Home

Context Navigation

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

Last change on this file since 1313 was 1310, checked in by raasch, 11 years ago
update of GPL copyright
Property svn:keywords set to `Id`
File size: 19.3 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	!
[1310]	17	! Copyright 1997-2014 Leibniz Universitaet Hannover
[1036]	18	!--------------------------------------------------------------------------------!
	19	!
[484]	20	! Current revisions:
[1]	21	! -----------------
	22	!
[1258]	23	!
[1]	24	! Former revisions:
	25	! -----------------
[3]	26	! $Id: diffusion_w.f90 1310 2014-03-14 08:01:56Z fricke $
[39]	27	!
[1258]	28	! 1257 2013-11-08 15:18:40Z raasch
	29	! openacc loop and loop vector clauses removed, declare create moved after
	30	! the FORTRAN declaration statement
	31	!
[1132]	32	! 1128 2013-04-12 06:19:32Z raasch
	33	! loop index bounds in accelerator version replaced by i_left, i_right, j_south,
	34	! j_north
	35	!
[1037]	36	! 1036 2012-10-22 13:43:42Z raasch
	37	! code put under GPL (PALM 3.9)
	38	!
[1017]	39	! 1015 2012-09-27 09:23:24Z raasch
	40	! accelerator version (*_acc) added
	41	!
[1002]	42	! 1001 2012-09-13 14:08:46Z raasch
	43	! arrays comunicated by module instead of parameter list
	44	!
[979]	45	! 978 2012-08-09 08:28:32Z fricke
	46	! outflow damping layer removed
	47	! kmxm_x/_z and kmxp_x/_z change to kmxm and kmxp
	48	! kmym_y/_z and kmyp_y/_z change to kmym and kmyp
	49	!
[668]	50	! 667 2010-12-23 12:06:00Z suehring/gryschka
	51	! nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng
	52	!
[392]	53	! 366 2009-08-25 08:06:27Z raasch
	54	! bc_lr/bc_ns replaced by bc_lr_cyc/bc_ns_cyc
	55	!
[77]	56	! 75 2007-03-22 09:54:05Z raasch
	57	! Wall functions now include diabatic conditions, call of routine wall_fluxes,
	58	! z0 removed from argument list
	59	!
[39]	60	! 20 2007-02-26 00:12:32Z raasch
	61	! Bugfix: ddzw dimensioned 1:nzt"+1"
	62	!
[3]	63	! RCS Log replace by Id keyword, revision history cleaned up
	64	!
[1]	65	! Revision 1.12 2006/02/23 10:38:03 raasch
	66	! nzb_2d replaced by nzb_w_outer, wall functions added for all vertical walls,
	67	! +z0 in argument list
	68	! WARNING: loops containing the MAX function are still not properly vectorized!
	69	!
	70	! Revision 1.1 1997/09/12 06:24:11 raasch
	71	! Initial revision
	72	!
	73	!
	74	! Description:
	75	! ------------
	76	! Diffusion term of the w-component
	77	!------------------------------------------------------------------------------!
	78
[56]	79	USE wall_fluxes_mod
	80
[1]	81	PRIVATE
[1015]	82	PUBLIC diffusion_w, diffusion_w_acc
[1]	83
	84	INTERFACE diffusion_w
	85	MODULE PROCEDURE diffusion_w
	86	MODULE PROCEDURE diffusion_w_ij
	87	END INTERFACE diffusion_w
	88
[1015]	89	INTERFACE diffusion_w_acc
	90	MODULE PROCEDURE diffusion_w_acc
	91	END INTERFACE diffusion_w_acc
	92
[1]	93	CONTAINS
	94
	95
	96	!------------------------------------------------------------------------------!
	97	! Call for all grid points
	98	!------------------------------------------------------------------------------!
[1001]	99	SUBROUTINE diffusion_w
[1]	100
[1001]	101	USE arrays_3d
[1]	102	USE control_parameters
	103	USE grid_variables
	104	USE indices
	105
	106	IMPLICIT NONE
	107
	108	INTEGER :: i, j, k
[978]	109	REAL :: kmxm, kmxp, kmym, kmyp
[1001]	110
[56]	111	REAL, DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) :: wsus, wsvs
[1]	112
	113
[56]	114	!
	115	!-- First calculate horizontal momentum flux w'u' and/or w'v' at vertical
	116	!-- walls, if neccessary
	117	IF ( topography /= 'flat' ) THEN
[75]	118	CALL wall_fluxes( wsus, 0.0, 0.0, 0.0, 1.0, nzb_w_inner, &
[56]	119	nzb_w_outer, wall_w_x )
[75]	120	CALL wall_fluxes( wsvs, 0.0, 0.0, 1.0, 0.0, nzb_w_inner, &
[56]	121	nzb_w_outer, wall_w_y )
	122	ENDIF
	123
[1]	124	DO i = nxl, nxr
	125	DO j = nys, nyn
	126	DO k = nzb_w_outer(j,i)+1, nzt-1
	127	!
	128	!-- Interpolate eddy diffusivities on staggered gridpoints
[978]	129	kmxp = 0.25 * &
	130	( km(k,j,i)+km(k,j,i+1)+km(k+1,j,i)+km(k+1,j,i+1) )
	131	kmxm = 0.25 * &
	132	( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) )
	133	kmyp = 0.25 * &
	134	( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) )
	135	kmym = 0.25 * &
	136	( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
[1]	137
	138	tend(k,j,i) = tend(k,j,i) &
[978]	139	& + ( kmxp * ( w(k,j,i+1) - w(k,j,i) ) * ddx &
	140	& + kmxp * ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1) &
	141	& - kmxm * ( w(k,j,i) - w(k,j,i-1) ) * ddx &
	142	& - kmxm * ( u(k+1,j,i) - u(k,j,i) ) * ddzu(k+1) &
[1]	143	& ) * ddx &
[978]	144	& + ( kmyp * ( w(k,j+1,i) - w(k,j,i) ) * ddy &
	145	& + kmyp * ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1) &
	146	& - kmym * ( w(k,j,i) - w(k,j-1,i) ) * ddy &
	147	& - kmym * ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1) &
[1]	148	& ) * ddy &
	149	& + 2.0 * ( &
	150	& km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1) &
	151	& - km(k,j,i) * ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) &
	152	& ) * ddzu(k+1)
	153	ENDDO
	154
	155	!
	156	!-- Wall functions at all vertical walls, where necessary
	157	IF ( wall_w_x(j,i) /= 0.0 .OR. wall_w_y(j,i) /= 0.0 ) THEN
[51]	158
[1]	159	DO k = nzb_w_inner(j,i)+1, nzb_w_outer(j,i)
	160	!
	161	!-- Interpolate eddy diffusivities on staggered gridpoints
[978]	162	kmxp = 0.25 * &
	163	( km(k,j,i)+km(k,j,i+1)+km(k+1,j,i)+km(k+1,j,i+1) )
	164	kmxm = 0.25 * &
	165	( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) )
	166	kmyp = 0.25 * &
	167	( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) )
	168	kmym = 0.25 * &
	169	( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
[1]	170
	171	tend(k,j,i) = tend(k,j,i) &
	172	+ ( fwxp(j,i) * ( &
[978]	173	kmxp * ( w(k,j,i+1) - w(k,j,i) ) * ddx &
	174	+ kmxp * ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1) &
[1]	175	) &
	176	- fwxm(j,i) * ( &
[978]	177	kmxm * ( w(k,j,i) - w(k,j,i-1) ) * ddx &
	178	+ kmxm * ( u(k+1,j,i) - u(k,j,i) ) * ddzu(k+1) &
[1]	179	) &
[56]	180	+ wall_w_x(j,i) * wsus(k,j,i) &
[1]	181	) * ddx &
	182	+ ( fwyp(j,i) * ( &
[978]	183	kmyp * ( w(k,j+1,i) - w(k,j,i) ) * ddy &
	184	+ kmyp * ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1) &
[1]	185	) &
	186	- fwym(j,i) * ( &
[978]	187	kmym * ( w(k,j,i) - w(k,j-1,i) ) * ddy &
	188	+ kmym * ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1) &
[1]	189	) &
[56]	190	+ wall_w_y(j,i) * wsvs(k,j,i) &
[1]	191	) * ddy &
	192	+ 2.0 * ( &
	193	km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1) &
	194	- km(k,j,i) * ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) &
	195	) * ddzu(k+1)
	196	ENDDO
	197	ENDIF
	198
	199	ENDDO
	200	ENDDO
	201
	202	END SUBROUTINE diffusion_w
	203
	204
	205	!------------------------------------------------------------------------------!
[1015]	206	! Call for all grid points - accelerator version
	207	!------------------------------------------------------------------------------!
	208	SUBROUTINE diffusion_w_acc
	209
	210	USE arrays_3d
	211	USE control_parameters
	212	USE grid_variables
	213	USE indices
	214
	215	IMPLICIT NONE
	216
	217	INTEGER :: i, j, k
	218	REAL :: kmxm, kmxp, kmym, kmyp
	219
[1257]	220	REAL, DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) :: wsus, wsvs
[1015]	221	!$acc declare create ( wsus, wsvs )
	222
	223	!
	224	!-- First calculate horizontal momentum flux w'u' and/or w'v' at vertical
	225	!-- walls, if neccessary
	226	IF ( topography /= 'flat' ) THEN
	227	CALL wall_fluxes_acc( wsus, 0.0, 0.0, 0.0, 1.0, nzb_w_inner, &
	228	nzb_w_outer, wall_w_x )
	229	CALL wall_fluxes_acc( wsvs, 0.0, 0.0, 1.0, 0.0, nzb_w_inner, &
	230	nzb_w_outer, wall_w_y )
	231	ENDIF
	232
	233	!$acc kernels present ( u, v, w, km, tend, vsws, vswst ) &
	234	!$acc present ( ddzu, ddzw, fwxm, fwxp, fwym, fwyp, wall_w_x, wall_w_y ) &
	235	!$acc present ( nzb_w_inner, nzb_w_outer )
[1128]	236	DO i = i_left, i_right
	237	DO j = j_south, j_north
[1015]	238	DO k = 1, nzt
	239	IF ( k > nzb_w_outer(j,i) ) THEN
	240	!
	241	!-- Interpolate eddy diffusivities on staggered gridpoints
	242	kmxp = 0.25 * &
	243	( km(k,j,i)+km(k,j,i+1)+km(k+1,j,i)+km(k+1,j,i+1) )
	244	kmxm = 0.25 * &
	245	( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) )
	246	kmyp = 0.25 * &
	247	( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) )
	248	kmym = 0.25 * &
	249	( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
	250
	251	tend(k,j,i) = tend(k,j,i) &
	252	& + ( kmxp * ( w(k,j,i+1) - w(k,j,i) ) * ddx &
	253	& + kmxp * ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1) &
	254	& - kmxm * ( w(k,j,i) - w(k,j,i-1) ) * ddx &
	255	& - kmxm * ( u(k+1,j,i) - u(k,j,i) ) * ddzu(k+1) &
	256	& ) * ddx &
	257	& + ( kmyp * ( w(k,j+1,i) - w(k,j,i) ) * ddy &
	258	& + kmyp * ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1) &
	259	& - kmym * ( w(k,j,i) - w(k,j-1,i) ) * ddy &
	260	& - kmym * ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1) &
	261	& ) * ddy &
	262	& + 2.0 * ( &
	263	& km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1) &
	264	& - km(k,j,i) * ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) &
	265	& ) * ddzu(k+1)
	266	ENDIF
	267	ENDDO
	268
	269	!
	270	!-- Wall functions at all vertical walls, where necessary
	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 |