Context Navigation

source: branches/publications/ORCHIDEE_gmd-2018-261/src_sechiba/thermosoil.f90 @ 6892

Last change on this file since 6892 was 4333, checked in by josefine.ghattas, 8 years ago
Merge with trunk for revision r3623 - r3977 Some corrections for gfortran
Property svn:keywords set to `Revision Date HeadURL Date Author Revision`
File size: 101.1 KB

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

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

Download in other formats: