Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

thermosoil.f90 @ 7746

Last change on this file since 7746 was 5745, checked in by thomas.verbeke, 6 years ago

Update GW version of ORCHIDEE-GW branche:
1) Addition of these trunk changesets:
https://forge.ipsl.jussieu.fr/orchidee/changeset/5433/trunk/ORCHIDEE
https://forge.ipsl.jussieu.fr/orchidee/changeset/5536/trunk/ORCHIDEE
https://forge.ipsl.jussieu.fr/orchidee/changeset/5573/trunk/ORCHIDEE
https://forge.ipsl.jussieu.fr/orchidee/changeset/5614/trunk/ORCHIDEE

2) Modification of wtd calculation in hydrol.f90
3) Modification of .xml files

Property svn:keywords set to Revision Date HeadURL Date Author Revision

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

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

Context Navigation

source: branches/publications/ORCHIDEE-Hillslope-r6515/src_sechiba/thermosoil.f90 @ 7746

Download in other formats: