MODULE subsurface_ice
!-----------------------------------------------------------------------
! NAME
!     subsurface_ice
!
! DESCRIPTION
!     Retreat and growth of subsurface ice on Mars with constant orbital
!     elements.
!
! AUTHORS & DATE
!     N. Schorghofer (MSIM dynamical program)
!     E. Vos, 2025
!     JB Clement, 2025
!
! NOTES
!    Based on Norbert Schorgofer's code. Updated for PEM.
!-----------------------------------------------------------------------

! DECLARATION
! -----------
implicit none

! MODULE PARAMETERS
! -----------------
real(8), parameter :: dt = 0.02  ! in units of Mars solar days
!real(8), parameter :: Fgeotherm = 0.
real(8), parameter :: Fgeotherm = 0 !0.028  ! [W/m^2]
real(8), parameter :: Lco2frost=6.0e5, co2albedo=0.60, co2emiss=1.
real(8), parameter :: emiss0 = 1.  ! emissivity of dry surface
integer, parameter :: EQUILTIME =15 ! [Mars years]

integer, parameter :: NMAX = 1000

contains
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++


!=======================================================================
!============================ dyn_ss_ice_m =============================
!=======================================================================


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 phys_constants,        only: pi
  ! DECLARATION
! -----------
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(nz,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,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
!***********************************************************************
  ! DECLARATION
! -----------
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]
      ! DECLARATION
! -----------
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


!=======================================================================
!=========================== fast_subs_mars ============================
!=======================================================================


!*****************************************************
! Subroutines for fast method
! written by Norbert Schorghofer 2007-2011
!*****************************************************


subroutine icelayer_mars(bigstep,nz,NP,thIn,rhoc,z,porosity,pfrost, &
     & Tb,zdepthF,zdepthE,porefill,Tmean1,Tmean3,zdepthG, &
     & albedo,p0,icefrac,zdepthT,avrho1, &
     & avrho1prescribed)
!*************************************************************************
! bigstep = time step [Earth years]
! latitude  [degree]
!*************************************************************************
  use ice_table,      only: rho_ice
  !use omp_lib
  ! DECLARATION
! -----------
implicit none
  integer, intent(IN) :: nz, NP
  real(8), intent(IN) :: bigstep
  real(8), intent(IN) :: thIn(NP), rhoc(NP), z(NMAX), porosity, pfrost(NP)
  real(8), intent(INOUT) :: porefill(nz,NP), zdepthF(NP), zdepthT(NP)
  real(8), intent(INOUT) :: Tb(nz)
  real(8), intent(OUT), dimension(NP) :: zdepthE, Tmean1, Tmean3, zdepthG
  real(8), intent(IN), dimension(NP) ::  albedo, p0
  real(8), intent(IN) ::  icefrac
  real(8), intent(OUT) :: avrho1(NP)
  real(8), intent(IN), optional :: avrho1prescribed(NP)  ! <0 means absent

  integer k, typeF, typeG, typeT, j, jump, typeP
  real(8) fracIR, fracDust, ti(NMAX), rhocv(NMAX)
  real(8) Diff, ypp(nz), avdrho(NP), avdrhoP(NP), B, deltaz
  real(8), SAVE :: avdrho_old(100), zdepth_old(100)  ! NP<=100
  logical mode2

  avdrho_old = 1. ! initialization

  !$omp parallel &
  !$omp    private(Diff,fracIR,fracDust,B,typeT,j,ti,rhocv,typeF,jump,typeG)
  !$omp do
  do k=1,NP   ! big loop

     Diff = 4e-4*600./p0(k)
     fracIR = 0.04*p0(k)/600.; fracDust = 0.02*p0(k)/600.
     B = Diff*bigstep*86400.*365.24/(porosity*927.)
     !B = Diff*bigstep*86400.*365.24/(porosity*rho_ice(Tb(),'h2o'))

     typeT = -9
     if (zdepthT(k)>=0. .and. zdepthT(k)<z(nz)) then
        do j=1,nz
           if (z(j)>zdepthT(k)) then ! ice
              typeT = j  ! first point in ice
              exit
           endif
        enddo
     endif

   !  call assignthermalproperties(nz,thIn(k),rhoc(k),ti,rhocv, &
   !       & typeT,icefrac,porosity,porefill(:,k))

     !----run thermal model
     call ajsub_mars(typeT, albedo(k), pfrost(k), nz, z, &
          &     ti, rhocv, fracIR, fracDust, p0(k), &
          &     avdrho(k), avdrhoP(k), avrho1(k), Tb(k), zdepthE(k), typeF, &
          &     zdepthF(k), ypp, porefill(:,k), Tmean1(k), Tmean3(k), B, &
          &     typeG, zdepthG(k), avrho1prescribed(k))

     if (typeF*zdepthF(k)<0.) stop 'error in icelayer_mars'
     ! diagnose
     if (zdepthT(k)>=0.) then
        jump = 0
        do j=1,nz
           if (zdepth_old(k)<z(j) .and. zdepthT(k)>z(j)) jump = jump+1
        enddo
     else
        jump = -9
     endif
     if (zdepthT(k)>=0. .and. avdrho(k)*avdrho_old(k)<0.) then
        write(34,*) '# zdepth arrested'
        if (jump>1) write(34,*) '# previous step was too large',jump
     endif
!     write(34,'(f8.3,1x,f6.2,1x,f11.5,2(1x,g11.4),1x,i3)') &
!          &        bigstep,latitude(k),zdepthT(k),avdrho(k),avdrhoP(k),jump
     zdepth_old(k) = zdepthT(k)
     avdrho_old(k) = avdrho(k)

!----mode 2 growth
     typeP = -9;  mode2 = .FALSE.
     do j=1,nz
        if (porefill(j,k)>0.) then
           typeP = j  ! first point with ice
           exit
        endif
     enddo
     if (typeT>0 .and. typeP>2 .and. zdepthE(k)>0.) then
        if (porefill(typeP,k)>=1. .and. porefill(typeP-1,k)==0. .and. &
             & zdepthE(k)<z(typeP) .and. &
             & z(typeP)-zdepthE(k)>2*(z(typeP)-z(typeP-1))) then  ! trick that avoids oscillations
           deltaz = -avdrhoP(k)/z(typeP)*18./8314.*B  ! conservation of mass
           if (deltaz>z(typeP)-z(typeP-1)) then  ! also implies avdrhoP<0.
              mode2 = .TRUE.
           endif
        endif
     endif

     !call icechanges_poreonly(nz,z,typeF,typeG,avdrhoP(k),ypp,B,porefill(:,k))
     call icechanges(nz,z(:),typeF,avdrho(k),avdrhoP(k),ypp(:), &
          & Diff,porosity,icefrac,bigstep,zdepthT(k),porefill(:,k),typeG)
     if (typeP>2) then
     if (mode2 .and. porefill(typeP,k)>=1. .and. porefill(typeP-1,k)==0.) then  ! nothing changed
        porefill(typeP-1,k)=1.  ! paste a layer
   !     write(34,*) '# mode 2 growth occurred',typeP,typeF,typeT
     endif
     endif

  enddo  ! end of big loop
  !$omp end do
  !$omp end parallel
end subroutine icelayer_mars



subroutine ajsub_mars(typeT, albedo0, pfrost, nz, z, ti, rhocv, &
     &     fracIR, fracDust, patm, avdrho, avdrhoP, avrho1, &
     &     Tb, zdepthE, typeF, zdepthF, ypp, porefill, Tmean1, Tmean3, &
     &     B, typeG, zdepthG, avrho1prescribed)
!***********************************************************************
!  A 1D thermal model that returns various averaged quantities
!
!  Tb<0 initializes temperatures
!  Tb>0 initializes temperatures with Tb
!***********************************************************************
  use constants_marspem_mod, only: sols_per_my, sec_per_sol
  ! DECLARATION
! -----------
implicit none
  integer, intent(IN) :: nz, typeT
!  real(8), intent(IN) :: latitude  ! in radians
  real(8), intent(IN) :: albedo0, pfrost, z(NMAX)
  real(8), intent(IN) :: ti(NMAX), rhocv(NMAX), fracIR, fracDust, patm
  real(8), intent(IN) ::  porefill(nz)
  real(8), intent(OUT) :: avdrho, avdrhoP  ! difference in annual mean vapor density
  real(8), intent(OUT) :: avrho1  ! mean vapor density on surface
  real(8), intent(INOUT) :: Tmean1
  real(8), intent(INOUT) :: Tb(nz)
  integer, intent(OUT) :: typeF  ! index of depth below which filling occurs
  real(8), intent(OUT) :: zdepthE, zdepthF
  real(8), intent(OUT) :: ypp(nz) ! (d rho/dt)/Diff
  real(8), intent(OUT) :: Tmean3, zdepthG
  real(8), intent(IN) :: B  ! just passing through
  integer, intent(OUT) :: typeG
  real(8), intent(IN), optional :: avrho1prescribed
  real(8), parameter :: a = 1.52366 ! Mars semimajor axis in A.U.
  integer nsteps, n, i, nm, typeP
  real(8) tmax, time, Qn, Qnp1, tdays
  real(8) marsR, marsLs, marsDec, HA
  real(8) Tsurf, Tco2frost, albedo, Fsurf, m, dE, emiss, T(NMAX)
  real(8) Told(nz), Fsurfold, Tsurfold, Tmean0, avrho2
  real(8) rhosatav0, rhosatav(nz), rlow

  tmax = EQUILTIME*sols_per_my
  nsteps = int(tmax/dt)     ! calculate total number of timesteps

  Tco2frost = tfrostco2(patm)

  !if (Tb<=0.) then  ! initialize
     !Tmean0 = 210.15       ! black-body temperature of planet
     !Tmean0 = (589.*(1.-albedo0)*cos(latitude)/(pi*emiss0*sigSB))**0.25 ! better estimate
     !Tmean0 = Tmean0-5.
     !write(34,*) '# initialized with temperature estimate at',latitude/d2r,'of',Tmean0,'K'
     !T(1:nz) = Tmean0
  !else
     T(1:nz) = Tb
     ! not so good when Fgeotherm is on
  !endif

  albedo = albedo0
  emiss = emiss0
  do i=1,nz
     if (T(i)<Tco2frost) T(i)=Tco2frost
  enddo
  Tsurf = T(1)
  m=0.; Fsurf=0.

  nm=0
  avrho1=0.; avrho2=0.
  Tmean1=0.; Tmean3=0.
  rhosatav0 = 0.
  rhosatav(:) = 0.

  time=0.
  !call generalorbit(0.d0,a,ecc,omega,eps,marsLs,marsDec,marsR)
  !HA = 2.*pi*time            ! hour angle
!  Qn=flux(marsR,marsDec,latitude,HA,albedo,fracir,fracdust,0.d0,0.d0)
  !Qn = flux_mars77(marsR,marsDec,latitude,HA,albedo,fracir,fracdust)
  !----loop over time steps
  do n=0,nsteps-1
     time = (n+1)*dt         !   time at n+1
   !  tdays = time*(sec_per_sol/earthDay) ! parenthesis may improve roundoff
   !  call generalorbit(tdays,a,ecc,omega,eps,marsLs,marsDec,marsR)
   !  HA = 2.*pi*mod(time,1.d0)  ! hour angle
!     Qnp1=flux(marsR,marsDec,latitude,HA,albedo,fracir,fracdust,0.d0,0.d0)
     !Qnp1 = flux_mars77(marsR,marsDec,latitude,HA,albedo,fracir,fracdust)

     Tsurfold = Tsurf
     Fsurfold = Fsurf
     Told(1:nz) = T(1:nz)
     if (m<=0. .or. Tsurf>Tco2frost) then
       ! call conductionQ(nz,z,dt*sec_per_sol,Qn,Qnp1,T,ti,rhocv,emiss, &
       !      &           Tsurf,Fgeotherm,Fsurf)
     endif
     if (Tsurf<Tco2frost .or. m>0.) then ! CO2 condensation
        T(1:nz) = Told(1:nz)
        !call conductionT(nz,z,dt*sec_per_sol,T,Tsurfold,Tco2frost,ti, &
             !&              rhocv,Fgeotherm,Fsurf)
        Tsurf = Tco2frost
    !    dE = (- Qn - Qnp1 + Fsurfold + Fsurf + &
    !         &      emiss*sigSB*(Tsurfold**4+Tsurf**4)/2.
        m = m + dt*sec_per_sol*dE/Lco2frost
     endif
     if (Tsurf>Tco2frost .or. m<=0.) then
        albedo = albedo0
        emiss = emiss0
     else
        albedo = co2albedo
        emiss = co2emiss
     endif
     !Qn=Qnp1

     if (time>=tmax-sols_per_my) then
        Tmean1 = Tmean1 + Tsurf
        Tmean3 = Tmean3 + T(nz)
        avrho1 = avrho1 + min(psv(Tsurf),pfrost)/Tsurf
        rhosatav0 = rhosatav0 + psv(Tsurf)/Tsurf
        do i=1,nz
           rhosatav(i) = rhosatav(i) + psv(T(i))/T(i)
        enddo
        nm = nm+1
     endif

  enddo  ! end of time loop

  avrho1 = avrho1/nm
  if (present(avrho1prescribed)) then
     if (avrho1prescribed>=0.) avrho1=avrho1prescribed
  endif
  Tmean1 = Tmean1/nm; Tmean3 = Tmean3/nm
  rhosatav0 = rhosatav0/nm; rhosatav(:)=rhosatav(:)/nm
  if (typeT>0) then
     avrho2 = rhosatav(typeT)
  else
     avrho2 = rhosatav(nz) ! no ice
  endif
  avdrho = avrho2-avrho1
  typeP = -9
  do i=1,nz
     if (porefill(i)>0.) then
        typeP = i  ! first point with ice
        exit
     endif
  enddo
  if (typeP>0) then
     avdrhoP = rhosatav(typeP) - avrho1
  else
     avdrhoP = -9999.
  end if

  zdepthE = equildepth(nz, z, rhosatav, rhosatav0, avrho1)

  if (Fgeotherm>0.) then
     !Tb = Tmean1
     typeG = 1   ! will be overwritten by depths_avmeth
     rlow = 2*rhosatav(nz)-rhosatav(nz-1)
  else
     !Tb = T(nz)
     typeG = -9
     rlow = rhosatav(nz-1)
  endif
  call depths_avmeth(typeT, nz, z, rhosatav(:), rhosatav0, rlow, avrho1,  &
       & porefill(:), typeF, zdepthF, B, ypp(:), typeG, zdepthG)

end subroutine ajsub_mars



real*8 function tfrostco2(p)
!     the inverse of function psvco2
!     input is pressure in Pascal, output is temperature in Kelvin
! DECLARATION
! -----------
implicit none
real*8 p
tfrostco2 = 3182.48/(23.3494+log(100./p))
end function

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

      real*8 function psv(T)
!     saturation vapor pressure of H2O ice [Pascal]
!     input is temperature [Kelvin]
      ! DECLARATION
! -----------
implicit none
      real*8 T

!-----parametrization 1
!      real*8 DHmelt,DHvap,DHsub,R,pt,Tt,C
!      parameter (DHmelt=6008.,DHvap=45050.)
!      parameter (DHsub=DHmelt+DHvap) ! sublimation enthalpy [J/mol]
!      parameter (R=8.314,pt=6.11e2,Tt=273.16)
!      C = (DHsub/R)*(1./T - 1./Tt)
!      psv = pt*exp(-C)

!-----parametrization 2
!     eq. (2) in Murphy & Koop, Q. J. R. Meteor. Soc. 131, 1539 (2005)
!     differs from parametrization 1 by only 0.1%
      real*8 A,B
      parameter (A=-6143.7, B=28.9074)
      psv = exp(A/T+B)  ! Clapeyron

!-----parametrization 3
!     eq. (7) in Murphy & Koop, Q. J. R. Meteor. Soc. 131, 1539 (2005)
!     psv = exp(9.550426 - 5723.265/T + 3.53068*log(T) - 0.00728332*T)

      end function psv


!=======================================================================
!=========================== fast_subs_univ ============================
!=======================================================================


!*****************************************************
! Commonly used subroutines for fast method
! written by Norbert Schorghofer 2007-2011
!*****************************************************

pure function zint(y1,y2,z1,z2)
  ! interpolate between two values, y1*y2<0
  ! DECLARATION
! -----------
implicit none
  real(8), intent(IN) :: y1,y2,z1,z2
  real(8) zint
  zint = (y1*z2 - y2*z1)/(y1-y2)
end function zint



pure function equildepth(nz, z, rhosatav, rhosatav0, avrho1)
!***********************************************************************
!  returns equilibrium depth for given ice content
!  this is not the true (final) equilibrium depth
!***********************************************************************
  ! DECLARATION
! -----------
implicit none
  integer, intent(IN) :: nz
  real(8), intent(IN) :: z(nz), rhosatav(nz)
  real(8), intent(IN) :: rhosatav0, avrho1
  integer i, typeE
  real(8) equildepth  ! =zdepthE
  !real(8), external :: zint  ! defined in allinterfaces.mod

  typeE = -9; equildepth = -9999.
  do i=1,nz
     if (rhosatav(i) <= avrho1) then
        typeE=i
        exit
     endif
  enddo
  if (typeE>1) then ! interpolate
     equildepth=zint(avrho1-rhosatav(typeE-1),avrho1-rhosatav(typeE),z(typeE-1),z(typeE))
  end if
  if (typeE==1) equildepth=zint(avrho1-rhosatav0,avrho1-rhosatav(1),0.d0,z(1))
  if (equildepth>z(nz)) equildepth=-9999.   ! happens very rarely
end function equildepth



subroutine depths_avmeth(typeT, nz, z, rhosatav, rhosatav0, rlow, avrho1,  &
     & porefill, typeF, zdepthF, B, ypp, typeG, zdepthG)
!***********************************************************************
!  returns interface depth and ypp
!  also returns lower geothermally limited boundary, if applicable
!
!  this is a crucial subroutine
!
!  B = Diff*bigstep/(porosity*icedensity)  [SI units]
!***********************************************************************
  use math_toolkit, only: deriv2_simple, deriv1_onesided, deriv1, colint
  use ice_table,    only: constriction
  ! DECLARATION
! -----------
implicit none
  integer, intent(IN) :: nz, typeT
  real(8), intent(IN), dimension(nz) :: z, rhosatav, porefill
  real(8), intent(IN) :: rhosatav0, rlow, avrho1
  integer, intent(OUT) :: typeF  ! index of depth below which filling occurs
  real(8), intent(INOUT) :: zdepthF
  real(8), intent(IN) :: B
  real(8), intent(OUT) :: ypp(nz), zdepthG
  integer, intent(INOUT) :: typeG  ! positive on input when Fgeotherm>0

  integer i, typeP, nlb, newtypeG
  real(8) eta(nz), Jpump1, help(nz), yp(nz), zdepthFold, ap_one, ap(nz)
  real(8) leak, cumfill, cumfillabove

  if (typeT<0) then
     nlb = nz
     do i=1,nz
        eta(i) = constriction(porefill(i))
     enddo
  else
     !nlb = typeT-1
     nlb = typeT ! changed 2010-09-29
     do i=1,typeT-1
        eta(i) = constriction(porefill(i))
     enddo
     do i=typeT,nz
        eta(i)=0.
     enddo
  end if

!-fill depth
  zdepthFold = zdepthF
  typeF = -9;  zdepthF = -9999.
  call deriv1(z,nz,rhosatav,rhosatav0,rlow,yp)  ! yp also used below
  do i=1,nlb
     Jpump1 = (rhosatav(i)-avrho1)/z(i)  ! <0 when stable
     ! yp is always <0
     help(i) = Jpump1 - eta(i)*yp(i)
     leak = porefill(i)/B*(z(i)-zdepthFold)/(18./8314.)
     !help(i) = help(i)-leak   ! optional
     if (help(i) <= 0.) then
        typeF=i
        !print *,'#',typeF,Jpump1,eta(typeF)*yp(typeF),leak
        exit
     endif
  enddo
  if (typeF>1) zdepthF = zint(help(typeF-1),help(typeF),z(typeF-1),z(typeF))
  if (typeF==1) zdepthF=z(1)


!-depth to shallowest perennial ice
  typeP = -9
  do i=1,nz
     if (porefill(i)>0.) then
        typeP = i  ! first point with ice
        exit
     endif
  enddo

!-calculate ypp
  !call deriv1(z,nz,rhosatav,rhosatav0,rlow,yp)
  call deriv1(z,nz,eta(:),1.d0,eta(nz-1),ap)
  if (typeP>0 .and. typeP<nz-2) then
     call deriv1_onesided(typeP,z,nz,eta(:),ap_one)
     ! print *,typeP,ap(typeP),ap_one
     ap(typeP)=ap_one
  endif
  call deriv2_simple(z,nz,rhosatav(1:nz),rhosatav0,rlow,ypp(:))
  !call deriv2_full(z,nz,eta(:),rhosatav(:),1.d0,rhosatav0,rhosatav(nz-1),ypp(:))
  !write(40,*) rhosatav
  !write(41,*) yp
  !write(42,*) ypp
  ypp(:) = ap(:)*yp(1:)+eta(:)*ypp(:)
  !write(43,*) ypp
  !write(44,*) eta(1:nz)
  !write(45,*) (rhosatav(:)-avrho1)/z(:)
  ypp(:) = ypp(:)*18./8314.
  ! ypp values beyond nlb should not be used

!-geothermal stuff
  if (typeT>0) typeG=-9
  if (typeG<0) zdepthG=-9999.
  if (typeG>0 .and. typeT<0) then
     typeG=-9
     do i=2,nz
        if (yp(i)>0.) then  ! first point with reversed flux
           typeG=i
           zdepthG=zint(yp(i-1),yp(i),z(i-1),z(i))
           !zdepthG=z(typeG)
           exit
        endif
     enddo
  else
     typeG = -9
  endif
  if (typeG>0 .and. typeT<0) then
     call colint(porefill(:)/eta(:),z,nz,typeG-1,nz,cumfillabove)
     newtypeG = -9
     do i=typeG,nz
        if (minval(eta(i:nz))<=0.) cycle  ! eta=0 means completely full
        call colint(porefill(:)/eta(:),z,nz,i,nz,cumfill)
        if (cumfill<yp(i)*18./8314.*B) then  ! usually executes on i=typeG
           if (i>typeG) then
              write(34,*) '# adjustment to geotherm depth by',i-typeG
              zdepthG = zint(yp(i-1)*18./8314.*B-cumfillabove, &  ! no good
                   &        yp(i)*18./8314.*B-cumfill,z(i-1),z(i))
              if (zdepthG>z(i) .or. zdepthG<z(i-1)) write(34,*) &
                   & '# WARNING: zdepthG interpolation failed',i,z(i-1),zdepthG,z(i)
              newtypeG=i
           endif
           ! otherwise leave zdepthG untouched
           exit
        endif
        cumfillabove = cumfill
     enddo
     if (newtypeG>0) typeG=newtypeG
  end if
  ! if typeG>0, then all ice at and below typeG should be erased
end subroutine depths_avmeth



pure subroutine icechanges(nz,z,typeF,avdrho,avdrhoP,ypp, &
     & Diff,porosity,icefrac,bigstep,zdepthT,porefill,typeG)
!***********************************************************
! advances ice table, advances interface, grows pore ice
!
! a crucial subroutine
!***********************************************************
  use math_toolkit, only: colint
  use ice_table,    only: rho_ice
  ! DECLARATION
! -----------
implicit none
  integer, intent(IN) :: nz, typeF, typeG
  real(8), intent(IN) :: z(nz), ypp(nz), avdrho, avdrhoP
  real(8), intent(IN) :: Diff, porosity, icefrac, bigstep
  real(8), intent(INOUT) :: zdepthT, porefill(nz)
  integer j, erase, newtypeP, ub, typeP, typeT
  real(8) B, beta, integ

  B = Diff*bigstep*86400.*365.24/(porosity*927.)
  !B = Diff*bigstep*86400.*365.24/(porosity*rho_ice(T,'h2o'))

  ! advance ice table, avdrho>0 is retreat
  if (zdepthT>=0. .and. avdrho>0.) then
     typeP=-9999; typeT=-9999
     do j=1,nz
        if (z(j)>zdepthT) then
           typeT=j
           exit
        endif
     enddo
     do j=1,nz
        if (porefill(j)>0.) then
           typeP=j
           exit
        endif
     enddo
     if (typeP==typeT) then   ! new 2011-09-01
        beta = (1-icefrac)/(1-porosity)/icefrac
        beta = Diff*bigstep*beta*86400*365.24/927.
        !beta = Diff*bigstep*beta*86400*365.24/rho_ice(T,'h2o')
        zdepthT = sqrt(2*beta*avdrho*18./8314. + zdepthT**2)
     endif
  endif
  if (zdepthT>z(nz)) zdepthT=-9999.

  ! advance interface, avdrhoP>0 is loss from zdepthP
  if (avdrhoP>0.) then
     erase=0
     do j=1,nz
        if (typeF>0 .and. j>=typeF) exit ! don't retreat beyond typeF
        if (zdepthT>=0. .and. z(j)>zdepthT) exit
        call colint(porefill(1:nz)*z(1:nz),z(1:nz),nz,1,j,integ)
        erase = j
        if (integ>B*avdrhoP*18./8314.) exit
     end do
     if (erase>0) porefill(1:erase)=0.
  endif

  ! new depth
  newtypeP = -9
  do j=1,nz
     if (zdepthT>=0. .and. z(j)>zdepthT) exit
     if (porefill(j)>0.) then
        newtypeP = j  ! first point with pore ice
        exit
     endif
  enddo

  ! diffusive filling
  ub = typeF
  if (newtypeP>0 .and. typeF>0 .and. newtypeP<ub) ub=newtypeP
  if (ub>0) then
     do j=ub,nz
        porefill(j) = porefill(j) + B*ypp(j)
        if (porefill(j)<0.) porefill(j)=0.
        if (porefill(j)>1.) porefill(j)=1.
        if (zdepthT>=0. .and. z(j)>zdepthT) exit
     enddo
  end if

  ! below icetable
  if (zdepthT>=0.) then
     do j=1,nz
        if (z(j)>zdepthT) porefill(j) = icefrac/porosity
     enddo
  else
     ! geothermal lower boundary
     if (typeG>0) porefill(typeG:nz)=0.
  end if
end subroutine icechanges

END MODULE subsurface_ice