Home

Context Navigation

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

Last change on this file since 4830 was 4828, checked in by Giersch, 3 years ago
Copyright updated to year 2021, interface pmc_sort removed to accelarate the nesting code
Property svn:keywords set to `Id`
File size: 129.2 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |