Home

Context Navigation

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

Last change on this file since 1004 was 1002, checked in by raasch, 12 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 15.4 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |