Home

Context Navigation

source: palm/trunk/SOURCE/wall_fluxes.f90 @ 1834

Last change on this file since 1834 was 1818, checked in by maronga, 8 years ago
last commit documented / copyright update
Property svn:keywords set to `Id`
File size: 43.3 KB

Rev	Line
[1682]	1	!> @file wall_fluxes.f90
[1036]	2	!--------------------------------------------------------------------------------!
	3	! This file is part of PALM.
	4	!
	5	! PALM is free software: you can redistribute it and/or modify it under the terms
	6	! of the GNU General Public License as published by the Free Software Foundation,
	7	! either version 3 of the License, or (at your option) any later version.
	8	!
	9	! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
	10	! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
	11	! A PARTICULAR PURPOSE. See the GNU General Public License for more details.
	12	!
	13	! You should have received a copy of the GNU General Public License along with
	14	! PALM. If not, see <http://www.gnu.org/licenses/>.
	15	!
[1818]	16	! Copyright 1997-2016 Leibniz Universitaet Hannover
[1036]	17	!--------------------------------------------------------------------------------!
	18	!
[484]	19	! Current revisions:
[52]	20	! -----------------
[1354]	21	!
[1692]	22	!
[1321]	23	! Former revisions:
	24	! -----------------
	25	! $Id: wall_fluxes.f90 1818 2016-04-06 15:53:27Z raasch $
	26	!
[1692]	27	! 1691 2015-10-26 16:17:44Z maronga
	28	! Renamed rif_min and rif_max with zeta_min and zeta_max, respectively.
	29	!
[1683]	30	! 1682 2015-10-07 23:56:08Z knoop
	31	! Code annotations made doxygen readable
	32	!
[1375]	33	! 1374 2014-04-25 12:55:07Z raasch
	34	! pt removed from acc-present-list
	35	!
[1354]	36	! 1353 2014-04-08 15:21:23Z heinze
	37	! REAL constants provided with KIND-attribute
	38	!
[1321]	39	! 1320 2014-03-20 08:40:49Z raasch
[1320]	40	! ONLY-attribute added to USE-statements,
	41	! kind-parameters added to all INTEGER and REAL declaration statements,
	42	! kinds are defined in new module kinds,
	43	! old module precision_kind is removed,
	44	! revision history before 2012 removed,
	45	! comment fields (!:) to be used for variable explanations added to
	46	! all variable declaration statements
[198]	47	!
[1258]	48	! 1257 2013-11-08 15:18:40Z raasch
	49	! openacc loop and loop vector clauses removed
	50	!
[1154]	51	! 1153 2013-05-10 14:33:08Z raasch
	52	! code adjustments of accelerator version required by PGI 12.3 / CUDA 5.0
	53	!
[1132]	54	! 1128 2013-04-12 06:19:32Z raasch
	55	! loop index bounds in accelerator version replaced by i_left, i_right, j_south,
	56	! j_north
	57	!
[1037]	58	! 1036 2012-10-22 13:43:42Z raasch
	59	! code put under GPL (PALM 3.9)
	60	!
[1017]	61	! 1015 2012-09-27 09:23:24Z raasch
	62	! accelerator version (*_acc) added
	63	!
[52]	64	! Initial version (2007/03/07)
	65	!
	66	! Description:
	67	! ------------
[1682]	68	!> Calculates momentum fluxes at vertical walls assuming Monin-Obukhov
	69	!> similarity.
	70	!> Indices: usvs a=1, vsus b=1, wsvs c1=1, wsus c2=1 (other=0).
	71	!> The all-gridpoint version of wall_fluxes_e is not used so far, because
	72	!> it gives slightly different results from the ij-version for some unknown
	73	!> reason.
[1691]	74	!>
	75	!> @todo Rename rif to zeta throughout the routine
[52]	76	!------------------------------------------------------------------------------!
[1682]	77	MODULE wall_fluxes_mod
	78
[56]	79	PRIVATE
[1015]	80	PUBLIC wall_fluxes, wall_fluxes_acc, wall_fluxes_e, wall_fluxes_e_acc
[56]	81
	82	INTERFACE wall_fluxes
	83	MODULE PROCEDURE wall_fluxes
	84	MODULE PROCEDURE wall_fluxes_ij
	85	END INTERFACE wall_fluxes
	86
[1015]	87	INTERFACE wall_fluxes_acc
	88	MODULE PROCEDURE wall_fluxes_acc
	89	END INTERFACE wall_fluxes_acc
	90
[56]	91	INTERFACE wall_fluxes_e
	92	MODULE PROCEDURE wall_fluxes_e
	93	MODULE PROCEDURE wall_fluxes_e_ij
	94	END INTERFACE wall_fluxes_e
	95
[1015]	96	INTERFACE wall_fluxes_e_acc
	97	MODULE PROCEDURE wall_fluxes_e_acc
	98	END INTERFACE wall_fluxes_e_acc
	99
[56]	100	CONTAINS
[52]	101
[56]	102	!------------------------------------------------------------------------------!
[1682]	103	! Description:
	104	! ------------
	105	!> Call for all grid points
[56]	106	!------------------------------------------------------------------------------!
[1320]	107	SUBROUTINE wall_fluxes( wall_flux, a, b, c1, c2, nzb_uvw_inner, &
[56]	108	nzb_uvw_outer, wall )
[52]	109
[1320]	110	USE arrays_3d, &
	111	ONLY: rif_wall, u, v, w, z0, pt
	112
	113	USE control_parameters, &
[1691]	114	ONLY: g, kappa, zeta_max, zeta_min
[1320]	115
	116	USE grid_variables, &
	117	ONLY: dx, dy
	118
	119	USE indices, &
	120	ONLY: nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb, nzt
	121
	122	USE kinds
	123
	124	USE statistics, &
	125	ONLY: hom
[52]	126
[56]	127	IMPLICIT NONE
[52]	128
[1682]	129	INTEGER(iwp) :: i !<
	130	INTEGER(iwp) :: j !<
	131	INTEGER(iwp) :: k !<
	132	INTEGER(iwp) :: wall_index !<
[52]	133
[1320]	134	INTEGER(iwp), &
	135	DIMENSION(nysg:nyng,nxlg:nxrg) :: &
[1682]	136	nzb_uvw_inner !<
[1320]	137	INTEGER(iwp), &
	138	DIMENSION(nysg:nyng,nxlg:nxrg) :: &
[1682]	139	nzb_uvw_outer !<
[1320]	140
[1682]	141	REAL(wp) :: a !<
	142	REAL(wp) :: b !<
	143	REAL(wp) :: c1 !<
	144	REAL(wp) :: c2 !<
	145	REAL(wp) :: h1 !<
	146	REAL(wp) :: h2 !<
	147	REAL(wp) :: zp !<
	148	REAL(wp) :: pts !<
	149	REAL(wp) :: pt_i !<
	150	REAL(wp) :: rifs !<
	151	REAL(wp) :: u_i !<
	152	REAL(wp) :: v_i !<
	153	REAL(wp) :: us_wall !<
	154	REAL(wp) :: vel_total !<
	155	REAL(wp) :: ws !<
	156	REAL(wp) :: wspts !<
[52]	157
[1320]	158	REAL(wp), &
	159	DIMENSION(nysg:nyng,nxlg:nxrg) :: &
[1682]	160	wall !<
[1320]	161
	162	REAL(wp), &
	163	DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) :: &
[1682]	164	wall_flux !<
[52]	165
	166
[1353]	167	zp = 0.5_wp * ( (a+c1) * dy + (b+c2) * dx )
	168	wall_flux = 0.0_wp
[56]	169	wall_index = NINT( a+ 2b + 3c1 + 4*c2 )
	170
[75]	171	DO i = nxl, nxr
	172	DO j = nys, nyn
[56]	173
[1353]	174	IF ( wall(j,i) /= 0.0_wp ) THEN
[52]	175	!
[56]	176	!-- All subsequent variables are computed for the respective
[187]	177	!-- location where the respective flux is defined.
[56]	178	DO k = nzb_uvw_inner(j,i)+1, nzb_uvw_outer(j,i)
[53]	179
[52]	180	!
[56]	181	!-- (1) Compute rifs, u_i, v_i, ws, pt' and w'pt'
	182	rifs = rif_wall(k,j,i,wall_index)
[53]	183
[1353]	184	u_i = a * u(k,j,i) + c1 * 0.25_wp * &
[56]	185	( u(k+1,j,i+1) + u(k+1,j,i) + u(k,j,i+1) + u(k,j,i) )
[53]	186
[1353]	187	v_i = b * v(k,j,i) + c2 * 0.25_wp * &
[56]	188	( v(k+1,j+1,i) + v(k+1,j,i) + v(k,j+1,i) + v(k,j,i) )
[53]	189
[1353]	190	ws = ( c1 + c2 ) * w(k,j,i) + 0.25_wp * ( &
[56]	191	a * ( w(k-1,j,i-1) + w(k-1,j,i) + w(k,j,i-1) + w(k,j,i) ) &
	192	+ b * ( w(k-1,j-1,i) + w(k-1,j,i) + w(k,j-1,i) + w(k,j,i) ) &
[1353]	193	)
	194	pt_i = 0.5_wp * ( pt(k,j,i) + a * pt(k,j,i-1) + &
[56]	195	b * pt(k,j-1,i) + ( c1 + c2 ) * pt(k+1,j,i) )
[53]	196
[56]	197	pts = pt_i - hom(k,1,4,0)
	198	wspts = ws * pts
[53]	199
[52]	200	!
[56]	201	!-- (2) Compute wall-parallel absolute velocity vel_total
	202	vel_total = SQRT( ws*2 + (a+c1) u_i*2 + (b+c2) v_i**2 )
[53]	203
[52]	204	!
[56]	205	!-- (3) Compute wall friction velocity us_wall
[1353]	206	IF ( rifs >= 0.0_wp ) THEN
[53]	207
[52]	208	!
[56]	209	!-- Stable stratification (and neutral)
	210	us_wall = kappa * vel_total / ( LOG( zp / z0(j,i) ) + &
[1353]	211	5.0_wp * rifs * ( zp - z0(j,i) ) / zp &
[56]	212	)
	213	ELSE
[53]	214
[52]	215	!
[56]	216	!-- Unstable stratification
[1353]	217	h1 = SQRT( SQRT( 1.0_wp - 16.0_wp * rifs ) )
	218	h2 = SQRT( SQRT( 1.0_wp - 16.0_wp * rifs * z0(j,i) / zp ) )
[53]	219
[187]	220	us_wall = kappa * vel_total / ( &
	221	LOG( zp / z0(j,i) ) - &
[1353]	222	LOG( ( 1.0_wp + h1 )*2 ( 1.0_wp + h1**2 ) / ( &
	223	( 1.0_wp + h2 )*2 ( 1.0_wp + h2**2 ) ) ) +&
	224	2.0_wp * ( ATAN( h1 ) - ATAN( h2 ) ) &
[187]	225	)
[56]	226	ENDIF
[53]	227
[52]	228	!
[56]	229	!-- (4) Compute zp/L (corresponds to neutral Richardson flux
	230	!-- number rifs)
[1353]	231	rifs = -1.0_wp * zp * kappa * g * wspts / &
	232	( pt_i * ( us_wall**3 + 1E-30 ) )
[53]	233
[52]	234	!
[56]	235	!-- Limit the value range of the Richardson numbers.
	236	!-- This is necessary for very small velocities (u,w --> 0),
	237	!-- because the absolute value of rif can then become very
	238	!-- large, which in consequence would result in very large
	239	!-- shear stresses and very small momentum fluxes (both are
	240	!-- generally unrealistic).
[1691]	241	IF ( rifs < zeta_min ) rifs = zeta_min
	242	IF ( rifs > zeta_max ) rifs = zeta_max
[53]	243
[52]	244	!
[56]	245	!-- (5) Compute wall_flux (u'v', v'u', w'v', or w'u')
[1353]	246	IF ( rifs >= 0.0_wp ) THEN
[53]	247
[52]	248	!
[56]	249	!-- Stable stratification (and neutral)
	250	wall_flux(k,j,i) = kappa * &
	251	( au(k,j,i) + bv(k,j,i) + (c1+c2)*w(k,j,i) ) / &
	252	( LOG( zp / z0(j,i) ) + &
[1353]	253	5.0_wp * rifs * ( zp - z0(j,i) ) / zp &
[56]	254	)
	255	ELSE
[53]	256
[52]	257	!
[56]	258	!-- Unstable stratification
[1353]	259	h1 = SQRT( SQRT( 1.0_wp - 16.0_wp * rifs ) )
	260	h2 = SQRT( SQRT( 1.0_wp - 16.0_wp * rifs * z0(j,i) / zp ) )
[53]	261
[187]	262	wall_flux(k,j,i) = kappa * &
	263	( au(k,j,i) + bv(k,j,i) + (c1+c2)*w(k,j,i) ) / ( &
	264	LOG( zp / z0(j,i) ) - &
[1353]	265	LOG( ( 1.0_wp + h1 )*2 ( 1.0_wp + h1**2 ) / ( &
	266	( 1.0_wp + h2 )*2 ( 1.0_wp + h2**2 ) ) ) +&
	267	2.0_wp * ( ATAN( h1 ) - ATAN( h2 ) ) &
[187]	268	)
[56]	269	ENDIF
[187]	270	wall_flux(k,j,i) = -wall_flux(k,j,i) * us_wall
[56]	271
	272	!
	273	!-- store rifs for next time step
	274	rif_wall(k,j,i,wall_index) = rifs
	275
	276	ENDDO
	277
	278	ENDIF
	279
	280	ENDDO
	281	ENDDO
	282
	283	END SUBROUTINE wall_fluxes
	284
	285
[1015]	286	!------------------------------------------------------------------------------!
[1682]	287	! Description:
	288	! ------------
	289	!> Call for all grid points - accelerator version
[1015]	290	!------------------------------------------------------------------------------!
[1320]	291	SUBROUTINE wall_fluxes_acc( wall_flux, a, b, c1, c2, nzb_uvw_inner, &
[1015]	292	nzb_uvw_outer, wall )
[56]	293
[1320]	294	USE arrays_3d, &
	295	ONLY: rif_wall, pt, u, v, w, z0
	296
	297	USE control_parameters, &
[1691]	298	ONLY: g, kappa, zeta_max, zeta_min
[1320]	299
	300	USE grid_variables, &
	301	ONLY: dx, dy
	302
	303	USE indices, &
	304	ONLY: i_left, i_right, j_north, j_south, nxl, nxlg, nxr, nxrg, &
	305	nyn, nyng, nys, nysg, nzb, nzt
	306
	307	USE kinds
	308
	309	USE statistics, &
	310	ONLY: hom
[1015]	311
	312	IMPLICIT NONE
	313
[1682]	314	INTEGER(iwp) :: i !<
	315	INTEGER(iwp) :: j !<
	316	INTEGER(iwp) :: k !<
	317	INTEGER(iwp) :: max_outer !<
	318	INTEGER(iwp) :: min_inner !<
	319	INTEGER(iwp) :: wall_index !<
[1015]	320
[1320]	321	INTEGER(iwp), &
	322	DIMENSION(nysg:nyng,nxlg:nxrg) :: &
[1682]	323	nzb_uvw_inner !<
[1320]	324	INTEGER(iwp), &
	325	DIMENSION(nysg:nyng,nxlg:nxrg) :: &
[1682]	326	nzb_uvw_outer !<
[1320]	327
[1682]	328	REAL(wp) :: a !<
	329	REAL(wp) :: b !<
	330	REAL(wp) :: c1 !<
	331	REAL(wp) :: c2 !<
	332	REAL(wp) :: h1 !<
	333	REAL(wp) :: h2 !<
	334	REAL(wp) :: zp !<
	335	REAL(wp) :: pts !<
	336	REAL(wp) :: pt_i !<
	337	REAL(wp) :: rifs !<
	338	REAL(wp) :: u_i !<
	339	REAL(wp) :: v_i !<
	340	REAL(wp) :: us_wall !<
	341	REAL(wp) :: vel_total !<
	342	REAL(wp) :: ws !<
	343	REAL(wp) :: wspts !<
[1015]	344
[1320]	345	REAL(wp), &
	346	DIMENSION(nysg:nyng,nxlg:nxrg) :: &
[1682]	347	wall !<
[1320]	348
	349	REAL(wp), &
	350	DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) :: &
[1682]	351	wall_flux !<
[1015]	352
	353
[1353]	354	zp = 0.5_wp * ( (a+c1) * dy + (b+c2) * dx )
	355	wall_flux = 0.0_wp
[1015]	356	wall_index = NINT( a+ 2b + 3c1 + 4*c2 )
	357
	358	min_inner = MINVAL( nzb_uvw_inner(nys:nyn,nxl:nxr) ) + 1
	359	max_outer = MINVAL( nzb_uvw_outer(nys:nyn,nxl:nxr) )
	360
	361	!$acc kernels present( hom, nzb_uvw_inner, nzb_uvw_outer, pt, rif_wall ) &
	362	!$acc present( u, v, w, wall, wall_flux, z0 )
[1153]	363	!$acc loop independent
[1128]	364	DO i = i_left, i_right
	365	DO j = j_south, j_north
[1153]	366
[1353]	367	IF ( wall(j,i) /= 0.0_wp ) THEN
[1015]	368	!
	369	!-- All subsequent variables are computed for the respective
	370	!-- location where the respective flux is defined.
[1257]	371	!$acc loop independent
[1153]	372	DO k = nzb_uvw_inner(j,i)+1, nzb_uvw_outer(j,i)
	373
[1015]	374	!
	375	!-- (1) Compute rifs, u_i, v_i, ws, pt' and w'pt'
	376	rifs = rif_wall(k,j,i,wall_index)
	377
[1353]	378	u_i = a * u(k,j,i) + c1 * 0.25_wp * &
[1015]	379	( u(k+1,j,i+1) + u(k+1,j,i) + u(k,j,i+1) + u(k,j,i) )
	380
[1353]	381	v_i = b * v(k,j,i) + c2 * 0.25_wp * &
[1015]	382	( v(k+1,j+1,i) + v(k+1,j,i) + v(k,j+1,i) + v(k,j,i) )
	383
[1353]	384	ws = ( c1 + c2 ) * w(k,j,i) + 0.25_wp * ( &
[1015]	385	a * ( w(k-1,j,i-1) + w(k-1,j,i) + w(k,j,i-1) + w(k,j,i) ) &
	386	+ b * ( w(k-1,j-1,i) + w(k-1,j,i) + w(k,j-1,i) + w(k,j,i) ) &
[1353]	387	)
	388	pt_i = 0.5_wp * ( pt(k,j,i) + a * pt(k,j,i-1) + &
[1015]	389	b * pt(k,j-1,i) + ( c1 + c2 ) * pt(k+1,j,i) )
	390
	391	pts = pt_i - hom(k,1,4,0)
	392	wspts = ws * pts
	393
	394	!
	395	!-- (2) Compute wall-parallel absolute velocity vel_total
	396	vel_total = SQRT( ws*2 + (a+c1) u_i*2 + (b+c2) v_i**2 )
	397
	398	!
	399	!-- (3) Compute wall friction velocity us_wall
[1353]	400	IF ( rifs >= 0.0_wp ) THEN
[1015]	401
	402	!
	403	!-- Stable stratification (and neutral)
	404	us_wall = kappa * vel_total / ( LOG( zp / z0(j,i) ) + &
[1353]	405	5.0_wp * rifs * ( zp - z0(j,i) ) / zp &
[1015]	406	)
	407	ELSE
	408
	409	!
	410	!-- Unstable stratification
[1353]	411	h1 = SQRT( SQRT( 1.0_wp - 16.0_wp * rifs ) )
	412	h2 = SQRT( SQRT( 1.0_wp - 16.0_wp * rifs * z0(j,i) / zp ) )
[1015]	413
	414	us_wall = kappa * vel_total / ( &
	415	LOG( zp / z0(j,i) ) - &
[1353]	416	LOG( ( 1.0_wp + h1 )*2 ( 1.0_wp + h1**2 ) / ( &
	417	( 1.0_wp + h2 )*2 ( 1.0_wp + h2**2 ) ) ) +&
	418	2.0_wp * ( ATAN( h1 ) - ATAN( h2 ) ) &
[1015]	419	)
	420	ENDIF
	421
	422	!
	423	!-- (4) Compute zp/L (corresponds to neutral Richardson flux
	424	!-- number rifs)
[1353]	425	rifs = -1.0_wp * zp * kappa * g * wspts / &
	426	( pt_i * ( us_wall**3 + 1E-30 ) )
[1015]	427
	428	!
	429	!-- Limit the value range of the Richardson numbers.
	430	!-- This is necessary for very small velocities (u,w --> 0),
	431	!-- because the absolute value of rif can then become very
	432	!-- large, which in consequence would result in very large
	433	!-- shear stresses and very small momentum fluxes (both are
	434	!-- generally unrealistic).
[1691]	435	IF ( rifs < zeta_min ) rifs = zeta_min
	436	IF ( rifs > zeta_max ) rifs = zeta_max
[1015]	437
	438	!
	439	!-- (5) Compute wall_flux (u'v', v'u', w'v', or w'u')
[1353]	440	IF ( rifs >= 0.0_wp ) THEN
[1015]	441
	442	!
	443	!-- Stable stratification (and neutral)
	444	wall_flux(k,j,i) = kappa * &
	445	( au(k,j,i) + bv(k,j,i) + (c1+c2)*w(k,j,i) ) / &
	446	( LOG( zp / z0(j,i) ) + &
[1353]	447	5.0_wp * rifs * ( zp - z0(j,i) ) / zp &
[1015]	448	)
	449	ELSE
	450
	451	!
	452	!-- Unstable stratification
[1353]	453	h1 = SQRT( SQRT( 1.0_wp - 16.0_wp * rifs ) )
	454	h2 = SQRT( SQRT( 1.0_wp - 16.0_wp * rifs * z0(j,i) / zp ) )
[1015]	455
	456	wall_flux(k,j,i) = kappa * &
	457	( au(k,j,i) + bv(k,j,i) + (c1+c2)*w(k,j,i) ) / ( &
	458	LOG( zp / z0(j,i) ) - &
[1353]	459	LOG( ( 1.0_wp + h1 )*2 ( 1.0_wp + h1**2 ) / ( &
	460	( 1.0_wp + h2 )*2 ( 1.0_wp + h2**2 ) ) ) +&
	461	2.0_wp * ( ATAN( h1 ) - ATAN( h2 ) ) &
[1015]	462	)
	463	ENDIF
	464	wall_flux(k,j,i) = -wall_flux(k,j,i) * us_wall
	465
	466	!
	467	!-- store rifs for next time step
	468	rif_wall(k,j,i,wall_index) = rifs
	469
[1153]	470	ENDDO
	471
	472	ENDIF
	473
[1015]	474	ENDDO
	475	ENDDO
	476	!$acc end kernels
	477
	478	END SUBROUTINE wall_fluxes_acc
	479
	480
[56]	481	!------------------------------------------------------------------------------!
[1682]	482	! Description:
	483	! ------------
	484	!> Call for all grid point i,j
[56]	485	!------------------------------------------------------------------------------!
	486	SUBROUTINE wall_fluxes_ij( i, j, nzb_w, nzt_w, wall_flux, a, b, c1, c2 )
	487
[1320]	488	USE arrays_3d, &
	489	ONLY: rif_wall, pt, u, v, w, z0
	490
	491	USE control_parameters, &
[1691]	492	ONLY: g, kappa, zeta_max, zeta_min
[1320]	493
	494	USE grid_variables, &
	495	ONLY: dx, dy
	496
	497	USE indices, &
	498	ONLY: nzb, nzt
	499
	500	USE kinds
	501
	502	USE statistics, &
	503	ONLY: hom
[56]	504
	505	IMPLICIT NONE
	506
[1682]	507	INTEGER(iwp) :: i !<
	508	INTEGER(iwp) :: j !<
	509	INTEGER(iwp) :: k !<
	510	INTEGER(iwp) :: nzb_w !<
	511	INTEGER(iwp) :: nzt_w !<
	512	INTEGER(iwp) :: wall_index !<
[1320]	513
[1682]	514	REAL(wp) :: a !<
	515	REAL(wp) :: b !<
	516	REAL(wp) :: c1 !<
	517	REAL(wp) :: c2 !<
	518	REAL(wp) :: h1 !<
	519	REAL(wp) :: h2 !<
	520	REAL(wp) :: zp !<
	521	REAL(wp) :: pts !<
	522	REAL(wp) :: pt_i !<
	523	REAL(wp) :: rifs !<
	524	REAL(wp) :: u_i !<
	525	REAL(wp) :: v_i !<
	526	REAL(wp) :: us_wall !<
	527	REAL(wp) :: vel_total !<
	528	REAL(wp) :: ws !<
	529	REAL(wp) :: wspts !<
[56]	530
[1682]	531	REAL(wp), DIMENSION(nzb:nzt+1) :: wall_flux !<
[56]	532
	533
[1353]	534	zp = 0.5_wp * ( (a+c1) * dy + (b+c2) * dx )
	535	wall_flux = 0.0_wp
[56]	536	wall_index = NINT( a+ 2b + 3c1 + 4*c2 )
	537
	538	!
	539	!-- All subsequent variables are computed for the respective location where
[187]	540	!-- the respective flux is defined.
[56]	541	DO k = nzb_w, nzt_w
	542
	543	!
	544	!-- (1) Compute rifs, u_i, v_i, ws, pt' and w'pt'
	545	rifs = rif_wall(k,j,i,wall_index)
	546
[1353]	547	u_i = a * u(k,j,i) + c1 * 0.25_wp * &
[56]	548	( u(k+1,j,i+1) + u(k+1,j,i) + u(k,j,i+1) + u(k,j,i) )
	549
[1353]	550	v_i = b * v(k,j,i) + c2 * 0.25_wp * &
[56]	551	( v(k+1,j+1,i) + v(k+1,j,i) + v(k,j+1,i) + v(k,j,i) )
	552
[1353]	553	ws = ( c1 + c2 ) * w(k,j,i) + 0.25_wp * ( &
[56]	554	a * ( w(k-1,j,i-1) + w(k-1,j,i) + w(k,j,i-1) + w(k,j,i) ) &
	555	+ b * ( w(k-1,j-1,i) + w(k-1,j,i) + w(k,j-1,i) + w(k,j,i) ) &
[1353]	556	)
	557	pt_i = 0.5_wp * ( pt(k,j,i) + a * pt(k,j,i-1) + b * pt(k,j-1,i) &
[56]	558	+ ( c1 + c2 ) * pt(k+1,j,i) )
	559
	560	pts = pt_i - hom(k,1,4,0)
	561	wspts = ws * pts
	562
	563	!
	564	!-- (2) Compute wall-parallel absolute velocity vel_total
	565	vel_total = SQRT( ws*2 + ( a+c1 ) u_i*2 + ( b+c2 ) v_i**2 )
	566
	567	!
	568	!-- (3) Compute wall friction velocity us_wall
[1353]	569	IF ( rifs >= 0.0_wp ) THEN
[56]	570
	571	!
	572	!-- Stable stratification (and neutral)
	573	us_wall = kappa * vel_total / ( LOG( zp / z0(j,i) ) + &
[1353]	574	5.0_wp * rifs * ( zp - z0(j,i) ) / zp &
[56]	575	)
	576	ELSE
	577
	578	!
	579	!-- Unstable stratification
[1353]	580	h1 = SQRT( SQRT( 1.0_wp - 16.0_wp * rifs ) )
	581	h2 = SQRT( SQRT( 1.0_wp - 16.0_wp * rifs * z0(j,i) / zp ) )
[56]	582
[1320]	583	us_wall = kappa * vel_total / ( &
	584	LOG( zp / z0(j,i) ) - &
[1353]	585	LOG( ( 1.0_wp + h1 )*2 ( 1.0_wp + h1**2 ) / ( &
	586	( 1.0_wp + h2 )*2 ( 1.0_wp + h2**2 ) ) ) + &
	587	2.0_wp * ( ATAN( h1 ) - ATAN( h2 ) ) &
[187]	588	)
[56]	589	ENDIF
	590
	591	!
	592	!-- (4) Compute zp/L (corresponds to neutral Richardson flux number
	593	!-- rifs)
[1353]	594	rifs = -1.0_wp * zp * kappa * g * wspts / &
	595	( pt_i * (us_wall**3 + 1E-30) )
[56]	596
	597	!
	598	!-- Limit the value range of the Richardson numbers.
	599	!-- This is necessary for very small velocities (u,w --> 0), because
	600	!-- the absolute value of rif can then become very large, which in
	601	!-- consequence would result in very large shear stresses and very
	602	!-- small momentum fluxes (both are generally unrealistic).
[1691]	603	IF ( rifs < zeta_min ) rifs = zeta_min
	604	IF ( rifs > zeta_max ) rifs = zeta_max
[56]	605
	606	!
	607	!-- (5) Compute wall_flux (u'v', v'u', w'v', or w'u')
[1353]	608	IF ( rifs >= 0.0_wp ) THEN
[56]	609
	610	!
	611	!-- Stable stratification (and neutral)
[1320]	612	wall_flux(k) = kappa * &
	613	( au(k,j,i) + bv(k,j,i) + (c1+c2)*w(k,j,i) ) / &
	614	( LOG( zp / z0(j,i) ) + &
[1353]	615	5.0_wp * rifs * ( zp - z0(j,i) ) / zp &
[53]	616	)
[52]	617	ELSE
[53]	618
[56]	619	!
	620	!-- Unstable stratification
[1353]	621	h1 = SQRT( SQRT( 1.0_wp - 16.0_wp * rifs ) )
	622	h2 = SQRT( SQRT( 1.0_wp - 16.0_wp * rifs * z0(j,i) / zp ) )
[52]	623
[1320]	624	wall_flux(k) = kappa * &
	625	( au(k,j,i) + bv(k,j,i) + (c1+c2)*w(k,j,i) ) / ( &
	626	LOG( zp / z0(j,i) ) - &
[1353]	627	LOG( ( 1.0_wp + h1 )*2 ( 1.0_wp + h1**2 ) / ( &
	628	( 1.0_wp + h2 )*2 ( 1.0_wp + h2**2 ) ) ) + &
	629	2.0_wp * ( ATAN( h1 ) - ATAN( h2 ) ) &
[187]	630	)
[56]	631	ENDIF
[187]	632	wall_flux(k) = -wall_flux(k) * us_wall
[53]	633
[56]	634	!
	635	!-- store rifs for next time step
	636	rif_wall(k,j,i,wall_index) = rifs
[53]	637
[56]	638	ENDDO
[53]	639
[56]	640	END SUBROUTINE wall_fluxes_ij
[53]	641
[56]	642
	643
[53]	644	!------------------------------------------------------------------------------!
	645	! Description:
	646	! ------------
[1682]	647	!> Call for all grid points
	648	!> Calculates momentum fluxes at vertical walls for routine production_e
	649	!> assuming Monin-Obukhov similarity.
	650	!> Indices: usvs a=1, vsus b=1, wsvs c1=1, wsus c2=1 (other=0).
[53]	651	!------------------------------------------------------------------------------!
	652
[1682]	653	SUBROUTINE wall_fluxes_e( wall_flux, a, b, c1, c2, wall )
	654
	655
[1320]	656	USE arrays_3d, &
	657	ONLY: rif_wall, u, v, w, z0
	658
	659	USE control_parameters, &
	660	ONLY: kappa
	661
	662	USE grid_variables, &
	663	ONLY: dx, dy
	664
	665	USE indices, &
	666	ONLY: nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb, &
	667	nzb_diff_s_inner, nzb_diff_s_outer, nzt
	668
	669	USE kinds
[53]	670
[56]	671	IMPLICIT NONE
[53]	672
[1682]	673	INTEGER(iwp) :: i !<
	674	INTEGER(iwp) :: j !<
	675	INTEGER(iwp) :: k !<
	676	INTEGER(iwp) :: kk !<
	677	INTEGER(iwp) :: wall_index !<
[1320]	678
[1682]	679	REAL(wp) :: a !<
	680	REAL(wp) :: b !<
	681	REAL(wp) :: c1 !<
	682	REAL(wp) :: c2 !<
	683	REAL(wp) :: h1 !<
	684	REAL(wp) :: h2 !<
	685	REAL(wp) :: u_i !<
	686	REAL(wp) :: v_i !<
	687	REAL(wp) :: us_wall !<
	688	REAL(wp) :: vel_total !<
	689	REAL(wp) :: vel_zp !<
	690	REAL(wp) :: ws !<
	691	REAL(wp) :: zp !<
	692	REAL(wp) :: rifs !<
[53]	693
[1320]	694	REAL(wp), &
	695	DIMENSION(nysg:nyng,nxlg:nxrg) :: &
[1682]	696	wall !<
[1320]	697
	698	REAL(wp), &
	699	DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) :: &
[1682]	700	wall_flux !<
[53]	701
	702
[1353]	703	zp = 0.5_wp * ( (a+c1) * dy + (b+c2) * dx )
	704	wall_flux = 0.0_wp
[56]	705	wall_index = NINT( a+ 2b + 3c1 + 4*c2 )
[53]	706
[56]	707	DO i = nxl, nxr
	708	DO j = nys, nyn
	709
[1353]	710	IF ( wall(j,i) /= 0.0_wp ) THEN
[53]	711	!
[187]	712	!-- All subsequent variables are computed for scalar locations.
[56]	713	DO k = nzb_diff_s_inner(j,i)-1, nzb_diff_s_outer(j,i)-2
[53]	714	!
[187]	715	!-- (1) Compute rifs, u_i, v_i, and ws
[56]	716	IF ( k == nzb_diff_s_inner(j,i)-1 ) THEN
	717	kk = nzb_diff_s_inner(j,i)-1
	718	ELSE
	719	kk = k-1
	720	ENDIF
[1353]	721	rifs = 0.5_wp * ( rif_wall(k,j,i,wall_index) + &
	722	a * rif_wall(k,j,i+1,1) + &
	723	b * rif_wall(k,j+1,i,2) + &
	724	c1 * rif_wall(kk,j,i,3) + &
	725	c2 * rif_wall(kk,j,i,4) &
	726	)
[53]	727
[1353]	728	u_i = 0.5_wp * ( u(k,j,i) + u(k,j,i+1) )
	729	v_i = 0.5_wp * ( v(k,j,i) + v(k,j+1,i) )
	730	ws = 0.5_wp * ( w(k,j,i) + w(k-1,j,i) )
[53]	731	!
[187]	732	!-- (2) Compute wall-parallel absolute velocity vel_total and
	733	!-- interpolate appropriate velocity component vel_zp.
	734	vel_total = SQRT( ws*2 + (a+c1) u_i*2 + (b+c2) v_i**2 )
[1353]	735	vel_zp = 0.5_wp * ( a * u_i + b * v_i + (c1+c2) * ws )
[187]	736	!
	737	!-- (3) Compute wall friction velocity us_wall
[1353]	738	IF ( rifs >= 0.0_wp ) THEN
[53]	739
	740	!
[187]	741	!-- Stable stratification (and neutral)
	742	us_wall = kappa * vel_total / ( LOG( zp / z0(j,i) ) + &
[1353]	743	5.0_wp * rifs * ( zp - z0(j,i) ) / zp &
[187]	744	)
	745	ELSE
	746
	747	!
	748	!-- Unstable stratification
[1353]	749	h1 = SQRT( SQRT( 1.0_wp - 16.0_wp * rifs ) )
	750	h2 = SQRT( SQRT( 1.0_wp - 16.0_wp * rifs * z0(j,i) / zp ) )
[187]	751
	752	us_wall = kappa * vel_total / ( &
	753	LOG( zp / z0(j,i) ) - &
[1353]	754	LOG( ( 1.0_wp + h1 )*2 ( 1.0_wp + h1**2 ) / ( &
	755	( 1.0_wp + h2 )*2 ( 1.0_wp + h2**2 ) ) ) +&
	756	2.0_wp * ( ATAN( h1 ) - ATAN( h2 ) ) &
[187]	757	)
	758	ENDIF
	759
	760	!
	761	!-- Skip step (4) of wall_fluxes, because here rifs is already
	762	!-- available from (1)
	763	!
[56]	764	!-- (5) Compute wall_flux (u'v', v'u', w'v', or w'u')
[55]	765
[1353]	766	IF ( rifs >= 0.0_wp ) THEN
[53]	767
	768	!
[56]	769	!-- Stable stratification (and neutral)
[1353]	770	wall_flux(k,j,i) = kappa * vel_zp / ( LOG( zp/z0(j,i) ) +&
	771	5.0_wp * rifs * ( zp-z0(j,i) ) / zp )
[56]	772	ELSE
[53]	773
	774	!
[56]	775	!-- Unstable stratification
[1353]	776	h1 = SQRT( SQRT( 1.0_wp - 16.0_wp * rifs ) )
	777	h2 = SQRT( SQRT( 1.0_wp - 16.0_wp * rifs * z0(j,i) / zp ) )
[53]	778
[187]	779	wall_flux(k,j,i) = kappa * vel_zp / ( &
	780	LOG( zp / z0(j,i) ) - &
[1353]	781	LOG( ( 1.0_wp + h1 )*2 ( 1.0_wp + h1**2 ) / ( &
	782	( 1.0_wp + h2 )*2 ( 1.0_wp + h2**2 ) ) ) +&
	783	2.0_wp * ( ATAN( h1 ) - ATAN( h2 ) ) &
[187]	784	)
[56]	785	ENDIF
[187]	786	wall_flux(k,j,i) = - wall_flux(k,j,i) * us_wall
[56]	787
	788	ENDDO
	789
	790	ENDIF
	791
	792	ENDDO
	793	ENDDO
	794
	795	END SUBROUTINE wall_fluxes_e
	796
	797
[1015]	798	!------------------------------------------------------------------------------!
	799	! Description:
	800	! ------------
[1682]	801	!> Call for all grid points - accelerator version
	802	!> Calculates momentum fluxes at vertical walls for routine production_e
	803	!> assuming Monin-Obukhov similarity.
	804	!> Indices: usvs a=1, vsus b=1, wsvs c1=1, wsus c2=1 (other=0).
[1015]	805	!------------------------------------------------------------------------------!
[1682]	806	SUBROUTINE wall_fluxes_e_acc( wall_flux, a, b, c1, c2, wall )
[1015]	807
[1682]	808
[1320]	809	USE arrays_3d, &
	810	ONLY: rif_wall, u, v, w, z0
	811
	812	USE control_parameters, &
	813	ONLY: kappa
	814
	815	USE grid_variables, &
	816	ONLY: dx, dy
	817
	818	USE indices, &
	819	ONLY: i_left, i_right, j_north, j_south, nxl, nxlg, nxr, nxrg, &
	820	nyn, nyng, nys, nysg, nzb, nzb_diff_s_inner, &
	821	nzb_diff_s_outer, nzt
	822
	823	USE kinds
[1015]	824
	825	IMPLICIT NONE
	826
[1682]	827	INTEGER(iwp) :: i !<
	828	INTEGER(iwp) :: j !<
	829	INTEGER(iwp) :: k !<
	830	INTEGER(iwp) :: kk !<
	831	INTEGER(iwp) :: max_outer !<
	832	INTEGER(iwp) :: min_inner !<
	833	INTEGER(iwp) :: wall_index !<
[1320]	834
[1682]	835	REAL(wp) :: a !<
	836	REAL(wp) :: b !<
	837	REAL(wp) :: c1 !<
	838	REAL(wp) :: c2 !<
	839	REAL(wp) :: h1 !<
	840	REAL(wp) :: h2 !<
	841	REAL(wp) :: u_i !<
	842	REAL(wp) :: v_i !<
	843	REAL(wp) :: us_wall !<
	844	REAL(wp) :: vel_total !<
	845	REAL(wp) :: vel_zp !<
	846	REAL(wp) :: ws !<
	847	REAL(wp) :: zp !<
	848	REAL(wp) :: rifs !<
[1015]	849
[1320]	850	REAL(wp), &
	851	DIMENSION(nysg:nyng,nxlg:nxrg) :: &
[1682]	852	wall !<
[1320]	853
	854	REAL(wp), &
	855	DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) :: &
[1682]	856	wall_flux !<
[1015]	857
	858
[1353]	859	zp = 0.5_wp * ( (a+c1) * dy + (b+c2) * dx )
	860	wall_flux = 0.0_wp
[1015]	861	wall_index = NINT( a+ 2b + 3c1 + 4*c2 )
	862
	863	min_inner = MINVAL( nzb_diff_s_inner(nys:nyn,nxl:nxr) ) - 1
	864	max_outer = MAXVAL( nzb_diff_s_outer(nys:nyn,nxl:nxr) ) - 2
	865
[1374]	866	!$acc kernels present( nzb_diff_s_inner, nzb_diff_s_outer, rif_wall ) &
[1015]	867	!$acc present( u, v, w, wall, wall_flux, z0 )
[1128]	868	DO i = i_left, i_right
	869	DO j = j_south, j_north
[1015]	870	DO k = min_inner, max_outer
	871	!
	872	!-- All subsequent variables are computed for scalar locations
[1320]	873	IF ( k >= nzb_diff_s_inner(j,i)-1 .AND. &
[1353]	874	k <= nzb_diff_s_outer(j,i)-2 .AND. &
	875	wall(j,i) /= 0.0_wp ) THEN
[1015]	876	!
	877	!-- (1) Compute rifs, u_i, v_i, and ws
	878	IF ( k == nzb_diff_s_inner(j,i)-1 ) THEN
	879	kk = nzb_diff_s_inner(j,i)-1
	880	ELSE
	881	kk = k-1
	882	ENDIF
[1353]	883	rifs = 0.5_wp * ( rif_wall(k,j,i,wall_index) + &
	884	a * rif_wall(k,j,i+1,1) + &
	885	b * rif_wall(k,j+1,i,2) + &
	886	c1 * rif_wall(kk,j,i,3) + &
	887	c2 * rif_wall(kk,j,i,4) &
	888	)
[1015]	889
[1353]	890	u_i = 0.5_wp * ( u(k,j,i) + u(k,j,i+1) )
	891	v_i = 0.5_wp * ( v(k,j,i) + v(k,j+1,i) )
	892	ws = 0.5_wp * ( w(k,j,i) + w(k-1,j,i) )
[1015]	893	!
	894	!-- (2) Compute wall-parallel absolute velocity vel_total and
	895	!-- interpolate appropriate velocity component vel_zp.
	896	vel_total = SQRT( ws*2 + (a+c1) u_i*2 + (b+c2) v_i**2 )
[1353]	897	vel_zp = 0.5_wp * ( a * u_i + b * v_i + (c1+c2) * ws )
[1015]	898	!
	899	!-- (3) Compute wall friction velocity us_wall
[1353]	900	IF ( rifs >= 0.0_wp ) THEN
[1015]	901
	902	!
	903	!-- Stable stratification (and neutral)
	904	us_wall = kappa * vel_total / ( LOG( zp / z0(j,i) ) + &
[1353]	905	5.0_wp * rifs * ( zp - z0(j,i) ) / zp &
[1015]	906	)
	907	ELSE
	908
	909	!
	910	!-- Unstable stratification
[1353]	911	h1 = SQRT( SQRT( 1.0_wp - 16.0_wp * rifs ) )
	912	h2 = SQRT( SQRT( 1.0_wp - 16.0_wp * rifs * z0(j,i) / zp ) )
[1015]	913
	914	us_wall = kappa * vel_total / ( &
	915	LOG( zp / z0(j,i) ) - &
[1353]	916	LOG( ( 1.0_wp + h1 )*2 ( 1.0_wp + h1**2 ) / ( &
	917	( 1.0_wp + h2 )*2 ( 1.0_wp + h2**2 ) ) ) +&
	918	2.0_wp * ( ATAN( h1 ) - ATAN( h2 ) ) &
[1015]	919	)
	920	ENDIF
	921
	922	!
	923	!-- Skip step (4) of wall_fluxes, because here rifs is already
	924	!-- available from (1)
	925	!
	926	!-- (5) Compute wall_flux (u'v', v'u', w'v', or w'u')
	927
[1353]	928	IF ( rifs >= 0.0_wp ) THEN
[1015]	929
	930	!
	931	!-- Stable stratification (and neutral)
[1353]	932	wall_flux(k,j,i) = kappa * vel_zp / ( &
	933	LOG( zp/z0(j,i) ) + &
	934	5.0_wp * rifs * ( zp-z0(j,i) ) / zp &
	935	)
[1015]	936	ELSE
	937
	938	!
	939	!-- Unstable stratification
[1353]	940	h1 = SQRT( SQRT( 1.0_wp - 16.0_wp * rifs ) )
	941	h2 = SQRT( SQRT( 1.0_wp - 16.0_wp * rifs * z0(j,i) / zp ) )
[1015]	942
	943	wall_flux(k,j,i) = kappa * vel_zp / ( &
	944	LOG( zp / z0(j,i) ) - &
[1353]	945	LOG( ( 1.0_wp + h1 )*2 ( 1.0_wp + h1**2 ) / ( &
	946	( 1.0_wp + h2 )*2 ( 1.0_wp + h2**2 ) ) ) +&
	947	2.0_wp * ( ATAN( h1 ) - ATAN( h2 ) ) &
[1015]	948	)
	949	ENDIF
	950	wall_flux(k,j,i) = - wall_flux(k,j,i) * us_wall
	951
	952	ENDIF
	953
	954	ENDDO
	955	ENDDO
	956	ENDDO
	957	!$acc end kernels
	958
	959	END SUBROUTINE wall_fluxes_e_acc
	960
	961
	962	!------------------------------------------------------------------------------!
[1682]	963	! Description:
	964	! ------------
	965	!> Call for grid point i,j
[56]	966	!------------------------------------------------------------------------------!
	967	SUBROUTINE wall_fluxes_e_ij( i, j, nzb_w, nzt_w, wall_flux, a, b, c1, c2 )
	968
[1320]	969	USE arrays_3d, &
	970	ONLY: rif_wall, u, v, w, z0
	971
	972	USE control_parameters, &
	973	ONLY: kappa
	974
	975	USE grid_variables, &
	976	ONLY: dx, dy
	977
	978	USE indices, &
	979	ONLY: nzb, nzt
	980
	981	USE kinds
[56]	982
	983	IMPLICIT NONE
	984
[1682]	985	INTEGER(iwp) :: i !<
	986	INTEGER(iwp) :: j !<
	987	INTEGER(iwp) :: k !<
	988	INTEGER(iwp) :: kk !<
	989	INTEGER(iwp) :: nzb_w !<
	990	INTEGER(iwp) :: nzt_w !<
	991	INTEGER(iwp) :: wall_index !<
[1320]	992
[1682]	993	REAL(wp) :: a !<
	994	REAL(wp) :: b !<
	995	REAL(wp) :: c1 !<
	996	REAL(wp) :: c2 !<
	997	REAL(wp) :: h1 !<
	998	REAL(wp) :: h2 !<
	999	REAL(wp) :: u_i !<
	1000	REAL(wp) :: v_i !<
	1001	REAL(wp) :: us_wall !<
	1002	REAL(wp) :: vel_total !<
	1003	REAL(wp) :: vel_zp !<
	1004	REAL(wp) :: ws !<
	1005	REAL(wp) :: zp !<
	1006	REAL(wp) :: rifs !<
[56]	1007
[1682]	1008	REAL(wp), DIMENSION(nzb:nzt+1) :: wall_flux !<
[56]	1009
	1010
[1353]	1011	zp = 0.5_wp * ( (a+c1) * dy + (b+c2) * dx )
	1012	wall_flux = 0.0_wp
[56]	1013	wall_index = NINT( a+ 2b + 3c1 + 4*c2 )
	1014
	1015	!
[187]	1016	!-- All subsequent variables are computed for scalar locations.
[56]	1017	DO k = nzb_w, nzt_w
	1018
	1019	!
[187]	1020	!-- (1) Compute rifs, u_i, v_i, and ws
[56]	1021	IF ( k == nzb_w ) THEN
	1022	kk = nzb_w
[53]	1023	ELSE
[56]	1024	kk = k-1
	1025	ENDIF
[1353]	1026	rifs = 0.5_wp * ( rif_wall(k,j,i,wall_index) + &
	1027	a * rif_wall(k,j,i+1,1) + &
	1028	b * rif_wall(k,j+1,i,2) + &
	1029	c1 * rif_wall(kk,j,i,3) + &
	1030	c2 * rif_wall(kk,j,i,4) &
	1031	)
[56]	1032
[1353]	1033	u_i = 0.5_wp * ( u(k,j,i) + u(k,j,i+1) )
	1034	v_i = 0.5_wp * ( v(k,j,i) + v(k,j+1,i) )
	1035	ws = 0.5_wp * ( w(k,j,i) + w(k-1,j,i) )
[56]	1036	!
[187]	1037	!-- (2) Compute wall-parallel absolute velocity vel_total and
	1038	!-- interpolate appropriate velocity component vel_zp.
	1039	vel_total = SQRT( ws*2 + (a+c1) u_i*2 + (b+c2) v_i**2 )
[1353]	1040	vel_zp = 0.5_wp * ( a * u_i + b * v_i + (c1+c2) * ws )
[187]	1041	!
	1042	!-- (3) Compute wall friction velocity us_wall
[1353]	1043	IF ( rifs >= 0.0_wp ) THEN
[56]	1044
	1045	!
[187]	1046	!-- Stable stratification (and neutral)
	1047	us_wall = kappa * vel_total / ( LOG( zp / z0(j,i) ) + &
[1353]	1048	5.0_wp * rifs * ( zp - z0(j,i) ) / zp &
[187]	1049	)
	1050	ELSE
	1051
	1052	!
	1053	!-- Unstable stratification
[1353]	1054	h1 = SQRT( SQRT( 1.0_wp - 16.0_wp * rifs ) )
	1055	h2 = SQRT( SQRT( 1.0_wp - 16.0_wp * rifs * z0(j,i) / zp ) )
[187]	1056
[1320]	1057	us_wall = kappa * vel_total / ( &
	1058	LOG( zp / z0(j,i) ) - &
[1353]	1059	LOG( ( 1.0_wp + h1 )*2 ( 1.0_wp + h1**2 ) / ( &
	1060	( 1.0_wp + h2 )*2 ( 1.0_wp + h2**2 ) ) ) + &
	1061	2.0_wp * ( ATAN( h1 ) - ATAN( h2 ) ) &
[187]	1062	)
	1063	ENDIF
	1064
	1065	!
	1066	!-- Skip step (4) of wall_fluxes, because here rifs is already
	1067	!-- available from (1)
	1068	!
[56]	1069	!-- (5) Compute wall_flux (u'v', v'u', w'v', or w'u')
[187]	1070	!-- First interpolate the velocity (this is different from
	1071	!-- subroutine wall_fluxes because fluxes in subroutine
	1072	!-- wall_fluxes_e are defined at scalar locations).
[1353]	1073	vel_zp = 0.5_wp * ( a * ( u(k,j,i) + u(k,j,i+1) ) + &
	1074	b * ( v(k,j,i) + v(k,j+1,i) ) + &
	1075	(c1+c2) * ( w(k,j,i) + w(k-1,j,i) ) &
	1076	)
[56]	1077
[1353]	1078	IF ( rifs >= 0.0_wp ) THEN
[56]	1079
	1080	!
	1081	!-- Stable stratification (and neutral)
[1320]	1082	wall_flux(k) = kappa * vel_zp / &
[1353]	1083	( LOG( zp/z0(j,i) ) + 5.0_wp * rifs * ( zp-z0(j,i) ) / zp )
[56]	1084	ELSE
	1085
	1086	!
	1087	!-- Unstable stratification
[1353]	1088	h1 = SQRT( SQRT( 1.0_wp - 16.0_wp * rifs ) )
	1089	h2 = SQRT( SQRT( 1.0_wp - 16.0_wp * rifs * z0(j,i) / zp ) )
[56]	1090
[1320]	1091	wall_flux(k) = kappa * vel_zp / ( &
	1092	LOG( zp / z0(j,i) ) - &
[1353]	1093	LOG( ( 1.0_wp + h1 )*2 ( 1.0_wp + h1**2 ) / ( &
	1094	( 1.0_wp + h2 )*2 ( 1.0_wp + h2**2 ) ) ) + &
	1095	2.0_wp * ( ATAN( h1 ) - ATAN( h2 ) ) &
	1096	)
[53]	1097	ENDIF
[187]	1098	wall_flux(k) = - wall_flux(k) * us_wall
[53]	1099
[56]	1100	ENDDO
[53]	1101
[56]	1102	END SUBROUTINE wall_fluxes_e_ij
	1103
	1104	END MODULE wall_fluxes_mod

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |