Context Navigation

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

Last change on this file since 7259 was 7258, checked in by agnes.ducharne, 3 years ago
As done in r6136 for SPREAD_PREC_CONT, and in r6326 for interpolation of regional maps.
File size: 180.4 KB

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

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

Download in other formats: