PROGRAM forcesoil
  !---------------------------------------------------------------------
  ! This program reads a forcing file for STOMATE's soil carbon routine
  ! which was created with STOMATE.
  ! It then integrates STOMATE's soil carbon parameterizations
  ! for a given number of years using this forcing.
  !---------------------------------------------------------------------
  !- IPSL (2006)
  !-  This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC
  USE netcdf
  !-
  USE defprec
  USE constantes
  USE constantes_mtc
  USE pft_parameters 
  USE stomate_data
  USE ioipsl
  USE stomate_soilcarbon
  USE parallel
  !-
  IMPLICIT NONE
  !-
  CHARACTER(LEN=80)          :: sto_restname_in,sto_restname_out
  INTEGER(i_std)             :: iim,jjm

  INTEGER(i_std),PARAMETER   :: llm = 1
  INTEGER(i_std)             :: kjpindex

  INTEGER(i_std)             :: itau_dep,itau_len
  CHARACTER(LEN=30)          :: time_str
  REAL(r_std)                :: dt_files
  REAL(r_std)                :: date0
  INTEGER(i_std)             :: rest_id_sto
  CHARACTER(LEN=20), SAVE    :: thecalendar = 'noleap'
  !-
  CHARACTER(LEN=100)         :: Cforcing_name
  INTEGER                    :: Cforcing_id
  INTEGER                    :: v_id
  REAL(r_std)                :: dt_forcesoil
  INTEGER                    :: nparan, nbyear

  INTEGER(i_std),DIMENSION(:),ALLOCATABLE    :: indices
  INTEGER(i_std),DIMENSION(:),ALLOCATABLE    :: indices_g
  REAL(r_std),DIMENSION(:),ALLOCATABLE       :: x_indices_g
  REAL(r_std),DIMENSION(:,:),ALLOCATABLE     :: lon, lat
  REAL(r_std),DIMENSION(llm)                 :: lev


  INTEGER                                    :: i,m,iatt,iv,iyear

  CHARACTER(LEN=80)                          :: var_name
  CHARACTER(LEN=800)                         :: taboo_vars
  REAL(r_std),DIMENSION(1)                   :: xtmp
  INTEGER(i_std),PARAMETER                   :: nbvarmax=300
  INTEGER(i_std)                             :: nbvar
  CHARACTER(LEN=50),DIMENSION(nbvarmax)      :: varnames
  INTEGER(i_std)                             :: varnbdim
  INTEGER(i_std),PARAMETER                   :: varnbdim_max=20
  INTEGER,DIMENSION(varnbdim_max)            :: vardims
  LOGICAL                                    :: l1d
  REAL(r_std),DIMENSION(:,:),ALLOCATABLE     :: var_3d
  REAL(r_std)                                :: x_tmp
  ! string suffix indicating an index
  CHARACTER(LEN=10)                          :: part_str
  !
  ! clay fraction
  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:)  :: clay_g
  REAL(r_std),DIMENSION(:,:,:,:),ALLOCATABLE :: soilcarbon_input_g
  REAL(r_std),DIMENSION(:,:,:),ALLOCATABLE   :: control_temp_g
  REAL(r_std),DIMENSION(:,:,:),ALLOCATABLE   :: control_moist_g
  REAL(r_std),DIMENSION(:,:,:),ALLOCATABLE   :: carbon_g

  REAL(r_std),ALLOCATABLE :: clay(:)
  REAL(r_std),ALLOCATABLE :: soilcarbon_input(:,:,:,:)
  REAL(r_std),ALLOCATABLE :: control_temp(:,:,:)
  REAL(r_std),ALLOCATABLE :: control_moist(:,:,:)
  REAL(r_std),ALLOCATABLE :: carbon(:,:,:)
  REAL(r_std),DIMENSION(:,:),ALLOCATABLE     :: resp_hetero_soil

  INTEGER(i_std)                             :: ier,iret

  CHARACTER(LEN=30) :: temp_name 
  LOGICAL :: debug
  LOGICAL :: l_error
 
  CALL Init_para(.FALSE.) 
  CALL init_timer

!-
! Configure the number of PFTS 
!-
  
  ! 1. Read the number of PFTs
  !
  !Config Key   = NVM
  !Config Desc  = number of PFTs  
  !Config If    = OK_SECHIBA or OK_STOMATE
  !Config Def   = 13
  !Config Help  = The number of vegetation types define by the user
  !Config Units = [-]
  CALL getin_p('NVM',nvm)

  ! 2. Allocation
  l_error = .FALSE.
  ALLOCATE(pft_to_mtc(nvm),stat=ier)
  l_error = l_error .OR. (ier .NE. 0)
  IF (l_error) THEN
     STOP 'pft_to_mtc (forcesoil only) : error in memory allocation'
  ENDIF

  ! 3. Initialisation of the correspondance table
  pft_to_mtc(:) = undef_int
  
  ! 4.Reading of the conrrespondance table in the .def file
  !
  !Config Key   = PFT_TO_MTC
  !Config Desc  = correspondance array linking a PFT to MTC
  !Config if    = OK_SECHIBA or OK_STOMATE
  !Config Def   = 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13
  !Config Help  =
  !Config Units = [-]
  CALL getin_p('PFT_TO_MTC',pft_to_mtc)

  ! 4.1 if nothing is found, we use the standard configuration
  IF(nvm .EQ. nvmc ) THEN
     IF(pft_to_mtc(1) .EQ. undef_int) THEN
        WRITE(numout,*) 'Note to the user : we will use ORCHIDEE to its standard configuration'
        pft_to_mtc(:) = (/ 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13 /)
     ENDIF
  ELSE   
     IF(pft_to_mtc(1) .EQ. undef_int) THEN
        WRITE(numout,*)' The array PFT_TO_MTC is empty : we stop'
     ENDIF
  ENDIF
  
  ! 4.2 What happened if pft_to_mtc(j) > nvmc (if the mtc doesn't exist)?
  DO i = 1, nvm
     IF(pft_to_mtc(i) .GT. nvmc) THEN
        WRITE(numout,*) "the MTC you chose doesn't exist"
        STOP 'we stop reading pft_to_mtc'
     ENDIF
  ENDDO
  
  ! 4.3 Check if pft_to_mtc(1) = 1 
  IF(pft_to_mtc(1) .NE. 1) THEN
     WRITE(numout,*) 'the first pft has to be the bare soil'
     STOP 'we stop reading next values of pft_to_mtc'
  ELSE
     DO i = 2,nvm
        IF(pft_to_mtc(i) .EQ.1) THEN
           WRITE(numout,*) 'only pft_to_mtc(1) has to be the bare soil'
           STOP 'we stop reading pft_to_mtc'
        ENDIF
     ENDDO
  ENDIF
  
  ! 5. Allocate and initialize natural ans is_c4
  
  ! 5.1 Memory allocation
  l_error = .FALSE.
  ALLOCATE(natural(nvm),stat=ier)
  l_error = l_error .OR. (ier .NE. 0)
  ALLOCATE(is_c4(nvm),stat=ier)

  IF (l_error) THEN
     STOP 'natural or is_c4 (forcesoil only) : error in memory allocation'
  ENDIF

  ! 5.2 Initialisation
  DO i = 1, nvm
     natural(i) = natural_mtc(pft_to_mtc(i))
     is_c4(i) = is_c4_mtc(pft_to_mtc(i))
  ENDDO

!-
!-
! set debug to have more information
!-
  !Config Key   = DEBUG_INFO
  !Config Desc  = Flag for debug information
  !Config If    = 
  !Config Def   = n
  !Config Help  = This option allows to switch on the output of debug
  !Config         information without recompiling the code.
  !Config Units = [FLAG]
!-
  debug = .FALSE.
  CALL getin_p('DEBUG_INFO',debug)
  !-
  ! Stomate's restart files
  !-
  IF (is_root_prc) THEN
     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_g, jjm.
     ! Get them from the restart files themselves.
     !-
     iret = NF90_OPEN (sto_restname_in, NF90_NOWRITE, rest_id_sto)
     iret = NF90_INQUIRE_DIMENSION (rest_id_sto,1,len=iim_g)
     iret = NF90_INQUIRE_DIMENSION (rest_id_sto,2,len=jjm_g)
     iret = NF90_INQ_VARID (rest_id_sto, "time", iv)
     iret = NF90_GET_ATT (rest_id_sto, iv, 'calendar',thecalendar)
     iret = NF90_CLOSE (rest_id_sto)
     i=INDEX(thecalendar,ACHAR(0))
     IF ( i > 0 ) THEN
        thecalendar(i:20)=' '
     ENDIF
     !-
     ! Allocate longitudes and latitudes
     !-
     ALLOCATE (lon(iim_g,jjm_g))
     ALLOCATE (lat(iim_g,jjm_g))
     lon(:,:) = zero
     lat(:,:) = zero
     lev(1)   = zero
     !-
     CALL restini &
          & (sto_restname_in, iim_g, jjm_g, lon, lat, llm, lev, &
          &  sto_restname_out, itau_dep, date0, dt_files, rest_id_sto)
  ENDIF
  CALL bcast(date0)
  CALL bcast(thecalendar)
  WRITE(numout,*) "calendar = ",thecalendar
  !-
  ! calendar
  !-
  CALL ioconf_calendar (thecalendar)
  CALL ioget_calendar  (one_year,one_day)
  CALL ioconf_startdate(date0)
  !
  IF (is_root_prc) THEN
     !-
     ! open FORCESOIL's forcing file to read some basic info
     !-
     Cforcing_name = 'NONE'
     CALL getin ('STOMATE_CFORCING_NAME',Cforcing_name)
     !-
     iret = NF90_OPEN (TRIM(Cforcing_name),NF90_NOWRITE,Cforcing_id)
     IF (iret /= NF90_NOERR) THEN
        CALL ipslerr (3,'forcesoil', &
             &        'Could not open file : ', &
             &          Cforcing_name,'(Do you have forget it ?)')
     ENDIF
     !-
     ier = NF90_GET_ATT (Cforcing_id,NF90_GLOBAL,'kjpindex',x_tmp)
     nbp_glo = NINT(x_tmp)
     ier = NF90_GET_ATT (Cforcing_id,NF90_GLOBAL,'nparan',x_tmp)
     nparan = NINT(x_tmp)
     ier = NF90_GET_ATT (Cforcing_id,NF90_GLOBAL,'nbyear',x_tmp)
     nbyear = NINT(x_tmp)
     !-
     ALLOCATE (indices_g(nbp_glo))
     ALLOCATE (clay_g(nbp_glo))
     !-
     ALLOCATE (x_indices_g(nbp_glo),stat=ier)
     ier = NF90_INQ_VARID (Cforcing_id,'index',v_id)
     ier = NF90_GET_VAR   (Cforcing_id,v_id,x_indices_g)
     indices_g(:) = NINT(x_indices_g(:))
     WRITE(numout,*) mpi_rank,"indices globaux : ",indices_g
     DEALLOCATE (x_indices_g)
     !-
     ier = NF90_INQ_VARID (Cforcing_id,'clay',v_id)
     ier = NF90_GET_VAR   (Cforcing_id,v_id,clay_g)
     !-
     ! time step of forcesoil
     !-
     dt_forcesoil = one_year / FLOAT(nparan)
     WRITE(numout,*) 'time step (d): ',dt_forcesoil
     WRITE(numout,*) 'nparan: ',nparan
     WRITE(numout,*) 'nbyear: ',nbyear    
     !-
     ! read (and partially write) the restart file
     !-
     ! read and write the variables we do not need
     !-
     taboo_vars ='$lon$ $lat$ $lev$ $nav_lon$ $nav_lat$ $nav_lev$ $time$ $time_steps$ '// &
          &             '$day_counter$ $dt_days$ $date$ '
     !-
     DO m = 1,nvm
        WRITE(part_str,'(I2)') m
        IF (m < 10) part_str(1:1) = '0'
        temp_name = '$carbon_'//part_str(1:LEN_TRIM(part_str))//'$'
        taboo_vars = TRIM(taboo_vars)//' '//TRIM(temp_name)
     ENDDO
     !-
     CALL ioget_vname(rest_id_sto, nbvar, varnames)
     !-
     ! read and write some special variables (1D or variables that we need)
     !-
     var_name = 'day_counter'
     CALL restget (rest_id_sto, var_name, 1, 1, 1, itau_dep, .TRUE., xtmp)
     CALL restput (rest_id_sto, var_name, 1, 1, 1, itau_dep, xtmp)
     !-
     var_name = 'dt_days'
     CALL restget (rest_id_sto, var_name, 1, 1, 1, itau_dep, .TRUE., xtmp)
     CALL restput (rest_id_sto, var_name, 1, 1, 1, itau_dep, xtmp)
     !-
     var_name = 'date'
     CALL restget (rest_id_sto, var_name, 1, 1, 1, itau_dep, .TRUE., xtmp)
     CALL restput (rest_id_sto, var_name, 1, 1, 1, itau_dep, xtmp)
     !-
     DO iv=1,nbvar
        !-- check if the variable is to be written here
        IF (INDEX(taboo_vars,'$'//TRIM(varnames(iv))//'$') == 0 ) THEN
           !---- get variable dimensions, especially 3rd dimension
           CALL ioget_vdim &
                &      (rest_id_sto, varnames(iv), varnbdim_max, varnbdim, vardims)
           l1d = ALL(vardims(1:varnbdim) == 1)
           !---- read it
           IF (l1d) THEN
              CALL restget &
                   &        (rest_id_sto, TRIM(varnames(iv)), 1, vardims(3), &
                   &         1, itau_dep, .TRUE., xtmp)
           ELSE
              ALLOCATE( var_3d(nbp_glo,vardims(3)), stat=ier)
              IF (ier /= 0) STOP 'ALLOCATION PROBLEM'
              !----
              CALL restget &
                   &        (rest_id_sto, 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_sto, TRIM(varnames(iv)), 1, vardims(3), &
                   &         1, itau_dep, xtmp)
           ELSE
              CALL restput &
                   &        (rest_id_sto, TRIM(varnames(iv)), nbp_glo, vardims(3), &
                   &         1, itau_dep, var_3d, 'scatter',  nbp_glo, indices_g)
              !----
              DEALLOCATE(var_3d)
           ENDIF
        ENDIF
     ENDDO
     !-
     ! read soil carbon
     !-
     ALLOCATE(carbon_g(nbp_glo,ncarb,nvm))
     carbon_g(:,:,:) = val_exp
     DO m = 1, nvm
        WRITE (part_str, '(I2)') m
        IF (m<10) part_str(1:1)='0'
        var_name = 'carbon_'//part_str(1:LEN_TRIM(part_str))
        CALL restget &
             &    (rest_id_sto, var_name, nbp_glo, ncarb , 1, itau_dep, &
             &     .TRUE., carbon_g(:,:,m), 'gather', nbp_glo, indices_g)
        IF (ALL(carbon_g(:,:,m) == val_exp)) carbon_g(:,:,m) = zero
        !-- do not write this variable: it will be modified.
     ENDDO
     WRITE(numout,*) "date0 : ",date0, itau_dep
     !-
     ! Length of run
     !-
     WRITE(time_str,'(a)') '10000Y'
     CALL getin('TIME_LENGTH', time_str)
     write(numout,*) 'Number of years for carbon spinup : ',time_str
     ! transform into itau
     CALL tlen2itau(time_str, dt_forcesoil*one_day, date0, itau_len)
     write(numout,*) 'Number of time steps to do: ',itau_len
     !-
     ! read the rest of the forcing file and store forcing in an array.
     ! We read an average year.
     !-
     ALLOCATE(soilcarbon_input_g(nbp_glo,ncarb,nvm,nparan*nbyear))
     ALLOCATE(control_temp_g(nbp_glo,nlevs,nparan*nbyear))
     ALLOCATE(control_moist_g(nbp_glo,nlevs,nparan*nbyear))
     !-
     ier = NF90_INQ_VARID (Cforcing_id,'soilcarbon_input',v_id)
     ier = NF90_GET_VAR   (Cforcing_id,v_id,soilcarbon_input_g)
     ier = NF90_INQ_VARID (Cforcing_id,   'control_moist',v_id)
     ier = NF90_GET_VAR   (Cforcing_id,v_id,control_moist_g)
     ier = NF90_INQ_VARID (Cforcing_id,    'control_temp',v_id)
     ier = NF90_GET_VAR   (Cforcing_id,v_id,control_temp_g)
     !-
     ier = NF90_CLOSE (Cforcing_id)
     !-
  ENDIF
  CALL bcast(nparan)
  CALL bcast(nbyear)
  CALL bcast(dt_forcesoil)
  CALL bcast(iim_g)
  CALL bcast(jjm_g)
  CALL bcast(nbp_glo)
  CALL bcast(itau_dep)
  CALL bcast(itau_len)
  !
  ! We must initialize data_para :
     !
     !
  CALL init_data_para(iim_g,jjm_g,nbp_glo,indices_g)

  kjpindex=nbp_loc
  jjm=jj_nb
  iim=iim_g
  IF (debug) WRITE(numout,*) "Local grid : ",kjpindex,iim,jjm

  !---
  !--- Create the index table
  !---
  !--- This job return a LOCAL kindex
  !---
  ALLOCATE (indices(kjpindex),stat=ier)
  CALL scatter(indices_g,indices)
  indices(1:kjpindex)=indices(1:kjpindex)-(jj_begin-1)*iim_g
  IF (debug) WRITE(numout,*) mpi_rank,"indices locaux = ",indices(1:kjpindex)
  !-
  !-
  ! there we go: time loop
  !-
  ALLOCATE(clay(kjpindex))
  ALLOCATE(soilcarbon_input(kjpindex,ncarb,nvm,nparan*nbyear))
  ALLOCATE(control_temp(kjpindex,nlevs,nparan*nbyear))
  ALLOCATE(control_moist(kjpindex,nlevs,nparan*nbyear))
  ALLOCATE(carbon(kjpindex,ncarb,nvm))
  ALLOCATE(resp_hetero_soil(kjpindex,nvm))
  iatt = 0

  CALL Scatter(clay_g,clay)
  CALL Scatter(soilcarbon_input_g,soilcarbon_input)
  CALL Scatter(control_temp_g,control_temp)
  CALL Scatter(control_moist_g,control_moist)
  CALL Scatter(carbon_g,carbon)

!-
! Configuration of the parameters
!-

  !Config Key   = FRAC_CARB_AP
  !Config Desc  = frac carb coefficients from active pool: depends on clay content
  !Config if    = OK_STOMATE 
  !Config Def   = 0.004
  !Config Help  = fraction of the active pool going to the passive pool
  !Config Units = [-]
  CALL getin_p('FRAC_CARB_AP',frac_carb_ap)  
  !
  !Config Key   = FRAC_CARB_SA
  !Config Desc  = frac_carb_coefficients from slow pool
  !Config if    = OK_STOMATE 
  !Config Def   = 0.42
  !Config Help  = fraction of the slow pool going to the active pool
  !Config Units = [-] 
  CALL getin_p('FRAC_CARB_SA',frac_carb_sa)
  !
  !Config Key   = FRAC_CARB_SP
  !Config Desc  = frac_carb_coefficients from slow pool
  !Config if    = OK_STOMATE 
  !Config Def   = 0.03
  !Config Help  = fraction of the slow pool going to the passive pool
  !Config Units = [-] 
  CALL getin_p('FRAC_CARB_SP',frac_carb_sp)
  !
  !Config Key   = FRAC_CARB_PA
  !Config Desc  = frac_carb_coefficients from passive pool
  !Config if    = OK_STOMATE 
  !Config Def   = 0.45
  !Config Help  = fraction of the passive pool going to the passive pool
  !Config Units = [-] 
  CALL getin_p('FRAC_CARB_PA',frac_carb_pa)
  !
  !Config Key   = FRAC_CARB_PS
  !Config Desc  = frac_carb_coefficients from passive pool
  !Config if    = OK_STOMATE 
  !Config Def   = 0.0
  !Config Help  = fraction of the passive pool going to the passive pool
  !Config Units = [-]
  CALL getin_p('FRAC_CARB_PS',frac_carb_ps)
  !
  !Config Key   = ACTIVE_TO_PASS_CLAY_FRAC
  !Config Desc  = 
  !Config if    = OK_STOMATE 
  !Config Def   =  .68  
  !Config Help  =
  !Config Units = [-]
  CALL getin_p('ACTIVE_TO_PASS_CLAY_FRAC',active_to_pass_clay_frac)
  !
  !Config Key   = CARBON_TAU_IACTIVE
  !Config Desc  = residence times in carbon pools
  !Config if    = OK_STOMATE 
  !Config Def   = 0.149
  !Config Help  =
  !Config Units = [days] 
  CALL getin_p('CARBON_TAU_IACTIVE',carbon_tau_iactive)
  !
  !Config Key   = CARBON_TAU_ISLOW
  !Config Desc  = residence times in carbon pools
  !Config if    = OK_STOMATE 
  !Config Def   = 5.48
  !Config Help  =
  !Config Units = [days]
  CALL getin_p('CARBON_TAU_ISLOW',carbon_tau_islow)
  !
  !Config Key   = CARBON_TAU_IPASSIVE
  !Config Desc  = residence times in carbon pools
  !Config if    = OK_STOMATE 
  !Config Def   = 241.
  !Config Help  =
  !Config Units = [days]
  CALL getin_p('CARBON_TAU_IPASSIVE',carbon_tau_ipassive)
  !
  !Config Key   = FLUX_TOT_COEFF
  !Config Desc  =
  !Config if    = OK_STOMATE 
  !Config Def   = 1.2, 1.4,.75
  !Config Help  =
  !Config Units = [days]
  CALL getin_p('FLUX_TOT_COEFF',flux_tot_coeff)

  iyear=1
  DO i=1,itau_len
     iatt = iatt+1
     IF (iatt > nparan*nbyear) THEN
        IF (debug) WRITE(numout,*) iyear
        iatt = 1
        iyear=iyear+1
     ENDIF
     CALL soilcarbon &
          &    (kjpindex, dt_forcesoil, clay, &
          &     soilcarbon_input(:,:,:,iatt), &
          &     control_temp(:,:,iatt), control_moist(:,:,iatt), &
          &     carbon, resp_hetero_soil)
  ENDDO
  WRITE(numout,*) "End of soilcarbon LOOP."
  CALL Gather(carbon,carbon_g)
  !-
  ! write new carbon into restart file
  !-
  IF (is_root_prc) THEN
     DO m=1,nvm
        WRITE (part_str, '(I2)') m
        IF (m<10) part_str(1:1)='0'
        var_name = 'carbon_'//part_str(1:LEN_TRIM(part_str))
        CALL restput &
             &    (rest_id_sto, var_name, nbp_glo, ncarb , 1, itau_dep, &
             &     carbon_g(:,:,m), 'scatter', nbp_glo, indices_g)
     ENDDO
     !-
     CALL getin_dump
     CALL restclo
  ENDIF
#ifdef CPP_PARA
  CALL MPI_FINALIZE(ier)
#endif
  WRITE(numout,*) "End of forcesoil."
  !--------------------
END PROGRAM forcesoil