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

Last change on this file since 4074 was 4074, checked in by jbclement, 11 days ago

PEM:

  • Correct management of H2O ice tendency in 1D when there is not enough ice anymore.
  • Clean initialization of allocatable module arrays (especially needed when no slope)
  • One more renaming for consistency + few small updates thoughout the code.

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