Home

Context Navigation

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

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