Home

Context Navigation

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

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |