Context Navigation

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

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

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

Download in other formats: