Context Navigation

source: branches/publications/ORCHIDEE_DFv1.0_site/src_driver/forcing_tools.f90 @ 6715

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

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

Download in other formats: