Context Navigation

forcing_tools.f90 @ 6593

Last change on this file since 6593 was 5599, checked in by josefine.ghattas, 6 years ago

Integrated changeset [5453] and [5598] done on branche ORCCHIDEE-ROUTING for interpolation of daily forcing files with orchideedriver:

A module was added to the new driver in order to generate diurnal cycles.
This allows to run with ISIMIP forcing data.
The forcing is expected to be in the same format qs 3h forcing, but two additional variables are expected : Tairmin and Tairmax.

Jan Polcher

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

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

Context Navigation

source: tags/ORCHIDEE_2_1/ORCHIDEE/src_driver/forcing_tools.f90 @ 6593

Download in other formats: