Home

Context Navigation

source: palm/trunk/SOURCE/wall_fluxes_mod.f90 @ 1850

Last change on this file since 1850 was 1850, checked in by maronga, 8 years ago
added _mod string to several filenames to meet the naming convection for modules
Property svn:keywords set to `Id`
File size: 43.3 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |