Context Navigation

source: branches/publications/ORCHIDEE_CAN_r2290/src_sechiba/thermosoil.f90 @ 7475

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

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

Download in other formats: