source: trunk/LMDZ.COMMON/libf/evolution/compute_soiltemp_mod.F90 @ 3705

Last change on this file since 3705 was 3666, checked in by jbclement, 9 months ago

PEM:

  • Adjusments of few scripts to handle more situations.
  • Bug correction related to deallocation in case of "soil_pem = .false." or "layering_algo = .true.".

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