Home

Context Navigation

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

Last change on this file since 2332 was 2332, checked in by maronga, 7 years ago
bugfix for previous revision
Property svn:keywords set to `Id`
File size: 185.8 KB

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