PROGRAM teststomate
!---------------------------------------------------------------------
! This program reads a forcing file for STOMATE
! which was created with STOMATE coupled to SECHIBA.
! It then integrates STOMATE for a
! given number of years using this forcing.
! The GPP is scaled to the new calculated LAI...
! Nothing is done for soil moisture which should
! actually evolve with the vegetation.
!---------------------------------------------------------------------
!- IPSL (2006)
!-  This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC
  USE netcdf
!-
  USE defprec
  USE constantes
  USE pft_parameters
  USE stomate_data
  USE ioipsl
  USE grid
  USE slowproc
  USE stomate
  USE intersurf, ONLY: stom_define_history, stom_ipcc_define_history, intsurf_time, l_first_intersurf, check_time
  USE parallel
!-
  IMPLICIT NONE
!-
! Declarations
!-
  INTEGER(i_std)                            :: kjpij,kjpindex
  REAL(r_std)                               :: dtradia,dt_force

  INTEGER(i_std),DIMENSION(:),ALLOCATABLE   :: indices
  INTEGER(i_std),DIMENSION(:),ALLOCATABLE   :: indexveg
  REAL(r_std),DIMENSION(:,:),ALLOCATABLE    :: soiltype, veget_x
  REAL(r_std),DIMENSION(:),ALLOCATABLE      :: totfrac_nobio
  REAL(r_std),DIMENSION(:,:),ALLOCATABLE    :: frac_nobio
  REAL(r_std),DIMENSION(:,:),ALLOCATABLE    :: veget_max_x
  REAL(r_std),DIMENSION(:,:),ALLOCATABLE    :: lai_x
  REAL(r_std),DIMENSION(:,:),ALLOCATABLE    :: veget_force_x
  REAL(r_std),DIMENSION(:,:),ALLOCATABLE    :: veget_max_force_x
  REAL(r_std),DIMENSION(:,:),ALLOCATABLE    :: lai_force_x
  REAL(r_std),DIMENSION(:),ALLOCATABLE      :: t2m,t2m_min,temp_sol
  REAL(r_std),DIMENSION(:,:),ALLOCATABLE    :: soiltemp,soilhum
  REAL(r_std),DIMENSION(:,:),ALLOCATABLE    :: humrel_x
  REAL(r_std),DIMENSION(:),ALLOCATABLE      :: litterhum
  REAL(r_std),DIMENSION(:),ALLOCATABLE      :: precip_rain,precip_snow
  REAL(r_std),DIMENSION(:,:),ALLOCATABLE    :: gpp_x
  REAL(r_std),DIMENSION(:),ALLOCATABLE      :: deadleaf_cover
  REAL(r_std),DIMENSION(:,:,:),ALLOCATABLE  :: assim_param_x
  REAL(r_std),DIMENSION(:,:),ALLOCATABLE    :: height_x
  REAL(r_std),DIMENSION(:,:),ALLOCATABLE    :: qsintmax_x
  REAL(r_std),DIMENSION(:,:),ALLOCATABLE    :: co2_flux
  REAL(r_std),DIMENSION(:),ALLOCATABLE      :: fco2_lu

  INTEGER(i_std),DIMENSION(:),ALLOCATABLE   :: indices_g
  REAL(r_std),DIMENSION(:),ALLOCATABLE   :: x_indices_g
  REAL(r_std),DIMENSION(:,:),ALLOCATABLE :: x_neighbours_g

  INTEGER    :: ier,iret
  INTEGER    :: ncfid
  CHARACTER(LEN=20),SAVE                      :: thecalendar='noleap'

  LOGICAL    :: a_er
  CHARACTER(LEN=80) :: &
 &  dri_restname_in,dri_restname_out, &
 &  sec_restname_in,sec_restname_out, &
 &  sto_restname_in,sto_restname_out, &
 &  stom_histname, stom_ipcc_histname
  INTEGER(i_std)                    :: iim,jjm,llm
  REAL, ALLOCATABLE, DIMENSION(:,:)  :: lon, lat
  REAL, ALLOCATABLE, DIMENSION(:)    :: lev
  LOGICAL                            :: rectilinear
  REAL, ALLOCATABLE, DIMENSION(:)    :: lon_rect, lat_rect
  REAL(r_std)                       :: dt
  INTEGER(i_std)                    :: itau_dep,itau,itau_len,itau_step
  REAL(r_std)                       :: date0
  INTEGER(i_std)                    :: rest_id_sec,rest_id_sto
  INTEGER(i_std)                    :: hist_id_sec,hist_id_sec2,hist_id_stom,hist_id_stom_IPCC
  CHARACTER(LEN=30)                :: time_str
  REAL(r_std)                       :: dt_slow_
  REAL                             :: hist_days_stom,hist_days_stom_ipcc,hist_dt_stom,hist_dt_stom_ipcc
!!$ >>  DS for the externalized version
  REAL(r_std),ALLOCATABLE,DIMENSION(:) :: hist_PFTaxis
!!$ >>  DS
  REAL(r_std),DIMENSION(10)         :: hist_pool_10axis     
  REAL(r_std),DIMENSION(100)        :: hist_pool_100axis     
  REAL(r_std),DIMENSION(11)         :: hist_pool_11axis     
  REAL(r_std),DIMENSION(101)        :: hist_pool_101axis     
  INTEGER                          :: hist_PFTaxis_id,hist_IPCC_PFTaxis_id,hori_id
  INTEGER                          :: hist_pool_10axis_id
  INTEGER                          :: hist_pool_100axis_id
  INTEGER                          :: hist_pool_11axis_id
  INTEGER                          :: hist_pool_101axis_id
  INTEGER(i_std)                    :: i,j,iv
  LOGICAL                          :: ldrestart_read,ldrestart_write
  LOGICAL                          :: l1d
  INTEGER(i_std),PARAMETER          :: nbvarmax=200
  INTEGER(i_std)                    :: nbvar
  CHARACTER*50,DIMENSION(nbvarmax) :: varnames
  INTEGER(i_std)                    :: varnbdim
  INTEGER(i_std),PARAMETER          :: varnbdim_max=20
  INTEGER,DIMENSION(varnbdim_max)  :: vardims
  CHARACTER*200                    :: taboo_vars
  REAL(r_std),DIMENSION(1)         :: xtmp
  INTEGER                          :: nsfm,nsft
  INTEGER                          :: iisf,iiisf
  INTEGER(i_std)                    :: max_totsize,totsize_1step,totsize_tmp

  INTEGER                          :: vid
  CHARACTER(LEN=100)               :: forcing_name
  REAL                             :: x

  REAL(r_std),DIMENSION(:,:),ALLOCATABLE :: var_3d

!-
  REAL(r_std)                             :: time_sec,time_step_sec
  REAL(r_std)                             :: time_dri,time_step_dri
  REAL(r_std),DIMENSION(1)                :: r1d
  REAL(r_std)                             :: julian,djulian

  INTEGER(i_std)                                :: ji,jv,l

  LOGICAL :: debug

!---------------------------------------------------------------------

  CALL init_para(.FALSE.)
  CALL init_timer

  IF (is_root_prc) THEN
     !-
     ! open STOMATE's forcing file to read some basic info
     !-
     forcing_name = 'stomate_forcing.nc'
     CALL getin ('STOMATE_FORCING_NAME',forcing_name)
     iret = NF90_OPEN (TRIM(forcing_name),NF90_NOWRITE,forcing_id)
     IF (iret /= NF90_NOERR) THEN
        CALL ipslerr (3,'teststomate', &
             &        'Could not open file : ', &
             &          forcing_name,'(Do you have forget it ?)')
     ENDIF
     ier = NF90_GET_ATT (forcing_id,NF90_GLOBAL,'dtradia',dtradia)
     ier = NF90_GET_ATT (forcing_id,NF90_GLOBAL,'dt_slow',dt_force)
     ier = NF90_GET_ATT (forcing_id,NF90_GLOBAL,'nsft',x)
     nsft = NINT(x)
     ier = NF90_GET_ATT (forcing_id,NF90_GLOBAL,'kjpij',x)
     kjpij = NINT(x)
     ier = NF90_GET_ATT (forcing_id,NF90_GLOBAL,'kjpindex',x)
     nbp_glo = NINT(x)
  ENDIF
  CALL bcast(dtradia)
  CALL bcast(dt_force)
  CALL bcast(nsft)
  CALL bcast(nbp_glo)
  !-
  write(numout,*) 'ATTENTION',dtradia,dt_force
  !-
  ! read info about land points
  !-
  IF (is_root_prc) THEN
     a_er=.FALSE.
     ALLOCATE (indices_g(nbp_glo),stat=ier)
     a_er = a_er .OR. (ier.NE.0)
     IF (a_er) THEN
        CALL ipslerr (3,'teststomate', &
             &        'PROBLEM WITH ALLOCATION', &
             &        'for local variables 1','')
     ENDIF
     !
     ALLOCATE (x_indices_g(nbp_glo),stat=ier)
     a_er = a_er .OR. (ier.NE.0)
     IF (a_er) THEN
        CALL ipslerr (3,'teststomate', &
             &        'PROBLEM WITH ALLOCATION', &
             &        'for global variables 1','')
     ENDIF
     ier = NF90_INQ_VARID (forcing_id,'index',vid)
     IF (ier .NE. 0) THEN
        CALL ipslerr (3,'teststomate', &
             &        'PROBLEM WITH ALLOCATION', &
             &        'for global variables 1','')
     ENDIF
     ier = NF90_GET_VAR   (forcing_id,vid,x_indices_g)
     IF (iret /= NF90_NOERR) THEN
        CALL ipslerr (3,'teststomate', &
             &        'PROBLEM WITH variable "index" in file ', &
             &        forcing_name,'(check this file)')
     ENDIF
     indices_g(:) = NINT(x_indices_g(:))
     DEALLOCATE (x_indices_g)
  ELSE
     ALLOCATE (indices_g(0))
  ENDIF
!---------------------------------------------------------------------
!-
! 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_p('DEBUG_INFO',debug)
  !
  !Config Key  = LONGPRINT
  !Config Desc = ORCHIDEE will print more messages
  !Config Def  = n
  !Config Help = This flag permits to print more debug messages in the run.
  !
  long_print = .FALSE.
  CALL getin_p('LONGPRINT',long_print)

!!$ >> DS added for the externalised version
!!$ we need to know the pft parameters values used for stomate

  ! 1. Read the number of PFTs
  !
  !Config Key  = NVM
  !Config Desc = number of PFTs  
  !Config  if  = ANYTIME
  !Config  Def  = 13
  !Config  Help = The number of vegetation types define by the user
  !Config  Units = NONE
  CALL getin_p('NVM',nvm)

  ! 2. Allocate and read the pft parameters stomate specific
  CALL pft_parameters_main
!!$ DS

  !-
  ! activate CO2, STOMATE, but not sechiba
  !-
  control%river_routing = .FALSE.
  control%hydrol_cwrr = .FALSE.
  control%ok_sechiba = .FALSE.
  !
  control%stomate_watchout = .TRUE.
  control%ok_co2 = .TRUE.
  control%ok_stomate = .TRUE.
  !-
  ! is DGVM activated?
  !-
  control%ok_dgvm = .FALSE.
  CALL getin_p('STOMATE_OK_DGVM',control%ok_dgvm)
  WRITE(numout,*) 'LPJ is activated: ',control%ok_dgvm

!!$ >> DS now we search the values for the scalar parameters in  some .def files

  ! 07/2011
  CALL activate_sub_models(control%ok_sechiba,control%river_routing,control%ok_stomate)
  CALL veget_config
  !
  CALL getin_stomate_pft_parameters
  CALL getin_co2_parameters
   CALL getin_stomate_parameters
  IF (control%ok_dgvm) THEN 
     CALL getin_dgvm_parameters   
  ENDIF
!!$ >> DS 

  !-
  ! restart files
  !-
  IF (is_root_prc) THEN
     ! Sechiba's restart files
     sec_restname_in = 'sechiba_start.nc'
     CALL getin('SECHIBA_restart_in',sec_restname_in)
     WRITE(numout,*) 'SECHIBA INPUT RESTART_FILE: ',TRIM(sec_restname_in)
     IF ( TRIM(sec_restname_in) .EQ. 'NONE' ) THEN
        STOP 'Need a restart file for Sechiba'
     ENDIF
     sec_restname_out = 'sechiba_rest_out.nc'
     CALL getin('SECHIBA_rest_out',sec_restname_out)
     WRITE(numout,*) 'SECHIBA OUTPUT RESTART_FILE: ',TRIM(sec_restname_out)
     ! Stomate's restart files
     sto_restname_in = 'stomate_start.nc'
     CALL getin('STOMATE_RESTART_FILEIN',sto_restname_in)
     WRITE(numout,*) 'STOMATE INPUT RESTART_FILE: ',TRIM(sto_restname_in)
     sto_restname_out = 'stomate_rest_out.nc'
     CALL getin('STOMATE_RESTART_FILEOUT',sto_restname_out)
     WRITE(numout,*) 'STOMATE OUTPUT RESTART_FILE: ',TRIM(sto_restname_out)

     !-
     ! We need to know iim, jjm.
     ! Get them from the restart files themselves.
     !-
     iret = NF90_OPEN (sec_restname_in,NF90_NOWRITE,ncfid)
     IF (iret /= NF90_NOERR) THEN
        CALL ipslerr (3,'teststomate', &
             &        'Could not open file : ', &
             &          sec_restname_in,'(Do you have forget it ?)')
     ENDIF
     iret = NF90_INQUIRE_DIMENSION (ncfid,1,len=iim_g)
     iret = NF90_INQUIRE_DIMENSION (ncfid,2,len=jjm_g)
     iret = NF90_INQ_VARID (ncfid, "time", iv)
     iret = NF90_GET_ATT (ncfid, iv, 'calendar',thecalendar)
     iret = NF90_CLOSE (ncfid)
     i=INDEX(thecalendar,ACHAR(0))
     IF ( i > 0 ) THEN
        thecalendar(i:20)=' '
     ENDIF
  ENDIF
  CALL bcast(iim_g)
  CALL bcast(jjm_g)
  CALL bcast(thecalendar)
  !-
  ! calendar
  !-
  CALL ioconf_calendar (thecalendar)
  CALL ioget_calendar  (one_year,one_day)
  !
  ! Parallelization :
  !
  CALL init_data_para(iim_g,jjm_g,nbp_glo,indices_g)
  kjpindex=nbp_loc
  jjm=jj_nb
  iim=iim_g
  kjpij=iim*jjm
  IF (debug) WRITE(numout,*) "Local grid : ",kjpindex,iim,jjm
  !-
  !-
  ! read info about grids
  !-
  !-
  llm=1
  ALLOCATE(lev(llm))
  IF (is_root_prc) THEN
     !-
     ier = NF90_INQ_VARID (forcing_id,'lalo',vid)
     ier = NF90_GET_VAR   (forcing_id,vid,lalo_g)
     !-
     ALLOCATE (x_neighbours_g(nbp_glo,8),stat=ier)
     ier = NF90_INQ_VARID (forcing_id,'neighbours',vid)
     ier = NF90_GET_VAR   (forcing_id,vid,x_neighbours_g)
     neighbours_g(:,:) = NINT(x_neighbours_g(:,:))
     DEALLOCATE (x_neighbours_g)
     !-
     ier = NF90_INQ_VARID (forcing_id,'resolution',vid)
     ier = NF90_GET_VAR   (forcing_id,vid,resolution_g)
     !-
     ier = NF90_INQ_VARID (forcing_id,'contfrac',vid)
     ier = NF90_GET_VAR   (forcing_id,vid,contfrac_g)

     lon_g(:,:) = 0.0
     lat_g(:,:) = 0.0
     lev(1)   = 0.0
     !-
     CALL restini &
          & (sec_restname_in, iim_g, jjm_g, lon_g, lat_g, llm, lev, &
          &  sec_restname_out, itau_dep, date0, dt, rest_id_sec)
     !-
     IF ( dt .NE. dtradia ) THEN
        WRITE(numout,*) 'dt',dt
        WRITE(numout,*) 'dtradia',dtradia
        CALL ipslerr (3,'teststomate', &
             &        'PROBLEM with time steps.', &
             &          sec_restname_in,'(dt .NE. dtradia)')
     ENDIF
     !-
     CALL restini &
          & (sto_restname_in, iim_g, jjm_g, lon_g, lat_g, llm, lev, &
          &  sto_restname_out, itau_dep, date0, dt, rest_id_sto)
     !-
     IF ( dt .NE. dtradia ) THEN
        WRITE(numout,*) 'dt',dt
        WRITE(numout,*) 'dtradia',dtradia
        CALL ipslerr (3,'teststomate', &
             &        'PROBLEM with time steps.', &
             &          sto_restname_in,'(dt .NE. dtradia)')
     ENDIF
  ENDIF
  CALL bcast(rest_id_sec)
  CALL bcast(rest_id_sto)
  CALL bcast(date0)
  CALL bcast(dt)
  CALL bcast(lev)
  !---
  !--- Create the index table
  !---
  !--- This job return a LOCAL kindex
  !---
  ALLOCATE (indices(kjpindex),stat=ier)
  IF (debug .AND. is_root_prc) WRITE(numout,*) "indices_g = ",indices_g(1:nbp_glo)
  CALL scatter(indices_g,indices)
  indices(1:kjpindex)=indices(1:kjpindex)-(jj_begin-1)*iim_g
  IF (debug) WRITE(numout,*) "indices = ",indices(1:kjpindex)

  !---
  !--- initialize global grid
  !---
  CALL init_grid( kjpindex ) 
  CALL grid_stuff (nbp_glo, iim_g, jjm_g, lon_g, lat_g, indices_g)

  !---
  !--- initialize local grid
  !---
  jlandindex = (((indices(1:kjpindex)-1)/iim) + 1)
  if (debug) WRITE(numout,*) "jlandindex = ",jlandindex(1:kjpindex)
  ilandindex = (indices(1:kjpindex) - (jlandindex(1:kjpindex)-1)*iim)
  IF (debug) WRITE(numout,*) "ilandindex = ",ilandindex(1:kjpindex)
  ALLOCATE(lon(iim,jjm))
  ALLOCATE(lat(iim,jjm))
  CALL scatter2D(lon_g,lon)
  CALL scatter2D(lat_g,lat)
  IF (debug) THEN
     IF (is_root_prc) THEN
        WRITE(numout,*) mpi_rank, lat_g
        WRITE(numout,*) mpi_rank, lon_g
     ENDIF
     WRITE(numout,*) mpi_rank, lat
     WRITE(numout,*) mpi_rank, lon
  ENDIF

  DO ji=1,kjpindex

     j = jlandindex(ji)
     i = ilandindex(ji)

     !- Create the internal coordinate table
!-
     lalo(ji,1) = lat(i,j)
     lalo(ji,2) = lon(i,j)
  ENDDO
  CALL scatter(neighbours_g,neighbours)
  CALL scatter(resolution_g,resolution)
  CALL scatter(contfrac_g,contfrac)
  CALL scatter(area_g,area)
  !-
  !- Check if we have by any change a rectilinear grid. This would allow us to 
  !- simplify the output files.
  !
  rectilinear = .FALSE.
  IF (is_root_prc) THEN
     IF ( ALL(lon_g(:,:) == SPREAD(lon_g(:,1), 2, SIZE(lon_g,2))) .AND. &
       & ALL(lat_g(:,:) == SPREAD(lat_g(1,:), 1, SIZE(lat_g,1))) ) THEN
        rectilinear = .TRUE.
     ENDIF
  ENDIF
  CALL bcast(rectilinear)
  IF (rectilinear) THEN
     ALLOCATE(lon_rect(iim),stat=ier)
     IF (ier .NE. 0) THEN
        WRITE (numout,*) ' error in lon_rect allocation. We stop. We need iim words = ',iim
        STOP 'intersurf_history'
     ENDIF
     ALLOCATE(lat_rect(jjm),stat=ier)
     IF (ier .NE. 0) THEN
        WRITE (numout,*) ' error in lat_rect allocation. We stop. We need jjm words = ',jjm
        STOP 'intersurf_history'
     ENDIF
     lon_rect(:) = lon(:,1)
     lat_rect(:) = lat(1,:)
  ENDIF
  !-
  ! allocate arrays
  !-
  !
  a_er = .FALSE.
!!$ >> DS added for the externalised version
  ALLOCATE (hist_PFTaxis(nvm),stat=ier)
  a_er = a_er .OR. (ier.NE.0)
! 01/07/2011 allocate fco2_lu
  ALLOCATE (fco2_lu(kjpindex),stat=ier)
  a_er = a_er .OR. (ier.NE.0)  
!!$ end add DS
  ALLOCATE (indexveg(kjpindex*nvm), stat=ier)
  a_er = a_er .OR. (ier.NE.0)
  ALLOCATE (soiltype(kjpindex,nstm),stat=ier)
  a_er = a_er .OR. (ier.NE.0)
  ALLOCATE (veget_x(kjpindex,nvm),stat=ier)
  a_er = a_er .OR. (ier.NE.0)
  ALLOCATE (totfrac_nobio(kjpindex),stat=ier)
  a_er = a_er .OR. (ier.NE.0)
  ALLOCATE (frac_nobio(kjpindex,nnobio),stat=ier)
  a_er = a_er .OR. (ier.NE.0)
  ALLOCATE (veget_max_x(kjpindex,nvm),stat=ier)
  a_er = a_er .OR. (ier.NE.0)
  ALLOCATE (lai_x(kjpindex,nvm),stat=ier)
  a_er = a_er .OR. (ier.NE.0)
  ALLOCATE (veget_force_x(kjpindex,nvm),stat=ier)
  a_er = a_er .OR. (ier.NE.0)
  ALLOCATE (veget_max_force_x(kjpindex,nvm),stat=ier)
  a_er = a_er .OR. (ier.NE.0)
  ALLOCATE (lai_force_x(kjpindex,nvm),stat=ier)
  a_er = a_er .OR. (ier.NE.0)
  ALLOCATE (t2m(kjpindex),stat=ier)
  a_er = a_er .OR. (ier.NE.0)
  ALLOCATE (t2m_min(kjpindex),stat=ier)
  a_er = a_er .OR. (ier.NE.0)
  ALLOCATE (temp_sol(kjpindex),stat=ier)
  a_er = a_er .OR. (ier.NE.0)
  ALLOCATE (soiltemp(kjpindex,nbdl),stat=ier)
  a_er = a_er .OR. (ier.NE.0)
  ALLOCATE (soilhum(kjpindex,nbdl),stat=ier)
  a_er = a_er .OR. (ier.NE.0)
  ALLOCATE (humrel_x(kjpindex,nvm),stat=ier)
  a_er = a_er .OR. (ier.NE.0)
  ALLOCATE (litterhum(kjpindex),stat=ier)
  a_er = a_er .OR. (ier.NE.0)
  ALLOCATE (precip_rain(kjpindex),stat=ier)
  a_er = a_er .OR. (ier.NE.0)
  ALLOCATE (precip_snow(kjpindex),stat=ier)
  a_er = a_er .OR. (ier.NE.0)
  ALLOCATE (gpp_x(kjpindex,nvm),stat=ier)
  a_er = a_er .OR. (ier.NE.0)
  ALLOCATE (deadleaf_cover(kjpindex),stat=ier)
  a_er = a_er .OR. (ier.NE.0)
  ALLOCATE (assim_param_x(kjpindex,nvm,npco2),stat=ier)
  a_er = a_er .OR. (ier.NE.0)
  ALLOCATE (height_x(kjpindex,nvm),stat=ier)
  a_er = a_er .OR. (ier.NE.0)
  ALLOCATE (qsintmax_x(kjpindex,nvm),stat=ier)
  a_er = a_er .OR. (ier.NE.0)
  ALLOCATE (co2_flux(kjpindex,nvm),stat=ier)
  a_er = a_er .OR. (ier.NE.0)
  ALLOCATE (fco2_lu(kjpindex),stat=ier)
  a_er = a_er .OR. (ier.NE.0)
  IF (a_er) THEN
     CALL ipslerr (3,'teststomate', &
          &        'PROBLEM WITH ALLOCATION', &
          &        'for local variables 1','')
  ENDIF
  !
  ! prepare forcing
  !
  max_totsize = 50
  CALL getin_p ('STOMATE_FORCING_MEMSIZE',max_totsize)
  max_totsize = max_totsize * 1000000

  totsize_1step = SIZE(soiltype(:,3))*KIND(soiltype(:,3)) + &
       SIZE(humrel_x)*KIND(humrel_x) + &
       SIZE(litterhum)*KIND(litterhum) + &
       SIZE(t2m)*KIND(t2m) + &
       SIZE(t2m_min)*KIND(t2m_min) + &
       SIZE(temp_sol)*KIND(temp_sol) + &
       SIZE(soiltemp)*KIND(soiltemp) + &
       SIZE(soilhum)*KIND(soilhum) + &
       SIZE(precip_rain)*KIND(precip_rain) + &
       SIZE(precip_snow)*KIND(precip_snow) + &
       SIZE(gpp_x)*KIND(gpp_x) + &
       SIZE(veget_force_x)*KIND(veget_force_x) + &
       SIZE(veget_max_force_x)*KIND(veget_max_force_x) + &
       SIZE(lai_force_x)*KIND(lai_force_x)
  CALL reduce_sum(totsize_1step,totsize_tmp)
  CALL bcast(totsize_tmp)
  totsize_1step=totsize_tmp 

! total number of forcing steps
  IF ( nsft .NE. INT(one_year/(dt_force/one_day)) ) THEN
     CALL ipslerr (3,'teststomate', &
          &        'stomate: error with total number of forcing steps', &
          &        'nsft','teststomate computation different with forcing file value.')
  ENDIF
! number of forcing steps in memory
  nsfm = MIN(nsft, &
       &       MAX(1,NINT( REAL(max_totsize,r_std) &
       &                  /REAL(totsize_1step,r_std))))
!-
  WRITE(numout,*) 'Offline forcing of Stomate:'
  WRITE(numout,*) '  Total number of forcing steps:',nsft
  WRITE(numout,*) '  Number of forcing steps in memory:',nsfm
!-
  CALL init_forcing(kjpindex,nsfm,nsft)
  !-
! ensure that we read all new forcing states
  iisf = nsfm
! initialize that table that contains the indices
! of the forcing states that will be in memory
  isf(:) = (/ (i,i=1,nsfm) /)

  nf_written(:) = .FALSE.
  nf_written(isf(:)) = .TRUE.

!-
! a time step for STOMATE corresponds to itau_step timesteps in SECHIBA
!-
  itau_step = NINT(dt_force/dtradia)
  IF (debug) WRITE(numout,*) "dtradia, dt_rest, dt_force, itau_step",dtradia, dt, dt_force, itau_step
  !
  CALL ioconf_startdate(date0)
  CALL intsurf_time( itau_dep, date0, dtradia )
!-
! Length of integration
!-
  WRITE(time_str,'(a)') '1Y'
  CALL getin ('TIME_LENGTH', time_str)
! transform into itau
  CALL tlen2itau(time_str, dt, date0, itau_len)
! itau_len*dtradia must be a multiple of dt_force
  itau_len = NINT( MAX(1.,FLOAT(NINT(itau_len*dtradia/dt_force))) &
       &                *dt_force/dtradia)
  !-
  ! set up STOMATE history file
  !-
  !Config  Key  = STOMATE_OUTPUT_FILE
  !Config  Desc = Name of file in which STOMATE's output is going
  !Config         to be written
  !Config  Def  = stomate_history.nc
  !Config  Help = This file is going to be created by the model
  !Config         and will contain the output from the model.
  !Config         This file is a truly COADS compliant netCDF file.
  !Config         It will be generated by the hist software from
  !Config         the IOIPSL package.
  !-
  stom_histname='stomate_history.nc'
  CALL getin_p ('STOMATE_OUTPUT_FILE', stom_histname)
  WRITE(numout,*) 'STOMATE_OUTPUT_FILE', TRIM(stom_histname)
  !-
  !Config  Key  = STOMATE_HIST_DT
  !Config  Desc = STOMATE history time step (d)
  !Config  Def  = 10.
  !Config  Help = Time step of the STOMATE history file
  !-
  hist_days_stom = 10.
  CALL getin_p ('STOMATE_HIST_DT', hist_days_stom)
  IF ( hist_days_stom == -1. ) THEN
     hist_dt_stom = -1.
     WRITE(numout,*) 'output frequency for STOMATE history file (d): one month.'
  ELSE
     hist_dt_stom = NINT( hist_days_stom ) * one_day
     WRITE(numout,*) 'output frequency for STOMATE history file (d): ', &
             hist_dt_stom/one_day
  ENDIF
!-
  !-
  !- initialize
  WRITE(numout,*) "before histbeg : ",date0,dt
  IF ( rectilinear ) THEN
#ifdef CPP_PARA
     CALL histbeg(stom_histname, iim, lon_rect, jjm, lat_rect,  1, iim, 1, jjm, &
          &     itau_dep, date0, dt, hori_id, hist_id_stom, domain_id=orch_domain_id)
#else
     CALL histbeg(stom_histname, iim, lon_rect, jjm, lat_rect,  1, iim, 1, jjm, &
          &     itau_dep, date0, dt, hori_id, hist_id_stom)
#endif
  ELSE
#ifdef CPP_PARA
     CALL histbeg(stom_histname, iim, lon, jjm, lat,  1, iim, 1, jjm, &
          &     itau_dep, date0, dt, hori_id, hist_id_stom, domain_id=orch_domain_id)
#else
     CALL histbeg(stom_histname, iim, lon, jjm, lat,  1, iim, 1, jjm, &
          &     itau_dep, date0, dt, hori_id, hist_id_stom)
#endif
  ENDIF
  !- define PFT axis
  hist_PFTaxis = (/ ( REAL(i,r_std), i=1,nvm ) /)
  !- declare this axis
  CALL histvert (hist_id_stom, 'PFT', 'Plant functional type', &
       & '1', nvm, hist_PFTaxis, hist_PFTaxis_id)
!- define Pool_10 axis
   hist_pool_10axis = (/ ( REAL(i,r_std), i=1,10 ) /)
!- declare this axis
  CALL histvert (hist_id_stom, 'P10', 'Pool 10 years', &
       & '1', 10, hist_pool_10axis, hist_pool_10axis_id)

!- define Pool_100 axis
   hist_pool_100axis = (/ ( REAL(i,r_std), i=1,100 ) /)
!- declare this axis
  CALL histvert (hist_id_stom, 'P100', 'Pool 100 years', &
       & '1', 100, hist_pool_100axis, hist_pool_100axis_id)

!- define Pool_11 axis
   hist_pool_11axis = (/ ( REAL(i,r_std), i=1,11 ) /)
!- declare this axis
  CALL histvert (hist_id_stom, 'P11', 'Pool 10 years + 1', &
       & '1', 11, hist_pool_11axis, hist_pool_11axis_id)

!- define Pool_101 axis
   hist_pool_101axis = (/ ( REAL(i,r_std), i=1,101 ) /)
!- declare this axis
  CALL histvert (hist_id_stom, 'P101', 'Pool 100 years + 1', &
       & '1', 101, hist_pool_101axis, hist_pool_101axis_id)

  !- define STOMATE history file
  CALL stom_define_history (hist_id_stom, nvm, iim, jjm, &
       & dt, hist_dt_stom, hori_id, hist_PFTaxis_id, &
 & hist_pool_10axis_id, hist_pool_100axis_id, &
 & hist_pool_11axis_id, hist_pool_101axis_id)

  !- end definition
  CALL histend(hist_id_stom)
  !-
  !-
  ! STOMATE IPCC OUTPUTS IS ACTIVATED
  !-
  !Config  Key  = STOMATE_IPCC_OUTPUT_FILE
  !Config  Desc = Name of file in which STOMATE's output is going
  !Config         to be written
  !Config  Def  = stomate_ipcc_history.nc
  !Config  Help = This file is going to be created by the model
  !Config         and will contain the output from the model.
  !Config         This file is a truly COADS compliant netCDF file.
  !Config         It will be generated by the hist software from
  !Config         the IOIPSL package.
  !-
  stom_ipcc_histname='stomate_ipcc_history.nc'
  CALL getin_p('STOMATE_IPCC_OUTPUT_FILE', stom_ipcc_histname)       
  WRITE(numout,*) 'STOMATE_IPCC_OUTPUT_FILE', TRIM(stom_ipcc_histname)
  !-
  !Config  Key  = STOMATE_IPCC_HIST_DT
  !Config  Desc = STOMATE IPCC history time step (d)
  !Config  Def  = 0.
  !Config  Help = Time step of the STOMATE IPCC history file
  !-
  hist_days_stom_ipcc = zero
  CALL getin_p('STOMATE_IPCC_HIST_DT', hist_days_stom_ipcc)       
  IF ( hist_days_stom_ipcc == moins_un ) THEN
     hist_dt_stom_ipcc = moins_un
     WRITE(numout,*) 'output frequency for STOMATE IPCC history file (d): one month.'
  ELSE
     hist_dt_stom_ipcc = NINT( hist_days_stom_ipcc ) * one_day
     WRITE(numout,*) 'output frequency for STOMATE IPCC history file (d): ', &
          hist_dt_stom_ipcc/one_day
  ENDIF

  ! test consistency between STOMATE_IPCC_HIST_DT and DT_SLOW parameters
  dt_slow_ = one_day
  CALL getin_p('DT_SLOW', dt_slow_)
  IF ( hist_days_stom_ipcc > zero ) THEN
     IF (dt_slow_ > hist_dt_stom_ipcc) THEN
        WRITE(numout,*) "DT_SLOW = ",dt_slow_,"  , STOMATE_IPCC_HIST_DT = ",hist_dt_stom_ipcc
        CALL ipslerr (3,'intsurf_history', &
             &          'Problem with DT_SLOW > STOMATE_IPCC_HIST_DT','', &
             &          '(must be less or equal)')
     ENDIF
  ENDIF

  IF ( hist_dt_stom_ipcc == 0 ) THEN
     hist_id_stom_ipcc = -1
  ELSE
     !-
     !- initialize
     IF ( rectilinear ) THEN
#ifdef CPP_PARA
        CALL histbeg(stom_ipcc_histname, iim, lon_rect, jjm, lat_rect,  1, iim, 1, jjm, &
             &     itau_dep, date0, dt, hori_id, hist_id_stom_ipcc,domain_id=orch_domain_id)
#else
        CALL histbeg(stom_ipcc_histname, iim, lon_rect, jjm, lat_rect,  1, iim, 1, jjm, &
             &     itau_dep, date0, dt, hori_id, hist_id_stom_ipcc)
#endif
     ELSE
#ifdef CPP_PARA
        CALL histbeg(stom_ipcc_histname, iim, lon, jjm, lat,  1, iim, 1, jjm, &
             &     itau_dep, date0, dt, hori_id, hist_id_stom_ipcc,domain_id=orch_domain_id)
#else
        CALL histbeg(stom_ipcc_histname, iim, lon, jjm, lat,  1, iim, 1, jjm, &
             &     itau_dep, date0, dt, hori_id, hist_id_stom_ipcc)
#endif
     ENDIF
     !- declare this axis
     CALL histvert (hist_id_stom_IPCC, 'PFT', 'Plant functional type', &
          & '1', nvm, hist_PFTaxis, hist_IPCC_PFTaxis_id)

     !- define STOMATE history file
     CALL stom_IPCC_define_history (hist_id_stom_IPCC, nvm, iim, jjm, &
          & dt, hist_dt_stom_ipcc, hori_id, hist_IPCC_PFTaxis_id)

     !- end definition
     CALL histend(hist_id_stom_IPCC)

  ENDIF
  !
  CALL histwrite(hist_id_stom, 'Areas',  itau_dep+itau_step, area, kjpindex, indices)
  IF ( hist_id_stom_IPCC > 0 ) THEN
     CALL histwrite(hist_id_stom_IPCC, 'Areas',  itau_dep+itau_step, area, kjpindex, indices)
  ENDIF
  !
  hist_id_sec = -1
  hist_id_sec2 = -1
!-
! first call of slowproc to initialize variables
!-
  itau = itau_dep
  !  
  DO ji=1,kjpindex
     DO jv=1,nvm
        indexveg((jv-1)*kjpindex + ji) = indices(ji) + (jv-1)*kjpij
     ENDDO
  ENDDO
  !-
  !MM Problem here with dpu which depends on soil type           
  DO l = 1, nbdl-1
     ! first 2.0 is dpu 
     ! second 2.0 is average
     diaglev(l) = dpu_cste/(2**(nbdl-1) -1) * ( ( 2**(l-1) -1) + ( 2**(l) -1) ) / 2.0
  ENDDO
  diaglev(nbdl) = dpu_cste
  !

  ! >> DS : getin the value of slowproc parameters
  !
  !Config Key  = CLAYFRACTION_DEFAULT
  !Config Desc = 
  !Config If   = OK_SECHIBA 
  !Config Def  = 0.2 
  !Config Help = 
  !Config Units = NONE   
  CALL getin_p('CLAYFRACTION_DEFAULT',clayfraction_default)
  !
  !Config Key  = MIN_VEGFRAC 
  !Config Desc = Minimal fraction of mesh a vegetation type can occupy 
  !Config If   = OK_SECHIBA 
  !Config Def  = 0.001 
  !Config Help = 
  !Config Units = NONE  
  CALL getin_p('MIN_VEGFRAC',min_vegfrac)
  !
  !Config Key  = STEMPDIAG_BID 
  !Config Desc = only needed for an initial LAI if there is no restart file
  !Config If   = OK_SECHIBA 
  !Config Def  = 280.
  !Config Help = 
  !Config Units = 
  CALL getin_p('STEMPDIAG_BID',stempdiag_bid)
  !
  !Config Key  = SOILTYPE_DEFAULT
  !Config Desc = Default soil texture distribution in the following order : sand, loam and clay
  !Config If   = OK_SECHIBA 
  !Config Def  = 0.0, 1.0, 0.0
  !Config Help = 
  !Config Units = NONE   
  CALL getin('SOILTYPE_DEFAULT',soiltype_default)
  ! >> DS

  CALL ioget_expval(val_exp)
  ldrestart_read = .TRUE.
  ldrestart_write = .FALSE.
  !-
  ! read some variables we need from SECHIBA's restart file
  !-
  CALL slowproc_main &
 &  (itau, kjpij, kjpindex, dt_force, date0, &
 &   ldrestart_read, ldrestart_write, .FALSE., .TRUE., &
 &   indices, indexveg, lalo, neighbours, resolution, contfrac, soiltype, &
 &   t2m, t2m_min, temp_sol, soiltemp, &
 &   humrel_x, soilhum, litterhum, precip_rain, precip_snow, gpp_x, &
 &   deadleaf_cover, assim_param_x, lai_x, height_x, veget_x, &
 &   frac_nobio, veget_max_x, totfrac_nobio, qsintmax_x, &
 &   rest_id_sec, hist_id_sec, hist_id_sec2, rest_id_sto, hist_id_stom, hist_id_stom_IPCC, co2_flux, fco2_lu)
  ! correct date
  day_counter = one_day - dt_force
  date=1
!-
! time loop
!-
  IF (debug) check_time=.TRUE.
  CALL intsurf_time( itau_dep+itau_step, date0, dtradia )
  l_first_intersurf=.FALSE.
  !
  DO itau = itau_dep+itau_step,itau_dep+itau_len,itau_step
     !
!-- next forcing state
    iisf = iisf+1
!---
    IF (iisf .GT. nsfm) THEN
!-----
!---- we have to read new forcing states from the file
!-----
!---- determine blocks of forcing states that are contiguous in memory
!-----
        CALL forcing_read(forcing_id,nsfm)

!--------------------------

!-----
!---- determine which forcing states must be read next time
!-----
      isf(1) = isf(nsfm)+1
      IF ( isf(1) .GT. nsft ) isf(1) = 1
        DO iiisf = 2, nsfm
           isf(iiisf) = isf(iiisf-1)+1
           IF ( isf(iiisf) .GT. nsft ) isf(iiisf) = 1
        ENDDO
        nf_written(isf(:)) = .TRUE.
!---- start again at first forcing state
        iisf = 1
     ENDIF
     ! Bug here ! soiltype(:,3) != clay
!     soiltype(:,3) = clay_fm(:,iisf)
     humrel_x(:,:) = humrel_daily_fm(:,:,iisf)
     litterhum(:) = litterhum_daily_fm(:,iisf)
     t2m(:) = t2m_daily_fm(:,iisf)
     t2m_min(:) = t2m_min_daily_fm(:,iisf)
     temp_sol(:) = tsurf_daily_fm(:,iisf)
     soiltemp(:,:) = tsoil_daily_fm(:,:,iisf)
     soilhum(:,:) = soilhum_daily_fm(:,:,iisf)
     precip_rain(:) = precip_fm(:,iisf)
     gpp_x(:,:) = gpp_daily_fm(:,:,iisf)
     veget_force_x(:,:) = veget_fm(:,:,iisf)
     veget_max_force_x(:,:) = veget_max_fm(:,:,iisf)
     lai_force_x(:,:) = lai_fm(:,:,iisf)
     WHERE ( t2m(:) .LT. ZeroCelsius )
        precip_snow(:) = precip_rain(:)
        precip_rain(:) = 0.
     ELSEWHERE
        precip_snow(:) = 0.
     ENDWHERE
!---
!-- scale GPP to new lai and veget_max
!---
     WHERE ( lai_x(:,:) .EQ. 0.0 ) gpp_x(:,:) = 0.0
!-- scale GPP to new LAI
     WHERE (lai_force_x(:,:) .GT. 0.0 )
        gpp_x(:,:) = gpp_x(:,:)*ATAN(2.*lai_x(:,:)) &
 &                           /ATAN( 2.*MAX(lai_force_x(:,:),0.01))
     ENDWHERE
!-- scale GPP to new veget_max
     WHERE (veget_max_force_x(:,:) .GT. 0.0 )
        gpp_x(:,:) = gpp_x(:,:)*veget_max_x(:,:)/veget_max_force_x(:,:)
     ENDWHERE
!---
!-- number crunching
!---
     ldrestart_read = .FALSE.
     ldrestart_write = .FALSE.
     CALL slowproc_main &
 &    (itau, kjpij, kjpindex, dt_force, date0, &
 &     ldrestart_read, ldrestart_write, .FALSE., .TRUE., &
 &     indices, indexveg, lalo, neighbours, resolution, contfrac, soiltype, &
 &     t2m, t2m_min, temp_sol, soiltemp, &
 &     humrel_x, soilhum, litterhum, precip_rain, precip_snow, gpp_x, &
 &     deadleaf_cover, assim_param_x, lai_x, height_x, veget_x, &
 &     frac_nobio, veget_max_x, totfrac_nobio, qsintmax_x, &
 &     rest_id_sec, hist_id_sec, hist_id_sec2, rest_id_sto, hist_id_stom, hist_id_stom_IPCC, co2_flux, fco2_lu)
     day_counter = one_day - dt_force
     !
     CALL intsurf_time( itau+itau_step, date0, dtradia )
  ENDDO ! end of the time loop
!-
  itau = itau -itau_step
!-
! write restart files
!-
  IF (is_root_prc) THEN
! first, read and write variables that are not managed otherwise
     taboo_vars = &
          &  '$lat$ $lon$ $lev$ $veget_year$ '// &
          &  '$height$ $day_counter$ $veget$ $veget_max$ $frac_nobio$ '// &
          &  '$lai$ $soiltype_frac$ $clay_frac$ '// &
          &  '$nav_lon$ $nav_lat$ $nav_lev$ $time$ $time_steps$'
!-
     CALL ioget_vname(rest_id_sec, nbvar, varnames)
!-
     DO iv = 1, nbvar
!-- check if the variable is to be written here
        IF (INDEX( taboo_vars,'$'//TRIM(varnames(iv))//'$') .EQ. 0 ) THEN
!---- get variable dimensions, especially 3rd dimension
           CALL ioget_vdim &
                &      (rest_id_sec, varnames(iv), varnbdim_max, varnbdim, vardims)
           l1d = ALL(vardims(1:varnbdim) .EQ. 1)
           ALLOCATE(var_3d(nbp_glo,vardims(3)),stat=ier)
           IF (ier .NE. 0) STOP 'ALLOCATION PROBLEM'
!---- read it
           IF (l1d) THEN
              CALL restget (rest_id_sec,TRIM(varnames(iv)), &
             1,vardims(3),1,itau_dep,.TRUE.,var_3d)
           ELSE
              CALL restget (rest_id_sec,TRIM(varnames(iv)), &
                   nbp_glo,vardims(3),1,itau_dep,.TRUE.,var_3d, &
                   "gather",nbp_glo,indices_g)
           ENDIF
!---- write it
           IF (l1d) THEN
              CALL restput (rest_id_sec,TRIM(varnames(iv)), &
                   1,vardims(3),1,itau,var_3d)
           ELSE
              CALL restput (rest_id_sec,TRIM(varnames(iv)), &
                   nbp_glo,vardims(3),1,itau,var_3d, &
                   'scatter',nbp_glo,indices_g)
           ENDIF
           DEALLOCATE (var_3d)
        ENDIF
     ENDDO
  ENDIF
  CALL barrier_para()

! call slowproc to write restart files
  ldrestart_read = .FALSE.
  ldrestart_write = .TRUE.
!-
  CALL slowproc_main &
 &  (itau, kjpij, kjpindex, dt_force, date0, &
 &   ldrestart_read, ldrestart_write, .FALSE., .TRUE., &
 &   indices, indexveg, lalo, neighbours, resolution, contfrac, soiltype, &
 &   t2m, t2m_min, temp_sol, soiltemp, &
 &   humrel_x, soilhum, litterhum, precip_rain, precip_snow, gpp_x, &
 &   deadleaf_cover, assim_param_x, lai_x, height_x, veget_x, &
 &   frac_nobio, veget_max_x, totfrac_nobio, qsintmax_x, &
 &   rest_id_sec, hist_id_sec, hist_id_sec2, rest_id_sto, hist_id_stom, hist_id_stom_IPCC, co2_flux, fco2_lu)
!-
! close files
!-
  IF (is_root_prc) &
       CALL restclo
  CALL histclo

  IF (is_root_prc) &
       ier = NF90_CLOSE (forcing_id)

  IF (is_root_prc) THEN
     CALL getin_dump()
     IF ( debug )  WRITE(numout,*) 'REST CLOSED'
  ENDIF
#ifdef CPP_PARA
  CALL MPI_FINALIZE(ier)
#endif
  WRITE(numout,*) "End of teststomate."

!---------------
END PROGRAM teststomate
!
!===
!