Home

Context Navigation

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

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