Home

Context Navigation

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

Last change on this file since 1025 was 1017, checked in by raasch, 12 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 23.4 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |