Home

Context Navigation

source: palm/trunk/SOURCE/wall_fluxes_mod.f90 @ 1850

Last change on this file since 1850 was 1850, checked in by maronga, 8 years ago
added _mod string to several filenames to meet the naming convection for modules
Property svn:keywords set to `Id`
File size: 43.3 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |