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

Last change on this file since 4147 was 4136, checked in by jbclement, 4 weeks ago

PEM:

  • Temporary memory load reduced in "xios_data".
  • Rename local variables in procedures to avoid masking variables in parent scope.

JBC
File size: 12.8 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! DEPENDENCIES
18! ------------
19use numerics, only: dp, di, k4
20
21! DECLARATION
22! -----------
23implicit none
24
25! PARAMETERS
26! ----------
27real(dp), parameter :: rho_regolith      = 2000._dp ! Density of the martian regolith [kg/m^3], Zent et al. 1995; Buhler and Piqueux 2021
28real(dp), parameter :: regolith_porosity = 0.4_dp   ! Porosity of the martian regolith for a random loose packing of monodisperse sphere, Scott 1960
29real(dp), parameter :: breccia_porosity  = 0.1_dp   ! Porosity of breccia
30
31logical(k4),                             protected          :: do_soil            ! To run with subsurface physics
32logical(k4),                             protected          :: reg_thprop_dependp ! Regolith thermal properties depend on surface pressure
33real(dp),                                protected          :: flux_geo           ! Geothermal flux [W/m^2]
34real(dp),                                protected          :: depth_breccia      ! Depth of breccia [m]
35real(dp),                                protected          :: depth_bedrock      ! Depth of bedrock [m]
36real(dp), allocatable, dimension(:,:,:), protected, private :: TI_PCM             ! Soil thermal inertia in the "startfi.nc" at the beginning [SI]
37
38! VARIABLES
39! ---------
40real(dp), allocatable, dimension(:)     :: layer              ! Soil layer depths [m]
41real(dp), allocatable, dimension(:)     :: mlayer             ! Soil mid-layer depths [m]
42real(dp), allocatable, dimension(:,:,:) :: TI                 ! Soil thermal inertia [SI]
43real(dp), allocatable, dimension(:,:)   :: inertiedat         ! Soil thermal inertia saved as reference for current climate [SI]
44real(dp), allocatable, dimension(:,:)   :: inertiedat_PCM     ! Soil thermal inertia saved as reference for current climate in the PCM [SI]
45real(dp), allocatable, dimension(:,:)   :: mthermdiff         ! Mid-layer thermal diffusivity [SI]
46real(dp), allocatable, dimension(:,:)   :: thermdiff          ! Inter-layer thermal diffusivity [SI]
47real(dp), allocatable, dimension(:)     :: coefq              ! q_{k+1/2} coefficients [SI]
48real(dp), allocatable, dimension(:,:)   :: coefd              ! d_k coefficients [SI]
49real(dp), allocatable, dimension(:,:)   :: alph               ! alpha_k coefficients [SI]
50real(dp), allocatable, dimension(:,:)   :: beta               ! beta_k coefficients [SI]
51real(dp)                                :: mu                 ! mu coefficient [SI]
52integer(di)                             :: index_breccia      ! Last index of the depth grid before having breccia
53integer(di)                             :: index_bedrock      ! Last index of the depth grid before having bedrock
54real(dp)                                :: volcapa            ! Soil volumetric heat capacity
55
56contains
57!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
58
59!=======================================================================
60SUBROUTINE ini_soil()
61!-----------------------------------------------------------------------
62! NAME
63!     ini_soil
64!
65! DESCRIPTION
66!     Allocate soil arrays based on grid dimensions.
67!
68! AUTHORS & DATE
69!     L. Lange, 2023
70!     JB Clement, 2023-2025
71!
72! NOTES
73!
74!-----------------------------------------------------------------------
75
76! DEPENDENCIES
77! ------------
78use geometry, only: ngrid, nslope, nsoil, nsoil_PCM
79
80! DECLARATION
81! -----------
82implicit none
83
84! CODE
85! ----
86allocate(layer(nsoil))
87allocate(mlayer(0:nsoil - 1))
88allocate(TI(ngrid,nsoil,nslope))
89allocate(mthermdiff(ngrid,0:nsoil - 1))
90allocate(thermdiff(ngrid,nsoil - 1))
91allocate(coefq(0:nsoil - 1))
92allocate(coefd(ngrid,nsoil - 1))
93allocate(alph(ngrid,nsoil - 1))
94allocate(beta(ngrid,nsoil - 1))
95allocate(inertiedat(ngrid,nsoil))
96allocate(inertiedat_PCM(ngrid,nsoil_PCM))
97allocate(TI_PCM(ngrid,nsoil_PCM,nslope))
98inertiedat_PCM(:,:) = 0._dp
99TI_PCM(:,:,:) = 0._dp
100
101END SUBROUTINE ini_soil
102!=======================================================================
103
104!=======================================================================
105SUBROUTINE end_soil
106!-----------------------------------------------------------------------
107! NAME
108!     end_soil
109!
110! DESCRIPTION
111!     Deallocate all soil arrays.
112!
113! AUTHORS & DATE
114!     L. Lange, 2023
115!     JB Clement, 2023-2025
116!
117! NOTES
118!
119!-----------------------------------------------------------------------
120
121! DECLARATION
122! -----------
123implicit none
124
125! CODE
126! ----
127if (allocated(layer)) deallocate(layer)
128if (allocated(mlayer)) deallocate(mlayer)
129if (allocated(TI)) deallocate(TI)
130if (allocated(mthermdiff)) deallocate(mthermdiff)
131if (allocated(thermdiff)) deallocate(thermdiff)
132if (allocated(coefq)) deallocate(coefq)
133if (allocated(coefd)) deallocate(coefd)
134if (allocated(alph)) deallocate(alph)
135if (allocated(beta)) deallocate(beta)
136if (allocated(inertiedat)) deallocate(inertiedat)
137if (allocated(inertiedat_PCM)) deallocate(inertiedat_PCM)
138if (allocated(TI_PCM)) deallocate(TI_PCM)
139
140END SUBROUTINE end_soil
141!=======================================================================
142
143!=======================================================================
144SUBROUTINE set_soil_config(do_soil_in,reg_thprop_dependp_in,flux_geo_in,depth_breccia_in,depth_bedrock_in)
145!-----------------------------------------------------------------------
146! NAME
147!     set_soil_config
148!
149! DESCRIPTION
150!     Setter for 'soil' configuration parameters.
151!
152! AUTHORS & DATE
153!     JB Clement, 02/2026
154!
155! NOTES
156!
157!-----------------------------------------------------------------------
158
159! DEPENDENCIES
160! ------------
161use utility,  only: real2str, bool2str
162use display,  only: print_msg, LVL_NFO
163use stoppage, only: stop_clean
164
165! DECLARATION
166! -----------
167implicit none
168
169! ARGUMENTS
170! ---------
171logical(k4), intent(in) :: do_soil_in, reg_thprop_dependp_in
172real(dp),    intent(in) :: flux_geo_in, depth_breccia_in, depth_bedrock_in
173
174! CODE
175! ----
176do_soil = do_soil_in
177reg_thprop_dependp = reg_thprop_dependp_in
178flux_geo = flux_geo_in
179depth_breccia = depth_breccia_in
180depth_bedrock = depth_bedrock_in
181call print_msg('do_soil            = '//bool2str(do_soil),LVL_NFO)
182call print_msg('reg_thprop_dependp = '//bool2str(reg_thprop_dependp),LVL_NFO)
183call print_msg('flux_geo           = '//real2str(flux_geo),LVL_NFO)
184call print_msg('depth_breccia      = '//real2str(depth_breccia),LVL_NFO)
185call print_msg('depth_bedrock      = '//real2str(depth_bedrock),LVL_NFO)
186if (depth_breccia < 0._dp) call stop_clean(__FILE__,__LINE__,'''depth_breccia'' must be positive!',1)
187if (depth_bedrock < 0._dp) call stop_clean(__FILE__,__LINE__,'''depth_bedrock'' must be positive!',1)
188if (depth_bedrock <= depth_breccia) call stop_clean(__FILE__,__LINE__,'''depth_bedrock'' must be greater than ''depth_breccia''!',1)
189
190END SUBROUTINE set_soil_config
191!=======================================================================
192
193!=======================================================================
194SUBROUTINE set_TI_PCM(TI_PCM_in)
195!-----------------------------------------------------------------------
196! NAME
197!     set_TI_PCM
198!
199! DESCRIPTION
200!     Setter for 'TI_PCM'.
201!
202! AUTHORS & DATE
203!     JB Clement, 01/2026
204!
205! NOTES
206!
207!-----------------------------------------------------------------------
208
209! DECLARATION
210! -----------
211implicit none
212
213! ARGUMENTS
214! ---------
215real(dp), dimension(:,:,:), intent(in) :: TI_PCM_in
216
217! CODE
218! ----
219TI_PCM(:,:,:) = TI_PCM_in(:,:,:)
220
221END SUBROUTINE set_TI_PCM
222!=======================================================================
223
224!=======================================================================
225SUBROUTINE set_inertiedat_PCM(inertiedat_PCM_in)
226!-----------------------------------------------------------------------
227! NAME
228!     set_inertiedat_PCM
229!
230! DESCRIPTION
231!     Setter for 'inertiedat_PCM'.
232!
233! AUTHORS & DATE
234!     JB Clement, 01/2026
235!
236! NOTES
237!
238!-----------------------------------------------------------------------
239
240! DECLARATION
241! -----------
242implicit none
243
244! ARGUMENTS
245! ---------
246real(dp), dimension(:,:), intent(in) :: inertiedat_PCM_in
247
248! CODE
249! ----
250inertiedat_PCM(:,:) = inertiedat_PCM_in(:,:)
251
252END SUBROUTINE set_inertiedat_PCM
253!=======================================================================
254
255!=======================================================================
256SUBROUTINE set_soil(therm_initia)
257!-----------------------------------------------------------------------
258! NAME
259!     set_soil
260!
261! DESCRIPTION
262!     Initialize soil depths and properties. Builds layer depths using
263!     exponential then power-law distribution; interpolates thermal inertia
264!     from PCM to PEM grid.
265!
266! AUTHORS & DATE
267!     E. Millour, 07/2006
268!     L. Lange, 2023
269!     JB Clement, 2023-2025
270!
271! NOTES
272!     Modifications:
273!          Aug.2010 EM: use NetCDF90 to load variables (enables using
274!                       r4 or r8 restarts independently of having compiled
275!                       the GCM in r4 or r8)
276!          June 2013 TN: Possibility to read files with a time axis
277!     The various actions and variable read/initialized are:
278!     1. Read/build layer (and midlayer) depth
279!     2. Interpolate thermal inertia and temperature on the new grid, if necessary
280!-----------------------------------------------------------------------
281
282! DEPENDENCIES
283! ------------
284use geometry, only: ngrid, nslope, nsoil, nsoil_PCM
285
286! DECLARATION
287! -----------
288implicit none
289
290! ARGUMENTS
291! ---------
292real(dp), dimension(:,:,:), intent(inout) :: therm_initia ! Thermal inertia in the PEM [SI]
293
294! LOCAL VARIABLES
295! ---------------
296integer(di) :: ig, iloop, islope ! loop counters
297real(dp)    :: alpha, lay1       ! coefficients for building layers
298real(dp)    :: index_powerlaw    ! coefficient to match the power low grid with the exponential one
299
300! CODE
301! ----
302! 1. Depth coordinate
303! -------------------
304! mlayer distribution: For the grid in common with the PCM: lay1*(1+k^(2.9)*(1-exp(-k/20))
305! Then we use a power low : lay1*alpha**(k-1/2)
306lay1 = 2.e-4
307alpha = 2
308do iloop = 0,nsoil_PCM - 1
309    mlayer(iloop) = lay1*(1 + iloop**2.9*(1._dp - exp(-real(iloop,dp)/20._dp)))
310end do
311
312do iloop = 0,nsoil - 1
313    if (lay1*(alpha**(iloop - 0.5_dp)) > mlayer(nsoil_PCM - 1)) then
314        index_powerlaw = iloop
315        exit
316    end if
317end do
318
319do iloop = nsoil_PCM,nsoil - 1
320    mlayer(iloop) = lay1*(alpha**(index_powerlaw + (iloop - nsoil_PCM) - 0.5_dp))
321end do
322
323! Build layer()
324do iloop = 1,nsoil - 1
325    layer(iloop) = (mlayer(iloop) + mlayer(iloop - 1))/2._dp
326end do
327layer(nsoil) = 2._dp*mlayer(nsoil - 1) - mlayer(nsoil - 2)
328
329
330! 2. Thermal inertia (note: it is declared in comsoil_h)
331! ------------------
332do ig = 1,ngrid
333    do islope = 1,nslope
334        do iloop = 1,nsoil_PCM
335            therm_initia(ig,iloop,islope) = TI_PCM(ig,iloop,islope)
336        end do
337        if (nsoil > nsoil_PCM) then
338            do iloop = nsoil_PCM + 1,nsoil
339                therm_initia(ig,iloop,islope) = TI_PCM(ig,nsoil_PCM,islope)
340            end do
341        end if
342    end do
343end do
344
345! 3. Index for subsurface layering
346! ------------------
347index_breccia = 1
348do iloop = 1,nsoil - 1
349    if (depth_breccia >= layer(iloop)) then
350        index_breccia = iloop
351    else
352        exit
353    end if
354end do
355
356index_bedrock = 1
357do iloop = 1,nsoil - 1
358    if (depth_bedrock >= layer(iloop)) then
359        index_bedrock = iloop
360    else
361        exit
362    end if
363end do
364
365END SUBROUTINE set_soil
366!=======================================================================
367
368!=======================================================================
369SUBROUTINE build4PCM_soil(tsoil_avg,tsoil_dev,inertiesoil4PCM,tsoil4PCM,flux_geo4PCM)
370!-----------------------------------------------------------------------
371! NAME
372!     build4PCM_soil
373!
374! DESCRIPTION
375!     Reconstructs the soil temperature profile, thermal inertia and
376!     the geothermal flux for the PCM.
377!
378! AUTHORS & DATE
379!     JB Clement, 12/2025
380!     C. Metz, 01/2026
381!
382! NOTES
383!
384!-----------------------------------------------------------------------
385
386! DEPENDENCIES
387! ------------
388use geometry, only: nsoil_PCM
389use display,  only: print_msg, LVL_NFO
390
391! DECLARATION
392! -----------
393implicit none
394
395! ARGUMENTS
396! ---------
397real(dp), dimension(:,:,:), intent(in)  :: tsoil_avg, tsoil_dev
398real(dp), dimension(:,:,:), intent(out) :: tsoil4PCM, inertiesoil4PCM
399real(dp), dimension(:,:),   intent(out) :: flux_geo4PCM
400
401! CODE
402! ----
403call print_msg('> Building soil thermal inertia for the PCM',LVL_NFO)
404inertiesoil4PCM(:,:,:) = TI(:,1:nsoil_PCM,:)
405
406call print_msg('> Building soil temperature profile for the PCM',LVL_NFO)
407tsoil4PCM(:,:,:) = tsoil_avg(:,1:nsoil_PCM,:) + tsoil_dev(:,:,:)
408
409call print_msg('> Building geothermal flux for the PCM',LVL_NFO)
410flux_geo4PCM(:,:) = (tsoil_avg(:,nsoil_PCM + 1,:) - tsoil4PCM(:,nsoil_PCM,:))/(mlayer(nsoil_PCM + 1) - mlayer(nsoil_PCM))
411
412END SUBROUTINE build4PCM_soil
413!=======================================================================
414
415END MODULE soil
Note: See TracBrowser for help on using the repository browser.