Home

Context Navigation

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

Last change on this file since 4515 was 4515, checked in by suehring, 4 years ago
Rename error number again since this was recently given in -r 4511
Property svn:keywords set to `Id`
File size: 114.2 KB

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