source: trunk/LMDZ.COMMON/libf/evolution/pemetat0.F90 @ 3129

MODULE pemetat0_mod

implicit none

!=======================================================================
contains
!=======================================================================

SUBROUTINE pemetat0(filename,ngrid,nsoil_GCM,nsoil_PEM,nslope,timelen,timestep,TI_PEM,tsoil_PEM,ice_table, ice_table_thickness,   &
                    tsurf_ave_yr1,tsurf_ave_yr2,q_co2,q_h2o,ps_inst,tsoil_inst,tend_h2oglaciers,tend_co2glaciers,co2ice,waterice, &
                    global_ave_pressure,watersurf_ave,watersoil_ave, m_co2_regolith_phys,deltam_co2_regolith_phys,                &
                    m_h2o_regolith_phys,deltam_h2o_regolith_phys, water_reservoir)

use iostart_PEM,                only: open_startphy, close_startphy, get_field, get_var
use comsoil_h_PEM,              only: soil_pem, layer_PEM, mlayer_PEM, fluxgeo, inertiedat_PEM, water_reservoir_nom, depth_breccia, depth_bedrock, index_breccia, index_bedrock
use comsoil_h,                  only: volcapa,inertiedat
use adsorption_mod,             only: regolith_adsorption, adsorption_pem
use ice_table_mod,              only: computeice_table_equilibrium
use constants_marspem_mod,      only: alpha_clap_h2o, beta_clap_h2o, TI_breccia, TI_bedrock
use soil_thermalproperties_mod, only: update_soil_thermalproperties
use tracer_mod,                 only: mmol, igcm_h2o_vap ! tracer names and molar masses
use abort_pem_mod,              only: abort_pem
use soil_pem_compute_mod,       only: soil_pem_compute
use soil_pem_ini_mod,           only: soil_pem_ini

#ifndef CPP_STD
    use comcstfi_h,   only: r, mugaz
    use surfdat_h,    only: watercaptag
#else
    use comcstfi_mod, only: r, mugaz
#endif

implicit none

include "callkeys.h"

character(len=*),               intent(in) :: filename            ! name of the startfi_PEM.nc
integer,                        intent(in) :: ngrid               ! # of physical grid points
integer,                        intent(in) :: nsoil_GCM           ! # of vertical grid points in the GCM
integer,                        intent(in) :: nsoil_PEM           ! # of vertical grid points in the PEM
integer,                        intent(in) :: nslope              ! # of sub-grid slopes
integer,                        intent(in) :: timelen             ! # time samples
real,                           intent(in) :: timestep            ! time step [s]
real,                           intent(in) :: global_ave_pressure ! mean average pressure on the planet [Pa]
real, dimension(ngrid,nslope),  intent(in) :: tsurf_ave_yr1       ! surface temperature at the first year of GCM call [K]
real, dimension(ngrid,nslope),  intent(in) :: tsurf_ave_yr2       ! surface temperature at the second  year of GCM call [K]
real, dimension(ngrid,timelen), intent(in) :: q_co2               ! MMR tracer co2 [kg/kg]
real, dimension(ngrid,timelen), intent(in) :: q_h2o               ! MMR tracer h2o [kg/kg]
real, dimension(ngrid,timelen), intent(in) :: ps_inst             ! surface pressure [Pa]
real, dimension(ngrid,nslope),  intent(in) :: tend_h2oglaciers    ! tendencies on h2o glaciers
real, dimension(ngrid,nslope),  intent(in) :: tend_co2glaciers    ! tendencies on co2 glaciers
real, dimension(ngrid,nslope),  intent(in) :: co2ice              ! co2 ice amount [kg/m^2]
real, dimension(ngrid,nslope),  intent(in) :: waterice            ! water ice amount [kg/m^2]
real, dimension(ngrid,nslope),  intent(in) :: watersurf_ave       ! surface water ice density, yearly averaged  (kg/m^3)
! outputs
real, dimension(ngrid),                          intent(out) :: deltam_co2_regolith_phys ! mass of co2 that is exchanged due to adsorption desorption [kg/m^2]
real, dimension(ngrid),                          intent(out) :: deltam_h2o_regolith_phys ! mass of h2o that is exchanged due to adsorption desorption [kg/m^2]
real, dimension(ngrid,nsoil_PEM,nslope),         intent(inout) :: TI_PEM              ! soil (mid-layer) thermal inertia in the PEM grid [SI]
real, dimension(ngrid,nsoil_PEM,nslope),         intent(inout) :: tsoil_PEM           ! soil (mid-layer) temperature [K]
real, dimension(ngrid,nslope),                   intent(inout) :: ice_table           ! Ice table depth [m]
real, dimension(ngrid,nslope),                   intent(inout) :: ice_table_thickness ! Ice table thickness [m]
real, dimension(ngrid,nsoil_PEM,nslope,timelen), intent(inout) :: tsoil_inst          ! instantaneous soil (mid-layer) temperature [K]
real, dimension(ngrid,nsoil_PEM,nslope),         intent(inout) :: m_co2_regolith_phys ! mass of co2 adsorbed [kg/m^2]
real, dimension(ngrid,nsoil_PEM,nslope),         intent(inout) :: m_h2o_regolith_phys ! mass of h2o adsorbed [kg/m^2]
real, dimension(ngrid),                          intent(inout) :: water_reservoir     ! mass of h2o that is exchanged due to adsorption desorption [kg/m^2]
real, dimension(ngrid,nsoil_PEM,nslope),         intent(inout) :: watersoil_ave       ! surface water ice density, yearly averaged (kg/m^3)
! local
real, dimension(ngrid,nsoil_PEM,nslope) :: tsoil_startPEM               ! soil temperature saved in the start [K]
real, dimension(ngrid,nsoil_PEM,nslope) :: TI_startPEM                  ! soil thermal inertia saved in the start [SI]
logical                                 :: found                        ! check if variables are found in the start
logical                                 :: found2                       ! check if variables are found in the start
integer                                 :: iloop, ig, islope, it, isoil ! index for loops
real                                    :: kcond                        ! Thermal conductivity, intermediate variable [SI]
real                                    :: delta                        ! Depth of the interface regolith-breccia, breccia -bedrock [m]
character(2)                            :: num                          ! intermediate string to read PEM start sloped variables
real, dimension(ngrid,nsoil_PEM)        :: tsoil_saved                  ! saved soil temperature [K]
real, dimension(ngrid,nsoil_PEM,nslope) :: tsoil_tmp_yr1                ! intermediate soil temperature during yr1[K]
real, dimension(ngrid,nsoil_PEM,nslope) :: tsoil_tmp_yr2                ! intermediate soil temperature during yr 2 [K]
real, dimension(ngrid,nsoil_PEM - 1)    :: alph_tmp                     ! Intermediate for tsoil computation []
real, dimension(ngrid,nsoil_PEM - 1)    :: beta_tmp                     ! Intermediate for tsoil computatio []
logical                                 :: startpem_file                ! boolean to check if we read the startfile or not

#ifdef CPP_STD
   logical, dimension(ngrid) :: watercaptag
   watercaptag(:) = .false.
#endif

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!
!!! Purpose: read start_pem. Need a specific iostart_PEM
!!!
!!! Order: 0. Previous year of the PEM run
!!!        1. Thermal Inertia
!!!        2. Soil Temperature
!!!        3. Ice table
!!!        4. Mass of CO2 & H2O adsorbed
!!!        5. Water reservoir
!!!
!!! /!\ This order must be respected !
!!! Author: LL
!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!0. Check if the start_PEM exist.

inquire(file=filename,exist =  startpem_file)

write(*,*)'Is start PEM?',startpem_file

!1. Run
if (startpem_file) then
    ! open pem initial state file:
    call open_startphy(filename)

    if (soil_pem) then

!1. Thermal Inertia
! a. General case
    do islope = 1,nslope
        write(num,'(i2.2)') islope
        call get_field("TI_PEM_slope"//num,TI_startPEM(:,:,islope),found)
        if (.not. found) then
            write(*,*)'PEM settings: failed loading <TI_PEM_slope'//num//'>'
            write(*,*)'will reconstruct the values of TI_PEM'

            do ig = 1,ngrid
                if (TI_PEM(ig,index_breccia,islope) < TI_breccia) then
                    !!! transition
                    delta = depth_breccia
                    TI_PEM(ig,index_breccia+1,islope) = sqrt((layer_PEM(index_breccia + 1) - layer_PEM(index_breccia))/ &
                                                        (((delta - layer_PEM(index_breccia))/(TI_PEM(ig,index_breccia,islope)**2)) + &
                                                        ((layer_PEM(index_breccia + 1) - delta)/(TI_breccia**2))))
                    do iloop=index_breccia+2,index_bedrock
                        TI_PEM(ig,iloop,islope) = TI_breccia
                    enddo
                else ! we keep the high ti values
                    do iloop = index_breccia + 1,index_bedrock
                        TI_PEM(ig,iloop,islope) = TI_PEM(ig,index_breccia,islope)
                    enddo
                endif ! TI PEM and breccia comparison
                !! transition
                delta = depth_bedrock
                TI_PEM(ig,index_bedrock + 1,islope) = sqrt((layer_PEM(index_bedrock + 1) - layer_PEM(index_bedrock))/ &
                                                        (((delta - layer_PEM(index_bedrock))/(TI_PEM(ig,index_bedrock,islope)**2)) + &
                                                        ((layer_PEM(index_bedrock + 1) - delta)/(TI_bedrock**2))))
                do iloop = index_bedrock + 2,nsoil_PEM
                    TI_PEM(ig,iloop,islope) = TI_bedrock
                enddo
            enddo
        else ! found
            do iloop = nsoil_GCM + 1,nsoil_PEM
                TI_PEM(:,iloop,islope) = TI_startPEM(:,iloop,islope)  ! ! 1st layers can change because of the presence of ice at the surface, so we don't change it here.
            enddo
        endif ! not found
    enddo ! islope

    write(*,*) 'PEMETAT0: THERMAL INERTIA done'

! b. Special case for inertiedat, inertiedat_PEM
call get_field("inertiedat_PEM",inertiedat_PEM,found)
if(.not.found) then
 do iloop = 1,nsoil_GCM
   inertiedat_PEM(:,iloop) = inertiedat(:,iloop)
 enddo
!!! zone de transition
delta = depth_breccia
do ig = 1,ngrid
if(inertiedat_PEM(ig,index_breccia).lt.TI_breccia) then
inertiedat_PEM(ig,index_breccia+1) = sqrt((layer_PEM(index_breccia+1)-layer_PEM(index_breccia))/ &
            (((delta-layer_PEM(index_breccia))/(inertiedat(ig,index_breccia)**2))+ &
                      ((layer_PEM(index_breccia+1)-delta)/(TI_breccia**2))))

 do iloop = index_breccia+2,index_bedrock
       inertiedat_PEM(ig,iloop) = TI_breccia
  enddo

else
   do iloop=index_breccia+1,index_bedrock
      inertiedat_PEM(ig,iloop) = inertiedat_PEM(ig,nsoil_GCM)
   enddo
endif ! comparison ti breccia
enddo!ig

!!! zone de transition
delta = depth_bedrock
do ig = 1,ngrid
inertiedat_PEM(ig,index_bedrock+1) = sqrt((layer_PEM(index_bedrock+1)-layer_PEM(index_bedrock))/ &
            (((delta-layer_PEM(index_bedrock))/(inertiedat_PEM(ig,index_bedrock)**2))+ &
                      ((layer_PEM(index_bedrock+1)-delta)/(TI_bedrock**2))))
enddo

 do iloop = index_bedrock+2, nsoil_PEM
   do ig = 1,ngrid
      inertiedat_PEM(ig,iloop) = TI_bedrock
   enddo
 enddo
endif ! not found

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!2. Soil Temperature
do islope=1,nslope
  write(num,fmt='(i2.2)') islope
   call get_field("tsoil_PEM_slope"//num,tsoil_startPEM(:,:,islope),found)
  if(.not.found) then
      write(*,*)'PEM settings: failed loading <tsoil_PEM_slope'//num//'>'
      write(*,*)'will reconstruct the values of Tsoil'
!      do ig = 1,ngrid
!        kcond = (TI_PEM(ig,index_breccia+1,islope)*TI_PEM(ig,index_breccia+1,islope))/volcapa
!        tsoil_PEM(ig,index_breccia+1,islope) = tsoil_PEM(ig,index_breccia,islope) + fluxgeo/kcond*(mlayer_PEM(index_breccia)-mlayer_PEM(index_breccia-1))
!       do iloop=index_breccia+2,index_bedrock
!            kcond = (TI_PEM(ig,iloop,islope)*TI_PEM(ig,iloop,islope))/volcapa
!            tsoil_PEM(ig,iloop,islope) = tsoil_PEM(ig,index_breccia+1,islope) + fluxgeo/kcond*(mlayer_PEM(iloop-1)-mlayer_PEM(index_breccia))
!      enddo
!        kcond = (TI_PEM(ig,index_bedrock+1,islope)*TI_PEM(ig,index_bedrock+1,islope))/volcapa
!        tsoil_PEM(ig,index_bedrock+1,islope) = tsoil_PEM(ig,index_bedrock,islope) + fluxgeo/kcond*(mlayer_PEM(index_bedrock)-mlayer_PEM(index_bedrock-1))
!
!      do iloop=index_bedrock+2,nsoil_PEM
!            kcond = (TI_PEM(ig,iloop,islope)*TI_PEM(ig,iloop,islope))/volcapa
!            tsoil_PEM(ig,iloop,islope) = tsoil_PEM(ig,index_bedrock+1,islope) + fluxgeo/kcond*(mlayer_PEM(iloop-1)-mlayer_PEM(index_bedrock))
!      enddo
!      enddo


     call soil_pem_ini(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_ave_yr2(:,islope),tsoil_PEM(:,:,islope))
     call soil_pem_compute(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_ave_yr2(:,islope),tsoil_PEM(:,:,islope))


   else
! predictor corrector: restart from year 1 of the GCM and build the evolution of
! tsoil at depth

    tsoil_tmp_yr1(:,:,islope) = tsoil_startPEM(:,:,islope)
    call soil_pem_compute(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_ave_yr1(:,islope),tsoil_tmp_yr1(:,:,islope))
    call soil_pem_compute(ngrid,nsoil_PEM,.false.,TI_PEM(:,:,islope),timestep,tsurf_ave_yr1(:,islope),tsoil_tmp_yr1(:,:,islope))
    tsoil_tmp_yr2(:,:,islope) = tsoil_tmp_yr1(:,:,islope)
    call soil_pem_compute(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_ave_yr2(:,islope),tsoil_tmp_yr2(:,:,islope))


     do iloop = nsoil_GCM+1,nsoil_PEM
       tsoil_PEM(:,iloop,islope) = tsoil_tmp_yr2(:,iloop,islope)
     enddo
   endif !found

    do it = 1,timelen
        do isoil = nsoil_GCM+1,nsoil_PEM
        tsoil_inst(:,isoil,islope,it) = tsoil_PEM(:,isoil,islope)
        enddo
     enddo
      do isoil = nsoil_GCM+1,nsoil_PEM
        do ig = 1,ngrid
        watersoil_ave(ig,isoil,islope) = exp(beta_clap_h2o/tsoil_PEM(ig,isoil,islope) + alpha_clap_h2o)/tsoil_PEM(ig,isoil,islope)*mmol(igcm_h2o_vap)/(mugaz*r)
        enddo
      enddo
enddo ! islope

write(*,*) 'PEMETAT0: SOIL TEMP done'

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!3. Ice Table
  call get_field("ice_table",ice_table,found)
   if(.not.found) then
      write(*,*)'PEM settings: failed loading <ice_table>'
      write(*,*)'will reconstruct the values of the ice table given the current state'
      ice_table(:,:) = -1  ! by default, no ice table  
      ice_table_thickness(:,:) = -1  ! by default, no ice table
     call computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,watersurf_ave,watersoil_ave, TI_PEM(:,1,:),ice_table,ice_table_thickness)
     call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,tend_h2oglaciers,waterice,global_ave_pressure,ice_table,ice_table_thickness,TI_PEM)
     do islope = 1,nslope
     call soil_pem_ini(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_ave_yr2(:,islope),tsoil_PEM(:,:,islope))
     enddo
   endif

write(*,*) 'PEMETAT0: ICE TABLE done'

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!4. CO2 & H2O Adsorption
 if(adsorption_pem) then
  do islope=1,nslope
   write(num,fmt='(i2.2)') islope
   call get_field("mco2_reg_ads_slope"//num,m_co2_regolith_phys(:,:,islope),found)
    if((.not.found)) then
       m_co2_regolith_phys(:,:,:) = 0.
       exit
    endif

  enddo

  do islope=1,nslope
   write(num,fmt='(i2.2)') islope
   call get_field("mh2o_reg_ads_slope"//num,m_h2o_regolith_phys(:,:,islope),found2)
    if((.not.found2)) then
       m_h2o_regolith_phys(:,:,:) = 0.
      exit
    endif

  enddo

    call regolith_adsorption(ngrid,nslope,nsoil_PEM,timelen,tend_h2oglaciers,tend_co2glaciers,waterice,co2ice,tsoil_PEM,TI_PEM,ps_inst,q_co2,q_h2o, &
                                m_h2o_regolith_phys,deltam_h2o_regolith_phys, m_co2_regolith_phys,deltam_co2_regolith_phys)

    if((.not.found)) then
      deltam_co2_regolith_phys(:) = 0.
    endif
    if((.not.found2)) then
       deltam_h2o_regolith_phys(:) = 0.
    endif
write(*,*) 'PEMETAT0: CO2 & H2O adsorption done'
    endif
endif ! soil_pem

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!. 5 water reservoir
#ifndef CPP_STD
   call get_field("water_reservoir",water_reservoir,found)
   if(.not.found) then
      write(*,*)'Pemetat0: failed loading <water_reservoir>'
      write(*,*)'will reconstruct the values from watercaptag'
      do ig=1,ngrid
        if(watercaptag(ig)) then
           water_reservoir(ig)=water_reservoir_nom
        else
           water_reservoir(ig)=0.
        endif
      enddo
   endif
#endif

call close_startphy

else !No startfi, let's build all by hand

    if (soil_pem) then

!a) Thermal inertia
   do islope = 1,nslope
      do ig = 1,ngrid

          if(TI_PEM(ig,index_breccia,islope).lt.TI_breccia) then
!!! transition
             delta = depth_breccia
             TI_PEM(ig,index_breccia+1,islope) =sqrt((layer_PEM(index_breccia+1)-layer_PEM(index_breccia))/ &
            (((delta-layer_PEM(index_breccia))/(TI_PEM(ig,index_breccia,islope)**2))+ &
                      ((layer_PEM(index_breccia+1)-delta)/(TI_breccia**2))))

          do iloop=index_breccia+2,index_bedrock
            TI_PEM(ig,iloop,islope) = TI_breccia
         enddo
         else ! we keep the high ti values
           do iloop=index_breccia+1,index_bedrock
                  TI_PEM(ig,iloop,islope) = TI_PEM(ig,index_breccia,islope)
           enddo
         endif

!! transition
             delta = depth_bedrock
             TI_PEM(ig,index_bedrock+1,islope) = sqrt((layer_PEM(index_bedrock+1)-layer_PEM(index_bedrock))/ &
            (((delta-layer_PEM(index_bedrock))/(TI_PEM(ig,index_bedrock,islope)**2))+ &
                      ((layer_PEM(index_bedrock+1)-delta)/(TI_breccia**2))))
          do iloop=index_bedrock+2,nsoil_PEM
            TI_PEM(ig,iloop,islope) = TI_bedrock
         enddo
      enddo
enddo

 do iloop = 1,nsoil_GCM
   inertiedat_PEM(:,iloop) = inertiedat(:,iloop)
 enddo
!!! zone de transition
delta = depth_breccia
do ig = 1,ngrid
      if(inertiedat_PEM(ig,index_breccia).lt.TI_breccia) then
inertiedat_PEM(ig,index_breccia+1) = sqrt((layer_PEM(index_breccia+1)-layer_PEM(index_breccia))/ &
            (((delta-layer_PEM(index_breccia))/(inertiedat(ig,index_breccia)**2))+ &
                      ((layer_PEM(index_breccia+1)-delta)/(TI_breccia**2))))


 do iloop = index_breccia+2,index_bedrock

       inertiedat_PEM(ig,iloop) = TI_breccia

 enddo
else
   do iloop = index_breccia+1,index_bedrock
       inertiedat_PEM(ig,iloop) = inertiedat_PEM(ig,index_breccia)
    enddo

endif
enddo

!!! zone de transition
delta = depth_bedrock
do ig = 1,ngrid
inertiedat_PEM(ig,index_bedrock+1) = sqrt((layer_PEM(index_bedrock+1)-layer_PEM(index_bedrock))/ &
            (((delta-layer_PEM(index_bedrock))/(inertiedat_PEM(ig,index_bedrock)**2))+ &
                      ((layer_PEM(index_bedrock+1)-delta)/(TI_bedrock**2))))
enddo



 do iloop = index_bedrock+2, nsoil_PEM
   do ig = 1,ngrid
      inertiedat_PEM(ig,iloop) = TI_bedrock
   enddo
 enddo

write(*,*) 'PEMETAT0: TI done'

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!b) Soil temperature
    do islope = 1,nslope
!     do ig = 1,ngrid
!        kcond = (TI_PEM(ig,index_breccia+1,islope)*TI_PEM(ig,index_breccia+1,islope))/volcapa
!        tsoil_PEM(ig,index_breccia+1,islope) = tsoil_PEM(ig,index_breccia,islope) + fluxgeo/kcond*(mlayer_PEM(index_breccia)-mlayer_PEM(index_breccia-1))
!
!       do iloop=index_breccia+2,index_bedrock
!            kcond = (TI_PEM(ig,iloop,islope)*TI_PEM(ig,iloop,islope))/volcapa
!            tsoil_PEM(ig,iloop,islope) = tsoil_PEM(ig,index_breccia+1,islope) + fluxgeo/kcond*(mlayer_PEM(iloop-1)-mlayer_PEM(index_breccia))
!      enddo
!        kcond = (TI_PEM(ig,index_bedrock+1,islope)*TI_PEM(ig,index_bedrock+1,islope))/volcapa
!        tsoil_PEM(ig,index_bedrock+1,islope) = tsoil_PEM(ig,index_bedrock,islope) + fluxgeo/kcond*(mlayer_PEM(index_bedrock)-mlayer_PEM(index_bedrock-1))

!       do iloop=index_bedrock+2,nsoil_PEM
!            kcond = (TI_PEM(ig,iloop,islope)*TI_PEM(ig,iloop,islope))/volcapa
!            tsoil_PEM(ig,iloop,islope) = tsoil_PEM(ig,index_bedrock+1,islope) + fluxgeo/kcond*(mlayer_PEM(iloop-1)-mlayer_PEM(index_bedrock))
!      enddo

!       enddo

      call soil_pem_ini(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_ave_yr2(:,islope),tsoil_PEM(:,:,islope))
      call soil_pem_compute(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_ave_yr2(:,islope),tsoil_PEM(:,:,islope))
  
! First raw initialization 
      do it = 1,timelen
        do isoil = nsoil_GCM+1,nsoil_PEM
          tsoil_inst(:,isoil,islope,it) = tsoil_PEM(:,isoil,islope)
        enddo
      enddo

      do it = 1,timelen
        do isoil = nsoil_GCM+1,nsoil_PEM
          call soil_pem_ini(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_ave_yr2(:,islope),tsoil_inst(:,:,islope,it))
        enddo
       enddo

       do isoil = nsoil_GCM+1,nsoil_PEM
         do ig = 1,ngrid
           watersoil_ave(ig,isoil,islope) = exp(beta_clap_h2o/tsoil_PEM(ig,isoil,islope) + alpha_clap_h2o)/tsoil_PEM(ig,isoil,islope)*mmol(igcm_h2o_vap)/(mugaz*r)
         enddo
       enddo
enddo !islope
write(*,*) 'PEMETAT0: TSOIL done'

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!c) Ice table
       ice_table(:,:) = -1.  ! by default, no ice table
       ice_table_thickness(:,:) = -1.
       call computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,watersurf_ave,watersoil_ave,TI_PEM(:,1,:),ice_table,ice_table_thickness)
       call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,tend_h2oglaciers,waterice,global_ave_pressure,ice_table,ice_table_thickness,TI_PEM)
       do islope = 1,nslope
         call soil_pem_ini(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_ave_yr2(:,islope),tsoil_PEM(:,:,islope))
       enddo
       write(*,*) 'PEMETAT0: Ice table done'

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!d) Regolith adsorbed
if (adsorption_pem) then
    m_co2_regolith_phys(:,:,:) = 0.
    m_h2o_regolith_phys(:,:,:) = 0.

    call regolith_adsorption(ngrid,nslope,nsoil_PEM,timelen,tend_h2oglaciers,tend_co2glaciers,waterice,co2ice,tsoil_PEM,TI_PEM,ps_inst,q_co2,q_h2o, &
                             m_h2o_regolith_phys,deltam_h2o_regolith_phys, m_co2_regolith_phys,deltam_co2_regolith_phys)

    deltam_co2_regolith_phys(:) = 0.
    deltam_h2o_regolith_phys(:) = 0.
endif

write(*,*) 'PEMETAT0: CO2 adsorption done'
endif !soil_pem

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!. e) water reservoir
#ifndef CPP_STD
    do ig = 1,ngrid
        if (watercaptag(ig)) then
            water_reservoir(ig) = water_reservoir_nom
        else
            water_reservoir(ig) = 0.
        endif
    enddo
#endif

endif ! of if (startphy_file)

if (soil_pem) then ! Sanity check
    do ig = 1,ngrid
        do islope = 1,nslope
            do iloop = 1,nsoil_PEM
                if (isnan(tsoil_PEM(ig,iloop,islope))) call abort_pem("PEM - pemetat0","NaN detected in Tsoil",1)
            enddo
        enddo
    enddo
endif !soil_pem

END SUBROUTINE pemetat0

END MODULE pemetat0_mod