Home

Context Navigation

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

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