Home

Context Navigation

diffusion_w.f90 @ 4901

Last change on this file since 4901 was 4828, checked in by Giersch, 4 years ago
Copyright updated to year 2021, interface pmc_sort removed to accelarate the nesting code
Property svn:keywords set to `Id`
File size: 19.3 KB

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

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

Context Navigation

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

Download in other formats: