Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |