Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

thermosoil.f90 @ 7346

Last change on this file since 7346 was 7297, checked in by wei.li, 3 years ago
updated code for publication, 2021,9,25
Property svn:keywords set to `Revision Date HeadURL Date Author Revision`
File size: 204.2 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	USE ioipsl_para
78	USE xios_orchidee
79	USE constantes
80	USE time, ONLY : one_day, dt_sechiba
81	USE constantes_soil
82	USE sechiba_io_p
83	USE grid
84	USE pft_parameters_var
85	USE vertical_soil
86	USE constantes_var
87	USE interpol_help
88
89	IMPLICIT NONE
90
91	!private and public routines :
92	PRIVATE
93	PUBLIC :: thermosoil_main, thermosoil_clear, thermosoil_initialize, thermosoil_finalize, thermosoil_rotation_update
94
95	REAL(r_std), SAVE :: lambda !! See Module description
96	!$OMP THREADPRIVATE(lambda)
97	REAL(r_std), SAVE :: fz1 !! usefull constants for diverse use
98	!$OMP THREADPRIVATE(fz1)
99	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:,:) :: ptn !! vertically discretized
100	!$OMP THREADPRIVATE(ptn)
101	!! soil temperatures @tex ($K$) @endtex.
102	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: ptn_pftmean !! Different levels soil temperature, mean across all pfts
103	!$OMP THREADPRIVATE(ptn_pftmean)
104	REAL(r_std), ALLOCATABLE,SAVE, DIMENSION (:) :: dz1 !! numerical constant used in the thermal numerical
105	!! scheme @tex ($m^{-1}$) @endtex. ; it corresponds
106	!! to the coefficient @tex $d_k$ @endtex of equation
107	!! (A.12) in F. Hourdin PhD thesis.
108	!$OMP THREADPRIVATE(dz1)
109	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: z1 !! constant of the numerical scheme; it is an
110	!! intermediate buffer for the calculation of the
111	!! integration coefficients cgrnd and dgrnd.
112	!$OMP THREADPRIVATE(z1)
113	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:,:) :: cgrnd !! integration coefficient for the numerical scheme,
114	!! see eq.1
115	!$OMP THREADPRIVATE(cgrnd)
116	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:,:) :: dgrnd !! integration coefficient for the numerical scheme,
117	!! see eq.1
118	!$OMP THREADPRIVATE(dgrnd)
119	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:,:) :: pcapa !! volumetric vertically discretized soil heat
120	!$OMP THREADPRIVATE(pcapa)
121
122	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:,:) :: pkappa !! vertically discretized soil thermal conductivity
123	!! @tex ($W K^{-1} m^{-1}$) @endtex. Same as pcapa.
124	!$OMP THREADPRIVATE(pkappa)
125	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:,:) :: pcapa_en !! heat capacity used for surfheat_incr and
126	!$OMP THREADPRIVATE(pcapa_en)
127	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: pcapa_snow !! volumetric vertically discretized snow heat
128	!! capacity @tex ($J K^{-1} m^{-3}$) @endtex.
129	!$OMP THREADPRIVATE(pcapa_snow)
130	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: pkappa_snow !! vertically discretized snow thermal conductivity
131	!! @tex ($W K^{-1} m^{-1}$) @endtex.
132	!$OMP THREADPRIVATE(pkappa_snow)
133
134	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:,:) :: ptn_beg !! ptn as it is after thermosoil_profile but before thermosoil_coef,
135	!! used in thermosoil_readjust
136	!$OMP THREADPRIVATE(ptn_beg)
137	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: temp_sol_beg !! Surface temperature at previous timestep (K)
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	REAL(r_std), SAVE :: so_cond = 1.5396 !! Thermix soil layer discretization constant
148	!$OMP THREADPRIVATE(so_cond)
149	REAL(r_std), SAVE :: so_capa = 2.0514e+6 !! Thermix soil layer discretization constant
150	!$OMP THREADPRIVATE(so_capa)
151
152	! Variables related to soil freezing
153	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:,:) :: profil_froz !! Frozen fraction of the soil on hydrological levels (-)
154	!$OMP THREADPRIVATE(profil_froz)
155	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:,:) :: shum_ngrnd_permalong !! Long-term soil humidity (for permafrost) if ok_freeze_thermix ; shum_ngrnd_perma sinon.
156	!$OMP THREADPRIVATE(shum_ngrnd_permalong)
157	LOGICAL, SAVE :: ok_shum_ngrnd_permalong
158	!$OMP THREADPRIVATE(ok_shum_ngrnd_permalong)
159
160	REAL(r_std),ALLOCATABLE, SAVE, DIMENSION (:,:,:) :: pcappa_supp !! Additional heat capacity due to soil freezing for each soil layer (J/K)
161	!$OMP THREADPRIVATE(pcappa_supp)
162	REAL(r_std),ALLOCATABLE, SAVE, DIMENSION (:,:) :: e_soil_lat !! Accumulated latent heat for the whole soil (J)
163	!$OMP THREADPRIVATE(e_soil_lat)
164	REAL(r_std), ALLOCATABLE, SAVE,DIMENSION(:) :: overburden !! Information read from IPA map for option read_permafrost_mapn
165	!$OMP THREADPRIVATE(overburden)
166	REAL(r_std), ALLOCATABLE, SAVE,DIMENSION(:) :: excess_ice !! Information read from IPA map for option read_permafrost_map
167	!$OMP THREADPRIVATE(excess_ice)
168	REAL(r_std), ALLOCATABLE, SAVE,DIMENSION(:) :: permafrost !! Information read from IPA map for option read_permafrost_map
169	!$OMP THREADPRIVATE(permafrost)
170	REAL(r_std), ALLOCATABLE, SAVE,DIMENSION(:,:) :: reftemp !! Flag to initialize soil temperature using climatological temperature
171	!$OMP THREADPRIVATE(reftemp)
172	REAL(r_std), ALLOCATABLE, SAVE,DIMENSION(:,:) :: refSOC !! initialize soil organic carbon only used to calculate thermal insulating effect
173	!$OMP THREADPRIVATE(refSOC)
174	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: dz5 !! Used for numerical calculation [-]
175	!$OMP THREADPRIVATE(dz5)
176	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: mcs !! Saturation humidity [m3/m3]
177	!$OMP THREADPRIVATE(mcs)
178	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: SMCMAX !! Soil porosity [m3/m3]
179	!$OMP THREADPRIVATE(SMCMAX)
180	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: QZ !! quartz content [-]
181	!$OMP THREADPRIVATE(QZ)
182	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: so_capa_dry_ns !! Dry soil Heat capacity of soils,J.m^{-3}.K^{-1}
183	!$OMP THREADPRIVATE(so_capa_dry_ns)
184	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: mc_layt !! Volumetric soil moisture (liquid+ice) (m3/m3) on the thermodynamical levels at interface
185	!$OMP THREADPRIVATE(mc_layt)
186	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: mcl_layt !! Volumetric soil moisture (liquid) (m3/m3) on the thermodynamical levels at interface
187	!$OMP THREADPRIVATE(mcl_layt)
188	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: tmc_layt !! Total soil moisture content for each layer (liquid+ice) (mm) on the thermodynamical levels
189	!$OMP THREADPRIVATE(tmc_layt)
190	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:,:) :: mc_layt_pft !! Volumetric soil moisture (liquid+ice) (m3/m3) on the thermodynamical levels at interface
191	!$OMP THREADPRIVATE(mc_layt_pft)
192	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:,:) :: mcl_layt_pft !! Volumetric soil moisture (liquid) (m3/m3) on the thermodynamical levels at interface
193	!$OMP THREADPRIVATE(mcl_layt_pft)
194	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:,:) :: tmc_layt_pft !! Total soil moisture content for each layer (liquid+ice) (mm) on the thermodynamical levels
195	!$OMP THREADPRIVATE(tmc_layt_pft)
196	INTEGER(i_std), SAVE :: brk_flag = 0 !! Flag to consider bedrock: 0.no; 1.yes
197	!$OMP THREADPRIVATE(brk_flag)
198
199	!Vertical Permafrost Carbon
200	LOGICAL, SAVE :: use_toporganiclayer_tempdiff = .FALSE.
201	!$OMP THREADPRIVATE(use_toporganiclayer_tempdiff)
202	LOGICAL, SAVE :: use_soilc_tempdiff = .TRUE.
203	!$OMP THREADPRIVATE(use_soilc_tempdiff)
204	LOGICAL, SAVE :: use_refSOC = .TRUE. !! which SOC to use in thermix:refSOC or modeled SOC
205	!$OMP THREADPRIVATE(use_refSOC)
206	INTEGER(i_std), PARAMETER :: SOILC_METHOD_ARITHMETIC = 1
207	INTEGER(i_std), PARAMETER :: SOILC_METHOD_GEOMETRIC = 2
208	INTEGER(i_std), SAVE :: use_soilc_method = SOILC_METHOD_ARITHMETIC
209	!$OMP THREADPRIVATE(use_soilc_method)
210	!! how to average thermal conductivity of mineral soil and organic soil:
211	!! 1=arithmetic mean ; 2=geometric mean
212
213	INTEGER(i_std), PARAMETER :: SNOW_COND_METHOD_DEFAULT = 1
214	INTEGER(i_std), PARAMETER :: SNOW_COND_METHOD_DECHARME16 = 2
215	INTEGER(i_std), SAVE :: snow_cond_method = SNOW_COND_METHOD_DEFAULT !! 1: original 2: follows Decharme et al 2016
216	!$OMP THREADPRIVATE(snow_cond_method)
217
218	LOGICAL, SAVE :: satsoil = .FALSE.
219	!$OMP THREADPRIVATE(satsoil)
220
221
222
223	CONTAINS
224	!! =============================================================================================================================
225	!! SUBROUTINE : thermosoil_initialize
226	!!
227	!>\BRIEF Allocate module variables, read from restart file or initialize with default values
228	!!
229	!! DESCRIPTION : Allocate module variables, read from restart file or initialize with default values.
230	!! Call thermosoil_var_init to calculate physical constants.
231	!! Call thermosoil_coef to calculate thermal soil properties.
232	!!
233	!! RECENT CHANGE(S) : None
234	!!
235	!! REFERENCE(S) : None
236	!!
237	!! FLOWCHART : None
238	!! \n
239	!_ ==============================================================================================================================
240	SUBROUTINE thermosoil_initialize(kjit, kjpindex, lalo,neighbours, resolution,contfrac, rest_id, veget_max, &
241	shumdiag_perma, snow, thawed_humidity, soilc_total, &
242	temp_sol_new, temp_sol_new_pft, &
243	organic_layer_thick, stempdiag, soilcap, soilcap_pft, soilflx, soilflx_pft, &
244	gtemp, &
245	mc_layh, mcl_layh, tmc_layh, mc_layh_pft, mcl_layh_pft, tmc_layh_pft, njsc, &
246	frac_snow_veg,frac_snow_nobio,totfrac_nobio, &
247	snowdz, snowrho, snowtemp, lambda_snow, cgrnd_snow, dgrnd_snow, pb)
248
249	!! 0. Variable and parameter declaration
250	!! 0.1 Input variables
251	REAL(r_std),DIMENSION (kjpindex,nslm), INTENT (in) :: mc_layh !! Volumetric soil moisture content (liquid+ice) for hydrological layers, at node (m3/m3)
252	REAL(r_std),DIMENSION (kjpindex,nslm), INTENT (in) :: mcl_layh !! Volumetric soil moisture content (liquid) for hydrological layers, at node (m3/m3)
253	REAL(r_std),DIMENSION (kjpindex,nslm), INTENT (in) :: tmc_layh !! Total soil moisture content(liquid+ice) for hydrological layers (mm)
254	INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: njsc !! Index of the dominant soil textural class in the grid cell (1-nscm, unitless)
255	REAL(r_std),DIMENSION (kjpindex,nslm,nvm), INTENT (in) :: mc_layh_pft !! Volumetric soil moisture content (liquid+ice) for hydrological layers, at node (m3/m3)
256	REAL(r_std),DIMENSION (kjpindex,nslm,nvm), INTENT (in) :: mcl_layh_pft !! Volumetric soil moisture content (liquid) for hydrological layers, at node (m3/m3)
257	REAL(r_std),DIMENSION (kjpindex,nslm,nvm), INTENT (in) :: tmc_layh_pft !! Total soil moisture content(liquid+ice) for hydrological layers (mm)
258	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: frac_snow_veg !! Snow cover fraction on vegeted area
259	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT(in) :: frac_snow_nobio !! Snow cover fraction on non-vegeted area
260	REAL(r_std),DIMENSION (kjpindex),INTENT(in) :: totfrac_nobio !! Total fraction of continental ice+lakes+cities+...
261	!! (unitless,0-1)
262	REAL(r_std),DIMENSION (kjpindex,nsnow),INTENT (in) :: snowdz !! Snow depth
263	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(in) :: snowrho !! Snow density
264	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(in) :: snowtemp !! Snow temperature (K)
265	REAL(r_std), DIMENSION (kjpindex), INTENT (in) :: pb !! Surface presure (hPa)
266
267	INTEGER(i_std), INTENT (in) :: kjit !! Time step number (unitless)
268	INTEGER(i_std), INTENT (in) :: kjpindex !! Domain size (unitless)
269	REAL(r_std), DIMENSION (kjpindex,2), INTENT(in) :: lalo !! coordinates
270	INTEGER(i_std),DIMENSION (kjpindex,NbNeighb), INTENT(in):: neighbours !! Neighbouring land grid cell
271	REAL(r_std), DIMENSION (kjpindex,2), INTENT(in) :: resolution !! Size of grid in x and y direction (m)
272	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: contfrac !! Fraction of land in each grid box
273	INTEGER(i_std), INTENT (in) :: rest_id !! Restart file identifier (unitless)
274	REAL(r_std), DIMENSION (kjpindex,nvm), INTENT (in) :: veget_max !! Fraction of vegetation type
275	REAL(r_std), DIMENSION (kjpindex,nslm), INTENT (in) :: shumdiag_perma !! Soil saturation degree on the diagnostic axis (0-1, unitless)
276	REAL(r_std), DIMENSION (kjpindex), INTENT (in) :: snow !! Snow mass @tex ($kg$) @endtex.
277	REAL(r_std), DIMENSION(kjpindex), INTENT (in) :: thawed_humidity !! specified humidity of thawed soil
278	REAL(r_std), DIMENSION(kjpindex,ndeep,nvm), INTENT (in) :: soilc_total !! total soil carbon for use in thermal calcs
279	REAL(r_std), DIMENSION (kjpindex), INTENT (in) :: temp_sol_new !! Surface temperature at the present time-step,
280	REAL(r_std), DIMENSION (kjpindex, nvm), INTENT (in) :: temp_sol_new_pft !! Surface temperature at the present time-step,
281
282	!! 0.2 Output variables
283
284
285	!! 0.3 Modified variables
286	REAL(r_std), DIMENSION(kjpindex), INTENT (inout) :: organic_layer_thick!! how deep is the organic soil?
287	REAL(r_std),DIMENSION (kjpindex,nslm), INTENT (out) :: stempdiag !! temperature profile on the levels in hydrol(K)
288	REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: soilflx !! apparent soil heat flux @tex ($W m^{-2}$) @endtex
289	REAL(r_std),DIMENSION (kjpindex, nvm), INTENT (inout) :: soilflx_pft !! apparent soil heat flux @tex ($W m^{-2}$) @endtex
290	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: soilcap !! apparent surface heat capacity considering snow and soil surface (J m-2 K-1)
291	REAL(r_std),DIMENSION (kjpindex, nvm), INTENT (inout) :: soilcap_pft !! apparent soil heat flux @tex ($W m^{-2}$) @endtex
292	REAL(r_std),DIMENSION (kjpindex),INTENT(out) :: gtemp !! First soil layer temperature
293
294	!! 0.3 Modified variables
295	REAL(r_std), DIMENSION (kjpindex), INTENT(inout) :: lambda_snow !! Coefficient of the linear extrapolation of surface temperature
296	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT (inout):: cgrnd_snow !! Integration coefficient for snow numerical scheme
297	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT (inout):: dgrnd_snow !! Integration coefficient for snow numerical scheme
298
299	!! 0.4 Local variables
300	REAL(r_std),DIMENSION (kjpindex,ngrnd,nvm) :: reftemp_3d
301	INTEGER(i_std) :: ier, i, m
302	INTEGER(i_std) :: jv, jg
303
304	CHARACTER(LEN=80) :: var_name !! To store variables names for I/O
305	REAL(r_std), DIMENSION (kjpindex,nvm) :: veget_max_bg
306
307	LOGICAL, SAVE :: ok_zimov
308	!$OMP THREADPRIVATE(ok_zimov)
309	REAL(r_std),DIMENSION (kjpindex,ngrnd) :: temp !! buffer
310	REAL(r_std),DIMENSION (kjpindex,ngrnd-1) :: temp1 !! buffer
311	REAL(r_std),DIMENSION (kjpindex) :: temp2 !! buffer
312	LOGICAL :: calculate_coef !! Local flag to initialize variables by call to thermosoil_coef
313	!_ ================================================================================================================================
314
315	IF (printlev >= 3) WRITE (numout,*) 'Start thermosoil_initialize '
316
317	!! 1. Initialisation
318
319	!
320	! !! Flag to consider bedrock at deeper layers
321	! !! It affects heat capacity and thermal conductivity (energy balance).
322	!
323	!Config Key = BEDROCK_FLAG
324	!Config Desc = Flag to consider bedrock at deeper layers.
325	!Config If = ok_freeze_thermix
326	!Config Def = 0
327	!Config Help = 0, no, 1, yes.
328	!Config Units = [FLAG]
329	brk_flag = 0
330	CALL getin_p('BEDROCK_FLAG', brk_flag)
331
332	!Config Key = OK_WETDIAGLONG
333	!Config Desc = Long-term soil humidity (for permafrost)
334	!Config If = ok_freeze_thermix
335	!Config Def =
336	!Config Help =
337	!Config Units = [FLAG]
338	ok_shum_ngrnd_permalong = .FALSE.
339	CALL getin_p ('OK_WETDIAGLONG',ok_shum_ngrnd_permalong)
340
341	IF (ok_freeze_thermix .AND. ok_pc) THEN
342	ok_shum_ngrnd_permalong = .TRUE.
343	ENDIF
344
345	!Config Key = satsoil
346	!Config Desc =
347	!Config If =
348	!Config Def =
349	!Config Help =
350	!Config Units = [FLAG]
351	CALL getin_p('satsoil', satsoil)
352	IF (ok_freeze_thermix .AND. ok_pc) THEN
353	use_toporganiclayer_tempdiff = .false.
354	!Config Key = USE_TOPORGANICLAYER_TEMPDIFF
355	!Config Desc =
356	!Config If =
357	!Config Def =
358	!Config Help =
359	!Config Units = [FLAG]
360	CALL getin_p('USE_TOPORGANICLAYER_TEMPDIFF',use_toporganiclayer_tempdiff)
361
362	use_soilc_tempdiff = .false.
363	!Config Key = USE_SOILC_TEMPDIFF
364	!Config Desc =
365	!Config If =
366	!Config Def =
367	!Config Help =
368	!Config Units = [FLAG]
369	CALL getin_p('USE_SOILC_TEMPDIFF', use_soilc_tempdiff)
370	IF (use_toporganiclayer_tempdiff .AND. use_soilc_tempdiff) THEN
371	WRITE(,) 'warning: thermosoil_getdiff: cant have both use_toporganiclayer_tempdiff and'
372	WRITE(,) 'use_soilc_tempdiff set to .true.. using only use_soilc_tempdiff.'
373	use_toporganiclayer_tempdiff = .FALSE.
374	ENDIF
375
376	IF (use_soilc_tempdiff) THEN
377	use_refSOC = .TRUE.
378	!Config Key = use_refSOC
379	!Config Desc =
380	!Config If =
381	!Config Def =
382	!Config Help =
383	!Config Units = [FLAG]
384	CALL getin_p('use_refSOC',use_refSOC)
385	ENDIF
386	ENDIF
387
388	!Config Key = USE_SOILC_METHOD
389	!Config Desc = Flag to control the way to average thermal conductivity of mineral soil and organic soil
390	!Config If =
391	!Config Def = 1
392	!Config Help = 1=arithmetic mean ; 2=geometric mean
393	!Config Units = [FLAG]
394	use_soilc_method = SOILC_METHOD_ARITHMETIC
395	CALL getin_p('USE_SOILC_METHOD', use_soilc_method)
396
397	!Config Key = SNOW_COND_METHOD
398	!Config Desc = Flag to choose the way to calculate snow thermal conductivity
399	!Config If =
400	!Config Def = 1
401	!Config Help = 1: original 2: follows Decharme et al 2016
402	!Config Units = [FLAG]
403	snow_cond_method = SNOW_COND_METHOD_DEFAULT
404	CALL getin_p('SNOW_COND_METHOD', snow_cond_method)
405
406	!! 2. Arrays allocations
407	ALLOCATE (reftemp(kjpindex,ngrnd),stat=ier)
408	IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of reftemp','','')
409	reftemp(:,:) = 0
410
411	ALLOCATE (refSOC(kjpindex,ngrnd),stat=ier)
412	IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of refSOC','','')
413	refSOC(:,:) = 0
414
415	ALLOCATE (ptn(kjpindex,ngrnd,nvm),stat=ier)
416	IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of ptn','','')
417	ptn(:,:,:) = 0
418
419	ALLOCATE (ptn_pftmean(kjpindex,ngrnd),stat=ier)
420	IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of ptn_pftmean','','')
421	ptn_pftmean(:,:) = 0
422
423	ALLOCATE (dz1(ngrnd),stat=ier)
424	IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of dz1','','')
425	dz1(:) = 0
426
427	ALLOCATE (z1(kjpindex),stat=ier)
428	IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of z1','','')
429	z1(:) = 0
430
431	ALLOCATE (cgrnd(kjpindex,ngrnd-1,nvm),stat=ier)
432	IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of cgrnd','','')
433	cgrnd(:,:,:) = 0
434
435	ALLOCATE (dgrnd(kjpindex,ngrnd-1,nvm),stat=ier)
436	IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of dgrnd','','')
437	dgrnd(:,:,:) = 0
438
439	ALLOCATE (pcapa(kjpindex,ngrnd,nvm),stat=ier)
440	IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of pcapa','','')
441	pcapa(:,:,:) = 0
442
443	ALLOCATE (pkappa(kjpindex,ngrnd,nvm),stat=ier)
444	IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of pkappa','','')
445	pkappa(:,:,:) = 0
446
447	ALLOCATE (pcapa_snow(kjpindex,nsnow),stat=ier)
448	IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of pcapa_snow','','')
449	pcapa_snow(:,:) = 0
450
451	ALLOCATE (pkappa_snow(kjpindex,nsnow),stat=ier)
452	IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of pkappa_snow','','')
453	pkappa_snow(:,:) = 0
454
455	ALLOCATE (surfheat_incr(kjpindex),stat=ier)
456	IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of surfheat_incr','','')
457
458	ALLOCATE (coldcont_incr(kjpindex),stat=ier)
459	IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of coldcont_incr','','')
460
461	ALLOCATE (pcapa_en(kjpindex,ngrnd,nvm),stat=ier)
462	IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of pcapa_en','','')
463	! Initialization to zero used at first time step in thermosoil_energy_diag, only for diagnostic variables coldcont_incr and surfheat_incr
464	pcapa_en = 0.
465
466	ALLOCATE (ptn_beg(kjpindex,ngrnd,nvm),stat=ier)
467	IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of ptn_beg','','')
468
469	ALLOCATE (temp_sol_beg(kjpindex),stat=ier)
470	IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of temp_sol_beg','','')
471
472	ALLOCATE (shum_ngrnd_perma(kjpindex,ngrnd,nvm),stat=ier)
473	IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of shum_ngrnd_perma','','')
474	shum_ngrnd_perma(:,:,:)=val_exp
475
476	ALLOCATE (shum_ngrnd_permalong(kjpindex,ngrnd,nvm),stat=ier)
477	IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of shum_ngrnd_permalong','','')
478	shum_ngrnd_permalong = val_exp
479
480	ALLOCATE (profil_froz(kjpindex,ngrnd,nvm),stat=ier)
481	IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of profil_froz','','')
482
483	IF (ok_freeze_thermix) THEN
484	ALLOCATE (pcappa_supp(kjpindex,ngrnd,nvm),stat=ier)
485	IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of pcapa_supp','','')
486	ELSE
487	ALLOCATE(pcappa_supp(1,1,1),stat=ier)
488	IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of pcapa_supp','','')
489	END IF
490
491	IF (ok_Ecorr) THEN
492	ALLOCATE (e_soil_lat(kjpindex,nvm),stat=ier)
493	IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of e_soil_lat','','')
494	END IF
495
496	ALLOCATE (dz5(ngrnd),stat=ier)
497	IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of dz5','','')
498
499	ALLOCATE (mc_layt(kjpindex,ngrnd),stat=ier)
500	IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of mc_layt','','')
501
502	ALLOCATE (mcl_layt(kjpindex,ngrnd),stat=ier)
503	IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of mcl_layt','','')
504
505	ALLOCATE (tmc_layt(kjpindex,ngrnd),stat=ier)
506	IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of tmc_layt','','')
507
508	ALLOCATE (mc_layt_pft(kjpindex,ngrnd,nvm),stat=ier)
509	IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of mc_layt_pft','','')
510
511	ALLOCATE (mcl_layt_pft(kjpindex,ngrnd,nvm),stat=ier)
512	IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of mcl_layt_pft','','')
513
514	ALLOCATE (tmc_layt_pft(kjpindex,ngrnd,nvm),stat=ier)
515	IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of tmc_layt_pft','','')
516
517	ALLOCATE (mcs(nscm),stat=ier)
518	IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of mcs','','')
519
520	ALLOCATE (SMCMAX(nscm),stat=ier)
521	IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of SMCMAX','','')
522
523	ALLOCATE (QZ(nscm),stat=ier)
524	IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of QZ','','')
525
526	ALLOCATE (so_capa_dry_ns(nscm),stat=ier)
527	IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of so_capa_dry_ns','','')
528
529
530	!! Soil texture choose
531	SELECTCASE (nscm)
532	CASE (3)
533	SMCMAX(:) = SMCMAX_fao(:)
534	QZ(:) = QZ_fao(:)
535	so_capa_dry_ns(:) = so_capa_dry_ns_fao(:)
536	mcs(:) = mcs_fao(:)
537	CASE (12)
538	SMCMAX(:) = SMCMAX_usda(:)
539	QZ(:) = QZ_usda(:)
540	so_capa_dry_ns(:) = so_capa_dry_ns_usda(:)
541	mcs(:) = mcs_usda(:)
542	CASE DEFAULT
543	WRITE (numout,*) 'Unsupported soil type classification. Choose between zobler, fao and usda according to the map'
544	STOP 'thermosoil_initialize'
545	ENDSELECT
546
547
548	!! 2. Initialize variable from restart file or with default values
549
550	!! Reads restart files for soil temperatures only. If no restart file is
551	!! found, the initial soil temperature is by default set to 280K at all depths. The user
552	!! can decide to initialize soil temperatures at an other value, in which case he should set the flag THERMOSOIL_TPRO
553	!! to this specific value in the run.def.
554	IF (printlev>=3) WRITE (numout,*) ' we have to READ a restart file for THERMOSOIL variables'
555
556	CALL ioconf_setatt_p('UNITS', 'K')
557	CALL ioconf_setatt_p('LONG_NAME','Soil Temperature profile')
558	CALL restget_p (rest_id, 'ptn', nbp_glo, ngrnd, nvm, kjit, .TRUE., ptn, "gather", nbp_glo, index_g) !need to add veg dim
559	CALL restget_p (rest_id, 'refSOC', nbp_glo, ngrnd, 1, kjit, .TRUE., refSOC, "gather", nbp_glo, index_g)
560
561	! Initialize ptn if it was not found in restart file
562	IF (ALL(ptn(:,:,:) == val_exp)) THEN
563	! ptn was not found in restart file
564
565	IF (read_reftemp) THEN
566	! Read variable ptn from file
567	CALL thermosoil_read_reftempfile(kjpindex,lalo,reftemp)
568	DO jv = 1,nvm
569	ptn(:,:,jv)=reftemp(:,:)
570	ENDDO ! jv = 1,nvm
571	ELSE
572	! Initialize ptn with a constant value which can be set in run.def
573
574	!Config Key = THERMOSOIL_TPRO
575	!Config Desc = Initial soil temperature profile if not found in restart
576	!Config Def = 280.
577	!Config If = OK_SECHIBA
578	!Config Help = The initial value of the temperature profile in the soil if
579	!Config its value is not found in the restart file. This should only
580	!Config be used if the model is started without a restart file. Here
581	!Config we only require one value as we will assume a constant
582	!Config throughout the column.
583	!Config Units = Kelvin [K]
584	CALL setvar_p (ptn, val_exp,'THERMOSOIL_TPRO',280._r_std)
585	ENDIF
586	ENDIF
587
588	IF (ALL(refSOC(:,:) == val_exp)) THEN
589	IF (use_refSOC) THEN
590	CALL read_refSOCfile(kjpindex,lalo,neighbours, resolution, contfrac)
591	ENDIF
592	ENDIF
593
594	! Initialize ptn_beg (variable needed in thermosoil_readadjust called from thermosoil_coef)
595	ptn_beg(:,:,:) = ptn(:,:,:)
596
597	! Initialize temp_sol_beg with values from previous time-step
598	temp_sol_beg(:) = temp_sol_new(:)
599
600	IF (ok_shum_ngrnd_permalong) THEN
601	shum_ngrnd_permalong(:,:,:) = val_exp
602	CALL ioconf_setatt_p('UNITS', '-')
603	CALL ioconf_setatt_p('LONG_NAME','Long-term soil humidity')
604	CALL restget_p (rest_id, 'shum_ngrnd_prmlng', nbp_glo, ngrnd, nvm, kjit, .TRUE.,shum_ngrnd_permalong, "gather", nbp_glo, index_g)
605	IF ( ALL(ABS(shum_ngrnd_permalong(:,:,:)-val_exp).LT.EPSILON(val_exp)) ) THEN
606	shum_ngrnd_permalong(:,:,:) = 1.
607	ENDIF
608	ENDIF
609
610	shum_ngrnd_perma(:,:,:) = val_exp
611	CALL ioconf_setatt_p('UNITS', '-')
612	CALL ioconf_setatt_p('LONG_NAME','soil humidity')
613	CALL restget_p (rest_id, 'shum_ngrnd_perma', nbp_glo, ngrnd, nvm, kjit, .TRUE.,shum_ngrnd_perma, "gather", nbp_glo, index_g)
614
615	IF (printlev>=3) WRITE(numout,*) 'before comparing epsilon'
616	IF ( ALL(ABS(shum_ngrnd_perma(:,:,:)-val_exp).LT.EPSILON(val_exp)) ) THEN
617	shum_ngrnd_perma(:,:,:) = 1.
618	ENDIF
619
620	IF (ok_Ecorr) THEN
621	CALL restget_p (rest_id, 'e_soil_lat', nbp_glo, nvm, 1, kjit,.TRUE.,e_soil_lat, "gather", nbp_glo, index_g)
622	CALL setvar_p (e_soil_lat, val_exp,'NO_KEYWORD',zero)
623	ENDIF
624
625	IF (printlev>=3) WRITE (numout,*) ' thermosoil_init done '
626
627	veget_max_bg(:,2:nvm) = veget_max(:,2:nvm)
628	veget_max_bg(:,1) = MAX((un - SUM(veget_max(:,2:nvm), 2)), zero)
629
630	! Initialize ptn_pftmean
631	DO m=1,nvm
632	DO jg = 1, ngrnd
633	ptn_pftmean(:,jg) = ptn_pftmean(:,jg) + ptn(:,jg,m) * veget_max_bg(:,m)
634	ENDDO ! jg = 1, ngrnd
635	ENDDO ! m=1,nvm
636
637	CALL getin_p('OK_ZIMOV',ok_zimov)
638
639	!! 1.1. Allocate and initialize soil temperatures variables
640	!! by reading restart files or using default values.
641	! CALL thermosoil_init (kjit, ldrestart_read, kjpindex, index, lalo, rest_id, &
642	! & snowdz)
643
644	!! 1.2.Computes physical constants and arrays; initializes soil thermal properties; produces the first stempdiag
645	!! Computes some physical constants and arrays depending on the soil vertical discretization
646	!! (lskin, cstgrnd, zz, zlt, dz1, dz2); get the vertical humidity onto the thermal levels, and
647	!! initializes soil thermal properties (pkappa, pcapa); produces the first temperature diagnostic stempdiag.
648
649	CALL thermosoil_var_init (kjpindex, &
650	& shumdiag_perma, stempdiag, profil_froz, snowdz, &
651	& thawed_humidity, organic_layer_thick, soilc_total, veget_max_bg, njsc, &
652	& mc_layh, mcl_layh, tmc_layh, mc_layh_pft, mcl_layh_pft, tmc_layh_pft, &
653	& snowrho, snowtemp, pb)
654	!
655	!! 1.3. Computes cgrd, dgrd, soilflx and soilcap coefficients from restart values or initialisation values.
656	! computes cgrd and dgrd coefficient from previous time step (restart)
657	!
658	CALL thermosoil_coef (&
659	kjpindex, temp_sol_new, temp_sol_new_pft, snow, &
660	soilcap, soilcap_pft, soilflx, soilflx_pft, njsc, &
661	cgrnd, dgrnd, profil_froz, pcappa_supp, &
662	organic_layer_thick, soilc_total, veget_max, snowdz, &
663	snowrho, snowtemp, pb, &
664	frac_snow_veg, frac_snow_nobio, totfrac_nobio, &
665	lambda_snow, cgrnd_snow, dgrnd_snow)
666
667	! make vegetation masks so that we don't bother to calculated pfts on
668	! gridcells where they don's exist
669	veget_max_bg(:,2:nvm) = veget_max(:,2:nvm)
670	veget_max_bg(:,1) = MAX((un - SUM(veget_max(:,2:nvm), 2)), zero)
671
672	! Read gtemp from restart file
673	CALL restget_p (rest_id, 'gtemp', nbp_glo, 1, 1, kjit, .TRUE., &
674	gtemp, "gather", nbp_glo, index_g)
675	CALL setvar_p (gtemp, val_exp,'NO_KEYWORD',zero)
676
677
678	! Read variables calculated in thermosoil_coef from restart file
679	! If the variables were not found in the restart file, the logical
680	! calculate_coef will be true and thermosoil_coef will be called further below.
681	! These variables need to be in the restart file to avoid a time shift that
682	! would be done using thermosoil_coef at this stage.
683	calculate_coef=.FALSE.
684	CALL ioconf_setatt_p('UNITS', 'J m-2 K-1')
685	CALL ioconf_setatt_p('LONG_NAME','Apparent surface heat capacity')
686	CALL restget_p (rest_id, 'soilcap', nbp_glo, 1, 1, kjit, .TRUE., &
687	soilcap, "gather", nbp_glo, index_g)
688	IF (ALL(soilcap(:)==val_exp)) calculate_coef=.TRUE.
689
690	CALL restget_p (rest_id, 'soilcap_pft', nbp_glo, nvm, 1, kjit, .TRUE., &
691	soilcap_pft, "gather", nbp_glo, index_g)
692	IF (ALL(soilcap_pft(:,:) == val_exp)) calculate_coef=.TRUE.
693
694	CALL ioconf_setatt_p('UNITS', 'W m-2')
695	CALL ioconf_setatt_p('LONG_NAME','Apparent soil heat flux')
696	CALL restget_p (rest_id, 'soilflx', nbp_glo, 1, 1, kjit, .TRUE., &
697	soilflx, "gather", nbp_glo, index_g)
698	IF (ALL(soilflx(:)==val_exp)) calculate_coef=.TRUE.
699
700	CALL restget_p (rest_id, 'soilflx_pft', nbp_glo, nvm, 1, kjit, .TRUE., &
701	soilflx_pft, "gather", nbp_glo, index_g)
702	IF (ALL(soilflx_pft(:,:) == val_exp)) calculate_coef=.TRUE.
703
704	CALL ioconf_setatt_p('UNITS', 'J m-2 K-1')
705	CALL ioconf_setatt_p('LONG_NAME','Integration coefficient for the numerical scheme')
706	CALL restget_p (rest_id, 'cgrnd', nbp_glo, ngrnd-1, nvm, kjit, .TRUE., &
707	cgrnd, "gather", nbp_glo, index_g)
708	IF (ALL(cgrnd(:,:,:)==val_exp)) calculate_coef=.TRUE.
709
710	CALL ioconf_setatt_p('UNITS', '')
711	CALL ioconf_setatt_p('LONG_NAME','Integration coefficient for the numerical scheme')
712	CALL restget_p (rest_id, 'dgrnd', nbp_glo, ngrnd-1, nvm, kjit, .TRUE., &
713	dgrnd, "gather", nbp_glo, index_g)
714	IF (ALL(dgrnd(:,:,:)==val_exp)) calculate_coef=.TRUE.
715
716	CALL ioconf_setatt_p('UNITS', '')
717	CALL ioconf_setatt_p('LONG_NAME','Integration coefficient for the numerical scheme')
718	CALL restget_p (rest_id, 'cgrnd_snow', nbp_glo, nsnow, 1, kjit, .TRUE., &
719	cgrnd_snow, "gather", nbp_glo, index_g)
720	IF (ALL(cgrnd_snow(:,:)==val_exp)) calculate_coef=.TRUE.
721
722	CALL ioconf_setatt_p('UNITS', '')
723	CALL ioconf_setatt_p('LONG_NAME','Integration coefficient for the numerical scheme')
724	CALL restget_p (rest_id, 'dgrnd_snow', nbp_glo, nsnow, 1, kjit, .TRUE., &
725	dgrnd_snow, "gather", nbp_glo, index_g)
726	IF (ALL(dgrnd_snow(:,:)==val_exp)) calculate_coef=.TRUE.
727
728	CALL ioconf_setatt_p('UNITS', '')
729	CALL ioconf_setatt_p('LONG_NAME','Coefficient of the linear extrapolation of surface temperature')
730	CALL restget_p (rest_id, 'lambda_snow', nbp_glo, 1, 1, kjit, .TRUE., &
731	lambda_snow, "gather", nbp_glo, index_g)
732	IF (ALL(lambda_snow(:)==val_exp)) calculate_coef=.TRUE.
733
734	!! 2.2.Computes physical constants and arrays; initializes soil thermal properties; produces the first stempdiag
735	!! Computes some physical constants and arrays depending on the soil vertical discretization
736	!! (lskin, cstgrnd, zz, zlt, dz1, dz2); get the vertical humidity onto the thermal levels
737	CALL thermosoil_var_init (kjpindex, &
738	shumdiag_perma, stempdiag, profil_froz, snowdz, &
739	thawed_humidity, organic_layer_thick, soilc_total, veget_max, njsc, &
740	mc_layh, mcl_layh, tmc_layh, mc_layh_pft, mcl_layh_pft, tmc_layh_pft, &
741	snowrho, snowtemp, pb)
742
743	!! 2.3. Computes cgrnd, dgrnd, soilflx and soilcap coefficients only if they were not found in restart file.
744	IF (calculate_coef) THEN
745	! Interpolate variables needed by thermosoil_coef to the thermal levels
746	CALL thermosoil_humlev(kjpindex, shumdiag_perma, thawed_humidity, mc_layh, mcl_layh, tmc_layh, &
747	mc_layh_pft, mcl_layh_pft, tmc_layh_pft, &
748	mc_layt, mcl_layt, tmc_layt, & ! out
749	mc_layt_pft, mcl_layt_pft, tmc_layt_pft, & ! out
750	shum_ngrnd_perma ) ! out
751
752	IF (printlev>=3) WRITE (numout,*) 'thermosoil_coef will be called in the intialization phase'
753
754	CALL thermosoil_coef (&
755	kjpindex, temp_sol_new, temp_sol_new_pft, snow, &
756	soilcap, soilcap_pft, soilflx, soilflx_pft, njsc, &
757	cgrnd, dgrnd, profil_froz, pcappa_supp, &
758	organic_layer_thick, soilc_total, veget_max, snowdz, &
759	snowrho, snowtemp, pb, &
760	frac_snow_veg,frac_snow_nobio,totfrac_nobio, &
761	lambda_snow, cgrnd_snow, dgrnd_snow)
762
763	END IF
764
765	END SUBROUTINE thermosoil_initialize
766
767
768	!! ================================================================================================================================
769	!! SUBROUTINE : thermosoil_main
770	!!
771	!>\BRIEF Thermosoil_main computes the soil thermal properties and dynamics, ie solves
772	!! the heat diffusion equation within the soil.
773	!!
774	!! DESCRIPTION : The resolution of the soil heat diffusion equation
775	!! relies on a numerical finite-difference implicit scheme
776	!! fully described in the reference and in the header of the thermosoil module.
777	!! - The dependency of the previous timestep hidden in the
778	!! integration coefficients cgrnd and dgrnd (EQ1), calculated in thermosoil_coef, and
779	!! called at the end of the routine to prepare for the next timestep.
780	!! - The effective computation of the new soil temperatures is performed in thermosoil_profile.
781	!!
782	!! - thermosoil_coef calculates the coefficients for the numerical scheme for the very first iteration of thermosoil;
783	!! after that, thermosoil_coef is called only at the end of the module to calculate the coefficients for the next timestep.
784	!! - thermosoil_profile solves the numerical scheme.\n
785	!!
786	!! - Flags : one unique flag : THERMOSOIL_TPRO (to be set to the desired initial soil in-depth temperature in K; by default 280K)
787	!!
788	!! RECENT CHANGE(S) : Change vertical discretization (consistent with hydrology layers) and soil thermal properties (taking into account soil texture effects).
789	!!
790	!! MAIN OUTPUT VARIABLE(S): vertically discretized soil temperatures ptn, soil
791	!! thermal properties (pcapa, pkappa), apparent surface heat capacity (soilcap)
792	!! and heat flux (soilflx) to be used in enerbil at the next timestep to solve
793	!! the surface energy balance.
794	!!
795	!! REFERENCE(S) :
796	!! - Hourdin, F. (1992). Study and numerical simulation of the general circulation of planetary atmospheres,
797	!! Ph.D. thesis, Paris VII University. Remark: the part of F. Hourdin' s PhD thesis relative to the thermal
798	!! integration scheme has been scanned and is provided along with the documentation, with name :
799	!! Hourdin_1992_PhD_thermal_scheme.pdf
800	!!
801	!! FLOWCHART :
802	!! \latexonly
803	!! \includegraphics[scale = 1]{thermosoil_flowchart.png}
804	!! \endlatexonly
805	!!
806	!! \n
807	!_ ================================================================================================================================
808
809	SUBROUTINE thermosoil_main (kjit, kjpindex, &
810	index, indexgrnd, &
811	temp_sol_new, temp_sol_new_pft, snow, soilcap, soilcap_pft, soilflx, soilflx_pft, &
812	shumdiag_perma, stempdiag, ptnlev1, hist_id, hist2_id, &
813	snowdz,snowrho, snowtemp,gtemp, pb, &
814	mc_layh, mcl_layh, tmc_layh, mc_layh_pft, mcl_layh_pft, tmc_layh_pft, njsc, &
815	thawed_humidity, organic_layer_thick, heat_Zimov, deeptemp_prof, deephum_prof,&
816	soilc_total, veget_max, &
817	frac_snow_veg,frac_snow_nobio,totfrac_nobio, temp_sol_add, &
818	lambda_snow, cgrnd_snow, dgrnd_snow)
819
820
821	!! 0. Variable and parameter declaration
822
823	!! 0.1 Input variables
824
825	INTEGER(i_std), INTENT(in) :: kjit !! Time step number (unitless)
826	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (unitless)
827	INTEGER(i_std),INTENT (in) :: hist_id !! Restart_ history file identifier
828	!! (unitless)
829	INTEGER(i_std),INTENT (in) :: hist2_id !! history file 2 identifier (unitless)
830	INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: index !! Indeces of the points on the map (unitless)
831	INTEGER(i_std),DIMENSION (kjpindex*ngrnd), INTENT (in):: indexgrnd !! Indeces of the points on the 3D map (vertical
832	!! dimension towards the ground) (unitless)
833	REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: temp_sol_new !! Surface temperature at the present time-step,
834	!! temp_sol_new is only modified for the case ok_explicitsnow
835	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: temp_sol_new_pft !! Surface temperature at the present time-step,
836	!! Ts @tex ($K$) @endtex
837	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: snow !! Snow mass @tex ($kg$) @endtex.
838	!! Caveat: when there is snow on the
839	!! ground, the snow is integrated into the soil for
840	!! the calculation of the thermal dynamics. It means
841	!! that the uppermost soil layers can completely or
842	!! partially consist in snow. In the second case, zx1
843	!! and zx2 are the fraction of the soil layer
844	!! consisting in snow and 'normal' soil, respectively
845	!! This is calculated in thermosoil_coef.
846	REAL(r_std),DIMENSION (kjpindex,nslm), INTENT (in) :: shumdiag_perma !! Soil saturation degree (0-1, unitless)
847	REAL(r_std), DIMENSION(kjpindex), INTENT (in) :: thawed_humidity !! specified humidity of thawed soil
848	REAL(r_std), DIMENSION(kjpindex,ndeep,nvm), INTENT (in) :: heat_Zimov !! heating associated with decomposition
849	REAL(r_std), DIMENSION(kjpindex,ndeep,nvm), INTENT (in) :: soilc_total !! total soil carbon for use in thermal calcs
850	REAL(r_std), DIMENSION (kjpindex,nvm), INTENT (in) :: veget_max !! Fraction of vegetation type
851	REAL(r_std),DIMENSION (kjpindex,nsnow),INTENT (in) :: snowdz !! Snow depth
852	REAL(r_std), DIMENSION(kjpindex),INTENT (in) :: organic_layer_thick !! how deep is the organic soil?
853	REAL(r_std),DIMENSION (kjpindex,nsnow),INTENT (in) :: snowrho !! Snow density
854	REAL(r_std),DIMENSION (kjpindex,nsnow),INTENT (inout) :: snowtemp !! Snow temperature (K)
855	REAL(r_std), DIMENSION (kjpindex),INTENT (in) :: pb !! Surface presure (hPa)
856	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)
857	REAL(r_std),DIMENSION (kjpindex,nslm), INTENT (in) :: mcl_layh !! Volumetric soil moisture content for each layer in hydrol at nodes(liquid) (m3/m3)
858	REAL(r_std),DIMENSION (kjpindex,nslm), INTENT (in) :: tmc_layh !! Total soil moisture content for each layer in hydrol(liquid + ice) (mm)
859	REAL(r_std),DIMENSION (kjpindex,nslm,nvm), INTENT (in) :: mc_layh_pft !! Volumetric soil moisture content for each layer in hydrol at nodes(liquid + ice) (m3/m3)
860	REAL(r_std),DIMENSION (kjpindex,nslm,nvm), INTENT (in) :: mcl_layh_pft !! Volumetric soil moisture content for each layer in hydrol at nodes(liquid) (m3/m3)
861	REAL(r_std),DIMENSION (kjpindex,nslm,nvm), INTENT (in) :: tmc_layh_pft !! Total soil moisture content for each layer in hydrol(liquid + ice) (mm)
862	INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: njsc !! Index of the dominant soil textural class in the grid cell (1-nscm, unitless)
863	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: frac_snow_veg !! Snow cover fraction on vegeted area
864	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT(in) :: frac_snow_nobio !! Snow cover fraction on non-vegeted area
865	REAL(r_std),DIMENSION (kjpindex),INTENT(in) :: totfrac_nobio !! Total fraction of continental ice+lakes+cities+...
866	!!(unitless,0-1)
867	REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: temp_sol_add !! additional surface temperature due to the melt of first layer
868	!! at the present time-step @tex ($K$) @endtex
869
870	!! 0.2 Output variables
871
872	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: ptnlev1 !! 1st level soil temperature
873	REAL(r_std), DIMENSION (kjpindex,ndeep,nvm), INTENT (out) :: deephum_prof !! moisture on a deep thermodynamic profile for permafrost calcs
874	REAL(r_std), DIMENSION (kjpindex,ndeep,nvm), INTENT (out) :: deeptemp_prof!! temp on a deep thermodynamic profile for permafrost calcs
875	REAL(r_std),DIMENSION (kjpindex),INTENT(out) :: gtemp !! First soil layer temperature
876
877	!! 0.3 Modified variables
878
879	REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: soilcap !! apparent surface heat capacity considering snow and soil surface
880	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (inout) :: soilcap_pft !! apparent surface heat capacity
881	!! @tex ($J m^{-2} K^{-1}$) @endtex
882	REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: soilflx !! apparent soil heat flux considering snow and soil surface
883	!! @tex ($W m^{-2}$) @endtex
884	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (inout) :: soilflx_pft !! apparent soil heat flux @tex ($W m^{-2}$) @endtex
885	!! , positive
886	!! towards the soil, writen as Qg (ground heat flux)
887	!! in the history files, and computed at the end of
888	!! thermosoil for the calculation of Ts in enerbil,
889	!! see EQ3.
890	REAL(r_std),DIMENSION (kjpindex,nslm), INTENT (out) :: stempdiag !! temperature profile @tex ($K$) @endtex
891	REAL(r_std),DIMENSION (kjpindex), INTENT(inout) :: lambda_snow !! Coefficient of the linear extrapolation of surface temperature
892	REAL(r_std),DIMENSION (kjpindex,nsnow), INTENT (inout):: cgrnd_snow !! Integration coefficient for snow numerical scheme
893	REAL(r_std),DIMENSION (kjpindex,nsnow), INTENT (inout):: dgrnd_snow !! Integration coefficient for snow numerical scheme
894
895	!! 0.4 Local variables
896
897	REAL(r_std), DIMENSION (kjpindex,nvm) :: veget_max_bg !! Fraction of vegetation type
898	LOGICAL, SAVE :: ok_zimov
899	!$OMP THREADPRIVATE(ok_zimov)
900	REAL(r_std),DIMENSION (kjpindex,ngrnd) :: pkappa_pftmean
901	INTEGER(i_std) :: jv,ji,m,jg
902	CHARACTER(LEN=10) :: part_str !! string suffix indicating an index
903	REAL(r_std), DIMENSION (kjpindex,ngrnd) :: tmparray !! Temporary array
904	REAL(r_std),DIMENSION (kjpindex) :: snowtemp_weighted!! Snow temperature weighted by snow density, only for diag (K)
905	REAL(r_std),DIMENSION (kjpindex, nsnow) :: pkappa_snow_diag !! Only for diag, containing xios_default_val
906	REAL(r_std),DIMENSION (kjpindex, nsnow) :: pcapa_snow_diag !! Only for diag, containing xios_default_val
907	REAL(r_std),DIMENSION (kjpindex, nsnow) :: snowtemp_diag !! Only for diag, containing xios_default_val
908
909
910	!_ ================================================================================================================================
911
912	veget_max_bg(:,2:nvm) = veget_max(:,2:nvm)
913	veget_max_bg(:,1) = MAX((un - SUM(veget_max(:,2:nvm), 2)), zero)
914
915	!! 3. Put the soil wetness diagnostic on the levels of the soil temperature
916
917	!!?? this could logically be put just before the last call to
918	!!thermosoil_coef, as the results are used there...
919	CALL thermosoil_humlev(kjpindex, shumdiag_perma, thawed_humidity, mc_layh, mcl_layh, tmc_layh, &
920	mc_layh_pft, mcl_layh_pft, tmc_layh_pft, &
921	mc_layt, mcl_layt, tmc_layt, & ! out
922	mc_layt_pft, mcl_layt_pft, tmc_layt_pft, & ! out
923	shum_ngrnd_perma ) ! out
924
925	! Compute long-term soil humidity (for permafrost)
926	!
927	IF (ok_shum_ngrnd_permalong) THEN
928	CALL thermosoil_wlupdate( kjpindex, shum_ngrnd_perma, shum_ngrnd_permalong )
929	ELSE
930	shum_ngrnd_permalong(:,:,:)=shum_ngrnd_perma(:,:,:)
931	ENDIF
932	!! 4. Effective computation of the soil temperatures profile.
933	!! cgrnd and dgrnd have been calculated in thermosoil_coef at the previous time step
934	!! but they are correct for the actual time-step.
935	CALL thermosoil_profile (kjpindex, temp_sol_new, temp_sol_new_pft, ptn, stempdiag, &
936	snowtemp, veget_max, &
937	frac_snow_veg, frac_snow_nobio, totfrac_nobio, &
938	cgrnd_snow, dgrnd_snow )
939
940
941	!! 5. Call to thermosoil_energy_diag for calculation of diagnostic variables
942	CALL thermosoil_energy_diag (kjpindex, temp_sol_new, soilcap, veget_max_bg)
943	ptn_pftmean(:,:) = zero
944	DO m=1,nvm
945	DO jg = 1, ngrnd
946	ptn_pftmean(:,jg) = ptn_pftmean(:,jg) + ptn(:,jg,m) * veget_max_bg(:,m)
947	ENDDO ! jg = 1, ngrnd
948	ENDDO ! m=1,nvm
949
950	!! Save ptn at current stage, to be used in thermosoil_readjust
951	ptn_beg(:,:,:) = ptn(:,:,:)
952
953	!! 6. Writing the history files according to the ALMA standards (or not..)
954
955	! Add XIOS default value where no snow
956	DO ji=1,kjpindex
957	IF (snow(ji) .GT. zero) THEN
958	pkappa_snow_diag(ji,:) = pkappa_snow(ji,:)
959	pcapa_snow_diag(ji,:) = pcapa_snow(ji,:)
960	snowtemp_diag(ji,:) = snowtemp(ji,:)
961	ELSE
962	pkappa_snow_diag(ji,:) = xios_default_val
963	pcapa_snow_diag(ji,:) = xios_default_val
964	snowtemp_diag(ji,:) = xios_default_val
965	END IF
966	END DO
967
968	IF (ok_explicitsnow) THEN
969	DO ji=1,kjpindex
970	! Use min_sechiba instead of zero to avoid problem with division by zero
971	IF (snow(ji) .GT. min_sechiba) THEN
972	snowtemp_weighted(ji) = SUM(snowtemp(ji,:)*snowrho(ji,:))/SUM(snowrho(ji,:))
973	ELSE
974	snowtemp_weighted(ji) = xios_default_val
975	END IF
976	END DO
977	CALL xios_orchidee_send_field("snowtemp_weighted",snowtemp_weighted)
978	END IF
979
980	CALL xios_orchidee_send_field("ptn",ptn)
981	CALL xios_orchidee_send_field("ptn_pftmean",ptn_pftmean)
982	CALL xios_orchidee_send_field("soilflx",soilflx)
983	CALL xios_orchidee_send_field("surfheat_incr",surfheat_incr)
984	CALL xios_orchidee_send_field("coldcont_incr",coldcont_incr)
985	CALL xios_orchidee_send_field("pkappa",pkappa)
986	CALL xios_orchidee_send_field("pkappa_snow",pkappa_snow_diag)
987	CALL xios_orchidee_send_field("pcapa",pcapa)
988	CALL xios_orchidee_send_field("pcapa_snow",pcapa_snow_diag)
989	CALL xios_orchidee_send_field("snowtemp",snowtemp_diag)
990
991	IF (ok_freeze_thermix) CALL xios_orchidee_send_field("pcappa_supp",pcappa_supp)
992	CALL xios_orchidee_send_field("shum_ngrnd_perma", shum_ngrnd_perma)
993	IF (ok_shum_ngrnd_permalong) CALL xios_orchidee_send_field("shum_ngrnd_prmlng", shum_ngrnd_permalong)
994	CALL xios_orchidee_send_field("ptn_beg", ptn_beg)
995
996	IF ( .NOT. almaoutput ) THEN
997	!!need to write with PFT dimension
998	CALL histwrite_p(hist_id, 'ptn', kjit, ptn, kjpindex*ngrnd, indexgrnd)
999	CALL histwrite_p(hist_id, 'ptn_pftmean', kjit, ptn_pftmean, kjpindex*ngrnd, indexgrnd)
1000	IF (hydrol_cwrr) THEN
1001	DO jv = 1, nvm
1002	IF (ok_freeze_thermix .AND. permafrost_veg_exists(jv)) THEN
1003	WRITE(part_str,'(I2)') jv
1004	IF (jv < 10) part_str(1:1) = '0'
1005	tmparray = pcapa(:,:,jv)
1006	CALL histwrite_p(hist_id, 'pcapa_'//part_str(1:LEN_TRIM(part_str)), &
1007	kjit, pcapa(:,:,jv), kjpindex*ngrnd, indexgrnd)
1008	tmparray = pcappa_supp(:,:,jv)
1009	CALL histwrite_p(hist_id, 'pcappa_supp_'//part_str(1:LEN_TRIM(part_str)), &
1010	kjit, pcappa_supp(:,:,jv), kjpindex*ngrnd, indexgrnd)
1011	tmparray = pkappa(:,:,jv)
1012	CALL histwrite_p(hist_id, 'pkappa_'//part_str(1:LEN_TRIM(part_str)), &
1013	kjit, pkappa(:,:,jv), kjpindex*ngrnd, indexgrnd)
1014	tmparray = ptn_beg(:,:,jv)
1015	CALL histwrite_p(hist_id, 'ptn_beg_'//part_str(1:LEN_TRIM(part_str)), &
1016	kjit, ptn_beg(:,:,jv), kjpindex*ngrnd, indexgrnd)
1017
1018	ENDIF
1019	ENDDO
1020
1021	CALL histwrite_p(hist_id, 'shum_ngrnd_perma', kjit, shum_ngrnd_perma, kjpindex*ngrnd, indexgrnd)
1022	IF (ok_shum_ngrnd_permalong) CALL histwrite_p(hist_id,'shum_ngrnd_prmlng', kjit, shum_ngrnd_permalong, kjpindex*ngrnd, indexgrnd)
1023	!CALL histwrite_p(hist_id,'profil_froz_'//part_str(1:LEN_TRIM(part_str)), &
1024	! kjit, profil_froz(:,:,jv), kjpindex*ngrnd, indexgrnd)
1025	!CALL histwrite_p(hist_id, 'shumdiag_perma', kjit, shumdiag_perma, kjpindex*nslm, indexnslm)
1026	END IF
1027	CALL histwrite_p(hist_id, 'Qg', kjit, soilflx, kjpindex, index)
1028
1029	ELSE !IF ( .NOT. almaoutput ) THEN
1030	CALL histwrite_p(hist_id, 'SoilTemp', kjit, ptn, kjpindex*ngrnd, indexgrnd)
1031	CALL histwrite_p(hist_id, 'Qg', kjit, soilflx, kjpindex, index)
1032	CALL histwrite_p(hist_id, 'DelSurfHeat', kjit, surfheat_incr, kjpindex, index)
1033	CALL histwrite_p(hist_id, 'DelColdCont', kjit, coldcont_incr, kjpindex, index)
1034	ENDIF !IF ( .NOT. almaoutput ) THEN
1035	IF ( hist2_id > 0 ) THEN
1036	IF ( .NOT. almaoutput ) THEN
1037	CALL histwrite_p(hist2_id, 'ptn_pftmean', kjit, ptn_pftmean, kjpindex*ngrnd, indexgrnd)
1038	ELSE
1039	CALL histwrite_p(hist2_id, 'SoilTemp', kjit, ptn, kjpindex*ngrnd, indexgrnd)
1040	CALL histwrite_p(hist2_id, 'Qg', kjit, soilflx, kjpindex, index)
1041	CALL histwrite_p(hist2_id, 'DelSurfHeat', kjit, surfheat_incr, kjpindex, index)
1042	CALL histwrite_p(hist2_id, 'DelColdCont', kjit, coldcont_incr, kjpindex, index)
1043	ENDIF
1044	ENDIF
1045
1046	!! 7. Considering the heat released by microbial respiration
1047	IF (ok_zimov) THEN
1048	CALL add_heat_Zimov(kjpindex, veget_max_bg, ptn, heat_zimov)
1049	END IF
1050
1051	!! 8. A last final call to thermosoil_coef
1052
1053	!! A last final call to thermosoil_coef, which calculates the different
1054	!!coefficients (cgrnd, dgrnd, soilcap, soilflx) from this time step to be
1055	!!used at the next time step, either in the surface temperature calculation
1056	!!(soilcap, soilflx) or in the soil thermal numerical scheme.
1057	!
1058	CALL thermosoil_coef (&
1059	kjpindex, temp_sol_new, temp_sol_new_pft, snow, &
1060	soilcap, soilcap_pft, soilflx, soilflx_pft, njsc, &
1061	cgrnd, dgrnd, profil_froz, pcappa_supp, &
1062	organic_layer_thick, soilc_total, veget_max_bg, snowdz, &
1063	snowrho, snowtemp, pb, &
1064	frac_snow_veg, frac_snow_nobio, totfrac_nobio, &
1065	lambda_snow, cgrnd_snow, dgrnd_snow)
1066
1067	!save some useful variables for new snow model
1068	ptn_pftmean(:,:) = zero
1069	pkappa_pftmean(:,:) = zero
1070	DO m=1,nvm
1071	DO jg = 1, ngrnd
1072	ptn_pftmean(:,jg) = ptn_pftmean(:,jg) + ptn(:,jg,m) * veget_max_bg(:,m)
1073	pkappa_pftmean(:,jg) = pkappa_pftmean(:,jg) + pkappa(:,jg,m) * veget_max_bg(:,m)
1074	END DO
1075	END DO
1076
1077	! Save variables for explicit snow model
1078
1079	DO ji=1,kjpindex
1080	gtemp(ji) = ptn_pftmean(ji,1)
1081	ENDDO
1082
1083	ptnlev1(:) = ptn_pftmean(:,1)
1084
1085	!++cdk prep updated temp and moisture fields so they can be sent to stomate
1086	!permafrost calcs
1087	deephum_prof = shum_ngrnd_permalong
1088	deeptemp_prof = ptn
1089	!--cdk
1090
1091	#ifdef STRICT_CHECK
1092	IF (ANY(deeptemp_prof > 500)) THEN
1093	CALL ipslerr_p(3, 'thermosoil_main', 'deeptemp_prof is bigger than 500 degrees kelvin', '', '')
1094	ENDIF
1095	#endif
1096
1097	!! Surface temperature is forced to zero celcius if its value is larger than melting point, only for explicit snow scheme
1098	IF (ok_explicitsnow) THEN
1099	DO ji=1,kjpindex
1100	IF (SUM(snowdz(ji,:)) .GT. 0.0) THEN
1101	IF (temp_sol_new(ji) .GE. tp_00) THEN
1102	temp_sol_new(ji) = tp_00
1103	ENDIF
1104	END IF
1105	END DO
1106	ENDIF
1107
1108	IF (printlev>=3) WRITE (numout,*) ' thermosoil_main done '
1109
1110	END SUBROUTINE thermosoil_main
1111
1112	!! =============================================================================================================================
1113	!! SUBROUTINE : thermosoil_finalize
1114	!!
1115	!>\BRIEF Write to restart file
1116	!!
1117	!! DESCRIPTION : This subroutine writes the module variables and variables calculated in thermosoil
1118	!! to restart file
1119	!! \n
1120	!_ ==============================================================================================================================
1121	SUBROUTINE thermosoil_finalize (kjit, kjpindex, rest_id, gtemp, &
1122	soilcap, soilcap_pft, soilflx, soilflx_pft, lambda_snow, cgrnd_snow, dgrnd_snow)
1123
1124	!! 0. Variable and parameter declaration
1125	!! 0.1 Input variables
1126	INTEGER(i_std), INTENT(in) :: kjit !! Time step number (unitless)
1127	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (unitless)
1128	INTEGER(i_std),INTENT (in) :: rest_id !! Restart file identifier(unitless)
1129	!! 0.2 Modified variables
1130	REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: soilcap !! apparent surface heat capacity
1131	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (inout) :: soilcap_pft !! apparent surface heat capacity
1132	!! @tex ($J m^{-2} K^{-1}$) @endtex
1133	REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: soilflx !! apparent soil heat flux @tex ($W m^{-2}$) @endtex
1134	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (inout) :: soilflx_pft !! apparent soil heat flux @tex ($W m^{-2}$) @endtex
1135	!! , positive
1136	!! towards the soil, writen as Qg (ground heat flux)
1137	!! in the history files, and computed at the end of
1138	!! thermosoil for the calculation of Ts in enerbil,
1139	!! see EQ3.
1140	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: lambda_snow !! Coefficient of the linear extrapolation of surface temperature
1141	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT (in) :: cgrnd_snow !! Integration coefficient for snow numerical scheme
1142	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT (in) :: dgrnd_snow !! Integration coefficient for snow numerical scheme
1143
1144	REAL(r_std),DIMENSION (kjpindex),INTENT(in) :: gtemp !! the first soil layer temperature
1145	!! 0.3 Local variables
1146	INTEGER(i_std) :: m
1147	CHARACTER(LEN=80) :: var_name !! To store variables names for I/O
1148
1149
1150	!_ ================================================================================================================================
1151
1152	!! 1. Write variables to restart file to be used for the next simulation
1153	IF (printlev>=3) WRITE (numout,*) 'Write restart file with THERMOSOIL variables'
1154
1155	!! 2. Prepares the restart files for the next simulation
1156
1157	IF (printlev>=3) WRITE (numout,*) ' we have to complete restart file with THERMOSOIL variables'
1158
1159	CALL restput_p (rest_id, 'ptn', nbp_glo, ngrnd, nvm, kjit, ptn, 'scatter', nbp_glo, index_g)
1160
1161	CALL restput_p (rest_id, 'refSOC', nbp_glo, ngrnd, 1, kjit, refSOC, 'scatter', nbp_glo, index_g)
1162
1163	IF (ok_shum_ngrnd_permalong) THEN
1164	CALL restput_p (rest_id, 'shum_ngrnd_prmlng', nbp_glo, ngrnd, nvm, kjit,shum_ngrnd_permalong, 'scatter', nbp_glo, index_g) !need to add veg dim
1165	END IF
1166
1167	CALL restput_p (rest_id, 'shum_ngrnd_perma', nbp_glo, ngrnd, nvm, kjit, shum_ngrnd_perma, 'scatter', nbp_glo, index_g) !need to add veg dim
1168
1169	IF (ok_Ecorr) THEN
1170	var_name = 'e_soil_lat'
1171	CALL restput_p (rest_id, var_name, nbp_glo, nvm, 1, kjit, e_soil_lat, 'scatter', nbp_glo, index_g)
1172	END IF
1173
1174	CALL restput_p (rest_id, 'cgrnd', nbp_glo, ngrnd-1, nvm, kjit, cgrnd, 'scatter', nbp_glo, index_g)
1175	CALL restput_p (rest_id, 'dgrnd', nbp_glo, ngrnd-1, nvm, kjit, dgrnd, 'scatter', nbp_glo, index_g)
1176	CALL restput_p (rest_id, 'z1', nbp_glo, 1, 1, kjit, z1, 'scatter', nbp_glo, index_g)
1177	CALL restput_p (rest_id, 'pcapa', nbp_glo, ngrnd, nvm, kjit, pcapa, 'scatter', nbp_glo, index_g)
1178	CALL restput_p (rest_id, 'pcapa_en', nbp_glo, ngrnd, nvm, kjit, pcapa_en, 'scatter', nbp_glo, index_g)
1179	CALL restput_p (rest_id, 'pkappa', nbp_glo, ngrnd, nvm, kjit, pkappa, 'scatter', nbp_glo, index_g)
1180
1181	var_name= 'temp_sol_beg'
1182	CALL restput_p(rest_id, var_name, nbp_glo, 1, 1, kjit, temp_sol_beg, 'scatter', nbp_glo, index_g)
1183
1184	CALL restput_p(rest_id, 'gtemp', nbp_glo, 1, 1, kjit, gtemp, 'scatter', nbp_glo, index_g)
1185
1186	var_name= 'soilcap'
1187	CALL restput_p(rest_id, var_name, nbp_glo, 1, 1, kjit, soilcap, 'scatter', nbp_glo, index_g)
1188	CALL restput_p(rest_id, 'soilcap_pft', nbp_glo, nvm, 1, kjit, soilcap_pft, 'scatter', nbp_glo, index_g)
1189
1190	var_name= 'soilflx'
1191	CALL restput_p(rest_id, var_name, nbp_glo, 1, 1, kjit, soilflx, 'scatter', nbp_glo, index_g)
1192	CALL restput_p(rest_id, 'soilflx_pft', nbp_glo, nvm, 1, kjit, soilflx_pft, 'scatter', nbp_glo, index_g)
1193
1194	CALL restput_p(rest_id, 'cgrnd_snow', nbp_glo, nsnow, 1, kjit, cgrnd_snow, 'scatter', nbp_glo, index_g)
1195	CALL restput_p(rest_id, 'dgrnd_snow', nbp_glo, nsnow, 1, kjit, dgrnd_snow, 'scatter', nbp_glo, index_g)
1196	CALL restput_p(rest_id, 'lambda_snow', nbp_glo, 1, 1, kjit, lambda_snow, 'scatter', nbp_glo, index_g)
1197
1198	END SUBROUTINE thermosoil_finalize
1199
1200	!! ================================================================================================================================
1201	!! SUBROUTINE : thermosoil_clear
1202	!!
1203	!>\BRIEF Deallocates the allocated arrays.
1204	!! The call of thermosoil_clear originates from sechiba_clear but the calling sequence and
1205	!! its purpose require further investigation.
1206	!!
1207	!! DESCRIPTION : None
1208	!!
1209	!! RECENT CHANGE(S) : None
1210	!!
1211	!! MAIN OUTPUT VARIABLE(S): None
1212	!!
1213	!! REFERENCE(S) : None
1214	!!
1215	!! FLOWCHART : None
1216	!! \n
1217	!_ ================================================================================================================================
1218
1219	SUBROUTINE thermosoil_clear()
1220
1221	IF ( ALLOCATED (ptn)) DEALLOCATE (ptn)
1222	IF ( ALLOCATED (ptn_pftmean)) DEALLOCATE (ptn_pftmean)
1223	IF ( ALLOCATED (z1)) DEALLOCATE (z1)
1224	IF ( ALLOCATED (cgrnd)) DEALLOCATE (cgrnd)
1225	IF ( ALLOCATED (dgrnd)) DEALLOCATE (dgrnd)
1226	IF ( ALLOCATED (pcapa)) DEALLOCATE (pcapa)
1227	IF ( ALLOCATED (pkappa)) DEALLOCATE (pkappa)
1228	IF ( ALLOCATED (pcapa_snow)) DEALLOCATE (pcapa_snow)
1229	IF ( ALLOCATED (pkappa_snow)) DEALLOCATE (pkappa_snow)
1230	IF ( ALLOCATED (pcapa_en)) DEALLOCATE (pcapa_en)
1231	IF ( ALLOCATED (ptn_beg)) DEALLOCATE (ptn_beg)
1232	IF ( ALLOCATED (temp_sol_beg)) DEALLOCATE (temp_sol_beg)
1233	IF ( ALLOCATED (surfheat_incr)) DEALLOCATE (surfheat_incr)
1234	IF ( ALLOCATED (coldcont_incr)) DEALLOCATE (coldcont_incr)
1235	IF ( ALLOCATED (shum_ngrnd_perma)) DEALLOCATE (shum_ngrnd_perma)
1236	IF ( ALLOCATED (profil_froz)) DEALLOCATE (profil_froz)
1237	IF ( ALLOCATED (shum_ngrnd_permalong)) DEALLOCATE (shum_ngrnd_permalong)
1238	IF ( ALLOCATED (mc_layt)) DEALLOCATE (mc_layt)
1239	IF ( ALLOCATED (mcl_layt)) DEALLOCATE (mcl_layt)
1240	IF ( ALLOCATED (tmc_layt)) DEALLOCATE (tmc_layt)
1241	IF ( ALLOCATED (mc_layt_pft)) DEALLOCATE (mc_layt_pft)
1242	IF ( ALLOCATED (mcl_layt_pft)) DEALLOCATE (mcl_layt_pft)
1243	IF ( ALLOCATED (tmc_layt_pft)) DEALLOCATE (tmc_layt_pft)
1244	IF ( ALLOCATED (reftemp)) DEALLOCATE (reftemp)
1245	IF ( ALLOCATED (refSOC)) DEALLOCATE (refSOC)
1246	END SUBROUTINE thermosoil_clear
1247
1248
1249	!! ================================================================================================================================
1250	!! SUBROUTINE : thermosoil_var_init
1251	!!
1252	!>\BRIEF Define and initializes the soil thermal parameters
1253	!!
1254	!! DESCRIPTION : This routine\n
1255	!! 1. Defines the parameters ruling the vertical grid of the thermal scheme (fz1, zalpha).\n
1256	!! 2. Defines the scaling coefficients for adimensional depths (lskin, cstgrnd, see explanation in the
1257	!! variables description of thermosoil_main). \n
1258	!! 3. Calculates the vertical discretization of the soil (zz, zlt, dz2) and the constants used
1259	!! in the numerical scheme and which depend only on the discretization (dz1, lambda).\n
1260	!! 4. Initializes the soil thermal parameters (capacity, conductivity) based on initial soil moisture content.\n
1261	!! 5. Produces a first temperature diagnostic based on temperature initialization.\n
1262	!!
1263	!! The scheme comprizes ngrnd=7 layers by default.
1264	!! The layer' s boundaries depths (zlt) follow a geometric series of ratio zalph=2 and first term fz1.\n
1265	!! zlt(jg)=fz1.(1-zalph^jg)/(1-zalph) \n
1266	!! The layers' boudaries depths are first calculated 'adimensionally', ie with a
1267	!! discretization adapted to EQ5. This discretization is chosen for its ability at
1268	!! reproducing a thermal signal with periods ranging from days to centuries. (see
1269	!! Hourdin, 1992). Typically, fz1 is chosen as : fz1=0.3*cstgrnd (with cstgrnd the
1270	!! adimensional attenuation depth). \n
1271	!! The factor lskin/cstgrnd is then used to go from adimensional depths to
1272	!! depths in m.\n
1273	!! zz(real)=lskin/cstgrnd*zz(adimensional)\n
1274	!! Similarly, the depths of the numerical nodes are first calculated
1275	!! adimensionally, then the conversion factor is applied.\n
1276	!! the numerical nodes (zz) are not exactly the layers' centers : their depths are calculated as follows:\n
1277	!! zz(jg)=fz1.(1-zalph^(jg-1/2))/(1-zalph)\n
1278	!! The values of zz and zlt used in the default thermal discretization are in the following table.
1279	!! \latexonly
1280	!! \includegraphics{thermosoil_var_init1.jpg}
1281	!! \endlatexonly\n
1282	!!
1283	!! RECENT CHANGE(S) : None
1284	!!
1285	!! MAIN OUTPUT VARIABLE(S) : None
1286	!!
1287	!! REFERENCE(S) :
1288	!! - Hourdin, F. (1992). Study and numerical simulation of the general circulation of
1289	!! planetary atmospheres, Ph.D. thesis, Paris VII University.
1290	!!
1291	!! FLOWCHART : None
1292	!! \n
1293	!_ ================================================================================================================================
1294
1295	SUBROUTINE thermosoil_var_init(kjpindex, &
1296	& shumdiag_perma, stempdiag, profil_froz,snowdz, &
1297	& thawed_humidity,organic_layer_thick, soilc_total, veget_max, njsc, &
1298	& mc_layh, mcl_layh, tmc_layh, mc_layh_pft, mcl_layh_pft, tmc_layh_pft, &
1299	& snowrho, snowtemp, pb )
1300
1301	!! 0. Variables and parameter declaration
1302
1303	!! 0.1 Input variables
1304
1305	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (unitless)
1306	REAL(r_std), DIMENSION (kjpindex,nslm), INTENT (in) :: shumdiag_perma !! Relative soil humidity on the diagnostic axis
1307	!! (unitless), [0,1]. (see description of the
1308	!! variables of thermosoil_main for more
1309	!! explanations)
1310
1311	REAL(r_std), DIMENSION (kjpindex,nvm), INTENT (in) :: veget_max !! Fraction of vegetation type
1312	REAL(r_std), DIMENSION (kjpindex), INTENT (in) :: thawed_humidity !! specified humidity of thawed soil
1313
1314	REAL(r_std), DIMENSION (kjpindex,nsnow),INTENT(in) :: snowdz
1315	REAL(r_std), DIMENSION (kjpindex), INTENT (in) :: organic_layer_thick !! how deep is the organic soil?
1316	REAL(r_std), DIMENSION (kjpindex,ndeep,nvm), INTENT (in) :: soilc_total !! total soil carbon for use in thermal calcs
1317	INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: njsc !! Index of the dominant soil textural class in the grid cell (1-nscm, unitless)
1318	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]
1319	REAL(r_std),DIMENSION (kjpindex,nslm), INTENT (in) :: mcl_layh !! Volumetric soil moisture content for each layer in hydrol at nodes(liquid) [m/s]
1320	REAL(r_std),DIMENSION (kjpindex,nslm), INTENT (in) :: tmc_layh !! Total soil moisture content for each layer in hydrol(liquid+ice) [mm]
1321	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(in) :: snowrho !!Snow density
1322	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(in) :: snowtemp !! Snow temperature
1323	REAL(r_std), DIMENSION (kjpindex), INTENT (in) :: pb !! Surface presure (hPa)
1324	REAL(r_std),DIMENSION (kjpindex,nslm,nvm), INTENT (in) :: mc_layh_pft !! Volumetric soil moisture content for each layer in hydrol at nodes(liquid+ice) [m/s]
1325	REAL(r_std),DIMENSION (kjpindex,nslm,nvm), INTENT (in) :: mcl_layh_pft !! Volumetric soil moisture content for each layer in hydrol at nodes(liquid) [m/s]
1326	REAL(r_std),DIMENSION (kjpindex,nslm,nvm), INTENT (in) :: tmc_layh_pft !! Total soil moisture content for each layer in hydrol(liquid+ice) [mm]
1327
1328	!! 0.2 Output variables
1329
1330	REAL(r_std), DIMENSION (kjpindex,nslm), INTENT (out) :: stempdiag !! Diagnostic temperature profile @tex ($K$)
1331	!! @endtex
1332	REAL(r_std), DIMENSION (kjpindex,ngrnd,nvm), INTENT(out) :: profil_froz
1333
1334	!! 0.3 Modified variables
1335
1336	!! 0.4 Local variables
1337
1338	INTEGER(i_std) :: ji,jv,jg !! Index (unitless)
1339
1340	!! 1. Initialization of the parameters of the vertical discretization and of the attenuation depths
1341
1342	dz5(:) = 0.0
1343	!! Calculate so_capa_ice
1344	so_capa_ice = so_capa_dry + poroscapa_icerho_ice
1345	WRITE(numout,*) 'Calculation of so_capa_ice=', so_capa_ice,' using poros=',poros,' and capa_ice=',capa_ice
1346
1347	!! 2. Get the depth of the thermal levels (numerical nodes) and the layers boundaries from vertical module
1348
1349	!! 2.2 Computing some usefull constants for the numerical scheme
1350	DO jg=1,ngrnd-1
1351	dz1(jg) = un / (znt(jg+1) - znt(jg))
1352	dz5(jg) = (zlt(jg) - znt(jg)) * dz1(jg)
1353	ENDDO
1354	lambda = znt(1) * dz1(1)
1355
1356	!! 2.3 Get the wetness profile on the thermal vertical grid from the diagnostic axis
1357	CALL thermosoil_humlev(kjpindex, shumdiag_perma, thawed_humidity, mc_layh, mcl_layh, tmc_layh, &
1358	mc_layh_pft, mcl_layh_pft, tmc_layh_pft, &
1359	mc_layt, mcl_layt, tmc_layt, & ! out
1360	mc_layt_pft, mcl_layt_pft, tmc_layt_pft, & ! out
1361	shum_ngrnd_perma ) ! out
1362	!
1363
1364	!! 2.4 Thermal conductivity at all levels
1365	if (ok_explicitsnow) then
1366	CALL thermosoil_getdiff( kjpindex, ptn, njsc, veget_max, shum_ngrnd_permalong, &
1367	profil_froz, pcappa_supp, organic_layer_thick, soilc_total, snowrho, &
1368	snowtemp, pb, mc_layt, mc_layt_pft, tmc_layt_pft, pcapa, pcapa_en, pkappa)
1369
1370	! this is for the thin snow in order to prevent the warm surface
1371	CALL thermosoil_getdiff_thinsnow (kjpindex, ptn, shum_ngrnd_permalong, snowdz, profil_froz)
1372	! else
1373	!if (ok_thermix_trunc) then
1374	! ! pour convergence avec le trunc
1375	! CALL thermosoil_getdiff_old_thermix_trunc2( kjpindex, pkappa, pcapa, pcapa_en )
1376	!else
1377	! CALL thermosoil_getdiff_old_thermix_with_snow( kjpindex, ptn, wetdiaglong, snow, pkappa, pcapa, pcapa_en,profil_froz )
1378	!endif
1379	endif
1380
1381	!! 3. Compute a first diagnostic temperature profile
1382
1383	CALL thermosoil_diaglev(kjpindex, stempdiag, veget_max)
1384
1385	IF (printlev>=3) WRITE (numout,*) ' thermosoil_var_init done '
1386
1387	END SUBROUTINE thermosoil_var_init
1388
1389
1390	!! ================================================================================================================================
1391	!! SUBROUTINE : thermosoil_coef
1392	!!
1393	!>\BRIEF Calculate soil thermal properties, integration coefficients, apparent heat flux,
1394	!! surface heat capacity,
1395	!!
1396	!! DESCRIPTION : This routine computes : \n
1397	!! 1. the soil thermal properties. \n
1398	!! 2. the integration coefficients of the thermal numerical scheme, cgrnd and dgrnd,
1399	!! which depend on the vertical grid and on soil properties, and are used at the next
1400	!! timestep.\n
1401	!! 3. the soil apparent heat flux and surface heat capacity (soilflx
1402	!! and soilcap), used by enerbil to compute the surface temperature at the next
1403	!! timestep.\n
1404	!! - The soil thermal properties depend on water content (shum_ngrnd_perma, shumdiag_perma,
1405	!! mc_layt, mcl_layt, tmc_layt), dominant soil texture(njsc), and on the presence
1406	!! of snow : snow is integrated into the soil for the thermal calculations, ie if there
1407	!! is snow on the ground, the first thermal layer(s) consist in snow, depending on the
1408	!! snow-depth. If a layer consists out of snow and soil, wheighed soil properties are
1409	!! calculated\n
1410	!! - The coefficients cgrnd and dgrnd are the integration
1411	!! coefficients for the thermal scheme \n
1412	!! T(k+1)=cgrnd(k)+dgrnd(k)*T(k) \n
1413	!! -- EQ1 -- \n
1414	!! They correspond respectively to $\beta$ and $\alpha$ from F. Hourdin\'s thesis and
1415	!! their expression can be found in this document (eq A19 and A20)
1416	!! - soilcap and soilflx are the apparent surface heat capacity and flux
1417	!! used in enerbil at the next timestep to solve the surface
1418	!! balance for Ts (EQ3); they correspond to $C_s$ and $F_s$ in F.
1419	!! Hourdin\'s PhD thesis and are expressed in eq. A30 and A31. \n
1420	!! soilcap*(Ts(t)-Ts(t-1))/dt=soilflx+otherfluxes(Ts(t)) \n
1421	!! -- EQ3 --\n
1422	!!
1423	!! RECENT CHANGE(S) : None
1424	!!
1425	!! MAIN OUTPUT VARIABLE(S): cgrnd, dgrnd, pcapa, pkappa, soilcap, soilflx
1426	!!
1427	!! REFERENCE(S) :
1428	!! - Hourdin, F. (1992). Study and numerical simulation of the general circulation of planetary atmospheres,
1429	!! Ph.D. thesis, Paris VII University. Remark: the part of F. Hourdin's PhD thesis relative to the thermal
1430	!! integration scheme has been scanned and is provided along with the documentation, with name :
1431	!! Hourdin_1992_PhD_thermal_scheme.pdf
1432	!!
1433	!! FLOWCHART : None
1434	!! \n
1435	!_ ================================================================================================================================
1436
1437	SUBROUTINE thermosoil_coef (kjpindex, temp_sol_new, temp_sol_new_pft, snow, &
1438	soilcap, soilcap_pft, soilflx, soilflx_pft, njsc, &
1439	cgrnd, dgrnd, profil_froz, pcappa_supp, &
1440	organic_layer_thick, soilc_total, veget_max, snowdz, &
1441	snowrho, snowtemp, pb, &
1442	frac_snow_veg, frac_snow_nobio, totfrac_nobio, &
1443	lambda_snow, cgrnd_snow, dgrnd_snow)
1444
1445	!! 0. Variables and parameter declaration
1446
1447	!! 0.1 Input variables
1448
1449	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (unitless)
1450	REAL(r_std), DIMENSION (kjpindex), INTENT (in) :: temp_sol_new !! soil surface temperature @tex ($K$) @endtex
1451	REAL(r_std), DIMENSION (kjpindex,nvm), INTENT (in) :: temp_sol_new_pft !! soil surface temperature @tex ($K$) @endtex
1452	REAL(r_std), DIMENSION (kjpindex), INTENT (in) :: snow !! snow mass @tex ($Kg$) @endtex
1453	INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: njsc !! Index of the dominant soil textural class in the grid cell (1-nscm, unitless)
1454	REAL(r_std), DIMENSION (kjpindex,nvm), INTENT (in) :: veget_max !!Fraction of vegetation type
1455	REAL(r_std), DIMENSION(kjpindex), INTENT (in) :: organic_layer_thick !! how deep is the organic soil?
1456	REAL(r_std), DIMENSION(kjpindex,ndeep,nvm), INTENT (in) :: soilc_total !! total soil carbon for use in thermal calcs
1457	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: frac_snow_veg !! Snow cover fraction on vegeted area
1458	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT(in) :: frac_snow_nobio !! Snow cover fraction on non-vegeted area
1459	REAL(r_std),DIMENSION (kjpindex),INTENT(in) :: totfrac_nobio !! Total fraction of continental ice+lakes+cities+...
1460	!!(unitless,0-1)
1461	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(in) :: snowdz !! Snow depth (m)
1462	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(in) :: snowrho !! Snow density
1463	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(in) :: snowtemp !! Snow temperature (K)
1464	REAL(r_std), DIMENSION (kjpindex), INTENT (in) :: pb !! Surface presure (hPa)
1465
1466	!! 0.2 Output variables
1467
1468	REAL(r_std), DIMENSION (kjpindex,ngrnd-1,nvm), INTENT(out):: cgrnd !! matrix coefficient for the computation of soil
1469	!! temperatures (beta in F. Hourdin thesis)
1470	REAL(r_std), DIMENSION (kjpindex,ngrnd-1,nvm), INTENT(out):: dgrnd !! matrix coefficient for the computation of soil
1471	!! temperatures (alpha in F. Hourdin thesis)
1472	REAL(r_std), DIMENSION (kjpindex,ngrnd,nvm), INTENT(out) :: profil_froz
1473	REAL(r_std),DIMENSION(kjpindex,ngrnd,nvm),INTENT(out) :: pcappa_supp
1474	REAL(r_std), DIMENSION (kjpindex), INTENT (out) :: soilcap !! surface heat capacity considering snow and soil surface
1475	REAL(r_std), DIMENSION (kjpindex,nvm), INTENT (out) :: soilcap_pft !! surface heat capacity
1476	!! @tex ($J m^{-2} K^{-1}$) @endtex
1477	REAL(r_std), DIMENSION (kjpindex), INTENT (out) :: soilflx !! surface heat flux considering snow and soil surface @tex ($W m^{-2}$) @endtex,
1478	REAL(r_std), DIMENSION (kjpindex,nvm), INTENT (out) :: soilflx_pft !! surface heat flux @tex ($W m^{-2}$) @endtex,
1479	!! positive towards the
1480	!! soil, writen as Qg (ground heat flux) in the history
1481	!! files.
1482
1483	!! 0.3 Modified variable
1484
1485	REAL(r_std), DIMENSION (kjpindex), INTENT(inout) :: lambda_snow !! Coefficient of the linear extrapolation of surface temperature
1486	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT (inout):: cgrnd_snow !! Integration coefficient for snow numerical scheme
1487	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT (inout):: dgrnd_snow !! Integration coefficient for snow numerical scheme
1488
1489	!! 0.4 Local variables
1490
1491
1492	INTEGER(i_std) :: ji, jg,jv
1493	REAL(r_std), DIMENSION (kjpindex,ngrnd-1,nvm) :: zdz1
1494
1495	REAL(r_std), DIMENSION (kjpindex,ngrnd,nvm) :: zdz2
1496	REAL(r_std), DIMENSION (kjpindex) :: z1 !! numerical constant @tex ($W m^{-1} K^{-1}$) @endtex
1497
1498
1499	REAL(r_std), DIMENSION (kjpindex) :: soilcap_nosnow!! surface heat capacity
1500	!! @tex ($J m^{-2} K^{-1}$)
1501	!! @endtex
1502	REAL(r_std), DIMENSION (kjpindex) :: soilflx_nosnow!! surface heat flux @tex ($W m^{-2}$) @endtex,
1503	!! positive towards the soil, written as Qg
1504	!!(ground heat flux in the history files).
1505
1506	REAL(r_std), DIMENSION (kjpindex) :: cgrnd_soil !! surface soil layer
1507	REAL(r_std), DIMENSION (kjpindex) :: dgrnd_soil !! surface soil layer
1508	REAL(r_std), DIMENSION (kjpindex) :: zdz1_soil !! surface soil layer
1509	REAL(r_std), DIMENSION (kjpindex) :: zdz2_soil !! surface soil layer
1510
1511	!_ ================================================================================================================================
1512
1513	!! Initialisation of local variables
1514	z1(:) = zero
1515	zdz2(:,:,:) = zero
1516	zdz1(:,:,:) = zero
1517	soilcap_nosnow(:) = zero
1518	soilflx_nosnow(:) = zero
1519	cgrnd_soil(:) = zero
1520	dgrnd_soil(:) = zero
1521	zdz1_soil(:) = zero
1522	zdz2_soil(:) = zero
1523	soilcap(:) = zero
1524	soilflx(:) = zero
1525	soilcap_pft(:,:) = zero
1526	soilflx_pft(:,:) = zero
1527
1528	!! 1. Computation of the soil thermal properties
1529
1530	! Computation of the soil thermal properties; snow properties are also accounted for
1531	IF (ok_explicitsnow) THEN
1532	CALL thermosoil_getdiff( kjpindex, ptn, njsc, veget_max, shum_ngrnd_permalong, &
1533	profil_froz, pcappa_supp, organic_layer_thick, soilc_total, snowrho, &
1534	snowtemp, pb, mc_layt, mc_layt_pft, tmc_layt_pft, pcapa, pcapa_en, pkappa)
1535
1536	! this is for the thin snow in order to prevent the warm surface
1537	! CALL thermosoil_getdiff_thinsnow (kjpindex, ptn, shum_ngrnd_permalong, snowdz,profil_froz)
1538	ELSE
1539	CALL thermosoil_getdiff_old_thermix_with_snow( kjpindex, snow, njsc )
1540	ENDIF
1541
1542	! ok_freeze_thermix must be true
1543	! IF (ok_Ecorr) THEN
1544	! CALL thermosoil_readjust(kjpindex, ptn)
1545	! ENDIF
1546
1547	!! 2. Computation of the coefficients of the numerical integration scheme for the soil layers
1548
1549	!! 2.1 Calculate numerical coefficients zdz1 and zdz2
1550
1551	DO jv=1,nvm
1552	DO jg=1,ngrnd
1553	zdz2(:,jg,jv)=pcapa(:,jg,jv) * dlt(jg)/dt_sechiba
1554	ENDDO ! DO jg=1,ngrnd
1555
1556	DO jg=1,ngrnd-1
1557	zdz1(:,jg,jv) = dz1(jg) * pkappa(:,jg,jv)
1558	ENDDO !DO jg=1,ngrnd-1
1559
1560
1561	!! 2.2 Calculate coefficients cgrnd and dgrnd for soil
1562	z1(:) = zdz2(:,ngrnd,jv) + zdz1(:,ngrnd-1,jv)
1563	cgrnd(:,ngrnd-1,jv) = (phigeoth + zdz2(:,ngrnd,jv) * ptn(:,ngrnd,jv)) / z1(:)
1564	dgrnd(:,ngrnd-1,jv) = zdz1(:,ngrnd-1,jv) / z1(:)
1565	DO jg = ngrnd-1,2,-1
1566	z1(:) = un / (zdz2(:,jg,jv) + zdz1(:,jg-1,jv) + zdz1(:,jg,jv) * (un - dgrnd(:,jg,jv)))
1567	cgrnd(:,jg-1,jv) = (ptn(:,jg,jv) * zdz2(:,jg,jv) + zdz1(:,jg,jv) * cgrnd(:,jg,jv)) * z1(:)
1568	dgrnd(:,jg-1,jv) = zdz1(:,jg-1,jv) * z1(:)
1569	ENDDO ! jg = ngrnd-1,2,-1
1570
1571	!! 3. Computation of the apparent ground heat flux
1572
1573	!! Computation of the apparent ground heat flux (> towards the soil) and
1574	!! apparent surface heat capacity, used at the next timestep by enerbil to
1575	!! compute the surface temperature.
1576
1577	!! no snow involved
1578	IF (ok_LAIdev(jv)) THEN
1579	!!! calculating soil heat flux with PFT specific surface temperature
1580	soilflx_pft(:,jv) = zdz1(:,1,jv) * (cgrnd(:,1,jv) + (dgrnd(:,1,jv)-1.) * ptn(:,1,jv))
1581	soilcap_pft(:,jv) = (zdz2(:,1,jv) * dt_sechiba + dt_sechiba * (un - dgrnd(:,1,jv)) * zdz1(:,1,jv))
1582	z1(:) = lambda * (un - dgrnd(:,1,jv)) + un
1583	soilcap_pft(:,jv) = soilcap_pft(:,jv) / z1(:)
1584	soilflx_pft(:,jv) = soilflx_pft(:,jv) + &
1585	& soilcap_pft(:,jv) * (ptn(:,1,jv) * z1(:) - lambda * cgrnd(:,1,jv) - temp_sol_new_pft(:,jv)) / dt_sechiba
1586	ELSE
1587	soilflx_pft(:,jv) = zdz1(:,1,jv) * (cgrnd(:,1,jv) + (dgrnd(:,1,jv)-1.) * ptn(:,1,jv))
1588	soilcap_pft(:,jv) = (zdz2(:,1,jv) * dt_sechiba + dt_sechiba * (un - dgrnd(:,1,jv)) * zdz1(:,1,jv))
1589	z1(:) = lambda * (un - dgrnd(:,1,jv)) + un
1590	soilcap_pft(:,jv) = soilcap_pft(:,jv) / z1(:)
1591	soilflx_pft(:,jv) = soilflx_pft(:,jv) + &
1592	& soilcap_pft(:,jv) * (ptn(:,1,jv) * z1(:) - lambda * cgrnd(:,1,jv) - temp_sol_new(:)) / dt_sechiba
1593	ENDIF
1594	ENDDO ! jv=1,nvm
1595
1596
1597	! 4 here is where I normalize to take the weighted means of each of the
1598	! PFTs for surface energy fluxes
1599
1600	DO ji = 1,kjpindex
1601	DO jv=1,nvm !pft
1602	soilflx_nosnow(ji) = soilflx_nosnow(ji) + (soilflx_pft(ji,jv)*veget_max(ji,jv))
1603	soilcap_nosnow(ji) = soilcap_nosnow(ji) + (soilcap_pft(ji,jv)*veget_max(ji,jv))
1604	cgrnd_soil(ji) = cgrnd_soil(ji) + (cgrnd(ji,1,jv)*veget_max(ji,jv))
1605	dgrnd_soil(ji) = dgrnd_soil(ji) + (dgrnd(ji,1,jv)*veget_max(ji,jv))
1606	zdz1_soil(ji) = zdz1_soil(ji) + (zdz1(ji,1,jv)*veget_max(ji,jv))
1607	zdz2_soil(ji) = zdz2_soil(ji) + (zdz2(ji,1,jv)*veget_max(ji,jv))
1608	END DO
1609	END DO
1610
1611	!! 3. Computation of the coefficients of the numerical integration scheme for the snow layers
1612	IF (ok_explicitsnow) THEN
1613	CALL thermosoil_coef_snow(kjpindex, temp_sol_new, snow, &
1614	soilcap_nosnow, soilflx_nosnow,soilcap, soilflx, &
1615	cgrnd, dgrnd, cgrnd_soil, dgrnd_soil, &
1616	snowdz, &
1617	snowtemp, zdz1_soil, zdz2_soil, &
1618	frac_snow_veg, frac_snow_nobio, totfrac_nobio, &
1619	lambda_snow, cgrnd_snow, dgrnd_snow)
1620	ELSE
1621	lambda_snow(:) = lambda
1622	cgrnd_snow(:,:) = zero
1623	dgrnd_snow(:,:) = zero
1624	soilcap(:)=soilcap_nosnow(:)
1625	soilflx(:)=soilflx_nosnow(:)
1626	ENDIF
1627
1628	IF (printlev>=3) WRITE (numout,*) ' thermosoil_coef done '
1629
1630	END SUBROUTINE thermosoil_coef
1631
1632
1633
1634	!! ================================================================================================================================
1635	!! SUBROUTINE : thermosoil_coef_snow
1636	!!
1637	!>\BRIEF Calculate soil thermal snow properties
1638	!!
1639	!! DESCRIPTION
1640	!!
1641	!! RECENT CHANGE(S) : None
1642	!!
1643	!! MAIN OUTPUT VARIABLE(S): cgrnd_snow, dgrnd_snow, soilcap, soilflx
1644	!!
1645	!! REFERENCE(S) :
1646	!! - Hourdin, F. (1992). Study and numerical simulation of the general circulation of planetary atmospheres,
1647	!! Ph.D. thesis, Paris VII University. Remark: the part of F. Hourdin's PhD thesis relative to the thermal
1648	!! integration scheme has been scanned and is provided along with the documentation, with name :
1649	!! Hourdin_1992_PhD_thermal_scheme.pdf
1650	!!
1651	!! FLOWCHART : None
1652	!! \n
1653	!_ ================================================================================================================================
1654	SUBROUTINE thermosoil_coef_snow(kjpindex, temp_sol_new, snow, &
1655	soilcap_nosnow, soilflx_nosnow, soilcap, soilflx, &
1656	cgrnd, dgrnd, cgrnd_soil, dgrnd_soil, &
1657	snowdz, &
1658	snowtemp, zdz1_soil, zdz2_soil, &
1659	frac_snow_veg, frac_snow_nobio, totfrac_nobio, &
1660	lambda_snow, cgrnd_snow, dgrnd_snow)
1661
1662	!! 0. Variables and parameter declaration
1663
1664	!! 0.1 Input variables
1665
1666	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (unitless)
1667	REAL(r_std), DIMENSION (kjpindex), INTENT (in) :: temp_sol_new !! soil surface temperature @tex ($K$) @endtex
1668	REAL(r_std), DIMENSION (kjpindex), INTENT (in) :: snow !! snow mass @tex ($Kg$) @endtex
1669	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: frac_snow_veg !! Snow cover fraction on vegeted area
1670	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT(in) :: frac_snow_nobio !! Snow cover fraction on non-vegeted area
1671	REAL(r_std),DIMENSION (kjpindex),INTENT(in) :: totfrac_nobio !! Total fraction of continental ice+lakes+cities+...
1672	!!(unitless,0-1)
1673	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(in) :: snowdz !! Snow depth (m)
1674	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(in) :: snowtemp !! Snow temperature (K)
1675
1676	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: zdz1_soil !! surface soil layer
1677	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: zdz2_soil !! surface soil layer
1678
1679	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: soilcap_nosnow!! surface heat capacity
1680	!! @tex ($J m^{-2} K^{-1}$)
1681	!! @endtex
1682	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: soilflx_nosnow!! surface heat flux @tex ($W m^{-2}$) @endtex,
1683	!! positive towards the soil, written as Qg
1684	!!(ground heat flux in the history files).
1685	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: cgrnd_soil !! surface soil layer
1686	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: dgrnd_soil !! surface soil layer
1687	!! 0.2 Output variables
1688
1689	REAL(r_std), DIMENSION (kjpindex,ngrnd-1,nvm), INTENT(in):: cgrnd !! matrix coefficient for the computation of soil
1690	!! temperatures (beta in F. Hourdin thesis)
1691	REAL(r_std), DIMENSION (kjpindex,ngrnd-1,nvm), INTENT(in):: dgrnd !! matrix coefficient for the computation of soil
1692	!! temperatures (alpha in F. Hourdin thesis)
1693	REAL(r_std), DIMENSION (kjpindex), INTENT (out) :: soilcap !! surface heat capacity considering snow and soil surface
1694	REAL(r_std), DIMENSION (kjpindex), INTENT (out) :: soilflx !! surface heat flux considering snow and soil surface @tex ($W m^{-2}$) @endtex,
1695
1696	!! 0.3 Modified variable
1697
1698	REAL(r_std), DIMENSION (kjpindex), INTENT(out) :: lambda_snow !! Coefficient of the linear extrapolation of surface temperature
1699	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT (out):: cgrnd_snow !! Integration coefficient for snow numerical scheme
1700	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT (out):: dgrnd_snow !! Integration coefficient for snow numerical scheme
1701
1702	!! 0.4 Local variables
1703
1704	REAL(r_std), DIMENSION (kjpindex) :: snowcap !! apparent snow heat capacity @tex ($J m^{-2} K^{-1}$)
1705	REAL(r_std), DIMENSION (kjpindex) :: snowflx !! apparent snow-atmosphere heat flux @tex ($W m^{-2}$) @endtex
1706	REAL(r_std), DIMENSION (kjpindex,nsnow) :: dz1_snow
1707	REAL(r_std), DIMENSION (kjpindex,nsnow) :: ZSNOWDZM
1708	REAL(r_std), DIMENSION (kjpindex,nsnow) :: dz2_snow
1709	REAL(r_std), DIMENSION (kjpindex,nsnow) :: zdz1_snow
1710	REAL(r_std), DIMENSION (kjpindex,nsnow) :: zdz2_snow
1711	REAL(r_std), DIMENSION (kjpindex) :: z1_snow
1712
1713	INTEGER(i_std) :: ji, jg,jv
1714	REAL(r_std), DIMENSION (kjpindex) :: z1 !! numerical constant @tex ($W m^{-1} K^{-1}$) @endtex
1715
1716	!_ ================================================================================================================================
1717
1718	!! Initialisation of local variables
1719	snowcap(:) = zero
1720	snowflx(:) = zero
1721	dz1_snow(:,:) = zero
1722	ZSNOWDZM(:,:) = zero
1723	dz2_snow(:,:) = zero
1724	zdz1_snow(:,:) = zero
1725	zdz2_snow(:,:) = zero
1726	z1_snow(:) = zero
1727
1728	IF (printlev>=3) WRITE (numout,*) ' thermosoil_coef_snow start '
1729
1730	IF (.NOT. ok_explicitsnow) THEN
1731	CALL ipslerr_p(3, 'thermosoi_coef_snow', 'ok_explicitsnow must be enabled to call this subroutine', '', '')
1732	ENDIF
1733
1734	!! 1. Computation of the coefficients of the numerical integration scheme for the snow layers
1735
1736	!! 1.1 Calculate numerical coefficients zdz1_snow, zdz2_snow and lambda_snow
1737	DO ji = 1, kjpindex
1738	! Calculate internal values
1739	DO jg = 1, nsnow
1740	ZSNOWDZM(ji,jg) = MAX(snowdz(ji,jg),psnowdzmin)
1741	ENDDO
1742	dz2_snow(ji,:)=ZSNOWDZM(ji,:)
1743
1744	DO jg = 1, nsnow-1
1745	dz1_snow(ji,jg) = 2.0 / (dz2_snow(ji,jg+1)+dz2_snow(ji,jg))
1746	ENDDO
1747
1748	lambda_snow(ji) = dz2_snow(ji,1)/2.0 * dz1_snow(ji,1)
1749
1750	DO jg=1,nsnow
1751	zdz2_snow(ji,jg)=pcapa_snow(ji,jg) * dz2_snow(ji,jg)/dt_sechiba
1752	ENDDO
1753
1754	DO jg=1,nsnow-1
1755	zdz1_snow(ji,jg) = dz1_snow(ji,jg) * pkappa_snow(ji,jg)
1756	ENDDO
1757
1758	! the bottom snow
1759	zdz1_snow(ji,nsnow) = pkappa_snow(ji,nsnow) / ( zlt(1) + dz2_snow(ji,nsnow)/2 )
1760
1761	ENDDO
1762
1763	!! 1.2 Calculate coefficients cgrnd_snow and dgrnd_snow for snow
1764	DO ji = 1,kjpindex
1765	! bottom level
1766	z1_snow(ji) = zdz2_soil(ji)+(un-dgrnd_soil(ji))*zdz1_soil(ji)+zdz1_snow(ji,nsnow)
1767	cgrnd_snow(ji,nsnow) = (zdz2_soil(ji) * ptn_pftmean(ji,1) + zdz1_soil(ji) * cgrnd_soil(ji) ) / z1_snow(ji)
1768	dgrnd_snow(ji,nsnow) = zdz1_snow(ji,nsnow) / z1_snow(ji)
1769
1770	! next-to-bottom level
1771	z1_snow(ji) = zdz2_snow(ji,nsnow)+(un-dgrnd_snow(ji,nsnow))*zdz1_snow(ji,nsnow)+zdz1_snow(ji,nsnow-1)
1772	cgrnd_snow(ji,nsnow-1) = (zdz2_snow(ji,nsnow)*snowtemp(ji,nsnow)+&
1773	zdz1_snow(ji,nsnow)*cgrnd_snow(ji,nsnow))/z1_snow(ji)
1774	dgrnd_snow(ji,nsnow-1) = zdz1_snow(ji,nsnow-1) / z1_snow(ji)
1775
1776	DO jg = nsnow-1,2,-1
1777	z1_snow(ji) = un / (zdz2_snow(ji,jg) + zdz1_snow(ji,jg-1) + zdz1_snow(ji,jg) * (un - dgrnd_snow(ji,jg)))
1778	cgrnd_snow(ji,jg-1) = (snowtemp(ji,jg) * zdz2_snow(ji,jg) + zdz1_snow(ji,jg) * cgrnd_snow(ji,jg)) * z1_snow(ji)
1779	dgrnd_snow(ji,jg-1) = zdz1_snow(ji,jg-1) * z1_snow(ji)
1780	ENDDO
1781	ENDDO
1782
1783	!! 2. Computation of the apparent ground heat flux
1784	!! Computation of apparent snow-atmosphere flux
1785	DO ji = 1,kjpindex
1786	snowflx(ji) = zdz1_snow(ji,1) * (cgrnd_snow(ji,1) + (dgrnd_snow(ji,1)-1.) * snowtemp(ji,1))
1787	snowcap(ji) = (zdz2_snow(ji,1) * dt_sechiba + dt_sechiba * (un - dgrnd_snow(ji,1)) * zdz1_snow(ji,1))
1788	z1_snow(ji) = lambda_snow(ji) * (un - dgrnd_snow(ji,1)) + un
1789	snowcap(ji) = snowcap(ji) / z1_snow(ji)
1790	snowflx(ji) = snowflx(ji) + &
1791	& snowcap(ji) * (snowtemp(ji,1) * z1_snow(ji) - lambda_snow(ji) * cgrnd_snow(ji,1) - temp_sol_new(ji)) / dt_sechiba
1792	ENDDO
1793
1794	!! Add snow fraction
1795	! Using an effective heat capacity and heat flux by a simple pondering of snow and soil fraction
1796	DO ji = 1, kjpindex
1797	soilcap(ji) = snowcap(ji)frac_snow_veg(ji)(1-totfrac_nobio(ji))+ & ! weights related to snow cover fraction on vegetation
1798	soilcap_nosnow(ji)SUM(frac_snow_nobio(ji,:))totfrac_nobio(ji)+ & ! weights related to SCF on nobio
1799	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
1800	soilflx(ji) = snowflx(ji)frac_snow_veg(ji)(1-totfrac_nobio(ji))+ & ! weights related to snow cover fraction on vegetation
1801	soilflx_nosnow(ji)SUM(frac_snow_nobio(ji,:))totfrac_nobio(ji)+ & ! weights related to SCF on nobio
1802	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
1803	ENDDO
1804
1805	IF (printlev>=3) WRITE (numout,*) ' thermosoil_coef_snow done '
1806	END SUBROUTINE thermosoil_coef_snow
1807
1808
1809	!! ================================================================================================================================
1810	!! SUBROUTINE : thermosoil_profile
1811	!!
1812	!>\BRIEF In this routine solves the numerical soil thermal scheme, ie calculates the new soil temperature profile;
1813	!!
1814	!!
1815	!! DESCRIPTION : The calculation of the new soil temperature profile is based on
1816	!! the cgrnd and dgrnd values from the previous timestep and the surface temperature Ts aka temp_sol_new. (see detailed
1817	!! explanation in the header of the thermosoil module or in the reference).\n
1818	!! T(k+1)=cgrnd(k)+dgrnd(k)*T(k)\n
1819	!! -- EQ1 --\n
1820	!! Ts=(1+lambda)T(1) -lambdaT(2)\n
1821	!! -- EQ2--\n
1822	!!
1823	!! RECENT CHANGE(S) : None
1824	!!
1825	!! MAIN OUTPUT VARIABLE(S): ptn (soil temperature profile on the thermal axis),
1826	!! stempdiag (soil temperature profile on the diagnostic axis)
1827	!!
1828	!! REFERENCE(S) :
1829	!! - Hourdin, F. (1992). Study and numerical simulation of the general circulation of planetary atmospheres,
1830	!! Ph.D. thesis, Paris VII University. Remark: the part of F. Hourdin's PhD thesis relative to the thermal
1831	!! integration scheme has been scanned and is provided along with the documentation, with name :
1832	!! Hourdin_1992_PhD_thermal_scheme.pdf
1833	!!
1834	!! FLOWCHART : None
1835	!! \n
1836	!_ ================================================================================================================================
1837
1838	SUBROUTINE thermosoil_profile (kjpindex, temp_sol_new, temp_sol_new_pft, ptn, stempdiag, &
1839	snowtemp, veget_max, &
1840	frac_snow_veg, frac_snow_nobio,totfrac_nobio, &
1841	cgrnd_snow, dgrnd_snow)
1842
1843	!! 0. Variables and parameter declaration
1844
1845	!! 0.1 Input variables
1846
1847	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (unitless)
1848	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_sol_new !! Surface temperature at the present time-step
1849	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: temp_sol_new_pft !! Surface temperature at the present time-step
1850	!! @tex ($K$) @endtex
1851	REAL(r_std), DIMENSION (kjpindex,nvm), INTENT (in) :: veget_max !! Fraction of vegetation type
1852	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: frac_snow_veg !! Snow cover fraction on vegeted area
1853	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT(in) :: frac_snow_nobio!! Snow cover fraction on non-vegeted area
1854	REAL(r_std),DIMENSION (kjpindex),INTENT(in) :: totfrac_nobio !! Total fraction of continental ice+lakes+cities+...
1855	REAL(r_std),DIMENSION (kjpindex,nsnow), INTENT(in) :: snowtemp !! Snow temperature (K)
1856	REAL(r_std),DIMENSION (kjpindex,nsnow), INTENT(in) :: cgrnd_snow !! Integration coefficient for snow numerical scheme
1857	REAL(r_std),DIMENSION (kjpindex,nsnow), INTENT(in) :: dgrnd_snow !! Integration coefficient for snow numerical scheme
1858
1859	!! 0.3 Modified variables
1860
1861	!! 0.2 Output variables
1862
1863	REAL(r_std),DIMENSION (kjpindex,nslm), INTENT (out) :: stempdiag !! diagnostic temperature profile
1864	!! @tex ($K$) @endtex
1865
1866	!! 0.3 Modified variables
1867
1868	REAL(r_std),DIMENSION (kjpindex,ngrnd,nvm), INTENT (inout) :: ptn !! vertically discretized soil temperatures
1869	!! @tex ($K$) @endtex
1870
1871	!! 0.4 Local variables
1872
1873	INTEGER(i_std) :: ji, jg, jv
1874	REAL(r_std) :: temp_sol_eff !! effective surface temperature including snow and soil
1875
1876	!_ ================================================================================================================================
1877
1878	!! 1. Computes the snow temperatures
1879
1880	!! 1.1. ptn(jg=1) using EQ1 and EQ2
1881	DO jv = 1,nvm
1882	DO ji = 1,kjpindex
1883	IF (ok_explicitsnow ) THEN !!! again, explicit snow not compatible with crop column temperature, xuhui
1884	! Using an effective surface temperature by a simple pondering
1885	temp_sol_eff=snowtemp(ji,nsnow)frac_snow_veg(ji)(1-totfrac_nobio(ji))+ & ! weights related to snow cover fraction on vegetation
1886	temp_sol_new(ji)SUM(frac_snow_nobio(ji,:))totfrac_nobio(ji)+ & ! weights related to SCF on nobio
1887	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
1888	! Soil temperature calculation with explicit snow if there is snow on the ground
1889	ptn(ji,1,jv) = cgrnd_snow(ji,nsnow) + dgrnd_snow(ji,nsnow) * temp_sol_eff
1890	ELSE
1891	! Standard soil temperature calculation
1892	IF (ok_LAIdev(jv) .AND. veget_max(ji,jv)>0) THEN
1893	ptn(ji,1,jv) = (lambda * cgrnd(ji,1,jv) + temp_sol_new_pft(ji,jv)) / (lambda *(un - dgrnd(ji,1,jv)) + un)
1894	ELSE
1895	ptn(ji,1,jv) = (lambda * cgrnd(ji,1,jv) + temp_sol_new(ji)) / (lambda *(un - dgrnd(ji,1,jv)) + un)
1896	ENDIF
1897	ENDIF
1898	ENDDO
1899
1900	!! 1.2. ptn(jg=2:ngrnd) using EQ1.
1901	DO jg = 1,ngrnd-1
1902	DO ji = 1,kjpindex
1903	ptn(ji,jg+1,jv) = cgrnd(ji,jg,jv) + dgrnd(ji,jg,jv) * ptn(ji,jg,jv)
1904	ENDDO
1905	ENDDO
1906	ENDDO
1907
1908	!! 2. Assigne the soil temperature to the output variable. It is already on the right axis.
1909	CALL thermosoil_diaglev(kjpindex, stempdiag, veget_max)
1910
1911	IF (printlev>=3) WRITE (numout,*) ' thermosoil_profile done '
1912
1913	END SUBROUTINE thermosoil_profile
1914
1915
1916	!================================================================================================================================
1917	!! SUBROUTINE : thermosoil_cond
1918	!!
1919	!>\BRIEF Calculate soil thermal conductivity from Orchidee trunk
1920	!!
1921	!! DESCRIPTION : This routine computes soil thermal conductivity
1922	!! Code introduced from NOAH LSM. Used in Orchidee trunk.
1923	!!
1924	!! RECENT CHANGE(S) : None
1925	!!
1926	!! MAIN OUTPUT VARIABLE(S): cnd
1927	!!
1928	!! REFERENCE(S) :
1929	!! Farouki, O.T.,1986: Thermal Properties of Soils. Series on Rock
1930	!! and Soil Mechanics, Vol. 11, Trans Tech, 136 PP.
1931	!! Johansen, O., 1975: Thermal Conductivity of Soils. Ph.D. Thesis,
1932	!! University of Trondheim,
1933	!! Peters-Lidard, C. D., Blackburn, E., Liang, X., & Wood, E. F.,
1934	!! 1998: The effect of soil thermal conductivity
1935	!! Parameterization on Surface Energy fluxes
1936	!! and Temperatures. J. of The Atmospheric Sciences,
1937	!! Vol. 55, pp. 1209-1224.
1938	!! Modify histroy:
1939	!!
1940	!! FLOWCHART : None
1941	!! \n
1942	!_
1943	!================================================================================================================================
1944
1945	SUBROUTINE thermosoil_cond (kjpindex, njsc, smc, qz, smcmax, sh2o, cnd)
1946
1947	!! 0. Variables and parameter declaration
1948
1949	!! 0.1 Input variables
1950	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (unitless)
1951	INTEGER(i_std), DIMENSION (kjpindex), INTENT (in) :: njsc !! Index of the dominant soil textural class in the grid cell (1-nscm, unitless)
1952	REAL(r_std), DIMENSION (kjpindex,ngrnd), INTENT(IN) :: smc !! Volumetric Soil Moisture Content (m3/m3)
1953	REAL(r_std), DIMENSION (nscm), INTENT(IN) :: qz !! Quartz Content (Soil Type Dependent) (0-1)
1954	REAL(r_std), DIMENSION (nscm), INTENT(IN) :: smcmax !! Soil Porosity (0-1)
1955	REAL(r_std), DIMENSION (kjpindex,ngrnd), INTENT(IN) :: sh2o !! Unfrozen Soil Moisture Content; Frozen Soil Moisture = smc - sh2o
1956
1957	!! 0.2 Output variables
1958	REAL(r_std), DIMENSION (kjpindex,ngrnd), INTENT(OUT) :: cnd !! Soil Thermal Conductivity (W/m/k)
1959
1960	!! 0.3 Modified variables
1961
1962	!! 0.4 Local variables
1963	REAL(r_std), DIMENSION (kjpindex,ngrnd) :: ake !! Kersten Number (unitless)
1964	REAL(r_std), DIMENSION (kjpindex,ngrnd) :: thksat !! Saturated Thermal Conductivity (W/m/k)
1965	REAL(r_std), DIMENSION (kjpindex,ngrnd) :: satratio !! Degree of Saturation (0-1)
1966	REAL(r_std), DIMENSION (kjpindex,ngrnd) :: xu !! Unfrozen Volume For Saturation (0-1)
1967	REAL(r_std), DIMENSION (kjpindex,ngrnd) :: xunfroz !! Unfrozon Volume Fraction (0-1)
1968	REAL(r_std) :: thko !! Thermal Conductivity for Other Ssoil Components (W/m/k)
1969	REAL(r_std) :: gammd !! Dry Dendity (kg/m3)
1970	REAL(r_std) :: thkdry !! Dry Thermal Conductivity (W/m/k)
1971	REAL(r_std) :: thks !! Thermal Conductivity for the Solids Combined (Quartz + Other) (W/m/k)
1972	INTEGER(i_std) :: ji, jg, jst
1973
1974	!_================================================================================================================================
1975
1976	!! 1. Dry and Saturated Thermal Conductivity.
1977
1978	DO ji = 1,kjpindex
1979	jst = njsc(ji)
1980
1981	!! 1.1. Dry density (Kg/m3) and Dry thermal conductivity (W.M-1.K-1)
1982	gammd = (1. - smcmax(jst))*2700.
1983	thkdry = (0.135* gammd+ 64.7)/ (2700. - 0.947* gammd)
1984
1985	!! 1.2. thermal conductivity of "other" soil components
1986	IF (qz(jst) > 0.2) THEN
1987	thko = 2.0
1988	ELSEIF (qz(jst) <= 0.2) THEN
1989	thko = 3.0
1990	ENDIF
1991
1992	!! 1.3. Thermal conductivity of solids
1993	thks = (THKQTZ ** qz(jst))* (thko ** (1. - qz(jst)))
1994
1995	DO jg = 1,ngrnd
1996	!! 1.4. saturation ratio
1997	satratio(ji,jg) = smc(ji,jg) / smcmax(jst)
1998
1999	!! 1.5. Saturated Thermal Conductivity (thksat)
2000	IF ( smc(ji,jg) > min_sechiba ) THEN
2001	xunfroz(ji,jg) = sh2o(ji,jg) / smc(ji,jg) ! Unfrozen Fraction (From i.e., 100%Liquid, to 0. (100% Frozen))
2002	xu(ji,jg) = xunfroz(ji,jg) * smcmax(jst) ! Unfrozen volume for saturation (porosity*xunfroz)
2003	thksat(ji,jg) = thks ** (1. - smcmax(jst))* THKICE ** (smcmax(jst) - xu(ji,jg))* THKW ** (xu(ji,jg))
2004	ELSE
2005	! this value will not be used since ake=0 for this case
2006	thksat(ji,jg)=0
2007	END IF
2008	END DO ! DO jg = 1,ngrnd
2009
2010	!! 2. Kersten Number (ake)
2011	DO jg = 1,ngrnd
2012	IF ( (sh2o(ji,jg) + 0.0005) < smc(ji,jg) ) THEN
2013	! Frozen
2014	ake(ji,jg) = satratio(ji,jg)
2015	ELSE
2016	! Unfrozen
2017	! Eq 11 in Peters-Lidard et al., 1998
2018	IF ( satratio(ji,jg) > 0.1 ) THEN
2019	IF ((jst < 4 .AND. soil_classif == 'usda') .OR. (jst == 1 .AND. soil_classif == 'zobler') ) THEN
2020	! Coarse
2021	ake(ji,jg) = 0.7 * LOG10 (SATRATIO(ji,jg)) + 1.0
2022	ELSE
2023	! Fine
2024	ake(ji,jg) = LOG10 (satratio(ji,jg)) + 1.0
2025	ENDIF
2026	ELSEIF ( satratio(ji,jg) > 0.05 .AND. satratio(ji,jg) <= 0.1 ) THEN
2027	IF ((jst < 4 .AND. soil_classif == 'usda') .OR. (jst == 1 .AND. soil_classif == 'zobler') ) THEN
2028	! Coarse
2029	ake(ji,jg) = 0.7 * LOG10 (satratio(ji,jg)) + 1.0
2030	ELSE
2031	! Fine
2032	ake(ji,jg) = 0.0
2033	ENDIF
2034	ELSE
2035	ake(ji,jg) = 0.0 ! use k = kdry
2036	END IF
2037	END IF
2038	END DO ! DO jg = 1,ngrnd
2039
2040	!! 3. Thermal conductivity (cnd)
2041	DO jg = 1,ngrnd
2042	cnd(ji,jg) = ake(ji,jg) * (thksat(ji,jg) - thkdry) + thkdry
2043	END DO ! DO jg = 1,ngrnd
2044
2045	END DO !DO ji = 1,kjpindex
2046
2047	END SUBROUTINE thermosoil_cond
2048
2049	!================================================================================================================================
2050	!! SUBROUTINE : thermosoil_cond_pft
2051	!!
2052	!>\BRIEF Calculate soil thermal conductivity.
2053	!!
2054	!! DESCRIPTION : This routine computes soil thermal conductivity
2055	!! but considers the fact that soil organic carbon can decrease
2056	!! conductivity
2057	!!
2058	!! thermosoil_cond_pft: original implementation for MICT(pft based). Useful for CROP(3d)
2059	!!
2060	!! thermosoil_cond_nopft: clean version of thermosoil_cond_pft. Useful when
2061	!! ther is not need for PFTs. Results are exactly the same
2062	!! as thermosoil_cond_pft when prescribed SOC (currently 2d) is used.
2063	!! There is a substantial performance improvement.
2064	!!
2065	!! RECENT CHANGE(S) : None
2066	!!
2067	!! MAIN OUTPUT VARIABLE(S): cnd
2068	!!
2069	!! REFERENCE(S) :
2070	!! Farouki, O.T.,1986: Thermal Properties of Soils. Series on Rock
2071	!! and Soil Mechanics, Vol. 11, Trans Tech, 136 PP.
2072	!! Johansen, O., 1975: Thermal Conductivity of Soils. Ph.D. Thesis,
2073	!! University of Trondheim,
2074	!! Peters-Lidard, C. D., Blackburn, E., Liang, X., & Wood, E. F.,
2075	!! 1998: The effect of soil thermal conductivity
2076	!! Parameterization on Surface Energy fluxes
2077	!! and Temperatures. J. of The Atmospheric Sciences,
2078	!! Vol. 55, pp. 1209-1224.
2079	!! Lawrence and Slater,2008: Incorporating organic soil into a global climate
2080	!! model
2081	!! Modify histroy:
2082	!!
2083	!! FLOWCHART : None
2084	!! \n
2085	!_
2086	!================================================================================================================================
2087
2088	SUBROUTINE thermosoil_cond_pft (kjpindex, njsc, smc, qz, smcmax, sh2o,zx1,zx2,porosnet,cnd)
2089
2090	!! 0. Variables and parameter declaration
2091
2092	!! 0.1 Input variables
2093	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (unitless)
2094	INTEGER(i_std), DIMENSION (kjpindex), INTENT (in) :: njsc !! Index of the dominant soil textural class in the grid cell (1-nscm, unitless)
2095	REAL(r_std), DIMENSION (kjpindex,ngrnd), INTENT(IN) :: smc !! Volumetric Soil Moisture Content (m3/m3)
2096	!!! xuhui: smc should be a 3-D variable in order to have PFT specific column of
2097	!soil moisture
2098	! REAL(r_std), DIMENSION (kjpindex,ngrnd,nvm), INTENT(IN) :: smc !! Volumetric Soil Moisture Content (m3/m3)
2099	REAL(r_std), DIMENSION (nscm), INTENT(IN) :: qz !! Quartz Content (Soil Type Dependent) (0-1)
2100	REAL(r_std), DIMENSION (nscm), INTENT(IN) :: smcmax !! Soil Porosity (0-1)
2101	REAL(r_std), DIMENSION (kjpindex,ngrnd,nvm), INTENT(IN) :: porosnet !! Soil Porosity (0-1)
2102	REAL(r_std), DIMENSION (kjpindex,ngrnd), INTENT(IN) :: sh2o !! Unfrozen Soil Moisture Content; Frozen Soil Moisture = smc - sh2o
2103	! xuhui: sh2o should be a 3-d variable in order to have PFT-specific column of
2104	! thermal conductivity
2105	! REAL(r_std), DIMENSION (kjpindex,ngrnd,nvm), INTENT(IN) :: sh2o !! Unfrozen Soil Moisture Content; Frozen Soil Moisture = smc - sh2o
2106	REAL(r_std), DIMENSION(kjpindex,ngrnd,nvm), INTENT(IN) :: zx1, zx2 !! proportion of organic and mineral soil
2107
2108	!! 0.2 Output variables
2109	REAL(r_std), DIMENSION (kjpindex,ngrnd,nvm), INTENT(OUT) :: cnd !! Soil Thermal Conductivity (W/m/k)
2110
2111	!! 0.3 Modified variables
2112
2113	!! 0.4 Local variables
2114	REAL(r_std), DIMENSION (kjpindex,ngrnd,nvm) :: ake !! Kerston Number (unitless)
2115	REAL(r_std), DIMENSION (kjpindex,ngrnd,nvm) :: thksat !! Saturated Thermal Conductivity (W/m/k)
2116	REAL(r_std), DIMENSION (kjpindex,ngrnd,nvm) :: satratio !! Degree of Saturation (0-1)
2117	REAL(r_std), DIMENSION (kjpindex,ngrnd,nvm) :: xu !! Unfrozen Volume For Saturation (0-1)
2118	REAL(r_std), DIMENSION (kjpindex,ngrnd,nvm) :: xunfroz !! Unfrozon Volume Fraction (0-1)
2119	REAL(r_std) :: thko !! Thermal Conductivity for Other Ssoil Components (W/m/k)
2120	REAL(r_std), DIMENSION (kjpindex) :: gammd !! Dry Density (kg/m3)
2121	REAL(r_std), DIMENSION (kjpindex) :: thkdry_min !! Dry Thermal Conductivity for mineral soil (W/m/k)
2122	REAL(r_std), DIMENSION (kjpindex,ngrnd,nvm) :: thkdry !! Dry Thermal Conductivity considering organic carbon (W/m/k)
2123	REAL(r_std), DIMENSION (kjpindex) :: thks_min !! Thermal Conductivity for the Solids Combined (Quartz + Other) (W/m/k)
2124	REAL(r_std), DIMENSION (kjpindex,ngrnd,nvm) :: thks !! Thermal Conductivity considering organic carbon (W/m/k)
2125	INTEGER(i_std) :: ji, jg, jst, jv
2126
2127	!_================================================================================================================================
2128
2129	thksat(:,:,:)=0
2130	!! 1. Dry and Saturated Thermal Conductivity.
2131
2132	DO ji = 1,kjpindex
2133	jst = njsc(ji)
2134
2135	!! 1.1. Dry density (Kg/m3) and Dry thermal conductivity (W.M-1.K-1)
2136	gammd(ji) = (1. - smcmax(jst))*2700.
2137	thkdry_min(ji) = (0.135* gammd(ji) + 64.7)/ (2700. - 0.947* gammd(ji))
2138
2139
2140	!! 1.2. thermal conductivity of "other" soil components
2141	IF (qz(jst) > 0.2) THEN
2142	thko = 2.0
2143	ELSEIF (qz(jst) <= 0.2) THEN
2144	thko = 3.0
2145	ENDIF
2146
2147	!! 1.3. Thermal conductivity of solids
2148	thks_min(ji) = (THKQTZ ** qz(jst))* (thko ** (1. - qz(jst)))
2149	ENDDO
2150
2151	DO jv = 1,nvm
2152
2153	SELECTCASE (use_soilc_method)
2154	CASE (SOILC_METHOD_ARITHMETIC)
2155	DO jg = 1,ngrnd
2156	DO ji = 1,kjpindex
2157	thks(ji,jg,jv) = zx1(ji,jg,jv) * cond_solid_org + zx2(ji,jg,jv) * thks_min(ji)
2158	ENDDO
2159	ENDDO
2160	CASE (SOILC_METHOD_GEOMETRIC)
2161	DO jg = 1,ngrnd
2162	DO ji = 1,kjpindex
2163	! use geometric mean rather than arithmetic mean (Decharme et al 2016)
2164	thks(ji,jg,jv) =(cond_solid_org*zx1(ji,jg,jv)) (thks_min(ji)**zx2(ji,jg,jv))
2165	ENDDO
2166	ENDDO
2167	ENDSELECT
2168
2169	DO jg = 1,ngrnd
2170	!! 1.4. saturation ratio
2171	DO ji = 1,kjpindex
2172	satratio(ji,jg,jv) = smc(ji,jg) / porosnet(ji, jg, jv)
2173	ENDDO
2174
2175	!! 1.5. Saturated Thermal Conductivity (thksat)
2176	DO ji = 1,kjpindex
2177	IF ( smc(ji,jg) > min_sechiba ) THEN
2178	xunfroz(ji,jg,jv) = sh2o(ji,jg) / smc(ji,jg) ! Unfrozen Fraction (From i.e., 100%Liquid, to 0. (100% Frozen))
2179	xu(ji,jg,jv) = xunfroz(ji,jg,jv) * porosnet(ji, jg, jv) ! Unfrozen volume for saturation (porosity*xunfroz)
2180	thksat(ji,jg,jv) = thks(ji,jg,jv) ** (1. - porosnet(ji, jg, jv))* THKICE ** (porosnet(ji, jg, jv) - xu(ji,jg,jv)) * THKW ** xu(ji,jg,jv)
2181	END IF
2182	ENDDO
2183	END DO ! DO jg = 1,ngrnd
2184
2185	!! 2. Kerston Number (ake)
2186	DO jg = 1,ngrnd
2187	DO ji = 1,kjpindex
2188	jst = njsc(ji)
2189	IF ( (sh2o(ji,jg) + 0.0005) < smc(ji,jg) ) THEN
2190	! Frozen
2191	ake(ji,jg,jv) = satratio(ji,jg,jv)
2192	ELSE
2193	! Unfrozen
2194	IF ( satratio(ji,jg,jv) > 0.1 ) THEN
2195	IF ((jst < 4 .AND. soil_classif == 'usda') .OR. (jst == 1 .AND. soil_classif == 'zobler') ) THEN
2196	! Coarse
2197	ake(ji,jg,jv) = 0.7 * LOG10 (satratio(ji,jg,jv)) + 1.0
2198	ELSE
2199	! Fine
2200	ake(ji,jg,jv) = LOG10 (satratio(ji,jg,jv)) + 1.0
2201	ENDIF
2202	ELSEIF ( satratio(ji,jg,jv) > 0.05 .AND. satratio(ji,jg,jv) <= 0.1 ) THEN
2203	IF ((jst < 4 .AND. soil_classif == 'usda') .OR. (jst == 1 .AND. soil_classif == 'zobler') ) THEN
2204	! Coarse
2205	ake(ji,jg,jv) = 0.7 * LOG10 (satratio(ji,jg,jv)) + 1.0
2206	ELSE
2207	! Fine
2208	ake(ji,jg,jv) = 0.0
2209	ENDIF
2210	ELSE
2211	ake(ji,jg,jv) = 0.0 ! use k = kdry
2212	END IF
2213	END IF
2214	ENDDO
2215	END DO ! DO jg = 1,ngrnd
2216
2217	SELECTCASE (use_soilc_method)
2218	CASE (SOILC_METHOD_ARITHMETIC)
2219	DO jg = 1,ngrnd
2220	DO ji = 1,kjpindex
2221	thkdry(ji,jg,jv) = zx1(ji,jg,jv) * cond_dry_org + zx2(ji,jg,jv) * thkdry_min(ji)
2222	ENDDO
2223	ENDDO
2224	CASE (SOILC_METHOD_GEOMETRIC)
2225	DO jg = 1,ngrnd
2226	DO ji = 1,kjpindex
2227	! use geometric mean rather than arithmetic mean (Decharme et al 2016)
2228	thkdry(ji,jg,jv) =(cond_dry_org*zx1(ji,jg,jv)) (thkdry_min(ji)**zx2(ji,jg,jv))
2229	ENDDO
2230	ENDDO
2231	CASE DEFAULT
2232	CALL ipslerr_p(3,'thermosoil_cond_pft','Unsupported USE_SOILC_METHOD','','')
2233	ENDSELECT
2234
2235	!! 3. Thermal conductivity (cnd)
2236	DO jg = 1,ngrnd
2237	DO ji = 1,kjpindex
2238	cnd(ji,jg,jv) = ake(ji,jg,jv) * (thksat(ji,jg,jv) - thkdry(ji, jg, jv)) + thkdry(ji, jg, jv)
2239	ENDDO
2240	END DO
2241
2242	END DO !DO jv = 1,nvm
2243
2244	END SUBROUTINE thermosoil_cond_pft
2245
2246
2247	SUBROUTINE thermosoil_cond_nopft (kjpindex, njsc, smc, qz, smcmax, sh2o,zx1,zx2,porosnet,cnd)
2248
2249	!! 0. Variables and parameter declaration
2250
2251	!! 0.1 Input variables
2252	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (unitless)
2253	INTEGER(i_std), DIMENSION (:), INTENT (in) :: njsc !! Index of the dominant soil textural class in the grid cell (1-nscm, unitless)
2254	REAL(r_std), DIMENSION (:,:), INTENT(IN) :: smc !! Volumetric Soil Moisture Content (m3/m3)
2255	REAL(r_std), DIMENSION (:), INTENT(IN) :: qz !! Quartz Content (Soil Type Dependent) (0-1)
2256	REAL(r_std), DIMENSION (:), INTENT(IN) :: smcmax !! Soil Porosity (0-1)
2257	REAL(r_std), DIMENSION (:,:), INTENT(IN) :: porosnet !! Soil Porosity (0-1)
2258	REAL(r_std), DIMENSION (:,:), INTENT(IN) :: sh2o !! Unfrozen Soil Moisture Content; Frozen Soil Moisture = smc - sh2o
2259	REAL(r_std), DIMENSION(:,:), INTENT(IN) :: zx1, zx2 !! proportion of organic and mineral soil
2260
2261	!! 0.2 Output variables
2262	REAL(r_std), DIMENSION (:,:), INTENT(OUT) :: cnd !! Soil Thermal Conductivity (W/m/k)
2263
2264	!! 0.3 Modified variables
2265
2266	!! 0.4 Local variables
2267	REAL(r_std), DIMENSION (kjpindex,ngrnd) :: ake !! Kerston Number (unitless)
2268	REAL(r_std), DIMENSION (kjpindex,ngrnd) :: thksat !! Saturated Thermal Conductivity (W/m/k)
2269	REAL(r_std), DIMENSION (kjpindex,ngrnd) :: satratio !! Degree of Saturation (0-1)
2270	REAL(r_std), DIMENSION (kjpindex,ngrnd) :: xu !! Unfrozen Volume For Saturation (0-1)
2271	REAL(r_std), DIMENSION (kjpindex,ngrnd) :: xunfroz !! Unfrozon Volume Fraction (0-1)
2272	REAL(r_std) :: thko !! Thermal Conductivity for Other Ssoil Components (W/m/k)
2273	REAL(r_std), DIMENSION (kjpindex) :: gammd !! Dry Density (kg/m3)
2274	REAL(r_std), DIMENSION (kjpindex) :: thkdry_min !! Dry Thermal Conductivity for mineral soil (W/m/k)
2275	REAL(r_std), DIMENSION (kjpindex,ngrnd) :: thkdry !! Dry Thermal Conductivity considering organic carbon (W/m/k)
2276	REAL(r_std), DIMENSION (kjpindex) :: thks_min !! Thermal Conductivity for the Solids Combined (Quartz + Other) (W/m/k)
2277	REAL(r_std), DIMENSION (kjpindex,ngrnd) :: thks !! Thermal Conductivity considering organic carbon (W/m/k)
2278	INTEGER(i_std) :: ji, jg, jst
2279
2280	!_================================================================================================================================
2281
2282	thksat(:,:)=0
2283	!! 1. Dry and Saturated Thermal Conductivity.
2284
2285	DO ji = 1,kjpindex
2286	jst = njsc(ji)
2287
2288	!! 1.1. Dry density (Kg/m3) and Dry thermal conductivity (W.M-1.K-1)
2289	gammd(ji) = (1. - smcmax(jst))*2700.
2290	thkdry_min(ji) = (0.135* gammd(ji) + 64.7)/ (2700. - 0.947* gammd(ji))
2291
2292
2293	!! 1.2. thermal conductivity of "other" soil components
2294	IF (qz(jst) > 0.2) THEN
2295	thko = 2.0
2296	ELSEIF (qz(jst) <= 0.2) THEN
2297	thko = 3.0
2298	ENDIF
2299
2300	!! 1.3. Thermal conductivity of solids
2301	thks_min(ji) = (THKQTZ ** qz(jst))* (thko ** (1. - qz(jst)))
2302	ENDDO
2303
2304	SELECTCASE (use_soilc_method)
2305	CASE (SOILC_METHOD_ARITHMETIC)
2306	DO jg = 1,ngrnd
2307	DO ji = 1,kjpindex
2308	thks(ji,jg) = zx1(ji,jg) * cond_solid_org + zx2(ji,jg) * thks_min(ji)
2309	ENDDO
2310	ENDDO
2311	CASE (SOILC_METHOD_GEOMETRIC)
2312	DO jg = 1,ngrnd
2313	DO ji = 1,kjpindex
2314	! use geometric mean rather than arithmetic mean (Decharme et al 2016)
2315	thks(ji,jg) =(cond_solid_org*zx1(ji,jg)) (thks_min(ji)**zx2(ji,jg))
2316	ENDDO
2317	ENDDO
2318	ENDSELECT
2319
2320	DO jg = 1,ngrnd
2321	!! 1.4. saturation ratio
2322	DO ji = 1,kjpindex
2323	satratio(ji,jg) = smc(ji,jg) / porosnet(ji, jg)
2324	ENDDO
2325
2326	!! 1.5. Saturated Thermal Conductivity (thksat)
2327	DO ji = 1,kjpindex
2328	IF ( smc(ji,jg) > min_sechiba ) THEN
2329	xunfroz(ji,jg) = sh2o(ji,jg) / smc(ji,jg) ! Unfrozen Fraction (From i.e., 100%Liquid, to 0. (100% Frozen))
2330	xu(ji,jg) = xunfroz(ji,jg) * porosnet(ji, jg) ! Unfrozen volume for saturation (porosity*xunfroz)
2331	thksat(ji,jg) = thks(ji,jg) ** (1. - porosnet(ji, jg))* THKICE ** (porosnet(ji, jg) - xu(ji,jg)) * THKW ** xu(ji,jg)
2332	END IF
2333	ENDDO
2334	END DO ! DO jg = 1,ngrnd
2335
2336	!! 2. Kerston Number (ake)
2337	DO jg = 1,ngrnd
2338	DO ji = 1,kjpindex
2339	IF ( (sh2o(ji,jg) + 0.0005) < smc(ji,jg) ) THEN
2340	! Frozen
2341	ake(ji,jg) = satratio(ji,jg)
2342	ELSE
2343	! Unfrozen
2344	IF ( satratio(ji,jg) > 0.1 ) THEN
2345	ake(ji,jg) = LOG10 (satratio(ji,jg)) + 1.0
2346	ELSEIF ( satratio(ji,jg) > 0.05 .AND. satratio(ji,jg) <= 0.1 ) THEN
2347	ake(ji,jg) = 0.7 * LOG10 (satratio(ji,jg)) + 1.0
2348	ELSE
2349	ake(ji,jg) = 0.0 ! use k = kdry
2350	END IF
2351	END IF
2352	ENDDO
2353	END DO ! DO jg = 1,ngrnd
2354
2355	SELECTCASE (use_soilc_method)
2356	CASE (SOILC_METHOD_ARITHMETIC)
2357	DO jg = 1,ngrnd
2358	DO ji = 1,kjpindex
2359	thkdry(ji,jg) = zx1(ji,jg) * cond_dry_org + zx2(ji,jg) * thkdry_min(ji)
2360	ENDDO
2361	ENDDO
2362	CASE (SOILC_METHOD_GEOMETRIC)
2363	DO jg = 1,ngrnd
2364	DO ji = 1,kjpindex
2365	! use geometric mean rather than arithmetic mean (Decharme et al 2016)
2366	thkdry(ji,jg) =(cond_dry_org*zx1(ji,jg)) (thkdry_min(ji)**zx2(ji,jg))
2367	ENDDO
2368	ENDDO
2369	CASE DEFAULT
2370	CALL ipslerr_p(3,'thermosoil_cond_nopft','Unsupported USE_SOILC_METHOD','','')
2371	ENDSELECT
2372
2373	!! 3. Thermal conductivity (cnd)
2374	DO jg = 1,ngrnd
2375	DO ji = 1,kjpindex
2376	cnd(ji,jg) = ake(ji,jg) * (thksat(ji,jg) - thkdry(ji, jg)) + thkdry(ji, jg)
2377	ENDDO
2378	END DO
2379
2380	END SUBROUTINE thermosoil_cond_nopft
2381
2382	!! ================================================================================================================================
2383	!! SUBROUTINE : thermosoil_humlev
2384	!!
2385	!>\BRIEF Interpolates the diagnostic soil humidity profile shumdiag_perma(nslm, diagnostic axis) onto
2386	!! the thermal axis, which gives shum_ngrnd_perma(ngrnd, thermal axis).
2387	!!
2388	!! DESCRIPTION : Interpolate the volumetric soil moisture content from the node to the interface of the layer.
2389	!! The values for the deep layers in thermosoil where hydrology is not existing are constant.
2390	!! No interpolation is needed for the total soil moisture content and for the soil saturation degree.
2391	!! The depths of the diagnostic levels are diaglev(1:nslm), computed in slowproc.f90.
2392	!! Recall that when the 11-layer hydrology is used,
2393	!! shum_ngrnd_perma and shumdiag_perma are with reference to the moisture content (mc)
2394	!! at the wilting point mcw : shum_ngrnd_perma=(mc-mcw)/(mcs-mcw).
2395	!! with mcs the saturated soil moisture content.
2396	!!
2397	!! RECENT CHANGE(S) : None
2398	!!
2399	!! MAIN OUTPUT VARIABLE(S): mc_layt, mcl_layt, tmc_layt, shum_ngrnd_perma (soil humidity profile on the thermal axis)
2400	!!
2401	!! REFERENCE(S) : None
2402	!!
2403	!! FLOWCHART : None
2404	!! \n
2405	!_ ================================================================================================================================
2406	SUBROUTINE thermosoil_humlev(kjpindex, shumdiag_perma, thawed_humidity,& ! in
2407	mc_layh, mcl_layh, tmc_layh, & ! in
2408	mc_layh_pft, mcl_layh_pft, tmc_layh_pft, & ! in
2409	mc_layt, mcl_layt, tmc_layt, & ! out
2410	mc_layt_pft, mcl_layt_pft, tmc_layt_pft, & ! out
2411	shum_ngrnd_perma ) ! out
2412
2413	!! 0. Variables and parameter declaration
2414
2415	!! 0.1 Input variables
2416
2417	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (unitless)
2418	REAL(r_std),DIMENSION (kjpindex,nslm), INTENT (in) :: shumdiag_perma !! Relative soil humidity on the diagnostic axis.
2419	!! (0-1, unitless). Caveats : when "hydrol" (the 11-layers
2420	!! hydrology) is used, this humidity is calculated with
2421	!! respect to the wilting point :
2422	!! shumdiag_perma= (mc-mcw)/(mcs-mcw), with mc : moisture
2423	!! content; mcs : saturated soil moisture content; mcw:
2424	!! soil moisture content at the wilting point. when the 2-layers
2425	!! hydrology "hydrolc" is used, shumdiag_perma is just
2426	!! a diagnostic humidity index, with no real physical
2427	!! meaning.
2428	REAL(r_std), DIMENSION(kjpindex), INTENT (in) :: thawed_humidity !! specified humidity of thawed soil
2429	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]
2430	REAL(r_std),DIMENSION (kjpindex,nslm), INTENT (in) :: mcl_layh !! Volumetric soil moisture content for each layer in hydrol at nodes(liquid) [m/s]
2431	REAL(r_std),DIMENSION (kjpindex,nslm), INTENT (in) :: tmc_layh !! Total soil moisture content for each layer in hydrol(liquid+ice) [mm]
2432	REAL(r_std),DIMENSION (kjpindex,nslm,nvm), INTENT (in) :: mc_layh_pft !! Volumetric soil moisture content for each layer in hydrol at nodes(liquid+ice) [m/s]
2433	REAL(r_std),DIMENSION (kjpindex,nslm,nvm), INTENT (in) :: mcl_layh_pft !! Volumetric soil moisture content for each layer in hydrol at nodes(liquid) [m/s]
2434	REAL(r_std),DIMENSION (kjpindex,nslm,nvm), INTENT (in) :: tmc_layh_pft !! Total soil moisture content for each layer in hydrol(liquid+ice) [mm]
2435
2436	!! 0.2 Output variables
2437	REAL(r_std), DIMENSION(:,:), INTENT(out):: mc_layt !! Volumetric soil moisture (liquid+ice) (m3/m3) on the thermodynamical levels at interface
2438	REAL(r_std), DIMENSION(:,:), INTENT(out):: mcl_layt !! Volumetric soil moisture (liquid) (m3/m3) on the thermodynamical levels at interface
2439	REAL(r_std), DIMENSION(:,:), INTENT(out):: tmc_layt !! Total soil moisture content for each layer (liquid+ice) (mm) on the thermodynamical levels
2440	REAL(r_std), DIMENSION(:,:,:), INTENT(out):: mc_layt_pft !! Volumetric soil moisture (liquid+ice) (m3/m3) on the thermodynamical levels at interface
2441	REAL(r_std), DIMENSION(:,:,:), INTENT(out):: mcl_layt_pft !! Volumetric soil moisture (liquid) (m3/m3) on the thermodynamical levels at interface
2442	REAL(r_std), DIMENSION(:,:,:), INTENT(out):: tmc_layt_pft !! Total soil moisture content for each layer (liquid+ice) (mm) on the thermodynamical levels
2443	REAL(r_std), DIMENSION(:,:,:), INTENT(out):: shum_ngrnd_perma !! Saturation degree on the thermal axes (0-1, dimensionless)
2444
2445	!! 0.3 Modified variables
2446
2447	!! 0.4 Local variables
2448
2449	INTEGER(i_std) :: ji, jd, jv
2450
2451	!_ ================================================================================================================================
2452
2453	shum_ngrnd_perma(:,:,:) = zero
2454	IF (printlev >= 4) WRITE(numout,*) 'Start thermosoil_humlev'
2455
2456
2457	!!!! *layt should be PFT specified, xuhui
2458	! The values for the deep layers in thermosoil where hydrology is not existing are constant.
2459	! For exemple if thermosoil uses 8m, and hydrol uses 2m vertical discretization,
2460	! the values between 2m and 8m are constant.
2461	! The moisture computed in hydrol is at the nodes (except for the
2462	! top and bottom layer which are at interfaces)
2463	! A linear interpolation is applied to obtain the moisture values at
2464	! the interfaces (mc_layt), from the mc_layh at the nodes
2465
2466	DO ji=1,kjpindex
2467	DO jd = 1, nslm
2468	IF(jd == 1) THEN ! the moisture at the 1st interface mc_layh(1) is at the surface, no interpolation
2469	mc_layt(ji,jd) = mc_layh(ji,jd)
2470	mcl_layt(ji,jd) = mcl_layh(ji,jd)
2471	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
2472	mc_layt(ji, jd) = mc_layh(ji,jd-1)*(znt(jd)-zlt(jd-1))/(znt(jd)-0.0) + &
2473	mc_layh(ji, jd)*(zlt(jd-1)-0.0)/(znt(jd)-0.0)
2474	mcl_layt(ji, jd) = mcl_layh(ji,jd-1)*(znt(jd)-zlt(jd-1))/(znt(jd)-0.0) + &
2475	mcl_layh(ji, jd)*(zlt(jd-1)-0.0)/(znt(jd)-0.0)
2476	ELSEIF(jd == nslm) THEN ! the mc_layt at the nslm interface is interpolated using mc_layh(nslm) and mc_layh(nslm-1)
2477	mc_layt(ji, jd) = mc_layh(ji,jd-1)*(zlt(jd)-zlt(jd-1))/(zlt(jd)-znt(jd-1)) + &
2478	mc_layh(ji,jd)*(zlt(jd-1)-znt(jd-1))/(zlt(jd)-znt(jd-1))
2479	mcl_layt(ji, jd) = mcl_layh(ji,jd-1)*(zlt(jd)-zlt(jd-1))/(zlt(jd)-znt(jd-1)) + &
2480	mcl_layh(ji,jd)*(zlt(jd-1)-znt(jd-1))/(zlt(jd)-znt(jd-1))
2481	ELSE ! the mc_layt at the other interfaces are interpolated using mc_layh at adjacent nodes.
2482	mc_layt(ji, jd) = mc_layh(ji, jd-1)(1-dz5(jd-1)) + mc_layh(ji,jd)dz5(jd-1)
2483	mcl_layt(ji, jd) = mcl_layh(ji, jd-1)(1-dz5(jd-1)) + mcl_layh(ji,jd)dz5(jd-1)
2484	ENDIF
2485	tmc_layt(ji,jd) = tmc_layh(ji,jd)
2486	ENDDO !jd
2487
2488	! The deep layers in thermosoil where hydro is not existing
2489	DO jd = nslm+1, ngrnd
2490	mc_layt(ji,jd) = mc_layh(ji,nslm)
2491	mcl_layt(ji,jd) = mcl_layh(ji,nslm)
2492	tmc_layt(ji,jd) = tmc_layh(ji,nslm)/dlt(nslm) *dlt(jd)
2493	ENDDO
2494	ENDDO
2495
2496	DO ji=1,kjpindex
2497	DO jd = 1, nslm
2498	DO jv = 1,nvm
2499	IF(jd == 1) THEN
2500	mc_layt_pft(ji,jd,jv) = MAX(mc_layh_pft(ji,jd,jv), min_sechiba)
2501	mcl_layt_pft(ji,jd,jv) = MAX(mcl_layh_pft(ji,jd,jv), min_sechiba)
2502	ELSEIF(jd == 2) THEN
2503	mc_layt_pft(ji,jd,jv) = MAX(mc_layh_pft(ji,jd-1,jv)*(znt(jd)-zlt(jd-1))/(znt(jd)-0.0) + &
2504	mc_layh_pft(ji,jd,jv)*(zlt(jd-1)-0.0)/(znt(jd)-0.0), min_sechiba)
2505	mcl_layt_pft(ji,jd,jv) = MAX(mcl_layh_pft(ji,jd-1,jv)*(znt(jd)-zlt(jd-1))/(znt(jd)-0.0) + &
2506	mcl_layh_pft(ji,jd,jv)*(zlt(jd-1)-0.0)/(znt(jd)-0.0), min_sechiba)
2507	ELSEIF(jd == nslm) THEN
2508	mc_layt_pft(ji,jd,jv) = MAX(mc_layh_pft(ji,jd-1,jv)*(zlt(jd)-zlt(jd-1))/(zlt(jd)-znt(jd-1)) + &
2509	mc_layh_pft(ji,jd,jv)*(zlt(jd-1)-znt(jd-1))/(zlt(jd)-znt(jd-1)),min_sechiba)
2510	mcl_layt_pft(ji,jd,jv) = MAX(mcl_layh_pft(ji,jd-1,jv)*(zlt(jd)-zlt(jd-1))/(zlt(jd)-znt(jd-1)) + &
2511	mcl_layh_pft(ji,jd,jv)*(zlt(jd-1)-znt(jd-1))/(zlt(jd)-znt(jd-1)),min_sechiba)
2512	ELSE
2513	mc_layt_pft(ji,jd,jv) = MAX(mc_layh_pft(ji,jd-1,jv)(1-dz5(jd-1)) + mc_layh_pft(ji,jd,jv)dz5(jd-1), min_sechiba)
2514	mcl_layt_pft(ji,jd,jv) = MAX(mcl_layh_pft(ji,jd-1,jv)(1-dz5(jd-1)) + mcl_layh_pft(ji,jd,jv)dz5(jd-1), min_sechiba)
2515	ENDIF
2516	ENDDO ! jv
2517	tmc_layt_pft(ji,jd,:) = tmc_layh_pft(ji,jd,:)
2518	ENDDO !jd
2519
2520	! The deep layers in thermosoil where hydro is not existing
2521	DO jd = nslm+1, ngrnd
2522	mc_layt_pft(ji,jd,:) = mc_layh_pft(ji,nslm,:)
2523	mcl_layt_pft(ji,jd,:) = mcl_layh_pft(ji,nslm,:)
2524	tmc_layt_pft(ji,jd,:) = tmc_layh_pft(ji,nslm,:)
2525	ENDDO
2526	ENDDO
2527
2528	IF (.NOT. satsoil ) THEN
2529
2530	DO jv = 1, nvm
2531
2532	! The values for the deep layers in thermosoil where hydro is not existing are constant.
2533	! For exemple if thermosoil uses 8m, and hydrol uses 2m vertical discretization,
2534	! the values between 2m and 8m are constant.
2535
2536	DO jd = 1, nslm
2537	shum_ngrnd_perma(:,jd,jv) = shumdiag_perma(:,jd)
2538	END DO
2539	DO jd = nslm+1,ngrnd
2540	shum_ngrnd_perma(:,jd,jv) = shumdiag_perma(:,nslm)
2541	!Former version of MICT before the new soil vertical discretization
2542	!see update_deep_soil_moisture
2543	! IF ( (ptn(ji,jd,jv) .GT. (ZeroCelsius+fr_dT/2.)) THEN
2544	! shum_ngrnd_perma(ji,jd,jv) = thawed_humidity(ji)
2545	! ENDIF
2546	!No else defined ??
2547	!-> Right now we stay with the TRUNK version. Possibility to add a flag later to reactivate this part as an option.
2548	END DO
2549	END DO
2550
2551	!now update the deep permafrost soil moisture separately
2552	!CALL update_deep_soil_moisture(kjpindex, shumdiag_perma,proglevel_bottomdiaglev, proglevel_zdeep, &
2553	! thawed_humidity)
2554
2555	ELSE
2556	!This is a weird option, what about the coherence with shumdiag_perma ans the hydrology in general?
2557	shum_ngrnd_perma(:,:,:) = 1.
2558	ENDIF
2559
2560	IF (printlev >= 4) WRITE(numout,*) 'thermosoil_humlev done'
2561
2562	END SUBROUTINE thermosoil_humlev
2563
2564
2565	!! ================================================================================================================================
2566	!! SUBROUTINE : thermosoil_energy_diag
2567	!!
2568	!>\BRIEF Calculate diagnostics
2569	!!
2570	!! DESCRIPTION : Calculate diagnostic variables coldcont_incr and coldcont_incr
2571	!!
2572	!! RECENT CHANGE(S) : None
2573	!!
2574	!! MAIN OUTPUT VARIABLE(S) :
2575	!!
2576	!! REFERENCE(S) : None
2577	!!
2578	!! FLOWCHART : None
2579	!! \n
2580	!_ ================================================================================================================================
2581
2582	SUBROUTINE thermosoil_energy_diag(kjpindex, temp_sol_new, soilcap, veget_max)
2583
2584	!! 0. Variables and parameter declaration
2585
2586	!! 0.1 Input variables
2587	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size
2588	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_sol_new!! New soil temperature
2589	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: soilcap !! Soil capacity
2590	REAL(r_std), DIMENSION (kjpindex,nvm), INTENT (in) :: veget_max !! Fraction of vegetation type
2591
2592	!! 0.2 Local variables
2593	INTEGER(i_std) :: ji, jg
2594
2595	!! 0.3 Modified variables
2596
2597	!! 0.4 Local variables
2598
2599	!_ ================================================================================================================================
2600	!
2601	! Sum up the energy content of all layers in the soil.
2602	!
2603	DO ji = 1, kjpindex
2604	!
2605	IF (SUM(pcapa_en(ji,1,:)*veget_max(ji,:)) .LE. sn_capa) THEN
2606	!
2607	! Verify the energy conservation in the surface layer
2608	!
2609	coldcont_incr(ji) = soilcap(ji) * (temp_sol_new(ji) - temp_sol_beg(ji))
2610	surfheat_incr(ji) = zero
2611	ELSE
2612	!
2613	! Verify the energy conservation in the surface layer
2614	!
2615	surfheat_incr(ji) = soilcap(ji) * (temp_sol_new(ji) - temp_sol_beg(ji))
2616	coldcont_incr(ji) = zero
2617	ENDIF
2618	ENDDO
2619
2620	! Save temp_sol_new to be used at next timestep
2621	temp_sol_beg(:) = temp_sol_new(:)
2622
2623	END SUBROUTINE thermosoil_energy_diag
2624
2625
2626
2627	!! ================================================================================================================================
2628	!! SUBROUTINE : thermosoil_readjust
2629	!!
2630	!>\BRIEF
2631	!!
2632	!! DESCRIPTION : Energy conservation : Correction to make sure that the same latent heat is released and
2633	!! consumed during freezing and thawing
2634	!!
2635	!! RECENT CHANGE(S) : None
2636	!!
2637	!! MAIN OUTPUT VARIABLE(S): ptn (soil temperature profile on the thermal axis),
2638	!!
2639	!! REFERENCE(S) :
2640	!!
2641	!! FLOWCHART : None
2642	!! \n
2643	!_ ================================================================================================================================
2644
2645	SUBROUTINE thermosoil_readjust(kjpindex, ptn)
2646
2647	!! 0. Variables and parameter declaration
2648
2649	!! 0.1 Input variables
2650	INTEGER(i_std), INTENT(in) :: kjpindex
2651
2652	!! 0.2 Modified variables
2653	REAL(r_std),DIMENSION(kjpindex,ngrnd,nvm),INTENT(inout) :: ptn
2654
2655	!! 0.3 Local variables
2656	INTEGER(i_std) :: ji, jg, jv
2657	INTEGER(i_std) :: lev3m !! Closest interface level to 3m
2658	REAL(r_std) :: ptn_tmp
2659
2660	! The energy is spread over the layers down to approximatly 3m
2661	! Find the closest level to 3m. It can be below or above 3m.
2662	lev3m=MINLOC(ABS(zlt(:)-3.0),dim=1)
2663	IF (printlev >= 3) WRITE(numout,*) 'In thermosoil_adjust: lev3m=',lev3m, ' zlt(lev3m)=', zlt(lev3m)
2664
2665	DO jv = 1,nvm
2666	DO jg=1, ngrnd
2667	DO ji=1, kjpindex
2668	! All soil latent energy is put into e_soil_lat(ji, 1)
2669	! because the variable soil layers make it difficult to keep track of all
2670	! layers in this version
2671	! NOTE : pcapa has unit J/K/m3 and pcappa_supp has J/K
2672	e_soil_lat(ji, jv)=e_soil_lat(ji, jv)+pcappa_supp(ji,jg,jv)*(ptn(ji,jg,jv)-ptn_beg(ji,jg,jv))
2673	ENDDO ! ji=1, kjpindex
2674	ENDDO ! jg=1, ngrnd
2675	ENDDO ! jv = 1,nvm
2676
2677	DO jv = 1,nvm
2678	DO ji=1, kjpindex
2679	IF (e_soil_lat(ji,jv).GT.min_sechiba.AND.MINVAL(ptn(ji,:,jv)).GT.ZeroCelsius+fr_dT/2.) THEN
2680	! The soil is thawed: we spread the excess of energy over the uppermost 6 levels e.g. 2.7m
2681	! Here we increase the temperatures
2682	DO jg=1,lev3m
2683	ptn_tmp=ptn(ji,jg,jv)
2684
2685	ptn(ji,jg,jv)=ptn(ji,jg,jv)+MIN(e_soil_lat(ji,jv)/pcapa(ji,jg,jv)/zlt(lev3m), 0.5)
2686	e_soil_lat(ji,jv)=e_soil_lat(ji,jv)-(ptn(ji,jg,jv)-ptn_tmp)pcapa(ji,jg,jv)dlt(jg)
2687	ENDDO ! jg=1,lev3m
2688	ELSE IF (e_soil_lat(ji,jv).LT.-min_sechiba.AND.MINVAL(ptn(ji,:,jv)).GT.ZeroCelsius+fr_dT/2.) THEN
2689	! The soil is thawed
2690	! Here we decrease the temperatures
2691	DO jg=1,lev3m
2692	ptn_tmp=ptn(ji,jg,jv)
2693	ptn(ji,jg,jv)=MAX(ZeroCelsius+fr_dT/2., ptn_tmp+e_soil_lat(ji,jv)/pcapa(ji,jg,jv)/zlt(lev3m))
2694	e_soil_lat(ji,jv)=e_soil_lat(ji,jv)+(ptn_tmp-ptn(ji,jg,jv))pcapa(ji,jg,jv)dlt(jg)
2695	ENDDO ! jg=1,6
2696	ENDIF
2697	ENDDO ! ji=1, kjpindex
2698	ENDDO ! jv = 1,nvm
2699
2700	END SUBROUTINE thermosoil_readjust
2701
2702	!-------------------------------------------------------------------
2703
2704
2705
2706	!! ================================================================================================================================
2707	!! SUBROUTINE : thermosoil_getdiff
2708	!!
2709	!>\BRIEF Computes soil and snow heat capacity and conductivity
2710	!!
2711	!! DESCRIPTION : Computation of the soil thermal properties; snow properties are also accounted for
2712	!!
2713	!! RECENT CHANGE(S) : None
2714	!!
2715	!! MAIN OUTPUT VARIABLE(S):
2716	!!
2717	!! REFERENCE(S) :
2718	!!
2719	!! FLOWCHART : None
2720	!! \n
2721	!_ ================================================================================================================================
2722	SUBROUTINE thermosoil_getdiff( kjpindex, ptn, njsc, veget_max, shum_ngrnd_permalong, &
2723	profil_froz, pcappa_supp, organic_layer_thick, soilc_total, snowrho, &
2724	snowtemp, pb, mc_layt, mc_layt_pft, tmc_layt_pft, pcapa, pcapa_en, pkappa)
2725
2726	!! 0. Variables and parameter declaration
2727
2728	!! 0.1 Input variables
2729	INTEGER(i_std),INTENT(in) :: kjpindex
2730	REAL(r_std), DIMENSION (kjpindex,nvm), INTENT (in) :: veget_max !! Fraction of vegetation type
2731	REAL(r_std),DIMENSION(kjpindex,ngrnd,nvm),INTENT(in) :: shum_ngrnd_permalong
2732	REAL(r_std), DIMENSION(kjpindex), INTENT (in) :: organic_layer_thick !! how deep is the organic soil?
2733	REAL(r_std), DIMENSION(kjpindex,ndeep,nvm), INTENT (in) :: soilc_total !! total soil carbon for use in thermal calcs
2734	INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: njsc !! Index of the dominant soil textural class in the grid cell (1-nscm, unitless)
2735	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(in) :: snowrho !! Snow density
2736	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(in) :: snowtemp !! Snow temperature (K)
2737	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: pb !! Surface presure (hPa)
2738
2739	REAL(r_std), DIMENSION (kjpindex,ngrnd,nvm), INTENT(in) :: tmc_layt_pft !! Total soil moisture content for each layer(liquid+ice) (mm)
2740	REAL(r_std), DIMENSION (kjpindex,ngrnd), INTENT(in) :: mc_layt !!
2741	REAL(r_std), DIMENSION (kjpindex,ngrnd,nvm), INTENT(in) :: mc_layt_pft !!
2742
2743	!! 0.2 Modified variables
2744	REAL(r_std),DIMENSION(kjpindex,ngrnd,nvm),INTENT(in) :: ptn !! Soil temperature profile
2745
2746	!! 0.3 Output variables
2747	REAL(r_std),DIMENSION(kjpindex,ngrnd,nvm),INTENT(out) :: pcappa_supp
2748	REAL(r_std),DIMENSION(kjpindex,ngrnd,nvm),INTENT(out) :: profil_froz
2749	REAL(r_std),DIMENSION(kjpindex,ngrnd,nvm),INTENT(out) :: pcapa, pcapa_en, pkappa
2750
2751	!! 0.3 Local variables
2752	REAL(r_std) :: xx !! Unfrozen fraction of the soil
2753	REAL(r_std) :: p
2754	REAL(r_std) :: cap_iw !! Heat capacity of ice/water mixture
2755	REAL(r_std), DIMENSION(kjpindex,ngrnd,nvm) :: so_capa_dry_net
2756	REAL(r_std) :: cond_solid_net
2757	REAL(r_std) :: so_cond_dry_net
2758	INTEGER(i_std) :: ji,jg,jv
2759	INTEGER(i_std) :: jst
2760	REAL(r_std), DIMENSION(kjpindex,ngrnd,nvm) :: poros_net
2761	REAL(r_std), DIMENSION(kjpindex,ngrnd,nvm) :: zx1, zx2
2762	REAL(r_std), DIMENSION(kjpindex,ngrnd) :: profil_froz_mean, tmp, zx1_tmp
2763
2764	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: zx1_iface, zx2_iface
2765	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: poros_net_iface
2766	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: pkappa_iface
2767	INTEGER(i_std) :: ier
2768
2769	! Organic and anorgaic layer fraction
2770	!
2771	! Default: organic layer not taken into account
2772	zx1(:,:,:) = 0.0
2773	zx2(:,:,:) = 0.0
2774	poros_net(:,:,:) = 0.0
2775	tmp(:,:) = 0.0
2776	!
2777	IF ( use_toporganiclayer_tempdiff ) THEN
2778	CALL thermosoil_toporganiclayer_tempdiff(organic_layer_thick, zlt, zx1_tmp)
2779	! Extend values to all pft's
2780	DO jv = 1, nvm
2781	DO jg = 1, ngrnd !- 2
2782	DO ji = 1,kjpindex
2783	zx1(ji,jg,jv) = zx1_tmp(ji,jg)
2784	ENDDO
2785	ENDDO
2786	ENDDO
2787
2788	ELSEIF ( use_soilc_tempdiff ) THEN
2789	!
2790	IF (use_refSOC) THEN
2791	DO jv = 1,nvm
2792	DO jg = 1, ngrnd
2793	DO ji = 1,kjpindex
2794	zx1(ji,jg,jv) = refSOC(ji,jg)/soilc_max !after lawrence and slater
2795	ENDDO
2796	ENDDO
2797	ENDDO
2798	ELSE ! use the simulated SOC(summed over PFTs)
2799	DO jv = 1,nvm
2800	DO jg = 1, ngrnd
2801	DO ji = 1,kjpindex
2802	tmp(ji,jg) = tmp(ji,jg) + soilc_total(ji,jg,jv)*veget_max(ji,jv)/soilc_max !after lawrence and slater
2803	ENDDO
2804	ENDDO
2805	ENDDO
2806	DO jg = 1, ngrnd
2807	DO ji = 1,kjpindex
2808	zx1(ji,jg,:) = tmp(ji,jg)
2809	ENDDO
2810	ENDDO
2811
2812	ENDIF
2813	!
2814	WHERE (zx1 > 1) zx1 = 1
2815	!
2816	ENDIF ! ( use_soilc_tempdiff ) THEN
2817	!
2818	zx2(:,:,:) = 1.-zx1(:,:,:)
2819
2820	DO jv = 1,nvm
2821	DO jg = 1, ngrnd
2822	DO ji = 1,kjpindex
2823	jst = njsc(ji)
2824	!
2825	! 1. Calculate dry heat capacity and conductivity, taking
2826	! into account the organic and mineral fractions in the layer
2827	!
2828	! Former MICT version
2829	!Here we take into account the new dependance of the soil heat capacity from the soil type.
2830	so_capa_dry_net(ji,jg,jv) = zx1(ji,jg,jv) * SO_CAPA_DRY_ORG + zx2(ji,jg,jv) * so_capa_dry_ns(jst)
2831
2832	!cond_solid_net = un / ( zx1(ji,jg,jv) / cond_solid_org + zx2(ji,jg,jv) / cond_solid ) ! TO DELETE
2833	!Here we take into account the new dependance of the porosity from the soil type.
2834	poros_net(ji,jg,jv) = zx1(ji,jg,jv) * poros_org + zx2(ji,jg,jv) * SMCMAX(jst)
2835	!
2836	!so_cond_dry_net = un / ( zx1(ji,jg,jv) / cond_dry_org + zx2(ji,jg,jv) / so_cond_dry ) ! TO DELETE
2837	!
2838	! 2. Calculate heat capacity with allowance for permafrost
2839	ENDDO
2840	ENDDO
2841	ENDDO
2842	!
2843	IF (ok_freeze_thermix) THEN
2844
2845	#ifdef STRICT_CHECK
2846	IF (ANY(tmc_layt_pft < min_sechiba)) CALL ipslerr_p(3, "thermosoil_getdiff", "tmc_layt_pft has negative values", "", "") ! prec issues
2847	#endif
2848	CALL thermosoil_freeze_thermix(kjpindex, ngrnd, nvm, njsc, ptn, shum_ngrnd_permalong, &
2849	mc_layt, mc_layt_pft, tmc_layt_pft, so_capa_dry_net, dlt, &
2850	profil_froz, pcapa, pcappa_supp) ! out
2851	ELSE !++cdk this is physically wrong and only to be used to test the influence of latent heat
2852	DO jv = 1,nvm
2853	DO jg = 1, ngrnd
2854	DO ji = 1,kjpindex
2855	profil_froz(ji,jg,jv) = 0.
2856
2857	IF (ok_LAIdev(jv)) THEN
2858	pcapa(ji,jg,jv) = so_capa_dry_net(ji,jg,jv) + water_capa * mc_layt_pft(ji,jg,jv)
2859	ELSE
2860	pcapa(ji,jg,jv) = so_capa_dry_net(ji,jg,jv) + water_capa * mc_layt(ji,jg)
2861	ENDIF
2862
2863	IF (brk_flag == 1) THEN
2864	! Bedrock flag is activated
2865	pcapa(ji,ngrnd-1:ngrnd,jv) = brk_capa
2866	pkappa(ji,ngrnd-1:ngrnd,jv) = brk_cond
2867	ENDIF
2868	ENDDO
2869	ENDDO
2870	ENDDO
2871	ENDIF
2872
2873	pcapa_en(:,:,:) = pcapa(:,:,:)
2874	!
2875	! 3. Calculate the heat conductivity with allowance for permafrost
2876	! Note: mc_layt has no PFT dimention,so we calculate here profil_froz_mean. Actually, profil_froz along the PFT dimention currently has no difference for each PFT.
2877	IF ( ANY(MAXVAL(profil_froz,DIM=3)>MINVAL(profil_froz,DIM=3)) ) THEN
2878	CALL ipslerr_p(3,'thermosoil_getdiff','profil_froz_mean wrong','','')
2879	ENDIF
2880	profil_froz_mean=MINVAL(profil_froz,DIM=3)
2881	tmp = mc_layt*(1.- profil_froz_mean)
2882
2883	IF (ANY(ok_LAIdev)) THEN ! CROP module
2884	CALL thermosoil_cond_pft (kjpindex, njsc, mc_layt, QZ, SMCMAX, tmp, zx1,zx2, poros_net,pkappa)
2885	ELSE
2886	ALLOCATE(zx1_iface(kjpindex,ngrnd), zx2_iface(kjpindex,ngrnd), stat=ier)
2887	IF (ier /= 0) CALL ipslerr_p(3, 'thermosoil_cond_pft', 'Allocation error for variables', 'zx1_iface and zx2_iface', '')
2888	ALLOCATE(poros_net_iface(kjpindex,ngrnd), pkappa_iface(kjpindex,ngrnd), stat=ier)
2889	IF (ier /= 0) CALL ipslerr_p(3, 'thermosoil_cond_pft', 'Allocation error for variables', 'zx1_iface and pkappa_iface', '')
2890	zx1_iface = zero
2891	zx2_iface = zero
2892	poros_net_iface = zero
2893	pkappa_iface = zero
2894
2895	! transform arrays from 3D to 2D
2896	DO jg = 1, ngrnd
2897	DO ji = 1,kjpindex
2898	zx1_iface(ji,jg) = zx1(ji,jg,1)
2899	zx2_iface(ji,jg) = zx2(ji,jg,1)
2900	poros_net_iface(ji,jg) = poros_net(ji,jg,1)
2901	ENDDO
2902	ENDDO
2903
2904	CALL thermosoil_cond_nopft (kjpindex, njsc, mc_layt, qz, smcmax, tmp, &
2905	zx1_iface, zx2_iface, poros_net_iface, pkappa_iface)
2906
2907	! Put values back to its original array
2908	DO jg = 1, ngrnd
2909	DO ji = 1,kjpindex
2910	pkappa(ji,jg,:) = pkappa_iface(ji,jg)
2911	ENDDO
2912	ENDDO
2913
2914	DEALLOCATE(zx1_iface)
2915	DEALLOCATE(zx2_iface)
2916	DEALLOCATE(poros_net_iface)
2917	DEALLOCATE(pkappa_iface)
2918	ENDIF
2919
2920	! CALL thermosoil_cond_pft (kjpindex, njsc, mc_layt_pft_tmp, QZ, poros_net, mcl_layt_pft_tmp, pkappa)
2921	!! xuhui: the above line should be activated if soil moisture budget is PFT specific
2922	! DO jv = 1,nvm
2923	! CALL thermosoil_cond_pft(kjpindex, njst, mc_layt_pft_tmp, QZ, poros_net, mcl_layt_pft_tmp, pkappa_pft(:,:,jv))
2924	! ENDDO
2925
2926	!! Computes snow heat capacity and conductivity
2927	CALL thermosoil_snowheat(kjpindex, pb, snowrho, snowtemp, pkappa_snow, pcapa_snow)
2928
2929	END SUBROUTINE thermosoil_getdiff
2930
2931	!! ================================================================================================================================
2932	!! SUBROUTINE : thermosoil_snowheat
2933	!!
2934	!>\BRIEF
2935	!!
2936	!! DESCRIPTION : Computes snow heat capacity and conductivity
2937	!!
2938	!! RECENT CHANGE(S) : None
2939	!!
2940	!! MAIN OUTPUT VARIABLE(S):
2941	!!
2942	!! REFERENCE(S) :
2943	!!
2944	!! FLOWCHART : None
2945	!! \n
2946	!_ ================================================================================================================================
2947	SUBROUTINE thermosoil_snowheat(kjpindex, pb, snowrho, snowtemp, pkappa_snow_out, pcapa_snow_out)
2948	!! 0. Variables and parameter declaration
2949
2950	!! 0.1 Input variables
2951	INTEGER(i_std),INTENT(in) :: kjpindex
2952	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(in) :: snowrho !! Snow density
2953	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(in) :: snowtemp !! Snow temperature (K)
2954	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: pb !! Surface presure (hPa)
2955
2956	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(out) :: pkappa_snow_out !!
2957	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(out) :: pcapa_snow_out !!
2958
2959	INTEGER(i_std) :: ji
2960
2961	!! Computes snow heat capacity and conductivity
2962	DO ji = 1,kjpindex
2963	pcapa_snow_out(ji,:) = snowrho(ji,:) * xci
2964
2965	SELECTCASE (snow_cond_method)
2966	CASE (SNOW_COND_METHOD_DEFAULT)
2967	pkappa_snow_out(ji,:) = (ZSNOWTHRMCOND1 + ZSNOWTHRMCOND2snowrho(ji,:)snowrho(ji,:)) + &
2968	MAX(0.0,(ZSNOWTHRMCOND_AVAP+(ZSNOWTHRMCOND_BVAP/(snowtemp(ji,:)+ &
2969	ZSNOWTHRMCOND_CVAP)))(XP00/(pb(ji)100.)))
2970	CASE (SNOW_COND_METHOD_DECHARME16)
2971	pkappa_snow_out(ji,:) = 2.2((snowrho(ji,:)/1000.)*2.0) + &
2972	MAX(0.0,(ZSNOWTHRMCOND_AVAP+(ZSNOWTHRMCOND_BVAP/(snowtemp(ji,:)+ &
2973	ZSNOWTHRMCOND_CVAP)))(XP00/(pb(ji)100.)))
2974	CASE DEFAULT
2975	CALL ipslerr_p(3,'thermosoil_getdiff','Unsupported SNOW_COND_METHOD', &
2976	'Currently supported methods are ','default(1) or ducharme16(2)')
2977	ENDSELECT
2978
2979	END DO
2980	END SUBROUTINE thermosoil_snowheat
2981
2982	!! ================================================================================================================================
2983	!! SUBROUTINE : thermosoil_freeze_thermix
2984	!!
2985	!>\BRIEF
2986	!!
2987	!! DESCRIPTION :
2988	!!
2989	!!
2990	!! RECENT CHANGE(S) : None
2991	!!
2992	!! MAIN OUTPUT VARIABLE(S): pcappa_supp, profil_froz, pcapa
2993	!!
2994	!! REFERENCE(S) :
2995	!!
2996	!! FLOWCHART : None
2997	!! \n
2998	!_ ================================================================================================================================
2999	SUBROUTINE thermosoil_freeze_thermix(kjpindex, ngrnd, nvm, njsc, ptn, shum_ngrnd_permalong, &
3000	mc_layt, mc_layt_pft, tmc_layt_pft, so_capa_dry_net, dlt, &
3001	profil_froz, pcapa, pcappa_supp) ! out
3002
3003	INTEGER(i_std), INTENT(in) :: kjpindex, ngrnd, nvm
3004	INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: njsc !! Index of the dominant soil textural class in the grid cell (1-nscm, unitless)
3005	REAL(r_std), DIMENSION(kjpindex,ngrnd,nvm),INTENT(in) :: ptn !! Soil temperature profile
3006	REAL(r_std), DIMENSION(kjpindex,ngrnd,nvm),INTENT(in) :: shum_ngrnd_permalong
3007	REAL(r_std), DIMENSION(ngrnd), INTENT(in) :: dlt !!
3008	REAL(r_std), DIMENSION(kjpindex,ngrnd), INTENT(in) :: mc_layt !!
3009	REAL(r_std), DIMENSION(kjpindex,ngrnd,nvm), INTENT(in) :: mc_layt_pft !!
3010	REAL(r_std), DIMENSION(kjpindex,ngrnd,nvm), INTENT(in) :: tmc_layt_pft !!
3011	REAL(r_std), DIMENSION(kjpindex,ngrnd,nvm) :: so_capa_dry_net
3012
3013	!! 0.3 Output variables
3014	REAL(r_std),DIMENSION(kjpindex,ngrnd,nvm),INTENT(out) :: pcappa_supp
3015	REAL(r_std),DIMENSION(kjpindex,ngrnd,nvm),INTENT(out) :: profil_froz
3016	REAL(r_std),DIMENSION(kjpindex,ngrnd,nvm),INTENT(out) :: pcapa
3017
3018	!! 0.4 Local variables
3019	INTEGER(i_std) :: ji,jg,jv
3020
3021	REAL(r_std), DIMENSION(kjpindex,ngrnd) :: pcapa_spec !! SPECIFIC soil heat capacity (J/kg/K)
3022	REAL(r_std) :: rho_tot !! Soil density (kg/m3)
3023	REAL(r_std) :: xx !! Unfrozen fraction of the soil
3024	REAL(r_std), DIMENSION(kjpindex) :: mcs_index !! Convert mcs(nscm) to mcs(kjpindex)
3025
3026	pcapa_spec = zero
3027
3028	!! Precalculate mcs_index
3029	DO ji = 1,kjpindex
3030	mcs_index(ji) = mcs(njsc(ji))
3031	ENDDO
3032
3033	DO jv = 1,nvm
3034	DO jg = 1, ngrnd
3035	DO ji = 1,kjpindex
3036	! 2.1. soil heat capacity depending on temperature and humidity
3037	IF (ptn(ji,jg,jv) .LT. ZeroCelsius-fr_dT/2.) THEN
3038	! frozen soil
3039	profil_froz(ji,jg,jv) = 1.
3040	pcappa_supp(ji,jg, jv)= 0.
3041	!! this is from Koven's version: pcapa(ji,jg,jv) = so_capa_dry_net + shum_ngrnd_permalong(ji,jg,jv)poros_net(ji,jg,jv)capa_ice*rho_ice
3042	IF (ok_LAIdev(jv)) THEN
3043	pcapa(ji,jg,jv) = so_capa_dry_net(ji,jg,jv) + so_capa_ice * mc_layt_pft(ji,jg,jv)
3044	ELSE
3045	pcapa(ji,jg,jv) = so_capa_dry_net(ji,jg,jv) * (1-mcs_index(ji)) + so_capa_ice * mc_layt(ji,jg)
3046	ENDIF
3047	rho_tot = rho_soil * (1-mcs_index(ji)) + rho_ice * tmc_layt(ji,jg) / mille / dlt(jg)
3048	pcapa_spec(ji, jg) = pcapa_spec(ji, jg) + pcapa(ji, jg, jv) / rho_tot
3049	ELSEIF (ptn(ji,jg,jv) .GT. ZeroCelsius+fr_dT/2.) THEN
3050	! unfrozen soil
3051	profil_froz(ji,jg,jv) = 0.
3052	pcappa_supp(ji,jg,jv)= 0.
3053	!! this is from Koven's version: pcapa(ji,jg,jv) = so_capa_dry_net + shum_ngrnd_permalong(ji,jg,jv)poros_net(ji,jg,jv)capa_water*rho_water
3054	IF (ok_LAIdev(jv)) THEN
3055	pcapa(ji,jg,jv) = so_capa_dry_net(ji,jg,jv) + water_capa * mc_layt_pft(ji,jg,jv)
3056	ELSE
3057	pcapa(ji,jg,jv) = so_capa_dry_net(ji,jg,jv) * (1-mcs_index(ji)) + water_capa * mc_layt(ji,jg)
3058	ENDIF
3059	rho_tot = rho_soil * (1-mcs_index(ji)) + rho_water * tmc_layt(ji,jg)/mille/dlt(jg)
3060	pcapa_spec(ji, jg) = pcapa_spec(ji, jg) + pcapa(ji, jg, jv) / rho_tot
3061	ELSE
3062
3063	pcappa_supp(ji,jg,jv)= shum_ngrnd_permalong(ji,jg,jv)lhfrho_water/fr_dT
3064	IF (jg .GT. nslm) pcappa_supp(ji,jg,jv)= 0.
3065
3066	! x is the unfrozen fraction of soil water
3067	xx = (ptn(ji,jg,jv)-(ZeroCelsius-fr_dT/2.)) / fr_dT
3068	profil_froz(ji,jg,jv) = (1. - xx)
3069	! net heat capacity of the ice/water mixture
3070	IF (ok_LAIdev(jv)) THEN
3071	pcapa(ji,jg,jv) = so_capa_dry_net(ji,jg,jv) + &
3072	& water_capa * tmc_layt_pft(ji,jg,jv)/ mille / dlt(jg) * xx + so_capa_ice * tmc_layt_pft(ji,jg,jv) / mille/dlt(jg) * (1.-xx)
3073	ELSE
3074	pcapa(ji,jg,jv) = so_capa_dry_net(ji,jg,jv) * (1-mcs_index(ji)) + &
3075	& water_capa * mc_layt(ji,jg) * xx + so_capa_ice * mc_layt(ji,jg) * (1.-xx) + pcappa_supp(ji,jg,jv)
3076	ENDIF
3077	rho_tot = rho_soil* (1-mcs_index(ji)) + &
3078	rho_water * tmc_layt(ji,jg)/mille / dlt(jg) * xx + &
3079	rho_ice * tmc_layt(ji,jg) / mille/dlt(jg) * (1.-xx)
3080	pcapa_spec(ji, jg) = pcapa_spec(ji, jg) + pcapa(ji, jg, jv) / rho_tot
3081	ENDIF
3082	ENDDO
3083	ENDDO
3084	ENDDO
3085
3086	! Output the specific heat capcaity for SP-MIP
3087	CALL xios_orchidee_send_field("pcapa_spec",pcapa_spec)
3088
3089	END SUBROUTINE thermosoil_freeze_thermix
3090
3091
3092	!! ================================================================================================================================
3093	!! SUBROUTINE : thermosoil_toporganiclayer_tempdiff
3094	!!
3095	!>\BRIEF
3096	!!
3097	!! DESCRIPTION :
3098	!!
3099	!!
3100	!! RECENT CHANGE(S) : None
3101	!!
3102	!! MAIN OUTPUT VARIABLE(S): zx1
3103	!!
3104	!! REFERENCE(S) :
3105	!!
3106	!! FLOWCHART : None
3107	!! \n
3108	!_ ================================================================================================================================
3109	SUBROUTINE thermosoil_toporganiclayer_tempdiff(organic_layer_thick, zlt, zx1)
3110	! Arguments
3111	REAL(r_std), DIMENSION(:), INTENT(in) :: zlt ! ngrnd
3112	REAL(r_std), DIMENSION(:), INTENT (in) :: organic_layer_thick !! how deep is the organic soil? kjpindex
3113	REAL(r_std), DIMENSION(:,:), INTENT(OUT) :: zx1 ! kpjindex, ngrnd
3114
3115	! Local
3116	INTEGER(i_std) :: jg, ji, kjpindex
3117
3118	kjpindex = SIZE(zx1, DIM=1)
3119	ngrnd = SIZE(zx1, DIM=2)
3120	!
3121	! level 1
3122	!
3123	DO ji = 1,kjpindex
3124	IF ( organic_layer_thick(ji) .GT. zlt(1) ) THEN
3125	!! the 1st level is in the organic => the 1st layer is entirely organic
3126	zx1(ji,1) = 1. !!zx1 being the fraction of each level that is organic, zx2 is the remainder
3127	ELSE IF ( organic_layer_thick(ji) .GT. zero ) THEN
3128	!! the 1st level is beyond the organic and the organic is present
3129	zx1(ji,1) = organic_layer_thick(ji) / zlt(1)
3130	ELSE
3131	! there is no organic at all
3132	zx1(ji,1) = 0.
3133	ENDIF
3134	ENDDO
3135	!
3136	! other levels
3137	!
3138	DO jg = 2, ngrnd !- 2
3139	DO ji = 1,kjpindex
3140	IF ( organic_layer_thick(ji) .GT. zlt(jg) ) THEN
3141	! the current level is in the organic => the current layer is
3142	! entirely organic
3143	zx1(ji,jg) = 1.
3144	ELSE IF ( organic_layer_thick(ji) .GT. zlt(jg-1) ) THEN
3145	! the current layer is partially organic
3146	zx1(ji,jg) = (organic_layer_thick(ji) - zlt(jg-1)) / (zlt(jg) - zlt(jg-1))
3147	ELSE
3148	! both levels are out of organic => the current layer is entirely
3149	! mineral soil
3150	zx1(ji,jg) = 0.
3151	ENDIF
3152	ENDDO
3153	ENDDO
3154
3155	END SUBROUTINE thermosoil_toporganiclayer_tempdiff
3156
3157	!! ================================================================================================================================
3158	!! SUBROUTINE : thermosoil_getdiff_old_thermix_with_snow
3159	!!
3160	!>\BRIEF Computes soil heat capacity and conductivity
3161	!!
3162	!! DESCRIPTION : Computes soil heat capacity and conductivity
3163	!! Special case with old snow without soil freezing
3164	!!
3165	!! RECENT CHANGE(S) : None
3166	!!
3167	!! MAIN OUTPUT VARIABLE(S):
3168	!!
3169	!! REFERENCE(S) :
3170	!!
3171	!! FLOWCHART : None
3172	!! \n
3173	!_ ================================================================================================================================
3174
3175
3176	SUBROUTINE thermosoil_getdiff_old_thermix_with_snow( kjpindex, snow, njsc )
3177
3178
3179	!! 0. Variables and parameter declaration
3180
3181	!! 0.1 Input variables
3182	INTEGER(i_std), INTENT(in) :: kjpindex
3183	REAL(r_std),DIMENSION(kjpindex),INTENT (in) :: snow
3184	INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: njsc !! Index of the dominant soil textural class in the grid cell (1-nscm, unitless)
3185
3186
3187	!! 0.2 Local variables
3188	INTEGER :: ji,jg,jv
3189	REAL(r_std) :: snow_h !! snow_h is the snow height @tex ($m$) @endtex
3190	REAL(r_std) :: zx1, zx2 !! zx1 and zx2 are the layer fraction consisting in snow and soil respectively.
3191	INTEGER :: jst
3192	REAL(r_std), DIMENSION (kjpindex,ngrnd) :: pkappa_tmp !! soil thermal conductivity (W/m/K)
3193	REAL(r_std), DIMENSION (kjpindex,ngrnd) :: pcapa_tmp !! soil heat capacity (J/m3/K)
3194	REAL(r_std), DIMENSION (kjpindex,ngrnd,nvm) :: pkappa_pft_tmp !! soil thermal conductivity (W/m/K)
3195	REAL(r_std), DIMENSION (kjpindex,ngrnd,nvm) :: pcapa_pft_tmp !! soil heat capacity (J/m3/K)
3196	REAL(r_std), DIMENSION (kjpindex,ngrnd) :: pkappa_wet !! wet soil thermal conductivity (W/m/K)
3197	REAL(r_std), DIMENSION (kjpindex,ngrnd) :: pcapa_wet !! wet soil heat capacity (J/m3/K)
3198	REAL(r_std), DIMENSION (kjpindex,ngrnd) :: mc_layt_tmp !! volumetric soil moisture (liquid+ice) (m/m3)
3199	REAL(r_std), DIMENSION (kjpindex,ngrnd) :: mcl_layt_tmp !! volumetric soil moisture (liquid) (m/m3)
3200	REAL(r_std), DIMENSION (kjpindex,ngrnd) :: tmc_layt_tmp !! total soil moisture content for each layer, mm
3201	REAL(r_std), DIMENSION (kjpindex,ngrnd,nvm) :: mc_layt_pft_tmp !! volumetric soil moisture (liquid+ice) (m/m3)
3202	REAL(r_std), DIMENSION (kjpindex,ngrnd,nvm) :: mcl_layt_pft_tmp !! volumetric soil moisture (liquid) (m/m3)
3203	REAL(r_std), DIMENSION (kjpindex,ngrnd,nvm) :: tmc_layt_pft_tmp !! total soil moisture content for each layer, mm
3204
3205	REAL(r_std), DIMENSION (kjpindex,ngrnd) :: mcs_tmp !! Saturated soil moisture (liquid+ice) (m3/m3)
3206
3207	! Computation of the soil thermal properties; snow properties are also accounted for
3208
3209	pkappa_tmp(:,:) = 0.0
3210	pcapa_tmp(:,:) = 0.0
3211	pkappa_pft_tmp(:,:,:) = 0.0
3212	pcapa_pft_tmp(:,:,:) = 0.0
3213	pkappa_wet(:,:) = 0.0
3214	pcapa_wet(:,:) = 0.0
3215	DO ji = 1, kjpindex
3216	jst = njsc(ji)
3217	mcs_tmp(ji,:) = mcs(jst)
3218	DO jg = 1, ngrnd
3219	pcapa_tmp(ji, jg) = so_capa_dry_ns(jst) * (1-mcs(jst)) + water_capa * tmc_layt(ji,jg)/mille/dlt(jg)
3220	pcapa_wet(ji, jg) = so_capa_dry_ns(jst) * (1-mcs(jst)) + water_capa * mcs(jst)
3221	DO jv = 1,nvm
3222	mc_layt_pft_tmp(ji, jg, jv) = MAX(mc_layt_pft(ji,jg,jv), min_sechiba)
3223	mcl_layt_pft_tmp(ji, jg, jv) = MAX(mcl_layt_pft(ji,jg,jv), min_sechiba)
3224	! CALL thermosoil_cond(kjpindex, njsc, mc_layt_pft_tmp(:,:, jv), QZ, &
3225	! SMCMAX, mcl_layt_pft_tmp(:, :, jv), pkappa_pft_tmp(:, :, jv))
3226	tmc_layt_pft_tmp(ji, jg, jv) = MAX(tmc_layt_pft(ji,jg,jv), min_sechiba)
3227	IF (ok_LAIdev(jv)) THEN
3228	pcapa_pft_tmp(ji, jg, jv) = so_capa_dry_ns(jst) + water_capa * tmc_layt_pft_tmp(ji,jg,jv)/mille/dlt(jg)
3229	ELSE
3230	pcapa_pft_tmp(ji, jg, jv) = pcapa_tmp(ji,jg)
3231	ENDIF
3232	ENDDO
3233	ENDDO
3234	ENDDO
3235
3236	DO jv = 1,nvm
3237	CALL thermosoil_cond(kjpindex, njsc, mc_layt_pft_tmp(:,:, jv), QZ, &
3238	SMCMAX, mcl_layt_pft_tmp(:, :, jv), pkappa_pft_tmp(:, :, jv))
3239	ENDDO
3240	CALL thermosoil_cond(kjpindex, njsc, mc_layt, QZ, SMCMAX, mcl_layt, pkappa_tmp)
3241	CALL thermosoil_cond(kjpindex, njsc, mcs_tmp, QZ, SMCMAX, mcs_tmp, pkappa_wet)
3242
3243	DO ji = 1,kjpindex
3244	snow_h = snow(ji) / sn_dens
3245
3246	! First layer
3247	IF ( snow_h .GT. zlt(1) ) THEN
3248	pcapa(ji,1,:) = sn_capa
3249	pcapa_en(ji,1,:) = sn_capa
3250	pkappa(ji,1,:) = sn_cond
3251	ELSE IF ( snow_h .GT. zero ) THEN
3252	pcapa_en(ji,1,:) = sn_capa
3253	zx1 = snow_h / zlt(1)
3254	zx2 = ( zlt(1) - snow_h) / zlt(1)
3255	pcapa(ji,1,:) = zx1 * sn_capa + zx2 * pcapa_wet(ji,1)
3256	pkappa(ji,1,:) = un / ( zx1 / sn_cond + zx2 / (pkappa_wet(ji,1)) )
3257	ELSE
3258	DO jv = 1,nvm
3259	IF (ok_LAIdev(jv)) THEN
3260	! pkappa(ji,1,:) = pkappa_tmp(ji,1)
3261	! pcapa(ji,1,:) = pcapa_tmp(ji,1)
3262	! pcapa_en(ji,1,:) = pcapa_tmp(ji,1)
3263	pkappa(ji,1,jv) = pkappa_pft_tmp(ji,1,jv)
3264	pcapa(ji,1,jv) = pcapa_pft_tmp(ji,1,jv)
3265	pcapa_en(ji,1,jv) = pcapa_pft_tmp(ji,1,jv)
3266	ELSE
3267	pkappa(ji,1,jv) = pkappa_tmp(ji,1)
3268	pcapa(ji,1,jv) = pcapa_tmp(ji,1)
3269	pcapa_en(ji,1,jv) = pcapa_tmp(ji,1)
3270	ENDIF
3271	ENDDO
3272	ENDIF
3273
3274	! Mid layers
3275	DO jg = 2, ngrnd - 2
3276	IF ( snow_h .GT. zlt(jg) ) THEN
3277	pcapa(ji,jg,:) = sn_capa
3278	pkappa(ji,jg,:) = sn_cond
3279	pcapa_en(ji,jg,:) = sn_capa
3280	ELSE IF ( snow_h .GT. zlt(jg-1) ) THEN
3281	zx1 = (snow_h - zlt(jg-1)) / (zlt(jg) - zlt(jg-1))
3282	zx2 = ( zlt(jg) - snow_h) / (zlt(jg) - zlt(jg-1))
3283	pcapa_en(ji,jg,:) = sn_capa
3284	pcapa(ji, jg,:) = zx1 * sn_capa + zx2 * pcapa_wet(ji,jg)
3285	pkappa(ji,jg,:) = un / ( zx1 / sn_cond + zx2 / (pkappa_wet(ji,jg)))
3286	ELSE
3287	DO jv = 1,nvm
3288	IF (ok_LAIdev(jv)) THEN
3289	! pcapa(ji,jg,:) = pcapa_tmp(ji, jg)
3290	! pkappa(ji,jg,:) = pkappa_tmp(ji,jg)
3291	! pcapa_en(ji,jg,:) = pcapa_tmp(ji, jg)
3292	pcapa(ji,jg,jv) = pcapa_pft_tmp(ji, jg, jv)
3293	pkappa(ji,jg,jv) = pkappa_pft_tmp(ji,jg, jv)
3294	pcapa_en(ji,jg,jv) = pcapa_pft_tmp(ji, jg, jv)
3295	ELSE
3296	pcapa(ji,jg,jv) = pcapa_tmp(ji, jg)
3297	pkappa(ji,jg,jv) = pkappa_tmp(ji,jg)
3298	pcapa_en(ji,jg,jv) = pcapa_tmp(ji, jg)
3299	ENDIF
3300	ENDDO
3301	ENDIF
3302	ENDDO
3303
3304	! Last two layers: These layers can not be filled with snow
3305	DO jg = ngrnd - 1, ngrnd
3306	pcapa(ji,jg,:) = so_capa_dry
3307	pkappa(ji,jg,:) = so_cond_dry
3308	pcapa_en(ji,jg,:) = so_capa_dry
3309	END DO
3310
3311	IF (brk_flag == 1) THEN
3312	! Bedrock flag is activated
3313	DO jg = ngrnd-1,ngrnd
3314	pcapa(ji,jg,:) = brk_capa
3315	pcapa_en(ji,jg,:) = brk_capa
3316	pkappa(ji,jg,:) = brk_cond
3317	ENDDO
3318	ENDIF
3319
3320	ENDDO ! DO ji = 1,kjpindex
3321
3322
3323	END SUBROUTINE thermosoil_getdiff_old_thermix_with_snow
3324
3325
3326	!! ================================================================================================================================
3327	!! SUBROUTINE : thermosoil_read_reftempfile
3328	!!
3329	!>\BRIEF
3330	!!
3331	!! DESCRIPTION : Read file with longterm soil temperature
3332	!!
3333	!!
3334	!! RECENT CHANGE(S) : None
3335	!!
3336	!! MAIN OUTPUT VARIABLE(S): reftemp : Reference temerature
3337	!!
3338	!! REFERENCE(S) :
3339	!!
3340	!! FLOWCHART : None
3341	!! \n
3342	!_ ================================================================================================================================
3343	SUBROUTINE thermosoil_read_reftempfile(kjpindex,lalo,reftemp)
3344
3345	USE interpweight
3346
3347	IMPLICIT NONE
3348
3349	!! 0. Variables and parameter declaration
3350
3351	!! 0.1 Input variables
3352	INTEGER(i_std), INTENT(in) :: kjpindex
3353	REAL(r_std), DIMENSION(kjpindex,2), INTENT(in) :: lalo
3354
3355	!! 0.2 Output variables
3356	REAL(r_std), DIMENSION(kjpindex, ngrnd), INTENT(out) :: reftemp
3357
3358	!! 0.3 Local variables
3359	INTEGER(i_std) :: ib
3360	CHARACTER(LEN=80) :: filename
3361	REAL(r_std),DIMENSION(kjpindex) :: reftemp_file !! Horizontal temperature field interpolated from file [C]
3362	INTEGER(i_std),DIMENSION(kjpindex,8) :: neighbours
3363	REAL(r_std) :: vmin, vmax !! min/max values to use for the
3364	!! renormalization
3365	REAL(r_std), DIMENSION(kjpindex) :: areftemp !! Availability of data for the interpolation
3366	CHARACTER(LEN=80) :: variablename !! Variable to interpolate
3367	!! the file
3368	CHARACTER(LEN=80) :: lonname, latname !! lon, lat names in input file
3369	REAL(r_std), DIMENSION(:), ALLOCATABLE :: variabletypevals !! Values for all the types of the variable
3370	!! (variabletypevals(1) = -un, not used)
3371	CHARACTER(LEN=50) :: fractype !! method of calculation of fraction
3372	!! 'XYKindTime': Input values are kinds
3373	!! of something with a temporal
3374	!! evolution on the dx*dy matrix'
3375	LOGICAL :: nonegative !! whether negative values should be removed
3376	CHARACTER(LEN=50) :: maskingtype !! Type of masking
3377	!! 'nomask': no-mask is applied
3378	!! 'mbelow': take values below maskvals(1)
3379	!! 'mabove': take values above maskvals(1)
3380	!! 'msumrange': take values within 2 ranges;
3381	!! maskvals(2) <= SUM(vals(k)) <= maskvals(1)
3382	!! maskvals(1) < SUM(vals(k)) <= maskvals(3)
3383	!! (normalized by maskvals(3))
3384	!! 'var': mask values are taken from a
3385	!! variable inside the file (>0)
3386	REAL(r_std), DIMENSION(3) :: maskvals !! values to use to mask (according to
3387	!! `maskingtype')
3388	CHARACTER(LEN=250) :: namemaskvar !! name of the variable to use to mask
3389	REAL(r_std) :: reftemp_norefinf
3390	REAL(r_std) :: reftemp_default !! Default value
3391
3392
3393	!Config Key = SOIL_REFTEMP_FILE
3394	!Config Desc = File with climatological soil temperature
3395	!Config If = READ_REFTEMP
3396	!Config Def = reftemp.nc
3397	!Config Help =
3398	!Config Units = [FILE]
3399	filename = 'reftemp.nc'
3400	CALL getin_p('REFTEMP_FILE',filename)
3401
3402	variablename = 'temperature'
3403
3404	IF (printlev >= 1) WRITE(numout,*) "thermosoil_read_reftempfile: Read and interpolate file " &
3405	// TRIM(filename) //" for variable " //TRIM(variablename)
3406
3407	! For this case there are not types/categories. We have 'only' a continuos field
3408	! Assigning values to vmin, vmax
3409
3410	vmin = 0.
3411	vmax = 9999.
3412
3413	! For this file we do not need neightbours!
3414	neighbours = 0
3415
3416	!! Variables for interpweight
3417	! Type of calculation of cell fractions
3418	fractype = 'default'
3419	! Name of the longitude and latitude in the input file
3420	lonname = 'nav_lon'
3421	latname = 'nav_lat'
3422	! Default value when no value is get from input file
3423	reftemp_default = 1.
3424	! Reference value when no value is get from input file
3425	reftemp_norefinf = 1.
3426	! Should negative values be set to zero from input file?
3427	nonegative = .FALSE.
3428	! Type of mask to apply to the input data (see header for more details)
3429	maskingtype = 'nomask'
3430	! Values to use for the masking (here not used)
3431	maskvals = (/ undef_sechiba, undef_sechiba, undef_sechiba /)
3432	! Name of the variable with the values for the mask in the input file (only if maskkingtype='var') (here not used)
3433	namemaskvar = ''
3434
3435	CALL interpweight_2Dcont(kjpindex, 0, 0, lalo, resolution, neighbours, &
3436	contfrac, filename, variablename, lonname, latname, vmin, vmax, nonegative, maskingtype, &
3437	maskvals, namemaskvar, -1, fractype, reftemp_default, reftemp_norefinf, &
3438	reftemp_file, areftemp)
3439	IF (printlev >= 5) WRITE(numout,*)' thermosoil_read_reftempfile after interpweight_2Dcont'
3440
3441	! Copy reftemp_file temperature to all ground levels and transform into Kelvin
3442	DO ib=1, kjpindex
3443	reftemp(ib, :) = reftemp_file(ib)+ZeroCelsius
3444	END DO
3445
3446	! Write diagnostics
3447	CALL xios_orchidee_send_field("areftemp",areftemp)
3448
3449	END SUBROUTINE thermosoil_read_reftempfile
3450
3451	!! ================================================================================================================================
3452	!! SUBROUTINE : read_refSOCfile
3453	!!
3454	!>\BRIEF
3455	!!
3456	!! DESCRIPTION : Read file of soil organic carbon to be used in thermix
3457	!! (insulating effect)
3458	!!
3459	!!
3460	!! RECENT CHANGE(S) : None
3461	!!
3462	!! MAIN OUTPUT VARIABLE(S): refSOC : soil organic carbon from data
3463	!!
3464	!! REFERENCE(S) :
3465	!!
3466	!! FLOWCHART : None
3467	!! \n
3468	!_ ================================================================================================================================
3469
3470	SUBROUTINE read_refSOCfile(nbpt, lalo, neighbours, resolution, contfrac)
3471
3472	!! 0. Variable and parameter declaration
3473
3474	!! 0.1 Input variables
3475
3476	INTEGER(i_std), INTENT(in) :: nbpt !! Number of points for which the data needs to be interpolated (unitless)
3477	REAL(r_std), INTENT(in) :: lalo(nbpt,2) !! Vector of latitude and longitudes (degree)
3478	INTEGER(i_std), INTENT(in) :: neighbours(nbpt,NbNeighb)!! Vector of neighbours for each grid point (1=N,2=E,3=S,4=W)
3479	REAL(r_std), INTENT(in) :: resolution(nbpt,2) !! The size of each grid cell in X and Y (km)
3480	REAL(r_std), INTENT(in) :: contfrac(nbpt) !! Fraction of land in each grid cell (unitless)
3481
3482	!! 0.4 Local variables
3483
3484	INTEGER(i_std) :: nbvmax !! nbvmax for interpolation (unitless)
3485	CHARACTER(LEN=80) :: filename
3486	INTEGER(i_std) :: iml, jml, lml, tml !! Indices
3487	INTEGER(i_std) :: fid, ib, ip, jp, fopt !! Indices
3488	INTEGER(i_std) :: ilf, ks !! Indices
3489	REAL(r_std) :: totarea !! Help variable to compute average SOC
3490	REAL(r_std), ALLOCATABLE, DIMENSION(:) :: lat_lu, lon_lu !! Latitudes and longitudes read from input file
3491	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: lat_rel, lon_rel !! Help variable to read file data and allocate memory
3492	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: mask_lu !! Help variable to read file data and allocate memory
3493	INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:) :: mask
3494	REAL(r_std), ALLOCATABLE, DIMENSION(:,:,:) :: refSOC_file !! Help variable to read file data and allocate memory
3495	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: sub_area !! Help variable to read file data and allocate memory
3496	INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:,:) :: sub_index !! Help variable to read file data and allocate memory
3497	CHARACTER(LEN=30) :: callsign !! Help variable to read file data and allocate memory
3498	CHARACTER(LEN=100) :: str !! Temporary string var
3499	LOGICAL :: ok_interpol !! Optional return of aggregate_2d
3500	INTEGER :: ALLOC_ERR !! Help varialbe to count allocation error
3501	!_
3502	!================================================================================================================================
3503
3504	!! 1. Open file and allocate memory
3505
3506	! Open file with SOC map
3507
3508	!Config Key = SOIL_REFSOC_FILE
3509	!Config Desc = File with climatological soil temperature
3510	!Config If = READ_REFTEMP
3511	!Config Def = reftemp.nc
3512	!Config Help =
3513	!Config Units = [FILE]
3514	filename = 'refSOC.nc'
3515	CALL getin_p('SOIL_REFSOC_FILE',filename)
3516
3517	! Read data from file
3518	IF (is_root_prc) CALL flininfo(filename, iml, jml, lml, tml, fid)
3519	CALL bcast(iml)
3520	CALL bcast(jml)
3521	CALL bcast(lml)
3522	CALL bcast(tml)
3523
3524	IF (lml .NE. ngrnd) THEN
3525	WRITE(str, *) 'ngrnd=', ngrnd, ', depth found in file=', lml
3526	CALL ipslerr_p(3, 'read_refSOCfile', &
3527	'depth from the file must be the same as ngrnd', &
3528	str, &
3529	filename )
3530	ENDIF
3531
3532	ALLOCATE(lon_lu(iml), STAT=ALLOC_ERR)
3533	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'read_refSOCfile','Problem in allocation of variable lon_lu','','')
3534
3535	ALLOCATE(lat_lu(jml), STAT=ALLOC_ERR)
3536	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'read_refSOCfile','Problem in allocation of variable lat_lu','','')
3537
3538	ALLOCATE(mask_lu(iml,jml), STAT=ALLOC_ERR)
3539	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'read_refSOCfile','Pb in allocation for mask_lu','','')
3540
3541	ALLOCATE(refSOC_file(iml,jml,lml), STAT=ALLOC_ERR)
3542	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'read_refSOCfile','Pb in allocation for refSOC_file','','')
3543
3544	IF (is_root_prc) THEN
3545	CALL flinget(fid, 'longitude', iml, 0, 0, 0, 1, 1, lon_lu)
3546	CALL flinget(fid, 'latitude', jml, 0, 0, 0, 1, 1, lat_lu)
3547	CALL flinget(fid, 'mask', iml, jml, 0, 0, 1, 1, mask_lu)
3548	CALL flinget(fid, 'soil_organic_carbon', iml, jml, lml, tml, 1, 1, refSOC_file)
3549
3550	CALL flinclo(fid)
3551	ENDIF
3552
3553	CALL bcast(lon_lu)
3554	CALL bcast(lat_lu)
3555	CALL bcast(mask_lu)
3556	CALL bcast(refSOC_file)
3557
3558	! Check for Nan values
3559	IF (ANY(refSOC_file .NE. refSOC_file)) THEN
3560	CALL ipslerr_p(3,'read_refSOCfile','Filename:'//filename, &
3561	'variable: soil_organic_carbon', &
3562	'Nan values not allowed. Check the input data')
3563	ENDIF
3564
3565	ALLOCATE(lon_rel(iml,jml), STAT=ALLOC_ERR)
3566	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'read_refSOCfile','Pb in allocation for lon_rel','','')
3567
3568	ALLOCATE(lat_rel(iml,jml), STAT=ALLOC_ERR)
3569	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'read_refSOCfile','Pb in allocation for lat_rel','','')
3570
3571	ALLOCATE(mask(iml,jml), STAT=ALLOC_ERR)
3572	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'read_refSOCfile','Problem in allocation of variable mask','','')
3573
3574	DO jp=1,jml
3575	lon_rel(:,jp) = lon_lu(:)
3576	ENDDO
3577	DO ip=1,iml
3578	lat_rel(ip,:) = lat_lu(:)
3579	ENDDO
3580
3581	mask(:,:) = zero
3582	WHERE (mask_lu(:,:) > zero )
3583	mask(:,:) = un
3584	ENDWHERE
3585
3586	! Set nbvmax to 200 for interpolation
3587	! This number is the dimension of the variables in which we store
3588	! the list of points of the source grid which fit into one grid box of the
3589	! target.
3590	nbvmax = 16
3591	callsign = 'soil organic carbon'
3592
3593	! Start interpolation
3594	ok_interpol=.FALSE.
3595	DO WHILE ( .NOT. ok_interpol )
3596	WRITE(numout,*) "Projection arrays for ",callsign," : "
3597	WRITE(numout,*) "nbvmax = ",nbvmax
3598
3599	ALLOCATE(sub_area(nbpt,nbvmax), STAT=ALLOC_ERR)
3600	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'read_refSOCfile','Pb in allocation for sub_area','','')
3601	sub_area(:,:)=zero
3602
3603	ALLOCATE(sub_index(nbpt,nbvmax,2), STAT=ALLOC_ERR)
3604	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'read_refSOCfile','Pb in allocation for sub_index','','')
3605	sub_index(:,:,:)=0
3606
3607	CALL aggregate_p(nbpt, lalo, neighbours, resolution, contfrac, &
3608	iml, jml, lon_rel, lat_rel, mask, callsign, &
3609	nbvmax, sub_index, sub_area, ok_interpol)
3610
3611	IF ( .NOT. ok_interpol ) THEN
3612	DEALLOCATE(sub_area)
3613	DEALLOCATE(sub_index)
3614	nbvmax = nbvmax * 2
3615	ENDIF
3616	ENDDO
3617
3618	! Compute the average
3619	refSOC(:,:) = zero
3620	DO ib = 1, nbpt
3621	fopt = COUNT(sub_area(ib,:) > zero)
3622	IF ( fopt > 0 ) THEN
3623	totarea = zero
3624	DO ilf = 1, fopt
3625	ip = sub_index(ib,ilf,1)
3626	jp = sub_index(ib,ilf,2)
3627	refSOC(ib,:) = refSOC(ib,:) + refSOC_file(ip,jp,:) * sub_area(ib,ilf)
3628	totarea = totarea + sub_area(ib,ilf)
3629	ENDDO
3630	! Normalize
3631	refSOC(ib,:) = refSOC(ib,:)/totarea
3632	ELSE
3633	! Set defalut value for points where the interpolation fail
3634	WRITE(numout,*) 'On point ', ib, ' no points were found for interpolation data. Mean value is used.'
3635	WRITE(numout,*) 'Location : ', lalo(ib,2), lalo(ib,1)
3636	refSOC(ib,:) = 0.
3637	ENDIF
3638	ENDDO
3639
3640	DEALLOCATE (lat_lu)
3641	DEALLOCATE (lat_rel)
3642	DEALLOCATE (lon_lu)
3643	DEALLOCATE (lon_rel)
3644	DEALLOCATE (mask_lu)
3645	DEALLOCATE (mask)
3646	DEALLOCATE (refSOC_file)
3647	DEALLOCATE (sub_area)
3648	DEALLOCATE (sub_index)
3649
3650	END SUBROUTINE read_refSOCfile
3651
3652
3653	!!
3654	!================================================================================================================================
3655	!! SUBROUTINE : add_heat_Zimov
3656	!!
3657	!>\BRIEF heat
3658	!!
3659	!! DESCRIPTION :
3660	!!
3661	!! RECENT CHANGE(S) : None
3662	!!
3663	!! MAIN OUTPUT VARIABLE(S):
3664	!!
3665	!! REFERENCE(S) :
3666	!!
3667	!! FLOWCHART : None
3668	!! \n
3669	!_
3670	!================================================================================================================================
3671	SUBROUTINE add_heat_Zimov(kjpindex, veget_max_bg, ptn, heat_Zimov)
3672	!! 0. Variables and parameter declaration
3673
3674	!! 0.1 Input variables
3675	INTEGER(i_std),INTENT(in) :: kjpindex
3676	REAL(r_std), DIMENSION (kjpindex,nvm), INTENT(in) :: veget_max_bg !! Fraction of vegetation type
3677
3678	REAL(r_std), DIMENSION(kjpindex,ndeep,nvm), INTENT (in) :: heat_Zimov !! heating associated with decomposition
3679
3680	!! 0.2 Modified variables
3681	REAL(r_std),DIMENSION(kjpindex,ngrnd,nvm),INTENT(inout) :: ptn
3682
3683	!! 0.3 Local variables
3684	INTEGER(r_std) :: ji, jg, jv
3685
3686	IF (printlev>=3) WRITE (numout,*) 'entering add_heat_Zimov'
3687
3688	DO ji = 1, kjpindex
3689	DO jv = 1,nvm
3690	DO jg = 1, ngrnd
3691	ptn(ji,jg,jv) = ptn(ji,jg,jv) + heat_zimov(ji,jg,jv) * dt_sechiba / ( pcapa(ji,jg,jv) * dlt(jg) )
3692	END DO
3693	END DO
3694	END DO
3695
3696	! ptn_pftmean needs to be updated to ensure consistency
3697	ptn_pftmean(:,:) = zero
3698	DO jv=1,nvm
3699	DO jg = 1, ngrnd
3700	ptn_pftmean(:,jg) = ptn_pftmean(:,jg) + ptn(:,jg,jv) * veget_max_bg(:,jv)
3701	ENDDO ! jg = 1, ngrnd
3702	ENDDO ! m=1,nvm
3703
3704	IF (printlev>=3) WRITE (numout,*) ' add_heat_Zimov done'
3705
3706	END SUBROUTINE add_heat_Zimov
3707
3708	!!
3709	!! ================================================================================================================================
3710	!! SUBROUTINE : thermosoil_diaglev
3711	!!
3712	!>\BRIEF Interpolation of the soil in-depth temperatures onto the diagnostic profile.
3713	!!
3714	!! DESCRIPTION : This is a very easy linear interpolation, with intfact(jd, jg) the fraction
3715	!! the thermal layer jg comprised within the diagnostic layer jd. The depths of
3716	!! the diagnostic levels are diaglev(1:nslm), computed in slowproc.f90.
3717	!!
3718	!! RECENT CHANGE(S) : None
3719	!!
3720	!! MAIN OUTPUT VARIABLE(S): stempdiag (soil temperature profile on the diagnostic axis)
3721	!!
3722	!! REFERENCE(S) : None
3723	!!
3724	!! FLOWCHART : None
3725	!! \n
3726	!_ ================================================================================================================================
3727	SUBROUTINE thermosoil_diaglev(kjpindex, stempdiag, veget_max)
3728
3729	!! 0. Variables and parameter declaration
3730
3731	!! 0.1 Input variables
3732
3733	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (unitless)
3734	REAL(r_std), DIMENSION (kjpindex,nvm), INTENT (in) :: veget_max !! Fraction of vegetation type
3735	!! 0.2 Output variables
3736
3737	REAL(r_std),DIMENSION (kjpindex,nslm), INTENT (out) :: stempdiag !! Diagnostoc soil temperature profile @tex ($K$) @endtex
3738
3739	!! 0.3 Modified variables
3740
3741	!! 0.4 Local variables
3742	INTEGER(i_std) :: jg,jv
3743	REAL(r_std), DIMENSION (kjpindex,nvm) :: veget_max_bg !! Fraction of vegetation type
3744
3745	!_ ================================================================================================================================
3746	veget_max_bg(:,2:nvm) = veget_max(:,2:nvm)
3747	veget_max_bg(:,1) = MAX((un - SUM(veget_max(:,2:nvm), 2)), zero)
3748
3749	stempdiag(:,:) = 0.
3750	DO jg = 1, nslm
3751	DO jv = 1, nvm
3752	stempdiag(:,jg) = stempdiag(:,jg) + ptn(:,jg,jv)*veget_max_bg(:,jv)
3753	ENDDO
3754	ENDDO
3755
3756	END SUBROUTINE thermosoil_diaglev
3757
3758	!================================================================================================================================
3759	!! SUBROUTINE : update_deep_soil_moisture
3760	!!
3761	!>\BRIEF updating deep soil moisture
3762	!!
3763	!! DESCRIPTION :
3764	!!
3765	!! RECENT CHANGE(S) : None
3766	!!
3767	!! MAIN OUTPUT VARIABLE(S):
3768	!!
3769	!! REFERENCE(S) : None
3770	!!
3771	!! FLOWCHART : None
3772	!! \n
3773	!_
3774	!================================================================================================================================
3775	SUBROUTINE update_deep_soil_moisture (kjpindex, shumdiag_perma, proglevel_bottomdiaglev, &
3776	proglevel_zdeep, thawed_humidity)
3777
3778	!! 0. Variables and parameter declaration
3779
3780	!! 0.1 Input variables
3781	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size
3782	REAL(r_std),DIMENSION (kjpindex,nslm), INTENT (in) :: shumdiag_perma !! Diagnostoc profile
3783	INTEGER(i_std), INTENT (in) :: proglevel_bottomdiaglev !! for keeping track of where the base of the diagnostic level meets the prognostic level
3784	INTEGER(i_std), INTENT (in) :: proglevel_zdeep !! for keeping track of where the prognostic levels meet zdeep
3785	REAL(r_std), DIMENSION(kjpindex), INTENT (in) :: thawed_humidity !! specified humidity of thawed soil
3786
3787	!! 0.2 Modified variables
3788
3789	!! 0.3 Output variables
3790
3791	!! 0.4 Local variables
3792	INTEGER(i_std) :: ji, jd, jv
3793
3794	IF (printlev>=3) WRITE (numout,*) 'entering update_deep_soil_misture'
3795
3796
3797	DO ji = 1, kjpindex
3798	DO jv = 1,nvm
3799	DO jd = proglevel_zdeep, ngrnd
3800	IF ( (ptn(ji,jd,jv) .GT. (ZeroCelsius+fr_dT/2.)) ) THEN
3801	shum_ngrnd_perma(ji,jd,jv) = thawed_humidity(ji)
3802	END IF
3803	END DO
3804	END DO
3805	END DO
3806
3807	DO jd = proglevel_bottomdiaglev, proglevel_zdeep-1
3808	DO ji = 1, kjpindex
3809	DO jv = 1,nvm
3810	CALL lint (diaglev(nslm), shumdiag_perma(ji,nslm), z_deepsoil,shum_ngrnd_perma(ji,proglevel_zdeep,jv), &
3811	znt(jd), shum_ngrnd_perma(ji,jd,jv), 1)
3812	END DO
3813	END DO
3814	END DO
3815
3816	IF (printlev>=3) WRITE (numout,*) ' update_deep_soil_misture done'
3817
3818	END SUBROUTINE update_deep_soil_moisture
3819
3820	!!
3821	!================================================================================================================================
3822	!! SUBROUTINE : lint
3823	!!
3824	!>\BRIEF Simple interpolation
3825	!!
3826	!! DESCRIPTION : ! Interpolation linÃÂ©aire entre des points (x1,y1) et(x2,y2))
3827	!! Ces commentaires en mauvais franÃÂ§ais permettent savoir qui a
3828	!! ecrit la subroutine :-) - DK
3829	!!
3830	!! RECENT CHANGE(S) : None
3831	!!
3832	!! MAIN OUTPUT VARIABLE(S):
3833	!!
3834	!! REFERENCE(S) : None
3835	!!
3836	!! FLOWCHART : None
3837	!! \n
3838	!_
3839	!================================================================================================================================
3840	SUBROUTINE lint(x1,y1,x2,y2,x,y,NY)
3841	!! 0. Variables and parameter declaration
3842
3843	!! 0.1 Input variables
3844
3845	REAL, INTENT(in) :: x1,x2,y1,y2,x
3846	INTEGER, INTENT(in) :: NY
3847
3848	!! 0.2 Modified variables
3849	REAL, DIMENSION(NY), INTENT(inout) :: y
3850
3851	!! 0.3 Local variables
3852	REAL, PARAMETER :: EPSILON = 1.E-10
3853
3854	IF (ABS(x1 - x2) .LT. EPSILON) THEN
3855	PRINT *, 'ERROR IN lint(x1,y1,x2,y2,y,NY) : x1==x2!'
3856	PRINT *, 'x1=',x1,' x2=',x2
3857	PRINT *, 'y1=',y1,' y2=',y2
3858	STOP
3859	END IF
3860
3861	IF (x1 .LE. x .AND. x .LE. x2) THEN
3862	y = x(y2-y1)/(x2-x1) + (y1x2 - y2*x1)/(x2-x1)
3863	! ELSE
3864	! y = UNDEF
3865	END IF
3866
3867	END SUBROUTINE lint
3868
3869	!!
3870	!================================================================================================================================
3871	!! SUBROUTINE : thermosoil_wlupdate
3872	!!
3873	!>\BRIEF Updates the long-term soil humidity
3874	!!
3875	!! DESCRIPTION :
3876	!!
3877	!! RECENT CHANGE(S) : None
3878	!!
3879	!! MAIN OUTPUT VARIABLE(S):
3880	!!
3881	!! REFERENCE(S) :
3882	!!
3883	!! FLOWCHART : None
3884	!! \n
3885	!_
3886	!================================================================================================================================
3887	SUBROUTINE thermosoil_wlupdate( kjpindex, hsd, hsdlong )
3888	!! 0. Variables and parameter declaration
3889
3890	!! 0.1 Input variables
3891	INTEGER(i_std),INTENT(in) :: kjpindex
3892	REAL(r_std),DIMENSION(kjpindex,ngrnd,nvm),INTENT(in) :: hsd
3893
3894	!! 0.2 Modified variables
3895	REAL(r_std),DIMENSION(kjpindex,ngrnd,nvm),INTENT(inout) :: hsdlong
3896
3897	!! 0.3 Local variables
3898	INTEGER(i_std) :: il
3899	REAL(r_std), PARAMETER :: tau_freezesoil = 30.*86400.
3900
3901	!
3902	! DO il = 1, ndeep
3903	! WHERE ( ( ptn(:,il,:) .GT. ZeroCelsius + fr_dT/2. ))
3904	! hsdlong(:,il,:) = ( hsd(:,il,:) * dt_sechiba + hsdlong(:,il,:) *(tau_freezesoil-dt_sechiba) ) / tau_freezesoil
3905	! ENDWHERE
3906	! END DO
3907	hsdlong(:,:,:) = ( hsd(:,:,:) * dt_sechiba + hsdlong(:,:,:) *(tau_freezesoil-dt_sechiba) ) / tau_freezesoil
3908
3909	IF (printlev>=3) WRITE (numout,*) 'entering thermosoil_wlupdate'
3910
3911	END SUBROUTINE thermosoil_wlupdate
3912
3913	!!
3914	!================================================================================================================================
3915	!! SUBROUTINE : thermosoil_getdiff_thinsnow
3916	!!
3917	!>\BRIEF Computes soil heat capacity and conductivity
3918	!!
3919	!! DESCRIPTION : Computation of the soil thermal properties; snow properties are also accounted for
3920	!!
3921	!! RECENT CHANGE(S) : None
3922	!!
3923	!! MAIN OUTPUT VARIABLE(S):
3924	!!
3925	!! REFERENCE(S) :
3926	!!
3927	!! FLOWCHART : None
3928	!! \n
3929	!_
3930	!================================================================================================================================
3931	SUBROUTINE thermosoil_getdiff_thinsnow (kjpindex, ptn, shum_ngrnd_permalong, snowdz, profil_froz)
3932
3933	!! 0. Variables and parameter declaration
3934
3935	!! 0.1 Input variables
3936	INTEGER(i_std),INTENT(in) :: kjpindex
3937	REAL(r_std),DIMENSION(kjpindex,ngrnd,nvm),INTENT(in) :: ptn
3938	REAL(r_std),DIMENSION(kjpindex,ngrnd,nvm),INTENT(in) :: shum_ngrnd_permalong
3939	REAL(r_std),DIMENSION(kjpindex,nsnow),INTENT (in) :: snowdz
3940
3941	!! 0.2 Output variables
3942	REAL(r_std),DIMENSION(kjpindex,ngrnd,nvm),INTENT(out) :: profil_froz
3943
3944	!! 0.3 Local variables
3945	REAL(r_std) :: x
3946	REAL(r_std), DIMENSION(kjpindex) :: snow_h
3947	REAL(r_std), DIMENSION(kjpindex,ngrnd) :: zx1, zx2
3948	INTEGER(i_std) :: ji,jg,jv
3949
3950	zx1 = 0
3951	zx2 = 0
3952
3953	DO ji = 1,kjpindex
3954
3955	! 1. Determine the fractions of snow and soil
3956
3957	snow_h(ji) = SUM(snowdz(ji,:))
3958
3959	IF (snow_h(ji) .LE. 0.01 .AND. snow_h(ji) .GT. 0 ) THEN
3960
3961	!
3962	! 1.1. The first level
3963	!
3964	IF ( snow_h(ji) .GT. zlt(1) ) THEN
3965
3966	! the 1st level is in the snow => the 1st layer is entirely snow
3967	zx1(ji,1) = 1.
3968	zx2(ji,1) = 0.
3969
3970	ELSE IF ( snow_h(ji) .GT. zero ) THEN
3971
3972	! the 1st level is beyond the snow and the snow is present
3973	zx1(ji,1) = snow_h(ji) / zlt(1)
3974	zx2(ji,1) = ( zlt(1) - snow_h(ji)) / zlt(1)
3975	ENDIF
3976
3977	!
3978	DO jv = 1,nvm
3979	DO jg = 1, 1
3980	!
3981	! 2. Calculate frozen profile for hydrolc.f90
3982	!
3983	IF (ptn(ji,jg,jv) .LT. ZeroCelsius-fr_dT/2.) THEN
3984	profil_froz(ji,jg,jv) = 1.
3985
3986	ELSEIF (ptn(ji,jg,jv) .GT. ZeroCelsius+fr_dT/2.) THEN
3987	profil_froz(ji,jg,jv) = 0.
3988	ELSE
3989
3990	! x is the unfrozen fraction of soil water
3991	x = (ptn(ji,jg,jv)-(ZeroCelsius-fr_dT/2.)) / fr_dT
3992	profil_froz(ji,jg,jv) = (1. - x)
3993
3994	ENDIF
3995
3996	! 3. heat capacity calculation
3997	!
3998	! 3.0 old heat capacity calculation
3999	pcapa(ji,jg,jv) = so_capa_dry + shum_ngrnd_permalong(ji,jg,jv)*(so_capa_wet - so_capa_dry)
4000
4001	! 3.1. Still some improvement from the old_version : Take into account the snow and soil fractions in the layer
4002
4003	pcapa(ji,jg,jv) = zx1(ji,jg) * sn_capa + zx2(ji,jg) * pcapa(ji,jg,jv)
4004
4005	! 3.2. Calculate the heat capacity for energy conservation check
4006	IF ( zx1(ji,jg).GT.0. ) THEN
4007	pcapa_en(ji,jg,jv) = sn_capa
4008	ELSE
4009	pcapa_en(ji,jg,jv) = pcapa(ji,jg,jv)
4010	ENDIF
4011	!
4012	!4. heat conductivity calculation
4013	!
4014	!4.0 old heat conductivity calculation
4015	pkappa(ji,jg,jv) = so_cond_dry + shum_ngrnd_permalong(ji,jg,jv)*(so_cond_wet - so_cond_dry)
4016
4017	!4.0 Still some improvement from the old_version : Take into account the snow and soil fractions in the layer
4018
4019	pkappa(ji,jg,jv) = un / ( zx1(ji,jg) / sn_cond + zx2(ji,jg) / pkappa(ji,jg,jv) )
4020
4021	END DO
4022	END DO
4023	ENDIF
4024	ENDDO
4025
4026
4027	END SUBROUTINE thermosoil_getdiff_thinsnow
4028
4029	SUBROUTINE thermosoil_rotation_update(ji, kjpindex, matrix_rot, old_veget_max)
4030	!! 0.1 Input variables
4031	INTEGER(i_std), INTENT(in) :: ji, kjpindex !! domain size
4032	REAL(r_std),DIMENSION (nvm), INTENT (in) :: old_veget_max !! max fraction of vegetation type
4033	REAL(r_std), DIMENSION (nvm, nvm), INTENT(in) :: matrix_rot !! rotation matrix
4034
4035	!! 0.4 Local variables
4036	INTEGER(i_std) :: jv, jsrc, jtar, ii
4037
4038	REAL(r_std),DIMENSION(ngrnd,nvm) :: ptn_old, dilu_ptn
4039	REAL(r_std),DIMENSION(ngrnd-1,nvm) :: cgrnd_old, dgrnd_old, dilu_cgrnd, dilu_dgrnd
4040	REAL(r_std), DIMENSION(nvm) :: maxfrac, maxfrac_new
4041
4042	!!!!----------------------------------------------------------------------------------------------
4043
4044	maxfrac = old_veget_max
4045	maxfrac_new = old_veget_max(:)
4046	DO jsrc = 1,nvm
4047	DO jtar = 1,nvm
4048	IF (matrix_rot(jsrc,jtar) .GT. 0.0) THEN
4049	maxfrac_new(jtar) = maxfrac_new(jtar) + maxfrac(jsrc) * matrix_rot(jsrc,jtar)
4050	maxfrac_new(jsrc) = maxfrac_new(jsrc) - maxfrac(jsrc) * matrix_rot(jsrc,jtar)
4051	ENDIF
4052	ENDDO
4053	ENDDO
4054
4055
4056	ptn_old = ptn(ji,:,:)
4057	cgrnd_old = cgrnd(ji,:,:)
4058	dgrnd_old = dgrnd(ji,:,:)
4059	DO jtar = 1,nvm
4060	dilu_ptn(:,:) = zero
4061	dilu_cgrnd(:,:) = zero
4062	dilu_dgrnd(:,:) = zero
4063	IF ( SUM(matrix_rot(:,jtar)) .GT. min_sechiba ) THEN
4064	DO jsrc = 1,nvm
4065	IF ( matrix_rot(jsrc,jtar) .GT. min_sechiba ) THEN
4066	dilu_ptn(:,jsrc) = ptn_old(:,jsrc)
4067	dilu_cgrnd(:,jsrc) = cgrnd_old(:,jsrc)
4068	dilu_dgrnd(:,jsrc) = dgrnd_old(:,jsrc)
4069	ENDIF
4070	ENDDO
4071	ptn(ji,:,jtar) = ptn_old(:,jtar) * maxfrac(jtar) * (1.0 - SUM(matrix_rot(jtar,:)))
4072	cgrnd(ji,:,jtar) = cgrnd_old(:,jtar) * maxfrac(jtar) * (1.0 - SUM(matrix_rot(jtar,:)))
4073	dgrnd(ji,:,jtar) = dgrnd_old(:,jtar) * maxfrac(jtar) * (1.0 - SUM(matrix_rot(jtar,:)))
4074	DO jsrc = 1,nvm
4075	ptn(ji,:,jtar) = ptn(ji,:,jtar) + maxfrac(jsrc) * matrix_rot(jsrc,jtar) * dilu_ptn(:,jsrc)
4076	cgrnd(ji,:,jtar) = cgrnd(ji,:,jtar) + maxfrac(jsrc) * matrix_rot(jsrc,jtar) * dilu_cgrnd(:,jsrc)
4077	dgrnd(ji,:,jtar) = dgrnd(ji,:,jtar) + maxfrac(jsrc) * matrix_rot(jsrc,jtar) * dilu_dgrnd(:,jsrc)
4078	ENDDO
4079	ptn(ji,:,jtar) = ptn(ji,:,jtar) / maxfrac_new(jtar)
4080	cgrnd(ji,:,jtar) = cgrnd(ji,:,jtar) / maxfrac_new(jtar)
4081	dgrnd(ji,:,jtar) = dgrnd(ji,:,jtar) / maxfrac_new(jtar)
4082	ENDIF
4083	ENDDO
4084
4085	IF (printlev>=4) THEN
4086	WRITE(numout,*) 'xuhui: debug for thermosoil rotation, ji:',ji
4087	DO ii=1,2
4088	WRITE(numout,*) 'checking first 2 layers:'
4089	WRITE(numout,*) 'ii, ptn_old(ii,:)', ii, ptn_old(ii,:)
4090	WRITE(numout,*) 'ii, ptn(ji,ii,:)', ii, ptn(ji,ii,:)
4091	WRITE(numout,*) 'ii, cgrnd_old(ii,:)', ii, cgrnd_old(ii,:)
4092	WRITE(numout,*) 'ii, cgrnd(ji,ii,:)', ii, cgrnd(ji,ii,:)
4093	WRITE(numout,*) 'ii, dgrnd_old(ii,:)', ii, dgrnd_old(ii,:)
4094	WRITE(numout,*) 'ii, dgrnd(ji,ii,:)', ii, dgrnd(ji,ii,:)
4095	ENDDO
4096	ENDIF
4097	END SUBROUTINE thermosoil_rotation_update
4098
4099	END MODULE thermosoil

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

Context Navigation

source: branches/publications/ORCHIDEE-MICT-BIOENERGY_r7298/src_sechiba/thermosoil.f90 @ 7346

Download in other formats: