Context Navigation

source: tags/ORCHIDEE_2_0/ORCHIDEE/src_sechiba/thermosoilc.f90 @ 5164

Last change on this file since 5164 was 5096, checked in by josefine.ghattas, 6 years ago
Protected calculation and output of snowtemp_weighted done only for explicitsnow.
Property svn:keywords set to `Revision Date HeadURL Date Author Revision`
File size: 105.0 KB

Line
1	! =================================================================================================================================
2	! MODULE : thermosoilc
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 Choisnel hydrology.
11	!!
12	!!\n DESCRIPTION : General important informations about the numerical scheme and
13	!! the soil vertical discretization:\n
14	!! - the soil is divided into "ngrnd" (=7 by default) layers, reaching to as
15	!! deep as 5.5m down within the soil, with thiscknesses
16	!! following a geometric series of ration 2.\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=soilflux+otherfluxes(Ts(t))\n
46	!! -- EQ3 --\n
47	!!
48	!! and the dependency from the previous timestep is hidden in
49	!! soilcap and soilflux (apparent surface heat capacity and heat
50	!! flux respectively). Soilcap and soilflux are therefore
51	!! calculated at the previsou 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 = (zz_coeff(1))/((zz_coef(2)-zz_coef(1))) \n
60	!!
61	!! - cstgrnd is the attenuation depth of the diurnal temperature signal
62	!! (period : one_day) as a result of the heat conduction equation
63	!! with no coefficients :
64	!!\latexonly
65	!!\input{thermosoil_var_init0.tex}
66	!!\endlatexonly
67	!! -- EQ4 --\n
68	!! This equation results from the change of variables :
69	!! z' =z*sqrt(Cp/K) where z' is the new depth (homogeneous
70	!! to sqrt(time) ), z the real depth (in m), Cp and K the soil heat
71	!! capacity and conductivity respectively.\n
72	!!
73	!! the attenuation depth of a diurnal thermal signal for EQ4 is therefore homogeneous to sqrt(time) and
74	!! equals : \n
75	!! cstgrnd = sqrt(oneday/Pi)
76	!!
77	!! - lskin is the attenuation depth of the diurnal temperature signal
78	!! (period : one_day) within the soil for the complete heat conduction equation
79	!! (ie : with coefficients)
80	!!\latexonly
81	!!\input{thermosoil_var_init00.tex}
82	!!\endlatexonly
83	!! -- EQ5 -- \n
84	!! it can be retrieved from cstgrnd using the change of variable z' =z*sqrt(Cp/K):\n
85	!! lskin = sqrt(K/Cp)cstgrnd = sqrt(K/Cp)sqrt(oneday//Pi)\n
86	!!
87	!! In thermosoil, the ratio lskin/cstgrnd is frequently used as the
88	!! multiplicative factor to go from
89	!!'adimensional' depths (like z' ) to real depths (z). z' is not really
90	!! adimensional but is reffered to like this in the code.
91	!!
92	!!
93	!! RECENT CHANGE(S) : thermosoilc module is a copy of thermosoil before the vertical discretization was changed.
94	!! This module is used only for the Choisnel hydrology.
95	!!
96	!! REFERENCE(S) : None
97	!!
98	!! SVN :
99	!! $HeadURL$
100	!! $Date$
101	!! $Revision$
102	!! \n
103	!_ ================================================================================================================================
104
105	MODULE thermosoilc
106
107	! modules used :
108	USE ioipsl
109	USE ioipsl_para
110	USE xios_orchidee
111	USE constantes
112	USE time, ONLY : one_day, dt_sechiba
113	USE constantes_soil
114	USE sechiba_io_p
115	USE grid
116
117
118	IMPLICIT NONE
119
120	!private and public routines :
121	PRIVATE
122	PUBLIC :: thermosoilc_main, thermosoilc_clear, thermosoilc_levels, thermosoilc_initialize, thermosoilc_finalize
123
124	REAL(r_std), SAVE :: lambda, cstgrnd, lskin !! See Module description
125	!$OMP THREADPRIVATE(lambda, cstgrnd, lskin)
126
127	REAL(r_std), SAVE :: fz1, zalph !! usefull constants for diverse use
128	!$OMP THREADPRIVATE(fz1, zalph)
129	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: ptn !! vertically discretized
130	!! soil temperatures @tex ($K$) @endtex.
131	!$OMP THREADPRIVATE(ptn)
132	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: zz !! depths of the soil thermal numerical nodes.
133	!! Caveats: they are not exactly the centers of the
134	!! thermal layers, see the calculation in
135	!! ::thermosoilc_var_init @tex ($m$) @endtex.
136	!$OMP THREADPRIVATE(zz)
137	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: zz_coef !! depths of the boundaries of the thermal layers,
138	!! see the calculation in
139	!! thermosoilc_var_init @tex ($m$) @endtex.
140	!$OMP THREADPRIVATE(zz_coef)
141	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: dz1 !! numerical constant used in the thermal numerical
142	!! scheme @tex ($m^{-1}$) @endtex. ; it corresponds
143	!! to the coefficient @tex $d_k$ @endtex of equation
144	!! (A.12) in F. Hourdin PhD thesis.
145	!$OMP THREADPRIVATE(dz1)
146	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: dz2 !! thicknesses of the thermal layers @tex ($m$)
147	!! @endtex; typically:
148	!! dz2(jg)=zz_coef(jg+1)-zz_coef(jg); calculated once
149	!! and for all in thermosoilc_var_init
150	!$OMP THREADPRIVATE(dz2)
151	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: cgrnd !! integration coefficient for the numerical scheme,
152	!! see eq.1
153	!$OMP THREADPRIVATE(cgrnd)
154	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: dgrnd !! integration coefficient for the numerical scheme,
155	!! see eq.1
156	!$OMP THREADPRIVATE(dgrnd)
157
158	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: pcapa !! volumetric vertically discretized soil heat
159	!! capacity @tex ($J K^{-1} m^{-3}$) @endtex.
160	!! It depends on the soil
161	!! moisture content (shum_ngrnd_perma) and is calculated at
162	!! each time step in thermosoilc_coef.
163	!$OMP THREADPRIVATE(pcapa)
164	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: pkappa !! vertically discretized soil thermal conductivity
165	!! @tex ($W K^{-1} m^{-1}$) @endtex. Same as pcapa.
166	!$OMP THREADPRIVATE(pkappa)
167	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: pcapa_snow !! volumetric vertically discretized snow heat
168	!! capacity @tex ($J K^{-1} m^{-3}$) @endtex.
169	!$OMP THREADPRIVATE(pcapa_snow)
170	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: pkappa_snow !! vertically discretized snow thermal conductivity
171	!! @tex ($W K^{-1} m^{-1}$) @endtex.
172	!$OMP THREADPRIVATE(pkappa_snow)
173
174	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: pcapa_en !! heat capacity used for surfheat_incr and
175	!! coldcont_incr
176	!$OMP THREADPRIVATE(pcapa_en)
177	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: ptn_beg !! vertically discretized temperature at the
178	!! beginning of the time step @tex ($K$) @endtex;
179	!! is used in
180	!! thermosoilc_energy for energy-related diagnostic of
181	!! the routine.
182	!$OMP THREADPRIVATE(ptn_beg)
183	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: temp_sol_beg !! Surface temperature at the beginning of the
184	!! timestep @tex ($K$) @endtex
185	!$OMP THREADPRIVATE(temp_sol_beg)
186	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: surfheat_incr !! Change in soil heat content during the timestep
187	!! @tex ($J$) @endtex.
188	!$OMP THREADPRIVATE(surfheat_incr)
189	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: coldcont_incr !! Change in snow heat content @tex ($J$) @endtex.
190	!$OMP THREADPRIVATE(coldcont_incr)
191	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: shum_ngrnd_perma !! Saturation degree on the thermal axes (0-1, dimensionless)
192	!$OMP THREADPRIVATE(shum_ngrnd_perma)
193
194	! Variables related to soil freezing
195	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: profil_froz !! Frozen fraction of the soil on hydrological levels (-)
196	!$OMP THREADPRIVATE(profil_froz)
197	REAL(r_std),ALLOCATABLE, SAVE, DIMENSION (:) :: e_soil_lat !! Accumulated latent heat for the whole soil (J)
198	!$OMP THREADPRIVATE(e_soil_lat)
199	REAL(r_std),ALLOCATABLE, SAVE, DIMENSION (:,:) :: pcappa_supp !! Additional heat capacity due to soil freezing for each soil layer (J/K)
200	!$OMP THREADPRIVATE(pcappa_supp)
201
202
203	CONTAINS
204	!! =============================================================================================================================
205	!! SUBROUTINE : thermosoilc_initialize
206	!!
207	!>\BRIEF Allocate module variables, read from restart file or initialize with default values
208	!!
209	!! DESCRIPTION : Allocate module variables, read from restart file or initialize with default values.
210	!! Call thermosoilc_var_init to calculate physical constants.
211	!! Call thermosoilc_coef to calculate thermal soil properties.
212	!!
213	!! RECENT CHANGE(S) : None
214	!!
215	!! REFERENCE(S) : None
216	!!
217	!! FLOWCHART : None
218	!! \n
219	!_ ==============================================================================================================================
220	SUBROUTINE thermosoilc_initialize (kjit, kjpindex, rest_id, &
221	temp_sol_new, snow, shumdiag_perma, &
222	soilcap, soilflx, stempdiag, &
223	gtemp, &
224	frac_snow_veg,frac_snow_nobio,totfrac_nobio, &
225	snowdz, snowrho, snowtemp, lambda_snow, cgrnd_snow, dgrnd_snow, pb)
226
227	!! 0. Variable and parameter declaration
228	!! 0.1 Input variables
229	INTEGER(i_std), INTENT(in) :: kjit !! Time step number (unitless)
230	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (unitless)
231	INTEGER(i_std),INTENT (in) :: rest_id !! Restart file identifier (unitless)
232	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_sol_new !! Surface temperature at the present time-step,
233	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: snow !! Snow mass (kg)
234	REAL(r_std),DIMENSION (kjpindex,nslm), INTENT (in) :: shumdiag_perma !! Soil saturation degree on the diagnostic axis (0-1, unitless)
235	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: frac_snow_veg !! Snow cover fraction on vegeted area
236	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT(in) :: frac_snow_nobio !! Snow cover fraction on non-vegeted area
237	REAL(r_std),DIMENSION (kjpindex),INTENT(in) :: totfrac_nobio !! Total fraction of continental ice+lakes+cities+...
238	!! (unitless,0-1)
239	REAL(r_std),DIMENSION (kjpindex,nsnow),INTENT (in) :: snowdz !! Snow depth
240	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(in) :: snowrho !! Snow density
241	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(in) :: snowtemp !! Snow temperature
242	REAL(r_std), DIMENSION (kjpindex), INTENT (in) :: pb !! Surface presure (hPa)
243
244
245	!! 0.2 Output variables
246	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: soilcap !! apparent surface heat capacity (J m-2 K-1)
247	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: soilflx !! apparent soil heat flux (W m-2)
248	REAL(r_std),DIMENSION (kjpindex,nslm), INTENT (out) :: stempdiag !! diagnostic temperature profile (K)
249	REAL(r_std),DIMENSION (kjpindex),INTENT(out) :: gtemp !! First soil layer temperature
250
251	!! 0.3 Modified variables
252	REAL(r_std), DIMENSION (kjpindex), INTENT(inout) :: lambda_snow !! Coefficient of the linear extrapolation of surface temperature
253	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT (inout):: cgrnd_snow !! Integration coefficient for snow numerical scheme
254	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT (inout):: dgrnd_snow !! Integration coefficient for snow numerical scheme
255
256	!! 0.4 Local variables
257	INTEGER(i_std) :: ier, i
258	LOGICAL :: calculate_coef !! Local flag to initialize variables by call to thermosoilc_coef
259	!_ ================================================================================================================================
260
261	IF (printlev >= 3) WRITE (numout,*) 'Start thermosoilc_initialize '
262
263	!! 1. Allocate soil temperatures variables
264	ALLOCATE (ptn(kjpindex,ngrnd),stat=ier)
265	IF (ier /= 0) CALL ipslerr_p(3,'thermosoilc_initialize', 'Error in allocation of ptn','','')
266
267	ALLOCATE (zz(ngrnd),stat=ier)
268	IF (ier /= 0) CALL ipslerr_p(3,'thermosoilc_initialize', 'Error in allocation of zz','','')
269
270	ALLOCATE (zz_coef(ngrnd),stat=ier)
271	IF (ier /= 0) CALL ipslerr_p(3,'thermosoilc_initialize', 'Error in allocation of zz_coef','','')
272
273	ALLOCATE (dz1(ngrnd),stat=ier)
274	IF (ier /= 0) CALL ipslerr_p(3,'thermosoilc_initialize', 'Error in allocation of dz1','','')
275
276	ALLOCATE (dz2(ngrnd),stat=ier)
277	IF (ier /= 0) CALL ipslerr_p(3,'thermosoilc_initialize', 'Error in allocation of dz2','','')
278
279	ALLOCATE (cgrnd(kjpindex,ngrnd-1),stat=ier)
280	IF (ier /= 0) CALL ipslerr_p(3,'thermosoilc_initialize', 'Error in allocation of cgrnd','','')
281
282	ALLOCATE (dgrnd(kjpindex,ngrnd-1),stat=ier)
283	IF (ier /= 0) CALL ipslerr_p(3,'thermosoilc_initialize', 'Error in allocation of dgrnd','','')
284
285	ALLOCATE (pkappa_snow(kjpindex,nsnow),stat=ier)
286	IF (ier /= 0) CALL ipslerr_p(3,'thermosoilc_initialize', 'Error in allocation of pkappa_snow','','')
287
288	ALLOCATE (pcapa_snow(kjpindex,nsnow),stat=ier)
289	IF (ier /= 0) CALL ipslerr_p(3,'thermosoilc_initialize', 'Error in allocation of pcapa_snow','','')
290
291	ALLOCATE (pcapa(kjpindex,ngrnd),stat=ier)
292	IF (ier /= 0) CALL ipslerr_p(3,'thermosoilc_initialize', 'Error in allocation of pcapa','','')
293
294	ALLOCATE (pkappa(kjpindex,ngrnd),stat=ier)
295	IF (ier /= 0) CALL ipslerr_p(3,'thermosoilc_initialize', 'Error in allocation of pkappa','','')
296
297	! Temporary fix: Initialize following variable because they are output to xios before the first calculation
298	pcapa = 0
299	pkappa = 0
300	pcapa_snow = 0
301	pkappa_snow = 0
302
303	ALLOCATE (surfheat_incr(kjpindex),stat=ier)
304	IF (ier /= 0) CALL ipslerr_p(3,'thermosoilc_initialize', 'Error in allocation of surfheat_incr','','')
305
306	ALLOCATE (coldcont_incr(kjpindex),stat=ier)
307	IF (ier /= 0) CALL ipslerr_p(3,'thermosoilc_initialize', 'Error in allocation of coldcont_incr','','')
308
309	ALLOCATE (pcapa_en(kjpindex,ngrnd),stat=ier)
310	IF (ier /= 0) CALL ipslerr_p(3,'thermosoilc_initialize', 'Error in allocation of pcapa_en','','')
311
312	ALLOCATE (ptn_beg(kjpindex,ngrnd),stat=ier)
313	IF (ier /= 0) CALL ipslerr_p(3,'thermosoilc_initialize', 'Error in allocation of ptn_beg','','')
314
315	ALLOCATE (temp_sol_beg(kjpindex),stat=ier)
316	IF (ier /= 0) CALL ipslerr_p(3,'thermosoilc_initialize', 'Error in allocation of temp_sol_beg','','')
317
318	ALLOCATE (shum_ngrnd_perma(kjpindex,ngrnd),stat=ier)
319	IF (ier /= 0) CALL ipslerr_p(3,'thermosoilc_initialize', 'Error in allocation of shum_ngrnd_perma','','')
320
321	ALLOCATE (profil_froz(kjpindex,ngrnd),stat=ier)
322	IF (ier /= 0) CALL ipslerr_p(3,'thermosoilc_initialize', 'Error in allocation of profil_froz','','')
323
324	IF (ok_freeze_thermix) THEN
325	ALLOCATE (pcappa_supp(kjpindex,ngrnd),stat=ier)
326	IF (ier /= 0) CALL ipslerr_p(3,'thermosoilc_initialize', 'Error in allocation of ok_freeze_termix','','')
327	END IF
328	IF (ok_Ecorr) THEN
329	ALLOCATE (e_soil_lat(kjpindex),stat=ier)
330	IF (ier /= 0) CALL ipslerr_p(3,'thermosoilc_initialize', 'Error in allocation of e_soil_lat','','')
331	END IF
332
333	!! 2. Initialize variable from restart file or with default values
334
335	!! Reads restart files for soil temperatures only. If no restart file is
336	!! found, the initial soil temperature is by default set to 280K at all depths. The user
337	!! can decide to initialize soil temperatures at an other value, in which case he should set the flag THERMOSOIL_TPRO
338	!! to this specific value in the run.def.
339	IF (printlev>=3) WRITE (numout,*) 'Read restart file for THERMOSOIL variables'
340
341	CALL ioconf_setatt_p('UNITS', 'K')
342	CALL ioconf_setatt_p('LONG_NAME','Soil Temperature profile')
343	CALL restget_p (rest_id, 'ptn', nbp_glo, ngrnd, 1, kjit, .TRUE., ptn, "gather", nbp_glo, index_g)
344
345	! Initialize ptn if it was not found in restart file
346	IF (ALL(ptn(:,:)==val_exp)) THEN
347	! ptn was not found in restart file
348
349	IF (read_reftemp) THEN
350	! Read variable ptn from file
351	CALL thermosoilc_read_reftempfile(kjpindex,lalo,ptn)
352	ELSE
353	! Initialize ptn with a constant value which can be set in run.def
354
355	!Config Key = THERMOSOIL_TPRO
356	!Config Desc = Initial soil temperature profile if not found in restart
357	!Config Def = 280.
358	!Config If = OK_SECHIBA
359	!Config Help = The initial value of the temperature profile in the soil if
360	!Config its value is not found in the restart file. Here
361	!Config we only require one value as we will assume a constant
362	!Config throughout the column.
363	!Config Units = Kelvin [K]
364	CALL setvar_p (ptn, val_exp,'THERMOSOIL_TPRO',280._r_std)
365	END IF
366	END IF
367
368	! Initialize ptn_beg (variable needed in thermosoilc_coef before calucation in thermosoilc_energy)
369	ptn_beg(:,:) = ptn(:,:)
370
371	! Initialize temp_sol_beg with values from previous time-step
372	temp_sol_beg(:) = temp_sol_new(:)
373
374	! Read e_soil_lat from restart file or initialize
375	IF (ok_Ecorr) THEN
376	CALL restget_p (rest_id, 'e_soil_lat', nbp_glo, 1, 1, kjit, .TRUE., &
377	e_soil_lat, "gather", nbp_glo, index_g)
378	CALL setvar_p (e_soil_lat, val_exp,'NO_KEYWORD',zero)
379	END IF
380
381	! Read gtemp from restart file
382	CALL restget_p (rest_id, 'gtemp', nbp_glo, 1, 1, kjit, .TRUE., &
383	gtemp, "gather", nbp_glo, index_g)
384	CALL setvar_p (gtemp, val_exp,'NO_KEYWORD',zero)
385
386	! Read variables calculated in thermosoilc_coef from restart file
387	! If the variables were not found in the restart file, the logical
388	! calculate_coef will be true and thermosoilc_coef will be called further below.
389	! These variables need to be in the restart file to avoid a time shift that
390	! would be done using thermosoilc_coef at this stage.
391	calculate_coef=.FALSE.
392	CALL ioconf_setatt_p('UNITS', 'J m-2 K-1')
393	CALL ioconf_setatt_p('LONG_NAME','Apparent surface heat capacity')
394	CALL restget_p (rest_id, 'soilcap', nbp_glo, 1, 1, kjit, .TRUE., &
395	soilcap, "gather", nbp_glo, index_g)
396	IF (ALL(soilcap(:)==val_exp)) calculate_coef=.TRUE.
397
398	CALL ioconf_setatt_p('UNITS', 'W m-2')
399	CALL ioconf_setatt_p('LONG_NAME','Apparent soil heat flux')
400	CALL restget_p (rest_id, 'soilflx', nbp_glo, 1, 1, kjit, .TRUE., &
401	soilflx, "gather", nbp_glo, index_g)
402	IF (ALL(soilflx(:)==val_exp)) calculate_coef=.TRUE.
403
404	CALL ioconf_setatt_p('UNITS', '')
405	CALL ioconf_setatt_p('LONG_NAME','Integration coefficient for the numerical scheme')
406	CALL restget_p (rest_id, 'cgrnd', nbp_glo, ngrnd-1, 1, kjit, .TRUE., &
407	cgrnd, "gather", nbp_glo, index_g)
408	IF (ALL(cgrnd(:,:)==val_exp)) calculate_coef=.TRUE.
409
410	CALL ioconf_setatt_p('UNITS', '')
411	CALL ioconf_setatt_p('LONG_NAME','Integration coefficient for the numerical scheme')
412	CALL restget_p (rest_id, 'dgrnd', nbp_glo, ngrnd-1, 1, kjit, .TRUE., &
413	dgrnd, "gather", nbp_glo, index_g)
414	IF (ALL(dgrnd(:,:)==val_exp)) calculate_coef=.TRUE.
415
416	CALL ioconf_setatt_p('UNITS', '')
417	CALL ioconf_setatt_p('LONG_NAME','Integration coefficient for the numerical scheme')
418	CALL restget_p (rest_id, 'cgrnd_snow', nbp_glo, nsnow, 1, kjit, .TRUE., &
419	cgrnd_snow, "gather", nbp_glo, index_g)
420	IF (ALL(cgrnd_snow(:,:)==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, 'dgrnd_snow', nbp_glo, nsnow, 1, kjit, .TRUE., &
425	dgrnd_snow, "gather", nbp_glo, index_g)
426	IF (ALL(dgrnd_snow(:,:)==val_exp)) calculate_coef=.TRUE.
427
428	CALL ioconf_setatt_p('UNITS', '')
429	CALL ioconf_setatt_p('LONG_NAME','Coefficient of the linear extrapolation of surface temperature')
430	CALL restget_p (rest_id, 'lambda_snow', nbp_glo, 1, 1, kjit, .TRUE., &
431	lambda_snow, "gather", nbp_glo, index_g)
432	IF (ALL(lambda_snow(:)==val_exp)) calculate_coef=.TRUE.
433
434
435	!! 2.2.Computes physical constants and arrays; initializes soil thermal properties; produces the first stempdiag
436	!! Computes some physical constants and arrays depending on the soil vertical discretization
437	!! (lskin, cstgrnd, zz, zz_coef, dz1, dz2); get the vertical humidity onto the thermal levels
438	CALL thermosoilc_var_init (kjpindex, shumdiag_perma, snow, &
439	zz, zz_coef, dz1, dz2, stempdiag)
440
441	!! 2.3. Computes cgrd, dgrd, soilflx and soilcap coefficients from restart values or initialisation values.
442	!! This is done only if they were not found in restart file.
443	IF (calculate_coef) THEN
444	IF (printlev>=3) WRITE (numout,*) 'thermosoilc_coef will be called in the intialization phase'
445	CALL thermosoilc_coef (&
446	kjpindex, temp_sol_new, snow, &
447	ptn, soilcap, soilflx, zz, &
448	dz1, dz2, &
449	cgrnd, dgrnd, &
450	frac_snow_veg,frac_snow_nobio,totfrac_nobio, &
451	snowdz,snowrho, snowtemp, pb, &
452	lambda_snow, cgrnd_snow, dgrnd_snow)
453	END IF
454
455	END SUBROUTINE thermosoilc_initialize
456
457
458	!! ================================================================================================================================
459	!! SUBROUTINE : thermosoilc_main
460	!!
461	!>\BRIEF Thermosoilc_main computes the soil thermal properties and dynamics, ie solves
462	!! the heat diffusion equation within the soil. The soil temperature profile is
463	!! then interpolated onto the diagnostic axis.
464	!!
465	!! DESCRIPTION : The resolution of the soil heat diffusion equation
466	!! relies on a numerical finite-difference implicit scheme
467	!! fully described in the reference and in the header of the thermosoilc module.
468	!! - The dependency of the previous timestep hidden in the
469	!! integration coefficients cgrnd and dgrnd (EQ1), calculated in thermosoilc_coef, and
470	!! called at the end of the routine to prepare for the next timestep.
471	!! - The effective computation of the new soil temperatures is performed in thermosoilc_profile.
472	!!
473	!! - thermosoilc_coef calculates the coefficients for the numerical scheme for the very first iteration of thermosoilc;
474	!! after that, thermosoilc_coef is called only at the end of the module to calculate the coefficients for the next timestep.
475	!! - thermosoilc_profile solves the numerical scheme.\n
476	!!
477	!! - Flags : one unique flag : THERMOSOIL_TPRO (to be set to the desired initial soil in-depth temperature in K; by default 280K)
478	!!
479	!! RECENT CHANGE(S) : None
480	!!
481	!! MAIN OUTPUT VARIABLE(S): vertically discretized soil temperatures ptn, soil
482	!! thermal properties (pcapa, pkappa), apparent surface heat capacity (soilcap)
483	!! and heat flux (soilflux) to be used in enerbil at the next timestep to solve
484	!! the surface energy balance.
485	!!
486	!! REFERENCE(S) :
487	!! - Hourdin, F. (1992). Study and numerical simulation of the general circulation of planetary atmospheres,
488	!! Ph.D. thesis, Paris VII University. Remark: the part of F. Hourdin' s PhD thesis relative to the thermal
489	!! integration scheme has been scanned and is provided along with the documentation, with name :
490	!! Hourdin_1992_PhD_thermal_scheme.pdf
491	!!
492	!! FLOWCHART :
493	!! \latexonly
494	!! \includegraphics[scale = 1]{thermosoil_flowchart.png}
495	!! \endlatexonly
496	!!
497	!! \n
498	!_ ================================================================================================================================
499
500	SUBROUTINE thermosoilc_main (kjit, kjpindex, &
501	index, indexgrnd, &
502	indexnslm, &
503	temp_sol_new, snow, soilcap, soilflx, &
504	shumdiag_perma, stempdiag, ptnlev1, rest_id, hist_id, hist2_id, &
505	snowdz,snowrho,snowtemp,gtemp,pb,&
506	frac_snow_veg,frac_snow_nobio,totfrac_nobio,temp_sol_add, &
507	lambda_snow, cgrnd_snow, dgrnd_snow)
508
509	!! 0. Variable and parameter declaration
510
511	!! 0.1 Input variables
512
513	INTEGER(i_std), INTENT(in) :: kjit !! Time step number (unitless)
514	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (unitless)
515	INTEGER(i_std),INTENT (in) :: rest_id,hist_id !! Restart_ file and history file identifier
516	!! (unitless)
517	INTEGER(i_std),INTENT (in) :: hist2_id !! history file 2 identifier (unitless)
518	INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: index !! Indeces of the points on the map (unitless)
519	INTEGER(i_std),DIMENSION (kjpindex*ngrnd), INTENT (in):: indexgrnd !! Indeces of the points on the 3D map (vertical
520	!! dimension towards the ground) (unitless)
521	REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: temp_sol_new !! Surface temperature at the present time-step,
522	!! temp_sol_new is only modified for the case ok_explicitsnow
523	!! Ts @tex ($K$) @endtex
524	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: snow !! Snow mass @tex ($kg$) @endtex.
525	!! Caveat: when there is snow on the
526	!! ground, the snow is integrated into the soil for
527	!! the calculation of the thermal dynamics. It means
528	!! that the uppermost soil layers can completely or
529	!! partially consist in snow. In the second case, zx1
530	!! and zx2 are the fraction of the soil layer
531	!! consisting in snow and 'normal' soil, respectively
532	!! This is calculated in thermosoilc_coef.
533	REAL(r_std),DIMENSION (kjpindex,nslm), INTENT (in) :: shumdiag_perma !! Soil saturation degree on the diagnostic axis (0-1, unitless)
534	INTEGER(i_std),DIMENSION (kjpindex*nslm), INTENT (in) :: indexnslm !! Indeces of the points on the 3D map
535	REAL(r_std),DIMENSION (kjpindex,nsnow),INTENT (in) :: snowdz !! Snow depth
536	REAL(r_std),DIMENSION (kjpindex,nsnow),INTENT (in) :: snowrho !! Snow density
537	REAL(r_std), DIMENSION (kjpindex), INTENT (in) :: pb !! Surface presure (hPa)
538
539	REAL(r_std),DIMENSION (kjpindex,nsnow),INTENT (inout) :: snowtemp !! Snow temperature
540	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: frac_snow_veg !! Snow cover fraction on vegeted area
541	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT(in) :: frac_snow_nobio !! Snow cover fraction on non-vegeted area
542	REAL(r_std),DIMENSION (kjpindex),INTENT(in) :: totfrac_nobio !! Total fraction of continental ice+lakes+cities+...
543	!!(unitless,0-1)
544	REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: temp_sol_add !! additional surface temperature due to the melt of first layer
545	!!at the present time-step @tex ($K$) @endtex
546
547	!! 0.2 Output variables
548
549	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: ptnlev1 !! 1st level soil temperature
550	REAL(r_std),DIMENSION (kjpindex),INTENT(out) :: gtemp !! First soil layer temperature
551
552
553	!! 0.3 Modified variables
554
555	REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: soilcap !! apparent surface heat capacity
556	!! @tex ($J m^{-2} K^{-1}$) @endtex
557	REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: soilflx !! apparent soil heat flux @tex ($W m^{-2}$) @endtex
558	!! , positive
559	!! towards the soil, writen as Qg (ground heat flux)
560	!! in the history files, and computed at the end of
561	!! thermosoilc for the calculation of Ts in enerbil,
562	!! see EQ3.
563	REAL(r_std),DIMENSION (kjpindex,nslm), INTENT (inout) :: stempdiag !! diagnostic temperature profile @tex ($K$) @endtex
564	!! , eg on the
565	!! diagnostic axis (levels:1:nslm). The soil
566	!! temperature is put on this diagnostic axis to be
567	!! used by other modules (slowproc.f90; routing.f90;
568	!! hydrol or hydrolc when a frozen soil
569	!! parametrization is used..)
570	REAL(r_std), DIMENSION (kjpindex), INTENT(inout) :: lambda_snow !! Coefficient of the linear extrapolation of surface temperature
571	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT (inout):: cgrnd_snow !! Integration coefficient for snow numerical scheme
572	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT (inout):: dgrnd_snow !! Integration coefficient for snow numerical scheme
573
574	!! 0.4 Local variables
575
576	INTEGER(i_std) :: jv,ji,ii
577	REAL(r_std),DIMENSION (kjpindex) :: snowtemp_weighted!! Snow temperature weighted by snow density (K)
578
579	!_ ================================================================================================================================
580
581	!! 3. Put the soil wetness diagnostic on the levels of the soil temperature
582
583	!!?? this could logically be put just before the last call to
584	!!thermosoil_coef, as the results are used there...
585	CALL thermosoilc_humlev(kjpindex, shumdiag_perma, snow)
586
587
588	!! 4. Effective computation of the soil temperatures profile, using the cgrd and dgrd coefficients from previsou tstep.
589
590	CALL thermosoilc_profile (kjpindex, temp_sol_new, ptn, &
591	stempdiag,snowtemp, frac_snow_veg, frac_snow_nobio,totfrac_nobio, &
592	cgrnd_snow, dgrnd_snow)
593
594	!! 5. Call to thermosoilc_energy, still to be clarified..
595
596	CALL thermosoilc_energy (kjpindex, temp_sol_new, soilcap)
597
598	!! 6. Writing the history files according to the ALMA standards (or not..)
599
600	!in only one file (hist2_id <=0) or in 2 different files (hist2_id >0).
601	CALL xios_orchidee_send_field("ptn",ptn)
602	CALL xios_orchidee_send_field("soilflx",soilflx)
603	CALL xios_orchidee_send_field("surfheat_incr",surfheat_incr)
604	CALL xios_orchidee_send_field("coldcont_incr",coldcont_incr)
605	CALL xios_orchidee_send_field("pkappa",pkappa)
606	CALL xios_orchidee_send_field("pkappa_snow",pkappa_snow)
607	CALL xios_orchidee_send_field("pcapa",pcapa)
608	CALL xios_orchidee_send_field("pcapa_snow",pcapa_snow)
609	CALL xios_orchidee_send_field("snowtemp",snowtemp)
610
611	IF (ok_explicitsnow) THEN
612	DO ji=1,kjpindex
613	! Use min_sechiba instead of zero to avoid problem with division by zero
614	IF (snow(ji) .GT. min_sechiba) THEN
615	snowtemp_weighted(ji) = SUM(snowtemp(ji,:)*snowrho(ji,:))/SUM(snowrho(ji,:))
616	ELSE
617	snowtemp_weighted(ji) = xios_default_val
618	END IF
619	END DO
620	CALL xios_orchidee_send_field("snowtemp_weighted",snowtemp_weighted)
621	END IF
622
623	IF ( .NOT. almaoutput ) THEN
624	CALL histwrite_p(hist_id, 'ptn', kjit, ptn, kjpindex*ngrnd, indexgrnd)
625	CALL histwrite_p(hist_id, 'Qg', kjit, soilflx, kjpindex, index)
626	ELSE
627	CALL histwrite_p(hist_id, 'SoilTemp', kjit, ptn, kjpindex*ngrnd, indexgrnd)
628	CALL histwrite_p(hist_id, 'Qg', kjit, soilflx, kjpindex, index)
629	CALL histwrite_p(hist_id, 'DelSurfHeat', kjit, surfheat_incr, kjpindex, index)
630	CALL histwrite_p(hist_id, 'DelColdCont', kjit, coldcont_incr, kjpindex, index)
631	ENDIF
632	IF ( hist2_id > 0 ) THEN
633	IF ( .NOT. almaoutput ) THEN
634	CALL histwrite_p(hist2_id, 'ptn', kjit, ptn, kjpindex*ngrnd, indexgrnd)
635	ELSE
636	CALL histwrite_p(hist2_id, 'SoilTemp', kjit, ptn, kjpindex*ngrnd, indexgrnd)
637	CALL histwrite_p(hist2_id, 'Qg', kjit, soilflx, kjpindex, index)
638	CALL histwrite_p(hist2_id, 'DelSurfHeat', kjit, surfheat_incr, kjpindex, index)
639	CALL histwrite_p(hist2_id, 'DelColdCont', kjit, coldcont_incr, kjpindex, index)
640	ENDIF
641	ENDIF
642
643	!! 7. A last final call to thermosoilc_coef
644
645	!! A last final call to thermosoilc_coef, which calculates the different
646	!!coefficients (cgrnd, dgrnd, dz1, soilcap, soilflx) from this time step to be
647	!!used at the next time step, either in the surface temperature calculation
648	!!(soilcap, soilflx) or in the soil thermal numerical scheme.
649	CALL thermosoilc_coef (&
650	kjpindex, temp_sol_new, snow, &
651	ptn, soilcap, soilflx, zz, &
652	dz1, dz2, &
653	cgrnd, dgrnd, &
654	frac_snow_veg,frac_snow_nobio,totfrac_nobio,&
655	snowdz,snowrho, snowtemp, pb, &
656	lambda_snow, cgrnd_snow, dgrnd_snow)
657
658	! Save variables for explicit snow model
659	gtemp(:) = ptn(:,1)
660
661	!! Initialize output arguments to be used in sechiba
662	ptnlev1(:) = ptn(:,1)
663
664
665	!! Surface temperature is forced to zero celcius if its value is larger than melting point, only for explicit snow scheme
666	IF (ok_explicitsnow) THEN
667	DO ji=1,kjpindex
668	IF (SUM(snowdz(ji,:)) .GT. 0.0) THEN
669	IF (temp_sol_new(ji) .GE. tp_00) THEN
670	temp_sol_new(ji) = tp_00
671	ENDIF
672	END IF
673	END DO
674	ENDIF
675
676
677	IF (printlev>=3) WRITE (numout,*) ' thermosoilc_main done '
678
679	END SUBROUTINE thermosoilc_main
680
681	!! =============================================================================================================================
682	!! SUBROUTINE : thermosoilc_finalize
683	!!
684	!>\BRIEF Write to restart file
685	!!
686	!! DESCRIPTION : This subroutine writes the module variables and variables calculated in thermosoilc
687	!! to restart file
688	!! \n
689	!_ ==============================================================================================================================
690	SUBROUTINE thermosoilc_finalize (kjit, kjpindex, rest_id, gtemp, &
691	soilcap, soilflx, lambda_snow, cgrnd_snow, dgrnd_snow)
692
693	!! 0. Variable and parameter declaration
694	!! 0.1 Input variables
695	INTEGER(i_std), INTENT(in) :: kjit !! Time step number (unitless)
696	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (unitless)
697	INTEGER(i_std),INTENT (in) :: rest_id !! Restart file identifier(unitless)
698	REAL(r_std),DIMENSION (kjpindex),INTENT(in) :: gtemp !! First soil layer temperature
699	REAL(r_std),DIMENSION (kjpindex),INTENT(in) :: soilcap
700	REAL(r_std),DIMENSION (kjpindex),INTENT(in) :: soilflx
701	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: lambda_snow !! Coefficient of the linear extrapolation of surface temperature
702	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT (in) :: cgrnd_snow !! Integration coefficient for snow numerical scheme
703	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT (in) :: dgrnd_snow !! Integration coefficient for snow numerical scheme
704
705	!_ ================================================================================================================================
706
707	!! 1. Write variables to restart file to be used for the next simulation
708	IF (printlev>=3) WRITE (numout,*) 'Write restart file with THERMOSOIL variables'
709
710	CALL restput_p(rest_id, 'ptn', nbp_glo, ngrnd, 1, kjit, ptn, 'scatter', nbp_glo, index_g)
711
712	IF (ok_Ecorr) THEN
713	CALL restput_p(rest_id, 'e_soil_lat', nbp_glo, 1 , 1, kjit, e_soil_lat, 'scatter', nbp_glo, index_g)
714	END IF
715
716	CALL restput_p(rest_id, 'gtemp', nbp_glo, 1, 1, kjit, gtemp, 'scatter', nbp_glo, index_g)
717	CALL restput_p(rest_id, 'soilcap', nbp_glo, 1, 1, kjit, soilcap, 'scatter', nbp_glo, index_g)
718	CALL restput_p(rest_id, 'soilflx', nbp_glo, 1, 1, kjit, soilflx, 'scatter', nbp_glo, index_g)
719	CALL restput_p(rest_id, 'cgrnd', nbp_glo, ngrnd-1, 1, kjit, cgrnd, 'scatter', nbp_glo, index_g)
720	CALL restput_p(rest_id, 'dgrnd', nbp_glo, ngrnd-1, 1, kjit, dgrnd, 'scatter', nbp_glo, index_g)
721	CALL restput_p(rest_id, 'cgrnd_snow', nbp_glo, nsnow, 1, kjit, cgrnd_snow, 'scatter', nbp_glo, index_g)
722	CALL restput_p(rest_id, 'dgrnd_snow', nbp_glo, nsnow, 1, kjit, dgrnd_snow, 'scatter', nbp_glo, index_g)
723	CALL restput_p(rest_id, 'lambda_snow', nbp_glo, 1, 1, kjit, lambda_snow, 'scatter', nbp_glo, index_g)
724
725	END SUBROUTINE thermosoilc_finalize
726
727
728	!! ================================================================================================================================
729	!! SUBROUTINE : thermosoilc_clear
730	!!
731	!>\BRIEF Desallocates the allocated arrays.
732	!! The call of thermosoilc_clear originates from sechiba_clear but the calling sequence and
733	!! its purpose require further investigation.
734	!!
735	!! DESCRIPTION : None
736	!!
737	!! RECENT CHANGE(S) : None
738	!!
739	!! MAIN OUTPUT VARIABLE(S): None
740	!!
741	!! REFERENCE(S) : None
742	!!
743	!! FLOWCHART : None
744	!! \n
745	!_ ================================================================================================================================
746
747	SUBROUTINE thermosoilc_clear()
748
749	IF ( ALLOCATED (ptn)) DEALLOCATE (ptn)
750	IF ( ALLOCATED (cgrnd)) DEALLOCATE (cgrnd)
751	IF ( ALLOCATED (dgrnd)) DEALLOCATE (dgrnd)
752	IF ( ALLOCATED (pcapa)) DEALLOCATE (pcapa)
753	IF ( ALLOCATED (pkappa)) DEALLOCATE (pkappa)
754	IF ( ALLOCATED (pcapa_en)) DEALLOCATE (pcapa_en)
755	IF ( ALLOCATED (ptn_beg)) DEALLOCATE (ptn_beg)
756	IF ( ALLOCATED (temp_sol_beg)) DEALLOCATE (temp_sol_beg)
757	IF ( ALLOCATED (surfheat_incr)) DEALLOCATE (surfheat_incr)
758	IF ( ALLOCATED (coldcont_incr)) DEALLOCATE (coldcont_incr)
759	IF ( ALLOCATED (shum_ngrnd_perma)) DEALLOCATE (shum_ngrnd_perma)
760	IF ( ALLOCATED (profil_froz)) DEALLOCATE (profil_froz)
761
762	END SUBROUTINE thermosoilc_clear
763
764
765	!! ================================================================================================================================
766	!! FUNCTION : fz
767	!!
768	!>\BRIEF fz(rk), the function's result, is the "rk"th element of a geometric series
769	!! with first element fz1 and ration zalph.
770	!!
771	!! DESCRIPTION : This function is used to calculate the depths of the boudaries of the thermal layers (zz_coef) and
772	!! of the numerical nodes (zz) of the thermal scheme. Formulae to get the adimensional depths are followings :
773	!! zz(jg) = fz(REAL(jg,r_std) - undemi); \n
774	!! zz_coef(jg) = fz(REAL(jg,r_std))
775	!!
776	!! RECENT CHANGE(S) : None
777	!!
778	!! RETURN VALUE : fz(rk)
779	!!
780	!! REFERENCE(S) : None
781	!!
782	!! FLOWCHART : None
783	!! \n
784	!_ ================================================================================================================================
785
786	FUNCTION fz(rk) RESULT (fz_result)
787
788	!! 0. Variables and parameter declaration
789
790	!! 0.1 Input variables
791
792	REAL(r_std), INTENT(in) :: rk
793
794	!! 0.2 Output variables
795
796	REAL(r_std) :: fz_result
797
798	!! 0.3 Modified variables
799
800	!! 0.4 Local variables
801
802	!_ ================================================================================================================================
803
804	fz_result = fz1 * (zalph ** rk - un) / (zalph - un)
805
806	END FUNCTION fz
807
808
809	!! ================================================================================================================================
810	!! FUNCTION : thermosoilc_levels
811	!!
812	!>\BRIEF Depth of nodes for the thermal layers in meters.
813	!!
814	!! DESCRIPTION : Calculate and return the depth in meters of the nodes of the soil layers. This calculation is the same
815	!! as done in thermosoilc_var_init for zz. See thermosoilc_var_init for more details.
816	!!
817	!! RECENT CHANGE(S) : None
818	!!
819	!! RETURN VALUE : Vector of soil depth for the nodes in meters
820	!!
821	!! REFERENCE(S) : None
822	!!
823	!! FLOWCHART : None
824	!! \n
825	!_ ================================================================================================================================
826
827	FUNCTION thermosoilc_levels() RESULT (zz_out)
828
829	!! 0.1 Return variable
830
831	REAL(r_std), DIMENSION (ngrnd) :: zz_out !! Depth of soil layers in meters
832
833	!! 0.2 Local variables
834	INTEGER(i_std) :: jg
835	REAL(r_std) :: so_capa
836	REAL(r_std) :: so_cond
837
838	!_ ================================================================================================================================
839
840	!! 1. Define some parameters
841	so_capa = (so_capa_dry + so_capa_wet)/deux
842	so_cond = (so_cond_dry + so_cond_wet)/deux
843
844	cstgrnd=SQRT(one_day / pi)
845	lskin = SQRT(so_cond / so_capa * one_day / pi)
846
847	!! Parameters needed by fz function
848	fz1 = 0.3_r_std * cstgrnd
849	zalph = deux
850
851	!! 2. Get adimentional depth of the numerical nodes
852	DO jg=1,ngrnd
853	zz_out(jg) = fz(REAL(jg,r_std) - undemi)
854	ENDDO
855
856	!! 3. Convert to meters
857	DO jg=1,ngrnd
858	zz_out(jg) = zz_out(jg) / cstgrnd * lskin
859	END DO
860
861	END FUNCTION thermosoilc_levels
862
863
864	!! ================================================================================================================================
865	!! SUBROUTINE : thermosoilc_var_init
866	!!
867	!>\BRIEF Define and initializes the soil thermal parameters
868	!!
869	!! DESCRIPTION : This routine\n
870	!! 1. Defines the parameters ruling the vertical grid of the thermal scheme (fz1, zalpha).\n
871	!! 2. Defines the scaling coefficients for adimensional depths (lskin, cstgrnd, see explanation in the
872	!! variables description of thermosoilc_main). \n
873	!! 3. Calculates the vertical discretization of the soil (zz, zz_coef, dz2) and the constants used
874	!! in the numerical scheme and which depend only on the discretization (dz1, lambda).\n
875	!! 4. Initializes the soil thermal parameters (capacity, conductivity) based on initial soil moisture content.\n
876	!! 5. Produces a first temperature diagnostic based on temperature initialization.\n
877	!!
878	!! The scheme comprizes ngrnd=7 layers by default.
879	!! The layer' s boundaries depths (zz_coef) follow a geometric series of ratio zalph=2 and first term fz1.\n
880	!! zz_coef(jg)=fz1.(1-zalph^jg)/(1-zalph) \n
881	!! The layers' boudaries depths are first calculated 'adimensionally', ie with a
882	!! discretization adapted to EQ5. This discretization is chosen for its ability at
883	!! reproducing a thermal signal with periods ranging from days to centuries. (see
884	!! Hourdin, 1992). Typically, fz1 is chosen as : fz1=0.3*cstgrnd (with cstgrnd the
885	!! adimensional attenuation depth). \n
886	!! The factor lskin/cstgrnd is then used to go from adimensional depths to
887	!! depths in m.\n
888	!! zz(real)=lskin/cstgrnd*zz(adimensional)\n
889	!! Similarly, the depths of the numerical nodes are first calculated
890	!! adimensionally, then the conversion factor is applied.\n
891	!! the numerical nodes (zz) are not exactly the layers' centers : their depths are calculated as follows:\n
892	!! zz(jg)=fz1.(1-zalph^(jg-1/2))/(1-zalph)\n
893	!! The values of zz and zz_coef used in the default thermal discretization are in the following table.
894	!! \latexonly
895	!! \includegraphics{thermosoilc_var_init1.jpg}
896	!! \endlatexonly\n
897	!!
898	!! RECENT CHANGE(S) : None
899	!!
900	!! MAIN OUTPUT VARIABLE(S) : None
901	!!
902	!! REFERENCE(S) :
903	!! - Hourdin, F. (1992). Study and numerical simulation of the general circulation of
904	!! planetary atmospheres, Ph.D. thesis, Paris VII University.
905	!!
906	!! FLOWCHART : None
907	!! \n
908	!_ ================================================================================================================================
909
910	SUBROUTINE thermosoilc_var_init(kjpindex, shumdiag_perma, snow, &
911	zz, zz_coef, dz1, dz2, stempdiag)
912
913	!! 0. Variables and parameter declaration
914
915	!! 0.1 Input variables
916
917	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (unitless)
918	REAL(r_std), DIMENSION (kjpindex,nslm), INTENT (in) :: shumdiag_perma !! Relative soil humidity on the diagnostic axis
919	!! (unitless), [0,1]. (see description of the
920	!! variables of thermosoilc_main for more
921	!! explanations)
922	REAL(r_std), DIMENSION (kjpindex), INTENT (in) :: snow !! Snow quantity
923
924	!! 0.2 Output variables
925
926	REAL(r_std), DIMENSION (ngrnd), INTENT(out) :: zz !! depths of the layers'numerical nodes
927	!! @tex ($m$)@endtex
928	REAL(r_std), DIMENSION (ngrnd), INTENT(out) :: zz_coef !! depths of the layers'boundaries
929	!! @tex ($m$)@endtex
930	REAL(r_std), DIMENSION (ngrnd), INTENT(out) :: dz1 !! numerical constant depending on the vertical
931	!! thermal grid only @tex ($m^{-1}$) @endtex.
932	!! (see description
933	!! of the variables of thermosoilc_main for more
934	!! explanations)
935	REAL(r_std), DIMENSION (ngrnd), INTENT(out) :: dz2 !! thicknesses of the soil thermal layers
936	!! @tex ($m$) @endtex
937	REAL(r_std), DIMENSION (kjpindex,nslm), INTENT (out) :: stempdiag !! Diagnostic temperature profile @tex ($K$)
938	!! @endtex
939
940	!! 0.3 Modified variables
941
942	!! 0.4 Local variables
943
944	INTEGER(i_std) :: ier, ji, jg !! Index (unitless)
945	REAL(r_std) :: sum !! Temporary variable
946	REAL(r_std) :: so_capa !! Average Thermal Conductivity of soils
947	!! @tex $(W.m^{-2}.K^{-1})$ @endtex
948	REAL(r_std) :: so_cond !! Average Thermal Conductivity of soils
949	!! @tex $(W.m^{-2}.K^{-1})$ @endtex
950
951	!_ ================================================================================================================================
952
953	!! 1. Initialization of the parameters of the vertical discretization and of the attenuation depths
954
955	!! so_capa and so_cond are temporary variables which contain the average values of soil conductivity
956	!! and soil conductivity and are only needed in thermosoilc_var_init to set the vertical layering.
957	so_capa = (so_capa_dry + so_capa_wet)/deux
958	so_cond = (so_cond_dry + so_cond_wet)/deux
959
960	cstgrnd=SQRT(one_day / pi)
961	lskin = SQRT(so_cond / so_capa * one_day / pi)
962	fz1 = 0.3_r_std * cstgrnd
963	zalph = deux
964
965	!! Calculate so_capa_ice
966	so_capa_ice = so_capa_dry + poroscapa_icerho_ice
967	WRITE(numout,*) 'Calculation of so_capa_ice=', so_capa_ice,' using poros=',poros,' and capa_ice=',capa_ice
968
969	!! 2. Computing the depth of the thermal levels (numerical nodes) and the layers boundaries
970
971	!! Computing the depth of the thermal levels (numerical nodes) and
972	!! the layers boundariesusing the so-called
973	!! adimentional variable z' = z/lskin*cstgrnd (with z in m)
974
975	!! 2.1 adimensional thicknesses of the layers
976	DO jg=1,ngrnd
977
978	!!?? code simplification hopefully possible here with up-to-date compilers !
979	!!! This needs to be solved soon. Either we allow CPP options in SECHIBA or the VPP
980	!!! fixes its compiler
981	!!!#ifdef VPP5000
982	dz2(jg) = fz(REAL(jg,r_std)-undemi+undemi) - fz(REAL(jg-1,r_std)-undemi+undemi)
983	!!!#else
984	!!! dz2(jg) = fz(REAL(jg,r_std)) - fz(REAL(jg-1,r_std))
985	!!!#endif
986	ENDDO
987
988	!! 2.2 adimentional depth of the numerical nodes and layers' boudaries
989	DO jg=1,ngrnd
990	zz(jg) = fz(REAL(jg,r_std) - undemi)
991	zz_coef(jg) = fz(REAL(jg,r_std)-undemi+undemi)
992	ENDDO
993
994	!! 2.3 Converting to meters
995	DO jg=1,ngrnd
996	zz(jg) = zz(jg) / cstgrnd * lskin
997	zz_coef(jg) = zz_coef(jg) / cstgrnd * lskin
998	dz2(jg) = dz2(jg) / cstgrnd * lskin
999	ENDDO
1000
1001	!! 2.4 Computing some usefull constants for the numerical scheme
1002	DO jg=1,ngrnd-1
1003	dz1(jg) = un / (zz(jg+1) - zz(jg))
1004	ENDDO
1005	lambda = zz(1) * dz1(1)
1006
1007	!! 2.5 Get the wetness profile on the thermal vertical grid from the diagnostic axis
1008	CALL thermosoilc_humlev(kjpindex, shumdiag_perma, snow)
1009
1010	!! 3. Diagnostics : consistency checks on the vertical grid.
1011
1012	WRITE (numout,*) 'diagnostic des niveaux dans le sol' !!?? to be changed,
1013	WRITE (numout,*) 'niveaux intermediaires et pleins'
1014	sum = zero
1015	DO jg=1,ngrnd
1016	sum = sum + dz2(jg)
1017	WRITE (numout,*) zz(jg),sum
1018	ENDDO
1019
1020
1021	!! 4. Compute a first diagnostic temperature profile
1022
1023	CALL thermosoilc_diaglev(kjpindex, stempdiag)
1024
1025	IF (printlev>=3) WRITE (numout,*) ' thermosoilc_var_init done '
1026
1027	END SUBROUTINE thermosoilc_var_init
1028
1029
1030	!! ================================================================================================================================
1031	!! SUBROUTINE : thermosoilc_coef
1032	!!
1033	!>\BRIEF Calculate soil thermal properties, integration coefficients, apparent heat flux,
1034	!! surface heat capacity,
1035	!!
1036	!! DESCRIPTION : This routine computes : \n
1037	!! 1. the soil thermal properties. \n
1038	!! 2. the integration coefficients of the thermal numerical scheme, cgrnd and dgrnd,
1039	!! which depend on the vertical grid and on soil properties, and are used at the next
1040	!! timestep.\n
1041	!! 3. the soil apparent heat flux and surface heat capacity soilflux
1042	!! and soilcap, used by enerbil to compute the surface temperature at the next
1043	!! timestep.\n
1044	!! - The soil thermal properties depend on water content (shum_ngrnd_perma, shumdiag_perma) and on the presence
1045	!! of snow : snow is integrated into the soil for the thermal calculations, ie if there
1046	!! is snow on the ground, the first thermal layer(s) consist in snow, depending on the
1047	!! snow-depth. If a layer consists out of snow and soil, wheighed soil properties are
1048	!! calculated\n
1049	!! - The coefficients cgrnd and dgrnd are the integration
1050	!! coefficients for the thermal scheme \n
1051	!! T(k+1)=cgrnd(k)+dgrnd(k)*T(k) \n
1052	!! -- EQ1 -- \n
1053	!! They correspond respectively to $\beta$ and $\alpha$ from F. Hourdin\'s thesis and
1054	!! their expression can be found in this document (eq A19 and A20)
1055	!! - soilcap and soilflux are the apparent surface heat capacity and flux
1056	!! used in enerbil at the next timestep to solve the surface
1057	!! balance for Ts (EQ3); they correspond to $C_s$ and $F_s$ in F.
1058	!! Hourdin\'s PhD thesis and are expressed in eq. A30 and A31. \n
1059	!! soilcap*(Ts(t)-Ts(t-1))/dt=soilflux+otherfluxes(Ts(t)) \n
1060	!! -- EQ3 --\n
1061	!!
1062	!! RECENT CHANGE(S) : None
1063	!!
1064	!! MAIN OUTPUT VARIABLE(S): cgrnd, dgrnd, pcapa, pkappa, soilcap, soilflx
1065	!!
1066	!! REFERENCE(S) :
1067	!! - Hourdin, F. (1992). Study and numerical simulation of the general circulation of planetary atmospheres,
1068	!! Ph.D. thesis, Paris VII University. Remark: the part of F. Hourdin's PhD thesis relative to the thermal
1069	!! integration scheme has been scanned and is provided along with the documentation, with name :
1070	!! Hourdin_1992_PhD_thermal_scheme.pdf
1071	!!
1072	!! FLOWCHART : None
1073	!! \n
1074	!_ ================================================================================================================================
1075
1076	SUBROUTINE thermosoilc_coef (kjpindex, temp_sol_new, snow, &
1077	ptn, soilcap, soilflx, zz, &
1078	dz1, dz2, &
1079	cgrnd, dgrnd, &
1080	frac_snow_veg,frac_snow_nobio,totfrac_nobio, &
1081	snowdz,snowrho, snowtemp, pb, &
1082	lambda_snow, cgrnd_snow, dgrnd_snow)
1083
1084	!! 0. Variables and parameter declaration
1085
1086	!! 0.1 Input variables
1087
1088	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (unitless)
1089	REAL(r_std), DIMENSION (kjpindex), INTENT (in) :: temp_sol_new !! soil surface temperature @tex ($K$) @endtex
1090	REAL(r_std), DIMENSION (kjpindex), INTENT (in) :: snow !! snow mass @tex ($Kg$) @endtex
1091	REAL(r_std), DIMENSION (ngrnd), INTENT(in) :: zz !! depths of the soil thermal numerical nodes
1092	!! @tex ($m$) @endtex
1093	REAL(r_std), DIMENSION (ngrnd), INTENT(in) :: dz1 !! numerical constant depending on the vertical
1094	!! thermal grid only @tex ($m^{-1}$) @endtex
1095	REAL(r_std), DIMENSION (ngrnd), INTENT(in) :: dz2 !! thicknesses of the soil thermal layers
1096	!! @tex ($m$) @endtex
1097	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: frac_snow_veg !! Snow cover fraction on vegeted area
1098	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT(in) :: frac_snow_nobio !! Snow cover fraction on non-vegeted area
1099	REAL(r_std),DIMENSION (kjpindex),INTENT(in) :: totfrac_nobio !! Total fraction of continental ice+lakes+cities+...
1100	!!(unitless,0-1)
1101	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(in) :: snowdz !! Snow depth
1102	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(in) :: snowrho !! Snow density
1103	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(in) :: snowtemp !! Snow temperature
1104	REAL(r_std), DIMENSION (kjpindex), INTENT (in) :: pb !! Surface presure (hPa)
1105
1106	!! 0.2 Output variables
1107
1108	REAL(r_std), DIMENSION (kjpindex), INTENT (out) :: soilcap !! surface heat capacity
1109	!! @tex ($J m^{-2} K^{-1}$) @endtex
1110	REAL(r_std), DIMENSION (kjpindex), INTENT (out) :: soilflx !! surface heat flux @tex ($W m^{-2}$) @endtex,
1111	!! positive towards the
1112	!! soil, writen as Qg (ground heat flux) in the history
1113	!! files.
1114	REAL(r_std), DIMENSION (kjpindex,ngrnd-1), INTENT(out) :: cgrnd !! matrix coefficient for the computation of soil
1115	!! temperatures (beta in F. Hourdin thesis)
1116	REAL(r_std), DIMENSION (kjpindex,ngrnd-1), INTENT(out) :: dgrnd !! matrix coefficient for the computation of soil
1117	!! temperatures (alpha in F. Hourdin thesis)
1118
1119
1120	!! 0.3 Modified variable
1121
1122	REAL(r_std), DIMENSION (kjpindex,ngrnd), INTENT (inout):: ptn !! vertically discretized soil temperatures
1123	!! @tex ($K$) @endtex
1124	REAL(r_std), DIMENSION (kjpindex), INTENT(inout) :: lambda_snow !! Coefficient of the linear extrapolation of surface temperature
1125	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT (inout):: cgrnd_snow !! Integration coefficient for snow numerical scheme
1126	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT (inout):: dgrnd_snow !! Integration coefficient for snow numerical scheme
1127
1128	!! 0.4 Local variables
1129
1130	INTEGER(i_std) :: ji, jg
1131	REAL(r_std), DIMENSION (kjpindex,ngrnd-1) :: zdz1 !! numerical (buffer) constant
1132	!! @tex ($W m^{-1} K^{-1}$) @endtex
1133	REAL(r_std), DIMENSION (kjpindex,ngrnd) :: zdz2 !! numerical (buffer) constant
1134	!! @tex ($W m^{-1} K^{-1}$) @endtex
1135	REAL(r_std), DIMENSION (kjpindex) :: z1 !! numerical constant @tex ($W m^{-1} K^{-1}$) @endtex
1136	REAL(r_std), DIMENSION (kjpindex) :: soilcap_nosnow !! surface heat capacity
1137	!! @tex ($J m^{-2} K^{-1}$) @endtex
1138	REAL(r_std), DIMENSION (kjpindex) :: soilflx_nosnow !! surface heat flux @tex ($W m^{-2}$) @endtex,
1139	!! positive towards the soil, written as Qg (ground heat flux in the history files).
1140	REAL(r_std), DIMENSION (kjpindex) :: snowcap !! apparent snow heat capacity @tex ($J m^{-2} K^{-1}$)
1141	REAL(r_std), DIMENSION (kjpindex) :: snowflx !! apparent snow-atmosphere heat flux @tex ($W m^{-2}$) @endtex
1142	REAL(r_std), DIMENSION (kjpindex,nsnow) :: dz1_snow
1143	REAL(r_std), DIMENSION (kjpindex,nsnow) :: ZSNOWDZM
1144	REAL(r_std), DIMENSION (kjpindex,nsnow) :: dz2_snow
1145	REAL(r_std), DIMENSION (kjpindex,nsnow) :: zdz1_snow
1146	REAL(r_std), DIMENSION (kjpindex,nsnow) :: zdz2_snow
1147	REAL(r_std), DIMENSION (kjpindex) :: z1_snow
1148
1149	!_ ================================================================================================================================
1150
1151	!! 1. Computation of the soil thermal properties
1152
1153	! Computation of the soil thermal properties; snow properties are also accounted for
1154	IF (ok_freeze_thermix) THEN
1155	CALL thermosoilc_getdiff( kjpindex, snow, ptn, snowrho, snowtemp, pb )
1156	ELSE IF (ok_explicitsnow) THEN
1157	! Special case with explicit snow model without soil freezing
1158	CALL thermosoilc_getdiff_old_thermix_without_snow( kjpindex, snowrho, snowtemp, pb )
1159	ELSE
1160	! Special case with old snow without soil freezing
1161	CALL thermosoilc_getdiff_old_thermix_with_snow( kjpindex, snow )
1162	ENDIF
1163
1164	! Energy conservation : Correction to make sure that the same latent heat is released and
1165	! consumed during freezing and thawing
1166	IF (ok_Ecorr) THEN
1167	CALL thermosoilc_readjust(kjpindex, ptn)
1168	ENDIF
1169
1170
1171	!! 2. Computation of the coefficients of the numerical integration scheme for the soil layers
1172
1173	!! 2.1 Calculate numerical coefficients zdz1 and zdz2
1174	DO jg=1,ngrnd
1175	DO ji=1,kjpindex
1176	zdz2(ji,jg)=pcapa(ji,jg) * dz2(jg)/dt_sechiba
1177	ENDDO
1178	ENDDO
1179
1180	DO jg=1,ngrnd-1
1181	DO ji=1,kjpindex
1182	zdz1(ji,jg) = dz1(jg) * pkappa(ji,jg)
1183	ENDDO
1184	ENDDO
1185
1186	!! 2.2 Calculate coefficients cgrnd and dgrnd for soil
1187	DO ji = 1,kjpindex
1188	z1(ji) = zdz2(ji,ngrnd) + zdz1(ji,ngrnd-1)
1189	cgrnd(ji,ngrnd-1) = zdz2(ji,ngrnd) * ptn(ji,ngrnd) / z1(ji)
1190	dgrnd(ji,ngrnd-1) = zdz1(ji,ngrnd-1) / z1(ji)
1191	ENDDO
1192
1193	DO jg = ngrnd-1,2,-1
1194	DO ji = 1,kjpindex
1195	z1(ji) = un / (zdz2(ji,jg) + zdz1(ji,jg-1) + zdz1(ji,jg) * (un - dgrnd(ji,jg)))
1196	cgrnd(ji,jg-1) = (ptn(ji,jg) * zdz2(ji,jg) + zdz1(ji,jg) * cgrnd(ji,jg)) * z1(ji)
1197	dgrnd(ji,jg-1) = zdz1(ji,jg-1) * z1(ji)
1198	ENDDO
1199	ENDDO
1200
1201
1202	!! 3. Computation of the coefficients of the numerical integration scheme for the snow layers
1203
1204	!! 3.1 Calculate numerical coefficients zdz1_snow, zdz2_snow and lambda_snow
1205	DO ji = 1, kjpindex
1206
1207	dz2_snow(ji,:) = 0
1208
1209	IF ( ok_explicitsnow ) THEN
1210
1211	! Calculate internal values
1212	DO jg = 1, nsnow
1213	ZSNOWDZM(ji,jg) = MAX(snowdz(ji,jg),psnowdzmin)
1214	ENDDO
1215	dz2_snow(ji,:)=ZSNOWDZM(ji,:)
1216
1217	DO jg = 1, nsnow-1
1218	dz1_snow(ji,jg) = 2.0 / (dz2_snow(ji,jg+1)+dz2_snow(ji,jg))
1219	ENDDO
1220
1221	lambda_snow(ji) = dz2_snow(ji,1)/2.0 * dz1_snow(ji,1)
1222
1223	DO jg=1,nsnow
1224	zdz2_snow(ji,jg)=pcapa_snow(ji,jg) * dz2_snow(ji,jg)/dt_sechiba
1225	ENDDO
1226
1227	DO jg=1,nsnow-1
1228	zdz1_snow(ji,jg) = dz1_snow(ji,jg) * pkappa_snow(ji,jg)
1229	ENDDO
1230
1231	! the bottom snow
1232	zdz1_snow(ji,nsnow) = pkappa_snow(ji,nsnow) / ( zz_coef(1) + dz2_snow(ji,nsnow)/2 )
1233
1234	ELSE
1235	! There is no snow, these values will never be used
1236	lambda_snow(ji) = lambda
1237	ENDIF
1238
1239	ENDDO
1240
1241	!! 3.2 Calculate coefficients cgrnd_snow and dgrnd_snow for snow
1242	DO ji = 1,kjpindex
1243	IF ( ok_explicitsnow ) THEN
1244	! bottom level
1245	z1_snow(ji) = zdz2(ji,1)+(un-dgrnd(ji,1))*zdz1(ji,1)+zdz1_snow(ji,nsnow)
1246	cgrnd_snow(ji,nsnow) = (zdz2(ji,1) * ptn(ji,1) + zdz1(ji,1) * cgrnd(ji,1) ) / z1_snow(ji)
1247	dgrnd_snow(ji,nsnow) = zdz1_snow(ji,nsnow) / z1_snow(ji)
1248
1249	! next-to-bottom level
1250	z1_snow(ji) = zdz2_snow(ji,nsnow)+(un-dgrnd_snow(ji,nsnow))*zdz1_snow(ji,nsnow)+zdz1_snow(ji,nsnow-1)
1251	cgrnd_snow(ji,nsnow-1) = (zdz2_snow(ji,nsnow)*snowtemp(ji,nsnow)+&
1252	zdz1_snow(ji,nsnow)*cgrnd_snow(ji,nsnow))/z1_snow(ji)
1253	dgrnd_snow(ji,nsnow-1) = zdz1_snow(ji,nsnow-1) / z1_snow(ji)
1254
1255	DO jg = nsnow-1,2,-1
1256	z1_snow(ji) = un / (zdz2_snow(ji,jg) + zdz1_snow(ji,jg-1) + zdz1_snow(ji,jg) * (un - dgrnd_snow(ji,jg)))
1257	cgrnd_snow(ji,jg-1) = (snowtemp(ji,jg) * zdz2_snow(ji,jg) + zdz1_snow(ji,jg) * cgrnd_snow(ji,jg)) * z1_snow(ji)
1258	dgrnd_snow(ji,jg-1) = zdz1_snow(ji,jg-1) * z1_snow(ji)
1259	ENDDO
1260	ELSE
1261	! There is no snow, these values will never be used
1262	cgrnd_snow(ji,:) = cgrnd(ji,1)
1263	dgrnd_snow(ji,:) = dgrnd(ji,1)
1264	ENDIF
1265	ENDDO
1266
1267
1268	!! 4. Computation of the apparent ground heat flux
1269	!! Computation of apparent snow-atmosphere flux
1270	DO ji = 1,kjpindex
1271	IF ( ok_explicitsnow ) THEN
1272	snowflx(ji) = zdz1_snow(ji,1) * (cgrnd_snow(ji,1) + (dgrnd_snow(ji,1)-1.) * snowtemp(ji,1))
1273	snowcap(ji) = (zdz2_snow(ji,1) * dt_sechiba + dt_sechiba * (un - dgrnd_snow(ji,1)) * zdz1_snow(ji,1))
1274	z1_snow(ji) = lambda_snow(ji) * (un - dgrnd_snow(ji,1)) + un
1275	snowcap(ji) = snowcap(ji) / z1_snow(ji)
1276	snowflx(ji) = snowflx(ji) + &
1277	& snowcap(ji) * (snowtemp(ji,1) * z1_snow(ji) - lambda_snow(ji) * cgrnd_snow(ji,1) - temp_sol_new(ji)) / dt_sechiba
1278	ELSE
1279	snowflx(ji) = zero
1280	snowcap(ji) = zero
1281	ENDIF
1282	ENDDO
1283
1284
1285	!! Computation of the apparent ground heat flux (> towards the soil) and
1286	!! apparent surface heat capacity, used at the next timestep by enerbil to
1287	!! compute the surface temperature.
1288	DO ji = 1,kjpindex
1289	soilflx_nosnow(ji) = zdz1(ji,1) * (cgrnd(ji,1) + (dgrnd(ji,1)-1.) * ptn(ji,1))
1290	soilcap_nosnow(ji) = (zdz2(ji,1) * dt_sechiba + dt_sechiba * (un - dgrnd(ji,1)) * zdz1(ji,1))
1291	z1(ji) = lambda * (un - dgrnd(ji,1)) + un
1292	soilcap_nosnow(ji) = soilcap_nosnow(ji) / z1(ji)
1293	soilflx_nosnow(ji) = soilflx_nosnow(ji) + &
1294	& soilcap_nosnow(ji) * (ptn(ji,1) * z1(ji) - lambda * cgrnd(ji,1) - temp_sol_new(ji)) / dt_sechiba
1295	ENDDO
1296
1297	!! Add snow fraction
1298	IF ( ok_explicitsnow ) THEN
1299	! Using an effective heat capacity and heat flux by a simple pondering of snow and soil fraction
1300	DO ji = 1, kjpindex
1301	soilcap(ji) = snowcap(ji)frac_snow_veg(ji)(1-totfrac_nobio(ji))+ & ! weights related to snow cover fraction on vegetation
1302	soilcap_nosnow(ji)SUM(frac_snow_nobio(ji,:))totfrac_nobio(ji)+ & ! weights related to SCF on nobio
1303	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
1304	soilflx(ji) = snowflx(ji)frac_snow_veg(ji)(1-totfrac_nobio(ji))+ & ! weights related to snow cover fraction on vegetation
1305	soilflx_nosnow(ji)SUM(frac_snow_nobio(ji,:))totfrac_nobio(ji)+ & ! weights related to SCF on nobio
1306	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
1307	ENDDO
1308	ELSE
1309	! Do not consider snow fraction
1310	soilcap(:)=soilcap_nosnow(:)
1311	soilflx(:)=soilflx_nosnow(:)
1312	END IF
1313
1314	IF (printlev>=3) WRITE (numout,*) ' thermosoilc_coef done '
1315
1316	END SUBROUTINE thermosoilc_coef
1317
1318
1319	!! ================================================================================================================================
1320	!! SUBROUTINE : thermosoilc_profile
1321	!!
1322	!>\BRIEF In this routine solves the numerical soil thermal scheme, ie calculates the new soil temperature profile;
1323	!! This profile is then exported onto the diagnostic axis (call thermosoilc_diaglev)
1324	!!
1325	!! DESCRIPTION : The calculation of the new soil temperature profile is based on
1326	!! the cgrnd and dgrnd values from the previous timestep and the surface temperature Ts aka temp_sol_new. (see detailed
1327	!! explanation in the header of the thermosoilc module or in the reference).\n
1328	!! T(k+1)=cgrnd(k)+dgrnd(k)*T(k)\n
1329	!! -- EQ1 --\n
1330	!! Ts=(1-lambda)T(1) -lambdaT(2)\n
1331	!! -- EQ2--\n
1332	!!
1333	!! RECENT CHANGE(S) : None
1334	!!
1335	!! MAIN OUTPUT VARIABLE(S): ptn (soil temperature profile on the thermal axis),
1336	!! stempdiag (soil temperature profile on the diagnostic axis)
1337	!!
1338	!! REFERENCE(S) :
1339	!! - Hourdin, F. (1992). Study and numerical simulation of the general circulation of planetary atmospheres,
1340	!! Ph.D. thesis, Paris VII University. Remark: the part of F. Hourdin's PhD thesis relative to the thermal
1341	!! integration scheme has been scanned and is provided along with the documentation, with name :
1342	!! Hourdin_1992_PhD_thermal_scheme.pdf
1343	!!
1344	!! FLOWCHART : None
1345	!! \n
1346	!_ ================================================================================================================================
1347
1348	SUBROUTINE thermosoilc_profile (kjpindex, temp_sol_new, ptn, stempdiag, &
1349	snowtemp, frac_snow_veg,frac_snow_nobio,totfrac_nobio, &
1350	cgrnd_snow,dgrnd_snow)
1351
1352	!! 0. Variables and parameter declaration
1353
1354	!! 0.1 Input variables
1355
1356	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (unitless)
1357	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_sol_new !! Surface temperature at the present time-step
1358	!! @tex ($K$) @endtex
1359	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: frac_snow_veg !! Snow cover fraction on vegeted area
1360	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT(in) :: frac_snow_nobio!! Snow cover fraction on non-vegeted area
1361	REAL(r_std),DIMENSION (kjpindex),INTENT(in) :: totfrac_nobio !! Total fraction of continental ice+lakes+cities+...(unitless,0-1)
1362	REAL(r_std),DIMENSION (kjpindex,nsnow), INTENT(in) :: snowtemp !! Snow temperature
1363	REAL(r_std),DIMENSION (kjpindex,nsnow), INTENT(in) :: cgrnd_snow !! Integration coefficient for snow numerical scheme
1364	REAL(r_std),DIMENSION (kjpindex,nsnow), INTENT(in) :: dgrnd_snow !! Integration coefficient for snow numerical scheme
1365
1366
1367	!! 0.2 Output variables
1368	REAL(r_std),DIMENSION (kjpindex,nslm), INTENT (out) :: stempdiag !! diagnostic temperature profile
1369	!! @tex ($K$) @endtex
1370	REAL(r_std),DIMENSION (kjpindex,ngrnd), INTENT (out) :: ptn !! vertically discretized soil temperatures
1371	!! @tex ($K$) @endtex
1372
1373	!! 0.3 Modified variables
1374
1375
1376	!! 0.4 Local variables
1377
1378	INTEGER(i_std) :: ji, jg
1379	REAL(r_std) :: temp_sol_eff
1380
1381	!_ ================================================================================================================================
1382
1383	!! 1. Computes the soil temperatures ptn.
1384
1385	!! 1.1. ptn(jg=1) using EQ1 and EQ2
1386	DO ji = 1,kjpindex
1387
1388	IF ( ok_explicitsnow ) THEN
1389	! Soil temperature calculation with explicit snow if there is snow on the ground
1390	! using an effective surface temperature by a simple pondering
1391	temp_sol_eff=snowtemp(ji,nsnow)frac_snow_veg(ji)(1-totfrac_nobio(ji))+ & ! weights related to snow cover fraction on vegetation
1392	temp_sol_new(ji)SUM(frac_snow_nobio(ji,:))totfrac_nobio(ji)+ & ! weights related to SCF on nobio
1393	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
1394	! Soil temperature calculation with explicit snow if there is snow on the ground
1395	ptn(ji,1) = cgrnd_snow(ji,nsnow) + dgrnd_snow(ji,nsnow) * temp_sol_eff
1396	ELSE
1397	! Standard soil temperature calculation
1398	ptn(ji,1) = (lambda * cgrnd(ji,1) + temp_sol_new(ji)) / (lambda *(un - dgrnd(ji,1)) + un)
1399	ENDIF
1400	ENDDO
1401
1402	!! 1.2. ptn(jg=2:ngrnd) using EQ1.
1403	DO jg = 1,ngrnd-1
1404	DO ji = 1,kjpindex
1405	ptn(ji,jg+1) = cgrnd(ji,jg) + dgrnd(ji,jg) * ptn(ji,jg)
1406	ENDDO
1407	ENDDO
1408
1409	!! 2. Put the soil temperatures onto the diagnostic axis
1410
1411	!! Put the soil temperatures onto the diagnostic axis for convenient
1412	!! use in other routines (stomate..)
1413	CALL thermosoilc_diaglev(kjpindex, stempdiag)
1414
1415	IF (printlev>=3) WRITE (numout,*) ' thermosoilc_profile done '
1416
1417	END SUBROUTINE thermosoilc_profile
1418
1419
1420	!! ================================================================================================================================
1421	!! SUBROUTINE : thermosoilc_diaglev
1422	!!
1423	!>\BRIEF Interpolation of the soil in-depth temperatures onto the diagnostic profile.
1424	!!
1425	!! DESCRIPTION : This is a very easy linear interpolation, with intfact(jsl, jg) the fraction
1426	!! the thermal layer jg comprised within the diagnostic layer jsl. The depths of
1427	!! the diagnostic levels are diaglev(1:nslm), computed in slowproc.f90.
1428	!!
1429	!! RECENT CHANGE(S) : None
1430	!!
1431	!! MAIN OUTPUT VARIABLE(S): stempdiag (soil temperature profile on the diagnostic axis)
1432	!!
1433	!! REFERENCE(S) : None
1434	!!
1435	!! FLOWCHART : None
1436	!! \n
1437	!_ ================================================================================================================================
1438
1439	SUBROUTINE thermosoilc_diaglev(kjpindex, stempdiag)
1440
1441	!! 0. Variables and parameter declaration
1442
1443	!! 0.1 Input variables
1444
1445	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (unitless)
1446
1447	!! 0.2 Output variables
1448
1449	REAL(r_std),DIMENSION (kjpindex,nslm), INTENT (out) :: stempdiag !! Diagnostoc soil temperature profile @tex ($K$) @endtex
1450
1451	!! 0.3 Modified variables
1452
1453	!! 0.4 Local variables
1454
1455	INTEGER(i_std) :: ji, jsl, jg
1456	REAL(r_std) :: lev_diag, prev_diag, lev_prog, prev_prog
1457	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: intfact
1458	!$OMP THREADPRIVATE(intfact)
1459	LOGICAL, PARAMETER :: check=.FALSE.
1460	!_ ================================================================================================================================
1461
1462	!! 1. Computes intfact(jsl, jg)
1463
1464	!! Computes intfact(jsl, jg), the fraction
1465	!! the thermal layer jg comprised within the diagnostic layer jsl.
1466	IF ( .NOT. ALLOCATED(intfact)) THEN
1467
1468	ALLOCATE(intfact(nslm, ngrnd))
1469
1470	prev_diag = zero
1471	DO jsl = 1, nslm
1472	lev_diag = diaglev(jsl)
1473	prev_prog = zero
1474	DO jg = 1, ngrnd
1475	IF ( jg == ngrnd .AND. (prev_prog + dz2(jg)) < lev_diag ) THEN
1476	lev_prog = lev_diag
1477	ELSE
1478	lev_prog = prev_prog + dz2(jg)
1479	ENDIF
1480	intfact(jsl,jg) = MAX(MIN(lev_diag,lev_prog)-MAX(prev_diag, prev_prog), zero)/(lev_diag-prev_diag)
1481	prev_prog = lev_prog
1482	ENDDO
1483	prev_diag = lev_diag
1484	ENDDO
1485
1486	IF ( check ) THEN
1487	WRITE(numout,*) 'thermosoilc_diagev -- thermosoilc_diaglev -- thermosoilc_diaglev --'
1488	DO jsl = 1, nslm
1489	WRITE(numout,*) jsl, '-', intfact(jsl,1:ngrnd)
1490	ENDDO
1491	WRITE(numout,*) "SUM -- SUM -- SUM SUM -- SUM -- SUM"
1492	DO jsl = 1, nslm
1493	WRITE(numout,*) jsl, '-', SUM(intfact(jsl,1:ngrnd))
1494	ENDDO
1495	WRITE(numout,*) 'thermosoilc_diaglev -- thermosoilc_diaglev -- thermosoilc_diaglev --'
1496	ENDIF
1497
1498	ENDIF
1499
1500	!! 2. does the interpolation
1501
1502	stempdiag(:,:) = zero
1503	DO jg = 1, ngrnd
1504	DO jsl = 1, nslm
1505	DO ji = 1, kjpindex
1506	stempdiag(ji,jsl) = stempdiag(ji,jsl) + ptn(ji,jg)*intfact(jsl,jg)
1507	ENDDO
1508	ENDDO
1509	ENDDO
1510
1511	END SUBROUTINE thermosoilc_diaglev
1512
1513
1514	!! ================================================================================================================================
1515	!! SUBROUTINE : thermosoilc_humlev
1516	!!
1517	!>\BRIEF Interpolates the diagnostic soil humidity profile shumdiag_perma(nslm, diagnostic axis) onto
1518	!! the thermal axis, which gives shum_ngrnd_perma(ngrnd, thermal axis).
1519	!!
1520	!! DESCRIPTION : Same as in thermosoilc_diaglev : This is a very easy linear interpolation, with intfactw(jsl, jg) the fraction
1521	!! the thermal layer jsl comprised within the diagnostic layer jg.
1522	!!
1523	!! The depths of the diagnostic levels are diaglev(1:nslm), computed in slowproc.f90.
1524	!!
1525	!! RECENT CHANGE(S) : None
1526	!!
1527	!! MAIN OUTPUT VARIABLE(S): shum_ngrnd_perma (soil saturation degree on the thermal axis)
1528	!!
1529	!! REFERENCE(S) : None
1530	!!
1531	!! FLOWCHART : None
1532	!! \n
1533	!_ ================================================================================================================================
1534	SUBROUTINE thermosoilc_humlev(kjpindex, shumdiag_perma, snow)
1535
1536	!! 0. Variables and parameter declaration
1537
1538	!! 0.1 Input variables
1539
1540	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (unitless)
1541	REAL(r_std),DIMENSION (kjpindex,nslm), INTENT (in) :: shumdiag_perma
1542	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: snow
1543
1544	!! 0.2 Output variables
1545
1546	!! 0.3 Modified variables
1547
1548	!! 0.4 Local variables
1549	INTEGER(i_std) :: ji, jsl, jg
1550	REAL(r_std) :: lev_diag, prev_diag, lev_prog, prev_prog
1551	REAL(r_std), DIMENSION(ngrnd,nslm) :: intfactw !! fraction of each diagnostic layer (jsl) comprized within
1552	!! a given thermal layer (jg)(0-1, unitless)
1553	REAL(r_std), DIMENSION(kjpindex) :: snow_h
1554	LOGICAL, PARAMETER :: check=.FALSE.
1555
1556	!_ ================================================================================================================================
1557
1558	!! 1. computes intfactw(jsl,jg), the fraction of each diagnostic layer (jg) comprized within a given thermal layer (jsl)
1559	IF ( check ) &
1560	WRITE(numout,*) 'thermosoilc_humlev --'
1561
1562	! Snow height
1563	snow_h(:)=snow(:)/sn_dens
1564	!
1565	shum_ngrnd_perma(:,:) = zero
1566	DO ji=1,kjpindex
1567	prev_diag = zero
1568	DO jsl = 1, ngrnd
1569	lev_diag = prev_diag + dz2(jsl)
1570	prev_prog = snow_h(ji)
1571	DO jg = 1, nslm
1572	IF ( jg == nslm .AND. diaglev(jg)+snow_h(ji) < lev_diag ) THEN
1573	lev_prog = lev_diag+snow_h(ji)
1574	ELSE
1575	lev_prog = diaglev(jg)+snow_h(ji)
1576	ENDIF
1577	intfactw(jsl,jg) = MAX(MIN(lev_diag,lev_prog)-MAX(prev_diag, prev_prog), zero)/(lev_diag-prev_diag)
1578	prev_prog = lev_prog
1579	ENDDO
1580	prev_diag = lev_diag
1581	ENDDO
1582
1583	DO jsl = 1, nslm
1584	DO jg = 1, ngrnd
1585	shum_ngrnd_perma(ji,jg) = shum_ngrnd_perma(ji,jg) + shumdiag_perma(ji,jsl)*intfactw(jg,jsl)
1586	ENDDO
1587	ENDDO
1588	!
1589	IF ( check ) THEN
1590	DO jg = 1, ngrnd
1591	WRITE(numout,*) ji,jg, '-', intfactw(jg,1:nslm),'-sum-', SUM(intfactw(jg,1:nslm))
1592	ENDDO
1593	ENDIF
1594	ENDDO
1595
1596	IF ( check ) &
1597	WRITE(numout,*) 'thermosoilc_humlev --'
1598
1599	END SUBROUTINE thermosoilc_humlev
1600
1601
1602	!! ================================================================================================================================
1603	!! SUBROUTINE : thermosoilc_energy
1604	!!
1605	!>\BRIEF Energy check-up.
1606	!!
1607	!! DESCRIPTION : I didn\'t comment this routine since at do not understand its use, please
1608	!! ask initial designers (Jan ? Nathalie ?).
1609	!!
1610	!! RECENT CHANGE(S) : None
1611	!!
1612	!! MAIN OUTPUT VARIABLE(S) : ??
1613	!!
1614	!! REFERENCE(S) : None
1615	!!
1616	!! FLOWCHART : None
1617	!! \n
1618	!_ ================================================================================================================================
1619
1620	SUBROUTINE thermosoilc_energy(kjpindex, temp_sol_new, soilcap)
1621
1622	!! 0. Variables and parameter declaration
1623
1624	!! 0.1 Input variables
1625
1626	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (unitless)
1627	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_sol_new !! Surface temperature at the present time-step, Ts
1628	!! @tex ($K$) @endtex
1629	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: soilcap !! Apparent surface heat capacity
1630	!! @tex ($J m^{-2} K^{-1}$) @endtex,
1631	!! see eq. A29 of F. Hourdin\'s PhD thesis.
1632
1633	!! 0.2 Output variables
1634
1635	!! 0.3 Modified variables
1636
1637	!! 0.4 Local variables
1638
1639	INTEGER(i_std) :: ji, jg
1640	!_ ================================================================================================================================
1641
1642
1643	DO ji = 1, kjpindex
1644	surfheat_incr(ji) = zero
1645	coldcont_incr(ji) = zero
1646	ENDDO
1647
1648	! Sum up the energy content of all layers in the soil.
1649	DO ji = 1, kjpindex
1650
1651	IF (pcapa_en(ji,1) .LE. sn_capa) THEN
1652
1653	! Verify the energy conservation in the surface layer
1654	coldcont_incr(ji) = soilcap(ji) * (temp_sol_new(ji) - temp_sol_beg(ji))
1655	surfheat_incr(ji) = zero
1656	ELSE
1657
1658	! Verify the energy conservation in the surface layer
1659	surfheat_incr(ji) = soilcap(ji) * (temp_sol_new(ji) - temp_sol_beg(ji))
1660	coldcont_incr(ji) = zero
1661	ENDIF
1662	ENDDO
1663
1664	ptn_beg(:,:) = ptn(:,:)
1665	temp_sol_beg(:) = temp_sol_new(:)
1666
1667	END SUBROUTINE thermosoilc_energy
1668
1669
1670
1671	!! ================================================================================================================================
1672	!! SUBROUTINE : thermosoilc_readjust
1673	!!
1674	!>\BRIEF
1675	!!
1676	!! DESCRIPTION : Energy conservation : Correction to make sure that the same latent heat is released and
1677	!! consumed during freezing and thawing
1678	!!
1679	!! RECENT CHANGE(S) : None
1680	!!
1681	!! MAIN OUTPUT VARIABLE(S): ptn (soil temperature profile on the thermal axis),
1682	!!
1683	!! REFERENCE(S) :
1684	!!
1685	!! FLOWCHART : None
1686	!! \n
1687	!_ ================================================================================================================================
1688
1689	SUBROUTINE thermosoilc_readjust(kjpindex, ptn)
1690
1691	!! 0. Variables and parameter declaration
1692
1693	!! 0.1 Input variables
1694	INTEGER(i_std), INTENT(in) :: kjpindex
1695
1696	!! 0.2 Modified variables
1697	REAL(r_std),DIMENSION(kjpindex,ngrnd),INTENT(inout) :: ptn
1698
1699	!! 0.3 Local variables
1700	INTEGER(i_std) :: ji, jg
1701	REAL(r_std) :: ptn_tmp
1702
1703	DO jg=1, ngrnd
1704	DO ji=1, kjpindex
1705	! All soil latent energy is put into e_soil_lat(ji)
1706	! because the variable soil layers make it difficult to keep track of all
1707	! layers in this version
1708	! NOTE : pcapa has unit J/K/m3 and pcappa_supp has J/K
1709	e_soil_lat(ji)=e_soil_lat(ji)+pcappa_supp(ji,jg)*(ptn(ji,jg)-ptn_beg(ji,jg))
1710	END DO
1711	END DO
1712
1713	DO ji=1, kjpindex
1714	IF (e_soil_lat(ji).GT.min_sechiba.AND.MINVAL(ptn(ji,:)).GT.ZeroCelsius+fr_dT/2.) THEN
1715	! The soil is thawed: we spread the excess of energy over the uppermost 6 levels e.g. 2.7m
1716	! Here we increase the temperatures
1717	DO jg=1,6
1718	ptn_tmp=ptn(ji,jg)
1719
1720	ptn(ji,jg)=ptn(ji,jg)+MIN(e_soil_lat(ji)/pcapa(ji,jg)/zz_coef(6), 0.5)
1721	e_soil_lat(ji)=e_soil_lat(ji)-(ptn(ji,jg)-ptn_tmp)pcapa(ji,jg)dz2(jg)
1722	ENDDO
1723	ELSE IF (e_soil_lat(ji).LT.-min_sechiba.AND.MINVAL(ptn(ji,:)).GT.ZeroCelsius+fr_dT/2.) THEN
1724	! The soil is thawed
1725	! Here we decrease the temperatures
1726	DO jg=1,6
1727	ptn_tmp=ptn(ji,jg)
1728	ptn(ji,jg)=MAX(ZeroCelsius+fr_dT/2., ptn_tmp+e_soil_lat(ji)/pcapa(ji,jg)/zz_coef(6))
1729	e_soil_lat(ji)=e_soil_lat(ji)+(ptn_tmp-ptn(ji,jg))pcapa(ji,jg)dz2(jg)
1730	END DO
1731	END IF
1732	END DO
1733
1734	END SUBROUTINE thermosoilc_readjust
1735
1736	!-------------------------------------------------------------------
1737
1738
1739
1740	!! ================================================================================================================================
1741	!! SUBROUTINE : thermosoilc_getdiff
1742	!!
1743	!>\BRIEF Computes soil and snow heat capacity and conductivity
1744	!!
1745	!! DESCRIPTION : Computation of the soil and snow thermal properties; snow properties are also accounted for
1746	!!
1747	!! RECENT CHANGE(S) : None
1748	!!
1749	!! MAIN OUTPUT VARIABLE(S):
1750	!!
1751	!! REFERENCE(S) :
1752	!!
1753	!! FLOWCHART : None
1754	!! \n
1755	!_ ================================================================================================================================
1756
1757	SUBROUTINE thermosoilc_getdiff( kjpindex, snow, ptn, snowrho, snowtemp, pb )
1758
1759	!! 0. Variables and parameter declaration
1760
1761	!! 0.1 Input variables
1762	INTEGER(i_std),INTENT(in) :: kjpindex
1763	REAL(r_std),DIMENSION(kjpindex),INTENT (in) :: snow !! Snow mass
1764	REAL(r_std),DIMENSION(kjpindex,nsnow),INTENT(in) :: snowrho !! Snow density
1765	REAL(r_std),DIMENSION(kjpindex,nsnow),INTENT(in) :: snowtemp !! Snow temperature
1766	REAL(r_std),DIMENSION(kjpindex),INTENT(in) :: pb !! Surface pressure
1767	REAL(r_std),DIMENSION(kjpindex,ngrnd),INTENT(in) :: ptn !! Soil temperature profile
1768
1769	!! 0.3 Local variables
1770	REAL :: xx !! Unfrozen fraction of the soil
1771	REAL(r_std), DIMENSION(kjpindex) :: snow_h
1772	REAL(r_std), DIMENSION(kjpindex,ngrnd) :: zx1, zx2
1773	REAL :: cap_iw !! Heat capacity of ice/water mixture
1774	REAL :: csat !! Thermal conductivity for saturated soil
1775	INTEGER :: ji,jg
1776
1777
1778	DO ji = 1,kjpindex
1779	IF (.NOT. ok_explicitsnow) THEN
1780	! 1. Determine the fractions of snow and soil
1781	snow_h(ji) = snow(ji) / sn_dens
1782
1783	!
1784	! 1.1. The first level
1785	!
1786	IF ( snow_h(ji) .GT. zz_coef(1) ) THEN
1787	! the 1st level is in the snow => the 1st layer is entirely snow
1788	zx1(ji,1) = 1.
1789	zx2(ji,1) = 0.
1790	ELSE IF ( snow_h(ji) .GT. zero ) THEN
1791	! the 1st level is beyond the snow and the snow is present
1792	zx1(ji,1) = snow_h(ji) / zz_coef(1)
1793	zx2(ji,1) = ( zz_coef(1) - snow_h(ji)) / zz_coef(1)
1794	ELSE
1795	! there is no snow at all, quoi ;-)
1796	zx1(ji,1) = 0.
1797	zx2(ji,1) = 1.
1798	ENDIF
1799
1800	!
1801	! 1.2. The other levels
1802	!
1803	DO jg = 2, ngrnd
1804	IF ( snow_h(ji) .GT. zz_coef(jg) ) THEN
1805	! the current level is in the snow => the current layer is entirely snow
1806	zx1(ji,jg) = 1.
1807	zx2(ji,jg) = 0.
1808	ELSE IF ( snow_h(ji) .GT. zz_coef(jg-1) ) THEN
1809	! the current layer is partially snow
1810	zx1(ji,jg) = (snow_h(ji) - zz_coef(jg-1)) / (zz_coef(jg) - zz_coef(jg-1))
1811	zx2(ji,jg) = ( zz_coef(jg) - snow_h(ji)) / (zz_coef(jg) - zz_coef(jg-1))
1812	ELSE
1813	! both levels are out of snow => the current layer is entirely soil
1814	zx1(ji,jg) = 0.
1815	zx2(ji,jg) = 1.
1816	ENDIF
1817	ENDDO
1818	ELSE
1819	zx1(ji,:) = 0.
1820	zx2(ji,:) = 1.
1821	END IF
1822
1823	DO jg = 1, ngrnd
1824	!
1825	! 2. Calculate heat capacity with allowance for permafrost
1826	! 2.1. soil heat capacity depending on temperature and humidity
1827
1828	IF (ptn(ji,jg) .LT. ZeroCelsius-fr_dT/2.) THEN
1829	! frozen soil
1830	pcapa(ji,jg) = so_capa_dry + shum_ngrnd_perma(ji,jg)*(so_capa_ice - so_capa_dry)!Isa : old version, proved to be correct
1831	profil_froz(ji,jg) = 1.
1832	pcappa_supp(ji,jg)= 0.
1833
1834	ELSEIF (ptn(ji,jg) .GT. ZeroCelsius+fr_dT/2.) THEN
1835	! unfrozen soil
1836	pcapa(ji,jg) = so_capa_dry + shum_ngrnd_perma(ji,jg)*(so_capa_wet - so_capa_dry)!Isa : old version, proved to be correct
1837	profil_froz(ji,jg) = 0.
1838	pcappa_supp(ji,jg)= 0.
1839	ELSE
1840	! xx is the unfrozen fraction of soil water
1841	xx = (ptn(ji,jg)-(ZeroCelsius-fr_dT/2.)) / fr_dT
1842	profil_froz(ji,jg) = (1. - xx)
1843
1844	! net heat capacity of the ice/water mixture
1845	cap_iw = xx * so_capa_wet + (1.-xx) * so_capa_ice
1846	pcapa(ji,jg) = so_capa_dry + shum_ngrnd_perma(ji,jg)(cap_iw-so_capa_dry) + shum_ngrnd_perma(ji,jg)poroslhfrho_water/fr_dT
1847	pcappa_supp(ji,jg)= shum_ngrnd_perma(ji,jg)poroslhfrho_water/fr_dTzx2(ji,jg)*dz2(jg)
1848
1849	ENDIF
1850
1851	!
1852	! 2.2. Take into account the snow and soil fractions in the layer
1853	!
1854	pcapa(ji,jg) = zx1(ji,jg) * sn_capa + zx2(ji,jg) * pcapa(ji,jg)
1855
1856	!
1857	! 2.3. Calculate the heat capacity for energy conservation check
1858	IF ( zx1(ji,jg).GT.0. ) THEN
1859	pcapa_en(ji,jg) = sn_capa
1860	ELSE
1861	pcapa_en(ji,jg) = pcapa(ji,jg)
1862	ENDIF
1863
1864	!
1865	! 3. Calculate the heat conductivity with allowance for permafrost (Farouki, 1981, Cold Reg. Sci. Technol.)
1866	!
1867	! 3.1. unfrozen fraction
1868	xx = (ptn(ji,jg)-(ZeroCelsius-fr_dT/2.)) / fr_dT * poros
1869	xx = MIN( poros, MAX( 0., xx ) )
1870
1871	! 3.2. saturated conductivity
1872	csat = cond_solid*(1.-poros) cond_ice*(poros-xx) cond_water**xx
1873
1874	! 3.3. unsaturated conductivity
1875	pkappa(ji,jg) =(csat - so_cond_dry)*shum_ngrnd_perma(ji,jg) + so_cond_dry
1876
1877	!
1878	! 3.4. Take into account the snow and soil fractions in the layer
1879	pkappa(ji,jg) = un / ( zx1(ji,jg) / sn_cond + zx2(ji,jg) / pkappa(ji,jg) )
1880	END DO
1881	ENDDO
1882
1883	! 4. Calculate snow heat capacity and conductivity
1884	DO ji = 1,kjpindex
1885	pcapa_snow(ji,:) = snowrho(ji,:) * xci
1886	pkappa_snow(ji,:) = (ZSNOWTHRMCOND1 + ZSNOWTHRMCOND2snowrho(ji,:)snowrho(ji,:)) + &
1887	MAX(0.0,(ZSNOWTHRMCOND_AVAP+(ZSNOWTHRMCOND_BVAP/(snowtemp(ji,:)+ &
1888	ZSNOWTHRMCOND_CVAP)))(XP00/(pb(ji)100.)))
1889	END DO
1890
1891
1892
1893	END SUBROUTINE thermosoilc_getdiff
1894
1895	!! ================================================================================================================================
1896	!! SUBROUTINE : thermosoilc_getdiff_old_thermix_with_snow
1897	!!
1898	!>\BRIEF Computes soil heat capacity and conductivity
1899	!!
1900	!! DESCRIPTION : Computes soil heat capacity and conductivity
1901	!! Special case with old snow without soil freezing
1902	!!
1903	!! RECENT CHANGE(S) : None
1904	!!
1905	!! MAIN OUTPUT VARIABLE(S):
1906	!!
1907	!! REFERENCE(S) :
1908	!!
1909	!! FLOWCHART : None
1910	!! \n
1911	!_ ================================================================================================================================
1912
1913
1914	SUBROUTINE thermosoilc_getdiff_old_thermix_with_snow( kjpindex, snow )
1915
1916
1917	!! 0. Variables and parameter declaration
1918
1919	!! 0.1 Input variables
1920	INTEGER(i_std), INTENT(in) :: kjpindex
1921	REAL(r_std),DIMENSION(kjpindex),INTENT (in) :: snow
1922
1923	!! 0.2 Local variables
1924	INTEGER :: ji,jg
1925	REAL(r_std) :: snow_h !! snow_h is the snow height @tex ($m$) @endtex
1926	REAL(r_std) :: zx1, zx2 !! zx1 and zx2 are the layer fraction consisting in snow and soil respectively.
1927
1928
1929	! Computation of the soil thermal properties; snow properties are also accounted for
1930	DO ji = 1,kjpindex
1931	snow_h = snow(ji) / sn_dens
1932
1933	! First layer
1934	IF ( snow_h .GT. zz_coef(1) ) THEN
1935	pcapa(ji,1) = sn_capa
1936	pcapa_en(ji,1) = sn_capa
1937	pkappa(ji,1) = sn_cond
1938	ELSE IF ( snow_h .GT. zero ) THEN
1939	pcapa_en(ji,1) = sn_capa
1940	zx1 = snow_h / zz_coef(1)
1941	zx2 = ( zz_coef(1) - snow_h) / zz_coef(1)
1942	pcapa(ji,1) = zx1 * sn_capa + zx2 * so_capa_wet
1943	pkappa(ji,1) = un / ( zx1 / sn_cond + zx2 / so_cond_wet )
1944	ELSE
1945	pcapa(ji,1) = so_capa_dry + shum_ngrnd_perma(ji,1)*(so_capa_wet - so_capa_dry)
1946	pkappa(ji,1) = so_cond_dry + shum_ngrnd_perma(ji,1)*(so_cond_wet - so_cond_dry)
1947	pcapa_en(ji,1) = so_capa_dry + shum_ngrnd_perma(ji,1)*(so_capa_wet - so_capa_dry)
1948	ENDIF
1949
1950	! Mid layers
1951	DO jg = 2, ngrnd - 2
1952	IF ( snow_h .GT. zz_coef(jg) ) THEN
1953	pcapa(ji,jg) = sn_capa
1954	pkappa(ji,jg) = sn_cond
1955	pcapa_en(ji,jg) = sn_capa
1956	ELSE IF ( snow_h .GT. zz_coef(jg-1) ) THEN
1957	zx1 = (snow_h - zz_coef(jg-1)) / (zz_coef(jg) - zz_coef(jg-1))
1958	zx2 = ( zz_coef(jg) - snow_h) / (zz_coef(jg) - zz_coef(jg-1))
1959	pcapa(ji,jg) = zx1 * sn_capa + zx2 * so_capa_wet
1960	pkappa(ji,jg) = un / ( zx1 / sn_cond + zx2 / so_cond_wet )
1961	pcapa_en(ji,jg) = sn_capa
1962	ELSE
1963	pcapa(ji,jg) = so_capa_dry + shum_ngrnd_perma(ji,jg)*(so_capa_wet - so_capa_dry)
1964	pkappa(ji,jg) = so_cond_dry + shum_ngrnd_perma(ji,jg)*(so_cond_wet - so_cond_dry)
1965	pcapa_en(ji,jg) = so_capa_dry + shum_ngrnd_perma(ji,jg)*(so_capa_wet - so_capa_dry)
1966	ENDIF
1967	ENDDO
1968
1969	! Last two layers: These layers can not be filled with snow
1970	DO jg = ngrnd - 1, ngrnd
1971	pcapa(ji,jg) = so_capa_dry
1972	pkappa(ji,jg) = so_cond_dry
1973	pcapa_en(ji,jg) = so_capa_dry
1974	END DO
1975
1976	ENDDO ! DO ji = 1,kjpindex
1977
1978
1979	END SUBROUTINE thermosoilc_getdiff_old_thermix_with_snow
1980
1981
1982
1983	!! ================================================================================================================================
1984	!! SUBROUTINE : thermosoilc_getdiff_old_thermix_without_snow
1985	!!
1986	!>\BRIEF Computes soil and snow heat capacity and conductivity
1987	!!
1988	!! DESCRIPTION : Calculations of soil and snow thermal properties without effect of freezing. This subroutine is only called
1989	!! when explicitsnow is activated.
1990	!!
1991	!!
1992	!! RECENT CHANGE(S) : None
1993	!!
1994	!! MAIN OUTPUT VARIABLE(S):
1995	!!
1996	!! REFERENCE(S) :
1997	!!
1998	!! FLOWCHART : None
1999	!! \n
2000	!_ ================================================================================================================================
2001
2002	SUBROUTINE thermosoilc_getdiff_old_thermix_without_snow( kjpindex,snowrho, snowtemp, pb )
2003
2004	!! 0. Variables and parameter declaration
2005
2006	!! 0.1 Input variables
2007	INTEGER(i_std), INTENT(in) :: kjpindex
2008	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(in) :: snowrho !! Snow density
2009	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(in) :: snowtemp !! Snow temperature
2010	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: pb !! Surface presure (hPa)
2011
2012
2013	!! 0.1 Local variables
2014	INTEGER :: ji,jg
2015
2016	! soil
2017	DO jg = 1,ngrnd
2018	DO ji = 1,kjpindex
2019	pkappa(ji,jg) = so_cond_dry + shum_ngrnd_perma(ji,jg)*(so_cond_wet - so_cond_dry)
2020	pcapa(ji,jg) = so_capa_dry + shum_ngrnd_perma(ji,jg)*(so_capa_wet - so_capa_dry)
2021	pcapa_en(ji,jg) = so_capa_dry + shum_ngrnd_perma(ji,jg)*(so_capa_wet - so_capa_dry)
2022	ENDDO
2023	ENDDO
2024
2025	! snow
2026	DO ji = 1,kjpindex
2027	pcapa_snow(ji,:) = snowrho(ji,:) * xci
2028	pkappa_snow(ji,:) = (ZSNOWTHRMCOND1 + ZSNOWTHRMCOND2snowrho(ji,:)snowrho(ji,:)) + &
2029	MAX(0.0,(ZSNOWTHRMCOND_AVAP+(ZSNOWTHRMCOND_BVAP/(snowtemp(ji,:)+ &
2030	ZSNOWTHRMCOND_CVAP)))(XP00/(pb(ji)100.)))
2031	END DO
2032
2033
2034	END SUBROUTINE thermosoilc_getdiff_old_thermix_without_snow
2035
2036
2037	!! ================================================================================================================================
2038	!! SUBROUTINE : thermosoilc_read_reftempfile
2039	!!
2040	!>\BRIEF
2041	!!
2042	!! DESCRIPTION : Read file with longterm temperature
2043	!!
2044	!!
2045	!! RECENT CHANGE(S) : None
2046	!!
2047	!! MAIN OUTPUT VARIABLE(S): reftemp : Reference temerature
2048	!!
2049	!! REFERENCE(S) :
2050	!!
2051	!! FLOWCHART : None
2052	!! \n
2053	!_ ================================================================================================================================
2054
2055	SUBROUTINE thermosoilc_read_reftempfile(kjpindex,lalo,reftemp)
2056
2057	USE interpweight
2058
2059	IMPLICIT NONE
2060	!! 0. Variables and parameter declaration
2061
2062	!! 0.1 Input variables
2063	INTEGER(i_std), INTENT(in) :: kjpindex
2064	REAL(r_std), DIMENSION(kjpindex,2), INTENT(in) :: lalo
2065
2066	!! 0.2 Output variables
2067	REAL(r_std), DIMENSION(kjpindex, ngrnd), INTENT(out) :: reftemp
2068	REAL(r_std), DIMENSION(kjpindex) :: areftemp !! Availability of data for the interpolation
2069
2070	!! 0.3 Local variables
2071	INTEGER(i_std) :: ib
2072	CHARACTER(LEN=80) :: filename
2073	REAL(r_std),DIMENSION(kjpindex) :: reftemp_file !! Horizontal temperature field interpolated from file [C]
2074	REAL(r_std) :: vmin, vmax !! min/max values to use for the
2075	!! renormalization
2076	CHARACTER(LEN=80) :: variablename !! Variable to interpolate
2077	!! the file
2078	CHARACTER(LEN=80) :: lonname, latname !! lon, lat names in input file
2079	REAL(r_std), DIMENSION(:), ALLOCATABLE :: variabletypevals !! Values for all the types of the variable
2080	!! (variabletypevals(1) = -un, not used)
2081	CHARACTER(LEN=50) :: fractype !! method of calculation of fraction
2082	!! 'XYKindTime': Input values are kinds
2083	!! of something with a temporal
2084	!! evolution on the dx*dy matrix'
2085	LOGICAL :: nonegative !! whether negative values should be removed
2086	CHARACTER(LEN=50) :: maskingtype !! Type of masking
2087	!! 'nomask': no-mask is applied
2088	!! 'mbelow': take values below maskvals(1)
2089	!! 'mabove': take values above maskvals(1)
2090	!! 'msumrange': take values within 2 ranges;
2091	!! maskvals(2) <= SUM(vals(k)) <= maskvals(1)
2092	!! maskvals(1) < SUM(vals(k)) <= maskvals(3)
2093	!! (normalized by maskvals(3))
2094	!! 'var': mask values are taken from a
2095	!! variable inside the file (>0)
2096	REAL(r_std), DIMENSION(3) :: maskvals !! values to use to mask (according to
2097	!! `maskingtype')
2098	CHARACTER(LEN=250) :: namemaskvar !! name of the variable to use to mask
2099	REAL(r_std) :: reftemp_norefinf
2100	REAL(r_std) :: reftemp_default !! Default value
2101
2102
2103	!Config Key = SOIL_REFTEMP_FILE
2104	!Config Desc = File with climatological soil temperature
2105	!Config If = READ_REFTEMP
2106	!Config Def = reftemp.nc
2107	!Config Help =
2108	!Config Units = [FILE]
2109	filename = 'reftemp.nc'
2110	CALL getin_p('REFTEMP_FILE',filename)
2111
2112	variablename = 'temperature'
2113
2114	IF (printlev >= 1) WRITE(numout,*) "thermosoilc_read_reftempfile: Read and interpolate file " &
2115	// TRIM(filename) //" for variable " //TRIM(variablename)
2116
2117	! For this case there are not types/categories. We have 'only' a continuos
2118	! field
2119	! Assigning values to vmin, vmax
2120
2121	vmin = 0.
2122	vmax = 9999.
2123
2124	! For this file we do not need neightbours!
2125	neighbours = 0
2126
2127	!! Variables for interpweight
2128	! Type of calculation of cell fractions
2129	fractype = 'default'
2130	! Name of the longitude and latitude in the input file
2131	lonname = 'nav_lon'
2132	latname = 'nav_lat'
2133	! Default value when no value is get from input file
2134	reftemp_default = 1.
2135	! Reference value when no value is get from input file
2136	reftemp_norefinf = 1.
2137	! Should negative values be set to zero from input file?
2138	nonegative = .FALSE.
2139	! Type of mask to apply to the input data (see header for more details)
2140	maskingtype = 'nomask'
2141	! Values to use for the masking (here not used)
2142	maskvals = (/ undef_sechiba, undef_sechiba, undef_sechiba /)
2143	! Name of the variable with the values for the mask in the input file (only if maskkingtype='var') (here not used)
2144	namemaskvar = ''
2145
2146	CALL interpweight_2Dcont(kjpindex, 0, 0, lalo, resolution, neighbours, &
2147	contfrac, filename, variablename, lonname, latname, vmin, vmax, nonegative, maskingtype, &
2148	maskvals, namemaskvar, -1, fractype, reftemp_default, reftemp_norefinf, &
2149	reftemp_file, areftemp)
2150	IF (printlev >= 5) WRITE(numout,*)' thermosoilc_read_reftempfile after interpweight2D_cont'
2151
2152	! Copy reftemp_file temperature to all ground levels and transform into Kelvin
2153	DO ib=1, kjpindex
2154	reftemp(ib, :) = reftemp_file(ib)+ZeroCelsius
2155	END DO
2156
2157	! Write diagnostics
2158	CALL xios_orchidee_send_field("areftemp",areftemp)
2159
2160	END SUBROUTINE thermosoilc_read_reftempfile
2161
2162
2163	END MODULE thermosoilc

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

Download in other formats: