Home

Context Navigation

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

Last change on this file since 2076 was 2001, checked in by knoop, 8 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 43.5 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |