module surfdat_h

  ! arrays are allocated in conf_phys
  real,save,allocatable :: albedodat(:) ! albedo of bare ground
  real,save,allocatable :: phisfi(:) ! geopotential at ground level
  real,save :: albedice(2) ! default albedo for ice (1: North H. 2: South H.)
  real,save :: emisice(2) ! ice emissivity; 1:Northern hemisphere 2:Southern hemisphere
  real,save :: emissiv ! emissivity of bare ground
  logical,save :: TESicealbedo ! use TES ice cap albedoes (if set to .true.)
  logical,save,allocatable :: watercaptag(:) ! flag for water ice surface
  real, save, allocatable :: dryness(:)

!$OMP THREADPRIVATE(albedodat, phisfi,albedice,emisice,emissiv,TESicealbedo,     &
!$OMP                watercaptag,dryness)
      
  logical,save :: temptag !temp tag for water caps

!$OMP THREADPRIVATE(temptag)
      
  real,save :: albedo_h2o_cap ! water cap albedo
  real,save :: albedo_h2o_frost ! water frost albedo
  real,save :: inert_h2o_ice ! water ice thermal inertia
  real,save :: frost_albedo_threshold ! water frost thickness on the ground (kg.m^-2, ie mm)
  real,save :: frost_metam_threshold ! water frost threshold before conversion to ice (kg.m^-2, ie mm)
  real,save :: TESice_Ncoef ! coefficient for TES ice albedo in Northern hemisphere
  real,save :: TESice_Scoef ! coefficient for TES ice albedo in Southern hemisphere
  real,save :: iceradius(2) , dtemisice(2)
  real,save,allocatable :: zmea(:),zstd(:),zsig(:),zgam(:),zthe(:)
  real,save,allocatable :: hmons(:),summit(:),base(:)
  real,save,allocatable :: z0(:) ! surface roughness length (m)
  real,save :: z0_default ! default (constant over planet) surface roughness (m)

!$OMP THREADPRIVATE(albedo_h2o_cap,albedo_h2o_frost,inert_h2o_ice,               &
!$OMP                frost_albedo_threshold,frost_metam_threshold,TESice_Ncoef,  &
!$OMP                TESice_Scoef,iceradius,dtemisice,                           &
!$OMP                zmea,zstd,zsig,zgam,zthe,hmons,summit,base,z0,z0_default )

  !! mountain top dust flows
  REAL,SAVE,ALLOCATABLE :: alpha_hmons(:) ! sub-grid scale mountain mesh fraction
  REAL,SAVE,ALLOCATABLE :: hsummit(:) ! mountain height above the GCM surface
  LOGICAL,SAVE,ALLOCATABLE :: contains_mons(:) ! is there a mountain in the grid mesh ?
          
!$OMP THREADPRIVATE(alpha_hmons,hsummit,contains_mons)

  !! variables
  REAL,SAVE,ALLOCATABLE :: tsurf(:,:)   ! Surface temperature (K)
  REAL,SAVE,ALLOCATABLE :: emis(:,:)    ! Thermal IR surface emissivity
  REAL,SAVE,ALLOCATABLE :: capcal(:,:) ! surface heat capacity (J m-2 K-1)
  REAL,SAVE,ALLOCATABLE :: fluxgrd(:,:) ! surface conduction flux (W.m-2)
  REAL,ALLOCATABLE,SAVE :: qsurf(:,:,:) ! tracer on surface (e.g. kg.m-2)
  REAL,SAVE,ALLOCATABLE :: watercap(:,:) ! Surface water ice (kg.m-2)
  REAL,SAVE,ALLOCATABLE :: perenial_co2ice(:,:) ! Perenial CO2 ice (kg.m-2)
!$OMP THREADPRIVATE(tsurf,emis,capcal,fluxgrd,qsurf,watercap,perenial_co2ice)

contains

  subroutine ini_surfdat_h(ngrid,nq,nslope)
  
  implicit none
  integer,intent(in) :: ngrid ! number of atmospheric columns
  integer,intent(in) :: nq ! number of tracers  
  integer,intent(in) :: nslope ! number of sub-grid scale slope  
    allocate(albedodat(ngrid))
    allocate(phisfi(ngrid))
    allocate(watercaptag(ngrid))
    allocate(dryness(ngrid))
    allocate(zmea(ngrid))
    allocate(zstd(ngrid))
    allocate(zsig(ngrid))
    allocate(zgam(ngrid))
    allocate(zthe(ngrid))
    allocate(z0(ngrid))
    allocate(qsurf(ngrid,nq,nslope))
    allocate(tsurf(ngrid,nslope))
    allocate(watercap(ngrid,nslope))
    allocate(perenial_co2ice(ngrid,nslope))
    allocate(emis(ngrid,nslope))
    allocate(capcal(ngrid,nslope))
    allocate(fluxgrd(ngrid,nslope))
    allocate(hmons(ngrid))
    allocate(summit(ngrid))
    allocate(base(ngrid))
    allocate(alpha_hmons(ngrid))
    allocate(hsummit(ngrid))
    allocate(contains_mons(ngrid))
       
  end subroutine ini_surfdat_h


  subroutine end_surfdat_h

  implicit none

    if (allocated(albedodat))        deallocate(albedodat)
    if (allocated(phisfi))           deallocate(phisfi)
    if (allocated(watercaptag))      deallocate(watercaptag)
    if (allocated(dryness))          deallocate(dryness)
    if (allocated(zmea))             deallocate(zmea)
    if (allocated(zstd))             deallocate(zstd)
    if (allocated(zsig))             deallocate(zsig)
    if (allocated(zgam))             deallocate(zgam)
    if (allocated(zthe))             deallocate(zthe)
    if (allocated(z0))               deallocate(z0)
    if (allocated(qsurf))            deallocate(qsurf)
    if (allocated(tsurf))            deallocate(tsurf)
    if (allocated(watercap))         deallocate(watercap)
    if (allocated(perenial_co2ice))  deallocate(perenial_co2ice)
    if (allocated(emis))             deallocate(emis)
    if (allocated(capcal))           deallocate(capcal)
    if (allocated(fluxgrd))          deallocate(fluxgrd)
    if (allocated(hmons))            deallocate(hmons)
    if (allocated(summit))           deallocate(summit)
    if (allocated(base))             deallocate(base)
    if (allocated(alpha_hmons))      deallocate(alpha_hmons)
    if (allocated(hsummit))          deallocate(hsummit)
    if (allocated(contains_mons))    deallocate(contains_mons)
    
  end subroutine end_surfdat_h

  subroutine ini_surfdat_h_slope_var(ngrid,nq,nslope)
  
  implicit none
  integer,intent(in) :: ngrid ! number of atmospheric columns
  integer,intent(in) :: nq ! number of tracers  
  integer,intent(in) :: nslope ! number of sub-grid scale slope  
    allocate(qsurf(ngrid,nq,nslope))
    allocate(tsurf(ngrid,nslope))
    allocate(watercap(ngrid,nslope))
    allocate(perenial_co2ice(ngrid,nslope))
    allocate(emis(ngrid,nslope))
    allocate(capcal(ngrid,nslope))
    allocate(fluxgrd(ngrid,nslope))
       
  end subroutine ini_surfdat_h_slope_var

  subroutine end_surfdat_h_slope_var

  implicit none

    if (allocated(qsurf))            deallocate(qsurf)
    if (allocated(tsurf))            deallocate(tsurf)
    if (allocated(watercap))         deallocate(watercap)
    if (allocated(perenial_co2ice))  deallocate(perenial_co2ice)
    if (allocated(emis))             deallocate(emis)
    if (allocated(capcal))           deallocate(capcal)
    if (allocated(fluxgrd))          deallocate(fluxgrd)
    
  end subroutine end_surfdat_h_slope_var

end module surfdat_h