MODULE planet
!-----------------------------------------------------------------------
! NAME
!     planet
!
! DESCRIPTION
!     Centralized climate-state storage and lifecycle helpers for PEM.
!
! AUTHORS & DATE
!     JB Clement, 03/2026
!
! NOTES
!     Ownership criterion for declarations:
!       - Declare in module "planet" variables that are persistent climate
!         state or persistent references used across multiple time steps.
!       - Declare in program "pem.F90" transient workflow/control varaibles,
!         one-shot initialization buffers and per-iteration working buffers.
!-----------------------------------------------------------------------

! DEPENDENCIES
! ------------
use numerics,         only: dp, qp, k4
use layered_deposits, only: layering, del_layering

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

! PARAMETERS
! ----------
! Pressure-related:
real(dp), dimension(:),   allocatable :: ps_avg          ! Average surface pressure [Pa]
real(dp), dimension(:,:), allocatable :: ps_ts           ! Surface pressure timeseries [Pa]
real(dp), dimension(:),   allocatable :: ps_dev          ! Deviation of surface pressure [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]

! Ice-related:
real(dp),    dimension(:,:), allocatable :: h2o_ice                    ! H2O ice [kg/m2]
real(dp),    dimension(:,:), allocatable :: co2_ice                    ! CO2 ice [kg/m2]
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_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_dev           ! Deviation of surface temperature [K]
real(dp), dimension(:,:), allocatable :: h2o_surfdensity_avg ! Average water surface density [kg/m^3]

! Soil-related:
real(dp), dimension(:,:,:),   allocatable :: tsoil_avg    ! Average soil temperature [K]
real(dp), dimension(:,:,:),   allocatable :: tsoil_dev    ! Deviation of 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]

! Layering-related:
type(layering), dimension(:,:), allocatable :: layerings_map ! Layering for each grid point and slope

! Sorption-related:
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(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]

! Ice table-related:
real(dp), dimension(:,:),   allocatable :: icetable_depth     ! Depth of the ice table [m]
real(dp), dimension(:,:),   allocatable :: icetable_thickness ! Thickness of the ice table [m]
real(dp), dimension(:,:,:), allocatable :: ice_porefilling    ! Amount of porefilling [m^3/m^3]
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 exchanged with the ice table [kg]
real(dp), dimension(:,:),   allocatable :: flux_ssice_avg     ! Average of total flux exchanged with subsurface ice [kg/m2/y]

! 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]

! 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]

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

!=======================================================================
SUBROUTINE ini_planet()
!-----------------------------------------------------------------------
! NAME
!     ini_planet
!
! DESCRIPTION
!     Initialize scalar state values of module planet.
!
! AUTHORS & DATE
!     JB Clement, 03/2026
!
! NOTES
!
!-----------------------------------------------------------------------

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

! CODE
! ----
ps_avg_glob_ini = 0._dp
ps_avg_glob_old = 0._dp
ps_avg_glob = 0._dp
h2oice_sublim_coverage_ini = 0._dp
co2ice_sublim_coverage_ini = 0._dp

END SUBROUTINE ini_planet
!=======================================================================

!=======================================================================
SUBROUTINE allocate_xios_state()
!-----------------------------------------------------------------------
! NAME
!     allocate_xios_state
!
! DESCRIPTION
!     Allocate arrays needed before loading PCM/XIOS data.
!
! AUTHORS & DATE
!     JB Clement, 03/2026
!
! NOTES
!
!-----------------------------------------------------------------------

! DEPENDENCIES
! ------------
use geometry, only: ngrid, nslope, nday, nsoil

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

! CODE
! ----
allocate(ps_avg(ngrid))
allocate(ps_ts(ngrid,nday))
allocate(tsurf_avg(ngrid,nslope))
allocate(h2o_surfdensity_avg(ngrid,nslope))
allocate(tsoil_avg(ngrid,nsoil,nslope))
allocate(tsoil_ts(ngrid,nsoil,nslope,nday))
allocate(h2o_soildensity_avg(ngrid,nsoil,nslope))
allocate(q_h2o_ts(ngrid,nday))
allocate(q_co2_ts(ngrid,nday))
allocate(flux_ssice_avg(ngrid,nslope))

ps_avg(:) = 0._dp
ps_ts(:,:) = 0._dp
tsurf_avg(:,:) = 0._dp
h2o_surfdensity_avg(:,:) = 0._dp
tsoil_avg(:,:,:) = 0._dp
tsoil_ts(:,:,:,:) = 0._dp
h2o_soildensity_avg(:,:,:) = 0._dp
q_h2o_ts(:,:) = 0._dp
q_co2_ts(:,:) = 0._dp
flux_ssice_avg(:,:) = 0._dp

END SUBROUTINE allocate_xios_state
!=======================================================================

!=======================================================================
SUBROUTINE allocate_deviation_state()
!-----------------------------------------------------------------------
! NAME
!     allocate_deviation_state
!
! DESCRIPTION
!     Allocate persistent deviation fields derived from PCM averages.
!
! AUTHORS & DATE
!     JB Clement, 03/2026
!
! NOTES
!
!-----------------------------------------------------------------------

! DEPENDENCIES
! ------------
use geometry, only: ngrid, nslope, nsoil_PCM

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

! CODE
! ----
allocate(ps_dev(ngrid))
allocate(tsurf_dev(ngrid,nslope))
allocate(tsoil_dev(ngrid,nsoil_PCM,nslope))

ps_dev(:) = 0._dp
tsurf_dev(:,:) = 0._dp
tsoil_dev(:,:,:) = 0._dp

END SUBROUTINE allocate_deviation_state
!=======================================================================

!=======================================================================
SUBROUTINE allocate_startevo_state()
!-----------------------------------------------------------------------
! NAME
!     allocate_startevo_state
!
! DESCRIPTION
!     Allocate arrays loaded from startevo or evolving afterward.
!
! AUTHORS & DATE
!     JB Clement, 03/2026
!
! NOTES
!
!-----------------------------------------------------------------------

! DEPENDENCIES
! ------------
use geometry, only: ngrid, nslope, nsoil

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

! LOCAL VARIABLES
! ---------------
integer :: i, islope

! CODE
! ----
allocate(h2o_ice(ngrid,nslope))
allocate(co2_ice(ngrid,nslope))
allocate(icetable_depth(ngrid,nslope))
allocate(icetable_thickness(ngrid,nslope))
allocate(ice_porefilling(ngrid,nsoil,nslope))
allocate(h2o_ads_reg(ngrid,nsoil,nslope))
allocate(co2_ads_reg(ngrid,nsoil,nslope))
allocate(delta_h2o_ads(ngrid))
allocate(delta_co2_ads(ngrid))
allocate(layerings_map(ngrid,nslope))

h2o_ice(:,:) = 0._dp
co2_ice(:,:) = 0._dp
icetable_depth(:,:) = 0._dp
icetable_thickness(:,:) = 0._dp
ice_porefilling(:,:,:) = 0._dp
h2o_ads_reg(:,:,:) = 0._dp
co2_ads_reg(:,:,:) = 0._dp
delta_h2o_ads(:) = 0._dp
delta_co2_ads(:) = 0._dp
do islope = 1,nslope
    do i = 1,ngrid
        layerings_map(i,islope)%nb_str = 0
        nullify(layerings_map(i,islope)%bottom,layerings_map(i,islope)%top)
    end do
end do

END SUBROUTINE allocate_startevo_state
!=======================================================================

!=======================================================================
SUBROUTINE allocate_tendencies()
!-----------------------------------------------------------------------
! NAME
!     allocate_tendencies
!
! DESCRIPTION
!     Allocate perennial ice tendency arrays.
!
! AUTHORS & DATE
!     JB Clement, 03/2026
!
! NOTES
!
!-----------------------------------------------------------------------

! DEPENDENCIES
! ------------
use geometry, only: ngrid, nslope

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

! CODE
! ----
allocate(d_h2oice(ngrid,nslope))
allocate(d_co2ice(ngrid,nslope))

d_h2oice(:,:) = 0._dp
d_co2ice(:,:) = 0._dp

END SUBROUTINE allocate_tendencies
!=======================================================================

!=======================================================================
SUBROUTINE allocate_initial_snapshots()
!-----------------------------------------------------------------------
! NAME
!     allocate_initial_snapshots
!
! DESCRIPTION
!     Allocate arrays storing initial-state snapshots.
!
! AUTHORS & DATE
!     JB Clement, 03/2026
!
! NOTES
!
!-----------------------------------------------------------------------

! DEPENDENCIES
! ------------
use geometry, only: ngrid, nslope, nday

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

! CODE
! ----
allocate(d_co2ice_ini(ngrid,nslope))
allocate(q_co2_ts_ini(ngrid,nday))
allocate(is_h2oice_ini(ngrid,nslope))
allocate(is_co2ice_ini(ngrid,nslope))

d_co2ice_ini(:,:) = 0._dp
q_co2_ts_ini(:,:) = 0._dp
is_h2oice_ini(:,:) = .false.
is_co2ice_ini(:,:) = .false.

END SUBROUTINE allocate_initial_snapshots
!=======================================================================

!=======================================================================
SUBROUTINE allocate_loop_state()
!-----------------------------------------------------------------------
! NAME
!     allocate_loop_state
!
! DESCRIPTION
!     Allocate arrays that start being used during the main PEM loop.
!
! AUTHORS & DATE
!     JB Clement, 03/2026
!
! NOTES
!
!-----------------------------------------------------------------------

! DEPENDENCIES
! ------------
use geometry, only: ngrid, nslope, nsoil, nday

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

! CODE
! ----
allocate(delta_icetable(ngrid))
allocate(icetable_depth_old(ngrid,nslope))
allocate(is_co2ice_disappeared(ngrid,nslope))
allocate(tsoil_ts_old(ngrid,nsoil,nslope,nday))

delta_icetable(:) = 0._dp
icetable_depth_old(:,:) = 0._dp
is_co2ice_disappeared(:,:) = .false.
tsoil_ts_old(:,:,:,:) = 0._dp

END SUBROUTINE allocate_loop_state
!=======================================================================

!=======================================================================
SUBROUTINE end_planet()
!-----------------------------------------------------------------------
! NAME
!     end_planet
!
! DESCRIPTION
!     Finalize module planet.
!
! AUTHORS & DATE
!     JB Clement, 03/2026
!
! NOTES
!     Deallocate all allocatable state arrays of module planet.
!-----------------------------------------------------------------------

! DEPENDENCIES
! ------------
use geometry, only: ngrid, nslope

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

! LOCAL VARIABLES
! ---------------
integer :: i, islope

! CODE
! ----
if (allocated(ps_avg)) deallocate(ps_avg)
if (allocated(ps_ts)) deallocate(ps_ts)
if (allocated(ps_dev)) deallocate(ps_dev)
if (allocated(h2o_ice)) deallocate(h2o_ice)
if (allocated(co2_ice)) deallocate(co2_ice)
if (allocated(is_h2oice_ini)) deallocate(is_h2oice_ini)
if (allocated(is_co2ice_ini)) deallocate(is_co2ice_ini)
if (allocated(is_co2ice_disappeared)) deallocate(is_co2ice_disappeared)
if (allocated(tsurf_avg)) deallocate(tsurf_avg)
if (allocated(tsurf_dev)) deallocate(tsurf_dev)
if (allocated(h2o_surfdensity_avg)) deallocate(h2o_surfdensity_avg)
if (allocated(tsoil_avg)) deallocate(tsoil_avg)
if (allocated(tsoil_dev)) deallocate(tsoil_dev)
if (allocated(tsoil_ts)) deallocate(tsoil_ts)
if (allocated(tsoil_ts_old)) deallocate(tsoil_ts_old)
if (allocated(layerings_map)) then
    do islope = 1,nslope
        do i = 1,ngrid
            call del_layering(layerings_map(i,islope))
        end do
    end do
    deallocate(layerings_map)
end if
if (allocated(h2o_soildensity_avg)) deallocate(h2o_soildensity_avg)
if (allocated(delta_co2_ads)) deallocate(delta_co2_ads)
if (allocated(delta_h2o_ads)) deallocate(delta_h2o_ads)
if (allocated(h2o_ads_reg)) deallocate(h2o_ads_reg)
if (allocated(co2_ads_reg)) deallocate(co2_ads_reg)
if (allocated(icetable_depth)) deallocate(icetable_depth)
if (allocated(icetable_thickness)) deallocate(icetable_thickness)
if (allocated(ice_porefilling)) deallocate(ice_porefilling)
if (allocated(icetable_depth_old)) deallocate(icetable_depth_old)
if (allocated(delta_icetable)) deallocate(delta_icetable)
if (allocated(q_co2_ts)) deallocate(q_co2_ts)
if (allocated(q_co2_ts_ini)) deallocate(q_co2_ts_ini)
if (allocated(q_h2o_ts)) deallocate(q_h2o_ts)
if (allocated(d_co2ice)) deallocate(d_co2ice)
if (allocated(d_co2ice_ini)) deallocate(d_co2ice_ini)
if (allocated(d_h2oice)) deallocate(d_h2oice)
if (allocated(flux_ssice_avg)) deallocate(flux_ssice_avg)

END SUBROUTINE end_planet
!=======================================================================

END MODULE planet