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

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

code has been put under the GNU General Public License (v3)

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