source: trunk/LMDZ.COMMON/libf/evolution/clim_state_init.F90 @ 4066

MODULE clim_state_init
!-----------------------------------------------------------------------
! NAME
!     clim_state_init
!
! DESCRIPTION
!     Read the start files to initialize the climate state.
!
! AUTHORS & DATE
!     JB Clement, 12/2025
!
! NOTES
!
!-----------------------------------------------------------------------

! DEPENDENCIES
! ------------
use numerics,  only: dp, di, k4
use io_netcdf, only: open_nc, close_nc, get_var_nc, get_dim_nc, put_var_nc, start_name, start1D_name, startfi_name, startpem_name

! DECLARATION
! -----------
implicit none

contains
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

!=======================================================================
SUBROUTINE read_start()
!-----------------------------------------------------------------------
! NAME
!     read_start
!
! DESCRIPTION
!     Read the file "start.nc".
!
! AUTHORS & DATE
!     JB Clement, 12/2025
!
! NOTES
!
!-----------------------------------------------------------------------

! DEPENDENCIES
! ------------
use geometry,   only: nlayer, nlon, nlat, ngrid
use atmosphere, only: set_ap, set_bp, set_ps_PCM, set_teta_PCM
use tracers,    only: nq, qnames, set_q_PCM
use stoppage,   only: stop_clean
use display,    only: print_msg

! DECLARATION
! -----------
implicit none

! LOCAL VARIABLES
! ---------------
real(dp), dimension(nlayer + 1)           :: tmp1d
real(dp), dimension(nlon + 1,nlat)        :: tmp2d
real(dp), dimension(nlon + 1,nlat,nlayer) :: tmp3d
logical(k4)                               :: here
integer(di)                               :: i

! CODE
! ----
! In case of 1D
! ~~~~~~~~~~~~~
if (ngrid == 1) then
    call read_start1D()
    return
end if

! In case of 3D
! ~~~~~~~~~~~~~
! Open
inquire(file = start_name,exist = here)
if (.not. here) call stop_clean(__FILE__,__LINE__,'cannot find required file "'//start_name//'"!',1)
call print_msg('> Reading "'//start_name//'"')
call open_nc(start_name,'read')

! Get the variables
call get_var_nc('ap',tmp1d)
call set_ap(tmp1d)

call get_var_nc('bp',tmp1d)
call set_bp(tmp1d)

call get_var_nc('ps',tmp2d)
call set_ps_PCM(tmp2d)

call get_var_nc('teta',tmp3d)
call set_teta_PCM(tmp3d)

do i = 1,nq
    call get_var_nc(trim(qnames(i)),tmp3d)
    call set_q_PCM(tmp3d,i)
end do

! Close
call close_nc(start_name)

END SUBROUTINE read_start
!=======================================================================

!=======================================================================
SUBROUTINE read_start1D()
!-----------------------------------------------------------------------
! NAME
!     read_start1D
!
! DESCRIPTION
!     Read the file "start1D.txt".
!
! AUTHORS & DATE
!     JB Clement, 12/2025
!
! NOTES
!
!-----------------------------------------------------------------------

! DEPENDENCIES
! ------------
use geometry,   only: nlayer, nslope
use atmosphere, only: set_ap, set_bp, set_ps_PCM, set_teta_PCM, set_u_PCM, set_v_PCM, compute_alt_coord, set_ap, set_bp
use tracers,    only: nq, set_q_PCM
use stoppage,   only: stop_clean
use config,     only: read_callphys
use display,    only: print_msg

! DECLARATION
! -----------
implicit none

! LOCAL VARIABLES
! ---------------
integer(di)                     :: i, j, k, ierr, funit
character(30)                   :: header
real(dp), dimension(1,1)        :: tmp
real(dp), dimension(nslope)     :: tmp_1d
real(dp), dimension(1,1,nlayer) :: q_tmp, teta_tmp, wind_tmp
real(dp), dimension(nlayer + 1) :: ap_tmp, bp_tmp
real(dp)                        :: pa, preff
logical(k4)                     :: here

! CODE
! ----
! Open "start1D.txt"
inquire(file = start1D_name,exist = here)
if (.not. here) call stop_clean(__FILE__,__LINE__,'cannot find required file "'//start1D_name//'"!',1)
call print_msg('> Reading "'//start1D_name//'"')
open(newunit = funit,file = start1D_name,status = "old",action = "read",iostat = ierr)
if (ierr /= 0) call stop_clean(__FILE__,__LINE__,'error opening file "'//start1D_name//'"!',ierr)

! Get the variables
read(funit,*) header, tmp, pa, preff
call set_ps_PCM(tmp)

do i = 1,nq
    read(funit,*,iostat = ierr) header, (tmp_1d(j),j = 1,nslope), (q_tmp(1,1,k),k = 1,nlayer)
    if (ierr /= 0) call stop_clean(__FILE__,__LINE__,'not enough atmospheric layers defined in the file "'//start1D_name//'" for the tracer "'//trim(header)//'"!',1)
    call set_q_PCM(q_tmp,i)
end do
read(funit,*,iostat = ierr) header, (wind_tmp(1,1,k),k = 1,nlayer)
call set_u_PCM(wind_tmp)
read(funit,*,iostat = ierr) header, (wind_tmp(1,1,k),k = 1,nlayer)
call set_v_PCM(wind_tmp)
read(funit,*,iostat = ierr) header, (tmp_1d(j),j = 1,nslope), (teta_tmp(1,1,k),k = 1,nlayer)
call set_teta_PCM(teta_tmp)

! Close
close(funit)

! Construct altitude coordinates (not stored in "start1D.txt")
call compute_alt_coord(pa,preff,ap_tmp,bp_tmp)
call set_ap(ap_tmp)
call set_bp(bp_tmp)

! Read the "callphys.def"
call read_callphys() ! To get 'CO2cond_ps'

END SUBROUTINE read_start1D
!=======================================================================

!=======================================================================
SUBROUTINE read_startfi()
!-----------------------------------------------------------------------
! NAME
!     read_startfi
!
! DESCRIPTION
!     Read the file "startfi.nc".
!
! AUTHORS & DATE
!     JB Clement, 12/2025
!
! NOTES
!
!-----------------------------------------------------------------------

! DEPENDENCIES
! ------------
use geometry,  only: ngrid, nslope, nsoil_PCM
use stoppage,  only: stop_clean
use config,    only: read_controldata
use slopes,    only: set_def_slope_mean, set_subslope_dist, set_iflat
use surface,   only: set_albedodat_PCM, set_albedo_PCM, set_emissivity_PCM
use surf_temp, only: set_tsurf_PCM
use soil_temp, only: set_tsoil_PCM, set_flux_geo_PCM
use frost,     only: set_h2ofrost_PCM, set_co2frost_PCM
use surf_ice,  only: set_is_h2o_perice_PCM, set_co2_perice_PCM
use soil,      only: set_TI_PCM, set_inertiedat_PCM
use display,   only: print_msg

! DECLARATION
! -----------
implicit none

! LOCAL VARIABLES
! ---------------
real(dp), dimension(:),   allocatable       :: tmp1d
real(dp), dimension(:,:), allocatable       :: tmp2d
real(dp), dimension(ngrid,nsoil_PCM,nslope) :: tmp3d
logical(k4)                                 :: here

! CODE
! ----
inquire(file = startfi_name,exist = here)
if (.not. here) call stop_clean(__FILE__,__LINE__,'cannot find required file "'//startfi_name//'"!',1)

! Allocate the array to store the variables
call print_msg('> Reading "'//startfi_name//'"')
allocate(tmp1d(nslope + 1),tmp2d(ngrid,nslope))

! Get control data
call read_controldata()

! Open
call open_nc(startfi_name,'read')

! Get the variables
if (nslope > 1) then
    call get_var_nc('def_slope_mean',tmp1d)
    call set_def_slope_mean(tmp1d)

    call get_var_nc('subslope_dist',tmp2d)
    call set_subslope_dist(tmp2d)
end if

call get_var_nc('flux_geo',tmp2d)
call set_flux_geo_PCM(tmp2d)

call get_var_nc('h2o_ice',tmp2d)
where (tmp2d < 0.) tmp2d = 0._dp ! Remove unphysical values of surface tracer
call set_h2ofrost_PCM(tmp2d)

call get_var_nc('co2',tmp2d)
where (tmp2d < 0.) tmp2d = 0._dp ! Remove unphysical values of surface tracer
call set_co2frost_PCM(tmp2d)

call get_var_nc('perennial_co2ice',tmp2d)
where (tmp2d < 0.) tmp2d = 0._dp ! Remove unphysical values of surface tracer
call set_co2_perice_PCM(tmp2d)

deallocate(tmp1d); allocate(tmp1d(ngrid))
call get_var_nc('watercaptag',tmp1d)
call set_is_h2o_perice_PCM(tmp1d)

call get_var_nc('albedodat',tmp1d)
call set_albedodat_PCM(tmp1d)

call get_var_nc('albedo',tmp2d)
call set_albedo_PCM(tmp2d)

call get_var_nc('emis',tmp2d)
call set_emissivity_PCM(tmp2d)

call get_var_nc('tsurf',tmp2d)
call set_tsurf_PCM(tmp2d)

call get_var_nc('tsoil',tmp3d)
call set_tsoil_PCM(tmp3d)

deallocate(tmp2d); allocate(tmp2d(ngrid,nsoil_PCM))
call get_var_nc('inertiedat',tmp2d)
call set_inertiedat_PCM(tmp2d)

call get_var_nc('inertiesoil',tmp3d)
call set_TI_PCM(tmp3d)

! To do?
!   h2oice_depth
!   d_coef
!   lag_co2_ice

! Close/Deallocate
call close_nc(startfi_name)
deallocate(tmp1d,tmp2d)

END SUBROUTINE read_startfi
!=======================================================================

!=======================================================================
SUBROUTINE read_startpem(tsurf_avg_yr1,tsurf_avg_yr2,ps_avg_glob,ps_ts,q_co2_ts,q_h2o_ts,h2o_surfdensity_avg,d_h2oice,d_co2ice,h2o_ice,co2_ice, &
                         tsoil_avg,h2o_soildensity_avg,icetable_depth,icetable_thickness,ice_porefilling,layerings_map,                         &
                         h2o_ads_reg,co2_ads_reg,delta_h2o_ads,delta_co2_ads)
!-----------------------------------------------------------------------
! NAME
!     read_startpem
!
! DESCRIPTION
!     Read the file "startpem.nc" which stores the PEM state.
!
! AUTHORS & DATE
!     L. Lange, 2023
!     JB Clement, 2023-2026
!
! NOTES
!
!-----------------------------------------------------------------------

! DEPENDENCIES
! ------------
use stoppage,           only: stop_clean
use geometry,           only: ngrid, nslope, nsoil, nsoil_PCM
use evolution,          only: dt
use physics,            only: mugaz, r
use planet,             only: TI_breccia, TI_bedrock, alpha_clap_h2o, beta_clap_h2o
use frost,              only: h2ofrost_PCM, h2o_frost4PCM, co2frost_PCM, co2_frost4PCM
use surf_ice,           only: is_h2o_perice_PCM, h2oice_huge_ini, co2_perice_PCM
use soil,               only: do_soil, TI, depth_breccia, depth_bedrock, index_breccia, index_bedrock, inertiedat, layer, inertiedat_PCM
use soil_temp,          only: ini_tsoil_pem, compute_tsoil
use soil_therm_inertia, only: update_soil_TI
use ice_table,          only: icetable_equilibrium, icetable_dynamic, computeice_table_equilibrium
use sorption,           only: do_sorption, evolve_regolith_adsorption
use tracers,            only: mmol, iPCM_qh2o
use layered_deposits,   only: layering, do_layering, array2map, ini_layering, add_stratum
use surf_ice,           only: rho_co2ice, rho_h2oice
use display,            only: print_msg
use utility,            only: int2str

! DECLARATION
! -----------
implicit none

! ARGUMENTS
! ---------
real(dp), dimension(:,:),       intent(in)    :: tsurf_avg_yr1, tsurf_avg_yr2 ! Surface temperature at the last but one PCM run and at the last PCM run
real(dp),                       intent(in)    :: ps_avg_glob                  ! Mean average pressure on the planet
real(dp), dimension(:,:),       intent(in)    :: ps_ts                        ! Surface pressure timeseries
real(dp), dimension(:,:),       intent(in)    :: q_co2_ts, q_h2o_ts           ! MMR of CO2 and H2O
real(dp), dimension(:,:),       intent(in)    :: h2o_surfdensity_avg          ! Average of surface water density
real(dp), dimension(:,:),       intent(in)    :: d_h2oice, d_co2ice           ! Tendencies on ice
real(dp), dimension(:,:),       intent(inout) :: h2o_ice, co2_ice             ! Surface ice
real(dp), dimension(:,:,:),     intent(inout) :: tsoil_avg                    ! Soil temperature
real(dp), dimension(:,:,:),     intent(inout) :: h2o_soildensity_avg          ! Average of soil water soil density
real(dp), dimension(:,:),       intent(inout) :: icetable_depth               ! Ice table depth
real(dp), dimension(:,:),       intent(inout) :: icetable_thickness           ! Ice table thickness
real(dp), dimension(:,:,:),     intent(inout) :: ice_porefilling              ! Subsurface ice pore filling
real(dp), dimension(:,:,:),     intent(inout) :: h2o_ads_reg, co2_ads_reg     ! Mass of H2O and CO2 adsorbed
type(layering), dimension(:,:), intent(inout) :: layerings_map                ! Layerings
real(dp), dimension(:),         intent(inout) :: delta_h2o_ads, delta_co2_ads ! Mass of H2O and CO2 exchanged due to adsorption/desorption

! LOCAL VARIABLES
! ---------------
logical(k4)                               :: here
integer(di)                               :: i, islope, k, nb_str_max, nsoil_startpem
real(dp)                                  :: delta           ! Depth of the interface regolith/breccia, breccia/bedrock [m]
real(dp), dimension(ngrid,nsoil,nslope)   :: TI_startpem     ! Soil thermal inertia saved in the startpem [SI]
real(dp), dimension(ngrid,nsoil,nslope)   :: tsoil_startpem  ! Soil temperature saved in the startpem [K]
real(dp), dimension(:,:,:,:), allocatable :: layerings_array ! Array for layerings

! CODE
! ----
! Check if the file exists
inquire(file = startpem_name,exist = here)

! If the file is not here
! ~~~~~~~~~~~~~~~~~~~~~~~
if (.not. here) then
    ! It is possibly because it is the very first PEM run so everything needs to be initalized
    call print_msg('> The file "'//startpem_name//'" was not found (possibly because this is the first PEM run)')

    ! H2O ice
    call print_msg("'h2o_ice' is initialized with default value 'h2oice_huge_ini' where 'watercaptag' is true.")
    h2o_ice(:,:) = 0._dp
    do i = 1,ngrid
        if (is_h2o_perice_PCM(i)) h2o_ice(i,:) = h2oice_huge_ini + h2ofrost_PCM(i,:) - h2o_frost4PCM(i,:)
    end do

    ! CO2 ice
    call print_msg("'co2_ice' is initialized with 'perennial_co2ice' found in the PCM.")
    co2_ice(:,:) = co2_perice_PCM(:,:) + co2frost_PCM(:,:) - co2_frost4PCM(:,:)

    if (do_soil) then
        ! Thermal inertia
        do islope = 1,nslope
            do i = 1,ngrid
                if (TI(i,index_breccia,islope) < TI_breccia) then
                    !!! transition
                    delta = depth_breccia
                    TI(i,index_breccia + 1,islope) = sqrt((layer(index_breccia + 1) - layer(index_breccia))/                       &
                                                               (((delta - layer(index_breccia))/(TI(i,index_breccia,islope)**2)) + &
                                                               ((layer(index_breccia + 1) - delta)/(TI_breccia**2))))
                    do k = index_breccia + 2,index_bedrock
                        TI(i,k,islope) = TI_breccia
                    end do
                else ! we keep the high ti values
                    do k = index_breccia + 1,index_bedrock
                        TI(i,k,islope) = TI(i,index_breccia,islope)
                    end do
                end if
                !! transition
                delta = depth_bedrock
                TI(i,index_bedrock + 1,islope) = sqrt((layer(index_bedrock + 1) - layer(index_bedrock))/                   &
                                                       (((delta - layer(index_bedrock))/(TI(i,index_bedrock,islope)**2)) + &
                                                       ((layer(index_bedrock + 1) - delta)/(TI_breccia**2))))
                do k = index_bedrock + 2,nsoil
                    TI(i,k,islope) = TI_bedrock
                end do
            end do
        end do

        do k = 1,nsoil_PCM
            inertiedat(:,k) = inertiedat_PCM(:,k)
        end do
        !!! zone de transition
        delta = depth_breccia
        do i = 1,ngrid
            if (inertiedat(i,index_breccia) < TI_breccia) then
                inertiedat(i,index_breccia + 1) = sqrt((layer(index_breccia + 1) - layer(index_breccia))/                    &
                                                        (((delta - layer(index_breccia))/(inertiedat(i,index_breccia)**2)) + &
                                                        ((layer(index_breccia + 1) - delta)/(TI_breccia**2))))
                do k = index_breccia + 2,index_bedrock
                    inertiedat(i,k) = TI_breccia
                end do
            else
                do k = index_breccia + 1,index_bedrock
                    inertiedat(i,k) = inertiedat(i,index_breccia)
                end do
            end if
        end do
        !!! zone de transition
        delta = depth_bedrock
        do i = 1,ngrid
            inertiedat(i,index_bedrock + 1) = sqrt((layer(index_bedrock + 1) - layer(index_bedrock))/                    &
                                                    (((delta - layer(index_bedrock))/(inertiedat(i,index_bedrock)**2)) + &
                                                    ((layer(index_bedrock + 1) - delta)/(TI_bedrock**2))))
        end do

        do k = index_bedrock + 2,nsoil
            do i = 1,ngrid
                inertiedat(i,k) = TI_bedrock
            end do
        end do

        ! Soil temperature
        do islope = 1,nslope
            call ini_tsoil_pem(ngrid,nsoil,TI(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_avg(:,:,islope))
            call compute_tsoil(ngrid,nsoil,.true.,TI(:,:,islope),dt,tsurf_avg_yr2(:,islope),tsoil_avg(:,:,islope))

            ! First raw initialization
            h2o_soildensity_avg(:,nsoil_PCM + 1:nsoil,islope) = exp(beta_clap_h2o/tsoil_avg(:,nsoil_PCM + 1:nsoil,islope) + alpha_clap_h2o)/tsoil_avg(:,nsoil_PCM + 1:nsoil,islope)**mmol(iPCM_qh2o)/(mugaz*r)
        end do

        ! Ice table
        if (icetable_equilibrium) then
            call computeice_table_equilibrium(is_h2o_perice_PCM,h2o_surfdensity_avg,h2o_soildensity_avg,TI(:,1,:),icetable_depth,icetable_thickness)
            call update_soil_TI(d_h2oice,h2o_ice,ps_avg_glob,icetable_depth,icetable_thickness,ice_porefilling,icetable_equilibrium,icetable_dynamic,TI)
            do islope = 1,nslope
                call ini_tsoil_pem(ngrid,nsoil,TI(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_avg(:,:,islope))
            end do
        else if (icetable_dynamic) then
            call update_soil_TI(d_h2oice,h2o_ice,ps_avg_glob,icetable_depth,icetable_thickness,ice_porefilling,icetable_equilibrium,icetable_dynamic,TI)
            do islope = 1,nslope
                call ini_tsoil_pem(ngrid,nsoil,TI(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_avg(:,:,islope))
            end do
        end if

        ! Absorption in regolith
        if (do_sorption) then
            h2o_ads_reg = 0._dp
            co2_ads_reg = 0._dp
            call evolve_regolith_adsorption(d_h2oice,d_co2ice,h2o_ice,co2_ice,tsoil_avg,TI,ps_ts,q_h2o_ts,q_co2_ts,h2o_ads_reg,co2_ads_reg,delta_h2o_ads,delta_co2_ads)
            delta_co2_ads = 0._dp
            delta_h2o_ads = 0._dp
        end if ! do_sorption
    end if ! do_soil

    ! Layering
    if (do_layering) then
        call print_msg('layerings_map is initialized with sub-surface strata.')
        call print_msg("Ice is added with 'h2oice_huge_ini' where 'watercaptag' is true and otherwise with 'perennial_co2ice' found in the PCM.")
        do i = 1,ngrid
            if (is_h2o_perice_PCM(i)) then
                do islope = 1,nslope
                    call ini_layering(layerings_map(i,islope))
                    call add_stratum(layerings_map(i,islope),1.05_dp*h2oice_huge_ini/rho_h2oice,0._dp,h2oice_huge_ini/rho_h2oice,0.05_dp*h2oice_huge_ini/rho_h2oice,0._dp,0._dp)
                end do
            else
                do islope = 1,nslope
                    call ini_layering(layerings_map(i,islope))
                    if (co2_perice_PCM(i,islope) > 0.) call add_stratum(layerings_map(i,islope),1.05_dp*co2_perice_PCM(i,islope)/rho_co2ice,co2_perice_PCM(i,islope)/rho_co2ice,0._dp,0.05_dp*co2_perice_PCM(i,islope)/rho_co2ice,0._dp,0._dp)
                end do
            end if
        end do
    end if ! do_layering

    return
end if

! If the file is present
! ~~~~~~~~~~~~~~~~~~~~~~
call print_msg('> Reading "'//startpem_name//'"')
call open_nc(startpem_name,'read')

! H2O ice
h2o_ice(:,:) = 0._dp
call get_var_nc('h2o_ice',h2o_ice(:,:))

! CO2 ice
co2_ice(:,:) = 0._dp
call get_var_nc('co2_ice',co2_ice(:,:))

if (do_soil) then
    ! Check if the number of soil layers is compatible
    call get_dim_nc('subsurface_layers',nsoil_startpem)
    if (nsoil_startpem /= nsoil) then
        call print_msg('nsoil (PEM) = '//int2str(nsoil)//' | nsoil ("'//startpem_name//'") = '//int2str(nsoil_startpem))
        call stop_clean(__FILE__,__LINE__,'nsoil defined in the PEM is different from the one read in "'//startpem_name//'"!',1)
    end if

    do islope = 1,nslope
        ! Soil temperature
        call get_var_nc('tsoil',tsoil_startpem(:,:,islope))
        ! Predictor-corrector: due to changes of surface temperature in the PCM, the tsoil_startpem is adapted firstly
        !                      for PCM year 1 and then for PCM year 2 in order to build the evolution of the profile at depth
        call compute_tsoil(ngrid,nsoil,.true.,TI(:,:,islope),dt,tsurf_avg_yr1(:,islope),tsoil_startpem(:,:,islope))
        call compute_tsoil(ngrid,nsoil,.false.,TI(:,:,islope),dt,tsurf_avg_yr1(:,islope),tsoil_startpem(:,:,islope))
        call compute_tsoil(ngrid,nsoil,.false.,TI(:,:,islope),dt,tsurf_avg_yr2(:,islope),tsoil_startpem(:,:,islope))
    end do
    tsoil_avg(:,nsoil_PCM + 1:nsoil,:) = tsoil_startpem(:,nsoil_PCM + 1:nsoil,:)
    if (any(isnan(tsoil_avg(:,:,:)))) call stop_clean(__FILE__,__LINE__,"NaN detected in 'tsoil_avg'",1)

    ! Thermal inertia
    call get_var_nc('TI',TI_startpem(:,:,:))
    TI(:,nsoil_PCM + 1:nsoil,:) = TI_startpem(:,nsoil_PCM + 1:nsoil,:) ! 1st layers can change because of the presence of ice at the surface, so we don't change it here
    call get_var_nc('inertiedat',inertiedat(:,:))

    ! Ice table
    call get_var_nc('icetable_depth',icetable_depth(:,:))
    if (icetable_dynamic) call get_var_nc('ice_porefilling',ice_porefilling(:,:,:))
    if (icetable_equilibrium) then
        call get_var_nc('icetable_thickness',icetable_thickness(:,:))
    else if (icetable_dynamic) then
        call get_var_nc('ice_porefilling',ice_porefilling(:,:,:))
    end if

    ! Absorption in regolith
    if (do_sorption) then
        call get_var_nc('co2_ads_reg',co2_ads_reg(:,:,:))
        call get_var_nc('h2o_ads_reg',h2o_ads_reg(:,:,:))
        call evolve_regolith_adsorption(d_h2oice,d_co2ice,h2o_ice,co2_ice,tsoil_avg,TI,ps_ts,q_h2o_ts,q_co2_ts,h2o_ads_reg,co2_ads_reg,delta_h2o_ads,delta_co2_ads)
    end if ! do_sorption
end if ! do_soil

! Layering
if (do_layering) then
    call get_dim_nc('nb_str_max',nb_str_max)
    allocate(layerings_array(ngrid,nslope,nb_str_max,6))
    layerings_array = 0.
    call get_var_nc('stratif_top_elevation',layerings_array(:,:,:,1))
    call get_var_nc('stratif_h_co2ice',layerings_array(:,:,:,2))
    call get_var_nc('stratif_h_h2oice',layerings_array(:,:,:,3))
    call get_var_nc('stratif_h_dust',layerings_array(:,:,:,4))
    call get_var_nc('stratif_h_pore',layerings_array(:,:,:,5))
    call get_var_nc('stratif_poreice_volfrac',layerings_array(:,:,:,6))
    call array2map(layerings_array,layerings_map)
    deallocate(layerings_array)
end if

! Close
call close_nc(startpem_name)

END SUBROUTINE read_startpem
!=======================================================================

END MODULE clim_state_init