Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |