Home

Context Navigation

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

Last change on this file since 2630 was 2608, checked in by schwenkel, 7 years ago
Inital revision of diagnostic_quantities_mod allows unified calculation of magnus equation and saturion mixing ratio
Property svn:keywords set to `Id`
File size: 195.0 KB

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