Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |