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: layer_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) = -10.

         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
               z1 = (diff_rho(isoil) - diff_rho(isoil+1))/(layer_PEM(isoil) - layer_PEM(isoil+1))
               z2 = -layer_PEM(isoil+1)*z1 +  diff_rho(isoil+1)
               ice_table(ig,islope) = -z2/z1
               exit
             endif !diff_rho(z) <0 & diff_rho(z+1) > 0
           enddo
          endif ! diff_rho(1) > 0
        enddo
       endif ! watercaptag
      enddo
!=======================================================================
      RETURN
#endif
      END



end module