Home

Context Navigation

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

Last change on this file since 974 was 668, checked in by suehring, 14 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 17.8 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |