Home

Context Navigation

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

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