source: trunk/LMDZ.MARS/libf/phymars/phyredem.F90 @ 3026

module phyredem

implicit none

contains

subroutine physdem0(filename,lonfi,latfi,nsoil,ngrid,nlay,nq, &
                         phystep,day_ini,time,airefi, &
                         alb,ith,def_slope, &
                         subslope_dist)
! create physics restart file and write time-independent variables
  use comsoil_h, only: inertiedat, volcapa, mlayer
  use geometry_mod, only: cell_area
  use surfdat_h, only: zmea, zstd, zsig, zgam, zthe, &
                       z0_default, albedice, emisice, emissiv, &
                       iceradius, dtemisice, phisfi, z0,   &
                       hmons,summit,base,watercaptag
  use dimradmars_mod, only: tauvis
  use iostart, only : open_restartphy, close_restartphy, & 
                      put_var, put_field, length
  use mod_grid_phy_lmdz, only : klon_glo
  use planete_h, only: aphelie, emin_turb, lmixmin, obliquit, &
                       peri_day, periheli, year_day
  use comcstfi_h, only: g, mugaz, omeg, rad, rcp
  use time_phylmdz_mod, only: daysec
  use comslope_mod, ONLY: nslope
  USE paleoclimate_mod, ONLY: paleoclimate, lag_h2o_ice, lag_co2_ice  
  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) ! thermal inertia for present day 
  real, intent(in) :: def_slope(nslope+1) !boundaries for bining of the slopes
  real, intent(in) :: subslope_dist(ngrid,nslope) !undermesh statistics

  real :: tab_cntrl(length) ! nb "length=100" defined in iostart module
  integer :: ig
  real :: watercaptag_tmp(ngrid)

  ! Create physics start file
  call open_restartphy(filename)
  
  ! Build tab_cntrl(:) array
  tab_cntrl(:)=0.0
  !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
  ! Fill control array tab_cntrl(:) with parameters for this run
  !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
  ! Informations on the physics grid
  tab_cntrl(1) = float(klon_glo)  ! total number of nodes on physics grid
  tab_cntrl(2) = float(nlay) ! number of atmospheric layers
  tab_cntrl(3) = day_ini + int(time)      ! initial day 
  tab_cntrl(4) = time -int(time)          ! initial time of day

  ! Informations about Mars, used by dynamics and physics
  tab_cntrl(5) = rad      ! radius of Mars (m) ~3397200
  tab_cntrl(6) = omeg     ! rotation rate (rad.s-1)
  tab_cntrl(7) = g        ! gravity (m.s-2) ~3.72
  tab_cntrl(8) = mugaz    ! Molar mass of the atmosphere (g.mol-1) ~43.49
  tab_cntrl(9) = rcp      !  = r/cp  ~0.256793 (=kappa dans dynamique)
  tab_cntrl(10) = daysec  ! length of a sol (s)  ~88775

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

  ! Informations about Mars, only for physics
  tab_cntrl(14) = year_day  ! length of year (sols) ~668.6
  tab_cntrl(15) = periheli  ! min. Sun-Mars distance (Mkm) ~206.66
  tab_cntrl(16) = aphelie   ! max. SUn-Mars distance (Mkm) ~249.22
  tab_cntrl(17) = peri_day  ! date of perihelion (sols since N. spring)
  tab_cntrl(18) = obliquit  ! Obliquity of the planet (deg) ~23.98

  ! Boundary layer and turbulence
  tab_cntrl(19) = z0_default   ! default surface roughness (m) ~0.01
  tab_cntrl(20) = lmixmin   ! mixing length ~100
  tab_cntrl(21) = emin_turb ! minimal energy ~1.e-8

  ! Optical properties of polar caps and ground emissivity
  tab_cntrl(22) = albedice(1)  ! Albedo of northern cap ~0.5
  tab_cntrl(23) = albedice(2)  ! Albedo of southern cap ~0.5
  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)

  ! dust aerosol properties
  tab_cntrl(27) = tauvis      ! mean visible optical depth

  ! Soil properties:
  tab_cntrl(35) = volcapa ! soil volumetric heat capacity

  ! Write the controle array
  call put_var("controle","Control parameters",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)
  call put_field("hmons","Relief: hmons parameter (summit - base)",hmons)
  call put_field("summit","Relief: altitude from the aeroid of the summit of the highest subgrid topography",summit)
  call put_field("base","Relief: altitude from the aeroid of the base of the highest subgrid topography",base)
     
  ! Soil thermal inertia
  call put_field("inertiedat","Soil thermal inertia - present day",ith)
  
  ! Surface roughness length
  call put_field("z0","Surface roughness length",z0)
  
  do ig=1,ngrid
    if(watercaptag(ig)) then
       watercaptag_tmp(ig)=1
    else
       watercaptag_tmp(ig)=0
    endif
  enddo

  call put_field("watercaptag","Infinite water reservoir",watercaptag_tmp)

  ! Sub grid scale slope parametrization
  call put_var("def_slope","slope criterium stages",def_slope)
  call put_field("subslope_dist","under mesh slope distribution",subslope_dist)

  ! Paleoclimate outputs
  if(paleoclimate) then
    call put_field("lag_h2o_ice","Depth of the H2O lags",lag_h2o_ice)
    call put_field("lag_co2_ice","Depth of the CO2 lags",lag_co2_ice)
  endif
  ! Close file
  call close_restartphy
  
end subroutine physdem0

subroutine physdem1(filename,nsoil,ngrid,nlay,nq, &
                    phystep,time,tsurf,tsoil,inertiesoil, &
                    albedo,emis,q2,qsurf,&
                    tauscaling,totcloudfrac,wstar, &
                    watercap,perenial_co2ice)
  ! write time-dependent variable to restart file
  use iostart, only : open_restartphy, close_restartphy, & 
                      put_var, put_field
  use tracer_mod, only: noms ! tracer names
  use nonoro_gwd_ran_mod, only: du_nonoro_gwd, dv_nonoro_gwd
  use compute_dtau_mod, only: dtau
  use dust_rad_adjust_mod, only: dust_rad_adjust_prev,dust_rad_adjust_next
  use dust_param_mod, only: dustscaling_mode
  use comsoil_h,only: flux_geo
  use comslope_mod, only: nslope
  use paleoclimate_mod, only: paleoclimate
  implicit none
  
  include "callkeys.h"
  
  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,nslope)
  real,intent(in) :: tsoil(ngrid,nsoil,nslope)
  real,intent(in) :: inertiesoil(ngrid,nsoil,nslope)
  real,intent(in) :: albedo(ngrid,2,nslope)
  real,intent(in) :: emis(ngrid,nslope)
  real,intent(in) :: q2(ngrid,nlay+1)
  real,intent(in) :: qsurf(ngrid,nq,nslope)
  real,intent(in) :: tauscaling(ngrid)
  real,intent(in) :: totcloudfrac(ngrid)
  real,intent(in) :: wstar(ngrid)
  real,intent(in) :: watercap(ngrid,nslope)
  real,intent(in) :: perenial_co2ice(ngrid,nslope)

  integer :: iq
  character(len=30) :: txt ! to store some text
  ! indexes of water vapour & water ice tracers (if any):
  integer :: i_h2o_vap=0
  integer :: i_h2o_ice=0
  integer :: i_hdo_vap=0
  integer :: i_hdo_ice=0

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

  ! Water ice layer
  call put_field("watercap","H2O ice cover",watercap,time)
 
 ! Perenial CO2 ice layer
  if(paleoclimate) then
    call put_field("perenial_co2ice","CO2 ice cover",perenial_co2ice,time)
  endif
  ! Surface temperature
  call put_field("tsurf","Surface temperature",tsurf,time)

  ! Soil temperature
  call put_field("inertiesoil","Soil thermal inertia",inertiesoil,time)
  
  ! Soil temperature
  call put_field("tsoil","Soil temperature",tsoil,time)
  
  ! Albedo of the surface
  call put_field("albedo","Surface albedo",albedo(:,1,:),time)
  
  ! Emissivity of the surface
  call put_field("emis","Surface emissivity",emis,time)
  
  ! Planetary Boundary Layer
  call put_field("q2","pbl wind variance",q2,time)

  ! Sub-grid cloud fraction
  call put_field("totcloudfrac","Total cloud fraction",totcloudfrac,time) 
 
  ! Dust conversion factor
  ! Only to be read by newstart to convert to actual dust values
  ! Or by any user who wants to reconstruct dust, opacity from the start files.
  call put_field("tauscaling","dust conversion factor",tauscaling,time)

  ! Radiative scaling coefficients
  if (dustscaling_mode==2) then
    call put_field("dust_rad_adjust_prev", &
                   "radiative dust adjustement factor prev. sol", &
                   dust_rad_adjust_prev,time)
    call put_field("dust_rad_adjust_next", &
                   "radiative dust adjustement factor next sol", &
                   dust_rad_adjust_next,time)
  endif
  
  if (dustinjection.gt.0) then
    call put_field("dtau","dust opacity difference between GCM and scenario",&
                   dtau,time)
  endif

  if (calltherm) then
    call put_field("wstar","Max vertical velocity in thermals",wstar,time)
  endif

  ! Tracers on the surface
  ! preliminary stuff: look for water vapour & water ice tracers (if any)
  do iq=1,nq
    if (noms(iq).eq."h2o_vap") then
      i_h2o_vap=iq
    endif
    if (noms(iq).eq."h2o_ice") then
      i_h2o_ice=iq
    endif

  ! look for HDO vapour & HDO ice tracers (if any)
    if (noms(iq).eq."hdo_vap") then
      i_hdo_vap=iq
    endif
    if (noms(iq).eq."hdo_ice") then
      i_hdo_ice=iq
    endif
  enddo

  
  if (nq.gt.0) then
    do iq=1,nq
      txt=noms(iq)
      ! Exception: there is no water vapour surface tracer
      if (txt.eq."h2o_vap") then
        write(*,*)"physdem1: skipping water vapour tracer"
        if (i_h2o_ice.eq.i_h2o_vap) then
          ! then there is no "water ice" tracer; but still
          ! there is some water ice on the surface
          write(*,*)"          writing water ice instead"
          txt="h2o_ice"
        else
          ! there is a "water ice" tracer which has been / will be
          ! delt with in due time
          cycle
        endif ! of if (igcm_h2o_ice.eq.igcm_h2o_vap)
      endif ! of if (txt.eq."h2o_vap")

      if (txt.eq."hdo_vap") then
        write(*,*)"physdem1: skipping HDO vapour tracer"
        if (i_hdo_ice.eq.i_hdo_vap) then
          ! then there is no "water ice" tracer; but still
          ! there is some water ice on the surface
          write(*,*)"          writing HDO ice instead"
          txt="hdo_ice"
        else
          ! there is a "water ice" tracer which has been / will be
          ! delt with in due time
          cycle
        endif ! of if (igcm_hdo_ice.eq.igcm_hdo_vap)
      endif ! of if (txt.eq."hdo_vap")

      ! co2_ice has been added to qsurf(:,igcm_co2) in co2condens4micro
      if (txt.eq."co2_ice") then
        write(*,*)"physdem1: skipping co2_ice tracer"
        cycle
      end if     

      call put_field(trim(txt),"tracer on surface",qsurf(:,iq,:),time)
    enddo
  endif

  ! Non-orographic gavity waves
  if (calllott_nonoro) then
     call put_field("du_nonoro_gwd","Zonal wind tendency due to GW",du_nonoro_gwd,time)
     call put_field("dv_nonoro_gwd","Meridional wind tendency due to GW",dv_nonoro_gwd,time)
  endif


  ! Geothermal Flux
     call put_field('flux_geo','Geothermal flux',flux_geo,time)
  ! Close file
  call close_restartphy
  
end subroutine physdem1

end module phyredem