Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |