Home

Context Navigation

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

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

Rev	Line
[1]	1	MODULE diffusion_w_mod
	2
	3	!------------------------------------------------------------------------------!
[484]	4	! Current revisions:
[1]	5	! -----------------
	6	!
[1017]	7	!
[1]	8	! Former revisions:
	9	! -----------------
[3]	10	! $Id: diffusion_w.f90 1017 2012-09-27 11:28:50Z letzel $
[39]	11	!
[1017]	12	! 1015 2012-09-27 09:23:24Z raasch
	13	! accelerator version (*_acc) added
	14	!
[1002]	15	! 1001 2012-09-13 14:08:46Z raasch
	16	! arrays comunicated by module instead of parameter list
	17	!
[979]	18	! 978 2012-08-09 08:28:32Z fricke
	19	! outflow damping layer removed
	20	! kmxm_x/_z and kmxp_x/_z change to kmxm and kmxp
	21	! kmym_y/_z and kmyp_y/_z change to kmym and kmyp
	22	!
[668]	23	! 667 2010-12-23 12:06:00Z suehring/gryschka
	24	! nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng
	25	!
[392]	26	! 366 2009-08-25 08:06:27Z raasch
	27	! bc_lr/bc_ns replaced by bc_lr_cyc/bc_ns_cyc
	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
	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.12 2006/02/23 10:38:03 raasch
	39	! nzb_2d replaced by nzb_w_outer, wall functions added for all vertical walls,
	40	! +z0 in argument list
	41	! WARNING: loops containing the MAX function are still not properly vectorized!
	42	!
	43	! Revision 1.1 1997/09/12 06:24:11 raasch
	44	! Initial revision
	45	!
	46	!
	47	! Description:
	48	! ------------
	49	! Diffusion term of the w-component
	50	!------------------------------------------------------------------------------!
	51
[56]	52	USE wall_fluxes_mod
	53
[1]	54	PRIVATE
[1015]	55	PUBLIC diffusion_w, diffusion_w_acc
[1]	56
	57	INTERFACE diffusion_w
	58	MODULE PROCEDURE diffusion_w
	59	MODULE PROCEDURE diffusion_w_ij
	60	END INTERFACE diffusion_w
	61
[1015]	62	INTERFACE diffusion_w_acc
	63	MODULE PROCEDURE diffusion_w_acc
	64	END INTERFACE diffusion_w_acc
	65
[1]	66	CONTAINS
	67
	68
	69	!------------------------------------------------------------------------------!
	70	! Call for all grid points
	71	!------------------------------------------------------------------------------!
[1001]	72	SUBROUTINE diffusion_w
[1]	73
[1001]	74	USE arrays_3d
[1]	75	USE control_parameters
	76	USE grid_variables
	77	USE indices
	78
	79	IMPLICIT NONE
	80
	81	INTEGER :: i, j, k
[978]	82	REAL :: kmxm, kmxp, kmym, kmyp
[1001]	83
[56]	84	REAL, DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) :: wsus, wsvs
[1]	85
	86
[56]	87	!
	88	!-- First calculate horizontal momentum flux w'u' and/or w'v' at vertical
	89	!-- walls, if neccessary
	90	IF ( topography /= 'flat' ) THEN
[75]	91	CALL wall_fluxes( wsus, 0.0, 0.0, 0.0, 1.0, nzb_w_inner, &
[56]	92	nzb_w_outer, wall_w_x )
[75]	93	CALL wall_fluxes( wsvs, 0.0, 0.0, 1.0, 0.0, nzb_w_inner, &
[56]	94	nzb_w_outer, wall_w_y )
	95	ENDIF
	96
[1]	97	DO i = nxl, nxr
	98	DO j = nys, nyn
	99	DO k = nzb_w_outer(j,i)+1, nzt-1
	100	!
	101	!-- Interpolate eddy diffusivities on staggered gridpoints
[978]	102	kmxp = 0.25 * &
	103	( km(k,j,i)+km(k,j,i+1)+km(k+1,j,i)+km(k+1,j,i+1) )
	104	kmxm = 0.25 * &
	105	( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) )
	106	kmyp = 0.25 * &
	107	( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) )
	108	kmym = 0.25 * &
	109	( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
[1]	110
	111	tend(k,j,i) = tend(k,j,i) &
[978]	112	& + ( kmxp * ( w(k,j,i+1) - w(k,j,i) ) * ddx &
	113	& + kmxp * ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1) &
	114	& - kmxm * ( w(k,j,i) - w(k,j,i-1) ) * ddx &
	115	& - kmxm * ( u(k+1,j,i) - u(k,j,i) ) * ddzu(k+1) &
[1]	116	& ) * ddx &
[978]	117	& + ( kmyp * ( w(k,j+1,i) - w(k,j,i) ) * ddy &
	118	& + kmyp * ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1) &
	119	& - kmym * ( w(k,j,i) - w(k,j-1,i) ) * ddy &
	120	& - kmym * ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1) &
[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
[978]	135	kmxp = 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 = 0.25 * &
	138	( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) )
	139	kmyp = 0.25 * &
	140	( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) )
	141	kmym = 0.25 * &
	142	( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
[1]	143
	144	tend(k,j,i) = tend(k,j,i) &
	145	+ ( fwxp(j,i) * ( &
[978]	146	kmxp * ( w(k,j,i+1) - w(k,j,i) ) * ddx &
	147	+ kmxp * ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1) &
[1]	148	) &
	149	- fwxm(j,i) * ( &
[978]	150	kmxm * ( w(k,j,i) - w(k,j,i-1) ) * ddx &
	151	+ kmxm * ( u(k+1,j,i) - u(k,j,i) ) * ddzu(k+1) &
[1]	152	) &
[56]	153	+ wall_w_x(j,i) * wsus(k,j,i) &
[1]	154	) * ddx &
	155	+ ( fwyp(j,i) * ( &
[978]	156	kmyp * ( w(k,j+1,i) - w(k,j,i) ) * ddy &
	157	+ kmyp * ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1) &
[1]	158	) &
	159	- fwym(j,i) * ( &
[978]	160	kmym * ( w(k,j,i) - w(k,j-1,i) ) * ddy &
	161	+ kmym * ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1) &
[1]	162	) &
[56]	163	+ wall_w_y(j,i) * wsvs(k,j,i) &
[1]	164	) * ddy &
	165	+ 2.0 * ( &
	166	km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1) &
	167	- km(k,j,i) * ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) &
	168	) * ddzu(k+1)
	169	ENDDO
	170	ENDIF
	171
	172	ENDDO
	173	ENDDO
	174
	175	END SUBROUTINE diffusion_w
	176
	177
	178	!------------------------------------------------------------------------------!
[1015]	179	! Call for all grid points - accelerator version
	180	!------------------------------------------------------------------------------!
	181	SUBROUTINE diffusion_w_acc
	182
	183	USE arrays_3d
	184	USE control_parameters
	185	USE grid_variables
	186	USE indices
	187
	188	IMPLICIT NONE
	189
	190	INTEGER :: i, j, k
	191	REAL :: kmxm, kmxp, kmym, kmyp
	192
	193	!$acc declare create ( wsus, wsvs )
	194	REAL, DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) :: wsus, wsvs
	195
	196
	197	!
	198	!-- First calculate horizontal momentum flux w'u' and/or w'v' at vertical
	199	!-- walls, if neccessary
	200	IF ( topography /= 'flat' ) THEN
	201	CALL wall_fluxes_acc( wsus, 0.0, 0.0, 0.0, 1.0, nzb_w_inner, &
	202	nzb_w_outer, wall_w_x )
	203	CALL wall_fluxes_acc( wsvs, 0.0, 0.0, 1.0, 0.0, nzb_w_inner, &
	204	nzb_w_outer, wall_w_y )
	205	ENDIF
	206
	207	!$acc kernels present ( u, v, w, km, tend, vsws, vswst ) &
	208	!$acc present ( ddzu, ddzw, fwxm, fwxp, fwym, fwyp, wall_w_x, wall_w_y ) &
	209	!$acc present ( nzb_w_inner, nzb_w_outer )
	210	!$acc loop
	211	DO i = nxl, nxr
	212	DO j = nys, nyn
	213	!$acc loop vector( 32 )
	214	DO k = 1, nzt
	215	IF ( k > nzb_w_outer(j,i) ) THEN
	216	!
	217	!-- Interpolate eddy diffusivities on staggered gridpoints
	218	kmxp = 0.25 * &
	219	( km(k,j,i)+km(k,j,i+1)+km(k+1,j,i)+km(k+1,j,i+1) )
	220	kmxm = 0.25 * &
	221	( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) )
	222	kmyp = 0.25 * &
	223	( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) )
	224	kmym = 0.25 * &
	225	( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
	226
	227	tend(k,j,i) = tend(k,j,i) &
	228	& + ( kmxp * ( w(k,j,i+1) - w(k,j,i) ) * ddx &
	229	& + kmxp * ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1) &
	230	& - kmxm * ( w(k,j,i) - w(k,j,i-1) ) * ddx &
	231	& - kmxm * ( u(k+1,j,i) - u(k,j,i) ) * ddzu(k+1) &
	232	& ) * ddx &
	233	& + ( kmyp * ( w(k,j+1,i) - w(k,j,i) ) * ddy &
	234	& + kmyp * ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1) &
	235	& - kmym * ( w(k,j,i) - w(k,j-1,i) ) * ddy &
	236	& - kmym * ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1) &
	237	& ) * ddy &
	238	& + 2.0 * ( &
	239	& km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1) &
	240	& - km(k,j,i) * ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) &
	241	& ) * ddzu(k+1)
	242	ENDIF
	243	ENDDO
	244
	245	!
	246	!-- Wall functions at all vertical walls, where necessary
	247	!$acc loop vector( 32 )
	248	DO k = 1,nzt
	249
	250	IF ( k > nzb_w_inner(j,i) .AND. k <= nzb_w_outer(j,i) .AND. &
	251	wall_w_x(j,i) /= 0.0 .AND. wall_w_y(j,i) /= 0.0 ) THEN
	252	!
	253	!-- Interpolate eddy diffusivities on staggered gridpoints
	254	kmxp = 0.25 * &
	255	( km(k,j,i)+km(k,j,i+1)+km(k+1,j,i)+km(k+1,j,i+1) )
	256	kmxm = 0.25 * &
	257	( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) )
	258	kmyp = 0.25 * &
	259	( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) )
	260	kmym = 0.25 * &
	261	( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
	262
	263	tend(k,j,i) = tend(k,j,i) &
	264	+ ( fwxp(j,i) * ( &
	265	kmxp * ( w(k,j,i+1) - w(k,j,i) ) * ddx &
	266	+ kmxp * ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1) &
	267	) &
	268	- fwxm(j,i) * ( &
	269	kmxm * ( w(k,j,i) - w(k,j,i-1) ) * ddx &
	270	+ kmxm * ( u(k+1,j,i) - u(k,j,i) ) * ddzu(k+1) &
	271	) &
	272	+ wall_w_x(j,i) * wsus(k,j,i) &
	273	) * ddx &
	274	+ ( fwyp(j,i) * ( &
	275	kmyp * ( w(k,j+1,i) - w(k,j,i) ) * ddy &
	276	+ kmyp * ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1) &
	277	) &
	278	- fwym(j,i) * ( &
	279	kmym * ( w(k,j,i) - w(k,j-1,i) ) * ddy &
	280	+ kmym * ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1) &
	281	) &
	282	+ wall_w_y(j,i) * wsvs(k,j,i) &
	283	) * ddy &
	284	+ 2.0 * ( &
	285	km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1) &
	286	- km(k,j,i) * ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) &
	287	) * ddzu(k+1)
	288	ENDIF
	289	ENDDO
	290
	291	ENDDO
	292	ENDDO
	293	!$acc end kernels
	294
	295	END SUBROUTINE diffusion_w_acc
	296
	297
	298	!------------------------------------------------------------------------------!
[1]	299	! Call for grid point i,j
	300	!------------------------------------------------------------------------------!
[1001]	301	SUBROUTINE diffusion_w_ij( i, j )
[1]	302
[1001]	303	USE arrays_3d
[1]	304	USE control_parameters
	305	USE grid_variables
	306	USE indices
	307
	308	IMPLICIT NONE
	309
	310	INTEGER :: i, j, k
[978]	311	REAL :: kmxm, kmxp, kmym, kmyp
[1]	312
[1001]	313	REAL, DIMENSION(nzb:nzt+1) :: wsus, wsvs
[1]	314
[1001]	315
[1]	316	DO k = nzb_w_outer(j,i)+1, nzt-1
	317	!
	318	!-- Interpolate eddy diffusivities on staggered gridpoints
[978]	319	kmxp = 0.25 * ( km(k,j,i)+km(k,j,i+1)+km(k+1,j,i)+km(k+1,j,i+1) )
	320	kmxm = 0.25 * ( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) )
	321	kmyp = 0.25 * ( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) )
	322	kmym = 0.25 * ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
[1]	323
	324	tend(k,j,i) = tend(k,j,i) &
[978]	325	& + ( kmxp * ( w(k,j,i+1) - w(k,j,i) ) * ddx &
	326	& + kmxp * ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1) &
	327	& - kmxm * ( w(k,j,i) - w(k,j,i-1) ) * ddx &
	328	& - kmxm * ( u(k+1,j,i) - u(k,j,i) ) * ddzu(k+1) &
[1]	329	& ) * ddx &
[978]	330	& + ( kmyp * ( w(k,j+1,i) - w(k,j,i) ) * ddy &
	331	& + kmyp * ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1) &
	332	& - kmym * ( w(k,j,i) - w(k,j-1,i) ) * ddy &
	333	& - kmym * ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1) &
[1]	334	& ) * ddy &
	335	& + 2.0 * ( &
	336	& km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1) &
	337	& - km(k,j,i) * ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) &
	338	& ) * ddzu(k+1)
	339	ENDDO
	340
	341	!
	342	!-- Wall functions at all vertical walls, where necessary
	343	IF ( wall_w_x(j,i) /= 0.0 .OR. wall_w_y(j,i) /= 0.0 ) THEN
[51]	344
	345	!
	346	!-- Calculate the horizontal momentum fluxes w'u' and/or w'v'
	347	IF ( wall_w_x(j,i) /= 0.0 ) THEN
	348	CALL wall_fluxes( i, j, nzb_w_inner(j,i)+1, nzb_w_outer(j,i), &
	349	wsus, 0.0, 0.0, 0.0, 1.0 )
	350	ELSE
	351	wsus = 0.0
	352	ENDIF
	353
	354	IF ( wall_w_y(j,i) /= 0.0 ) THEN
	355	CALL wall_fluxes( i, j, nzb_w_inner(j,i)+1, nzb_w_outer(j,i), &
	356	wsvs, 0.0, 0.0, 1.0, 0.0 )
	357	ELSE
	358	wsvs = 0.0
	359	ENDIF
	360
[1]	361	DO k = nzb_w_inner(j,i)+1, nzb_w_outer(j,i)
	362	!
	363	!-- Interpolate eddy diffusivities on staggered gridpoints
[978]	364	kmxp = 0.25 * ( km(k,j,i)+km(k,j,i+1)+km(k+1,j,i)+km(k+1,j,i+1) )
	365	kmxm = 0.25 * ( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) )
	366	kmyp = 0.25 * ( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) )
	367	kmym = 0.25 * ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
[1]	368
	369	tend(k,j,i) = tend(k,j,i) &
	370	+ ( fwxp(j,i) * ( &
[978]	371	kmxp * ( w(k,j,i+1) - w(k,j,i) ) * ddx &
	372	+ kmxp * ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1) &
[1]	373	) &
	374	- fwxm(j,i) * ( &
[978]	375	kmxm * ( w(k,j,i) - w(k,j,i-1) ) * ddx &
	376	+ kmxm * ( u(k+1,j,i) - u(k,j,i) ) * ddzu(k+1) &
[1]	377	) &
[51]	378	+ wall_w_x(j,i) * wsus(k) &
[1]	379	) * ddx &
	380	+ ( fwyp(j,i) * ( &
[978]	381	kmyp * ( w(k,j+1,i) - w(k,j,i) ) * ddy &
	382	+ kmyp * ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1) &
[1]	383	) &
	384	- fwym(j,i) * ( &
[978]	385	kmym * ( w(k,j,i) - w(k,j-1,i) ) * ddy &
	386	+ kmym * ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1) &
[1]	387	) &
[51]	388	+ wall_w_y(j,i) * wsvs(k) &
[1]	389	) * ddy &
	390	+ 2.0 * ( &
	391	km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1) &
	392	- km(k,j,i) * ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) &
	393	) * ddzu(k+1)
	394	ENDDO
	395	ENDIF
	396
	397	END SUBROUTINE diffusion_w_ij
	398
	399	END MODULE diffusion_w_mod

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |