Home

Context Navigation

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

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