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