Home

Context Navigation

source: palm/trunk/SOURCE/plant_canopy_model_mod.f90 @ 4526

Last change on this file since 4526 was 4525, checked in by raasch, 5 years ago
salsa: added reading/writing of global restart data + reading/writing restart data with MPI-IO, bugfix for MPI-IO in plant_canopy_model_mod, tutorial user-interface dispersion_eularian_and_lpm_extended updated
Property svn:keywords set to `Id`
File size: 117.7 KB

Line
1	!> @file plant_canopy_model_mod.f90
2	!------------------------------------------------------------------------------!
3	! This file is part of the PALM model system.
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-2020 Leibniz Universitaet Hannover
18	! Copyright 2017-2020 Institute of Computer Science of the
19	! Czech Academy of Sciences, Prague
20	!------------------------------------------------------------------------------!
21	!
22	! Current revisions:
23	! ------------------
24	!
25	!
26	! Former revisions:
27	! -----------------
28	! $Id: plant_canopy_model_mod.f90 4525 2020-05-10 17:05:07Z oliver.maas $
29	! bugfix for reading/writing pcm_...rate_av with MPI-IO
30	!
31	! 4517 2020-05-03 14:29:30Z raasch
32	! added restart with MPI-IO for reading local arrays
33	!
34	! 4515 2020-04-30 16:37:18Z suehring
35	! Rename error number again since this was recently given in -r 4511
36	!
37	! 4514 2020-04-30 16:29:59Z suehring
38	! - Bugfix in output of pcm_heatrate_av in a restart run. In order to fix this,
39	! pch_index is now output for a restart run. Therefore, define global restart
40	! routines.
41	! - Error message number renamed and check for PA0505 revised in order to also
42	! consider natural surfaces with plant-canopy.
43	!
44	! 4495 2020-04-13 20:11:20Z raasch
45	! restart data handling with MPI-IO added
46	!
47	! 4457 2020-03-11 14:20:43Z raasch
48	!
49	! use statement for exchange horiz added
50	! (salim) removed the error message PA0672 to consider PC 3d data via ascii file
51	!
52	! 4392 2020-01-31 16:14:57Z pavelkrc (resler)
53	! Make pcm_heatrate_av, pcm_latentrate_av public to allow calculation
54	! of averaged Bowen ratio in the user procedure
55	!
56	! 4381 2020-01-20 13:51:46Z suehring
57	! Give error message 313 only once
58	!
59	! 4363 2020-01-07 18:11:28Z suehring
60	! Fix for last commit
61	!
62	! 4362 2020-01-07 17:15:02Z suehring
63	! Input of plant canopy variables from static driver moved to plant-canopy
64	! model
65	!
66	! 4361 2020-01-07 12:22:38Z suehring
67	! - Remove unused arrays in pmc_rrd_local
68	! - Remove one exchange of ghost points
69	!
70	! 4360 2020-01-07 11:25:50Z suehring
71	! - Bugfix, read restart data for time-averaged pcm output quantities
72	! - Output of plant-canopy quantities will fill values
73	!
74	! 4356 2019-12-20 17:09:33Z suehring
75	! Correct single message call, local check must be given by the respective
76	! mpi rank.
77	!
78	! 4346 2019-12-18 11:55:56Z motisi
79	! Introduction of wall_flags_total_0, which currently sets bits based on static
80	! topography information used in wall_flags_static_0
81	!
82	! 4342 2019-12-16 13:49:14Z Giersch
83	! Use statements moved to module level, ocean dependency removed, redundant
84	! variables removed
85	!
86	! 4341 2019-12-16 10:43:49Z motisi
87	! - Unification of variable names: pc_-variables now pcm_-variables
88	! (pc_latent_rate, pc_heating_rate, pc_transpiration_rate)
89	! - Removal of pcm_bowenratio output
90	! - Renamed canopy-mode 'block' to 'homogeneous'
91	! - Renamed value 'read_from_file_3d' to 'read_from_file'
92	! - Removal of confusing comment lines
93	! - Replacement of k_wall by topo_top_ind
94	! - Removal of Else-Statement in tendency-calculation
95	!
96	! 4335 2019-12-12 16:39:05Z suehring
97	! Fix for LAD at building edges also implemented in vector branch.
98	!
99	! 4331 2019-12-10 18:25:02Z suehring
100	! Typo corrected
101	!
102	! 4329 2019-12-10 15:46:36Z motisi
103	! Renamed wall_flags_0 to wall_flags_static_0
104	!
105	! 4314 2019-11-29 10:29:20Z suehring
106	! - Bugfix, plant canopy was still considered at building edges on for the u-
107	! and v-component.
108	! - Relax restriction of LAD on building tops. LAD is only omitted at
109	! locations where building grid points emerged artificially by the
110	! topography filtering.
111	!
112	! 4309 2019-11-26 18:49:59Z suehring
113	! Typo
114	!
115	! 4302 2019-11-22 13:15:56Z suehring
116	! Omit tall canopy mapped on top of buildings
117	!
118	! 4279 2019-10-29 08:48:17Z scharf
119	! unused variables removed
120	!
121	! 4258 2019-10-07 13:29:08Z scharf
122	! changed check for static driver and fixed bugs in initialization and header
123	!
124	! 4258 2019-10-07 13:29:08Z suehring
125	! Check if any LAD is prescribed when plant-canopy model is applied.
126	!
127	! 4226 2019-09-10 17:03:24Z suehring
128	! Bugfix, missing initialization of heating rate
129	!
130	! 4221 2019-09-09 08:50:35Z suehring
131	! Further bugfix in 3d data output for plant canopy
132	!
133	! 4216 2019-09-04 09:09:03Z suehring
134	! Bugfixes in 3d data output
135	!
136	! 4205 2019-08-30 13:25:00Z suehring
137	! Missing working precision + bugfix in calculation of wind speed
138	!
139	! 4188 2019-08-26 14:15:47Z suehring
140	! Minor adjustment in error number
141	!
142	! 4187 2019-08-26 12:43:15Z suehring
143	! Give specific error numbers instead of PA0999
144	!
145	! 4182 2019-08-22 15:20:23Z scharf
146	! Corrected "Former revisions" section
147	!
148	! 4168 2019-08-16 13:50:17Z suehring
149	! Replace function get_topography_top_index by topo_top_ind
150	!
151	! 4127 2019-07-30 14:47:10Z suehring
152	! Output of 3D plant canopy variables changed. It is now relative to the local
153	! terrain rather than located at the acutal vertical level in the model. This
154	! way, the vertical dimension of the output can be significantly reduced.
155	! (merge from branch resler)
156	!
157	! 3885 2019-04-11 11:29:34Z kanani
158	! Changes related to global restructuring of location messages and introduction
159	! of additional debug messages
160	!
161	! 3864 2019-04-05 09:01:56Z monakurppa
162	! unsed variables removed
163	!
164	! 3745 2019-02-15 18:57:56Z suehring
165	! Bugfix in transpiration, floating invalid when temperature
166	! becomes > 40 degrees
167	!
168	! 3744 2019-02-15 18:38:58Z suehring
169	! Some interface calls moved to module_interface + cleanup
170	!
171	! 3655 2019-01-07 16:51:22Z knoop
172	! unused variables removed
173	!
174	! 138 2007-11-28 10:03:58Z letzel
175	! Initial revision
176	!
177	! Description:
178	! ------------
179	!> 1) Initialization of the canopy model, e.g. construction of leaf area density
180	!> profile (subroutine pcm_init).
181	!> 2) Calculation of sinks and sources of momentum, heat and scalar concentration
182	!> due to canopy elements (subroutine pcm_tendency).
183	!
184	! @todo - precalculate constant terms in pcm_calc_transpiration_rate
185	! @todo - unify variable names (pcm_, pc_, ...)
186	! @todo - get rid-off dependency on radiation model
187	!------------------------------------------------------------------------------!
188	MODULE plant_canopy_model_mod
189
190	USE arrays_3d, &
191	ONLY: dzu, dzw, e, exner, hyp, pt, q, s, tend, u, v, w, zu, zw
192
193	USE basic_constants_and_equations_mod, &
194	ONLY: c_p, degc_to_k, l_v, lv_d_cp, r_d, rd_d_rv
195
196	USE bulk_cloud_model_mod, &
197	ONLY: bulk_cloud_model, microphysics_seifert
198
199	USE control_parameters, &
200	ONLY: average_count_3d, &
201	coupling_char, &
202	debug_output, &
203	dt_3d, &
204	dz, &
205	humidity, &
206	land_surface, &
207	length, &
208	message_string, &
209	ocean_mode, &
210	passive_scalar, &
211	plant_canopy, &
212	restart_data_format_output, &
213	restart_string, &
214	urban_surface
215
216	USE grid_variables, &
217	ONLY: dx, dy
218
219	USE indices, &
220	ONLY: nbgp, nxl, nxlg, nxlu, nxr, nxrg, nyn, nyng, nys, nysg, nysv, &
221	nz, nzb, nzt, topo_top_ind, wall_flags_total_0
222
223	USE kinds
224
225	USE netcdf_data_input_mod, &
226	ONLY: input_pids_static, &
227	char_fill, &
228	check_existence, &
229	close_input_file, &
230	get_attribute, &
231	get_dimension_length, &
232	get_variable, &
233	inquire_num_variables, &
234	inquire_variable_names, &
235	input_file_static, &
236	num_var_pids, &
237	open_read_file, &
238	pids_id, &
239	real_3d, &
240	vars_pids
241
242	USE pegrid
243
244	USE restart_data_mpi_io_mod, &
245	ONLY: rd_mpi_io_check_array, &
246	rrd_mpi_io, &
247	wrd_mpi_io
248
249	USE surface_mod, &
250	ONLY: surf_def_h, surf_lsm_h, surf_usm_h
251
252
253	IMPLICIT NONE
254
255	CHARACTER (LEN=30) :: canopy_mode = 'homogeneous' !< canopy coverage
256	LOGICAL :: plant_canopy_transpiration = .FALSE. !< flag to switch calculation of transpiration and corresponding latent heat
257	!< for resolved plant canopy inside radiation model
258	!< (calls subroutine pcm_calc_transpiration_rate from module plant_canopy_mod)
259
260	INTEGER(iwp) :: pch_index = 0 !< plant canopy height/top index
261	INTEGER(iwp) :: lad_vertical_gradient_level_ind(10) = -9999 !< lad-profile levels (index)
262
263	INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: pch_index_ji !< local plant canopy top
264
265	LOGICAL :: calc_beta_lad_profile = .FALSE. !< switch for calc. of lad from beta func.
266
267	REAL(wp) :: alpha_lad = 9999999.9_wp !< coefficient for lad calculation
268	REAL(wp) :: beta_lad = 9999999.9_wp !< coefficient for lad calculation
269	REAL(wp) :: canopy_drag_coeff = 0.0_wp !< canopy drag coefficient (parameter)
270	REAL(wp) :: cthf = 0.0_wp !< canopy top heat flux
271	REAL(wp) :: dt_plant_canopy = 0.0_wp !< timestep account. for canopy drag
272	REAL(wp) :: ext_coef = 0.6_wp !< extinction coefficient
273	REAL(wp) :: lad_surface = 0.0_wp !< lad surface value
274	REAL(wp) :: lai_beta = 0.0_wp !< leaf area index (lai) for lad calc.
275	REAL(wp) :: leaf_scalar_exch_coeff = 0.0_wp !< canopy scalar exchange coeff.
276	REAL(wp) :: leaf_surface_conc = 0.0_wp !< leaf surface concentration
277
278	REAL(wp) :: lad_vertical_gradient(10) = 0.0_wp !< lad gradient
279	REAL(wp) :: lad_vertical_gradient_level(10) = -9999999.9_wp !< lad-prof. levels (in m)
280
281	REAL(wp) :: lad_type_coef(0:10) = 1.0_wp !< multiplicative coeficients for particular types
282	!< of plant canopy (e.g. deciduous tree during winter)
283
284	REAL(wp), DIMENSION(:), ALLOCATABLE :: lad !< leaf area density
285	REAL(wp), DIMENSION(:), ALLOCATABLE :: pre_lad !< preliminary lad
286
287	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: cum_lai_hf !< cumulative lai for heatflux calc.
288	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: lad_s !< lad on scalar-grid
289	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: pcm_heating_rate !< plant canopy heating rate
290	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: pcm_transpiration_rate !< plant canopy transpiration rate
291	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: pcm_latent_rate !< plant canopy latent heating rate
292
293	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: pcm_heatrate_av !< array for averaging plant canopy sensible heating rate
294	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: pcm_latentrate_av !< array for averaging plant canopy latent heating rate
295	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: pcm_transpirationrate_av !< array for averaging plant canopy transpiration rate
296
297	TYPE(real_3d) :: basal_area_density_f !< input variable for basal area density - resolved vegetation
298	TYPE(real_3d) :: leaf_area_density_f !< input variable for leaf area density - resolved vegetation
299	TYPE(real_3d) :: root_area_density_lad_f !< input variable for root area density - resolved vegetation
300
301	SAVE
302
303	PRIVATE
304
305	!
306	!-- Public functions
307	PUBLIC pcm_calc_transpiration_rate, &
308	pcm_check_data_output, &
309	pcm_check_parameters, &
310	pcm_3d_data_averaging, &
311	pcm_data_output_3d, &
312	pcm_define_netcdf_grid, &
313	pcm_header, &
314	pcm_init, &
315	pcm_parin, &
316	pcm_rrd_global, &
317	pcm_rrd_local, &
318	pcm_tendency, &
319	pcm_wrd_global, &
320	pcm_wrd_local
321
322	!
323	!-- Public variables and constants
324	PUBLIC canopy_drag_coeff, pcm_heating_rate, pcm_transpiration_rate, &
325	pcm_latent_rate, canopy_mode, cthf, dt_plant_canopy, lad, lad_s, &
326	pch_index, plant_canopy_transpiration, &
327	pcm_heatrate_av, pcm_latentrate_av
328
329	INTERFACE pcm_calc_transpiration_rate
330	MODULE PROCEDURE pcm_calc_transpiration_rate
331	END INTERFACE pcm_calc_transpiration_rate
332
333	INTERFACE pcm_check_data_output
334	MODULE PROCEDURE pcm_check_data_output
335	END INTERFACE pcm_check_data_output
336
337	INTERFACE pcm_check_parameters
338	MODULE PROCEDURE pcm_check_parameters
339	END INTERFACE pcm_check_parameters
340
341	INTERFACE pcm_3d_data_averaging
342	MODULE PROCEDURE pcm_3d_data_averaging
343	END INTERFACE pcm_3d_data_averaging
344
345	INTERFACE pcm_data_output_3d
346	MODULE PROCEDURE pcm_data_output_3d
347	END INTERFACE pcm_data_output_3d
348
349	INTERFACE pcm_define_netcdf_grid
350	MODULE PROCEDURE pcm_define_netcdf_grid
351	END INTERFACE pcm_define_netcdf_grid
352
353	INTERFACE pcm_header
354	MODULE PROCEDURE pcm_header
355	END INTERFACE pcm_header
356
357	INTERFACE pcm_init
358	MODULE PROCEDURE pcm_init
359	END INTERFACE pcm_init
360
361	INTERFACE pcm_parin
362	MODULE PROCEDURE pcm_parin
363	END INTERFACE pcm_parin
364
365	INTERFACE pcm_read_plant_canopy_3d
366	MODULE PROCEDURE pcm_read_plant_canopy_3d
367	END INTERFACE pcm_read_plant_canopy_3d
368
369	INTERFACE pcm_rrd_local
370	MODULE PROCEDURE pcm_rrd_local_ftn
371	MODULE PROCEDURE pcm_rrd_local_mpi
372	END INTERFACE pcm_rrd_local
373
374	INTERFACE pcm_rrd_global
375	MODULE PROCEDURE pcm_rrd_global_ftn
376	MODULE PROCEDURE pcm_rrd_global_mpi
377	END INTERFACE pcm_rrd_global
378
379	INTERFACE pcm_tendency
380	MODULE PROCEDURE pcm_tendency
381	MODULE PROCEDURE pcm_tendency_ij
382	END INTERFACE pcm_tendency
383
384	INTERFACE pcm_wrd_local
385	MODULE PROCEDURE pcm_wrd_local
386	END INTERFACE pcm_wrd_local
387
388	INTERFACE pcm_wrd_global
389	MODULE PROCEDURE pcm_wrd_global
390	END INTERFACE pcm_wrd_global
391
392
393	CONTAINS
394
395
396	!------------------------------------------------------------------------------!
397	! Description:
398	! ------------
399	!> Calculation of the plant canopy transpiration rate based on the Jarvis-Stewart
400	!> with parametrizations described in Daudet et al. (1999; Agricult. and Forest
401	!> Meteorol. 97) and Ngao, Adam and Saudreau (2017; Agricult. and Forest Meteorol
402	!> 237-238). Model functions f1-f4 were adapted from Stewart (1998; Agric.
403	!> and Forest. Meteorol. 43) instead, because they are valid for broader intervals
404	!> of values. Funcion f4 used in form present in van Wijk et al. (1998;
405	!> Tree Physiology 20).
406	!>
407	!> This subroutine is called from subroutine radiation_interaction
408	!> after the calculation of radiation in plant canopy boxes.
409	!> (arrays pcbinsw and pcbinlw).
410	!>
411	!------------------------------------------------------------------------------!
412	SUBROUTINE pcm_calc_transpiration_rate(i, j, k, kk, pcbsw, pcblw, pcbtr, pcblh)
413
414	!
415	!-- input parameters
416	INTEGER(iwp), INTENT(IN) :: i, j, k, kk !< indices of the pc gridbox
417	REAL(wp), INTENT(IN) :: pcbsw !< sw radiation in gridbox (W)
418	REAL(wp), INTENT(IN) :: pcblw !< lw radiation in gridbox (W)
419	REAL(wp), INTENT(OUT) :: pcbtr !< transpiration rate dq/dt (kg/kg/s)
420	REAL(wp), INTENT(OUT) :: pcblh !< latent heat from transpiration dT/dt (K/s)
421
422	!-- variables and parameters for calculation of transpiration rate
423	REAL(wp) :: sat_press, sat_press_d, temp, v_lad
424	REAL(wp) :: d_fact, g_b, g_s, wind_speed, evapor_rate
425	REAL(wp) :: f1, f2, f3, f4, vpd, rswc, e_eq, e_imp, rad
426	REAL(wp), PARAMETER :: gama_psychr = 66.0_wp !< psychrometric constant (Pa/K)
427	REAL(wp), PARAMETER :: g_s_max = 0.01 !< maximum stomatal conductivity (m/s)
428	REAL(wp), PARAMETER :: m_soil = 0.4_wp !< soil water content (needs to adjust or take from LSM)
429	REAL(wp), PARAMETER :: m_wilt = 0.01_wp !< wilting point soil water content (needs to adjust or take from LSM)
430	REAL(wp), PARAMETER :: m_sat = 0.51_wp !< saturation soil water content (needs to adjust or take from LSM)
431	REAL(wp), PARAMETER :: t2_min = 0.0_wp !< minimal temperature for calculation of f2
432	REAL(wp), PARAMETER :: t2_max = 40.0_wp !< maximal temperature for calculation of f2
433
434
435	!-- Temperature (deg C)
436	temp = pt(k,j,i) * exner(k) - degc_to_k
437	!-- Coefficient for conversion of radiation to grid to radiation to unit leaves surface
438	v_lad = 1.0_wp / ( MAX( lad_s(kk,j,i), 1.0E-10_wp ) * dx * dy * dz(1) )
439	!-- Magnus formula for the saturation pressure (see Ngao, Adam and Saudreau (2017) eq. 1)
440	!-- There are updated formulas available, kept consistent with the rest of the parametrization
441	sat_press = 610.8_wp * exp(17.27_wp * temp/(temp + 237.3_wp))
442	!-- Saturation pressure derivative (derivative of the above)
443	sat_press_d = sat_press * 17.27_wp * 237.3_wp / (temp + 237.3_wp)**2
444	!-- Wind speed
445	wind_speed = SQRT( ( 0.5_wp * ( u(k,j,i) + u(k,j,i+1) ) )**2 + &
446	( 0.5_wp * ( v(k,j,i) + v(k,j+1,i) ) )**2 + &
447	( 0.5_wp * ( w(k,j,i) + w(k-1,j,i) ) )**2 )
448	!-- Aerodynamic conductivity (Daudet et al. (1999) eq. 14
449	g_b = 0.01_wp * wind_speed + 0.0071_wp
450	!-- Radiation flux per leaf surface unit
451	rad = pcbsw * v_lad
452	!-- First function for calculation of stomatal conductivity (radiation dependency)
453	!-- Stewart (1988; Agric. and Forest. Meteorol. 43) eq. 17
454	f1 = rad * (1000.0_wp+42.1_wp) / 1000.0_wp / (rad+42.1_wp)
455	!-- Second function for calculation of stomatal conductivity (temperature dependency)
456	!-- Stewart (1988; Agric. and Forest. Meteorol. 43) eq. 21
457	f2 = MAX(t2_min, (temp-t2_min) * MAX(0.0_wp,t2_max-temp)**((t2_max-16.9_wp)/(16.9_wp-t2_min)) / &
458	((16.9_wp-t2_min) * (t2_max-16.9_wp)**((t2_max-16.9_wp)/(16.9_wp-t2_min))) )
459	!-- Water pressure deficit
460	!-- Ngao, Adam and Saudreau (2017) eq. 6 but with water vapour partial pressure
461	vpd = max( sat_press - q(k,j,i) * hyp(k) / rd_d_rv, 0._wp )
462	!-- Third function for calculation of stomatal conductivity (water pressure deficit dependency)
463	!-- Ngao, Adam and Saudreau (2017) Table 1, limited from below according to Stewart (1988)
464	!-- The coefficients of the linear dependence should better correspond to broad-leaved trees
465	!-- than the coefficients from Stewart (1988) which correspond to conifer trees.
466	vpd = MIN(MAX(vpd,770.0_wp),3820.0_wp)
467	f3 = -2E-4_wp * vpd + 1.154_wp
468	!-- Fourth function for calculation of stomatal conductivity (soil moisture dependency)
469	!-- Residual soil water content
470	!-- van Wijk et al. (1998; Tree Physiology 20) eq. 7
471	!-- TODO - over LSM surface might be calculated from LSM parameters
472	rswc = ( m_sat - m_soil ) / ( m_sat - m_wilt )
473	!-- van Wijk et al. (1998; Tree Physiology 20) eq. 5-6 (it is a reformulation of eq. 22-23 of Stewart(1988))
474	f4 = MAX(0.0_wp, MIN(1.0_wp - 0.041_wp * EXP(3.2_wp * rswc), 1.0_wp - 0.041_wp))
475	!-- Stomatal conductivity
476	!-- Stewart (1988; Agric. and Forest. Meteorol. 43) eq. 12
477	!-- (notation according to Ngao, Adam and Saudreau (2017) and others)
478	g_s = g_s_max * f1 * f2 * f3 * f4 + 1.0E-10_wp
479	!-- Decoupling factor
480	!-- Daudet et al. (1999) eq. 6
481	d_fact = (sat_press_d / gama_psychr + 2.0_wp ) / &
482	(sat_press_d / gama_psychr + 2.0_wp + 2.0_wp * g_b / g_s )
483	!-- Equilibrium evaporation rate
484	!-- Daudet et al. (1999) eq. 4
485	e_eq = (pcbsw + pcblw) * v_lad * sat_press_d / &
486	gama_psychr /( sat_press_d / gama_psychr + 2.0_wp ) / l_v
487	!-- Imposed evaporation rate
488	!-- Daudet et al. (1999) eq. 5
489	e_imp = r_d * pt(k,j,i) * exner(k) / hyp(k) * c_p * g_s * vpd / gama_psychr / l_v
490	!-- Evaporation rate
491	!-- Daudet et al. (1999) eq. 3
492	!-- (evaporation rate is limited to non-negative values)
493	evapor_rate = MAX(d_fact * e_eq + ( 1.0_wp - d_fact ) * e_imp, 0.0_wp)
494	!-- Conversion of evaporation rate to q tendency in gridbox
495	!-- dq/dt = E * LAD * V_g / (rho_air * V_g)
496	pcbtr = evapor_rate * r_d * pt(k,j,i) * exner(k) * lad_s(kk,j,i) / hyp(k) !-- = dq/dt
497	!-- latent heat from evaporation
498	pcblh = pcbtr * lv_d_cp !-- = - dT/dt
499
500	END SUBROUTINE pcm_calc_transpiration_rate
501
502
503	!------------------------------------------------------------------------------!
504	! Description:
505	! ------------
506	!> Check data output for plant canopy model
507	!------------------------------------------------------------------------------!
508	SUBROUTINE pcm_check_data_output( var, unit )
509
510	CHARACTER (LEN=*) :: unit !<
511	CHARACTER (LEN=*) :: var !<
512
513
514	SELECT CASE ( TRIM( var ) )
515
516	CASE ( 'pcm_heatrate' )
517	!
518	!-- Output of heatrate can be only done if it is explicitely set by cthf,
519	!-- or parametrized by absorption of radiation. The latter, however, is
520	!-- only available if radiation_interactions are on. Note, these are
521	!-- enabled if land-surface or urban-surface is switched-on. Using
522	!-- radiation_interactions_on directly is not possible since it belongs
523	!-- to the radition_model, which in turn depends on the plant-canopy model,
524	!-- creating circular dependencies.
525	IF ( cthf == 0.0_wp .AND. ( .NOT. urban_surface .AND. &
526	.NOT. land_surface ) ) THEN
527	message_string = 'output of "' // TRIM( var ) // '" requi' // &
528	'res setting of parameter cthf /= 0.0'
529	CALL message( 'pcm_check_data_output', 'PA0718', 1, 2, 0, 6, 0 )
530	ENDIF
531	unit = 'K s-1'
532
533	CASE ( 'pcm_transpirationrate' )
534	unit = 'kg kg-1 s-1'
535
536	CASE ( 'pcm_latentrate' )
537	unit = 'K s-1'
538
539	CASE ( 'pcm_lad' )
540	unit = 'm2 m-3'
541
542
543	CASE DEFAULT
544	unit = 'illegal'
545
546	END SELECT
547
548
549	END SUBROUTINE pcm_check_data_output
550
551
552	!------------------------------------------------------------------------------!
553	! Description:
554	! ------------
555	!> Check parameters routine for plant canopy model
556	!------------------------------------------------------------------------------!
557	SUBROUTINE pcm_check_parameters
558
559	IF ( ocean_mode ) THEN
560	message_string = 'plant_canopy = .TRUE. is not allowed in the '// &
561	'ocean'
562	CALL message( 'pcm_check_parameters', 'PA0696', 1, 2, 0, 6, 0 )
563	ENDIF
564
565	IF ( canopy_drag_coeff == 0.0_wp ) THEN
566	message_string = 'plant_canopy = .TRUE. requires a non-zero drag '// &
567	'coefficient & given value is canopy_drag_coeff = 0.0'
568	CALL message( 'pcm_check_parameters', 'PA0041', 1, 2, 0, 6, 0 )
569	ENDIF
570
571	IF ( ( alpha_lad /= 9999999.9_wp .AND. beta_lad == 9999999.9_wp ) .OR.&
572	beta_lad /= 9999999.9_wp .AND. alpha_lad == 9999999.9_wp ) THEN
573	message_string = 'using the beta function for the construction ' // &
574	'of the leaf area density profile requires ' // &
575	'both alpha_lad and beta_lad to be /= 9999999.9'
576	CALL message( 'pcm_check_parameters', 'PA0118', 1, 2, 0, 6, 0 )
577	ENDIF
578
579	IF ( calc_beta_lad_profile .AND. lai_beta == 0.0_wp ) THEN
580	message_string = 'using the beta function for the construction ' // &
581	'of the leaf area density profile requires ' // &
582	'a non-zero lai_beta, but given value is ' // &
583	'lai_beta = 0.0'
584	CALL message( 'pcm_check_parameters', 'PA0119', 1, 2, 0, 6, 0 )
585	ENDIF
586
587	IF ( calc_beta_lad_profile .AND. lad_surface /= 0.0_wp ) THEN
588	message_string = 'simultaneous setting of alpha_lad /= 9999999.9 '// &
589	'combined with beta_lad /= 9999999.9 ' // &
590	'and lad_surface /= 0.0 is not possible, ' // &
591	'use either vertical gradients or the beta ' // &
592	'function for the construction of the leaf area '// &
593	'density profile'
594	CALL message( 'pcm_check_parameters', 'PA0120', 1, 2, 0, 6, 0 )
595	ENDIF
596
597	IF ( bulk_cloud_model .AND. microphysics_seifert ) THEN
598	message_string = 'plant_canopy = .TRUE. requires cloud_scheme /=' // &
599	' seifert_beheng'
600	CALL message( 'pcm_check_parameters', 'PA0360', 1, 2, 0, 6, 0 )
601	ENDIF
602
603	END SUBROUTINE pcm_check_parameters
604
605
606	!------------------------------------------------------------------------------!
607	!
608	! Description:
609	! ------------
610	!> Subroutine for averaging 3D data
611	!------------------------------------------------------------------------------!
612	SUBROUTINE pcm_3d_data_averaging( mode, variable )
613
614	CHARACTER (LEN=*) :: mode !<
615	CHARACTER (LEN=*) :: variable !<
616
617	INTEGER(iwp) :: i !<
618	INTEGER(iwp) :: j !<
619	INTEGER(iwp) :: k !<
620
621
622	IF ( mode == 'allocate' ) THEN
623
624	SELECT CASE ( TRIM( variable ) )
625
626	CASE ( 'pcm_heatrate' )
627	IF ( .NOT. ALLOCATED( pcm_heatrate_av ) ) THEN
628	ALLOCATE( pcm_heatrate_av(0:pch_index,nysg:nyng,nxlg:nxrg) )
629	ENDIF
630	pcm_heatrate_av = 0.0_wp
631
632
633	CASE ( 'pcm_latentrate' )
634	IF ( .NOT. ALLOCATED( pcm_latentrate_av ) ) THEN
635	ALLOCATE( pcm_latentrate_av(0:pch_index,nysg:nyng,nxlg:nxrg) )
636	ENDIF
637	pcm_latentrate_av = 0.0_wp
638
639
640	CASE ( 'pcm_transpirationrate' )
641	IF ( .NOT. ALLOCATED( pcm_transpirationrate_av ) ) THEN
642	ALLOCATE( pcm_transpirationrate_av(0:pch_index,nysg:nyng,nxlg:nxrg) )
643	ENDIF
644	pcm_transpirationrate_av = 0.0_wp
645
646	CASE DEFAULT
647	CONTINUE
648
649	END SELECT
650
651	ELSEIF ( mode == 'sum' ) THEN
652
653	SELECT CASE ( TRIM( variable ) )
654
655	CASE ( 'pcm_heatrate' )
656	IF ( ALLOCATED( pcm_heatrate_av ) ) THEN
657	DO i = nxl, nxr
658	DO j = nys, nyn
659	IF ( pch_index_ji(j,i) /= 0 ) THEN
660	DO k = 0, pch_index_ji(j,i)
661	pcm_heatrate_av(k,j,i) = pcm_heatrate_av(k,j,i) + pcm_heating_rate(k,j,i)
662	ENDDO
663	ENDIF
664	ENDDO
665	ENDDO
666	ENDIF
667
668
669	CASE ( 'pcm_latentrate' )
670	IF ( ALLOCATED( pcm_latentrate_av ) ) THEN
671	DO i = nxl, nxr
672	DO j = nys, nyn
673	IF ( pch_index_ji(j,i) /= 0 ) THEN
674	DO k = 0, pch_index_ji(j,i)
675	pcm_latentrate_av(k,j,i) = pcm_latentrate_av(k,j,i) + pcm_latent_rate(k,j,i)
676	ENDDO
677	ENDIF
678	ENDDO
679	ENDDO
680	ENDIF
681
682
683	CASE ( 'pcm_transpirationrate' )
684	IF ( ALLOCATED( pcm_transpirationrate_av ) ) THEN
685	DO i = nxl, nxr
686	DO j = nys, nyn
687	IF ( pch_index_ji(j,i) /= 0 ) THEN
688	DO k = 0, pch_index_ji(j,i)
689	pcm_transpirationrate_av(k,j,i) = pcm_transpirationrate_av(k,j,i) + pcm_transpiration_rate(k,j,i)
690	ENDDO
691	ENDIF
692	ENDDO
693	ENDDO
694	ENDIF
695
696	CASE DEFAULT
697	CONTINUE
698
699	END SELECT
700
701	ELSEIF ( mode == 'average' ) THEN
702
703	SELECT CASE ( TRIM( variable ) )
704
705	CASE ( 'pcm_heatrate' )
706	IF ( ALLOCATED( pcm_heatrate_av ) ) THEN
707	DO i = nxlg, nxrg
708	DO j = nysg, nyng
709	IF ( pch_index_ji(j,i) /= 0 ) THEN
710	DO k = 0, pch_index_ji(j,i)
711	pcm_heatrate_av(k,j,i) = pcm_heatrate_av(k,j,i) &
712	/ REAL( average_count_3d, KIND=wp )
713	ENDDO
714	ENDIF
715	ENDDO
716	ENDDO
717	ENDIF
718
719
720	CASE ( 'pcm_latentrate' )
721	IF ( ALLOCATED( pcm_latentrate_av ) ) THEN
722	DO i = nxlg, nxrg
723	DO j = nysg, nyng
724	IF ( pch_index_ji(j,i) /= 0 ) THEN
725	DO k = 0, pch_index_ji(j,i)
726	pcm_latentrate_av(k,j,i) = pcm_latentrate_av(k,j,i) &
727	/ REAL( average_count_3d, KIND=wp )
728	ENDDO
729	ENDIF
730	ENDDO
731	ENDDO
732	ENDIF
733
734
735	CASE ( 'pcm_transpirationrate' )
736	IF ( ALLOCATED( pcm_transpirationrate_av ) ) THEN
737	DO i = nxlg, nxrg
738	DO j = nysg, nyng
739	IF ( pch_index_ji(j,i) /= 0 ) THEN
740	DO k = 0, pch_index_ji(j,i)
741	pcm_transpirationrate_av(k,j,i) = pcm_transpirationrate_av(k,j,i) &
742	/ REAL( average_count_3d, KIND=wp )
743	ENDDO
744	ENDIF
745	ENDDO
746	ENDDO
747	ENDIF
748
749	END SELECT
750
751	ENDIF
752
753	END SUBROUTINE pcm_3d_data_averaging
754
755	!------------------------------------------------------------------------------!
756	!
757	! Description:
758	! ------------
759	!> Subroutine defining 3D output variables.
760	!> Note, 3D plant-canopy output has it's own vertical output dimension, meaning
761	!> that 3D output is relative to the model surface now rather than at the actual
762	!> grid point where the plant canopy is located.
763	!------------------------------------------------------------------------------!
764	SUBROUTINE pcm_data_output_3d( av, variable, found, local_pf, fill_value, &
765	nzb_do, nzt_do )
766
767	CHARACTER (LEN=*) :: variable !< treated variable
768
769	INTEGER(iwp) :: av !< flag indicating instantaneous or averaged data output
770	INTEGER(iwp) :: i !< grid index x-direction
771	INTEGER(iwp) :: j !< grid index y-direction
772	INTEGER(iwp) :: k !< grid index z-direction
773	INTEGER(iwp) :: nzb_do !< lower limit of the data output (usually 0)
774	INTEGER(iwp) :: nzt_do !< vertical upper limit of the data output (usually nz_do3d)
775
776	LOGICAL :: found !< flag indicating if variable is found
777
778	REAL(wp) :: fill_value !< fill value
779	REAL(sp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) :: local_pf !< data output array
780
781
782	found = .TRUE.
783
784	local_pf = REAL( fill_value, KIND = 4 )
785
786	SELECT CASE ( TRIM( variable ) )
787	!
788	!-- Note, to save memory arrays for heating are allocated from 0:pch_index.
789	!-- Thus, output must be relative to these array indices. Further, check
790	!-- whether the output is within the vertical output range,
791	!-- i.e. nzb_do:nzt_do, which is necessary as local_pf is only allocated
792	!-- for this index space. Note, plant-canopy output has a separate
793	!-- vertical output coordinate zlad, so that output is mapped down to the
794	!-- surface.
795	CASE ( 'pcm_heatrate' )
796	IF ( av == 0 ) THEN
797	DO i = nxl, nxr
798	DO j = nys, nyn
799	DO k = MAX( 1, nzb_do ), MIN( pch_index_ji(j,i), nzt_do )
800	local_pf(i,j,k) = pcm_heating_rate(k,j,i)
801	ENDDO
802	ENDDO
803	ENDDO
804	ELSE
805	DO i = nxl, nxr
806	DO j = nys, nyn
807	DO k = MAX( 1, nzb_do ), MIN( pch_index_ji(j,i), nzt_do )
808	local_pf(i,j,k) = pcm_heatrate_av(k,j,i)
809	ENDDO
810	ENDDO
811	ENDDO
812	ENDIF
813
814	CASE ( 'pcm_latentrate' )
815	IF ( av == 0 ) THEN
816	DO i = nxl, nxr
817	DO j = nys, nyn
818	DO k = MAX( 1, nzb_do ), MIN( pch_index_ji(j,i), nzt_do )
819	local_pf(i,j,k) = pcm_latent_rate(k,j,i)
820	ENDDO
821	ENDDO
822	ENDDO
823	ELSE
824	DO i = nxl, nxr
825	DO j = nys, nyn
826	DO k = MAX( 1, nzb_do ), MIN( pch_index_ji(j,i), nzt_do )
827	local_pf(i,j,k) = pcm_latentrate_av(k,j,i)
828	ENDDO
829	ENDDO
830	ENDDO
831	ENDIF
832
833	CASE ( 'pcm_transpirationrate' )
834	IF ( av == 0 ) THEN
835	DO i = nxl, nxr
836	DO j = nys, nyn
837	DO k = MAX( 1, nzb_do ), MIN( pch_index_ji(j,i), nzt_do )
838	local_pf(i,j,k) = pcm_transpiration_rate(k,j,i)
839	ENDDO
840	ENDDO
841	ENDDO
842	ELSE
843	DO i = nxl, nxr
844	DO j = nys, nyn
845	DO k = MAX( 1, nzb_do ), MIN( pch_index_ji(j,i), nzt_do )
846	local_pf(i,j,k) = pcm_transpirationrate_av(k,j,i)
847	ENDDO
848	ENDDO
849	ENDDO
850	ENDIF
851
852	CASE ( 'pcm_lad' )
853	IF ( av == 0 ) THEN
854	DO i = nxl, nxr
855	DO j = nys, nyn
856	DO k = MAX( 1, nzb_do ), MIN( pch_index_ji(j,i), nzt_do )
857	local_pf(i,j,k) = lad_s(k,j,i)
858	ENDDO
859	ENDDO
860	ENDDO
861	ENDIF
862
863	CASE DEFAULT
864	found = .FALSE.
865
866	END SELECT
867
868	END SUBROUTINE pcm_data_output_3d
869
870	!------------------------------------------------------------------------------!
871	!
872	! Description:
873	! ------------
874	!> Subroutine defining appropriate grid for netcdf variables.
875	!> It is called from subroutine netcdf.
876	!------------------------------------------------------------------------------!
877	SUBROUTINE pcm_define_netcdf_grid( var, found, grid_x, grid_y, grid_z )
878
879	CHARACTER (LEN=*), INTENT(IN) :: var !<
880	LOGICAL, INTENT(OUT) :: found !<
881	CHARACTER (LEN=*), INTENT(OUT) :: grid_x !<
882	CHARACTER (LEN=*), INTENT(OUT) :: grid_y !<
883	CHARACTER (LEN=*), INTENT(OUT) :: grid_z !<
884
885	found = .TRUE.
886
887	!
888	!-- Check for the grid. zpc is zu(nzb:nzb+pch_index)
889	SELECT CASE ( TRIM( var ) )
890
891	CASE ( 'pcm_heatrate', 'pcm_lad', 'pcm_transpirationrate', 'pcm_latentrate')
892	grid_x = 'x'
893	grid_y = 'y'
894	grid_z = 'zpc'
895
896	CASE DEFAULT
897	found = .FALSE.
898	grid_x = 'none'
899	grid_y = 'none'
900	grid_z = 'none'
901	END SELECT
902
903	END SUBROUTINE pcm_define_netcdf_grid
904
905
906	!------------------------------------------------------------------------------!
907	! Description:
908	! ------------
909	!> Header output for plant canopy model
910	!------------------------------------------------------------------------------!
911	SUBROUTINE pcm_header ( io )
912
913	CHARACTER (LEN=10) :: coor_chr !<
914
915	CHARACTER (LEN=86) :: coordinates !<
916	CHARACTER (LEN=86) :: gradients !<
917	CHARACTER (LEN=86) :: leaf_area_density !<
918	CHARACTER (LEN=86) :: slices !<
919
920	INTEGER(iwp) :: i !<
921	INTEGER(iwp), INTENT(IN) :: io !< Unit of the output file
922	INTEGER(iwp) :: k !<
923
924	REAL(wp) :: canopy_height !< canopy height (in m)
925
926	canopy_height = zw(pch_index)
927
928	WRITE ( io, 1 ) canopy_mode, canopy_height, pch_index, &
929	canopy_drag_coeff
930	IF ( passive_scalar ) THEN
931	WRITE ( io, 2 ) leaf_scalar_exch_coeff, &
932	leaf_surface_conc
933	ENDIF
934
935	!
936	! Heat flux at the top of vegetation
937	WRITE ( io, 3 ) cthf
938
939	!
940	! Leaf area density profile, calculated either from given vertical
941	! gradients or from beta probability density function.
942	IF ( .NOT. calc_beta_lad_profile ) THEN
943
944	! Building output strings, starting with surface value
945	WRITE ( leaf_area_density, '(F7.4)' ) lad_surface
946	gradients = '------'
947	slices = ' 0'
948	coordinates = ' 0.0'
949	DO i = 1, UBOUND(lad_vertical_gradient_level_ind, DIM=1)
950	IF ( lad_vertical_gradient_level_ind(i) /= -9999 ) THEN
951
952	WRITE (coor_chr,'(F7.2)') lad(lad_vertical_gradient_level_ind(i))
953	leaf_area_density = TRIM( leaf_area_density ) // ' ' // TRIM( coor_chr )
954
955	WRITE (coor_chr,'(F7.2)') lad_vertical_gradient(i)
956	gradients = TRIM( gradients ) // ' ' // TRIM( coor_chr )
957
958	WRITE (coor_chr,'(I7)') lad_vertical_gradient_level_ind(i)
959	slices = TRIM( slices ) // ' ' // TRIM( coor_chr )
960
961	WRITE (coor_chr,'(F7.1)') lad_vertical_gradient_level(i)
962	coordinates = TRIM( coordinates ) // ' ' // TRIM( coor_chr )
963	ELSE
964	EXIT
965	ENDIF
966	ENDDO
967
968	WRITE ( io, 4 ) TRIM( coordinates ), TRIM( leaf_area_density ), &
969	TRIM( gradients ), TRIM( slices )
970
971	ELSE
972
973	WRITE ( leaf_area_density, '(F7.4)' ) lad_surface
974	coordinates = ' 0.0'
975
976	DO k = 1, pch_index
977
978	WRITE (coor_chr,'(F7.2)') lad(k)
979	leaf_area_density = TRIM( leaf_area_density ) // ' ' // &
980	TRIM( coor_chr )
981
982	WRITE (coor_chr,'(F7.1)') zu(k)
983	coordinates = TRIM( coordinates ) // ' ' // TRIM( coor_chr )
984
985	ENDDO
986
987	WRITE ( io, 5 ) TRIM( coordinates ), TRIM( leaf_area_density ), &
988	alpha_lad, beta_lad, lai_beta
989
990	ENDIF
991
992	1 FORMAT (//' Vegetation canopy (drag) model:'/ &
993	' ------------------------------'// &
994	' Canopy mode: ', A / &
995	' Canopy height: ',F6.2,'m (',I4,' grid points)' / &
996	' Leaf drag coefficient: ',F6.2 /)
997	2 FORMAT (/ ' Scalar exchange coefficient: ',F6.2 / &
998	' Scalar concentration at leaf surfaces in kg/m**3: ',F6.2 /)
999	3 FORMAT (' Predefined constant heatflux at the top of the vegetation: ', &
1000	F6.2, ' K m/s')
1001	4 FORMAT (/ ' Characteristic levels of the leaf area density:'// &
1002	' Height: ',A,' m'/ &
1003	' Leaf area density: ',A,' m2/m3'/ &
1004	' Gradient: ',A,' m2/m4'/ &
1005	' Gridpoint: ',A)
1006	5 FORMAT (//' Characteristic levels of the leaf area density and coefficients:'&
1007	// ' Height: ',A,' m'/ &
1008	' Leaf area density: ',A,' m2/m3'/ &
1009	' Coefficient alpha: ',F6.2 / &
1010	' Coefficient beta: ',F6.2 / &
1011	' Leaf area index: ',F6.2,' m2/m2' /)
1012
1013	END SUBROUTINE pcm_header
1014
1015
1016	!------------------------------------------------------------------------------!
1017	! Description:
1018	! ------------
1019	!> Initialization of the plant canopy model
1020	!------------------------------------------------------------------------------!
1021	SUBROUTINE pcm_init
1022
1023	USE exchange_horiz_mod, &
1024	ONLY: exchange_horiz
1025
1026	INTEGER(iwp) :: i !< running index
1027	INTEGER(iwp) :: j !< running index
1028	INTEGER(iwp) :: k !< running index
1029	INTEGER(iwp) :: m !< running index
1030
1031	LOGICAL :: lad_on_top = .FALSE. !< dummy flag to indicate that LAD is defined on a building roof
1032
1033	REAL(wp) :: canopy_height !< canopy height for lad-profile construction
1034	REAL(wp) :: gradient !< gradient for lad-profile construction
1035	REAL(wp) :: int_bpdf !< vertical integral for lad-profile construction
1036	REAL(wp) :: lad_max !< maximum LAD value in the model domain, used to perform a check
1037
1038	IF ( debug_output ) CALL debug_message( 'pcm_init', 'start' )
1039	!
1040	!-- Allocate one-dimensional arrays for the computation of the
1041	!-- leaf area density (lad) profile
1042	ALLOCATE( lad(0:nz+1), pre_lad(0:nz+1) )
1043	lad = 0.0_wp
1044	pre_lad = 0.0_wp
1045
1046	!
1047	!-- Set flag that indicates that the lad-profile shall be calculated by using
1048	!-- a beta probability density function
1049	IF ( alpha_lad /= 9999999.9_wp .AND. beta_lad /= 9999999.9_wp ) THEN
1050	calc_beta_lad_profile = .TRUE.
1051	ENDIF
1052
1053
1054	!
1055	!-- Compute the profile of leaf area density used in the plant
1056	!-- canopy model. The profile can either be constructed from
1057	!-- prescribed vertical gradients of the leaf area density or by
1058	!-- using a beta probability density function (see e.g. Markkanen et al.,
1059	!-- 2003: Boundary-Layer Meteorology, 106, 437-459)
1060	IF ( .NOT. calc_beta_lad_profile ) THEN
1061
1062	!
1063	!-- Use vertical gradients for lad-profile construction
1064	i = 1
1065	gradient = 0.0_wp
1066
1067	lad(0) = lad_surface
1068	lad_vertical_gradient_level_ind(1) = 0
1069
1070	DO k = 1, pch_index
1071	IF ( i < 11 ) THEN
1072	IF ( lad_vertical_gradient_level(i) < zu(k) .AND. &
1073	lad_vertical_gradient_level(i) >= 0.0_wp ) THEN
1074	gradient = lad_vertical_gradient(i)
1075	lad_vertical_gradient_level_ind(i) = k - 1
1076	i = i + 1
1077	ENDIF
1078	ENDIF
1079	IF ( gradient /= 0.0_wp ) THEN
1080	IF ( k /= 1 ) THEN
1081	lad(k) = lad(k-1) + dzu(k) * gradient
1082	ELSE
1083	lad(k) = lad_surface + dzu(k) * gradient
1084	ENDIF
1085	ELSE
1086	lad(k) = lad(k-1)
1087	ENDIF
1088	ENDDO
1089
1090	!
1091	!-- In case of no given leaf area density gradients, choose a vanishing
1092	!-- gradient. This information is used for the HEADER and the RUN_CONTROL
1093	!-- file.
1094	IF ( lad_vertical_gradient_level(1) == -9999999.9_wp ) THEN
1095	lad_vertical_gradient_level(1) = 0.0_wp
1096	ENDIF
1097
1098	ELSE
1099
1100	!
1101	!-- Use beta function for lad-profile construction
1102	int_bpdf = 0.0_wp
1103	canopy_height = zw(pch_index)
1104
1105	DO k = 0, pch_index
1106	int_bpdf = int_bpdf + &
1107	( ( ( zw(k) / canopy_height )*( alpha_lad-1.0_wp ) ) &
1108	( ( 1.0_wp - ( zw(k) / canopy_height ) )**( &
1109	beta_lad-1.0_wp ) ) &
1110	* ( ( zw(k+1)-zw(k) ) / canopy_height ) )
1111	ENDDO
1112
1113	!
1114	!-- Preliminary lad profile (defined on w-grid)
1115	DO k = 0, pch_index
1116	pre_lad(k) = lai_beta * &
1117	( ( ( zw(k) / canopy_height )**( alpha_lad-1.0_wp ) ) &
1118	* ( ( 1.0_wp - ( zw(k) / canopy_height ) )**( &
1119	beta_lad-1.0_wp ) ) / int_bpdf &
1120	) / canopy_height
1121	ENDDO
1122
1123	!
1124	!-- Final lad profile (defined on scalar-grid level, since most prognostic
1125	!-- quantities are defined there, hence, less interpolation is required
1126	!-- when calculating the canopy tendencies)
1127	lad(0) = pre_lad(0)
1128	DO k = 1, pch_index
1129	lad(k) = 0.5 * ( pre_lad(k-1) + pre_lad(k) )
1130	ENDDO
1131
1132	ENDIF
1133
1134	!
1135	!-- Allocate 3D-array for the leaf area density (lad_s).
1136	ALLOCATE( lad_s(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1137
1138	!
1139	!-- Initialization of the canopy coverage in the model domain:
1140	!-- Setting the parameter canopy_mode = 'homogeneous' initializes a canopy, which
1141	!-- fully covers the domain surface
1142	SELECT CASE ( TRIM( canopy_mode ) )
1143
1144	CASE( 'homogeneous' )
1145
1146	DO i = nxlg, nxrg
1147	DO j = nysg, nyng
1148	lad_s(:,j,i) = lad(:)
1149	ENDDO
1150	ENDDO
1151
1152	CASE ( 'read_from_file' )
1153	!
1154	!-- Read plant canopy
1155	IF ( input_pids_static ) THEN
1156	!
1157	!-- Open the static input file
1158	#if defined( __netcdf )
1159	CALL open_read_file( TRIM( input_file_static ) // &
1160	TRIM( coupling_char ), &
1161	pids_id )
1162
1163	CALL inquire_num_variables( pids_id, num_var_pids )
1164	!
1165	!-- Allocate memory to store variable names and read them
1166	ALLOCATE( vars_pids(1:num_var_pids) )
1167	CALL inquire_variable_names( pids_id, vars_pids )
1168	!
1169	!-- Read leaf area density - resolved vegetation
1170	IF ( check_existence( vars_pids, 'lad' ) ) THEN
1171	leaf_area_density_f%from_file = .TRUE.
1172	CALL get_attribute( pids_id, char_fill, &
1173	leaf_area_density_f%fill, &
1174	.FALSE., 'lad' )
1175	!
1176	!-- Inquire number of vertical vegetation layer
1177	CALL get_dimension_length( pids_id, &
1178	leaf_area_density_f%nz, &
1179	'zlad' )
1180	!
1181	!-- Allocate variable for leaf-area density
1182	ALLOCATE( leaf_area_density_f%var &
1183	(0:leaf_area_density_f%nz-1, &
1184	nys:nyn,nxl:nxr) )
1185
1186	CALL get_variable( pids_id, 'lad', leaf_area_density_f%var, &
1187	nxl, nxr, nys, nyn, &
1188	0, leaf_area_density_f%nz-1 )
1189
1190	ELSE
1191	leaf_area_density_f%from_file = .FALSE.
1192	ENDIF
1193	!
1194	!-- Read basal area density - resolved vegetation
1195	IF ( check_existence( vars_pids, 'bad' ) ) THEN
1196	basal_area_density_f%from_file = .TRUE.
1197	CALL get_attribute( pids_id, char_fill, &
1198	basal_area_density_f%fill, &
1199	.FALSE., 'bad' )
1200	!
1201	!-- Inquire number of vertical vegetation layer
1202	CALL get_dimension_length( pids_id, &
1203	basal_area_density_f%nz, &
1204	'zlad' )
1205	!
1206	!-- Allocate variable
1207	ALLOCATE( basal_area_density_f%var &
1208	(0:basal_area_density_f%nz-1, &
1209	nys:nyn,nxl:nxr) )
1210
1211	CALL get_variable( pids_id, 'bad', basal_area_density_f%var,&
1212	nxl, nxr, nys, nyn, &
1213	0, basal_area_density_f%nz-1 )
1214	ELSE
1215	basal_area_density_f%from_file = .FALSE.
1216	ENDIF
1217	!
1218	!-- Read root area density - resolved vegetation
1219	IF ( check_existence( vars_pids, 'root_area_dens_r' ) ) THEN
1220	root_area_density_lad_f%from_file = .TRUE.
1221	CALL get_attribute( pids_id, char_fill, &
1222	root_area_density_lad_f%fill, &
1223	.FALSE., 'root_area_dens_r' )
1224	!
1225	!-- Inquire number of vertical soil layers
1226	CALL get_dimension_length( pids_id, &
1227	root_area_density_lad_f%nz, &
1228	'zsoil' )
1229	!
1230	!-- Allocate variable
1231	ALLOCATE( root_area_density_lad_f%var &
1232	(0:root_area_density_lad_f%nz-1,&
1233	nys:nyn,nxl:nxr) )
1234
1235	CALL get_variable( pids_id, 'root_area_dens_r', &
1236	root_area_density_lad_f%var, &
1237	nxl, nxr, nys, nyn, &
1238	0, root_area_density_lad_f%nz-1 )
1239	ELSE
1240	root_area_density_lad_f%from_file = .FALSE.
1241	ENDIF
1242
1243	DEALLOCATE( vars_pids )
1244	!
1245	!-- Finally, close the input file and deallocate temporary array
1246	CALL close_input_file( pids_id )
1247	#endif
1248	ENDIF
1249
1250	!
1251	!-- Initialize LAD with data from file. If LAD is given in NetCDF file,
1252	!-- use these values, else take LAD profiles from ASCII file.
1253	!-- Please note, in NetCDF file LAD is only given up to the maximum
1254	!-- canopy top, indicated by leaf_area_density_f%nz.
1255	lad_s = 0.0_wp
1256	IF ( leaf_area_density_f%from_file ) THEN
1257	!
1258	!-- Set also pch_index, used to be the upper bound of the vertical
1259	!-- loops. Therefore, use the global top of the canopy layer.
1260	pch_index = leaf_area_density_f%nz - 1
1261
1262	DO i = nxl, nxr
1263	DO j = nys, nyn
1264	DO k = 0, leaf_area_density_f%nz - 1
1265	IF ( leaf_area_density_f%var(k,j,i) /= &
1266	leaf_area_density_f%fill ) &
1267	lad_s(k,j,i) = leaf_area_density_f%var(k,j,i)
1268	ENDDO
1269	!
1270	!-- Check if resolved vegetation is mapped onto buildings.
1271	!-- In general, this is allowed and also meaningful, e.g.
1272	!-- when trees carry across roofs. However, due to the
1273	!-- topography filtering, new building grid points can emerge
1274	!-- at location where also plant canopy is defined. As a
1275	!-- result, plant canopy is mapped on top of roofs, with
1276	!-- siginficant impact on the downstream flow field and the
1277	!-- nearby surface radiation. In order to avoid that
1278	!-- plant canopy is mistakenly mapped onto building roofs,
1279	!-- check for building grid points (bit 6) that emerge from
1280	!-- the filtering (bit 4) and set LAD to zero at these
1281	!-- artificially created building grid points. This case,
1282	!-- an informative message is given.
1283	IF ( ANY( lad_s(:,j,i) /= 0.0_wp ) .AND. &
1284	ANY( BTEST( wall_flags_total_0(:,j,i), 6 ) ) .AND. &
1285	ANY( BTEST( wall_flags_total_0(:,j,i), 4 ) ) ) THEN
1286	lad_s(:,j,i) = 0.0_wp
1287	lad_on_top = .TRUE.
1288	ENDIF
1289	ENDDO
1290	ENDDO
1291	#if defined( __parallel )
1292	CALL MPI_ALLREDUCE( MPI_IN_PLACE, lad_on_top, 1, MPI_LOGICAL, &
1293	MPI_LOR, comm2d, ierr)
1294	#endif
1295	IF ( lad_on_top ) THEN
1296	WRITE( message_string, * ) &
1297	'Resolved plant-canopy is ' // &
1298	'defined on top of an artificially '// &
1299	'created building grid point(s) ' // &
1300	'(emerged from the filtering) - ' // &
1301	'LAD profile is omitted at this / ' // &
1302	'these grid point(s).'
1303	CALL message( 'pcm_init', 'PA0313', 0, 0, 0, 6, 0 )
1304	ENDIF
1305	CALL exchange_horiz( lad_s, nbgp )
1306	!
1307	! ASCII file
1308	!-- Initialize canopy parameters canopy_drag_coeff,
1309	!-- leaf_scalar_exch_coeff, leaf_surface_conc
1310	!-- from file which contains complete 3D data (separate vertical profiles for
1311	!-- each location).
1312	ELSE
1313	CALL pcm_read_plant_canopy_3d
1314	ENDIF
1315
1316	CASE DEFAULT
1317	!
1318	!-- The DEFAULT case is reached either if the parameter
1319	!-- canopy mode contains a wrong character string or if the
1320	!-- user has coded a special case in the user interface.
1321	!-- There, the subroutine user_init_plant_canopy checks
1322	!-- which of these two conditions applies.
1323	CALL user_init_plant_canopy
1324
1325	END SELECT
1326	!
1327	!-- Check that at least one grid point has an LAD /= 0, else this may
1328	!-- cause errors in the radiation model.
1329	lad_max = MAXVAL( lad_s )
1330	#if defined( __parallel )
1331	CALL MPI_ALLREDUCE( MPI_IN_PLACE, lad_max, 1, MPI_REAL, MPI_MAX, &
1332	comm2d, ierr)
1333	#endif
1334	IF ( lad_max <= 0.0_wp ) THEN
1335	message_string = 'Plant-canopy model is switched-on but no ' // &
1336	'plant canopy is present in the model domain.'
1337	CALL message( 'pcm_init', 'PA0685', 1, 2, 0, 6, 0 )
1338	ENDIF
1339
1340	!
1341	!-- Initialize 2D index array indicating canopy top index.
1342	ALLOCATE( pch_index_ji(nysg:nyng,nxlg:nxrg) )
1343	pch_index_ji = 0
1344
1345	DO i = nxlg, nxrg
1346	DO j = nysg, nyng
1347	DO k = 0, pch_index
1348	IF ( lad_s(k,j,i) /= 0 ) pch_index_ji(j,i) = k
1349	ENDDO
1350	!
1351	!-- Check whether topography and local vegetation on top exceed
1352	!-- height of the model domain.
1353	IF ( topo_top_ind(j,i,0) + pch_index_ji(j,i) >= nzt + 1 ) THEN
1354	message_string = 'Local vegetation height on top of ' // &
1355	'topography exceeds height of model domain.'
1356	CALL message( 'pcm_init', 'PA0674', 2, 2, myid, 6, 0 )
1357	ENDIF
1358
1359	ENDDO
1360	ENDDO
1361	!
1362	!-- Calculate global pch_index value (index of top of plant canopy from ground)
1363	pch_index = MAXVAL( pch_index_ji )
1364	!
1365	!-- Exchange pch_index from all processors
1366	#if defined( __parallel )
1367	CALL MPI_ALLREDUCE( MPI_IN_PLACE, pch_index, 1, MPI_INTEGER, &
1368	MPI_MAX, comm2d, ierr)
1369	#endif
1370	!
1371	!-- Allocation of arrays pcm_heating_rate, pcm_transpiration_rate and pcm_latent_rate
1372	ALLOCATE( pcm_heating_rate(0:pch_index,nysg:nyng,nxlg:nxrg) )
1373	pcm_heating_rate = 0.0_wp
1374
1375	IF ( humidity ) THEN
1376	ALLOCATE( pcm_transpiration_rate(0:pch_index,nysg:nyng,nxlg:nxrg) )
1377	pcm_transpiration_rate = 0.0_wp
1378	ALLOCATE( pcm_latent_rate(0:pch_index,nysg:nyng,nxlg:nxrg) )
1379	pcm_latent_rate = 0.0_wp
1380	ENDIF
1381	!
1382	!-- Initialization of the canopy heat source distribution due to heating
1383	!-- of the canopy layers by incoming solar radiation, in case that a non-zero
1384	!-- value is set for the canopy top heat flux (cthf), which equals the
1385	!-- available net radiation at canopy top.
1386	!-- The heat source distribution is calculated by a decaying exponential
1387	!-- function of the downward cumulative leaf area index (cum_lai_hf),
1388	!-- assuming that the foliage inside the plant canopy is heated by solar
1389	!-- radiation penetrating the canopy layers according to the distribution
1390	!-- of net radiation as suggested by Brown & Covey (1966; Agric. Meteorol. 3,
1391	!-- 73â96). This approach has been applied e.g. by Shaw & Schumann (1992;
1392	!-- Bound.-Layer Meteorol. 61, 47â64).
1393	!-- When using the radiation_interactions, canopy heating (pcm_heating_rate)
1394	!-- and plant canopy transpiration (pcm_transpiration_rate, pcm_latent_rate)
1395	!-- are calculated in the RTM after the calculation of radiation.
1396	IF ( cthf /= 0.0_wp ) THEN
1397
1398	ALLOCATE( cum_lai_hf(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1399	!
1400	!-- Piecewise calculation of the cumulative leaf area index by vertical
1401	!-- integration of the leaf area density
1402	cum_lai_hf(:,:,:) = 0.0_wp
1403	DO i = nxlg, nxrg
1404	DO j = nysg, nyng
1405	DO k = pch_index_ji(j,i)-1, 0, -1
1406	IF ( k == pch_index_ji(j,i)-1 ) THEN
1407	cum_lai_hf(k,j,i) = cum_lai_hf(k+1,j,i) + &
1408	( 0.5_wp * lad_s(k+1,j,i) * &
1409	( zw(k+1) - zu(k+1) ) ) + &
1410	( 0.5_wp * ( 0.5_wp * ( lad_s(k+1,j,i) + &
1411	lad_s(k,j,i) ) + &
1412	lad_s(k+1,j,i) ) * &
1413	( zu(k+1) - zw(k) ) )
1414	ELSE
1415	cum_lai_hf(k,j,i) = cum_lai_hf(k+1,j,i) + &
1416	( 0.5_wp * ( 0.5_wp * ( lad_s(k+2,j,i) + &
1417	lad_s(k+1,j,i) ) + &
1418	lad_s(k+1,j,i) ) * &
1419	( zw(k+1) - zu(k+1) ) ) + &
1420	( 0.5_wp * ( 0.5_wp * ( lad_s(k+1,j,i) + &
1421	lad_s(k,j,i) ) + &
1422	lad_s(k+1,j,i) ) * &
1423	( zu(k+1) - zw(k) ) )
1424	ENDIF
1425	ENDDO
1426	ENDDO
1427	ENDDO
1428
1429	!
1430	!-- In areas with canopy the surface value of the canopy heat
1431	!-- flux distribution overrides the surface heat flux (shf)
1432	!-- Start with default surface type
1433	DO m = 1, surf_def_h(0)%ns
1434	i = surf_def_h(0)%i(m)
1435	j = surf_def_h(0)%j(m)
1436	IF ( cum_lai_hf(0,j,i) /= 0.0_wp ) &
1437	surf_def_h(0)%shf(m) = cthf * exp( -ext_coef * cum_lai_hf(0,j,i) )
1438	ENDDO
1439	!
1440	!-- Natural surfaces
1441	DO m = 1, surf_lsm_h%ns
1442	i = surf_lsm_h%i(m)
1443	j = surf_lsm_h%j(m)
1444	IF ( cum_lai_hf(0,j,i) /= 0.0_wp ) &
1445	surf_lsm_h%shf(m) = cthf * exp( -ext_coef * cum_lai_hf(0,j,i) )
1446	ENDDO
1447	!
1448	!-- Urban surfaces
1449	DO m = 1, surf_usm_h%ns
1450	i = surf_usm_h%i(m)
1451	j = surf_usm_h%j(m)
1452	IF ( cum_lai_hf(0,j,i) /= 0.0_wp ) &
1453	surf_usm_h%shf(m) = cthf * exp( -ext_coef * cum_lai_hf(0,j,i) )
1454	ENDDO
1455	!
1456	!
1457	!-- Calculation of the heating rate (K/s) within the different layers of
1458	!-- the plant canopy. Calculation is only necessary in areas covered with
1459	!-- canopy.
1460	!-- Within the different canopy layers the plant-canopy heating
1461	!-- rate (pcm_heating_rate) is calculated as the vertical
1462	!-- divergence of the canopy heat fluxes at the top and bottom
1463	!-- of the respective layer
1464	DO i = nxlg, nxrg
1465	DO j = nysg, nyng
1466	DO k = 1, pch_index_ji(j,i)
1467	IF ( cum_lai_hf(0,j,i) /= 0.0_wp ) THEN
1468	pcm_heating_rate(k,j,i) = cthf * &
1469	( exp(-ext_coef*cum_lai_hf(k,j,i)) - &
1470	exp(-ext_coef*cum_lai_hf(k-1,j,i) ) ) / dzw(k)
1471	ENDIF
1472	ENDDO
1473	ENDDO
1474	ENDDO
1475
1476	ENDIF
1477
1478	IF ( debug_output ) CALL debug_message( 'pcm_init', 'end' )
1479
1480	END SUBROUTINE pcm_init
1481
1482
1483	!------------------------------------------------------------------------------!
1484	! Description:
1485	! ------------
1486	!> Parin for &plant_canopy_parameters for plant canopy model
1487	!------------------------------------------------------------------------------!
1488	SUBROUTINE pcm_parin
1489
1490	CHARACTER (LEN=80) :: line !< dummy string that contains the current line of the parameter file
1491
1492	NAMELIST /plant_canopy_parameters/ &
1493	alpha_lad, beta_lad, canopy_drag_coeff, &
1494	canopy_mode, cthf, &
1495	lad_surface, lad_type_coef, &
1496	lad_vertical_gradient, &
1497	lad_vertical_gradient_level, &
1498	lai_beta, &
1499	leaf_scalar_exch_coeff, &
1500	leaf_surface_conc, pch_index, &
1501	plant_canopy_transpiration
1502
1503	NAMELIST /canopy_par/ alpha_lad, beta_lad, canopy_drag_coeff, &
1504	canopy_mode, cthf, &
1505	lad_surface, lad_type_coef, &
1506	lad_vertical_gradient, &
1507	lad_vertical_gradient_level, &
1508	lai_beta, &
1509	leaf_scalar_exch_coeff, &
1510	leaf_surface_conc, pch_index, &
1511	plant_canopy_transpiration
1512
1513	line = ' '
1514
1515	!
1516	!-- Try to find plant-canopy model package
1517	REWIND ( 11 )
1518	line = ' '
1519	DO WHILE ( INDEX( line, '&plant_canopy_parameters' ) == 0 )
1520	READ ( 11, '(A)', END=12 ) line
1521	ENDDO
1522	BACKSPACE ( 11 )
1523
1524	!
1525	!-- Read user-defined namelist
1526	READ ( 11, plant_canopy_parameters, ERR = 10 )
1527
1528	!
1529	!-- Set flag that indicates that the plant-canopy model is switched on
1530	plant_canopy = .TRUE.
1531
1532	GOTO 14
1533
1534	10 BACKSPACE( 11 )
1535	READ( 11 , '(A)') line
1536	CALL parin_fail_message( 'plant_canopy_parameters', line )
1537	!
1538	!-- Try to find old namelist
1539	12 REWIND ( 11 )
1540	line = ' '
1541	DO WHILE ( INDEX( line, '&canopy_par' ) == 0 )
1542	READ ( 11, '(A)', END=14 ) line
1543	ENDDO
1544	BACKSPACE ( 11 )
1545
1546	!
1547	!-- Read user-defined namelist
1548	READ ( 11, canopy_par, ERR = 13, END = 14 )
1549
1550	message_string = 'namelist canopy_par is deprecated and will be ' // &
1551	'removed in near future. Please use namelist ' // &
1552	'plant_canopy_parameters instead'
1553	CALL message( 'pcm_parin', 'PA0487', 0, 1, 0, 6, 0 )
1554
1555	!
1556	!-- Set flag that indicates that the plant-canopy model is switched on
1557	plant_canopy = .TRUE.
1558
1559	GOTO 14
1560
1561	13 BACKSPACE( 11 )
1562	READ( 11 , '(A)') line
1563	CALL parin_fail_message( 'canopy_par', line )
1564
1565	14 CONTINUE
1566
1567	END SUBROUTINE pcm_parin
1568
1569
1570	!------------------------------------------------------------------------------!
1571	! Description:
1572	! ------------
1573	!
1574	!> Loads 3D plant canopy data from file. File format is as follows:
1575	!>
1576	!> num_levels
1577	!> dtype,x,y,pctype,value(nzb),value(nzb+1), ... ,value(nzb+num_levels-1)
1578	!> dtype,x,y,pctype,value(nzb),value(nzb+1), ... ,value(nzb+num_levels-1)
1579	!> dtype,x,y,pctype,value(nzb),value(nzb+1), ... ,value(nzb+num_levels-1)
1580	!> ...
1581	!>
1582	!> i.e. first line determines number of levels and further lines represent plant
1583	!> canopy data, one line per column and variable. In each data line,
1584	!> dtype represents variable to be set:
1585	!>
1586	!> dtype=1: leaf area density (lad_s)
1587	!> dtype=2....n: some additional plant canopy input data quantity
1588	!>
1589	!> Zeros are added automatically above num_levels until top of domain. Any
1590	!> non-specified (x,y) columns have zero values as default.
1591	!------------------------------------------------------------------------------!
1592	SUBROUTINE pcm_read_plant_canopy_3d
1593
1594	USE exchange_horiz_mod, &
1595	ONLY: exchange_horiz
1596
1597	INTEGER(iwp) :: dtype !< type of input data (1=lad)
1598	INTEGER(iwp) :: pctype !< type of plant canopy (deciduous,non-deciduous,...)
1599	INTEGER(iwp) :: i, j !< running index
1600	INTEGER(iwp) :: nzp !< number of vertical layers of plant canopy
1601	INTEGER(iwp) :: nzpltop !<
1602	INTEGER(iwp) :: nzpl !<
1603	INTEGER(iwp) :: kk !<
1604
1605	REAL(wp), DIMENSION(:), ALLOCATABLE :: col !< vertical column of input data
1606
1607	!
1608	!-- Initialize lad_s array
1609	lad_s = 0.0_wp
1610
1611	!
1612	!-- Open and read plant canopy input data
1613	OPEN(152, FILE='PLANT_CANOPY_DATA_3D' // TRIM( coupling_char ), &
1614	ACCESS='SEQUENTIAL', ACTION='READ', STATUS='OLD', &
1615	FORM='FORMATTED', ERR=515)
1616	READ(152, *, ERR=516, END=517) nzp !< read first line = number of vertical layers
1617	nzpltop = MIN(nzt+1, nzb+nzp-1)
1618	nzpl = nzpltop - nzb + 1 !< no. of layers to assign
1619	ALLOCATE( col(0:nzp-1) )
1620
1621	DO
1622	READ(152, *, ERR=516, END=517) dtype, i, j, pctype, col(:)
1623	IF ( i < nxlg .OR. i > nxrg .OR. j < nysg .OR. j > nyng ) CYCLE
1624
1625	SELECT CASE (dtype)
1626	CASE( 1 ) !< leaf area density
1627	!
1628	!-- This is just the pure canopy layer assumed to be grounded to
1629	!-- a flat domain surface. At locations where plant canopy sits
1630	!-- on top of any kind of topography, the vertical plant column
1631	!-- must be "lifted", which is done in SUBROUTINE pcm_tendency.
1632	IF ( pctype < 0 .OR. pctype > 10 ) THEN !< incorrect plant canopy type
1633	WRITE( message_string, * ) 'Incorrect type of plant canopy. ' // &
1634	'Allowed values 0 <= pctype <= 10, ' // &
1635	'but pctype is ', pctype
1636	CALL message( 'pcm_read_plant_canopy_3d', 'PA0349', 1, 2, 0, 6, 0 )
1637	ENDIF
1638	kk = topo_top_ind(j,i,0)
1639	lad_s(nzb:nzpltop-kk, j, i) = col(kk:nzpl-1)*lad_type_coef(pctype)
1640	CASE DEFAULT
1641	WRITE(message_string, '(a,i2,a)') &
1642	'Unknown record type in file PLANT_CANOPY_DATA_3D: "', dtype, '"'
1643	CALL message( 'pcm_read_plant_canopy_3d', 'PA0530', 1, 2, 0, 6, 0 )
1644	END SELECT
1645	ENDDO
1646
1647	515 message_string = 'error opening file PLANT_CANOPY_DATA_3D'
1648	CALL message( 'pcm_read_plant_canopy_3d', 'PA0531', 1, 2, 0, 6, 0 )
1649
1650	516 message_string = 'error reading file PLANT_CANOPY_DATA_3D'
1651	CALL message( 'pcm_read_plant_canopy_3d', 'PA0532', 1, 2, 0, 6, 0 )
1652
1653	517 CLOSE(152)
1654	DEALLOCATE( col )
1655
1656	CALL exchange_horiz( lad_s, nbgp )
1657
1658	END SUBROUTINE pcm_read_plant_canopy_3d
1659
1660	!------------------------------------------------------------------------------!
1661	! Description:
1662	! ------------
1663	!> Read module-specific global restart data (Fortran binary format).
1664	!------------------------------------------------------------------------------!
1665	SUBROUTINE pcm_rrd_global_ftn( found )
1666
1667	LOGICAL, INTENT(OUT) :: found
1668
1669	found = .TRUE.
1670
1671	SELECT CASE ( restart_string(1:length) )
1672
1673	CASE ( 'pch_index' )
1674	READ ( 13 ) pch_index
1675
1676	CASE DEFAULT
1677
1678	found = .FALSE.
1679
1680	END SELECT
1681
1682	END SUBROUTINE pcm_rrd_global_ftn
1683
1684	!------------------------------------------------------------------------------!
1685	! Description:
1686	! ------------
1687	!> Read module-specific global restart data (MPI-IO).
1688	!------------------------------------------------------------------------------!
1689	SUBROUTINE pcm_rrd_global_mpi
1690
1691	CALL rrd_mpi_io( 'pch_index', pch_index )
1692
1693	END SUBROUTINE pcm_rrd_global_mpi
1694
1695
1696	!------------------------------------------------------------------------------!
1697	! Description:
1698	! ------------
1699	!> Read module-specific local restart data arrays (Fortran binary format).
1700	!------------------------------------------------------------------------------!
1701	SUBROUTINE pcm_rrd_local_ftn( k, nxlf, nxlc, nxl_on_file, nxrf, nxrc, &
1702	nxr_on_file, nynf, nync, nyn_on_file, nysf, &
1703	nysc, nys_on_file, found )
1704
1705	INTEGER(iwp) :: k !<
1706	INTEGER(iwp) :: nxlc !<
1707	INTEGER(iwp) :: nxlf !<
1708	INTEGER(iwp) :: nxl_on_file !<
1709	INTEGER(iwp) :: nxrc !<
1710	INTEGER(iwp) :: nxrf !<
1711	INTEGER(iwp) :: nxr_on_file !<
1712	INTEGER(iwp) :: nync !<
1713	INTEGER(iwp) :: nynf !<
1714	INTEGER(iwp) :: nyn_on_file !<
1715	INTEGER(iwp) :: nysc !<
1716	INTEGER(iwp) :: nysf !<
1717	INTEGER(iwp) :: nys_on_file !<
1718
1719	LOGICAL, INTENT(OUT) :: found
1720
1721	REAL(wp), DIMENSION(0:pch_index, &
1722	nys_on_file-nbgp:nyn_on_file+nbgp, &
1723	nxl_on_file-nbgp:nxr_on_file+nbgp) :: tmp_3d2 !< temporary 3D array for entire vertical
1724	!< extension of canopy layer
1725	found = .TRUE.
1726
1727	SELECT CASE ( restart_string(1:length) )
1728
1729	CASE ( 'pcm_heatrate_av' )
1730	IF ( .NOT. ALLOCATED( pcm_heatrate_av ) ) THEN
1731	ALLOCATE( pcm_heatrate_av(0:pch_index,nysg:nyng,nxlg:nxrg) )
1732	pcm_heatrate_av = 0.0_wp
1733	ENDIF
1734	IF ( k == 1 ) READ ( 13 ) tmp_3d2
1735	pcm_heatrate_av(0:pch_index,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
1736	tmp_3d2(0:pch_index,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
1737
1738	CASE ( 'pcm_latentrate_av' )
1739	IF ( .NOT. ALLOCATED( pcm_latentrate_av ) ) THEN
1740	ALLOCATE( pcm_latentrate_av(0:pch_index,nysg:nyng,nxlg:nxrg) )
1741	pcm_latentrate_av = 0.0_wp
1742	ENDIF
1743	IF ( k == 1 ) READ ( 13 ) tmp_3d2
1744	pcm_latentrate_av(0:pch_index,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
1745	tmp_3d2(0:pch_index,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
1746
1747	CASE ( 'pcm_transpirationrate_av' )
1748	IF ( .NOT. ALLOCATED( pcm_transpirationrate_av ) ) THEN
1749	ALLOCATE( pcm_transpirationrate_av(0:pch_index,nysg:nyng,nxlg:nxrg) )
1750	pcm_transpirationrate_av = 0.0_wp
1751	ENDIF
1752	IF ( k == 1 ) READ ( 13 ) tmp_3d2
1753	pcm_transpirationrate_av(0:pch_index,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
1754	tmp_3d2(0:pch_index,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
1755
1756	CASE DEFAULT
1757
1758	found = .FALSE.
1759
1760	END SELECT
1761
1762	END SUBROUTINE pcm_rrd_local_ftn
1763
1764
1765	!------------------------------------------------------------------------------!
1766	! Description:
1767	! ------------
1768	!> Read module-specific local restart data arrays (MPI-IO).
1769	!------------------------------------------------------------------------------!
1770	SUBROUTINE pcm_rrd_local_mpi
1771
1772	IMPLICIT NONE
1773
1774	LOGICAL :: array_found !<
1775
1776	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: tmp_3d !< temporary array to store pcm data with
1777	!< non-standard vertical index bounds
1778
1779	!
1780	!-- Plant canopy arrays have non standard reduced vertical index bounds. They are stored with
1781	!-- full vertical bounds (bzb:nzt+1) in the restart file and must be re-stored after reading.
1782	CALL rd_mpi_io_check_array( 'pcm_heatrate_av' , found = array_found )
1783	IF ( array_found ) THEN
1784	IF ( .NOT. ALLOCATED( pcm_heatrate_av ) ) THEN
1785	ALLOCATE( pcm_heatrate_av(nzb:pch_index,nysg:nyng,nxlg:nxrg) )
1786	ENDIF
1787	ALLOCATE( tmp_3d(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1788	CALL rrd_mpi_io( 'pcm_heatrate_av', tmp_3d )
1789	pcm_heatrate_av = tmp_3d(nzb:pch_index,:,:)
1790	DEALLOCATE( tmp_3d )
1791	ENDIF
1792
1793	CALL rd_mpi_io_check_array( 'pcm_latentrate_av' , found = array_found )
1794	IF ( array_found ) THEN
1795	IF ( .NOT. ALLOCATED( pcm_latentrate_av ) ) THEN
1796	ALLOCATE( pcm_latentrate_av(nzb:pch_index,nysg:nyng,nxlg:nxrg) )
1797	ENDIF
1798	ALLOCATE( tmp_3d(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1799	CALL rrd_mpi_io( 'pcm_latentrate_av', tmp_3d )
1800	pcm_latentrate_av = tmp_3d(nzb:pch_index,:,:)
1801	DEALLOCATE( tmp_3d )
1802	ENDIF
1803
1804	CALL rd_mpi_io_check_array( 'pcm_transpirationrate_av' , found = array_found )
1805	IF ( array_found ) THEN
1806	IF ( .NOT. ALLOCATED( pcm_transpirationrate_av ) ) THEN
1807	ALLOCATE( pcm_transpirationrate_av(nzb:pch_index,nysg:nyng,nxlg:nxrg) )
1808	ENDIF
1809	ALLOCATE( tmp_3d(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1810	CALL rrd_mpi_io( 'pcm_transpirationrate_av', tmp_3d )
1811	pcm_transpirationrate_av = tmp_3d(nzb:pch_index,:,:)
1812	DEALLOCATE( tmp_3d )
1813	ENDIF
1814
1815	END SUBROUTINE pcm_rrd_local_mpi
1816
1817
1818	!------------------------------------------------------------------------------!
1819	! Description:
1820	! ------------
1821	!> Calculation of the tendency terms, accounting for the effect of the plant
1822	!> canopy on momentum and scalar quantities.
1823	!>
1824	!> The canopy is located where the leaf area density lad_s(k,j,i) > 0.0
1825	!> (defined on scalar grid), as initialized in subroutine pcm_init.
1826	!> The lad on the w-grid is vertically interpolated from the surrounding
1827	!> lad_s. The upper boundary of the canopy is defined on the w-grid at
1828	!> k = pch_index. Here, the lad is zero.
1829	!>
1830	!> The canopy drag must be limited (previously accounted for by calculation of
1831	!> a limiting canopy timestep for the determination of the maximum LES timestep
1832	!> in subroutine timestep), since it is physically impossible that the canopy
1833	!> drag alone can locally change the sign of a velocity component. This
1834	!> limitation is realized by calculating preliminary tendencies and velocities.
1835	!> It is subsequently checked if the preliminary new velocity has a different
1836	!> sign than the current velocity. If so, the tendency is limited in a way that
1837	!> the velocity can at maximum be reduced to zero by the canopy drag.
1838	!>
1839	!>
1840	!> Call for all grid points
1841	!------------------------------------------------------------------------------!
1842	SUBROUTINE pcm_tendency( component )
1843
1844	INTEGER(iwp) :: component !< prognostic variable (u,v,w,pt,q,e)
1845	INTEGER(iwp) :: i !< running index
1846	INTEGER(iwp) :: j !< running index
1847	INTEGER(iwp) :: k !< running index
1848	INTEGER(iwp) :: kk !< running index for flat lad arrays
1849
1850	LOGICAL :: building_edge_e !< control flag indicating an eastward-facing building edge
1851	LOGICAL :: building_edge_n !< control flag indicating a north-facing building edge
1852	LOGICAL :: building_edge_s !< control flag indicating a south-facing building edge
1853	LOGICAL :: building_edge_w !< control flag indicating a westward-facing building edge
1854
1855	REAL(wp) :: ddt_3d !< inverse of the LES timestep (dt_3d)
1856	REAL(wp) :: lad_local !< local lad value
1857	REAL(wp) :: pre_tend !< preliminary tendency
1858	REAL(wp) :: pre_u !< preliminary u-value
1859	REAL(wp) :: pre_v !< preliminary v-value
1860	REAL(wp) :: pre_w !< preliminary w-value
1861
1862
1863	ddt_3d = 1.0_wp / dt_3d
1864
1865	!
1866	!-- Compute drag for the three velocity components and the SGS-TKE:
1867	SELECT CASE ( component )
1868
1869	!
1870	!-- u-component
1871	CASE ( 1 )
1872	DO i = nxlu, nxr
1873	DO j = nys, nyn
1874	!
1875	!-- Set control flags indicating east- and westward-orientated
1876	!-- building edges. Note, building_egde_w is set from the perspective
1877	!-- of the potential rooftop grid point, while building_edge_e is
1878	!-- is set from the perspective of the non-building grid point.
1879	building_edge_w = ANY( BTEST( wall_flags_total_0(:,j,i), 6 ) )&
1880	.AND. .NOT. ANY( BTEST( wall_flags_total_0(:,j,i-1), 6 ) )
1881	building_edge_e = ANY( BTEST( wall_flags_total_0(:,j,i-1), 6 ) )&
1882	.AND. .NOT. ANY( BTEST( wall_flags_total_0(:,j,i), 6 ) )
1883	!
1884	!-- Determine topography-top index on u-grid
1885	DO k = topo_top_ind(j,i,1)+1, topo_top_ind(j,i,1) + pch_index_ji(j,i)
1886
1887	kk = k - topo_top_ind(j,i,1) !- lad arrays are defined flat
1888	!
1889	!-- In order to create sharp boundaries of the plant canopy,
1890	!-- the lad on the u-grid at index (k,j,i) is equal to lad_s(k,j,i),
1891	!-- rather than being interpolated from the surrounding lad_s,
1892	!-- because this would yield smaller lad at the canopy boundaries
1893	!-- than inside of the canopy.
1894	!-- For the same reason, the lad at the rightmost(i+1)canopy
1895	!-- boundary on the u-grid equals lad_s(k,j,i), which is considered
1896	!-- in the next if-statement. Note, at left-sided building edges
1897	!-- this is not applied, here the LAD is equals the LAD at grid
1898	!-- point (k,j,i), in order to avoid that LAD is mistakenly mapped
1899	!-- on top of a roof where (usually) is no LAD is defined.
1900	lad_local = lad_s(kk,j,i)
1901	IF ( lad_local == 0.0_wp .AND. lad_s(kk,j,i-1) > 0.0_wp &
1902	.AND. .NOT. building_edge_w ) lad_local = lad_s(kk,j,i-1)
1903	!
1904	!-- In order to avoid that LAD is mistakenly considered at right-
1905	!-- sided building edges (here the topography-top index for the
1906	!-- u-component at index j,i is still on the building while the
1907	!-- topography top for the scalar isn't), LAD is taken from grid
1908	!-- point (j,i-1).
1909	IF ( lad_local > 0.0_wp .AND. lad_s(kk,j,i-1) == 0.0_wp &
1910	.AND. building_edge_e ) lad_local = lad_s(kk,j,i-1)
1911
1912	pre_tend = 0.0_wp
1913	pre_u = 0.0_wp
1914	!
1915	!-- Calculate preliminary value (pre_tend) of the tendency
1916	pre_tend = - canopy_drag_coeff * &
1917	lad_local * &
1918	SQRT( u(k,j,i)**2 + &
1919	( 0.25_wp * ( v(k,j,i-1) + &
1920	v(k,j,i) + &
1921	v(k,j+1,i) + &
1922	v(k,j+1,i-1) ) &
1923	)**2 + &
1924	( 0.25_wp * ( w(k-1,j,i-1) + &
1925	w(k-1,j,i) + &
1926	w(k,j,i-1) + &
1927	w(k,j,i) ) &
1928	)**2 &
1929	) * &
1930	u(k,j,i)
1931
1932	!
1933	!-- Calculate preliminary new velocity, based on pre_tend
1934	pre_u = u(k,j,i) + dt_3d * pre_tend
1935	!
1936	!-- Compare sign of old velocity and new preliminary velocity,
1937	!-- and in case the signs are different, limit the tendency
1938	IF ( SIGN(pre_u,u(k,j,i)) /= pre_u ) THEN
1939	pre_tend = - u(k,j,i) * ddt_3d
1940	ENDIF
1941	!
1942	!-- Calculate final tendency
1943	tend(k,j,i) = tend(k,j,i) + pre_tend
1944
1945	ENDDO
1946	ENDDO
1947	ENDDO
1948
1949	!
1950	!-- v-component
1951	CASE ( 2 )
1952	DO i = nxl, nxr
1953	DO j = nysv, nyn
1954	!
1955	!-- Set control flags indicating north- and southward-orientated
1956	!-- building edges. Note, building_egde_s is set from the perspective
1957	!-- of the potential rooftop grid point, while building_edge_n is
1958	!-- is set from the perspective of the non-building grid point.
1959	building_edge_s = ANY( BTEST( wall_flags_total_0(:,j,i), 6 ) )&
1960	.AND. .NOT. ANY( BTEST( wall_flags_total_0(:,j-1,i), 6 ) )
1961	building_edge_n = ANY( BTEST( wall_flags_total_0(:,j-1,i), 6 ) )&
1962	.AND. .NOT. ANY( BTEST( wall_flags_total_0(:,j,i), 6 ) )
1963	!
1964	!-- Determine topography-top index on v-grid
1965	DO k = topo_top_ind(j,i,2)+1, topo_top_ind(j,i,2) + pch_index_ji(j,i)
1966
1967	kk = k - topo_top_ind(j,i,2) !- lad arrays are defined flat
1968	!
1969	!-- In order to create sharp boundaries of the plant canopy,
1970	!-- the lad on the v-grid at index (k,j,i) is equal to lad_s(k,j,i),
1971	!-- rather than being interpolated from the surrounding lad_s,
1972	!-- because this would yield smaller lad at the canopy boundaries
1973	!-- than inside of the canopy.
1974	!-- For the same reason, the lad at the northmost(j+1)canopy
1975	!-- boundary on the v-grid equals lad_s(k,j,i), which is considered
1976	!-- in the next if-statement. Note, at left-sided building edges
1977	!-- this is not applied, here the LAD is equals the LAD at grid
1978	!-- point (k,j,i), in order to avoid that LAD is mistakenly mapped
1979	!-- on top of a roof where (usually) is no LAD is defined.
1980	lad_local = lad_s(kk,j,i)
1981	IF ( lad_local == 0.0_wp .AND. lad_s(kk,j-1,i) > 0.0_wp &
1982	.AND. .NOT. building_edge_s ) lad_local = lad_s(kk,j-1,i)
1983	!
1984	!-- In order to avoid that LAD is mistakenly considered at right-
1985	!-- sided building edges (here the topography-top index for the
1986	!-- u-component at index j,i is still on the building while the
1987	!-- topography top for the scalar isn't), LAD is taken from grid
1988	!-- point (j,i-1).
1989	IF ( lad_local > 0.0_wp .AND. lad_s(kk,j-1,i) == 0.0_wp &
1990	.AND. building_edge_n ) lad_local = lad_s(kk,j-1,i)
1991
1992	pre_tend = 0.0_wp
1993	pre_v = 0.0_wp
1994	!
1995	!-- Calculate preliminary value (pre_tend) of the tendency
1996	pre_tend = - canopy_drag_coeff * &
1997	lad_local * &
1998	SQRT( ( 0.25_wp * ( u(k,j-1,i) + &
1999	u(k,j-1,i+1) + &
2000	u(k,j,i) + &
2001	u(k,j,i+1) ) &
2002	)**2 + &
2003	v(k,j,i)**2 + &
2004	( 0.25_wp * ( w(k-1,j-1,i) + &
2005	w(k-1,j,i) + &
2006	w(k,j-1,i) + &
2007	w(k,j,i) ) &
2008	)**2 &
2009	) * &
2010	v(k,j,i)
2011
2012	!
2013	!-- Calculate preliminary new velocity, based on pre_tend
2014	pre_v = v(k,j,i) + dt_3d * pre_tend
2015	!
2016	!-- Compare sign of old velocity and new preliminary velocity,
2017	!-- and in case the signs are different, limit the tendency
2018	IF ( SIGN(pre_v,v(k,j,i)) /= pre_v ) THEN
2019	pre_tend = - v(k,j,i) * ddt_3d
2020	ELSE
2021	pre_tend = pre_tend
2022	ENDIF
2023	!
2024	!-- Calculate final tendency
2025	tend(k,j,i) = tend(k,j,i) + pre_tend
2026
2027	ENDDO
2028	ENDDO
2029	ENDDO
2030
2031	!
2032	!-- w-component
2033	CASE ( 3 )
2034	DO i = nxl, nxr
2035	DO j = nys, nyn
2036	!
2037	!-- Determine topography-top index on w-grid
2038	DO k = topo_top_ind(j,i,3)+1, topo_top_ind(j,i,3) + pch_index_ji(j,i) - 1
2039
2040	kk = k - topo_top_ind(j,i,3) !- lad arrays are defined flat
2041
2042	pre_tend = 0.0_wp
2043	pre_w = 0.0_wp
2044	!
2045	!-- Calculate preliminary value (pre_tend) of the tendency
2046	pre_tend = - canopy_drag_coeff * &
2047	(0.5_wp * &
2048	( lad_s(kk+1,j,i) + lad_s(kk,j,i) )) * &
2049	SQRT( ( 0.25_wp * ( u(k,j,i) + &
2050	u(k,j,i+1) + &
2051	u(k+1,j,i) + &
2052	u(k+1,j,i+1) ) &
2053	)**2 + &
2054	( 0.25_wp * ( v(k,j,i) + &
2055	v(k,j+1,i) + &
2056	v(k+1,j,i) + &
2057	v(k+1,j+1,i) ) &
2058	)**2 + &
2059	w(k,j,i)**2 &
2060	) * &
2061	w(k,j,i)
2062	!
2063	!-- Calculate preliminary new velocity, based on pre_tend
2064	pre_w = w(k,j,i) + dt_3d * pre_tend
2065	!
2066	!-- Compare sign of old velocity and new preliminary velocity,
2067	!-- and in case the signs are different, limit the tendency
2068	IF ( SIGN(pre_w,w(k,j,i)) /= pre_w ) THEN
2069	pre_tend = - w(k,j,i) * ddt_3d
2070	ELSE
2071	pre_tend = pre_tend
2072	ENDIF
2073	!
2074	!-- Calculate final tendency
2075	tend(k,j,i) = tend(k,j,i) + pre_tend
2076
2077	ENDDO
2078	ENDDO
2079	ENDDO
2080
2081	!
2082	!-- potential temperature
2083	CASE ( 4 )
2084	IF ( humidity ) THEN
2085	DO i = nxl, nxr
2086	DO j = nys, nyn
2087	!-- Determine topography-top index on scalar-grid
2088	DO k = topo_top_ind(j,i,0)+1, topo_top_ind(j,i,0) + pch_index_ji(j,i)
2089	kk = k - topo_top_ind(j,i,0) !- lad arrays are defined flat
2090	tend(k,j,i) = tend(k,j,i) + pcm_heating_rate(kk,j,i) - pcm_latent_rate(kk,j,i)
2091	ENDDO
2092	ENDDO
2093	ENDDO
2094	ELSE
2095	DO i = nxl, nxr
2096	DO j = nys, nyn
2097	!-- Determine topography-top index on scalar-grid
2098	DO k = topo_top_ind(j,i,0)+1, topo_top_ind(j,i,0) + pch_index_ji(j,i)
2099	kk = k - topo_top_ind(j,i,0) !- lad arrays are defined flat
2100	tend(k,j,i) = tend(k,j,i) + pcm_heating_rate(kk,j,i)
2101	ENDDO
2102	ENDDO
2103	ENDDO
2104	ENDIF
2105
2106	!
2107	!-- humidity
2108	CASE ( 5 )
2109	DO i = nxl, nxr
2110	DO j = nys, nyn
2111	!
2112	!-- Determine topography-top index on scalar-grid
2113	DO k = topo_top_ind(j,i,0)+1, topo_top_ind(j,i,0) + pch_index_ji(j,i)
2114
2115	kk = k - topo_top_ind(j,i,0) !- lad arrays are defined flat
2116
2117	IF ( .NOT. plant_canopy_transpiration ) THEN
2118	! pcm_transpiration_rate is calculated in radiation model
2119	! in case of plant_canopy_transpiration = .T.
2120	! to include also the dependecy to the radiation
2121	! in the plant canopy box
2122	pcm_transpiration_rate(kk,j,i) = - leaf_scalar_exch_coeff &
2123	* lad_s(kk,j,i) * &
2124	SQRT( ( 0.5_wp * ( u(k,j,i) + &
2125	u(k,j,i+1) ) &
2126	)**2 + &
2127	( 0.5_wp * ( v(k,j,i) + &
2128	v(k,j+1,i) ) &
2129	)**2 + &
2130	( 0.5_wp * ( w(k-1,j,i) + &
2131	w(k,j,i) ) &
2132	)**2 &
2133	) * &
2134	( q(k,j,i) - leaf_surface_conc )
2135	ENDIF
2136
2137	tend(k,j,i) = tend(k,j,i) + pcm_transpiration_rate(kk,j,i)
2138	ENDDO
2139	ENDDO
2140	ENDDO
2141
2142	!
2143	!-- sgs-tke
2144	CASE ( 6 )
2145	DO i = nxl, nxr
2146	DO j = nys, nyn
2147	!
2148	!-- Determine topography-top index on scalar-grid
2149	DO k = topo_top_ind(j,i,0)+1, topo_top_ind(j,i,0) + pch_index_ji(j,i)
2150
2151	kk = k - topo_top_ind(j,i,0) !- lad arrays are defined flat
2152	tend(k,j,i) = tend(k,j,i) - &
2153	2.0_wp * canopy_drag_coeff * &
2154	lad_s(kk,j,i) * &
2155	SQRT( ( 0.5_wp * ( u(k,j,i) + &
2156	u(k,j,i+1) ) &
2157	)**2 + &
2158	( 0.5_wp * ( v(k,j,i) + &
2159	v(k,j+1,i) ) &
2160	)**2 + &
2161	( 0.5_wp * ( w(k,j,i) + &
2162	w(k+1,j,i) ) &
2163	)**2 &
2164	) * &
2165	e(k,j,i)
2166	ENDDO
2167	ENDDO
2168	ENDDO
2169	!
2170	!-- scalar concentration
2171	CASE ( 7 )
2172	DO i = nxl, nxr
2173	DO j = nys, nyn
2174	!
2175	!-- Determine topography-top index on scalar-grid
2176	DO k = topo_top_ind(j,i,0)+1, topo_top_ind(j,i,0) + pch_index_ji(j,i)
2177
2178	kk = k - topo_top_ind(j,i,0) !- lad arrays are defined flat
2179	tend(k,j,i) = tend(k,j,i) - &
2180	leaf_scalar_exch_coeff * &
2181	lad_s(kk,j,i) * &
2182	SQRT( ( 0.5_wp * ( u(k,j,i) + &
2183	u(k,j,i+1) ) &
2184	)**2 + &
2185	( 0.5_wp * ( v(k,j,i) + &
2186	v(k,j+1,i) ) &
2187	)**2 + &
2188	( 0.5_wp * ( w(k-1,j,i) + &
2189	w(k,j,i) ) &
2190	)**2 &
2191	) * &
2192	( s(k,j,i) - leaf_surface_conc )
2193	ENDDO
2194	ENDDO
2195	ENDDO
2196
2197
2198
2199	CASE DEFAULT
2200
2201	WRITE( message_string, * ) 'wrong component: ', component
2202	CALL message( 'pcm_tendency', 'PA0279', 1, 2, 0, 6, 0 )
2203
2204	END SELECT
2205
2206	END SUBROUTINE pcm_tendency
2207
2208
2209	!------------------------------------------------------------------------------!
2210	! Description:
2211	! ------------
2212	!> Calculation of the tendency terms, accounting for the effect of the plant
2213	!> canopy on momentum and scalar quantities.
2214	!>
2215	!> The canopy is located where the leaf area density lad_s(k,j,i) > 0.0
2216	!> (defined on scalar grid), as initialized in subroutine pcm_init.
2217	!> The lad on the w-grid is vertically interpolated from the surrounding
2218	!> lad_s. The upper boundary of the canopy is defined on the w-grid at
2219	!> k = pch_index. Here, the lad is zero.
2220	!>
2221	!> The canopy drag must be limited (previously accounted for by calculation of
2222	!> a limiting canopy timestep for the determination of the maximum LES timestep
2223	!> in subroutine timestep), since it is physically impossible that the canopy
2224	!> drag alone can locally change the sign of a velocity component. This
2225	!> limitation is realized by calculating preliminary tendencies and velocities.
2226	!> It is subsequently checked if the preliminary new velocity has a different
2227	!> sign than the current velocity. If so, the tendency is limited in a way that
2228	!> the velocity can at maximum be reduced to zero by the canopy drag.
2229	!>
2230	!>
2231	!> Call for grid point i,j
2232	!------------------------------------------------------------------------------!
2233	SUBROUTINE pcm_tendency_ij( i, j, component )
2234
2235	INTEGER(iwp) :: component !< prognostic variable (u,v,w,pt,q,e)
2236	INTEGER(iwp) :: i !< running index
2237	INTEGER(iwp) :: j !< running index
2238	INTEGER(iwp) :: k !< running index
2239	INTEGER(iwp) :: kk !< running index for flat lad arrays
2240
2241	LOGICAL :: building_edge_e !< control flag indicating an eastward-facing building edge
2242	LOGICAL :: building_edge_n !< control flag indicating a north-facing building edge
2243	LOGICAL :: building_edge_s !< control flag indicating a south-facing building edge
2244	LOGICAL :: building_edge_w !< control flag indicating a westward-facing building edge
2245
2246	REAL(wp) :: ddt_3d !< inverse of the LES timestep (dt_3d)
2247	REAL(wp) :: lad_local !< local lad value
2248	REAL(wp) :: pre_tend !< preliminary tendency
2249	REAL(wp) :: pre_u !< preliminary u-value
2250	REAL(wp) :: pre_v !< preliminary v-value
2251	REAL(wp) :: pre_w !< preliminary w-value
2252
2253
2254	ddt_3d = 1.0_wp / dt_3d
2255	!
2256	!-- Compute drag for the three velocity components and the SGS-TKE
2257	SELECT CASE ( component )
2258
2259	!
2260	!-- u-component
2261	CASE ( 1 )
2262	!
2263	!-- Set control flags indicating east- and westward-orientated
2264	!-- building edges. Note, building_egde_w is set from the perspective
2265	!-- of the potential rooftop grid point, while building_edge_e is
2266	!-- is set from the perspective of the non-building grid point.
2267	building_edge_w = ANY( BTEST( wall_flags_total_0(:,j,i), 6 ) ) .AND. &
2268	.NOT. ANY( BTEST( wall_flags_total_0(:,j,i-1), 6 ) )
2269	building_edge_e = ANY( BTEST( wall_flags_total_0(:,j,i-1), 6 ) ) .AND. &
2270	.NOT. ANY( BTEST( wall_flags_total_0(:,j,i), 6 ) )
2271	!
2272	!-- Determine topography-top index on u-grid
2273	DO k = topo_top_ind(j,i,1) + 1, topo_top_ind(j,i,1) + pch_index_ji(j,i)
2274
2275	kk = k - topo_top_ind(j,i,1) !- lad arrays are defined flat
2276
2277	!
2278	!-- In order to create sharp boundaries of the plant canopy,
2279	!-- the lad on the u-grid at index (k,j,i) is equal to lad_s(k,j,i),
2280	!-- rather than being interpolated from the surrounding lad_s,
2281	!-- because this would yield smaller lad at the canopy boundaries
2282	!-- than inside of the canopy.
2283	!-- For the same reason, the lad at the rightmost(i+1)canopy
2284	!-- boundary on the u-grid equals lad_s(k,j,i), which is considered
2285	!-- in the next if-statement. Note, at left-sided building edges
2286	!-- this is not applied, here the LAD is equals the LAD at grid
2287	!-- point (k,j,i), in order to avoid that LAD is mistakenly mapped
2288	!-- on top of a roof where (usually) is no LAD is defined.
2289	lad_local = lad_s(kk,j,i)
2290	IF ( lad_local == 0.0_wp .AND. lad_s(kk,j,i-1) > 0.0_wp .AND. &
2291	.NOT. building_edge_w ) lad_local = lad_s(kk,j,i-1)
2292	!
2293	!-- In order to avoid that LAD is mistakenly considered at right-
2294	!-- sided building edges (here the topography-top index for the
2295	!-- u-component at index j,i is still on the building while the
2296	!-- topography top for the scalar isn't), LAD is taken from grid
2297	!-- point (j,i-1).
2298	IF ( lad_local > 0.0_wp .AND. lad_s(kk,j,i-1) == 0.0_wp .AND. &
2299	building_edge_e ) lad_local = lad_s(kk,j,i-1)
2300
2301	pre_tend = 0.0_wp
2302	pre_u = 0.0_wp
2303	!
2304	!-- Calculate preliminary value (pre_tend) of the tendency
2305	pre_tend = - canopy_drag_coeff * &
2306	lad_local * &
2307	SQRT( u(k,j,i)**2 + &
2308	( 0.25_wp * ( v(k,j,i-1) + &
2309	v(k,j,i) + &
2310	v(k,j+1,i) + &
2311	v(k,j+1,i-1) ) &
2312	)**2 + &
2313	( 0.25_wp * ( w(k-1,j,i-1) + &
2314	w(k-1,j,i) + &
2315	w(k,j,i-1) + &
2316	w(k,j,i) ) &
2317	)**2 &
2318	) * &
2319	u(k,j,i)
2320
2321	!
2322	!-- Calculate preliminary new velocity, based on pre_tend
2323	pre_u = u(k,j,i) + dt_3d * pre_tend
2324	!
2325	!-- Compare sign of old velocity and new preliminary velocity,
2326	!-- and in case the signs are different, limit the tendency
2327	IF ( SIGN(pre_u,u(k,j,i)) /= pre_u ) THEN
2328	pre_tend = - u(k,j,i) * ddt_3d
2329	ELSE
2330	pre_tend = pre_tend
2331	ENDIF
2332	!
2333	!-- Calculate final tendency
2334	tend(k,j,i) = tend(k,j,i) + pre_tend
2335	ENDDO
2336
2337
2338	!
2339	!-- v-component
2340	CASE ( 2 )
2341	!
2342	!-- Set control flags indicating north- and southward-orientated
2343	!-- building edges. Note, building_egde_s is set from the perspective
2344	!-- of the potential rooftop grid point, while building_edge_n is
2345	!-- is set from the perspective of the non-building grid point.
2346	building_edge_s = ANY( BTEST( wall_flags_total_0(:,j,i), 6 ) ) .AND. &
2347	.NOT. ANY( BTEST( wall_flags_total_0(:,j-1,i), 6 ) )
2348	building_edge_n = ANY( BTEST( wall_flags_total_0(:,j-1,i), 6 ) ) .AND. &
2349	.NOT. ANY( BTEST( wall_flags_total_0(:,j,i), 6 ) )
2350	!
2351	!-- Determine topography-top index on v-grid
2352	DO k = topo_top_ind(j,i,2) + 1, topo_top_ind(j,i,2) + pch_index_ji(j,i)
2353
2354	kk = k - topo_top_ind(j,i,2) !- lad arrays are defined flat
2355	!
2356	!-- In order to create sharp boundaries of the plant canopy,
2357	!-- the lad on the v-grid at index (k,j,i) is equal to lad_s(k,j,i),
2358	!-- rather than being interpolated from the surrounding lad_s,
2359	!-- because this would yield smaller lad at the canopy boundaries
2360	!-- than inside of the canopy.
2361	!-- For the same reason, the lad at the northmost(j+1)canopy
2362	!-- boundary on the v-grid equals lad_s(k,j,i), which is considered
2363	!-- in the next if-statement. Note, at left-sided building edges
2364	!-- this is not applied, here the LAD is equals the LAD at grid
2365	!-- point (k,j,i), in order to avoid that LAD is mistakenly mapped
2366	!-- on top of a roof where (usually) is no LAD is defined.
2367	lad_local = lad_s(kk,j,i)
2368	IF ( lad_local == 0.0_wp .AND. lad_s(kk,j-1,i) > 0.0_wp .AND. &
2369	.NOT. building_edge_s ) lad_local = lad_s(kk,j-1,i)
2370	!
2371	!-- In order to avoid that LAD is mistakenly considered at right-
2372	!-- sided building edges (here the topography-top index for the
2373	!-- u-component at index j,i is still on the building while the
2374	!-- topography top for the scalar isn't), LAD is taken from grid
2375	!-- point (j,i-1).
2376	IF ( lad_local > 0.0_wp .AND. lad_s(kk,j-1,i) == 0.0_wp .AND. &
2377	building_edge_n ) lad_local = lad_s(kk,j-1,i)
2378
2379	pre_tend = 0.0_wp
2380	pre_v = 0.0_wp
2381	!
2382	!-- Calculate preliminary value (pre_tend) of the tendency
2383	pre_tend = - canopy_drag_coeff * &
2384	lad_local * &
2385	SQRT( ( 0.25_wp * ( u(k,j-1,i) + &
2386	u(k,j-1,i+1) + &
2387	u(k,j,i) + &
2388	u(k,j,i+1) ) &
2389	)**2 + &
2390	v(k,j,i)**2 + &
2391	( 0.25_wp * ( w(k-1,j-1,i) + &
2392	w(k-1,j,i) + &
2393	w(k,j-1,i) + &
2394	w(k,j,i) ) &
2395	)**2 &
2396	) * &
2397	v(k,j,i)
2398
2399	!
2400	!-- Calculate preliminary new velocity, based on pre_tend
2401	pre_v = v(k,j,i) + dt_3d * pre_tend
2402	!
2403	!-- Compare sign of old velocity and new preliminary velocity,
2404	!-- and in case the signs are different, limit the tendency
2405	IF ( SIGN(pre_v,v(k,j,i)) /= pre_v ) THEN
2406	pre_tend = - v(k,j,i) * ddt_3d
2407	ELSE
2408	pre_tend = pre_tend
2409	ENDIF
2410	!
2411	!-- Calculate final tendency
2412	tend(k,j,i) = tend(k,j,i) + pre_tend
2413	ENDDO
2414
2415
2416	!
2417	!-- w-component
2418	CASE ( 3 )
2419	!
2420	!-- Determine topography-top index on w-grid
2421	DO k = topo_top_ind(j,i,3) + 1, topo_top_ind(j,i,3) + pch_index_ji(j,i) - 1
2422
2423	kk = k - topo_top_ind(j,i,3) !- lad arrays are defined flat
2424
2425	pre_tend = 0.0_wp
2426	pre_w = 0.0_wp
2427	!
2428	!-- Calculate preliminary value (pre_tend) of the tendency
2429	pre_tend = - canopy_drag_coeff * &
2430	(0.5_wp * &
2431	( lad_s(kk+1,j,i) + lad_s(kk,j,i) )) * &
2432	SQRT( ( 0.25_wp * ( u(k,j,i) + &
2433	u(k,j,i+1) + &
2434	u(k+1,j,i) + &
2435	u(k+1,j,i+1) ) &
2436	)**2 + &
2437	( 0.25_wp * ( v(k,j,i) + &
2438	v(k,j+1,i) + &
2439	v(k+1,j,i) + &
2440	v(k+1,j+1,i) ) &
2441	)**2 + &
2442	w(k,j,i)**2 &
2443	) * &
2444	w(k,j,i)
2445	!
2446	!-- Calculate preliminary new velocity, based on pre_tend
2447	pre_w = w(k,j,i) + dt_3d * pre_tend
2448	!
2449	!-- Compare sign of old velocity and new preliminary velocity,
2450	!-- and in case the signs are different, limit the tendency
2451	IF ( SIGN(pre_w,w(k,j,i)) /= pre_w ) THEN
2452	pre_tend = - w(k,j,i) * ddt_3d
2453	ELSE
2454	pre_tend = pre_tend
2455	ENDIF
2456	!
2457	!-- Calculate final tendency
2458	tend(k,j,i) = tend(k,j,i) + pre_tend
2459	ENDDO
2460
2461	!
2462	!-- potential temperature
2463	CASE ( 4 )
2464	!
2465	!-- Determine topography-top index on scalar grid
2466	IF ( humidity ) THEN
2467	DO k = topo_top_ind(j,i,0) + 1, topo_top_ind(j,i,0) + pch_index_ji(j,i)
2468	kk = k - topo_top_ind(j,i,0) !- lad arrays are defined flat
2469	tend(k,j,i) = tend(k,j,i) + pcm_heating_rate(kk,j,i) - &
2470	pcm_latent_rate(kk,j,i)
2471	ENDDO
2472	ELSE
2473	DO k = topo_top_ind(j,i,0) + 1, topo_top_ind(j,i,0) + pch_index_ji(j,i)
2474	kk = k - topo_top_ind(j,i,0) !- lad arrays are defined flat
2475	tend(k,j,i) = tend(k,j,i) + pcm_heating_rate(kk,j,i)
2476	ENDDO
2477	ENDIF
2478
2479	!
2480	!-- humidity
2481	CASE ( 5 )
2482	!
2483	!-- Determine topography-top index on scalar grid
2484	DO k = topo_top_ind(j,i,0) + 1, topo_top_ind(j,i,0) + pch_index_ji(j,i)
2485	kk = k - topo_top_ind(j,i,0) !- lad arrays are defined flat
2486	IF ( .NOT. plant_canopy_transpiration ) THEN
2487	! pcm_transpiration_rate is calculated in radiation model
2488	! in case of plant_canopy_transpiration = .T.
2489	! to include also the dependecy to the radiation
2490	! in the plant canopy box
2491	pcm_transpiration_rate(kk,j,i) = - leaf_scalar_exch_coeff &
2492	* lad_s(kk,j,i) * &
2493	SQRT( ( 0.5_wp * ( u(k,j,i) + &
2494	u(k,j,i+1) ) &
2495	)**2 + &
2496	( 0.5_wp * ( v(k,j,i) + &
2497	v(k,j+1,i) ) &
2498	)**2 + &
2499	( 0.5_wp * ( w(k-1,j,i) + &
2500	w(k,j,i) ) &
2501	)**2 &
2502	) * &
2503	( q(k,j,i) - leaf_surface_conc )
2504	ENDIF
2505
2506	tend(k,j,i) = tend(k,j,i) + pcm_transpiration_rate(kk,j,i)
2507
2508	ENDDO
2509
2510	!
2511	!-- sgs-tke
2512	CASE ( 6 )
2513	!
2514	!-- Determine topography-top index on scalar grid
2515	DO k = topo_top_ind(j,i,0) + 1, topo_top_ind(j,i,0) + pch_index_ji(j,i)
2516
2517	kk = k - topo_top_ind(j,i,0)
2518	tend(k,j,i) = tend(k,j,i) - &
2519	2.0_wp * canopy_drag_coeff * &
2520	lad_s(kk,j,i) * &
2521	SQRT( ( 0.5_wp * ( u(k,j,i) + &
2522	u(k,j,i+1) ) &
2523	)**2 + &
2524	( 0.5_wp * ( v(k,j,i) + &
2525	v(k,j+1,i) ) &
2526	)**2 + &
2527	( 0.5_wp * ( w(k,j,i) + &
2528	w(k+1,j,i) ) &
2529	)**2 &
2530	) * &
2531	e(k,j,i)
2532	ENDDO
2533
2534	!
2535	!-- scalar concentration
2536	CASE ( 7 )
2537	!
2538	!-- Determine topography-top index on scalar grid
2539	DO k = topo_top_ind(j,i,0) + 1, topo_top_ind(j,i,0) + pch_index_ji(j,i)
2540
2541	kk = k - topo_top_ind(j,i,0)
2542	tend(k,j,i) = tend(k,j,i) - &
2543	leaf_scalar_exch_coeff * &
2544	lad_s(kk,j,i) * &
2545	SQRT( ( 0.5_wp * ( u(k,j,i) + &
2546	u(k,j,i+1) ) &
2547	)**2 + &
2548	( 0.5_wp * ( v(k,j,i) + &
2549	v(k,j+1,i) ) &
2550	)**2 + &
2551	( 0.5_wp * ( w(k-1,j,i) + &
2552	w(k,j,i) ) &
2553	)**2 &
2554	) * &
2555	( s(k,j,i) - leaf_surface_conc )
2556	ENDDO
2557
2558	CASE DEFAULT
2559
2560	WRITE( message_string, * ) 'wrong component: ', component
2561	CALL message( 'pcm_tendency', 'PA0279', 1, 2, 0, 6, 0 )
2562
2563	END SELECT
2564
2565	END SUBROUTINE pcm_tendency_ij
2566
2567	!------------------------------------------------------------------------------!
2568	! Description:
2569	! ------------
2570	!> Subroutine writes global restart data
2571	!------------------------------------------------------------------------------!
2572	SUBROUTINE pcm_wrd_global
2573
2574	IF ( TRIM( restart_data_format_output ) == 'fortran_binary' ) THEN
2575
2576	CALL wrd_write_string( 'pch_index' )
2577	WRITE ( 14 ) pch_index
2578
2579	ELSEIF ( TRIM( restart_data_format_output ) == 'mpi' ) THEN
2580
2581	CALL wrd_mpi_io( 'pch_index', pch_index )
2582
2583	ENDIF
2584
2585	END SUBROUTINE pcm_wrd_global
2586
2587	!------------------------------------------------------------------------------!
2588	! Description:
2589	! ------------
2590	!> Subroutine writes local (subdomain) restart data
2591	!------------------------------------------------------------------------------!
2592	SUBROUTINE pcm_wrd_local
2593
2594	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: tmp_3d !< temporary array to store pcm data with
2595	!< non-standard vertical index bounds
2596
2597	IF ( TRIM( restart_data_format_output ) == 'fortran_binary' ) THEN
2598
2599	IF ( ALLOCATED( pcm_heatrate_av ) ) THEN
2600	CALL wrd_write_string( 'pcm_heatrate_av' )
2601	WRITE ( 14 ) pcm_heatrate_av
2602	ENDIF
2603
2604	IF ( ALLOCATED( pcm_latentrate_av ) ) THEN
2605	CALL wrd_write_string( 'pcm_latentrate_av' )
2606	WRITE ( 14 ) pcm_latentrate_av
2607	ENDIF
2608
2609	IF ( ALLOCATED( pcm_transpirationrate_av ) ) THEN
2610	CALL wrd_write_string( 'pcm_transpirationrate_av' )
2611	WRITE ( 14 ) pcm_transpirationrate_av
2612	ENDIF
2613
2614	ELSEIF ( TRIM( restart_data_format_output ) == 'mpi' ) THEN
2615
2616	!
2617	!-- Plant canopy arrays have non standard reduced vertical index bounds. They are stored with
2618	!-- full vertical bounds (bzb:nzt+1) in the restart file and must be re-stored before writing.
2619	IF ( ALLOCATED( pcm_heatrate_av ) ) THEN
2620	ALLOCATE( tmp_3d(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2621	tmp_3d(nzb:pch_index,:,:) = pcm_heatrate_av
2622	tmp_3d(pch_index+1:nzt+1,:,:) = 0.0_wp
2623	CALL wrd_mpi_io( 'pcm_heatrate_av', tmp_3d )
2624	DEALLOCATE( tmp_3d )
2625	ENDIF
2626	IF ( ALLOCATED( pcm_latentrate_av ) ) THEN
2627	ALLOCATE( tmp_3d(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2628	tmp_3d(nzb:pch_index,:,:) = pcm_latentrate_av
2629	tmp_3d(pch_index+1:nzt+1,:,:) = 0.0_wp
2630	CALL wrd_mpi_io( 'pcm_latentrate_av', tmp_3d )
2631	DEALLOCATE( tmp_3d )
2632	ENDIF
2633	IF ( ALLOCATED( pcm_transpirationrate_av ) ) THEN
2634	ALLOCATE( tmp_3d(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2635	tmp_3d(nzb:pch_index,:,:) = pcm_transpirationrate_av
2636	tmp_3d(pch_index+1:nzt+1,:,:) = 0.0_wp
2637	CALL wrd_mpi_io( 'pcm_transpirationrate_av', tmp_3d )
2638	DEALLOCATE( tmp_3d )
2639	ENDIF
2640
2641	ENDIF
2642
2643	END SUBROUTINE pcm_wrd_local
2644
2645
2646	END MODULE plant_canopy_model_mod

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |