Context Navigation

source: branches/publications/ORCHIDEE_DFv1.0_site/src_sechiba/thermosoil.f90 @ 7346

Last change on this file since 7346 was 5433, checked in by josefine.ghattas, 6 years ago
Initialize pcappa_supp to zero at first time step to avoid crash with debug options as intialization to nan. This variable is output with histwrite(no xios send_field for the moment) at first time step before it is calculated. This output variable is therfor not correct restarted. The variable is always calculated before it is used in readadjust.
Property svn:keywords set to `Revision Date HeadURL Date Author Revision`
File size: 107.8 KB

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

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

Download in other formats: