Context Navigation

source: tags/ORCHIDEE_2_0/ORCHIDEE/src_sechiba/thermosoil.f90

Last change on this file was 5096, checked in by josefine.ghattas, 6 years ago
Protected calculation and output of snowtemp_weighted done only for explicitsnow.
Property svn:keywords set to `Revision Date HeadURL Date Author Revision`
File size: 104.6 KB

Line
1	! =================================================================================================================================
2	! MODULE : thermosoil
3	!
4	! CONTACT : orchidee-help _at_ listes.ipsl.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. This module is only used with CWRR hydrology.
11	!!
12	!!\n DESCRIPTION : General important informations about the numerical scheme and
13	!! the soil vertical discretization:\n
14	!! - the soil is zmaxt deep (by default 10m) and divided into "ngrnd" layers.
15	!! From 0-zmaxh(default 2m), the discretization is the same as for hydrology.
16	!! From zmaxh(2m) and below, the depth increase linearly (by default) or geometrically. \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=soilflx+otherfluxes(Ts(t))\n
46	!! -- EQ3 --\n
47	!!
48	!! and the dependency from the previous timestep is hidden in
49	!! soilcap and soilflx (apparent surface heat capacity and heat
50	!! flux respectively). Soilcap and soilflx are therefore
51	!! calculated at the previous 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 = (zlt(1))/((zlt(2)-zlt(1))) \n
60	!!
61	!! RECENT CHANGE(S) : - Change soil thermal properties to consider also soil texture, rev 2922.
62	!! - Change vertical discretization, rev 2917. Note: In the revised thermosoil,
63	!! cstgrnd and lskin are not needed any more. The depth znt, zlt and dlt
64	!! are computed in vertical_soil and are in meter
65	!!
66	!! REFERENCE(S) : None
67	!!
68	!! SVN :
69	!! $HeadURL$
70	!! $Date$
71	!! $Revision$
72	!! \n
73	!_ ================================================================================================================================
74
75	MODULE thermosoil
76
77	! modules used :
78	USE ioipsl
79	USE ioipsl_para
80	USE xios_orchidee
81	USE constantes
82	USE time, ONLY : one_day, dt_sechiba
83	USE constantes_soil
84	USE sechiba_io_p
85	USE grid
86
87	IMPLICIT NONE
88
89	!private and public routines :
90	PRIVATE
91	PUBLIC :: thermosoil_main, thermosoil_clear, thermosoil_initialize, thermosoil_finalize
92
93	REAL(r_std), SAVE :: lambda !! See Module description
94	!$OMP THREADPRIVATE(lambda)
95	REAL(r_std), SAVE :: fz1, zalph !! usefull constants for diverse use
96	!$OMP THREADPRIVATE(fz1, zalph)
97	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: ptn !! vertically discretized
98	!! soil temperatures @tex ($K$) @endtex.
99	!$OMP THREADPRIVATE(ptn)
100	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: dz1 !! numerical constant used in the thermal numerical
101	!! scheme @tex ($m^{-1}$) @endtex. ; it corresponds
102	!! to the coefficient @tex $d_k$ @endtex of equation
103	!! (A.12) in F. Hourdin PhD thesis.
104	!$OMP THREADPRIVATE(dz1)
105	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: cgrnd !! integration coefficient for the numerical scheme,
106	!! see eq.1
107	!$OMP THREADPRIVATE(cgrnd)
108	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: dgrnd !! integration coefficient for the numerical scheme,
109	!! see eq.1
110	!$OMP THREADPRIVATE(dgrnd)
111
112	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: pcapa !! volumetric vertically discretized soil heat
113	!! capacity @tex ($J K^{-1} m^{-3}$) @endtex.
114	!! It depends on the soil
115	!! moisture content (shum_ngrnd_perma) and is calculated at
116	!! each time step in thermosoil_coef.
117	!$OMP THREADPRIVATE(pcapa)
118	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: pkappa !! vertically discretized soil thermal conductivity
119	!! @tex ($W K^{-1} m^{-1}$) @endtex. Same as pcapa.
120	!$OMP THREADPRIVATE(pkappa)
121	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: pcapa_snow !! volumetric vertically discretized snow heat
122	!! capacity @tex ($J K^{-1} m^{-3}$) @endtex.
123	!$OMP THREADPRIVATE(pcapa_snow)
124	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: pkappa_snow !! vertically discretized snow thermal conductivity
125	!! @tex ($W K^{-1} m^{-1}$) @endtex.
126	!$OMP THREADPRIVATE(pkappa_snow)
127
128	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: pcapa_en !! heat capacity used for surfheat_incr and
129	!! coldcont_incr
130	!$OMP THREADPRIVATE(pcapa_en)
131	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: ptn_beg !! ptn as it is after thermosoil_profile but before thermosoil_coef,
132	!! used in thermosoil_readjust
133	!$OMP THREADPRIVATE(ptn_beg)
134	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: temp_sol_beg !! Surface temperature at previous timestep (K)
135	!$OMP THREADPRIVATE(temp_sol_beg)
136	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: surfheat_incr !! Change in soil heat content during the timestep
137	!! @tex ($J$) @endtex.
138	!$OMP THREADPRIVATE(surfheat_incr)
139	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: coldcont_incr !! Change in snow heat content @tex ($J$) @endtex.
140	!$OMP THREADPRIVATE(coldcont_incr)
141	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: shum_ngrnd_perma !! Saturation degree on the thermal axes (0-1, dimensionless)
142	!$OMP THREADPRIVATE(shum_ngrnd_perma)
143
144	! Variables related to soil freezing
145	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: profil_froz !! Frozen fraction of the soil on hydrological levels (-)
146	!$OMP THREADPRIVATE(profil_froz)
147	REAL(r_std),ALLOCATABLE, SAVE, DIMENSION (:) :: e_soil_lat !! Accumulated latent heat for the whole soil (J)
148	!$OMP THREADPRIVATE(e_soil_lat)
149	REAL(r_std),ALLOCATABLE, SAVE, DIMENSION (:,:) :: pcappa_supp !! Additional heat capacity due to soil freezing for each soil layer (J/K)
150	!$OMP THREADPRIVATE(pcappa_supp)
151	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: dz5 !! Used for numerical calculation [-]
152	!$OMP THREADPRIVATE(dz5)
153	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: mcs !! Saturation humidity [m3/m3]
154	!$OMP THREADPRIVATE(mcs)
155	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: SMCMAX !! Soil porosity [m3/m3]
156	!$OMP THREADPRIVATE(SMCMAX)
157	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: QZ !! quartz content [-]
158	!$OMP THREADPRIVATE(QZ)
159	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: so_capa_dry_ns !! Dry soil Heat capacity of soils,J.m^{-3}.K^{-1}
160	!$OMP THREADPRIVATE(so_capa_dry_ns)
161	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: mc_layt !! Volumetric soil moisture (liquid+ice) (m3/m3) on the thermodynamical levels at interface
162	!$OMP THREADPRIVATE(mc_layt)
163	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: mcl_layt !! Volumetric soil moisture (liquid) (m3/m3) on the thermodynamical levels at interface
164	!$OMP THREADPRIVATE(mcl_layt)
165	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: tmc_layt !! Total soil moisture content for each layer (liquid+ice) (mm) on the thermodynamical levels
166	!$OMP THREADPRIVATE(tmc_layt)
167	INTEGER(i_std), SAVE :: brk_flag = 0 !! Flag to consider bedrock: 0.no; 1.yes
168	!$OMP THREADPRIVATE(brk_flag)
169
170
171	CONTAINS
172	!! =============================================================================================================================
173	!! SUBROUTINE : thermosoil_initialize
174	!!
175	!>\BRIEF Allocate module variables, read from restart file or initialize with default values
176	!!
177	!! DESCRIPTION : Allocate module variables, read from restart file or initialize with default values.
178	!! Call thermosoil_var_init to calculate physical constants.
179	!! Call thermosoil_coef to calculate thermal soil properties.
180	!!
181	!! RECENT CHANGE(S) : None
182	!!
183	!! REFERENCE(S) : None
184	!!
185	!! FLOWCHART : None
186	!! \n
187	!_ ==============================================================================================================================
188	SUBROUTINE thermosoil_initialize (kjit, kjpindex, rest_id, &
189	temp_sol_new, snow, shumdiag_perma, &
190	soilcap, soilflx, stempdiag, &
191	gtemp, &
192	mc_layh, mcl_layh, tmc_layh, njsc, &
193	frac_snow_veg,frac_snow_nobio,totfrac_nobio, &
194	snowdz, snowrho, snowtemp, lambda_snow, cgrnd_snow, dgrnd_snow, pb)
195
196	!! 0. Variable and parameter declaration
197	!! 0.1 Input variables
198	INTEGER(i_std), INTENT(in) :: kjit !! Time step number (unitless)
199	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (unitless)
200	INTEGER(i_std),INTENT (in) :: rest_id !! Restart file identifier (unitless)
201	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_sol_new !! Surface temperature at the present time-step,
202	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: snow !! Snow mass (kg)
203	REAL(r_std),DIMENSION (kjpindex,nslm), INTENT (in) :: shumdiag_perma !! Soil saturation degree (0-1, unitless)
204	REAL(r_std),DIMENSION (kjpindex,nslm), INTENT (in) :: mc_layh !! Volumetric soil moisture content (liquid+ice) for hydrological layers, at node (m3/m3)
205	REAL(r_std),DIMENSION (kjpindex,nslm), INTENT (in) :: mcl_layh !! Volumetric soil moisture content (liquid) for hydrological layers, at node (m3/m3)
206	REAL(r_std),DIMENSION (kjpindex,nslm), INTENT (in) :: tmc_layh !! Total soil moisture content(liquid+ice) for hydrological layers (mm)
207	INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: njsc !! Index of the dominant soil textural class in the grid cell (1-nscm, unitless)
208	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: frac_snow_veg !! Snow cover fraction on vegeted area
209	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT(in) :: frac_snow_nobio !! Snow cover fraction on non-vegeted area
210	REAL(r_std),DIMENSION (kjpindex),INTENT(in) :: totfrac_nobio !! Total fraction of continental ice+lakes+cities+...
211	!! (unitless,0-1)
212	REAL(r_std),DIMENSION (kjpindex,nsnow),INTENT (in) :: snowdz !! Snow depth
213	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(in) :: snowrho !! Snow density
214	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(in) :: snowtemp !! Snow temperature (K)
215	REAL(r_std), DIMENSION (kjpindex), INTENT (in) :: pb !! Surface presure (hPa)
216
217	!! 0.2 Output variables
218	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: soilcap !! apparent surface heat capacity considering snow and soil surface (J m-2 K-1)
219	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: soilflx !! apparent soil heat flux considering snow and soil surface (W m-2)
220	REAL(r_std),DIMENSION (kjpindex,nslm), INTENT (out) :: stempdiag !! temperature profile on the levels in hydrol(K)
221	REAL(r_std),DIMENSION (kjpindex),INTENT(out) :: gtemp !! First soil layer temperature
222
223	!! 0.3 Modified variables
224	REAL(r_std), DIMENSION (kjpindex), INTENT(inout) :: lambda_snow !! Coefficient of the linear extrapolation of surface temperature
225	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT (inout):: cgrnd_snow !! Integration coefficient for snow numerical scheme
226	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT (inout):: dgrnd_snow !! Integration coefficient for snow numerical scheme
227
228	!! 0.4 Local variables
229	INTEGER(i_std) :: ier, i, jg
230	LOGICAL :: calculate_coef !! Local flag to initialize variables by call to thermosoil_coef
231	!_ ================================================================================================================================
232
233
234	!
235	! !! Flag to consider bedrock at deeper layers
236	! !! It affects heat capacity and thermal conductivity (energy balance).
237	!
238	!Config Key = BEDROCK_FLAG
239	!Config Desc = Flag to consider bedrock at deeper layers.
240	!Config If =
241	!Config Def = 0
242	!Config Help = 0, no, 1, yes.
243	!Config Units = [FLAG]
244	brk_flag = 0
245	CALL getin_p('BEDROCK_FLAG', brk_flag)
246
247	IF (printlev >= 3) WRITE (numout,*) 'Start thermosoil_initialize '
248
249	!! 1. Allocate soil temperatures variables
250	ALLOCATE (ptn(kjpindex,ngrnd),stat=ier)
251	IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of ptn','','')
252
253	ALLOCATE (dz1(ngrnd),stat=ier)
254	IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of dz1','','')
255
256	ALLOCATE (cgrnd(kjpindex,ngrnd-1),stat=ier)
257	IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of cgrnd','','')
258
259	ALLOCATE (dgrnd(kjpindex,ngrnd-1),stat=ier)
260	IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of dgrnd','','')
261
262	ALLOCATE (pcapa(kjpindex,ngrnd),stat=ier)
263	IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of pcapa','','')
264
265	ALLOCATE (pkappa(kjpindex,ngrnd),stat=ier)
266	IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of pkappa','','')
267
268	ALLOCATE (pcapa_snow(kjpindex,nsnow),stat=ier)
269	IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of pcapa_snow','','')
270
271	ALLOCATE (pkappa_snow(kjpindex,nsnow),stat=ier)
272	IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of pkappa_snow','','')
273
274	! Temporary fix: Initialize following variable because they are output to xios before the first calculation
275	pcapa = 0
276	pkappa = 0
277	pcapa_snow = 0
278	pkappa_snow = 0
279
280	ALLOCATE (surfheat_incr(kjpindex),stat=ier)
281	IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of surfheat_incr','','')
282
283	ALLOCATE (coldcont_incr(kjpindex),stat=ier)
284	IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of coldcont_incr','','')
285
286	ALLOCATE (pcapa_en(kjpindex,ngrnd),stat=ier)
287	IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of pcapa_en','','')
288
289	ALLOCATE (ptn_beg(kjpindex,ngrnd),stat=ier)
290	IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of ptn_beg','','')
291
292	ALLOCATE (temp_sol_beg(kjpindex),stat=ier)
293	IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of temp_sol_beg','','')
294
295	ALLOCATE (shum_ngrnd_perma(kjpindex,ngrnd),stat=ier)
296	IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of shum_ngrnd_perma','','')
297
298	ALLOCATE (profil_froz(kjpindex,ngrnd),stat=ier)
299	IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of profil_froz','','')
300
301	IF (ok_freeze_thermix) THEN
302	ALLOCATE (pcappa_supp(kjpindex,ngrnd),stat=ier)
303	IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of ok_freeze_termix','','')
304	END IF
305	IF (ok_Ecorr) THEN
306	ALLOCATE (e_soil_lat(kjpindex),stat=ier)
307	IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of e_soil_lat','','')
308	END IF
309
310	ALLOCATE (dz5(ngrnd),stat=ier)
311	IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of dz5','','')
312
313	ALLOCATE (mc_layt(kjpindex,ngrnd),stat=ier)
314	IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of mc_layt','','')
315
316	ALLOCATE (mcl_layt(kjpindex,ngrnd),stat=ier)
317	IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of mcl_layt','','')
318
319	ALLOCATE (tmc_layt(kjpindex,ngrnd),stat=ier)
320	IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of tmc_layt','','')
321
322	ALLOCATE (mcs(nscm),stat=ier)
323	IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of mcs','','')
324
325	ALLOCATE (SMCMAX(nscm),stat=ier)
326	IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of SMCMAX','','')
327
328	ALLOCATE (QZ(nscm),stat=ier)
329	IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of QZ','','')
330
331	ALLOCATE (so_capa_dry_ns(nscm),stat=ier)
332	IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of so_capa_dry_ns','','')
333
334
335	!! Soil texture choose
336	SELECTCASE (nscm)
337	CASE (3)
338	SMCMAX(:) = SMCMAX_fao(:)
339	QZ(:) = QZ_fao(:)
340	so_capa_dry_ns(:) = so_capa_dry_ns_fao(:)
341	mcs(:) = mcs_fao(:)
342	CASE (12)
343	SMCMAX(:) = SMCMAX_usda(:)
344	QZ(:) = QZ_usda(:)
345	so_capa_dry_ns(:) = so_capa_dry_ns_usda(:)
346	mcs(:) = mcs_usda(:)
347	CASE DEFAULT
348	WRITE (numout,*) 'Unsupported soil type classification. Choose between zobler, fao and usda according to the map'
349	STOP 'thermosoil_initialize'
350	ENDSELECT
351
352
353	!! 2. Initialize variable from restart file or with default values
354
355	!! Reads restart files for soil temperatures only. If no restart file is
356	!! found, the initial soil temperature is by default set to 280K at all depths. The user
357	!! can decide to initialize soil temperatures at an other value, in which case he should set the flag THERMOSOIL_TPRO
358	!! to this specific value in the run.def.
359	IF (printlev>=3) WRITE (numout,*) 'Read restart file for THERMOSOIL variables'
360
361	CALL ioconf_setatt_p('UNITS', 'K')
362	CALL ioconf_setatt_p('LONG_NAME','Soil Temperature profile')
363	CALL restget_p (rest_id, 'ptn', nbp_glo, ngrnd, 1, kjit, .TRUE., ptn, "gather", nbp_glo, index_g)
364
365	! Initialize ptn if it was not found in restart file
366	IF (ALL(ptn(:,:)==val_exp)) THEN
367	! ptn was not found in restart file
368
369	IF (read_reftemp) THEN
370	! Read variable ptn from file
371	CALL thermosoil_read_reftempfile(kjpindex,lalo,ptn)
372	ELSE
373	! Initialize ptn with a constant value which can be set in run.def
374
375	!Config Key = THERMOSOIL_TPRO
376	!Config Desc = Initial soil temperature profile if not found in restart
377	!Config Def = 280.
378	!Config If = OK_SECHIBA
379	!Config Help = The initial value of the temperature profile in the soil if
380	!Config its value is not found in the restart file. Here
381	!Config we only require one value as we will assume a constant
382	!Config throughout the column.
383	!Config Units = Kelvin [K]
384	CALL setvar_p (ptn, val_exp,'THERMOSOIL_TPRO',280._r_std)
385	END IF
386	END IF
387
388	! Initialize ptn_beg (variable needed in thermosoil_readadjust called from thermosoil_coef)
389	ptn_beg(:,:) = ptn(:,:)
390
391	! Initialize temp_sol_beg with values from previous time-step
392	temp_sol_beg(:) = temp_sol_new(:)
393
394	! Read e_soil_lat from restart file or initialize
395	IF (ok_Ecorr) THEN
396	CALL restget_p (rest_id, 'e_soil_lat', nbp_glo, 1, 1, kjit, .TRUE., &
397	e_soil_lat, "gather", nbp_glo, index_g)
398	CALL setvar_p (e_soil_lat, val_exp,'NO_KEYWORD',zero)
399	END IF
400
401	! Read gtemp from restart file
402	CALL restget_p (rest_id, 'gtemp', nbp_glo, 1, 1, kjit, .TRUE., &
403	gtemp, "gather", nbp_glo, index_g)
404	CALL setvar_p (gtemp, val_exp,'NO_KEYWORD',zero)
405
406
407	! Read variables calculated in thermosoil_coef from restart file
408	! If the variables were not found in the restart file, the logical
409	! calculate_coef will be true and thermosoil_coef will be called further below.
410	! These variables need to be in the restart file to avoid a time shift that
411	! would be done using thermosoil_coef at this stage.
412	calculate_coef=.FALSE.
413	CALL ioconf_setatt_p('UNITS', 'J m-2 K-1')
414	CALL ioconf_setatt_p('LONG_NAME','Apparent surface heat capacity')
415	CALL restget_p (rest_id, 'soilcap', nbp_glo, 1, 1, kjit, .TRUE., &
416	soilcap, "gather", nbp_glo, index_g)
417	IF (ALL(soilcap(:)==val_exp)) calculate_coef=.TRUE.
418
419	CALL ioconf_setatt_p('UNITS', 'W m-2')
420	CALL ioconf_setatt_p('LONG_NAME','Apparent soil heat flux')
421	CALL restget_p (rest_id, 'soilflx', nbp_glo, 1, 1, kjit, .TRUE., &
422	soilflx, "gather", nbp_glo, index_g)
423	IF (ALL(soilflx(:)==val_exp)) calculate_coef=.TRUE.
424
425	CALL ioconf_setatt_p('UNITS', '')
426	CALL ioconf_setatt_p('LONG_NAME','Integration coefficient for the numerical scheme')
427	CALL restget_p (rest_id, 'cgrnd', nbp_glo, ngrnd-1, 1, kjit, .TRUE., &
428	cgrnd, "gather", nbp_glo, index_g)
429	IF (ALL(cgrnd(:,:)==val_exp)) calculate_coef=.TRUE.
430
431	CALL ioconf_setatt_p('UNITS', '')
432	CALL ioconf_setatt_p('LONG_NAME','Integration coefficient for the numerical scheme')
433	CALL restget_p (rest_id, 'dgrnd', nbp_glo, ngrnd-1, 1, kjit, .TRUE., &
434	dgrnd, "gather", nbp_glo, index_g)
435	IF (ALL(dgrnd(:,:)==val_exp)) calculate_coef=.TRUE.
436
437	CALL ioconf_setatt_p('UNITS', '')
438	CALL ioconf_setatt_p('LONG_NAME','Integration coefficient for the numerical scheme')
439	CALL restget_p (rest_id, 'cgrnd_snow', nbp_glo, nsnow, 1, kjit, .TRUE., &
440	cgrnd_snow, "gather", nbp_glo, index_g)
441	IF (ALL(cgrnd_snow(:,:)==val_exp)) calculate_coef=.TRUE.
442
443	CALL ioconf_setatt_p('UNITS', '')
444	CALL ioconf_setatt_p('LONG_NAME','Integration coefficient for the numerical scheme')
445	CALL restget_p (rest_id, 'dgrnd_snow', nbp_glo, nsnow, 1, kjit, .TRUE., &
446	dgrnd_snow, "gather", nbp_glo, index_g)
447	IF (ALL(dgrnd_snow(:,:)==val_exp)) calculate_coef=.TRUE.
448
449	CALL ioconf_setatt_p('UNITS', '')
450	CALL ioconf_setatt_p('LONG_NAME','Coefficient of the linear extrapolation of surface temperature')
451	CALL restget_p (rest_id, 'lambda_snow', nbp_glo, 1, 1, kjit, .TRUE., &
452	lambda_snow, "gather", nbp_glo, index_g)
453	IF (ALL(lambda_snow(:)==val_exp)) calculate_coef=.TRUE.
454
455	!! 2.2 Computes some physical constants and arrays depending on the soil vertical discretization
456
457	! Calculate so_capa_ice
458	so_capa_ice = so_capa_dry + poroscapa_icerho_ice
459	IF (printlev>=2) WRITE(numout,*) 'Calculation of so_capa_ice=', so_capa_ice,' using poros=',poros,' and capa_ice=',capa_ice
460
461	! Computing some usefull constants for the numerical scheme
462	! Use znt(depth of nodes) and zlt(depth of deeper layer interface) from vertical_soil module.
463	DO jg=1,ngrnd-1
464	dz1(jg) = un / (znt(jg+1) - znt(jg))
465	dz5(jg) = (zlt(jg) - znt(jg)) * dz1(jg)
466	ENDDO
467	dz5(ngrnd) = 0.0
468	lambda = znt(1) * dz1(1)
469
470	! Send out the temperature profile on the first nslm levels(the levels treated in hydrol)
471	stempdiag(:,:) = ptn(:,1:nslm)
472
473
474	!! 2.3. Computes cgrnd, dgrnd, soilflx and soilcap coefficients only if they were not found in restart file.
475	IF (calculate_coef) THEN
476	! Interpolate variables needed by thermosoil_coef to the thermal levels
477	CALL thermosoil_humlev(kjpindex, shumdiag_perma, mc_layh, mcl_layh, tmc_layh)
478
479	IF (printlev>=3) WRITE (numout,*) 'thermosoil_coef will be called in the intialization phase'
480	CALL thermosoil_coef (&
481	kjpindex, temp_sol_new, snow, njsc, &
482	frac_snow_veg, frac_snow_nobio, totfrac_nobio, &
483	snowdz, snowrho, snowtemp, pb, &
484	ptn, &
485	soilcap, soilflx, cgrnd, dgrnd,&
486	lambda_snow, cgrnd_snow, dgrnd_snow)
487	END IF
488
489	END SUBROUTINE thermosoil_initialize
490
491
492	!! ================================================================================================================================
493	!! SUBROUTINE : thermosoil_main
494	!!
495	!>\BRIEF Thermosoil_main computes the soil thermal properties and dynamics, ie solves
496	!! the heat diffusion equation within the soil.
497	!!
498	!! DESCRIPTION : The resolution of the soil heat diffusion equation
499	!! relies on a numerical finite-difference implicit scheme
500	!! fully described in the reference and in the header of the thermosoil module.
501	!! - The dependency of the previous timestep hidden in the
502	!! integration coefficients cgrnd and dgrnd (EQ1), calculated in thermosoil_coef, and
503	!! called at the end of the routine to prepare for the next timestep.
504	!! - The effective computation of the new soil temperatures is performed in thermosoil_profile.
505	!!
506	!! - thermosoil_coef calculates the coefficients for the numerical scheme for the very first iteration of thermosoil;
507	!! after that, thermosoil_coef is called only at the end of the module to calculate the coefficients for the next timestep.
508	!! - thermosoil_profile solves the numerical scheme.\n
509	!!
510	!! - Flags : one unique flag : THERMOSOIL_TPRO (to be set to the desired initial soil in-depth temperature in K; by default 280K)
511	!!
512	!! RECENT CHANGE(S) : Change vertical discretization (consistent with hydrology layers) and soil thermal properties (taking into account soil texture effects).
513	!!
514	!! MAIN OUTPUT VARIABLE(S): vertically discretized soil temperatures ptn, soil
515	!! thermal properties (pcapa, pkappa), apparent surface heat capacity (soilcap)
516	!! and heat flux (soilflx) to be used in enerbil at the next timestep to solve
517	!! the surface energy balance.
518	!!
519	!! REFERENCE(S) :
520	!! - Hourdin, F. (1992). Study and numerical simulation of the general circulation of planetary atmospheres,
521	!! Ph.D. thesis, Paris VII University. Remark: the part of F. Hourdin' s PhD thesis relative to the thermal
522	!! integration scheme has been scanned and is provided along with the documentation, with name :
523	!! Hourdin_1992_PhD_thermal_scheme.pdf
524	!!
525	!! FLOWCHART :
526	!! \latexonly
527	!! \includegraphics[scale = 1]{thermosoil_flowchart.png}
528	!! \endlatexonly
529	!!
530	!! \n
531	!_ ================================================================================================================================
532
533	SUBROUTINE thermosoil_main (kjit, kjpindex, &
534	index, indexgrnd, &
535	temp_sol_new, snow, soilcap, soilflx, &
536	shumdiag_perma, stempdiag, ptnlev1, rest_id, hist_id, hist2_id, &
537	snowdz,snowrho,snowtemp,gtemp,pb,&
538	mc_layh, mcl_layh, tmc_layh, njsc, frac_snow_veg,frac_snow_nobio,totfrac_nobio,temp_sol_add, &
539	lambda_snow, cgrnd_snow, dgrnd_snow)
540
541	!! 0. Variable and parameter declaration
542
543	!! 0.1 Input variables
544
545	INTEGER(i_std), INTENT(in) :: kjit !! Time step number (unitless)
546	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (unitless)
547	INTEGER(i_std),INTENT (in) :: rest_id,hist_id !! Restart_ file and history file identifier
548	!! (unitless)
549	INTEGER(i_std),INTENT (in) :: hist2_id !! history file 2 identifier (unitless)
550	INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: index !! Indeces of the points on the map (unitless)
551	INTEGER(i_std),DIMENSION (kjpindex*ngrnd), INTENT (in):: indexgrnd !! Indeces of the points on the 3D map (vertical
552	!! dimension towards the ground) (unitless)
553	REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: temp_sol_new !! Surface temperature at the present time-step,
554	!! temp_sol_new is only modified for the case ok_explicitsnow
555	!! Ts @tex ($K$) @endtex
556	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: snow !! Snow mass @tex ($kg$) @endtex.
557	!! Caveat: when there is snow on the
558	!! ground, the snow is integrated into the soil for
559	!! the calculation of the thermal dynamics. It means
560	!! that the uppermost soil layers can completely or
561	!! partially consist in snow. In the second case, zx1
562	!! and zx2 are the fraction of the soil layer
563	!! consisting in snow and 'normal' soil, respectively
564	!! This is calculated in thermosoil_coef.
565	REAL(r_std),DIMENSION (kjpindex,nslm), INTENT (in) :: shumdiag_perma !! Soil saturation degree (0-1, unitless)
566	REAL(r_std),DIMENSION (kjpindex,nsnow),INTENT (in) :: snowdz !! Snow depth
567	REAL(r_std),DIMENSION (kjpindex,nsnow),INTENT (in) :: snowrho !! Snow density
568	REAL(r_std),DIMENSION (kjpindex,nsnow),INTENT (inout) :: snowtemp !! Snow temperature (K)
569	REAL(r_std), DIMENSION (kjpindex),INTENT (in) :: pb !! Surface presure (hPa)
570	REAL(r_std),DIMENSION (kjpindex,nslm), INTENT (in) :: mc_layh !! Volumetric soil moisture content for each layer in hydrol at nodes(liquid + ice) (m3/m3)
571	REAL(r_std),DIMENSION (kjpindex,nslm), INTENT (in) :: mcl_layh !! Volumetric soil moisture content for each layer in hydrol at nodes(liquid) (m3/m3)
572	REAL(r_std),DIMENSION (kjpindex,nslm), INTENT (in) :: tmc_layh !! Total soil moisture content for each layer in hydrol(liquid + ice) (mm)
573	INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: njsc !! Index of the dominant soil textural class in the grid cell (1-nscm, unitless)
574	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: frac_snow_veg !! Snow cover fraction on vegeted area
575	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT(in) :: frac_snow_nobio !! Snow cover fraction on non-vegeted area
576	REAL(r_std),DIMENSION (kjpindex),INTENT(in) :: totfrac_nobio !! Total fraction of continental ice+lakes+cities+...
577	!!(unitless,0-1)
578	REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: temp_sol_add !! additional surface temperature due to the melt of first layer
579	!! at the present time-step @tex ($K$) @endtex
580
581	!! 0.2 Output variables
582
583	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: ptnlev1 !! 1st level soil temperature
584	REAL(r_std),DIMENSION (kjpindex),INTENT(out) :: gtemp !! First soil layer temperature
585
586
587	!! 0.3 Modified variables
588
589	REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: soilcap !! apparent surface heat capacity considering snow and soil surface
590	!! @tex ($J m^{-2} K^{-1}$) @endtex
591	REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: soilflx !! apparent soil heat flux considering snow and soil surface
592	!! @tex ($W m^{-2}$) @endtex
593	!! , positive
594	!! towards the soil, writen as Qg (ground heat flux)
595	!! in the history files, and computed at the end of
596	!! thermosoil for the calculation of Ts in enerbil,
597	!! see EQ3.
598	REAL(r_std),DIMENSION (kjpindex,nslm), INTENT (out) :: stempdiag !! temperature profile @tex ($K$) @endtex
599	REAL(r_std),DIMENSION (kjpindex), INTENT(inout) :: lambda_snow !! Coefficient of the linear extrapolation of surface temperature
600	REAL(r_std),DIMENSION (kjpindex,nsnow), INTENT (inout):: cgrnd_snow !! Integration coefficient for snow numerical scheme
601	REAL(r_std),DIMENSION (kjpindex,nsnow), INTENT (inout):: dgrnd_snow !! Integration coefficient for snow numerical scheme
602
603	!! 0.4 Local variables
604
605	INTEGER(i_std) :: jv,ji,ii
606	REAL(r_std),DIMENSION (kjpindex) :: snowtemp_weighted!! Snow temperature weighted by snow density, only for diag (K)
607	REAL(r_std),DIMENSION (kjpindex, nsnow) :: pkappa_snow_diag !! Only for diag, containing xios_default_val
608	REAL(r_std),DIMENSION (kjpindex, nsnow) :: pcapa_snow_diag !! Only for diag, containing xios_default_val
609	REAL(r_std),DIMENSION (kjpindex, nsnow) :: snowtemp_diag !! Only for diag, containing xios_default_val
610
611	!_ ================================================================================================================================
612
613	!! 3. Put the soil wetness diagnostic on the levels of the soil temperature
614
615	!!?? this could logically be put just before the last call to
616	!!thermosoil_coef, as the results are used there...
617	CALL thermosoil_humlev(kjpindex, shumdiag_perma, mc_layh, mcl_layh, tmc_layh)
618
619
620	!! 4. Effective computation of the soil temperatures profile.
621	!! cgrnd and dgrnd have been calculated in thermosoil_coef at the previous time step
622	!! but they are correct for the actual time-step.
623	CALL thermosoil_profile (kjpindex, temp_sol_new, &
624	frac_snow_veg, frac_snow_nobio, totfrac_nobio, &
625	ptn, stempdiag, snowtemp, &
626	cgrnd_snow, dgrnd_snow)
627
628
629	!! 5. Call to thermosoil_energy_diag for calculation of diagnostic variables
630	CALL thermosoil_energy_diag(kjpindex, temp_sol_new, soilcap)
631
632	!! Save ptn at current stage, to be used in thermosoil_readjust
633	ptn_beg(:,:) = ptn(:,:)
634
635	!! 6. Writing the history files according to the ALMA standards (or not..)
636
637	! Add XIOS default value where no snow
638	DO ji=1,kjpindex
639	IF (snow(ji) .GT. zero) THEN
640	pkappa_snow_diag(ji,:) = pkappa_snow(ji,:)
641	pcapa_snow_diag(ji,:) = pcapa_snow(ji,:)
642	snowtemp_diag(ji,:) = snowtemp(ji,:)
643	ELSE
644	pkappa_snow_diag(ji,:) = xios_default_val
645	pcapa_snow_diag(ji,:) = xios_default_val
646	snowtemp_diag(ji,:) = xios_default_val
647	END IF
648	END DO
649
650	IF (ok_explicitsnow) THEN
651	DO ji=1,kjpindex
652	! Use min_sechiba instead of zero to avoid problem with division by zero
653	IF (snow(ji) .GT. min_sechiba) THEN
654	snowtemp_weighted(ji) = SUM(snowtemp(ji,:)*snowrho(ji,:))/SUM(snowrho(ji,:))
655	ELSE
656	snowtemp_weighted(ji) = xios_default_val
657	END IF
658	END DO
659	CALL xios_orchidee_send_field("snowtemp_weighted",snowtemp_weighted)
660	END IF
661
662	CALL xios_orchidee_send_field("ptn",ptn)
663	CALL xios_orchidee_send_field("soilflx",soilflx)
664	CALL xios_orchidee_send_field("surfheat_incr",surfheat_incr)
665	CALL xios_orchidee_send_field("coldcont_incr",coldcont_incr)
666	CALL xios_orchidee_send_field("pkappa",pkappa)
667	CALL xios_orchidee_send_field("pkappa_snow",pkappa_snow_diag)
668	CALL xios_orchidee_send_field("pcapa",pcapa)
669	CALL xios_orchidee_send_field("pcapa_snow",pcapa_snow_diag)
670	CALL xios_orchidee_send_field("snowtemp",snowtemp_diag)
671
672
673	IF ( .NOT. almaoutput ) THEN
674	CALL histwrite_p(hist_id, 'ptn', kjit, ptn, kjpindex*ngrnd, indexgrnd)
675	CALL histwrite_p(hist_id, 'Qg', kjit, soilflx, kjpindex, index)
676
677	IF ( hydrol_cwrr ) THEN
678	CALL histwrite_p(hist_id, 'ptn_beg', kjit, ptn_beg, kjpindex*ngrnd, indexgrnd)
679	CALL histwrite_p(hist_id, 'pkappa', kjit, pkappa, kjpindex*ngrnd, indexgrnd)
680	CALL histwrite_p(hist_id, 'pcapa', kjit, pcapa, kjpindex*ngrnd, indexgrnd)
681
682	IF (ok_freeze_thermix) THEN
683	CALL histwrite_p(hist_id, 'profil_froz', kjit, profil_froz, kjpindex*ngrnd, indexgrnd)
684	CALL histwrite_p(hist_id, 'pcappa_supp', kjit, pcappa_supp, kjpindex*ngrnd, indexgrnd)
685	END IF
686	CALL histwrite_p(hist_id, 'shum_ngrnd_perma', kjit, shum_ngrnd_perma(:,:), kjpindex*ngrnd, indexgrnd)
687	END IF
688
689	ELSE
690	CALL histwrite_p(hist_id, 'SoilTemp', kjit, ptn, kjpindex*ngrnd, indexgrnd)
691	CALL histwrite_p(hist_id, 'Qg', kjit, soilflx, kjpindex, index)
692	CALL histwrite_p(hist_id, 'DelSurfHeat', kjit, surfheat_incr, kjpindex, index)
693	CALL histwrite_p(hist_id, 'DelColdCont', kjit, coldcont_incr, kjpindex, index)
694	ENDIF
695	IF ( hist2_id > 0 ) THEN
696	IF ( .NOT. almaoutput ) THEN
697	CALL histwrite_p(hist2_id, 'ptn', kjit, ptn, kjpindex*ngrnd, indexgrnd)
698	ELSE
699	CALL histwrite_p(hist2_id, 'SoilTemp', kjit, ptn, kjpindex*ngrnd, indexgrnd)
700	CALL histwrite_p(hist2_id, 'Qg', kjit, soilflx, kjpindex, index)
701	CALL histwrite_p(hist2_id, 'DelSurfHeat', kjit, surfheat_incr, kjpindex, index)
702	CALL histwrite_p(hist2_id, 'DelColdCont', kjit, coldcont_incr, kjpindex, index)
703	ENDIF
704	ENDIF
705
706	!! 7. A last final call to thermosoil_coef
707
708	!! A last final call to thermosoil_coef, which calculates the different
709	!!coefficients (cgrnd, dgrnd, soilcap, soilflx) from this time step to be
710	!!used at the next time step, either in the surface temperature calculation
711	!!(soilcap, soilflx) or in the soil thermal numerical scheme.
712	CALL thermosoil_coef (&
713	kjpindex, temp_sol_new, snow, njsc, &
714	frac_snow_veg, frac_snow_nobio, totfrac_nobio, &
715	snowdz, snowrho, snowtemp, pb, &
716	ptn, &
717	soilcap, soilflx, cgrnd, dgrnd,&
718	lambda_snow, cgrnd_snow, dgrnd_snow)
719
720
721	! Save variables for explicit snow model
722	gtemp(:) = ptn(:,1)
723
724	!! Initialize output arguments to be used in sechiba
725	ptnlev1(:) = ptn(:,1)
726
727	!! Surface temperature is forced to zero celcius if its value is larger than melting point, only for explicit snow scheme
728	IF (ok_explicitsnow) THEN
729	DO ji=1,kjpindex
730	IF (SUM(snowdz(ji,:)) .GT. 0.0) THEN
731	IF (temp_sol_new(ji) .GE. tp_00) THEN
732	temp_sol_new(ji) = tp_00
733	ENDIF
734	END IF
735	END DO
736	ENDIF
737
738	IF (printlev>=3) WRITE (numout,*) ' thermosoil_main done '
739
740	END SUBROUTINE thermosoil_main
741
742	!! =============================================================================================================================
743	!! SUBROUTINE : thermosoil_finalize
744	!!
745	!>\BRIEF Write to restart file
746	!!
747	!! DESCRIPTION : This subroutine writes the module variables and variables calculated in thermosoil
748	!! to restart file
749	!! \n
750	!_ ==============================================================================================================================
751	SUBROUTINE thermosoil_finalize (kjit, kjpindex, rest_id, gtemp, &
752	soilcap, soilflx, lambda_snow, cgrnd_snow, dgrnd_snow)
753
754	!! 0. Variable and parameter declaration
755	!! 0.1 Input variables
756	INTEGER(i_std), INTENT(in) :: kjit !! Time step number (unitless)
757	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (unitless)
758	INTEGER(i_std),INTENT (in) :: rest_id !! Restart file identifier(unitless)
759	REAL(r_std),DIMENSION (kjpindex),INTENT(in) :: gtemp !! First soil layer temperature
760	REAL(r_std),DIMENSION (kjpindex),INTENT(in) :: soilcap
761	REAL(r_std),DIMENSION (kjpindex),INTENT(in) :: soilflx
762	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: lambda_snow !! Coefficient of the linear extrapolation of surface temperature
763	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT (in) :: cgrnd_snow !! Integration coefficient for snow numerical scheme
764	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT (in) :: dgrnd_snow !! Integration coefficient for snow numerical scheme
765
766	!_ ================================================================================================================================
767
768	!! 1. Write variables to restart file to be used for the next simulation
769	IF (printlev>=3) WRITE (numout,*) 'Write restart file with THERMOSOIL variables'
770
771	CALL restput_p(rest_id, 'ptn', nbp_glo, ngrnd, 1, kjit, ptn, 'scatter', nbp_glo, index_g)
772
773	IF (ok_Ecorr) THEN
774	CALL restput_p(rest_id, 'e_soil_lat', nbp_glo, 1 , 1, kjit, e_soil_lat, 'scatter', nbp_glo, index_g)
775	END IF
776
777	CALL restput_p(rest_id, 'gtemp', nbp_glo, 1, 1, kjit, gtemp, 'scatter', nbp_glo, index_g)
778	CALL restput_p(rest_id, 'soilcap', nbp_glo, 1, 1, kjit, soilcap, 'scatter', nbp_glo, index_g)
779	CALL restput_p(rest_id, 'soilflx', nbp_glo, 1, 1, kjit, soilflx, 'scatter', nbp_glo, index_g)
780	CALL restput_p(rest_id, 'cgrnd', nbp_glo, ngrnd-1, 1, kjit, cgrnd, 'scatter', nbp_glo, index_g)
781	CALL restput_p(rest_id, 'dgrnd', nbp_glo, ngrnd-1, 1, kjit, dgrnd, 'scatter', nbp_glo, index_g)
782	CALL restput_p(rest_id, 'cgrnd_snow', nbp_glo, nsnow, 1, kjit, cgrnd_snow, 'scatter', nbp_glo, index_g)
783	CALL restput_p(rest_id, 'dgrnd_snow', nbp_glo, nsnow, 1, kjit, dgrnd_snow, 'scatter', nbp_glo, index_g)
784	CALL restput_p(rest_id, 'lambda_snow', nbp_glo, 1, 1, kjit, lambda_snow, 'scatter', nbp_glo, index_g)
785
786	END SUBROUTINE thermosoil_finalize
787
788
789	!! ================================================================================================================================
790	!! SUBROUTINE : thermosoil_clear
791	!!
792	!>\BRIEF Deallocates the allocated arrays.
793	!! The call of thermosoil_clear originates from sechiba_clear but the calling sequence and
794	!! its purpose require further investigation.
795	!!
796	!! DESCRIPTION : None
797	!!
798	!! RECENT CHANGE(S) : None
799	!!
800	!! MAIN OUTPUT VARIABLE(S): None
801	!!
802	!! REFERENCE(S) : None
803	!!
804	!! FLOWCHART : None
805	!! \n
806	!_ ================================================================================================================================
807
808	SUBROUTINE thermosoil_clear()
809
810	IF ( ALLOCATED (ptn)) DEALLOCATE (ptn)
811	IF ( ALLOCATED (cgrnd)) DEALLOCATE (cgrnd)
812	IF ( ALLOCATED (dgrnd)) DEALLOCATE (dgrnd)
813	IF ( ALLOCATED (pcapa)) DEALLOCATE (pcapa)
814	IF ( ALLOCATED (pkappa)) DEALLOCATE (pkappa)
815	IF ( ALLOCATED (pcapa_snow)) DEALLOCATE (pcapa_snow)
816	IF ( ALLOCATED (pkappa_snow)) DEALLOCATE (pkappa_snow)
817	IF ( ALLOCATED (pcapa_en)) DEALLOCATE (pcapa_en)
818	IF ( ALLOCATED (ptn_beg)) DEALLOCATE (ptn_beg)
819	IF ( ALLOCATED (temp_sol_beg)) DEALLOCATE (temp_sol_beg)
820	IF ( ALLOCATED (surfheat_incr)) DEALLOCATE (surfheat_incr)
821	IF ( ALLOCATED (coldcont_incr)) DEALLOCATE (coldcont_incr)
822	IF ( ALLOCATED (shum_ngrnd_perma)) DEALLOCATE (shum_ngrnd_perma)
823	IF ( ALLOCATED (profil_froz)) DEALLOCATE (profil_froz)
824	IF ( ALLOCATED (mc_layt)) DEALLOCATE (mc_layt)
825	IF ( ALLOCATED (mcl_layt)) DEALLOCATE (mcl_layt)
826	IF ( ALLOCATED (tmc_layt)) DEALLOCATE (tmc_layt)
827	END SUBROUTINE thermosoil_clear
828
829
830	!! ================================================================================================================================
831	!! SUBROUTINE : thermosoil_coef
832	!!
833	!>\BRIEF Calculate soil thermal properties, integration coefficients, apparent heat flux,
834	!! surface heat capacity,
835	!!
836	!! DESCRIPTION : This routine computes : \n
837	!! 1. the soil thermal properties. \n
838	!! 2. the integration coefficients of the thermal numerical scheme, cgrnd and dgrnd,
839	!! which depend on the vertical grid and on soil properties, and are used at the next
840	!! timestep.\n
841	!! 3. the soil apparent heat flux and surface heat capacity (soilflx
842	!! and soilcap), used by enerbil to compute the surface temperature at the next
843	!! timestep.\n
844	!! - The soil thermal properties depend on water content (shum_ngrnd_perma, shumdiag_perma,
845	!! mc_layt, mcl_layt, tmc_layt), dominant soil texture(njsc), and on the presence
846	!! of snow : snow is integrated into the soil for the thermal calculations, ie if there
847	!! is snow on the ground, the first thermal layer(s) consist in snow, depending on the
848	!! snow-depth. If a layer consists out of snow and soil, wheighed soil properties are
849	!! calculated\n
850	!! - The coefficients cgrnd and dgrnd are the integration
851	!! coefficients for the thermal scheme \n
852	!! T(k+1)=cgrnd(k)+dgrnd(k)*T(k) \n
853	!! -- EQ1 -- \n
854	!! They correspond respectively to $\beta$ and $\alpha$ from F. Hourdin\'s thesis and
855	!! their expression can be found in this document (eq A19 and A20)
856	!! - soilcap and soilflx are the apparent surface heat capacity and flux
857	!! used in enerbil at the next timestep to solve the surface
858	!! balance for Ts (EQ3); they correspond to $C_s$ and $F_s$ in F.
859	!! Hourdin\'s PhD thesis and are expressed in eq. A30 and A31. \n
860	!! soilcap*(Ts(t)-Ts(t-1))/dt=soilflx+otherfluxes(Ts(t)) \n
861	!! -- EQ3 --\n
862	!!
863	!! RECENT CHANGE(S) : None
864	!!
865	!! MAIN OUTPUT VARIABLE(S): cgrnd, dgrnd, pcapa, pkappa, soilcap, soilflx
866	!!
867	!! REFERENCE(S) :
868	!! - Hourdin, F. (1992). Study and numerical simulation of the general circulation of planetary atmospheres,
869	!! Ph.D. thesis, Paris VII University. Remark: the part of F. Hourdin's PhD thesis relative to the thermal
870	!! integration scheme has been scanned and is provided along with the documentation, with name :
871	!! Hourdin_1992_PhD_thermal_scheme.pdf
872	!!
873	!! FLOWCHART : None
874	!! \n
875	!_ ================================================================================================================================
876
877	SUBROUTINE thermosoil_coef (kjpindex, temp_sol_new, snow, njsc, &
878	frac_snow_veg, frac_snow_nobio, totfrac_nobio, &
879	snowdz, snowrho, snowtemp, pb, &
880	ptn, &
881	soilcap, soilflx, cgrnd, dgrnd,&
882	lambda_snow, cgrnd_snow, dgrnd_snow)
883
884	!! 0. Variables and parameter declaration
885
886	!! 0.1 Input variables
887
888	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (unitless)
889	REAL(r_std), DIMENSION (kjpindex), INTENT (in) :: temp_sol_new !! soil surface temperature @tex ($K$) @endtex
890	REAL(r_std), DIMENSION (kjpindex), INTENT (in) :: snow !! snow mass @tex ($Kg$) @endtex
891	INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: njsc !! Index of the dominant soil textural class in the grid cell (1-nscm, unitless)
892	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: frac_snow_veg !! Snow cover fraction on vegeted area
893	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT(in) :: frac_snow_nobio !! Snow cover fraction on non-vegeted area
894	REAL(r_std),DIMENSION (kjpindex),INTENT(in) :: totfrac_nobio !! Total fraction of continental ice+lakes+cities+...
895	!!(unitless,0-1)
896	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(in) :: snowdz !! Snow depth (m)
897	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(in) :: snowrho !! Snow density
898	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(in) :: snowtemp !! Snow temperature (K)
899	REAL(r_std), DIMENSION (kjpindex), INTENT (in) :: pb !! Surface presure (hPa)
900
901	!! 0.2 Output variables
902
903	REAL(r_std), DIMENSION (kjpindex), INTENT (out) :: soilcap !! surface heat capacity considering snow and soil surface
904	!! @tex ($J m^{-2} K^{-1}$) @endtex
905	REAL(r_std), DIMENSION (kjpindex), INTENT (out) :: soilflx !! surface heat flux considering snow and soil surface @tex ($W m^{-2}$) @endtex,
906	!! positive towards the
907	!! soil, writen as Qg (ground heat flux) in the history
908	!! files.
909	REAL(r_std), DIMENSION (kjpindex,ngrnd-1), INTENT(out) :: cgrnd !! matrix coefficient for the computation of soil
910	!! temperatures (beta in F. Hourdin thesis)
911	REAL(r_std), DIMENSION (kjpindex,ngrnd-1), INTENT(out) :: dgrnd !! matrix coefficient for the computation of soil
912	!! temperatures (alpha in F. Hourdin thesis)
913
914
915	!! 0.3 Modified variable
916
917	REAL(r_std), DIMENSION (kjpindex,ngrnd), INTENT (inout):: ptn !! vertically discretized soil temperatures. ptn is only modified if ok_Ecorr.
918	REAL(r_std), DIMENSION (kjpindex), INTENT(inout) :: lambda_snow !! Coefficient of the linear extrapolation of surface temperature
919	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT (inout):: cgrnd_snow !! Integration coefficient for snow numerical scheme
920	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT (inout):: dgrnd_snow !! Integration coefficient for snow numerical scheme
921
922	!! 0.4 Local variables
923
924	INTEGER(i_std) :: ji, jg
925	REAL(r_std), DIMENSION (kjpindex,ngrnd-1) :: zdz1 !! numerical (buffer) constant
926	!! @tex ($W m^{-1} K^{-1}$) @endtex
927	REAL(r_std), DIMENSION (kjpindex,ngrnd) :: zdz2 !! numerical (buffer) constant
928	!! @tex ($W m^{-1} K^{-1}$) @endtex
929	REAL(r_std), DIMENSION (kjpindex) :: z1 !! numerical constant @tex ($W m^{-1} K^{-1}$) @endtex
930	REAL(r_std), DIMENSION (kjpindex) :: soilcap_nosnow !! surface heat capacity
931	!! @tex ($J m^{-2} K^{-1}$)
932	!! @endtex
933	REAL(r_std), DIMENSION (kjpindex) :: soilflx_nosnow !! surface heat flux @tex ($W m^{-2}$) @endtex,
934	!! positive towards the soil, written as Qg
935	!!(ground heat flux in the history files).
936	REAL(r_std), DIMENSION (kjpindex) :: snowcap !! apparent snow heat capacity @tex ($J m^{-2} K^{-1}$)
937	REAL(r_std), DIMENSION (kjpindex) :: snowflx !! apparent snow-atmosphere heat flux @tex ($W m^{-2}$) @endtex
938	REAL(r_std), DIMENSION (kjpindex,nsnow) :: dz1_snow
939	REAL(r_std), DIMENSION (kjpindex,nsnow) :: ZSNOWDZM
940	REAL(r_std), DIMENSION (kjpindex,nsnow) :: dz2_snow
941	REAL(r_std), DIMENSION (kjpindex,nsnow) :: zdz1_snow
942	REAL(r_std), DIMENSION (kjpindex,nsnow) :: zdz2_snow
943	REAL(r_std), DIMENSION (kjpindex) :: z1_snow
944	REAL(r_std), DIMENSION (kjpindex) :: snowflxtot !! Total snow flux (including snow on vegetated and bare soil and nobio areas)
945	!! @tex ($W m^{-2}$) @endtex
946	!! positive towards the soil
947
948	!_ ================================================================================================================================
949
950	!! 1. Computation of the soil thermal properties
951
952	! Computation of the soil thermal properties; snow properties are also accounted for
953	IF (ok_freeze_thermix) THEN
954	CALL thermosoil_getdiff( kjpindex, snow, ptn, njsc, snowrho, snowtemp, pb )
955	ELSE IF (ok_explicitsnow) THEN
956	! Special case with explicit snow model without soil freezing
957	CALL thermosoil_getdiff_old_thermix_without_snow( kjpindex, njsc, snowrho, snowtemp, pb )
958	ELSE
959	! Special case with old snow without soil freezing
960	CALL thermosoil_getdiff_old_thermix_with_snow( kjpindex, snow, njsc )
961	ENDIF
962
963	! Energy conservation : Correction to make sure that the same latent heat is released and
964	! consumed during freezing and thawing
965	IF (ok_Ecorr) THEN
966	CALL thermosoil_readjust(kjpindex, ptn)
967	ENDIF
968
969
970	!! 2. Computation of the coefficients of the numerical integration scheme for the soil layers
971
972	!! 2.1 Calculate numerical coefficients zdz1 and zdz2
973	DO jg=1,ngrnd
974	DO ji=1,kjpindex
975	zdz2(ji,jg)=pcapa(ji,jg) * dlt(jg)/dt_sechiba
976	ENDDO
977	ENDDO
978
979	DO jg=1,ngrnd-1
980	DO ji=1,kjpindex
981	zdz1(ji,jg) = dz1(jg) * pkappa(ji,jg)
982	ENDDO
983	ENDDO
984
985	!! 2.2 Calculate coefficients cgrnd and dgrnd for soil
986	DO ji = 1,kjpindex
987	z1(ji) = zdz2(ji,ngrnd) + zdz1(ji,ngrnd-1)
988	cgrnd(ji,ngrnd-1) = zdz2(ji,ngrnd) * ptn(ji,ngrnd) / z1(ji)
989	dgrnd(ji,ngrnd-1) = zdz1(ji,ngrnd-1) / z1(ji)
990	ENDDO
991
992	DO jg = ngrnd-1,2,-1
993	DO ji = 1,kjpindex
994	z1(ji) = un / (zdz2(ji,jg) + zdz1(ji,jg-1) + zdz1(ji,jg) * (un - dgrnd(ji,jg)))
995	cgrnd(ji,jg-1) = (ptn(ji,jg) * zdz2(ji,jg) + zdz1(ji,jg) * cgrnd(ji,jg)) * z1(ji)
996	dgrnd(ji,jg-1) = zdz1(ji,jg-1) * z1(ji)
997	ENDDO
998	ENDDO
999
1000
1001	!! 3. Computation of the coefficients of the numerical integration scheme for the snow layers
1002
1003	!! 3.1 Calculate numerical coefficients zdz1_snow, zdz2_snow and lambda_snow
1004	DO ji = 1, kjpindex
1005
1006	IF ( ok_explicitsnow ) THEN
1007
1008	! Calculate internal values
1009	DO jg = 1, nsnow
1010	ZSNOWDZM(ji,jg) = MAX(snowdz(ji,jg),psnowdzmin)
1011	ENDDO
1012	dz2_snow(ji,:)=ZSNOWDZM(ji,:)
1013
1014	DO jg = 1, nsnow-1
1015	dz1_snow(ji,jg) = 2.0 / (dz2_snow(ji,jg+1)+dz2_snow(ji,jg))
1016	ENDDO
1017
1018	lambda_snow(ji) = dz2_snow(ji,1)/2.0 * dz1_snow(ji,1)
1019
1020	DO jg=1,nsnow
1021	zdz2_snow(ji,jg)=pcapa_snow(ji,jg) * dz2_snow(ji,jg)/dt_sechiba
1022	ENDDO
1023
1024	DO jg=1,nsnow-1
1025	zdz1_snow(ji,jg) = dz1_snow(ji,jg) * pkappa_snow(ji,jg)
1026	ENDDO
1027
1028	! the bottom snow
1029	zdz1_snow(ji,nsnow) = pkappa_snow(ji,nsnow) / ( zlt(1) + dz2_snow(ji,nsnow)/2 )
1030
1031	ELSE
1032	! Without explict snow
1033	lambda_snow(ji) = lambda
1034	ENDIF
1035
1036	ENDDO
1037
1038	!! 3.2 Calculate coefficients cgrnd_snow and dgrnd_snow for snow
1039	DO ji = 1,kjpindex
1040	IF ( ok_explicitsnow ) THEN
1041	! bottom level
1042	z1_snow(ji) = zdz2(ji,1)+(un-dgrnd(ji,1))*zdz1(ji,1)+zdz1_snow(ji,nsnow)
1043	cgrnd_snow(ji,nsnow) = (zdz2(ji,1) * ptn(ji,1) + zdz1(ji,1) * cgrnd(ji,1) ) / z1_snow(ji)
1044	dgrnd_snow(ji,nsnow) = zdz1_snow(ji,nsnow) / z1_snow(ji)
1045
1046	! next-to-bottom level
1047	z1_snow(ji) = zdz2_snow(ji,nsnow)+(un-dgrnd_snow(ji,nsnow))*zdz1_snow(ji,nsnow)+zdz1_snow(ji,nsnow-1)
1048	cgrnd_snow(ji,nsnow-1) = (zdz2_snow(ji,nsnow)*snowtemp(ji,nsnow)+&
1049	zdz1_snow(ji,nsnow)*cgrnd_snow(ji,nsnow))/z1_snow(ji)
1050	dgrnd_snow(ji,nsnow-1) = zdz1_snow(ji,nsnow-1) / z1_snow(ji)
1051
1052	DO jg = nsnow-1,2,-1
1053	z1_snow(ji) = un / (zdz2_snow(ji,jg) + zdz1_snow(ji,jg-1) + zdz1_snow(ji,jg) * (un - dgrnd_snow(ji,jg)))
1054	cgrnd_snow(ji,jg-1) = (snowtemp(ji,jg) * zdz2_snow(ji,jg) + zdz1_snow(ji,jg) * cgrnd_snow(ji,jg)) * z1_snow(ji)
1055	dgrnd_snow(ji,jg-1) = zdz1_snow(ji,jg-1) * z1_snow(ji)
1056	ENDDO
1057	ELSE
1058	! Without explict snow
1059	cgrnd_snow(ji,:) = cgrnd(ji,1)
1060	dgrnd_snow(ji,:) = dgrnd(ji,1)
1061	ENDIF
1062	ENDDO
1063
1064
1065
1066	!! 4. Computation of the apparent ground heat flux
1067	!! Computation of apparent snow-atmosphere flux
1068	DO ji = 1,kjpindex
1069	IF ( ok_explicitsnow ) THEN
1070	snowflx(ji) = zdz1_snow(ji,1) * (cgrnd_snow(ji,1) + (dgrnd_snow(ji,1)-1.) * snowtemp(ji,1))
1071	snowcap(ji) = (zdz2_snow(ji,1) * dt_sechiba + dt_sechiba * (un - dgrnd_snow(ji,1)) * zdz1_snow(ji,1))
1072	z1_snow(ji) = lambda_snow(ji) * (un - dgrnd_snow(ji,1)) + un
1073	snowcap(ji) = snowcap(ji) / z1_snow(ji)
1074	snowflx(ji) = snowflx(ji) + &
1075	& snowcap(ji) * (snowtemp(ji,1) * z1_snow(ji) - lambda_snow(ji) * cgrnd_snow(ji,1) - temp_sol_new(ji)) / dt_sechiba
1076	ELSE
1077	snowflx(ji) = zero
1078	snowcap(ji) = zero
1079	ENDIF
1080	ENDDO
1081
1082
1083	!! Computation of the apparent ground heat flux (> towards the soil) and
1084	!! apparent surface heat capacity, used at the next timestep by enerbil to
1085	!! compute the surface temperature.
1086	DO ji = 1,kjpindex
1087	soilflx_nosnow(ji) = zdz1(ji,1) * (cgrnd(ji,1) + (dgrnd(ji,1)-1.) * ptn(ji,1))
1088	soilcap_nosnow(ji) = (zdz2(ji,1) * dt_sechiba + dt_sechiba * (un - dgrnd(ji,1)) * zdz1(ji,1))
1089	z1(ji) = lambda * (un - dgrnd(ji,1)) + un
1090	soilcap_nosnow(ji) = soilcap_nosnow(ji) / z1(ji)
1091	soilflx_nosnow(ji) = soilflx_nosnow(ji) + &
1092	& soilcap_nosnow(ji) * (ptn(ji,1) * z1(ji) - lambda * cgrnd(ji,1) - temp_sol_new(ji)) / dt_sechiba
1093	ENDDO
1094
1095	!! Add snow fraction
1096	IF ( ok_explicitsnow ) THEN
1097	! Using an effective heat capacity and heat flux by a simple pondering of snow and soil fraction
1098	DO ji = 1, kjpindex
1099	soilcap(ji) = snowcap(ji)frac_snow_veg(ji)(1-totfrac_nobio(ji))+ & ! weights related to snow cover fraction on vegetation
1100	soilcap_nosnow(ji)SUM(frac_snow_nobio(ji,:))totfrac_nobio(ji)+ & ! weights related to SCF on nobio
1101	soilcap_nosnow(ji)(1-(frac_snow_veg(ji)(1-totfrac_nobio(ji))+SUM(frac_snow_nobio(ji,:))*totfrac_nobio(ji))) ! weights related to non snow fraction
1102	soilflx(ji) = snowflx(ji)frac_snow_veg(ji)(1-totfrac_nobio(ji))+ & ! weights related to snow cover fraction on vegetation
1103	soilflx_nosnow(ji)SUM(frac_snow_nobio(ji,:))totfrac_nobio(ji)+ & ! weights related to SCF on nobio
1104	soilflx_nosnow(ji)(1-(frac_snow_veg(ji)(1-totfrac_nobio(ji))+SUM(frac_snow_nobio(ji,:))*totfrac_nobio(ji))) ! weights related to non snow fraction
1105	ENDDO
1106	ELSE
1107	! Do not consider snow fraction
1108	soilcap(:)=soilcap_nosnow(:)
1109	soilflx(:)=soilflx_nosnow(:)
1110	END IF
1111
1112	! Total snow flux (including snow on vegetated and bare soil and nobio areas)
1113	DO ji = 1, kjpindex
1114	snowflxtot(ji) = snowflx(ji)frac_snow_veg(ji)(1-totfrac_nobio(ji)) + &
1115	soilflx_nosnow(ji)SUM(frac_snow_nobio(ji,:))totfrac_nobio(ji)
1116	ENDDO
1117	CALL xios_orchidee_send_field("snowflxtot",snowflxtot(:))
1118
1119	IF (printlev>=3) WRITE (numout,*) ' thermosoil_coef done '
1120
1121	END SUBROUTINE thermosoil_coef
1122
1123
1124	!! ================================================================================================================================
1125	!! SUBROUTINE : thermosoil_profile
1126	!!
1127	!>\BRIEF In this routine solves the numerical soil thermal scheme, ie calculates the new soil temperature profile;
1128	!!
1129	!!
1130	!! DESCRIPTION : The calculation of the new soil temperature profile is based on
1131	!! the cgrnd and dgrnd values from the previous timestep and the surface temperature Ts aka temp_sol_new. (see detailed
1132	!! explanation in the header of the thermosoil module or in the reference).\n
1133	!! T(k+1)=cgrnd(k)+dgrnd(k)*T(k)\n
1134	!! -- EQ1 --\n
1135	!! Ts=(1+lambda)T(1) -lambdaT(2)\n
1136	!! -- EQ2--\n
1137	!!
1138	!! RECENT CHANGE(S) : None
1139	!!
1140	!! MAIN OUTPUT VARIABLE(S): ptn (soil temperature profile on the thermal axis),
1141	!! stempdiag (soil temperature profile on the diagnostic axis)
1142	!!
1143	!! REFERENCE(S) :
1144	!! - Hourdin, F. (1992). Study and numerical simulation of the general circulation of planetary atmospheres,
1145	!! Ph.D. thesis, Paris VII University. Remark: the part of F. Hourdin's PhD thesis relative to the thermal
1146	!! integration scheme has been scanned and is provided along with the documentation, with name :
1147	!! Hourdin_1992_PhD_thermal_scheme.pdf
1148	!!
1149	!! FLOWCHART : None
1150	!! \n
1151	!_ ================================================================================================================================
1152
1153	SUBROUTINE thermosoil_profile (kjpindex, temp_sol_new, &
1154	frac_snow_veg, frac_snow_nobio, totfrac_nobio, &
1155	ptn, stempdiag, snowtemp, &
1156	cgrnd_snow, dgrnd_snow)
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),DIMENSION (kjpindex), INTENT (in) :: temp_sol_new !! Surface temperature at the present time-step
1164	!! @tex ($K$) @endtex
1165	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: frac_snow_veg !! Snow cover fraction on vegeted area
1166	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT(in) :: frac_snow_nobio!! Snow cover fraction on non-vegeted area
1167	REAL(r_std),DIMENSION (kjpindex),INTENT(in) :: totfrac_nobio !! Total fraction of continental ice+lakes+cities+...
1168	!! (unitless,0-1)
1169	REAL(r_std),DIMENSION (kjpindex,nsnow), INTENT(in) :: snowtemp !! Snow temperature (K)
1170	REAL(r_std),DIMENSION (kjpindex,nsnow), INTENT(in) :: cgrnd_snow !! Integration coefficient for snow numerical scheme
1171	REAL(r_std),DIMENSION (kjpindex,nsnow), INTENT(in) :: dgrnd_snow !! Integration coefficient for snow numerical scheme
1172
1173	!! 0.3 Modified variables
1174
1175
1176	!! 0.2 Output variables
1177	REAL(r_std),DIMENSION (kjpindex,ngrnd), INTENT (out) :: ptn !! vertically discretized soil temperatures
1178	!! @tex ($K$) @endtex
1179	REAL(r_std),DIMENSION (kjpindex,nslm), INTENT (out) :: stempdiag !! diagnostic temperature profile
1180	!! @tex ($K$) @endtex
1181
1182	!! 0.4 Local variables
1183
1184	INTEGER(i_std) :: ji, jg
1185	REAL(r_std) :: temp_sol_eff !! effective surface temperature including snow and soil
1186
1187	!_ ================================================================================================================================
1188
1189	!! 1. Computes the soil temperatures ptn.
1190
1191	!! 1.1. ptn(jg=1) using EQ1 and EQ2
1192	DO ji = 1,kjpindex
1193
1194	IF ( ok_explicitsnow ) THEN
1195	! Using an effective surface temperature by a simple pondering
1196	temp_sol_eff=snowtemp(ji,nsnow)frac_snow_veg(ji)(1-totfrac_nobio(ji))+ & ! weights related to snow cover fraction on vegetation
1197	temp_sol_new(ji)SUM(frac_snow_nobio(ji,:))totfrac_nobio(ji)+ & ! weights related to SCF on nobio
1198	temp_sol_new(ji)(1-(frac_snow_veg(ji)(1-totfrac_nobio(ji))+SUM(frac_snow_nobio(ji,:))*totfrac_nobio(ji))) ! weights related to non snow fraction
1199	! Soil temperature calculation with explicit snow if there is snow on the ground
1200	ptn(ji,1) = cgrnd_snow(ji,nsnow) + dgrnd_snow(ji,nsnow) * temp_sol_eff
1201	ELSE
1202	! Standard soil temperature calculation
1203	ptn(ji,1) = (lambda * cgrnd(ji,1) + temp_sol_new(ji)) / (lambda *(un - dgrnd(ji,1)) + un)
1204	ENDIF
1205	ENDDO
1206
1207	!! 1.2. ptn(jg=2:ngrnd) using EQ1.
1208	DO jg = 1,ngrnd-1
1209	DO ji = 1,kjpindex
1210	ptn(ji,jg+1) = cgrnd(ji,jg) + dgrnd(ji,jg) * ptn(ji,jg)
1211	ENDDO
1212	ENDDO
1213
1214	!! 2. Assigne the soil temperature to the output variable. It is already on the right axis.
1215	stempdiag(:,:) = ptn(:,1:nslm)
1216
1217	IF (printlev>=3) WRITE (numout,*) ' thermosoil_profile done '
1218
1219	END SUBROUTINE thermosoil_profile
1220
1221	!================================================================================================================================
1222	!! SUBROUTINE : thermosoil_cond
1223	!!
1224	!>\BRIEF Calculate soil thermal conductivity.
1225	!!
1226	!! DESCRIPTION : This routine computes soil thermal conductivity
1227	!! Code introduced from NOAH LSM.
1228	!!
1229	!! RECENT CHANGE(S) : None
1230	!!
1231	!! MAIN OUTPUT VARIABLE(S): cnd
1232	!!
1233	!! REFERENCE(S) :
1234	!! Farouki, O.T.,1986: Thermal Properties of Soils. Series on Rock
1235	!! and Soil Mechanics, Vol. 11, Trans Tech, 136 PP.
1236	!! Johansen, O., 1975: Thermal Conductivity of Soils. Ph.D. Thesis,
1237	!! University of Trondheim,
1238	!! Peters-Lidard, C. D., Blackburn, E., Liang, X., & Wood, E. F.,
1239	!! 1998: The effect of soil thermal conductivity
1240	!! Parameterization on Surface Energy fluxes
1241	!! and Temperatures. J. of The Atmospheric Sciences,
1242	!! Vol. 55, pp. 1209-1224.
1243	!! Modify histroy:
1244	!!
1245	!! FLOWCHART : None
1246	!! \n
1247	!_
1248	!================================================================================================================================
1249
1250	SUBROUTINE thermosoil_cond (kjpindex, njsc, smc, qz, smcmax, sh2o, cnd)
1251
1252	!! 0. Variables and parameter declaration
1253
1254	!! 0.1 Input variables
1255	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (unitless)
1256	INTEGER(i_std), DIMENSION (kjpindex), INTENT (in) :: njsc !! Index of the dominant soil textural class in the grid cell (1-nscm, unitless)
1257	REAL(r_std), DIMENSION (kjpindex,ngrnd), INTENT(IN) :: smc !! Volumetric Soil Moisture Content (m3/m3)
1258	REAL(r_std), DIMENSION (nscm), INTENT(IN) :: qz !! Quartz Content (Soil Type Dependent) (0-1)
1259	REAL(r_std), DIMENSION (nscm), INTENT(IN) :: smcmax !! Soil Porosity (0-1)
1260	REAL(r_std), DIMENSION (kjpindex,ngrnd), INTENT(IN) :: sh2o !! Unfrozen Soil Moisture Content; Frozen Soil Moisture = smc - sh2o
1261
1262	!! 0.2 Output variables
1263	REAL(r_std), DIMENSION (kjpindex,ngrnd), INTENT(OUT) :: cnd !! Soil Thermal Conductivity (W/m/k)
1264
1265	!! 0.3 Modified variables
1266
1267	!! 0.4 Local variables
1268	REAL(r_std), DIMENSION (kjpindex,ngrnd) :: ake !! Kersten Number (unitless)
1269	REAL(r_std), DIMENSION (kjpindex,ngrnd) :: thksat !! Saturated Thermal Conductivity (W/m/k)
1270	REAL(r_std), DIMENSION (kjpindex,ngrnd) :: satratio !! Degree of Saturation (0-1)
1271	REAL(r_std), DIMENSION (kjpindex,ngrnd) :: xu !! Unfrozen Volume For Saturation (0-1)
1272	REAL(r_std), DIMENSION (kjpindex,ngrnd) :: xunfroz !! Unfrozon Volume Fraction (0-1)
1273	REAL(r_std) :: thko !! Thermal Conductivity for Other Ssoil Components (W/m/k)
1274	REAL(r_std) :: gammd !! Dry Dendity (kg/m3)
1275	REAL(r_std) :: thkdry !! Dry Thermal Conductivity (W/m/k)
1276	REAL(r_std) :: thks !! Thermal Conductivity for the Solids Combined (Quartz + Other) (W/m/k)
1277	REAL(r_std), PARAMETER :: THKICE = 2.2 !! Ice Thermal Conductivity (W/m/k)
1278	REAL(r_std), PARAMETER :: THKQTZ = 7.7 !! Thermal Conductivity for Quartz (W/m/k)
1279	REAL(r_std), PARAMETER :: THKW = 0.57 !! Water Thermal Conductivity (W/m/k)
1280	INTEGER(i_std) :: ji, jg, jst
1281
1282	!_================================================================================================================================
1283
1284	!! 1. Dry and Saturated Thermal Conductivity.
1285
1286	DO ji = 1,kjpindex
1287	jst = njsc(ji)
1288
1289	!! 1.1. Dry density (Kg/m3) and Dry thermal conductivity (W.M-1.K-1)
1290	gammd = (1. - smcmax(jst))*2700.
1291	thkdry = (0.135* gammd+ 64.7)/ (2700. - 0.947* gammd)
1292
1293	!! 1.2. thermal conductivity of "other" soil components
1294	IF (qz(jst) > 0.2) THEN
1295	thko = 2.0
1296	ELSEIF (qz(jst) <= 0.2) THEN
1297	thko = 3.0
1298	ENDIF
1299
1300	!! 1.3. Thermal conductivity of solids
1301	thks = (THKQTZ ** qz(jst))* (thko ** (1. - qz(jst)))
1302
1303	DO jg = 1,ngrnd
1304	!! 1.4. saturation ratio
1305	satratio(ji,jg) = smc(ji,jg) / smcmax(jst)
1306
1307	!! 1.5. Saturated Thermal Conductivity (thksat)
1308	IF ( smc(ji,jg) > min_sechiba ) THEN
1309	xunfroz(ji,jg) = sh2o(ji,jg) / smc(ji,jg) ! Unfrozen Fraction (From i.e., 100%Liquid, to 0. (100% Frozen))
1310	xu(ji,jg) = xunfroz(ji,jg) * smcmax(jst) ! Unfrozen volume for saturation (porosity*xunfroz)
1311	thksat(ji,jg) = thks ** (1. - smcmax(jst))* THKICE ** (smcmax(jst) - xu(ji,jg))* THKW ** (xu(ji,jg))
1312	ELSE
1313	! this value will not be used since ake=0 for this case
1314	thksat(ji,jg)=0
1315	END IF
1316	END DO ! DO jg = 1,ngrnd
1317
1318	!! 2. Kersten Number (ake)
1319	DO jg = 1,ngrnd
1320	IF ( (sh2o(ji,jg) + 0.0005) < smc(ji,jg) ) THEN
1321	! Frozen
1322	ake(ji,jg) = satratio(ji,jg)
1323	ELSE
1324	! Unfrozen
1325	! Eq 11 in Peters-Lidard et al., 1998
1326	IF ( satratio(ji,jg) > 0.1 ) THEN
1327	IF ((jst < 4 .AND. soil_classif == 'usda') .OR. (jst == 1 .AND. soil_classif == 'zobler') ) THEN
1328	! Coarse
1329	ake(ji,jg) = 0.7 * LOG10 (SATRATIO(ji,jg)) + 1.0
1330	ELSE
1331	! Fine
1332	ake(ji,jg) = LOG10 (satratio(ji,jg)) + 1.0
1333	ENDIF
1334	ELSEIF ( satratio(ji,jg) > 0.05 .AND. satratio(ji,jg) <= 0.1 ) THEN
1335	IF ((jst < 4 .AND. soil_classif == 'usda') .OR. (jst == 1 .AND. soil_classif == 'zobler') ) THEN
1336	! Coarse
1337	ake(ji,jg) = 0.7 * LOG10 (satratio(ji,jg)) + 1.0
1338	ELSE
1339	! Fine
1340	ake(ji,jg) = 0.0
1341	ENDIF
1342	ELSE
1343	ake(ji,jg) = 0.0 ! use k = kdry
1344	END IF
1345	END IF
1346	END DO ! DO jg = 1,ngrnd
1347
1348	!! 3. Thermal conductivity (cnd)
1349	DO jg = 1,ngrnd
1350	cnd(ji,jg) = ake(ji,jg) * (thksat(ji,jg) - thkdry) + thkdry
1351	END DO ! DO jg = 1,ngrnd
1352
1353	END DO !DO ji = 1,kjpindex
1354
1355	END SUBROUTINE thermosoil_cond
1356
1357
1358	!! ================================================================================================================================
1359	!! SUBROUTINE : thermosoil_humlev
1360	!!
1361	!>\BRIEF Interpolate variables from the hydrology layers to the thermodynamic layers
1362	!!
1363	!! DESCRIPTION : Interpolate the volumetric soil moisture content from the node to the interface of the layer.
1364	!! The values for the deep layers in thermosoil where hydrology is not existing are constant.
1365	!! No interpolation is needed for the total soil moisture content and for the soil saturation degree.
1366	!!
1367	!! RECENT CHANGE(S) : None
1368	!!
1369	!! MAIN OUTPUT VARIABLE(S): mc_layt, mcl_layt, tmc_layt, shum_ngrnd_perma
1370	!!
1371	!! REFERENCE(S) : None
1372	!!
1373	!! FLOWCHART : None
1374	!! \n
1375	!_ ================================================================================================================================
1376	SUBROUTINE thermosoil_humlev(kjpindex, shumdiag_perma, mc_layh, mcl_layh, tmc_layh)
1377
1378	!! 0. Variables and parameter declaration
1379
1380	!! 0.1 Input variables
1381
1382	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (unitless)
1383	REAL(r_std),DIMENSION (kjpindex,nslm), INTENT (in) :: shumdiag_perma !! Soil saturation degree on the diagnostic axis (0-1, unitless)
1384	REAL(r_std),DIMENSION (kjpindex,nslm), INTENT (in) :: mc_layh !! Volumetric soil moisture content for each layer in hydrol at nodes(liquid+ice) [m/s]
1385	REAL(r_std),DIMENSION (kjpindex,nslm), INTENT (in) :: mcl_layh !! Volumetric soil moisture content for each layer in hydrol at nodes(liquid) [m/s]
1386	REAL(r_std),DIMENSION (kjpindex,nslm), INTENT (in) :: tmc_layh !! Total soil moisture content for each layer in hydrol(liquid+ice) [mm]
1387
1388	!! 0.2 Output variables
1389
1390	!! 0.3 Modified variables
1391
1392	!! 0.4 Local variables
1393	INTEGER(i_std) :: ji, jd
1394
1395	!_ ================================================================================================================================
1396
1397	IF (printlev >= 4) WRITE(numout,*) 'Start thermosoil_humlev'
1398
1399	! The values for the deep layers in thermosoil where hydrology is not existing are constant.
1400	! For exemple if thermosoil uses 8m, and hydrol uses 2m vertical discretization,
1401	! the values between 2m and 8m are constant.
1402	! The moisture computed in hydrol is at the nodes (except for the
1403	! top and bottom layer which are at interfaces)
1404	! A linear interpolation is applied to obtain the moisture values at
1405	! the interfaces (mc_layt), from the mc_layh at the nodes
1406
1407	DO ji=1,kjpindex
1408	DO jd = 1, nslm
1409	IF(jd == 1) THEN ! the moisture at the 1st interface mc_layh(1) is at the surface, no interpolation
1410	mc_layt(ji,jd) = mc_layh(ji,jd)
1411	mcl_layt(ji,jd) = mcl_layh(ji,jd)
1412	ELSEIF(jd == 2) THEN !! the mc_layt at the 2nd interface is interpolated using mc_layh(1) at surface and mc_layh(2) at the node
1413	mc_layt(ji, jd) = mc_layh(ji,jd-1)*(znt(jd)-zlt(jd-1))/(znt(jd)-0.0) + &
1414	mc_layh(ji, jd)*(zlt(jd-1)-0.0)/(znt(jd)-0.0)
1415	mcl_layt(ji, jd) = mcl_layh(ji,jd-1)*(znt(jd)-zlt(jd-1))/(znt(jd)-0.0) + &
1416	mcl_layh(ji, jd)*(zlt(jd-1)-0.0)/(znt(jd)-0.0)
1417	ELSEIF(jd == nslm) THEN ! the mc_layt at the nslm interface is interpolated using mc_layh(nslm) and mc_layh(nslm-1)
1418	mc_layt(ji, jd) = mc_layh(ji,jd-1)*(zlt(jd)-zlt(jd-1))/(zlt(jd)-znt(jd-1)) + &
1419	mc_layh(ji,jd)*(zlt(jd-1)-znt(jd-1))/(zlt(jd)-znt(jd-1))
1420	mcl_layt(ji, jd) = mcl_layh(ji,jd-1)*(zlt(jd)-zlt(jd-1))/(zlt(jd)-znt(jd-1)) + &
1421	mcl_layh(ji,jd)*(zlt(jd-1)-znt(jd-1))/(zlt(jd)-znt(jd-1))
1422	ELSE ! the mc_layt at the other interfaces are interpolated using mc_layh at adjacent nodes.
1423	mc_layt(ji, jd) = mc_layh(ji, jd-1)(1-dz5(jd-1)) + mc_layh(ji,jd)dz5(jd-1)
1424	mcl_layt(ji, jd) = mcl_layh(ji, jd-1)(1-dz5(jd-1)) + mcl_layh(ji,jd)dz5(jd-1)
1425	ENDIF
1426
1427	shum_ngrnd_perma(ji,jd) = shumdiag_perma(ji,jd)
1428	tmc_layt(ji,jd) = tmc_layh(ji,jd)
1429	ENDDO
1430
1431	! The deep layers in thermosoil where hydro is not existing
1432	DO jd = nslm+1, ngrnd
1433	shum_ngrnd_perma(ji,jd) = shumdiag_perma(ji,nslm)
1434	mc_layt(ji,jd) = mc_layh(ji,nslm)
1435	mcl_layt(ji,jd) = mcl_layh(ji,nslm)
1436	tmc_layt(ji,jd) = tmc_layh(ji,nslm)/dlt(nslm) *dlt(jd)
1437	ENDDO
1438	ENDDO
1439
1440	IF (printlev >= 4) WRITE(numout,*) 'thermosoil_humlev done'
1441
1442	END SUBROUTINE thermosoil_humlev
1443
1444
1445	!! ================================================================================================================================
1446	!! SUBROUTINE : thermosoil_energy_diag
1447	!!
1448	!>\BRIEF Calculate diagnostics
1449	!!
1450	!! DESCRIPTION : Calculate diagnostic variables coldcont_incr and coldcont_incr
1451	!!
1452	!! RECENT CHANGE(S) : None
1453	!!
1454	!! MAIN OUTPUT VARIABLE(S) :
1455	!!
1456	!! REFERENCE(S) : None
1457	!!
1458	!! FLOWCHART : None
1459	!! \n
1460	!_ ================================================================================================================================
1461
1462	SUBROUTINE thermosoil_energy_diag(kjpindex, temp_sol_new, soilcap)
1463
1464	!! 0. Variables and parameter declaration
1465
1466	!! 0.1 Input variables
1467
1468	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (unitless)
1469	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_sol_new !! Surface temperature at the present time-step, Ts
1470	!! @tex ($K$) @endtex
1471	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: soilcap !! Apparent surface heat capacity
1472	!! @tex ($J m^{-2} K^{-1}$) @endtex,
1473	!! see eq. A29 of F. Hourdin\'s PhD thesis.
1474
1475	!! 0.2 Output variables
1476
1477	!! 0.3 Modified variables
1478
1479	!! 0.4 Local variables
1480
1481	INTEGER(i_std) :: ji, jg
1482	!_ ================================================================================================================================
1483
1484	! Sum up the energy content of all layers in the soil.
1485	DO ji = 1, kjpindex
1486
1487	IF (pcapa_en(ji,1) .LE. sn_capa) THEN
1488
1489	! Verify the energy conservation in the surface layer
1490	coldcont_incr(ji) = soilcap(ji) * (temp_sol_new(ji) - temp_sol_beg(ji))
1491	surfheat_incr(ji) = zero
1492	ELSE
1493
1494	! Verify the energy conservation in the surface layer
1495	surfheat_incr(ji) = soilcap(ji) * (temp_sol_new(ji) - temp_sol_beg(ji))
1496	coldcont_incr(ji) = zero
1497	ENDIF
1498	ENDDO
1499
1500	! Save temp_sol_new to be used at next timestep
1501	temp_sol_beg(:) = temp_sol_new(:)
1502
1503	END SUBROUTINE thermosoil_energy_diag
1504
1505
1506
1507	!! ================================================================================================================================
1508	!! SUBROUTINE : thermosoil_readjust
1509	!!
1510	!>\BRIEF
1511	!!
1512	!! DESCRIPTION : Energy conservation : Correction to make sure that the same latent heat is released and
1513	!! consumed during freezing and thawing
1514	!!
1515	!! RECENT CHANGE(S) : None
1516	!!
1517	!! MAIN OUTPUT VARIABLE(S): ptn (soil temperature profile on the thermal axis),
1518	!!
1519	!! REFERENCE(S) :
1520	!!
1521	!! FLOWCHART : None
1522	!! \n
1523	!_ ================================================================================================================================
1524
1525	SUBROUTINE thermosoil_readjust(kjpindex, ptn)
1526
1527	!! 0. Variables and parameter declaration
1528
1529	!! 0.1 Input variables
1530	INTEGER(i_std), INTENT(in) :: kjpindex
1531
1532	!! 0.2 Modified variables
1533	REAL(r_std),DIMENSION(kjpindex,ngrnd),INTENT(inout) :: ptn
1534
1535	!! 0.3 Local variables
1536	INTEGER(i_std) :: ji, jg
1537	INTEGER(i_std) :: lev3m !! Closest interface level to 3m
1538	REAL(r_std) :: ptn_tmp
1539
1540	! The energy is spread over the layers down to approximatly 3m
1541	! Find the closest level to 3m. It can be below or above 3m.
1542	lev3m=MINLOC(ABS(zlt(:)-3.0),dim=1)
1543	IF (printlev >= 3) WRITE(numout,*) 'In thermosoil_adjust: lev3m=',lev3m, ' zlt(lev3m)=', zlt(lev3m)
1544
1545	DO jg=1, ngrnd
1546	DO ji=1, kjpindex
1547	! All soil latent energy is put into e_soil_lat(ji)
1548	! because the variable soil layers make it difficult to keep track of all
1549	! layers in this version
1550	! NOTE : pcapa has unit J/K/m3 and pcappa_supp has J/K
1551	e_soil_lat(ji)=e_soil_lat(ji)+pcappa_supp(ji,jg)*(ptn(ji,jg)-ptn_beg(ji,jg))
1552	END DO
1553	END DO
1554
1555	DO ji=1, kjpindex
1556	IF (e_soil_lat(ji).GT.min_sechiba.AND.MINVAL(ptn(ji,:)).GT.ZeroCelsius+fr_dT/2.) THEN
1557	! The soil is thawed: we spread the excess of energy over the uppermost lev3m levels
1558	! Here we increase the temperatures
1559	DO jg=1, lev3m
1560	ptn_tmp=ptn(ji,jg)
1561
1562	ptn(ji,jg)=ptn(ji,jg)+MIN(e_soil_lat(ji)/pcapa(ji,jg)/zlt(lev3m), 0.5)
1563	e_soil_lat(ji)=e_soil_lat(ji)-(ptn(ji,jg)-ptn_tmp)pcapa(ji,jg)dlt(jg)
1564	ENDDO
1565	ELSE IF (e_soil_lat(ji).LT.-min_sechiba.AND.MINVAL(ptn(ji,:)).GT.ZeroCelsius+fr_dT/2.) THEN
1566	! The soil is thawed
1567	! Here we decrease the temperatures
1568	DO jg=1, lev3m
1569	ptn_tmp=ptn(ji,jg)
1570	ptn(ji,jg)=MAX(ZeroCelsius+fr_dT/2., ptn_tmp+e_soil_lat(ji)/pcapa(ji,jg)/zlt(lev3m))
1571	e_soil_lat(ji)=e_soil_lat(ji)+(ptn_tmp-ptn(ji,jg))pcapa(ji,jg)dlt(jg)
1572	END DO
1573	END IF
1574	END DO
1575
1576	END SUBROUTINE thermosoil_readjust
1577
1578	!-------------------------------------------------------------------
1579
1580
1581
1582	!! ================================================================================================================================
1583	!! SUBROUTINE : thermosoil_getdiff
1584	!!
1585	!>\BRIEF Computes soil and snow heat capacity and conductivity
1586	!!
1587	!! DESCRIPTION : Computation of the soil thermal properties; snow properties are also accounted for
1588	!!
1589	!! RECENT CHANGE(S) : None
1590	!!
1591	!! MAIN OUTPUT VARIABLE(S):
1592	!!
1593	!! REFERENCE(S) :
1594	!!
1595	!! FLOWCHART : None
1596	!! \n
1597	!_ ================================================================================================================================
1598
1599	SUBROUTINE thermosoil_getdiff( kjpindex, snow, ptn, njsc, snowrho, snowtemp, pb )
1600
1601	!! 0. Variables and parameter declaration
1602
1603	!! 0.1 Input variables
1604	INTEGER(i_std),INTENT(in) :: kjpindex
1605	REAL(r_std),DIMENSION(kjpindex),INTENT (in) :: snow !! Snow mass
1606	INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: njsc !! Index of the dominant soil textural class in the grid cell (1-nscm, unitless)
1607	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(in) :: snowrho !! Snow density
1608	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(in) :: snowtemp !! Snow temperature (K)
1609	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: pb !! Surface presure (hPa)
1610	REAL(r_std),DIMENSION(kjpindex,ngrnd),INTENT(in) :: ptn !! Soil temperature profile
1611
1612	!! 0.3 Local variables
1613	REAL :: xx !! Unfrozen fraction of the soil
1614	REAL(r_std), DIMENSION(kjpindex) :: snow_h
1615	REAL(r_std), DIMENSION(kjpindex,ngrnd) :: zx1, zx2
1616	REAL :: cap_iw !! Heat capacity of ice/water mixture
1617	INTEGER :: ji,jg
1618	INTEGER :: jst
1619	REAL(r_std), DIMENSION (kjpindex,ngrnd) :: pcapa_tmp !! soil heat capacity (J/m3/K)
1620	REAL(r_std), DIMENSION (kjpindex,ngrnd) :: pcapa_spec !! SPECIFIC soil heat capacity (J/kg/K)
1621	REAL(r_std) :: rho_tot !! Soil density (kg/m3)
1622
1623	pcapa_tmp(:,:) = 0.0
1624
1625	!! Computes soil heat capacity and conductivity
1626	DO ji = 1,kjpindex
1627	IF (.NOT. ok_explicitsnow) THEN
1628	! 1. Determine the fractions of snow and soil
1629	snow_h(ji) = snow(ji) / sn_dens
1630
1631	!
1632	! 1.1. The first level
1633	!
1634	IF ( snow_h(ji) .GT. zlt(1) ) THEN
1635	! the 1st level is in the snow => the 1st layer is entirely snow
1636	zx1(ji,1) = 1.
1637	zx2(ji,1) = 0.
1638	ELSE IF ( snow_h(ji) .GT. zero ) THEN
1639	! the 1st level is beyond the snow and the snow is present
1640	zx1(ji,1) = snow_h(ji) / zlt(1)
1641	zx2(ji,1) = ( zlt(1) - snow_h(ji)) / zlt(1)
1642	ELSE
1643	! there is no snow at all, quoi ;-)
1644	zx1(ji,1) = 0.
1645	zx2(ji,1) = 1.
1646	ENDIF
1647
1648	!
1649	! 1.2. The other levels
1650	!
1651	DO jg = 2, ngrnd
1652	IF ( snow_h(ji) .GT. zlt(jg) ) THEN
1653	! the current level is in the snow => the current layer is entirely snow
1654	zx1(ji,jg) = 1.
1655	zx2(ji,jg) = 0.
1656	ELSE IF ( snow_h(ji) .GT. zlt(jg-1) ) THEN
1657	! the current layer is partially snow
1658	zx1(ji,jg) = (snow_h(ji) - zlt(jg-1)) / (zlt(jg) - zlt(jg-1))
1659	zx2(ji,jg) = ( zlt(jg) - snow_h(ji)) / (zlt(jg) - zlt(jg-1))
1660	ELSE
1661	! both levels are out of snow => the current layer is entirely soil
1662	zx1(ji,jg) = 0.
1663	zx2(ji,jg) = 1.
1664	ENDIF
1665	ENDDO
1666	ELSE
1667	zx1(ji,:) = 0.
1668	zx2(ji,:) = 1.
1669	END IF
1670
1671	DO jg = 1, ngrnd
1672	jst = njsc(ji)
1673	pcapa_tmp(ji, jg) = so_capa_dry_ns(jst) * (1-mcs(jst)) + water_capa * tmc_layt(ji,jg)/mille/dlt(jg)
1674	!
1675	! 2. Calculate volumetric heat capacity with allowance for permafrost
1676	! 2.1. soil heat capacity depending on temperature and humidity
1677	! For SP6MIP we also diagnose a specific heat capacity (pcapa_spec),
1678	! which requires to calculate the total density of the soil (rho_tot), always >> 0
1679
1680	IF (ptn(ji,jg) .LT. ZeroCelsius-fr_dT/2.) THEN
1681	! frozen soil
1682	profil_froz(ji,jg) = 1.
1683	pcappa_supp(ji,jg)= 0.
1684	pcapa(ji, jg) = so_capa_dry_ns(jst) * (1-mcs(jst)) + so_capa_ice * tmc_layt(ji,jg) / mille / dlt(jg)
1685	rho_tot = rho_soil * (1-mcs(jst)) + rho_ice * tmc_layt(ji,jg) / mille / dlt(jg)
1686	pcapa_spec(ji, jg) = pcapa(ji, jg) / rho_tot
1687
1688	ELSEIF (ptn(ji,jg) .GT. ZeroCelsius+fr_dT/2.) THEN
1689	! unfrozen soil
1690	pcapa(ji, jg) = pcapa_tmp(ji, jg)
1691	profil_froz(ji,jg) = 0.
1692	pcappa_supp(ji,jg)= 0.
1693	rho_tot = rho_soil * (1-mcs(jst)) + rho_water * tmc_layt(ji,jg)/mille/dlt(jg)
1694	pcapa_spec(ji, jg) = pcapa(ji, jg) / rho_tot
1695	ELSE
1696	! xx is the unfrozen fraction of soil water
1697	xx = (ptn(ji,jg)-(ZeroCelsius-fr_dT/2.)) / fr_dT
1698	profil_froz(ji,jg) = (1. - xx)
1699
1700	! net heat capacity of the ice/water mixture
1701	cap_iw = xx * so_capa_wet + (1.-xx) * so_capa_ice
1702	pcapa(ji, jg) = so_capa_dry_ns(jst) * (1-mcs(jst)) + &
1703	water_capa * tmc_layt(ji,jg)/mille / dlt(jg) * xx + &
1704	so_capa_ice * tmc_layt(ji,jg) / mille/dlt(jg) * (1.-xx)
1705	rho_tot = rho_soil* (1-mcs(jst)) + &
1706	rho_water * tmc_layt(ji,jg)/mille / dlt(jg) * xx + &
1707	rho_ice * tmc_layt(ji,jg) / mille/dlt(jg) * (1.-xx)
1708	pcapa_spec(ji, jg) = pcapa(ji, jg) / rho_tot
1709
1710	pcappa_supp(ji,jg)= shum_ngrnd_perma(ji,jg)poroslhfrho_water/fr_dTzx2(ji,jg)*dlt(jg)
1711
1712	ENDIF
1713
1714	!
1715	! 2.2. Take into account the snow and soil fractions in the layer
1716	!
1717	pcapa(ji,jg) = zx1(ji,jg) * sn_capa + zx2(ji,jg) * pcapa(ji,jg)
1718
1719	!
1720	! 2.3. Calculate the heat capacity for energy conservation check
1721	IF ( zx1(ji,jg).GT.0. ) THEN
1722	pcapa_en(ji,jg) = sn_capa
1723	ELSE
1724	pcapa_en(ji,jg) = pcapa(ji,jg)
1725	ENDIF
1726
1727	END DO
1728	ENDDO
1729
1730	! Output the specific heat capcaity for SP-MIP
1731	CALL xios_orchidee_send_field("pcapa_spec",pcapa_spec)
1732
1733	!
1734	! 3. Calculate the heat conductivity with allowance for permafrost
1735	!
1736	CALL thermosoil_cond (kjpindex, njsc, mc_layt, QZ, SMCMAX, mcl_layt, pkappa)
1737
1738	!! Computes snow heat capacity and conductivity
1739	DO ji = 1,kjpindex
1740	pcapa_snow(ji,:) = snowrho(ji,:) * xci
1741	pkappa_snow(ji,:) = (ZSNOWTHRMCOND1 + ZSNOWTHRMCOND2snowrho(ji,:)snowrho(ji,:)) + &
1742	MAX(0.0,(ZSNOWTHRMCOND_AVAP+(ZSNOWTHRMCOND_BVAP/(snowtemp(ji,:)+ &
1743	ZSNOWTHRMCOND_CVAP)))(XP00/(pb(ji)100.)))
1744	END DO
1745
1746	END SUBROUTINE thermosoil_getdiff
1747
1748	!! ================================================================================================================================
1749	!! SUBROUTINE : thermosoil_getdiff_old_thermix_with_snow
1750	!!
1751	!>\BRIEF Computes soil heat capacity and conductivity
1752	!!
1753	!! DESCRIPTION : Computes soil heat capacity and conductivity
1754	!! Special case with old snow without soil freezing
1755	!!
1756	!! RECENT CHANGE(S) : None
1757	!!
1758	!! MAIN OUTPUT VARIABLE(S):
1759	!!
1760	!! REFERENCE(S) :
1761	!!
1762	!! FLOWCHART : None
1763	!! \n
1764	!_ ================================================================================================================================
1765
1766
1767	SUBROUTINE thermosoil_getdiff_old_thermix_with_snow( kjpindex, snow, njsc )
1768
1769
1770	!! 0. Variables and parameter declaration
1771
1772	!! 0.1 Input variables
1773	INTEGER(i_std), INTENT(in) :: kjpindex
1774	REAL(r_std),DIMENSION(kjpindex),INTENT (in) :: snow
1775	INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: njsc !! Index of the dominant soil textural class in the grid cell (1-nscm, unitless)
1776
1777
1778	!! 0.2 Local variables
1779	INTEGER :: ji,jg
1780	REAL(r_std) :: snow_h !! snow_h is the snow height @tex ($m$) @endtex
1781	REAL(r_std) :: zx1, zx2 !! zx1 and zx2 are the layer fraction consisting in snow and soil respectively.
1782	INTEGER :: jst
1783	REAL(r_std), DIMENSION (kjpindex,ngrnd) :: pkappa_tmp !! soil thermal conductivity (W/m/K)
1784	REAL(r_std), DIMENSION (kjpindex,ngrnd) :: pcapa_tmp !! soil heat capacity (J/m3/K)
1785	REAL(r_std), DIMENSION (kjpindex,ngrnd) :: pkappa_wet !! wet soil thermal conductivity (W/m/K)
1786	REAL(r_std), DIMENSION (kjpindex,ngrnd) :: pcapa_wet !! wet soil heat capacity (J/m3/K)
1787	REAL(r_std), DIMENSION (kjpindex,ngrnd) :: mcs_tmp !! Saturated soil moisture (liquid+ice) (m3/m3)
1788
1789	! Computation of the soil thermal properties; snow properties are also accounted for
1790
1791	pkappa_tmp(:,:) = 0.0
1792	pcapa_tmp(:,:) = 0.0
1793	pkappa_wet(:,:) = 0.0
1794	pcapa_wet(:,:) = 0.0
1795	DO ji = 1, kjpindex
1796	jst = njsc(ji)
1797	mcs_tmp(ji,:) = mcs(jst)
1798	DO jg = 1, ngrnd
1799	pcapa_tmp(ji, jg) = so_capa_dry_ns(jst) * (1-mcs(jst)) + water_capa * tmc_layt(ji,jg)/mille/dlt(jg)
1800	pcapa_wet(ji, jg) = so_capa_dry_ns(jst) * (1-mcs(jst)) + water_capa * mcs(jst)
1801	ENDDO
1802	ENDDO
1803
1804	CALL thermosoil_cond(kjpindex, njsc, mc_layt, QZ, SMCMAX, mcl_layt, pkappa_tmp)
1805	CALL thermosoil_cond(kjpindex, njsc, mcs_tmp, QZ, SMCMAX, mcs_tmp, pkappa_wet)
1806
1807	DO ji = 1,kjpindex
1808	snow_h = snow(ji) / sn_dens
1809
1810	! First layer
1811	IF ( snow_h .GT. zlt(1) ) THEN
1812	pcapa(ji,1) = sn_capa
1813	pcapa_en(ji,1) = sn_capa
1814	pkappa(ji,1) = sn_cond
1815	ELSE IF ( snow_h .GT. zero ) THEN
1816	pcapa_en(ji,1) = sn_capa
1817	zx1 = snow_h / zlt(1)
1818	zx2 = ( zlt(1) - snow_h) / zlt(1)
1819	pcapa(ji, 1) = zx1 * sn_capa + zx2 * pcapa_wet(ji,1)
1820	pkappa(ji,1) = un / ( zx1 / sn_cond + zx2 / (pkappa_wet(ji,1)) )
1821	ELSE
1822	pkappa(ji,1) = pkappa_tmp(ji,1)
1823	pcapa(ji, 1) = pcapa_tmp(ji, 1)
1824	pcapa_en(ji, 1) = pcapa_tmp(ji, 1)
1825	ENDIF
1826
1827	! Mid layers
1828	DO jg = 2, ngrnd - 2
1829	IF ( snow_h .GT. zlt(jg) ) THEN
1830	pcapa(ji,jg) = sn_capa
1831	pkappa(ji,jg) = sn_cond
1832	pcapa_en(ji,jg) = sn_capa
1833	ELSE IF ( snow_h .GT. zlt(jg-1) ) THEN
1834	zx1 = (snow_h - zlt(jg-1)) / (zlt(jg) - zlt(jg-1))
1835	zx2 = ( zlt(jg) - snow_h) / (zlt(jg) - zlt(jg-1))
1836	pcapa_en(ji,jg) = sn_capa
1837	pcapa(ji, jg) = zx1 * sn_capa + zx2 * pcapa_wet(ji,jg)
1838	pkappa(ji,jg) = un / ( zx1 / sn_cond + zx2 / (pkappa_wet(ji,jg)))
1839	ELSE
1840	pcapa(ji, jg) = pcapa_tmp(ji, jg)
1841	pkappa(ji,jg) = pkappa_tmp(ji,jg)
1842	pcapa_en(ji, jg) = pcapa_tmp(ji, jg)
1843	ENDIF
1844	ENDDO
1845
1846	! Last two layers: These layers can not be filled with snow
1847	DO jg = ngrnd - 1, ngrnd
1848	pcapa(ji,jg) = so_capa_dry
1849	pkappa(ji,jg) = so_cond_dry
1850	pcapa_en(ji,jg) = so_capa_dry
1851	END DO
1852
1853	IF (brk_flag == 1) THEN
1854	! Bedrock flag is activated
1855	DO jg = ngrnd-1,ngrnd
1856	pcapa(ji,jg) = brk_capa
1857	pcapa_en(ji,jg) = brk_capa
1858	pkappa(ji,jg) = brk_cond
1859	ENDDO
1860	ENDIF
1861
1862	ENDDO ! DO ji = 1,kjpindex
1863
1864
1865	END SUBROUTINE thermosoil_getdiff_old_thermix_with_snow
1866
1867
1868
1869	!! ================================================================================================================================
1870	!! SUBROUTINE : thermosoil_getdiff_old_thermix_without_snow
1871	!!
1872	!>\BRIEF Computes soil and snow heat capacity and conductivity
1873	!!
1874	!! DESCRIPTION : Calculations of soil and snow thermal properties without effect of freezing. This subroutine is only
1875	!! call for the case with explictsnow activated.
1876	!!
1877	!!
1878	!! RECENT CHANGE(S) : None
1879	!!
1880	!! MAIN OUTPUT VARIABLE(S):
1881	!!
1882	!! REFERENCE(S) :
1883	!!
1884	!! FLOWCHART : None
1885	!! \n
1886	!_ ================================================================================================================================
1887
1888	SUBROUTINE thermosoil_getdiff_old_thermix_without_snow( kjpindex, njsc, snowrho, snowtemp, pb )
1889
1890	!! 0. Variables and parameter declaration
1891
1892	!! 0.1 Input variables
1893	INTEGER(i_std), INTENT(in) :: kjpindex
1894	INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: njsc !! Index of the dominant soil textural class in the grid cell (1-nscm, unitless)
1895	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(in) :: snowrho !! Snow density
1896	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(in) :: snowtemp !! Snow temperature (K)
1897	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: pb !! Surface presure (hPa)
1898
1899
1900	!! 0.1 Local variables
1901	INTEGER(i_std) :: ji,jg, jst !! Index
1902	REAL(r_std), DIMENSION (kjpindex,ngrnd) :: pcapa_tmp !! Soil heat capacity (J/m3/K)
1903
1904	!! Computes soil heat capacity and conductivity
1905	DO jg = 1,ngrnd
1906	DO ji = 1,kjpindex
1907	jst = njsc(ji)
1908	pcapa_tmp(ji, jg) = so_capa_dry_ns(jst) * (1-mcs(jst)) + water_capa * tmc_layt(ji,jg)/mille/dlt(jg)
1909	pcapa(ji,jg) = pcapa_tmp(ji, jg)
1910	pcapa_en(ji,jg) = pcapa_tmp(ji, jg)
1911	ENDDO
1912	ENDDO
1913
1914	CALL thermosoil_cond (kjpindex, njsc, mc_layt, QZ, SMCMAX, mcl_layt, pkappa)
1915
1916	IF (brk_flag == 1) THEN
1917	! Bedrock flag is activated
1918	DO jg = ngrnd-1,ngrnd
1919	DO ji = 1,kjpindex
1920	pcapa(ji,jg) = brk_capa
1921	pcapa_en(ji,jg) = brk_capa
1922	pkappa(ji,jg) = brk_cond
1923	ENDDO
1924	ENDDO
1925	ENDIF
1926
1927	!! Computes snow heat capacity and conductivity
1928	DO ji = 1,kjpindex
1929	pcapa_snow(ji,:) = snowrho(ji,:) * xci
1930	pkappa_snow(ji,:) = (ZSNOWTHRMCOND1 + ZSNOWTHRMCOND2snowrho(ji,:)snowrho(ji,:)) + &
1931	MAX(0.0,(ZSNOWTHRMCOND_AVAP+(ZSNOWTHRMCOND_BVAP/(snowtemp(ji,:)+ &
1932	ZSNOWTHRMCOND_CVAP)))(XP00/(pb(ji)100.)))
1933	END DO
1934
1935	END SUBROUTINE thermosoil_getdiff_old_thermix_without_snow
1936
1937
1938	!! ================================================================================================================================
1939	!! SUBROUTINE : thermosoil_read_reftempfile
1940	!!
1941	!>\BRIEF
1942	!!
1943	!! DESCRIPTION : Read file with longterm soil temperature
1944	!!
1945	!!
1946	!! RECENT CHANGE(S) : None
1947	!!
1948	!! MAIN OUTPUT VARIABLE(S): reftemp : Reference temerature
1949	!!
1950	!! REFERENCE(S) :
1951	!!
1952	!! FLOWCHART : None
1953	!! \n
1954	!_ ================================================================================================================================
1955	SUBROUTINE thermosoil_read_reftempfile(kjpindex,lalo,reftemp)
1956
1957	USE interpweight
1958
1959	IMPLICIT NONE
1960
1961	!! 0. Variables and parameter declaration
1962
1963	!! 0.1 Input variables
1964	INTEGER(i_std), INTENT(in) :: kjpindex
1965	REAL(r_std), DIMENSION(kjpindex,2), INTENT(in) :: lalo
1966
1967	!! 0.2 Output variables
1968	REAL(r_std), DIMENSION(kjpindex, ngrnd), INTENT(out) :: reftemp
1969
1970	!! 0.3 Local variables
1971	INTEGER(i_std) :: ib
1972	CHARACTER(LEN=80) :: filename
1973	REAL(r_std),DIMENSION(kjpindex) :: reftemp_file !! Horizontal temperature field interpolated from file [C]
1974	INTEGER(i_std),DIMENSION(kjpindex,8) :: neighbours
1975	REAL(r_std) :: vmin, vmax !! min/max values to use for the
1976	!! renormalization
1977	REAL(r_std), DIMENSION(kjpindex) :: areftemp !! Availability of data for the interpolation
1978	CHARACTER(LEN=80) :: variablename !! Variable to interpolate
1979	!! the file
1980	CHARACTER(LEN=80) :: lonname, latname !! lon, lat names in input file
1981	REAL(r_std), DIMENSION(:), ALLOCATABLE :: variabletypevals !! Values for all the types of the variable
1982	!! (variabletypevals(1) = -un, not used)
1983	CHARACTER(LEN=50) :: fractype !! method of calculation of fraction
1984	!! 'XYKindTime': Input values are kinds
1985	!! of something with a temporal
1986	!! evolution on the dx*dy matrix'
1987	LOGICAL :: nonegative !! whether negative values should be removed
1988	CHARACTER(LEN=50) :: maskingtype !! Type of masking
1989	!! 'nomask': no-mask is applied
1990	!! 'mbelow': take values below maskvals(1)
1991	!! 'mabove': take values above maskvals(1)
1992	!! 'msumrange': take values within 2 ranges;
1993	!! maskvals(2) <= SUM(vals(k)) <= maskvals(1)
1994	!! maskvals(1) < SUM(vals(k)) <= maskvals(3)
1995	!! (normalized by maskvals(3))
1996	!! 'var': mask values are taken from a
1997	!! variable inside the file (>0)
1998	REAL(r_std), DIMENSION(3) :: maskvals !! values to use to mask (according to
1999	!! `maskingtype')
2000	CHARACTER(LEN=250) :: namemaskvar !! name of the variable to use to mask
2001	REAL(r_std) :: reftemp_norefinf
2002	REAL(r_std) :: reftemp_default !! Default value
2003
2004
2005	!Config Key = SOIL_REFTEMP_FILE
2006	!Config Desc = File with climatological soil temperature
2007	!Config If = READ_REFTEMP
2008	!Config Def = reftemp.nc
2009	!Config Help =
2010	!Config Units = [FILE]
2011	filename = 'reftemp.nc'
2012	CALL getin_p('REFTEMP_FILE',filename)
2013
2014	variablename = 'temperature'
2015
2016	IF (printlev >= 1) WRITE(numout,*) "thermosoil_read_reftempfile: Read and interpolate file " &
2017	// TRIM(filename) //" for variable " //TRIM(variablename)
2018
2019	! For this case there are not types/categories. We have 'only' a continuos field
2020	! Assigning values to vmin, vmax
2021
2022	vmin = 0.
2023	vmax = 9999.
2024
2025	! For this file we do not need neightbours!
2026	neighbours = 0
2027
2028	!! Variables for interpweight
2029	! Type of calculation of cell fractions
2030	fractype = 'default'
2031	! Name of the longitude and latitude in the input file
2032	lonname = 'nav_lon'
2033	latname = 'nav_lat'
2034	! Default value when no value is get from input file
2035	reftemp_default = 1.
2036	! Reference value when no value is get from input file
2037	reftemp_norefinf = 1.
2038	! Should negative values be set to zero from input file?
2039	nonegative = .FALSE.
2040	! Type of mask to apply to the input data (see header for more details)
2041	maskingtype = 'nomask'
2042	! Values to use for the masking (here not used)
2043	maskvals = (/ undef_sechiba, undef_sechiba, undef_sechiba /)
2044	! Name of the variable with the values for the mask in the input file (only if maskkingtype='var') (here not used)
2045	namemaskvar = ''
2046
2047	CALL interpweight_2Dcont(kjpindex, 0, 0, lalo, resolution, neighbours, &
2048	contfrac, filename, variablename, lonname, latname, vmin, vmax, nonegative, maskingtype, &
2049	maskvals, namemaskvar, -1, fractype, reftemp_default, reftemp_norefinf, &
2050	reftemp_file, areftemp)
2051	IF (printlev >= 5) WRITE(numout,*)' thermosoil_read_reftempfile after interpweight_2Dcont'
2052
2053	! Copy reftemp_file temperature to all ground levels and transform into Kelvin
2054	DO ib=1, kjpindex
2055	reftemp(ib, :) = reftemp_file(ib)+ZeroCelsius
2056	END DO
2057
2058	! Write diagnostics
2059	CALL xios_orchidee_send_field("areftemp",areftemp)
2060
2061	END SUBROUTINE thermosoil_read_reftempfile
2062
2063	END MODULE thermosoil

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

Download in other formats: