Home

Context Navigation

land_surface_model_mod.f90 @ 2608

Last change on this file since 2608 was 2608, checked in by schwenkel, 5 years ago
Inital revision of diagnostic_quantities_mod allows unified calculation of magnus equation and saturion mixing ratio
Property svn:keywords set to `Id`
File size: 195.0 KB

Line
1	!> @file land_surface_model_mod.f90
2	!------------------------------------------------------------------------------!
3	! This file is part of PALM.
4	!
5	! PALM is free software: you can redistribute it and/or modify it under the
6	! terms of the GNU General Public License as published by the Free Software
7	! Foundation, either version 3 of the License, or (at your option) any later
8	! version.
9	!
10	! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11	! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12	! A PARTICULAR PURPOSE. See the GNU General Public License for more details.
13	!
14	! You should have received a copy of the GNU General Public License along with
15	! PALM. If not, see <http://www.gnu.org/licenses/>.
16	!
17	! Copyright 1997-2017 Leibniz Universitaet Hannover
18	!------------------------------------------------------------------------------!
19	!
20	! Current revisions:
21	! -----------------
22	!
23	!
24	! Former revisions:
25	! -----------------
26	! $Id: land_surface_model_mod.f90 2608 2017-11-13 14:04:26Z schwenkel $
27	! Calculation of magnus equation in external module (diagnostic_quantities_mod).
28	! Adjust calculation of vapor pressure and saturation mixing ratio that it is
29	! consistent with formulations in other parts of PALM.
30	!
31	! 2575 2017-10-24 09:57:58Z maronga
32	! Pavement parameterization revised
33	!
34	! 2573 2017-10-20 15:57:49Z scharf
35	! bugfixes in last_actions
36	!
37	! 2548 2017-10-16 13:18:20Z suehring
38	! extended by cloud_droplets option
39	!
40	! 2532 2017-10-11 16:00:46Z scharf
41	! bugfixes in data_output_3d
42	!
43	! 2516 2017-10-04 11:03:04Z suehring
44	! Remove tabs
45	!
46	! 2514 2017-10-04 09:52:37Z suehring
47	! upper bounds of cross section and 3d output changed from nx+1,ny+1 to nx,ny
48	! no output of ghost layer data
49	!
50	! 2504 2017-09-27 10:36:13Z maronga
51	! Support roots and water under pavement. Added several pavement types.
52	!
53	! 2476 2017-09-18 07:54:32Z maronga
54	! Bugfix for last commit
55	!
56	! 2475 2017-09-18 07:42:36Z maronga
57	! Bugfix: setting of vegetation_pars for bare soil corrected.
58	!
59	! 2354 2017-08-17 10:49:36Z schwenkel
60	! minor bugfixes
61	!
62	! 2340 2017-08-07 17:11:13Z maronga
63	! Revised root_distribution tabel and implemented a pseudo-generic root fraction
64	! calculation
65	!
66	! 2333 2017-08-04 09:08:26Z maronga
67	! minor bugfixes
68	!
69	! 2332 2017-08-03 21:15:22Z maronga
70	! bugfix in pavement_pars
71	!
72	! 2328 2017-08-03 12:34:22Z maronga
73	! Revised skin layer concept.
74	! Bugfix for runs with pavement surface and humidity
75	! Revised some standard values in vegetation_pars
76	! Added emissivity and default albedo_type as variable to tables
77	! Changed default surface type to vegetation
78	! Revised input of soil layer configuration
79	!
80	! 2307 2017-07-07 11:32:10Z suehring
81	! Bugfix, variable names corrected
82	!
83	! 2299 2017-06-29 10:14:38Z maronga
84	! Removed pt_p from USE statement. Adjusted call to lsm_soil_model to allow
85	! spinups without soil moisture prediction
86	!
87	! 2298 2017-06-29 09:28:18Z raasch
88	! type of write_binary changed from CHARACTER to LOGICAL
89	!
90	! 2296 2017-06-28 07:53:56Z maronga
91	! Bugfix in calculation of bare soil heat capacity.
92	! Bugfix in calculation of shf
93	! Added support for spinups
94	!
95	! 2282 2017-06-13 11:38:46Z schwenkel
96	! Bugfix for check of saturation moisture
97	!
98	! 2273 2017-06-09 12:46:06Z sward
99	! Error number changed
100	!
101	! 2270 2017-06-09 12:18:47Z maronga
102	! Revised parameterization of heat conductivity between skin layer and soil.
103	! Temperature and moisture are now defined at the center of the layers.
104	! Renamed veg_type to vegetation_type and pave_type to pavement_type_name
105	! Renamed and reduced the number of look-up tables (vegetation_pars, soil_pars)
106	! Revised land surface model initialization
107	! Removed output of shf_eb and qsws_eb and removed _eb throughout code
108	! Removed Clapp & Hornberger parameterization
109	!
110	! 2249 2017-06-06 13:58:01Z sward
111	!
112	! 2248 2017-06-06 13:52:54Z sward $
113	! Error no changed
114	!
115	! 2246 2017-06-06 13:09:34Z sward
116	! Error no changed
117	!
118	! Changed soil configuration to 8 layers. The number of soil layers is now
119	! freely adjustable via the NAMELIST.
120	!
121	! 2237 2017-05-31 10:34:53Z suehring
122	! Bugfix in write restart data
123	!
124	! 2233 2017-05-30 18:08:54Z suehring
125	!
126	! 2232 2017-05-30 17:47:52Z suehring
127	! Adjustments to new topography and surface concept
128	! - now, also vertical walls are possible
129	! - for vertical walls, parametrization of r_a (aerodynamic resisistance) is
130	! implemented.
131	!
132	! Add check for soil moisture, it must not exceed its saturation value.
133	!
134	! 2149 2017-02-09 16:57:03Z scharf
135	! Land surface parameters II corrected for vegetation_type 18 and 19
136	!
137	! 2031 2016-10-21 15:11:58Z knoop
138	! renamed variable rho to rho_ocean
139	!
140	! 2000 2016-08-20 18:09:15Z knoop
141	! Forced header and separation lines into 80 columns
142	!
143	! 1978 2016-07-29 12:08:31Z maronga
144	! Bugfix: initial values of pave_surface and water_surface were not set.
145	!
146	! 1976 2016-07-27 13:28:04Z maronga
147	! Parts of the code have been reformatted. Use of radiation model output is
148	! generalized and simplified. Added more output quantities due to modularization
149	!
150	! 1972 2016-07-26 07:52:02Z maronga
151	! Further modularization: output of cross sections and 3D data is now done in this
152	! module. Moreover, restart data is written and read directly within this module.
153	!
154	!
155	! 1966 2016-07-18 11:54:18Z maronga
156	! Bugfix: calculation of m_total in soil model was not set to zero at model start
157	!
158	! 1949 2016-06-17 07:19:16Z maronga
159	! Bugfix: calculation of qsws_soil_eb with precipitation = .TRUE. gave
160	! qsws_soil_eb = 0 due to a typo
161	!
162	! 1856 2016-04-13 12:56:17Z maronga
163	! Bugfix: for water surfaces, the initial water surface temperature is set equal
164	! to the intital skin temperature. Moreover, the minimum value of r_a is now
165	! 1.0 to avoid too large fluxes at the first model time step
166	!
167	! 1849 2016-04-08 11:33:18Z hoffmann
168	! prr moved to arrays_3d
169	!
170	! 1826 2016-04-07 12:01:39Z maronga
171	! Cleanup after modularization
172	!
173	! 1817 2016-04-06 15:44:20Z maronga
174	! Added interface for lsm_init_arrays. Added subroutines for check_parameters,
175	! header, and parin. Renamed some subroutines.
176	!
177	! 1788 2016-03-10 11:01:04Z maronga
178	! Bugfix: calculate lambda_surface based on temperature gradient between skin
179	! layer and soil layer instead of Obukhov length
180	! Changed: moved calculation of surface specific humidity to energy balance solver
181	! New: water surfaces are available by using a fixed sea surface temperature.
182	! The roughness lengths are calculated dynamically using the Charnock
183	! parameterization. This involves the new roughness length for moisture z0q.
184	! New: modified solution of the energy balance solver and soil model for
185	! paved surfaces (i.e. asphalt concrete).
186	! Syntax layout improved.
187	! Changed: parameter dewfall removed.
188	!
189	! 1783 2016-03-06 18:36:17Z raasch
190	! netcdf variables moved to netcdf module
191	!
192	! 1757 2016-02-22 15:49:32Z maronga
193	! Bugfix: set tm_soil_m to zero after allocation. Added parameter
194	! unscheduled_radiation_calls to control calls of the radiation model based on
195	! the skin temperature change during one time step (preliminary version). Set
196	! qsws_soil_eb to zero at model start (previously set to qsws_eb). Removed MAX
197	! function as it cannot be vectorized.
198	!
199	! 1709 2015-11-04 14:47:01Z maronga
200	! Renamed pt_1 and qv_1 to pt1 and qv1.
201	! Bugfix: set initial values for t_surface_p in case of restart runs
202	! Bugfix: zero resistance caused crash when using radiation_scheme = 'clear-sky'
203	! Bugfix: calculation of rad_net when using radiation_scheme = 'clear-sky'
204	! Added todo action
205	!
206	! 1697 2015-10-28 17:14:10Z raasch
207	! bugfix: misplaced cpp-directive
208	!
209	! 1695 2015-10-27 10:03:11Z maronga
210	! Bugfix: REAL constants provided with KIND-attribute in call of
211	! Replaced rif with ol
212	!
213	! 1691 2015-10-26 16:17:44Z maronga
214	! Added skip_time_do_lsm to allow for spin-ups without LSM. Various bugfixes:
215	! Soil temperatures are now defined at the edges of the layers, calculation of
216	! shb_eb corrected, prognostic equation for skin temperature corrected. Surface
217	! fluxes are now directly transfered to atmosphere
218	!
219	! 1682 2015-10-07 23:56:08Z knoop
220	! Code annotations made doxygen readable
221	!
222	! 1590 2015-05-08 13:56:27Z maronga
223	! Bugfix: definition of character strings requires same length for all elements
224	!
225	! 1585 2015-04-30 07:05:52Z maronga
226	! Modifications for RRTMG. Changed tables to PARAMETER type.
227	!
228	! 1571 2015-03-12 16:12:49Z maronga
229	! Removed upper-case variable names. Corrected distribution of precipitation to
230	! the liquid water reservoir and the bare soil fractions.
231	!
232	! 1555 2015-03-04 17:44:27Z maronga
233	! Added output of r_a and r_s
234	!
235	! 1553 2015-03-03 17:33:54Z maronga
236	! Improved better treatment of roughness lengths. Added default soil temperature
237	! profile
238	!
239	! 1551 2015-03-03 14:18:16Z maronga
240	! Flux calculation is now done in prandtl_fluxes. Added support for data output.
241	! Vertical indices have been replaced. Restart runs are now possible. Some
242	! variables have beem renamed. Bugfix in the prognostic equation for the surface
243	! temperature. Introduced z0_eb and z0h_eb, which overwrite the setting of
244	! roughness_length and z0_factor. Added Clapp & Hornberger parametrization for
245	! the hydraulic conductivity. Bugfix for root fraction and extraction
246	! calculation
247	!
248	! intrinsic function MAX and MIN
249	!
250	! 1500 2014-12-03 17:42:41Z maronga
251	! Corrected calculation of aerodynamic resistance (r_a).
252	! Precipitation is now added to liquid water reservoir using LE_liq.
253	! Added support for dry runs.
254	!
255	! 1496 2014-12-02 17:25:50Z maronga
256	! Initial revision
257	!
258	!
259	! Description:
260	! ------------
261	!> Land surface model, consisting of a solver for the energy balance at the
262	!> surface and a multi layer soil scheme. The scheme is similar to the TESSEL
263	!> scheme implemented in the ECMWF IFS model, with modifications according to
264	!> H-TESSEL. The implementation is based on the formulation implemented in the
265	!> DALES and UCLA-LES models.
266	!>
267	!> @todo Restart data for vertical natural land-surfaces
268	!> @todo Extensive verification energy-balance solver for vertical surfaces,
269	!> e.g. parametrization of r_a
270	!> @todo Revise single land-surface processes for vertical surfaces, e.g.
271	!> treatment of humidity, etc.
272	!> @todo Consider partial absorption of the net shortwave radiation by the
273	!> skin layer.
274	!> @todo Improve surface water parameterization
275	!> @todo Invert indices (running from -3 to 0. Currently: nzb_soil=0,
276	!> nzt_soil=3)).
277	!> @todo Implement surface runoff model (required when performing long-term LES
278	!> with considerable precipitation.
279	!> @todo Revise calculation of f2 when wilting point is non-constant in the
280	!> soil
281	!> @todo Allow for zero soil moisture (currently, it is set to wilting point)
282	!> @note No time step criterion is required as long as the soil layers do not
283	!> become too thin.
284	!------------------------------------------------------------------------------!
285	MODULE land_surface_model_mod
286
287	USE arrays_3d, &
288	ONLY: hyp, pt, prr, q, q_p, ql, vpt, u, v, w
289
290	USE cloud_parameters, &
291	ONLY: cp, hyrho, l_d_cp, l_d_r, l_v, pt_d_t, rho_l, r_d, r_v
292
293	USE control_parameters, &
294	ONLY: cloud_droplets, cloud_physics, coupling_start_time, dt_3d, &
295	end_time, humidity, intermediate_timestep_count, &
296	initializing_actions, intermediate_timestep_count_max, &
297	land_surface, max_masks, precipitation, pt_surface, &
298	rho_surface, spinup, spinup_pt_mean, spinup_time, &
299	surface_pressure, timestep_scheme, tsc, &
300	time_since_reference_point
301
302	USE indices, &
303	ONLY: nbgp, nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb
304
305	USE kinds
306
307	USE pegrid
308
309	USE radiation_model_mod, &
310	ONLY: albedo, albedo_type, emissivity, force_radiation_call, rad_net, &
311	rad_sw_in, rad_lw_out, rad_lw_out_change_0, radiation_scheme, &
312	unscheduled_radiation_calls
313
314	USE statistics, &
315	ONLY: hom, statistic_regions
316
317	USE surface_mod, &
318	ONLY : surf_lsm_h, surf_lsm_v, surf_type
319
320	IMPLICIT NONE
321
322	TYPE surf_type_lsm
323	REAL(wp), DIMENSION(:), ALLOCATABLE :: var_1D !< 1D prognostic variable
324	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: var_2D !< 2D prognostic variable
325	END TYPE surf_type_lsm
326
327	!
328	!-- LSM model constants
329
330	REAL(wp), PARAMETER :: &
331	b_ch = 6.04_wp, & ! Clapp & Hornberger exponent
332	lambda_h_dry = 0.19_wp, & ! heat conductivity for dry soil (W/m/K)
333	lambda_h_sm = 3.44_wp, & ! heat conductivity of the soil matrix (W/m/K)
334	lambda_h_water = 0.57_wp, & ! heat conductivity of water (W/m/K)
335	psi_sat = -0.388_wp, & ! soil matrix potential at saturation
336	rho_c_soil = 2.19E6_wp, & ! volumetric heat capacity of soil (J/m3/K)
337	rho_c_water = 4.20E6_wp, & ! volumetric heat capacity of water (J/m3/K)
338	m_max_depth = 0.0002_wp ! Maximum capacity of the water reservoir (m)
339
340
341	REAL(wp), DIMENSION(0:7), PARAMETER :: dz_soil_default = & ! default soil layer configuration
342	(/ 0.01_wp, 0.02_wp, 0.04_wp, &
343	0.06_wp, 0.14_wp, 0.26_wp, &
344	0.54_wp, 1.86_wp/)
345
346	REAL(wp), DIMENSION(0:3), PARAMETER :: dz_soil_ref = & ! reference four layer soil configuration used for estimating the root fractions
347	(/ 0.07_wp, 0.21_wp, 0.72_wp, &
348	1.89_wp /)
349
350	REAL(wp), DIMENSION(0:3), PARAMETER :: zs_ref = & ! reference four layer soil configuration used for estimating the root fractions
351	(/ 0.07_wp, 0.28_wp, 1.0_wp, &
352	2.89_wp /)
353
354
355	!
356	!-- LSM variables
357	CHARACTER(10) :: surface_type = 'vegetation' !< general classification. Allowed are:
358	!< 'vegetation', 'pavement', ('building'),
359	!< 'water', and 'netcdf'
360
361
362
363	INTEGER(iwp) :: nzb_soil = 0, & !< bottom of the soil model (Earth's surface)
364	nzt_soil = 7, & !< top of the soil model
365	nzt_pavement = 0, & !< top of the pavement within the soil
366	nzs = 8, & !< number of soil layers
367	pavement_depth_level = 0, & !< default NAMELIST nzt_pavement
368	pavement_type = 1, & !< default NAMELIST pavement_type
369	soil_type = 3, & !< default NAMELIST soil_type
370	vegetation_type = 2, & !< default NAMELIST vegetation_type
371	water_type = 1 !< default NAMELISt water_type
372
373
374
375	LOGICAL :: conserve_water_content = .TRUE., & !< open or closed bottom surface for the soil model
376	constant_roughness = .FALSE., & !< use fixed/dynamic roughness lengths for water surfaces
377	force_radiation_call_l = .FALSE., & !< flag to force calling of radiation routine
378	aero_resist_kray = .TRUE. !< flag to control parametrization of aerodynamic resistance at vertical surface elements
379
380	! value 9999999.9_wp -> generic available or user-defined value must be set
381	! otherwise -> no generic variable and user setting is optional
382	REAL(wp) :: alpha_vangenuchten = 9999999.9_wp, & !< NAMELIST alpha_vg
383	canopy_resistance_coefficient = 9999999.9_wp, & !< NAMELIST g_d
384	c_surface = 9999999.9_wp, & !< Surface (skin) heat capacity (J/m2/K)
385	drho_l_lv, & !< (rho_l * l_v)**-1
386	exn, & !< value of the Exner function
387	e_s = 0.0_wp, & !< saturation water vapour pressure
388	field_capacity = 9999999.9_wp, & !< NAMELIST m_fc
389	f_shortwave_incoming = 9999999.9_wp, & !< NAMELIST f_sw_in
390	hydraulic_conductivity = 9999999.9_wp, & !< NAMELIST gamma_w_sat
391	ke = 0.0_wp, & !< Kersten number
392	lambda_h_sat = 0.0_wp, & !< heat conductivity for saturated soil (W/m/K)
393	lambda_surface_stable = 9999999.9_wp, & !< NAMELIST lambda_surface_s (W/m2/K)
394	lambda_surface_unstable = 9999999.9_wp, & !< NAMELIST lambda_surface_u (W/m2/K)
395	leaf_area_index = 9999999.9_wp, & !< NAMELIST lai
396	l_vangenuchten = 9999999.9_wp, & !< NAMELIST l_vg
397	min_canopy_resistance = 9999999.9_wp, & !< NAMELIST r_canopy_min
398	min_soil_resistance = 50.0_wp, & !< NAMELIST r_soil_min
399	m_total = 0.0_wp, & !< weighted total water content of the soil (m3/m3)
400	n_vangenuchten = 9999999.9_wp, & !< NAMELIST n_vg
401	pavement_heat_capacity = 9999999.9_wp, & !< volumetric heat capacity of pavement (e.g. roads) (J/m3/K)
402	pavement_heat_conduct = 9999999.9_wp, & !< heat conductivity for pavements (e.g. roads) (W/m/K)
403	q_s = 0.0_wp, & !< saturation specific humidity
404	residual_moisture = 9999999.9_wp, & !< NAMELIST m_res
405	rho_cp, & !< rho_surface * cp
406	rho_lv, & !< rho_ocean * l_v
407	rd_d_rv, & !< r_d / r_v
408	saturation_moisture = 9999999.9_wp, & !< NAMELIST m_sat
409	skip_time_do_lsm = 0.0_wp, & !< LSM is not called before this time
410	vegetation_coverage = 9999999.9_wp, & !< NAMELIST c_veg
411	water_temperature = 9999999.9_wp, & !< water temperature
412	wilting_point = 9999999.9_wp, & !< NAMELIST m_wilt
413	z0_vegetation = 9999999.9_wp, & !< NAMELIST z0 (lsm_par)
414	z0h_vegetation = 9999999.9_wp, & !< NAMELIST z0h (lsm_par)
415	z0q_vegetation = 9999999.9_wp, & !< NAMELIST z0q (lsm_par)
416	z0_pavement = 9999999.9_wp, & !< NAMELIST z0 (lsm_par)
417	z0h_pavement = 9999999.9_wp, & !< NAMELIST z0h (lsm_par)
418	z0q_pavement = 9999999.9_wp, & !< NAMELIST z0q (lsm_par)
419	z0_water = 9999999.9_wp, & !< NAMELIST z0 (lsm_par)
420	z0h_water = 9999999.9_wp, & !< NAMELIST z0h (lsm_par)
421	z0q_water = 9999999.9_wp !< NAMELIST z0q (lsm_par)
422
423
424	REAL(wp), DIMENSION(:), ALLOCATABLE :: ddz_soil_center, & !< 1/dz_soil_center
425	ddz_soil, & !< 1/dz_soil
426	dz_soil_center, & !< soil grid spacing (center-center)
427	zs, & !< depth of the temperature/moisute levels
428	root_extr !< root extraction
429
430
431
432	REAL(wp), DIMENSION(0:20) :: root_fraction = 9999999.9_wp, & !< (NAMELIST) distribution of root surface area to the individual soil layers
433	soil_moisture = 0.0_wp, & !< NAMELIST soil moisture content (m3/m3)
434	soil_temperature = 300.0_wp, & !< NAMELIST soil temperature (K) +1
435	dz_soil = 9999999.9_wp, & !< (NAMELIST) soil layer depths (spacing)
436	zs_layer = 9999999.9_wp !< soil layer depths (edge)
437
438	#if defined( __nopointer )
439	TYPE(surf_type_lsm), TARGET :: t_soil_h, & !< Soil temperature (K), horizontal surface elements
440	t_soil_h_p, & !< Prog. soil temperature (K), horizontal surface elements
441	m_soil_h, & !< Soil moisture (m3/m3), horizontal surface elements
442	m_soil_h_p !< Prog. soil moisture (m3/m3), horizontal surface elements
443
444	TYPE(surf_type_lsm), DIMENSION(0:3), TARGET :: &
445	t_soil_v, & !< Soil temperature (K), vertical surface elements
446	t_soil_v_p, & !< Prog. soil temperature (K), vertical surface elements
447	m_soil_v, & !< Soil moisture (m3/m3), vertical surface elements
448	m_soil_v_p !< Prog. soil moisture (m3/m3), vertical surface elements
449	#else
450	TYPE(surf_type_lsm), POINTER :: t_soil_h, & !< Soil temperature (K), horizontal surface elements
451	t_soil_h_p, & !< Prog. soil temperature (K), horizontal surface elements
452	m_soil_h, & !< Soil moisture (m3/m3), horizontal surface elements
453	m_soil_h_p !< Prog. soil moisture (m3/m3), horizontal surface elements
454
455	TYPE(surf_type_lsm), TARGET :: t_soil_h_1, & !<
456	t_soil_h_2, & !<
457	m_soil_h_1, & !<
458	m_soil_h_2 !<
459
460	TYPE(surf_type_lsm), DIMENSION(:), POINTER :: &
461	t_soil_v, & !< Soil temperature (K), vertical surface elements
462	t_soil_v_p, & !< Prog. soil temperature (K), vertical surface elements
463	m_soil_v, & !< Soil moisture (m3/m3), vertical surface elements
464	m_soil_v_p !< Prog. soil moisture (m3/m3), vertical surface elements
465
466	TYPE(surf_type_lsm), DIMENSION(0:3), TARGET ::&
467	t_soil_v_1, & !<
468	t_soil_v_2, & !<
469	m_soil_v_1, & !<
470	m_soil_v_2 !<
471	#endif
472
473	#if defined( __nopointer )
474	TYPE(surf_type_lsm), TARGET :: t_surface_h, & !< surface temperature (K), horizontal surface elements
475	t_surface_h_p, & !< progn. surface temperature (K), horizontal surface elements
476	m_liq_h, & !< liquid water reservoir (m), horizontal surface elements
477	m_liq_h_p !< progn. liquid water reservoir (m), horizontal surface elements
478
479	TYPE(surf_type_lsm), DIMENSION(0:3), TARGET :: &
480	t_surface_v, & !< surface temperature (K), vertical surface elements
481	t_surface_v_p, & !< progn. surface temperature (K), vertical surface elements
482	m_liq_v, & !< liquid water reservoir (m), vertical surface elements
483	m_liq_v_p !< progn. liquid water reservoir (m), vertical surface elements
484	#else
485	TYPE(surf_type_lsm), POINTER :: t_surface_h, & !< surface temperature (K), horizontal surface elements
486	t_surface_h_p, & !< progn. surface temperature (K), horizontal surface elements
487	m_liq_h, & !< liquid water reservoir (m), horizontal surface elements
488	m_liq_h_p !< progn. liquid water reservoir (m), horizontal surface elements
489
490	TYPE(surf_type_lsm), TARGET :: t_surface_h_1, & !<
491	t_surface_h_2, & !<
492	m_liq_h_1, & !<
493	m_liq_h_2 !<
494
495	TYPE(surf_type_lsm), DIMENSION(:), POINTER :: &
496	t_surface_v, & !< surface temperature (K), vertical surface elements
497	t_surface_v_p, & !< progn. surface temperature (K), vertical surface elements
498	m_liq_v, & !< liquid water reservoir (m), vertical surface elements
499	m_liq_v_p !< progn. liquid water reservoir (m), vertical surface elements
500
501	TYPE(surf_type_lsm), DIMENSION(0:3), TARGET :: &
502	t_surface_v_1, & !<
503	t_surface_v_2, & !<
504	m_liq_v_1, & !<
505	m_liq_v_2 !<
506	#endif
507
508	#if defined( __nopointer )
509	REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET :: m_liq_av
510	#else
511	REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET :: m_liq_av
512	#endif
513
514	#if defined( __nopointer )
515	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: t_soil_av, & !< Average of t_soil
516	m_soil_av !< Average of m_soil
517	#else
518	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: t_soil_av, & !< Average of t_soil
519	m_soil_av !< Average of m_soil
520	#endif
521
522	TYPE(surf_type_lsm), TARGET :: tm_liq_h_m !< liquid water reservoir tendency (m), horizontal surface elements
523	TYPE(surf_type_lsm), TARGET :: tt_surface_h_m !< surface temperature tendency (K), horizontal surface elements
524	TYPE(surf_type_lsm), TARGET :: tt_soil_h_m !< t_soil storage array, horizontal surface elements
525	TYPE(surf_type_lsm), TARGET :: tm_soil_h_m !< m_soil storage array, horizontal surface elements
526
527	TYPE(surf_type_lsm), DIMENSION(0:3), TARGET :: tm_liq_v_m !< liquid water reservoir tendency (m), vertical surface elements
528	TYPE(surf_type_lsm), DIMENSION(0:3), TARGET :: tt_surface_v_m !< surface temperature tendency (K), vertical surface elements
529	TYPE(surf_type_lsm), DIMENSION(0:3), TARGET :: tt_soil_v_m !< t_soil storage array, vertical surface elements
530	TYPE(surf_type_lsm), DIMENSION(0:3), TARGET :: tm_soil_v_m !< m_soil storage array, vertical surface elements
531
532	!
533	!-- Energy balance variables
534	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: &
535	c_liq_av, & !< average of c_liq
536	c_soil_av, & !< average of c_soil
537	c_veg_av, & !< average of c_veg
538	ghf_av, & !< average of ghf
539	lai_av, & !< average of lai
540	qsws_liq_av, & !< average of qsws_liq
541	qsws_soil_av, & !< average of qsws_soil
542	qsws_veg_av, & !< average of qsws_veg
543	r_a_av, & !< average of r_a
544	r_s_av !< average of r_s
545
546
547	!
548	!-- Predefined Land surface classes (vegetation_type)
549	CHARACTER(26), DIMENSION(0:18), PARAMETER :: vegetation_type_name = (/ &
550	'user defined ', & ! 0
551	'bare soil ', & ! 1
552	'crops, mixed farming ', & ! 2
553	'short grass ', & ! 3
554	'evergreen needleleaf trees', & ! 4
555	'deciduous needleleaf trees', & ! 5
556	'evergreen broadleaf trees ', & ! 6
557	'deciduous broadleaf trees ', & ! 7
558	'tall grass ', & ! 8
559	'desert ', & ! 9
560	'tundra ', & ! 10
561	'irrigated crops ', & ! 11
562	'semidesert ', & ! 12
563	'ice caps and glaciers ', & ! 13
564	'bogs and marshes ', & ! 14
565	'evergreen shrubs ', & ! 15
566	'deciduous shrubs ', & ! 16
567	'mixed forest/woodland ', & ! 17
568	'interrupted forest ' & ! 18
569	/)
570
571	!
572	!-- Soil model classes (soil_type)
573	CHARACTER(12), DIMENSION(0:6), PARAMETER :: soil_type_name = (/ &
574	'user defined', & ! 0
575	'coarse ', & ! 1
576	'medium ', & ! 2
577	'medium-fine ', & ! 3
578	'fine ', & ! 4
579	'very fine ', & ! 5
580	'organic ' & ! 6
581	/)
582
583	!
584	!-- Pavement classes
585	CHARACTER(29), DIMENSION(0:15), PARAMETER :: pavement_type_name = (/ &
586	'user defined ', & ! 0
587	'asphalt/concrete mix ', & ! 1
588	'asphalt (asphalt concrete) ', & ! 2
589	'concrete (Portland concrete) ', & ! 3
590	'sett ', & ! 4
591	'paving stones ', & ! 5
592	'cobblestone ', & ! 6
593	'metal ', & ! 7
594	'wood ', & ! 8
595	'gravel ', & ! 9
596	'fine gravel ', & ! 10
597	'pebblestone ', & ! 11
598	'woodchips ', & ! 12
599	'tartan (sports) ', & ! 13
600	'artifical turf (sports) ', & ! 14
601	'clay (sports) ' & ! 15
602	/)
603
604
605
606
607
608	!
609	!-- Water classes
610	CHARACTER(12), DIMENSION(0:5), PARAMETER :: water_type_name = (/ &
611	'user defined', & ! 0
612	'lake ', & ! 1
613	'river ', & ! 2
614	'ocean ', & ! 3
615	'pond ', & ! 4
616	'fountain ' & ! 5
617	/)
618
619	!
620	!-- Land surface parameters according to the respective classes (vegetation_type)
621
622	!
623	!-- Land surface parameters
624	!-- r_canopy_min, lai, c_veg, g_d z0, z0h, lambda_s_s, lambda_s_u, f_sw_in, c_surface, albedo_type, emissivity
625	REAL(wp), DIMENSION(0:11,1:18), PARAMETER :: vegetation_pars = RESHAPE( (/ &
626	0.0_wp, 0.00_wp, 0.00_wp, 0.00_wp, 0.005_wp, 0.5E-4_wp, 0.0_wp, 0.0_wp, 0.00_wp, 0.00_wp, 17.0_wp, 0.94_wp, & ! 1
627	180.0_wp, 3.00_wp, 1.00_wp, 0.00_wp, 0.10_wp, 0.001_wp, 10.0_wp, 10.0_wp, 0.05_wp, 0.00_wp, 2.0_wp, 0.95_wp, & ! 2
628	110.0_wp, 2.00_wp, 1.00_wp, 0.00_wp, 0.03_wp, 0.3E-4_wp, 10.0_wp, 10.0_wp, 0.05_wp, 0.00_wp, 2.0_wp, 0.95_wp, & ! 3
629	500.0_wp, 5.00_wp, 1.00_wp, 0.03_wp, 2.00_wp, 2.00_wp, 20.0_wp, 15.0_wp, 0.03_wp, 0.00_wp, 5.0_wp, 0.97_wp, & ! 4
630	500.0_wp, 5.00_wp, 1.00_wp, 0.03_wp, 2.00_wp, 2.00_wp, 20.0_wp, 15.0_wp, 0.03_wp, 0.00_wp, 6.0_wp, 0.97_wp, & ! 5
631	175.0_wp, 5.00_wp, 1.00_wp, 0.03_wp, 2.00_wp, 2.00_wp, 20.0_wp, 15.0_wp, 0.03_wp, 0.00_wp, 8.0_wp, 0.97_wp, & ! 6
632	240.0_wp, 6.00_wp, 0.99_wp, 0.13_wp, 2.00_wp, 2.00_wp, 20.0_wp, 15.0_wp, 0.03_wp, 0.00_wp, 9.0_wp, 0.97_wp, & ! 7
633	100.0_wp, 2.00_wp, 0.70_wp, 0.00_wp, 0.47_wp, 0.47E-2_wp, 10.0_wp, 10.0_wp, 0.05_wp, 0.00_wp, 8.0_wp, 0.97_wp, & ! 8
634	250.0_wp, 0.05_wp, 0.00_wp, 0.00_wp, 0.013_wp, 0.013E-2_wp, 15.0_wp, 15.0_wp, 0.00_wp, 0.00_wp, 3.0_wp, 0.94_wp, & ! 9
635	80.0_wp, 1.00_wp, 0.50_wp, 0.00_wp, 0.034_wp, 0.034E-2_wp, 10.0_wp, 10.0_wp, 0.05_wp, 0.00_wp, 11.0_wp, 0.97_wp, & ! 10
636	180.0_wp, 3.00_wp, 1.00_wp, 0.00_wp, 0.5_wp, 0.50E-2_wp, 10.0_wp, 10.0_wp, 0.05_wp, 0.00_wp, 13.0_wp, 0.97_wp, & ! 11
637	150.0_wp, 0.50_wp, 0.10_wp, 0.00_wp, 0.17_wp, 0.17E-2_wp, 10.0_wp, 10.0_wp, 0.05_wp, 0.00_wp, 2.0_wp, 0.97_wp, & ! 12
638	0.0_wp, 0.00_wp, 0.00_wp, 0.00_wp, 1.3E-3_wp, 1.3E-4_wp, 58.0_wp, 58.0_wp, 0.00_wp, 0.00_wp, 11.0_wp, 0.97_wp, & ! 13
639	240.0_wp, 4.00_wp, 0.60_wp, 0.00_wp, 0.83_wp, 0.83E-2_wp, 10.0_wp, 10.0_wp, 0.05_wp, 0.00_wp, 4.0_wp, 0.97_wp, & ! 14
640	225.0_wp, 3.00_wp, 0.50_wp, 0.00_wp, 0.10_wp, 0.10E-2_wp, 10.0_wp, 10.0_wp, 0.05_wp, 0.00_wp, 4.0_wp, 0.97_wp, & ! 15
641	225.0_wp, 1.50_wp, 0.50_wp, 0.00_wp, 0.25_wp, 0.25E-2_wp, 10.0_wp, 10.0_wp, 0.05_wp, 0.00_wp, 4.0_wp, 0.97_wp, & ! 16
642	250.0_wp, 5.00_wp, 1.00_wp, 0.03_wp, 2.00_wp, 2.00_wp, 20.0_wp, 15.0_wp, 0.03_wp, 0.00_wp, 7.0_wp, 0.97_wp, & ! 17
643	175.0_wp, 2.50_wp, 1.00_wp, 0.03_wp, 1.10_wp, 1.10_wp, 20.0_wp, 15.0_wp, 0.03_wp, 0.00_wp, 8.0_wp, 0.97_wp & ! 18
644	/), (/ 12, 18 /) )
645
646
647	!
648	!-- Root distribution for default soil layer configuration (sum = 1)
649	!-- level 1 - level 4 according to zs_ref
650	REAL(wp), DIMENSION(0:3,1:18), PARAMETER :: root_distribution = RESHAPE( (/ &
651	1.00_wp, 0.00_wp, 0.00_wp, 0.00_wp, & ! 1
652	0.24_wp, 0.41_wp, 0.31_wp, 0.04_wp, & ! 2
653	0.35_wp, 0.38_wp, 0.23_wp, 0.04_wp, & ! 3
654	0.26_wp, 0.39_wp, 0.29_wp, 0.06_wp, & ! 4
655	0.26_wp, 0.38_wp, 0.29_wp, 0.07_wp, & ! 5
656	0.24_wp, 0.38_wp, 0.31_wp, 0.07_wp, & ! 6
657	0.25_wp, 0.34_wp, 0.27_wp, 0.14_wp, & ! 7
658	0.27_wp, 0.27_wp, 0.27_wp, 0.09_wp, & ! 8
659	1.00_wp, 0.00_wp, 0.00_wp, 0.00_wp, & ! 9
660	0.47_wp, 0.45_wp, 0.08_wp, 0.00_wp, & ! 10
661	0.24_wp, 0.41_wp, 0.31_wp, 0.04_wp, & ! 11
662	0.17_wp, 0.31_wp, 0.33_wp, 0.19_wp, & ! 12
663	0.00_wp, 0.00_wp, 0.00_wp, 0.00_wp, & ! 13
664	0.25_wp, 0.34_wp, 0.27_wp, 0.11_wp, & ! 14
665	0.23_wp, 0.36_wp, 0.30_wp, 0.11_wp, & ! 15
666	0.23_wp, 0.36_wp, 0.30_wp, 0.11_wp, & ! 16
667	0.19_wp, 0.35_wp, 0.36_wp, 0.10_wp, & ! 17
668	0.19_wp, 0.35_wp, 0.36_wp, 0.10_wp & ! 18
669	/), (/ 4, 18 /) )
670
671	!
672	!-- Soil parameters according to the following porosity classes (soil_type)
673
674	!
675	!-- Soil parameters alpha_vg, l_vg, n_vg, gamma_w_sat, m_sat, m_fc, m_wilt, m_res
676	REAL(wp), DIMENSION(0:7,1:6), PARAMETER :: soil_pars = RESHAPE( (/ &
677	3.83_wp, 1.250_wp, 1.38_wp, 6.94E-6_wp, 0.403_wp, 0.244_wp, 0.059_wp, 0.025_wp,& ! 1
678	3.14_wp, -2.342_wp, 1.28_wp, 1.16E-6_wp, 0.439_wp, 0.347_wp, 0.151_wp, 0.010_wp,& ! 2
679	0.83_wp, -0.588_wp, 1.25_wp, 0.26E-6_wp, 0.430_wp, 0.383_wp, 0.133_wp, 0.010_wp,& ! 3
680	3.67_wp, -1.977_wp, 1.10_wp, 2.87E-6_wp, 0.520_wp, 0.448_wp, 0.279_wp, 0.010_wp,& ! 4
681	2.65_wp, 2.500_wp, 1.10_wp, 1.74E-6_wp, 0.614_wp, 0.541_wp, 0.335_wp, 0.010_wp,& ! 5
682	1.30_wp, 0.400_wp, 1.20_wp, 0.93E-6_wp, 0.766_wp, 0.663_wp, 0.267_wp, 0.010_wp & ! 6
683	/), (/ 8, 6 /) )
684
685	!
686	!-- TO BE FILLED
687	!-- Pavement parameters z0, z0h, albedo_type, emissivity
688	REAL(wp), DIMENSION(0:3,1:15), PARAMETER :: pavement_pars = RESHAPE( (/ &
689	1.0E-4_wp, 1.0E-5_wp, 18.0_wp, 0.97_wp, & ! 1
690	1.0E-4_wp, 1.0E-5_wp, 19.0_wp, 0.94_wp, & ! 2
691	1.0E-4_wp, 1.0E-5_wp, 20.0_wp, 0.98_wp, & ! 3
692	1.0E-4_wp, 1.0E-5_wp, 21.0_wp, 0.93_wp, & ! 4
693	1.0E-4_wp, 1.0E-5_wp, 22.0_wp, 0.97_wp, & ! 5
694	1.0E-4_wp, 1.0E-5_wp, 23.0_wp, 0.97_wp, & ! 6
695	1.0E-4_wp, 1.0E-5_wp, 24.0_wp, 0.97_wp, & ! 7
696	1.0E-4_wp, 1.0E-5_wp, 25.0_wp, 0.94_wp, & ! 8
697	1.0E-4_wp, 1.0E-5_wp, 26.0_wp, 0.98_wp, & ! 9
698	1.0E-4_wp, 1.0E-5_wp, 27.0_wp, 0.93_wp, & ! 10
699	1.0E-4_wp, 1.0E-5_wp, 28.0_wp, 0.97_wp, & ! 11
700	1.0E-4_wp, 1.0E-5_wp, 29.0_wp, 0.97_wp, & ! 12
701	1.0E-4_wp, 1.0E-5_wp, 30.0_wp, 0.97_wp, & ! 13
702	1.0E-4_wp, 1.0E-5_wp, 31.0_wp, 0.94_wp, & ! 14
703	1.0E-4_wp, 1.0E-5_wp, 32.0_wp, 0.98_wp & ! 15
704	/), (/ 4, 15 /) )
705	!
706	!-- Pavement subsurface parameters part 1: thermal conductivity (W/m/K)
707	!-- 0.0-0.01, 0.01-0.03, 0.03-0.07, 0.07-0.15, 0.15-0.30, 0.30-0.50, 0.50-1.25, 1.25-3.00
708	REAL(wp), DIMENSION(0:7,1:15), PARAMETER :: pavement_subsurface_pars_1 = RESHAPE( (/ &
709	1.00_wp, 1.00_wp, 1.00_wp, 1.00_wp, 1.00_wp, 1.00_wp, 9999999.9_wp, 9999999.9_wp, & ! 1
710	1.00_wp, 1.00_wp, 1.00_wp, 1.00_wp, 1.00_wp, 1.00_wp, 9999999.9_wp, 9999999.9_wp, & ! 2
711	1.00_wp, 1.00_wp, 1.00_wp, 1.00_wp, 1.00_wp, 1.00_wp, 9999999.9_wp, 9999999.9_wp, & ! 3
712	1.00_wp, 1.00_wp, 1.00_wp, 1.00_wp, 1.00_wp, 1.00_wp, 9999999.9_wp, 9999999.9_wp, & ! 4
713	1.00_wp, 1.00_wp, 1.00_wp, 1.00_wp, 1.00_wp, 1.00_wp, 9999999.9_wp, 9999999.9_wp, & ! 5
714	1.00_wp, 1.00_wp, 1.00_wp, 1.00_wp, 1.00_wp, 1.00_wp, 9999999.9_wp, 9999999.9_wp, & ! 6
715	1.00_wp, 1.00_wp, 1.00_wp, 1.00_wp, 1.00_wp, 1.00_wp, 9999999.9_wp, 9999999.9_wp, & ! 7
716	1.00_wp, 1.00_wp, 1.00_wp, 1.00_wp, 1.00_wp, 1.00_wp, 9999999.9_wp, 9999999.9_wp, & ! 8
717	1.00_wp, 1.00_wp, 1.00_wp, 1.00_wp, 1.00_wp, 1.00_wp, 9999999.9_wp, 9999999.9_wp, & ! 9
718	1.00_wp, 1.00_wp, 1.00_wp, 1.00_wp, 1.00_wp, 1.00_wp, 9999999.9_wp, 9999999.9_wp, & ! 10
719	1.00_wp, 1.00_wp, 1.00_wp, 1.00_wp, 1.00_wp, 1.00_wp, 9999999.9_wp, 9999999.9_wp, & ! 11
720	1.00_wp, 1.00_wp, 1.00_wp, 1.00_wp, 1.00_wp, 1.00_wp, 9999999.9_wp, 9999999.9_wp, & ! 12
721	1.00_wp, 1.00_wp, 1.00_wp, 1.00_wp, 1.00_wp, 1.00_wp, 9999999.9_wp, 9999999.9_wp, & ! 13
722	1.00_wp, 1.00_wp, 1.00_wp, 1.00_wp, 1.00_wp, 1.00_wp, 9999999.9_wp, 9999999.9_wp, & ! 14
723	1.00_wp, 1.00_wp, 1.00_wp, 1.00_wp, 1.00_wp, 1.00_wp, 9999999.9_wp, 9999999.9_wp & ! 15
724	/), (/ 8, 15 /) )
725
726	!
727	!-- Pavement subsurface parameters part 2: volumetric heat capacity (J/m3/K)
728	!-- 0.0-0.01, 0.01-0.03, 0.03-0.07, 0.07-0.15, 0.15-0.30, 0.30-0.50, 0.50-1.25, 1.25-3.00
729	REAL(wp), DIMENSION(0:7,1:15), PARAMETER :: pavement_subsurface_pars_2 = RESHAPE( (/ &
730	1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 9999999.9_wp, 9999999.9_wp, & ! 1
731	1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 9999999.9_wp, 9999999.9_wp, & ! 2
732	1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 9999999.9_wp, 9999999.9_wp, & ! 3
733	1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 9999999.9_wp, 9999999.9_wp, & ! 4
734	1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 9999999.9_wp, 9999999.9_wp, & ! 5
735	1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 9999999.9_wp, 9999999.9_wp, & ! 6
736	1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 9999999.9_wp, 9999999.9_wp, & ! 7
737	1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 9999999.9_wp, 9999999.9_wp, & ! 8
738	1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 9999999.9_wp, 9999999.9_wp, & ! 9
739	1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 9999999.9_wp, 9999999.9_wp, & ! 10
740	1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 9999999.9_wp, 9999999.9_wp, & ! 11
741	1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 9999999.9_wp, 9999999.9_wp, & ! 12
742	1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 9999999.9_wp, 9999999.9_wp, & ! 13
743	1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 9999999.9_wp, 9999999.9_wp, & ! 14
744	1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 9999999.9_wp, 9999999.9_wp & ! 15
745	/), (/ 8, 15 /) )
746
747	!
748	!-- TO BE FILLED
749	!-- Water parameters temperature, z0, z0h, albedo_type, emissivity,
750	REAL(wp), DIMENSION(0:4,1:5), PARAMETER :: water_pars = RESHAPE( (/ &
751	283.0_wp, 0.01_wp, 0.001_wp, 1.0_wp, 0.99_wp, & ! 1
752	283.0_wp, 0.01_wp, 0.001_wp, 1.0_wp, 0.99_wp, & ! 2
753	283.0_wp, 0.01_wp, 0.001_wp, 1.0_wp, 0.99_wp, & ! 3
754	283.0_wp, 0.01_wp, 0.001_wp, 1.0_wp, 0.99_wp, & ! 4
755	283.0_wp, 0.01_wp, 0.001_wp, 1.0_wp, 0.99_wp & ! 5
756	/), (/ 5, 5 /) )
757
758	SAVE
759
760
761	PRIVATE
762
763
764	!
765	!-- Public functions
766	PUBLIC lsm_check_data_output, lsm_check_data_output_pr, &
767	lsm_check_parameters, lsm_define_netcdf_grid, lsm_3d_data_averaging,&
768	lsm_data_output_2d, lsm_data_output_3d, lsm_energy_balance, &
769	lsm_header, lsm_init, lsm_init_arrays, lsm_parin, lsm_soil_model, &
770	lsm_swap_timelevel, lsm_read_restart_data, lsm_last_actions
771	! !vegetat
772	!-- Public parameters, constants and initial values
773	PUBLIC aero_resist_kray, skip_time_do_lsm
774
775	!
776	!-- Public grid variables
777	PUBLIC nzb_soil, nzs, nzt_soil, zs
778
779	!
780	!-- Public prognostic variables
781	PUBLIC m_soil_h, t_soil_h
782
783
784	INTERFACE lsm_check_data_output
785	MODULE PROCEDURE lsm_check_data_output
786	END INTERFACE lsm_check_data_output
787
788	INTERFACE lsm_check_data_output_pr
789	MODULE PROCEDURE lsm_check_data_output_pr
790	END INTERFACE lsm_check_data_output_pr
791
792	INTERFACE lsm_check_parameters
793	MODULE PROCEDURE lsm_check_parameters
794	END INTERFACE lsm_check_parameters
795
796	INTERFACE lsm_3d_data_averaging
797	MODULE PROCEDURE lsm_3d_data_averaging
798	END INTERFACE lsm_3d_data_averaging
799
800	INTERFACE lsm_data_output_2d
801	MODULE PROCEDURE lsm_data_output_2d
802	END INTERFACE lsm_data_output_2d
803
804	INTERFACE lsm_data_output_3d
805	MODULE PROCEDURE lsm_data_output_3d
806	END INTERFACE lsm_data_output_3d
807
808	INTERFACE lsm_define_netcdf_grid
809	MODULE PROCEDURE lsm_define_netcdf_grid
810	END INTERFACE lsm_define_netcdf_grid
811
812	INTERFACE lsm_energy_balance
813	MODULE PROCEDURE lsm_energy_balance
814	END INTERFACE lsm_energy_balance
815
816	INTERFACE lsm_header
817	MODULE PROCEDURE lsm_header
818	END INTERFACE lsm_header
819
820	INTERFACE lsm_init
821	MODULE PROCEDURE lsm_init
822	END INTERFACE lsm_init
823
824	INTERFACE lsm_init_arrays
825	MODULE PROCEDURE lsm_init_arrays
826	END INTERFACE lsm_init_arrays
827
828	INTERFACE lsm_parin
829	MODULE PROCEDURE lsm_parin
830	END INTERFACE lsm_parin
831
832	INTERFACE lsm_soil_model
833	MODULE PROCEDURE lsm_soil_model
834	END INTERFACE lsm_soil_model
835
836	INTERFACE lsm_swap_timelevel
837	MODULE PROCEDURE lsm_swap_timelevel
838	END INTERFACE lsm_swap_timelevel
839
840	INTERFACE lsm_read_restart_data
841	MODULE PROCEDURE lsm_read_restart_data
842	END INTERFACE lsm_read_restart_data
843
844	INTERFACE lsm_last_actions
845	MODULE PROCEDURE lsm_last_actions
846	END INTERFACE lsm_last_actions
847
848	CONTAINS
849
850	!------------------------------------------------------------------------------!
851	! Description:
852	! ------------
853	!> Check data output for land surface model
854	!------------------------------------------------------------------------------!
855	SUBROUTINE lsm_check_data_output( var, unit, i, ilen, k )
856
857
858	USE control_parameters, &
859	ONLY: data_output, message_string
860
861	IMPLICIT NONE
862
863	CHARACTER (LEN=*) :: unit !<
864	CHARACTER (LEN=*) :: var !<
865
866	INTEGER(iwp) :: i
867	INTEGER(iwp) :: ilen
868	INTEGER(iwp) :: k
869
870	SELECT CASE ( TRIM( var ) )
871
872	CASE ( 'm_soil' )
873	IF ( .NOT. land_surface ) THEN
874	message_string = 'output of "' // TRIM( var ) // '" requi' // &
875	'res land_surface = .TRUE.'
876	CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
877	ENDIF
878	unit = 'm3/m3'
879
880	CASE ( 't_soil' )
881	IF ( .NOT. land_surface ) THEN
882	message_string = 'output of "' // TRIM( var ) // '" requi' // &
883	'res land_surface = .TRUE.'
884	CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
885	ENDIF
886	unit = 'K'
887
888	CASE ( 'lai', 'c_liq', 'c_soil', 'c_veg', 'ghf', 'm_liq', &
889	'qsws_liq', 'qsws_soil', 'qsws_veg*', &
890	'r_a', 'r_s' )
891	IF ( k == 0 .OR. data_output(i)(ilen-2:ilen) /= '_xy' ) THEN
892	message_string = 'illegal value for data_output: "' // &
893	TRIM( var ) // '" & only 2d-horizontal ' // &
894	'cross sections are allowed for this value'
895	CALL message( 'check_parameters', 'PA0111', 1, 2, 0, 6, 0 )
896	ENDIF
897	IF ( TRIM( var ) == 'lai*' .AND. .NOT. land_surface ) THEN
898	message_string = 'output of "' // TRIM( var ) // '" requi' // &
899	'res land_surface = .TRUE.'
900	CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
901	ENDIF
902	IF ( TRIM( var ) == 'c_liq*' .AND. .NOT. land_surface ) THEN
903	message_string = 'output of "' // TRIM( var ) // '" requi' // &
904	'res land_surface = .TRUE.'
905	CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
906	ENDIF
907	IF ( TRIM( var ) == 'c_soil*' .AND. .NOT. land_surface ) THEN
908	message_string = 'output of "' // TRIM( var ) // '" requi' // &
909	'res land_surface = .TRUE.'
910	CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
911	ENDIF
912	IF ( TRIM( var ) == 'c_veg*' .AND. .NOT. land_surface ) THEN
913	message_string = 'output of "' // TRIM( var ) // '" requi' // &
914	'res land_surface = .TRUE.'
915	CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
916	ENDIF
917	IF ( TRIM( var ) == 'ghf*' .AND. .NOT. land_surface ) THEN
918	message_string = 'output of "' // TRIM( var ) // '" requi' // &
919	'res land_surface = .TRUE.'
920	CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
921	ENDIF
922	IF ( TRIM( var ) == 'm_liq*' .AND. .NOT. land_surface ) THEN
923	message_string = 'output of "' // TRIM( var ) // '" requi' // &
924	'res land_surface = .TRUE.'
925	CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
926	ENDIF
927	IF ( TRIM( var ) == 'qsws_liq*' .AND. .NOT. land_surface ) &
928	THEN
929	message_string = 'output of "' // TRIM( var ) // '" requi' // &
930	'res land_surface = .TRUE.'
931	CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
932	ENDIF
933	IF ( TRIM( var ) == 'qsws_soil*' .AND. .NOT. land_surface ) &
934	THEN
935	message_string = 'output of "' // TRIM( var ) // '" requi' // &
936	'res land_surface = .TRUE.'
937	CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
938	ENDIF
939	IF ( TRIM( var ) == 'qsws_veg*' .AND. .NOT. land_surface ) &
940	THEN
941	message_string = 'output of "' // TRIM( var ) // '" requi' // &
942	'res land_surface = .TRUE.'
943	CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
944	ENDIF
945	IF ( TRIM( var ) == 'r_a*' .AND. .NOT. land_surface ) &
946	THEN
947	message_string = 'output of "' // TRIM( var ) // '" requi' // &
948	'res land_surface = .TRUE.'
949	CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
950	ENDIF
951	IF ( TRIM( var ) == 'r_s*' .AND. .NOT. land_surface ) &
952	THEN
953	message_string = 'output of "' // TRIM( var ) // '" requi' // &
954	'res land_surface = .TRUE.'
955	CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
956	ENDIF
957
958	IF ( TRIM( var ) == 'lai*' ) unit = 'none'
959	IF ( TRIM( var ) == 'c_liq*' ) unit = 'none'
960	IF ( TRIM( var ) == 'c_soil*') unit = 'none'
961	IF ( TRIM( var ) == 'c_veg*' ) unit = 'none'
962	IF ( TRIM( var ) == 'ghf*') unit = 'W/m2'
963	IF ( TRIM( var ) == 'm_liq*' ) unit = 'm'
964	IF ( TRIM( var ) == 'qsws_liq*' ) unit = 'W/m2'
965	IF ( TRIM( var ) == 'qsws_soil*' ) unit = 'W/m2'
966	IF ( TRIM( var ) == 'qsws_veg*' ) unit = 'W/m2'
967	IF ( TRIM( var ) == 'r_a*') unit = 's/m'
968	IF ( TRIM( var ) == 'r_s*') unit = 's/m'
969
970	CASE DEFAULT
971	unit = 'illegal'
972
973	END SELECT
974
975
976	END SUBROUTINE lsm_check_data_output
977
978
979	!------------------------------------------------------------------------------!
980	! Description:
981	! ------------
982	!> Check data output of profiles for land surface model
983	!------------------------------------------------------------------------------!
984	SUBROUTINE lsm_check_data_output_pr( variable, var_count, unit, dopr_unit )
985
986	USE control_parameters, &
987	ONLY: data_output_pr, message_string
988
989	USE indices
990
991	USE profil_parameter
992
993	USE statistics
994
995	IMPLICIT NONE
996
997	CHARACTER (LEN=*) :: unit !<
998	CHARACTER (LEN=*) :: variable !<
999	CHARACTER (LEN=*) :: dopr_unit !< local value of dopr_unit
1000
1001	INTEGER(iwp) :: user_pr_index !<
1002	INTEGER(iwp) :: var_count !<
1003
1004	SELECT CASE ( TRIM( variable ) )
1005
1006	CASE ( 't_soil', '#t_soil' )
1007	IF ( .NOT. land_surface ) THEN
1008	message_string = 'data_output_pr = ' // &
1009	TRIM( data_output_pr(var_count) ) // ' is' // &
1010	'not implemented for land_surface = .FALSE.'
1011	CALL message( 'check_parameters', 'PA0402', 1, 2, 0, 6, 0 )
1012	ELSE
1013	dopr_index(var_count) = 89
1014	dopr_unit = 'K'
1015	hom(0:nzs-1,2,89,:) = SPREAD( - zs(nzb_soil:nzt_soil), 2, statistic_regions+1 )
1016	IF ( data_output_pr(var_count)(1:1) == '#' ) THEN
1017	dopr_initial_index(var_count) = 90
1018	hom(0:nzs-1,2,90,:) = SPREAD( - zs(nzb_soil:nzt_soil), 2, statistic_regions+1 )
1019	data_output_pr(var_count) = data_output_pr(var_count)(2:)
1020	ENDIF
1021	unit = dopr_unit
1022	ENDIF
1023
1024	CASE ( 'm_soil', '#m_soil' )
1025	IF ( .NOT. land_surface ) THEN
1026	message_string = 'data_output_pr = ' // &
1027	TRIM( data_output_pr(var_count) ) // ' is' // &
1028	' not implemented for land_surface = .FALSE.'
1029	CALL message( 'check_parameters', 'PA0402', 1, 2, 0, 6, 0 )
1030	ELSE
1031	dopr_index(var_count) = 91
1032	dopr_unit = 'm3/m3'
1033	hom(0:nzs-1,2,91,:) = SPREAD( - zs(nzb_soil:nzt_soil), 2, statistic_regions+1 )
1034	IF ( data_output_pr(var_count)(1:1) == '#' ) THEN
1035	dopr_initial_index(var_count) = 92
1036	hom(0:nzs-1,2,92,:) = SPREAD( - zs(nzb_soil:nzt_soil), 2, statistic_regions+1 )
1037	data_output_pr(var_count) = data_output_pr(var_count)(2:)
1038	ENDIF
1039	unit = dopr_unit
1040	ENDIF
1041
1042
1043	CASE DEFAULT
1044	unit = 'illegal'
1045
1046	END SELECT
1047
1048
1049	END SUBROUTINE lsm_check_data_output_pr
1050
1051
1052	!------------------------------------------------------------------------------!
1053	! Description:
1054	! ------------
1055	!> Check parameters routine for land surface model
1056	!------------------------------------------------------------------------------!
1057	SUBROUTINE lsm_check_parameters
1058
1059	USE control_parameters, &
1060	ONLY: bc_pt_b, bc_q_b, constant_flux_layer, message_string, &
1061	most_method
1062
1063	USE radiation_model_mod, &
1064	ONLY: radiation
1065
1066
1067	IMPLICIT NONE
1068
1069	INTEGER(iwp) :: k !< running index, z-dimension
1070
1071	!
1072	!-- Check for a valid setting of surface_type. The default value is 'netcdf'.
1073	!-- In that case, the surface types are read from NetCDF file
1074	IF ( TRIM( surface_type ) /= 'vegetation' .AND. &
1075	TRIM( surface_type ) /= 'pavement' .AND. &
1076	TRIM( surface_type ) /= 'water' .AND. &
1077	TRIM( surface_type ) /= 'netcdf' ) THEN
1078	message_string = 'unknown surface type surface_type = "' // &
1079	TRIM( surface_type ) // '"'
1080	CALL message( 'check_parameters', 'PA0019', 1, 2, 0, 6, 0 )
1081	ENDIF
1082
1083	!
1084	!-- Dirichlet boundary conditions are required as the surface fluxes are
1085	!-- calculated from the temperature/humidity gradients in the land surface
1086	!-- model
1087	IF ( bc_pt_b == 'neumann' .OR. bc_q_b == 'neumann' ) THEN
1088	message_string = 'lsm requires setting of'// &
1089	'bc_pt_b = "dirichlet" and '// &
1090	'bc_q_b = "dirichlet"'
1091	CALL message( 'check_parameters', 'PA0399', 1, 2, 0, 6, 0 )
1092	ENDIF
1093
1094	IF ( .NOT. constant_flux_layer ) THEN
1095	message_string = 'lsm requires '// &
1096	'constant_flux_layer = .T.'
1097	CALL message( 'check_parameters', 'PA0400', 1, 2, 0, 6, 0 )
1098	ENDIF
1099
1100	IF ( TRIM( surface_type ) == 'vegetation' ) THEN
1101
1102	IF ( vegetation_type == 0 ) THEN
1103	IF ( min_canopy_resistance == 9999999.9_wp ) THEN
1104	message_string = 'vegetation_type = 0 (user defined)'// &
1105	'requires setting of min_canopy_resistance'// &
1106	'/= 9999999.9'
1107	CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
1108	ENDIF
1109
1110	IF ( leaf_area_index == 9999999.9_wp ) THEN
1111	message_string = 'vegetation_type = 0 (user_defined)'// &
1112	'requires setting of leaf_area_index'// &
1113	'/= 9999999.9'
1114	CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
1115	ENDIF
1116
1117	IF ( vegetation_coverage == 9999999.9_wp ) THEN
1118	message_string = 'vegetation_type = 0 (user_defined)'// &
1119	'requires setting of vegetation_coverage'// &
1120	'/= 9999999.9'
1121	CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
1122	ENDIF
1123
1124	IF ( canopy_resistance_coefficient == 9999999.9_wp) THEN
1125	message_string = 'vegetation_type = 0 (user_defined)'// &
1126	'requires setting of'// &
1127	'canopy_resistance_coefficient /= 9999999.9'
1128	CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
1129	ENDIF
1130
1131	IF ( lambda_surface_stable == 9999999.9_wp ) THEN
1132	message_string = 'vegetation_type = 0 (user_defined)'// &
1133	'requires setting of lambda_surface_stable'// &
1134	'/= 9999999.9'
1135	CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
1136	ENDIF
1137
1138	IF ( lambda_surface_unstable == 9999999.9_wp ) THEN
1139	message_string = 'vegetation_type = 0 (user_defined)'// &
1140	'requires setting of lambda_surface_unstable'// &
1141	'/= 9999999.9'
1142	CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
1143	ENDIF
1144
1145	IF ( f_shortwave_incoming == 9999999.9_wp ) THEN
1146	message_string = 'vegetation_type = 0 (user_defined)'// &
1147	'requires setting of f_shortwave_incoming'// &
1148	'/= 9999999.9'
1149	CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
1150	ENDIF
1151
1152	IF ( z0_vegetation == 9999999.9_wp ) THEN
1153	message_string = 'vegetation_type = 0 (user_defined)'// &
1154	'requires setting of z0_vegetation'// &
1155	'/= 9999999.9'
1156	CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
1157	ENDIF
1158
1159	IF ( z0h_vegetation == 9999999.9_wp ) THEN
1160	message_string = 'vegetation_type = 0 (user_defined)'// &
1161	'requires setting of z0h_vegetation'// &
1162	'/= 9999999.9'
1163	CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
1164	ENDIF
1165	ENDIF
1166
1167	IF ( vegetation_type == 1 ) THEN
1168	IF ( vegetation_coverage /= 9999999.9_wp .AND. vegetation_coverage &
1169	/= 0.0_wp ) THEN
1170	message_string = 'vegetation_type = 1 (bare soil)'// &
1171	' requires vegetation_coverage = 0'
1172	CALL message( 'check_parameters', 'PA0471', 1, 2, 0, 6, 0 )
1173	ENDIF
1174	ENDIF
1175
1176	ENDIF
1177
1178	IF ( TRIM( surface_type ) == 'water' ) THEN
1179
1180	IF ( TRIM( most_method ) == 'lookup' ) THEN
1181	WRITE( message_string, * ) 'surface_type = ', surface_type, &
1182	' is not allowed in combination with ', &
1183	'most_method = ', most_method
1184	CALL message( 'check_parameters', 'PA0417', 1, 2, 0, 6, 0 )
1185	ENDIF
1186
1187	IF ( water_type == 0 ) THEN
1188
1189	IF ( z0_water == 9999999.9_wp ) THEN
1190	message_string = 'water_type = 0 (user_defined)'// &
1191	'requires setting of z0_water'// &
1192	'/= 9999999.9'
1193	CALL message( 'check_parameters', 'PA0415', 1, 2, 0, 6, 0 )
1194	ENDIF
1195
1196	IF ( z0h_water == 9999999.9_wp ) THEN
1197	message_string = 'water_type = 0 (user_defined)'// &
1198	'requires setting of z0h_water'// &
1199	'/= 9999999.9'
1200	CALL message( 'check_parameters', 'PA0392', 1, 2, 0, 6, 0 )
1201	ENDIF
1202
1203	IF ( water_temperature == 9999999.9_wp ) THEN
1204	message_string = 'water_type = 0 (user_defined)'// &
1205	'requires setting of water_temperature'// &
1206	'/= 9999999.9'
1207	CALL message( 'check_parameters', 'PA0379', 1, 2, 0, 6, 0 )
1208	ENDIF
1209
1210	ENDIF
1211
1212	ENDIF
1213
1214	IF ( TRIM( surface_type ) == 'pavement' ) THEN
1215
1216	IF ( ANY( dz_soil /= 9999999.9_wp ) .AND. pavement_type /= 0 ) THEN
1217	message_string = 'non-default setting of dz_soil '// &
1218	'does not allow to use pavement_type /= 0)'
1219	CALL message( 'check_parameters', 'PA0341', 1, 2, 0, 6, 0 )
1220	ENDIF
1221
1222	IF ( pavement_type == 0 ) THEN
1223
1224	IF ( z0_pavement == 9999999.9_wp ) THEN
1225	message_string = 'pavement_type = 0 (user_defined)'// &
1226	'requires setting of z0_pavement'// &
1227	'/= 9999999.9'
1228	CALL message( 'check_parameters', 'PA0352', 1, 2, 0, 6, 0 )
1229	ENDIF
1230
1231	IF ( z0h_pavement == 9999999.9_wp ) THEN
1232	message_string = 'pavement_type = 0 (user_defined)'// &
1233	'requires setting of z0h_pavement'// &
1234	'/= 9999999.9'
1235	CALL message( 'check_parameters', 'PA0353', 1, 2, 0, 6, 0 )
1236	ENDIF
1237
1238	IF ( pavement_heat_conduct == 9999999.9_wp ) THEN
1239	message_string = 'pavement_type = 0 (user_defined)'// &
1240	'requires setting of pavement_heat_conduct'// &
1241	'/= 9999999.9'
1242	CALL message( 'check_parameters', 'PA0342', 1, 2, 0, 6, 0 )
1243	ENDIF
1244
1245	IF ( pavement_heat_capacity == 9999999.9_wp ) THEN
1246	message_string = 'pavement_type = 0 (user_defined)'// &
1247	'requires setting of pavement_heat_capacity'// &
1248	'/= 9999999.9'
1249	CALL message( 'check_parameters', 'PA0139', 1, 2, 0, 6, 0 )
1250	ENDIF
1251
1252	IF ( pavement_depth_level == 0 ) THEN
1253	message_string = 'pavement_type = 0 (user_defined)'// &
1254	'requires setting of pavement_depth_level'// &
1255	'/= 0'
1256	CALL message( 'check_parameters', 'PA0474', 1, 2, 0, 6, 0 )
1257	ENDIF
1258
1259	ENDIF
1260
1261	ENDIF
1262
1263	!
1264	!-- Temporary message as long as NetCDF input is not available
1265	IF ( TRIM( surface_type ) == 'netcdf' ) THEN
1266	message_string = 'surface_type = netcdf '// &
1267	'is not supported at the moment'
1268	CALL message( 'check_parameters', 'PA0465', 1, 2, 0, 6, 0 )
1269	ENDIF
1270
1271	IF ( soil_type == 0 ) THEN
1272
1273	IF ( alpha_vangenuchten == 9999999.9_wp ) THEN
1274	message_string = 'soil_type = 0 (user_defined)'// &
1275	'requires setting of alpha_vangenuchten'// &
1276	'/= 9999999.9'
1277	CALL message( 'check_parameters', 'PA0403', 1, 2, 0, 6, 0 )
1278	ENDIF
1279
1280	IF ( l_vangenuchten == 9999999.9_wp ) THEN
1281	message_string = 'soil_type = 0 (user_defined)'// &
1282	'requires setting of l_vangenuchten'// &
1283	'/= 9999999.9'
1284	CALL message( 'check_parameters', 'PA0403', 1, 2, 0, 6, 0 )
1285	ENDIF
1286
1287	IF ( n_vangenuchten == 9999999.9_wp ) THEN
1288	message_string = 'soil_type = 0 (user_defined)'// &
1289	'requires setting of n_vangenuchten'// &
1290	'/= 9999999.9'
1291	CALL message( 'check_parameters', 'PA0403', 1, 2, 0, 6, 0 )
1292	ENDIF
1293
1294	IF ( hydraulic_conductivity == 9999999.9_wp ) THEN
1295	message_string = 'soil_type = 0 (user_defined)'// &
1296	'requires setting of hydraulic_conductivity'// &
1297	'/= 9999999.9'
1298	CALL message( 'check_parameters', 'PA0403', 1, 2, 0, 6, 0 )
1299	ENDIF
1300
1301	IF ( saturation_moisture == 9999999.9_wp ) THEN
1302	message_string = 'soil_type = 0 (user_defined)'// &
1303	'requires setting of saturation_moisture'// &
1304	'/= 9999999.9'
1305	CALL message( 'check_parameters', 'PA0403', 1, 2, 0, 6, 0 )
1306	ENDIF
1307
1308	IF ( field_capacity == 9999999.9_wp ) THEN
1309	message_string = 'soil_type = 0 (user_defined)'// &
1310	'requires setting of field_capacity'// &
1311	'/= 9999999.9'
1312	CALL message( 'check_parameters', 'PA0403', 1, 2, 0, 6, 0 )
1313	ENDIF
1314
1315	IF ( wilting_point == 9999999.9_wp ) THEN
1316	message_string = 'soil_type = 0 (user_defined)'// &
1317	'requires setting of wilting_point'// &
1318	'/= 9999999.9'
1319	CALL message( 'check_parameters', 'PA0403', 1, 2, 0, 6, 0 )
1320	ENDIF
1321
1322	IF ( residual_moisture == 9999999.9_wp ) THEN
1323	message_string = 'soil_type = 0 (user_defined)'// &
1324	'requires setting of residual_moisture'// &
1325	'/= 9999999.9'
1326	CALL message( 'check_parameters', 'PA0403', 1, 2, 0, 6, 0 )
1327	ENDIF
1328
1329	ENDIF
1330
1331	!
1332	!-- Determine number of soil layers to be used and check whether an appropriate
1333	!-- root fraction is prescribed
1334	nzb_soil = 0
1335	nzt_soil = -1
1336	IF ( ALL( dz_soil == 9999999.9_wp ) ) THEN
1337	nzt_soil = 7
1338	dz_soil(nzb_soil:nzt_soil) = dz_soil_default
1339	ELSE
1340	DO k = 0, 19
1341	IF ( dz_soil(k) /= 9999999.9_wp ) THEN
1342	nzt_soil = nzt_soil + 1
1343	ENDIF
1344	ENDDO
1345	ENDIF
1346	nzs = nzt_soil + 1
1347
1348
1349	IF ( vegetation_type == 0 ) THEN
1350	IF ( SUM( root_fraction(nzb_soil:nzt_soil) ) /= 1.0_wp ) THEN
1351	message_string = 'vegetation_type = 0 (user_defined)'// &
1352	'requires setting of root_fraction'// &
1353	'/= 9999999.9 and SUM(root_fraction) = 1'
1354	CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
1355	ENDIF
1356	ENDIF
1357
1358
1359	!
1360	!-- Check for proper setting of soil moisture, must not be larger than its
1361	!-- saturation value.
1362	DO k = nzb_soil, nzt_soil
1363	IF ( soil_moisture(k) > saturation_moisture ) THEN
1364	message_string = 'soil_moisture must not exceed its saturation' // &
1365	' value'
1366	CALL message( 'check_parameters', 'PA0458', 1, 2, 0, 6, 0 )
1367	ENDIF
1368	ENDDO
1369
1370	!
1371	!-- Calculate grid spacings. Temperature and moisture are defined at
1372	!-- the center of the soil layers, whereas gradients/fluxes are
1373	!-- defined at the edges (_layer)
1374	!
1375	!-- Allocate global 1D arrays
1376	ALLOCATE ( ddz_soil_center(nzb_soil:nzt_soil) )
1377	ALLOCATE ( ddz_soil(nzb_soil:nzt_soil+1) )
1378	ALLOCATE ( dz_soil_center(nzb_soil:nzt_soil) )
1379	ALLOCATE ( zs(nzb_soil:nzt_soil+1) )
1380
1381	zs(nzb_soil) = 0.5_wp * dz_soil(nzb_soil)
1382	zs_layer(nzb_soil) = dz_soil(nzb_soil)
1383
1384	DO k = nzb_soil+1, nzt_soil
1385	zs_layer(k) = zs_layer(k-1) + dz_soil(k)
1386	zs(k) = (zs_layer(k) + zs_layer(k-1)) * 0.5_wp
1387	ENDDO
1388
1389	dz_soil(nzt_soil+1) = zs_layer(nzt_soil) + dz_soil(nzt_soil)
1390	zs(nzt_soil+1) = zs_layer(nzt_soil) + 0.5_wp * dz_soil(nzt_soil)
1391
1392	DO k = nzb_soil, nzt_soil-1
1393	dz_soil_center(k) = zs(k+1) - zs(k)
1394	IF ( dz_soil_center(k) == 0.0_wp ) THEN
1395	message_string = 'invalid soil layer configuration found ' // &
1396	'(dz_soil_center(k) = 0.0)'
1397	CALL message( 'lsm_read_restart_data', 'PA0140', 1, 2, 0, 6, 0 )
1398	ENDIF
1399	ENDDO
1400
1401	dz_soil_center(nzt_soil) = zs_layer(k-1) + dz_soil(k) - zs(nzt_soil)
1402
1403	ddz_soil_center = 1.0_wp / dz_soil_center
1404	ddz_soil(nzb_soil:nzt_soil) = 1.0_wp / dz_soil(nzb_soil:nzt_soil)
1405
1406
1407
1408	END SUBROUTINE lsm_check_parameters
1409
1410	!------------------------------------------------------------------------------!
1411	! Description:
1412	! ------------
1413	!> Solver for the energy balance at the surface.
1414	!------------------------------------------------------------------------------!
1415	SUBROUTINE lsm_energy_balance( horizontal, l )
1416
1417	USE diagnostic_quantities_mod, &
1418	ONLY: magnus
1419
1420	USE pegrid
1421
1422	IMPLICIT NONE
1423
1424	INTEGER(iwp) :: i !< running index
1425	INTEGER(iwp) :: i_off !< offset to determine index of surface element, seen from atmospheric grid point, for x
1426	INTEGER(iwp) :: j !< running index
1427	INTEGER(iwp) :: j_off !< offset to determine index of surface element, seen from atmospheric grid point, for y
1428	INTEGER(iwp) :: k !< running index
1429	INTEGER(iwp) :: k_off !< offset to determine index of surface element, seen from atmospheric grid point, for z
1430	INTEGER(iwp) :: k_rad !< index to access radiation array
1431	INTEGER(iwp) :: ks !< running index
1432	INTEGER(iwp) :: l !< surface-facing index
1433	INTEGER(iwp) :: m !< running index concerning wall elements
1434
1435	LOGICAL :: horizontal !< Flag indicating horizontal or vertical surfaces
1436
1437	REAL(wp) :: c_surface_tmp,& !< temporary variable for storing the volumetric heat capacity of the surface
1438	f1, & !< resistance correction term 1
1439	f2, & !< resistance correction term 2
1440	f3, & !< resistance correction term 3
1441	m_min, & !< minimum soil moisture
1442	e, & !< water vapour pressure
1443	e_s, & !< water vapour saturation pressure
1444	e_s_dt, & !< derivate of e_s with respect to T
1445	tend, & !< tendency
1446	dq_s_dt, & !< derivate of q_s with respect to T
1447	coef_1, & !< coef. for prognostic equation
1448	coef_2, & !< coef. for prognostic equation
1449	f_qsws, & !< factor for qsws
1450	f_qsws_veg, & !< factor for qsws_veg
1451	f_qsws_soil, & !< factor for qsws_soil
1452	f_qsws_liq, & !< factor for qsws_liq
1453	f_shf, & !< factor for shf
1454	lambda_soil, & !< Thermal conductivity of the uppermost soil layer (W/m2/K)
1455	lambda_surface, & !< Current value of lambda_surface (W/m2/K)
1456	m_liq_max, & !< maxmimum value of the liq. water reservoir
1457	pt1, & !< potential temperature at first grid level
1458	qv1 !< specific humidity at first grid level
1459
1460	TYPE(surf_type_lsm), POINTER :: surf_t_surface
1461	TYPE(surf_type_lsm), POINTER :: surf_t_surface_p
1462	TYPE(surf_type_lsm), POINTER :: surf_tt_surface_m
1463	TYPE(surf_type_lsm), POINTER :: surf_m_liq
1464	TYPE(surf_type_lsm), POINTER :: surf_m_liq_p
1465	TYPE(surf_type_lsm), POINTER :: surf_tm_liq_m
1466
1467	TYPE(surf_type_lsm), POINTER :: surf_m_soil
1468	TYPE(surf_type_lsm), POINTER :: surf_t_soil
1469
1470	TYPE(surf_type), POINTER :: surf !< surface-date type variable
1471
1472	IF ( horizontal ) THEN
1473	surf => surf_lsm_h
1474	surf_t_surface => t_surface_h
1475	surf_t_surface_p => t_surface_h_p
1476	surf_tt_surface_m => tt_surface_h_m
1477	surf_m_liq => m_liq_h
1478	surf_m_liq_p => m_liq_h_p
1479	surf_tm_liq_m => tm_liq_h_m
1480	surf_m_soil => m_soil_h
1481	surf_t_soil => t_soil_h
1482
1483	k_off = -1
1484	j_off = 0
1485	i_off = 0
1486	ELSE
1487	surf => surf_lsm_v(l)
1488	surf_t_surface => t_surface_v(l)
1489	surf_t_surface_p => t_surface_v_p(l)
1490	surf_tt_surface_m => tt_surface_v_m(l)
1491	surf_m_liq => m_liq_v(l)
1492	surf_m_liq_p => m_liq_v_p(l)
1493	surf_tm_liq_m => tm_liq_v_m(l)
1494	surf_m_soil => m_soil_v(l)
1495	surf_t_soil => t_soil_v(l)
1496
1497	k_off = 0
1498	IF ( l == 0 ) THEN
1499	j_off = -1
1500	i_off = 0
1501	ELSEIF ( l == 1 ) THEN
1502	j_off = 1
1503	i_off = 0
1504	ELSEIF ( l == 2 ) THEN
1505	j_off = 0
1506	i_off = -1
1507	ELSEIF ( l == 3 ) THEN
1508	j_off = 0
1509	i_off = 1
1510	ENDIF
1511	ENDIF
1512
1513	!
1514	!-- Calculate the exner function for the current time step
1515	exn = ( surface_pressure / 1000.0_wp )**0.286_wp
1516
1517	DO m = 1, surf%ns
1518
1519	i = surf%i(m)
1520	j = surf%j(m)
1521	k = surf%k(m)
1522	!
1523	!-- Determine height index for radiation. Note, in clear-sky case radiation
1524	!-- arrays have rank 0 in first dimensions, so index must be zero. In case
1525	!-- of RRTMG radiation arrays have non-zero rank in first dimension, so that
1526	!-- radiation can be obtained at respective height level
1527	k_rad = MERGE( 0, k + k_off, radiation_scheme /= 'rrtmg' )
1528
1529	!
1530	!-- Define heat conductivity between surface and soil depending on surface
1531	!-- type. For vegetation, a skin layer parameterization is used. The new
1532	!-- parameterization uses a combination of two conductivities: a constant
1533	!-- conductivity for the skin layer, and a conductivity according to the
1534	!-- uppermost soil layer. For bare soil and pavements, no skin layer is
1535	!-- applied. In these cases, the temperature is assumed to be constant
1536	!-- between the surface and the first soil layer. The heat conductivity is
1537	!-- then derived from the soil/pavement properties.
1538	!-- For water surfaces, the conductivity is already set to 1E10.
1539	!-- Moreover, the heat capacity is set. For bare soil the heat capacity is
1540	!-- the capacity of the uppermost soil layer, for pavement it is that of
1541	!-- the material involved.
1542
1543	!
1544	!-- for vegetation type surfaces, the thermal conductivity of the soil is
1545	!-- needed
1546
1547	IF ( surf%vegetation_surface(m) ) THEN
1548
1549	lambda_h_sat = lambda_h_sm*(1.0_wp - surf%m_sat(nzb_soil,m)) &
1550	lambda_h_water ** surf_m_soil%var_2d(nzb_soil,m)
1551
1552	ke = 1.0_wp + LOG10( MAX( 0.1_wp, surf_m_soil%var_2d(nzb_soil,m) / &
1553	surf%m_sat(nzb_soil,m) ) )
1554
1555	lambda_soil = (ke * (lambda_h_sat - lambda_h_dry) + lambda_h_dry ) &
1556	* ddz_soil(nzb_soil) * 2.0_wp
1557
1558	!
1559	!-- When bare soil is set without a thermal conductivity (no skin layer),
1560	!-- a heat capacity is that of the soil layer, otherwise it is a
1561	!-- combination of the conductivities from the skin and the soil layer
1562	IF ( surf%lambda_surface_s(m) == 0.0_wp ) THEN
1563	surf%c_surface(m) = (rho_c_soil * (1.0_wp - surf%m_sat(nzb_soil,m))&
1564	+ rho_c_water * surf_m_soil%var_2d(nzb_soil,m) ) &
1565	* dz_soil(nzb_soil) * 0.5_wp
1566	lambda_surface = lambda_soil
1567
1568	ELSE IF ( surf_t_surface%var_1d(m) >= surf_t_soil%var_2d(nzb_soil,m))&
1569	THEN
1570	lambda_surface = surf%lambda_surface_s(m) * lambda_soil &
1571	/ ( surf%lambda_surface_s(m) + lambda_soil )
1572	ELSE
1573
1574	lambda_surface = surf%lambda_surface_u(m) * lambda_soil &
1575	/ ( surf%lambda_surface_u(m) + lambda_soil )
1576	ENDIF
1577	ELSE
1578	lambda_surface = surf%lambda_surface_s(m)
1579	ENDIF
1580
1581	!
1582	!-- Set heat capacity of the skin/surface. It is ususally zero when a skin
1583	!-- layer is used, and non-zero otherwise.
1584	c_surface_tmp = surf%c_surface(m)
1585
1586	!
1587	!-- First step: calculate aerodyamic resistance. As pt, us, ts
1588	!-- are not available for the prognostic time step, data from the last
1589	!-- time step is used here. Note that this formulation is the
1590	!-- equivalent to the ECMWF formulation using drag coefficients
1591	IF ( cloud_physics ) THEN
1592	pt1 = pt(k,j,i) + l_d_cp * pt_d_t(k) * ql(k,j,i)
1593	qv1 = q(k,j,i) - ql(k,j,i)
1594	ELSEIF ( cloud_droplets ) THEN
1595	pt1 = pt(k,j,i) + l_d_cp * pt_d_t(k) * ql(k,j,i)
1596	qv1 = q(k,j,i)
1597	ELSE
1598	pt1 = pt(k,j,i)
1599	IF ( humidity ) THEN
1600	qv1 = q(k,j,i)
1601	ELSE
1602	qv1 = 0.0_wp
1603	ENDIF
1604	ENDIF
1605	!
1606	!-- Calculate aerodynamical resistance. For horizontal and vertical
1607	!-- surfaces MOST is applied. Moreover, for vertical surfaces, resistance
1608	!-- can be obtain via parameterization of Mason (2000) /
1609	!-- Krayenhoff and Voogt (2006).
1610	!-- To do: detailed investigation which approach is better!
1611	IF ( horizontal .OR. .NOT. aero_resist_kray ) THEN
1612	surf%r_a(m) = ( pt1 - surf_lsm_h%pt_surface(m) ) / &
1613	( surf%ts(m) * surf%us(m) + 1.0E-20_wp )
1614	ELSE
1615	surf%r_a(m) = 1.0_wp / ( 11.8_wp + 4.2_wp * &
1616	SQRT( ( ( u(k,j,i) + u(k,j,i+1) ) * 0.5_wp )**2 + &
1617	( ( v(k,j,i) + v(k,j+1,i) ) * 0.5_wp )**2 + &
1618	( ( w(k,j,i) + w(k-1,j,i) ) * 0.5_wp )**2 ) &
1619	)
1620	ENDIF
1621	!
1622	!-- Make sure that the resistance does not drop to zero for neutral
1623	!-- stratification
1624	IF ( ABS( surf%r_a(m) ) < 1.0_wp ) surf%r_a(m) = 1.0_wp
1625	!
1626	!-- Second step: calculate canopy resistance r_canopy
1627	!-- f1-f3 here are defined as 1/f1-f3 as in ECMWF documentation
1628
1629	!-- f1: correction for incoming shortwave radiation (stomata close at
1630	!-- night)
1631	f1 = MIN( 1.0_wp, ( 0.004_wp * rad_sw_in(k_rad,j+j_off,i+i_off) &
1632	+ 0.05_wp ) / (0.81_wp * (0.004_wp &
1633	* rad_sw_in(k_rad,j+j_off,i+i_off) + 1.0_wp)) )
1634	!
1635	!-- f2: correction for soil moisture availability to plants (the
1636	!-- integrated soil moisture must thus be considered here)
1637	!-- f2 = 0 for very dry soils
1638	m_total = 0.0_wp
1639	DO ks = nzb_soil, nzt_soil
1640	m_total = m_total + surf%root_fr(ks,m) &
1641	* MAX( surf_m_soil%var_2d(ks,m), surf%m_wilt(ks,m) )
1642	ENDDO
1643
1644	!
1645	!-- The calculation of f2 is based on only one wilting point value for all
1646	!-- soil layers. The value at k=nzb_soil is used here as a proxy but might
1647	!-- need refinement in the future.
1648	IF ( m_total > surf%m_wilt(nzb_soil,m) .AND. &
1649	m_total < surf%m_fc(nzb_soil,m) ) THEN
1650	f2 = ( m_total - surf%m_wilt(nzb_soil,m) ) / &
1651	( surf%m_fc(nzb_soil,m) - surf%m_wilt(nzb_soil,m) )
1652	ELSEIF ( m_total >= surf%m_fc(nzb_soil,m) ) THEN
1653	f2 = 1.0_wp
1654	ELSE
1655	f2 = 1.0E-20_wp
1656	ENDIF
1657
1658	!
1659	!-- Calculate water vapour pressure at saturation and convert to hPa
1660	e_s = 0.01_wp * magnus( surf_t_surface%var_1d(m) )
1661
1662	!
1663	!-- f3: correction for vapour pressure deficit
1664	IF ( surf%g_d(m) /= 0.0_wp ) THEN
1665	!
1666	!-- Calculate vapour pressure
1667	e = qv1 * surface_pressure / ( qv1 + 0.622_wp )
1668	f3 = EXP ( - surf%g_d(m) * (e_s - e) )
1669	ELSE
1670	f3 = 1.0_wp
1671	ENDIF
1672	!
1673	!-- Calculate canopy resistance. In case that c_veg is 0 (bare soils),
1674	!-- this calculation is obsolete, as r_canopy is not used below.
1675	!-- To do: check for very dry soil -> r_canopy goes to infinity
1676	surf%r_canopy(m) = surf%r_canopy_min(m) / &
1677	( surf%lai(m) * f1 * f2 * f3 + 1.0E-20_wp )
1678	!
1679	!-- Third step: calculate bare soil resistance r_soil.
1680	m_min = surf%c_veg(m) * surf%m_wilt(nzb_soil,m) + &
1681	( 1.0_wp - surf%c_veg(m) ) * surf%m_res(nzb_soil,m)
1682
1683
1684	f2 = ( surf_m_soil%var_2d(nzb_soil,m) - m_min ) &
1685	/ ( surf%m_fc(nzb_soil,m) - m_min )
1686	f2 = MAX( f2, 1.0E-20_wp )
1687	f2 = MIN( f2, 1.0_wp )
1688
1689	surf%r_soil(m) = surf%r_soil_min(m) / f2
1690
1691	!
1692	!-- Calculate the maximum possible liquid water amount on plants and
1693	!-- bare surface. For vegetated surfaces, a maximum depth of 0.2 mm is
1694	!-- assumed, while paved surfaces might hold up 1 mm of water. The
1695	!-- liquid water fraction for paved surfaces is calculated after
1696	!-- Noilhan & Planton (1989), while the ECMWF formulation is used for
1697	!-- vegetated surfaces and bare soils.
1698	IF ( surf%pavement_surface(m) ) THEN
1699	m_liq_max = m_max_depth * 5.0_wp
1700	surf%c_liq(m) = MIN( 1.0_wp, (surf_m_liq%var_1d(m) / m_liq_max)**0.67 )
1701	ELSE
1702	m_liq_max = m_max_depth * ( surf%c_veg(m) * surf%lai(m) &
1703	+ ( 1.0_wp - surf%c_veg(m) ) )
1704	surf%c_liq(m) = MIN( 1.0_wp, surf_m_liq%var_1d(m) / m_liq_max )
1705	ENDIF
1706	!
1707	!-- Calculate saturation specific humidity
1708	q_s = 0.622_wp * e_s / ( surface_pressure - e_s )
1709	!
1710	!-- In case of dewfall, set evapotranspiration to zero
1711	!-- All super-saturated water is then removed from the air
1712	IF ( humidity .AND. q_s <= qv1 ) THEN
1713	surf%r_canopy(m) = 0.0_wp
1714	surf%r_soil(m) = 0.0_wp
1715	ENDIF
1716
1717	!
1718	!-- Calculate coefficients for the total evapotranspiration
1719	!-- In case of water surface, set vegetation and soil fluxes to zero.
1720	!-- For pavements, only evaporation of liquid water is possible.
1721	IF ( surf%water_surface(m) ) THEN
1722	f_qsws_veg = 0.0_wp
1723	f_qsws_soil = 0.0_wp
1724	f_qsws_liq = rho_lv / surf%r_a(m)
1725	ELSEIF ( surf%pavement_surface (m) ) THEN
1726	f_qsws_veg = 0.0_wp
1727	f_qsws_soil = 0.0_wp
1728	f_qsws_liq = rho_lv * surf%c_liq(m) / surf%r_a(m)
1729	ELSE
1730	f_qsws_veg = rho_lv * surf%c_veg(m) * &
1731	( 1.0_wp - surf%c_liq(m) ) / &
1732	( surf%r_a(m) + surf%r_canopy(m) )
1733	f_qsws_soil = rho_lv * (1.0_wp - surf%c_veg(m) ) / &
1734	( surf%r_a(m) + surf%r_soil(m) )
1735	f_qsws_liq = rho_lv * surf%c_veg(m) * surf%c_liq(m) / &
1736	surf%r_a(m)
1737	ENDIF
1738
1739	f_shf = rho_cp / surf%r_a(m)
1740	f_qsws = f_qsws_veg + f_qsws_soil + f_qsws_liq
1741	!
1742	!-- Calculate derivative of q_s for Taylor series expansion
1743	e_s_dt = e_s * ( 17.62_wp / ( surf_t_surface%var_1d(m) - 29.65_wp) - &
1744	17.62_wp*( surf_t_surface%var_1d(m) - 273.15_wp) &
1745	/ ( surf_t_surface%var_1d(m) - 29.65_wp)**2 )
1746
1747	dq_s_dt = 0.622_wp * e_s_dt / ( surface_pressure - e_s_dt )
1748	!
1749	!-- Add LW up so that it can be removed in prognostic equation
1750	surf%rad_net_l(m) = rad_net(j,i) + rad_lw_out(nzb,j,i)
1751	!
1752	!-- Calculate new skin temperature
1753	IF ( humidity ) THEN
1754	!
1755	!-- Numerator of the prognostic equation
1756	coef_1 = surf%rad_net_l(m) + rad_lw_out_change_0(j,i) &
1757	* surf_t_surface%var_1d(m) - rad_lw_out(nzb,j,i) &
1758	+ f_shf * pt1 + f_qsws * ( qv1 - q_s &
1759	+ dq_s_dt * surf_t_surface%var_1d(m) ) + lambda_surface &
1760	* surf_t_soil%var_2d(nzb_soil,m)
1761	!
1762	!-- Denominator of the prognostic equation
1763	coef_2 = rad_lw_out_change_0(j,i) + f_qsws * dq_s_dt &
1764	+ lambda_surface + f_shf / exn
1765	ELSE
1766	!
1767	!-- Numerator of the prognostic equation
1768	coef_1 = surf%rad_net_l(m) + rad_lw_out_change_0(j,i) &
1769	* surf_t_surface%var_1d(m) - rad_lw_out(nzb,j,i) &
1770	+ f_shf * pt1 + lambda_surface &
1771	* surf_t_soil%var_2d(nzb_soil,m)
1772	!
1773	!-- Denominator of the prognostic equation
1774	coef_2 = rad_lw_out_change_0(j,i) + lambda_surface + f_shf / exn
1775
1776	ENDIF
1777
1778	tend = 0.0_wp
1779
1780	!
1781	!-- Implicit solution when the surface layer has no heat capacity,
1782	!-- otherwise use RK3 scheme.
1783	surf_t_surface_p%var_1d(m) = ( coef_1 * dt_3d * tsc(2) + c_surface_tmp *&
1784	surf_t_surface%var_1d(m) ) / ( c_surface_tmp + coef_2&
1785	* dt_3d * tsc(2) )
1786
1787	!
1788	!-- Add RK3 term
1789	IF ( c_surface_tmp /= 0.0_wp ) THEN
1790
1791	surf_t_surface_p%var_1d(m) = surf_t_surface_p%var_1d(m) + dt_3d * &
1792	tsc(3) * surf_tt_surface_m%var_1d(m)
1793
1794	!
1795	!-- Calculate true tendency
1796	tend = ( surf_t_surface_p%var_1d(m) - surf_t_surface%var_1d(m) - &
1797	dt_3d * tsc(3) * surf_tt_surface_m%var_1d(m)) / (dt_3d * tsc(2))
1798	!
1799	!-- Calculate t_surface tendencies for the next Runge-Kutta step
1800	IF ( timestep_scheme(1:5) == 'runge' ) THEN
1801	IF ( intermediate_timestep_count == 1 ) THEN
1802	surf_tt_surface_m%var_1d(m) = tend
1803	ELSEIF ( intermediate_timestep_count < &
1804	intermediate_timestep_count_max ) THEN
1805	surf_tt_surface_m%var_1d(m) = -9.5625_wp * tend + &
1806	5.3125_wp * surf_tt_surface_m%var_1d(m)
1807	ENDIF
1808	ENDIF
1809	ENDIF
1810
1811	!
1812	!-- In case of fast changes in the skin temperature, it is possible to
1813	!-- update the radiative fluxes independently from the prescribed
1814	!-- radiation call frequency. This effectively prevents oscillations,
1815	!-- especially when setting skip_time_do_radiation /= 0. The threshold
1816	!-- value of 0.2 used here is just a first guess. This method should be
1817	!-- revised in the future as tests have shown that the threshold is
1818	!-- often reached, when no oscillations would occur (causes immense
1819	!-- computing time for the radiation code).
1820	IF ( ABS( surf_t_surface_p%var_1d(m) - surf_t_surface%var_1d(m) ) &
1821	> 0.2_wp .AND. unscheduled_radiation_calls ) THEN
1822	force_radiation_call_l = .TRUE.
1823	ENDIF
1824
1825	pt(k+k_off,j+j_off,i+i_off) = surf_t_surface_p%var_1d(m) / exn !is actually no air temperature
1826	surf%pt_surface(m) = surf_t_surface_p%var_1d(m) / exn
1827
1828	!
1829	!-- Calculate fluxes
1830	surf%rad_net_l(m) = surf%rad_net_l(m) + &
1831	rad_lw_out_change_0(j,i) &
1832	* surf_t_surface%var_1d(m) - rad_lw_out(nzb,j,i) &
1833	- rad_lw_out_change_0(j,i) * surf_t_surface_p%var_1d(m)
1834
1835	rad_net(j,i) = surf%rad_net_l(m)
1836	rad_lw_out(nzb,j,i) = rad_lw_out(nzb,j,i) + rad_lw_out_change_0(j,i) * &
1837	( surf_t_surface_p%var_1d(m) - surf_t_surface%var_1d(m) )
1838
1839	surf%ghf(m) = lambda_surface * ( surf_t_surface_p%var_1d(m) &
1840	- surf_t_soil%var_2d(nzb_soil,m) )
1841
1842	surf%shf(m) = - f_shf * ( pt1 - surf%pt_surface(m) ) / cp
1843
1844
1845	IF ( humidity ) THEN
1846	surf%qsws(m) = - f_qsws * ( qv1 - q_s + dq_s_dt &
1847	* surf_t_surface%var_1d(m) - dq_s_dt * &
1848	surf_t_surface_p%var_1d(m) )
1849
1850	surf%qsws_veg(m) = - f_qsws_veg * ( qv1 - q_s &
1851	+ dq_s_dt * surf_t_surface%var_1d(m) - dq_s_dt &
1852	* surf_t_surface_p%var_1d(m) )
1853
1854	surf%qsws_soil(m) = - f_qsws_soil * ( qv1 - q_s &
1855	+ dq_s_dt * surf_t_surface%var_1d(m) - dq_s_dt &
1856	* surf_t_surface_p%var_1d(m) )
1857
1858	surf%qsws_liq(m) = - f_qsws_liq * ( qv1 - q_s &
1859	+ dq_s_dt * surf_t_surface%var_1d(m) - dq_s_dt &
1860	* surf_t_surface_p%var_1d(m) )
1861	ENDIF
1862
1863	!
1864	!-- Calculate the true surface resistance
1865	IF ( .NOT. humidity ) THEN
1866	surf%r_s(m) = 1.0E10_wp
1867	ELSE
1868	surf%r_s(m) = - rho_lv * ( qv1 - q_s + dq_s_dt &
1869	* surf_t_surface%var_1d(m) - dq_s_dt * &
1870	surf_t_surface_p%var_1d(m) ) / &
1871	(surf%qsws(m) + 1.0E-20) - surf%r_a(m)
1872	ENDIF
1873
1874	!
1875	!-- Calculate change in liquid water reservoir due to dew fall or
1876	!-- evaporation of liquid water
1877	IF ( humidity ) THEN
1878	!
1879	!-- If precipitation is activated, add rain water to qsws_liq
1880	!-- and qsws_soil according the the vegetation coverage.
1881	!-- precipitation_rate is given in mm.
1882	IF ( precipitation ) THEN
1883
1884	!
1885	!-- Add precipitation to liquid water reservoir, if possible.
1886	!-- Otherwise, add the water to soil. In case of
1887	!-- pavements, the exceeding water amount is implicitely removed
1888	!-- as runoff as qsws_soil is then not used in the soil model
1889	IF ( surf_m_liq%var_1d(m) /= m_liq_max ) THEN
1890	surf%qsws_liq(m) = surf%qsws_liq(m) &
1891	+ surf%c_veg(m) * prr(k+k_off,j+j_off,i+i_off)&
1892	* hyrho(k+k_off) &
1893	* 0.001_wp * rho_l * l_v
1894	ELSE
1895	surf%qsws_soil(m) = surf%qsws_soil(m) &
1896	+ surf%c_veg(m) * prr(k+k_off,j+j_off,i+i_off)&
1897	* hyrho(k+k_off) &
1898	* 0.001_wp * rho_l * l_v
1899	ENDIF
1900
1901	!-- Add precipitation to bare soil according to the bare soil
1902	!-- coverage.
1903	surf%qsws_soil(m) = surf%qsws_soil(m) + ( 1.0_wp &
1904	- surf%c_veg(m) ) * prr(k+k_off,j+j_off,i+i_off)&
1905	* hyrho(k+k_off) &
1906	* 0.001_wp * rho_l * l_v
1907	ENDIF
1908
1909	!
1910	!-- If the air is saturated, check the reservoir water level
1911	IF ( surf%qsws(m) < 0.0_wp ) THEN
1912	!
1913	!-- Check if reservoir is full (avoid values > m_liq_max)
1914	!-- In that case, qsws_liq goes to qsws_soil. In this
1915	!-- case qsws_veg is zero anyway (because c_liq = 1),
1916	!-- so that tend is zero and no further check is needed
1917	IF ( surf_m_liq%var_1d(m) == m_liq_max ) THEN
1918	surf%qsws_soil(m) = surf%qsws_soil(m) + surf%qsws_liq(m)
1919
1920	surf%qsws_liq(m) = 0.0_wp
1921	ENDIF
1922
1923	!
1924	!-- In case qsws_veg becomes negative (unphysical behavior),
1925	!-- let the water enter the liquid water reservoir as dew on the
1926	!-- plant
1927	IF ( surf%qsws_veg(m) < 0.0_wp ) THEN
1928	surf%qsws_liq(m) = surf%qsws_liq(m) + surf%qsws_veg(m)
1929	surf%qsws_veg(m) = 0.0_wp
1930	ENDIF
1931	ENDIF
1932
1933	surf%qsws(m) = surf%qsws(m) / l_v
1934
1935	tend = - surf%qsws_liq(m) * drho_l_lv
1936	surf_m_liq_p%var_1d(m) = surf_m_liq%var_1d(m) + dt_3d * &
1937	( tsc(2) * tend + &
1938	tsc(3) * surf_tm_liq_m%var_1d(m) )
1939	!
1940	!-- Check if reservoir is overfull -> reduce to maximum
1941	!-- (conservation of water is violated here)
1942	surf_m_liq_p%var_1d(m) = MIN( surf_m_liq_p%var_1d(m),m_liq_max )
1943
1944	!
1945	!-- Check if reservoir is empty (avoid values < 0.0)
1946	!-- (conservation of water is violated here)
1947	surf_m_liq_p%var_1d(m) = MAX( surf_m_liq_p%var_1d(m), 0.0_wp )
1948	!
1949	!-- Calculate m_liq tendencies for the next Runge-Kutta step
1950	IF ( timestep_scheme(1:5) == 'runge' ) THEN
1951	IF ( intermediate_timestep_count == 1 ) THEN
1952	surf_tm_liq_m%var_1d(m) = tend
1953	ELSEIF ( intermediate_timestep_count < &
1954	intermediate_timestep_count_max ) THEN
1955	surf_tm_liq_m%var_1d(m) = -9.5625_wp * tend + &
1956	5.3125_wp * surf_tm_liq_m%var_1d(m)
1957	ENDIF
1958	ENDIF
1959
1960	ENDIF
1961
1962	ENDDO
1963
1964	!
1965	!-- Make a logical OR for all processes. Force radiation call if at
1966	!-- least one processor reached the threshold change in skin temperature
1967	IF ( unscheduled_radiation_calls .AND. intermediate_timestep_count &
1968	== intermediate_timestep_count_max-1 ) THEN
1969	#if defined( __parallel )
1970	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
1971	CALL MPI_ALLREDUCE( force_radiation_call_l, force_radiation_call, &
1972	1, MPI_LOGICAL, MPI_LOR, comm2d, ierr )
1973	#else
1974	force_radiation_call = force_radiation_call_l
1975	#endif
1976	force_radiation_call_l = .FALSE.
1977	ENDIF
1978
1979	!
1980	!-- Calculate surface specific humidity
1981	IF ( humidity ) THEN
1982	CALL calc_q_surface
1983	ENDIF
1984
1985	!
1986	!-- Calculate new roughness lengths (for water surfaces only)
1987	IF ( horizontal .AND. .NOT. constant_roughness ) CALL calc_z0_water_surface
1988
1989	CONTAINS
1990	!------------------------------------------------------------------------------!
1991	! Description:
1992	! ------------
1993	!> Calculation of specific humidity of the skin layer (surface). It is assumend
1994	!> that the skin is always saturated.
1995	!------------------------------------------------------------------------------!
1996	SUBROUTINE calc_q_surface
1997
1998	USE diagnostic_quantities_mod
1999
2000	IMPLICIT NONE
2001
2002	REAL(wp) :: resistance !< aerodynamic and soil resistance term
2003
2004	DO m = 1, surf%ns
2005
2006	i = surf%i(m)
2007	j = surf%j(m)
2008	k = surf%k(m)
2009	!
2010	!-- Calculate water vapour pressure at saturation and convert to hPa
2011	e_s = 0.01_wp * magnus( surf_t_surface_p%var_1d(m) )
2012
2013	!
2014	!-- Calculate specific humidity at saturation
2015	q_s = 0.622_wp * e_s / ( surface_pressure - e_s )
2016
2017	resistance = surf%r_a(m) / ( surf%r_a(m) + surf%r_s(m) )
2018
2019	!
2020	!-- Calculate specific humidity at surface
2021	IF ( cloud_physics ) THEN
2022	q(k+k_off,j+j_off,i+i_off) = resistance * q_s + &
2023	( 1.0_wp - resistance ) * &
2024	( q(k,j,i) - ql(k,j,i) )
2025	ELSE
2026	q(k+k_off,j+j_off,i+i_off) = resistance * q_s + &
2027	( 1.0_wp - resistance ) * &
2028	q(k,j,i)
2029	ENDIF
2030
2031	!
2032	!-- Update virtual potential temperature
2033	vpt(k+k_off,j+j_off,i+i_off) = pt(k+k_off,j+j_off,i+i_off) * &
2034	( 1.0_wp + 0.61_wp * q(k+k_off,j+j_off,i+i_off) )
2035
2036	ENDDO
2037
2038	END SUBROUTINE calc_q_surface
2039
2040
2041
2042	END SUBROUTINE lsm_energy_balance
2043
2044
2045	!------------------------------------------------------------------------------!
2046	! Description:
2047	! ------------
2048	!> Header output for land surface model
2049	!------------------------------------------------------------------------------!
2050	SUBROUTINE lsm_header ( io )
2051
2052
2053	IMPLICIT NONE
2054
2055	CHARACTER (LEN=86) :: t_soil_chr !< String for soil temperature profile
2056	CHARACTER (LEN=86) :: roots_chr !< String for root profile
2057	CHARACTER (LEN=86) :: vertical_index_chr !< String for the vertical index
2058	CHARACTER (LEN=86) :: m_soil_chr !< String for soil moisture
2059	CHARACTER (LEN=86) :: soil_depth_chr !< String for soil depth
2060	CHARACTER (LEN=10) :: coor_chr !< Temporary string
2061
2062	INTEGER(iwp) :: i !< Loop index over soil layers
2063
2064	INTEGER(iwp), INTENT(IN) :: io !< Unit of the output file
2065
2066	t_soil_chr = ''
2067	m_soil_chr = ''
2068	soil_depth_chr = ''
2069	roots_chr = ''
2070	vertical_index_chr = ''
2071
2072	i = 1
2073	DO i = nzb_soil, nzt_soil
2074	WRITE (coor_chr,'(F10.2,7X)') soil_temperature(i)
2075	t_soil_chr = TRIM( t_soil_chr ) // ' ' // TRIM( coor_chr )
2076
2077	WRITE (coor_chr,'(F10.2,7X)') soil_moisture(i)
2078	m_soil_chr = TRIM( m_soil_chr ) // ' ' // TRIM( coor_chr )
2079
2080	WRITE (coor_chr,'(F10.2,7X)') - zs(i)
2081	soil_depth_chr = TRIM( soil_depth_chr ) // ' ' // TRIM( coor_chr )
2082
2083	WRITE (coor_chr,'(F10.2,7X)') root_fraction(i)
2084	roots_chr = TRIM( roots_chr ) // ' ' // TRIM( coor_chr )
2085
2086	WRITE (coor_chr,'(I10,7X)') i
2087	vertical_index_chr = TRIM( vertical_index_chr ) // ' ' // &
2088	TRIM( coor_chr )
2089	ENDDO
2090
2091	!
2092	!-- Write land surface model header
2093	WRITE( io, 1 )
2094	IF ( conserve_water_content ) THEN
2095	WRITE( io, 2 )
2096	ELSE
2097	WRITE( io, 3 )
2098	ENDIF
2099
2100	WRITE( io, 4 ) TRIM( vegetation_type_name(vegetation_type) ), &
2101	TRIM (soil_type_name(soil_type) )
2102	WRITE( io, 5 ) TRIM( soil_depth_chr ), TRIM( t_soil_chr ), &
2103	TRIM( m_soil_chr ), TRIM( roots_chr ), &
2104	TRIM( vertical_index_chr )
2105
2106	1 FORMAT (//' Land surface model information:'/ &
2107	' ------------------------------'/)
2108	2 FORMAT (' --> Soil bottom is closed (water content is conserved', &
2109	', default)')
2110	3 FORMAT (' --> Soil bottom is open (water content is not conserved)')
2111	4 FORMAT (' --> Land surface type : ',A,/ &
2112	' --> Soil porosity type : ',A)
2113	5 FORMAT (/' Initial soil temperature and moisture profile:'// &
2114	' Height: ',A,' m'/ &
2115	' Temperature: ',A,' K'/ &
2116	' Moisture: ',A,' m3/m3'/ &
2117	' Root fraction: ',A,' '/ &
2118	' Grid point: ',A)
2119
2120	END SUBROUTINE lsm_header
2121
2122
2123	!------------------------------------------------------------------------------!
2124	! Description:
2125	! ------------
2126	!> Initialization of the land surface model
2127	!------------------------------------------------------------------------------!
2128	SUBROUTINE lsm_init
2129
2130	USE control_parameters, &
2131	ONLY: message_string
2132
2133	IMPLICIT NONE
2134
2135	INTEGER(iwp) :: i !< running index
2136	INTEGER(iwp) :: i_off !< index offset of surface element, seen from atmospheric grid point
2137	INTEGER(iwp) :: j !< running index
2138	INTEGER(iwp) :: j_off !< index offset of surface element, seen from atmospheric grid point
2139	INTEGER(iwp) :: k !< running index
2140	INTEGER(iwp) :: kn !< running index
2141	INTEGER(iwp) :: ko !< running index
2142	INTEGER(iwp) :: kroot !< running index
2143	INTEGER(iwp) :: kzs !< running index
2144	INTEGER(iwp) :: l !< running index surface facing
2145	INTEGER(iwp) :: m !< running index
2146	INTEGER(iwp) :: n_soil_layers_total !< temperature variable, stores the total number of soil layers + 4
2147
2148
2149	REAL(wp) :: pt1 !< potential temperature at first grid level
2150
2151	REAL(wp), DIMENSION(:), ALLOCATABLE :: bound, bound_root_fr !< temporary arrays for storing index bounds
2152
2153	!
2154	!-- Calculate Exner function
2155	exn = ( surface_pressure / 1000.0_wp )**0.286_wp
2156	!
2157	!-- If no cloud physics is used, rho_surface has not been calculated before
2158	IF ( .NOT. cloud_physics ) THEN
2159	rho_surface = surface_pressure * 100.0_wp / ( r_d * pt_surface * exn )
2160	ENDIF
2161
2162	!
2163	!-- Calculate frequently used parameters
2164	rho_cp = cp * rho_surface
2165	rd_d_rv = r_d / r_v
2166	rho_lv = rho_surface * l_v
2167	drho_l_lv = 1.0_wp / (rho_l * l_v)
2168
2169	!
2170	!-- Set initial values for prognostic quantities
2171	!-- Horizontal surfaces
2172	tt_surface_h_m%var_1d = 0.0_wp
2173	tt_soil_h_m%var_2d = 0.0_wp
2174	tm_soil_h_m%var_2d = 0.0_wp
2175	tm_liq_h_m%var_1d = 0.0_wp
2176	surf_lsm_h%c_liq = 0.0_wp
2177
2178	surf_lsm_h%ghf = 0.0_wp
2179
2180	surf_lsm_h%qsws_liq = 0.0_wp
2181	surf_lsm_h%qsws_soil = 0.0_wp
2182	surf_lsm_h%qsws_veg = 0.0_wp
2183
2184	surf_lsm_h%r_a = 50.0_wp
2185	surf_lsm_h%r_s = 50.0_wp
2186	surf_lsm_h%r_canopy = 0.0_wp
2187	surf_lsm_h%r_soil = 0.0_wp
2188	!
2189	!-- Do the same for vertical surfaces
2190	DO l = 0, 3
2191	tt_surface_v_m(l)%var_1d = 0.0_wp
2192	tt_soil_v_m(l)%var_2d = 0.0_wp
2193	tm_soil_v_m(l)%var_2d = 0.0_wp
2194	tm_liq_v_m(l)%var_1d = 0.0_wp
2195	surf_lsm_v(l)%c_liq = 0.0_wp
2196
2197	surf_lsm_v(l)%ghf = 0.0_wp
2198
2199	surf_lsm_v(l)%qsws_liq = 0.0_wp
2200	surf_lsm_v(l)%qsws_soil = 0.0_wp
2201	surf_lsm_v(l)%qsws_veg = 0.0_wp
2202
2203	surf_lsm_v(l)%r_a = 50.0_wp
2204	surf_lsm_v(l)%r_s = 50.0_wp
2205	surf_lsm_v(l)%r_canopy = 0.0_wp
2206	surf_lsm_v(l)%r_soil = 0.0_wp
2207	ENDDO
2208
2209	!
2210	!-- Allocate 3D soil model arrays
2211	!-- First, for horizontal surfaces
2212	ALLOCATE ( surf_lsm_h%alpha_vg(nzb_soil:nzt_soil,1:surf_lsm_h%ns) )
2213	ALLOCATE ( surf_lsm_h%gamma_w_sat(nzb_soil:nzt_soil,1:surf_lsm_h%ns) )
2214	ALLOCATE ( surf_lsm_h%lambda_h(nzb_soil:nzt_soil,1:surf_lsm_h%ns) )
2215	ALLOCATE ( surf_lsm_h%lambda_h_def(nzb_soil:nzt_soil,1:surf_lsm_h%ns))
2216	ALLOCATE ( surf_lsm_h%l_vg(nzb_soil:nzt_soil,1:surf_lsm_h%ns) )
2217	ALLOCATE ( surf_lsm_h%m_fc(nzb_soil:nzt_soil,1:surf_lsm_h%ns) )
2218	ALLOCATE ( surf_lsm_h%m_res(nzb_soil:nzt_soil,1:surf_lsm_h%ns) )
2219	ALLOCATE ( surf_lsm_h%m_sat(nzb_soil:nzt_soil,1:surf_lsm_h%ns) )
2220	ALLOCATE ( surf_lsm_h%m_wilt(nzb_soil:nzt_soil,1:surf_lsm_h%ns) )
2221	ALLOCATE ( surf_lsm_h%n_vg(nzb_soil:nzt_soil,1:surf_lsm_h%ns) )
2222	ALLOCATE ( surf_lsm_h%rho_c_total(nzb_soil:nzt_soil,1:surf_lsm_h%ns) )
2223	ALLOCATE ( surf_lsm_h%rho_c_total_def(nzb_soil:nzt_soil,1:surf_lsm_h%ns) )
2224	ALLOCATE ( surf_lsm_h%root_fr(nzb_soil:nzt_soil,1:surf_lsm_h%ns) )
2225
2226	surf_lsm_h%lambda_h = 0.0_wp
2227	!
2228	!-- If required, allocate humidity-related variables for the soil model
2229	IF ( humidity ) THEN
2230	ALLOCATE ( surf_lsm_h%lambda_w(nzb_soil:nzt_soil,1:surf_lsm_h%ns) )
2231	ALLOCATE ( surf_lsm_h%gamma_w(nzb_soil:nzt_soil,1:surf_lsm_h%ns) )
2232
2233	surf_lsm_h%lambda_w = 0.0_wp
2234	ENDIF
2235	!
2236	!-- For vertical surfaces
2237	DO l = 0, 3
2238	ALLOCATE ( surf_lsm_v(l)%alpha_vg(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns) )
2239	ALLOCATE ( surf_lsm_v(l)%gamma_w_sat(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns) )
2240	ALLOCATE ( surf_lsm_v(l)%lambda_h(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns) )
2241	ALLOCATE ( surf_lsm_v(l)%lambda_h_def(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns))
2242	ALLOCATE ( surf_lsm_v(l)%l_vg(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns) )
2243	ALLOCATE ( surf_lsm_v(l)%m_fc(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns) )
2244	ALLOCATE ( surf_lsm_v(l)%m_res(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns) )
2245	ALLOCATE ( surf_lsm_v(l)%m_sat(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns) )
2246	ALLOCATE ( surf_lsm_v(l)%m_wilt(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns) )
2247	ALLOCATE ( surf_lsm_v(l)%n_vg(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns) )
2248	ALLOCATE ( surf_lsm_v(l)%rho_c_total(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns) )
2249	ALLOCATE ( surf_lsm_v(l)%rho_c_total_def(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns) )
2250	ALLOCATE ( surf_lsm_v(l)%root_fr(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns) )
2251
2252	surf_lsm_v(l)%lambda_h = 0.0_wp
2253
2254	!
2255	!-- If required, allocate humidity-related variables for the soil model
2256	IF ( humidity ) THEN
2257	ALLOCATE ( surf_lsm_v(l)%lambda_w(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns) )
2258	ALLOCATE ( surf_lsm_v(l)%gamma_w(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns) )
2259
2260	surf_lsm_v(l)%lambda_w = 0.0_wp
2261	ENDIF
2262	ENDDO
2263
2264	!
2265	!-- Set flag parameter for the prescribed surface type depending on user
2266	!-- input.
2267	DO m = 1, surf_lsm_h%ns
2268
2269	SELECT CASE ( TRIM( surface_type ) )
2270
2271	CASE ( 'vegetation' )
2272
2273	surf_lsm_h%vegetation_surface(m) = .TRUE.
2274
2275	CASE ( 'water' )
2276
2277	surf_lsm_h%water_surface(m) = .TRUE.
2278
2279	CASE ( 'pavement' )
2280
2281	surf_lsm_h%pavement_surface(m) = .TRUE.
2282
2283	END SELECT
2284
2285	ENDDO
2286
2287	DO l = 0, 3
2288	SELECT CASE ( TRIM( surface_type ) )
2289
2290	CASE ( 'vegetation' )
2291
2292	surf_lsm_v(l)%vegetation_surface = .TRUE.
2293
2294	CASE ( 'water' )
2295
2296	surf_lsm_v(l)%water_surface = .TRUE.
2297
2298	CASE ( 'pavement' )
2299
2300	surf_lsm_h%pavement_surface = .TRUE.
2301
2302	END SELECT
2303
2304	ENDDO
2305	!
2306	!-- Initialize standard soil types. It is possible to overwrite each
2307	!-- parameter by setting the respecticy NAMELIST variable to a
2308	!-- value /= 9999999.9.
2309	IF ( soil_type /= 0 ) THEN
2310
2311	IF ( alpha_vangenuchten == 9999999.9_wp ) THEN
2312	alpha_vangenuchten = soil_pars(0,soil_type)
2313	ENDIF
2314
2315	IF ( l_vangenuchten == 9999999.9_wp ) THEN
2316	l_vangenuchten = soil_pars(1,soil_type)
2317	ENDIF
2318
2319	IF ( n_vangenuchten == 9999999.9_wp ) THEN
2320	n_vangenuchten = soil_pars(2,soil_type)
2321	ENDIF
2322
2323	IF ( hydraulic_conductivity == 9999999.9_wp ) THEN
2324	hydraulic_conductivity = soil_pars(3,soil_type)
2325	ENDIF
2326
2327	IF ( saturation_moisture == 9999999.9_wp ) THEN
2328	saturation_moisture = soil_pars(4,soil_type)
2329	ENDIF
2330
2331	IF ( field_capacity == 9999999.9_wp ) THEN
2332	field_capacity = soil_pars(5,soil_type)
2333	ENDIF
2334
2335	IF ( wilting_point == 9999999.9_wp ) THEN
2336	wilting_point = soil_pars(6,soil_type)
2337	ENDIF
2338
2339	IF ( residual_moisture == 9999999.9_wp ) THEN
2340	residual_moisture = soil_pars(7,soil_type)
2341	ENDIF
2342
2343	ENDIF
2344
2345	!
2346	!-- Check whether parameters from the lookup tables are to be used
2347	IF ( vegetation_type /= 0 ) THEN
2348
2349	IF ( min_canopy_resistance == 9999999.9_wp ) THEN
2350	min_canopy_resistance = vegetation_pars(0,vegetation_type)
2351	ENDIF
2352
2353	IF ( leaf_area_index == 9999999.9_wp ) THEN
2354	leaf_area_index = vegetation_pars(1,vegetation_type)
2355	ENDIF
2356
2357	IF ( vegetation_coverage == 9999999.9_wp ) THEN
2358	vegetation_coverage = vegetation_pars(2,vegetation_type)
2359	ENDIF
2360
2361	IF ( canopy_resistance_coefficient == 9999999.9_wp ) THEN
2362	canopy_resistance_coefficient= vegetation_pars(3,vegetation_type)
2363	ENDIF
2364
2365	IF ( z0_vegetation == 9999999.9_wp ) THEN
2366	z0_vegetation = vegetation_pars(4,vegetation_type)
2367	ENDIF
2368
2369	IF ( z0h_vegetation == 9999999.9_wp ) THEN
2370	z0h_vegetation = vegetation_pars(5,vegetation_type)
2371	ENDIF
2372
2373	IF ( lambda_surface_stable == 9999999.9_wp ) THEN
2374	lambda_surface_stable = vegetation_pars(6,vegetation_type)
2375	ENDIF
2376
2377	IF ( lambda_surface_unstable == 9999999.9_wp ) THEN
2378	lambda_surface_unstable = vegetation_pars(7,vegetation_type)
2379	ENDIF
2380
2381	IF ( f_shortwave_incoming == 9999999.9_wp ) THEN
2382	f_shortwave_incoming = vegetation_pars(8,vegetation_type)
2383	ENDIF
2384
2385	IF ( c_surface == 9999999.9_wp ) THEN
2386	c_surface = vegetation_pars(9,vegetation_type)
2387	ENDIF
2388
2389	IF ( albedo_type == 9999999 .AND. albedo == 9999999.9_wp ) THEN
2390	albedo_type = INT(vegetation_pars(10,vegetation_type))
2391	ENDIF
2392
2393	IF ( emissivity == 9999999.9_wp ) THEN
2394	emissivity = vegetation_pars(11,vegetation_type)
2395	ENDIF
2396
2397	ENDIF
2398
2399	IF ( water_type /= 0 ) THEN
2400
2401	IF ( water_temperature == 9999999.9_wp ) THEN
2402	water_temperature = water_pars(0,water_type)
2403	ENDIF
2404
2405	IF ( z0_water == 9999999.9_wp ) THEN
2406	z0_water = water_pars(1,water_type)
2407	ENDIF
2408
2409	IF ( z0h_water == 9999999.9_wp ) THEN
2410	z0h_water = water_pars(2,water_type)
2411	ENDIF
2412
2413	IF ( albedo_type == 9999999 .AND. albedo == 9999999.9_wp ) THEN
2414	albedo_type = INT(water_pars(3,water_type))
2415	ENDIF
2416
2417	IF ( emissivity == 9999999.9_wp ) THEN
2418	emissivity = water_pars(4,water_type)
2419	ENDIF
2420
2421	ENDIF
2422
2423	IF ( pavement_type /= 0 ) THEN
2424
2425	!
2426	!-- When a pavement_type is used, overwrite a possible setting of
2427	!-- the pavement depth as it is already defined by the pavement type
2428	pavement_depth_level = 0
2429
2430	IF ( z0_pavement == 9999999.9_wp ) THEN
2431	z0_pavement = pavement_pars(0,pavement_type)
2432	ENDIF
2433
2434	IF ( z0h_pavement == 9999999.9_wp ) THEN
2435	z0h_pavement = pavement_pars(1,pavement_type)
2436	ENDIF
2437
2438	IF ( pavement_heat_conduct == 9999999.9_wp ) THEN
2439	pavement_heat_conduct = pavement_subsurface_pars_1(0,pavement_type)
2440	ENDIF
2441
2442	IF ( pavement_heat_capacity == 9999999.9_wp ) THEN
2443	pavement_heat_capacity = pavement_subsurface_pars_2(0,pavement_type)
2444	ENDIF
2445
2446	IF ( albedo_type == 9999999 .AND. albedo == 9999999.9_wp ) THEN
2447	albedo_type = INT(pavement_pars(2,pavement_type))
2448	ENDIF
2449
2450	IF ( emissivity == 9999999.9_wp ) THEN
2451	emissivity = pavement_pars(3,pavement_type)
2452	ENDIF
2453
2454	!
2455	!-- If the depth level of the pavement is not set, determine it from
2456	!-- lookup table.
2457	IF ( pavement_depth_level == 0 ) THEN
2458	DO k = nzb_soil, nzt_soil
2459	IF ( pavement_subsurface_pars_1(k,pavement_type) == 9999999.9_wp &
2460	.OR. pavement_subsurface_pars_2(k,pavement_type) &
2461	== 9999999.9_wp) THEN
2462	nzt_pavement = k-1
2463	EXIT
2464	ENDIF
2465	ENDDO
2466	ELSE
2467	nzt_pavement = pavement_depth_level
2468	ENDIF
2469
2470	ENDIF
2471	!
2472	!-- Map values to the respective 2D/3D arrays
2473	!-- Horizontal surfaces
2474	surf_lsm_h%alpha_vg = alpha_vangenuchten
2475	surf_lsm_h%l_vg = l_vangenuchten
2476	surf_lsm_h%n_vg = n_vangenuchten
2477	surf_lsm_h%gamma_w_sat = hydraulic_conductivity
2478	surf_lsm_h%m_sat = saturation_moisture
2479	surf_lsm_h%m_fc = field_capacity
2480	surf_lsm_h%m_wilt = wilting_point
2481	surf_lsm_h%m_res = residual_moisture
2482	surf_lsm_h%r_soil_min = min_soil_resistance
2483	!
2484	!-- Vertical surfaces
2485	DO l = 0, 3
2486	surf_lsm_v(l)%alpha_vg = alpha_vangenuchten
2487	surf_lsm_v(l)%l_vg = l_vangenuchten
2488	surf_lsm_v(l)%n_vg = n_vangenuchten
2489	surf_lsm_v(l)%gamma_w_sat = hydraulic_conductivity
2490	surf_lsm_v(l)%m_sat = saturation_moisture
2491	surf_lsm_v(l)%m_fc = field_capacity
2492	surf_lsm_v(l)%m_wilt = wilting_point
2493	surf_lsm_v(l)%m_res = residual_moisture
2494	surf_lsm_v(l)%r_soil_min = min_soil_resistance
2495	ENDDO
2496
2497	!
2498	!-- Initial run actions
2499	IF ( TRIM( initializing_actions ) /= 'read_restart_data' ) THEN
2500	!
2501	!-- First, for horizontal surfaces
2502	!
2503	!-- Set artifical values for ts and us so that r_a has its initial value
2504	!-- for the first time step. Only for interior core domain, not for ghost points
2505	DO m = 1, surf_lsm_h%ns
2506
2507	t_soil_h%var_2d(:,m) = 0.0_wp
2508	m_soil_h%var_2d(:,m) = 0.0_wp
2509	m_liq_h%var_1d(m) = 0.0_wp
2510
2511	!-- Map user settings of T and q for each soil layer
2512	!-- (make sure that the soil moisture does not drop below the permanent
2513	!-- wilting point) -> problems with devision by zero)
2514	IF ( surf_lsm_h%water_surface(m) ) THEN
2515
2516	surf_lsm_h%z0(m) = z0_water
2517	surf_lsm_h%z0h(m) = z0h_water
2518	surf_lsm_h%z0q(m) = z0h_water
2519	t_soil_h%var_2d(:,m) = water_temperature
2520	surf_lsm_h%lambda_surface_s(m) = 1.0E10_wp
2521	surf_lsm_h%lambda_surface_u(m) = 1.0E10_wp
2522	surf_lsm_h%c_surface(m) = 0.0_wp
2523
2524	ELSEIF ( surf_lsm_h%pavement_surface(m) ) THEN
2525
2526	surf_lsm_h%z0(m) = z0_pavement
2527	surf_lsm_h%z0h(m) = z0h_pavement
2528	surf_lsm_h%z0q(m) = z0h_pavement
2529	DO k = nzb_soil, nzt_soil
2530	t_soil_h%var_2d(k,m) = soil_temperature(k)
2531	m_soil_h%var_2d(k,m) = soil_moisture(k)
2532	ENDDO
2533	t_soil_h%var_2d(nzt_soil+1,m) = soil_temperature(nzt_soil+1)
2534	surf_lsm_h%lambda_surface_s(m) = pavement_heat_conduct &
2535	* ddz_soil(nzb_soil) &
2536	* 2.0_wp
2537	surf_lsm_h%lambda_surface_u(m) = surf_lsm_h%lambda_surface_s(m)
2538	surf_lsm_h%c_surface(m) = pavement_heat_capacity &
2539	* dz_soil(nzb_soil) &
2540	* 0.25_wp
2541
2542	IF ( pavement_type /= 0 ) THEN
2543	DO k = nzb_soil, nzt_soil
2544	surf_lsm_h%lambda_h_def(k,m) = pavement_subsurface_pars_1(k,pavement_type)
2545	surf_lsm_h%rho_c_total_def(k,m) = pavement_subsurface_pars_2(k,pavement_type)
2546	ENDDO
2547	ELSE
2548	surf_lsm_h%lambda_h_def(:,m) = pavement_heat_conduct
2549	surf_lsm_h%rho_c_total_def(:,m) = pavement_heat_capacity
2550	ENDIF
2551
2552	ELSEIF ( surf_lsm_h%vegetation_surface(m) ) THEN
2553
2554	surf_lsm_h%z0(m) = z0_vegetation
2555	surf_lsm_h%z0h(m) = z0h_vegetation
2556	surf_lsm_h%z0q(m) = z0h_vegetation
2557	DO k = nzb_soil, nzt_soil
2558	t_soil_h%var_2d(k,m) = soil_temperature(k)
2559	m_soil_h%var_2d(k,m) = soil_moisture(k)
2560	ENDDO
2561	t_soil_h%var_2d(nzt_soil+1,m) = soil_temperature(nzt_soil+1)
2562	surf_lsm_h%lambda_surface_s(m) = lambda_surface_stable
2563	surf_lsm_h%lambda_surface_u(m) = lambda_surface_unstable
2564	surf_lsm_h%c_surface(m) = c_surface
2565
2566	ENDIF
2567
2568	i = surf_lsm_h%i(m)
2569	j = surf_lsm_h%j(m)
2570	k = surf_lsm_h%k(m)
2571	!
2572	!-- Calculate surface temperature. In case of bare soil, the surface
2573	!-- temperature must be reset to the soil temperature in the first soil
2574	!-- layer
2575	IF ( surf_lsm_h%lambda_surface_s(m) == 0.0_wp ) THEN
2576	t_surface_h%var_1d(:) = t_soil_h%var_2d(nzb_soil,m)
2577	surf_lsm_h%pt_surface(:) = t_soil_h%var_2d(nzb_soil,m) / exn
2578	ELSE
2579	t_surface_h%var_1d(:) = pt_surface * exn
2580	surf_lsm_h%pt_surface(:) = pt_surface
2581	ENDIF
2582
2583	IF ( cloud_physics .OR. cloud_droplets ) THEN
2584	pt1 = pt(k,j,i) + l_d_cp * pt_d_t(k) * ql(k,j,i)
2585	ELSE
2586	pt1 = pt(k,j,i)
2587	ENDIF
2588
2589	!
2590	!-- Assure that r_a cannot be zero at model start
2591	IF ( pt1 == surf_lsm_h%pt_surface(m) ) pt1 = pt1 + 1.0E-10_wp
2592
2593	surf_lsm_h%us(m) = 0.1_wp
2594	surf_lsm_h%ts(m) = ( pt1 - surf_lsm_h%pt_surface(m) ) &
2595	/ surf_lsm_h%r_a(m)
2596	surf_lsm_h%shf(m) = - surf_lsm_h%us(m) * surf_lsm_h%ts(m) &
2597	* rho_surface
2598
2599	ENDDO
2600
2601	!
2602	!-- Vertical surfaces
2603	DO l = 0, 3
2604	DO m = 1, surf_lsm_v(l)%ns
2605
2606	t_soil_v(l)%var_2d = 0.0_wp
2607	m_soil_v(l)%var_2d = 0.0_wp
2608	m_liq_v(l)%var_1d = 0.0_wp
2609
2610
2611	!-- Map user settings of T and q for each soil layer
2612	!-- (make sure that the soil moisture does not drop below the permanent
2613	!-- wilting point) -> problems with devision by zero)
2614	IF ( surf_lsm_v(l)%water_surface(m) ) THEN
2615	surf_lsm_v(l)%z0(m) = z0_water
2616	surf_lsm_v(l)%z0h(m) = z0h_water
2617	surf_lsm_v(l)%z0q(m) = z0h_water
2618	t_soil_v(l)%var_2d(:,m) = water_temperature
2619	surf_lsm_v(l)%lambda_surface_s(m) = 1.0E10_wp
2620	surf_lsm_v(l)%lambda_surface_u(m) = 1.0E10_wp
2621	surf_lsm_v(l)%c_surface(m) = 0.0_wp
2622
2623	ELSEIF ( surf_lsm_v(l)%pavement_surface(m) ) THEN
2624	surf_lsm_v(l)%z0(m) = z0_pavement
2625	surf_lsm_v(l)%z0h(m) = z0h_pavement
2626	surf_lsm_v(l)%z0q(m) = z0h_pavement
2627	DO k = nzb_soil, nzt_soil
2628	t_soil_v(l)%var_2d(k,m) = soil_temperature(k)
2629	m_soil_v(l)%var_2d(k,m) = soil_moisture(k)
2630	ENDDO
2631	t_soil_v(l)%var_2d(nzt_soil+1,m) = soil_temperature(nzt_soil+1)
2632	surf_lsm_v(l)%lambda_surface_s(m) = pavement_heat_conduct &
2633	* ddz_soil(nzb_soil) &
2634	* 2.0_wp
2635	surf_lsm_v(l)%lambda_surface_u(m) = surf_lsm_v(l)%lambda_surface_s(m)
2636	surf_lsm_v(l)%c_surface(m) = pavement_heat_capacity &
2637	* dz_soil(nzb_soil) &
2638	* 0.25_wp
2639	IF ( pavement_type /= 0 ) THEN
2640	DO k = nzb_soil, nzt_soil
2641	surf_lsm_v(l)%lambda_h_def(k,m) = pavement_subsurface_pars_1(k,pavement_type)
2642	surf_lsm_v(l)%rho_c_total_def(k,m) = pavement_subsurface_pars_2(k,pavement_type)
2643	ENDDO
2644	ELSE
2645	surf_lsm_v(l)%lambda_h_def(:,m) = pavement_heat_conduct
2646	surf_lsm_v(l)%rho_c_total_def(:,m) = pavement_heat_capacity
2647	ENDIF
2648
2649
2650	ELSEIF ( surf_lsm_v(l)%vegetation_surface(m) ) THEN
2651
2652	surf_lsm_v(l)%z0(m) = z0_vegetation
2653	surf_lsm_v(l)%z0h(m) = z0h_vegetation
2654	surf_lsm_v(l)%z0q(m) = z0h_vegetation
2655	DO k = nzb_soil, nzt_soil
2656	t_soil_v(l)%var_2d(k,m) = soil_temperature(k)
2657	m_soil_v(l)%var_2d(k,m) = soil_moisture(k)
2658	ENDDO
2659	t_soil_v(l)%var_2d(nzt_soil+1,m) = soil_temperature(nzt_soil+1)
2660	surf_lsm_v(l)%lambda_surface_s(m) = lambda_surface_stable
2661	surf_lsm_v(l)%lambda_surface_u(m) = lambda_surface_unstable
2662	surf_lsm_v(l)%c_surface(m) = c_surface
2663	ENDIF
2664
2665	!
2666	!-- Calculate surface temperature. In case of bare soil, the surface
2667	!-- temperature must be reset to the soil temperature in the first soil
2668	!-- layer
2669	IF ( surf_lsm_v(l)%lambda_surface_s(m) == 0.0_wp ) THEN
2670	t_surface_v(l)%var_1d(:) = t_soil_v(l)%var_2d(nzb_soil,m)
2671	surf_lsm_v(l)%pt_surface(:) = t_soil_v(l)%var_2d(nzb_soil,m) / exn
2672	ELSE
2673	t_surface_v(l)%var_1d(:) = pt_surface * exn
2674	surf_lsm_v(l)%pt_surface(:) = pt_surface
2675	ENDIF
2676
2677	!
2678	!-- Set artifical values for ts and us so that r_a has its initial value
2679	!-- for the first time step. Only for interior core domain, not for ghost points
2680
2681	i = surf_lsm_v(l)%i(m)
2682	j = surf_lsm_v(l)%j(m)
2683	k = surf_lsm_v(l)%k(m)
2684
2685	IF ( cloud_physics .OR. cloud_droplets ) THEN
2686	pt1 = pt(k,j,i) + l_d_cp * pt_d_t(k) * ql(k,j,i)
2687	ELSE
2688	pt1 = pt(k,j,i)
2689	ENDIF
2690
2691	!
2692	!-- Assure that r_a cannot be zero at model start
2693	IF ( pt1 == surf_lsm_v(l)%pt_surface(m) ) pt1 = pt1 + 1.0E-10_wp
2694
2695	surf_lsm_v(l)%us(m) = 0.1_wp
2696	surf_lsm_v(l)%ts(m) = ( pt1 - surf_lsm_v(l)%pt_surface(m) ) / surf_lsm_v(l)%r_a(m)
2697	surf_lsm_v(l)%shf(m) = - surf_lsm_v(l)%us(m) * surf_lsm_v(l)%ts(m) * rho_surface
2698	ENDDO
2699	ENDDO
2700
2701	!
2702	!-- Actions for restart runs
2703	ELSE
2704	!
2705	!-- Horizontal surfaces
2706	DO m = 1, surf_lsm_h%ns
2707	i = surf_lsm_h%i(m)
2708	j = surf_lsm_h%j(m)
2709	k = surf_lsm_h%k(m)
2710	t_surface_h%var_1d(m) = pt(k-1,j,i) * exn
2711	ENDDO
2712	!
2713	!-- Vertical surfaces
2714	DO l = 0, 3
2715	!
2716	!-- Set index offset of surface element, seen from atmospheric grid point
2717	IF ( l == 0 ) THEN
2718	j_off = -1
2719	i_off = 0
2720	ELSEIF ( l == 1 ) THEN
2721	j_off = 1
2722	i_off = 0
2723	ELSEIF ( l == 2 ) THEN
2724	j_off = 0
2725	i_off = -1
2726	ELSEIF ( l == 3 ) THEN
2727	j_off = 0
2728	i_off = 1
2729	ENDIF
2730	DO m = 1, surf_lsm_v(l)%ns
2731	i = surf_lsm_v(l)%i(m)
2732	j = surf_lsm_v(l)%j(m)
2733	k = surf_lsm_v(l)%k(m)
2734	t_surface_v(l)%var_1d(m) = pt(k,j+j_off,i+i_off) * exn
2735	ENDDO
2736	ENDDO
2737
2738	ENDIF
2739	!
2740	!-- Initialize root fraction
2741	!-- Horizontal surfaces
2742	DO m = 1, surf_lsm_h%ns
2743	i = surf_lsm_h%i(m)
2744	j = surf_lsm_h%j(m)
2745
2746	DO k = nzb_soil, nzt_soil
2747	surf_lsm_h%root_fr(k,m) = root_fraction(k)
2748	ENDDO
2749	ENDDO
2750	!
2751	!-- Vertical surfaces
2752	DO l = 0, 3
2753	DO m = 1, surf_lsm_v(l)%ns
2754	i = surf_lsm_v(l)%i(m)
2755	j = surf_lsm_v(l)%j(m)
2756
2757	DO k = nzb_soil, nzt_soil
2758	surf_lsm_v(l)%root_fr(k,m) = root_fraction(k)
2759	ENDDO
2760	ENDDO
2761	ENDDO
2762
2763	!
2764	!-- When no root distribution is given by the user, use look-up table to prescribe
2765	!-- the root fraction in the individual soil layers
2766	IF ( ALL( root_fraction == 9999999.9_wp ) ) THEN
2767
2768	!
2769	!-- First, calculate the index bounds for integration
2770	n_soil_layers_total = nzt_soil - nzb_soil + 6
2771	ALLOCATE ( bound(0:n_soil_layers_total) )
2772	ALLOCATE ( bound_root_fr(0:n_soil_layers_total) )
2773
2774	kn = 0
2775	ko = 0
2776	bound(0) = 0.0_wp
2777	DO k = 1, n_soil_layers_total-1
2778	IF ( zs_layer(kn) <= zs_ref(ko) ) THEN
2779	bound(k) = zs_layer(kn)
2780	bound_root_fr(k) = ko
2781	kn = kn + 1
2782	IF ( kn > nzt_soil+1 ) THEN
2783	kn = nzt_soil
2784	ENDIF
2785	ELSE
2786	bound(k) = zs_ref(ko)
2787	bound_root_fr(k) = ko
2788	ko = ko + 1
2789	IF ( ko > 3 ) THEN
2790	ko = 3
2791	ENDIF
2792	ENDIF
2793
2794	ENDDO
2795
2796	!
2797	!-- Integrate over all soil layers based on the four-layer root fraction
2798	kzs = 1
2799	root_fraction = 0.0_wp
2800	DO k = 0, n_soil_layers_total-2
2801	kroot = bound_root_fr(k+1)
2802	root_fraction(kzs-1) = root_fraction(kzs-1) &
2803	+ root_distribution(kroot,vegetation_type) &
2804	/ dz_soil_ref(kroot) * ( bound(k+1) - bound(k) )
2805
2806	IF ( bound(k+1) == zs_layer(kzs-1) ) THEN
2807	kzs = kzs+1
2808	ENDIF
2809	ENDDO
2810
2811
2812	!
2813	!-- Normalize so that the sum of all fractions equals one
2814	root_fraction = root_fraction / SUM(root_fraction)
2815
2816	DEALLOCATE ( bound )
2817	DEALLOCATE ( bound_root_fr )
2818
2819	!
2820	!-- Map calculated root fractions
2821	DO m = 1, surf_lsm_h%ns
2822	i = surf_lsm_h%i(m)
2823	j = surf_lsm_h%j(m)
2824	DO k = nzb_soil, nzt_soil
2825	IF ( surf_lsm_h%pavement_surface(m) .AND. &
2826	k <= nzt_pavement ) THEN
2827	surf_lsm_h%root_fr(k,m) = 0.0_wp
2828	ELSE
2829	surf_lsm_h%root_fr(k,m) = root_fraction(k)
2830	ENDIF
2831	ENDDO
2832
2833	!
2834	!-- Normalize so that the sum = 1. Only relevant when the root
2835	!-- distribution was set to zero due to pavement at some layers.
2836	IF ( SUM( surf_lsm_h%root_fr(:,m) ) > 0.0_wp ) THEN
2837	DO k = nzb_soil, nzt_soil
2838	surf_lsm_h%root_fr(k,m) = surf_lsm_h%root_fr(k,m) &
2839	/ SUM( surf_lsm_h%root_fr(:,m) )
2840	ENDDO
2841	ENDIF
2842
2843	ENDDO
2844	DO l = 0, 3
2845	DO m = 1, surf_lsm_v(l)%ns
2846	i = surf_lsm_v(l)%i(m)
2847	j = surf_lsm_v(l)%j(m)
2848
2849	DO k = nzb_soil, nzt_soil
2850	IF ( surf_lsm_v(l)%pavement_surface(m) .AND. &
2851	k <= nzt_pavement ) THEN
2852	surf_lsm_v(l)%root_fr(k,m) = 0.0_wp
2853	ELSE
2854	surf_lsm_v(l)%root_fr(k,m) = root_fraction(k)
2855	ENDIF
2856	ENDDO
2857	!
2858	!-- Normalize so that the sum = 1. Only relevant when the root
2859	!-- distribution was set to zero due to pavement at some layers.
2860	IF ( SUM( surf_lsm_v(l)%root_fr(:,m) ) > 0.0_wp ) THEN
2861	DO k = nzb_soil, nzt_soil
2862	surf_lsm_v(l)%root_fr(k,m) = surf_lsm_v(l)%root_fr(k,m) &
2863	/ SUM( surf_lsm_v(l)%root_fr(:,m) )
2864	ENDDO
2865	ENDIF
2866
2867	ENDDO
2868	ENDDO
2869	ENDIF
2870	!
2871	!-- Map vegetation parameters to the respective 2D arrays
2872	surf_lsm_h%r_canopy_min = min_canopy_resistance
2873	surf_lsm_h%lai = leaf_area_index
2874	surf_lsm_h%c_veg = vegetation_coverage
2875	surf_lsm_h%g_d = canopy_resistance_coefficient
2876	surf_lsm_h%f_sw_in = f_shortwave_incoming
2877
2878	!-- Vertical surfaces
2879	DO l = 0, 3
2880
2881	!
2882	!-- Map vegetation parameters to the respective 2D arrays
2883	surf_lsm_v(l)%r_canopy_min = min_canopy_resistance
2884	surf_lsm_v(l)%lai = leaf_area_index
2885	surf_lsm_v(l)%c_veg = vegetation_coverage
2886	surf_lsm_v(l)%g_d = canopy_resistance_coefficient
2887	surf_lsm_v(l)%f_sw_in = f_shortwave_incoming
2888	ENDDO
2889
2890	!
2891	!-- Possibly do user-defined actions (e.g. define heterogeneous land surface)
2892	CALL user_init_land_surface
2893
2894
2895	!
2896	!-- Calculate new roughness lengths (for water surfaces only, i.e. only
2897	!- horizontal surfaces)
2898	CALL calc_z0_water_surface
2899
2900	t_soil_h_p = t_soil_h
2901	m_soil_h_p = m_soil_h
2902	m_liq_h_p = m_liq_h
2903	t_surface_h_p = t_surface_h
2904
2905	t_soil_v_p = t_soil_v
2906	m_soil_v_p = m_soil_v
2907	m_liq_v_p = m_liq_v
2908	t_surface_v_p = t_surface_v
2909
2910
2911
2912	!-- Store initial profiles of t_soil and m_soil (assuming they are
2913	!-- horizontally homogeneous on this PE)
2914	hom(nzb_soil:nzt_soil,1,90,:) = SPREAD( t_soil_h%var_2d(nzb_soil:nzt_soil,1), &
2915	2, statistic_regions+1 )
2916	hom(nzb_soil:nzt_soil,1,92,:) = SPREAD( m_soil_h%var_2d(nzb_soil:nzt_soil,1), &
2917	2, statistic_regions+1 )
2918
2919	END SUBROUTINE lsm_init
2920
2921
2922	!------------------------------------------------------------------------------!
2923	! Description:
2924	! ------------
2925	!> Allocate land surface model arrays and define pointers
2926	!------------------------------------------------------------------------------!
2927	SUBROUTINE lsm_init_arrays
2928
2929
2930	IMPLICIT NONE
2931
2932	INTEGER(iwp) :: l !< index indicating facing of surface array
2933
2934	ALLOCATE ( root_extr(nzb_soil:nzt_soil) )
2935	root_extr = 0.0_wp
2936
2937	!
2938	!-- Allocate surface and soil temperature / humidity. Please note,
2939	!-- these arrays are allocated according to surface-data structure,
2940	!-- even if they do not belong to the data type due to the
2941	!-- pointer arithmetric (TARGET attribute is not allowed in a data-type).
2942	#if defined( __nopointer )
2943	!
2944	!-- Horizontal surfaces
2945	ALLOCATE ( m_liq_h%var_1d(1:surf_lsm_h%ns) )
2946	ALLOCATE ( m_liq_h_p%var_1d(1:surf_lsm_h%ns) )
2947	ALLOCATE ( t_surface_h%var_1d(1:surf_lsm_h%ns) )
2948	ALLOCATE ( t_surface_h_p%var_1d(1:surf_lsm_h%ns) )
2949	ALLOCATE ( m_soil_h%var_2d(nzb_soil:nzt_soil,1:surf_lsm_h%ns) )
2950	ALLOCATE ( m_soil_h_p%var_2d(nzb_soil:nzt_soil,1:surf_lsm_h%ns) )
2951	ALLOCATE ( t_soil_h%var_2d(nzb_soil:nzt_soil+1,1:surf_lsm_h%ns) )
2952	ALLOCATE ( t_soil_h_p%var_2d(nzb_soil:nzt_soil+1,1:surf_lsm_h%ns) )
2953
2954	!
2955	!-- Vertical surfaces
2956	DO l = 0, 3
2957	ALLOCATE ( m_liq_v(l)%var_1d(1:surf_lsm_v(l)%ns) )
2958	ALLOCATE ( m_liq_v_p(l)%var_1d(1:surf_lsm_v(l)%ns) )
2959	ALLOCATE ( t_surface_v(l)%var_1d(1:surf_lsm_v(l)%ns) )
2960	ALLOCATE ( t_surface_v_p(l)%var_1d(1:surf_lsm_v(l)%ns) )
2961	ALLOCATE ( m_soil_v(l)%var_2d(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns) )
2962	ALLOCATE ( m_soil_v_p(l)%var_2d(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns) )
2963	ALLOCATE ( t_soil_v(l)%var_2d(nzb_soil:nzt_soil+1,1:surf_lsm_v(l)%ns) )
2964	ALLOCATE ( t_soil_v_p(l)%var_2d(nzb_soil:nzt_soil+1,1:surf_lsm_v(l)%ns) )
2965	ENDDO
2966	#else
2967	!
2968	!-- Horizontal surfaces
2969	ALLOCATE ( m_liq_h_1%var_1d(1:surf_lsm_h%ns) )
2970	ALLOCATE ( m_liq_h_2%var_1d(1:surf_lsm_h%ns) )
2971	ALLOCATE ( t_surface_h_1%var_1d(1:surf_lsm_h%ns) )
2972	ALLOCATE ( t_surface_h_2%var_1d(1:surf_lsm_h%ns) )
2973	ALLOCATE ( m_soil_h_1%var_2d(nzb_soil:nzt_soil,1:surf_lsm_h%ns) )
2974	ALLOCATE ( m_soil_h_2%var_2d(nzb_soil:nzt_soil,1:surf_lsm_h%ns) )
2975	ALLOCATE ( t_soil_h_1%var_2d(nzb_soil:nzt_soil+1,1:surf_lsm_h%ns) )
2976	ALLOCATE ( t_soil_h_2%var_2d(nzb_soil:nzt_soil+1,1:surf_lsm_h%ns) )
2977	!
2978	!-- Vertical surfaces
2979	DO l = 0, 3
2980	ALLOCATE ( m_liq_v_1(l)%var_1d(1:surf_lsm_v(l)%ns) )
2981	ALLOCATE ( m_liq_v_2(l)%var_1d(1:surf_lsm_v(l)%ns) )
2982	ALLOCATE ( t_surface_v_1(l)%var_1d(1:surf_lsm_v(l)%ns) )
2983	ALLOCATE ( t_surface_v_2(l)%var_1d(1:surf_lsm_v(l)%ns) )
2984	ALLOCATE ( m_soil_v_1(l)%var_2d(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns) )
2985	ALLOCATE ( m_soil_v_2(l)%var_2d(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns) )
2986	ALLOCATE ( t_soil_v_1(l)%var_2d(nzb_soil:nzt_soil+1,1:surf_lsm_v(l)%ns) )
2987	ALLOCATE ( t_soil_v_2(l)%var_2d(nzb_soil:nzt_soil+1,1:surf_lsm_v(l)%ns) )
2988	ENDDO
2989	#endif
2990
2991	!
2992	!-- Allocate intermediate timestep arrays
2993	!-- Horizontal surfaces
2994	ALLOCATE ( tm_liq_h_m%var_1d(1:surf_lsm_h%ns) )
2995	ALLOCATE ( tt_surface_h_m%var_1d(1:surf_lsm_h%ns) )
2996	ALLOCATE ( tm_soil_h_m%var_2d(nzb_soil:nzt_soil,1:surf_lsm_h%ns) )
2997	ALLOCATE ( tt_soil_h_m%var_2d(nzb_soil:nzt_soil,1:surf_lsm_h%ns) )
2998	!
2999	!-- Horizontal surfaces
3000	DO l = 0, 3
3001	ALLOCATE ( tm_liq_v_m(l)%var_1d(1:surf_lsm_v(l)%ns) )
3002	ALLOCATE ( tt_surface_v_m(l)%var_1d(1:surf_lsm_v(l)%ns) )
3003	ALLOCATE ( tm_soil_v_m(l)%var_2d(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns) )
3004	ALLOCATE ( tt_soil_v_m(l)%var_2d(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns) )
3005	ENDDO
3006	!
3007	!-- Allocate skin-surface temperature
3008	ALLOCATE ( surf_lsm_h%pt_surface(1:surf_lsm_h%ns) )
3009	DO l = 0, 3
3010	ALLOCATE ( surf_lsm_v(l)%pt_surface(1:surf_lsm_v(l)%ns) )
3011	ENDDO
3012
3013	!
3014	!-- Allocate 2D vegetation model arrays
3015	!-- Horizontal surfaces
3016	ALLOCATE ( surf_lsm_h%building_surface(1:surf_lsm_h%ns) )
3017	ALLOCATE ( surf_lsm_h%c_liq(1:surf_lsm_h%ns) )
3018	ALLOCATE ( surf_lsm_h%c_surface(1:surf_lsm_h%ns) )
3019	ALLOCATE ( surf_lsm_h%c_veg(1:surf_lsm_h%ns) )
3020	ALLOCATE ( surf_lsm_h%f_sw_in(1:surf_lsm_h%ns) )
3021	ALLOCATE ( surf_lsm_h%ghf(1:surf_lsm_h%ns) )
3022	ALLOCATE ( surf_lsm_h%g_d(1:surf_lsm_h%ns) )
3023	ALLOCATE ( surf_lsm_h%lai(1:surf_lsm_h%ns) )
3024	ALLOCATE ( surf_lsm_h%lambda_surface_u(1:surf_lsm_h%ns) )
3025	ALLOCATE ( surf_lsm_h%lambda_surface_s(1:surf_lsm_h%ns) )
3026	ALLOCATE ( surf_lsm_h%vegetation_surface(1:surf_lsm_h%ns) )
3027	ALLOCATE ( surf_lsm_h%pavement_surface(1:surf_lsm_h%ns) )
3028	ALLOCATE ( surf_lsm_h%qsws_soil(1:surf_lsm_h%ns) )
3029	ALLOCATE ( surf_lsm_h%qsws_liq(1:surf_lsm_h%ns) )
3030	ALLOCATE ( surf_lsm_h%qsws_veg(1:surf_lsm_h%ns) )
3031	ALLOCATE ( surf_lsm_h%rad_net_l(1:surf_lsm_h%ns) )
3032	ALLOCATE ( surf_lsm_h%r_a(1:surf_lsm_h%ns) )
3033	ALLOCATE ( surf_lsm_h%r_canopy(1:surf_lsm_h%ns) )
3034	ALLOCATE ( surf_lsm_h%r_soil(1:surf_lsm_h%ns) )
3035	ALLOCATE ( surf_lsm_h%r_soil_min(1:surf_lsm_h%ns) )
3036	ALLOCATE ( surf_lsm_h%r_s(1:surf_lsm_h%ns) )
3037	ALLOCATE ( surf_lsm_h%r_canopy_min(1:surf_lsm_h%ns) )
3038	ALLOCATE ( surf_lsm_h%water_surface(1:surf_lsm_h%ns) )
3039
3040	surf_lsm_h%water_surface = .FALSE.
3041	surf_lsm_h%pavement_surface = .FALSE.
3042	surf_lsm_h%vegetation_surface = .FALSE.
3043	!
3044	!-- Vertical surfaces
3045	DO l = 0, 3
3046	ALLOCATE ( surf_lsm_v(l)%building_surface(1:surf_lsm_v(l)%ns) )
3047	ALLOCATE ( surf_lsm_v(l)%c_liq(1:surf_lsm_v(l)%ns) )
3048	ALLOCATE ( surf_lsm_v(l)%c_surface(1:surf_lsm_v(l)%ns) )
3049	ALLOCATE ( surf_lsm_v(l)%c_veg(1:surf_lsm_v(l)%ns) )
3050	ALLOCATE ( surf_lsm_v(l)%f_sw_in(1:surf_lsm_v(l)%ns) )
3051	ALLOCATE ( surf_lsm_v(l)%ghf(1:surf_lsm_v(l)%ns) )
3052	ALLOCATE ( surf_lsm_v(l)%g_d(1:surf_lsm_v(l)%ns) )
3053	ALLOCATE ( surf_lsm_v(l)%lai(1:surf_lsm_v(l)%ns) )
3054	ALLOCATE ( surf_lsm_v(l)%lambda_surface_u(1:surf_lsm_v(l)%ns) )
3055	ALLOCATE ( surf_lsm_v(l)%lambda_surface_s(1:surf_lsm_v(l)%ns) )
3056	ALLOCATE ( surf_lsm_v(l)%vegetation_surface(1:surf_lsm_v(l)%ns) )
3057	ALLOCATE ( surf_lsm_v(l)%pavement_surface(1:surf_lsm_v(l)%ns) )
3058	ALLOCATE ( surf_lsm_v(l)%qsws_soil(1:surf_lsm_v(l)%ns) )
3059	ALLOCATE ( surf_lsm_v(l)%qsws_liq(1:surf_lsm_v(l)%ns) )
3060	ALLOCATE ( surf_lsm_v(l)%qsws_veg(1:surf_lsm_v(l)%ns) )
3061	ALLOCATE ( surf_lsm_v(l)%rad_net_l(1:surf_lsm_v(l)%ns) )
3062	ALLOCATE ( surf_lsm_v(l)%r_a(1:surf_lsm_v(l)%ns) )
3063	ALLOCATE ( surf_lsm_v(l)%r_canopy(1:surf_lsm_v(l)%ns) )
3064	ALLOCATE ( surf_lsm_v(l)%r_soil(1:surf_lsm_v(l)%ns) )
3065	ALLOCATE ( surf_lsm_v(l)%r_soil_min(1:surf_lsm_v(l)%ns) )
3066	ALLOCATE ( surf_lsm_v(l)%r_s(1:surf_lsm_v(l)%ns) )
3067	ALLOCATE ( surf_lsm_v(l)%r_canopy_min(1:surf_lsm_v(l)%ns) )
3068	ALLOCATE ( surf_lsm_v(l)%water_surface(1:surf_lsm_v(l)%ns) )
3069
3070	surf_lsm_v(l)%water_surface = .FALSE.
3071	surf_lsm_v(l)%pavement_surface = .FALSE.
3072	surf_lsm_v(l)%vegetation_surface = .FALSE.
3073
3074	ENDDO
3075
3076	#if ! defined( __nopointer )
3077	!
3078	!-- Initial assignment of the pointers
3079	!-- Horizontal surfaces
3080	t_soil_h => t_soil_h_1; t_soil_h_p => t_soil_h_2
3081	t_surface_h => t_surface_h_1; t_surface_h_p => t_surface_h_2
3082	m_soil_h => m_soil_h_1; m_soil_h_p => m_soil_h_2
3083	m_liq_h => m_liq_h_1; m_liq_h_p => m_liq_h_2
3084	!
3085	!-- Vertical surfaces
3086	t_soil_v => t_soil_v_1; t_soil_v_p => t_soil_v_2
3087	t_surface_v => t_surface_v_1; t_surface_v_p => t_surface_v_2
3088	m_soil_v => m_soil_v_1; m_soil_v_p => m_soil_v_2
3089	m_liq_v => m_liq_v_1; m_liq_v_p => m_liq_v_2
3090	#endif
3091
3092
3093	END SUBROUTINE lsm_init_arrays
3094
3095
3096	!------------------------------------------------------------------------------!
3097	! Description:
3098	! ------------
3099	!> Parin for &lsmpar for land surface model
3100	!------------------------------------------------------------------------------!
3101	SUBROUTINE lsm_parin
3102
3103
3104	IMPLICIT NONE
3105
3106	CHARACTER (LEN=80) :: line !< dummy string that contains the current line of the parameter file
3107
3108	NAMELIST /lsm_par/ alpha_vangenuchten, c_surface, &
3109	canopy_resistance_coefficient, &
3110	constant_roughness, &
3111	conserve_water_content, &
3112	dz_soil, &
3113	f_shortwave_incoming, field_capacity, &
3114	aero_resist_kray, hydraulic_conductivity, &
3115	lambda_surface_stable, &
3116	lambda_surface_unstable, leaf_area_index, &
3117	l_vangenuchten, min_canopy_resistance, &
3118	min_soil_resistance, n_vangenuchten, &
3119	pavement_depth_level, &
3120	pavement_heat_capacity, &
3121	pavement_heat_conduct, pavement_type, &
3122	residual_moisture, root_fraction, &
3123	saturation_moisture, skip_time_do_lsm, &
3124	soil_moisture, soil_temperature, soil_type, &
3125	surface_type, &
3126	vegetation_coverage, vegetation_type, &
3127	water_temperature, water_type, &
3128	wilting_point, z0_vegetation, &
3129	z0h_vegetation, z0q_vegetation, z0_water, &
3130	z0h_water, z0q_water, z0_pavement, &
3131	z0h_pavement, z0q_pavement
3132
3133	line = ' '
3134
3135	!
3136	!-- Try to find land surface model package
3137	REWIND ( 11 )
3138	line = ' '
3139	DO WHILE ( INDEX( line, '&lsm_par' ) == 0 )
3140	READ ( 11, '(A)', END=10 ) line
3141	ENDDO
3142	BACKSPACE ( 11 )
3143
3144	!
3145	!-- Read user-defined namelist
3146	READ ( 11, lsm_par )
3147
3148	!
3149	!-- Set flag that indicates that the land surface model is switched on
3150	land_surface = .TRUE.
3151
3152	!
3153	!-- Activate spinup
3154	IF ( spinup_time > 0.0_wp ) THEN
3155	coupling_start_time = spinup_time
3156	end_time = end_time + spinup_time
3157	IF ( spinup_pt_mean == 9999999.9_wp ) THEN
3158	spinup_pt_mean = pt_surface
3159	ENDIF
3160	spinup = .TRUE.
3161	ENDIF
3162
3163
3164	10 CONTINUE
3165
3166
3167	END SUBROUTINE lsm_parin
3168
3169
3170	!------------------------------------------------------------------------------!
3171	! Description:
3172	! ------------
3173	!> Soil model as part of the land surface model. The model predicts soil
3174	!> temperature and water content.
3175	!------------------------------------------------------------------------------!
3176	SUBROUTINE lsm_soil_model( horizontal, l, calc_soil_moisture )
3177
3178
3179	IMPLICIT NONE
3180
3181	INTEGER(iwp) :: k !< running index
3182	INTEGER(iwp) :: l !< surface-data type index indication facing
3183	INTEGER(iwp) :: m !< running index
3184
3185	LOGICAL, INTENT(IN) :: calc_soil_moisture !< flag indicating whether soil moisture shall be calculated or not.
3186
3187	LOGICAL :: horizontal !< flag indication horizontal wall, required to set pointer accordingly
3188
3189	REAL(wp) :: h_vg !< Van Genuchten coef. h
3190
3191	REAL(wp), DIMENSION(nzb_soil:nzt_soil) :: gamma_temp, & !< temp. gamma
3192	lambda_temp, & !< temp. lambda
3193	tend !< tendency
3194
3195	TYPE(surf_type_lsm), POINTER :: surf_m_soil
3196	TYPE(surf_type_lsm), POINTER :: surf_m_soil_p
3197	TYPE(surf_type_lsm), POINTER :: surf_t_soil
3198	TYPE(surf_type_lsm), POINTER :: surf_t_soil_p
3199	TYPE(surf_type_lsm), POINTER :: surf_tm_soil_m
3200	TYPE(surf_type_lsm), POINTER :: surf_tt_soil_m
3201
3202	TYPE(surf_type), POINTER :: surf !< surface-date type variable
3203
3204	IF ( horizontal ) THEN
3205	surf => surf_lsm_h
3206	surf_m_soil => m_soil_h
3207	surf_m_soil_p => m_soil_h_p
3208	surf_t_soil => t_soil_h
3209	surf_t_soil_p => t_soil_h_p
3210	surf_tm_soil_m => tm_soil_h_m
3211	surf_tt_soil_m => tt_soil_h_m
3212	ELSE
3213	surf => surf_lsm_v(l)
3214	surf_m_soil => m_soil_v(l)
3215	surf_m_soil_p => m_soil_v_p(l)
3216	surf_t_soil => t_soil_v(l)
3217	surf_t_soil_p => t_soil_v_p(l)
3218	surf_tm_soil_m => tm_soil_v_m(l)
3219	surf_tt_soil_m => tt_soil_v_m(l)
3220	ENDIF
3221
3222	DO m = 1, surf%ns
3223
3224	IF ( .NOT. surf%water_surface(m) ) THEN
3225	DO k = nzb_soil, nzt_soil
3226
3227	IF ( surf%pavement_surface(m) .AND. k <= nzt_pavement ) THEN
3228
3229	surf%rho_c_total(k,m) = surf%rho_c_total_def(k,m)
3230	lambda_temp(k) = surf%lambda_h_def(k,m)
3231
3232	ELSE
3233	!
3234	!-- Calculate volumetric heat capacity of the soil, taking
3235	!-- into account water content
3236	surf%rho_c_total(k,m) = (rho_c_soil * &
3237	( 1.0_wp - surf%m_sat(k,m) ) &
3238	+ rho_c_water * surf_m_soil%var_2d(k,m) )
3239
3240	!
3241	!-- Calculate soil heat conductivity at the center of the soil
3242	!-- layers
3243	lambda_h_sat = lambda_h_sm*(1.0_wp - surf%m_sat(k,m)) &
3244	lambda_h_water ** surf_m_soil%var_2d(k,m)
3245
3246	ke = 1.0_wp + LOG10( MAX( 0.1_wp, surf_m_soil%var_2d(k,m) / &
3247	surf%m_sat(k,m) ) )
3248
3249	lambda_temp(k) = ke * (lambda_h_sat - lambda_h_dry) + &
3250	lambda_h_dry
3251	ENDIF
3252	ENDDO
3253
3254	!
3255	!-- Calculate soil heat conductivity (lambda_h) at the _layer level
3256	!-- using linear interpolation. For pavement surface, the
3257	!-- true pavement depth is considered
3258	DO k = nzb_soil, nzt_soil-1
3259	surf%lambda_h(k,m) = ( lambda_temp(k+1) + lambda_temp(k) ) &
3260	* 0.5_wp
3261	ENDDO
3262	surf%lambda_h(nzt_soil,m) = lambda_temp(nzt_soil)
3263
3264	!
3265	!-- Prognostic equation for soil temperature t_soil
3266	tend(:) = 0.0_wp
3267
3268	tend(nzb_soil) = ( 1.0_wp / surf%rho_c_total(nzb_soil,m) ) * &
3269	( surf%lambda_h(nzb_soil,m) &
3270	* ( surf_t_soil%var_2d(nzb_soil+1,m) &
3271	- surf_t_soil%var_2d(nzb_soil,m) ) &
3272	* ddz_soil_center(nzb_soil) + surf%ghf(m) ) &
3273	* ddz_soil(nzb_soil)
3274
3275	DO k = nzb_soil+1, nzt_soil
3276	tend(k) = ( 1.0_wp / surf%rho_c_total(k,m) ) &
3277	* ( surf%lambda_h(k,m) &
3278	* ( surf_t_soil%var_2d(k+1,m) - surf_t_soil%var_2d(k,m) ) &
3279	* ddz_soil_center(k) - surf%lambda_h(k-1,m) &
3280	* ( surf_t_soil%var_2d(k,m) - surf_t_soil%var_2d(k-1,m) ) &
3281	* ddz_soil_center(k-1) ) * ddz_soil(k)
3282
3283	ENDDO
3284
3285	surf_t_soil_p%var_2d(nzb_soil:nzt_soil,m) = &
3286	surf_t_soil%var_2d(nzb_soil:nzt_soil,m) &
3287	+ dt_3d * ( tsc(2) &
3288	* tend(nzb_soil:nzt_soil) + tsc(3) &
3289	* surf_tt_soil_m%var_2d(nzb_soil:nzt_soil,m) )
3290
3291	!
3292	!-- Calculate t_soil tendencies for the next Runge-Kutta step
3293	IF ( timestep_scheme(1:5) == 'runge' ) THEN
3294	IF ( intermediate_timestep_count == 1 ) THEN
3295	DO k = nzb_soil, nzt_soil
3296	surf_tt_soil_m%var_2d(k,m) = tend(k)
3297	ENDDO
3298	ELSEIF ( intermediate_timestep_count < &
3299	intermediate_timestep_count_max ) THEN
3300	DO k = nzb_soil, nzt_soil
3301	surf_tt_soil_m%var_2d(k,m) = -9.5625_wp * tend(k) &
3302	+ 5.3125_wp &
3303	* surf_tt_soil_m%var_2d(k,m)
3304	ENDDO
3305	ENDIF
3306	ENDIF
3307
3308
3309	DO k = nzb_soil, nzt_soil
3310
3311	!
3312	!-- In order to prevent water tranport through paved surfaces,
3313	!-- conductivity and diffusivity are set to zero
3314	IF ( surf%pavement_surface(m) .AND. k <= nzt_pavement ) THEN
3315
3316	lambda_temp(k) = 0.0_wp
3317	gamma_temp(k) = 0.0_wp
3318
3319	ELSE
3320
3321	!
3322	!-- Calculate soil diffusivity at the center of the soil layers
3323	lambda_temp(k) = (- b_ch * surf%gamma_w_sat(k,m) * psi_sat &
3324	/ surf%m_sat(k,m) ) * ( &
3325	MAX( surf_m_soil%var_2d(k,m), &
3326	surf%m_wilt(k,m) ) / surf%m_sat(k,m) )**( &
3327	b_ch + 2.0_wp )
3328
3329	!
3330	!-- Parametrization of Van Genuchten
3331	!-- Calculate the hydraulic conductivity after Van Genuchten (1980)
3332	h_vg = ( ( ( surf%m_res(k,m) - surf%m_sat(k,m) ) / &
3333	( surf%m_res(k,m) - &
3334	MAX(surf_m_soil%var_2d(k,m), surf%m_wilt(k,m)) &
3335	) &
3336	)**( surf%n_vg(k,m) / ( surf%n_vg(k,m) - 1.0_wp ) &
3337	) - 1.0_wp &
3338	)**( 1.0_wp / surf%n_vg(k,m) ) / surf%alpha_vg(k,m)
3339
3340	gamma_temp(k) = surf%gamma_w_sat(k,m) * ( ( ( 1.0_wp + &
3341	( surf%alpha_vg(k,m) * h_vg )**surf%n_vg(k,m) &
3342	)**( 1.0_wp - 1.0_wp / surf%n_vg(k,m) ) &
3343	- ( surf%alpha_vg(k,m) * h_vg )**( surf%n_vg(k,m) &
3344	- 1.0_wp) )**2 ) / ( ( 1.0_wp + ( surf%alpha_vg(k,m) &
3345	* h_vg )surf%n_vg(k,m) )(( 1.0_wp - 1.0_wp &
3346	/ surf%n_vg(k,m) ) * ( surf%l_vg(k,m) + 2.0_wp) ))
3347	ENDIF
3348
3349	ENDDO
3350
3351	ENDIF
3352
3353	ENDDO
3354
3355
3356	DO m = 1, surf%ns
3357
3358	IF ( .NOT. surf%water_surface(m) .AND. calc_soil_moisture ) THEN
3359
3360
3361	!
3362	!-- Prognostic equation for soil moisture content. Only performed,
3363	!-- when humidity is enabled in the atmosphere.
3364	IF ( humidity ) THEN
3365	!
3366	!-- Calculate soil diffusivity (lambda_w) at the _layer level
3367	!-- using linear interpolation. To do: replace this with
3368	!-- ECMWF-IFS Eq. 8.81
3369	DO k = nzb_soil, nzt_soil-1
3370	surf%lambda_w(k,m) = ( lambda_temp(k+1) + lambda_temp(k) ) &
3371	* 0.5_wp
3372	surf%gamma_w(k,m) = ( gamma_temp(k+1) + gamma_temp(k) ) &
3373	* 0.5_wp
3374	ENDDO
3375	!
3376	!
3377	!-- In case of a closed bottom (= water content is conserved),
3378	!-- set hydraulic conductivity to zero to that no water will be
3379	!-- lost in the bottom layer. As gamma_w is always a positive value,
3380	!-- it cannot be set to zero in case of purely dry soil since this
3381	!-- would cause accumulation of (non-existing) water in the lowest
3382	!-- soil layer
3383	IF ( conserve_water_content .AND. surf_m_soil%var_2d(nzt_soil,m) /= 0.0_wp ) THEN
3384	surf%gamma_w(nzt_soil,m) = 0.0_wp
3385	ELSE
3386	surf%gamma_w(nzt_soil,m) = gamma_temp(nzt_soil)
3387	ENDIF
3388
3389	!-- The root extraction (= root_extr * qsws_veg / (rho_l
3390	!-- * l_v)) ensures the mass conservation for water. The
3391	!-- transpiration of plants equals the cumulative withdrawals by
3392	!-- the roots in the soil. The scheme takes into account the
3393	!-- availability of water in the soil layers as well as the root
3394	!-- fraction in the respective layer. Layer with moisture below
3395	!-- wilting point will not contribute, which reflects the
3396	!-- preference of plants to take water from moister layers.
3397	!
3398	!-- Calculate the root extraction (ECMWF 7.69, the sum of
3399	!-- root_extr = 1). The energy balance solver guarantees a
3400	!-- positive transpiration, so that there is no need for an
3401	!-- additional check.
3402	m_total = 0.0_wp
3403	DO k = nzb_soil, nzt_soil
3404	IF ( surf_m_soil%var_2d(k,m) > surf%m_wilt(k,m) ) THEN
3405	m_total = m_total + surf%root_fr(k,m) &
3406	* surf_m_soil%var_2d(k,m)
3407	ENDIF
3408	ENDDO
3409	IF ( m_total > 0.0_wp ) THEN
3410	DO k = nzb_soil, nzt_soil
3411	IF ( surf_m_soil%var_2d(k,m) > surf%m_wilt(k,m) ) THEN
3412	root_extr(k) = surf%root_fr(k,m) &
3413	* surf_m_soil%var_2d(k,m) / m_total
3414	ELSE
3415	root_extr(k) = 0.0_wp
3416	ENDIF
3417	ENDDO
3418	ENDIF
3419	!
3420	!-- Prognostic equation for soil water content m_soil_h.
3421	tend(:) = 0.0_wp
3422
3423	tend(nzb_soil) = ( surf%lambda_w(nzb_soil,m) * ( &
3424	surf_m_soil%var_2d(nzb_soil+1,m) &
3425	- surf_m_soil%var_2d(nzb_soil,m) ) &
3426	* ddz_soil_center(nzb_soil) - surf%gamma_w(nzb_soil,m)&
3427	- ( root_extr(nzb_soil) * surf%qsws_veg(m) &
3428	+ surf%qsws_soil(m) ) * drho_l_lv ) &
3429	* ddz_soil(nzb_soil)
3430
3431	DO k = nzb_soil+1, nzt_soil-1
3432	tend(k) = ( surf%lambda_w(k,m) * ( surf_m_soil%var_2d(k+1,m) &
3433	- surf_m_soil%var_2d(k,m) ) * ddz_soil_center(k) &
3434	- surf%gamma_w(k,m) &
3435	- surf%lambda_w(k-1,m) * ( surf_m_soil%var_2d(k,m) &
3436	- surf_m_soil%var_2d(k-1,m)) * ddz_soil_center(k-1) &
3437	+ surf%gamma_w(k-1,m) - (root_extr(k) &
3438	* surf%qsws_veg(m) * drho_l_lv) &
3439	) * ddz_soil(k)
3440	ENDDO
3441	tend(nzt_soil) = ( - surf%gamma_w(nzt_soil,m) &
3442	- surf%lambda_w(nzt_soil-1,m) &
3443	* ( surf_m_soil%var_2d(nzt_soil,m) &
3444	- surf_m_soil%var_2d(nzt_soil-1,m)) &
3445	* ddz_soil_center(nzt_soil-1) &
3446	+ surf%gamma_w(nzt_soil-1,m) - ( &
3447	root_extr(nzt_soil) &
3448	* surf%qsws_veg(m) * drho_l_lv ) &
3449	) * ddz_soil(nzt_soil)
3450
3451	surf_m_soil_p%var_2d(nzb_soil:nzt_soil,m) = &
3452	surf_m_soil%var_2d(nzb_soil:nzt_soil,m) &
3453	+ dt_3d * ( tsc(2) * tend(:) &
3454	+ tsc(3) * surf_tm_soil_m%var_2d(:,m) )
3455	!
3456	!-- Account for dry soils (find a better solution here!)
3457	DO k = nzb_soil, nzt_soil
3458	IF ( surf_m_soil_p%var_2d(k,m) < 0.0_wp ) surf_m_soil_p%var_2d(k,m) = 0.0_wp
3459	ENDDO
3460
3461	!
3462	!-- Calculate m_soil tendencies for the next Runge-Kutta step
3463	IF ( timestep_scheme(1:5) == 'runge' ) THEN
3464	IF ( intermediate_timestep_count == 1 ) THEN
3465	DO k = nzb_soil, nzt_soil
3466	surf_tm_soil_m%var_2d(k,m) = tend(k)
3467	ENDDO
3468	ELSEIF ( intermediate_timestep_count < &
3469	intermediate_timestep_count_max ) THEN
3470	DO k = nzb_soil, nzt_soil
3471	surf_tm_soil_m%var_2d(k,m) = -9.5625_wp * tend(k) &
3472	+ 5.3125_wp &
3473	* surf_tm_soil_m%var_2d(k,m)
3474	ENDDO
3475	ENDIF
3476	ENDIF
3477	ENDIF
3478
3479	ENDIF
3480
3481	ENDDO
3482
3483	END SUBROUTINE lsm_soil_model
3484
3485
3486	!------------------------------------------------------------------------------!
3487	! Description:
3488	! ------------
3489	!> Swapping of timelevels
3490	!------------------------------------------------------------------------------!
3491	SUBROUTINE lsm_swap_timelevel ( mod_count )
3492
3493	IMPLICIT NONE
3494
3495	INTEGER, INTENT(IN) :: mod_count
3496
3497	#if defined( __nopointer )
3498	!
3499	!-- Horizontal surfaces
3500	t_surface_h = t_surface_h_p
3501	t_soil_h = t_soil_h_p
3502	IF ( humidity ) THEN
3503	m_soil_h = m_soil_h_p
3504	m_liq_h = m_liq_h_p
3505	ENDIF
3506	!
3507	!-- Vertical surfaces
3508	t_surface_v = t_surface_v_p
3509	t_soil_v = t_soil_v_p
3510	IF ( humidity ) THEN
3511	m_soil_v = m_soil_v_p
3512	m_liq_v = m_liq_v_p
3513	ENDIF
3514
3515	#else
3516
3517	SELECT CASE ( mod_count )
3518
3519	CASE ( 0 )
3520	!
3521	!-- Horizontal surfaces
3522	t_surface_h => t_surface_h_1; t_surface_h_p => t_surface_h_2
3523	t_soil_h => t_soil_h_1; t_soil_h_p => t_soil_h_2
3524	IF ( humidity ) THEN
3525	m_soil_h => m_soil_h_1; m_soil_h_p => m_soil_h_2
3526	m_liq_h => m_liq_h_1; m_liq_h_p => m_liq_h_2
3527	ENDIF
3528	!
3529	!-- Vertical surfaces
3530	t_surface_v => t_surface_v_1; t_surface_v_p => t_surface_v_2
3531	t_soil_v => t_soil_v_1; t_soil_v_p => t_soil_v_2
3532	IF ( humidity ) THEN
3533	m_soil_v => m_soil_v_1; m_soil_v_p => m_soil_v_2
3534	m_liq_v => m_liq_v_1; m_liq_v_p => m_liq_v_2
3535	ENDIF
3536
3537
3538
3539	CASE ( 1 )
3540	!
3541	!-- Horizontal surfaces
3542	t_surface_h => t_surface_h_2; t_surface_h_p => t_surface_h_1
3543	t_soil_h => t_soil_h_2; t_soil_h_p => t_soil_h_1
3544	IF ( humidity ) THEN
3545	m_soil_h => m_soil_h_2; m_soil_h_p => m_soil_h_1
3546	m_liq_h => m_liq_h_2; m_liq_h_p => m_liq_h_1
3547	ENDIF
3548	!
3549	!-- Vertical surfaces
3550	t_surface_v => t_surface_v_2; t_surface_v_p => t_surface_v_1
3551	t_soil_v => t_soil_v_2; t_soil_v_p => t_soil_v_1
3552	IF ( humidity ) THEN
3553	m_soil_v => m_soil_v_2; m_soil_v_p => m_soil_v_1
3554	m_liq_v => m_liq_v_2; m_liq_v_p => m_liq_v_1
3555	ENDIF
3556
3557	END SELECT
3558	#endif
3559
3560	END SUBROUTINE lsm_swap_timelevel
3561
3562
3563
3564
3565	!------------------------------------------------------------------------------!
3566	!
3567	! Description:
3568	! ------------
3569	!> Subroutine for averaging 3D data
3570	!------------------------------------------------------------------------------!
3571	SUBROUTINE lsm_3d_data_averaging( mode, variable )
3572
3573
3574	USE control_parameters
3575
3576	USE indices
3577
3578	USE kinds
3579
3580	IMPLICIT NONE
3581
3582	CHARACTER (LEN=*) :: mode !<
3583	CHARACTER (LEN=*) :: variable !<
3584
3585	INTEGER(iwp) :: i !<
3586	INTEGER(iwp) :: j !<
3587	INTEGER(iwp) :: k !<
3588	INTEGER(iwp) :: m !< running index
3589
3590	IF ( mode == 'allocate' ) THEN
3591
3592	SELECT CASE ( TRIM( variable ) )
3593
3594	CASE ( 'c_liq*' )
3595	IF ( .NOT. ALLOCATED( c_liq_av ) ) THEN
3596	ALLOCATE( c_liq_av(nysg:nyng,nxlg:nxrg) )
3597	ENDIF
3598	c_liq_av = 0.0_wp
3599
3600	CASE ( 'c_soil*' )
3601	IF ( .NOT. ALLOCATED( c_soil_av ) ) THEN
3602	ALLOCATE( c_soil_av(nysg:nyng,nxlg:nxrg) )
3603	ENDIF
3604	c_soil_av = 0.0_wp
3605
3606	CASE ( 'c_veg*' )
3607	IF ( .NOT. ALLOCATED( c_veg_av ) ) THEN
3608	ALLOCATE( c_veg_av(nysg:nyng,nxlg:nxrg) )
3609	ENDIF
3610	c_veg_av = 0.0_wp
3611
3612	CASE ( 'ghf*' )
3613	IF ( .NOT. ALLOCATED( ghf_av ) ) THEN
3614	ALLOCATE( ghf_av(nysg:nyng,nxlg:nxrg) )
3615	ENDIF
3616	ghf_av = 0.0_wp
3617
3618	CASE ( 'lai*' )
3619	IF ( .NOT. ALLOCATED( lai_av ) ) THEN
3620	ALLOCATE( lai_av(nysg:nyng,nxlg:nxrg) )
3621	ENDIF
3622	lai_av = 0.0_wp
3623
3624	CASE ( 'm_liq*' )
3625	IF ( .NOT. ALLOCATED( m_liq_av ) ) THEN
3626	ALLOCATE( m_liq_av(nysg:nyng,nxlg:nxrg) )
3627	ENDIF
3628	m_liq_av = 0.0_wp
3629
3630	CASE ( 'm_soil' )
3631	IF ( .NOT. ALLOCATED( m_soil_av ) ) THEN
3632	ALLOCATE( m_soil_av(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
3633	ENDIF
3634	m_soil_av = 0.0_wp
3635
3636	CASE ( 'qsws_liq*' )
3637	IF ( .NOT. ALLOCATED( qsws_liq_av ) ) THEN
3638	ALLOCATE( qsws_liq_av(nysg:nyng,nxlg:nxrg) )
3639	ENDIF
3640	qsws_liq_av = 0.0_wp
3641
3642	CASE ( 'qsws_soil*' )
3643	IF ( .NOT. ALLOCATED( qsws_soil_av ) ) THEN
3644	ALLOCATE( qsws_soil_av(nysg:nyng,nxlg:nxrg) )
3645	ENDIF
3646	qsws_soil_av = 0.0_wp
3647
3648	CASE ( 'qsws_veg*' )
3649	IF ( .NOT. ALLOCATED( qsws_veg_av ) ) THEN
3650	ALLOCATE( qsws_veg_av(nysg:nyng,nxlg:nxrg) )
3651	ENDIF
3652	qsws_veg_av = 0.0_wp
3653
3654	CASE ( 'r_a*' )
3655	IF ( .NOT. ALLOCATED( r_a_av ) ) THEN
3656	ALLOCATE( r_a_av(nysg:nyng,nxlg:nxrg) )
3657	ENDIF
3658	r_a_av = 0.0_wp
3659
3660	CASE ( 'r_s*' )
3661	IF ( .NOT. ALLOCATED( r_s_av ) ) THEN
3662	ALLOCATE( r_s_av(nysg:nyng,nxlg:nxrg) )
3663	ENDIF
3664	r_s_av = 0.0_wp
3665
3666	CASE ( 't_soil' )
3667	IF ( .NOT. ALLOCATED( t_soil_av ) ) THEN
3668	ALLOCATE( t_soil_av(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
3669	ENDIF
3670	t_soil_av = 0.0_wp
3671
3672	CASE DEFAULT
3673	CONTINUE
3674
3675	END SELECT
3676
3677	ELSEIF ( mode == 'sum' ) THEN
3678
3679	SELECT CASE ( TRIM( variable ) )
3680
3681	CASE ( 'c_liq*' )
3682	DO m = 1, surf_lsm_h%ns
3683	i = surf_lsm_h%i(m)
3684	j = surf_lsm_h%j(m)
3685	c_liq_av(j,i) = c_liq_av(j,i) + surf_lsm_h%c_liq(m)
3686	ENDDO
3687
3688	CASE ( 'c_soil*' )
3689	DO m = 1, surf_lsm_h%ns
3690	i = surf_lsm_h%i(m)
3691	j = surf_lsm_h%j(m)
3692	c_soil_av(j,i) = c_soil_av(j,i) + (1.0 - surf_lsm_h%c_veg(m))
3693	ENDDO
3694
3695	CASE ( 'c_veg*' )
3696	DO m = 1, surf_lsm_h%ns
3697	i = surf_lsm_h%i(m)
3698	j = surf_lsm_h%j(m)
3699	c_veg_av(j,i) = c_veg_av(j,i) + surf_lsm_h%c_veg(m)
3700	ENDDO
3701
3702	CASE ( 'ghf*' )
3703	DO m = 1, surf_lsm_h%ns
3704	i = surf_lsm_h%i(m)
3705	j = surf_lsm_h%j(m)
3706	ghf_av(j,i) = ghf_av(j,i) + surf_lsm_h%ghf(m)
3707	ENDDO
3708
3709	CASE ( 'lai*' )
3710	DO m = 1, surf_lsm_h%ns
3711	i = surf_lsm_h%i(m)
3712	j = surf_lsm_h%j(m)
3713	lai_av(j,i) = lai_av(j,i) + surf_lsm_h%lai(m)
3714	ENDDO
3715
3716	CASE ( 'm_liq*' )
3717	DO m = 1, surf_lsm_h%ns
3718	i = surf_lsm_h%i(m)
3719	j = surf_lsm_h%j(m)
3720	m_liq_av(j,i) = m_liq_av(j,i) + m_liq_h%var_1d(m)
3721	ENDDO
3722
3723	CASE ( 'm_soil' )
3724	DO m = 1, surf_lsm_h%ns
3725	i = surf_lsm_h%i(m)
3726	j = surf_lsm_h%j(m)
3727	DO k = nzb_soil, nzt_soil
3728	m_soil_av(k,j,i) = m_soil_av(k,j,i) + m_soil_h%var_2d(k,m)
3729	ENDDO
3730	ENDDO
3731
3732	CASE ( 'qsws_liq*' )
3733	DO m = 1, surf_lsm_h%ns
3734	i = surf_lsm_h%i(m)
3735	j = surf_lsm_h%j(m)
3736	qsws_liq_av(j,i) = qsws_liq_av(j,i) + &
3737	surf_lsm_h%qsws_liq(m)
3738	ENDDO
3739
3740	CASE ( 'qsws_soil*' )
3741	DO m = 1, surf_lsm_h%ns
3742	i = surf_lsm_h%i(m)
3743	j = surf_lsm_h%j(m)
3744	qsws_soil_av(j,i) = qsws_soil_av(j,i) + &
3745	surf_lsm_h%qsws_soil(m)
3746	ENDDO
3747
3748	CASE ( 'qsws_veg*' )
3749	DO m = 1, surf_lsm_h%ns
3750	i = surf_lsm_h%i(m)
3751	j = surf_lsm_h%j(m)
3752	qsws_veg_av(j,i) = qsws_veg_av(j,i) + &
3753	surf_lsm_h%qsws_veg(m)
3754	ENDDO
3755
3756	CASE ( 'r_a*' )
3757	DO m = 1, surf_lsm_h%ns
3758	i = surf_lsm_h%i(m)
3759	j = surf_lsm_h%j(m)
3760	r_a_av(j,i) = r_a_av(j,i) + surf_lsm_h%r_a(m)
3761	ENDDO
3762
3763	CASE ( 'r_s*' )
3764	DO m = 1, surf_lsm_h%ns
3765	i = surf_lsm_h%i(m)
3766	j = surf_lsm_h%j(m)
3767	r_s_av(j,i) = r_s_av(j,i) + surf_lsm_h%r_s(m)
3768	ENDDO
3769
3770	CASE ( 't_soil' )
3771	DO m = 1, surf_lsm_h%ns
3772	i = surf_lsm_h%i(m)
3773	j = surf_lsm_h%j(m)
3774	DO k = nzb_soil, nzt_soil
3775	t_soil_av(k,j,i) = t_soil_av(k,j,i) + t_soil_h%var_2d(k,m)
3776	ENDDO
3777	ENDDO
3778
3779	CASE DEFAULT
3780	CONTINUE
3781
3782	END SELECT
3783
3784	ELSEIF ( mode == 'average' ) THEN
3785
3786	SELECT CASE ( TRIM( variable ) )
3787
3788	CASE ( 'c_liq*' )
3789	DO i = nxl, nxr
3790	DO j = nys, nyn
3791	c_liq_av(j,i) = c_liq_av(j,i) &
3792	/ REAL( average_count_3d, KIND=wp )
3793	ENDDO
3794	ENDDO
3795
3796	CASE ( 'c_soil*' )
3797	DO i = nxl, nxr
3798	DO j = nys, nyn
3799	c_soil_av(j,i) = c_soil_av(j,i) &
3800	/ REAL( average_count_3d, KIND=wp )
3801	ENDDO
3802	ENDDO
3803
3804	CASE ( 'c_veg*' )
3805	DO i = nxl, nxr
3806	DO j = nys, nyn
3807	c_veg_av(j,i) = c_veg_av(j,i) &
3808	/ REAL( average_count_3d, KIND=wp )
3809	ENDDO
3810	ENDDO
3811
3812	CASE ( 'ghf*' )
3813	DO i = nxl, nxr
3814	DO j = nys, nyn
3815	ghf_av(j,i) = ghf_av(j,i) &
3816	/ REAL( average_count_3d, KIND=wp )
3817	ENDDO
3818	ENDDO
3819
3820	CASE ( 'lai*' )
3821	DO i = nxl, nxr
3822	DO j = nys, nyn
3823	lai_av(j,i) = lai_av(j,i) &
3824	/ REAL( average_count_3d, KIND=wp )
3825	ENDDO
3826	ENDDO
3827
3828	CASE ( 'm_liq*' )
3829	DO i = nxl, nxr
3830	DO j = nys, nyn
3831	m_liq_av(j,i) = m_liq_av(j,i) &
3832	/ REAL( average_count_3d, KIND=wp )
3833	ENDDO
3834	ENDDO
3835
3836	CASE ( 'm_soil' )
3837	DO i = nxl, nxr
3838	DO j = nys, nyn
3839	DO k = nzb_soil, nzt_soil
3840	m_soil_av(k,j,i) = m_soil_av(k,j,i) &
3841	/ REAL( average_count_3d, KIND=wp )
3842	ENDDO
3843	ENDDO
3844	ENDDO
3845
3846	CASE ( 'qsws_liq*' )
3847	DO i = nxl, nxr
3848	DO j = nys, nyn
3849	qsws_liq_av(j,i) = qsws_liq_av(j,i) &
3850	/ REAL( average_count_3d, KIND=wp )
3851	ENDDO
3852	ENDDO
3853
3854	CASE ( 'qsws_soil*' )
3855	DO i = nxl, nxr
3856	DO j = nys, nyn
3857	qsws_soil_av(j,i) = qsws_soil_av(j,i) &
3858	/ REAL( average_count_3d, KIND=wp )
3859	ENDDO
3860	ENDDO
3861
3862	CASE ( 'qsws_veg*' )
3863	DO i = nxl, nxr
3864	DO j = nys, nyn
3865	qsws_veg_av(j,i) = qsws_veg_av(j,i) &
3866	/ REAL( average_count_3d, KIND=wp )
3867	ENDDO
3868	ENDDO
3869
3870	CASE ( 'r_a*' )
3871	DO i = nxl, nxr
3872	DO j = nys, nyn
3873	r_a_av(j,i) = r_a_av(j,i) / REAL( average_count_3d, KIND=wp )
3874	ENDDO
3875	ENDDO
3876
3877	CASE ( 'r_s*' )
3878	DO i = nxl, nxr
3879	DO j = nys, nyn
3880	r_s_av(j,i) = r_s_av(j,i) / REAL( average_count_3d, KIND=wp )
3881	ENDDO
3882	ENDDO
3883
3884	CASE ( 't_soil' )
3885	DO i = nxl, nxr
3886	DO j = nys, nyn
3887	DO k = nzb_soil, nzt_soil
3888	t_soil_av(k,j,i) = t_soil_av(k,j,i) &
3889	/ REAL( average_count_3d, KIND=wp )
3890	ENDDO
3891	ENDDO
3892	ENDDO
3893	!
3894	!--
3895
3896	END SELECT
3897
3898	ENDIF
3899
3900	END SUBROUTINE lsm_3d_data_averaging
3901
3902
3903	!------------------------------------------------------------------------------!
3904	!
3905	! Description:
3906	! ------------
3907	!> Subroutine defining appropriate grid for netcdf variables.
3908	!> It is called out from subroutine netcdf.
3909	!------------------------------------------------------------------------------!
3910	SUBROUTINE lsm_define_netcdf_grid( var, found, grid_x, grid_y, grid_z )
3911
3912	IMPLICIT NONE
3913
3914	CHARACTER (LEN=*), INTENT(IN) :: var !<
3915	LOGICAL, INTENT(OUT) :: found !<
3916	CHARACTER (LEN=*), INTENT(OUT) :: grid_x !<
3917	CHARACTER (LEN=*), INTENT(OUT) :: grid_y !<
3918	CHARACTER (LEN=*), INTENT(OUT) :: grid_z !<
3919
3920	found = .TRUE.
3921
3922	!
3923	!-- Check for the grid
3924	SELECT CASE ( TRIM( var ) )
3925
3926	CASE ( 'm_soil', 't_soil', 'm_soil_xy', 't_soil_xy', 'm_soil_xz', &
3927	't_soil_xz', 'm_soil_yz', 't_soil_yz' )
3928	grid_x = 'x'
3929	grid_y = 'y'
3930	grid_z = 'zs'
3931
3932	CASE DEFAULT
3933	found = .FALSE.
3934	grid_x = 'none'
3935	grid_y = 'none'
3936	grid_z = 'none'
3937	END SELECT
3938
3939	END SUBROUTINE lsm_define_netcdf_grid
3940
3941	!------------------------------------------------------------------------------!
3942	!
3943	! Description:
3944	! ------------
3945	!> Subroutine defining 3D output variables
3946	!------------------------------------------------------------------------------!
3947	SUBROUTINE lsm_data_output_2d( av, variable, found, grid, mode, local_pf, &
3948	two_d, nzb_do, nzt_do )
3949
3950	USE indices
3951
3952	USE kinds
3953
3954
3955	IMPLICIT NONE
3956
3957	CHARACTER (LEN=*) :: grid !<
3958	CHARACTER (LEN=*) :: mode !<
3959	CHARACTER (LEN=*) :: variable !<
3960
3961	INTEGER(iwp) :: av !<
3962	INTEGER(iwp) :: i !< running index
3963	INTEGER(iwp) :: j !< running index
3964	INTEGER(iwp) :: k !< running index
3965	INTEGER(iwp) :: m !< running index
3966	INTEGER(iwp) :: nzb_do !<
3967	INTEGER(iwp) :: nzt_do !<
3968
3969	LOGICAL :: found !<
3970	LOGICAL :: two_d !< flag parameter that indicates 2D variables (horizontal cross sections)
3971
3972	REAL(wp), DIMENSION(nxl:nxr,nys:nyn,nzb:nzt+1) :: local_pf !<
3973
3974
3975	found = .TRUE.
3976
3977	SELECT CASE ( TRIM( variable ) )
3978	!
3979	!-- Before data is transfered to local_pf, transfer is it 2D dummy variable and exchange ghost points therein.
3980	!-- However, at this point this is only required for instantaneous arrays, time-averaged quantities are already exchanged.
3981	CASE ( 'c_liq*_xy' ) ! 2d-array
3982	IF ( av == 0 ) THEN
3983	DO m = 1, surf_lsm_h%ns
3984	i = surf_lsm_h%i(m)
3985	j = surf_lsm_h%j(m)
3986	local_pf(i,j,nzb+1) = surf_lsm_h%c_liq(m) * surf_lsm_h%c_veg(m)
3987	ENDDO
3988	ELSE
3989	DO i = nxl, nxr
3990	DO j = nys, nyn
3991	local_pf(i,j,nzb+1) = c_liq_av(j,i)
3992	ENDDO
3993	ENDDO
3994	ENDIF
3995
3996	two_d = .TRUE.
3997	grid = 'zu1'
3998
3999	CASE ( 'c_soil*_xy' ) ! 2d-array
4000	IF ( av == 0 ) THEN
4001	DO m = 1, surf_lsm_h%ns
4002	i = surf_lsm_h%i(m)
4003	j = surf_lsm_h%j(m)
4004	local_pf(i,j,nzb+1) = 1.0_wp - surf_lsm_h%c_veg(m)
4005	ENDDO
4006	ELSE
4007	DO i = nxl, nxr
4008	DO j = nys, nyn
4009	local_pf(i,j,nzb+1) = c_soil_av(j,i)
4010	ENDDO
4011	ENDDO
4012	ENDIF
4013
4014	two_d = .TRUE.
4015	grid = 'zu1'
4016
4017	CASE ( 'c_veg*_xy' ) ! 2d-array
4018	IF ( av == 0 ) THEN
4019	DO m = 1, surf_lsm_h%ns
4020	i = surf_lsm_h%i(m)
4021	j = surf_lsm_h%j(m)
4022	local_pf(i,j,nzb+1) = surf_lsm_h%c_veg(m)
4023	ENDDO
4024	ELSE
4025	DO i = nxl, nxr
4026	DO j = nys, nyn
4027	local_pf(i,j,nzb+1) = c_veg_av(j,i)
4028	ENDDO
4029	ENDDO
4030	ENDIF
4031
4032	two_d = .TRUE.
4033	grid = 'zu1'
4034
4035	CASE ( 'ghf*_xy' ) ! 2d-array
4036	IF ( av == 0 ) THEN
4037	DO m = 1, surf_lsm_h%ns
4038	i = surf_lsm_h%i(m)
4039	j = surf_lsm_h%j(m)
4040	local_pf(i,j,nzb+1) = surf_lsm_h%ghf(m)
4041	ENDDO
4042	ELSE
4043	DO i = nxl, nxr
4044	DO j = nys, nyn
4045	local_pf(i,j,nzb+1) = ghf_av(j,i)
4046	ENDDO
4047	ENDDO
4048	ENDIF
4049
4050	two_d = .TRUE.
4051	grid = 'zu1'
4052
4053	CASE ( 'lai*_xy' ) ! 2d-array
4054	IF ( av == 0 ) THEN
4055	DO m = 1, surf_lsm_h%ns
4056	i = surf_lsm_h%i(m)
4057	j = surf_lsm_h%j(m)
4058	�