Context Navigation

source: branches/ORCHIDEE_2_2/ORCHIDEE/src_driver/forcing_tools.f90 @ 7262

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

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

Download in other formats: