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

PROGRAM pem
!-----------------------------------------------------------------------
! NAME
!     pem
!
! DESCRIPTION
!     Main entry point for the Planetary Evolution Model (PEM).
!
! AUTHORS & DATE
!     R. Vandemeulebrouck, 22/07/2022
!     L. Lange, 22/07/2022
!     JB Clement, 2023-2025
!
! NOTES
!
!-----------------------------------------------------------------------

! DEPENDENCIES
! ------------
! Common modules
use job_timelimit_mod, only: timelimit, antetime, timewall
use parse_args_mod,    only: parse_args
! PEM modules
use allocation,         only: ini_allocation, end_allocation
use atmosphere,         only: ps_PCM, evolve_pressure, CO2cond_ps_PCM, build4PCM_atmosphere
use clim_state_init,    only: read_start, read_startfi, read_startpem
use clim_state_rec,     only: write_restart, write_restartfi, write_restartpem
use config,             only: read_rundef
use display,            only: print_ini, print_end, print_msg
use evolution,          only: n_yr_run, n_yr_sim, nmax_yr_sim, nmax_yr_run, dt, idt, r_plnt2earth_yr, pem_ini_date
use geometry,           only: ngrid, nslope, nday, nsoil_PCM, nsoil, cell_area, total_surface, nlayer
use glaciers,           only: h2oice_flow, co2ice_flow, flow_co2glaciers, flow_h2oglaciers
use ice_table,          only: icetable_equilibrium, icetable_dynamic, icetable_depth, icetable_thickness, ice_porefilling, evolve_ice_table
use info,               only: read_info, update_info
use layered_deposits,   only: layering, do_layering, del_layering, evolve_layering, ptrarray, layering2surfice, surfice2layering
use maths,              only: pi
use numerics,           only: dp, qp, di, li, k4, minieps, minieps_qp
use orbit,              only: evol_orbit, read_orbitpm, compute_maxyr_orbit, update_orbit
use output,             only: write_diagevol, dim_ngrid, dim_nsoil
use physics,            only: g
use slopes,             only: subslope_dist, def_slope_mean
use soil,               only: do_soil, set_soil, TI, build4PCM_soil
use soil_temp,          only: tsoil_PCM, shift_tsoil2surf, evolve_soil_temp
use soil_therm_inertia, only: update_soil_TI
use sorption,           only: do_sorption, compute_totmass_adsorbed, evolve_regolith_adsorption
use stopping_crit,      only: stopFlags, stopping_crit_pressure, stopping_crit_h2o, stopping_crit_h2oice, stopping_crit_co2ice
use surface,            only: emissivity_PCM, build4PCM_surf_rad_prop
use surf_ice,           only: evolve_co2ice, evolve_h2oice, balance_h2oice_reservoirs, build4PCM_perice
use surf_temp,          only: tsurf_PCM, adapt_tsurf2disappearedice, build4PCM_tsurf
use tendencies,         only: evolve_tend_co2, evolve_tend_h2o
use tracers,            only: adapt_tracers2pressure, build4PCM_tracers, nq
use utility,            only: real2str
use xios_data,          only: load_xios_data

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

! LOCAL VARIABLES
! ---------------
! Utility-related:
integer(li)               :: cr     ! Number of clock ticks per second (count rate)
integer(li)               :: c1, c2 ! Counts of processor clock
character(:), allocatable :: num    ! To write slope variables
integer(di)               :: i, islope
! Pressure-related:
real(dp), dimension(:),   allocatable :: ps_avg          ! Average surface pressure [Pa]
real(dp), dimension(:),   allocatable :: ps_dev          ! Deviation of surface pressure [Pa]
real(dp), dimension(:,:), allocatable :: ps_ts           ! Surface pressure timeseries [Pa]
real(dp)                              :: ps_avg_glob_ini ! Global average pressure at initialization [Pa]
real(dp)                              :: ps_avg_glob_old ! Global average pressure of previous time step [Pa]
real(dp)                              :: ps_avg_glob     ! Global average pressure of current time step [Pa]
real(dp), dimension(:),   allocatable :: ps4PCM          ! Surface pressure reconstruction to feed back into PCM [Pa]
real(dp), dimension(:,:), allocatable :: teta4PCM        ! Potential temperature reconstruction to feed back into PCM [K]
real(dp), dimension(:,:), allocatable :: air_mass4PCM    ! Air mass reconstruction to feed back into PCM [kg]
real(dp)                              :: pa4PCM          ! Transition pressure (for hybrid coord.) [Pa]
real(dp)                              :: preff4PCM       ! Reference pressure [Pa]
! Ice-related:
real(dp),    dimension(:,:), allocatable :: h2o_ice                    ! H2O ice [kg.m-2]
real(dp),    dimension(:,:), allocatable :: co2_ice                    ! CO2 ice [kg.m-2]
real(dp)                                 :: h2oice_sublim_coverage_ini ! Initial surface area of sublimating H2O ice [m2]
real(dp)                                 :: co2ice_sublim_coverage_ini ! Initial surface area of sublimating CO2 ice [m2]
logical(k4), dimension(:,:), allocatable :: is_h2oice_ini              ! Initial location of H2O ice
logical(k4), dimension(:,:), allocatable :: is_co2ice_ini              ! Initial location of CO2 ice
logical(k4), dimension(:,:), allocatable :: is_co2ice_flow             ! Flag for location of CO2 glacier flow
logical(k4), dimension(:,:), allocatable :: is_h2oice_flow             ! Flag for location of H2O glacier flow
real(dp),    dimension(:,:), allocatable :: h2o_ice4PCM                ! H2O ice reconstruction to feed back into PCM [kg.m-2]
real(dp),    dimension(:,:), allocatable :: co2_ice4PCM                ! CO2 ice reconstruction to feed back into PCM [kg.m-2]
logical(k4), dimension(:),   allocatable :: is_h2o_perice              ! Location of H2O infinite reservoir, called 'watercaptag' in PCM
logical(k4), dimension(:,:), allocatable :: is_co2ice_disappeared      ! Flag to check if CO2 ice disappeared at the previous timestep
! Surface-related:
real(dp), dimension(:,:), allocatable :: tsurf_avg           ! Average surface temperature [K]
real(dp), dimension(:,:), allocatable :: tsurf_avg_yr1       ! Average surface temperature of the second to last PCM run [K]
real(dp), dimension(:,:), allocatable :: tsurf_dev           ! Deviation of surface temperature [K]
real(dp), dimension(:,:), allocatable :: h2o_surfdensity_avg ! Average water surface density [kg/m^3]
real(dp), dimension(:,:), allocatable :: zshift_surf         ! Elevation shift for the surface [m]
real(dp), dimension(:,:), allocatable :: zlag                ! Newly built lag thickness [m]
real(dp), dimension(:,:), allocatable :: albedo4PCM          ! Albedo reconstruction to feed back into PCM
real(dp), dimension(:,:), allocatable :: emissivity4PCM      ! Emissivity reconstruction to feed back into PCM
real(dp), dimension(:,:), allocatable :: tsurf4PCM           ! Surface temperature reconstruction to feed back into PCM [K]
! Soil-related:
real(dp), dimension(:,:,:),   allocatable :: tsoil_avg                ! Average soil temperature [K]
real(dp), dimension(:,:,:),   allocatable :: tsoil_dev                ! Deviation pf soil temperature [K]
real(dp), dimension(:,:,:,:), allocatable :: tsoil_ts                 ! Soil temperature timeseries [K]
real(dp), dimension(:,:,:,:), allocatable :: tsoil_ts_old             ! Soil temperature timeseries at the previous time step [K]
real(dp), dimension(:,:,:),   allocatable :: h2o_soildensity_avg      ! Average of soil water soil density [kg/m^3]
real(dp), dimension(:),       allocatable :: delta_co2_ads            ! Quantity of CO2 exchanged due to adsorption/desorption [kg/m^2]
real(dp), dimension(:),       allocatable :: delta_h2o_ads            ! Quantity of H2O exchanged due to adsorption/desorption [kg/m^2]
real(qp)                                  :: totmass_adsh2o           ! Total mass of H2O exchanged because of adsorption/desorption [kg]
real(dp), dimension(:,:,:),   allocatable :: inertiesoil4PCM          ! Soil thermal inertia reconstruction to feed back into PCM [SI]
real(dp), dimension(:,:,:),   allocatable :: tsoil4PCM                ! Average soil temperature reconstruction to feed back into PCM [K]
real(dp), dimension(:,:),     allocatable :: flux_geo4PCM             ! Geothermal flux reconstruction to feed back into PCM [W/m2]
real(dp), dimension(:,:,:),   allocatable :: h2o_ads_reg              ! H2O adsorbed in the regolith [kg/m^2]
real(dp), dimension(:,:,:),   allocatable :: co2_ads_reg              ! CO2 adsorbed in the regolith [kg/m^2]
real(dp), dimension(:,:),     allocatable :: icetable_depth_old       ! Old depth of the ice table [m]
real(dp), dimension(:),       allocatable :: delta_icetable ! Total mass of the H2O that has sublimated / condenses from the ice table [kg]
! Tracer-related:
real(dp), dimension(:,:),   allocatable :: q_co2_ts     ! CO2 mass mixing ratio in the first layer [kg/kg]
real(dp), dimension(:,:),   allocatable :: q_co2_ts_ini ! Initial CO2 mass mixing ratio in the first layer [kg/kg]
real(dp), dimension(:,:),   allocatable :: q_h2o_ts     ! H2O mass mixing ratio in the first layer [kg/kg]
real(dp), dimension(:,:,:), allocatable :: q4PCM        ! Tracers reconstruction to feed back into PCM [kg/kg]
! Tendency-related:
real(dp), dimension(:,:), allocatable :: d_co2ice     ! Tendency of perennial CO2 ice [kg/m2/y]
real(dp), dimension(:,:), allocatable :: d_co2ice_ini ! Tendency of perennial CO2 ice at the beginning [kg/m2/y]
real(dp), dimension(:,:), allocatable :: d_h2oice     ! Tendency of perennial H2O ice [kg/m2/y]
real(dp), dimension(:,:), allocatable :: d_h2oice_new ! Adjusted tendency of perennial H2O ice to keep balance between donor and recipient [kg/m2/y]
! Layering-related:
type(layering), dimension(:,:), allocatable :: layerings_map    ! Layering for each grid point and slope
type(ptrarray), dimension(:,:), allocatable :: current          ! Current active stratum in the layering
real(dp),       dimension(:,:), allocatable :: h2oice_depth     ! Depth of subsurface ice layer
real(dp),       dimension(:,:), allocatable :: h2oice_depth_old ! Old depth of subsurface ice layer
logical(k4),    dimension(:,:), allocatable :: new_str, new_lag ! Flags for the layering algorithm
! Evolution-related:
real(dp)        :: nmax_yr_runorb ! Maximum number of years for the run due to orbital parameters
type(stopFlags) :: stopcrit       ! Stopping criteria
! Balance-related
real(qp) :: totmass_co2ice, totmass_atmco2, totmass_adsco2             ! Current total CO2 masses (surface ice|atmospheric|adsorbed)
real(qp) :: totmass_co2ice_ini, totmass_atmco2_ini, totmass_adsco2_ini ! Initial total CO2 masses (surface ice|atmospheric|adsorbed)
real(qp) :: totmass_ini                                                ! Initial total CO2 mass [kg]
real(qp) :: S_atm_2_h2o, S_h2o_2_atm, S_atm_2_h2oice, S_h2oice_2_atm   ! Variables to balance H2O ice reservoirs

! CODE
! ----
! Pre-processing step
!~~~~~~~~~~~~~~~~~~~~
! Elapsed time with system clock
call system_clock(count_rate = cr)
call system_clock(c1)

! Parse command-line options
call parse_args()

! Initialization
! ~~~~~~~~~~~~~~
! Header
call print_ini()

! Allocate module arrays
call ini_allocation()

! Read the duration information of the simulation
call read_info()

! Read the PEM parameters
call read_rundef()

! Read the orbital parameters
call read_orbitpm(n_yr_sim)

! Read the "start.nc"
call read_start()

! Read the "startfi.nc"
call read_startfi()

! Read the PCM data given by XIOS
allocate(ps_avg(ngrid),ps_ts(ngrid,nday))
allocate(tsurf_avg(ngrid,nslope),tsurf_avg_yr1(ngrid,nslope),h2o_surfdensity_avg(ngrid,nslope))
allocate(tsoil_avg(ngrid,nsoil,nslope),tsoil_ts(ngrid,nsoil,nslope,nday),h2o_soildensity_avg(ngrid,nsoil,nslope))
allocate(q_h2o_ts(ngrid,nday),q_co2_ts(ngrid,nday))
allocate(d_h2oice(ngrid,nslope),d_co2ice(ngrid,nslope))

call load_xios_data(ps_avg,ps_ts,tsurf_avg,tsurf_avg_yr1,tsoil_avg,tsoil_ts,h2o_surfdensity_avg,h2o_soildensity_avg, &
                    q_h2o_ts,q_co2_ts,d_h2oice,d_co2ice)

! Initiate soil settings and TI
if (do_soil) call set_soil(TI)

! Compute the deviation from the average
allocate(ps_dev(ngrid),tsurf_dev(ngrid,nslope),tsoil_dev(ngrid,nsoil_PCM,nslope))
ps_dev(:) = ps_PCM(:) - ps_avg(:)
tsurf_dev(:,:) = tsurf_PCM(:,:) - tsurf_avg(:,:)
tsoil_dev(:,:,:) = tsoil_PCM(:,:,:) - tsoil_avg(:,1:nsoil_PCM,:)

! Compute global surface pressure
ps_avg_glob = sum(cell_area*ps_avg)/total_surface

! Compute the accepted maximum number of years due to orbital parameters (if needed)
call compute_maxyr_orbit(n_yr_sim,nmax_yr_runorb)

! Read the "startpem.nc"
allocate(h2o_ice(ngrid,nslope),co2_ice(ngrid,nslope))
allocate(h2o_ads_reg(ngrid,nsoil,nslope),co2_ads_reg(ngrid,nsoil,nslope),delta_h2o_ads(ngrid),delta_co2_ads(ngrid))
allocate(layerings_map(ngrid,nslope))
call read_startpem(tsurf_avg_yr1,tsurf_avg,ps_avg_glob_ini,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)
deallocate(tsurf_avg_yr1)

! Save initial set-up useful for the next computations
call print_msg('> Saving some initial climate state variables')
allocate(d_co2ice_ini(ngrid,nslope),q_co2_ts_ini(ngrid,nday))
allocate(is_h2oice_ini(ngrid,nslope),is_co2ice_ini(ngrid,nslope))
ps_avg_glob_ini = ps_avg_glob
d_co2ice_ini(:,:) = d_co2ice(:,:)
q_co2_ts_ini(:,:) = q_co2_ts(:,:)
h2oice_sublim_coverage_ini = 0._dp
co2ice_sublim_coverage_ini = 0._dp
is_h2oice_ini(:,:) = .false.
is_co2ice_ini(:,:) = .false.
totmass_co2ice_ini = 0._qp
totmass_atmco2_ini = 0._qp
totmass_adsco2_ini = 0._qp
totmass_adsh2o = 0._qp
do i = 1,ngrid
    totmass_atmco2_ini = totmass_atmco2_ini + cell_area(i)*ps_avg(i)/g
    do islope = 1,nslope
        totmass_co2ice_ini = totmass_co2ice_ini + co2_ice(i,islope)*cell_area(i)*subslope_dist(i,islope)/cos(pi*def_slope_mean(islope)/180._dp)
        if (co2_ice(i,islope) > 0._dp) then
            is_co2ice_ini(i,islope) = .true.
            if (d_co2ice(i,islope) < 0._dp) co2ice_sublim_coverage_ini = co2ice_sublim_coverage_ini + cell_area(i)*subslope_dist(i,islope)/cos(pi*def_slope_mean(islope)/180._dp)
        end if
        if (h2o_ice(i,islope) > 0._dp) then
            is_h2oice_ini(i,islope) = .true.
            if (d_co2ice(i,islope) < 0._dp) h2oice_sublim_coverage_ini = h2oice_sublim_coverage_ini + cell_area(i)*subslope_dist(i,islope)/cos(pi*def_slope_mean(islope)/180._dp)
        end if
    end do
end do
call print_msg("Initial global average pressure [Pa] = "//real2str(ps_avg_glob_ini))
call print_msg("Initial surface area of sublimating CO2 ice [m2] = "//real2str(co2ice_sublim_coverage_ini))
call print_msg("Initial surface area of sublimating H2O ice [m2] = "//real2str(h2oice_sublim_coverage_ini))
if (do_sorption) call compute_totmass_adsorbed(h2o_ads_reg,co2_ads_reg,totmass_adsco2_ini,totmass_adsh2o)

! Main evolution loop
! ~~~~~~~~~~~~~~~~~~~
call print_msg('')
call print_msg('********* Evolution *********')
if (do_layering) then
    allocate(h2oice_depth(ngrid,nslope),h2oice_depth_old(ngrid,nslope),new_str(ngrid,nslope),new_lag(ngrid,nslope),current(ngrid,nslope))
    new_str = .true.
    new_lag = .true.
    do islope = 1,nslope
        do i = 1,ngrid
            current(i,islope)%p => layerings_map(i,islope)%top
        end do
    end do
    ! Conversion to surface ice
    call layering2surfice(layerings_map,h2o_ice,co2_ice,h2oice_depth)
    ! We put the sublimating tendency coming from subsurface ice into the overall tendency
    !where (h2oice_depth > 0. .and. zdqsdif_ssi_tot < -minieps) d_h2oice = zdqsdif_ssi_tot
end if
if (nslope == 1) then ! No slope
    allocate(character(0) :: num)
else ! Using slopes
    allocate(character(8) :: num)
end if
allocate(delta_icetable(ngrid),icetable_depth_old(ngrid,nslope),is_co2ice_disappeared(ngrid,nslope))
is_co2ice_disappeared(:,:) = .false.
delta_icetable(:) = 0._dp
n_yr_run = 0
idt = 0
do while (n_yr_run < min(nmax_yr_runorb,nmax_yr_run) .and. n_yr_sim < nmax_yr_sim)
    call print_msg('**** Iteration of the PEM leg (Planetary years): '//real2str(n_yr_run + dt))
    ! Evolve global surface pressure
    call evolve_pressure(d_co2ice,delta_co2_ads,do_sorption,ps_avg_glob_old,ps_avg_glob,ps_avg)

    ! Adapt tracers according to global surface pressure
    call adapt_tracers2pressure(ps_avg_glob_old,ps_avg_glob,ps_ts,q_h2o_ts,q_co2_ts)

    ! Evolve surface ice
    allocate(zshift_surf(ngrid,nslope),zlag(ngrid,nslope))
    if (do_layering) then
        h2oice_depth_old(:,:) = h2oice_depth(:,:)
        do islope = 1,nslope
            do i = 1,ngrid
                call evolve_layering(layerings_map(i,islope),d_co2ice(i,islope),d_h2oice_new(i,islope),new_str(i,islope),zshift_surf(i,islope),new_lag(i,islope),zlag(i,islope),current(i,islope)%p)
                !call print_layering(layerings_map(i,islope))
            end do
        end do
        ! Conversion to surface ice
        call layering2surfice(layerings_map,h2o_ice,co2_ice,h2oice_depth)

    else
        zlag(:,:) = 0._dp
        call evolve_h2oice(delta_h2o_ads,delta_icetable,h2o_ice,d_h2oice,zshift_surf,stopcrit)
        call evolve_co2ice(co2_ice,d_co2ice,zshift_surf)
    end if

    if (do_layering) then
        allocate(d_h2oice_new(ngrid,nslope))
        call stopping_crit_h2o(delta_h2o_ads,delta_icetable,h2o_ice,d_h2oice,S_atm_2_h2o,S_h2o_2_atm,S_atm_2_h2oice,S_h2oice_2_atm,stopcrit)
        call balance_h2oice_reservoirs(S_atm_2_h2o,S_h2o_2_atm,S_atm_2_h2oice,S_h2oice_2_atm,h2o_ice,d_h2oice,d_h2oice_new)
        deallocate(d_h2oice_new)
    end if

    ! Flow glaciers according to surface ice
    if (co2ice_flow) then
        allocate(is_co2ice_flow(ngrid,nslope))
        call flow_co2glaciers(q_co2_ts,ps_ts,ps_avg_glob_old,ps_avg_glob,co2_ice,is_co2ice_flow)
    end if
    if (h2oice_flow) then
        allocate(is_h2oice_flow(ngrid,nslope))
        call flow_h2oglaciers(tsurf_avg,h2o_ice,is_h2oice_flow)
    end if
    if (do_layering) call surfice2layering(h2o_ice,co2_ice,layerings_map)

    ! Adapt surface temperature if surface ice disappeared
    call adapt_tsurf2disappearedice(co2_ice,is_co2ice_ini,is_co2ice_disappeared,tsurf_avg)
    !call adapt_tsurf2disappearedice(h2o_ice,is_h2oice_ini,tsurf_avg)

    if (do_soil) then
        ! Shift soil temperature to surface
        call shift_tsoil2surf(ngrid,nsoil,nslope,zshift_surf,zlag,tsurf_avg,tsoil_avg)

        ! Evolve soil temperature
        allocate(tsoil_ts_old(ngrid,nsoil,nslope,nday))
        call evolve_soil_temp(tsoil_avg,tsurf_avg,tsoil_ts,tsoil_ts_old,h2o_soildensity_avg)

        ! Evolve ice table
        call evolve_ice_table(h2o_ice,h2o_surfdensity_avg,h2o_soildensity_avg,tsoil_avg,tsurf_avg,delta_icetable,q_h2o_ts, &
                              ps_avg,icetable_depth,icetable_thickness,ice_porefilling,icetable_depth_old)

        ! Update soil thermal properties according to ice table and soil temperature
        call update_soil_TI(d_h2oice,h2o_ice,ps_avg_glob,icetable_depth,icetable_thickness,ice_porefilling,icetable_equilibrium,icetable_dynamic,TI)

        ! Evolve adsorbed species
        if (do_sorption) then
            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)
            call compute_totmass_adsorbed(h2o_ads_reg,co2_ads_reg,totmass_adsco2,totmass_adsh2o)
        else
            totmass_adsh2o = 0._dp
            totmass_adsco2 = 0._dp
        end if
    end if ! do_soil
    deallocate(zshift_surf,zlag)

    ! Output the results
    call print_msg('> Standard outputs')
    call write_diagevol('ps_avg_glob','Global average pressure','Pa',ps_avg_glob)
    do islope = 1,nslope
        if (nslope /= 1) then
            num = '  '
            write(num,'(i2.2)') islope
            num = '_slope'//num
        end if
        call write_diagevol('h2oice'//num,'H2O ice','kg.m-2',h2o_ice(:,islope),(/dim_ngrid/))
        call write_diagevol('co2ice'//num,'CO2 ice','kg.m-2',co2_ice(:,islope),(/dim_ngrid/))
        call write_diagevol('d_h2oice'//num,'H2O ice tendency','kg.m-2.yr-1',d_h2oice(:,islope),(/dim_ngrid/))
        call write_diagevol('d_co2ice'//num,'CO2 ice tendency','kg.m-2.yr-1',d_co2ice(:,islope),(/dim_ngrid/))
        if (co2ice_flow) then
            call write_diagevol('Flow_co2ice'//num,'CO2 ice flow location','T/F',merge(1._dp,0._dp,is_co2ice_flow(:,islope)),(/dim_ngrid/))
            deallocate(is_co2ice_flow)
        end if
        if (h2oice_flow) then
            call write_diagevol('Flow_h2oice'//num,'H2O ice flow location','T/F',merge(1._dp,0._dp,is_h2oice_flow(:,islope)),(/dim_ngrid/))
            deallocate(is_h2oice_flow)
        end if
        call write_diagevol('tsurf'//num,'Surface temperature','K',tsurf_avg(:,islope),(/dim_ngrid/))
        if (do_soil) then
            if (icetable_equilibrium) then
                call write_diagevol('icetable_depth'//num,'Ice table depth','m',icetable_depth(:,islope),(/dim_ngrid/))
                call write_diagevol('icetable_thick'//num,'Ice table thickness','m',icetable_thickness(:,islope),(/dim_ngrid/))
            else if (icetable_dynamic) then
                call write_diagevol('icetable_depth'//num,'Ice table depth','m',icetable_depth(:,islope),(/dim_ngrid/))
                call write_diagevol('ice_porefilling'//num,'Ice pore filling','-',ice_porefilling(:,:,islope),(/dim_ngrid,dim_nsoil/))
            end if
            call write_diagevol('tsoil_avg'//num,'Soil temperature','K',tsoil_avg(:,:,islope),(/dim_ngrid,dim_nsoil/))
            call write_diagevol('inertiesoil'//num,'Thermal inertia','SI',TI(:,:,islope),(/dim_ngrid,dim_nsoil/))
            if (do_sorption) then
                call write_diagevol('co2_ads_reg'//num,'CO2 adsorbed in regolith','kg.m-2',co2_ads_reg(:,:,islope),(/dim_ngrid,dim_nsoil/))
                call write_diagevol('h2o_ads_reg'//num,'H2O adsorbed in regolith','kg.m-2',h2o_ads_reg(:,:,islope),(/dim_ngrid,dim_nsoil/))
            end if
        end if
    end do

    ! Checking mass balance for CO2
    if (abs(CO2cond_ps_PCM - 1._dp) < minieps) then
        totmass_co2ice = 0._dp
        totmass_atmco2 = 0._dp
        do i = 1,ngrid
            totmass_atmco2 = totmass_atmco2 + cell_area(i)*ps_avg(i)/g
            do islope = 1,nslope
                totmass_co2ice = totmass_co2ice + co2_ice(i,islope)*cell_area(i)*subslope_dist(i,islope)/cos(pi*def_slope_mean(islope)/180.)
            end do
        end do
        totmass_ini = max(totmass_atmco2_ini + totmass_co2ice_ini + totmass_adsco2_ini,minieps_qp) ! To avoid division by 0
        call print_msg(" > Relative total CO2 mass balance = "//real2str(100._qp*(totmass_atmco2 + totmass_co2ice + totmass_adsco2 - totmass_atmco2_ini - totmass_co2ice_ini - totmass_adsco2_ini)/totmass_ini)//' %')
        if (abs((totmass_atmco2 + totmass_co2ice + totmass_adsco2 - totmass_atmco2_ini - totmass_co2ice_ini - totmass_adsco2_ini)/totmass_ini) > 0.01_qp) then
            call print_msg('  /!\ Warning: mass balance is not conserved!')
            totmass_ini = max(totmass_atmco2_ini,minieps_qp) ! To avoid division by 0
            call print_msg('       Atmospheric CO2 mass balance = '//real2str(100._qp*(totmass_atmco2 - totmass_atmco2_ini)/totmass_ini)//' %')
            totmass_ini = max(totmass_co2ice_ini,minieps_qp) ! To avoid division by 0
            call print_msg('       CO2 ice mass balance         = '//real2str(100._qp*(totmass_co2ice - totmass_co2ice_ini)/totmass_ini)//' %')
            totmass_ini = max(totmass_adsco2_ini,minieps_qp) ! To avoid division by 0
            call print_msg('       Adsorbed CO2 mass balance    = '//real2str(100._qp*(totmass_adsco2 - totmass_adsco2_ini)/totmass_ini)//' %')
        end if
    end if

    ! Evolve the tendencies
    call evolve_tend_co2(d_co2ice_ini,co2_ice,emissivity_PCM,q_co2_ts_ini,q_co2_ts,ps_ts,ps_avg_glob_ini,ps_avg_glob,d_co2ice)
!~     call print_msg("> Updating the H2O sub-surface ice tendency due to lag layer")
!~     if (do_layering) then
!~         do i = 1,ngrid
!~             do islope = 1,nslope
!~                 if (is_h2oice_sublim_ini(i,islope) .and. h2oice_depth(i,islope) > 0._dp) call evolve_tend_h2o(h2oice_depth_old(i,islope),h2oice_depth(i,islope),tsurf_avg(i,islope),tsoil_ts_old(i,:,islope,:),tsoil_ts(i,:,islope,:),d_h2oice(i,islope))
!~             end do
!~         end do
!    else
!        do i = 1,ngrid
!            do islope = 1,nslope
!                call evolve_tend_h2o(icetable_depth_old(i,islope),icetable_depth(i,islope),tsurf_avg(i,islope),tsoil_ts_old(i,:,islope,:),tsoil_ts(i,:,islope,:),d_h2oice(i,islope))
!            end do
!        end do
!~     end if

    ! Increment time
    n_yr_run = n_yr_run + dt
    n_yr_sim = n_yr_sim + dt
    idt = idt + 1

    ! Check the stopping criteria
    call print_msg("> Checking the stopping criteria")
    call stopping_crit_pressure(ps_avg_glob_ini,ps_avg_glob,stopcrit)
    call stopping_crit_h2oice(h2oice_sublim_coverage_ini,h2o_ice,d_h2oice,stopcrit)
    call stopping_crit_co2ice(co2ice_sublim_coverage_ini,co2_ice,d_co2ice,stopcrit)
    call system_clock(c2)
    if (stopcrit%stop_code() == 0 .and. n_yr_run >= nmax_yr_run) stopcrit%nmax_yr_run_reached = .true.
    if (stopcrit%stop_code() == 0 .and. n_yr_run >= nmax_yr_runorb) stopcrit%nmax_yr_runorb_reached = .true.
    if (stopcrit%stop_code() == 0 .and. n_yr_sim >= nmax_yr_sim) stopcrit%nmax_yr_sim_reached = .true.
    if (stopcrit%stop_code() == 0 .and. timewall .and. real(c2 - c1,dp)/real(cr,dp) >= timelimit - antetime) stopcrit%time_limit_reached = .true.
    if (stopcrit%is_any_set()) then
        call print_msg(stopcrit%stop_message())
        exit
    else
        call print_msg('**** This run has achieved '//real2str(n_yr_run)//' Planetary years.')
        call print_msg('**** The chained simulation has achieved '//real2str(n_yr_sim)//' Planetary years.')
        call print_msg('**** This run is continuing!')
        call print_msg('****')
    end if
end do ! End of the evolution loop
deallocate(is_co2ice_disappeared)

! Finalization
! ~~~~~~~~~~~~
call print_msg('')
call print_msg('********* Finalization *********')
! Update orbital parameters
if (evol_orbit) call update_orbit(n_yr_sim,n_yr_run)

! Build ice for the PCM
allocate(h2o_ice4PCM(ngrid,nslope),co2_ice4PCM(ngrid,nslope),is_h2o_perice(ngrid))
call build4PCM_perice(h2o_ice,co2_ice,is_h2o_perice,h2o_ice4PCM,co2_ice4PCM)

! Build surface temperature for the PCM
allocate(tsurf4PCM(ngrid,nslope))
call build4PCM_tsurf(tsurf_avg,tsurf_dev,tsurf4PCM)

! Build soil for the PCM
if (do_soil) then
    allocate(tsoil4PCM(ngrid,nsoil_PCM,nslope),inertiesoil4PCM(ngrid,nsoil_PCM,nslope),flux_geo4PCM(ngrid,nslope))
    call build4PCM_soil(tsoil_avg,tsoil_dev,inertiesoil4PCM,tsoil4PCM,flux_geo4PCM)
end if

! Build atmosphere for the PCM
allocate(ps4PCM(ngrid),teta4PCM(ngrid,nlayer),air_mass4PCM(ngrid,nlayer))
call build4PCM_atmosphere(ps_avg,ps_dev,ps_avg_glob,ps_avg_glob_ini,ps4PCM,pa4PCM,preff4PCM,teta4PCM,air_mass4PCM)

! Build tracers for the PCM
allocate(q4PCM(ngrid,nlayer,nq))
call build4PCM_tracers(ps4PCM,q4PCM)

! Build surface radiative properties state for the PCM
allocate(albedo4PCM(ngrid,nslope),emissivity4PCM(ngrid,nslope))
call build4PCM_surf_rad_prop(h2o_ice,co2_ice,albedo4PCM,emissivity4PCM)

! Write the "startpem.nc"
call write_restartpem(h2o_ice,co2_ice,tsoil_avg,TI,icetable_depth,icetable_thickness,ice_porefilling,h2o_ads_reg,co2_ads_reg,layerings_map)

! Write the "startfi.nc"
call write_restartfi(is_h2o_perice,h2o_ice4PCM,co2_ice4PCM,tsurf4PCM,tsoil4PCM,inertiesoil4PCM,albedo4PCM,emissivity4PCM,flux_geo4PCM)

! Write the "start.nc"
call write_restart(ps4PCM,pa4PCM,preff4PCM,tsurf4PCM,q4PCM,teta4PCM,air_mass4PCM)

! Update the duration information of the simulation
call update_info(n_yr_run,stopcrit%stop_code(),n_yr_sim,nmax_yr_sim)

! Deallocation
deallocate(emissivity4PCM,albedo4PCM)
deallocate(q4PCM)
deallocate(ps4PCM,teta4PCM,air_mass4PCM)
if (do_soil) deallocate(tsoil4PCM,inertiesoil4PCM,flux_geo4PCM)
deallocate(tsurf4PCM)
deallocate(co2_ice4PCM,h2o_ice4PCM,is_h2o_perice)
deallocate(num)
if (do_layering) then
    deallocate(h2oice_depth,h2oice_depth_old,new_str,new_lag,current)
    do islope = 1,nslope
        do i = 1,ngrid
            call del_layering(layerings_map(i,islope))
        end do
    end do
end if
deallocate(layerings_map)
deallocate(h2o_ads_reg,co2_ads_reg)
deallocate(h2o_ice,co2_ice,is_h2oice_ini,is_co2ice_ini)
deallocate(ps_avg,ps_dev,ps_ts)
deallocate(tsurf_avg,tsurf_dev,h2o_surfdensity_avg)
deallocate(tsoil_avg,tsoil_dev,tsoil_ts,tsoil_ts_old,h2o_soildensity_avg)
deallocate(q_h2o_ts,q_co2_ts,q_co2_ts_ini)
deallocate(d_h2oice,d_co2ice,d_co2ice_ini)
deallocate(delta_h2o_ads,delta_co2_ads,delta_icetable,icetable_depth_old)
call end_allocation()

! Footer
call system_clock(c2)
call print_end(n_yr_run,real((c2 - c1),dp)/real(cr,dp),r_plnt2earth_yr,pem_ini_date,n_yr_sim)

END PROGRAM pem