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 pft_parameters 
  USE stomate_data
  USE ioipsl
  USE stomate_soilcarbon
  USE parallel
  !-
  IMPLICIT NONE
  !-
  CHARACTER(LEN=80) :: sto_restname_in,sto_restname_out,var_name
  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
  INTEGER(i_std)                             :: ier,iret
  REAL(r_std)                                :: dt_files
  REAL(r_std)                                :: date0
  INTEGER(i_std)                             :: rest_id_sto
  INTEGER(i_std)                             :: ncfid
  REAL(r_std)                                :: dt_force,dt_forcesoil
  INTEGER                                   :: nparan
  INTEGER,PARAMETER                         :: nparanmax=36
  REAL(r_std)                                :: xbid1,xbid2
  INTEGER(i_std)                             :: ibid
  INTEGER(i_std),DIMENSION(:),ALLOCATABLE    :: indices
  REAL(r_std),DIMENSION(llm)                 :: lev
  REAL(r_std),DIMENSION(:,:,:,:),ALLOCATABLE :: soilcarbon_input
  REAL(r_std),DIMENSION(:,:,:),ALLOCATABLE   :: &
       &  carbon,control_moist,control_temp
  REAL(r_std),DIMENSION(:,:),ALLOCATABLE     :: &
       &  lon,lat,resp_hetero_soil,var_3d
  REAL(r_std),DIMENSION(:),ALLOCATABLE       :: &
       &  x_indices
  REAL(r_std)                                :: time
  INTEGER                                   :: i,j,m,iatt,iv
  CHARACTER(LEN=400)                        :: 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)                                :: x_tmp
  ! clay fraction
  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:)  :: clay
  !-
  ! string suffix indicating an index
  CHARACTER(LEN=10)  :: part_str
  !
  CHARACTER(LEN=100) :: Cforcing_name
  INTEGER            :: Cforcing_id
  INTEGER            :: v_id

  REAL(r_std),ALLOCATABLE :: clay_loc(:)
  REAL(r_std),ALLOCATABLE :: soilcarbon_input_loc(:,:,:,:)
  REAL(r_std),ALLOCATABLE :: control_temp_loc(:,:,:)
  REAL(r_std),ALLOCATABLE :: control_moist_loc(:,:,:)
  REAL(r_std),ALLOCATABLE :: carbon_loc(:,:,:)
  INTEGER :: ierr
  !>> DS add for externalization
  LOGICAL  :: l_error
  ! >> DS

  CALL Init_para(.FALSE.) 

  !
  ! DS : For externalization cause we decoupled forcesoil from ORCHIDEE
  !
  
  ! 1. Read the number of PFTs
  CALL getin('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_integer
  
  ! 4.Reading of the conrrespondance table in the .def file
  CALL getin('PFT_TO_MTC',pft_to_mtc)

  ! 4.1 if nothing is found, we use the standard configuration
  IF(nvm .EQ. 13 ) THEN
     IF(pft_to_mtc(1) .EQ. undef_integer) 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_integer) THEN
        WRITE(numout,*)' The array PFT_TO_MTC is empty : we stop'
     ENDIF
  ENDIF
  
  ! 5. 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
  
  
  ! 6. 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
  

  ! 7. Allocate and initialize natural ans is_c4
  
  ! 7.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

  ! 7.2 Initialisation
  DO j= 1, nvm
     natural(j) = natural_mtc(pft_to_mtc(j))
     is_c4(j) = is_c4_mtc(pft_to_mtc(j))
  ENDDO

  !-
  ! Stomate's restart files
  !-
  IF (is_root_prc) THEN
     sto_restname_in = 'stomate_start.nc'
     CALL getin ('STOMATE_RESTART_FILEIN',sto_restname_in)
     WRITE(*,*) 'STOMATE INPUT RESTART_FILE: ',TRIM(sto_restname_in)
     sto_restname_out = 'stomate_restart.nc'
     CALL getin ('STOMATE_RESTART_FILEOUT',sto_restname_out)
     WRITE(*,*) 'STOMATE OUTPUT RESTART_FILE: ',TRIM(sto_restname_out)
     !-
     ! We need to know iim, jjm.
     ! Get them from the restart files themselves.
     !-
     iret = NF90_OPEN (sto_restname_in, NF90_NOWRITE, ncfid)
     iret = NF90_INQUIRE_DIMENSION (ncfid,1,len=iim)
     iret = NF90_INQUIRE_DIMENSION (ncfid,2,len=jjm)
     iret = NF90_CLOSE (ncfid)
     !-
     ! Allocate longitudes and latitudes
     !-
     ALLOCATE (lon(iim,jjm))
     ALLOCATE (lat(iim,jjm))
     lon(:,:) = 0.0
     lat(:,:) = 0.0
     lev(1)   = 0.0
     !-
     CALL restini &
          & (sto_restname_in, iim, jjm, lon, lat, llm, lev, &
          &  sto_restname_out, itau_dep, date0, dt_files, rest_id_sto)
  ENDIF
  !-
  ! calendar
  !-
  CALL bcast(date0)
!!! MM : à revoir : choix du calendrier dans forcesoil ?? Il est dans le restart de stomate !
  !  CALL ioconf_calendar ('noleap')
  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 = 'stomate_Cforcing.nc'
     CALL getin ('STOMATE_CFORCING_NAME',Cforcing_name)
     !-
     ier = NF90_OPEN (TRIM(Cforcing_name),NF90_NOWRITE,Cforcing_id)
     !-
     ier = NF90_GET_ATT (Cforcing_id,NF90_GLOBAL,'kjpindex',x_tmp)
     kjpindex = NINT(x_tmp)
     ier = NF90_GET_ATT (Cforcing_id,NF90_GLOBAL,'nparan',x_tmp)
     nparan = NINT(x_tmp)
     !-
     ALLOCATE (indices(kjpindex))
     ALLOCATE (clay(kjpindex))
     !-
     ALLOCATE (x_indices(kjpindex),stat=ier)
     ier = NF90_INQ_VARID (Cforcing_id,'index',v_id)
     ier = NF90_GET_VAR   (Cforcing_id,v_id,x_indices)
     indices(:) = NINT(x_indices(:))
     DEALLOCATE (x_indices)
     !-
     ier = NF90_INQ_VARID (Cforcing_id,'clay',v_id)
     ier = NF90_GET_VAR   (Cforcing_id,v_id,clay)
     !-
     ! time step of forcesoil
     !-
     dt_forcesoil = one_year / FLOAT(nparan)
     WRITE(*,*) 'time step (d): ',dt_forcesoil
     !-
     ! 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$ '// &
          &             '$carbon_01$ $carbon_02$ $carbon_03$ $carbon_04$ $carbon_05$'// &
          &             '$carbon_06$ $carbon_07$ $carbon_08$ $carbon_09$ $carbon_10$'// &
          &             '$carbon_11$ $carbon_12$ $carbon_13$'
     !-
     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)
           !----
           ALLOCATE( var_3d(kjpindex,vardims(3)), stat=ier)
           IF (ier /= 0) STOP 'ALLOCATION PROBLEM'
           !---- read it
           IF (l1d) THEN
              CALL restget &
                   &        (rest_id_sto, TRIM(varnames(iv)), 1, vardims(3), &
                   &         1, itau_dep, .TRUE., xtmp)
           ELSE
              CALL restget &
                   &        (rest_id_sto, TRIM(varnames(iv)), kjpindex, vardims(3), &
                   &         1, itau_dep, .TRUE., var_3d, "gather", kjpindex, indices)
           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)), kjpindex, vardims(3), &
                   &         1, itau_dep, var_3d, 'scatter',  kjpindex, indices)
           ENDIF
           !----
           DEALLOCATE(var_3d)
        ENDIF
     ENDDO
     !-
     ! read soil carbon
     !-
     ALLOCATE(carbon(kjpindex,ncarb,nvm))
     carbon(:,:,:) = 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_p &
             &    (rest_id_sto, var_name, kjpindex, ncarb , 1, itau_dep, &
             &     .TRUE., carbon(:,:,m), 'gather', kjpindex, indices)
        IF (ALL(carbon(:,:,m) == val_exp)) carbon(:,:,m) = zero
        !-- do not write this variable: it will be modified.
     ENDDO
     !-
     ! Length of run
     !-
     WRITE(time_str,'(a)') '10000Y'
     CALL getin('TIME_LENGTH', time_str)
     ! transform into itau
     CALL tlen2itau(time_str, dt_forcesoil*one_year, date0, itau_len)
     write(*,*) '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(kjpindex,ncarb,nvm,nparan))
     ALLOCATE(control_temp(kjpindex,nlevs,nparan))
     ALLOCATE(control_moist(kjpindex,nlevs,nparan))
     !-
     ier = NF90_INQ_VARID (Cforcing_id,'soilcarbon_input',v_id)
     ier = NF90_GET_VAR   (Cforcing_id,v_id,soilcarbon_input)
     ier = NF90_INQ_VARID (Cforcing_id,   'control_moist',v_id)
     ier = NF90_GET_VAR   (Cforcing_id,v_id,control_moist)
     ier = NF90_INQ_VARID (Cforcing_id,    'control_temp',v_id)
     ier = NF90_GET_VAR   (Cforcing_id,v_id,control_temp)
     !-
     ier = NF90_CLOSE (Cforcing_id)
     !-
     !MM Problem here with dpu which depends on soil type           
     DO iv = 1, nbdl-1
        ! first 2.0 is dpu 
        ! second 2.0 is average
        diaglev(iv) = 2.0/(2**(nbdl-1) -1) * ( ( 2**(iv-1) -1) + ( 2**(iv) -1) ) / 2.0
     ENDDO
     diaglev(nbdl) = 2.0
     !-
     ! For sequential use only, we must initialize data_para :
     !
     !
  ENDIF

  CALL bcast(iim)
  CALL bcast(jjm)
  call bcast(kjpindex)
  CALL init_data_para(iim,jjm,kjpindex,indices)

  !-
  !-
  ! there we go: time loop
  !-
  CALL bcast(nparan)
  ALLOCATE(clay_loc(nbp_loc))
  ALLOCATE(soilcarbon_input_loc(nbp_loc,ncarb,nvm,nparan))
  ALLOCATE(control_temp_loc(nbp_loc,nlevs,nparan))
  ALLOCATE(control_moist_loc(nbp_loc,nlevs,nparan))
  ALLOCATE(carbon_loc(nbp_loc,ncarb,nvm))
  ALLOCATE(resp_hetero_soil(nbp_loc,nvm))
  iatt = 0

  CALL bcast(itau_len)
  CALL bcast(nparan)
  CALL bcast(dt_forcesoil)
  CALL Scatter(clay,clay_loc)
  CALL Scatter(soilcarbon_input,soilcarbon_input_loc)
  CALL Scatter(control_temp,control_temp_loc)
  CALL Scatter(control_moist,control_moist_loc)
  CALL Scatter(carbon,carbon_loc)

!!$ DS 16/06/2011 : calling the new_values of soilcarbon parameters before loop
     !
     CALL getin_p('FRAC_CARB_AA',frac_carb_aa)
     CALL getin_p('FRAC_CARB_AP',frac_carb_ap)   
     CALL getin_p('FRAC_CARB_SS',frac_carb_ss)
     CALL getin_p('FRAC_CARB_SA',frac_carb_sa)
     CALL getin_p('FRAC_CARB_SP',frac_carb_sp)
     CALL getin_p('FRAC_CARB_PP',frac_carb_pp)
     CALL getin_p('FRAC_CARB_PA',frac_carb_pa)
     CALL getin_p('FRAC_CARB_PS',frac_carb_ps)
     !
     CALL getin_p('ACTIVE_TO_PASS_CLAY_FRAC',active_to_pass_clay_frac)
     CALL getin_p('CARBON_TAU_IACTIVE',carbon_tau_iactive)
     CALL getin_p('CARBON_TAU_ISLOW',carbon_tau_islow)
     CALL getin_p('CARBON_TAU_IPASSIVE',carbon_tau_ipassive)
     CALL getin_p('FLUX_TOT_COEFF',flux_tot_coeff)

  DO i=1,itau_len
     iatt = iatt+1
     IF (iatt > nparan) iatt = 1
     CALL soilcarbon &
          &    (nbp_loc, dt_forcesoil, clay_loc, &
          &     soilcarbon_input_loc(:,:,:,iatt), &
          &     control_temp_loc(:,:,iatt), control_moist_loc(:,:,iatt), &
          &     carbon_loc, resp_hetero_soil)
  ENDDO

  CALL Gather(carbon_loc,carbon)
  !-
  ! 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_p &
             &    (rest_id_sto, var_name, kjpindex, ncarb , 1, itau_dep, &
             &     carbon(:,:,m), 'scatter', kjpindex, indices)
     ENDDO
     !-
     CALL getin_dump
     CALL restclo
  ENDIF
#ifdef CPP_PARA
  CALL MPI_FINALIZE(ierr)
#endif
  !--------------------
END PROGRAM forcesoil