Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |