MODULE dyn_ss_ice_m_mod

implicit none

CONTAINS

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!
!!! Purpose:  Retreat and growth of subsurface ice on Mars
!!!           orbital elements remain constant
!!!
!!!
!!! Author: EV, updated NS MSIM dynamical program for the PEM
!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

SUBROUTINE dyn_ss_ice_m(ssi_depth_in,T1,Tb,nz,thIn,p0,pfrost,porefill_in,porefill,ssi_depth)

!***********************************************************************
! Retreat and growth of subsurface ice on Mars
! orbital elements remain constant
!***********************************************************************
  use constants_marspem_mod, only: sec_per_sol
  use fast_subs_mars, only: psv, icelayer_mars, NMAX
#ifndef CPP_STD
    use comcstfi_h,   only: pi
#else
    use comcstfi_mod, only: pi
#endif
  implicit none
  integer, parameter :: NP=1   ! # of sites
  integer nz, i, k, iloop
  real(8)  zmax, delta, z(NMAX), icetime, porosity, icefrac
  real(8), dimension(NP) ::  albedo, thIn, rhoc
  real(8), dimension(NP) :: pfrost, p0
  real(8) newti, stretch, newrhoc, ecc, omega, eps, timestep
  real(8)  ssi_depth_in, ssi_depth, T1
  real(8), dimension(NP) :: zdepthF, zdepthE, zdepthT, zdepthG
  real(8), dimension(NMAX,NP) :: porefill, porefill_in
  real(8), dimension(nz) :: Tb
  real(8), dimension(NP) :: Tmean1, Tmean3, avrho1
  real(8) tmax, tlast, avrho1prescribed(NP), l1
  real(8), external :: smartzfac

  !if (iargc() /= 1) then
  !   stop 'USAGE: icages ext'
  !endif
  !call getarg( 1, ext )

  if (NP>100) stop 'subroutine icelayer_mars cannot handle this many sites'

  ! parameters that never ever change
  porosity = 0.4d0  ! porosity of till
  !rhoc(:) = 1500.*800.  ! will be overwritten
  icefrac = 0.98
  tmax = 1
  tlast = 0.
  avrho1prescribed(:) = pfrost/T1  ! <0 means absent
  albedo=0.23
  !avrho1prescribed(:) = 0.16/200.  ! units are Pa/K

  !open(unit=21,file='lats.'//ext,action='read',status='old',iostat=ierr)
  !if (ierr /= 0) then
  !   print *,'File lats.'//ext,'not found'
  !   stop
  !endif
  do k=1,NP
     !read(21,*) latitude(k),albedo(k),thIn(k),htopo(k)
     ! empirical relation from Mellon & Jakosky
     rhoc(k) = 800.*(150.+100.*sqrt(34.2+0.714*thIn(k)))
  enddo
  !close(21)

  ! set eternal grid
  zmax = 25.
  !zfac = smartzfac(nz,zmax,6,0.032d0)
  !call setgrid(nz,z,zmax,zfac)
  l1=2.e-4
  do iloop=0,nz - 1
    z(iloop + 1) = l1*(1+iloop**2.9*(1-exp(-real(iloop)/20.)))
  enddo


  !open(unit=30,file='z.'//ext,action='write',status='unknown')
  !write(30,'(999(f8.5,1x))') z(1:nz)
  !close(30)

  !ecc = ecc_in;  eps = obl_in*d2r;  omega = Lp_in*d2r   ! today
  ! total atmospheric pressure
  !p0(:) = 600.
  ! presently 520 Pa at zero elevation (Smith & Zuber, 1998)
 ! do k=1,NP
 !    p0(k)=520*exp(-htopo(k)/10800.)
 ! enddo
  timestep = 1 ! must be integer fraction of 1 ka
  icetime = -tmax-timestep  ! earth years

  ! initializations
  !Tb = -9999.
  zdepthF(:) = -9999.

  !zdepthT(1:NP) = -9999.   ! reset again below
 ! zdepthT(1:NP) = 0.

 ! print *,'RUNNING MARS_FAST'
 ! print *,'Global model parameters:'
 ! print *,'nz=',nz,' zfac=',zfac,'zmax=',zmax
 ! print *,'porosity=',porosity
 ! print *,'starting at time',icetime,'years'
 ! print *,'time step=',timestep,'years'
 ! print *,'eps=',eps/d2r,'ecc=',ecc,'omega=',omega/d2r
 ! print *,'number of sites=',NP
 ! print *,'Site specific parameters:'
  do k=1,NP
     if (NP>1) print *,'  Site ',k
 !    print *,'  latitude (deg)',latitude(k),' rho*c (J/m^3/K)',rhoc(k),' thIn=',thIn(k)
 !    print *,'  total pressure=',p0(k),'partial pressure=',pfrost(k)
     delta = thIn(k)/rhoc(k)*sqrt(sec_per_sol/pi)
 !    print *,'  skin depths (m)',delta,delta*sqrt(solsperyear)
     call soilthprop(porosity,1.d0,rhoc(k),thIn(k),1,newrhoc,newti,icefrac)
     stretch = (newti/thIn(k))*(rhoc(k)/newrhoc)
     do i=1,nz
        if (z(i)<delta) cycle
 !       print *,'  ',i-1,' grid points within diurnal skin depth'
        exit
     enddo
 !    print *,'  ',zmax/(sqrt(solsperyear)*delta),'times seasonal dry skin depth'
 !    print *,'  ',zmax/(sqrt(solsperyear)*delta*stretch),'times seasonal filled skin depth'
 !    print *,'  Initial ice depth=',zdepthT(k)
 !    print *
  enddo
 ! call outputmoduleparameters
 ! print *

  ! open and name all output files
!  open(unit=34,file='subout.'//ext,action='write',status='unknown')
!  open(unit=36,file='depthF.'//ext,action='write',status='unknown')
!  open(unit=37,file='depths.'//ext,action='write',status='unknown')

 ! print *,'Equilibrating initial temperature'
 ! do i=1,4
 !    call icelayer_mars(0d0,nz,NP,thIn,rhoc,z,porosity,pfrost,Tb,zdepthF, &
 !      &  zdepthE,porefill(1:nz,:),Tmean1,Tmean3,zdepthG, &
 !      &  latitude,albedo,p0,ecc,omega,eps,icefrac,zdepthT,avrho1, &
 !      &  avrho1prescribed)
 ! enddo

  !print *,'History begins here'
 porefill(1:nz,1:NP) =  porefill_in(1:nz,1:NP)
  zdepthT(1:NP) = ssi_depth_in
  do
  !print *,'Zt0=  ',ZdepthT
     call icelayer_mars(timestep,nz,NP,thIn,rhoc,z,porosity,pfrost,Tb,zdepthF, &
          & zdepthE,porefill(1:nz,:),Tmean1,Tmean3,zdepthG, &
          & albedo,p0,icefrac,zdepthT,avrho1, &
          & avrho1prescribed)
     icetime = icetime+timestep
   !  print *,'T_after=  ',Tb(:)
 !    print *,'z=  ',z(:)
 !    print *,'Zt=  ',ZdepthT
     ssi_depth=ZdepthT(1)
    ! if (abs(mod(icetime/100.,1.d0))<1.e-3) then ! output every 1000 years
    !    do k=1,NP
         !write(36,*) icetime,latitude(k),zdepthF(k),porefill(1:nz,k)
           ! compact output format
 !          write(36,'(f10.0,2x,f7.3,1x,f11.5,1x)',advance='no') &
 !               & icetime,latitude(k),zdepthF(k)
 !          call compactoutput(36,porefill(:,k),nz)
 !          write(37,501) icetime,latitude(k),zdepthT(k), &
 !               & Tmean1(k),Tmean3(k),zdepthG(k),avrho1(k)
  !      enddo
  !   endif
!     print *,icetime
     if (icetime>=tlast) exit
  enddo

 ! close(34)
 ! close(36); close(37)

!501 format (f10.0,2x,f7.3,2x,f10.4,2(2x,f6.2),2x,f9.3,2x,g11.4)

end subroutine dyn_ss_ice_m

!=======================================================================

subroutine soilthprop(porosity,fill,rhocobs,tiobs,layertype, &
     &     newrhoc,newti,icefrac)
!***********************************************************************
! soilthprop: assign thermal properties of icy soil or dirty ice
!
!     porositiy = void space / total volume
!     rhof = density of free ice in space not occupied by regolith [kg/m^3]
!     fill = rhof/icedensity <=1 (only relevant for layertype 1)
!     rhocobs = heat capacity per volume of dry regolith [J/m^3]
!     tiobs = thermal inertia of dry regolith [SI-units]
!     layertype: 1=interstitial ice, 2=pure ice or ice with dirt
!                3=pure ice, 4=ice-cemented soil, 5=custom values
!     icefrac: fraction of ice in icelayer (only relevant for layertype 2)
!     output are newti and newrhoc
!***********************************************************************
  implicit none
  integer, intent(IN) :: layertype
  real(8), intent(IN) :: porosity, fill, rhocobs, tiobs
  real(8), intent(OUT) :: newti, newrhoc
  real(8), intent(IN) :: icefrac
  real(8) kobs, cice, icedensity, kice
  !parameter (cice=2000.d0, icedensity=926.d0, kice=2.4d0) ! unaffected by scaling
  parameter (cice=1540.d0, icedensity=927.d0, kice=3.2d0) ! at 198 Kelvin
  real(8) fA, ki0, ki, k
  real(8), parameter :: kw=3. ! Mellon et al., JGR 102, 19357 (1997)

  kobs = tiobs**2/rhocobs
  ! k, rhoc, and ti are defined in between grid points
  ! rhof and T are defined on grid points

  newrhoc = -9999.
  newti  = -9999.

  select case (layertype)
  case (1) ! interstitial ice
     newrhoc = rhocobs + porosity*fill*icedensity*cice
     if (fill>0.) then
        !--linear addition (option A)
        k = porosity*fill*kice + kobs
        !--Mellon et al. 1997 (option B)
        ki0 = porosity/(1/kobs-(1-porosity)/kw)
        fA = sqrt(fill)
        ki = (1-fA)*ki0 + fA*kice
        !k = kw*ki/((1-porosity)*ki+porosity*kw)
     else
        k = kobs
     endif
     newti = sqrt(newrhoc*k)

  case (2)  ! massive ice (pure or dirty ice)
     newrhoc = rhocobs*(1.-icefrac)/(1.-porosity) + icefrac*icedensity*cice
     k = icefrac*kice + (1.-icefrac)*kw
     newti = sqrt(newrhoc*k)

  case (3)  ! all ice, special case of layertype 2, which doesn't use porosity
     newrhoc = icedensity*cice
     k = kice
     newti = sqrt(newrhoc*k)

  case (4)  ! pores completely filled with ice, special case of layertype 1
     newrhoc = rhocobs + porosity*icedensity*cice
     k = porosity*kice + kobs ! option A, end-member case of type 1, option A
     !k = kw*kice/((1-porosity)*kice+porosity*kw) ! option B, harmonic average
     newti = sqrt(newrhoc*k)

  case (5)  ! custom values
     ! values from Mellon et al. (2004) for ice-cemented soil
     newrhoc = 2018.*1040.
     k = 2.5
     newti = sqrt(newrhoc*k)

  case default
     error stop 'invalid layer type'

  end select

end subroutine soilthprop


!=======================================================================

      real*8 function frostpoint(p)
!     inverse of psv
!     input is partial pressure [Pascal]
!     output is temperature [Kelvin]
      implicit none
      real*8 p

!-----inverse of parametrization 1
!      real*8 DHmelt,DHvap,DHsub,R,pt,Tt
!      parameter (DHmelt=6008.,DHvap=45050.)
!      parameter (DHsub=DHmelt+DHvap)
!      parameter (R=8.314,pt=6.11e2,Tt=273.16)
!      frostpoint = 1./(1./Tt-R/DHsub*log(p/pt))

!-----inverse of parametrization 2
!     inverse of eq. (2) in Murphy & Koop (2005)
      real*8 A,B
      parameter (A=-6143.7, B=28.9074)
      frostpoint = A / (log(p) - B)

!-----approximate inverse of parametrization 3
!     eq. (8) in Murphy & Koop (2005)
!      frostpoint = (1.814625*log(p) + 6190.134)/(29.120 - log(p))

      end function frostpoint

END MODULE dyn_ss_ice_m_mod