Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |