Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

forcing_tools.f90

Last change on this file was 6939, checked in by jinfeng.chang, 4 years ago
Update ORCHIDEE-GMv3.2 for publication
File size: 175.1 KB

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

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

Context Navigation

source: branches/publications/ORCHIDEE-GMv3.2/ORCHIDEE/src_driver/forcing_tools.f90

Download in other formats: