Home

Context Navigation

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

Last change on this file since 2118 was 2118, checked in by raasch, 7 years ago
all OpenACC directives and related parts removed from the code
Property svn:keywords set to `Id`
File size: 28.4 KB

Rev	Line
[1873]	1	!> @file wall_fluxes.f90
[2000]	2	!------------------------------------------------------------------------------!
[1036]	3	! This file is part of PALM.
	4	!
[2000]	5	! PALM is free software: you can redistribute it and/or modify it under the
	6	! terms of the GNU General Public License as published by the Free Software
	7	! Foundation, either version 3 of the License, or (at your option) any later
	8	! version.
[1036]	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	!
[2101]	17	! Copyright 1997-2017 Leibniz Universitaet Hannover
[2000]	18	!------------------------------------------------------------------------------!
[1036]	19	!
[484]	20	! Current revisions:
[52]	21	! -----------------
[2118]	22	! OpenACC versions of subroutines removed
[1354]	23	!
[1321]	24	! Former revisions:
	25	! -----------------
	26	! $Id: wall_fluxes.f90 2118 2017-01-17 16:38:49Z raasch $
	27	!
[2001]	28	! 2000 2016-08-20 18:09:15Z knoop
	29	! Forced header and separation lines into 80 columns
	30	!
[1874]	31	! 1873 2016-04-18 14:50:06Z maronga
	32	! Module renamed (removed _mod)
	33	!
	34	!
[1851]	35	! 1850 2016-04-08 13:29:27Z maronga
	36	! Module renamed
	37	!
	38	!
[1692]	39	! 1691 2015-10-26 16:17:44Z maronga
	40	! Renamed rif_min and rif_max with zeta_min and zeta_max, respectively.
	41	!
[1683]	42	! 1682 2015-10-07 23:56:08Z knoop
	43	! Code annotations made doxygen readable
	44	!
[1375]	45	! 1374 2014-04-25 12:55:07Z raasch
	46	! pt removed from acc-present-list
	47	!
[1354]	48	! 1353 2014-04-08 15:21:23Z heinze
	49	! REAL constants provided with KIND-attribute
	50	!
[1321]	51	! 1320 2014-03-20 08:40:49Z raasch
[1320]	52	! ONLY-attribute added to USE-statements,
	53	! kind-parameters added to all INTEGER and REAL declaration statements,
	54	! kinds are defined in new module kinds,
	55	! old module precision_kind is removed,
	56	! revision history before 2012 removed,
	57	! comment fields (!:) to be used for variable explanations added to
	58	! all variable declaration statements
[198]	59	!
[1258]	60	! 1257 2013-11-08 15:18:40Z raasch
	61	! openacc loop and loop vector clauses removed
	62	!
[1154]	63	! 1153 2013-05-10 14:33:08Z raasch
	64	! code adjustments of accelerator version required by PGI 12.3 / CUDA 5.0
	65	!
[1132]	66	! 1128 2013-04-12 06:19:32Z raasch
	67	! loop index bounds in accelerator version replaced by i_left, i_right, j_south,
	68	! j_north
	69	!
[1037]	70	! 1036 2012-10-22 13:43:42Z raasch
	71	! code put under GPL (PALM 3.9)
	72	!
[1017]	73	! 1015 2012-09-27 09:23:24Z raasch
	74	! accelerator version (*_acc) added
	75	!
[52]	76	! Initial version (2007/03/07)
	77	!
	78	! Description:
	79	! ------------
[1682]	80	!> Calculates momentum fluxes at vertical walls assuming Monin-Obukhov
	81	!> similarity.
	82	!> Indices: usvs a=1, vsus b=1, wsvs c1=1, wsus c2=1 (other=0).
	83	!> The all-gridpoint version of wall_fluxes_e is not used so far, because
	84	!> it gives slightly different results from the ij-version for some unknown
	85	!> reason.
[1691]	86	!>
	87	!> @todo Rename rif to zeta throughout the routine
[52]	88	!------------------------------------------------------------------------------!
[1682]	89	MODULE wall_fluxes_mod
	90
[56]	91	PRIVATE
[2118]	92	PUBLIC wall_fluxes, wall_fluxes_e
[56]	93
	94	INTERFACE wall_fluxes
	95	MODULE PROCEDURE wall_fluxes
	96	MODULE PROCEDURE wall_fluxes_ij
	97	END INTERFACE wall_fluxes
	98
	99	INTERFACE wall_fluxes_e
	100	MODULE PROCEDURE wall_fluxes_e
	101	MODULE PROCEDURE wall_fluxes_e_ij
	102	END INTERFACE wall_fluxes_e
	103
	104	CONTAINS
[52]	105
[56]	106	!------------------------------------------------------------------------------!
[1682]	107	! Description:
	108	! ------------
	109	!> Call for all grid points
[56]	110	!------------------------------------------------------------------------------!
[1320]	111	SUBROUTINE wall_fluxes( wall_flux, a, b, c1, c2, nzb_uvw_inner, &
[56]	112	nzb_uvw_outer, wall )
[52]	113
[1320]	114	USE arrays_3d, &
	115	ONLY: rif_wall, u, v, w, z0, pt
	116
	117	USE control_parameters, &
[1691]	118	ONLY: g, kappa, zeta_max, zeta_min
[1320]	119
	120	USE grid_variables, &
	121	ONLY: dx, dy
	122
	123	USE indices, &
	124	ONLY: nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb, nzt
	125
	126	USE kinds
	127
	128	USE statistics, &
	129	ONLY: hom
[52]	130
[56]	131	IMPLICIT NONE
[52]	132
[1682]	133	INTEGER(iwp) :: i !<
	134	INTEGER(iwp) :: j !<
	135	INTEGER(iwp) :: k !<
	136	INTEGER(iwp) :: wall_index !<
[52]	137
[1320]	138	INTEGER(iwp), &
	139	DIMENSION(nysg:nyng,nxlg:nxrg) :: &
[1682]	140	nzb_uvw_inner !<
[1320]	141	INTEGER(iwp), &
	142	DIMENSION(nysg:nyng,nxlg:nxrg) :: &
[1682]	143	nzb_uvw_outer !<
[1320]	144
[1682]	145	REAL(wp) :: a !<
	146	REAL(wp) :: b !<
	147	REAL(wp) :: c1 !<
	148	REAL(wp) :: c2 !<
	149	REAL(wp) :: h1 !<
	150	REAL(wp) :: h2 !<
	151	REAL(wp) :: zp !<
	152	REAL(wp) :: pts !<
	153	REAL(wp) :: pt_i !<
	154	REAL(wp) :: rifs !<
	155	REAL(wp) :: u_i !<
	156	REAL(wp) :: v_i !<
	157	REAL(wp) :: us_wall !<
	158	REAL(wp) :: vel_total !<
	159	REAL(wp) :: ws !<
	160	REAL(wp) :: wspts !<
[52]	161
[1320]	162	REAL(wp), &
	163	DIMENSION(nysg:nyng,nxlg:nxrg) :: &
[1682]	164	wall !<
[1320]	165
	166	REAL(wp), &
	167	DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) :: &
[1682]	168	wall_flux !<
[52]	169
	170
[1353]	171	zp = 0.5_wp * ( (a+c1) * dy + (b+c2) * dx )
	172	wall_flux = 0.0_wp
[56]	173	wall_index = NINT( a+ 2b + 3c1 + 4*c2 )
	174
[75]	175	DO i = nxl, nxr
	176	DO j = nys, nyn
[56]	177
[1353]	178	IF ( wall(j,i) /= 0.0_wp ) THEN
[52]	179	!
[56]	180	!-- All subsequent variables are computed for the respective
[187]	181	!-- location where the respective flux is defined.
[56]	182	DO k = nzb_uvw_inner(j,i)+1, nzb_uvw_outer(j,i)
[53]	183
[52]	184	!
[56]	185	!-- (1) Compute rifs, u_i, v_i, ws, pt' and w'pt'
	186	rifs = rif_wall(k,j,i,wall_index)
[53]	187
[1353]	188	u_i = a * u(k,j,i) + c1 * 0.25_wp * &
[56]	189	( u(k+1,j,i+1) + u(k+1,j,i) + u(k,j,i+1) + u(k,j,i) )
[53]	190
[1353]	191	v_i = b * v(k,j,i) + c2 * 0.25_wp * &
[56]	192	( v(k+1,j+1,i) + v(k+1,j,i) + v(k,j+1,i) + v(k,j,i) )
[53]	193
[1353]	194	ws = ( c1 + c2 ) * w(k,j,i) + 0.25_wp * ( &
[56]	195	a * ( w(k-1,j,i-1) + w(k-1,j,i) + w(k,j,i-1) + w(k,j,i) ) &
	196	+ b * ( w(k-1,j-1,i) + w(k-1,j,i) + w(k,j-1,i) + w(k,j,i) ) &
[1353]	197	)
	198	pt_i = 0.5_wp * ( pt(k,j,i) + a * pt(k,j,i-1) + &
[56]	199	b * pt(k,j-1,i) + ( c1 + c2 ) * pt(k+1,j,i) )
[53]	200
[56]	201	pts = pt_i - hom(k,1,4,0)
	202	wspts = ws * pts
[53]	203
[52]	204	!
[56]	205	!-- (2) Compute wall-parallel absolute velocity vel_total
	206	vel_total = SQRT( ws*2 + (a+c1) u_i*2 + (b+c2) v_i**2 )
[53]	207
[52]	208	!
[56]	209	!-- (3) Compute wall friction velocity us_wall
[1353]	210	IF ( rifs >= 0.0_wp ) THEN
[53]	211
[52]	212	!
[56]	213	!-- Stable stratification (and neutral)
	214	us_wall = kappa * vel_total / ( LOG( zp / z0(j,i) ) + &
[1353]	215	5.0_wp * rifs * ( zp - z0(j,i) ) / zp &
[56]	216	)
	217	ELSE
[53]	218
[52]	219	!
[56]	220	!-- Unstable stratification
[1353]	221	h1 = SQRT( SQRT( 1.0_wp - 16.0_wp * rifs ) )
	222	h2 = SQRT( SQRT( 1.0_wp - 16.0_wp * rifs * z0(j,i) / zp ) )
[53]	223
[187]	224	us_wall = kappa * vel_total / ( &
	225	LOG( zp / z0(j,i) ) - &
[1353]	226	LOG( ( 1.0_wp + h1 )*2 ( 1.0_wp + h1**2 ) / ( &
	227	( 1.0_wp + h2 )*2 ( 1.0_wp + h2**2 ) ) ) +&
	228	2.0_wp * ( ATAN( h1 ) - ATAN( h2 ) ) &
[187]	229	)
[56]	230	ENDIF
[53]	231
[52]	232	!
[56]	233	!-- (4) Compute zp/L (corresponds to neutral Richardson flux
	234	!-- number rifs)
[1353]	235	rifs = -1.0_wp * zp * kappa * g * wspts / &
	236	( pt_i * ( us_wall**3 + 1E-30 ) )
[53]	237
[52]	238	!
[56]	239	!-- Limit the value range of the Richardson numbers.
	240	!-- This is necessary for very small velocities (u,w --> 0),
	241	!-- because the absolute value of rif can then become very
	242	!-- large, which in consequence would result in very large
	243	!-- shear stresses and very small momentum fluxes (both are
	244	!-- generally unrealistic).
[1691]	245	IF ( rifs < zeta_min ) rifs = zeta_min
	246	IF ( rifs > zeta_max ) rifs = zeta_max
[53]	247
[52]	248	!
[56]	249	!-- (5) Compute wall_flux (u'v', v'u', w'v', or w'u')
[1353]	250	IF ( rifs >= 0.0_wp ) THEN
[53]	251
[52]	252	!
[56]	253	!-- Stable stratification (and neutral)
	254	wall_flux(k,j,i) = kappa * &
	255	( au(k,j,i) + bv(k,j,i) + (c1+c2)*w(k,j,i) ) / &
	256	( LOG( zp / z0(j,i) ) + &
[1353]	257	5.0_wp * rifs * ( zp - z0(j,i) ) / zp &
[56]	258	)
	259	ELSE
[53]	260
[52]	261	!
[56]	262	!-- Unstable stratification
[1353]	263	h1 = SQRT( SQRT( 1.0_wp - 16.0_wp * rifs ) )
	264	h2 = SQRT( SQRT( 1.0_wp - 16.0_wp * rifs * z0(j,i) / zp ) )
[53]	265
[187]	266	wall_flux(k,j,i) = kappa * &
	267	( au(k,j,i) + bv(k,j,i) + (c1+c2)*w(k,j,i) ) / ( &
	268	LOG( zp / z0(j,i) ) - &
[1353]	269	LOG( ( 1.0_wp + h1 )*2 ( 1.0_wp + h1**2 ) / ( &
	270	( 1.0_wp + h2 )*2 ( 1.0_wp + h2**2 ) ) ) +&
	271	2.0_wp * ( ATAN( h1 ) - ATAN( h2 ) ) &
[187]	272	)
[56]	273	ENDIF
[187]	274	wall_flux(k,j,i) = -wall_flux(k,j,i) * us_wall
[56]	275
	276	!
	277	!-- store rifs for next time step
	278	rif_wall(k,j,i,wall_index) = rifs
	279
	280	ENDDO
	281
	282	ENDIF
	283
	284	ENDDO
	285	ENDDO
	286
	287	END SUBROUTINE wall_fluxes
	288
	289
[1015]	290	!------------------------------------------------------------------------------!
[1682]	291	! Description:
	292	! ------------
	293	!> Call for all grid point i,j
[56]	294	!------------------------------------------------------------------------------!
	295	SUBROUTINE wall_fluxes_ij( i, j, nzb_w, nzt_w, wall_flux, a, b, c1, c2 )
	296
[1320]	297	USE arrays_3d, &
	298	ONLY: rif_wall, pt, u, v, w, z0
	299
	300	USE control_parameters, &
[1691]	301	ONLY: g, kappa, zeta_max, zeta_min
[1320]	302
	303	USE grid_variables, &
	304	ONLY: dx, dy
	305
	306	USE indices, &
	307	ONLY: nzb, nzt
	308
	309	USE kinds
	310
	311	USE statistics, &
	312	ONLY: hom
[56]	313
	314	IMPLICIT NONE
	315
[1682]	316	INTEGER(iwp) :: i !<
	317	INTEGER(iwp) :: j !<
	318	INTEGER(iwp) :: k !<
	319	INTEGER(iwp) :: nzb_w !<
	320	INTEGER(iwp) :: nzt_w !<
	321	INTEGER(iwp) :: wall_index !<
[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 !<
[56]	339
[1682]	340	REAL(wp), DIMENSION(nzb:nzt+1) :: wall_flux !<
[56]	341
	342
[1353]	343	zp = 0.5_wp * ( (a+c1) * dy + (b+c2) * dx )
	344	wall_flux = 0.0_wp
[56]	345	wall_index = NINT( a+ 2b + 3c1 + 4*c2 )
	346
	347	!
	348	!-- All subsequent variables are computed for the respective location where
[187]	349	!-- the respective flux is defined.
[56]	350	DO k = nzb_w, nzt_w
	351
	352	!
	353	!-- (1) Compute rifs, u_i, v_i, ws, pt' and w'pt'
	354	rifs = rif_wall(k,j,i,wall_index)
	355
[1353]	356	u_i = a * u(k,j,i) + c1 * 0.25_wp * &
[56]	357	( u(k+1,j,i+1) + u(k+1,j,i) + u(k,j,i+1) + u(k,j,i) )
	358
[1353]	359	v_i = b * v(k,j,i) + c2 * 0.25_wp * &
[56]	360	( v(k+1,j+1,i) + v(k+1,j,i) + v(k,j+1,i) + v(k,j,i) )
	361
[1353]	362	ws = ( c1 + c2 ) * w(k,j,i) + 0.25_wp * ( &
[56]	363	a * ( w(k-1,j,i-1) + w(k-1,j,i) + w(k,j,i-1) + w(k,j,i) ) &
	364	+ b * ( w(k-1,j-1,i) + w(k-1,j,i) + w(k,j-1,i) + w(k,j,i) ) &
[1353]	365	)
	366	pt_i = 0.5_wp * ( pt(k,j,i) + a * pt(k,j,i-1) + b * pt(k,j-1,i) &
[56]	367	+ ( c1 + c2 ) * pt(k+1,j,i) )
	368
	369	pts = pt_i - hom(k,1,4,0)
	370	wspts = ws * pts
	371
	372	!
	373	!-- (2) Compute wall-parallel absolute velocity vel_total
	374	vel_total = SQRT( ws*2 + ( a+c1 ) u_i*2 + ( b+c2 ) v_i**2 )
	375
	376	!
	377	!-- (3) Compute wall friction velocity us_wall
[1353]	378	IF ( rifs >= 0.0_wp ) THEN
[56]	379
	380	!
	381	!-- Stable stratification (and neutral)
	382	us_wall = kappa * vel_total / ( LOG( zp / z0(j,i) ) + &
[1353]	383	5.0_wp * rifs * ( zp - z0(j,i) ) / zp &
[56]	384	)
	385	ELSE
	386
	387	!
	388	!-- Unstable stratification
[1353]	389	h1 = SQRT( SQRT( 1.0_wp - 16.0_wp * rifs ) )
	390	h2 = SQRT( SQRT( 1.0_wp - 16.0_wp * rifs * z0(j,i) / zp ) )
[56]	391
[1320]	392	us_wall = kappa * vel_total / ( &
	393	LOG( zp / z0(j,i) ) - &
[1353]	394	LOG( ( 1.0_wp + h1 )*2 ( 1.0_wp + h1**2 ) / ( &
	395	( 1.0_wp + h2 )*2 ( 1.0_wp + h2**2 ) ) ) + &
	396	2.0_wp * ( ATAN( h1 ) - ATAN( h2 ) ) &
[187]	397	)
[56]	398	ENDIF
	399
	400	!
	401	!-- (4) Compute zp/L (corresponds to neutral Richardson flux number
	402	!-- rifs)
[1353]	403	rifs = -1.0_wp * zp * kappa * g * wspts / &
	404	( pt_i * (us_wall**3 + 1E-30) )
[56]	405
	406	!
	407	!-- Limit the value range of the Richardson numbers.
	408	!-- This is necessary for very small velocities (u,w --> 0), because
	409	!-- the absolute value of rif can then become very large, which in
	410	!-- consequence would result in very large shear stresses and very
	411	!-- small momentum fluxes (both are generally unrealistic).
[1691]	412	IF ( rifs < zeta_min ) rifs = zeta_min
	413	IF ( rifs > zeta_max ) rifs = zeta_max
[56]	414
	415	!
	416	!-- (5) Compute wall_flux (u'v', v'u', w'v', or w'u')
[1353]	417	IF ( rifs >= 0.0_wp ) THEN
[56]	418
	419	!
	420	!-- Stable stratification (and neutral)
[1320]	421	wall_flux(k) = kappa * &
	422	( au(k,j,i) + bv(k,j,i) + (c1+c2)*w(k,j,i) ) / &
	423	( LOG( zp / z0(j,i) ) + &
[1353]	424	5.0_wp * rifs * ( zp - z0(j,i) ) / zp &
[53]	425	)
[52]	426	ELSE
[53]	427
[56]	428	!
	429	!-- Unstable stratification
[1353]	430	h1 = SQRT( SQRT( 1.0_wp - 16.0_wp * rifs ) )
	431	h2 = SQRT( SQRT( 1.0_wp - 16.0_wp * rifs * z0(j,i) / zp ) )
[52]	432
[1320]	433	wall_flux(k) = kappa * &
	434	( au(k,j,i) + bv(k,j,i) + (c1+c2)*w(k,j,i) ) / ( &
	435	LOG( zp / z0(j,i) ) - &
[1353]	436	LOG( ( 1.0_wp + h1 )*2 ( 1.0_wp + h1**2 ) / ( &
	437	( 1.0_wp + h2 )*2 ( 1.0_wp + h2**2 ) ) ) + &
	438	2.0_wp * ( ATAN( h1 ) - ATAN( h2 ) ) &
[187]	439	)
[56]	440	ENDIF
[187]	441	wall_flux(k) = -wall_flux(k) * us_wall
[53]	442
[56]	443	!
	444	!-- store rifs for next time step
	445	rif_wall(k,j,i,wall_index) = rifs
[53]	446
[56]	447	ENDDO
[53]	448
[56]	449	END SUBROUTINE wall_fluxes_ij
[53]	450
[56]	451
	452
[53]	453	!------------------------------------------------------------------------------!
	454	! Description:
	455	! ------------
[1682]	456	!> Call for all grid points
	457	!> Calculates momentum fluxes at vertical walls for routine production_e
	458	!> assuming Monin-Obukhov similarity.
	459	!> Indices: usvs a=1, vsus b=1, wsvs c1=1, wsus c2=1 (other=0).
[53]	460	!------------------------------------------------------------------------------!
	461
[1682]	462	SUBROUTINE wall_fluxes_e( wall_flux, a, b, c1, c2, wall )
	463
	464
[1320]	465	USE arrays_3d, &
	466	ONLY: rif_wall, u, v, w, z0
	467
	468	USE control_parameters, &
	469	ONLY: kappa
	470
	471	USE grid_variables, &
	472	ONLY: dx, dy
	473
	474	USE indices, &
	475	ONLY: nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb, &
	476	nzb_diff_s_inner, nzb_diff_s_outer, nzt
	477
	478	USE kinds
[53]	479
[56]	480	IMPLICIT NONE
[53]	481
[1682]	482	INTEGER(iwp) :: i !<
	483	INTEGER(iwp) :: j !<
	484	INTEGER(iwp) :: k !<
	485	INTEGER(iwp) :: kk !<
	486	INTEGER(iwp) :: wall_index !<
[1320]	487
[1682]	488	REAL(wp) :: a !<
	489	REAL(wp) :: b !<
	490	REAL(wp) :: c1 !<
	491	REAL(wp) :: c2 !<
	492	REAL(wp) :: h1 !<
	493	REAL(wp) :: h2 !<
	494	REAL(wp) :: u_i !<
	495	REAL(wp) :: v_i !<
	496	REAL(wp) :: us_wall !<
	497	REAL(wp) :: vel_total !<
	498	REAL(wp) :: vel_zp !<
	499	REAL(wp) :: ws !<
	500	REAL(wp) :: zp !<
	501	REAL(wp) :: rifs !<
[53]	502
[1320]	503	REAL(wp), &
	504	DIMENSION(nysg:nyng,nxlg:nxrg) :: &
[1682]	505	wall !<
[1320]	506
	507	REAL(wp), &
	508	DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) :: &
[1682]	509	wall_flux !<
[53]	510
	511
[1353]	512	zp = 0.5_wp * ( (a+c1) * dy + (b+c2) * dx )
	513	wall_flux = 0.0_wp
[56]	514	wall_index = NINT( a+ 2b + 3c1 + 4*c2 )
[53]	515
[56]	516	DO i = nxl, nxr
	517	DO j = nys, nyn
	518
[1353]	519	IF ( wall(j,i) /= 0.0_wp ) THEN
[53]	520	!
[187]	521	!-- All subsequent variables are computed for scalar locations.
[56]	522	DO k = nzb_diff_s_inner(j,i)-1, nzb_diff_s_outer(j,i)-2
[53]	523	!
[187]	524	!-- (1) Compute rifs, u_i, v_i, and ws
[56]	525	IF ( k == nzb_diff_s_inner(j,i)-1 ) THEN
	526	kk = nzb_diff_s_inner(j,i)-1
	527	ELSE
	528	kk = k-1
	529	ENDIF
[1353]	530	rifs = 0.5_wp * ( rif_wall(k,j,i,wall_index) + &
	531	a * rif_wall(k,j,i+1,1) + &
	532	b * rif_wall(k,j+1,i,2) + &
	533	c1 * rif_wall(kk,j,i,3) + &
	534	c2 * rif_wall(kk,j,i,4) &
	535	)
[53]	536
[1353]	537	u_i = 0.5_wp * ( u(k,j,i) + u(k,j,i+1) )
	538	v_i = 0.5_wp * ( v(k,j,i) + v(k,j+1,i) )
	539	ws = 0.5_wp * ( w(k,j,i) + w(k-1,j,i) )
[53]	540	!
[187]	541	!-- (2) Compute wall-parallel absolute velocity vel_total and
	542	!-- interpolate appropriate velocity component vel_zp.
	543	vel_total = SQRT( ws*2 + (a+c1) u_i*2 + (b+c2) v_i**2 )
[1353]	544	vel_zp = 0.5_wp * ( a * u_i + b * v_i + (c1+c2) * ws )
[187]	545	!
	546	!-- (3) Compute wall friction velocity us_wall
[1353]	547	IF ( rifs >= 0.0_wp ) THEN
[53]	548
	549	!
[187]	550	!-- Stable stratification (and neutral)
	551	us_wall = kappa * vel_total / ( LOG( zp / z0(j,i) ) + &
[1353]	552	5.0_wp * rifs * ( zp - z0(j,i) ) / zp &
[187]	553	)
	554	ELSE
	555
	556	!
	557	!-- Unstable stratification
[1353]	558	h1 = SQRT( SQRT( 1.0_wp - 16.0_wp * rifs ) )
	559	h2 = SQRT( SQRT( 1.0_wp - 16.0_wp * rifs * z0(j,i) / zp ) )
[187]	560
	561	us_wall = kappa * vel_total / ( &
	562	LOG( zp / z0(j,i) ) - &
[1353]	563	LOG( ( 1.0_wp + h1 )*2 ( 1.0_wp + h1**2 ) / ( &
	564	( 1.0_wp + h2 )*2 ( 1.0_wp + h2**2 ) ) ) +&
	565	2.0_wp * ( ATAN( h1 ) - ATAN( h2 ) ) &
[187]	566	)
	567	ENDIF
	568
	569	!
	570	!-- Skip step (4) of wall_fluxes, because here rifs is already
	571	!-- available from (1)
	572	!
[56]	573	!-- (5) Compute wall_flux (u'v', v'u', w'v', or w'u')
[55]	574
[1353]	575	IF ( rifs >= 0.0_wp ) THEN
[53]	576
	577	!
[56]	578	!-- Stable stratification (and neutral)
[1353]	579	wall_flux(k,j,i) = kappa * vel_zp / ( LOG( zp/z0(j,i) ) +&
	580	5.0_wp * rifs * ( zp-z0(j,i) ) / zp )
[56]	581	ELSE
[53]	582
	583	!
[56]	584	!-- Unstable stratification
[1353]	585	h1 = SQRT( SQRT( 1.0_wp - 16.0_wp * rifs ) )
	586	h2 = SQRT( SQRT( 1.0_wp - 16.0_wp * rifs * z0(j,i) / zp ) )
[53]	587
[187]	588	wall_flux(k,j,i) = kappa * vel_zp / ( &
	589	LOG( zp / z0(j,i) ) - &
[1353]	590	LOG( ( 1.0_wp + h1 )*2 ( 1.0_wp + h1**2 ) / ( &
	591	( 1.0_wp + h2 )*2 ( 1.0_wp + h2**2 ) ) ) +&
	592	2.0_wp * ( ATAN( h1 ) - ATAN( h2 ) ) &
[187]	593	)
[56]	594	ENDIF
[187]	595	wall_flux(k,j,i) = - wall_flux(k,j,i) * us_wall
[56]	596
	597	ENDDO
	598
	599	ENDIF
	600
	601	ENDDO
	602	ENDDO
	603
	604	END SUBROUTINE wall_fluxes_e
	605
	606
[1015]	607	!------------------------------------------------------------------------------!
	608	! Description:
	609	! ------------
[1682]	610	!> Call for grid point i,j
[56]	611	!------------------------------------------------------------------------------!
	612	SUBROUTINE wall_fluxes_e_ij( i, j, nzb_w, nzt_w, wall_flux, a, b, c1, c2 )
	613
[1320]	614	USE arrays_3d, &
	615	ONLY: rif_wall, u, v, w, z0
	616
	617	USE control_parameters, &
	618	ONLY: kappa
	619
	620	USE grid_variables, &
	621	ONLY: dx, dy
	622
	623	USE indices, &
	624	ONLY: nzb, nzt
	625
	626	USE kinds
[56]	627
	628	IMPLICIT NONE
	629
[1682]	630	INTEGER(iwp) :: i !<
	631	INTEGER(iwp) :: j !<
	632	INTEGER(iwp) :: k !<
	633	INTEGER(iwp) :: kk !<
	634	INTEGER(iwp) :: nzb_w !<
	635	INTEGER(iwp) :: nzt_w !<
	636	INTEGER(iwp) :: wall_index !<
[1320]	637
[1682]	638	REAL(wp) :: a !<
	639	REAL(wp) :: b !<
	640	REAL(wp) :: c1 !<
	641	REAL(wp) :: c2 !<
	642	REAL(wp) :: h1 !<
	643	REAL(wp) :: h2 !<
	644	REAL(wp) :: u_i !<
	645	REAL(wp) :: v_i !<
	646	REAL(wp) :: us_wall !<
	647	REAL(wp) :: vel_total !<
	648	REAL(wp) :: vel_zp !<
	649	REAL(wp) :: ws !<
	650	REAL(wp) :: zp !<
	651	REAL(wp) :: rifs !<
[56]	652
[1682]	653	REAL(wp), DIMENSION(nzb:nzt+1) :: wall_flux !<
[56]	654
	655
[1353]	656	zp = 0.5_wp * ( (a+c1) * dy + (b+c2) * dx )
	657	wall_flux = 0.0_wp
[56]	658	wall_index = NINT( a+ 2b + 3c1 + 4*c2 )
	659
	660	!
[187]	661	!-- All subsequent variables are computed for scalar locations.
[56]	662	DO k = nzb_w, nzt_w
	663
	664	!
[187]	665	!-- (1) Compute rifs, u_i, v_i, and ws
[56]	666	IF ( k == nzb_w ) THEN
	667	kk = nzb_w
[53]	668	ELSE
[56]	669	kk = k-1
	670	ENDIF
[1353]	671	rifs = 0.5_wp * ( rif_wall(k,j,i,wall_index) + &
	672	a * rif_wall(k,j,i+1,1) + &
	673	b * rif_wall(k,j+1,i,2) + &
	674	c1 * rif_wall(kk,j,i,3) + &
	675	c2 * rif_wall(kk,j,i,4) &
	676	)
[56]	677
[1353]	678	u_i = 0.5_wp * ( u(k,j,i) + u(k,j,i+1) )
	679	v_i = 0.5_wp * ( v(k,j,i) + v(k,j+1,i) )
	680	ws = 0.5_wp * ( w(k,j,i) + w(k-1,j,i) )
[56]	681	!
[187]	682	!-- (2) Compute wall-parallel absolute velocity vel_total and
	683	!-- interpolate appropriate velocity component vel_zp.
	684	vel_total = SQRT( ws*2 + (a+c1) u_i*2 + (b+c2) v_i**2 )
[1353]	685	vel_zp = 0.5_wp * ( a * u_i + b * v_i + (c1+c2) * ws )
[187]	686	!
	687	!-- (3) Compute wall friction velocity us_wall
[1353]	688	IF ( rifs >= 0.0_wp ) THEN
[56]	689
	690	!
[187]	691	!-- Stable stratification (and neutral)
	692	us_wall = kappa * vel_total / ( LOG( zp / z0(j,i) ) + &
[1353]	693	5.0_wp * rifs * ( zp - z0(j,i) ) / zp &
[187]	694	)
	695	ELSE
	696
	697	!
	698	!-- Unstable stratification
[1353]	699	h1 = SQRT( SQRT( 1.0_wp - 16.0_wp * rifs ) )
	700	h2 = SQRT( SQRT( 1.0_wp - 16.0_wp * rifs * z0(j,i) / zp ) )
[187]	701
[1320]	702	us_wall = kappa * vel_total / ( &
	703	LOG( zp / z0(j,i) ) - &
[1353]	704	LOG( ( 1.0_wp + h1 )*2 ( 1.0_wp + h1**2 ) / ( &
	705	( 1.0_wp + h2 )*2 ( 1.0_wp + h2**2 ) ) ) + &
	706	2.0_wp * ( ATAN( h1 ) - ATAN( h2 ) ) &
[187]	707	)
	708	ENDIF
	709
	710	!
	711	!-- Skip step (4) of wall_fluxes, because here rifs is already
	712	!-- available from (1)
	713	!
[56]	714	!-- (5) Compute wall_flux (u'v', v'u', w'v', or w'u')
[187]	715	!-- First interpolate the velocity (this is different from
	716	!-- subroutine wall_fluxes because fluxes in subroutine
	717	!-- wall_fluxes_e are defined at scalar locations).
[1353]	718	vel_zp = 0.5_wp * ( a * ( u(k,j,i) + u(k,j,i+1) ) + &
	719	b * ( v(k,j,i) + v(k,j+1,i) ) + &
	720	(c1+c2) * ( w(k,j,i) + w(k-1,j,i) ) &
	721	)
[56]	722
[1353]	723	IF ( rifs >= 0.0_wp ) THEN
[56]	724
	725	!
	726	!-- Stable stratification (and neutral)
[1320]	727	wall_flux(k) = kappa * vel_zp / &
[1353]	728	( LOG( zp/z0(j,i) ) + 5.0_wp * rifs * ( zp-z0(j,i) ) / zp )
[56]	729	ELSE
	730
	731	!
	732	!-- Unstable stratification
[1353]	733	h1 = SQRT( SQRT( 1.0_wp - 16.0_wp * rifs ) )
	734	h2 = SQRT( SQRT( 1.0_wp - 16.0_wp * rifs * z0(j,i) / zp ) )
[56]	735
[1320]	736	wall_flux(k) = kappa * vel_zp / ( &
	737	LOG( zp / z0(j,i) ) - &
[1353]	738	LOG( ( 1.0_wp + h1 )*2 ( 1.0_wp + h1**2 ) / ( &
	739	( 1.0_wp + h2 )*2 ( 1.0_wp + h2**2 ) ) ) + &
	740	2.0_wp * ( ATAN( h1 ) - ATAN( h2 ) ) &
	741	)
[53]	742	ENDIF
[187]	743	wall_flux(k) = - wall_flux(k) * us_wall
[53]	744
[56]	745	ENDDO
[53]	746
[56]	747	END SUBROUTINE wall_fluxes_e_ij
	748
	749	END MODULE wall_fluxes_mod

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |