Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

thermosoilc.f90 @ 5143

Last change on this file since 5143 was 3965, checked in by jan.polcher, 8 years ago

Merge with trunk at revision3959.
This includes all the developments made for CMIP6 and passage to XIOS2.
All conflicts are resolved and the code compiles.

But it still does not link because of an "undefined reference to `_intel_fast_memmove'"

Property svn:keywords set to Revision Date HeadURL Date Author Revision

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

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

Context Navigation

source: branches/publications/ORCHIDEE_gmd-2018-57/src_sechiba/thermosoilc.f90 @ 5143

Download in other formats: