Context Navigation

thermosoil.f90 @ 8557

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

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

Context Navigation

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

Download in other formats: