Home

Context Navigation

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

Last change on this file since 1876 was 1874, checked in by maronga, 9 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 43.4 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |