Home

Context Navigation

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

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