Context Navigation

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

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

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

Download in other formats: