Context Navigation

source: branches/publications/ORCHIDEE_gmd-2018-57/src_driver/forcing_tools.f90 @ 5143

Last change on this file since 5143 was 4433, checked in by jan.polcher, 7 years ago
Small corrections noted when doing tests on ADA
File size: 164.1 KB

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

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

Download in other formats: