Changeset 1557


Ignore:
Timestamp:
Mar 5, 2015 4:43:04 PM (7 years ago)
Author:
suehring
Message:

Enable monotone advection for scalars in combination with fifth-order scheme using monotonic limiter.

Location:
palm/trunk/SOURCE
Files:
7 edited

Legend:

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

    r1375 r1557  
    2020! Current revisions:
    2121! ------------------
    22 !
     22! Enable monotone advection for scalars using monotonic limiter
    2323!
    2424! Former revisions:
     
    194194       USE control_parameters,                                                 &
    195195           ONLY:  cloud_physics, humidity, icloud_scheme, loop_optimization,   &
    196                   passive_scalar, precipitation, ocean, ws_scheme_mom,         &
    197                   ws_scheme_sca
     196                  monotonic_adjustment, passive_scalar, precipitation, ocean,  &
     197                  ws_scheme_mom, ws_scheme_sca
    198198
    199199       USE indices,                                                            &
     
    386386
    387387       USE control_parameters,                                                 &
    388            ONLY:  intermediate_timestep_count, u_gtrans, v_gtrans
     388           ONLY:  intermediate_timestep_count, monotonic_adjustment, u_gtrans, &
     389                  v_gtrans
    389390
    390391       USE grid_variables,                                                     &
     
    421422       INTEGER(iwp) ::  k     !:
    422423       INTEGER(iwp) ::  k_mm  !:
     424       INTEGER(iwp) ::  k_mmm !:
    423425       INTEGER(iwp) ::  k_pp  !:
    424426       INTEGER(iwp) ::  k_ppp !:
     
    428430       REAL(wp)     ::  div    !:
    429431       REAL(wp)     ::  flux_d !:
     432       REAL(wp)     ::  fd_1   !:
     433       REAL(wp)     ::  fl_1   !:
     434       REAL(wp)     ::  fn_1   !:
     435       REAL(wp)     ::  fr_1   !:
     436       REAL(wp)     ::  fs_1   !:
     437       REAL(wp)     ::  ft_1   !:
     438       REAL(wp)     ::  phi_d  !:
     439       REAL(wp)     ::  phi_l  !:
     440       REAL(wp)     ::  phi_n  !:
     441       REAL(wp)     ::  phi_r  !:
     442       REAL(wp)     ::  phi_s  !:
     443       REAL(wp)     ::  phi_t  !:
     444       REAL(wp)     ::  rd     !:
     445       REAL(wp)     ::  rl     !:
     446       REAL(wp)     ::  rn     !:
     447       REAL(wp)     ::  rr     !:
     448       REAL(wp)     ::  rs     !:
     449       REAL(wp)     ::  rt     !:
    430450       REAL(wp)     ::  u_comp !:
    431451       REAL(wp)     ::  v_comp !:
     
    697717                                         )
    698718!
     719!--       Apply monotonic adjustment.
     720          IF ( monotonic_adjustment )  THEN
     721!
     722!--          At first, calculate first order fluxes.
     723             u_comp = u(k,j,i) - u_gtrans
     724             fl_1   =  ( u_comp   * ( sk(k,j,i) + sk(k,j,i-1) )               &
     725                   -ABS( u_comp ) * ( sk(k,j,i) - sk(k,j,i-1) )               &
     726                       ) * adv_sca_1
     727
     728             u_comp = u(k,j,i+1) - u_gtrans
     729             fr_1   =  ( u_comp   * ( sk(k,j,i+1) + sk(k,j,i) )               &
     730                   -ABS( u_comp ) * ( sk(k,j,i+1) - sk(k,j,i) )               &
     731                       ) * adv_sca_1
     732
     733             v_comp = v(k,j,i) - v_gtrans
     734             fs_1   =  ( v_comp   * ( sk(k,j,i) + sk(k,j-1,i) )               &
     735                   -ABS( v_comp ) * ( sk(k,j,i) - sk(k,j-1,i) )               &
     736                       ) * adv_sca_1
     737
     738             v_comp = v(k,j+1,i) - v_gtrans
     739             fn_1   =  ( v_comp   * ( sk(k,j+1,i) + sk(k,j,i) )               &
     740                   -ABS( v_comp ) * ( sk(k,j+1,i) - sk(k,j,i) )               &
     741                       ) * adv_sca_1
     742
     743             fd_1   =  ( w(k-1,j,i)   * ( sk(k,j,i) + sk(k-1,j,i) )           &
     744                   -ABS( w(k-1,j,i) ) * ( sk(k,j,i) - sk(k-1,j,i) )           &
     745                       ) * adv_sca_1
     746
     747             ft_1   =  ( w(k,j,i)   * ( sk(k+1,j,i) + sk(k,j,i) )             &
     748                   -ABS( w(k,j,i) ) * ( sk(k+1,j,i) - sk(k,j,i) )             &
     749                      ) * adv_sca_1
     750!
     751!--          Calculate ratio of upwind gradients. Note, Min/Max is just to
     752!--          avoid if statements.
     753             rl     = ( MAX( 0.0, u(k,j,i) - u_gtrans ) *                     &
     754                           ABS( ( sk(k,j,i-1) - sk(k,j,i-2)            ) /    &
     755                                ( sk(k,j,i)   - sk(k,j,i-1) + 1E-20_wp )      &
     756                              ) +                                             &
     757                        MIN( 0.0, u(k,j,i) - u_gtrans ) *                     &
     758                           ABS( ( sk(k,j,i)   - sk(k,j,i+1)            ) /    &
     759                                ( sk(k,j,i-1) - sk(k,j,i)   + 1E-20_wp )      &
     760                              )                                               &
     761                      ) / ABS( u(k,j,i) - u_gtrans + 1E-20_wp )
     762
     763             rr     = ( MAX( 0.0, u(k,j,i+1) - u_gtrans ) *                   &
     764                           ABS( ( sk(k,j,i)   - sk(k,j,i-1)            ) /    &
     765                                ( sk(k,j,i+1) - sk(k,j,i)   + 1E-20_wp )      &
     766                              ) +                                             &
     767                        MIN( 0.0, u(k,j,i+1) - u_gtrans ) *                   &
     768                           ABS( ( sk(k,j,i+1) - sk(k,j,i+2)            ) /    &
     769                                ( sk(k,j,i)   - sk(k,j,i+1) + 1E-20_wp )      &
     770                              )                                               &
     771                      ) / ABS( u(k,j,i+1) - u_gtrans + 1E-20_wp )
     772
     773             rs     = ( MAX( 0.0, v(k,j,i) - v_gtrans ) *                     &
     774                           ABS( ( sk(k,j-1,i) - sk(k,j-2,i)            ) /    &
     775                                ( sk(k,j,i)   - sk(k,j-1,i) + 1E-20_wp )      &
     776                              ) +                                             &
     777                        MIN( 0.0, v(k,j,i) - v_gtrans ) *                     &
     778                           ABS( ( sk(k,j,i)   - sk(k,j+1,i)            ) /    &
     779                                ( sk(k,j-1,i) - sk(k,j,i)   + 1E-20_wp )      &
     780                              )                                               &
     781                      ) / ABS( v(k,j,i) - v_gtrans + 1E-20_wp )
     782
     783             rn     = ( MAX( 0.0, v(k,j+1,i) - v_gtrans ) *                   &
     784                           ABS( ( sk(k,j,i)   - sk(k,j-1,i)            ) /    &
     785                                ( sk(k,j+1,i) - sk(k,j,i)   + 1E-20_wp )      &
     786                              ) +                                             &
     787                       MIN( 0.0, v(k,j+1,i) - v_gtrans ) *                    &
     788                           ABS( ( sk(k,j+1,i) - sk(k,j+2,i)            ) /    &
     789                                ( sk(k,j,i)   - sk(k,j+1,i) + 1E-20_wp )      &
     790                              )                                               &
     791                      ) / ABS( v(k,j+1,i) - v_gtrans + 1E-20_wp )     
     792!   
     793!--          Reuse k_mm and compute k_mmm for the vertical gradient ratios.
     794!--          Note, for vertical advection below the third grid point above
     795!--          surface ( or below the model top) rd and rt are set to 0, i.e.
     796!--          use of first order scheme is enforced.
     797             k_mmm  = k - 3 * ibit8
     798
     799             rd     = ( MAX( 0.0, w(k-1,j,i) ) *                              &
     800                           ABS( ( sk(k_mm,j,i) - sk(k_mmm,j,i)           ) /  &
     801                                ( sk(k-1,j,i)  - sk(k_mm,j,i) + 1E-20_wp )    &
     802                              ) +                                             &
     803                        MIN( 0.0, w(k-1,j,i) ) *                              &
     804                           ABS( ( sk(k-1,j,i) - sk(k,j,i)            ) /      &
     805                                ( sk(k_mm,j,i) - sk(k-1,j,i)   + 1E-20_wp )   &
     806                              )                                               &
     807                      ) * ibit8 / ABS( w(k-1,j,i) + 1E-20_wp )
     808 
     809             rt     = ( MAX( 0.0, w(k,j,i) ) *                                &
     810                           ABS( ( sk(k,j,i)   - sk(k-1,j,i)            ) /    &
     811                                ( sk(k+1,j,i) - sk(k,j,i)   + 1E-20_wp )      &
     812                              ) +                                             &
     813                        MIN( 0.0, w(k,j,i) ) *                                &
     814                           ABS( ( sk(k+1,j,i) - sk(k_pp,j,i)           ) /    &
     815                                ( sk(k,j,i)   - sk(k+1,j,i) + 1E-20_wp )      &
     816                              )                                               &
     817                      ) * ibit8 / ABS( w(k,j,i) + 1E-20_wp )
     818!
     819!--           Calculate empirical limiter function (van Albada2 limiter).
     820              phi_l = MIN( 1.0_wp, ( 2.0_wp * rl ) /                          &
     821                                         ( rl**2 + 1.0_wp ) )
     822              phi_r = MIN( 1.0_wp, ( 2.0_wp * rr ) /                          &
     823                                         ( rr**2 + 1.0_wp ) )
     824              phi_s = MIN( 1.0_wp, ( 2.0_wp * rs ) /                          &
     825                                         ( rs**2 + 1.0_wp ) )
     826              phi_n = MIN( 1.0_wp, ( 2.0_wp * rn ) /                          &
     827                                         ( rn**2 + 1.0_wp ) )
     828              phi_d = MIN( 1.0_wp, ( 2.0_wp * rd ) /                          &
     829                                         ( rd**2 + 1.0_wp ) )
     830              phi_t = MIN( 1.0_wp, ( 2.0_wp * rt ) /                          &
     831                                         ( rt**2 + 1.0_wp ) )
     832!
     833!--           Calculate the resulting monotone flux.
     834              swap_flux_x_local(k,j,tn) = fl_1 - phi_l *                      &
     835                                        ( fl_1 - swap_flux_x_local(k,j,tn)  )
     836              flux_r(k)                 = fr_1 - phi_r *                      &
     837                                        ( fr_1 - flux_r(k)                  )
     838              swap_flux_y_local(k,tn)   = fs_1 - phi_s *                      &
     839                                        ( fs_1 - swap_flux_y_local(k,tn)    )
     840              flux_n(k)                 = fn_1 - phi_n *                      &
     841                                        ( fn_1 - flux_n(k)                  )
     842              flux_d                    = fd_1  !- phi_d *                      &
     843!                                         ( fd_1 - flux_d                     )
     844              flux_t(k)                 = ft_1 !- phi_t *                      &
     845!                                         ( ft_1 - flux_t(k)                  )
     846!
     847!--          Moreover, modify dissipation flux according to the limiter.
     848             swap_diss_x_local(k,j,tn) = swap_diss_x_local(k,j,tn) * phi_l
     849             diss_r(k)                 = diss_r(k)                 * phi_r
     850             swap_diss_y_local(k,tn)   = swap_diss_y_local(k,tn)   * phi_s
     851             diss_n(k)                 = diss_n(k)                 * phi_n
     852             diss_d                    = 0.0 !diss_d                    * phi_d
     853             diss_t(k)                 = 0.0 !diss_t(k)                 * phi_t
     854
     855          ENDIF
     856!
    699857!--       Calculate the divergence of the velocity field. A respective
    700858!--       correction is needed to overcome numerical instabilities caused
     
    756914          k_pp  = k + 2 * ( 1 - ibit6  )
    757915          k_mm  = k - 2 * ibit8
    758 
    759916
    760917          flux_t(k) = w(k,j,i) * (                                            &
     
    786943                             ( sk(k_ppp,j,i) - sk(k_mm,j,i) )                 &
    787944                                         )
     945
     946
     947!
     948!--       Apply monotonic adjustment.
     949          IF ( monotonic_adjustment )  THEN
     950!
     951!--          At first, calculate first order fluxes.
     952             u_comp = u(k,j,i) - u_gtrans
     953             fl_1   =  ( u_comp   * ( sk(k,j,i) + sk(k,j,i-1) )               &
     954                   -ABS( u_comp ) * ( sk(k,j,i) - sk(k,j,i-1) )               &
     955                       ) * adv_sca_1
     956
     957             u_comp = u(k,j,i+1) - u_gtrans
     958             fr_1   =  ( u_comp   * ( sk(k,j,i+1) + sk(k,j,i) )               &
     959                   -ABS( u_comp ) * ( sk(k,j,i+1) - sk(k,j,i) )               &
     960                       ) * adv_sca_1
     961
     962             v_comp = v(k,j,i) - v_gtrans
     963             fs_1   =  ( v_comp   * ( sk(k,j,i) + sk(k,j-1,i) )               &
     964                   -ABS( v_comp ) * ( sk(k,j,i) - sk(k,j-1,i) )               &
     965                       ) * adv_sca_1
     966
     967             v_comp = v(k,j+1,i) - v_gtrans
     968             fn_1   =  ( v_comp   * ( sk(k,j+1,i) + sk(k,j,i) )               &
     969                   -ABS( v_comp ) * ( sk(k,j+1,i) - sk(k,j,i) )               &
     970                       ) * adv_sca_1
     971
     972             fd_1   =  ( w(k-1,j,i)   * ( sk(k,j,i) + sk(k-1,j,i) )           &
     973                   -ABS( w(k-1,j,i) ) * ( sk(k,j,i) - sk(k-1,j,i) )           &
     974                       ) * adv_sca_1
     975
     976             ft_1   =  ( w(k,j,i)   * ( sk(k+1,j,i) + sk(k,j,i) )             &
     977                   -ABS( w(k,j,i) ) * ( sk(k+1,j,i) - sk(k,j,i) )             &
     978                       ) * adv_sca_1
     979!
     980!--          Calculate ratio of upwind gradients. Note, Min/Max is just to
     981!--          avoid if statements.
     982             rl     = ( MAX( 0.0, u(k,j,i) - u_gtrans ) *                     &
     983                           ABS( ( sk(k,j,i-1) - sk(k,j,i-2)            ) /    &
     984                                ( sk(k,j,i)   - sk(k,j,i-1) + 1E-20_wp )      &
     985                              ) +                                             &
     986                        MIN( 0.0, u(k,j,i) - u_gtrans ) *                     &
     987                           ABS( ( sk(k,j,i)   - sk(k,j,i+1)            ) /    &
     988                                ( sk(k,j,i-1) - sk(k,j,i)   + 1E-20_wp )      &
     989                              )                                               &
     990                      ) / ABS( u(k,j,i) - u_gtrans + 1E-20_wp )
     991
     992             rr     = ( MAX( 0.0, u(k,j,i+1) - u_gtrans ) *                   &
     993                           ABS( ( sk(k,j,i)   - sk(k,j,i-1)            ) /    &
     994                                ( sk(k,j,i+1) - sk(k,j,i)   + 1E-20_wp )      &
     995                              ) +                                             &
     996                        MIN( 0.0, u(k,j,i+1) - u_gtrans ) *                   &
     997                           ABS( ( sk(k,j,i+1) - sk(k,j,i+2)            ) /    &
     998                                ( sk(k,j,i)   - sk(k,j,i+1) + 1E-20_wp )      &
     999                              )                                               &
     1000                      ) / ABS( u(k,j,i+1) - u_gtrans + 1E-20_wp )
     1001
     1002             rs     = ( MAX( 0.0, v(k,j,i) - v_gtrans ) *                     &
     1003                           ABS( ( sk(k,j-1,i) - sk(k,j-2,i)            ) /    &
     1004                                ( sk(k,j,i)   - sk(k,j-1,i) + 1E-20_wp )      &
     1005                              ) +                                             &
     1006                        MIN( 0.0, v(k,j,i) - v_gtrans ) *                     &
     1007                           ABS( ( sk(k,j,i)   - sk(k,j+1,i)            ) /    &
     1008                                ( sk(k,j-1,i) - sk(k,j,i)   + 1E-20_wp )      &
     1009                              )                                               &
     1010                      ) / ABS( v(k,j,i) - v_gtrans + 1E-20_wp )
     1011
     1012             rn     = ( MAX( 0.0, v(k,j+1,i) - v_gtrans ) *                   &
     1013                           ABS( ( sk(k,j,i)   - sk(k,j-1,i)            ) /    &
     1014                                ( sk(k,j+1,i) - sk(k,j,i)   + 1E-20_wp )      &
     1015                              ) +                                             &
     1016                       MIN( 0.0, v(k,j+1,i) - v_gtrans ) *                    &
     1017                           ABS( ( sk(k,j+1,i) - sk(k,j+2,i)            ) /    &
     1018                                ( sk(k,j,i)   - sk(k,j+1,i) + 1E-20_wp )      &
     1019                              )                                               &
     1020                      ) / ABS( v(k,j+1,i) - v_gtrans + 1E-20_wp )     
     1021!
     1022!--          Reuse k_mm and compute k_mmm for the vertical gradient ratios.
     1023!--          Note, for vertical advection below the third grid point above
     1024!--          surface ( or below the model top) rd and rt are set to 0, i.e.
     1025!--          use of first order scheme is enforced.
     1026             k_mmm  = k - 3 * ibit8
     1027
     1028             rd     = ( MAX( 0.0, w(k-1,j,i) ) *                              &
     1029                           ABS( ( sk(k_mm,j,i) - sk(k_mmm,j,i)           ) /  &
     1030                                ( sk(k-1,j,i)  - sk(k_mm,j,i) + 1E-20_wp )    &
     1031                              ) +                                             &
     1032                        MIN( 0.0, w(k-1,j,i) ) *                              &
     1033                           ABS( ( sk(k-1,j,i) - sk(k,j,i)            ) /      &
     1034                                ( sk(k_mm,j,i) - sk(k-1,j,i)   + 1E-20_wp )   &
     1035                              )                                               &
     1036                      ) * ibit8 / ABS( w(k-1,j,i) + 1E-20_wp )
     1037 
     1038             rt     = ( MAX( 0.0, w(k,j,i) ) *                                &
     1039                           ABS( ( sk(k,j,i)   - sk(k-1,j,i)            ) /    &
     1040                                ( sk(k+1,j,i) - sk(k,j,i)   + 1E-20_wp )      &
     1041                              ) +                                             &
     1042                        MIN( 0.0, w(k,j,i) ) *                                &
     1043                           ABS( ( sk(k+1,j,i) - sk(k_pp,j,i)           ) /    &
     1044                                ( sk(k,j,i)   - sk(k+1,j,i) + 1E-20_wp )      &
     1045                              )                                               &
     1046                      ) * ibit8 / ABS( w(k,j,i) + 1E-20_wp )
     1047!
     1048!--           Calculate empirical limiter function (van Albada2 limiter).
     1049              phi_l = MIN( 1.0_wp, ( 2.0_wp * rl ) /                          &
     1050                                         ( rl**2 + 1.0_wp ) )
     1051              phi_r = MIN( 1.0_wp, ( 2.0_wp * rr ) /                          &
     1052                                         ( rr**2 + 1.0_wp ) )
     1053              phi_s = MIN( 1.0_wp, ( 2.0_wp * rs ) /                          &
     1054                                         ( rs**2 + 1.0_wp ) )
     1055              phi_n = MIN( 1.0_wp, ( 2.0_wp * rn ) /                          &
     1056                                         ( rn**2 + 1.0_wp ) )
     1057              phi_d = MIN( 1.0_wp, ( 2.0_wp * rd ) /                          &
     1058                                         ( rd**2 + 1.0_wp ) )
     1059              phi_t = MIN( 1.0_wp, ( 2.0_wp * rt ) /                          &
     1060                                         ( rt**2 + 1.0_wp ) )
     1061!
     1062!--           Calculate the resulting monotone flux.
     1063              swap_flux_x_local(k,j,tn) = fl_1 - phi_l *                      &
     1064                                        ( fl_1 - swap_flux_x_local(k,j,tn)  )
     1065              flux_r(k)                 = fr_1 - phi_r *                      &
     1066                                        ( fr_1 - flux_r(k)                  )
     1067              swap_flux_y_local(k,tn)   = fs_1 - phi_s *                      &
     1068                                        ( fs_1 - swap_flux_y_local(k,tn)    )
     1069              flux_n(k)                 = fn_1 - phi_n *                      &
     1070                                        ( fn_1 - flux_n(k)                  )
     1071              flux_d                    = fd_1 - phi_d *                      &
     1072                                        ( fd_1 - flux_d                     )
     1073              flux_t(k)                 = ft_1 - phi_t *                      &
     1074                                        ( ft_1 - flux_t(k)                  )
     1075!
     1076!--          Moreover, modify dissipation flux according to the limiter.
     1077             swap_diss_x_local(k,j,tn) = swap_diss_x_local(k,j,tn) * phi_l
     1078             diss_r(k)                 = diss_r(k)                 * phi_r
     1079             swap_diss_y_local(k,tn)   = swap_diss_y_local(k,tn)   * phi_s
     1080             diss_n(k)                 = diss_n(k)                 * phi_n
     1081             diss_d                    = diss_d                    * phi_d
     1082             diss_t(k)                 = diss_t(k)                 * phi_t
     1083
     1084          ENDIF
    7881085!
    7891086!--       Calculate the divergence of the velocity field. A respective
     
    8031100                                      ) + sk(k,j,i) * div
    8041101
     1102
    8051103          swap_flux_y_local(k,tn)   = flux_n(k)
    8061104          swap_diss_y_local(k,tn)   = diss_n(k)
     
    8131111
    8141112!
    815 !--    Evaluation of statistics
    816        SELECT CASE ( sk_char )
    817 
    818           CASE ( 'pt' )
    819 
    820              DO  k = nzb, nzt
    821                 sums_wspts_ws_l(k,tn) = sums_wspts_ws_l(k,tn) +               &
    822                                        ( flux_t(k) + diss_t(k) )              &
    823                                  * weight_substep(intermediate_timestep_count)
     1113!--    Evaluation of statistics. Please note, in case of using monotone limiter
     1114!--    vertical fluxes will be not calculated by the advective fluxes.
     1115       IF ( .NOT. monotonic_adjustment )  THEN
     1116
     1117          SELECT CASE ( sk_char )
     1118
     1119             CASE ( 'pt' )
     1120
     1121                DO  k = nzb, nzt
     1122                   sums_wspts_ws_l(k,tn) = sums_wspts_ws_l(k,tn) +            &
     1123                                          ( flux_t(k) + diss_t(k) )           &
     1124                                  * weight_substep(intermediate_timestep_count)
    8241125             ENDDO
    8251126             
    826           CASE ( 'sa' )
    827 
    828              DO  k = nzb, nzt
    829                 sums_wssas_ws_l(k,tn) = sums_wssas_ws_l(k,tn) +               &
    830                                        ( flux_t(k) + diss_t(k) )              &
    831                                  * weight_substep(intermediate_timestep_count)
    832              ENDDO
     1127             CASE ( 'sa' )
     1128
     1129                DO  k = nzb, nzt
     1130                   sums_wssas_ws_l(k,tn) = sums_wssas_ws_l(k,tn) +            &
     1131                                          ( flux_t(k) + diss_t(k) )           &
     1132                                  * weight_substep(intermediate_timestep_count)
     1133                ENDDO
    8331134             
    834           CASE ( 'q' )
    835 
    836              DO  k = nzb, nzt
    837                 sums_wsqs_ws_l(k,tn)  = sums_wsqs_ws_l(k,tn) +                &
    838                                       ( flux_t(k) + diss_t(k) )               &
    839                                  * weight_substep(intermediate_timestep_count)
    840              ENDDO
    841 
    842           CASE ( 'qr' )
    843 
    844              DO  k = nzb, nzt
    845                 sums_wsqrs_ws_l(k,tn)  = sums_wsqrs_ws_l(k,tn) +              &
    846                                       ( flux_t(k) + diss_t(k) )               &
    847                                  * weight_substep(intermediate_timestep_count)
    848              ENDDO
    849 
    850           CASE ( 'nr' )
    851 
    852              DO  k = nzb, nzt
    853                 sums_wsnrs_ws_l(k,tn)  = sums_wsnrs_ws_l(k,tn) +              &
    854                                       ( flux_t(k) + diss_t(k) )               &
    855                                  * weight_substep(intermediate_timestep_count)
    856              ENDDO
    857 
    858          END SELECT
     1135             CASE ( 'q' )
     1136
     1137                DO  k = nzb, nzt
     1138                   sums_wsqs_ws_l(k,tn)  = sums_wsqs_ws_l(k,tn) +             &
     1139                                         ( flux_t(k) + diss_t(k) )            &
     1140                                  * weight_substep(intermediate_timestep_count)
     1141                ENDDO
     1142
     1143             CASE ( 'qr' )
     1144
     1145                DO  k = nzb, nzt
     1146                   sums_wsqrs_ws_l(k,tn)  = sums_wsqrs_ws_l(k,tn) +           &
     1147                                         ( flux_t(k) + diss_t(k) )            &
     1148                                  * weight_substep(intermediate_timestep_count)
     1149                ENDDO
     1150
     1151             CASE ( 'nr' )
     1152
     1153                DO  k = nzb, nzt
     1154                   sums_wsnrs_ws_l(k,tn)  = sums_wsnrs_ws_l(k,tn) +           &
     1155                                         ( flux_t(k) + diss_t(k) )            &
     1156                                  * weight_substep(intermediate_timestep_count)
     1157                ENDDO
     1158
     1159            END SELECT
     1160         ENDIF
    8591161
    8601162    END SUBROUTINE advec_s_ws_ij
     
    22282530
    22292531       USE control_parameters,                                                &
    2230            ONLY:  intermediate_timestep_count, u_gtrans, v_gtrans
     2532           ONLY:  intermediate_timestep_count, monotonic_adjustment, u_gtrans,&
     2533                  v_gtrans
    22312534
    22322535       USE grid_variables,                                                    &
     
    22472550       CHARACTER (LEN = *), INTENT(IN)    ::  sk_char !:
    22482551       
    2249        INTEGER(iwp) ::  i      !:
    2250        INTEGER(iwp) ::  ibit0  !:
    2251        INTEGER(iwp) ::  ibit1  !:
    2252        INTEGER(iwp) ::  ibit2  !:
    2253        INTEGER(iwp) ::  ibit3  !:
    2254        INTEGER(iwp) ::  ibit4  !:
    2255        INTEGER(iwp) ::  ibit5  !:
    2256        INTEGER(iwp) ::  ibit6  !:
    2257        INTEGER(iwp) ::  ibit7  !:
    2258        INTEGER(iwp) ::  ibit8  !:
    2259        INTEGER(iwp) ::  j      !:
    2260        INTEGER(iwp) ::  k      !:
    2261        INTEGER(iwp) ::  k_mm   !:
    2262        INTEGER(iwp) ::  k_pp   !:
    2263        INTEGER(iwp) ::  k_ppp  !:
    2264        INTEGER(iwp) ::  tn = 0 !:
    2265        
    2266 #if defined( __nopointer )
    2267        REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  sk !:
    2268 #else
    2269        REAL(wp), DIMENSION(:,:,:), POINTER ::  sk !:
    2270 #endif
    2271 
    2272        REAL(wp) ::  diss_d !:
    2273        REAL(wp) ::  div    !:
    2274        REAL(wp) ::  flux_d !:
    2275        REAL(wp) ::  u_comp !:
    2276        REAL(wp) ::  v_comp !:
    2277        
    2278        REAL(wp), DIMENSION(nzb:nzt)   ::  diss_n !:
    2279        REAL(wp), DIMENSION(nzb:nzt)   ::  diss_r !:
    2280        REAL(wp), DIMENSION(nzb:nzt)   ::  diss_t !:
    2281        REAL(wp), DIMENSION(nzb:nzt)   ::  flux_n !:
    2282        REAL(wp), DIMENSION(nzb:nzt)   ::  flux_r !:
    2283        REAL(wp), DIMENSION(nzb:nzt)   ::  flux_t !:
    2284        
    2285        REAL(wp), DIMENSION(nzb+1:nzt) ::  swap_diss_y_local !:
    2286        REAL(wp), DIMENSION(nzb+1:nzt) ::  swap_flux_y_local !:
    2287        
    2288        REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) ::  swap_diss_x_local !:
    2289        REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) ::  swap_flux_x_local !:
    2290        
    2291 
    2292 !
    2293 !--    Compute the fluxes for the whole left boundary of the processor domain.
    2294        i = nxl
    2295        DO  j = nys, nyn
    2296 
    2297           DO  k = nzb+1, nzb_max
    2298 
    2299              ibit2 = IBITS(wall_flags_0(k,j,i),2,1)
    2300              ibit1 = IBITS(wall_flags_0(k,j,i),1,1)
    2301              ibit0 = IBITS(wall_flags_0(k,j,i),0,1)
    2302 
    2303              u_comp                 = u(k,j,i) - u_gtrans
    2304              swap_flux_x_local(k,j) = u_comp * (                              &
    2305                                              ( 37.0_wp * ibit2 * adv_sca_5    &
    2306                                           +     7.0_wp * ibit1 * adv_sca_3    &
    2307                                           +              ibit0 * adv_sca_1    &
    2308                                              ) *                              &
    2309                                           ( sk(k,j,i)   + sk(k,j,i-1)    )    &
    2310                                       -      (  8.0_wp * ibit2 * adv_sca_5    &
    2311                                           +              ibit1 * adv_sca_3    &
    2312                                              ) *                              &
    2313                                           ( sk(k,j,i+1) + sk(k,j,i-2)    )    &
    2314                                       +      (           ibit2 * adv_sca_5    &
    2315                                              ) *                              &
    2316                                           ( sk(k,j,i+2) + sk(k,j,i-3)    )    &
    2317                                                )
    2318 
    2319               swap_diss_x_local(k,j) = -ABS( u_comp ) * (                     &
    2320                                              ( 10.0_wp * ibit2 * adv_sca_5    &
    2321                                           +     3.0_wp * ibit1 * adv_sca_3    &
    2322                                           +              ibit0 * adv_sca_1    &
    2323                                              ) *                              &
    2324                                           ( sk(k,j,i)   - sk(k,j,i-1) )       &
    2325                                       -      (  5.0_wp * ibit2 * adv_sca_5    &
    2326                                           +              ibit1 * adv_sca_3    &
    2327                                              ) *                              &
    2328                                          ( sk(k,j,i+1) - sk(k,j,i-2)  )       &
    2329                                       +      (           ibit2 * adv_sca_5    &
    2330                                              ) *                              &
    2331                                           ( sk(k,j,i+2) - sk(k,j,i-3) )       &
    2332                                                         )
    2333 
    2334           ENDDO
    2335 
    2336           DO  k = nzb_max+1, nzt
    2337 
    2338              u_comp                 = u(k,j,i) - u_gtrans
    2339              swap_flux_x_local(k,j) = u_comp * (                              &
    2340                                       37.0_wp * ( sk(k,j,i)   + sk(k,j,i-1) ) &
    2341                                     -  8.0_wp * ( sk(k,j,i+1) + sk(k,j,i-2) ) &
    2342                                     +           ( sk(k,j,i+2) + sk(k,j,i-3) ) &
    2343                                                ) * adv_sca_5
    2344 
    2345              swap_diss_x_local(k,j) = -ABS( u_comp ) * (                      &
    2346                                       10.0_wp * ( sk(k,j,i)   - sk(k,j,i-1) ) &
    2347                                     -  5.0_wp * ( sk(k,j,i+1) - sk(k,j,i-2) ) &
    2348                                     +           ( sk(k,j,i+2) - sk(k,j,i-3) ) &
    2349                                                        ) * adv_sca_5
    2350 
    2351           ENDDO
    2352 
    2353        ENDDO
    2354 
    2355        DO  i = nxl, nxr
    2356 
    2357           j = nys
    2358           DO  k = nzb+1, nzb_max
    2359 
    2360              ibit5 = IBITS(wall_flags_0(k,j,i),5,1)
    2361              ibit4 = IBITS(wall_flags_0(k,j,i),4,1)
    2362              ibit3 = IBITS(wall_flags_0(k,j,i),3,1)
    2363 
    2364              v_comp               = v(k,j,i) - v_gtrans
    2365              swap_flux_y_local(k) = v_comp * (                                &
    2366                                              ( 37.0_wp * ibit5 * adv_sca_5    &
    2367                                           +     7.0_wp * ibit4 * adv_sca_3    &
    2368                                           +              ibit3 * adv_sca_1    &
    2369                                              ) *                              &
    2370                                          ( sk(k,j,i)  + sk(k,j-1,i)     )     &
    2371                                        -     (  8.0_wp * ibit5 * adv_sca_5    &
    2372                                           +              ibit4 * adv_sca_3    &
    2373                                               ) *                             &
    2374                                          ( sk(k,j+1,i) + sk(k,j-2,i)    )     &
    2375                                       +      (           ibit5 * adv_sca_5    &
    2376                                              ) *                              &
    2377                                         ( sk(k,j+2,i) + sk(k,j-3,i)     )     &
    2378                                              )
    2379 
    2380              swap_diss_y_local(k) = -ABS( v_comp ) * (                        &
    2381                                              ( 10.0_wp * ibit5 * adv_sca_5    &
    2382                                           +     3.0_wp * ibit4 * adv_sca_3    &
    2383                                           +              ibit3 * adv_sca_1    &
    2384                                              ) *                              &
    2385                                           ( sk(k,j,i)   - sk(k,j-1,i)    )    &
    2386                                       -      (  5.0_wp * ibit5 * adv_sca_5    &
    2387                                           +              ibit4 * adv_sca_3    &
    2388                                              ) *                              &
    2389                                           ( sk(k,j+1,i) - sk(k,j-2,i)    )    &
    2390                                       +      (           ibit5 * adv_sca_5    &
    2391                                              ) *                              &
    2392                                           ( sk(k,j+2,i) - sk(k,j-3,i)    )    &
    2393                                                      )
    2394 
    2395           ENDDO
    2396 !
    2397 !--       Above to the top of the highest topography. No degradation necessary.
    2398           DO  k = nzb_max+1, nzt
    2399 
    2400              v_comp               = v(k,j,i) - v_gtrans
    2401              swap_flux_y_local(k) = v_comp * (                               &
    2402                                     37.0_wp * ( sk(k,j,i)   + sk(k,j-1,i) )  &
    2403                                   -  8.0_wp * ( sk(k,j+1,i) + sk(k,j-2,i) )  &
    2404                                   +           ( sk(k,j+2,i) + sk(k,j-3,i) )  &
    2405                                              ) * adv_sca_5
    2406               swap_diss_y_local(k) = -ABS( v_comp ) * (                      &
    2407                                     10.0_wp * ( sk(k,j,i)   - sk(k,j-1,i) )  &
    2408                                   -  5.0_wp * ( sk(k,j+1,i) - sk(k,j-2,i) )  &
    2409                                   +             sk(k,j+2,i) - sk(k,j-3,i)    &
    2410                                                       ) * adv_sca_5
    2411 
    2412           ENDDO
    2413 
    2414           DO  j = nys, nyn
    2415 
    2416              flux_t(0) = 0.0_wp
    2417              diss_t(0) = 0.0_wp
    2418              flux_d    = 0.0_wp
    2419              diss_d    = 0.0_wp
    2420 
    2421              DO  k = nzb+1, nzb_max
    2422 
    2423                 ibit2 = IBITS(wall_flags_0(k,j,i),2,1)
    2424                 ibit1 = IBITS(wall_flags_0(k,j,i),1,1)
    2425                 ibit0 = IBITS(wall_flags_0(k,j,i),0,1)
    2426 
    2427                 u_comp    = u(k,j,i+1) - u_gtrans
    2428                 flux_r(k) = u_comp * (                                        &
    2429                           ( 37.0_wp * ibit2 * adv_sca_5                       &
    2430                       +      7.0_wp * ibit1 * adv_sca_3                       &
    2431                       +               ibit0 * adv_sca_1                       &
    2432                           ) *                                                 &
    2433                              ( sk(k,j,i+1) + sk(k,j,i)   )                    &
    2434                    -      (  8.0_wp * ibit2 * adv_sca_5                       &
    2435                        +              ibit1 * adv_sca_3                       &
    2436                           ) *                                                 &
    2437                              ( sk(k,j,i+2) + sk(k,j,i-1) )                    &
    2438                    +      (           ibit2 * adv_sca_5                       &
    2439                           ) *                                                 &
    2440                              ( sk(k,j,i+3) + sk(k,j,i-2) )                    &
    2441                                      )
    2442 
    2443                 diss_r(k) = -ABS( u_comp ) * (                                &
    2444                           ( 10.0_wp * ibit2 * adv_sca_5                       &
    2445                        +     3.0_wp * ibit1 * adv_sca_3                       &
    2446                        +              ibit0 * adv_sca_1                       &
    2447                           ) *                                                 &
    2448                              ( sk(k,j,i+1) - sk(k,j,i)   )                    &
    2449                    -      (  5.0_wp * ibit2 * adv_sca_5                       &
    2450                        +              ibit1 * adv_sca_3                       &
    2451                           ) *                                                 &
    2452                              ( sk(k,j,i+2) - sk(k,j,i-1) )                    &
    2453                    +      (           ibit2 * adv_sca_5                       &
    2454                           ) *                                                 &
    2455                              ( sk(k,j,i+3) - sk(k,j,i-2) )                    &
    2456                                              )
    2457 
    2458                 ibit5 = IBITS(wall_flags_0(k,j,i),5,1)
    2459                 ibit4 = IBITS(wall_flags_0(k,j,i),4,1)
    2460                 ibit3 = IBITS(wall_flags_0(k,j,i),3,1)
    2461 
    2462                 v_comp    = v(k,j+1,i) - v_gtrans
    2463                 flux_n(k) = v_comp * (                                        &
    2464                           ( 37.0_wp * ibit5 * adv_sca_5                       &
    2465                        +     7.0_wp * ibit4 * adv_sca_3                       &
    2466                        +              ibit3 * adv_sca_1                       &
    2467                           ) *                                                 &
    2468                              ( sk(k,j+1,i) + sk(k,j,i)   )                    &
    2469                    -      (  8.0_wp * ibit5 * adv_sca_5                       &
    2470                        +              ibit4 * adv_sca_3                       &
    2471                           ) *                                                 &
    2472                              ( sk(k,j+2,i) + sk(k,j-1,i) )                    &
    2473                    +      (           ibit5 * adv_sca_5                       &
    2474                           ) *                                                 &
    2475                              ( sk(k,j+3,i) + sk(k,j-2,i) )                    &
    2476                                      )
    2477 
    2478                 diss_n(k) = -ABS( v_comp ) * (                                &
    2479                           ( 10.0_wp * ibit5 * adv_sca_5                       &
    2480                        +     3.0_wp * ibit4 * adv_sca_3                       &
    2481                        +              ibit3 * adv_sca_1                       &
    2482                           ) *                                                 &
    2483                              ( sk(k,j+1,i) - sk(k,j,i)    )                   &
    2484                    -      (  5.0_wp * ibit5 * adv_sca_5                       &
    2485                        +              ibit4 * adv_sca_3                       &
    2486                           ) *                                                 &
    2487                              ( sk(k,j+2,i) - sk(k,j-1,i)  )                   &
    2488                    +      (           ibit5 * adv_sca_5                       &
    2489                           ) *                                                 &
    2490                              ( sk(k,j+3,i) - sk(k,j-2,i) )                    &
    2491                                              )
    2492 !
    2493 !--             k index has to be modified near bottom and top, else array
    2494 !--             subscripts will be exceeded.
    2495                 ibit8 = IBITS(wall_flags_0(k,j,i),8,1)
    2496                 ibit7 = IBITS(wall_flags_0(k,j,i),7,1)
    2497                 ibit6 = IBITS(wall_flags_0(k,j,i),6,1)
    2498 
    2499                 k_ppp = k + 3 * ibit8
    2500                 k_pp  = k + 2 * ( 1 - ibit6  )
    2501                 k_mm  = k - 2 * ibit8
    2502 
    2503 
    2504                 flux_t(k) = w(k,j,i) * (                                      &
    2505                            ( 37.0_wp * ibit8 * adv_sca_5                      &
    2506                         +     7.0_wp * ibit7 * adv_sca_3                      &
    2507                         +           ibit6 * adv_sca_1                         &
    2508                            ) *                                                &
    2509                                    ( sk(k+1,j,i)  + sk(k,j,i)    )            &
    2510                     -      (  8.0_wp * ibit8 * adv_sca_5                      &
    2511                         +              ibit7 * adv_sca_3                      &
    2512                            ) *                                                &
    2513                                    ( sk(k_pp,j,i) + sk(k-1,j,i)  )            &
    2514                     +      (           ibit8 * adv_sca_5                      &
    2515                            ) *     ( sk(k_ppp,j,i)+ sk(k_mm,j,i) )            &
    2516                                        )
    2517 
    2518                 diss_t(k) = -ABS( w(k,j,i) ) * (                              &
    2519                            ( 10.0_wp * ibit8 * adv_sca_5                      &
    2520                         +     3.0_wp * ibit7 * adv_sca_3                      &
    2521                         +              ibit6 * adv_sca_1                      &
    2522                            ) *                                                &
    2523                                    ( sk(k+1,j,i)   - sk(k,j,i)    )           &
    2524                     -      (  5.0_wp * ibit8 * adv_sca_5                      &
    2525                         +              ibit7 * adv_sca_3                      &
    2526                            ) *                                                &
    2527                                    ( sk(k_pp,j,i)  - sk(k-1,j,i)  )           &
    2528                     +      (           ibit8 * adv_sca_5                      &
    2529                            ) *                                                &
    2530                                    ( sk(k_ppp,j,i) - sk(k_mm,j,i) )           &
    2531                                                )
    2532 !
    2533 !--             Calculate the divergence of the velocity field. A respective
    2534 !--             correction is needed to overcome numerical instabilities caused
    2535 !--             by a not sufficient reduction of divergences near topography.
    2536                 div         =   ( u(k,j,i+1) - u(k,j,i)   ) * ddx             &
    2537                               + ( v(k,j+1,i) - v(k,j,i)   ) * ddy             &
    2538                               + ( w(k,j,i)   - w(k-1,j,i) ) * ddzw(k)
    2539 
    2540                 tend(k,j,i) = tend(k,j,i) - (                                 &
    2541                         ( flux_r(k) + diss_r(k) - swap_flux_x_local(k,j) -    &
    2542                           swap_diss_x_local(k,j)            ) * ddx           &
    2543                       + ( flux_n(k) + diss_n(k) - swap_flux_y_local(k)   -    &
    2544                           swap_diss_y_local(k)              ) * ddy           &
    2545                       + ( flux_t(k) + diss_t(k) - flux_d - diss_d             &
    2546                                                                ) * ddzw(k)    &
    2547                                             ) + sk(k,j,i) * div
    2548 
    2549                 swap_flux_y_local(k)   = flux_n(k)
    2550                 swap_diss_y_local(k)   = diss_n(k)
    2551                 swap_flux_x_local(k,j) = flux_r(k)
    2552                 swap_diss_x_local(k,j) = diss_r(k)
    2553                 flux_d                 = flux_t(k)
    2554                 diss_d                 = diss_t(k)
    2555 
    2556              ENDDO
    2557 
    2558              DO  k = nzb_max+1, nzt
    2559 
    2560                 u_comp    = u(k,j,i+1) - u_gtrans
    2561                 flux_r(k) = u_comp * (                                        &
    2562                       37.0_wp * ( sk(k,j,i+1) + sk(k,j,i)   )                 &
    2563                     -  8.0_wp * ( sk(k,j,i+2) + sk(k,j,i-1) )                 &
    2564                     +           ( sk(k,j,i+3) + sk(k,j,i-2) ) ) * adv_sca_5
    2565                 diss_r(k) = -ABS( u_comp ) * (                                &
    2566                       10.0_wp * ( sk(k,j,i+1) - sk(k,j,i)   )                 &
    2567                     -  5.0_wp * ( sk(k,j,i+2) - sk(k,j,i-1) )                 &
    2568                     +           ( sk(k,j,i+3) - sk(k,j,i-2) ) ) * adv_sca_5
    2569 
    2570                 v_comp    = v(k,j+1,i) - v_gtrans
    2571                 flux_n(k) = v_comp * (                                        &
    2572                       37.0_wp * ( sk(k,j+1,i) + sk(k,j,i)   )                 &
    2573                     -  8.0_wp * ( sk(k,j+2,i) + sk(k,j-1,i) )                 &
    2574                     +           ( sk(k,j+3,i) + sk(k,j-2,i) ) ) * adv_sca_5
    2575                 diss_n(k) = -ABS( v_comp ) * (                                &
    2576                       10.0_wp * ( sk(k,j+1,i) - sk(k,j,i)   )                 &
    2577                     -  5.0_wp * ( sk(k,j+2,i) - sk(k,j-1,i) )                 &
    2578                     +           ( sk(k,j+3,i) - sk(k,j-2,i) ) ) * adv_sca_5
    2579 !
    2580 !--             k index has to be modified near bottom and top, else array
    2581 !--             subscripts will be exceeded.
    2582                 ibit8 = IBITS(wall_flags_0(k,j,i),8,1)
    2583                 ibit7 = IBITS(wall_flags_0(k,j,i),7,1)
    2584                 ibit6 = IBITS(wall_flags_0(k,j,i),6,1)
    2585 
    2586                 k_ppp = k + 3 * ibit8
    2587                 k_pp  = k + 2 * ( 1 - ibit6  )
    2588                 k_mm  = k - 2 * ibit8
    2589 
    2590 
    2591                 flux_t(k) = w(k,j,i) * (                                      &
    2592                            ( 37.0_wp * ibit8 * adv_sca_5                      &
    2593                         +     7.0_wp * ibit7 * adv_sca_3                      &
    2594                         +              ibit6 * adv_sca_1                      &
    2595                            ) *                                                &
    2596                                    ( sk(k+1,j,i)  + sk(k,j,i)     )           &
    2597                     -      (  8.0_wp * ibit8 * adv_sca_5                      &
    2598                         +              ibit7 * adv_sca_3                      &
    2599                            ) *                                                &
    2600                                    ( sk(k_pp,j,i) + sk(k-1,j,i)   )           &
    2601                     +      (           ibit8 * adv_sca_5                      &
    2602                            ) *     ( sk(k_ppp,j,i)+ sk(k_mm,j,i)  )           &
    2603                                        )
    2604 
    2605                 diss_t(k) = -ABS( w(k,j,i) ) * (                              &
    2606                            ( 10.0_wp * ibit8 * adv_sca_5                      &
    2607                         +     3.0_wp * ibit7 * adv_sca_3                      &
    2608                         +              ibit6 * adv_sca_1                      &
    2609                            ) *                                                &
    2610                                    ( sk(k+1,j,i)   - sk(k,j,i)    )           &
    2611                     -      (  5.0_wp * ibit8 * adv_sca_5                      &
    2612                         +              ibit7 * adv_sca_3                      &
    2613                            ) *                                                &
    2614                                    ( sk(k_pp,j,i)  - sk(k-1,j,i)  )           &
    2615                     +      (           ibit8 * adv_sca_5                      &
    2616                            ) *                                                &
    2617                                    ( sk(k_ppp,j,i) - sk(k_mm,j,i) )           &
    2618                                                )
    2619 !
    2620 !--             Calculate the divergence of the velocity field. A respective
    2621 !--             correction is needed to overcome numerical instabilities introduced
    2622 !--             by a not sufficient reduction of divergences near topography.
    2623                 div         =   ( u(k,j,i+1) - u(k,j,i)   ) * ddx             &
    2624                               + ( v(k,j+1,i) - v(k,j,i)   ) * ddy             &
    2625                               + ( w(k,j,i)   - w(k-1,j,i) ) * ddzw(k)
    2626 
    2627                 tend(k,j,i) = tend(k,j,i) - (                                 &
    2628                         ( flux_r(k) + diss_r(k) - swap_flux_x_local(k,j) -    &
    2629                           swap_diss_x_local(k,j)            ) * ddx           &
    2630                       + ( flux_n(k) + diss_n(k) - swap_flux_y_local(k)   -    &
    2631                           swap_diss_y_local(k)              ) * ddy           &
    2632                       + ( flux_t(k) + diss_t(k) - flux_d - diss_d             &
    2633                                                                ) * ddzw(k)    &
    2634                                             ) + sk(k,j,i) * div
    2635 
    2636                 swap_flux_y_local(k)   = flux_n(k)
    2637                 swap_diss_y_local(k)   = diss_n(k)
    2638                 swap_flux_x_local(k,j) = flux_r(k)
    2639                 swap_diss_x_local(k,j) = diss_r(k)
    2640                 flux_d                 = flux_t(k)
    2641                 diss_d                 = diss_t(k)
    2642 
    2643              ENDDO
    2644 !
    2645 !--          evaluation of statistics
    2646              SELECT CASE ( sk_char )
    2647 
    2648                  CASE ( 'pt' )
    2649                     DO  k = nzb, nzt
    2650                        sums_wspts_ws_l(k,tn) = sums_wspts_ws_l(k,tn)         &
    2651                         + ( flux_t(k) + diss_t(k) )                          &
    2652                         *   weight_substep(intermediate_timestep_count)
    2653                     ENDDO
    2654                  CASE ( 'sa' )
    2655                     DO  k = nzb, nzt
    2656                        sums_wssas_ws_l(k,tn) = sums_wssas_ws_l(k,tn)         &
    2657                         + ( flux_t(k) + diss_t(k) )                          &
    2658                         *   weight_substep(intermediate_timestep_count)
    2659                     ENDDO
    2660                  CASE ( 'q' )
    2661                     DO  k = nzb, nzt
    2662                        sums_wsqs_ws_l(k,tn) = sums_wsqs_ws_l(k,tn)           &
    2663                         + ( flux_t(k) + diss_t(k) )                          &
    2664                         *   weight_substep(intermediate_timestep_count)
    2665                     ENDDO
    2666                  CASE ( 'qr' )
    2667                     DO  k = nzb, nzt
    2668                        sums_wsqrs_ws_l(k,tn) = sums_wsqrs_ws_l(k,tn)         &
    2669                         + ( flux_t(k) + diss_t(k) )                          &
    2670                         *   weight_substep(intermediate_timestep_count)
    2671                     ENDDO
    2672                  CASE ( 'nr' )
    2673                     DO  k = nzb, nzt
    2674                        sums_wsnrs_ws_l(k,tn) = sums_wsnrs_ws_l(k,tn)         &
    2675                         + ( flux_t(k) + diss_t(k) )                          &
    2676                         *   weight_substep(intermediate_timestep_count)
    2677                     ENDDO
    2678 
    2679               END SELECT
    2680 
    2681          ENDDO
    2682       ENDDO
    2683 
    2684     END SUBROUTINE advec_s_ws
    2685 
    2686 
    2687 !------------------------------------------------------------------------------!
    2688 ! Scalar advection - Call for all grid points - accelerator version
    2689 !------------------------------------------------------------------------------!
    2690     SUBROUTINE advec_s_ws_acc ( sk, sk_char )
    2691 
    2692        USE arrays_3d,                                                         &
    2693            ONLY:  ddzw, tend, u, v, w
    2694 
    2695        USE constants,                                                         &
    2696            ONLY:  adv_sca_1, adv_sca_3, adv_sca_5
    2697 
    2698        USE control_parameters,                                                &
    2699            ONLY:  intermediate_timestep_count, u_gtrans, v_gtrans
    2700 
    2701        USE grid_variables,                                                    &
    2702            ONLY:  ddx, ddy
    2703 
    2704        USE indices,                                                           &
    2705            ONLY:  i_left, i_right, j_north, j_south, nxlg, nxrg, nyng, nysg,  &
    2706                   nzb, nzb_max, nzt, wall_flags_0
    2707 
    2708        USE kinds
    2709        
    2710 !        USE statistics,                                                       &
    2711 !            ONLY:  sums_wspts_ws_l, sums_wsqs_ws_l, sums_wssas_ws_l,          &
    2712 !                   sums_wsqrs_ws_l, sums_wsnrs_ws_l, weight_substep
    2713 
    2714        IMPLICIT NONE
    2715 
    2716        CHARACTER (LEN = *), INTENT(IN)    :: sk_char !:
    2717 
    27182552       INTEGER(iwp) ::  i      !:
    27192553       INTEGER(iwp) ::  ibit0  !:
     
    27332567       INTEGER(iwp) ::  k_ppp  !:
    27342568       INTEGER(iwp) ::  tn = 0 !:
     2569       
     2570#if defined( __nopointer )
     2571       REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  sk !:
     2572#else
     2573       REAL(wp), DIMENSION(:,:,:), POINTER ::  sk !:
     2574#endif
     2575
     2576       REAL(wp) ::  diss_d !:
     2577       REAL(wp) ::  div    !:
     2578       REAL(wp) ::  flux_d !:
     2579       REAL(wp) ::  fd_1   !:
     2580       REAL(wp) ::  fl_1   !:
     2581       REAL(wp) ::  fn_1   !:
     2582       REAL(wp) ::  fr_1   !:
     2583       REAL(wp) ::  fs_1   !:
     2584       REAL(wp) ::  ft_1   !:
     2585       REAL(wp) ::  phi_d  !:
     2586       REAL(wp) ::  phi_l  !:
     2587       REAL(wp) ::  phi_n  !:
     2588       REAL(wp) ::  phi_r  !:
     2589       REAL(wp) ::  phi_s  !:
     2590       REAL(wp) ::  phi_t  !:
     2591       REAL(wp) ::  rd     !:
     2592       REAL(wp) ::  rl     !:
     2593       REAL(wp) ::  rn     !:
     2594       REAL(wp) ::  rr     !:
     2595       REAL(wp) ::  rs     !:
     2596       REAL(wp) ::  rt     !:
     2597       REAL(wp) ::  u_comp !:
     2598       REAL(wp) ::  v_comp !:
     2599       
     2600       REAL(wp), DIMENSION(nzb:nzt)   ::  diss_n !:
     2601       REAL(wp), DIMENSION(nzb:nzt)   ::  diss_r !:
     2602       REAL(wp), DIMENSION(nzb:nzt)   ::  diss_t !:
     2603       REAL(wp), DIMENSION(nzb:nzt)   ::  flux_n !:
     2604       REAL(wp), DIMENSION(nzb:nzt)   ::  flux_r !:
     2605       REAL(wp), DIMENSION(nzb:nzt)   ::  flux_t !:
     2606       
     2607       REAL(wp), DIMENSION(nzb+1:nzt) ::  swap_diss_y_local !:
     2608       REAL(wp), DIMENSION(nzb+1:nzt) ::  swap_flux_y_local !:
     2609       
     2610       REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) ::  swap_diss_x_local !:
     2611       REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) ::  swap_flux_x_local !:
     2612       
     2613
     2614!
     2615!--    Compute the fluxes for the whole left boundary of the processor domain.
     2616       i = nxl
     2617       DO  j = nys, nyn
     2618
     2619          DO  k = nzb+1, nzb_max
     2620
     2621             ibit2 = IBITS(wall_flags_0(k,j,i),2,1)
     2622             ibit1 = IBITS(wall_flags_0(k,j,i),1,1)
     2623             ibit0 = IBITS(wall_flags_0(k,j,i),0,1)
     2624
     2625             u_comp                 = u(k,j,i) - u_gtrans
     2626             swap_flux_x_local(k,j) = u_comp * (                              &
     2627                                             ( 37.0_wp * ibit2 * adv_sca_5    &
     2628                                          +     7.0_wp * ibit1 * adv_sca_3    &
     2629                                          +              ibit0 * adv_sca_1    &
     2630                                             ) *                              &
     2631                                          ( sk(k,j,i)   + sk(k,j,i-1)    )    &
     2632                                      -      (  8.0_wp * ibit2 * adv_sca_5    &
     2633                                          +              ibit1 * adv_sca_3    &
     2634                                             ) *                              &
     2635                                          ( sk(k,j,i+1) + sk(k,j,i-2)    )    &
     2636                                      +      (           ibit2 * adv_sca_5    &
     2637                                             ) *                              &
     2638                                          ( sk(k,j,i+2) + sk(k,j,i-3)    )    &
     2639                                               )
     2640
     2641              swap_diss_x_local(k,j) = -ABS( u_comp ) * (                     &
     2642                                             ( 10.0_wp * ibit2 * adv_sca_5    &
     2643                                          +     3.0_wp * ibit1 * adv_sca_3    &
     2644                                          +              ibit0 * adv_sca_1    &
     2645                                             ) *                              &
     2646                                          ( sk(k,j,i)   - sk(k,j,i-1) )       &
     2647                                      -      (  5.0_wp * ibit2 * adv_sca_5    &
     2648                                          +              ibit1 * adv_sca_3    &
     2649                                             ) *                              &
     2650                                         ( sk(k,j,i+1) - sk(k,j,i-2)  )       &
     2651                                      +      (           ibit2 * adv_sca_5    &
     2652                                             ) *                              &
     2653                                          ( sk(k,j,i+2) - sk(k,j,i-3) )       &
     2654                                                        )
     2655
     2656          ENDDO
     2657
     2658          DO  k = nzb_max+1, nzt
     2659
     2660             u_comp                 = u(k,j,i) - u_gtrans
     2661             swap_flux_x_local(k,j) = u_comp * (                              &
     2662                                      37.0_wp * ( sk(k,j,i)   + sk(k,j,i-1) ) &
     2663                                    -  8.0_wp * ( sk(k,j,i+1) + sk(k,j,i-2) ) &
     2664                                    +           ( sk(k,j,i+2) + sk(k,j,i-3) ) &
     2665                                               ) * adv_sca_5
     2666
     2667             swap_diss_x_local(k,j) = -ABS( u_comp ) * (                      &
     2668                                      10.0_wp * ( sk(k,j,i)   - sk(k,j,i-1) ) &
     2669                                    -  5.0_wp * ( sk(k,j,i+1) - sk(k,j,i-2) ) &
     2670                                    +           ( sk(k,j,i+2) - sk(k,j,i-3) ) &
     2671                                                       ) * adv_sca_5
     2672
     2673          ENDDO
     2674
     2675       ENDDO
     2676
     2677       DO  i = nxl, nxr
     2678
     2679          j = nys
     2680          DO  k = nzb+1, nzb_max
     2681
     2682             ibit5 = IBITS(wall_flags_0(k,j,i),5,1)
     2683             ibit4 = IBITS(wall_flags_0(k,j,i),4,1)
     2684             ibit3 = IBITS(wall_flags_0(k,j,i),3,1)
     2685
     2686             v_comp               = v(k,j,i) - v_gtrans
     2687             swap_flux_y_local(k) = v_comp * (                                &
     2688                                             ( 37.0_wp * ibit5 * adv_sca_5    &
     2689                                          +     7.0_wp * ibit4 * adv_sca_3    &
     2690                                          +              ibit3 * adv_sca_1    &
     2691                                             ) *                              &
     2692                                         ( sk(k,j,i)  + sk(k,j-1,i)     )     &
     2693                                       -     (  8.0_wp * ibit5 * adv_sca_5    &
     2694                                          +              ibit4 * adv_sca_3    &
     2695                                              ) *                             &
     2696                                         ( sk(k,j+1,i) + sk(k,j-2,i)    )     &
     2697                                      +      (           ibit5 * adv_sca_5    &
     2698                                             ) *                              &
     2699                                        ( sk(k,j+2,i) + sk(k,j-3,i)     )     &
     2700                                             )
     2701
     2702             swap_diss_y_local(k) = -ABS( v_comp ) * (                        &
     2703                                             ( 10.0_wp * ibit5 * adv_sca_5    &
     2704                                          +     3.0_wp * ibit4 * adv_sca_3    &
     2705                                          +              ibit3 * adv_sca_1    &
     2706                                             ) *                              &
     2707                                          ( sk(k,j,i)   - sk(k,j-1,i)    )    &
     2708                                      -      (  5.0_wp * ibit5 * adv_sca_5    &
     2709                                          +              ibit4 * adv_sca_3    &
     2710                                             ) *                              &
     2711                                          ( sk(k,j+1,i) - sk(k,j-2,i)    )    &
     2712                                      +      (           ibit5 * adv_sca_5    &
     2713                                             ) *                              &
     2714                                          ( sk(k,j+2,i) - sk(k,j-3,i)    )    &
     2715                                                     )
     2716
     2717          ENDDO
     2718!
     2719!--       Above to the top of the highest topography. No degradation necessary.
     2720          DO  k = nzb_max+1, nzt
     2721
     2722             v_comp               = v(k,j,i) - v_gtrans
     2723             swap_flux_y_local(k) = v_comp * (                               &
     2724                                    37.0_wp * ( sk(k,j,i)   + sk(k,j-1,i) )  &
     2725                                  -  8.0_wp * ( sk(k,j+1,i) + sk(k,j-2,i) )  &
     2726                                  +           ( sk(k,j+2,i) + sk(k,j-3,i) )  &
     2727                                             ) * adv_sca_5
     2728              swap_diss_y_local(k) = -ABS( v_comp ) * (                      &
     2729                                    10.0_wp * ( sk(k,j,i)   - sk(k,j-1,i) )  &
     2730                                  -  5.0_wp * ( sk(k,j+1,i) - sk(k,j-2,i) )  &
     2731                                  +             sk(k,j+2,i) - sk(k,j-3,i)    &
     2732                                                      ) * adv_sca_5
     2733
     2734          ENDDO
     2735
     2736          DO  j = nys, nyn
     2737
     2738             flux_t(0) = 0.0_wp
     2739             diss_t(0) = 0.0_wp
     2740             flux_d    = 0.0_wp
     2741             diss_d    = 0.0_wp
     2742
     2743             DO  k = nzb+1, nzb_max
     2744
     2745                ibit2 = IBITS(wall_flags_0(k,j,i),2,1)
     2746                ibit1 = IBITS(wall_flags_0(k,j,i),1,1)
     2747                ibit0 = IBITS(wall_flags_0(k,j,i),0,1)
     2748
     2749                u_comp    = u(k,j,i+1) - u_gtrans
     2750                flux_r(k) = u_comp * (                                        &
     2751                          ( 37.0_wp * ibit2 * adv_sca_5                       &
     2752                      +      7.0_wp * ibit1 * adv_sca_3                       &
     2753                      +               ibit0 * adv_sca_1                       &
     2754                          ) *                                                 &
     2755                             ( sk(k,j,i+1) + sk(k,j,i)   )                    &
     2756                   -      (  8.0_wp * ibit2 * adv_sca_5                       &
     2757                       +              ibit1 * adv_sca_3                       &
     2758                          ) *                                                 &
     2759                             ( sk(k,j,i+2) + sk(k,j,i-1) )                    &
     2760                   +      (           ibit2 * adv_sca_5                       &
     2761                          ) *                                                 &
     2762                             ( sk(k,j,i+3) + sk(k,j,i-2) )                    &
     2763                                     )
     2764
     2765                diss_r(k) = -ABS( u_comp ) * (                                &
     2766                          ( 10.0_wp * ibit2 * adv_sca_5                       &
     2767                       +     3.0_wp * ibit1 * adv_sca_3                       &
     2768                       +              ibit0 * adv_sca_1                       &
     2769                          ) *                                                 &
     2770                             ( sk(k,j,i+1) - sk(k,j,i)   )                    &
     2771                   -      (  5.0_wp * ibit2 * adv_sca_5                       &
     2772                       +              ibit1 * adv_sca_3                       &
     2773                          ) *                                                 &
     2774                             ( sk(k,j,i+2) - sk(k,j,i-1) )                    &
     2775                   +      (           ibit2 * adv_sca_5                       &
     2776                          ) *                                                 &
     2777                             ( sk(k,j,i+3) - sk(k,j,i-2) )                    &
     2778                                             )
     2779
     2780                ibit5 = IBITS(wall_flags_0(k,j,i),5,1)
     2781                ibit4 = IBITS(wall_flags_0(k,j,i),4,1)
     2782                ibit3 = IBITS(wall_flags_0(k,j,i),3,1)
     2783
     2784                v_comp    = v(k,j+1,i) - v_gtrans
     2785                flux_n(k) = v_comp * (                                        &
     2786                          ( 37.0_wp * ibit5 * adv_sca_5                       &
     2787                       +     7.0_wp * ibit4 * adv_sca_3                       &
     2788                       +              ibit3 * adv_sca_1                       &
     2789                          ) *                                                 &
     2790                             ( sk(k,j+1,i) + sk(k,j,i)   )                    &
     2791                   -      (  8.0_wp * ibit5 * adv_sca_5                       &
     2792                       +              ibit4 * adv_sca_3                       &
     2793                          ) *                                                 &
     2794                             ( sk(k,j+2,i) + sk(k,j-1,i) )                    &
     2795                   +      (           ibit5 * adv_sca_5                       &
     2796                          ) *                                                 &
     2797                             ( sk(k,j+3,i) + sk(k,j-2,i) )                    &
     2798                                     )
     2799
     2800                diss_n(k) = -ABS( v_comp ) * (                                &
     2801                          ( 10.0_wp * ibit5 * adv_sca_5                       &
     2802                       +     3.0_wp * ibit4 * adv_sca_3                       &
     2803                       +              ibit3 * adv_sca_1                       &
     2804                          ) *                                                 &
     2805                             ( sk(k,j+1,i) - sk(k,j,i)    )                   &
     2806                   -      (  5.0_wp * ibit5 * adv_sca_5                       &
     2807                       +              ibit4 * adv_sca_3                       &
     2808                          ) *                                                 &
     2809                             ( sk(k,j+2,i) - sk(k,j-1,i)  )                   &
     2810                   +      (           ibit5 * adv_sca_5                       &
     2811                          ) *                                                 &
     2812                             ( sk(k,j+3,i) - sk(k,j-2,i) )                    &
     2813                                             )
     2814!
     2815!--             k index has to be modified near bottom and top, else array
     2816!--             subscripts will be exceeded.
     2817                ibit8 = IBITS(wall_flags_0(k,j,i),8,1)
     2818                ibit7 = IBITS(wall_flags_0(k,j,i),7,1)
     2819                ibit6 = IBITS(wall_flags_0(k,j,i),6,1)
     2820
     2821                k_ppp = k + 3 * ibit8
     2822                k_pp  = k + 2 * ( 1 - ibit6  )
     2823                k_mm  = k - 2 * ibit8
     2824
     2825
     2826                flux_t(k) = w(k,j,i) * (                                      &
     2827                           ( 37.0_wp * ibit8 * adv_sca_5                      &
     2828                        +     7.0_wp * ibit7 * adv_sca_3                      &
     2829                        +           ibit6 * adv_sca_1                         &
     2830                           ) *                                                &
     2831                                   ( sk(k+1,j,i)  + sk(k,j,i)    )            &
     2832                    -      (  8.0_wp * ibit8 * adv_sca_5                      &
     2833                        +              ibit7 * adv_sca_3                      &
     2834                           ) *                                                &
     2835                                   ( sk(k_pp,j,i) + sk(k-1,j,i)  )            &
     2836                    +      (           ibit8 * adv_sca_5                      &
     2837                           ) *     ( sk(k_ppp,j,i)+ sk(k_mm,j,i) )            &
     2838                                       )
     2839
     2840                diss_t(k) = -ABS( w(k,j,i) ) * (                              &
     2841                           ( 10.0_wp * ibit8 * adv_sca_5                      &
     2842                        +     3.0_wp * ibit7 * adv_sca_3                      &
     2843                        +              ibit6 * adv_sca_1                      &
     2844                           ) *                                                &
     2845                                   ( sk(k+1,j,i)   - sk(k,j,i)    )           &
     2846                    -      (  5.0_wp * ibit8 * adv_sca_5                      &
     2847                        +              ibit7 * adv_sca_3                      &
     2848                           ) *                                                &
     2849                                   ( sk(k_pp,j,i)  - sk(k-1,j,i)  )           &
     2850                    +      (           ibit8 * adv_sca_5                      &
     2851                           ) *                                                &
     2852                                   ( sk(k_ppp,j,i) - sk(k_mm,j,i) )           &
     2853                                               )
     2854!
     2855!--             Apply monotonic adjustment.
     2856                IF ( monotonic_adjustment )  THEN
     2857!
     2858!--                At first, calculate first order fluxes.
     2859                   u_comp = u(k,j,i) - u_gtrans
     2860                   fl_1   =  ( u_comp   * ( sk(k,j,i) + sk(k,j,i-1) )         &
     2861                         -ABS( u_comp ) * ( sk(k,j,i) - sk(k,j,i-1) )         &
     2862                             ) * adv_sca_1
     2863
     2864                   u_comp = u(k,j,i+1) - u_gtrans
     2865                   fr_1   =  ( u_comp   * ( sk(k,j,i+1) + sk(k,j,i) )         &
     2866                         -ABS( u_comp ) * ( sk(k,j,i+1) - sk(k,j,i) )         &
     2867                             ) * adv_sca_1
     2868
     2869                   v_comp = v(k,j,i) - v_gtrans
     2870                   fs_1   =  ( v_comp   * ( sk(k,j,i) + sk(k,j-1,i) )         &
     2871                         -ABS( v_comp ) * ( sk(k,j,i) - sk(k,j-1,i) )         &
     2872                             ) * adv_sca_1
     2873
     2874                   v_comp = v(k,j+1,i) - v_gtrans
     2875                   fn_1   =  ( v_comp   * ( sk(k,j+1,i) + sk(k,j,i) )         &
     2876                         -ABS( v_comp ) * ( sk(k,j+1,i) - sk(k,j,i) )         &
     2877                             ) * adv_sca_1
     2878
     2879                   fd_1   = (  w(k-1,j,i)   * ( sk(k,j,i) + sk(k-1,j,i) )     &
     2880                         -ABS( w(k-1,j,i) ) * ( sk(k,j,i) - sk(k-1,j,i) )     &
     2881                            ) * adv_sca_1
     2882
     2883                   ft_1   = (  w(k,j,i)   * ( sk(k+1,j,i) + sk(k,j,i) )       &
     2884                         -ABS( w(k,j,i) ) * ( sk(k+1,j,i) - sk(k,j,i) )       &
     2885                            ) * adv_sca_1
     2886!
     2887!--                Calculate ratio of upwind gradients. Note, Min/Max is just
     2888!--                to avoid if statements.
     2889                   rl     = ( MAX( 0.0, u(k,j,i) - u_gtrans ) *               &
     2890                               ABS( ( sk(k,j,i-1) - sk(k,j,i-2)            ) /&
     2891                                    ( sk(k,j,i)   - sk(k,j,i-1) + 1E-20_wp )  &
     2892                                  ) +                                         &
     2893                              MIN( 0.0, u(k,j,i) - u_gtrans ) *               &
     2894                               ABS( ( sk(k,j,i)   - sk(k,j,i+1)            ) /&
     2895                                    ( sk(k,j,i-1) - sk(k,j,i)   + 1E-20_wp )  &
     2896                                  )                                           &
     2897                            ) / ABS( u(k,j,i) - u_gtrans + 1E-20_wp )
     2898
     2899                   rr     = ( MAX( 0.0, u(k,j,i+1) - u_gtrans ) *             &
     2900                               ABS( ( sk(k,j,i)   - sk(k,j,i-1)            ) /&
     2901                                    ( sk(k,j,i+1) - sk(k,j,i)   + 1E-20_wp )  &
     2902                                  ) +                                         &
     2903                              MIN( 0.0, u(k,j,i+1) - u_gtrans ) *             &
     2904                               ABS( ( sk(k,j,i+1) - sk(k,j,i+2)            ) /&
     2905                                    ( sk(k,j,i)   - sk(k,j,i+1) + 1E-20_wp )  &
     2906                                  )                                           &
     2907                            ) / ABS( u(k,j,i+1) - u_gtrans + 1E-20_wp )
     2908
     2909                   rs     = ( MAX( 0.0, v(k,j,i) - v_gtrans ) *               &
     2910                               ABS( ( sk(k,j-1,i) - sk(k,j-2,i)            ) /&
     2911                                    ( sk(k,j,i)   - sk(k,j-1,i) + 1E-20_wp )  &
     2912                                  ) +                                         &
     2913                              MIN( 0.0, v(k,j,i) - v_gtrans ) *               &
     2914                               ABS( ( sk(k,j,i)   - sk(k,j+1,i)            ) /&
     2915                                    ( sk(k,j-1,i) - sk(k,j,i)   + 1E-20_wp )  &
     2916                                  )                                           &
     2917                            ) / ABS( v(k,j,i) - v_gtrans + 1E-20_wp )
     2918
     2919                   rn     = ( MAX( 0.0, v(k,j+1,i) - v_gtrans ) *             &
     2920                               ABS( ( sk(k,j,i)   - sk(k,j-1,i)            ) /&
     2921                                    ( sk(k,j+1,i) - sk(k,j,i)   + 1E-20_wp )  &
     2922                                  ) +                                         &
     2923                              MIN( 0.0, v(k,j+1,i) - v_gtrans ) *             &
     2924                               ABS( ( sk(k,j+1,i) - sk(k,j+2,i)            ) /&
     2925                                    ( sk(k,j,i)   - sk(k,j+1,i) + 1E-20_wp )  &
     2926                                  )                                           &
     2927                            ) / ABS( v(k,j+1,i) - v_gtrans + 1E-20_wp )     
     2928!   
     2929!--                Reuse k_mm and compute k_mmm for the vertical gradient ratios.
     2930!--                Note, for vertical advection below the third grid point above
     2931!--                surface ( or below the model top) rd and rt are set to 0, i.e.
     2932!--                use of first order scheme is enforced.
     2933                   k_mmm  = k - 3 * ibit8
     2934
     2935                   rd     = ( MAX( 0.0, w(k-1,j,i) ) *                        &
     2936                            ABS( ( sk(k_mm,j,i) - sk(k_mmm,j,i)           ) / &
     2937                                 ( sk(k-1,j,i)  - sk(k_mm,j,i) + 1E-20_wp )   &
     2938                               ) +                                            &
     2939                              MIN( 0.0, w(k-1,j,i) ) *                        &
     2940                            ABS( ( sk(k-1,j,i) - sk(k,j,i)            ) /     &
     2941                                 ( sk(k_mm,j,i) - sk(k-1,j,i)   + 1E-20_wp )  &
     2942                               )                                              &
     2943                            ) * ibit8 / ABS( w(k-1,j,i) + 1E-20_wp )
     2944 
     2945                   rt     = ( MAX( 0.0, w(k,j,i) ) *                          &
     2946                            ABS( ( sk(k,j,i)   - sk(k-1,j,i)            ) /   &
     2947                                 ( sk(k+1,j,i) - sk(k,j,i)   + 1E-20_wp )     &
     2948                               ) +                                            &
     2949                              MIN( 0.0, w(k,j,i) ) *                          &
     2950                            ABS( ( sk(k+1,j,i) - sk(k_pp,j,i)           ) /   &
     2951                                 ( sk(k,j,i)   - sk(k+1,j,i) + 1E-20_wp )     &
     2952                               )                                              &
     2953                            ) * ibit8 / ABS( w(k,j,i) + 1E-20_wp )
     2954!
     2955!--                Calculate empirical limiter function (van Albada2 limiter).
     2956                   phi_l = MIN( 1.0_wp, ( 2.0_wp * rl ) /                     &
     2957                                        ( rl**2 + 1.0_wp ) )
     2958                   phi_r = MIN( 1.0_wp, ( 2.0_wp * rr ) /                     &
     2959                                        ( rr**2 + 1.0_wp ) )
     2960                   phi_s = MIN( 1.0_wp, ( 2.0_wp * rs ) /                     &
     2961                                        ( rs**2 + 1.0_wp ) )
     2962                   phi_n = MIN( 1.0_wp, ( 2.0_wp * rn ) /                     &
     2963                                        ( rn**2 + 1.0_wp ) )
     2964                   phi_d = MIN( 1.0_wp, ( 2.0_wp * rd ) /                     &
     2965                                        ( rd**2 + 1.0_wp ) )
     2966                   phi_t = MIN( 1.0_wp, ( 2.0_wp * rt ) /                     &
     2967                                        ( rt**2 + 1.0_wp ) )
     2968!
     2969!--                Calculate the resulting monotone flux.
     2970                   swap_flux_x_local(k,j) = fl_1 - phi_l *                    &
     2971                                          ( fl_1 - swap_flux_x_local(k,j) )
     2972                   flux_r(k)              = fr_1 - phi_r *                    &
     2973                                          ( fr_1 - flux_r(k)              )
     2974                   swap_flux_y_local(k)   = fs_1 - phi_s *                    &
     2975                                          ( fs_1 - swap_flux_y_local(k)   )
     2976                   flux_n(k)              = fn_1 - phi_n *                    &
     2977                                          ( fn_1 - flux_n(k)              )
     2978                   flux_d                 = fd_1 - phi_d *                    &
     2979                                          ( fd_1 - flux_d                 )
     2980                   flux_t(k)              = ft_1 - phi_t *                    &
     2981                                          ( ft_1 - flux_t(k)              )
     2982!
     2983!--                Moreover, modify dissipation flux according to the limiter.
     2984                   swap_diss_x_local(k,j) = swap_diss_x_local(k,j) * phi_l
     2985                   diss_r(k)              = diss_r(k)              * phi_r
     2986                   swap_diss_y_local(k)   = swap_diss_y_local(k)   * phi_s
     2987                   diss_n(k)              = diss_n(k)              * phi_n
     2988                   diss_d                 = diss_d                 * phi_d
     2989                   diss_t(k)              = diss_t(k)              * phi_t
     2990
     2991                ENDIF
     2992!
     2993!--             Calculate the divergence of the velocity field. A respective
     2994!--             correction is needed to overcome numerical instabilities caused
     2995!--             by a not sufficient reduction of divergences near topography.
     2996                div         =   ( u(k,j,i+1) - u(k,j,i)   ) * ddx             &
     2997                              + ( v(k,j+1,i) - v(k,j,i)   ) * ddy             &
     2998                              + ( w(k,j,i)   - w(k-1,j,i) ) * ddzw(k)
     2999
     3000                tend(k,j,i) = tend(k,j,i) - (                                 &
     3001                        ( flux_r(k) + diss_r(k) - swap_flux_x_local(k,j) -    &
     3002                          swap_diss_x_local(k,j)            ) * ddx           &
     3003                      + ( flux_n(k) + diss_n(k) - swap_flux_y_local(k)   -    &
     3004                          swap_diss_y_local(k)              ) * ddy           &
     3005                      + ( flux_t(k) + diss_t(k) - flux_d - diss_d             &
     3006                                                               ) * ddzw(k)    &
     3007                                            ) + sk(k,j,i) * div
     3008
     3009                swap_flux_y_local(k)   = flux_n(k)
     3010                swap_diss_y_local(k)   = diss_n(k)
     3011                swap_flux_x_local(k,j) = flux_r(k)
     3012                swap_diss_x_local(k,j) = diss_r(k)
     3013                flux_d                 = flux_t(k)
     3014                diss_d                 = diss_t(k)
     3015
     3016             ENDDO
     3017
     3018             DO  k = nzb_max+1, nzt
     3019
     3020                u_comp    = u(k,j,i+1) - u_gtrans
     3021                flux_r(k) = u_comp * (                                        &
     3022                      37.0_wp * ( sk(k,j,i+1) + sk(k,j,i)   )                 &
     3023                    -  8.0_wp * ( sk(k,j,i+2) + sk(k,j,i-1) )                 &
     3024                    +           ( sk(k,j,i+3) + sk(k,j,i-2) ) ) * adv_sca_5
     3025                diss_r(k) = -ABS( u_comp ) * (                                &
     3026                      10.0_wp * ( sk(k,j,i+1) - sk(k,j,i)   )                 &
     3027                    -  5.0_wp * ( sk(k,j,i+2) - sk(k,j,i-1) )                 &
     3028                    +           ( sk(k,j,i+3) - sk(k,j,i-2) ) ) * adv_sca_5
     3029
     3030                v_comp    = v(k,j+1,i) - v_gtrans
     3031                flux_n(k) = v_comp * (                                        &
     3032                      37.0_wp * ( sk(k,j+1,i) + sk(k,j,i)   )                 &
     3033                    -  8.0_wp * ( sk(k,j+2,i) + sk(k,j-1,i) )                 &
     3034                    +           ( sk(k,j+3,i) + sk(k,j-2,i) ) ) * adv_sca_5
     3035                diss_n(k) = -ABS( v_comp ) * (                                &
     3036                      10.0_wp * ( sk(k,j+1,i) - sk(k,j,i)   )                 &
     3037                    -  5.0_wp * ( sk(k,j+2,i) - sk(k,j-1,i) )                 &
     3038                    +           ( sk(k,j+3,i) - sk(k,j-2,i) ) ) * adv_sca_5
     3039!
     3040!--             k index has to be modified near bottom and top, else array
     3041!--             subscripts will be exceeded.
     3042                ibit8 = IBITS(wall_flags_0(k,j,i),8,1)
     3043                ibit7 = IBITS(wall_flags_0(k,j,i),7,1)
     3044                ibit6 = IBITS(wall_flags_0(k,j,i),6,1)
     3045
     3046                k_ppp = k + 3 * ibit8
     3047                k_pp  = k + 2 * ( 1 - ibit6  )
     3048                k_mm  = k - 2 * ibit8
     3049
     3050
     3051                flux_t(k) = w(k,j,i) * (                                      &
     3052                           ( 37.0_wp * ibit8 * adv_sca_5                      &
     3053                        +     7.0_wp * ibit7 * adv_sca_3                      &
     3054                        +              ibit6 * adv_sca_1                      &
     3055                           ) *                                                &
     3056                                   ( sk(k+1,j,i)  + sk(k,j,i)     )           &
     3057                    -      (  8.0_wp * ibit8 * adv_sca_5                      &
     3058                        +              ibit7 * adv_sca_3                      &
     3059                           ) *                                                &
     3060                                   ( sk(k_pp,j,i) + sk(k-1,j,i)   )           &
     3061                    +      (           ibit8 * adv_sca_5                      &
     3062                           ) *     ( sk(k_ppp,j,i)+ sk(k_mm,j,i)  )           &
     3063                                       )
     3064
     3065                diss_t(k) = -ABS( w(k,j,i) ) * (                              &
     3066                           ( 10.0_wp * ibit8 * adv_sca_5                      &
     3067                        +     3.0_wp * ibit7 * adv_sca_3                      &
     3068                        +              ibit6 * adv_sca_1                      &
     3069                           ) *                                                &
     3070                                   ( sk(k+1,j,i)   - sk(k,j,i)    )           &
     3071                    -      (  5.0_wp * ibit8 * adv_sca_5                      &
     3072                        +              ibit7 * adv_sca_3                      &
     3073                           ) *                                                &
     3074                                   ( sk(k_pp,j,i)  - sk(k-1,j,i)  )           &
     3075                    +      (           ibit8 * adv_sca_5                      &
     3076                           ) *                                                &
     3077                                   ( sk(k_ppp,j,i) - sk(k_mm,j,i) )           &
     3078                                               )
     3079!
     3080!--             Apply monotonic adjustment.
     3081                IF ( monotonic_adjustment )  THEN
     3082!
     3083!--                At first, calculate first order fluxes.
     3084                   u_comp = u(k,j,i) - u_gtrans
     3085                   fl_1   =  ( u_comp   * ( sk(k,j,i) + sk(k,j,i-1) )         &
     3086                         -ABS( u_comp ) * ( sk(k,j,i) - sk(k,j,i-1) )         &
     3087                             ) * adv_sca_1
     3088
     3089                   u_comp = u(k,j,i+1) - u_gtrans
     3090                   fr_1   =  ( u_comp   * ( sk(k,j,i+1) + sk(k,j,i) )         &
     3091                         -ABS( u_comp ) * ( sk(k,j,i+1) - sk(k,j,i) )         &
     3092                             ) * adv_sca_1
     3093
     3094                   v_comp = v(k,j,i) - v_gtrans
     3095                   fs_1   =  ( v_comp   * ( sk(k,j,i) + sk(k,j-1,i) )         &
     3096                         -ABS( v_comp ) * ( sk(k,j,i) - sk(k,j-1,i) )         &
     3097                             ) * adv_sca_1
     3098
     3099                   v_comp = v(k,j+1,i) - v_gtrans
     3100                   fn_1   =  ( v_comp   * ( sk(k,j+1,i) + sk(k,j,i) )         &
     3101                         -ABS( v_comp ) * ( sk(k,j+1,i) - sk(k,j,i) )         &
     3102                             ) * adv_sca_1
     3103
     3104                   fd_1   = (  w(k-1,j,i)   * ( sk(k,j,i) + sk(k-1,j,i) )     &
     3105                         -ABS( w(k-1,j,i) ) * ( sk(k,j,i) - sk(k-1,j,i) )     &
     3106                            ) * adv_sca_1
     3107
     3108                   ft_1   = (  w(k,j,i)   * ( sk(k+1,j,i) + sk(k,j,i) )       &
     3109                         -ABS( w(k,j,i) ) * ( sk(k+1,j,i) - sk(k,j,i) )       &
     3110                            ) * adv_sca_1
     3111!
     3112!--                Calculate ratio of upwind gradients. Note, Min/Max is just
     3113!--                to avoid if statements.
     3114                   rl     = ( MAX( 0.0, u(k,j,i) - u_gtrans ) *               &
     3115                               ABS( ( sk(k,j,i-1) - sk(k,j,i-2)            ) /&
     3116                                    ( sk(k,j,i)   - sk(k,j,i-1) + 1E-20_wp )  &
     3117                                  ) +                                         &
     3118                              MIN( 0.0, u(k,j,i) - u_gtrans ) *               &
     3119                               ABS( ( sk(k,j,i)   - sk(k,j,i+1)            ) /&
     3120                                    ( sk(k,j,i-1) - sk(k,j,i)   + 1E-20_wp )  &
     3121                                  )                                           &
     3122                            ) / ABS( u(k,j,i) - u_gtrans + 1E-20_wp )
     3123
     3124                   rr     = ( MAX( 0.0, u(k,j,i+1) - u_gtrans ) *             &
     3125                               ABS( ( sk(k,j,i)   - sk(k,j,i-1)            ) /&
     3126                                    ( sk(k,j,i+1) - sk(k,j,i)   + 1E-20_wp )  &
     3127                                  ) +                                         &
     3128                              MIN( 0.0, u(k,j,i+1) - u_gtrans ) *             &
     3129                               ABS( ( sk(k,j,i+1) - sk(k,j,i+2)            ) /&
     3130                                    ( sk(k,j,i)   - sk(k,j,i+1) + 1E-20_wp )  &
     3131                                  )                                           &
     3132                            ) / ABS( u(k,j,i+1) - u_gtrans + 1E-20_wp )
     3133
     3134                   rs     = ( MAX( 0.0, v(k,j,i) - v_gtrans ) *               &
     3135                               ABS( ( sk(k,j-1,i) - sk(k,j-2,i)            ) /&
     3136                                    ( sk(k,j,i)   - sk(k,j-1,i) + 1E-20_wp )  &
     3137                                  ) +                                         &
     3138                              MIN( 0.0, v(k,j,i) - v_gtrans ) *               &
     3139                               ABS( ( sk(k,j,i)   - sk(k,j+1,i)            ) /&
     3140                                    ( sk(k,j-1,i) - sk(k,j,i)   + 1E-20_wp )  &
     3141                                  )                                           &
     3142                            ) / ABS( v(k,j,i) - v_gtrans + 1E-20_wp )
     3143
     3144                   rn     = ( MAX( 0.0, v(k,j+1,i) - v_gtrans ) *             &
     3145                               ABS( ( sk(k,j,i)   - sk(k,j-1,i)            ) /&
     3146                                    ( sk(k,j+1,i) - sk(k,j,i)   + 1E-20_wp )  &
     3147                                  ) +                                         &
     3148                              MIN( 0.0, v(k,j+1,i) - v_gtrans ) *             &
     3149                               ABS( ( sk(k,j+1,i) - sk(k,j+2,i)            ) /&
     3150                                    ( sk(k,j,i)   - sk(k,j+1,i) + 1E-20_wp )  &
     3151                                  )                                           &
     3152                            ) / ABS( v(k,j+1,i) - v_gtrans + 1E-20_wp )     
     3153!   
     3154!--                Reuse k_mm and compute k_mmm for the vertical gradient ratios.
     3155!--                Note, for vertical advection below the third grid point above
     3156!--                surface ( or below the model top) rd and rt are set to 0, i.e.
     3157!--                use of first order scheme is enforced.
     3158                   k_mmm  = k - 3 * ibit8
     3159
     3160                   rd     = ( MAX( 0.0, w(k-1,j,i) ) *                        &
     3161                            ABS( ( sk(k_mm,j,i) - sk(k_mmm,j,i)           ) / &
     3162                                 ( sk(k-1,j,i)  - sk(k_mm,j,i) + 1E-20_wp )   &
     3163                               ) +                                            &
     3164                              MIN( 0.0, w(k-1,j,i) ) *                        &
     3165                            ABS( ( sk(k-1,j,i) - sk(k,j,i)            ) /     &
     3166                                 ( sk(k_mm,j,i) - sk(k-1,j,i)   + 1E-20_wp )  &
     3167                               )                                              &
     3168                            ) * ibit8 / ABS( w(k-1,j,i) + 1E-20_wp )
     3169 
     3170                   rt     = ( MAX( 0.0, w(k,j,i) ) *                          &
     3171                            ABS( ( sk(k,j,i)   - sk(k-1,j,i)            ) /   &
     3172                                 ( sk(k+1,j,i) - sk(k,j,i)   + 1E-20_wp )     &
     3173                               ) +                                            &
     3174                              MIN( 0.0, w(k,j,i) ) *                          &
     3175                            ABS( ( sk(k+1,j,i) - sk(k_pp,j,i)           ) /   &
     3176                                 ( sk(k,j,i)   - sk(k+1,j,i) + 1E-20_wp )     &
     3177                               )                                              &
     3178                            ) * ibit8 / ABS( w(k,j,i) + 1E-20_wp )
     3179!
     3180!--                Calculate empirical limiter function (van Albada2 limiter).
     3181                   phi_l = MIN( 1.0_wp, ( 2.0_wp * rl ) /                     &
     3182                                        ( rl**2 + 1.0_wp ) )
     3183                   phi_r = MIN( 1.0_wp, ( 2.0_wp * rr ) /                     &
     3184                                        ( rr**2 + 1.0_wp ) )
     3185                   phi_s = MIN( 1.0_wp, ( 2.0_wp * rs ) /                     &
     3186                                        ( rs**2 + 1.0_wp ) )
     3187                   phi_n = MIN( 1.0_wp, ( 2.0_wp * rn ) /                     &
     3188                                        ( rn**2 + 1.0_wp ) )
     3189                   phi_d = MIN( 1.0_wp, ( 2.0_wp * rd ) /                     &
     3190                                        ( rd**2 + 1.0_wp ) )
     3191                   phi_t = MIN( 1.0_wp, ( 2.0_wp * rt ) /                     &
     3192                                        ( rt**2 + 1.0_wp ) )
     3193!
     3194!--                Calculate the resulting monotone flux.
     3195                   swap_flux_x_local(k,j) = fl_1 - phi_l *                    &
     3196                                          ( fl_1 - swap_flux_x_local(k,j) )
     3197                   flux_r(k)              = fr_1 - phi_r *                    &
     3198                                          ( fr_1 - flux_r(k)              )
     3199                   swap_flux_y_local(k)   = fs_1 - phi_s *                    &
     3200                                          ( fs_1 - swap_flux_y_local(k)   )
     3201                   flux_n(k)              = fn_1 - phi_n *                    &
     3202                                          ( fn_1 - flux_n(k)              )
     3203                   flux_d                 = fd_1 - phi_d *                    &
     3204                                          ( fd_1 - flux_d                 )
     3205                   flux_t(k)              = ft_1 - phi_t *                    &
     3206                                          ( ft_1 - flux_t(k)              )
     3207!
     3208!--                Moreover, modify dissipation flux according to the limiter.
     3209                   swap_diss_x_local(k,j) = swap_diss_x_local(k,j) * phi_l
     3210                   diss_r(k)              = diss_r(k)              * phi_r
     3211                   swap_diss_y_local(k)   = swap_diss_y_local(k)   * phi_s
     3212                   diss_n(k)              = diss_n(k)              * phi_n
     3213                   diss_d                 = diss_d                 * phi_d
     3214                   diss_t(k)              = diss_t(k)              * phi_t
     3215
     3216                ENDIF
     3217!
     3218!--             Calculate the divergence of the velocity field. A respective
     3219!--             correction is needed to overcome numerical instabilities introduced
     3220!--             by a not sufficient reduction of divergences near topography.
     3221                div         =   ( u(k,j,i+1) - u(k,j,i)   ) * ddx             &
     3222                              + ( v(k,j+1,i) - v(k,j,i)   ) * ddy             &
     3223                              + ( w(k,j,i)   - w(k-1,j,i) ) * ddzw(k)
     3224
     3225                tend(k,j,i) = tend(k,j,i) - (                                 &
     3226                        ( flux_r(k) + diss_r(k) - swap_flux_x_local(k,j) -    &
     3227                          swap_diss_x_local(k,j)            ) * ddx           &
     3228                      + ( flux_n(k) + diss_n(k) - swap_flux_y_local(k)   -    &
     3229                          swap_diss_y_local(k)              ) * ddy           &
     3230                      + ( flux_t(k) + diss_t(k) - flux_d - diss_d             &
     3231                                                               ) * ddzw(k)    &
     3232                                            ) + sk(k,j,i) * div
     3233
     3234                swap_flux_y_local(k)   = flux_n(k)
     3235                swap_diss_y_local(k)   = diss_n(k)
     3236                swap_flux_x_local(k,j) = flux_r(k)
     3237                swap_diss_x_local(k,j) = diss_r(k)
     3238                flux_d                 = flux_t(k)
     3239                diss_d                 = diss_t(k)
     3240
     3241             ENDDO
     3242!
     3243!--          Evaluation of statistics. Please note, in case of using monotone
     3244!--          limiter vertical fluxes will be not calculated by the advective
     3245!--          fluxes.
     3246             IF ( .NOT. monotonic_adjustment )  THEN
     3247
     3248                SELECT CASE ( sk_char )
     3249
     3250                    CASE ( 'pt' )
     3251                       DO  k = nzb, nzt
     3252                          sums_wspts_ws_l(k,tn) = sums_wspts_ws_l(k,tn)       &
     3253                           + ( flux_t(k) + diss_t(k) )                        &
     3254                           *   weight_substep(intermediate_timestep_count)
     3255                       ENDDO
     3256                    CASE ( 'sa' )
     3257                       DO  k = nzb, nzt
     3258                          sums_wssas_ws_l(k,tn) = sums_wssas_ws_l(k,tn)       &
     3259                           + ( flux_t(k) + diss_t(k) )                        &
     3260                           *   weight_substep(intermediate_timestep_count)
     3261                       ENDDO
     3262                    CASE ( 'q' )
     3263                       DO  k = nzb, nzt
     3264                          sums_wsqs_ws_l(k,tn) = sums_wsqs_ws_l(k,tn)         &
     3265                           + ( flux_t(k) + diss_t(k) )                        &
     3266                           *   weight_substep(intermediate_timestep_count)
     3267                       ENDDO
     3268                    CASE ( 'qr' )
     3269                       DO  k = nzb, nzt
     3270                          sums_wsqrs_ws_l(k,tn) = sums_wsqrs_ws_l(k,tn)       &
     3271                           + ( flux_t(k) + diss_t(k) )                        &
     3272                           *   weight_substep(intermediate_timestep_count)
     3273                       ENDDO
     3274                    CASE ( 'nr' )
     3275                       DO  k = nzb, nzt
     3276                          sums_wsnrs_ws_l(k,tn) = sums_wsnrs_ws_l(k,tn)       &
     3277                           + ( flux_t(k) + diss_t(k) )                        &
     3278                           *   weight_substep(intermediate_timestep_count)
     3279                       ENDDO
     3280
     3281                 END SELECT
     3282
     3283              ENDIF
     3284
     3285         ENDDO
     3286      ENDDO
     3287
     3288    END SUBROUTINE advec_s_ws
     3289
     3290
     3291!------------------------------------------------------------------------------!
     3292! Scalar advection - Call for all grid points - accelerator version
     3293!------------------------------------------------------------------------------!
     3294    SUBROUTINE advec_s_ws_acc ( sk, sk_char )
     3295
     3296       USE arrays_3d,                                                         &
     3297           ONLY:  ddzw, tend, u, v, w
     3298
     3299       USE constants,                                                         &
     3300           ONLY:  adv_sca_1, adv_sca_3, adv_sca_5
     3301
     3302       USE control_parameters,                                                &
     3303           ONLY:  intermediate_timestep_count, monotonic_adjustment, u_gtrans,&
     3304                  v_gtrans
     3305
     3306       USE grid_variables,                                                    &
     3307           ONLY:  ddx, ddy
     3308
     3309       USE indices,                                                           &
     3310           ONLY:  i_left, i_right, j_north, j_south, nxlg, nxrg, nyng, nysg,  &
     3311                  nzb, nzb_max, nzt, wall_flags_0
     3312
     3313       USE kinds
     3314       
     3315!        USE statistics,                                                       &
     3316!            ONLY:  sums_wspts_ws_l, sums_wsqs_ws_l, sums_wssas_ws_l,          &
     3317!                   sums_wsqrs_ws_l, sums_wsnrs_ws_l, weight_substep
     3318
     3319       IMPLICIT NONE
     3320
     3321       CHARACTER (LEN = *), INTENT(IN)    :: sk_char !:
     3322
     3323       INTEGER(iwp) ::  i      !:
     3324       INTEGER(iwp) ::  ibit0  !:
     3325       INTEGER(iwp) ::  ibit1  !:
     3326       INTEGER(iwp) ::  ibit2  !:
     3327       INTEGER(iwp) ::  ibit3  !:
     3328       INTEGER(iwp) ::  ibit4  !:
     3329       INTEGER(iwp) ::  ibit5  !:
     3330       INTEGER(iwp) ::  ibit6  !:
     3331       INTEGER(iwp) ::  ibit7  !:
     3332       INTEGER(iwp) ::  ibit8  !:
     3333       INTEGER(iwp) ::  j      !:
     3334       INTEGER(iwp) ::  k      !:
     3335       INTEGER(iwp) ::  k_mm   !:
     3336       INTEGER(iwp) ::  k_mmm  !:
     3337       INTEGER(iwp) ::  k_pp   !:
     3338       INTEGER(iwp) ::  k_ppp  !:
     3339       INTEGER(iwp) ::  tn = 0 !:
    27353340
    27363341       REAL(wp)    ::  diss_d !:
     
    27473352       REAL(wp)    ::  flux_s !:
    27483353       REAL(wp)    ::  flux_t !:
     3354       REAL(wp)    ::  fd_1   !:
     3355       REAL(wp)    ::  fl_1   !:
     3356       REAL(wp)    ::  fn_1   !:
     3357       REAL(wp)    ::  fr_1   !:
     3358       REAL(wp)    ::  fs_1   !:
     3359       REAL(wp)    ::  ft_1   !:
     3360       REAL(wp)    ::  phi_d  !:
     3361       REAL(wp)    ::  phi_l  !:
     3362       REAL(wp)    ::  phi_n  !:
     3363       REAL(wp)    ::  phi_r  !:
     3364       REAL(wp)    ::  phi_s  !:
     3365       REAL(wp)    ::  phi_t  !:
     3366       REAL(wp)    ::  rd     !:
     3367       REAL(wp)    ::  rl     !:
     3368       REAL(wp)    ::  rn     !:
     3369       REAL(wp)    ::  rr     !:
     3370       REAL(wp)    ::  rs     !:
     3371       REAL(wp)    ::  rt     !:
    27493372       REAL(wp)    ::  u_comp !:
    27503373       REAL(wp)    ::  v_comp !:
     
    29673590                                   ( sk(k_ppp,j,i) - sk(k_mm,j,i) )           &
    29683591                                         )
     3592!
     3593!--             Apply monotonic adjustment.
     3594                IF ( monotonic_adjustment )  THEN
     3595!
     3596!--                At first, calculate first order fluxes.
     3597                   u_comp = u(k,j,i) - u_gtrans
     3598                   fl_1   =  ( u_comp   * ( sk(k,j,i) + sk(k,j,i-1) )         &
     3599                         -ABS( u_comp ) * ( sk(k,j,i) - sk(k,j,i-1) )         &
     3600                             ) * adv_sca_1
     3601
     3602                   u_comp = u(k,j,i+1) - u_gtrans
     3603                   fr_1   =  ( u_comp   * ( sk(k,j,i+1) + sk(k,j,i) )         &
     3604                         -ABS( u_comp ) * ( sk(k,j,i+1) - sk(k,j,i) )         &
     3605                             ) * adv_sca_1
     3606
     3607                   v_comp = v(k,j,i) - v_gtrans
     3608                   fs_1   =  ( v_comp   * ( sk(k,j,i) + sk(k,j-1,i) )         &
     3609                         -ABS( v_comp ) * ( sk(k,j,i) - sk(k,j-1,i) )         &
     3610                             ) * adv_sca_1
     3611
     3612                   v_comp = v(k,j+1,i) - v_gtrans
     3613                   fn_1   =  ( v_comp   * ( sk(k,j+1,i) + sk(k,j,i) )         &
     3614                         -ABS( v_comp ) * ( sk(k,j+1,i) - sk(k,j,i) )         &
     3615                             ) * adv_sca_1
     3616
     3617                   fd_1   = (  w(k-1,j,i)   * ( sk(k,j,i) + sk(k-1,j,i) )     &
     3618                        -ABS( w(k-1,j,i) ) * ( sk(k,j,i) - sk(k-1,j,i) )      &
     3619                            ) * adv_sca_1
     3620
     3621                   ft_1   = (  w(k,j,i)   * ( sk(k+1,j,i) + sk(k,j,i) )       &
     3622                        -ABS( w(k,j,i) ) * ( sk(k+1,j,i) - sk(k,j,i) )        &
     3623                            ) * adv_sca_1
     3624!
     3625!--                Calculate ratio of upwind gradients. Note, Min/Max is just
     3626!--                to avoid if statements.
     3627                   rl     = ( MAX( 0.0, u(k,j,i) - u_gtrans ) *               &
     3628                               ABS( ( sk(k,j,i-1) - sk(k,j,i-2)            ) /&
     3629                                    ( sk(k,j,i)   - sk(k,j,i-1) + 1E-20_wp )  &
     3630                                  ) +                                         &
     3631                              MIN( 0.0, u(k,j,i) - u_gtrans ) *               &
     3632                               ABS( ( sk(k,j,i)   - sk(k,j,i+1)            ) /&
     3633                                    ( sk(k,j,i-1) - sk(k,j,i)   + 1E-20_wp )  &
     3634                                  )                                           &
     3635                            ) / ABS( u(k,j,i) - u_gtrans + 1E-20_wp )
     3636
     3637                   rr     = ( MAX( 0.0, u(k,j,i+1) - u_gtrans ) *             &
     3638                               ABS( ( sk(k,j,i)   - sk(k,j,i-1)            ) /&
     3639                                    ( sk(k,j,i+1) - sk(k,j,i)   + 1E-20_wp )  &
     3640                                  ) +                                         &
     3641                              MIN( 0.0, u(k,j,i+1) - u_gtrans ) *             &
     3642                               ABS( ( sk(k,j,i+1) - sk(k,j,i+2)            ) /&
     3643                                    ( sk(k,j,i)   - sk(k,j,i+1) + 1E-20_wp )  &
     3644                                  )                                           &
     3645                            ) / ABS( u(k,j,i+1) - u_gtrans + 1E-20_wp )
     3646
     3647                   rs     = ( MAX( 0.0, v(k,j,i) - v_gtrans ) *               &
     3648                               ABS( ( sk(k,j-1,i) - sk(k,j-2,i)            ) /&
     3649                                    ( sk(k,j,i)   - sk(k,j-1,i) + 1E-20_wp )  &
     3650                                  ) +                                         &
     3651                              MIN( 0.0, v(k,j,i) - v_gtrans ) *               &
     3652                               ABS( ( sk(k,j,i)   - sk(k,j+1,i)            ) /&
     3653                                    ( sk(k,j-1,i) - sk(k,j,i)   + 1E-20_wp )  &
     3654                                  )                                           &
     3655                            ) / ABS( v(k,j,i) - v_gtrans + 1E-20_wp )
     3656
     3657                   rn     = ( MAX( 0.0, v(k,j+1,i) - v_gtrans ) *             &
     3658                               ABS( ( sk(k,j,i)   - sk(k,j-1,i)            ) /&
     3659                                    ( sk(k,j+1,i) - sk(k,j,i)   + 1E-20_wp )  &
     3660                                  ) +                                         &
     3661                              MIN( 0.0, v(k,j+1,i) - v_gtrans ) *             &
     3662                               ABS( ( sk(k,j+1,i) - sk(k,j+2,i)            ) /&
     3663                                    ( sk(k,j,i)   - sk(k,j+1,i) + 1E-20_wp )  &
     3664                                  )                                           &
     3665                            ) / ABS( v(k,j+1,i) - v_gtrans + 1E-20_wp )     
     3666!   
     3667!--                Reuse k_mm and compute k_mmm for the vertical gradient ratios.
     3668!--                Note, for vertical advection below the third grid point above
     3669!--                surface ( or below the model top) rd and rt are set to 0, i.e.
     3670!--                use of first order scheme is enforced.
     3671                   k_mmm  = k - 3 * ibit8
     3672
     3673                   rd     = ( MAX( 0.0, w(k-1,j,i) ) *                        &
     3674                            ABS( ( sk(k_mm,j,i) - sk(k_mmm,j,i)           ) / &
     3675                                 ( sk(k-1,j,i)  - sk(k_mm,j,i) + 1E-20_wp )   &
     3676                               ) +                                            &
     3677                              MIN( 0.0, w(k-1,j,i) ) *                        &
     3678                            ABS( ( sk(k-1,j,i) - sk(k,j,i)            ) /     &
     3679                                 ( sk(k_mm,j,i) - sk(k-1,j,i)   + 1E-20_wp )  &
     3680                               )                                              &
     3681                            ) * ibit8 / ABS( w(k-1,j,i) + 1E-20_wp )
     3682 
     3683                   rt     = ( MAX( 0.0, w(k,j,i) ) *                          &
     3684                            ABS( ( sk(k,j,i)   - sk(k-1,j,i)            ) /   &
     3685                                 ( sk(k+1,j,i) - sk(k,j,i)   + 1E-20_wp )     &
     3686                               ) +                                            &
     3687                              MIN( 0.0, w(k,j,i) ) *                          &
     3688                            ABS( ( sk(k+1,j,i) - sk(k_pp,j,i)           ) /   &
     3689                                 ( sk(k,j,i)   - sk(k+1,j,i) + 1E-20_wp )     &
     3690                               )                                              &
     3691                            ) * ibit8 / ABS( w(k,j,i) + 1E-20_wp )
     3692!
     3693!--                Calculate empirical limiter function (van Albada2 limiter).
     3694                   phi_l = MIN( 1.0_wp, ( 2.0_wp * rl ) /                     &
     3695                                        ( rl**2 + 1.0_wp ) )
     3696                   phi_r = MIN( 1.0_wp, ( 2.0_wp * rr ) /                     &
     3697                                        ( rr**2 + 1.0_wp ) )
     3698                   phi_s = MIN( 1.0_wp, ( 2.0_wp * rs ) /                     &
     3699                                        ( rs**2 + 1.0_wp ) )
     3700                   phi_n = MIN( 1.0_wp, ( 2.0_wp * rn ) /                     &
     3701                                        ( rn**2 + 1.0_wp ) )
     3702                   phi_d = MIN( 1.0_wp, ( 2.0_wp * rd ) /                     &
     3703                                        ( rd**2 + 1.0_wp ) )
     3704                   phi_t = MIN( 1.0_wp, ( 2.0_wp * rt ) /                     &
     3705                                        ( rt**2 + 1.0_wp ) )
     3706!
     3707!--                Calculate the resulting monotone flux.
     3708                   flux_l = fl_1 - phi_l * ( fl_1 - flux_l )
     3709                   flux_r = fr_1 - phi_r * ( fr_1 - flux_r )
     3710                   flux_s = fs_1 - phi_s * ( fs_1 - flux_s )
     3711                   flux_n = fn_1 - phi_n * ( fn_1 - flux_n )
     3712                   flux_d = fd_1 - phi_d * ( fd_1 - flux_d )
     3713                   flux_t = ft_1 - phi_t * ( ft_1 - flux_t )
     3714!
     3715!--                Moreover, modify dissipation flux according to the limiter.
     3716                   diss_l = diss_l * phi_l
     3717                   diss_r = diss_r * phi_r
     3718                   diss_s = diss_s * phi_s
     3719                   diss_n = diss_n * phi_n
     3720                   diss_d = diss_d * phi_d
     3721                   diss_t = diss_t * phi_t
     3722
     3723                ENDIF
    29693724!
    29703725!--             Calculate the divergence of the velocity field. A respective
  • palm/trunk/SOURCE/check_parameters.f90

    r1556 r1557  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! Added checks for monotonic limiter
    2323!
    2424! Former revisions:
     
    605605    IF ( topography /= 'flat' )  THEN
    606606       action = ' '
    607        IF ( scalar_advec /= 'pw-scheme' .AND. scalar_advec /= 'ws-scheme')  THEN
     607       IF ( scalar_advec /= 'pw-scheme' .AND. scalar_advec /= 'ws-scheme'     &
     608      .AND. scalar_advec /= 'ws-scheme-mono' )  THEN
    608609          WRITE( action, '(A,A)' )  'scalar_advec = ', scalar_advec
    609610       ENDIF
     
    751752       CALL message( 'check_parameters', 'PA0022', 1, 2, 0, 6, 0 )
    752753    ENDIF
    753     IF ( ( momentum_advec == 'ws-scheme' .OR.  scalar_advec == 'ws-scheme' )   &
     754    IF ( ( momentum_advec == 'ws-scheme' .OR.  scalar_advec == 'ws-scheme'     &
     755           .OR. scalar_advec == 'ws-scheme-mono' )                             &
    754756           .AND. ( timestep_scheme == 'euler' .OR.                             &
    755757                   timestep_scheme == 'runge-kutta-2' ) )                      &
     
    761763    ENDIF
    762764    IF ( scalar_advec /= 'pw-scheme'  .AND.  scalar_advec /= 'ws-scheme' .AND. &
    763          scalar_advec /= 'bc-scheme' )                                         &
     765         scalar_advec /= 'ws-scheme-mono' .AND. scalar_advec /= 'bc-scheme' )  &
    764766    THEN
    765767       message_string = 'unknown advection scheme: scalar_advec = "' // &
     
    776778
    777779    IF ( use_sgs_for_particles  .AND.  .NOT. use_upstream_for_tke .AND.        &
    778          scalar_advec /= 'ws-scheme' )  THEN
     780         ( scalar_advec /= 'ws-scheme' .OR.                                    &
     781           scalar_advec /= 'ws-scheme-mono' )                                  &
     782       )  THEN
    779783       use_upstream_for_tke = .TRUE.
    780784       message_string = 'use_upstream_for_tke set .TRUE. because ' //          &
     
    792796!
    793797!-- Set LOGICAL switches to enhance performance
    794     IF ( momentum_advec == 'ws-scheme' )    ws_scheme_mom = .TRUE.
    795     IF ( scalar_advec   == 'ws-scheme'   )  ws_scheme_sca = .TRUE.
     798    IF ( momentum_advec == 'ws-scheme' )       ws_scheme_mom = .TRUE.
     799    IF ( scalar_advec   == 'ws-scheme' .OR.                                   &
     800         scalar_advec   == 'ws-scheme-mono' )  ws_scheme_sca = .TRUE.
     801    IF ( scalar_advec   == 'ws-scheme-mono' )  monotonic_adjustment = .TRUE.
     802
    796803
    797804!
     
    17071714          CALL message( 'check_parameters', 'PA0051', 1, 2, 0, 6, 0 )
    17081715       ENDIF
    1709        IF ( momentum_advec /= 'pw-scheme' .AND. &
    1710             momentum_advec /= 'ws-scheme')  THEN
     1716       IF ( scalar_advec /= 'pw-scheme' .AND.                                 &
     1717          ( scalar_advec /= 'ws-scheme' .OR.                                  &
     1718            scalar_advec /= 'ws-scheme-mono' )                                &
     1719          )  THEN
     1720
    17111721          message_string = 'non-cyclic lateral boundaries do not allow ' // &
    17121722                           'momentum_advec = "' // TRIM( momentum_advec ) // '"'
    17131723          CALL message( 'check_parameters', 'PA0052', 1, 2, 0, 6, 0 )
    17141724       ENDIF
    1715        IF ( scalar_advec /= 'pw-scheme' .AND. &
    1716             scalar_advec /= 'ws-scheme' )  THEN
     1725       IF ( scalar_advec /= 'pw-scheme' .AND.                                  &
     1726          ( scalar_advec /= 'ws-scheme' .OR.                                   &
     1727            scalar_advec /= 'ws-scheme-mono' )                                 &
     1728          )  THEN
    17171729          message_string = 'non-cyclic lateral boundaries do not allow ' // &
    17181730                           'scalar_advec = "' // TRIM( scalar_advec ) // '"'
  • palm/trunk/SOURCE/flow_statistics.f90

    r1556 r1557  
    2121! Current revisions:
    2222! -----------------
    23 !
     23! Adjustments for monotonic limiter
    2424!
    2525! Former revisions:
     
    143143        ONLY:   average_count_pr, cloud_droplets, cloud_physics, do_sum,       &
    144144                dt_3d, g, humidity, icloud_scheme, kappa, large_scale_forcing, &
    145                 large_scale_subsidence, max_pr_user, message_string, ocean,    &
     145                large_scale_subsidence, max_pr_user, message_string,           &
     146                monotonic_adjustment, ocean,                                   &
    146147                passive_scalar, precipitation, simulated_time,                 &
    147148                use_subsidence_tendencies, use_surface_fluxes, use_top_fluxes, &
     
    281282       ENDIF
    282283
    283        IF ( ws_scheme_sca .AND. sr == 0 )  THEN
     284       IF ( ws_scheme_sca .AND. .NOT. monotonic_adjustment                    &
     285           .AND. sr == 0 )  THEN
    284286
    285287          DO  i = 0, threads_per_task-1
     
    802804!--             but so far there is no other suitable place to calculate)
    803805                IF ( ocean )  THEN
    804                    IF( .NOT. ws_scheme_sca .OR. sr /= 0 )  THEN
     806                   IF( .NOT. ws_scheme_sca .OR. monotonic_adjustment .OR.      &
     807                             sr /= 0 )  THEN
    805808                      pts = 0.5_wp * ( sa(k,j,i)   - hom(k,1,23,sr) +          &
    806809                                       sa(k+1,j,i) - hom(k+1,1,23,sr) )
     
    853856                      ENDIF
    854857                   ELSE
    855                       IF( .NOT. ws_scheme_sca .OR. sr /= 0 )  THEN
     858                      IF( .NOT. ws_scheme_sca .OR. monotonic_adjustment .OR.   &
     859                                sr /= 0 )  THEN
    856860                         pts = 0.5_wp * ( vpt(k,j,i)   - hom(k,1,44,sr) +      &
    857861                                          vpt(k+1,j,i) - hom(k+1,1,44,sr) )
    858862                         sums_l(k,46,tn) = sums_l(k,46,tn) + pts * w(k,j,i) *  &
    859863                                                             rmask(j,i,sr)
    860                       ELSE IF ( ws_scheme_sca .AND. sr == 0 )  THEN
     864                      ELSE IF ( ws_scheme_sca .AND. .NOT. monotonic_adjustment &
     865                                .AND. sr == 0 )  THEN
    861866                         sums_l(k,46,tn) = ( 1.0_wp + 0.61_wp *                &
    862867                                             hom(k,1,41,sr) ) *                &
     
    870875!--             Passive scalar flux
    871876                IF ( passive_scalar .AND. ( .NOT. ws_scheme_sca                &
    872                      .OR. sr /= 0 ) )  THEN
     877                     .OR. monotonic_adjustment .OR. sr /= 0 ) )  THEN
    873878                   pts = 0.5_wp * ( q(k,j,i)   - hom(k,1,41,sr) +              &
    874879                                    q(k+1,j,i) - hom(k+1,1,41,sr) )
     
    914919
    915920       ENDIF
    916        IF ( .NOT. ws_scheme_sca .OR. sr /= 0 )  THEN
     921       IF ( .NOT. ws_scheme_sca .OR. monotonic_adjustment .OR. sr /= 0 )  THEN
    917922         !$OMP DO
    918923         DO  i = nxl, nxr
     
    15251530        ONLY :  average_count_pr, cloud_droplets, cloud_physics, do_sum,       &
    15261531                dt_3d, g, humidity, icloud_scheme, kappa, large_scale_forcing, &
    1527                 large_scale_subsidence, max_pr_user, message_string, ocean,    &
     1532                large_scale_subsidence, max_pr_user, message_string,           &
     1533                monotonic_adjustment, ocean,                                   &
    15281534                passive_scalar, precipitation, simulated_time,                 &
    15291535                use_subsidence_tendencies, use_surface_fluxes, use_top_fluxes, &
     
    16631669       ENDIF
    16641670
    1665        IF ( ws_scheme_sca .AND. sr == 0 )  THEN
     1671       IF ( ws_scheme_sca .AND. .NOT. monotonic_adjustment .AND. sr == 0 )  THEN
    16661672
    16671673          DO  i = 0, threads_per_task-1
     
    24862492       IF ( ocean )  THEN
    24872493
    2488           IF( .NOT. ws_scheme_sca .OR. sr /= 0 )  THEN
     2494          IF( .NOT. ws_scheme_sca .OR. monotonic_adjustment .OR. sr /= 0 )  THEN
    24892495
    24902496             !$acc parallel loop gang present( hom, rflags_invers, rmask, sa, sums_l, w ) create( s1 )
     
    26282634          ELSE
    26292635
    2630              IF( .NOT. ws_scheme_sca .OR.  sr /= 0 )  THEN
     2636             IF( .NOT. ws_scheme_sca .OR. monotonic_adjustment .OR.  sr /= 0 )  THEN
    26312637
    26322638                !$acc parallel loop gang present( hom, rflags_invers, rmask, sums_l, vpt, w ) create( s1 )
     
    26602666!
    26612667!--    Passive scalar flux
    2662        IF ( passive_scalar  .AND.  ( .NOT. ws_scheme_sca  .OR.  sr /= 0 ) )  THEN
     2668       IF ( passive_scalar  .AND.  ( .NOT. ws_scheme_sca .OR. monotonic_adjustment &
     2669                                     .OR.  sr /= 0 ) )  THEN
    26632670
    26642671          !$acc parallel loop gang present( hom, q, rflags_invers, rmask, sums_l, w ) create( s1 )
     
    27142721       ENDIF
    27152722
    2716        IF ( .NOT. ws_scheme_sca .OR.  sr /= 0 )  THEN
     2723       IF ( .NOT. ws_scheme_sca .OR. monotonic_adjustment .OR.  sr /= 0 )  THEN
    27172724
    27182725          !$OMP DO
  • palm/trunk/SOURCE/header.f90

    r1552 r1557  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! output for monotonic limiter
    2323!
    2424! Former revisions:
     
    452452    ELSEIF ( scalar_advec == 'ws-scheme' )  THEN
    453453       WRITE ( io, 504 )
     454    ELSEIF ( scalar_advec == 'ws-scheme-mono' )  THEN
     455       WRITE ( io, 513 )
    454456    ELSE
    455457       WRITE ( io, 118 )
     
    23302332            ' Time:                 ',A8,6X,'Run-No.:   ',I2.2/     &
    23312333            ' Run on host:        ',A10,6X,'En-No.:    ',I2.2)
     2334513 FORMAT (' --> Scalar advection via Wicker-Skamarock-Scheme 5th order ' // &
     2335            '+ monotonic adjustment')
     2336
    23322337
    23332338 END SUBROUTINE header
  • palm/trunk/SOURCE/init_grid.f90

    r1419 r1557  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! Adjustment for monotoinic limiter
    2323!
    2424! Former revisions:
     
    11691169    wall_flags_00 = 0
    11701170
    1171     IF ( scalar_advec == 'ws-scheme' )  THEN
     1171    IF ( scalar_advec == 'ws-scheme' .OR.                                     &
     1172         scalar_advec == 'ws-scheme-mono' )  THEN
    11721173!
    11731174!--    Set flags to steer the degradation of the advection scheme in advec_ws
  • palm/trunk/SOURCE/init_pegrid.f90

    r1469 r1557  
    2020! Current revisions:
    2121! ------------------
    22 !
     22! Adjustment for monotonic limiter
    2323!
    2424! Former revisions:
     
    591591!
    592592!-- Determine the number of ghost point layers
    593     IF ( scalar_advec == 'ws-scheme' .OR. momentum_advec == 'ws-scheme' )  THEN
     593    IF ( scalar_advec == 'ws-scheme'      .OR.                                &
     594         scalar_advec == 'ws-scheme-mono' .OR.                                &
     595         momentum_advec == 'ws-scheme' )  THEN
    594596       nbgp = 3
    595597    ELSE
  • palm/trunk/SOURCE/modules.f90

    r1552 r1557  
    2020! Current revisions:
    2121! ------------------
    22 !
     22! +monotonic_adjustment
    2323!
    2424! Former revisions:
     
    659659                lunudge = .FALSE., lvnudge = .FALSE., lwnudge = .FALSE., &
    660660                masking_method = .FALSE., mg_switch_to_pe0 = .FALSE., &
     661                monotonic_adjustment = .FALSE., &
    661662                neutral = .FALSE., nudging = .FALSE., &
    662663                ocean = .FALSE., on_device = .FALSE., &
Note: See TracChangeset for help on using the changeset viewer.