Home

Context Navigation

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

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