Home

Context Navigation

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

Last change on this file since 597 was 484, checked in by raasch, 15 years ago
typo in file headers removed
Property svn:keywords set to `Id`
File size: 17.7 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |