Home

Context Navigation

land_surface_model_mod.f90 @ 2475

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