Changeset 4468 for palm


Ignore:
Timestamp:
Mar 23, 2020 1:49:05 PM (15 months ago)
Author:
suehring
Message:

bugfix for last commit in openacc branch; some loop bounds revised (only to be consistent with cache version, no effect on solution); setting of nzb_max_l for advection of the w-component revised (no effect on the solution)

File:
1 edited

Legend:

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

    r4466 r4468  
    2020! Current revisions:
    2121! ------------------
    22 !
    23 !
     22! 
     23! 
    2424! Former revisions:
    2525! -----------------
    2626! $Id$
     27! - bugfix for last commit in openacc branch
     28! - some loop bounds revised (only to be consistent with cache version)
     29! - setting of nzb_max_l for advection of the w-component revised
     30!
     31! 4466 2020-03-20 16:14:41Z suehring
    2732! - vector branch further optimized (linear dependencies along z removed and
    2833!   loops are splitted)
     
    18221827                                         )                                     &
    18231828                          ) * ddy                                              &
    1824                         + ( w(k,j,i) * rho_air_zw(k) *                         &
     1829                        + ( w(k,j,i)   * rho_air_zw(k) *                       &
    18251830                                         ( ibit6 + ibit7 + ibit8 )             &
    18261831                          - w(k-1,j,i) * rho_air_zw(k-1) *                     &
     
    18601865          div         =   ( u(k,j,i+1) - u(k,j,i) ) * ddx                      &
    18611866                        + ( v(k,j+1,i) - v(k,j,i) ) * ddy                      &
    1862                         + ( w(k,j,i) * rho_air_zw(k)                           &
     1867                        + ( w(k,j,i)   * rho_air_zw(k)                         &
    18631868                          - w(k-1,j,i) * rho_air_zw(k-1)                       &
    18641869                                                  ) * drho_air(k) * ddzw(k)
     
    31883193           ( bc_dirichlet_s  .OR.  bc_radiation_s )  .AND.  j <= nys + 2  .OR. &
    31893194           ( bc_dirichlet_n  .OR.  bc_radiation_n )  .AND.  j >= nyn - 2 )  THEN
    3190           nzb_max_l = nzt
     3195          nzb_max_l = nzt - 1
    31913196       ELSE
    31923197          nzb_max_l = nzb_max
     
    37053710       INTEGER(iwp) ::  nzb_max_l !< index indicating upper bound for order degradation of horizontal advection terms
    37063711       INTEGER(iwp) ::  sk_num    !< integer identifier, used for assign fluxes to the correct dimension in the analysis array
    3707        INTEGER(iwp) ::  tn = 0    !< number of OpenMP thread
     3712       INTEGER(iwp) ::  tn = 0    !< number of OpenMP thread (is always zero here)
    37083713
    37093714       INTEGER(iwp), INTENT(IN), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::   &
     
    37553760       REAL(wp), DIMENSION(nzb:nzt+1)         ::  flux_r     !< discretized 6th-order flux at rightward-side of the grid box
    37563761       REAL(wp), DIMENSION(nzb:nzt+1)         ::  flux_t     !< discretized 6th-order flux at top of the grid box
    3757 
    3758        REAL(wp), DIMENSION(nzb+1:nzt,0:threads_per_task-1) ::  swap_diss_y_local !< discretized artificial dissipation at southward-side
    3759        REAL(wp), DIMENSION(nzb+1:nzt,0:threads_per_task-1) ::  swap_flux_y_local !< discretized 6th-order flux at northward-side
    3760 
    3761        REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn,0:threads_per_task-1) ::  swap_diss_x_local !< discretized artificial dissipation at leftward-side
    3762        REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn,0:threads_per_task-1) ::  swap_flux_x_local !< discretized 6th-order flux at leftward-side
     3762#ifndef _OPENACC
     3763       REAL(wp), DIMENSION(nzb+1:nzt)         ::  swap_diss_y_local !< discretized artificial dissipation at southward-side
     3764       REAL(wp), DIMENSION(nzb+1:nzt)         ::  swap_flux_y_local !< discretized 6th-order flux at northward-side
     3765#endif
     3766#ifdef _OPENACC
     3767       REAL(wp), DIMENSION(nzb+1:nzt)         ::  diss_l            !< discretized artificial dissipation at leftward-side
     3768       REAL(wp), DIMENSION(nzb+1:nzt)         ::  diss_s            !< discretized artificial dissipation at southward-side
     3769       REAL(wp), DIMENSION(nzb+1:nzt)         ::  flux_l            !< discretized 6th-order flux at leftward-side
     3770       REAL(wp), DIMENSION(nzb+1:nzt)         ::  flux_s            !< discretized 6th-order flux at southward-side
     3771#endif
     3772
     3773       REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) ::  swap_diss_x_local !< discretized artificial dissipation at leftward-side
     3774       REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) ::  swap_flux_x_local !< discretized 6th-order flux at leftward-side
    37633775
    37643776       CALL cpu_log( log_point_s(49), 'advec_s_ws', 'start' )
     
    37953807       !$ACC PRIVATE(ibit3_s, ibit4_s, ibit5_s) &
    37963808       !$ACC PRIVATE(ibit6, ibit7, ibit8) &
    3797        !$ACC PRIVATE(flux_r, diss_r) &
    3798        !$ACC PRIVATE(flux_n, diss_n) &
    3799        !$ACC PRIVATE(swap_diss_y_local, swap_flux_y_local) &
     3809       !$ACC PRIVATE(nzb_max_l) &
     3810       !$ACC PRIVATE(diss_l, diss_r, flux_l, flux_r) &
     3811       !$ACC PRIVATE(diss_n, diss_s, flux_s, flux_n) &
    38003812       !$ACC PRIVATE(flux_t, diss_t, flux_d, diss_d) &
    38013813       !$ACC PRIVATE(div, u_comp, u_comp_l, v_comp, v_comp_s) &
     
    38383850
    38393851                   u_comp                    = u(k,j,i) - u_gtrans + u_stokes_zu(k)
    3840                    swap_flux_x_local(k,j,tn) = u_comp * (                      &
     3852                   swap_flux_x_local(k,j)    = u_comp * (                      &
    38413853                                               ( 37.0_wp * ibit2 * adv_sca_5   &
    38423854                                            +     7.0_wp * ibit1 * adv_sca_3   &
     
    38533865                                                        )
    38543866
    3855                    swap_diss_x_local(k,j,tn) = -ABS( u_comp ) * (              &
     3867                   swap_diss_x_local(k,j)    = -ABS( u_comp ) * (              &
    38563868                                                ( 10.0_wp * ibit2 * adv_sca_5  &
    38573869                                             +     3.0_wp * ibit1 * adv_sca_3  &
     
    38733885
    38743886                   u_comp                    = u(k,j,i) - u_gtrans + u_stokes_zu(k)
    3875                    swap_flux_x_local(k,j,tn) = u_comp * (                      &
     3887                   swap_flux_x_local(k,j)    = u_comp * (                      &
    38763888                                        37.0_wp * ( sk(k,j,i)   + sk(k,j,i-1) )&
    38773889                                      -  8.0_wp * ( sk(k,j,i+1) + sk(k,j,i-2) )&
     
    38793891                                                        ) * adv_sca_5
    38803892
    3881                    swap_diss_x_local(k,j,tn) = -ABS( u_comp ) * (              &
     3893                   swap_diss_x_local(k,j)    = -ABS( u_comp ) * (              &
    38823894                                        10.0_wp * ( sk(k,j,i)   - sk(k,j,i-1) )&
    38833895                                      -  5.0_wp * ( sk(k,j,i+1) - sk(k,j,i-2) )&
     
    39003912
    39013913                   v_comp                  = v(k,j,i) - v_gtrans + v_stokes_zu(k)
    3902                    swap_flux_y_local(k,tn) = v_comp *         (                &
     3914                   swap_flux_y_local(k)    = v_comp *         (                &
    39033915                                                ( 37.0_wp * ibit5 * adv_sca_5  &
    39043916                                             +     7.0_wp * ibit4 * adv_sca_3  &
     
    39153927                                                              )
    39163928
    3917                    swap_diss_y_local(k,tn) = -ABS( v_comp ) * (                &
     3929                   swap_diss_y_local(k)    = -ABS( v_comp ) * (                &
    39183930                                                ( 10.0_wp * ibit5 * adv_sca_5  &
    39193931                                             +     3.0_wp * ibit4 * adv_sca_3  &
     
    39363948
    39373949                   v_comp                  = v(k,j,i) - v_gtrans + v_stokes_zu(k)
    3938                    swap_flux_y_local(k,tn) = v_comp * (                        &
     3950                   swap_flux_y_local(k)    = v_comp * (                        &
    39393951                                     37.0_wp * ( sk(k,j,i)   + sk(k,j-1,i) )   &
    39403952                                   -  8.0_wp * ( sk(k,j+1,i) + sk(k,j-2,i) )   &
    39413953                                   +           ( sk(k,j+2,i) + sk(k,j-3,i) )   &
    39423954                                                      ) * adv_sca_5
    3943                    swap_diss_y_local(k,tn) = -ABS( v_comp ) * (                &
     3955                   swap_diss_y_local(k)    = -ABS( v_comp ) * (                &
    39443956                                     10.0_wp * ( sk(k,j,i)   - sk(k,j-1,i) )   &
    39453957                                   -  5.0_wp * ( sk(k,j+1,i) - sk(k,j-2,i) )   &
     
    40004012                ibit0_l = REAL( IBITS(advc_flag(k,j,i-1),0,1), KIND = wp )
    40014013
    4002                 u_comp_l = u(k,j,i) - u_gtrans + u_stokes_zu(k)
    4003                 swap_flux_x_local(k,j,tn) = u_comp_l * (                       &
     4014                u_comp_l                  = u(k,j,i) - u_gtrans + u_stokes_zu(k)
     4015                flux_l(k)                = u_comp_l * (                       &
    40044016                                             ( 37.0_wp * ibit2_l * adv_sca_5   &
    40054017                                          +     7.0_wp * ibit1_l * adv_sca_3   &
     
    40164028                                                       )
    40174029
    4018                 swap_diss_x_local(k,j,tn) = -ABS( u_comp_l ) * (               &
     4030                diss_l(k)                = -ABS( u_comp_l ) * (               &
    40194031                                            ( 10.0_wp * ibit2_l * adv_sca_5    &
    40204032                                         +     3.0_wp * ibit1_l * adv_sca_3    &
     
    40254037                                         +              ibit1_l * adv_sca_3    &
    40264038                                            ) *                                &
    4027                                         ( sk(k,j,i+1) - sk(k,j,i-2) )         &
     4039                                         ( sk(k,j,i+1) - sk(k,j,i-2) )         &
    40284040                                     +      (           ibit2_l * adv_sca_5    &
    40294041                                            ) *                                &
    40304042                                         ( sk(k,j,i+2) - sk(k,j,i-3) )         &
    40314043                                                               )
     4044#else
     4045                flux_l(k) = swap_flux_x_local(k,j)
     4046                diss_l(k) = swap_diss_x_local(k,j)
    40324047#endif
    40334048                ibit5 = REAL( IBITS(advc_flag(k,j,i),5,1), KIND = wp )
     
    40734088
    40744089                v_comp_s                = v(k,j,i) - v_gtrans + v_stokes_zu(k)
    4075                 swap_flux_y_local(k,tn) = v_comp_s * (                         &
     4090                flux_s(k)              = v_comp_s * (                         &
    40764091                                             ( 37.0_wp * ibit5_s * adv_sca_5   &
    40774092                                          +     7.0_wp * ibit4_s * adv_sca_3   &
    40784093                                          +              ibit3_s * adv_sca_1   &
    40794094                                             ) *                               &
    4080                                          ( sk(k,j,i)  + sk(k,j-1,i)     )      &
     4095                                         ( sk(k,j,i)   + sk(k,j-1,i)     )     &
    40814096                                       -     (  8.0_wp * ibit5_s * adv_sca_5   &
    40824097                                          +              ibit4_s * adv_sca_3   &
     
    40854100                                      +      (           ibit5_s * adv_sca_5   &
    40864101                                             ) *                               &
    4087                                         ( sk(k,j+2,i) + sk(k,j-3,i)     )      &
     4102                                         ( sk(k,j+2,i) + sk(k,j-3,i)    )      &
    40884103                                                     )
    40894104
    4090                 swap_flux_y_local(k,tn) = -ABS( v_comp_s ) * (                 &
     4105                diss_s(k)              = -ABS( v_comp_s ) * (                 &
    40914106                                             ( 10.0_wp * ibit5_s * adv_sca_5   &
    40924107                                          +     3.0_wp * ibit4_s * adv_sca_3   &
     
    41024117                                          ( sk(k,j+2,i) - sk(k,j-3,i)    )     &
    41034118                                                             )
     4119#else
     4120                flux_s(k) = swap_flux_y_local(k)
     4121                diss_s(k) = swap_diss_y_local(k)
    41044122#endif
    41054123             ENDDO
     
    41234141!--             Recompute the left fluxes.
    41244142                u_comp_l                  = u(k,j,i) - u_gtrans + u_stokes_zu(k)
    4125                 swap_flux_x_local(k,j,tn) = u_comp_l * (                       &
     4143                flux_l(k)                = u_comp_l * (                       &
    41264144                                      37.0_wp * ( sk(k,j,i)   + sk(k,j,i-1) )  &
    41274145                                    -  8.0_wp * ( sk(k,j,i+1) + sk(k,j,i-2) )  &
     
    41294147                                                       ) * adv_sca_5
    41304148
    4131                 swap_diss_x_local(k,j,tn) = -ABS( u_comp_l ) * (               &
     4149                diss_l(k)                = -ABS( u_comp_l ) * (               &
    41324150                                      10.0_wp * ( sk(k,j,i)   - sk(k,j,i-1) )  &
    41334151                                    -  5.0_wp * ( sk(k,j,i+1) - sk(k,j,i-2) )  &
    41344152                                    +           ( sk(k,j,i+2) - sk(k,j,i-3) )  &
    41354153                                                               ) * adv_sca_5
     4154#else
     4155                flux_l(k) = swap_flux_x_local(k,j)
     4156                diss_l(k) = swap_diss_x_local(k,j)
     4157
    41364158#endif
    41374159
     
    41494171!--             Recompute the south fluxes.
    41504172                v_comp_s                = v(k,j,i) - v_gtrans + v_stokes_zu(k)
    4151                 swap_flux_y_local(k,tn) = v_comp_s * (                         &
     4173                flux_s(k)              = v_comp_s * (                         &
    41524174                                    37.0_wp * ( sk(k,j,i)   + sk(k,j-1,i) )    &
    41534175                                  -  8.0_wp * ( sk(k,j+1,i) + sk(k,j-2,i) )    &
    41544176                                  +           ( sk(k,j+2,i) + sk(k,j-3,i) )    &
    41554177                                                     ) * adv_sca_5
    4156                 swap_diss_y_local(k,tn) = -ABS( v_comp_s ) * (                 &
     4178                diss_s(k)              = -ABS( v_comp_s ) * (                 &
    41574179                                    10.0_wp * ( sk(k,j,i)   - sk(k,j-1,i) )    &
    41584180                                  -  5.0_wp * ( sk(k,j+1,i) - sk(k,j-2,i) )    &
    41594181                                  +             sk(k,j+2,i) - sk(k,j-3,i)      &
    41604182                                                             ) * adv_sca_5
     4183#else
     4184                flux_s(k) = swap_flux_y_local(k)
     4185                diss_s(k) = swap_diss_y_local(k)
    41614186#endif
    41624187             ENDDO
     
    43044329             DO  k = nzb+1, nzb_max_l
    43054330
    4306                 flux_d    = flux_t(k-1)
    4307                 diss_d    = diss_t(k-1)
     4331                flux_d = flux_t(k-1)
     4332                diss_d = diss_t(k-1)
    43084333
    43094334                ibit2 = REAL( IBITS(advc_flag(k,j,i),2,1), KIND = wp )
     
    43364361                                               )                               &
    43374362                                ) * ddy                                        &
    4338                               + ( w(k,j,i) * rho_air_zw(k) *                   &
     4363                              + ( w(k,j,i)   * rho_air_zw(k)   *               &
    43394364                                               ( ibit6 + ibit7 + ibit8 )       &
    43404365                                - w(k-1,j,i) * rho_air_zw(k-1) *               &
     
    43474372
    43484373                tend(k,j,i) = tend(k,j,i) - (                                  &
    4349                          ( flux_r(k) + diss_r(k) - swap_flux_x_local(k,j,tn) - &
    4350                            swap_diss_x_local(k,j,tn)            ) * ddx        &
    4351                        + ( flux_n(k) + diss_n(k) - swap_flux_y_local(k,tn)   - &
    4352                            swap_diss_y_local(k,tn)              ) * ddy        &
    4353                        + ( ( flux_t(k) + diss_t(k) ) -                         &
    4354                            ( flux_d    + diss_d    )                           &
    4355                                                                 )              &
    4356                           * drho_air(k) * ddzw(k) &
     4374                         ( flux_r(k) + diss_r(k) -                             &
     4375                           flux_l(k) - diss_l(k) ) * ddx                       &
     4376                       + ( flux_n(k) + diss_n(k) -                             &
     4377                           flux_s(k) - diss_s(k) ) * ddy                       &
     4378                       + ( flux_t(k) + diss_t(k) -                             &
     4379                           flux_d    - diss_d    ) * drho_air(k) * ddzw(k)     &
    43574380                                            ) + sk(k,j,i) * div
    43584381
    43594382#ifndef _OPENACC
    4360                 swap_flux_y_local(k,tn)   = flux_n(k)
    4361                 swap_diss_y_local(k,tn)   = diss_n(k)
    4362                 swap_flux_x_local(k,j,tn) = flux_r(k)
    4363                 swap_diss_x_local(k,j,tn) = diss_r(k)
     4383                swap_flux_y_local(k)   = flux_n(k)
     4384                swap_diss_y_local(k)   = diss_n(k)
     4385                swap_flux_x_local(k,j) = flux_r(k)
     4386                swap_diss_x_local(k,j) = diss_r(k)
    43644387#endif
    43654388
     
    43684391             DO  k = nzb_max_l+1, nzt
    43694392
    4370                 flux_d    = flux_t(k-1)
    4371                 diss_d    = diss_t(k-1)
     4393                flux_d = flux_t(k-1)
     4394                diss_d = diss_t(k-1)
    43724395!
    43734396!--             Calculate the divergence of the velocity field. A respective
     
    43764399                div         =   ( u(k,j,i+1) - u(k,j,i) ) * ddx                &
    43774400                              + ( v(k,j+1,i) - v(k,j,i) ) * ddy                &
    4378                               + ( w(k,j,i) * rho_air_zw(k)                     &
     4401                              + ( w(k,j,i)   * rho_air_zw(k)                   &
    43794402                                - w(k-1,j,i) * rho_air_zw(k-1)                 &
    43804403                                                        ) * drho_air(k) * ddzw(k)
    43814404
    43824405                tend(k,j,i) = tend(k,j,i) - (                                  &
    4383                          ( flux_r(k) + diss_r(k) - swap_flux_x_local(k,j,tn) - &
    4384                            swap_diss_x_local(k,j,tn)            ) * ddx        &
    4385                        + ( flux_n(k) + diss_n(k) - swap_flux_y_local(k,tn)   - &
    4386                            swap_diss_y_local(k,tn)              ) * ddy        &
    4387                        + ( ( flux_t(k) + diss_t(k) ) -                         &
    4388                            ( flux_d    + diss_d    )            )              &
    4389                            * drho_air(k) * ddzw(k)                             &
     4406                         ( flux_r(k) + diss_r(k) -                             &
     4407                           flux_l(k) - diss_l(k) ) * ddx                       &
     4408                       + ( flux_n(k) + diss_n(k) -                             &
     4409                           flux_s(k) - diss_s(k) ) * ddy                       &
     4410                       + ( flux_t(k) + diss_t(k) -                             &
     4411                           flux_d    - diss_d    ) * drho_air(k) * ddzw(k)     &
    43904412                                            ) + sk(k,j,i) * div
    43914413
    43924414#ifndef _OPENACC
    4393                 swap_flux_y_local(k,tn)   = flux_n(k)
    4394                 swap_diss_y_local(k,tn)   = diss_n(k)
    4395                 swap_flux_x_local(k,j,tn) = flux_r(k)
    4396                 swap_diss_x_local(k,j,tn) = diss_r(k)
     4415                swap_flux_y_local(k)   = flux_n(k)
     4416                swap_diss_y_local(k)   = diss_n(k)
     4417                swap_flux_x_local(k,j) = flux_r(k)
     4418                swap_diss_x_local(k,j) = diss_r(k)
    43974419#endif
    43984420            ENDDO
     
    45694591       !$ACC PRIVATE(ibit0_l, ibit1_l, ibit2_l) &
    45704592       !$ACC PRIVATE(ibit3_s, ibit4_s, ibit5_s) &
     4593       !$ACC PRIVATE(nzb_max_l) &
    45714594       !$ACC PRIVATE(ibit6, ibit7, ibit8) &
    45724595       !$ACC PRIVATE(flux_r, diss_r) &
     
    49314954             diss_t(nzb) = 0.0_wp
    49324955             w_comp(nzb) = 0.0_wp
    4933              DO  k = nzb+1, nzb+2
     4956             DO  k = nzb+1, nzb+1
    49344957!
    49354958!--             k index has to be modified near bottom and top, else array
     
    49754998             ENDDO
    49764999
    4977              DO  k = nzb+3, nzt-2
     5000             DO  k = nzb+2, nzt-2
    49785001
    49795002                ibit8 = REAL( IBITS(advc_flags_m(k,j,i),8,1), KIND = wp )
     
    52825305       !$ACC PRIVATE(ibit9_l, ibit10_l, ibit11_l) &
    52835306       !$ACC PRIVATE(ibit12_s, ibit13_s, ibit14_s) &
     5307       !$ACC PRIVATE(nzb_max_l) &
    52845308       !$ACC PRIVATE(flux_r, diss_r) &
    52855309       !$ACC PRIVATE(flux_n, diss_n) &
     
    56375661             diss_t(nzb) = 0.0_wp
    56385662             w_comp(nzb) = 0.0_wp
    5639              DO  k = nzb+1, nzb+2
     5663             DO  k = nzb+1, nzb+1
    56405664!
    56415665!--             k index has to be modified near bottom and top, else array
     
    56815705             ENDDO
    56825706
    5683              DO  k = nzb+3, nzt-2
     5707             DO  k = nzb+2, nzt-2
    56845708                ibit17 = REAL( IBITS(advc_flags_m(k,j,i),17,1), KIND = wp )
    56855709                ibit16 = REAL( IBITS(advc_flags_m(k,j,i),16,1), KIND = wp )
     
    60046028       !$ACC PRIVATE(ibit18_l, ibit19_l, ibit20_l) &
    60056029       !$ACC PRIVATE(ibit21_s, ibit22_s, ibit23_s) &
     6030       !$ACC PRIVATE(nzb_max_l) &
    60066031       !$ACC PRIVATE(flux_r, diss_r) &
    60076032       !$ACC PRIVATE(flux_n, diss_n) &
     
    60276052                 ( bc_dirichlet_s  .OR.  bc_radiation_s )  .AND.  j <= nys + 2  .OR. &
    60286053                 ( bc_dirichlet_n  .OR.  bc_radiation_n )  .AND.  j >= nyn - 2 )  THEN
    6029                 nzb_max_l = nzt
     6054                nzb_max_l = nzt - 1
    60306055             ELSE
    60316056                nzb_max_l = nzb_max
Note: See TracChangeset for help on using the changeset viewer.