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

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

code adjustment of accelerator version for PGI 12.3 / CUDA 5.0

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