1 | ! ================================================================================================================================= |
---|
2 | ! MODULE : thermosoil |
---|
3 | ! |
---|
4 | ! CONTACT : orchidee-help _at_ ipsl.jussieu.fr |
---|
5 | ! |
---|
6 | ! LICENCE : IPSL (2006) |
---|
7 | ! This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC |
---|
8 | ! |
---|
9 | !>\BRIEF Calculates the soil temperatures by solving the heat |
---|
10 | !! diffusion equation within the soil |
---|
11 | !! |
---|
12 | !!\n DESCRIPTION : General important informations about the numerical scheme and |
---|
13 | !! the soil vertical discretization:\n |
---|
14 | !! - the soil is divided into "ngrnd" (=7 by default) layers, reaching to as |
---|
15 | !! deep as 5.5m down within the soil, with thiscknesses |
---|
16 | !! following a geometric series of ration 2.\n |
---|
17 | !! - "jg" is usually used as the index going from 1 to ngrnd to describe the |
---|
18 | !! layers, from top (jg=1) to bottom (jg=ngrnd)\n |
---|
19 | !! - the thermal numerical scheme is implicit finite differences.\n |
---|
20 | !! -- When it is resolved in thermosoil_profile at the present timestep t, the |
---|
21 | !! dependancy from the previous timestep (t-1) is hidden in the |
---|
22 | !! integration coefficients cgrnd and dgrnd, which are therefore |
---|
23 | !! calculated at the very end of thermosoil_main (call to |
---|
24 | !! thermosoil_coef) for use in the next timestep.\n |
---|
25 | !! -- At timestep t, the system becomes :\n |
---|
26 | !! |
---|
27 | !! T(k+1)=cgrnd(k)+dgrnd(k)*T(k) \n |
---|
28 | !! -- EQ1 -- \n |
---|
29 | !! |
---|
30 | !! (the bottom boundary condition has been used to obtained this equation).\n |
---|
31 | !! To solve it, the uppermost soil temperature T(1) is required. |
---|
32 | !! It is obtained from the surface temperature Ts, which is |
---|
33 | !! considered a linear extrapolation of T(1) and T(2)\n |
---|
34 | !! |
---|
35 | !! Ts=(1-lambda)*T(1) -lambda*T(2) \n |
---|
36 | !! -- EQ2--\n |
---|
37 | !! |
---|
38 | !! -- caveat 1 : Ts is called 'temp_soil_new' in this routine, |
---|
39 | !! don' t act.\n |
---|
40 | !! -- caveat 2 : actually, the surface temperature at time t Ts |
---|
41 | !! depends on the soil temperature at time t through the |
---|
42 | !! ground heat flux. This is again implicitly solved, with Ts(t) |
---|
43 | !! expressed as :\n |
---|
44 | !! |
---|
45 | !! soilcap*(Ts(t)-Ts(t-1))/dt=soilflux+otherfluxes(Ts(t))\n |
---|
46 | !! -- EQ3 --\n |
---|
47 | !! |
---|
48 | !! and the dependency from the previous timestep is hidden in |
---|
49 | !! soilcap and soilflux (apparent surface heat capacity and heat |
---|
50 | !! flux respectively). Soilcap and soilflux are therefore |
---|
51 | !! calculated at the previsou timestep, at the very end of thermosoil |
---|
52 | !! (final call to thermosoil_coef) and stored to be used at the next time step. |
---|
53 | !! At timestep t, EQ3 is solved for Ts in enerbil, and Ts |
---|
54 | !! is used in thermosoil to get T(1) and solve EQ1.\n |
---|
55 | !! |
---|
56 | !! - lambda is the @tex $\mu$ @endtex of F. Hourdin' s PhD thesis, equation (A28); ie the |
---|
57 | !! coefficient of the linear extrapolation of Ts (surface temperature) from T1 and T2 (ptn(jg=1) and ptn(jg=2)), so that:\n |
---|
58 | !! Ts= (1+lambda)*T(1)-lambda*T(2) --EQ2-- \n |
---|
59 | !! lambda = (zz_coeff(1))/((zz_coef(2)-zz_coef(1))) \n |
---|
60 | !! |
---|
61 | !! - cstgrnd is the attenuation depth of the diurnal temperature signal |
---|
62 | !! (period : one_day) as a result of the heat conduction equation |
---|
63 | !! with no coefficients : |
---|
64 | !!\latexonly |
---|
65 | !!\input{thermosoil_var_init0.tex} |
---|
66 | !!\endlatexonly |
---|
67 | !! -- EQ4 --\n |
---|
68 | !! This equation results from the change of variables : |
---|
69 | !! z' =z*sqrt(Cp/K) where z' is the new depth (homogeneous |
---|
70 | !! to sqrt(time) ), z the real depth (in m), Cp and K the soil heat |
---|
71 | !! capacity and conductivity respectively.\n |
---|
72 | !! |
---|
73 | !! the attenuation depth of a diurnal thermal signal for EQ4 is therefore homogeneous to sqrt(time) and |
---|
74 | !! equals : \n |
---|
75 | !! cstgrnd = sqrt(oneday/Pi) |
---|
76 | !! |
---|
77 | !! - lskin is the attenuation depth of the diurnal temperature signal |
---|
78 | !! (period : one_day) within the soil for the complete heat conduction equation |
---|
79 | !! (ie : with coefficients) |
---|
80 | !!\latexonly |
---|
81 | !!\input{thermosoil_var_init00.tex} |
---|
82 | !!\endlatexonly |
---|
83 | !! -- EQ5 -- \n |
---|
84 | !! it can be retrieved from cstgrnd using the change of variable z' =z*sqrt(Cp/K):\n |
---|
85 | !! lskin = sqrt(K/Cp)*cstgrnd = sqrt(K/Cp)*sqrt(oneday//Pi)\n |
---|
86 | !! |
---|
87 | !! In thermosoil, the ratio lskin/cstgrnd is frequently used as the |
---|
88 | !! multiplicative factor to go from |
---|
89 | !!'adimensional' depths (like z' ) to real depths (z). z' is not really |
---|
90 | !! adimensional but is reffered to like this in the code. |
---|
91 | !! |
---|
92 | !! |
---|
93 | !! RECENT CHANGE(S) : None |
---|
94 | !! |
---|
95 | !! REFERENCE(S) : None |
---|
96 | !! |
---|
97 | !! SVN : |
---|
98 | !! $HeadURL$ |
---|
99 | !! $Date$ |
---|
100 | !! $Revision$ |
---|
101 | !! \n |
---|
102 | !_ ================================================================================================================================ |
---|
103 | |
---|
104 | MODULE thermosoil |
---|
105 | |
---|
106 | USE ioipsl |
---|
107 | |
---|
108 | ! modules used : |
---|
109 | USE constantes |
---|
110 | USE constantes_soil |
---|
111 | USE sechiba_io |
---|
112 | USE grid |
---|
113 | USE parallel |
---|
114 | |
---|
115 | IMPLICIT NONE |
---|
116 | |
---|
117 | !private and public routines : |
---|
118 | PRIVATE |
---|
119 | PUBLIC :: thermosoil_main,thermosoil_clear |
---|
120 | |
---|
121 | LOGICAL, SAVE :: l_first_thermosoil=.TRUE.!! does the initialisation of the routine |
---|
122 | !! (true/false) |
---|
123 | CHARACTER(LEN=80) , SAVE :: var_name !! To store variables names for the |
---|
124 | !! input-outputs dealt with by IOIPSL |
---|
125 | REAL(r_std), SAVE :: lambda, cstgrnd, lskin !! See Module description |
---|
126 | REAL(r_std), SAVE :: fz1, zalph !! usefull constants for diverse use |
---|
127 | REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: ptn !! vertically discretized |
---|
128 | !! soil temperatures @tex ($K$) @endtex. |
---|
129 | REAL(r_std), SAVE, DIMENSION (ngrnd) :: zz !! depths of the soil thermal numerical nodes. |
---|
130 | !! Caveats: they are not exactly the centers of the |
---|
131 | !! thermal layers, see the calculation in |
---|
132 | !! ::thermosoil_var_init @tex ($m$) @endtex. |
---|
133 | REAL(r_std), SAVE, DIMENSION (ngrnd) :: zz_coef !! depths of the boundaries of the thermal layers, |
---|
134 | !! see the calculation in |
---|
135 | !! thermosoil_var_init @tex ($m$) @endtex. |
---|
136 | REAL(r_std), SAVE, DIMENSION (ngrnd) :: dz1 !! numerical constant used in the thermal numerical |
---|
137 | !! scheme @tex ($m^{-1}$) @endtex. ; it corresponds |
---|
138 | !! to the coefficient @tex $d_k$ @endtex of equation |
---|
139 | !! (A.12) in F. Hourdin PhD thesis. |
---|
140 | REAL(r_std), SAVE, DIMENSION (ngrnd) :: dz2 !! thicknesses of the thermal layers @tex ($m$) |
---|
141 | !! @endtex; typically: |
---|
142 | !! dz2(jg)=zz_coef(jg+1)-zz_coef(jg); calculated once |
---|
143 | !! and for all in thermosoil_var_init |
---|
144 | REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: z1 !! constant of the numerical scheme; it is an |
---|
145 | !! intermediate buffer for the calculation of the |
---|
146 | !! integration coefficients cgrnd and dgrnd. |
---|
147 | REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: cgrnd !! integration coefficient for the numerical scheme, |
---|
148 | !! see eq.1 |
---|
149 | REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: dgrnd !! integration coefficient for the numerical scheme, |
---|
150 | !! see eq.1 |
---|
151 | REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: pcapa !! volumetric vertically discretized soil heat |
---|
152 | !! capacity @tex ($J K^{-1} m^{-3}$) @endtex. |
---|
153 | !! It depends on the soil |
---|
154 | !! moisture content (wetdiag) and is calculated at |
---|
155 | !! each time step in thermosoil_coef. |
---|
156 | REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: pkappa !! vertically discretized soil thermal conductivity |
---|
157 | !! @tex ($W K^{-1} m^{-1}$) @endtex. Same as pcapa. |
---|
158 | REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: zdz1 !! numerical constant of the numerical scheme; it is |
---|
159 | !! an intermediate buffer for the calculation of the |
---|
160 | !! integration coefficients cgrnd and dgrnd |
---|
161 | !! @tex ($W K^{-1} m^{-1}$) @endtex |
---|
162 | REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: zdz2 !! numerical constant of the numerical scheme; it is |
---|
163 | !! an intermediate buffer for the calculation of the |
---|
164 | !! integration coefficients cgrnd and dgrnd |
---|
165 | !! @tex ($W K^{-1} m^{-1}$) @endtex |
---|
166 | REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: pcapa_en !! heat capacity used for surfheat_incr and |
---|
167 | !! coldcont_incr |
---|
168 | REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: ptn_beg !! vertically discretized temperature at the |
---|
169 | !! beginning of the time step @tex ($K$) @endtex; |
---|
170 | !! is used in |
---|
171 | !! thermosoil_energy for energy-related diagnostic of |
---|
172 | !! the routine. |
---|
173 | REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: temp_sol_beg !! Surface temperature at the beginning of the |
---|
174 | !! timestep @tex ($K$) @endtex |
---|
175 | REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: surfheat_incr !! Change in soil heat content during the timestep |
---|
176 | !! @tex ($J$) @endtex. |
---|
177 | REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: coldcont_incr !! Change in snow heat content @tex ($J$) @endtex. |
---|
178 | REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: wetdiag !! Soil wetness on the thermodynamical levels |
---|
179 | !! (1, ngrnd) (0-1, dimensionless). corresponds to the |
---|
180 | !! relative soil humidity to the wilting point when |
---|
181 | !! the 11-layers hydrology (hydrol) is used, see more |
---|
182 | !! precisions in thermosoil_humlev. |
---|
183 | CONTAINS |
---|
184 | |
---|
185 | |
---|
186 | !! ================================================================================================================================ |
---|
187 | !! SUBROUTINE : thermosoil_main |
---|
188 | !! |
---|
189 | !>\BRIEF Thermosoil_main computes the soil thermal properties and dynamics, ie solves |
---|
190 | !! the heat diffusion equation within the soil. The soil temperature profile is |
---|
191 | !! then interpolated onto the diagnostic axis. |
---|
192 | !! |
---|
193 | !! DESCRIPTION : The resolution of the soil heat diffusion equation |
---|
194 | !! relies on a numerical finite-difference implicit scheme |
---|
195 | !! fully described in the reference and in the header of the thermosoil module. |
---|
196 | !! - The dependency of the previous timestep hidden in the |
---|
197 | !! integration coefficients cgrnd and dgrnd (EQ1), calculated in thermosoil_coef, and |
---|
198 | !! called at the end of the routine to prepare for the next timestep. |
---|
199 | !! - The effective computation of the new soil temperatures is performed in thermosoil_profile. |
---|
200 | !! |
---|
201 | !! - The calling sequence of thermosoil_main is summarized in the flowchart below. |
---|
202 | !! - Thermosoil_init and thermosoil_var_init initialize the variables from |
---|
203 | !! restart files or with default values; they also set up |
---|
204 | !! the vertical discretization for the numerical scheme. |
---|
205 | !! - thermosoil_coef calculates the coefficients for the numerical scheme for the very first iteration of thermosoil; |
---|
206 | !! after that, thermosoil_coef is called only at the end of the module to calculate the coefficients for the next timestep. |
---|
207 | !! - thermosoil_profile solves the numerical scheme.\n |
---|
208 | !! |
---|
209 | !! - Flags : one unique flag : THERMOSOIL_TPRO (to be set to the desired initial soil in-depth temperature in K; by default 280K) |
---|
210 | !! |
---|
211 | !! RECENT CHANGE(S) : None |
---|
212 | !! |
---|
213 | !! MAIN OUTPUT VARIABLE(S): vertically discretized soil temperatures ptn, soil |
---|
214 | !! thermal properties (pcapa, pkappa), apparent surface heat capacity (soilcap) |
---|
215 | !! and heat flux (soilflux) to be used in enerbil at the next timestep to solve |
---|
216 | !! the surface energy balance. |
---|
217 | !! |
---|
218 | !! REFERENCE(S) : |
---|
219 | !! - Hourdin, F. (1992). Study and numerical simulation of the general circulation of planetary atmospheres, |
---|
220 | !! Ph.D. thesis, Paris VII University. Remark: the part of F. Hourdin' s PhD thesis relative to the thermal |
---|
221 | !! integration scheme has been scanned and is provided along with the documentation, with name : |
---|
222 | !! Hourdin_1992_PhD_thermal_scheme.pdf |
---|
223 | !! |
---|
224 | !! FLOWCHART : |
---|
225 | !! \latexonly |
---|
226 | !! \includegraphics[scale = 1]{thermosoil_flowchart.png} |
---|
227 | !! \endlatexonly |
---|
228 | !! |
---|
229 | !! \n |
---|
230 | !_ ================================================================================================================================ |
---|
231 | |
---|
232 | SUBROUTINE thermosoil_main (kjit, kjpindex, dtradia, ldrestart_read, ldrestart_write, index, indexgrnd, & |
---|
233 | & temp_sol_new, snow, soilcap, soilflx, shumdiag, stempdiag, rest_id, hist_id, hist2_id) |
---|
234 | |
---|
235 | !! 0. Variable and parameter declaration |
---|
236 | |
---|
237 | !! 0.1 Input variables |
---|
238 | |
---|
239 | INTEGER(i_std), INTENT(in) :: kjit !! Time step number (unitless) |
---|
240 | INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (unitless) |
---|
241 | INTEGER(i_std),INTENT (in) :: rest_id,hist_id !! Restart_ file and history file identifier |
---|
242 | !! (unitless) |
---|
243 | INTEGER(i_std),INTENT (in) :: hist2_id !! history file 2 identifier (unitless) |
---|
244 | REAL(r_std), INTENT (in) :: dtradia !! model iteration time step in seconds (s) |
---|
245 | LOGICAL, INTENT(in) :: ldrestart_read !! Logical for restart files to be read |
---|
246 | !! (true/false) |
---|
247 | LOGICAL, INTENT(in) :: ldrestart_write !! Logical for restart files to be writen |
---|
248 | !! (true/false) |
---|
249 | INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: index !! Indeces of the points on the map (unitless) |
---|
250 | INTEGER(i_std),DIMENSION (kjpindex*ngrnd), INTENT (in):: indexgrnd !! Indeces of the points on the 3D map (vertical |
---|
251 | !! dimension towards the ground) (unitless) |
---|
252 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_sol_new !! Surface temperature at the present time-step, |
---|
253 | !! Ts @tex ($K$) @endtex |
---|
254 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: snow !! Snow mass @tex ($kg$) @endtex. |
---|
255 | !! Caveat: when there is snow on the |
---|
256 | !! ground, the snow is integrated into the soil for |
---|
257 | !! the calculation of the thermal dynamics. It means |
---|
258 | !! that the uppermost soil layers can completely or |
---|
259 | !! partially consist in snow. In the second case, zx1 |
---|
260 | !! and zx2 are the fraction of the soil layer |
---|
261 | !! consisting in snow and 'normal' soil, respectively |
---|
262 | !! This is calculated in thermosoil_coef. |
---|
263 | REAL(r_std),DIMENSION (kjpindex,nbdl), INTENT (in) :: shumdiag !! Relative soil humidity on the diagnostic axis |
---|
264 | !! (0-1, unitless). Caveats: when "hydrol" |
---|
265 | !! (the 11-layers hydrology) |
---|
266 | !! is used, this humidity is |
---|
267 | !! calculated with respect to the wilting point: |
---|
268 | !! shumdiag= (mc-mcw)/(mcs-mcw), with mc : moisture |
---|
269 | !! content; mcs : saturated soil moisture content; |
---|
270 | !! mcw: soil moisture content at the wilting point. |
---|
271 | !! When the 2-layers hydrology "hydrolc" is used, |
---|
272 | !! shumdiag is just a soil wetness index, from 0 to 1 |
---|
273 | !! but cannot direcly be linked to a soil moisture |
---|
274 | !! content. |
---|
275 | |
---|
276 | !! 0.2 Output variables |
---|
277 | |
---|
278 | REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: soilcap !! apparent surface heat capacity |
---|
279 | !! @tex ($J m^{-2} K^{-1}$) @endtex |
---|
280 | REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: soilflx !! apparent soil heat flux @tex ($W m^{-2}$) @endtex |
---|
281 | !! , positive |
---|
282 | !! towards the soil, writen as Qg (ground heat flux) |
---|
283 | !! in the history files, and computed at the end of |
---|
284 | !! thermosoil for the calculation of Ts in enerbil, |
---|
285 | !! see EQ3. |
---|
286 | REAL(r_std),DIMENSION (kjpindex,nbdl), INTENT (inout) :: stempdiag !! diagnostic temperature profile @tex ($K$) @endtex |
---|
287 | !! , eg on the |
---|
288 | !! diagnostic axis (levels:1:nbdl). The soil |
---|
289 | !! temperature is put on this diagnostic axis to be |
---|
290 | !! used by other modules (slowproc.f90; routing.f90; |
---|
291 | !! hydrol or hydrolc when a frozen soil |
---|
292 | !! parametrization is used..) |
---|
293 | |
---|
294 | !! 0.3 Modified variables |
---|
295 | |
---|
296 | !! 0.4 Local variables |
---|
297 | |
---|
298 | REAL(r_std),DIMENSION (kjpindex,ngrnd) :: temp !! buffer |
---|
299 | REAL(r_std),DIMENSION (kjpindex,ngrnd-1) :: temp1 !! buffer |
---|
300 | REAL(r_std),DIMENSION (kjpindex) :: temp2 !! buffer |
---|
301 | !_ ================================================================================================================================ |
---|
302 | |
---|
303 | !! 1. do initialisation |
---|
304 | |
---|
305 | IF (l_first_thermosoil) THEN |
---|
306 | |
---|
307 | IF (long_print) WRITE (numout,*) ' l_first_thermosoil : call thermosoil_init ' |
---|
308 | |
---|
309 | |
---|
310 | !! 1.1. Allocate and initialize soil temperatures variables |
---|
311 | !! by reading restart files or using default values. |
---|
312 | CALL thermosoil_init (kjit, ldrestart_read, kjpindex, index, rest_id) |
---|
313 | |
---|
314 | |
---|
315 | !! 1.2.Computes physical constants and arrays; initializes soil thermal properties; produces the first stempdiag |
---|
316 | !! Computes some physical constants and arrays depending on the soil vertical discretization |
---|
317 | !! (lskin, cstgrnd, zz, zz_coef, dz1, dz2); get the vertical humidity onto the thermal levels, and |
---|
318 | !! initializes soil thermal properties (pkappa, pcapa); produces the first temperature diagnostic stempdiag. |
---|
319 | CALL thermosoil_var_init (kjpindex, zz, zz_coef, dz1, dz2, pkappa, pcapa, pcapa_en, shumdiag, stempdiag) |
---|
320 | |
---|
321 | |
---|
322 | !! 1.3. Computes cgrd, dgrd, soilflx and soilcap coefficients from restart values or initialisation values. |
---|
323 | CALL thermosoil_coef (kjpindex, dtradia, temp_sol_new, snow, ptn, soilcap, soilflx, zz, dz1, dz2, z1, zdz1,& |
---|
324 | & zdz2, cgrnd, dgrnd, pcapa, pcapa_en, pkappa) |
---|
325 | |
---|
326 | !! 1.4. call to thermosoil_energy, if you wish to perform some checks (?) |
---|
327 | !!?? the usefulness of this routine seems questionable. |
---|
328 | CALL thermosoil_energy (kjpindex, temp_sol_new, soilcap, .TRUE.) |
---|
329 | |
---|
330 | !! 1.5. read restart files for other variables than ptn. |
---|
331 | !!?? mind the use of ok_var here. |
---|
332 | !!?? ok_var is a function of sechiba_io_p.f90, documented as follows : |
---|
333 | !!!! pour déclancher les restarts rajoutés avec un paramÚtre externe |
---|
334 | !!FUNCTION ok_var ( varname ) |
---|
335 | !!CHARACTER(LEN=*), INTENT(IN) :: varname |
---|
336 | !!LOGICAL ok_var |
---|
337 | !!ok_var=.FALSE. |
---|
338 | !!CALL getin_p(varname, ok_var) |
---|
339 | !!END FUNCTION ok_var |
---|
340 | !! |
---|
341 | !! from what we understand, it looks for the chain varname in |
---|
342 | !!run.def; if absent, returns .FALSE., and the variable named |
---|
343 | !!'varname' is not searched for in the restart. This looks like a |
---|
344 | !!trick to read variables in restart files when they are not read |
---|
345 | !!there by default. For all variables in the following sequence, ok_var |
---|
346 | !!is by default false, so don' t bother about this. |
---|
347 | !! this is also logical as those variables have been initialized |
---|
348 | !!above. |
---|
349 | !!?? so maybe this part of the code could be deleted to add clarity. |
---|
350 | IF (ldrestart_read) THEN |
---|
351 | IF (long_print) WRITE (numout,*) ' we have to READ a restart file for THERMOSOIL variables' |
---|
352 | |
---|
353 | var_name= 'cgrnd' |
---|
354 | CALL ioconf_setatt('UNITS', '-') |
---|
355 | CALL ioconf_setatt('LONG_NAME','Cgrnd coefficient.') |
---|
356 | IF ( ok_var(var_name) ) THEN |
---|
357 | CALL restget_p (rest_id, var_name, nbp_glo, ngrnd-1, 1, kjit, .TRUE., temp1, "gather", nbp_glo, index_g) |
---|
358 | IF (MINVAL(temp1) < MAXVAL(temp1) .OR. MAXVAL(temp1) .NE. val_exp) THEN |
---|
359 | cgrnd(:,:)=temp1(:,:) |
---|
360 | ENDIF |
---|
361 | ENDIF |
---|
362 | |
---|
363 | var_name= 'dgrnd' |
---|
364 | CALL ioconf_setatt('UNITS', '-') |
---|
365 | CALL ioconf_setatt('LONG_NAME','Dgrnd coefficient.') |
---|
366 | IF ( ok_var(var_name) ) THEN |
---|
367 | CALL restget_p (rest_id, var_name, nbp_glo, ngrnd-1, 1, kjit, .TRUE., temp1, "gather", nbp_glo, index_g) |
---|
368 | IF (MINVAL(temp1) < MAXVAL(temp1) .OR. MAXVAL(temp1) .NE. val_exp) THEN |
---|
369 | dgrnd(:,:)=temp1(:,:) |
---|
370 | ENDIF |
---|
371 | ENDIF |
---|
372 | |
---|
373 | var_name= 'z1' |
---|
374 | CALL ioconf_setatt('UNITS', '-') |
---|
375 | CALL ioconf_setatt('LONG_NAME','?.') |
---|
376 | IF ( ok_var(var_name) ) THEN |
---|
377 | CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., temp2, "gather", nbp_glo, index_g) |
---|
378 | IF (MINVAL(temp2) < MAXVAL(temp2) .OR. MAXVAL(temp2) .NE. val_exp) THEN |
---|
379 | z1(:)=temp2(:) |
---|
380 | ENDIF |
---|
381 | ENDIF |
---|
382 | |
---|
383 | var_name= 'pcapa' |
---|
384 | CALL ioconf_setatt('UNITS', '-') |
---|
385 | CALL ioconf_setatt('LONG_NAME','?.') |
---|
386 | IF ( ok_var(var_name) ) THEN |
---|
387 | CALL restget_p (rest_id, var_name, nbp_glo, ngrnd, 1, kjit, .TRUE., temp, "gather", nbp_glo, index_g) |
---|
388 | IF (MINVAL(temp) < MAXVAL(temp) .OR. MAXVAL(temp) .NE. val_exp) THEN |
---|
389 | pcapa(:,:)=temp(:,:) |
---|
390 | ENDIF |
---|
391 | ENDIF |
---|
392 | |
---|
393 | var_name= 'pcapa_en' |
---|
394 | CALL ioconf_setatt('UNITS', '-') |
---|
395 | CALL ioconf_setatt('LONG_NAME','?.') |
---|
396 | IF ( ok_var(var_name) ) THEN |
---|
397 | CALL restget_p (rest_id, var_name, nbp_glo, ngrnd, 1, kjit, .TRUE., temp, "gather", nbp_glo, index_g) |
---|
398 | IF (MINVAL(temp) < MAXVAL(temp) .OR. MAXVAL(temp) .NE. val_exp) THEN |
---|
399 | pcapa_en(:,:)=temp(:,:) |
---|
400 | ENDIF |
---|
401 | ENDIF |
---|
402 | |
---|
403 | var_name= 'pkappa' |
---|
404 | CALL ioconf_setatt('UNITS', '-') |
---|
405 | CALL ioconf_setatt('LONG_NAME','?.') |
---|
406 | IF ( ok_var(var_name) ) THEN |
---|
407 | CALL restget_p (rest_id, var_name, nbp_glo, ngrnd, 1, kjit, .TRUE., temp, "gather", nbp_glo, index_g) |
---|
408 | IF (MINVAL(temp) < MAXVAL(temp) .OR. MAXVAL(temp) .NE. val_exp) THEN |
---|
409 | pkappa(:,:)=temp(:,:) |
---|
410 | ENDIF |
---|
411 | ENDIF |
---|
412 | |
---|
413 | var_name= 'zdz1' |
---|
414 | CALL ioconf_setatt('UNITS', '-') |
---|
415 | CALL ioconf_setatt('LONG_NAME','?.') |
---|
416 | IF ( ok_var(var_name) ) THEN |
---|
417 | CALL restget_p (rest_id, var_name, nbp_glo, ngrnd-1, 1, kjit, .TRUE., temp1, "gather", nbp_glo, index_g) |
---|
418 | IF (MINVAL(temp1) < MAXVAL(temp1) .OR. MAXVAL(temp1) .NE. val_exp) THEN |
---|
419 | zdz1(:,:)=temp1(:,:) |
---|
420 | ENDIF |
---|
421 | ENDIF |
---|
422 | |
---|
423 | var_name= 'zdz2' |
---|
424 | CALL ioconf_setatt('UNITS', '-') |
---|
425 | CALL ioconf_setatt('LONG_NAME','?.') |
---|
426 | IF ( ok_var(var_name) ) THEN |
---|
427 | CALL restget_p (rest_id, var_name, nbp_glo, ngrnd, 1, kjit, .TRUE., temp, "gather", nbp_glo, index_g) |
---|
428 | IF (MINVAL(temp) < MAXVAL(temp) .OR. MAXVAL(temp) .NE. val_exp) THEN |
---|
429 | zdz2(:,:)=temp(:,:) |
---|
430 | ENDIF |
---|
431 | ENDIF |
---|
432 | |
---|
433 | var_name='temp_sol_beg' |
---|
434 | CALL ioconf_setatt('UNITS', 'K') |
---|
435 | CALL ioconf_setatt('LONG_NAME','Old Surface temperature') |
---|
436 | IF ( ok_var(var_name) ) THEN |
---|
437 | CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., temp2, "gather", nbp_glo, index_g) |
---|
438 | IF (MINVAL(temp2) < MAXVAL(temp2) .OR. MAXVAL(temp2) .NE. val_exp) THEN |
---|
439 | temp_sol_beg(:) = temp2(:) |
---|
440 | ENDIF |
---|
441 | ENDIF |
---|
442 | |
---|
443 | ENDIF !ldrestart_read |
---|
444 | |
---|
445 | RETURN |
---|
446 | |
---|
447 | ENDIF !l_first_thermosoil |
---|
448 | |
---|
449 | |
---|
450 | !! 2. Prepares the restart files for the next simulation |
---|
451 | |
---|
452 | !!?? do all the coefficients (cgrnd, dgrnd...) be put in the restart file |
---|
453 | !! as they are by default not read there, but calculated in |
---|
454 | !!thermosoil_var_init from the restart or initial temperature ? |
---|
455 | !! exceptions are soilcap and soilflx, used in enerbil, and of course ptn. |
---|
456 | IF (ldrestart_write) THEN |
---|
457 | |
---|
458 | IF (long_print) WRITE (numout,*) ' we have to complete restart file with THERMOSOIL variables' |
---|
459 | |
---|
460 | var_name= 'ptn' |
---|
461 | CALL restput_p(rest_id, var_name, nbp_glo, ngrnd, 1, kjit, ptn, 'scatter', nbp_glo, index_g) |
---|
462 | |
---|
463 | var_name= 'cgrnd' |
---|
464 | CALL restput_p(rest_id, var_name, nbp_glo, ngrnd-1, 1, kjit, cgrnd, 'scatter', nbp_glo, index_g) |
---|
465 | var_name= 'dgrnd' |
---|
466 | CALL restput_p(rest_id, var_name, nbp_glo, ngrnd-1, 1, kjit, dgrnd, 'scatter', nbp_glo, index_g) |
---|
467 | |
---|
468 | var_name= 'z1' |
---|
469 | CALL restput_p(rest_id, var_name, nbp_glo, 1, 1, kjit, z1, 'scatter', nbp_glo, index_g) |
---|
470 | |
---|
471 | var_name= 'pcapa' |
---|
472 | CALL restput_p(rest_id, var_name, nbp_glo, ngrnd, 1, kjit, pcapa, 'scatter', nbp_glo, index_g) |
---|
473 | |
---|
474 | var_name= 'pcapa_en' |
---|
475 | CALL restput_p(rest_id, var_name, nbp_glo, ngrnd, 1, kjit, pcapa_en, 'scatter', nbp_glo, index_g) |
---|
476 | |
---|
477 | var_name= 'pkappa' |
---|
478 | CALL restput_p(rest_id, var_name, nbp_glo, ngrnd, 1, kjit, pkappa, 'scatter', nbp_glo, index_g) |
---|
479 | |
---|
480 | var_name= 'zdz1' |
---|
481 | CALL restput_p(rest_id, var_name, nbp_glo, ngrnd-1, 1, kjit, zdz1, 'scatter', nbp_glo, index_g) |
---|
482 | |
---|
483 | var_name= 'zdz2' |
---|
484 | CALL restput_p(rest_id, var_name, nbp_glo, ngrnd, 1, kjit, zdz2, 'scatter', nbp_glo, index_g) |
---|
485 | |
---|
486 | var_name= 'temp_sol_beg' |
---|
487 | CALL restput_p(rest_id, var_name, nbp_glo, 1, 1, kjit, temp_sol_beg, 'scatter', nbp_glo, index_g) |
---|
488 | |
---|
489 | var_name= 'soilcap' |
---|
490 | CALL restput_p(rest_id, var_name, nbp_glo, 1, 1, kjit, soilcap, 'scatter', nbp_glo, index_g) |
---|
491 | |
---|
492 | var_name= 'soilflx' |
---|
493 | CALL restput_p(rest_id, var_name, nbp_glo, 1, 1, kjit, soilflx, 'scatter', nbp_glo, index_g) |
---|
494 | |
---|
495 | ! read in enerbil |
---|
496 | var_name= 'temp_sol_new' |
---|
497 | CALL restput_p(rest_id, var_name, nbp_glo, 1, 1, kjit, temp_sol_new, 'scatter', nbp_glo, index_g) |
---|
498 | |
---|
499 | RETURN |
---|
500 | |
---|
501 | END IF !ldrestart_write |
---|
502 | |
---|
503 | !! 3. Put the soil wetness diagnostic on the levels of the soil temperature |
---|
504 | |
---|
505 | !!?? this could logically be put just before the last call to |
---|
506 | !!thermosoil_coef, as the results are used there... |
---|
507 | CALL thermosoil_humlev(kjpindex, shumdiag) |
---|
508 | |
---|
509 | |
---|
510 | !! 4. Effective computation of the soil temperatures profile, using the cgrd and dgrd coefficients from previsou tstep. |
---|
511 | |
---|
512 | CALL thermosoil_profile (kjpindex, temp_sol_new, ptn, stempdiag) |
---|
513 | |
---|
514 | !! 5. Call to thermosoil_energy, still to be clarified.. |
---|
515 | |
---|
516 | CALL thermosoil_energy (kjpindex, temp_sol_new, soilcap, .FALSE.) |
---|
517 | |
---|
518 | !! 6. Writing the history files according to the ALMA standards (or not..) |
---|
519 | |
---|
520 | !in only one file (hist2_id <=0) or in 2 different files (hist2_id >0). |
---|
521 | IF ( .NOT. almaoutput ) THEN |
---|
522 | CALL histwrite(hist_id, 'ptn', kjit, ptn, kjpindex*ngrnd, indexgrnd) |
---|
523 | ELSE |
---|
524 | CALL histwrite(hist_id, 'SoilTemp', kjit, ptn, kjpindex*ngrnd, indexgrnd) |
---|
525 | CALL histwrite(hist_id, 'Qg', kjit, soilflx, kjpindex, index) |
---|
526 | CALL histwrite(hist_id, 'DelSurfHeat', kjit, surfheat_incr, kjpindex, index) |
---|
527 | CALL histwrite(hist_id, 'DelColdCont', kjit, coldcont_incr, kjpindex, index) |
---|
528 | ENDIF |
---|
529 | IF ( hist2_id > 0 ) THEN |
---|
530 | IF ( .NOT. almaoutput ) THEN |
---|
531 | CALL histwrite(hist2_id, 'ptn', kjit, ptn, kjpindex*ngrnd, indexgrnd) |
---|
532 | ELSE |
---|
533 | CALL histwrite(hist2_id, 'SoilTemp', kjit, ptn, kjpindex*ngrnd, indexgrnd) |
---|
534 | CALL histwrite(hist2_id, 'Qg', kjit, soilflx, kjpindex, index) |
---|
535 | CALL histwrite(hist2_id, 'DelSurfHeat', kjit, surfheat_incr, kjpindex, index) |
---|
536 | CALL histwrite(hist2_id, 'DelColdCont', kjit, coldcont_incr, kjpindex, index) |
---|
537 | ENDIF |
---|
538 | ENDIF |
---|
539 | |
---|
540 | !! 7. A last final call to thermosoil_coef |
---|
541 | |
---|
542 | !! A last final call to thermosoil_coef, which calculates the different |
---|
543 | !!coefficients (cgrnd, dgrnd, dz1, z1, zdz2, soilcap, soilflx) from this time step to be |
---|
544 | !!used at the next time step, either in the surface temperature calculation |
---|
545 | !!(soilcap, soilflx) or in the soil thermal numerical scheme. |
---|
546 | CALL thermosoil_coef (kjpindex, dtradia, temp_sol_new, snow, ptn, soilcap, soilflx, zz, dz1, dz2, z1, zdz1,& |
---|
547 | & zdz2, cgrnd, dgrnd, pcapa, pcapa_en, pkappa) |
---|
548 | |
---|
549 | IF (long_print) WRITE (numout,*) ' thermosoil_main done ' |
---|
550 | |
---|
551 | END SUBROUTINE thermosoil_main |
---|
552 | |
---|
553 | |
---|
554 | !! ================================================================================================================================ |
---|
555 | !! SUBROUTINE : thermosoil_init |
---|
556 | !! |
---|
557 | !>\BRIEF Allocates local and global arrays; initializes soil temperatures using either restart files |
---|
558 | !! or a fixed value set by the flag THERMOSOIL_TPRO. |
---|
559 | !! |
---|
560 | !! DESCRIPTION : flag : THERMOSOIL_TPRO (to be set to the desired initial temperature in K; by default 280K). |
---|
561 | !! |
---|
562 | !! RECENT CHANGE(S) : None |
---|
563 | !! |
---|
564 | !! MAIN OUTPUT VARIABLE(S): None |
---|
565 | !! |
---|
566 | !! REFERENCE(S) : None |
---|
567 | !! |
---|
568 | !! FLOWCHART : None |
---|
569 | !! \n |
---|
570 | !_ ================================================================================================================================ |
---|
571 | |
---|
572 | SUBROUTINE thermosoil_init(kjit, ldrestart_read, kjpindex, index, rest_id) |
---|
573 | |
---|
574 | !! 0. Variable and parameter declaration |
---|
575 | |
---|
576 | !! 0.1 Input variables |
---|
577 | |
---|
578 | INTEGER(i_std), INTENT (in) :: kjit !! Time step number (unitless) |
---|
579 | LOGICAL,INTENT (in) :: ldrestart_read !! Logical for restart file to read (true/false) |
---|
580 | INTEGER(i_std), INTENT (in) :: kjpindex !! Domain size (unitless) |
---|
581 | INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: index !! Indeces of the points on the map (unitless) |
---|
582 | INTEGER(i_std), INTENT (in) :: rest_id !! Restart file identifier (unitless) |
---|
583 | |
---|
584 | !! 0.2 Output variables |
---|
585 | |
---|
586 | !! 0.3 Modified variables |
---|
587 | |
---|
588 | !! 0.4 Local variables |
---|
589 | |
---|
590 | INTEGER(i_std) :: ier |
---|
591 | !_ ================================================================================================================================ |
---|
592 | |
---|
593 | !! 1. Initialisation |
---|
594 | |
---|
595 | !! Initialisation has to be done only one time, so the logical |
---|
596 | !! logical l_first_thermosoil has to be set to .FALSE. now.. |
---|
597 | IF (l_first_thermosoil) THEN |
---|
598 | l_first_thermosoil=.FALSE. |
---|
599 | ELSE |
---|
600 | WRITE (numout,*) ' l_first_thermosoil false . we stop ' |
---|
601 | STOP 'thermosoil_init' |
---|
602 | ENDIF |
---|
603 | |
---|
604 | !! 2. Arrays allocations |
---|
605 | |
---|
606 | ALLOCATE (ptn(kjpindex,ngrnd),stat=ier) |
---|
607 | IF (ier.NE.0) THEN |
---|
608 | WRITE (numout,*) ' error in ptn allocation. We stop. We need ',kjpindex,' fois ',ngrnd,' words = '& |
---|
609 | & , kjpindex*ngrnd |
---|
610 | STOP 'thermosoil_init' |
---|
611 | END IF |
---|
612 | |
---|
613 | ALLOCATE (z1(kjpindex),stat=ier) |
---|
614 | IF (ier.NE.0) THEN |
---|
615 | WRITE (numout,*) ' error in z1 allocation. We STOP. We need ',kjpindex,' words ' |
---|
616 | STOP 'thermosoil_init' |
---|
617 | END IF |
---|
618 | |
---|
619 | ALLOCATE (cgrnd(kjpindex,ngrnd-1),stat=ier) |
---|
620 | IF (ier.NE.0) THEN |
---|
621 | WRITE (numout,*) ' error in cgrnd allocation. We STOP. We need ',kjpindex,' fois ',ngrnd-1 ,' words = '& |
---|
622 | & , kjpindex*(ngrnd-1) |
---|
623 | STOP 'thermosoil_init' |
---|
624 | END IF |
---|
625 | |
---|
626 | ALLOCATE (dgrnd(kjpindex,ngrnd-1),stat=ier) |
---|
627 | IF (ier.NE.0) THEN |
---|
628 | WRITE (numout,*) ' error in dgrnd allocation. We STOP. We need ',kjpindex,' fois ',ngrnd-1 ,' words = '& |
---|
629 | & , kjpindex*(ngrnd-1) |
---|
630 | STOP 'thermosoil_init' |
---|
631 | END IF |
---|
632 | |
---|
633 | ALLOCATE (pcapa(kjpindex,ngrnd),stat=ier) |
---|
634 | IF (ier.NE.0) THEN |
---|
635 | WRITE (numout,*) ' error in pcapa allocation. We STOP. We need ',kjpindex,' fois ',ngrnd ,' words = '& |
---|
636 | & , kjpindex*ngrnd |
---|
637 | STOP 'thermosoil_init' |
---|
638 | END IF |
---|
639 | |
---|
640 | ALLOCATE (pkappa(kjpindex,ngrnd),stat=ier) |
---|
641 | IF (ier.NE.0) THEN |
---|
642 | WRITE (numout,*) ' error in pkappa allocation. We STOP. We need ',kjpindex,' fois ',ngrnd ,' words = '& |
---|
643 | & , kjpindex*ngrnd |
---|
644 | STOP 'thermosoil_init' |
---|
645 | END IF |
---|
646 | |
---|
647 | ALLOCATE (zdz1(kjpindex,ngrnd-1),stat=ier) |
---|
648 | IF (ier.NE.0) THEN |
---|
649 | WRITE (numout,*) ' error in zdz1 allocation. We STOP. We need ',kjpindex,' fois ',ngrnd-1 ,' words = '& |
---|
650 | & , kjpindex*(ngrnd-1) |
---|
651 | STOP 'thermosoil_init' |
---|
652 | END IF |
---|
653 | |
---|
654 | ALLOCATE (zdz2(kjpindex,ngrnd),stat=ier) |
---|
655 | IF (ier.NE.0) THEN |
---|
656 | WRITE (numout,*) ' error in zdz2 allocation. We STOP. We need ',kjpindex,' fois ',ngrnd ,' words = '& |
---|
657 | & , kjpindex*ngrnd |
---|
658 | STOP 'thermosoil_init' |
---|
659 | END IF |
---|
660 | |
---|
661 | ALLOCATE (surfheat_incr(kjpindex),stat=ier) |
---|
662 | IF (ier.NE.0) THEN |
---|
663 | WRITE (numout,*) ' error in surfheat_incr allocation. We STOP. We need ',kjpindex,' words = '& |
---|
664 | & , kjpindex |
---|
665 | STOP 'thermosoil_init' |
---|
666 | END IF |
---|
667 | |
---|
668 | ALLOCATE (coldcont_incr(kjpindex),stat=ier) |
---|
669 | IF (ier.NE.0) THEN |
---|
670 | WRITE (numout,*) ' error in coldcont_incr allocation. We STOP. We need ',kjpindex,' words = '& |
---|
671 | & , kjpindex |
---|
672 | STOP 'thermosoil_init' |
---|
673 | END IF |
---|
674 | |
---|
675 | ALLOCATE (pcapa_en(kjpindex,ngrnd),stat=ier) |
---|
676 | IF (ier.NE.0) THEN |
---|
677 | WRITE (numout,*) ' error in pcapa_en allocation. We STOP. We need ',kjpindex,' fois ',ngrnd ,' words = '& |
---|
678 | & , kjpindex*ngrnd |
---|
679 | STOP 'thermosoil_init' |
---|
680 | END IF |
---|
681 | |
---|
682 | ALLOCATE (ptn_beg(kjpindex,ngrnd),stat=ier) |
---|
683 | IF (ier.NE.0) THEN |
---|
684 | WRITE (numout,*) ' error in ptn_beg allocation. We STOP. We need ',kjpindex,' fois ',ngrnd ,' words = '& |
---|
685 | & , kjpindex*ngrnd |
---|
686 | STOP 'thermosoil_init' |
---|
687 | END IF |
---|
688 | |
---|
689 | ALLOCATE (temp_sol_beg(kjpindex),stat=ier) |
---|
690 | IF (ier.NE.0) THEN |
---|
691 | WRITE (numout,*) ' error in temp_sol_beg allocation. We STOP. We need ',kjpindex,' words = '& |
---|
692 | & , kjpindex |
---|
693 | STOP 'thermosoil_init' |
---|
694 | END IF |
---|
695 | |
---|
696 | ALLOCATE (wetdiag(kjpindex,ngrnd),stat=ier) |
---|
697 | IF (ier.NE.0) THEN |
---|
698 | WRITE (numout,*) ' error in wetdiag allocation. We STOP. We need ',kjpindex,' fois ',ngrnd ,' words = '& |
---|
699 | & , kjpindex*ngrnd |
---|
700 | STOP 'thermosoil_init' |
---|
701 | END IF |
---|
702 | |
---|
703 | !! 3. Reads restart files for soil temperatures only |
---|
704 | |
---|
705 | !! Reads restart files for soil temperatures only. If no restart file is |
---|
706 | !! found, the initial soil temperature is by default set to 280K at all depths. The user |
---|
707 | !! can decide to initialize soil temperatures at an other value, in which case he should set the flag THERMOSOIL_TPRO |
---|
708 | !! to this specific value in the run.def. |
---|
709 | IF (ldrestart_read) THEN |
---|
710 | IF (long_print) WRITE (numout,*) ' we have to READ a restart file for THERMOSOIL variables' |
---|
711 | |
---|
712 | var_name= 'ptn' |
---|
713 | CALL ioconf_setatt('UNITS', 'K') |
---|
714 | CALL ioconf_setatt('LONG_NAME','Soil Temperature profile') |
---|
715 | CALL restget_p (rest_id, var_name, nbp_glo, ngrnd, 1, kjit, .TRUE., ptn, "gather", nbp_glo, index_g) |
---|
716 | ! |
---|
717 | !Config Key = THERMOSOIL_TPRO |
---|
718 | !Config Desc = Initial soil temperature profile if not found in restart |
---|
719 | !Config Def = 280. |
---|
720 | !Config Help = The initial value of the temperature profile in the soil if |
---|
721 | !Config its value is not found in the restart file. This should only |
---|
722 | !Config be used if the model is started without a restart file. Here |
---|
723 | !Config we only require one value as we will assume a constant |
---|
724 | !Config throughout the column. |
---|
725 | ! |
---|
726 | CALL setvar_p (ptn, val_exp,'THERMOSOIL_TPRO',280._r_std) |
---|
727 | |
---|
728 | ENDIF |
---|
729 | |
---|
730 | IF (long_print) WRITE (numout,*) ' thermosoil_init done ' |
---|
731 | |
---|
732 | END SUBROUTINE thermosoil_init |
---|
733 | |
---|
734 | |
---|
735 | !! ================================================================================================================================ |
---|
736 | !! SUBROUTINE : thermosoil_clear |
---|
737 | !! |
---|
738 | !>\BRIEF Sets the flag l_first_thermosoil to true and desallocates the allocated arrays. |
---|
739 | !! ??!! the call of thermosoil_clear originates from sechiba_clear but the calling sequence and |
---|
740 | !! its purpose require further investigation. |
---|
741 | !! |
---|
742 | !! DESCRIPTION : None |
---|
743 | !! |
---|
744 | !! RECENT CHANGE(S) : None |
---|
745 | !! |
---|
746 | !! MAIN OUTPUT VARIABLE(S): None |
---|
747 | !! |
---|
748 | !! REFERENCE(S) : None |
---|
749 | !! |
---|
750 | !! FLOWCHART : None |
---|
751 | !! \n |
---|
752 | !_ ================================================================================================================================ |
---|
753 | |
---|
754 | SUBROUTINE thermosoil_clear() |
---|
755 | |
---|
756 | l_first_thermosoil=.TRUE. |
---|
757 | |
---|
758 | IF ( ALLOCATED (ptn)) DEALLOCATE (ptn) |
---|
759 | IF ( ALLOCATED (z1)) DEALLOCATE (z1) |
---|
760 | IF ( ALLOCATED (cgrnd)) DEALLOCATE (cgrnd) |
---|
761 | IF ( ALLOCATED (dgrnd)) DEALLOCATE (dgrnd) |
---|
762 | IF ( ALLOCATED (pcapa)) DEALLOCATE (pcapa) |
---|
763 | IF ( ALLOCATED (pkappa)) DEALLOCATE (pkappa) |
---|
764 | IF ( ALLOCATED (zdz1)) DEALLOCATE (zdz1) |
---|
765 | IF ( ALLOCATED (zdz2)) DEALLOCATE (zdz2) |
---|
766 | IF ( ALLOCATED (pcapa_en)) DEALLOCATE (pcapa_en) |
---|
767 | IF ( ALLOCATED (ptn_beg)) DEALLOCATE (ptn_beg) |
---|
768 | IF ( ALLOCATED (temp_sol_beg)) DEALLOCATE (temp_sol_beg) |
---|
769 | IF ( ALLOCATED (surfheat_incr)) DEALLOCATE (surfheat_incr) |
---|
770 | IF ( ALLOCATED (coldcont_incr)) DEALLOCATE (coldcont_incr) |
---|
771 | IF ( ALLOCATED (wetdiag)) DEALLOCATE (wetdiag) |
---|
772 | |
---|
773 | END SUBROUTINE thermosoil_clear |
---|
774 | |
---|
775 | |
---|
776 | !! ================================================================================================================================ |
---|
777 | !! FUNCTION : fz |
---|
778 | !! |
---|
779 | !>\BRIEF fz(rk), the function's result, is the "rk"th element of a geometric series |
---|
780 | !! with first element fz1 and ration zalph. |
---|
781 | !! |
---|
782 | !! DESCRIPTION : This function is used to calculate the depths of the boudaries of the thermal layers (zz_coef) and |
---|
783 | !! of the numerical nodes (zz) of the thermal scheme. Formulae to get the adimensional depths are followings : |
---|
784 | !! zz(jg) = fz(REAL(jg,r_std) - undemi); \n |
---|
785 | !! zz_coef(jg) = fz(REAL(jg,r_std)) |
---|
786 | !! |
---|
787 | !! RECENT CHANGE(S) : None |
---|
788 | !! |
---|
789 | !! RETURN VALUE : fz(rk) |
---|
790 | !! |
---|
791 | !! REFERENCE(S) : None |
---|
792 | !! |
---|
793 | !! FLOWCHART : None |
---|
794 | !! \n |
---|
795 | !_ ================================================================================================================================ |
---|
796 | |
---|
797 | FUNCTION fz(rk) RESULT (fz_result) |
---|
798 | |
---|
799 | !! 0. Variables and parameter declaration |
---|
800 | |
---|
801 | !! 0.1 Input variables |
---|
802 | |
---|
803 | REAL(r_std), INTENT(in) :: rk |
---|
804 | |
---|
805 | !! 0.2 Output variables |
---|
806 | |
---|
807 | REAL(r_std) :: fz_result |
---|
808 | |
---|
809 | !! 0.3 Modified variables |
---|
810 | |
---|
811 | !! 0.4 Local variables |
---|
812 | |
---|
813 | !_ ================================================================================================================================ |
---|
814 | |
---|
815 | fz_result = fz1 * (zalph ** rk - un) / (zalph - un) |
---|
816 | |
---|
817 | END FUNCTION fz |
---|
818 | |
---|
819 | |
---|
820 | !! ================================================================================================================================ |
---|
821 | !! SUBROUTINE : thermosoil_var_init |
---|
822 | !! |
---|
823 | !>\BRIEF Define and initializes the soil thermal parameters |
---|
824 | !! |
---|
825 | !! DESCRIPTION : This routine\n |
---|
826 | !! 1. Defines the parameters ruling the vertical grid of the thermal scheme (fz1, zalpha).\n |
---|
827 | !! 2. Defines the scaling coefficients for adimensional depths (lskin, cstgrnd, see explanation in the |
---|
828 | !! variables description of thermosoil_main). \n |
---|
829 | !! 3. Calculates the vertical discretization of the soil (zz, zz_coef, dz2) and the constants used |
---|
830 | !! in the numerical scheme and which depend only on the discretization (dz1, lambda).\n |
---|
831 | !! 4. Initializes the soil thermal parameters (capacity, conductivity) based on initial soil moisture content.\n |
---|
832 | !! 5. Produces a first temperature diagnostic based on temperature initialization.\n |
---|
833 | !! |
---|
834 | !! The scheme comprizes ngrnd=7 layers by default. |
---|
835 | !! The layer' s boundaries depths (zz_coef) follow a geometric series of ratio zalph=2 and first term fz1.\n |
---|
836 | !! zz_coef(jg)=fz1.(1-zalph^jg)/(1-zalph) \n |
---|
837 | !! The layers' boudaries depths are first calculated 'adimensionally', ie with a |
---|
838 | !! discretization adapted to EQ5. This discretization is chosen for its ability at |
---|
839 | !! reproducing a thermal signal with periods ranging from days to centuries. (see |
---|
840 | !! Hourdin, 1992). Typically, fz1 is chosen as : fz1=0.3*cstgrnd (with cstgrnd the |
---|
841 | !! adimensional attenuation depth). \n |
---|
842 | !! The factor lskin/cstgrnd is then used to go from adimensional depths to |
---|
843 | !! depths in m.\n |
---|
844 | !! zz(real)=lskin/cstgrnd*zz(adimensional)\n |
---|
845 | !! Similarly, the depths of the numerical nodes are first calculated |
---|
846 | !! adimensionally, then the conversion factor is applied.\n |
---|
847 | !! the numerical nodes (zz) are not exactly the layers' centers : their depths are calculated as follows:\n |
---|
848 | !! zz(jg)=fz1.(1-zalph^(jg-1/2))/(1-zalph)\n |
---|
849 | !! The values of zz and zz_coef used in the default thermal discretization are in the following table. |
---|
850 | !! \latexonly |
---|
851 | !! \includegraphics{thermosoil_var_init1.jpg} |
---|
852 | !! \endlatexonly\n |
---|
853 | !! |
---|
854 | !! RECENT CHANGE(S) : None |
---|
855 | !! |
---|
856 | !! MAIN OUTPUT VARIABLE(S) : None |
---|
857 | !! |
---|
858 | !! REFERENCE(S) : |
---|
859 | !! - Hourdin, F. (1992). Study and numerical simulation of the general circulation of |
---|
860 | !! planetary atmospheres, Ph.D. thesis, Paris VII University. |
---|
861 | !! |
---|
862 | !! FLOWCHART : None |
---|
863 | !! \n |
---|
864 | !_ ================================================================================================================================ |
---|
865 | |
---|
866 | SUBROUTINE thermosoil_var_init(kjpindex, zz, zz_coef, dz1, dz2, pkappa, pcapa, pcapa_en, shumdiag, stempdiag) |
---|
867 | |
---|
868 | !! 0. Variables and parameter declaration |
---|
869 | |
---|
870 | !! 0.1 Input variables |
---|
871 | |
---|
872 | INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (unitless) |
---|
873 | REAL(r_std), DIMENSION (kjpindex,nbdl), INTENT (in) :: shumdiag !! Relative soil humidity on the diagnostic axis |
---|
874 | !! (unitless), [0,1]. (see description of the |
---|
875 | !! variables of thermosoil_main for more |
---|
876 | !! explanations) |
---|
877 | |
---|
878 | !! 0.2 Output variables |
---|
879 | |
---|
880 | REAL(r_std), DIMENSION (ngrnd), INTENT(out) :: zz !! depths of the layers'numerical nodes |
---|
881 | !! @tex ($m$)@endtex |
---|
882 | REAL(r_std), DIMENSION (ngrnd), INTENT(out) :: zz_coef !! depths of the layers'boundaries |
---|
883 | !! @tex ($m$)@endtex |
---|
884 | REAL(r_std), DIMENSION (ngrnd), INTENT(out) :: dz1 !! numerical constant depending on the vertical |
---|
885 | !! thermal grid only @tex ($m^{-1}$) @endtex. |
---|
886 | !! (see description |
---|
887 | !! of the variables of thermosoil_main for more |
---|
888 | !! explanations) |
---|
889 | REAL(r_std), DIMENSION (ngrnd), INTENT(out) :: dz2 !! thicknesses of the soil thermal layers |
---|
890 | !! @tex ($m$) @endtex |
---|
891 | REAL(r_std), DIMENSION (kjpindex,ngrnd), INTENT(out) :: pcapa !! volumetric vertically discretized soil heat |
---|
892 | !! capacity @tex ($J K^{-1} m^{-3}$) @endtex |
---|
893 | REAL(r_std), DIMENSION (kjpindex,ngrnd), INTENT(out) :: pcapa_en !! volumetric vertically discretized heat |
---|
894 | !! capacity used in thermosoil_energy |
---|
895 | !! @tex ($J K^{-1} m^{-3}$) @endtex ; |
---|
896 | !! usefulness still to be clarified. |
---|
897 | REAL(r_std), DIMENSION (kjpindex,ngrnd), INTENT(out) :: pkappa !! vertically discretized soil thermal |
---|
898 | !! conductivity @tex ($W m^{-1} K^{-1}$) @endtex |
---|
899 | REAL(r_std), DIMENSION (kjpindex,nbdl), INTENT (out) :: stempdiag !! Diagnostic temperature profile @tex ($K$) |
---|
900 | !! @endtex |
---|
901 | |
---|
902 | !! 0.3 Modified variables |
---|
903 | |
---|
904 | !! 0.4 Local variables |
---|
905 | |
---|
906 | INTEGER(i_std) :: ier, ji, jg |
---|
907 | REAL(r_std) :: sum |
---|
908 | !_ ================================================================================================================================ |
---|
909 | |
---|
910 | !! 1. Initialization of the parameters of the vertical discretization and of the attenuation depths |
---|
911 | |
---|
912 | cstgrnd=SQRT(one_day / pi) |
---|
913 | lskin = SQRT(so_cond / so_capa * one_day / pi) |
---|
914 | fz1 = 0.3_r_std * cstgrnd |
---|
915 | zalph = deux |
---|
916 | |
---|
917 | !! 2. Computing the depth of the thermal levels (numerical nodes) and the layers boundaries |
---|
918 | |
---|
919 | !! Computing the depth of the thermal levels (numerical nodes) and |
---|
920 | !! the layers boundariesusing the so-called |
---|
921 | !! adimentional variable z' = z/lskin*cstgrnd (with z in m) |
---|
922 | |
---|
923 | !! 2.1 adimensional thicknesses of the layers |
---|
924 | DO jg=1,ngrnd |
---|
925 | |
---|
926 | !!?? code simplification hopefully possible here with up-to-date compilers ! |
---|
927 | !!! This needs to be solved soon. Either we allow CPP options in SECHIBA or the VPP |
---|
928 | !!! fixes its compiler |
---|
929 | !!!#ifdef VPP5000 |
---|
930 | dz2(jg) = fz(REAL(jg,r_std)-undemi+undemi) - fz(REAL(jg-1,r_std)-undemi+undemi) |
---|
931 | !!!#else |
---|
932 | !!! dz2(jg) = fz(REAL(jg,r_std)) - fz(REAL(jg-1,r_std)) |
---|
933 | !!!#endif |
---|
934 | ENDDO |
---|
935 | |
---|
936 | !! 2.2 adimentional depth of the numerical nodes and layers' boudaries |
---|
937 | DO jg=1,ngrnd |
---|
938 | zz(jg) = fz(REAL(jg,r_std) - undemi) |
---|
939 | zz_coef(jg) = fz(REAL(jg,r_std)-undemi+undemi) |
---|
940 | ENDDO |
---|
941 | |
---|
942 | !! 2.3 Converting to meters |
---|
943 | DO jg=1,ngrnd |
---|
944 | zz(jg) = zz(jg) / cstgrnd * lskin |
---|
945 | zz_coef(jg) = zz_coef(jg) / cstgrnd * lskin |
---|
946 | dz2(jg) = dz2(jg) / cstgrnd * lskin |
---|
947 | ENDDO |
---|
948 | |
---|
949 | !! 2.4 Computing some usefull constants for the numerical scheme |
---|
950 | DO jg=1,ngrnd-1 |
---|
951 | dz1(jg) = un / (zz(jg+1) - zz(jg)) |
---|
952 | ENDDO |
---|
953 | lambda = zz(1) * dz1(1) |
---|
954 | |
---|
955 | !! 2.5 Get the wetness profile on the thermal vertical grid from the diagnostic axis |
---|
956 | CALL thermosoil_humlev(kjpindex, shumdiag) |
---|
957 | |
---|
958 | !! 2.6 Thermal conductivity at all levels |
---|
959 | DO jg = 1,ngrnd |
---|
960 | DO ji = 1,kjpindex |
---|
961 | pkappa(ji,jg) = so_cond_dry + wetdiag(ji,jg)*(so_cond_wet - so_cond_dry) |
---|
962 | pcapa(ji,jg) = so_capa_dry + wetdiag(ji,jg)*(so_capa_wet - so_capa_dry) |
---|
963 | pcapa_en(ji,jg) = so_capa_dry + wetdiag(ji,jg)*(so_capa_wet - so_capa_dry) |
---|
964 | ENDDO |
---|
965 | ENDDO |
---|
966 | |
---|
967 | !! 3. Diagnostics : consistency checks on the vertical grid. |
---|
968 | |
---|
969 | WRITE (numout,*) 'diagnostic des niveaux dans le sol' !!?? to be changed, |
---|
970 | WRITE (numout,*) 'niveaux intermediaires et pleins' |
---|
971 | sum = zero |
---|
972 | DO jg=1,ngrnd |
---|
973 | sum = sum + dz2(jg) |
---|
974 | WRITE (numout,*) zz(jg),sum |
---|
975 | ENDDO |
---|
976 | |
---|
977 | |
---|
978 | !! 4. Compute a first diagnostic temperature profile |
---|
979 | |
---|
980 | CALL thermosoil_diaglev(kjpindex, stempdiag) |
---|
981 | |
---|
982 | IF (long_print) WRITE (numout,*) ' thermosoil_var_init done ' |
---|
983 | |
---|
984 | END SUBROUTINE thermosoil_var_init |
---|
985 | |
---|
986 | |
---|
987 | !! ================================================================================================================================ |
---|
988 | !! SUBROUTINE : thermosoil_coef |
---|
989 | !! |
---|
990 | !>\BRIEF Calculate soil thermal properties, integration coefficients, apparent heat flux, |
---|
991 | !! surface heat capacity, |
---|
992 | !! |
---|
993 | !! DESCRIPTION : This routine computes : \n |
---|
994 | !! 1. the soil thermal properties. \n |
---|
995 | !! 2. the integration coefficients of the thermal numerical scheme, cgrnd and dgrnd, |
---|
996 | !! which depend on the vertical grid and on soil properties, and are used at the next |
---|
997 | !! timestep.\n |
---|
998 | !! 3. the soil apparent heat flux and surface heat capacity soilflux |
---|
999 | !! and soilcap, used by enerbil to compute the surface temperature at the next |
---|
1000 | !! timestep.\n |
---|
1001 | !! - The soil thermal properties depend on water content (wetdiag) and on the presence |
---|
1002 | !! of snow : snow is integrated into the soil for the thermal calculations, ie if there |
---|
1003 | !! is snow on the ground, the first thermal layer(s) consist in snow, depending on the |
---|
1004 | !! snow-depth. If a layer consists out of snow and soil, wheighed soil properties are |
---|
1005 | !! calculated\n |
---|
1006 | !! - The coefficients cgrnd and dgrnd are the integration |
---|
1007 | !! coefficients for the thermal scheme \n |
---|
1008 | !! T(k+1)=cgrnd(k)+dgrnd(k)*T(k) \n |
---|
1009 | !! -- EQ1 -- \n |
---|
1010 | !! They correspond respectively to $\beta$ and $\alpha$ from F. Hourdin\'s thesis and |
---|
1011 | !! their expression can be found in this document (eq A19 and A20) |
---|
1012 | !! - soilcap and soilflux are the apparent surface heat capacity and flux |
---|
1013 | !! used in enerbil at the next timestep to solve the surface |
---|
1014 | !! balance for Ts (EQ3); they correspond to $C_s$ and $F_s$ in F. |
---|
1015 | !! Hourdin\'s PhD thesis and are expressed in eq. A30 and A31. \n |
---|
1016 | !! soilcap*(Ts(t)-Ts(t-1))/dt=soilflux+otherfluxes(Ts(t)) \n |
---|
1017 | !! -- EQ3 --\n |
---|
1018 | !! |
---|
1019 | !! RECENT CHANGE(S) : None |
---|
1020 | !! |
---|
1021 | !! MAIN OUTPUT VARIABLE(S): cgrnd, dgrnd, pcapa, pkappa, soilcap, soilflx |
---|
1022 | !! |
---|
1023 | !! REFERENCE(S) : |
---|
1024 | !! - Hourdin, F. (1992). Study and numerical simulation of the general circulation of planetary atmospheres, |
---|
1025 | !! Ph.D. thesis, Paris VII University. Remark: the part of F. Hourdin's PhD thesis relative to the thermal |
---|
1026 | !! integration scheme has been scanned and is provided along with the documentation, with name : |
---|
1027 | !! Hourdin_1992_PhD_thermal_scheme.pdf |
---|
1028 | !! |
---|
1029 | !! FLOWCHART : None |
---|
1030 | !! \n |
---|
1031 | !_ ================================================================================================================================ |
---|
1032 | |
---|
1033 | SUBROUTINE thermosoil_coef (kjpindex, dtradia, temp_sol_new, snow, ptn, soilcap, soilflx, zz, dz1, dz2, z1, zdz1,& |
---|
1034 | & zdz2, cgrnd, dgrnd, pcapa, pcapa_en, pkappa) |
---|
1035 | |
---|
1036 | !! 0. Variables and parameter declaration |
---|
1037 | |
---|
1038 | !! 0.1 Input variables |
---|
1039 | |
---|
1040 | INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (unitless) |
---|
1041 | REAL(r_std), INTENT (in) :: dtradia !! Time step in seconds @tex ($s$) @endtex |
---|
1042 | REAL(r_std), DIMENSION (kjpindex), INTENT (in) :: temp_sol_new !! soil surface temperature @tex ($K$) @endtex |
---|
1043 | REAL(r_std), DIMENSION (kjpindex), INTENT (in) :: snow !! snow mass @tex ($Kg$) @endtex |
---|
1044 | REAL(r_std), DIMENSION (ngrnd), INTENT(in) :: zz !! depths of the soil thermal numerical nodes |
---|
1045 | !! @tex ($m$) @endtex |
---|
1046 | REAL(r_std), DIMENSION (ngrnd), INTENT(in) :: dz1 !! numerical constant depending on the vertical |
---|
1047 | !! thermal grid only @tex ($m^{-1}$) @endtex |
---|
1048 | REAL(r_std), DIMENSION (ngrnd), INTENT(in) :: dz2 !! thicknesses of the soil thermal layers |
---|
1049 | !! @tex ($m$) @endtex |
---|
1050 | REAL(r_std), DIMENSION (kjpindex,ngrnd), INTENT (in) :: ptn !! vertically discretized soil temperatures |
---|
1051 | !! @tex ($K$) @endtex |
---|
1052 | |
---|
1053 | !! 0.2 Output variables |
---|
1054 | |
---|
1055 | REAL(r_std), DIMENSION (kjpindex), INTENT (out) :: soilcap !! surface heat capacity |
---|
1056 | !! @tex ($J m^{-2} K^{-1}$) @endtex |
---|
1057 | REAL(r_std), DIMENSION (kjpindex), INTENT (out) :: soilflx !! surface heat flux @tex ($W m^{-2}$) @endtex, |
---|
1058 | !! positive towards the |
---|
1059 | !! soil, writen as Qg (ground heat flux) in the history |
---|
1060 | !! files. |
---|
1061 | REAL(r_std), DIMENSION (kjpindex), INTENT (out) :: z1 !! numerical constant @tex ($W m^{-1} K^{-1}$) @endtex |
---|
1062 | |
---|
1063 | REAL(r_std), DIMENSION (kjpindex,ngrnd-1), INTENT(out) :: cgrnd !! matrix coefficient for the computation of soil |
---|
1064 | !! temperatures (beta in F. Hourdin thesis) |
---|
1065 | REAL(r_std), DIMENSION (kjpindex,ngrnd-1), INTENT(out) :: dgrnd !! matrix coefficient for the computation of soil |
---|
1066 | !! temperatures (alpha in F. Hourdin thesis) |
---|
1067 | REAL(r_std), DIMENSION (kjpindex,ngrnd-1), INTENT(out) :: zdz1 !! numerical (buffer) constant |
---|
1068 | !! @tex ($W m^{-1} K^{-1}$) @endtex |
---|
1069 | |
---|
1070 | REAL(r_std), DIMENSION (kjpindex,ngrnd), INTENT(out) :: zdz2 !! numerical (buffer) constant |
---|
1071 | !! @tex ($W m^{-1} K^{-1}$) @endtex |
---|
1072 | |
---|
1073 | |
---|
1074 | !! 0.3 Modified variable |
---|
1075 | |
---|
1076 | REAL(r_std), DIMENSION (kjpindex,ngrnd), INTENT(inout) :: pcapa !! volumetric vertically discretized soil heat capacity |
---|
1077 | !! @tex ($J K^{-1} m^{-3}$) @endtex |
---|
1078 | REAL(r_std), DIMENSION (kjpindex,ngrnd), INTENT(inout) :: pcapa_en !! volumetric vertically discretized heat capacity used |
---|
1079 | !! to calculate surfheat_incr |
---|
1080 | !! @tex ($J K^{-1} m^{-3}$) @endtex |
---|
1081 | REAL(r_std), DIMENSION (kjpindex,ngrnd), INTENT(inout) :: pkappa !! vertically discretized soil thermal conductivity |
---|
1082 | !! @tex ($W m^{-1} K^{-1}$) @endtex |
---|
1083 | |
---|
1084 | !! 0.4 Local variables |
---|
1085 | |
---|
1086 | INTEGER(i_std) :: ji, jg |
---|
1087 | REAL(r_std), DIMENSION(kjpindex) :: snow_h !! snow_h is the snow height @tex ($m$) @endtex |
---|
1088 | REAL(r_std), DIMENSION(kjpindex) :: zx1, zx2 !! zx1 and zx2 are the layer fraction consisting in snow |
---|
1089 | !! and soil respectively. |
---|
1090 | !_ ================================================================================================================================ |
---|
1091 | |
---|
1092 | !! 1. Computation of the soil thermal properties |
---|
1093 | |
---|
1094 | ! Computation of the soil thermal properties; snow properties are also accounted for |
---|
1095 | DO ji = 1,kjpindex |
---|
1096 | snow_h(ji) = snow(ji) / sn_dens |
---|
1097 | |
---|
1098 | IF ( snow_h(ji) .GT. zz_coef(1) ) THEN |
---|
1099 | pcapa(ji,1) = sn_capa |
---|
1100 | pcapa_en(ji,1) = sn_capa |
---|
1101 | pkappa(ji,1) = sn_cond |
---|
1102 | ELSE IF ( snow_h(ji) .GT. zero ) THEN |
---|
1103 | pcapa_en(ji,1) = sn_capa |
---|
1104 | zx1(ji) = snow_h(ji) / zz_coef(1) |
---|
1105 | zx2(ji) = ( zz_coef(1) - snow_h(ji)) / zz_coef(1) |
---|
1106 | pcapa(ji,1) = zx1(ji) * sn_capa + zx2(ji) * so_capa_wet |
---|
1107 | pkappa(ji,1) = un / ( zx1(ji) / sn_cond + zx2(ji) / so_cond_wet ) |
---|
1108 | ELSE |
---|
1109 | pcapa(ji,1) = so_capa_dry + wetdiag(ji,1)*(so_capa_wet - so_capa_dry) |
---|
1110 | pkappa(ji,1) = so_cond_dry + wetdiag(ji,1)*(so_cond_wet - so_cond_dry) |
---|
1111 | pcapa_en(ji,1) = so_capa_dry + wetdiag(ji,1)*(so_capa_wet - so_capa_dry) |
---|
1112 | ENDIF |
---|
1113 | ! |
---|
1114 | DO jg = 2, ngrnd - 2 |
---|
1115 | IF ( snow_h(ji) .GT. zz_coef(jg) ) THEN |
---|
1116 | pcapa(ji,jg) = sn_capa |
---|
1117 | pkappa(ji,jg) = sn_cond |
---|
1118 | pcapa_en(ji,jg) = sn_capa |
---|
1119 | ELSE IF ( snow_h(ji) .GT. zz_coef(jg-1) ) THEN |
---|
1120 | zx1(ji) = (snow_h(ji) - zz_coef(jg-1)) / (zz_coef(jg) - zz_coef(jg-1)) |
---|
1121 | zx2(ji) = ( zz_coef(jg) - snow_h(ji)) / (zz_coef(jg) - zz_coef(jg-1)) |
---|
1122 | pcapa(ji,jg) = zx1(ji) * sn_capa + zx2(ji) * so_capa_wet |
---|
1123 | pkappa(ji,jg) = un / ( zx1(ji) / sn_cond + zx2(ji) / so_cond_wet ) |
---|
1124 | pcapa_en(ji,jg) = sn_capa |
---|
1125 | ELSE |
---|
1126 | pcapa(ji,jg) = so_capa_dry + wetdiag(ji,jg)*(so_capa_wet - so_capa_dry) |
---|
1127 | pkappa(ji,jg) = so_cond_dry + wetdiag(ji,jg)*(so_cond_wet - so_cond_dry) |
---|
1128 | pcapa_en(ji,jg) = so_capa_dry + wetdiag(ji,jg)*(so_capa_wet - so_capa_dry) |
---|
1129 | ENDIF |
---|
1130 | ENDDO |
---|
1131 | |
---|
1132 | ENDDO |
---|
1133 | |
---|
1134 | !! 2. computation of the coefficients of the numerical integration scheme |
---|
1135 | |
---|
1136 | ! cgrnd, dgrnd |
---|
1137 | |
---|
1138 | !! 2.1. some "buffer" values |
---|
1139 | DO jg=1,ngrnd |
---|
1140 | DO ji=1,kjpindex |
---|
1141 | zdz2(ji,jg)=pcapa(ji,jg) * dz2(jg)/dtradia |
---|
1142 | ENDDO |
---|
1143 | ENDDO |
---|
1144 | |
---|
1145 | DO jg=1,ngrnd-1 |
---|
1146 | DO ji=1,kjpindex |
---|
1147 | zdz1(ji,jg) = dz1(jg) * pkappa(ji,jg) |
---|
1148 | ENDDO |
---|
1149 | ENDDO |
---|
1150 | |
---|
1151 | !! 2.2. the coefficients ! |
---|
1152 | DO ji = 1,kjpindex |
---|
1153 | z1(ji) = zdz2(ji,ngrnd) + zdz1(ji,ngrnd-1) |
---|
1154 | cgrnd(ji,ngrnd-1) = zdz2(ji,ngrnd) * ptn(ji,ngrnd) / z1(ji) |
---|
1155 | dgrnd(ji,ngrnd-1) = zdz1(ji,ngrnd-1) / z1(ji) |
---|
1156 | ENDDO |
---|
1157 | |
---|
1158 | DO jg = ngrnd-1,2,-1 |
---|
1159 | DO ji = 1,kjpindex |
---|
1160 | z1(ji) = un / (zdz2(ji,jg) + zdz1(ji,jg-1) + zdz1(ji,jg) * (un - dgrnd(ji,jg))) |
---|
1161 | cgrnd(ji,jg-1) = (ptn(ji,jg) * zdz2(ji,jg) + zdz1(ji,jg) * cgrnd(ji,jg)) * z1(ji) |
---|
1162 | dgrnd(ji,jg-1) = zdz1(ji,jg-1) * z1(ji) |
---|
1163 | ENDDO |
---|
1164 | ENDDO |
---|
1165 | |
---|
1166 | !! 3. Computation of the apparent ground heat flux |
---|
1167 | |
---|
1168 | !! Computation of the apparent ground heat flux (> towards the soil) and |
---|
1169 | !! apparent surface heat capacity, used at the next timestep by enerbil to |
---|
1170 | !! compute the surface temperature. |
---|
1171 | DO ji = 1,kjpindex |
---|
1172 | soilflx(ji) = zdz1(ji,1) * (cgrnd(ji,1) + (dgrnd(ji,1)-1.) * ptn(ji,1)) |
---|
1173 | soilcap(ji) = (zdz2(ji,1) * dtradia + dtradia * (un - dgrnd(ji,1)) * zdz1(ji,1)) |
---|
1174 | z1(ji) = lambda * (un - dgrnd(ji,1)) + un |
---|
1175 | soilcap(ji) = soilcap(ji) / z1(ji) |
---|
1176 | soilflx(ji) = soilflx(ji) + & |
---|
1177 | & soilcap(ji) * (ptn(ji,1) * z1(ji) - lambda * cgrnd(ji,1) - temp_sol_new(ji)) / dtradia |
---|
1178 | ENDDO |
---|
1179 | |
---|
1180 | IF (long_print) WRITE (numout,*) ' thermosoil_coef done ' |
---|
1181 | |
---|
1182 | END SUBROUTINE thermosoil_coef |
---|
1183 | |
---|
1184 | |
---|
1185 | !! ================================================================================================================================ |
---|
1186 | !! SUBROUTINE : thermosoil_profile |
---|
1187 | !! |
---|
1188 | !>\BRIEF In this routine solves the numerical soil thermal scheme, ie calculates the new soil temperature profile; |
---|
1189 | !! This profile is then exported onto the diagnostic axis (call thermosoil_diaglev) |
---|
1190 | !! |
---|
1191 | !! DESCRIPTION : The calculation of the new soil temperature profile is based on |
---|
1192 | !! the cgrnd and dgrnd values from the previous timestep and the surface temperature Ts aka temp_sol_new. (see detailed |
---|
1193 | !! explanation in the header of the thermosoil module or in the reference).\n |
---|
1194 | !! T(k+1)=cgrnd(k)+dgrnd(k)*T(k)\n |
---|
1195 | !! -- EQ1 --\n |
---|
1196 | !! Ts=(1-lambda)*T(1) -lambda*T(2)\n |
---|
1197 | !! -- EQ2--\n |
---|
1198 | !! |
---|
1199 | !! RECENT CHANGE(S) : None |
---|
1200 | !! |
---|
1201 | !! MAIN OUTPUT VARIABLE(S): ptn (soil temperature profile on the thermal axis), |
---|
1202 | !! stempdiag (soil temperature profile on the diagnostic axis) |
---|
1203 | !! |
---|
1204 | !! REFERENCE(S) : |
---|
1205 | !! - Hourdin, F. (1992). Study and numerical simulation of the general circulation of planetary atmospheres, |
---|
1206 | !! Ph.D. thesis, Paris VII University. Remark: the part of F. Hourdin's PhD thesis relative to the thermal |
---|
1207 | !! integration scheme has been scanned and is provided along with the documentation, with name : |
---|
1208 | !! Hourdin_1992_PhD_thermal_scheme.pdf |
---|
1209 | !! |
---|
1210 | !! FLOWCHART : None |
---|
1211 | !! \n |
---|
1212 | !_ ================================================================================================================================ |
---|
1213 | |
---|
1214 | SUBROUTINE thermosoil_profile (kjpindex, temp_sol_new, ptn, stempdiag) |
---|
1215 | |
---|
1216 | !! 0. Variables and parameter declaration |
---|
1217 | |
---|
1218 | !! 0.1 Input variables |
---|
1219 | |
---|
1220 | INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (unitless) |
---|
1221 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_sol_new !! Surface temperature at the present time-step |
---|
1222 | !! @tex ($K$) @endtex |
---|
1223 | |
---|
1224 | !! 0.2 Output variables |
---|
1225 | REAL(r_std),DIMENSION (kjpindex,nbdl), INTENT (out) :: stempdiag !! diagnostic temperature profile |
---|
1226 | !! @tex ($K$) @endtex |
---|
1227 | |
---|
1228 | !! 0.3 Modified variables |
---|
1229 | |
---|
1230 | REAL(r_std),DIMENSION (kjpindex,ngrnd), INTENT (inout) :: ptn !! vertically discretized soil temperatures |
---|
1231 | !! @tex ($K$) @endtex |
---|
1232 | |
---|
1233 | |
---|
1234 | !! 0.4 Local variables |
---|
1235 | |
---|
1236 | INTEGER(i_std) :: ji, jg |
---|
1237 | !_ ================================================================================================================================ |
---|
1238 | |
---|
1239 | !! 1. Computes the soil temperatures ptn. |
---|
1240 | |
---|
1241 | !! 1.1. ptn(jg=1) using EQ1 and EQ2 |
---|
1242 | DO ji = 1,kjpindex |
---|
1243 | ptn(ji,1) = (lambda * cgrnd(ji,1) + temp_sol_new(ji)) / (lambda * (un - dgrnd(ji,1)) + un) |
---|
1244 | ENDDO |
---|
1245 | |
---|
1246 | !! 1.2. ptn(jg=2:ngrnd) using EQ1. |
---|
1247 | DO jg = 1,ngrnd-1 |
---|
1248 | DO ji = 1,kjpindex |
---|
1249 | ptn(ji,jg+1) = cgrnd(ji,jg) + dgrnd(ji,jg) * ptn(ji,jg) |
---|
1250 | ENDDO |
---|
1251 | ENDDO |
---|
1252 | |
---|
1253 | !! 2. Put the soil temperatures onto the diagnostic axis |
---|
1254 | |
---|
1255 | !! Put the soil temperatures onto the diagnostic axis for convenient |
---|
1256 | !! use in other routines (stomate..) |
---|
1257 | CALL thermosoil_diaglev(kjpindex, stempdiag) |
---|
1258 | |
---|
1259 | IF (long_print) WRITE (numout,*) ' thermosoil_profile done ' |
---|
1260 | |
---|
1261 | END SUBROUTINE thermosoil_profile |
---|
1262 | |
---|
1263 | |
---|
1264 | !! ================================================================================================================================ |
---|
1265 | !! SUBROUTINE : thermosoil_diaglev |
---|
1266 | !! |
---|
1267 | !>\BRIEF Interpolation of the soil in-depth temperatures onto the diagnostic profile. |
---|
1268 | !! |
---|
1269 | !! DESCRIPTION : This is a very easy linear interpolation, with intfact(jd, jg) the fraction |
---|
1270 | !! the thermal layer jg comprised within the diagnostic layer jd. The depths of |
---|
1271 | !! the diagnostic levels are diaglev(1:nbdl), computed in slowproc.f90. |
---|
1272 | !! |
---|
1273 | !! RECENT CHANGE(S) : None |
---|
1274 | !! |
---|
1275 | !! MAIN OUTPUT VARIABLE(S): stempdiag (soil temperature profile on the diagnostic axis) |
---|
1276 | !! |
---|
1277 | !! REFERENCE(S) : None |
---|
1278 | !! |
---|
1279 | !! FLOWCHART : None |
---|
1280 | !! \n |
---|
1281 | !_ ================================================================================================================================ |
---|
1282 | |
---|
1283 | SUBROUTINE thermosoil_diaglev(kjpindex, stempdiag) |
---|
1284 | |
---|
1285 | !! 0. Variables and parameter declaration |
---|
1286 | |
---|
1287 | !! 0.1 Input variables |
---|
1288 | |
---|
1289 | INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (unitless) |
---|
1290 | |
---|
1291 | !! 0.2 Output variables |
---|
1292 | |
---|
1293 | REAL(r_std),DIMENSION (kjpindex,nbdl), INTENT (out) :: stempdiag !! Diagnostoc soil temperature profile @tex ($K$) @endtex |
---|
1294 | |
---|
1295 | !! 0.3 Modified variables |
---|
1296 | |
---|
1297 | !! 0.4 Local variables |
---|
1298 | |
---|
1299 | INTEGER(i_std) :: ji, jd, jg |
---|
1300 | REAL(r_std) :: lev_diag, prev_diag, lev_prog, prev_prog |
---|
1301 | REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: intfact |
---|
1302 | LOGICAL, PARAMETER :: check=.FALSE. |
---|
1303 | !_ ================================================================================================================================ |
---|
1304 | |
---|
1305 | !! 1. Computes intfact(jd, jg) |
---|
1306 | |
---|
1307 | !! Computes intfact(jd, jg), the fraction |
---|
1308 | !! the thermal layer jg comprised within the diagnostic layer jd. |
---|
1309 | IF ( .NOT. ALLOCATED(intfact)) THEN |
---|
1310 | |
---|
1311 | ALLOCATE(intfact(nbdl, ngrnd)) |
---|
1312 | |
---|
1313 | prev_diag = zero |
---|
1314 | DO jd = 1, nbdl |
---|
1315 | lev_diag = diaglev(jd) |
---|
1316 | prev_prog = zero |
---|
1317 | DO jg = 1, ngrnd |
---|
1318 | IF ( jg == ngrnd .AND. (prev_prog + dz2(jg)) < lev_diag ) THEN |
---|
1319 | lev_prog = lev_diag |
---|
1320 | ELSE |
---|
1321 | lev_prog = prev_prog + dz2(jg) |
---|
1322 | ENDIF |
---|
1323 | intfact(jd,jg) = MAX(MIN(lev_diag,lev_prog)-MAX(prev_diag, prev_prog), zero)/(lev_diag-prev_diag) |
---|
1324 | prev_prog = lev_prog |
---|
1325 | ENDDO |
---|
1326 | prev_diag = lev_diag |
---|
1327 | ENDDO |
---|
1328 | |
---|
1329 | IF ( check ) THEN |
---|
1330 | WRITE(numout,*) 'thermosoil_diagev -- thermosoil_diaglev -- thermosoil_diaglev --' |
---|
1331 | DO jd = 1, nbdl |
---|
1332 | WRITE(numout,*) jd, '-', intfact(jd,1:ngrnd) |
---|
1333 | ENDDO |
---|
1334 | WRITE(numout,*) "SUM -- SUM -- SUM SUM -- SUM -- SUM" |
---|
1335 | DO jd = 1, nbdl |
---|
1336 | WRITE(numout,*) jd, '-', SUM(intfact(jd,1:ngrnd)) |
---|
1337 | ENDDO |
---|
1338 | WRITE(numout,*) 'thermosoil_diaglev -- thermosoil_diaglev -- thermosoil_diaglev --' |
---|
1339 | ENDIF |
---|
1340 | |
---|
1341 | ENDIF |
---|
1342 | |
---|
1343 | !! 2. does the interpolation |
---|
1344 | |
---|
1345 | stempdiag(:,:) = zero |
---|
1346 | DO jg = 1, ngrnd |
---|
1347 | DO jd = 1, nbdl |
---|
1348 | DO ji = 1, kjpindex |
---|
1349 | stempdiag(ji,jd) = stempdiag(ji,jd) + ptn(ji,jg)*intfact(jd,jg) |
---|
1350 | ENDDO |
---|
1351 | ENDDO |
---|
1352 | ENDDO |
---|
1353 | |
---|
1354 | END SUBROUTINE thermosoil_diaglev |
---|
1355 | |
---|
1356 | |
---|
1357 | !! ================================================================================================================================ |
---|
1358 | !! SUBROUTINE : thermosoil_humlev |
---|
1359 | !! |
---|
1360 | !>\BRIEF Interpolates the diagnostic soil humidity profile shumdiag(nbdl, diagnostic axis) onto |
---|
1361 | !! the thermal axis, which gives wetdiag(ngrnd, thermal axis). |
---|
1362 | !! |
---|
1363 | !! DESCRIPTION : Same as in thermosoil_diaglev : This is a very easy linear interpolation, with intfactw(jd, jg) the fraction |
---|
1364 | !! the thermal layer jd comprised within the diagnostic layer jg. |
---|
1365 | !!?? I would think wise to change the indeces here, to keep jD for Diagnostic |
---|
1366 | !!?? and jG for Ground thermal levels... |
---|
1367 | !! |
---|
1368 | !! The depths of the diagnostic levels are diaglev(1:nbdl), computed in slowproc.f90. |
---|
1369 | !! Recall that when the 11-layer hydrology is used, |
---|
1370 | !! wetdiag and shumdiag are with reference to the moisture content (mc) |
---|
1371 | !! at the wilting point mcw : wetdiag=(mc-mcw)/(mcs-mcw). |
---|
1372 | !! with mcs the saturated soil moisture content. |
---|
1373 | !! |
---|
1374 | !! RECENT CHANGE(S) : None |
---|
1375 | !! |
---|
1376 | !! MAIN OUTPUT VARIABLE(S): wetdiag (soil soil humidity profile on the thermal axis) |
---|
1377 | !! |
---|
1378 | !! REFERENCE(S) : None |
---|
1379 | !! |
---|
1380 | !! FLOWCHART : None |
---|
1381 | !! \n |
---|
1382 | !_ ================================================================================================================================ |
---|
1383 | |
---|
1384 | SUBROUTINE thermosoil_humlev(kjpindex, shumdiag) |
---|
1385 | |
---|
1386 | !! 0. Variables and parameter declaration |
---|
1387 | |
---|
1388 | !! 0.1 Input variables |
---|
1389 | |
---|
1390 | INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (unitless) |
---|
1391 | REAL(r_std),DIMENSION (kjpindex,nbdl), INTENT (in) :: shumdiag !! Relative soil humidity on the diagnostic axis. |
---|
1392 | !! (0-1, unitless). Caveats : when "hydrol" (the 11-layers |
---|
1393 | !! hydrology) is used, this humidity is calculated with |
---|
1394 | !! respect to the wilting point : |
---|
1395 | !! shumdiag= (mc-mcw)/(mcs-mcw), with mc : moisture |
---|
1396 | !! content; mcs : saturated soil moisture content; mcw: |
---|
1397 | !! soil moisture content at the wilting point. when the 2-layers |
---|
1398 | !! hydrology "hydrolc" is used, shumdiag is just |
---|
1399 | !! a diagnostic humidity index, with no real physical |
---|
1400 | !! meaning. |
---|
1401 | |
---|
1402 | !! 0.2 Output variables |
---|
1403 | |
---|
1404 | !! 0.3 Modified variables |
---|
1405 | |
---|
1406 | !! 0.4 Local variables |
---|
1407 | INTEGER(i_std) :: ji, jd, jg |
---|
1408 | REAL(r_std) :: lev_diag, prev_diag, lev_prog, prev_prog |
---|
1409 | REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: intfactw !! fraction of each diagnostic layer (jd) comprized within |
---|
1410 | !! a given thermal layer (jg)(0-1, unitless) |
---|
1411 | LOGICAL, PARAMETER :: check=.FALSE. |
---|
1412 | !_ ================================================================================================================================ |
---|
1413 | |
---|
1414 | !! 1. computes intfactw(jd,jg), the fraction of each diagnostic layer (jg) comprized within a given thermal layer (jd) |
---|
1415 | |
---|
1416 | IF ( .NOT. ALLOCATED(intfactw)) THEN |
---|
1417 | |
---|
1418 | ALLOCATE(intfactw(ngrnd, nbdl)) |
---|
1419 | |
---|
1420 | prev_diag = zero |
---|
1421 | DO jd = 1, ngrnd |
---|
1422 | lev_diag = prev_diag + dz2(jd) |
---|
1423 | prev_prog = zero |
---|
1424 | DO jg = 1, nbdl |
---|
1425 | IF ( jg == nbdl .AND. diaglev(jg) < lev_diag ) THEN |
---|
1426 | lev_prog = lev_diag |
---|
1427 | ELSE |
---|
1428 | lev_prog = diaglev(jg) |
---|
1429 | ENDIF |
---|
1430 | intfactw(jd,jg) = MAX(MIN(lev_diag,lev_prog)-MAX(prev_diag, prev_prog), zero)/(lev_diag-prev_diag) |
---|
1431 | prev_prog = lev_prog |
---|
1432 | ENDDO |
---|
1433 | prev_diag = lev_diag |
---|
1434 | ENDDO |
---|
1435 | |
---|
1436 | IF ( check ) THEN |
---|
1437 | WRITE(numout,*) 'thermosoil_humlev -- thermosoil_humlev -- thermosoil_humlev --' |
---|
1438 | DO jd = 1, ngrnd |
---|
1439 | WRITE(numout,*) jd, '-', intfactw(jd,1:nbdl) |
---|
1440 | ENDDO |
---|
1441 | WRITE(numout,*) "SUM -- SUM -- SUM SUM -- SUM -- SUM" |
---|
1442 | DO jd = 1, ngrnd |
---|
1443 | WRITE(numout,*) jd, '-', SUM(intfactw(jd,1:nbdl)) |
---|
1444 | ENDDO |
---|
1445 | WRITE(numout,*) 'thermosoil_humlev -- thermosoil_humlev -- thermosoil_humlev --' |
---|
1446 | ENDIF |
---|
1447 | |
---|
1448 | ENDIF |
---|
1449 | |
---|
1450 | !! 2. does the interpolation |
---|
1451 | |
---|
1452 | wetdiag(:,:) = zero |
---|
1453 | DO jg = 1, nbdl |
---|
1454 | DO jd = 1, ngrnd |
---|
1455 | DO ji = 1, kjpindex |
---|
1456 | wetdiag(ji,jd) = wetdiag(ji,jd) + shumdiag(ji,jg)*intfactw(jd,jg) |
---|
1457 | ENDDO |
---|
1458 | ENDDO |
---|
1459 | ENDDO |
---|
1460 | |
---|
1461 | END SUBROUTINE thermosoil_humlev |
---|
1462 | |
---|
1463 | |
---|
1464 | !! ================================================================================================================================ |
---|
1465 | !! SUBROUTINE : thermosoil_energy |
---|
1466 | !! |
---|
1467 | !>\BRIEF Energy check-up. |
---|
1468 | !! |
---|
1469 | !! DESCRIPTION : I didn\'t comment this routine since at do not understand its use, please |
---|
1470 | !! ask initial designers (Jan ? Nathalie ?). |
---|
1471 | !! |
---|
1472 | !! RECENT CHANGE(S) : None |
---|
1473 | !! |
---|
1474 | !! MAIN OUTPUT VARIABLE(S) : ?? |
---|
1475 | !! |
---|
1476 | !! REFERENCE(S) : None |
---|
1477 | !! |
---|
1478 | !! FLOWCHART : None |
---|
1479 | !! \n |
---|
1480 | !_ ================================================================================================================================ |
---|
1481 | |
---|
1482 | SUBROUTINE thermosoil_energy(kjpindex, temp_sol_new, soilcap, first_call) |
---|
1483 | |
---|
1484 | !! 0. Variables and parameter declaration |
---|
1485 | |
---|
1486 | !! 0.1 Input variables |
---|
1487 | |
---|
1488 | INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (unitless) |
---|
1489 | LOGICAL, INTENT (in) :: first_call !! First call (true/false) |
---|
1490 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_sol_new !! Surface temperature at the present time-step, Ts |
---|
1491 | !! @tex ($K$) @endtex |
---|
1492 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: soilcap !! Apparent surface heat capacity |
---|
1493 | !! @tex ($J m^{-2} K^{-1}$) @endtex, |
---|
1494 | !! see eq. A29 of F. Hourdin\'s PhD thesis. |
---|
1495 | |
---|
1496 | !! 0.2 Output variables |
---|
1497 | |
---|
1498 | !! 0.3 Modified variables |
---|
1499 | |
---|
1500 | !! 0.4 Local variables |
---|
1501 | |
---|
1502 | INTEGER(i_std) :: ji, jg |
---|
1503 | !_ ================================================================================================================================ |
---|
1504 | |
---|
1505 | IF (first_call) THEN |
---|
1506 | |
---|
1507 | DO ji = 1, kjpindex |
---|
1508 | surfheat_incr(ji) = zero |
---|
1509 | coldcont_incr(ji) = zero |
---|
1510 | temp_sol_beg(ji) = temp_sol_new(ji) |
---|
1511 | |
---|
1512 | DO jg = 1, ngrnd |
---|
1513 | ptn_beg(ji,jg) = ptn(ji,jg) |
---|
1514 | ENDDO |
---|
1515 | |
---|
1516 | ENDDO |
---|
1517 | |
---|
1518 | RETURN |
---|
1519 | |
---|
1520 | ENDIF |
---|
1521 | |
---|
1522 | DO ji = 1, kjpindex |
---|
1523 | surfheat_incr(ji) = zero |
---|
1524 | coldcont_incr(ji) = zero |
---|
1525 | ENDDO |
---|
1526 | |
---|
1527 | ! Sum up the energy content of all layers in the soil. |
---|
1528 | DO ji = 1, kjpindex |
---|
1529 | |
---|
1530 | IF (pcapa_en(ji,1) .LE. sn_capa) THEN |
---|
1531 | |
---|
1532 | ! Verify the energy conservation in the surface layer |
---|
1533 | coldcont_incr(ji) = soilcap(ji) * (temp_sol_new(ji) - temp_sol_beg(ji)) |
---|
1534 | surfheat_incr(ji) = zero |
---|
1535 | ELSE |
---|
1536 | |
---|
1537 | ! Verify the energy conservation in the surface layer |
---|
1538 | surfheat_incr(ji) = soilcap(ji) * (temp_sol_new(ji) - temp_sol_beg(ji)) |
---|
1539 | coldcont_incr(ji) = zero |
---|
1540 | ENDIF |
---|
1541 | ENDDO |
---|
1542 | |
---|
1543 | ptn_beg(:,:) = ptn(:,:) |
---|
1544 | temp_sol_beg(:) = temp_sol_new(:) |
---|
1545 | |
---|
1546 | END SUBROUTINE thermosoil_energy |
---|
1547 | |
---|
1548 | END MODULE thermosoil |
---|