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

Last change on this file since 3552 was 3532, checked in by jbclement, 3 weeks ago

PEM:
Removing unecessary module/subroutine "interpol_TI_PEM2PCM.F90" + Few small corrections/cleanings.
JBC
File size: 7.9 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
218END MODULE compute_soiltemp_mod
Note: See TracBrowser for help on using the repository browser.