Home

Context Navigation

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

Last change on this file since 2504 was 2504, checked in by maronga, 7 years ago
lsm now allows for moist soil and roots below paved surfaces
Property svn:keywords set to `Id`
File size: 190.6 KB

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