source: trunk/LMDZ.COMMON/libf/evolution/soil_temp.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: 14.8 KB
Line 
1MODULE soil_temp
2!-----------------------------------------------------------------------
3! NAME
4!     soil_temp
5!
6! DESCRIPTION
7!     Routines to compute soil temperature evolution and initialization.
8!
9! AUTHORS & DATE
10!     L. Lange, 2023
11!     JB Clement, 2023-2025
12!
13! NOTES
14!
15!-----------------------------------------------------------------------
16
17! DECLARATION
18! -----------
19implicit none
20
21contains
22!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
23
24!=======================================================================
25SUBROUTINE compute_tsoil(ngrid,nsoil,firstcall,therm_i,timestep,tsurf,tsoil)
26!-----------------------------------------------------------------------
27! NAME
28!     compute_tsoil
29!
30! DESCRIPTION
31!     Compute soil temperature using an implicit 1st order scheme.
32!
33! AUTHORS & DATE
34!     L. Lange, 2023
35!     JB Clement, 2023-2025
36!
37! NOTES
38!
39!-----------------------------------------------------------------------
40
41! DEPENDENCIES
42! ------------
43use soil,      only: layer_PEM, mlayer_PEM, mthermdiff_PEM, thermdiff_PEM, coefq_PEM, coefd_PEM, mu_PEM, alph_PEM, beta_PEM, fluxgeo
44use comsoil_h, only: volcapa
45
46! DECLARATION
47! -----------
48implicit none
49
50#include "dimensions.h"
51
52! ARGUMENTS
53! ---------
54integer,                      intent(in)    :: ngrid     ! number of (horizontal) grid-points
55integer,                      intent(in)    :: nsoil     ! number of soil layers
56logical,                      intent(in)    :: firstcall ! identifier for initialization call
57real, dimension(ngrid,nsoil), intent(in)    :: therm_i   ! thermal inertia [SI]
58real,                         intent(in)    :: timestep  ! time step [s]
59real, dimension(ngrid),       intent(in)    :: tsurf     ! surface temperature [K]
60real, dimension(ngrid,nsoil), intent(inout) :: tsoil     ! soil (mid-layer) temperature [K]
61
62! LOCAL VARIABLES
63! ---------------
64integer :: ig, ik
65
66! CODE
67! ----
68! 0. Initialisations and preprocessing step
69 if (firstcall) then
70! 0.1 Build mthermdiff_PEM(:), the mid-layer thermal diffusivities
71    do ig = 1,ngrid
72        do ik = 0,nsoil - 1
73            mthermdiff_PEM(ig,ik) = therm_i(ig,ik + 1)*therm_i(ig,ik + 1)/volcapa
74        enddo
75    enddo
76
77! 0.2 Build thermdiff(:), the "interlayer" thermal diffusivities
78    do ig = 1,ngrid
79        do ik = 1,nsoil - 1
80            thermdiff_PEM(ig,ik) = ((layer_PEM(ik) - mlayer_PEM(ik - 1))*mthermdiff_PEM(ig,ik) &
81                                   + (mlayer_PEM(ik) - layer_PEM(ik))*mthermdiff_PEM(ig,ik - 1))/(mlayer_PEM(ik) - mlayer_PEM(ik - 1))
82        enddo
83    enddo
84
85! 0.3 Build coefficients mu_PEM, q_{k+1/2}, d_k, alph_PEM
86    ! mu_PEM
87    mu_PEM = mlayer_PEM(0)/(mlayer_PEM(1) - mlayer_PEM(0))
88
89    ! q_{1/2}
90    coefq_PEM(0) = volcapa*layer_PEM(1)/timestep
91    ! q_{k+1/2}
92    do ik = 1,nsoil - 1
93        coefq_PEM(ik) = volcapa*(layer_PEM(ik + 1) - layer_PEM(ik))/timestep
94    enddo
95
96    do ig = 1,ngrid
97        ! d_k
98        do ik = 1,nsoil - 1
99            coefd_PEM(ig,ik) = thermdiff_PEM(ig,ik)/(mlayer_PEM(ik)-mlayer_PEM(ik - 1))
100        enddo
101
102        ! alph_PEM_{N-1}
103        alph_PEM(ig,nsoil - 1) = coefd_PEM(ig,nsoil-1)/(coefq_PEM(nsoil - 1) + coefd_PEM(ig,nsoil - 1))
104        ! alph_PEM_k
105        do ik = nsoil - 2,1,-1
106            alph_PEM(ig,ik) = coefd_PEM(ig,ik)/(coefq_PEM(ik) + coefd_PEM(ig,ik + 1)*(1. - alph_PEM(ig,ik + 1)) + coefd_PEM(ig,ik))
107        enddo
108      enddo ! of do ig=1,ngrid
109endif ! of if (firstcall)
110
111if (.not. firstcall) THEN
112! 2. Compute soil temperatures
113! First layer:
114    do ig = 1,ngrid
115        tsoil(ig,1) = (tsurf(ig) + mu_PEM*beta_PEM(ig,1)*thermdiff_PEM(ig,1)/mthermdiff_PEM(ig,0))/ &
116                      (1. + mu_PEM*(1. - alph_PEM(ig,1))*thermdiff_PEM(ig,1)/mthermdiff_PEM(ig,0))
117
118! Other layers:
119        do ik = 1,nsoil - 1
120                tsoil(ig,ik + 1) = alph_PEM(ig,ik)*tsoil(ig,ik) + beta_PEM(ig,ik)
121        enddo
122    enddo
123endif
124
125!  2. Compute beta_PEM coefficients (preprocessing for next time step)
126! Bottom layer, beta_PEM_{N-1}
127do ig = 1,ngrid
128    beta_PEM(ig,nsoil - 1) = coefq_PEM(nsoil - 1)*tsoil(ig,nsoil)/(coefq_PEM(nsoil - 1) + coefd_PEM(ig,nsoil - 1)) &
129                             + fluxgeo/(coefq_PEM(nsoil - 1) + coefd_PEM(ig,nsoil - 1))
130enddo
131! Other layers
132do ik = nsoil-2,1,-1
133    do ig = 1,ngrid
134        beta_PEM(ig,ik) = (coefq_PEM(ik)*tsoil(ig,ik + 1) + coefd_PEM(ig,ik + 1)*beta_PEM(ig,ik + 1))/ &
135                          (coefq_PEM(ik) + coefd_PEM(ig,ik + 1)*(1. - alph_PEM(ig,ik + 1)) + coefd_PEM(ig,ik))
136    enddo
137enddo
138
139END SUBROUTINE compute_tsoil
140!=======================================================================
141
142!=======================================================================
143SUBROUTINE ini_tsoil_pem(ngrid,nsoil,therm_i,tsurf,tsoil)
144!-----------------------------------------------------------------------
145! NAME
146!     ini_tsoil_pem
147!
148! DESCRIPTION
149!     Initialize soil with the solution of the stationary heat conduction problem.
150!     Boundary conditions: Tsurf averaged from PCM; geothermal flux at bottom.
151!
152! AUTHORS & DATE
153!     L. Lang, 2023
154!     JB Clement, 2023-2025
155!
156! NOTES
157!
158!-----------------------------------------------------------------------
159
160! DEPENDENCIES
161! ------------
162use soil,      only: layer_PEM, mlayer_PEM, mthermdiff_PEM, thermdiff_PEM, coefq_PEM, coefd_PEM, mu_PEM, alph_PEM, beta_PEM, fluxgeo
163use comsoil_h, only: volcapa
164
165! DECLARATION
166! -----------
167implicit none
168
169#include "dimensions.h"
170
171! ARGUMENTS
172! ---------
173integer,                      intent(in)    :: ngrid   ! number of (horizontal) grid-points
174integer,                      intent(in)    :: nsoil   ! number of soil layers
175real, dimension(ngrid,nsoil), intent(in)    :: therm_i ! thermal inertia [SI]
176real, dimension(ngrid),       intent(in)    :: tsurf   ! surface temperature [K]
177real, dimension(ngrid,nsoil), intent(inout) :: tsoil   ! soil (mid-layer) temperature [K]
178
179! LOCAL VARIABLES
180! ---------------
181integer :: ig, ik, iloop
182
183! CODE
184! ----
185! 0. Initialisations and preprocessing step
186! 0.1 Build mthermdiff_PEM(:), the mid-layer thermal diffusivities
187do ig = 1,ngrid
188    do ik = 0,nsoil - 1
189        mthermdiff_PEM(ig,ik) = therm_i(ig,ik + 1)*therm_i(ig,ik + 1)/volcapa
190    enddo
191enddo
192
193! 0.2 Build thermdiff(:), the "interlayer" thermal diffusivities
194do ig = 1,ngrid
195    do ik = 1,nsoil - 1
196        thermdiff_PEM(ig,ik) = ((layer_PEM(ik) - mlayer_PEM(ik - 1))*mthermdiff_PEM(ig,ik) &
197                             + (mlayer_PEM(ik) - layer_PEM(ik))*mthermdiff_PEM(ig,ik - 1))/(mlayer_PEM(ik) - mlayer_PEM(ik - 1))
198    enddo
199enddo
200
201! 0.3 Build coefficients mu_PEM, q_{k+1/2}, d_k, alph_PEM
202! mu_PEM
203mu_PEM = mlayer_PEM(0)/(mlayer_PEM(1) - mlayer_PEM(0))
204
205! q_{1/2}
206coefq_PEM(:) = 0.
207! q_{k+1/2}
208
209do ig = 1,ngrid
210    ! d_k
211    do ik = 1,nsoil - 1
212        coefd_PEM(ig,ik) = thermdiff_PEM(ig,ik)/(mlayer_PEM(ik) - mlayer_PEM(ik - 1))
213    enddo
214
215    ! alph_PEM_{N-1}
216    alph_PEM(ig,nsoil - 1) = coefd_PEM(ig,nsoil - 1)/(coefq_PEM(nsoil - 1) + coefd_PEM(ig,nsoil - 1))
217    ! alph_PEM_k
218    do ik = nsoil - 2,1,-1
219        alph_PEM(ig,ik) = coefd_PEM(ig,ik)/(coefq_PEM(ik) + coefd_PEM(ig,ik + 1)*(1. - alph_PEM(ig,ik + 1)) + coefd_PEM(ig,ik))
220    enddo
221enddo ! of do ig=1,ngrid
222
223!  1. Compute beta_PEM coefficients
224! Bottom layer, beta_PEM_{N-1}
225do ig = 1,ngrid
226        beta_PEM(ig,nsoil - 1) = coefq_PEM(nsoil - 1)*tsoil(ig,nsoil)/(coefq_PEM(nsoil - 1) + coefd_PEM(ig,nsoil - 1)) &
227                                 + fluxgeo/(coefq_PEM(nsoil - 1) + coefd_PEM(ig,nsoil - 1))
228enddo
229! Other layers
230do ik = nsoil - 2,1,-1
231    do ig = 1,ngrid
232        beta_PEM(ig,ik) = (coefq_PEM(ik)*tsoil(ig,ik + 1) + coefd_PEM(ig,ik+1)*beta_PEM(ig,ik + 1))/ &
233                          (coefq_PEM(ik) + coefd_PEM(ig,ik + 1)*(1. - alph_PEM(ig,ik + 1)) + coefd_PEM(ig,ik))
234    enddo
235enddo
236
237!  2. Compute soil temperatures
238do iloop = 1,10 !just convergence
239    do ig = 1,ngrid
240        ! First layer:
241        tsoil(ig,1) = (tsurf(ig) + mu_PEM*beta_PEM(ig,1)*thermdiff_PEM(ig,1)/mthermdiff_PEM(ig,0))/ &
242                      (1. + mu_PEM*(1. - alph_PEM(ig,1))*thermdiff_PEM(ig,1)/mthermdiff_PEM(ig,0))
243        ! Other layers:
244        do ik = 1,nsoil - 1
245            tsoil(ig,ik + 1) = alph_PEM(ig,ik)*tsoil(ig,ik) + beta_PEM(ig,ik)
246        enddo
247    enddo
248enddo ! iloop
249
250END SUBROUTINE ini_tsoil_pem
251!=======================================================================
252
253!=======================================================================
254SUBROUTINE shift_tsoil2surf(ngrid,nsoil,nslope,zshift_surf,zlag,tsurf,tsoil)
255!-----------------------------------------------------------------------
256! NAME
257!     shift_tsoil2surf
258!
259! DESCRIPTION
260!     Shift soil temperature profile to follow surface evolution due to ice
261!     condensation/sublimation.
262!
263! AUTHORS & DATE
264!     JB Clement, 2025
265!
266! NOTES
267!
268!-----------------------------------------------------------------------
269
270! DEPENDENCIES
271! ------------
272use soil,         only: layer_PEM, mlayer_PEM, fluxgeo, thermdiff_PEM, mthermdiff_PEM
273use math_toolkit, only: solve_steady_heat
274
275! DECLARATION
276! -----------
277implicit none
278
279! ARGUMENTS
280! ---------
281integer,                             intent(in)    :: ngrid       ! number of (horizontal) grid-points
282integer,                             intent(in)    :: nsoil       ! number of soil layers
283integer,                             intent(in)    :: nslope      ! number of sub-slopes
284real, dimension(ngrid,nslope),       intent(in)    :: zshift_surf ! elevation shift for the surface [m]
285real, dimension(ngrid,nslope),       intent(in)    :: zlag        ! newly built lag thickness [m]
286real, dimension(ngrid,nslope),       intent(in)    :: tsurf       ! surface temperature [K]
287real, dimension(ngrid,nsoil,nslope), intent(inout) :: tsoil       ! soil (mid-layer) temperature [K]
288
289! LOCAL VARIABLES
290! ---------------
291integer                             :: ig, isoil, islope, ishift, ilag, iz
292real                                :: a, z, zshift_surfloc, tsoil_minus, mlayer_minus
293real, dimension(ngrid,nsoil,nslope) :: tsoil_old
294
295! CODE
296! ----
297write(*,*) "> Shifting soil temperature profile to match surface evolution"
298tsoil_old = tsoil
299
300do ig = 1,ngrid
301    do islope = 1,nslope
302        zshift_surfloc = zshift_surf(ig,islope)
303
304        if (zshift_surfloc >= 0.) then ! In case of the surface is higher than initially
305            if (zshift_surfloc < mlayer_PEM(0)) then ! Surface change is too small to be taken into account
306                ! Nothing to do; we keep the soil temperature profile
307            else if (zshift_surfloc >= mlayer_PEM(nsoil - 1)) then ! Surface change is much larger than the discretization
308                tsoil(ig,:,islope) = tsurf(ig,islope)
309            else
310                ishift = 0
311                do while (mlayer_PEM(ishift) <= zshift_surfloc)
312                    ishift = ishift + 1 ! mlayer indices begin at 0 so this the good index for tsoil!
313                enddo
314                ! The "new soil" temperature is set to tsurf
315                tsoil(ig,:ishift,islope) = tsurf(ig,islope)
316                do isoil = ishift + 1,nsoil
317                    ! Position in the old discretization of the depth
318                    z = mlayer_PEM(isoil - 1) - zshift_surfloc
319                    ! Interpolation of the temperature profile from the old discretization
320                    tsoil(ig,isoil,islope) = itp_tsoil(tsoil_old(ig,:,islope),tsurf(ig,islope),z)
321                enddo
322            endif
323        else ! In case of the surface is lower than initially
324            if (abs(zshift_surfloc) < mlayer_PEM(0)) then ! Surface change is too small to be taken into account
325                ! Nothing to do; we keep the soil temperature profile
326            else if (abs(zshift_surfloc) >= mlayer_PEM(nsoil - 1)) then ! Surface change is much larger than the discretization
327                call solve_steady_heat(nsoil,mlayer_PEM,layer_PEM,mthermdiff_PEM(ig,:),thermdiff_PEM(ig,:),tsurf(ig,islope),fluxgeo,tsoil(ig,:,islope))
328            else
329                if (zlag(ig,islope) < mlayer_PEM(0)) then ! The lag is too thin to be taken into account
330                    ilag = 0
331                else
332                    ilag = 0
333                    do while (mlayer_PEM(ilag) <= zlag(ig,islope))
334                        ilag = ilag + 1 ! mlayer indices begin at 0 so this the good index for tsoil!
335                    enddo
336                    ! Position of the lag bottom in the old discretization of the depth
337                    z = zlag(ig,islope) - zshift_surfloc
338                    ! The "new lag" temperature is set to the ice temperature just below
339                    tsoil(ig,:ilag,islope) = itp_tsoil(tsoil_old(ig,:,islope),tsurf(ig,islope),z)
340                endif
341
342                ishift = nsoil - 1
343                z = mlayer_PEM(nsoil - 1) + zshift_surfloc
344                do while (mlayer_PEM(ishift) >= z)
345                    ishift = ishift - 1
346                enddo
347                ishift = ishift + 1 ! Adding 1 is needed to match the good index for tsoil!
348                do isoil = ilag + 1,ishift
349                    ! Position in the old discretization of the depth
350                    z = mlayer_PEM(isoil - 1) - zshift_surfloc
351                    ! Interpolation of the temperature profile from the old discretization
352                    tsoil(ig,isoil,islope) = itp_tsoil(tsoil_old(ig,:,islope),tsurf(ig,islope),z)
353                enddo
354
355                ! The "new deepest layers" temperature is set by solving the steady heat equation
356                call solve_steady_heat(nsoil - ishift + 1,mlayer_PEM(ishift - 1:),layer_PEM(ishift:),mthermdiff_PEM(ig,ishift - 1:),thermdiff_PEM(ig,ishift:),tsoil(ig,ishift,islope),fluxgeo,tsoil(ig,ishift:,islope))
357            endif
358        endif
359    enddo
360enddo
361
362END SUBROUTINE shift_tsoil2surf
363!=======================================================================
364
365!=======================================================================
366FUNCTION itp_tsoil(tsoil,tsurf,z) RESULT(tsoil_z)
367!-----------------------------------------------------------------------
368! NAME
369!     itp_tsoil
370!
371! DESCRIPTION
372!     Interpolate soil temperature profile.
373!
374! AUTHORS & DATE
375!     JB Clement, 2025
376!
377! NOTES
378!
379!-----------------------------------------------------------------------
380
381! DEPENDENCIES
382! ------------
383use soil, only: mlayer_PEM
384
385! DECLARATION
386! -----------
387implicit none
388
389! ARGUMENTS
390! ---------
391real, dimension(:), intent(in) :: tsoil
392real,               intent(in) :: z, tsurf
393
394! LOCAL VARIABLES
395! ---------------
396real    :: tsoil_z, tsoil_minus, mlayer_minus, a
397integer :: iz
398
399! CODE
400! ----
401! Find the interval [mlayer_PEM(iz - 1),mlayer_PEM(iz)[ where the position z belongs
402iz = 0
403do while (mlayer_PEM(iz) <= z)
404    iz = iz + 1
405enddo
406if (iz == 0) then
407    tsoil_minus = tsurf
408    mlayer_minus = 0.
409else
410    tsoil_minus = tsoil(iz)
411    mlayer_minus = mlayer_PEM(iz - 1)
412endif
413! Interpolation of the temperature profile from the old discretization
414a = (tsoil(iz + 1) - tsoil_minus)/(mlayer_PEM(iz) - mlayer_minus)
415tsoil_z = a*(z - mlayer_minus) + tsoil_minus
416
417END FUNCTION itp_tsoil
418!=======================================================================
419
420END MODULE soil_temp
Note: See TracBrowser for help on using the repository browser.