MODULE soil !----------------------------------------------------------------------- ! NAME ! soil ! ! DESCRIPTION ! Soil state and properties for PEM. ! ! AUTHORS & DATE ! L. Lange, 04/2023 ! JB Clement, 2023-2025 ! ! NOTES ! !----------------------------------------------------------------------- ! DECLARATION ! ----------- implicit none ! MODULE PARAMETERS ! ----------------- integer, parameter :: nsoilmx_PEM = 68 ! number of layers in the PEM ! MODULE VARIABLES ! ----------------- real, allocatable, dimension(:) :: layer_PEM ! soil layer depths [m] real, allocatable, dimension(:) :: mlayer_PEM ! soil mid-layer depths [m] real, allocatable, dimension(:,:,:) :: TI_PEM ! soil thermal inertia [SI] real, allocatable, dimension(:,:) :: inertiedat_PEM ! soil thermal inertia saved as reference for current climate [SI] real, allocatable, dimension(:,:,:) :: tsoil_PEM ! sub-surface temperatures [K] real, allocatable, dimension(:,:) :: mthermdiff_PEM ! mid-layer thermal diffusivity [SI] real, allocatable, dimension(:,:) :: thermdiff_PEM ! inter-layer thermal diffusivity [SI] real, allocatable, dimension(:) :: coefq_PEM ! q_{k+1/2} coefficients [SI] real, allocatable, dimension(:,:) :: coefd_PEM ! d_k coefficients [SI] real, allocatable, dimension(:,:) :: alph_PEM ! alpha_k coefficients [SI] real, allocatable, dimension(:,:) :: beta_PEM ! beta_k coefficients [SI] real :: mu_PEM ! mu coefficient [SI] real :: fluxgeo ! Geothermal flux [W/m^2] real :: depth_breccia ! Depth at which we have breccia [m] real :: depth_bedrock ! Depth at which we have bedrock [m] integer :: index_breccia ! last index of the depth grid before having breccia integer :: index_bedrock ! last index of the depth grid before having bedrock logical :: do_soil ! True by default, to run with the subsurface physic. Read in run_PEM.def real :: h2oice_huge_ini ! Initial value for huge reservoirs of H2O ice [kg/m^2] logical :: reg_thprop_dependp ! thermal properites of the regolith vary with the surface pressure contains !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ !======================================================================= SUBROUTINE ini_soil(ngrid,nslope) !----------------------------------------------------------------------- ! NAME ! ini_soil ! ! DESCRIPTION ! Allocate soil arrays based on grid dimensions. ! ! AUTHORS & DATE ! L. Lange, 2023 ! JB Clement, 2023-2025 ! ! NOTES ! !----------------------------------------------------------------------- ! DECLARATION ! ----------- implicit none ! ARGUMENTS ! --------- integer, intent(in) :: ngrid ! number of atmospheric columns integer, intent(in) :: nslope ! number of slope within a mesh ! CODE ! ---- allocate(layer_PEM(nsoilmx_PEM)) allocate(mlayer_PEM(0:nsoilmx_PEM - 1)) allocate(TI_PEM(ngrid,nsoilmx_PEM,nslope)) allocate(tsoil_PEM(ngrid,nsoilmx_PEM,nslope)) allocate(mthermdiff_PEM(ngrid,0:nsoilmx_PEM - 1)) allocate(thermdiff_PEM(ngrid,nsoilmx_PEM - 1)) allocate(coefq_PEM(0:nsoilmx_PEM - 1)) allocate(coefd_PEM(ngrid,nsoilmx_PEM - 1)) allocate(alph_PEM(ngrid,nsoilmx_PEM - 1)) allocate(beta_PEM(ngrid,nsoilmx_PEM - 1)) allocate(inertiedat_PEM(ngrid,nsoilmx_PEM)) 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_PEM)) deallocate(layer_PEM) if (allocated(mlayer_PEM)) deallocate(mlayer_PEM) if (allocated(TI_PEM)) deallocate(TI_PEM) if (allocated(tsoil_PEM)) deallocate(tsoil_PEM) if (allocated(mthermdiff_PEM)) deallocate(mthermdiff_PEM) if (allocated(thermdiff_PEM)) deallocate(thermdiff_PEM) if (allocated(coefq_PEM)) deallocate(coefq_PEM) if (allocated(coefd_PEM)) deallocate(coefd_PEM) if (allocated(alph_PEM)) deallocate(alph_PEM) if (allocated(beta_PEM)) deallocate(beta_PEM) if (allocated(inertiedat_PEM)) deallocate(inertiedat_PEM) END SUBROUTINE end_soil !======================================================================= !======================================================================= SUBROUTINE set_soil(ngrid,nslope,nsoil_PEM,nsoil_PCM,TI_PCM,TI_PEM) !----------------------------------------------------------------------- ! 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 iostart, only: inquire_field_ndims, get_var, get_field, inquire_field, inquire_dimension_length ! DECLARATION ! ----------- implicit none ! ARGUMENTS ! --------- integer, intent(in) :: ngrid ! # of horizontal grid points integer, intent(in) :: nslope ! # of subslope wihtin the mesh integer, intent(in) :: nsoil_PEM ! # of soil layers in the PEM integer, intent(in) :: nsoil_PCM ! # of soil layers in the PCM real, dimension(ngrid,nsoil_PCM,nslope), intent(in) :: TI_PCM ! Thermal inertia in the PCM [SI] real, dimension(ngrid,nsoil_PEM,nslope), intent(inout) :: TI_PEM ! Thermal inertia in the PEM [SI] ! LOCAL VARIABLES ! --------------- integer :: ig, iloop, islope ! loop counters real :: alpha, lay1 ! coefficients for building layers real :: index_powerlaw ! coefficient to match the power low grid with the exponential one ! CODE ! ---- ! 1. Depth coordinate ! ------------------- ! mlayer_PEM 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_PEM(iloop) = lay1*(1 + iloop**2.9*(1 - exp(-real(iloop)/20.))) enddo do iloop = 0,nsoil_PEM - 1 if (lay1*(alpha**(iloop-0.5)) > mlayer_PEM(nsoil_PCM - 1)) then index_powerlaw = iloop exit endif enddo do iloop = nsoil_PCM,nsoil_PEM - 1 mlayer_PEM(iloop) = lay1*(alpha**(index_powerlaw + (iloop - nsoil_PCM) - 0.5)) enddo ! Build layer_PEM() do iloop = 1,nsoil_PEM - 1 layer_PEM(iloop) = (mlayer_PEM(iloop) + mlayer_PEM(iloop - 1))/2. enddo layer_PEM(nsoil_PEM) = 2*mlayer_PEM(nsoil_PEM - 1) - mlayer_PEM(nsoil_PEM - 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_PEM(ig,iloop,islope) = TI_PCM(ig,iloop,islope) enddo if (nsoil_PEM > nsoil_PCM) then do iloop = nsoil_PCM + 1,nsoil_PEM TI_PEM(ig,iloop,islope) = TI_PCM(ig,nsoil_PCM,islope) enddo endif enddo enddo ! 3. Index for subsurface layering ! ------------------ index_breccia = 1 do iloop = 1,nsoil_PEM - 1 if (depth_breccia >= layer_PEM(iloop)) then index_breccia = iloop else exit endif enddo index_bedrock = 1 do iloop = 1,nsoil_PEM - 1 if (depth_bedrock >= layer_PEM(iloop)) then index_bedrock = iloop else exit endif enddo END SUBROUTINE set_soil !======================================================================= END MODULE soil