Context Navigation

forcing_tools.f90

Last change on this file was 4646, checked in by josefine.ghattas, 7 years ago

Ticket #185 :
Created new module containing time variables and subroutines to calculate them. Variables in the time module are public and can be used elswhere in the model. They should not be modified outside the time module.

Note that now each time step have a time interval. xx_start or xx_end values for the date can be used. Previously the date values correspondend to the end of the interval. See description in module time for more information.

These changes do not make any difference in results for couled mode with LMDZ or dim2_driver offline driver.

These modifications also corrects the problem with do_slow related to orchideedriver, ticket #327

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

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

Context Navigation

source: branches/publications/ORCHIDEE_2.0_gmd_Africa/ORCHIDEE/src_driver/forcing_tools.f90

Download in other formats: