Home

Context Navigation

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

Last change on this file since 2270 was 2270, checked in by maronga, 7 years ago
major revisions in land surface model
Property svn:keywords set to `Id`
File size: 180.8 KB

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