Home

Context Navigation

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

Last change on this file since 1350 was 1321, checked in by raasch, 11 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 42.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	! -----------------
[1321]	22	!
	23	!
	24	! Former revisions:
	25	! -----------------
	26	! $Id: wall_fluxes.f90 1321 2014-03-20 09:40:40Z maronga $
	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
[56]	150	zp = 0.5 * ( (a+c1) * dy + (b+c2) * dx )
	151	wall_flux = 0.0
	152	wall_index = NINT( a+ 2b + 3c1 + 4*c2 )
	153
[75]	154	DO i = nxl, nxr
	155	DO j = nys, nyn
[56]	156
	157	IF ( wall(j,i) /= 0.0 ) 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
[1320]	167	u_i = a * u(k,j,i) + c1 * 0.25 * &
[56]	168	( u(k+1,j,i+1) + u(k+1,j,i) + u(k,j,i+1) + u(k,j,i) )
[53]	169
[1320]	170	v_i = b * v(k,j,i) + c2 * 0.25 * &
[56]	171	( v(k+1,j+1,i) + v(k+1,j,i) + v(k,j+1,i) + v(k,j,i) )
[53]	172
[56]	173	ws = ( c1 + c2 ) * w(k,j,i) + 0.25 * ( &
	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) ) &
	176	)
[1320]	177	pt_i = 0.5 * ( 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
	189	IF ( rifs >= 0.0 ) THEN
[53]	190
[52]	191	!
[56]	192	!-- Stable stratification (and neutral)
	193	us_wall = kappa * vel_total / ( LOG( zp / z0(j,i) ) + &
	194	5.0 * rifs * ( zp - z0(j,i) ) / zp &
	195	)
	196	ELSE
[53]	197
[52]	198	!
[56]	199	!-- Unstable stratification
[187]	200	h1 = SQRT( SQRT( 1.0 - 16.0 * rifs ) )
	201	h2 = SQRT( SQRT( 1.0 - 16.0 * rifs * z0(j,i) / zp ) )
[53]	202
[187]	203	us_wall = kappa * vel_total / ( &
	204	LOG( zp / z0(j,i) ) - &
	205	LOG( ( 1.0 + h1 )*2 ( 1.0 + h1**2 ) / ( &
	206	( 1.0 + h2 )*2 ( 1.0 + h2**2 ) ) ) + &
	207	2.0 * ( ATAN( h1 ) - ATAN( h2 ) ) &
	208	)
[56]	209	ENDIF
[53]	210
[52]	211	!
[56]	212	!-- (4) Compute zp/L (corresponds to neutral Richardson flux
	213	!-- number rifs)
[1320]	214	rifs = -1.0 * zp * kappa * g * wspts / ( pt_i * &
[56]	215	( 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')
	229	IF ( rifs >= 0.0 ) 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) ) + &
	236	5.0 * rifs * ( zp - z0(j,i) ) / zp &
	237	)
	238	ELSE
[53]	239
[52]	240	!
[56]	241	!-- Unstable stratification
[187]	242	h1 = SQRT( SQRT( 1.0 - 16.0 * rifs ) )
	243	h2 = SQRT( SQRT( 1.0 - 16.0 * 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) ) - &
	248	LOG( ( 1.0 + h1 )*2 ( 1.0 + h1**2 ) / ( &
	249	( 1.0 + h2 )*2 ( 1.0 + h2**2 ) ) ) + &
	250	2.0 * ( ATAN( h1 ) - ATAN( h2 ) ) &
	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
	335	zp = 0.5 * ( (a+c1) * dy + (b+c2) * dx )
	336	wall_flux = 0.0
	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
	348	IF ( wall(j,i) /= 0.0 ) 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
[1320]	359	u_i = a * u(k,j,i) + c1 * 0.25 * &
[1015]	360	( u(k+1,j,i+1) + u(k+1,j,i) + u(k,j,i+1) + u(k,j,i) )
	361
[1320]	362	v_i = b * v(k,j,i) + c2 * 0.25 * &
[1015]	363	( v(k+1,j+1,i) + v(k+1,j,i) + v(k,j+1,i) + v(k,j,i) )
	364
	365	ws = ( c1 + c2 ) * w(k,j,i) + 0.25 * ( &
	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) ) &
	368	)
[1320]	369	pt_i = 0.5 * ( 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
	381	IF ( rifs >= 0.0 ) THEN
	382
	383	!
	384	!-- Stable stratification (and neutral)
	385	us_wall = kappa * vel_total / ( LOG( zp / z0(j,i) ) + &
	386	5.0 * rifs * ( zp - z0(j,i) ) / zp &
	387	)
	388	ELSE
	389
	390	!
	391	!-- Unstable stratification
	392	h1 = SQRT( SQRT( 1.0 - 16.0 * rifs ) )
	393	h2 = SQRT( SQRT( 1.0 - 16.0 * rifs * z0(j,i) / zp ) )
	394
	395	us_wall = kappa * vel_total / ( &
	396	LOG( zp / z0(j,i) ) - &
	397	LOG( ( 1.0 + h1 )*2 ( 1.0 + h1**2 ) / ( &
	398	( 1.0 + h2 )*2 ( 1.0 + h2**2 ) ) ) + &
	399	2.0 * ( ATAN( h1 ) - ATAN( h2 ) ) &
	400	)
	401	ENDIF
	402
	403	!
	404	!-- (4) Compute zp/L (corresponds to neutral Richardson flux
	405	!-- number rifs)
[1320]	406	rifs = -1.0 * zp * kappa * g * wspts / ( pt_i * &
[1015]	407	( us_wall**3 + 1E-30 ) )
	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')
	421	IF ( rifs >= 0.0 ) THEN
	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) ) + &
	428	5.0 * rifs * ( zp - z0(j,i) ) / zp &
	429	)
	430	ELSE
	431
	432	!
	433	!-- Unstable stratification
	434	h1 = SQRT( SQRT( 1.0 - 16.0 * rifs ) )
	435	h2 = SQRT( SQRT( 1.0 - 16.0 * rifs * z0(j,i) / zp ) )
	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) ) - &
	440	LOG( ( 1.0 + h1 )*2 ( 1.0 + h1**2 ) / ( &
	441	( 1.0 + h2 )*2 ( 1.0 + h2**2 ) ) ) + &
	442	2.0 * ( ATAN( h1 ) - ATAN( h2 ) ) &
	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
	513	zp = 0.5 * ( (a+c1) * dy + (b+c2) * dx )
	514	wall_flux = 0.0
	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
[1320]	526	u_i = a * u(k,j,i) + c1 * 0.25 * &
[56]	527	( u(k+1,j,i+1) + u(k+1,j,i) + u(k,j,i+1) + u(k,j,i) )
	528
[1320]	529	v_i = b * v(k,j,i) + c2 * 0.25 * &
[56]	530	( v(k+1,j+1,i) + v(k+1,j,i) + v(k,j+1,i) + v(k,j,i) )
	531
	532	ws = ( c1 + c2 ) * w(k,j,i) + 0.25 * ( &
	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) ) &
	535	)
[1320]	536	pt_i = 0.5 * ( 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
	548	IF ( rifs >= 0.0 ) THEN
	549
	550	!
	551	!-- Stable stratification (and neutral)
	552	us_wall = kappa * vel_total / ( LOG( zp / z0(j,i) ) + &
	553	5.0 * rifs * ( zp - z0(j,i) ) / zp &
	554	)
	555	ELSE
	556
	557	!
	558	!-- Unstable stratification
[187]	559	h1 = SQRT( SQRT( 1.0 - 16.0 * rifs ) )
	560	h2 = SQRT( SQRT( 1.0 - 16.0 * rifs * z0(j,i) / zp ) )
[56]	561
[1320]	562	us_wall = kappa * vel_total / ( &
	563	LOG( zp / z0(j,i) ) - &
	564	LOG( ( 1.0 + h1 )*2 ( 1.0 + h1**2 ) / ( &
	565	( 1.0 + h2 )*2 ( 1.0 + h2**2 ) ) ) + &
	566	2.0 * ( ATAN( h1 ) - ATAN( h2 ) ) &
[187]	567	)
[56]	568	ENDIF
	569
	570	!
	571	!-- (4) Compute zp/L (corresponds to neutral Richardson flux number
	572	!-- rifs)
	573	rifs = -1.0 * zp * kappa * g * wspts / ( pt_i * (us_wall**3 + 1E-30) )
	574
	575	!
	576	!-- Limit the value range of the Richardson numbers.
	577	!-- This is necessary for very small velocities (u,w --> 0), because
	578	!-- the absolute value of rif can then become very large, which in
	579	!-- consequence would result in very large shear stresses and very
	580	!-- small momentum fluxes (both are generally unrealistic).
	581	IF ( rifs < rif_min ) rifs = rif_min
	582	IF ( rifs > rif_max ) rifs = rif_max
	583
	584	!
	585	!-- (5) Compute wall_flux (u'v', v'u', w'v', or w'u')
	586	IF ( rifs >= 0.0 ) THEN
	587
	588	!
	589	!-- Stable stratification (and neutral)
[1320]	590	wall_flux(k) = kappa * &
	591	( au(k,j,i) + bv(k,j,i) + (c1+c2)*w(k,j,i) ) / &
	592	( LOG( zp / z0(j,i) ) + &
	593	5.0 * rifs * ( zp - z0(j,i) ) / zp &
[53]	594	)
[52]	595	ELSE
[53]	596
[56]	597	!
	598	!-- Unstable stratification
[187]	599	h1 = SQRT( SQRT( 1.0 - 16.0 * rifs ) )
	600	h2 = SQRT( SQRT( 1.0 - 16.0 * rifs * z0(j,i) / zp ) )
[52]	601
[1320]	602	wall_flux(k) = kappa * &
	603	( au(k,j,i) + bv(k,j,i) + (c1+c2)*w(k,j,i) ) / ( &
	604	LOG( zp / z0(j,i) ) - &
	605	LOG( ( 1.0 + h1 )*2 ( 1.0 + h1**2 ) / ( &
	606	( 1.0 + h2 )*2 ( 1.0 + h2**2 ) ) ) + &
	607	2.0 * ( ATAN( h1 ) - ATAN( h2 ) ) &
[187]	608	)
[56]	609	ENDIF
[187]	610	wall_flux(k) = -wall_flux(k) * us_wall
[53]	611
[56]	612	!
	613	!-- store rifs for next time step
	614	rif_wall(k,j,i,wall_index) = rifs
[53]	615
[56]	616	ENDDO
[53]	617
[56]	618	END SUBROUTINE wall_fluxes_ij
[53]	619
[56]	620
	621
[53]	622	!------------------------------------------------------------------------------!
[56]	623	! Call for all grid points
	624	!------------------------------------------------------------------------------!
	625	SUBROUTINE wall_fluxes_e( wall_flux, a, b, c1, c2, wall )
	626
	627	!------------------------------------------------------------------------------!
[53]	628	! Description:
	629	! ------------
	630	! Calculates momentum fluxes at vertical walls for routine production_e
	631	! assuming Monin-Obukhov similarity.
	632	! Indices: usvs a=1, vsus b=1, wsvs c1=1, wsus c2=1 (other=0).
	633	!------------------------------------------------------------------------------!
	634
[1320]	635	USE arrays_3d, &
	636	ONLY: rif_wall, u, v, w, z0
	637
	638	USE control_parameters, &
	639	ONLY: kappa
	640
	641	USE grid_variables, &
	642	ONLY: dx, dy
	643
	644	USE indices, &
	645	ONLY: nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb, &
	646	nzb_diff_s_inner, nzb_diff_s_outer, nzt
	647
	648	USE kinds
[53]	649
[56]	650	IMPLICIT NONE
[53]	651
[1320]	652	INTEGER(iwp) :: i !:
	653	INTEGER(iwp) :: j !:
	654	INTEGER(iwp) :: k !:
	655	INTEGER(iwp) :: kk !:
	656	INTEGER(iwp) :: wall_index !:
	657
	658	REAL(wp) :: a !:
	659	REAL(wp) :: b !:
	660	REAL(wp) :: c1 !:
	661	REAL(wp) :: c2 !:
	662	REAL(wp) :: h1 !:
	663	REAL(wp) :: h2 !:
	664	REAL(wp) :: u_i !:
	665	REAL(wp) :: v_i !:
	666	REAL(wp) :: us_wall !:
	667	REAL(wp) :: vel_total !:
	668	REAL(wp) :: vel_zp !:
	669	REAL(wp) :: ws !:
	670	REAL(wp) :: zp !:
	671	REAL(wp) :: rifs !:
[53]	672
[1320]	673	REAL(wp), &
	674	DIMENSION(nysg:nyng,nxlg:nxrg) :: &
	675	wall !:
	676
	677	REAL(wp), &
	678	DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) :: &
	679	wall_flux !:
[53]	680
	681
[56]	682	zp = 0.5 * ( (a+c1) * dy + (b+c2) * dx )
	683	wall_flux = 0.0
	684	wall_index = NINT( a+ 2b + 3c1 + 4*c2 )
[53]	685
[56]	686	DO i = nxl, nxr
	687	DO j = nys, nyn
	688
	689	IF ( wall(j,i) /= 0.0 ) THEN
[53]	690	!
[187]	691	!-- All subsequent variables are computed for scalar locations.
[56]	692	DO k = nzb_diff_s_inner(j,i)-1, nzb_diff_s_outer(j,i)-2
[53]	693	!
[187]	694	!-- (1) Compute rifs, u_i, v_i, and ws
[56]	695	IF ( k == nzb_diff_s_inner(j,i)-1 ) THEN
	696	kk = nzb_diff_s_inner(j,i)-1
	697	ELSE
	698	kk = k-1
	699	ENDIF
	700	rifs = 0.5 * ( rif_wall(k,j,i,wall_index) + &
	701	a * rif_wall(k,j,i+1,1) + b * rif_wall(k,j+1,i,2) + &
	702	c1 * rif_wall(kk,j,i,3) + c2 * rif_wall(kk,j,i,4) &
	703	)
[53]	704
[187]	705	u_i = 0.5 * ( u(k,j,i) + u(k,j,i+1) )
	706	v_i = 0.5 * ( v(k,j,i) + v(k,j+1,i) )
	707	ws = 0.5 * ( w(k,j,i) + w(k-1,j,i) )
[53]	708	!
[187]	709	!-- (2) Compute wall-parallel absolute velocity vel_total and
	710	!-- interpolate appropriate velocity component vel_zp.
	711	vel_total = SQRT( ws*2 + (a+c1) u_i*2 + (b+c2) v_i**2 )
	712	vel_zp = 0.5 * ( a * u_i + b * v_i + (c1+c2) * ws )
	713	!
	714	!-- (3) Compute wall friction velocity us_wall
	715	IF ( rifs >= 0.0 ) THEN
[53]	716
	717	!
[187]	718	!-- Stable stratification (and neutral)
	719	us_wall = kappa * vel_total / ( LOG( zp / z0(j,i) ) + &
	720	5.0 * rifs * ( zp - z0(j,i) ) / zp &
	721	)
	722	ELSE
	723
	724	!
	725	!-- Unstable stratification
	726	h1 = SQRT( SQRT( 1.0 - 16.0 * rifs ) )
	727	h2 = SQRT( SQRT( 1.0 - 16.0 * rifs * z0(j,i) / zp ) )
	728
	729	us_wall = kappa * vel_total / ( &
	730	LOG( zp / z0(j,i) ) - &
	731	LOG( ( 1.0 + h1 )*2 ( 1.0 + h1**2 ) / ( &
	732	( 1.0 + h2 )*2 ( 1.0 + h2**2 ) ) ) + &
	733	2.0 * ( ATAN( h1 ) - ATAN( h2 ) ) &
	734	)
	735	ENDIF
	736
	737	!
	738	!-- Skip step (4) of wall_fluxes, because here rifs is already
	739	!-- available from (1)
	740	!
[56]	741	!-- (5) Compute wall_flux (u'v', v'u', w'v', or w'u')
[55]	742
[56]	743	IF ( rifs >= 0.0 ) THEN
[53]	744
	745	!
[56]	746	!-- Stable stratification (and neutral)
[1320]	747	wall_flux(k,j,i) = kappa * vel_zp / &
[56]	748	( LOG( zp/z0(j,i) ) + 5.0rifs ( zp-z0(j,i) ) / zp )
	749	ELSE
[53]	750
	751	!
[56]	752	!-- Unstable stratification
[187]	753	h1 = SQRT( SQRT( 1.0 - 16.0 * rifs ) )
	754	h2 = SQRT( SQRT( 1.0 - 16.0 * rifs * z0(j,i) / zp ) )
[53]	755
[187]	756	wall_flux(k,j,i) = kappa * vel_zp / ( &
	757	LOG( zp / z0(j,i) ) - &
	758	LOG( ( 1.0 + h1 )*2 ( 1.0 + h1**2 ) / ( &
	759	( 1.0 + h2 )*2 ( 1.0 + h2**2 ) ) ) + &
	760	2.0 * ( ATAN( h1 ) - ATAN( h2 ) ) &
	761	)
[56]	762	ENDIF
[187]	763	wall_flux(k,j,i) = - wall_flux(k,j,i) * us_wall
[56]	764
	765	ENDDO
	766
	767	ENDIF
	768
	769	ENDDO
	770	ENDDO
	771
	772	END SUBROUTINE wall_fluxes_e
	773
	774
[1015]	775	!------------------------------------------------------------------------------!
	776	! Call for all grid points - accelerator version
	777	!------------------------------------------------------------------------------!
	778	SUBROUTINE wall_fluxes_e_acc( wall_flux, a, b, c1, c2, wall )
[56]	779
	780	!------------------------------------------------------------------------------!
[1015]	781	! Description:
	782	! ------------
	783	! Calculates momentum fluxes at vertical walls for routine production_e
	784	! assuming Monin-Obukhov similarity.
	785	! Indices: usvs a=1, vsus b=1, wsvs c1=1, wsus c2=1 (other=0).
	786	!------------------------------------------------------------------------------!
	787
[1320]	788	USE arrays_3d, &
	789	ONLY: rif_wall, u, v, w, z0
	790
	791	USE control_parameters, &
	792	ONLY: kappa
	793
	794	USE grid_variables, &
	795	ONLY: dx, dy
	796
	797	USE indices, &
	798	ONLY: i_left, i_right, j_north, j_south, nxl, nxlg, nxr, nxrg, &
	799	nyn, nyng, nys, nysg, nzb, nzb_diff_s_inner, &
	800	nzb_diff_s_outer, nzt
	801
	802	USE kinds
[1015]	803
	804	IMPLICIT NONE
	805
[1320]	806	INTEGER(iwp) :: i !:
	807	INTEGER(iwp) :: j !:
	808	INTEGER(iwp) :: k !:
	809	INTEGER(iwp) :: kk !:
	810	INTEGER(iwp) :: max_outer !:
	811	INTEGER(iwp) :: min_inner !:
	812	INTEGER(iwp) :: wall_index !:
	813
	814	REAL(wp) :: a !:
	815	REAL(wp) :: b !:
	816	REAL(wp) :: c1 !:
	817	REAL(wp) :: c2 !:
	818	REAL(wp) :: h1 !:
	819	REAL(wp) :: h2 !:
	820	REAL(wp) :: u_i !:
	821	REAL(wp) :: v_i !:
	822	REAL(wp) :: us_wall !:
	823	REAL(wp) :: vel_total !:
	824	REAL(wp) :: vel_zp !:
	825	REAL(wp) :: ws !:
	826	REAL(wp) :: zp !:
	827	REAL(wp) :: rifs !:
[1015]	828
[1320]	829	REAL(wp), &
	830	DIMENSION(nysg:nyng,nxlg:nxrg) :: &
	831	wall !:
	832
	833	REAL(wp), &
	834	DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) :: &
	835	wall_flux !:
[1015]	836
	837
	838	zp = 0.5 * ( (a+c1) * dy + (b+c2) * dx )
	839	wall_flux = 0.0
	840	wall_index = NINT( a+ 2b + 3c1 + 4*c2 )
	841
	842	min_inner = MINVAL( nzb_diff_s_inner(nys:nyn,nxl:nxr) ) - 1
	843	max_outer = MAXVAL( nzb_diff_s_outer(nys:nyn,nxl:nxr) ) - 2
	844
	845	!$acc kernels present( nzb_diff_s_inner, nzb_diff_s_outer, pt, rif_wall ) &
	846	!$acc present( u, v, w, wall, wall_flux, z0 )
[1128]	847	DO i = i_left, i_right
	848	DO j = j_south, j_north
[1015]	849	DO k = min_inner, max_outer
	850	!
	851	!-- All subsequent variables are computed for scalar locations
[1320]	852	IF ( k >= nzb_diff_s_inner(j,i)-1 .AND. &
[1015]	853	k <= nzb_diff_s_outer(j,i)-2 .AND. wall(j,i) /= 0.0 ) THEN
	854	!
	855	!-- (1) Compute rifs, u_i, v_i, and ws
	856	IF ( k == nzb_diff_s_inner(j,i)-1 ) THEN
	857	kk = nzb_diff_s_inner(j,i)-1
	858	ELSE
	859	kk = k-1
	860	ENDIF
	861	rifs = 0.5 * ( rif_wall(k,j,i,wall_index) + &
	862	a * rif_wall(k,j,i+1,1) + b * rif_wall(k,j+1,i,2) + &
	863	c1 * rif_wall(kk,j,i,3) + c2 * rif_wall(kk,j,i,4) &
	864	)
	865
	866	u_i = 0.5 * ( u(k,j,i) + u(k,j,i+1) )
	867	v_i = 0.5 * ( v(k,j,i) + v(k,j+1,i) )
	868	ws = 0.5 * ( w(k,j,i) + w(k-1,j,i) )
	869	!
	870	!-- (2) Compute wall-parallel absolute velocity vel_total and
	871	!-- interpolate appropriate velocity component vel_zp.
	872	vel_total = SQRT( ws*2 + (a+c1) u_i*2 + (b+c2) v_i**2 )
	873	vel_zp = 0.5 * ( a * u_i + b * v_i + (c1+c2) * ws )
	874	!
	875	!-- (3) Compute wall friction velocity us_wall
	876	IF ( rifs >= 0.0 ) THEN
	877
	878	!
	879	!-- Stable stratification (and neutral)
	880	us_wall = kappa * vel_total / ( LOG( zp / z0(j,i) ) + &
	881	5.0 * rifs * ( zp - z0(j,i) ) / zp &
	882	)
	883	ELSE
	884
	885	!
	886	!-- Unstable stratification
	887	h1 = SQRT( SQRT( 1.0 - 16.0 * rifs ) )
	888	h2 = SQRT( SQRT( 1.0 - 16.0 * rifs * z0(j,i) / zp ) )
	889
	890	us_wall = kappa * vel_total / ( &
	891	LOG( zp / z0(j,i) ) - &
	892	LOG( ( 1.0 + h1 )*2 ( 1.0 + h1**2 ) / ( &
	893	( 1.0 + h2 )*2 ( 1.0 + h2**2 ) ) ) + &
	894	2.0 * ( ATAN( h1 ) - ATAN( h2 ) ) &
	895	)
	896	ENDIF
	897
	898	!
	899	!-- Skip step (4) of wall_fluxes, because here rifs is already
	900	!-- available from (1)
	901	!
	902	!-- (5) Compute wall_flux (u'v', v'u', w'v', or w'u')
	903
	904	IF ( rifs >= 0.0 ) THEN
	905
	906	!
	907	!-- Stable stratification (and neutral)
[1320]	908	wall_flux(k,j,i) = kappa * vel_zp / &
[1015]	909	( LOG( zp/z0(j,i) ) + 5.0rifs ( zp-z0(j,i) ) / zp )
	910	ELSE
	911
	912	!
	913	!-- Unstable stratification
	914	h1 = SQRT( SQRT( 1.0 - 16.0 * rifs ) )
	915	h2 = SQRT( SQRT( 1.0 - 16.0 * rifs * z0(j,i) / zp ) )
	916
	917	wall_flux(k,j,i) = kappa * vel_zp / ( &
	918	LOG( zp / z0(j,i) ) - &
	919	LOG( ( 1.0 + h1 )*2 ( 1.0 + h1**2 ) / ( &
	920	( 1.0 + h2 )*2 ( 1.0 + h2**2 ) ) ) + &
	921	2.0 * ( ATAN( h1 ) - ATAN( h2 ) ) &
	922	)
	923	ENDIF
	924	wall_flux(k,j,i) = - wall_flux(k,j,i) * us_wall
	925
	926	ENDIF
	927
	928	ENDDO
	929	ENDDO
	930	ENDDO
	931	!$acc end kernels
	932
	933	END SUBROUTINE wall_fluxes_e_acc
	934
	935
	936	!------------------------------------------------------------------------------!
[56]	937	! Call for grid point i,j
	938	!------------------------------------------------------------------------------!
	939	SUBROUTINE wall_fluxes_e_ij( i, j, nzb_w, nzt_w, wall_flux, a, b, c1, c2 )
	940
[1320]	941	USE arrays_3d, &
	942	ONLY: rif_wall, u, v, w, z0
	943
	944	USE control_parameters, &
	945	ONLY: kappa
	946
	947	USE grid_variables, &
	948	ONLY: dx, dy
	949
	950	USE indices, &
	951	ONLY: nzb, nzt
	952
	953	USE kinds
[56]	954
	955	IMPLICIT NONE
	956
[1320]	957	INTEGER(iwp) :: i !:
	958	INTEGER(iwp) :: j !:
	959	INTEGER(iwp) :: k !:
	960	INTEGER(iwp) :: kk !:
	961	INTEGER(iwp) :: nzb_w !:
	962	INTEGER(iwp) :: nzt_w !:
	963	INTEGER(iwp) :: wall_index !:
	964
	965	REAL(wp) :: a !:
	966	REAL(wp) :: b !:
	967	REAL(wp) :: c1 !:
	968	REAL(wp) :: c2 !:
	969	REAL(wp) :: h1 !:
	970	REAL(wp) :: h2 !:
	971	REAL(wp) :: u_i !:
	972	REAL(wp) :: v_i !:
	973	REAL(wp) :: us_wall !:
	974	REAL(wp) :: vel_total !:
	975	REAL(wp) :: vel_zp !:
	976	REAL(wp) :: ws !:
	977	REAL(wp) :: zp !:
	978	REAL(wp) :: rifs !:
[56]	979
[1320]	980	REAL(wp), DIMENSION(nzb:nzt+1) :: wall_flux !:
[56]	981
	982
	983	zp = 0.5 * ( (a+c1) * dy + (b+c2) * dx )
	984	wall_flux = 0.0
	985	wall_index = NINT( a+ 2b + 3c1 + 4*c2 )
	986
	987	!
[187]	988	!-- All subsequent variables are computed for scalar locations.
[56]	989	DO k = nzb_w, nzt_w
	990
	991	!
[187]	992	!-- (1) Compute rifs, u_i, v_i, and ws
[56]	993	IF ( k == nzb_w ) THEN
	994	kk = nzb_w
[53]	995	ELSE
[56]	996	kk = k-1
	997	ENDIF
	998	rifs = 0.5 * ( rif_wall(k,j,i,wall_index) + &
	999	a * rif_wall(k,j,i+1,1) + b * rif_wall(k,j+1,i,2) + &
	1000	c1 * rif_wall(kk,j,i,3) + c2 * rif_wall(kk,j,i,4) &
	1001	)
	1002
[187]	1003	u_i = 0.5 * ( u(k,j,i) + u(k,j,i+1) )
	1004	v_i = 0.5 * ( v(k,j,i) + v(k,j+1,i) )
	1005	ws = 0.5 * ( w(k,j,i) + w(k-1,j,i) )
[56]	1006	!
[187]	1007	!-- (2) Compute wall-parallel absolute velocity vel_total and
	1008	!-- interpolate appropriate velocity component vel_zp.
	1009	vel_total = SQRT( ws*2 + (a+c1) u_i*2 + (b+c2) v_i**2 )
	1010	vel_zp = 0.5 * ( a * u_i + b * v_i + (c1+c2) * ws )
	1011	!
	1012	!-- (3) Compute wall friction velocity us_wall
	1013	IF ( rifs >= 0.0 ) THEN
[56]	1014
	1015	!
[187]	1016	!-- Stable stratification (and neutral)
	1017	us_wall = kappa * vel_total / ( LOG( zp / z0(j,i) ) + &
	1018	5.0 * rifs * ( zp - z0(j,i) ) / zp &
	1019	)
	1020	ELSE
	1021
	1022	!
	1023	!-- Unstable stratification
	1024	h1 = SQRT( SQRT( 1.0 - 16.0 * rifs ) )
	1025	h2 = SQRT( SQRT( 1.0 - 16.0 * rifs * z0(j,i) / zp ) )
	1026
[1320]	1027	us_wall = kappa * vel_total / ( &
	1028	LOG( zp / z0(j,i) ) - &
	1029	LOG( ( 1.0 + h1 )*2 ( 1.0 + h1**2 ) / ( &
	1030	( 1.0 + h2 )*2 ( 1.0 + h2**2 ) ) ) + &
	1031	2.0 * ( ATAN( h1 ) - ATAN( h2 ) ) &
[187]	1032	)
	1033	ENDIF
	1034
	1035	!
	1036	!-- Skip step (4) of wall_fluxes, because here rifs is already
	1037	!-- available from (1)
	1038	!
[56]	1039	!-- (5) Compute wall_flux (u'v', v'u', w'v', or w'u')
[187]	1040	!-- First interpolate the velocity (this is different from
	1041	!-- subroutine wall_fluxes because fluxes in subroutine
	1042	!-- wall_fluxes_e are defined at scalar locations).
[1320]	1043	vel_zp = 0.5 * ( a * ( u(k,j,i) + u(k,j,i+1) ) + &
	1044	b * ( v(k,j,i) + v(k,j+1,i) ) + &
	1045	(c1+c2) * ( w(k,j,i) + w(k-1,j,i) ) &
[56]	1046	)
	1047
	1048	IF ( rifs >= 0.0 ) THEN
	1049
	1050	!
	1051	!-- Stable stratification (and neutral)
[1320]	1052	wall_flux(k) = kappa * vel_zp / &
[56]	1053	( LOG( zp/z0(j,i) ) + 5.0rifs ( zp-z0(j,i) ) / zp )
	1054	ELSE
	1055
	1056	!
	1057	!-- Unstable stratification
[187]	1058	h1 = SQRT( SQRT( 1.0 - 16.0 * rifs ) )
	1059	h2 = SQRT( SQRT( 1.0 - 16.0 * rifs * z0(j,i) / zp ) )
[56]	1060
[1320]	1061	wall_flux(k) = kappa * vel_zp / ( &
	1062	LOG( zp / z0(j,i) ) - &
	1063	LOG( ( 1.0 + h1 )*2 ( 1.0 + h1**2 ) / ( &
	1064	( 1.0 + h2 )*2 ( 1.0 + h2**2 ) ) ) + &
	1065	2.0 * ( ATAN( h1 ) - ATAN( h2 ) ) &
[187]	1066	)
[53]	1067	ENDIF
[187]	1068	wall_flux(k) = - wall_flux(k) * us_wall
[53]	1069
[56]	1070	ENDDO
[53]	1071
[56]	1072	END SUBROUTINE wall_fluxes_e_ij
	1073
	1074	END MODULE wall_fluxes_mod

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |