Context Navigation

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

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

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

Download in other formats: