Home

Context Navigation

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

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