PROGRAM main
!---------------------------------------------------------------------
!- This PROGRAM is the driver of the dim2 version of SECHIBA. It's
!- main use is for PILPS type experiments. To run it you need to have
!- the following software :
!- - SECHIBA
!- - ioipsl
!- - netcdf
!- - F90 compiler
!- - tk (optional but very useful for configuring the model)
!- Obviously also the associated Makefiles.
!-
!- The forcing data needs to be in netCDF format and should
!- contain the following variables :
!- - Incoming SW radiation
!- - Incoming LW radiation
!- - Precipitation
!- - Air temperature at a reference level
!- - Air humidity at the same level
!- - wind at the same level
!- - surface pressure
!-
!- Once you have all this and compiled the model you can configure it.
!- Type make config and set all the values presented in the menu.
!- This tool will create a run.def which will be used by the model.
!- The run.def can also be edited by hand but it is more tedious.
!- Once this run.def is created the model is ready to run.
!-
!- $Header: /home/ssipsl/CVSREP/sechiba/dim2_driver.f90,v 1.38 2000/05/11 14:48:27 ssipsl Exp $
!---------------------------------------------------------------------
   USE ioipsl
   USE intersurf,  ONLY : intersurf_main
   USE readdim2
!-
   IMPLICIT NONE
!-
   REAL,PARAMETER :: cp_air = 1004.675
   REAL,PARAMETER :: cte_grav = 9.80665
   REAL,PARAMETER :: tp_00 = 273.15
   REAL,PARAMETER :: cte_molr = 287.05
!-
   REAL,SAVE  :: val_exp
!-
   INTEGER :: iim, jjm, llm
   INTEGER :: im, jm, lm, tm, tmp_tm, is, force_id, itest, jtest
   REAL :: dt, dt_force, dt_rest, dw, date0, date0_rest
   REAL :: zlflu, un_jour
   REAL :: un_an, relax, alpha
!-
   REAL, ALLOCATABLE, DIMENSION(:,:) :: &
   & swdown, precip_rain, precip_snow, tair_obs, u, v, qair_obs, pb, lwdown, &
   & eair_obs, zlev_vec, zlflu_vec
!- Variables which are forcings for SECHIBA
   REAL, ALLOCATABLE, DIMENSION(:,:) :: &
   & petAcoef, peqAcoef, petBcoef, peqBcoef, for_cdrag, &
   & for_u, for_v, for_swnet, for_swdown, for_lwdown, &
   & for_psurf, for_qair, for_tair, for_eair, &
   & for_ccanopy, for_rau, tmp_eair, tmp_tair, &
   & tmp_qair, tmp_pb
!-
   REAL, ALLOCATABLE, DIMENSION(:,:) :: &
   & old_qair, old_eair, tsol_rad, vevapp, &
   & temp_sol_NEW, temp_sol_old, dtdt, run_off_tot, &
   & fluxsens, fluxlat, emis, z0, epot_sol
   REAL, ALLOCATABLE, DIMENSION(:,:,:) :: albedo
!-
   INTEGER, ALLOCATABLE, DIMENSION(:) :: kindex
   REAL, ALLOCATABLE, DIMENSION(:,:) :: lon, lat
!-
   REAL :: old_tair
   INTEGER :: nbindex
   REAL,DIMENSION(:), ALLOCATABLE :: lev
   REAL :: julian, ss
   INTEGER :: yy, mm, dd
!-
   CHARACTER(LEN=10) :: coupling
!-
!-
   CHARACTER(LEN=80) :: filename, restname_in, restname_out
   CHARACTER(LEN=30) :: time_str, var_name
!-
   INTEGER :: i, j, it, istp, istp_old, rest_id, it_for
!-
   INTEGER :: split, split_start, nb_spread, for_offset
   INTEGER :: itau_dep, itau_fin, itau_skip, itau_len
   INTEGER :: it_force, it_forcep1
   INTEGER,DIMENSION(2) :: ml
!-
   LOGICAL :: last_CALL, first_CALL, debug, IMPLICIT
   LOGICAL :: no_inter, inter_lin, netrad_cons
!
!---------------------------------------------------------------------
!-
   CALL ioconf_calendar('noleap')
   CALL ioget_calendar(un_an, un_jour)
!-
!- set debug to have more information
!-
   !Config  Key  = DEBUG_INFO
   !Config  Desc = Flag for debug information
   !Config  Def  = n
   !Config  Help = This option allows to switch on the output of debug
   !Config         information without recompiling the code.
!-
   debug = .FALSE.
   CALL getin('DEBUG_INFO',debug)
!=====================================================================
!- 1.0 This section defines the general set-up of the experiment :
!-   - Forcing data to be used with its time step
!-   - Restart file to be used
!-   - The time step that will be used
!-   - The starting date of the simulation
!-   - Length of integration
!-   - Basic flags for SSIPSL
!=====================================================================
!- 1.1 Initialize the driving variables. It essentialy gets the mode
!-     used and the size of the driving variables.
!=====================================================================
!- 
!-
!- 
   IF (debug) WRITE(*,*) 'Reading name of the forcing file'
   !Config  Key  = FORCING_FILE
   !Config  Desc = Name of file containing the forcing data
   !Config  Def  = islscp_for.nc
   !Config  Help = This is the name of the file which should be opened
   !Config         for reading the forcing data of the dim0 model.
   !Config         The format of the file has to be netCDF and COADS
   !Config         compliant.
!- 
   filename='islscp_for.nc'
   CALL getin('FORCING_FILE',filename)
!-
   IF (debug) WRITE(*,*) 'Opening forcing file'
!-
!- We call flininfo to obtain the dimensions
!- of iim, jjm and others.
!- This information will allow us to allocate all the space needed.
!-
!-
   CALL forcing_info(filename, iim, jjm, llm, tm, date0, dt_force, &
        & force_id) 
!-
   IF (debug) WRITE(*,*) 'Allocate memeory for the driver :', iim, jjm, llm
!-
   ALLOCATE( lev(llm) )
   ALLOCATE &
   & (swdown(iim,jjm), precip_rain(iim,jjm), precip_snow(iim,jjm), tair_obs(iim,jjm), &
   & u(iim,jjm), v(iim,jjm), qair_obs(iim,jjm), pb(iim,jjm), lwdown(iim,jjm), &
   & eair_obs(iim,jjm), zlev_vec(iim,jjm), zlflu_vec(iim,jjm))
!- 
   ALLOCATE &
   & (petAcoef(iim,jjm), peqAcoef(iim,jjm), petBcoef(iim,jjm), peqBcoef(iim,jjm),&
   & for_cdrag(iim,jjm), for_u(iim,jjm), for_v(iim,jjm), for_swnet(iim,jjm),&
   & for_swdown(iim,jjm), for_lwdown(iim,jjm), &
   & for_psurf(iim,jjm), for_qair(iim,jjm), for_tair(iim,jjm), for_eair(iim,jjm), &
   & for_ccanopy(iim,jjm), for_rau(iim,jjm), tmp_eair(iim,jjm),&
   & tmp_tair(iim,jjm), tmp_qair(iim,jjm), tmp_pb(iim,jjm))
!-
   ALLOCATE &
   & (old_qair(iim,jjm), old_eair(iim,jjm), tsol_rad(iim,jjm), vevapp(iim,jjm), &
   & temp_sol_NEW(iim,jjm), temp_sol_old(iim,jjm), dtdt(iim,jjm), run_off_tot(iim,jjm), &
   & fluxsens(iim,jjm), fluxlat(iim,jjm), emis(iim,jjm), z0(iim,jjm), epot_sol(iim,jjm))

   ALLOCATE(albedo(iim,jjm,2))
!-
   ALLOCATE(kindex(iim*jjm))
   ALLOCATE(lon(iim,jjm), lat(iim,jjm))
!-
!- We need to know the grid. Then we can initialize the restart files, and then we
!- can give read the restart files in the forcing subroutine.
!-
   CALL forcing_grid (iim,jjm,llm,lon,lat,lev)
!-
!=====================================================================
!- 1.2 Time step to be used.
!-     This is relative to the time step of the  forcing data
!=====================================================================
   !Config  Key  = SPLIT_DT
   !Config  Desc = splits the timestep imposed by the forcing
   !Config  Def  = 12
   !Config  Help = With this value the time step of the forcing
   !Config         will be devided. In principle this can be run
   !Config         in explicit mode but it is strongly suggested
   !Config         to use the implicit method so that the
   !Config         atmospheric forcing has a smooth evolution.
!-
   split = 12
   CALL getin('SPLIT_DT', split)
!-
!- The time step which is going to be used is computed
!- from the forcing time step and the user input.
!-
   IF (split > 0) THEN
     dt = dt_force/split
   ELSE
     dt = dt_force
   ENDIF
!-
!=====================================================================
!- 1.3 Initialize the restart file for the driver
!=====================================================================
!-
   !
   !Config  Key  = RESTART_FILEIN
   !Config  Desc = Name of restart to READ for initial conditions
   !Config  Def  = NONE
   !Config  Help = This is the name of the file which will be opened
   !Config         to extract the initial values of all prognostic
   !Config         values of the model. This has to be a netCDF file.
   !Config         Not truly COADS compliant. NONE will mean that
   !Config         no restart file is to be expected.
!-
   restname_in = 'NONE'
   CALL getin('RESTART_FILEIN',restname_in)
   if ( debug ) WRITE(*,*) 'INPUT RESTART_FILE: ', restname_in(1:len_trim(restname_in))
!-
   !Config Key  = RESTART_FILEOUT
   !Config Desc = Name of restart files to be created by the driver
   !Config Def  = driver_rest_out.nc
   !Config Help = This variable give the  name for
   !Config        the restart files. The restart software within
   !Config        IOIPSL will add .nc if needed
!-
   restname_out = 'driver_rest_out.nc'
   CALL getin('RESTART_FILEOUT', restname_out)
   if ( debug ) WRITE(*,*) 'OUTPUT RESTART_FILE: ', restname_out(1:len_trim(restname_out))
!-
!- We choose some default values for the chronology
!- of the restart file.
!-
   dt_rest = dt
   date0_rest = date0
!-
!- Set default values for the start and end of the simulation
!- in the forcing chronology.
!-
   itau_dep = 0
   itau_fin = tm
!-
   CALL restini &
   & (restname_in, iim, jjm, lon, lat, llm, lev, &
   &  restname_out, itau_dep, date0_rest, dt_rest, rest_id)
!-
   WRITE(*,*) 'Restart file initialized : itau_dep, date0, dt_rest',&
        & itau_dep, date0_rest, dt_rest
!-
!=====================================================================
!- 1.4 Here we do the first real reading of the driving. It only
!-     gets a few variables.
!=====================================================================
!-
!- prepares kindex table from the information
!- obtained from the forcing data and reads some initial values for
!- temperature, etc.
!- 
   kindex(1) = 1
!-
!-
   CALL forcing_READ &
        & (filename, rest_id, .TRUE., .FALSE., &
        &  0, itau_dep, 0, split, nb_spread, netrad_cons, date0, dt_force,     &
        &  iim, jjm, lon, lat, tm, &
        &  swdown, precip_rain, precip_snow, tair_obs, &
        &  u, v, qair_obs, pb, lwdown, &
        &  kindex, nbindex, force_id)
!-
!-
   im = iim; jm = jjm; lm = llm;
   IF ( iim .GT. 1 .AND. jjm .GT. 1) THEN
      jtest = INT(kindex(INT(nbindex/2))/iim)+1
      itest = kindex(INT(nbindex/2)) - (jtest-1)*iim
   ELSE
      jtest = 1
      itest = 1
   ENDIF
!-
   IF ((im /= iim) .AND. (jm /= jjm) .AND. (lm /= llm))  THEN
     WRITE (*,*) ' dimensions are not good. Verfify FILE :'
     WRITE (*,*) ' filename = ',filename
     WRITE (*,*) ' im, jm, lm lus         = ', im, jm, lm
     WRITE (*,*) ' iim, jjm, llm demandes = ', iim, jjm, llm
     STOP 'dim2_driver'
   ENDIF
!-
!=====================================================================
!- 1.5  Configures the time-steps and other parameters 
!-      of the run.
!=====================================================================
!-
   istp_old = itau_dep
   istp = itau_dep+1
!-
!- If the time steping of the restart is different from the one
!- of the forcing we need to convert the itau_dep into the
!- chronology of the forcing. This ensures that the forcing will
!- start at the date of the restart file. Obviously the the itau_fin
!- needs to be shifted as well !
!-
   IF ( dt_rest /= dt_force .AND. itau_dep > 1) THEN
     itau_dep = (itau_dep* dt_rest )/dt_force
     itau_fin = itau_dep + tm
   ENDIF
!-
!- We need to select the right position of the splited time steps.
!-
   IF ( (itau_dep-1)*split /= istp) THEN
     split_start = istp-(itau_dep-1)*split
   ELSE
     split_start = 1
   ENDIF
!-
!- Same things if the starting dates are not the same.
!- Everything should look as if we had only one forcing file !
!-
   IF ( date0 /= date0_rest) THEN
     WRITE(*,*) 'date0_rest , date0 : ',date0_rest , date0
     for_offset = (date0_rest-date0)*un_jour/dt_force
     date0 = date0_rest
   ELSE
     for_offset = 0
   ENDIF
   WRITE(*,*) 'OFFSET FOR THE data read :', for_offset
!-
   !Config  Key  = TIME_SKIP
   !Config  Desc = Time in the forcing file at which the model
   !Config         is started.
   !Config  Def  = 0
   !Config  Help = This time give the point in time at which the model
   !Config         should be started. The FORMAT of this date can be 
   !Config         either of the following :
   !Config  n   : time step n within the forcing file
   !Config  nS  : n seconds after the first time-step in the file
   !Config  nD  : n days after the first time-step
   !Config  nM  : n month after the first time-step (year of 365 days)
   !Config  nY  : n years after the first time-step (year of 365 days)
   !Config         Or combinations :
   !Config  nYmM: n years and m month
!-
   itau_skip = 0
   WRITE(time_str,'(I10)') itau_skip
   if ( debug ) CALL getin('TIME_SKIP', time_str)
!-
!- Transform into itau
!-
  CALL tlen2itau(time_str, dt_force, date0, itau_skip)
!-
  itau_dep = itau_dep+itau_skip
  it_force = MAX(MOD(itau_dep,INT(un_jour*un_an/dt_force)),1)+itau_skip
  WRITE(*,*) 'itau_dep, it_force :', itau_dep, it_force
!-
  itau_len = (itau_fin-1)-itau_dep
!-
   !Config  Key  = TIME_LENGTH
   !Config  Desc = Length of the integration in time.
   !Config  Def  = DEF
   !Config  Help = Length of integration. By default the entire length
   !Config         of the forcing is used. The FORMAT of this date can
   !Config         be either of the following :
   !Config  n   : time step n within the forcing file
   !Config  nS  : n seconds after the first time-step in the file
   !Config  nD  : n days after the first time-step
   !Config  nM  : n month after the first time-step (year of 365 days)
   !Config  nY  : n years after the first time-step (year of 365 days)
   !Config         Or combinations :
   !Config  nYmM: n years and m month
!-
   WRITE(time_str,'(I10)') itau_len
   CALL getin('TIME_LENGTH', time_str)
!-
!- Transform into itau
!-
   CALL tlen2itau(time_str, dt_force, date0, itau_len)
!-
   itau_fin = itau_dep+itau_len
!-
   WRITE(*,*) '>> Time origine in the forcing file :', date0
   WRITE(*,*) '>> Time origine in the restat file :', date0_rest
   WRITE(*,*) '>> Simulate starts at forcing time-step : ', itau_dep
   WRITE(*,*) '>> The splited time-steps start at (Sets the '
   WRITE(*,*) '>>     chronology for the history and restart files):', istp
   WRITE(*,*) '>> The time spliting starts at :', split_start
   WRITE(*,*) '>> Simulation ends at forcing time-step: ', itau_fin
   WRITE(*,*) '>> Length of the simulation is thus :', itau_len
   WRITE(*,*) '>> Length of the forcing data is in time-steps : ', tm

!=====================================================================
!- 2.0 This section is going to define the details by which
!-     the input data is going to be used to force the
!-     land-surface scheme. The tasks are the following :
!-   - Is the coupling going to be explicit or implicit
!-   - Type of interpolation to be used.
!-   - At which height are the atmospheric forcings going to be used ?
!-   - Is a relaxation method going to be used on the forcing
!-   - Does net radiation in the interpolated data need to be conserved
!-   - How do we distribute precipitation.
!=====================================================================
   !Config  Key  = IMPLICIT
   !Config  Desc = Coupling to the atmospheric variables is IMPLICIT
   !Config  Def  = n
   !Config  Help = Allows to configure the method used to couple
   !Config         the surface to the atmospheric variables.
   !Config         One can either impose the atmospheric temperature
   !Config         and moisture directly (IMPLICIT=n, i.e. explicit
   !Config         coupling) or relax towards it (IMPLICIT=y). The 
   !Config         first method in the standard one for PILPS-TYPE
   !Config         exercices but may lead to instabilities as it does 
   !Config         not take into account the fluxes computed by the SVAT.
   !Config         The second method is actually ONLY experimental
   !Config         but will hopefully lead to a way of forcing SVATS
   !Config         which takes into account the computed fluxes and
   !Config         produce better results. It may certainly help in
   !Config         cases WHERE the DATA is of bad quality (as was the
   !Config         CASE in PILPS 1c) or not frequent enough (more
   !Config         than 1 hour).
!-
   IMPLICIT = .FALSE.
   CALL getin('IMPLICIT', IMPLICIT)
!-
   IF (IMPLICIT) THEN
     coupling = 'implicit'
   ELSE
     coupling = 'explicit'
   ENDIF
!-
   !Config  Key  = RELAX_A
   !Config  Desc = Capacity of the relaxation layer
   !Config  If   = IMPLICIT
   !Config  Def  = 1.0
   !Config  Help = The capacity associated to the atmospheric
   !Config         conditions which are going to be computed
   !Config         in the implicit coupling. The largest stable
   !Config         value is 1.0 after that is looses all its
   !Config         meaning as you should USE the explicit coupling.
   !Config         Values below 1.0 reduce the weight of the
   !Config         relaxation and smoothes the atmospheric forcing.
   !Config         =>   0. < RELAX_A < 1.
!-
   alpha = 1.0
   CALL getin('RELAX_A', alpha)
!-
   relax = alpha
   IF ( IMPLICIT ) THEN
      zlflu = dt/relax
   ELSE
      zlflu = lev(1)/2.
   ENDIF
!-
   !Config  Key  = NO_INTER||INTER_LIN
   !Config  Desc = Interpolation or not IF split is larger than 1
   !Config  Def  = NO_INTER
   !Config  Help = Choose IF you wish to interpolate linearly or not.
!-
   no_inter = .TRUE.
   inter_lin = .FALSE.
   CALL getin('NO_INTER', no_inter)
   CALL getin('INTER_LIN', inter_lin)
!-
   IF (inter_lin) no_inter = .FALSE.
!-
   IF (split == 1) THEN
     no_inter = .TRUE.
     inter_lin = .FALSE.
   ENDIF
!-
   !Config  Key  = SPRED_PREC
   !Config  Desc = Spread the precipitation.
   !Config  Def  = 1
   !Config  Help = Spread the the precipitaiton over n steps of the
   !Config         splited time step. This ONLY applied if the time
   !Config         steps have been splited.
!-
   nb_spread = 1
   CALL getin('SPRED_PREC', nb_spread)
   IF (nb_spread > split) THEN
     WRITE(*,*) 'WARNING : nb_spread is too large it will be '
     WRITE(*,*) '          set to the value of split'
     nb_spread = split
   ELSE IF (split == 1) THEN
     nb_spread = 1
   ENDIF
!-
   !Config  Key  = NETRAD_CONS
   !Config  Desc = Conserve net radiation in the forcing
   !Config  Def  = y
   !Config  If  = INTER_LIN
   !Config  Help = When the interpolation is used the net radiation
   !Config         provided by the forcing is not conserved anymore.
   !Config         This should be avoided and thus this option should
   !Config         be TRUE (y).
   !Config         This option is not used for short-wave if the time-step of 
   !Config         the forcing in longer than an hour. It does not make sense to try
   !Config         and reconstruct a diurnal cycle and at the same time conserve the 
   !Config         incoming solar radiation.
!-
   netrad_cons = .TRUE.
   CALL getin('NETRAD_CONS', netrad_cons)

!=====================================================================
!- 3.0 Finaly we can prepare all the variables before launching
!-     the simulation !
!=====================================================================
!-
!- Initialize LOGICAL and the length of the integration
!-
   last_CALL = .FALSE.
   first_CALL = .TRUE.
!-
   temp_sol_NEW(:,:) = tp_00
   zlev_vec(:,:) = lev(1)
!-
!- Preparing for the implicit scheme.
!- This means loading the prognostic varibales from the restart file.
!-
   CALL ioget_expval(val_exp)
!-
   var_name= 'fluxsens'
   CALL restget (rest_id, var_name, iim, jjm, 1, istp_old, &
   &             .TRUE., fluxsens)
   IF (ALL(fluxsens(:,:) == val_exp)) THEN
     fluxsens(:,:) = 0.
   ENDIF
!-
   var_name= 'vevapp'
   CALL restget (rest_id, var_name, iim, jjm, 1, istp_old, &
   &             .TRUE., vevapp)
   IF (ALL(vevapp(:,:) == val_exp)) THEN
     vevapp(:,:) = 0.
   ENDIF
!-
   var_name= 'qair_old'
   CALL restget (rest_id, var_name, iim, jjm, 1, istp_old, &
   &             .TRUE., old_qair)
   IF (ALL(old_qair(:,:) == val_exp)) THEN
     old_qair(:,:) = qair_obs(:,:)
   ENDIF
!-
   var_name= 'eair_old'
   CALL restget (rest_id, var_name, iim, jjm, 1, istp_old, &
   &             .TRUE., old_eair)
   IF (ALL(old_eair(:,:) == val_exp)) THEN
     old_eair(:,:) = &
   &  cp_air*(tair_obs(:,:))+cte_grav*zlev_vec(:,:)
   ENDIF
!- old density is also needed because we do not yet have the right pb
   var_name= 'rau_old'
   CALL restget (rest_id, var_name, iim, jjm, 1, istp_old, &
   &             .TRUE., for_rau)
   IF (ALL(for_rau(:,:) == val_exp)) THEN
     for_rau(:,:) = pb(:,:)/(cte_molr*(tair_obs(:,:)))
   ENDIF
!- For this variable the restart is extracted by SECHIBA
   temp_sol_NEW(:,:) = tair_obs(:,:)
!-
!- And other variables which need initial variables
!-
   albedo(:,:,:) = 0.13
!=====================================================================
!- 4.0 START THE TIME LOOP
!=====================================================================
   DO it = itau_dep, itau_fin
!---
     it_forcep1 = it_force+1
     IF (it_forcep1 > tm) it_forcep1 = it_force
!---
     DO is=split_start, split
       julian = itau2date(istp, date0, dt)
!-----
       it_for = it+for_offset
       IF (it_for < 0) THEN
         WRITE(*,*) 'The day is not in the forcing file :', it_for, it, for_offset
         STOP 'dim2_driver'
       ENDIF
       IF (it_for > tm) THEN
          WRITE(*,*) 'ERROR : more time-steps than data'
          WRITE(*,*) 'it_force : ', it_force, ' tm : ', tm
          STOP 'dim2_driver'
       ENDIF 
!----- 
       IF ( (it == itau_fin).AND.(is == split) ) THEN
         last_CALL = .TRUE.
       ENDIF
!-------
       IF ( debug ) WRITE(*,*) 'Into forcing_read'
!-------
       CALL forcing_READ &
            &  (filename, rest_id, .FALSE., last_call, &
            &   it_for, istp, is, split, nb_spread, netrad_cons, date0, dt_force, &
            &   iim, jjm, lon, lat, tm, &
            &   swdown, precip_rain, precip_snow, tair_obs, &
            &   for_u, for_v, qair_obs, pb, for_lwdown, &
            &   kindex, nbindex, force_id)
!
!-----
!----- Prepare : tmp_qair, tmp_eair, tmp_tair, tmp_pb
!----- and     : for_u, for_v, for_lwdown, for_swnet, for_swdown
!----- All the work is done in forcing_read
!----- We do not need the no_inter options as we will
!----- allways interpolate
!----- 
       IF ( debug ) WRITE(*,*) 'Prepare the atmospheric forcing'
!-----
       eair_obs(:,:) = cp_air*tair_obs(:,:)+cte_grav*zlev_vec(:,:)
       precip_rain(:,:) = precip_rain(:,:)*dt_force
       precip_snow(:,:) = precip_snow(:,:)*dt_force
!----- 
       for_swnet(:,:)  = (1.-(albedo(:,:,1)+albedo(:,:,2))/2.)*swdown(:,:)
       for_swdown(:,:) = swdown(:,:)
!----- 
!----- Computing the buffer zone !
!----- 
       IF (INDEX(coupling,'implicit') > 0) THEN
         for_qair(:,:) = old_qair(:,:) &
         & +dt/zlflu*(relax*(qair_obs(:,:)-old_qair(:,:)) &
         &            +vevapp(:,:)/(dt*for_rau(:,:)))
         for_eair(:,:) = old_eair(:,:) &
         & +dt/zlflu*(relax*(eair_obs(:,:)-old_eair(:,:)) &
         &            +fluxsens(:,:)/for_rau(:,:))
!-------
         for_tair(:,:) = (for_eair(:,:)-cte_grav*zlev_vec(:,:))/cp_air
!-------
         epot_sol(:,:) = cp_air*temp_sol_NEW(:,:)
!-------
         for_rau(:,:) = pb(:,:)/(cte_molr*for_tair(:,:))
!-------
         peqAcoef(:,:) = 0.0
         peqBcoef(:,:) = relax*(qair_obs(:,:)-for_qair(:,:))
         petAcoef(:,:) = 0.0
         petBcoef(:,:) = relax*(eair_obs(:,:)-for_eair(:,:))
       ELSE
         for_qair(:,:) = qair_obs(:,:)
         for_eair(:,:) = eair_obs(:,:)
         for_tair(:,:) = tair_obs(:,:)
!-------
         epot_sol(:,:) =  cp_air*temp_sol_NEW(:,:)
!-------
         for_rau(:,:) = pb(:,:)/(cte_molr*for_tair(:,:))
!-------
         petAcoef(:,:) = 0.0
         peqAcoef(:,:) = 0.0
         petBcoef(:,:) = 0.0
         peqBcoef(:,:) = 0.0
       ENDIF
!-----
       for_ccanopy(:,:) = 350.
       for_cdrag(:,:)  = 0.0
!-----
!-----  SECHIBA expects surface pressure in hPa
!-----
       for_psurf(:,:)  = pb(:,:)/100.
       zlflu_vec(:,:) = zlflu
!-----
!----- Prepare the other variables WITH the special CASE
!----- of splited time steps
!-----
!----- PRINT input value for debug
!-----
       IF (debug) THEN
         WRITE (*,*) ' >>>>>> time step it_force = ', it_force
         WRITE (*,*) ' tair, qair, eair = ', &
         & for_tair(itest,jtest), for_qair(itest,jtest), for_eair(itest,jtest)
         WRITE (*,*) ' OBS : tair, qair, eair = ', &
         & tair_obs(itest,jtest), qair_obs(itest,jtest), eair_obs(itest,jtest)
         WRITE (*,*) ' u et v = ', for_u(itest,jtest), for_v(itest,jtest)
         WRITE (*,*) ' precip rain et snow = ', &
         & precip_rain(itest,jtest), precip_snow(itest,jtest)
         WRITE (*,*) ' lwdown et swnet = ', &
         & for_lwdown(itest,jtest), for_swnet(itest,jtest)
         WRITE (*,*) ' petAcoef et peqAcoef = ', &
         & petAcoef(itest,jtest),  peqAcoef(itest,jtest)
         WRITE (*,*) ' petBcoef et peqAcoef = ', &
         & petBcoef(itest,jtest), peqBcoef(itest,jtest)
         WRITE (*,*) ' zlev = ', zlev_vec(itest,jtest)
       ENDIF
!-----
       IF (first_CALL) THEN
!-------
!------- CALL sechiba to initialize fields
!------- and have some initial results: emis, albedo, z0
!-------
         CALL intersurf_main &
         & (istp_old, iim, jjm, nbindex, kindex, &
         & dt, first_CALL, .FALSE., coupling, lon, lat, date0, &
         ! first level conditions
         & zlev_vec, zlflu_vec, for_u, for_v, &
         & for_qair, for_tair, for_eair, for_ccanopy, &
         ! Variables for the implicit coupling
         & for_cdrag, petAcoef, peqAcoef, petBcoef, peqBcoef, &
         ! Rain, snow, radiation and surface pressure
         & precip_rain, precip_snow, &
         & for_lwdown, for_swnet, for_swdown, for_psurf, &
         ! Output : Fluxes
         & vevapp, fluxsens, fluxlat, run_off_tot,  &
         ! Surface temperatures and surface properties
         & tsol_rad, temp_sol_NEW, albedo, emis, z0)
!-------
         first_CALL = .FALSE.
!-------
         temp_sol_old(:,:) = temp_sol_NEW(:,:)
         for_swnet(:,:)  = ( 1.- (albedo(:,:,1)+albedo(:,:,2))/2. ) * swdown(:,:)
         for_swdown(:,:)  = swdown(:,:)
!-------
         IF (debug) THEN
           WRITE (*,*) ' OUT rest : z0, albedoVis, albedoNIR, emis = ', &
           & z0(itest,jtest),albedo(itest,jtest,1),albedo(itest,jtest,2),emis(itest,jtest)
           WRITE (*,*) ' OUT rest : run_off_tot = ', &
           & run_off_tot(itest,jtest)
           WRITE (*,*) ' OUT rest : tsol_rad, vevapp = ', &
           & tsol_rad(itest,jtest), vevapp(itest,jtest)
           WRITE (*,*) ' OUT rest : temp_sol_new =', &
           & temp_sol_NEW(itest,jtest)
         ENDIF
       ENDIF
!-----
!----- Calling SECHIBA and doing the number crunching.
!----- Note that for the first time step SECHIBA is called twice.
!-----
       CALL intersurf_main &
       & (istp, iim, jjm, nbindex, kindex, &
       & dt, first_CALL, last_CALL, coupling, lon, lat, date0, &
       ! first level conditions
       &  zlev_vec, zlflu_vec, for_u, for_v, &
       &  for_qair, for_tair, for_eair, for_ccanopy, &
       ! Variables for the implicit coupling
       &  for_cdrag, petAcoef, peqAcoef, petBcoef, peqBcoef, &
       ! Rain, snow, radiation and surface pressure
       &  precip_rain, precip_snow, &
       &  for_lwdown, for_swnet, for_swdown, for_psurf, &
       ! Output : Fluxes
       &  vevapp, fluxsens, fluxlat, run_off_tot,  &
       ! Surface temperatures and surface properties
       &  tsol_rad, temp_sol_NEW, albedo, emis, z0)
!-----
       dtdt(:,:) = ABS(temp_sol_NEW(:,:)-temp_sol_old(:,:))/dt
!-----
!----- Test if the point with the largest change has more than 5K per dt
!-----
       IF (MAXVAL(dtdt(:,:)) > 5./dt) THEN
         ml = MAXLOC(dtdt)
         CALL ju2ymds(julian, yy, mm, dd, ss)
         WRITE(*,"('ATT :',I5,' big temperature jumps on ', &
         &          I4,'-',I2.2,'-',I2.2,':',F8.4)") &
         &         COUNT( dtdt(:,:) .GT. 5./dt ),yy, mm, dd, ss/3600.
         WRITE(*,*) &
         & 'Maximum change of surface temperature located at :', &
         & lon(ml(1),ml(2)), lat(ml(1),ml(2))
         WRITE(*,*) 'Coordinates in grid space: ',ml(1),ml(2)
         WRITE(*,*) 'Change from ', temp_sol_old(ml(1),ml(2)), &
         & ' to ', temp_sol_new(ml(1),ml(2)),&
         & 'with sw_in = ', for_swnet(ml(1),ml(2))
         old_tair = &
         & (old_eair(ml(1),ml(2))-cte_grav*zlev_vec(ml(1),ml(2)))/cp_air
         WRITE(*,*) 'Air temperature change from ', old_tair, &
         &          ' to ', for_tair(ml(1),ml(2))
         WRITE(*,*) 'Max of dtdt : ', dtdt(ml(1),ml(2)), &
         &          ' with dt = ', dt
       ENDIF
!-----
       temp_sol_old(:,:) = temp_sol_NEW(:,:)
!-----
!----- PRINT output value for debug
!-----
       IF (debug) THEN
         WRITE (*,*) ' OUT : z0, albedoVIS, albedoNIR, emis = ', &
         & z0(itest,jtest),albedo(itest,jtest,1),albedo(itest,jtest,2),emis(itest,jtest)
         WRITE (*,*) ' OUT : run_off_tot = ',  run_off_tot(itest,jtest)
         WRITE (*,*) ' OUT : tsol_rad, vevapp = ', &
         & tsol_rad(itest,jtest), vevapp(itest,jtest)
         WRITE (*,*) ' OUT : temp_sol_new =', temp_sol_NEW(itest,jtest)
       ENDIF
!-----
!----- Give some variables to the output package
!----- for saving on the history tape
!-----

!-----
       IF ( debug ) WRITE(*,*) 'history written for ', istp
!-----
       istp_old = istp
       istp = istp+1
!-----
       old_qair(:,:) = for_qair(:,:)
       old_eair(:,:) = for_eair(:,:)
     ENDDO
     split_start = 1
     it_force =it_force+1
     !!!CALL histsync
   ENDDO
!-
!- Archive in restart file the prognostic variables
!-
   IF ( debug ) WRITE(*,*) 'Write the restart for the driver', istp_old
!-
   var_name= 'fluxsens'
   CALL restput(rest_id, var_name, iim, jjm, 1, istp_old, fluxsens)
   var_name= 'vevapp'
   CALL restput(rest_id, var_name, iim, jjm, 1, istp_old, vevapp)
   var_name= 'qair_old'
   CALL restput(rest_id, var_name, iim, jjm, 1, istp_old, old_qair)
   var_name= 'eair_old'
   CALL restput(rest_id, var_name, iim, jjm, 1, istp_old, old_eair)
   var_name= 'rau_old'
   CALL restput(rest_id, var_name, iim, jjm, 1, istp_old, for_rau)
!-
   IF ( debug ) WRITE(*,*) 'Restart for the driver written'
!-
!-
!-
!=====================================================================
!- 5.0 Closing all files
!=====================================================================
  CALL flinclo(force_id)
  IF ( debug )  WRITE(*,*) 'FLIN CLOSED'
  CALL histclo
  IF ( debug )   WRITE(*,*) 'HIST CLOSED'
  CALL restclo
  IF ( debug )  WRITE(*,*) 'REST CLOSED'
!---------------
  STOP 'END of dim2_driver'
!---------------
END PROGRAM main
!-
