source: trunk/LMDZ.COMMON/libf/evolution/soil.F90 @ 3989

Last change on this file since 3989 was 3989, checked in by jbclement, 5 weeks ago

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 size: 7.1 KB
Line 
1MODULE soil
2
3implicit none
4
5integer, parameter                  :: nsoilmx_PEM = 68   ! number of layers in the PEM
6real, allocatable, dimension(:)     :: layer_PEM          ! soil layer depths [m]
7real, allocatable, dimension(:)     :: mlayer_PEM         ! soil mid-layer depths [m]
8real, allocatable, dimension(:,:,:) :: TI_PEM             ! soil thermal inertia [SI]
9real, allocatable, dimension(:,:)   :: inertiedat_PEM     ! soil thermal inertia saved as reference for current climate [SI]
10real, allocatable, dimension(:,:,:) :: tsoil_PEM          ! sub-surface temperatures [K]
11real, allocatable, dimension(:,:)   :: mthermdiff_PEM     ! mid-layer thermal diffusivity [SI]
12real, allocatable, dimension(:,:)   :: thermdiff_PEM      ! inter-layer thermal diffusivity [SI]
13real, allocatable, dimension(:)     :: coefq_PEM          ! q_{k+1/2} coefficients [SI]
14real, allocatable, dimension(:,:)   :: coefd_PEM          ! d_k coefficients [SI]
15real, allocatable, dimension(:,:)   :: alph_PEM           ! alpha_k coefficients [SI]
16real, allocatable, dimension(:,:)   :: beta_PEM           ! beta_k coefficients [SI]
17real                                :: mu_PEM             ! mu coefficient [SI]
18real                                :: fluxgeo            ! Geothermal flux [W/m^2]
19real                                :: depth_breccia      ! Depth at which we have breccia [m]
20real                                :: depth_bedrock      ! Depth at which we have bedrock [m]
21integer                             :: index_breccia      ! last index of the depth grid before having breccia
22integer                             :: index_bedrock      ! last index of the depth grid before having bedrock
23logical                             :: do_soil            ! True by default, to run with the subsurface physic. Read in run_PEM.def
24real                                :: h2oice_huge_ini    ! Initial value for huge reservoirs of H2O ice [kg/m^2]
25logical                             :: reg_thprop_dependp ! thermal properites of the regolith vary with the surface pressure
26
27!=======================================================================
28contains
29!=======================================================================
30
31SUBROUTINE ini_soil(ngrid,nslope)
32
33implicit none
34
35integer, intent(in) :: ngrid  ! number of atmospheric columns
36integer, intent(in) :: nslope ! number of slope within a mesh
37
38allocate(layer_PEM(nsoilmx_PEM))
39allocate(mlayer_PEM(0:nsoilmx_PEM - 1))
40allocate(TI_PEM(ngrid,nsoilmx_PEM,nslope))
41allocate(tsoil_PEM(ngrid,nsoilmx_PEM,nslope))
42allocate(mthermdiff_PEM(ngrid,0:nsoilmx_PEM - 1))
43allocate(thermdiff_PEM(ngrid,nsoilmx_PEM - 1))
44allocate(coefq_PEM(0:nsoilmx_PEM - 1))
45allocate(coefd_PEM(ngrid,nsoilmx_PEM - 1))
46allocate(alph_PEM(ngrid,nsoilmx_PEM - 1))
47allocate(beta_PEM(ngrid,nsoilmx_PEM - 1))
48allocate(inertiedat_PEM(ngrid,nsoilmx_PEM))
49
50END SUBROUTINE ini_soil
51
52!=======================================================================
53
54SUBROUTINE end_soil
55
56implicit none
57
58if (allocated(layer_PEM)) deallocate(layer_PEM)
59if (allocated(mlayer_PEM)) deallocate(mlayer_PEM)
60if (allocated(TI_PEM)) deallocate(TI_PEM)
61if (allocated(tsoil_PEM)) deallocate(tsoil_PEM)
62if (allocated(mthermdiff_PEM)) deallocate(mthermdiff_PEM)
63if (allocated(thermdiff_PEM)) deallocate(thermdiff_PEM)
64if (allocated(coefq_PEM)) deallocate(coefq_PEM)
65if (allocated(coefd_PEM)) deallocate(coefd_PEM)
66if (allocated(alph_PEM)) deallocate(alph_PEM)
67if (allocated(beta_PEM)) deallocate(beta_PEM)
68if (allocated(inertiedat_PEM)) deallocate(inertiedat_PEM)
69
70END SUBROUTINE end_soil
71
72!=======================================================================
73SUBROUTINE set_soil(ngrid,nslope,nsoil_PEM,nsoil_PCM,TI_PCM,TI_PEM)
74
75use iostart, only: inquire_field_ndims, get_var, get_field, inquire_field, inquire_dimension_length
76
77implicit none
78
79!=======================================================================
80!  Author: LL, based on work by Ehouarn Millour (07/2006)
81!
82!  Purpose: Read and/or initialise soil depths and properties
83!
84! Modifications: Aug.2010 EM: use NetCDF90 to load variables (enables using
85!                             r4 or r8 restarts independently of having compiled
86!                             the GCM in r4 or r8)
87!                June 2013 TN: Possibility to read files with a time axis
88!
89!  The various actions and variable read/initialized are:
90!  1. Read/build layer (and midlayer) depth
91!  2. Interpolate thermal inertia and temperature on the new grid, if necessary
92!=======================================================================
93
94!=======================================================================
95!  arguments
96!  ---------
97!  inputs:
98integer,                                 intent(in) :: ngrid     ! # of horizontal grid points
99integer,                                 intent(in) :: nslope    ! # of subslope wihtin the mesh
100integer,                                 intent(in) :: nsoil_PEM ! # of soil layers in the PEM
101integer,                                 intent(in) :: nsoil_PCM ! # of soil layers in the PCM
102real, dimension(ngrid,nsoil_PCM,nslope), intent(in) :: TI_PCM    ! Thermal inertia  in the PCM [SI]
103
104real, dimension(ngrid,nsoil_PEM,nslope), intent(inout) :: TI_PEM ! Thermal inertia   in the PEM [SI]
105!=======================================================================
106! local variables:
107integer :: ig, iloop, islope ! loop counters
108real    :: alpha, lay1       ! coefficients for building layers
109real    :: index_powerlaw    ! coefficient to match the power low grid with the exponential one
110!=======================================================================
111! 1. Depth coordinate
112! -------------------
113! mlayer_PEM distribution: For the grid in common with the PCM: lay1*(1+k^(2.9)*(1-exp(-k/20))
114! Then we use a power low : lay1*alpha**(k-1/2)
115lay1 = 2.e-4
116alpha = 2
117do iloop = 0,nsoil_PCM - 1
118    mlayer_PEM(iloop) = lay1*(1 + iloop**2.9*(1 - exp(-real(iloop)/20.)))
119enddo
120
121do iloop = 0,nsoil_PEM - 1
122    if (lay1*(alpha**(iloop-0.5)) > mlayer_PEM(nsoil_PCM - 1)) then
123        index_powerlaw = iloop
124        exit
125    endif
126enddo
127
128do iloop = nsoil_PCM,nsoil_PEM - 1
129    mlayer_PEM(iloop) = lay1*(alpha**(index_powerlaw + (iloop - nsoil_PCM) - 0.5))
130enddo
131
132! Build layer_PEM()
133do iloop = 1,nsoil_PEM - 1
134    layer_PEM(iloop) = (mlayer_PEM(iloop) + mlayer_PEM(iloop - 1))/2.
135enddo
136layer_PEM(nsoil_PEM) = 2*mlayer_PEM(nsoil_PEM - 1) - mlayer_PEM(nsoil_PEM - 2)
137
138
139! 2. Thermal inertia (note: it is declared in comsoil_h)
140! ------------------
141do ig = 1,ngrid
142    do islope = 1,nslope
143        do iloop = 1,nsoil_PCM
144            TI_PEM(ig,iloop,islope) = TI_PCM(ig,iloop,islope)
145        enddo
146        if (nsoil_PEM > nsoil_PCM) then
147            do iloop = nsoil_PCM + 1,nsoil_PEM
148                TI_PEM(ig,iloop,islope) = TI_PCM(ig,nsoil_PCM,islope)
149            enddo
150        endif
151    enddo
152enddo
153
154! 3. Index for subsurface layering
155! ------------------
156index_breccia = 1
157do iloop = 1,nsoil_PEM - 1
158    if (depth_breccia >= layer_PEM(iloop)) then
159        index_breccia = iloop
160    else
161        exit
162    endif
163enddo
164
165index_bedrock = 1
166do iloop = 1,nsoil_PEM - 1
167    if (depth_bedrock >= layer_PEM(iloop)) then
168        index_bedrock = iloop
169    else
170        exit
171    endif
172enddo
173
174END SUBROUTINE set_soil
175
176END MODULE soil
Note: See TracBrowser for help on using the repository browser.