Context Navigation

source: branches/publications/ORCHIDEE_SOM_13C_r4859/src_sechiba/thermosoilc.f90 @ 7390

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

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

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

Download in other formats: