Home

Context Navigation

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

Last change on this file since 4514 was 4514, checked in by suehring, 4 years ago
Bugfix in plant-canopy model for output of averaged transpiration rate after a restart; Revise check for output for plant heating rate and rename error message number; Surface-data output: enable output of mixing ratio and passive scalar concentration at the surface; Surface-data input: Add possibility to prescribe surface sensible and latent heat fluxes via static input file
Property svn:keywords set to `Id`
File size: 114.1 KB

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