Home

Context Navigation

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

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