Home

Context Navigation

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

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |