Home

Context Navigation

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

Last change on this file since 2062 was 2025, checked in by kanani, 8 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 59.4 KB

Rev	Line
[1826]	1	!> @file plant_canopy_model_mod.f90
[2000]	2	!------------------------------------------------------------------------------!
[1036]	3	! This file is part of PALM.
	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	!
[1818]	17	! Copyright 1997-2016 Leibniz Universitaet Hannover
[2000]	18	!------------------------------------------------------------------------------!
[1036]	19	!
[257]	20	! Current revisions:
[138]	21	! -----------------
[2008]	22	!
[2025]	23	!
[2008]	24	! Former revisions:
	25	! -----------------
	26	! $Id: plant_canopy_model_mod.f90 2025 2016-10-12 16:44:03Z maronga $
	27	!
[2025]	28	! 2024 2016-10-12 16:42:37Z kanani
	29	! Added missing lad_s initialization
	30	!
[2012]	31	! 2011 2016-09-19 17:29:57Z kanani
	32	! Renamed canopy_heat_flux to pc_heating_rate, since the original meaning/
	33	! calculation of the quantity has changed, related to the urban surface model
	34	! and similar future applications.
	35	!
[2008]	36	! 2007 2016-08-24 15:47:17Z kanani
[2007]	37	! Added SUBROUTINE pcm_read_plant_canopy_3d for reading 3d plant canopy data
	38	! from file (new case canopy_mode=read_from_file_3d) in the course of
	39	! introduction of urban surface model,
	40	! introduced variable ext_coef,
	41	! resorted SUBROUTINEs to alphabetical order
[1827]	42	!
[2001]	43	! 2000 2016-08-20 18:09:15Z knoop
	44	! Forced header and separation lines into 80 columns
	45	!
[1961]	46	! 1960 2016-07-12 16:34:24Z suehring
	47	! Separate humidity and passive scalar
	48	!
[1954]	49	! 1953 2016-06-21 09:28:42Z suehring
	50	! Bugfix, lad_s and lad must be public
	51	!
[1827]	52	! 1826 2016-04-07 12:01:39Z maronga
	53	! Further modularization
	54	!
[1722]	55	! 1721 2015-11-16 12:56:48Z raasch
	56	! bugfixes: shf is reduced in areas covered with canopy only,
	57	! canopy is set on top of topography
	58	!
[1683]	59	! 1682 2015-10-07 23:56:08Z knoop
	60	! Code annotations made doxygen readable
	61	!
[1485]	62	! 1484 2014-10-21 10:53:05Z kanani
[1484]	63	! Changes due to new module structure of the plant canopy model:
	64	! module plant_canopy_model_mod now contains a subroutine for the
[1826]	65	! initialization of the canopy model (pcm_init),
[1484]	66	! limitation of the canopy drag (previously accounted for by calculation of
	67	! a limiting canopy timestep for the determination of the maximum LES timestep
	68	! in subroutine timestep) is now realized by the calculation of pre-tendencies
[1826]	69	! and preliminary velocities in subroutine pcm_tendency,
	70	! some redundant MPI communication removed in subroutine pcm_init
[1484]	71	! (was previously in init_3d_model),
	72	! unnecessary 3d-arrays lad_u, lad_v, lad_w removed - lad information on the
	73	! respective grid is now provided only by lad_s (e.g. in the calculation of
	74	! the tendency terms or of cum_lai_hf),
	75	! drag_coefficient, lai, leaf_surface_concentration,
	76	! scalar_exchange_coefficient, sec and sls renamed to canopy_drag_coeff,
	77	! cum_lai_hf, leaf_surface_conc, leaf_scalar_exch_coeff, lsec and lsc,
	78	! respectively,
	79	! unnecessary 3d-arrays cdc, lsc and lsec now defined as single-value constants,
	80	! USE-statements and ONLY-lists modified accordingly
[1341]	81	!
	82	! 1340 2014-03-25 19:45:13Z kanani
	83	! REAL constants defined as wp-kind
	84	!
[1321]	85	! 1320 2014-03-20 08:40:49Z raasch
[1320]	86	! ONLY-attribute added to USE-statements,
	87	! kind-parameters added to all INTEGER and REAL declaration statements,
	88	! kinds are defined in new module kinds,
	89	! old module precision_kind is removed,
	90	! revision history before 2012 removed,
	91	! comment fields (!:) to be used for variable explanations added to
	92	! all variable declaration statements
[153]	93	!
[1037]	94	! 1036 2012-10-22 13:43:42Z raasch
	95	! code put under GPL (PALM 3.9)
	96	!
[139]	97	! 138 2007-11-28 10:03:58Z letzel
	98	! Initial revision
	99	!
[138]	100	! Description:
	101	! ------------
[1682]	102	!> 1) Initialization of the canopy model, e.g. construction of leaf area density
[1826]	103	!> profile (subroutine pcm_init).
[1682]	104	!> 2) Calculation of sinks and sources of momentum, heat and scalar concentration
[1826]	105	!> due to canopy elements (subroutine pcm_tendency).
[138]	106	!------------------------------------------------------------------------------!
[1682]	107	MODULE plant_canopy_model_mod
	108
[1484]	109	USE arrays_3d, &
[1960]	110	ONLY: dzu, dzw, e, q, s, shf, tend, u, v, w, zu, zw
[138]	111
[1484]	112	USE indices, &
	113	ONLY: nbgp, nxl, nxlg, nxlu, nxr, nxrg, nyn, nyng, nys, nysg, nysv, &
	114	nz, nzb, nzb_s_inner, nzb_u_inner, nzb_v_inner, nzb_w_inner, nzt
	115
	116	USE kinds
	117
	118
	119	IMPLICIT NONE
	120
	121
[1682]	122	CHARACTER (LEN=20) :: canopy_mode = 'block' !< canopy coverage
[1484]	123
[1682]	124	INTEGER(iwp) :: pch_index = 0 !< plant canopy height/top index
[1484]	125	INTEGER(iwp) :: &
[1682]	126	lad_vertical_gradient_level_ind(10) = -9999 !< lad-profile levels (index)
[1484]	127
[1682]	128	LOGICAL :: calc_beta_lad_profile = .FALSE. !< switch for calc. of lad from beta func.
	129	LOGICAL :: plant_canopy = .FALSE. !< switch for use of canopy model
[1484]	130
[1682]	131	REAL(wp) :: alpha_lad = 9999999.9_wp !< coefficient for lad calculation
	132	REAL(wp) :: beta_lad = 9999999.9_wp !< coefficient for lad calculation
	133	REAL(wp) :: canopy_drag_coeff = 0.0_wp !< canopy drag coefficient (parameter)
	134	REAL(wp) :: cdc = 0.0_wp !< canopy drag coeff. (abbreviation used in equations)
	135	REAL(wp) :: cthf = 0.0_wp !< canopy top heat flux
	136	REAL(wp) :: dt_plant_canopy = 0.0_wp !< timestep account. for canopy drag
[2007]	137	REAL(wp) :: ext_coef = 0.6_wp !< extinction coefficient
[1682]	138	REAL(wp) :: lad_surface = 0.0_wp !< lad surface value
	139	REAL(wp) :: lai_beta = 0.0_wp !< leaf area index (lai) for lad calc.
[1484]	140	REAL(wp) :: &
[1682]	141	leaf_scalar_exch_coeff = 0.0_wp !< canopy scalar exchange coeff.
[1484]	142	REAL(wp) :: &
[1682]	143	leaf_surface_conc = 0.0_wp !< leaf surface concentration
	144	REAL(wp) :: lsec = 0.0_wp !< leaf scalar exchange coeff.
	145	REAL(wp) :: lsc = 0.0_wp !< leaf surface concentration
[1484]	146
	147	REAL(wp) :: &
[1682]	148	lad_vertical_gradient(10) = 0.0_wp !< lad gradient
[1484]	149	REAL(wp) :: &
[1682]	150	lad_vertical_gradient_level(10) = -9999999.9_wp !< lad-prof. levels (in m)
[1484]	151
[1682]	152	REAL(wp), DIMENSION(:), ALLOCATABLE :: lad !< leaf area density
	153	REAL(wp), DIMENSION(:), ALLOCATABLE :: pre_lad !< preliminary lad
[1484]	154
	155	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
[2011]	156	pc_heating_rate !< plant canopy heating rate
[1682]	157	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: cum_lai_hf !< cumulative lai for heatflux calc.
	158	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: lad_s !< lad on scalar-grid
[1484]	159
	160
	161	SAVE
	162
	163
[138]	164	PRIVATE
[1826]	165
	166	!
	167	!-- Public functions
	168	PUBLIC pcm_check_parameters, pcm_header, pcm_init, pcm_parin, pcm_tendency
[138]	169
[1826]	170	!
	171	!-- Public variables and constants
[2011]	172	PUBLIC pc_heating_rate, canopy_mode, cthf, dt_plant_canopy, lad, lad_s, &
[2007]	173	pch_index, plant_canopy
	174
[1484]	175
	176
[1826]	177	INTERFACE pcm_check_parameters
	178	MODULE PROCEDURE pcm_check_parameters
	179	END INTERFACE pcm_check_parameters
	180
	181	INTERFACE pcm_header
	182	MODULE PROCEDURE pcm_header
	183	END INTERFACE pcm_header
	184
	185	INTERFACE pcm_init
	186	MODULE PROCEDURE pcm_init
	187	END INTERFACE pcm_init
[138]	188
[1826]	189	INTERFACE pcm_parin
	190	MODULE PROCEDURE pcm_parin
[2007]	191	END INTERFACE pcm_parin
	192
	193	INTERFACE pcm_read_plant_canopy_3d
	194	MODULE PROCEDURE pcm_read_plant_canopy_3d
	195	END INTERFACE pcm_read_plant_canopy_3d
[1826]	196
	197	INTERFACE pcm_tendency
	198	MODULE PROCEDURE pcm_tendency
	199	MODULE PROCEDURE pcm_tendency_ij
	200	END INTERFACE pcm_tendency
[1484]	201
	202
[138]	203	CONTAINS
	204
[1826]	205
	206	!------------------------------------------------------------------------------!
	207	! Description:
	208	! ------------
	209	!> Check parameters routine for plant canopy model
	210	!------------------------------------------------------------------------------!
	211	SUBROUTINE pcm_check_parameters
[138]	212
[1826]	213	USE control_parameters, &
	214	ONLY: cloud_physics, message_string, microphysics_seifert
	215
	216
	217	IMPLICIT NONE
	218
	219
	220	IF ( canopy_drag_coeff == 0.0_wp ) THEN
	221	message_string = 'plant_canopy = .TRUE. requires a non-zero drag '// &
	222	'coefficient & given value is canopy_drag_coeff = 0.0'
	223	CALL message( 'check_parameters', 'PA0041', 1, 2, 0, 6, 0 )
	224	ENDIF
	225
	226	IF ( ( alpha_lad /= 9999999.9_wp .AND. beta_lad == 9999999.9_wp ) .OR.&
	227	beta_lad /= 9999999.9_wp .AND. alpha_lad == 9999999.9_wp ) THEN
	228	message_string = 'using the beta function for the construction ' // &
	229	'of the leaf area density profile requires ' // &
	230	'both alpha_lad and beta_lad to be /= 9999999.9'
	231	CALL message( 'check_parameters', 'PA0118', 1, 2, 0, 6, 0 )
	232	ENDIF
	233
	234	IF ( calc_beta_lad_profile .AND. lai_beta == 0.0_wp ) THEN
	235	message_string = 'using the beta function for the construction ' // &
	236	'of the leaf area density profile requires ' // &
	237	'a non-zero lai_beta, but given value is ' // &
	238	'lai_beta = 0.0'
	239	CALL message( 'check_parameters', 'PA0119', 1, 2, 0, 6, 0 )
	240	ENDIF
	241
	242	IF ( calc_beta_lad_profile .AND. lad_surface /= 0.0_wp ) THEN
	243	message_string = 'simultaneous setting of alpha_lad /= 9999999.9' // &
	244	'and lad_surface /= 0.0 is not possible, ' // &
	245	'use either vertical gradients or the beta ' // &
	246	'function for the construction of the leaf area '// &
	247	'density profile'
	248	CALL message( 'check_parameters', 'PA0120', 1, 2, 0, 6, 0 )
	249	ENDIF
	250
	251	IF ( cloud_physics .AND. microphysics_seifert ) THEN
	252	message_string = 'plant_canopy = .TRUE. requires cloud_scheme /=' // &
	253	' seifert_beheng'
	254	CALL message( 'check_parameters', 'PA0360', 1, 2, 0, 6, 0 )
	255	ENDIF
	256
	257
	258	END SUBROUTINE pcm_check_parameters
	259
	260
[138]	261	!------------------------------------------------------------------------------!
[1484]	262	! Description:
	263	! ------------
[1826]	264	!> Header output for plant canopy model
	265	!------------------------------------------------------------------------------!
	266	SUBROUTINE pcm_header ( io )
	267
	268	USE control_parameters, &
	269	ONLY: dz, passive_scalar
	270
	271
	272	IMPLICIT NONE
	273
	274	CHARACTER (LEN=10) :: coor_chr !<
	275
	276	CHARACTER (LEN=86) :: coordinates !<
	277	CHARACTER (LEN=86) :: gradients !<
	278	CHARACTER (LEN=86) :: leaf_area_density !<
	279	CHARACTER (LEN=86) :: slices !<
	280
	281	INTEGER(iwp) :: i !<
	282	INTEGER(iwp), INTENT(IN) :: io !< Unit of the output file
	283	INTEGER(iwp) :: k !<
	284
	285	REAL(wp) :: canopy_height !< canopy height (in m)
	286
	287	canopy_height = pch_index * dz
	288
	289	WRITE ( io, 1 ) canopy_mode, canopy_height, pch_index, &
	290	canopy_drag_coeff
	291	IF ( passive_scalar ) THEN
	292	WRITE ( io, 2 ) leaf_scalar_exch_coeff, &
	293	leaf_surface_conc
	294	ENDIF
	295
	296	!
	297	!-- Heat flux at the top of vegetation
	298	WRITE ( io, 3 ) cthf
	299
	300	!
	301	!-- Leaf area density profile, calculated either from given vertical
	302	!-- gradients or from beta probability density function.
	303	IF ( .NOT. calc_beta_lad_profile ) THEN
	304
	305	!-- Building output strings, starting with surface value
	306	WRITE ( leaf_area_density, '(F7.4)' ) lad_surface
	307	gradients = '------'
	308	slices = ' 0'
	309	coordinates = ' 0.0'
	310	i = 1
	311	DO WHILE ( i < 11 .AND. lad_vertical_gradient_level_ind(i) &
	312	/= -9999 )
	313
	314	WRITE (coor_chr,'(F7.2)') lad(lad_vertical_gradient_level_ind(i))
	315	leaf_area_density = TRIM( leaf_area_density ) // ' ' // &
	316	TRIM( coor_chr )
	317
	318	WRITE (coor_chr,'(F7.2)') lad_vertical_gradient(i)
	319	gradients = TRIM( gradients ) // ' ' // TRIM( coor_chr )
	320
	321	WRITE (coor_chr,'(I7)') lad_vertical_gradient_level_ind(i)
	322	slices = TRIM( slices ) // ' ' // TRIM( coor_chr )
	323
	324	WRITE (coor_chr,'(F7.1)') lad_vertical_gradient_level(i)
	325	coordinates = TRIM( coordinates ) // ' ' // TRIM( coor_chr )
	326
	327	i = i + 1
	328	ENDDO
	329
	330	WRITE ( io, 4 ) TRIM( coordinates ), TRIM( leaf_area_density ), &
	331	TRIM( gradients ), TRIM( slices )
	332
	333	ELSE
	334
	335	WRITE ( leaf_area_density, '(F7.4)' ) lad_surface
	336	coordinates = ' 0.0'
	337
	338	DO k = 1, pch_index
	339
	340	WRITE (coor_chr,'(F7.2)') lad(k)
	341	leaf_area_density = TRIM( leaf_area_density ) // ' ' // &
	342	TRIM( coor_chr )
	343
	344	WRITE (coor_chr,'(F7.1)') zu(k)
	345	coordinates = TRIM( coordinates ) // ' ' // TRIM( coor_chr )
	346
	347	ENDDO
	348
	349	WRITE ( io, 5 ) TRIM( coordinates ), TRIM( leaf_area_density ), &
	350	alpha_lad, beta_lad, lai_beta
	351
	352	ENDIF
	353
	354	1 FORMAT (//' Vegetation canopy (drag) model:'/ &
	355	' ------------------------------'// &
	356	' Canopy mode: ', A / &
	357	' Canopy height: ',F6.2,'m (',I4,' grid points)' / &
	358	' Leaf drag coefficient: ',F6.2 /)
	359	2 FORMAT (/ ' Scalar exchange coefficient: ',F6.2 / &
	360	' Scalar concentration at leaf surfaces in kg/m**3: ',F6.2 /)
	361	3 FORMAT (' Predefined constant heatflux at the top of the vegetation: ',F6.2, &
	362	' K m/s')
	363	4 FORMAT (/ ' Characteristic levels of the leaf area density:'// &
	364	' Height: ',A,' m'/ &
	365	' Leaf area density: ',A,' m2/m3'/ &
	366	' Gradient: ',A,' m2/m4'/ &
	367	' Gridpoint: ',A)
	368	5 FORMAT (//' Characteristic levels of the leaf area density and coefficients:'&
	369	// ' Height: ',A,' m'/ &
	370	' Leaf area density: ',A,' m2/m3'/ &
	371	' Coefficient alpha: ',F6.2 / &
	372	' Coefficient beta: ',F6.2 / &
	373	' Leaf area index: ',F6.2,' m2/m2' /)
	374
	375	END SUBROUTINE pcm_header
	376
	377
	378	!------------------------------------------------------------------------------!
	379	! Description:
	380	! ------------
[1682]	381	!> Initialization of the plant canopy model
[138]	382	!------------------------------------------------------------------------------!
[1826]	383	SUBROUTINE pcm_init
[1484]	384
	385
	386	USE control_parameters, &
[2007]	387	ONLY: coupling_char, dz, humidity, io_blocks, io_group, &
	388	message_string, ocean, passive_scalar
[1484]	389
	390
	391	IMPLICIT NONE
	392
[2007]	393	CHARACTER(10) :: pct
	394
	395	INTEGER(iwp) :: i !< running index
	396	INTEGER(iwp) :: ii !< index
	397	INTEGER(iwp) :: j !< running index
	398	INTEGER(iwp) :: k !< running index
[1484]	399
[2007]	400	REAL(wp) :: int_bpdf !< vertical integral for lad-profile construction
	401	REAL(wp) :: dzh !< vertical grid spacing in units of canopy height
	402	REAL(wp) :: gradient !< gradient for lad-profile construction
	403	REAL(wp) :: canopy_height !< canopy height for lad-profile construction
	404	REAL(wp) :: pcv(nzb:nzt+1) !<
	405
[1484]	406	!
	407	!-- Allocate one-dimensional arrays for the computation of the
	408	!-- leaf area density (lad) profile
	409	ALLOCATE( lad(0:nz+1), pre_lad(0:nz+1) )
	410	lad = 0.0_wp
	411	pre_lad = 0.0_wp
	412
	413	!
[1826]	414	!-- Set flag that indicates that the lad-profile shall be calculated by using
	415	!-- a beta probability density function
	416	IF ( alpha_lad /= 9999999.9_wp .AND. beta_lad /= 9999999.9_wp ) THEN
	417	calc_beta_lad_profile = .TRUE.
	418	ENDIF
	419
	420
	421	!
[1484]	422	!-- Compute the profile of leaf area density used in the plant
	423	!-- canopy model. The profile can either be constructed from
	424	!-- prescribed vertical gradients of the leaf area density or by
	425	!-- using a beta probability density function (see e.g. Markkanen et al.,
	426	!-- 2003: Boundary-Layer Meteorology, 106, 437-459)
	427	IF ( .NOT. calc_beta_lad_profile ) THEN
	428
	429	!
	430	!-- Use vertical gradients for lad-profile construction
	431	i = 1
	432	gradient = 0.0_wp
	433
	434	IF ( .NOT. ocean ) THEN
	435
	436	lad(0) = lad_surface
	437	lad_vertical_gradient_level_ind(1) = 0
	438
	439	DO k = 1, pch_index
	440	IF ( i < 11 ) THEN
	441	IF ( lad_vertical_gradient_level(i) < zu(k) .AND. &
	442	lad_vertical_gradient_level(i) >= 0.0_wp ) THEN
	443	gradient = lad_vertical_gradient(i)
	444	lad_vertical_gradient_level_ind(i) = k - 1
	445	i = i + 1
	446	ENDIF
	447	ENDIF
	448	IF ( gradient /= 0.0_wp ) THEN
	449	IF ( k /= 1 ) THEN
	450	lad(k) = lad(k-1) + dzu(k) * gradient
	451	ELSE
	452	lad(k) = lad_surface + dzu(k) * gradient
	453	ENDIF
	454	ELSE
	455	lad(k) = lad(k-1)
	456	ENDIF
	457	ENDDO
	458
	459	ENDIF
	460
	461	!
	462	!-- In case of no given leaf area density gradients, choose a vanishing
	463	!-- gradient. This information is used for the HEADER and the RUN_CONTROL
	464	!-- file.
	465	IF ( lad_vertical_gradient_level(1) == -9999999.9_wp ) THEN
	466	lad_vertical_gradient_level(1) = 0.0_wp
	467	ENDIF
	468
	469	ELSE
	470
	471	!
	472	!-- Use beta function for lad-profile construction
	473	int_bpdf = 0.0_wp
	474	canopy_height = pch_index * dz
	475
	476	DO k = nzb, pch_index
	477	int_bpdf = int_bpdf + &
[1826]	478	( ( ( zw(k) / canopy_height )*( alpha_lad-1.0_wp ) ) &
	479	( ( 1.0_wp - ( zw(k) / canopy_height ) )**( &
	480	beta_lad-1.0_wp ) ) &
	481	* ( ( zw(k+1)-zw(k) ) / canopy_height ) )
[1484]	482	ENDDO
	483
	484	!
	485	!-- Preliminary lad profile (defined on w-grid)
	486	DO k = nzb, pch_index
[1826]	487	pre_lad(k) = lai_beta * &
	488	( ( ( zw(k) / canopy_height )**( alpha_lad-1.0_wp ) ) &
	489	* ( ( 1.0_wp - ( zw(k) / canopy_height ) )**( &
	490	beta_lad-1.0_wp ) ) / int_bpdf &
	491	) / canopy_height
[1484]	492	ENDDO
	493
	494	!
	495	!-- Final lad profile (defined on scalar-grid level, since most prognostic
	496	!-- quantities are defined there, hence, less interpolation is required
	497	!-- when calculating the canopy tendencies)
	498	lad(0) = pre_lad(0)
	499	DO k = nzb+1, pch_index
	500	lad(k) = 0.5 * ( pre_lad(k-1) + pre_lad(k) )
	501	ENDDO
	502
	503	ENDIF
	504
	505	!
	506	!-- Allocate 3D-array for the leaf area density (lad_s). In case of a
	507	!-- prescribed canopy-top heat flux (cthf), allocate 3D-arrays for
	508	!-- the cumulative leaf area index (cum_lai_hf) and the canopy heat flux.
	509	ALLOCATE( lad_s(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
	510
	511	IF ( cthf /= 0.0_wp ) THEN
	512	ALLOCATE( cum_lai_hf(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
[2011]	513	pc_heating_rate(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
[1484]	514	ENDIF
	515
	516	!
	517	!-- Initialize canopy parameters cdc (canopy drag coefficient),
	518	!-- lsec (leaf scalar exchange coefficient), lsc (leaf surface concentration)
	519	!-- with the prescribed values
	520	cdc = canopy_drag_coeff
	521	lsec = leaf_scalar_exch_coeff
	522	lsc = leaf_surface_conc
	523
	524	!
	525	!-- Initialization of the canopy coverage in the model domain:
	526	!-- Setting the parameter canopy_mode = 'block' initializes a canopy, which
	527	!-- fully covers the domain surface
	528	SELECT CASE ( TRIM( canopy_mode ) )
	529
	530	CASE( 'block' )
	531
	532	DO i = nxlg, nxrg
	533	DO j = nysg, nyng
	534	lad_s(:,j,i) = lad(:)
	535	ENDDO
	536	ENDDO
	537
[2007]	538	CASE ( 'read_from_file_3d' )
	539	!
	540	!-- Initialize canopy parameters cdc (canopy drag coefficient),
	541	!-- lsec (leaf scalar exchange coefficient), lsc (leaf surface concentration)
	542	!-- from file which contains complete 3D data (separate vertical profiles for
	543	!-- each location).
	544	CALL pcm_read_plant_canopy_3d
	545
[1484]	546	CASE DEFAULT
	547	!
[2007]	548	!-- The DEFAULT case is reached either if the parameter
	549	!-- canopy mode contains a wrong character string or if the
	550	!-- user has coded a special case in the user interface.
	551	!-- There, the subroutine user_init_plant_canopy checks
	552	!-- which of these two conditions applies.
	553	CALL user_init_plant_canopy
[1484]	554
	555	END SELECT
	556
	557	!
[2011]	558	!-- Initialization of the canopy heat source distribution due to heating
	559	!-- of the canopy layers by incoming solar radiation, in case that a non-zero
	560	!-- value is set for the canopy top heat flux (cthf), which equals the
	561	!-- available net radiation at canopy top.
	562	!-- The heat source distribution is calculated by a decaying exponential
	563	!-- function of the downward cumulative leaf area index (cum_lai_hf),
	564	!-- assuming that the foliage inside the plant canopy is heated by solar
	565	!-- radiation penetrating the canopy layers according to the distribution
	566	!-- of net radiation as suggested by Brown & Covey (1966; Agric. Meteorol. 3,
	567	!-- 73â96). This approach has been applied e.g. by Shaw & Schumann (1992;
	568	!-- Bound.-Layer Meteorol. 61, 47â64)
[1484]	569	IF ( cthf /= 0.0_wp ) THEN
	570	!
[2011]	571	!-- Piecewise calculation of the cumulative leaf area index by vertical
[1484]	572	!-- integration of the leaf area density
	573	cum_lai_hf(:,:,:) = 0.0_wp
	574	DO i = nxlg, nxrg
	575	DO j = nysg, nyng
	576	DO k = pch_index-1, 0, -1
	577	IF ( k == pch_index-1 ) THEN
	578	cum_lai_hf(k,j,i) = cum_lai_hf(k+1,j,i) + &
	579	( 0.5_wp * lad_s(k+1,j,i) * &
	580	( zw(k+1) - zu(k+1) ) ) + &
	581	( 0.5_wp * ( 0.5_wp * ( lad_s(k+1,j,i) + &
	582	lad_s(k,j,i) ) + &
	583	lad_s(k+1,j,i) ) * &
	584	( zu(k+1) - zw(k) ) )
	585	ELSE
	586	cum_lai_hf(k,j,i) = cum_lai_hf(k+1,j,i) + &
	587	( 0.5_wp * ( 0.5_wp * ( lad_s(k+2,j,i) + &
	588	lad_s(k+1,j,i) ) + &
	589	lad_s(k+1,j,i) ) * &
	590	( zw(k+1) - zu(k+1) ) ) + &
	591	( 0.5_wp * ( 0.5_wp * ( lad_s(k+1,j,i) + &
	592	lad_s(k,j,i) ) + &
	593	lad_s(k+1,j,i) ) * &
	594	( zu(k+1) - zw(k) ) )
	595	ENDIF
	596	ENDDO
	597	ENDDO
	598	ENDDO
	599
	600	!
[2011]	601	!-- Calculation of the heating rate (K/s) within the different layers of
	602	!-- the plant canopy
[1484]	603	DO i = nxlg, nxrg
	604	DO j = nysg, nyng
	605	!
[2011]	606	!-- Calculation only necessary in areas covered with canopy
	607	IF ( cum_lai_hf(0,j,i) /= 0.0_wp ) THEN
	608	!--
	609	!-- In areas with canopy the surface value of the canopy heat
	610	!-- flux distribution overrides the surface heat flux (shf)
	611	shf(j,i) = cthf * exp( -ext_coef * cum_lai_hf(0,j,i) )
	612	!
	613	!-- Within the different canopy layers the plant-canopy heating
	614	!-- rate (pc_heating_rate) is calculated as the vertical
	615	!-- divergence of the canopy heat fluxes at the top and bottom
	616	!-- of the respective layer
	617	DO k = 1, pch_index
	618	pc_heating_rate(k,j,i) = cthf * &
	619	( exp(-ext_coef*cum_lai_hf(k,j,i)) - &
	620	exp(-ext_coef*cum_lai_hf(k-1,j,i)) ) / dzw(k)
	621	ENDDO
[1721]	622	ENDIF
	623	ENDDO
	624	ENDDO
[1484]	625
	626	ENDIF
	627
	628
	629
[1826]	630	END SUBROUTINE pcm_init
[1484]	631
	632
[2007]	633	!------------------------------------------------------------------------------!
	634	! Description:
	635	! ------------
	636	!> Parin for &canopy_par for plant canopy model
	637	!------------------------------------------------------------------------------!
	638	SUBROUTINE pcm_parin
[1484]	639
[2007]	640
	641	IMPLICIT NONE
	642
	643	CHARACTER (LEN=80) :: line !< dummy string that contains the current line of the parameter file
	644
	645	NAMELIST /canopy_par/ alpha_lad, beta_lad, canopy_drag_coeff, &
	646	canopy_mode, cthf, &
	647	lad_surface, &
	648	lad_vertical_gradient, &
	649	lad_vertical_gradient_level, &
	650	lai_beta, &
	651	leaf_scalar_exch_coeff, &
	652	leaf_surface_conc, pch_index
	653
	654	line = ' '
	655
	656	!
	657	!-- Try to find radiation model package
	658	REWIND ( 11 )
	659	line = ' '
	660	DO WHILE ( INDEX( line, '&canopy_par' ) == 0 )
	661	READ ( 11, '(A)', END=10 ) line
	662	ENDDO
	663	BACKSPACE ( 11 )
	664
	665	!
	666	!-- Read user-defined namelist
	667	READ ( 11, canopy_par )
	668
	669	!
	670	!-- Set flag that indicates that the radiation model is switched on
	671	plant_canopy = .TRUE.
	672
	673	10 CONTINUE
	674
	675
	676	END SUBROUTINE pcm_parin
	677
	678
	679
[1484]	680	!------------------------------------------------------------------------------!
	681	! Description:
	682	! ------------
[2007]	683	!
	684	!> Loads 3D plant canopy data from file. File format is as follows:
	685	!>
	686	!> num_levels
	687	!> dtype,x,y,value(nzb),value(nzb+1), ... ,value(nzb+num_levels-1)
	688	!> dtype,x,y,value(nzb),value(nzb+1), ... ,value(nzb+num_levels-1)
	689	!> dtype,x,y,value(nzb),value(nzb+1), ... ,value(nzb+num_levels-1)
	690	!> ...
	691	!>
	692	!> i.e. first line determines number of levels and further lines represent plant
	693	!> canopy data, one line per column and variable. In each data line,
	694	!> dtype represents variable to be set:
	695	!>
	696	!> dtype=1: leaf area density (lad_s)
	697	!> dtype=2: canopy drag coefficient (cdc)
	698	!> dtype=3: leaf scalar exchange coefficient (lsec)
	699	!> dtype=4: leaf surface concentration (lsc)
	700	!>
	701	!> Zeros are added automatically above num_levels until top of domain. Any
	702	!> non-specified (x,y) columns have zero values as default.
	703	!------------------------------------------------------------------------------!
	704	SUBROUTINE pcm_read_plant_canopy_3d
	705	USE control_parameters, &
	706	ONLY: passive_scalar, message_string
	707	IMPLICIT NONE
	708
	709	INTEGER(iwp) :: i, j, dtype, nzp, nzpltop, nzpl, kk
	710	REAL(wp), DIMENSION(:), ALLOCATABLE :: col
	711
[2024]	712	lad_s = 0.0_wp
[2007]	713	OPEN(152, file='PLANT_CANOPY_DATA_3D', access='SEQUENTIAL', &
	714	action='READ', status='OLD', form='FORMATTED', err=515)
	715	READ(152, *, err=516, end=517) nzp !< read first line = number of vertical layers
	716	ALLOCATE(col(1:nzp))
	717	nzpltop = MIN(nzt+1, nzb+nzp-1)
	718	nzpl = nzpltop - nzb + 1 !< no. of layers to assign
	719
	720	DO
	721	READ(152, *, err=516, end=517) dtype, i, j, col(:)
	722	IF ( i < nxlg .or. i > nxrg .or. j < nysg .or. j > nyng ) CYCLE
	723
	724	SELECT CASE (dtype)
	725	CASE( 1 ) !< leaf area density
	726	!-- only lad_s has flat z-coordinate, others have regular
	727	kk = nzb_s_inner(j, i)
	728	lad_s(nzb:nzpltop-kk, j, i) = col(1+kk:nzpl)
	729	! CASE( 2 ) !< canopy drag coefficient
	730	! cdc(nzb:nzpltop, j, i) = col(1:nzpl)
	731	! CASE( 3 ) !< leaf scalar exchange coefficient
	732	! lsec(nzb:nzpltop, j, i) = col(1:nzpl)
	733	! CASE( 4 ) !< leaf surface concentration
	734	! lsc(nzb:nzpltop, j, i) = col(1:nzpl)
	735	CASE DEFAULT
	736	write(message_string, '(a,i2,a)') &
	737	'Unknown record type in file PLANT_CANOPY_DATA_3D: "', dtype, '"'
	738	CALL message( 'pcm_read_plant_canopy_3d', 'PA0530', 1, 2, 0, 6, 0 )
	739	END SELECT
	740	ENDDO
	741
	742	515 message_string = 'error opening file PLANT_CANOPY_DATA_3D'
	743	CALL message( 'pcm_read_plant_canopy_3d', 'PA0531', 1, 2, 0, 6, 0 )
	744
	745	516 message_string = 'error reading file PLANT_CANOPY_DATA_3D'
	746	CALL message( 'pcm_read_plant_canopy_3d', 'PA0532', 1, 2, 0, 6, 0 )
	747
	748	517 CLOSE(152)
	749	DEALLOCATE(col)
	750
	751	END SUBROUTINE pcm_read_plant_canopy_3d
	752
	753
	754
	755	!------------------------------------------------------------------------------!
	756	! Description:
	757	! ------------
[1682]	758	!> Calculation of the tendency terms, accounting for the effect of the plant
	759	!> canopy on momentum and scalar quantities.
	760	!>
	761	!> The canopy is located where the leaf area density lad_s(k,j,i) > 0.0
[1826]	762	!> (defined on scalar grid), as initialized in subroutine pcm_init.
[1682]	763	!> The lad on the w-grid is vertically interpolated from the surrounding
	764	!> lad_s. The upper boundary of the canopy is defined on the w-grid at
	765	!> k = pch_index. Here, the lad is zero.
	766	!>
	767	!> The canopy drag must be limited (previously accounted for by calculation of
	768	!> a limiting canopy timestep for the determination of the maximum LES timestep
	769	!> in subroutine timestep), since it is physically impossible that the canopy
	770	!> drag alone can locally change the sign of a velocity component. This
	771	!> limitation is realized by calculating preliminary tendencies and velocities.
	772	!> It is subsequently checked if the preliminary new velocity has a different
	773	!> sign than the current velocity. If so, the tendency is limited in a way that
	774	!> the velocity can at maximum be reduced to zero by the canopy drag.
	775	!>
	776	!>
	777	!> Call for all grid points
[1484]	778	!------------------------------------------------------------------------------!
[1826]	779	SUBROUTINE pcm_tendency( component )
[138]	780
	781
[1320]	782	USE control_parameters, &
[1484]	783	ONLY: dt_3d, message_string
[1320]	784
	785	USE kinds
	786
[138]	787	IMPLICIT NONE
	788
[1682]	789	INTEGER(iwp) :: component !< prognostic variable (u,v,w,pt,q,e)
	790	INTEGER(iwp) :: i !< running index
	791	INTEGER(iwp) :: j !< running index
	792	INTEGER(iwp) :: k !< running index
[1721]	793	INTEGER(iwp) :: kk !< running index for flat lad arrays
[1484]	794
[1682]	795	REAL(wp) :: ddt_3d !< inverse of the LES timestep (dt_3d)
	796	REAL(wp) :: lad_local !< local lad value
	797	REAL(wp) :: pre_tend !< preliminary tendency
	798	REAL(wp) :: pre_u !< preliminary u-value
	799	REAL(wp) :: pre_v !< preliminary v-value
	800	REAL(wp) :: pre_w !< preliminary w-value
[1484]	801
	802
	803	ddt_3d = 1.0_wp / dt_3d
[138]	804
	805	!
[1484]	806	!-- Compute drag for the three velocity components and the SGS-TKE:
[138]	807	SELECT CASE ( component )
	808
	809	!
	810	!-- u-component
	811	CASE ( 1 )
	812	DO i = nxlu, nxr
	813	DO j = nys, nyn
[1721]	814	DO k = nzb_u_inner(j,i)+1, nzb_u_inner(j,i)+pch_index
[1484]	815
[1721]	816	kk = k - nzb_u_inner(j,i) !- lad arrays are defined flat
[1484]	817	!
	818	!-- In order to create sharp boundaries of the plant canopy,
	819	!-- the lad on the u-grid at index (k,j,i) is equal to
	820	!-- lad_s(k,j,i), rather than being interpolated from the
	821	!-- surrounding lad_s, because this would yield smaller lad
	822	!-- at the canopy boundaries than inside of the canopy.
	823	!-- For the same reason, the lad at the rightmost(i+1)canopy
	824	!-- boundary on the u-grid equals lad_s(k,j,i).
[1721]	825	lad_local = lad_s(kk,j,i)
	826	IF ( lad_local == 0.0_wp .AND. lad_s(kk,j,i-1) > 0.0_wp )&
	827	THEN
	828	lad_local = lad_s(kk,j,i-1)
[1484]	829	ENDIF
	830
	831	pre_tend = 0.0_wp
	832	pre_u = 0.0_wp
	833	!
	834	!-- Calculate preliminary value (pre_tend) of the tendency
	835	pre_tend = - cdc * &
	836	lad_local * &
	837	SQRT( u(k,j,i)**2 + &
	838	( 0.25_wp * ( v(k,j,i-1) + &
	839	v(k,j,i) + &
	840	v(k,j+1,i) + &
	841	v(k,j+1,i-1) ) &
	842	)**2 + &
	843	( 0.25_wp * ( w(k-1,j,i-1) + &
	844	w(k-1,j,i) + &
	845	w(k,j,i-1) + &
	846	w(k,j,i) ) &
	847	)**2 &
	848	) * &
	849	u(k,j,i)
	850
	851	!
	852	!-- Calculate preliminary new velocity, based on pre_tend
	853	pre_u = u(k,j,i) + dt_3d * pre_tend
	854	!
	855	!-- Compare sign of old velocity and new preliminary velocity,
	856	!-- and in case the signs are different, limit the tendency
	857	IF ( SIGN(pre_u,u(k,j,i)) /= pre_u ) THEN
	858	pre_tend = - u(k,j,i) * ddt_3d
	859	ELSE
	860	pre_tend = pre_tend
	861	ENDIF
	862	!
	863	!-- Calculate final tendency
	864	tend(k,j,i) = tend(k,j,i) + pre_tend
	865
[138]	866	ENDDO
	867	ENDDO
	868	ENDDO
	869
	870	!
	871	!-- v-component
	872	CASE ( 2 )
	873	DO i = nxl, nxr
	874	DO j = nysv, nyn
[1721]	875	DO k = nzb_v_inner(j,i)+1, nzb_v_inner(j,i)+pch_index
[1484]	876
[1721]	877	kk = k - nzb_v_inner(j,i) !- lad arrays are defined flat
[1484]	878	!
	879	!-- In order to create sharp boundaries of the plant canopy,
	880	!-- the lad on the v-grid at index (k,j,i) is equal to
	881	!-- lad_s(k,j,i), rather than being interpolated from the
	882	!-- surrounding lad_s, because this would yield smaller lad
	883	!-- at the canopy boundaries than inside of the canopy.
	884	!-- For the same reason, the lad at the northmost(j+1) canopy
	885	!-- boundary on the v-grid equals lad_s(k,j,i).
[1721]	886	lad_local = lad_s(kk,j,i)
	887	IF ( lad_local == 0.0_wp .AND. lad_s(kk,j-1,i) > 0.0_wp )&
	888	THEN
	889	lad_local = lad_s(kk,j-1,i)
[1484]	890	ENDIF
	891
	892	pre_tend = 0.0_wp
	893	pre_v = 0.0_wp
	894	!
	895	!-- Calculate preliminary value (pre_tend) of the tendency
	896	pre_tend = - cdc * &
	897	lad_local * &
	898	SQRT( ( 0.25_wp * ( u(k,j-1,i) + &
	899	u(k,j-1,i+1) + &
	900	u(k,j,i) + &
	901	u(k,j,i+1) ) &
	902	)**2 + &
	903	v(k,j,i)**2 + &
	904	( 0.25_wp * ( w(k-1,j-1,i) + &
	905	w(k-1,j,i) + &
	906	w(k,j-1,i) + &
	907	w(k,j,i) ) &
	908	)**2 &
	909	) * &
	910	v(k,j,i)
	911
	912	!
	913	!-- Calculate preliminary new velocity, based on pre_tend
	914	pre_v = v(k,j,i) + dt_3d * pre_tend
	915	!
	916	!-- Compare sign of old velocity and new preliminary velocity,
	917	!-- and in case the signs are different, limit the tendency
	918	IF ( SIGN(pre_v,v(k,j,i)) /= pre_v ) THEN
	919	pre_tend = - v(k,j,i) * ddt_3d
	920	ELSE
	921	pre_tend = pre_tend
	922	ENDIF
	923	!
	924	!-- Calculate final tendency
	925	tend(k,j,i) = tend(k,j,i) + pre_tend
	926
[138]	927	ENDDO
	928	ENDDO
	929	ENDDO
	930
	931	!
	932	!-- w-component
	933	CASE ( 3 )
	934	DO i = nxl, nxr
	935	DO j = nys, nyn
[1721]	936	DO k = nzb_w_inner(j,i)+1, nzb_w_inner(j,i)+pch_index-1
[1484]	937
[1721]	938	kk = k - nzb_w_inner(j,i) !- lad arrays are defined flat
	939
[1484]	940	pre_tend = 0.0_wp
	941	pre_w = 0.0_wp
	942	!
	943	!-- Calculate preliminary value (pre_tend) of the tendency
	944	pre_tend = - cdc * &
	945	(0.5_wp * &
[1721]	946	( lad_s(kk+1,j,i) + lad_s(kk,j,i) )) * &
[1484]	947	SQRT( ( 0.25_wp * ( u(k,j,i) + &
	948	u(k,j,i+1) + &
	949	u(k+1,j,i) + &
	950	u(k+1,j,i+1) ) &
	951	)**2 + &
	952	( 0.25_wp * ( v(k,j,i) + &
	953	v(k,j+1,i) + &
	954	v(k+1,j,i) + &
	955	v(k+1,j+1,i) ) &
	956	)**2 + &
	957	w(k,j,i)**2 &
	958	) * &
	959	w(k,j,i)
	960	!
	961	!-- Calculate preliminary new velocity, based on pre_tend
	962	pre_w = w(k,j,i) + dt_3d * pre_tend
	963	!
	964	!-- Compare sign of old velocity and new preliminary velocity,
	965	!-- and in case the signs are different, limit the tendency
	966	IF ( SIGN(pre_w,w(k,j,i)) /= pre_w ) THEN
	967	pre_tend = - w(k,j,i) * ddt_3d
	968	ELSE
	969	pre_tend = pre_tend
	970	ENDIF
	971	!
	972	!-- Calculate final tendency
	973	tend(k,j,i) = tend(k,j,i) + pre_tend
	974
[138]	975	ENDDO
	976	ENDDO
	977	ENDDO
	978
	979	!
[153]	980	!-- potential temperature
[138]	981	CASE ( 4 )
	982	DO i = nxl, nxr
	983	DO j = nys, nyn
[1721]	984	DO k = nzb_s_inner(j,i)+1, nzb_s_inner(j,i)+pch_index
	985	kk = k - nzb_s_inner(j,i) !- lad arrays are defined flat
[2011]	986	tend(k,j,i) = tend(k,j,i) + pc_heating_rate(kk,j,i)
[153]	987	ENDDO
	988	ENDDO
	989	ENDDO
	990
	991	!
[1960]	992	!-- humidity
[153]	993	CASE ( 5 )
	994	DO i = nxl, nxr
	995	DO j = nys, nyn
[1721]	996	DO k = nzb_s_inner(j,i)+1, nzb_s_inner(j,i)+pch_index
	997	kk = k - nzb_s_inner(j,i) !- lad arrays are defined flat
[1484]	998	tend(k,j,i) = tend(k,j,i) - &
	999	lsec * &
[1721]	1000	lad_s(kk,j,i) * &
[1484]	1001	SQRT( ( 0.5_wp * ( u(k,j,i) + &
	1002	u(k,j,i+1) ) &
	1003	)**2 + &
	1004	( 0.5_wp * ( v(k,j,i) + &
	1005	v(k,j+1,i) ) &
	1006	)**2 + &
	1007	( 0.5_wp * ( w(k-1,j,i) + &
	1008	w(k,j,i) ) &
	1009	)**2 &
	1010	) * &
	1011	( q(k,j,i) - lsc )
[153]	1012	ENDDO
	1013	ENDDO
	1014	ENDDO
	1015
	1016	!
	1017	!-- sgs-tke
	1018	CASE ( 6 )
	1019	DO i = nxl, nxr
	1020	DO j = nys, nyn
[1721]	1021	DO k = nzb_s_inner(j,i)+1, nzb_s_inner(j,i)+pch_index
	1022	kk = k - nzb_s_inner(j,i) !- lad arrays are defined flat
[1484]	1023	tend(k,j,i) = tend(k,j,i) - &
	1024	2.0_wp * cdc * &
[1721]	1025	lad_s(kk,j,i) * &
[1484]	1026	SQRT( ( 0.5_wp * ( u(k,j,i) + &
	1027	u(k,j,i+1) ) &
	1028	)**2 + &
	1029	( 0.5_wp * ( v(k,j,i) + &
	1030	v(k,j+1,i) ) &
	1031	)**2 + &
	1032	( 0.5_wp * ( w(k,j,i) + &
	1033	w(k+1,j,i) ) &
	1034	)**2 &
	1035	) * &
	1036	e(k,j,i)
[138]	1037	ENDDO
	1038	ENDDO
	1039	ENDDO
[1960]	1040	!
	1041	!-- scalar concentration
	1042	CASE ( 7 )
	1043	DO i = nxl, nxr
	1044	DO j = nys, nyn
	1045	DO k = nzb_s_inner(j,i)+1, nzb_s_inner(j,i)+pch_index
	1046	kk = k - nzb_s_inner(j,i) !- lad arrays are defined flat
	1047	tend(k,j,i) = tend(k,j,i) - &
	1048	lsec * &
	1049	lad_s(kk,j,i) * &
	1050	SQRT( ( 0.5_wp * ( u(k,j,i) + &
	1051	u(k,j,i+1) ) &
	1052	)**2 + &
	1053	( 0.5_wp * ( v(k,j,i) + &
	1054	v(k,j+1,i) ) &
	1055	)**2 + &
	1056	( 0.5_wp * ( w(k-1,j,i) + &
	1057	w(k,j,i) ) &
	1058	)**2 &
	1059	) * &
	1060	( s(k,j,i) - lsc )
	1061	ENDDO
	1062	ENDDO
	1063	ENDDO
[1484]	1064
	1065
[1960]	1066
[138]	1067	CASE DEFAULT
	1068
[257]	1069	WRITE( message_string, * ) 'wrong component: ', component
[1826]	1070	CALL message( 'pcm_tendency', 'PA0279', 1, 2, 0, 6, 0 )
[138]	1071
	1072	END SELECT
	1073
[1826]	1074	END SUBROUTINE pcm_tendency
[138]	1075
	1076
	1077	!------------------------------------------------------------------------------!
[1484]	1078	! Description:
	1079	! ------------
[1682]	1080	!> Calculation of the tendency terms, accounting for the effect of the plant
	1081	!> canopy on momentum and scalar quantities.
	1082	!>
	1083	!> The canopy is located where the leaf area density lad_s(k,j,i) > 0.0
[1826]	1084	!> (defined on scalar grid), as initialized in subroutine pcm_init.
[1682]	1085	!> The lad on the w-grid is vertically interpolated from the surrounding
	1086	!> lad_s. The upper boundary of the canopy is defined on the w-grid at
	1087	!> k = pch_index. Here, the lad is zero.
	1088	!>
	1089	!> The canopy drag must be limited (previously accounted for by calculation of
	1090	!> a limiting canopy timestep for the determination of the maximum LES timestep
	1091	!> in subroutine timestep), since it is physically impossible that the canopy
	1092	!> drag alone can locally change the sign of a velocity component. This
	1093	!> limitation is realized by calculating preliminary tendencies and velocities.
	1094	!> It is subsequently checked if the preliminary new velocity has a different
	1095	!> sign than the current velocity. If so, the tendency is limited in a way that
	1096	!> the velocity can at maximum be reduced to zero by the canopy drag.
	1097	!>
	1098	!>
	1099	!> Call for grid point i,j
[138]	1100	!------------------------------------------------------------------------------!
[1826]	1101	SUBROUTINE pcm_tendency_ij( i, j, component )
[138]	1102
	1103
[1320]	1104	USE control_parameters, &
[1484]	1105	ONLY: dt_3d, message_string
[1320]	1106
	1107	USE kinds
	1108
[138]	1109	IMPLICIT NONE
	1110
[1682]	1111	INTEGER(iwp) :: component !< prognostic variable (u,v,w,pt,q,e)
	1112	INTEGER(iwp) :: i !< running index
	1113	INTEGER(iwp) :: j !< running index
	1114	INTEGER(iwp) :: k !< running index
[1721]	1115	INTEGER(iwp) :: kk !< running index for flat lad arrays
[138]	1116
[1682]	1117	REAL(wp) :: ddt_3d !< inverse of the LES timestep (dt_3d)
	1118	REAL(wp) :: lad_local !< local lad value
	1119	REAL(wp) :: pre_tend !< preliminary tendency
	1120	REAL(wp) :: pre_u !< preliminary u-value
	1121	REAL(wp) :: pre_v !< preliminary v-value
	1122	REAL(wp) :: pre_w !< preliminary w-value
[1484]	1123
	1124
	1125	ddt_3d = 1.0_wp / dt_3d
	1126
[138]	1127	!
[1484]	1128	!-- Compute drag for the three velocity components and the SGS-TKE
[142]	1129	SELECT CASE ( component )
[138]	1130
	1131	!
[142]	1132	!-- u-component
[1484]	1133	CASE ( 1 )
[1721]	1134	DO k = nzb_u_inner(j,i)+1, nzb_u_inner(j,i)+pch_index
[138]	1135
[1721]	1136	kk = k - nzb_u_inner(j,i) !- lad arrays are defined flat
[138]	1137	!
[1484]	1138	!-- In order to create sharp boundaries of the plant canopy,
	1139	!-- the lad on the u-grid at index (k,j,i) is equal to lad_s(k,j,i),
	1140	!-- rather than being interpolated from the surrounding lad_s,
	1141	!-- because this would yield smaller lad at the canopy boundaries
	1142	!-- than inside of the canopy.
	1143	!-- For the same reason, the lad at the rightmost(i+1)canopy
	1144	!-- boundary on the u-grid equals lad_s(k,j,i).
[1721]	1145	lad_local = lad_s(kk,j,i)
	1146	IF ( lad_local == 0.0_wp .AND. lad_s(kk,j,i-1) > 0.0_wp ) THEN
	1147	lad_local = lad_s(kk,j,i-1)
[1484]	1148	ENDIF
	1149
	1150	pre_tend = 0.0_wp
	1151	pre_u = 0.0_wp
	1152	!
	1153	!-- Calculate preliminary value (pre_tend) of the tendency
	1154	pre_tend = - cdc * &
	1155	lad_local * &
	1156	SQRT( u(k,j,i)**2 + &
	1157	( 0.25_wp * ( v(k,j,i-1) + &
	1158	v(k,j,i) + &
	1159	v(k,j+1,i) + &
	1160	v(k,j+1,i-1) ) &
	1161	)**2 + &
	1162	( 0.25_wp * ( w(k-1,j,i-1) + &
	1163	w(k-1,j,i) + &
	1164	w(k,j,i-1) + &
	1165	w(k,j,i) ) &
	1166	)**2 &
	1167	) * &
	1168	u(k,j,i)
	1169
	1170	!
	1171	!-- Calculate preliminary new velocity, based on pre_tend
	1172	pre_u = u(k,j,i) + dt_3d * pre_tend
	1173	!
	1174	!-- Compare sign of old velocity and new preliminary velocity,
	1175	!-- and in case the signs are different, limit the tendency
	1176	IF ( SIGN(pre_u,u(k,j,i)) /= pre_u ) THEN
	1177	pre_tend = - u(k,j,i) * ddt_3d
	1178	ELSE
	1179	pre_tend = pre_tend
	1180	ENDIF
	1181	!
	1182	!-- Calculate final tendency
	1183	tend(k,j,i) = tend(k,j,i) + pre_tend
	1184	ENDDO
	1185
	1186
	1187	!
[142]	1188	!-- v-component
[1484]	1189	CASE ( 2 )
[1721]	1190	DO k = nzb_v_inner(j,i)+1, nzb_v_inner(j,i)+pch_index
[138]	1191
[1721]	1192	kk = k - nzb_v_inner(j,i) !- lad arrays are defined flat
[138]	1193	!
[1484]	1194	!-- In order to create sharp boundaries of the plant canopy,
	1195	!-- the lad on the v-grid at index (k,j,i) is equal to lad_s(k,j,i),
	1196	!-- rather than being interpolated from the surrounding lad_s,
	1197	!-- because this would yield smaller lad at the canopy boundaries
	1198	!-- than inside of the canopy.
	1199	!-- For the same reason, the lad at the northmost(j+1)canopy
	1200	!-- boundary on the v-grid equals lad_s(k,j,i).
[1721]	1201	lad_local = lad_s(kk,j,i)
	1202	IF ( lad_local == 0.0_wp .AND. lad_s(kk,j-1,i) > 0.0_wp ) THEN
	1203	lad_local = lad_s(kk,j-1,i)
[1484]	1204	ENDIF
	1205
	1206	pre_tend = 0.0_wp
	1207	pre_v = 0.0_wp
	1208	!
	1209	!-- Calculate preliminary value (pre_tend) of the tendency
	1210	pre_tend = - cdc * &
	1211	lad_local * &
	1212	SQRT( ( 0.25_wp * ( u(k,j-1,i) + &
	1213	u(k,j-1,i+1) + &
	1214	u(k,j,i) + &
	1215	u(k,j,i+1) ) &
	1216	)**2 + &
	1217	v(k,j,i)**2 + &
	1218	( 0.25_wp * ( w(k-1,j-1,i) + &
	1219	w(k-1,j,i) + &
	1220	w(k,j-1,i) + &
	1221	w(k,j,i) ) &
	1222	)**2 &
	1223	) * &
	1224	v(k,j,i)
	1225
	1226	!
	1227	!-- Calculate preliminary new velocity, based on pre_tend
	1228	pre_v = v(k,j,i) + dt_3d * pre_tend
	1229	!
	1230	!-- Compare sign of old velocity and new preliminary velocity,
	1231	!-- and in case the signs are different, limit the tendency
	1232	IF ( SIGN(pre_v,v(k,j,i)) /= pre_v ) THEN
	1233	pre_tend = - v(k,j,i) * ddt_3d
	1234	ELSE
	1235	pre_tend = pre_tend
	1236	ENDIF
	1237	!
	1238	!-- Calculate final tendency
	1239	tend(k,j,i) = tend(k,j,i) + pre_tend
	1240	ENDDO
	1241
	1242
	1243	!
[142]	1244	!-- w-component
[1484]	1245	CASE ( 3 )
[1721]	1246	DO k = nzb_w_inner(j,i)+1, nzb_w_inner(j,i)+pch_index-1
[138]	1247
[1721]	1248	kk = k - nzb_w_inner(j,i) !- lad arrays are defined flat
	1249
[1484]	1250	pre_tend = 0.0_wp
	1251	pre_w = 0.0_wp
[138]	1252	!
[1484]	1253	!-- Calculate preliminary value (pre_tend) of the tendency
	1254	pre_tend = - cdc * &
	1255	(0.5_wp * &
[1721]	1256	( lad_s(kk+1,j,i) + lad_s(kk,j,i) )) * &
[1484]	1257	SQRT( ( 0.25_wp * ( u(k,j,i) + &
	1258	u(k,j,i+1) + &
	1259	u(k+1,j,i) + &
	1260	u(k+1,j,i+1) ) &
	1261	)**2 + &
	1262	( 0.25_wp * ( v(k,j,i) + &
	1263	v(k,j+1,i) + &
	1264	v(k+1,j,i) + &
	1265	v(k+1,j+1,i) ) &
	1266	)**2 + &
	1267	w(k,j,i)**2 &
	1268	) * &
	1269	w(k,j,i)
	1270	!
	1271	!-- Calculate preliminary new velocity, based on pre_tend
	1272	pre_w = w(k,j,i) + dt_3d * pre_tend
	1273	!
	1274	!-- Compare sign of old velocity and new preliminary velocity,
	1275	!-- and in case the signs are different, limit the tendency
	1276	IF ( SIGN(pre_w,w(k,j,i)) /= pre_w ) THEN
	1277	pre_tend = - w(k,j,i) * ddt_3d
	1278	ELSE
	1279	pre_tend = pre_tend
	1280	ENDIF
	1281	!
	1282	!-- Calculate final tendency
	1283	tend(k,j,i) = tend(k,j,i) + pre_tend
	1284	ENDDO
	1285
	1286	!
[153]	1287	!-- potential temperature
	1288	CASE ( 4 )
[1721]	1289	DO k = nzb_s_inner(j,i)+1, nzb_s_inner(j,i)+pch_index
	1290	kk = k - nzb_s_inner(j,i) !- lad arrays are defined flat
[2011]	1291	tend(k,j,i) = tend(k,j,i) + pc_heating_rate(kk,j,i)
[153]	1292	ENDDO
	1293
	1294
	1295	!
[1960]	1296	!-- humidity
[153]	1297	CASE ( 5 )
[1721]	1298	DO k = nzb_s_inner(j,i)+1, nzb_s_inner(j,i)+pch_index
	1299	kk = k - nzb_s_inner(j,i) !- lad arrays are defined flat
[1484]	1300	tend(k,j,i) = tend(k,j,i) - &
	1301	lsec * &
[1721]	1302	lad_s(kk,j,i) * &
[1484]	1303	SQRT( ( 0.5_wp * ( u(k,j,i) + &
	1304	u(k,j,i+1) ) &
	1305	)**2 + &
	1306	( 0.5_wp * ( v(k,j,i) + &
	1307	v(k,j+1,i) ) &
	1308	)**2 + &
	1309	( 0.5_wp * ( w(k-1,j,i) + &
	1310	w(k,j,i) ) &
	1311	)**2 &
	1312	) * &
	1313	( q(k,j,i) - lsc )
[153]	1314	ENDDO
	1315
	1316	!
[142]	1317	!-- sgs-tke
[1484]	1318	CASE ( 6 )
[1721]	1319	DO k = nzb_s_inner(j,i)+1, nzb_s_inner(j,i)+pch_index
	1320	kk = k - nzb_s_inner(j,i) !- lad arrays are defined flat
[1484]	1321	tend(k,j,i) = tend(k,j,i) - &
	1322	2.0_wp * cdc * &
[1721]	1323	lad_s(kk,j,i) * &
[1484]	1324	SQRT( ( 0.5_wp * ( u(k,j,i) + &
	1325	u(k,j,i+1) ) &
	1326	)**2 + &
	1327	( 0.5_wp * ( v(k,j,i) + &
	1328	v(k,j+1,i) ) &
	1329	)**2 + &
	1330	( 0.5_wp * ( w(k,j,i) + &
	1331	w(k+1,j,i) ) &
	1332	)**2 &
	1333	) * &
	1334	e(k,j,i)
	1335	ENDDO
[1960]	1336
	1337	!
	1338	!-- scalar concentration
	1339	CASE ( 7 )
	1340	DO k = nzb_s_inner(j,i)+1, nzb_s_inner(j,i)+pch_index
	1341	kk = k - nzb_s_inner(j,i) !- lad arrays are defined flat
	1342	tend(k,j,i) = tend(k,j,i) - &
	1343	lsec * &
	1344	lad_s(kk,j,i) * &
	1345	SQRT( ( 0.5_wp * ( u(k,j,i) + &
	1346	u(k,j,i+1) ) &
	1347	)**2 + &
	1348	( 0.5_wp * ( v(k,j,i) + &
	1349	v(k,j+1,i) ) &
	1350	)**2 + &
	1351	( 0.5_wp * ( w(k-1,j,i) + &
	1352	w(k,j,i) ) &
	1353	)**2 &
	1354	) * &
	1355	( s(k,j,i) - lsc )
	1356	ENDDO
[138]	1357
[142]	1358	CASE DEFAULT
[138]	1359
[257]	1360	WRITE( message_string, * ) 'wrong component: ', component
[1826]	1361	CALL message( 'pcm_tendency', 'PA0279', 1, 2, 0, 6, 0 )
[138]	1362
[142]	1363	END SELECT
[138]	1364
[1826]	1365	END SUBROUTINE pcm_tendency_ij
[138]	1366
[2007]	1367
	1368
[138]	1369	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 |