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(ngrid + 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,m_co2_regolith,m_h2o_regolith,h2o_ice)

! 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

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) :: m_co2_regolith, m_h2o_regolith
real, dimension(ngrid,nslope),           intent(in) :: h2o_ice

integer       :: islope
character(2)  :: num
real          :: Year ! Year of the simulation

! 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 (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("inertiedat_PEM","Thermal inertie of PEM ",inertiedat_PEM,Year)
endif ! soil_pem

! Close file
call close_restartphy

END SUBROUTINE pemdem1

END MODULE pemredem