Context Navigation

source: branches/publications/ORCHIDEE_CAMEO_gmd_2022/src_sechiba/thermosoil.f90 @ 7693

Last change on this file since 7693 was 6499, checked in by josefine.ghattas, 4 years ago
Removed variable cap_iw never used, see ticket #604
Property svn:keywords set to `Revision Date HeadURL Date Author Revision`
File size: 151.5 KB

Line
1	! =================================================================================================================================
2	! MODULE : thermosoil
3	!
4	! CONTACT : orchidee-help _at_ listes.ipsl.fr
5	!
6	! LICENCE : IPSL (2006)
7	! This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC
8	!
9	!>\BRIEF Calculates the soil temperatures by solving the heat
10	!! diffusion equation within the soil. This module is only used with CWRR hydrology.
11	!!
12	!!\n DESCRIPTION : General important informations about the numerical scheme and
13	!! the soil vertical discretization:\n
14	!! - the soil is zmaxt deep (by default 10m) and divided into "ngrnd" layers.
15	!! From 0-zmaxh(default 2m), the discretization is the same as for hydrology.
16	!! From zmaxh(2m) and below, the depth increase linearly (by default) or geometrically. \n
17	!! - "jg" is usually used as the index going from 1 to ngrnd to describe the
18	!! layers, from top (jg=1) to bottom (jg=ngrnd)\n
19	!! - the thermal numerical scheme is implicit finite differences.\n
20	!! -- When it is resolved in thermosoil_profile at the present timestep t, the
21	!! dependancy from the previous timestep (t-1) is hidden in the
22	!! integration coefficients cgrnd and dgrnd, which are therefore
23	!! calculated at the very end of thermosoil_main (call to
24	!! thermosoil_coef) for use in the next timestep.\n
25	!! -- At timestep t, the system becomes :\n
26	!!
27	!! T(k+1)=cgrnd(k)+dgrnd(k)*T(k) \n
28	!! -- EQ1 -- \n
29	!!
30	!! (the bottom boundary condition has been used to obtained this equation).\n
31	!! To solve it, the uppermost soil temperature T(1) is required.
32	!! It is obtained from the surface temperature Ts, which is
33	!! considered a linear extrapolation of T(1) and T(2)\n
34	!!
35	!! Ts=(1+lambda)T(1) -lambdaT(2) \n
36	!! -- EQ2--\n
37	!!
38	!! -- caveat 1 : Ts is called 'temp_soil_new' in this routine,
39	!! don' t act.\n
40	!! -- caveat 2 : actually, the surface temperature at time t Ts
41	!! depends on the soil temperature at time t through the
42	!! ground heat flux. This is again implicitly solved, with Ts(t)
43	!! expressed as :\n
44	!!
45	!! soilcap*(Ts(t)-Ts(t-1))/dt=soilflx+otherfluxes(Ts(t))\n
46	!! -- EQ3 --\n
47	!!
48	!! and the dependency from the previous timestep is hidden in
49	!! soilcap and soilflx (apparent surface heat capacity and heat
50	!! flux respectively). Soilcap and soilflx are therefore
51	!! calculated at the previous timestep, at the very end of thermosoil
52	!! (final call to thermosoil_coef) and stored to be used at the next time step.
53	!! At timestep t, EQ3 is solved for Ts in enerbil, and Ts
54	!! is used in thermosoil to get T(1) and solve EQ1.\n
55	!!
56	!! - lambda is the @tex $\mu$ @endtex of F. Hourdin' s PhD thesis, equation (A28); ie the
57	!! coefficient of the linear extrapolation of Ts (surface temperature) from T1 and T2 (ptn(jg=1) and ptn(jg=2)), so that:\n
58	!! Ts= (1+lambda)T(1)-lambdaT(2) --EQ2-- \n
59	!! lambda = (zlt(1))/((zlt(2)-zlt(1))) \n
60	!!
61	!! RECENT CHANGE(S) : - Change soil thermal properties to consider also soil texture, rev 2922.
62	!! - Change vertical discretization, rev 2917. Note: In the revised thermosoil,
63	!! cstgrnd and lskin are not needed any more. The depth znt, zlt and dlt
64	!! are computed in vertical_soil and are in meter
65	!!
66	!! REFERENCE(S) : None
67	!!
68	!! SVN :
69	!! $HeadURL$
70	!! $Date$
71	!! $Revision$
72	!! \n
73	!_ ================================================================================================================================
74
75	MODULE thermosoil
76
77	! modules used :
78	USE ioipsl
79	USE ioipsl_para
80	USE xios_orchidee
81	USE constantes
82	USE time, ONLY : one_day, dt_sechiba
83	USE constantes_soil
84	USE sechiba_io_p
85	USE grid
86	USE pft_parameters
87	USE interpol_help
88	USE vertical_soil
89
90	IMPLICIT NONE
91
92	!private and public routines :
93	PRIVATE
94	PUBLIC :: thermosoil_main, thermosoil_clear, thermosoil_initialize, thermosoil_finalize, thermosoil_xios_initialize
95
96	REAL(r_std), SAVE :: lambda !! See Module description
97	!$OMP THREADPRIVATE(lambda)
98	REAL(r_std), SAVE :: fz1, zalph !! usefull constants for diverse use
99	!$OMP THREADPRIVATE(fz1, zalph)
100	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:,:) :: ptn !! Vertically discretized soil temperature, per soil layer and per pft (K)
101	!$OMP THREADPRIVATE(ptn)
102	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: ptn_pftmean !! Vertically discretized soil temperature, mean across all pfts
103	!$OMP THREADPRIVATE(ptn_pftmean)
104
105	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: dz1 !! numerical constant used in the thermal numerical
106	!! scheme @tex ($m^{-1}$) @endtex. ; it corresponds
107	!! to the coefficient @tex $d_k$ @endtex of equation
108	!! (A.12) in F. Hourdin PhD thesis.
109	!$OMP THREADPRIVATE(dz1)
110	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: cgrnd !! integration coefficient for the numerical scheme,
111	!! see eq.1
112	!$OMP THREADPRIVATE(cgrnd)
113	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: dgrnd !! integration coefficient for the numerical scheme,
114	!! see eq.1
115	!$OMP THREADPRIVATE(dgrnd)
116
117	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: pcapa !! volumetric vertically discretized soil heat
118	!! capacity @tex ($J K^{-1} m^{-3}$) @endtex.
119	!! It depends on the soil
120	!! moisture content (shum_ngrnd_perma) and is calculated at
121	!! each time step in thermosoil_coef.
122	!$OMP THREADPRIVATE(pcapa)
123	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:,:) :: pcapa_per_pft !! volumetric vertically discretized soil heat per pft
124	!! capacity @tex ($J K^{-1} m^{-3}$) @endtex.
125	!! It depends on the soil
126	!! moisture content (shum_ngrnd_perma) and is calculated at
127	!! each time step in thermosoil_coef.
128	!$OMP THREADPRIVATE(pcapa_per_pft)
129	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: pkappa !! vertically discretized soil thermal conductivity
130	!! @tex ($W K^{-1} m^{-1}$) @endtex. Same as pcapa.
131	!$OMP THREADPRIVATE(pkappa)
132	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:,:) :: pkappa_per_pft !! vertically discretized soil thermal conductivity per pft
133	!! @tex ($W K^{-1} m^{-1}$) @endtex. Same as pcapa.
134	!$OMP THREADPRIVATE(pkappa_per_pft)
135	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: pcapa_snow !! volumetric vertically discretized snow heat
136	!! capacity @tex ($J K^{-1} m^{-3}$) @endtex.
137	!$OMP THREADPRIVATE(pcapa_snow)
138	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: pkappa_snow !! vertically discretized snow thermal conductivity
139	!! @tex ($W K^{-1} m^{-1}$) @endtex.
140	!$OMP THREADPRIVATE(pkappa_snow)
141
142	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: pcapa_en !! heat capacity used for surfheat_incr and
143	!! coldcont_incr
144	!$OMP THREADPRIVATE(pcapa_en)
145	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:,:) :: pcapa_en_per_pft !! heat capacity used for surfheat_incr and per pft
146	!! coldcont_incr
147	!$OMP THREADPRIVATE(pcapa_en_per_pft)
148	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:,:) :: ptn_beg !! ptn as it is after thermosoil_profile but before thermosoil_coef,
149	!! used in thermosoil_readjust
150	!$OMP THREADPRIVATE(ptn_beg)
151	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: temp_sol_beg !! Surface temperature at previous timestep (K)
152	!$OMP THREADPRIVATE(temp_sol_beg)
153	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: surfheat_incr !! Change in soil heat content during the timestep
154	!! @tex ($J$) @endtex.
155	!$OMP THREADPRIVATE(surfheat_incr)
156	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: coldcont_incr !! Change in snow heat content @tex ($J$) @endtex.
157	!$OMP THREADPRIVATE(coldcont_incr)
158	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: shum_ngrnd_perma !! Water saturation degree on the soil layers define by the thermic (0-1, dimensionless)
159	!$OMP THREADPRIVATE(shum_ngrnd_perma)
160
161	! Variables related to soil freezing
162	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: profil_froz !! Frozen fraction of the soil on hydrological levels (-)
163	!$OMP THREADPRIVATE(profil_froz)
164	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:,:,:):: profil_froz_pft !! Frozen fraction of the soil on hydrological levels per pft (-)
165	!$OMP THREADPRIVATE(profil_froz_pft)
166	REAL(r_std), ALLOCATABLE, SAVE,DIMENSION(:,:) :: refsoc !! Initialize soil organic carbon only used to calculate thermal insulating effect [kgC/m2]
167	!$OMP THREADPRIVATE(refsoc)
168	REAL(r_std),ALLOCATABLE, SAVE, DIMENSION (:) :: e_soil_lat !! Latent heat released or consumed in the freezing/thawing processes summed vertically
169	!! for the whole soil (J/m2) and on the whole simulation to check/correct energy conservation
170	!$OMP THREADPRIVATE(e_soil_lat)
171	REAL(r_std),ALLOCATABLE, SAVE, DIMENSION (:,:) :: pcappa_supp !! Additional surfacic heat capacity due to soil freezing for each soil layer (J/K/m2)
172	!$OMP THREADPRIVATE(pcappa_supp)
173	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: dz5 !! Used for numerical calculation [-]
174	!$OMP THREADPRIVATE(dz5)
175	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: mcs !! Saturation humidity [m3/m3]
176	!$OMP THREADPRIVATE(mcs)
177	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: QZ !! quartz content [-]
178	!$OMP THREADPRIVATE(QZ)
179	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: so_capa_dry_ns !! Dry soil Heat capacity of soils,J.m^{-3}.K^{-1}
180	!$OMP THREADPRIVATE(so_capa_dry_ns)
181	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: so_capa_ice !! Heat capacity of saturated frozen soil (J/K/m3)
182	!$OMP THREADPRIVATE(so_capa_ice)
183	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: mc_layt !! Volumetric soil moisture (liquid+ice) (m3/m3) on the thermodynamical levels at interface
184	!$OMP THREADPRIVATE(mc_layt)
185	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: mcl_layt !! Volumetric soil moisture (liquid) (m3/m3) on the thermodynamical levels at interface
186	!$OMP THREADPRIVATE(mcl_layt)
187	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: tmc_layt !! Total soil moisture content for each layer (liquid+ice) (mm) on the thermodynamical levels
188	!$OMP THREADPRIVATE(tmc_layt)
189	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:,:) :: mc_layt_pft !! Volumetric soil moisture (liquid+ice) (m3/m3) on the thermodynamical levels at interface per pft
190	!$OMP THREADPRIVATE(mc_layt_pft)
191	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:,:) :: mcl_layt_pft !! Volumetric soil moisture (liquid) (m3/m3) on the thermodynamical levels at interface per pft
192	!$OMP THREADPRIVATE(mcl_layt_pft)
193	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:,:) :: tmc_layt_pft !! Total soil moisture content for each layer (liquid+ice) (mm) on the thermodynamical levels per pft
194	!$OMP THREADPRIVATE(tmc_layt_pft)
195	INTEGER(i_std), SAVE :: brk_flag = 0 !! Flag to consider bedrock: 0.no; 1.yes
196	!$OMP THREADPRIVATE(brk_flag)
197
198	!Vertical Permafrost Carbon
199	LOGICAL, SAVE :: use_soilc_tempdiff = .FALSE. !! Do we want to activate the soil C effect on thermix
200	!$OMP THREADPRIVATE(use_soilc_tempdiff)
201	LOGICAL, SAVE :: use_refsoc = .FALSE. !! which SOC to use in thermix: a map that we read or the SOC calculated by the model.
202	!$OMP THREADPRIVATE(use_refsoc)
203	INTEGER(i_std), SAVE :: use_soilc_method !! Method to average thermal conductivity of mineral soil and organic soil:
204	!! Possible values are SOILC_METHOD_ARITHMETIC or SOILC_METHOD_GEOMETRIC
205	!$OMP THREADPRIVATE(use_soilc_method)
206	INTEGER(i_std), PARAMETER :: SOILC_METHOD_ARITHMETIC = 1 !! Index to use arithmetic mean
207	INTEGER(i_std), PARAMETER :: SOILC_METHOD_GEOMETRIC = 2 !! Index to use geometric mean
208
209	INTEGER(i_std), SAVE :: snow_cond_method !! Method to calculate snow thermal conductivity
210	!! Possible values are SNOW_COND_METHOD_DEFAULT and SNOW_COND_METHOD_DECHARME16
211	!$OMP THREADPRIVATE(snow_cond_method)
212	INTEGER(i_std), PARAMETER :: SNOW_COND_METHOD_DEFAULT = 1 !! Index for original method
213	INTEGER(i_std), PARAMETER :: SNOW_COND_METHOD_DECHARME16 = 2 !! Index to follow the method by Decharme et al 2016
214
215	CONTAINS
216
217
218	!! =============================================================================================================================
219	!! SUBROUTINE: thermosoil_xios_initialize
220	!!
221	!>\BRIEF Initialize xios dependant defintion before closing context defintion
222	!!
223	!! DESCRIPTION: Initialize xios dependant defintion before closing context defintion
224	!! Reading is deactivated if the sechiba restart file exists because the
225	!! variable should be in the restart file already.
226	!!
227	!! \n
228	!_ ==============================================================================================================================
229
230	SUBROUTINE thermosoil_xios_initialize
231
232	CHARACTER(LEN=255) :: filename, name
233
234	filename = 'reftemp.nc'
235	CALL getin_p('REFTEMP_FILE',filename)
236
237	name = filename(1:LEN_TRIM(FILENAME)-3)
238	CALL xios_orchidee_set_file_attr("reftemp_file",name=name)
239
240	! Check if the reftemp file will be read by XIOS, by IOIPSL or not at all
241	IF (xios_interpolation .AND. read_reftemp .AND. restname_in=='NONE') THEN
242	! The reftemp file will be read using XIOS
243	IF (printlev>=2) WRITE(numout,*) 'Reading of reftemp file will be done later using XIOS. The filename is ', filename
244	ELSE
245	IF (.NOT. read_reftemp) THEN
246	IF (printlev>=2) WRITE (numout,*) 'No reading of reftemp will be done because read_reftemp=FALSE'
247	ELSE IF (restname_in=='NONE') THEN
248	IF (printlev>=2) WRITE (numout,*) 'The reftemp file will be read later by IOIPSL'
249	ELSE
250	IF (printlev>=2) WRITE (numout,*) 'The reftemp file will not be read because the restart file exists.'
251	END IF
252
253	! The reftemp file will not be read by XIOS. Now deactivate albedo for XIOS.
254	CALL xios_orchidee_set_file_attr("reftemp_file",enabled=.FALSE.)
255	CALL xios_orchidee_set_field_attr("reftemp_interp",enabled=.FALSE.)
256	ENDIF
257
258	END SUBROUTINE thermosoil_xios_initialize
259
260	!! =============================================================================================================================
261	!! SUBROUTINE : thermosoil_initialize
262	!!
263	!>\BRIEF Allocate module variables, read from restart file or initialize with default values
264	!!
265	!! DESCRIPTION : Allocate module variables, read from restart file or initialize with default values.
266	!! Call thermosoil_var_init to calculate physical constants.
267	!! Call thermosoil_coef to calculate thermal soil properties.
268	!!
269	!! RECENT CHANGE(S) : None
270	!!
271	!! REFERENCE(S) : None
272	!!
273	!! FLOWCHART : None
274	!! \n
275	!_ ==============================================================================================================================
276	SUBROUTINE thermosoil_initialize (kjit, kjpindex, rest_id, &
277	temp_sol_new, snow, shumdiag_perma, &
278	soilcap, soilflx, depth_organic_soil, &
279	stempdiag, gtemp, mc_layh, &
280	mcl_layh, tmc_layh, njsc, &
281	frac_snow_veg,frac_snow_nobio,totfrac_nobio, &
282	snowdz, snowrho, snowtemp, lambda_snow, cgrnd_snow, &
283	dgrnd_snow, pb, &
284	som_total, veget_max, mc_layh_pft, mcl_layh_pft, tmc_layh_pft)
285
286	!! 0. Variable and parameter declaration
287	!! 0.1 Input variables
288	INTEGER(i_std), INTENT(in) :: kjit !! Time step number (unitless)
289	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (unitless)
290	INTEGER(i_std),INTENT (in) :: rest_id !! Restart file identifier (unitless)
291	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_sol_new !! Surface temperature at the present time-step,
292	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: snow !! Snow mass (kg)
293	REAL(r_std),DIMENSION (kjpindex,nslm), INTENT (in) :: shumdiag_perma !! Soil saturation degree (0-1, unitless)
294	REAL(r_std),DIMENSION (kjpindex,nslm), INTENT (in) :: mc_layh !! Volumetric soil moisture content (liquid+ice) for hydrological layers, at node (m3/m3)
295	REAL(r_std),DIMENSION (kjpindex,nslm), INTENT (in) :: mcl_layh !! Volumetric soil moisture content (liquid) for hydrological layers, at node (m3/m3)
296	REAL(r_std),DIMENSION (kjpindex,nslm), INTENT (in) :: tmc_layh !! Total soil moisture content(liquid+ice) for hydrological layers (mm)
297	INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: njsc !! Index of the dominant soil textural class in the grid cell (1-nscm, unitless)
298	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: frac_snow_veg !! Snow cover fraction on vegeted area
299	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT(in) :: frac_snow_nobio !! Snow cover fraction on non-vegeted area
300	REAL(r_std),DIMENSION (kjpindex),INTENT(in) :: totfrac_nobio !! Total fraction of continental ice+lakes+cities+...
301	!! (unitless,0-1)
302	REAL(r_std),DIMENSION (kjpindex,nsnow),INTENT (in) :: snowdz !! Snow depth [m]
303	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(in) :: snowrho !! Snow density (Kg/m^3)
304	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(in) :: snowtemp !! Snow temperature (K)
305	REAL(r_std), DIMENSION (kjpindex), INTENT (in) :: pb !! Surface presure (hPa)
306	REAL(r_std), DIMENSION(kjpindex,ngrnd,nvm,nelements), INTENT (in) :: som_total !! total soil organic matter for use in thermal calcs (g/m**3)
307	REAL(r_std), DIMENSION (kjpindex,nvm), INTENT(in) :: veget_max !! Fraction of vegetation type
308	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]
309	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]
310	REAL(r_std),DIMENSION (kjpindex,nslm,nvm), INTENT (in):: tmc_layh_pft !! Total soil moisture content for each layer in hydrol(liquid+ice) [mm]
311
312	!! 0.2 Output variables
313	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: soilcap !! apparent surface heat capacity considering snow and soil surface (J m-2 K-1)
314	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: soilflx !! apparent soil heat flux considering snow and soil surface (W m-2)
315	REAL(r_std),DIMENSION (kjpindex,nslm), INTENT (out) :: stempdiag !! temperature profile on the levels in hydrol(K)
316	REAL(r_std),DIMENSION (kjpindex),INTENT(out) :: gtemp !! First soil layer temperature
317
318	!! 0.3 Modified variables
319	REAL(r_std), DIMENSION( kjpindex), INTENT (inout) :: depth_organic_soil!! Depth at which there is still organic matter (m)
320	REAL(r_std), DIMENSION (kjpindex), INTENT(inout) :: lambda_snow !! Coefficient of the linear extrapolation of surface temperature
321	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT (inout):: cgrnd_snow !! Integration coefficient for snow numerical scheme
322	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT (inout):: dgrnd_snow !! Integration coefficient for snow numerical scheme
323
324	!! 0.4 Local variables
325	INTEGER(i_std) :: ier, i, jg, iv, m, jv, jsc
326	LOGICAL :: calculate_coef !! Local flag to initialize variables by call to thermosoil_coef
327	CHARACTER(LEN=10) :: part_str !! string suffix indicating an index
328	CHARACTER(LEN=80) :: var_name
329	!_ ================================================================================================================================
330
331
332	!
333	! !! Flag to consider bedrock at deeper layers
334	! !! It affects heat capacity and thermal conductivity (energy balance).
335	!
336	!Config Key = BEDROCK_FLAG
337	!Config Desc = Flag to consider bedrock at deeper layers.
338	!Config If =
339	!Config Def = 0
340	!Config Help = 0, no, 1, yes.
341	!Config Units = [FLAG]
342	brk_flag = 0
343	CALL getin_p('BEDROCK_FLAG', brk_flag)
344
345	IF (ok_freeze_thermix .AND. ok_soil_carbon_discretization) THEN
346	use_soilc_tempdiff = .false.
347	!Config Key = USE_SOILC_TEMPDIFF
348	!Config Desc = insolation effect of the organic top soil layer
349	!Config If = OK_SOIL_CARBON_DISCRETIZATION
350	!Config Def = FALSE
351	!Config Help =
352	!Config Units = [FLAG]
353	CALL getin_p('USE_SOILC_TEMPDIFF', use_soilc_tempdiff)
354
355	IF (use_soilc_tempdiff) THEN
356	USE_REFSOC = .TRUE.
357	!Config Key = USE_REFSOC
358	!Config Desc = Read a SOC map to perform the insolation effect
359	!Config If = USE_SOILC_TEMPDIFF
360	!Config Def = TRUE
361	!Config Help =
362	!Config Units = [FLAG]
363	CALL getin_p('USE_REFSOC',use_refsoc)
364	ENDIF
365	ENDIF
366	!Config Key = USE_SOILC_METHOD
367	!Config Desc = Flag to control the way to average thermal conductivity of mineral soil and organic soil
368	!Config If = OK_SOIL_CARBON_DISCRETIZATION
369	!Config Def = 1
370	!Config Help = 1=arithmetic mean ; 2=geometric mean
371	!Config Units = [FLAG]
372	use_soilc_method = SOILC_METHOD_ARITHMETIC
373	CALL getin_p('USE_SOILC_METHOD', use_soilc_method)
374	IF ( (use_soilc_method /= SOILC_METHOD_ARITHMETIC) .AND. (use_soilc_method /= SOILC_METHOD_GEOMETRIC) ) THEN
375	CALL ipslerr_p(3,'thermosoil_initialize', 'Error in variable use_soilc_method','USE_SOILC_METHOD must be equal 1 or 2','')
376	END IF
377
378
379	!Config Key = SNOW_COND_METHOD
380	!Config Desc = Flag to choose the way to calculate snow thermal conductivity
381	!Config If = OK_SOIL_CARBON_DISCRETIZATION
382	!Config Def = 1
383	!Config Help = 1: original 2: follows Decharme et al 2016
384	!Config Units = [1=original method, 2=method by Decharme et al 2016]
385	snow_cond_method = SNOW_COND_METHOD_DEFAULT
386	CALL getin_p('SNOW_COND_METHOD', snow_cond_method)
387
388	IF (printlev >= 3) WRITE (numout,*) 'Start thermosoil_initialize '
389
390	!! 1. Allocate soil temperatures variables
391	ALLOCATE (refsoc(kjpindex,ngrnd),stat=ier)
392	IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of refsoc','','')
393
394	ALLOCATE (ptn(kjpindex,ngrnd,nvm),stat=ier)
395	IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of ptn','','')
396
397	ALLOCATE (ptn_pftmean(kjpindex,ngrnd),stat=ier)
398	IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of ptn_pftmean','','')
399
400	ALLOCATE (dz1(ngrnd),stat=ier)
401	IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of dz1','','')
402
403	ALLOCATE (cgrnd(kjpindex,ngrnd-1),stat=ier)
404	IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of cgrnd','','')
405
406	ALLOCATE (dgrnd(kjpindex,ngrnd-1),stat=ier)
407	IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of dgrnd','','')
408
409	ALLOCATE (pcapa(kjpindex,ngrnd),stat=ier)
410	IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of pcapa','','')
411
412	ALLOCATE (pkappa(kjpindex,ngrnd),stat=ier)
413	IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of pkappa','','')
414
415	ALLOCATE (pkappa_per_pft(kjpindex,ngrnd,nvm),stat=ier)
416	IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of pkappa_per_pft','','')
417
418	ALLOCATE (pcapa_per_pft(kjpindex,ngrnd,nvm),stat=ier)
419	IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of pcapa_per_pft','','')
420
421	ALLOCATE (pcapa_en_per_pft(kjpindex,ngrnd,nvm),stat=ier)
422	IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of pcapa_en_per_pft','','')
423
424	ALLOCATE (pcapa_snow(kjpindex,nsnow),stat=ier)
425	IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of pcapa_snow','','')
426
427	ALLOCATE (pkappa_snow(kjpindex,nsnow),stat=ier)
428	IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of pkappa_snow','','')
429
430	! Temporary fix: Initialize following variable because they are output to xios before the first calculation
431	pcapa = 0
432	pkappa = 0
433	pkappa_per_pft = 0
434	pcapa_per_pft = 0
435	pcapa_en_per_pft = 0
436	pcapa_snow = 0
437	pkappa_snow = 0
438
439	ALLOCATE (surfheat_incr(kjpindex),stat=ier)
440	IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of surfheat_incr','','')
441
442	ALLOCATE (coldcont_incr(kjpindex),stat=ier)
443	IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of coldcont_incr','','')
444
445	ALLOCATE (pcapa_en(kjpindex,ngrnd),stat=ier)
446	IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of pcapa_en','','')
447	! Initialization to zero used at first time step in thermosoil_energy_diag, only for diagnostic variables coldcont_incr and surfheat_incr
448	pcapa_en(:,:) = 0.
449
450	ALLOCATE (ptn_beg(kjpindex,ngrnd,nvm),stat=ier)
451	IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of ptn_beg','','')
452
453	ALLOCATE (temp_sol_beg(kjpindex),stat=ier)
454	IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of temp_sol_beg','','')
455
456	ALLOCATE (shum_ngrnd_perma(kjpindex,ngrnd),stat=ier)
457	IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of shum_ngrnd_perma','','')
458	shum_ngrnd_perma(:,:) = 0
459
460	ALLOCATE (profil_froz(kjpindex,ngrnd),stat=ier)
461	IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of profil_froz','','')
462
463	ALLOCATE (profil_froz_pft(kjpindex,ngrnd,nvm),stat=ier)
464	IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of profil_froz_pft','','')
465
466	IF (ok_freeze_thermix) THEN
467	ALLOCATE (pcappa_supp(kjpindex,ngrnd),stat=ier)
468	IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of ok_freeze_termix','','')
469	! Initialization to zero used at first time step only for diagnostic output.
470	! This variable is only used in thermosoil_readajust and always calculated before in thermosoil_getdiff.
471	pcappa_supp(:,:) = 0.
472	END IF
473	IF (ok_Ecorr) THEN
474	ALLOCATE (e_soil_lat(kjpindex),stat=ier)
475	IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of e_soil_lat','','')
476	END IF
477
478	ALLOCATE (dz5(ngrnd),stat=ier)
479	IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of dz5','','')
480
481	ALLOCATE (mc_layt(kjpindex,ngrnd),stat=ier)
482	IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of mc_layt','','')
483
484	ALLOCATE (mcl_layt(kjpindex,ngrnd),stat=ier)
485	IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of mcl_layt','','')
486
487	ALLOCATE (tmc_layt(kjpindex,ngrnd),stat=ier)
488	IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of tmc_layt','','')
489
490	ALLOCATE (mc_layt_pft(kjpindex,ngrnd,nvm),stat=ier)
491	IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of mc_layt_pft','','')
492
493	ALLOCATE (mcl_layt_pft(kjpindex,ngrnd,nvm),stat=ier)
494	IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of mcl_layt_pft','','')
495
496	ALLOCATE (tmc_layt_pft(kjpindex,ngrnd,nvm),stat=ier)
497	IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of tmc_layt_pft','','')
498
499	ALLOCATE (mcs(nscm),stat=ier)
500	IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of mcs','','')
501
502	ALLOCATE (QZ(nscm),stat=ier)
503	IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of QZ','','')
504
505	ALLOCATE (so_capa_dry_ns(nscm),stat=ier)
506	IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of so_capa_dry_ns','','')
507
508	ALLOCATE (so_capa_ice(nscm),stat=ier)
509	IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of so_capa_ice','','')
510
511
512	!! Soil texture choose
513	SELECTCASE (nscm)
514	CASE (3)
515	QZ(:) = QZ_fao(:)
516	so_capa_dry_ns(:) = so_capa_dry_ns_fao(:)
517	mcs(:) = mcs_fao(:)
518	CASE (12)
519	QZ(:) = QZ_usda(:)
520	so_capa_dry_ns(:) = so_capa_dry_ns_usda(:)
521	mcs(:) = mcs_usda(:)
522	CASE DEFAULT
523	WRITE (numout,*) 'Unsupported soil type classification. Choose between zobler, fao and usda according to the map'
524	STOP 'thermosoil_initialize'
525	ENDSELECT
526
527
528	!! 2. Initialize variable from restart file or with default values
529
530	!! Reads restart files for soil temperatures only. If no restart file is
531	!! found, the initial soil temperature is by default set to 280K at all depths. The user
532	!! can decide to initialize soil temperatures at an other value, in which case he should set the flag THERMOSOIL_TPRO
533	!! to this specific value in the run.def.
534	IF (printlev>=3) WRITE (numout,*) 'Read restart file for THERMOSOIL variables'
535
536	CALL ioconf_setatt_p('UNITS', 'K')
537	CALL ioconf_setatt_p('LONG_NAME','Soil Temperature profile')
538	CALL restget_p (rest_id, 'ptn', nbp_glo, ngrnd, nvm, kjit, .TRUE., ptn, "gather", nbp_glo, index_g)
539
540	! Initialize ptn if it was not found in restart file
541	IF (ALL(ptn(:,:,:)==val_exp)) THEN
542	! ptn was not found in restart file
543
544	IF (read_reftemp) THEN
545	! Read variable ptn from file
546	CALL thermosoil_read_reftempfile(kjpindex,lalo,ptn(:,:,1))
547	DO jv = 2,nvm
548	ptn(:,:,jv)=ptn(:,:,1)
549	ENDDO ! jv = 1,nvm
550	ELSE
551	! Initialize ptn with a constant value which can be set in run.def
552
553	!Config Key = THERMOSOIL_TPRO
554	!Config Desc = Initial soil temperature profile if not found in restart
555	!Config Def = 280.
556	!Config If = OK_SECHIBA
557	!Config Help = The initial value of the temperature profile in the soil if
558	!Config its value is not found in the restart file. Here
559	!Config we only require one value as we will assume a constant
560	!Config throughout the column.
561	!Config Units = Kelvin [K]
562	CALL setvar_p (ptn, val_exp,'THERMOSOIL_TPRO',280._r_std)
563	END IF
564	END IF
565
566	CALL restget_p (rest_id, 'refsoc', nbp_glo, ngrnd, 1, kjit, .TRUE., refsoc, "gather", nbp_glo, index_g)
567	IF (ALL(refsoc(:,:) == val_exp)) THEN
568	IF (use_refsoc) THEN
569	CALL thermosoil_read_refsoc_file(kjpindex,lalo,neighbours, resolution, contfrac)
570	ENDIF
571	ENDIF
572
573	! Initialize ptn_beg (variable needed in thermosoil_readadjust called from thermosoil_coef)
574	ptn_beg(:,:,:) = ptn(:,:,:)
575
576	! Initialize temp_sol_beg with values from previous time-step
577	temp_sol_beg(:) = temp_sol_new(:)
578
579	! Read e_soil_lat from restart file or initialize
580	IF (ok_Ecorr) THEN
581	CALL restget_p (rest_id, 'e_soil_lat', nbp_glo, 1, 1, kjit, .TRUE., &
582	e_soil_lat, "gather", nbp_glo, index_g)
583	CALL setvar_p (e_soil_lat, val_exp,'NO_KEYWORD',zero)
584	END IF
585
586	! Read gtemp from restart file
587	CALL restget_p (rest_id, 'gtemp', nbp_glo, 1, 1, kjit, .TRUE., &
588	gtemp, "gather", nbp_glo, index_g)
589	CALL setvar_p (gtemp, val_exp,'NO_KEYWORD',zero)
590
591
592	! Read variables calculated in thermosoil_coef from restart file
593	! If the variables were not found in the restart file, the logical
594	! calculate_coef will be true and thermosoil_coef will be called further below.
595	! These variables need to be in the restart file to avoid a time shift that
596	! would be done using thermosoil_coef at this stage.
597	calculate_coef=.FALSE.
598	CALL ioconf_setatt_p('UNITS', 'J m-2 K-1')
599	CALL ioconf_setatt_p('LONG_NAME','Apparent surface heat capacity')
600	CALL restget_p (rest_id, 'soilcap', nbp_glo, 1, 1, kjit, .TRUE., &
601	soilcap, "gather", nbp_glo, index_g)
602	IF (ALL(soilcap(:)==val_exp)) calculate_coef=.TRUE.
603
604	CALL ioconf_setatt_p('UNITS', 'W m-2')
605	CALL ioconf_setatt_p('LONG_NAME','Apparent soil heat flux')
606	CALL restget_p (rest_id, 'soilflx', nbp_glo, 1, 1, kjit, .TRUE., &
607	soilflx, "gather", nbp_glo, index_g)
608	IF (ALL(soilflx(:)==val_exp)) calculate_coef=.TRUE.
609
610	CALL ioconf_setatt_p('UNITS', '')
611	CALL ioconf_setatt_p('LONG_NAME','Integration coefficient for the numerical scheme')
612	CALL restget_p (rest_id, 'cgrnd', nbp_glo, ngrnd-1, 1, kjit, .TRUE., &
613	cgrnd, "gather", nbp_glo, index_g)
614	IF (ALL(cgrnd(:,:)==val_exp)) calculate_coef=.TRUE.
615
616	CALL ioconf_setatt_p('UNITS', '')
617	CALL ioconf_setatt_p('LONG_NAME','Integration coefficient for the numerical scheme')
618	CALL restget_p (rest_id, 'dgrnd', nbp_glo, ngrnd-1, 1, kjit, .TRUE., &
619	dgrnd, "gather", nbp_glo, index_g)
620	IF (ALL(dgrnd(:,:)==val_exp)) calculate_coef=.TRUE.
621
622	CALL ioconf_setatt_p('UNITS', '')
623	CALL ioconf_setatt_p('LONG_NAME','Integration coefficient for the numerical scheme')
624	CALL restget_p (rest_id, 'cgrnd_snow', nbp_glo, nsnow, 1, kjit, .TRUE., &
625	cgrnd_snow, "gather", nbp_glo, index_g)
626	IF (ALL(cgrnd_snow(:,:)==val_exp)) calculate_coef=.TRUE.
627
628	CALL ioconf_setatt_p('UNITS', '')
629	CALL ioconf_setatt_p('LONG_NAME','Integration coefficient for the numerical scheme')
630	CALL restget_p (rest_id, 'dgrnd_snow', nbp_glo, nsnow, 1, kjit, .TRUE., &
631	dgrnd_snow, "gather", nbp_glo, index_g)
632	IF (ALL(dgrnd_snow(:,:)==val_exp)) calculate_coef=.TRUE.
633
634	CALL ioconf_setatt_p('UNITS', '')
635	CALL ioconf_setatt_p('LONG_NAME','Coefficient of the linear extrapolation of surface temperature')
636	CALL restget_p (rest_id, 'lambda_snow', nbp_glo, 1, 1, kjit, .TRUE., &
637	lambda_snow, "gather", nbp_glo, index_g)
638	IF (ALL(lambda_snow(:)==val_exp)) calculate_coef=.TRUE.
639
640	!! 2.2 Computes some physical constants and arrays depending on the soil vertical discretization
641
642	! Calculate so_capa_ice
643	DO jsc = 1, nscm
644	so_capa_ice(jsc) = so_capa_dry + mcs(jsc)capa_icerho_ice
645	END DO
646	IF (printlev>=2) WRITE(numout,*) 'Calculation of so_capa_ice(:)=', so_capa_ice(:),' and capa_ice=',capa_ice
647
648	! Computing some usefull constants for the numerical scheme
649	! Use znt(depth of nodes) and zlt(depth of deeper layer interface) from vertical_soil module.
650	DO jg=1,ngrnd-1
651	dz1(jg) = un / (znt(jg+1) - znt(jg))
652	dz5(jg) = (zlt(jg) - znt(jg)) * dz1(jg)
653	ENDDO
654	dz1(ngrnd) = 0.0
655	dz5(ngrnd) = 0.0
656	lambda = znt(1) * dz1(1)
657
658
659	CALL ioconf_setatt_p('UNITS', 'K')
660	CALL ioconf_setatt_p('LONG_NAME','Soil Temperature profile mean over pft')
661	CALL restget_p (rest_id, 'ptn_pftmean', nbp_glo, ngrnd, 1, kjit, .TRUE., ptn_pftmean, "gather", nbp_glo, index_g)
662	IF (ALL(ptn_pftmean(:,:)==val_exp)) THEN
663	! Initialize ptn_pftmean if not found in restart file
664	IF (ok_soil_carbon_discretization) THEN
665	ptn_pftmean(:,:)=0.0
666	DO m=1,nvm
667	DO jg = 1, ngrnd
668	ptn_pftmean(:,jg) = ptn_pftmean(:,jg) + ptn(:,jg,m) * veget_max(:,m)
669	ENDDO
670	ENDDO
671	ELSE
672	! For this case, ptn is constant over all pfts. Use here PFT=1 to initialize ptn_pftmean.
673	ptn_pftmean(:,:) = ptn(:,:,1)
674	END IF
675	END IF
676
677	! Send out the temperature profile on the first nslm levels(the levels treated in hydrol)
678	stempdiag(:,:) = ptn_pftmean(:,1:nslm)
679
680
681	!! 2.3. Computes cgrnd, dgrnd, soilflx and soilcap coefficients only if they were not found in restart file.
682	IF (calculate_coef) THEN
683	! Interpolate variables needed by thermosoil_coef to the thermal levels
684	CALL thermosoil_humlev(kjpindex, shumdiag_perma, mc_layh, mcl_layh, tmc_layh, &
685	mc_layh_pft, mcl_layh_pft, tmc_layh_pft)
686
687	IF (printlev>=3) WRITE (numout,*) 'thermosoil_coef will be called in the intialization phase'
688	CALL thermosoil_coef (&
689	kjpindex, temp_sol_new, snow, njsc, &
690	frac_snow_veg, frac_snow_nobio, totfrac_nobio, &
691	snowdz, snowrho, snowtemp, pb, &
692	veget_max, som_total, depth_organic_soil, &
693	ptn, ptn_pftmean, &
694	soilcap, soilflx, cgrnd, dgrnd,&
695	lambda_snow, cgrnd_snow, dgrnd_snow)
696
697	END IF
698
699	END SUBROUTINE thermosoil_initialize
700
701
702	!! ================================================================================================================================
703	!! SUBROUTINE : thermosoil_main
704	!!
705	!>\BRIEF Thermosoil_main computes the soil thermal properties and dynamics, ie solves
706	!! the heat diffusion equation within the soil.
707	!!
708	!! DESCRIPTION : The resolution of the soil heat diffusion equation
709	!! relies on a numerical finite-difference implicit scheme
710	!! fully described in the reference and in the header of the thermosoil module.
711	!! - The dependency of the previous timestep hidden in the
712	!! integration coefficients cgrnd and dgrnd (EQ1), calculated in thermosoil_coef, and
713	!! called at the end of the routine to prepare for the next timestep.
714	!! - The effective computation of the new soil temperatures is performed in thermosoil_profile.
715	!!
716	!! - thermosoil_coef calculates the coefficients for the numerical scheme for the very first iteration of thermosoil;
717	!! after that, thermosoil_coef is called only at the end of the module to calculate the coefficients for the next timestep.
718	!! - thermosoil_profile solves the numerical scheme.\n
719	!!
720	!! - Flags : one unique flag : THERMOSOIL_TPRO (to be set to the desired initial soil in-depth temperature in K; by default 280K)
721	!!
722	!! RECENT CHANGE(S) : Change vertical discretization (consistent with hydrology layers) and soil thermal properties (taking into account soil texture effects).
723	!!
724	!! MAIN OUTPUT VARIABLE(S): vertically discretized soil temperatures ptn, soil
725	!! thermal properties (pcapa, pkappa), apparent surface heat capacity (soilcap)
726	!! and heat flux (soilflx) to be used in enerbil at the next timestep to solve
727	!! the surface energy balance.
728	!!
729	!! REFERENCE(S) :
730	!! - Hourdin, F. (1992). Study and numerical simulation of the general circulation of planetary atmospheres,
731	!! Ph.D. thesis, Paris VII University. Remark: the part of F. Hourdin' s PhD thesis relative to the thermal
732	!! integration scheme has been scanned and is provided along with the documentation, with name :
733	!! Hourdin_1992_PhD_thermal_scheme.pdf
734	!!
735	!! FLOWCHART :
736	!! \latexonly
737	!! \includegraphics[scale = 1]{thermosoil_flowchart.png}
738	!! \endlatexonly
739	!!
740	!! \ A flag to activate the heat production by soil microbial activityn
741	!_ ================================================================================================================================
742
743	SUBROUTINE thermosoil_main (kjit, kjpindex, &
744	index, indexgrnd, &
745	temp_sol_new, snow, soilcap, soilflx, &
746	shumdiag_perma, stempdiag, ptnlev1, rest_id, hist_id, hist2_id, &
747	snowdz,snowrho,snowtemp,gtemp,pb,&
748	mc_layh, mcl_layh, tmc_layh, mc_layh_pft, mcl_layh_pft, tmc_layh_pft, njsc, &
749	depth_organic_soil, heat_Zimov, deeptemp_prof, deephum_prof,&
750	som_total, veget_max, frac_snow_veg,frac_snow_nobio,totfrac_nobio,temp_sol_add, &
751	lambda_snow, cgrnd_snow, dgrnd_snow)
752
753	!! 0. Variable and parameter declaration
754
755	!! 0.1 Input variables
756
757	INTEGER(i_std), INTENT(in) :: kjit !! Time step number (unitless)
758	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (unitless)
759	INTEGER(i_std),INTENT (in) :: rest_id,hist_id !! Restart_ file and history file identifier
760	!! (unitless)
761	INTEGER(i_std),INTENT (in) :: hist2_id !! history file 2 identifier (unitless)
762	INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: index !! Indeces of the points on the map (unitless)
763	INTEGER(i_std),DIMENSION (kjpindex*ngrnd), INTENT (in):: indexgrnd !! Indeces of the points on the 3D map (vertical
764	!! dimension towards the ground) (unitless)
765	REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: temp_sol_new !! Surface temperature at the present time-step,
766	!! Ts @tex ($K$) @endtex
767	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: snow !! Snow mass @tex ($kg$) @endtex.
768	!! Caveat: when there is snow on the
769	!! ground, the snow is integrated into the soil for
770	!! the calculation of the thermal dynamics. It means
771	!! that the uppermost soil layers can completely or
772	!! partially consist in snow. In the second case, zx1
773	!! and zx2 are the fraction of the soil layer
774	!! consisting in snow and 'normal' soil, respectively
775	!! This is calculated in thermosoil_coef.
776	REAL(r_std),DIMENSION (kjpindex,nslm), INTENT (in) :: shumdiag_perma !! Soil saturation degree (0-1, unitless)
777	REAL(r_std), DIMENSION(kjpindex),INTENT (in) :: depth_organic_soil !! Depth at which there is still organic matter (m)
778	REAL(r_std), DIMENSION(kjpindex,ngrnd,nvm), INTENT (in) :: heat_Zimov !! Heating associated with decomposition [W/m**3 soil]
779	REAL(r_std), DIMENSION(kjpindex,ngrnd,nvm,nelements), INTENT (in) :: som_total !! Total soil organic matter for use in thermal calcs (g/m**3)
780	REAL(r_std), DIMENSION (kjpindex,nvm), INTENT (in) :: veget_max !! Fraction of vegetation type
781	REAL(r_std),DIMENSION (kjpindex,nsnow),INTENT (in) :: snowdz !! Snow depth [m]
782	REAL(r_std),DIMENSION (kjpindex,nsnow),INTENT (in) :: snowrho !! Snow density (Kg/m^3)
783	REAL(r_std),DIMENSION (kjpindex,nsnow),INTENT (inout) :: snowtemp !! Snow temperature (K)
784	REAL(r_std), DIMENSION (kjpindex),INTENT (in) :: pb !! Surface presure (hPa)
785	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)
786	REAL(r_std),DIMENSION (kjpindex,nslm), INTENT (in) :: mcl_layh !! Volumetric soil moisture content for each layer in hydrol at nodes(liquid) (m3/m3)
787	REAL(r_std),DIMENSION (kjpindex,nslm), INTENT (in) :: tmc_layh !! Total soil moisture content for each layer in hydrol(liquid + ice) (mm)
788	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)
789	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)
790	REAL(r_std),DIMENSION (kjpindex,nslm,nvm), INTENT (in):: tmc_layh_pft !! Total soil moisture content for each layer in hydrol(liquid + ice) (mm)
791	INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: njsc !! Index of the dominant soil textural class in the grid cell (1-nscm, unitless)
792	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: frac_snow_veg !! Snow cover fraction on vegeted area
793	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT(in) :: frac_snow_nobio !! Snow cover fraction on non-vegeted area
794	REAL(r_std),DIMENSION (kjpindex),INTENT(in) :: totfrac_nobio !! Total fraction of continental ice+lakes+cities+...
795	!!(unitless,0-1)
796	REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: temp_sol_add !! additional surface temperature due to the melt of first layer
797	!! at the present time-step @tex ($K$) @endtex
798
799	!! 0.2 Output variables
800
801	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: ptnlev1 !! 1st level soil temperature
802	REAL(r_std),DIMENSION (kjpindex),INTENT(out) :: gtemp !! First soil layer temperature
803	REAL(r_std), DIMENSION (kjpindex,ngrnd,nvm), INTENT (out) :: deephum_prof !! moisture on a deep thermodynamic profile for permafrost calcs
804	REAL(r_std), DIMENSION (kjpindex,ngrnd,nvm), INTENT (out) :: deeptemp_prof!! temp on a deep thermodynamic profile for permafrost calcs
805
806	!! 0.3 Modified variables
807
808	REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: soilcap !! apparent surface heat capacity considering snow and soil surface
809	!! @tex ($J m^{-2} K^{-1}$) @endtex
810	REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: soilflx !! apparent soil heat flux considering snow and soil surface
811	!! @tex ($W m^{-2}$) @endtex
812	!! , positive
813	!! towards the soil, writen as Qg (ground heat flux)
814	!! in the history files, and computed at the end of
815	!! thermosoil for the calculation of Ts in enerbil,
816	!! see EQ3.
817	REAL(r_std),DIMENSION (kjpindex,nslm), INTENT (out) :: stempdiag !! temperature profile @tex ($K$) @endtex
818	REAL(r_std),DIMENSION (kjpindex), INTENT(inout) :: lambda_snow !! Coefficient of the linear extrapolation of surface temperature
819	REAL(r_std),DIMENSION (kjpindex,nsnow), INTENT (inout):: cgrnd_snow !! Integration coefficient for snow numerical scheme
820	REAL(r_std),DIMENSION (kjpindex,nsnow), INTENT (inout):: dgrnd_snow !! Integration coefficient for snow numerical scheme
821
822	!! 0.4 Local variables
823
824	INTEGER(i_std) :: jv,ji,ii, jg, m
825	REAL(r_std),DIMENSION (kjpindex) :: snowtemp_weighted!! Snow temperature weighted by snow density, only for diag (K)
826	REAL(r_std),DIMENSION (kjpindex, nsnow) :: pkappa_snow_diag !! Only for diag, containing xios_default_val
827	REAL(r_std),DIMENSION (kjpindex, nsnow) :: pcapa_snow_diag !! Only for diag, containing xios_default_val
828	REAL(r_std),DIMENSION (kjpindex, nsnow) :: snowtemp_diag !! Only for diag, containing xios_default_val
829	LOGICAL :: ok_zimov !! A flag to activate the heat production by soil microbial activity
830	!_ ================================================================================================================================
831
832	!! 3. Put the soil wetness diagnostic on the levels of the soil temperature
833
834	!!?? this could logically be put just before the last call to
835	!!thermosoil_coef, as the results are used there...
836	CALL thermosoil_humlev(kjpindex, shumdiag_perma, mc_layh, mcl_layh, tmc_layh, &
837	mc_layh_pft, mcl_layh_pft, tmc_layh_pft)
838
839	!! 4. Effective computation of the soil temperatures profile.
840	!! cgrnd and dgrnd have been calculated in thermosoil_coef at the previous time step
841	!! but they are correct for the actual time-step.
842	CALL thermosoil_profile (kjpindex, temp_sol_new, &
843	frac_snow_veg, frac_snow_nobio, totfrac_nobio, &
844	ptn, ptn_pftmean, stempdiag, snowtemp, &
845	cgrnd_snow, dgrnd_snow)
846
847
848	!! 5. Call to thermosoil_energy_diag for calculation of diagnostic variables
849	CALL thermosoil_energy_diag(kjpindex, temp_sol_new, soilcap)
850
851	!! Save ptn at current stage, to be used in thermosoil_readjust
852	ptn_beg(:,:,:) = ptn(:,:,:)
853
854	!! 6. Writing the history files according to the ALMA standards (or not..)
855
856	! Add XIOS default value where no snow
857	DO ji=1,kjpindex
858	IF (snow(ji) .GT. zero) THEN
859	pkappa_snow_diag(ji,:) = pkappa_snow(ji,:)
860	pcapa_snow_diag(ji,:) = pcapa_snow(ji,:)
861	snowtemp_diag(ji,:) = snowtemp(ji,:)
862	ELSE
863	pkappa_snow_diag(ji,:) = xios_default_val
864	pcapa_snow_diag(ji,:) = xios_default_val
865	snowtemp_diag(ji,:) = xios_default_val
866	END IF
867	END DO
868
869	DO ji=1,kjpindex
870	! Use min_sechiba instead of zero to avoid problem with division by zero
871	IF (snow(ji) .GT. min_sechiba) THEN
872	snowtemp_weighted(ji) = SUM(snowtemp(ji,:)*snowrho(ji,:))/SUM(snowrho(ji,:))
873	ELSE
874	snowtemp_weighted(ji) = xios_default_val
875	END IF
876	END DO
877	CALL xios_orchidee_send_field("snowtemp_weighted",snowtemp_weighted)
878	CALL xios_orchidee_send_field("ptn_pftmean",ptn_pftmean)
879	CALL xios_orchidee_send_field("soilflx",soilflx)
880	CALL xios_orchidee_send_field("surfheat_incr",surfheat_incr)
881	CALL xios_orchidee_send_field("coldcont_incr",coldcont_incr)
882	CALL xios_orchidee_send_field("pkappa",pkappa)
883	CALL xios_orchidee_send_field("pkappa_snow",pkappa_snow_diag)
884	CALL xios_orchidee_send_field("pcapa",pcapa)
885	CALL xios_orchidee_send_field("pcapa_snow",pcapa_snow_diag)
886	CALL xios_orchidee_send_field("snowtemp",snowtemp_diag)
887	IF (ok_freeze_thermix) CALL xios_orchidee_send_field("shum_ngrnd_perma", shum_ngrnd_perma)
888
889	IF ( .NOT. almaoutput ) THEN
890	CALL histwrite_p(hist_id, 'ptn', kjit, ptn_pftmean, kjpindex*ngrnd, indexgrnd)
891	CALL histwrite_p(hist_id, 'Qg', kjit, soilflx, kjpindex, index)
892	CALL histwrite_p(hist_id, 'pkappa', kjit, pkappa, kjpindex*ngrnd, indexgrnd)
893	CALL histwrite_p(hist_id, 'pcapa', kjit, pcapa, kjpindex*ngrnd, indexgrnd)
894
895	IF (ok_freeze_thermix) THEN
896	CALL histwrite_p(hist_id, 'profil_froz', kjit, profil_froz, kjpindex*ngrnd, indexgrnd)
897	CALL histwrite_p(hist_id, 'pcappa_supp', kjit, pcappa_supp, kjpindex*ngrnd, indexgrnd)
898	END IF
899	CALL histwrite_p(hist_id, 'shum_ngrnd_perma', kjit, shum_ngrnd_perma(:,:), kjpindex*ngrnd, indexgrnd)
900
901	ELSE
902	CALL histwrite_p(hist_id, 'SoilTemp', kjit, ptn_pftmean, kjpindex*ngrnd, indexgrnd)
903	CALL histwrite_p(hist_id, 'Qg', kjit, soilflx, kjpindex, index)
904	CALL histwrite_p(hist_id, 'DelSurfHeat', kjit, surfheat_incr, kjpindex, index)
905	CALL histwrite_p(hist_id, 'DelColdCont', kjit, coldcont_incr, kjpindex, index)
906	ENDIF
907	IF ( hist2_id > 0 ) THEN
908	IF ( .NOT. almaoutput ) THEN
909	CALL histwrite_p(hist2_id, 'ptn', kjit, ptn_pftmean, kjpindex*ngrnd, indexgrnd)
910	ELSE
911	CALL histwrite_p(hist2_id, 'SoilTemp', kjit, ptn_pftmean, kjpindex*ngrnd, indexgrnd)
912	CALL histwrite_p(hist2_id, 'Qg', kjit, soilflx, kjpindex, index)
913	CALL histwrite_p(hist2_id, 'DelSurfHeat', kjit, surfheat_incr, kjpindex, index)
914	CALL histwrite_p(hist2_id, 'DelColdCont', kjit, coldcont_incr, kjpindex, index)
915	ENDIF
916	ENDIF
917
918	!! 7. Considering the heat released by microbial respiration
919	ok_zimov=.FALSE.
920	IF (ok_zimov) THEN
921	CALL thermosoil_add_heat_zimov(kjpindex, veget_max, ptn, heat_zimov)
922	END IF
923
924	!! 8. A last final call to thermosoil_coef
925
926	!! A last final call to thermosoil_coef, which calculates the different
927	!!coefficients (cgrnd, dgrnd, soilcap, soilflx) from this time step to be
928	!!used at the next time step, either in the surface temperature calculation
929	!!(soilcap, soilflx) or in the soil thermal numerical scheme.
930	CALL thermosoil_coef (&
931	kjpindex, temp_sol_new, snow, njsc, &
932	frac_snow_veg, frac_snow_nobio, totfrac_nobio, &
933	snowdz, snowrho, snowtemp, pb, &
934	veget_max, som_total, depth_organic_soil,&
935	ptn, ptn_pftmean, &
936	soilcap, soilflx, cgrnd, dgrnd,&
937	lambda_snow, cgrnd_snow, dgrnd_snow)
938
939	! Save variables for explicit snow model
940	gtemp(:) = ptn_pftmean(:,1)
941
942	!! Initialize output arguments to be used in sechiba
943	ptnlev1(:) = ptn_pftmean(:,1)
944
945	!! Initialize output arguments to be returned to sechiba and further used in stomate
946	DO jv= 1, nvm
947	deephum_prof(:,:,jv) = shum_ngrnd_perma(:,:)
948	END DO
949	deeptemp_prof = ptn
950
951	!! Surface temperature is forced to zero celcius if its value is larger than melting point, only for explicit snow scheme
952	DO ji=1,kjpindex
953	IF (SUM(snowdz(ji,:)) .GT. 0.0) THEN
954	IF (temp_sol_new(ji) .GE. tp_00) THEN
955	temp_sol_new(ji) = tp_00
956	ENDIF
957	END IF
958	END DO
959
960	IF (printlev>=3) WRITE (numout,*) ' thermosoil_main done '
961
962	END SUBROUTINE thermosoil_main
963
964	!! =============================================================================================================================
965	!! SUBROUTINE : thermosoil_finalize
966	!!
967	!>\BRIEF Write to restart file
968	!!
969	!! DESCRIPTION : This subroutine writes the module variables and variables calculated in thermosoil
970	!! to restart file
971	!! \n
972	!_ ==============================================================================================================================
973	SUBROUTINE thermosoil_finalize (kjit, kjpindex, rest_id, gtemp, &
974	soilcap, soilflx, lambda_snow, cgrnd_snow, dgrnd_snow)
975
976	!! 0. Variable and parameter declaration
977	!! 0.1 Input variables
978	INTEGER(i_std), INTENT(in) :: kjit !! Time step number (unitless)
979	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (unitless)
980	INTEGER(i_std),INTENT (in) :: rest_id !! Restart file identifier(unitless)
981	REAL(r_std),DIMENSION (kjpindex),INTENT(in) :: gtemp !! First soil layer temperature
982	REAL(r_std),DIMENSION (kjpindex),INTENT(in) :: soilcap
983	REAL(r_std),DIMENSION (kjpindex),INTENT(in) :: soilflx
984	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: lambda_snow !! Coefficient of the linear extrapolation of surface temperature
985	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT (in) :: cgrnd_snow !! Integration coefficient for snow numerical scheme
986	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT (in) :: dgrnd_snow !! Integration coefficient for snow numerical scheme
987
988	!! 0.2 Local variables
989	! To store variables names for I/O
990	CHARACTER(LEN=80) :: var_name
991	CHARACTER(LEN=10) :: part_str
992	INTEGER(i_std) :: iv
993	!_ ================================================================================================================================
994
995	!! 1. Write variables to restart file to be used for the next simulation
996	IF (printlev>=3) WRITE (numout,*) 'Write restart file with THERMOSOIL variables'
997
998	CALL restput_p(rest_id, 'ptn', nbp_glo, ngrnd, nvm, kjit, ptn, 'scatter', nbp_glo, index_g)
999
1000	CALL restput_p(rest_id, 'ptn_pftmean', nbp_glo, ngrnd, 1, kjit, ptn_pftmean, 'scatter', nbp_glo, index_g)
1001	CALL restput_p (rest_id, 'refsoc', nbp_glo, ngrnd, 1, kjit, refsoc, 'scatter', nbp_glo, index_g)
1002
1003	IF (ok_Ecorr) THEN
1004	CALL restput_p(rest_id, 'e_soil_lat', nbp_glo, 1 , 1, kjit, e_soil_lat, 'scatter', nbp_glo, index_g)
1005	END IF
1006
1007	CALL restput_p(rest_id, 'gtemp', nbp_glo, 1, 1, kjit, gtemp, 'scatter', nbp_glo, index_g)
1008	CALL restput_p(rest_id, 'soilcap', nbp_glo, 1, 1, kjit, soilcap, 'scatter', nbp_glo, index_g)
1009	CALL restput_p(rest_id, 'soilflx', nbp_glo, 1, 1, kjit, soilflx, 'scatter', nbp_glo, index_g)
1010	CALL restput_p(rest_id, 'cgrnd', nbp_glo, ngrnd-1, 1, kjit, cgrnd, 'scatter', nbp_glo, index_g)
1011	CALL restput_p(rest_id, 'dgrnd', nbp_glo, ngrnd-1, 1, kjit, dgrnd, 'scatter', nbp_glo, index_g)
1012	CALL restput_p(rest_id, 'cgrnd_snow', nbp_glo, nsnow, 1, kjit, cgrnd_snow, 'scatter', nbp_glo, index_g)
1013	CALL restput_p(rest_id, 'dgrnd_snow', nbp_glo, nsnow, 1, kjit, dgrnd_snow, 'scatter', nbp_glo, index_g)
1014	CALL restput_p(rest_id, 'lambda_snow', nbp_glo, 1, 1, kjit, lambda_snow, 'scatter', nbp_glo, index_g)
1015
1016	END SUBROUTINE thermosoil_finalize
1017
1018
1019	!! ================================================================================================================================
1020	!! SUBROUTINE : thermosoil_clear
1021	!!
1022	!>\BRIEF Deallocates the allocated arrays.
1023	!! The call of thermosoil_clear originates from sechiba_clear but the calling sequence and
1024	!! its purpose require further investigation.
1025	!!
1026	!! DESCRIPTION : None
1027	!!
1028	!! RECENT CHANGE(S) : None
1029	!!
1030	!! MAIN OUTPUT VARIABLE(S): None
1031	!!
1032	!! REFERENCE(S) : None
1033	!!
1034	!! FLOWCHART : None
1035	!! \n
1036	!_ ================================================================================================================================
1037
1038	SUBROUTINE thermosoil_clear()
1039
1040	IF ( ALLOCATED (ptn)) DEALLOCATE (ptn)
1041	IF ( ALLOCATED (ptn_pftmean)) DEALLOCATE (ptn_pftmean)
1042	IF ( ALLOCATED (cgrnd)) DEALLOCATE (cgrnd)
1043	IF ( ALLOCATED (dgrnd)) DEALLOCATE (dgrnd)
1044	IF ( ALLOCATED (pcapa)) DEALLOCATE (pcapa)
1045	IF ( ALLOCATED (pkappa)) DEALLOCATE (pkappa)
1046	IF ( ALLOCATED (pkappa_per_pft)) DEALLOCATE (pkappa_per_pft)
1047	IF ( ALLOCATED (pcapa_per_pft)) DEALLOCATE (pcapa_per_pft)
1048	IF ( ALLOCATED (pcapa_en_per_pft)) DEALLOCATE (pcapa_en_per_pft)
1049	IF ( ALLOCATED (pcapa_snow)) DEALLOCATE (pcapa_snow)
1050	IF ( ALLOCATED (pkappa_snow)) DEALLOCATE (pkappa_snow)
1051	IF ( ALLOCATED (pcapa_en)) DEALLOCATE (pcapa_en)
1052	IF ( ALLOCATED (ptn_beg)) DEALLOCATE (ptn_beg)
1053	IF ( ALLOCATED (temp_sol_beg)) DEALLOCATE (temp_sol_beg)
1054	IF ( ALLOCATED (surfheat_incr)) DEALLOCATE (surfheat_incr)
1055	IF ( ALLOCATED (coldcont_incr)) DEALLOCATE (coldcont_incr)
1056	IF ( ALLOCATED (shum_ngrnd_perma)) DEALLOCATE (shum_ngrnd_perma)
1057	IF ( ALLOCATED (profil_froz)) DEALLOCATE (profil_froz)
1058	IF ( ALLOCATED (mc_layt)) DEALLOCATE (mc_layt)
1059	IF ( ALLOCATED (mcl_layt)) DEALLOCATE (mcl_layt)
1060	IF ( ALLOCATED (tmc_layt)) DEALLOCATE (tmc_layt)
1061	IF ( ALLOCATED (mc_layt_pft)) DEALLOCATE (mc_layt_pft)
1062	IF ( ALLOCATED (mcl_layt_pft)) DEALLOCATE (mcl_layt_pft)
1063	IF ( ALLOCATED (tmc_layt_pft)) DEALLOCATE (tmc_layt_pft)
1064	IF ( ALLOCATED (profil_froz_pft)) DEALLOCATE (profil_froz_pft)
1065	IF ( ALLOCATED (refsoc)) DEALLOCATE (refsoc)
1066	END SUBROUTINE thermosoil_clear
1067
1068
1069	!! ================================================================================================================================
1070	!! SUBROUTINE : thermosoil_coef
1071	!!
1072	!>\BRIEF Calculate soil thermal properties, integration coefficients, apparent heat flux,
1073	!! surface heat capacity,
1074	!!
1075	!! DESCRIPTION : This routine computes : \n
1076	!! 1. the soil thermal properties. \n
1077	!! 2. the integration coefficients of the thermal numerical scheme, cgrnd and dgrnd,
1078	!! which depend on the vertical grid and on soil properties, and are used at the next
1079	!! timestep.\n
1080	!! 3. the soil apparent heat flux and surface heat capacity (soilflx
1081	!! and soilcap), used by enerbil to compute the surface temperature at the next
1082	!! timestep.\n
1083	!! - The soil thermal properties depend on water content (shum_ngrnd_perma, shumdiag_perma,
1084	!! mc_layt, mcl_layt, tmc_layt), dominant soil texture(njsc), and on the presence
1085	!! of snow : snow is integrated into the soil for the thermal calculations, ie if there
1086	!! is snow on the ground, the first thermal layer(s) consist in snow, depending on the
1087	!! snow-depth. If a layer consists out of snow and soil, wheighed soil properties are
1088	!! calculated\n
1089	!! - The coefficients cgrnd and dgrnd are the integration
1090	!! coefficients for the thermal scheme \n
1091	!! T(k+1)=cgrnd(k)+dgrnd(k)*T(k) \n
1092	!! -- EQ1 -- \n
1093	!! They correspond respectively to $\beta$ and $\alpha$ from F. Hourdin\'s thesis and
1094	!! their expression can be found in this document (eq A19 and A20)
1095	!! - soilcap and soilflx are the apparent surface heat capacity and flux
1096	!! used in enerbil at the next timestep to solve the surface
1097	!! balance for Ts (EQ3); they correspond to $C_s$ and $F_s$ in F.
1098	!! Hourdin\'s PhD thesis and are expressed in eq. A30 and A31. \n
1099	!! soilcap*(Ts(t)-Ts(t-1))/dt=soilflx+otherfluxes(Ts(t)) \n
1100	!! -- EQ3 --\n
1101	!!
1102	!! RECENT CHANGE(S) : None
1103	!!
1104	!! MAIN OUTPUT VARIABLE(S): cgrnd, dgrnd, pcapa, pkappa, soilcap, soilflx
1105	!!
1106	!! REFERENCE(S) :
1107	!! - Hourdin, F. (1992). Study and numerical simulation of the general circulation of planetary atmospheres,
1108	!! Ph.D. thesis, Paris VII University. Remark: the part of F. Hourdin's PhD thesis relative to the thermal
1109	!! integration scheme has been scanned and is provided along with the documentation, with name :
1110	!! Hourdin_1992_PhD_thermal_scheme.pdf
1111	!!
1112	!! FLOWCHART : None
1113	!! \n
1114	!_ ================================================================================================================================
1115
1116	SUBROUTINE thermosoil_coef (kjpindex, temp_sol_new, snow, njsc, &
1117	frac_snow_veg, frac_snow_nobio, totfrac_nobio, &
1118	snowdz, snowrho, snowtemp, pb, &
1119	veget_max, som_total, depth_organic_soil, &
1120	ptn, ptn_pftmean, &
1121	soilcap, soilflx, cgrnd, dgrnd,&
1122	lambda_snow, cgrnd_snow, dgrnd_snow)
1123
1124	!! 0. Variables and parameter declaration
1125
1126	!! 0.1 Input variables
1127
1128	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (unitless)
1129	REAL(r_std), DIMENSION (kjpindex), INTENT (in) :: temp_sol_new !! soil surface temperature @tex ($K$) @endtex
1130	REAL(r_std), DIMENSION (kjpindex), INTENT (in) :: snow !! snow mass @tex ($Kg$) @endtex
1131	INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: njsc !! Index of the dominant soil textural class in the grid cell (1-nscm, unitless)
1132	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: frac_snow_veg !! Snow cover fraction on vegeted area
1133	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT(in) :: frac_snow_nobio !! Snow cover fraction on non-vegeted area
1134	REAL(r_std),DIMENSION (kjpindex),INTENT(in) :: totfrac_nobio !! Total fraction of continental ice+lakes+cities+...
1135	!!(unitless,0-1)
1136	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(in) :: snowdz !! Snow depth (m)
1137	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(in) :: snowrho !! Snow density (Kg/m^3)
1138	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(in) :: snowtemp !! Snow temperature (K)
1139	REAL(r_std), DIMENSION (kjpindex), INTENT (in) :: pb !! Surface presure (hPa)
1140	REAL(r_std), DIMENSION (kjpindex,nvm), INTENT(in) :: veget_max !! Fraction of vegetation type
1141	REAL(r_std), DIMENSION (kjpindex,ngrnd,nvm,nelements), INTENT (in) :: som_total !! total soil carbon for use in thermal calcs (g/m**3)
1142	REAL(r_std), DIMENSION (kjpindex), INTENT (in) :: depth_organic_soil !! Depth at which there is still organic matter (m)
1143	!! 0.2 Output variables
1144
1145	REAL(r_std), DIMENSION (kjpindex), INTENT (out) :: soilcap !! surface heat capacity considering snow and soil surface
1146	!! @tex ($J m^{-2} K^{-1}$) @endtex
1147	REAL(r_std), DIMENSION (kjpindex), INTENT (out) :: soilflx !! surface heat flux considering snow and soil surface @tex ($W m^{-2}$) @endtex,
1148	!! positive towards the
1149	!! soil, writen as Qg (ground heat flux) in the history
1150	!! files.
1151	REAL(r_std), DIMENSION (kjpindex,ngrnd-1), INTENT(out) :: cgrnd !! matrix coefficient for the computation of soil
1152	!! temperatures (beta in F. Hourdin thesis)
1153	REAL(r_std), DIMENSION (kjpindex,ngrnd-1), INTENT(out) :: dgrnd !! matrix coefficient for the computation of soil
1154	!! temperatures (alpha in F. Hourdin thesis)
1155
1156
1157	!! 0.3 Modified variable
1158
1159	REAL(r_std), DIMENSION (kjpindex,ngrnd,nvm), INTENT (inout):: ptn !! vertically discretized soil temperatures per pft. ptn is only modified if ok_Ecorr.
1160	REAL(r_std), DIMENSION (kjpindex,ngrnd), INTENT (inout):: ptn_pftmean !! vertically discretized soil temperatures. ptn is only modified if ok_Ecorr.
1161	REAL(r_std), DIMENSION (kjpindex), INTENT(inout) :: lambda_snow !! Coefficient of the linear extrapolation of surface temperature
1162	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT (inout):: cgrnd_snow !! Integration coefficient for snow numerical scheme
1163	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT (inout):: dgrnd_snow !! Integration coefficient for snow numerical scheme
1164
1165	!! 0.4 Local variables
1166
1167	INTEGER(i_std) :: ji, jg, m
1168	REAL(r_std), DIMENSION (kjpindex,ngrnd-1) :: zdz1 !! numerical (buffer) constant
1169	!! @tex ($W m^{-1} K^{-1}$) @endtex
1170	REAL(r_std), DIMENSION (kjpindex,ngrnd) :: zdz2 !! numerical (buffer) constant
1171	!! @tex ($W m^{-1} K^{-1}$) @endtex
1172	REAL(r_std), DIMENSION (kjpindex) :: z1 !! numerical constant @tex ($W m^{-1} K^{-1}$) @endtex
1173	REAL(r_std), DIMENSION (kjpindex) :: soilcap_nosnow !! surface heat capacity
1174	!! @tex ($J m^{-2} K^{-1}$)
1175	!! @endtex
1176	REAL(r_std), DIMENSION (kjpindex) :: soilflx_nosnow !! surface heat flux @tex ($W m^{-2}$) @endtex,
1177	!! positive towards the soil, written as Qg
1178	!!(ground heat flux in the history files).
1179	REAL(r_std), DIMENSION (kjpindex) :: snowcap !! apparent snow heat capacity @tex ($J m^{-2} K^{-1}$)
1180	REAL(r_std), DIMENSION (kjpindex) :: snowflx !! apparent snow-atmosphere heat flux @tex ($W m^{-2}$) @endtex
1181	REAL(r_std), DIMENSION (kjpindex,nsnow) :: dz1_snow
1182	REAL(r_std), DIMENSION (kjpindex,nsnow) :: ZSNOWDZM
1183	REAL(r_std), DIMENSION (kjpindex,nsnow) :: dz2_snow
1184	REAL(r_std), DIMENSION (kjpindex,nsnow) :: zdz1_snow
1185	REAL(r_std), DIMENSION (kjpindex,nsnow) :: zdz2_snow
1186	REAL(r_std), DIMENSION (kjpindex) :: z1_snow
1187	REAL(r_std), DIMENSION (kjpindex) :: snowflxtot !! Total snow flux (including snow on vegetated and bare soil and nobio areas)
1188	!! @tex ($W m^{-2}$) @endtex
1189	!! positive towards the soil
1190
1191	!_ ================================================================================================================================
1192
1193	!! 1. Computation of the soil thermal properties
1194
1195	! Computation of the soil thermal properties; snow properties are also accounted for
1196	IF (ok_freeze_thermix) THEN
1197	IF (ok_soil_carbon_discretization) THEN
1198	pcapa(:,:) = zero
1199	pcapa_en(:,:) = zero
1200	pkappa(:,:) = zero
1201
1202	CALL thermosoil_getdiff_pft( kjpindex, ptn, njsc, shum_ngrnd_perma, &
1203	depth_organic_soil, som_total, snowrho, snowtemp, pb, veget_max)
1204
1205	DO m=1,nvm
1206	DO jg = 1, ngrnd
1207	pcapa(:,jg) = pcapa(:,jg) + pcapa_per_pft(:,jg,m) * veget_max(:,m)
1208	pcapa_en(:,jg) = pcapa_en(:,jg) + pcapa_en_per_pft(:,jg,m) * veget_max(:,m)
1209	pkappa(:,jg) = pkappa(:,jg) + pkappa_per_pft(:,jg,m) * veget_max(:,m)
1210	END DO
1211	END DO
1212	ELSE
1213	CALL thermosoil_getdiff( kjpindex, snow, ptn_pftmean, njsc, snowrho, snowtemp, pb )
1214	ENDIF
1215	ELSE
1216	! Special case without soil freezing
1217	CALL thermosoil_getdiff_old_thermix_without_snow( kjpindex, njsc, snowrho, snowtemp, pb )
1218	ENDIF
1219
1220	! Energy conservation : Correction to make sure that the same latent heat is released and
1221	! consumed during freezing and thawing
1222	IF (ok_Ecorr) THEN
1223	CALL thermosoil_readjust(kjpindex, ptn, ptn_pftmean, veget_max)
1224	ENDIF
1225
1226
1227	!! 2. Computation of the coefficients of the numerical integration scheme for the soil layers
1228
1229	!! 2.1 Calculate numerical coefficients zdz1 and zdz2
1230	DO jg=1,ngrnd
1231	DO ji=1,kjpindex
1232	zdz2(ji,jg)=pcapa(ji,jg) * dlt(jg)/dt_sechiba
1233	ENDDO
1234	ENDDO
1235
1236	DO jg=1,ngrnd-1
1237	DO ji=1,kjpindex
1238	zdz1(ji,jg) = dz1(jg) * pkappa(ji,jg)
1239	ENDDO
1240	ENDDO
1241
1242	!! 2.2 Calculate coefficients cgrnd and dgrnd for soil
1243	DO ji = 1,kjpindex
1244	z1(ji) = zdz2(ji,ngrnd) + zdz1(ji,ngrnd-1)
1245	cgrnd(ji,ngrnd-1) = zdz2(ji,ngrnd) * ptn_pftmean(ji,ngrnd) / z1(ji)
1246	dgrnd(ji,ngrnd-1) = zdz1(ji,ngrnd-1) / z1(ji)
1247	ENDDO
1248
1249	DO jg = ngrnd-1,2,-1
1250	DO ji = 1,kjpindex
1251	z1(ji) = un / (zdz2(ji,jg) + zdz1(ji,jg-1) + zdz1(ji,jg) * (un - dgrnd(ji,jg)))
1252	cgrnd(ji,jg-1) = (ptn_pftmean(ji,jg) * zdz2(ji,jg) + zdz1(ji,jg) * cgrnd(ji,jg)) * z1(ji)
1253	dgrnd(ji,jg-1) = zdz1(ji,jg-1) * z1(ji)
1254	ENDDO
1255	ENDDO
1256
1257
1258	!! 3. Computation of the coefficients of the numerical integration scheme for the snow layers
1259
1260	!! 3.1 Calculate numerical coefficients zdz1_snow, zdz2_snow and lambda_snow
1261	DO ji = 1, kjpindex
1262
1263	! Calculate internal values
1264	DO jg = 1, nsnow
1265	ZSNOWDZM(ji,jg) = MAX(snowdz(ji,jg),psnowdzmin)
1266	ENDDO
1267	dz2_snow(ji,:)=ZSNOWDZM(ji,:)
1268
1269	DO jg = 1, nsnow-1
1270	dz1_snow(ji,jg) = 2.0 / (dz2_snow(ji,jg+1)+dz2_snow(ji,jg))
1271	ENDDO
1272
1273	lambda_snow(ji) = dz2_snow(ji,1)/2.0 * dz1_snow(ji,1)
1274
1275	DO jg=1,nsnow
1276	zdz2_snow(ji,jg)=pcapa_snow(ji,jg) * dz2_snow(ji,jg)/dt_sechiba
1277	ENDDO
1278
1279	DO jg=1,nsnow-1
1280	zdz1_snow(ji,jg) = dz1_snow(ji,jg) * pkappa_snow(ji,jg)
1281	ENDDO
1282
1283	! the bottom snow
1284	zdz1_snow(ji,nsnow) = pkappa_snow(ji,nsnow) / ( zlt(1) + dz2_snow(ji,nsnow)/2 )
1285
1286	ENDDO
1287
1288
1289	!! 3.2 Calculate coefficients cgrnd_snow and dgrnd_snow for snow
1290	DO ji = 1,kjpindex
1291	! bottom level
1292	z1_snow(ji) = zdz2(ji,1)+(un-dgrnd(ji,1))*zdz1(ji,1)+zdz1_snow(ji,nsnow)
1293	cgrnd_snow(ji,nsnow) = (zdz2(ji,1) * ptn_pftmean(ji,1) + zdz1(ji,1) * cgrnd(ji,1) ) / z1_snow(ji)
1294	dgrnd_snow(ji,nsnow) = zdz1_snow(ji,nsnow) / z1_snow(ji)
1295
1296	! next-to-bottom level
1297	z1_snow(ji) = zdz2_snow(ji,nsnow)+(un-dgrnd_snow(ji,nsnow))*zdz1_snow(ji,nsnow)+zdz1_snow(ji,nsnow-1)
1298	cgrnd_snow(ji,nsnow-1) = (zdz2_snow(ji,nsnow)*snowtemp(ji,nsnow)+&
1299	zdz1_snow(ji,nsnow)*cgrnd_snow(ji,nsnow))/z1_snow(ji)
1300	dgrnd_snow(ji,nsnow-1) = zdz1_snow(ji,nsnow-1) / z1_snow(ji)
1301
1302	DO jg = nsnow-1,2,-1
1303	z1_snow(ji) = un / (zdz2_snow(ji,jg) + zdz1_snow(ji,jg-1) + zdz1_snow(ji,jg) * (un - dgrnd_snow(ji,jg)))
1304	cgrnd_snow(ji,jg-1) = (snowtemp(ji,jg) * zdz2_snow(ji,jg) + zdz1_snow(ji,jg) * cgrnd_snow(ji,jg)) * z1_snow(ji)
1305	dgrnd_snow(ji,jg-1) = zdz1_snow(ji,jg-1) * z1_snow(ji)
1306	ENDDO
1307	ENDDO
1308
1309
1310
1311	!! 4. Computation of the apparent ground heat flux
1312	!! Computation of apparent snow-atmosphere flux
1313	DO ji = 1,kjpindex
1314	snowflx(ji) = zdz1_snow(ji,1) * (cgrnd_snow(ji,1) + (dgrnd_snow(ji,1)-1.) * snowtemp(ji,1))
1315	snowcap(ji) = (zdz2_snow(ji,1) * dt_sechiba + dt_sechiba * (un - dgrnd_snow(ji,1)) * zdz1_snow(ji,1))
1316	z1_snow(ji) = lambda_snow(ji) * (un - dgrnd_snow(ji,1)) + un
1317	snowcap(ji) = snowcap(ji) / z1_snow(ji)
1318	snowflx(ji) = snowflx(ji) + &
1319	snowcap(ji) * (snowtemp(ji,1) * z1_snow(ji) - lambda_snow(ji) * cgrnd_snow(ji,1) - temp_sol_new(ji)) / dt_sechiba
1320	ENDDO
1321
1322
1323	!! Computation of the apparent ground heat flux (> towards the soil) and
1324	!! apparent surface heat capacity, used at the next timestep by enerbil to
1325	!! compute the surface temperature.
1326	DO ji = 1,kjpindex
1327	soilflx_nosnow(ji) = zdz1(ji,1) * (cgrnd(ji,1) + (dgrnd(ji,1)-1.) * ptn_pftmean(ji,1))
1328	soilcap_nosnow(ji) = (zdz2(ji,1) * dt_sechiba + dt_sechiba * (un - dgrnd(ji,1)) * zdz1(ji,1))
1329	z1(ji) = lambda * (un - dgrnd(ji,1)) + un
1330	soilcap_nosnow(ji) = soilcap_nosnow(ji) / z1(ji)
1331	soilflx_nosnow(ji) = soilflx_nosnow(ji) + &
1332	& soilcap_nosnow(ji) * (ptn_pftmean(ji,1) * z1(ji) - lambda * cgrnd(ji,1) - temp_sol_new(ji)) / dt_sechiba
1333
1334	ENDDO
1335
1336	!! Add snow fraction
1337	! Using an effective heat capacity and heat flux by a simple pondering of snow and soil fraction
1338	DO ji = 1, kjpindex
1339	soilcap(ji) = snowcap(ji)frac_snow_veg(ji)(1-totfrac_nobio(ji))+ & ! weights related to snow cover fraction on vegetation
1340	soilcap_nosnow(ji)SUM(frac_snow_nobio(ji,:))totfrac_nobio(ji)+ & ! weights related to SCF on nobio
1341	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
1342	soilflx(ji) = snowflx(ji)frac_snow_veg(ji)(1-totfrac_nobio(ji))+ & ! weights related to snow cover fraction on vegetation
1343	soilflx_nosnow(ji)SUM(frac_snow_nobio(ji,:))totfrac_nobio(ji)+ & ! weights related to SCF on nobio
1344	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
1345	ENDDO
1346
1347	! Total snow flux (including snow on vegetated and bare soil and nobio areas)
1348	DO ji = 1, kjpindex
1349	snowflxtot(ji) = snowflx(ji)frac_snow_veg(ji)(1-totfrac_nobio(ji)) + &
1350	soilflx_nosnow(ji)SUM(frac_snow_nobio(ji,:))totfrac_nobio(ji)
1351	ENDDO
1352
1353	CALL xios_orchidee_send_field("snowflxtot",snowflxtot(:))
1354
1355	IF (printlev>=3) WRITE (numout,*) ' thermosoil_coef done '
1356
1357	END SUBROUTINE thermosoil_coef
1358
1359
1360	!! ================================================================================================================================
1361	!! SUBROUTINE : thermosoil_profile
1362	!!
1363	!>\BRIEF In this routine solves the numerical soil thermal scheme, ie calculates the new soil temperature profile;
1364	!!
1365	!!
1366	!! DESCRIPTION : The calculation of the new soil temperature profile is based on
1367	!! the cgrnd and dgrnd values from the previous timestep and the surface temperature Ts aka temp_sol_new. (see detailed
1368	!! explanation in the header of the thermosoil module or in the reference).\n
1369	!! T(k+1)=cgrnd(k)+dgrnd(k)*T(k)\n
1370	!! -- EQ1 --\n
1371	!! Ts=(1+lambda)T(1) -lambdaT(2)\n
1372	!! -- EQ2--\n
1373	!!
1374	!! RECENT CHANGE(S) : None
1375	!!
1376	!! MAIN OUTPUT VARIABLE(S): ptn (soil temperature profile on the thermal axis),
1377	!! stempdiag (soil temperature profile on the diagnostic axis)
1378	!!
1379	!! REFERENCE(S) :
1380	!! - Hourdin, F. (1992). Study and numerical simulation of the general circulation of planetary atmospheres,
1381	!! Ph.D. thesis, Paris VII University. Remark: the part of F. Hourdin's PhD thesis relative to the thermal
1382	!! integration scheme has been scanned and is provided along with the documentation, with name :
1383	!! Hourdin_1992_PhD_thermal_scheme.pdf
1384	!!
1385	!! FLOWCHART : None
1386	!! \n
1387	!_ ================================================================================================================================
1388
1389	SUBROUTINE thermosoil_profile (kjpindex, temp_sol_new, &
1390	frac_snow_veg, frac_snow_nobio, totfrac_nobio, &
1391	ptn, ptn_pftmean, stempdiag, snowtemp, &
1392	cgrnd_snow, dgrnd_snow)
1393
1394	!! 0. Variables and parameter declaration
1395
1396	!! 0.1 Input variables
1397
1398	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (unitless)
1399	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_sol_new !! Surface temperature at the present time-step
1400	!! @tex ($K$) @endtex
1401	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: frac_snow_veg !! Snow cover fraction on vegeted area
1402	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT(in) :: frac_snow_nobio!! Snow cover fraction on non-vegeted area
1403	REAL(r_std),DIMENSION (kjpindex),INTENT(in) :: totfrac_nobio !! Total fraction of continental ice+lakes+cities+...
1404	!! (unitless,0-1)
1405	REAL(r_std),DIMENSION (kjpindex,nsnow), INTENT(in) :: snowtemp !! Snow temperature (K)
1406	REAL(r_std),DIMENSION (kjpindex,nsnow), INTENT(in) :: cgrnd_snow !! Integration coefficient for snow numerical scheme
1407	REAL(r_std),DIMENSION (kjpindex,nsnow), INTENT(in) :: dgrnd_snow !! Integration coefficient for snow numerical scheme
1408
1409	!! 0.3 Modified variables
1410
1411
1412	!! 0.2 Output variables
1413	REAL(r_std),DIMENSION (kjpindex,ngrnd,nvm), INTENT (out) :: ptn !! vertically discretized soil temperatures per pft (K)
1414	!! @tex ($K$) @endtex
1415	REAL(r_std),DIMENSION (kjpindex,ngrnd), INTENT (out) :: ptn_pftmean !! vertically discretized soil temperatures (K)
1416	REAL(r_std),DIMENSION (kjpindex,nslm), INTENT (out) :: stempdiag !! diagnostic temperature profile
1417	!! @tex ($K$) @endtex
1418
1419	!! 0.4 Local variables
1420	INTEGER(i_std) :: ji, jg, jv
1421	REAL(r_std) :: temp_sol_eff !! effective surface temperature including snow and soil
1422
1423	!_ ================================================================================================================================
1424
1425	!! 1. Computes the soil temperatures ptn.
1426
1427	!! 1.1. ptn(jg=1) using EQ1 and EQ2
1428	DO ji = 1,kjpindex
1429
1430	! Using an effective surface temperature by a simple pondering
1431	temp_sol_eff=snowtemp(ji,nsnow)frac_snow_veg(ji)(1-totfrac_nobio(ji))+ & ! weights related to snow cover fraction on vegetation
1432	temp_sol_new(ji)SUM(frac_snow_nobio(ji,:))totfrac_nobio(ji)+ & ! weights related to SCF on nobio
1433	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
1434	! Soil temperature calculation with explicit snow if there is snow on the ground
1435	ptn_pftmean(ji,1) = cgrnd_snow(ji,nsnow) + dgrnd_snow(ji,nsnow) * temp_sol_eff
1436	ENDDO
1437
1438	!! 1.2. ptn_pftmean(jg=2:ngrnd) using EQ1.
1439	DO jg = 1,ngrnd-1
1440	DO ji = 1,kjpindex
1441	ptn_pftmean(ji,jg+1) = cgrnd(ji,jg) + dgrnd(ji,jg) * ptn_pftmean(ji,jg)
1442	ENDDO
1443	ENDDO
1444
1445	!! Calculate ptn per pft
1446	! Here ptn is the same at all pfts when ok_soil_carbon_discretization is activated or not.
1447	! ptn is modified in thermosoil_add_heat_zimov if ok_zimov=TRUE.
1448	DO jv=1,nvm
1449	ptn(:,:,jv) = ptn_pftmean(:,:)
1450	END DO
1451
1452
1453	!! 2. Assigne the soil temperature to the output variable. It is already on the right axis.
1454	stempdiag(:,:) = ptn_pftmean(:,1:nslm)
1455
1456	IF (printlev>=3) WRITE (numout,*) ' thermosoil_profile done '
1457
1458	END SUBROUTINE thermosoil_profile
1459
1460	!================================================================================================================================
1461	!! SUBROUTINE : thermosoil_cond
1462	!!
1463	!>\BRIEF Calculate soil thermal conductivity.
1464	!!
1465	!! DESCRIPTION : This routine computes soil thermal conductivity
1466	!! Code introduced from NOAH LSM.
1467	!!
1468	!! RECENT CHANGE(S) : None
1469	!!
1470	!! MAIN OUTPUT VARIABLE(S): cnd
1471	!!
1472	!! REFERENCE(S) :
1473	!! Farouki, O.T.,1986: Thermal Properties of Soils. Series on Rock
1474	!! and Soil Mechanics, Vol. 11, Trans Tech, 136 PP.
1475	!! Johansen, O., 1975: Thermal Conductivity of Soils. Ph.D. Thesis,
1476	!! University of Trondheim,
1477	!! Peters-Lidard, C. D., Blackburn, E., Liang, X., & Wood, E. F.,
1478	!! 1998: The effect of soil thermal conductivity
1479	!! Parameterization on Surface Energy fluxes
1480	!! and Temperatures. J. of The Atmospheric Sciences,
1481	!! Vol. 55, pp. 1209-1224.
1482	!! Modify histroy:
1483	!!
1484	!! FLOWCHART : None
1485	!! \n
1486	!_
1487	!================================================================================================================================
1488
1489	SUBROUTINE thermosoil_cond (kjpindex, njsc, smc, qz, sh2o, cnd)
1490
1491	!! 0. Variables and parameter declaration
1492
1493	!! 0.1 Input variables
1494	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (unitless)
1495	INTEGER(i_std), DIMENSION (kjpindex), INTENT (in) :: njsc !! Index of the dominant soil textural class in the grid cell (1-nscm, unitless)
1496	REAL(r_std), DIMENSION (kjpindex,ngrnd), INTENT(IN) :: smc !! Volumetric Soil Moisture Content (m3/m3)
1497	REAL(r_std), DIMENSION (nscm), INTENT(IN) :: qz !! Quartz Content (Soil Type Dependent) (0-1)
1498	REAL(r_std), DIMENSION (kjpindex,ngrnd), INTENT(IN) :: sh2o !! Unfrozen Soil Moisture Content; Frozen Soil Moisture = smc - sh2o
1499
1500	!! 0.2 Output variables
1501	REAL(r_std), DIMENSION (kjpindex,ngrnd), INTENT(OUT) :: cnd !! Soil Thermal Conductivity (W/m/k)
1502
1503	!! 0.3 Modified variables
1504
1505	!! 0.4 Local variables
1506	REAL(r_std), DIMENSION (kjpindex,ngrnd) :: ake !! Kersten Number (unitless)
1507	REAL(r_std), DIMENSION (kjpindex,ngrnd) :: thksat !! Saturated Thermal Conductivity (W/m/k)
1508	REAL(r_std), DIMENSION (kjpindex,ngrnd) :: satratio !! Degree of Saturation (0-1)
1509	REAL(r_std), DIMENSION (kjpindex,ngrnd) :: xu !! Unfrozen Volume For Saturation (0-1)
1510	REAL(r_std), DIMENSION (kjpindex,ngrnd) :: xunfroz !! Unfrozon Volume Fraction (0-1)
1511	REAL(r_std) :: thko !! Thermal Conductivity for Other Ssoil Components (W/m/k)
1512	REAL(r_std) :: gammd !! Dry Dendity (kg/m3)
1513	REAL(r_std) :: thkdry !! Dry Thermal Conductivity (W/m/k)
1514	REAL(r_std) :: thks !! Thermal Conductivity for the Solids Combined (Quartz + Other) (W/m/k)
1515	REAL(r_std), PARAMETER :: THKICE = 2.2 !! Ice Thermal Conductivity (W/m/k)
1516	REAL(r_std), PARAMETER :: THKQTZ = 7.7 !! Thermal Conductivity for Quartz (W/m/k)
1517	REAL(r_std), PARAMETER :: THKW = 0.57 !! Water Thermal Conductivity (W/m/k)
1518	INTEGER(i_std) :: ji, jg, jst
1519
1520	!_================================================================================================================================
1521
1522	!! 1. Dry and Saturated Thermal Conductivity.
1523
1524	DO ji = 1,kjpindex
1525	jst = njsc(ji)
1526
1527	!! 1.1. Dry density (Kg/m3) and Dry thermal conductivity (W.M-1.K-1)
1528	gammd = (1. - mcs(jst))*2700.
1529	thkdry = (0.135* gammd+ 64.7)/ (2700. - 0.947* gammd)
1530
1531	!! 1.2. thermal conductivity of "other" soil components
1532	IF (qz(jst) > 0.2) THEN
1533	thko = 2.0
1534	ELSEIF (qz(jst) <= 0.2) THEN
1535	thko = 3.0
1536	ENDIF
1537
1538	!! 1.3. Thermal conductivity of solids
1539	thks = (THKQTZ ** qz(jst))* (thko ** (1. - qz(jst)))
1540
1541	DO jg = 1,ngrnd
1542	!! 1.4. saturation ratio
1543	satratio(ji,jg) = smc(ji,jg) / mcs(jst)
1544
1545	!! 1.5. Saturated Thermal Conductivity (thksat)
1546	IF ( smc(ji,jg) > min_sechiba ) THEN
1547	xunfroz(ji,jg) = sh2o(ji,jg) / smc(ji,jg) ! Unfrozen Fraction (From i.e., 100%Liquid, to 0. (100% Frozen))
1548	xu(ji,jg) = xunfroz(ji,jg) * mcs(jst) ! Unfrozen volume for saturation (porosity*xunfroz)
1549	thksat(ji,jg) = thks ** (1. - mcs(jst))* THKICE ** (mcs(jst) - xu(ji,jg))* THKW ** (xu(ji,jg))
1550	ELSE
1551	! this value will not be used since ake=0 for this case
1552	thksat(ji,jg)=0
1553	END IF
1554	END DO ! DO jg = 1,ngrnd
1555
1556	!! 2. Kersten Number (ake)
1557	DO jg = 1,ngrnd
1558	IF ( (sh2o(ji,jg) + 0.0005) < smc(ji,jg) ) THEN
1559	! Frozen
1560	ake(ji,jg) = satratio(ji,jg)
1561	ELSE
1562	! Unfrozen
1563	! Eq 11 in Peters-Lidard et al., 1998
1564	IF ( satratio(ji,jg) > 0.1 ) THEN
1565	IF ((jst < 4 .AND. soil_classif == 'usda') .OR. (jst == 1 .AND. soil_classif == 'zobler') ) THEN
1566	! Coarse
1567	ake(ji,jg) = 0.7 * LOG10 (SATRATIO(ji,jg)) + 1.0
1568	ELSE
1569	! Fine
1570	ake(ji,jg) = LOG10 (satratio(ji,jg)) + 1.0
1571	ENDIF
1572	ELSEIF ( satratio(ji,jg) > 0.05 .AND. satratio(ji,jg) <= 0.1 ) THEN
1573	IF ((jst < 4 .AND. soil_classif == 'usda') .OR. (jst == 1 .AND. soil_classif == 'zobler') ) THEN
1574	! Coarse
1575	ake(ji,jg) = 0.7 * LOG10 (satratio(ji,jg)) + 1.0
1576	ELSE
1577	! Fine
1578	ake(ji,jg) = 0.0
1579	ENDIF
1580	ELSE
1581	ake(ji,jg) = 0.0 ! use k = kdry
1582	END IF
1583	END IF
1584	END DO ! DO jg = 1,ngrnd
1585
1586	!! 3. Thermal conductivity (cnd)
1587	DO jg = 1,ngrnd
1588	cnd(ji,jg) = ake(ji,jg) * (thksat(ji,jg) - thkdry) + thkdry
1589	END DO ! DO jg = 1,ngrnd
1590
1591	END DO !DO ji = 1,kjpindex
1592
1593	END SUBROUTINE thermosoil_cond
1594
1595	!================================================================================================================================
1596	!! SUBROUTINE : thermosoil_cond_pft
1597	!!
1598	!>\BRIEF Calculate soil thermal conductivity.
1599	!!
1600	!! DESCRIPTION : This routine computes soil thermal conductivity
1601	!! but considers the fact that soil organic carbon can decrease
1602	!! conductivity
1603	!!
1604	!! RECENT CHANGE(S) : None
1605	!!
1606	!! MAIN OUTPUT VARIABLE(S): cnd
1607	!!
1608	!! REFERENCE(S) :
1609	!! Farouki, O.T.,1986: Thermal Properties of Soils. Series on Rock
1610	!! and Soil Mechanics, Vol. 11, Trans Tech, 136 PP.
1611	!! Johansen, O., 1975: Thermal Conductivity of Soils. Ph.D. Thesis,
1612	!! University of Trondheim,
1613	!! Peters-Lidard, C. D., Blackburn, E., Liang, X., & Wood, E. F.,
1614	!! 1998: The effect of soil thermal conductivity
1615	!! Parameterization on Surface Energy fluxes
1616	!! and Temperatures. J. of The Atmospheric Sciences,
1617	!! Vol. 55, pp. 1209-1224.
1618	!! Lawrence and Slater,2008: Incorporating organic soil into a global climate
1619	!! model
1620	!! Modify histroy:
1621	!!
1622	!! FLOWCHART : None
1623	!! \n
1624	!_
1625	!================================================================================================================================
1626
1627
1628	SUBROUTINE thermosoil_cond_pft (kjpindex, njsc, smc, qz, sh2o,zx1,zx2,porosnet,cnd)
1629
1630	!! 0. Variables and parameter declaration
1631
1632	!! 0.1 Input variables
1633	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (unitless)
1634	INTEGER(i_std), DIMENSION (kjpindex), INTENT (in) :: njsc !! Index of the dominant soil textural class in the grid cell (1-nscm, unitless)
1635	REAL(r_std), DIMENSION (kjpindex,ngrnd), INTENT(IN) :: smc !! Volumetric Soil Moisture Content (m3/m3)
1636	REAL(r_std), DIMENSION (nscm), INTENT(IN) :: qz !! Quartz Content (Soil Type Dependent) (0-1)
1637	REAL(r_std), DIMENSION (kjpindex,ngrnd,nvm), INTENT(IN) :: porosnet !! Soil Porosity (0-1)
1638	REAL(r_std), DIMENSION (kjpindex,ngrnd), INTENT(IN) :: sh2o !! Unfrozen Soil Moisture Content; Frozen Soil Moisture = smc - sh2o
1639	REAL(r_std), DIMENSION(kjpindex,ngrnd,nvm), INTENT(IN) :: zx1, zx2 !! proportion of organic and mineral soil
1640
1641	!! 0.2 Output variables
1642	REAL(r_std), DIMENSION (kjpindex,ngrnd,nvm), INTENT(OUT) :: cnd !! Soil Thermal Conductivity (W/m/k)
1643
1644	!! 0.3 Modified variables
1645
1646	!! 0.4 Local variables
1647	REAL(r_std), DIMENSION (kjpindex,ngrnd,nvm) :: ake !! Kerston Number (unitless)
1648	REAL(r_std), DIMENSION (kjpindex,ngrnd,nvm) :: thksat !! Saturated Thermal Conductivity (W/m/k)
1649	REAL(r_std), DIMENSION (kjpindex,ngrnd,nvm) :: satratio !! Degree of Saturation (0-1)
1650	REAL(r_std), DIMENSION (kjpindex,ngrnd,nvm) :: xu !! Unfrozen Volume For Saturation (0-1)
1651	REAL(r_std), DIMENSION (kjpindex,ngrnd,nvm) :: xunfroz !! Unfrozon Volume Fraction (0-1)
1652	REAL(r_std) :: thko !! Thermal Conductivity for Other Ssoil Components (W/m/k)
1653	REAL(r_std), DIMENSION (kjpindex) :: gammd !! Dry Density (kg/m3)
1654	REAL(r_std), DIMENSION (kjpindex) :: thkdry_min !! Dry Thermal Conductivity for mineral soil (W/m/k)
1655	REAL(r_std), DIMENSION (kjpindex,ngrnd,nvm) :: thkdry !! Dry Thermal Conductivity considering organic carbon (W/m/k)
1656	REAL(r_std), DIMENSION (kjpindex) :: thks_min !! Thermal Conductivity for the Solids Combined (Quartz + Other) (W/m/k)
1657	REAL(r_std), DIMENSION (kjpindex,ngrnd,nvm) :: thks !! Thermal Conductivity considering organic carbon (W/m/k)
1658	INTEGER(i_std) :: ji, jg, jst, jv
1659	REAL(r_std), PARAMETER :: THKICE = 2.2 !! Ice Thermal Conductivity (W/m/k)
1660	REAL(r_std), PARAMETER :: THKQTZ = 7.7 !! Thermal Conductivity for Quartz (W/m/k)
1661	REAL(r_std), PARAMETER :: THKW = 0.57 !! Water Thermal Conductivity (W/m/k)
1662	!_================================================================================================================================
1663
1664	thksat(:,:,:)=0
1665	!! 1. Dry and Saturated Thermal Conductivity.
1666
1667	DO ji = 1,kjpindex
1668	jst = njsc(ji)
1669
1670	!! 1.1. Dry density (Kg/m3) and Dry thermal conductivity (W.M-1.K-1)
1671	gammd(ji) = (1. - mcs(jst))*2700.
1672	thkdry_min(ji) = (0.135* gammd(ji) + 64.7)/ (2700. - 0.947* gammd(ji))
1673
1674
1675	!! 1.2. thermal conductivity of "other" soil components
1676	IF (qz(jst) > 0.2) THEN
1677	thko = 2.0
1678	ELSEIF (qz(jst) <= 0.2) THEN
1679	thko = 3.0
1680	ENDIF
1681
1682	!! 1.3. Thermal conductivity of solids
1683	thks_min(ji) = (THKQTZ ** qz(jst))* (thko ** (1. - qz(jst)))
1684	ENDDO
1685
1686	DO jv = 1,nvm
1687
1688	SELECTCASE (use_soilc_method)
1689	CASE (SOILC_METHOD_ARITHMETIC)
1690	DO jg = 1,ngrnd
1691	DO ji = 1,kjpindex
1692	thks(ji,jg,jv) = zx1(ji,jg,jv) * cond_solid_org + zx2(ji,jg,jv) * thks_min(ji)
1693	ENDDO
1694	ENDDO
1695	CASE (SOILC_METHOD_GEOMETRIC)
1696	DO jg = 1,ngrnd
1697	DO ji = 1,kjpindex
1698	! use geometric mean rather than arithmetic mean (Decharme et al 2016)
1699	thks(ji,jg,jv) =(cond_solid_org*zx1(ji,jg,jv)) (thks_min(ji)**zx2(ji,jg,jv))
1700	ENDDO
1701	ENDDO
1702	ENDSELECT
1703
1704	DO jg = 1,ngrnd
1705	!! 1.4. saturation ratio
1706	DO ji = 1,kjpindex
1707	satratio(ji,jg,jv) = smc(ji,jg) / porosnet(ji, jg, jv)
1708	ENDDO
1709
1710	!! 1.5. Saturated Thermal Conductivity (thksat)
1711	DO ji = 1,kjpindex
1712	IF ( smc(ji,jg) > min_sechiba ) THEN
1713	xunfroz(ji,jg,jv) = sh2o(ji,jg) / smc(ji,jg) ! Unfrozen Fraction (From i.e., 100%Liquid, to 0. (100% Frozen))
1714	xu(ji,jg,jv) = xunfroz(ji,jg,jv) * porosnet(ji, jg, jv) ! Unfrozen volume for saturation (porosity*xunfroz)
1715	thksat(ji,jg,jv) = thks(ji,jg,jv) ** (1. - porosnet(ji, jg, jv)) &
1716	* THKICE ** (porosnet(ji, jg, jv) - xu(ji,jg,jv)) * THKW ** xu(ji,jg,jv)
1717	ELSE
1718	! this value will not be used since ake=0 for this case
1719	thksat(ji,jg,jv)=0
1720	END IF
1721	ENDDO
1722	END DO ! DO jg = 1,ngrnd
1723
1724	!! 2. Kerston Number (ake)
1725	DO jg = 1,ngrnd
1726	DO ji = 1,kjpindex
1727	IF ( (sh2o(ji,jg) + 0.0005) < smc(ji,jg) ) THEN
1728	! Frozen
1729	ake(ji,jg,jv) = satratio(ji,jg,jv)
1730	ELSE
1731	! Unfrozen
1732	IF ( satratio(ji,jg,jv) > 0.1 ) THEN
1733	IF ((jst < 4 .AND. soil_classif == 'usda') .OR. (jst == 1 .AND. soil_classif == 'zobler') ) THEN
1734	! Coarse
1735	ake(ji,jg,jv) = 0.7 * LOG10 (satratio(ji,jg,jv)) + 1.0
1736	ELSE
1737	! Fine
1738	ake(ji,jg,jv) = LOG10 (satratio(ji,jg,jv)) + 1.0
1739	ENDIF
1740	ELSEIF ( satratio(ji,jg,jv) > 0.05 .AND. satratio(ji,jg,jv) <= 0.1 ) THEN
1741	IF ((jst < 4 .AND. soil_classif == 'usda') .OR. (jst == 1 .AND. soil_classif == 'zobler') ) THEN
1742	! Coarse
1743	ake(ji,jg,jv) = 0.7 * LOG10 (satratio(ji,jg,jv)) + 1.0
1744	ELSE
1745	! Fine
1746	ake(ji,jg,jv) = 0.0
1747	ENDIF
1748	ELSE
1749	ake(ji,jg,jv) = 0.0 ! use k = kdry
1750	END IF
1751	END IF
1752	ENDDO
1753	END DO ! DO jg = 1,ngrnd
1754
1755	SELECTCASE (use_soilc_method)
1756	CASE (SOILC_METHOD_ARITHMETIC)
1757	DO jg = 1,ngrnd
1758	DO ji = 1,kjpindex
1759	thkdry(ji,jg,jv) = zx1(ji,jg,jv) * cond_dry_org + zx2(ji,jg,jv) * thkdry_min(ji)
1760	ENDDO
1761	ENDDO
1762	CASE (SOILC_METHOD_GEOMETRIC)
1763	DO jg = 1,ngrnd
1764	DO ji = 1,kjpindex
1765	! use geometric mean rather than arithmetic mean (Decharme et al 2016)
1766	thkdry(ji,jg,jv) =(cond_dry_org*zx1(ji,jg,jv)) (thkdry_min(ji)**zx2(ji,jg,jv))
1767	ENDDO
1768	ENDDO
1769	CASE DEFAULT
1770	CALL ipslerr_p(3,'thermosoil_cond_pft','Unsupported USE_SOILC_METHOD','','')
1771	ENDSELECT
1772
1773	!! 3. Thermal conductivity (cnd)
1774	DO jg = 1,ngrnd
1775	DO ji = 1,kjpindex
1776	cnd(ji,jg,jv) = ake(ji,jg,jv) * (thksat(ji,jg,jv) - thkdry(ji, jg, jv)) + thkdry(ji, jg, jv)
1777	ENDDO
1778	END DO
1779
1780	END DO !DO jv = 1,nvm
1781
1782
1783	END SUBROUTINE thermosoil_cond_pft
1784
1785	!! ================================================================================================================================
1786	!! SUBROUTINE : thermosoil_humlev
1787	!!
1788	!>\BRIEF Interpolates the diagnostic soil humidity profile shumdiag_perma(nslm, diagnostic axis) onto
1789	!! the thermal axis, which gives shum_ngrnd_perma(ngrnd, thermal axis).
1790	!!
1791	!! DESCRIPTION : Interpolate the volumetric soil moisture content from the node to the interface of the layer.
1792	!! The values for the deep layers in thermosoil where hydrology is not existing are constant.
1793	!! No interpolation is needed for the total soil moisture content and for the soil saturation degree.
1794	!! The depths of the diagnostic levels are diaglev(1:nslm), computed in slowproc.f90.
1795	!! Recall that when the 11-layer hydrology is used,
1796	!! shum_ngrnd_perma and shumdiag_perma are with reference to the moisture content (mc)
1797	!! at the wilting point mcw : shum_ngrnd_perma=(mc-mcw)/(mcs-mcw).
1798	!! with mcs the saturated soil moisture content.
1799	!!
1800	!! RECENT CHANGE(S) : None
1801	!!
1802	!! MAIN OUTPUT VARIABLE(S): mc_layt, mcl_layt, tmc_layt, shum_ngrnd_perma (soil humidity profile on the thermal axis)
1803	!!
1804	!! REFERENCE(S) : None
1805	!!
1806	!! FLOWCHART : None
1807	!! \n
1808	!_ ================================================================================================================================
1809	SUBROUTINE thermosoil_humlev(kjpindex, shumdiag_perma, mc_layh, mcl_layh, tmc_layh, &
1810	mc_layh_pft, mcl_layh_pft, tmc_layh_pft)
1811
1812	!! 0. Variables and parameter declaration
1813
1814	!! 0.1 Input variables
1815
1816	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (unitless)
1817	REAL(r_std),DIMENSION (kjpindex,nslm), INTENT (in) :: shumdiag_perma !! Relative soil humidity on the diagnostic axis.
1818	!! (0-1, unitless). Caveats : when "hydrol" (the 11-layers
1819	!! hydrology) is used, this humidity is calculated with
1820	!! respect to the wilting point :
1821	!! shumdiag_perma= (mc-mcw)/(mcs-mcw), with mc : moisture
1822	!! content; mcs : saturated soil moisture content; mcw:
1823	!! soil moisture content at the wilting point. when the 2-layers
1824	!! hydrology "hydrolc" is used, shumdiag_perma is just
1825	!! a diagnostic humidity index, with no real physical
1826	!! meaning.
1827	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]
1828	REAL(r_std),DIMENSION (kjpindex,nslm), INTENT (in) :: mcl_layh !! Volumetric soil moisture content for each layer in hydrol at nodes(liquid) [m/s]
1829	REAL(r_std),DIMENSION (kjpindex,nslm), INTENT (in) :: tmc_layh !! Total soil moisture content for each layer in hydrol(liquid+ice) [mm]
1830	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]
1831	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]
1832	REAL(r_std),DIMENSION (kjpindex,nslm,nvm), INTENT (in) :: tmc_layh_pft !! Total soil moisture content for each layer in hydrol(liquid+ice) [mm]
1833
1834	!! 0.2 Output variables
1835
1836	!! 0.3 Modified variables
1837
1838	!! 0.4 Local variables
1839
1840	INTEGER(i_std) :: ji, jd, jv
1841
1842	!_ ================================================================================================================================
1843
1844	IF (printlev >= 4) WRITE(numout,*) 'Start thermosoil_humlev'
1845
1846	! The values for the deep layers in thermosoil where hydrology is not existing are constant.
1847	! For exemple if thermosoil uses 8m, and hydrol uses 2m vertical discretization,
1848	! the values between 2m and 8m are constant.
1849	! The moisture computed in hydrol is at the nodes (except for the
1850	! top and bottom layer which are at interfaces)
1851	! A linear interpolation is applied to obtain the moisture values at
1852	! the interfaces (mc_layt), from the mc_layh at the nodes
1853
1854	DO ji=1,kjpindex
1855	DO jd = 1, nslm
1856	IF(jd == 1) THEN ! the moisture at the 1st interface mc_layh(1) is at the surface, no interpolation
1857	mc_layt(ji,jd) = mc_layh(ji,jd)
1858	mcl_layt(ji,jd) = mcl_layh(ji,jd)
1859	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
1860	mc_layt(ji, jd) = mc_layh(ji,jd-1)*(znt(jd)-zlt(jd-1))/(znt(jd)-0.0) + &
1861	mc_layh(ji, jd)*(zlt(jd-1)-0.0)/(znt(jd)-0.0)
1862	mcl_layt(ji, jd) = mcl_layh(ji,jd-1)*(znt(jd)-zlt(jd-1))/(znt(jd)-0.0) + &
1863	mcl_layh(ji, jd)*(zlt(jd-1)-0.0)/(znt(jd)-0.0)
1864	ELSEIF(jd == nslm) THEN ! the mc_layt at the nslm interface is interpolated using mc_layh(nslm) and mc_layh(nslm-1)
1865	mc_layt(ji, jd) = mc_layh(ji,jd-1)*(zlt(jd)-zlt(jd-1))/(zlt(jd)-znt(jd-1)) + &
1866	mc_layh(ji,jd)*(zlt(jd-1)-znt(jd-1))/(zlt(jd)-znt(jd-1))
1867	mcl_layt(ji, jd) = mcl_layh(ji,jd-1)*(zlt(jd)-zlt(jd-1))/(zlt(jd)-znt(jd-1)) + &
1868	mcl_layh(ji,jd)*(zlt(jd-1)-znt(jd-1))/(zlt(jd)-znt(jd-1))
1869	ELSE ! the mc_layt at the other interfaces are interpolated using mc_layh at adjacent nodes.
1870	mc_layt(ji, jd) = mc_layh(ji, jd-1)(1-dz5(jd-1)) + mc_layh(ji,jd)dz5(jd-1)
1871	mcl_layt(ji, jd) = mcl_layh(ji, jd-1)(1-dz5(jd-1)) + mcl_layh(ji,jd)dz5(jd-1)
1872	ENDIF
1873	DO jv = 1,nvm
1874	IF(jd == 1) THEN
1875	mc_layt_pft(ji,jd,jv) = MAX(mc_layh_pft(ji,jd,jv), min_sechiba)
1876	mcl_layt_pft(ji,jd,jv) = MAX(mcl_layh_pft(ji,jd,jv), min_sechiba)
1877	ELSEIF(jd == 2) THEN
1878	mc_layt_pft(ji,jd,jv) = MAX(mc_layh_pft(ji,jd-1,jv)*(znt(jd)-zlt(jd-1))/(znt(jd)-0.0) + &
1879	mc_layh_pft(ji,jd,jv)*(zlt(jd-1)-0.0)/(znt(jd)-0.0), min_sechiba)
1880	mcl_layt_pft(ji,jd,jv) = MAX(mcl_layh_pft(ji,jd-1,jv)*(znt(jd)-zlt(jd-1))/(znt(jd)-0.0) + &
1881	mcl_layh_pft(ji,jd,jv)*(zlt(jd-1)-0.0)/(znt(jd)-0.0), min_sechiba)
1882	ELSEIF(jd == nslm) THEN
1883	mc_layt_pft(ji,jd,jv) = MAX(mc_layh_pft(ji,jd-1,jv)*(zlt(jd)-zlt(jd-1))/(zlt(jd)-znt(jd-1)) + &
1884	mc_layh_pft(ji,jd,jv)*(zlt(jd-1)-znt(jd-1))/(zlt(jd)-znt(jd-1)),min_sechiba)
1885	mcl_layt_pft(ji,jd,jv) = MAX(mcl_layh_pft(ji,jd-1,jv)*(zlt(jd)-zlt(jd-1))/(zlt(jd)-znt(jd-1)) + &
1886	mcl_layh_pft(ji,jd,jv)*(zlt(jd-1)-znt(jd-1))/(zlt(jd)-znt(jd-1)),min_sechiba)
1887	ELSE
1888	mc_layt_pft(ji,jd,jv) = MAX(mc_layh_pft(ji,jd-1,jv)*(1-dz5(jd-1)) + &
1889	mc_layh_pft(ji,jd,jv)*dz5(jd-1), min_sechiba)
1890	mcl_layt_pft(ji,jd,jv) = MAX(mcl_layh_pft(ji,jd-1,jv)*(1-dz5(jd-1)) + &
1891	mcl_layh_pft(ji,jd,jv)*dz5(jd-1), min_sechiba)
1892	ENDIF
1893	ENDDO ! jv
1894	tmc_layt(ji,jd) = tmc_layh(ji,jd)
1895	tmc_layt_pft(ji,jd,:) = tmc_layh_pft(ji,jd,:)
1896	ENDDO !jd
1897
1898	! The deep layers in thermosoil where hydro is not existing
1899	DO jd = nslm+1, ngrnd
1900	mc_layt(ji,jd) = mc_layh(ji,nslm)
1901	mcl_layt(ji,jd) = mcl_layh(ji,nslm)
1902	tmc_layt(ji,jd) = tmc_layh(ji,nslm)/dlt(nslm) *dlt(jd)
1903	mc_layt_pft(ji,jd,:) = mc_layh_pft(ji,nslm,:)
1904	mcl_layt_pft(ji,jd,:) = mcl_layh_pft(ji,nslm,:)
1905	tmc_layt_pft(ji,jd,:) = tmc_layh_pft(ji,nslm,:)/dlt(nslm) *dlt(jd)
1906	ENDDO
1907	ENDDO
1908
1909
1910	! The values for the deep layers in thermosoil where hydro is not existing are constant.
1911	! For exemple if thermosoil uses 8m, and hydrol uses 2m vertical discretization,
1912	! the values between 2m and 8m are constant.
1913
1914	DO jd = 1, nslm
1915	shum_ngrnd_perma(:,jd) = shumdiag_perma(:,jd)
1916	END DO
1917	DO jd = nslm+1,ngrnd
1918	shum_ngrnd_perma(:,jd) = shumdiag_perma(:,nslm)
1919	END DO
1920
1921	IF (printlev >= 4) WRITE(numout,*) 'thermosoil_humlev done'
1922
1923	END SUBROUTINE thermosoil_humlev
1924
1925
1926
1927	!! ================================================================================================================================
1928	!! SUBROUTINE : thermosoil_energy_diag
1929	!!
1930	!>\BRIEF Calculate diagnostics
1931	!!
1932	!! DESCRIPTION : Calculate diagnostic variables coldcont_incr and coldcont_incr
1933	!!
1934	!! RECENT CHANGE(S) : None
1935	!!
1936	!! MAIN OUTPUT VARIABLE(S) :
1937	!!
1938	!! REFERENCE(S) : None
1939	!!
1940	!! FLOWCHART : None
1941	!! \n
1942	!_ ================================================================================================================================
1943
1944	SUBROUTINE thermosoil_energy_diag(kjpindex, temp_sol_new, soilcap)
1945
1946	!! 0. Variables and parameter declaration
1947
1948	!! 0.1 Input variables
1949
1950	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (unitless)
1951	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_sol_new !! Surface temperature at the present time-step, Ts
1952	!! @tex ($K$) @endtex
1953	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: soilcap !! Apparent surface heat capacity
1954	!! @tex ($J m^{-2} K^{-1}$) @endtex,
1955	!! see eq. A29 of F. Hourdin\'s PhD thesis.
1956
1957	!! 0.2 Output variables
1958
1959	!! 0.3 Modified variables
1960
1961	!! 0.4 Local variables
1962
1963	INTEGER(i_std) :: ji, jg
1964	!_ ================================================================================================================================
1965
1966	! Sum up the energy content of all layers in the soil.
1967	DO ji = 1, kjpindex
1968
1969	IF (pcapa_en(ji,1) .LE. sn_capa) THEN
1970
1971	! Verify the energy conservation in the surface layer
1972	coldcont_incr(ji) = soilcap(ji) * (temp_sol_new(ji) - temp_sol_beg(ji))
1973	surfheat_incr(ji) = zero
1974	ELSE
1975
1976	! Verify the energy conservation in the surface layer
1977	surfheat_incr(ji) = soilcap(ji) * (temp_sol_new(ji) - temp_sol_beg(ji))
1978	coldcont_incr(ji) = zero
1979	ENDIF
1980	ENDDO
1981
1982	! Save temp_sol_new to be used at next timestep
1983	temp_sol_beg(:) = temp_sol_new(:)
1984
1985	END SUBROUTINE thermosoil_energy_diag
1986
1987
1988
1989	!! ================================================================================================================================
1990	!! SUBROUTINE : thermosoil_readjust
1991	!!
1992	!>\BRIEF
1993	!!
1994	!! DESCRIPTION : Energy conservation : Correction to make sure that the same latent heat is released and
1995	!! consumed during freezing and thawing
1996	!!
1997	!! RECENT CHANGE(S) : None
1998	!!
1999	!! MAIN OUTPUT VARIABLE(S): ptn (soil temperature profile on the thermal axis per pft) and ptn_pftmean
2000	!!
2001	!! REFERENCE(S) :
2002	!!
2003	!! FLOWCHART : None
2004	!! \n
2005	!_ ================================================================================================================================
2006
2007	SUBROUTINE thermosoil_readjust(kjpindex, ptn, ptn_pftmean, veget_max)
2008
2009	!! 0. Variables and parameter declaration
2010
2011	!! 0.1 Input variables
2012	INTEGER(i_std), INTENT(in) :: kjpindex
2013	REAL(r_std), DIMENSION (kjpindex,nvm), INTENT(in) :: veget_max !! Fraction of vegetation type
2014
2015	!! 0.2 Modified variables
2016	REAL(r_std),DIMENSION(kjpindex,ngrnd,nvm),INTENT(inout) :: ptn
2017	REAL(r_std),DIMENSION(kjpindex,ngrnd),INTENT(inout) :: ptn_pftmean
2018
2019	!! 0.3 Local variables
2020	INTEGER(i_std) :: ji, jg, m, jv
2021	INTEGER(i_std) :: lev3m !! Closest interface level to 3m
2022	REAL(r_std) :: ptn_tmp
2023	REAL(r_std), DIMENSION (kjpindex,ngrnd) :: ptn_beg_pftmean
2024
2025	! The energy is spread over the layers down to approximatly 3m
2026	! Find the closest level to 3m. It can be below or above 3m.
2027	lev3m=MINLOC(ABS(zlt(:)-3.0),dim=1)
2028	IF (printlev >= 3) WRITE(numout,*) 'In thermosoil_adjust: lev3m=',lev3m, ' zlt(lev3m)=', zlt(lev3m)
2029
2030	! Calculate ptn_beg for mean PFT
2031	IF (ok_soil_carbon_discretization) THEN
2032	ptn_beg_pftmean(:,:) = zero
2033	DO m=1,nvm
2034	DO jg = 1, ngrnd
2035	ptn_beg_pftmean(:,jg) = ptn_beg_pftmean(:,jg) + ptn_beg(:,jg,m) * veget_max(:,m)
2036	END DO
2037	END DO
2038	ELSE
2039	! ptn is not depending on pft
2040	ptn_beg_pftmean(:,:) = ptn_beg(:,:,1)
2041	END IF
2042
2043	DO jg=1, ngrnd
2044	DO ji=1, kjpindex
2045	! All soil latent energy is put into e_soil_lat(ji)
2046	! because the variable soil layers make it difficult to keep track of all
2047	! layers in this version
2048	! NOTE: pcapa has unit J/K/m3 and pcappa_supp has J/K/m2
2049	! NOTE: Right now this is not effective as ptn_beg_pftmean equal ptn_pftmean
2050	e_soil_lat(ji)=e_soil_lat(ji)+pcappa_supp(ji,jg)*(ptn_pftmean(ji,jg)-ptn_beg_pftmean(ji,jg))
2051	END DO
2052	END DO
2053
2054	DO ji=1, kjpindex
2055	IF (e_soil_lat(ji).GT.min_sechiba.AND.MINVAL(ptn_pftmean(ji,:)).GT.ZeroCelsius+fr_dT/2.) THEN
2056	! The soil is thawed: we spread the excess of energy over the uppermost lev3m levels
2057	! Here we increase the temperatures
2058	DO jg=1, lev3m
2059	ptn_tmp=ptn_pftmean(ji,jg)
2060
2061	ptn_pftmean(ji,jg)=ptn_pftmean(ji,jg)+MIN(e_soil_lat(ji)/pcapa(ji,jg)/zlt(lev3m), 0.5)
2062	e_soil_lat(ji)=e_soil_lat(ji)-(ptn_pftmean(ji,jg)-ptn_tmp)pcapa(ji,jg)dlt(jg)
2063	ENDDO
2064	ELSE IF (e_soil_lat(ji).LT.-min_sechiba.AND.MINVAL(ptn_pftmean(ji,:)).GT.ZeroCelsius+fr_dT/2.) THEN
2065	! The soil is thawed
2066	! Here we decrease the temperatures
2067	DO jg=1, lev3m
2068	ptn_tmp=ptn_pftmean(ji,jg)
2069	ptn_pftmean(ji,jg)=MAX(ZeroCelsius+fr_dT/2., ptn_tmp+e_soil_lat(ji)/pcapa(ji,jg)/zlt(lev3m))
2070	e_soil_lat(ji)=e_soil_lat(ji)+(ptn_tmp-ptn_pftmean(ji,jg))pcapa(ji,jg)dlt(jg)
2071	END DO
2072	END IF
2073	END DO
2074
2075
2076	!! Calculate ptn per pft
2077	! Here ptn is the same at all pfts when ok_soil_carbon_discretization is activated or not.
2078	! ptn is modified in thermosoil_add_heat_zimov if ok_zimov=TRUE.
2079	DO jv=1,nvm
2080	ptn(:,:,jv) = ptn_pftmean(:,:)
2081	END DO
2082
2083
2084	END SUBROUTINE thermosoil_readjust
2085
2086	!-------------------------------------------------------------------
2087
2088
2089
2090	!! ================================================================================================================================
2091	!! SUBROUTINE : thermosoil_getdiff
2092	!!
2093	!>\BRIEF Computes soil and snow heat capacity and conductivity
2094	!!
2095	!! DESCRIPTION : Computation of the soil thermal properties; snow properties are also accounted for
2096	!!
2097	!! RECENT CHANGE(S) : None
2098	!!
2099	!! MAIN OUTPUT VARIABLE(S):
2100	!!
2101	!! REFERENCE(S) :
2102	!!
2103	!! FLOWCHART : None
2104	!! \n
2105	!_ ================================================================================================================================
2106
2107	SUBROUTINE thermosoil_getdiff( kjpindex, snow, ptn, njsc, snowrho, snowtemp, pb )
2108
2109	!! 0. Variables and parameter declaration
2110
2111	!! 0.1 Input variables
2112	INTEGER(i_std),INTENT(in) :: kjpindex
2113	REAL(r_std),DIMENSION(kjpindex),INTENT (in) :: snow !! Snow mass
2114	INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: njsc !! Index of the dominant soil textural class in the grid cell (1-nscm, unitless)
2115	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(in) :: snowrho !! Snow density (Kg/m^3)
2116	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(in) :: snowtemp !! Snow temperature (K)
2117	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: pb !! Surface presure (hPa)
2118	REAL(r_std),DIMENSION(kjpindex,ngrnd),INTENT(in) :: ptn !! Soil temperature profile
2119
2120	!! 0.3 Local variables
2121	REAL :: xx !! Unfrozen fraction of the soil
2122	REAL(r_std), DIMENSION(kjpindex) :: snow_h
2123	REAL(r_std), DIMENSION(kjpindex,ngrnd) :: zx1, zx2
2124	INTEGER :: ji,jg
2125	INTEGER :: jst
2126	REAL(r_std), DIMENSION (kjpindex,ngrnd) :: pcapa_tmp !! soil heat capacity (J/m3/K)
2127	REAL(r_std), DIMENSION (kjpindex,ngrnd) :: pcapa_spec !! SPECIFIC soil heat capacity (J/kg/K)
2128	REAL(r_std) :: rho_tot !! Soil density (kg/m3)
2129
2130	pcapa_tmp(:,:) = 0.0
2131
2132	!! Computes soil heat capacity and conductivity
2133	DO ji = 1,kjpindex
2134
2135	! Since explicitsnow module is implemented set zx1=0 and zx2=1
2136	zx1(ji,:) = 0.
2137	zx2(ji,:) = 1.
2138
2139	DO jg = 1, ngrnd
2140	jst = njsc(ji)
2141	pcapa_tmp(ji, jg) = so_capa_dry_ns(jst) * (1-mcs(jst)) + water_capa * tmc_layt(ji,jg)/mille/dlt(jg)
2142	!
2143	! 2. Calculate volumetric heat capacity with allowance for permafrost
2144	! 2.1. soil heat capacity depending on temperature and humidity
2145	! For SP6MIP we also diagnose a specific heat capacity (pcapa_spec),
2146	! which requires to calculate the total density of the soil (rho_tot), always >> 0
2147
2148	IF (ptn(ji,jg) .LT. ZeroCelsius-fr_dT/2.) THEN
2149	! frozen soil
2150	profil_froz(ji,jg) = 1.
2151	pcappa_supp(ji,jg)= 0.
2152	pcapa(ji, jg) = so_capa_dry_ns(jst) * (1-mcs(jst)) + so_capa_ice(jst) * tmc_layt(ji,jg) / mille / dlt(jg)
2153	rho_tot = rho_soil * (1-mcs(jst)) + rho_ice * tmc_layt(ji,jg) / mille / dlt(jg)
2154	pcapa_spec(ji, jg) = pcapa(ji, jg) / rho_tot
2155
2156	ELSEIF (ptn(ji,jg) .GT. ZeroCelsius+fr_dT/2.) THEN
2157	! unfrozen soil
2158	pcapa(ji, jg) = pcapa_tmp(ji, jg)
2159	profil_froz(ji,jg) = 0.
2160	pcappa_supp(ji,jg)= 0.
2161	rho_tot = rho_soil * (1-mcs(jst)) + rho_water * tmc_layt(ji,jg)/mille/dlt(jg)
2162	pcapa_spec(ji, jg) = pcapa(ji, jg) / rho_tot
2163	ELSE
2164	! xx is the unfrozen fraction of soil water
2165	xx = (ptn(ji,jg)-(ZeroCelsius-fr_dT/2.)) / fr_dT
2166	profil_froz(ji,jg) = (1. - xx)
2167
2168	IF (ok_freeze_thaw_latent_heat) THEN
2169	pcapa(ji, jg) = so_capa_dry_ns(jst) * (1-mcs(jst)) + &
2170	water_capa * tmc_layt(ji,jg)/mille / dlt(jg) * xx + &
2171	so_capa_ice(jst) * tmc_layt(ji,jg) / mille/dlt(jg) * (1.-xx) + &
2172	shum_ngrnd_perma(ji,jg)mcs(jst)lhf*rho_water/fr_dT
2173	ELSE
2174	pcapa(ji, jg) = so_capa_dry_ns(jst) * (1-mcs(jst)) + &
2175	water_capa * tmc_layt(ji,jg)/mille / dlt(jg) * xx + &
2176	so_capa_ice(jst) * tmc_layt(ji,jg) / mille/dlt(jg) * (1.-xx)
2177	ENDIF
2178
2179	rho_tot = rho_soil* (1-mcs(jst)) + &
2180	rho_water * tmc_layt(ji,jg)/mille / dlt(jg) * xx + &
2181	rho_ice * tmc_layt(ji,jg) / mille/dlt(jg) * (1.-xx)
2182	pcapa_spec(ji, jg) = pcapa(ji, jg) / rho_tot
2183
2184	! NOTE: pcappa_supp is transformed from volumetric heat capacity into surfacic heat capacity (J/K/m2)
2185	pcappa_supp(ji,jg)= shum_ngrnd_perma(ji,jg)mcs(jst)lhfrho_water/fr_dTzx2(ji,jg)*dlt(jg)
2186
2187	ENDIF
2188	!BG:is this 2.2 and 2.3 necessary since zx1=0?
2189	!
2190	! 2.2. Take into account the snow and soil fractions in the layer
2191	!
2192	pcapa(ji,jg) = zx1(ji,jg) * sn_capa + zx2(ji,jg) * pcapa(ji,jg)
2193
2194	!
2195	! 2.3. Calculate the heat capacity for energy conservation check
2196	IF ( zx1(ji,jg).GT.0. ) THEN
2197	pcapa_en(ji,jg) = sn_capa
2198	ELSE
2199	pcapa_en(ji,jg) = pcapa(ji,jg)
2200	ENDIF
2201
2202	END DO
2203	ENDDO
2204
2205	! Output the specific heat capcaity for SP-MIP
2206	CALL xios_orchidee_send_field("pcapa_spec",pcapa_spec)
2207
2208	!
2209	! 3. Calculate the heat conductivity with allowance for permafrost
2210	!
2211	IF (ok_freeze_thaw_latent_heat) THEN
2212	CALL thermosoil_cond (kjpindex, njsc, mc_layt, QZ, mcl_layt*(1-profil_froz), pkappa)
2213	ELSE
2214	CALL thermosoil_cond (kjpindex, njsc, mc_layt, QZ, mcl_layt, pkappa)
2215	ENDIF
2216
2217	!! Computes snow heat capacity and conductivity
2218	DO ji = 1,kjpindex
2219	pcapa_snow(ji,:) = snowrho(ji,:) * xci
2220	pkappa_snow(ji,:) = (ZSNOWTHRMCOND1 + ZSNOWTHRMCOND2snowrho(ji,:)snowrho(ji,:)) + &
2221	MAX(0.0,(ZSNOWTHRMCOND_AVAP+(ZSNOWTHRMCOND_BVAP/(snowtemp(ji,:)+ &
2222	ZSNOWTHRMCOND_CVAP)))(XP00/(pb(ji)100.)))
2223	END DO
2224
2225
2226	END SUBROUTINE thermosoil_getdiff
2227
2228
2229
2230	!! ================================================================================================================================
2231	!! SUBROUTINE : thermosoil_getdiff_pft
2232	!!
2233	!>\BRIEF Computes soil and snow heat capacity and conductivity
2234	!!
2235	!! DESCRIPTION : Computation of the soil thermal properties; snow properties are also accounted for
2236	!!
2237	!! RECENT CHANGE(S) : None
2238	!!
2239	!! MAIN OUTPUT VARIABLE(S): Following module variables are calculted in this subroutine:
2240	!! profil_froz_pft, pcapa_per_pft, pkappa_per_pft, pcapa_en_per_pft, tmc_layt_pft
2241	!!
2242	!! REFERENCE(S) :
2243	!!
2244	!! FLOWCHART : None
2245	!! \n
2246	!_ ================================================================================================================================
2247
2248	SUBROUTINE thermosoil_getdiff_pft( kjpindex, ptn, njsc, shum_ngrnd_perma, &
2249	depth_organic_soil, som_total, snowrho, snowtemp, pb, veget_max)
2250
2251	!! 0. Variables and parameter declaration
2252
2253	!! 0.1 Input variables
2254	INTEGER(i_std),INTENT(in) :: kjpindex
2255	REAL(r_std),DIMENSION(kjpindex,ngrnd),INTENT(in) :: shum_ngrnd_perma !! Water saturation degree on the soil depth axes define by the thermic (0-1, dimensionless)
2256	REAL(r_std), DIMENSION(kjpindex), INTENT (in) :: depth_organic_soil !! Depth at which there is still organic matter (m)
2257	REAL(r_std), DIMENSION(kjpindex,ngrnd,nvm,nelements), INTENT (in) :: som_total !! total soil carbon for use in thermal calcs (g/m**3)
2258	INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: njsc !! Index of the dominant soil textural class in the grid cell (1-nscm, unitless)
2259	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(in) :: snowrho !! Snow density (Kg/m^3)
2260	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(in) :: snowtemp !! Snow temperature (K)
2261	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: pb !! Surface presure (hPa)
2262	REAL(r_std),DIMENSION(kjpindex,ngrnd,nvm),INTENT(in) :: ptn !! Soil temperature profile
2263	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(in) :: veget_max !! Fraction of vegetation type
2264	!! 0.2 Modified variables
2265
2266	!! 0.3 Output variables
2267
2268	!! 0.3 Local variables
2269	REAL(r_std) :: xx !! Unfrozen fraction of the soil
2270	REAL(r_std) :: p
2271	REAL(r_std), DIMENSION(kjpindex,ngrnd,nvm) :: so_capa_dry_net
2272	REAL(r_std) :: cond_solid_net
2273	REAL(r_std) :: so_cond_dry_net
2274	INTEGER(i_std) :: ji,jg,jv
2275	INTEGER(i_std) :: jst
2276	REAL(r_std), DIMENSION(kjpindex,ngrnd,nvm) :: poros_net
2277	REAL(r_std), DIMENSION(kjpindex,ngrnd,nvm) :: zx1, zx2 !! BE CAREFUL here zx1 and zx2 split between the organic and mineral soils not
2278	!! not snow vs. liquid water
2279	REAL(r_std), DIMENSION(kjpindex,ngrnd) :: profil_froz_mean, tmp
2280	REAL(r_std), DIMENSION(kjpindex,ngrnd,nvm) :: pcappa_supp_pft
2281	REAL(r_std), DIMENSION(kjpindex,ngrnd,nvm) :: pcapa_pft_tmp !! soil heat capacity (J/m3/K)
2282	REAL(r_std), DIMENSION(kjpindex,ngrnd,nvm) :: pcapa_spec_pft !! SPECIFIC soil heat capacity (J/kg/K)
2283	REAL(r_std) :: rho_tot !! Soil density (kg/m3)
2284	REAL(r_std), DIMENSION(kjpindex,ngrnd) :: pcapa_spec !! SPECIFIC soil heat capacity (J/kg/K)
2285	! Organic and anorgaic layer fraction
2286	!
2287	! Default: organic layer not taken into account
2288	zx1(:,:,:) = 0.0
2289	zx2(:,:,:) = 0.0
2290	poros_net(:,:,:) = 0.0
2291	pcapa_spec_pft(:,:,:) = 0.0
2292	pcapa_pft_tmp(:,:,:) = 0.0
2293	pcapa_spec(:,:) = 0.0
2294	!
2295	IF ( use_soilc_tempdiff ) THEN
2296	!
2297	IF (use_refsoc) THEN
2298	DO jv = 1,nvm
2299	DO jg = 1, ngrnd
2300	DO ji = 1,kjpindex
2301	zx1(ji,jg,jv) = refsoc(ji,jg)/soilc_max !after lawrence and slater
2302	ENDDO
2303	ENDDO
2304	ENDDO
2305	ELSE
2306	DO jv = 1,nvm
2307	DO jg = 1, ngrnd
2308	DO ji = 1,kjpindex
2309	zx1(ji,jg,jv) = som_total(ji,jg,jv,icarbon)/soilc_max !after lawrence and slater
2310	ENDDO
2311	ENDDO
2312	ENDDO
2313	ENDIF
2314	!
2315	WHERE (zx1 > 1) zx1 = 1
2316	!
2317	ELSE
2318	!
2319	! level 1
2320	!
2321	DO jv = 1,nvm
2322	DO ji = 1,kjpindex
2323	IF ( depth_organic_soil(ji) .GT. zlt(1) ) THEN
2324	!! the 1st level is in the organic => the 1st layer is entirely organic
2325	zx1(ji,1,jv) = 1. !!zx1 being the fraction of each level that is organic, zx2 is the remainder
2326	ELSE IF ( depth_organic_soil(ji) .GT. zero ) THEN
2327	!! the 1st level is beyond the organic and the organic is present
2328	zx1(ji,1,jv) = depth_organic_soil(ji) / zlt(1)
2329	ELSE
2330	! there is no organic at all
2331	zx1(ji,1,jv) = 0.
2332	ENDIF
2333	ENDDO
2334	ENDDO
2335	!
2336	! other levels
2337	!
2338	DO jg = 2, ngrnd !- 2
2339	DO ji = 1,kjpindex
2340	IF ( depth_organic_soil(ji) .GT. zlt(jg) ) THEN
2341	! the current level is in the organic => the current layer is
2342	! entirely organic
2343	zx1(ji,jg,1) = 1.
2344	ELSE IF ( depth_organic_soil(ji) .GT. zlt(jg-1) ) THEN
2345	! the current layer is partially organic
2346	zx1(ji,jg,1) = (depth_organic_soil(ji) - zlt(jg-1)) / (zlt(jg) - zlt(jg-1))
2347	ELSE
2348	! both levels are out of organic => the current layer is entirely
2349	! mineral soil
2350	zx1(ji,jg,1) = 0.
2351	ENDIF
2352	ENDDO
2353	ENDDO
2354	DO jv = 2, nvm
2355	zx1(:,:,jv) = zx1(:,:,1)
2356	ENDDO
2357	ENDIF ! ( use_soilc_tempdiff )
2358	!
2359	zx2(:,:,:) = 1.-zx1(:,:,:)
2360
2361	#ifdef STRICT_CHECK
2362	IF (ANY(tmc_layt_pft < min_sechiba)) CALL ipslerr_p(3, "thermosoil_getdiff_pft", "tmc_layt_pft has negative values", "", "") ! prec issues
2363	#endif
2364
2365	DO jv = 1,nvm
2366	DO jg = 1, ngrnd
2367	DO ji = 1,kjpindex
2368	jst = njsc(ji)
2369	!
2370	! 1. Calculate dry heat capacity and conductivity, taking
2371	! into account the organic and mineral fractions in the layer
2372	!
2373	! Former MICT version
2374	!Here we take into account the new dependance of the soil heat capacity from the soil type.
2375	so_capa_dry_net(ji,jg,jv) = zx1(ji,jg,jv) * SO_CAPA_DRY_ORG + zx2(ji,jg,jv) * so_capa_dry_ns(jst)
2376
2377	!cond_solid_net = un / ( zx1(ji,jg,jv) / cond_solid_org + zx2(ji,jg,jv) / cond_solid ) ! TO DELETE
2378	!Here we take into account the new dependance of the porosity from the soil type.
2379	poros_net(ji,jg,jv) = zx1(ji,jg,jv) * poros_org + zx2(ji,jg,jv) * mcs(jst)
2380	!
2381	!so_cond_dry_net = un / ( zx1(ji,jg,jv) / cond_dry_org + zx2(ji,jg,jv) / so_cond_dry ) ! TO DELETE
2382	!
2383	! 2. Calculate heat capacity with allowance for permafrost
2384	ENDDO
2385	ENDDO
2386	ENDDO
2387
2388	DO jv = 1,nvm
2389	DO jg = 1, ngrnd
2390	DO ji = 1,kjpindex
2391	jst = njsc(ji)
2392	pcapa_pft_tmp(ji, jg,jv) = so_capa_dry_net(ji,jg,jv) * (1-mcs(jst)) + water_capa * tmc_layt(ji,jg)/mille/dlt(jg)
2393	! 2.1. soil heat capacity depending on temperature and humidity
2394	IF (ptn(ji,jg,jv) .LT. ZeroCelsius-fr_dT/2.) THEN
2395	! frozen soil
2396	profil_froz_pft(ji,jg,jv) = 1.
2397	pcappa_supp_pft(ji,jg, jv)= 0.
2398	pcapa_per_pft(ji,jg,jv) = so_capa_dry_net(ji,jg,jv)* (1-mcs(jst)) + &
2399	so_capa_ice(jst) * tmc_layt(ji,jg) / mille / dlt(jg)
2400	rho_tot = rho_soil * (1-poros_net(ji,jg,jv)) + rho_ice * tmc_layt(ji,jg) / mille / dlt(jg)
2401	pcapa_spec_pft(ji, jg, jv) = pcapa_per_pft(ji, jg, jv) / rho_tot
2402	ELSEIF (ptn(ji,jg,jv) .GT. ZeroCelsius+fr_dT/2.) THEN
2403	! unfrozen soil
2404	profil_froz_pft(ji,jg,jv) = 0.
2405	pcappa_supp_pft(ji,jg,jv)= 0.
2406	pcapa_per_pft(ji,jg,jv) = pcapa_pft_tmp(ji, jg,jv)
2407	rho_tot = rho_soil * (1-poros_net(ji,jg,jv)) + rho_water * tmc_layt(ji,jg)/mille/dlt(jg)
2408	pcapa_spec_pft(ji, jg, jv) = pcapa_per_pft(ji, jg, jv) / rho_tot
2409	ELSE
2410
2411	pcappa_supp_pft(ji,jg,jv)= shum_ngrnd_perma(ji,jg)lhfrho_water/fr_dT
2412	IF (jg .GT. nslm) pcappa_supp_pft(ji,jg,jv)= 0.
2413
2414	! x is the unfrozen fraction of soil water
2415	xx = (ptn(ji,jg,jv)-(ZeroCelsius-fr_dT/2.)) / fr_dT
2416	profil_froz_pft(ji,jg,jv) = (1. - xx)
2417
2418	IF (ok_freeze_thaw_latent_heat) THEN
2419	pcapa_per_pft(ji,jg,jv) = so_capa_dry_net(ji,jg,jv)*(1-mcs(jst)) + &
2420	water_capa * tmc_layt(ji,jg)/mille / dlt(jg) * xx + &
2421	so_capa_ice(jst) * tmc_layt(ji,jg) / mille/dlt(jg) * (1.-xx) + &
2422	shum_ngrnd_perma(ji,jg)mcs(jst)lhf*rho_water/fr_dT
2423	ELSE
2424	pcapa_per_pft(ji,jg,jv) = so_capa_dry_net(ji,jg,jv)*(1-mcs(jst)) + &
2425	water_capa * tmc_layt(ji,jg)/mille / dlt(jg) * xx + &
2426	so_capa_ice(jst) * tmc_layt(ji,jg) / mille/dlt(jg) * (1.-xx)
2427	ENDIF
2428	rho_tot = rho_soil* (1-poros_net(ji,jg,jv)) + &
2429	rho_water * tmc_layt(ji,jg)/mille / dlt(jg) * xx + &
2430	rho_ice * tmc_layt(ji,jg) / mille/dlt(jg) * (1.-xx)
2431	pcapa_spec_pft(ji, jg,jv) = pcapa_per_pft(ji, jg,jv) / rho_tot
2432
2433	! NOTE: pcappa_supp is transformed from volumetric heat capacity into surfacic heat capacity (J/K/m2)
2434	pcappa_supp_pft(ji,jg,jv)= shum_ngrnd_perma(ji,jg)mcs(jst)lhfrho_water/fr_dTzx2(ji,jg,jv)*dlt(jg)
2435	ENDIF
2436
2437	ENDDO
2438	ENDDO
2439	ENDDO
2440
2441	! Output the specific heat capcaity for SP-MIP
2442	DO jv = 1,nvm
2443	DO jg = 1, ngrnd
2444	DO ji = 1,kjpindex
2445	pcapa_spec(ji,jg)=pcapa_spec(ji,jg) + pcapa_spec_pft(ji,jg,jv)*veget_max(ji,jv)
2446	ENDDO
2447	ENDDO
2448	ENDDO
2449
2450	CALL xios_orchidee_send_field("pcapa_spec",pcapa_spec)
2451	!
2452	! 3. Calculate the heat conductivity with allowance for permafrost
2453	! Note: mc_layt has no PFT dimention,so we calculate here profil_froz_mean. Actually, profil_froz_pft along the PFT dimention currently has no difference for each PFT.
2454	IF ( ANY(MAXVAL(profil_froz_pft,DIM=3)>MINVAL(profil_froz_pft,DIM=3)) ) THEN
2455	CALL ipslerr_p(3,'thermosoil_getdiff_pft','profil_froz_mean wrong','','')
2456	ENDIF
2457	profil_froz_mean=MINVAL(profil_froz_pft,DIM=3)
2458	tmp = mc_layt*(1.- profil_froz_mean)
2459	CALL thermosoil_cond_pft (kjpindex, njsc, mc_layt, QZ, tmp, zx1,zx2, poros_net,pkappa_per_pft)
2460
2461	!! Computes snow heat capacity and conductivity
2462	DO ji = 1,kjpindex
2463	pcapa_snow(ji,:) = snowrho(ji,:) * xci
2464
2465	SELECTCASE (snow_cond_method)
2466	CASE (SNOW_COND_METHOD_DEFAULT)
2467	pkappa_snow(ji,:) = (ZSNOWTHRMCOND1 + ZSNOWTHRMCOND2snowrho(ji,:)snowrho(ji,:)) + &
2468	MAX(0.0,(ZSNOWTHRMCOND_AVAP+(ZSNOWTHRMCOND_BVAP/(snowtemp(ji,:)+ &
2469	ZSNOWTHRMCOND_CVAP)))(XP00/(pb(ji)100.)))
2470	CASE (SNOW_COND_METHOD_DECHARME16)
2471	pkappa_snow(ji,:) = 2.2((snowrho(ji,:)/1000.)*2.0) + &
2472	MAX(0.0,(ZSNOWTHRMCOND_AVAP+(ZSNOWTHRMCOND_BVAP/(snowtemp(ji,:)+ &
2473	ZSNOWTHRMCOND_CVAP)))(XP00/(pb(ji)100.)))
2474	CASE DEFAULT
2475	CALL ipslerr_p(3,'thermosoil_getdiff_pft','Unsupported SNOW_COND_METHOD', &
2476	'Currently supported methods are ','default(1) or ducharme16(2)')
2477	ENDSELECT
2478
2479	END DO
2480	END SUBROUTINE thermosoil_getdiff_pft
2481
2482
2483	!! ================================================================================================================================
2484	!! SUBROUTINE : thermosoil_getdiff_old_thermix_without_snow
2485	!!
2486	!>\BRIEF Computes soil and snow heat capacity and conductivity
2487	!!
2488	!! DESCRIPTION : Calculations of soil and snow thermal properties without effect of freezing.
2489	!!
2490	!!
2491	!! RECENT CHANGE(S) : None
2492	!!
2493	!! MAIN OUTPUT VARIABLE(S):
2494	!!
2495	!! REFERENCE(S) :
2496	!!
2497	!! FLOWCHART : None
2498	!! \n
2499	!_ ================================================================================================================================
2500
2501	SUBROUTINE thermosoil_getdiff_old_thermix_without_snow( kjpindex, njsc, snowrho, snowtemp, pb )
2502
2503	!! 0. Variables and parameter declaration
2504
2505	!! 0.1 Input variables
2506	INTEGER(i_std), INTENT(in) :: kjpindex
2507	INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: njsc !! Index of the dominant soil textural class in the grid cell (1-nscm, unitless)
2508	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(in) :: snowrho !! Snow density (Kg/m^3)
2509	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(in) :: snowtemp !! Snow temperature (K)
2510	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: pb !! Surface presure (hPa)
2511
2512
2513	!! 0.1 Local variables
2514	INTEGER(i_std) :: ji,jg, jst !! Index
2515	REAL(r_std), DIMENSION (kjpindex,ngrnd) :: pcapa_tmp !! Soil heat capacity (J/m3/K)
2516
2517	!! Computes soil heat capacity and conductivity
2518	DO jg = 1,ngrnd
2519	DO ji = 1,kjpindex
2520	jst = njsc(ji)
2521	pcapa_tmp(ji, jg) = so_capa_dry_ns(jst) * (1-mcs(jst)) + water_capa * tmc_layt(ji,jg)/mille/dlt(jg)
2522	pcapa(ji,jg) = pcapa_tmp(ji, jg)
2523	pcapa_en(ji,jg) = pcapa_tmp(ji, jg)
2524	ENDDO
2525	ENDDO
2526
2527	CALL thermosoil_cond (kjpindex, njsc, mc_layt, QZ, mcl_layt, pkappa)
2528
2529	IF (brk_flag == 1) THEN
2530	! Bedrock flag is activated
2531	DO jg = ngrnd-1,ngrnd
2532	DO ji = 1,kjpindex
2533	pcapa(ji,jg) = brk_capa
2534	pcapa_en(ji,jg) = brk_capa
2535	pkappa(ji,jg) = brk_cond
2536	ENDDO
2537	ENDDO
2538	ENDIF
2539
2540	!! Computes snow heat capacity and conductivity
2541	DO ji = 1,kjpindex
2542	pcapa_snow(ji,:) = snowrho(ji,:) * xci
2543	pkappa_snow(ji,:) = (ZSNOWTHRMCOND1 + ZSNOWTHRMCOND2snowrho(ji,:)snowrho(ji,:)) + &
2544	MAX(0.0,(ZSNOWTHRMCOND_AVAP+(ZSNOWTHRMCOND_BVAP/(snowtemp(ji,:)+ &
2545	ZSNOWTHRMCOND_CVAP)))(XP00/(pb(ji)100.)))
2546	END DO
2547
2548	END SUBROUTINE thermosoil_getdiff_old_thermix_without_snow
2549
2550
2551	!! ================================================================================================================================
2552	!! SUBROUTINE : thermosoil_read_reftempfile
2553	!!
2554	!>\BRIEF
2555	!!
2556	!! DESCRIPTION : Read file with longterm soil temperature
2557	!!
2558	!!
2559	!! RECENT CHANGE(S) : None
2560	!!
2561	!! MAIN OUTPUT VARIABLE(S): reftemp : Reference temerature
2562	!!
2563	!! REFERENCE(S) :
2564	!!
2565	!! FLOWCHART : None
2566	!! \n
2567	!_ ================================================================================================================================
2568	SUBROUTINE thermosoil_read_reftempfile(kjpindex,lalo,reftemp)
2569
2570	USE interpweight
2571
2572	IMPLICIT NONE
2573
2574	!! 0. Variables and parameter declaration
2575
2576	!! 0.1 Input variables
2577	INTEGER(i_std), INTENT(in) :: kjpindex
2578	REAL(r_std), DIMENSION(kjpindex,2), INTENT(in) :: lalo
2579
2580	!! 0.2 Output variables
2581	REAL(r_std), DIMENSION(kjpindex, ngrnd), INTENT(out) :: reftemp
2582
2583	!! 0.3 Local variables
2584	INTEGER(i_std) :: ib
2585	CHARACTER(LEN=80) :: filename
2586	REAL(r_std),DIMENSION(kjpindex) :: reftemp_file !! Horizontal temperature field interpolated from file [C]
2587	INTEGER(i_std),DIMENSION(kjpindex,8) :: neighbours
2588	REAL(r_std) :: vmin, vmax !! min/max values to use for the
2589	!! renormalization
2590	REAL(r_std), DIMENSION(kjpindex) :: areftemp !! Availability of data for the interpolation
2591	CHARACTER(LEN=80) :: variablename !! Variable to interpolate
2592	!! the file
2593	CHARACTER(LEN=80) :: lonname, latname !! lon, lat names in input file
2594	REAL(r_std), DIMENSION(:), ALLOCATABLE :: variabletypevals !! Values for all the types of the variable
2595	!! (variabletypevals(1) = -un, not used)
2596	CHARACTER(LEN=50) :: fractype !! method of calculation of fraction
2597	!! 'XYKindTime': Input values are kinds
2598	!! of something with a temporal
2599	!! evolution on the dx*dy matrix'
2600	LOGICAL :: nonegative !! whether negative values should be removed
2601	CHARACTER(LEN=50) :: maskingtype !! Type of masking
2602	!! 'nomask': no-mask is applied
2603	!! 'mbelow': take values below maskvals(1)
2604	!! 'mabove': take values above maskvals(1)
2605	!! 'msumrange': take values within 2 ranges;
2606	!! maskvals(2) <= SUM(vals(k)) <= maskvals(1)
2607	!! maskvals(1) < SUM(vals(k)) <= maskvals(3)
2608	!! (normalized by maskvals(3))
2609	!! 'var': mask values are taken from a
2610	!! variable inside the file (>0)
2611	REAL(r_std), DIMENSION(3) :: maskvals !! values to use to mask (according to
2612	!! `maskingtype')
2613	CHARACTER(LEN=250) :: namemaskvar !! name of the variable to use to mask
2614	REAL(r_std) :: reftemp_norefinf
2615	REAL(r_std) :: reftemp_default !! Default value
2616
2617	!Config Key = SOIL_REFTEMP_FILE
2618	!Config Desc = File with climatological soil temperature
2619	!Config If = READ_REFTEMP
2620	!Config Def = reftemp.nc
2621	!Config Help =
2622	!Config Units = [FILE]
2623	filename = 'reftemp.nc'
2624	CALL getin_p('REFTEMP_FILE',filename)
2625
2626	variablename = 'temperature'
2627
2628	IF (printlev >= 1) WRITE(numout,*) "thermosoil_read_reftempfile: Read and interpolate file " &
2629	// TRIM(filename) //" for variable " //TRIM(variablename)
2630
2631	IF (xios_interpolation) THEN
2632
2633	CALL xios_orchidee_recv_field('reftemp_interp',reftemp(:,1))
2634
2635	DO ib=1, kjpindex
2636	reftemp(ib,:) = reftemp(ib,1) + ZeroCelsius
2637	END DO
2638	areftemp = 1.0
2639	ELSE
2640
2641
2642	! For this case there are not types/categories. We have 'only' a continuos field
2643	! Assigning values to vmin, vmax
2644
2645	vmin = 0.
2646	vmax = 9999.
2647
2648	! For this file we do not need neightbours!
2649	neighbours = 0
2650
2651	!! Variables for interpweight
2652	! Type of calculation of cell fractions
2653	fractype = 'default'
2654	! Name of the longitude and latitude in the input file
2655	lonname = 'nav_lon'
2656	latname = 'nav_lat'
2657	! Default value when no value is get from input file
2658	reftemp_default = 1.
2659	! Reference value when no value is get from input file
2660	reftemp_norefinf = 1.
2661	! Should negative values be set to zero from input file?
2662	nonegative = .FALSE.
2663	! Type of mask to apply to the input data (see header for more details)
2664	maskingtype = 'nomask'
2665	! Values to use for the masking (here not used)
2666	maskvals = (/ undef_sechiba, undef_sechiba, undef_sechiba /)
2667	! Name of the variable with the values for the mask in the input file (only if maskkingtype='var') (here not used)
2668	namemaskvar = ''
2669
2670	CALL interpweight_2Dcont(kjpindex, 0, 0, lalo, resolution, neighbours, &
2671	contfrac, filename, variablename, lonname, latname, vmin, vmax, nonegative, maskingtype, &
2672	maskvals, namemaskvar, -1, fractype, reftemp_default, reftemp_norefinf, &
2673	reftemp_file, areftemp)
2674	IF (printlev >= 5) WRITE(numout,*)' thermosoil_read_reftempfile after interpweight_2Dcont'
2675
2676	! Copy reftemp_file temperature to all ground levels and transform into Kelvin
2677	DO ib=1, kjpindex
2678	reftemp(ib, :) = reftemp_file(ib)+ZeroCelsius
2679	END DO
2680
2681	END IF
2682
2683	! Write diagnostics
2684	CALL xios_orchidee_send_field("areftemp",areftemp)
2685
2686	END SUBROUTINE thermosoil_read_reftempfile
2687
2688	!! ================================================================================================================================
2689	!! SUBROUTINE : thermosoil_read_refsoc_file
2690	!!
2691	!>\BRIEF
2692	!!
2693	!! DESCRIPTION : Read file of soil organic carbon to be used in thermix
2694	!! (insulating effect)
2695	!!
2696	!!
2697	!! RECENT CHANGE(S) : None
2698	!!
2699	!! MAIN OUTPUT VARIABLE(S): refsoc : soil organic carbon from data
2700	!!
2701	!! REFERENCE(S) :
2702	!!
2703	!! FLOWCHART : None
2704	!! \n
2705	!_ ================================================================================================================================
2706
2707	SUBROUTINE thermosoil_read_refsoc_file(nbpt, lalo, neighbours, resolution, contfrac)
2708
2709	!! 0. Variable and parameter declaration
2710
2711	!! 0.1 Input variables
2712
2713	INTEGER(i_std), INTENT(in) :: nbpt !! Number of points for which the data needs to be interpolated (unitless)
2714	REAL(r_std), INTENT(in) :: lalo(nbpt,2) !! Vector of latitude and longitudes (degree)
2715	INTEGER(i_std), INTENT(in) :: neighbours(nbpt,NbNeighb)!! Vector of neighbours for each grid point (1=N,2=E,3=S,4=W)
2716	REAL(r_std), INTENT(in) :: resolution(nbpt,2) !! The size of each grid cell in X and Y (km)
2717	REAL(r_std), INTENT(in) :: contfrac(nbpt) !! Fraction of land in each grid cell (unitless)
2718
2719	!! 0.4 Local variables
2720
2721	INTEGER(i_std) :: nbvmax !! nbvmax for interpolation (unitless)
2722	CHARACTER(LEN=80) :: filename
2723	INTEGER(i_std) :: iml, jml, lml, tml !! Indices
2724	INTEGER(i_std) :: fid, ib, ip, jp, fopt !! Indices
2725	INTEGER(i_std) :: ilf, ks !! Indices
2726	REAL(r_std) :: totarea !! Help variable to compute average SOC
2727	REAL(r_std), ALLOCATABLE, DIMENSION(:) :: lat_lu, lon_lu !! Latitudes and longitudes read from input file
2728	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: lat_rel, lon_rel !! Help variable to read file data and allocate memory
2729	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: mask_lu !! Help variable to read file data and allocate memory
2730	INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:) :: mask
2731	REAL(r_std), ALLOCATABLE, DIMENSION(:,:,:) :: refsoc_file !! Help variable to read file data and allocate memory
2732	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: sub_area !! Help variable to read file data and allocate memory
2733	INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:,:) :: sub_index !! Help variable to read file data and allocate memory
2734	CHARACTER(LEN=30) :: callsign !! Help variable to read file data and allocate memory
2735	CHARACTER(LEN=100) :: str !! Temporary string var
2736	LOGICAL :: ok_interpol !! Optional return of aggregate_2d
2737	INTEGER :: ALLOC_ERR !! Help varialbe to count allocation error
2738	!_
2739	!================================================================================================================================
2740
2741	!! 1. Open file and allocate memory
2742
2743	! Open file with SOC map
2744
2745	!Config Key = SOIL_REFSOC_FILE
2746	!Config Desc = File with soil carbon stocks
2747	!Config If = OK_SOIL_CARBON_DISCRETIZATION, USE_REFSOC, SOIL_CTEMPDIFF
2748	!Config Def = refSOC.nc
2749	!Config Help =
2750	!Config Units = [FILE]
2751	filename = 'refSOC.nc'
2752	CALL getin_p('SOIL_REFSOC_FILE',filename)
2753
2754	! Read data from file
2755	IF (is_root_prc) CALL flininfo(filename, iml, jml, lml, tml, fid)
2756	CALL bcast(iml)
2757	CALL bcast(jml)
2758	CALL bcast(lml)
2759	CALL bcast(tml)
2760
2761	IF (lml .NE. ngrnd) THEN
2762	WRITE(str, *) 'ngrnd=', ngrnd, ', depth found in file=', lml
2763	CALL ipslerr_p(3, 'thermosoil_read_refsoc_file', &
2764	'depth from the file must be the same as ngrnd', &
2765	str, &
2766	filename )
2767	ENDIF
2768
2769	ALLOCATE(lon_lu(iml), STAT=ALLOC_ERR)
2770	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'thermosoil_read_refsoc_file','Problem in allocation of variable lon_lu','','')
2771
2772	ALLOCATE(lat_lu(jml), STAT=ALLOC_ERR)
2773	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'thermosoil_read_refsoc_file','Problem in allocation of variable lat_lu','','')
2774
2775	ALLOCATE(mask_lu(iml,jml), STAT=ALLOC_ERR)
2776	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'thermosoil_read_refsoc_file','Pb in allocation for mask_lu','','')
2777
2778	ALLOCATE(refsoc_file(iml,jml,lml), STAT=ALLOC_ERR)
2779	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'thermosoil_read_refsoc_file','Pb in allocation for refsoc_file','','')
2780
2781	IF (is_root_prc) THEN
2782	CALL flinget(fid, 'longitude', iml, 0, 0, 0, 1, 1, lon_lu)
2783	CALL flinget(fid, 'latitude', jml, 0, 0, 0, 1, 1, lat_lu)
2784	CALL flinget(fid, 'mask', iml, jml, 0, 0, 1, 1, mask_lu)
2785	CALL flinget(fid, 'soil_organic_carbon', iml, jml, lml, tml, 1, 1, refsoc_file)
2786
2787	CALL flinclo(fid)
2788	ENDIF
2789
2790	CALL bcast(lon_lu)
2791	CALL bcast(lat_lu)
2792	CALL bcast(mask_lu)
2793	CALL bcast(refsoc_file)
2794
2795	ALLOCATE(lon_rel(iml,jml), STAT=ALLOC_ERR)
2796	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'thermosoil_read_refsoc_file','Pb in allocation for lon_rel','','')
2797
2798	ALLOCATE(lat_rel(iml,jml), STAT=ALLOC_ERR)
2799	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'thermosoil_read_refsoc_file','Pb in allocation for lat_rel','','')
2800
2801	ALLOCATE(mask(iml,jml), STAT=ALLOC_ERR)
2802	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'thermosoil_read_refsoc_file','Problem in allocation of variable mask','','')
2803
2804	DO jp=1,jml
2805	lon_rel(:,jp) = lon_lu(:)
2806	ENDDO
2807	DO ip=1,iml
2808	lat_rel(ip,:) = lat_lu(:)
2809	ENDDO
2810
2811	mask(:,:) = zero
2812	WHERE (mask_lu(:,:) > zero )
2813	mask(:,:) = un
2814	ENDWHERE
2815
2816	! Set nbvmax to 200 for interpolation
2817	! This number is the dimension of the variables in which we store
2818	! the list of points of the source grid which fit into one grid box of the
2819	! target.
2820	nbvmax = 16
2821	callsign = 'soil organic carbon'
2822
2823	! Start interpolation
2824	ok_interpol=.FALSE.
2825	DO WHILE ( .NOT. ok_interpol )
2826	WRITE(numout,*) "Projection arrays for ",callsign," : "
2827	WRITE(numout,*) "nbvmax = ",nbvmax
2828
2829	ALLOCATE(sub_area(nbpt,nbvmax), STAT=ALLOC_ERR)
2830	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'thermosoil_read_refsoc_file','Pb in allocation for sub_area','','')
2831	sub_area(:,:)=zero
2832
2833	ALLOCATE(sub_index(nbpt,nbvmax,2), STAT=ALLOC_ERR)
2834	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'thermosoil_read_refsoc_file','Pb in allocation for sub_index','','')
2835	sub_index(:,:,:)=0
2836
2837	CALL aggregate_p(nbpt, lalo, neighbours, resolution, contfrac, &
2838	iml, jml, lon_rel, lat_rel, mask, callsign, &
2839	nbvmax, sub_index, sub_area, ok_interpol)
2840
2841	IF ( .NOT. ok_interpol ) THEN
2842	DEALLOCATE(sub_area)
2843	DEALLOCATE(sub_index)
2844	nbvmax = nbvmax * 2
2845	ENDIF
2846	ENDDO
2847
2848	! Compute the average
2849	refsoc(:,:) = zero
2850	DO ib = 1, nbpt
2851	fopt = COUNT(sub_area(ib,:) > zero)
2852	IF ( fopt > 0 ) THEN
2853	totarea = zero
2854	DO ilf = 1, fopt
2855	ip = sub_index(ib,ilf,1)
2856	jp = sub_index(ib,ilf,2)
2857	refsoc(ib,:) = refsoc(ib,:) + refsoc_file(ip,jp,:) * sub_area(ib,ilf)
2858	totarea = totarea + sub_area(ib,ilf)
2859	ENDDO
2860	! Normalize
2861	refsoc(ib,:) = refsoc(ib,:)/totarea
2862	ELSE
2863	! Set defalut value for points where the interpolation fail
2864	WRITE(numout,*) 'On point ', ib, ' no points were found for interpolation data. Mean value is used.'
2865	WRITE(numout,*) 'Location : ', lalo(ib,2), lalo(ib,1)
2866	refsoc(ib,:) = 0.
2867	ENDIF
2868	ENDDO
2869
2870	DEALLOCATE (lat_lu)
2871	DEALLOCATE (lat_rel)
2872	DEALLOCATE (lon_lu)
2873	DEALLOCATE (lon_rel)
2874	DEALLOCATE (mask_lu)
2875	DEALLOCATE (mask)
2876	DEALLOCATE (refsoc_file)
2877	DEALLOCATE (sub_area)
2878	DEALLOCATE (sub_index)
2879
2880	END SUBROUTINE thermosoil_read_refsoc_file
2881
2882
2883	!!
2884	!================================================================================================================================
2885	!! SUBROUTINE : thermosoil_add_heat_zimov
2886	!!
2887	!>\BRIEF heat
2888	!!
2889	!! DESCRIPTION :
2890	!!
2891	!! RECENT CHANGE(S) : None
2892	!!
2893	!! MAIN OUTPUT VARIABLE(S):
2894	!!
2895	!! REFERENCE(S) :
2896	!!
2897	!! FLOWCHART : None
2898	!! \n
2899	!_
2900	!================================================================================================================================
2901	SUBROUTINE thermosoil_add_heat_zimov(kjpindex, veget_max, ptn, heat_Zimov)
2902	!! 0. Variables and parameter declaration
2903
2904	!! 0.1 Input variables
2905	INTEGER(i_std),INTENT(in) :: kjpindex
2906	REAL(r_std), DIMENSION (kjpindex,nvm), INTENT(in) :: veget_max !! Fraction of vegetation type
2907
2908	REAL(r_std), DIMENSION(kjpindex,ngrnd,nvm), INTENT (in) :: heat_Zimov !! heating associated with decomposition [W/m**3 soil]
2909
2910	!! 0.2 Modified variables
2911	REAL(r_std),DIMENSION(kjpindex,ngrnd,nvm),INTENT(inout) :: ptn !! Soil temperature profile (K)
2912
2913	!! 0.3 Local variables
2914	INTEGER(r_std) :: ji, jg, jv
2915
2916	IF (printlev>=3) WRITE (numout,*) 'entering thermosoil_add_heat_zimov'
2917
2918	DO ji = 1, kjpindex
2919	DO jv = 1,nvm
2920	DO jg = 1, ngrnd
2921	! BG: should we use pcapa without pft dimension or pkappa_per_pft since they seem similar?
2922	ptn(ji,jg,jv) = ptn(ji,jg,jv) + heat_zimov(ji,jg,jv) * dt_sechiba / ( pcapa(ji,jg) * dlt(jg) )
2923	END DO
2924	END DO
2925	END DO
2926
2927	! ptn_pftmean needs to be updated to ensure consistency
2928	ptn_pftmean(:,:) = zero
2929	DO jv=1,nvm
2930	DO jg = 1, ngrnd
2931	ptn_pftmean(:,jg) = ptn_pftmean(:,jg) + ptn(:,jg,jv) * veget_max(:,jv)
2932	ENDDO ! jg = 1, ngrnd
2933	ENDDO ! m=1,nvm
2934
2935	IF (printlev>=3) WRITE (numout,*) ' thermosoil_add_heat_zimov done'
2936
2937	END SUBROUTINE thermosoil_add_heat_zimov
2938
2939
2940	END MODULE thermosoil

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

Download in other formats: