Home

Context Navigation

wall_fluxes_mod.f90 @ 1859

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

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

Context Navigation

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

Download in other formats: