Home

Context Navigation

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

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