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

Last change on this file since 3525 was 3525, checked in by jbclement, 2 days ago

PEM:
Computation of <soil thermal intertia> and <H2O mass subsurface/surface exchange> according to the presence of subsurface ice provided by the (Norbert's) dynamic method + few cleanings.
JBC
File size: 8.1 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
14
15
16SUBROUTINE compute_tsoil_pem(ngrid,nsoil,firstcall,therm_i,timestep,tsurf,tsoil)
17
18use comsoil_h_PEM, only: layer_PEM, mlayer_PEM, mthermdiff_PEM, thermdiff_PEM, coefq_PEM, coefd_PEM, mu_PEM, alph_PEM, beta_PEM, fluxgeo
19use comsoil_h,     only: volcapa
20
21implicit none
22
23!-----------------------------------------------------------------------
24!  Author: LL
25!  Purpose: Compute soil temperature using an implict 1st order scheme
26! 
27!  Note: depths of layers and mid-layers, soil thermal inertia and
28!        heat capacity are commons in comsoil_PEM.h
29!-----------------------------------------------------------------------
30
31#include "dimensions.h"
32
33!-----------------------------------------------------------------------
34!  arguments
35!  ---------
36!  inputs:
37integer,                      intent(in) :: ngrid     ! number of (horizontal) grid-points
38integer,                      intent(in) :: nsoil     ! number of soil layers 
39logical,                      intent(in) :: firstcall ! identifier for initialization call
40real, dimension(ngrid,nsoil), intent(in) :: therm_i   ! thermal inertia [SI]
41real,                         intent(in) :: timestep  ! time step [s]
42real, dimension(ngrid),       intent(in) :: tsurf     ! surface temperature [K]
43 
44! outputs:
45real, dimension(ngrid,nsoil), intent(inout) :: tsoil ! soil (mid-layer) temperature [K]
46! local variables:
47integer :: ig, ik   
48
49! 0. Initialisations and preprocessing step
50 if (firstcall) then
51! 0.1 Build mthermdiff_PEM(:), the mid-layer thermal diffusivities
52    do ig = 1,ngrid
53        do ik = 0,nsoil - 1
54            mthermdiff_PEM(ig,ik) = therm_i(ig,ik + 1)*therm_i(ig,ik + 1)/volcapa   
55        enddo
56    enddo
57
58! 0.2 Build thermdiff(:), the "interlayer" thermal diffusivities
59    do ig = 1,ngrid
60        do ik = 1,nsoil - 1
61            thermdiff_PEM(ig,ik) = ((layer_PEM(ik) - mlayer_PEM(ik - 1))*mthermdiff_PEM(ig,ik) &
62                                   + (mlayer_PEM(ik) - layer_PEM(ik))*mthermdiff_PEM(ig,ik - 1))/(mlayer_PEM(ik) - mlayer_PEM(ik - 1))
63        enddo
64    enddo
65
66! 0.3 Build coefficients mu_PEM, q_{k+1/2}, d_k, alph_PEMa_k and capcal
67    ! mu_PEM
68    mu_PEM = mlayer_PEM(0)/(mlayer_PEM(1) - mlayer_PEM(0))
69
70    ! q_{1/2}
71    coefq_PEM(0) = volcapa*layer_PEM(1)/timestep
72    ! q_{k+1/2}
73    do ik = 1,nsoil - 1
74        coefq_PEM(ik) = volcapa*(layer_PEM(ik + 1) - layer_PEM(ik))/timestep
75    enddo
76
77    do ig = 1,ngrid
78        ! d_k
79        do ik = 1,nsoil - 1
80            coefd_PEM(ig,ik) = thermdiff_PEM(ig,ik)/(mlayer_PEM(ik)-mlayer_PEM(ik - 1))
81        enddo
82
83        ! alph_PEM_{N-1}
84        alph_PEM(ig,nsoil - 1) = coefd_PEM(ig,nsoil-1)/(coefq_PEM(nsoil - 1) + coefd_PEM(ig,nsoil - 1))
85        ! alph_PEM_k
86        do ik = nsoil - 2,1,-1
87            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))
88        enddo
89      enddo ! of do ig=1,ngrid
90endif ! of if (firstcall)
91
92if (.not. firstcall) THEN
93! 2. Compute soil temperatures
94! First layer:
95    do ig = 1,ngrid
96        tsoil(ig,1) = (tsurf(ig) + mu_PEM*beta_PEM(ig,1)*thermdiff_PEM(ig,1)/mthermdiff_PEM(ig,0))/ &
97                      (1. + mu_PEM*(1. - alph_PEM(ig,1))*thermdiff_PEM(ig,1)/mthermdiff_PEM(ig,0))
98
99! Other layers:
100        do ik = 1,nsoil - 1
101                tsoil(ig,ik + 1) = alph_PEM(ig,ik)*tsoil(ig,ik) + beta_PEM(ig,ik)
102        enddo
103    enddo
104endif
105
106!  2. Compute beta_PEM coefficients (preprocessing for next time step)
107! Bottom layer, beta_PEM_{N-1}
108do ig = 1,ngrid
109    beta_PEM(ig,nsoil - 1) = coefq_PEM(nsoil - 1)*tsoil(ig,nsoil)/(coefq_PEM(nsoil - 1) + coefd_PEM(ig,nsoil - 1)) &
110                             + fluxgeo/(coefq_PEM(nsoil - 1) + coefd_PEM(ig,nsoil - 1))
111enddo
112! Other layers
113do ik = nsoil-2,1,-1
114    do ig = 1,ngrid
115        beta_PEM(ig,ik) = (coefq_PEM(ik)*tsoil(ig,ik + 1) + coefd_PEM(ig,ik + 1)*beta_PEM(ig,ik + 1))/ &
116                          (coefq_PEM(ik) + coefd_PEM(ig,ik + 1)*(1. - alph_PEM(ig,ik + 1)) + coefd_PEM(ig,ik))
117    enddo
118enddo
119
120END SUBROUTINE compute_tsoil_pem
121
122!=======================================================================
123
124SUBROUTINE ini_tsoil_pem(ngrid,nsoil,therm_i,tsurf,tsoil)
125
126use comsoil_h_PEM, only: layer_PEM, mlayer_PEM, mthermdiff_PEM, thermdiff_PEM, coefq_PEM, coefd_PEM, mu_PEM, alph_PEM, beta_PEM, fluxgeo
127use comsoil_h,     only: volcapa
128
129implicit none
130
131!-----------------------------------------------------------------------
132!  Author: LL
133!  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
134!
135!  Note: depths of layers and mid-layers, soil thermal inertia and
136!        heat capacity are commons in comsoil_PEM.h
137!-----------------------------------------------------------------------
138
139#include "dimensions.h"
140
141!-----------------------------------------------------------------------
142!  arguments
143!  ---------
144!  inputs:
145integer,                      intent(in) :: ngrid   ! number of (horizontal) grid-points
146integer,                      intent(in) :: nsoil   ! number of soil layers
147real, dimension(ngrid,nsoil), intent(in) :: therm_i ! thermal inertia [SI]
148real, dimension(ngrid),       intent(in) :: tsurf   ! surface temperature [K]
149
150! outputs:
151real, dimension(ngrid,nsoil), intent(inout) :: tsoil ! soil (mid-layer) temperature [K]
152! local variables:
153integer :: ig, ik, iloop
154
155! 0. Initialisations and preprocessing step
156! 0.1 Build mthermdiff_PEM(:), the mid-layer thermal diffusivities
157do ig = 1,ngrid
158    do ik = 0,nsoil - 1
159        mthermdiff_PEM(ig,ik) = therm_i(ig,ik + 1)*therm_i(ig,ik + 1)/volcapa
160    enddo
161enddo
162
163! 0.2 Build thermdiff(:), the "interlayer" thermal diffusivities
164do ig = 1,ngrid
165    do ik = 1,nsoil - 1
166        thermdiff_PEM(ig,ik) = ((layer_PEM(ik) - mlayer_PEM(ik - 1))*mthermdiff_PEM(ig,ik) &
167                             + (mlayer_PEM(ik) - layer_PEM(ik))*mthermdiff_PEM(ig,ik - 1))/(mlayer_PEM(ik) - mlayer_PEM(ik - 1))
168    enddo
169enddo
170
171! 0.3 Build coefficients mu_PEM, q_{k+1/2}, d_k, alph_PEMa_k and capcal
172! mu_PEM
173mu_PEM = mlayer_PEM(0)/(mlayer_PEM(1) - mlayer_PEM(0))
174
175! q_{1/2}
176coefq_PEM(:) = 0.
177! q_{k+1/2}
178
179do ig = 1,ngrid
180    ! d_k
181    do ik = 1,nsoil-1
182        coefd_PEM(ig,ik) = thermdiff_PEM(ig,ik)/(mlayer_PEM(ik) - mlayer_PEM(ik - 1))
183    enddo
184
185    ! alph_PEM_{N-1}
186    alph_PEM(ig,nsoil - 1) = coefd_PEM(ig,nsoil - 1)/(coefq_PEM(nsoil - 1) + coefd_PEM(ig,nsoil - 1))
187    ! alph_PEM_k
188    do ik = nsoil - 2,1,-1
189        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))
190    enddo
191enddo ! of do ig=1,ngrid
192
193!  1. Compute beta_PEM coefficients
194! Bottom layer, beta_PEM_{N-1}
195do ig = 1,ngrid
196        beta_PEM(ig,nsoil - 1) = coefq_PEM(nsoil - 1)*tsoil(ig,nsoil)/(coefq_PEM(nsoil - 1) + coefd_PEM(ig,nsoil - 1)) &
197                                 + fluxgeo/(coefq_PEM(nsoil - 1) + coefd_PEM(ig,nsoil - 1))
198enddo
199! Other layers
200do ik = nsoil - 2,1,-1
201    do ig = 1,ngrid
202        beta_PEM(ig,ik) = (coefq_PEM(ik)*tsoil(ig,ik + 1) + coefd_PEM(ig,ik+1)*beta_PEM(ig,ik + 1))/ &
203                          (coefq_PEM(ik) + coefd_PEM(ig,ik + 1)*(1. - alph_PEM(ig,ik + 1)) + coefd_PEM(ig,ik))
204    enddo
205enddo
206
207!  2. Compute soil temperatures
208do iloop = 1,10 !just convergence
209    ! First layer:
210    do ig = 1,ngrid
211        tsoil(ig,1)=(tsurf(ig) + mu_PEM*beta_PEM(ig,1)*thermdiff_PEM(ig,1)/mthermdiff_PEM(ig,0))/ &
212                   (1. + mu_PEM*(1. - alph_PEM(ig,1))*thermdiff_PEM(ig,1)/mthermdiff_PEM(ig,0))
213    ! Other layers:
214    do ik = 1,nsoil - 1
215        tsoil(ig,ik + 1) = alph_PEM(ig,ik)*tsoil(ig,ik) + beta_PEM(ig,ik)
216    enddo
217enddo
218
219enddo ! iloop
220
221END SUBROUTINE ini_tsoil_pem
222
223END MODULE compute_soiltemp_mod
Note: See TracBrowser for help on using the repository browser.