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

Last change on this file since 1132 was 1132, checked in by raasch, 11 years ago

r1028 documented

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