Home

Context Navigation

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

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