Home

Context Navigation

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

Last change on this file since 1634 was 1341, checked in by kanani, 11 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 27.0 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |