Home

Context Navigation

source: palm/trunk/SOURCE/plant_canopy_model.f90 @ 1625

Last change on this file since 1625 was 1485, checked in by kanani, 10 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 40.7 KB

Rev	Line
[138]	1	MODULE plant_canopy_model_mod
	2
[1036]	3	!--------------------------------------------------------------------------------!
	4	! This file is part of PALM.
	5	!
	6	! PALM is free software: you can redistribute it and/or modify it under the terms
	7	! of the GNU General Public License as published by the Free Software Foundation,
	8	! either version 3 of the License, or (at your option) any later version.
	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	!
[1310]	17	! Copyright 1997-2014 Leibniz Universitaet Hannover
[1036]	18	!--------------------------------------------------------------------------------!
	19	!
[257]	20	! Current revisions:
[138]	21	! -----------------
[1485]	22	!
	23	!
	24	! Former revisions:
	25	! -----------------
	26	! $Id: plant_canopy_model.f90 1485 2014-10-21 11:09:54Z heinze $
	27	!
	28	! 1484 2014-10-21 10:53:05Z kanani
[1484]	29	! Changes due to new module structure of the plant canopy model:
	30	! module plant_canopy_model_mod now contains a subroutine for the
	31	! initialization of the canopy model (init_plant_canopy),
	32	! limitation of the canopy drag (previously accounted for by calculation of
	33	! a limiting canopy timestep for the determination of the maximum LES timestep
	34	! in subroutine timestep) is now realized by the calculation of pre-tendencies
	35	! and preliminary velocities in subroutine plant_canopy_model,
	36	! some redundant MPI communication removed in subroutine init_plant_canopy
	37	! (was previously in init_3d_model),
	38	! unnecessary 3d-arrays lad_u, lad_v, lad_w removed - lad information on the
	39	! respective grid is now provided only by lad_s (e.g. in the calculation of
	40	! the tendency terms or of cum_lai_hf),
	41	! drag_coefficient, lai, leaf_surface_concentration,
	42	! scalar_exchange_coefficient, sec and sls renamed to canopy_drag_coeff,
	43	! cum_lai_hf, leaf_surface_conc, leaf_scalar_exch_coeff, lsec and lsc,
	44	! respectively,
	45	! unnecessary 3d-arrays cdc, lsc and lsec now defined as single-value constants,
	46	! USE-statements and ONLY-lists modified accordingly
[1341]	47	!
	48	! 1340 2014-03-25 19:45:13Z kanani
	49	! REAL constants defined as wp-kind
	50	!
[1321]	51	! 1320 2014-03-20 08:40:49Z raasch
[1320]	52	! ONLY-attribute added to USE-statements,
	53	! kind-parameters added to all INTEGER and REAL declaration statements,
	54	! kinds are defined in new module kinds,
	55	! old module precision_kind is removed,
	56	! revision history before 2012 removed,
	57	! comment fields (!:) to be used for variable explanations added to
	58	! all variable declaration statements
[153]	59	!
[1037]	60	! 1036 2012-10-22 13:43:42Z raasch
	61	! code put under GPL (PALM 3.9)
	62	!
[139]	63	! 138 2007-11-28 10:03:58Z letzel
	64	! Initial revision
	65	!
[138]	66	! Description:
	67	! ------------
[1484]	68	! 1) Initialization of the canopy model, e.g. construction of leaf area density
	69	! profile (subroutine init_plant_canopy).
	70	! 2) Calculation of sinks and sources of momentum, heat and scalar concentration
	71	! due to canopy elements (subroutine plant_canopy_model).
[138]	72	!------------------------------------------------------------------------------!
[1484]	73	USE arrays_3d, &
	74	ONLY: dzu, dzw, e, q, shf, tend, u, v, w, zu, zw
[138]	75
[1484]	76	USE indices, &
	77	ONLY: nbgp, nxl, nxlg, nxlu, nxr, nxrg, nyn, nyng, nys, nysg, nysv, &
	78	nz, nzb, nzb_s_inner, nzb_u_inner, nzb_v_inner, nzb_w_inner, nzt
	79
	80	USE kinds
	81
	82
	83	IMPLICIT NONE
	84
	85
	86	CHARACTER (LEN=20) :: canopy_mode = 'block' !: canopy coverage
	87
	88	INTEGER(iwp) :: pch_index = 0 !: plant canopy height/top index
	89	INTEGER(iwp) :: &
	90	lad_vertical_gradient_level_ind(10) = -9999 !: lad-profile levels (index)
	91
	92	LOGICAL :: calc_beta_lad_profile = .FALSE. !: switch for calc. of lad from beta func.
	93	LOGICAL :: plant_canopy = .FALSE. !: switch for use of canopy model
	94
	95	REAL(wp) :: alpha_lad = 9999999.9_wp !: coefficient for lad calculation
	96	REAL(wp) :: beta_lad = 9999999.9_wp !: coefficient for lad calculation
	97	REAL(wp) :: canopy_drag_coeff = 0.0_wp !: canopy drag coefficient (parameter)
	98	REAL(wp) :: cdc = 0.0_wp !: canopy drag coeff. (abbreviation used in equations)
	99	REAL(wp) :: cthf = 0.0_wp !: canopy top heat flux
	100	REAL(wp) :: dt_plant_canopy = 0.0_wp !: timestep account. for canopy drag
	101	REAL(wp) :: lad_surface = 0.0_wp !: lad surface value
	102	REAL(wp) :: lai_beta = 0.0_wp !: leaf area index (lai) for lad calc.
	103	REAL(wp) :: &
	104	leaf_scalar_exch_coeff = 0.0_wp !: canopy scalar exchange coeff.
	105	REAL(wp) :: &
	106	leaf_surface_conc = 0.0_wp !: leaf surface concentration
	107	REAL(wp) :: lsec = 0.0_wp !: leaf scalar exchange coeff.
	108	REAL(wp) :: lsc = 0.0_wp !: leaf surface concentration
	109
	110	REAL(wp) :: &
	111	lad_vertical_gradient(10) = 0.0_wp !: lad gradient
	112	REAL(wp) :: &
	113	lad_vertical_gradient_level(10) = -9999999.9_wp !: lad-prof. levels (in m)
	114
	115	REAL(wp), DIMENSION(:), ALLOCATABLE :: lad !: leaf area density
	116	REAL(wp), DIMENSION(:), ALLOCATABLE :: pre_lad !: preliminary lad
	117
	118	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
	119	canopy_heat_flux !: canopy heat flux
	120	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: cum_lai_hf !: cumulative lai for heatflux calc.
	121	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: lad_s !: lad on scalar-grid
	122
	123
	124	SAVE
	125
	126
[138]	127	PRIVATE
[1484]	128	PUBLIC alpha_lad, beta_lad, calc_beta_lad_profile, canopy_drag_coeff, &
	129	canopy_mode, cdc, cthf, dt_plant_canopy, init_plant_canopy, lad, &
	130	lad_s, lad_surface, lad_vertical_gradient, &
	131	lad_vertical_gradient_level, lad_vertical_gradient_level_ind, &
	132	lai_beta, leaf_scalar_exch_coeff, leaf_surface_conc, lsc, lsec, &
	133	pch_index, plant_canopy, plant_canopy_model
[138]	134
[1484]	135
	136	INTERFACE init_plant_canopy
	137	MODULE PROCEDURE init_plant_canopy
	138	END INTERFACE init_plant_canopy
	139
[138]	140	INTERFACE plant_canopy_model
	141	MODULE PROCEDURE plant_canopy_model
	142	MODULE PROCEDURE plant_canopy_model_ij
	143	END INTERFACE plant_canopy_model
	144
[1484]	145
	146
	147
[138]	148	CONTAINS
	149
	150
	151	!------------------------------------------------------------------------------!
[1484]	152	! Description:
	153	! ------------
	154	!-- Initialization of the plant canopy model
[138]	155	!------------------------------------------------------------------------------!
[1484]	156	SUBROUTINE init_plant_canopy
	157
	158
	159	USE control_parameters, &
	160	ONLY: dz, ocean, passive_scalar
	161
	162
	163	IMPLICIT NONE
	164
	165	INTEGER(iwp) :: i !: running index
	166	INTEGER(iwp) :: j !: running index
	167	INTEGER(iwp) :: k !: running index
	168
	169	REAL(wp) :: int_bpdf !: vertical integral for lad-profile construction
	170	REAL(wp) :: dzh !: vertical grid spacing in units of canopy height
	171	REAL(wp) :: gradient !: gradient for lad-profile construction
	172	REAL(wp) :: canopy_height !: canopy height for lad-profile construction
	173
	174	!
	175	!-- Allocate one-dimensional arrays for the computation of the
	176	!-- leaf area density (lad) profile
	177	ALLOCATE( lad(0:nz+1), pre_lad(0:nz+1) )
	178	lad = 0.0_wp
	179	pre_lad = 0.0_wp
	180
	181	!
	182	!-- Compute the profile of leaf area density used in the plant
	183	!-- canopy model. The profile can either be constructed from
	184	!-- prescribed vertical gradients of the leaf area density or by
	185	!-- using a beta probability density function (see e.g. Markkanen et al.,
	186	!-- 2003: Boundary-Layer Meteorology, 106, 437-459)
	187	IF ( .NOT. calc_beta_lad_profile ) THEN
	188
	189	!
	190	!-- Use vertical gradients for lad-profile construction
	191	i = 1
	192	gradient = 0.0_wp
	193
	194	IF ( .NOT. ocean ) THEN
	195
	196	lad(0) = lad_surface
	197	lad_vertical_gradient_level_ind(1) = 0
	198
	199	DO k = 1, pch_index
	200	IF ( i < 11 ) THEN
	201	IF ( lad_vertical_gradient_level(i) < zu(k) .AND. &
	202	lad_vertical_gradient_level(i) >= 0.0_wp ) THEN
	203	gradient = lad_vertical_gradient(i)
	204	lad_vertical_gradient_level_ind(i) = k - 1
	205	i = i + 1
	206	ENDIF
	207	ENDIF
	208	IF ( gradient /= 0.0_wp ) THEN
	209	IF ( k /= 1 ) THEN
	210	lad(k) = lad(k-1) + dzu(k) * gradient
	211	ELSE
	212	lad(k) = lad_surface + dzu(k) * gradient
	213	ENDIF
	214	ELSE
	215	lad(k) = lad(k-1)
	216	ENDIF
	217	ENDDO
	218
	219	ENDIF
	220
	221	!
	222	!-- In case of no given leaf area density gradients, choose a vanishing
	223	!-- gradient. This information is used for the HEADER and the RUN_CONTROL
	224	!-- file.
	225	IF ( lad_vertical_gradient_level(1) == -9999999.9_wp ) THEN
	226	lad_vertical_gradient_level(1) = 0.0_wp
	227	ENDIF
	228
	229	ELSE
	230
	231	!
	232	!-- Use beta function for lad-profile construction
	233	int_bpdf = 0.0_wp
	234	canopy_height = pch_index * dz
	235
	236	DO k = nzb, pch_index
	237	int_bpdf = int_bpdf + &
	238	( ( ( zw(k) / canopy_height )*( alpha_lad-1.0_wp ) ) &
	239	( ( 1.0_wp - ( zw(k) / canopy_height ) )*( beta_lad-1.0_wp ) ) &
	240	( ( zw(k+1)-zw(k) ) / canopy_height ) )
	241	ENDDO
	242
	243	!
	244	!-- Preliminary lad profile (defined on w-grid)
	245	DO k = nzb, pch_index
	246	pre_lad(k) = lai_beta * &
	247	( ( ( zw(k) / canopy_height )*( alpha_lad-1.0_wp ) ) &
	248	( ( 1.0_wp - ( zw(k) / canopy_height ) )**( beta_lad-1.0_wp ) ) / &
	249	int_bpdf &
	250	) / canopy_height
	251	ENDDO
	252
	253	!
	254	!-- Final lad profile (defined on scalar-grid level, since most prognostic
	255	!-- quantities are defined there, hence, less interpolation is required
	256	!-- when calculating the canopy tendencies)
	257	lad(0) = pre_lad(0)
	258	DO k = nzb+1, pch_index
	259	lad(k) = 0.5 * ( pre_lad(k-1) + pre_lad(k) )
	260	ENDDO
	261
	262	ENDIF
	263
	264	!
	265	!-- Allocate 3D-array for the leaf area density (lad_s). In case of a
	266	!-- prescribed canopy-top heat flux (cthf), allocate 3D-arrays for
	267	!-- the cumulative leaf area index (cum_lai_hf) and the canopy heat flux.
	268	ALLOCATE( lad_s(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
	269
	270	IF ( cthf /= 0.0_wp ) THEN
	271	ALLOCATE( cum_lai_hf(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
	272	canopy_heat_flux(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
	273	ENDIF
	274
	275	!
	276	!-- Initialize canopy parameters cdc (canopy drag coefficient),
	277	!-- lsec (leaf scalar exchange coefficient), lsc (leaf surface concentration)
	278	!-- with the prescribed values
	279	cdc = canopy_drag_coeff
	280	lsec = leaf_scalar_exch_coeff
	281	lsc = leaf_surface_conc
	282
	283	!
	284	!-- Initialization of the canopy coverage in the model domain:
	285	!-- Setting the parameter canopy_mode = 'block' initializes a canopy, which
	286	!-- fully covers the domain surface
	287	SELECT CASE ( TRIM( canopy_mode ) )
	288
	289	CASE( 'block' )
	290
	291	DO i = nxlg, nxrg
	292	DO j = nysg, nyng
	293	lad_s(:,j,i) = lad(:)
	294	ENDDO
	295	ENDDO
	296
	297	CASE DEFAULT
	298
	299	!
	300	!-- The DEFAULT case is reached either if the parameter
	301	!-- canopy mode contains a wrong character string or if the
	302	!-- user has coded a special case in the user interface.
	303	!-- There, the subroutine user_init_plant_canopy checks
	304	!-- which of these two conditions applies.
	305	CALL user_init_plant_canopy
	306
	307	END SELECT
	308
	309	!
	310	!-- Initialization of the canopy heat source distribution
	311	IF ( cthf /= 0.0_wp ) THEN
	312	!
	313	!-- Piecewise calculation of the leaf area index by vertical
	314	!-- integration of the leaf area density
	315	cum_lai_hf(:,:,:) = 0.0_wp
	316	DO i = nxlg, nxrg
	317	DO j = nysg, nyng
	318	DO k = pch_index-1, 0, -1
	319	IF ( k == pch_index-1 ) THEN
	320	cum_lai_hf(k,j,i) = cum_lai_hf(k+1,j,i) + &
	321	( 0.5_wp * lad_s(k+1,j,i) * &
	322	( zw(k+1) - zu(k+1) ) ) + &
	323	( 0.5_wp * ( 0.5_wp * ( lad_s(k+1,j,i) + &
	324	lad_s(k,j,i) ) + &
	325	lad_s(k+1,j,i) ) * &
	326	( zu(k+1) - zw(k) ) )
	327	ELSE
	328	cum_lai_hf(k,j,i) = cum_lai_hf(k+1,j,i) + &
	329	( 0.5_wp * ( 0.5_wp * ( lad_s(k+2,j,i) + &
	330	lad_s(k+1,j,i) ) + &
	331	lad_s(k+1,j,i) ) * &
	332	( zw(k+1) - zu(k+1) ) ) + &
	333	( 0.5_wp * ( 0.5_wp * ( lad_s(k+1,j,i) + &
	334	lad_s(k,j,i) ) + &
	335	lad_s(k+1,j,i) ) * &
	336	( zu(k+1) - zw(k) ) )
	337	ENDIF
	338	ENDDO
	339	ENDDO
	340	ENDDO
	341
	342	!
	343	!-- Calculation of the upward kinematic vertical heat flux within the
	344	!-- canopy
	345	DO i = nxlg, nxrg
	346	DO j = nysg, nyng
	347	DO k = 0, pch_index
	348	canopy_heat_flux(k,j,i) = cthf * &
	349	exp( -0.6_wp * cum_lai_hf(k,j,i) )
	350	ENDDO
	351	ENDDO
	352	ENDDO
	353
	354	!
	355	!-- The surface heat flux is set to the surface value of the calculated
	356	!-- in-canopy heat flux distribution
	357	shf(:,:) = canopy_heat_flux(0,:,:)
	358
	359	ENDIF
	360
	361
	362
	363	END SUBROUTINE init_plant_canopy
	364
	365
	366
	367	!------------------------------------------------------------------------------!
	368	! Description:
	369	! ------------
	370	!-- Calculation of the tendency terms, accounting for the effect of the plant
	371	!-- canopy on momentum and scalar quantities.
	372	!--
	373	!-- The canopy is located where the leaf area density lad_s(k,j,i) > 0.0
	374	!-- (defined on scalar grid), as initialized in subroutine init_plant_canopy.
	375	!-- The lad on the w-grid is vertically interpolated from the surrounding
	376	!-- lad_s. The upper boundary of the canopy is defined on the w-grid at
	377	!-- k = pch_index. Here, the lad is zero.
	378	!--
	379	!-- The canopy drag must be limited (previously accounted for by calculation of
	380	!-- a limiting canopy timestep for the determination of the maximum LES timestep
	381	!-- in subroutine timestep), since it is physically impossible that the canopy
	382	!-- drag alone can locally change the sign of a velocity component. This
	383	!-- limitation is realized by calculating preliminary tendencies and velocities.
	384	!-- It is subsequently checked if the preliminary new velocity has a different
	385	!-- sign than the current velocity. If so, the tendency is limited in a way that
	386	!-- the velocity can at maximum be reduced to zero by the canopy drag.
	387	!--
	388	!--
	389	!-- Call for all grid points
	390	!------------------------------------------------------------------------------!
[138]	391	SUBROUTINE plant_canopy_model( component )
	392
	393
[1320]	394	USE control_parameters, &
[1484]	395	ONLY: dt_3d, message_string
[1320]	396
	397	USE kinds
	398
[138]	399	IMPLICIT NONE
	400
[1484]	401	INTEGER(iwp) :: component !: prognostic variable (u,v,w,pt,q,e)
	402	INTEGER(iwp) :: i !: running index
	403	INTEGER(iwp) :: j !: running index
	404	INTEGER(iwp) :: k !: running index
	405
	406	REAL(wp) :: ddt_3d !: inverse of the LES timestep (dt_3d)
	407	REAL(wp) :: lad_local !: local lad value
	408	REAL(wp) :: pre_tend !: preliminary tendency
	409	REAL(wp) :: pre_u !: preliminary u-value
	410	REAL(wp) :: pre_v !: preliminary v-value
	411	REAL(wp) :: pre_w !: preliminary w-value
	412
	413
	414	ddt_3d = 1.0_wp / dt_3d
[138]	415
	416	!
[1484]	417	!-- Compute drag for the three velocity components and the SGS-TKE:
[138]	418	SELECT CASE ( component )
	419
	420	!
	421	!-- u-component
	422	CASE ( 1 )
	423	DO i = nxlu, nxr
	424	DO j = nys, nyn
	425	DO k = nzb_u_inner(j,i)+1, pch_index
[1484]	426
	427	!
	428	!-- In order to create sharp boundaries of the plant canopy,
	429	!-- the lad on the u-grid at index (k,j,i) is equal to
	430	!-- lad_s(k,j,i), rather than being interpolated from the
	431	!-- surrounding lad_s, because this would yield smaller lad
	432	!-- at the canopy boundaries than inside of the canopy.
	433	!-- For the same reason, the lad at the rightmost(i+1)canopy
	434	!-- boundary on the u-grid equals lad_s(k,j,i).
	435	lad_local = lad_s(k,j,i)
	436	IF ( lad_local == 0.0_wp .AND. &
	437	lad_s(k,j,i-1) > 0.0_wp ) THEN
	438	lad_local = lad_s(k,j,i-1)
	439	ENDIF
	440
	441	pre_tend = 0.0_wp
	442	pre_u = 0.0_wp
	443	!
	444	!-- Calculate preliminary value (pre_tend) of the tendency
	445	pre_tend = - cdc * &
	446	lad_local * &
	447	SQRT( u(k,j,i)**2 + &
	448	( 0.25_wp * ( v(k,j,i-1) + &
	449	v(k,j,i) + &
	450	v(k,j+1,i) + &
	451	v(k,j+1,i-1) ) &
	452	)**2 + &
	453	( 0.25_wp * ( w(k-1,j,i-1) + &
	454	w(k-1,j,i) + &
	455	w(k,j,i-1) + &
	456	w(k,j,i) ) &
	457	)**2 &
	458	) * &
	459	u(k,j,i)
	460
	461	!
	462	!-- Calculate preliminary new velocity, based on pre_tend
	463	pre_u = u(k,j,i) + dt_3d * pre_tend
	464	!
	465	!-- Compare sign of old velocity and new preliminary velocity,
	466	!-- and in case the signs are different, limit the tendency
	467	IF ( SIGN(pre_u,u(k,j,i)) /= pre_u ) THEN
	468	pre_tend = - u(k,j,i) * ddt_3d
	469	ELSE
	470	pre_tend = pre_tend
	471	ENDIF
	472	!
	473	!-- Calculate final tendency
	474	tend(k,j,i) = tend(k,j,i) + pre_tend
	475
[138]	476	ENDDO
	477	ENDDO
	478	ENDDO
	479
	480	!
	481	!-- v-component
	482	CASE ( 2 )
	483	DO i = nxl, nxr
	484	DO j = nysv, nyn
	485	DO k = nzb_v_inner(j,i)+1, pch_index
[1484]	486
	487	!
	488	!-- In order to create sharp boundaries of the plant canopy,
	489	!-- the lad on the v-grid at index (k,j,i) is equal to
	490	!-- lad_s(k,j,i), rather than being interpolated from the
	491	!-- surrounding lad_s, because this would yield smaller lad
	492	!-- at the canopy boundaries than inside of the canopy.
	493	!-- For the same reason, the lad at the northmost(j+1) canopy
	494	!-- boundary on the v-grid equals lad_s(k,j,i).
	495	lad_local = lad_s(k,j,i)
	496	IF ( lad_local == 0.0_wp .AND. &
	497	lad_s(k,j-1,i) > 0.0_wp ) THEN
	498	lad_local = lad_s(k,j-1,i)
	499	ENDIF
	500
	501	pre_tend = 0.0_wp
	502	pre_v = 0.0_wp
	503	!
	504	!-- Calculate preliminary value (pre_tend) of the tendency
	505	pre_tend = - cdc * &
	506	lad_local * &
	507	SQRT( ( 0.25_wp * ( u(k,j-1,i) + &
	508	u(k,j-1,i+1) + &
	509	u(k,j,i) + &
	510	u(k,j,i+1) ) &
	511	)**2 + &
	512	v(k,j,i)**2 + &
	513	( 0.25_wp * ( w(k-1,j-1,i) + &
	514	w(k-1,j,i) + &
	515	w(k,j-1,i) + &
	516	w(k,j,i) ) &
	517	)**2 &
	518	) * &
	519	v(k,j,i)
	520
	521	!
	522	!-- Calculate preliminary new velocity, based on pre_tend
	523	pre_v = v(k,j,i) + dt_3d * pre_tend
	524	!
	525	!-- Compare sign of old velocity and new preliminary velocity,
	526	!-- and in case the signs are different, limit the tendency
	527	IF ( SIGN(pre_v,v(k,j,i)) /= pre_v ) THEN
	528	pre_tend = - v(k,j,i) * ddt_3d
	529	ELSE
	530	pre_tend = pre_tend
	531	ENDIF
	532	!
	533	!-- Calculate final tendency
	534	tend(k,j,i) = tend(k,j,i) + pre_tend
	535
[138]	536	ENDDO
	537	ENDDO
	538	ENDDO
	539
	540	!
	541	!-- w-component
	542	CASE ( 3 )
	543	DO i = nxl, nxr
	544	DO j = nys, nyn
[1484]	545	DO k = nzb_w_inner(j,i)+1, pch_index-1
	546
	547	pre_tend = 0.0_wp
	548	pre_w = 0.0_wp
	549	!
	550	!-- Calculate preliminary value (pre_tend) of the tendency
	551	pre_tend = - cdc * &
	552	(0.5_wp * &
	553	( lad_s(k+1,j,i) + lad_s(k,j,i) )) * &
	554	SQRT( ( 0.25_wp * ( u(k,j,i) + &
	555	u(k,j,i+1) + &
	556	u(k+1,j,i) + &
	557	u(k+1,j,i+1) ) &
	558	)**2 + &
	559	( 0.25_wp * ( v(k,j,i) + &
	560	v(k,j+1,i) + &
	561	v(k+1,j,i) + &
	562	v(k+1,j+1,i) ) &
	563	)**2 + &
	564	w(k,j,i)**2 &
	565	) * &
	566	w(k,j,i)
	567	!
	568	!-- Calculate preliminary new velocity, based on pre_tend
	569	pre_w = w(k,j,i) + dt_3d * pre_tend
	570	!
	571	!-- Compare sign of old velocity and new preliminary velocity,
	572	!-- and in case the signs are different, limit the tendency
	573	IF ( SIGN(pre_w,w(k,j,i)) /= pre_w ) THEN
	574	pre_tend = - w(k,j,i) * ddt_3d
	575	ELSE
	576	pre_tend = pre_tend
	577	ENDIF
	578	!
	579	!-- Calculate final tendency
	580	tend(k,j,i) = tend(k,j,i) + pre_tend
	581
[138]	582	ENDDO
	583	ENDDO
	584	ENDDO
	585
	586	!
[153]	587	!-- potential temperature
[138]	588	CASE ( 4 )
	589	DO i = nxl, nxr
	590	DO j = nys, nyn
	591	DO k = nzb_s_inner(j,i)+1, pch_index
[1484]	592	tend(k,j,i) = tend(k,j,i) + &
	593	( canopy_heat_flux(k,j,i) - &
	594	canopy_heat_flux(k-1,j,i) ) / dzw(k)
[153]	595	ENDDO
	596	ENDDO
	597	ENDDO
	598
	599	!
	600	!-- scalar concentration
	601	CASE ( 5 )
	602	DO i = nxl, nxr
	603	DO j = nys, nyn
	604	DO k = nzb_s_inner(j,i)+1, pch_index
[1484]	605	tend(k,j,i) = tend(k,j,i) - &
	606	lsec * &
	607	lad_s(k,j,i) * &
	608	SQRT( ( 0.5_wp * ( u(k,j,i) + &
	609	u(k,j,i+1) ) &
	610	)**2 + &
	611	( 0.5_wp * ( v(k,j,i) + &
	612	v(k,j+1,i) ) &
	613	)**2 + &
	614	( 0.5_wp * ( w(k-1,j,i) + &
	615	w(k,j,i) ) &
	616	)**2 &
	617	) * &
	618	( q(k,j,i) - lsc )
[153]	619	ENDDO
	620	ENDDO
	621	ENDDO
	622
	623	!
	624	!-- sgs-tke
	625	CASE ( 6 )
	626	DO i = nxl, nxr
	627	DO j = nys, nyn
	628	DO k = nzb_s_inner(j,i)+1, pch_index
[1484]	629	tend(k,j,i) = tend(k,j,i) - &
	630	2.0_wp * cdc * &
	631	lad_s(k,j,i) * &
	632	SQRT( ( 0.5_wp * ( u(k,j,i) + &
	633	u(k,j,i+1) ) &
	634	)**2 + &
	635	( 0.5_wp * ( v(k,j,i) + &
	636	v(k,j+1,i) ) &
	637	)**2 + &
	638	( 0.5_wp * ( w(k,j,i) + &
	639	w(k+1,j,i) ) &
	640	)**2 &
	641	) * &
	642	e(k,j,i)
[138]	643	ENDDO
	644	ENDDO
	645	ENDDO
[1484]	646
	647
[138]	648	CASE DEFAULT
	649
[257]	650	WRITE( message_string, * ) 'wrong component: ', component
	651	CALL message( 'plant_canopy_model', 'PA0279', 1, 2, 0, 6, 0 )
[138]	652
	653	END SELECT
	654
	655	END SUBROUTINE plant_canopy_model
	656
	657
	658	!------------------------------------------------------------------------------!
[1484]	659	! Description:
	660	! ------------
	661	!-- Calculation of the tendency terms, accounting for the effect of the plant
	662	!-- canopy on momentum and scalar quantities.
	663	!--
	664	!-- The canopy is located where the leaf area density lad_s(k,j,i) > 0.0
	665	!-- (defined on scalar grid), as initialized in subroutine init_plant_canopy.
	666	!-- The lad on the w-grid is vertically interpolated from the surrounding
	667	!-- lad_s. The upper boundary of the canopy is defined on the w-grid at
	668	!-- k = pch_index. Here, the lad is zero.
	669	!--
	670	!-- The canopy drag must be limited (previously accounted for by calculation of
	671	!-- a limiting canopy timestep for the determination of the maximum LES timestep
	672	!-- in subroutine timestep), since it is physically impossible that the canopy
	673	!-- drag alone can locally change the sign of a velocity component. This
	674	!-- limitation is realized by calculating preliminary tendencies and velocities.
	675	!-- It is subsequently checked if the preliminary new velocity has a different
	676	!-- sign than the current velocity. If so, the tendency is limited in a way that
	677	!-- the velocity can at maximum be reduced to zero by the canopy drag.
	678	!--
	679	!--
	680	!-- Call for grid point i,j
[138]	681	!------------------------------------------------------------------------------!
	682	SUBROUTINE plant_canopy_model_ij( i, j, component )
	683
	684
[1320]	685	USE control_parameters, &
[1484]	686	ONLY: dt_3d, message_string
[1320]	687
	688	USE kinds
	689
[138]	690	IMPLICIT NONE
	691
[1484]	692	INTEGER(iwp) :: component !: prognostic variable (u,v,w,pt,q,e)
	693	INTEGER(iwp) :: i !: running index
	694	INTEGER(iwp) :: j !: running index
	695	INTEGER(iwp) :: k !: running index
[138]	696
[1484]	697	REAL(wp) :: ddt_3d !: inverse of the LES timestep (dt_3d)
	698	REAL(wp) :: lad_local !: local lad value
	699	REAL(wp) :: pre_tend !: preliminary tendency
	700	REAL(wp) :: pre_u !: preliminary u-value
	701	REAL(wp) :: pre_v !: preliminary v-value
	702	REAL(wp) :: pre_w !: preliminary w-value
	703
	704
	705	ddt_3d = 1.0_wp / dt_3d
	706
[138]	707	!
[1484]	708	!-- Compute drag for the three velocity components and the SGS-TKE
[142]	709	SELECT CASE ( component )
[138]	710
	711	!
[142]	712	!-- u-component
[1484]	713	CASE ( 1 )
	714	DO k = nzb_u_inner(j,i)+1, pch_index
[138]	715
	716	!
[1484]	717	!-- In order to create sharp boundaries of the plant canopy,
	718	!-- the lad on the u-grid at index (k,j,i) is equal to lad_s(k,j,i),
	719	!-- rather than being interpolated from the surrounding lad_s,
	720	!-- because this would yield smaller lad at the canopy boundaries
	721	!-- than inside of the canopy.
	722	!-- For the same reason, the lad at the rightmost(i+1)canopy
	723	!-- boundary on the u-grid equals lad_s(k,j,i).
	724	lad_local = lad_s(k,j,i)
	725	IF ( lad_local == 0.0_wp .AND. &
	726	lad_s(k,j,i-1) > 0.0_wp ) THEN
	727	lad_local = lad_s(k,j,i-1)
	728	ENDIF
	729
	730	pre_tend = 0.0_wp
	731	pre_u = 0.0_wp
	732	!
	733	!-- Calculate preliminary value (pre_tend) of the tendency
	734	pre_tend = - cdc * &
	735	lad_local * &
	736	SQRT( u(k,j,i)**2 + &
	737	( 0.25_wp * ( v(k,j,i-1) + &
	738	v(k,j,i) + &
	739	v(k,j+1,i) + &
	740	v(k,j+1,i-1) ) &
	741	)**2 + &
	742	( 0.25_wp * ( w(k-1,j,i-1) + &
	743	w(k-1,j,i) + &
	744	w(k,j,i-1) + &
	745	w(k,j,i) ) &
	746	)**2 &
	747	) * &
	748	u(k,j,i)
	749
	750	!
	751	!-- Calculate preliminary new velocity, based on pre_tend
	752	pre_u = u(k,j,i) + dt_3d * pre_tend
	753	!
	754	!-- Compare sign of old velocity and new preliminary velocity,
	755	!-- and in case the signs are different, limit the tendency
	756	IF ( SIGN(pre_u,u(k,j,i)) /= pre_u ) THEN
	757	pre_tend = - u(k,j,i) * ddt_3d
	758	ELSE
	759	pre_tend = pre_tend
	760	ENDIF
	761	!
	762	!-- Calculate final tendency
	763	tend(k,j,i) = tend(k,j,i) + pre_tend
	764	ENDDO
	765
	766
	767	!
[142]	768	!-- v-component
[1484]	769	CASE ( 2 )
	770	DO k = nzb_v_inner(j,i)+1, pch_index
[138]	771
	772	!
[1484]	773	!-- In order to create sharp boundaries of the plant canopy,
	774	!-- the lad on the v-grid at index (k,j,i) is equal to lad_s(k,j,i),
	775	!-- rather than being interpolated from the surrounding lad_s,
	776	!-- because this would yield smaller lad at the canopy boundaries
	777	!-- than inside of the canopy.
	778	!-- For the same reason, the lad at the northmost(j+1)canopy
	779	!-- boundary on the v-grid equals lad_s(k,j,i).
	780	lad_local = lad_s(k,j,i)
	781	IF ( lad_local == 0.0_wp .AND. &
	782	lad_s(k,j-1,i) > 0.0_wp ) THEN
	783	lad_local = lad_s(k,j-1,i)
	784	ENDIF
	785
	786	pre_tend = 0.0_wp
	787	pre_v = 0.0_wp
	788	!
	789	!-- Calculate preliminary value (pre_tend) of the tendency
	790	pre_tend = - cdc * &
	791	lad_local * &
	792	SQRT( ( 0.25_wp * ( u(k,j-1,i) + &
	793	u(k,j-1,i+1) + &
	794	u(k,j,i) + &
	795	u(k,j,i+1) ) &
	796	)**2 + &
	797	v(k,j,i)**2 + &
	798	( 0.25_wp * ( w(k-1,j-1,i) + &
	799	w(k-1,j,i) + &
	800	w(k,j-1,i) + &
	801	w(k,j,i) ) &
	802	)**2 &
	803	) * &
	804	v(k,j,i)
	805
	806	!
	807	!-- Calculate preliminary new velocity, based on pre_tend
	808	pre_v = v(k,j,i) + dt_3d * pre_tend
	809	!
	810	!-- Compare sign of old velocity and new preliminary velocity,
	811	!-- and in case the signs are different, limit the tendency
	812	IF ( SIGN(pre_v,v(k,j,i)) /= pre_v ) THEN
	813	pre_tend = - v(k,j,i) * ddt_3d
	814	ELSE
	815	pre_tend = pre_tend
	816	ENDIF
	817	!
	818	!-- Calculate final tendency
	819	tend(k,j,i) = tend(k,j,i) + pre_tend
	820	ENDDO
	821
	822
	823	!
[142]	824	!-- w-component
[1484]	825	CASE ( 3 )
	826	DO k = nzb_w_inner(j,i)+1, pch_index-1
[138]	827
[1484]	828	pre_tend = 0.0_wp
	829	pre_w = 0.0_wp
[138]	830	!
[1484]	831	!-- Calculate preliminary value (pre_tend) of the tendency
	832	pre_tend = - cdc * &
	833	(0.5_wp * &
	834	( lad_s(k+1,j,i) + lad_s(k,j,i) )) * &
	835	SQRT( ( 0.25_wp * ( u(k,j,i) + &
	836	u(k,j,i+1) + &
	837	u(k+1,j,i) + &
	838	u(k+1,j,i+1) ) &
	839	)**2 + &
	840	( 0.25_wp * ( v(k,j,i) + &
	841	v(k,j+1,i) + &
	842	v(k+1,j,i) + &
	843	v(k+1,j+1,i) ) &
	844	)**2 + &
	845	w(k,j,i)**2 &
	846	) * &
	847	w(k,j,i)
	848	!
	849	!-- Calculate preliminary new velocity, based on pre_tend
	850	pre_w = w(k,j,i) + dt_3d * pre_tend
	851	!
	852	!-- Compare sign of old velocity and new preliminary velocity,
	853	!-- and in case the signs are different, limit the tendency
	854	IF ( SIGN(pre_w,w(k,j,i)) /= pre_w ) THEN
	855	pre_tend = - w(k,j,i) * ddt_3d
	856	ELSE
	857	pre_tend = pre_tend
	858	ENDIF
	859	!
	860	!-- Calculate final tendency
	861	tend(k,j,i) = tend(k,j,i) + pre_tend
	862	ENDDO
	863
	864	!
[153]	865	!-- potential temperature
	866	CASE ( 4 )
	867	DO k = nzb_s_inner(j,i)+1, pch_index
[1484]	868	tend(k,j,i) = tend(k,j,i) + &
	869	( canopy_heat_flux(k,j,i) - &
	870	canopy_heat_flux(k-1,j,i) ) / dzw(k)
[153]	871	ENDDO
	872
	873
	874	!
	875	!-- scalar concentration
	876	CASE ( 5 )
	877	DO k = nzb_s_inner(j,i)+1, pch_index
[1484]	878	tend(k,j,i) = tend(k,j,i) - &
	879	lsec * &
	880	lad_s(k,j,i) * &
	881	SQRT( ( 0.5_wp * ( u(k,j,i) + &
	882	u(k,j,i+1) ) &
	883	)**2 + &
	884	( 0.5_wp * ( v(k,j,i) + &
	885	v(k,j+1,i) ) &
	886	)**2 + &
	887	( 0.5_wp * ( w(k-1,j,i) + &
	888	w(k,j,i) ) &
	889	)**2 &
	890	) * &
	891	( q(k,j,i) - lsc )
[153]	892	ENDDO
	893
	894	!
[142]	895	!-- sgs-tke
[1484]	896	CASE ( 6 )
	897	DO k = nzb_s_inner(j,i)+1, pch_index
	898	tend(k,j,i) = tend(k,j,i) - &
	899	2.0_wp * cdc * &
	900	lad_s(k,j,i) * &
	901	SQRT( ( 0.5_wp * ( u(k,j,i) + &
	902	u(k,j,i+1) ) &
	903	)**2 + &
	904	( 0.5_wp * ( v(k,j,i) + &
	905	v(k,j+1,i) ) &
	906	)**2 + &
	907	( 0.5_wp * ( w(k,j,i) + &
	908	w(k+1,j,i) ) &
	909	)**2 &
	910	) * &
	911	e(k,j,i)
	912	ENDDO
[138]	913
[142]	914	CASE DEFAULT
[138]	915
[257]	916	WRITE( message_string, * ) 'wrong component: ', component
	917	CALL message( 'plant_canopy_model', 'PA0279', 1, 2, 0, 6, 0 )
[138]	918
[142]	919	END SELECT
[138]	920
	921	END SUBROUTINE plant_canopy_model_ij
	922
	923	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 |