module phyredem

implicit none

contains

subroutine physdem0(filename,lonfi,latfi,nsoil,ngrid,nlay,nq, &
                         phystep,day_ini,time,airefi, &
                         alb,ith,pzmea,pzstd,pzsig,pzgam,pzthe)
! create physics restart file and write time-independent variables
  use comsoil_h, only: volcapa, mlayer
  use geometry_mod, only: cell_area
  use surfdat_h, only: zmea, zstd, zsig, zgam, zthe, &
                       emisice, emissiv,             &
                       iceradius, dtemisice, phisfi
  use iostart, only : create_restartphy, close_restartphy, & 
                      put_var, put_field, length, ldscrpt, ndscrpt
  use mod_grid_phy_lmdz, only : klon_glo
  use planete_mod, only: year_day, periastr, apoastr, peri_day, &
                         obliquit, z0
  use comcstfi_mod, only: rad, omeg, g, mugaz, rcp
  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
  integer,intent(in) :: ngrid
  integer,intent(in) :: nlay
  integer,intent(in) :: nq
  real,intent(in) :: phystep
  real,intent(in) :: day_ini
  real,intent(in) :: time
  real,intent(in) :: airefi(ngrid)
  real,intent(in) :: alb(ngrid)
  real,intent(in) :: ith(ngrid,nsoil)
  real,intent(in) :: pzmea(ngrid)
  real,intent(in) :: pzstd(ngrid)
  real,intent(in) :: pzsig(ngrid)
  real,intent(in) :: pzgam(ngrid)
  real,intent(in) :: pzthe(ngrid)
 
  character(ndscrpt), dimension(ldscrpt), parameter :: dscrpt_tab_cntrl = (/ &
      "(1)  Number of atmospheric columns in physics     ", &
      "(2)  Number of atmospheric layers                 ", &
      "(3)  Final day                                    ", &
      "(4)  Final time of day                            ", &
      "(5)  Planet radius (m)                            ", &
      "(6)  Rotation rate (rad.s-1)                      ", &
      "(7)  Gravity (m.s-2)                              ", &
      "(8)  Molar mass of the atmosphere (g.mol-1)       ", &
      "(9)  = r/Cp           (=kappa in the dynamics)    ", &
      "(10) Length of a solar day (s)                    ", &
      "(11) Physics time step (s)                        ", &
      "(12) -                                            ", &
      "(13) -                                            ", &
      "(14) Length of year (in solar days)               ", &
      "(15) Minimum star-planet distance (AU)            ", &
      "(16) Maximum star-planet distance (AU)            ", &
      "(17) Date of periastro (sols since N. spring)     ", &
      "(18) Obliquity of the planet (deg)                ", &
      "(19) Default surface roughness (m)                ", &
      "(20) -                                            ", &
      "(21) -                                            ", &
      "(22) -                                            ", &
      "(23) -                                            ", &
      "(24) Emissivity of northern cap ~0.95             ", &
      "(25) Emissivity of southern cap ~0.95             ", &
      "(26) Emissivity of martian soil ~.95              ", &
      "(27) -                                            ", &
      "(28) -                                            ", &
      "(29) -                                            ", &
      "(30) -                                            ", &
      "(31) Mean scat radius of CO2 snow (north)         ", &
      "(32) Mean scat radius of CO2 snow (south)         ", &
      "(33) Time scale for snow metamorphism (north)     ", &
      "(34) Time scale for snow metamorphism (south)     ", &
      "(35) Soil volumetric heat capacity                "/)
  real :: tab_cntrl(length) ! nb "length=100" defined in iostart module
  
  ! Create physics start file
  call create_restartphy(filename)

! tab_cntrl() contains run parameters
  tab_cntrl(:)=0 ! initialization
!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
! Fill control array tab_cntrl(:) with paramleters for this run
!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
! Informations on the physics grid
  tab_cntrl(1) = float(klon_glo)  ! number of atmospheric columns in physics
  tab_cntrl(2) = float(nlay) ! number of atmospheric layers
  tab_cntrl(3) = day_ini + int(time)         ! final day 
  tab_cntrl(4) = time -int(time)            ! final time of day

! Informations about the planet, used by dynamics and physics
  tab_cntrl(5) = rad      ! planet radius (m) 
  tab_cntrl(6) = omeg     ! rotation rate (rad.s-1)
  tab_cntrl(7) = g        ! gravity (m.s-2) 
  tab_cntrl(8) = mugaz    ! Molar mass of the atmosphere (g.mol-1)
  tab_cntrl(9) = rcp      !  = r/cp  (=kappa in the dynamics)
  tab_cntrl(10) = daysec  ! length of a solar day (s) 

  tab_cntrl(11) = phystep  ! physics time step (s)
  tab_cntrl(12) = 0.
  tab_cntrl(13) = 0.

! Informations about the planet, only for physics
  tab_cntrl(14) = year_day  ! length of year (in solar days)
  tab_cntrl(15) = periastr  ! minimum star-planet distance (AU)
  tab_cntrl(16) = apoastr   ! maximum star-planet distance (AU) 
  tab_cntrl(17) = peri_day  ! date of periastron (sols since N. spring)
  tab_cntrl(18) = obliquit  ! Obliquity of the planet (deg)

! Boundary layer and turbulence
  tab_cntrl(19) = z0        ! surface roughness (m)
!  tab_cntrl(20) = lmixmin   ! mixing length (m)
!  tab_cntrl(21) = emin_turb ! minimal energy 

! Optical properties of polar caps and ground emissivity
  tab_cntrl(24) = emisice(1)   ! Emissivity of northern cap ~0.95
  tab_cntrl(25) = emisice(2)   ! Emissivity of southern cap ~0.95
  tab_cntrl(26) = emissiv      ! Emissivity of martian soil ~.95
  tab_cntrl(31) = iceradius(1) ! mean scat radius of CO2 snow (north)
  tab_cntrl(32) = iceradius(2) ! mean scat radius of CO2 snow (south)
  tab_cntrl(33) = dtemisice(1) ! time scale for snow metamorphism (north)
  tab_cntrl(34) = dtemisice(2) ! time scale for snow metamorphism (south)

  tab_cntrl(28) = 0. 
  tab_cntrl(29) = 0.
  tab_cntrl(30) = 0.
! Soil properties:
  tab_cntrl(35) = volcapa ! soil volumetric heat capacity

  ! Write the controle array
  call put_var("controle","Control parameters",tab_cntrl)
  
  ! Write the controle array descriptor
  call put_var("controle_descriptor",&
               "Description of control parameters",dscrpt_tab_cntrl)
  
  ! Write the mid-layer depths
  call put_var("soildepth","Soil mid-layer depth",mlayer)
  
  ! 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)
  
  ! Write surface geopotential
  call put_field("phisfi","Geopotential at the surface",phisfi)
  
  ! Write surface albedo
  call put_field("albedodat","Albedo of bare ground",alb)
  
  ! Subgrid topogaphy variables
  call put_field("ZMEA","Relief: mean relief",zmea)
  call put_field("ZSTD","Relief: standard deviation",zstd)
  call put_field("ZSIG","Relief: sigma parameter",zsig)
  call put_field("ZGAM","Relief: gamma parameter",zgam)
  call put_field("ZTHE","Relief: theta parameter",zthe)
  
  ! Soil thermal inertia
  call put_field("inertiedat","Soil thermal inertia",ith)
  
  ! Close file
  call close_restartphy
  
end subroutine physdem0

subroutine physdem1(filename,nsoil,ngrid,nlay,nq, &
                    phystep,time,tsurf,tsoil,emis,albedo,q2,qsurf, &
                    cloudfrac,totcloudfrac,hice, &
                    rnat,pctsrf_sic,tslab,tsea_ice,tice,sea_ice)
  ! write time-dependent variable to restart file
  use iostart, only : open_restartphy, close_restartphy, & 
                      put_var, put_field
  use tracer_h, only: noms
  USE radinc_h, ONLY: L_NSPECTV
!  use slab_ice_h, only: noceanmx
  USE ocean_slab_mod, ONLY: nslay
  use callkeys_mod, only: ok_slab_ocean, calllott_nonoro
  use nonoro_gwd_ran_mod, only: du_nonoro_gwd, dv_nonoro_gwd, &
                                east_gwstress, west_gwstress

  implicit none

  character(len=*),intent(in) :: filename
  integer,intent(in) :: nsoil
  integer,intent(in) :: ngrid
  integer,intent(in) :: nlay
  integer,intent(in) :: nq
  real,intent(in) :: phystep
  real,intent(in) :: time
  real,intent(in) :: tsurf(ngrid)
  real,intent(in) :: tsoil(ngrid,nsoil)
  real,intent(in) :: emis(ngrid)
  real,intent(in) :: albedo(ngrid,L_NSPECTV)
  real,intent(in) :: q2(ngrid,nlay+1)
  real,intent(in) :: qsurf(ngrid,nq)
  real,intent(in) :: cloudfrac(ngrid,nlay)
  real,intent(in) :: totcloudfrac(ngrid)
  real,intent(in) :: hice(ngrid)
  real,intent(in) :: rnat(ngrid)
  real,intent(in) :: pctsrf_sic(ngrid)
  real,intent(in) :: tslab(ngrid,nslay)
  real,intent(in) :: tsea_ice(ngrid)
  real,intent(in) :: sea_ice(ngrid)
  real,intent(in) :: tice(ngrid)

  integer :: iq
  
  ! 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","Temps de simulation",time)
  
  ! Surface temperature
  call put_field("tsurf","Surface temperature",tsurf)
  
  ! Soil temperature
  call put_field("tsoil","Soil temperature",tsoil)
  
  ! Emissivity of the surface
  call put_field("emis","Surface emissivity",emis)
  
  ! Surface albedo (assume homegeneous spectral albedo for now)
  call put_field("albedo","Surface albedo",albedo(:,1))
  
  ! Planetary Boundary Layer
  call put_field("q2","pbl wind variance",q2)
  
! cloud fraction and sea ice (NB: these should be optional... to be improved)
  call put_field("cloudfrac","Cloud fraction",cloudfrac)
  call put_field("totcloudfrac","Total fraction",totcloudfrac)
  call put_field("hice","Height of oceanic ice",hice)

 !Slab ocean 
 if(ok_slab_ocean) then
   call put_field("rnat","Nature of surface",rnat)
   call put_field("pctsrf_sic","Pourcentage sea ice",pctsrf_sic)
   call put_field("tslab","Temperature slab ocean",tslab)
   call put_field("tsea_ice","Temperature of top layer (seaice or snow)",tsea_ice)
   call put_field("tice","Temperature of seaice",tice)
   call put_field("sea_ice","Sea ice mass",sea_ice)
 endif!(ok_slab_ocean)

! tracers
  if (nq>0) then
    do iq=1,nq
      call put_field(noms(iq),"tracer on surface",qsurf(:,iq))
    enddo
  endif ! of if (nq>0)

! Non-orographic gavity waves
if (calllott_nonoro) then
   call put_field("du_nonoro_gwd","Zonal wind tendency due to GW",du_nonoro_gwd)
   call put_field("dv_nonoro_gwd","Meridional wind tendency due to GW",dv_nonoro_gwd)
   call put_field("east_gwstress","Eastward stress profile due to GW",east_gwstress)
   call put_field("west_gwstress","Westward stress profile due to GW",west_gwstress)
endif

! close file
      CALL close_restartphy
!$OMP BARRIER

end subroutine physdem1

end module phyredem