Context Navigation

source: branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/thermosoil.f90 @ 7710

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

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

Download in other formats: