Home

Context Navigation

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

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