Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

forcing_tools.f90

Last change on this file was 8606, checked in by josefine.ghattas, 3 days ago
Correction on commit [8595]. Unvisible bad caracteres where induced. Ticket #836
File size: 195.1 KB

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