Context Navigation

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

Last change on this file since 7328 was 7328, checked in by josefine.ghattas, 3 years ago
Integrated commit [7315] from the trunk: Variables should always be allocated when passed as argument. Otherwise strict debug options might fail. This was the case with the case at jean-zay.
File size: 185.7 KB

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

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

Download in other formats: