Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

thermosoil.f90 @ 7541

Last change on this file since 7541 was 7541, checked in by fabienne.maignan, 2 years ago
Zhang publication on coupling factor
File size: 100.9 KB

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

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

Context Navigation

source: branches/publications/ORCHIDEE_2.2_r7266/ORCHIDEE/src_sechiba/thermosoil.f90 @ 7541

Download in other formats: