Home

Context Navigation

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

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |