Home

Context Navigation

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

Last change on this file since 980 was 979, checked in by fricke, 12 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 15.8 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |