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

implicit none

contains

subroutine pemdem0(filename,lonfi,latfi,cell_area,nsoil_PEM,ngrid,day_ini,time,nslope,def_slope,subslope_dist)
! create physics restart file and write time-independent variables
  use comsoil_h_PEM,     only: soil_pem,mlayer_PEM,fluxgeo,inertiedat_PEM
  use iostart_PEM,       only: open_restartphy, close_restartphy, put_var, put_field, length
  use mod_grid_phy_lmdz, only: klon_glo
  use time_phylmdz_mod,  only: daysec
#ifndef CPP_STD
  use planete_h,    only: aphelie, emin_turb, lmixmin, obliquit, peri_day, periheli, year_day
  use comcstfi_h,   only: g, mugaz, omeg, rad, rcp
#else
  use planete_mod,  only: apoastr, emin_turb, lmixmin, obliquit, peri_day, periastr, year_day
  use comcstfi_mod, only: g, mugaz, omeg, rad, rcp
#endif

  implicit none
 
  character(len=*), intent(in) :: filename
  real, intent(in)             :: lonfi(ngrid)
  real, intent(in)             :: latfi(ngrid)
  integer, intent(in)          :: nsoil_PEM
  integer, intent(in)          :: ngrid
  real, intent(in)             :: day_ini
  real, intent(in)             :: time
  real, intent(in)             :: cell_area(ngrid) !boundaries for bining of the slopes
  integer, intent(in)          :: nslope
  real, intent(in)             :: def_slope(nslope+1) !boundaries for bining of the slopes
  real, intent(in)             :: subslope_dist(ngrid,nslope) !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 (nogcm 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,water_reservoir)


  ! write time-dependent variable to restart file
  use iostart_PEM,   only: open_restartphy, close_restartphy, put_var, put_field
  use comsoil_h_PEM, only: mlayer_PEM,fluxgeo, inertiedat_PEM,soil_pem
  use time_evol_mod, only: year_bp_ini, convert_years
                
  implicit none
  
#ifndef CPP_STD
  include "callkeys.h"
#endif
  
  character(len=*), intent(in) :: filename
  integer, intent(in)          :: nsoil_PEM, ngrid, nslope, i_myear
  real, intent(in)             :: tsoil_slope_PEM(ngrid,nsoil_PEM,nslope) !under mesh bining according to slope
  real, intent(in)             :: inertiesoil_slope_PEM(ngrid,nsoil_PEM,nslope) !under mesh bining according to slope
  real, intent(in)             :: ice_table_depth(ngrid,nslope) !under mesh bining according to slope
  real, intent(in)             :: ice_table_thickness(ngrid,nslope) !under mesh bining according to slope
  real, intent(in)             :: m_co2_regolith(ngrid,nsoil_PEM,nslope)
  real, intent(in)             :: m_h2o_regolith(ngrid,nsoil_PEM,nslope)
  real, intent(in)             :: water_reservoir(ngrid)

  integer           :: islope
  character*2       :: num
  integer           :: iq
  character(len=30) :: txt  ! To store some text
  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('water_reservoir','water_reservoir',water_reservoir,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