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

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

Makes all keywords related to SPREAD_PREC include an A in spreAd, cf ticket #730, comment 5 and on. All forms with SPRED_PREC get deprecated. The file orchidee.default is adapted, and now includes info on SPREAD_PREC_SEC and SPREAD_PREC_CONT (only used by new driver). The configurations, however, have not yet been updated (the only committed configuration affected is FG3nd with SPRED_PREC_SEC = 5400).

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