Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |