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

Last change on this file since 3995 was 3991, checked in by jbclement, 4 weeks ago

PEM:
Apply documentation template everywhere: standardized headers format with short description, separators between functions/subroutines, normalized code sections, aligned dependencies/arguments/variables declaration.
JBC
File size: 8.5 KB
Line 
1MODULE soil
2!-----------------------------------------------------------------------
3! NAME
4!     soil
5!
6! DESCRIPTION
7!     Soil state and properties for PEM.
8!
9! AUTHORS & DATE
10!     L. Lange, 04/2023
11!     JB Clement, 2023-2025
12!
13! NOTES
14!
15!-----------------------------------------------------------------------
16
17! DECLARATION
18! -----------
19implicit none
20
21! MODULE PARAMETERS
22! -----------------
23integer, parameter :: nsoilmx_PEM = 68 ! number of layers in the PEM
24
25! MODULE VARIABLES
26! -----------------
27real, allocatable, dimension(:)     :: layer_PEM          ! soil layer depths [m]
28real, allocatable, dimension(:)     :: mlayer_PEM         ! soil mid-layer depths [m]
29real, allocatable, dimension(:,:,:) :: TI_PEM             ! soil thermal inertia [SI]
30real, allocatable, dimension(:,:)   :: inertiedat_PEM     ! soil thermal inertia saved as reference for current climate [SI]
31real, allocatable, dimension(:,:,:) :: tsoil_PEM          ! sub-surface temperatures [K]
32real, allocatable, dimension(:,:)   :: mthermdiff_PEM     ! mid-layer thermal diffusivity [SI]
33real, allocatable, dimension(:,:)   :: thermdiff_PEM      ! inter-layer thermal diffusivity [SI]
34real, allocatable, dimension(:)     :: coefq_PEM          ! q_{k+1/2} coefficients [SI]
35real, allocatable, dimension(:,:)   :: coefd_PEM          ! d_k coefficients [SI]
36real, allocatable, dimension(:,:)   :: alph_PEM           ! alpha_k coefficients [SI]
37real, allocatable, dimension(:,:)   :: beta_PEM           ! beta_k coefficients [SI]
38real                                :: mu_PEM             ! mu coefficient [SI]
39real                                :: fluxgeo            ! Geothermal flux [W/m^2]
40real                                :: depth_breccia      ! Depth at which we have breccia [m]
41real                                :: depth_bedrock      ! Depth at which we have bedrock [m]
42integer                             :: index_breccia      ! last index of the depth grid before having breccia
43integer                             :: index_bedrock      ! last index of the depth grid before having bedrock
44logical                             :: do_soil            ! True by default, to run with the subsurface physic. Read in run_PEM.def
45real                                :: h2oice_huge_ini    ! Initial value for huge reservoirs of H2O ice [kg/m^2]
46logical                             :: reg_thprop_dependp ! thermal properites of the regolith vary with the surface pressure
47
48contains
49!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
50
51!=======================================================================
52SUBROUTINE ini_soil(ngrid,nslope)
53!-----------------------------------------------------------------------
54! NAME
55!     ini_soil
56!
57! DESCRIPTION
58!     Allocate soil arrays based on grid dimensions.
59!
60! AUTHORS & DATE
61!     L. Lange, 2023
62!     JB Clement, 2023-2025
63!
64! NOTES
65!
66!-----------------------------------------------------------------------
67
68! DECLARATION
69! -----------
70implicit none
71
72! ARGUMENTS
73! ---------
74integer, intent(in) :: ngrid  ! number of atmospheric columns
75integer, intent(in) :: nslope ! number of slope within a mesh
76
77! CODE
78! ----
79allocate(layer_PEM(nsoilmx_PEM))
80allocate(mlayer_PEM(0:nsoilmx_PEM - 1))
81allocate(TI_PEM(ngrid,nsoilmx_PEM,nslope))
82allocate(tsoil_PEM(ngrid,nsoilmx_PEM,nslope))
83allocate(mthermdiff_PEM(ngrid,0:nsoilmx_PEM - 1))
84allocate(thermdiff_PEM(ngrid,nsoilmx_PEM - 1))
85allocate(coefq_PEM(0:nsoilmx_PEM - 1))
86allocate(coefd_PEM(ngrid,nsoilmx_PEM - 1))
87allocate(alph_PEM(ngrid,nsoilmx_PEM - 1))
88allocate(beta_PEM(ngrid,nsoilmx_PEM - 1))
89allocate(inertiedat_PEM(ngrid,nsoilmx_PEM))
90
91END SUBROUTINE ini_soil
92!=======================================================================
93
94!=======================================================================
95SUBROUTINE end_soil
96!-----------------------------------------------------------------------
97! NAME
98!     end_soil
99!
100! DESCRIPTION
101!     Deallocate all soil arrays.
102!
103! AUTHORS & DATE
104!     L. Lange, 2023
105!     JB Clement, 2023-2025
106!
107! NOTES
108!
109!-----------------------------------------------------------------------
110
111! DECLARATION
112! -----------
113implicit none
114
115! CODE
116! ----
117if (allocated(layer_PEM)) deallocate(layer_PEM)
118if (allocated(mlayer_PEM)) deallocate(mlayer_PEM)
119if (allocated(TI_PEM)) deallocate(TI_PEM)
120if (allocated(tsoil_PEM)) deallocate(tsoil_PEM)
121if (allocated(mthermdiff_PEM)) deallocate(mthermdiff_PEM)
122if (allocated(thermdiff_PEM)) deallocate(thermdiff_PEM)
123if (allocated(coefq_PEM)) deallocate(coefq_PEM)
124if (allocated(coefd_PEM)) deallocate(coefd_PEM)
125if (allocated(alph_PEM)) deallocate(alph_PEM)
126if (allocated(beta_PEM)) deallocate(beta_PEM)
127if (allocated(inertiedat_PEM)) deallocate(inertiedat_PEM)
128
129END SUBROUTINE end_soil
130!=======================================================================
131
132!=======================================================================
133SUBROUTINE set_soil(ngrid,nslope,nsoil_PEM,nsoil_PCM,TI_PCM,TI_PEM)
134!-----------------------------------------------------------------------
135! NAME
136!     set_soil
137!
138! DESCRIPTION
139!     Initialize soil depths and properties. Builds layer depths using
140!     exponential then power-law distribution; interpolates thermal inertia
141!     from PCM to PEM grid.
142!
143! AUTHORS & DATE
144!     E. Millour, 07/2006
145!     L. Lange, 2023
146!     JB Clement, 2023-2025
147!
148! NOTES
149!     Modifications:
150!          Aug.2010 EM: use NetCDF90 to load variables (enables using
151!                       r4 or r8 restarts independently of having compiled
152!                       the GCM in r4 or r8)
153!          June 2013 TN: Possibility to read files with a time axis
154!     The various actions and variable read/initialized are:
155!     1. Read/build layer (and midlayer) depth
156!     2. Interpolate thermal inertia and temperature on the new grid, if necessary
157!-----------------------------------------------------------------------
158
159! DEPENDENCIES
160! ------------
161use iostart, only: inquire_field_ndims, get_var, get_field, inquire_field, inquire_dimension_length
162
163! DECLARATION
164! -----------
165implicit none
166
167! ARGUMENTS
168! ---------
169integer,                                 intent(in)    :: ngrid     ! # of horizontal grid points
170integer,                                 intent(in)    :: nslope    ! # of subslope wihtin the mesh
171integer,                                 intent(in)    :: nsoil_PEM ! # of soil layers in the PEM
172integer,                                 intent(in)    :: nsoil_PCM ! # of soil layers in the PCM
173real, dimension(ngrid,nsoil_PCM,nslope), intent(in)    :: TI_PCM    ! Thermal inertia in the PCM [SI]
174real, dimension(ngrid,nsoil_PEM,nslope), intent(inout) :: TI_PEM    ! Thermal inertia in the PEM [SI]
175
176! LOCAL VARIABLES
177! ---------------
178integer :: ig, iloop, islope ! loop counters
179real    :: alpha, lay1       ! coefficients for building layers
180real    :: index_powerlaw    ! coefficient to match the power low grid with the exponential one
181
182! CODE
183! ----
184! 1. Depth coordinate
185! -------------------
186! mlayer_PEM distribution: For the grid in common with the PCM: lay1*(1+k^(2.9)*(1-exp(-k/20))
187! Then we use a power low : lay1*alpha**(k-1/2)
188lay1 = 2.e-4
189alpha = 2
190do iloop = 0,nsoil_PCM - 1
191    mlayer_PEM(iloop) = lay1*(1 + iloop**2.9*(1 - exp(-real(iloop)/20.)))
192enddo
193
194do iloop = 0,nsoil_PEM - 1
195    if (lay1*(alpha**(iloop-0.5)) > mlayer_PEM(nsoil_PCM - 1)) then
196        index_powerlaw = iloop
197        exit
198    endif
199enddo
200
201do iloop = nsoil_PCM,nsoil_PEM - 1
202    mlayer_PEM(iloop) = lay1*(alpha**(index_powerlaw + (iloop - nsoil_PCM) - 0.5))
203enddo
204
205! Build layer_PEM()
206do iloop = 1,nsoil_PEM - 1
207    layer_PEM(iloop) = (mlayer_PEM(iloop) + mlayer_PEM(iloop - 1))/2.
208enddo
209layer_PEM(nsoil_PEM) = 2*mlayer_PEM(nsoil_PEM - 1) - mlayer_PEM(nsoil_PEM - 2)
210
211
212! 2. Thermal inertia (note: it is declared in comsoil_h)
213! ------------------
214do ig = 1,ngrid
215    do islope = 1,nslope
216        do iloop = 1,nsoil_PCM
217            TI_PEM(ig,iloop,islope) = TI_PCM(ig,iloop,islope)
218        enddo
219        if (nsoil_PEM > nsoil_PCM) then
220            do iloop = nsoil_PCM + 1,nsoil_PEM
221                TI_PEM(ig,iloop,islope) = TI_PCM(ig,nsoil_PCM,islope)
222            enddo
223        endif
224    enddo
225enddo
226
227! 3. Index for subsurface layering
228! ------------------
229index_breccia = 1
230do iloop = 1,nsoil_PEM - 1
231    if (depth_breccia >= layer_PEM(iloop)) then
232        index_breccia = iloop
233    else
234        exit
235    endif
236enddo
237
238index_bedrock = 1
239do iloop = 1,nsoil_PEM - 1
240    if (depth_bedrock >= layer_PEM(iloop)) then
241        index_bedrock = iloop
242    else
243        exit
244    endif
245enddo
246
247END SUBROUTINE set_soil
248!=======================================================================
249
250END MODULE soil
Note: See TracBrowser for help on using the repository browser.