Home

Context Navigation

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

Last change on this file since 2516 was 2516, checked in by suehring, 7 years ago
document changes
Property svn:keywords set to `Id`
File size: 190.8 KB

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