Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |