Home

Context Navigation

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

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