Changeset 3989 for trunk/LMDZ.COMMON/libf/evolution/soil.F90

Timestamp:

Dec 11, 2025, 12:56:05 PM (5 weeks ago)

Author:

jbclement

Message:

PEM:
Massive structural refactor of the PEM codebase for improved readability, consistency and maintainability. The goal is to modernize, standardize and consolidate the code while removing legacy complexity. In detail, this change:

• Performs large-scale cleanup removing unused code, obsolete routines, duplicated functionality and deprecated initialization logic;
• Removes "*_mod" wrappers;
• Replaces mixed naming conventions with clearer variable names, domain-based file/module names and purpose-based routine names;
• Adds native reading/writing capabilities to the PEM removing the cumbersome dependency on Mars PCM subroutines;
• Centralizes soil/ice/orbital/PEM configuration logic into dedicated modules;
• Simplifies routines and updates calls/interfaces to match the new structure.

This refactor significantly clarifies the codebase and provides a cleaner foundation for forthcoming developments.
JBC

File:

• trunk/LMDZ.COMMON/libf/evolution/soil.F90 (moved) (moved from trunk/LMDZ.COMMON/libf/evolution/comsoil_h_PEM.F90) (6 diffs)

trunk/LMDZ.COMMON/libf/evolution/soil.F90

@@ -1 @@
-MODULE comsoil_h_PEM
+MODULE soil
 
 implicit none
@@ -8 +8 @@
 real, allocatable, dimension(:,:,:) :: TI_PEM             ! soil thermal inertia [SI]
 real, allocatable, dimension(:,:)   :: inertiedat_PEM     ! soil thermal inertia saved as reference for current climate [SI]
-! Variables (FC: built in firstcall in soil.F)
 real, allocatable, dimension(:,:,:) :: tsoil_PEM          ! sub-surface temperatures [K]
-real, allocatable, dimension(:,:)   :: mthermdiff_PEM     ! (FC) mid-layer thermal diffusivity [SI]
-real, allocatable, dimension(:,:)   :: thermdiff_PEM      ! (FC) inter-layer thermal diffusivity [SI]
-real, allocatable, dimension(:)     :: coefq_PEM          ! (FC) q_{k+1/2} coefficients [SI]
-real, allocatable, dimension(:,:)   :: coefd_PEM          ! (FC) d_k coefficients [SI]
-real, allocatable, dimension(:,:)   :: alph_PEM           ! (FC) alpha_k coefficients [SI]
+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]
@@ -22 +21 @@
 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                             :: soil_pem           ! True by default, to run with the subsurface physic. Read in run_PEM.def
-real                                :: ini_huge_h2oice    ! Initial value for huge reservoirs of H2O ice [kg/m^2]
+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
 
@@ -30 +29 @@
 !=======================================================================
 
-SUBROUTINE ini_comsoil_h_PEM(ngrid,nslope)
+SUBROUTINE ini_soil(ngrid,nslope)
 
 implicit none
@@ -49 +48 @@
 allocate(inertiedat_PEM(ngrid,nsoilmx_PEM))
 
-END SUBROUTINE ini_comsoil_h_PEM
+END SUBROUTINE ini_soil
 
 !=======================================================================
 
-SUBROUTINE end_comsoil_h_PEM
+SUBROUTINE end_soil
 
 implicit none
@@ -69 +68 @@
 if (allocated(inertiedat_PEM)) deallocate(inertiedat_PEM)
 
-END SUBROUTINE end_comsoil_h_PEM
+END SUBROUTINE end_soil
 
-END MODULE comsoil_h_PEM
+!=======================================================================
+SUBROUTINE set_soil(ngrid,nslope,nsoil_PEM,nsoil_PCM,TI_PCM,TI_PEM)
+
+use iostart, only: inquire_field_ndims, get_var, get_field, inquire_field, inquire_dimension_length
+
+implicit none
+
+!=======================================================================
+!  Author: LL, based on work by Ehouarn Millour (07/2006)
+!
+!  Purpose: Read and/or initialise soil depths and properties
+!
+! 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
+!=======================================================================
+
+!=======================================================================
+!  arguments
+!  ---------
+!  inputs:
+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
+!=======================================================================
+! 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