Home

Context Navigation

source: palm/trunk/SOURCE/diffusion_w.f90 @ 513

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

Rev	Line
[1]	1	MODULE diffusion_w_mod
	2
	3	!------------------------------------------------------------------------------!
[484]	4	! Current revisions:
[1]	5	! -----------------
[392]	6	!
[1]	7	!
	8	! Former revisions:
	9	! -----------------
[3]	10	! $Id: diffusion_w.f90 484 2010-02-05 07:36:54Z heinze $
[39]	11	!
[392]	12	! 366 2009-08-25 08:06:27Z raasch
	13	! bc_lr/bc_ns replaced by bc_lr_cyc/bc_ns_cyc
	14	!
[77]	15	! 75 2007-03-22 09:54:05Z raasch
	16	! Wall functions now include diabatic conditions, call of routine wall_fluxes,
	17	! z0 removed from argument list
	18	!
[39]	19	! 20 2007-02-26 00:12:32Z raasch
	20	! Bugfix: ddzw dimensioned 1:nzt"+1"
	21	!
[3]	22	! RCS Log replace by Id keyword, revision history cleaned up
	23	!
[1]	24	! Revision 1.12 2006/02/23 10:38:03 raasch
	25	! nzb_2d replaced by nzb_w_outer, wall functions added for all vertical walls,
	26	! +z0 in argument list
	27	! WARNING: loops containing the MAX function are still not properly vectorized!
	28	!
	29	! Revision 1.1 1997/09/12 06:24:11 raasch
	30	! Initial revision
	31	!
	32	!
	33	! Description:
	34	! ------------
	35	! Diffusion term of the w-component
	36	!------------------------------------------------------------------------------!
	37
[56]	38	USE wall_fluxes_mod
	39
[1]	40	PRIVATE
	41	PUBLIC diffusion_w
	42
	43	INTERFACE diffusion_w
	44	MODULE PROCEDURE diffusion_w
	45	MODULE PROCEDURE diffusion_w_ij
	46	END INTERFACE diffusion_w
	47
	48	CONTAINS
	49
	50
	51	!------------------------------------------------------------------------------!
	52	! Call for all grid points
	53	!------------------------------------------------------------------------------!
	54	SUBROUTINE diffusion_w( ddzu, ddzw, km, km_damp_x, km_damp_y, tend, u, v, &
[57]	55	w )
[1]	56
	57	USE control_parameters
	58	USE grid_variables
	59	USE indices
	60
	61	IMPLICIT NONE
	62
	63	INTEGER :: i, j, k
	64	REAL :: kmxm_x, kmxm_z, kmxp_x, kmxp_z, kmym_y, kmym_z, kmyp_y, &
[51]	65	kmyp_z
[20]	66	REAL :: ddzu(1:nzt+1), ddzw(1:nzt+1), km_damp_x(nxl-1:nxr+1), &
[1]	67	km_damp_y(nys-1:nyn+1)
	68	REAL :: tend(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1)
	69	REAL, DIMENSION(:,:,:), POINTER :: km, u, v, w
[56]	70	REAL, DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) :: wsus, wsvs
[1]	71
	72
[56]	73	!
	74	!-- First calculate horizontal momentum flux w'u' and/or w'v' at vertical
	75	!-- walls, if neccessary
	76	IF ( topography /= 'flat' ) THEN
[75]	77	CALL wall_fluxes( wsus, 0.0, 0.0, 0.0, 1.0, nzb_w_inner, &
[56]	78	nzb_w_outer, wall_w_x )
[75]	79	CALL wall_fluxes( wsvs, 0.0, 0.0, 1.0, 0.0, nzb_w_inner, &
[56]	80	nzb_w_outer, wall_w_y )
	81	ENDIF
	82
[1]	83	DO i = nxl, nxr
	84	DO j = nys, nyn
	85	DO k = nzb_w_outer(j,i)+1, nzt-1
	86	!
	87	!-- Interpolate eddy diffusivities on staggered gridpoints
	88	kmxp_x = 0.25 * &
	89	( km(k,j,i)+km(k,j,i+1)+km(k+1,j,i)+km(k+1,j,i+1) )
	90	kmxm_x = 0.25 * &
	91	( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) )
	92	kmxp_z = kmxp_x
	93	kmxm_z = kmxm_x
	94	kmyp_y = 0.25 * &
	95	( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) )
	96	kmym_y = 0.25 * &
	97	( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
	98	kmyp_z = kmyp_y
	99	kmym_z = kmym_y
	100	!
	101	!-- Increase diffusion at the outflow boundary in case of
	102	!-- non-cyclic lateral boundaries. Damping is only needed for
	103	!-- velocity components parallel to the outflow boundary in
	104	!-- the direction normal to the outflow boundary.
[366]	105	IF ( .NOT. bc_lr_cyc ) THEN
[1]	106	kmxp_x = MAX( kmxp_x, km_damp_x(i) )
	107	kmxm_x = MAX( kmxm_x, km_damp_x(i) )
	108	ENDIF
[366]	109	IF ( .NOT. bc_ns_cyc ) THEN
[1]	110	kmyp_y = MAX( kmyp_y, km_damp_y(j) )
	111	kmym_y = MAX( kmym_y, km_damp_y(j) )
	112	ENDIF
	113
	114	tend(k,j,i) = tend(k,j,i) &
	115	& + ( kmxp_x * ( w(k,j,i+1) - w(k,j,i) ) * ddx &
	116	& + kmxp_z * ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1) &
	117	& - kmxm_x * ( w(k,j,i) - w(k,j,i-1) ) * ddx &
	118	& - kmxm_z * ( u(k+1,j,i) - u(k,j,i) ) * ddzu(k+1) &
	119	& ) * ddx &
	120	& + ( kmyp_y * ( w(k,j+1,i) - w(k,j,i) ) * ddy &
	121	& + kmyp_z * ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1) &
	122	& - kmym_y * ( w(k,j,i) - w(k,j-1,i) ) * ddy &
	123	& - kmym_z * ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1) &
	124	& ) * ddy &
	125	& + 2.0 * ( &
	126	& km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1) &
	127	& - km(k,j,i) * ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) &
	128	& ) * ddzu(k+1)
	129	ENDDO
	130
	131	!
	132	!-- Wall functions at all vertical walls, where necessary
	133	IF ( wall_w_x(j,i) /= 0.0 .OR. wall_w_y(j,i) /= 0.0 ) THEN
[51]	134
[1]	135	DO k = nzb_w_inner(j,i)+1, nzb_w_outer(j,i)
	136	!
	137	!-- Interpolate eddy diffusivities on staggered gridpoints
	138	kmxp_x = 0.25 * &
	139	( km(k,j,i)+km(k,j,i+1)+km(k+1,j,i)+km(k+1,j,i+1) )
	140	kmxm_x = 0.25 * &
	141	( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) )
	142	kmxp_z = kmxp_x
	143	kmxm_z = kmxm_x
	144	kmyp_y = 0.25 * &
	145	( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) )
	146	kmym_y = 0.25 * &
	147	( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
	148	kmyp_z = kmyp_y
	149	kmym_z = kmym_y
	150	!
	151	!-- Increase diffusion at the outflow boundary in case of
	152	!-- non-cyclic lateral boundaries. Damping is only needed for
	153	!-- velocity components parallel to the outflow boundary in
	154	!-- the direction normal to the outflow boundary.
[366]	155	IF ( .NOT. bc_lr_cyc ) THEN
[1]	156	kmxp_x = MAX( kmxp_x, km_damp_x(i) )
	157	kmxm_x = MAX( kmxm_x, km_damp_x(i) )
	158	ENDIF
[366]	159	IF ( .NOT. bc_ns_cyc ) THEN
[1]	160	kmyp_y = MAX( kmyp_y, km_damp_y(j) )
	161	kmym_y = MAX( kmym_y, km_damp_y(j) )
	162	ENDIF
	163
	164	tend(k,j,i) = tend(k,j,i) &
	165	+ ( fwxp(j,i) * ( &
	166	kmxp_x * ( w(k,j,i+1) - w(k,j,i) ) * ddx &
	167	+ kmxp_z * ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1) &
	168	) &
	169	- fwxm(j,i) * ( &
	170	kmxm_x * ( w(k,j,i) - w(k,j,i-1) ) * ddx &
	171	+ kmxm_z * ( u(k+1,j,i) - u(k,j,i) ) * ddzu(k+1) &
	172	) &
[56]	173	+ wall_w_x(j,i) * wsus(k,j,i) &
[1]	174	) * ddx &
	175	+ ( fwyp(j,i) * ( &
	176	kmyp_y * ( w(k,j+1,i) - w(k,j,i) ) * ddy &
	177	+ kmyp_z * ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1) &
	178	) &
	179	- fwym(j,i) * ( &
	180	kmym_y * ( w(k,j,i) - w(k,j-1,i) ) * ddy &
	181	+ kmym_z * ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1) &
	182	) &
[56]	183	+ wall_w_y(j,i) * wsvs(k,j,i) &
[1]	184	) * ddy &
	185	+ 2.0 * ( &
	186	km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1) &
	187	- km(k,j,i) * ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) &
	188	) * ddzu(k+1)
	189	ENDDO
	190	ENDIF
	191
	192	ENDDO
	193	ENDDO
	194
	195	END SUBROUTINE diffusion_w
	196
	197
	198	!------------------------------------------------------------------------------!
	199	! Call for grid point i,j
	200	!------------------------------------------------------------------------------!
	201	SUBROUTINE diffusion_w_ij( i, j, ddzu, ddzw, km, km_damp_x, km_damp_y, &
[57]	202	tend, u, v, w )
[1]	203
	204	USE control_parameters
	205	USE grid_variables
	206	USE indices
	207
	208	IMPLICIT NONE
	209
	210	INTEGER :: i, j, k
	211	REAL :: kmxm_x, kmxm_z, kmxp_x, kmxp_z, kmym_y, kmym_z, kmyp_y, &
[51]	212	kmyp_z
[20]	213	REAL :: ddzu(1:nzt+1), ddzw(1:nzt+1), km_damp_x(nxl-1:nxr+1), &
[1]	214	km_damp_y(nys-1:nyn+1)
	215	REAL :: tend(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1)
[51]	216	REAL, DIMENSION(nzb:nzt+1) :: wsus, wsvs
[1]	217	REAL, DIMENSION(:,:,:), POINTER :: km, u, v, w
	218
	219
	220	DO k = nzb_w_outer(j,i)+1, nzt-1
	221	!
	222	!-- Interpolate eddy diffusivities on staggered gridpoints
	223	kmxp_x = 0.25 * ( km(k,j,i)+km(k,j,i+1)+km(k+1,j,i)+km(k+1,j,i+1) )
	224	kmxm_x = 0.25 * ( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) )
	225	kmxp_z = kmxp_x
	226	kmxm_z = kmxm_x
	227	kmyp_y = 0.25 * ( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) )
	228	kmym_y = 0.25 * ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
	229	kmyp_z = kmyp_y
	230	kmym_z = kmym_y
	231	!
	232	!-- Increase diffusion at the outflow boundary in case of non-cyclic
	233	!-- lateral boundaries. Damping is only needed for velocity components
	234	!-- parallel to the outflow boundary in the direction normal to the
	235	!-- outflow boundary.
[366]	236	IF ( .NOT. bc_lr_cyc ) THEN
[1]	237	kmxp_x = MAX( kmxp_x, km_damp_x(i) )
	238	kmxm_x = MAX( kmxm_x, km_damp_x(i) )
	239	ENDIF
[366]	240	IF ( .NOT. bc_ns_cyc ) THEN
[1]	241	kmyp_y = MAX( kmyp_y, km_damp_y(j) )
	242	kmym_y = MAX( kmym_y, km_damp_y(j) )
	243	ENDIF
	244
	245	tend(k,j,i) = tend(k,j,i) &
	246	& + ( kmxp_x * ( w(k,j,i+1) - w(k,j,i) ) * ddx &
	247	& + kmxp_z * ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1) &
	248	& - kmxm_x * ( w(k,j,i) - w(k,j,i-1) ) * ddx &
	249	& - kmxm_z * ( u(k+1,j,i) - u(k,j,i) ) * ddzu(k+1) &
	250	& ) * ddx &
	251	& + ( kmyp_y * ( w(k,j+1,i) - w(k,j,i) ) * ddy &
	252	& + kmyp_z * ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1) &
	253	& - kmym_y * ( w(k,j,i) - w(k,j-1,i) ) * ddy &
	254	& - kmym_z * ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1) &
	255	& ) * ddy &
	256	& + 2.0 * ( &
	257	& km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1) &
	258	& - km(k,j,i) * ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) &
	259	& ) * ddzu(k+1)
	260	ENDDO
	261
	262	!
	263	!-- Wall functions at all vertical walls, where necessary
	264	IF ( wall_w_x(j,i) /= 0.0 .OR. wall_w_y(j,i) /= 0.0 ) THEN
[51]	265
	266	!
	267	!-- Calculate the horizontal momentum fluxes w'u' and/or w'v'
	268	IF ( wall_w_x(j,i) /= 0.0 ) THEN
	269	CALL wall_fluxes( i, j, nzb_w_inner(j,i)+1, nzb_w_outer(j,i), &
	270	wsus, 0.0, 0.0, 0.0, 1.0 )
	271	ELSE
	272	wsus = 0.0
	273	ENDIF
	274
	275	IF ( wall_w_y(j,i) /= 0.0 ) THEN
	276	CALL wall_fluxes( i, j, nzb_w_inner(j,i)+1, nzb_w_outer(j,i), &
	277	wsvs, 0.0, 0.0, 1.0, 0.0 )
	278	ELSE
	279	wsvs = 0.0
	280	ENDIF
	281
[1]	282	DO k = nzb_w_inner(j,i)+1, nzb_w_outer(j,i)
	283	!
	284	!-- Interpolate eddy diffusivities on staggered gridpoints
	285	kmxp_x = 0.25 * ( km(k,j,i)+km(k,j,i+1)+km(k+1,j,i)+km(k+1,j,i+1) )
	286	kmxm_x = 0.25 * ( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) )
	287	kmxp_z = kmxp_x
	288	kmxm_z = kmxm_x
	289	kmyp_y = 0.25 * ( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) )
	290	kmym_y = 0.25 * ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
	291	kmyp_z = kmyp_y
	292	kmym_z = kmym_y
	293	!
	294	!-- Increase diffusion at the outflow boundary in case of
	295	!-- non-cyclic lateral boundaries. Damping is only needed for
	296	!-- velocity components parallel to the outflow boundary in
	297	!-- the direction normal to the outflow boundary.
[366]	298	IF ( .NOT. bc_lr_cyc ) THEN
[1]	299	kmxp_x = MAX( kmxp_x, km_damp_x(i) )
	300	kmxm_x = MAX( kmxm_x, km_damp_x(i) )
	301	ENDIF
[366]	302	IF ( .NOT. bc_ns_cyc ) THEN
[1]	303	kmyp_y = MAX( kmyp_y, km_damp_y(j) )
	304	kmym_y = MAX( kmym_y, km_damp_y(j) )
	305	ENDIF
	306
	307	tend(k,j,i) = tend(k,j,i) &
	308	+ ( fwxp(j,i) * ( &
	309	kmxp_x * ( w(k,j,i+1) - w(k,j,i) ) * ddx &
	310	+ kmxp_z * ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1) &
	311	) &
	312	- fwxm(j,i) * ( &
	313	kmxm_x * ( w(k,j,i) - w(k,j,i-1) ) * ddx &
	314	+ kmxm_z * ( u(k+1,j,i) - u(k,j,i) ) * ddzu(k+1) &
	315	) &
[51]	316	+ wall_w_x(j,i) * wsus(k) &
[1]	317	) * ddx &
	318	+ ( fwyp(j,i) * ( &
	319	kmyp_y * ( w(k,j+1,i) - w(k,j,i) ) * ddy &
	320	+ kmyp_z * ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1) &
	321	) &
	322	- fwym(j,i) * ( &
	323	kmym_y * ( w(k,j,i) - w(k,j-1,i) ) * ddy &
	324	+ kmym_z * ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1) &
	325	) &
[51]	326	+ wall_w_y(j,i) * wsvs(k) &
[1]	327	) * ddy &
	328	+ 2.0 * ( &
	329	km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1) &
	330	- km(k,j,i) * ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) &
	331	) * ddzu(k+1)
	332	ENDDO
	333	ENDIF
	334
	335	END SUBROUTINE diffusion_w_ij
	336
	337	END MODULE diffusion_w_mod

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |