Ignore:
Timestamp:
Oct 26, 2016 11:15:40 AM (5 years ago)
Author:
knoop
Message:

Anelastic approximation implemented

File:
1 edited

Legend:

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

    r2001 r2037  
    2020! Current revisions:
    2121! ------------------
    22 !
     22! Anelastic approximation implemented
    2323!
    2424! Former revisions:
     
    902902
    903903       USE arrays_3d,                                                          &
    904            ONLY:  ddzw, tend, u, v, w
     904           ONLY:  ddzw, drho_air, tend, u, v, w, rho_air, rho_air_zw
    905905
    906906       USE constants,                                                          &
     
    12111211
    12121212
    1213           flux_t(k) = w(k,j,i) * (                                            &
     1213          flux_t(k) = w(k,j,i) * rho_air_zw(k) * (                            &
    12141214                     ( 37.0_wp * ibit8 * adv_sca_5                            &
    12151215                  +     7.0_wp * ibit7 * adv_sca_3                            &
     
    12251225                                 )
    12261226
    1227           diss_t(k) = -ABS( w(k,j,i) ) * (                                    &
     1227          diss_t(k) = -ABS( w(k,j,i) ) * rho_air_zw(k) * (                    &
    12281228                     ( 10.0_wp * ibit8 * adv_sca_5                            &
    12291229                  +     3.0_wp * ibit7 * adv_sca_3                            &
     
    12661266             fd_1   =  ( w(k-1,j,i)   * ( sk(k,j,i) + sk(k-1,j,i) )           &
    12671267                   -ABS( w(k-1,j,i) ) * ( sk(k,j,i) - sk(k-1,j,i) )           &
    1268                        ) * adv_sca_1
     1268                       ) * adv_sca_1 * rho_air_zw(k)
    12691269
    12701270             ft_1   =  ( w(k,j,i)   * ( sk(k+1,j,i) + sk(k,j,i) )             &
    12711271                   -ABS( w(k,j,i) ) * ( sk(k+1,j,i) - sk(k,j,i) )             &
    1272                       ) * adv_sca_1
     1272                      ) * adv_sca_1 * rho_air_zw(k)
    12731273!
    12741274!--          Calculate ratio of upwind gradients. Note, Min/Max is just to
     
    13861386                                         + IBITS(wall_flags_0(k,j,i-1),2,1)    &
    13871387                                         )                                     &
    1388                           ) * ddx                                              &
     1388                          ) * rho_air(k) * ddx                                 &
    13891389                        + ( v(k,j+1,i) * ( ibit3 + ibit4 + ibit5 )             &
    13901390                          - v(k,j,i)   * ( IBITS(wall_flags_0(k,j-1,i),3,1)    &
     
    13921392                                         + IBITS(wall_flags_0(k,j-1,i),5,1)    &
    13931393                                         )                                     &
    1394                           ) * ddy                                              &
    1395                         + ( w(k,j,i)   * ( ibit6 + ibit7 + ibit8 )             &
    1396                           - w(k-1,j,i) * ( IBITS(wall_flags_0(k-1,j,i),6,1)    &
     1394                          ) * rho_air(k) * ddy                                 &
     1395                        + ( w(k,j,i) * rho_air_zw(k) *                         &
     1396                                         ( ibit6 + ibit7 + ibit8 )             &
     1397                          - w(k-1,j,i) * rho_air_zw(k-1) *                     &
     1398                                         ( IBITS(wall_flags_0(k-1,j,i),6,1)    &
    13971399                                         + IBITS(wall_flags_0(k-1,j,i),7,1)    &
    13981400                                         + IBITS(wall_flags_0(k-1,j,i),8,1)    &
     
    14061408                      + ( flux_n(k) + diss_n(k) - swap_flux_y_local(k,tn)   - &
    14071409                          swap_diss_y_local(k,tn)              ) * ddy        &
    1408                       + ( flux_t(k) + diss_t(k) - flux_d - diss_d             &
    1409                                                                ) * ddzw(k)    &
     1410                      + ( ( flux_t(k) + diss_t(k) ) -                         &
     1411                          ( flux_d    + diss_d    )                           &
     1412                                                    ) * drho_air(k) * ddzw(k) &
    14101413                                      ) + sk(k,j,i) * div
    14111414
     
    14541457          k_mm  = k - 2 * ibit8
    14551458
    1456           flux_t(k) = w(k,j,i) * (                                            &
     1459          flux_t(k) = w(k,j,i) * rho_air_zw(k) * (                            &
    14571460                    ( 37.0_wp * ibit8 * adv_sca_5                             &
    14581461                 +     7.0_wp * ibit7 * adv_sca_3                             &
     
    14681471                                 )
    14691472
    1470           diss_t(k) = -ABS( w(k,j,i) ) * (                                    &
     1473          diss_t(k) = -ABS( w(k,j,i) ) * rho_air_zw(k) * (                    &
    14711474                    ( 10.0_wp * ibit8 * adv_sca_5                             &
    14721475                 +     3.0_wp * ibit7 * adv_sca_3                             &
     
    15111514             fd_1   =  ( w(k-1,j,i)   * ( sk(k,j,i) + sk(k-1,j,i) )           &
    15121515                   -ABS( w(k-1,j,i) ) * ( sk(k,j,i) - sk(k-1,j,i) )           &
    1513                        ) * adv_sca_1
     1516                       ) * adv_sca_1  * rho_air_zw(k)
    15141517
    15151518             ft_1   =  ( w(k,j,i)   * ( sk(k+1,j,i) + sk(k,j,i) )             &
    15161519                   -ABS( w(k,j,i) ) * ( sk(k+1,j,i) - sk(k,j,i) )             &
    1517                        ) * adv_sca_1
     1520                       ) * adv_sca_1  * rho_air_zw(k)
    15181521!
    15191522!--          Calculate ratio of upwind gradients. Note, Min/Max is just to
     
    16261629!--       correction is needed to overcome numerical instabilities introduced
    16271630!--       by a not sufficient reduction of divergences near topography.
    1628           div         =   ( u(k,j,i+1) - u(k,j,i)   ) * ddx                   &
    1629                         + ( v(k,j+1,i) - v(k,j,i)   ) * ddy                   &
    1630                         + ( w(k,j,i)   - w(k-1,j,i) ) * ddzw(k)
     1631          div         =   ( u(k,j,i+1) - u(k,j,i)   ) * rho_air(k) * ddx      &
     1632                        + ( v(k,j+1,i) - v(k,j,i)   ) * rho_air(k) * ddy      &
     1633                        + ( w(k,j,i)   * rho_air_zw(k) -                      &
     1634                            w(k-1,j,i) * rho_air_zw(k-1) ) * ddzw(k)
    16311635
    16321636          tend(k,j,i) = tend(k,j,i) - (                                       &
     
    16351639                      + ( flux_n(k) + diss_n(k) - swap_flux_y_local(k,tn)   - &
    16361640                          swap_diss_y_local(k,tn)              ) * ddy        &
    1637                       + ( flux_t(k) + diss_t(k) - flux_d - diss_d             &
    1638                                                                ) * ddzw(k)    &
     1641                      + ( ( flux_t(k) + diss_t(k) ) -                         &
     1642                          ( flux_d    + diss_d    )                           &
     1643                                                    ) * drho_air(k) * ddzw(k) &
    16391644                                      ) + sk(k,j,i) * div
    16401645
     
    17341739
    17351740       USE arrays_3d,                                                         &
    1736            ONLY:  ddzw, diss_l_u, diss_s_u, flux_l_u, flux_s_u, tend, u, v, w
     1741           ONLY:  ddzw, diss_l_u, diss_s_u, flux_l_u, flux_s_u, tend, u, v, w,&
     1742                  drho_air, rho_air, rho_air_zw
    17371743
    17381744       USE constants,                                                         &
     
    19992005
    20002006          w_comp    = w(k,j,i) + w(k,j,i-1)
    2001           flux_t(k) = w_comp  * (                                             &
     2007          flux_t(k) = w_comp * rho_air_zw(k) * (                              &
    20022008                     ( 37.0_wp * ibit17 * adv_mom_5                           &
    20032009                  +     7.0_wp * ibit16 * adv_mom_3                           &
     
    20142020                                 )
    20152021
    2016           diss_t(k) = - ABS( w_comp ) * (                                     &
     2022          diss_t(k) = - ABS( w_comp ) * rho_air_zw(k) * (                     &
    20172023                     ( 10.0_wp * ibit17 * adv_mom_5                           &
    20182024                  +     3.0_wp * ibit16 * adv_mom_3                           &
     
    20372043                                      + IBITS(wall_flags_0(k,j,i-1),10,1)     &
    20382044                                      + IBITS(wall_flags_0(k,j,i-1),11,1)     &
    2039                                       )                                       &   
    2040                   ) * ddx                                                     &
     2045                                      )                                       &
     2046                  ) * rho_air(k) * ddx                                        &
    20412047               +  ( ( v_comp + gv ) * ( ibit12 + ibit13 + ibit14 )            &
    20422048                  - ( v(k,j,i)   + v(k,j,i-1 )  )                             &
     
    20452051                                      + IBITS(wall_flags_0(k,j-1,i),14,1)     &
    20462052                                      )                                       &
    2047                   ) * ddy                                                     &
    2048                +  ( w_comp          * ( ibit15 + ibit16 + ibit17 )            &
    2049                 - ( w(k-1,j,i) + w(k-1,j,i-1) )                               &
     2053                  ) * rho_air(k) * ddy                                        &
     2054               +  ( w_comp * rho_air_zw(k) * ( ibit15 + ibit16 + ibit17 )     &
     2055                - ( w(k-1,j,i) + w(k-1,j,i-1) ) * rho_air_zw(k-1)             &
    20502056                                    * ( IBITS(wall_flags_0(k-1,j,i),15,1)     &
    20512057                                      + IBITS(wall_flags_0(k-1,j,i),16,1)     &
     
    20612067                          + ( flux_n(k) + diss_n(k)                           &
    20622068                          -   flux_s_u(k,tn) - diss_s_u(k,tn)     ) * ddy     &
    2063                           + ( flux_t(k) + diss_t(k)                           &
    2064                           -   flux_d    - diss_d                  ) * ddzw(k) &
     2069                          + ( ( flux_t(k) + diss_t(k) )                       &
     2070                          -   ( flux_d    + diss_d )                          &
     2071                                                    ) * drho_air(k) * ddzw(k) &
    20652072                                       ) + div * u(k,j,i)
    20662073
     
    21272134
    21282135          w_comp    = w(k,j,i) + w(k,j,i-1)
    2129           flux_t(k) = w_comp  * (                                             &
     2136          flux_t(k) = w_comp * rho_air_zw(k) * (                              &
    21302137                     ( 37.0_wp * ibit17 * adv_mom_5                           &
    21312138                  +     7.0_wp * ibit16 * adv_mom_3                           &
     
    21422149                                 )
    21432150
    2144           diss_t(k) = - ABS( w_comp ) * (                                     &
     2151          diss_t(k) = - ABS( w_comp ) * rho_air_zw(k) * (                     &
    21452152                     ( 10.0_wp * ibit17 * adv_mom_5                           &
    21462153                  +     3.0_wp * ibit16 * adv_mom_3                           &
     
    21612168!--       by a not sufficient reduction of divergences near topography.
    21622169          div = ( ( u_comp(k)   - ( u(k,j,i)   + u(k,j,i-1)   ) ) * ddx       &
     2170                                                                  * rho_air(k)&
    21632171               +  ( v_comp + gv - ( v(k,j,i)   + v(k,j,i-1 )  ) ) * ddy       &
    2164                +  ( w_comp      - ( w(k-1,j,i) + w(k-1,j,i-1) ) ) * ddzw(k)   &
     2172                                                                  * rho_air(k)&
     2173               +  (   w_comp                      * rho_air_zw(k)   -         &
     2174                    ( w(k-1,j,i) + w(k-1,j,i-1) ) * rho_air_zw(k-1)           &
     2175                  ) * ddzw(k)                                                 &
    21652176                ) * 0.5_wp
    21662177
     
    21702181                          + ( flux_n(k) + diss_n(k)                           &
    21712182                          -   flux_s_u(k,tn) - diss_s_u(k,tn)     ) * ddy     &
    2172                           + ( flux_t(k) + diss_t(k)                           &
    2173                           -   flux_d    - diss_d                  ) * ddzw(k) &
     2183                          + ( ( flux_t(k) + diss_t(k) )                       &
     2184                          -   ( flux_d    + diss_d    )                       &
     2185                                                    ) * drho_air(k) * ddzw(k) &
    21742186                                       ) + div * u(k,j,i)
    21752187
     
    22192231
    22202232       USE arrays_3d,                                                          &
    2221            ONLY:  ddzw, diss_l_v, diss_s_v, flux_l_v, flux_s_v, tend, u, v, w
     2233           ONLY:  ddzw, diss_l_v, diss_s_v, flux_l_v, flux_s_v, tend, u, v, w, &
     2234                  drho_air, rho_air, rho_air_zw
    22222235
    22232236       USE constants,                                                          &
     
    24862499
    24872500          w_comp    = w(k,j-1,i) + w(k,j,i)
    2488           flux_t(k) = w_comp  * (                                             &
     2501          flux_t(k) = w_comp * rho_air_zw(k) * (                              &
    24892502                     ( 37.0_wp * ibit26 * adv_mom_5                           &
    24902503                  +     7.0_wp * ibit25 * adv_mom_3                           &
     
    25012514                                 )
    25022515
    2503           diss_t(k) = - ABS( w_comp ) * (                                     &
     2516          diss_t(k) = - ABS( w_comp ) * rho_air_zw(k) * (                     &
    25042517                     ( 10.0_wp * ibit26 * adv_mom_5                           &
    25052518                  +     3.0_wp * ibit25 * adv_mom_3                           &
     
    25252538                                         + IBITS(wall_flags_0(k,j,i-1),19,1)  &
    25262539                                         + IBITS(wall_flags_0(k,j,i-1),20,1)  &
    2527                                          )                                    &   
    2528                   ) * ddx                                                     &
     2540                                         )                                    &
     2541                  ) * rho_air(k) * ddx                                        &
    25292542               +  ( v_comp(k)                                                 &
    25302543                                       * ( ibit21 + ibit22 + ibit23 )         &
     
    25332546                                         + IBITS(wall_flags_0(k,j-1,i),22,1)  &
    25342547                                         + IBITS(wall_flags_0(k,j-1,i),23,1)  &
    2535                                          )                                    &   
    2536                   ) * ddy                                                     &
    2537                +  ( w_comp                                                    &
    2538                                        * ( ibit24 + ibit25 + ibit26 )         &
    2539                 - ( w(k-1,j-1,i) + w(k-1,j,i) )                               &
     2548                                         )                                    &
     2549                  ) * rho_air(k) * ddy                                        &
     2550               +  ( w_comp * rho_air_zw(k) * ( ibit24 + ibit25 + ibit26 )     &
     2551                - ( w(k-1,j-1,i) + w(k-1,j,i) ) * rho_air_zw(k-1)             &
    25402552                                       * ( IBITS(wall_flags_0(k-1,j,i),24,1)  &
    25412553                                         + IBITS(wall_flags_0(k-1,j,i),25,1)  &
     
    25512563                       + ( flux_n(k) + diss_n(k)                              &
    25522564                       -   flux_s_v(k,tn) - diss_s_v(k,tn)       ) * ddy      &
    2553                        + ( flux_t(k) + diss_t(k)                              &
    2554                        -   flux_d    - diss_d                    ) * ddzw(k)  &
     2565                       + ( ( flux_t(k) + diss_t(k) )                          &
     2566                       -   ( flux_d    + diss_d    )                          &
     2567                                                   ) * drho_air(k) * ddzw(k)  &
    25552568                                      ) + v(k,j,i) * div
    25562569
     
    26222635
    26232636          w_comp    = w(k,j-1,i) + w(k,j,i)
    2624           flux_t(k) = w_comp  * (                                             &
     2637          flux_t(k) = w_comp * rho_air_zw(k) * (                              &
    26252638                     ( 37.0_wp * ibit26 * adv_mom_5                           &
    26262639                  +     7.0_wp * ibit25 * adv_mom_3                           &
     
    26372650                                 )
    26382651
    2639           diss_t(k) = - ABS( w_comp ) * (                                     &
     2652          diss_t(k) = - ABS( w_comp ) * rho_air_zw(k) * (                     &
    26402653                     ( 10.0_wp * ibit26 * adv_mom_5                           &
    26412654                  +     3.0_wp * ibit25 * adv_mom_3                           &
     
    26562669!--       by a not sufficient reduction of divergences near topography.
    26572670          div = ( ( u_comp + gu - ( u(k,j-1,i)   + u(k,j,i)   ) ) * ddx       &
     2671                                                                  * rho_air(k)&
    26582672               +  ( v_comp(k)   - ( v(k,j,i)     + v(k,j-1,i) ) ) * ddy       &
    2659                +  ( w_comp      - ( w(k-1,j-1,i) + w(k-1,j,i) ) ) * ddzw(k)   &
     2673                                                                  * rho_air(k)&
     2674               +  (   w_comp                      * rho_air_zw(k)   -         &
     2675                    ( w(k-1,j-1,i) + w(k-1,j,i) ) * rho_air_zw(k-1)           &
     2676                  ) * ddzw(k)                                                 &
    26602677                ) * 0.5_wp
    26612678
     
    26652682                       + ( flux_n(k) + diss_n(k)                              &
    26662683                       -   flux_s_v(k,tn) - diss_s_v(k,tn)       ) * ddy      &
    2667                        + ( flux_t(k) + diss_t(k)                              &
    2668                        -   flux_d    - diss_d                    ) * ddzw(k)  &
     2684                       + ( ( flux_t(k) + diss_t(k) )                          &
     2685                       -   ( flux_d    + diss_d    )                          &
     2686                                                   ) * drho_air(k) * ddzw(k)  &
    26692687                                      ) + v(k,j,i) * div
    26702688
     
    27142732
    27152733       USE arrays_3d,                                                         &
    2716            ONLY:  ddzu, diss_l_w, diss_s_w, flux_l_w, flux_s_w, tend, u, v, w
     2734           ONLY:  ddzu, diss_l_w, diss_s_w, flux_l_w, flux_s_w, tend, u, v, w,&
     2735                  drho_air_zw, rho_air, rho_air_zw
    27172736
    27182737       USE constants,                                                         &
     
    29843003
    29853004          w_comp    = w(k+1,j,i) + w(k,j,i)
    2986           flux_t(k) = w_comp  * (                                             &
     3005          flux_t(k) = w_comp * rho_air(k+1) * (                               &
    29873006                     ( 37.0_wp * ibit35 * adv_mom_5                           &
    29883007                  +     7.0_wp * ibit34 * adv_mom_3                           &
     
    29993018                                )
    30003019
    3001           diss_t(k) = - ABS( w_comp ) * (                                     &
     3020          diss_t(k) = - ABS( w_comp ) * rho_air(k+1) * (                      &
    30023021                     ( 10.0_wp * ibit35 * adv_mom_5                           &
    30033022                  +     3.0_wp * ibit34 * adv_mom_3                           &
     
    30243043                                      + IBITS(wall_flags_0(k,j,i-1),29,1)     &
    30253044                                      )                                       &
    3026                   ) * ddx                                                     &
     3045                  ) * rho_air_zw(k) * ddx                                     &
    30273046              +   ( ( v_comp + gv ) * ( ibit30 + ibit31 + ibit32 )            &
    30283047                  - ( v(k+1,j,i) + v(k,j,i)   )                               &
     
    30313050                                      + IBITS(wall_flags_00(k,j-1,i),0,1)     &
    30323051                                      )                                       &
    3033                   ) * ddy                                                     &
    3034               +   ( w_comp          * ( ibit33 + ibit34 + ibit35 )            &
    3035                 - ( w(k,j,i)   + w(k-1,j,i)   )                               &
     3052                  ) * rho_air_zw(k) * ddy                                     &
     3053              +   ( w_comp * rho_air(k+1) * ( ibit33 + ibit34 + ibit35 )      &
     3054                - ( w(k,j,i)   + w(k-1,j,i)   ) * rho_air(k)                  &
    30363055                                    * ( IBITS(wall_flags_00(k-1,j,i),1,1)     &
    30373056                                      + IBITS(wall_flags_00(k-1,j,i),2,1)     &
     
    30473066                    + ( flux_n(k) + diss_n(k)                                 &
    30483067                    -   flux_s_w(k,tn) - diss_s_w(k,tn)       ) * ddy         &
    3049                     + ( flux_t(k) + diss_t(k)                                 &
    3050                     -   flux_d    - diss_d                    ) * ddzu(k+1)   &
     3068                    + ( ( flux_t(k) + diss_t(k) )                             &
     3069                    -   ( flux_d    + diss_d    )                             &
     3070                                              ) * drho_air_zw(k) * ddzu(k+1)  &
    30513071                                      ) + div * w(k,j,i)
    30523072
     
    31053125
    31063126          w_comp    = w(k+1,j,i) + w(k,j,i)
    3107           flux_t(k) = w_comp  * (                                             &
     3127          flux_t(k) = w_comp * rho_air(k+1) * (                               &
    31083128                     ( 37.0_wp * ibit35 * adv_mom_5                           &
    31093129                  +     7.0_wp * ibit34 * adv_mom_3                           &
     
    31203140                                )
    31213141
    3122           diss_t(k) = - ABS( w_comp ) * (                                     &
     3142          diss_t(k) = - ABS( w_comp ) * rho_air(k+1) * (                      &
    31233143                     ( 10.0_wp * ibit35 * adv_mom_5                           &
    31243144                  +     3.0_wp * ibit34 * adv_mom_3                           &
     
    31393159!--       by a not sufficient reduction of divergences near topography.
    31403160          div = ( ( u_comp + gu - ( u(k+1,j,i) + u(k,j,i)   ) ) * ddx         &
     3161                                                              * rho_air_zw(k) &
    31413162              +   ( v_comp + gv - ( v(k+1,j,i) + v(k,j,i)   ) ) * ddy         &
    3142               +   ( w_comp      - ( w(k,j,i)   + w(k-1,j,i) ) ) * ddzu(k+1)   &
     3163                                                              * rho_air_zw(k) &
     3164              +   (   w_comp                    * rho_air(k+1) -              &
     3165                    ( w(k,j,i)   + w(k-1,j,i) ) * rho_air(k)                  &
     3166                  ) * ddzu(k+1)                                               &
    31433167                ) * 0.5_wp
    31443168
     
    31483172                    + ( flux_n(k) + diss_n(k)                                 &
    31493173                    -   flux_s_w(k,tn) - diss_s_w(k,tn)       ) * ddy         &
    3150                     + ( flux_t(k) + diss_t(k)                                 &
    3151                     -   flux_d    - diss_d                    ) * ddzu(k+1)   &
     3174                    + ( ( flux_t(k) + diss_t(k) )                             &
     3175                    -   ( flux_d    + diss_d    )                             &
     3176                                              ) * drho_air_zw(k) * ddzu(k+1)  &
    31523177                                      ) + div * w(k,j,i)
    31533178
     
    31833208
    31843209       USE arrays_3d,                                                         &
    3185            ONLY:  ddzw, tend, u, v, w
     3210           ONLY:  ddzw, drho_air, tend, u, v, w, rho_air, rho_air_zw
    31863211
    31873212       USE constants,                                                         &
     
    34843509
    34853510
    3486                 flux_t(k) = w(k,j,i) * (                                      &
     3511                flux_t(k) = w(k,j,i) * rho_air_zw(k) * (                      &
    34873512                           ( 37.0_wp * ibit8 * adv_sca_5                      &
    34883513                        +     7.0_wp * ibit7 * adv_sca_3                      &
     
    34983523                                       )
    34993524
    3500                 diss_t(k) = -ABS( w(k,j,i) ) * (                              &
     3525                diss_t(k) = -ABS( w(k,j,i) ) * rho_air_zw(k) * (              &
    35013526                           ( 10.0_wp * ibit8 * adv_sca_5                      &
    35023527                        +     3.0_wp * ibit7 * adv_sca_3                      &
     
    35393564                   fd_1   = (  w(k-1,j,i)   * ( sk(k,j,i) + sk(k-1,j,i) )     &
    35403565                         -ABS( w(k-1,j,i) ) * ( sk(k,j,i) - sk(k-1,j,i) )     &
    3541                             ) * adv_sca_1
     3566                            ) * adv_sca_1 * rho_air_zw(k)
    35423567
    35433568                   ft_1   = (  w(k,j,i)   * ( sk(k+1,j,i) + sk(k,j,i) )       &
    35443569                         -ABS( w(k,j,i) ) * ( sk(k+1,j,i) - sk(k,j,i) )       &
    3545                             ) * adv_sca_1
     3570                            ) * adv_sca_1 * rho_air_zw(k)
    35463571!
    35473572!--                Calculate ratio of upwind gradients. Note, Min/Max is just
     
    36593684                                         + IBITS(wall_flags_0(k,j,i-1),2,1)    &
    36603685                                         )                                     &
    3661                           ) * ddx                                              &
     3686                          ) * rho_air(k) * ddx                                 &
    36623687                        + ( v(k,j+1,i) * ( ibit3 + ibit4 + ibit5 )             &
    36633688                          - v(k,j,i)   * ( IBITS(wall_flags_0(k,j-1,i),3,1)    &
     
    36653690                                         + IBITS(wall_flags_0(k,j-1,i),5,1)    &
    36663691                                         )                                     &
    3667                           ) * ddy                                              &
    3668                         + ( w(k,j,i)   * ( ibit6 + ibit7 + ibit8 )             &
    3669                           - w(k-1,j,i) * ( IBITS(wall_flags_0(k-1,j,i),6,1)    &
     3692                          ) * rho_air(k) * ddy                                 &
     3693                        + ( w(k,j,i) * rho_air_zw(k) *                         &
     3694                                         ( ibit6 + ibit7 + ibit8 )             &
     3695                          - w(k-1,j,i) * rho_air_zw(k-1) *                     &
     3696                                         ( IBITS(wall_flags_0(k-1,j,i),6,1)    &
    36703697                                         + IBITS(wall_flags_0(k-1,j,i),7,1)    &
    36713698                                         + IBITS(wall_flags_0(k-1,j,i),8,1)    &
     
    36793706                      + ( flux_n(k) + diss_n(k) - swap_flux_y_local(k)   -    &
    36803707                          swap_diss_y_local(k)              ) * ddy           &
    3681                       + ( flux_t(k) + diss_t(k) - flux_d - diss_d             &
    3682                                                                ) * ddzw(k)    &
     3708                      + ( ( flux_t(k) + diss_t(k) ) -                         &
     3709                          ( flux_d    + diss_d    )                           &
     3710                                                    ) * drho_air(k) * ddzw(k) &
    36833711                                            ) + sk(k,j,i) * div
    36843712
     
    37253753
    37263754
    3727                 flux_t(k) = w(k,j,i) * (                                      &
     3755                flux_t(k) = w(k,j,i) * rho_air_zw(k) * (                      &
    37283756                           ( 37.0_wp * ibit8 * adv_sca_5                      &
    37293757                        +     7.0_wp * ibit7 * adv_sca_3                      &
     
    37393767                                       )
    37403768
    3741                 diss_t(k) = -ABS( w(k,j,i) ) * (                              &
     3769                diss_t(k) = -ABS( w(k,j,i) ) * rho_air_zw(k) * (              &
    37423770                           ( 10.0_wp * ibit8 * adv_sca_5                      &
    37433771                        +     3.0_wp * ibit7 * adv_sca_3                      &
     
    37803808                   fd_1   = (  w(k-1,j,i)   * ( sk(k,j,i) + sk(k-1,j,i) )     &
    37813809                         -ABS( w(k-1,j,i) ) * ( sk(k,j,i) - sk(k-1,j,i) )     &
    3782                             ) * adv_sca_1
     3810                            ) * adv_sca_1 * rho_air_zw(k)
    37833811
    37843812                   ft_1   = (  w(k,j,i)   * ( sk(k+1,j,i) + sk(k,j,i) )       &
    37853813                         -ABS( w(k,j,i) ) * ( sk(k+1,j,i) - sk(k,j,i) )       &
    3786                             ) * adv_sca_1
     3814                            ) * adv_sca_1 * rho_air_zw(k)
    37873815!
    37883816!--                Calculate ratio of upwind gradients. Note, Min/Max is just
     
    38953923!--             correction is needed to overcome numerical instabilities introduced
    38963924!--             by a not sufficient reduction of divergences near topography.
    3897                 div         =   ( u(k,j,i+1) - u(k,j,i)   ) * ddx             &
    3898                               + ( v(k,j+1,i) - v(k,j,i)   ) * ddy             &
    3899                               + ( w(k,j,i)   - w(k-1,j,i) ) * ddzw(k)
     3925                div         =   ( u(k,j,i+1) - u(k,j,i)   ) * rho_air(k) * ddx &
     3926                              + ( v(k,j+1,i) - v(k,j,i)   ) * rho_air(k) * ddy &
     3927                              + ( w(k,j,i)   * rho_air_zw(k) -                 &
     3928                                  w(k-1,j,i) * rho_air_zw(k-1) ) * ddzw(k)
    39003929
    39013930                tend(k,j,i) = tend(k,j,i) - (                                 &
     
    39043933                      + ( flux_n(k) + diss_n(k) - swap_flux_y_local(k)   -    &
    39053934                          swap_diss_y_local(k)              ) * ddy           &
    3906                       + ( flux_t(k) + diss_t(k) - flux_d - diss_d             &
    3907                                                                ) * ddzw(k)    &
     3935                      + ( ( flux_t(k) + diss_t(k) ) -                         &
     3936                          ( flux_d    + diss_d    )                           &
     3937                                                    ) * drho_air(k) * ddzw(k) &
    39083938                                            ) + sk(k,j,i) * div
    39093939
     
    40044034
    40054035       USE arrays_3d,                                                         &
    4006            ONLY:  ddzw, tend, u, v, w
     4036           ONLY:  ddzw, drho_air, tend, u, v, w, rho_air, rho_air_zw
    40074037
    40084038       USE constants,                                                         &
     
    42784308                k_mm  = k - 2 * ibit8
    42794309
    4280                 flux_t    = w(k,j,i) * (                                      &
     4310                flux_t    = w(k,j,i) * rho_air_zw(k) * (                      &
    42814311                           ( 37.0_wp * ibit8 * adv_sca_5                      &
    42824312                        +     7.0_wp * ibit7 * adv_sca_3                      &
     
    42924322                                       )
    42934323
    4294                 diss_t    = -ABS( w(k,j,i) ) * (                              &
     4324                diss_t    = -ABS( w(k,j,i) ) * rho_air_zw(k) * (              &
    42954325                           ( 10.0_wp * ibit8 * adv_sca_5                      &
    42964326                        +     3.0_wp * ibit7 * adv_sca_3                      &
     
    43334363                   fd_1   = (  w(k-1,j,i)   * ( sk(k,j,i) + sk(k-1,j,i) )     &
    43344364                        -ABS( w(k-1,j,i) ) * ( sk(k,j,i) - sk(k-1,j,i) )      &
    4335                             ) * adv_sca_1
     4365                            ) * adv_sca_1 * rho_air_zw(k)
    43364366
    43374367                   ft_1   = (  w(k,j,i)   * ( sk(k+1,j,i) + sk(k,j,i) )       &
    43384368                        -ABS( w(k,j,i) ) * ( sk(k+1,j,i) - sk(k,j,i) )        &
    4339                             ) * adv_sca_1
     4369                            ) * adv_sca_1 * rho_air_zw(k)
    43404370!
    43414371!--                Calculate ratio of upwind gradients. Note, Min/Max is just
     
    44474477                                         + IBITS(wall_flags_0(k,j,i-1),2,1)    &
    44484478                                         )                                     &
    4449                           ) * ddx                                              &
     4479                          ) * rho_air(k) * ddx                                 &
    44504480                        + ( v(k,j+1,i) * ( ibit3 + ibit4 + ibit5 )             &
    44514481                          - v(k,j,i)   * ( IBITS(wall_flags_0(k,j-1,i),3,1)    &
     
    44534483                                         + IBITS(wall_flags_0(k,j-1,i),5,1)    &
    44544484                                         )                                     &
    4455                           ) * ddy                                              &
    4456                         + ( w(k,j,i)   * ( ibit6 + ibit7 + ibit8 )             &
    4457                           - w(k-1,j,i) * ( IBITS(wall_flags_0(k-1,j,i),6,1)    &
     4485                          ) * rho_air(k) * ddy                                 &
     4486                        + ( w(k,j,i) * rho_air_zw(k) *                         &
     4487                                         ( ibit6 + ibit7 + ibit8 )             &
     4488                          - w(k-1,j,i) * rho_air_zw(k-1) *                     &
     4489                                         ( IBITS(wall_flags_0(k-1,j,i),6,1)    &
    44584490                                         + IBITS(wall_flags_0(k-1,j,i),7,1)    &
    44594491                                         + IBITS(wall_flags_0(k-1,j,i),8,1)    &
     
    44654497                               ( flux_r + diss_r - flux_l - diss_l ) * ddx    &
    44664498                             + ( flux_n + diss_n - flux_s - diss_s ) * ddy    &
    4467                              + ( flux_t + diss_t - flux_d - diss_d ) * ddzw(k)&
     4499                             + ( ( flux_t + diss_t ) -                        &
     4500                                 ( flux_d + diss_d )                          &
     4501                                                    ) * drho_air(k) * ddzw(k) &
    44684502                                ) + div * sk(k,j,i)
    44694503
     
    45114545
    45124546       USE arrays_3d,                                                          &
    4513            ONLY:  ddzw, tend, u, v, w
     4547           ONLY:  ddzw, drho_air, tend, u, v, w, rho_air, rho_air_zw
    45144548
    45154549       USE constants,                                                          &
     
    47804814
    47814815                w_comp    = w(k,j,i) + w(k,j,i-1)
    4782                 flux_t(k) = w_comp  * (                                      &
     4816                flux_t(k) = w_comp * rho_air_zw(k) * (                       &
    47834817                          ( 37.0_wp * ibit17 * adv_mom_5                        &
    47844818                       +     7.0_wp * ibit16 * adv_mom_3                        &
     
    47954829                                      )
    47964830
    4797                 diss_t(k) = - ABS( w_comp ) * (                              &
     4831                diss_t(k) = - ABS( w_comp ) * rho_air_zw(k) * (              &
    47984832                          ( 10.0_wp * ibit17 * adv_mom_5                        &
    47994833                       +     3.0_wp * ibit16 * adv_mom_3                        &
     
    48184852                                      + IBITS(wall_flags_0(k,j,i-1),10,1)     &
    48194853                                      + IBITS(wall_flags_0(k,j,i-1),11,1)     &
    4820                                       )                                       &   
    4821                   ) * ddx                                                     &
     4854                                      )                                       &
     4855                  ) * rho_air(k) * ddx                                        &
    48224856               +  ( ( v_comp + gv ) * ( ibit12 + ibit13 + ibit14 )            &
    48234857                  - ( v(k,j,i)   + v(k,j,i-1 )  )                             &
     
    48264860                                      + IBITS(wall_flags_0(k,j-1,i),14,1)     &
    48274861                                      )                                       &
    4828                   ) * ddy                                                     &
    4829                +  ( w_comp          * ( ibit15 + ibit16 + ibit17 )            &
    4830                 - ( w(k-1,j,i) + w(k-1,j,i-1) )                               &
     4862                  ) * rho_air(k) * ddy                                        &
     4863               +  ( w_comp * rho_air_zw(k) * ( ibit15 + ibit16 + ibit17 )     &
     4864                - ( w(k-1,j,i) + w(k-1,j,i-1) ) * rho_air_zw(k-1)             &
    48314865                                    * ( IBITS(wall_flags_0(k-1,j,i),15,1)     &
    48324866                                      + IBITS(wall_flags_0(k-1,j,i),16,1)     &
     
    48434877               + ( flux_n(k) + diss_n(k)                                       &
    48444878               -   swap_flux_y_local_u(k)   - swap_diss_y_local_u(k)   ) * ddy &
    4845                + ( flux_t(k) + diss_t(k)                                       &
    4846                -   flux_d    - diss_d                                          &
    4847                                                                   ) * ddzw(k)  &
     4879               + ( ( flux_t(k) + diss_t(k) )                                   &
     4880               -   ( flux_d    + diss_d    )                                   &
     4881                                                    ) * drho_air(k) * ddzw(k)  &
    48484882                                           ) + div * u(k,j,i)
    48494883
     
    49114945
    49124946                w_comp    = w(k,j,i) + w(k,j,i-1)
    4913                 flux_t(k) = w_comp  * (                                      &
     4947                flux_t(k) = w_comp * rho_air_zw(k) * (                       &
    49144948                          ( 37.0_wp * ibit17 * adv_mom_5                        &
    49154949                       +     7.0_wp * ibit16 * adv_mom_3                        &
     
    49264960                                      )
    49274961
    4928                 diss_t(k) = - ABS( w_comp ) * (                              &
     4962                diss_t(k) = - ABS( w_comp ) * rho_air_zw(k) * (              &
    49294963                          ( 10.0_wp * ibit17 * adv_mom_5                        &
    49304964                       +     3.0_wp * ibit16 * adv_mom_3                        &
     
    49454979!--             by a not sufficient reduction of divergences near topography.
    49464980                div = ( ( u_comp(k)   - ( u(k,j,i)   + u(k,j,i-1)   ) ) * ddx &
     4981                                                                  * rho_air(k)&
    49474982                     +  ( v_comp + gv - ( v(k,j,i)   + v(k,j,i-1 )  ) ) * ddy &
    4948                      +  ( w_comp      - ( w(k-1,j,i) + w(k-1,j,i-1) ) )       &
    4949                                                                     * ddzw(k) &
     4983                                                                  * rho_air(k)&
     4984                     +  (   w_comp                      * rho_air_zw(k) -     &
     4985                          ( w(k-1,j,i) + w(k-1,j,i-1) ) * rho_air_zw(k-1)     &
     4986                        ) * ddzw(k)                                           &
    49504987                      ) * 0.5_wp
    49514988
     
    49554992               + ( flux_n(k) + diss_n(k)                                       &
    49564993               -   swap_flux_y_local_u(k)   - swap_diss_y_local_u(k)   ) * ddy &
    4957                + ( flux_t(k) + diss_t(k)                                       &
    4958                -   flux_d    - diss_d                                          &
    4959                                                                   ) * ddzw(k)  &
     4994               + ( ( flux_t(k) + diss_t(k) )                                   &
     4995               -   ( flux_d    + diss_d    )                                   &
     4996                                                    ) * drho_air(k) * ddzw(k)  &
    49604997                                           ) + div * u(k,j,i)
    49614998
     
    50045041
    50055042       USE arrays_3d,                                                          &
    5006            ONLY:  ddzw, tend, u, v, w
     5043           ONLY:  ddzw, drho_air, tend, u, v, w, rho_air, rho_air_zw
    50075044
    50085045       USE constants,                                                          &
     
    52675304
    52685305                w_comp    = w(k,j,i) + w(k,j,i-1)
    5269                 flux_t    = w_comp  * (                                      &
     5306                flux_t    = w_comp * rho_air_zw(k) * (                       &
    52705307                          ( 37.0_wp * ibit17 * adv_mom_5                        &
    52715308                       +     7.0_wp * ibit16 * adv_mom_3                        &
     
    52825319                                      )
    52835320
    5284                 diss_t    = - ABS( w_comp ) * (                              &
     5321                diss_t    = - ABS( w_comp ) * rho_air_zw(k) * (              &
    52855322                          ( 10.0_wp * ibit17 * adv_mom_5                        &
    52865323                       +     3.0_wp * ibit16 * adv_mom_3                        &
     
    53055342                                      + IBITS(wall_flags_0(k,j,i-1),10,1)     &
    53065343                                      + IBITS(wall_flags_0(k,j,i-1),11,1)     &
    5307                                       )                                       &   
    5308                   ) * ddx                                                     &
     5344                                      )                                       &
     5345                  ) * rho_air(k) * ddx                                        &
    53095346               +  ( ( v_comp + gv ) * ( ibit12 + ibit13 + ibit14 )            &
    53105347                  - ( v(k,j,i)   + v(k,j,i-1 )  )                             &
     
    53135350                                      + IBITS(wall_flags_0(k,j-1,i),14,1)     &
    53145351                                      )                                       &
    5315                   ) * ddy                                                     &
    5316                +  ( w_comp          * ( ibit15 + ibit16 + ibit17 )            &
    5317                 - ( w(k-1,j,i) + w(k-1,j,i-1) )                               &
     5352                  ) * rho_air(k) * ddy                                        &
     5353               +  ( w_comp * rho_air_zw(k) * ( ibit15 + ibit16 + ibit17 )     &
     5354                - ( w(k-1,j,i) + w(k-1,j,i-1) ) * rho_air_zw(k-1)             &
    53185355                                    * ( IBITS(wall_flags_0(k-1,j,i),15,1)     &
    53195356                                      + IBITS(wall_flags_0(k-1,j,i),16,1)     &
    53205357                                      + IBITS(wall_flags_0(k-1,j,i),17,1)     &
    5321                                       )                                       & 
     5358                                      )                                       &
    53225359                  ) * ddzw(k)   &
    53235360                ) * 0.5_wp
     
    53275364                               ( flux_r + diss_r - flux_l - diss_l ) * ddx     &
    53285365                             + ( flux_n + diss_n - flux_s - diss_s ) * ddy     &
    5329                              + ( flux_t + diss_t - flux_d - diss_d ) * ddzw(k) &
     5366                             + ( ( flux_t + diss_t ) -                         &
     5367                                 ( flux_d + diss_d )                           &
     5368                                                     ) * drho_air(k) * ddzw(k) &
    53305369                                ) + div * u(k,j,i)
    53315370
     
    53655404
    53665405       USE arrays_3d,                                                          &
    5367            ONLY:  ddzw, tend, u, v, w
     5406           ONLY:  ddzw, drho_air, tend, u, v, w, rho_air, rho_air_zw
    53685407
    53695408       USE constants,                                                          &
     
    56335672
    56345673                w_comp    = w(k,j-1,i) + w(k,j,i)
    5635                 flux_t(k) = w_comp  * (                                      &
     5674                flux_t(k) = w_comp * rho_air_zw(k) * (                       &
    56365675                          ( 37.0_wp * ibit26 * adv_mom_5                        &
    56375676                       +     7.0_wp * ibit25 * adv_mom_3                        &
     
    56485687                                      )
    56495688
    5650                 diss_t(k) = - ABS( w_comp ) * (                              &
     5689                diss_t(k) = - ABS( w_comp ) * rho_air_zw(k) * (              &
    56515690                          ( 10.0_wp * ibit26 * adv_mom_5                        &
    56525691                       +     3.0_wp * ibit25 * adv_mom_3                        &
     
    56725711                                         + IBITS(wall_flags_0(k,j,i-1),19,1)  &
    56735712                                         + IBITS(wall_flags_0(k,j,i-1),20,1)  &
    5674                                          )                                    &   
    5675                   ) * ddx                                                     &
     5713                                         )                                    &
     5714                  ) * rho_air(k) * ddx                                        &
    56765715               +  ( v_comp(k)                                                 &
    56775716                                       * ( ibit21 + ibit22 + ibit23 )         &
     
    56805719                                         + IBITS(wall_flags_0(k,j-1,i),22,1)  &
    56815720                                         + IBITS(wall_flags_0(k,j-1,i),23,1)  &
    5682                                          )                                    &   
    5683                   ) * ddy                                                     &
    5684                +  ( w_comp                                                    &
     5721                                         )                                    &
     5722                  ) * rho_air(k) * ddy                                        &
     5723               +  ( w_comp * rho_air_zw(k)                                    &
    56855724                                       * ( ibit24 + ibit25 + ibit26 )         &
    5686                 - ( w(k-1,j-1,i) + w(k-1,j,i) )                               &
     5725                - ( w(k-1,j-1,i) + w(k-1,j,i) ) * rho_air_zw(k-1)             &
    56875726                                       * ( IBITS(wall_flags_0(k-1,j,i),24,1)  &
    56885727                                         + IBITS(wall_flags_0(k-1,j,i),25,1)  &
     
    57005739                     -   swap_flux_y_local_v(k) - swap_diss_y_local_v(k)      &
    57015740                       ) * ddy                                                &
    5702                      + ( flux_t(k) + diss_t(k)                                &
    5703                      -   flux_d    - diss_d                                   &
    5704                        ) * ddzw(k)                                            &
     5741                     + ( ( flux_t(k) + diss_t(k) )                            &
     5742                     -   ( flux_d    + diss_d    )                            &
     5743                       ) * drho_air(k) * ddzw(k)                              &
    57055744                                            )  + v(k,j,i) * div
    57065745
     
    57725811
    57735812                w_comp    = w(k,j-1,i) + w(k,j,i)
    5774                 flux_t(k) = w_comp  * (                                      &
     5813                flux_t(k) = w_comp * rho_air_zw(k) * (                       &
    57755814                          ( 37.0_wp * ibit26 * adv_mom_5                        &
    57765815                       +     7.0_wp * ibit25 * adv_mom_3                        &
     
    57875826                                      )
    57885827
    5789                 diss_t(k) = - ABS( w_comp ) * (                              &
     5828                diss_t(k) = - ABS( w_comp ) * rho_air_zw(k) * (              &
    57905829                          ( 10.0_wp * ibit26 * adv_mom_5                        &
    57915830                       +     3.0_wp * ibit25 * adv_mom_3                        &
     
    58065845!--             by a not sufficient reduction of divergences near topography.
    58075846                div = ( ( u_comp + gu - ( u(k,j-1,i)   + u(k,j,i)   ) ) * ddx &
     5847                                                                  * rho_air(k)&
    58085848                     +  ( v_comp(k)   - ( v(k,j,i)     + v(k,j-1,i) ) ) * ddy &
    5809                      +  ( w_comp      - ( w(k-1,j-1,i) + w(k-1,j,i) ) )       &
    5810                                                                     * ddzw(k) &
     5849                                                                  * rho_air(k)&
     5850                     +  (   w_comp                      * rho_air_zw(k) -     &
     5851                          ( w(k-1,j-1,i) + w(k-1,j,i) ) * rho_air_zw(k-1)     &
     5852                        ) * ddzw(k)                                           &
    58115853                      ) * 0.5_wp
    58125854 
     
    58185860                     -   swap_flux_y_local_v(k) - swap_diss_y_local_v(k)      &
    58195861                       ) * ddy                                                &
    5820                      + ( flux_t(k) + diss_t(k)                                &
    5821                      -   flux_d    - diss_d                                   &
    5822                        ) * ddzw(k)                                            &
     5862                     + ( ( flux_t(k) + diss_t(k) )                            &
     5863                     -   ( flux_d    + diss_d    )                            &
     5864                       ) * drho_air(k) * ddzw(k)                              &
    58235865                                            )  + v(k,j,i) * div
    58245866
     
    58695911
    58705912       USE arrays_3d,                                                          &
    5871            ONLY:  ddzw, tend, u, v, w
     5913           ONLY:  ddzw, drho_air, tend, u, v, w, rho_air, rho_air_zw
    58725914
    58735915       USE constants,                                                          &
     
    61326174
    61336175                w_comp    = w(k,j-1,i) + w(k,j,i)
    6134                 flux_t    = w_comp  * (                                      &
     6176                flux_t    = w_comp * rho_air_zw(k) * (                       &
    61356177                          ( 37.0_wp * ibit26 * adv_mom_5                        &
    61366178                       +     7.0_wp * ibit25 * adv_mom_3                        &
     
    61476189                                      )
    61486190
    6149                 diss_t    = - ABS( w_comp ) * (                              &
     6191                diss_t    = - ABS( w_comp ) * rho_air_zw(k) * (              &
    61506192                          ( 10.0_wp * ibit26 * adv_mom_5                        &
    61516193                       +     3.0_wp * ibit25 * adv_mom_3                        &
     
    61716213                                         + IBITS(wall_flags_0(k,j,i-1),19,1)  &
    61726214                                         + IBITS(wall_flags_0(k,j,i-1),20,1)  &
    6173                                          )                                    &   
    6174                   ) * ddx                                                     &
     6215                                         )                                    &
     6216                  ) * rho_air(k) * ddx                                        &
    61756217               +  ( v_comp                                                    &
    61766218                                       * ( ibit21 + ibit22 + ibit23 )         &
     
    61796221                                         + IBITS(wall_flags_0(k,j-1,i),22,1)  &
    61806222                                         + IBITS(wall_flags_0(k,j-1,i),23,1)  &
    6181                                          )                                    &   
    6182                   ) * ddy                                                     &
    6183                +  ( w_comp                                                    &
     6223                                         )                                    &
     6224                  ) * rho_air(k) * ddy                                        &
     6225               +  ( w_comp * rho_air_zw(k)                                    &
    61846226                                       * ( ibit24 + ibit25 + ibit26 )         &
    6185                 - ( w(k-1,j-1,i) + w(k-1,j,i) )                               &
     6227                - ( w(k-1,j-1,i) + w(k-1,j,i) ) * rho_air_zw(k-1)             &
    61866228                                       * ( IBITS(wall_flags_0(k-1,j,i),24,1)  &
    61876229                                         + IBITS(wall_flags_0(k-1,j,i),25,1)  &
     
    61956237                               ( flux_r + diss_r - flux_l - diss_l ) * ddx     &
    61966238                             + ( flux_n + diss_n - flux_s - diss_s ) * ddy     &
    6197                              + ( flux_t + diss_t - flux_d - diss_d ) * ddzw(k) &
     6239                             + ( ( flux_t + diss_t ) -                         &
     6240                                 ( flux_d + diss_d )                           &
     6241                               ) * drho_air(k) * ddzw(k)                       &
    61986242                                ) + div * v(k,j,i)
    61996243
     
    62356279
    62366280       USE arrays_3d,                                                          &
    6237            ONLY:  ddzu, tend, u, v, w
     6281           ONLY:  ddzu, drho_air_zw, tend, u, v, w, rho_air, rho_air_zw
    62386282
    62396283       USE constants,                                                          &
     
    65096553
    65106554                w_comp    = w(k+1,j,i) + w(k,j,i)
    6511                 flux_t(k) = w_comp  * (                                      &
     6555                flux_t(k) = w_comp * rho_air(k+1) * (                        &
    65126556                          ( 37.0_wp * ibit35 * adv_mom_5                        &
    65136557                       +     7.0_wp * ibit34 * adv_mom_3                        &
     
    65246568                                       )
    65256569
    6526                 diss_t(k) = - ABS( w_comp ) * (                              &
     6570                diss_t(k) = - ABS( w_comp ) * rho_air(k+1) * (               &
    65276571                          ( 10.0_wp * ibit35 * adv_mom_5                        &
    65286572                       +     3.0_wp * ibit34 * adv_mom_3                        &
     
    65486592                                      + IBITS(wall_flags_0(k,j,i-1),29,1)     &
    65496593                                      )                                       &
    6550                   ) * ddx                                                     &
    6551               +   ( ( v_comp + gv ) * ( ibit30 + ibit31 + ibit32 )            & 
     6594                  ) * rho_air_zw(k) * ddx                                     &
     6595              +   ( ( v_comp + gv ) * ( ibit30 + ibit31 + ibit32 )            &
    65526596                  - ( v(k+1,j,i) + v(k,j,i)   )                               &
    65536597                                    * ( IBITS(wall_flags_0(k,j-1,i),30,1)     &
     
    65556599                                      + IBITS(wall_flags_00(k,j-1,i),0,1)     &
    65566600                                      )                                       &
    6557                   ) * ddy                                                     &
    6558               +   ( w_comp          * ( ibit33 + ibit34 + ibit35 )            &
    6559                 - ( w(k,j,i)   + w(k-1,j,i)   )                               &
     6601                  ) * rho_air_zw(k) * ddy                                     &
     6602              +   ( w_comp * rho_air(k+1) * ( ibit33 + ibit34 + ibit35 )      &
     6603                - ( w(k,j,i)   + w(k-1,j,i)   ) * rho_air(k)                  &
    65606604                                    * ( IBITS(wall_flags_00(k-1,j,i),1,1)     &
    65616605                                      + IBITS(wall_flags_00(k-1,j,i),2,1)     &
     
    65746618                    -   swap_flux_y_local_w(k)   - swap_diss_y_local_w(k)     &
    65756619                      ) * ddy                                                 &
    6576                     + ( flux_t(k) + diss_t(k)                                 &
    6577                     -   flux_d    - diss_d                                    &
    6578                       ) * ddzu(k+1)                                           &
     6620                    + ( ( flux_t(k) + diss_t(k) )                             &
     6621                    -   ( flux_d    + diss_d    )                             &
     6622                      ) * drho_air_zw(k) * ddzu(k+1)                          &
    65796623                                            )  + div * w(k,j,i)
    65806624
     
    66326676
    66336677                w_comp    = w(k+1,j,i) + w(k,j,i)
    6634                 flux_t(k) = w_comp  * (                                      &
     6678                flux_t(k) = w_comp * rho_air(k+1) * (                        &
    66356679                          ( 37.0_wp * ibit35 * adv_mom_5                        &
    66366680                       +     7.0_wp * ibit34 * adv_mom_3                        &
     
    66476691                                       )
    66486692
    6649                 diss_t(k) = - ABS( w_comp ) * (                              &
     6693                diss_t(k) = - ABS( w_comp ) * rho_air(k+1) * (               &
    66506694                          ( 10.0_wp * ibit35 * adv_mom_5                        &
    66516695                       +     3.0_wp * ibit34 * adv_mom_3                        &
     
    66666710!--             by a not sufficient reduction of divergences near topography.
    66676711                div = ( ( u_comp + gu - ( u(k+1,j,i) + u(k,j,i)   ) ) * ddx  &
     6712                                                             * rho_air_zw(k) &
    66686713                    +   ( v_comp + gv - ( v(k+1,j,i) + v(k,j,i)   ) ) * ddy  &
    6669                     +   ( w_comp      - ( w(k,j,i)   + w(k-1,j,i) ) )        &
    6670                                                                  * ddzu(k+1) &
     6714                                                             * rho_air_zw(k) &
     6715                    +   (   w_comp                    * rho_air(k+1) -       &
     6716                          ( w(k,j,i)   + w(k-1,j,i) ) * rho_air(k)           &
     6717                        ) * ddzu(k+1)                                        &
    66716718                      ) * 0.5_wp
    66726719
     
    66786725                    -   swap_flux_y_local_w(k)   - swap_diss_y_local_w(k)     &
    66796726                      ) * ddy                                                 &
    6680                     + ( flux_t(k) + diss_t(k)                                 &
    6681                     -   flux_d    - diss_d                                    &
    6682                       ) * ddzu(k+1)                                           &
     6727                    + ( ( flux_t(k) + diss_t(k) )                             &
     6728                    -   ( flux_d    + diss_d    )                             &
     6729                      ) * drho_air_zw(k) * ddzu(k+1)                          &
    66836730                                            )  + div * w(k,j,i)
    66846731
     
    67146761
    67156762       USE arrays_3d,                                                          &
    6716            ONLY:  ddzu, tend, u, v, w
     6763           ONLY:  ddzu, drho_air_zw, tend, u, v, w, rho_air, rho_air_zw
    67176764
    67186765       USE constants,                                                          &
     
    69767023
    69777024                w_comp    = w(k+1,j,i) + w(k,j,i)
    6978                 flux_t    = w_comp  * (                                      &
     7025                flux_t    = w_comp * rho_air(k+1) * (                        &
    69797026                          ( 37.0_wp * ibit35 * adv_mom_5                        &
    69807027                       +     7.0_wp * ibit34 * adv_mom_3                        &
     
    69917038                                       )
    69927039
    6993                 diss_t    = - ABS( w_comp ) * (                              &
     7040                diss_t    = - ABS( w_comp ) * rho_air(k+1) * (               &
    69947041                          ( 10.0_wp * ibit35 * adv_mom_5                        &
    69957042                       +     3.0_wp * ibit34 * adv_mom_3                        &
     
    70157062                                      + IBITS(wall_flags_0(k,j,i-1),29,1)     &
    70167063                                      )                                       &
    7017                   ) * ddx                                                     &
    7018               +   ( ( v_comp + gv ) * ( ibit30 + ibit31 + ibit32 )            & 
     7064                  ) * rho_air_zw(k) * ddx                                     &
     7065              +   ( ( v_comp + gv ) * ( ibit30 + ibit31 + ibit32 )            &
    70197066                  - ( v(k+1,j,i) + v(k,j,i)   )                               &
    70207067                                    * ( IBITS(wall_flags_0(k,j-1,i),30,1)     &
     
    70227069                                      + IBITS(wall_flags_00(k,j-1,i),0,1)     &
    70237070                                      )                                       &
    7024                   ) * ddy                                                     &
    7025               +   ( w_comp          * ( ibit33 + ibit34 + ibit35 )            &
    7026                 - ( w(k,j,i)   + w(k-1,j,i)   )                               &
     7071                  ) * rho_air_zw(k) * ddy                                     &
     7072              +   ( w_comp * rho_air(k+1) * ( ibit33 + ibit34 + ibit35 )      &
     7073                - ( w(k,j,i)   + w(k-1,j,i)   ) * rho_air(k)                  &
    70277074                                    * ( IBITS(wall_flags_00(k-1,j,i),1,1)     &
    70287075                                      + IBITS(wall_flags_00(k-1,j,i),2,1)     &
     
    70367083                               ( flux_r + diss_r - flux_l - diss_l ) * ddx       &
    70377084                             + ( flux_n + diss_n - flux_s - diss_s ) * ddy       &
    7038                              + ( flux_t + diss_t - flux_d - diss_d ) * ddzu(k+1) &
     7085                             + ( ( flux_t + diss_t ) -                           &
     7086                                 ( flux_d + diss_d ) * rho_air(k)                &
     7087                               ) * drho_air_zw(k) * ddzu(k+1)                    &
    70397088                                 ) + div * w(k,j,i)
    70407089
Note: See TracChangeset for help on using the changeset viewer.