Context Navigation

source: branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/thermosoil.f90 @ 7199

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

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

Download in other formats: