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
#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 comconst_mod, ONLY: g, omeg, rad
  use comcstfi_mod, only: mugaz, rcp
#endif
  use time_phylmdz_mod, only: daysec
  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)

  call put_var("FluxGeo","Geothermal Flux (W/m^2)",fluxgeo)

  ! Close file
  call close_restartphy
  
end subroutine pemdem0

subroutine pemdem1(filename,year_iter,nsoil_PEM,ngrid,nslope, &
                    tsoil_slope_PEM,inertiesoil_slope_PEM,ice_table, &
                    m_co2_regolith,m_h2o_regolith)
  ! 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 temps_mod_evol, ONLY: year_PEM
                
  implicit none
  
#ifndef CPP_STD
  include "callkeys.h"
#endif
  
  character(len=*),intent(in) :: filename
  integer,intent(in) :: nsoil_PEM
  integer,intent(in) :: ngrid
  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(ngrid,nslope) !under mesh bining according to slope
  integer,intent(IN) :: nslope
  real,intent(in) :: m_co2_regolith(ngrid,nsoil_PEM,nslope)
  real,intent(in) :: m_h2o_regolith(ngrid,nsoil_PEM,nslope)
  integer :: islope
  CHARACTER*2 :: num  
  integer, intent(IN) :: year_iter

  real :: year_tot
 
  integer :: iq
  character(len=30) :: txt ! to store some text
  ! indexes of water vapour & water ice tracers (if any):

  year_tot=real(year_iter+year_PEM)
  print *, "Year of simulation=", year_tot  

  ! Open file
  call open_restartphy(filename)
  
  ! First variable to write must be "Time", in order to correctly
  ! set time counter in file
  call put_var("Time","Year of simulation",year_tot)

    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_tot)
  call put_field("TI_PEM_slope"//num,"Soil Thermal Inertia by slope type", &
                 inertiesoil_slope_PEM(:,:,islope),year_tot)

 call put_field("mco2_reg_ads_slope"//num, "Mass of co2 adsorbded in the regolith", &
                    m_co2_regolith(:,:,islope), year_tot)
  call put_field("mh2o_reg_ads_slope"//num, "Mass of h2o adsorbded in the regolith", &
                    m_h2o_regolith(:,:,islope), year_tot)

ENDDO

  call put_field("ice_table","Depth of ice table", &
                 ice_table,year_tot)

  call put_field("inertiedat_PEM","Thermal inertie of PEM ", &
                 inertiedat_PEM,year_tot)

   endif !soil_pem

  ! Close file
  call close_restartphy
  
end subroutine pemdem1

end module pemredem