Context Navigation

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

Last change on this file since 7429 was 7429, checked in by josefine.ghattas, 2 years ago

Correction to run with the new driver using compressed forcing data (tested with crujra half deg). Already done in the trunk rev [7281].

Done by Xiaoni Wang-Faivre

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

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

Download in other formats: