Home

Context Navigation

source: palm/trunk/SOURCE/land_surface_model_mod.f90 @ 2514

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |