Home

Context Navigation

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

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