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) #ifndef CPP_STD USE comsoil_h, only: inertiedat, volcapa USE comsoil_h_PEM, only: layer_PEM,inertiedat_PEM,depth_breccia,depth_bedrock,index_breccia,index_bedrock,reg_thprop_dependp USE vertical_layers_mod, ONLY: ap,bp 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 :: inertie_thresold = 800. ! Above this thermal inertia, it's ice [SI] REAL :: ice_inertia ! Inertia of water ice [SI] REAL :: P610 = 610.0 ! current average pressure on Mars [Pa] ! Local variables: INTEGER :: ig,islope,iloop,iref,k,iend REAL :: regolith_inertia(ngrid,nslope) ! TI of the regolith REAL :: delta REAL :: TI_breccia_new REAL :: ice_bottom_depth ! 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.inertie_thresold) then regolith_inertia(ig,islope) = regolith_inertia(ig,islope)*(p_avg_new/P610)**0.3 endif TI_breccia_new = TI_breccia*(p_avg_new/P610)**0.3 else TI_breccia_new = TI_breccia endif enddo enddo ! 2. Build new Thermal Inertia do islope=1,nslope do ig = 1,ngrid do iloop = 1,index_breccia TI_PEM(ig,iloop,islope) = regolith_inertia(ig,islope) enddo if(regolith_inertia(ig,islope).lt.TI_breccia_new) 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_new**2)))) do iloop=index_breccia+2,index_bedrock TI_PEM(ig,iloop,islope) = TI_breccia_new enddo else ! we keep the high ti values do iloop=index_breccia+1,index_bedrock TI_PEM(ig,iloop,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 iloop=index_bedrock+2,nsoil_PEM TI_PEM(ig,iloop,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 perenial ice if (ice_depth(ig,islope).lt. 1e-10) then call ice_thermal_properties(.true.,1.,0.,ice_inertia) do iloop = 1,nsoil_PEM TI_PEM(ig,iloop,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 k=1,nsoil_PEM ! loop on layers to find the beginning of the ice table if (ice_depth(ig,islope).ge.layer_PEM(k)) then iref=k ! 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 k=1,nsoil_PEM ! loop on layers to find the end of the ice table if (ice_bottom_depth.ge.layer_PEM(k)) then iend=k ! pure regolith layer up to here else ! correct iref was obtained in previous cycle exit endif enddo ! 3.2 Build the new ti do iloop=1,iref TI_PEM(ig,iloop,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 iloop=iref+2,iend TI_PEM(ig,iloop,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 !======================================================================= RETURN #endif END subroutine update_soil_thermalproperties end module soil_thermalproperties_mod