Home

Context Navigation

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

Last change on this file since 1361 was 1354, checked in by heinze, 10 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 43.1 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |