Home

Context Navigation

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

Last change on this file since 2381 was 2354, checked in by schwenkel, 7 years ago
bugfix for lsm data output
Property svn:keywords set to `Id`
File size: 185.6 KB

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