MODULE soil_temp !----------------------------------------------------------------------- ! NAME ! soil_temp ! ! DESCRIPTION ! Routines to compute soil temperature evolution and initialization. ! ! AUTHORS & DATE ! L. Lange, 2023 ! JB Clement, 2023-2025 ! ! NOTES ! !----------------------------------------------------------------------- ! DECLARATION ! ----------- implicit none contains !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ !======================================================================= SUBROUTINE compute_tsoil(ngrid,nsoil,firstcall,therm_i,timestep,tsurf,tsoil) !----------------------------------------------------------------------- ! NAME ! compute_tsoil ! ! DESCRIPTION ! Compute soil temperature using an implicit 1st order scheme. ! ! AUTHORS & DATE ! L. Lange, 2023 ! JB Clement, 2023-2025 ! ! NOTES ! !----------------------------------------------------------------------- ! DEPENDENCIES ! ------------ use soil, only: layer_PEM, mlayer_PEM, mthermdiff_PEM, thermdiff_PEM, coefq_PEM, coefd_PEM, mu_PEM, alph_PEM, beta_PEM, fluxgeo use comsoil_h, only: volcapa ! DECLARATION ! ----------- implicit none #include "dimensions.h" ! ARGUMENTS ! --------- integer, intent(in) :: ngrid ! number of (horizontal) grid-points integer, intent(in) :: nsoil ! number of soil layers logical, intent(in) :: firstcall ! identifier for initialization call real, dimension(ngrid,nsoil), intent(in) :: therm_i ! thermal inertia [SI] real, intent(in) :: timestep ! time step [s] real, dimension(ngrid), intent(in) :: tsurf ! surface temperature [K] real, dimension(ngrid,nsoil), intent(inout) :: tsoil ! soil (mid-layer) temperature [K] ! LOCAL VARIABLES ! --------------- integer :: ig, ik ! CODE ! ---- ! 0. Initialisations and preprocessing step if (firstcall) then ! 0.1 Build mthermdiff_PEM(:), the mid-layer thermal diffusivities do ig = 1,ngrid do ik = 0,nsoil - 1 mthermdiff_PEM(ig,ik) = therm_i(ig,ik + 1)*therm_i(ig,ik + 1)/volcapa enddo enddo ! 0.2 Build thermdiff(:), the "interlayer" thermal diffusivities do ig = 1,ngrid do ik = 1,nsoil - 1 thermdiff_PEM(ig,ik) = ((layer_PEM(ik) - mlayer_PEM(ik - 1))*mthermdiff_PEM(ig,ik) & + (mlayer_PEM(ik) - layer_PEM(ik))*mthermdiff_PEM(ig,ik - 1))/(mlayer_PEM(ik) - mlayer_PEM(ik - 1)) enddo enddo ! 0.3 Build coefficients mu_PEM, q_{k+1/2}, d_k, alph_PEM ! mu_PEM mu_PEM = mlayer_PEM(0)/(mlayer_PEM(1) - mlayer_PEM(0)) ! q_{1/2} coefq_PEM(0) = volcapa*layer_PEM(1)/timestep ! q_{k+1/2} do ik = 1,nsoil - 1 coefq_PEM(ik) = volcapa*(layer_PEM(ik + 1) - layer_PEM(ik))/timestep enddo do ig = 1,ngrid ! d_k do ik = 1,nsoil - 1 coefd_PEM(ig,ik) = thermdiff_PEM(ig,ik)/(mlayer_PEM(ik)-mlayer_PEM(ik - 1)) enddo ! alph_PEM_{N-1} alph_PEM(ig,nsoil - 1) = coefd_PEM(ig,nsoil-1)/(coefq_PEM(nsoil - 1) + coefd_PEM(ig,nsoil - 1)) ! alph_PEM_k do ik = nsoil - 2,1,-1 alph_PEM(ig,ik) = coefd_PEM(ig,ik)/(coefq_PEM(ik) + coefd_PEM(ig,ik + 1)*(1. - alph_PEM(ig,ik + 1)) + coefd_PEM(ig,ik)) enddo enddo ! of do ig=1,ngrid endif ! of if (firstcall) if (.not. firstcall) THEN ! 2. Compute soil temperatures ! First layer: do ig = 1,ngrid tsoil(ig,1) = (tsurf(ig) + mu_PEM*beta_PEM(ig,1)*thermdiff_PEM(ig,1)/mthermdiff_PEM(ig,0))/ & (1. + mu_PEM*(1. - alph_PEM(ig,1))*thermdiff_PEM(ig,1)/mthermdiff_PEM(ig,0)) ! Other layers: do ik = 1,nsoil - 1 tsoil(ig,ik + 1) = alph_PEM(ig,ik)*tsoil(ig,ik) + beta_PEM(ig,ik) enddo enddo endif ! 2. Compute beta_PEM coefficients (preprocessing for next time step) ! Bottom layer, beta_PEM_{N-1} do ig = 1,ngrid beta_PEM(ig,nsoil - 1) = coefq_PEM(nsoil - 1)*tsoil(ig,nsoil)/(coefq_PEM(nsoil - 1) + coefd_PEM(ig,nsoil - 1)) & + fluxgeo/(coefq_PEM(nsoil - 1) + coefd_PEM(ig,nsoil - 1)) enddo ! Other layers do ik = nsoil-2,1,-1 do ig = 1,ngrid beta_PEM(ig,ik) = (coefq_PEM(ik)*tsoil(ig,ik + 1) + coefd_PEM(ig,ik + 1)*beta_PEM(ig,ik + 1))/ & (coefq_PEM(ik) + coefd_PEM(ig,ik + 1)*(1. - alph_PEM(ig,ik + 1)) + coefd_PEM(ig,ik)) enddo enddo END SUBROUTINE compute_tsoil !======================================================================= !======================================================================= SUBROUTINE ini_tsoil_pem(ngrid,nsoil,therm_i,tsurf,tsoil) !----------------------------------------------------------------------- ! NAME ! ini_tsoil_pem ! ! DESCRIPTION ! Initialize soil with the solution of the stationary heat conduction problem. ! Boundary conditions: Tsurf averaged from PCM; geothermal flux at bottom. ! ! AUTHORS & DATE ! L. Lang, 2023 ! JB Clement, 2023-2025 ! ! NOTES ! !----------------------------------------------------------------------- ! DEPENDENCIES ! ------------ use soil, only: layer_PEM, mlayer_PEM, mthermdiff_PEM, thermdiff_PEM, coefq_PEM, coefd_PEM, mu_PEM, alph_PEM, beta_PEM, fluxgeo use comsoil_h, only: volcapa ! DECLARATION ! ----------- implicit none #include "dimensions.h" ! ARGUMENTS ! --------- integer, intent(in) :: ngrid ! number of (horizontal) grid-points integer, intent(in) :: nsoil ! number of soil layers real, dimension(ngrid,nsoil), intent(in) :: therm_i ! thermal inertia [SI] real, dimension(ngrid), intent(in) :: tsurf ! surface temperature [K] real, dimension(ngrid,nsoil), intent(inout) :: tsoil ! soil (mid-layer) temperature [K] ! LOCAL VARIABLES ! --------------- integer :: ig, ik, iloop ! CODE ! ---- ! 0. Initialisations and preprocessing step ! 0.1 Build mthermdiff_PEM(:), the mid-layer thermal diffusivities do ig = 1,ngrid do ik = 0,nsoil - 1 mthermdiff_PEM(ig,ik) = therm_i(ig,ik + 1)*therm_i(ig,ik + 1)/volcapa enddo enddo ! 0.2 Build thermdiff(:), the "interlayer" thermal diffusivities do ig = 1,ngrid do ik = 1,nsoil - 1 thermdiff_PEM(ig,ik) = ((layer_PEM(ik) - mlayer_PEM(ik - 1))*mthermdiff_PEM(ig,ik) & + (mlayer_PEM(ik) - layer_PEM(ik))*mthermdiff_PEM(ig,ik - 1))/(mlayer_PEM(ik) - mlayer_PEM(ik - 1)) enddo enddo ! 0.3 Build coefficients mu_PEM, q_{k+1/2}, d_k, alph_PEM ! mu_PEM mu_PEM = mlayer_PEM(0)/(mlayer_PEM(1) - mlayer_PEM(0)) ! q_{1/2} coefq_PEM(:) = 0. ! q_{k+1/2} do ig = 1,ngrid ! d_k do ik = 1,nsoil - 1 coefd_PEM(ig,ik) = thermdiff_PEM(ig,ik)/(mlayer_PEM(ik) - mlayer_PEM(ik - 1)) enddo ! alph_PEM_{N-1} alph_PEM(ig,nsoil - 1) = coefd_PEM(ig,nsoil - 1)/(coefq_PEM(nsoil - 1) + coefd_PEM(ig,nsoil - 1)) ! alph_PEM_k do ik = nsoil - 2,1,-1 alph_PEM(ig,ik) = coefd_PEM(ig,ik)/(coefq_PEM(ik) + coefd_PEM(ig,ik + 1)*(1. - alph_PEM(ig,ik + 1)) + coefd_PEM(ig,ik)) enddo enddo ! of do ig=1,ngrid ! 1. Compute beta_PEM coefficients ! Bottom layer, beta_PEM_{N-1} do ig = 1,ngrid beta_PEM(ig,nsoil - 1) = coefq_PEM(nsoil - 1)*tsoil(ig,nsoil)/(coefq_PEM(nsoil - 1) + coefd_PEM(ig,nsoil - 1)) & + fluxgeo/(coefq_PEM(nsoil - 1) + coefd_PEM(ig,nsoil - 1)) enddo ! Other layers do ik = nsoil - 2,1,-1 do ig = 1,ngrid beta_PEM(ig,ik) = (coefq_PEM(ik)*tsoil(ig,ik + 1) + coefd_PEM(ig,ik+1)*beta_PEM(ig,ik + 1))/ & (coefq_PEM(ik) + coefd_PEM(ig,ik + 1)*(1. - alph_PEM(ig,ik + 1)) + coefd_PEM(ig,ik)) enddo enddo ! 2. Compute soil temperatures do iloop = 1,10 !just convergence do ig = 1,ngrid ! First layer: tsoil(ig,1) = (tsurf(ig) + mu_PEM*beta_PEM(ig,1)*thermdiff_PEM(ig,1)/mthermdiff_PEM(ig,0))/ & (1. + mu_PEM*(1. - alph_PEM(ig,1))*thermdiff_PEM(ig,1)/mthermdiff_PEM(ig,0)) ! Other layers: do ik = 1,nsoil - 1 tsoil(ig,ik + 1) = alph_PEM(ig,ik)*tsoil(ig,ik) + beta_PEM(ig,ik) enddo enddo enddo ! iloop END SUBROUTINE ini_tsoil_pem !======================================================================= !======================================================================= SUBROUTINE shift_tsoil2surf(ngrid,nsoil,nslope,zshift_surf,zlag,tsurf,tsoil) !----------------------------------------------------------------------- ! NAME ! shift_tsoil2surf ! ! DESCRIPTION ! Shift soil temperature profile to follow surface evolution due to ice ! condensation/sublimation. ! ! AUTHORS & DATE ! JB Clement, 2025 ! ! NOTES ! !----------------------------------------------------------------------- ! DEPENDENCIES ! ------------ use soil, only: layer_PEM, mlayer_PEM, fluxgeo, thermdiff_PEM, mthermdiff_PEM use math_toolkit, only: solve_steady_heat ! DECLARATION ! ----------- implicit none ! ARGUMENTS ! --------- integer, intent(in) :: ngrid ! number of (horizontal) grid-points integer, intent(in) :: nsoil ! number of soil layers integer, intent(in) :: nslope ! number of sub-slopes real, dimension(ngrid,nslope), intent(in) :: zshift_surf ! elevation shift for the surface [m] real, dimension(ngrid,nslope), intent(in) :: zlag ! newly built lag thickness [m] real, dimension(ngrid,nslope), intent(in) :: tsurf ! surface temperature [K] real, dimension(ngrid,nsoil,nslope), intent(inout) :: tsoil ! soil (mid-layer) temperature [K] ! LOCAL VARIABLES ! --------------- integer :: ig, isoil, islope, ishift, ilag, iz real :: a, z, zshift_surfloc, tsoil_minus, mlayer_minus real, dimension(ngrid,nsoil,nslope) :: tsoil_old ! CODE ! ---- write(*,*) "> Shifting soil temperature profile to match surface evolution" tsoil_old = tsoil do ig = 1,ngrid do islope = 1,nslope zshift_surfloc = zshift_surf(ig,islope) if (zshift_surfloc >= 0.) then ! In case of the surface is higher than initially if (zshift_surfloc < mlayer_PEM(0)) then ! Surface change is too small to be taken into account ! Nothing to do; we keep the soil temperature profile else if (zshift_surfloc >= mlayer_PEM(nsoil - 1)) then ! Surface change is much larger than the discretization tsoil(ig,:,islope) = tsurf(ig,islope) else ishift = 0 do while (mlayer_PEM(ishift) <= zshift_surfloc) ishift = ishift + 1 ! mlayer indices begin at 0 so this the good index for tsoil! enddo ! The "new soil" temperature is set to tsurf tsoil(ig,:ishift,islope) = tsurf(ig,islope) do isoil = ishift + 1,nsoil ! Position in the old discretization of the depth z = mlayer_PEM(isoil - 1) - zshift_surfloc ! Interpolation of the temperature profile from the old discretization tsoil(ig,isoil,islope) = itp_tsoil(tsoil_old(ig,:,islope),tsurf(ig,islope),z) enddo endif else ! In case of the surface is lower than initially if (abs(zshift_surfloc) < mlayer_PEM(0)) then ! Surface change is too small to be taken into account ! Nothing to do; we keep the soil temperature profile else if (abs(zshift_surfloc) >= mlayer_PEM(nsoil - 1)) then ! Surface change is much larger than the discretization call solve_steady_heat(nsoil,mlayer_PEM,layer_PEM,mthermdiff_PEM(ig,:),thermdiff_PEM(ig,:),tsurf(ig,islope),fluxgeo,tsoil(ig,:,islope)) else if (zlag(ig,islope) < mlayer_PEM(0)) then ! The lag is too thin to be taken into account ilag = 0 else ilag = 0 do while (mlayer_PEM(ilag) <= zlag(ig,islope)) ilag = ilag + 1 ! mlayer indices begin at 0 so this the good index for tsoil! enddo ! Position of the lag bottom in the old discretization of the depth z = zlag(ig,islope) - zshift_surfloc ! The "new lag" temperature is set to the ice temperature just below tsoil(ig,:ilag,islope) = itp_tsoil(tsoil_old(ig,:,islope),tsurf(ig,islope),z) endif ishift = nsoil - 1 z = mlayer_PEM(nsoil - 1) + zshift_surfloc do while (mlayer_PEM(ishift) >= z) ishift = ishift - 1 enddo ishift = ishift + 1 ! Adding 1 is needed to match the good index for tsoil! do isoil = ilag + 1,ishift ! Position in the old discretization of the depth z = mlayer_PEM(isoil - 1) - zshift_surfloc ! Interpolation of the temperature profile from the old discretization tsoil(ig,isoil,islope) = itp_tsoil(tsoil_old(ig,:,islope),tsurf(ig,islope),z) enddo ! The "new deepest layers" temperature is set by solving the steady heat equation call solve_steady_heat(nsoil - ishift + 1,mlayer_PEM(ishift - 1:),layer_PEM(ishift:),mthermdiff_PEM(ig,ishift - 1:),thermdiff_PEM(ig,ishift:),tsoil(ig,ishift,islope),fluxgeo,tsoil(ig,ishift:,islope)) endif endif enddo enddo END SUBROUTINE shift_tsoil2surf !======================================================================= !======================================================================= FUNCTION itp_tsoil(tsoil,tsurf,z) RESULT(tsoil_z) !----------------------------------------------------------------------- ! NAME ! itp_tsoil ! ! DESCRIPTION ! Interpolate soil temperature profile. ! ! AUTHORS & DATE ! JB Clement, 2025 ! ! NOTES ! !----------------------------------------------------------------------- ! DEPENDENCIES ! ------------ use soil, only: mlayer_PEM ! DECLARATION ! ----------- implicit none ! ARGUMENTS ! --------- real, dimension(:), intent(in) :: tsoil real, intent(in) :: z, tsurf ! LOCAL VARIABLES ! --------------- real :: tsoil_z, tsoil_minus, mlayer_minus, a integer :: iz ! CODE ! ---- ! Find the interval [mlayer_PEM(iz - 1),mlayer_PEM(iz)[ where the position z belongs iz = 0 do while (mlayer_PEM(iz) <= z) iz = iz + 1 enddo if (iz == 0) then tsoil_minus = tsurf mlayer_minus = 0. else tsoil_minus = tsoil(iz) mlayer_minus = mlayer_PEM(iz - 1) endif ! Interpolation of the temperature profile from the old discretization a = (tsoil(iz + 1) - tsoil_minus)/(mlayer_PEM(iz) - mlayer_minus) tsoil_z = a*(z - mlayer_minus) + tsoil_minus END FUNCTION itp_tsoil !======================================================================= END MODULE soil_temp