MODULE waterice_tifeedback_mod
 
!======================================================================================================================!
! Subject:
!---------
!   Module used to compute the thermal inertia of an icy soil (either pore filling or massive ice)
!----------------------------------------------------------------------------------------------------------------------!
! Reference:
!-----------
!  For pure ice on the surface: J.-B. Madeleine, F. Forget, James W. Head, B. Levrard, F. Montmessin, E. Millour,
!  Amazonian northern mid-latitude glaciation on Mars: A proposed climate scenario, Icarus (10.1016/j.icarus.2009.04.037).
!
!  For pore-filling ice: Siegler, M., O. Aharonson, E. Carey, M. Choukroun, T. Hudson, N. Schorghofer, and S. Xu (2012), Measurements of thermal properties of icy Mars regolith analogs, JGR (10.1029/2011JE003938). 
!
! Originally written by JBM for pure ice (2008-2012), moved in to a module and .F90 by LL (2024). Pore filling ice included by LL (2024)
!
!======================================================================================================================!



IMPLICIT NONE

CONTAINS

SUBROUTINE waterice_tifeedback(ngrid,nsoil,nslope,icecover,poreice,newtherm_i)

      use tracer_mod, only: rho_ice
      use comsoil_h, only: layer, inertiedat, porosity_reg
      use surfdat_h, only: watercaptag, inert_h2o_ice
      IMPLICIT NONE
      include "callkeys.h"
!=======================================================================
!   Description :
!       Surface water ice / pore filling ice thermal inertia feedback.
!
!   When surface water-ice is thick enough (flag surfaceice_tifeedback), this routine creates a new
!   soil thermal inertia with three different layers :
!   - One layer of surface water ice (the thickness is given
!     by the variable icecover (in kg of ice per m2) and the thermal
!     inertia is prescribed by inert_h2o_ice (see surfdat_h)); 
!   - A transitional layer of mixed thermal inertia;
!   - A last layer of regolith below the ice cover whose thermal inertia
!     is equal to inertiedat.
!
!   To use the model :
!       SET THE surfaceice_tifeedback LOGICAL TO ".true." in callphys.def.
!
!   When pore filling ice is present, surface ice is computed as in Siegler et al., 2012
!   \sqrt(surf_thermalinertia**2 + porosity*pore_filling*inertie_purewaterice**2) 
!
!   For now, the code can not run with both options
!
!
!   Author: J.-B. Madeleine Mars 2008 - Updated November 2012; Added porous ice by LL February 2024
!=======================================================================

! Local variables
! ---------------

      INTEGER :: ig                       ! Grid point (ngrid)
      INTEGER :: islope                   ! SubGrid point (nslope)
      INTEGER :: ik                       ! Grid point (nsoil)
      INTEGER :: iref                     ! Ice/Regolith boundary index
      REAL :: icedepth                    ! Ice cover thickness (m)
      REAL :: inertie_purewaterice = 2100 ! Thermal inertia of pure water ice (J/m^2/K/s^1/2)

! Inputs
! ------
      INTEGER :: ngrid                 ! Number of horizontal grid points 
      INTEGER :: nsoil                 ! Number of soil layers 
      INTEGER :: nslope                ! Number of subgrid slopes 
      REAL icecover(ngrid,nslope)      ! water iceon the surface (kg.m-2)
      REAL poreice(ngrid,nsoil,nslope) ! pore ice filling fraction (1)
! Outputs
! -------

      REAL newtherm_i(ngrid,nsoil,nslope)  ! New soil thermal inertia

! Initialization
! --------------

  newtherm_i(1:ngrid,1:nsoil,1:nslope) = 0
  
  IF (surfaceice_tifeedback) THEN
      
! Creating the new soil thermal inertia table because of the massive surface ice
! ------------------------------------------------------------------------------
     DO islope = 1,nslope
        DO ig=1,ngrid
!       Calculating the ice cover thickness
           icedepth=icecover(ig,islope)/rho_ice
!       If the ice cover is too thick or watercaptag=true,
!       the entire column is changed :
           IF ((icedepth.ge.layer(nsoil)).or.(watercaptag(ig))) THEN
              DO ik=1,nsoil
                 newtherm_i(ig,ik,islope)=inert_h2o_ice
              ENDDO 
!          We neglect the effect of a very thin ice cover :
           ELSE IF (icedepth.lt.layer(1)) THEN
              DO ik=1,nsoil
                 newtherm_i(ig,ik,islope)=inertiedat(ig,ik)
              ENDDO 
           ELSE
!          Ice/regolith boundary index :
              iref=1
!          Otherwise, we find the ice/regolith boundary:
              DO ik=1,nsoil-1
                 IF ((icedepth.ge.layer(ik)).and.icedepth.lt.layer(ik+1)) THEN
                    iref=ik+1
                 EXIT
                 ENDIF
              ENDDO
!             And we change the thermal inertia:
              DO ik=1,iref-1
                 newtherm_i(ig,ik,islope)=inert_h2o_ice
              ENDDO
!             Transition (based on the equations of thermal conduction):
              newtherm_i(ig,iref,islope)=sqrt( (layer(iref)-layer(iref-1))/(((icedepth-layer(iref-1))/inert_h2o_ice**2) + &
              ((layer(iref)-icedepth)/inertiedat(ig,ik)**2) ) )
!              Underlying regolith:
              DO ik=iref+1,nsoil
                 newtherm_i(ig,ik,islope)=inertiedat(ig,ik)
              ENDDO
           ENDIF ! icedepth
        ENDDO ! ig
     ENDDO ! islope
     
  ELSE IF (poreice_tifeedback) THEN
  
! Updating soil thermal properties to includes the effect of porous ice
! ------------------------------------------------------------------------------ 
     DO islope = 1,nslope
        newtherm_i(:,:,islope) = sqrt(inertiedat(:,:)**2 + porosity_reg*poreice(:,:,islope)*inertie_purewaterice**2)
     ENDDO
  ENDIF

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

END SUBROUTINE waterice_tifeedback
END MODULE waterice_tifeedback_mod