Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |