Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |