Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

forcing_tools.f90

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

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

Context Navigation

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

Download in other formats: