Home

Context Navigation

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

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