Home

Context Navigation

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

Last change on this file since 1357 was 1354, checked in by heinze, 11 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 43.1 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |