Context Navigation

thermosoil.f90 @ 854

Last change on this file since 854 was 816, checked in by fabienne.maignan, 12 years ago
Correction of a single !
Property svn:keywords set to `Revision Date HeadURL Date Author Revision`
File size: 76.2 KB

Line
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) -lambdaT(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)-lambdaT(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) -lambdaT(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

Note: See TracBrowser for help on using the repository browser.

Context Navigation

source: perso/abdelouhab.djerrah/ORCHIDEE/src_sechiba/thermosoil.f90 @ 854

Download in other formats: