Home

Context Navigation

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

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