Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |