Home

Context Navigation

diffusion_u.f90 @ 4851

Last change on this file since 4851 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: 25.5 KB

Rev	Line
[1873]	1	!> @file diffusion_u.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	! -----------------
[4671]	21	!
	22	!
[1321]	23	! Former revisions:
	24	! -----------------
	25	! $Id: diffusion_u.f90 4828 2021-01-05 11:21:41Z gronemeier $
[4674]	26	! Update ACC directives for downward facing USM and LSM surfaces
	27	!
	28	! 4671 2020-09-09 20:27:58Z pavelkrc
[4671]	29	! Implementation of downward facing USM and LSM surfaces
	30	!
	31	! 4583 2020-06-29 12:36:47Z raasch
[4583]	32	! file re-formatted to follow the PALM coding standard
	33	!
	34	! 4360 2020-01-07 11:25:50Z suehring
	35	! Introduction of wall_flags_total_0, which currently sets bits based on static topography
	36	! information used in wall_flags_static_0
	37	!
[4346]	38	! 4329 2019-12-10 15:46:36Z motisi
[4329]	39	! Renamed wall_flags_0 to wall_flags_static_0
[4583]	40	!
[4329]	41	! 4182 2019-08-22 15:20:23Z scharf
[4182]	42	! Corrected "Former revisions" section
[4583]	43	!
[4182]	44	! 3655 2019-01-07 16:51:22Z knoop
[3634]	45	! OpenACC port for SPEC
[2716]	46	!
[4182]	47	! Revision 1.1 1997/09/12 06:23:51 raasch
	48	! Initial revision
	49	!
	50	!
[1]	51	! Description:
	52	! ------------
[1682]	53	!> Diffusion term of the u-component
[4583]	54	!> @todo additional damping (needed for non-cyclic bc) causes bad vectorization and slows down the
	55	! speed on NEC about 5-10%
	56	!--------------------------------------------------------------------------------------------------!
[1682]	57	MODULE diffusion_u_mod
[1]	58
[4583]	59
[1]	60	PRIVATE
[2118]	61	PUBLIC diffusion_u
[1]	62
	63	INTERFACE diffusion_u
	64	MODULE PROCEDURE diffusion_u
	65	MODULE PROCEDURE diffusion_u_ij
	66	END INTERFACE diffusion_u
	67
	68	CONTAINS
	69
	70
[4583]	71	!--------------------------------------------------------------------------------------------------!
[1682]	72	! Description:
	73	! ------------
	74	!> Call for all grid points
[4583]	75	!--------------------------------------------------------------------------------------------------!
[1001]	76	SUBROUTINE diffusion_u
[1]	77
[4583]	78	USE arrays_3d, &
	79	ONLY: ddzu, ddzw, drho_air, km, rho_air_zw, tend, u, v, w
	80
	81	USE control_parameters, &
	82	ONLY: constant_top_momentumflux, use_surface_fluxes, use_top_fluxes
	83
	84	USE grid_variables, &
[2232]	85	ONLY: ddx, ddx2, ddy
[4583]	86
	87	USE indices, &
[4346]	88	ONLY: nxlu, nxr, nyn, nys, nzb, nzt, wall_flags_total_0
[4583]	89
[1320]	90	USE kinds
[1]	91
[4583]	92	USE surface_mod, &
	93	ONLY : surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, surf_usm_v
[2232]	94
[1]	95	IMPLICIT NONE
	96
[2232]	97	INTEGER(iwp) :: i !< running index x direction
	98	INTEGER(iwp) :: j !< running index y direction
	99	INTEGER(iwp) :: k !< running index z direction
	100	INTEGER(iwp) :: l !< running index of surface type, south- or north-facing wall
	101	INTEGER(iwp) :: m !< running index surface elements
[3547]	102	INTEGER(iwp) :: surf_e !< end index of surface elements at (j,i)-gridpoint
	103	INTEGER(iwp) :: surf_s !< start index of surface elements at (j,i)-gridpoint
[1001]	104
[2232]	105	REAL(wp) :: flag !< flag to mask topography grid points
[3547]	106	REAL(wp) :: kmym !< diffusion coefficient on southward side of the u-gridbox - interpolated onto xu-yv grid
	107	REAL(wp) :: kmyp !< diffusion coefficient on northward side of the u-gridbox - interpolated onto xu-yv grid
	108	REAL(wp) :: kmzm !< diffusion coefficient on bottom of the gridbox - interpolated onto xu-zw grid
	109	REAL(wp) :: kmzp !< diffusion coefficient on top of the gridbox - interpolated onto xu-zw grid
[4583]	110	REAL(wp) :: mask_bottom !< flag to mask vertical upward-facing surface
	111	REAL(wp) :: mask_north !< flag to mask vertical surface north of the grid point
	112	REAL(wp) :: mask_south !< flag to mask vertical surface south of the grid point
	113	REAL(wp) :: mask_top !< flag to mask vertical downward-facing surface
[1]	114
[56]	115
[2232]	116
[3634]	117	!$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i, j, k, l, m) &
	118	!$ACC PRIVATE(surf_e, surf_s, flag, kmym, kmyp, kmzm, kmzp) &
	119	!$ACC PRIVATE(mask_bottom, mask_north, mask_south, mask_top) &
[4346]	120	!$ACC PRESENT(wall_flags_total_0, km) &
[3634]	121	!$ACC PRESENT(u, v, w) &
	122	!$ACC PRESENT(ddzu, ddzw, drho_air, rho_air_zw) &
	123	!$ACC PRESENT(surf_def_h(0:2), surf_def_v(0:1)) &
[4674]	124	!$ACC PRESENT(surf_lsm_h(0:1), surf_lsm_v(0:1)) &
	125	!$ACC PRESENT(surf_usm_h(0:1), surf_usm_v(0:1)) &
[3634]	126	!$ACC PRESENT(tend)
[106]	127	DO i = nxlu, nxr
[1001]	128	DO j = nys, nyn
[1]	129	!
	130	!-- Compute horizontal diffusion
[2232]	131	DO k = nzb+1, nzt
[1]	132	!
[4583]	133	!-- Predetermine flag to mask topography and wall-bounded grid points.
	134	!-- It is sufficient to masked only north- and south-facing surfaces, which need special
	135	!-- treatment for the u-component.
	136	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 1 ) )
[4346]	137	mask_south = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j-1,i), 1 ) )
	138	mask_north = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j+1,i), 1 ) )
[2232]	139	!
[1]	140	!-- Interpolate eddy diffusivities on staggered gridpoints
[4583]	141	kmyp = 0.25_wp * ( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) )
	142	kmym = 0.25_wp * ( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) )
[1]	143
[4583]	144	tend(k,j,i) = tend(k,j,i) &
	145	+ 2.0_wp * ( &
	146	km(k,j,i) * ( u(k,j,i+1) - u(k,j,i) ) &
	147	- km(k,j,i-1) * ( u(k,j,i) - u(k,j,i-1) ) &
	148	) * ddx2 * flag &
	149	+ ( mask_north * ( &
	150	kmyp * ( u(k,j+1,i) - u(k,j,i) ) * ddy &
	151	+ kmyp * ( v(k,j+1,i) - v(k,j+1,i-1) ) * ddx &
	152	) &
	153	- mask_south * ( &
	154	kmym * ( u(k,j,i) - u(k,j-1,i) ) * ddy &
	155	+ kmym * ( v(k,j,i) - v(k,j,i-1) ) * ddx &
	156	) &
	157	) * ddy * flag
[1]	158	ENDDO
	159	!
[4583]	160	!-- Add horizontal momentum flux u'v' at north- (l=0) and south-facing (l=1) surfaces.
	161	!-- Note, in the the flat case, loops won't be entered as start_index > end_index.
	162	!-- Furtermore, note, no vertical natural surfaces so far.
[2232]	163	!-- Default-type surfaces
	164	DO l = 0, 1
	165	surf_s = surf_def_v(l)%start_index(j,i)
	166	surf_e = surf_def_v(l)%end_index(j,i)
	167	DO m = surf_s, surf_e
	168	k = surf_def_v(l)%k(m)
[4583]	169	tend(k,j,i) = tend(k,j,i) + surf_def_v(l)%mom_flux_uv(m) * ddy
	170	ENDDO
[2232]	171	ENDDO
	172	!
	173	!-- Natural-type surfaces
	174	DO l = 0, 1
	175	surf_s = surf_lsm_v(l)%start_index(j,i)
	176	surf_e = surf_lsm_v(l)%end_index(j,i)
	177	DO m = surf_s, surf_e
	178	k = surf_lsm_v(l)%k(m)
[4583]	179	tend(k,j,i) = tend(k,j,i) + surf_lsm_v(l)%mom_flux_uv(m) * ddy
	180	ENDDO
[2232]	181	ENDDO
	182	!
	183	!-- Urban-type surfaces
	184	DO l = 0, 1
	185	surf_s = surf_usm_v(l)%start_index(j,i)
	186	surf_e = surf_usm_v(l)%end_index(j,i)
	187	DO m = surf_s, surf_e
	188	k = surf_usm_v(l)%k(m)
[4583]	189	tend(k,j,i) = tend(k,j,i) + surf_usm_v(l)%mom_flux_uv(m) * ddy
	190	ENDDO
[2232]	191	ENDDO
[51]	192
[1]	193	!
[4583]	194	!-- Compute vertical diffusion. In case of simulating a surface layer, respective grid
	195	!-- diffusive fluxes are masked (flag 8) within this loop, and added further below, else,
	196	!-- simple gradient approach is applied. Model top is also mask if top-momentum flux is
	197	!-- given.
[2232]	198	DO k = nzb+1, nzt
[1]	199	!
[4583]	200	!-- Determine flags to mask topography below and above. Flag 1 is used to mask
	201	!-- topography in general, and flag 8 implies information about use_surface_fluxes.
	202	!-- Flag 9 is used to control momentum flux at model top.
	203	mask_bottom = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k-1,j,i), 8 ) )
	204	mask_top = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k+1,j,i), 8 ) ) * &
	205	MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k+1,j,i), 9 ) )
	206	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 1 ) )
[2232]	207	!
[1]	208	!-- Interpolate eddy diffusivities on staggered gridpoints
[4583]	209	kmzp = 0.25_wp * ( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) )
	210	kmzm = 0.25_wp * ( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
[1]	211
[4583]	212	tend(k,j,i) = tend(k,j,i) &
	213	+ ( kmzp * ( ( u(k+1,j,i) - u(k,j,i) ) * ddzu(k+1) &
	214	+ ( w(k,j,i) - w(k,j,i-1) ) * ddx &
	215	) * rho_air_zw(k) * mask_top &
	216	- kmzm * ( ( u(k,j,i) - u(k-1,j,i) ) * ddzu(k) &
	217	+ ( w(k-1,j,i) - w(k-1,j,i-1) ) * ddx &
	218	) * rho_air_zw(k-1) * mask_bottom &
	219	) * ddzw(k) * drho_air(k) * flag
[1]	220	ENDDO
	221
	222	!
[4583]	223	!-- Vertical diffusion at the first grid point above the surface, if the momentum flux at
	224	!-- the bottom is given by the Prandtl law or if it is prescribed by the user.
	225	!-- Difference quotient of the momentum flux is not formed over half of the grid spacing
	226	!-- (2.0*ddzw(k)) any more, since the comparison with other (LES) models showed that the
	227	!-- values of the momentum flux becomes too large in this case.
	228	!-- The term containing w(k-1,..) (see above equation) is removed here because the vertical
	229	!-- velocity is assumed to be zero at the surface.
[1]	230	IF ( use_surface_fluxes ) THEN
	231	!
[2232]	232	!-- Default-type surfaces, upward-facing
	233	surf_s = surf_def_h(0)%start_index(j,i)
	234	surf_e = surf_def_h(0)%end_index(j,i)
	235	DO m = surf_s, surf_e
[1]	236
[2232]	237	k = surf_def_h(0)%k(m)
[1]	238
[4583]	239	tend(k,j,i) = tend(k,j,i) &
	240	+ ( - ( - surf_def_h(0)%usws(m) ) ) * ddzw(k) * drho_air(k)
[2232]	241	ENDDO
[102]	242	!
[2232]	243	!-- Default-type surfaces, dowward-facing
	244	surf_s = surf_def_h(1)%start_index(j,i)
	245	surf_e = surf_def_h(1)%end_index(j,i)
	246	DO m = surf_s, surf_e
	247
	248	k = surf_def_h(1)%k(m)
	249
[4583]	250	tend(k,j,i) = tend(k,j,i) &
	251	+ ( - surf_def_h(1)%usws(m) ) * ddzw(k) * drho_air(k)
[2232]	252	ENDDO
[102]	253	!
[2232]	254	!-- Natural-type surfaces, upward-facing
[4671]	255	surf_s = surf_lsm_h(0)%start_index(j,i)
	256	surf_e = surf_lsm_h(0)%end_index(j,i)
[2232]	257	DO m = surf_s, surf_e
[102]	258
[4671]	259	k = surf_lsm_h(0)%k(m)
[2232]	260
[4583]	261	tend(k,j,i) = tend(k,j,i) &
[4671]	262	+ ( - ( - surf_lsm_h(0)%usws(m) ) ) * ddzw(k) * drho_air(k)
[2232]	263	ENDDO
	264	!
[4671]	265	!-- Natural-type surfaces, downward-facing
	266	surf_s = surf_lsm_h(1)%start_index(j,i)
	267	surf_e = surf_lsm_h(1)%end_index(j,i)
	268	DO m = surf_s, surf_e
	269
	270	k = surf_lsm_h(1)%k(m)
	271
	272	tend(k,j,i) = tend(k,j,i) &
	273	+ ( - surf_lsm_h(1)%usws(m) ) * ddzw(k) * drho_air(k)
	274	ENDDO
	275	!
[2232]	276	!-- Urban-type surfaces, upward-facing
[4671]	277	surf_s = surf_usm_h(0)%start_index(j,i)
	278	surf_e = surf_usm_h(0)%end_index(j,i)
[2232]	279	DO m = surf_s, surf_e
	280
[4671]	281	k = surf_usm_h(0)%k(m)
[2232]	282
[4583]	283	tend(k,j,i) = tend(k,j,i) &
[4671]	284	+ ( - ( - surf_usm_h(0)%usws(m) ) ) * ddzw(k) * drho_air(k)
[2232]	285	ENDDO
[4671]	286	!
	287	!-- Urban-type surfaces, downward-facing
	288	surf_s = surf_usm_h(1)%start_index(j,i)
	289	surf_e = surf_usm_h(1)%end_index(j,i)
	290	DO m = surf_s, surf_e
[2232]	291
[4671]	292	k = surf_usm_h(1)%k(m)
	293
	294	tend(k,j,i) = tend(k,j,i) &
	295	+ ( - surf_usm_h(1)%usws(m) ) * ddzw(k) * drho_air(k)
	296	ENDDO
	297
[102]	298	ENDIF
[2232]	299	!
	300	!-- Add momentum flux at model top
[2638]	301	IF ( use_top_fluxes .AND. constant_top_momentumflux ) THEN
[2232]	302	surf_s = surf_def_h(2)%start_index(j,i)
	303	surf_e = surf_def_h(2)%end_index(j,i)
	304	DO m = surf_s, surf_e
[102]	305
[2232]	306	k = surf_def_h(2)%k(m)
	307
[4583]	308	tend(k,j,i) = tend(k,j,i) &
	309	+ ( - surf_def_h(2)%usws(m) ) * ddzw(k) * drho_air(k)
[2232]	310	ENDDO
	311	ENDIF
	312
[1]	313	ENDDO
	314	ENDDO
	315
	316	END SUBROUTINE diffusion_u
	317
	318
[4583]	319	!--------------------------------------------------------------------------------------------------!
[1682]	320	! Description:
	321	! ------------
	322	!> Call for grid point i,j
[4583]	323	!--------------------------------------------------------------------------------------------------!
[1001]	324	SUBROUTINE diffusion_u_ij( i, j )
[1]	325
[4583]	326	USE arrays_3d, &
	327	ONLY: ddzu, ddzw, drho_air, km, tend, u, v, w, rho_air_zw
	328
	329	USE control_parameters, &
	330	ONLY: constant_top_momentumflux, use_surface_fluxes, use_top_fluxes
	331
	332	USE grid_variables, &
[2232]	333	ONLY: ddx, ddx2, ddy
[4583]	334
	335	USE indices, &
[4346]	336	ONLY: nzb, nzt, wall_flags_total_0
[4583]	337
[1320]	338	USE kinds
[1]	339
[4583]	340	USE surface_mod, &
	341	ONLY : surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, surf_usm_v
[2232]	342
[1]	343	IMPLICIT NONE
	344
[2232]	345	INTEGER(iwp) :: i !< running index x direction
	346	INTEGER(iwp) :: j !< running index y direction
	347	INTEGER(iwp) :: k !< running index z direction
	348	INTEGER(iwp) :: l !< running index of surface type, south- or north-facing wall
	349	INTEGER(iwp) :: m !< running index surface elements
	350	INTEGER(iwp) :: surf_e !< End index of surface elements at (j,i)-gridpoint
	351	INTEGER(iwp) :: surf_s !< Start index of surface elements at (j,i)-gridpoint
[1]	352
[2232]	353	REAL(wp) :: flag !< flag to mask topography grid points
[3547]	354	REAL(wp) :: kmym !< diffusion coefficient on southward side of the u-gridbox - interpolated onto xu-yv grid
	355	REAL(wp) :: kmyp !<diffusion coefficient on northward side of the u-gridbox - interpolated onto xu-yv grid
	356	REAL(wp) :: kmzm !< diffusion coefficient on bottom of the gridbox - interpolated onto xu-zw grid
	357	REAL(wp) :: kmzp !< diffusion coefficient on top of the gridbox - interpolated onto xu-zw grid
[4583]	358	REAL(wp) :: mask_bottom !< flag to mask vertical upward-facing surface
	359	REAL(wp) :: mask_north !< flag to mask vertical surface north of the grid point
	360	REAL(wp) :: mask_south !< flag to mask vertical surface south of the grid point
	361	REAL(wp) :: mask_top !< flag to mask vertical downward-facing surface
	362	!
[1]	363	!-- Compute horizontal diffusion
[2232]	364	DO k = nzb+1, nzt
[1]	365	!
[4583]	366	!-- Predetermine flag to mask topography and wall-bounded grid points.
	367	!-- It is sufficient to masked only north- and south-facing surfaces, which need special
	368	!-- treatment for the u-component.
	369	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 1 ) )
[4346]	370	mask_south = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j-1,i), 1 ) )
	371	mask_north = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j+1,i), 1 ) )
[2232]	372	!
[1]	373	!-- Interpolate eddy diffusivities on staggered gridpoints
[1340]	374	kmyp = 0.25_wp * ( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) )
	375	kmym = 0.25_wp * ( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) )
[1]	376
[4583]	377	tend(k,j,i) = tend(k,j,i) &
	378	+ 2.0_wp * ( &
	379	km(k,j,i) * ( u(k,j,i+1) - u(k,j,i) ) &
	380	- km(k,j,i-1) * ( u(k,j,i) - u(k,j,i-1) ) &
	381	) * ddx2 * flag &
	382	+ ( &
	383	mask_north * kmyp * ( ( u(k,j+1,i) - u(k,j,i) ) * ddy &
	384	+ ( v(k,j+1,i) - v(k,j+1,i-1) ) * ddx &
	385	) &
	386	- mask_south * kmym * ( ( u(k,j,i) - u(k,j-1,i) ) * ddy &
	387	+ ( v(k,j,i) - v(k,j,i-1) ) * ddx &
	388	) &
	389	) * ddy * flag
[1]	390	ENDDO
	391
	392	!
[4583]	393	!-- Add horizontal momentum flux u'v' at north- (l=0) and south-facing (l=1) surfaces. Note, in
	394	!-- the the flat case, loops won't be entered as start_index > end_index. Furtermore, note, no
	395	!-- vertical natural surfaces so far.
[2232]	396	!-- Default-type surfaces
	397	DO l = 0, 1
	398	surf_s = surf_def_v(l)%start_index(j,i)
	399	surf_e = surf_def_v(l)%end_index(j,i)
	400	DO m = surf_s, surf_e
	401	k = surf_def_v(l)%k(m)
	402	tend(k,j,i) = tend(k,j,i) + surf_def_v(l)%mom_flux_uv(m) * ddy
[4583]	403	ENDDO
[2232]	404	ENDDO
[51]	405	!
[2232]	406	!-- Natural-type surfaces
	407	DO l = 0, 1
	408	surf_s = surf_lsm_v(l)%start_index(j,i)
	409	surf_e = surf_lsm_v(l)%end_index(j,i)
	410	DO m = surf_s, surf_e
	411	k = surf_lsm_v(l)%k(m)
	412	tend(k,j,i) = tend(k,j,i) + surf_lsm_v(l)%mom_flux_uv(m) * ddy
[4583]	413	ENDDO
[2232]	414	ENDDO
[1]	415	!
[2232]	416	!-- Urban-type surfaces
	417	DO l = 0, 1
	418	surf_s = surf_usm_v(l)%start_index(j,i)
	419	surf_e = surf_usm_v(l)%end_index(j,i)
	420	DO m = surf_s, surf_e
	421	k = surf_usm_v(l)%k(m)
	422	tend(k,j,i) = tend(k,j,i) + surf_usm_v(l)%mom_flux_uv(m) * ddy
[4583]	423	ENDDO
[2232]	424	ENDDO
[1]	425	!
[4583]	426	!-- Compute vertical diffusion. In case of simulating a surface layer, respective grid diffusive
	427	!-- fluxes are masked (flag 8) within this loop, and added further below, else, simple gradient
	428	!-- approach is applied. Model top is also mask if top-momentum flux is given.
[2232]	429	DO k = nzb+1, nzt
	430	!
[4583]	431	!-- Determine flags to mask topography below and above. Flag 1 is used to mask topography in
	432	!-- general, and flag 8 implies information about use_surface_fluxes. Flag 9 is used to
	433	!-- control momentum flux at model top.
	434	mask_bottom = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k-1,j,i), 8 ) )
	435	mask_top = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k+1,j,i), 8 ) ) * &
	436	MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k+1,j,i), 9 ) )
	437	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 1 ) )
[2232]	438	!
[1]	439	!-- Interpolate eddy diffusivities on staggered gridpoints
[1340]	440	kmzp = 0.25_wp * ( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) )
	441	kmzm = 0.25_wp * ( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
[1]	442
[4583]	443	tend(k,j,i) = tend(k,j,i) &
	444	+ ( kmzp * ( ( u(k+1,j,i) - u(k,j,i) ) * ddzu(k+1) &
	445	+ ( w(k,j,i) - w(k,j,i-1) ) * ddx &
	446	) * rho_air_zw(k) * mask_top &
	447	- kmzm * ( ( u(k,j,i) - u(k-1,j,i) ) * ddzu(k) &
	448	+ ( w(k-1,j,i) - w(k-1,j,i-1) ) * ddx &
	449	) * rho_air_zw(k-1) * mask_bottom &
[2232]	450	) * ddzw(k) * drho_air(k) * flag
[1]	451	ENDDO
	452
	453	!
[4583]	454	!-- Vertical diffusion at the first surface grid points, if the momentum flux at the bottom is
	455	!-- given by the Prandtl law or if it is prescribed by the user.
	456	!-- Difference quotient of the momentum flux is not formed over half of the grid spacing
	457	!-- (2.0*ddzw(k)) any more, since the comparison with other (LES) models showed that the values
	458	!-- of the momentum flux becomes too large in this case.
[1]	459	IF ( use_surface_fluxes ) THEN
	460	!
[2232]	461	!-- Default-type surfaces, upward-facing
	462	surf_s = surf_def_h(0)%start_index(j,i)
	463	surf_e = surf_def_h(0)%end_index(j,i)
	464	DO m = surf_s, surf_e
[1]	465
[2232]	466	k = surf_def_h(0)%k(m)
[1]	467
[4583]	468	tend(k,j,i) = tend(k,j,i) + ( - ( - surf_def_h(0)%usws(m) ) ) * ddzw(k) * drho_air(k)
[2232]	469	ENDDO
[102]	470	!
[2232]	471	!-- Default-type surfaces, dowward-facing (except for model-top fluxes)
	472	surf_s = surf_def_h(1)%start_index(j,i)
	473	surf_e = surf_def_h(1)%end_index(j,i)
	474	DO m = surf_s, surf_e
	475
	476	k = surf_def_h(1)%k(m)
	477
[4583]	478	tend(k,j,i) = tend(k,j,i) + ( - surf_def_h(1)%usws(m) ) * ddzw(k) * drho_air(k)
[2232]	479	ENDDO
[102]	480	!
[2232]	481	!-- Natural-type surfaces, upward-facing
[4671]	482	surf_s = surf_lsm_h(0)%start_index(j,i)
	483	surf_e = surf_lsm_h(0)%end_index(j,i)
[2232]	484	DO m = surf_s, surf_e
[102]	485
[4671]	486	k = surf_lsm_h(0)%k(m)
[2232]	487
[4671]	488	tend(k,j,i) = tend(k,j,i) + ( - ( - surf_lsm_h(0)%usws(m) ) ) * ddzw(k) * drho_air(k)
[2232]	489	ENDDO
	490	!
[4671]	491	!-- Natural-type surfaces, downward-facing
	492	surf_s = surf_lsm_h(1)%start_index(j,i)
	493	surf_e = surf_lsm_h(1)%end_index(j,i)
	494	DO m = surf_s, surf_e
	495
	496	k = surf_lsm_h(1)%k(m)
	497
	498	tend(k,j,i) = tend(k,j,i) + ( - surf_lsm_h(1)%usws(m) ) * ddzw(k) * drho_air(k)
	499	ENDDO
	500	!
[2232]	501	!-- Urban-type surfaces, upward-facing
[4671]	502	surf_s = surf_usm_h(0)%start_index(j,i)
	503	surf_e = surf_usm_h(0)%end_index(j,i)
[2232]	504	DO m = surf_s, surf_e
	505
[4671]	506	k = surf_usm_h(0)%k(m)
[2232]	507
[4671]	508	tend(k,j,i) = tend(k,j,i) + ( - ( - surf_usm_h(0)%usws(m) ) ) * ddzw(k) * drho_air(k)
[2232]	509	ENDDO
[4671]	510	!
	511	!-- Urban-type surfaces, downward-facing
	512	surf_s = surf_usm_h(1)%start_index(j,i)
	513	surf_e = surf_usm_h(1)%end_index(j,i)
	514	DO m = surf_s, surf_e
[2232]	515
[4671]	516	k = surf_usm_h(1)%k(m)
	517
	518	tend(k,j,i) = tend(k,j,i) + ( - surf_usm_h(1)%usws(m) ) * ddzw(k) * drho_air(k)
	519	ENDDO
	520
[102]	521	ENDIF
[2232]	522	!
	523	!-- Add momentum flux at model top
[2638]	524	IF ( use_top_fluxes .AND. constant_top_momentumflux ) THEN
[2232]	525	surf_s = surf_def_h(2)%start_index(j,i)
	526	surf_e = surf_def_h(2)%end_index(j,i)
	527	DO m = surf_s, surf_e
[102]	528
[2232]	529	k = surf_def_h(2)%k(m)
	530
[4583]	531	tend(k,j,i) = tend(k,j,i) + ( - surf_def_h(2)%usws(m) ) * ddzw(k) * drho_air(k)
[2232]	532	ENDDO
	533	ENDIF
	534
	535
[1]	536	END SUBROUTINE diffusion_u_ij
	537
	538	END MODULE diffusion_u_mod

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

Context Navigation

source: palm/trunk/SOURCE/diffusion_u.f90 @ 4851

Download in other formats: