Home

Context Navigation

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

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