Home

Context Navigation

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

Last change on this file since 1172 was 1132, checked in by raasch, 12 years ago
r1028 documented
Property svn:keywords set to `Id`
File size: 19.2 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |