Context Navigation

source: branches/publications/ORCHIDEE_GLUC_r6545/src_driver/forcing_tools.f90 @ 6737

Last change on this file since 6737 was 4719, checked in by albert.jornet, 7 years ago

Merge: from revisions [4491:4695/trunk/ORCHIDEE]

Merge done in [4671:4718/perso/albert.jornet/MICT_MERGE]

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

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

Download in other formats: