Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |