Home

Context Navigation

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

Last change on this file since 4526 was 4525, checked in by raasch, 5 years ago
salsa: added reading/writing of global restart data + reading/writing restart data with MPI-IO, bugfix for MPI-IO in plant_canopy_model_mod, tutorial user-interface dispersion_eularian_and_lpm_extended updated
Property svn:keywords set to `Id`
File size: 117.7 KB

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