module ice_table_mod

  implicit none

  contains


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!
!!! Purpose: Compute the ice table in two ways: dynamic and at equilibrium
!!! Author:  LL, 02/2023
!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!



   SUBROUTINE computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,rhowatersurf_ave,rhowatersoil_ave,ice_table)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!
!!! Purpose: Compute the ice table depth knowing the yearly average water
!!! density at the surface and at depth.
!!! Computations are made following the methods in Schorgofer et al., 2005
!!! This subroutine only gives the ice table at equilibrium and does not consider exchange with the atmosphere
!!! 
!!! Author: LL
!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

#ifndef CPP_STD
    USE comsoil_h_PEM, only: mlayer_PEM                             ! Depth of the vertical grid 
    implicit none 

    integer,intent(in) :: ngrid,nslope,nsoil_PEM                   ! Size of the physical grid, number of subslope, number of soil layer in the PEM 
    logical,intent(in) :: watercaptag(ngrid)                       ! Boolean to check the presence of a perennial glacier
    real,intent(in) :: rhowatersurf_ave(ngrid,nslope)              ! Water density at the surface, yearly averaged [kg/m^3]
    real,intent(in) :: rhowatersoil_ave(ngrid,nsoil_PEM,nslope)    ! Water density at depth, computed from clapeyron law's (Murchy and Koop 2005), yearly averaged  [kg/m^3]
    real,intent(inout) :: ice_table(ngrid,nslope)                  ! ice table depth [m]
    real :: z1,z2                                                  ! intermediate variables used when doing a linear interpolation between two depths to find the root 
    integer ig, islope,isoil                                       ! loop variables
    real :: diff_rho(nsoil_PEM)                                    ! difference of water vapor density between the surface and at depth [kg/m^3]


     do ig = 1,ngrid
      if(watercaptag(ig)) then
         ice_table(ig,:) = 0.
      else
        do islope = 1,nslope
           ice_table(ig,islope) = -1.

         do isoil = 1,nsoil_PEM
           diff_rho(isoil) = rhowatersurf_ave(ig,islope) - rhowatersoil_ave(ig,isoil,islope)

         enddo
         if(diff_rho(1) > 0) then                                   ! ice is at the surface
           ice_table(ig,islope) = 0.
         else
           do isoil = 1,nsoil_PEM -1                                ! general case, we find the ice table depth by doing a linear approximation between the two depth, and then solve the first degree equation to find the root
             if((diff_rho(isoil).lt.0).and.(diff_rho(isoil+1).gt.0.)) then
               call findroot(diff_rho(isoil),diff_rho(isoil+1),mlayer_PEM(isoil),mlayer_PEM(isoil+1),ice_table(ig,islope))
               exit
             endif !diff_rho(z) <0 & diff_rho(z+1) > 0
           enddo
          endif ! diff_rho(1) > 0
        enddo
       endif ! watercaptag
      enddo
!=======================================================================
      RETURN
#endif
      END

   SUBROUTINE find_layering_icetable(porefill,psat_soil,psat_surf,pwat_surf,psat_bottom,B,index_IS,depth_filling,index_filling,index_geothermal,depth_geothermal,dz_etadz_rho)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!
!!! Purpose: Compute layering between dry soil, pore filling ice, and ice sheet based on Schorgofer, Icarus (2010). Adapted from NS MSIM 
!!! 
!!! Author: LL
!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
use comsoil_h_PEM,only: nsoilmx_PEM,mlayer_PEM 
implicit none
! inputs
! ------

real,intent(in) :: porefill(nsoilmx_PEM) ! Fraction of pore space filled with ice [Unitless] 0 <= f <= 1 for pore ice
real,intent(in) :: psat_soil(nsoilmx_PEM) ! Soil water pressure at saturation, yearly averaged [Pa]
real,intent(in) :: psat_surf ! surface water pressure at saturation, yearly averaged [Pa]
real,intent(in) :: pwat_surf ! Water vapor pressure  at the surface, not necesseraly at saturation, yearly averaged [Pa]
real,intent(in) :: psat_bottom ! Boundary conditions for soil vapor pressure [Pa]
real,intent(in) :: B ! constant (Eq. 8 from  Schorgofer, Icarus (2010).) 
integer, intent(in) :: index_IS ! index of the soil layer where the ice sheet begins [1]
real, intent(inout) :: depth_filling ! depth where pore filling begins [m]

! outputs
! -------
integer,intent(out) :: index_filling ! index where the pore filling begins [1]
integer, intent(out) :: index_geothermal ! index where the ice table stops [1]
real, intent(out) :: depth_geothermal ! depth where the ice table stops [m]
real, intent(out) :: dz_etadz_rho(nsoilmx_PEM) ! \partial z(eta \partial z rho), eta is the constriction, used later for pore filling increase

! local

real :: eta(nsoilmx_PEM)  ! constriction
integer :: ilay ! index for loop
real :: old_depth_filling ! depth_filling saved
real :: dz_psat(nsoilmx_PEM) ! first derivative of the vapor pressure at saturationn
integer :: index_tmp ! for loop 
real :: Jdry ! flux trought the dry layer
real :: Jsat ! flux trought the ice layer
real :: Jdry_prevlay,Jsat_prevlay ! same but for the previous ice layer
integer :: index_firstice ! first index where ice appears (i.e., f > 0)
real :: dz_eta(nsoilmx_PEM) ! \partial z \eta
real :: dz_eta_low ! same but evaluated at the interface for ice
real :: dzz_psat(nsoilmx_PEM) ! \partial \partial psat
real :: massfillabove,massfillafter ! h2O mass above and after index_geothermal

! constant
real :: pvap2rho = 18.e-3/8.314
!




  
! 0. Compute constriction over the layer
! Within the ice sheet, constriction is set to 0. Elsewhere, constriction =  (1-porefilling)**2
if (index_IS.lt.0) then
   index_tmp = nsoilmx_PEM
   do ilay = 1,nsoilmx_PEM
       call constriction(porefill(ilay),eta(ilay))
   enddo
else
   index_tmp = index_IS
   do ilay = 1,index_IS-1
       call constriction(porefill(ilay),eta(ilay))
   enddo
   do ilay = index_IS,nsoilmx_PEM
       eta(ilay) = 0.
   enddo
endif

! 1. Depth at which pore filling occurs. We solve Eq. 9 from  Schorgofer, Icarus  (2010)    

old_depth_filling = depth_filling

call deriv1(mlayer_PEM,nsoilmx_PEM,psat_soil,psat_surf,psat_bottom,dz_psat)  

do ilay = 1,index_tmp
  Jdry = (psat_soil(ilay) - pwat_surf)/mlayer_PEM(ilay) ! left member of Eq. 9 from Schorgofer, Icarus  (2010)
  Jsat =  eta(ilay)*dz_psat(ilay) !right member of Eq. 9 from Schorgofer, Icarus (2010)
  if((Jdry - Jsat).le.0) then
     index_filling = ilay
     exit
   endif
enddo

if(index_filling.eq.1)  depth_filling = mlayer_PEM(1) 

if(index_filling.gt.1) then
    Jdry_prevlay = (psat_soil(index_filling-1) - pwat_surf)/mlayer_PEM(index_filling-1)
    Jsat_prevlay = eta(index_filling-1)*dz_psat(index_filling-1)
   call findroot(Jdry-Jsat,Jdry_prevlay-Jsat_prevlay,mlayer_PEM(index_filling),mlayer_PEM(index_filling-1),depth_filling)
endif


! 2. Compute d_z (eta* d_z(rho)) (last term in Eq. 13 of Schorgofer, Icarus (2010))


! 2.0 preliminary: depth to shallowest  ice (discontinuity at interface)

index_firstice = -1
do ilay = 1,nsoilmx_PEM
   if (porefill(ilay).le.0.) then
        index_firstice = ilay  ! first point with ice
        exit
    endif
enddo

! 2.1: now we can computeCompute d_z (eta* d_z(rho))

call deriv1(mlayer_PEM,nsoilmx_PEM,eta,1.,eta(nsoilmx_PEM-1),dz_eta)

if ((index_firstice.gt.0).and.(index_firstice.lt.nsoilmx_PEM-2)) then
     call deriv1_onesided(index_firstice,mlayer_PEM,nsoilmx_PEM,eta,dz_eta(index_firstice))
endif

call deriv2_simple(mlayer_PEM,nsoilmx_PEM,psat_soil,psat_surf,psat_bottom,dzz_psat)
dz_etadz_rho(:) = pvap2rho*(dz_eta(:)*dz_psat(:) + eta(:)*dzz_psat(:))

! 3. Ice table boundary due to geothermal heating

if(index_IS.gt.0) index_geothermal = -1
if(index_geothermal.lt.0) depth_geothermal = -1.
if((index_geothermal.gt.0).and.(index_IS.lt.0)) then ! Eq. 21 from Schorfoger, Icarus (2010)
   index_geothermal = -1
   do ilay=2,nsoilmx_PEM
        if (dz_psat(ilay).gt.0.) then  ! first point with reversed flux
           index_geothermal=ilay
           call findroot(dz_psat(ilay-1),dz_psat(ilay),mlayer_PEM(ilay-1),mlayer_PEM(ilay),depth_geothermal)
           exit
        endif
     enddo
  else
     index_geothermal = -1
  endif
 if ((index_geothermal.gt.0).and.(index_IS.lt.0)) then !  Eq. 24 from Schorgofer, Icarus (2010)
     call  colint(porefill(:)/eta(:),mlayer_PEM,nsoilmx_PEM,index_geothermal-1,nsoilmx_PEM,massfillabove) 
     index_tmp = -1
     do ilay=index_geothermal,nsoilmx_PEM
        if (minval(eta(ilay:nsoilmx_PEM)).le.0.) cycle  ! eta=0 means completely full
         call colint(porefill(:)/eta(:),mlayer_PEM,nsoilmx_PEM,ilay,nsoilmx_PEM,massfillafter)
        if (massfillafter<dz_psat(ilay)*pvap2rho*B) then  ! usually executes on i=typeG
           if (ilay>index_geothermal) then
!              write(34,*) '# adjustment to geotherm depth by',ilay-index_geothermal
              call findroot(dz_psat(ilay-1)*pvap2rho*B-massfillabove, &  
                   dz_psat(ilay)*pvap2rho*B-massfillafter,mlayer_PEM(ilay-1),mlayer_PEM(ilay),depth_geothermal)
!               if (depth_geothermal.gt.mlayer_PEM(ilay) .or. depth_geothermal.lt.<mlayer_PEM(ilay-1)) write(34,*) 
!                     '# WARNING: zdepthG interpolation failed',ilay,mlayer_PEM(ilay-1),depth_geothermal,mlayer_PEM(ilay)
              index_tmp=ilay
           endif
           ! otherwise leave depth_geothermal unchanged
           exit
        endif
        massfillabove = massfillafter
     enddo
     if (index_tmp.gt.0) index_geothermal = index_tmp
  end if
  return
end subroutine 





SUBROUTINE findroot(y1,y2,z1,z2,zr)
        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        !!!
        !!! Purpose: Compute the root zr, between two values y1 and y2 at depth z1,z2
        !!! 
        !!! Author: LL
        !!!
        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

        implicit none
        real,intent(in) :: y1,y2 ! difference between surface water density and at depth  [kg/m^3]
        real,intent(in) :: z1,z2 ! depth [m]
        real,intent(out) :: zr ! depth at which we have zero
        zr = (y1*z2 - y2*z1)/(y1-y2)
        RETURN
        end 

SUBROUTINE constriction(porefill,eta)

        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        !!!
        !!! Purpose: Compute the  constriction of vapor flux by pore ice 
        !!! 
        !!! Author: LL
        !!!
        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        implicit none
        real,intent(in) :: porefill ! pore filling fraction
        real,intent(out) :: eta ! constriction

        !!!


        if (porefill.le.0.) eta = 1.
        if ((porefill.gt.0.) .and.(porefill.lt.1.)) then
        eta = (1-porefill)**2  ! Hudson et al., JGR, 2009
        endif
        if (porefill.le.1.) eta = 0.
        return
        end



subroutine deriv1(z,nz,y,y0,ybot,dzY)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!
!!! Purpose: Compute the first derivative of a function y(z) on an irregular grid 
!!!           
!!! Author: From N.S (N.S, Icarus 2010), impletented here by LL
!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        implicit none

! first derivative of a function y(z) on irregular grid
! upper boundary conditions: y(0)=y0
! lower boundary condition.: yp = ybottom
  integer, intent(IN) :: nz ! number of layer
  real, intent(IN) :: z(nz) ! depth layer
  real, intent(IN) :: y(nz) ! function which needs to be derived
  real, intent(IN) :: y0,ybot ! boundary conditions
  real, intent(OUT) :: dzY(nz) ! derivative of y w.r.t depth
! local
  integer :: j
  real :: hm,hp,c1,c2,c3

  hp = z(2)-z(1)
  c1 = z(1)/(hp*z(2))
  c2 = 1/z(1) - 1/(z(2)-z(1))
  c3 = -hp/(z(1)*z(2))
  dzY(1) = c1*y(2) + c2*y(1) + c3*y0
  do j=2,nz-1
     hp = z(j+1)-z(j)
     hm = z(j)-z(j-1)
     c1 = +hm/(hp*(z(j+1)-z(j-1)))
     c2 = 1/hm - 1/hp
     c3 = -hp/(hm*(z(j+1)-z(j-1)))
     dzY(j) = c1*y(j+1) + c2*y(j) + c3*y(j-1)
  enddo
  dzY(nz) = (ybot - y(nz-1))/(2.*(z(nz)-z(nz-1)))
 return
end subroutine deriv1



subroutine deriv2_simple(z,nz,y,y0,yNp1,yp2)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!
!!! Purpose: Compute the second derivative of a function y(z) on an irregular grid 
!!!           
!!! Author: N.S (raw copy/paste from MSIM)
!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  ! second derivative y_zz on irregular grid
  ! boundary conditions: y(0)=y0, y(nz+1)=yNp1
  implicit none
  integer, intent(IN) :: nz
  real, intent(IN) :: z(nz),y(nz),y0,yNp1
  real, intent(OUT) :: yp2(nz)
  integer j
  real hm,hp,c1,c2,c3

  c1 = +2./((z(2)-z(1))*z(2))
  c2 = -2./((z(2)-z(1))*z(1))
  c3 = +2./(z(1)*z(2))
  yp2(1) = c1*y(2) + c2*y(1) + c3*y0

  do j=2,nz-1
     hp = z(j+1)-z(j)
     hm = z(j)-z(j-1)
     c1 = +2./(hp*(z(j+1)-z(j-1)))
     c2 = -2./(hp*hm)
     c3 = +2./(hm*(z(j+1)-z(j-1)))
     yp2(j) = c1*y(j+1) + c2*y(j) + c3*y(j-1)
  enddo
  yp2(nz) = (yNp1 - 2*y(nz) + y(nz-1))/(z(nz)-z(nz-1))**2
  return
end subroutine deriv2_simple



subroutine  deriv1_onesided(j,z,nz,y,dy_zj)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!
!!! Purpose:   First derivative of function y(z) at z(j)  one-sided derivative on irregular grid
!!!           
!!! Author: N.S (raw copy/paste from MSIM)
!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  implicit none
  integer, intent(IN) :: nz,j
  real, intent(IN) :: z(nz),y(nz)
  real, intent(out) :: dy_zj
  real h1,h2,c1,c2,c3
  if (j<1 .or. j>nz-2) then
     dy_zj = -1.
  else
     h1 = z(j+1)-z(j)
     h2 = z(j+2)-z(j+1)
     c1 = -(2*h1+h2)/(h1*(h1+h2))
     c2 =  (h1+h2)/(h1*h2)
     c3 = -h1/(h2*(h1+h2))
     dy_zj = c1*y(j) + c2*y(j+1) + c3*y(j+2)
  endif
  return
end subroutine deriv1_onesided


subroutine  colint(y,z,nz,i1,i2,integral)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!
!!! Purpose:  Column integrates y on irregular grid
!!!           
!!! Author: N.S (raw copy/paste from MSIM)
!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  
  implicit none
  integer, intent(IN) :: nz, i1, i2
  real, intent(IN) :: y(nz), z(nz)
  real,intent(out) :: integral
  integer i
  real dz(nz)
  
  dz(1) = (z(2)-0.)/2
  do i=2,nz-1
     dz(i) = (z(i+1)-z(i-1))/2.
  enddo
  dz(nz) = z(nz)-z(nz-1)
  integral = sum(y(i1:i2)*dz(i1:i2))
end subroutine colint






        end module