Home

Context Navigation

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

Last change on this file since 2092 was 2032, checked in by knoop, 8 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 141.8 KB

Line
1	!> @file land_surface_model_mod.f90
2	!------------------------------------------------------------------------------!
3	! This file is part of PALM.
4	!
5	! PALM is free software: you can redistribute it and/or modify it under the
6	! terms of the GNU General Public License as published by the Free Software
7	! Foundation, either version 3 of the License, or (at your option) any later
8	! version.
9	!
10	! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11	! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12	! A PARTICULAR PURPOSE. See the GNU General Public License for more details.
13	!
14	! You should have received a copy of the GNU General Public License along with
15	! PALM. If not, see <http://www.gnu.org/licenses/>.
16	!
17	! Copyright 1997-2016 Leibniz Universitaet Hannover
18	!------------------------------------------------------------------------------!
19	!
20	! Current revisions:
21	! -----------------
22	!
23	!
24	! Former revisions:
25	! -----------------
26	! $Id: land_surface_model_mod.f90 2032 2016-10-21 15:13:51Z suehring $
27	!
28	! 2031 2016-10-21 15:11:58Z knoop
29	! renamed variable rho to rho_ocean
30	!
31	! 2000 2016-08-20 18:09:15Z knoop
32	! Forced header and separation lines into 80 columns
33	!
34	! 1978 2016-07-29 12:08:31Z maronga
35	! Bugfix: initial values of pave_surface and water_surface were not set.
36	!
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	!
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	!
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	!
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	!
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	!
58	! 1849 2016-04-08 11:33:18Z hoffmann
59	! prr moved to arrays_3d
60	!
61	! 1826 2016-04-07 12:01:39Z maronga
62	! Cleanup after modularization
63	!
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	!
68	! 1788 2016-03-10 11:01:04Z maronga
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.
79	!
80	! 1783 2016-03-06 18:36:17Z raasch
81	! netcdf variables moved to netcdf module
82	!
83	! 1757 2016-02-22 15:49:32Z maronga
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.
89	!
90	! 1709 2015-11-04 14:47:01Z maronga
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
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
101	! Bugfix: REAL constants provided with KIND-attribute in call of
102	! Replaced rif with ol
103	!
104	! 1691 2015-10-26 16:17:44Z maronga
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
109	!
110	! 1682 2015-10-07 23:56:08Z knoop
111	! Code annotations made doxygen readable
112	!
113	! 1590 2015-05-08 13:56:27Z maronga
114	! Bugfix: definition of character strings requires same length for all elements
115	!
116	! 1585 2015-04-30 07:05:52Z maronga
117	! Modifications for RRTMG. Changed tables to PARAMETER type.
118	!
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	!
123	! 1555 2015-03-04 17:44:27Z maronga
124	! Added output of r_a and r_s
125	!
126	! 1553 2015-03-03 17:33:54Z maronga
127	! Improved better treatment of roughness lengths. Added default soil temperature
128	! profile
129	!
130	! 1551 2015-03-03 14:18:16Z maronga
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
138	!
139	! intrinsic function MAX and MIN
140	!
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	!
146	! 1496 2014-12-02 17:25:50Z maronga
147	! Initial revision
148	!
149	!
150	! Description:
151	! ------------
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	!>
158	!> @todo Consider partial absorption of the net shortwave radiation by the
159	!> skin layer.
160	!> @todo Improve surface water parameterization
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.
165	!> @todo Fix crashes with radiation_scheme == 'constant'
166	!>
167	!> @note No time step criterion is required as long as the soil layers do not
168	!> become too thin.
169	!------------------------------------------------------------------------------!
170	MODULE land_surface_model_mod
171
172	USE arrays_3d, &
173	ONLY: hyp, ol, pt, pt_p, prr, q, q_p, ql, qsws, shf, ts, us, vpt, z0, &
174	z0h, z0q
175
176	USE cloud_parameters, &
177	ONLY: cp, hyrho, l_d_cp, l_d_r, l_v, pt_d_t, rho_l, r_d, r_v
178
179	USE control_parameters, &
180	ONLY: cloud_physics, dt_3d, humidity, intermediate_timestep_count, &
181	initializing_actions, intermediate_timestep_count_max, &
182	max_masks, precipitation, pt_surface, &
183	rho_surface, roughness_length, surface_pressure, &
184	timestep_scheme, tsc, z0h_factor, time_since_reference_point
185
186	USE indices, &
187	ONLY: nbgp, nxlg, nxrg, nyng, nysg, nzb, nzb_s_inner
188
189	USE kinds
190
191	USE pegrid
192
193	USE radiation_model_mod, &
194	ONLY: force_radiation_call, rad_net, rad_sw_in, rad_lw_out, &
195	rad_lw_out_change_0, unscheduled_radiation_calls
196
197	USE statistics, &
198	ONLY: hom, statistic_regions
199
200	IMPLICIT NONE
201
202	!
203	!-- LSM model constants
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)
207
208	REAL(wp), PARAMETER :: &
209	b_ch = 6.04_wp, & ! Clapp & Hornberger exponent
210	lambda_h_dry = 0.19_wp, & ! heat conductivity for dry soil
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
214	rho_c_soil = 2.19E6_wp, & ! volumetric heat capacity of soil
215	rho_c_water = 4.20E6_wp, & ! volumetric heat capacity of water
216	m_max_depth = 0.0002_wp ! Maximum capacity of the water reservoir (m)
217
218
219	!
220	!-- LSM variables
221	INTEGER(iwp) :: veg_type = 2, & !< NAMELIST veg_type_2d
222	soil_type = 3 !< NAMELIST soil_type_2d
223
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
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
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
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
247	lambda_h_sat = 0.0_wp, & !< heat conductivity for saturated soil
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
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)
259	q_s = 0.0_wp, & !< saturation specific humidity
260	residual_moisture = 9999999.9_wp, & !< NAMELIST m_res
261	rho_cp, & !< rho_surface * cp
262	rho_lv, & !< rho_ocean * l_v
263	rd_d_rv, & !< r_d / r_v
264	saturation_moisture = 9999999.9_wp, & !< NAMELIST m_sat
265	skip_time_do_lsm = 0.0_wp, & !< LSM is not called before this time
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)
269	z0h_eb = 9999999.9_wp, & !< NAMELIST z0h (lsm_par)
270	z0q_eb = 9999999.9_wp !< NAMELIST z0q (lsm_par)
271
272	REAL(wp), DIMENSION(nzb_soil:nzt_soil) :: &
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
276	root_fraction = (/9999999.9_wp, 9999999.9_wp, &
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)
280
281	REAL(wp), DIMENSION(nzb_soil:nzt_soil+1) :: &
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)
286
287	#if defined( __nopointer )
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)
293	#else
294	REAL(wp), DIMENSION(:,:), POINTER :: t_surface, &
295	t_surface_p, &
296	m_liq_eb, &
297	m_liq_eb_p
298
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
302	#endif
303
304	!
305	!-- Temporal tendencies for time stepping
306	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: tt_surface_m, & !< surface temperature tendency (K)
307	tm_liq_eb_m !< liquid water reservoir tendency (m)
308
309	!
310	!-- Energy balance variables
311	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: &
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
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)
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
343	r_a_av, & !< average of r_a
344	r_canopy, & !< canopy resistance
345	r_soil, & !< soil resistance
346	r_soil_min, & !< minimum soil resistance
347	r_s, & !< total surface resistance (combination of r_soil and r_canopy)
348	r_s_av, & !< average of r_s
349	r_canopy_min, & !< minimum canopy (stomatal) resistance
350	shf_eb, & !< surface flux of sensible heat
351	shf_eb_av !< average of shf_eb
352
353
354	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
355	lambda_h, & !< heat conductivity of soil (W/m/K)
356	lambda_w, & !< hydraulic diffusivity of soil (?)
357	gamma_w, & !< hydraulic conductivity of soil (W/m/K)
358	rho_c_total !< volumetric heat capacity of the actual soil matrix (?)
359
360	#if defined( __nopointer )
361	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: &
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)
368	#else
369	REAL(wp), DIMENSION(:,:,:), POINTER :: &
370	t_soil, t_soil_p, &
371	m_soil, m_soil_p
372
373	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: &
374	t_soil_av, t_soil_1, t_soil_2, &
375	m_soil_av, m_soil_1, m_soil_2
376	#endif
377
378
379	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
380	tt_soil_m, & !< t_soil storage array
381	tm_soil_m, & !< m_soil storage array
382	root_fr !< root fraction (sum=1)
383
384
385	!
386	!-- Predefined Land surface classes (veg_type)
387	CHARACTER(26), DIMENSION(0:20), PARAMETER :: veg_type_name = (/ &
388	'user defined ', & ! 0
389	'crops, mixed farming ', & ! 1
390	'short grass ', & ! 2
391	'evergreen needleleaf trees', & ! 3
392	'deciduous needleleaf trees', & ! 4
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
407	'interrupted forest ', & ! 19
408	'pavements/roads ' & ! 20
409	/)
410
411	!
412	!-- Soil model classes (soil_type)
413	CHARACTER(12), DIMENSION(0:7), PARAMETER :: soil_type_name = (/ &
414	'user defined', & ! 0
415	'coarse ', & ! 1
416	'medium ', & ! 2
417	'medium-fine ', & ! 3
418	'fine ', & ! 4
419	'very fine ', & ! 5
420	'organic ', & ! 6
421	'loamy (CH) ' & ! 7
422	/)
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
429	REAL(wp), DIMENSION(0:3,1:20), PARAMETER :: veg_pars = RESHAPE( (/ &
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
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 /) )
451
452	!
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 /) )
476
477	!
478	!-- Land surface parameters III lambda_surface_s, lambda_surface_u, f_sw_in
479	REAL(wp), DIMENSION(0:2,1:20), PARAMETER :: surface_pars = RESHAPE( (/ &
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
493	1.0E10_wp, 1.0E10_wp, 0.00_wp, & ! 14
494	1.0E10_wp, 1.0E10_wp, 0.00_wp, & ! 15
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
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 /) )
501
502	!
503	!-- Root distribution (sum = 1) level 1, level 2, level 3, level 4,
504	REAL(wp), DIMENSION(0:3,1:20), PARAMETER :: root_distribution = RESHAPE( (/ &
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
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 /) )
526
527	!
528	!-- Soil parameters according to the following porosity classes (soil_type)
529
530	!
531	!-- Soil parameters I alpha_vg, l_vg, n_vg, gamma_w_sat
532	REAL(wp), DIMENSION(0:3,1:7), PARAMETER :: soil_pars = RESHAPE( (/ &
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
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 /) )
541
542	!
543	!-- Soil parameters II m_sat, m_fc, m_wilt, m_res
544	REAL(wp), DIMENSION(0:3,1:7), PARAMETER :: m_soil_pars = RESHAPE( (/ &
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
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 /) )
553
554
555	SAVE
556
557
558	PRIVATE
559
560
561	!
562	!-- Public functions
563	PUBLIC lsm_check_data_output, lsm_check_data_output_pr, &
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
568	!
569	!-- Public parameters, constants and initial values
570	PUBLIC land_surface, skip_time_do_lsm
571
572	!
573	!-- Public grid variables
574	PUBLIC nzb_soil, nzs, nzt_soil, zs
575
576	!
577	!-- Public 2D output variables
578	PUBLIC ghf_eb, qsws_eb, qsws_liq_eb, qsws_soil_eb,qsws_veg_eb, r_a, r_s, &
579	shf_eb
580
581	!
582	!-- Public prognostic variables
583	PUBLIC m_soil, t_soil
584
585
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
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
614	INTERFACE lsm_energy_balance
615	MODULE PROCEDURE lsm_energy_balance
616	END INTERFACE lsm_energy_balance
617
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
634	INTERFACE lsm_soil_model
635	MODULE PROCEDURE lsm_soil_model
636	END INTERFACE lsm_soil_model
637
638	INTERFACE lsm_swap_timelevel
639	MODULE PROCEDURE lsm_swap_timelevel
640	END INTERFACE lsm_swap_timelevel
641
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
650	CONTAINS
651
652	!------------------------------------------------------------------------------!
653	! Description:
654	! ------------
655	!> Check data output for land surface model
656	!------------------------------------------------------------------------------!
657	SUBROUTINE lsm_check_data_output( var, unit, i, ilen, k )
658
659
660	USE control_parameters, &
661	ONLY: data_output, message_string
662
663	IMPLICIT NONE
664
665	CHARACTER (LEN=*) :: unit !<
666	CHARACTER (LEN=*) :: var !<
667
668	INTEGER(iwp) :: i
669	INTEGER(iwp) :: ilen
670	INTEGER(iwp) :: k
671
672	SELECT CASE ( TRIM( var ) )
673
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'
689
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
764
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'
778
779	CASE DEFAULT
780	unit = 'illegal'
781
782	END SELECT
783
784
785	END SUBROUTINE lsm_check_data_output
786
787
788	!------------------------------------------------------------------------------!
789	! Description:
790	! ------------
791	!> Check data output of profiles for land surface model
792	!------------------------------------------------------------------------------!
793	SUBROUTINE lsm_check_data_output_pr( variable, var_count, unit, dopr_unit )
794
795	USE control_parameters, &
796	ONLY: data_output_pr, message_string
797
798	USE indices
799
800	USE profil_parameter
801
802	USE statistics
803
804	IMPLICIT NONE
805
806	CHARACTER (LEN=*) :: unit !<
807	CHARACTER (LEN=*) :: variable !<
808	CHARACTER (LEN=*) :: dopr_unit !< local value of dopr_unit
809
810	INTEGER(iwp) :: user_pr_index !<
811	INTEGER(iwp) :: var_count !<
812
813	SELECT CASE ( TRIM( variable ) )
814
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:)
829	ENDIF
830	unit = dopr_unit
831	ENDIF
832
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:)
847	ENDIF
848	unit = dopr_unit
849	ENDIF
850
851
852	CASE DEFAULT
853	unit = 'illegal'
854
855	END SELECT
856
857
858	END SUBROUTINE lsm_check_data_output_pr
859
860
861	!------------------------------------------------------------------------------!
862	! Description:
863	! ------------
864	!> Check parameters routine for land surface model
865	!------------------------------------------------------------------------------!
866	SUBROUTINE lsm_check_parameters
867
868	USE control_parameters, &
869	ONLY: bc_pt_b, bc_q_b, constant_flux_layer, message_string, &
870	most_method, topography
871
872	USE radiation_model_mod, &
873	ONLY: radiation
874
875
876	IMPLICIT NONE
877
878
879	!
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
889
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
895
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
901
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 )
916	ENDIF
917
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
924
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
931
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'
936	CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
937	ENDIF
938
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
945
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
952
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
959
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
966
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
973
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 )
979	ENDIF
980
981
982	ENDIF
983
984	IF ( soil_type == 0 ) THEN
985
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
992
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
999
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
1006
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
1013
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
1020
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
1027
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 )
1033	ENDIF
1034
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 )
1040	ENDIF
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
1049
1050
1051	END SUBROUTINE lsm_check_parameters
1052
1053	!------------------------------------------------------------------------------!
1054	! Description:
1055	! ------------
1056	!> Solver for the energy balance at the surface.
1057	!------------------------------------------------------------------------------!
1058	SUBROUTINE lsm_energy_balance
1059
1060
1061	IMPLICIT NONE
1062
1063	INTEGER(iwp) :: i !< running index
1064	INTEGER(iwp) :: j !< running index
1065	INTEGER(iwp) :: k, ks !< running index
1066
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
1088
1089	!
1090	!-- Calculate the exner function for the current time step
1091	exn = ( surface_pressure / 1000.0_wp )**0.286_wp
1092
1093	DO i = nxlg, nxrg
1094	DO j = nysg, nyng
1095	k = nzb_s_inner(j,i)
1096
1097	!
1098	!-- Set lambda_surface according to stratification between skin layer and soil
1099	IF ( .NOT. pave_surface(j,i) ) THEN
1100
1101	c_surface_tmp = c_surface
1102
1103	IF ( t_surface(j,i) >= t_soil(nzb_soil,j,i)) THEN
1104	lambda_surface = lambda_surface_s(j,i)
1105	ELSE
1106	lambda_surface = lambda_surface_u(j,i)
1107	ENDIF
1108	ELSE
1109
1110	c_surface_tmp = pave_heat_capacity * dz_soil(nzb_soil) * 0.5_wp
1111	lambda_surface = pave_heat_conductivity * ddz_soil(nzb_soil)
1112
1113	ENDIF
1114
1115	!
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
1127
1128	r_a(j,i) = (pt1 - pt(k,j,i)) / (ts(j,i) * us(j,i) + 1.0E-20_wp)
1129
1130	!
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
1134
1135	!
1136	!-- Second step: calculate canopy resistance r_canopy
1137	!-- f1-f3 here are defined as 1/f1-f3 as in ECMWF documentation
1138
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)) )
1144
1145
1146
1147	!
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
1156
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
1164
1165	!
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 ) )
1169
1170	!
1171	!-- f3: correction for vapour pressure deficit
1172	IF ( g_d(j,i) /= 0.0_wp ) THEN
1173	!
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
1180
1181	!
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)
1187
1188	!
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
1197
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)
1201
1202	r_soil(j,i) = r_soil_min(j,i) / f2
1203
1204	!
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
1219
1220	!
1221	!-- Calculate saturation specific humidity
1222	q_s = 0.622_wp * e_s / surface_pressure
1223
1224	!
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
1231
1232	!
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
1251	!
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
1257
1258	f_shf = rho_cp / r_a(j,i)
1259	f_qsws = f_qsws_veg + f_qsws_soil + f_qsws_liq
1260
1261	!
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 )
1266
1267	dq_s_dt = 0.622_wp * e_s_dt / surface_pressure
1268
1269	!
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)
1272
1273	!
1274	!-- Calculate new skin temperature
1275	IF ( humidity ) THEN
1276
1277	!
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)
1284
1285	!
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
1290
1291	!
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)
1297
1298	!
1299	!-- Denominator of the prognostic equation
1300	coef_2 = rad_lw_out_change_0(j,i) + lambda_surface + f_shf / exn
1301
1302	ENDIF
1303
1304	tend = 0.0_wp
1305
1306	!
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 &
1311	* dt_3d * tsc(2) )
1312
1313	!
1314	!-- Add RK3 term
1315	IF ( c_surface_tmp /= 0.0_wp ) THEN
1316
1317	t_surface_p(j,i) = t_surface_p(j,i) + dt_3d * tsc(3) &
1318	* tt_surface_m(j,i)
1319
1320	!
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))
1324	!
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)
1333	ENDIF
1334	ENDIF
1335	ENDIF
1336
1337	!
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
1350
1351	pt(k,j,i) = t_surface_p(j,i) / exn
1352
1353	!
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)
1358
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) )
1362
1363	ghf_eb(j,i) = lambda_surface * (t_surface_p(j,i) &
1364	- t_soil(nzb_soil,j,i))
1365
1366	shf_eb(j,i) = - f_shf * ( pt1 - pt(k,j,i) )
1367
1368	shf(j,i) = shf_eb(j,i) / rho_cp
1369
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) )
1373
1374	qsws(j,i) = qsws_eb(j,i) / rho_lv
1375
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) )
1379
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) )
1383
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
1388
1389	!
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
1398
1399	!
1400	!-- Calculate change in liquid water reservoir due to dew fall or
1401	!-- evaporation of liquid water
1402	IF ( humidity ) THEN
1403	!
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
1408
1409	!
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
1422	ENDIF
1423
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
1431	!
1432	!-- If the air is saturated, check the reservoir water level
1433	IF ( qsws_eb(j,i) < 0.0_wp ) THEN
1434
1435	!
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)
1443
1444	qsws_liq_eb(j,i) = 0.0_wp
1445	ENDIF
1446
1447	!
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
1456
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) )
1460
1461	!
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)
1465
1466	!
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)
1470
1471
1472	!
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)
1481	ENDIF
1482	ENDIF
1483
1484	ENDIF
1485
1486	ENDDO
1487	ENDDO
1488
1489	!
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
1494	#if defined( __parallel )
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 )
1498	#else
1499	force_radiation_call = force_radiation_call_l
1500	#endif
1501	force_radiation_call_l = .FALSE.
1502	ENDIF
1503
1504	!
1505	!-- Calculate surface specific humidity
1506	IF ( humidity ) THEN
1507	CALL calc_q_surface
1508	ENDIF
1509
1510	!
1511	!-- Calculate new roughness lengths (for water surfaces only)
1512	CALL calc_z0_water_surface
1513
1514
1515	END SUBROUTINE lsm_energy_balance
1516
1517
1518	!------------------------------------------------------------------------------!
1519	! Description:
1520	! ------------
1521	!> Header output for land surface model
1522	!------------------------------------------------------------------------------!
1523	SUBROUTINE lsm_header ( io )
1524
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	' ------------------------------'/)
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)
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,' '/ &
1591	' Grid point: ',A)
1592
1593	END SUBROUTINE lsm_header
1594
1595
1596	!------------------------------------------------------------------------------!
1597	! Description:
1598	! ------------
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
1643	qsws_eb = rho_lv * qsws
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.
1903	t_soil(:,j,i) = t_surface(j,i)
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
2010	water_surface = .FALSE.
2011	pave_surface = .FALSE.
2012
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	! ------------
2083	!> Soil model as part of the land surface model. The model predicts soil
2084	!> temperature and water content.
2085	!------------------------------------------------------------------------------!
2086	SUBROUTINE lsm_soil_model
2087
2088
2089	IMPLICIT NONE
2090
2091	INTEGER(iwp) :: i !< running index
2092	INTEGER(iwp) :: j !< running index
2093	INTEGER(iwp) :: k !< running index
2094
2095	REAL(wp) :: h_vg !< Van Genuchten coef. h
2096
2097	REAL(wp), DIMENSION(nzb_soil:nzt_soil) :: gamma_temp, & !< temp. gamma
2098	lambda_temp, & !< temp. lambda
2099	tend !< tendency
2100
2101	DO i = nxlg, nxrg
2102	DO j = nysg, nyng
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
2119	!
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))
2124
2125	!
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)
2130
2131	ke = 1.0_wp + LOG10(MAX(0.1_wp,m_soil(k,j,i) &
2132	/ m_sat(j,i)))
2133
2134	lambda_temp(k) = ke * (lambda_h_sat - lambda_h_dry) + &
2135	lambda_h_dry
2136	ENDIF
2137
2138	ENDDO
2139
2140	!
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)
2157
2158
2159
2160
2161	!
2162	!-- Prognostic equation for soil temperature t_soil
2163	tend(:) = 0.0_wp
2164
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)
2169
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)
2179
2180	ENDDO
2181
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) )
2187
2188	!
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 &
2199	* tt_soil_m(k,j,i)
2200	ENDDO
2201	ENDIF
2202	ENDIF
2203
2204
2205	DO k = nzb_soil, nzt_soil
2206
2207	!
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 )
2213
2214	!
2215	!-- Parametrization of Van Genuchten
2216	IF ( soil_type /= 7 ) THEN
2217	!
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)
2224
2225
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) ) )
2234
2235	!
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
2241
2242	ENDDO
2243
2244	!
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
2254
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
2259
2260	ENDDO
2261
2262	!
2263	!
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.
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
2272
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.
2281
2282	!
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.
2287	m_total = 0.0_wp
2288	DO k = nzb_soil, nzt_soil
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
2293
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
2297	root_extr(k) = root_fr(k,j,i) * m_soil(k,j,i) &
2298	/ m_total
2299	ELSE
2300	root_extr(k) = 0.0_wp
2301	ENDIF
2302	ENDDO
2303	ENDIF
2304
2305	!
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) * ( &
2310	m_soil(nzb_soil+1,j,i) - m_soil(nzb_soil,j,i) ) &
2311	* ddz_soil(nzb_soil+1) - gamma_w(nzb_soil,j,i) - ( &
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)
2315
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)
2325
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)
2336
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
2341	!
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
2346
2347	!
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&
2358	* tm_soil_m(k,j,i)
2359	ENDDO
2360	ENDIF
2361	ENDIF
2362
2363	ENDIF
2364
2365	ENDIF
2366
2367	ENDDO
2368	ENDDO
2369
2370	END SUBROUTINE lsm_soil_model
2371
2372
2373	!------------------------------------------------------------------------------!
2374	! Description:
2375	! ------------
2376	!> Swapping of timelevels
2377	!------------------------------------------------------------------------------!
2378	SUBROUTINE lsm_swap_timelevel ( mod_count )
2379
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
2422
2423
2424	!------------------------------------------------------------------------------!
2425	!
2426	! Description:
2427	! ------------
2428	!> Subroutine for averaging 3D data
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	! ------------
2782	!> Subroutine defining appropriate grid for netcdf variables.
2783	!> It is called out from subroutine netcdf.
2784	!------------------------------------------------------------------------------!
2785	SUBROUTINE lsm_define_netcdf_grid( var, found, grid_x, grid_y, grid_z )
2786
2787	IMPLICIT NONE
2788
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 !<
2794
2795	found = .TRUE.
2796
2797	!
2798	!-- Check for the grid
2799	SELECT CASE ( TRIM( var ) )
2800
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'
2806
2807	CASE DEFAULT
2808	found = .FALSE.
2809	grid_x = 'none'
2810	grid_y = 'none'
2811	grid_z = 'none'
2812	END SELECT
2813
2814	END SUBROUTINE lsm_define_netcdf_grid
2815
2816	!------------------------------------------------------------------------------!
2817	!
2818	! Description:
2819	! ------------
2820	!> Subroutine defining 3D output variables
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	! ------------
3149	!> Subroutine defining 3D output variables
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	! ------------
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
3537	!
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 ) )
3548
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
3555
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 )
3563
3564	z0h(j,i) = 0.40_wp * molecular_viscosity / us(j,i)
3565	z0q(j,i) = 0.62_wp * molecular_viscosity / us(j,i)
3566
3567	ENDIF
3568	ENDDO
3569	ENDDO
3570
3571	END SUBROUTINE calc_z0_water_surface
3572
3573
3574	!------------------------------------------------------------------------------!
3575	! Description:
3576	! ------------
3577	!> Calculation of specific humidity of the skin layer (surface). It is assumend
3578	!> that the skin is always saturated.
3579	!------------------------------------------------------------------------------!
3580	SUBROUTINE calc_q_surface
3581
3582	IMPLICIT NONE
3583
3584	INTEGER :: i !< running index
3585	INTEGER :: j !< running index
3586	INTEGER :: k !< running index
3587
3588	REAL(wp) :: resistance !< aerodynamic and soil resistance term
3589
3590	DO i = nxlg, nxrg
3591	DO j = nysg, nyng
3592	k = nzb_s_inner(j,i)
3593
3594	!
3595	!-- Calculate water vapour pressure at saturation
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 ) )
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
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
3614
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
3619	ENDDO
3620	ENDDO
3621
3622	END SUBROUTINE calc_q_surface
3623
3624
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 |