Home

Context Navigation

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

Last change on this file since 2142 was 2119, checked in by raasch, 8 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 28.5 KB

Line
1	!> @file wall_fluxes.f90
2	!------------------------------------------------------------------------------!
3	! This file is part of PALM.
4	!
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.
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	!
17	! Copyright 1997-2017 Leibniz Universitaet Hannover
18	!------------------------------------------------------------------------------!
19	!
20	! Current revisions:
21	! -----------------
22	!
23	!
24	! Former revisions:
25	! -----------------
26	! $Id: wall_fluxes.f90 2119 2017-01-17 16:51:50Z gronemeier $
27	!
28	! 2118 2017-01-17 16:38:49Z raasch
29	! OpenACC versions of subroutines removed
30	!
31	! 2000 2016-08-20 18:09:15Z knoop
32	! Forced header and separation lines into 80 columns
33	!
34	! 1873 2016-04-18 14:50:06Z maronga
35	! Module renamed (removed _mod)
36	!
37	!
38	! 1850 2016-04-08 13:29:27Z maronga
39	! Module renamed
40	!
41	!
42	! 1691 2015-10-26 16:17:44Z maronga
43	! Renamed rif_min and rif_max with zeta_min and zeta_max, respectively.
44	!
45	! 1682 2015-10-07 23:56:08Z knoop
46	! Code annotations made doxygen readable
47	!
48	! 1374 2014-04-25 12:55:07Z raasch
49	! pt removed from acc-present-list
50	!
51	! 1353 2014-04-08 15:21:23Z heinze
52	! REAL constants provided with KIND-attribute
53	!
54	! 1320 2014-03-20 08:40:49Z raasch
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
62	!
63	! 1257 2013-11-08 15:18:40Z raasch
64	! openacc loop and loop vector clauses removed
65	!
66	! 1153 2013-05-10 14:33:08Z raasch
67	! code adjustments of accelerator version required by PGI 12.3 / CUDA 5.0
68	!
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	!
73	! 1036 2012-10-22 13:43:42Z raasch
74	! code put under GPL (PALM 3.9)
75	!
76	! 1015 2012-09-27 09:23:24Z raasch
77	! accelerator version (*_acc) added
78	!
79	! Initial version (2007/03/07)
80	!
81	! Description:
82	! ------------
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.
89	!>
90	!> @todo Rename rif to zeta throughout the routine
91	!------------------------------------------------------------------------------!
92	MODULE wall_fluxes_mod
93
94	PRIVATE
95	PUBLIC wall_fluxes, wall_fluxes_e
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
108
109	!------------------------------------------------------------------------------!
110	! Description:
111	! ------------
112	!> Call for all grid points
113	!------------------------------------------------------------------------------!
114	SUBROUTINE wall_fluxes( wall_flux, a, b, c1, c2, nzb_uvw_inner, &
115	nzb_uvw_outer, wall )
116
117	USE arrays_3d, &
118	ONLY: rif_wall, u, v, w, z0, pt
119
120	USE control_parameters, &
121	ONLY: g, kappa, zeta_max, zeta_min
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
133
134	IMPLICIT NONE
135
136	INTEGER(iwp) :: i !<
137	INTEGER(iwp) :: j !<
138	INTEGER(iwp) :: k !<
139	INTEGER(iwp) :: wall_index !<
140
141	INTEGER(iwp), &
142	DIMENSION(nysg:nyng,nxlg:nxrg) :: &
143	nzb_uvw_inner !<
144	INTEGER(iwp), &
145	DIMENSION(nysg:nyng,nxlg:nxrg) :: &
146	nzb_uvw_outer !<
147
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 !<
164
165	REAL(wp), &
166	DIMENSION(nysg:nyng,nxlg:nxrg) :: &
167	wall !<
168
169	REAL(wp), &
170	DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) :: &
171	wall_flux !<
172
173
174	zp = 0.5_wp * ( (a+c1) * dy + (b+c2) * dx )
175	wall_flux = 0.0_wp
176	wall_index = NINT( a+ 2b + 3c1 + 4*c2 )
177
178	DO i = nxl, nxr
179	DO j = nys, nyn
180
181	IF ( wall(j,i) /= 0.0_wp ) THEN
182	!
183	!-- All subsequent variables are computed for the respective
184	!-- location where the respective flux is defined.
185	DO k = nzb_uvw_inner(j,i)+1, nzb_uvw_outer(j,i)
186
187	!
188	!-- (1) Compute rifs, u_i, v_i, ws, pt' and w'pt'
189	rifs = rif_wall(k,j,i,wall_index)
190
191	u_i = a * u(k,j,i) + c1 * 0.25_wp * &
192	( u(k+1,j,i+1) + u(k+1,j,i) + u(k,j,i+1) + u(k,j,i) )
193
194	v_i = b * v(k,j,i) + c2 * 0.25_wp * &
195	( v(k+1,j+1,i) + v(k+1,j,i) + v(k,j+1,i) + v(k,j,i) )
196
197	ws = ( c1 + c2 ) * w(k,j,i) + 0.25_wp * ( &
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) ) &
200	)
201	pt_i = 0.5_wp * ( pt(k,j,i) + a * pt(k,j,i-1) + &
202	b * pt(k,j-1,i) + ( c1 + c2 ) * pt(k+1,j,i) )
203
204	pts = pt_i - hom(k,1,4,0)
205	wspts = ws * pts
206
207	!
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 )
210
211	!
212	!-- (3) Compute wall friction velocity us_wall
213	IF ( rifs >= 0.0_wp ) THEN
214
215	!
216	!-- Stable stratification (and neutral)
217	us_wall = kappa * vel_total / ( LOG( zp / z0(j,i) ) + &
218	5.0_wp * rifs * ( zp - z0(j,i) ) / zp &
219	)
220	ELSE
221
222	!
223	!-- Unstable stratification
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 ) )
226
227	us_wall = kappa * vel_total / ( &
228	LOG( zp / z0(j,i) ) - &
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 ) ) &
232	)
233	ENDIF
234
235	!
236	!-- (4) Compute zp/L (corresponds to neutral Richardson flux
237	!-- number rifs)
238	rifs = -1.0_wp * zp * kappa * g * wspts / &
239	( pt_i * ( us_wall**3 + 1E-30 ) )
240
241	!
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).
248	IF ( rifs < zeta_min ) rifs = zeta_min
249	IF ( rifs > zeta_max ) rifs = zeta_max
250
251	!
252	!-- (5) Compute wall_flux (u'v', v'u', w'v', or w'u')
253	IF ( rifs >= 0.0_wp ) THEN
254
255	!
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) ) + &
260	5.0_wp * rifs * ( zp - z0(j,i) ) / zp &
261	)
262	ELSE
263
264	!
265	!-- Unstable stratification
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 ) )
268
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) ) - &
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 ) ) &
275	)
276	ENDIF
277	wall_flux(k,j,i) = -wall_flux(k,j,i) * us_wall
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
293	!------------------------------------------------------------------------------!
294	! Description:
295	! ------------
296	!> Call for all grid point i,j
297	!------------------------------------------------------------------------------!
298	SUBROUTINE wall_fluxes_ij( i, j, nzb_w, nzt_w, wall_flux, a, b, c1, c2 )
299
300	USE arrays_3d, &
301	ONLY: rif_wall, pt, u, v, w, z0
302
303	USE control_parameters, &
304	ONLY: g, kappa, zeta_max, zeta_min
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
316
317	IMPLICIT NONE
318
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 !<
325
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 !<
342
343	REAL(wp), DIMENSION(nzb:nzt+1) :: wall_flux !<
344
345
346	zp = 0.5_wp * ( (a+c1) * dy + (b+c2) * dx )
347	wall_flux = 0.0_wp
348	wall_index = NINT( a+ 2b + 3c1 + 4*c2 )
349
350	!
351	!-- All subsequent variables are computed for the respective location where
352	!-- the respective flux is defined.
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
359	u_i = a * u(k,j,i) + c1 * 0.25_wp * &
360	( u(k+1,j,i+1) + u(k+1,j,i) + u(k,j,i+1) + u(k,j,i) )
361
362	v_i = b * v(k,j,i) + c2 * 0.25_wp * &
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_wp * ( &
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	)
369	pt_i = 0.5_wp * ( pt(k,j,i) + a * pt(k,j,i-1) + b * pt(k,j-1,i) &
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
381	IF ( rifs >= 0.0_wp ) THEN
382
383	!
384	!-- Stable stratification (and neutral)
385	us_wall = kappa * vel_total / ( LOG( zp / z0(j,i) ) + &
386	5.0_wp * rifs * ( zp - z0(j,i) ) / zp &
387	)
388	ELSE
389
390	!
391	!-- Unstable stratification
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 ) )
394
395	us_wall = kappa * vel_total / ( &
396	LOG( zp / z0(j,i) ) - &
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 ) ) &
400	)
401	ENDIF
402
403	!
404	!-- (4) Compute zp/L (corresponds to neutral Richardson flux number
405	!-- rifs)
406	rifs = -1.0_wp * zp * kappa * g * wspts / &
407	( pt_i * (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), 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).
415	IF ( rifs < zeta_min ) rifs = zeta_min
416	IF ( rifs > zeta_max ) rifs = zeta_max
417
418	!
419	!-- (5) Compute wall_flux (u'v', v'u', w'v', or w'u')
420	IF ( rifs >= 0.0_wp ) THEN
421
422	!
423	!-- Stable stratification (and neutral)
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) ) + &
427	5.0_wp * rifs * ( zp - z0(j,i) ) / zp &
428	)
429	ELSE
430
431	!
432	!-- Unstable stratification
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 ) )
435
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) ) - &
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 ) ) &
442	)
443	ENDIF
444	wall_flux(k) = -wall_flux(k) * us_wall
445
446	!
447	!-- store rifs for next time step
448	rif_wall(k,j,i,wall_index) = rifs
449
450	ENDDO
451
452	END SUBROUTINE wall_fluxes_ij
453
454
455
456	!------------------------------------------------------------------------------!
457	! Description:
458	! ------------
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).
463	!------------------------------------------------------------------------------!
464
465	SUBROUTINE wall_fluxes_e( wall_flux, a, b, c1, c2, wall )
466
467
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
482
483	IMPLICIT NONE
484
485	INTEGER(iwp) :: i !<
486	INTEGER(iwp) :: j !<
487	INTEGER(iwp) :: k !<
488	INTEGER(iwp) :: kk !<
489	INTEGER(iwp) :: wall_index !<
490
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 !<
505
506	REAL(wp), &
507	DIMENSION(nysg:nyng,nxlg:nxrg) :: &
508	wall !<
509
510	REAL(wp), &
511	DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) :: &
512	wall_flux !<
513
514
515	zp = 0.5_wp * ( (a+c1) * dy + (b+c2) * dx )
516	wall_flux = 0.0_wp
517	wall_index = NINT( a+ 2b + 3c1 + 4*c2 )
518
519	DO i = nxl, nxr
520	DO j = nys, nyn
521
522	IF ( wall(j,i) /= 0.0_wp ) THEN
523	!
524	!-- All subsequent variables are computed for scalar locations.
525	DO k = nzb_diff_s_inner(j,i)-1, nzb_diff_s_outer(j,i)-2
526	!
527	!-- (1) Compute rifs, u_i, v_i, and ws
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
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	)
539
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) )
543	!
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 )
547	vel_zp = 0.5_wp * ( a * u_i + b * v_i + (c1+c2) * ws )
548	!
549	!-- (3) Compute wall friction velocity us_wall
550	IF ( rifs >= 0.0_wp ) THEN
551
552	!
553	!-- Stable stratification (and neutral)
554	us_wall = kappa * vel_total / ( LOG( zp / z0(j,i) ) + &
555	5.0_wp * rifs * ( zp - z0(j,i) ) / zp &
556	)
557	ELSE
558
559	!
560	!-- Unstable stratification
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 ) )
563
564	us_wall = kappa * vel_total / ( &
565	LOG( zp / z0(j,i) ) - &
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 ) ) &
569	)
570	ENDIF
571
572	!
573	!-- Skip step (4) of wall_fluxes, because here rifs is already
574	!-- available from (1)
575	!
576	!-- (5) Compute wall_flux (u'v', v'u', w'v', or w'u')
577
578	IF ( rifs >= 0.0_wp ) THEN
579
580	!
581	!-- Stable stratification (and neutral)
582	wall_flux(k,j,i) = kappa * vel_zp / ( LOG( zp/z0(j,i) ) +&
583	5.0_wp * rifs * ( zp-z0(j,i) ) / zp )
584	ELSE
585
586	!
587	!-- Unstable stratification
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 ) )
590
591	wall_flux(k,j,i) = kappa * vel_zp / ( &
592	LOG( zp / z0(j,i) ) - &
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 ) ) &
596	)
597	ENDIF
598	wall_flux(k,j,i) = - wall_flux(k,j,i) * us_wall
599
600	ENDDO
601
602	ENDIF
603
604	ENDDO
605	ENDDO
606
607	END SUBROUTINE wall_fluxes_e
608
609
610	!------------------------------------------------------------------------------!
611	! Description:
612	! ------------
613	!> Call for grid point i,j
614	!------------------------------------------------------------------------------!
615	SUBROUTINE wall_fluxes_e_ij( i, j, nzb_w, nzt_w, wall_flux, a, b, c1, c2 )
616
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
630
631	IMPLICIT NONE
632
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 !<
640
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 !<
655
656	REAL(wp), DIMENSION(nzb:nzt+1) :: wall_flux !<
657
658
659	zp = 0.5_wp * ( (a+c1) * dy + (b+c2) * dx )
660	wall_flux = 0.0_wp
661	wall_index = NINT( a+ 2b + 3c1 + 4*c2 )
662
663	!
664	!-- All subsequent variables are computed for scalar locations.
665	DO k = nzb_w, nzt_w
666
667	!
668	!-- (1) Compute rifs, u_i, v_i, and ws
669	IF ( k == nzb_w ) THEN
670	kk = nzb_w
671	ELSE
672	kk = k-1
673	ENDIF
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	)
680
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) )
684	!
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 )
688	vel_zp = 0.5_wp * ( a * u_i + b * v_i + (c1+c2) * ws )
689	!
690	!-- (3) Compute wall friction velocity us_wall
691	IF ( rifs >= 0.0_wp ) THEN
692
693	!
694	!-- Stable stratification (and neutral)
695	us_wall = kappa * vel_total / ( LOG( zp / z0(j,i) ) + &
696	5.0_wp * rifs * ( zp - z0(j,i) ) / zp &
697	)
698	ELSE
699
700	!
701	!-- Unstable stratification
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 ) )
704
705	us_wall = kappa * vel_total / ( &
706	LOG( zp / z0(j,i) ) - &
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 ) ) &
710	)
711	ENDIF
712
713	!
714	!-- Skip step (4) of wall_fluxes, because here rifs is already
715	!-- available from (1)
716	!
717	!-- (5) Compute wall_flux (u'v', v'u', w'v', or w'u')
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).
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	)
725
726	IF ( rifs >= 0.0_wp ) THEN
727
728	!
729	!-- Stable stratification (and neutral)
730	wall_flux(k) = kappa * vel_zp / &
731	( LOG( zp/z0(j,i) ) + 5.0_wp * rifs * ( zp-z0(j,i) ) / zp )
732	ELSE
733
734	!
735	!-- Unstable stratification
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 ) )
738
739	wall_flux(k) = kappa * vel_zp / ( &
740	LOG( zp / z0(j,i) ) - &
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	)
745	ENDIF
746	wall_flux(k) = - wall_flux(k) * us_wall
747
748	ENDDO
749
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 |