Home

Context Navigation

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

Last change on this file since 2583 was 2575, checked in by maronga, 6 years ago
bugfix in radiation model and improvements in land surface scheme
Property svn:keywords set to `Id`
File size: 194.8 KB

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