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