Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

forcing_tools.f90 @ 7541

Last change on this file since 7541 was 7541, checked in by fabienne.maignan, 2 years ago
Zhang publication on coupling factor
File size: 185.6 KB

Line
1	! ==============================================================================================================================\n
2	! MODULE forcing_tools : This module concentrates on the temporal interpolation of the forcing for ORCHIDEE.
3	! It provides basic service for the grid when this is provided in the forcing file. The main
4	! work for the grid is done in glogrid.f90. The approach of forcing_tools to handle the time
5	! aspect of the forcing is to read as many time steps as possible in memory and then
6	! interpolate that to the time interval requested by the calling program.
7	! The data is read on root_proc but then distributed over all processors according to the
8	! domain decomposition of ORCHIDEE. This allows to be more efficient in the temporal interpolation.
9	! It is important to keep in mind that forcing_tools works on time intervals. So the request for data
10	! of ORCHIDEE as to be for an interval and the forcing file needs to have a description of the time interval
11	! over which the forcing is valid.
12	! The general description of how the attributes needed in the netCDF file for describing the cell_methods
13	! for time are provided in this document :
14	! https://forge.ipsl.jussieu.fr/orchidee/attachment/wiki/Documentation/Forcings/Description_Forcing_Files.pdf
15	!
16	! The most important routines of foring_tools are forcing_open and forcing_getvalues
17	!
18	! forcing_integration_time : Computes the interval over which the simulation should be carried out.
19	! forcing_open : Opens the forcing files and extract the main information.
20	! forcing_getvalues : Gets the forcing data for a time interval.
21	! forcing_close : Closes the forcing file
22	! forcing_printdate : A tool to print the dates in human readable form.
23	! forcing_printpoint : Print the values for a given point in time.
24	! forcing_givegridsize : Allows other programs to get the dimensions of the forcing grid.
25	! forcing_getglogrid : Allows other programs to get the spatial grid of the forcing.
26	! forcing_givegrid : Returns the description of the grid.
27	! forcing_zoomgrid : Extract a sub-region of the forcing grid.
28	!
29	! CONTACT : jan.polcher@lmd.jussieu.fr
30	!
31	! LICENCE : IPSL (2016)
32	! This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC
33	!
34	!>\BRIEF
35	!!
36	!! RECENT CHANGE(S): None
37	!!
38	!! REFERENCE(S) : None
39	!!
40	!_ ================================================================================================================================
41	!!
42	MODULE forcing_tools
43	!
44	USE defprec
45	USE netcdf
46	!
47	USE ioipsl
48	USE constantes
49	USE solar
50	!
51	USE mod_orchidee_para
52	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 SPREAD_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 = SPREAD_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 SPREAD_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('SPREAD_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 SPREAD_PREC_SEC is smaller than
1049	!Config the forcing time step. Should the user set SPREAD_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 'SPREAD_PREC' not be present in the run.def
1054	!
1055	time_to_spread = forcing_tstep_ave/2.0
1056	!
1057	CALL getin_p('SPREAD_PREC_SEC', time_to_spread)
1058	!
1059	ELSE
1060	time_to_spread = dt_sechiba_keep * nb_spread
1061	ENDIF
1062	!
1063	! Add the option to introduce a continuous distribution of precipitation
1064	!
1065	!Config Key = SPREAD_PREC_CONT
1066	!Config Desc = Take into account precipitation on the next forcing step for spreading it.
1067	!Config Def = FALSE
1068	!Config Help = This allows to extend the spreading of rainfall to the entire forcing
1069	!Config should it be raining on the following time step. This ensures that if it rains
1070	!Config in two consecutive forcing time step at least during the first period rainfall
1071	!Config will be continuous. This avoids the spiky nature of rainfall in periods of
1072	!Config prolonged rainfall.
1073	!Config Units = -
1074	!
1075	! This is the default should 'SPREAD_PREC_CONT' not be present in the run.def
1076	!
1077	spreadprec_cont = .FALSE.
1078	!
1079	CALL getin_p('SPREAD_PREC_CONT', spreadprec_cont)
1080	!
1081	! Do some verifications on the information read from run.def
1082	!
1083	IF ( time_to_spread > forcing_tstep_ave) THEN
1084	time_to_spread = forcing_tstep_ave
1085	ELSE IF ( time_to_spread <= 0 ) THEN
1086	time_to_spread = forcing_tstep_ave
1087	ENDIF
1088	!
1089	first_call_spreadprec = .FALSE.
1090	!
1091	ENDIF
1092	!
1093	! First test that we have the right time interval from the forcing to spread the precipitation
1094	!
1095	IF ( time_int(1) >= timebnd_central(1,1) .AND. time_int(2) <= timebnd_central(slab_size,2)) THEN
1096	!
1097	! Create a mask so that MINLOC does not look outside of the valid interval of time_central
1098	!
1099	mask(1:slab_size) = .TRUE.
1100	!
1101	! To get better precision on the time difference we get a common offset to substract
1102	!
1103	offset = INT(forcingstartdate)
1104	!
1105	! In principle 3 time steps can contribute to the time step closest to the center of the forcing interval
1106	!
1107	imin = MINLOC( ABS(time_central(1:slab_size)-(time_int(1)+time_int(2))/2.0), mask(1:slab_size) )
1108	tind(1) = MAX(imin(1)-1,1)
1109	tind(2) = imin(1)
1110	tind(3) = MIN(imin(1)+1,slab_size)
1111	IF (imin(1)+1 > slab_size) THEN
1112	WRITE(,) "We have a problem here imin(1)+1,slab_size ", imin(1)+1,slab_size
1113	WRITE(,) "Interval : ", time_int(1),time_int(2)
1114	ENDIF
1115	!
1116	tspread(:) = time_to_spread
1117	!
1118	DO i=1, nbpoint_proc
1119	!
1120	IF ( spreadprec_cont ) THEN
1121	!
1122	! If the next time step has also rainfall, then time_to_spread becomes the length of the forcing time step.
1123	!
1124	IF ( rainf_slab(i,tind(3)) > zero .OR. snowf_slab(i,tind(3)) > zero) THEN
1125	tspread(i) = forcing_tstep_ave
1126	ELSE
1127	tspread(i) = time_to_spread
1128	ENDIF
1129	ELSE
1130	! Default behavious
1131	tspread(i) = time_to_spread
1132	ENDIF
1133	!
1134	! Do we need to take some rain from the previous time step ?
1135	!
1136	!! Time computation is not better than 1/1000 seconds
1137	IF ( time_int(1) < timebnd_central(tind(2),1) .AND. preciptime_slab(i,tind(1)) < (tspread(i)-0.001) ) THEN
1138	dt = ((timebnd_central(tind(2),1)-offset)-(time_int(1)-offset))*one_day
1139	ft(1) = MIN(tspread(i) - preciptime_slab(i,tind(1)), dt)/tlen
1140	ELSE
1141	ft(1) = zero
1142	ENDIF
1143	!
1144	! Is there still some rain to spread from the current forcing time step ?
1145	!
1146	!! Time computation is not better than 1/1000 seconds
1147	IF (preciptime_slab(i,tind(2)) < (tspread(i)-0.001) ) THEN
1148	left = MAX(time_int(1), timebnd_central(tind(2),1))
1149	right = MIN(time_int(2),timebnd_central(tind(2),2))
1150	dt = ((right-offset)-(left-offset))*one_day
1151	ft(2) = MIN(tspread(i) - preciptime_slab(i,tind(2)), dt)/tlen
1152	ELSE
1153	ft(2) = zero
1154	ENDIF
1155	!
1156	! Do we need to take some rain from the next time step ?
1157	!
1158	!! Time computation is not better than 1/1000 seconds
1159	IF ( time_int(2) > timebnd_central(tind(2),2) .AND. preciptime_slab(i,tind(3)) < (tspread(i)-0.001) ) THEN
1160	dt = ((time_int(2)-offset)-(timebnd_central(tind(2),2)-offset))*one_day
1161	ft(3) = MIN(tspread(i) - preciptime_slab(i,tind(3)), dt)/tlen
1162	ELSE
1163	ft(3) = zero
1164	ENDIF
1165	!
1166	! Do the actual calculation
1167	!
1168	rainf(i) = (rainf_slab(i,tind(1)) * forcing_tstep_ave * ft(1) + &
1169	& rainf_slab(i,tind(2)) * forcing_tstep_ave * ft(2) + &
1170	& rainf_slab(i,tind(3)) * forcing_tstep_ave * ft(3))*tlen/tspread(i)
1171
1172	snowf(i) = (snowf_slab(i,tind(1)) * forcing_tstep_ave * ft(1) + &
1173	& snowf_slab(i,tind(2)) * forcing_tstep_ave * ft(2) + &
1174	& snowf_slab(i,tind(3)) * forcing_tstep_ave * ft(3))*tlen/tspread(i)
1175	!
1176	! Update the time over which we have already spread the rainf
1177	!
1178	preciptime_slab(i,tind(1)) = preciptime_slab(i,tind(1)) + tlen * ft(1)
1179	preciptime_slab(i,tind(2)) = preciptime_slab(i,tind(2)) + tlen * ft(2)
1180	preciptime_slab(i,tind(3)) = preciptime_slab(i,tind(3)) + tlen * ft(3)
1181	!
1182	ENDDO
1183	ELSE
1184	WRITE(numout,*) "Time interval toward which we will interpolate : ", time_int
1185	WRITE(numout,*) "Limits of the time slab we have : ", timebnd_central(1,1), timebnd_central(slab_size,2)
1186	CALL forcing_printdate(time_int(1), "Start of target time interval.")
1187	CALL forcing_printdate(time_int(2), "End of target time interval.")
1188	CALL forcing_printdate(timebnd_central(1,1), "Start of time slab we have.")
1189	CALL forcing_printdate(timebnd_central(slab_size,2), "End of time slab we have.")
1190	CALL ipslerr (3,'forcing_spreadprec', 'The sitation should not occur Why are we here ?',&
1191	& "","")
1192	ENDIF
1193
1194	END SUBROUTINE forcing_spreadprec
1195
1196	!! =============================================================================================================================
1197	!! SUBROUTINE: forcing_solarint
1198	!!
1199	!>\BRIEF Interpolates incoming solar radiation to the interval requested.
1200	!!
1201	!! DESCRIPTION: The interpolation here takes into account the variation of the solar zenith angle
1202	!! to ensure the diurnal cycle of solar radiation is as well represented as possible.
1203	!!
1204	!! \n
1205	!_ ==============================================================================================================================
1206	SUBROUTINE forcing_solarint(time_int_in, tlen, timebnd_in, time_cent_in, iim, jjm, lon, lat, swdown, solarangle)
1207	!
1208	! ARGUMENTS
1209	!
1210	REAL(r_std), INTENT(in) :: time_int_in(2) ! Time interval for which we will compute radiation
1211	REAL(r_std), INTENT(in) :: tlen ! size of time interval in seconds (time step !)
1212	REAL(r_std), INTENT(in) :: timebnd_in(:,:) ! Time interval over which the read data is valid
1213	REAL(r_std), INTENT(in) :: time_cent_in(:) ! Center of the time interval
1214	INTEGER(i_std), INTENT(in) :: iim, jjm ! Size of 2D domain
1215	REAL(r_std), INTENT(in) :: lon(iim,jjm), lat(iim,jjm) ! Longitude and latitude
1216	REAL(r_std), INTENT(out) :: swdown(:), solarangle(:) ! interpolated downward solar radiation and corresponding
1217	! ! solar angle.
1218	!
1219	! LOCAL SAVED
1220	!
1221	LOGICAL, SAVE :: first_call_solarint=.TRUE.
1222	REAL(r_std), SAVE :: solaryearstart
1223	INTEGER(i_std), SAVE :: split, split_max
1224	REAL(r_std), SAVE :: last_time
1225	!
1226	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: mean_sinang
1227	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:,:) :: sinangles
1228	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: time_angles
1229	! Dusk-dawn management
1230	REAL(r_std), SAVE :: dusk_angle
1231	!
1232	! LOCAL - temporary
1233	!
1234	REAL(r_std) :: time_int(2)
1235	REAL(r_std) :: timebnd(slab_size_max,2)
1236	REAL(r_std) :: time_cent(slab_size_max)
1237	INTEGER(i_std) :: year, month, day, hours, minutes
1238	REAL(r_std) :: sec
1239	!
1240	REAL(r_std) :: mean_sol, split_time
1241	REAL(r_std) :: julian, julian_tmp
1242	REAL(r_std) :: sinang(iim,jjm)
1243	INTEGER(i_std) :: is, i, ii, jj, imin(1), tmin(1)
1244	LOGICAL, ALLOCATABLE, DIMENSION(:) :: mask
1245	!
1246	!
1247	IF ( .NOT. ALLOCATED(mask) ) THEN
1248	ALLOCATE(mask(slab_size_max))
1249	mask(:) = .FALSE.
1250	ENDIF
1251
1252	IF ( first_call_solarint ) THEN
1253	!
1254	! Ensure the offset is on the 1st of Januray of the current years so that we do not
1255	! perturbe the solar angle calculations.
1256	!
1257	CALL ju2ymds (startdate, year, month, day, sec)
1258	CALL ymds2ju (year, 1, 1, 0.0, solaryearstart)
1259	!
1260	last_time = -9999.0
1261	!
1262	ALLOCATE(mean_sinang(iim,jjm))
1263	mean_sinang(:,:) = 0.0
1264	!
1265	split = NINT(forcing_tstep_ave/tlen)
1266	!
1267	! Allow for more space than estimated with the size of the first time step.
1268	!
1269	ALLOCATE(time_angles(split*2))
1270	time_angles(:) = 0.0
1271	!
1272	ALLOCATE(sinangles(iim,jjm,split*2))
1273	sinangles(:,:,:) = 0.0
1274	!
1275	dusk_angle=0.01
1276	!
1277	first_call_solarint = .FALSE.
1278	!
1279	split = 0
1280	!
1281	ENDIF
1282	!
1283	! Shift the input dates in order to gain in precision for the time calculations
1284	!
1285	time_int(:) = time_int_in(:)-INT(solaryearstart)
1286	time_cent(1:slab_size) = time_cent_in(1:slab_size)-INT(solaryearstart)
1287	timebnd(1:slab_size,1) = timebnd_in(1:slab_size,1)-INT(solaryearstart)
1288	timebnd(1:slab_size,2) = timebnd_in(1:slab_size,2)-INT(solaryearstart)
1289	!
1290	! Create a mask so that MINLOC does not look outside of the valid interval of time_central
1291	!
1292	mask(1:slab_size) = .TRUE.
1293	!
1294	! Locate the time step in the SLAB at hand
1295	!
1296	imin = MINLOC( ABS(time_cent(1:slab_size)-(time_int(1)+time_int(2))/2.0), mask(1:slab_size) )
1297	!
1298	! Compute all the angels we will encounter for the current forcing interval
1299	!
1300	IF ( last_time .NE. timebnd(imin(1),1) ) THEN
1301	!
1302	! Verify that we have used all the angles of the previous decomposition of the forcing
1303	! time step.
1304	!
1305	IF ( split .NE. 0 ) THEN
1306	!
1307	WRITE(numout,*) "The forcing has a time step of : ", forcing_tstep_ave
1308	WRITE(numout,*) "The model is configured to run with a time step of : ", tlen
1309	WRITE(numout,*) "We are left with split = ", split, " starting from ", split_max
1310	!
1311	CALL ipslerr (3,'forcing_solarint',"The decomposition of solar downward radiation of the forcing file over the model",&
1312	& "has failed. This means the average of the solar radiation over the forcing time step is not conserved.",&
1313	& "This can be caused by a time step repeated twice.")
1314	ENDIF
1315	!
1316	! Compute the number of time steps the model will put in the current interval of forcing data.
1317	!
1318	split = 0
1319	julian_tmp = (time_int(1)+time_int(2))/2.0
1320	split_time = julian_tmp+split*tlen/one_day
1321	tmin = MINLOC( ABS(time_cent(1:slab_size) - split_time), mask(1:slab_size))
1322	DO WHILE ( tmin(1) .EQ. imin(1) .AND. split_time .LE. timebnd(slab_size,2) )
1323	split = split + 1
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	ENDDO
1327	!
1328	mean_sinang(:,:) = 0.0
1329	time_angles(:) = 0.0
1330	!
1331	DO is = 1,split
1332	!
1333	julian = julian_tmp + (is-1)*tlen/one_day
1334	!
1335	! This call should be better at it allows to compute the difference between the
1336	! current date and a reference time to higher precision. But it produces noisy
1337	! SWdown fluxes !
1338	!! CALL solarang (julian, solaryearstart, iim, jjm, lon, lat, sinang)
1339	! The old approach.
1340	CALL solarang (julian+INT(solaryearstart), solaryearstart, iim, jjm, lon, lat, sinang)
1341	!
1342	! During the dusk,dawn period maintain a minimum angle to take into account the
1343	! diffuse radiation which starts before the sun is over the horizon.
1344	!
1345	DO ii=1,iim
1346	DO jj=1,jjm
1347	IF ( sinang(ii,jj) > zero .AND. sinang(ii,jj) < dusk_angle ) THEN
1348	sinang(ii,jj) = dusk_angle
1349	ENDIF
1350	mean_sinang(ii,jj) = mean_sinang(ii,jj)+sinang(ii,jj)
1351	ENDDO
1352	ENDDO
1353	!
1354	! Save the solar angle information for later use. That is when we have the target time we will
1355	! look in this table the angle we have forseen.
1356	!
1357	time_angles(is) = julian
1358	sinangles(:,:,is) = sinang(:,:)
1359	!
1360	ENDDO
1361	!
1362	mean_sinang(:,:) = mean_sinang(:,:)/split
1363	last_time = timebnd(imin(1),1)
1364	split_max = split
1365	!
1366	ENDIF
1367	!
1368	! For the current timle step get the time of the closest pre-computed solar angle.
1369	!
1370	julian = (time_int(1)+time_int(2))/2.0
1371	tmin = MINLOC(ABS(julian-time_angles(1:split_max)), mask)
1372	sinang(:,:) = sinangles(:,:,tmin(1))
1373	! Remember that we have taken one value of the table for later verification
1374	split = split - 1
1375	!
1376	DO i=1, nbpoint_proc
1377	!
1378	jj = ((glolindex_proc(i)-1)/iim)+1
1379	ii = (glolindex_proc(i)-(jj-1)*iim)
1380	!
1381	IF ( mean_sinang(ii,jj) > zero ) THEN
1382	swdown(i) = swdown_slab(i,imin(1))*sinang(ii,jj)/mean_sinang(ii,jj)
1383	ELSE
1384	swdown(i) = zero
1385	ENDIF
1386	!
1387	! Why is this ??? Should we not take the angle corresponding to this time step ?? (Jan)
1388	!
1389	solarangle(i) = mean_sinang(ii,jj)
1390	!
1391	ENDDO
1392	!
1393	END SUBROUTINE forcing_solarint
1394	!!
1395	!! =============================================================================================================================
1396	!! SUBROUTINE: forcing_readslab
1397	!!
1398	!>\BRIEF Interface routine to read the data. This routine prepares the memory on each procesor and scatters the read data.
1399	!!
1400	!! DESCRIPTION:
1401	!!
1402	!! \n
1403	!_ ==============================================================================================================================
1404	SUBROUTINE forcing_readslab(time_int)
1405	!
1406	! This routine serves to interface with forcing_readslab_root and ensure that
1407	! the data is distributed correctly on all processors.
1408	!
1409	REAL(r_std), INTENT(in) :: time_int(2) !! The time interval over which the forcing is needed.
1410	!
1411	! Local
1412	!
1413	INTEGER(i_std) :: is
1414	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: tair_full, qair_full
1415	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: tairmin_full, tairmax_full
1416	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: rainf_full, snowf_full
1417	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: swdown_full, lwdown_full
1418	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: u_full, v_full
1419	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: ps_full, ztq_full, zuv_full
1420	REAL(r_std), ALLOCATABLE, DIMENSION(:) :: preciptime_tmp
1421	!
1422	! 1.0 Verify that for the slabs the memory is allocated for the variable
1423	! as well as its time axis.
1424	!
1425	IF ( .NOT. ALLOCATED(tair_slab) ) ALLOCATE(tair_slab(nbpoint_proc,slab_size))
1426	IF ( .NOT. ALLOCATED(time_tair) ) ALLOCATE(time_tair(slab_size))
1427	IF ( .NOT. ALLOCATED(timebnd_tair) ) ALLOCATE(timebnd_tair(slab_size,2))
1428	IF ( forcing_tstep_ave >= one_day/3.0) THEN
1429	IF ( .NOT. ALLOCATED(tairmax_slab) ) ALLOCATE(tairmax_slab(nbpoint_proc,slab_size))
1430	IF ( .NOT. ALLOCATED(tairmin_slab) ) ALLOCATE(tairmin_slab(nbpoint_proc,slab_size))
1431	ENDIF
1432	!
1433	IF ( .NOT. ALLOCATED(qair_slab) ) ALLOCATE(qair_slab(nbpoint_proc,slab_size))
1434	IF ( .NOT. ALLOCATED(time_qair) ) ALLOCATE(time_qair(slab_size))
1435	IF ( .NOT. ALLOCATED(timebnd_qair) ) ALLOCATE(timebnd_qair(slab_size,2))
1436	!
1437	IF ( .NOT. ALLOCATED(rainf_slab) ) ALLOCATE(rainf_slab(nbpoint_proc,slab_size))
1438	IF ( .NOT. ALLOCATED(snowf_slab) ) ALLOCATE(snowf_slab(nbpoint_proc,slab_size))
1439	IF ( .NOT. ALLOCATED(time_precip) ) ALLOCATE(time_precip(slab_size))
1440	IF ( .NOT. ALLOCATED(timebnd_precip) ) ALLOCATE(timebnd_precip(slab_size,2))
1441	IF ( .NOT. ALLOCATED(preciptime_slab) ) ALLOCATE(preciptime_slab(nbpoint_proc,slab_size))
1442	IF ( .NOT. ALLOCATED(preciptime_tmp) ) ALLOCATE(preciptime_tmp(slab_size))
1443	!
1444	IF ( .NOT. ALLOCATED(swdown_slab) ) ALLOCATE(swdown_slab(nbpoint_proc,slab_size))
1445	IF ( .NOT. ALLOCATED(time_swdown) ) ALLOCATE(time_swdown(slab_size))
1446	IF ( .NOT. ALLOCATED(timebnd_swdown) ) ALLOCATE(timebnd_swdown(slab_size,2))
1447	!
1448	IF ( .NOT. ALLOCATED(lwdown_slab) ) ALLOCATE(lwdown_slab(nbpoint_proc,slab_size))
1449	IF ( .NOT. ALLOCATED(time_lwdown) ) ALLOCATE(time_lwdown(slab_size))
1450	IF ( .NOT. ALLOCATED(timebnd_lwdown) ) ALLOCATE(timebnd_lwdown(slab_size,2))
1451	!
1452	IF ( .NOT. ALLOCATED(u_slab) ) ALLOCATE(u_slab(nbpoint_proc,slab_size))
1453	IF ( .NOT. ALLOCATED(time_u) ) ALLOCATE(time_u(slab_size))
1454	IF ( .NOT. ALLOCATED(timebnd_u) ) ALLOCATE(timebnd_u(slab_size,2))
1455	!
1456	IF ( .NOT. ALLOCATED(v_slab) ) ALLOCATE(v_slab(nbpoint_proc,slab_size))
1457	IF ( .NOT. ALLOCATED(time_v) ) ALLOCATE(time_v(slab_size))
1458	IF ( .NOT. ALLOCATED(timebnd_v) ) ALLOCATE(timebnd_v(slab_size,2))
1459	!
1460	IF ( .NOT. ALLOCATED(ps_slab) ) ALLOCATE(ps_slab(nbpoint_proc,slab_size))
1461	IF ( .NOT. ALLOCATED(time_ps) ) ALLOCATE(time_ps(slab_size))
1462	IF ( .NOT. ALLOCATED(timebnd_ps) ) ALLOCATE(timebnd_ps(slab_size,2))
1463	!
1464	IF ( .NOT. ALLOCATED(ztq_slab) ) ALLOCATE(ztq_slab(nbpoint_proc,slab_size))
1465	IF ( .NOT. ALLOCATED(zuv_slab) ) ALLOCATE(zuv_slab(nbpoint_proc,slab_size))
1466	!
1467	!
1468	IF ( is_root_prc) THEN
1469	!
1470	! Allocate ont he root processor the memory to receive the data from the file
1471	!
1472	ALLOCATE(tair_full(nbpoint_loc,slab_size))
1473	IF ( forcing_tstep_ave >= one_day/3.0) THEN
1474	ALLOCATE(tairmax_full(nbpoint_loc,slab_size))
1475	ALLOCATE(tairmin_full(nbpoint_loc,slab_size))
1476	ENDIF
1477	ALLOCATE(qair_full(nbpoint_loc,slab_size))
1478	ALLOCATE(rainf_full(nbpoint_loc,slab_size))
1479	ALLOCATE(snowf_full(nbpoint_loc,slab_size))
1480	ALLOCATE(swdown_full(nbpoint_loc,slab_size))
1481	ALLOCATE(lwdown_full(nbpoint_loc,slab_size))
1482	ALLOCATE(u_full(nbpoint_loc,slab_size))
1483	ALLOCATE(v_full(nbpoint_loc,slab_size))
1484	ALLOCATE(ps_full(nbpoint_loc,slab_size))
1485	ALLOCATE(ztq_full(nbpoint_loc,slab_size))
1486	ALLOCATE(zuv_full(nbpoint_loc,slab_size))
1487	!
1488	CALL forcing_readslab_root(time_int, &
1489	& tair_full, tairmax_full, tairmin_full, time_tair, timebnd_tair, &
1490	& qair_full, time_qair, timebnd_qair, &
1491	& rainf_full, snowf_full, time_precip, timebnd_precip, preciptime_tmp, &
1492	& swdown_full, time_swdown, timebnd_swdown, &
1493	& lwdown_full, time_lwdown, timebnd_lwdown, &
1494	& u_full, time_u, timebnd_u, &
1495	& v_full, time_v, timebnd_v, &
1496	& ps_full, time_ps, timebnd_ps, &
1497	& ztq_full, zuv_full)
1498	!
1499	ELSE
1500	!
1501	ALLOCATE(tair_full(1,1))
1502	IF ( forcing_tstep_ave >= one_day/3.0) THEN
1503	ALLOCATE(tairmax_full(1,1))
1504	ALLOCATE(tairmin_full(1,1))
1505	ENDIF
1506	ALLOCATE(qair_full(1,1))
1507	ALLOCATE(rainf_full(1,1))
1508	ALLOCATE(snowf_full(1,1))
1509	ALLOCATE(swdown_full(1,1))
1510	ALLOCATE(lwdown_full(1,1))
1511	ALLOCATE(u_full(1,1))
1512	ALLOCATE(v_full(1,1))
1513	ALLOCATE(ps_full(1,1))
1514	ALLOCATE(ztq_full(1,1))
1515	ALLOCATE(zuv_full(1,1))
1516	!
1517	ENDIF
1518	!
1519	! Broadcast the time information to all procs.
1520	!
1521	CALL bcast(slab_size)
1522	CALL bcast(time_tair)
1523	CALL bcast(timebnd_tair)
1524	CALL bcast(time_qair)
1525	CALL bcast(timebnd_qair)
1526	CALL bcast(time_precip)
1527	CALL bcast(timebnd_precip)
1528	CALL bcast(preciptime_tmp)
1529	CALL bcast(time_swdown)
1530	CALL bcast(timebnd_swdown)
1531	CALL bcast(time_lwdown)
1532	CALL bcast(timebnd_lwdown)
1533	CALL bcast(time_u)
1534	CALL bcast(timebnd_u)
1535	CALL bcast(time_v)
1536	CALL bcast(timebnd_v)
1537	CALL bcast(time_ps)
1538	CALL bcast(timebnd_ps)
1539	!
1540	! Scatter the slabs of data to all processors
1541	!
1542	IF ( landonly ) THEN
1543	CALL scatter(tair_full, tair_slab)
1544	IF ( forcing_tstep_ave >= one_day/3.0) THEN
1545	CALL scatter(tairmax_full, tairmax_slab)
1546	CALL scatter(tairmin_full, tairmin_slab)
1547	ENDIF
1548	CALL scatter(qair_full, qair_slab)
1549	CALL scatter(rainf_full, rainf_slab)
1550	CALL scatter(snowf_full, snowf_slab)
1551	CALL scatter(swdown_full, swdown_slab)
1552	CALL scatter(lwdown_full, lwdown_slab)
1553	CALL scatter(u_full, u_slab)
1554	CALL scatter(v_full, v_slab)
1555	CALL scatter(ps_full, ps_slab)
1556	CALL scatter(ztq_full, ztq_slab)
1557	CALL scatter(zuv_full, zuv_slab)
1558	ELSE
1559	tair_slab(:,:) = tair_full(:,:)
1560	IF ( forcing_tstep_ave >= one_day/3.0) THEN
1561	tairmax_slab(:,:) = tairmax_full(:,:)
1562	tairmin_slab(:,:) = tairmin_full(:,:)
1563	ENDIF
1564	qair_slab(:,:) = qair_full(:,:)
1565	rainf_slab(:,:) = rainf_full(:,:)
1566	snowf_slab(:,:) = snowf_full(:,:)
1567	swdown_slab(:,:) = swdown_full(:,:)
1568	lwdown_slab(:,:) = lwdown_full(:,:)
1569	u_slab(:,:) = u_full(:,:)
1570	v_slab(:,:) = v_full(:,:)
1571	ps_slab(:,:) = ps_full(:,:)
1572	ztq_slab(:,:) = ztq_full(:,:)
1573	zuv_slab(:,:) = zuv_full(:,:)
1574	ENDIF
1575	!
1576	!
1577	!
1578	DO is=1,nbpoint_proc
1579	preciptime_slab(is,:) = preciptime_tmp(:)
1580	ENDDO
1581	!
1582	! Clean-up to free the memory on the root processor.
1583	!
1584	IF ( ALLOCATED(tair_full) ) DEALLOCATE(tair_full)
1585	IF ( ALLOCATED(tairmax_full) ) DEALLOCATE(tairmax_full)
1586	IF ( ALLOCATED(tairmin_full) ) DEALLOCATE(tairmin_full)
1587	IF ( ALLOCATED(qair_full) ) DEALLOCATE(qair_full)
1588	IF ( ALLOCATED(rainf_full) ) DEALLOCATE(rainf_full)
1589	IF ( ALLOCATED(snowf_full) ) DEALLOCATE(snowf_full)
1590	IF ( ALLOCATED(swdown_full) ) DEALLOCATE(swdown_full)
1591	IF ( ALLOCATED(lwdown_full) ) DEALLOCATE(lwdown_full)
1592	IF ( ALLOCATED(u_full) ) DEALLOCATE(u_full)
1593	IF ( ALLOCATED(v_full) ) DEALLOCATE(v_full)
1594	IF ( ALLOCATED(ps_full) ) DEALLOCATE(ps_full)
1595	IF ( ALLOCATED(ztq_full) ) DEALLOCATE(ztq_full)
1596	IF ( ALLOCATED(zuv_full) ) DEALLOCATE(zuv_full)
1597	!
1598	END SUBROUTINE forcing_readslab
1599	!!
1600	!! =============================================================================================================================
1601	!! SUBROUTINE: forcing_readslab_root
1602	!!
1603	!>\BRIEF Routine which reads a slab of data from the netCDF file and writes it onto the memory.
1604	!!
1605	!! DESCRIPTION: It is important to read the next slab of data while still keeping an overlap so that
1606	!! interpolation can continue.
1607	!! It also attributes a time axis to each variable.
1608	!!
1609	!! \n
1610	!_ ==============================================================================================================================
1611
1612	SUBROUTINE forcing_readslab_root(time_int, &
1613	& tair, tairmax, tairmin, t_tair, tbnd_tair, &
1614	& qair, t_qair, tbnd_qair, &
1615	& rainf, snowf, t_prec, tbnd_prec, prectime, &
1616	& swdown, t_swdown, tbnd_swdown, &
1617	& lwdown, t_lwdown, tbnd_lwdown, &
1618	& u, t_u, tbnd_u, &
1619	& v, t_v, tbnd_v, &
1620	& ps, t_ps, tbnd_ps, &
1621	& ztq, zuv)
1622	!
1623	! Arguments
1624	!
1625	REAL(r_std), INTENT(in) :: time_int(2) !! The time interval over which the forcing is needed.
1626	!
1627	REAL(r_std), INTENT(out) :: tair(:,:), tairmax(:,:), tairmin(:,:)
1628	REAL(r_std), INTENT(out) :: t_tair(:)
1629	REAL(r_std), INTENT(out) :: tbnd_tair(:,:)
1630	!
1631	REAL(r_std), INTENT(out) :: qair(:,:)
1632	REAL(r_std), INTENT(out) :: t_qair(:)
1633	REAL(r_std), INTENT(out) :: tbnd_qair(:,:)
1634	!
1635	REAL(r_std), INTENT(out) :: rainf(:,:)
1636	REAL(r_std), INTENT(out) :: snowf(:,:)
1637	REAL(r_std), INTENT(out) :: t_prec(:)
1638	REAL(r_std), INTENT(out) :: tbnd_prec(:,:)
1639	REAL(r_std), INTENT(out) :: prectime(:)
1640	!
1641	REAL(r_std), INTENT(out) :: swdown(:,:)
1642	REAL(r_std), INTENT(out) :: t_swdown(:)
1643	REAL(r_std), INTENT(out) :: tbnd_swdown(:,:)
1644	!
1645	REAL(r_std), INTENT(out) :: lwdown(:,:)
1646	REAL(r_std), INTENT(out) :: t_lwdown(:)
1647	REAL(r_std), INTENT(out) :: tbnd_lwdown(:,:)
1648	!
1649	REAL(r_std), INTENT(out) :: u(:,:)
1650	REAL(r_std), INTENT(out) :: t_u(:)
1651	REAL(r_std), INTENT(out) :: tbnd_u(:,:)
1652	!
1653	REAL(r_std), INTENT(out) :: v(:,:)
1654	REAL(r_std), INTENT(out) :: t_v(:)
1655	REAL(r_std), INTENT(out) :: tbnd_v(:,:)
1656	!
1657	REAL(r_std), INTENT(out) :: ps(:,:)
1658	REAL(r_std), INTENT(out) :: t_ps(:)
1659	REAL(r_std), INTENT(out) :: tbnd_ps(:,:)
1660	!
1661	REAL(r_std), INTENT(out) :: ztq(:,:)
1662	REAL(r_std), INTENT(out) :: zuv(:,:)
1663	!
1664	! Local
1665	!
1666	INTEGER(i_std) :: iret, varid
1667	INTEGER(i_std) :: if, it
1668	INTEGER(i_std) :: tstart(3), tcount(3)
1669	INTEGER(i_std) :: imin(1), imax(1), firstif(1)
1670	INTEGER(i_std) :: nctstart, nctcount, inslabpos
1671	INTEGER(i_std) :: start_globtime, end_globtime
1672	INTEGER(i_std) :: timeid_tair, timeid_qair, timeid_precip, timeid_swdown
1673	INTEGER(i_std) :: timeid_lwdown, timeid_u, timeid_v, timeid_ps, timeid_tmp
1674	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: time_tmp
1675	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: rau
1676	CHARACTER(LEN=80) :: cellmethod
1677	!
1678	LOGICAL, SAVE :: first_call_readslab=.TRUE.
1679	!
1680	ALLOCATE(time_tmp(slab_size,nbtax))
1681	ALLOCATE(rau(nbpoint_loc,slab_size))
1682	!
1683	!
1684	! Catch any stupid utilisation of this routine.
1685	!
1686	IF ( .NOT. is_root_prc) THEN
1687	CALL ipslerr (3,'forcing_readslab_root',"Cannot run this routine o other procs than root.",&
1688	& "All the information on the forcing files is only lated on the root processor."," ")
1689	ENDIF
1690	!
1691	!Set some variables to zero
1692	!
1693	IF ( first_call_readslab ) THEN
1694	!
1695	preciptime(:) = 0
1696	!
1697	! If the first file is only there to provide the last time step of the previous year, we
1698	! do not need to read all. We will start 2 forcing time steps before the start of the first
1699	! time interval requested.
1700	!
1701	imin=MINLOC(ABS(time(:,1)-time_int(1)))
1702	current_offset = MAX(imin(1)-2,1)
1703	!
1704	first_call_readslab = .FALSE.
1705	write(numout, *) "first_call_readslab in forcing_readslab_root"
1706	!
1707	ELSE
1708	!
1709	! Put back the cummulated time of rainfall into the global array
1710	!
1711	preciptime(position_slab(1):position_slab(2)) = preciptime(position_slab(1):position_slab(2)) + &
1712	& prectime(1:slab_size)
1713	!
1714	! Compute new offset
1715	!
1716	current_offset = position_slab(2)-2
1717	write(numout, *) "first_call_readslab in forcing_readslab_root 22"
1718	!
1719	ENDIF
1720
1721	!
1722	! Check that the slab size is not too large
1723	!
1724	IF ( (current_offset-1)+slab_size > nb_forcing_steps) THEN
1725	slab_size = nb_forcing_steps - (current_offset-1)
1726	ENDIF
1727	!
1728	! 1.1 Check that the slab we have to read exists
1729	!
1730	IF ( slab_size > 0 ) THEN
1731	!
1732	start_globtime = current_offset
1733	end_globtime = (current_offset-1)+slab_size
1734	inslabpos=1
1735	WRITE(,) ">> Reading from global position ", start_globtime, "up to ", end_globtime
1736	write(numout,*) time_sourcefile
1737	!
1738	DO if=MINVAL(time_sourcefile(start_globtime:end_globtime)),MAXVAL(time_sourcefile(start_globtime:end_globtime))
1739	!
1740	! Get position of the part of the global time axis we need to read from this file
1741	!
1742	firstif = MINLOC(ABS(time_sourcefile-if))
1743	! start = distance of start_globtime or nothing + 1 to follow netCDF convention.
1744	nctstart = MAX((start_globtime-firstif(1)), 0)+1
1745	! count = end index - start index + 1
1746	nctcount = MIN((firstif(1)-1)+nbtime_perfile(if),end_globtime)-MAX(firstif(1),start_globtime)+1
1747	!
1748	!
1749	! Read time over the indexes computed above in order to position the slab correctly
1750	!
1751	WRITE(,) ">> From file ", if," reading from position ", nctstart, "up to ", (nctstart-1)+nctcount
1752	!
1753	DO it =1,nbtax
1754	tstart(1) = nctstart
1755	tcount(1) = nctcount
1756	iret = NF90_GET_VAR(force_id(if), time_id(if,it), time_tmp(inslabpos:inslabpos+nctcount-1,it), tstart, tcount)
1757	IF (iret /= NF90_NOERR) THEN
1758	WRITE(,) TRIM(NF90_STRERROR(iret))
1759	WRITE(,) "Working on file ", IF," starting at ",tstart(1)," counting ",tcount(1)
1760	WRITE(,) "The data was to be written in to section ", inslabpos,":",inslabpos+nctcount-1," of time_tmp"
1761	CALL ipslerr (3,'forcing_readslab_root',"Could not read the time for the current interval."," "," ")
1762	ENDIF
1763	time_tmp(inslabpos:inslabpos+nctcount-1,it) = date0_file(if,it) + &
1764	time_tmp(inslabpos:inslabpos+nctcount-1,it)*convtosec(if)/one_day
1765	ENDDO
1766	!
1767	! 2.0 Find and read variables.
1768	!
1769	! 2.1 Deal with air temperature and humidity as the fist and basic case
1770	!
1771	!
1772	!
1773	IF ( forcing_tstep_ave >= one_day/3.0) THEN
1774	CALL forcing_varforslab(if, "Tairmax", nctstart, nctcount, inslabpos, tairmax, cellmethod)
1775	CALL forcing_varforslab(if, "Tairmin", nctstart, nctcount, inslabpos, tairmin, cellmethod)
1776	ENDIF
1777
1778
1779	CALL forcing_varforslab(if, "Tair", nctstart, nctcount, inslabpos, tair, cellmethod)
1780	CALL forcing_attributetimeaxe(cellmethod, timeid_tair)
1781	!
1782	CALL forcing_varforslab(if, "Qair", nctstart, nctcount, inslabpos, qair, cellmethod)
1783	CALL forcing_attributetimeaxe(cellmethod, timeid_qair)
1784	!
1785	! 2.2 Deal with rainfall and snowfall.
1786	!
1787	CALL forcing_varforslab(if, "Rainf", nctstart, nctcount, inslabpos, rainf, cellmethod)
1788	CALL forcing_attributetimeaxe(cellmethod, timeid_precip)
1789	!
1790	CALL forcing_varforslab(if, "Snowf", nctstart, nctcount, inslabpos, snowf, cellmethod)
1791	CALL forcing_attributetimeaxe(cellmethod, timeid_tmp)
1792	IF ( timeid_precip .NE. timeid_tmp) THEN
1793	CALL ipslerr(3, 'forcing_readslab_root','Rainf and Snwof have different time axes.', &
1794	& 'Please check the forcing file to ensure both variable have the same cell_method.','')
1795	ENDIF
1796	!
1797	!
1798	! 2.3 Deal with downward shorwave and longwave radiation
1799	! The SW radiation can have 2 names SWdown_aerosol or SWdown. The first one is
1800	! given priority
1801	!
1802	CALL forcing_varforslab(if, "SWdown", nctstart, nctcount, inslabpos, swdown, cellmethod)
1803	CALL forcing_attributetimeaxe(cellmethod, timeid_swdown)
1804	!
1805	CALL forcing_varforslab(if, "LWdown", nctstart, nctcount, inslabpos, lwdown, cellmethod)
1806	CALL forcing_attributetimeaxe(cellmethod, timeid_lwdown)
1807	!
1808	!
1809	! 2.4 Deal with wind and pressure
1810	!
1811	CALL forcing_varforslab(if, "PSurf", nctstart, nctcount, inslabpos, ps, cellmethod)
1812	CALL forcing_attributetimeaxe(cellmethod, timeid_ps)
1813	!
1814	CALL forcing_varforslab(if, "Wind_E", nctstart, nctcount, inslabpos, u, cellmethod)
1815	CALL forcing_attributetimeaxe(cellmethod, timeid_u)
1816	!
1817	CALL forcing_varforslab(if, "Wind_N", nctstart, nctcount, inslabpos, v, cellmethod)
1818	CALL forcing_attributetimeaxe(cellmethod, timeid_v)
1819	!
1820	! Verify on Tair that we have a credible field.
1821	!
1822	IF ( MINVAL(tair(:,inslabpos:inslabpos+nctcount-1)) < 100.0 .OR. &
1823	& MAXVAL(tair(:,inslabpos:inslabpos+nctcount-1)) > 500.0 ) THEN
1824	WRITE(,) "ERROR on range of Tair : ", MINVAL(tair(:,inslabpos:inslabpos+nctcount-1)), &
1825	& MAXVAL(tair(:,inslabpos:inslabpos+nctcount-1))
1826	CALL ipslerr(3, 'forcing_readslab_root','The air temperature is not in a credible range.', &
1827	& 'Please verify your forcing file.','Are variables for all points to be simulated ?')
1828	ENDIF
1829	!
1830	! Do the height of the variables T&Q and U&V
1831	!
1832	! First the T&Q level
1833	!
1834	IF ( zheight ) THEN
1835	ztq(:,:) = zlev_fixed
1836	ELSE IF ( zsigma .OR. zhybrid ) THEN
1837	DO it=inslabpos,inslabpos+nctcount-1
1838	rau(:,it) = ps(:,it)/(cte_molr*tair(:,it))
1839	ztq(:,it) = (ps(:,it) - (zhybrid_a + zhybrid_bps(:,it)))/(rau(:,it) cte_grav)
1840	ENDDO
1841	ELSE IF ( zlevels ) THEN
1842	CALL forcing_varforslab(IF, "Levels", nctstart, nctcount, inslabpos, ztq, cellmethod)
1843	ELSE
1844	CALL ipslerr(3, 'forcing_readslab_root','No case for the vertical levels was specified.', &
1845	& 'We cannot determine the height for T and Q.','')
1846	ENDIF
1847	!
1848	! Now the U&V level
1849	!
1850	IF ( zsamelev_uv ) THEN
1851	zuv(:,:) = ztq(:,:)
1852	ELSE
1853	IF ( zheight ) THEN
1854	zuv(:,:) = zlevuv_fixed
1855	ELSE IF ( zsigma .OR. zhybrid ) THEN
1856	DO it=inslabpos,inslabpos+nctcount-1
1857	rau(:,it) = ps(:,it)/(cte_molr*tair(:,it))
1858	zuv(:,it) = (ps(:,it) - (zhybriduv_a + zhybriduv_bps(:,it)))/(rau(:,it) cte_grav)
1859	ENDDO
1860	ELSE IF ( zlevels ) THEN
1861	CALL forcing_varforslab(IF, "Levels_uv", nctstart, nctcount, inslabpos, zuv, cellmethod)
1862	ELSE
1863	CALL ipslerr(3, 'forcing_readslab_root','No case for the vertical levels was specified.', &
1864	& 'We cannot determine the height for U and V.','stop readdim2')
1865	ENDIF
1866	ENDIF
1867
1868	inslabpos = inslabpos+nctcount
1869
1870	ENDDO
1871	!
1872	! Use the read time of the slab to place it in the global variables. We consider
1873	! that we can do that on the first axis.
1874	!
1875	imin = MINLOC(ABS(time(:,1)-time_tmp(1,1)))
1876	position_slab(1) = imin(1)
1877	imax = MINLOC(ABS(time(:,1)-time_tmp(slab_size,1)))
1878	position_slab(2) = imax(1)
1879	!
1880	!
1881	IF ( position_slab(2)-position_slab(1) .GT. slab_size ) THEN
1882	DO it =1,nbtax
1883	WRITE(,) "Checking time_tmp on idex : ", it
1884	WRITE(,) "Time_tmp start and end : ",time_tmp(1,it), time_tmp(slab_size,it)
1885	imin = MINLOC(ABS(time(:,1)-time_tmp(1,it)))
1886	imax = MINLOC(ABS(time(:,1)-time_tmp(slab_size,it)))
1887	WRITE(,) "Interval read : ", imax(1)-imin(1)+1
1888	ENDDO
1889	CALL ipslerr (3,'forcing_readslab_root',"The time slab read does not fit the number of variables read.",&
1890	& "Could there be an error in the time axis ?"," ")
1891	ENDIF
1892	!
1893	! Transfer the global time axis into the time variables approriate for each variable. This way
1894	! the time axis for each variable will be centered on the interval of validity. This is an essential assumption
1895	! the time interpolation to be done later.
1896	!
1897	WRITE(,) "We have found the following axes : ", time_axename(:)
1898	WRITE(,) "For Tair we are using time axis '",TRIM(time_axename(timeid_tair)),&
1899	& "' with cell method ",TRIM(time_cellmethod(timeid_tair)),"."
1900	t_tair(1:slab_size) = time(position_slab(1):position_slab(2), timeid_tair)
1901	tbnd_tair(1:slab_size,:) = time_bounds(position_slab(1):position_slab(2),timeid_tair,:)
1902	!
1903	WRITE(,) "For Qair we are using time axis '",TRIM(time_axename(timeid_qair)),&
1904	& "' with cell method ",TRIM(time_cellmethod(timeid_qair)),"."
1905	t_qair(1:slab_size) = time(position_slab(1):position_slab(2), timeid_qair)
1906	tbnd_qair(1:slab_size,:) = time_bounds(position_slab(1):position_slab(2),timeid_qair,:)
1907	!
1908	WRITE(,) "For Rainf and Snowf we are using time axis '",TRIM(time_axename(timeid_precip)),&
1909	& "' with cell method ",TRIM(time_cellmethod(timeid_precip)),"."
1910	t_prec(1:slab_size) = time(position_slab(1):position_slab(2), timeid_precip)
1911	tbnd_prec(1:slab_size,:) = time_bounds(position_slab(1):position_slab(2),timeid_precip,:)
1912	prectime(1:slab_size) = preciptime(position_slab(1):position_slab(2))
1913	!
1914	WRITE(,) "For SWdown we are using time axis '",TRIM(time_axename(timeid_swdown)),&
1915	& "' with cell method ",TRIM(time_cellmethod(timeid_swdown)),"."
1916	t_swdown(1:slab_size) = time(position_slab(1):position_slab(2), timeid_swdown)
1917	tbnd_swdown(1:slab_size,:) = time_bounds(position_slab(1):position_slab(2),timeid_swdown,:)
1918	!
1919	WRITE(,) "For LWdown we are using time axis '",TRIM(time_axename(timeid_lwdown)),&
1920	& "' with cell method ",TRIM(time_cellmethod(timeid_lwdown)),"."
1921	t_lwdown(1:slab_size) = time(position_slab(1):position_slab(2), timeid_lwdown)
1922	tbnd_lwdown(1:slab_size,:) = time_bounds(position_slab(1):position_slab(2),timeid_lwdown,:)
1923	!
1924	WRITE(,) "For Wind_E we are using time axis '",TRIM(time_axename(timeid_u)),&
1925	& "' with cell method ",TRIM(time_cellmethod(timeid_u)),"."
1926	t_u(1:slab_size) = time(position_slab(1):position_slab(2), timeid_u)
1927	tbnd_u(1:slab_size,:) = time_bounds(position_slab(1):position_slab(2),timeid_u,:)
1928	!
1929	WRITE(,) "For Wind_N we are using time axis '",TRIM(time_axename(timeid_v)),&
1930	& "' with cell method ",TRIM(time_cellmethod(timeid_v)),"."
1931	t_v(1:slab_size) = time(position_slab(1):position_slab(2), timeid_v)
1932	tbnd_v(1:slab_size,:) = time_bounds(position_slab(1):position_slab(2),timeid_v,:)
1933	!
1934	WRITE(,) "For PSurf we are using time axis '",TRIM(time_axename(timeid_ps)),&
1935	& "' with cell method ",TRIM(time_cellmethod(timeid_ps)),"."
1936	t_ps(1:slab_size) = time(position_slab(1):position_slab(2), timeid_ps)
1937	tbnd_ps(1:slab_size,:) = time_bounds(position_slab(1):position_slab(2),timeid_ps,:)
1938	!
1939	ELSE
1940	CALL ipslerr (2,'forcing_readslab_root',"We have reached the end of the slabs we can read.",&
1941	& "The calling program needs to manage this situation"," ")
1942	ENDIF
1943	!
1944	! Have we read to the end of the files ?
1945	!
1946	IF ( current_offset+slab_size >= nb_forcing_steps ) THEN
1947	end_of_file = .TRUE.
1948	ELSE
1949	end_of_file = .FALSE.
1950	ENDIF
1951	!
1952	IF ( ALLOCATED(rau) ) DEALLOCATE(rau)
1953	IF ( ALLOCATED(time_tmp) ) DEALLOCATE(time_tmp)
1954	!
1955	END SUBROUTINE forcing_readslab_root
1956
1957	!! =============================================================================================================================
1958	!! SUBROUTINE: forcing_reindex3d
1959	!!
1960	!>\BRIEF
1961	!!
1962	!! DESCRIPTION:
1963	!!
1964	!! \n
1965	!_ ==============================================================================================================================
1966	SUBROUTINE forcing_reindex3d(nbi, nbj, tin, slab_in, nbout, tout, slab_out, tstart, reindex)
1967	!
1968	! ARGUMENTS
1969	!
1970	INTEGER(i_std), INTENT(in) :: nbi, nbj, tin, nbout, tout
1971	REAL(r_std), INTENT(in) :: slab_in(nbi,nbj,tin)
1972	REAL(r_std), INTENT(out) :: slab_out(nbout,tout)
1973	INTEGER(i_std), INTENT(in) :: tstart
1974	INTEGER(i_std), INTENT(in) :: reindex(nbout,2)
1975	!
1976	! LOCAL
1977	!
1978	INTEGER(i_std) :: is, in
1979	!
1980	DO is=1,tin
1981	DO in=1,nbout
1982	slab_out(in,tstart+(is-1)) = slab_in(reindex(in,1),reindex(in,2),is)
1983	ENDDO
1984	ENDDO
1985	!
1986	END SUBROUTINE forcing_reindex3d
1987
1988	!! =============================================================================================================================
1989	!! SUBROUTINE: forcing_reindex2d
1990	!!
1991	!>\BRIEF
1992	!!
1993	!! DESCRIPTION:
1994	!!
1995	!! \n
1996	!_ ==============================================================================================================================
1997	SUBROUTINE forcing_reindex2d(nbi, nbj, slab_in, nbout, slab_out, reindex)
1998	!
1999	! ARGUMENTS
2000	!
2001	INTEGER(i_std), INTENT(in) :: nbi, nbj, nbout
2002	REAL(r_std), INTENT(in) :: slab_in(nbi,nbj)
2003	REAL(r_std), INTENT(out) :: slab_out(nbout)
2004	INTEGER(i_std), INTENT(in) :: reindex(nbout,2)
2005	!
2006	! LOCAL
2007	!
2008	INTEGER(i_std) :: in
2009	!
2010	DO in=1,nbout
2011	slab_out(in) = slab_in(reindex(in,1),reindex(in,2))
2012	ENDDO
2013	!
2014	END SUBROUTINE forcing_reindex2d
2015	!!
2016	!! =============================================================================================================================
2017	!! SUBROUTINE: forcing_reindex2dt
2018	!!
2019	!>\BRIEF
2020	!!
2021	!! DESCRIPTION:
2022	!!
2023	!! \n
2024	!_ ==============================================================================================================================
2025
2026	SUBROUTINE forcing_reindex2dt(nbin, tin, slab_in, nbout, tout, slab_out, tstart, reindex)
2027	!
2028	! ARGUMENTS
2029	!
2030	INTEGER(i_std), INTENT(in) :: nbin, tin, nbout, tout
2031	REAL(r_std), INTENT(in) :: slab_in(nbin,tin)
2032	REAL(r_std), INTENT(out) :: slab_out(nbout,tout)
2033	INTEGER(i_std), INTENT(in) :: tstart
2034	INTEGER(i_std), INTENT(in) :: reindex(nbout)
2035	!
2036	! LOCAL
2037	!
2038	INTEGER(i_std) :: is, in
2039	!
2040	DO is=1,tin
2041	DO in=1,nbout
2042	slab_out(in,tstart+(is-1)) = slab_in(reindex(in),is)
2043	ENDDO
2044	ENDDO
2045	!
2046	END SUBROUTINE forcing_reindex2dt
2047
2048	!! =============================================================================================================================
2049	!! SUBROUTINE: forcing_reindex1d
2050	!!
2051	!>\BRIEF
2052	!!
2053	!! DESCRIPTION:
2054	!!
2055	!! \n
2056	!_ ==============================================================================================================================
2057
2058	SUBROUTINE forcing_reindex1d(nbin, slab_in, nbout, slab_out, reindex)
2059	!
2060	! ARGUMENTS
2061	!
2062	INTEGER(i_std), INTENT(in) :: nbin, nbout
2063	REAL(r_std), INTENT(in) :: slab_in(nbin)
2064	REAL(r_std), INTENT(out) :: slab_out(nbout)
2065	INTEGER(i_std), INTENT(in) :: reindex(nbout)
2066	!
2067	! LOCAL
2068	!
2069	INTEGER(i_std) :: is
2070	!
2071	DO is=1,nbout
2072	slab_out(is) = slab_in(reindex(is))
2073	ENDDO
2074	!
2075	END SUBROUTINE forcing_reindex1d
2076	!!
2077	!! =============================================================================================================================
2078	!! SUBROUTINE: forcing_reindex2to1
2079	!!
2080	!>\BRIEF
2081	!!
2082	!! DESCRIPTION:
2083	!!
2084	!! \n
2085	!_ ==============================================================================================================================
2086
2087	SUBROUTINE forcing_reindex2to1(nbi, nbj, slab_in, nbout, slab_out, reindex)
2088	!
2089	! ARGUMENTS
2090	!
2091	INTEGER(i_std), INTENT(in) :: nbi, nbj, nbout
2092	REAL(r_std), INTENT(in) :: slab_in(nbi,nbj)
2093	REAL(r_std), INTENT(out) :: slab_out(nbout)
2094	INTEGER(i_std), INTENT(in) :: reindex(nbout)
2095	!
2096	! LOCAL
2097	!
2098	INTEGER(i_std) :: i, j, is
2099	!
2100	DO is=1,nbout
2101	j = INT((reindex(is)-1)/nbi)+1
2102	i = (reindex(is)-(j-1)*nbi)
2103	slab_out(is) = slab_in(i,j)
2104	ENDDO
2105	!
2106	END SUBROUTINE forcing_reindex2to1
2107
2108	!! =============================================================================================================================
2109	!! SUBROUTINE: forcing_reindex1to2
2110	!!
2111	!>\BRIEF
2112	!!
2113	!! DESCRIPTION:
2114	!!
2115	!! \n
2116	!_ ==============================================================================================================================
2117
2118	SUBROUTINE forcing_reindex1to2(nbin, slab_in, nbi, nbj, slab_out, reindex)
2119	!
2120	! ARGUMENTS
2121	!
2122	INTEGER(i_std), INTENT(in) :: nbin, nbi, nbj
2123	REAL(r_std), INTENT(in) :: slab_in(nbin)
2124	REAL(r_std), INTENT(out) :: slab_out(nbi, nbj)
2125	INTEGER(i_std), INTENT(in) :: reindex(nbin)
2126	!
2127	! LOCAL
2128	!
2129	INTEGER(i_std) :: i, j, is
2130	!
2131	DO is=1,nbin
2132	j = INT((reindex(is)-1)/nbi)+1
2133	i = (reindex(is)-(j-1)*nbi)
2134	slab_out(i,j) = slab_in(is)
2135	ENDDO
2136	!
2137	END SUBROUTINE forcing_reindex1to2
2138
2139	!! =============================================================================================================================
2140	!! SUBROUTINE: forcing_close
2141	!!
2142	!>\BRIEF Close all forcing files
2143	!!
2144	!! DESCRIPTION:
2145	!!
2146	!! \n
2147	!_ ==============================================================================================================================
2148	SUBROUTINE forcing_close()
2149
2150	INTEGER(i_std) :: ierr, if
2151
2152	DO if=1,nb_forcefile
2153	ierr = NF90_CLOSE(force_id(if))
2154	ENDDO
2155
2156	END SUBROUTINE forcing_close
2157
2158	!! =============================================================================================================================
2159	!! SUBROUTINE: forcing_printdate
2160	!!
2161	!>\BRIEF Print the date in a human readable format.
2162	!!
2163	!! DESCRIPTION:
2164	!!
2165	!! \n
2166	!_ ==============================================================================================================================
2167
2168	SUBROUTINE forcing_printdate(julian_day, message, wunit)
2169	!
2170	REAL(r_std), INTENT(in) :: julian_day
2171	CHARACTER(len=*), INTENT(in) :: message
2172	INTEGER, INTENT(in), OPTIONAL :: wunit
2173	!
2174	!
2175	!
2176	INTEGER(i_std) :: year, month, day, hours, minutes, seci
2177	REAL(r_std) :: sec
2178	!
2179	CALL ju2ymds (julian_day, year, month, day, sec)
2180	hours = INT(sec/3600)
2181	sec = sec - 3600 * hours
2182	minutes = INT(sec / 60)
2183	sec = sec - 60 * minutes
2184	seci = INT(sec)
2185	!
2186	IF (PRESENT(wunit)) THEN
2187	WRITE(wunit,'(I4.4,"-",I2.2,"-",I2.2," ",I2.2,":",I2.2,":",I2.2," > ", A60)') &
2188	& year, month, day, hours, minutes, seci, message
2189	ELSE
2190	WRITE(*,'(I4.4,"-",I2.2,"-",I2.2," ",I2.2,":",I2.2,":",I2.2," > ", A60)') &
2191	& year, month, day, hours, minutes, seci, message
2192	ENDIF
2193	!
2194	END SUBROUTINE forcing_printdate
2195
2196	!! =============================================================================================================================
2197	!! SUBROUTINE: forcing_printpoint_forgrid
2198	!!
2199	!>\BRIEF Together with the date print some sample values. Useful for checking the forcing.
2200	!!
2201	!! DESCRIPTION:
2202	!!
2203	!! \n
2204	!_ ==============================================================================================================================
2205
2206	SUBROUTINE forcing_printpoint_forgrid(julian_day, lon_pt, lat_pt, var, message)
2207	!
2208	REAL(r_std), INTENT(in) :: julian_day
2209	REAL(r_std), INTENT(in) :: lon_pt, lat_pt
2210	REAL(r_std), INTENT(in) :: var(:)
2211	CHARACTER(len=*), INTENT(in) :: message
2212	!
2213	!
2214	!
2215	INTEGER(i_std) :: year, month, day, hours, minutes, seci
2216	REAL(r_std) :: sec
2217	INTEGER(i_std) :: lon_ind, lat_ind, ind
2218	INTEGER(i_std), DIMENSION(1) :: i,j,k
2219	!
2220	! Check if there is anything to be done
2221	!
2222	IF ( MAX(lon_pt, lat_pt) > 360.0 ) THEN
2223	RETURN
2224	ENDIF
2225	!
2226	! Convert time first
2227	!
2228	CALL ju2ymds (julian_day, year, month, day, sec)
2229	hours = INT(sec/3600)
2230	sec = sec - 3600 * hours
2231	minutes = INT(sec / 60)
2232	sec = sec - 60 * minutes
2233	seci = INT(sec)
2234	!
2235	! Get the local to be analysed
2236	!
2237	i = MINLOC(ABS(lon_loc(:,1)-lon_pt))
2238	j = MINLOC(ABS(lat_loc(1,:)-lat_pt))
2239	ind = (j(1)-1)*iim_loc+i(1)
2240	k = MINLOC(ABS(lindex_loc(:)-ind))
2241	!
2242	WRITE(*,"(I2.2,':',I2.2,':',I2.2,' Loc : ', F5.1,',', F5.1,'(i=',I6,') Value = ',F12.4,A40)") &
2243	& hours, minutes, seci, lon_loc(i(1),1), lat_loc(1,j(1)), k(1), var(k(1)), message
2244
2245	END SUBROUTINE forcing_printpoint_forgrid
2246	!! =============================================================================================================================
2247	!! SUBROUTINE: forcing_printpoint_forgrid2d
2248	!!
2249	!>\BRIEF Together with the date print some sample values. Useful for checking the forcing.
2250	!!
2251	!! DESCRIPTION:
2252	!!
2253	!! \n
2254	!_ ==============================================================================================================================
2255
2256	SUBROUTINE forcing_printpoint_forgrid2d(julian_day, lon_pt, lat_pt, var, message)
2257	!
2258	REAL(r_std), INTENT(in) :: julian_day
2259	REAL(r_std), INTENT(in) :: lon_pt, lat_pt
2260	REAL(r_std), INTENT(in) :: var(:,:)
2261	CHARACTER(len=*), INTENT(in) :: message
2262	!
2263	!
2264	!
2265	INTEGER(i_std) :: year, month, day, hours, minutes, seci
2266	REAL(r_std) :: sec
2267	INTEGER(i_std) :: lon_ind, lat_ind
2268	INTEGER(i_std), DIMENSION(1) :: i,j
2269	!
2270	! Check if there is anything to be done
2271	!
2272	IF ( MAX(lon_pt, lat_pt) > 360.0 ) THEN
2273	RETURN
2274	ENDIF
2275	!
2276	! Convert time first
2277	!
2278	CALL ju2ymds (julian_day, year, month, day, sec)
2279	hours = INT(sec/3600)
2280	sec = sec - 3600 * hours
2281	minutes = INT(sec / 60)
2282	sec = sec - 60 * minutes
2283	seci = INT(sec)
2284	!
2285	! Get the local to be analysed
2286	!
2287	i = MINLOC(ABS(lon_loc(:,1)-lon_pt))
2288	j = MINLOC(ABS(lat_loc(1,:)-lat_pt))
2289	!
2290	WRITE(*,"(I2.2,':',I2.2,':',I2.2,' Loc : ', F5.1,',', F5.1,'(i=',I6,') Value = ',F12.4,A40)") &
2291	& hours, minutes, seci, lon_loc(i(1),1), lat_loc(1,j(1)), i(1), j(1), var(i(1),j(1)), message
2292
2293	END SUBROUTINE forcing_printpoint_forgrid2d
2294
2295	!! =============================================================================================================================
2296	!! SUBROUTINE: forcing_printpoint_gen
2297	!!
2298	!>\BRIEF Together with the date print some sample values. Useful for checking the forcing.
2299	!!
2300	!! DESCRIPTION:
2301	!!
2302	!! \n
2303	!_ ==============================================================================================================================
2304
2305	SUBROUTINE forcing_printpoint_gen(julian_day, lon_pt, lat_pt, nbind, lalo_in, var, message, ktest)
2306	!
2307	REAL(r_std), INTENT(in) :: julian_day
2308	REAL(r_std), INTENT(in) :: lon_pt, lat_pt
2309	INTEGER(i_std), INTENT(in) :: nbind
2310	REAL(r_std), INTENT(in) :: lalo_in(:,:)
2311	REAL(r_std), INTENT(in) :: var(:)
2312	CHARACTER(len=*), INTENT(in) :: message
2313	INTEGER(i_std), OPTIONAL, INTENT(out) :: ktest
2314	!
2315	!
2316	!
2317	INTEGER(i_std) :: year, month, day, hours, minutes, seci
2318	REAL(r_std) :: sec, mindist
2319	REAL(r_std), ALLOCATABLE, DIMENSION(:) :: dist, refdist
2320	INTEGER(i_std) :: lon_ind, lat_ind, ind
2321	INTEGER(i_std) :: i, imin(1)
2322	REAL(r_std), PARAMETER :: mincos = 0.0001
2323	REAL(r_std), PARAMETER :: pi = 3.141592653589793238
2324	REAL(r_std), PARAMETER :: R_Earth = 6378000.
2325	!
2326	! Check if there is anything to be done
2327	!
2328	IF ( MAX(lon_pt, lat_pt) > 360.0 ) THEN
2329	IF ( PRESENT(ktest) ) ktest = -1
2330	RETURN
2331	ENDIF
2332	!
2333	! Allocate memory
2334	!
2335	ALLOCATE(dist(nbind))
2336	ALLOCATE(refdist(nbind))
2337	!
2338	! Convert time first
2339	!
2340	CALL ju2ymds (julian_day, year, month, day, sec)
2341	hours = INT(sec/3600)
2342	sec = sec - 3600 * hours
2343	minutes = INT(sec / 60)
2344	sec = sec - 60 * minutes
2345	seci = INT(sec)
2346	!
2347	! Get the location to be analysed
2348	!
2349	DO i=1,nbind
2350	dist(i) = acos( sin(lat_ptpi/180)sin(lalo_in(i,1)*pi/180) + &
2351	& cos(lat_ptpi/180)cos(lalo_in(i,1)pi/180)&
2352	& cos((lalo_in(i,2)-lon_pt)pi/180) ) R_Earth
2353	ENDDO
2354	!
2355	! Look for the next grid point closest to the one with the smalest distance.
2356	!
2357	imin = MINLOC(dist)
2358	DO i=1,nbind
2359	refdist(i) = acos( sin(lalo_in(imin(1),1)pi/180)sin(lalo_in(i,1)*pi/180) + &
2360	& cos(lalo_in(imin(1),1)pi/180)cos(lalo_in(i,1)pi/180) &
2361	& cos((lalo_in(i,2)-lalo_in(imin(1),2))pi/180) ) R_Earth
2362	ENDDO
2363	refdist(imin(1)) = MAXVAL(refdist)
2364	mindist = MINVAL(refdist)
2365	!
2366	! Are we closer than the closest points ?
2367	!
2368	IF ( PRESENT(ktest) ) ktest = -1
2369	IF ( dist(imin(1)) <= mindist ) THEN
2370	!
2371	WRITE(*,"(I2.2,':',I2.2,':',I2.2,' Loc : ', F6.1,',', F6.1,'(i=',I6,') Value = ',F12.4,A38)") &
2372	& hours, minutes, seci, lalo_in(imin(1),2), lalo_in(imin(1),1), imin(1), var(imin(1)), message
2373	!
2374	IF ( PRESENT(ktest) ) ktest = imin(1)
2375	ENDIF
2376	!
2377	END SUBROUTINE forcing_printpoint_gen
2378
2379	!! =============================================================================================================================
2380	!! SUBROUTINE: forcing_getglogrid
2381	!!
2382	!>\BRIEF Routine to read the full grid of the forcing file.
2383	!!
2384	!! DESCRIPTION: The data is stored in the saved variables of the module and serve all other routines.
2385	!!
2386	!! \n
2387	!_ ==============================================================================================================================
2388
2389	SUBROUTINE forcing_getglogrid (nbfiles, filename, iim_tmp, jjm_tmp, nbpoint_tmp, closefile, landonly_arg)
2390	!
2391	! This routine reads the global grid information from the forcing file and generates all the
2392	! information needed for this grid.
2393	!
2394	! ARGUMENTS
2395	!
2396	INTEGER(i_std), INTENT(in) :: nbfiles
2397	CHARACTER(LEN=*), INTENT(in) :: filename(:)
2398	INTEGER(i_std), INTENT(out) :: iim_tmp, jjm_tmp, nbpoint_tmp
2399	LOGICAL, INTENT(in) :: closefile
2400	LOGICAL, OPTIONAL, INTENT(in) :: landonly_arg
2401	!
2402	! LOCAL
2403	!
2404	INTEGER(i_std) :: iret, iv, if, lll
2405	CHARACTER(LEN=20) :: dimname, varname
2406	CHARACTER(LEN=60) :: lon_units, lat_units, units
2407	INTEGER(i_std), ALLOCATABLE, DIMENSION(:) :: dimids, londim_ids, latdim_ids
2408	INTEGER(i_std) :: lon_id, lat_id, land_id, lon_nbdims, lat_nbdims, land_nbdims
2409	INTEGER(i_std) :: lonvar_id, latvar_id, landvar_id
2410	LOGICAL :: dump_mask
2411	INTEGER(i_std) :: ik, i, j, iff, ndimsvar
2412	! Read a test variabe
2413	CHARACTER(len=8) :: testvarname='Tair'
2414	INTEGER(i_std) :: testvar_id, contfrac_id
2415	REAL(r_std) :: testvar_missing, contfrac_missing
2416	REAL(r_std), ALLOCATABLE, DIMENSION(:) :: testvar
2417	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: testvar2d, contfrac2d
2418	!
2419	! 0.0 Check variables are allocated
2420	!
2421	IF ( .NOT. ALLOCATED(force_id)) ALLOCATE(force_id(nbfiles))
2422	IF ( .NOT. ALLOCATED(id_unlim)) ALLOCATE(id_unlim(nbfiles))
2423	IF ( .NOT. ALLOCATED(nb_atts)) ALLOCATE(nb_atts(nbfiles))
2424	IF ( .NOT. ALLOCATED(ndims)) ALLOCATE(ndims(nbfiles))
2425	IF ( .NOT. ALLOCATED(nvars)) ALLOCATE( nvars(nbfiles))
2426	!
2427	! 0.1 Are we one the root proc ?
2428	!
2429	IF ( .NOT. is_root_prc ) THEN
2430	CALL ipslerr (3,'forcing_getglogrid'," This routine can only be called on the root processor.", " ", " ")
2431	ENDIF
2432	!
2433	! The default behaviour is to provide only land points to the calling program.
2434	! But for forcing ocean model there is also the option to pass on land and ocean values.
2435	! When the grid is initialized landonly_tmp=.FALSE. has to be set to obtian this behaviour.
2436	!
2437	IF ( PRESENT(landonly_arg) ) THEN
2438	landonly=landonly_arg
2439	ELSE
2440	landonly=.TRUE.
2441	ENDIF
2442	!
2443	! 1.0 Open the netCDF file and get basic dimensions
2444	!
2445	DO iff=1,nbfiles
2446	!
2447	iret = NF90_OPEN(filename(iff), NF90_NOWRITE, force_id(iff))
2448	IF (iret /= NF90_NOERR) THEN
2449	CALL ipslerr (3,'forcing_getglogrid',"Error opening the forcing file :", filename(iff), " ")
2450	ENDIF
2451	!
2452	iret = NF90_INQUIRE (force_id(iff), nDimensions=ndims(iff), nVariables=nvars(iff), &
2453	nAttributes=nb_atts(iff), unlimitedDimId=id_unlim(iff))
2454	!
2455	!
2456	! 2.0 Read the dimensions found in the forcing file. Only deal with the spatial dimension as
2457	! time is an unlimited dimension.
2458	!
2459	DO iv=1,ndims(iff)
2460	!
2461	iret = NF90_INQUIRE_DIMENSION (force_id(iff), iv, name=dimname, len=lll)
2462	IF (iret /= NF90_NOERR) THEN
2463	CALL ipslerr (3,'forcing_getglogrid',"Could not get size of dimensions in file : ", filename(iff), " ")
2464	ENDIF
2465	!
2466	SELECT CASE(lowercase(dimname))
2467	!
2468	CASE("west_east")
2469	CALL forcing_checkdim(iff, filename, iim_glo, lon_id, lll, iv)
2470	CASE("south_north")
2471	CALL forcing_checkdim(iff, filename, jjm_glo, lat_id, lll, iv)
2472	CASE("longitude")
2473	CALL forcing_checkdim(iff, filename, iim_glo, lon_id, lll, iv)
2474	CASE("latitude")
2475	CALL forcing_checkdim(iff, filename, jjm_glo, lat_id, lll, iv)
2476	CASE("lon")
2477	CALL forcing_checkdim(iff, filename, iim_glo, lon_id, lll, iv)
2478	CASE("lat")
2479	CALL forcing_checkdim(iff, filename, jjm_glo, lat_id, lll, iv)
2480	CASE("nav_lon")
2481	CALL forcing_checkdim(iff, filename, iim_glo, lon_id, lll, iv)
2482	CASE("nav_lat")
2483	CALL forcing_checkdim(iff, filename, jjm_glo, lat_id, lll, iv)
2484	CASE("x")
2485	CALL forcing_checkdim(iff, filename, iim_glo, lon_id, lll, iv)
2486	CASE("y")
2487	CALL forcing_checkdim(iff, filename, jjm_glo, lat_id, lll, iv)
2488	CASE("land")
2489	CALL forcing_checkdim(iff, filename, nbland_glo, land_id, lll, iv)
2490	END SELECT
2491	!
2492	ENDDO
2493	IF ( iim_glo == 0 .AND. jjm_glo == 0 ) THEN
2494	CALL ipslerr (3,'forcing_getglogrid',"Did not recognize any dimensions in : ", filename(iff), " ")
2495	ENDIF
2496	ENDDO
2497
2498	!
2499	! 3.0 Read the spatial coordinate variables found in the first file.
2500	!
2501	ALLOCATE(dimids(NF90_MAX_VAR_DIMS), londim_ids(NF90_MAX_VAR_DIMS), latdim_ids(NF90_MAX_VAR_DIMS))
2502	lonvar_id = -1
2503	latvar_id = -1
2504	landvar_id = -1
2505	testvar_id = -1
2506	contfrac_id = -1
2507	! Count the number of time axis we have
2508	nbtax = 0
2509	!
2510	DO iv=1,nvars(1)
2511	!
2512	iret = NF90_INQUIRE_VARIABLE(force_id(1), iv, name=varname, ndims=ndimsvar, dimids=dimids)
2513	iret = NF90_GET_ATT(force_id(1), iv, 'units', units)
2514	!
2515	! Check that we have the longitude
2516	!
2517	IF ( INDEX(lowercase(varname), 'lon') > 0 .AND. INDEX(lowercase(units), 'east') > 0 ) THEN
2518	lonvar_id = iv
2519	lon_units=units
2520	lon_nbdims = ndimsvar
2521	londim_ids = dimids
2522	ENDIF
2523	!
2524	! Check that we have the latitude
2525	!
2526	IF ( INDEX(lowercase(varname), 'lat') > 0 .AND. INDEX(lowercase(units), 'north') > 0) THEN
2527	latvar_id = iv
2528	lat_units=units
2529	lat_nbdims = ndimsvar
2530	latdim_ids = dimids
2531	ENDIF
2532	!
2533	! Check that we have the land re-indexing table
2534	!
2535	IF ( INDEX(lowercase(varname), 'land') > 0 ) THEN
2536	landvar_id = iv
2537	land_nbdims = ndimsvar
2538	latdim_ids = dimids
2539	ENDIF
2540	!
2541	! Check if we have the contfrac variable
2542	!
2543	IF ( INDEX(lowercase(varname), 'contfrac') > 0 ) THEN
2544	contfrac_id = iv
2545	iret = NF90_GET_ATT(force_id(1), iv, "missing_value", contfrac_missing)
2546	IF (iret /= NF90_NOERR) THEN
2547	! No missing_value found, try to read _FillValue instead
2548	iret = NF90_GET_ATT(force_id(1), iv, "_FillValue", contfrac_missing)
2549	IF (iret /= NF90_NOERR) THEN
2550	WRITE(,) TRIM(nf90_strerror(iret))
2551	WRITE(,) " >> No _FillValue or missing_value found for contfrac"
2552	contfrac_missing=0.0
2553	END IF
2554	ENDIF
2555	ENDIF
2556	!
2557	! Find the test variable
2558	!
2559	IF ( INDEX(lowercase(varname), TRIM(lowercase(testvarname))) > 0 .AND. &
2560	& LEN_TRIM(varname) == LEN_TRIM(testvarname)) THEN
2561	testvar_id = iv
2562	iret = NF90_GET_ATT(force_id(1), iv, "missing_value", testvar_missing)
2563	IF (iret /= NF90_NOERR) THEN
2564	! No missing_value found, try to read _FillValue instead
2565	iret = NF90_GET_ATT(force_id(1), iv, "_FillValue", testvar_missing)
2566	IF (iret /= NF90_NOERR) THEN
2567	WRITE(,) TRIM(nf90_strerror(iret))
2568	WRITE(,) " >> No _FillValue or missing_value found for variable=",varname
2569	testvar_missing=-1
2570	END IF
2571	ENDIF
2572	ENDIF
2573	!
2574	! If we come across time get the calendar and archive it.
2575	!
2576	IF ( INDEX(lowercase(units),'seconds since') > 0 .OR. &
2577	& INDEX(lowercase(units),'minutes since') > 0 .OR. &
2578	& INDEX(lowercase(units),'hours since') > 0 .OR. &
2579	& INDEX(lowercase(units),'days since') > 0) THEN
2580	!
2581	! Get calendar used for the time axis
2582	!
2583	iret = NF90_GET_ATT(force_id(1), iv, "calendar", calendar)
2584	IF (iret == NF90_NOERR) THEN
2585	WRITE(,) ">> Setting the calendar to ",calendar
2586	ELSE
2587	WRITE(,) ">> Keeping proleptic Gregorian calendar"
2588	calendar="proleptic_gregorian"
2589	ENDIF
2590	!
2591	nbtax = nbtax + 1
2592	!
2593	ENDIF
2594	ENDDO
2595	!
2596	! 4.0 Verification that we have found both coordinate variables and the land point index
2597	!
2598	IF ( ( lonvar_id < 0 ) .AND. ( INDEX(lowercase(lon_units), 'east') <= 0 ) ) THEN
2599	CALL ipslerr (3,'forcing_getglogrid',"Have not found a valid longitude. Units = ", lon_units, " ")
2600	ENDIF
2601	IF ( ( latvar_id < 0 ) .AND. ( INDEX(lowercase(lat_units), 'north') <= 0 ) ) THEN
2602	CALL ipslerr (3,'forcing_getglogrid',"Have not found a valid latitude. Units = : ", lat_units, " ")
2603	ENDIF
2604	IF ( landvar_id < 0 ) THEN
2605	CALL ipslerr (1,'forcing_getglogrid',"No reindexing table was found. ", &
2606	& "This forcing file is not compressed by gathering.", " ")
2607	ENDIF
2608	!
2609	! 5.0 Allocate the space for the global variables and read them.
2610	!
2611	IF ( .NOT. ALLOCATED(lon_glo)) ALLOCATE(lon_glo(iim_glo, jjm_glo))
2612	IF ( .NOT. ALLOCATED(lat_glo)) ALLOCATE(lat_glo(iim_glo, jjm_glo))
2613	!
2614	IF ( lon_nbdims == 2 .AND. lat_nbdims == 2 ) THEN
2615	iret = NF90_GET_VAR(force_id(1), lonvar_id, lon_glo)
2616	iret = NF90_GET_VAR(force_id(1), latvar_id, lat_glo)
2617	ELSE IF ( lon_nbdims == 1 .AND. lat_nbdims == 1 ) THEN
2618	DO iv=1,jjm_glo
2619	iret = NF90_GET_VAR(force_id(1), lonvar_id, lon_glo(:,iv))
2620	ENDDO
2621	DO iv=1,iim_glo
2622	iret = NF90_GET_VAR(force_id(1), latvar_id, lat_glo(iv,:))
2623	ENDDO
2624	ELSE
2625	WRITE(,) "Found dimensions for lon and lat :", lon_nbdims, lat_nbdims
2626	CALL ipslerr (3,'forcing_getglogrid',"Longitude and Latitude variables do not have the right dimensions.", " ", " ")
2627	ENDIF
2628	!
2629	! If we have a indexing variable then the data is compressed by gathering, else we have to construct it.
2630	!
2631	compressed = .FALSE.
2632	IF ( landvar_id > 0 ) THEN
2633	IF ( .NOT. ALLOCATED(lindex_glo)) ALLOCATE(lindex_glo(nbland_glo))
2634	iret = NF90_GET_VAR(force_id(1), landvar_id, lindex_glo)
2635	compressed = .TRUE.
2636	ENDIF
2637	!
2638	IF ( .NOT. ALLOCATED(mask_glo)) ALLOCATE(mask_glo(iim_glo, jjm_glo))
2639	!
2640	! Get the land/sea mask and contfrac based on a test variable, if contfract is not available. Else
2641	! we use the contfrac variable.
2642	!
2643	IF ( compressed ) THEN
2644	IF ( landonly ) THEN
2645	IF ( .NOT. ALLOCATED(contfrac_glo)) ALLOCATE(contfrac_glo(nbland_glo))
2646	IF ( .NOT. ALLOCATED(testvar)) ALLOCATE(testvar(nbland_glo))
2647	CALL forcing_contfrac1d(force_id(1), testvar_id, contfrac_id, testvar)
2648	nbpoint_glo = nbland_glo
2649	ELSE
2650	WRITE(,) "forcing_tools cannot provide data over ocean points as the"
2651	WRITE(,) "data in the file is compressed by gathering land points."
2652	WRITE(,) "Fatal error"
2653	CALL ipslerr (3,'forcing_getglogrid',"forcing_tools cannot provide data over ocean points as the", &
2654	& "data in the file is compressed by gathering land points.", " ")
2655	ENDIF
2656	ELSE
2657	IF ( .NOT. ALLOCATED(testvar2d)) ALLOCATE(testvar2d(iim_glo, jjm_glo))
2658	IF ( .NOT. ALLOCATED(contfrac2d)) ALLOCATE(contfrac2d(iim_glo, jjm_glo))
2659	CALL forcing_contfrac2d(force_id(1), testvar_id, contfrac_id, testvar2d, contfrac2d, &
2660	& testvar_missing, contfrac_missing, nbland_glo)
2661
2662
2663	!
2664	! We have found a variable on which we can count the number of land points. We can build
2665	! the indexing table and gather the information needed.
2666	!
2667	IF ( landonly ) THEN
2668	nbpoint_glo = nbland_glo
2669
2670
2671	IF ( .NOT. ALLOCATED(lindex_glo)) ALLOCATE(lindex_glo(nbpoint_glo))
2672	IF ( .NOT. ALLOCATED(contfrac_glo)) ALLOCATE(contfrac_glo(nbpoint_glo))
2673	IF ( .NOT. ALLOCATED(testvar)) ALLOCATE(testvar(nbpoint_glo))
2674
2675
2676	IF ( contfrac_id > 0 ) THEN
2677	CALL forcing_buildindex(contfrac2d, contfrac_missing, landonly, lindex_glo, contfrac_glo)
2678	CALL forcing_reindex(iim_glo, jjm_glo, testvar2d, nbland_glo, testvar, lindex_glo)
2679
2680
2681	ELSE
2682	CALL forcing_buildindex(testvar2d, testvar_missing, landonly, lindex_glo, testvar)
2683	contfrac_glo(:) = 1.0
2684
2685	ENDIF
2686	ELSE
2687	nbpoint_glo = iim_glo*jjm_glo
2688	IF ( .NOT. ALLOCATED(lindex_glo)) ALLOCATE(lindex_glo(nbpoint_glo))
2689	IF ( .NOT. ALLOCATED(contfrac_glo)) ALLOCATE(contfrac_glo(nbpoint_glo))
2690	IF ( .NOT. ALLOCATED(testvar)) ALLOCATE(testvar(nbpoint_glo))
2691	IF ( contfrac_id > 0 ) THEN
2692	CALL forcing_buildindex(contfrac2d, contfrac_missing, landonly, lindex_glo, contfrac_glo)
2693	CALL forcing_reindex(iim_glo, jjm_glo, testvar2d, nbland_glo, testvar, lindex_glo)
2694	ELSE
2695	CALL forcing_buildindex(testvar2d, testvar_missing, landonly, lindex_glo, testvar)
2696	contfrac_glo(:) = 1.0
2697	ENDIF
2698	ENDIF
2699	!
2700	ENDIF
2701	!
2702	!
2703	! We assume that if the forcing file is closed at the end of this subroutine this is a test
2704	! or initialisation of the grids. So we will dump the mask in a netCDF file for the user to
2705	! check.
2706	!
2707	dump_mask = closefile
2708	CALL forcing_checkindex(dump_mask, testvarname, testvar)
2709
2710	!
2711	!
2712	! 8.0 Close file if needed
2713	!
2714	IF ( closefile ) THEN
2715	CALL forcing_close()
2716	ENDIF
2717	!
2718	! Prepare variables to be returnned to calling subroutine.
2719	!
2720	iim_tmp = iim_glo
2721	jjm_tmp = jjm_glo
2722	nbpoint_tmp = nbpoint_glo
2723
2724	!
2725	! Clean up !
2726	!
2727	IF ( ALLOCATED(dimids) ) DEALLOCATE(dimids)
2728	IF ( ALLOCATED(londim_ids) ) DEALLOCATE(londim_ids)
2729	IF ( ALLOCATED(latdim_ids) ) DEALLOCATE(latdim_ids)
2730	IF ( ALLOCATED(testvar) ) DEALLOCATE(testvar)
2731	IF ( ALLOCATED(testvar2d) ) DEALLOCATE(testvar2d)
2732	IF ( ALLOCATED(contfrac2d) ) DEALLOCATE(contfrac2d)
2733	!
2734	END SUBROUTINE forcing_getglogrid
2735
2736	!! =============================================================================================================================
2737	!! SUBROUTINE: forcing_buildindex
2738	!!
2739	!>\BRIEF
2740	!!
2741	!! DESCRIPTION: When the forcing file does not contain compressed variables we need
2742	!! to build the land index variable from the mask defined by missing variables in
2743	!! a test variable.
2744	!!
2745	!! \n
2746	!_ ==============================================================================================================================
2747
2748	SUBROUTINE forcing_buildindex(var2d, var_missing, landonly, lindex, var)
2749	!
2750	! When the forcing file does not contain compressed variables we need
2751	! to build the land index variable from the mask defined by missing variables in
2752	! a test variable.
2753	!
2754	! Arguments
2755	!
2756	REAL(r_std), INTENT(in) :: var2d(:,:)
2757	REAL(r_std), INTENT(in) :: var_missing
2758	LOGICAL, INTENT(in) :: landonly
2759	INTEGER(i_std), INTENT(out) :: lindex(:)
2760	REAL(r_std), INTENT(out) :: var(:)
2761	!
2762	! Local
2763	!
2764	INTEGER(i_std) :: i,j,k
2765
2766	IF ( MAXVAL(var2d) >= var_missing ) THEN
2767	! Case when we have missing values in the vard2d
2768	k=0
2769	DO i=1,iim_glo
2770	DO j=1,jjm_glo
2771	IF ( landonly ) THEN
2772	IF ( var2d(i,j) /= var_missing .AND. var2d(i,j) > 0.0 ) THEN
2773	k = k + 1
2774	lindex(k) = (j-1)*iim_glo+i
2775	var(k) = var2d(i,j)
2776	ENDIF
2777	ELSE
2778	! When we take all point, no test is performed.
2779	k = k + 1
2780	lindex(k) = (j-1)*iim_glo+i
2781	var(k) = var2d(i,j)
2782	ENDIF
2783	ENDDO
2784	ENDDO
2785	ELSE
2786	! We suppose that this is land fraction variable
2787	k=0
2788	DO i=1,iim_glo
2789	DO j=1,jjm_glo
2790	IF ( landonly ) THEN
2791	IF ( var2d(i,j) > 0.0 ) THEN
2792	k = k + 1
2793	lindex(k) = (j-1)*iim_glo+i
2794	var(k) = var2d(i,j)
2795	ENDIF
2796	ELSE
2797	! When we take all point, no test is performed.
2798	k = k + 1
2799	lindex(k) = (j-1)*iim_glo+i
2800	var(k) = var2d(i,j)
2801	ENDIF
2802	ENDDO
2803	ENDDO
2804	ENDIF
2805	!
2806	!
2807	END SUBROUTINE forcing_buildindex
2808
2809	!! =============================================================================================================================
2810	!! SUBROUTINE: forcing_contfrac1d
2811	!!
2812	!>\BRIEF
2813	!!
2814	!! DESCRIPTION: This routine build the land/mask if needed and gets the contfrac variable from forcing file.
2815	!! Here we treat the case where the variables are compressed by gathering. Thus only
2816	!! land points are available in the file.
2817	!!
2818	!! \n
2819	!_ ==============================================================================================================================
2820
2821	SUBROUTINE forcing_contfrac1d(ifile, testvar_id, contfrac_id, testvar)
2822	!
2823	! This routine build the land/mask if needed and gets the contfrac variable from forcing file.
2824	! Here we treat the case where the variables are compressed by gathering. Thus only
2825	! land points are available in the file.
2826	!
2827	! ARGUMENTS
2828	!
2829	INTEGER(i_std), INTENT(in) :: ifile
2830	INTEGER(i_std), INTENT(in) :: testvar_id, contfrac_id
2831	REAL(r_std), DIMENSION(:), INTENT(inout) :: testvar
2832	!
2833	! LOCAL
2834	!
2835	INTEGER(i_std) :: it, iret
2836	INTEGER(i_std), DIMENSION(3) :: start, count
2837	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: contfrac2d
2838	!
2839	! First determine the contfrac variable
2840	!
2841
2842	IF ( contfrac_id > 0 ) THEN
2843	iret = NF90_INQUIRE_VARIABLE(ifile, contfrac_id, ndims=it)
2844	IF ( it == 1 ) THEN
2845	start = (/1,1,0/)
2846	count = (/nbpoint_glo,1,0/)
2847	iret = NF90_GET_VAR(ifile, contfrac_id, contfrac_glo, start, count)
2848	IF (iret /= NF90_NOERR) THEN
2849	WRITE(,) TRIM(nf90_strerror(iret))
2850	CALL ipslerr (3,'forcing_contfrac1d',"Error reading contfrac ", " ", " ")
2851	ENDIF
2852	ELSE IF ( it == 2 ) THEN
2853	ALLOCATE(contfrac2d(iim_glo,jjm_glo))
2854	start = (/1,1,0/)
2855	count = (/iim_glo,jjm_glo,0/)
2856	iret = NF90_GET_VAR(ifile, contfrac_id, contfrac2d)
2857	IF (iret /= NF90_NOERR) THEN
2858	WRITE(,) TRIM(nf90_strerror(iret))
2859	CALL ipslerr (3,'forcing_contfrac1d',"Error reading contfrac ", " ", " ")
2860	ENDIF
2861	CALL forcing_reindex(iim_glo, jjm_glo, contfrac2d, nbpoint_glo, contfrac_glo, lindex_glo)
2862	DEALLOCATE(contfrac2d)
2863	ELSE
2864	CALL ipslerr (3,'forcing_contfrac1d',"Contfrac has a dimension larger than 2. ", &
2865	"We do not know how to handle this.", " ")
2866	ENDIF
2867	ELSE
2868	contfrac_glo(:) = 1.0
2869	ENDIF
2870	!
2871	! Read our test variable
2872	!
2873	iret = NF90_INQUIRE_VARIABLE(ifile, testvar_id, ndims=it)
2874	IF ( it == 2 ) THEN
2875	start = (/1,1,0/)
2876	count = (/nbpoint_glo,1,0/)
2877	ELSE IF ( it == 3 ) THEN
2878	start = (/1,1,1/)
2879	count = (/nbpoint_glo,1,1/)
2880	ELSE
2881	CALL ipslerr (3,'forcing_contfrac1d',"Testvar has a dimension larger than 3.", &
2882	"We do not know how to handle this", " ")
2883	ENDIF
2884	iret = NF90_GET_VAR(ifile, testvar_id, testvar, start, count)
2885	IF (iret /= NF90_NOERR) THEN
2886	WRITE(,) TRIM(nf90_strerror(iret))
2887	CALL ipslerr (3,'forcing_contfrac1d',"Error reading testvar.", " ", " ")
2888	ENDIF
2889	!
2890	END SUBROUTINE forcing_contfrac1d
2891
2892	!! =============================================================================================================================
2893	!! SUBROUTINE: forcing_contfrac2d
2894	!!
2895	!>\BRIEF
2896	!!
2897	!! DESCRIPTION: This routine build the land/mask if needed and gets the contfrac variable from forcing file.
2898	!! Here we treat the case where the variables is 2D. Thus we also need to identify the land points.
2899	!! For this we can either use the contfrac variable or the test variable.
2900	!!
2901	!! \n
2902	!_ ==============================================================================================================================
2903
2904	SUBROUTINE forcing_contfrac2d(ifile, testvar_id, contfrac_id, testvar, contfrac, testvar_missing, contfrac_missing, nbland)
2905	!
2906	! This routine build the land/mask if needed and gets the contfrac variable from forcing file.
2907	! Here we treat the case where the variables is 2D. Thus we also need to identify the land points.
2908	! For this we can either use the contfrac variable or the test variable.
2909	!
2910	! ARGUMENTS
2911	!
2912	INTEGER(i_std), INTENT(in) :: ifile
2913	INTEGER(i_std), INTENT(in) :: testvar_id, contfrac_id
2914	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: testvar
2915	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: contfrac
2916	REAL(r_std), INTENT(in) :: testvar_missing
2917	REAL(r_std), INTENT(in) :: contfrac_missing
2918	INTEGER(i_std), INTENT(out) :: nbland
2919	!
2920	! LOCAL
2921	!
2922	INTEGER(i_std) :: i, j, it, iret
2923	INTEGER(i_std), DIMENSION(4) :: start, count
2924	!
2925	!
2926	nbland = 0
2927	!
2928	IF ( contfrac_id > 0 ) THEN
2929	!
2930	iret = NF90_INQUIRE_VARIABLE(ifile, contfrac_id, ndims=it)
2931	IF ( it == 2 ) THEN
2932	start = (/1,1,0,0/)
2933	count = (/iim_glo,jjm_glo,0,0/)
2934	iret = NF90_GET_VAR(ifile, contfrac_id, contfrac, start, count)
2935	IF (iret /= NF90_NOERR) THEN
2936	WRITE(,) TRIM(nf90_strerror(iret))
2937	CALL ipslerr (3,'forcing_contfrac2d',"Error reading contfrac.", " ", " ")
2938	ENDIF
2939	ELSE
2940	CALL ipslerr (3,'forcing_contfrac2d',"Contfrac has a dimension different of 1.", &
2941	"We do not know how to handle this.", " ")
2942	ENDIF
2943
2944	IF ( MAXVAL(contfrac) >= contfrac_missing ) THEN
2945	! We have missing values in contfrac and we use it to count number of land points
2946	DO i=1,iim_glo
2947	DO j=1,jjm_glo
2948	IF ( contfrac(i,j) /= contfrac_missing .AND. contfrac(i,j) > 0.0 ) THEN
2949	nbland = nbland + 1
2950	ENDIF
2951	ENDDO
2952	ENDDO
2953
2954	ELSE
2955	! Then ocean is fully contfrc=0 !
2956	DO i=1,iim_glo
2957	DO j=1,jjm_glo
2958	IF ( contfrac(i,j) > 0.0 ) THEN
2959	nbland = nbland + 1
2960	ENDIF
2961	ENDDO
2962	ENDDO
2963
2964	ENDIF
2965
2966	! If we did not find any land points on the map (i.e. iim_glo > 1 and jjm_glo > 1) then we
2967	! look for fractions larger then 0.
2968	!
2969	IF ( iim_glo > 1 .AND. jjm_glo > 1 .AND. nbland < 1 ) THEN
2970	DO i=1,iim_glo
2971	DO j=1,jjm_glo
2972	IF ( contfrac(i,j) > 0.0 ) THEN
2973	nbland = nbland + 1
2974	ENDIF
2975	ENDDO
2976	ENDDO
2977	ENDIF
2978
2979	! Did we get a result ?
2980	!
2981	IF ( iim_glo > 1 .AND. jjm_glo > 1 .AND. nbland < 1 ) THEN
2982	CALL ipslerr (3,'forcing_contfrac2d',"Contfrac was used to count the number of land points.", &
2983	& "We still have not found any land points when we looked for contfrac > 0.", " ")
2984	ENDIF
2985	!
2986	ELSE
2987	! Just so that we have no un-initialized variable
2988	contfrac(:,:) = 0.0
2989	ENDIF
2990	!
2991	IF ( testvar_id > 0 ) THEN
2992	!
2993	iret = NF90_INQUIRE_VARIABLE(ifile, testvar_id, ndims=it)
2994	IF ( it == 2 ) THEN
2995	start = (/1,1,0,0/)
2996	count = (/iim_glo,jjm_glo,0,0/)
2997	ELSE IF ( it == 3 ) THEN
2998	start = (/1,1,1,0/)
2999	count = (/iim_glo,jjm_glo,1,0/)
3000	ELSE IF ( it == 4 ) THEN
3001	start = (/1,1,1,1/)
3002	count = (/iim_glo,jjm_glo,1,1/)
3003	ELSE
3004	CALL ipslerr (3,'forcing_contfrac2d',"testvar has a dimension of 1 or larger than 4.", &
3005	"We do not know how to handle this.", " ")
3006	ENDIF
3007	iret = NF90_GET_VAR(ifile, testvar_id, testvar, start, count)
3008	IF (iret /= NF90_NOERR) THEN
3009	WRITE(,) TRIM(nf90_strerror(iret))
3010	CALL ipslerr (3,'forcing_contfrac2d',"Error reading testvar.", " ", " ")
3011	ENDIF
3012	!
3013	! IF with count frac we did not get the number of land points, let us try it here
3014	!
3015	IF ( nbland < 1 ) THEN
3016	DO i=1,iim_glo
3017	DO j=1,jjm_glo
3018	IF ( testvar(i,j) < testvar_missing ) THEN
3019	nbland = nbland + 1
3020	! Add infor to contfrac
3021	IF ( contfrac_id < 0 ) THEN
3022	contfrac(i,j) = 1.0
3023	ENDIF
3024	ENDIF
3025	ENDDO
3026	ENDDO
3027	ENDIF
3028	!
3029	!
3030	! Did we get a result here ?
3031	!
3032	IF ( iim_glo > 1 .AND. jjm_glo > 1 .AND. nbland < 1 ) THEN
3033	CALL ipslerr (3,'forcing_contfrac2d',"Contfrac and testvar were used to count the number", &
3034	& "of land points. We have not found any land points.", " ")
3035	ENDIF
3036	!
3037	ENDIF
3038	!
3039	END SUBROUTINE forcing_contfrac2d
3040
3041	!! =============================================================================================================================
3042	!! SUBROUTINE: forcing_checkindex
3043	!!
3044	!>\BRIEF
3045	!!
3046	!! DESCRIPTION: For ORCHIDEE its paralelisation requires that the land points are ordered
3047	!! in such a way that the longitude runs fastest. That means that we go over the
3048	!! globle filling one line after the other.
3049	!! As this might not be the case in a compressed vector of land points, we need to
3050	!! put all the points on the 2d grid and then scan them in the right order.
3051	!! The reindexing is prepared here.
3052	!!
3053	!! \n
3054	!_ ==============================================================================================================================
3055
3056	SUBROUTINE forcing_checkindex(dump_mask, testvarname, testvar)
3057	!
3058	! For ORCHIDEE its paralelisation requires that the land points are ordered
3059	! in such a way that the longitude runs fastest. That means that we go over the
3060	! globle filling one line after the other.
3061	! As this might not be the case in a compressed vector of land points, we need to
3062	! put all the points on the 2d grid and then scan them in the right order.
3063	! The reindexing is prepared here.
3064	!
3065	LOGICAL :: dump_mask
3066	CHARACTER(LEN=*) :: testvarname
3067	REAL(r_std) :: testvar(:)
3068	!
3069	INTEGER(i_std) :: j, i, ik
3070	REAL(r_std), ALLOCATABLE, DIMENSION(:) :: testvar_reind
3071	!
3072	!
3073	!
3074	! Get the indices of the land points in the focing file
3075	!
3076	IF ( .NOT. ALLOCATED(reindex_glo)) ALLOCATE(reindex_glo(nbpoint_glo))
3077	IF ( .NOT. ALLOCATED(origind)) ALLOCATE(origind(iim_glo,jjm_glo))
3078	!
3079	! Find the origine of each point in the gathered vector on the xy grid.
3080	!
3081	origind(:,:) = -1
3082	mask_glo(:,:) = 0
3083	DO ik=1,nbpoint_glo
3084	j = INT((lindex_glo(ik)-1)/iim_glo)+1
3085	i = (lindex_glo(ik)-(j-1)*iim_glo)
3086	origind(i,j) = ik
3087	mask_glo(i,j) = 1
3088	ENDDO
3089	!
3090	! Prepare a reindexing array so that the vector goes in the right order : longitude runs
3091	! faster than the latitude. Put then also the right information into lindex_glo.
3092	!
3093	ik=0
3094	DO j=1,jjm_glo
3095	DO i=1,iim_glo
3096	IF ( origind(i,j) > zero ) THEN
3097	ik = ik+1
3098	reindex_glo(ik) = origind(i,j)
3099	lindex_glo(ik) = (j-1)*iim_glo+i
3100	ENDIF
3101	ENDDO
3102	ENDDO
3103	!
3104	!
3105	! Write the mask and a test variable to a file so that the user can check that all is OK
3106	!
3107	IF ( dump_mask) THEN
3108	!
3109	! Scatter the test variable and save it in the file
3110	!
3111	WRITE(,) MINVAL(testvar), "<< test variable ", TRIM(testvarname), " <<", MAXVAL(testvar)
3112	ALLOCATE(testvar_reind(nbpoint_glo))
3113	!
3114	CALL forcing_reindex(nbpoint_glo, testvar, nbpoint_glo, testvar_reind, reindex_glo)
3115	!
3116
3117	CALL forcing_writetestvar("forcing_mask_glo.nc", iim_glo, jjm_glo, nbpoint_glo, &
3118	& lon_glo(:,1), lat_glo(1,:), lindex_glo, mask_glo, &
3119	& testvarname, testvar_reind)
3120	!
3121	ENDIF
3122	!
3123	! Clean up !
3124	!
3125	IF ( ALLOCATED(testvar_reind) ) DEALLOCATE(testvar_reind)
3126	!
3127	END SUBROUTINE forcing_checkindex
3128
3129	!! =============================================================================================================================
3130	!! SUBROUTINE: forcing_writetestvar
3131	!!
3132	!>\BRIEF Write the mask and a test variable to a netCDF file.
3133	!!
3134	!! DESCRIPTION: This routine allows to test if the variables read from the forcing files is well read.
3135	!! Typically the forcing is compressed by gathering and thus it is a safe practice
3136	!! to verify that the un-compression is done correctly and that all points end-up in the
3137	!! right place on the global lat/lon grid.
3138	!!
3139	!! \n
3140	!_ ==============================================================================================================================
3141
3142	SUBROUTINE forcing_writetestvar(ncdffile, iim, jjm, nbland, lon, lat, lindex, mask, varname, var)
3143	!
3144	! Write the mask and a test variable to a netCDF file
3145	!
3146	! ARGUMENTS
3147	!
3148	CHARACTER(LEN=*), INTENT(in) :: ncdffile
3149	INTEGER(i_std), INTENT(in) :: iim, jjm, nbland
3150	REAL(r_std), INTENT(in) :: lon(iim), lat(jjm)
3151	INTEGER(i_std), INTENT(in) :: lindex(nbland)
3152	INTEGER(i_std), INTENT(in) :: mask(iim,jjm)
3153	CHARACTER(LEN=*), INTENT(in) :: varname
3154	REAL(r_std), INTENT(in) :: var(nbland)
3155	!
3156	! Local
3157	!
3158	INTEGER(i_std) :: ik, i, j
3159	INTEGER(i_std) :: iret, nlonid, nlatid, varid, fid, ierr, iland
3160	INTEGER(i_std) :: testid
3161	INTEGER(i_std), DIMENSION(2) :: lolaid
3162	REAL(r_std) :: test_scat(iim,jjm)
3163	!
3164	!
3165	test_scat(:,:) = NF90_FILL_REAL
3166	CALL forcing_reindex(nbland, var, iim, jjm, test_scat, lindex)
3167	!
3168	iret = NF90_CREATE(ncdffile, NF90_WRITE, fid)
3169	IF (iret /= NF90_NOERR) THEN
3170	CALL ipslerr (3,'forcing_writetestvar',"Error opening the output file : ", ncdffile, " ")
3171	ENDIF
3172	!
3173	! Define dimensions
3174	!
3175	iret = NF90_DEF_DIM(fid,'lon',iim,lolaid(1))
3176	iret = NF90_DEF_DIM(fid,'lat',jjm,lolaid(2))
3177	!
3178	!
3179	iret = NF90_DEF_VAR(fid,"lon",NF90_REAL4,lolaid(1),nlonid)
3180	iret = NF90_PUT_ATT(fid,nlonid,'axis',"X")
3181	iret = NF90_PUT_ATT(fid,nlonid,'standard_name',"longitude")
3182	iret = NF90_PUT_ATT(fid,nlonid,'units',"degree_east")
3183	iret = NF90_PUT_ATT(fid,nlonid,'valid_min',MINVAL(lon_glo))
3184	iret = NF90_PUT_ATT(fid,nlonid,'valid_max',MAXVAL(lon_glo))
3185	iret = NF90_PUT_ATT(fid,nlonid,'long_name',"Longitude")
3186	!
3187	iret = NF90_DEF_VAR(fid,"lat",NF90_REAL4,lolaid(2),nlatid)
3188	iret = NF90_PUT_ATT(fid,nlatid,'axis',"Y")
3189	iret = NF90_PUT_ATT(fid,nlatid,'standard_name',"latitude")
3190	iret = NF90_PUT_ATT(fid,nlatid,'units',"degree_north")
3191	iret = NF90_PUT_ATT(fid,nlatid,'valid_min',MINVAL(lat_glo))
3192	iret = NF90_PUT_ATT(fid,nlatid,'valid_max',MAXVAL(lat_glo))
3193	iret = NF90_PUT_ATT(fid,nlatid,'long_name',"Latitude")
3194	!
3195	iret = NF90_DEF_VAR(fid,"mask",NF90_REAL4,lolaid,varid)
3196	!
3197	iret = NF90_DEF_VAR(fid,TRIM(varname),NF90_REAL4,lolaid,testid)
3198	iret = NF90_PUT_ATT(fid,testid,'_FillValue',NF90_FILL_REAL)
3199	iret = NF90_PUT_ATT(fid,testid,'missing_value',NF90_FILL_REAL)
3200	!
3201	iret = NF90_ENDDEF (fid)
3202	IF (iret /= NF90_NOERR) THEN
3203	WRITE(,) TRIM(nf90_strerror(iret))
3204	CALL ipslerr (3,'forcing_writetestvar',"Error ending definitions in file : ", ncdffile, " ")
3205	ENDIF
3206	!
3207	! Write variables
3208	!
3209	iret = NF90_PUT_VAR(fid,nlonid,lon)
3210	iret = NF90_PUT_VAR(fid,nlatid,lat)
3211	iret = NF90_PUT_VAR(fid,varid,REAL(mask))
3212	iret = NF90_PUT_VAR(fid,testid,test_scat)
3213	!
3214	! Close file
3215	!
3216	iret = NF90_CLOSE(fid)
3217	IF (iret /= NF90_NOERR) THEN
3218	CALL ipslerr (3,'forcing_writetestvar',"Error closing the output file : ", ncdffile, " ")
3219	ENDIF
3220	!
3221	END SUBROUTINE forcing_writetestvar
3222
3223	!! =============================================================================================================================
3224	!! SUBROUTINE: forcing_zoomgrid
3225	!!
3226	!>\BRIEF We zoom into the region requested by the user.
3227	!!
3228	!! DESCRIPTION: Get the area to be zoomed and the sizes of arrays we will need.
3229	!! This subroutine fills the *_loc variables.
3230	!! If requested it will dump a test vraible into a netCDF file.
3231	!!
3232	!! \n
3233	!_ ==============================================================================================================================
3234
3235	SUBROUTINE forcing_zoomgrid (zoom_lon, zoom_lat, filename, dumpncdf)
3236	!
3237	! Get the area to be zoomed and the sizes of arrays we will need.
3238	! This subroutine fills the *_loc variables.
3239	! If requested it will dump a test vraible into a netCDF file.
3240	!
3241	! ARGUMENTS
3242	!
3243	REAL(r_std), DIMENSION(2), INTENT(in) :: zoom_lon, zoom_lat
3244	CHARACTER(LEN=*), INTENT(in) :: filename
3245	LOGICAL, INTENT(in) :: dumpncdf
3246	!
3247	! LOCAL
3248	!
3249	INTEGER(i_std) :: i, j, ic, jc, ik, ig
3250	REAL(r_std) :: dx, dy, coslat
3251	REAL(r_std) :: lon_tmp(iim_glo), lat_tmp(jjm_glo)
3252	REAL(r_std) :: longlo_tmp(iim_glo,jjm_glo)
3253	REAL(r_std), ALLOCATABLE, DIMENSION(:) :: lon_val, lat_val
3254	INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:,:) :: zoom_index
3255	!
3256	INTEGER(i_std) :: iret, force_id, iv
3257	INTEGER(i_std), DIMENSION(1) :: imin, jmin
3258	INTEGER(i_std), DIMENSION(2) :: start, count
3259	INTEGER(i_std), DIMENSION(3) :: start2d, count2d
3260	REAL(r_std), ALLOCATABLE, DIMENSION(:) :: readvar, zoomedvar
3261	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: readvar2d
3262	INTEGER(i_std), ALLOCATABLE, DIMENSION(:) :: index_glotoloc
3263	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: lalo
3264	CHARACTER(LEN=8) :: testvarname="Tair"
3265	!
3266	! 0.0 Verify we are on the root processor
3267	!
3268	IF ( .NOT. is_root_prc ) THEN
3269	CALL ipslerr (3,'forcing_zoomgrid'," This routine can only be called on the root processor.", " ", " ")
3270	ENDIF
3271	!
3272	! 0.1 Inform the user
3273	!
3274	WRITE(,) "Zoom forcing : lon = ", zoom_lon
3275	WRITE(,) "Zoom forcing : lat = ", zoom_lat
3276	!
3277	! Some forcing files have longitudes going from 0 to 360. This code works on the
3278	! -180 to 180 longitude grid. So if needed we transform the longitudes of the global grid.
3279	!
3280	IF ( MAXVAL(lon_glo) <= 180.0 ) THEN
3281	longlo_tmp=lon_glo
3282	ELSE
3283	DO i=1,iim_glo
3284	DO j=1,jjm_glo
3285	IF ( lon_glo(i,j) <= 180.0 ) THEN
3286	longlo_tmp(i,j) = lon_glo(i,j)
3287	ELSE
3288	longlo_tmp(i,j) = lon_glo(i,j)-360
3289	ENDIF
3290	ENDDO
3291	ENDDO
3292	ENDIF
3293	!
3294	! See if we need to zoom
3295	!
3296	IF (MINVAL(zoom_lon) > MINVAL(longlo_tmp) .OR. MAXVAL(zoom_lon) < MAXVAL(longlo_tmp) .OR.&
3297	& MINVAL(zoom_lat) > MINVAL(lat_glo) .OR. MAXVAL(zoom_lat) < MAXVAL(lat_glo) ) THEN
3298	zoom_forcing = .TRUE.
3299	ENDIF
3300	!
3301	! Determine the size in x and y of the zoom
3302	!
3303	IF ( zoom_forcing ) THEN
3304	!
3305	! Working with the hypothesis it is a regular global grid and bring it back to the -180 to 180 interval
3306	! if needed.
3307	!
3308	lon_tmp(:) = longlo_tmp(:,1)
3309	lat_tmp(:) = lat_glo(1,:)
3310	!
3311	DO i=1,iim_glo
3312	IF ( lon_tmp(i) <= MINVAL(zoom_lon) .OR. lon_tmp(i) >= MAXVAL(zoom_lon) ) THEN
3313	lon_tmp(i) = 0.0
3314	ELSE
3315	lon_tmp(i) = 1.0
3316	ENDIF
3317	ENDDO
3318	DO j=1,jjm_glo
3319	IF ( lat_tmp(j) <= MINVAL(zoom_lat) .OR. lat_tmp(j) >= MAXVAL(zoom_lat) ) THEN
3320	lat_tmp(j) = 0.0
3321	ELSE
3322	lat_tmp(j) = 1.0
3323	ENDIF
3324	ENDDO
3325	iim_loc = NINT(SUM(lon_tmp))
3326	jjm_loc = NINT(SUM(lat_tmp))
3327	ELSE
3328	iim_loc = iim_glo
3329	jjm_loc = jjm_glo
3330	lon_tmp(:) = 1.0
3331	lat_tmp(:) = 1.0
3332	ENDIF
3333	!
3334	! Determine the number of land points in the zoom
3335	!
3336	IF ( .NOT. ALLOCATED(lon_loc) ) ALLOCATE(lon_loc(iim_loc,jjm_loc))
3337	IF ( .NOT. ALLOCATED(lat_loc) ) ALLOCATE(lat_loc(iim_loc,jjm_loc))
3338	IF ( .NOT. ALLOCATED(mask_loc) ) ALLOCATE(mask_loc(iim_loc,jjm_loc))
3339	IF ( .NOT. ALLOCATED(zoom_index) ) ALLOCATE(zoom_index(iim_loc,jjm_loc,2))
3340	!
3341	IF ( .NOT. ALLOCATED(lon_val)) ALLOCATE(lon_val(iim_loc))
3342	IF ( .NOT. ALLOCATED(lat_val)) ALLOCATE(lat_val(jjm_loc))
3343	!
3344	! Create our new lat/lon system which is in the -180/180 range and South to North and West to East
3345	!
3346	ic=0
3347	DO i=1,iim_glo
3348	IF ( lon_tmp(i) > 0 ) THEN
3349	ic = ic+1
3350	lon_val(ic) = longlo_tmp(i,1)
3351	ENDIF
3352	ENDDO
3353	jc=0
3354	DO j=1,jjm_glo
3355	IF ( lat_tmp(j) > 0 ) THEN
3356	jc = jc+1
3357	lat_val(jc) = lat_glo(1,j)
3358	ENDIF
3359	ENDDO
3360	CALL sort(lon_val, iim_loc)
3361	CALL sort(lat_val, jjm_loc)
3362	!
3363	! Now find the correspondance between the zoomed & re-ordered grid and the global one.
3364	!
3365	DO i=1,iim_loc
3366	DO j=1,jjm_loc
3367	!
3368	imin=MINLOC(ABS(longlo_tmp(:,1)-lon_val(i)))
3369	jmin=MINLOC(ABS(lat_glo(1,:)-lat_val(j)))
3370	!
3371	lon_loc(i,j) = longlo_tmp(imin(1),jmin(1))
3372	lat_loc(i,j) = lat_glo(imin(1),jmin(1))
3373	mask_loc(i,j) = mask_glo(imin(1),jmin(1))
3374	!
3375	zoom_index(i,j,1) = imin(1)
3376	zoom_index(i,j,2) = jmin(1)
3377	!
3378	ENDDO
3379	ENDDO
3380	!
3381	nbpoint_loc = SUM(mask_loc)
3382
3383
3384	IF ( .NOT. zoom_forcing .AND. nbpoint_loc .NE. nbpoint_glo) THEN
3385	WRITE(,) "We have not zoomed into the forcing file still we get a number of"
3386	WRITE(,) "land points that is different from what we have in the forcing file."
3387	STOP "forcing_zoomgrid"
3388	ENDIF
3389	!
3390	IF ( .NOT. ALLOCATED(lindex_loc)) ALLOCATE(lindex_loc(nbpoint_loc))
3391	IF ( .NOT. ALLOCATED(reindex_loc)) ALLOCATE(reindex_loc(nbpoint_loc))
3392	IF ( .NOT. ALLOCATED(contfrac_loc)) ALLOCATE(contfrac_loc(nbpoint_loc))
3393	!
3394	IF ( .NOT. ALLOCATED(reindex2d_loc)) ALLOCATE(reindex2d_loc(nbpoint_loc,2))
3395	IF ( .NOT. ALLOCATED(index_glotoloc)) ALLOCATE(index_glotoloc(nbpoint_glo))
3396	IF ( .NOT. ALLOCATED(lalo)) ALLOCATE(lalo(nbpoint_loc,2))
3397	!
3398	! Do the actual zoom on the grid
3399	!
3400	! Set indices of all points as non existant so that we can fill in as we zoom the
3401	! indices of the points which exist.
3402	index_glotoloc(:) = -1
3403	!
3404	ik = 0
3405	!
3406	! Loop only over the zoomed grid
3407	!
3408	! Why does the inner loop need to be ic for the pralalisation ????
3409	!
3410	DO jc=1,jjm_loc
3411	DO ic=1,iim_loc
3412	!
3413	! Point back from the local to the original global i*j grid
3414	!
3415	i = zoom_index(ic,jc,1)
3416	j = zoom_index(ic,jc,2)
3417	!
3418	IF ( origind(i,j) > 0 ) THEN
3419	ik = ik+1
3420	! index of the points in the local grid
3421	lindex_loc(ik) = (jc-1)*iim_loc+ic
3422	!
3423	! For land points, the index of global grid is saved for the this point on the local grid
3424	reindex_loc(ik) = origind(i,j)
3425	!
3426	! Keep also the i and j of the global grid for this land point on the local grid
3427	reindex2d_loc(ik,1) = i
3428	reindex2d_loc(ik,2) = j
3429	!
3430	! Keep the reverse : on the global grid the location where we put the value of the local grid.
3431	index_glotoloc(origind(i,j)) = ik
3432	!
3433	contfrac_loc(ik) = contfrac_glo(origind(i,j))
3434	!
3435	lalo(ik,1) = lat_glo(i,j)
3436	lalo(ik,2) = longlo_tmp(i,j)
3437	!
3438	ENDIF
3439	ENDDO
3440	ENDDO
3441	!
3442	!
3443	nbland_loc = 0
3444	DO ik=1, SIZE(contfrac_loc)
3445	IF (contfrac_loc(ik) > 0.0) THEN
3446	nbland_loc = nbland_loc + 1.0
3447	ENDIF
3448	ENDDO
3449	!
3450	!
3451	ncdfstart = MINVAL(reindex_loc)
3452	reindex_loc(:) = reindex_loc(:)-ncdfstart+1
3453	ncdfcount = MAXVAL(reindex_loc)
3454	!
3455	! Compute the areas and the corners on the grid over which we will run ORCHIDEE.
3456	! As this module is dedicated for regular lat/lon forcing we know that we can compute these
3457	! variables without further worries.
3458	!
3459	IF ( .NOT. ALLOCATED(area_loc)) ALLOCATE(area_loc(iim_loc,jjm_loc))
3460	IF ( .NOT. ALLOCATED(corners_loc)) ALLOCATE(corners_loc(iim_loc,jjm_loc,4,2))
3461	!
3462	! Treat first the longitudes
3463	!
3464	DO j=1,jjm_loc
3465	dx = zero
3466	DO i=1,iim_loc-1
3467	dx = dx+ABS(lon_loc(i,j)-lon_loc(i+1,j))
3468	ENDDO
3469	dx = dx/(iim_loc-1)
3470	DO i=1,iim_loc
3471	corners_loc(i,j,1,1) = lon_loc(i,j)-dx/2.0
3472	corners_loc(i,j,2,1) = lon_loc(i,j)+dx/2.0
3473	corners_loc(i,j,3,1) = lon_loc(i,j)+dx/2.0
3474	corners_loc(i,j,4,1) = lon_loc(i,j)-dx/2.0
3475	ENDDO
3476	ENDDO
3477	!
3478	! Now treat the latitudes
3479	!
3480	DO i=1,iim_loc
3481	dy = zero
3482	DO j=1,jjm_loc-1
3483	dy = dy + ABS(lat_loc(i,j)-lat_loc(i,j+1))
3484	ENDDO
3485	dy = dy/(jjm_loc-1)
3486	DO j=1,jjm_loc
3487	corners_loc(i,j,1,2) = lat_loc(i,j)+dy/2.0
3488	corners_loc(i,j,2,2) = lat_loc(i,j)+dy/2.0
3489	corners_loc(i,j,3,2) = lat_loc(i,j)-dy/2.0
3490	corners_loc(i,j,4,2) = lat_loc(i,j)-dy/2.0
3491	ENDDO
3492	ENDDO
3493	!
3494	! Compute the areas of the grid (using the simplification that the grid is regular in lon/lat).
3495	!
3496	DO i=1,iim_loc
3497	DO j=1,jjm_loc
3498	coslat = MAX( COS(lat_loc(i,j) * pi/180. ), mincos )
3499	dx = ABS(corners_loc(i,j,2,1) - corners_loc(i,j,1,1)) * pi/180. * R_Earth * coslat
3500	dy = ABS(corners_loc(i,j,1,2) - corners_loc(i,j,3,2)) * pi/180. * R_Earth
3501	area_loc(i,j) = dx*dy
3502	ENDDO
3503	ENDDO
3504	!
3505	! If requested we read a variable, zoomin and dump the result into a test file.
3506	!
3507	IF ( dumpncdf ) THEN
3508	iret = NF90_OPEN (filename, NF90_NOWRITE, force_id)
3509	IF (iret /= NF90_NOERR) THEN
3510	CALL ipslerr (3,'forcing_zoomgrid',"Error opening the forcing file :", filename, " ")
3511	ENDIF
3512	!
3513	ALLOCATE(readvar(ncdfcount), readvar2d(iim_glo,jjm_glo), zoomedvar(nbpoint_loc))
3514	!
3515	iret = NF90_INQ_VARID(force_id, TRIM(testvarname), iv)
3516	IF (iret /= NF90_NOERR) THEN
3517	CALL ipslerr (3,'forcing_zoomgrid',"Could not find variable Tair in file."," "," ")
3518	ENDIF
3519
3520	IF ( compressed ) THEN
3521	!
3522	start(1) = ncdfstart
3523	start(2) = 1
3524	count(1) = ncdfcount
3525	count(2) = 1
3526	!
3527	iret = NF90_GET_VAR(force_id, iv, readvar, start, count)
3528	IF (iret /= NF90_NOERR) THEN
3529	CALL ipslerr (3,'forcing_zoomgrid',"Could not read compressed variable Tair from file."," "," ")
3530	ENDIF
3531	CALL forcing_reindex(ncdfcount, readvar, nbpoint_loc, zoomedvar, reindex_loc)
3532	!
3533	ELSE
3534	!
3535	start2d(1) = 1
3536	start2d(2) = 1
3537	start2d(3) = 1
3538	count2d(1) = iim_glo
3539	count2d(2) = jjm_glo
3540	count2d(3) = 1
3541	!
3542	iret = NF90_GET_VAR(force_id, iv, readvar2d, start2d, count2d)
3543	IF (iret /= NF90_NOERR) THEN
3544	CALL ipslerr (3,'forcing_zoomgrid',"Could not read variable Tair from file."," "," ")
3545	ENDIF
3546	CALL forcing_reindex(iim_glo, jjm_glo, readvar2d, nbpoint_loc, zoomedvar, reindex2d_loc)
3547	!
3548	ENDIF
3549	!
3550	CALL forcing_writetestvar("forcing_mask_loc.nc", iim_loc, jjm_loc, nbpoint_loc, &
3551	& lon_loc(:,1), lat_loc(1,:), lindex_loc, mask_loc, &
3552	& TRIM(testvarname), zoomedvar)
3553	!
3554	ENDIF
3555	!
3556	! Clean up
3557	!
3558	IF ( ALLOCATED(readvar) ) DEALLOCATE(readvar)
3559	IF ( ALLOCATED(readvar2d) ) DEALLOCATE(readvar2d)
3560	IF ( ALLOCATED(zoomedvar) ) DEALLOCATE(zoomedvar)
3561	IF ( ALLOCATED(index_glotoloc) ) DEALLOCATE(index_glotoloc)
3562	IF ( ALLOCATED(lalo) ) DEALLOCATE(lalo)
3563	!
3564	END SUBROUTINE forcing_zoomgrid
3565
3566	!! =============================================================================================================================
3567	!! SUBROUTINE: forcing_givegridsize
3568	!!
3569	!>\BRIEF Routine which exports the size of the grid on which the model will run, i.e. the zoomed grid.
3570	!!
3571	!! DESCRIPTION: This is needed to transfer the grid information from this module to the glogrid.f90 module.
3572	!!
3573	!! \n
3574	!_ ==============================================================================================================================
3575
3576	SUBROUTINE forcing_givegridsize (iim, jjm, nblindex)
3577	!
3578	! Provides the size of the grid to be used to the calling program
3579	!
3580	! Size of the x and y direction of the zoomed area
3581	INTEGER(i_std), INTENT(out) :: iim, jjm
3582	! Number of land points in the zoomed area
3583	INTEGER(i_std), INTENT(out) :: nblindex
3584	!
3585	IF ( .NOT. is_root_prc ) THEN
3586	CALL ipslerr (3,'forcing_givegridsize'," This routine can only be called on the root processor.", &
3587	& "The information requested is only available on root processor.", " ")
3588	ENDIF
3589	!
3590	iim = iim_loc
3591	jjm = jjm_loc
3592	nblindex = nbland_loc
3593	!
3594	END SUBROUTINE forcing_givegridsize
3595
3596	!! =============================================================================================================================
3597	!! SUBROUTINE: forcing_
3598	!!
3599	!>\BRIEF
3600	!!
3601	!! DESCRIPTION:
3602	!!
3603	!! \n
3604	!_ ==============================================================================================================================
3605
3606	SUBROUTINE forcing_vertical(force_id)
3607	!
3608	!- This subroutine explores the forcing file to decide what information is available
3609	!- on the vertical coordinate.
3610	!
3611	INTEGER, INTENT(IN) :: force_id
3612	!
3613	INTEGER(i_std) :: iret, ireta, iretb
3614	!
3615	INTEGER(i_std) :: sigma_id = -1, sigma_uv_id = -1
3616	INTEGER(i_std) :: hybsiga_id = -1, hybsiga_uv_id = -1
3617	INTEGER(i_std) :: hybsigb_id = -1, hybsigb_uv_id = -1
3618	INTEGER(i_std) :: levels_id = -1, levels_uv_id = -1
3619	INTEGER(i_std) :: height_id = -1, height_uv_id = -1
3620	INTEGER(i_std) :: lev_id = -1
3621	!
3622	LOGICAL :: var_exists, vara_exists, varb_exists, varuv_exists
3623	LOGICAL :: foundvar
3624	LOGICAL :: levlegacy
3625	!
3626	!- Set all the defaults
3627	!
3628	zfixed=.FALSE.
3629	zsigma=.FALSE.
3630	zhybrid=.FALSE.
3631	zlevels=.FALSE.
3632	zheight=.FALSE.
3633	zsamelev_uv = .TRUE.
3634	levlegacy = .FALSE.
3635	!
3636	foundvar = .FALSE.
3637	!
3638	!- We have a forcing file to explore so let us see if we find any of the conventions
3639	!- which allow us to find the height of T,Q,U and V.
3640	!
3641	IF ( force_id > 0 ) THEN
3642	!
3643	! Case for sigma levels
3644	!
3645	IF ( .NOT. foundvar ) THEN
3646	ireta = NF90_INQ_VARID(force_id, 'Sigma', sigma_id)
3647	IF ( (sigma_id >= 0) .AND. (ireta == NF90_NOERR) ) THEN
3648	foundvar = .TRUE.
3649	zsigma = .TRUE.
3650	iretb = NF90_INQ_VARID(force_id, 'Sigma_uv', sigma_uv_id)
3651	IF ( (sigma_uv_id >= 0) .OR. (iretb == NF90_NOERR) ) zsamelev_uv = .FALSE.
3652	ENDIF
3653	ENDIF
3654	!
3655	! Case for Hybrid levels
3656	!
3657	IF ( .NOT. foundvar ) THEN
3658	var_exists = .FALSE.
3659	varuv_exists = .FALSE.
3660	ireta = NF90_INQ_VARID(force_id, 'HybSigA', hybsiga_id)
3661	IF ( (hybsiga_id >= 0 ) .AND. (ireta == NF90_NOERR) ) THEN
3662	iretb = NF90_INQ_VARID(force_id, 'HybSigB', hybsigb_id)
3663	IF ( (hybsigb_id >= 0 ) .AND. (iretb == NF90_NOERR) ) THEN
3664	var_exists=.TRUE.
3665	ELSE
3666	CALL ipslerr ( 3, 'forcing_vertical','Missing the B coefficient for', &
3667	& 'Hybrid vertical levels for T and Q','stop')
3668	ENDIF
3669	ENDIF
3670	ireta = NF90_INQ_VARID(force_id, 'HybSigA_uv', hybsiga_uv_id)
3671	IF ( (hybsiga_uv_id >= 0 ) .AND. (ireta == NF90_NOERR) ) THEN
3672	iretb = NF90_INQ_VARID(force_id, 'HybSigB_uv', hybsigb_uv_id)
3673	IF ( (hybsigb_uv_id >= 0 ) .AND. (iretb == NF90_NOERR) ) THEN
3674	varuv_exists=.TRUE.
3675	ELSE
3676	CALL ipslerr ( 3, 'forcing_vertical','Missing the B coefficient for', &
3677	& 'Hybrid vertical levels for U and V','stop')
3678	ENDIF
3679	ENDIF
3680	IF ( var_exists ) THEN
3681	foundvar = .TRUE.
3682	zhybrid = .TRUE.
3683	IF ( varuv_exists ) zsamelev_uv = .FALSE.
3684	ENDIF
3685	ENDIF
3686	!
3687	! Case for levels (i.e. a 2d time varying field with the height in meters)
3688	!
3689	IF ( .NOT. foundvar ) THEN
3690	ireta = NF90_INQ_VARID(force_id, 'Levels', levels_id)
3691	IF ( (levels_id >= 0 ) .AND. (ireta == NF90_NOERR) ) THEN
3692	foundvar = .TRUE.
3693	zlevels = .TRUE.
3694	iretb = NF90_INQ_VARID(force_id, 'Levels_uv', levels_uv_id)
3695	IF ( (levels_uv_id >= 0 ) .AND. (iretb == NF90_NOERR) ) zsamelev_uv = .FALSE.
3696	ENDIF
3697	ENDIF
3698	!
3699	! Case where a fixed height is provided in meters
3700	!
3701	IF ( .NOT. foundvar ) THEN
3702	ireta = NF90_INQ_VARID(force_id, 'Height_Lev1', height_id)
3703	IF ( (height_id >= 0 ) .AND. (ireta == NF90_NOERR) ) THEN
3704	foundvar = .TRUE.
3705	zheight = .TRUE.
3706	iretb = NF90_INQ_VARID(force_id, 'Height_Levuv', height_uv_id)
3707	IF ( (height_uv_id >= 0 ) .AND. (iretb == NF90_NOERR) ) zsamelev_uv = .FALSE.
3708	ENDIF
3709	ENDIF
3710	!
3711	! Case where a fixed height is provided in meters in the lev variable
3712	!
3713	IF ( .NOT. foundvar ) THEN
3714	ireta = NF90_INQ_VARID(force_id, 'lev', lev_id)
3715	IF ( (lev_id >= 0 ) .AND. (ireta == NF90_NOERR) ) THEN
3716	foundvar = .TRUE.
3717	zheight = .TRUE.
3718	levlegacy = .TRUE.
3719	ENDIF
3720	ENDIF
3721	!
3722	ENDIF
3723	!
3724	! We found forcing variables so we need to extract the values if we are dealing with constant values (i.e. all
3725	! except the case zlevels
3726	!
3727	IF ( foundvar .AND. .NOT. zlevels ) THEN
3728	!
3729	IF ( zheight ) THEN
3730	!
3731	! Constant height
3732	!
3733	IF ( levlegacy ) THEN
3734	iret = NF90_GET_VAR(force_id, lev_id, zlev_fixed)
3735	IF ( iret /= NF90_NOERR ) THEN
3736	CALL ipslerr ( 3, 'forcing_vertical','Attempted to read variable lev from forcing file in legacy mode', &
3737	& 'NF90_GET_VAR failed.','stop')
3738	ENDIF
3739	ELSE
3740	iret = NF90_GET_VAR(force_id, height_id, zlev_fixed)
3741	IF ( iret /= NF90_NOERR ) THEN
3742	CALL ipslerr ( 3, 'forcing_vertical','Attempted to read variable Height_Lev1 from forcing file', &
3743	& 'NF90_GET_VAR failed.','stop')
3744	ENDIF
3745	IF ( .NOT. zsamelev_uv ) THEN
3746	iret = NF90_GET_VAR(force_id, height_uv_id, zlevuv_fixed)
3747	IF ( iret /= NF90_NOERR ) THEN
3748	CALL ipslerr ( 3, 'forcing_vertical','Attempted to read variable Height_Levuv from forcing file', &
3749	& 'NF90_GET_VAR failed.','stop')
3750	ENDIF
3751	ENDIF
3752	ENDIF
3753	WRITE(,) "forcing_vertical : case ZLEV : Read from forcing file :", zlev_fixed, zlevuv_fixed
3754	!
3755	ELSE IF ( zsigma .OR. zhybrid ) THEN
3756	!
3757	! Sigma or hybrid levels
3758	!
3759	IF ( zsigma ) THEN
3760	iret = NF90_GET_VAR(force_id, sigma_id, zhybrid_b)
3761	zhybrid_a = zero
3762	IF ( .NOT. zsamelev_uv ) THEN
3763	iret = NF90_GET_VAR(force_id, sigma_uv_id, zhybriduv_b)
3764	zhybriduv_a = zero
3765	ENDIF
3766	ELSE
3767	ireta = NF90_GET_VAR(force_id, hybsigb_id, zhybrid_b)
3768	iretb = NF90_GET_VAR(force_id, hybsiga_id, zhybrid_a)
3769	IF ( ireta /= NF90_NOERR .OR. iretb /= NF90_NOERR) THEN
3770	CALL ipslerr ( 3, 'forcing_vertical','Attempted to read variable HybSigA and HybSigB from forcing file', &
3771	& 'NF90_GET_VAR failed.','stop')
3772	ENDIF
3773	IF ( .NOT. zsamelev_uv ) THEN
3774	ireta = NF90_GET_VAR(force_id, hybsigb_uv_id, zhybriduv_b)
3775	iretb = NF90_GET_VAR(force_id, hybsiga_uv_id, zhybriduv_a)
3776	IF ( ireta /= NF90_NOERR .OR. iretb /= NF90_NOERR) THEN
3777	CALL ipslerr ( 3, 'forcing_vertical',&
3778	& 'Attempted to read variable HybSigA_uv and HybSigB_uv from forcing file', &
3779	& 'NF90_GET_VAR failed.','stop')
3780	ENDIF
3781	ENDIF
3782	ENDIF
3783	WRITE(,) "forcing_vertical : case Pressure coordinates : "
3784	WRITE(,) "Read from forcing file :", zhybrid_b, zhybrid_a, zhybriduv_b, zhybriduv_a
3785	ELSE
3786	!
3787	! Why are we here ???
3788	!
3789	CALL ipslerr ( 3, 'forcing_vertical','What is the option used to describe the height of', &
3790	& 'the atmospheric forcing ?','Please check your forcing file.')
3791	ENDIF
3792	ENDIF
3793	!
3794	!- We have no forcing file to explore or we did not find anything. So revert back to the run.def and
3795	!- read what has been specified by the user.
3796	!
3797	IF ( force_id < 0 .OR. .NOT. foundvar ) THEN
3798	!
3799	!-
3800	!Config Key = HEIGHT_LEV1
3801	!Config Desc = Height at which T and Q are given
3802	!Config Def = 2.0
3803	!Config Help = The atmospheric variables (temperature and specific
3804	!Config humidity) are measured at a specific level.
3805	!Config The height of this level is needed to compute
3806	!Config correctly the turbulent transfer coefficients.
3807	!Config Look at the description of the forcing
3808	!Config DATA for the correct value.
3809	!-
3810	zlev_fixed = 2.0
3811	CALL getin('HEIGHT_LEV1', zlev_fixed)
3812	!-
3813	!Config Key = HEIGHT_LEVW
3814	!Config Desc = Height at which the wind is given
3815	!Config Def = 10.0
3816	!Config Help = The height at which wind is needed to compute
3817	!Config correctly the turbulent transfer coefficients.
3818	!-
3819	zlevuv_fixed = 10.0
3820	CALL getin('HEIGHT_LEVW', zlevuv_fixed)
3821
3822	zheight = .TRUE.
3823
3824	IF ( ABS(zlevuv_fixed-zlev_fixed) > EPSILON(zlev_fixed)) THEN
3825	zsamelev_uv = .FALSE.
3826	ENDIF
3827
3828	CALL ipslerr ( 2, 'forcing_vertical','The height of the atmospheric forcing variables', &
3829	& 'was not found in the netCDF file.','Thus the values in run.def were used ... or their defaults.')
3830	ENDIF
3831
3832	END SUBROUTINE forcing_vertical
3833
3834	!! =============================================================================================================================
3835	!! SUBROUTINE: forcing_givegrid
3836	!!
3837	!>\BRIEF Routine which exports the grid (longitude, latitude, land indices) on which the model will run, i.e. the zoomed grid.
3838	!!
3839	!! DESCRIPTION: This is needed to transfer the grid information from this module to the glogrid.f90 module.
3840	!!
3841	!!
3842	!! \n
3843	!_ ==============================================================================================================================
3844
3845	SUBROUTINE forcing_givegrid (lon, lat, mask, area, corners, lindex, contfrac, calendar_tmp)
3846	!
3847	! This subroutine will return to the caller the grid which has been extracted from the
3848	! the forcing file. It is assumed that the caller has called forcing_givegridsize before
3849	! and knows the dimensions of the fields and thus has done the correct allocations.
3850	!
3851	!
3852	REAL(r_std), INTENT(out) :: lon(iim_loc,jjm_loc), lat(iim_loc,jjm_loc)
3853	REAL(r_std), INTENT(out) :: mask(iim_loc,jjm_loc)
3854	REAL(r_std), INTENT(out) :: area(iim_loc,jjm_loc)
3855	REAL(r_std), INTENT(out) :: corners(iim_loc,jjm_loc,4,2)
3856	INTEGER(i_std), INTENT(out) :: lindex(nbpoint_loc)
3857	REAL(r_std), INTENT(out) :: contfrac(nbpoint_loc)
3858	CHARACTER(LEN=20), INTENT(out) :: calendar_tmp
3859	!
3860	IF ( .NOT. is_root_prc ) THEN
3861	CALL ipslerr (3,'forcing_givegrid'," This routine can only be called on the root processor.", &
3862	& "The information requested is only available on root processor.", " ")
3863	ENDIF
3864	!
3865	IF (nbpoint_loc .NE. nbland_loc) THEN
3866	WRITE(numout, *) "forcing_givegrid:: nbpoint_loc=", nbpoint_loc
3867	WRITE(numout, *) "forcing_givegrid:: nbland_loc=", nbland_loc
3868	CALL ipslerr(3,'forcing_givegrid','nbpoint_loc and nbland_loc do match', &
3869	'The calculation of land points is not correct','')
3870	ENDIF
3871	!
3872	lon(:,:) = lon_loc(:,:)
3873	lat(:,:) = lat_loc(:,:)
3874	!
3875	mask(:,:) = mask_loc(:,:)
3876	area(:,:) = area_loc(:,:)
3877	corners(:,:,:,:) = corners_loc(:,:,:,:)
3878	!
3879	!
3880	lindex(:) = lindex_loc(:)
3881	contfrac(:) = contfrac_loc(:)
3882	!
3883	calendar_tmp = calendar
3884	!
3885	END SUBROUTINE forcing_givegrid
3886
3887	!! =============================================================================================================================
3888	!! SUBROUTINE: forcing_checkdim
3889	!!
3890	!>\BRIEF
3891	!!
3892	!! DESCRIPTION: Save the dimension or check that it is equal to the previous value.
3893	!! Should one of the spatial dimensions be different between 2 files, then we have a big problem.
3894	!! They probably do not belong to the same set of forcing files.
3895	!!
3896	!! \n
3897	!_ ==============================================================================================================================
3898
3899	SUBROUTINE forcing_checkdim(ifile, filenames, out_dim, out_id, in_dim, in_id)
3900	!
3901	! Save the dimension or check that it is equal to the previous value.
3902	! Should one of the spatial dimensions be different between 2 files, then we have a big problem.
3903	! They probably do not belong to the same set of forcing files.
3904	!
3905	INTEGER(i_std), INTENT(in) :: ifile
3906	CHARACTER(LEN=*), INTENT(in) :: filenames(:)
3907	INTEGER(i_std), INTENT(out) :: out_dim, out_id
3908	INTEGER(i_std), INTENT(in) :: in_dim, in_id
3909	!
3910	IF ( ifile == 1 ) THEN
3911	out_dim = in_dim
3912	out_id = in_id
3913	ELSE
3914	IF ( out_dim /= in_dim ) THEN
3915	CALL ipslerr (3,'forcing_ocheckdim', 'The dimension of the file is not the same of the first file opened.', &
3916	& 'The offending file is : ', filenames(ifile))
3917	ENDIF
3918	ENDIF
3919	!
3920	END SUBROUTINE forcing_checkdim
3921
3922	!! =============================================================================================================================
3923	!! SUBROUTINE: forcing_time
3924	!!
3925	!>\BRIEF Read the time from each file and create the time axis to be the basis for the simulation.
3926	!!
3927	!! DESCRIPTION: This is an important routine which analyses the time axis of the forcing files and
3928	!! stores the main information in the SAVED variables of this routine.
3929	!! As this module manages a list of forcing files we also need to check that the time
3930	!! axis of all these files is continuous and homogeneous.
3931	!! The bounds are also build for all the time axes so that we know how to interpret the
3932	!! various variables.
3933	!!
3934	!! \n
3935	!_ ==============================================================================================================================
3936
3937	SUBROUTINE forcing_time(nbfiles, filenames)
3938	!
3939	! Read the time from each file and create the time axis to be the basis
3940	! for the simulation.
3941	!
3942	INTEGER(i_std) :: nbfiles
3943	CHARACTER(LEN=*) :: filenames(nbfiles)
3944	!
3945	INTEGER(i_std) :: iv, it, iff, tcnt, itbase, itbasetmp, ittmp
3946	INTEGER(i_std) :: tstart, maxtime_infile
3947	REAL(r_std), ALLOCATABLE, DIMENSION(:) :: timeint, time_read
3948	REAL(r_std), ALLOCATABLE, DIMENSION(:) :: time_infiles
3949	CHARACTER(LEN=20) :: axname, calendar, timevarname
3950	CHARACTER(LEN=60) :: timestamp, tmpatt, ymd, hms
3951	INTEGER(i_std) :: tncstart(3), tnccount(3)
3952	!
3953	INTEGER(i_std) :: iret, year0, month0, day0, hours0, minutes0, seci
3954	INTEGER(i_std), DIMENSION(1) :: imax, imin
3955	REAL(r_std) :: sec0, date_int, date0_tmp
3956	CHARACTER :: strc
3957	LOGICAL :: check=.FALSE.
3958	!
3959	! Check that we are working on the root proc.
3960	!
3961	IF ( .NOT. is_root_prc) THEN
3962	CALL ipslerr (3,'forcing_time',"Cannot run this routine o other procs than root.",&
3963	& "All the information on the forcing files is only lated on the root processor."," ")
3964	ENDIF
3965	!
3966	! Size of unlimited dimension added up through the files. If variable not allocated before by another
3967	! subroutine, it needs to be done here.
3968	!
3969	IF ( .NOT. ALLOCATED(nbtime_perfile) ) ALLOCATE(nbtime_perfile(nbfiles))
3970	IF ( .NOT. ALLOCATED(date0_file) ) ALLOCATE(date0_file(nbfiles,nbtax))
3971	!
3972	! Go through all files in the list in order to get the total number of time steps we have
3973	! in the nbfiles files to be read
3974	!
3975	nb_forcing_steps = 0
3976	maxtime_infile = 0
3977	DO iff=1,nbfiles
3978	!
3979	iret = NF90_INQUIRE_DIMENSION(force_id(iff), id_unlim(iff), name=axname, len=nbtime_perfile(iff))
3980	IF (iret /= NF90_NOERR) THEN
3981	CALL ipslerr (3,'forcing_time',"Could not get size of dimension of unlimited axis"," "," ")
3982	ENDIF
3983	nb_forcing_steps = nb_forcing_steps + nbtime_perfile(iff)
3984	IF ( nbtime_perfile(iff) > maxtime_infile ) maxtime_infile = nbtime_perfile(iff)
3985	ENDDO
3986	!
3987	! Allocate the variables needed with the time length just calculated.
3988	! These variables are saved in the module
3989	!
3990	ALLOCATE(time_infiles(nb_forcing_steps))
3991	ALLOCATE(time(nb_forcing_steps, nbtaxnbtmethods), time_bounds(nb_forcing_steps,nbtaxnbtmethods,2))
3992	ALLOCATE(time_axename(nbtaxnbtmethods), time_cellmethod(nbtaxnbtmethods))
3993	ALLOCATE(preciptime(nb_forcing_steps))
3994	ALLOCATE(time_sourcefile(nb_forcing_steps))
3995	ALLOCATE(time_id(nb_forcing_steps, nbtax))
3996	! Allocate local variables
3997	ALLOCATE(time_read(nb_forcing_steps))
3998	ALLOCATE(timeint(nb_forcing_steps))
3999	!
4000	! Get through all variables to find time_id
4001	! The key variables to filled up here are time (the time stamps read in the file) and
4002	! time_bounds which give the validity interval for the variables.
4003	!
4004	tstart=0
4005	!
4006	IF ( check ) WRITE(,) "forcing_time : going through ", nbfiles, " files to get the time."
4007	!
4008	DO iff=1,nbfiles
4009	!
4010	time_id(iff,:)=-1
4011	!
4012	! Go through the variables in the file and find the one which is a time axis.
4013	!
4014	tcnt=1
4015	DO iv=1,nvars(iff)
4016	iret = NF90_GET_ATT(force_id(iff), iv, "units", tmpatt)
4017	IF ( INDEX(lowercase(tmpatt),'seconds since') > 0) THEN
4018	time_id(iff,tcnt)=iv
4019	tcnt=tcnt+1
4020	convtosec(iff)=1.0
4021	ELSE IF ( INDEX(lowercase(tmpatt),'minutes since') > 0) THEN
4022	time_id(iff,tcnt)=iv
4023	tcnt=tcnt+1
4024	convtosec(iff)=60.0
4025	ELSE IF ( INDEX(lowercase(tmpatt),'hours since') > 0) THEN
4026	time_id(iff,tcnt)=iv
4027	tcnt=tcnt+1
4028	convtosec(iff)=3600.0
4029	ELSE IF ( INDEX(lowercase(tmpatt),'days since') > 0) THEN
4030	time_id(iff,tcnt)=iv
4031	tcnt=tcnt+1
4032	convtosec(iff)=one_day
4033	ENDIF
4034	ENDDO
4035	IF ( ANY(time_id(iff,:) < 0) ) THEN
4036	CALL ipslerr (3,'forcing_time',"Incorrect numer of time axes. A time axis is missing ",&
4037	& "in file :", filenames(iff))
4038	ENDIF
4039	!
4040	IF ( check ) WRITE(,) "forcing_time : Looking at time axis for file ", force_id(iff)
4041	!
4042	! Looping through the time axes and read them.
4043	!
4044	DO tcnt=1,nbtax
4045	!
4046	iret = NF90_INQUIRE_VARIABLE(force_id(iff), time_id(iff,tcnt), name=timevarname)
4047	IF ( check ) WRITE(,) "forcing_time : in ", iff, " found variable ", timevarname
4048	!
4049	! Get units of time axis
4050	!
4051	iret = NF90_GET_ATT(force_id(iff), time_id(iff,tcnt), "units", timestamp)
4052	IF ( check ) WRITE(,) "forcing_time : has time stamp ", timestamp
4053	!
4054	! Transform the start date of the netCDF file into a julian date for the model
4055	!
4056	timestamp = TRIM(timestamp(INDEX(timestamp,'since')+6:LEN_TRIM(timestamp)))
4057	ymd=TRIM(timestamp(1:INDEX(timestamp,' ')))
4058	hms=TRIM(timestamp(INDEX(timestamp,' ')+1:LEN_TRIM(timestamp)))
4059	!
4060	! First extral the information from the date string
4061	!
4062	READ(ymd(1:INDEX(ymd,'-')-1),'(I4)') year0
4063	ymd=TRIM(ymd(INDEX(ymd,'-')+1:LEN_TRIM(ymd)))
4064	READ(ymd(1:INDEX(ymd,'-')-1),'(I2)') month0
4065	ymd=TRIM(ymd(INDEX(ymd,'-')+1:LEN_TRIM(ymd)))
4066	READ(ymd,'(I2)') day0
4067	!
4068	! Now extract from the time string
4069	!
4070	READ(hms(1:INDEX(hms,':')-1),'(I2)') hours0
4071	hms=TRIM(hms(INDEX(hms,':')+1:LEN_TRIM(hms)))
4072	READ(hms(1:INDEX(hms,':')-1),'(I2)') minutes0
4073	hms=TRIM(hms(INDEX(hms,':')+1:LEN_TRIM(hms)))
4074	READ(hms,'(I2)') seci
4075	!
4076	sec0 = hours03600. + minutes060. + seci
4077	CALL ymds2ju (year0, month0, day0, sec0, date0_tmp)
4078	date0_file(iff,tcnt) = date0_tmp
4079	!
4080	! Now get the actual dates
4081	!
4082	tncstart(1) = 1
4083	tnccount(1) = nbtime_perfile(iff)
4084	IF ( check ) WRITE(,) "forcing_time : number of values read : ", tnccount(1)
4085	iret = NF90_GET_VAR(force_id(iff), time_id(iff,tcnt), time_read, tncstart, tnccount)
4086	IF (iret /= NF90_NOERR) THEN
4087	CALL ipslerr (3,'forcing_time',"An error occured while reading time from the file."," "," ")
4088	ENDIF
4089	!
4090	! Convert the variable time from seconds since to julian days
4091	!
4092	DO it=1,nbtime_perfile(iff)
4093	!!time_infiles(tstart+it) = date0_file(iff,tcnt) + time_read(it)*convtosec(iff)/one_day
4094	IF ( convtosec(iff) < one_day ) THEN
4095	time_infiles(tstart+it) = date0_file(iff,tcnt) + time_read(it)*convtosec(iff)/one_day
4096	ELSE
4097	! In the case of daily forcing the start time has to be 00UTC in Julian days.
4098	time_infiles(tstart+it) = date0_file(iff,tcnt) + INT(time_read(it))
4099	ENDIF
4100	ENDDO
4101	if ( check ) WRITE(,) "File ", iff, "goes from ", time_infiles(tstart+1), " to ", &
4102	time_infiles(tstart+nbtime_perfile(iff))
4103	!
4104	! Estimate the bounds as this information is not yet in the forcing file.
4105	!
4106	date_int = (time_infiles(tstart+nbtime_perfile(iff)) - time_infiles(tstart+1))/(nbtime_perfile(iff)-1)
4107	forcing_tstep_ave = date_int*one_day
4108	!
4109	! If this is the first file then simply keep the name of the time axis. Else search the same name
4110	! in what has already been read
4111	!
4112	IF ( iff == 1 ) THEN
4113	itbase = (tcnt-1)*nbtmethods
4114	time_axename(itbase+1:itbase+4) = timevarname
4115	time_cellmethod(itbase+1) = "reference"
4116	time_cellmethod(itbase+2) = "start"
4117	time_cellmethod(itbase+3) = "cent"
4118	time_cellmethod(itbase+4) = "end"
4119	ELSE
4120	!
4121	! If this is not the first file then we need to find the correct axis to add to.
4122	! All information have already been saved with the first file.
4123	!
4124	DO ittmp=1,nbtax
4125	itbasetmp=(ittmp-1)*nbtmethods
4126	IF ( time_axename(itbasetmp+1) == timevarname ) THEN
4127	itbase = itbasetmp
4128	ENDIF
4129	ENDDO
4130
4131	ENDIF
4132	!
4133	!
4134	! Keep for future usage the various information on the time axis we have just read. This time axis can
4135	! be understood in 3 different ways and each variable might use a different cell method for this time
4136	! axis.
4137	!
4138	! time(:,(tcnt-1)*nbtmethods+1) : corresponds to the reference time axis as it has been read from the file
4139	! time(:,(tcnt-1)*nbtmethods+2) : is the time axis with a cell method which place the value at the
4140	! beginning of the time interval
4141	! time(:,(tcnt-1)*nbtmethods+3) : is the time axis corresponding to variables placed at the center of the
4142	! time interval
4143	! time(:,(tcnt-1)*nbtmethods+4) : for variables put at the end of the time interval over which they aere
4144	! for instance averaged.
4145	!
4146	! In variable time_cellmethod we will write the type of cell method as descirbed above so that the selection
4147	! of the right axis for each variable can be made automaticaly.
4148	!
4149	! We also keep the name of the time axis read in preparation of file where we might have to read more than one
4150	! time axis.
4151	!
4152	DO it=tstart+1,tstart+nbtime_perfile(iff)
4153	!
4154	! Reference time
4155	!
4156	time(it,itbase+1) = time_infiles(it)
4157	time_bounds(it,itbase+1,1) = time_infiles(it)-date_int/2.0
4158	time_bounds(it,itbase+1,2) = time_infiles(it)+date_int/2.0
4159	!
4160	! Start cell method
4161	time(it,itbase+2) = time_infiles(it)+date_int/2.0
4162	time_bounds(it,itbase+2,1) = time_infiles(it)
4163	time_bounds(it,itbase+2,2) = time_infiles(it)+date_int
4164	!
4165	! Centered cell method
4166	time(it,itbase+3) = time_infiles(it)
4167	time_bounds(it,itbase+3,1) = time_infiles(it)-date_int/2.0
4168	time_bounds(it,itbase+3,2) = time_infiles(it)+date_int/2.0
4169	!
4170	! End cell method
4171	time(it,itbase+4) = time_infiles(it)-date_int/2.0
4172	time_bounds(it,itbase+4,1) = time_infiles(it)-date_int
4173	time_bounds(it,itbase+4,2) = time_infiles(it)
4174	!
4175	ENDDO
4176	!
4177	! Keep the number of the file from which we read this time.
4178	!
4179	time_sourcefile(tstart+1:tstart+nbtime_perfile(iff))=iff
4180	!
4181	IF ( check ) WRITE(,) "forcing_time : finished file ", iff
4182	!
4183	ENDDO
4184	!
4185	! Before moving to the next file advance the pointer in the time arrays.
4186	!
4187	tstart=tstart+nbtime_perfile(iff)
4188	!
4189	ENDDO
4190	!
4191	IF ( check ) WRITE(,) "forcing_time : All files have been processed"
4192	!
4193	! Verify that the forcing comes in regular time intervals. If not, many of the
4194	! interpolation schemes will fail.
4195	! This is only done on the first time axis ... is it enough ?
4196	!
4197	DO ittmp=1,nbtax
4198	itbase=(ittmp-1)*nbtmethods
4199	!
4200	date_int = (time(nb_forcing_steps,itbase+1) - time(1,itbase+1))/(nb_forcing_steps-1)
4201	forcing_tstep_ave = date_int*one_day
4202	!
4203	timeint(:) = 0
4204	DO it=1, nb_forcing_steps-1
4205	timeint(it) = time(it+1,itbase+1)-time(it,itbase+1)
4206	ENDDO
4207	!
4208	IF ( MAXVAL(timeint(1:nb_forcing_steps-1)) > date_int+0.1*date_int .OR.&
4209	& MINVAL(timeint(1:nb_forcing_steps-1)) < date_int-0.1*date_int) THEN
4210	WRITE(,) "The time steping of the forcing files does not seem to be regular on axis",time_axename(itbase+1),":"
4211	WRITE(,) "Average time step : ", date_int, "days = ", date_int*one_day, "sec."
4212	imax = MAXLOC(timeint(1:nb_forcing_steps-1))
4213	imin = MINLOC(timeint(1:nb_forcing_steps-1))
4214	WRITE(,) "Maximum time step : ", MAXVAL(timeint(1:nb_forcing_steps-1)), " at ", imax(1)
4215	WRITE(,) "Minimum time step : ", MINVAL(timeint(1:nb_forcing_steps-1)), " at ", imin(1)
4216	WRITE(,) "++++ Values around Maximum"
4217	DO it=MAX(imax(1)-5,1),MIN(imax(1)+5,nb_forcing_steps)
4218	WRITE(,) it, " from file ", time_sourcefile(it), " Value ", time(it,itbase+1)
4219	CALL forcing_printdate(time(it,itbase+1), "Time values.")
4220	ENDDO
4221	WRITE(,) "++++ Values around Minimum"
4222	DO it=MAX(imin(1)-5,1),MIN(imin(1)+5,nb_forcing_steps)
4223	WRITE(,) it, " from file ", time_sourcefile(it), " Value ", time(it,itbase+1)
4224	CALL forcing_printdate(time(it,itbase+1), "Time values.")
4225	ENDDO
4226	CALL ipslerr (3,'forcing_time', 'The time handling could be improved to allow the driver',&
4227	& "to cope with irregular forcing files."," ")
4228	ENDIF
4229	ENDDO
4230	!
4231	! Print some test values
4232	!
4233	DO ittmp=1,nbtax
4234	itbase=(ittmp-1)*nbtmethods
4235	!
4236	WRITE(,) "Bounds for axis ",time_axename(itbase+1)," :"
4237	!
4238	CALL forcing_printdate(time_bounds(1,itbase+1,1), "Start time of first forcing interval.")
4239	CALL forcing_printdate(time_bounds(1,itbase+1,2), "End time of first forcing interval.")
4240	CALL forcing_printdate(time_bounds(nb_forcing_steps,itbase+1,1), "Start time of last forcing interval.")
4241	CALL forcing_printdate(time_bounds(nb_forcing_steps,itbase+1,2), "End time of last forcing interval.")
4242	ENDDO
4243	!
4244	! Set to zero the variable for the cummulated time for rainfall
4245	!
4246	preciptime(:) = zero
4247	!
4248	! Keep the very first date of the time axis for future use
4249	!
4250	forcingstartdate = time(1,1)
4251	!
4252	! Clean-up
4253	!
4254	DEALLOCATE(timeint, time_read)
4255	!
4256	END SUBROUTINE forcing_time
4257
4258	!! =============================================================================================================================
4259	!! SUBROUTINE: forcing_varforslab
4260	!!
4261	!>\BRIEF
4262	!!
4263	!! DESCRIPTION: 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	!! \n
4267	!_ ==============================================================================================================================
4268
4269	SUBROUTINE forcing_varforslab(fileindex, varname, timestart, timecount, inslabpos, data, cellmethod)
4270	!
4271	! This subroutine will read the named variable and put it in the right place in the
4272	! slab of data kept in the memory of the driver.
4273	!
4274	INTEGER(i_std), INTENT(in) :: fileindex
4275	CHARACTER(LEN=*), INTENT(in) :: varname
4276	INTEGER(i_std), INTENT(in) :: timestart, timecount, inslabpos
4277	REAL(r_std), INTENT(inout) :: data(nbpoint_loc,slab_size)
4278	CHARACTER(LEN=*), INTENT(out) :: cellmethod
4279	!
4280	! Local variables
4281	!
4282	INTEGER(i_std) :: varid, windid, windndims, it, ig, iv
4283	INTEGER(i_std) :: iret, ndims
4284	INTEGER(i_std), DIMENSION(4) :: start, count
4285	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: tmp_slab
4286	REAL(r_std), ALLOCATABLE, DIMENSION(:,:,:) :: tmp_slab2d
4287	CHARACTER(LEN=80) :: name
4288	LOGICAL :: windzero
4289	!
4290	! Allocate the temporary data array if not already available
4291	!
4292	IF ( compressed ) THEN
4293	IF ( .NOT. ALLOCATED(tmp_slab) ) ALLOCATE(tmp_slab(ncdfcount,slab_size))
4294	ELSE
4295	IF ( .NOT. ALLOCATED(tmp_slab2d) ) ALLOCATE(tmp_slab2d(iim_glo,jjm_glo,slab_size))
4296	ENDIF
4297	!
4298	! Reset the counters and flags to forget the past !
4299	!
4300	varid=-1
4301	windid=-1
4302	windzero=.FALSE.
4303	!
4304	! Find the variable in the file
4305	!
4306	DO iv=1,nvars(fileindex)
4307	!
4308	iret = NF90_INQUIRE_VARIABLE(force_id(fileindex), iv, name=name, ndims=it)
4309	!
4310	IF ( INDEX(name, varname) > 0 ) THEN
4311	varid = iv
4312	ndims = it
4313	ENDIF
4314	IF ( (INDEX(name, "Wind") > 0) .AND. (LEN_TRIM(name) == LEN_TRIM("Wind")) ) THEN
4315	windid = iv
4316	windndims = it
4317	ENDIF
4318	!
4319	ENDDO
4320	!
4321	! Treat some special cases and catch errors
4322	!
4323	IF ( varid < 0 ) THEN
4324	!
4325	! If we requested a wind component but did not find it, it could be that only the module is available.
4326	! If that is the case, then we use the module (windid) for the U component and set V top zero.
4327	!
4328	IF ( INDEX(varname, "Wind_E") > 0 ) THEN
4329	varid = windid
4330	ndims = windndims
4331	windzero = .FALSE.
4332	ELSE IF ( INDEX(varname, "Wind_N") > 0 ) THEN
4333	windzero = .TRUE.
4334	ELSE
4335	CALL ipslerr (3,'forcing_varforslab',"Could not find variable",varname," in file.")
4336	ENDIF
4337	ENDIF
4338	!
4339	! If there is some variable to be read then do it
4340	!
4341	IF ( .NOT. windzero ) THEN
4342	!
4343	! Get the attributes needed for intepretating the data
4344	!
4345	! First get the cell method used for this variable
4346	iret = NF90_GET_ATT(force_id(fileindex), varid, 'cell_methods', cellmethod)
4347	IF (iret /= NF90_NOERR) THEN
4348	! If the attribute is not found then we set a reasonable default : instantaneous and centered.
4349	cellmethod="time: instantaneous"
4350	ENDIF
4351	!
4352	!
4353	! Getsize of data to be read from the netCDF file
4354	!
4355	!
4356	IF ( compressed ) THEN
4357	!
4358	IF ( ndims == 2 ) THEN
4359	start = (/ncdfstart,timestart,0,0/)
4360	count = (/ncdfcount,timecount,0,0/)
4361	ELSE IF ( ndims == 3 ) THEN
4362	start = (/ncdfstart,1,timestart,0/)
4363	count = (/ncdfcount,1,timecount,0/)
4364	ELSE
4365	CALL ipslerr (3,'forcing_varforslab',"Compressed variable : ",varname,&
4366	& " does not have a compatible number of dimensions.")
4367	ENDIF
4368	!
4369	iret = NF90_GET_VAR(force_id(fileindex), varid, tmp_slab, start, count)
4370	IF (iret /= NF90_NOERR) THEN
4371	CALL ipslerr (3,'forcing_varforslab',"Could not read from file variable : ",varname," Compressed in the file.")
4372	ENDIF
4373	!
4374	! Zoom into the data and put it in the right place in the slab of data.
4375	!
4376	CALL forcing_reindex(ncdfcount, timecount, tmp_slab, nbpoint_loc, slab_size, data, inslabpos, reindex_loc)
4377	ELSE
4378	!
4379	IF ( ndims == 3 ) THEN
4380	start = (/1,1,timestart,0/)
4381	count = (/iim_glo,jjm_glo,timecount,0/)
4382	ELSE IF (ndims == 4 ) THEN
4383	start = (/1,1,1,timestart/)
4384	count = (/iim_glo,jjm_glo,1,timecount/)
4385	ELSE
4386	CALL ipslerr (3,'forcing_varforslab',"Full lat Lon variable : ",varname,&
4387	& " does not have a compatible number of dimensions.")
4388	ENDIF
4389	!
4390	iret = NF90_GET_VAR(force_id(fileindex), varid, tmp_slab2d, start, count)
4391	IF (iret /= NF90_NOERR) THEN
4392	WRITE(,) TRIM(NF90_STRERROR(iret))
4393	WRITE(,) "File =", fileindex, "Size =", SIZE(tmp_slab2d,DIM=1), SIZE(tmp_slab2d,DIM=2), SIZE(tmp_slab2d,DIM=3)
4394	WRITE(,) "Start :", start(1:3)
4395	WRITE(,) "Count :", count(1:3)
4396	CALL ipslerr (3,'forcing_varforslab',"Could not read from file variable : ",varname," Not compressed.")
4397	ENDIF
4398	!
4399	! Zoom into the data and put it in the right place in the slab of data.
4400	!
4401
4402	CALL forcing_reindex(iim_glo, jjm_glo, timecount, tmp_slab2d, nbpoint_loc, slab_size, data, inslabpos, reindex2d_loc)
4403
4404	ENDIF
4405	ELSE
4406	cellmethod="time: instantaneous"
4407	DO it=0,timecount-1
4408	data(:,inslabpos+it) = zero
4409	ENDDO
4410	ENDIF
4411	!
4412	END SUBROUTINE forcing_varforslab
4413
4414	!! =============================================================================================================================
4415	!! SUBROUTINE: forcing_attributetimeaxe
4416	!!
4417	!>\BRIEF Find the time axis which corresponds to the variable at hand.
4418	!!
4419	!! DESCRIPTION: We interpret the cell_method provided in the netCDF file so that
4420	!! we can determine how we need to understand the values we read.
4421	!!
4422	!! \n
4423	!_ ==============================================================================================================================
4424
4425	SUBROUTINE forcing_attributetimeaxe(cellmethod, timeindex)
4426	!
4427	! We will analyse the time axis of the cell method found in the NetCDF file in order to
4428	! attribute the correct time axis to this variable.
4429	!
4430	CHARACTER(LEN=*), INTENT(in) :: cellmethod
4431	INTEGER(i_std), INTENT(out) :: timeindex
4432	!
4433	INTEGER(i_std) :: itax, timepos, pos, lentime, itbase, im
4434	CHARACTER(LEN=20) :: TARGET, tmpstr
4435	CHARACTER(LEN=80) :: method
4436	!
4437	! Clean the string to delete spaces in front of ":" and "("
4438	!
4439	method = cellmethod
4440	DO WHILE ( INDEX(method," :") > 0 )
4441	pos = INDEX(method," :")
4442	method = method(1:pos-1)//method(pos+1:LEN_TRIM(method))
4443	ENDDO
4444	DO WHILE ( INDEX(method,"( ") > 0 )
4445	pos = INDEX(method,"( ")
4446	method = method(1:pos)//method(pos+2:LEN_TRIM(method))
4447	ENDDO
4448	!
4449	! Go through all the time axes we have to find the right one.
4450	!
4451	timeindex=0
4452	DO itax=1,nbtax
4453	!
4454	itbase=(itax-1)*nbtmethods
4455	! find the time axis name in the cell method
4456	TARGET = TRIM(time_axename(itbase+1))//":"
4457	timepos = INDEX(method,TRIM(TARGET))
4458	!
4459	! If we found the time axis then we look for the method with a time position description
4460	! which is expected to be between parenthesis. For instance : mean(end)
4461	!
4462	IF ( timepos > 0 ) THEN
4463	!
4464	lentime=LEN_TRIM(time_axename(itbase+1))
4465	tmpstr = method(lentime+2:LEN_TRIM(method))
4466	!
4467	! If there is ":" then there is information for another axis which needs to be deleted
4468	!
4469	IF ( INDEX(tmpstr,":") > 0 ) THEN
4470	tmpstr = tmpstr(1:INDEX(tmpstr,":")-1)
4471	ENDIF
4472	!
4473	! Now that we have found a time axis see if we have between parenthesis a position
4474	! on that time avis.
4475	!
4476	! Look for a "("
4477	IF ( INDEX(tmpstr, "(") > 0 ) THEN
4478	DO im=1,nbtmethods
4479	TARGET = "("//TRIM(time_cellmethod(itbase+im))
4480	timepos = INDEX(tmpstr,TRIM(TARGET))
4481	IF ( timepos > 0 ) THEN
4482	timeindex = itbase+im
4483	ENDIF
4484	ENDDO
4485	!
4486	! If there is no "(" then we have to find the centered axis.
4487	ELSE
4488	DO im=1,nbtmethods
4489	IF ( INDEX(time_cellmethod(itbase+im), "cent") > 0 ) THEN
4490	timeindex = itbase+im
4491	ENDIF
4492	ENDDO
4493	ENDIF
4494	!
4495	! The name of the time axis was found bu no method could be identified
4496	!
4497	IF ( timeindex < 1 ) THEN
4498	CALL ipslerr (3,'forcing_attributetimeaxe',"Found a time axis name but could not identify method.", &
4499	"This should not happen !"," ")
4500	ENDIF
4501	!
4502	ELSE
4503	! Continue in loop over nbtax
4504	ENDIF
4505	ENDDO
4506	!
4507	! Should no corresponding time axis name be found,
4508	! then we use the first centered one.
4509	!
4510	itax=1
4511	DO WHILE ( timeindex < 1 )
4512	IF ( INDEX(time_cellmethod(itax), "cent") > 0 ) THEN
4513	timeindex = itax
4514	ELSE
4515	itax = itax + 1
4516	ENDIF
4517	ENDDO
4518	!
4519	END SUBROUTINE forcing_attributetimeaxe
4520
4521	!! =============================================================================================================================
4522	!! SUBROUTINE: forcing_filenamecheck
4523	!!
4524	!>\BRIEF Check the list of files obtained from the calling program.
4525	!!
4526	!! DESCRIPTION: A small code to check the forcing files. They have to be NetCDF (i.e. .nc termination) and
4527	!! we dont keep files that appear more than once in the list.
4528	!!
4529	!! \n
4530	!_ ==============================================================================================================================
4531
4532	SUBROUTINE forcing_filenamecheck(filenames_in, nb_files)
4533	!
4534	! A small code to check the forcing files. They have to
4535	! be NetCDF (i.e. .nc termination) and we dont keep files
4536	! that appear more than once in the list.
4537	!
4538	!
4539	! INPUT
4540	!
4541	CHARACTER(LEN=*), DIMENSION(:), INTENT(in) :: filenames_in
4542	INTEGER(i_std), INTENT(out) :: nb_files
4543	!
4544	! LOCAL
4545	!
4546	INTEGER(i_std) :: ii, is, sizein
4547	LOGICAL :: notfound
4548	!
4549	sizein = SIZE(filenames_in)
4550	IF ( sizein > 0 ) THEN
4551	IF ( ALLOCATED(forfilename) ) THEN
4552	DEALLOCATE(forfilename)
4553	ENDIF
4554	ALLOCATE(forfilename(sizein))
4555	nb_files=0
4556	ELSE
4557	CALL ipslerr (3,'forcing_filenamecheck',"List of forcing files is empty.","Please check your run.def file."," ")
4558	ENDIF
4559	!
4560	DO ii=1,sizein
4561	IF ( INDEX(filenames_in(ii), '.nc') > 0 ) THEN
4562	IF ( nb_files == 0 ) THEN
4563	nb_files = nb_files+1
4564	forfilename(nb_files)=TRIM(ADJUSTL(filenames_in(ii)))
4565	ELSE
4566	notfound=.TRUE.
4567	DO is=1,nb_files
4568	IF ( INDEX(TRIM(filenames_in(ii)), TRIM(ADJUSTL(forfilename(is)))) > 0 ) notfound=.FALSE.
4569	ENDDO
4570	IF ( notfound ) THEN
4571	nb_files = nb_files+1
4572	forfilename(nb_files)=TRIM(adjustl(filenames_in(ii)))
4573	ENDIF
4574	ENDIF
4575	ELSE
4576	!!! This is not a NetCDF file, so we ignore it
4577	ENDIF
4578	ENDDO
4579	!
4580	!
4581	END SUBROUTINE forcing_filenamecheck
4582
4583	!! =============================================================================================================================
4584	!! FUNCTION: lowercase, FindMinimum, Swap
4585	!!
4586	!>\BRIEF Help functions for the forcing_tools module.
4587	!!
4588	!! DESCRIPTION:
4589	!!
4590	!! \n
4591	!_ ==============================================================================================================================
4592
4593	FUNCTION lowercase(strIn) RESULT(strOut)
4594	! Adapted from http://www.star.le.ac.uk/~cgp/fortran.html (25 May 2012)
4595
4596	IMPLICIT NONE
4597
4598	CHARACTER(len=*), INTENT(in) :: strIn
4599	CHARACTER(len=LEN(strIn)) :: strOut
4600	INTEGER :: i,j
4601
4602	DO i = 1, LEN(strIn)
4603	j = IACHAR(strIn(i:i))
4604	IF (j>= IACHAR("A") .AND. j<=IACHAR("Z") ) THEN
4605	strOut(i:i) = ACHAR(IACHAR(strIn(i:i))+32)
4606	ELSE
4607	strOut(i:i) = strIn(i:i)
4608	END IF
4609	END DO
4610
4611	END FUNCTION lowercase
4612	!
4613	! Some help function found on Internet : http://www.cs.mtu.edu/~shene/COURSES/cs201/NOTES/chap08/sorting.f90
4614	!
4615	! --------------------------------------------------------------------
4616	! INTEGER FUNCTION FindMinimum():
4617	! This function returns the location of the minimum in the section
4618	! between Start and End.
4619	! --------------------------------------------------------------------
4620	INTEGER FUNCTION FindMinimum(x, Start, END)
4621	IMPLICIT NONE
4622	REAL(r_std), DIMENSION(1:), INTENT(IN) :: x
4623	INTEGER(i_std), INTENT(IN) :: Start, END
4624	REAL(r_std) :: Minimum
4625	INTEGER(i_std) :: Location
4626	INTEGER(i_std) :: i
4627
4628	Minimum = x(Start) ! assume the first is the min
4629	Location = Start ! record its position
4630	DO i = Start+1, END ! start with next elements
4631	IF (x(i) < Minimum) THEN ! if x(i) less than the min?
4632	Minimum = x(i) ! Yes, a new minimum found
4633	Location = i ! record its position
4634	ENDIF
4635	END DO
4636	FindMinimum = Location ! return the position
4637	END FUNCTION FindMinimum
4638	! --------------------------------------------------------------------
4639	! SUBROUTINE Swap():
4640	! This subroutine swaps the values of its two formal arguments.
4641	! --------------------------------------------------------------------
4642	SUBROUTINE Swap(a, b)
4643	IMPLICIT NONE
4644	REAL(r_std), INTENT(INOUT) :: a, b
4645	REAL(r_std) :: Temp
4646
4647	Temp = a
4648	a = b
4649	b = Temp
4650	END SUBROUTINE Swap
4651	! --------------------------------------------------------------------
4652	! SUBROUTINE Sort():
4653	! This subroutine receives an array x() and sorts it into ascending
4654	! order.
4655	! --------------------------------------------------------------------
4656	SUBROUTINE Sort(x, Size)
4657	IMPLICIT NONE
4658	REAL(r_std), DIMENSION(1:), INTENT(INOUT) :: x
4659	INTEGER(i_std), INTENT(IN) :: Size
4660	INTEGER(i_std) :: i
4661	INTEGER(i_std) :: Location
4662
4663	DO i = 1, Size-1 ! except for the last
4664	Location = FindMinimum(x, i, Size) ! find min from this to last
4665	CALL Swap(x(i), x(Location)) ! swap this and the minimum
4666	END DO
4667	END SUBROUTINE Sort
4668
4669	END MODULE forcing_tools

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

Context Navigation

source: branches/publications/ORCHIDEE_2.2_r7266/ORCHIDEE/src_driver/forcing_tools.f90 @ 7541

Download in other formats: