MODULE pemredem

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!
!!! Purpose: Write specific netcdf restart for the PEM
!!!
!!!
!!! Author: LL, inspired by phyredem from the PCM
!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

implicit none

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

SUBROUTINE pemdem0(filename,lonfi,latfi,cell_area,ngrid,nslope,def_slope,subslope_dist)

! create physics restart file and write time-independent variables
use comsoil_h_PEM, only: mlayer_PEM
use iostart_PEM,   only: open_restartphy, close_restartphy, put_var, put_field, length

implicit none

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

! Create physics start file
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, &
                   ice_table_depth,ice_table_thickness,ice_porefilling,m_co2_regolith,m_h2o_regolith,h2o_ice,stratif)

! write time-dependent variable to restart file
use iostart_PEM,   only: open_restartphy, close_restartphy, put_var, put_field
use comsoil_h_PEM, only: inertiedat_PEM, soil_pem
use time_evol_mod, only: year_bp_ini, convert_years
use layering_mod,  only: layering, nb_str_max, stratif2array, print_layering, layering_algo

implicit none

#ifndef CPP_STD
    include "callkeys.h"
#endif

character(*),                            intent(in) :: filename
integer,                                 intent(in) :: nsoil_PEM, ngrid, nslope, 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) :: ice_table_depth       ! under mesh bining according to slope
real, dimension(ngrid,nslope),           intent(in) :: ice_table_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
type(layering), dimension(ngrid,nslope), intent(in) :: stratif ! Stratification (layerings)

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

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

if (layering_algo) then
    allocate(stratif_array(ngrid,nslope,nb_str_max,6))
    call stratif2array(stratif,ngrid,nslope,stratif_array)
    do islope = 1,nslope
        write(num,fmt='(i2.2)') islope
        call put_field('stratif_slope'//num//'_thickness','Layering thickness',stratif_array(:,islope,:,1),Year)
        call put_field('stratif_slope'//num//'_top_elevation','Layering top elevation',stratif_array(:,islope,:,2),Year)
        call put_field('stratif_slope'//num//'_co2ice_volfrac','Layering CO2 ice volume fraction',stratif_array(:,islope,:,3),Year)
        call put_field('stratif_slope'//num//'_h2oice_volfrac','Layering H2O ice volume fraction',stratif_array(:,islope,:,4),Year)
        call put_field('stratif_slope'//num//'_dust_volfrac','Layering dust volume fraction',stratif_array(:,islope,:,5),Year)
        call put_field('stratif_slope'//num//'_air_volfrac','Layering air volume fraction',stratif_array(:,islope,:,6),Year)
    enddo
    deallocate(stratif_array)
endif

if (soil_pem) 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 adsorbded in the regolith",m_co2_regolith(:,:,islope),Year)
        call put_field("mh2o_reg_ads_slope"//num, "Mass of h2o adsorbded in the regolith",m_h2o_regolith(:,:,islope),Year)
    enddo
    call put_field("ice_table_depth","Depth of ice table",ice_table_depth,Year)
    call put_field("ice_table_thickness","Depth of ice table",ice_table_thickness,Year)
    call put_field("ice_porefilling","Subsurface ice pore filling",ice_porefilling,Year)
    call put_field("inertiedat_PEM","Thermal inertie of PEM ",inertiedat_PEM,Year)
endif ! soil_pem

! Close file
call close_restartphy

END SUBROUTINE pemdem1

END MODULE pemredem