Context Navigation

source: branches/publications/ORCHIDEE-LEAK-r5919/src_sechiba/thermosoil.f90 @ 7346

Last change on this file since 7346 was 3345, checked in by josefine.ghattas, 8 years ago
Set READ_REFTEMP=y as default if OK_FREEZE Only call read_reftemp if the variable ptn was not set in restart file. Before the file was always read but the variable ptn was only initialized if it was not found in the restart file. No change in results. Remove option READ_PERMAFROST_MAP and subroutine read_permafrostmap : this option is not usefull as it only read a file without using the variables.
Property svn:keywords set to `Revision Date HeadURL Date Author Revision`
File size: 96.3 KB

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

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

Download in other formats: