Home

Context Navigation

source: palm/trunk/SOURCE/diffusion_u.f90 @ 1185

Last change on this file since 1185 was 1132, checked in by raasch, 12 years ago
r1028 documented
Property svn:keywords set to `Id`
File size: 24.4 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |