Home

Context Navigation

source: palm/trunk/SOURCE/diffusion_v.f90 @ 2005

Last change on this file since 2005 was 2001, checked in by knoop, 8 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 26.8 KB

Rev	Line
[1873]	1	!> @file diffusion_v.f90
[2000]	2	!------------------------------------------------------------------------------!
[1036]	3	! This file is part of PALM.
	4	!
[2000]	5	! PALM is free software: you can redistribute it and/or modify it under the
	6	! terms of the GNU General Public License as published by the Free Software
	7	! Foundation, either version 3 of the License, or (at your option) any later
	8	! version.
[1036]	9	!
	10	! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
	11	! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
	12	! A PARTICULAR PURPOSE. See the GNU General Public License for more details.
	13	!
	14	! You should have received a copy of the GNU General Public License along with
	15	! PALM. If not, see <http://www.gnu.org/licenses/>.
	16	!
[1818]	17	! Copyright 1997-2016 Leibniz Universitaet Hannover
[2000]	18	!------------------------------------------------------------------------------!
[1036]	19	!
[484]	20	! Current revisions:
[1]	21	! -----------------
[1341]	22	!
[2001]	23	!
[1321]	24	! Former revisions:
	25	! -----------------
	26	! $Id: diffusion_v.f90 2001 2016-08-20 18:41:22Z suehring $
	27	!
[2001]	28	! 2000 2016-08-20 18:09:15Z knoop
	29	! Forced header and separation lines into 80 columns
	30	!
[1874]	31	! 1873 2016-04-18 14:50:06Z maronga
	32	! Module renamed (removed _mod)
	33	!
	34	!
[1851]	35	! 1850 2016-04-08 13:29:27Z maronga
	36	! Module renamed
	37	!
	38	!
[1741]	39	! 1740 2016-01-13 08:19:40Z raasch
	40	! unnecessary calculations of kmzm and kmzp in wall bounded parts removed
	41	!
[1683]	42	! 1682 2015-10-07 23:56:08Z knoop
	43	! Code annotations made doxygen readable
	44	!
[1341]	45	! 1340 2014-03-25 19:45:13Z kanani
	46	! REAL constants defined as wp-kind
	47	!
[1321]	48	! 1320 2014-03-20 08:40:49Z raasch
[1320]	49	! ONLY-attribute added to USE-statements,
	50	! kind-parameters added to all INTEGER and REAL declaration statements,
	51	! kinds are defined in new module kinds,
	52	! revision history before 2012 removed,
	53	! comment fields (!:) to be used for variable explanations added to
	54	! all variable declaration statements
[1321]	55	!
[1258]	56	! 1257 2013-11-08 15:18:40Z raasch
	57	! openacc loop and loop vector clauses removed, declare create moved after
	58	! the FORTRAN declaration statement
	59	!
[1132]	60	! 1128 2013-04-12 06:19:32Z raasch
	61	! loop index bounds in accelerator version replaced by i_left, i_right, j_south,
	62	! j_north
	63	!
[1037]	64	! 1036 2012-10-22 13:43:42Z raasch
	65	! code put under GPL (PALM 3.9)
	66	!
[1017]	67	! 1015 2012-09-27 09:23:24Z raasch
	68	! accelerator version (*_acc) added
	69	!
[1002]	70	! 1001 2012-09-13 14:08:46Z raasch
	71	! arrays comunicated by module instead of parameter list
	72	!
[979]	73	! 978 2012-08-09 08:28:32Z fricke
	74	! outflow damping layer removed
	75	! kmxm_x/_y and kmxp_x/_y change to kmxm and kmxp
	76	!
[1]	77	! Revision 1.1 1997/09/12 06:24:01 raasch
	78	! Initial revision
	79	!
	80	!
	81	! Description:
	82	! ------------
[1682]	83	!> Diffusion term of the v-component
[1]	84	!------------------------------------------------------------------------------!
[1682]	85	MODULE diffusion_v_mod
	86
[1]	87
[56]	88	USE wall_fluxes_mod
	89
[1]	90	PRIVATE
[1015]	91	PUBLIC diffusion_v, diffusion_v_acc
[1]	92
	93	INTERFACE diffusion_v
	94	MODULE PROCEDURE diffusion_v
	95	MODULE PROCEDURE diffusion_v_ij
	96	END INTERFACE diffusion_v
	97
[1015]	98	INTERFACE diffusion_v_acc
	99	MODULE PROCEDURE diffusion_v_acc
	100	END INTERFACE diffusion_v_acc
	101
[1]	102	CONTAINS
	103
	104
	105	!------------------------------------------------------------------------------!
[1682]	106	! Description:
	107	! ------------
	108	!> Call for all grid points
[1]	109	!------------------------------------------------------------------------------!
[1001]	110	SUBROUTINE diffusion_v
[1]	111
[1320]	112	USE arrays_3d, &
	113	ONLY: ddzu, ddzw, km, tend, u, v, vsws, vswst, w
	114
	115	USE control_parameters, &
	116	ONLY: constant_top_momentumflux, topography, use_surface_fluxes, &
	117	use_top_fluxes
	118
	119	USE grid_variables, &
	120	ONLY: ddx, ddy, ddy2, fxm, fxp, wall_v
	121
	122	USE indices, &
	123	ONLY: nxl, nxr, nyn, nys, nysv, nzb, nzb_diff_v, nzb_v_inner, &
	124	nzb_v_outer, nzt, nzt_diff
	125
	126	USE kinds
[1]	127
	128	IMPLICIT NONE
	129
[1682]	130	INTEGER(iwp) :: i !<
	131	INTEGER(iwp) :: j !<
	132	INTEGER(iwp) :: k !<
	133	REAL(wp) :: kmxm !<
	134	REAL(wp) :: kmxp !<
	135	REAL(wp) :: kmzm !<
	136	REAL(wp) :: kmzp !<
[1001]	137
[1682]	138	REAL(wp), DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) :: vsus !<
[1]	139
[56]	140	!
	141	!-- First calculate horizontal momentum flux v'u' at vertical walls,
	142	!-- if neccessary
	143	IF ( topography /= 'flat' ) THEN
[1320]	144	CALL wall_fluxes( vsus, 0.0_wp, 1.0_wp, 0.0_wp, 0.0_wp, nzb_v_inner, &
[56]	145	nzb_v_outer, wall_v )
	146	ENDIF
	147
[1]	148	DO i = nxl, nxr
[106]	149	DO j = nysv, nyn
[1]	150	!
	151	!-- Compute horizontal diffusion
	152	DO k = nzb_v_outer(j,i)+1, nzt
	153	!
	154	!-- Interpolate eddy diffusivities on staggered gridpoints
[1340]	155	kmxp = 0.25_wp * &
[978]	156	( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) )
[1340]	157	kmxm = 0.25_wp * &
[978]	158	( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) )
[1]	159
[1320]	160	tend(k,j,i) = tend(k,j,i) &
	161	& + ( kmxp * ( v(k,j,i+1) - v(k,j,i) ) * ddx &
	162	& + kmxp * ( u(k,j,i+1) - u(k,j-1,i+1) ) * ddy &
	163	& - kmxm * ( v(k,j,i) - v(k,j,i-1) ) * ddx &
	164	& - kmxm * ( u(k,j,i) - u(k,j-1,i) ) * ddy &
	165	& ) * ddx &
[1340]	166	& + 2.0_wp * ( &
[1320]	167	& km(k,j,i) * ( v(k,j+1,i) - v(k,j,i) ) &
	168	& - km(k,j-1,i) * ( v(k,j,i) - v(k,j-1,i) ) &
[1340]	169	& ) * ddy2
[1]	170	ENDDO
	171
	172	!
	173	!-- Wall functions at the left and right walls, respectively
[1340]	174	IF ( wall_v(j,i) /= 0.0_wp ) THEN
[51]	175
[1]	176	DO k = nzb_v_inner(j,i)+1, nzb_v_outer(j,i)
[1340]	177	kmxp = 0.25_wp * &
[978]	178	( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) )
[1340]	179	kmxm = 0.25_wp * &
[978]	180	( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) )
	181
[1]	182	tend(k,j,i) = tend(k,j,i) &
[1340]	183	+ 2.0_wp * ( &
[1]	184	km(k,j,i) * ( v(k,j+1,i) - v(k,j,i) ) &
	185	- km(k,j-1,i) * ( v(k,j,i) - v(k,j-1,i) ) &
[1340]	186	) * ddy2 &
[1]	187	+ ( fxp(j,i) * ( &
[978]	188	kmxp * ( v(k,j,i+1) - v(k,j,i) ) * ddx &
	189	+ kmxp * ( u(k,j,i+1) - u(k,j-1,i+1) ) * ddy &
[1]	190	) &
	191	- fxm(j,i) * ( &
[978]	192	kmxm * ( v(k,j,i) - v(k,j,i-1) ) * ddx &
	193	+ kmxm * ( u(k,j,i) - u(k,j-1,i) ) * ddy &
[1]	194	) &
[56]	195	+ wall_v(j,i) * vsus(k,j,i) &
[1]	196	) * ddx
	197	ENDDO
	198	ENDIF
	199
	200	!
	201	!-- Compute vertical diffusion. In case of simulating a Prandtl
	202	!-- layer, index k starts at nzb_v_inner+2.
[102]	203	DO k = nzb_diff_v(j,i), nzt_diff
[1]	204	!
	205	!-- Interpolate eddy diffusivities on staggered gridpoints
[1340]	206	kmzp = 0.25_wp * &
[1]	207	( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
[1340]	208	kmzm = 0.25_wp * &
[1]	209	( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) )
	210
[1320]	211	tend(k,j,i) = tend(k,j,i) &
	212	& + ( kmzp * ( ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1) &
	213	& + ( w(k,j,i) - w(k,j-1,i) ) * ddy &
	214	& ) &
	215	& - kmzm * ( ( v(k,j,i) - v(k-1,j,i) ) * ddzu(k) &
	216	& + ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy &
	217	& ) &
[1]	218	& ) * ddzw(k)
	219	ENDDO
	220
	221	!
	222	!-- Vertical diffusion at the first grid point above the surface,
	223	!-- if the momentum flux at the bottom is given by the Prandtl law
	224	!-- or if it is prescribed by the user.
	225	!-- Difference quotient of the momentum flux is not formed over
	226	!-- half of the grid spacing (2.0*ddzw(k)) any more, since the
[1320]	227	!-- comparison with other (LES) models showed that the values of
[1]	228	!-- the momentum flux becomes too large in this case.
	229	!-- The term containing w(k-1,..) (see above equation) is removed here
	230	!-- because the vertical velocity is assumed to be zero at the surface.
	231	IF ( use_surface_fluxes ) THEN
	232	k = nzb_v_inner(j,i)+1
	233	!
	234	!-- Interpolate eddy diffusivities on staggered gridpoints
[1340]	235	kmzp = 0.25_wp * &
[1]	236	( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
	237
[1320]	238	tend(k,j,i) = tend(k,j,i) &
	239	& + ( kmzp * ( w(k,j,i) - w(k,j-1,i) ) * ddy &
	240	& ) * ddzw(k) &
	241	& + ( kmzp * ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1) &
	242	& + vsws(j,i) &
[1]	243	& ) * ddzw(k)
	244	ENDIF
	245
[102]	246	!
	247	!-- Vertical diffusion at the first gridpoint below the top boundary,
	248	!-- if the momentum flux at the top is prescribed by the user
[103]	249	IF ( use_top_fluxes .AND. constant_top_momentumflux ) THEN
[102]	250	k = nzt
	251	!
	252	!-- Interpolate eddy diffusivities on staggered gridpoints
[1340]	253	kmzm = 0.25_wp * &
[102]	254	( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) )
	255
[1320]	256	tend(k,j,i) = tend(k,j,i) &
	257	& - ( kmzm * ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy &
	258	& ) * ddzw(k) &
	259	& + ( -vswst(j,i) &
	260	& - kmzm * ( v(k,j,i) - v(k-1,j,i) ) * ddzu(k) &
[102]	261	& ) * ddzw(k)
	262	ENDIF
	263
[1]	264	ENDDO
	265	ENDDO
	266
	267	END SUBROUTINE diffusion_v
	268
	269
	270	!------------------------------------------------------------------------------!
[1682]	271	! Description:
	272	! ------------
	273	!> Call for all grid points - accelerator version
[1015]	274	!------------------------------------------------------------------------------!
	275	SUBROUTINE diffusion_v_acc
	276
[1320]	277	USE arrays_3d, &
	278	ONLY: ddzu, ddzw, km, tend, u, v, vsws, vswst, w
	279
	280	USE control_parameters, &
	281	ONLY: constant_top_momentumflux, topography, use_surface_fluxes, &
	282	use_top_fluxes
	283
	284	USE grid_variables, &
	285	ONLY: ddx, ddy, ddy2, fxm, fxp, wall_v
	286
	287	USE indices, &
	288	ONLY: i_left, i_right, j_north, j_south, nxl, nxr, nyn, nys, nzb, &
	289	nzb_diff_v, nzb_v_inner, nzb_v_outer, nzt, nzt_diff
	290
	291	USE kinds
[1015]	292
	293	IMPLICIT NONE
	294
[1682]	295	INTEGER(iwp) :: i !<
	296	INTEGER(iwp) :: j !<
	297	INTEGER(iwp) :: k !<
	298	REAL(wp) :: kmxm !<
	299	REAL(wp) :: kmxp !<
	300	REAL(wp) :: kmzm !<
	301	REAL(wp) :: kmzp !<
[1015]	302
[1682]	303	REAL(wp), DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) :: vsus !<
[1015]	304	!$acc declare create ( vsus )
	305
	306	!
	307	!-- First calculate horizontal momentum flux v'u' at vertical walls,
	308	!-- if neccessary
	309	IF ( topography /= 'flat' ) THEN
[1320]	310	CALL wall_fluxes_acc( vsus, 0.0_wp, 1.0_wp, 0.0_wp, 0.0_wp, &
	311	nzb_v_inner, nzb_v_outer, wall_v )
[1015]	312	ENDIF
	313
[1320]	314	!$acc kernels present ( u, v, w, km, tend, vsws, vswst ) &
	315	!$acc present ( ddzu, ddzw, fxm, fxp, wall_v ) &
[1015]	316	!$acc present ( nzb_v_inner, nzb_v_outer, nzb_diff_v )
[1128]	317	DO i = i_left, i_right
	318	DO j = j_south, j_north
[1015]	319	!
	320	!-- Compute horizontal diffusion
	321	DO k = 1, nzt
	322	IF ( k > nzb_v_outer(j,i) ) THEN
	323	!
	324	!-- Interpolate eddy diffusivities on staggered gridpoints
[1340]	325	kmxp = 0.25_wp * &
[1015]	326	( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) )
[1340]	327	kmxm = 0.25_wp * &
[1015]	328	( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) )
	329
	330	tend(k,j,i) = tend(k,j,i) &
	331	& + ( kmxp * ( v(k,j,i+1) - v(k,j,i) ) * ddx &
	332	& + kmxp * ( u(k,j,i+1) - u(k,j-1,i+1) ) * ddy &
	333	& - kmxm * ( v(k,j,i) - v(k,j,i-1) ) * ddx &
	334	& - kmxm * ( u(k,j,i) - u(k,j-1,i) ) * ddy &
	335	& ) * ddx &
[1340]	336	& + 2.0_wp * ( &
[1015]	337	& km(k,j,i) * ( v(k,j+1,i) - v(k,j,i) ) &
	338	& - km(k,j-1,i) * ( v(k,j,i) - v(k,j-1,i) ) &
[1340]	339	& ) * ddy2
[1015]	340	ENDIF
	341	ENDDO
	342
	343	!
	344	!-- Wall functions at the left and right walls, respectively
	345	DO k = 1, nzt
[1320]	346	IF( k > nzb_v_inner(j,i) .AND. k <= nzb_v_outer(j,i) .AND. &
[1340]	347	wall_v(j,i) /= 0.0_wp ) THEN
[1015]	348
[1340]	349	kmxp = 0.25_wp * &
[1015]	350	( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) )
[1340]	351	kmxm = 0.25_wp * &
[1015]	352	( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) )
	353
	354	tend(k,j,i) = tend(k,j,i) &
[1340]	355	+ 2.0_wp * ( &
[1015]	356	km(k,j,i) * ( v(k,j+1,i) - v(k,j,i) ) &
	357	- km(k,j-1,i) * ( v(k,j,i) - v(k,j-1,i) ) &
[1340]	358	) * ddy2 &
[1015]	359	+ ( fxp(j,i) * ( &
	360	kmxp * ( v(k,j,i+1) - v(k,j,i) ) * ddx &
	361	+ kmxp * ( u(k,j,i+1) - u(k,j-1,i+1) ) * ddy &
	362	) &
	363	- fxm(j,i) * ( &
	364	kmxm * ( v(k,j,i) - v(k,j,i-1) ) * ddx &
	365	+ kmxm * ( u(k,j,i) - u(k,j-1,i) ) * ddy &
	366	) &
	367	+ wall_v(j,i) * vsus(k,j,i) &
	368	) * ddx
	369	ENDIF
	370	ENDDO
	371
	372	!
	373	!-- Compute vertical diffusion. In case of simulating a Prandtl
	374	!-- layer, index k starts at nzb_v_inner+2.
	375	DO k = 1, nzt_diff
	376	IF ( k >= nzb_diff_v(j,i) ) THEN
	377	!
	378	!-- Interpolate eddy diffusivities on staggered gridpoints
[1340]	379	kmzp = 0.25_wp * &
[1015]	380	( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
[1340]	381	kmzm = 0.25_wp * &
[1015]	382	( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) )
	383
	384	tend(k,j,i) = tend(k,j,i) &
	385	& + ( kmzp * ( ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1)&
	386	& + ( w(k,j,i) - w(k,j-1,i) ) * ddy &
	387	& ) &
	388	& - kmzm * ( ( v(k,j,i) - v(k-1,j,i) ) * ddzu(k)&
	389	& + ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy &
	390	& ) &
	391	& ) * ddzw(k)
	392	ENDIF
	393	ENDDO
	394
	395	ENDDO
	396	ENDDO
	397
	398	!
	399	!-- Vertical diffusion at the first grid point above the surface,
	400	!-- if the momentum flux at the bottom is given by the Prandtl law
	401	!-- or if it is prescribed by the user.
	402	!-- Difference quotient of the momentum flux is not formed over
	403	!-- half of the grid spacing (2.0*ddzw(k)) any more, since the
[1320]	404	!-- comparison with other (LES) models showed that the values of
[1015]	405	!-- the momentum flux becomes too large in this case.
	406	!-- The term containing w(k-1,..) (see above equation) is removed here
	407	!-- because the vertical velocity is assumed to be zero at the surface.
	408	IF ( use_surface_fluxes ) THEN
	409
[1128]	410	DO i = i_left, i_right
	411	DO j = j_south, j_north
[1015]	412
	413	k = nzb_v_inner(j,i)+1
	414	!
	415	!-- Interpolate eddy diffusivities on staggered gridpoints
[1340]	416	kmzp = 0.25_wp * &
[1015]	417	( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
	418
[1320]	419	tend(k,j,i) = tend(k,j,i) &
	420	& + ( kmzp * ( w(k,j,i) - w(k,j-1,i) ) * ddy &
	421	& ) * ddzw(k) &
	422	& + ( kmzp * ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1) &
	423	& + vsws(j,i) &
[1015]	424	& ) * ddzw(k)
	425	ENDDO
	426	ENDDO
	427
	428	ENDIF
	429
	430	!
	431	!-- Vertical diffusion at the first gridpoint below the top boundary,
	432	!-- if the momentum flux at the top is prescribed by the user
	433	IF ( use_top_fluxes .AND. constant_top_momentumflux ) THEN
	434
	435	k = nzt
	436
[1128]	437	DO i = i_left, i_right
	438	DO j = j_south, j_north
[1015]	439
	440	!
	441	!-- Interpolate eddy diffusivities on staggered gridpoints
[1340]	442	kmzm = 0.25_wp * &
[1015]	443	( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) )
	444
[1320]	445	tend(k,j,i) = tend(k,j,i) &
	446	& - ( kmzm * ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy &
	447	& ) * ddzw(k) &
	448	& + ( -vswst(j,i) &
	449	& - kmzm * ( v(k,j,i) - v(k-1,j,i) ) * ddzu(k) &
[1015]	450	& ) * ddzw(k)
	451	ENDDO
	452	ENDDO
	453
	454	ENDIF
	455	!$acc end kernels
	456
	457	END SUBROUTINE diffusion_v_acc
	458
	459
	460	!------------------------------------------------------------------------------!
[1682]	461	! Description:
	462	! ------------
	463	!> Call for grid point i,j
[1]	464	!------------------------------------------------------------------------------!
[1001]	465	SUBROUTINE diffusion_v_ij( i, j )
[1]	466
[1320]	467	USE arrays_3d, &
	468	ONLY: ddzu, ddzw, km, tend, u, v, vsws, vswst, w
	469
	470	USE control_parameters, &
	471	ONLY: constant_top_momentumflux, use_surface_fluxes, use_top_fluxes
	472
	473	USE grid_variables, &
	474	ONLY: ddx, ddy, ddy2, fxm, fxp, wall_v
	475
	476	USE indices, &
	477	ONLY: nzb, nzb_diff_v, nzb_v_inner, nzb_v_outer, nzt, nzt_diff
	478
	479	USE kinds
[1]	480
	481	IMPLICIT NONE
	482
[1682]	483	INTEGER(iwp) :: i !<
	484	INTEGER(iwp) :: j !<
	485	INTEGER(iwp) :: k !<
	486	REAL(wp) :: kmxm !<
	487	REAL(wp) :: kmxp !<
	488	REAL(wp) :: kmzm !<
	489	REAL(wp) :: kmzp !<
[1]	490
[1682]	491	REAL(wp), DIMENSION(nzb:nzt+1) :: vsus !<
[1001]	492
[1]	493	!
	494	!-- Compute horizontal diffusion
	495	DO k = nzb_v_outer(j,i)+1, nzt
	496	!
	497	!-- Interpolate eddy diffusivities on staggered gridpoints
[1340]	498	kmxp = 0.25_wp * ( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) )
	499	kmxm = 0.25_wp * ( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) )
[1]	500
[1320]	501	tend(k,j,i) = tend(k,j,i) &
	502	& + ( kmxp * ( v(k,j,i+1) - v(k,j,i) ) * ddx &
	503	& + kmxp * ( u(k,j,i+1) - u(k,j-1,i+1) ) * ddy &
	504	& - kmxm * ( v(k,j,i) - v(k,j,i-1) ) * ddx &
	505	& - kmxm * ( u(k,j,i) - u(k,j-1,i) ) * ddy &
	506	& ) * ddx &
[1340]	507	& + 2.0_wp * ( &
[1320]	508	& km(k,j,i) * ( v(k,j+1,i) - v(k,j,i) ) &
	509	& - km(k,j-1,i) * ( v(k,j,i) - v(k,j-1,i) ) &
[1340]	510	& ) * ddy2
[1]	511	ENDDO
	512
	513	!
	514	!-- Wall functions at the left and right walls, respectively
[1340]	515	IF ( wall_v(j,i) /= 0.0_wp ) THEN
[51]	516
	517	!
	518	!-- Calculate the horizontal momentum flux v'u'
[1320]	519	CALL wall_fluxes( i, j, nzb_v_inner(j,i)+1, nzb_v_outer(j,i), &
	520	vsus, 0.0_wp, 1.0_wp, 0.0_wp, 0.0_wp )
[51]	521
[1]	522	DO k = nzb_v_inner(j,i)+1, nzb_v_outer(j,i)
[1340]	523	kmxp = 0.25_wp * &
[978]	524	( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) )
[1340]	525	kmxm = 0.25_wp * &
[978]	526	( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) )
[1]	527
	528	tend(k,j,i) = tend(k,j,i) &
[1340]	529	+ 2.0_wp * ( &
[1]	530	km(k,j,i) * ( v(k,j+1,i) - v(k,j,i) ) &
	531	- km(k,j-1,i) * ( v(k,j,i) - v(k,j-1,i) ) &
[1340]	532	) * ddy2 &
[1]	533	+ ( fxp(j,i) * ( &
[978]	534	kmxp * ( v(k,j,i+1) - v(k,j,i) ) * ddx &
	535	+ kmxp * ( u(k,j,i+1) - u(k,j-1,i+1) ) * ddy &
[1]	536	) &
	537	- fxm(j,i) * ( &
[978]	538	kmxm * ( v(k,j,i) - v(k,j,i-1) ) * ddx &
	539	+ kmxm * ( u(k,j,i) - u(k,j-1,i) ) * ddy &
[1]	540	) &
[51]	541	+ wall_v(j,i) * vsus(k) &
[1]	542	) * ddx
	543	ENDDO
	544	ENDIF
	545
	546	!
	547	!-- Compute vertical diffusion. In case of simulating a Prandtl layer,
	548	!-- index k starts at nzb_v_inner+2.
[102]	549	DO k = nzb_diff_v(j,i), nzt_diff
[1]	550	!
	551	!-- Interpolate eddy diffusivities on staggered gridpoints
[1340]	552	kmzp = 0.25_wp * ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
	553	kmzm = 0.25_wp * ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) )
[1]	554
[1320]	555	tend(k,j,i) = tend(k,j,i) &
	556	& + ( kmzp * ( ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1) &
	557	& + ( w(k,j,i) - w(k,j-1,i) ) * ddy &
	558	& ) &
	559	& - kmzm * ( ( v(k,j,i) - v(k-1,j,i) ) * ddzu(k) &
	560	& + ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy &
	561	& ) &
[1]	562	& ) * ddzw(k)
	563	ENDDO
	564
	565	!
	566	!-- Vertical diffusion at the first grid point above the surface, if the
	567	!-- momentum flux at the bottom is given by the Prandtl law or if it is
	568	!-- prescribed by the user.
	569	!-- Difference quotient of the momentum flux is not formed over half of
	570	!-- the grid spacing (2.0*ddzw(k)) any more, since the comparison with
[1320]	571	!-- other (LES) models showed that the values of the momentum flux becomes
[1]	572	!-- too large in this case.
	573	!-- The term containing w(k-1,..) (see above equation) is removed here
	574	!-- because the vertical velocity is assumed to be zero at the surface.
	575	IF ( use_surface_fluxes ) THEN
	576	k = nzb_v_inner(j,i)+1
	577	!
	578	!-- Interpolate eddy diffusivities on staggered gridpoints
[1340]	579	kmzp = 0.25_wp * ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
[1]	580
[1320]	581	tend(k,j,i) = tend(k,j,i) &
	582	& + ( kmzp * ( w(k,j,i) - w(k,j-1,i) ) * ddy &
	583	& ) * ddzw(k) &
	584	& + ( kmzp * ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1) &
	585	& + vsws(j,i) &
[1]	586	& ) * ddzw(k)
	587	ENDIF
	588
[102]	589	!
	590	!-- Vertical diffusion at the first gridpoint below the top boundary,
	591	!-- if the momentum flux at the top is prescribed by the user
[103]	592	IF ( use_top_fluxes .AND. constant_top_momentumflux ) THEN
[102]	593	k = nzt
	594	!
	595	!-- Interpolate eddy diffusivities on staggered gridpoints
[1340]	596	kmzm = 0.25_wp * &
[102]	597	( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) )
	598
[1320]	599	tend(k,j,i) = tend(k,j,i) &
	600	& - ( kmzm * ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy &
	601	& ) * ddzw(k) &
	602	& + ( -vswst(j,i) &
	603	& - kmzm * ( v(k,j,i) - v(k-1,j,i) ) * ddzu(k) &
[102]	604	& ) * ddzw(k)
	605	ENDIF
	606
[1]	607	END SUBROUTINE diffusion_v_ij
	608
[1321]	609	END MODULE diffusion_v_mod

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |