Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |