Home

Context Navigation

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

Last change on this file since 1016 was 1015, checked in by raasch, 12 years ago
Starting with changes required for GPU optimization. OpenACC statements for using NVIDIA GPUs added. Adjustment of mixing length to the Prandtl mixing length at first grid point above ground removed. mask array is set zero for ghost boundaries
Property svn:keywords set to `Id`
File size: 18.2 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |