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

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

As in r7065 and r76158 for the new driver. Runs OK on jean-zay with the old driver, but not yet tested with the new driver. The corresponding configurations need to be added following r7081, r7082, and r7719.

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