MODULE soil !----------------------------------------------------------------------- ! NAME ! soil ! ! DESCRIPTION ! Soil state and properties for PEM. ! ! AUTHORS & DATE ! L. Lange, 04/2023 ! JB Clement, 2023-2025 ! ! NOTES ! !----------------------------------------------------------------------- ! DEPENDENCIES ! ------------ use numerics, only: dp, di, k4 ! DECLARATION ! ----------- implicit none ! PARAMETERS ! ---------- logical(k4), protected :: do_soil ! To run with subsurface physics logical(k4), protected :: reg_thprop_dependp ! Regolith thermal properties depend on surface pressure real(dp), protected :: flux_geo ! Geothermal flux [W/m^2] real(dp), protected :: depth_breccia ! Depth of breccia [m] real(dp), protected :: depth_bedrock ! Depth of bedrock [m] real(dp), allocatable, dimension(:,:,:), protected :: TI_PCM ! Soil thermal inertia in the "startfi.nc" at the beginning [SI] ! VARIABLES ! --------- real(dp), allocatable, dimension(:) :: layer ! Soil layer depths [m] real(dp), allocatable, dimension(:) :: mlayer ! Soil mid-layer depths [m] real(dp), allocatable, dimension(:,:,:) :: TI ! Soil thermal inertia [SI] real(dp), allocatable, dimension(:,:) :: inertiedat ! Soil thermal inertia saved as reference for current climate [SI] real(dp), allocatable, dimension(:,:) :: inertiedat_PCM ! Soil thermal inertia saved as reference for current climate in the PCM [SI] real(dp), allocatable, dimension(:,:) :: mthermdiff ! Mid-layer thermal diffusivity [SI] real(dp), allocatable, dimension(:,:) :: thermdiff ! Inter-layer thermal diffusivity [SI] real(dp), allocatable, dimension(:) :: coefq ! q_{k+1/2} coefficients [SI] real(dp), allocatable, dimension(:,:) :: coefd ! d_k coefficients [SI] real(dp), allocatable, dimension(:,:) :: alph ! alpha_k coefficients [SI] real(dp), allocatable, dimension(:,:) :: beta ! beta_k coefficients [SI] real(dp) :: mu ! mu coefficient [SI] integer(di) :: index_breccia ! Last index of the depth grid before having breccia integer(di) :: index_bedrock ! Last index of the depth grid before having bedrock real(dp) :: volcapa ! Soil volumetric heat capacity contains !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ !======================================================================= SUBROUTINE ini_soil() !----------------------------------------------------------------------- ! NAME ! ini_soil ! ! DESCRIPTION ! Allocate soil arrays based on grid dimensions. ! ! AUTHORS & DATE ! L. Lange, 2023 ! JB Clement, 2023-2025 ! ! NOTES ! !----------------------------------------------------------------------- ! DEPENDENCIES ! ------------ use geometry, only: ngrid, nslope, nsoil, nsoil_PCM ! DECLARATION ! ----------- implicit none ! CODE ! ---- allocate(layer(nsoil)) allocate(mlayer(0:nsoil - 1)) allocate(TI(ngrid,nsoil,nslope)) allocate(mthermdiff(ngrid,0:nsoil - 1)) allocate(thermdiff(ngrid,nsoil - 1)) allocate(coefq(0:nsoil - 1)) allocate(coefd(ngrid,nsoil - 1)) allocate(alph(ngrid,nsoil - 1)) allocate(beta(ngrid,nsoil - 1)) allocate(inertiedat(ngrid,nsoil)) if (.not. allocated(inertiedat_PCM)) allocate(inertiedat_PCM(ngrid,nsoil_PCM)) if (.not. allocated(TI_PCM)) allocate(TI_PCM(ngrid,nsoil_PCM,nslope)) END SUBROUTINE ini_soil !======================================================================= !======================================================================= SUBROUTINE end_soil !----------------------------------------------------------------------- ! NAME ! end_soil ! ! DESCRIPTION ! Deallocate all soil arrays. ! ! AUTHORS & DATE ! L. Lange, 2023 ! JB Clement, 2023-2025 ! ! NOTES ! !----------------------------------------------------------------------- ! DECLARATION ! ----------- implicit none ! CODE ! ---- if (allocated(layer)) deallocate(layer) if (allocated(mlayer)) deallocate(mlayer) if (allocated(TI)) deallocate(TI) if (allocated(mthermdiff)) deallocate(mthermdiff) if (allocated(thermdiff)) deallocate(thermdiff) if (allocated(coefq)) deallocate(coefq) if (allocated(coefd)) deallocate(coefd) if (allocated(alph)) deallocate(alph) if (allocated(beta)) deallocate(beta) if (allocated(inertiedat)) deallocate(inertiedat) if (allocated(inertiedat_PCM)) deallocate(inertiedat_PCM) if (allocated(TI_PCM)) deallocate(TI_PCM) END SUBROUTINE end_soil !======================================================================= !======================================================================= SUBROUTINE set_soil_config(do_soil_in,reg_thprop_dependp_in,flux_geo_in,depth_breccia_in,depth_bedrock_in) !----------------------------------------------------------------------- ! NAME ! set_soil_config ! ! DESCRIPTION ! Setter for 'soil' configuration parameters. ! ! AUTHORS & DATE ! JB Clement, 02/2026 ! ! NOTES ! !----------------------------------------------------------------------- ! DEPENDENCIES ! ------------ use utility, only: real2str, bool2str use display, only: print_msg ! DECLARATION ! ----------- implicit none ! ARGUMENTS ! --------- logical(k4), intent(in) :: do_soil_in, reg_thprop_dependp_in real(dp), intent(in) :: flux_geo_in, depth_breccia_in, depth_bedrock_in ! CODE ! ---- do_soil = do_soil_in reg_thprop_dependp = reg_thprop_dependp_in flux_geo = flux_geo_in depth_breccia = depth_breccia_in depth_bedrock = depth_bedrock_in call print_msg('do_soil = '//bool2str(do_soil)) call print_msg('reg_thprop_dependp = '//bool2str(reg_thprop_dependp)) call print_msg('flux_geo = '//real2str(flux_geo)) call print_msg('depth_breccia = '//real2str(depth_breccia)) call print_msg('depth_bedrock = '//real2str(depth_bedrock)) END SUBROUTINE set_soil_config !======================================================================= !======================================================================= SUBROUTINE set_TI_PCM(TI_PCM_in) !----------------------------------------------------------------------- ! NAME ! set_TI_PCM ! ! DESCRIPTION ! Setter for 'TI_PCM'. ! ! AUTHORS & DATE ! JB Clement, 01/2026 ! ! NOTES ! !----------------------------------------------------------------------- ! DECLARATION ! ----------- implicit none ! ARGUMENTS ! --------- real(dp), dimension(:,:,:), intent(in) :: TI_PCM_in ! CODE ! ---- TI_PCM(:,:,:) = TI_PCM_in(:,:,:) END SUBROUTINE set_TI_PCM !======================================================================= !======================================================================= SUBROUTINE set_inertiedat_PCM(inertiedat_PCM_in) !----------------------------------------------------------------------- ! NAME ! set_inertiedat_PCM ! ! DESCRIPTION ! Setter for 'inertiedat_PCM'. ! ! AUTHORS & DATE ! JB Clement, 01/2026 ! ! NOTES ! !----------------------------------------------------------------------- ! DECLARATION ! ----------- implicit none ! ARGUMENTS ! --------- real(dp), dimension(:,:), intent(in) :: inertiedat_PCM_in ! CODE ! ---- inertiedat_PCM(:,:) = inertiedat_PCM_in(:,:) END SUBROUTINE set_inertiedat_PCM !======================================================================= !======================================================================= SUBROUTINE set_soil(TI) !----------------------------------------------------------------------- ! NAME ! set_soil ! ! DESCRIPTION ! Initialize soil depths and properties. Builds layer depths using ! exponential then power-law distribution; interpolates thermal inertia ! from PCM to PEM grid. ! ! AUTHORS & DATE ! E. Millour, 07/2006 ! L. Lange, 2023 ! JB Clement, 2023-2025 ! ! NOTES ! Modifications: ! Aug.2010 EM: use NetCDF90 to load variables (enables using ! r4 or r8 restarts independently of having compiled ! the GCM in r4 or r8) ! June 2013 TN: Possibility to read files with a time axis ! The various actions and variable read/initialized are: ! 1. Read/build layer (and midlayer) depth ! 2. Interpolate thermal inertia and temperature on the new grid, if necessary !----------------------------------------------------------------------- ! DEPENDENCIES ! ------------ use geometry, only: ngrid, nslope, nsoil, nsoil_PCM ! DECLARATION ! ----------- implicit none ! ARGUMENTS ! --------- real(dp), dimension(:,:,:), intent(inout) :: TI ! Thermal inertia in the PEM [SI] ! LOCAL VARIABLES ! --------------- integer(di) :: ig, iloop, islope ! loop counters real(dp) :: alpha, lay1 ! coefficients for building layers real(dp) :: index_powerlaw ! coefficient to match the power low grid with the exponential one ! CODE ! ---- ! 1. Depth coordinate ! ------------------- ! mlayer distribution: For the grid in common with the PCM: lay1*(1+k^(2.9)*(1-exp(-k/20)) ! Then we use a power low : lay1*alpha**(k-1/2) lay1 = 2.e-4 alpha = 2 do iloop = 0,nsoil_PCM - 1 mlayer(iloop) = lay1*(1 + iloop**2.9*(1._dp - exp(-real(iloop,dp)/20._dp))) end do do iloop = 0,nsoil - 1 if (lay1*(alpha**(iloop - 0.5_dp)) > mlayer(nsoil_PCM - 1)) then index_powerlaw = iloop exit end if end do do iloop = nsoil_PCM,nsoil - 1 mlayer(iloop) = lay1*(alpha**(index_powerlaw + (iloop - nsoil_PCM) - 0.5_dp)) end do ! Build layer() do iloop = 1,nsoil - 1 layer(iloop) = (mlayer(iloop) + mlayer(iloop - 1))/2._dp end do layer(nsoil) = 2._dp*mlayer(nsoil - 1) - mlayer(nsoil - 2) ! 2. Thermal inertia (note: it is declared in comsoil_h) ! ------------------ do ig = 1,ngrid do islope = 1,nslope do iloop = 1,nsoil_PCM TI(ig,iloop,islope) = TI_PCM(ig,iloop,islope) end do if (nsoil > nsoil_PCM) then do iloop = nsoil_PCM + 1,nsoil TI(ig,iloop,islope) = TI_PCM(ig,nsoil_PCM,islope) end do end if end do end do ! 3. Index for subsurface layering ! ------------------ index_breccia = 1 do iloop = 1,nsoil - 1 if (depth_breccia >= layer(iloop)) then index_breccia = iloop else exit end if end do index_bedrock = 1 do iloop = 1,nsoil - 1 if (depth_bedrock >= layer(iloop)) then index_bedrock = iloop else exit end if end do END SUBROUTINE set_soil !======================================================================= !======================================================================= SUBROUTINE build4PCM_soil(tsoil_avg,tsoil_dev,inertiesoil4PCM,tsoil4PCM,flux_geo4PCM) !----------------------------------------------------------------------- ! NAME ! build4PCM_soil ! ! DESCRIPTION ! Reconstructs the soil temperature profile, thermal inertia and ! the geothermal flux for the PCM. ! ! AUTHORS & DATE ! JB Clement, 12/2025 ! C. Metz, 01/2026 ! ! NOTES ! !----------------------------------------------------------------------- ! DEPENDENCIES ! ------------ use geometry, only: nsoil_PCM use display, only: print_msg ! DECLARATION ! ----------- implicit none ! ARGUMENTS ! --------- real(dp), dimension(:,:,:), intent(in) :: tsoil_avg, tsoil_dev real(dp), dimension(:,:,:), intent(out) :: tsoil4PCM, inertiesoil4PCM real(dp), dimension(:,:), intent(out) :: flux_geo4PCM ! CODE ! ---- call print_msg('> Building soil thermal inertia for the PCM') inertiesoil4PCM(:,:,:) = TI(:,1:nsoil_PCM,:) call print_msg('> Building soil temperature profile for the PCM') tsoil4PCM(:,:,:) = tsoil_avg(:,1:nsoil_PCM,:) + tsoil_dev(:,:,:) call print_msg('> Building geothermal flux for the PCM') flux_geo4PCM(:,:) = (tsoil_avg(:,nsoil_PCM + 1,:) - tsoil4PCM(:,nsoil_PCM,:))/(mlayer(nsoil_PCM + 1) - mlayer(nsoil_PCM)) END SUBROUTINE build4PCM_soil !======================================================================= END MODULE soil