Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |