MODULE compute_soiltemp_mod implicit none !----------------------------------------------------------------------- ! Author: LL ! Purpose: This module gathers the different routines used in the PEM to compute the soil temperature evolution and initialisation. ! ! Note: depths of layers and mid-layers, soil thermal inertia and ! heat capacity are commons in comsoil_PEM.h !----------------------------------------------------------------------- contains !======================================================================= SUBROUTINE compute_tsoil_pem(ngrid,nsoil,firstcall,therm_i,timestep,tsurf,tsoil) use comsoil_h_PEM, 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 implicit none !----------------------------------------------------------------------- ! Author: LL ! Purpose: Compute soil temperature using an implict 1st order scheme ! ! Note: depths of layers and mid-layers, soil thermal inertia and ! heat capacity are commons in comsoil_PEM.h !----------------------------------------------------------------------- #include "dimensions.h" ! Inputs: ! ------- 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] ! Outputs: !--------- real, dimension(ngrid,nsoil), intent(inout) :: tsoil ! soil (mid-layer) temperature [K] ! Local: !------- integer :: ig, ik ! 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_pem !======================================================================= SUBROUTINE ini_tsoil_pem(ngrid,nsoil,therm_i,tsurf,tsoil) use comsoil_h_PEM, 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 implicit none !----------------------------------------------------------------------- ! Author: LL ! Purpose: Initialize the soil with the solution of the stationnary problem of Heat Conduction. Boundarry conditions: Tsurf averaged from the PCM; Geothermal flux at the bottom layer ! ! Note: depths of layers and mid-layers, soil thermal inertia and ! heat capacity are commons in comsoil_PEM.h !----------------------------------------------------------------------- #include "dimensions.h" ! Inputs: !-------- 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] ! Outputs: !--------- real, dimension(ngrid,nsoil), intent(inout) :: tsoil ! soil (mid-layer) temperature [K] ! Local: !------- integer :: ig, ik, iloop ! 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) use comsoil_h_PEM, only: layer_PEM, mlayer_PEM, fluxgeo, thermdiff_PEM, mthermdiff_PEM use math_mod, only: solve_steady_heat implicit none !----------------------------------------------------------------------- ! Author: JBC ! Purpose: Shifting the soil temperature profile to follow the surface evolution due to ice condensation/sublimation !----------------------------------------------------------------------- ! Inputs: ! ------- 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] ! Outputs: ! -------- real, dimension(ngrid,nsoil,nslope), intent(inout) :: tsoil ! soil (mid-layer) temperature [K] ! Local: ! ------ integer :: ig, isoil, islope, ishift, ilag, iz real :: a, z, zshift_surfloc real, dimension(ngrid,nsoil,nslope) :: tsoil_old write(*,*)"Shifting soil temperature to surface" 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 ! 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 ! Interpolation of the temperature profile from the old discretization a = (tsoil_old(ig,iz + 1,islope) - tsoil_old(ig,iz,islope))/(mlayer_PEM(iz) - mlayer_PEM(iz - 1)) tsoil(ig,isoil,islope) = a*(z - mlayer_PEM(iz - 1)) + tsoil_old(ig,iz,islope) 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 ! 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 ! The "new lag" temperature is set to the ice temperature just below a = (tsoil_old(ig,iz + 1,islope) - tsoil_old(ig,iz,islope))/(mlayer_PEM(iz) - mlayer_PEM(iz - 1)) tsoil(ig,:ilag,islope) = a*(z - mlayer_PEM(iz - 1)) + tsoil_old(ig,iz,islope) 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 ! 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 ! Interpolation of the temperature profile from the old discretization a = (tsoil_old(ig,iz + 1,islope) - tsoil_old(ig,iz,islope))/(mlayer_PEM(iz) - mlayer_PEM(iz - 1)) tsoil(ig,isoil,islope) = a*(z - mlayer_PEM(iz - 1)) + tsoil_old(ig,iz,islope) 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 END MODULE compute_soiltemp_mod