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

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 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	!
16	! Copyright 1997-2014 Leibniz Universitaet Hannover
17	!--------------------------------------------------------------------------------!
18	!
19	! Current revisions:
20	! -----------------
21	!
22	!
23	! Former revisions:
24	! -----------------
25	! $Id: wall_fluxes.f90 1683 2015-10-07 23:57:51Z knoop $
26	!
27	! 1682 2015-10-07 23:56:08Z knoop
28	! Code annotations made doxygen readable
29	!
30	! 1374 2014-04-25 12:55:07Z raasch
31	! pt removed from acc-present-list
32	!
33	! 1353 2014-04-08 15:21:23Z heinze
34	! REAL constants provided with KIND-attribute
35	!
36	! 1320 2014-03-20 08:40:49Z raasch
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
44	!
45	! 1257 2013-11-08 15:18:40Z raasch
46	! openacc loop and loop vector clauses removed
47	!
48	! 1153 2013-05-10 14:33:08Z raasch
49	! code adjustments of accelerator version required by PGI 12.3 / CUDA 5.0
50	!
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	!
55	! 1036 2012-10-22 13:43:42Z raasch
56	! code put under GPL (PALM 3.9)
57	!
58	! 1015 2012-09-27 09:23:24Z raasch
59	! accelerator version (*_acc) added
60	!
61	! Initial version (2007/03/07)
62	!
63	! Description:
64	! ------------
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.
71	!------------------------------------------------------------------------------!
72	MODULE wall_fluxes_mod
73
74	PRIVATE
75	PUBLIC wall_fluxes, wall_fluxes_acc, wall_fluxes_e, wall_fluxes_e_acc
76
77	INTERFACE wall_fluxes
78	MODULE PROCEDURE wall_fluxes
79	MODULE PROCEDURE wall_fluxes_ij
80	END INTERFACE wall_fluxes
81
82	INTERFACE wall_fluxes_acc
83	MODULE PROCEDURE wall_fluxes_acc
84	END INTERFACE wall_fluxes_acc
85
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
91	INTERFACE wall_fluxes_e_acc
92	MODULE PROCEDURE wall_fluxes_e_acc
93	END INTERFACE wall_fluxes_e_acc
94
95	CONTAINS
96
97	!------------------------------------------------------------------------------!
98	! Description:
99	! ------------
100	!> Call for all grid points
101	!------------------------------------------------------------------------------!
102	SUBROUTINE wall_fluxes( wall_flux, a, b, c1, c2, nzb_uvw_inner, &
103	nzb_uvw_outer, wall )
104
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
121
122	IMPLICIT NONE
123
124	INTEGER(iwp) :: i !<
125	INTEGER(iwp) :: j !<
126	INTEGER(iwp) :: k !<
127	INTEGER(iwp) :: wall_index !<
128
129	INTEGER(iwp), &
130	DIMENSION(nysg:nyng,nxlg:nxrg) :: &
131	nzb_uvw_inner !<
132	INTEGER(iwp), &
133	DIMENSION(nysg:nyng,nxlg:nxrg) :: &
134	nzb_uvw_outer !<
135
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 !<
152
153	REAL(wp), &
154	DIMENSION(nysg:nyng,nxlg:nxrg) :: &
155	wall !<
156
157	REAL(wp), &
158	DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) :: &
159	wall_flux !<
160
161
162	zp = 0.5_wp * ( (a+c1) * dy + (b+c2) * dx )
163	wall_flux = 0.0_wp
164	wall_index = NINT( a+ 2b + 3c1 + 4*c2 )
165
166	DO i = nxl, nxr
167	DO j = nys, nyn
168
169	IF ( wall(j,i) /= 0.0_wp ) THEN
170	!
171	!-- All subsequent variables are computed for the respective
172	!-- location where the respective flux is defined.
173	DO k = nzb_uvw_inner(j,i)+1, nzb_uvw_outer(j,i)
174
175	!
176	!-- (1) Compute rifs, u_i, v_i, ws, pt' and w'pt'
177	rifs = rif_wall(k,j,i,wall_index)
178
179	u_i = a * u(k,j,i) + c1 * 0.25_wp * &
180	( u(k+1,j,i+1) + u(k+1,j,i) + u(k,j,i+1) + u(k,j,i) )
181
182	v_i = b * v(k,j,i) + c2 * 0.25_wp * &
183	( v(k+1,j+1,i) + v(k+1,j,i) + v(k,j+1,i) + v(k,j,i) )
184
185	ws = ( c1 + c2 ) * w(k,j,i) + 0.25_wp * ( &
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) ) &
188	)
189	pt_i = 0.5_wp * ( pt(k,j,i) + a * pt(k,j,i-1) + &
190	b * pt(k,j-1,i) + ( c1 + c2 ) * pt(k+1,j,i) )
191
192	pts = pt_i - hom(k,1,4,0)
193	wspts = ws * pts
194
195	!
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 )
198
199	!
200	!-- (3) Compute wall friction velocity us_wall
201	IF ( rifs >= 0.0_wp ) THEN
202
203	!
204	!-- Stable stratification (and neutral)
205	us_wall = kappa * vel_total / ( LOG( zp / z0(j,i) ) + &
206	5.0_wp * rifs * ( zp - z0(j,i) ) / zp &
207	)
208	ELSE
209
210	!
211	!-- Unstable stratification
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 ) )
214
215	us_wall = kappa * vel_total / ( &
216	LOG( zp / z0(j,i) ) - &
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 ) ) &
220	)
221	ENDIF
222
223	!
224	!-- (4) Compute zp/L (corresponds to neutral Richardson flux
225	!-- number rifs)
226	rifs = -1.0_wp * zp * kappa * g * wspts / &
227	( pt_i * ( us_wall**3 + 1E-30 ) )
228
229	!
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
238
239	!
240	!-- (5) Compute wall_flux (u'v', v'u', w'v', or w'u')
241	IF ( rifs >= 0.0_wp ) THEN
242
243	!
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) ) + &
248	5.0_wp * rifs * ( zp - z0(j,i) ) / zp &
249	)
250	ELSE
251
252	!
253	!-- Unstable stratification
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 ) )
256
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	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 ) ) &
263	)
264	ENDIF
265	wall_flux(k,j,i) = -wall_flux(k,j,i) * us_wall
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
281	!------------------------------------------------------------------------------!
282	! Description:
283	! ------------
284	!> Call for all grid points - accelerator version
285	!------------------------------------------------------------------------------!
286	SUBROUTINE wall_fluxes_acc( wall_flux, a, b, c1, c2, nzb_uvw_inner, &
287	nzb_uvw_outer, wall )
288
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
306
307	IMPLICIT NONE
308
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 !<
315
316	INTEGER(iwp), &
317	DIMENSION(nysg:nyng,nxlg:nxrg) :: &
318	nzb_uvw_inner !<
319	INTEGER(iwp), &
320	DIMENSION(nysg:nyng,nxlg:nxrg) :: &
321	nzb_uvw_outer !<
322
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 !<
339
340	REAL(wp), &
341	DIMENSION(nysg:nyng,nxlg:nxrg) :: &
342	wall !<
343
344	REAL(wp), &
345	DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) :: &
346	wall_flux !<
347
348
349	zp = 0.5_wp * ( (a+c1) * dy + (b+c2) * dx )
350	wall_flux = 0.0_wp
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 )
358	!$acc loop independent
359	DO i = i_left, i_right
360	DO j = j_south, j_north
361
362	IF ( wall(j,i) /= 0.0_wp ) THEN
363	!
364	!-- All subsequent variables are computed for the respective
365	!-- location where the respective flux is defined.
366	!$acc loop independent
367	DO k = nzb_uvw_inner(j,i)+1, nzb_uvw_outer(j,i)
368
369	!
370	!-- (1) Compute rifs, u_i, v_i, ws, pt' and w'pt'
371	rifs = rif_wall(k,j,i,wall_index)
372
373	u_i = a * u(k,j,i) + c1 * 0.25_wp * &
374	( u(k+1,j,i+1) + u(k+1,j,i) + u(k,j,i+1) + u(k,j,i) )
375
376	v_i = b * v(k,j,i) + c2 * 0.25_wp * &
377	( v(k+1,j+1,i) + v(k+1,j,i) + v(k,j+1,i) + v(k,j,i) )
378
379	ws = ( c1 + c2 ) * w(k,j,i) + 0.25_wp * ( &
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) ) &
382	)
383	pt_i = 0.5_wp * ( pt(k,j,i) + a * pt(k,j,i-1) + &
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
395	IF ( rifs >= 0.0_wp ) THEN
396
397	!
398	!-- Stable stratification (and neutral)
399	us_wall = kappa * vel_total / ( LOG( zp / z0(j,i) ) + &
400	5.0_wp * rifs * ( zp - z0(j,i) ) / zp &
401	)
402	ELSE
403
404	!
405	!-- Unstable stratification
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 ) )
408
409	us_wall = kappa * vel_total / ( &
410	LOG( zp / z0(j,i) ) - &
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 ) ) &
414	)
415	ENDIF
416
417	!
418	!-- (4) Compute zp/L (corresponds to neutral Richardson flux
419	!-- number rifs)
420	rifs = -1.0_wp * zp * kappa * g * wspts / &
421	( pt_i * ( us_wall**3 + 1E-30 ) )
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')
435	IF ( rifs >= 0.0_wp ) THEN
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) ) + &
442	5.0_wp * rifs * ( zp - z0(j,i) ) / zp &
443	)
444	ELSE
445
446	!
447	!-- Unstable stratification
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 ) )
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) ) - &
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 ) ) &
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
465	ENDDO
466
467	ENDIF
468
469	ENDDO
470	ENDDO
471	!$acc end kernels
472
473	END SUBROUTINE wall_fluxes_acc
474
475
476	!------------------------------------------------------------------------------!
477	! Description:
478	! ------------
479	!> Call for all grid point i,j
480	!------------------------------------------------------------------------------!
481	SUBROUTINE wall_fluxes_ij( i, j, nzb_w, nzt_w, wall_flux, a, b, c1, c2 )
482
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
499
500	IMPLICIT NONE
501
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 !<
508
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 !<
525
526	REAL(wp), DIMENSION(nzb:nzt+1) :: wall_flux !<
527
528
529	zp = 0.5_wp * ( (a+c1) * dy + (b+c2) * dx )
530	wall_flux = 0.0_wp
531	wall_index = NINT( a+ 2b + 3c1 + 4*c2 )
532
533	!
534	!-- All subsequent variables are computed for the respective location where
535	!-- the respective flux is defined.
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
542	u_i = a * u(k,j,i) + c1 * 0.25_wp * &
543	( u(k+1,j,i+1) + u(k+1,j,i) + u(k,j,i+1) + u(k,j,i) )
544
545	v_i = b * v(k,j,i) + c2 * 0.25_wp * &
546	( v(k+1,j+1,i) + v(k+1,j,i) + v(k,j+1,i) + v(k,j,i) )
547
548	ws = ( c1 + c2 ) * w(k,j,i) + 0.25_wp * ( &
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) ) &
551	)
552	pt_i = 0.5_wp * ( pt(k,j,i) + a * pt(k,j,i-1) + b * pt(k,j-1,i) &
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
564	IF ( rifs >= 0.0_wp ) THEN
565
566	!
567	!-- Stable stratification (and neutral)
568	us_wall = kappa * vel_total / ( LOG( zp / z0(j,i) ) + &
569	5.0_wp * rifs * ( zp - z0(j,i) ) / zp &
570	)
571	ELSE
572
573	!
574	!-- Unstable stratification
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 ) )
577
578	us_wall = kappa * vel_total / ( &
579	LOG( zp / z0(j,i) ) - &
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 ) ) &
583	)
584	ENDIF
585
586	!
587	!-- (4) Compute zp/L (corresponds to neutral Richardson flux number
588	!-- rifs)
589	rifs = -1.0_wp * zp * kappa * g * wspts / &
590	( pt_i * (us_wall**3 + 1E-30) )
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')
603	IF ( rifs >= 0.0_wp ) THEN
604
605	!
606	!-- Stable stratification (and neutral)
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) ) + &
610	5.0_wp * rifs * ( zp - z0(j,i) ) / zp &
611	)
612	ELSE
613
614	!
615	!-- Unstable stratification
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 ) )
618
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) ) - &
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 ) ) &
625	)
626	ENDIF
627	wall_flux(k) = -wall_flux(k) * us_wall
628
629	!
630	!-- store rifs for next time step
631	rif_wall(k,j,i,wall_index) = rifs
632
633	ENDDO
634
635	END SUBROUTINE wall_fluxes_ij
636
637
638
639	!------------------------------------------------------------------------------!
640	! Description:
641	! ------------
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).
646	!------------------------------------------------------------------------------!
647
648	SUBROUTINE wall_fluxes_e( wall_flux, a, b, c1, c2, wall )
649
650
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
665
666	IMPLICIT NONE
667
668	INTEGER(iwp) :: i !<
669	INTEGER(iwp) :: j !<
670	INTEGER(iwp) :: k !<
671	INTEGER(iwp) :: kk !<
672	INTEGER(iwp) :: wall_index !<
673
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 !<
688
689	REAL(wp), &
690	DIMENSION(nysg:nyng,nxlg:nxrg) :: &
691	wall !<
692
693	REAL(wp), &
694	DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) :: &
695	wall_flux !<
696
697
698	zp = 0.5_wp * ( (a+c1) * dy + (b+c2) * dx )
699	wall_flux = 0.0_wp
700	wall_index = NINT( a+ 2b + 3c1 + 4*c2 )
701
702	DO i = nxl, nxr
703	DO j = nys, nyn
704
705	IF ( wall(j,i) /= 0.0_wp ) THEN
706	!
707	!-- All subsequent variables are computed for scalar locations.
708	DO k = nzb_diff_s_inner(j,i)-1, nzb_diff_s_outer(j,i)-2
709	!
710	!-- (1) Compute rifs, u_i, v_i, and ws
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
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	)
722
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) )
726	!
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 )
730	vel_zp = 0.5_wp * ( a * u_i + b * v_i + (c1+c2) * ws )
731	!
732	!-- (3) Compute wall friction velocity us_wall
733	IF ( rifs >= 0.0_wp ) THEN
734
735	!
736	!-- Stable stratification (and neutral)
737	us_wall = kappa * vel_total / ( LOG( zp / z0(j,i) ) + &
738	5.0_wp * rifs * ( zp - z0(j,i) ) / zp &
739	)
740	ELSE
741
742	!
743	!-- Unstable stratification
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 ) )
746
747	us_wall = kappa * vel_total / ( &
748	LOG( zp / z0(j,i) ) - &
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 ) ) &
752	)
753	ENDIF
754
755	!
756	!-- Skip step (4) of wall_fluxes, because here rifs is already
757	!-- available from (1)
758	!
759	!-- (5) Compute wall_flux (u'v', v'u', w'v', or w'u')
760
761	IF ( rifs >= 0.0_wp ) THEN
762
763	!
764	!-- Stable stratification (and neutral)
765	wall_flux(k,j,i) = kappa * vel_zp / ( LOG( zp/z0(j,i) ) +&
766	5.0_wp * rifs * ( zp-z0(j,i) ) / zp )
767	ELSE
768
769	!
770	!-- Unstable stratification
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 ) )
773
774	wall_flux(k,j,i) = kappa * vel_zp / ( &
775	LOG( zp / z0(j,i) ) - &
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 ) ) &
779	)
780	ENDIF
781	wall_flux(k,j,i) = - wall_flux(k,j,i) * us_wall
782
783	ENDDO
784
785	ENDIF
786
787	ENDDO
788	ENDDO
789
790	END SUBROUTINE wall_fluxes_e
791
792
793	!------------------------------------------------------------------------------!
794	! Description:
795	! ------------
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).
800	!------------------------------------------------------------------------------!
801	SUBROUTINE wall_fluxes_e_acc( wall_flux, a, b, c1, c2, wall )
802
803
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
819
820	IMPLICIT NONE
821
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 !<
829
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 !<
844
845	REAL(wp), &
846	DIMENSION(nysg:nyng,nxlg:nxrg) :: &
847	wall !<
848
849	REAL(wp), &
850	DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) :: &
851	wall_flux !<
852
853
854	zp = 0.5_wp * ( (a+c1) * dy + (b+c2) * dx )
855	wall_flux = 0.0_wp
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
861	!$acc kernels present( nzb_diff_s_inner, nzb_diff_s_outer, rif_wall ) &
862	!$acc present( u, v, w, wall, wall_flux, z0 )
863	DO i = i_left, i_right
864	DO j = j_south, j_north
865	DO k = min_inner, max_outer
866	!
867	!-- All subsequent variables are computed for scalar locations
868	IF ( k >= nzb_diff_s_inner(j,i)-1 .AND. &
869	k <= nzb_diff_s_outer(j,i)-2 .AND. &
870	wall(j,i) /= 0.0_wp ) THEN
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
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	)
884
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) )
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 )
892	vel_zp = 0.5_wp * ( a * u_i + b * v_i + (c1+c2) * ws )
893	!
894	!-- (3) Compute wall friction velocity us_wall
895	IF ( rifs >= 0.0_wp ) THEN
896
897	!
898	!-- Stable stratification (and neutral)
899	us_wall = kappa * vel_total / ( LOG( zp / z0(j,i) ) + &
900	5.0_wp * rifs * ( zp - z0(j,i) ) / zp &
901	)
902	ELSE
903
904	!
905	!-- Unstable stratification
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 ) )
908
909	us_wall = kappa * vel_total / ( &
910	LOG( zp / z0(j,i) ) - &
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 ) ) &
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
923	IF ( rifs >= 0.0_wp ) THEN
924
925	!
926	!-- Stable stratification (and neutral)
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	)
931	ELSE
932
933	!
934	!-- Unstable stratification
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 ) )
937
938	wall_flux(k,j,i) = kappa * vel_zp / ( &
939	LOG( zp / z0(j,i) ) - &
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 ) ) &
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	!------------------------------------------------------------------------------!
958	! Description:
959	! ------------
960	!> Call for grid point i,j
961	!------------------------------------------------------------------------------!
962	SUBROUTINE wall_fluxes_e_ij( i, j, nzb_w, nzt_w, wall_flux, a, b, c1, c2 )
963
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
977
978	IMPLICIT NONE
979
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 !<
987
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 !<
1002
1003	REAL(wp), DIMENSION(nzb:nzt+1) :: wall_flux !<
1004
1005
1006	zp = 0.5_wp * ( (a+c1) * dy + (b+c2) * dx )
1007	wall_flux = 0.0_wp
1008	wall_index = NINT( a+ 2b + 3c1 + 4*c2 )
1009
1010	!
1011	!-- All subsequent variables are computed for scalar locations.
1012	DO k = nzb_w, nzt_w
1013
1014	!
1015	!-- (1) Compute rifs, u_i, v_i, and ws
1016	IF ( k == nzb_w ) THEN
1017	kk = nzb_w
1018	ELSE
1019	kk = k-1
1020	ENDIF
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	)
1027
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) )
1031	!
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 )
1035	vel_zp = 0.5_wp * ( a * u_i + b * v_i + (c1+c2) * ws )
1036	!
1037	!-- (3) Compute wall friction velocity us_wall
1038	IF ( rifs >= 0.0_wp ) THEN
1039
1040	!
1041	!-- Stable stratification (and neutral)
1042	us_wall = kappa * vel_total / ( LOG( zp / z0(j,i) ) + &
1043	5.0_wp * rifs * ( zp - z0(j,i) ) / zp &
1044	)
1045	ELSE
1046
1047	!
1048	!-- Unstable stratification
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 ) )
1051
1052	us_wall = kappa * vel_total / ( &
1053	LOG( zp / z0(j,i) ) - &
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 ) ) &
1057	)
1058	ENDIF
1059
1060	!
1061	!-- Skip step (4) of wall_fluxes, because here rifs is already
1062	!-- available from (1)
1063	!
1064	!-- (5) Compute wall_flux (u'v', v'u', w'v', or w'u')
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).
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	)
1072
1073	IF ( rifs >= 0.0_wp ) THEN
1074
1075	!
1076	!-- Stable stratification (and neutral)
1077	wall_flux(k) = kappa * vel_zp / &
1078	( LOG( zp/z0(j,i) ) + 5.0_wp * rifs * ( zp-z0(j,i) ) / zp )
1079	ELSE
1080
1081	!
1082	!-- Unstable stratification
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 ) )
1085
1086	wall_flux(k) = kappa * vel_zp / ( &
1087	LOG( zp / z0(j,i) ) - &
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	)
1092	ENDIF
1093	wall_flux(k) = - wall_flux(k) * us_wall
1094
1095	ENDDO
1096
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 |