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

Last change on this file since 1128 was 1128, checked in by raasch, 8 years ago

asynchronous transfer of ghost point data for acc-optimized version

  • Property svn:keywords set to Id
File size: 33.4 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-2012  Leibniz University Hannover
18!--------------------------------------------------------------------------------!
19!
20! Current revisions:
21! -----------------
22! loop index bounds in accelerator version replaced by i_left, i_right, j_south,
23! j_north
24!
25! Former revisions:
26! -----------------
27! $Id: wall_fluxes.f90 1128 2013-04-12 06:19:32Z raasch $
28!
29! 1036 2012-10-22 13:43:42Z raasch
30! code put under GPL (PALM 3.9)
31!
32! 1015 2012-09-27 09:23:24Z raasch
33! accelerator version (*_acc) added
34!
35! 187 2008-08-06 16:25:09Z letzel
36! Bugfix: Modification of the evaluation of the vertical turbulent momentum
37! fluxes u'w' and v'w (see prandtl_fluxes), this requires the calculation of
38! us_wall (and vel_total, u_i, v_i, ws) also in wall_fluxes_e.
39! Bugfix: change definition of us_wall from 1D to 2D
40! Bugfix: storage of rifs to rifs_wall in wall_fluxes_e removed
41! Change: add 'minus' sign to fluxes produced by subroutine wall_fluxes_e for
42! consistency with subroutine wall_fluxes
43! Change: Modification of the integrated version of the profile function for
44! momentum for unstable stratification
45!
46! Initial version (2007/03/07)
47!
48! Description:
49! ------------
50! Calculates momentum fluxes at vertical walls assuming Monin-Obukhov
51! similarity.
52! Indices: usvs a=1, vsus b=1, wsvs c1=1, wsus c2=1 (other=0).
53! The all-gridpoint version of wall_fluxes_e is not used so far, because
54! it gives slightly different results from the ij-version for some unknown
55! reason.
56!------------------------------------------------------------------------------!
57    PRIVATE
58    PUBLIC wall_fluxes, wall_fluxes_acc, wall_fluxes_e, wall_fluxes_e_acc
59   
60    INTERFACE wall_fluxes
61       MODULE PROCEDURE wall_fluxes
62       MODULE PROCEDURE wall_fluxes_ij
63    END INTERFACE wall_fluxes
64   
65    INTERFACE wall_fluxes_acc
66       MODULE PROCEDURE wall_fluxes_acc
67    END INTERFACE wall_fluxes_acc
68
69    INTERFACE wall_fluxes_e
70       MODULE PROCEDURE wall_fluxes_e
71       MODULE PROCEDURE wall_fluxes_e_ij
72    END INTERFACE wall_fluxes_e
73 
74    INTERFACE wall_fluxes_e_acc
75       MODULE PROCEDURE wall_fluxes_e_acc
76    END INTERFACE wall_fluxes_e_acc
77
78 CONTAINS
79
80!------------------------------------------------------------------------------!
81! Call for all grid points
82!------------------------------------------------------------------------------!
83    SUBROUTINE wall_fluxes( wall_flux, a, b, c1, c2, nzb_uvw_inner, &
84                            nzb_uvw_outer, wall )
85
86       USE arrays_3d
87       USE control_parameters
88       USE grid_variables
89       USE indices
90       USE statistics
91
92       IMPLICIT NONE
93
94       INTEGER ::  i, j, k, wall_index
95
96       INTEGER, DIMENSION(nysg:nyng,nxlg:nxrg) ::  nzb_uvw_inner, &
97                                                       nzb_uvw_outer
98       REAL ::  a, b, c1, c2, h1, h2, zp
99       REAL ::  pts, pt_i, rifs, u_i, v_i, us_wall, vel_total, ws, wspts
100
101       REAL, DIMENSION(nysg:nyng,nxlg:nxrg)   ::  wall
102       REAL, DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) ::  wall_flux
103
104
105       zp         = 0.5 * ( (a+c1) * dy + (b+c2) * dx )
106       wall_flux  = 0.0
107       wall_index = NINT( a+ 2*b + 3*c1 + 4*c2 )
108
109       DO  i = nxl, nxr
110          DO  j = nys, nyn
111
112             IF ( wall(j,i) /= 0.0 )  THEN
113!
114!--             All subsequent variables are computed for the respective
115!--             location where the respective flux is defined.
116                DO  k = nzb_uvw_inner(j,i)+1, nzb_uvw_outer(j,i)
117
118!
119!--                (1) Compute rifs, u_i, v_i, ws, pt' and w'pt'
120                   rifs  = rif_wall(k,j,i,wall_index)
121
122                   u_i   = a * u(k,j,i) + c1 * 0.25 * &
123                           ( u(k+1,j,i+1) + u(k+1,j,i) + u(k,j,i+1) + u(k,j,i) )
124
125                   v_i   = b * v(k,j,i) + c2 * 0.25 * &
126                           ( v(k+1,j+1,i) + v(k+1,j,i) + v(k,j+1,i) + v(k,j,i) )
127
128                   ws    = ( c1 + c2 ) * w(k,j,i) + 0.25 * (                   &
129                     a * ( w(k-1,j,i-1) + w(k-1,j,i) + w(k,j,i-1) + w(k,j,i) ) &
130                   + b * ( w(k-1,j-1,i) + w(k-1,j,i) + w(k,j-1,i) + w(k,j,i) ) &
131                                                           )
132                   pt_i  = 0.5 * ( pt(k,j,i) + a *  pt(k,j,i-1) + &
133                                   b * pt(k,j-1,i) + ( c1 + c2 ) * pt(k+1,j,i) )
134
135                   pts   = pt_i - hom(k,1,4,0)
136                   wspts = ws * pts
137
138!
139!--                (2) Compute wall-parallel absolute velocity vel_total
140                   vel_total = SQRT( ws**2 + (a+c1) * u_i**2 + (b+c2) * v_i**2 )
141
142!
143!--                (3) Compute wall friction velocity us_wall
144                   IF ( rifs >= 0.0 )  THEN
145
146!
147!--                   Stable stratification (and neutral)
148                      us_wall = kappa * vel_total / ( LOG( zp / z0(j,i) ) +    &
149                                            5.0 * rifs * ( zp - z0(j,i) ) / zp &
150                                                    )
151                   ELSE
152
153!
154!--                   Unstable stratification
155                      h1 = SQRT( SQRT( 1.0 - 16.0 * rifs ) )
156                      h2 = SQRT( SQRT( 1.0 - 16.0 * rifs * z0(j,i) / zp ) )
157
158                      us_wall = kappa * vel_total / (                          &
159                           LOG( zp / z0(j,i) ) -                               &
160                           LOG( ( 1.0 + h1 )**2 * ( 1.0 + h1**2 ) / (          &
161                                ( 1.0 + h2 )**2 * ( 1.0 + h2**2 )   ) ) +      &
162                                2.0 * ( ATAN( h1 ) - ATAN( h2 ) )              &
163                                                    )
164                   ENDIF
165
166!
167!--                (4) Compute zp/L (corresponds to neutral Richardson flux
168!--                    number rifs)
169                   rifs = -1.0 * zp * kappa * g * wspts / ( pt_i * &
170                                                        ( us_wall**3 + 1E-30 ) )
171
172!
173!--                Limit the value range of the Richardson numbers.
174!--                This is necessary for very small velocities (u,w --> 0),
175!--                because the absolute value of rif can then become very
176!--                large, which in consequence would result in very large
177!--                shear stresses and very small momentum fluxes (both are
178!--                generally unrealistic).
179                   IF ( rifs < rif_min )  rifs = rif_min
180                   IF ( rifs > rif_max )  rifs = rif_max
181
182!
183!--                (5) Compute wall_flux (u'v', v'u', w'v', or w'u')
184                   IF ( rifs >= 0.0 )  THEN
185
186!
187!--                   Stable stratification (and neutral)
188                      wall_flux(k,j,i) = kappa *                               &
189                              ( a*u(k,j,i) + b*v(k,j,i) + (c1+c2)*w(k,j,i) ) / &
190                              (  LOG( zp / z0(j,i) ) +                         &
191                                 5.0 * rifs * ( zp - z0(j,i) ) / zp            &
192                              )
193                   ELSE
194
195!
196!--                   Unstable stratification
197                      h1 = SQRT( SQRT( 1.0 - 16.0 * rifs ) )
198                      h2 = SQRT( SQRT( 1.0 - 16.0 * rifs * z0(j,i) / zp ) )
199
200                      wall_flux(k,j,i) = kappa *                               &
201                           ( a*u(k,j,i) + b*v(k,j,i) + (c1+c2)*w(k,j,i) ) / (  &
202                           LOG( zp / z0(j,i) ) -                               &
203                           LOG( ( 1.0 + h1 )**2 * ( 1.0 + h1**2 ) / (          &
204                                ( 1.0 + h2 )**2 * ( 1.0 + h2**2 )   ) ) +      &
205                                2.0 * ( ATAN( h1 ) - ATAN( h2 ) )              &
206                                                                            )
207                   ENDIF
208                   wall_flux(k,j,i) = -wall_flux(k,j,i) * us_wall
209
210!
211!--                store rifs for next time step
212                   rif_wall(k,j,i,wall_index) = rifs
213
214                ENDDO
215
216             ENDIF
217
218          ENDDO
219       ENDDO
220
221    END SUBROUTINE wall_fluxes
222
223
224!------------------------------------------------------------------------------!
225! Call for all grid points - accelerator version
226!------------------------------------------------------------------------------!
227    SUBROUTINE wall_fluxes_acc( wall_flux, a, b, c1, c2, nzb_uvw_inner, &
228                                nzb_uvw_outer, wall )
229
230       USE arrays_3d
231       USE control_parameters
232       USE grid_variables
233       USE indices
234       USE statistics
235
236       IMPLICIT NONE
237
238       INTEGER ::  i, j, k, max_outer, min_inner, wall_index
239
240       INTEGER, DIMENSION(nysg:nyng,nxlg:nxrg) ::  nzb_uvw_inner, &
241                                                   nzb_uvw_outer
242       REAL ::  a, b, c1, c2, h1, h2, zp
243       REAL ::  pts, pt_i, rifs, u_i, v_i, us_wall, vel_total, ws, wspts
244
245       REAL, DIMENSION(nysg:nyng,nxlg:nxrg)   ::  wall
246       REAL, DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) ::  wall_flux
247
248
249       zp         = 0.5 * ( (a+c1) * dy + (b+c2) * dx )
250       wall_flux  = 0.0
251       wall_index = NINT( a+ 2*b + 3*c1 + 4*c2 )
252
253       min_inner = MINVAL( nzb_uvw_inner(nys:nyn,nxl:nxr) ) + 1
254       max_outer = MINVAL( nzb_uvw_outer(nys:nyn,nxl:nxr) )
255
256       !$acc kernels present( hom, nzb_uvw_inner, nzb_uvw_outer, pt, rif_wall ) &
257       !$acc         present( u, v, w, wall, wall_flux, z0 )
258       !$acc loop
259       DO  i = i_left, i_right
260          DO  j = j_south, j_north
261             !$acc loop vector( 32 )
262             DO  k = min_inner, max_outer
263!
264!--             All subsequent variables are computed for the respective
265!--             location where the respective flux is defined.
266                IF ( k >= nzb_uvw_inner(j,i)+1  .AND. &
267                     k <= nzb_uvw_outer(j,i)    .AND.  wall(j,i) /= 0.0 )  THEN
268!
269!--                (1) Compute rifs, u_i, v_i, ws, pt' and w'pt'
270                   rifs  = rif_wall(k,j,i,wall_index)
271
272                   u_i   = a * u(k,j,i) + c1 * 0.25 * &
273                           ( u(k+1,j,i+1) + u(k+1,j,i) + u(k,j,i+1) + u(k,j,i) )
274
275                   v_i   = b * v(k,j,i) + c2 * 0.25 * &
276                           ( v(k+1,j+1,i) + v(k+1,j,i) + v(k,j+1,i) + v(k,j,i) )
277
278                   ws    = ( c1 + c2 ) * w(k,j,i) + 0.25 * (                   &
279                     a * ( w(k-1,j,i-1) + w(k-1,j,i) + w(k,j,i-1) + w(k,j,i) ) &
280                   + b * ( w(k-1,j-1,i) + w(k-1,j,i) + w(k,j-1,i) + w(k,j,i) ) &
281                                                           )
282                   pt_i  = 0.5 * ( pt(k,j,i) + a *  pt(k,j,i-1) + &
283                                   b * pt(k,j-1,i) + ( c1 + c2 ) * pt(k+1,j,i) )
284
285                   pts   = pt_i - hom(k,1,4,0)
286                   wspts = ws * pts
287
288!
289!--                (2) Compute wall-parallel absolute velocity vel_total
290                   vel_total = SQRT( ws**2 + (a+c1) * u_i**2 + (b+c2) * v_i**2 )
291
292!
293!--                (3) Compute wall friction velocity us_wall
294                   IF ( rifs >= 0.0 )  THEN
295
296!
297!--                   Stable stratification (and neutral)
298                      us_wall = kappa * vel_total / ( LOG( zp / z0(j,i) ) +    &
299                                            5.0 * rifs * ( zp - z0(j,i) ) / zp &
300                                                    )
301                   ELSE
302
303!
304!--                   Unstable stratification
305                      h1 = SQRT( SQRT( 1.0 - 16.0 * rifs ) )
306                      h2 = SQRT( SQRT( 1.0 - 16.0 * rifs * z0(j,i) / zp ) )
307
308                      us_wall = kappa * vel_total / (                          &
309                           LOG( zp / z0(j,i) ) -                               &
310                           LOG( ( 1.0 + h1 )**2 * ( 1.0 + h1**2 ) / (          &
311                                ( 1.0 + h2 )**2 * ( 1.0 + h2**2 )   ) ) +      &
312                                2.0 * ( ATAN( h1 ) - ATAN( h2 ) )              &
313                                                    )
314                   ENDIF
315
316!
317!--                (4) Compute zp/L (corresponds to neutral Richardson flux
318!--                    number rifs)
319                   rifs = -1.0 * zp * kappa * g * wspts / ( pt_i * &
320                                                        ( us_wall**3 + 1E-30 ) )
321
322!
323!--                Limit the value range of the Richardson numbers.
324!--                This is necessary for very small velocities (u,w --> 0),
325!--                because the absolute value of rif can then become very
326!--                large, which in consequence would result in very large
327!--                shear stresses and very small momentum fluxes (both are
328!--                generally unrealistic).
329                   IF ( rifs < rif_min )  rifs = rif_min
330                   IF ( rifs > rif_max )  rifs = rif_max
331
332!
333!--                (5) Compute wall_flux (u'v', v'u', w'v', or w'u')
334                   IF ( rifs >= 0.0 )  THEN
335
336!
337!--                   Stable stratification (and neutral)
338                      wall_flux(k,j,i) = kappa *                               &
339                              ( a*u(k,j,i) + b*v(k,j,i) + (c1+c2)*w(k,j,i) ) / &
340                              (  LOG( zp / z0(j,i) ) +                         &
341                                 5.0 * rifs * ( zp - z0(j,i) ) / zp            &
342                              )
343                   ELSE
344
345!
346!--                   Unstable stratification
347                      h1 = SQRT( SQRT( 1.0 - 16.0 * rifs ) )
348                      h2 = SQRT( SQRT( 1.0 - 16.0 * rifs * z0(j,i) / zp ) )
349
350                      wall_flux(k,j,i) = kappa *                               &
351                           ( a*u(k,j,i) + b*v(k,j,i) + (c1+c2)*w(k,j,i) ) / (  &
352                           LOG( zp / z0(j,i) ) -                               &
353                           LOG( ( 1.0 + h1 )**2 * ( 1.0 + h1**2 ) / (          &
354                                ( 1.0 + h2 )**2 * ( 1.0 + h2**2 )   ) ) +      &
355                                2.0 * ( ATAN( h1 ) - ATAN( h2 ) )              &
356                                                                            )
357                   ENDIF
358                   wall_flux(k,j,i) = -wall_flux(k,j,i) * us_wall
359
360!
361!--                store rifs for next time step
362                   rif_wall(k,j,i,wall_index) = rifs
363
364                ENDIF
365
366             ENDDO
367          ENDDO
368       ENDDO
369       !$acc end kernels
370
371    END SUBROUTINE wall_fluxes_acc
372
373
374!------------------------------------------------------------------------------!
375! Call for all grid point i,j
376!------------------------------------------------------------------------------!
377    SUBROUTINE wall_fluxes_ij( i, j, nzb_w, nzt_w, wall_flux, a, b, c1, c2 )
378
379       USE arrays_3d
380       USE control_parameters
381       USE grid_variables
382       USE indices
383       USE statistics
384
385       IMPLICIT NONE
386
387       INTEGER ::  i, j, k, nzb_w, nzt_w, wall_index
388       REAL    ::  a, b, c1, c2, h1, h2, zp
389
390       REAL ::  pts, pt_i, rifs, u_i, v_i, us_wall, vel_total, ws, wspts
391
392       REAL, DIMENSION(nzb:nzt+1) ::  wall_flux
393
394
395       zp         = 0.5 * ( (a+c1) * dy + (b+c2) * dx )
396       wall_flux  = 0.0
397       wall_index = NINT( a+ 2*b + 3*c1 + 4*c2 )
398
399!
400!--    All subsequent variables are computed for the respective location where
401!--    the respective flux is defined.
402       DO  k = nzb_w, nzt_w
403
404!
405!--       (1) Compute rifs, u_i, v_i, ws, pt' and w'pt'
406          rifs  = rif_wall(k,j,i,wall_index)
407
408          u_i   = a * u(k,j,i) + c1 * 0.25 * &
409                  ( u(k+1,j,i+1) + u(k+1,j,i) + u(k,j,i+1) + u(k,j,i) )
410
411          v_i   = b * v(k,j,i) + c2 * 0.25 * &
412                  ( v(k+1,j+1,i) + v(k+1,j,i) + v(k,j+1,i) + v(k,j,i) )
413
414          ws    = ( c1 + c2 ) * w(k,j,i) + 0.25 * (                            &
415                     a * ( w(k-1,j,i-1) + w(k-1,j,i) + w(k,j,i-1) + w(k,j,i) ) &
416                   + b * ( w(k-1,j-1,i) + w(k-1,j,i) + w(k,j-1,i) + w(k,j,i) ) &
417                                                  )
418          pt_i  = 0.5 * ( pt(k,j,i) + a *  pt(k,j,i-1) + b * pt(k,j-1,i)  &
419                          + ( c1 + c2 ) * pt(k+1,j,i) )
420
421          pts   = pt_i - hom(k,1,4,0)
422          wspts = ws * pts
423
424!
425!--       (2) Compute wall-parallel absolute velocity vel_total
426          vel_total = SQRT( ws**2 + ( a+c1 ) * u_i**2 + ( b+c2 ) * v_i**2 )
427
428!
429!--       (3) Compute wall friction velocity us_wall
430          IF ( rifs >= 0.0 )  THEN
431
432!
433!--          Stable stratification (and neutral)
434             us_wall = kappa * vel_total / ( LOG( zp / z0(j,i) ) +             &
435                                            5.0 * rifs * ( zp - z0(j,i) ) / zp &
436                                           )
437          ELSE
438
439!
440!--          Unstable stratification
441             h1 = SQRT( SQRT( 1.0 - 16.0 * rifs ) )
442             h2 = SQRT( SQRT( 1.0 - 16.0 * rifs * z0(j,i) / zp ) )
443
444             us_wall = kappa * vel_total / (                          &
445                  LOG( zp / z0(j,i) ) -                               &
446                  LOG( ( 1.0 + h1 )**2 * ( 1.0 + h1**2 ) / (          &
447                       ( 1.0 + h2 )**2 * ( 1.0 + h2**2 )   ) ) +      &
448                       2.0 * ( ATAN( h1 ) - ATAN( h2 ) )              &
449                                           )
450          ENDIF
451
452!
453!--       (4) Compute zp/L (corresponds to neutral Richardson flux number
454!--           rifs)
455          rifs = -1.0 * zp * kappa * g * wspts / ( pt_i * (us_wall**3 + 1E-30) )
456
457!
458!--       Limit the value range of the Richardson numbers.
459!--       This is necessary for very small velocities (u,w --> 0), because
460!--       the absolute value of rif can then become very large, which in
461!--       consequence would result in very large shear stresses and very
462!--       small momentum fluxes (both are generally unrealistic).
463          IF ( rifs < rif_min )  rifs = rif_min
464          IF ( rifs > rif_max )  rifs = rif_max
465
466!
467!--       (5) Compute wall_flux (u'v', v'u', w'v', or w'u')
468          IF ( rifs >= 0.0 )  THEN
469
470!
471!--          Stable stratification (and neutral)
472             wall_flux(k) = kappa *                                          &
473                            ( a*u(k,j,i) + b*v(k,j,i) + (c1+c2)*w(k,j,i) ) / &
474                            (  LOG( zp / z0(j,i) ) +                         &
475                               5.0 * rifs * ( zp - z0(j,i) ) / zp            &
476                            )
477          ELSE
478
479!
480!--          Unstable stratification
481             h1 = SQRT( SQRT( 1.0 - 16.0 * rifs ) )
482             h2 = SQRT( SQRT( 1.0 - 16.0 * rifs * z0(j,i) / zp ) )
483
484             wall_flux(k) = kappa *                               &
485                  ( a*u(k,j,i) + b*v(k,j,i) + (c1+c2)*w(k,j,i) ) / (  &
486                  LOG( zp / z0(j,i) ) -                               &
487                  LOG( ( 1.0 + h1 )**2 * ( 1.0 + h1**2 ) / (          &
488                       ( 1.0 + h2 )**2 * ( 1.0 + h2**2 )   ) ) +      &
489                       2.0 * ( ATAN( h1 ) - ATAN( h2 ) )              &
490                                                                   )
491          ENDIF
492          wall_flux(k) = -wall_flux(k) * us_wall
493
494!
495!--       store rifs for next time step
496          rif_wall(k,j,i,wall_index) = rifs
497
498       ENDDO
499
500    END SUBROUTINE wall_fluxes_ij
501
502
503
504!------------------------------------------------------------------------------!
505! Call for all grid points
506!------------------------------------------------------------------------------!
507    SUBROUTINE wall_fluxes_e( wall_flux, a, b, c1, c2, wall )
508
509!------------------------------------------------------------------------------!
510! Description:
511! ------------
512! Calculates momentum fluxes at vertical walls for routine production_e
513! assuming Monin-Obukhov similarity.
514! Indices: usvs a=1, vsus b=1, wsvs c1=1, wsus c2=1 (other=0).
515!------------------------------------------------------------------------------!
516
517       USE arrays_3d
518       USE control_parameters
519       USE grid_variables
520       USE indices
521       USE statistics
522
523       IMPLICIT NONE
524
525       INTEGER ::  i, j, k, kk, wall_index
526       REAL    ::  a, b, c1, c2, h1, h2, u_i, v_i, us_wall, vel_total, vel_zp, &
527                   ws, zp
528
529       REAL ::  rifs
530
531       REAL, DIMENSION(nysg:nyng,nxlg:nxrg)   ::  wall
532       REAL, DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) ::  wall_flux
533
534
535       zp         = 0.5 * ( (a+c1) * dy + (b+c2) * dx )
536       wall_flux  = 0.0
537       wall_index = NINT( a+ 2*b + 3*c1 + 4*c2 )
538
539       DO  i = nxl, nxr
540          DO  j = nys, nyn
541
542             IF ( wall(j,i) /= 0.0 )  THEN
543!
544!--             All subsequent variables are computed for scalar locations.
545                DO  k = nzb_diff_s_inner(j,i)-1, nzb_diff_s_outer(j,i)-2
546!
547!--                (1) Compute rifs, u_i, v_i, and ws
548                   IF ( k == nzb_diff_s_inner(j,i)-1 )  THEN
549                      kk = nzb_diff_s_inner(j,i)-1
550                   ELSE
551                      kk = k-1
552                   ENDIF
553                   rifs  = 0.5 * ( rif_wall(k,j,i,wall_index) +                &
554                          a * rif_wall(k,j,i+1,1) +  b * rif_wall(k,j+1,i,2) + &
555                          c1 * rif_wall(kk,j,i,3) + c2 * rif_wall(kk,j,i,4)    &
556                                 )
557
558                   u_i   = 0.5 * ( u(k,j,i) + u(k,j,i+1) )
559                   v_i   = 0.5 * ( v(k,j,i) + v(k,j+1,i) )
560                   ws    = 0.5 * ( w(k,j,i) + w(k-1,j,i) )
561!
562!--                (2) Compute wall-parallel absolute velocity vel_total and
563!--                interpolate appropriate velocity component vel_zp.
564                   vel_total = SQRT( ws**2 + (a+c1) * u_i**2 + (b+c2) * v_i**2 )
565                   vel_zp = 0.5 * ( a * u_i + b * v_i + (c1+c2) * ws )
566!
567!--                (3) Compute wall friction velocity us_wall
568                   IF ( rifs >= 0.0 )  THEN
569
570!
571!--                   Stable stratification (and neutral)
572                      us_wall = kappa * vel_total / ( LOG( zp / z0(j,i) ) +    &
573                                            5.0 * rifs * ( zp - z0(j,i) ) / zp &
574                                                    )
575                   ELSE
576
577!
578!--                   Unstable stratification
579                      h1 = SQRT( SQRT( 1.0 - 16.0 * rifs ) )
580                      h2 = SQRT( SQRT( 1.0 - 16.0 * rifs * z0(j,i) / zp ) )
581
582                      us_wall = kappa * vel_total / (                          &
583                           LOG( zp / z0(j,i) ) -                               &
584                           LOG( ( 1.0 + h1 )**2 * ( 1.0 + h1**2 ) / (          &
585                                ( 1.0 + h2 )**2 * ( 1.0 + h2**2 )   ) ) +      &
586                                2.0 * ( ATAN( h1 ) - ATAN( h2 ) )              &
587                                                    )
588                   ENDIF
589
590!
591!--                Skip step (4) of wall_fluxes, because here rifs is already
592!--                available from (1)
593!
594!--                (5) Compute wall_flux (u'v', v'u', w'v', or w'u')
595
596                   IF ( rifs >= 0.0 )  THEN
597
598!
599!--                   Stable stratification (and neutral)
600                      wall_flux(k,j,i) = kappa *  vel_zp / &
601                          ( LOG( zp/z0(j,i) ) + 5.0*rifs * ( zp-z0(j,i) ) / zp )
602                   ELSE
603
604!
605!--                   Unstable stratification
606                      h1 = SQRT( SQRT( 1.0 - 16.0 * rifs ) )
607                      h2 = SQRT( SQRT( 1.0 - 16.0 * rifs * z0(j,i) / zp ) )
608
609                      wall_flux(k,j,i) = kappa * vel_zp / (                    &
610                           LOG( zp / z0(j,i) ) -                               &
611                           LOG( ( 1.0 + h1 )**2 * ( 1.0 + h1**2 ) / (          &
612                                ( 1.0 + h2 )**2 * ( 1.0 + h2**2 )   ) ) +      &
613                                2.0 * ( ATAN( h1 ) - ATAN( h2 ) )              &
614                                                          )
615                   ENDIF
616                   wall_flux(k,j,i) = - wall_flux(k,j,i) * us_wall
617
618                ENDDO
619
620             ENDIF
621
622          ENDDO
623       ENDDO
624
625    END SUBROUTINE wall_fluxes_e
626
627
628!------------------------------------------------------------------------------!
629! Call for all grid points - accelerator version
630!------------------------------------------------------------------------------!
631    SUBROUTINE wall_fluxes_e_acc( wall_flux, a, b, c1, c2, wall )
632
633!------------------------------------------------------------------------------!
634! Description:
635! ------------
636! Calculates momentum fluxes at vertical walls for routine production_e
637! assuming Monin-Obukhov similarity.
638! Indices: usvs a=1, vsus b=1, wsvs c1=1, wsus c2=1 (other=0).
639!------------------------------------------------------------------------------!
640
641       USE arrays_3d
642       USE control_parameters
643       USE grid_variables
644       USE indices
645       USE statistics
646
647       IMPLICIT NONE
648
649       INTEGER ::  i, j, k, kk, max_outer, min_inner, wall_index
650       REAL    ::  a, b, c1, c2, h1, h2, u_i, v_i, us_wall, vel_total, vel_zp, &
651                   ws, zp
652
653       REAL ::  rifs
654
655       REAL, DIMENSION(nysg:nyng,nxlg:nxrg)   ::  wall
656       REAL, DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) ::  wall_flux
657
658
659       zp         = 0.5 * ( (a+c1) * dy + (b+c2) * dx )
660       wall_flux  = 0.0
661       wall_index = NINT( a+ 2*b + 3*c1 + 4*c2 )
662
663       min_inner = MINVAL( nzb_diff_s_inner(nys:nyn,nxl:nxr) ) - 1
664       max_outer = MAXVAL( nzb_diff_s_outer(nys:nyn,nxl:nxr) ) - 2
665
666       !$acc kernels present( nzb_diff_s_inner, nzb_diff_s_outer, pt, rif_wall ) &
667       !$acc         present( u, v, w, wall, wall_flux, z0 )
668       !$acc loop
669       DO  i = i_left, i_right
670          DO  j = j_south, j_north
671             !$acc loop vector(32)
672             DO  k = min_inner, max_outer
673!
674!--             All subsequent variables are computed for scalar locations
675                IF ( k >= nzb_diff_s_inner(j,i)-1  .AND. &
676                     k <= nzb_diff_s_outer(j,i)-2  .AND.  wall(j,i) /= 0.0 )  THEN
677!
678!--                (1) Compute rifs, u_i, v_i, and ws
679                   IF ( k == nzb_diff_s_inner(j,i)-1 )  THEN
680                      kk = nzb_diff_s_inner(j,i)-1
681                   ELSE
682                      kk = k-1
683                   ENDIF
684                   rifs  = 0.5 * ( rif_wall(k,j,i,wall_index) +                &
685                          a * rif_wall(k,j,i+1,1) +  b * rif_wall(k,j+1,i,2) + &
686                          c1 * rif_wall(kk,j,i,3) + c2 * rif_wall(kk,j,i,4)    &
687                                 )
688
689                   u_i   = 0.5 * ( u(k,j,i) + u(k,j,i+1) )
690                   v_i   = 0.5 * ( v(k,j,i) + v(k,j+1,i) )
691                   ws    = 0.5 * ( w(k,j,i) + w(k-1,j,i) )
692!
693!--                (2) Compute wall-parallel absolute velocity vel_total and
694!--                interpolate appropriate velocity component vel_zp.
695                   vel_total = SQRT( ws**2 + (a+c1) * u_i**2 + (b+c2) * v_i**2 )
696                   vel_zp = 0.5 * ( a * u_i + b * v_i + (c1+c2) * ws )
697!
698!--                (3) Compute wall friction velocity us_wall
699                   IF ( rifs >= 0.0 )  THEN
700
701!
702!--                   Stable stratification (and neutral)
703                      us_wall = kappa * vel_total / ( LOG( zp / z0(j,i) ) +    &
704                                            5.0 * rifs * ( zp - z0(j,i) ) / zp &
705                                                    )
706                   ELSE
707
708!
709!--                   Unstable stratification
710                      h1 = SQRT( SQRT( 1.0 - 16.0 * rifs ) )
711                      h2 = SQRT( SQRT( 1.0 - 16.0 * rifs * z0(j,i) / zp ) )
712
713                      us_wall = kappa * vel_total / (                          &
714                           LOG( zp / z0(j,i) ) -                               &
715                           LOG( ( 1.0 + h1 )**2 * ( 1.0 + h1**2 ) / (          &
716                                ( 1.0 + h2 )**2 * ( 1.0 + h2**2 )   ) ) +      &
717                                2.0 * ( ATAN( h1 ) - ATAN( h2 ) )              &
718                                                    )
719                   ENDIF
720
721!
722!--                Skip step (4) of wall_fluxes, because here rifs is already
723!--                available from (1)
724!
725!--                (5) Compute wall_flux (u'v', v'u', w'v', or w'u')
726
727                   IF ( rifs >= 0.0 )  THEN
728
729!
730!--                   Stable stratification (and neutral)
731                      wall_flux(k,j,i) = kappa *  vel_zp / &
732                          ( LOG( zp/z0(j,i) ) + 5.0*rifs * ( zp-z0(j,i) ) / zp )
733                   ELSE
734
735!
736!--                   Unstable stratification
737                      h1 = SQRT( SQRT( 1.0 - 16.0 * rifs ) )
738                      h2 = SQRT( SQRT( 1.0 - 16.0 * rifs * z0(j,i) / zp ) )
739
740                      wall_flux(k,j,i) = kappa * vel_zp / (                    &
741                           LOG( zp / z0(j,i) ) -                               &
742                           LOG( ( 1.0 + h1 )**2 * ( 1.0 + h1**2 ) / (          &
743                                ( 1.0 + h2 )**2 * ( 1.0 + h2**2 )   ) ) +      &
744                                2.0 * ( ATAN( h1 ) - ATAN( h2 ) )              &
745                                                          )
746                   ENDIF
747                   wall_flux(k,j,i) = - wall_flux(k,j,i) * us_wall
748
749                ENDIF
750
751             ENDDO
752          ENDDO
753       ENDDO
754       !$acc end kernels
755
756    END SUBROUTINE wall_fluxes_e_acc
757
758
759!------------------------------------------------------------------------------!
760! Call for grid point i,j
761!------------------------------------------------------------------------------!
762    SUBROUTINE wall_fluxes_e_ij( i, j, nzb_w, nzt_w, wall_flux, a, b, c1, c2 )
763
764       USE arrays_3d
765       USE control_parameters
766       USE grid_variables
767       USE indices
768       USE statistics
769
770       IMPLICIT NONE
771
772       INTEGER ::  i, j, k, kk, nzb_w, nzt_w, wall_index
773       REAL    ::  a, b, c1, c2, h1, h2, u_i, v_i, us_wall, vel_total, vel_zp, &
774                   ws, zp
775
776       REAL ::  rifs
777
778       REAL, DIMENSION(nzb:nzt+1) ::  wall_flux
779
780
781       zp         = 0.5 * ( (a+c1) * dy + (b+c2) * dx )
782       wall_flux  = 0.0
783       wall_index = NINT( a+ 2*b + 3*c1 + 4*c2 )
784
785!
786!--    All subsequent variables are computed for scalar locations.
787       DO  k = nzb_w, nzt_w
788
789!
790!--       (1) Compute rifs, u_i, v_i, and ws
791          IF ( k == nzb_w )  THEN
792             kk = nzb_w
793          ELSE
794             kk = k-1
795          ENDIF
796          rifs  = 0.5 * ( rif_wall(k,j,i,wall_index) +                         &
797                          a * rif_wall(k,j,i+1,1) +  b * rif_wall(k,j+1,i,2) + &
798                          c1 * rif_wall(kk,j,i,3) + c2 * rif_wall(kk,j,i,4)    &
799                        )
800
801          u_i   = 0.5 * ( u(k,j,i) + u(k,j,i+1) )
802          v_i   = 0.5 * ( v(k,j,i) + v(k,j+1,i) )
803          ws    = 0.5 * ( w(k,j,i) + w(k-1,j,i) )
804!
805!--       (2) Compute wall-parallel absolute velocity vel_total and
806!--       interpolate appropriate velocity component vel_zp.
807          vel_total = SQRT( ws**2 + (a+c1) * u_i**2 + (b+c2) * v_i**2 )
808          vel_zp = 0.5 * ( a * u_i + b * v_i + (c1+c2) * ws )
809!
810!--       (3) Compute wall friction velocity us_wall
811          IF ( rifs >= 0.0 )  THEN
812
813!
814!--          Stable stratification (and neutral)
815             us_wall = kappa * vel_total / ( LOG( zp / z0(j,i) ) +             &
816                                            5.0 * rifs * ( zp - z0(j,i) ) / zp &
817                                           )
818          ELSE
819
820!
821!--          Unstable stratification
822             h1 = SQRT( SQRT( 1.0 - 16.0 * rifs ) )
823             h2 = SQRT( SQRT( 1.0 - 16.0 * rifs * z0(j,i) / zp ) )
824
825             us_wall = kappa * vel_total / (                          &
826                  LOG( zp / z0(j,i) ) -                               &
827                  LOG( ( 1.0 + h1 )**2 * ( 1.0 + h1**2 ) / (          &
828                       ( 1.0 + h2 )**2 * ( 1.0 + h2**2 )   ) ) +      &
829                       2.0 * ( ATAN( h1 ) - ATAN( h2 ) )              &
830                                           )
831          ENDIF
832
833!
834!--       Skip step (4) of wall_fluxes, because here rifs is already
835!--       available from (1)
836!
837!--       (5) Compute wall_flux (u'v', v'u', w'v', or w'u')
838!--       First interpolate the velocity (this is different from
839!--       subroutine wall_fluxes because fluxes in subroutine
840!--       wall_fluxes_e are defined at scalar locations).
841          vel_zp = 0.5 * (       a * ( u(k,j,i) + u(k,j,i+1) ) +  &
842                                 b * ( v(k,j,i) + v(k,j+1,i) ) +  &
843                           (c1+c2) * ( w(k,j,i) + w(k-1,j,i) )    &
844                         )
845
846          IF ( rifs >= 0.0 )  THEN
847
848!
849!--          Stable stratification (and neutral)
850             wall_flux(k) = kappa *  vel_zp / &
851                          ( LOG( zp/z0(j,i) ) + 5.0*rifs * ( zp-z0(j,i) ) / zp )
852          ELSE
853
854!
855!--          Unstable stratification
856             h1 = SQRT( SQRT( 1.0 - 16.0 * rifs ) )
857             h2 = SQRT( SQRT( 1.0 - 16.0 * rifs * z0(j,i) / zp ) )
858
859             wall_flux(k) = kappa * vel_zp / (                        &
860                  LOG( zp / z0(j,i) ) -                               &
861                  LOG( ( 1.0 + h1 )**2 * ( 1.0 + h1**2 ) / (          &
862                       ( 1.0 + h2 )**2 * ( 1.0 + h2**2 )   ) ) +      &
863                       2.0 * ( ATAN( h1 ) - ATAN( h2 ) )              &
864                                                 )
865          ENDIF
866          wall_flux(k) = - wall_flux(k) * us_wall
867
868       ENDDO
869
870    END SUBROUTINE wall_fluxes_e_ij
871
872 END MODULE wall_fluxes_mod
Note: See TracBrowser for help on using the repository browser.