Context Navigation

source: branches/publications/ORCHIDEE-ICE_SurfaceMassBalance/src_driver/forcing_tools.f90 @ 8398

Last change on this file since 8398 was 4287, checked in by josefine.ghattas, 7 years ago
Enhencement on clear functions. Added call to clear functions from all offline drivers. See ticket #232
File size: 164.6 KB

Line
1	! ==============================================================================================================================\n
2	! MODULE forcing_tools : This module concentrates on the temporal interpolation of the forcing for ORCHIDEE.
3	! It provides basic service for the grid when this is provided in the forcing file. The main
4	! work for the grid is done in glogrid.f90. The approach of forcing_tools to handle the time
5	! aspect of the forcing is to read as many time steps as possible in memory and then
6	! interpolate that to the time interval requested by the calling program.
7	! The data is read on root_proc but then distributed over all processors according to the
8	! domain decomposition of ORCHIDEE. This allows to be more efficient in the temporal interpolation.
9	! It is important to keep in mind that forcing_tools works on time intervals. So the request for data
10	! of ORCHIDEE as to be for an interval and the forcing file needs to have a description of the time interval
11	! over which the forcing is valid.
12	! The general description of how the attributes needed in the netCDF file for describing the cell_methods
13	! for time are provided in this document :
14	! https://forge.ipsl.jussieu.fr/orchidee/attachment/wiki/Documentation/Forcings/Description_Forcing_Files.pdf
15	!
16	! The most important routines of foring_tools are forcing_open and forcing_getvalues
17	!
18	! forcing_integration_time : Computes the interval over which the simulation should be carried out.
19	! forcing_open : Opens the forcing files and extract the main information.
20	! forcing_getvalues : Gets the forcing data for a time interval.
21	! forcing_close : Closes the forcing file
22	! forcing_printdate : A tool to print the dates in human readable form.
23	! forcing_printpoint : Print the values for a given point in time.
24	! forcing_givegridsize : Allows other programs to get the dimensions of the forcing grid.
25	! forcing_getglogrid : Allows other programs to get the spatial grid of the forcing.
26	! forcing_givegrid : Returns the description of the grid.
27	! forcing_zoomgrid : Extract a sub-region of the forcing grid.
28	!
29	! CONTACT : jan.polcher@lmd.jussieu.fr
30	!
31	! LICENCE : IPSL (2016)
32	! This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC
33	!
34	!>\BRIEF
35	!!
36	!! RECENT CHANGE(S): None
37	!!
38	!! REFERENCE(S) : None
39	!!
40	!_ ================================================================================================================================
41	!!
42	MODULE forcing_tools
43	!
44	USE defprec
45	USE netcdf
46	!
47	USE ioipsl
48	USE constantes
49	USE solar
50	!
51	USE mod_orchidee_para
52	!
53	IMPLICIT NONE
54	!
55	PRIVATE
56	PUBLIC :: forcing_open, forcing_close, forcing_printdate, forcing_getvalues, forcing_printpoint,&
57	& forcing_getglogrid, forcing_givegridsize, forcing_givegrid, forcing_zoomgrid, forcing_integration_time
58	PUBLIC :: forcing_tools_clear
59	!
60	!
61	!
62	INTERFACE forcing_reindex
63	MODULE PROCEDURE forcing_reindex3d, forcing_reindex2dt, forcing_reindex2d, forcing_reindex1d, &
64	& forcing_reindex2to1, forcing_reindex1to2
65	END INTERFACE forcing_reindex
66	!
67	INTERFACE forcing_printpoint
68	MODULE PROCEDURE forcing_printpoint_forgrid, forcing_printpoint_gen
69	END INTERFACE forcing_printpoint
70	!
71	! This PARAMETER essentially manages the memory usage of the module as it
72	! determines how much of the forcing will be uploaded from the netCDF file into
73	! memory.
74	!
75	INTEGER(i_std), PARAMETER :: slab_size_max=1500
76	!
77	! Time variables, all in Julian days
78	!
79	INTEGER(i_std), PARAMETER :: nbtmethods=4
80	INTEGER(i_std), SAVE :: nbtax
81	INTEGER(i_std), SAVE :: nb_forcing_steps
82	REAL(r_std), SAVE :: global_start_date, global_end_date, forcing_tstep_ave
83	REAL(r_std), SAVE :: dt_sechiba_keep
84	!
85	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: time
86	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:,:) :: time_bounds
87	CHARACTER(LEN=20), SAVE, ALLOCATABLE, DIMENSION(:) :: time_axename, time_cellmethod
88	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: preciptime
89	INTEGER(i_std), SAVE, ALLOCATABLE, DIMENSION(:) :: time_sourcefile
90	INTEGER(i_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: time_id
91	LOGICAL, SAVE :: end_of_file
92	!
93	! Forcing file information
94	!
95	INTEGER(i_std), SAVE :: nb_forcefile=0
96	CHARACTER(LEN=100), SAVE, ALLOCATABLE, DIMENSION(:) :: forfilename
97	INTEGER(i_std), SAVE, ALLOCATABLE, DIMENSION(:) :: force_id, id_unlim
98	INTEGER(i_std), SAVE, ALLOCATABLE, DIMENSION(:) :: nb_atts, ndims, nvars
99	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: convtosec
100	INTEGER(i_std), SAVE, ALLOCATABLE, DIMENSION(:) :: nbtime_perfile
101	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: date0_file
102	REAL(r_std), SAVE :: startdate, forcingstartdate
103	!
104	! Descrition of global grid
105	!
106	INTEGER(i_std), SAVE :: iim_glo, jjm_glo, nbland_glo
107	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: lon_glo, lat_glo
108	INTEGER(i_std), SAVE, ALLOCATABLE, DIMENSION(:,:):: mask_glo
109	INTEGER(i_std), SAVE, ALLOCATABLE, DIMENSION(:) :: lindex_glo
110	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: contfrac_glo
111	LOGICAL, SAVE :: compressed
112	!
113	! Descrition of zoomed grid
114	!
115	LOGICAL, SAVE :: zoom_forcing = .FALSE.
116	INTEGER(i_std), SAVE :: iim_loc, jjm_loc, nbland_loc
117	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: lon_loc, lat_loc
118	INTEGER(i_std), SAVE, ALLOCATABLE, DIMENSION(:) :: lindex_loc
119	INTEGER(i_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: mask_loc
120	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: area_loc
121	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: contfrac_loc
122	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:,:,:):: corners_loc
123	! Number of land points per proc
124	INTEGER(i_std), SAVE :: nbland_proc
125	INTEGER(i_std), SAVE, ALLOCATABLE, DIMENSION(:) :: glolindex_proc
126	!-
127	!- Heigh controls and data
128	!-
129	LOGICAL, SAVE :: zfixed, zsigma, zhybrid, zlevels, zheight
130	LOGICAL, SAVE :: zsamelev_uv
131	REAL, SAVE :: zlev_fixed, zlevuv_fixed
132	REAL, SAVE :: zhybrid_a, zhybrid_b
133	REAL, SAVE :: zhybriduv_a, zhybriduv_b
134	LOGICAL, SAVE :: lwdown_cons
135	!
136	! Forcing variables to be read and stored
137	!
138	! At 3000 we can fit in the slab an entire year of 3 hourly forcing.
139	INTEGER(i_std), SAVE :: slab_size=-1
140	INTEGER(i_std), SAVE :: current_offset=1
141	INTEGER(i_std), SAVE :: position_slab(2)
142	CHARACTER(LEN=20), SAVE :: calendar
143	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: tair_slab, qair_slab
144	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: time_tair, time_qair
145	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: timebnd_tair, timebnd_qair
146	!
147	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: rainf_slab, snowf_slab
148	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: time_precip
149	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: timebnd_precip
150	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: preciptime_slab !! Variable needed to keep track of how much rainfall was already distributed
151	!
152	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: swdown_slab, lwdown_slab
153	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: time_swdown, time_lwdown
154	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: timebnd_swdown, timebnd_lwdown
155	!
156	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: u_slab, v_slab, ps_slab
157	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: time_u, time_v, time_ps
158	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: timebnd_u, timebnd_v, timebnd_ps
159	!
160	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: ztq_slab, zuv_slab
161	INTEGER(i_std), SAVE, ALLOCATABLE, DIMENSION(:) :: reindex_glo, reindex_loc
162	INTEGER(i_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: reindex2d_loc
163	INTEGER(i_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: origind
164	!
165	INTEGER(i_std), SAVE :: ncdfstart, ncdfcount
166	!
167	! Flags to activate initialization phase
168	LOGICAL, SAVE :: first_call_readslab=.TRUE. !! Activate initialization phase in forcing_readslab_root
169	LOGICAL, SAVE :: first_call_solarint=.TRUE. !! Activate initialization phase in forcing_solarint
170	LOGICAL, SAVE :: first_call_spreadprec=.TRUE. !! Activate initialization phase in forcing_spreadprec
171
172
173	CONTAINS
174	!!
175	!! =============================================================================================================================
176	!! SUBROUTINE: forcing_integration_time
177	!!
178	!>\BRIEF Computes the interval over which the simulation should be carried out
179	!!
180	!! DESCRIPTION: This routing will get the following parameters from the run.def : 'START_DATE', 'END_DATE' and 'DT_SECHIBA'.
181	!! It allows to define the integration time of ORCHIDEE and later it will be used to verify that we have
182	!! the needed data in the forcing files to perform this simulation.
183	!!
184	!! \n
185	!_ ==============================================================================================================================
186	!!
187	SUBROUTINE forcing_integration_time(date_start, dt, nbdt)
188	!
189	!
190	! This subroutine gets the start date of the simulation, the time step and the number
191	! of time steps we need to do until the end of the simulations.
192	!
193	!
194	!
195	REAL(r_std), INTENT(out) :: date_start !! The date at which the simulation starts
196	REAL(r_std), INTENT(out) :: dt !! Time step length in seconds
197	INTEGER(i_std), INTENT(out) :: nbdt !! Number of timesteps to be executed
198	!
199	! Local
200	!
201	CHARACTER(LEN=20) :: str_sdate(2), str_edate(2), tmpstr
202	INTEGER(i_std) :: s_year, s_month, s_day, e_year, e_month, e_day
203	INTEGER(i_std) :: seci, hours, minutes
204	REAL(r_std) :: s_sec, e_sec, dateend, diff_sec, date_end
205	INTEGER(i_std) :: i, ic
206	!
207	!Config Key = START_DATE
208	!Config Desc = Date at which the simulation starts
209	!Config Def = NONE
210	!Config Help = The format is the same as in the CF convention : 1999-09-13 12:0:0
211	str_sdate = " "
212	CALL getin('START_DATE',str_sdate)
213	!
214	IF ( (INDEX(str_sdate(1),"-") .NE. INDEX(str_sdate(1),"-", .TRUE.)) .AND. &
215	& (INDEX(str_sdate(2),":") .NE. INDEX(str_sdate(2),":", .TRUE.)) ) THEN
216	DO i=1,2
217	tmpstr = str_sdate(1)
218	ic = INDEX(tmpstr,"-")
219	tmpstr(ic:ic) = " "
220	str_sdate(1) = tmpstr
221	tmpstr = str_sdate(2)
222	ic = INDEX(tmpstr,":")
223	tmpstr(ic:ic) = " "
224	str_sdate(2) = tmpstr
225	ENDDO
226	READ (str_sdate(1),*) s_year, s_month, s_day
227	READ (str_sdate(2),*) hours, minutes, seci
228	s_sec = hours3600. + minutes60. + seci
229	ELSE
230	CALL ipslerr(3, "forcing_integration_time", "START_DATE incorrectly specified in run.def", str_sdate(1), str_sdate(2))
231	ENDIF
232	CALL ymds2ju (s_year, s_month, s_day, s_sec, date_start)
233	CALL forcing_printdate(date_start, "This is after reading the start date")
234	!
235	!Config Key = END_DATE
236	!Config Desc = Date at which the simulation ends
237	!Config Def = NONE
238	!Config Help = The format is the same as in the CF convention : 1999-09-13 12:0:0
239	str_edate = " "
240	CALL getin('END_DATE',str_edate)
241	!
242	IF ( (INDEX(str_edate(1),"-") .NE. INDEX(str_edate(1),"-", .TRUE.)) .AND. &
243	& (INDEX(str_edate(2),":") .NE. INDEX(str_edate(2),":", .TRUE.)) ) THEN
244	DO i=1,2
245	tmpstr = str_edate(1)
246	ic = INDEX(tmpstr,"-")
247	tmpstr(ic:ic) = " "
248	str_edate(1) = tmpstr
249	tmpstr = str_edate(2)
250	ic = INDEX(tmpstr,":")
251	tmpstr(ic:ic) = " "
252	str_edate(2) = tmpstr
253	ENDDO
254	READ (str_edate(1),*) e_year, e_month, e_day
255	READ (str_edate(2),*) hours, minutes, seci
256	e_sec = hours3600. + minutes60. + seci
257	ELSE
258	CALL ipslerr(3, "forcing_integration_time", "END_DATE incorrectly specified in run.def", str_edate(1), str_edate(2))
259	ENDIF
260	CALL ymds2ju (e_year, e_month, e_day, e_sec, date_end)
261	!
262	CALL time_diff (s_year,s_month,s_day,s_sec,e_year,e_month,e_day,e_sec,diff_sec)
263	!
264	!Config Key = DT_SECHIBA
265	!Config Desc = Time step length in seconds for sechiba component
266	!Config Def = 1800
267	!Config Help =
268	!Config Units = [seconds]
269	dt = 1800
270	CALL getin('DT_SECHIBA', dt)
271	dt_sechiba_keep = dt
272	!
273	nbdt = NINT(diff_sec/dt)
274	!
275	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
276	!
277	! Read the configuration options for the time interpolations.
278	!
279	!Config Key = LWDOWN_CONS
280	!Config Desc = Conserve the longwave downward radiation of the forcing
281	!Config Def = n
282	!Config Help = This flag allows to conserve the downward longwave radiation
283	! provided in the forcing. It will do this by taking the closest
284	! neighbour in time from the forcing. This assumes that the forcing
285	! contains average fluxes. The default setting (LWDOWN_CONS=n) will
286	! force the model to perform a linear interpolation of the fluxes.
287	!Config Units = [FLAG]
288	!-
289	lwdown_cons = .FALSE.
290	CALL getin('LWDOWN_CONS', lwdown_cons)
291	!
292	END SUBROUTINE forcing_integration_time
293	!!
294	!! =============================================================================================================================
295	!! SUBROUTINE: forcing_open
296	!!
297	!>\BRIEF Opens the forcing files and extract the main information.
298	!!
299	!! DESCRIPTION: This routine opens all the forcing files provided in the list and verifies that the grid corresponds
300	!! to the coordinates provided (and which was obtained by the model from glogrid.f90.). It then zooms
301	!! into the forcing as requested by the user, extracts the vertical coordinates and final reads the time axis.
302	!! Some basic consistency checks are performed as for instance ensuring the that all the forcing data is available
303	!! to simulate the desired period.
304	!! All that information is also broadcasted to all processors.
305	!! Actual forcing data is not read at this stage.
306	!!
307	!! \n
308	!_ ==============================================================================================================================
309	!
310	SUBROUTINE forcing_open(filenames_in, iim, jjm, lon, lat, nbland_in, drvzoom_lon, drvzoom_lat, &
311	& kindex, nbindex_perproc, wunit)
312	!
313	! Opens the forcing file and reads some key information and stores them in the shared variables of the
314	! module.
315	!
316	! Lon, lat should come from the grid file we read before. This will give indication of the grid
317	! file is consistant with the forcing file and if we need to zoom into the forcing file.
318	!
319	! Time interval of the simulation is also determined.
320	!
321	! ARGUMENTS
322	!
323	CHARACTER(LEN=*), INTENT(in) :: filenames_in(:)
324	INTEGER(i_std), INTENT(in) :: iim, jjm, nbland_in
325	REAL(r_std), INTENT(in) :: lon(iim,jjm), lat(iim,jjm)
326	REAL(r_std), DIMENSION(2), INTENT(in) :: drvzoom_lon, drvzoom_lat
327	INTEGER(i_std), INTENT(in) :: kindex(nbland_in)
328	INTEGER(i_std), INTENT(in) :: nbindex_perproc
329	INTEGER(i_std), OPTIONAL :: wunit
330	!
331	! LOCAL
332	!
333	INTEGER(i_std) :: iim_tmp, jjm_tmp, nbland_tmp, nb_files
334	INTEGER(i_std) :: iv, it
335	INTEGER(i_std) :: inl, ii, jj, ik
336	INTEGER(i_std) :: land_id
337	REAL(r_std) :: dt
338	INTEGER(i_std) :: nbdt
339	!
340	! How many files do we have to open ?
341	!
342	! Number of points per processor
343	nbland_proc = nbindex_perproc
344	!
345	! All the meta information from the forcing file is ojnly needed on the root processor.
346	!
347	IF ( is_root_prc ) THEN
348	!
349	CALL forcing_filenamecheck(filenames_in, nb_files)
350	IF ( PRESENT(wunit) ) THEN
351	DO it=1,nb_files
352	WRITE(wunit,*) "Files to be used for forcing the simulation :", it, TRIM(forfilename(it))
353	ENDDO
354	ENDIF
355	!
356	! 0.0 Check if variables are allocated to the right size on root_proc
357	!
358	IF (nb_files > nb_forcefile) THEN
359	IF ( ALLOCATED(force_id) ) DEALLOCATE(force_id)
360	ALLOCATE(force_id(nb_files))
361	IF ( ALLOCATED(id_unlim) ) DEALLOCATE(id_unlim)
362	ALLOCATE(id_unlim(nb_files))
363	IF ( ALLOCATED(nb_atts) ) DEALLOCATE(nb_atts)
364	ALLOCATE(nb_atts(nb_files))
365	IF ( ALLOCATED(ndims) ) DEALLOCATE(ndims)
366	ALLOCATE(ndims(nb_files))
367	IF ( ALLOCATED(nvars) ) DEALLOCATE(nvars)
368	ALLOCATE( nvars(nb_files))
369	IF ( ALLOCATED(nbtime_perfile) ) DEALLOCATE(nbtime_perfile)
370	ALLOCATE(nbtime_perfile(nb_files))
371	IF ( ALLOCATED(convtosec) ) DEALLOCATE(convtosec)
372	ALLOCATE(convtosec(nb_files))
373	ENDIF
374	nb_forcefile = nb_files
375	!
376	! Get the global grid size from the forcing file. The output is in temporary variables as in this
377	! module the values are shared.
378	!
379	IF ( PRESENT(wunit) ) THEN
380	WRITE(wunit,*) "Getting global grid from ", nb_forcefile, "files."
381	CALL FLUSH(wunit)
382	ENDIF
383	CALL forcing_getglogrid(nb_forcefile, forfilename, iim_tmp, jjm_tmp, nbland_tmp, .FALSE.)
384	!
385	IF ( PRESENT(wunit) ) THEN
386	WRITE(wunit,*) "Getting the zoomed grid", nbland_tmp
387	CALL FLUSH(wunit)
388	ENDIF
389	CALL forcing_zoomgrid(drvzoom_lon, drvzoom_lat, forfilename(1), .FALSE.)
390	IF ( PRESENT(wunit) ) THEN
391	WRITE(wunit,*) "Out of the zoomed grid operation"
392	CALL FLUSH(wunit)
393	ENDIF
394	!
395	! Verification that the grid sizes coming from the calling program are consistant with what we get
396	! from the forcing file.
397	!
398	IF ( (iim_loc .NE. iim) .OR. (jjm_loc .NE. jjm) ) THEN
399	CALL ipslerr (3,'forcing_open',"At least one of the dimensions of the grid obtained from the",&
400	& "grid file is different from the one in the forcing file.",&
401	& "Run driver2oasis -init to generate a new grid file.")
402	ENDIF
403	! Special treatment for the number of land point, as we could have a case where the forcing
404	! file does not include the land/sea mask.
405	!
406	IF ( nbland_loc .NE. nbland_in ) THEN
407	! We trust the number of land points obtained from the gridfile. It has the land/sea mask.
408	nbland_loc = nbland_in
409	ENDIF
410	!
411	! Treat the time dimension now :
412	!
413	IF ( PRESENT(wunit) ) THEN
414	WRITE(wunit,*) "Getting forcing time"
415	CALL FLUSH(wunit)
416	ENDIF
417	CALL forcing_time(nb_forcefile, forfilename)
418	!
419	! Now that we know how much time steps are in the forcing we can set some realistic slab_size
420	!
421	slab_size=MIN(nb_forcing_steps, slab_size_max)
422	!
423	!
424	! Get the vertical information from the file
425	!
426	CALL forcing_vertical(force_id(1))
427	!
428	!
429	IF ( PRESENT(wunit) ) THEN
430	WRITE(wunit,*) "Getting integration time"
431	CALL FLUSH(wunit)
432	ENDIF
433	CALL forcing_integration_time(startdate, dt, nbdt)
434	!
435	! Test that the time interval requested by the user correspond to the time available in the
436	! forcing file.
437	!
438	IF ( startdate < time_bounds(1,1,1) .OR. startdate > time_bounds(nb_forcing_steps,1,2) ) THEN
439	CALL ipslerr (3,'forcing_open', 'Start time requested by the user is outside of the time interval',&
440	& "covered by the forcing file.","Please verify the configuration in the run.def file.")
441	ENDIF
442	!
443	IF ( startdate+(dt/one_day)*nbdt > time_bounds(nb_forcing_steps,1,2) .OR. &
444	& startdate+(dt/one_day)*nbdt < time_bounds(1,1,1)) THEN
445	CALL forcing_printdate(time_bounds(nb_forcing_steps,1,2), "Outer bound of forcing file.")
446	CALL forcing_printdate(startdate+(dt/one_day)*nbdt, "Last date to be simulated.")
447	WRITE(,) "ERROR : Final date of forcing needed is : ", startdate+(dt/one_day)*nbdt
448	WRITE(,) "ERROR : The outer bound of the last forcing time step is :", time_bounds(nb_forcing_steps,1,2)
449	CALL ipslerr (3,'forcing_open', 'End time requested by the user is outside of the time interval',&
450	& "covered by the forcing file.","Please verify the configuration in the run.def file.")
451	ENDIF
452	!
453	ENDIF
454	!
455	! Broadcast the local grid (i.e. the one resulting from the zoom) to all processors
456	!
457	CALL bcast(iim_loc)
458	CALL bcast(jjm_loc)
459	CALL bcast(nbland_loc)
460	! Time variables needed by all procs
461	CALL bcast(slab_size)
462	CALL bcast(startdate)
463	CALL bcast(forcingstartdate)
464	CALL bcast(forcing_tstep_ave)
465	!
466	! On the slave processes we need to allocate the memory for the data on root_prc to be bcast
467	! On the root_proc these allocations were done with CALL forcing_zoomgrid
468	!
469	ALLOCATE(glolindex_proc(nbland_proc))
470	IF ( .NOT. is_root_prc ) THEN
471	ALLOCATE(lon_loc(iim_loc,jjm_loc))
472	ALLOCATE(lat_loc(iim_loc,jjm_loc))
473	ALLOCATE(lindex_loc(nbland_loc))
474	ALLOCATE(mask_loc(iim_loc,jjm_loc))
475	ALLOCATE(area_loc(iim_loc,jjm_loc))
476	ALLOCATE(contfrac_loc(nbland_loc))
477	ALLOCATE(corners_loc(iim_loc,jjm_loc,4,2))
478	ENDIF
479	!
480	! Keep on each processor the index of each land point on the *_loc grid
481	!
482	CALL scatter(kindex, glolindex_proc)
483	!
484	CALL bcast(lon_loc)
485	CALL bcast(lat_loc)
486	CALL bcast(lindex_loc)
487	CALL bcast(mask_loc)
488	CALL bcast(area_loc)
489	CALL bcast(contfrac_loc)
490	CALL bcast(corners_loc)
491	!
492	END SUBROUTINE forcing_open
493	!!
494	!! =============================================================================================================================
495	!! SUBROUTINE: forcing_getvalues
496	!!
497	!>\BRIEF Gets the forcing data for a time interval.
498	!!
499	!! DESCRIPTION: The routine will get the forcing valid for the time interval provided by the caller.
500	!! First it will check that the data is already in memory for that time interval. If not
501	!! it will first read the data from the netCDF file.
502	!! Then the forcing date will be interpolated to the requested time interval.
503	!! The code calls linear interpolation for most variables except for SWdown and precipitation.
504	!! These temporal interpolations can be improved later.
505	!!
506	!! \n
507	!_ ==============================================================================================================================
508	SUBROUTINE forcing_getvalues(time_int, dt, zlev_tq, zlev_uv, tair, qair, rainf, snowf, &
509	& swdown, lwdown, solarang, u, v, ps)
510	!
511	! ARGUMENTS
512	!
513	REAL(r_std), INTENT(in) :: time_int(2) !! The time interval over which the forcing is needed.
514	REAL(r_std), INTENT(in) :: dt !! timestep, i.e. distance in seconds between time_int(1) and time_int(2)
515	REAL(r_std), INTENT(out) :: zlev_tq(:), zlev_uv(:)
516	REAL(r_std), INTENT(out) :: tair(:), qair(:), rainf(:), snowf(:)
517	REAL(r_std), INTENT(out) :: swdown(:), lwdown(:), solarang(:)
518	REAL(r_std), INTENT(out) :: u(:), v(:), ps(:)
519	!
520	! LOCAL
521	!
522	INTEGER(i_std) :: i
523	!
524	! Test that we have the time interval within our slab of data else we need to update it.
525	! Att : the tests are done here on time_tair as an exemple. This might need to have to be generalized.
526	!
527	! First case the time axis of the variable are not even yet allocated !
528	IF ( .NOT. ALLOCATED(time_tair) ) THEN
529	CALL forcing_readslab(time_int)
530	CALL forcing_printdate(timebnd_tair(1,1), "Start of time slab just read")
531	CALL forcing_printdate(timebnd_tair(slab_size,2), "End of time slab just read")
532	ELSE
533	! If we have time axis (for TAIR here) we test that it is long enough in time to allow for an interpolation.
534	!
535	IF ( time_int(2)+forcing_tstep_ave/one_day > time_tair(slab_size) .AND. (.NOT. end_of_file) ) THEN
536	CALL forcing_readslab(time_int)
537	CALL forcing_printdate(timebnd_tair(1,1), "Start of time slab just read")
538	CALL forcing_printdate(timebnd_tair(slab_size,2), "End of time slab just read")
539	ENDIF
540	ENDIF
541	!
542	! Interpolate the dynamical variables to the time step at which the driver is for the moment.
543	!
544	CALL forcing_interpol(time_int, dt, time_u, u_slab, u)
545	CALL forcing_interpol(time_int, dt, time_v, v_slab, v)
546	CALL forcing_interpol(time_int, dt, time_ps, ps_slab, ps)
547	!
548	! Compute the height of the first level (a routine will be needed for that !)
549	! ATT : we assume that the time axis for the height of the scalar variable is the one of TAIR
550	! and for the height of wind is the same as U.
551	CALL forcing_interpol(time_int, dt, time_tair, ztq_slab, zlev_tq)
552	CALL forcing_interpol(time_int, dt, time_u, zuv_slab, zlev_uv)
553	!
554	! Interpolate the state variables of the lower atmospheric level
555	!
556	CALL forcing_interpol(time_int, dt, time_tair, tair_slab, tair)
557	CALL forcing_interpol(time_int, dt, time_qair, qair_slab, qair)
558	!
559	! Spread the precipitation as requested by the user
560	!
561	CALL forcing_spreadprec(time_int, dt, timebnd_precip, time_precip, rainf, snowf)
562	!
563	! Deal with the interpolate of the radiative fluxes.
564	!
565	CALL forcing_solarint(time_int, dt, timebnd_swdown, time_swdown, iim_loc, jjm_loc, lon_loc, lat_loc, swdown, solarang)
566	!
567	! We have the option here to conserve LWdown by taking the closest point in the forcing.
568	! So no interpolation is done.
569	!
570	IF ( lwdown_cons ) THEN
571	CALL forcing_closest(time_int, dt, time_lwdown, lwdown_slab, lwdown)
572	ELSE
573	CALL forcing_interpol(time_int, dt, time_lwdown, lwdown_slab, lwdown)
574	ENDIF
575
576	END SUBROUTINE forcing_getvalues
577
578	!! =============================================================================================================================
579	!! SUBROUTINE: forcing_closest
580	!!
581	!>\BRIEF This routine does not interpolate and simply uses the closes value in time. It is useful for preserving
582	!! variables which are averaged in the forcing file.
583	!!
584	!! DESCRIPTION:
585	!!
586	!! \n
587	!_ ==============================================================================================================================
588	SUBROUTINE forcing_closest(time_int_in, dt, time_central_in, var_slab, var)
589	!
590	! ARGUMENTS
591	!
592	REAL(r_std), INTENT(in) :: time_int_in(2)
593	REAL(r_std), INTENT(in) :: dt
594	REAL(r_std), INTENT(in) :: time_central_in(:)
595	REAL(r_std), INTENT(in) :: var_slab(:,:)
596	REAL(r_std), INTENT(out) :: var(:)
597	!
598	! LOCAL
599	!
600	INTEGER(i_std) :: slabind_a, slabind_b, imin(1), i
601	REAL(r_std) :: time_int(2), time_central(slab_size_max)
602	REAL(r_std) :: mid_int, wa, wb, wt, wab, wae, tmp_mid_int
603	LOGICAL :: mask(slab_size_max)=.FALSE.
604	!
605	! Shift the input dates in order to gain in precision for the calculations
606	!
607	time_int(:) = time_int_in(:)-INT(forcingstartdate)
608	time_central(1:slab_size) = time_central_in(1:slab_size)-INT(forcingstartdate)
609	!
610	! Create a mask so that MINLOC does not look outside of the valid interval of time_central
611	!
612	mask(1:slab_size) = .TRUE.
613	!
614	! Select the forcing interval for which the center date is the closest to the time of
615	! the model.
616	!
617	mid_int = time_int(1) + (dt/2.0)/one_day
618	imin = MINLOC( ABS(time_central(1:slab_size) - mid_int), mask )
619	!
620	! Verify that this is a possible date
621	!
622	IF ( imin(1) > 0 .AND. imin(1) <= slab_size ) THEN
623	!
624	slabind_a = imin(1)
625	!
626	ELSE
627	WRITE(,) "imin(1) = ", imin(1), (time_int_in(1) + (dt/2.0)/one_day)
628	CALL forcing_printdate(time_int_in(1), "===> Start of target time interval.")
629	CALL forcing_printdate(time_int_in(2), "===> End of target time interval.")
630	CALL forcing_printdate(time_central_in(imin(1)), "===> Center of forcing time interval.")
631	CALL ipslerr (3,'forcing_closest', 'The target time interval has no acceptable closest',&
632	& "time in the forcing slab.","")
633	ENDIF
634	!
635	! Transfer the data from the sloest time of the forcing data slab.
636	!
637	DO i=1, nbland_proc
638	!
639	var(i) = var_slab(i,slabind_a)
640	!
641	ENDDO
642	!
643	!
644	END SUBROUTINE forcing_closest
645
646	!! =============================================================================================================================
647	!! SUBROUTINE: forcing_interpol
648	!!
649	!>\BRIEF Perform linear interpolation for the time interval requested.
650	!!
651	!! DESCRIPTION:
652	!! The code gets an interval over which the model will integrate (time_int_in) but only uses the centre. It also gets
653	!! the times representative of the forcing data for the variable at hand (time_central_in). Using this data we will
654	!! determine which 2 forcing times will need to be used for the interpolation. Once this is established the weights
655	!! are computed and used in order to interpolate the variable between the 2 times which bracket the model integration time.
656	!! \n
657	!_ ==============================================================================================================================
658	SUBROUTINE forcing_interpol(time_int_in, dt, time_central_in, var_slab, var)
659	!
660	! ARGUMENTS
661	!
662	REAL(r_std), INTENT(in) :: time_int_in(2) !! The time interval over which the forcing is needed by the model.
663	REAL(r_std), INTENT(in) :: dt !! Time step of the model
664	REAL(r_std), INTENT(in) :: time_central_in(:) !! Representative time for the interval of validity of the forcing data
665	REAL(r_std), INTENT(in) :: var_slab(:,:) !! The slab of forcing data read from the file.
666	REAL(r_std), INTENT(out) :: var(:) !! Result of the time interpolation.
667	!
668	! LOCAL
669	!
670	INTEGER(i_std) :: slabind_a, slabind_b, imin(1), i
671	REAL(r_std) :: time_int(2), time_central(slab_size_max)
672	REAL(r_std) :: mid_int, wa, wb, wt, wab, wae, tmp_mid_int
673	LOGICAL :: mask(slab_size_max)=.FALSE.
674	!
675	! Shift the input dates in order to gain in precision for the calculations
676	!
677	time_int(:) = time_int_in(:)-INT(forcingstartdate)
678	time_central(1:slab_size) = time_central_in(1:slab_size)-INT(forcingstartdate)
679	!
680	! Create a mask so that MINLOC does not look outside of the valid interval of time_central
681	!
682	mask(1:slab_size) = .TRUE.
683	!
684	! Select the type of interpolation to be done.
685	!
686	! Compute the central time of the model integration time.
687	!
688	mid_int = time_int(1) + (dt/2.0)/one_day
689	! Locate that time on the time axis of the forcing.
690	imin = MINLOC( ABS(time_central(1:slab_size) - mid_int), mask )
691	!
692	! Determine which indices are to the left (slabind_a) and right (slabind_b) of the model time and will be used
693	! for the linear interpolation.
694	!
695	IF ( imin(1) > 1 .AND. imin(1) < slab_size ) THEN
696	!
697	! Determine if the model time is to the left or right of the representative time
698	! of the forcing data. This allows to determine with which other position in the
699	! forcing data we need to interpolate.
700	!
701	IF ( mid_int < time_central(imin(1)) ) THEN
702	slabind_a = imin(1) - 1
703	slabind_b = imin(1)
704	ELSE
705	slabind_a = imin(1)
706	slabind_b = imin(1) + 1
707	ENDIF
708	!
709	ELSE IF ( imin(1) == 1 ) THEN
710	!
711	! If we are at the first time step of the forcing data we need to take care as there is
712	! no data earlier.
713	!
714	slabind_a = 1
715	slabind_b = 2
716	IF ( mid_int < time_central(slabind_a) ) THEN
717	IF ( time_int(2) < time_central(slabind_a) ) THEN
718	WRITE(,) "imin(1) = ", imin(1), (time_int_in(1) + (dt/2.0)/one_day)
719	CALL forcing_printdate(time_int_in(1), "===> Start of target time interval.")
720	CALL forcing_printdate(time_int_in(2), "===> End of target time interval.")
721	CALL forcing_printdate(time_central_in(slabind_a), "===> Center of forcing time interval.")
722	CALL ipslerr (3,'forcing_interpol', 'The target time interval lies before the first date of the slab.',&
723	& "","")
724	ELSE
725	mid_int = time_central(slabind_a)
726	ENDIF
727	ENDIF
728	ELSE IF ( imin(1) == slab_size ) THEN
729	!
730	! If we are at the end of the forcing data we need to pay attention as we have no data later in time.
731	!
732	slabind_a = slab_size - 1
733	slabind_b = slab_size
734	IF ( mid_int > time_central(slabind_b) ) THEN
735	IF ( time_int(1) > time_central(slabind_b) ) THEN
736	WRITE(,) "imin(1) = ", imin(1), (time_int_in(1) + (dt/2.0)/one_day)
737	CALL forcing_printdate(time_int_in(1), "===> Start of target time interval.")
738	CALL forcing_printdate(time_int_in(2), "===> End of target time interval.")
739	CALL forcing_printdate(time_central_in(slabind_b), "===> Center of forcing time interval.")
740	CALL ipslerr (3,'forcing_interpol', 'The target time interval lies after the last date of the slab.',&
741	& "","")
742	ELSE
743	mid_int = time_central(slabind_b)
744	ENDIF
745	ENDIF
746	ENDIF
747	!
748	! Compute the weights between the values at slabind_a and slabind_b. As with the time
749	! representation we are at the limit of precision we use 2 days to compute the distance
750	! in time between the first value (slabind_a) and the middle of the target interval.
751	!
752	wab = time_int(1) - time_central(slabind_a) + (dt/2.0)/one_day
753	wae = time_int(2) - time_central(slabind_a) - (dt/2.0)/one_day
754	wa = (wab+wae)/2.0
755	wb = time_central(slabind_b) - time_central(slabind_a)
756	wt = wa/wb
757	!
758	! Do the weighted average of all land points with the time indices and weights computed above.
759	!
760	DO i=1, nbland_proc
761	var(i) = var_slab(i,slabind_a) + wt*(var_slab(i,slabind_b) - var_slab(i,slabind_a))
762	ENDDO
763
764	END SUBROUTINE forcing_interpol
765
766
767	!! =============================================================================================================================
768	!! SUBROUTINE: forcing_spreadprec
769	!!
770	!>\BRIEF Spreads the precipitation over the interval chosen based on the interval chosen by the user.
771	!!
772	!! DESCRIPTION: The behaviour of this routine is controlled by the parameter SPRED_PREC_SEC in the run.def.
773	!! The time in second specified by the user will be the one over which the precipitation will last
774	!! where the forcing interval has rain or snow.
775	!!
776	!! \n
777	!_ ==============================================================================================================================
778	SUBROUTINE forcing_spreadprec(time_int, tlen, timebnd_central, time_central, rainf, snowf)
779	!
780	! ARGUMENTS
781	!
782	REAL(r_std), INTENT(in) :: time_int(2) ! Time interval to which we will spread precip
783	REAL(r_std), INTENT(in) :: tlen ! size of time interval in seconds (time step !)
784	REAL(r_std), INTENT(in) :: timebnd_central(:,:) ! Time interval over which the read data is valid
785	REAL(r_std), INTENT(in) :: time_central(:) ! Center of the time interval
786	REAL(r_std), INTENT(out) :: rainf(:), snowf(:)
787	!
788	! LOCAL
789	!
790	REAL(r_std), SAVE :: time_to_spread=3600.0
791	INTEGER(i_std) :: imin(1), i, tind(3)
792	REAL(r_std) :: ft(3), dt, left, right
793	INTEGER(i_std) :: offset, nb_spread
794	LOGICAL :: mask(slab_size_max)=.FALSE.
795	!
796	IF ( first_call_spreadprec ) THEN
797	!Config Key = SPRED_PREC
798	!Config Desc = Spread the precipitation.
799	!Config If = [-]
800	!Config Def = Half of the forcing time step or uniform, depending on dt_force and dt_sechiba
801	!Config Help = Spread the precipitation over SPRED_PREC steps of the splited forcing
802	!Config time step. This ONLY applied if the forcing time step has been splited.
803	!Config If the value indicated is greater than SPLIT_DT, SPLIT_DT is used for it.
804	!Config Units = [-]
805	!-
806	nb_spread = -1
807	CALL getin_p('SPRED_PREC', nb_spread)
808	!
809	! Test if we have read the number of time steps to spread in run.def
810	! If not, then probably the time was given in seconds.
811	!
812	IF ( nb_spread < 0 ) THEN
813	!Config Key = SPRED_PREC_SEC
814	!Config Desc = Spread the precipitation over an interval in seconds.
815	!Config Def = 3600
816	!Config Help = Spread the precipitation over n seconds of the forcing time step
817	!Config interval. This ONLY applies when the SPRED_PREC_SEC is smaller than
818	!Config the forcing time step. Should the user set SPRED_PREC_SEC=0 we will
819	!Config assume that the rainfall is uniformely distributed over the forcing interval.
820	!Config Units = seconds
821	!
822	! This is the default should 'SPRED_PREC' not be present in the run.def
823	!
824	time_to_spread = forcing_tstep_ave/2.0
825	!
826	CALL getin_p('SPRED_PREC_SEC', time_to_spread)
827	ELSE
828	time_to_spread = dt_sechiba_keep * nb_spread
829	ENDIF
830	!
831	! Do some verifications on the information read from run.def
832	!
833	IF ( time_to_spread > forcing_tstep_ave) THEN
834	time_to_spread = forcing_tstep_ave
835	ELSE IF ( time_to_spread <= 0 ) THEN
836	time_to_spread = forcing_tstep_ave
837	ENDIF
838	!
839	first_call_spreadprec = .FALSE.
840	!
841	ENDIF
842	!
843	! First test that we have the right time interval from the forcing to spread the precipitation
844	!
845	IF ( time_int(1) >= timebnd_central(1,1) .AND. time_int(2) <= timebnd_central(slab_size,2)) THEN
846	!
847	! Create a mask so that MINLOC does not look outside of the valid interval of time_central
848	!
849	mask(1:slab_size) = .TRUE.
850	!
851	! To get better precision on the time difference we get a common offset to substract
852	!
853	offset = INT(forcingstartdate)
854	!
855	! In principle 3 time steps can contribute to the time step closest to the center of the forcing interval
856	!
857	imin = MINLOC( ABS(time_central(1:slab_size)-(time_int(1)+time_int(2))/2.0), mask )
858	tind(1) = MAX(imin(1)-1,1)
859	tind(2) = imin(1)
860	tind(3) = MIN(imin(1)+1,slab_size)
861	IF (imin(1)+1 > slab_size) THEN
862	WRITE(,) "We have a problem here imin(1)+1,slab_size ", imin(1)+1,slab_size
863	WRITE(,) "Interval : ", time_int(1),time_int(2)
864	ENDIF
865	!
866	!
867	!
868	! Do we need to take some rain from the previous time step ?
869	!
870	!! Time computation is not better than 1/1000 seconds
871	IF ( time_int(1) < timebnd_central(tind(2),1) .AND. preciptime_slab(tind(1)) < (time_to_spread-0.001) ) THEN
872	dt = ((timebnd_central(tind(2),1)-offset)-(time_int(1)-offset))*one_day
873	ft(1) = MIN(time_to_spread - preciptime_slab(tind(1)), dt)/tlen
874	ELSE
875	ft(1) = zero
876	ENDIF
877	!
878	! Is there still some rain to spread from the current forcing time step ?
879	!
880	!! Time computation is not better than 1/1000 seconds
881	IF (preciptime_slab(tind(2)) < (time_to_spread-0.001) ) THEN
882	left = MAX(time_int(1), timebnd_central(tind(2),1))
883	right = MIN(time_int(2),timebnd_central(tind(2),2))
884	dt = ((right-offset)-(left-offset))*one_day
885	ft(2) = MIN(time_to_spread - preciptime_slab(tind(2)), dt)/tlen
886	ELSE
887	ft(2) = zero
888	ENDIF
889	!
890	! Do we need to take some rain from the next time step ?
891	!
892	!! Time computation is not better than 1/1000 seconds
893	IF ( time_int(2) > timebnd_central(tind(2),2) .AND. preciptime_slab(tind(3)) < (time_to_spread-0.001) ) THEN
894	dt = ((time_int(2)-offset)-(timebnd_central(tind(2),2)-offset))*one_day
895	ft(3) = MIN(time_to_spread - preciptime_slab(tind(3)), dt)/tlen
896	ELSE
897	ft(3) = zero
898	ENDIF
899	!
900	! Do the actual calculation
901	!
902	DO i=1, nbland_proc
903	rainf(i) = (rainf_slab(i,tind(1)) * forcing_tstep_ave * ft(1) + &
904	& rainf_slab(i,tind(2)) * forcing_tstep_ave * ft(2) + &
905	& rainf_slab(i,tind(3)) * forcing_tstep_ave * ft(3))*tlen/time_to_spread
906
907	snowf(i) = (snowf_slab(i,tind(1)) * forcing_tstep_ave * ft(1) + &
908	& snowf_slab(i,tind(2)) * forcing_tstep_ave * ft(2) + &
909	& snowf_slab(i,tind(3)) * forcing_tstep_ave * ft(3))*tlen/time_to_spread
910	ENDDO
911	!
912	! Update the time over which we have already spread the rainf
913	!
914	preciptime_slab(tind(1)) = preciptime_slab(tind(1)) + tlen * ft(1)
915	preciptime_slab(tind(2)) = preciptime_slab(tind(2)) + tlen * ft(2)
916	preciptime_slab(tind(3)) = preciptime_slab(tind(3)) + tlen * ft(3)
917	!
918	ELSE
919	WRITE(numout,*) "Time interval toward which we will interpolate : ", time_int
920	WRITE(numout,*) "Limits of the time slab we have : ", timebnd_central(1,1), timebnd_central(slab_size,2)
921	CALL forcing_printdate(time_int(1), "Start of target time interval.")
922	CALL forcing_printdate(time_int(2), "End of target time interval.")
923	CALL forcing_printdate(timebnd_central(1,1), "Start of time slab we have.")
924	CALL forcing_printdate(timebnd_central(slab_size,2), "End of time slab we have.")
925	CALL ipslerr (3,'forcing_spreadprec', 'The sitation should not occur Why are we here ?',&
926	& "","")
927	ENDIF
928
929	END SUBROUTINE forcing_spreadprec
930
931	!! =============================================================================================================================
932	!! SUBROUTINE: forcing_solarint
933	!!
934	!>\BRIEF Interpolates incoming solar radiation to the interval requested.
935	!!
936	!! DESCRIPTION: The interpolation here takes into account the variation of the solar zenith angle
937	!! to ensure the diurnal cycle of solar radiation is as well represented as possible.
938	!!
939	!! \n
940	!_ ==============================================================================================================================
941	SUBROUTINE forcing_solarint(time_int_in, tlen, timebnd_in, time_cent_in, iim, jjm, lon, lat, swdown, solarangle)
942	!
943	! ARGUMENTS
944	!
945	REAL(r_std), INTENT(in) :: time_int_in(2) ! Time interval for which we will compute radiation
946	REAL(r_std), INTENT(in) :: tlen ! size of time interval in seconds (time step !)
947	REAL(r_std), INTENT(in) :: timebnd_in(:,:) ! Time interval over which the read data is valid
948	REAL(r_std), INTENT(in) :: time_cent_in(:) ! Center of the time interval
949	INTEGER(i_std), INTENT(in) :: iim, jjm ! Size of 2D domain
950	REAL(r_std), INTENT(in) :: lon(iim,jjm), lat(iim,jjm) ! Longitude and latitude
951	REAL(r_std), INTENT(out) :: swdown(:), solarangle(:) ! interpolated downward solar radiation and corresponding
952	! ! solar angle.
953	!
954	! LOCAL SAVED
955	!
956	REAL(r_std), SAVE :: solaryearstart
957	INTEGER(i_std), SAVE :: split, split_max
958	REAL(r_std), SAVE :: last_time
959	!
960	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: mean_sinang
961	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:,:) :: sinangles
962	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: time_angles
963	! Dusk-dawn management
964	REAL(r_std), SAVE :: dusk_angle
965	!
966	! LOCAL - temporary
967	!
968	REAL(r_std) :: time_int(2)
969	REAL(r_std) :: timebnd(slab_size_max,2)
970	REAL(r_std) :: time_cent(slab_size_max)
971	INTEGER(i_std) :: year, month, day, hours, minutes
972	REAL(r_std) :: sec
973	!
974	REAL(r_std) :: mean_sol, split_time
975	REAL(r_std) :: julian, julian_tmp
976	REAL(r_std) :: sinang(iim,jjm)
977	INTEGER(i_std) :: is, i, ii, jj, imin(1), tmin(1)
978	LOGICAL :: mask(slab_size_max)=.FALSE.
979	!
980	IF ( first_call_solarint ) THEN
981	!
982	! Ensure the offset is on the 1st of Januray of the current years so that we do not
983	! perturbe the solar angle calculations.
984	!
985	CALL ju2ymds (startdate, year, month, day, sec)
986	CALL ymds2ju (year, 1, 1, 0.0, solaryearstart)
987	!
988	last_time = -9999.0
989	!
990	ALLOCATE(mean_sinang(iim,jjm))
991	mean_sinang(:,:) = 0.0
992	!
993	split = NINT(forcing_tstep_ave/tlen)
994	!
995	! Allow for more space than estimated with the size of the first time step.
996	!
997	ALLOCATE(time_angles(split*2))
998	time_angles(:) = 0.0
999	!
1000	ALLOCATE(sinangles(iim,jjm,split*2))
1001	sinangles(:,:,:) = 0.0
1002	!
1003	dusk_angle=0.01
1004	!
1005	first_call_solarint = .FALSE.
1006	split = 0
1007	!
1008	ENDIF
1009	!
1010	! Shift the input dates in order to gain in precision for the time calculations
1011	!
1012	time_int(:) = time_int_in(:)-INT(solaryearstart)
1013	time_cent(1:slab_size) = time_cent_in(1:slab_size)-INT(solaryearstart)
1014	timebnd(1:slab_size,1) = timebnd_in(1:slab_size,1)-INT(solaryearstart)
1015	timebnd(1:slab_size,2) = timebnd_in(1:slab_size,2)-INT(solaryearstart)
1016	!
1017	! Create a mask so that MINLOC does not look outside of the valid interval of time_central
1018	!
1019	mask(1:slab_size) = .TRUE.
1020	!
1021	! Locate the time step in the SLAB at hand
1022	!
1023	imin = MINLOC( ABS(time_cent(1:slab_size)-(time_int(1)+time_int(2))/2.0), mask )
1024	!
1025	! Compute all the angels we will encounter for the current forcing interval
1026	!
1027	IF ( last_time .NE. timebnd(imin(1),1) ) THEN
1028	!
1029	! Verify that we have used all the angles of the previous decomposition of the forcing
1030	! time step.
1031	!
1032	IF ( split .NE. 0 ) THEN
1033	!
1034	WRITE(numout,*) "The forcing has a time step of : ", forcing_tstep_ave
1035	WRITE(numout,*) "The model is configured to run with a time step of : ", tlen
1036	WRITE(numout,*) "We are left with split = ", split, " starting from ", split_max
1037	!
1038	CALL ipslerr (3,'forcing_solarint',"The decomposition of solar downward radiation of the forcing file over the model",&
1039	& "has failed. This means the average of the solar radiation over the forcing time step is not conserved.",&
1040	& "This can be caused by a time step repeated twice.")
1041	ENDIF
1042	!
1043	! Compute the number of time steps the model will put in the current interval of forcing data.
1044	!
1045	split = 0
1046	julian_tmp = (time_int(1)+time_int(2))/2.0
1047	split_time = julian_tmp+split*tlen/one_day
1048	tmin = MINLOC( ABS(time_cent(1:slab_size) - split_time), mask)
1049	DO WHILE ( tmin(1) .EQ. imin(1) .AND. split_time .LE. timebnd(slab_size,2) )
1050	split = split + 1
1051	split_time = julian_tmp+split*tlen/one_day
1052	tmin = MINLOC( ABS(time_cent(1:slab_size) - split_time), mask)
1053	ENDDO
1054	!
1055	mean_sinang(:,:) = 0.0
1056	time_angles(:) = 0.0
1057	!
1058	DO is = 1,split
1059	!
1060	julian = julian_tmp + (is-1)*tlen/one_day
1061	!
1062	! This call should be better at it allows to compute the difference between the
1063	! current date and a reference time to higher precision. But it produces noisy
1064	! SWdown fluxes !
1065	!! CALL solarang (julian, solaryearstart, iim, jjm, lon, lat, sinang)
1066	! The old approach.
1067	CALL solarang (julian+INT(solaryearstart), solaryearstart, iim, jjm, lon, lat, sinang)
1068	!
1069	! During the dusk,dawn period maintain a minimum angle to take into account the
1070	! diffuse radiation which starts before the sun is over the horizon.
1071	!
1072	DO ii=1,iim
1073	DO jj=1,jjm
1074	IF ( sinang(ii,jj) > zero .AND. sinang(ii,jj) < dusk_angle ) THEN
1075	sinang(ii,jj) = dusk_angle
1076	ENDIF
1077	mean_sinang(ii,jj) = mean_sinang(ii,jj)+sinang(ii,jj)
1078	ENDDO
1079	ENDDO
1080	!
1081	! Save the solar angle information for later use. That is when we have the target time we will
1082	! look in this table the angle we have forseen.
1083	!
1084	time_angles(is) = julian
1085	sinangles(:,:,is) = sinang(:,:)
1086	!
1087	ENDDO
1088	!
1089	mean_sinang(:,:) = mean_sinang(:,:)/split
1090	last_time = timebnd(imin(1),1)
1091	split_max = split
1092	!
1093	ENDIF
1094	!
1095	! For the current timle step get the time of the closest pre-computed solar angle.
1096	!
1097	julian = (time_int(1)+time_int(2))/2.0
1098	tmin = MINLOC(ABS(julian-time_angles(1:split_max)), mask)
1099	sinang(:,:) = sinangles(:,:,tmin(1))
1100	! Remember that we have taken one value of the table for later verification
1101	split = split - 1
1102	!
1103	DO i=1, nbland_proc
1104	!
1105	jj = ((glolindex_proc(i)-1)/iim)+1
1106	ii = (glolindex_proc(i)-(jj-1)*iim)
1107	!
1108	IF ( mean_sinang(ii,jj) > zero ) THEN
1109	swdown(i) = swdown_slab(i,imin(1))*sinang(ii,jj)/mean_sinang(ii,jj)
1110	ELSE
1111	swdown(i) = zero
1112	ENDIF
1113	!
1114	! Why is this ??? Should we not take the angle corresponding to this time step ?? (Jan)
1115	!
1116	solarangle(i) = mean_sinang(ii,jj)
1117	!
1118	ENDDO
1119	!
1120	END SUBROUTINE forcing_solarint
1121	!!
1122	!! =============================================================================================================================
1123	!! SUBROUTINE: forcing_readslab
1124	!!
1125	!>\BRIEF Interface routine to read the data. This routine prepares the memory on each procesor and scatters the read data.
1126	!!
1127	!! DESCRIPTION:
1128	!!
1129	!! \n
1130	!_ ==============================================================================================================================
1131	SUBROUTINE forcing_readslab(time_int)
1132	!
1133	! This routine serves to interface with forcing_readslab_root and ensure that
1134	! the data is distributed correctly on all processors.
1135	!
1136	REAL(r_std), INTENT(in) :: time_int(2) !! The time interval over which the forcing is needed.
1137	!
1138	! Local
1139	!
1140	INTEGER(i_std) :: is
1141	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: tair_full, qair_full
1142	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: rainf_full, snowf_full
1143	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: swdown_full, lwdown_full
1144	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: u_full, v_full
1145	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: ps_full, ztq_full, zuv_full
1146	!
1147	! 1.0 Verify that for the slabs the memory is allocated for the variable
1148	! as well as its time axis.
1149	!
1150	IF ( .NOT. ALLOCATED(tair_slab) ) ALLOCATE(tair_slab(nbland_proc,slab_size))
1151	IF ( .NOT. ALLOCATED(time_tair) ) ALLOCATE(time_tair(slab_size))
1152	IF ( .NOT. ALLOCATED(timebnd_tair) ) ALLOCATE(timebnd_tair(slab_size,2))
1153	!
1154	IF ( .NOT. ALLOCATED(qair_slab) ) ALLOCATE(qair_slab(nbland_proc,slab_size))
1155	IF ( .NOT. ALLOCATED(time_qair) ) ALLOCATE(time_qair(slab_size))
1156	IF ( .NOT. ALLOCATED(timebnd_qair) ) ALLOCATE(timebnd_qair(slab_size,2))
1157	!
1158	IF ( .NOT. ALLOCATED(rainf_slab) ) ALLOCATE(rainf_slab(nbland_proc,slab_size))
1159	IF ( .NOT. ALLOCATED(snowf_slab) ) ALLOCATE(snowf_slab(nbland_proc,slab_size))
1160	IF ( .NOT. ALLOCATED(time_precip) ) ALLOCATE(time_precip(slab_size))
1161	IF ( .NOT. ALLOCATED(timebnd_precip) ) ALLOCATE(timebnd_precip(slab_size,2))
1162	IF ( .NOT. ALLOCATED(preciptime_slab) ) ALLOCATE(preciptime_slab(slab_size))
1163	!
1164	IF ( .NOT. ALLOCATED(swdown_slab) ) ALLOCATE(swdown_slab(nbland_proc,slab_size))
1165	IF ( .NOT. ALLOCATED(time_swdown) ) ALLOCATE(time_swdown(slab_size))
1166	IF ( .NOT. ALLOCATED(timebnd_swdown) ) ALLOCATE(timebnd_swdown(slab_size,2))
1167	!
1168	IF ( .NOT. ALLOCATED(lwdown_slab) ) ALLOCATE(lwdown_slab(nbland_proc,slab_size))
1169	IF ( .NOT. ALLOCATED(time_lwdown) ) ALLOCATE(time_lwdown(slab_size))
1170	IF ( .NOT. ALLOCATED(timebnd_lwdown) ) ALLOCATE(timebnd_lwdown(slab_size,2))
1171	!
1172	IF ( .NOT. ALLOCATED(u_slab) ) ALLOCATE(u_slab(nbland_proc,slab_size))
1173	IF ( .NOT. ALLOCATED(time_u) ) ALLOCATE(time_u(slab_size))
1174	IF ( .NOT. ALLOCATED(timebnd_u) ) ALLOCATE(timebnd_u(slab_size,2))
1175	!
1176	IF ( .NOT. ALLOCATED(v_slab) ) ALLOCATE(v_slab(nbland_proc,slab_size))
1177	IF ( .NOT. ALLOCATED(time_v) ) ALLOCATE(time_v(slab_size))
1178	IF ( .NOT. ALLOCATED(timebnd_v) ) ALLOCATE(timebnd_v(slab_size,2))
1179	!
1180	IF ( .NOT. ALLOCATED(ps_slab) ) ALLOCATE(ps_slab(nbland_proc,slab_size))
1181	IF ( .NOT. ALLOCATED(time_ps) ) ALLOCATE(time_ps(slab_size))
1182	IF ( .NOT. ALLOCATED(timebnd_ps) ) ALLOCATE(timebnd_ps(slab_size,2))
1183	!
1184	IF ( .NOT. ALLOCATED(ztq_slab) ) ALLOCATE(ztq_slab(nbland_proc,slab_size))
1185	IF ( .NOT. ALLOCATED(zuv_slab) ) ALLOCATE(zuv_slab(nbland_proc,slab_size))
1186	!
1187	!
1188	IF ( is_root_prc) THEN
1189	!
1190	! Allocate ont he root processor the memory to receive the data from the file
1191	!
1192	ALLOCATE(tair_full(nbland_loc,slab_size))
1193	ALLOCATE(qair_full(nbland_loc,slab_size))
1194	ALLOCATE(rainf_full(nbland_loc,slab_size))
1195	ALLOCATE(snowf_full(nbland_loc,slab_size))
1196	ALLOCATE(swdown_full(nbland_loc,slab_size))
1197	ALLOCATE(lwdown_full(nbland_loc,slab_size))
1198	ALLOCATE(u_full(nbland_loc,slab_size))
1199	ALLOCATE(v_full(nbland_loc,slab_size))
1200	ALLOCATE(ps_full(nbland_loc,slab_size))
1201	ALLOCATE(ztq_full(nbland_loc,slab_size))
1202	ALLOCATE(zuv_full(nbland_loc,slab_size))
1203	!
1204	CALL forcing_readslab_root(time_int, &
1205	& tair_full, time_tair, timebnd_tair, &
1206	& qair_full, time_qair, timebnd_qair, &
1207	& rainf_full, snowf_full, time_precip, timebnd_precip, preciptime_slab, &
1208	& swdown_full, time_swdown, timebnd_swdown, &
1209	& lwdown_full, time_lwdown, timebnd_lwdown, &
1210	& u_full, time_u, timebnd_u, &
1211	& v_full, time_v, timebnd_v, &
1212	& ps_full, time_ps, timebnd_ps, &
1213	& ztq_full, zuv_full)
1214	!
1215	ELSE
1216	!
1217	ALLOCATE(tair_full(1,1))
1218	ALLOCATE(qair_full(1,1))
1219	ALLOCATE(rainf_full(1,1))
1220	ALLOCATE(snowf_full(1,1))
1221	ALLOCATE(swdown_full(1,1))
1222	ALLOCATE(lwdown_full(1,1))
1223	ALLOCATE(u_full(1,1))
1224	ALLOCATE(v_full(1,1))
1225	ALLOCATE(ps_full(1,1))
1226	ALLOCATE(ztq_full(1,1))
1227	ALLOCATE(zuv_full(1,1))
1228	!
1229	ENDIF
1230	!
1231	! Broadcast the time information to all procs.
1232	!
1233	CALL bcast(slab_size)
1234	CALL bcast(time_tair)
1235	CALL bcast(timebnd_tair)
1236	CALL bcast(time_qair)
1237	CALL bcast(timebnd_qair)
1238	CALL bcast(time_precip)
1239	CALL bcast(timebnd_precip)
1240	CALL bcast(preciptime_slab)
1241	CALL bcast(time_swdown)
1242	CALL bcast(timebnd_swdown)
1243	CALL bcast(time_lwdown)
1244	CALL bcast(timebnd_lwdown)
1245	CALL bcast(time_u)
1246	CALL bcast(timebnd_u)
1247	CALL bcast(time_v)
1248	CALL bcast(timebnd_v)
1249	CALL bcast(time_ps)
1250	CALL bcast(timebnd_ps)
1251	!
1252	! Scatter the slabs of data to all processors
1253	!
1254	CALL scatter(tair_full, tair_slab)
1255	CALL scatter(qair_full, qair_slab)
1256	CALL scatter(rainf_full, rainf_slab)
1257	CALL scatter(snowf_full, snowf_slab)
1258	CALL scatter(swdown_full, swdown_slab)
1259	CALL scatter(lwdown_full, lwdown_slab)
1260	CALL scatter(u_full, u_slab)
1261	CALL scatter(v_full, v_slab)
1262	CALL scatter(ps_full, ps_slab)
1263	CALL scatter(ztq_full, ztq_slab)
1264	CALL scatter(zuv_full, zuv_slab)
1265	!
1266	! Clean-up to free the memory on the root processor.
1267	!
1268	IF ( ALLOCATED(tair_full) ) DEALLOCATE(tair_full)
1269	IF ( ALLOCATED(qair_full) ) DEALLOCATE(qair_full)
1270	IF ( ALLOCATED(rainf_full) ) DEALLOCATE(rainf_full)
1271	IF ( ALLOCATED(snowf_full) ) DEALLOCATE(snowf_full)
1272	IF ( ALLOCATED(swdown_full) ) DEALLOCATE(swdown_full)
1273	IF ( ALLOCATED(lwdown_full) ) DEALLOCATE(lwdown_full)
1274	IF ( ALLOCATED(u_full) ) DEALLOCATE(u_full)
1275	IF ( ALLOCATED(v_full) ) DEALLOCATE(v_full)
1276	IF ( ALLOCATED(ps_full) ) DEALLOCATE(ps_full)
1277	IF ( ALLOCATED(ztq_full) ) DEALLOCATE(ztq_full)
1278	IF ( ALLOCATED(zuv_full) ) DEALLOCATE(zuv_full)
1279	!
1280	END SUBROUTINE forcing_readslab
1281	!!
1282	!! =============================================================================================================================
1283	!! SUBROUTINE: forcing_readslab_root
1284	!!
1285	!>\BRIEF Routine which reads a slab of data from the netCDF file and writes it onto the memory.
1286	!!
1287	!! DESCRIPTION: It is important to read the next slab of data while still keeping an overlap so that
1288	!! interpolation can continue.
1289	!! It also attributes a time axis to each variable.
1290	!!
1291	!! \n
1292	!_ ==============================================================================================================================
1293
1294	SUBROUTINE forcing_readslab_root(time_int, &
1295	& tair, t_tair, tbnd_tair, &
1296	& qair, t_qair, tbnd_qair, &
1297	& rainf, snowf, t_prec, tbnd_prec, prectime, &
1298	& swdown, t_swdown, tbnd_swdown, &
1299	& lwdown, t_lwdown, tbnd_lwdown, &
1300	& u, t_u, tbnd_u, &
1301	& v, t_v, tbnd_v, &
1302	& ps, t_ps, tbnd_ps, &
1303	& ztq, zuv)
1304	!
1305	! Arguments
1306	!
1307	REAL(r_std), INTENT(in) :: time_int(2) !! The time interval over which the forcing is needed.
1308	!
1309	REAL(r_std), INTENT(out) :: tair(:,:)
1310	REAL(r_std), INTENT(out) :: t_tair(:)
1311	REAL(r_std), INTENT(out) :: tbnd_tair(:,:)
1312	!
1313	REAL(r_std), INTENT(out) :: qair(:,:)
1314	REAL(r_std), INTENT(out) :: t_qair(:)
1315	REAL(r_std), INTENT(out) :: tbnd_qair(:,:)
1316	!
1317	REAL(r_std), INTENT(out) :: rainf(:,:)
1318	REAL(r_std), INTENT(out) :: snowf(:,:)
1319	REAL(r_std), INTENT(out) :: t_prec(:)
1320	REAL(r_std), INTENT(out) :: tbnd_prec(:,:)
1321	REAL(r_std), INTENT(out) :: prectime(:)
1322	!
1323	REAL(r_std), INTENT(out) :: swdown(:,:)
1324	REAL(r_std), INTENT(out) :: t_swdown(:)
1325	REAL(r_std), INTENT(out) :: tbnd_swdown(:,:)
1326	!
1327	REAL(r_std), INTENT(out) :: lwdown(:,:)
1328	REAL(r_std), INTENT(out) :: t_lwdown(:)
1329	REAL(r_std), INTENT(out) :: tbnd_lwdown(:,:)
1330	!
1331	REAL(r_std), INTENT(out) :: u(:,:)
1332	REAL(r_std), INTENT(out) :: t_u(:)
1333	REAL(r_std), INTENT(out) :: tbnd_u(:,:)
1334	!
1335	REAL(r_std), INTENT(out) :: v(:,:)
1336	REAL(r_std), INTENT(out) :: t_v(:)
1337	REAL(r_std), INTENT(out) :: tbnd_v(:,:)
1338	!
1339	REAL(r_std), INTENT(out) :: ps(:,:)
1340	REAL(r_std), INTENT(out) :: t_ps(:)
1341	REAL(r_std), INTENT(out) :: tbnd_ps(:,:)
1342	!
1343	REAL(r_std), INTENT(out) :: ztq(:,:)
1344	REAL(r_std), INTENT(out) :: zuv(:,:)
1345	!
1346	! Local
1347	!
1348	INTEGER(i_std) :: iret, varid
1349	INTEGER(i_std) :: if, it
1350	INTEGER(i_std) :: tstart(3), tcount(3)
1351	INTEGER(i_std) :: imin(1), imax(1), firstif(1)
1352	INTEGER(i_std) :: nctstart, nctcount, inslabpos
1353	INTEGER(i_std) :: start_globtime, end_globtime
1354	INTEGER(i_std) :: timeid_tair, timeid_qair, timeid_precip, timeid_swdown
1355	INTEGER(i_std) :: timeid_lwdown, timeid_u, timeid_v, timeid_ps, timeid_tmp
1356	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: time_tmp
1357	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: rau
1358	CHARACTER(LEN=80) :: cellmethod
1359	!
1360	!
1361	ALLOCATE(time_tmp(slab_size,nbtax))
1362	ALLOCATE(rau(nbland_loc,slab_size))
1363	!
1364	! Catch any stupid utilisation of this routine.
1365	!
1366	IF ( .NOT. is_root_prc) THEN
1367	CALL ipslerr (3,'forcing_readslab_root',"Cannot run this routine o other procs than root.",&
1368	& "All the information on the forcing files is only lated on the root processor."," ")
1369	ENDIF
1370	!
1371	!Set some variables to zero
1372	!
1373	IF ( first_call_readslab ) THEN
1374	!
1375	preciptime(:) = 0
1376	!
1377	! If the first file is only there to provide the last time step of the previous year, we
1378	! do not need to read all. We will start 2 forcing time steps before the start of the first
1379	! time interval requested.
1380	!
1381	imin=MINLOC(ABS(time(:,1)-time_int(1)))
1382	current_offset = MAX(imin(1)-2,1)
1383	!
1384	first_call_readslab = .FALSE.
1385	!
1386	ELSE
1387	!
1388	! Put back the cummulated time of rainfall into the global array
1389	!
1390	preciptime(position_slab(1):position_slab(2)) = preciptime(position_slab(1):position_slab(2)) + &
1391	& prectime(1:slab_size)
1392	!
1393	! Compute new offset
1394	!
1395	current_offset = position_slab(2)-2
1396	!
1397	ENDIF
1398	!
1399	! Check that the slab size is not too large
1400	!
1401	IF ( (current_offset-1)+slab_size > nb_forcing_steps) THEN
1402	slab_size = nb_forcing_steps - (current_offset-1)
1403	ENDIF
1404	!
1405	! 1.1 Check that the slab we have to read exists
1406	!
1407	IF ( slab_size > 0 ) THEN
1408	!
1409	start_globtime = current_offset
1410	end_globtime = (current_offset-1)+slab_size
1411	inslabpos=1
1412	WRITE(,) ">> Reading from global position ", start_globtime, "up to ", end_globtime
1413	!
1414	DO IF=MINVAL(time_sourcefile(start_globtime:end_globtime)),MAXVAL(time_sourcefile(start_globtime:end_globtime))
1415	!
1416	! Get position of the part of the global time axis we need to read from this file
1417	!
1418	firstif = MINLOC(ABS(time_sourcefile-if))
1419	! start = distance of start_globtime or nothing + 1 to follow netCDF convention.
1420	nctstart = MAX((start_globtime-firstif(1)), 0)+1
1421	! count = end index - start index + 1
1422	nctcount = MIN((firstif(1)-1)+nbtime_perfile(if),end_globtime)-MAX(firstif(1),start_globtime)+1
1423	!
1424	!
1425	! Read time over the indexes computed above in order to position the slab correctly
1426	!
1427	WRITE(,) ">> From file ", if," reading from position ", nctstart, "up to ", (nctstart-1)+nctcount
1428	!
1429	DO it =1,nbtax
1430	tstart(1) = nctstart
1431	tcount(1) = nctcount
1432	iret = NF90_GET_VAR(force_id(if), time_id(if,it), time_tmp(inslabpos:inslabpos+nctcount-1,it), tstart, tcount)
1433	IF (iret /= NF90_NOERR) THEN
1434	WRITE(,) TRIM(NF90_STRERROR(iret))
1435	WRITE(,) "Working on file ", IF," starting at ",tstart(1)," counting ",tcount(1)
1436	WRITE(,) "The data was to be written in to section ", inslabpos,":",inslabpos+nctcount-1," of time_tmp"
1437	CALL ipslerr (3,'forcing_readslab_root',"Could not read the time for the current interval."," "," ")
1438	ENDIF
1439	time_tmp(inslabpos:inslabpos+nctcount-1,it) = date0_file(if,it) + &
1440	time_tmp(inslabpos:inslabpos+nctcount-1,it)*convtosec(if)/one_day
1441	ENDDO
1442	!
1443	! 2.0 Find and read variables.
1444	!
1445	! 2.1 Deal with air temperature and humidity as the fist and basic case
1446	!
1447	!
1448	!
1449	CALL forcing_varforslab(if, "Tair", nctstart, nctcount, inslabpos, tair, cellmethod)
1450	CALL forcing_attributetimeaxe(cellmethod, timeid_tair)
1451	!
1452	CALL forcing_varforslab(if, "Qair", nctstart, nctcount, inslabpos, qair, cellmethod)
1453	CALL forcing_attributetimeaxe(cellmethod, timeid_qair)
1454	!
1455	! 2.2 Deal with rainfall and snowfall.
1456	!
1457	CALL forcing_varforslab(if, "Rainf", nctstart, nctcount, inslabpos, rainf, cellmethod)
1458	CALL forcing_attributetimeaxe(cellmethod, timeid_precip)
1459	!
1460	CALL forcing_varforslab(if, "Snowf", nctstart, nctcount, inslabpos, snowf, cellmethod)
1461	CALL forcing_attributetimeaxe(cellmethod, timeid_tmp)
1462	IF ( timeid_precip .NE. timeid_tmp) THEN
1463	CALL ipslerr(3, 'forcing_readslab_root','Rainf and Snwof have different time axes.', &
1464	& 'Please check the forcing file to ensure both variable have the same cell_method.','')
1465	ENDIF
1466	!
1467	!
1468	! 2.3 Deal with downward shorwave and longwave radiation
1469	! The SW radiation can have 2 names SWdown_aerosol or SWdown. The first one is
1470	! given priority
1471	!
1472	CALL forcing_varforslab(if, "SWdown", nctstart, nctcount, inslabpos, swdown, cellmethod)
1473	CALL forcing_attributetimeaxe(cellmethod, timeid_swdown)
1474	!
1475	CALL forcing_varforslab(if, "LWdown", nctstart, nctcount, inslabpos, lwdown, cellmethod)
1476	CALL forcing_attributetimeaxe(cellmethod, timeid_lwdown)
1477	!
1478	!
1479	! 2.4 Deal with wind and pressure
1480	!
1481	CALL forcing_varforslab(if, "PSurf", nctstart, nctcount, inslabpos, ps, cellmethod)
1482	CALL forcing_attributetimeaxe(cellmethod, timeid_ps)
1483	!
1484	CALL forcing_varforslab(if, "Wind_E", nctstart, nctcount, inslabpos, u, cellmethod)
1485	CALL forcing_attributetimeaxe(cellmethod, timeid_u)
1486	!
1487	CALL forcing_varforslab(if, "Wind_N", nctstart, nctcount, inslabpos, v, cellmethod)
1488	CALL forcing_attributetimeaxe(cellmethod, timeid_v)
1489	!
1490	!
1491	! Do the height of the variables T&Q and U&V
1492	!
1493	! First the T&Q level
1494	!
1495	IF ( zheight ) THEN
1496	ztq(:,:) = zlev_fixed
1497	ELSE IF ( zsigma .OR. zhybrid ) THEN
1498	rau(:,:) = ps(:,:)/(cte_molr*tair(:,:))
1499	ztq(:,:) = (ps(:,:) - (zhybrid_a + zhybrid_bps(:,:)))/(rau(:,:) cte_grav)
1500	ELSE IF ( zlevels ) THEN
1501	CALL forcing_varforslab(IF, "Levels", nctstart, nctcount, inslabpos, ztq, cellmethod)
1502	ELSE
1503	CALL ipslerr(3, 'forcing_readslab_root','No case for the vertical levels was specified.', &
1504	& 'We cannot determine the height for T and Q.','')
1505	ENDIF
1506	!
1507	! Now the U&V level
1508	!
1509	IF ( zsamelev_uv ) THEN
1510	zuv(:,:) = ztq(:,:)
1511	ELSE
1512	IF ( zheight ) THEN
1513	zuv(:,:) = zlevuv_fixed
1514	ELSE IF ( zsigma .OR. zhybrid ) THEN
1515	rau(:,:) = ps(:,:)/(cte_molr*tair(:,:))
1516	zuv(:,:) = (ps(:,:) - (zhybriduv_a + zhybriduv_bps(:,:)))/(rau(:,:) cte_grav)
1517	ELSE IF ( zlevels ) THEN
1518	CALL forcing_varforslab(IF, "Levels_uv", nctstart, nctcount, inslabpos, zuv, cellmethod)
1519	ELSE
1520	CALL ipslerr(3, 'forcing_readslab_root','No case for the vertical levels was specified.', &
1521	& 'We cannot determine the height for U and V.','stop readdim2')
1522	ENDIF
1523	ENDIF
1524
1525	inslabpos = inslabpos+nctcount
1526
1527	ENDDO
1528	!
1529	! Use the read time of the slab to place it in the global variables. We consider
1530	! that we can do that on the first axis.
1531	!
1532	imin = MINLOC(ABS(time(:,1)-time_tmp(1,1)))
1533	position_slab(1) = imin(1)
1534	imax = MINLOC(ABS(time(:,1)-time_tmp(slab_size,1)))
1535	position_slab(2) = imax(1)
1536	!
1537	!
1538	IF ( position_slab(2)-position_slab(1) .GT. slab_size ) THEN
1539	WRITE(,) "Postition_slab =", position_slab
1540	WRITE(,) "Interval read : ", position_slab(2)-position_slab(1)
1541	WRITE(,) "Time start and end : ", time(1,1), time(slab_size,1)
1542	DO it =1,nbtax
1543	WRITE(,) "Checking time_tmp on idex : ", it
1544	WRITE(,) "Time_tmp start and end : ",time_tmp(1,it), time_tmp(slab_size,it)
1545	imin = MINLOC(ABS(time(:,1)-time_tmp(1,it)))
1546	imax = MINLOC(ABS(time(:,1)-time_tmp(slab_size,it)))
1547	WRITE(,) "Interval read : ", imax(1)-imin(1)+1
1548	ENDDO
1549	WRITE(,) "current_offset, slab_size =", current_offset, slab_size
1550	CALL ipslerr (3,'forcing_readslab_root',"The time slab read does not fit the number of variables read.",&
1551	& "Could there be an error in the time axis ?"," ")
1552	ENDIF
1553	!
1554	! Transfer the global time axis into the time variables approriate for each variable. This way
1555	! the time axis for each variable will be centered on the interval of validity. This is an essential assumption
1556	! the time interpolation to be done later.
1557	!
1558	WRITE(,) "We have found the following axes : ", time_axename(:)
1559	WRITE(,) "For Tair we are using time axis '",TRIM(time_axename(timeid_tair)),&
1560	& "' with cell method ",TRIM(time_cellmethod(timeid_tair)),"."
1561	t_tair(1:slab_size) = time(position_slab(1):position_slab(2), timeid_tair)
1562	tbnd_tair(1:slab_size,:) = time_bounds(position_slab(1):position_slab(2),timeid_tair,:)
1563	!
1564	WRITE(,) "For Qair we are using time axis '",TRIM(time_axename(timeid_qair)),&
1565	& "' with cell method ",TRIM(time_cellmethod(timeid_qair)),"."
1566	t_qair(1:slab_size) = time(position_slab(1):position_slab(2), timeid_qair)
1567	tbnd_qair(1:slab_size,:) = time_bounds(position_slab(1):position_slab(2),timeid_qair,:)
1568	!
1569	WRITE(,) "For Rainf and Snowf we are using time axis '",TRIM(time_axename(timeid_precip)),&
1570	& "' with cell method ",TRIM(time_cellmethod(timeid_precip)),"."
1571	t_prec(1:slab_size) = time(position_slab(1):position_slab(2), timeid_precip)
1572	tbnd_prec(1:slab_size,:) = time_bounds(position_slab(1):position_slab(2),timeid_precip,:)
1573	prectime(1:slab_size) = preciptime(position_slab(1):position_slab(2))
1574	!
1575	WRITE(,) "For SWdown we are using time axis '",TRIM(time_axename(timeid_swdown)),&
1576	& "' with cell method ",TRIM(time_cellmethod(timeid_swdown)),"."
1577	t_swdown(1:slab_size) = time(position_slab(1):position_slab(2), timeid_swdown)
1578	tbnd_swdown(1:slab_size,:) = time_bounds(position_slab(1):position_slab(2),timeid_swdown,:)
1579	!
1580	WRITE(,) "For LWdown we are using time axis '",TRIM(time_axename(timeid_lwdown)),&
1581	& "' with cell method ",TRIM(time_cellmethod(timeid_lwdown)),"."
1582	t_lwdown(1:slab_size) = time(position_slab(1):position_slab(2), timeid_lwdown)
1583	tbnd_lwdown(1:slab_size,:) = time_bounds(position_slab(1):position_slab(2),timeid_lwdown,:)
1584	!
1585	WRITE(,) "For Wind_E we are using time axis '",TRIM(time_axename(timeid_u)),&
1586	& "' with cell method ",TRIM(time_cellmethod(timeid_u)),"."
1587	t_u(1:slab_size) = time(position_slab(1):position_slab(2), timeid_u)
1588	tbnd_u(1:slab_size,:) = time_bounds(position_slab(1):position_slab(2),timeid_u,:)
1589	!
1590	WRITE(,) "For Wind_N we are using time axis '",TRIM(time_axename(timeid_v)),&
1591	& "' with cell method ",TRIM(time_cellmethod(timeid_v)),"."
1592	t_v(1:slab_size) = time(position_slab(1):position_slab(2), timeid_v)
1593	tbnd_v(1:slab_size,:) = time_bounds(position_slab(1):position_slab(2),timeid_v,:)
1594	!
1595	WRITE(,) "For PSurf we are using time axis '",TRIM(time_axename(timeid_ps)),&
1596	& "' with cell method ",TRIM(time_cellmethod(timeid_ps)),"."
1597	t_ps(1:slab_size) = time(position_slab(1):position_slab(2), timeid_ps)
1598	tbnd_ps(1:slab_size,:) = time_bounds(position_slab(1):position_slab(2),timeid_ps,:)
1599	!
1600	ELSE
1601	CALL ipslerr (2,'forcing_readslab_root',"We have reached the end of the slabs we can read.",&
1602	& "The calling program needs to manage this situation"," ")
1603	ENDIF
1604	!
1605	! Have we read to the end of the files ?
1606	!
1607	IF ( current_offset+slab_size >= nb_forcing_steps ) THEN
1608	end_of_file = .TRUE.
1609	ELSE
1610	end_of_file = .FALSE.
1611	ENDIF
1612	!
1613	IF ( ALLOCATED(rau) ) DEALLOCATE(rau)
1614	IF ( ALLOCATED(time_tmp) ) DEALLOCATE(time_tmp)
1615	!
1616	END SUBROUTINE forcing_readslab_root
1617
1618	!! =============================================================================================================================
1619	!! SUBROUTINE: forcing_reindex3d
1620	!!
1621	!>\BRIEF
1622	!!
1623	!! DESCRIPTION:
1624	!!
1625	!! \n
1626	!_ ==============================================================================================================================
1627	SUBROUTINE forcing_reindex3d(nbi, nbj, tin, slab_in, nbout, tout, slab_out, tstart, reindex)
1628	!
1629	! ARGUMENTS
1630	!
1631	INTEGER(i_std), INTENT(in) :: nbi, nbj, tin, nbout, tout
1632	REAL(r_std), INTENT(in) :: slab_in(nbi,nbj,tin)
1633	REAL(r_std), INTENT(out) :: slab_out(nbout,tout)
1634	INTEGER(i_std), INTENT(in) :: tstart
1635	INTEGER(i_std), INTENT(in) :: reindex(nbout,2)
1636	!
1637	! LOCAL
1638	!
1639	INTEGER(i_std) :: is, in
1640	!
1641	DO is=1,tin
1642	DO in=1,nbout
1643	slab_out(in,tstart+(is-1)) = slab_in(reindex(in,1),reindex(in,2),is)
1644	ENDDO
1645	ENDDO
1646	!
1647	END SUBROUTINE forcing_reindex3d
1648
1649	!! =============================================================================================================================
1650	!! SUBROUTINE: forcing_reindex2d
1651	!!
1652	!>\BRIEF
1653	!!
1654	!! DESCRIPTION:
1655	!!
1656	!! \n
1657	!_ ==============================================================================================================================
1658	SUBROUTINE forcing_reindex2d(nbi, nbj, slab_in, nbout, slab_out, reindex)
1659	!
1660	! ARGUMENTS
1661	!
1662	INTEGER(i_std), INTENT(in) :: nbi, nbj, nbout
1663	REAL(r_std), INTENT(in) :: slab_in(nbi,nbj)
1664	REAL(r_std), INTENT(out) :: slab_out(nbout)
1665	INTEGER(i_std), INTENT(in) :: reindex(nbout,2)
1666	!
1667	! LOCAL
1668	!
1669	INTEGER(i_std) :: in
1670	!
1671	DO in=1,nbout
1672	slab_out(in) = slab_in(reindex(in,1),reindex(in,2))
1673	ENDDO
1674	!
1675	END SUBROUTINE forcing_reindex2d
1676	!!
1677	!! =============================================================================================================================
1678	!! SUBROUTINE: forcing_reindex2dt
1679	!!
1680	!>\BRIEF
1681	!!
1682	!! DESCRIPTION:
1683	!!
1684	!! \n
1685	!_ ==============================================================================================================================
1686
1687	SUBROUTINE forcing_reindex2dt(nbin, tin, slab_in, nbout, tout, slab_out, tstart, reindex)
1688	!
1689	! ARGUMENTS
1690	!
1691	INTEGER(i_std), INTENT(in) :: nbin, tin, nbout, tout
1692	REAL(r_std), INTENT(in) :: slab_in(nbin,tin)
1693	REAL(r_std), INTENT(out) :: slab_out(nbout,tout)
1694	INTEGER(i_std), INTENT(in) :: tstart
1695	INTEGER(i_std), INTENT(in) :: reindex(nbout)
1696	!
1697	! LOCAL
1698	!
1699	INTEGER(i_std) :: is, in
1700	!
1701	DO is=1,tin
1702	DO in=1,nbout
1703	slab_out(in,tstart+(is-1)) = slab_in(reindex(in),is)
1704	ENDDO
1705	ENDDO
1706	!
1707	END SUBROUTINE forcing_reindex2dt
1708
1709	!! =============================================================================================================================
1710	!! SUBROUTINE: forcing_reindex1d
1711	!!
1712	!>\BRIEF
1713	!!
1714	!! DESCRIPTION:
1715	!!
1716	!! \n
1717	!_ ==============================================================================================================================
1718
1719	SUBROUTINE forcing_reindex1d(nbin, slab_in, nbout, slab_out, reindex)
1720	!
1721	! ARGUMENTS
1722	!
1723	INTEGER(i_std), INTENT(in) :: nbin, nbout
1724	REAL(r_std), INTENT(in) :: slab_in(nbin)
1725	REAL(r_std), INTENT(out) :: slab_out(nbout)
1726	INTEGER(i_std), INTENT(in) :: reindex(nbout)
1727	!
1728	! LOCAL
1729	!
1730	INTEGER(i_std) :: is
1731	!
1732	DO is=1,nbout
1733	slab_out(is) = slab_in(reindex(is))
1734	ENDDO
1735	!
1736	END SUBROUTINE forcing_reindex1d
1737	!!
1738	!! =============================================================================================================================
1739	!! SUBROUTINE: forcing_reindex2to1
1740	!!
1741	!>\BRIEF
1742	!!
1743	!! DESCRIPTION:
1744	!!
1745	!! \n
1746	!_ ==============================================================================================================================
1747
1748	SUBROUTINE forcing_reindex2to1(nbi, nbj, slab_in, nbout, slab_out, reindex)
1749	!
1750	! ARGUMENTS
1751	!
1752	INTEGER(i_std), INTENT(in) :: nbi, nbj, nbout
1753	REAL(r_std), INTENT(in) :: slab_in(nbi,nbj)
1754	REAL(r_std), INTENT(out) :: slab_out(nbout)
1755	INTEGER(i_std), INTENT(in) :: reindex(nbout)
1756	!
1757	! LOCAL
1758	!
1759	INTEGER(i_std) :: i, j, is
1760	!
1761	DO is=1,nbout
1762	j = INT((reindex(is)-1)/nbi)+1
1763	i = (reindex(is)-(j-1)*nbi)
1764	slab_out(is) = slab_in(i,j)
1765	ENDDO
1766	!
1767	END SUBROUTINE forcing_reindex2to1
1768
1769	!! =============================================================================================================================
1770	!! SUBROUTINE: forcing_reindex1to2
1771	!!
1772	!>\BRIEF
1773	!!
1774	!! DESCRIPTION:
1775	!!
1776	!! \n
1777	!_ ==============================================================================================================================
1778
1779	SUBROUTINE forcing_reindex1to2(nbin, slab_in, nbi, nbj, slab_out, reindex)
1780	!
1781	! ARGUMENTS
1782	!
1783	INTEGER(i_std), INTENT(in) :: nbin, nbi, nbj
1784	REAL(r_std), INTENT(in) :: slab_in(nbin)
1785	REAL(r_std), INTENT(out) :: slab_out(nbi, nbj)
1786	INTEGER(i_std), INTENT(in) :: reindex(nbin)
1787	!
1788	! LOCAL
1789	!
1790	INTEGER(i_std) :: i, j, is
1791	!
1792	DO is=1,nbin
1793	j = INT((reindex(is)-1)/nbi)+1
1794	i = (reindex(is)-(j-1)*nbi)
1795	slab_out(i,j) = slab_in(is)
1796	ENDDO
1797	!
1798	END SUBROUTINE forcing_reindex1to2
1799
1800	!! =============================================================================================================================
1801	!! SUBROUTINE: forcing_close
1802	!!
1803	!>\BRIEF Close all forcing files
1804	!!
1805	!! DESCRIPTION:
1806	!!
1807	!! \n
1808	!_ ==============================================================================================================================
1809	SUBROUTINE forcing_close()
1810
1811	INTEGER(i_std) :: ierr, if
1812
1813	DO if=1,nb_forcefile
1814	ierr = NF90_CLOSE(force_id(if))
1815	ENDDO
1816
1817	END SUBROUTINE forcing_close
1818
1819	!! =============================================================================================================================
1820	!! SUBROUTINE: forcing_printdate
1821	!!
1822	!>\BRIEF Print the date in a human readable format.
1823	!!
1824	!! DESCRIPTION:
1825	!!
1826	!! \n
1827	!_ ==============================================================================================================================
1828
1829	SUBROUTINE forcing_printdate(julian_day, message, wunit)
1830	!
1831	REAL(r_std), INTENT(in) :: julian_day
1832	CHARACTER(len=*), INTENT(in) :: message
1833	INTEGER, INTENT(in), OPTIONAL :: wunit
1834	!
1835	!
1836	!
1837	INTEGER(i_std) :: year, month, day, hours, minutes, seci
1838	REAL(r_std) :: sec
1839	!
1840	CALL ju2ymds (julian_day, year, month, day, sec)
1841	hours = INT(sec/3600)
1842	sec = sec - 3600 * hours
1843	minutes = INT(sec / 60)
1844	sec = sec - 60 * minutes
1845	seci = INT(sec)
1846	!
1847	IF (PRESENT(wunit)) THEN
1848	WRITE(wunit,'(I4.4,"-",I2.2,"-",I2.2," ",I2.2,":",I2.2,":",I2.2," > ", A60)') &
1849	& year, month, day, hours, minutes, seci, message
1850	ELSE
1851	WRITE(*,'(I4.4,"-",I2.2,"-",I2.2," ",I2.2,":",I2.2,":",I2.2," > ", A60)') &
1852	& year, month, day, hours, minutes, seci, message
1853	ENDIF
1854	!
1855	END SUBROUTINE forcing_printdate
1856
1857	!! =============================================================================================================================
1858	!! SUBROUTINE: forcing_printpoint_forgrid
1859	!!
1860	!>\BRIEF Together with the date print some sample values. Useful for checking the forcing.
1861	!!
1862	!! DESCRIPTION:
1863	!!
1864	!! \n
1865	!_ ==============================================================================================================================
1866
1867	SUBROUTINE forcing_printpoint_forgrid(julian_day, lon_pt, lat_pt, var, message)
1868	!
1869	REAL(r_std), INTENT(in) :: julian_day
1870	REAL(r_std), INTENT(in) :: lon_pt, lat_pt
1871	REAL(r_std), INTENT(in) :: var(:)
1872	CHARACTER(len=*), INTENT(in) :: message
1873	!
1874	!
1875	!
1876	INTEGER(i_std) :: year, month, day, hours, minutes, seci
1877	REAL(r_std) :: sec
1878	INTEGER(i_std) :: lon_ind, lat_ind, ind
1879	INTEGER(i_std), DIMENSION(1) :: i,j,k
1880	!
1881	! Check if there is anything to be done
1882	!
1883	IF ( MAX(lon_pt, lat_pt) > 360.0 ) THEN
1884	RETURN
1885	ENDIF
1886	!
1887	! Convert time first
1888	!
1889	CALL ju2ymds (julian_day, year, month, day, sec)
1890	hours = INT(sec/3600)
1891	sec = sec - 3600 * hours
1892	minutes = INT(sec / 60)
1893	sec = sec - 60 * minutes
1894	seci = INT(sec)
1895	!
1896	! Get the local to be analysed
1897	!
1898	i = MINLOC(ABS(lon_loc(:,1)-lon_pt))
1899	j = MINLOC(ABS(lat_loc(1,:)-lat_pt))
1900	ind = (j(1)-1)*iim_loc+i(1)
1901	k = MINLOC(ABS(lindex_loc(:)-ind))
1902	!
1903	WRITE(*,"(I2.2,':',I2.2,':',I2.2,' Loc : ', F5.1,',', F5.1,'(i=',I6,') Value = ',F12.4,A40)") &
1904	& hours, minutes, seci, lon_loc(i(1),1), lat_loc(1,j(1)), k(1), var(k(1)), message
1905
1906	END SUBROUTINE forcing_printpoint_forgrid
1907
1908	!! =============================================================================================================================
1909	!! SUBROUTINE: forcing_printpoint_gen
1910	!!
1911	!>\BRIEF Together with the date print some sample values. Useful for checking the forcing.
1912	!!
1913	!! DESCRIPTION:
1914	!!
1915	!! \n
1916	!_ ==============================================================================================================================
1917
1918	SUBROUTINE forcing_printpoint_gen(julian_day, lon_pt, lat_pt, nbind, lalo_in, var, message, ktest)
1919	!
1920	REAL(r_std), INTENT(in) :: julian_day
1921	REAL(r_std), INTENT(in) :: lon_pt, lat_pt
1922	INTEGER(i_std), INTENT(in) :: nbind
1923	REAL(r_std), INTENT(in) :: lalo_in(:,:)
1924	REAL(r_std), INTENT(in) :: var(:)
1925	CHARACTER(len=*), INTENT(in) :: message
1926	INTEGER(i_std), OPTIONAL, INTENT(out) :: ktest
1927	!
1928	!
1929	!
1930	INTEGER(i_std) :: year, month, day, hours, minutes, seci
1931	REAL(r_std) :: sec, mindist
1932	REAL(r_std), ALLOCATABLE, DIMENSION(:) :: dist, refdist
1933	INTEGER(i_std) :: lon_ind, lat_ind, ind
1934	INTEGER(i_std) :: i, imin(1)
1935	REAL(r_std), PARAMETER :: mincos = 0.0001
1936	REAL(r_std), PARAMETER :: pi = 3.141592653589793238
1937	REAL(r_std), PARAMETER :: R_Earth = 6378000.
1938	!
1939	! Check if there is anything to be done
1940	!
1941	IF ( MAX(lon_pt, lat_pt) > 360.0 ) THEN
1942	IF ( PRESENT(ktest) ) ktest = -1
1943	RETURN
1944	ENDIF
1945	!
1946	! Allocate memory
1947	!
1948	ALLOCATE(dist(nbind))
1949	ALLOCATE(refdist(nbind))
1950	!
1951	! Convert time first
1952	!
1953	CALL ju2ymds (julian_day, year, month, day, sec)
1954	hours = INT(sec/3600)
1955	sec = sec - 3600 * hours
1956	minutes = INT(sec / 60)
1957	sec = sec - 60 * minutes
1958	seci = INT(sec)
1959	!
1960	! Get the location to be analysed
1961	!
1962	DO i=1,nbind
1963	dist(i) = acos( sin(lat_ptpi/180)sin(lalo_in(i,1)*pi/180) + &
1964	& cos(lat_ptpi/180)cos(lalo_in(i,1)pi/180)&
1965	& cos((lalo_in(i,2)-lon_pt)pi/180) ) R_Earth
1966	ENDDO
1967	!
1968	! Look for the next grid point closest to the one with the smalest distance.
1969	!
1970	imin = MINLOC(dist)
1971	DO i=1,nbind
1972	refdist(i) = acos( sin(lalo_in(imin(1),1)pi/180)sin(lalo_in(i,1)*pi/180) + &
1973	& cos(lalo_in(imin(1),1)pi/180)cos(lalo_in(i,1)pi/180) &
1974	& cos((lalo_in(i,2)-lalo_in(imin(1),2))pi/180) ) R_Earth
1975	ENDDO
1976	refdist(imin(1)) = MAXVAL(refdist)
1977	mindist = MINVAL(refdist)
1978	!
1979	! Are we closer than the closest points ?
1980	!
1981	IF ( PRESENT(ktest) ) ktest = -1
1982	IF ( dist(imin(1)) <= mindist ) THEN
1983	!
1984	WRITE(*,"(I2.2,':',I2.2,':',I2.2,' Loc : ', F6.1,',', F6.1,'(i=',I6,') Value = ',F12.4,A38)") &
1985	& hours, minutes, seci, lalo_in(imin(1),2), lalo_in(imin(1),1), imin(1), var(imin(1)), message
1986	!
1987	IF ( PRESENT(ktest) ) ktest = imin(1)
1988	ENDIF
1989	!
1990	END SUBROUTINE forcing_printpoint_gen
1991
1992	!! =============================================================================================================================
1993	!! SUBROUTINE: forcing_getglogrid
1994	!!
1995	!>\BRIEF Routine to read the full grid of the forcing file.
1996	!!
1997	!! DESCRIPTION: The data is stored in the saved variables of the module and serve all other routines.
1998	!!
1999	!! \n
2000	!_ ==============================================================================================================================
2001
2002	SUBROUTINE forcing_getglogrid (nbfiles, filename, iim_tmp, jjm_tmp, nbland_tmp, closefile)
2003	!
2004	! This routine reads the global grid information from the forcing file and generates all the
2005	! information needed for this grid.
2006	!
2007	! ARGUMENTS
2008	!
2009	INTEGER(i_std), INTENT(in) :: nbfiles
2010	CHARACTER(LEN=*), INTENT(in) :: filename(:)
2011	INTEGER(i_std), INTENT(out) :: iim_tmp, jjm_tmp, nbland_tmp
2012	LOGICAL, INTENT(in) :: closefile
2013	!
2014	! LOCAL
2015	!
2016	INTEGER(i_std) :: iret, iv, if, lll
2017	CHARACTER(LEN=20) :: dimname, varname
2018	CHARACTER(LEN=60) :: lon_units, lat_units, units
2019	INTEGER(i_std), ALLOCATABLE, DIMENSION(:) :: dimids, londim_ids, latdim_ids
2020	INTEGER(i_std) :: lon_id, lat_id, land_id, lon_nbdims, lat_nbdims, land_nbdims
2021	INTEGER(i_std) :: lonvar_id, latvar_id, landvar_id
2022	LOGICAL :: dump_mask
2023	INTEGER(i_std) :: ik, i, j, iff, ndimsvar
2024	! Read a test variabe
2025	CHARACTER(len=8) :: testvarname='tair'
2026	INTEGER(i_std) :: testvar_id, contfrac_id
2027	REAL(r_std) :: testvar_missing, contfrac_missing
2028	REAL(r_std), ALLOCATABLE, DIMENSION(:) :: testvar
2029	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: testvar2d, contfrac2d
2030	!
2031	! 0.0 Check variables are allocated
2032	!
2033	IF ( .NOT. ALLOCATED(force_id)) ALLOCATE(force_id(nbfiles))
2034	IF ( .NOT. ALLOCATED(id_unlim)) ALLOCATE(id_unlim(nbfiles))
2035	IF ( .NOT. ALLOCATED(nb_atts)) ALLOCATE(nb_atts(nbfiles))
2036	IF ( .NOT. ALLOCATED(ndims)) ALLOCATE(ndims(nbfiles))
2037	IF ( .NOT. ALLOCATED(nvars)) ALLOCATE( nvars(nbfiles))
2038	!
2039	! 0.1 Are we one the root proc ?
2040	!
2041	IF ( .NOT. is_root_prc ) THEN
2042	CALL ipslerr (3,'forcing_getglogrid'," This routine can only be called on the root processor.", " ", " ")
2043	ENDIF
2044	!
2045	! 1.0 Open the netCDF file and get basic dimensions
2046	!
2047	DO iff=1,nbfiles
2048	!
2049	iret = NF90_OPEN(filename(iff), NF90_NOWRITE, force_id(iff))
2050	IF (iret /= NF90_NOERR) THEN
2051	CALL ipslerr (3,'forcing_getglogrid',"Error opening the forcing file :", filename(iff), " ")
2052	ENDIF
2053	!
2054	iret = NF90_INQUIRE (force_id(iff), nDimensions=ndims(iff), nVariables=nvars(iff), &
2055	nAttributes=nb_atts(iff), unlimitedDimId=id_unlim(iff))
2056	!
2057	!
2058	! 2.0 Read the dimensions found in the forcing file. Only deal with the spatial dimension as
2059	! time is an unlimited dimension.
2060	!
2061	DO iv=1,ndims(iff)
2062	!
2063	iret = NF90_INQUIRE_DIMENSION (force_id(iff), iv, name=dimname, len=lll)
2064	IF (iret /= NF90_NOERR) THEN
2065	CALL ipslerr (3,'forcing_getglogrid',"Could not get size of dimensions in file : ", filename(iff), " ")
2066	ENDIF
2067	!
2068	SELECT CASE(dimname)
2069	!
2070	CASE("west_east")
2071	CALL forcing_checkdim(iff, filename, iim_glo, lon_id, lll, iv)
2072	CASE("south_north")
2073	CALL forcing_checkdim(iff, filename, jjm_glo, lat_id, lll, iv)
2074	CASE("lon")
2075	CALL forcing_checkdim(iff, filename, iim_glo, lon_id, lll, iv)
2076	CASE("lat")
2077	CALL forcing_checkdim(iff, filename, jjm_glo, lat_id, lll, iv)
2078	CASE("nav_lon")
2079	CALL forcing_checkdim(iff, filename, iim_glo, lon_id, lll, iv)
2080	CASE("nav_lat")
2081	CALL forcing_checkdim(iff, filename, jjm_glo, lat_id, lll, iv)
2082	CASE("x")
2083	CALL forcing_checkdim(iff, filename, iim_glo, lon_id, lll, iv)
2084	CASE("y")
2085	CALL forcing_checkdim(iff, filename, jjm_glo, lat_id, lll, iv)
2086	CASE("land")
2087	CALL forcing_checkdim(iff, filename, nbland_glo, land_id, lll, iv)
2088	END SELECT
2089	!
2090	ENDDO
2091	ENDDO
2092	!
2093	! 3.0 Read the spatial coordinate variables found in the first file.
2094	!
2095	ALLOCATE(dimids(NF90_MAX_VAR_DIMS), londim_ids(NF90_MAX_VAR_DIMS), latdim_ids(NF90_MAX_VAR_DIMS))
2096	lonvar_id = -1
2097	latvar_id = -1
2098	landvar_id = -1
2099	testvar_id = -1
2100	contfrac_id = -1
2101	! Count the number of time axis we have
2102	nbtax = 0
2103	!
2104	DO iv=1,nvars(1)
2105	!
2106	iret = NF90_INQUIRE_VARIABLE(force_id(1), iv, name=varname, ndims=ndimsvar, dimids=dimids)
2107	iret = NF90_GET_ATT(force_id(1), iv, 'units', units)
2108	!
2109	! Check that we have the longitude
2110	!
2111	IF ( INDEX(lowercase(varname), 'lon') > 0 .AND. INDEX(lowercase(units), 'east') > 0 ) THEN
2112	lonvar_id = iv
2113	lon_units=units
2114	lon_nbdims = ndimsvar
2115	londim_ids = dimids
2116	ENDIF
2117	!
2118	! Check that we have the latitude
2119	!
2120	IF ( INDEX(lowercase(varname), 'lat') > 0 .AND. INDEX(lowercase(units), 'north') > 0) THEN
2121	latvar_id = iv
2122	lat_units=units
2123	lat_nbdims = ndimsvar
2124	latdim_ids = dimids
2125	ENDIF
2126	!
2127	! Check that we have the land re-indexing table
2128	!
2129	IF ( INDEX(lowercase(varname), 'land') > 0 ) THEN
2130	landvar_id = iv
2131	land_nbdims = ndimsvar
2132	latdim_ids = dimids
2133	ENDIF
2134	!
2135	! Check if we have the contfrac variable
2136	!
2137	IF ( INDEX(lowercase(varname), 'contfrac') > 0 ) THEN
2138	contfrac_id = iv
2139	iret = NF90_GET_ATT(force_id(1), iv, "missing_value", contfrac_missing)
2140	IF (iret /= NF90_NOERR) THEN
2141	contfrac_missing=0.0
2142	ENDIF
2143	ENDIF
2144	!
2145	! Find the test variable
2146	!
2147	IF ( INDEX(lowercase(varname), TRIM(testvarname)) > 0 ) THEN
2148	testvar_id = iv
2149	iret = NF90_GET_ATT(force_id(1), iv, "missing_value", testvar_missing)
2150	IF (iret /= NF90_NOERR) THEN
2151	testvar_missing=-1
2152	ENDIF
2153	ENDIF
2154	!
2155	! If we come across time get the calendar and archive it.
2156	!
2157	IF ( INDEX(lowercase(units),'seconds since') > 0 .OR. &
2158	& INDEX(lowercase(units),'minutes since') > 0 .OR. &
2159	& INDEX(lowercase(units),'hours since') > 0) THEN
2160	!
2161	! Get calendar used for the time axis
2162	!
2163	iret = NF90_GET_ATT(force_id(1), iv, "calendar", calendar)
2164	IF (iret == NF90_NOERR) THEN
2165	WRITE(,) ">> Setting the calendar to ",calendar
2166	ELSE
2167	WRITE(,) ">> Keeping proleptic Gregorian calendar"
2168	calendar="proleptic_gregorian"
2169	ENDIF
2170	!
2171	nbtax = nbtax + 1
2172	!
2173	ENDIF
2174	ENDDO
2175	!
2176	! 4.0 Verification that we have found both coordinate variables and the land point index
2177	!
2178	IF ( ( lonvar_id < 0 ) .AND. ( INDEX(lowercase(lon_units), 'east') <= 0 ) ) THEN
2179	CALL ipslerr (3,'forcing_getglogrid',"Have not found a valid longitude. Units = ", lon_units, " ")
2180	ENDIF
2181	IF ( ( latvar_id < 0 ) .AND. ( INDEX(lowercase(lat_units), 'north') <= 0 ) ) THEN
2182	CALL ipslerr (3,'forcing_getglogrid',"Have not found a valid latitude. Units = : ", lat_units, " ")
2183	ENDIF
2184	IF ( landvar_id < 0 ) THEN
2185	CALL ipslerr (1,'forcing_getglogrid',"No reindexing table was found. ", &
2186	& "This forcing file is not compressed by gathering.", " ")
2187	ENDIF
2188	!
2189	! 5.0 Allocate the space for the global variables and read them.
2190	!
2191	IF ( .NOT. ALLOCATED(lon_glo)) ALLOCATE(lon_glo(iim_glo, jjm_glo))
2192	IF ( .NOT. ALLOCATED(lat_glo)) ALLOCATE(lat_glo(iim_glo, jjm_glo))
2193	!
2194	IF ( lon_nbdims == 2 .AND. lat_nbdims == 2 ) THEN
2195	iret = NF90_GET_VAR(force_id(1), lonvar_id, lon_glo)
2196	iret = NF90_GET_VAR(force_id(1), latvar_id, lat_glo)
2197	ELSE IF ( lon_nbdims == 1 .AND. lat_nbdims == 1 ) THEN
2198	DO iv=1,jjm_glo
2199	iret = NF90_GET_VAR(force_id(1), lonvar_id, lon_glo(:,iv))
2200	ENDDO
2201	DO iv=1,iim_glo
2202	iret = NF90_GET_VAR(force_id(1), latvar_id, lat_glo(iv,:))
2203	ENDDO
2204	ELSE
2205	WRITE(,) "Found dimensions for lon and lat :", lon_nbdims, lat_nbdims
2206	CALL ipslerr (3,'forcing_getglogrid',"Longitude and Latitude variables do not have the right dimensions.", " ", " ")
2207	ENDIF
2208	!
2209	! If we have a indexing variable then the data is compressed by gathering, else we have to construct it.
2210	!
2211	compressed = .FALSE.
2212	IF ( landvar_id > 0 ) THEN
2213	IF ( .NOT. ALLOCATED(lindex_glo)) ALLOCATE(lindex_glo(nbland_glo))
2214	iret = NF90_GET_VAR(force_id(1), landvar_id, lindex_glo)
2215	compressed = .TRUE.
2216	ENDIF
2217	!
2218	IF ( .NOT. ALLOCATED(mask_glo)) ALLOCATE(mask_glo(iim_glo, jjm_glo))
2219	!
2220	! Get the land/sea mask and contfrac based on a test variable, if contfract is not available. Else
2221	! we use the contfrac variable.
2222	!
2223	IF ( compressed ) THEN
2224	IF ( .NOT. ALLOCATED(contfrac_glo)) ALLOCATE(contfrac_glo(nbland_glo))
2225	IF ( .NOT. ALLOCATED(testvar)) ALLOCATE(testvar(nbland_glo))
2226	CALL forcing_contfrac1d(force_id(1), testvar_id, contfrac_id, testvar)
2227	ELSE
2228	IF ( .NOT. ALLOCATED(testvar2d)) ALLOCATE(testvar2d(iim_glo, jjm_glo))
2229	IF ( .NOT. ALLOCATED(contfrac2d)) ALLOCATE(contfrac2d(iim_glo, jjm_glo))
2230	CALL forcing_contfrac2d(force_id(1), testvar_id, contfrac_id, testvar2d, contfrac2d, &
2231	& testvar_missing, contfrac_missing, nbland_glo)
2232	!
2233	! We have found a variable on which we can count the number of land points. We can build
2234	! the indexing table and gather the information needed.
2235	!
2236	IF ( .NOT. ALLOCATED(lindex_glo)) ALLOCATE(lindex_glo(nbland_glo))
2237	IF ( .NOT. ALLOCATED(contfrac_glo)) ALLOCATE(contfrac_glo(nbland_glo))
2238	IF ( .NOT. ALLOCATED(testvar)) ALLOCATE(testvar(nbland_glo))
2239	IF ( contfrac_id > 0 ) THEN
2240	CALL forcing_buildindex(contfrac2d, contfrac_missing, lindex_glo, contfrac_glo)
2241	CALL forcing_reindex(iim_glo, jjm_glo, testvar2d, nbland_glo, testvar, lindex_glo)
2242	ELSE
2243	CALL forcing_buildindex(testvar2d, testvar_missing, lindex_glo, testvar)
2244	contfrac_glo(:) = 1.0
2245	ENDIF
2246	!
2247	ENDIF
2248	!
2249	! We assume that if the forcing file is closed at the end of this subroutine this is a test
2250	! or initialisation of the grids. So we will dump the mask in a netCDF file for the user to
2251	! check.
2252	!
2253	dump_mask = closefile
2254	CALL forcing_checkindex(dump_mask, testvarname, testvar)
2255	!
2256	!
2257	! 8.0 Close file if needed
2258	!
2259	IF ( closefile ) THEN
2260	CALL forcing_close()
2261	ENDIF
2262	!
2263	! Prepare variables to be returnned to calling subroutine.
2264	!
2265	iim_tmp = iim_glo
2266	jjm_tmp = jjm_glo
2267	nbland_tmp = nbland_glo
2268	!
2269	! Clean up !
2270	!
2271	IF ( ALLOCATED(dimids) ) DEALLOCATE(dimids)
2272	IF ( ALLOCATED(londim_ids) ) DEALLOCATE(londim_ids)
2273	IF ( ALLOCATED(latdim_ids) ) DEALLOCATE(latdim_ids)
2274	IF ( ALLOCATED(testvar) ) DEALLOCATE(testvar)
2275	IF ( ALLOCATED(testvar2d) ) DEALLOCATE(testvar2d)
2276	IF ( ALLOCATED(contfrac2d) ) DEALLOCATE(contfrac2d)
2277	!
2278	END SUBROUTINE forcing_getglogrid
2279
2280	!! =============================================================================================================================
2281	!! SUBROUTINE: forcing_buildindex
2282	!!
2283	!>\BRIEF
2284	!!
2285	!! DESCRIPTION: When the forcing file does not contain compressed variables we need
2286	!! to build the land index variable from the mask defined by missing variables in
2287	!! a test variable.
2288	!!
2289	!! \n
2290	!_ ==============================================================================================================================
2291
2292	SUBROUTINE forcing_buildindex(var2d, var_missing, lindex, var)
2293	!
2294	! When the forcing file does not contain compressed variables we need
2295	! to build the land index variable from the mask defined by missing variables in
2296	! a test variable.
2297	!
2298	! Arguments
2299	!
2300	REAL(r_std), INTENT(in) :: var2d(:,:)
2301	REAL(r_std), INTENT(in) :: var_missing
2302	INTEGER(i_std), INTENT(out) :: lindex(:)
2303	REAL(r_std), INTENT(out) :: var(:)
2304	!
2305	! Local
2306	!
2307	INTEGER(i_std) :: i,j,k
2308	!
2309	k=0
2310	DO i=1,iim_glo
2311	DO j=1,jjm_glo
2312	IF (var2d(i,j) /= var_missing ) THEN
2313	k = k + 1
2314	lindex(k) = (j-1)*iim_glo+i
2315	var(k) = var2d(i,j)
2316	ENDIF
2317	ENDDO
2318	ENDDO
2319	!
2320	END SUBROUTINE forcing_buildindex
2321
2322	!! =============================================================================================================================
2323	!! SUBROUTINE: forcing_contfrac1d
2324	!!
2325	!>\BRIEF
2326	!!
2327	!! DESCRIPTION: This routine build the land/mask if needed and gets the contfrac variable from forcing file.
2328	!! Here we treat the case where the variables are compressed by gathering. Thus only
2329	!! land points are available in the file.
2330	!!
2331	!! \n
2332	!_ ==============================================================================================================================
2333
2334	SUBROUTINE forcing_contfrac1d(ifile, testvar_id, contfrac_id, testvar)
2335	!
2336	! This routine build the land/mask if needed and gets the contfrac variable from forcing file.
2337	! Here we treat the case where the variables are compressed by gathering. Thus only
2338	! land points are available in the file.
2339	!
2340	! ARGUMENTS
2341	!
2342	INTEGER(i_std), INTENT(in) :: ifile
2343	INTEGER(i_std), INTENT(in) :: testvar_id, contfrac_id
2344	REAL(r_std), DIMENSION(:), INTENT(inout) :: testvar
2345	!
2346	! LOCAL
2347	!
2348	INTEGER(i_std) :: it, iret
2349	INTEGER(i_std), DIMENSION(3) :: start, count
2350	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: contfrac2d
2351	!
2352	! First determine the contfrac variable
2353	!
2354	IF ( contfrac_id > 0 ) THEN
2355	iret = NF90_INQUIRE_VARIABLE(ifile, contfrac_id, ndims=it)
2356	IF ( it == 1 ) THEN
2357	start = (/1,1,0/)
2358	count = (/nbland_glo,1,0/)
2359	iret = NF90_GET_VAR(ifile, contfrac_id, contfrac_glo, start, count)
2360	IF (iret /= NF90_NOERR) THEN
2361	WRITE(,) TRIM(nf90_strerror(iret))
2362	CALL ipslerr (3,'forcing_contfrac1d',"Error reading contfrac ", " ", " ")
2363	ENDIF
2364	ELSE IF ( it == 2 ) THEN
2365	ALLOCATE(contfrac2d(iim_glo,jjm_glo))
2366	start = (/1,1,0/)
2367	count = (/iim_glo,jjm_glo,0/)
2368	iret = NF90_GET_VAR(ifile, contfrac_id, contfrac2d)
2369	IF (iret /= NF90_NOERR) THEN
2370	WRITE(,) TRIM(nf90_strerror(iret))
2371	CALL ipslerr (3,'forcing_contfrac1d',"Error reading contfrac ", " ", " ")
2372	ENDIF
2373	CALL forcing_reindex(iim_glo, jjm_glo, contfrac2d, nbland_glo, contfrac_glo, lindex_glo)
2374	DEALLOCATE(contfrac2d)
2375	ELSE
2376	CALL ipslerr (3,'forcing_contfrac1d',"Contfrac has a dimension larger than 2. ", &
2377	"We do not know how to handle this.", " ")
2378	ENDIF
2379	ELSE
2380	contfrac_glo(:) = 1.0
2381	ENDIF
2382	!
2383	! Read our test variable
2384	!
2385	iret = NF90_INQUIRE_VARIABLE(ifile, testvar_id, ndims=it)
2386	IF ( it == 2 ) THEN
2387	start = (/1,1,0/)
2388	count = (/nbland_glo,1,0/)
2389	ELSE IF ( it == 3 ) THEN
2390	start = (/1,1,1/)
2391	count = (/nbland_glo,1,1/)
2392	ELSE
2393	CALL ipslerr (3,'forcing_contfrac1d',"Testvar has a dimension larger than 3.", &
2394	"We do not know how to handle this", " ")
2395	ENDIF
2396	iret = NF90_GET_VAR(ifile, testvar_id, testvar, start, count)
2397	IF (iret /= NF90_NOERR) THEN
2398	WRITE(,) TRIM(nf90_strerror(iret))
2399	CALL ipslerr (3,'forcing_contfrac1d',"Error reading testvar.", " ", " ")
2400	ENDIF
2401	!
2402	END SUBROUTINE forcing_contfrac1d
2403
2404	!! =============================================================================================================================
2405	!! SUBROUTINE: forcing_contfrac2d
2406	!!
2407	!>\BRIEF
2408	!!
2409	!! DESCRIPTION: This routine build the land/mask if needed and gets the contfrac variable from forcing file.
2410	!! Here we treat the case where the variables is 2D. Thus we also need to identify the land points.
2411	!! For this we can either use the contfrac variable or the test variable.
2412	!!
2413	!! \n
2414	!_ ==============================================================================================================================
2415
2416	SUBROUTINE forcing_contfrac2d(ifile, testvar_id, contfrac_id, testvar, contfrac, testvar_missing, contfrac_missing, nbland)
2417	!
2418	! This routine build the land/mask if needed and gets the contfrac variable from forcing file.
2419	! Here we treat the case where the variables is 2D. Thus we also need to identify the land points.
2420	! For this we can either use the contfrac variable or the test variable.
2421	!
2422	! ARGUMENTS
2423	!
2424	INTEGER(i_std), INTENT(in) :: ifile
2425	INTEGER(i_std), INTENT(in) :: testvar_id, contfrac_id
2426	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: testvar
2427	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: contfrac
2428	REAL(r_std), INTENT(in) :: testvar_missing
2429	REAL(r_std), INTENT(in) :: contfrac_missing
2430	INTEGER(i_std), INTENT(out) :: nbland
2431	!
2432	! LOCAL
2433	!
2434	INTEGER(i_std) :: i, j, it, iret
2435	INTEGER(i_std), DIMENSION(4) :: start, count
2436	!
2437	!
2438	nbland = 0
2439	!
2440	IF ( contfrac_id > 0 ) THEN
2441	!
2442	iret = NF90_INQUIRE_VARIABLE(ifile, contfrac_id, ndims=it)
2443	IF ( it == 2 ) THEN
2444	start = (/1,1,0,0/)
2445	count = (/iim_glo,jjm_glo,0,0/)
2446	iret = NF90_GET_VAR(ifile, contfrac_id, contfrac, start, count)
2447	IF (iret /= NF90_NOERR) THEN
2448	WRITE(,) TRIM(nf90_strerror(iret))
2449	CALL ipslerr (3,'forcing_contfrac2d',"Error reading contfrac.", " ", " ")
2450	ENDIF
2451	ELSE
2452	CALL ipslerr (3,'forcing_contfrac2d',"Contfrac has a dimension different of 1.", &
2453	"We do not know how to handle this.", " ")
2454	ENDIF
2455	!
2456	! Count the number of land points.
2457	!
2458	DO i=1,iim_glo
2459	DO j=1,jjm_glo
2460	IF ( contfrac(i,j) /= contfrac_missing ) THEN
2461	nbland = nbland + 1
2462	ENDIF
2463	ENDDO
2464	ENDDO
2465	!
2466	! If we did not find any land points on the map (i.e. iim_glo > 1 and jjm_glo > 1) then we
2467	! look for fractions larger then 0.
2468	!
2469	IF ( iim_glo > 1 .AND. jjm_glo > 1 .AND. nbland < 1 ) THEN
2470	DO i=1,iim_glo
2471	DO j=1,jjm_glo
2472	IF ( contfrac(i,j) > 0.0 ) THEN
2473	nbland = nbland + 1
2474	ENDIF
2475	ENDDO
2476	ENDDO
2477	ENDIF
2478	!
2479	! Did we get a result ?
2480	!
2481	IF ( iim_glo > 1 .AND. jjm_glo > 1 .AND. nbland < 1 ) THEN
2482	CALL ipslerr (3,'forcing_contfrac2d',"Contfrac was used to count the number of land points.", &
2483	& "We still have not found any land points when we looked for contfrac > 0.", " ")
2484	ENDIF
2485	!
2486	ELSE
2487	! Just so that we have no un-initialized variable
2488	contfrac(:,:) = 0.0
2489	ENDIF
2490	!
2491	IF ( testvar_id > 0 ) THEN
2492	!
2493	iret = NF90_INQUIRE_VARIABLE(ifile, testvar_id, ndims=it)
2494	IF ( it == 2 ) THEN
2495	start = (/1,1,0,0/)
2496	count = (/iim_glo,jjm_glo,0,0/)
2497	ELSE IF ( it == 3 ) THEN
2498	start = (/1,1,1,0/)
2499	count = (/iim_glo,jjm_glo,1,0/)
2500	ELSE IF ( it == 4 ) THEN
2501	start = (/1,1,1,1/)
2502	count = (/iim_glo,jjm_glo,1,1/)
2503	ELSE
2504	CALL ipslerr (3,'forcing_contfrac2d',"testvar has a dimension of 1 or larger than 4.", &
2505	"We do not know how to handle this.", " ")
2506	ENDIF
2507	iret = NF90_GET_VAR(ifile, testvar_id, testvar, start, count)
2508	IF (iret /= NF90_NOERR) THEN
2509	WRITE(,) TRIM(nf90_strerror(iret))
2510	CALL ipslerr (3,'forcing_contfrac2d',"Error reading testvar.", " ", " ")
2511	ENDIF
2512	!
2513	! IF with count frac we did not get the number of land points, let us try it here
2514	!
2515	IF ( nbland < 1 ) THEN
2516	DO i=1,iim_glo
2517	DO j=1,jjm_glo
2518	IF ( testvar(i,j) /= testvar_missing ) THEN
2519	nbland = nbland + 1
2520	! Add infor to contfrac
2521	IF ( contfrac_id < 0 ) THEN
2522	contfrac(i,j) = 1.0
2523	ENDIF
2524	ENDIF
2525	ENDDO
2526	ENDDO
2527	ENDIF
2528	!
2529	!
2530	! Did we get a result here ?
2531	!
2532	IF ( iim_glo > 1 .AND. jjm_glo > 1 .AND. nbland < 1 ) THEN
2533	CALL ipslerr (3,'forcing_contfrac2d',"Contfrac and testvar were used to count the number", &
2534	& "of land points. We have not found any land points.", " ")
2535	ENDIF
2536	!
2537	ENDIF
2538	!
2539	END SUBROUTINE forcing_contfrac2d
2540
2541	!! =============================================================================================================================
2542	!! SUBROUTINE: forcing_checkindex
2543	!!
2544	!>\BRIEF
2545	!!
2546	!! DESCRIPTION: For ORCHIDEE its paralelisation requires that the land points are ordered
2547	!! in such a way that the longitude runs fastest. That means that we go over the
2548	!! globle filling one line after the other.
2549	!! As this might not be the case in a compressed vector of land points, we need to
2550	!! put all the points on the 2d grid and then scan them in the right order.
2551	!! The reindexing is prepared here.
2552	!!
2553	!! \n
2554	!_ ==============================================================================================================================
2555
2556	SUBROUTINE forcing_checkindex(dump_mask, testvarname, testvar)
2557	!
2558	! For ORCHIDEE its paralelisation requires that the land points are ordered
2559	! in such a way that the longitude runs fastest. That means that we go over the
2560	! globle filling one line after the other.
2561	! As this might not be the case in a compressed vector of land points, we need to
2562	! put all the points on the 2d grid and then scan them in the right order.
2563	! The reindexing is prepared here.
2564	!
2565	LOGICAL :: dump_mask
2566	CHARACTER(LEN=*) :: testvarname
2567	REAL(r_std) :: testvar(:)
2568	!
2569	INTEGER(i_std) :: j, i, ik
2570	REAL(r_std), ALLOCATABLE, DIMENSION(:) :: testvar_reind
2571	!
2572	!
2573	!
2574	! Get the indices of the land points in the focing file
2575	!
2576	IF ( .NOT. ALLOCATED(reindex_glo)) ALLOCATE(reindex_glo(nbland_glo))
2577	IF ( .NOT. ALLOCATED(origind)) ALLOCATE(origind(iim_glo,jjm_glo))
2578	!
2579	! Find the origine of each point in the gathered vector on the xy grid.
2580	!
2581	origind(:,:) = -1
2582	mask_glo(:,:) = 0
2583	DO ik=1,nbland_glo
2584	j = INT((lindex_glo(ik)-1)/iim_glo)+1
2585	i = (lindex_glo(ik)-(j-1)*iim_glo)
2586	origind(i,j) = ik
2587	mask_glo(i,j) = 1
2588	ENDDO
2589	!
2590	! Prepare a reindexing array so that the vector goes in the right order : longitude runs
2591	! faster than the latitude. Put then also the right information into lindex_glo.
2592	!
2593	ik=0
2594	DO j=1,jjm_glo
2595	DO i=1,iim_glo
2596	IF ( origind(i,j) > zero ) THEN
2597	ik = ik+1
2598	reindex_glo(ik) = origind(i,j)
2599	lindex_glo(ik) = (j-1)*iim_glo+i
2600	ENDIF
2601	ENDDO
2602	ENDDO
2603	!
2604	!
2605	! Write the mask and a test variable to a file so that the user can check that all is OK
2606	!
2607	IF ( dump_mask) THEN
2608	!
2609	! Scatter the test variable and save it in the file
2610	!
2611	WRITE(,) MINVAL(testvar), "<< test variable ", TRIM(testvarname), " <<", MAXVAL(testvar)
2612	ALLOCATE(testvar_reind(nbland_glo))
2613	!
2614	CALL forcing_reindex(nbland_glo, testvar, nbland_glo, testvar_reind, reindex_glo)
2615	!
2616	CALL forcing_writetestvar("forcing_mask_glo.nc", iim_glo, jjm_glo, nbland_glo, &
2617	& lon_glo(:,1), lat_glo(1,:), lindex_glo, mask_glo, &
2618	& testvarname, testvar_reind)
2619	!
2620	ENDIF
2621	!
2622	! Clean up !
2623	!
2624	IF ( ALLOCATED(testvar_reind) ) DEALLOCATE(testvar_reind)
2625	!
2626	END SUBROUTINE forcing_checkindex
2627
2628	!! =============================================================================================================================
2629	!! SUBROUTINE: forcing_writetestvar
2630	!!
2631	!>\BRIEF Write the mask and a test variable to a netCDF file.
2632	!!
2633	!! DESCRIPTION: This routine allows to test if the variables read from the forcing files is well read.
2634	!! Typically the forcing is compressed by gathering and thus it is a safe practice
2635	!! to verify that the un-compression is done correctly and that all points end-up in the
2636	!! right place on the global lat/lon grid.
2637	!!
2638	!! \n
2639	!_ ==============================================================================================================================
2640
2641	SUBROUTINE forcing_writetestvar(ncdffile, iim, jjm, nbland, lon, lat, lindex, mask, varname, var)
2642	!
2643	! Write the mask and a test variable to a netCDF file
2644	!
2645	! ARGUMENTS
2646	!
2647	CHARACTER(LEN=*), INTENT(in) :: ncdffile
2648	INTEGER(i_std), INTENT(in) :: iim, jjm, nbland
2649	REAL(r_std), INTENT(in) :: lon(iim), lat(jjm)
2650	INTEGER(i_std), INTENT(in) :: lindex(nbland)
2651	INTEGER(i_std), INTENT(in) :: mask(iim,jjm)
2652	CHARACTER(LEN=*), INTENT(in) :: varname
2653	REAL(r_std), INTENT(in) :: var(nbland)
2654	!
2655	! Local
2656	!
2657	INTEGER(i_std) :: ik, i, j
2658	INTEGER(i_std) :: iret, nlonid, nlatid, varid, fid, ierr, iland
2659	INTEGER(i_std) :: testid
2660	INTEGER(i_std), DIMENSION(2) :: lolaid
2661	REAL(r_std) :: test_scat(iim,jjm)
2662	!
2663	WRITE(,) "Working on file ", TRIM(ncdffile)
2664	!
2665	test_scat(:,:) = NF90_FILL_REAL
2666	CALL forcing_reindex(nbland, var, iim, jjm, test_scat, lindex)
2667	!
2668	iret = NF90_CREATE(ncdffile, NF90_WRITE, fid)
2669	IF (iret /= NF90_NOERR) THEN
2670	CALL ipslerr (3,'forcing_writetestvar',"Error opening the output file : ", ncdffile, " ")
2671	ENDIF
2672	!
2673	! Define dimensions
2674	!
2675	iret = NF90_DEF_DIM(fid,'lon',iim,lolaid(1))
2676	iret = NF90_DEF_DIM(fid,'lat',jjm,lolaid(2))
2677	!
2678	!
2679	iret = NF90_DEF_VAR(fid,"lon",NF90_REAL4,lolaid(1),nlonid)
2680	iret = NF90_PUT_ATT(fid,nlonid,'axis',"X")
2681	iret = NF90_PUT_ATT(fid,nlonid,'standard_name',"longitude")
2682	iret = NF90_PUT_ATT(fid,nlonid,'units',"degree_east")
2683	iret = NF90_PUT_ATT(fid,nlonid,'valid_min',MINVAL(lon_glo))
2684	iret = NF90_PUT_ATT(fid,nlonid,'valid_max',MAXVAL(lon_glo))
2685	iret = NF90_PUT_ATT(fid,nlonid,'long_name',"Longitude")
2686	!
2687	iret = NF90_DEF_VAR(fid,"lat",NF90_REAL4,lolaid(2),nlatid)
2688	iret = NF90_PUT_ATT(fid,nlatid,'axis',"Y")
2689	iret = NF90_PUT_ATT(fid,nlatid,'standard_name',"latitude")
2690	iret = NF90_PUT_ATT(fid,nlatid,'units',"degree_north")
2691	iret = NF90_PUT_ATT(fid,nlatid,'valid_min',MINVAL(lat_glo))
2692	iret = NF90_PUT_ATT(fid,nlatid,'valid_max',MAXVAL(lat_glo))
2693	iret = NF90_PUT_ATT(fid,nlatid,'long_name',"Latitude")
2694	!
2695	iret = NF90_DEF_VAR(fid,"mask",NF90_REAL4,lolaid,varid)
2696	!
2697	iret = NF90_DEF_VAR(fid,TRIM(varname),NF90_REAL4,lolaid,testid)
2698	iret = NF90_PUT_ATT(fid,testid,'_FillValue',NF90_FILL_REAL)
2699	iret = NF90_PUT_ATT(fid,testid,'missing_value',NF90_FILL_REAL)
2700	!
2701	iret = NF90_ENDDEF (fid)
2702	IF (iret /= NF90_NOERR) THEN
2703	WRITE(,) TRIM(nf90_strerror(iret))
2704	CALL ipslerr (3,'forcing_writetestvar',"Error ending definitions in file : ", ncdffile, " ")
2705	ENDIF
2706	!
2707	! Write variables
2708	!
2709	iret = NF90_PUT_VAR(fid,nlonid,lon)
2710	iret = NF90_PUT_VAR(fid,nlatid,lat)
2711	iret = NF90_PUT_VAR(fid,varid,REAL(mask))
2712	iret = NF90_PUT_VAR(fid,testid,test_scat)
2713	!
2714	! Close file
2715	!
2716	iret = NF90_CLOSE(fid)
2717	IF (iret /= NF90_NOERR) THEN
2718	CALL ipslerr (3,'forcing_writetestvar',"Error closing the output file : ", ncdffile, " ")
2719	ENDIF
2720	!
2721	END SUBROUTINE forcing_writetestvar
2722
2723	!! =============================================================================================================================
2724	!! SUBROUTINE: forcing_zoomgrid
2725	!!
2726	!>\BRIEF We zoom into the region requested by the user.
2727	!!
2728	!! DESCRIPTION: Get the area to be zoomed and the sizes of arrays we will need.
2729	!! This subroutine fills the *_loc variables.
2730	!! If requested it will dump a test vraible into a netCDF file.
2731	!!
2732	!! \n
2733	!_ ==============================================================================================================================
2734
2735	SUBROUTINE forcing_zoomgrid (zoom_lon, zoom_lat, filename, dumpncdf)
2736	!
2737	! Get the area to be zoomed and the sizes of arrays we will need.
2738	! This subroutine fills the *_loc variables.
2739	! If requested it will dump a test vraible into a netCDF file.
2740	!
2741	! ARGUMENTS
2742	!
2743	REAL(r_std), DIMENSION(2), INTENT(in) :: zoom_lon, zoom_lat
2744	CHARACTER(LEN=*), INTENT(in) :: filename
2745	LOGICAL, INTENT(in) :: dumpncdf
2746	!
2747	! LOCAL
2748	!
2749	INTEGER(i_std) :: i, j, ic, jc, ik, ig
2750	REAL(r_std) :: dx, dy, coslat
2751	REAL(r_std) :: lon_tmp(iim_glo), lat_tmp(jjm_glo)
2752	REAL(r_std) :: longlo_tmp(iim_glo,jjm_glo)
2753	REAL(r_std), ALLOCATABLE, DIMENSION(:) :: lon_val, lat_val
2754	INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:,:) :: zoom_index
2755	!
2756	INTEGER(i_std) :: iret, force_id, iv
2757	INTEGER(i_std), DIMENSION(1) :: imin, jmin
2758	INTEGER(i_std), DIMENSION(2) :: start, count
2759	INTEGER(i_std), DIMENSION(3) :: start2d, count2d
2760	REAL(r_std), ALLOCATABLE, DIMENSION(:) :: readvar, zoomedvar
2761	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: readvar2d
2762	INTEGER(i_std), ALLOCATABLE, DIMENSION(:) :: index_glotoloc
2763	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: lalo
2764	CHARACTER(LEN=8) :: testvarname="Tair"
2765	!
2766	! 0.0 Verify we are on the root processor
2767	!
2768	IF ( .NOT. is_root_prc ) THEN
2769	CALL ipslerr (3,'forcing_zoomgrid'," This routine can only be called on the root processor.", " ", " ")
2770	ENDIF
2771	!
2772	! 0.1 Inform the user
2773	!
2774	WRITE(,) "Zoom forcing : lon = ", zoom_lon
2775	WRITE(,) "Zoom forcing : lat = ", zoom_lat
2776	!
2777	! Some forcing files have longitudes going from 0 to 360. This code works on the
2778	! -180 to 180 longitude grid. So if needed we transform the longitudes of the global grid.
2779	!
2780	IF ( MAXVAL(lon_glo) <= 180.0 ) THEN
2781	longlo_tmp=lon_glo
2782	ELSE
2783	DO i=1,iim_glo
2784	DO j=1,jjm_glo
2785	IF ( lon_glo(i,j) <= 180.0 ) THEN
2786	longlo_tmp(i,j) = lon_glo(i,j)
2787	ELSE
2788	longlo_tmp(i,j) = lon_glo(i,j)-360
2789	ENDIF
2790	ENDDO
2791	ENDDO
2792	ENDIF
2793	!
2794	! See if we need to zoom
2795	!
2796	IF (MINVAL(zoom_lon) > MINVAL(longlo_tmp) .OR. MAXVAL(zoom_lon) < MAXVAL(longlo_tmp) .OR.&
2797	& MINVAL(zoom_lat) > MINVAL(lat_glo) .OR. MAXVAL(zoom_lat) < MAXVAL(lat_glo) ) THEN
2798	zoom_forcing = .TRUE.
2799	ENDIF
2800	!
2801	! Determine the size in x and y of the zoom
2802	!
2803	IF ( zoom_forcing ) THEN
2804	!
2805	! Working with the hypothesis it is a regular global grid and bring it back to the -180 to 180 interval
2806	! if needed.
2807	!
2808	lon_tmp(:) = longlo_tmp(:,1)
2809	lat_tmp(:) = lat_glo(1,:)
2810	!
2811	DO i=1,iim_glo
2812	IF ( lon_tmp(i) <= MINVAL(zoom_lon) .OR. lon_tmp(i) >= MAXVAL(zoom_lon) ) THEN
2813	lon_tmp(i) = 0.0
2814	ELSE
2815	lon_tmp(i) = 1.0
2816	ENDIF
2817	ENDDO
2818	DO j=1,jjm_glo
2819	IF ( lat_tmp(j) <= MINVAL(zoom_lat) .OR. lat_tmp(j) >= MAXVAL(zoom_lat) ) THEN
2820	lat_tmp(j) = 0.0
2821	ELSE
2822	lat_tmp(j) = 1.0
2823	ENDIF
2824	ENDDO
2825	iim_loc = NINT(SUM(lon_tmp))
2826	jjm_loc = NINT(SUM(lat_tmp))
2827	ELSE
2828	iim_loc = iim_glo
2829	jjm_loc = jjm_glo
2830	lon_tmp(:) = 1.0
2831	lat_tmp(:) = 1.0
2832	ENDIF
2833	!
2834	! Determine the number of land points in the zoom
2835	!
2836	IF ( .NOT. ALLOCATED(lon_loc) ) ALLOCATE(lon_loc(iim_loc,jjm_loc))
2837	IF ( .NOT. ALLOCATED(lat_loc) ) ALLOCATE(lat_loc(iim_loc,jjm_loc))
2838	IF ( .NOT. ALLOCATED(mask_loc) ) ALLOCATE(mask_loc(iim_loc,jjm_loc))
2839	IF ( .NOT. ALLOCATED(zoom_index) ) ALLOCATE(zoom_index(iim_loc,jjm_loc,2))
2840	!
2841	IF ( .NOT. ALLOCATED(lon_val)) ALLOCATE(lon_val(iim_loc))
2842	IF ( .NOT. ALLOCATED(lat_val)) ALLOCATE(lat_val(jjm_loc))
2843	!
2844	! Create our new lat/lon system which is in the -180/180 range and South to North and West to East
2845	!
2846	ic=0
2847	DO i=1,iim_glo
2848	IF ( lon_tmp(i) > 0 ) THEN
2849	ic = ic+1
2850	lon_val(ic) = longlo_tmp(i,1)
2851	ENDIF
2852	ENDDO
2853	jc=0
2854	DO j=1,jjm_glo
2855	IF ( lat_tmp(j) > 0 ) THEN
2856	jc = jc+1
2857	lat_val(jc) = lat_glo(1,j)
2858	ENDIF
2859	ENDDO
2860	CALL sort(lon_val, iim_loc)
2861	CALL sort(lat_val, jjm_loc)
2862	!
2863	! Now find the correspondance between the zoomed & re-ordered grid and the global one.
2864	!
2865	DO i=1,iim_loc
2866	DO j=1,jjm_loc
2867	!
2868	imin=MINLOC(ABS(longlo_tmp(:,1)-lon_val(i)))
2869	jmin=MINLOC(ABS(lat_glo(1,:)-lat_val(j)))
2870	!
2871	lon_loc(i,j) = longlo_tmp(imin(1),jmin(1))
2872	lat_loc(i,j) = lat_glo(imin(1),jmin(1))
2873	mask_loc(i,j) = mask_glo(imin(1),jmin(1))
2874	!
2875	zoom_index(i,j,1) = imin(1)
2876	zoom_index(i,j,2) = jmin(1)
2877	!
2878	ENDDO
2879	ENDDO
2880	!
2881	nbland_loc = SUM(mask_loc)
2882	IF ( .NOT. zoom_forcing .AND. nbland_loc .NE. nbland_glo) THEN
2883	WRITE(,) "We have not zoomed into the forcing file still we get a number of"
2884	WRITE(,) "land points that is different from what we have in the forcing file."
2885	STOP "forcing_zoomgrid"
2886	ENDIF
2887	!
2888	IF ( .NOT. ALLOCATED(lindex_loc)) ALLOCATE(lindex_loc(nbland_loc))
2889	IF ( .NOT. ALLOCATED(reindex_loc)) ALLOCATE(reindex_loc(nbland_loc))
2890	IF ( .NOT. ALLOCATED(contfrac_loc)) ALLOCATE(contfrac_loc(nbland_loc))
2891	!
2892	IF ( .NOT. ALLOCATED(reindex2d_loc)) ALLOCATE(reindex2d_loc(nbland_loc,2))
2893	IF ( .NOT. ALLOCATED(index_glotoloc)) ALLOCATE(index_glotoloc(nbland_glo))
2894	IF ( .NOT. ALLOCATED(lalo)) ALLOCATE(lalo(nbland_loc,2))
2895	!
2896	! Do the actual zoom on the grid
2897	!
2898	! Set indices of all points as non existant so that we can fill in as we zoom the
2899	! indices of the points which exist.
2900	index_glotoloc(:) = -1
2901	!
2902	ik = 0
2903	!
2904	! Loop only over the zoomed grid
2905	!
2906	! Why does the inner loop need to be ic for the pralalisation ????
2907	!
2908	DO jc=1,jjm_loc
2909	DO ic=1,iim_loc
2910	!
2911	! Point back from the local to the original global i*j grid
2912	!
2913	i = zoom_index(ic,jc,1)
2914	j = zoom_index(ic,jc,2)
2915	!
2916	IF ( origind(i,j) > 0 ) THEN
2917	ik = ik+1
2918	! index of the points in the local grid
2919	lindex_loc(ik) = (jc-1)*iim_loc+ic
2920	!
2921	! For land points, the index of global grid is saved for the this point on the local grid
2922	reindex_loc(ik) = origind(i,j)
2923	!
2924	! Keep also the i and j of the global grid for this land point on the local grid
2925	reindex2d_loc(ik,1) = i
2926	reindex2d_loc(ik,2) = j
2927	!
2928	! Keep the reverse : on the global grid the location where we put the value of the local grid.
2929	index_glotoloc(origind(i,j)) = ik
2930	!
2931	contfrac_loc(ik) = contfrac_glo(origind(i,j))
2932	!
2933	lalo(ik,1) = lat_glo(i,j)
2934	lalo(ik,2) = longlo_tmp(i,j)
2935	!
2936	ENDIF
2937	ENDDO
2938	ENDDO
2939	!
2940	!
2941	!
2942	ncdfstart = MINVAL(reindex_loc)
2943	reindex_loc(:) = reindex_loc(:)-ncdfstart+1
2944	ncdfcount = MAXVAL(reindex_loc)
2945	!
2946	! Compute the areas and the corners on the grid over which we will run ORCHIDEE.
2947	! As this module is dedicated for regular lat/lon forcing we know that we can compute these
2948	! variables without further worries.
2949	!
2950	IF ( .NOT. ALLOCATED(area_loc)) ALLOCATE(area_loc(iim_loc,jjm_loc))
2951	IF ( .NOT. ALLOCATED(corners_loc)) ALLOCATE(corners_loc(iim_loc,jjm_loc,4,2))
2952	!
2953	! Treat first the longitudes
2954	!
2955	DO j=1,jjm_loc
2956	dx = zero
2957	DO i=1,iim_loc-1
2958	dx = dx+ABS(lon_loc(i,j)-lon_loc(i+1,j))
2959	ENDDO
2960	dx = dx/(iim_loc-1)
2961	DO i=1,iim_loc
2962	corners_loc(i,j,1,1) = lon_loc(i,j)-dx/2.0
2963	corners_loc(i,j,2,1) = lon_loc(i,j)+dx/2.0
2964	corners_loc(i,j,3,1) = lon_loc(i,j)+dx/2.0
2965	corners_loc(i,j,4,1) = lon_loc(i,j)-dx/2.0
2966	ENDDO
2967	ENDDO
2968	!
2969	! Now treat the latitudes
2970	!
2971	DO i=1,iim_loc
2972	dy = zero
2973	DO j=1,jjm_loc-1
2974	dy = dy + ABS(lat_loc(i,j)-lat_loc(i,j+1))
2975	ENDDO
2976	dy = dy/(jjm_loc-1)
2977	DO j=1,jjm_loc
2978	corners_loc(i,j,1,2) = lat_loc(i,j)+dy/2.0
2979	corners_loc(i,j,2,2) = lat_loc(i,j)+dy/2.0
2980	corners_loc(i,j,3,2) = lat_loc(i,j)-dy/2.0
2981	corners_loc(i,j,4,2) = lat_loc(i,j)-dy/2.0
2982	ENDDO
2983	ENDDO
2984	!
2985	! Compute the areas of the grid (using the simplification that the grid is regular in lon/lat).
2986	!
2987	DO i=1,iim_loc
2988	DO j=1,jjm_loc
2989	coslat = MAX( COS(lat_loc(i,j) * pi/180. ), mincos )
2990	dx = ABS(corners_loc(i,j,2,1) - corners_loc(i,j,1,1)) * pi/180. * R_Earth * coslat
2991	dy = ABS(corners_loc(i,j,1,2) - corners_loc(i,j,3,2)) * pi/180. * R_Earth
2992	area_loc(i,j) = dx*dy
2993	ENDDO
2994	ENDDO
2995	!
2996	! If requested we read a variable, zoomin and dump the result into a test file.
2997	!
2998	IF ( dumpncdf ) THEN
2999	iret = NF90_OPEN (filename, NF90_NOWRITE, force_id)
3000	IF (iret /= NF90_NOERR) THEN
3001	CALL ipslerr (3,'forcing_zoomgrid',"Error opening the forcing file :", filename, " ")
3002	ENDIF
3003	!
3004	ALLOCATE(readvar(ncdfcount), readvar2d(iim_glo,jjm_glo), zoomedvar(nbland_loc))
3005	!
3006	iret = NF90_INQ_VARID(force_id, TRIM(testvarname), iv)
3007	IF (iret /= NF90_NOERR) THEN
3008	CALL ipslerr (3,'forcing_zoomgrid',"Could not find variable Tair in file."," "," ")
3009	ENDIF
3010
3011	IF ( compressed ) THEN
3012	!
3013	start(1) = ncdfstart
3014	start(2) = 1
3015	count(1) = ncdfcount
3016	count(2) = 1
3017	!
3018	iret = NF90_GET_VAR(force_id, iv, readvar, start, count)
3019	IF (iret /= NF90_NOERR) THEN
3020	CALL ipslerr (3,'forcing_zoomgrid',"Could not read compressed variable Tair from file."," "," ")
3021	ENDIF
3022	CALL forcing_reindex(ncdfcount, readvar, nbland_loc, zoomedvar, reindex_loc)
3023	!
3024	ELSE
3025	!
3026	start2d(1) = 1
3027	start2d(2) = 1
3028	start2d(3) = 1
3029	count2d(1) = iim_glo
3030	count2d(2) = jjm_glo
3031	count2d(3) = 1
3032	!
3033	iret = NF90_GET_VAR(force_id, iv, readvar2d, start2d, count2d)
3034	IF (iret /= NF90_NOERR) THEN
3035	CALL ipslerr (3,'forcing_zoomgrid',"Could not read variable Tair from file."," "," ")
3036	ENDIF
3037	CALL forcing_reindex(iim_glo, jjm_glo, readvar2d, nbland_loc, zoomedvar, reindex2d_loc)
3038	!
3039	ENDIF
3040	!
3041	CALL forcing_writetestvar("forcing_mask_loc.nc", iim_loc, jjm_loc, nbland_loc, &
3042	& lon_loc(:,1), lat_loc(1,:), lindex_loc, mask_loc, &
3043	& TRIM(testvarname), zoomedvar)
3044	!
3045	ENDIF
3046	!
3047	! Clean up
3048	!
3049	IF ( ALLOCATED(readvar) ) DEALLOCATE(readvar)
3050	IF ( ALLOCATED(readvar2d) ) DEALLOCATE(readvar2d)
3051	IF ( ALLOCATED(zoomedvar) ) DEALLOCATE(zoomedvar)
3052	IF ( ALLOCATED(index_glotoloc) ) DEALLOCATE(index_glotoloc)
3053	IF ( ALLOCATED(lalo) ) DEALLOCATE(lalo)
3054	!
3055	END SUBROUTINE forcing_zoomgrid
3056
3057	!! =============================================================================================================================
3058	!! SUBROUTINE: forcing_givegridsize
3059	!!
3060	!>\BRIEF Routine which exports the size of the grid on which the model will run, i.e. the zoomed grid.
3061	!!
3062	!! DESCRIPTION: This is needed to transfer the grid information from this module to the glogrid.f90 module.
3063	!!
3064	!! \n
3065	!_ ==============================================================================================================================
3066
3067	SUBROUTINE forcing_givegridsize (iim, jjm, nblindex)
3068	!
3069	! Provides the size of the grid to be used to the calling program
3070	!
3071	! Size of the x and y direction of the zoomed area
3072	INTEGER(i_std), INTENT(out) :: iim, jjm
3073	! Number of land points in the zoomed area
3074	INTEGER(i_std), INTENT(out) :: nblindex
3075	!
3076	IF ( .NOT. is_root_prc ) THEN
3077	CALL ipslerr (3,'forcing_givegridsize'," This routine can only be called on the root processor.", &
3078	& "The information requested is only available on root processor.", " ")
3079	ENDIF
3080	!
3081	iim = iim_loc
3082	jjm = jjm_loc
3083	nblindex = nbland_loc
3084	!
3085	END SUBROUTINE forcing_givegridsize
3086
3087	!! =============================================================================================================================
3088	!! SUBROUTINE: forcing_
3089	!!
3090	!>\BRIEF
3091	!!
3092	!! DESCRIPTION:
3093	!!
3094	!! \n
3095	!_ ==============================================================================================================================
3096
3097	SUBROUTINE forcing_vertical(force_id)
3098	!
3099	!- This subroutine explores the forcing file to decide what information is available
3100	!- on the vertical coordinate.
3101	!
3102	INTEGER, INTENT(IN) :: force_id
3103	!
3104	INTEGER(i_std) :: iret, ireta, iretb
3105	!
3106	INTEGER(i_std) :: sigma_id = -1, sigma_uv_id = -1
3107	INTEGER(i_std) :: hybsiga_id = -1, hybsiga_uv_id = -1
3108	INTEGER(i_std) :: hybsigb_id = -1, hybsigb_uv_id = -1
3109	INTEGER(i_std) :: levels_id = -1, levels_uv_id = -1
3110	INTEGER(i_std) :: height_id = -1, height_uv_id = -1
3111	INTEGER(i_std) :: lev_id = -1
3112	!
3113	LOGICAL :: var_exists, vara_exists, varb_exists, varuv_exists
3114	LOGICAL :: foundvar
3115	LOGICAL :: levlegacy
3116	!
3117	!- Set all the defaults
3118	!
3119	zfixed=.FALSE.
3120	zsigma=.FALSE.
3121	zhybrid=.FALSE.
3122	zlevels=.FALSE.
3123	zheight=.FALSE.
3124	zsamelev_uv = .TRUE.
3125	levlegacy = .FALSE.
3126	!
3127	foundvar = .FALSE.
3128	!
3129	!- We have a forcing file to explore so let us see if we find any of the conventions
3130	!- which allow us to find the height of T,Q,U and V.
3131	!
3132	IF ( force_id > 0 ) THEN
3133	!
3134	! Case for sigma levels
3135	!
3136	IF ( .NOT. foundvar ) THEN
3137	ireta = NF90_INQ_VARID(force_id, 'Sigma', sigma_id)
3138	IF ( (sigma_id >= 0) .AND. (ireta == NF90_NOERR) ) THEN
3139	foundvar = .TRUE.
3140	zsigma = .TRUE.
3141	iretb = NF90_INQ_VARID(force_id, 'Sigma_uv', sigma_uv_id)
3142	IF ( (sigma_uv_id >= 0) .OR. (iretb == NF90_NOERR) ) zsamelev_uv = .FALSE.
3143	ENDIF
3144	ENDIF
3145	!
3146	! Case for Hybrid levels
3147	!
3148	IF ( .NOT. foundvar ) THEN
3149	var_exists = .FALSE.
3150	ireta = NF90_INQ_VARID(force_id, 'HybSigA', hybsiga_id)
3151	IF ( (hybsiga_id >= 0 ) .AND. (ireta == NF90_NOERR) ) THEN
3152	iretb = NF90_INQ_VARID(force_id, 'HybSigB', hybsigb_id)
3153	IF ( (hybsigb_id >= 0 ) .AND. (iretb == NF90_NOERR) ) THEN
3154	var_exists=.TRUE.
3155	ELSE
3156	CALL ipslerr ( 3, 'forcing_vertical','Missing the B coefficient for', &
3157	& 'Hybrid vertical levels for T and Q','stop')
3158	ENDIF
3159	ENDIF
3160	ireta = NF90_INQ_VARID(force_id, 'HybSigA_uv', hybsiga_uv_id)
3161	IF ( (hybsiga_uv_id >= 0 ) .AND. (ireta == NF90_NOERR) ) THEN
3162	iretb = NF90_INQ_VARID(force_id, 'HybSigB_uv', hybsigb_uv_id)
3163	IF ( (hybsigb_uv_id >= 0 ) .AND. (iretb == NF90_NOERR) ) THEN
3164	varuv_exists=.TRUE.
3165	ELSE
3166	CALL ipslerr ( 3, 'forcing_vertical','Missing the B coefficient for', &
3167	& 'Hybrid vertical levels for U and V','stop')
3168	ENDIF
3169	ENDIF
3170	IF ( var_exists ) THEN
3171	foundvar = .TRUE.
3172	zhybrid = .TRUE.
3173	IF ( varuv_exists ) zsamelev_uv = .FALSE.
3174	ENDIF
3175	ENDIF
3176	!
3177	! Case for levels (i.e. a 2d time varying field with the height in meters)
3178	!
3179	IF ( .NOT. foundvar ) THEN
3180	ireta = NF90_INQ_VARID(force_id, 'Levels', levels_id)
3181	IF ( (levels_id >= 0 ) .AND. (ireta == NF90_NOERR) ) THEN
3182	foundvar = .TRUE.
3183	zlevels = .TRUE.
3184	iretb = NF90_INQ_VARID(force_id, 'Levels_uv', levels_uv_id)
3185	IF ( (levels_uv_id >= 0 ) .AND. (iretb == NF90_NOERR) ) zsamelev_uv = .FALSE.
3186	ENDIF
3187	ENDIF
3188	!
3189	! Case where a fixed height is provided in meters
3190	!
3191	IF ( .NOT. foundvar ) THEN
3192	ireta = NF90_INQ_VARID(force_id, 'Height_Lev1', height_id)
3193	IF ( (height_id >= 0 ) .AND. (ireta == NF90_NOERR) ) THEN
3194	foundvar = .TRUE.
3195	zheight = .TRUE.
3196	iretb = NF90_INQ_VARID(force_id, 'Height_Levuv', height_uv_id)
3197	IF ( (height_uv_id >= 0 ) .AND. (iretb == NF90_NOERR) ) zsamelev_uv = .FALSE.
3198	ENDIF
3199	ENDIF
3200	!
3201	! Case where a fixed height is provided in meters in the lev variable
3202	!
3203	IF ( .NOT. foundvar ) THEN
3204	ireta = NF90_INQ_VARID(force_id, 'lev', lev_id)
3205	IF ( (lev_id >= 0 ) .AND. (ireta == NF90_NOERR) ) THEN
3206	foundvar = .TRUE.
3207	zheight = .TRUE.
3208	levlegacy = .TRUE.
3209	ENDIF
3210	ENDIF
3211	!
3212	ENDIF
3213	!
3214	! We found forcing variables so we need to extract the values if we are dealing with constant values (i.e. all
3215	! except the case zlevels
3216	!
3217	IF ( foundvar .AND. .NOT. zlevels ) THEN
3218	!
3219	IF ( zheight ) THEN
3220	!
3221	! Constant height
3222	!
3223	IF ( levlegacy ) THEN
3224	iret = NF90_GET_VAR(force_id, lev_id, zlev_fixed)
3225	IF ( iret /= NF90_NOERR ) THEN
3226	CALL ipslerr ( 3, 'forcing_vertical','Attempted to read variable lev from forcing file in legacy mode', &
3227	& 'NF90_GET_VAR failed.','stop')
3228	ENDIF
3229	ELSE
3230	iret = NF90_GET_VAR(force_id, height_id, zlev_fixed)
3231	IF ( iret /= NF90_NOERR ) THEN
3232	CALL ipslerr ( 3, 'forcing_vertical','Attempted to read variable Height_Lev1 from forcing file', &
3233	& 'NF90_GET_VAR failed.','stop')
3234	ENDIF
3235	IF ( .NOT. zsamelev_uv ) THEN
3236	iret = NF90_GET_VAR(force_id, height_uv_id, zlevuv_fixed)
3237	IF ( iret /= NF90_NOERR ) THEN
3238	CALL ipslerr ( 3, 'forcing_vertical','Attempted to read variable Height_Levuv from forcing file', &
3239	& 'NF90_GET_VAR failed.','stop')
3240	ENDIF
3241	ENDIF
3242	ENDIF
3243	WRITE(,) "forcing_vertical : case ZLEV : Read from forcing file :", zlev_fixed, zlevuv_fixed
3244	!
3245	ELSE IF ( zsigma .OR. zhybrid ) THEN
3246	!
3247	! Sigma or hybrid levels
3248	!
3249	IF ( zsigma ) THEN
3250	iret = NF90_GET_VAR(force_id, sigma_id, zhybrid_b)
3251	zhybrid_a = zero
3252	IF ( .NOT. zsamelev_uv ) THEN
3253	iret = NF90_GET_VAR(force_id, sigma_uv_id, zhybriduv_b)
3254	zhybriduv_a = zero
3255	ENDIF
3256	ELSE
3257	ireta = NF90_GET_VAR(force_id, hybsigb_id, zhybrid_b)
3258	iretb = NF90_GET_VAR(force_id, hybsiga_id, zhybrid_a)
3259	IF ( ireta /= NF90_NOERR .OR. iretb /= NF90_NOERR) THEN
3260	CALL ipslerr ( 3, 'forcing_vertical','Attempted to read variable HybSigA and HybSigB from forcing file', &
3261	& 'NF90_GET_VAR failed.','stop')
3262	ENDIF
3263	IF ( .NOT. zsamelev_uv ) THEN
3264	ireta = NF90_GET_VAR(force_id, hybsigb_uv_id, zhybriduv_b)
3265	iretb = NF90_GET_VAR(force_id, hybsiga_uv_id, zhybriduv_a)
3266	IF ( ireta /= NF90_NOERR .OR. iretb /= NF90_NOERR) THEN
3267	CALL ipslerr ( 3, 'forcing_vertical',&
3268	& 'Attempted to read variable HybSigA_uv and HybSigB_uv from forcing file', &
3269	& 'NF90_GET_VAR failed.','stop')
3270	ENDIF
3271	ENDIF
3272	ENDIF
3273	WRITE(,) "forcing_vertical : case Pressure coordinates : "
3274	WRITE(,) "Read from forcing file :", zhybrid_b, zhybrid_a, zhybriduv_b, zhybriduv_a
3275	ELSE
3276	!
3277	! Why are we here ???
3278	!
3279	CALL ipslerr ( 3, 'forcing_vertical','What is the option used to describe the height of', &
3280	& 'the atmospheric forcing ?','Please check your forcing file.')
3281	ENDIF
3282	ENDIF
3283	!
3284	!- We have no forcing file to explore or we did not find anything. So revert back to the run.def and
3285	!- read what has been specified by the user.
3286	!
3287	IF ( force_id < 0 .OR. .NOT. foundvar ) THEN
3288	!
3289	!-
3290	!Config Key = HEIGHT_LEV1
3291	!Config Desc = Height at which T and Q are given
3292	!Config Def = 2.0
3293	!Config Help = The atmospheric variables (temperature and specific
3294	!Config humidity) are measured at a specific level.
3295	!Config The height of this level is needed to compute
3296	!Config correctly the turbulent transfer coefficients.
3297	!Config Look at the description of the forcing
3298	!Config DATA for the correct value.
3299	!-
3300	zlev_fixed = 2.0
3301	CALL getin('HEIGHT_LEV1', zlev_fixed)
3302	!-
3303	!Config Key = HEIGHT_LEVW
3304	!Config Desc = Height at which the wind is given
3305	!Config Def = 10.0
3306	!Config Help = The height at which wind is needed to compute
3307	!Config correctly the turbulent transfer coefficients.
3308	!-
3309	zlevuv_fixed = 10.0
3310	CALL getin('HEIGHT_LEVW', zlevuv_fixed)
3311
3312	zheight = .TRUE.
3313
3314	IF ( ABS(zlevuv_fixed-zlev_fixed) > EPSILON(zlev_fixed)) THEN
3315	zsamelev_uv = .FALSE.
3316	ELSE
3317	zsamelev_uv = .TRUE.
3318	ENDIF
3319
3320	CALL ipslerr ( 2, 'forcing_vertical','The height of the atmospheric forcing variables', &
3321	& 'was not found in the netCDF file.','Thus the values in run.def were used ... or their defaults.')
3322	ENDIF
3323
3324	END SUBROUTINE forcing_vertical
3325
3326	!! =============================================================================================================================
3327	!! SUBROUTINE: forcing_givegrid
3328	!!
3329	!>\BRIEF Routine which exports the grid (longitude, latitude, land indices) on which the model will run, i.e. the zoomed grid.
3330	!!
3331	!! DESCRIPTION: This is needed to transfer the grid information from this module to the glogrid.f90 module.
3332	!!
3333	!!
3334	!! \n
3335	!_ ==============================================================================================================================
3336
3337	SUBROUTINE forcing_givegrid (lon, lat, mask, area, corners, lindex, contfrac, calendar_tmp)
3338	!
3339	! This subroutine will return to the caller the grid which has been extracted from the
3340	! the forcing file. It is assumed that the caller has called forcing_givegridsize before
3341	! and knows the dimensions of the fields and thus has done the correct allocations.
3342	!
3343	!
3344	REAL(r_std), INTENT(out) :: lon(iim_loc,jjm_loc), lat(iim_loc,jjm_loc)
3345	REAL(r_std), INTENT(out) :: mask(iim_loc,jjm_loc)
3346	REAL(r_std), INTENT(out) :: area(iim_loc,jjm_loc)
3347	REAL(r_std), INTENT(out) :: corners(iim_loc,jjm_loc,4,2)
3348	INTEGER(i_std), INTENT(out) :: lindex(nbland_loc)
3349	REAL(r_std), INTENT(out) :: contfrac(nbland_loc)
3350	CHARACTER(LEN=20), INTENT(out) :: calendar_tmp
3351	!
3352	IF ( .NOT. is_root_prc ) THEN
3353	CALL ipslerr (3,'forcing_givegrid'," This routine can only be called on the root processor.", &
3354	& "The information requested is only available on root processor.", " ")
3355	ENDIF
3356	!
3357	lon(:,:) = lon_loc(:,:)
3358	lat(:,:) = lat_loc(:,:)
3359	!
3360	mask(:,:) = mask_loc(:,:)
3361	area(:,:) = area_loc(:,:)
3362	corners(:,:,:,:) = corners_loc(:,:,:,:)
3363	!
3364	!
3365	lindex(:) = lindex_loc(:)
3366	contfrac(:) = contfrac_loc(:)
3367	!
3368	calendar_tmp = calendar
3369	!
3370	END SUBROUTINE forcing_givegrid
3371
3372	!! =============================================================================================================================
3373	!! SUBROUTINE: forcing_checkdim
3374	!!
3375	!>\BRIEF
3376	!!
3377	!! DESCRIPTION: Save the dimension or check that it is equal to the previous value.
3378	!! Should one of the spatial dimensions be different between 2 files, then we have a big problem.
3379	!! They probably do not belong to the same set of forcing files.
3380	!!
3381	!! \n
3382	!_ ==============================================================================================================================
3383
3384	SUBROUTINE forcing_checkdim(ifile, filenames, out_dim, out_id, in_dim, in_id)
3385	!
3386	! Save the dimension or check that it is equal to the previous value.
3387	! Should one of the spatial dimensions be different between 2 files, then we have a big problem.
3388	! They probably do not belong to the same set of forcing files.
3389	!
3390	INTEGER(i_std), INTENT(in) :: ifile
3391	CHARACTER(LEN=*), INTENT(in) :: filenames(:)
3392	INTEGER(i_std), INTENT(out) :: out_dim, out_id
3393	INTEGER(i_std), INTENT(in) :: in_dim, in_id
3394	!
3395	IF ( ifile == 1 ) THEN
3396	out_dim = in_dim
3397	out_id = in_id
3398	ELSE
3399	IF ( out_dim /= in_dim ) THEN
3400	CALL ipslerr (3,'forcing_ocheckdim', 'The dimension of the file is not the same of the first file opened.', &
3401	& 'The offending file is : ', filenames(ifile))
3402	ENDIF
3403	ENDIF
3404	!
3405	END SUBROUTINE forcing_checkdim
3406
3407	!! =============================================================================================================================
3408	!! SUBROUTINE: forcing_time
3409	!!
3410	!>\BRIEF Read the time from each file and create the time axis to be the basis for the simulation.
3411	!!
3412	!! DESCRIPTION: This is an important routine which analyses the time axis of the forcing files and
3413	!! stores the main information in the SAVED variables of this routine.
3414	!! As this module manages a list of forcing files we also need to check that the time
3415	!! axis of all these files is continuous and homogeneous.
3416	!! The bounds are also build for all the time axes so that we know how to interpret the
3417	!! various variables.
3418	!!
3419	!! \n
3420	!_ ==============================================================================================================================
3421
3422	SUBROUTINE forcing_time(nbfiles, filenames)
3423	!
3424	! Read the time from each file and create the time axis to be the basis
3425	! for the simulation.
3426	!
3427	INTEGER(i_std) :: nbfiles
3428	CHARACTER(LEN=*) :: filenames(nbfiles)
3429	!
3430	INTEGER(i_std) :: iv, it, iff, tcnt, itbase, itbasetmp, ittmp
3431	INTEGER(i_std) :: tstart, maxtime_infile
3432	REAL(r_std), ALLOCATABLE, DIMENSION(:) :: timeint, time_read
3433	REAL(r_std), ALLOCATABLE, DIMENSION(:) :: time_infiles
3434	CHARACTER(LEN=20) :: axname, calendar, timevarname
3435	CHARACTER(LEN=60) :: timestamp, tmpatt
3436	INTEGER(i_std) :: tncstart(3), tnccount(3)
3437	!
3438	INTEGER(i_std) :: iret, year0, month0, day0, hours0, minutes0, seci
3439	INTEGER(i_std), DIMENSION(1) :: imax, imin
3440	REAL(r_std) :: sec0, date_int, date0_tmp
3441	CHARACTER :: strc
3442	LOGICAL :: check=.FALSE.
3443	!
3444	! Check that we are working on the root proc.
3445	!
3446	IF ( .NOT. is_root_prc) THEN
3447	CALL ipslerr (3,'forcing_time',"Cannot run this routine o other procs than root.",&
3448	& "All the information on the forcing files is only lated on the root processor."," ")
3449	ENDIF
3450	!
3451	! Size of unlimited dimension added up through the files. If variable not allocated before by another
3452	! subroutine, it needs to be done here.
3453	!
3454	IF ( .NOT. ALLOCATED(nbtime_perfile) ) ALLOCATE(nbtime_perfile(nbfiles))
3455	IF ( .NOT. ALLOCATED(date0_file) ) ALLOCATE(date0_file(nbfiles,nbtax))
3456	!
3457	! Go through all files in the list in order to get the total number of time steps we have
3458	! in the nbfiles files to be read
3459	!
3460	nb_forcing_steps = 0
3461	maxtime_infile = 0
3462	DO iff=1,nbfiles
3463	!
3464	iret = NF90_INQUIRE_DIMENSION(force_id(iff), id_unlim(iff), name=axname, len=nbtime_perfile(iff))
3465	IF (iret /= NF90_NOERR) THEN
3466	CALL ipslerr (3,'forcing_time',"Could not get size of dimension of unlimited axis"," "," ")
3467	ENDIF
3468	nb_forcing_steps = nb_forcing_steps + nbtime_perfile(iff)
3469	IF ( nbtime_perfile(iff) > maxtime_infile ) maxtime_infile = nbtime_perfile(iff)
3470	ENDDO
3471	!
3472	! Allocate the variables needed with the time length just calculated.
3473	! These variables are saved in the module
3474	!
3475	ALLOCATE(time_infiles(nb_forcing_steps))
3476	ALLOCATE(time(nb_forcing_steps, nbtaxnbtmethods), time_bounds(nb_forcing_steps,nbtaxnbtmethods,2))
3477	ALLOCATE(time_axename(nbtaxnbtmethods), time_cellmethod(nbtaxnbtmethods))
3478	ALLOCATE(preciptime(nb_forcing_steps))
3479	ALLOCATE(time_sourcefile(nb_forcing_steps))
3480	ALLOCATE(time_id(nb_forcing_steps, nbtax))
3481	! Allocate local variables
3482	ALLOCATE(time_read(nb_forcing_steps))
3483	ALLOCATE(timeint(nb_forcing_steps))
3484	!
3485	! Get through all variables to find time_id
3486	! The key variables to filled up here are time (the time stamps read in the file) and
3487	! time_bounds which give the validity interval for the variables.
3488	!
3489	tstart=0
3490	!
3491	IF ( check ) WRITE(,) "forcing_time : going through ", nbfiles, " files to get the time."
3492	!
3493	DO iff=1,nbfiles
3494	!
3495	time_id(iff,:)=-1
3496	!
3497	! Go through the variables in the file and find the one which is a time axis.
3498	!
3499	tcnt=1
3500	DO iv=1,nvars(iff)
3501	iret = NF90_GET_ATT(force_id(iff), iv, "units", tmpatt)
3502	IF ( INDEX(lowercase(tmpatt),'seconds since') > 0) THEN
3503	time_id(iff,tcnt)=iv
3504	tcnt=tcnt+1
3505	convtosec(iff)=1.0
3506	ELSE IF ( INDEX(lowercase(tmpatt),'minutes since') > 0) THEN
3507	time_id(iff,tcnt)=iv
3508	tcnt=tcnt+1
3509	convtosec(iff)=60.0
3510	ELSE IF ( INDEX(lowercase(tmpatt),'hours since') > 0) THEN
3511	time_id(iff,tcnt)=iv
3512	tcnt=tcnt+1
3513	convtosec(iff)=3600.0
3514	ENDIF
3515	ENDDO
3516	IF ( ANY(time_id(iff,:) < 0) ) THEN
3517	CALL ipslerr (3,'forcing_time',"Incorrect numer of time axes. A time axis is missing ",&
3518	& "in file :", filenames(iff))
3519	ENDIF
3520	!
3521	IF ( check ) WRITE(,) "forcing_time : Looking at time axis for file ", force_id(iff)
3522	!
3523	! Looping through the time axes and read them.
3524	!
3525	DO tcnt=1,nbtax
3526	!
3527	iret = NF90_INQUIRE_VARIABLE(force_id(iff), time_id(iff,tcnt), name=timevarname)
3528	IF ( check ) WRITE(,) "forcing_time : in ", iff, " found variable ", timevarname
3529	!
3530	! Get units of time axis
3531	!
3532	iret = NF90_GET_ATT(force_id(iff), time_id(iff,tcnt), "units", timestamp)
3533	IF ( check ) WRITE(,) "forcing_time : has time stamp ", timestamp
3534	!
3535	! Transform the start date of the netCDF file into a julian date for the model
3536	!
3537	timestamp = TRIM(timestamp(INDEX(timestamp,'since')+6:LEN_TRIM(timestamp)))
3538	!
3539	! Temporary fix. We need a more general method to find the right format for reading
3540	! the elements of the start date.
3541	IF ( LEN_TRIM(timestamp) == 14 ) THEN
3542	READ (timestamp,'(I4.4,5(a,I1))') &
3543	year0, strc, month0, strc, day0, &
3544	strc, hours0, strc, minutes0, strc, seci
3545	ELSE
3546	READ (timestamp,'(I4.4,5(a,I2.2))') &
3547	year0, strc, month0, strc, day0, &
3548	strc, hours0, strc, minutes0, strc, seci
3549	ENDIF
3550	sec0 = hours03600. + minutes060. + seci
3551	CALL ymds2ju (year0, month0, day0, sec0, date0_tmp)
3552	date0_file(iff,tcnt) = date0_tmp
3553	!
3554	! Now get the actual dates
3555	!
3556	tncstart(1) = 1
3557	tnccount(1) = nbtime_perfile(iff)
3558	IF ( check ) WRITE(,) "forcing_time : number of values read : ", tnccount(1)
3559	iret = NF90_GET_VAR(force_id(iff), time_id(iff,tcnt), time_read, tncstart, tnccount)
3560	IF (iret /= NF90_NOERR) THEN
3561	CALL ipslerr (3,'forcing_time',"An error occured while reading time from the file."," "," ")
3562	ENDIF
3563	!
3564	! Convert the variable time from seconds since to julian days
3565	!
3566	DO it=1,nbtime_perfile(iff)
3567	time_infiles(tstart+it) = date0_file(iff,tcnt) + time_read(it)*convtosec(iff)/one_day
3568	ENDDO
3569	if ( check ) WRITE(,) "File ", iff, "goes from ", time_infiles(tstart+1), " to ", &
3570	time_infiles(tstart+nbtime_perfile(iff))
3571	!
3572	! Estimate the bounds as this information is not yet in the forcing file.
3573	!
3574	date_int = (time_infiles(tstart+nbtime_perfile(iff)) - time_infiles(tstart+1))/(nbtime_perfile(iff)-1)
3575	forcing_tstep_ave = date_int*one_day
3576	!
3577	! If this is the first file then simply keep the name of the time axis. Else search the same name
3578	! in what has already been read
3579	!
3580	IF ( iff == 1 ) THEN
3581	itbase = (tcnt-1)*nbtmethods
3582	time_axename(itbase+1:itbase+4) = timevarname
3583	time_cellmethod(itbase+1) = "reference"
3584	time_cellmethod(itbase+2) = "start"
3585	time_cellmethod(itbase+3) = "cent"
3586	time_cellmethod(itbase+4) = "end"
3587	ELSE
3588	!
3589	! If this is not the first file then we need to find the correct axis to add to.
3590	! All information have already been saved with the first file.
3591	!
3592	DO ittmp=1,nbtax
3593	itbasetmp=(ittmp-1)*nbtmethods
3594	IF ( time_axename(itbasetmp+1) == timevarname ) THEN
3595	itbase = itbasetmp
3596	ENDIF
3597	ENDDO
3598
3599	ENDIF
3600	!
3601	!
3602	! Keep for future usage the various information on the time axis we have just read. This time axis can
3603	! be understood in 3 different ways and each variable might use a different cell method for this time
3604	! axis.
3605	!
3606	! time(:,(tcnt-1)*nbtmethods+1) : corresponds to the reference time axis as it has been read from the file
3607	! time(:,(tcnt-1)*nbtmethods+2) : is the time axis with a cell method which place the value at the
3608	! beginning of the time interval
3609	! time(:,(tcnt-1)*nbtmethods+3) : is the time axis corresponding to variables placed at the center of the
3610	! time interval
3611	! time(:,(tcnt-1)*nbtmethods+4) : for variables put at the end of the time interval over which they aere
3612	! for instance averaged.
3613	!
3614	! In variable time_cellmethod we will write the type of cell method as descirbed above so that the selection
3615	! of the right axis for each variable can be made automaticaly.
3616	!
3617	! We also keep the name of the time axis read in preparation of file where we might have to read more than one
3618	! time axis.
3619	!
3620	DO it=tstart+1,tstart+nbtime_perfile(iff)
3621	!
3622	! Reference time
3623	!
3624	time(it,itbase+1) = time_infiles(it)
3625	time_bounds(it,itbase+1,1) = time_infiles(it)-date_int/2.0
3626	time_bounds(it,itbase+1,2) = time_infiles(it)+date_int/2.0
3627	!
3628	! Start cell method
3629	time(it,itbase+2) = time_infiles(it)+date_int/2.0
3630	time_bounds(it,itbase+2,1) = time_infiles(it)
3631	time_bounds(it,itbase+2,2) = time_infiles(it)+date_int
3632	!
3633	! Centered cell method
3634	time(it,itbase+3) = time_infiles(it)
3635	time_bounds(it,itbase+3,1) = time_infiles(it)+date_int/2.0
3636	time_bounds(it,itbase+3,2) = time_infiles(it)+date_int/2.0
3637	!
3638	! End cell method
3639	time(it,itbase+4) = time_infiles(it)-date_int/2.0
3640	time_bounds(it,itbase+4,1) = time_infiles(it)-date_int
3641	time_bounds(it,itbase+4,2) = time_infiles(it)
3642	!
3643	ENDDO
3644	!
3645	! Keep the number of the file from which we read this time.
3646	!
3647	time_sourcefile(tstart+1:tstart+nbtime_perfile(iff))=iff
3648	!
3649	IF ( check ) WRITE(,) "forcing_time : finished file ", iff
3650	!
3651	ENDDO
3652	!
3653	! Before moving to the next file advance the pointer in the time arrays.
3654	!
3655	tstart=tstart+nbtime_perfile(iff)
3656	!
3657	ENDDO
3658	!
3659	IF ( check ) WRITE(,) "forcing_time : All files have been processed"
3660	!
3661	! Verify that the forcing comes in regular time intervals. If not, many of the
3662	! interpolation schemes will fail.
3663	! This is only done on the first time axis ... is it enough ?
3664	!
3665	DO ittmp=1,nbtax
3666	itbase=(ittmp-1)*nbtmethods
3667	!
3668	date_int = (time(nb_forcing_steps,itbase+1) - time(1,itbase+1))/(nb_forcing_steps-1)
3669	forcing_tstep_ave = date_int*one_day
3670	!
3671	timeint(:) = 0
3672	DO it=1, nb_forcing_steps-1
3673	timeint(it) = time(it+1,itbase+1)-time(it,itbase+1)
3674	ENDDO
3675	!
3676	IF ( MAXVAL(timeint(1:nb_forcing_steps-1)) > date_int+0.1*date_int .OR.&
3677	& MINVAL(timeint(1:nb_forcing_steps-1)) < date_int-0.1*date_int) THEN
3678	WRITE(,) "The time steping of the forcing files does not seem to be regular on axis",time_axename(itbase+1),":"
3679	WRITE(,) "Average time step : ", date_int, "days = ", date_int*one_day, "sec."
3680	imax = MAXLOC(timeint(1:nb_forcing_steps-1))
3681	imin = MINLOC(timeint(1:nb_forcing_steps-1))
3682	WRITE(,) "Maximum time step : ", MAXVAL(timeint(1:nb_forcing_steps-1)), " at ", imax(1)
3683	WRITE(,) "Minimum time step : ", MINVAL(timeint(1:nb_forcing_steps-1)), " at ", imin(1)
3684	WRITE(,) "++++ Values around Maximum"
3685	DO it=MAX(imax(1)-5,1),MIN(imax(1)+5,nb_forcing_steps)
3686	WRITE(,) it, " from file ", time_sourcefile(it), " Value ", time(it,itbase+1)
3687	CALL forcing_printdate(time(it,itbase+1), "Time values.")
3688	ENDDO
3689	WRITE(,) "++++ Values around Minimum"
3690	DO it=MAX(imin(1)-5,1),MIN(imin(1)+5,nb_forcing_steps)
3691	WRITE(,) it, " from file ", time_sourcefile(it), " Value ", time(it,itbase+1)
3692	CALL forcing_printdate(time(it,itbase+1), "Time values.")
3693	ENDDO
3694	CALL ipslerr (3,'forcing_time', 'The time handling could be improved to allow the driver',&
3695	& "to cope with irregular forcing files."," ")
3696	ENDIF
3697	ENDDO
3698	!
3699	! Print some test values
3700	!
3701	DO ittmp=1,nbtax
3702	itbase=(ittmp-1)*nbtmethods
3703	!
3704	WRITE(,) "Bounds for axis ",time_axename(itbase+1)," :"
3705	!
3706	CALL forcing_printdate(time_bounds(1,itbase+1,1), "Start time of first forcing interval.")
3707	CALL forcing_printdate(time_bounds(1,itbase+1,2), "End time of first forcing interval.")
3708	CALL forcing_printdate(time_bounds(nb_forcing_steps,itbase+1,1), "Start time of last forcing interval.")
3709	CALL forcing_printdate(time_bounds(nb_forcing_steps,itbase+1,2), "End time of last forcing interval.")
3710	ENDDO
3711	!
3712	! Set to zero the variable for the cummulated time for rainfall
3713	!
3714	preciptime(:) = zero
3715	!
3716	! Keep the very first date of the time axis for future use
3717	!
3718	forcingstartdate = time(1,1)
3719	!
3720	! Clean-up
3721	!
3722	DEALLOCATE(timeint, time_read)
3723	!
3724	END SUBROUTINE forcing_time
3725
3726	!! =============================================================================================================================
3727	!! SUBROUTINE: forcing_varforslab
3728	!!
3729	!>\BRIEF
3730	!!
3731	!! DESCRIPTION: This subroutine will read the named variable and put it in the right place in the
3732	!! slab of data kept in the memory of the driver.
3733	!!
3734	!! \n
3735	!_ ==============================================================================================================================
3736
3737	SUBROUTINE forcing_varforslab(fileindex, varname, timestart, timecount, inslabpos, data, cellmethod)
3738	!
3739	! This subroutine will read the named variable and put it in the right place in the
3740	! slab of data kept in the memory of the driver.
3741	!
3742	INTEGER(i_std), INTENT(in) :: fileindex
3743	CHARACTER(LEN=*), INTENT(in) :: varname
3744	INTEGER(i_std), INTENT(in) :: timestart, timecount, inslabpos
3745	REAL(r_std), INTENT(inout) :: data(nbland_loc,slab_size)
3746	CHARACTER(LEN=*), INTENT(out) :: cellmethod
3747	!
3748	! Local variables
3749	!
3750	INTEGER(i_std) :: varid, windid, windndims, it, ig, iv
3751	INTEGER(i_std) :: iret, ndims
3752	INTEGER(i_std), DIMENSION(4) :: start, count
3753	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: tmp_slab
3754	REAL(r_std), ALLOCATABLE, DIMENSION(:,:,:) :: tmp_slab2d
3755	CHARACTER(LEN=80) :: name
3756	LOGICAL :: windzero
3757	!
3758	! Allocate the temporary data array if not already available
3759	!
3760	IF ( compressed ) THEN
3761	IF ( .NOT. ALLOCATED(tmp_slab) ) ALLOCATE(tmp_slab(ncdfcount,slab_size))
3762	ELSE
3763	IF ( .NOT. ALLOCATED(tmp_slab2d) ) ALLOCATE(tmp_slab2d(iim_glo,jjm_glo,slab_size))
3764	ENDIF
3765	!
3766	! Reset the counters and flags to forget the past !
3767	!
3768	varid=-1
3769	windid=-1
3770	windzero=.FALSE.
3771	!
3772	! Find the variable in the file
3773	!
3774	DO iv=1,nvars(fileindex)
3775	!
3776	iret = NF90_INQUIRE_VARIABLE(force_id(fileindex), iv, name=name, ndims=it)
3777	!
3778	IF ( INDEX(name, varname) > 0 ) THEN
3779	varid = iv
3780	ndims = it
3781	ENDIF
3782	IF ( (INDEX(name, "Wind") > 0) .AND. (LEN_TRIM(name) == LEN_TRIM("Wind")) ) THEN
3783	windid = iv
3784	windndims = it
3785	ENDIF
3786	!
3787	ENDDO
3788	!
3789	! Treat some special cases and catch errors
3790	!
3791	IF ( varid < 0 ) THEN
3792	!
3793	! If we requested a wind component but did not find it, it could be that only the module is available.
3794	! If that is the case, then we use the module (windid) for the U component and set V top zero.
3795	!
3796	IF ( INDEX(varname, "Wind_E") > 0 ) THEN
3797	varid = windid
3798	ndims = windndims
3799	windzero = .FALSE.
3800	ELSE IF ( INDEX(varname, "Wind_N") > 0 ) THEN
3801	windzero = .TRUE.
3802	ELSE
3803	CALL ipslerr (3,'forcing_varforslab',"Could not find variable",varname," in file.")
3804	ENDIF
3805	ENDIF
3806	!
3807	! If there is some variable to be read then do it
3808	!
3809	IF ( .NOT. windzero ) THEN
3810	!
3811	! Get the attributes needed for intepretating the data
3812	!
3813	! First get the cell method used for this variable
3814	iret = NF90_GET_ATT(force_id(fileindex), varid, 'cell_methods', cellmethod)
3815	IF (iret /= NF90_NOERR) THEN
3816	! If the attribute is not found then we set a reasonable default : instantaneous and centered.
3817	cellmethod="time: instantaneous"
3818	ENDIF
3819	!
3820	!
3821	! Getsize of data to be read from the netCDF file
3822	!
3823	!
3824	IF ( compressed ) THEN
3825	!
3826	IF ( ndims == 2 ) THEN
3827	start = (/ncdfstart,timestart,0,0/)
3828	count = (/ncdfcount,timecount,0,0/)
3829	ELSE IF ( ndims == 3 ) THEN
3830	start = (/ncdfstart,1,timestart,0/)
3831	count = (/ncdfcount,1,timecount,0/)
3832	ELSE
3833	CALL ipslerr (3,'forcing_varforslab',"Compressed variable : ",varname,&
3834	& " does not have a compatible number of dimensions.")
3835	ENDIF
3836	!
3837	iret = NF90_GET_VAR(force_id(fileindex), varid, tmp_slab, start, count)
3838	IF (iret /= NF90_NOERR) THEN
3839	CALL ipslerr (3,'forcing_varforslab',"Could not read from file variable : ",varname," Compressed in the file.")
3840	ENDIF
3841	!
3842	! Zoom into the data and put it in the right place in the slab of data.
3843	!
3844	CALL forcing_reindex(ncdfcount, timecount, tmp_slab, nbland_loc, slab_size, data, inslabpos, reindex_loc)
3845	ELSE
3846	!
3847	IF ( ndims == 3 ) THEN
3848	start = (/1,1,timestart,0/)
3849	count = (/iim_glo,jjm_glo,timecount,0/)
3850	ELSE IF (ndims == 4 ) THEN
3851	start = (/1,1,1,timestart/)
3852	count = (/iim_glo,jjm_glo,1,timecount/)
3853	ELSE
3854	CALL ipslerr (3,'forcing_varforslab',"Full lat Lon variable : ",varname,&
3855	& " does not have a compatible number of dimensions.")
3856	ENDIF
3857	!
3858	iret = NF90_GET_VAR(force_id(fileindex), varid, tmp_slab2d, start, count)
3859	IF (iret /= NF90_NOERR) THEN
3860	WRITE(,) TRIM(NF90_STRERROR(iret))
3861	WRITE(,) "File =", fileindex, "Size =", SIZE(tmp_slab2d,DIM=1), SIZE(tmp_slab2d,DIM=2), SIZE(tmp_slab2d,DIM=3)
3862	WRITE(,) "Start :", start(1:3)
3863	WRITE(,) "Count :", count(1:3)
3864	CALL ipslerr (3,'forcing_varforslab',"Could not read from file variable : ",varname," Not compressed.")
3865	ENDIF
3866	!
3867	! Zoom into the data and put it in the right place in the slab of data.
3868	!
3869	CALL forcing_reindex(iim_glo, jjm_glo, timecount, tmp_slab2d, nbland_loc, slab_size, data, inslabpos, reindex2d_loc)
3870	ENDIF
3871	ELSE
3872	cellmethod="time: instantaneous"
3873	DO it=0,timecount-1
3874	data(:,inslabpos+it) = zero
3875	ENDDO
3876	ENDIF
3877	!
3878	END SUBROUTINE forcing_varforslab
3879
3880	!! =============================================================================================================================
3881	!! SUBROUTINE: forcing_attributetimeaxe
3882	!!
3883	!>\BRIEF Find the time axis which corresponds to the variable at hand.
3884	!!
3885	!! DESCRIPTION: We interpret the cell_method provided in the netCDF file so that
3886	!! we can determine how we need to understand the values we read.
3887	!!
3888	!! \n
3889	!_ ==============================================================================================================================
3890
3891	SUBROUTINE forcing_attributetimeaxe(cellmethod, timeindex)
3892	!
3893	! We will analyse the time axis of the cell method found in the NetCDF file in order to
3894	! attribute the correct time axis to this variable.
3895	!
3896	CHARACTER(LEN=*), INTENT(in) :: cellmethod
3897	INTEGER(i_std), INTENT(out) :: timeindex
3898	!
3899	INTEGER(i_std) :: itax, timepos, pos, lentime, itbase, im
3900	CHARACTER(LEN=20) :: TARGET, tmpstr
3901	CHARACTER(LEN=80) :: method
3902	!
3903	! Clean the string to delete spaces in front of ":" and "("
3904	!
3905	method = cellmethod
3906	DO WHILE ( INDEX(method," :") > 0 )
3907	pos = INDEX(method," :")
3908	method = method(1:pos-1)//method(pos+1:LEN_TRIM(method))
3909	ENDDO
3910	DO WHILE ( INDEX(method,"( ") > 0 )
3911	pos = INDEX(method,"( ")
3912	method = method(1:pos)//method(pos+2:LEN_TRIM(method))
3913	ENDDO
3914	!
3915	! Go through all the time axes we have to find the right one.
3916	!
3917	timeindex=0
3918	DO itax=1,nbtax
3919	!
3920	itbase=(itax-1)*nbtmethods
3921	! find the time axis name in the cell method
3922	TARGET = TRIM(time_axename(itbase+1))//":"
3923	timepos = INDEX(method,TRIM(TARGET))
3924	!
3925	! If we found the time axis then we look for the method with a time position description
3926	! which is expected to be between parenthesis. For instance : mean(end)
3927	!
3928	IF ( timepos > 0 ) THEN
3929	!
3930	lentime=LEN_TRIM(time_axename(itbase+1))
3931	tmpstr = method(lentime+2:LEN_TRIM(method))
3932	!
3933	! If there is ":" then there is information for another axis which needs to be deleted
3934	!
3935	IF ( INDEX(tmpstr,":") > 0 ) THEN
3936	tmpstr = tmpstr(1:INDEX(tmpstr,":")-1)
3937	ENDIF
3938	!
3939	! Now that we have found a time axis see if we have between parenthesis a position
3940	! on that time avis.
3941	!
3942	! Look for a "("
3943	IF ( INDEX(tmpstr, "(") > 0 ) THEN
3944	DO im=1,nbtmethods
3945	TARGET = "("//TRIM(time_cellmethod(itbase+im))
3946	timepos = INDEX(tmpstr,TRIM(TARGET))
3947	IF ( timepos > 0 ) THEN
3948	timeindex = itbase+im
3949	ENDIF
3950	ENDDO
3951	!
3952	! If there is no "(" then we have to find the centered axis.
3953	ELSE
3954	DO im=1,nbtmethods
3955	IF ( INDEX(time_cellmethod(itbase+im), "cent") > 0 ) THEN
3956	timeindex = itbase+im
3957	ENDIF
3958	ENDDO
3959	ENDIF
3960	!
3961	! The name of the time axis was found bu no method could be identified
3962	!
3963	IF ( timeindex < 1 ) THEN
3964	CALL ipslerr (3,'forcing_attributetimeaxe',"Found a time axis name but could not identify method.", &
3965	"This should not happen !"," ")
3966	ENDIF
3967	!
3968	ELSE
3969	! Continue in loop over nbtax
3970	ENDIF
3971	ENDDO
3972	!
3973	! Should no corresponding time axis name be found,
3974	! then we use the first centered one.
3975	!
3976	itax=1
3977	DO WHILE ( timeindex < 1 )
3978	IF ( INDEX(time_cellmethod(itax), "cent") > 0 ) THEN
3979	timeindex = itax
3980	ELSE
3981	itax = itax + 1
3982	ENDIF
3983	ENDDO
3984	!
3985	END SUBROUTINE forcing_attributetimeaxe
3986
3987	!! =============================================================================================================================
3988	!! SUBROUTINE: forcing_filenamecheck
3989	!!
3990	!>\BRIEF Check the list of files obtained from the calling program.
3991	!!
3992	!! DESCRIPTION: A small code to check the forcing files. They have to be NetCDF (i.e. .nc termination) and
3993	!! we dont keep files that appear more than once in the list.
3994	!!
3995	!! \n
3996	!_ ==============================================================================================================================
3997
3998	SUBROUTINE forcing_filenamecheck(filenames_in, nb_files)
3999	!
4000	! A small code to check the forcing files. They have to
4001	! be NetCDF (i.e. .nc termination) and we dont keep files
4002	! that appear more than once in the list.
4003	!
4004	!
4005	! INPUT
4006	!
4007	CHARACTER(LEN=*), DIMENSION(:), INTENT(in) :: filenames_in
4008	INTEGER(i_std), INTENT(out) :: nb_files
4009	!
4010	! LOCAL
4011	!
4012	INTEGER(i_std) :: ii, is, sizein
4013	LOGICAL :: notfound
4014	!
4015	sizein = SIZE(filenames_in)
4016	IF ( sizein > 0 ) THEN
4017	IF ( ALLOCATED(forfilename) ) THEN
4018	DEALLOCATE(forfilename)
4019	ENDIF
4020	ALLOCATE(forfilename(sizein))
4021	nb_files=0
4022	ELSE
4023	CALL ipslerr (3,'forcing_filenamecheck',"List of forcing files is empty.","Please check your run.def file."," ")
4024	ENDIF
4025	!
4026	DO ii=1,sizein
4027	IF ( INDEX(filenames_in(ii), '.nc') > 0 ) THEN
4028	IF ( nb_files == 0 ) THEN
4029	nb_files = nb_files+1
4030	forfilename(nb_files)=TRIM(ADJUSTL(filenames_in(ii)))
4031	ELSE
4032	notfound=.TRUE.
4033	DO is=1,nb_files
4034	IF ( INDEX(TRIM(filenames_in(ii)), TRIM(ADJUSTL(forfilename(is)))) > 0 ) notfound=.FALSE.
4035	ENDDO
4036	IF ( notfound ) THEN
4037	nb_files = nb_files+1
4038	forfilename(nb_files)=TRIM(adjustl(filenames_in(ii)))
4039	ENDIF
4040	ENDIF
4041	ELSE
4042	!!! This is not a NetCDF file, so we ignore it
4043	ENDIF
4044	ENDDO
4045	!
4046	!
4047	END SUBROUTINE forcing_filenamecheck
4048
4049	!! =============================================================================================================================
4050	!! FUNCTION: lowercase, FindMinimum, Swap
4051	!!
4052	!>\BRIEF Help functions for the forcing_tools module.
4053	!!
4054	!! DESCRIPTION:
4055	!!
4056	!! \n
4057	!_ ==============================================================================================================================
4058
4059	FUNCTION lowercase(strIn) RESULT(strOut)
4060	! Adapted from http://www.star.le.ac.uk/~cgp/fortran.html (25 May 2012)
4061
4062	IMPLICIT NONE
4063
4064	CHARACTER(len=*), INTENT(in) :: strIn
4065	CHARACTER(len=LEN(strIn)) :: strOut
4066	INTEGER :: i,j
4067
4068	DO i = 1, LEN(strIn)
4069	j = IACHAR(strIn(i:i))
4070	IF (j>= IACHAR("A") .AND. j<=IACHAR("Z") ) THEN
4071	strOut(i:i) = ACHAR(IACHAR(strIn(i:i))+32)
4072	ELSE
4073	strOut(i:i) = strIn(i:i)
4074	END IF
4075	END DO
4076
4077	END FUNCTION lowercase
4078	!
4079	! Some help function found on Internet : http://www.cs.mtu.edu/~shene/COURSES/cs201/NOTES/chap08/sorting.f90
4080	!
4081	! --------------------------------------------------------------------
4082	! INTEGER FUNCTION FindMinimum():
4083	! This function returns the location of the minimum in the section
4084	! between Start and End.
4085	! --------------------------------------------------------------------
4086	INTEGER FUNCTION FindMinimum(x, Start, END)
4087	IMPLICIT NONE
4088	REAL(r_std), DIMENSION(1:), INTENT(IN) :: x
4089	INTEGER(i_std), INTENT(IN) :: Start, END
4090	REAL(r_std) :: Minimum
4091	INTEGER(i_std) :: Location
4092	INTEGER(i_std) :: i
4093
4094	Minimum = x(Start) ! assume the first is the min
4095	Location = Start ! record its position
4096	DO i = Start+1, END ! start with next elements
4097	IF (x(i) < Minimum) THEN ! if x(i) less than the min?
4098	Minimum = x(i) ! Yes, a new minimum found
4099	Location = i ! record its position
4100	ENDIF
4101	END DO
4102	FindMinimum = Location ! return the position
4103	END FUNCTION FindMinimum
4104	! --------------------------------------------------------------------
4105	! SUBROUTINE Swap():
4106	! This subroutine swaps the values of its two formal arguments.
4107	! --------------------------------------------------------------------
4108	SUBROUTINE Swap(a, b)
4109	IMPLICIT NONE
4110	REAL(r_std), INTENT(INOUT) :: a, b
4111	REAL(r_std) :: Temp
4112
4113	Temp = a
4114	a = b
4115	b = Temp
4116	END SUBROUTINE Swap
4117	! --------------------------------------------------------------------
4118	! SUBROUTINE Sort():
4119	! This subroutine receives an array x() and sorts it into ascending
4120	! order.
4121	! --------------------------------------------------------------------
4122	SUBROUTINE Sort(x, Size)
4123	IMPLICIT NONE
4124	REAL(r_std), DIMENSION(1:), INTENT(INOUT) :: x
4125	INTEGER(i_std), INTENT(IN) :: Size
4126	INTEGER(i_std) :: i
4127	INTEGER(i_std) :: Location
4128
4129	DO i = 1, Size-1 ! except for the last
4130	Location = FindMinimum(x, i, Size) ! find min from this to last
4131	CALL Swap(x(i), x(Location)) ! swap this and the minimum
4132	END DO
4133	END SUBROUTINE Sort
4134
4135
4136
4137	!! ================================================================================================================================
4138	!! SUBROUTINE : forcing_tools_clear
4139	!!
4140	!>\BRIEF Clear forcing_tools module
4141	!!
4142	!! DESCRIPTION : Deallocate memory and reset initialization variables to there original values
4143	!!
4144	!_ ================================================================================================================================
4145
4146	SUBROUTINE forcing_tools_clear
4147
4148	first_call_readslab = .TRUE.
4149	first_call_solarint = .TRUE.
4150	first_call_spreadprec = .TRUE.
4151
4152	END SUBROUTINE forcing_tools_clear
4153
4154
4155	END MODULE forcing_tools

Note: See TracBrowser for help on using the repository browser.

Download in other formats: