Home

Context Navigation

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

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