source: trunk/LMDZ.COMMON/libf/evolution/soil_thermalproperties_mod.F90 @ 3210

MODULE soil_thermalproperties_mod

implicit none

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!
!!! Purpose: Compute the soil thermal properties 
!!! Author:  LL, 04/2023
!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!=======================================================================
contains
!=======================================================================

SUBROUTINE ice_thermal_properties(ispureice,pore_filling,surf_thermalinertia,ice_thermalinertia)
!=======================================================================
!   subject: Compute ice thermal properties
!   --------
!
!   author: LL, 04/2023 
!   -------
!         
!=======================================================================

use constants_marspem_mod, only: porosity

implicit none

!-----------------------------------------------------------------------
!=======================================================================
!    Declarations :
!=======================================================================
!
!    Input/Output
!    ------------
logical, intent(in) :: ispureice           ! Boolean to check if ice is massive or just pore filling
real,    intent(in) :: pore_filling        ! ice pore filling in each layer (1)
real,    intent(in) :: surf_thermalinertia ! surface thermal inertia (J/m^2/K/s^1/2)  
real,    intent(out) :: ice_thermalinertia ! Thermal inertia of ice when present in the pore (J/m^2/K/s^1/2)

!    Local Variables
!    --------------
REAL :: inertie_purewaterice = 2100 ! 2050 is better according to my computations with the formula from Siegler et al., 2012, but let's follow Mellon et al. (2004) 
!=======================================================================
!   Beginning of the code
!=======================================================================

if (ispureice) then
    ice_thermalinertia = inertie_purewaterice
else
    ice_thermalinertia = sqrt(surf_thermalinertia**2 + porosity*pore_filling*inertie_purewaterice**2)  ! Siegler et al., 2012
endif

END SUBROUTINE 
!=======================================================================

SUBROUTINE update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,tendencies_waterice,waterice,p_avg_new,ice_depth,ice_thickness,TI_PEM)

use comsoil_h,             only: volcapa
use comsoil_h_PEM,         only: layer_PEM, inertiedat_PEM, depth_breccia, depth_bedrock, index_breccia, index_bedrock, reg_thprop_dependp
use constants_marspem_mod, only: TI_breccia, TI_bedrock, TI_regolith_avg

implicit none

! Input: 
integer, intent(in) ::  ngrid, nslope, nsoil_PEM      ! Shape of the arrays: physical grid, number of sub-grid slopes, number of layer in the soil
real, intent(in) :: p_avg_new                         ! Global average surface pressure [Pa]
real, intent(in) :: tendencies_waterice(ngrid,nslope) ! Tendencies on the water ice [kg/m^2/year]
real, intent(in) :: waterice(ngrid,nslope)            ! Surface Water ice [kg/m^2]
real, intent(in) :: ice_depth(ngrid,nslope)           ! Ice table depth [m]
real, intent(in) :: ice_thickness(ngrid,nslope)       ! Ice table thickness [m]
! Outputs:

real, intent(inout) :: TI_PEM(ngrid,nsoil_PEM,nslope) ! Soil Thermal Inertia [J/m^2/K/s^1/2]
  
! Constants:

 REAL ::  reg_inertie_thresold = 370.                      ! Above this thermal inertia, the regolith has too much cementation to be dependant on the pressure [J/m^2/K/s^1/2]
 REAL ::  reg_inertie_minvalue = 50.                       ! Minimum value of the Thermal Inertia at low pressure (Piqueux & Christensen 2009) [J/m^2/K/s^1/2]
 REAL ::  reg_inertie_maxvalue = 370.                      ! Maximum value of the Thermal Inertia at low pressure (Piqueux & Christensen 2009) [J/m^2/K/s^1/2]
 REAL ::  ice_inertia                                      ! Inertia of water ice [SI]
 REAL ::  P610 = 610.0                                     ! current average pressure on Mars [Pa]
 REAL ::  C = 0.0015                                       ! Constant to derive TI as a function of P, from Presley and Christensen 1997 [uniteless]
 REAL ::  K = 8.1*1e4                                      ! Constant to derive TI as a function of P, from Presley and Christensen 1997 [tor, or 133.3Pa]
 REAL ::  Pa2Tor = 1./133.3                                ! Conversion from Pa to tor [Pa/tor]


! Local variables:

 INTEGER :: ig,islope,isoil,iref,iend                      ! Loop variables
 REAL :: regolith_inertia(ngrid,nslope)                    ! Thermal inertia of the regolith (first layer in the GCM) [J/m^2/K/s^1/2]
 REAL :: delta                                             ! Difference of depth to compute the  thermal inertia in a mixed layer [m]
 REAL :: ice_bottom_depth                                  ! Bottom depth of the subsurface ice [m]
 REAL :: d_part                                            ! Regolith particle size [micrometer]
 

write(*,*) "Update soil propreties"
        
! 1. Modification of the regolith thermal inertia.
 do islope = 1,nslope
   regolith_inertia(:,islope) = inertiedat_PEM(:,1)
   do ig = 1,ngrid
      if((tendencies_waterice(ig,islope).lt.-1e-5).and.(waterice(ig,islope).eq.0)) then
              regolith_inertia(ig,islope) = TI_regolith_avg
       endif
      if(reg_thprop_dependp) then
          if(TI_PEM(ig,1,islope).lt.reg_inertie_thresold) then
            d_part = (regolith_inertia(ig,islope)**2/(volcapa*C*(P610*Pa2Tor)**(0.6)))**(-1./(0.11*log10(P610*Pa2Tor/K))) ! compute particle size, in micrometer
            regolith_inertia(ig,islope) = sqrt(volcapa*C*(p_avg_new*Pa2Tor)**(0.6)*d_part**(-0.11*log10(p_avg_new*Pa2Tor/K)))
            if(regolith_inertia(ig,islope).gt.reg_inertie_maxvalue) regolith_inertia(ig,islope) = reg_inertie_maxvalue
            if(regolith_inertia(ig,islope).lt.reg_inertie_minvalue) regolith_inertia(ig,islope) = reg_inertie_minvalue
          endif 
      endif
   enddo
 enddo

! 2. Build new Thermal Inertia 
do  islope=1,nslope
   do ig = 1,ngrid
     do isoil = 1,index_breccia
         TI_PEM(ig,isoil,islope) = regolith_inertia(ig,islope)
     enddo
     if(regolith_inertia(ig,islope).lt.TI_breccia) then 
!!! transition
             delta = depth_breccia
             TI_PEM(ig,index_breccia+1,islope) = sqrt((layer_PEM(index_breccia+1)-layer_PEM(index_breccia))/ &
            (((delta-layer_PEM(index_breccia))/(TI_PEM(ig,index_breccia,islope)**2))+ &
                      ((layer_PEM(index_breccia+1)-delta)/(TI_breccia**2))))
            do isoil=index_breccia+2,index_bedrock
              TI_PEM(ig,isoil,islope) = TI_breccia
            enddo      
      else ! we keep the high ti values
            do isoil=index_breccia+1,index_bedrock
              TI_PEM(ig,isoil,islope) = TI_PEM(ig,index_breccia,islope)
            enddo 
       endif ! TI PEM and breccia comparison
!! transition
       delta = depth_bedrock
       TI_PEM(ig,index_bedrock+1,islope) = sqrt((layer_PEM(index_bedrock+1)-layer_PEM(index_bedrock))/ &
            (((delta-layer_PEM(index_bedrock))/(TI_PEM(ig,index_bedrock,islope)**2))+ &
            ((layer_PEM(index_bedrock+1)-delta)/(TI_bedrock**2))))
       do isoil=index_bedrock+2,nsoil_PEM
            TI_PEM(ig,isoil,islope) = TI_bedrock
       enddo   
      enddo ! ig
ENDDO ! islope

!  3. Build new TI in case of ice table
  do ig=1,ngrid
    do islope=1,nslope
     if (ice_depth(ig,islope).gt.-1e-6) then
! 3.0 Case where it is perennial ice
      if (ice_depth(ig,islope).lt. 1e-10) then
       call ice_thermal_properties(.true.,1.,0.,ice_inertia)
       do isoil = 1,nsoil_PEM
         TI_PEM(ig,isoil,islope)=ice_inertia
       enddo
      else 

        call ice_thermal_properties(.false.,1.,regolith_inertia(ig,islope),ice_inertia)
        ! 3.1.1 find the index of the mixed layer
        iref=0 ! initialize iref
        do isoil=1,nsoil_PEM ! loop on layers to find the beginning of the ice table
          if (ice_depth(ig,islope).ge.layer_PEM(isoil)) then
            iref=isoil ! pure regolith layer up to here
          else
         ! correct iref was obtained in previous cycle
            exit
          endif
        enddo
        ! 3.1.2 find the index of the end of the ice table
        iend=0 ! initialize iend
        ice_bottom_depth = ice_depth(ig,islope)+ice_thickness(ig,islope)
        do isoil=1,nsoil_PEM ! loop on layers to find the end of the ice table
          if (ice_bottom_depth.ge.layer_PEM(isoil)) then
            iend=isoil ! pure regolith layer up to here
          else
         ! correct iref was obtained in previous cycle
            exit
          endif
        enddo
      ! 3.2 Build the new ti
        do isoil=1,iref
          TI_PEM(ig,isoil,islope) = TI_PEM(ig,1,islope)
        enddo
        if (iref.lt.nsoil_PEM) then
           if(iref.eq.iend) then
            ! Ice table begins and end in the same mixture Mixtures with three components 
            if (iref.ne.0) then ! 
            ! mixed layer
              TI_PEM(ig,iref+1,islope)=sqrt((layer_PEM(iref+1)-layer_PEM(iref))/ &
                      (((ice_depth(ig,islope)-layer_PEM(iref))/(TI_PEM(ig,iref,islope)**2))+ &
                      ((ice_bottom_depth-ice_depth(ig,islope))/(ice_inertia**2))+ &
                      ((layer_PEM(iref+1)-ice_bottom_depth)/(TI_PEM(ig,iref+1,islope)**2))))
            else ! first layer is already a mixed layer
              ! (ie: take layer(iref=0)=0)
              TI_PEM(ig,1,islope)=sqrt((layer_PEM(1))/ &
                      (((ice_depth(ig,islope))/(TI_PEM(ig,1,islope)**2))+ &
                      ((ice_bottom_depth-ice_depth(ig,islope))/(ice_inertia**2))+ &
                      ((layer_PEM(2)-ice_bottom_depth)/(TI_PEM(ig,2,islope)**2))))
            endif ! of if (iref.ne.0) 
           else
            if (iref.ne.0) then
            ! mixed layer
              TI_PEM(ig,iref+1,islope)=sqrt((layer_PEM(iref+1)-layer_PEM(iref))/ &
                      (((ice_depth(ig,islope)-layer_PEM(iref))/(TI_PEM(ig,iref,islope)**2))+ &
                      ((layer_PEM(iref+1)-ice_depth(ig,islope))/(ice_inertia**2))))
            else ! first layer is already a mixed layer
              ! (ie: take layer(iref=0)=0)
              TI_PEM(ig,1,islope)=sqrt((layer_PEM(1))/ &
                      (((ice_depth(ig,islope))/(TI_PEM(ig,1,islope)**2))+ &
                      ((layer_PEM(1)-ice_depth(ig,islope))/(ice_inertia**2))))
            endif ! of if (iref.ne.0)        
           endif ! iref.eq.iend
      ! 3.3 Build the new ti
          do isoil=iref+2,iend
            TI_PEM(ig,isoil,islope)=ice_inertia
          enddo


          if(iend.lt.nsoil_PEM) then
             TI_PEM(ig,iend+1,islope)=sqrt((layer_PEM(iend+1)-layer_PEM(iend))/ &
                      (((ice_bottom_depth-layer_PEM(iend))/(TI_PEM(ig,iend,islope)**2)) + &
                      ((layer_PEM(iend+1)-ice_bottom_depth)/(ice_inertia**2))))
          endif


        endif ! of if (iref.lt.(nsoilmx))
      endif ! permanent glaciers
     endif !ice_depth(ig,islope) > 0.
!     write(*,*) 'TI=', TI_PEM(ig,:,islope)
    enddo !islope
   enddo !ig

END SUBROUTINE update_soil_thermalproperties

END MODULE soil_thermalproperties_mod