MODULE soil_thermalproperties_mod implicit none !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!! !!! Purpose: Compute the soil thermal properties !!! Author: LL, 04/2023 !!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Constants: real, parameter :: 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, parameter :: reg_inertie_minvalue = 50. ! Minimum value of the Thermal Inertia at low pressure (Piqueux & Christensen 2009) [J/m^2/K/s^1/2] real, parameter :: reg_inertie_maxvalue = 370. ! Maximum value of the Thermal Inertia at low pressure (Piqueux & Christensen 2009) [J/m^2/K/s^1/2] real, parameter :: P610 = 610.0 ! current average pressure on Mars [Pa] real, parameter :: C = 0.0015 ! Constant to derive TI as a function of P, from Presley and Christensen 1997 [unitless] real, parameter :: K = 8.1*1e4 ! Constant to derive TI as a function of P, from Presley and Christensen 1997 [Torr, or 133.3Pa] real, parameter :: Pa2Torr = 1./133.3 ! Conversion from Pa to tor [Pa/Torr] !======================================================================= 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,tendencies_waterice,waterice,p_avg_new,icetable_depth,icetable_thickness,ice_porefilling,icetable_equilibrium,icetable_dynamic,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 ! 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, dimension(ngrid,nslope), intent(in) :: tendencies_waterice ! Tendencies on the water ice [kg/m^2/year] real, dimension(ngrid,nslope), intent(in) :: waterice ! Surface Water ice [kg/m^2] real, dimension(ngrid,nslope), intent(in) :: icetable_depth ! Ice table depth [m] real, dimension(ngrid,nslope), intent(in) :: icetable_thickness ! Ice table thickness [m] real, dimension(ngrid,nsoil,nslope), intent(in) :: ice_porefilling ! Ice porefilling [m^3/m^3] logical, intent(in) :: icetable_equilibrium, icetable_dynamic ! Computing method for ice table ! Outputs: real, dimension(ngrid,nsoil,nslope), intent(inout) :: TI_PEM ! Soil Thermal Inertia [J/m^2/K/s^1/2] ! Local variables: integer :: ig, islope, isoil, iref, iend ! Loop variables real, dimension(ngrid,nslope) :: regolith_inertia ! 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] real :: ice_inertia ! Inertia of water ice [SI] write(*,*) "Update soil properties" ! 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) < -1.e-5 .and. waterice(ig,islope) == 0) regolith_inertia(ig,islope) = TI_regolith_avg if (reg_thprop_dependp) then if (TI_PEM(ig,1,islope) < reg_inertie_thresold) then d_part = (regolith_inertia(ig,islope)**2/(volcapa*C*(P610*Pa2Torr)**(0.6)))**(-1./(0.11*log10(P610*Pa2Torr/K))) ! compute particle size, in micrometer regolith_inertia(ig,islope) = sqrt(volcapa*C*(p_avg_new*Pa2Torr)**(0.6)*d_part**(-0.11*log10(p_avg_new*Pa2Torr/K))) if (regolith_inertia(ig,islope) > reg_inertie_maxvalue) regolith_inertia(ig,islope) = reg_inertie_maxvalue if (regolith_inertia(ig,islope) < 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) < 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 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 (icetable_depth(ig,islope) > -1.e-6) then ! 3.0 Case where it is perennial ice if (icetable_depth(ig,islope) < 1.e-10) then call ice_thermal_properties(.true.,1.,0.,ice_inertia) do isoil = 1,nsoil TI_PEM(ig,isoil,islope) = ice_inertia enddo else if (icetable_equilibrium) then 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 ! loop on layers to find the beginning of the ice table if (icetable_depth(ig,islope) >= layer_PEM(isoil)) then iref = isoil ! pure regolith layer up to here else exit ! correct iref was obtained in previous cycle endif enddo ! 3.1.2 find the index of the end of the ice table iend = 0 ! initialize iend ice_bottom_depth = icetable_depth(ig,islope) + icetable_thickness(ig,islope) do isoil = 1,nsoil ! loop on layers to find the end of the ice table if (ice_bottom_depth >= layer_PEM(isoil)) then iend = isoil ! pure regolith layer up to here else exit ! correct iref was obtained in previous cycle endif enddo ! 3.2 Build the new ti if (iref < nsoil) then if (iref == iend) then ! Ice table begins and end in the same mixture with three components if (iref /= 0) then ! mixed layer TI_PEM(ig,iref + 1,islope) = sqrt((layer_PEM(iref + 1) - layer_PEM(iref))/ & (((icetable_depth(ig,islope) - layer_PEM(iref))/(TI_PEM(ig,iref,islope)**2)) + & ((ice_bottom_depth - icetable_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))/ & (((icetable_depth(ig,islope))/(TI_PEM(ig,1,islope)**2)) + & ((ice_bottom_depth - icetable_depth(ig,islope))/(ice_inertia**2)) + & ((layer_PEM(2) - ice_bottom_depth)/(TI_PEM(ig,2,islope)**2)))) endif ! of if (iref /= 0) else if (iref /= 0) then ! mixed layer TI_PEM(ig,iref + 1,islope) = sqrt((layer_PEM(iref + 1) - layer_PEM(iref))/ & (((icetable_depth(ig,islope) - layer_PEM(iref))/(TI_PEM(ig,iref,islope)**2)) + & ((layer_PEM(iref + 1) - icetable_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))/ & (((icetable_depth(ig,islope))/(TI_PEM(ig,1,islope)**2)) + & ((layer_PEM(1) - icetable_depth(ig,islope))/(ice_inertia**2)))) endif ! of if (iref /= 0) endif ! iref == iend TI_PEM(ig,iref + 2:iend,islope) = ice_inertia if (iend < nsoil) 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 < nsoilmx) else if (icetable_dynamic) then do isoil = 1,nsoil call ice_thermal_properties(.false.,ice_porefilling(ig,isoil,islope),regolith_inertia(ig,islope),TI_PEM(ig,isoil,islope)) enddo endif ! of if icetable_equilibrium endif ! permanent glaciers endif ! icetable_depth(ig,islope) > 0. enddo !islope enddo !ig END SUBROUTINE update_soil_thermalproperties END MODULE soil_thermalproperties_mod