Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

thermosoil.f90 @ 5501

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

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

Context Navigation

source: CONFIG_DEVT/IPSLCM6.5_work_ENSEMBLES/modeles/ORCHIDEE/src_sechiba/thermosoil.f90 @ 5501

Download in other formats: