source: branches/ORCHIDEE_2_2/ORCHIDEE/src_driver/forcing_tools.f90 @ 7258

Last change on this file since 7258 was 7258, checked in by agnes.ducharne, 6 months ago

As done in r6136 for SPREAD_PREC_CONT, and in r6326 for interpolation of regional maps.

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