Home

Context Navigation

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

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