Ignore:
Timestamp:
Sep 27, 2012 9:23:24 AM (9 years ago)
Author:
raasch
Message:

Starting with changes required for GPU optimization. OpenACC statements for using NVIDIA GPUs added.
Adjustment of mixing length to the Prandtl mixing length at first grid point above ground removed.
mask array is set zero for ghost boundaries

File:
1 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk/SOURCE/wall_fluxes.f90

    r667 r1015  
    33! Current revisions:
    44! -----------------
    5 !
     5! accelerator version (*_acc) added
    66!
    77! Former revisions:
     
    3232!------------------------------------------------------------------------------!
    3333    PRIVATE
    34     PUBLIC wall_fluxes, wall_fluxes_e
     34    PUBLIC wall_fluxes, wall_fluxes_acc, wall_fluxes_e, wall_fluxes_e_acc
    3535   
    3636    INTERFACE wall_fluxes
     
    3939    END INTERFACE wall_fluxes
    4040   
     41    INTERFACE wall_fluxes_acc
     42       MODULE PROCEDURE wall_fluxes_acc
     43    END INTERFACE wall_fluxes_acc
     44
    4145    INTERFACE wall_fluxes_e
    4246       MODULE PROCEDURE wall_fluxes_e
     
    4448    END INTERFACE wall_fluxes_e
    4549 
     50    INTERFACE wall_fluxes_e_acc
     51       MODULE PROCEDURE wall_fluxes_e_acc
     52    END INTERFACE wall_fluxes_e_acc
     53
    4654 CONTAINS
    4755
     
    189197    END SUBROUTINE wall_fluxes
    190198
     199
     200!------------------------------------------------------------------------------!
     201! Call for all grid points - accelerator version
     202!------------------------------------------------------------------------------!
     203    SUBROUTINE wall_fluxes_acc( wall_flux, a, b, c1, c2, nzb_uvw_inner, &
     204                                nzb_uvw_outer, wall )
     205
     206       USE arrays_3d
     207       USE control_parameters
     208       USE grid_variables
     209       USE indices
     210       USE statistics
     211
     212       IMPLICIT NONE
     213
     214       INTEGER ::  i, j, k, max_outer, min_inner, wall_index
     215
     216       INTEGER, DIMENSION(nysg:nyng,nxlg:nxrg) ::  nzb_uvw_inner, &
     217                                                   nzb_uvw_outer
     218       REAL ::  a, b, c1, c2, h1, h2, zp
     219       REAL ::  pts, pt_i, rifs, u_i, v_i, us_wall, vel_total, ws, wspts
     220
     221       REAL, DIMENSION(nysg:nyng,nxlg:nxrg)   ::  wall
     222       REAL, DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) ::  wall_flux
     223
     224
     225       zp         = 0.5 * ( (a+c1) * dy + (b+c2) * dx )
     226       wall_flux  = 0.0
     227       wall_index = NINT( a+ 2*b + 3*c1 + 4*c2 )
     228
     229       min_inner = MINVAL( nzb_uvw_inner(nys:nyn,nxl:nxr) ) + 1
     230       max_outer = MINVAL( nzb_uvw_outer(nys:nyn,nxl:nxr) )
     231
     232       !$acc kernels present( hom, nzb_uvw_inner, nzb_uvw_outer, pt, rif_wall ) &
     233       !$acc         present( u, v, w, wall, wall_flux, z0 )
     234       !$acc loop
     235       DO  i = nxl, nxr
     236          DO  j = nys, nyn
     237             !$acc loop vector( 32 )
     238             DO  k = min_inner, max_outer
     239!
     240!--             All subsequent variables are computed for the respective
     241!--             location where the respective flux is defined.
     242                IF ( k >= nzb_uvw_inner(j,i)+1  .AND. &
     243                     k <= nzb_uvw_outer(j,i)    .AND.  wall(j,i) /= 0.0 )  THEN
     244!
     245!--                (1) Compute rifs, u_i, v_i, ws, pt' and w'pt'
     246                   rifs  = rif_wall(k,j,i,wall_index)
     247
     248                   u_i   = a * u(k,j,i) + c1 * 0.25 * &
     249                           ( u(k+1,j,i+1) + u(k+1,j,i) + u(k,j,i+1) + u(k,j,i) )
     250
     251                   v_i   = b * v(k,j,i) + c2 * 0.25 * &
     252                           ( v(k+1,j+1,i) + v(k+1,j,i) + v(k,j+1,i) + v(k,j,i) )
     253
     254                   ws    = ( c1 + c2 ) * w(k,j,i) + 0.25 * (                   &
     255                     a * ( w(k-1,j,i-1) + w(k-1,j,i) + w(k,j,i-1) + w(k,j,i) ) &
     256                   + b * ( w(k-1,j-1,i) + w(k-1,j,i) + w(k,j-1,i) + w(k,j,i) ) &
     257                                                           )
     258                   pt_i  = 0.5 * ( pt(k,j,i) + a *  pt(k,j,i-1) + &
     259                                   b * pt(k,j-1,i) + ( c1 + c2 ) * pt(k+1,j,i) )
     260
     261                   pts   = pt_i - hom(k,1,4,0)
     262                   wspts = ws * pts
     263
     264!
     265!--                (2) Compute wall-parallel absolute velocity vel_total
     266                   vel_total = SQRT( ws**2 + (a+c1) * u_i**2 + (b+c2) * v_i**2 )
     267
     268!
     269!--                (3) Compute wall friction velocity us_wall
     270                   IF ( rifs >= 0.0 )  THEN
     271
     272!
     273!--                   Stable stratification (and neutral)
     274                      us_wall = kappa * vel_total / ( LOG( zp / z0(j,i) ) +    &
     275                                            5.0 * rifs * ( zp - z0(j,i) ) / zp &
     276                                                    )
     277                   ELSE
     278
     279!
     280!--                   Unstable stratification
     281                      h1 = SQRT( SQRT( 1.0 - 16.0 * rifs ) )
     282                      h2 = SQRT( SQRT( 1.0 - 16.0 * rifs * z0(j,i) / zp ) )
     283
     284                      us_wall = kappa * vel_total / (                          &
     285                           LOG( zp / z0(j,i) ) -                               &
     286                           LOG( ( 1.0 + h1 )**2 * ( 1.0 + h1**2 ) / (          &
     287                                ( 1.0 + h2 )**2 * ( 1.0 + h2**2 )   ) ) +      &
     288                                2.0 * ( ATAN( h1 ) - ATAN( h2 ) )              &
     289                                                    )
     290                   ENDIF
     291
     292!
     293!--                (4) Compute zp/L (corresponds to neutral Richardson flux
     294!--                    number rifs)
     295                   rifs = -1.0 * zp * kappa * g * wspts / ( pt_i * &
     296                                                        ( us_wall**3 + 1E-30 ) )
     297
     298!
     299!--                Limit the value range of the Richardson numbers.
     300!--                This is necessary for very small velocities (u,w --> 0),
     301!--                because the absolute value of rif can then become very
     302!--                large, which in consequence would result in very large
     303!--                shear stresses and very small momentum fluxes (both are
     304!--                generally unrealistic).
     305                   IF ( rifs < rif_min )  rifs = rif_min
     306                   IF ( rifs > rif_max )  rifs = rif_max
     307
     308!
     309!--                (5) Compute wall_flux (u'v', v'u', w'v', or w'u')
     310                   IF ( rifs >= 0.0 )  THEN
     311
     312!
     313!--                   Stable stratification (and neutral)
     314                      wall_flux(k,j,i) = kappa *                               &
     315                              ( a*u(k,j,i) + b*v(k,j,i) + (c1+c2)*w(k,j,i) ) / &
     316                              (  LOG( zp / z0(j,i) ) +                         &
     317                                 5.0 * rifs * ( zp - z0(j,i) ) / zp            &
     318                              )
     319                   ELSE
     320
     321!
     322!--                   Unstable stratification
     323                      h1 = SQRT( SQRT( 1.0 - 16.0 * rifs ) )
     324                      h2 = SQRT( SQRT( 1.0 - 16.0 * rifs * z0(j,i) / zp ) )
     325
     326                      wall_flux(k,j,i) = kappa *                               &
     327                           ( a*u(k,j,i) + b*v(k,j,i) + (c1+c2)*w(k,j,i) ) / (  &
     328                           LOG( zp / z0(j,i) ) -                               &
     329                           LOG( ( 1.0 + h1 )**2 * ( 1.0 + h1**2 ) / (          &
     330                                ( 1.0 + h2 )**2 * ( 1.0 + h2**2 )   ) ) +      &
     331                                2.0 * ( ATAN( h1 ) - ATAN( h2 ) )              &
     332                                                                            )
     333                   ENDIF
     334                   wall_flux(k,j,i) = -wall_flux(k,j,i) * us_wall
     335
     336!
     337!--                store rifs for next time step
     338                   rif_wall(k,j,i,wall_index) = rifs
     339
     340                ENDIF
     341
     342             ENDDO
     343          ENDDO
     344       ENDDO
     345       !$acc end kernels
     346
     347    END SUBROUTINE wall_fluxes_acc
    191348
    192349
     
    444601    END SUBROUTINE wall_fluxes_e
    445602
     603
     604!------------------------------------------------------------------------------!
     605! Call for all grid points - accelerator version
     606!------------------------------------------------------------------------------!
     607    SUBROUTINE wall_fluxes_e_acc( wall_flux, a, b, c1, c2, wall )
     608
     609!------------------------------------------------------------------------------!
     610! Description:
     611! ------------
     612! Calculates momentum fluxes at vertical walls for routine production_e
     613! assuming Monin-Obukhov similarity.
     614! Indices: usvs a=1, vsus b=1, wsvs c1=1, wsus c2=1 (other=0).
     615!------------------------------------------------------------------------------!
     616
     617       USE arrays_3d
     618       USE control_parameters
     619       USE grid_variables
     620       USE indices
     621       USE statistics
     622
     623       IMPLICIT NONE
     624
     625       INTEGER ::  i, j, k, kk, max_outer, min_inner, wall_index
     626       REAL    ::  a, b, c1, c2, h1, h2, u_i, v_i, us_wall, vel_total, vel_zp, &
     627                   ws, zp
     628
     629       REAL ::  rifs
     630
     631       REAL, DIMENSION(nysg:nyng,nxlg:nxrg)   ::  wall
     632       REAL, DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) ::  wall_flux
     633
     634
     635       zp         = 0.5 * ( (a+c1) * dy + (b+c2) * dx )
     636       wall_flux  = 0.0
     637       wall_index = NINT( a+ 2*b + 3*c1 + 4*c2 )
     638
     639       min_inner = MINVAL( nzb_diff_s_inner(nys:nyn,nxl:nxr) ) - 1
     640       max_outer = MAXVAL( nzb_diff_s_outer(nys:nyn,nxl:nxr) ) - 2
     641
     642       !$acc kernels present( nzb_diff_s_inner, nzb_diff_s_outer, pt, rif_wall ) &
     643       !$acc         present( u, v, w, wall, wall_flux, z0 )
     644       !$acc loop
     645       DO  i = nxl, nxr
     646          DO  j = nys, nyn
     647             !$acc loop vector(32)
     648             DO  k = min_inner, max_outer
     649!
     650!--             All subsequent variables are computed for scalar locations
     651                IF ( k >= nzb_diff_s_inner(j,i)-1  .AND. &
     652                     k <= nzb_diff_s_outer(j,i)-2  .AND.  wall(j,i) /= 0.0 )  THEN
     653!
     654!--                (1) Compute rifs, u_i, v_i, and ws
     655                   IF ( k == nzb_diff_s_inner(j,i)-1 )  THEN
     656                      kk = nzb_diff_s_inner(j,i)-1
     657                   ELSE
     658                      kk = k-1
     659                   ENDIF
     660                   rifs  = 0.5 * ( rif_wall(k,j,i,wall_index) +                &
     661                          a * rif_wall(k,j,i+1,1) +  b * rif_wall(k,j+1,i,2) + &
     662                          c1 * rif_wall(kk,j,i,3) + c2 * rif_wall(kk,j,i,4)    &
     663                                 )
     664
     665                   u_i   = 0.5 * ( u(k,j,i) + u(k,j,i+1) )
     666                   v_i   = 0.5 * ( v(k,j,i) + v(k,j+1,i) )
     667                   ws    = 0.5 * ( w(k,j,i) + w(k-1,j,i) )
     668!
     669!--                (2) Compute wall-parallel absolute velocity vel_total and
     670!--                interpolate appropriate velocity component vel_zp.
     671                   vel_total = SQRT( ws**2 + (a+c1) * u_i**2 + (b+c2) * v_i**2 )
     672                   vel_zp = 0.5 * ( a * u_i + b * v_i + (c1+c2) * ws )
     673!
     674!--                (3) Compute wall friction velocity us_wall
     675                   IF ( rifs >= 0.0 )  THEN
     676
     677!
     678!--                   Stable stratification (and neutral)
     679                      us_wall = kappa * vel_total / ( LOG( zp / z0(j,i) ) +    &
     680                                            5.0 * rifs * ( zp - z0(j,i) ) / zp &
     681                                                    )
     682                   ELSE
     683
     684!
     685!--                   Unstable stratification
     686                      h1 = SQRT( SQRT( 1.0 - 16.0 * rifs ) )
     687                      h2 = SQRT( SQRT( 1.0 - 16.0 * rifs * z0(j,i) / zp ) )
     688
     689                      us_wall = kappa * vel_total / (                          &
     690                           LOG( zp / z0(j,i) ) -                               &
     691                           LOG( ( 1.0 + h1 )**2 * ( 1.0 + h1**2 ) / (          &
     692                                ( 1.0 + h2 )**2 * ( 1.0 + h2**2 )   ) ) +      &
     693                                2.0 * ( ATAN( h1 ) - ATAN( h2 ) )              &
     694                                                    )
     695                   ENDIF
     696
     697!
     698!--                Skip step (4) of wall_fluxes, because here rifs is already
     699!--                available from (1)
     700!
     701!--                (5) Compute wall_flux (u'v', v'u', w'v', or w'u')
     702
     703                   IF ( rifs >= 0.0 )  THEN
     704
     705!
     706!--                   Stable stratification (and neutral)
     707                      wall_flux(k,j,i) = kappa *  vel_zp / &
     708                          ( LOG( zp/z0(j,i) ) + 5.0*rifs * ( zp-z0(j,i) ) / zp )
     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                      wall_flux(k,j,i) = kappa * vel_zp / (                    &
     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                   wall_flux(k,j,i) = - wall_flux(k,j,i) * us_wall
     724
     725                ENDIF
     726
     727             ENDDO
     728          ENDDO
     729       ENDDO
     730       !$acc end kernels
     731
     732    END SUBROUTINE wall_fluxes_e_acc
    446733
    447734
Note: See TracChangeset for help on using the changeset viewer.