Home

Context Navigation

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

Last change on this file since 2300 was 2299, checked in by maronga, 7 years ago
improvements for spinup mechanism
Property svn:keywords set to `Id`
File size: 181.6 KB

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