Context Navigation

source: tags/ORCHIDEE_4_1/ORCHIDEE/src_sechiba/thermosoil.f90 @ 8119

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

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

Download in other formats: