Home

Context Navigation

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

Last change on this file since 1120 was 1037, checked in by raasch, 12 years ago
last commit documented
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	!
	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 1037 2012-10-22 14:10:22Z raasch $
[39]	27	!
[1037]	28	! 1036 2012-10-22 13:43:42Z raasch
	29	! code put under GPL (PALM 3.9)
	30	!
[1017]	31	! 1015 2012-09-27 09:23:24Z raasch
	32	! accelerator version (*_acc) added
	33	!
[1002]	34	! 1001 2012-09-13 14:08:46Z raasch
	35	! arrays comunicated by module instead of parameter list
	36	!
[979]	37	! 978 2012-08-09 08:28:32Z fricke
	38	! outflow damping layer removed
	39	! kmym_x/_y and kmyp_x/_y change to kmym and kmyp
	40	!
[668]	41	! 667 2010-12-23 12:06:00Z suehring/gryschka
	42	! nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng
	43	!
[392]	44	! 366 2009-08-25 08:06:27Z raasch
	45	! bc_ns replaced by bc_ns_cyc
	46	!
[110]	47	! 106 2007-08-16 14:30:26Z raasch
	48	! Momentumflux at top (uswst) included as boundary condition,
	49	! i loop is starting from nxlu (needed for non-cyclic boundary conditions)
	50	!
[77]	51	! 75 2007-03-22 09:54:05Z raasch
	52	! Wall functions now include diabatic conditions, call of routine wall_fluxes,
	53	! z0 removed from argument list, uxrp eliminated
	54	!
[39]	55	! 20 2007-02-26 00:12:32Z raasch
	56	! Bugfix: ddzw dimensioned 1:nzt"+1"
	57	!
[3]	58	! RCS Log replace by Id keyword, revision history cleaned up
	59	!
[1]	60	! Revision 1.15 2006/02/23 10:35:35 raasch
	61	! nzb_2d replaced by nzb_u_outer in horizontal diffusion and by nzb_u_inner
	62	! or nzb_diff_u, respectively, in vertical diffusion,
	63	! wall functions added for north and south walls, +z0 in argument list,
	64	! terms containing w(k-1,..) are removed from the Prandtl-layer equation
	65	! because they cause errors at the edges of topography
	66	! WARNING: loops containing the MAX function are still not properly vectorized!
	67	!
	68	! Revision 1.1 1997/09/12 06:23:51 raasch
	69	! Initial revision
	70	!
	71	!
	72	! Description:
	73	! ------------
	74	! Diffusion term of the u-component
[51]	75	! To do: additional damping (needed for non-cyclic bc) causes bad vectorization
	76	! and slows down the speed on NEC about 5-10%
[1]	77	!------------------------------------------------------------------------------!
	78
[56]	79	USE wall_fluxes_mod
	80
[1]	81	PRIVATE
[1015]	82	PUBLIC diffusion_u, diffusion_u_acc
[1]	83
	84	INTERFACE diffusion_u
	85	MODULE PROCEDURE diffusion_u
	86	MODULE PROCEDURE diffusion_u_ij
	87	END INTERFACE diffusion_u
	88
[1015]	89	INTERFACE diffusion_u_acc
	90	MODULE PROCEDURE diffusion_u_acc
	91	END INTERFACE diffusion_u_acc
	92
[1]	93	CONTAINS
	94
	95
	96	!------------------------------------------------------------------------------!
	97	! Call for all grid points
	98	!------------------------------------------------------------------------------!
[1001]	99	SUBROUTINE diffusion_u
[1]	100
[1001]	101	USE arrays_3d
[1]	102	USE control_parameters
	103	USE grid_variables
	104	USE indices
	105
	106	IMPLICIT NONE
	107
	108	INTEGER :: i, j, k
[978]	109	REAL :: kmym, kmyp, kmzm, kmzp
[1001]	110
[75]	111	REAL, DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) :: usvs
[1]	112
[56]	113	!
	114	!-- First calculate horizontal momentum flux u'v' at vertical walls,
	115	!-- if neccessary
	116	IF ( topography /= 'flat' ) THEN
[75]	117	CALL wall_fluxes( usvs, 1.0, 0.0, 0.0, 0.0, nzb_u_inner, &
[56]	118	nzb_u_outer, wall_u )
	119	ENDIF
	120
[106]	121	DO i = nxlu, nxr
[1001]	122	DO j = nys, nyn
[1]	123	!
	124	!-- Compute horizontal diffusion
	125	DO k = nzb_u_outer(j,i)+1, nzt
	126	!
	127	!-- Interpolate eddy diffusivities on staggered gridpoints
[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 &
[978]	138	& + ( kmyp * ( u(k,j+1,i) - u(k,j,i) ) * ddy &
	139	& + kmyp * ( v(k,j+1,i) - v(k,j+1,i-1) ) * ddx &
	140	& - kmym * ( u(k,j,i) - u(k,j-1,i) ) * ddy &
	141	& - kmym * ( v(k,j,i) - v(k,j,i-1) ) * ddx &
[1]	142	& ) * ddy
	143	ENDDO
	144
	145	!
	146	!-- Wall functions at the north and south walls, respectively
	147	IF ( wall_u(j,i) /= 0.0 ) THEN
[51]	148
[1]	149	DO k = nzb_u_inner(j,i)+1, nzb_u_outer(j,i)
[978]	150	kmyp = 0.25 * &
	151	( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) )
	152	kmym = 0.25 * &
	153	( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) )
[1]	154
	155	tend(k,j,i) = tend(k,j,i) &
	156	+ 2.0 * ( &
	157	km(k,j,i) * ( u(k,j,i+1) - u(k,j,i) ) &
	158	- km(k,j,i-1) * ( u(k,j,i) - u(k,j,i-1) ) &
	159	) * ddx2 &
	160	+ ( fyp(j,i) * ( &
[978]	161	kmyp * ( u(k,j+1,i) - u(k,j,i) ) * ddy &
	162	+ kmyp * ( v(k,j+1,i) - v(k,j+1,i-1) ) * ddx &
[1]	163	) &
	164	- fym(j,i) * ( &
[978]	165	kmym * ( u(k,j,i) - u(k,j-1,i) ) * ddy &
	166	+ kmym * ( v(k,j,i) - v(k,j,i-1) ) * ddx &
[1]	167	) &
[56]	168	+ wall_u(j,i) * usvs(k,j,i) &
[1]	169	) * ddy
	170	ENDDO
	171	ENDIF
	172
	173	!
	174	!-- Compute vertical diffusion. In case of simulating a Prandtl layer,
	175	!-- index k starts at nzb_u_inner+2.
[102]	176	DO k = nzb_diff_u(j,i), nzt_diff
[1]	177	!
	178	!-- Interpolate eddy diffusivities on staggered gridpoints
	179	kmzp = 0.25 * &
	180	( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) )
	181	kmzm = 0.25 * &
	182	( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
	183
	184	tend(k,j,i) = tend(k,j,i) &
	185	& + ( kmzp * ( ( u(k+1,j,i) - u(k,j,i) ) * ddzu(k+1) &
	186	& + ( w(k,j,i) - w(k,j,i-1) ) * ddx &
	187	& ) &
	188	& - kmzm * ( ( u(k,j,i) - u(k-1,j,i) ) * ddzu(k) &
	189	& + ( w(k-1,j,i) - w(k-1,j,i-1) ) * ddx &
[667]	190	& ) &
[1]	191	& ) * ddzw(k)
	192	ENDDO
	193
	194	!
	195	!-- Vertical diffusion at the first grid point above the surface,
	196	!-- if the momentum flux at the bottom is given by the Prandtl law or
	197	!-- if it is prescribed by the user.
	198	!-- Difference quotient of the momentum flux is not formed over half
	199	!-- of the grid spacing (2.0*ddzw(k)) any more, since the comparison
	200	!-- with other (LES) modell showed that the values of the momentum
	201	!-- flux becomes too large in this case.
	202	!-- The term containing w(k-1,..) (see above equation) is removed here
	203	!-- because the vertical velocity is assumed to be zero at the surface.
	204	IF ( use_surface_fluxes ) THEN
	205	k = nzb_u_inner(j,i)+1
	206	!
	207	!-- Interpolate eddy diffusivities on staggered gridpoints
	208	kmzp = 0.25 * &
	209	( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) )
	210	kmzm = 0.25 * &
	211	( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
	212
	213	tend(k,j,i) = tend(k,j,i) &
	214	& + ( kmzp * ( w(k,j,i) - w(k,j,i-1) ) * ddx &
	215	& ) * ddzw(k) &
[102]	216	& + ( kmzp * ( u(k+1,j,i) - u(k,j,i) ) * ddzu(k+1) &
[1]	217	& + usws(j,i) &
	218	& ) * ddzw(k)
	219	ENDIF
	220
[102]	221	!
	222	!-- Vertical diffusion at the first gridpoint below the top boundary,
	223	!-- if the momentum flux at the top is prescribed by the user
[103]	224	IF ( use_top_fluxes .AND. constant_top_momentumflux ) THEN
[102]	225	k = nzt
	226	!
	227	!-- Interpolate eddy diffusivities on staggered gridpoints
	228	kmzp = 0.25 * &
	229	( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) )
	230	kmzm = 0.25 * &
	231	( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
	232
	233	tend(k,j,i) = tend(k,j,i) &
	234	& - ( kmzm * ( w(k-1,j,i) - w(k-1,j,i-1) ) * ddx &
	235	& ) * ddzw(k) &
	236	& + ( -uswst(j,i) &
	237	& - kmzm * ( u(k,j,i) - u(k-1,j,i) ) * ddzu(k) &
	238	& ) * ddzw(k)
	239	ENDIF
	240
[1]	241	ENDDO
	242	ENDDO
	243
	244	END SUBROUTINE diffusion_u
	245
	246
	247	!------------------------------------------------------------------------------!
[1015]	248	! Call for all grid points - accelerator version
	249	!------------------------------------------------------------------------------!
	250	SUBROUTINE diffusion_u_acc
	251
	252	USE arrays_3d
	253	USE control_parameters
	254	USE grid_variables
	255	USE indices
	256
	257	IMPLICIT NONE
	258
	259	INTEGER :: i, j, k
	260	REAL :: kmym, kmyp, kmzm, kmzp
	261
	262	!$acc declare create ( usvs )
	263	REAL, DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) :: usvs
	264
	265	!
	266	!-- First calculate horizontal momentum flux u'v' at vertical walls,
	267	!-- if neccessary
	268	IF ( topography /= 'flat' ) THEN
	269	CALL wall_fluxes_acc( usvs, 1.0, 0.0, 0.0, 0.0, nzb_u_inner, &
	270	nzb_u_outer, wall_u )
	271	ENDIF
	272
	273	!$acc kernels present ( u, v, w, km, tend, usws, uswst ) &
	274	!$acc present ( ddzu, ddzw, fym, fyp, wall_u ) &
	275	!$acc present ( nzb_u_inner, nzb_u_outer, nzb_diff_u )
	276	!$acc loop
	277	DO i = nxlu, nxr
	278	DO j = nys, nyn
	279	!
	280	!-- Compute horizontal diffusion
	281	!$acc loop vector(32)
	282	DO k = 1, nzt
	283	IF ( k > nzb_u_outer(j,i) ) THEN
	284	!
	285	!-- Interpolate eddy diffusivities on staggered gridpoints
	286	kmyp = 0.25 * &
	287	( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) )
	288	kmym = 0.25 * &
	289	( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) )
	290
	291	tend(k,j,i) = tend(k,j,i) &
	292	& + 2.0 * ( &
	293	& km(k,j,i) * ( u(k,j,i+1) - u(k,j,i) ) &
	294	& - km(k,j,i-1) * ( u(k,j,i) - u(k,j,i-1) ) &
	295	& ) * ddx2 &
	296	& + ( kmyp * ( u(k,j+1,i) - u(k,j,i) ) * ddy &
	297	& + kmyp * ( v(k,j+1,i) - v(k,j+1,i-1) ) * ddx &
	298	& - kmym * ( u(k,j,i) - u(k,j-1,i) ) * ddy &
	299	& - kmym * ( v(k,j,i) - v(k,j,i-1) ) * ddx &
	300	& ) * ddy
	301	ENDIF
	302	ENDDO
	303
	304	!
	305	!-- Wall functions at the north and south walls, respectively
	306	!$acc loop vector(32)
	307	DO k = 1, nzt
	308	IF( k > nzb_u_inner(j,i) .AND. k <= nzb_u_outer(j,i) .AND. &
	309	wall_u(j,i) /= 0.0 ) THEN
	310
	311	kmyp = 0.25 * &
	312	( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) )
	313	kmym = 0.25 * &
	314	( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) )
	315
	316	tend(k,j,i) = tend(k,j,i) &
	317	+ 2.0 * ( &
	318	km(k,j,i) * ( u(k,j,i+1) - u(k,j,i) ) &
	319	- km(k,j,i-1) * ( u(k,j,i) - u(k,j,i-1) ) &
	320	) * ddx2 &
	321	+ ( fyp(j,i) * ( &
	322	kmyp * ( u(k,j+1,i) - u(k,j,i) ) * ddy &
	323	+ kmyp * ( v(k,j+1,i) - v(k,j+1,i-1) ) * ddx &
	324	) &
	325	- fym(j,i) * ( &
	326	kmym * ( u(k,j,i) - u(k,j-1,i) ) * ddy &
	327	+ kmym * ( v(k,j,i) - v(k,j,i-1) ) * ddx &
	328	) &
	329	+ wall_u(j,i) * usvs(k,j,i) &
	330	) * ddy
	331	ENDIF
	332	ENDDO
	333
	334	!
	335	!-- Compute vertical diffusion. In case of simulating a Prandtl layer,
	336	!-- index k starts at nzb_u_inner+2.
	337	!$acc loop vector(32)
	338	DO k = 1, nzt_diff
	339	IF ( k >= nzb_diff_u(j,i) ) THEN
	340	!
	341	!-- Interpolate eddy diffusivities on staggered gridpoints
	342	kmzp = 0.25 * &
	343	( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) )
	344	kmzm = 0.25 * &
	345	( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
	346
	347	tend(k,j,i) = tend(k,j,i) &
	348	& + ( kmzp * ( ( u(k+1,j,i) - u(k,j,i) ) * ddzu(k+1)&
	349	& + ( w(k,j,i) - w(k,j,i-1) ) * ddx &
	350	& ) &
	351	& - kmzm * ( ( u(k,j,i) - u(k-1,j,i) ) * ddzu(k)&
	352	& + ( w(k-1,j,i) - w(k-1,j,i-1) ) * ddx &
	353	& ) &
	354	& ) * ddzw(k)
	355	ENDIF
	356	ENDDO
	357
	358	ENDDO
	359	ENDDO
	360
	361	!
	362	!-- Vertical diffusion at the first grid point above the surface,
	363	!-- if the momentum flux at the bottom is given by the Prandtl law or
	364	!-- if it is prescribed by the user.
	365	!-- Difference quotient of the momentum flux is not formed over half
	366	!-- of the grid spacing (2.0*ddzw(k)) any more, since the comparison
	367	!-- with other (LES) modell showed that the values of the momentum
	368	!-- flux becomes too large in this case.
	369	!-- The term containing w(k-1,..) (see above equation) is removed here
	370	!-- because the vertical velocity is assumed to be zero at the surface.
	371	IF ( use_surface_fluxes ) THEN
	372
	373	!$acc loop
	374	DO i = nxlu, nxr
	375	!$acc loop vector(32)
	376	DO j = nys, nyn
	377
	378	k = nzb_u_inner(j,i)+1
	379	!
	380	!-- Interpolate eddy diffusivities on staggered gridpoints
	381	kmzp = 0.25 * &
	382	( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) )
	383	kmzm = 0.25 * &
	384	( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
	385
	386	tend(k,j,i) = tend(k,j,i) &
	387	& + ( kmzp * ( w(k,j,i) - w(k,j,i-1) ) * ddx &
	388	& ) * ddzw(k) &
	389	& + ( kmzp * ( u(k+1,j,i) - u(k,j,i) ) * ddzu(k+1) &
	390	& + usws(j,i) &
	391	& ) * ddzw(k)
	392	ENDDO
	393	ENDDO
	394
	395	ENDIF
	396
	397	!
	398	!-- Vertical diffusion at the first gridpoint below the top boundary,
	399	!-- if the momentum flux at the top is prescribed by the user
	400	IF ( use_top_fluxes .AND. constant_top_momentumflux ) THEN
	401
	402	k = nzt
	403
	404	!$acc loop
	405	DO i = nxlu, nxr
	406	!$acc loop vector(32)
	407	DO j = nys, nyn
	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 |