Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |