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

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

Integrated commit [7314] from the trunk: Corrected error seen in debug mode with new driver: Treated case when variable CYCLIC_STARTDATE and CYCLIC_ENDDATE were not present in run.def.

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          ALLOCATE(hurs_full(1,1))
1495       ENDIF
1496       ALLOCATE(qair_full(nbpoint_loc,slab_size))
1497       ALLOCATE(rainf_full(nbpoint_loc,slab_size))
1498       ALLOCATE(snowf_full(nbpoint_loc,slab_size))
1499       ALLOCATE(swdown_full(nbpoint_loc,slab_size))
1500       ALLOCATE(lwdown_full(nbpoint_loc,slab_size))
1501       ALLOCATE(u_full(nbpoint_loc,slab_size))
1502       ALLOCATE(v_full(nbpoint_loc,slab_size))
1503       ALLOCATE(ps_full(nbpoint_loc,slab_size))
1504       ALLOCATE(ztq_full(nbpoint_loc,slab_size))
1505       ALLOCATE(zuv_full(nbpoint_loc,slab_size))
1506       !
1507       CALL forcing_readslab_root(time_int, &
1508            &                     tair_full, tairmax_full, tairmin_full, time_tair, timebnd_tair, &
1509            &                     qair_full, time_qair, timebnd_qair, &
1510            &                     rainf_full, snowf_full, time_precip, timebnd_precip, preciptime_tmp, &
1511            &                     swdown_full, time_swdown, timebnd_swdown, &
1512            &                     lwdown_full, time_lwdown, timebnd_lwdown, &
1513            &                     u_full, time_u, timebnd_u, &
1514            &                     v_full, time_v, timebnd_v, &
1515            &                     ps_full, time_ps, timebnd_ps, &
1516            &                     ztq_full, zuv_full)
1517       !
1518    ELSE
1519       !
1520       ALLOCATE(tair_full(1,1))
1521       IF ( forcing_tstep_ave >= one_day/3.0) THEN
1522          ALLOCATE(tairmax_full(1,1))
1523          ALLOCATE(tairmin_full(1,1))
1524       ENDIF
1525       ALLOCATE(qair_full(1,1))
1526       ALLOCATE(rainf_full(1,1))
1527       ALLOCATE(snowf_full(1,1))
1528       ALLOCATE(swdown_full(1,1))
1529       ALLOCATE(lwdown_full(1,1))
1530       ALLOCATE(u_full(1,1))
1531       ALLOCATE(v_full(1,1))
1532       ALLOCATE(ps_full(1,1))
1533       ALLOCATE(ztq_full(1,1))
1534       ALLOCATE(zuv_full(1,1))
1535       !
1536    ENDIF
1537    !
1538    ! Broadcast the time information to all procs.
1539    !
1540    CALL bcast(slab_size)
1541    CALL bcast(time_tair)
1542    CALL bcast(timebnd_tair)
1543    CALL bcast(time_qair)
1544    CALL bcast(timebnd_qair)
1545    CALL bcast(time_precip)
1546    CALL bcast(timebnd_precip)
1547    CALL bcast(preciptime_tmp)
1548    CALL bcast(time_swdown)
1549    CALL bcast(timebnd_swdown)
1550    CALL bcast(time_lwdown)
1551    CALL bcast(timebnd_lwdown)
1552    CALL bcast(time_u)
1553    CALL bcast(timebnd_u)
1554    CALL bcast(time_v)
1555    CALL bcast(timebnd_v)
1556    CALL bcast(time_ps)
1557    CALL bcast(timebnd_ps)
1558    !
1559    ! Scatter the slabs of data to all processors
1560    !
1561    IF ( landonly ) THEN
1562       CALL scatter(tair_full, tair_slab)
1563       IF ( forcing_tstep_ave >= one_day/3.0) THEN
1564          CALL scatter(tairmax_full, tairmax_slab)
1565          CALL scatter(tairmin_full, tairmin_slab)
1566       ENDIF
1567       CALL scatter(qair_full, qair_slab)
1568       CALL scatter(rainf_full, rainf_slab)
1569       CALL scatter(snowf_full, snowf_slab)
1570       CALL scatter(swdown_full, swdown_slab)
1571       CALL scatter(lwdown_full, lwdown_slab)
1572       CALL scatter(u_full, u_slab)
1573       CALL scatter(v_full, v_slab)
1574       CALL scatter(ps_full, ps_slab)
1575       CALL scatter(ztq_full, ztq_slab)
1576       CALL scatter(zuv_full, zuv_slab)
1577    ELSE
1578       tair_slab(:,:) = tair_full(:,:)
1579       IF ( forcing_tstep_ave >= one_day/3.0) THEN
1580          tairmax_slab(:,:) = tairmax_full(:,:)
1581          tairmin_slab(:,:) = tairmin_full(:,:)
1582       ENDIF
1583       qair_slab(:,:) = qair_full(:,:)
1584       rainf_slab(:,:) = rainf_full(:,:)
1585       snowf_slab(:,:) = snowf_full(:,:)
1586       swdown_slab(:,:) = swdown_full(:,:)
1587       lwdown_slab(:,:) = lwdown_full(:,:)
1588       u_slab(:,:) = u_full(:,:)
1589       v_slab(:,:) = v_full(:,:)
1590       ps_slab(:,:) = ps_full(:,:)
1591       ztq_slab(:,:) = ztq_full(:,:)
1592       zuv_slab(:,:) = zuv_full(:,:)
1593    ENDIF
1594    !
1595    !
1596    !
1597    DO is=1,nbpoint_proc
1598       preciptime_slab(is,:) = preciptime_tmp(:)
1599    ENDDO
1600    !
1601    ! Clean-up to free the memory on the root processor.
1602    !
1603    IF ( ALLOCATED(tair_full) ) DEALLOCATE(tair_full)
1604    IF ( ALLOCATED(tairmax_full) ) DEALLOCATE(tairmax_full)
1605    IF ( ALLOCATED(tairmin_full) ) DEALLOCATE(tairmin_full)
1606    IF ( ALLOCATED(qair_full) ) DEALLOCATE(qair_full)
1607    IF ( ALLOCATED(rainf_full) ) DEALLOCATE(rainf_full)
1608    IF ( ALLOCATED(snowf_full) ) DEALLOCATE(snowf_full)
1609    IF ( ALLOCATED(swdown_full) ) DEALLOCATE(swdown_full)
1610    IF ( ALLOCATED(lwdown_full) ) DEALLOCATE(lwdown_full)
1611    IF ( ALLOCATED(u_full) ) DEALLOCATE(u_full)
1612    IF ( ALLOCATED(v_full) ) DEALLOCATE(v_full)
1613    IF ( ALLOCATED(ps_full) ) DEALLOCATE(ps_full)
1614    IF ( ALLOCATED(ztq_full) ) DEALLOCATE(ztq_full)
1615    IF ( ALLOCATED(zuv_full) ) DEALLOCATE(zuv_full)
1616    !
1617  END SUBROUTINE forcing_readslab
1618!!
1619!!  =============================================================================================================================
1620!! SUBROUTINE: forcing_readslab_root
1621!!
1622!>\BRIEF Routine which reads a slab of data from the netCDF file and writes it onto the memory.
1623!!
1624!! DESCRIPTION: It is important to read the next slab of data while still keeping an overlap so that
1625!!              interpolation can continue.
1626!!              It also attributes a time axis to each variable.
1627!!
1628!! \n
1629!_ ==============================================================================================================================
1630
1631  SUBROUTINE forcing_readslab_root(time_int, &
1632            &                     tair, tairmax, tairmin, t_tair, tbnd_tair, &
1633            &                     qair, t_qair, tbnd_qair, &
1634            &                     rainf, snowf, t_prec, tbnd_prec, prectime, &
1635            &                     swdown, t_swdown, tbnd_swdown, &
1636            &                     lwdown, t_lwdown, tbnd_lwdown, &
1637            &                     u, t_u, tbnd_u, &
1638            &                     v, t_v, tbnd_v, &
1639            &                     ps, t_ps, tbnd_ps, &
1640            &                     ztq, zuv)
1641    !
1642    ! Arguments
1643    !
1644    REAL(r_std), INTENT(in)  :: time_int(2)                            !! The time interval over which the forcing is needed.
1645    !
1646    REAL(r_std), INTENT(out) :: tair(:,:), tairmax(:,:), tairmin(:,:)
1647    REAL(r_std), INTENT(out) :: t_tair(:)
1648    REAL(r_std), INTENT(out) :: tbnd_tair(:,:)
1649    !
1650    REAL(r_std), INTENT(out) :: qair(:,:)
1651    REAL(r_std), INTENT(out) :: t_qair(:)
1652    REAL(r_std), INTENT(out) :: tbnd_qair(:,:)
1653    !
1654    REAL(r_std), INTENT(out) :: rainf(:,:)
1655    REAL(r_std), INTENT(out) :: snowf(:,:)
1656    REAL(r_std), INTENT(out) :: t_prec(:)
1657    REAL(r_std), INTENT(out) :: tbnd_prec(:,:)
1658    REAL(r_std), INTENT(out) :: prectime(:)
1659    !
1660    REAL(r_std), INTENT(out) :: swdown(:,:)
1661    REAL(r_std), INTENT(out) :: t_swdown(:)
1662    REAL(r_std), INTENT(out) :: tbnd_swdown(:,:)
1663    !
1664    REAL(r_std), INTENT(out) :: lwdown(:,:)
1665    REAL(r_std), INTENT(out) :: t_lwdown(:)
1666    REAL(r_std), INTENT(out) :: tbnd_lwdown(:,:)
1667    !
1668    REAL(r_std), INTENT(out) :: u(:,:)
1669    REAL(r_std), INTENT(out) :: t_u(:)
1670    REAL(r_std), INTENT(out) :: tbnd_u(:,:)
1671    !
1672    REAL(r_std), INTENT(out) :: v(:,:)
1673    REAL(r_std), INTENT(out) :: t_v(:)
1674    REAL(r_std), INTENT(out) :: tbnd_v(:,:)
1675    !
1676    REAL(r_std), INTENT(out) :: ps(:,:)
1677    REAL(r_std), INTENT(out) :: t_ps(:)
1678    REAL(r_std), INTENT(out) :: tbnd_ps(:,:)
1679    !
1680    REAL(r_std), INTENT(out) :: ztq(:,:)
1681    REAL(r_std), INTENT(out) :: zuv(:,:)
1682    !
1683    ! Local
1684    !
1685    INTEGER(i_std) :: iret, varid
1686    INTEGER(i_std) :: if, it
1687    INTEGER(i_std) :: tstart(3), tcount(3)
1688    INTEGER(i_std) :: imin(1), imax(1), firstif(1)
1689    INTEGER(i_std) :: nctstart, nctcount, inslabpos
1690    INTEGER(i_std) :: start_globtime, end_globtime
1691    INTEGER(i_std) :: timeid_tair, timeid_qair, timeid_precip, timeid_swdown
1692    INTEGER(i_std) :: timeid_lwdown, timeid_u, timeid_v, timeid_ps, timeid_tmp
1693    REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: time_tmp
1694    REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: rau
1695    CHARACTER(LEN=80) :: cellmethod
1696    !
1697    LOGICAL, SAVE :: first_call_readslab=.TRUE.
1698    !
1699    ALLOCATE(time_tmp(slab_size,nbtax))
1700    ALLOCATE(rau(nbpoint_loc,slab_size))
1701    !
1702    !
1703    ! Catch any stupid utilisation of this routine.
1704    !
1705    IF ( .NOT. is_root_prc) THEN
1706       CALL ipslerr (3,'forcing_readslab_root',"Cannot run this routine o other procs than root.",&
1707            &        "All the information on the forcing files is only lated on the root processor."," ")
1708    ENDIF
1709    !
1710    !Set some variables to zero
1711    !
1712    IF ( first_call_readslab ) THEN
1713       !
1714       preciptime(:) = 0
1715       !
1716       ! If the first file is only there to provide the last time step of the previous year, we
1717       ! do not need to read all. We will start 2 forcing time steps before the start of the first
1718       ! time interval requested.
1719       !
1720       imin=MINLOC(ABS(time(:,1)-time_int(1)))
1721       current_offset = MAX(imin(1)-2,1)
1722       !
1723       first_call_readslab = .FALSE.
1724       write(numout, *) "first_call_readslab in forcing_readslab_root"
1725       !
1726    ELSE       
1727       !
1728       ! Put back the cummulated time of rainfall into the global array
1729       !
1730       preciptime(position_slab(1):position_slab(2)) = preciptime(position_slab(1):position_slab(2)) + &
1731            &    prectime(1:slab_size)
1732       !
1733       ! Compute new offset
1734       !
1735       current_offset = position_slab(2)-2
1736       write(numout, *) "first_call_readslab in forcing_readslab_root 22"
1737       !
1738    ENDIF
1739   
1740    !
1741    ! Check that the slab size is not too large
1742    !
1743    IF ( (current_offset-1)+slab_size > nb_forcing_steps) THEN
1744       slab_size = nb_forcing_steps - (current_offset-1)
1745    ENDIF
1746    !
1747    ! 1.1 Check that the slab we have to read exists
1748    !
1749    IF ( slab_size > 0 ) THEN
1750       !
1751       start_globtime = current_offset
1752       end_globtime = (current_offset-1)+slab_size
1753       inslabpos=1
1754       WRITE(*,*) ">> Reading from global position ", start_globtime, "up to ", end_globtime
1755       write(numout,*) time_sourcefile 
1756       !
1757       DO if=MINVAL(time_sourcefile(start_globtime:end_globtime)),MAXVAL(time_sourcefile(start_globtime:end_globtime))
1758          !
1759          ! Get position of the part of the global time axis we need to read from this file
1760          !
1761          firstif = MINLOC(ABS(time_sourcefile-if))
1762          ! start = distance of start_globtime or nothing + 1 to follow netCDF convention.
1763          nctstart =  MAX((start_globtime-firstif(1)), 0)+1
1764          ! count = end index - start index + 1
1765          nctcount = MIN((firstif(1)-1)+nbtime_perfile(if),end_globtime)-MAX(firstif(1),start_globtime)+1
1766          !
1767          !
1768          ! Read time over the indexes computed above in order to position the slab correctly
1769          !
1770          WRITE(*,*) ">> From file ", if," reading from position ", nctstart, "up to ", (nctstart-1)+nctcount
1771          !
1772          DO it =1,nbtax
1773             tstart(1) = nctstart
1774             tcount(1) = nctcount
1775             iret = NF90_GET_VAR(force_id(if), time_id(if,it), time_tmp(inslabpos:inslabpos+nctcount-1,it), tstart, tcount)
1776             IF (iret /= NF90_NOERR) THEN
1777                WRITE(*,*) TRIM(NF90_STRERROR(iret))
1778                WRITE(*,*) "Working on file ", IF," starting at ",tstart(1)," counting ",tcount(1)
1779                WRITE(*,*) "The data was to be written in to section ", inslabpos,":",inslabpos+nctcount-1," of time_tmp"
1780                CALL ipslerr (3,'forcing_readslab_root',"Could not read the time for the current interval."," "," ")
1781             ENDIF
1782             time_tmp(inslabpos:inslabpos+nctcount-1,it) = date0_file(if,it) + &
1783                  time_tmp(inslabpos:inslabpos+nctcount-1,it)*convtosec(if)/one_day
1784          ENDDO
1785          !
1786          ! 2.0 Find and read variables.
1787          !
1788          ! 2.1 Deal with air temperature and humidity as the fist and basic case
1789          !
1790          !
1791          !
1792          IF ( forcing_tstep_ave >= one_day/3.0) THEN
1793             CALL forcing_varforslab(if, "Tairmax", nctstart, nctcount, inslabpos, tairmax, cellmethod)
1794             CALL forcing_varforslab(if, "Tairmin", nctstart, nctcount, inslabpos, tairmin, cellmethod)
1795          ENDIF
1796         
1797         
1798          CALL forcing_varforslab(if, "Tair", nctstart, nctcount, inslabpos, tair, cellmethod)
1799          CALL forcing_attributetimeaxe(cellmethod, timeid_tair)
1800          !
1801          CALL forcing_varforslab(if, "Qair", nctstart, nctcount, inslabpos, qair, cellmethod)
1802          CALL forcing_attributetimeaxe(cellmethod, timeid_qair)
1803          !
1804          ! 2.2 Deal with rainfall and snowfall.
1805          !
1806          CALL forcing_varforslab(if, "Rainf", nctstart, nctcount, inslabpos, rainf, cellmethod)
1807          CALL forcing_attributetimeaxe(cellmethod, timeid_precip)
1808          !
1809          CALL forcing_varforslab(if, "Snowf", nctstart, nctcount, inslabpos, snowf, cellmethod)
1810          CALL forcing_attributetimeaxe(cellmethod, timeid_tmp)
1811          IF ( timeid_precip .NE. timeid_tmp) THEN
1812             CALL ipslerr(3, 'forcing_readslab_root','Rainf and Snwof have different time axes.', &
1813                  &         'Please check the forcing file to ensure both variable have the same cell_method.','')
1814          ENDIF
1815          !
1816          !
1817          ! 2.3 Deal with downward shorwave and longwave radiation
1818          !     The SW radiation can have 2 names SWdown_aerosol or SWdown. The first one is
1819          !     given priority
1820          !
1821          CALL forcing_varforslab(if, "SWdown", nctstart, nctcount, inslabpos, swdown, cellmethod)
1822          CALL forcing_attributetimeaxe(cellmethod, timeid_swdown)
1823          !
1824          CALL forcing_varforslab(if, "LWdown", nctstart, nctcount, inslabpos, lwdown, cellmethod)
1825          CALL forcing_attributetimeaxe(cellmethod, timeid_lwdown)
1826          !
1827          !
1828          ! 2.4 Deal with wind and pressure
1829          !
1830          CALL forcing_varforslab(if, "PSurf", nctstart, nctcount, inslabpos, ps, cellmethod)
1831          CALL forcing_attributetimeaxe(cellmethod, timeid_ps)
1832          !
1833          CALL forcing_varforslab(if, "Wind_E", nctstart, nctcount, inslabpos, u, cellmethod)
1834          CALL forcing_attributetimeaxe(cellmethod, timeid_u)
1835          !
1836          CALL forcing_varforslab(if, "Wind_N", nctstart, nctcount, inslabpos, v, cellmethod)
1837          CALL forcing_attributetimeaxe(cellmethod, timeid_v)
1838          !
1839          ! Verify on Tair that we have a credible field.
1840          !
1841          IF (   MINVAL(tair(:,inslabpos:inslabpos+nctcount-1)) < 100.0 .OR. &
1842               & MAXVAL(tair(:,inslabpos:inslabpos+nctcount-1)) > 500.0 ) THEN
1843             WRITE(*,*) "ERROR on range of Tair : ", MINVAL(tair(:,inslabpos:inslabpos+nctcount-1)), &
1844                  &                                  MAXVAL(tair(:,inslabpos:inslabpos+nctcount-1))
1845             CALL ipslerr(3, 'forcing_readslab_root','The air temperature is not in a credible range.', &
1846                  &         'Please verify your forcing file.','Are variables for all points to be simulated ?')
1847          ENDIF
1848          !
1849          ! Do the height of the variables T&Q and U&V
1850          !
1851          ! First the T&Q level
1852          !
1853          IF ( zheight ) THEN
1854             ztq(:,:) = zlev_fixed
1855          ELSE IF ( zsigma .OR. zhybrid ) THEN
1856             DO it=inslabpos,inslabpos+nctcount-1
1857                rau(:,it) = ps(:,it)/(cte_molr*tair(:,it))
1858                ztq(:,it) = (ps(:,it) - (zhybrid_a + zhybrid_b*ps(:,it)))/(rau(:,it) * cte_grav)
1859             ENDDO
1860          ELSE IF ( zlevels ) THEN
1861             CALL forcing_varforslab(IF, "Levels", nctstart, nctcount, inslabpos, ztq, cellmethod)
1862          ELSE
1863             CALL ipslerr(3, 'forcing_readslab_root','No case for the vertical levels was specified.', &
1864                  &         'We cannot determine the height for T and Q.','')
1865          ENDIF
1866          !
1867          ! Now the U&V level
1868          !
1869          IF ( zsamelev_uv ) THEN
1870             zuv(:,:) = ztq(:,:)
1871          ELSE
1872             IF ( zheight ) THEN
1873                zuv(:,:) = zlevuv_fixed
1874             ELSE IF ( zsigma .OR. zhybrid ) THEN
1875                DO it=inslabpos,inslabpos+nctcount-1
1876                   rau(:,it) = ps(:,it)/(cte_molr*tair(:,it))
1877                   zuv(:,it) = (ps(:,it) - (zhybriduv_a + zhybriduv_b*ps(:,it)))/(rau(:,it) * cte_grav)
1878                ENDDO
1879             ELSE IF ( zlevels ) THEN
1880                CALL forcing_varforslab(IF, "Levels_uv", nctstart, nctcount, inslabpos, zuv, cellmethod)
1881             ELSE
1882                CALL ipslerr(3, 'forcing_readslab_root','No case for the vertical levels was specified.', &
1883                     &         'We cannot determine the height for U and V.','stop readdim2')
1884             ENDIF
1885          ENDIF
1886         
1887          inslabpos = inslabpos+nctcount
1888         
1889       ENDDO
1890       !
1891       ! Use the read time of the slab to place it in the global variables. We consider
1892       ! that we can do that on the first axis.
1893       !
1894       imin = MINLOC(ABS(time(:,1)-time_tmp(1,1)))
1895       position_slab(1) = imin(1)
1896       imax = MINLOC(ABS(time(:,1)-time_tmp(slab_size,1)))
1897       position_slab(2) = imax(1)
1898       !
1899       !
1900       IF ( position_slab(2)-position_slab(1) .GT. slab_size ) THEN
1901          DO it =1,nbtax
1902             WRITE(*,*) "Checking time_tmp on idex : ", it
1903             WRITE(*,*) "Time_tmp start and end : ",time_tmp(1,it), time_tmp(slab_size,it)
1904             imin = MINLOC(ABS(time(:,1)-time_tmp(1,it)))
1905             imax = MINLOC(ABS(time(:,1)-time_tmp(slab_size,it)))
1906             WRITE(*,*) "Interval read : ", imax(1)-imin(1)+1
1907          ENDDO
1908          CALL ipslerr (3,'forcing_readslab_root',"The time slab read does not fit the number of variables read.",&
1909               &        "Could there be an error in the time axis ?"," ")
1910       ENDIF
1911       !
1912       ! Transfer the global time axis into the time variables approriate for each variable. This way
1913       ! the time axis for each variable will be centered on the interval of validity. This is an essential assumption
1914       ! the time interpolation to be done later.
1915       !
1916       WRITE(*,*) "We have found the following axes : ", time_axename(:)
1917       WRITE(*,*) "For Tair we are using time axis '",TRIM(time_axename(timeid_tair)),&
1918            &     "' with cell method ",TRIM(time_cellmethod(timeid_tair)),"."
1919       t_tair(1:slab_size) = time(position_slab(1):position_slab(2), timeid_tair)
1920       tbnd_tair(1:slab_size,:) = time_bounds(position_slab(1):position_slab(2),timeid_tair,:)
1921       !
1922       WRITE(*,*) "For Qair we are using time axis '",TRIM(time_axename(timeid_qair)),&
1923            &     "' with cell method ",TRIM(time_cellmethod(timeid_qair)),"."
1924       t_qair(1:slab_size) = time(position_slab(1):position_slab(2), timeid_qair)
1925       tbnd_qair(1:slab_size,:) = time_bounds(position_slab(1):position_slab(2),timeid_qair,:)
1926       !
1927       WRITE(*,*) "For Rainf and Snowf we are using time axis '",TRIM(time_axename(timeid_precip)),&
1928            &     "' with cell method ",TRIM(time_cellmethod(timeid_precip)),"."
1929       t_prec(1:slab_size) = time(position_slab(1):position_slab(2), timeid_precip)
1930       tbnd_prec(1:slab_size,:) = time_bounds(position_slab(1):position_slab(2),timeid_precip,:)
1931       prectime(1:slab_size) = preciptime(position_slab(1):position_slab(2))
1932       !
1933       WRITE(*,*) "For SWdown we are using time axis '",TRIM(time_axename(timeid_swdown)),&
1934            &     "' with cell method ",TRIM(time_cellmethod(timeid_swdown)),"."
1935       t_swdown(1:slab_size) = time(position_slab(1):position_slab(2), timeid_swdown)
1936       tbnd_swdown(1:slab_size,:) = time_bounds(position_slab(1):position_slab(2),timeid_swdown,:)
1937       !
1938       WRITE(*,*) "For LWdown we are using time axis '",TRIM(time_axename(timeid_lwdown)),&
1939            &     "' with cell method ",TRIM(time_cellmethod(timeid_lwdown)),"."
1940       t_lwdown(1:slab_size) = time(position_slab(1):position_slab(2), timeid_lwdown)
1941       tbnd_lwdown(1:slab_size,:) = time_bounds(position_slab(1):position_slab(2),timeid_lwdown,:)
1942       !
1943       WRITE(*,*) "For Wind_E we are using time axis '",TRIM(time_axename(timeid_u)),&
1944            &     "' with cell method ",TRIM(time_cellmethod(timeid_u)),"."
1945       t_u(1:slab_size) = time(position_slab(1):position_slab(2), timeid_u)
1946       tbnd_u(1:slab_size,:) = time_bounds(position_slab(1):position_slab(2),timeid_u,:)
1947       !
1948       WRITE(*,*) "For Wind_N we are using time axis '",TRIM(time_axename(timeid_v)),&
1949            &     "' with cell method ",TRIM(time_cellmethod(timeid_v)),"."
1950       t_v(1:slab_size) = time(position_slab(1):position_slab(2), timeid_v)
1951       tbnd_v(1:slab_size,:) = time_bounds(position_slab(1):position_slab(2),timeid_v,:)
1952       !
1953       WRITE(*,*) "For PSurf we are using time axis '",TRIM(time_axename(timeid_ps)),&
1954            &     "' with cell method ",TRIM(time_cellmethod(timeid_ps)),"."
1955       t_ps(1:slab_size) = time(position_slab(1):position_slab(2), timeid_ps)
1956       tbnd_ps(1:slab_size,:) = time_bounds(position_slab(1):position_slab(2),timeid_ps,:)
1957       !
1958    ELSE
1959       CALL ipslerr (2,'forcing_readslab_root',"We have reached the end of the slabs we can read.",&
1960            &          "The calling program needs to manage this situation"," ")
1961    ENDIF
1962    !
1963    ! Have we read to the end of the files ?
1964    !
1965    IF ( current_offset+slab_size >= nb_forcing_steps ) THEN
1966       end_of_file = .TRUE.
1967    ELSE
1968       end_of_file = .FALSE.
1969    ENDIF
1970    !
1971    IF ( ALLOCATED(rau) ) DEALLOCATE(rau)
1972    IF ( ALLOCATED(time_tmp) ) DEALLOCATE(time_tmp)
1973    !
1974  END SUBROUTINE forcing_readslab_root
1975
1976!!  =============================================================================================================================
1977!! SUBROUTINE: forcing_reindex3d
1978!!
1979!>\BRIEF     
1980!!
1981!! DESCRIPTION:   
1982!!
1983!! \n
1984!_ ==============================================================================================================================
1985  SUBROUTINE forcing_reindex3d(nbi, nbj, tin, slab_in, nbout, tout, slab_out, tstart, reindex)
1986    !
1987    ! ARGUMENTS
1988    !
1989    INTEGER(i_std), INTENT(in) :: nbi, nbj, tin, nbout, tout
1990    REAL(r_std), INTENT(in)    :: slab_in(nbi,nbj,tin)
1991    REAL(r_std), INTENT(out)   :: slab_out(nbout,tout)
1992    INTEGER(i_std), INTENT(in) :: tstart
1993    INTEGER(i_std), INTENT(in) :: reindex(nbout,2)
1994    !
1995    ! LOCAL
1996    !
1997    INTEGER(i_std) :: is, in
1998    !
1999    DO is=1,tin
2000       DO in=1,nbout
2001          slab_out(in,tstart+(is-1)) = slab_in(reindex(in,1),reindex(in,2),is)
2002       ENDDO
2003    ENDDO
2004    !
2005  END SUBROUTINE forcing_reindex3d
2006
2007!!  =============================================================================================================================
2008!! SUBROUTINE: forcing_reindex2d
2009!!
2010!>\BRIEF     
2011!!
2012!! DESCRIPTION:   
2013!!
2014!! \n
2015!_ ==============================================================================================================================
2016  SUBROUTINE forcing_reindex2d(nbi, nbj, slab_in, nbout, slab_out, reindex)
2017    !
2018    ! ARGUMENTS
2019    !
2020    INTEGER(i_std), INTENT(in) :: nbi, nbj, nbout
2021    REAL(r_std), INTENT(in)    :: slab_in(nbi,nbj)
2022    REAL(r_std), INTENT(out)   :: slab_out(nbout)
2023    INTEGER(i_std), INTENT(in) :: reindex(nbout,2)
2024    !
2025    ! LOCAL
2026    !
2027    INTEGER(i_std) :: in
2028    !
2029    DO in=1,nbout
2030       slab_out(in) = slab_in(reindex(in,1),reindex(in,2))
2031    ENDDO
2032    !
2033  END SUBROUTINE forcing_reindex2d
2034!!
2035!!  =============================================================================================================================
2036!! SUBROUTINE: forcing_reindex2dt
2037!!
2038!>\BRIEF     
2039!!
2040!! DESCRIPTION:   
2041!!
2042!! \n
2043!_ ==============================================================================================================================
2044
2045  SUBROUTINE forcing_reindex2dt(nbin, tin, slab_in, nbout, tout, slab_out, tstart, reindex)
2046    !
2047    ! ARGUMENTS
2048    !
2049    INTEGER(i_std), INTENT(in) :: nbin, tin, nbout, tout
2050    REAL(r_std), INTENT(in)    :: slab_in(nbin,tin)
2051    REAL(r_std), INTENT(out)   :: slab_out(nbout,tout)
2052    INTEGER(i_std), INTENT(in) :: tstart
2053    INTEGER(i_std), INTENT(in) :: reindex(nbout)
2054    !
2055    ! LOCAL
2056    !
2057    INTEGER(i_std) :: is, in
2058    !
2059    DO is=1,tin
2060       DO in=1,nbout
2061          slab_out(in,tstart+(is-1)) = slab_in(reindex(in),is)
2062       ENDDO
2063    ENDDO
2064    !
2065  END SUBROUTINE forcing_reindex2dt
2066
2067!!  =============================================================================================================================
2068!! SUBROUTINE: forcing_reindex1d
2069!!
2070!>\BRIEF     
2071!!
2072!! DESCRIPTION:   
2073!!
2074!! \n
2075!_ ==============================================================================================================================
2076
2077  SUBROUTINE forcing_reindex1d(nbin, slab_in, nbout, slab_out, reindex)
2078    !
2079    ! ARGUMENTS
2080    !
2081    INTEGER(i_std), INTENT(in) :: nbin, nbout
2082    REAL(r_std), INTENT(in)    :: slab_in(nbin)
2083    REAL(r_std), INTENT(out)   :: slab_out(nbout)
2084    INTEGER(i_std), INTENT(in) :: reindex(nbout)
2085    !
2086    ! LOCAL
2087    !
2088    INTEGER(i_std) :: is
2089    !
2090    DO is=1,nbout
2091       slab_out(is) = slab_in(reindex(is))
2092    ENDDO
2093    !
2094  END SUBROUTINE forcing_reindex1d
2095!!
2096!!  =============================================================================================================================
2097!! SUBROUTINE: forcing_reindex2to1
2098!!
2099!>\BRIEF     
2100!!
2101!! DESCRIPTION:   
2102!!
2103!! \n
2104!_ ==============================================================================================================================
2105
2106  SUBROUTINE forcing_reindex2to1(nbi, nbj, slab_in, nbout, slab_out, reindex)
2107    !
2108    ! ARGUMENTS
2109    !
2110    INTEGER(i_std), INTENT(in) :: nbi, nbj, nbout
2111    REAL(r_std), INTENT(in)    :: slab_in(nbi,nbj)
2112    REAL(r_std), INTENT(out)   :: slab_out(nbout)
2113    INTEGER(i_std), INTENT(in) :: reindex(nbout)
2114    !
2115    ! LOCAL
2116    !
2117    INTEGER(i_std) :: i, j, is
2118    !
2119    DO is=1,nbout
2120       j = INT((reindex(is)-1)/nbi)+1
2121       i = (reindex(is)-(j-1)*nbi)
2122       slab_out(is) = slab_in(i,j)
2123    ENDDO
2124    !
2125  END SUBROUTINE forcing_reindex2to1
2126
2127!!  =============================================================================================================================
2128!! SUBROUTINE: forcing_reindex1to2
2129!!
2130!>\BRIEF     
2131!!
2132!! DESCRIPTION:   
2133!!
2134!! \n
2135!_ ==============================================================================================================================
2136
2137  SUBROUTINE forcing_reindex1to2(nbin, slab_in, nbi, nbj, slab_out, reindex)
2138    !
2139    ! ARGUMENTS
2140    !
2141    INTEGER(i_std), INTENT(in) :: nbin, nbi, nbj
2142    REAL(r_std), INTENT(in)    :: slab_in(nbin)
2143    REAL(r_std), INTENT(out)   :: slab_out(nbi, nbj)
2144    INTEGER(i_std), INTENT(in) :: reindex(nbin)
2145    !
2146    ! LOCAL
2147    !
2148    INTEGER(i_std) :: i, j, is
2149    !
2150    DO is=1,nbin
2151       j = INT((reindex(is)-1)/nbi)+1
2152       i = (reindex(is)-(j-1)*nbi)
2153       slab_out(i,j) = slab_in(is)
2154    ENDDO
2155    !
2156  END SUBROUTINE forcing_reindex1to2
2157
2158!!  =============================================================================================================================
2159!! SUBROUTINE: forcing_close
2160!!
2161!>\BRIEF  Close all forcing files
2162!!
2163!! DESCRIPTION:   
2164!!
2165!! \n
2166!_ ==============================================================================================================================
2167  SUBROUTINE forcing_close()
2168
2169    INTEGER(i_std) :: ierr, if
2170
2171    DO if=1,nb_forcefile
2172       ierr = NF90_CLOSE(force_id(if))
2173    ENDDO
2174
2175  END SUBROUTINE forcing_close
2176
2177!!  =============================================================================================================================
2178!! SUBROUTINE: forcing_printdate
2179!!
2180!>\BRIEF    Print the date in a human readable format. 
2181!!
2182!! DESCRIPTION:   
2183!!
2184!! \n
2185!_ ==============================================================================================================================
2186
2187  SUBROUTINE forcing_printdate(julian_day, message, wunit)
2188    !
2189    REAL(r_std), INTENT(in) :: julian_day
2190    CHARACTER(len=*), INTENT(in) :: message
2191    INTEGER, INTENT(in), OPTIONAL :: wunit
2192    !
2193    !
2194    !
2195    INTEGER(i_std) :: year, month, day, hours, minutes, seci
2196    REAL(r_std) :: sec
2197    !
2198    CALL ju2ymds (julian_day, year, month, day, sec)
2199    hours = INT(sec/3600)
2200    sec = sec - 3600 * hours
2201    minutes = INT(sec / 60)
2202    sec = sec - 60 * minutes
2203    seci = INT(sec)
2204    !
2205    IF (PRESENT(wunit)) THEN
2206       WRITE(wunit,'(I4.4,"-",I2.2,"-",I2.2," ",I2.2,":",I2.2,":",I2.2," > ", A60)') &
2207            &            year, month, day, hours, minutes, seci, message
2208    ELSE
2209       WRITE(*,'(I4.4,"-",I2.2,"-",I2.2," ",I2.2,":",I2.2,":",I2.2," > ", A60)') &
2210            &            year, month, day, hours, minutes, seci, message
2211    ENDIF
2212    !
2213  END SUBROUTINE forcing_printdate
2214
2215!!  =============================================================================================================================
2216!! SUBROUTINE: forcing_printpoint_forgrid
2217!!
2218!>\BRIEF     Together with the date print some sample values. Useful for checking the forcing.
2219!!
2220!! DESCRIPTION:   
2221!!
2222!! \n
2223!_ ==============================================================================================================================
2224
2225  SUBROUTINE forcing_printpoint_forgrid(julian_day, lon_pt, lat_pt, var, message)
2226    !
2227    REAL(r_std), INTENT(in) :: julian_day
2228    REAL(r_std), INTENT(in) :: lon_pt, lat_pt
2229    REAL(r_std), INTENT(in) :: var(:)
2230    CHARACTER(len=*), INTENT(in) :: message
2231    !
2232    !
2233    !
2234    INTEGER(i_std) :: year, month, day, hours, minutes, seci
2235    REAL(r_std) :: sec
2236    INTEGER(i_std) :: lon_ind, lat_ind, ind
2237    INTEGER(i_std), DIMENSION(1) :: i,j,k
2238    !
2239    ! Check if there is anything to be done
2240    !
2241    IF ( MAX(lon_pt, lat_pt) > 360.0 ) THEN
2242       RETURN
2243    ENDIF
2244    !
2245    ! Convert time first
2246    !
2247    CALL ju2ymds (julian_day, year, month, day, sec)
2248    hours = INT(sec/3600)
2249    sec = sec - 3600 * hours
2250    minutes = INT(sec / 60)
2251    sec = sec - 60 * minutes
2252    seci = INT(sec)
2253    !
2254    ! Get the local to be analysed
2255    !
2256    i = MINLOC(ABS(lon_loc(:,1)-lon_pt))
2257    j = MINLOC(ABS(lat_loc(1,:)-lat_pt))
2258    ind = (j(1)-1)*iim_loc+i(1)
2259    k = MINLOC(ABS(lindex_loc(:)-ind))
2260    !
2261    WRITE(*,"(I2.2,':',I2.2,':',I2.2,' Loc : ', F5.1,',', F5.1,'(i=',I6,') Value = ',F12.4,A40)") &
2262         & hours, minutes, seci, lon_loc(i(1),1), lat_loc(1,j(1)), k(1), var(k(1)), message
2263   
2264  END SUBROUTINE forcing_printpoint_forgrid
2265  !!  =============================================================================================================================
2266!! SUBROUTINE: forcing_printpoint_forgrid2d
2267!!
2268!>\BRIEF     Together with the date print some sample values. Useful for checking the forcing.
2269!!
2270!! DESCRIPTION:   
2271!!
2272!! \n
2273!_ ==============================================================================================================================
2274
2275  SUBROUTINE forcing_printpoint_forgrid2d(julian_day, lon_pt, lat_pt, var, message)
2276    !
2277    REAL(r_std), INTENT(in) :: julian_day
2278    REAL(r_std), INTENT(in) :: lon_pt, lat_pt
2279    REAL(r_std), INTENT(in) :: var(:,:)
2280    CHARACTER(len=*), INTENT(in) :: message
2281    !
2282    !
2283    !
2284    INTEGER(i_std) :: year, month, day, hours, minutes, seci
2285    REAL(r_std) :: sec
2286    INTEGER(i_std) :: lon_ind, lat_ind
2287    INTEGER(i_std), DIMENSION(1) :: i,j
2288    !
2289    ! Check if there is anything to be done
2290    !
2291    IF ( MAX(lon_pt, lat_pt) > 360.0 ) THEN
2292       RETURN
2293    ENDIF
2294    !
2295    ! Convert time first
2296    !
2297    CALL ju2ymds (julian_day, year, month, day, sec)
2298    hours = INT(sec/3600)
2299    sec = sec - 3600 * hours
2300    minutes = INT(sec / 60)
2301    sec = sec - 60 * minutes
2302    seci = INT(sec)
2303    !
2304    ! Get the local to be analysed
2305    !
2306    i = MINLOC(ABS(lon_loc(:,1)-lon_pt))
2307    j = MINLOC(ABS(lat_loc(1,:)-lat_pt))
2308    !
2309    WRITE(*,"(I2.2,':',I2.2,':',I2.2,' Loc : ', F5.1,',', F5.1,'(i=',I6,') Value = ',F12.4,A40)") &
2310         & hours, minutes, seci, lon_loc(i(1),1), lat_loc(1,j(1)), i(1), j(1), var(i(1),j(1)), message
2311   
2312  END SUBROUTINE forcing_printpoint_forgrid2d
2313
2314!!  =============================================================================================================================
2315!! SUBROUTINE: forcing_printpoint_gen
2316!!
2317!>\BRIEF       Together with the date print some sample values. Useful for checking the forcing.
2318!!
2319!! DESCRIPTION:   
2320!!
2321!! \n
2322!_ ==============================================================================================================================
2323
2324  SUBROUTINE forcing_printpoint_gen(julian_day, lon_pt, lat_pt, nbind, lalo_in, var, message, ktest)
2325    !
2326    REAL(r_std), INTENT(in) :: julian_day
2327    REAL(r_std), INTENT(in) :: lon_pt, lat_pt
2328    INTEGER(i_std), INTENT(in) :: nbind
2329    REAL(r_std), INTENT(in) :: lalo_in(:,:)
2330    REAL(r_std), INTENT(in) :: var(:)
2331    CHARACTER(len=*), INTENT(in) :: message
2332    INTEGER(i_std), OPTIONAL, INTENT(out) :: ktest
2333    !
2334    !
2335    !
2336    INTEGER(i_std) :: year, month, day, hours, minutes, seci
2337    REAL(r_std) :: sec, mindist
2338    REAL(r_std), ALLOCATABLE, DIMENSION(:) :: dist, refdist
2339    INTEGER(i_std) :: lon_ind, lat_ind, ind
2340    INTEGER(i_std) :: i, imin(1)
2341    REAL(r_std), PARAMETER :: mincos  = 0.0001
2342    REAL(r_std), PARAMETER :: pi = 3.141592653589793238
2343    REAL(r_std), PARAMETER :: R_Earth = 6378000.
2344    !
2345    ! Check if there is anything to be done
2346    !
2347    IF ( MAX(lon_pt, lat_pt) > 360.0 ) THEN
2348       IF ( PRESENT(ktest) ) ktest = -1
2349       RETURN
2350    ENDIF
2351    !
2352    ! Allocate memory
2353    !
2354    ALLOCATE(dist(nbind))
2355    ALLOCATE(refdist(nbind))
2356    !
2357    ! Convert time first
2358    !
2359    CALL ju2ymds (julian_day, year, month, day, sec)
2360    hours = INT(sec/3600)
2361    sec = sec - 3600 * hours
2362    minutes = INT(sec / 60)
2363    sec = sec - 60 * minutes
2364    seci = INT(sec)
2365    !
2366    ! Get the location to be analysed
2367    !
2368    DO i=1,nbind
2369       dist(i) = acos( sin(lat_pt*pi/180)*sin(lalo_in(i,1)*pi/180) + &
2370            &    cos(lat_pt*pi/180)*cos(lalo_in(i,1)*pi/180)*&
2371            &    cos((lalo_in(i,2)-lon_pt)*pi/180) ) * R_Earth
2372    ENDDO
2373    !
2374    ! Look for the next grid point closest to the one with the smalest distance.
2375    !
2376    imin = MINLOC(dist)
2377    DO i=1,nbind
2378       refdist(i) = acos( sin(lalo_in(imin(1),1)*pi/180)*sin(lalo_in(i,1)*pi/180) + &
2379            &       cos(lalo_in(imin(1),1)*pi/180)*cos(lalo_in(i,1)*pi/180) * &
2380            &       cos((lalo_in(i,2)-lalo_in(imin(1),2))*pi/180) ) * R_Earth
2381    ENDDO
2382    refdist(imin(1)) =  MAXVAL(refdist)
2383    mindist = MINVAL(refdist)
2384    !
2385    ! Are we closer than the closest points ?
2386    !
2387    IF ( PRESENT(ktest) ) ktest = -1
2388    IF ( dist(imin(1)) <= mindist ) THEN
2389       !
2390       WRITE(*,"(I2.2,':',I2.2,':',I2.2,' Loc : ', F6.1,',', F6.1,'(i=',I6,') Value = ',F12.4,A38)") &
2391            & hours, minutes, seci, lalo_in(imin(1),2), lalo_in(imin(1),1), imin(1), var(imin(1)), message
2392       !
2393       IF ( PRESENT(ktest) ) ktest = imin(1)
2394    ENDIF
2395    !
2396  END SUBROUTINE forcing_printpoint_gen
2397
2398!!  =============================================================================================================================
2399!! SUBROUTINE: forcing_getglogrid
2400!!
2401!>\BRIEF       Routine to read the full grid of the forcing file.
2402!!
2403!! DESCRIPTION: The data is stored in the saved variables of the module and serve all other routines.     
2404!!
2405!! \n
2406!_ ==============================================================================================================================
2407
2408  SUBROUTINE forcing_getglogrid (nbfiles, filename, iim_tmp, jjm_tmp, nbpoint_tmp, closefile, landonly_arg)
2409    !
2410    ! This routine reads the global grid information from the forcing file and generates all the
2411    ! information needed for this grid.
2412    !
2413    ! ARGUMENTS
2414    !
2415    INTEGER(i_std), INTENT(in)    :: nbfiles
2416    CHARACTER(LEN=*), INTENT(in)  :: filename(:)
2417    INTEGER(i_std), INTENT(out)   :: iim_tmp, jjm_tmp, nbpoint_tmp
2418    LOGICAL, INTENT(in)           :: closefile
2419    LOGICAL, OPTIONAL, INTENT(in) :: landonly_arg
2420    !
2421    ! LOCAL
2422    !
2423    INTEGER(i_std) :: iret, iv, if, lll
2424    CHARACTER(LEN=20) :: dimname, varname
2425    CHARACTER(LEN=60) :: lon_units, lat_units, units
2426    INTEGER(i_std), ALLOCATABLE, DIMENSION(:) :: dimids, londim_ids, latdim_ids
2427    INTEGER(i_std) :: lon_id, lat_id, land_id, lon_nbdims, lat_nbdims, land_nbdims
2428    INTEGER(i_std) :: lonvar_id, latvar_id, landvar_id
2429    LOGICAL :: dump_mask
2430    INTEGER(i_std) :: ik, i, j, iff, ndimsvar
2431    ! Read a test variabe
2432    CHARACTER(len=8) :: testvarname='Tair'
2433    INTEGER(i_std)   :: testvar_id, contfrac_id
2434    REAL(r_std) :: testvar_missing, contfrac_missing
2435    REAL(r_std), ALLOCATABLE, DIMENSION(:)   :: testvar
2436    REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: testvar2d, contfrac2d
2437    !
2438    ! 0.0 Check variables are allocated
2439    !
2440    IF ( .NOT. ALLOCATED(force_id)) ALLOCATE(force_id(nbfiles))
2441    IF ( .NOT. ALLOCATED(id_unlim)) ALLOCATE(id_unlim(nbfiles))
2442    IF ( .NOT. ALLOCATED(nb_atts)) ALLOCATE(nb_atts(nbfiles))
2443    IF ( .NOT. ALLOCATED(ndims)) ALLOCATE(ndims(nbfiles))
2444    IF ( .NOT. ALLOCATED(nvars)) ALLOCATE( nvars(nbfiles))
2445    !
2446    ! 0.1 Are we one the root proc ?
2447    !
2448    IF ( .NOT. is_root_prc ) THEN
2449       CALL ipslerr (3,'forcing_getglogrid'," This routine can only be called on the root processor.", " ", " ")
2450    ENDIF
2451    !
2452    ! The default behaviour is to provide only land points to the calling program.
2453    ! But for forcing ocean model there is also the option to pass on land and ocean values.
2454    ! When the grid is initialized landonly_tmp=.FALSE. has to be set to obtian this behaviour.
2455    !
2456    IF ( PRESENT(landonly_arg) ) THEN
2457       landonly=landonly_arg
2458    ELSE
2459       landonly=.TRUE.
2460    ENDIF
2461    !
2462    ! 1.0 Open the netCDF file and get basic dimensions
2463    !
2464    DO iff=1,nbfiles
2465       !
2466       iret = NF90_OPEN(filename(iff), NF90_NOWRITE, force_id(iff))
2467       IF (iret /= NF90_NOERR) THEN
2468          CALL ipslerr (3,'forcing_getglogrid',"Error opening the forcing file :", filename(iff), " ")
2469       ENDIF
2470       !
2471       iret = NF90_INQUIRE (force_id(iff), nDimensions=ndims(iff), nVariables=nvars(iff), &
2472            nAttributes=nb_atts(iff), unlimitedDimId=id_unlim(iff))
2473       !
2474       !
2475       ! 2.0 Read the dimensions found in the forcing file. Only deal with the spatial dimension as
2476       !     time is an unlimited dimension.
2477       !
2478       DO iv=1,ndims(iff)
2479          !
2480          iret = NF90_INQUIRE_DIMENSION (force_id(iff), iv, name=dimname, len=lll)
2481          IF (iret /= NF90_NOERR) THEN
2482             CALL ipslerr (3,'forcing_getglogrid',"Could not get size of dimensions in file : ", filename(iff), " ")
2483          ENDIF
2484          !
2485          SELECT CASE(lowercase(dimname))
2486             !
2487          CASE("west_east")
2488             CALL forcing_checkdim(iff, filename, iim_glo, lon_id, lll, iv)
2489          CASE("south_north")
2490             CALL forcing_checkdim(iff, filename, jjm_glo, lat_id, lll, iv)
2491          CASE("longitude")
2492             CALL forcing_checkdim(iff, filename, iim_glo, lon_id, lll, iv)
2493          CASE("latitude")
2494             CALL forcing_checkdim(iff, filename, jjm_glo, lat_id, lll, iv)
2495          CASE("lon")
2496             CALL forcing_checkdim(iff, filename, iim_glo, lon_id, lll, iv)
2497          CASE("lat")
2498             CALL forcing_checkdim(iff, filename, jjm_glo, lat_id, lll, iv)
2499          CASE("nav_lon")
2500             CALL forcing_checkdim(iff, filename, iim_glo, lon_id, lll, iv)
2501          CASE("nav_lat")
2502             CALL forcing_checkdim(iff, filename, jjm_glo, lat_id, lll, iv)
2503          CASE("x")
2504             CALL forcing_checkdim(iff, filename, iim_glo, lon_id, lll, iv)
2505          CASE("y")
2506             CALL forcing_checkdim(iff, filename, jjm_glo, lat_id, lll, iv)
2507          CASE("land")
2508             CALL forcing_checkdim(iff, filename, nbland_glo, land_id, lll, iv)
2509          END SELECT
2510          !
2511       ENDDO
2512       IF ( iim_glo == 0 .AND. jjm_glo == 0 ) THEN
2513          CALL ipslerr (3,'forcing_getglogrid',"Did not recognize any dimensions in : ", filename(iff), " ")
2514       ENDIF
2515    ENDDO
2516   
2517    !
2518    ! 3.0 Read the spatial coordinate variables found in the first file.
2519    !
2520    ALLOCATE(dimids(NF90_MAX_VAR_DIMS), londim_ids(NF90_MAX_VAR_DIMS), latdim_ids(NF90_MAX_VAR_DIMS))
2521    lonvar_id = -1
2522    latvar_id = -1
2523    landvar_id = -1
2524    testvar_id = -1
2525    contfrac_id = -1
2526    ! Count the number of time axis we have
2527    nbtax = 0
2528    !
2529    DO iv=1,nvars(1)
2530       !
2531       iret = NF90_INQUIRE_VARIABLE(force_id(1), iv, name=varname, ndims=ndimsvar, dimids=dimids)
2532       iret = NF90_GET_ATT(force_id(1), iv, 'units', units)
2533       !
2534       ! Check that we have the longitude
2535       !
2536       IF ( INDEX(lowercase(varname), 'lon') > 0 .AND. INDEX(lowercase(units), 'east') > 0 ) THEN
2537          lonvar_id = iv
2538          lon_units=units
2539          lon_nbdims = ndimsvar
2540          londim_ids = dimids
2541       ENDIF
2542       !
2543       ! Check that we have the latitude
2544       !
2545       IF ( INDEX(lowercase(varname), 'lat') > 0 .AND. INDEX(lowercase(units), 'north') > 0) THEN
2546          latvar_id = iv
2547          lat_units=units
2548          lat_nbdims = ndimsvar
2549          latdim_ids = dimids
2550       ENDIF
2551       !
2552       ! Check that we have the land re-indexing table
2553       !
2554       IF ( INDEX(lowercase(varname), 'land') > 0 ) THEN
2555          landvar_id = iv
2556          land_nbdims = ndimsvar
2557          latdim_ids = dimids
2558       ENDIF
2559       !
2560       ! Check if we have the contfrac variable
2561       !
2562       IF ( INDEX(lowercase(varname), 'contfrac') > 0 ) THEN
2563          contfrac_id = iv
2564          iret = NF90_GET_ATT(force_id(1), iv, "missing_value", contfrac_missing)
2565          IF (iret /= NF90_NOERR) THEN
2566             ! No missing_value found, try to read _FillValue instead
2567             iret = NF90_GET_ATT(force_id(1), iv, "_FillValue", contfrac_missing)
2568             IF (iret /= NF90_NOERR) THEN
2569                WRITE(*,*) TRIM(nf90_strerror(iret))
2570                WRITE(*,*) " >> No _FillValue or missing_value found for contfrac"
2571                contfrac_missing=0.0
2572             END IF
2573          ENDIF
2574       ENDIF
2575       !
2576       ! Find the test variable
2577       !
2578       IF ( INDEX(lowercase(varname), TRIM(lowercase(testvarname))) > 0 .AND. &
2579            & LEN_TRIM(varname) == LEN_TRIM(testvarname)) THEN
2580          testvar_id = iv
2581          iret = NF90_GET_ATT(force_id(1), iv, "missing_value", testvar_missing)
2582          IF (iret /= NF90_NOERR) THEN
2583             ! No missing_value found, try to read _FillValue instead
2584             iret = NF90_GET_ATT(force_id(1), iv, "_FillValue", testvar_missing)
2585             IF (iret /= NF90_NOERR) THEN
2586                WRITE(*,*) TRIM(nf90_strerror(iret))
2587                WRITE(*,*) " >> No _FillValue or missing_value found for variable=",varname
2588                testvar_missing=-1
2589             END IF
2590          ENDIF
2591       ENDIF
2592       !
2593       ! If we come across time get the calendar and archive it.
2594       !
2595       IF ( INDEX(lowercase(units),'seconds since') > 0 .OR. &
2596            &  INDEX(lowercase(units),'minutes since') > 0 .OR. &
2597            &  INDEX(lowercase(units),'hours since') > 0 .OR. &
2598          &  INDEX(lowercase(units),'days since') > 0) THEN 
2599          !
2600          ! Get calendar used for the time axis
2601          !
2602          iret = NF90_GET_ATT(force_id(1), iv, "calendar", calendar)
2603          IF (iret == NF90_NOERR) THEN
2604             WRITE(*,*) ">> Setting the calendar to ",calendar
2605          ELSE
2606             WRITE(*,*) ">> Keeping proleptic Gregorian calendar" 
2607             calendar="proleptic_gregorian"
2608          ENDIF
2609          !
2610          nbtax = nbtax + 1
2611          !
2612       ENDIF
2613    ENDDO
2614    !
2615    ! 4.0 Verification that we have found both coordinate variables and the land point index
2616    !
2617    IF ( ( lonvar_id < 0 ) .AND. ( INDEX(lowercase(lon_units), 'east') <= 0 ) ) THEN
2618       CALL ipslerr (3,'forcing_getglogrid',"Have not found a valid longitude. Units = ", lon_units, " ")
2619    ENDIF
2620    IF ( ( latvar_id < 0 ) .AND. ( INDEX(lowercase(lat_units), 'north') <= 0 ) ) THEN
2621       CALL ipslerr (3,'forcing_getglogrid',"Have not found a valid latitude. Units = : ", lat_units, " ")
2622    ENDIF
2623    IF ( landvar_id < 0 ) THEN
2624       CALL ipslerr (1,'forcing_getglogrid',"No reindexing table was found. ", &
2625            &           "This forcing file is not compressed by gathering.", " ")
2626    ENDIF
2627    !
2628    ! 5.0 Allocate the space for the global variables and read them.
2629    !
2630    IF ( .NOT. ALLOCATED(lon_glo)) ALLOCATE(lon_glo(iim_glo, jjm_glo))
2631    IF ( .NOT. ALLOCATED(lat_glo)) ALLOCATE(lat_glo(iim_glo, jjm_glo))
2632    !
2633    IF ( lon_nbdims == 2 .AND. lat_nbdims == 2 ) THEN
2634       iret = NF90_GET_VAR(force_id(1), lonvar_id, lon_glo)
2635       iret = NF90_GET_VAR(force_id(1), latvar_id, lat_glo)
2636    ELSE IF ( lon_nbdims == 1 .AND. lat_nbdims == 1 ) THEN
2637       DO iv=1,jjm_glo
2638          iret = NF90_GET_VAR(force_id(1), lonvar_id, lon_glo(:,iv))
2639       ENDDO
2640       DO iv=1,iim_glo
2641          iret = NF90_GET_VAR(force_id(1), latvar_id, lat_glo(iv,:))
2642       ENDDO
2643    ELSE
2644       WRITE(*,*) "Found dimensions for lon and lat :", lon_nbdims, lat_nbdims
2645       CALL ipslerr (3,'forcing_getglogrid',"Longitude and Latitude variables do not have the right dimensions.", " ", " ")
2646    ENDIF
2647    !
2648    ! If we have a indexing variable then the data is compressed by gathering, else we have to construct it.
2649    !
2650    compressed = .FALSE.
2651    IF ( landvar_id > 0 ) THEN
2652       IF ( .NOT. ALLOCATED(lindex_glo)) ALLOCATE(lindex_glo(nbland_glo))
2653       iret = NF90_GET_VAR(force_id(1), landvar_id, lindex_glo)
2654       compressed = .TRUE.
2655    ENDIF
2656    !
2657    IF ( .NOT. ALLOCATED(mask_glo)) ALLOCATE(mask_glo(iim_glo, jjm_glo)) 
2658    !
2659    ! Get the land/sea mask and contfrac based on a test variable, if contfract is not available. Else
2660    ! we use the contfrac variable.
2661    !
2662    IF ( compressed ) THEN
2663       IF ( landonly ) THEN
2664          IF ( .NOT. ALLOCATED(contfrac_glo)) ALLOCATE(contfrac_glo(nbland_glo))
2665          IF ( .NOT. ALLOCATED(testvar)) ALLOCATE(testvar(nbland_glo))
2666          CALL forcing_contfrac1d(force_id(1), testvar_id, contfrac_id, testvar)
2667          nbpoint_glo = nbland_glo
2668       ELSE
2669          WRITE(*,*) "forcing_tools cannot provide data over ocean points as the"
2670          WRITE(*,*) "data in the file is compressed by gathering land points."
2671          WRITE(*,*) "Fatal error"
2672          CALL ipslerr (3,'forcing_getglogrid',"forcing_tools cannot provide data over ocean points as the", &
2673               &                               "data in the file is compressed by gathering land points.", " ")
2674       ENDIF
2675    ELSE
2676       IF ( .NOT. ALLOCATED(testvar2d)) ALLOCATE(testvar2d(iim_glo, jjm_glo))
2677       IF ( .NOT. ALLOCATED(contfrac2d)) ALLOCATE(contfrac2d(iim_glo, jjm_glo))
2678       CALL forcing_contfrac2d(force_id(1), testvar_id, contfrac_id, testvar2d, contfrac2d, &
2679            & testvar_missing, contfrac_missing, nbland_glo)
2680       
2681
2682       !
2683       ! We have found a variable on which we can count the number of land points. We can build
2684       ! the indexing table and gather the information needed.
2685       !
2686       IF ( landonly ) THEN
2687          nbpoint_glo = nbland_glo
2688         
2689   
2690          IF ( .NOT. ALLOCATED(lindex_glo)) ALLOCATE(lindex_glo(nbpoint_glo))
2691          IF ( .NOT. ALLOCATED(contfrac_glo)) ALLOCATE(contfrac_glo(nbpoint_glo))
2692          IF ( .NOT. ALLOCATED(testvar)) ALLOCATE(testvar(nbpoint_glo))
2693         
2694 
2695          IF ( contfrac_id > 0 ) THEN
2696             CALL forcing_buildindex(contfrac2d, contfrac_missing, landonly, lindex_glo, contfrac_glo)
2697             CALL forcing_reindex(iim_glo, jjm_glo, testvar2d, nbland_glo, testvar, lindex_glo)
2698             
2699
2700          ELSE
2701             CALL forcing_buildindex(testvar2d, testvar_missing, landonly, lindex_glo, testvar)
2702             contfrac_glo(:) = 1.0
2703             
2704          ENDIF
2705       ELSE
2706          nbpoint_glo = iim_glo*jjm_glo
2707          IF ( .NOT. ALLOCATED(lindex_glo)) ALLOCATE(lindex_glo(nbpoint_glo))
2708          IF ( .NOT. ALLOCATED(contfrac_glo)) ALLOCATE(contfrac_glo(nbpoint_glo))
2709          IF ( .NOT. ALLOCATED(testvar)) ALLOCATE(testvar(nbpoint_glo))
2710          IF ( contfrac_id > 0 ) THEN
2711             CALL forcing_buildindex(contfrac2d, contfrac_missing, landonly, lindex_glo, contfrac_glo)
2712             CALL forcing_reindex(iim_glo, jjm_glo, testvar2d, nbland_glo, testvar, lindex_glo)
2713          ELSE
2714             CALL forcing_buildindex(testvar2d, testvar_missing, landonly, lindex_glo, testvar)
2715             contfrac_glo(:) = 1.0
2716          ENDIF
2717       ENDIF
2718       !
2719    ENDIF
2720    !
2721    !
2722    ! We assume that if the forcing file is closed at the end of this subroutine this is a test
2723    ! or initialisation of the grids. So we will dump the mask in a netCDF file for the user to
2724    ! check.
2725    !
2726    dump_mask = closefile 
2727    CALL forcing_checkindex(dump_mask, testvarname, testvar)
2728   
2729    !
2730    !
2731    ! 8.0 Close file if needed
2732    !
2733    IF ( closefile ) THEN
2734       CALL forcing_close()
2735    ENDIF
2736    !
2737    ! Prepare variables to be returnned to calling subroutine.
2738    !
2739    iim_tmp = iim_glo
2740    jjm_tmp = jjm_glo
2741    nbpoint_tmp = nbpoint_glo
2742   
2743    !
2744    ! Clean up !
2745    !
2746    IF ( ALLOCATED(dimids) ) DEALLOCATE(dimids)
2747    IF ( ALLOCATED(londim_ids) ) DEALLOCATE(londim_ids)
2748    IF ( ALLOCATED(latdim_ids) ) DEALLOCATE(latdim_ids)
2749    IF ( ALLOCATED(testvar) ) DEALLOCATE(testvar)
2750    IF ( ALLOCATED(testvar2d) ) DEALLOCATE(testvar2d)
2751    IF ( ALLOCATED(contfrac2d) ) DEALLOCATE(contfrac2d)
2752    !
2753  END SUBROUTINE forcing_getglogrid
2754
2755!!  =============================================================================================================================
2756!! SUBROUTINE: forcing_buildindex
2757!!
2758!>\BRIEF     
2759!!
2760!! DESCRIPTION: When the forcing file does not contain compressed variables we need
2761!!              to build the land index variable from the mask defined by missing variables in
2762!!              a test variable. 
2763!!
2764!! \n
2765!_ ==============================================================================================================================
2766
2767  SUBROUTINE forcing_buildindex(var2d, var_missing, landonly, lindex, var)
2768    !
2769    ! When the forcing file does not contain compressed variables we need
2770    ! to build the land index variable from the mask defined by missing variables in
2771    ! a test variable.
2772    !
2773    ! Arguments
2774    !
2775    REAL(r_std), INTENT(in) :: var2d(:,:)
2776    REAL(r_std), INTENT(in) :: var_missing
2777    LOGICAL, INTENT(in)     :: landonly
2778    INTEGER(i_std), INTENT(out) :: lindex(:)
2779    REAL(r_std), INTENT(out) :: var(:)
2780    !
2781    ! Local
2782    !
2783    INTEGER(i_std) :: i,j,k
2784
2785    IF ( MAXVAL(var2d) >= var_missing ) THEN
2786       ! Case when we have missing values in the vard2d
2787       k=0
2788       DO i=1,iim_glo
2789          DO j=1,jjm_glo
2790             IF ( landonly ) THEN
2791                IF ( var2d(i,j) /= var_missing .AND. var2d(i,j) > 0.0 ) THEN
2792                   k = k + 1
2793                   lindex(k) = (j-1)*iim_glo+i 
2794                   var(k) = var2d(i,j)
2795                ENDIF
2796             ELSE
2797                ! When we take all point, no test is performed.
2798                k = k + 1
2799                lindex(k) = (j-1)*iim_glo+i 
2800                var(k) = var2d(i,j)
2801             ENDIF
2802          ENDDO
2803       ENDDO
2804    ELSE
2805       ! We suppose that this is land fraction variable
2806       k=0
2807       DO i=1,iim_glo
2808          DO j=1,jjm_glo
2809             IF ( landonly ) THEN
2810                IF ( var2d(i,j) > 0.0 ) THEN
2811                   k = k + 1
2812                   lindex(k) = (j-1)*iim_glo+i 
2813                   var(k) = var2d(i,j)
2814                ENDIF
2815             ELSE
2816                ! When we take all point, no test is performed.
2817                k = k + 1
2818                lindex(k) = (j-1)*iim_glo+i 
2819                var(k) = var2d(i,j)
2820             ENDIF
2821          ENDDO
2822       ENDDO
2823    ENDIF
2824    !
2825    !
2826  END SUBROUTINE forcing_buildindex
2827
2828!!  =============================================================================================================================
2829!! SUBROUTINE: forcing_contfrac1d
2830!!
2831!>\BRIEF     
2832!!
2833!! DESCRIPTION:  This routine build the land/mask if needed and gets the contfrac variable from forcing file.
2834!!               Here we treat the case where the variables are compressed by gathering. Thus only
2835!!               land points are available in the file.
2836!!
2837!! \n
2838!_ ==============================================================================================================================
2839
2840  SUBROUTINE forcing_contfrac1d(ifile, testvar_id, contfrac_id, testvar)
2841    !
2842    ! This routine build the land/mask if needed and gets the contfrac variable from forcing file.
2843    ! Here we treat the case where the variables are compressed by gathering. Thus only
2844    ! land points are available in the file.
2845    !
2846    ! ARGUMENTS
2847    !
2848    INTEGER(i_std), INTENT(in)               :: ifile
2849    INTEGER(i_std), INTENT(in)               :: testvar_id, contfrac_id
2850    REAL(r_std), DIMENSION(:), INTENT(inout) :: testvar 
2851    !
2852    ! LOCAL
2853    !
2854    INTEGER(i_std)                           :: it, iret
2855    INTEGER(i_std), DIMENSION(3)             :: start, count
2856    REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: contfrac2d
2857    !
2858    ! First determine the contfrac variable
2859    !
2860   
2861    IF ( contfrac_id > 0 ) THEN
2862       iret = NF90_INQUIRE_VARIABLE(ifile, contfrac_id, ndims=it)
2863       IF ( it == 1 ) THEN
2864          start = (/1,1,0/)
2865          count = (/nbpoint_glo,1,0/)
2866          iret = NF90_GET_VAR(ifile, contfrac_id, contfrac_glo, start, count)
2867          IF (iret /= NF90_NOERR) THEN
2868             WRITE(*,*) TRIM(nf90_strerror(iret))
2869             CALL ipslerr (3,'forcing_contfrac1d',"Error reading contfrac ", " ", " ")
2870          ENDIF
2871       ELSE IF ( it == 2 ) THEN
2872          ALLOCATE(contfrac2d(iim_glo,jjm_glo))
2873          start = (/1,1,0/)
2874          count = (/iim_glo,jjm_glo,0/)
2875          iret = NF90_GET_VAR(ifile, contfrac_id, contfrac2d)
2876          IF (iret /= NF90_NOERR) THEN
2877             WRITE(*,*) TRIM(nf90_strerror(iret))
2878             CALL ipslerr (3,'forcing_contfrac1d',"Error reading contfrac ", " ", " ")
2879          ENDIF
2880          CALL forcing_reindex(iim_glo, jjm_glo, contfrac2d, nbpoint_glo, contfrac_glo, lindex_glo)
2881          DEALLOCATE(contfrac2d)
2882       ELSE
2883          CALL ipslerr (3,'forcing_contfrac1d',"Contfrac has a dimension larger than 2. ", &
2884               "We do not know how to handle this.", " ")
2885       ENDIF
2886    ELSE
2887       contfrac_glo(:) = 1.0
2888    ENDIF
2889    !
2890    ! Read our test variable
2891    !
2892    iret = NF90_INQUIRE_VARIABLE(ifile, testvar_id, ndims=it)
2893    IF ( it == 2 ) THEN
2894       start = (/1,1,0/)
2895       count = (/nbpoint_glo,1,0/)
2896    ELSE IF ( it == 3 ) THEN
2897       start = (/1,1,1/)
2898       count = (/nbpoint_glo,1,1/)
2899    ELSE
2900       CALL ipslerr (3,'forcing_contfrac1d',"Testvar has a dimension larger than 3.", &
2901            "We do not know how to handle this", " ")
2902    ENDIF
2903    iret = NF90_GET_VAR(ifile, testvar_id, testvar, start, count)
2904    IF (iret /= NF90_NOERR) THEN
2905       WRITE(*,*) TRIM(nf90_strerror(iret))
2906       CALL ipslerr (3,'forcing_contfrac1d',"Error reading testvar.", " ", " ")
2907    ENDIF
2908    !
2909  END SUBROUTINE forcing_contfrac1d
2910
2911!!  =============================================================================================================================
2912!! SUBROUTINE: forcing_contfrac2d
2913!!
2914!>\BRIEF     
2915!!
2916!! DESCRIPTION: This routine build the land/mask if needed and gets the contfrac variable from forcing file.
2917!!              Here we treat the case where the variables is 2D. Thus we also need to identify the land points.
2918!!              For this we can either use the contfrac variable or the test variable.   
2919!!
2920!! \n
2921!_ ==============================================================================================================================
2922
2923  SUBROUTINE forcing_contfrac2d(ifile, testvar_id, contfrac_id, testvar, contfrac, testvar_missing, contfrac_missing, nbland)
2924    !
2925    ! This routine build the land/mask if needed and gets the contfrac variable from forcing file.
2926    ! Here we treat the case where the variables is 2D. Thus we also need to identify the land points.
2927    ! For this we can either use the contfrac variable or the test variable.
2928    !
2929    ! ARGUMENTS
2930    !
2931    INTEGER(i_std), INTENT(in)                 :: ifile
2932    INTEGER(i_std), INTENT(in)                 :: testvar_id, contfrac_id
2933    REAL(r_std), DIMENSION(:,:), INTENT(inout) :: testvar 
2934    REAL(r_std), DIMENSION(:,:), INTENT(inout) :: contfrac
2935    REAL(r_std), INTENT(in)                    :: testvar_missing
2936    REAL(r_std), INTENT(in)                    :: contfrac_missing
2937    INTEGER(i_std), INTENT(out)                :: nbland
2938    !
2939    ! LOCAL
2940    !
2941    INTEGER(i_std)                           :: i, j, it, iret
2942    INTEGER(i_std), DIMENSION(4)             :: start, count
2943    !
2944    !
2945    nbland = 0
2946    !
2947    IF ( contfrac_id > 0 ) THEN
2948       !
2949       iret = NF90_INQUIRE_VARIABLE(ifile, contfrac_id, ndims=it)
2950       IF ( it == 2 ) THEN
2951          start = (/1,1,0,0/)
2952          count = (/iim_glo,jjm_glo,0,0/)
2953          iret = NF90_GET_VAR(ifile, contfrac_id, contfrac, start, count)
2954          IF (iret /= NF90_NOERR) THEN
2955             WRITE(*,*) TRIM(nf90_strerror(iret))
2956             CALL ipslerr (3,'forcing_contfrac2d',"Error reading contfrac.", " ", " ")
2957          ENDIF
2958       ELSE
2959          CALL ipslerr (3,'forcing_contfrac2d',"Contfrac has a dimension different of 1.", &
2960               "We do not know how to handle this.", " ")
2961       ENDIF
2962       
2963       IF ( MAXVAL(contfrac) >= contfrac_missing ) THEN
2964          ! We have missing values in contfrac and we use it to count number of land points
2965          DO i=1,iim_glo
2966             DO j=1,jjm_glo
2967                IF ( contfrac(i,j) /= contfrac_missing .AND. contfrac(i,j) > 0.0 ) THEN
2968                   nbland = nbland + 1
2969                ENDIF
2970             ENDDO
2971          ENDDO
2972         
2973       ELSE
2974          ! Then ocean is fully contfrc=0 !
2975          DO i=1,iim_glo
2976             DO j=1,jjm_glo
2977                IF ( contfrac(i,j) > 0.0 ) THEN
2978                   nbland = nbland + 1
2979                ENDIF
2980             ENDDO
2981          ENDDO
2982         
2983       ENDIF
2984
2985       ! If we did not find any land points on the map (i.e. iim_glo > 1 and jjm_glo > 1) then we
2986       ! look for fractions larger then 0.
2987       !
2988       IF ( iim_glo > 1 .AND. jjm_glo > 1 .AND. nbland < 1 ) THEN
2989          DO i=1,iim_glo
2990             DO j=1,jjm_glo
2991                IF ( contfrac(i,j) > 0.0 ) THEN
2992                   nbland = nbland + 1
2993                ENDIF
2994             ENDDO
2995          ENDDO
2996       ENDIF
2997       
2998       ! Did we get a result ?
2999       !
3000       IF ( iim_glo > 1 .AND. jjm_glo > 1 .AND. nbland < 1 ) THEN
3001          CALL ipslerr (3,'forcing_contfrac2d',"Contfrac was used to count the number of land points.", &
3002               &          "We still have not found any land points when we looked for contfrac > 0.", " ")
3003       ENDIF
3004       !
3005    ELSE
3006       ! Just so that we have no un-initialized variable
3007       contfrac(:,:) = 0.0
3008    ENDIF
3009    !
3010    IF ( testvar_id > 0 ) THEN
3011       !
3012       iret = NF90_INQUIRE_VARIABLE(ifile, testvar_id, ndims=it)
3013       IF ( it == 2 ) THEN
3014          start = (/1,1,0,0/)
3015          count = (/iim_glo,jjm_glo,0,0/)
3016       ELSE IF ( it == 3 ) THEN
3017          start = (/1,1,1,0/)
3018          count = (/iim_glo,jjm_glo,1,0/)
3019       ELSE IF ( it == 4 ) THEN
3020          start = (/1,1,1,1/)
3021          count = (/iim_glo,jjm_glo,1,1/)
3022       ELSE
3023          CALL ipslerr (3,'forcing_contfrac2d',"testvar has a dimension of 1 or larger than 4.", &
3024               "We do not know how to handle this.", " ")
3025       ENDIF
3026       iret = NF90_GET_VAR(ifile, testvar_id, testvar, start, count)
3027       IF (iret /= NF90_NOERR) THEN
3028          WRITE(*,*) TRIM(nf90_strerror(iret))
3029          CALL ipslerr (3,'forcing_contfrac2d',"Error reading testvar.", " ", " ")
3030       ENDIF
3031       !
3032       ! IF with count frac we did not get the number of land points, let us try it here
3033       !
3034       IF ( nbland < 1 ) THEN
3035          DO i=1,iim_glo
3036             DO j=1,jjm_glo
3037                IF ( testvar(i,j) < testvar_missing ) THEN
3038                   nbland = nbland + 1
3039                   ! Add infor to contfrac
3040                   IF ( contfrac_id < 0 ) THEN
3041                      contfrac(i,j) = 1.0
3042                   ENDIF
3043                ENDIF
3044             ENDDO
3045          ENDDO
3046       ENDIF
3047       !
3048       !
3049       ! Did we get a result here ?
3050       !
3051       IF ( iim_glo > 1 .AND. jjm_glo > 1 .AND. nbland < 1 ) THEN
3052          CALL ipslerr (3,'forcing_contfrac2d',"Contfrac and testvar were used to count the number", &
3053               &          "of land points. We have not found any land points.", " ")
3054       ENDIF
3055       !
3056    ENDIF
3057    !
3058  END SUBROUTINE forcing_contfrac2d
3059
3060!!  =============================================================================================================================
3061!! SUBROUTINE: forcing_checkindex
3062!!
3063!>\BRIEF     
3064!!
3065!! DESCRIPTION:  For ORCHIDEE its paralelisation requires that the land points are ordered
3066!!               in such a way that the longitude runs fastest. That means that we go over the
3067!!               globle filling one line after the other.
3068!!               As this might not be the case in a compressed vector of land points, we need to
3069!!               put all the points on the 2d grid and then scan them in the right order.
3070!!               The reindexing is prepared here. 
3071!!
3072!! \n
3073!_ ==============================================================================================================================
3074
3075  SUBROUTINE forcing_checkindex(dump_mask, testvarname, testvar)
3076    !
3077    ! For ORCHIDEE its paralelisation requires that the land points are ordered
3078    ! in such a way that the longitude runs fastest. That means that we go over the
3079    ! globle filling one line after the other.
3080    ! As this might not be the case in a compressed vector of land points, we need to
3081    ! put all the points on the 2d grid and then scan them in the right order.
3082    ! The reindexing is prepared here.
3083    !
3084    LOGICAL          :: dump_mask
3085    CHARACTER(LEN=*) :: testvarname
3086    REAL(r_std)      :: testvar(:)
3087    !
3088    INTEGER(i_std) :: j, i, ik
3089    REAL(r_std), ALLOCATABLE, DIMENSION(:)   :: testvar_reind
3090    !
3091    !
3092    !
3093    ! Get the indices of the land points in the focing file
3094    !
3095    IF ( .NOT. ALLOCATED(reindex_glo)) ALLOCATE(reindex_glo(nbpoint_glo))
3096    IF ( .NOT. ALLOCATED(origind)) ALLOCATE(origind(iim_glo,jjm_glo))
3097    !
3098    ! Find the origine of each point in the gathered vector on the xy grid.
3099    !
3100    origind(:,:) = -1
3101    mask_glo(:,:) = 0
3102    DO ik=1,nbpoint_glo
3103       j = INT((lindex_glo(ik)-1)/iim_glo)+1
3104       i = (lindex_glo(ik)-(j-1)*iim_glo)
3105       origind(i,j) = ik
3106       mask_glo(i,j) = 1
3107    ENDDO
3108    !
3109    ! Prepare a reindexing array so that the vector goes in the right order : longitude runs
3110    ! faster than the latitude. Put then also the right information into lindex_glo.
3111    !
3112    ik=0
3113    DO j=1,jjm_glo
3114       DO i=1,iim_glo
3115          IF ( origind(i,j) > zero ) THEN
3116             ik = ik+1
3117             reindex_glo(ik) = origind(i,j)
3118             lindex_glo(ik) = (j-1)*iim_glo+i
3119          ENDIF
3120       ENDDO
3121    ENDDO
3122    !
3123    !
3124    ! Write the mask and a test variable to a file so that the user can check that all is OK
3125    !
3126    IF ( dump_mask) THEN
3127       !
3128       ! Scatter the test variable and save it in the file
3129       !
3130       WRITE(*,*) MINVAL(testvar), "<< test variable ", TRIM(testvarname), " <<", MAXVAL(testvar) 
3131       ALLOCATE(testvar_reind(nbpoint_glo))
3132       !
3133       CALL forcing_reindex(nbpoint_glo, testvar, nbpoint_glo, testvar_reind, reindex_glo)
3134       !
3135       
3136       CALL forcing_writetestvar("forcing_mask_glo.nc", iim_glo, jjm_glo, nbpoint_glo, &
3137            &                    lon_glo(:,1), lat_glo(1,:), lindex_glo, mask_glo, &
3138            &                    testvarname, testvar_reind)
3139       !
3140    ENDIF
3141    !
3142    ! Clean up !
3143    !
3144    IF ( ALLOCATED(testvar_reind) ) DEALLOCATE(testvar_reind)
3145    !
3146  END SUBROUTINE forcing_checkindex
3147
3148!!  =============================================================================================================================
3149!! SUBROUTINE: forcing_writetestvar
3150!!
3151!>\BRIEF Write the mask and a test variable to a netCDF file.     
3152!!
3153!! DESCRIPTION: This routine allows to test if the variables read from the forcing files is well read.
3154!!              Typically the forcing is compressed by gathering and thus it is a safe practice
3155!!              to verify that the un-compression is done correctly and that all points end-up in the
3156!!              right place on the global lat/lon grid.
3157!!
3158!! \n
3159!_ ==============================================================================================================================
3160
3161  SUBROUTINE forcing_writetestvar(ncdffile, iim, jjm, nbland, lon, lat, lindex, mask, varname, var)
3162    !
3163    ! Write the mask and a test variable to a netCDF file
3164    !
3165    ! ARGUMENTS
3166    !
3167    CHARACTER(LEN=*), INTENT(in) :: ncdffile
3168    INTEGER(i_std), INTENT(in)   :: iim, jjm, nbland
3169    REAL(r_std), INTENT(in)      :: lon(iim), lat(jjm)
3170    INTEGER(i_std), INTENT(in)   :: lindex(nbland)
3171    INTEGER(i_std), INTENT(in)   :: mask(iim,jjm)
3172    CHARACTER(LEN=*), INTENT(in) :: varname
3173    REAL(r_std), INTENT(in)      :: var(nbland)
3174    !
3175    ! Local
3176    !
3177    INTEGER(i_std) :: ik, i, j
3178    INTEGER(i_std) :: iret, nlonid, nlatid, varid, fid, ierr, iland
3179    INTEGER(i_std) :: testid
3180    INTEGER(i_std), DIMENSION(2) :: lolaid
3181    REAL(r_std) :: test_scat(iim,jjm)
3182    !
3183    !
3184    test_scat(:,:) = NF90_FILL_REAL
3185    CALL forcing_reindex(nbland, var, iim, jjm, test_scat, lindex)
3186    !
3187    iret = NF90_CREATE(ncdffile, NF90_WRITE, fid)
3188    IF (iret /= NF90_NOERR) THEN
3189       CALL ipslerr (3,'forcing_writetestvar',"Error opening the output file : ", ncdffile, " ")
3190    ENDIF
3191    !
3192    ! Define dimensions
3193    !
3194    iret = NF90_DEF_DIM(fid,'lon',iim,lolaid(1))
3195    iret = NF90_DEF_DIM(fid,'lat',jjm,lolaid(2))
3196    !
3197    !
3198    iret = NF90_DEF_VAR(fid,"lon",NF90_REAL4,lolaid(1),nlonid)
3199    iret = NF90_PUT_ATT(fid,nlonid,'axis',"X")
3200    iret = NF90_PUT_ATT(fid,nlonid,'standard_name',"longitude")
3201    iret = NF90_PUT_ATT(fid,nlonid,'units',"degree_east")
3202    iret = NF90_PUT_ATT(fid,nlonid,'valid_min',MINVAL(lon_glo))
3203    iret = NF90_PUT_ATT(fid,nlonid,'valid_max',MAXVAL(lon_glo))
3204    iret = NF90_PUT_ATT(fid,nlonid,'long_name',"Longitude")
3205    !
3206    iret = NF90_DEF_VAR(fid,"lat",NF90_REAL4,lolaid(2),nlatid)
3207    iret = NF90_PUT_ATT(fid,nlatid,'axis',"Y")
3208    iret = NF90_PUT_ATT(fid,nlatid,'standard_name',"latitude")
3209    iret = NF90_PUT_ATT(fid,nlatid,'units',"degree_north")
3210    iret = NF90_PUT_ATT(fid,nlatid,'valid_min',MINVAL(lat_glo))
3211    iret = NF90_PUT_ATT(fid,nlatid,'valid_max',MAXVAL(lat_glo))
3212    iret = NF90_PUT_ATT(fid,nlatid,'long_name',"Latitude")
3213    !
3214    iret = NF90_DEF_VAR(fid,"mask",NF90_REAL4,lolaid,varid)
3215    !
3216    iret = NF90_DEF_VAR(fid,TRIM(varname),NF90_REAL4,lolaid,testid)
3217    iret = NF90_PUT_ATT(fid,testid,'_FillValue',NF90_FILL_REAL)
3218    iret = NF90_PUT_ATT(fid,testid,'missing_value',NF90_FILL_REAL)
3219    !
3220    iret = NF90_ENDDEF (fid)
3221    IF (iret /= NF90_NOERR) THEN
3222       WRITE(*,*) TRIM(nf90_strerror(iret))
3223       CALL ipslerr (3,'forcing_writetestvar',"Error ending definitions in file : ", ncdffile, " ")
3224    ENDIF
3225    !
3226    ! Write variables
3227    !
3228    iret = NF90_PUT_VAR(fid,nlonid,lon)
3229    iret = NF90_PUT_VAR(fid,nlatid,lat)
3230    iret = NF90_PUT_VAR(fid,varid,REAL(mask))
3231    iret = NF90_PUT_VAR(fid,testid,test_scat)
3232    !
3233    ! Close file
3234    !
3235    iret = NF90_CLOSE(fid)
3236    IF (iret /= NF90_NOERR) THEN
3237       CALL ipslerr (3,'forcing_writetestvar',"Error closing the output file : ", ncdffile, " ")
3238    ENDIF
3239    !
3240  END SUBROUTINE forcing_writetestvar
3241
3242!!  =============================================================================================================================
3243!! SUBROUTINE: forcing_zoomgrid
3244!!
3245!>\BRIEF      We zoom into the region requested by the user.
3246!!
3247!! DESCRIPTION: Get the area to be zoomed and the sizes of arrays we will need.
3248!!              This subroutine fills the *_loc variables.
3249!!              If requested it will dump a test vraible into a netCDF file. 
3250!!
3251!! \n
3252!_ ==============================================================================================================================
3253
3254  SUBROUTINE forcing_zoomgrid (zoom_lon, zoom_lat, filename, dumpncdf)
3255    !
3256    ! Get the area to be zoomed and the sizes of arrays we will need.
3257    ! This subroutine fills the *_loc variables.
3258    ! If requested it will dump a test vraible into a netCDF file.
3259    !
3260    ! ARGUMENTS
3261    !
3262    REAL(r_std), DIMENSION(2), INTENT(in) :: zoom_lon, zoom_lat
3263    CHARACTER(LEN=*), INTENT(in) :: filename
3264    LOGICAL, INTENT(in) :: dumpncdf
3265    !
3266    ! LOCAL
3267    !
3268    INTEGER(i_std) :: i, j, ic, jc, ik, ig
3269    REAL(r_std) :: dx, dy, coslat
3270    REAL(r_std) :: lon_tmp(iim_glo), lat_tmp(jjm_glo)
3271    REAL(r_std) :: longlo_tmp(iim_glo,jjm_glo)
3272    REAL(r_std), ALLOCATABLE, DIMENSION(:) :: lon_val, lat_val
3273    INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:,:) :: zoom_index
3274    !
3275    INTEGER(i_std) :: iret, force_id, iv
3276    INTEGER(i_std), DIMENSION(1) :: imin, jmin
3277    INTEGER(i_std), DIMENSION(2) :: start, count
3278    INTEGER(i_std), DIMENSION(3) :: start2d, count2d
3279    REAL(r_std), ALLOCATABLE, DIMENSION(:) :: readvar, zoomedvar
3280     REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: readvar2d
3281    INTEGER(i_std), ALLOCATABLE, DIMENSION(:) :: index_glotoloc
3282    REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: lalo
3283    CHARACTER(LEN=8) :: testvarname="Tair"
3284    !
3285    ! 0.0 Verify we are on the root processor
3286    !
3287    IF ( .NOT. is_root_prc ) THEN
3288       CALL ipslerr (3,'forcing_zoomgrid'," This routine can only be called on the root processor.", " ", " ")
3289    ENDIF
3290    !
3291    ! 0.1 Inform the user
3292    !
3293    WRITE(*,*) "Zoom forcing : lon = ", zoom_lon
3294    WRITE(*,*) "Zoom forcing : lat = ", zoom_lat
3295    !
3296    ! Some forcing files have longitudes going from 0 to 360. This code works on the
3297    ! -180 to 180 longitude grid. So if needed we transform the longitudes of the global grid.
3298    !
3299    IF ( MAXVAL(lon_glo) <= 180.0 ) THEN
3300       longlo_tmp=lon_glo
3301    ELSE
3302       DO i=1,iim_glo
3303          DO j=1,jjm_glo
3304             IF ( lon_glo(i,j) <= 180.0 ) THEN
3305                longlo_tmp(i,j) = lon_glo(i,j)
3306             ELSE
3307                longlo_tmp(i,j) = lon_glo(i,j)-360
3308             ENDIF
3309          ENDDO
3310       ENDDO
3311    ENDIF
3312    !
3313    ! See if we need to zoom
3314    !
3315    IF (MINVAL(zoom_lon) > MINVAL(longlo_tmp) .OR. MAXVAL(zoom_lon) < MAXVAL(longlo_tmp) .OR.&
3316         & MINVAL(zoom_lat) > MINVAL(lat_glo) .OR. MAXVAL(zoom_lat) < MAXVAL(lat_glo) ) THEN
3317       zoom_forcing = .TRUE.
3318    ENDIF
3319    !
3320    ! Determine the size in x and y of the zoom
3321    !
3322    IF ( zoom_forcing ) THEN
3323       !
3324       ! Working with the hypothesis it is a regular global grid and bring it back to the -180 to 180 interval
3325       ! if needed.
3326       !
3327       lon_tmp(:) = longlo_tmp(:,1)
3328       lat_tmp(:) = lat_glo(1,:)
3329       !
3330       DO i=1,iim_glo
3331          IF ( lon_tmp(i) <= MINVAL(zoom_lon) .OR.  lon_tmp(i) >= MAXVAL(zoom_lon) ) THEN
3332             lon_tmp(i) = 0.0
3333          ELSE
3334             lon_tmp(i) = 1.0
3335          ENDIF
3336       ENDDO
3337       DO j=1,jjm_glo
3338          IF ( lat_tmp(j) <= MINVAL(zoom_lat) .OR. lat_tmp(j) >= MAXVAL(zoom_lat) ) THEN
3339             lat_tmp(j) = 0.0
3340          ELSE
3341             lat_tmp(j) = 1.0
3342          ENDIF
3343       ENDDO
3344       iim_loc = NINT(SUM(lon_tmp))
3345       jjm_loc = NINT(SUM(lat_tmp))
3346    ELSE
3347       iim_loc = iim_glo
3348       jjm_loc = jjm_glo
3349       lon_tmp(:) = 1.0
3350       lat_tmp(:) = 1.0
3351    ENDIF
3352    !
3353    ! Determine the number of land points in the zoom
3354    !
3355    IF ( .NOT. ALLOCATED(lon_loc) ) ALLOCATE(lon_loc(iim_loc,jjm_loc))
3356    IF ( .NOT. ALLOCATED(lat_loc) ) ALLOCATE(lat_loc(iim_loc,jjm_loc))
3357    IF ( .NOT. ALLOCATED(mask_loc) ) ALLOCATE(mask_loc(iim_loc,jjm_loc))
3358    IF ( .NOT. ALLOCATED(zoom_index) ) ALLOCATE(zoom_index(iim_loc,jjm_loc,2))
3359    !
3360    IF ( .NOT. ALLOCATED(lon_val)) ALLOCATE(lon_val(iim_loc))
3361    IF ( .NOT. ALLOCATED(lat_val)) ALLOCATE(lat_val(jjm_loc))
3362    !
3363    ! Create our new lat/lon system which is in the -180/180 range and South to North and West to East
3364    !
3365    ic=0
3366    DO i=1,iim_glo
3367       IF ( lon_tmp(i) > 0 ) THEN
3368          ic = ic+1
3369          lon_val(ic) = longlo_tmp(i,1)
3370       ENDIF
3371    ENDDO
3372    jc=0
3373    DO j=1,jjm_glo
3374       IF ( lat_tmp(j) > 0 ) THEN
3375          jc = jc+1
3376          lat_val(jc) = lat_glo(1,j)
3377       ENDIF
3378    ENDDO
3379    CALL sort(lon_val, iim_loc)
3380    CALL sort(lat_val, jjm_loc)
3381    !
3382    ! Now find the correspondance between the zoomed & re-ordered grid and the global one.
3383    !
3384    DO i=1,iim_loc
3385       DO j=1,jjm_loc
3386          !
3387          imin=MINLOC(ABS(longlo_tmp(:,1)-lon_val(i)))
3388          jmin=MINLOC(ABS(lat_glo(1,:)-lat_val(j)))
3389          !
3390          lon_loc(i,j) = longlo_tmp(imin(1),jmin(1))
3391          lat_loc(i,j) = lat_glo(imin(1),jmin(1))
3392          mask_loc(i,j) = mask_glo(imin(1),jmin(1))
3393          !
3394          zoom_index(i,j,1) = imin(1)
3395          zoom_index(i,j,2) = jmin(1)
3396          !
3397       ENDDO
3398    ENDDO
3399    !
3400    nbpoint_loc = SUM(mask_loc)
3401   
3402
3403    IF ( .NOT. zoom_forcing .AND. nbpoint_loc .NE. nbpoint_glo) THEN
3404       WRITE(*,*) "We have not zoomed into the forcing file still we get a number of"
3405       WRITE(*,*) "land points that is different from what we have in the forcing file."
3406       STOP "forcing_zoomgrid"
3407    ENDIF
3408    !
3409    IF ( .NOT. ALLOCATED(lindex_loc)) ALLOCATE(lindex_loc(nbpoint_loc))
3410    IF ( .NOT. ALLOCATED(reindex_loc)) ALLOCATE(reindex_loc(nbpoint_loc))
3411    IF ( .NOT. ALLOCATED(contfrac_loc)) ALLOCATE(contfrac_loc(nbpoint_loc))
3412    !
3413    IF ( .NOT. ALLOCATED(reindex2d_loc)) ALLOCATE(reindex2d_loc(nbpoint_loc,2))
3414    IF ( .NOT. ALLOCATED(index_glotoloc)) ALLOCATE(index_glotoloc(nbpoint_glo))
3415    IF ( .NOT. ALLOCATED(lalo)) ALLOCATE(lalo(nbpoint_loc,2))
3416    !
3417    ! Do the actual zoom on the grid
3418    !
3419    ! Set indices of all points as non existant so that we can fill in as we zoom the
3420    ! indices of the points which exist.
3421    index_glotoloc(:) = -1
3422    !
3423    ik = 0
3424    !
3425    ! Loop only over the zoomed grid
3426    !
3427    ! Why does the inner loop need to be ic for the pralalisation ????
3428    !
3429    DO jc=1,jjm_loc
3430       DO ic=1,iim_loc
3431          !
3432          ! Point back from the local to the original global i*j grid
3433          !
3434          i = zoom_index(ic,jc,1)
3435          j = zoom_index(ic,jc,2)
3436          !
3437          IF ( origind(i,j) > 0 ) THEN
3438             ik = ik+1
3439             ! index of the points in the local grid
3440             lindex_loc(ik) = (jc-1)*iim_loc+ic
3441             !
3442             ! For land points, the index of global grid is saved for the this point on the local grid
3443             reindex_loc(ik) = origind(i,j)
3444             !
3445             ! Keep also the i and j of the global grid for this land point on the local grid
3446             reindex2d_loc(ik,1) = i
3447             reindex2d_loc(ik,2) = j
3448             !
3449             ! Keep the reverse : on the global grid the location where we put the value of the local grid.
3450             index_glotoloc(origind(i,j)) = ik
3451             !
3452             contfrac_loc(ik) = contfrac_glo(origind(i,j))
3453             !
3454             lalo(ik,1) = lat_glo(i,j)
3455             lalo(ik,2) = longlo_tmp(i,j)
3456             !
3457          ENDIF
3458       ENDDO
3459    ENDDO
3460    !
3461    !
3462    nbland_loc = 0
3463    DO ik=1, SIZE(contfrac_loc)
3464       IF (contfrac_loc(ik) > 0.0) THEN
3465          nbland_loc = nbland_loc + 1.0
3466       ENDIF
3467    ENDDO
3468    !
3469    !
3470    ncdfstart = MINVAL(reindex_loc)
3471    reindex_loc(:) = reindex_loc(:)-ncdfstart+1
3472    ncdfcount =  MAXVAL(reindex_loc)
3473    !
3474    ! Compute the areas and the corners on the grid over which we will run ORCHIDEE.
3475    ! As this module is dedicated for regular lat/lon forcing we know that we can compute these
3476    ! variables without further worries.
3477    !
3478    IF ( .NOT. ALLOCATED(area_loc)) ALLOCATE(area_loc(iim_loc,jjm_loc))
3479    IF ( .NOT. ALLOCATED(corners_loc)) ALLOCATE(corners_loc(iim_loc,jjm_loc,4,2))
3480    !
3481    ! Treat first the longitudes
3482    !
3483    DO j=1,jjm_loc
3484       dx = zero
3485       DO i=1,iim_loc-1
3486          dx = dx+ABS(lon_loc(i,j)-lon_loc(i+1,j))
3487       ENDDO
3488       dx = dx/(iim_loc-1)
3489       DO i=1,iim_loc
3490          corners_loc(i,j,1,1) = lon_loc(i,j)-dx/2.0
3491          corners_loc(i,j,2,1) = lon_loc(i,j)+dx/2.0
3492          corners_loc(i,j,3,1) = lon_loc(i,j)+dx/2.0
3493          corners_loc(i,j,4,1) = lon_loc(i,j)-dx/2.0
3494       ENDDO
3495    ENDDO
3496    !
3497    ! Now treat the latitudes
3498    !
3499    DO i=1,iim_loc
3500       dy = zero
3501       DO j=1,jjm_loc-1
3502          dy = dy + ABS(lat_loc(i,j)-lat_loc(i,j+1))
3503       ENDDO
3504       dy = dy/(jjm_loc-1)
3505       DO j=1,jjm_loc
3506          corners_loc(i,j,1,2) = lat_loc(i,j)+dy/2.0
3507          corners_loc(i,j,2,2) = lat_loc(i,j)+dy/2.0
3508          corners_loc(i,j,3,2) = lat_loc(i,j)-dy/2.0
3509          corners_loc(i,j,4,2) = lat_loc(i,j)-dy/2.0
3510       ENDDO
3511    ENDDO
3512    !
3513    ! Compute the areas of the grid (using the simplification that the grid is regular in lon/lat).
3514    !
3515    DO i=1,iim_loc
3516       DO j=1,jjm_loc
3517          coslat = MAX( COS(lat_loc(i,j) * pi/180. ), mincos )
3518          dx = ABS(corners_loc(i,j,2,1) - corners_loc(i,j,1,1)) * pi/180. * R_Earth * coslat
3519          dy = ABS(corners_loc(i,j,1,2) - corners_loc(i,j,3,2)) * pi/180. * R_Earth
3520          area_loc(i,j) = dx*dy
3521       ENDDO
3522    ENDDO
3523    !
3524    ! If requested we read a variable, zoomin and dump the result into a test file.
3525    !
3526    IF ( dumpncdf ) THEN
3527       iret = NF90_OPEN (filename, NF90_NOWRITE, force_id)
3528       IF (iret /= NF90_NOERR) THEN
3529          CALL ipslerr (3,'forcing_zoomgrid',"Error opening the forcing file :", filename, " ")
3530       ENDIF
3531       !
3532       ALLOCATE(readvar(ncdfcount), readvar2d(iim_glo,jjm_glo), zoomedvar(nbpoint_loc))
3533       !
3534       iret = NF90_INQ_VARID(force_id, TRIM(testvarname), iv)
3535       IF (iret /= NF90_NOERR) THEN
3536          CALL ipslerr (3,'forcing_zoomgrid',"Could not find variable Tair in file."," "," ")
3537       ENDIF
3538
3539       IF ( compressed ) THEN
3540          !
3541          start(1) = ncdfstart
3542          start(2) = 1
3543          count(1) = ncdfcount
3544          count(2) = 1
3545          !
3546          iret = NF90_GET_VAR(force_id, iv, readvar, start, count)
3547          IF (iret /= NF90_NOERR) THEN
3548             CALL ipslerr (3,'forcing_zoomgrid',"Could not read compressed variable Tair from file."," "," ")
3549          ENDIF
3550          CALL forcing_reindex(ncdfcount, readvar, nbpoint_loc, zoomedvar, reindex_loc)
3551          !
3552       ELSE
3553          !
3554          start2d(1) = 1
3555          start2d(2) = 1
3556          start2d(3) = 1
3557          count2d(1) = iim_glo
3558          count2d(2) = jjm_glo
3559          count2d(3) = 1
3560          !
3561          iret = NF90_GET_VAR(force_id, iv, readvar2d, start2d, count2d)
3562          IF (iret /= NF90_NOERR) THEN
3563             CALL ipslerr (3,'forcing_zoomgrid',"Could not read variable Tair from file."," "," ")
3564          ENDIF
3565          CALL forcing_reindex(iim_glo, jjm_glo, readvar2d, nbpoint_loc, zoomedvar, reindex2d_loc)
3566          !
3567       ENDIF
3568       !
3569       CALL forcing_writetestvar("forcing_mask_loc.nc", iim_loc, jjm_loc, nbpoint_loc, &
3570            &                    lon_loc(:,1), lat_loc(1,:), lindex_loc, mask_loc, &
3571            &                    TRIM(testvarname), zoomedvar)
3572       !
3573    ENDIF
3574    !
3575    ! Clean up
3576    !
3577    IF ( ALLOCATED(readvar) ) DEALLOCATE(readvar)
3578    IF ( ALLOCATED(readvar2d) ) DEALLOCATE(readvar2d)
3579    IF ( ALLOCATED(zoomedvar) ) DEALLOCATE(zoomedvar)
3580    IF ( ALLOCATED(index_glotoloc) ) DEALLOCATE(index_glotoloc)
3581    IF ( ALLOCATED(lalo) ) DEALLOCATE(lalo)
3582    !
3583  END SUBROUTINE forcing_zoomgrid
3584
3585!!  =============================================================================================================================
3586!! SUBROUTINE: forcing_givegridsize
3587!!
3588!>\BRIEF      Routine which exports the size of the grid on which the model will run, i.e. the zoomed grid.
3589!!
3590!! DESCRIPTION: This is needed to transfer the grid information from this module to the glogrid.f90 module. 
3591!!
3592!! \n
3593!_ ==============================================================================================================================
3594
3595  SUBROUTINE forcing_givegridsize (iim, jjm, nblindex)
3596    !
3597    ! Provides the size of the grid to be used to the calling program
3598    !
3599    ! Size of the x and y direction of the zoomed area
3600    INTEGER(i_std), INTENT(out) :: iim, jjm
3601    ! Number of land points in the zoomed area
3602    INTEGER(i_std), INTENT(out) :: nblindex
3603    !
3604    IF ( .NOT. is_root_prc ) THEN
3605       CALL ipslerr (3,'forcing_givegridsize'," This routine can only be called on the root processor.", &
3606            &          "The information requested is only available on root processor.", " ")
3607    ENDIF
3608    !
3609    iim = iim_loc
3610    jjm = jjm_loc
3611    nblindex = nbland_loc
3612    !
3613  END SUBROUTINE forcing_givegridsize
3614
3615!!  =============================================================================================================================
3616!! SUBROUTINE: forcing_
3617!!
3618!>\BRIEF     
3619!!
3620!! DESCRIPTION:   
3621!!
3622!! \n
3623!_ ==============================================================================================================================
3624
3625  SUBROUTINE forcing_vertical(force_id)
3626    !
3627    !- This subroutine explores the forcing file to decide what information is available
3628    !- on the vertical coordinate.
3629    !
3630    INTEGER, INTENT(IN) :: force_id
3631    !
3632    INTEGER(i_std) :: iret, ireta, iretb
3633    !
3634    INTEGER(i_std) :: sigma_id = -1, sigma_uv_id = -1
3635    INTEGER(i_std) :: hybsiga_id = -1, hybsiga_uv_id = -1
3636    INTEGER(i_std) :: hybsigb_id = -1, hybsigb_uv_id = -1
3637    INTEGER(i_std) :: levels_id = -1, levels_uv_id = -1
3638    INTEGER(i_std) :: height_id = -1, height_uv_id = -1
3639    INTEGER(i_std) :: lev_id = -1
3640    !
3641    LOGICAL :: var_exists, vara_exists, varb_exists, varuv_exists
3642    LOGICAL :: foundvar
3643    LOGICAL :: levlegacy
3644    !
3645    !- Set all the defaults
3646    !
3647    zfixed=.FALSE.
3648    zsigma=.FALSE.
3649    zhybrid=.FALSE.
3650    zlevels=.FALSE.
3651    zheight=.FALSE.
3652    zsamelev_uv = .TRUE.
3653    levlegacy = .FALSE.
3654    !
3655    foundvar = .FALSE.
3656    !
3657    !- We have a forcing file to explore so let us see if we find any of the conventions
3658    !- which allow us to find the height of T,Q,U and V.
3659    !
3660    IF ( force_id > 0 ) THEN
3661       !
3662       ! Case for sigma levels
3663       !
3664       IF ( .NOT. foundvar ) THEN
3665          ireta = NF90_INQ_VARID(force_id, 'Sigma', sigma_id)
3666          IF ( (sigma_id >= 0) .AND. (ireta == NF90_NOERR) ) THEN
3667             foundvar = .TRUE.
3668             zsigma = .TRUE.
3669             iretb = NF90_INQ_VARID(force_id, 'Sigma_uv', sigma_uv_id)
3670             IF ( (sigma_uv_id >= 0) .OR. (iretb == NF90_NOERR) ) zsamelev_uv = .FALSE.
3671          ENDIF
3672       ENDIF
3673       !
3674       ! Case for Hybrid levels
3675       !
3676       IF ( .NOT. foundvar ) THEN
3677          var_exists = .FALSE.
3678          varuv_exists = .FALSE.
3679          ireta = NF90_INQ_VARID(force_id, 'HybSigA', hybsiga_id)
3680          IF ( (hybsiga_id >= 0 ) .AND. (ireta == NF90_NOERR) ) THEN
3681             iretb = NF90_INQ_VARID(force_id, 'HybSigB', hybsigb_id)
3682             IF ( (hybsigb_id >= 0 ) .AND. (iretb == NF90_NOERR) ) THEN
3683                var_exists=.TRUE.
3684             ELSE
3685                CALL ipslerr ( 3, 'forcing_vertical','Missing the B coefficient for', &
3686                     &         'Hybrid vertical levels for T and Q','stop')
3687             ENDIF
3688          ENDIF
3689          ireta = NF90_INQ_VARID(force_id, 'HybSigA_uv', hybsiga_uv_id)
3690          IF ( (hybsiga_uv_id >= 0 ) .AND. (ireta == NF90_NOERR) ) THEN
3691             iretb = NF90_INQ_VARID(force_id, 'HybSigB_uv', hybsigb_uv_id)
3692             IF ( (hybsigb_uv_id >= 0 ) .AND. (iretb == NF90_NOERR) ) THEN
3693                varuv_exists=.TRUE.
3694             ELSE
3695                CALL ipslerr ( 3, 'forcing_vertical','Missing the B coefficient for', &
3696                     &         'Hybrid vertical levels for U and V','stop')
3697             ENDIF
3698          ENDIF
3699          IF ( var_exists ) THEN
3700             foundvar = .TRUE.
3701             zhybrid = .TRUE.
3702             IF ( varuv_exists ) zsamelev_uv = .FALSE.
3703          ENDIF
3704       ENDIF
3705       !
3706       ! Case for levels (i.e. a 2d time varying field with the height in meters)
3707       !
3708       IF ( .NOT. foundvar ) THEN
3709          ireta = NF90_INQ_VARID(force_id, 'Levels', levels_id)
3710          IF ( (levels_id >= 0 ) .AND. (ireta == NF90_NOERR) ) THEN
3711             foundvar = .TRUE.
3712             zlevels = .TRUE.
3713             iretb = NF90_INQ_VARID(force_id, 'Levels_uv', levels_uv_id)
3714             IF ( (levels_uv_id >= 0 ) .AND. (iretb == NF90_NOERR) ) zsamelev_uv = .FALSE.
3715          ENDIF
3716       ENDIF
3717       !
3718       ! Case where a fixed height is provided in meters
3719       !
3720       IF ( .NOT. foundvar ) THEN
3721          ireta = NF90_INQ_VARID(force_id, 'Height_Lev1', height_id)
3722          IF ( (height_id >= 0 ) .AND. (ireta == NF90_NOERR) ) THEN
3723             foundvar = .TRUE.
3724             zheight = .TRUE.       
3725             iretb = NF90_INQ_VARID(force_id, 'Height_Levuv', height_uv_id)
3726             IF ( (height_uv_id >= 0 ) .AND. (iretb == NF90_NOERR) ) zsamelev_uv = .FALSE.
3727          ENDIF
3728       ENDIF
3729       !
3730       ! Case where a fixed height is provided in meters in the lev variable
3731       !
3732       IF ( .NOT. foundvar ) THEN
3733          ireta = NF90_INQ_VARID(force_id, 'lev', lev_id)
3734          IF ( (lev_id >= 0 ) .AND. (ireta == NF90_NOERR) ) THEN
3735             foundvar = .TRUE.
3736             zheight = .TRUE.
3737             levlegacy = .TRUE.
3738          ENDIF
3739       ENDIF
3740       !
3741    ENDIF
3742    !
3743    ! We found forcing variables so we need to extract the values if we are dealing with constant values (i.e. all
3744    ! except the case zlevels
3745    !
3746    IF ( foundvar .AND. .NOT. zlevels ) THEN
3747       !
3748       IF ( zheight ) THEN
3749          !
3750          ! Constant height
3751          !
3752          IF ( levlegacy ) THEN
3753             iret = NF90_GET_VAR(force_id, lev_id, zlev_fixed)
3754             IF ( iret /= NF90_NOERR ) THEN
3755                CALL ipslerr ( 3, 'forcing_vertical','Attempted to read variable lev from forcing file in legacy mode', &
3756                     &         'NF90_GET_VAR failed.','stop')
3757             ENDIF
3758          ELSE
3759             iret = NF90_GET_VAR(force_id, height_id, zlev_fixed)
3760             IF ( iret /= NF90_NOERR ) THEN
3761                CALL ipslerr ( 3, 'forcing_vertical','Attempted to read variable Height_Lev1 from forcing file', &
3762                     &         'NF90_GET_VAR failed.','stop')
3763             ENDIF
3764             IF ( .NOT. zsamelev_uv ) THEN
3765                iret = NF90_GET_VAR(force_id, height_uv_id, zlevuv_fixed)
3766                IF ( iret /= NF90_NOERR ) THEN
3767                   CALL ipslerr ( 3, 'forcing_vertical','Attempted to read variable Height_Levuv from forcing file', &
3768                        &         'NF90_GET_VAR failed.','stop')
3769                ENDIF
3770             ENDIF
3771          ENDIF
3772          WRITE(*,*) "forcing_vertical : case ZLEV : Read from forcing file :", zlev_fixed, zlevuv_fixed
3773          !
3774       ELSE IF ( zsigma .OR. zhybrid ) THEN
3775          !
3776          ! Sigma or hybrid levels
3777          !
3778          IF ( zsigma ) THEN
3779             iret = NF90_GET_VAR(force_id, sigma_id, zhybrid_b)
3780             zhybrid_a = zero
3781             IF ( .NOT. zsamelev_uv ) THEN
3782                iret = NF90_GET_VAR(force_id, sigma_uv_id, zhybriduv_b)
3783                zhybriduv_a = zero
3784             ENDIF
3785          ELSE
3786             ireta = NF90_GET_VAR(force_id, hybsigb_id, zhybrid_b)
3787             iretb = NF90_GET_VAR(force_id, hybsiga_id, zhybrid_a)
3788             IF ( ireta /= NF90_NOERR .OR. iretb /= NF90_NOERR) THEN
3789                CALL ipslerr ( 3, 'forcing_vertical','Attempted to read variable HybSigA and HybSigB from forcing file', &
3790                     &         'NF90_GET_VAR failed.','stop')
3791             ENDIF
3792             IF ( .NOT. zsamelev_uv ) THEN
3793                ireta = NF90_GET_VAR(force_id, hybsigb_uv_id, zhybriduv_b)
3794                iretb = NF90_GET_VAR(force_id, hybsiga_uv_id, zhybriduv_a)
3795                IF ( ireta /= NF90_NOERR .OR. iretb /= NF90_NOERR) THEN
3796                   CALL ipslerr ( 3, 'forcing_vertical',&
3797                        &        'Attempted to read variable HybSigA_uv and HybSigB_uv from forcing file', &
3798                        &        'NF90_GET_VAR failed.','stop')
3799                ENDIF
3800             ENDIF
3801          ENDIF
3802          WRITE(*,*) "forcing_vertical : case Pressure coordinates : "
3803          WRITE(*,*) "Read from forcing file :", zhybrid_b, zhybrid_a, zhybriduv_b, zhybriduv_a
3804       ELSE
3805          !
3806          ! Why are we here ???
3807          !
3808          CALL ipslerr ( 3, 'forcing_vertical','What is the option used to describe the height of', &
3809               &         'the atmospheric forcing ?','Please check your forcing file.')
3810       ENDIF
3811    ENDIF
3812    !
3813    !- We have no forcing file to explore or we did not find anything. So revert back to the run.def and
3814    !- read what has been specified by the user.
3815    !
3816    IF ( force_id < 0 .OR. .NOT. foundvar ) THEN
3817       !
3818       !-
3819       !Config  Key  = HEIGHT_LEV1
3820       !Config  Desc = Height at which T and Q are given
3821       !Config  Def  = 2.0
3822       !Config  Help = The atmospheric variables (temperature and specific
3823       !Config         humidity) are measured at a specific level.
3824       !Config         The height of this level is needed to compute
3825       !Config         correctly the turbulent transfer coefficients.
3826       !Config         Look at the description of the forcing
3827       !Config         DATA for the correct value.
3828       !-
3829       zlev_fixed = 2.0
3830       CALL getin('HEIGHT_LEV1', zlev_fixed)
3831       !-
3832       !Config  Key  = HEIGHT_LEVW
3833       !Config  Desc = Height at which the wind is given
3834       !Config  Def  = 10.0
3835       !Config  Help = The height at which wind is needed to compute
3836       !Config         correctly the turbulent transfer coefficients.
3837       !-
3838       zlevuv_fixed = 10.0
3839       CALL getin('HEIGHT_LEVW', zlevuv_fixed)
3840
3841       zheight = .TRUE.
3842
3843       IF ( ABS(zlevuv_fixed-zlev_fixed) > EPSILON(zlev_fixed)) THEN
3844          zsamelev_uv = .FALSE.
3845       ENDIF
3846
3847       CALL ipslerr ( 2, 'forcing_vertical','The height of the atmospheric forcing variables', &
3848            &         'was not found in the netCDF file.','Thus the values in run.def were used ... or their defaults.')
3849    ENDIF
3850
3851  END SUBROUTINE forcing_vertical
3852
3853!!  =============================================================================================================================
3854!! SUBROUTINE: forcing_givegrid
3855!!
3856!>\BRIEF      Routine which exports the grid (longitude, latitude, land indices) on which the model will run, i.e. the zoomed grid.
3857!!
3858!! DESCRIPTION: This is needed to transfer the grid information from this module to the glogrid.f90 module. 
3859!!
3860!!
3861!! \n
3862!_ ==============================================================================================================================
3863
3864  SUBROUTINE forcing_givegrid (lon, lat, mask, area, corners, lindex, contfrac, calendar_tmp)
3865    !
3866    ! This subroutine will return to the caller the grid which has been extracted from the
3867    ! the forcing file. It is assumed that the caller has called forcing_givegridsize before
3868    ! and knows the dimensions of the fields and thus has done the correct allocations.
3869    !
3870    !
3871    REAL(r_std), INTENT(out) :: lon(iim_loc,jjm_loc), lat(iim_loc,jjm_loc)
3872    REAL(r_std), INTENT(out) :: mask(iim_loc,jjm_loc)
3873    REAL(r_std), INTENT(out) :: area(iim_loc,jjm_loc)
3874    REAL(r_std), INTENT(out) :: corners(iim_loc,jjm_loc,4,2)
3875    INTEGER(i_std), INTENT(out) :: lindex(nbpoint_loc)
3876    REAL(r_std), INTENT(out) :: contfrac(nbpoint_loc)
3877    CHARACTER(LEN=20), INTENT(out) :: calendar_tmp
3878    !
3879    IF ( .NOT. is_root_prc ) THEN
3880       CALL ipslerr (3,'forcing_givegrid'," This routine can only be called on the root processor.", &
3881            &          "The information requested is only available on root processor.", " ")
3882    ENDIF
3883    !
3884    IF (nbpoint_loc .NE. nbland_loc) THEN
3885       WRITE(numout, *) "forcing_givegrid:: nbpoint_loc=", nbpoint_loc
3886       WRITE(numout, *) "forcing_givegrid:: nbland_loc=", nbland_loc
3887       CALL ipslerr(3,'forcing_givegrid','nbpoint_loc and nbland_loc do match', & 
3888                    'The calculation of land points is not correct','')
3889    ENDIF
3890    !
3891    lon(:,:) = lon_loc(:,:)
3892    lat(:,:) = lat_loc(:,:)
3893    !
3894    mask(:,:) = mask_loc(:,:)
3895    area(:,:) = area_loc(:,:)
3896    corners(:,:,:,:) = corners_loc(:,:,:,:)
3897    !
3898    !
3899    lindex(:) = lindex_loc(:)
3900    contfrac(:) = contfrac_loc(:)
3901    !
3902    calendar_tmp = calendar
3903    !
3904  END SUBROUTINE forcing_givegrid
3905
3906!!  =============================================================================================================================
3907!! SUBROUTINE: forcing_checkdim
3908!!
3909!>\BRIEF     
3910!!
3911!! DESCRIPTION: Save the dimension or check that it is equal to the previous value.
3912!!              Should one of the spatial dimensions be different between 2 files, then we have a big problem.
3913!!              They probably do not belong to the same set of forcing files. 
3914!!
3915!! \n
3916!_ ==============================================================================================================================
3917
3918SUBROUTINE forcing_checkdim(ifile, filenames, out_dim, out_id, in_dim, in_id)
3919  !
3920  ! Save the dimension or check that it is equal to the previous value.
3921  ! Should one of the spatial dimensions be different between 2 files, then we have a big problem.
3922  ! They probably do not belong to the same set of forcing files.
3923  !
3924  INTEGER(i_std), INTENT(in) :: ifile
3925  CHARACTER(LEN=*), INTENT(in) :: filenames(:)
3926  INTEGER(i_std), INTENT(out) :: out_dim, out_id
3927  INTEGER(i_std), INTENT(in) :: in_dim, in_id
3928  !
3929  IF ( ifile == 1 ) THEN
3930     out_dim = in_dim
3931     out_id = in_id
3932  ELSE
3933     IF ( out_dim /= in_dim ) THEN
3934        CALL ipslerr (3,'forcing_ocheckdim', 'The dimension of the file is not the same of the first file opened.', &
3935             &        'The offending file is : ', filenames(ifile))
3936     ENDIF
3937  ENDIF
3938  !
3939END SUBROUTINE forcing_checkdim
3940
3941!!  =============================================================================================================================
3942!! SUBROUTINE: forcing_time
3943!!
3944!>\BRIEF Read the time from each file and create the time axis to be the basis for the simulation.     
3945!!
3946!! DESCRIPTION: This is an important routine which analyses the time axis of the forcing files and
3947!!              stores the main information in the SAVED variables of this routine.
3948!!              As this module manages a list of forcing files we also need to check that the time
3949!!              axis of all these files is continuous and homogeneous.
3950!!              The bounds are also build for all the time axes so that we know how to interpret the
3951!!              various variables.
3952!!
3953!! \n
3954!_ ==============================================================================================================================
3955
3956SUBROUTINE forcing_time(nbfiles, filenames)
3957  !
3958  ! Read the time from each file and create the time axis to be the basis
3959  ! for the simulation.
3960  !
3961  INTEGER(i_std) :: nbfiles
3962  CHARACTER(LEN=*) :: filenames(nbfiles)
3963  !
3964  INTEGER(i_std) :: iv, it, iff, tcnt, itbase, itbasetmp, ittmp
3965  INTEGER(i_std) :: tstart, maxtime_infile
3966  REAL(r_std), ALLOCATABLE, DIMENSION(:) :: timeint, time_read
3967  REAL(r_std), ALLOCATABLE, DIMENSION(:)   :: time_infiles
3968  CHARACTER(LEN=20) :: axname, calendar, timevarname
3969  CHARACTER(LEN=60) :: timestamp, tmpatt, ymd, hms
3970  INTEGER(i_std) :: tncstart(3), tnccount(3)
3971  !
3972  INTEGER(i_std) :: iret, year0, month0, day0, hours0, minutes0, seci
3973  INTEGER(i_std), DIMENSION(1) :: imax, imin
3974  REAL(r_std) :: sec0, date_int, date0_tmp
3975  CHARACTER :: strc
3976  LOGICAL :: check=.FALSE.
3977  !
3978  ! Check that we are working on the root proc.
3979  !
3980  IF ( .NOT. is_root_prc) THEN
3981     CALL ipslerr (3,'forcing_time',"Cannot run this routine o other procs than root.",&
3982          &        "All the information on the forcing files is only lated on the root processor."," ")
3983  ENDIF
3984  !
3985  ! Size of unlimited dimension added up through the files. If variable not allocated before by another
3986  ! subroutine, it needs to be done here.
3987  !
3988  IF ( .NOT. ALLOCATED(nbtime_perfile) ) ALLOCATE(nbtime_perfile(nbfiles))
3989  IF ( .NOT. ALLOCATED(date0_file) ) ALLOCATE(date0_file(nbfiles,nbtax))
3990  !
3991  ! Go through all files in the list in order to get the total number of time steps we have
3992  ! in the nbfiles files to be read
3993  !
3994  nb_forcing_steps = 0
3995  maxtime_infile = 0
3996  DO iff=1,nbfiles
3997     !
3998     iret = NF90_INQUIRE_DIMENSION(force_id(iff), id_unlim(iff), name=axname, len=nbtime_perfile(iff))
3999     IF (iret /= NF90_NOERR) THEN
4000        CALL ipslerr (3,'forcing_time',"Could not get size of dimension of unlimited axis"," "," ")
4001     ENDIF
4002     nb_forcing_steps =  nb_forcing_steps + nbtime_perfile(iff)
4003     IF ( nbtime_perfile(iff) > maxtime_infile ) maxtime_infile = nbtime_perfile(iff)
4004  ENDDO
4005  !
4006  ! Allocate the variables needed with the time length just calculated.
4007  ! These variables are saved in the module
4008  !
4009  ALLOCATE(time_infiles(nb_forcing_steps))
4010  ALLOCATE(time(nb_forcing_steps, nbtax*nbtmethods), time_bounds(nb_forcing_steps,nbtax*nbtmethods,2))
4011  ALLOCATE(time_axename(nbtax*nbtmethods), time_cellmethod(nbtax*nbtmethods))
4012  ALLOCATE(preciptime(nb_forcing_steps))
4013  ALLOCATE(time_sourcefile(nb_forcing_steps))
4014  ALLOCATE(time_id(nb_forcing_steps, nbtax))
4015  ! Allocate local variables
4016  ALLOCATE(time_read(nb_forcing_steps))
4017  ALLOCATE(timeint(nb_forcing_steps))
4018  !
4019  ! Get through all variables to find time_id
4020  ! The key variables to filled up here are time (the time stamps read in the file) and
4021  ! time_bounds which give the validity interval for the variables.
4022  !
4023  tstart=0
4024  !
4025  IF ( check ) WRITE(*,*) "forcing_time : going through ", nbfiles, " files to get the time."
4026  !
4027  DO iff=1,nbfiles
4028     !
4029     time_id(iff,:)=-1
4030     !
4031     ! Go through the variables in the file and find the one which is a time axis.
4032     !
4033     tcnt=1
4034     DO iv=1,nvars(iff)
4035        iret = NF90_GET_ATT(force_id(iff), iv, "units", tmpatt)
4036        IF ( INDEX(lowercase(tmpatt),'seconds since') > 0) THEN
4037           time_id(iff,tcnt)=iv
4038           tcnt=tcnt+1
4039           convtosec(iff)=1.0
4040        ELSE IF ( INDEX(lowercase(tmpatt),'minutes since') > 0) THEN
4041           time_id(iff,tcnt)=iv
4042           tcnt=tcnt+1
4043           convtosec(iff)=60.0
4044        ELSE IF ( INDEX(lowercase(tmpatt),'hours since') > 0) THEN
4045           time_id(iff,tcnt)=iv
4046           tcnt=tcnt+1
4047           convtosec(iff)=3600.0
4048        ELSE IF ( INDEX(lowercase(tmpatt),'days since') > 0) THEN
4049           time_id(iff,tcnt)=iv
4050           tcnt=tcnt+1
4051           convtosec(iff)=one_day
4052        ENDIF
4053     ENDDO
4054     IF ( ANY(time_id(iff,:) < 0) ) THEN
4055        CALL ipslerr (3,'forcing_time',"Incorrect numer of time axes. A time axis is missing ",&
4056             &        "in file :", filenames(iff))
4057     ENDIF
4058     !
4059     IF ( check ) WRITE(*,*) "forcing_time : Looking at time axis for file ", force_id(iff)
4060     !
4061     ! Looping through the time axes and read them.
4062     !
4063     DO tcnt=1,nbtax
4064        !
4065        iret = NF90_INQUIRE_VARIABLE(force_id(iff), time_id(iff,tcnt), name=timevarname)
4066        IF ( check ) WRITE(*,*) "forcing_time : in ", iff, " found variable ", timevarname
4067        !
4068        ! Get units of time axis
4069        !
4070        iret = NF90_GET_ATT(force_id(iff), time_id(iff,tcnt), "units", timestamp) 
4071        IF ( check ) WRITE(*,*) "forcing_time : has time stamp ", timestamp
4072        !
4073        ! Transform the start date of the netCDF file into a julian date for the model
4074        !
4075        timestamp = TRIM(timestamp(INDEX(timestamp,'since')+6:LEN_TRIM(timestamp)))
4076        ymd=TRIM(timestamp(1:INDEX(timestamp,' ')))
4077        hms=TRIM(timestamp(INDEX(timestamp,' ')+1:LEN_TRIM(timestamp)))
4078        !
4079        ! First extral the information from the date string
4080        !
4081        READ(ymd(1:INDEX(ymd,'-')-1),'(I4)') year0
4082        ymd=TRIM(ymd(INDEX(ymd,'-')+1:LEN_TRIM(ymd)))
4083        READ(ymd(1:INDEX(ymd,'-')-1),'(I2)') month0
4084        ymd=TRIM(ymd(INDEX(ymd,'-')+1:LEN_TRIM(ymd)))
4085        READ(ymd,'(I2)') day0
4086        !
4087        ! Now extract from the time string
4088        !
4089        READ(hms(1:INDEX(hms,':')-1),'(I2)') hours0
4090        hms=TRIM(hms(INDEX(hms,':')+1:LEN_TRIM(hms)))
4091        READ(hms(1:INDEX(hms,':')-1),'(I2)') minutes0
4092        hms=TRIM(hms(INDEX(hms,':')+1:LEN_TRIM(hms)))
4093        READ(hms,'(I2)') seci
4094        !
4095        sec0 = hours0*3600. + minutes0*60. + seci
4096        CALL ymds2ju (year0, month0, day0, sec0, date0_tmp)
4097        date0_file(iff,tcnt) = date0_tmp
4098        !
4099        ! Now get the actual dates
4100        !
4101        tncstart(1) = 1
4102        tnccount(1) = nbtime_perfile(iff)
4103        IF ( check ) WRITE(*,*) "forcing_time : number of values read : ", tnccount(1)
4104        iret = NF90_GET_VAR(force_id(iff), time_id(iff,tcnt), time_read, tncstart, tnccount)
4105        IF (iret /= NF90_NOERR) THEN
4106           CALL ipslerr (3,'forcing_time',"An error occured while reading time from the file."," "," ")
4107        ENDIF
4108        !
4109        ! Convert the variable time from seconds since to julian days
4110        !
4111        DO it=1,nbtime_perfile(iff)
4112           !!time_infiles(tstart+it) = date0_file(iff,tcnt) + time_read(it)*convtosec(iff)/one_day
4113           IF ( convtosec(iff) < one_day ) THEN
4114              time_infiles(tstart+it) = date0_file(iff,tcnt) + time_read(it)*convtosec(iff)/one_day
4115           ELSE
4116              ! In the case of daily forcing the start time has to be 00UTC in Julian days.
4117              time_infiles(tstart+it) = date0_file(iff,tcnt) + INT(time_read(it))
4118           ENDIF
4119        ENDDO
4120        if ( check ) WRITE(*,*) "File ", iff, "goes from ",  time_infiles(tstart+1), " to ", &
4121             time_infiles(tstart+nbtime_perfile(iff))
4122        !
4123        ! Estimate the bounds as this information is not yet in the forcing file.
4124        !
4125        date_int = (time_infiles(tstart+nbtime_perfile(iff)) - time_infiles(tstart+1))/(nbtime_perfile(iff)-1)
4126        forcing_tstep_ave = date_int*one_day
4127        !
4128        ! If this is the first file then simply keep the name of the time axis. Else search the same name
4129        ! in what has already been read
4130        !
4131        IF ( iff == 1 ) THEN
4132           itbase = (tcnt-1)*nbtmethods
4133           time_axename(itbase+1:itbase+4) = timevarname
4134           time_cellmethod(itbase+1) = "reference"
4135           time_cellmethod(itbase+2) = "start"
4136           time_cellmethod(itbase+3) = "cent"
4137           time_cellmethod(itbase+4) = "end"
4138        ELSE
4139           !
4140           ! If this is not the first file then we need to find the correct axis to add to.
4141           ! All information have already been saved with the first file.
4142           !
4143           DO ittmp=1,nbtax
4144              itbasetmp=(ittmp-1)*nbtmethods
4145              IF ( time_axename(itbasetmp+1) == timevarname ) THEN
4146                 itbase = itbasetmp
4147              ENDIF
4148           ENDDO
4149
4150        ENDIF
4151        !
4152        !
4153        ! Keep for future usage the various information on the time axis we have just read. This time axis can
4154        ! be understood in 3 different ways and each variable might use a different cell method for this time
4155        ! axis.
4156        !
4157        ! time(:,(tcnt-1)*nbtmethods+1) : corresponds to the reference time axis as it has been read from the file
4158        ! time(:,(tcnt-1)*nbtmethods+2) : is the time axis with a cell method which place the value at the
4159        !                                beginning of the time interval
4160        ! time(:,(tcnt-1)*nbtmethods+3) : is the time axis corresponding to variables placed at the center of the
4161        !                                time interval
4162        ! time(:,(tcnt-1)*nbtmethods+4) : for variables put at the end of the time interval over which they aere
4163        !                                for instance averaged.
4164        !
4165        ! In variable time_cellmethod we will write the type of cell method as descirbed above so that the selection
4166        ! of the right axis for each variable can be made automaticaly.
4167        !
4168        ! We also keep the name of the time axis read in preparation of file where we might have to read more than one
4169        ! time axis.
4170        !
4171        DO it=tstart+1,tstart+nbtime_perfile(iff)
4172           !
4173           ! Reference time
4174           !
4175           time(it,itbase+1) = time_infiles(it)
4176           time_bounds(it,itbase+1,1) = time_infiles(it)-date_int/2.0
4177           time_bounds(it,itbase+1,2) = time_infiles(it)+date_int/2.0
4178           !
4179           ! Start cell method
4180           time(it,itbase+2) = time_infiles(it)+date_int/2.0
4181           time_bounds(it,itbase+2,1) = time_infiles(it)
4182           time_bounds(it,itbase+2,2) = time_infiles(it)+date_int
4183           !
4184           ! Centered cell method
4185           time(it,itbase+3) = time_infiles(it)
4186           time_bounds(it,itbase+3,1) = time_infiles(it)-date_int/2.0
4187           time_bounds(it,itbase+3,2) = time_infiles(it)+date_int/2.0
4188           !
4189           ! End cell method
4190           time(it,itbase+4) = time_infiles(it)-date_int/2.0
4191           time_bounds(it,itbase+4,1) = time_infiles(it)-date_int
4192           time_bounds(it,itbase+4,2) = time_infiles(it)
4193           !
4194        ENDDO
4195        !
4196        ! Keep the number of the file from which we read this time.
4197        !
4198        time_sourcefile(tstart+1:tstart+nbtime_perfile(iff))=iff
4199        !
4200        IF ( check ) WRITE(*,*) "forcing_time : finished file ", iff
4201        !
4202     ENDDO
4203     !
4204     ! Before moving to the next file advance the pointer in the time arrays.
4205     !
4206     tstart=tstart+nbtime_perfile(iff)
4207     !
4208  ENDDO
4209  !
4210  IF ( check ) WRITE(*,*) "forcing_time : All files have been processed"
4211  !
4212  ! Verify that the forcing comes in regular time intervals. If not, many of the
4213  ! interpolation schemes will fail.
4214  ! This is only done on the first time axis ... is it enough ?
4215  !
4216  DO ittmp=1,nbtax
4217     itbase=(ittmp-1)*nbtmethods
4218     !
4219     date_int = (time(nb_forcing_steps,itbase+1) - time(1,itbase+1))/(nb_forcing_steps-1)
4220     forcing_tstep_ave = date_int*one_day
4221     !
4222     timeint(:) = 0
4223     DO it=1, nb_forcing_steps-1
4224        timeint(it) = time(it+1,itbase+1)-time(it,itbase+1)
4225     ENDDO
4226     !
4227     IF (  MAXVAL(timeint(1:nb_forcing_steps-1)) > date_int+0.1*date_int .OR.&
4228          & MINVAL(timeint(1:nb_forcing_steps-1)) < date_int-0.1*date_int) THEN
4229        WRITE(*,*) "The time steping of the forcing files does not seem to be regular on axis",time_axename(itbase+1),":"
4230        WRITE(*,*) "Average time step : ", date_int, "days = ", date_int*one_day, "sec."
4231        imax = MAXLOC(timeint(1:nb_forcing_steps-1))
4232        imin = MINLOC(timeint(1:nb_forcing_steps-1))
4233        WRITE(*,*) "Maximum time step : ", MAXVAL(timeint(1:nb_forcing_steps-1)), " at ", imax(1)
4234        WRITE(*,*) "Minimum time step : ", MINVAL(timeint(1:nb_forcing_steps-1)), " at ", imin(1)
4235        WRITE(*,*) "++++ Values around Maximum"
4236        DO it=MAX(imax(1)-5,1),