Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |