Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

forcing_tools.f90 @ 8367

Last change on this file since 8367 was 6717, checked in by josefine.ghattas, 4 years ago

Integrated relevant corrections done in the trunk into branch ORCHIDEE_3, between creation of the branch until 6708. It concerns following changeset: [6700], [6695], [6657], [6652], [6651]

Following command has been used:

svn merge -c 6651,6652,6657,6695,6700 svn://forge.ipsl.jussieu.fr/orchidee/trunk/ORCHIDEE --accept postpone

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

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

Context Navigation

source: branches/ORCHIDEE_3_CMIP6/ORCHIDEE/src_driver/forcing_tools.f90 @ 8367

Download in other formats: