MODULE pemredem
!-----------------------------------------------------------------------
! NAME
!     pemredem
!
! DESCRIPTION
!     Write PEM-specific NetCDF restart files.
!
! AUTHORS & DATE
!     L. Lange
!     JB Clement, 2023-2025
!
! NOTES
!     Inspired by phyredem from the PCM. Handles time-independent and
!     time-dependent variables for restart functionality.
!-----------------------------------------------------------------------

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

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

!=======================================================================
SUBROUTINE pemdem0(filename,lonfi,latfi,cell_area,ngrid,nslope,def_slope,subslope_dist)
!-----------------------------------------------------------------------
! NAME
!     pemdem0
!
! DESCRIPTION
!     Create physics restart file and write time-independent variables.
!
! AUTHORS & DATE
!     L. Lange
!     JB Clement, 2023-2025
!
! NOTES
!
!-----------------------------------------------------------------------

! DEPENDENCIES
! ------------
use soil,        only: mlayer_PEM
use iostart_pem, only: open_restartphy, close_restartphy, put_var, put_field, length

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

! ARGUMENTS
! ---------
character(*),                  intent(in) :: filename
integer,                       intent(in) :: ngrid, nslope
real, dimension(ngrid),        intent(in) :: lonfi, latfi
real, dimension(ngrid),        intent(in) :: cell_area     ! boundaries for bining of the slopes
real, dimension(nslope + 1),   intent(in) :: def_slope     ! boundaries for bining of the slopes
real, dimension(ngrid,nslope), intent(in) :: subslope_dist ! undermesh statistics

! CODE
! ----
call open_restartphy(filename)

! Write the mid-layer depths
call put_var("soildepth","Soil mid-layer depth",mlayer_PEM)

! Write longitudes
call put_field("longitude","Longitudes of physics grid",lonfi)

! Write latitudes
call put_field("latitude","Latitudes of physics grid",latfi)

! Write mesh areas
call put_field("area","Mesh area",cell_area)

! Multidimensionnal variables (nopcm undermesh slope statistics)
call put_var("def_slope","slope criterium stages",def_slope)
call put_field("subslope_dist","under mesh slope distribution",subslope_dist)

! Close file
call close_restartphy

END SUBROUTINE pemdem0
!=======================================================================

!=======================================================================
SUBROUTINE pemdem1(filename,i_myear,nsoil_PEM,ngrid,nslope,tsoil_slope_PEM,inertiesoil_slope_PEM, &
                   icetable_depth,icetable_thickness,ice_porefilling,m_co2_regolith,m_h2o_regolith,h2o_ice,co2_ice,layerings_map)
!-----------------------------------------------------------------------
! NAME
!     pemdem1
!
! DESCRIPTION
!     Write time-dependent variables to restart file.
!
! AUTHORS & DATE
!     L. Lange
!     JB Clement, 2023-2025
!
! NOTES
!
!-----------------------------------------------------------------------

! DEPENDENCIES
! ------------
use iostart_pem,      only: open_restartphy, close_restartphy, put_var, put_field
use soil,             only: inertiedat_PEM, do_soil
use evolution,        only: year_bp_ini, convert_years
use layered_deposits, only: layering, nb_str_max, map2array, print_layering, layering_algo

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

! ARGUMENTS
! ---------
character(*),                            intent(in) :: filename
integer,                                 intent(in) :: nsoil_PEM, ngrid, nslope
real,                                    intent(in) :: i_myear
real, dimension(ngrid,nsoil_PEM,nslope), intent(in) :: tsoil_slope_PEM       ! under mesh bining according to slope
real, dimension(ngrid,nsoil_PEM,nslope), intent(in) :: inertiesoil_slope_PEM ! under mesh bining according to slope
real, dimension(ngrid,nslope),           intent(in) :: icetable_depth        ! under mesh bining according to slope
real, dimension(ngrid,nslope),           intent(in) :: icetable_thickness    ! under mesh bining according to slope
real, dimension(ngrid,nsoil_PEM,nslope), intent(in) :: ice_porefilling       ! under mesh bining according to slope
real, dimension(ngrid,nsoil_PEM,nslope), intent(in) :: m_co2_regolith, m_h2o_regolith
real, dimension(ngrid,nslope),           intent(in) :: h2o_ice, co2_ice
type(layering), dimension(ngrid,nslope), intent(in) :: layerings_map ! Layerings

! LOCAL VARIABLES
! ---------------
integer                               :: islope
character(2)                          :: num
real                                  :: Year          ! Year of the simulation
real, dimension(:,:,:,:), allocatable :: layerings_array ! Array for stratification (layerings)

! CODE
! ----
! Open file
call open_restartphy(filename)

! First variable to write must be "Time", in order to correctly
! set time counter in file
Year = (year_bp_ini + i_myear)*convert_years
call put_var("Time","Year of simulation",Year)
call put_field('h2o_ice','h2o_ice',h2o_ice,Year)
call put_field('co2_ice','co2_ice',co2_ice,Year)

if (layering_algo) then
    allocate(layerings_array(ngrid,nslope,nb_str_max,6))
    call map2array(layerings_map,ngrid,nslope,layerings_array)
    do islope = 1,nslope
        write(num,fmt='(i2.2)') islope
        call put_field('stratif_slope'//num//'_top_elevation','Layering top elevation',layerings_array(:,islope,:,1),Year)
        call put_field('stratif_slope'//num//'_h_co2ice','Layering CO2 ice height',layerings_array(:,islope,:,2),Year)
        call put_field('stratif_slope'//num//'_h_h2oice','Layering H2O ice height',layerings_array(:,islope,:,3),Year)
        call put_field('stratif_slope'//num//'_h_dust','Layering dust height',layerings_array(:,islope,:,4),Year)
        call put_field('stratif_slope'//num//'_h_pore','Layering pore height',layerings_array(:,islope,:,5),Year)
        call put_field('stratif_slope'//num//'_poreice_volfrac','Layering ice pore volume fraction',layerings_array(:,islope,:,6),Year)
    enddo
    deallocate(layerings_array)
endif

if (do_soil) then
  ! Multidimensionnal variables (undermesh slope statistics)
    do islope = 1,nslope
        write(num,fmt='(i2.2)') islope
        call put_field("tsoil_PEM_slope"//num,"Soil temperature by slope type",tsoil_slope_PEM(:,:,islope),Year)
        call put_field("TI_PEM_slope"//num,"Soil Thermal Inertia by slope type",inertiesoil_slope_PEM(:,:,islope),Year)
        call put_field("mco2_reg_ads_slope"//num, "Mass of CO2 adsorbed in the regolith",m_co2_regolith(:,:,islope),Year)
        call put_field("mh2o_reg_ads_slope"//num, "Mass of H2O adsorbed in the regolith",m_h2o_regolith(:,:,islope),Year)
        call put_field("ice_porefilling"//num,"Subsurface ice pore filling",ice_porefilling(:,:,islope),Year)
    enddo
    call put_field("icetable_depth","Depth of ice table",icetable_depth,Year)
    call put_field("icetable_thickness","Depth of ice table",icetable_thickness,Year)
    call put_field("inertiedat_PEM","Thermal inertie of PEM ",inertiedat_PEM,Year)
endif ! do_soil

! Close file
call close_restartphy

END SUBROUTINE pemdem1

END MODULE pemredem