Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |