Home

Context Navigation

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

Last change on this file since 4326 was 4182, checked in by scharf, 5 years ago
corrected "Former revisions" section minor formatting in "Former revisions" section added "Author" section
Property svn:keywords set to `Id`
File size: 19.0 KB

Rev	Line
[1873]	1	!> @file diffusion_w.f90
[2000]	2	!------------------------------------------------------------------------------!
[2696]	3	! This file is part of the PALM model system.
[1036]	4	!
[2000]	5	! PALM is free software: you can redistribute it and/or modify it under the
	6	! terms of the GNU General Public License as published by the Free Software
	7	! Foundation, either version 3 of the License, or (at your option) any later
	8	! version.
[1036]	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	!
[3655]	17	! Copyright 1997-2019 Leibniz Universitaet Hannover
[2000]	18	!------------------------------------------------------------------------------!
[1036]	19	!
[484]	20	! Current revisions:
[1]	21	! -----------------
[1341]	22	!
[2233]	23	!
[1321]	24	! Former revisions:
	25	! -----------------
	26	! $Id: diffusion_w.f90 4182 2019-08-22 15:20:23Z oliver.maas $
[4182]	27	! Corrected "Former revisions" section
	28	!
	29	! 3655 2019-01-07 16:51:22Z knoop
[3634]	30	! OpenACC port for SPEC
[1321]	31	!
[4182]	32	! Revision 1.1 1997/09/12 06:24:11 raasch
	33	! Initial revision
	34	!
	35	!
[1]	36	! Description:
	37	! ------------
[1682]	38	!> Diffusion term of the w-component
[1]	39	!------------------------------------------------------------------------------!
[1682]	40	MODULE diffusion_w_mod
	41
[1]	42
	43	PRIVATE
[2118]	44	PUBLIC diffusion_w
[1]	45
	46	INTERFACE diffusion_w
	47	MODULE PROCEDURE diffusion_w
	48	MODULE PROCEDURE diffusion_w_ij
	49	END INTERFACE diffusion_w
	50
	51	CONTAINS
	52
	53
	54	!------------------------------------------------------------------------------!
[1682]	55	! Description:
	56	! ------------
	57	!> Call for all grid points
[1]	58	!------------------------------------------------------------------------------!
[1001]	59	SUBROUTINE diffusion_w
[1]	60
[1320]	61	USE arrays_3d, &
[2037]	62	ONLY : ddzu, ddzw, km, tend, u, v, w, drho_air_zw, rho_air
[1320]	63
	64	USE grid_variables, &
[2232]	65	ONLY : ddx, ddy
[1320]	66
	67	USE indices, &
[2232]	68	ONLY : nxl, nxr, nyn, nys, nzb, nzt, wall_flags_0
[1320]	69
	70	USE kinds
[1]	71
[2232]	72	USE surface_mod, &
	73	ONLY : surf_def_v, surf_lsm_v, surf_usm_v
	74
[1]	75	IMPLICIT NONE
	76
[2232]	77	INTEGER(iwp) :: i !< running index x direction
	78	INTEGER(iwp) :: j !< running index y direction
	79	INTEGER(iwp) :: k !< running index z direction
	80	INTEGER(iwp) :: l !< running index of surface type, south- or north-facing wall
	81	INTEGER(iwp) :: m !< running index surface elements
	82	INTEGER(iwp) :: surf_e !< End index of surface elements at (j,i)-gridpoint
	83	INTEGER(iwp) :: surf_s !< Start index of surface elements at (j,i)-gridpoint
[1320]	84
[2232]	85	REAL(wp) :: flag !< flag to mask topography grid points
[3547]	86	REAL(wp) :: kmxm !< diffusion coefficient on leftward side of the w-gridbox - interpolated onto xu-y grid
	87	REAL(wp) :: kmxp !<diffusion coefficient on rightward side of the w-gridbox - interpolated onto xu-y grid
	88	REAL(wp) :: kmym !< diffusion coefficient on southward side of the w-gridbox - interpolated onto x-yv grid
	89	REAL(wp) :: kmyp !< diffusion coefficient on northward side of the w-gridbox - interpolated onto x-yv grid
[2232]	90	REAL(wp) :: mask_west !< flag to mask vertical wall west of the grid point
	91	REAL(wp) :: mask_east !< flag to mask vertical wall east of the grid point
	92	REAL(wp) :: mask_south !< flag to mask vertical wall south of the grid point
	93	REAL(wp) :: mask_north !< flag to mask vertical wall north of the grid point
[1001]	94
[1]	95
	96
[3634]	97	!$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i, j, k, l, m) &
	98	!$ACC PRIVATE(surf_e, surf_s, flag, kmxm, kmxp, kmym, kmyp) &
	99	!$ACC PRIVATE(mask_west, mask_east, mask_south, mask_north) &
	100	!$ACC PRESENT(wall_flags_0, km) &
	101	!$ACC PRESENT(u, v, w) &
	102	!$ACC PRESENT(ddzu, ddzw, rho_air, drho_air_zw) &
	103	!$ACC PRESENT(surf_def_v(0:3)) &
	104	!$ACC PRESENT(surf_lsm_v(0:3)) &
	105	!$ACC PRESENT(surf_usm_v(0:3)) &
	106	!$ACC PRESENT(tend)
[1]	107	DO i = nxl, nxr
	108	DO j = nys, nyn
[2232]	109	DO k = nzb+1, nzt-1
[1]	110	!
[2232]	111	!-- Predetermine flag to mask topography and wall-bounded grid points.
	112	flag = MERGE( 1.0_wp, 0.0_wp, &
	113	BTEST( wall_flags_0(k,j,i), 3 ) )
	114	mask_east = MERGE( 1.0_wp, 0.0_wp, &
	115	BTEST( wall_flags_0(k,j,i+1), 3 ) )
	116	mask_west = MERGE( 1.0_wp, 0.0_wp, &
	117	BTEST( wall_flags_0(k,j,i-1), 3 ) )
	118	mask_south = MERGE( 1.0_wp, 0.0_wp, &
	119	BTEST( wall_flags_0(k,j-1,i), 3 ) )
	120	mask_north = MERGE( 1.0_wp, 0.0_wp, &
	121	BTEST( wall_flags_0(k,j+1,i), 3 ) )
	122	!
[1]	123	!-- Interpolate eddy diffusivities on staggered gridpoints
[3547]	124	kmxp = 0.25_wp * ( km(k,j,i) + km(k,j,i+1) + &
[2232]	125	km(k+1,j,i) + km(k+1,j,i+1) )
	126	kmxm = 0.25_wp * ( km(k,j,i) + km(k,j,i-1) + &
	127	km(k+1,j,i) + km(k+1,j,i-1) )
	128	kmyp = 0.25_wp * ( km(k,j,i) + km(k+1,j,i) + &
	129	km(k,j+1,i) + km(k+1,j+1,i) )
	130	kmym = 0.25_wp * ( km(k,j,i) + km(k+1,j,i) + &
	131	km(k,j-1,i) + km(k+1,j-1,i) )
[1]	132
	133	tend(k,j,i) = tend(k,j,i) &
[2232]	134	+ ( mask_east * kmxp * ( &
	135	( w(k,j,i+1) - w(k,j,i) ) * ddx &
	136	+ ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1) &
	137	) &
	138	- mask_west * kmxm * ( &
	139	( w(k,j,i) - w(k,j,i-1) ) * ddx &
	140	+ ( u(k+1,j,i) - u(k,j,i) ) * ddzu(k+1) &
	141	) &
	142	) * ddx * flag &
	143	+ ( mask_north * kmyp * ( &
	144	( w(k,j+1,i) - w(k,j,i) ) * ddy &
	145	+ ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1) &
	146	) &
	147	- mask_south * kmym * ( &
	148	( w(k,j,i) - w(k,j-1,i) ) * ddy &
	149	+ ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1) &
	150	) &
	151	) * ddy * flag &
	152	+ 2.0_wp * ( &
	153	km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1) &
	154	* rho_air(k+1) &
	155	- km(k,j,i) * ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) &
	156	* rho_air(k) &
	157	) * ddzu(k+1) * drho_air_zw(k) * flag
[1]	158	ENDDO
	159
	160	!
[2232]	161	!-- Add horizontal momentum flux v'w' at north- (l=0) and south-facing (l=1)
	162	!-- surfaces. Note, in the the flat case, loops won't be entered as
	163	!-- start_index > end_index. Furtermore, note, no vertical natural surfaces
	164	!-- so far.
	165	!-- Default-type surfaces
	166	DO l = 0, 1
	167	surf_s = surf_def_v(l)%start_index(j,i)
	168	surf_e = surf_def_v(l)%end_index(j,i)
	169	DO m = surf_s, surf_e
	170	k = surf_def_v(l)%k(m)
	171	tend(k,j,i) = tend(k,j,i) + &
	172	surf_def_v(l)%mom_flux_w(m) * ddy
	173	ENDDO
	174	ENDDO
[1]	175	!
[2232]	176	!-- Natural-type surfaces
	177	DO l = 0, 1
	178	surf_s = surf_lsm_v(l)%start_index(j,i)
	179	surf_e = surf_lsm_v(l)%end_index(j,i)
	180	DO m = surf_s, surf_e
	181	k = surf_lsm_v(l)%k(m)
	182	tend(k,j,i) = tend(k,j,i) + &
	183	surf_lsm_v(l)%mom_flux_w(m) * ddy
	184	ENDDO
	185	ENDDO
	186	!
	187	!-- Urban-type surfaces
	188	DO l = 0, 1
	189	surf_s = surf_usm_v(l)%start_index(j,i)
	190	surf_e = surf_usm_v(l)%end_index(j,i)
	191	DO m = surf_s, surf_e
	192	k = surf_usm_v(l)%k(m)
	193	tend(k,j,i) = tend(k,j,i) + &
	194	surf_usm_v(l)%mom_flux_w(m) * ddy
	195	ENDDO
	196	ENDDO
	197	!
	198	!-- Add horizontal momentum flux u'w' at east- (l=2) and west-facing (l=3)
	199	!-- surface.
	200	!-- Default-type surfaces
	201	DO l = 2, 3
	202	surf_s = surf_def_v(l)%start_index(j,i)
	203	surf_e = surf_def_v(l)%end_index(j,i)
	204	DO m = surf_s, surf_e
	205	k = surf_def_v(l)%k(m)
	206	tend(k,j,i) = tend(k,j,i) + &
	207	surf_def_v(l)%mom_flux_w(m) * ddx
	208	ENDDO
	209	ENDDO
	210	!
	211	!-- Natural-type surfaces
	212	DO l = 2, 3
	213	surf_s = surf_lsm_v(l)%start_index(j,i)
	214	surf_e = surf_lsm_v(l)%end_index(j,i)
	215	DO m = surf_s, surf_e
	216	k = surf_lsm_v(l)%k(m)
	217	tend(k,j,i) = tend(k,j,i) + &
	218	surf_lsm_v(l)%mom_flux_w(m) * ddx
	219	ENDDO
	220	ENDDO
	221	!
	222	!-- Urban-type surfaces
	223	DO l = 2, 3
	224	surf_s = surf_usm_v(l)%start_index(j,i)
	225	surf_e = surf_usm_v(l)%end_index(j,i)
	226	DO m = surf_s, surf_e
	227	k = surf_usm_v(l)%k(m)
	228	tend(k,j,i) = tend(k,j,i) + &
	229	surf_usm_v(l)%mom_flux_w(m) * ddx
	230	ENDDO
	231	ENDDO
[1]	232
	233	ENDDO
	234	ENDDO
	235
	236	END SUBROUTINE diffusion_w
	237
	238
	239	!------------------------------------------------------------------------------!
[1682]	240	! Description:
	241	! ------------
	242	!> Call for grid point i,j
[1]	243	!------------------------------------------------------------------------------!
[1001]	244	SUBROUTINE diffusion_w_ij( i, j )
[1]	245
[1320]	246	USE arrays_3d, &
[2037]	247	ONLY : ddzu, ddzw, km, tend, u, v, w, drho_air_zw, rho_air
[1320]	248
	249	USE grid_variables, &
[2232]	250	ONLY : ddx, ddy
[1320]	251
	252	USE indices, &
[3241]	253	ONLY : nzb, nzt, wall_flags_0
[1320]	254
	255	USE kinds
[1]	256
[2232]	257	USE surface_mod, &
	258	ONLY : surf_def_v, surf_lsm_v, surf_usm_v
	259
[1]	260	IMPLICIT NONE
	261
[2232]	262
	263	INTEGER(iwp) :: i !< running index x direction
	264	INTEGER(iwp) :: j !< running index y direction
	265	INTEGER(iwp) :: k !< running index z direction
	266	INTEGER(iwp) :: l !< running index of surface type, south- or north-facing wall
	267	INTEGER(iwp) :: m !< running index surface elements
	268	INTEGER(iwp) :: surf_e !< End index of surface elements at (j,i)-gridpoint
	269	INTEGER(iwp) :: surf_s !< Start index of surface elements at (j,i)-gridpoint
[1320]	270
[2232]	271	REAL(wp) :: flag !< flag to mask topography grid points
[3547]	272	REAL(wp) :: kmxm !< diffusion coefficient on leftward side of the w-gridbox - interpolated onto xu-y grid
	273	REAL(wp) :: kmxp !< diffusion coefficient on rightward side of the w-gridbox - interpolated onto xu-y grid
	274	REAL(wp) :: kmym !< diffusion coefficient on southward side of the w-gridbox - interpolated onto x-yv grid
	275	REAL(wp) :: kmyp !< diffusion coefficient on northward side of the w-gridbox - interpolated onto x-yv grid
[2232]	276	REAL(wp) :: mask_west !< flag to mask vertical wall west of the grid point
	277	REAL(wp) :: mask_east !< flag to mask vertical wall east of the grid point
	278	REAL(wp) :: mask_south !< flag to mask vertical wall south of the grid point
	279	REAL(wp) :: mask_north !< flag to mask vertical wall north of the grid point
[1]	280
	281
[2232]	282	DO k = nzb+1, nzt-1
[1]	283	!
[2232]	284	!-- Predetermine flag to mask topography and wall-bounded grid points.
	285	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 3 ) )
	286	mask_east = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i+1), 3 ) )
	287	mask_west = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i-1), 3 ) )
	288	mask_south = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j-1,i), 3 ) )
	289	mask_north = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j+1,i), 3 ) )
	290	!
[1]	291	!-- Interpolate eddy diffusivities on staggered gridpoints
[1322]	292	kmxp = 0.25_wp * ( km(k,j,i)+km(k,j,i+1)+km(k+1,j,i)+km(k+1,j,i+1) )
	293	kmxm = 0.25_wp * ( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) )
	294	kmyp = 0.25_wp * ( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) )
	295	kmym = 0.25_wp * ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
[1]	296
	297	tend(k,j,i) = tend(k,j,i) &
[2232]	298	+ ( mask_east * kmxp * ( &
	299	( w(k,j,i+1) - w(k,j,i) ) * ddx &
	300	+ ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1) &
	301	) &
	302	- mask_west * kmxm * ( &
	303	( w(k,j,i) - w(k,j,i-1) ) * ddx &
	304	+ ( u(k+1,j,i) - u(k,j,i) ) * ddzu(k+1) &
	305	) &
	306	) * ddx * flag &
	307	+ ( mask_north * kmyp * ( &
	308	( w(k,j+1,i) - w(k,j,i) ) * ddy &
	309	+ ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1) &
	310	) &
	311	- mask_south * kmym * ( &
	312	( w(k,j,i) - w(k,j-1,i) ) * ddy &
	313	+ ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1) &
	314	) &
	315	) * ddy * flag &
	316	+ 2.0_wp * ( &
	317	km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1) &
	318	* rho_air(k+1) &
	319	- km(k,j,i) * ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) &
	320	* rho_air(k) &
	321	) * ddzu(k+1) * drho_air_zw(k) * flag
[1]	322	ENDDO
	323	!
[2232]	324	!-- Add horizontal momentum flux v'w' at north- (l=0) and south-facing (l=1)
	325	!-- surfaces. Note, in the the flat case, loops won't be entered as
	326	!-- start_index > end_index. Furtermore, note, no vertical natural surfaces
	327	!-- so far.
	328	!-- Default-type surfaces
	329	DO l = 0, 1
	330	surf_s = surf_def_v(l)%start_index(j,i)
	331	surf_e = surf_def_v(l)%end_index(j,i)
	332	DO m = surf_s, surf_e
	333	k = surf_def_v(l)%k(m)
	334	tend(k,j,i) = tend(k,j,i) + &
	335	surf_def_v(l)%mom_flux_w(m) * ddy
	336	ENDDO
	337	ENDDO
[51]	338	!
[2232]	339	!-- Natural-type surfaces
	340	DO l = 0, 1
	341	surf_s = surf_lsm_v(l)%start_index(j,i)
	342	surf_e = surf_lsm_v(l)%end_index(j,i)
	343	DO m = surf_s, surf_e
	344	k = surf_lsm_v(l)%k(m)
	345	tend(k,j,i) = tend(k,j,i) + &
	346	surf_lsm_v(l)%mom_flux_w(m) * ddy
	347	ENDDO
	348	ENDDO
[1]	349	!
[2232]	350	!-- Urban-type surfaces
	351	DO l = 0, 1
	352	surf_s = surf_usm_v(l)%start_index(j,i)
	353	surf_e = surf_usm_v(l)%end_index(j,i)
	354	DO m = surf_s, surf_e
	355	k = surf_usm_v(l)%k(m)
	356	tend(k,j,i) = tend(k,j,i) + &
	357	surf_usm_v(l)%mom_flux_w(m) * ddy
	358	ENDDO
	359	ENDDO
	360	!
	361	!-- Add horizontal momentum flux u'w' at east- (l=2) and west-facing (l=3)
	362	!-- surfaces.
	363	!-- Default-type surfaces
	364	DO l = 2, 3
	365	surf_s = surf_def_v(l)%start_index(j,i)
	366	surf_e = surf_def_v(l)%end_index(j,i)
	367	DO m = surf_s, surf_e
	368	k = surf_def_v(l)%k(m)
	369	tend(k,j,i) = tend(k,j,i) + &
	370	surf_def_v(l)%mom_flux_w(m) * ddx
	371	ENDDO
	372	ENDDO
	373	!
	374	!-- Natural-type surfaces
	375	DO l = 2, 3
	376	surf_s = surf_lsm_v(l)%start_index(j,i)
	377	surf_e = surf_lsm_v(l)%end_index(j,i)
	378	DO m = surf_s, surf_e
	379	k = surf_lsm_v(l)%k(m)
	380	tend(k,j,i) = tend(k,j,i) + &
	381	surf_lsm_v(l)%mom_flux_w(m) * ddx
	382	ENDDO
	383	ENDDO
	384	!
	385	!-- Urban-type surfaces
	386	DO l = 2, 3
	387	surf_s = surf_usm_v(l)%start_index(j,i)
	388	surf_e = surf_usm_v(l)%end_index(j,i)
	389	DO m = surf_s, surf_e
	390	k = surf_usm_v(l)%k(m)
	391	tend(k,j,i) = tend(k,j,i) + &
	392	surf_usm_v(l)%mom_flux_w(m) * ddx
	393	ENDDO
	394	ENDDO
[1]	395
	396
	397	END SUBROUTINE diffusion_w_ij
	398
[1323]	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 |