Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |