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

Last change on this file since 7330 was 7330, checked in by josefine.ghattas, 3 months ago

Corrected error introduced in previous commit [7328]

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