Context Navigation

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

Last change on this file since 7508 was 7508, checked in by josefine.ghattas, 2 years ago

Removed DRY_SOIL_HEAT_COND and related variable so_cond_dry
Replaced previous scalar variable so_capa_dry with vector variable so_capa_dry_ns and moved read from run.def to thermosoil. In the end, so_capa_dry_ns is renamed so_capa_dry.

See ticket #780

No change in results.

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

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

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

Download in other formats: