Ignore:
Timestamp:
Jan 8, 2019 6:22:50 PM (5 years ago)
Author:
suehring
Message:

Split loops in advec_ws in order to reduce bit queries; Introduce new parameter to better control order degradation of advection scheme at non-cyclic boundaries; Remove setting of Neumann conditions for horizontal velocity variances; Minor bugfix in divergence correction in advection scheme (only has implications at downward-facing wall surfaces)

File:
1 edited

Legend:

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

    r3655 r3661  
    2525! -----------------
    2626! $Id$
     27! - Minor bugfix in divergence correction (only has implications at
     28!   downward-facing wall surfaces)
     29! - Remove setting of Neumann condition for horizontal velocity variances
     30! - Split loops for tendency calculation and divergence correction in order to
     31!   reduce bit queries
     32! - Introduce new parameter nzb_max_l to better control order degradation at
     33!   non-cyclic boundaries
     34!
     35! 3655 2019-01-07 16:51:22Z knoop
    2736! OpenACC port for SPEC
    2837!
     
    247256!> degraded.
    248257!> A divergence correction is applied. It is necessary for topography, since
    249 !> the divergence is not sufficiently reduced, resulting in erroneous fluxes and
    250 !> partly numerical instabilities.
     258!> the divergence is not sufficiently reduced, resulting in erroneous fluxes
     259!> and could lead to numerical instabilities.
    251260!-----------------------------------------------------------------------------!
    252261 MODULE advec_ws
     
    11921201
    11931202       USE control_parameters,                                                 &
    1194            ONLY:  intermediate_timestep_count, u_gtrans, v_gtrans
     1203           ONLY:  bc_dirichlet_l, bc_dirichlet_n, bc_dirichlet_r,              &
     1204                  bc_dirichlet_s, bc_radiation_l, bc_radiation_n,              &
     1205                  bc_radiation_r, bc_radiation_s, intermediate_timestep_count, &
     1206                  u_gtrans, v_gtrans
    11951207
    11961208       USE grid_variables,                                                     &
     
    12181230       CHARACTER (LEN = *), INTENT(IN) ::  sk_char !< string identifier, used for assign fluxes to the correct dimension in the analysis array
    12191231       
    1220        INTEGER(iwp) ::  i     !< grid index along x-direction
    1221        INTEGER(iwp) ::  i_omp !< leftmost index on subdomain, or in case of OpenMP, on thread
    1222        INTEGER(iwp) ::  j     !< grid index along y-direction
    1223        INTEGER(iwp) ::  k     !< grid index along z-direction
    1224        INTEGER(iwp) ::  k_mm  !< k-2 index in disretization, can be modified to avoid segmentation faults
    1225        INTEGER(iwp) ::  k_pp  !< k+2 index in disretization, can be modified to avoid segmentation faults
    1226        INTEGER(iwp) ::  k_ppp !< k+3 index in disretization, can be modified to avoid segmentation faults
    1227        INTEGER(iwp) ::  tn    !< number of OpenMP thread
     1232       INTEGER(iwp) ::  i         !< grid index along x-direction
     1233       INTEGER(iwp) ::  i_omp     !< leftmost index on subdomain, or in case of OpenMP, on thread
     1234       INTEGER(iwp) ::  j         !< grid index along y-direction
     1235       INTEGER(iwp) ::  k         !< grid index along z-direction
     1236       INTEGER(iwp) ::  k_mm      !< k-2 index in disretization, can be modified to avoid segmentation faults
     1237       INTEGER(iwp) ::  k_pp      !< k+2 index in disretization, can be modified to avoid segmentation faults
     1238       INTEGER(iwp) ::  k_ppp     !< k+3 index in disretization, can be modified to avoid segmentation faults
     1239       INTEGER(iwp) ::  nzb_max_l !< index indicating upper bound for order degradation of horizontal advection terms 
     1240       INTEGER(iwp) ::  tn        !< number of OpenMP thread
    12281241       
    12291242       REAL(wp)     ::  ibit0  !< flag indicating 1st-order scheme along x-direction
     
    12611274       REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn,0:threads_per_task-1) ::  swap_diss_x_local !< discretized artificial dissipation at leftward-side of the grid box
    12621275       REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn,0:threads_per_task-1) ::  swap_flux_x_local !< discretized 6th-order flux at leftward-side of the grid box
    1263        
    1264 
     1276!
     1277!--    Used local modified copy of nzb_max (used to degrade order of
     1278!--    discretization) at non-cyclic boundaries. Modify only at relevant points
     1279!--    instead of the entire subdomain. This should lead to better
     1280!--    load balance between boundary and non-boundary PEs.
     1281       IF( ( bc_dirichlet_l  .OR.  bc_radiation_l )  .AND.  i <= nxl + 2  .OR. &
     1282           ( bc_dirichlet_r  .OR.  bc_radiation_r )  .AND.  i >= nxr - 2  .OR. &
     1283           ( bc_dirichlet_s  .OR.  bc_radiation_s )  .AND.  j <= nys + 2  .OR. &
     1284           ( bc_dirichlet_n  .OR.  bc_radiation_n )  .AND.  j >= nyn - 2 )  THEN
     1285          nzb_max_l = nzt
     1286       ELSE
     1287          nzb_max_l = nzb_max
     1288       END IF
    12651289!
    12661290!--    Compute southside fluxes of the respective PE bounds.
     
    12681292!
    12691293!--       Up to the top of the highest topography.
    1270           DO  k = nzb+1, nzb_max
     1294          DO  k = nzb+1, nzb_max_l
    12711295
    12721296             ibit5 = REAL( IBITS(advc_flags_1(k,j-1,i),5,1), KIND = wp )
     
    13291353       IF ( i == i_omp )  THEN
    13301354       
    1331           DO  k = nzb+1, nzb_max
     1355          DO  k = nzb+1, nzb_max_l
    13321356
    13331357             ibit2 = REAL( IBITS(advc_flags_1(k,j,i-1),2,1), KIND = wp )
     
    13681392          ENDDO
    13691393
    1370           DO  k = nzb_max+1, nzt
     1394          DO  k = nzb_max_l+1, nzt
    13711395
    13721396             u_comp                 = u(k,j,i) - u_gtrans + u_stokes_zu(k)
     
    13891413!--    Now compute the fluxes and tendency terms for the horizontal and
    13901414!--    vertical parts up to the top of the highest topography.
    1391        DO  k = nzb+1, nzb_max
     1415       DO  k = nzb+1, nzb_max_l
    13921416!
    13931417!--       Note: It is faster to conduct all multiplications explicitly, e.g.
     
    14691493!--    vertical parts above the top of the highest topography. No degradation
    14701494!--    for the horizontal parts, but for the vertical it is stell needed.
    1471        DO  k = nzb_max+1, nzt
     1495       DO  k = nzb_max_l+1, nzt
    14721496
    14731497          u_comp    = u(k,j,i+1) - u_gtrans + u_stokes_zu(k)
     
    16221646       ENDDO
    16231647
    1624        DO  k = nzb+1, nzt
     1648       DO  k = nzb+1, nzb_max_l
    16251649
    16261650          flux_d    = flux_t(k-1)
    16271651          diss_d    = diss_t(k-1)
    1628 
     1652         
     1653          ibit2 = REAL( IBITS(advc_flags_1(k,j,i),2,1), KIND = wp )
     1654          ibit1 = REAL( IBITS(advc_flags_1(k,j,i),1,1), KIND = wp )
     1655          ibit0 = REAL( IBITS(advc_flags_1(k,j,i),0,1), KIND = wp )
     1656         
     1657          ibit5 = REAL( IBITS(advc_flags_1(k,j,i),5,1), KIND = wp )
     1658          ibit4 = REAL( IBITS(advc_flags_1(k,j,i),4,1), KIND = wp )
     1659          ibit3 = REAL( IBITS(advc_flags_1(k,j,i),3,1), KIND = wp )
     1660         
     1661          ibit8 = REAL( IBITS(advc_flags_1(k,j,i),8,1), KIND = wp )
     1662          ibit7 = REAL( IBITS(advc_flags_1(k,j,i),7,1), KIND = wp )
     1663          ibit6 = REAL( IBITS(advc_flags_1(k,j,i),6,1), KIND = wp )
    16291664!
    16301665!--       Calculate the divergence of the velocity field. A respective
     
    16541689                                         )                                    &     
    16551690                          ) * drho_air(k) * ddzw(k)
     1691
     1692          tend(k,j,i) = tend(k,j,i) - (                                       &
     1693                        ( flux_r(k) + diss_r(k) - swap_flux_x_local(k,j,tn) - &
     1694                          swap_diss_x_local(k,j,tn)            ) * ddx        &
     1695                      + ( flux_n(k) + diss_n(k) - swap_flux_y_local(k,tn)   - &
     1696                          swap_diss_y_local(k,tn)              ) * ddy        &
     1697                      + ( ( flux_t(k) + diss_t(k) ) -                         &
     1698                          ( flux_d    + diss_d    )                           &
     1699                                                    ) * drho_air(k) * ddzw(k) &
     1700                                      ) + sk(k,j,i) * div
     1701
     1702
     1703          swap_flux_y_local(k,tn)   = flux_n(k)
     1704          swap_diss_y_local(k,tn)   = diss_n(k)
     1705          swap_flux_x_local(k,j,tn) = flux_r(k)
     1706          swap_diss_x_local(k,j,tn) = diss_r(k)
     1707
     1708       ENDDO
     1709       
     1710       DO  k = nzb_max_l+1, nzt
     1711
     1712          flux_d    = flux_t(k-1)
     1713          diss_d    = diss_t(k-1)
     1714!
     1715!--       Calculate the divergence of the velocity field. A respective
     1716!--       correction is needed to overcome numerical instabilities introduced
     1717!--       by a not sufficient reduction of divergences near topography.
     1718          div         =   ( u(k,j,i+1) - u(k,j,i) ) * ddx                     &
     1719                        + ( v(k,j+1,i) - v(k,j,i) ) * ddy                     &
     1720                        + ( w(k,j,i) * rho_air_zw(k)                          &
     1721                          - w(k-1,j,i) * rho_air_zw(k-1)                      &
     1722                                                  ) * drho_air(k) * ddzw(k)
    16561723
    16571724          tend(k,j,i) = tend(k,j,i) - (                                       &
     
    17951862    SUBROUTINE advec_u_ws_ij( i, j, i_omp, tn )
    17961863
    1797        USE arrays_3d,                                                         &
    1798            ONLY:  ddzw, diss_l_u, diss_s_u, flux_l_u, flux_s_u, tend, u, v, w,&
     1864       USE arrays_3d,                                                          &
     1865           ONLY:  ddzw, diss_l_u, diss_s_u, flux_l_u, flux_s_u, tend, u, v, w, &
    17991866                  drho_air, rho_air_zw
    18001867
    1801        USE control_parameters,                                                &
    1802            ONLY:  intermediate_timestep_count, u_gtrans, v_gtrans
    1803 
    1804        USE grid_variables,                                                    &
     1868       USE control_parameters,                                                 &
     1869           ONLY:  bc_dirichlet_l, bc_dirichlet_n, bc_dirichlet_r,              &
     1870                  bc_dirichlet_s, bc_radiation_l, bc_radiation_n,              &
     1871                  bc_radiation_r, bc_radiation_s, intermediate_timestep_count, &
     1872                  u_gtrans, v_gtrans
     1873
     1874       USE grid_variables,                                                     &
    18051875           ONLY:  ddx, ddy
    18061876
    1807        USE indices,                                                           &
    1808            ONLY:  nxlu, nys, nzb, nzb_max, nzt, advc_flags_1
     1877       USE indices,                                                            &
     1878           ONLY:  nyn, nys, nxl, nxlu, nxr, nzb, nzb_max, nzt, advc_flags_1
    18091879
    18101880       USE kinds
    18111881
    1812        USE statistics,                                                        &
     1882       USE statistics,                                                         &
    18131883           ONLY:  hom, sums_us2_ws_l, sums_wsus_ws_l, weight_substep
    18141884
    18151885       IMPLICIT NONE
    18161886
    1817        INTEGER(iwp) ::  i      !< grid index along x-direction
    1818        INTEGER(iwp) ::  i_omp  !< leftmost index on subdomain, or in case of OpenMP, on thread
    1819        INTEGER(iwp) ::  j      !< grid index along y-direction
    1820        INTEGER(iwp) ::  k      !< grid index along z-direction
    1821        INTEGER(iwp) ::  k_mm   !< k-2 index in disretization, can be modified to avoid segmentation faults
    1822        INTEGER(iwp) ::  k_pp   !< k+2 index in disretization, can be modified to avoid segmentation faults
    1823        INTEGER(iwp) ::  k_ppp  !< k+3 index in disretization, can be modified to avoid segmentation faults
    1824        INTEGER(iwp) ::  tn     !< number of OpenMP thread
     1887       INTEGER(iwp) ::  i         !< grid index along x-direction
     1888       INTEGER(iwp) ::  i_omp     !< leftmost index on subdomain, or in case of OpenMP, on thread
     1889       INTEGER(iwp) ::  j         !< grid index along y-direction
     1890       INTEGER(iwp) ::  k         !< grid index along z-direction
     1891       INTEGER(iwp) ::  k_mm      !< k-2 index in disretization, can be modified to avoid segmentation faults
     1892       INTEGER(iwp) ::  k_pp      !< k+2 index in disretization, can be modified to avoid segmentation faults
     1893       INTEGER(iwp) ::  k_ppp     !< k+3 index in disretization, can be modified to avoid segmentation faults
     1894       INTEGER(iwp) ::  nzb_max_l !< index indicating upper bound for order degradation of horizontal advection terms 
     1895       INTEGER(iwp) ::  tn        !< number of OpenMP thread
    18251896       
    18261897       REAL(wp)    ::  ibit9    !< flag indicating 1st-order scheme along x-direction
     
    18491920       REAL(wp), DIMENSION(nzb:nzt+1) ::  v_comp !< advection velocity along y
    18501921       REAL(wp), DIMENSION(nzb:nzt+1) ::  w_comp !< advection velocity along z
    1851 
     1922!
     1923!--    Used local modified copy of nzb_max (used to degrade order of
     1924!--    discretization) at non-cyclic boundaries. Modify only at relevant points
     1925!--    instead of the entire subdomain. This should lead to better
     1926!--    load balance between boundary and non-boundary PEs.
     1927       IF( ( bc_dirichlet_l  .OR.  bc_radiation_l )  .AND.  i <= nxl + 2  .OR. &
     1928           ( bc_dirichlet_r  .OR.  bc_radiation_r )  .AND.  i >= nxr - 2  .OR. &
     1929           ( bc_dirichlet_s  .OR.  bc_radiation_s )  .AND.  j <= nys + 2  .OR. &
     1930           ( bc_dirichlet_n  .OR.  bc_radiation_n )  .AND.  j >= nyn - 2 )  THEN
     1931          nzb_max_l = nzt
     1932       ELSE
     1933          nzb_max_l = nzb_max
     1934       END IF
     1935       
    18521936       gu = 2.0_wp * u_gtrans
    18531937       gv = 2.0_wp * v_gtrans
     
    18561940       IF ( j == nys  )  THEN
    18571941       
    1858           DO  k = nzb+1, nzb_max
     1942          DO  k = nzb+1, nzb_max_l
    18591943
    18601944             ibit14 = REAL( IBITS(advc_flags_1(k,j-1,i),14,1), KIND = wp )
     
    18951979          ENDDO
    18961980
    1897           DO  k = nzb_max+1, nzt
     1981          DO  k = nzb_max_l+1, nzt
    18981982
    18991983             v_comp(k)      = v(k,j,i) + v(k,j,i-1) - gv
     
    19141998       IF ( i == i_omp  .OR.  i == nxlu )  THEN
    19151999       
    1916           DO  k = nzb+1, nzb_max
     2000          DO  k = nzb+1, nzb_max_l
    19172001
    19182002             ibit11 = REAL( IBITS(advc_flags_1(k,j,i-1),11,1), KIND = wp )
     
    19532037          ENDDO
    19542038
    1955           DO  k = nzb_max+1, nzt
     2039          DO  k = nzb_max_l+1, nzt
    19562040
    19572041             u_comp_l         = u(k,j,i) + u(k,j,i-1) - gu
     
    19712055!--    Now compute the fluxes tendency terms for the horizontal and
    19722056!--    vertical parts.
    1973        DO  k = nzb+1, nzb_max
     2057       DO  k = nzb+1, nzb_max_l
    19742058
    19752059          ibit11 = REAL( IBITS(advc_flags_1(k,j,i),11,1), KIND = wp )
     
    20442128       ENDDO
    20452129
    2046        DO  k = nzb_max+1, nzt
     2130       DO  k = nzb_max_l+1, nzt
    20472131
    20482132          u_comp(k) = u(k,j,i+1) + u(k,j,i)
     
    22052289       ENDDO
    22062290       
    2207        DO  k = nzb+1, nzt
     2291       DO  k = nzb+1, nzb_max_l
    22082292
    22092293          flux_d    = flux_t(k-1)
    22102294          diss_d    = diss_t(k-1)
     2295         
     2296          ibit11 = REAL( IBITS(advc_flags_1(k,j,i),11,1), KIND = wp )
     2297          ibit10 = REAL( IBITS(advc_flags_1(k,j,i),10,1), KIND = wp )
     2298          ibit9  = REAL( IBITS(advc_flags_1(k,j,i),9,1),  KIND = wp )
     2299         
     2300          ibit14 = REAL( IBITS(advc_flags_1(k,j,i),14,1), KIND = wp )
     2301          ibit13 = REAL( IBITS(advc_flags_1(k,j,i),13,1), KIND = wp )
     2302          ibit12 = REAL( IBITS(advc_flags_1(k,j,i),12,1), KIND = wp )
     2303         
     2304          ibit17 = REAL( IBITS(advc_flags_1(k,j,i),17,1), KIND = wp )
     2305          ibit16 = REAL( IBITS(advc_flags_1(k,j,i),16,1), KIND = wp )
     2306          ibit15 = REAL( IBITS(advc_flags_1(k,j,i),15,1), KIND = wp )
    22112307!
    22122308!--       Calculate the divergence of the velocity field. A respective
     
    22392335                ) * 0.5_wp
    22402336
    2241            tend(k,j,i) = tend(k,j,i) - (                                      &
     2337          tend(k,j,i) = tend(k,j,i) - (                                       &
    22422338                            ( flux_r(k) + diss_r(k)                           &
    22432339                          -   flux_l_u(k,j,tn) - diss_l_u(k,j,tn) ) * ddx     &
     
    22492345                                       ) + div * u(k,j,i)
    22502346
    2251            flux_l_u(k,j,tn) = flux_r(k)
    2252            diss_l_u(k,j,tn) = diss_r(k)
    2253            flux_s_u(k,tn)   = flux_n(k)
    2254            diss_s_u(k,tn)   = diss_n(k)
    2255 !
    2256 !--        Statistical Evaluation of u'u'. The factor has to be applied for
    2257 !--        right evaluation when gallilei_trans = .T. .
    2258            sums_us2_ws_l(k,tn) = sums_us2_ws_l(k,tn)                           &
     2347          flux_l_u(k,j,tn) = flux_r(k)
     2348          diss_l_u(k,j,tn) = diss_r(k)
     2349          flux_s_u(k,tn)   = flux_n(k)
     2350          diss_s_u(k,tn)   = diss_n(k)
     2351!
     2352!--       Statistical Evaluation of u'u'. The factor has to be applied for
     2353!--       right evaluation when gallilei_trans = .T. .
     2354          sums_us2_ws_l(k,tn) = sums_us2_ws_l(k,tn)                            &
    22592355                + ( flux_r(k)                                                  &
    22602356                    * ( u_comp(k) - 2.0_wp * hom(k,1,1,0)                   )  &
     
    22652361                  ) *   weight_substep(intermediate_timestep_count)
    22662362!
    2267 !--        Statistical Evaluation of w'u'.
    2268            sums_wsus_ws_l(k,tn) = sums_wsus_ws_l(k,tn)                         &
     2363!--       Statistical Evaluation of w'u'.
     2364          sums_wsus_ws_l(k,tn) = sums_wsus_ws_l(k,tn)                          &
    22692365                + ( flux_t(k)                                                  &
    22702366                    * ( w_comp(k) - 2.0_wp * hom(k,1,3,0)                   )  &
     
    22752371                  ) *   weight_substep(intermediate_timestep_count)
    22762372       ENDDO
    2277 
    2278        sums_us2_ws_l(nzb,tn) = sums_us2_ws_l(nzb+1,tn)
     2373       
     2374       DO  k = nzb_max_l+1, nzt
     2375
     2376          flux_d    = flux_t(k-1)
     2377          diss_d    = diss_t(k-1)
     2378!
     2379!--       Calculate the divergence of the velocity field. A respective
     2380!--       correction is needed to overcome numerical instabilities introduced
     2381!--       by a not sufficient reduction of divergences near topography.
     2382          div = ( ( u_comp(k)       - ( u(k,j,i)   + u(k,j,i-1) ) ) * ddx      &
     2383               +  ( v_comp(k) + gv  - ( v(k,j,i)   + v(k,j,i-1) ) ) * ddy      &
     2384               +  ( w_comp(k)   * rho_air_zw(k)                                &
     2385                 -  w_comp(k-1) * rho_air_zw(k-1)                              & 
     2386                  ) * drho_air(k) * ddzw(k)                                    &
     2387                ) * 0.5_wp
     2388
     2389          tend(k,j,i) = tend(k,j,i) - (                                        &
     2390                            ( flux_r(k) + diss_r(k)                            &
     2391                          -   flux_l_u(k,j,tn) - diss_l_u(k,j,tn) ) * ddx      &
     2392                          + ( flux_n(k) + diss_n(k)                            &
     2393                          -   flux_s_u(k,tn) - diss_s_u(k,tn)     ) * ddy      &
     2394                          + ( ( flux_t(k) + diss_t(k) )                        &
     2395                          -   ( flux_d    + diss_d    )                        &
     2396                                                    ) * drho_air(k) * ddzw(k)  &
     2397                                       ) + div * u(k,j,i)
     2398
     2399          flux_l_u(k,j,tn) = flux_r(k)
     2400          diss_l_u(k,j,tn) = diss_r(k)
     2401          flux_s_u(k,tn)   = flux_n(k)
     2402          diss_s_u(k,tn)   = diss_n(k)
     2403!
     2404!--       Statistical Evaluation of u'u'. The factor has to be applied for
     2405!--       right evaluation when gallilei_trans = .T. .
     2406          sums_us2_ws_l(k,tn) = sums_us2_ws_l(k,tn)                            &
     2407                + ( flux_r(k)                                                  &
     2408                    * ( u_comp(k) - 2.0_wp * hom(k,1,1,0)                   )  &
     2409                    / ( u_comp(k) - gu + SIGN( 1.0E-20_wp, u_comp(k) - gu ) )  &
     2410                  + diss_r(k)                                                  &
     2411                    *   ABS( u_comp(k) - 2.0_wp * hom(k,1,1,0)              )  &
     2412                    / ( ABS( u_comp(k) - gu ) + 1.0E-20_wp                  )  &
     2413                  ) *   weight_substep(intermediate_timestep_count)
     2414!
     2415!--       Statistical Evaluation of w'u'.
     2416          sums_wsus_ws_l(k,tn) = sums_wsus_ws_l(k,tn)                          &
     2417                + ( flux_t(k)                                                  &
     2418                    * ( w_comp(k) - 2.0_wp * hom(k,1,3,0)                   )  &
     2419                    / ( w_comp(k) + SIGN( 1.0E-20_wp, w_comp(k) )           )  &
     2420                  + diss_t(k)                                                  &
     2421                    *   ABS( w_comp(k) - 2.0_wp * hom(k,1,3,0)              )  &
     2422                    / ( ABS( w_comp(k) ) + 1.0E-20_wp                       )  &
     2423                  ) *   weight_substep(intermediate_timestep_count)
     2424       ENDDO
    22792425
    22802426
     
    22962442
    22972443       USE control_parameters,                                                 &
    2298            ONLY:  intermediate_timestep_count, u_gtrans, v_gtrans
     2444           ONLY:  bc_dirichlet_l, bc_dirichlet_n, bc_dirichlet_r,              &
     2445                  bc_dirichlet_s, bc_radiation_l, bc_radiation_n,              &
     2446                  bc_radiation_r, bc_radiation_s, intermediate_timestep_count, &
     2447                  u_gtrans, v_gtrans
    22992448
    23002449       USE grid_variables,                                                     &
     
    23022451
    23032452       USE indices,                                                            &
    2304            ONLY:  nysv, nzb, nzb_max, nzt, advc_flags_1
     2453           ONLY:  nyn, nys, nysv, nxl, nxr, nzb, nzb_max, nzt, advc_flags_1
    23052454
    23062455       USE kinds
     
    23112460       IMPLICIT NONE
    23122461
    2313        INTEGER(iwp)  ::  i      !< grid index along x-direction
    2314        INTEGER(iwp)  ::  i_omp  !< leftmost index on subdomain, or in case of OpenMP, on thread
    2315        INTEGER(iwp)  ::  j      !< grid index along y-direction
    2316        INTEGER(iwp)  ::  k      !< grid index along z-direction
    2317        INTEGER(iwp)  ::  k_mm   !< k-2 index in disretization, can be modified to avoid segmentation faults
    2318        INTEGER(iwp)  ::  k_pp   !< k+2 index in disretization, can be modified to avoid segmentation faults
    2319        INTEGER(iwp)  ::  k_ppp  !< k+3 index in disretization, can be modified to avoid segmentation faults
    2320        INTEGER(iwp)  ::  tn     !< number of OpenMP thread
     2462       INTEGER(iwp)  ::  i         !< grid index along x-direction
     2463       INTEGER(iwp)  ::  i_omp     !< leftmost index on subdomain, or in case of OpenMP, on thread
     2464       INTEGER(iwp)  ::  j         !< grid index along y-direction
     2465       INTEGER(iwp)  ::  k         !< grid index along z-direction
     2466       INTEGER(iwp)  ::  k_mm      !< k-2 index in disretization, can be modified to avoid segmentation faults
     2467       INTEGER(iwp)  ::  k_pp      !< k+2 index in disretization, can be modified to avoid segmentation faults
     2468       INTEGER(iwp)  ::  k_ppp     !< k+3 index in disretization, can be modified to avoid segmentation faults
     2469       INTEGER(iwp)  ::  nzb_max_l !< index indicating upper bound for order degradation of horizontal advection terms 
     2470       INTEGER(iwp)  ::  tn        !< number of OpenMP thread
    23212471       
    23222472       REAL(wp)      ::  ibit18   !< flag indicating 1st-order scheme along x-direction
     
    23452495       REAL(wp), DIMENSION(nzb:nzt+1)  ::  v_comp !< advection velocity along y
    23462496       REAL(wp), DIMENSION(nzb:nzt+1)  ::  w_comp !< advection velocity along z
    2347 
     2497!
     2498!--    Used local modified copy of nzb_max (used to degrade order of
     2499!--    discretization) at non-cyclic boundaries. Modify only at relevant points
     2500!--    instead of the entire subdomain. This should lead to better
     2501!--    load balance between boundary and non-boundary PEs.
     2502       IF( ( bc_dirichlet_l  .OR.  bc_radiation_l )  .AND.  i <= nxl + 2  .OR. &
     2503           ( bc_dirichlet_r  .OR.  bc_radiation_r )  .AND.  i >= nxr - 2  .OR. &
     2504           ( bc_dirichlet_s  .OR.  bc_radiation_s )  .AND.  j <= nys + 2  .OR. &
     2505           ( bc_dirichlet_n  .OR.  bc_radiation_n )  .AND.  j >= nyn - 2 )  THEN
     2506          nzb_max_l = nzt
     2507       ELSE
     2508          nzb_max_l = nzb_max
     2509       END IF
     2510       
    23482511       gu = 2.0_wp * u_gtrans
    23492512       gv = 2.0_wp * v_gtrans
     
    23532516       IF ( i == i_omp )  THEN
    23542517
    2355           DO  k = nzb+1, nzb_max
     2518          DO  k = nzb+1, nzb_max_l
    23562519
    23572520             ibit20 = REAL( IBITS(advc_flags_1(k,j,i-1),20,1), KIND = wp )
     
    23922555          ENDDO
    23932556
    2394           DO  k = nzb_max+1, nzt
     2557          DO  k = nzb_max_l+1, nzt
    23952558
    23962559             u_comp(k)           = u(k,j-1,i) + u(k,j,i) - gu
     
    24112574       IF ( j == nysv )  THEN
    24122575       
    2413           DO  k = nzb+1, nzb_max
     2576          DO  k = nzb+1, nzb_max_l
    24142577
    24152578             ibit23 = REAL( IBITS(advc_flags_1(k,j-1,i),23,1), KIND = wp )
     
    24502613          ENDDO
    24512614
    2452           DO  k = nzb_max+1, nzt
     2615          DO  k = nzb_max_l+1, nzt
    24532616
    24542617             v_comp_l       = v(k,j,i) + v(k,j-1,i) - gv
     
    24682631!--    Now compute the fluxes and tendency terms for the horizontal and
    24692632!--    verical parts.
    2470        DO  k = nzb+1, nzb_max
     2633       DO  k = nzb+1, nzb_max_l
    24712634
    24722635          ibit20 = REAL( IBITS(advc_flags_1(k,j,i),20,1), KIND = wp )
     
    25422705       ENDDO
    25432706
    2544        DO  k = nzb_max+1, nzt
     2707       DO  k = nzb_max_l+1, nzt
    25452708
    25462709          u_comp(k) = u(k,j-1,i+1) + u(k,j,i+1) - gu
     
    27062869       ENDDO
    27072870       
    2708        DO  k = nzb+1, nzt
     2871       DO  k = nzb+1, nzb_max_l
    27092872
    27102873          flux_d    = flux_t(k-1)
    27112874          diss_d    = diss_t(k-1)
     2875         
     2876          ibit20 = REAL( IBITS(advc_flags_1(k,j,i),20,1), KIND = wp )
     2877          ibit19 = REAL( IBITS(advc_flags_1(k,j,i),19,1), KIND = wp )
     2878          ibit18 = REAL( IBITS(advc_flags_1(k,j,i),18,1), KIND = wp )
     2879         
     2880          ibit23 = REAL( IBITS(advc_flags_1(k,j,i),23,1), KIND = wp )
     2881          ibit22 = REAL( IBITS(advc_flags_1(k,j,i),22,1), KIND = wp )
     2882          ibit21 = REAL( IBITS(advc_flags_1(k,j,i),21,1), KIND = wp )
     2883         
     2884          ibit26 = REAL( IBITS(advc_flags_1(k,j,i),26,1), KIND = wp )
     2885          ibit25 = REAL( IBITS(advc_flags_1(k,j,i),25,1), KIND = wp )
     2886          ibit24 = REAL( IBITS(advc_flags_1(k,j,i),24,1), KIND = wp )
    27122887!
    27132888!--       Calculate the divergence of the velocity field. A respective
     
    27522927                                      ) + v(k,j,i) * div
    27532928
    2754            flux_l_v(k,j,tn) = flux_r(k)
    2755            diss_l_v(k,j,tn) = diss_r(k)
    2756            flux_s_v(k,tn)   = flux_n(k)
    2757            diss_s_v(k,tn)   = diss_n(k)
    2758 !
    2759 !--        Statistical Evaluation of v'v'. The factor has to be applied for
    2760 !--        right evaluation when gallilei_trans = .T. .
    2761            sums_vs2_ws_l(k,tn) = sums_vs2_ws_l(k,tn)                           &
     2929          flux_l_v(k,j,tn) = flux_r(k)
     2930          diss_l_v(k,j,tn) = diss_r(k)
     2931          flux_s_v(k,tn)   = flux_n(k)
     2932          diss_s_v(k,tn)   = diss_n(k)
     2933!
     2934!--       Statistical Evaluation of v'v'. The factor has to be applied for
     2935!--       right evaluation when gallilei_trans = .T. .
     2936          sums_vs2_ws_l(k,tn) = sums_vs2_ws_l(k,tn)                            &
    27622937                + ( flux_n(k)                                                  &
    27632938                    * ( v_comp(k) - 2.0_wp * hom(k,1,2,0)                   )  &
     
    27682943                  ) *   weight_substep(intermediate_timestep_count)
    27692944!
    2770 !--        Statistical Evaluation of w'u'.
    2771            sums_wsvs_ws_l(k,tn) = sums_wsvs_ws_l(k,tn)                         &
     2945!--       Statistical Evaluation of w'u'.
     2946          sums_wsvs_ws_l(k,tn) = sums_wsvs_ws_l(k,tn)                          &
    27722947                + ( flux_t(k)                                                  &
    27732948                    * ( w_comp(k) - 2.0_wp * hom(k,1,3,0)                   )  &
     
    27792954
    27802955       ENDDO
    2781        sums_vs2_ws_l(nzb,tn) = sums_vs2_ws_l(nzb+1,tn)
     2956       
     2957       DO  k = nzb_max_l+1, nzt
     2958
     2959          flux_d    = flux_t(k-1)
     2960          diss_d    = diss_t(k-1)
     2961!
     2962!--       Calculate the divergence of the velocity field. A respective
     2963!--       correction is needed to overcome numerical instabilities introduced
     2964!--       by a not sufficient reduction of divergences near topography.
     2965          div = ( ( u_comp(k) + gu - ( u(k,j-1,i) + u(k,j,i)   ) ) * ddx       &
     2966               +  ( v_comp(k)      - ( v(k,j,i)   + v(k,j-1,i) ) ) * ddy       &
     2967               +  ( w_comp(k)   * rho_air_zw(k)                                &
     2968                 -  w_comp(k-1) * rho_air_zw(k-1)                              &
     2969                  ) * drho_air(k) * ddzw(k)                                    &
     2970                ) * 0.5_wp
     2971
     2972          tend(k,j,i) = tend(k,j,i) - (                                        &
     2973                         ( flux_r(k) + diss_r(k)                               &
     2974                       -   flux_l_v(k,j,tn) - diss_l_v(k,j,tn)   ) * ddx       &
     2975                       + ( flux_n(k) + diss_n(k)                               &
     2976                       -   flux_s_v(k,tn) - diss_s_v(k,tn)       ) * ddy       &
     2977                       + ( ( flux_t(k) + diss_t(k) )                           &
     2978                       -   ( flux_d    + diss_d    )                           &
     2979                                                   ) * drho_air(k) * ddzw(k)   &
     2980                                      ) + v(k,j,i) * div
     2981
     2982          flux_l_v(k,j,tn) = flux_r(k)
     2983          diss_l_v(k,j,tn) = diss_r(k)
     2984          flux_s_v(k,tn)   = flux_n(k)
     2985          diss_s_v(k,tn)   = diss_n(k)
     2986!
     2987!--       Statistical Evaluation of v'v'. The factor has to be applied for
     2988!--       right evaluation when gallilei_trans = .T. .
     2989          sums_vs2_ws_l(k,tn) = sums_vs2_ws_l(k,tn)                            &
     2990                + ( flux_n(k)                                                  &
     2991                    * ( v_comp(k) - 2.0_wp * hom(k,1,2,0)                   )  &
     2992                    / ( v_comp(k) - gv + SIGN( 1.0E-20_wp, v_comp(k) - gv ) )  &
     2993                  + diss_n(k)                                                  &
     2994                    *   ABS( v_comp(k) - 2.0_wp * hom(k,1,2,0)              )  &
     2995                    / ( ABS( v_comp(k) - gv ) + 1.0E-20_wp                  )  &
     2996                  ) *   weight_substep(intermediate_timestep_count)
     2997!
     2998!--       Statistical Evaluation of w'u'.
     2999          sums_wsvs_ws_l(k,tn) = sums_wsvs_ws_l(k,tn)                          &
     3000                + ( flux_t(k)                                                  &
     3001                    * ( w_comp(k) - 2.0_wp * hom(k,1,3,0)                   )  &
     3002                    / ( w_comp(k) + SIGN( 1.0E-20_wp, w_comp(k) )           )  &
     3003                  + diss_t(k)                                                  &
     3004                    *   ABS( w_comp(k) - 2.0_wp * hom(k,1,3,0)              )  &
     3005                    / ( ABS( w_comp(k) ) + 1.0E-20_wp                       )  &
     3006                  ) *   weight_substep(intermediate_timestep_count)
     3007
     3008       ENDDO
    27823009
    27833010
     
    27933020    SUBROUTINE advec_w_ws_ij( i, j, i_omp, tn )
    27943021
    2795        USE arrays_3d,                                                         &
    2796            ONLY:  ddzu, diss_l_w, diss_s_w, flux_l_w, flux_s_w, tend, u, v, w,&
     3022       USE arrays_3d,                                                          &
     3023           ONLY:  ddzu, diss_l_w, diss_s_w, flux_l_w, flux_s_w, tend, u, v, w, &
    27973024                  drho_air_zw, rho_air
    27983025
    2799        USE control_parameters,                                                &
    2800            ONLY:  intermediate_timestep_count, u_gtrans, v_gtrans
    2801 
    2802        USE grid_variables,                                                    &
     3026       USE control_parameters,                                                 &
     3027           ONLY:  bc_dirichlet_l, bc_dirichlet_n, bc_dirichlet_r,              &
     3028                  bc_dirichlet_s, bc_radiation_l, bc_radiation_n,              &
     3029                  bc_radiation_r, bc_radiation_s,intermediate_timestep_count,  &
     3030                  u_gtrans, v_gtrans
     3031
     3032       USE grid_variables,                                                     &
    28033033           ONLY:  ddx, ddy
    28043034
    2805        USE indices,                                                           &
    2806            ONLY:  nys, nzb, nzb_max, nzt, advc_flags_1, advc_flags_2
     3035       USE indices,                                                            &
     3036           ONLY:  nyn, nys, nxl, nxr, nzb, nzb_max, nzt, advc_flags_1,         &
     3037                  advc_flags_2
    28073038
    28083039       USE kinds
    28093040       
    2810        USE statistics,                                                        &
     3041       USE statistics,                                                         &
    28113042           ONLY:  hom, sums_ws2_ws_l, weight_substep
    28123043
    28133044       IMPLICIT NONE
    28143045
    2815        INTEGER(iwp) ::  i      !< grid index along x-direction
    2816        INTEGER(iwp) ::  i_omp  !< leftmost index on subdomain, or in case of OpenMP, on thread
    2817        INTEGER(iwp) ::  j      !< grid index along y-direction
    2818        INTEGER(iwp) ::  k      !< grid index along z-direction
    2819        INTEGER(iwp) ::  k_mm   !< k-2 index in disretization, can be modified to avoid segmentation faults
    2820        INTEGER(iwp) ::  k_pp   !< k+2 index in disretization, can be modified to avoid segmentation faults
    2821        INTEGER(iwp) ::  k_ppp  !< k+3 index in disretization, can be modified to avoid segmentation faults
    2822        INTEGER(iwp) ::  tn     !< number of OpenMP thread
     3046       INTEGER(iwp) ::  i         !< grid index along x-direction
     3047       INTEGER(iwp) ::  i_omp     !< leftmost index on subdomain, or in case of OpenMP, on thread
     3048       INTEGER(iwp) ::  j         !< grid index along y-direction
     3049       INTEGER(iwp) ::  k         !< grid index along z-direction
     3050       INTEGER(iwp) ::  k_mm      !< k-2 index in disretization, can be modified to avoid segmentation faults
     3051       INTEGER(iwp) ::  k_pp      !< k+2 index in disretization, can be modified to avoid segmentation faults
     3052       INTEGER(iwp) ::  k_ppp     !< k+3 index in disretization, can be modified to avoid segmentation faults
     3053       INTEGER(iwp) ::  nzb_max_l !< index indicating upper bound for order degradation of horizontal advection terms 
     3054       INTEGER(iwp) ::  tn        !< number of OpenMP thread
    28233055       
    28243056       REAL(wp)    ::  ibit27  !< flag indicating 1st-order scheme along x-direction
     
    28463078       REAL(wp), DIMENSION(nzb:nzt+1)  ::  v_comp !< advection velocity along y
    28473079       REAL(wp), DIMENSION(nzb:nzt+1)  ::  w_comp !< advection velocity along z
    2848 
     3080!
     3081!--    Used local modified copy of nzb_max (used to degrade order of
     3082!--    discretization) at non-cyclic boundaries. Modify only at relevant points
     3083!--    instead of the entire subdomain. This should lead to better
     3084!--    load balance between boundary and non-boundary PEs.
     3085       IF( ( bc_dirichlet_l  .OR.  bc_radiation_l )  .AND.  i <= nxl + 2  .OR. &
     3086           ( bc_dirichlet_r  .OR.  bc_radiation_r )  .AND.  i >= nxr - 2  .OR. &
     3087           ( bc_dirichlet_s  .OR.  bc_radiation_s )  .AND.  j <= nys + 2  .OR. &
     3088           ( bc_dirichlet_n  .OR.  bc_radiation_n )  .AND.  j >= nyn - 2 )  THEN
     3089          nzb_max_l = nzt
     3090       ELSE
     3091          nzb_max_l = nzb_max
     3092       END IF
     3093       
    28493094       gu = 2.0_wp * u_gtrans
    28503095       gv = 2.0_wp * v_gtrans
    2851 
    28523096!
    28533097!--    Compute southside fluxes for the respective boundary.
    28543098       IF ( j == nys )  THEN
    28553099
    2856           DO  k = nzb+1, nzb_max
     3100          DO  k = nzb+1, nzb_max_l
    28573101             ibit32 = REAL( IBITS(advc_flags_2(k,j-1,i),0,1),  KIND = wp )
    28583102             ibit31 = REAL( IBITS(advc_flags_1(k,j-1,i),31,1), KIND = wp )
     
    28923136          ENDDO
    28933137
    2894           DO  k = nzb_max+1, nzt
     3138          DO  k = nzb_max_l+1, nzt
    28953139
    28963140             v_comp(k)      = v(k+1,j,i) + v(k,j,i) - gv
     
    29113155       IF ( i == i_omp ) THEN
    29123156
    2913           DO  k = nzb+1, nzb_max
     3157          DO  k = nzb+1, nzb_max_l
    29143158
    29153159             ibit29 = REAL( IBITS(advc_flags_1(k,j,i-1),29,1), KIND = wp )
     
    29503194          ENDDO
    29513195
    2952           DO  k = nzb_max+1, nzt
     3196          DO  k = nzb_max_l+1, nzt
    29533197
    29543198             u_comp(k)        = u(k+1,j,i) + u(k,j,i) - gu
     
    29683212!--    Now compute the fluxes and tendency terms for the horizontal
    29693213!--    and vertical parts.
    2970        DO  k = nzb+1, nzb_max
     3214       DO  k = nzb+1, nzb_max_l
    29713215
    29723216          ibit29 = REAL( IBITS(advc_flags_1(k,j,i),29,1), KIND = wp )
     
    30413285       ENDDO
    30423286
    3043        DO  k = nzb_max+1, nzt
     3287       DO  k = nzb_max_l+1, nzt
    30443288
    30453289          u_comp(k) = u(k+1,j,i+1) + u(k,j,i+1) - gu
     
    32063450       ENDDO
    32073451       
    3208        DO  k = nzb+1, nzt
     3452       DO  k = nzb+1, nzb_max_l
    32093453
    32103454          flux_d    = flux_t(k-1)
    32113455          diss_d    = diss_t(k-1)
     3456         
     3457          ibit29 = REAL( IBITS(advc_flags_1(k,j,i),29,1), KIND = wp )
     3458          ibit28 = REAL( IBITS(advc_flags_1(k,j,i),28,1), KIND = wp )
     3459          ibit27 = REAL( IBITS(advc_flags_1(k,j,i),27,1), KIND = wp )
     3460         
     3461          ibit32 = REAL( IBITS(advc_flags_2(k,j,i),0,1),  KIND = wp )
     3462          ibit31 = REAL( IBITS(advc_flags_1(k,j,i),31,1), KIND = wp )
     3463          ibit30 = REAL( IBITS(advc_flags_1(k,j,i),30,1), KIND = wp )
     3464         
     3465          ibit35 = REAL( IBITS(advc_flags_2(k,j,i),3,1), KIND = wp )
     3466          ibit34 = REAL( IBITS(advc_flags_2(k,j,i),2,1), KIND = wp )
     3467          ibit33 = REAL( IBITS(advc_flags_2(k,j,i),1,1), KIND = wp )
    32123468!
    32133469!--       Calculate the divergence of the velocity field. A respective
     
    32673523
    32683524       ENDDO
    3269 
     3525       
     3526       DO  k = nzb_max_l+1, nzt
     3527
     3528          flux_d    = flux_t(k-1)
     3529          diss_d    = diss_t(k-1)
     3530!
     3531!--       Calculate the divergence of the velocity field. A respective
     3532!--       correction is needed to overcome numerical instabilities introduced
     3533!--       by a not sufficient reduction of divergences near topography.
     3534          div = ( ( u_comp(k) + gu - ( u(k+1,j,i) + u(k,j,i)   ) ) * ddx       &
     3535              +   ( v_comp(k) + gv - ( v(k+1,j,i) + v(k,j,i)   ) ) * ddy       &
     3536              +   ( w_comp(k)               * rho_air(k+1)                     &
     3537                - ( w(k,j,i) + w(k-1,j,i) ) * rho_air(k)                       &
     3538                  ) * drho_air_zw(k) * ddzu(k+1)                               &
     3539                ) * 0.5_wp
     3540
     3541          tend(k,j,i) = tend(k,j,i) - (                                        &
     3542                      ( flux_r(k) + diss_r(k)                                  &
     3543                    -   flux_l_w(k,j,tn) - diss_l_w(k,j,tn)   ) * ddx          &
     3544                    + ( flux_n(k) + diss_n(k)                                  &
     3545                    -   flux_s_w(k,tn) - diss_s_w(k,tn)       ) * ddy          &
     3546                    + ( ( flux_t(k) + diss_t(k) )                              &
     3547                    -   ( flux_d    + diss_d    )                              &
     3548                                              ) * drho_air_zw(k) * ddzu(k+1)   &
     3549                                      ) + div * w(k,j,i)
     3550
     3551          flux_l_w(k,j,tn) = flux_r(k)
     3552          diss_l_w(k,j,tn) = diss_r(k)
     3553          flux_s_w(k,tn)   = flux_n(k)
     3554          diss_s_w(k,tn)   = diss_n(k)
     3555!
     3556!--       Statistical Evaluation of w'w'.
     3557          sums_ws2_ws_l(k,tn)  = sums_ws2_ws_l(k,tn)                           &
     3558                      + ( flux_t(k)                                            &
     3559                       * ( w_comp(k) - 2.0_wp * hom(k,1,3,0)                )  &
     3560                       / ( w_comp(k) + SIGN( 1.0E-20_wp, w_comp(k) )        )  &
     3561                        + diss_t(k)                                            &
     3562                       *   ABS( w_comp(k) - 2.0_wp * hom(k,1,3,0)           )  &
     3563                       / ( ABS( w_comp(k) ) + 1.0E-20_wp                    )  &
     3564                        ) *   weight_substep(intermediate_timestep_count)
     3565
     3566       ENDDO
    32703567
    32713568    END SUBROUTINE advec_w_ws_ij
     
    32843581
    32853582       USE control_parameters,                                                 &
    3286            ONLY:  intermediate_timestep_count, u_gtrans, v_gtrans
     3583           ONLY:  bc_dirichlet_l, bc_dirichlet_n, bc_dirichlet_r,              &
     3584                  bc_dirichlet_s, bc_radiation_l, bc_radiation_n,              &
     3585                  bc_radiation_r, bc_radiation_s, intermediate_timestep_count, &
     3586                  u_gtrans, v_gtrans
    32873587
    32883588       USE grid_variables,                                                     &
     
    33103610       INTEGER(iwp) ::  sk_num !< integer identifier, used for assign fluxes to the correct dimension in the analysis array
    33113611       
    3312        INTEGER(iwp) ::  i      !< grid index along x-direction
    3313        INTEGER(iwp) ::  j      !< grid index along y-direction
    3314        INTEGER(iwp) ::  k      !< grid index along z-direction
    3315        INTEGER(iwp) ::  k_mm   !< k-2 index in disretization, can be modified to avoid segmentation faults
    3316        INTEGER(iwp) ::  k_pp   !< k+2 index in disretization, can be modified to avoid segmentation faults
    3317        INTEGER(iwp) ::  k_ppp  !< k+3 index in disretization, can be modified to avoid segmentation faults
    3318        INTEGER(iwp) ::  tn = 0 !< number of OpenMP thread
     3612       INTEGER(iwp) ::  i         !< grid index along x-direction
     3613       INTEGER(iwp) ::  j         !< grid index along y-direction
     3614       INTEGER(iwp) ::  k         !< grid index along z-direction
     3615       INTEGER(iwp) ::  k_mm      !< k-2 index in disretization, can be modified to avoid segmentation faults
     3616       INTEGER(iwp) ::  k_pp      !< k+2 index in disretization, can be modified to avoid segmentation faults
     3617       INTEGER(iwp) ::  k_ppp     !< k+3 index in disretization, can be modified to avoid segmentation faults
     3618       INTEGER(iwp) ::  nzb_max_l !< index indicating upper bound for order degradation of horizontal advection terms
     3619       INTEGER(iwp) ::  tn = 0    !< number of OpenMP thread
    33193620       
    33203621!       
     
    33783679#endif
    33793680       
     3681!
     3682!--    Set local version of nzb_max. At non-cyclic boundaries the order of the
     3683!--    advection need to be degraded near the boundary. Please note, in contrast
     3684!--    to the cache-optimized routines, nzb_max_l is set constantly for the
     3685!--    entire subdomain, in order to avoid unsymmetric loops which might be
     3686!--    an issue for GPUs.
     3687       IF( bc_dirichlet_l  .OR.  bc_radiation_l  .OR.                          &
     3688           bc_dirichlet_r  .OR.  bc_radiation_r  .OR.                          &
     3689           bc_dirichlet_s  .OR.  bc_radiation_s  .OR.                          &
     3690           bc_dirichlet_n  .OR.  bc_radiation_n )  THEN
     3691          nzb_max_l = nzt
     3692       ELSE
     3693          nzb_max_l = nzb_max
     3694       END IF
     3695       
    33803696       SELECT CASE ( sk_char )
    33813697
     
    34073723       DO  j = nys, nyn
    34083724
    3409           DO  k = nzb+1, nzb_max
     3725          DO  k = nzb+1, nzb_max_l
    34103726
    34113727             ibit2 = REAL( IBITS(advc_flags_1(k,j,i-1),2,1), KIND = wp )
     
    34463762          ENDDO
    34473763
    3448           DO  k = nzb_max+1, nzt
     3764          DO  k = nzb_max_l+1, nzt
    34493765
    34503766             u_comp                 = u(k,j,i) - u_gtrans + u_stokes_zu(k)
     
    34803796       !$ACC PRESENT(drho_air, rho_air_zw, ddzw) &
    34813797       !$ACC PRESENT(tend) &
    3482        !$ACC PRESENT(hom(nzb+1:nzb_max,1,1:3,0)) &
     3798       !$ACC PRESENT(hom(nzb+1:nzb_max_l,1,1:3,0)) &
    34833799       !$ACC PRESENT(weight_substep(intermediate_timestep_count)) &
    34843800       !$ACC PRESENT(sums_wspts_ws_l, sums_wssas_ws_l) &
     
    34913807#ifndef _OPENACC
    34923808          j = nys
    3493           DO  k = nzb+1, nzb_max
     3809          DO  k = nzb+1, nzb_max_l
    34943810
    34953811             ibit5 = REAL( IBITS(advc_flags_1(k,j-1,i),5,1), KIND = wp )
     
    35313847!
    35323848!--       Above to the top of the highest topography. No degradation necessary.
    3533           DO  k = nzb_max+1, nzt
     3849          DO  k = nzb_max_l+1, nzt
    35343850
    35353851             v_comp               = v(k,j,i) - v_gtrans + v_stokes_zu(k)
     
    35533869             diss_d    = 0.0_wp
    35543870
    3555              DO  k = nzb+1, nzb_max
     3871             DO  k = nzb+1, nzb_max_l
    35563872
    35573873                ibit2 = REAL( IBITS(advc_flags_1(k,j,i),2,1), KIND = wp )
     
    38954211             ENDDO
    38964212
    3897              DO  k = nzb_max+1, nzt
     4213             DO  k = nzb_max_l+1, nzt
    38984214
    38994215                u_comp = u(k,j,i+1) - u_gtrans + u_stokes_zu(k)
     
    39684284
    39694285
    3970                 flux_t = w(k,j,i) * rho_air_zw(k) * (                      &
     4286                flux_t = w(k,j,i) * rho_air_zw(k) * (                         &
    39714287                           ( 37.0_wp * ibit8 * adv_sca_5                      &
    39724288                        +     7.0_wp * ibit7 * adv_sca_3                      &
     
    39824298                                       )
    39834299
    3984                 diss_t = -ABS( w(k,j,i) ) * rho_air_zw(k) * (              &
     4300                diss_t = -ABS( w(k,j,i) ) * rho_air_zw(k) * (                 &
    39854301                           ( 10.0_wp * ibit8 * adv_sca_5                      &
    39864302                        +     3.0_wp * ibit7 * adv_sca_3                      &
     
    41404456
    41414457       USE control_parameters,                                                 &
    4142            ONLY:  intermediate_timestep_count, u_gtrans, v_gtrans
     4458           ONLY:  bc_dirichlet_l, bc_dirichlet_n, bc_dirichlet_r,              &
     4459                  bc_dirichlet_s, bc_radiation_l, bc_radiation_n,              &
     4460                  bc_radiation_r, bc_radiation_s, intermediate_timestep_count, &
     4461                  u_gtrans, v_gtrans
    41434462
    41444463       USE grid_variables,                                                     &
     
    41554474       IMPLICIT NONE
    41564475
    4157        INTEGER(iwp) ::  i      !< grid index along x-direction
    4158        INTEGER(iwp) ::  j      !< grid index along y-direction
    4159        INTEGER(iwp) ::  k      !< grid index along z-direction
    4160        INTEGER(iwp) ::  k_mm   !< k-2 index in disretization, can be modified to avoid segmentation faults
    4161        INTEGER(iwp) ::  k_pp   !< k+2 index in disretization, can be modified to avoid segmentation faults
    4162        INTEGER(iwp) ::  k_ppp  !< k+3 index in disretization, can be modified to avoid segmentation faults
    4163        INTEGER(iwp) ::  tn = 0 !< number of OpenMP thread
     4476       INTEGER(iwp) ::  i         !< grid index along x-direction
     4477       INTEGER(iwp) ::  j         !< grid index along y-direction
     4478       INTEGER(iwp) ::  k         !< grid index along z-direction
     4479       INTEGER(iwp) ::  k_mm      !< k-2 index in disretization, can be modified to avoid segmentation faults
     4480       INTEGER(iwp) ::  k_pp      !< k+2 index in disretization, can be modified to avoid segmentation faults
     4481       INTEGER(iwp) ::  k_ppp     !< k+3 index in disretization, can be modified to avoid segmentation faults
     4482       INTEGER(iwp) ::  nzb_max_l !< index indicating upper bound for order degradation of horizontal advection terms
     4483       INTEGER(iwp) ::  tn = 0    !< number of OpenMP thread
    41644484       
    41654485       REAL(wp)    ::  ibit9  !< flag indicating 1st-order scheme along x-direction
     
    42174537       REAL(wp)    ::  u_comp_l !<
    42184538#endif
     4539!
     4540!--    Set local version of nzb_max. At non-cyclic boundaries the order of the
     4541!--    advection need to be degraded near the boundary. Please note, in contrast
     4542!--    to the cache-optimized routines, nzb_max_l is set constantly for the
     4543!--    entire subdomain, in order to avoid unsymmetric loops which might be
     4544!--    an issue for GPUs.
     4545       IF( bc_dirichlet_l  .OR.  bc_radiation_l  .OR.                          &
     4546           bc_dirichlet_r  .OR.  bc_radiation_r  .OR.                          &
     4547           bc_dirichlet_s  .OR.  bc_radiation_s  .OR.                          &
     4548           bc_dirichlet_n  .OR.  bc_radiation_n )  THEN
     4549          nzb_max_l = nzt
     4550       ELSE
     4551          nzb_max_l = nzb_max
     4552       END IF
    42194553 
    42204554       gu = 2.0_wp * u_gtrans
     
    42264560       i = nxlu
    42274561       DO  j = nys, nyn
    4228           DO  k = nzb+1, nzb_max
     4562          DO  k = nzb+1, nzb_max_l
    42294563
    42304564             ibit11 = REAL( IBITS(advc_flags_1(k,j,i-1),11,1), KIND = wp )
     
    42654599          ENDDO
    42664600
    4267           DO  k = nzb_max+1, nzt
     4601          DO  k = nzb_max_l+1, nzt
    42684602
    42694603             u_comp            = u(k,j,i) + u(k,j,i-1) - gu
     
    42954629       !$ACC PRESENT(drho_air, rho_air_zw, ddzw) &
    42964630       !$ACC PRESENT(tend) &
    4297        !$ACC PRESENT(hom(nzb+1:nzb_max,1,1:3,0)) &
     4631       !$ACC PRESENT(hom(nzb+1:nzb_max_l,1,1:3,0)) &
    42984632       !$ACC PRESENT(weight_substep(intermediate_timestep_count)) &
    42994633       !$ACC PRESENT(sums_us2_ws_l, sums_wsus_ws_l)
     
    43034637!--       The following loop computes the fluxes for the south boundary points
    43044638          j = nys
    4305           DO  k = nzb+1, nzb_max
     4639          DO  k = nzb+1, nzb_max_l
    43064640
    43074641             ibit14 = REAL( IBITS(advc_flags_1(k,j-1,i),14,1), KIND = wp )
     
    43424676          ENDDO
    43434677
    4344           DO  k = nzb_max+1, nzt
     4678          DO  k = nzb_max_l+1, nzt
    43454679
    43464680             v_comp                 = v(k,j,i) + v(k,j,i-1) - gv
     
    43644698             diss_d    = 0.0_wp
    43654699
    4366              DO  k = nzb+1, nzb_max
     4700             DO  k = nzb+1, nzb_max_l
    43674701
    43684702                ibit11 = REAL( IBITS(advc_flags_1(k,j,i),11,1), KIND = wp )
     
    46364970             ENDDO
    46374971
    4638              DO  k = nzb_max+1, nzt
     4972             DO  k = nzb_max_l+1, nzt
    46394973
    46404974                u_comp = u(k,j,i+1) + u(k,j,i)
     
    47875121          ENDDO
    47885122       ENDDO
    4789        sums_us2_ws_l(nzb,tn) = sums_us2_ws_l(nzb+1,tn)
    4790 
    47915123
    47925124    END SUBROUTINE advec_u_ws
     
    48045136
    48055137       USE control_parameters,                                                 &
    4806            ONLY:  intermediate_timestep_count, u_gtrans, v_gtrans
     5138           ONLY:  bc_dirichlet_l, bc_dirichlet_n, bc_dirichlet_r,              &
     5139                  bc_dirichlet_s, bc_radiation_l, bc_radiation_n,              &
     5140                  bc_radiation_r, bc_radiation_s, intermediate_timestep_count, &
     5141                  u_gtrans, v_gtrans
    48075142
    48085143       USE grid_variables,                                                     &
     
    48205155
    48215156
    4822        INTEGER(iwp) ::  i      !< grid index along x-direction
    4823        INTEGER(iwp) ::  j      !< grid index along y-direction
    4824        INTEGER(iwp) ::  k      !< grid index along z-direction
    4825        INTEGER(iwp) ::  k_mm   !< k-2 index in disretization, can be modified to avoid segmentation faults
    4826        INTEGER(iwp) ::  k_pp   !< k+2 index in disretization, can be modified to avoid segmentation faults
    4827        INTEGER(iwp) ::  k_ppp  !< k+3 index in disretization, can be modified to avoid segmentation faults
    4828        INTEGER(iwp) ::  tn = 0 !< number of OpenMP thread
     5157       INTEGER(iwp) ::  i         !< grid index along x-direction
     5158       INTEGER(iwp) ::  j         !< grid index along y-direction
     5159       INTEGER(iwp) ::  k         !< grid index along z-direction
     5160       INTEGER(iwp) ::  k_mm      !< k-2 index in disretization, can be modified to avoid segmentation faults
     5161       INTEGER(iwp) ::  k_pp      !< k+2 index in disretization, can be modified to avoid segmentation faults
     5162       INTEGER(iwp) ::  k_ppp     !< k+3 index in disretization, can be modified to avoid segmentation faults
     5163       INTEGER(iwp) ::  nzb_max_l !< index indicating upper bound for order degradation of horizontal advection terms
     5164       INTEGER(iwp) ::  tn = 0    !< number of OpenMP thread
    48295165       
    48305166       REAL(wp)    ::  ibit18 !< flag indicating 1st-order scheme along x-direction
     
    48825218       REAL(wp)    ::  v_comp_s !<
    48835219#endif
    4884 
     5220!
     5221!--    Set local version of nzb_max. At non-cyclic boundaries the order of the
     5222!--    advection need to be degraded near the boundary. Please note, in contrast
     5223!--    to the cache-optimized routines, nzb_max_l is set constantly for the
     5224!--    entire subdomain, in order to avoid unsymmetric loops which might be
     5225!--    an issue for GPUs.
     5226       IF( bc_dirichlet_l  .OR.  bc_radiation_l  .OR.                          &
     5227           bc_dirichlet_r  .OR.  bc_radiation_r  .OR.                          &
     5228           bc_dirichlet_s  .OR.  bc_radiation_s  .OR.                          &
     5229           bc_dirichlet_n  .OR.  bc_radiation_n )  THEN
     5230          nzb_max_l = nzt
     5231       ELSE
     5232          nzb_max_l = nzb_max
     5233       END IF
     5234       
    48855235       gu = 2.0_wp * u_gtrans
    48865236       gv = 2.0_wp * v_gtrans
     
    48915241       i = nxl
    48925242       DO  j = nysv, nyn
    4893           DO  k = nzb+1, nzb_max
     5243          DO  k = nzb+1, nzb_max_l
    48945244
    48955245             ibit20 = REAL( IBITS(advc_flags_1(k,j,i-1),20,1), KIND = wp )
     
    49305280          ENDDO
    49315281
    4932           DO  k = nzb_max+1, nzt
     5282          DO  k = nzb_max_l+1, nzt
    49335283
    49345284             u_comp                   = u(k,j-1,i) + u(k,j,i) - gu
     
    49615311       !$ACC PRESENT(drho_air, rho_air_zw, ddzw) &
    49625312       !$ACC PRESENT(tend) &
    4963        !$ACC PRESENT(hom(nzb+1:nzb_max,1,2:3,0)) &
     5313       !$ACC PRESENT(hom(nzb+1:nzb_max_l,1,2:3,0)) &
    49645314       !$ACC PRESENT(weight_substep(intermediate_timestep_count)) &
    49655315       !$ACC PRESENT(sums_vs2_ws_l, sums_wsvs_ws_l)
     
    49685318#ifndef _OPENACC
    49695319          j = nysv
    4970           DO  k = nzb+1, nzb_max
     5320          DO  k = nzb+1, nzb_max_l
    49715321
    49725322             ibit23 = REAL( IBITS(advc_flags_1(k,j-1,i),23,1), KIND = wp )
     
    50075357          ENDDO
    50085358
    5009           DO  k = nzb_max+1, nzt
     5359          DO  k = nzb_max_l+1, nzt
    50105360
    50115361             v_comp                 = v(k,j,i) + v(k,j-1,i) - gv
     
    50275377             diss_d    = 0.0_wp
    50285378
    5029              DO  k = nzb+1, nzb_max
     5379             DO  k = nzb+1, nzb_max_l
    50305380
    50315381                ibit20 = REAL( IBITS(advc_flags_1(k,j,i),20,1), KIND = wp )
     
    53025652             ENDDO
    53035653
    5304              DO  k = nzb_max+1, nzt
     5654             DO  k = nzb_max_l+1, nzt
    53055655
    53065656                u_comp = u(k,j-1,i+1) + u(k,j,i+1) - gu
     
    54575807          ENDDO
    54585808       ENDDO
    5459 !$ACC UPDATE HOST(sums_vs2_ws_l(nzb+1,tn))
    5460        sums_vs2_ws_l(nzb,tn) = sums_vs2_ws_l(nzb+1,tn)
    5461 !$ACC UPDATE DEVICE(sums_vs2_ws_l(nzb,tn))
    5462 
    54635809
    54645810    END SUBROUTINE advec_v_ws
     
    54765822
    54775823       USE control_parameters,                                                 &
    5478            ONLY:  intermediate_timestep_count, u_gtrans, v_gtrans
     5824           ONLY:  bc_dirichlet_l, bc_dirichlet_n, bc_dirichlet_r,              &
     5825                  bc_dirichlet_s, bc_radiation_l, bc_radiation_n,              &
     5826                  bc_radiation_r, bc_radiation_s, intermediate_timestep_count, &
     5827                  u_gtrans, v_gtrans
    54795828
    54805829       USE grid_variables,                                                     &
     
    54925841       IMPLICIT NONE
    54935842
    5494        INTEGER(iwp) ::  i      !< grid index along x-direction
    5495        INTEGER(iwp) ::  j      !< grid index along y-direction
    5496        INTEGER(iwp) ::  k      !< grid index along z-direction
    5497        INTEGER(iwp) ::  k_mm   !< k-2 index in disretization, can be modified to avoid segmentation faults
    5498        INTEGER(iwp) ::  k_pp   !< k+2 index in disretization, can be modified to avoid segmentation faults
    5499        INTEGER(iwp) ::  k_ppp  !< k+3 index in disretization, can be modified to avoid segmentation faults
    5500        INTEGER(iwp) ::  tn = 0 !< number of OpenMP thread
     5843       INTEGER(iwp) ::  i         !< grid index along x-direction
     5844       INTEGER(iwp) ::  j         !< grid index along y-direction
     5845       INTEGER(iwp) ::  k         !< grid index along z-direction
     5846       INTEGER(iwp) ::  k_mm      !< k-2 index in disretization, can be modified to avoid segmentation faults
     5847       INTEGER(iwp) ::  k_pp      !< k+2 index in disretization, can be modified to avoid segmentation faults
     5848       INTEGER(iwp) ::  k_ppp     !< k+3 index in disretization, can be modified to avoid segmentation faults
     5849       INTEGER(iwp) ::  nzb_max_l !< index indicating upper bound for order degradation of horizontal advection terms 
     5850       INTEGER(iwp) ::  tn = 0    !< number of OpenMP thread
    55015851       
    55025852       REAL(wp)    ::  ibit27 !< flag indicating 1st-order scheme along x-direction
     
    55555905       REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) ::  swap_flux_x_local_w !< discretized 6th-order flux at leftward-side of the grid box
    55565906#endif
    5557  
     5907!
     5908!--    Set local version of nzb_max. At non-cyclic boundaries the order of the
     5909!--    advection need to be degraded near the boundary. Please note, in contrast
     5910!--    to the cache-optimized routines, nzb_max_l is set constantly for the
     5911!--    entire subdomain, in order to avoid unsymmetric loops which might be
     5912!--    an issue for GPUs.
     5913       IF( bc_dirichlet_l  .OR.  bc_radiation_l  .OR.                          &
     5914           bc_dirichlet_r  .OR.  bc_radiation_r  .OR.                          &
     5915           bc_dirichlet_s  .OR.  bc_radiation_s  .OR.                          &
     5916           bc_dirichlet_n  .OR.  bc_radiation_n )  THEN
     5917          nzb_max_l = nzt
     5918       ELSE
     5919          nzb_max_l = nzb_max
     5920       END IF
     5921       
    55585922       gu = 2.0_wp * u_gtrans
    55595923       gv = 2.0_wp * v_gtrans
     
    55645928       i = nxl
    55655929       DO  j = nys, nyn
    5566           DO  k = nzb+1, nzb_max
     5930          DO  k = nzb+1, nzb_max_l
    55675931
    55685932             ibit29 = REAL( IBITS(advc_flags_1(k,j,i-1),29,1), KIND = wp )
     
    56035967          ENDDO
    56045968
    5605           DO  k = nzb_max+1, nzt
     5969          DO  k = nzb_max_l+1, nzt
    56065970
    56075971             u_comp                   = u(k+1,j,i) + u(k,j,i) - gu
     
    56345998       !$ACC PRESENT(rho_air, drho_air_zw, ddzu) &
    56355999       !$ACC PRESENT(tend) &
    5636        !$ACC PRESENT(hom(nzb+1:nzb_max,1,3,0)) &
     6000       !$ACC PRESENT(hom(nzb+1:nzb_max_l,1,3,0)) &
    56376001       !$ACC PRESENT(weight_substep(intermediate_timestep_count)) &
    5638        !$ACC PRESENT(sums_ws2_ws_l(nzb+1:nzb_max,0))
     6002       !$ACC PRESENT(sums_ws2_ws_l(nzb+1:nzb_max_l,0))
    56396003       DO i = nxl, nxr
    56406004
    56416005#ifndef _OPENACC
    56426006          j = nys
    5643           DO  k = nzb+1, nzb_max
     6007          DO  k = nzb+1, nzb_max_l
    56446008
    56456009             ibit32 = REAL( IBITS(advc_flags_2(k,j-1,i),0,1),  KIND = wp )
     
    56806044          ENDDO
    56816045
    5682           DO  k = nzb_max+1, nzt
     6046          DO  k = nzb_max_l+1, nzt
    56836047
    56846048             v_comp                 = v(k+1,j,i) + v(k,j,i) - gv
     
    57066070             diss_d = -ABS(w_comp) * ( w(k,j,i) - w(k-1,j,i) ) * adv_mom_1
    57076071
    5708              DO  k = nzb+1, nzb_max
     6072             DO  k = nzb+1, nzb_max_l
    57096073
    57106074                ibit29 = REAL( IBITS(advc_flags_1(k,j,i),29,1), KIND = wp )
     
    59666330             ENDDO
    59676331
    5968              DO  k = nzb_max+1, nzt
     6332             DO  k = nzb_max_l+1, nzt
    59696333
    59706334                u_comp = u(k+1,j,i+1) + u(k,j,i+1) - gu
Note: See TracChangeset for help on using the changeset viewer.