Home

Context Navigation

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

Last change on this file since 377 was 366, checked in by raasch, 15 years ago
speed optomizations +bugfix in init_ocean
Property svn:keywords set to `Id`
File size: 15.3 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |