Changeset 2037


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

Anelastic approximation implemented

Location:
palm/trunk/SOURCE
Files:
22 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
  • palm/trunk/SOURCE/check_parameters.f90

    r2032 r2037  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! Anelastic approximation implemented
    2323!
    2424! Former revisions:
     
    432432    USE netcdf_interface,                                                      &
    433433        ONLY:  dopr_unit, do2d_unit, do3d_unit, netcdf_data_format,            &
    434                netcdf_data_format_string
     434               netcdf_data_format_string, dots_unit, heatflux_output_unit,     &
     435               waterflux_output_unit, momentumflux_output_unit
    435436
    436437    USE particle_attributes
     
    837838
    838839    ENDIF
     840
     841!
     842!-- Check approximation
     843    IF ( TRIM( approximation ) /= 'boussinesq'   .AND.   &
     844         TRIM( approximation ) /= 'anelastic' )  THEN
     845       message_string = 'unknown approximation: approximation = "' //    &
     846                        TRIM( approximation ) // '"'
     847       CALL message( 'check_parameters', 'PA0000', 1, 2, 0, 6, 0 )
     848    ENDIF
     849
     850!
     851!-- Check approximation requirements
     852    IF ( TRIM( approximation ) == 'anelastic'   .AND.   &
     853         TRIM( momentum_advec ) /= 'ws-scheme' )  THEN
     854       message_string = 'Anelastic approximation requires: ' //                &
     855                        'momentum_advec = "ws-scheme"'
     856       CALL message( 'check_parameters', 'PA0000', 1, 2, 0, 6, 0 )
     857    ENDIF
     858    IF ( TRIM( approximation ) == 'anelastic'   .AND.   &
     859         TRIM( psolver ) == 'multigrid' )  THEN
     860       message_string = 'Anelastic approximation currently only supports: ' // &
     861                        'psolver = "poisfft", ' // &
     862                        'psolver = "sor" and ' // &
     863                        'psolver = "multigrid_noopt"'
     864       CALL message( 'check_parameters', 'PA0000', 1, 2, 0, 6, 0 )
     865    ENDIF
     866    IF ( TRIM( approximation ) == 'anelastic'   .AND.   &
     867         conserve_volume_flow )  THEN
     868       message_string = 'Anelastic approximation is not allowed with:' //      &
     869                        'conserve_volume_flow = .TRUE.'
     870       CALL message( 'check_parameters', 'PA0000', 1, 2, 0, 6, 0 )
     871    ENDIF
     872
     873!
     874!-- Check flux input mode
     875    IF ( TRIM( flux_input_mode ) /= 'dynamic'    .AND.   &
     876         TRIM( flux_input_mode ) /= 'kinematic'  .AND.   &
     877         TRIM( flux_input_mode ) /= 'approximation-specific' )  THEN
     878       message_string = 'unknown flux input mode: flux_input_mode = "' //      &
     879                        TRIM( flux_input_mode ) // '"'
     880       CALL message( 'check_parameters', 'PA0000', 1, 2, 0, 6, 0 )
     881    ENDIF
     882!-- Set flux input mode according to approximation if applicable
     883    IF ( TRIM( flux_input_mode ) == 'approximation-specific' )  THEN
     884       IF ( TRIM( approximation ) == 'anelastic' )  THEN
     885          flux_input_mode = 'dynamic'
     886       ELSEIF ( TRIM( approximation ) == 'boussinesq' )  THEN
     887          flux_input_mode = 'kinematic'
     888       ENDIF
     889    ENDIF
     890
     891!
     892!-- Check flux output mode
     893    IF ( TRIM( flux_output_mode ) /= 'dynamic'    .AND.   &
     894         TRIM( flux_output_mode ) /= 'kinematic'  .AND.   &
     895         TRIM( flux_output_mode ) /= 'approximation-specific' )  THEN
     896       message_string = 'unknown flux output mode: flux_output_mode = "' //    &
     897                        TRIM( flux_output_mode ) // '"'
     898       CALL message( 'check_parameters', 'PA0000', 1, 2, 0, 6, 0 )
     899    ENDIF
     900!-- Set flux output mode according to approximation if applicable
     901    IF ( TRIM( flux_output_mode ) == 'approximation-specific' )  THEN
     902       IF ( TRIM( approximation ) == 'anelastic' )  THEN
     903          flux_output_mode = 'dynamic'
     904       ELSEIF ( TRIM( approximation ) == 'boussinesq' )  THEN
     905          flux_output_mode = 'kinematic'
     906       ENDIF
     907    ENDIF
     908
     909!
     910!-- set the flux output units according to flux_output_mode
     911    IF ( TRIM( flux_output_mode ) == 'kinematic' ) THEN
     912        heatflux_output_unit              = 'K m/s'
     913        waterflux_output_unit             = 'kg/kg m/s'
     914        momentumflux_output_unit          = 'm2/s2'
     915    ELSEIF ( TRIM( flux_output_mode ) == 'dynamic' ) THEN
     916        heatflux_output_unit              = 'W/m2'
     917        waterflux_output_unit             = 'W/m2'
     918        momentumflux_output_unit          = 'N/m2'
     919    ENDIF
     920
     921!-- set time series output units for fluxes
     922    dots_unit(14:16) = heatflux_output_unit
     923    dots_unit(21)    = waterflux_output_unit
     924    dots_unit(19:20) = momentumflux_output_unit
    839925
    840926!
     
    21152201          CASE ( 'w"u"' )
    21162202             dopr_index(i) = 12
    2117              dopr_unit(i)  = 'm2/s2'
     2203             dopr_unit(i)  = TRIM ( momentumflux_output_unit )
    21182204             hom(:,2,12,:) = SPREAD( zw, 2, statistic_regions+1 )
    21192205             IF ( constant_flux_layer )  hom(nzb,2,12,:) = zu(1)
     
    21212207          CASE ( 'w*u*' )
    21222208             dopr_index(i) = 13
    2123              dopr_unit(i)  = 'm2/s2'
     2209             dopr_unit(i)  = TRIM ( momentumflux_output_unit )
    21242210             hom(:,2,13,:) = SPREAD( zw, 2, statistic_regions+1 )
    21252211
    21262212          CASE ( 'w"v"' )
    21272213             dopr_index(i) = 14
    2128              dopr_unit(i)  = 'm2/s2'
     2214             dopr_unit(i)  = TRIM ( momentumflux_output_unit )
    21292215             hom(:,2,14,:) = SPREAD( zw, 2, statistic_regions+1 )
    21302216             IF ( constant_flux_layer )  hom(nzb,2,14,:) = zu(1)
     
    21322218          CASE ( 'w*v*' )
    21332219             dopr_index(i) = 15
    2134              dopr_unit(i)  = 'm2/s2'
     2220             dopr_unit(i)  = TRIM ( momentumflux_output_unit )
    21352221             hom(:,2,15,:) = SPREAD( zw, 2, statistic_regions+1 )
    21362222
    21372223          CASE ( 'w"pt"' )
    21382224             dopr_index(i) = 16
    2139              dopr_unit(i)  = 'K m/s'
     2225             dopr_unit(i)  = TRIM ( heatflux_output_unit )
    21402226             hom(:,2,16,:) = SPREAD( zw, 2, statistic_regions+1 )
    21412227
    21422228          CASE ( 'w*pt*' )
    21432229             dopr_index(i) = 17
    2144              dopr_unit(i)  = 'K m/s'
     2230             dopr_unit(i)  = TRIM ( heatflux_output_unit )
    21452231             hom(:,2,17,:) = SPREAD( zw, 2, statistic_regions+1 )
    21462232
    21472233          CASE ( 'wpt' )
    21482234             dopr_index(i) = 18
    2149              dopr_unit(i)  = 'K m/s'
     2235             dopr_unit(i)  = TRIM ( heatflux_output_unit )
    21502236             hom(:,2,18,:) = SPREAD( zw, 2, statistic_regions+1 )
    21512237
    21522238          CASE ( 'wu' )
    21532239             dopr_index(i) = 19
    2154              dopr_unit(i)  = 'm2/s2'
     2240             dopr_unit(i)  = TRIM ( momentumflux_output_unit )
    21552241             hom(:,2,19,:) = SPREAD( zw, 2, statistic_regions+1 )
    21562242             IF ( constant_flux_layer )  hom(nzb,2,19,:) = zu(1)
     
    21582244          CASE ( 'wv' )
    21592245             dopr_index(i) = 20
    2160              dopr_unit(i)  = 'm2/s2'
     2246             dopr_unit(i)  = TRIM ( momentumflux_output_unit )
    21612247             hom(:,2,20,:) = SPREAD( zw, 2, statistic_regions+1 )
    21622248             IF ( constant_flux_layer )  hom(nzb,2,20,:) = zu(1)
     
    21642250          CASE ( 'w*pt*BC' )
    21652251             dopr_index(i) = 21
    2166              dopr_unit(i)  = 'K m/s'
     2252             dopr_unit(i)  = TRIM ( heatflux_output_unit )
    21672253             hom(:,2,21,:) = SPREAD( zw, 2, statistic_regions+1 )
    21682254
    21692255          CASE ( 'wptBC' )
    21702256             dopr_index(i) = 22
    2171              dopr_unit(i)  = 'K m/s'
     2257             dopr_unit(i)  = TRIM ( heatflux_output_unit )
    21722258             hom(:,2,22,:) = SPREAD( zw, 2, statistic_regions+1 )
    21732259
     
    23352421          CASE ( 'w"vpt"' )
    23362422             dopr_index(i) = 45
    2337              dopr_unit(i)  = 'K m/s'
     2423             dopr_unit(i)  = TRIM ( heatflux_output_unit )
    23382424             hom(:,2,45,:) = SPREAD( zw, 2, statistic_regions+1 )
    23392425
    23402426          CASE ( 'w*vpt*' )
    23412427             dopr_index(i) = 46
    2342              dopr_unit(i)  = 'K m/s'
     2428             dopr_unit(i)  = TRIM ( heatflux_output_unit )
    23432429             hom(:,2,46,:) = SPREAD( zw, 2, statistic_regions+1 )
    23442430
    23452431          CASE ( 'wvpt' )
    23462432             dopr_index(i) = 47
    2347              dopr_unit(i)  = 'K m/s'
     2433             dopr_unit(i)  = TRIM ( heatflux_output_unit )
    23482434             hom(:,2,47,:) = SPREAD( zw, 2, statistic_regions+1 )
    23492435
     
    23562442             ELSE
    23572443                dopr_index(i) = 48
    2358                 dopr_unit(i)  = 'kg/kg m/s'
     2444                dopr_unit(i)  = TRIM ( waterflux_output_unit )
    23592445                hom(:,2,48,:) = SPREAD( zw, 2, statistic_regions+1 )
    23602446             ENDIF
     
    23682454             ELSE
    23692455                dopr_index(i) = 49
    2370                 dopr_unit(i)  = 'kg/kg m/s'
     2456                dopr_unit(i)  = TRIM ( waterflux_output_unit )
    23712457                hom(:,2,49,:) = SPREAD( zw, 2, statistic_regions+1 )
    23722458             ENDIF
     
    23802466             ELSE
    23812467                dopr_index(i) = 50
    2382                 dopr_unit(i)  = 'kg/kg m/s'
     2468                dopr_unit(i)  = TRIM ( waterflux_output_unit )
    23832469                hom(:,2,50,:) = SPREAD( zw, 2, statistic_regions+1 )
    23842470             ENDIF
     
    24232509             IF ( humidity  .AND.  .NOT.  cloud_physics )  THEN
    24242510                dopr_index(i) = 48
    2425                 dopr_unit(i)  = 'kg/kg m/s'
     2511                dopr_unit(i)  = TRIM ( waterflux_output_unit )
    24262512                hom(:,2,48,:) = SPREAD( zw, 2, statistic_regions+1 )
    24272513             ELSEIF ( humidity  .AND.  cloud_physics )  THEN
    24282514                dopr_index(i) = 51
    2429                 dopr_unit(i)  = 'kg/kg m/s'
     2515                dopr_unit(i)  = TRIM ( waterflux_output_unit )
    24302516                hom(:,2,51,:) = SPREAD( zw, 2, statistic_regions+1 )
    24312517             ELSE
     
    24412527             THEN
    24422528                dopr_index(i) = 49
    2443                 dopr_unit(i)  = 'kg/kg m/s'
     2529                dopr_unit(i)  = TRIM ( waterflux_output_unit )
    24442530                hom(:,2,49,:) = SPREAD( zw, 2, statistic_regions+1 )
    24452531             ELSEIF( humidity .AND. cloud_physics ) THEN
    24462532                dopr_index(i) = 52
    2447                 dopr_unit(i)  = 'kg/kg m/s'
     2533                dopr_unit(i)  = TRIM ( waterflux_output_unit )
    24482534                hom(:,2,52,:) = SPREAD( zw, 2, statistic_regions+1 )
    24492535             ELSE
     
    24582544             IF ( humidity  .AND.  .NOT.  cloud_physics )  THEN
    24592545                dopr_index(i) = 50
    2460                 dopr_unit(i)  = 'kg/kg m/s'
     2546                dopr_unit(i)  = TRIM ( waterflux_output_unit )
    24612547                hom(:,2,50,:) = SPREAD( zw, 2, statistic_regions+1 )
    24622548             ELSEIF ( humidity  .AND.  cloud_physics )  THEN
    24632549                dopr_index(i) = 53
    2464                 dopr_unit(i)  = 'kg/kg m/s'
     2550                dopr_unit(i)  = TRIM ( waterflux_output_unit )
    24652551                hom(:,2,53,:) = SPREAD( zw, 2, statistic_regions+1 )
    24662552             ELSE
     
    25032589          CASE ( 'u"pt"' )
    25042590             dopr_index(i) = 58
    2505              dopr_unit(i)  = 'K m/s'
     2591             dopr_unit(i)  = TRIM ( heatflux_output_unit )
    25062592             hom(:,2,58,:) = SPREAD( zu, 2, statistic_regions+1 )
    25072593
    25082594          CASE ( 'u*pt*' )
    25092595             dopr_index(i) = 59
    2510              dopr_unit(i)  = 'K m/s'
     2596             dopr_unit(i)  = TRIM ( heatflux_output_unit )
    25112597             hom(:,2,59,:) = SPREAD( zu, 2, statistic_regions+1 )
    25122598
    25132599          CASE ( 'upt_t' )
    25142600             dopr_index(i) = 60
    2515              dopr_unit(i)  = 'K m/s'
     2601             dopr_unit(i)  = TRIM ( heatflux_output_unit )
    25162602             hom(:,2,60,:) = SPREAD( zu, 2, statistic_regions+1 )
    25172603
    25182604          CASE ( 'v"pt"' )
    25192605             dopr_index(i) = 61
    2520              dopr_unit(i)  = 'K m/s'
     2606             dopr_unit(i)  = TRIM ( heatflux_output_unit )
    25212607             hom(:,2,61,:) = SPREAD( zu, 2, statistic_regions+1 )
    25222608             
    25232609          CASE ( 'v*pt*' )
    25242610             dopr_index(i) = 62
    2525              dopr_unit(i)  = 'K m/s'
     2611             dopr_unit(i)  = TRIM ( heatflux_output_unit )
    25262612             hom(:,2,62,:) = SPREAD( zu, 2, statistic_regions+1 )
    25272613
    25282614          CASE ( 'vpt_t' )
    25292615             dopr_index(i) = 63
    2530              dopr_unit(i)  = 'K m/s'
     2616             dopr_unit(i)  = TRIM ( heatflux_output_unit )
    25312617             hom(:,2,63,:) = SPREAD( zu, 2, statistic_regions+1 )
    25322618
     
    26212707          CASE ( 'hyp' )
    26222708             dopr_index(i) = 72
    2623              dopr_unit(i)  = 'dbar'
     2709             dopr_unit(i)  = 'hPa'
    26242710             hom(:,2,72,:) = SPREAD( zu, 2, statistic_regions+1 )
     2711
     2712          CASE ( 'rho_air' )
     2713             dopr_index(i)  = 121
     2714             dopr_unit(i)   = 'kg/m3'
     2715             hom(:,2,121,:) = SPREAD( zu, 2, statistic_regions+1 )
     2716
     2717          CASE ( 'rho_air_zw' )
     2718             dopr_index(i)  = 122
     2719             dopr_unit(i)   = 'kg/m3'
     2720             hom(:,2,122,:) = SPREAD( zw, 2, statistic_regions+1 )
    26252721
    26262722          CASE ( 'nr' )
  • palm/trunk/SOURCE/diffusion_e.f90

    r2001 r2037  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! Anelastic approximation implemented
    2323!
    2424! Former revisions:
     
    128128
    129129       USE arrays_3d,                                                          &
    130            ONLY:  dd2zu, ddzu, ddzw, diss, e, km, l_grid, tend, zu, zw
     130           ONLY:  dd2zu, ddzu, ddzw, diss, e, km, l_grid, tend, zu, zw,        &
     131                  drho_air, rho_air_zw
    131132           
    132133       USE control_parameters,                                                 &
     
    222223                                        + (                                    &
    223224               ( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1) &
     225                                                             * rho_air_zw(k)   &
    224226             - ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k)   &
    225                                           ) * ddzw(k)                          &
     227                                                             * rho_air_zw(k-1) &
     228                                          ) * ddzw(k) * drho_air(k)            &
    226229                             - dissipation(k,j)
    227230
     
    294297                                        + (                                    &
    295298               ( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1) &
     299                                                             * rho_air_zw(k)   &
    296300             - ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k)   &
    297                                           ) * ddzw(k)                          &
     301                                                             * rho_air_zw(k-1) &
     302                                          ) * ddzw(k) * drho_air(k)            &
    298303                             - dissipation(k,j)
    299304
     
    339344
    340345       USE arrays_3d,                                                          &
    341            ONLY:  dd2zu, ddzu, ddzw, diss, e, km, l_grid, tend, zu, zw
     346           ONLY:  dd2zu, ddzu, ddzw, diss, e, km, l_grid, tend, zu, zw,        &
     347                  drho_air, rho_air_zw
    342348         
    343349       USE control_parameters,                                                 &
     
    430436                                        + (                                    &
    431437               ( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1) &
     438                                                             * rho_air_zw(k)   &
    432439             - ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k)   &
    433                                           ) * ddzw(k)                          &
     440                                                             * rho_air_zw(k-1) &
     441                                          ) * ddzw(k) * drho_air(k)            &
    434442                                  - dissipation
    435443
     
    497505                                        + (                                    &
    498506               ( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1) &
     507                                                             * rho_air_zw(k)   &
    499508             - ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k)   &
    500                                           ) * ddzw(k)                          &
     509                                                             * rho_air_zw(k-1) &
     510                                          ) * ddzw(k) * drho_air(k)            &
    501511                                  - dissipation
    502512
     
    542552
    543553       USE arrays_3d,                                                          &
    544            ONLY:  dd2zu, ddzu, ddzw, diss, e, km, l_grid, tend, zu, zw
     554           ONLY:  dd2zu, ddzu, ddzw, diss, e, km, l_grid, tend, zu, zw,        &
     555                  drho_air, rho_air_zw
    545556         
    546557       USE control_parameters,                                                 &
     
    625636                                       + (                                    &
    626637              ( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1) &
     638                                                            * rho_air_zw(k)   &
    627639            - ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k)   &
    628                                          ) * ddzw(k)                          &
     640                                                            * rho_air_zw(k-1) &
     641                                         ) * ddzw(k) * drho_air(k)            &
    629642                                       - dissipation(k)
    630643
  • palm/trunk/SOURCE/diffusion_s.f90

    r2001 r2037  
    2020! Current revisions:
    2121! ------------------
    22 !
     22! Anelastic approximation implemented
    2323!
    2424! Former revisions:
     
    113113
    114114       USE arrays_3d,                                                          &
    115            ONLY:  ddzu, ddzw, kh, tend
     115           ONLY:  ddzu, ddzw, kh, tend, drho_air, rho_air_zw
    116116       
    117117       USE control_parameters,                                                 &
     
    190190                                       + 0.5_wp * (                            &
    191191            ( kh(k,j,i) + kh(k+1,j,i) ) * ( s(k+1,j,i)-s(k,j,i) ) * ddzu(k+1)  &
     192                                                            * rho_air_zw(k)    &
    192193          - ( kh(k,j,i) + kh(k-1,j,i) ) * ( s(k,j,i)-s(k-1,j,i) ) * ddzu(k)    &
    193                                                   ) * ddzw(k)
     194                                                            * rho_air_zw(k-1)  &
     195                                                  ) * ddzw(k) * drho_air(k)
    194196             ENDDO
    195197
     
    205207                                                  * ( s(k+1,j,i)-s(k,j,i) )    &
    206208                                                  * ddzu(k+1)                  &
     209                                                  * rho_air_zw(k)              &
    207210                                           + s_flux_b(j,i)                     &
    208                                          ) * ddzw(k)
     211                                         ) * ddzw(k) * drho_air(k)
    209212
    210213             ENDIF
     
    222225                                                    * ( s(k,j,i)-s(k-1,j,i) )  &
    223226                                                    * ddzu(k)                  &
    224                                          ) * ddzw(k)
     227                                                    * rho_air_zw(k-1)          &
     228                                         ) * ddzw(k) * drho_air(k)
    225229
    226230             ENDIF
     
    240244
    241245       USE arrays_3d,                                                          &
    242            ONLY:  ddzu, ddzw, kh, tend
     246           ONLY:  ddzu, ddzw, kh, tend, drho_air, rho_air_zw
    243247           
    244248       USE control_parameters,                                                 &
     
    324328                                       + 0.5_wp * (                            &
    325329            ( kh(k,j,i) + kh(k+1,j,i) ) * ( s(k+1,j,i)-s(k,j,i) ) * ddzu(k+1)  &
     330                                                            * rho_air_zw(k)    &
    326331          - ( kh(k,j,i) + kh(k-1,j,i) ) * ( s(k,j,i)-s(k-1,j,i) ) * ddzu(k)    &
    327                                                   ) * ddzw(k)
     332                                                            * rho_air_zw(k-1)  &
     333                                                  ) * ddzw(k) * drho_air(k)
    328334                ENDIF
    329335             ENDDO
     
    338344                                                     * ( s(k+1,j,i)-s(k,j,i) ) &
    339345                                                     * ddzu(k+1)               &
     346                                                     * rho_air_zw(k)           &
    340347                                              + s_flux_b(j,i)                  &
    341                                             ) * ddzw(k)
     348                                            ) * ddzw(k) * drho_air(k)
    342349                ENDIF
    343350
     
    351358                                                       * ( s(k,j,i)-s(k-1,j,i) )  &
    352359                                                       * ddzu(k)                  &
    353                                             ) * ddzw(k)
     360                                                       * rho_air_zw(k-1)          &
     361                                            ) * ddzw(k) * drho_air(k)
    354362                ENDIF
    355363             ENDDO
     
    370378
    371379       USE arrays_3d,                                                          &
    372            ONLY:  ddzu, ddzw, kh, tend
     380           ONLY:  ddzu, ddzw, kh, tend, drho_air, rho_air_zw
    373381           
    374382       USE control_parameters,                                                 &
     
    446454                                       + 0.5_wp * (                            &
    447455            ( kh(k,j,i) + kh(k+1,j,i) ) * ( s(k+1,j,i)-s(k,j,i) ) * ddzu(k+1)  &
     456                                                            * rho_air_zw(k)    &
    448457          - ( kh(k,j,i) + kh(k-1,j,i) ) * ( s(k,j,i)-s(k-1,j,i) ) * ddzu(k)    &
    449                                                   ) * ddzw(k)
     458                                                            * rho_air_zw(k-1)  &
     459                                                  ) * ddzw(k) * drho_air(k)
    450460       ENDDO
    451461
     
    459469                                               * ( s(k+1,j,i)-s(k,j,i) )       &
    460470                                               * ddzu(k+1)                     &
     471                                               * rho_air_zw(k)                 &
    461472                                        + s_flux_b(j,i)                        &
    462                                       ) * ddzw(k)
     473                                      ) * ddzw(k) * drho_air(k)
    463474
    464475       ENDIF
     
    474485                                               * ( s(k,j,i)-s(k-1,j,i) )       &
    475486                                               * ddzu(k)                       &
    476                                       ) * ddzw(k)
     487                                               * rho_air_zw(k-1)               &
     488                                      ) * ddzw(k) * drho_air(k)
    477489
    478490       ENDIF
  • palm/trunk/SOURCE/diffusion_u.f90

    r2001 r2037  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! Anelastic approximation implemented
    2323!
    2424! Former revisions:
     
    116116
    117117       USE arrays_3d,                                                          &
    118            ONLY:  ddzu, ddzw, km, tend, u, usws, uswst, v, w
     118           ONLY:  ddzu, ddzw, km, tend, u, usws, uswst, v, w,                  &
     119                  drho_air, rho_air_zw
    119120       
    120121       USE control_parameters,                                                 &
     
    217218                      & + ( kmzp * ( ( u(k+1,j,i) - u(k,j,i)   ) * ddzu(k+1)   &
    218219                      &            + ( w(k,j,i)   - w(k,j,i-1) ) * ddx         &
    219                       &            )                                           &
     220                      &            ) * rho_air_zw(k)                           &
    220221                      &   - kmzm * ( ( u(k,j,i)   - u(k-1,j,i)   ) * ddzu(k)   &
    221222                      &            + ( w(k-1,j,i) - w(k-1,j,i-1) ) * ddx       &
    222                       &            )                                           &
    223                       &   ) * ddzw(k)
     223                      &            ) * rho_air_zw(k-1)                         &
     224                      &   ) * ddzw(k) * drho_air(k)
    224225             ENDDO
    225226
     
    242243
    243244                tend(k,j,i) = tend(k,j,i)                                      &
    244                       & + ( kmzp * ( w(k,j,i)   - w(k,j,i-1)   ) * ddx         &
    245                       &   ) * ddzw(k)                                          &
    246                       & + ( kmzp * ( u(k+1,j,i) - u(k,j,i)     ) * ddzu(k+1)   &
    247                       &   + usws(j,i)                                          &
    248                       &   ) * ddzw(k)
     245                      & + ( kmzp * ( ( u(k+1,j,i) - u(k,j,i)   ) * ddzu(k+1)   &
     246                      &            + ( w(k,j,i)   - w(k,j,i-1) ) * ddx         &
     247                      &            ) * rho_air_zw(k)                           &
     248                      &   - ( -usws(j,i) )                                     &
     249                      &   ) * ddzw(k) * drho_air(k)
    249250             ENDIF
    250251
     
    260261
    261262                tend(k,j,i) = tend(k,j,i)                                      &
    262                       & - ( kmzm * ( w(k-1,j,i) - w(k-1,j,i-1) ) * ddx         &
    263                       &   ) * ddzw(k)                                          &
    264                       & + ( -uswst(j,i)                                        &
    265                       &   - kmzm * ( u(k,j,i)   - u(k-1,j,i)   ) * ddzu(k)     &
    266                       &   ) * ddzw(k)
     263                      & + ( ( -uswst(j,i) )                                    &
     264                      &   - kmzm * ( ( u(k,j,i)   - u(k-1,j,i)   ) * ddzu(k)   &
     265                      &            + ( w(k-1,j,i) - w(k-1,j,i-1) ) * ddx       &
     266                      &            ) * rho_air_zw(k-1)                         &
     267                      &   ) * ddzw(k) * drho_air(k)
    267268             ENDIF
    268269
     
    281282
    282283       USE arrays_3d,                                                          &
    283            ONLY:  ddzu, ddzw, km, tend, u, usws, uswst, v, w
     284           ONLY:  ddzu, ddzw, km, tend, u, usws, uswst, v, w,                  &
     285                  drho_air, rho_air_zw
    284286       
    285287       USE control_parameters,                                                 &
     
    390392                         & + ( kmzp * ( ( u(k+1,j,i) - u(k,j,i)   ) * ddzu(k+1)&
    391393                         &            + ( w(k,j,i)   - w(k,j,i-1) ) * ddx      &
    392                          &            )                                        &
     394                         &            ) * rho_air_zw(k)                        &
    393395                         &   - kmzm * ( ( u(k,j,i)   - u(k-1,j,i)   ) * ddzu(k)&
    394396                         &            + ( w(k-1,j,i) - w(k-1,j,i-1) ) * ddx    &
    395                          &            )                                        &
    396                          &   ) * ddzw(k)
     397                         &            ) * rho_air_zw(k-1)                      &
     398                         &   ) * ddzw(k) * drho_air(k)
    397399                ENDIF
    398400             ENDDO
     
    423425
    424426                tend(k,j,i) = tend(k,j,i)                                      &
    425                       & + ( kmzp * ( w(k,j,i)   - w(k,j,i-1)   ) * ddx         &
    426                       &   ) * ddzw(k)                                          &
    427                       & + ( kmzp * ( u(k+1,j,i) - u(k,j,i)     ) * ddzu(k+1)   &
    428                       &   + usws(j,i)                                          &
    429                       &   ) * ddzw(k)
     427                      & + ( kmzp * ( ( u(k+1,j,i) - u(k,j,i)   ) * ddzu(k+1)   &
     428                      &            + ( w(k,j,i)   - w(k,j,i-1) ) * ddx         &
     429                      &            ) * rho_air_zw(k)                           &
     430                      &   - ( -usws(j,i) )                                     &
     431                      &   ) * ddzw(k) * drho_air(k)
    430432             ENDDO
    431433          ENDDO
     
    449451
    450452                tend(k,j,i) = tend(k,j,i)                                      &
    451                       & - ( kmzm * ( w(k-1,j,i) - w(k-1,j,i-1) ) * ddx         &
    452                       &   ) * ddzw(k)                                          &
    453                       & + ( -uswst(j,i)                                        &
    454                       &   - kmzm * ( u(k,j,i)   - u(k-1,j,i)   ) * ddzu(k)     &
    455                       &   ) * ddzw(k)
     453                      & + ( ( -uswst(j,i) )                                    &
     454                      &   - kmzm * ( ( u(k,j,i)   - u(k-1,j,i)   ) * ddzu(k)   &
     455                      &            + ( w(k-1,j,i) - w(k-1,j,i-1) ) * ddx       &
     456                      &            ) * rho_air_zw(k-1)                         &
     457                      &   ) * ddzw(k) * drho_air(k)
    456458             ENDDO
    457459          ENDDO
     
    471473
    472474       USE arrays_3d,                                                          &
    473            ONLY:  ddzu, ddzw, km, tend, u, usws, uswst, v, w
     475           ONLY:  ddzu, ddzw, km, tend, u, usws, uswst, v, w,                  &
     476                  drho_air, rho_air_zw
    474477       
    475478       USE control_parameters,                                                 &
     
    559562                      & + ( kmzp * ( ( u(k+1,j,i) - u(k,j,i)   ) * ddzu(k+1)   &
    560563                      &            + ( w(k,j,i)   - w(k,j,i-1) ) * ddx         &
    561                       &            )                                           &
     564                      &            ) * rho_air_zw(k)                           &
    562565                      &   - kmzm * ( ( u(k,j,i)   - u(k-1,j,i)   ) * ddzu(k)   &
    563566                      &            + ( w(k-1,j,i) - w(k-1,j,i-1) ) * ddx       &
    564                       &            )                                           &
    565                       &   ) * ddzw(k)
     567                      &            ) * rho_air_zw(k-1)                         &
     568                      &   ) * ddzw(k) * drho_air(k)
    566569       ENDDO
    567570
     
    583586
    584587          tend(k,j,i) = tend(k,j,i)                                            &
    585                       & + ( kmzp * ( w(k,j,i)   - w(k,j,i-1)   ) * ddx         &
    586                       &   ) * ddzw(k)                                          &
    587                       & + ( kmzp * ( u(k+1,j,i) - u(k,j,i)     ) * ddzu(k+1)   &
    588                       &   + usws(j,i)                                          &
    589                       &   ) * ddzw(k)
     588                      & + ( kmzp * ( ( u(k+1,j,i) - u(k,j,i)   ) * ddzu(k+1)   &
     589                      &            + ( w(k,j,i)   - w(k,j,i-1) ) * ddx         &
     590                      &            ) * rho_air_zw(k)                           &
     591                      &   - ( -usws(j,i) )                                     &
     592                      &   ) * ddzw(k) * drho_air(k)
    590593       ENDIF
    591594
     
    597600!
    598601!--       Interpolate eddy diffusivities on staggered gridpoints
    599           kmzm = 0.25_wp *                                                     &
    600                  ( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
     602          kmzm = 0.25_wp * ( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
    601603
    602604          tend(k,j,i) = tend(k,j,i)                                            &
    603                       & - ( kmzm * ( w(k-1,j,i) - w(k-1,j,i-1) ) * ddx         &
    604                       &   ) * ddzw(k)                                          &
    605                       & + ( -uswst(j,i)                                        &
    606                       &   - kmzm * ( u(k,j,i)   - u(k-1,j,i)   ) * ddzu(k)     &
    607                       &   ) * ddzw(k)
     605                      & + ( ( -uswst(j,i) )                                    &
     606                      &   - kmzm * ( ( u(k,j,i)   - u(k-1,j,i)   ) * ddzu(k)   &
     607                      &            + ( w(k-1,j,i) - w(k-1,j,i-1) ) * ddx       &
     608                      &            ) * rho_air_zw(k-1)                         &
     609                      &   ) * ddzw(k) * drho_air(k)
    608610       ENDIF
    609611
  • palm/trunk/SOURCE/diffusion_v.f90

    r2001 r2037  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! Anelastic approximation implemented
    2323!
    2424! Former revisions:
     
    111111
    112112       USE arrays_3d,                                                          &
    113            ONLY:  ddzu, ddzw, km, tend, u, v, vsws, vswst, w
     113           ONLY:  ddzu, ddzw, km, tend, u, v, vsws, vswst, w,                  &
     114                  drho_air, rho_air_zw
    114115       
    115116       USE control_parameters,                                                 &
     
    212213                      & + ( kmzp * ( ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1)     &
    213214                      &            + ( w(k,j,i) - w(k,j-1,i) ) * ddy           &
    214                       &            )                                           &
     215                      &            ) * rho_air_zw(k)                           &
    215216                      &   - kmzm * ( ( v(k,j,i)   - v(k-1,j,i)   ) * ddzu(k)   &
    216217                      &            + ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy       &
    217                       &            )                                           &
    218                       &   ) * ddzw(k)
     218                      &            ) * rho_air_zw(k-1)                         &
     219                      &   ) * ddzw(k) * drho_air(k)
    219220             ENDDO
    220221
     
    237238
    238239                tend(k,j,i) = tend(k,j,i)                                      &
    239                       & + ( kmzp * ( w(k,j,i) - w(k,j-1,i)     ) * ddy         &
    240                       &   ) * ddzw(k)                                          &
    241                       & + ( kmzp * ( v(k+1,j,i) - v(k,j,i)     ) * ddzu(k+1)   &
    242                       &   + vsws(j,i)                                          &
    243                       &   ) * ddzw(k)
     240                      & + ( kmzp * ( ( v(k+1,j,i) - v(k,j,i)   ) * ddzu(k+1)   &
     241                      &            + ( w(k,j,i)   - w(k,j-1,i) ) * ddy         &
     242                      &            ) * rho_air_zw(k)                           &
     243                      &   - ( -vsws(j,i) )                                     &
     244                      &   ) * ddzw(k) * drho_air(k)
    244245             ENDIF
    245246
     
    255256
    256257                tend(k,j,i) = tend(k,j,i)                                      &
    257                       & - ( kmzm *  ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy        &
    258                       &   ) * ddzw(k)                                          &
    259                       & + ( -vswst(j,i)                                        &
    260                       &   - kmzm * ( v(k,j,i)   - v(k-1,j,i)    ) * ddzu(k)    &
    261                       &   ) * ddzw(k)
     258                      & + ( ( -vswst(j,i) )                                    &
     259                      &   - kmzm * ( ( v(k,j,i)   - v(k-1,j,i)   ) * ddzu(k)   &
     260                      &            + ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy       &
     261                      &            ) * rho_air_zw(k-1)                         &
     262                      &   ) * ddzw(k) * drho_air(k)
    262263             ENDIF
    263264
     
    276277
    277278       USE arrays_3d,                                                          &
    278            ONLY:  ddzu, ddzw, km, tend, u, v, vsws, vswst, w
     279           ONLY:  ddzu, ddzw, km, tend, u, v, vsws, vswst, w,                  &
     280                  drho_air, rho_air_zw
    279281       
    280282       USE control_parameters,                                                 &
     
    385387                         & + ( kmzp * ( ( v(k+1,j,i) - v(k,j,i)   ) * ddzu(k+1)&
    386388                         &            + ( w(k,j,i)   - w(k,j-1,i) ) * ddy      &
    387                          &            )                                        &
     389                         &            ) * rho_air_zw(k)                        &
    388390                         &   - kmzm * ( ( v(k,j,i)   - v(k-1,j,i)   ) * ddzu(k)&
    389391                         &            + ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy    &
    390                          &            )                                        &
    391                          &   ) * ddzw(k)
     392                         &            ) * rho_air_zw(k-1)                      &
     393                         &   ) * ddzw(k) * drho_air(k)
    392394                ENDIF
    393395             ENDDO
     
    418420
    419421                tend(k,j,i) = tend(k,j,i)                                      &
    420                       & + ( kmzp * ( w(k,j,i) - w(k,j-1,i)     ) * ddy         &
    421                       &   ) * ddzw(k)                                          &
    422                       & + ( kmzp * ( v(k+1,j,i) - v(k,j,i)     ) * ddzu(k+1)   &
    423                       &   + vsws(j,i)                                          &
    424                       &   ) * ddzw(k)
     422                      & + ( kmzp * ( ( v(k+1,j,i) - v(k,j,i)   ) * ddzu(k+1)   &
     423                      &            + ( w(k,j,i)   - w(k,j-1,i) ) * ddy         &
     424                      &            ) * rho_air_zw(k)                           &
     425                      &   - ( -vsws(j,i) )                                     &
     426                      &   ) * ddzw(k) * drho_air(k)
    425427             ENDDO
    426428          ENDDO
     
    444446
    445447                tend(k,j,i) = tend(k,j,i)                                      &
    446                       & - ( kmzm *  ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy        &
    447                       &   ) * ddzw(k)                                          &
    448                       & + ( -vswst(j,i)                                        &
    449                       &   - kmzm * ( v(k,j,i)   - v(k-1,j,i)    ) * ddzu(k)    &
    450                       &   ) * ddzw(k)
     448                      & + ( ( -vswst(j,i) )                                    &
     449                      &   - kmzm * ( ( v(k,j,i)   - v(k-1,j,i)   ) * ddzu(k)   &
     450                      &            + ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy       &
     451                      &            ) * rho_air_zw(k-1)                         &
     452                      &   ) * ddzw(k) * drho_air(k)
    451453             ENDDO
    452454          ENDDO
     
    466468
    467469       USE arrays_3d,                                                          &
    468            ONLY:  ddzu, ddzw, km, tend, u, v, vsws, vswst, w
     470           ONLY:  ddzu, ddzw, km, tend, u, v, vsws, vswst, w,                  &
     471                  drho_air, rho_air_zw
    469472       
    470473       USE control_parameters,                                                 &
     
    556559                      & + ( kmzp * ( ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1)     &
    557560                      &            + ( w(k,j,i) - w(k,j-1,i) ) * ddy           &
    558                       &            )                                           &
     561                      &            ) * rho_air_zw(k)                           &
    559562                      &   - kmzm * ( ( v(k,j,i)   - v(k-1,j,i)   ) * ddzu(k)   &
    560563                      &            + ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy       &
    561                       &            )                                           &
    562                       &   ) * ddzw(k)
     564                      &            ) * rho_air_zw(k-1)                         &
     565                      &   ) * ddzw(k) * drho_air(k)
    563566       ENDDO
    564567
     
    580583
    581584          tend(k,j,i) = tend(k,j,i)                                            &
    582                       & + ( kmzp * ( w(k,j,i) - w(k,j-1,i)     ) * ddy         &
    583                       &   ) * ddzw(k)                                          &
    584                       & + ( kmzp * ( v(k+1,j,i) - v(k,j,i)     ) * ddzu(k+1)   &
    585                       &   + vsws(j,i)                                          &
    586                       &   ) * ddzw(k)
     585                      & + ( kmzp * ( ( v(k+1,j,i) - v(k,j,i)   ) * ddzu(k+1)   &
     586                      &            + ( w(k,j,i)   - w(k,j-1,i) ) * ddy         &
     587                      &            ) * rho_air_zw(k)                           &
     588                      &   - ( -vsws(j,i) )                                     &
     589                      &   ) * ddzw(k) * drho_air(k)
    587590       ENDIF
    588591
     
    594597!
    595598!--       Interpolate eddy diffusivities on staggered gridpoints
    596           kmzm = 0.25_wp * &
    597                  ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) )
     599          kmzm = 0.25_wp * ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) )
    598600
    599601          tend(k,j,i) = tend(k,j,i)                                            &
    600                       & - ( kmzm *  ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy        &
    601                       &   ) * ddzw(k)                                          &
    602                       & + ( -vswst(j,i)                                        &
    603                       &   - kmzm * ( v(k,j,i)   - v(k-1,j,i)    ) * ddzu(k)    &
    604                       &   ) * ddzw(k)
     602                      & + ( ( -vswst(j,i) )                                    &
     603                      &   - kmzm * ( ( v(k,j,i)   - v(k-1,j,i)   ) * ddzu(k)   &
     604                      &            + ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy       &
     605                      &            ) * rho_air_zw(k-1)                         &
     606                      &   ) * ddzw(k) * drho_air(k)
    605607       ENDIF
    606608
  • palm/trunk/SOURCE/diffusion_w.f90

    r2001 r2037  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! Anelastic approximation implemented
    2323!
    2424! Former revisions:
     
    119119
    120120       USE arrays_3d,                                                          &         
    121            ONLY :  ddzu, ddzw, km, tend, u, v, w
     121           ONLY :  ddzu, ddzw, km, tend, u, v, w, drho_air_zw, rho_air
    122122           
    123123       USE control_parameters,                                                 &
     
    184184                      & + 2.0_wp * (                                           &
    185185                      &   km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1)  &
     186                      &               * rho_air(k+1)                           &
    186187                      & - km(k,j,i)   * ( w(k,j,i)   - w(k-1,j,i) ) * ddzw(k)  &
    187                       &            ) * ddzu(k+1)
     188                      &               * rho_air(k)                             &
     189                      &            ) * ddzu(k+1) * drho_air_zw(k)
    188190             ENDDO
    189191
     
    227229                                 + 2.0_wp * (                                  &
    228230                           km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1) &
     231                                       * rho_air(k+1)                          &
    229232                         - km(k,j,i)   * ( w(k,j,i)   - w(k-1,j,i) ) * ddzw(k) &
    230                                             ) * ddzu(k+1)
     233                                       * rho_air(k)                            &
     234                                            ) * ddzu(k+1) * drho_air_zw(k)
    231235                ENDDO
    232236             ENDIF
     
    246250
    247251       USE arrays_3d,                                                          &
    248            ONLY :  ddzu, ddzw, km, tend, u, v, w
     252           ONLY :  ddzu, ddzw, km, tend, u, v, w, drho_air_zw, rho_air
    249253           
    250254       USE control_parameters,                                                 &
     
    316320                         & + 2.0_wp * (                                          &
    317321                         &   km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1) &
     322                         &               * rho_air(k+1)                          &
    318323                         & - km(k,j,i)   * ( w(k,j,i)   - w(k-1,j,i) ) * ddzw(k) &
    319                          &            ) * ddzu(k+1)
     324                         &               * rho_air(k)                            &
     325                         &            ) * ddzu(k+1) * drho_air_zw(k)
    320326                ENDIF
    321327             ENDDO
     
    361367                                 + 2.0_wp * (                                  &
    362368                           km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1) &
     369                                       * rho_air(k+1)                          &
    363370                         - km(k,j,i)   * ( w(k,j,i)   - w(k-1,j,i) ) * ddzw(k) &
    364                                             ) * ddzu(k+1)
     371                                       * rho_air(k)                            &
     372                                            ) * ddzu(k+1) * drho_air_zw(k)
    365373                ENDIF
    366374             ENDDO
     
    381389
    382390       USE arrays_3d,                                                          &         
    383            ONLY :  ddzu, ddzw, km, tend, u, v, w
     391           ONLY :  ddzu, ddzw, km, tend, u, v, w, drho_air_zw, rho_air
    384392           
    385393       USE control_parameters,                                                 &
     
    430438                      & + 2.0_wp * (                                           &
    431439                      &   km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1)  &
     440                      &               * rho_air(k+1)                           &
    432441                      & - km(k,j,i)   * ( w(k,j,i)   - w(k-1,j,i) ) * ddzw(k)  &
    433                       &            ) * ddzu(k+1)
     442                      &               * rho_air(k)                             &
     443                      &            ) * ddzu(k+1) * drho_air_zw(k)
    434444       ENDDO
    435445
     
    485495                                 + 2.0_wp * (                                  &
    486496                           km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1) &
     497                                       * rho_air(k+1)                          &
    487498                         - km(k,j,i)   * ( w(k,j,i)   - w(k-1,j,i) ) * ddzw(k) &
    488                                             ) * ddzu(k+1)
     499                                       * rho_air(k)                            &
     500                                            ) * ddzu(k+1) * drho_air_zw(k)
    489501          ENDDO
    490502       ENDIF
  • palm/trunk/SOURCE/disturb_heatflux.f90

    r2001 r2037  
    2020! Current revisions:
    2121! ------------------
    22 !
     22! Anelastic approximation implemented
    2323!
    2424! Former revisions:
     
    6565
    6666    USE arrays_3d,                                                             &
    67         ONLY:  shf
     67        ONLY:  shf, heatflux_input_conversion
    6868       
    6969    USE control_parameters,                                                    &
     
    7676   
    7777    USE indices,                                                               &
    78         ONLY:  nx, nxl, nxr, ny, nyn, nys, nzb_s_inner
     78        ONLY:  nx, nxl, nxr, ny, nyn, nys, nzb, nzb_s_inner
    7979
    8080    IMPLICIT NONE
     
    9797          THEN
    9898             IF ( nzb_s_inner(j,i) == 0 )  THEN
    99                 shf(j,i) = randomnumber * surface_heatflux
     99                shf(j,i) = randomnumber * surface_heatflux                     &
     100                           * heatflux_input_conversion(nzb)
    100101             ELSE
    101102!
    102103!--             Over topography surface_heatflux is replaced by wall_heatflux(0)
    103                 shf(j,i) = randomnumber * wall_heatflux(0)
     104                shf(j,i) = randomnumber * wall_heatflux(0)                     &
     105                           * heatflux_input_conversion(nzb_s_inner(j,i))
    104106             ENDIF
    105107          ENDIF
  • palm/trunk/SOURCE/flow_statistics.f90

    r2032 r2037  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! Anelastic approximation implemented
    2323!
    2424! Former revisions:
     
    214214
    215215    USE arrays_3d,                                                             &
    216         ONLY:  ddzu, ddzw, e, hyp, km, kh, nr, ol, p, prho, prr, pt, q, qc, ql,&
    217                qr, qs, qsws, qswst, rho_ocean, s, sa, ss, ssws, sswst, saswsb,       &
    218                saswst, shf, td_lsa_lpt, td_lsa_q, td_sub_lpt, td_sub_q,        &
    219                time_vert, ts, tswst, u, ug, us, usws, uswst, vsws, v, vg, vpt, &
    220                vswst, w, w_subs, zw
     216        ONLY:  ddzu, ddzw, e, heatflux_output_conversion, hyp, km, kh,         &
     217               momentumflux_output_conversion, nr, ol, p, prho, prr, pt, q,    &
     218               qc, ql, qr, qs, qsws, qswst, rho_air, rho_air_zw, rho_ocean, s, &
     219               sa, ss, ssws, sswst, saswsb, saswst, shf, td_lsa_lpt, td_lsa_q, &
     220               td_sub_lpt, td_sub_q, time_vert, ts, tswst, u, ug, us, usws,    &
     221               uswst, vsws, v, vg, vpt, vswst, w, w_subs,                      &
     222               waterflux_output_conversion, zw
    221223       
    222224    USE cloud_parameters,                                                      &
     
    326328       sums_l(:,11,:) = sums_l_l(:,sr,:)      ! mixing length from diffusivities
    327329!--    WARNING: next line still has to be adjusted for OpenMP
    328        sums_l(:,21,0) = sums_wsts_bc_l(:,sr)  ! heat flux from advec_s_bc
     330       sums_l(:,21,0) = sums_wsts_bc_l(:,sr) *                                 &
     331                        heatflux_output_conversion  ! heat flux from advec_s_bc
    329332       sums_l(nzb+9,pr_palm,0)  = sums_divold_l(sr)  ! old divergence from pres
    330333       sums_l(nzb+10,pr_palm,0) = sums_divnew_l(sr)  ! new divergence from pres
     
    358361!
    359362!--          Swap the turbulent quantities evaluated in advec_ws.
    360              sums_l(:,13,i) = sums_wsus_ws_l(:,i)       ! w*u*
    361              sums_l(:,15,i) = sums_wsvs_ws_l(:,i)       ! w*v*
     363             sums_l(:,13,i) = sums_wsus_ws_l(:,i)                              &
     364                              * momentumflux_output_conversion ! w*u*
     365             sums_l(:,15,i) = sums_wsvs_ws_l(:,i)                              &
     366                              * momentumflux_output_conversion ! w*v*
    362367             sums_l(:,30,i) = sums_us2_ws_l(:,i)        ! u*2
    363368             sums_l(:,31,i) = sums_vs2_ws_l(:,i)        ! v*2
     
    373378
    374379          DO  i = 0, threads_per_task-1
    375              sums_l(:,17,i)                        = sums_wspts_ws_l(:,i) ! w*pt*
     380             sums_l(:,17,i)                        = sums_wspts_ws_l(:,i)      &
     381                                           * heatflux_output_conversion  ! w*pt*
    376382             IF ( ocean          ) sums_l(:,66,i)  = sums_wssas_ws_l(:,i) ! w*sa*
    377              IF ( humidity       ) sums_l(:,49,i)  = sums_wsqs_ws_l(:,i)  ! w*q*
     383             IF ( humidity       ) sums_l(:,49,i)  = sums_wsqs_ws_l(:,i)       &
     384                                           * waterflux_output_conversion  ! w*q*
    378385             IF ( passive_scalar ) sums_l(:,116,i) = sums_wsss_ws_l(:,i)  ! w*s*
    379386          ENDDO
     
    722729                                   ( u(k+1,j,i) - u(k,j,i)   ) * ddzu(k+1)     &
    723730                                 + ( w(k,j,i)   - w(k,j,i-1) ) * ddx           &
    724                                                                ) * rmask(j,i,sr)
     731                                                           ) * rmask(j,i,sr)   &
     732                                         * rho_air_zw(k)                       &
     733                                         * momentumflux_output_conversion(k)
    725734!
    726735!--             Momentum flux w"v"
     
    730739                                   ( v(k+1,j,i) - v(k,j,i)   ) * ddzu(k+1)     &
    731740                                 + ( w(k,j,i)   - w(k,j-1,i) ) * ddy           &
    732                                                                ) * rmask(j,i,sr)
     741                                                           ) * rmask(j,i,sr)   &
     742                                         * rho_air_zw(k)                       &
     743                                         * momentumflux_output_conversion(k)
    733744!
    734745!--             Heat flux w"pt"
     
    736747                                         - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )&
    737748                                               * ( pt(k+1,j,i) - pt(k,j,i) )   &
     749                                               * rho_air_zw(k)                 &
     750                                               * heatflux_output_conversion(k) &
    738751                                               * ddzu(k+1) * rmask(j,i,sr)
    739752
     
    754767                                         - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )&
    755768                                               * ( vpt(k+1,j,i) - vpt(k,j,i) ) &
     769                                               * rho_air_zw(k)                 &
     770                                               * heatflux_output_conversion(k) &
    756771                                               * ddzu(k+1) * rmask(j,i,sr)
    757772                   sums_l(k,48,tn) = sums_l(k,48,tn)                           &
    758773                                         - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )&
    759774                                               * ( q(k+1,j,i) - q(k,j,i) )     &
     775                                               * rho_air_zw(k)                 &
     776                                               * waterflux_output_conversion(k)&
    760777                                               * ddzu(k+1) * rmask(j,i,sr)
    761778
     
    765782                                               * ( ( q(k+1,j,i) - ql(k+1,j,i) )&
    766783                                                - ( q(k,j,i) - ql(k,j,i) ) )   &
     784                                               * rho_air_zw(k)                 &
     785                                               * waterflux_output_conversion(k)&
    767786                                               * ddzu(k+1) * rmask(j,i,sr)
    768787                   ENDIF
     
    784803             IF ( use_surface_fluxes )  THEN
    785804                sums_l(nzb,12,tn) = sums_l(nzb,12,tn) + &
     805                                    momentumflux_output_conversion(nzb) * &
    786806                                    usws(j,i) * rmask(j,i,sr)     ! w"u"
    787807                sums_l(nzb,14,tn) = sums_l(nzb,14,tn) + &
     808                                    momentumflux_output_conversion(nzb) * &
    788809                                    vsws(j,i) * rmask(j,i,sr)     ! w"v"
    789810                sums_l(nzb,16,tn) = sums_l(nzb,16,tn) + &
     811                                    heatflux_output_conversion(nzb) * &
    790812                                    shf(j,i)  * rmask(j,i,sr)     ! w"pt"
    791813                sums_l(nzb,58,tn) = sums_l(nzb,58,tn) + &
     
    799821                IF ( humidity )  THEN
    800822                   sums_l(nzb,48,tn) = sums_l(nzb,48,tn) +                     &
     823                                       waterflux_output_conversion(nzb) *      &
    801824                                       qsws(j,i) * rmask(j,i,sr)  ! w"q" (w"qv")
    802825                   sums_l(nzb,45,tn) = sums_l(nzb,45,tn) + (                   &
    803826                                       ( 1.0_wp + 0.61_wp * q(nzb,j,i) ) *     &
    804827                                       shf(j,i) + 0.61_wp * pt(nzb,j,i) *      &
    805                                                   qsws(j,i) )
     828                                                  qsws(j,i) )                  &
     829                                       * heatflux_output_conversion(nzb)
    806830                   IF ( cloud_droplets )  THEN
    807831                      sums_l(nzb,45,tn) = sums_l(nzb,45,tn) + (                &
    808832                                         ( 1.0_wp + 0.61_wp * q(nzb,j,i) -     &
    809833                                           ql(nzb,j,i) ) * shf(j,i) +          &
    810                                            0.61_wp * pt(nzb,j,i) * qsws(j,i) )
     834                                           0.61_wp * pt(nzb,j,i) * qsws(j,i) ) &
     835                                          * heatflux_output_conversion(nzb)
    811836                   ENDIF
    812837                   IF ( cloud_physics )  THEN
    813838!
    814839!--                   Formula does not work if ql(nzb) /= 0.0
    815                       sums_l(nzb,51,tn) = sums_l(nzb,51,tn) +                  &
     840                      sums_l(nzb,51,tn) = sums_l(nzb,51,tn) +                  &
     841                                          waterflux_output_conversion(nzb) *   &
    816842                                          qsws(j,i) * rmask(j,i,sr) ! w"q" (w"qv")
    817843                   ENDIF
     
    860886             IF ( use_top_fluxes )  THEN
    861887                sums_l(nzt:nzt+1,12,tn) = sums_l(nzt:nzt+1,12,tn) + &
     888                                    momentumflux_output_conversion(nzt:nzt+1) * &
    862889                                    uswst(j,i) * rmask(j,i,sr)    ! w"u"
    863890                sums_l(nzt:nzt+1,14,tn) = sums_l(nzt:nzt+1,14,tn) + &
     891                                    momentumflux_output_conversion(nzt:nzt+1) * &
    864892                                    vswst(j,i) * rmask(j,i,sr)    ! w"v"
    865893                sums_l(nzt:nzt+1,16,tn) = sums_l(nzt:nzt+1,16,tn) + &
     894                                    heatflux_output_conversion(nzt:nzt+1) * &
    866895                                    tswst(j,i)  * rmask(j,i,sr)   ! w"pt"
    867896                sums_l(nzt:nzt+1,58,tn) = sums_l(nzt:nzt+1,58,tn) + &
     
    876905                IF ( humidity )  THEN
    877906                   sums_l(nzt,48,tn) = sums_l(nzt,48,tn) +                     &
     907                                       waterflux_output_conversion(nzt) *      &
    878908                                       qswst(j,i) * rmask(j,i,sr) ! w"q" (w"qv")
    879909                   sums_l(nzt,45,tn) = sums_l(nzt,45,tn) + (                   &
    880910                                       ( 1.0_wp + 0.61_wp * q(nzt,j,i) ) *     &
    881911                                       tswst(j,i) + 0.61_wp * pt(nzt,j,i) *    &
    882                                                              qswst(j,i) )
     912                                                             qswst(j,i) )      &
     913                                       * heatflux_output_conversion(nzt)
    883914                   IF ( cloud_droplets )  THEN
    884915                      sums_l(nzt,45,tn) = sums_l(nzt,45,tn) + (                &
    885916                                          ( 1.0_wp + 0.61_wp * q(nzt,j,i) -    &
    886917                                            ql(nzt,j,i) ) * tswst(j,i) +       &
    887                                             0.61_wp * pt(nzt,j,i) * qswst(j,i) )
     918                                           0.61_wp * pt(nzt,j,i) * qswst(j,i) )&
     919                                           * heatflux_output_conversion(nzt)
    888920                   ENDIF
    889921                   IF ( cloud_physics )  THEN
     
    891923!--                   Formula does not work if ql(nzb) /= 0.0
    892924                      sums_l(nzt,51,tn) = sums_l(nzt,51,tn) + &   ! w"q" (w"qv")
     925                                          waterflux_output_conversion(nzt) *   &
    893926                                          qswst(j,i) * rmask(j,i,sr)
    894927                   ENDIF
     
    943976                                    vpt(k+1,j,i) - hom(k+1,1,44,sr) )
    944977                      sums_l(k,46,tn) = sums_l(k,46,tn) + pts * w(k,j,i) *     &
     978                                               heatflux_output_conversion(k) * &
    945979                                                          rmask(j,i,sr)
    946980                      sums_l(k,54,tn) = sums_l(k,54,tn) + ql(k,j,i) * rmask(j,i,sr)
     
    953987                              hom(k+1,1,42,sr) )
    954988                         sums_l(k,52,tn) = sums_l(k,52,tn) + pts * w(k,j,i) *  &
     989                                             waterflux_output_conversion(k) *  &
    955990                                                             rmask(j,i,sr)
    956991                         sums_l(k,75,tn) = sums_l(k,75,tn) + qc(k,j,i) *       &
     
    9711006                                          vpt(k+1,j,i) - hom(k+1,1,44,sr) )
    9721007                         sums_l(k,46,tn) = sums_l(k,46,tn) + pts * w(k,j,i) *  &
     1008                                              heatflux_output_conversion(k) *  &
    9731009                                                             rmask(j,i,sr)
    9741010                      ELSE IF ( ws_scheme_sca .AND. sr == 0 )  THEN
    975                          sums_l(k,46,tn) = ( 1.0_wp + 0.61_wp *                &
    976                                              hom(k,1,41,sr) ) *                &
    977                                            sums_l(k,17,tn) +                   &
    978                                            0.61_wp * hom(k,1,4,sr) *           &
    979                                            sums_l(k,49,tn)
     1011                         sums_l(k,46,tn) = ( ( 1.0_wp + 0.61_wp *              &
     1012                                               hom(k,1,41,sr) ) *              &
     1013                                             sums_l(k,17,tn) +                 &
     1014                                             0.61_wp * hom(k,1,4,sr) *         &
     1015                                             sums_l(k,49,tn)                   &
     1016                                           ) * heatflux_output_conversion(k)
    9801017                      END IF
    9811018                   END IF
     
    9961033                sums_l(k,37,tn) = sums_l(k,37,tn) + w(k,j,i) * 0.5_wp *        &
    9971034                                             ( ust**2 + vst**2 + w(k,j,i)**2 ) &
     1035                                           * momentumflux_output_conversion(k) &
    9981036                                             * rmask(j,i,sr)
    9991037             ENDDO
     
    10171055                  sums_l(k,13,tn) = sums_l(k,13,tn) + 0.5_wp *                 &
    10181056                                                    ( w(k,j,i-1) + w(k,j,i) )  &
     1057                                          * momentumflux_output_conversion(k)  &
    10191058                                                    * ust * rmask(j,i,sr)
    10201059!
     
    10221061                  sums_l(k,15,tn) = sums_l(k,15,tn) + 0.5_wp *                 &
    10231062                                                    ( w(k,j-1,i) + w(k,j,i) )  &
     1063                                          * momentumflux_output_conversion(k)  &
    10241064                                                    * vst * rmask(j,i,sr)
    10251065               ENDDO
     
    10381078                           ( pt(k,j,i)   - hom(k,1,4,sr) +                     &
    10391079                             pt(k+1,j,i) - hom(k+1,1,4,sr) )                   &
     1080                           * heatflux_output_conversion(k)                     &
    10401081                           * w(k,j,i) * rmask(j,i,sr)
    10411082                  IF ( humidity )  THEN
     
    10431084                                      q(k+1,j,i) - hom(k+1,1,41,sr) )
    10441085                     sums_l(k,49,tn) = sums_l(k,49,tn) + pts * w(k,j,i) *      &
     1086                                       waterflux_output_conversion(k) *        &
    10451087                                       rmask(j,i,sr)
    10461088                  ENDIF
     
    11471189                                                   ( kh(k,j,i) + kh(k,j,i-1) ) &
    11481190                                                 * ( pt(k,j,i-1) - pt(k,j,i) ) &
     1191                                               * rho_air_zw(k)                 &
     1192                                               * heatflux_output_conversion(k) &
    11491193                                                 * ddx * rmask(j,i,sr)
    11501194                   sums_l(k,61,tn) = sums_l(k,61,tn) - 0.5_wp *                &
    11511195                                                   ( kh(k,j,i) + kh(k,j-1,i) ) &
    11521196                                                 * ( pt(k,j-1,i) - pt(k,j,i) ) &
     1197                                               * rho_air_zw(k)                 &
     1198                                               * heatflux_output_conversion(k) &
    11531199                                                 * ddy * rmask(j,i,sr)
    11541200!
     
    11571203                                                  ( u(k,j,i) - hom(k,1,1,sr) ) &
    11581204                                    * 0.5_wp * ( pt(k,j,i-1) - hom(k,1,4,sr) + &
    1159                                                  pt(k,j,i)   - hom(k,1,4,sr) )
     1205                                                 pt(k,j,i)   - hom(k,1,4,sr) ) &
     1206                                               * heatflux_output_conversion(k)
    11601207                   pts = 0.5_wp * ( pt(k,j-1,i) - hom(k,1,4,sr) +              &
    11611208                                    pt(k,j,i)   - hom(k,1,4,sr) )
     
    11631210                                                  ( v(k,j,i) - hom(k,1,2,sr) ) &
    11641211                                    * 0.5_wp * ( pt(k,j-1,i) - hom(k,1,4,sr) + &
    1165                                                  pt(k,j,i)   - hom(k,1,4,sr) )
     1212                                                 pt(k,j,i)   - hom(k,1,4,sr) ) &
     1213                                               * heatflux_output_conversion(k)
    11661214                ENDDO
    11671215             ENDDO
     
    14891537       ENDIF
    14901538
     1539       hom(:,1,121,sr) = rho_air       ! rho_air in Kg/m^3
     1540       hom(:,1,122,sr) = rho_air_zw    ! rho_air_zw in Kg/m^3
     1541
    14911542       hom(:,1,pr_palm,sr) =   sums(:,pr_palm)
    14921543                                       ! u*, w'u', w'v', t* (in last profile)
     
    15971648       THEN
    15981649          hom(nzb+8,1,pr_palm,sr) = &
    1599              ( g / hom(k_surface_level+1,1,4,sr) * hom(k_surface_level,1,18,sr)&
     1650             ( g / hom(k_surface_level+1,1,4,sr) *                             &
     1651             ( hom(k_surface_level,1,18,sr) / heatflux_output_conversion(nzb) )&
    16001652             * ABS( z_i(1) - mean_surface_level_height(sr) ) )**0.333333333_wp
    16011653       ELSE
     
    17081760
    17091761    USE arrays_3d,                                                             &
    1710         ONLY:  ddzu, ddzw, e, hyp, km, kh, nr, p, prho, pt, q, qc, ql, qr, qs, &
    1711                qsws, qswst, rho_ocean, s, sa, saswsb, saswst, shf, ss, ssws, sswst,  &
    1712                td_lsa_lpt, td_lsa_q, td_sub_lpt, td_sub_q, time_vert, ts,      &
    1713                tswst, u, ug, us, usws, uswst, vsws, v, vg, vpt, vswst, w,      &
    1714                w_subs, zw
     1762        ONLY:  ddzu, ddzw, e, heatflux_output_conversion, hyp, km, kh,         &
     1763               momentumflux_output_conversion, nr, p, prho, pt, q, qc, ql, qr, &
     1764               qs, qsws, qswst, rho_air, rho_air_zw, rho_ocean, s, sa, saswsb, &
     1765               saswst, shf, ss, ssws, sswst, td_lsa_lpt, td_lsa_q, td_sub_lpt, &
     1766               td_sub_q, time_vert, ts, tswst, u, ug, us, usws, uswst, vsws,   &
     1767               v, vg, vpt, vswst, w, w_subs, waterflux_output_conversion, zw
    17151768                 
    17161769       
     
    18281881       sums_l(:,11,:) = sums_l_l(:,sr,:)      ! mixing length from diffusivities
    18291882!--    WARNING: next line still has to be adjusted for OpenMP
    1830        sums_l(:,21,0) = sums_wsts_bc_l(:,sr)  ! heat flux from advec_s_bc
     1883       sums_l(:,21,0) = sums_wsts_bc_l(:,sr) *                                 &
     1884                        heatflux_output_conversion  ! heat flux from advec_s_bc
    18311885       sums_l(nzb+9,pr_palm,0)  = sums_divold_l(sr)  ! old divergence from pres
    18321886       sums_l(nzb+10,pr_palm,0) = sums_divnew_l(sr)  ! new divergence from pres
     
    18601914!
    18611915!--          Swap the turbulent quantities evaluated in advec_ws.
    1862              sums_l(:,13,i) = sums_wsus_ws_l(:,i)       ! w*u*
    1863              sums_l(:,15,i) = sums_wsvs_ws_l(:,i)       ! w*v*
     1916             sums_l(:,13,i) = sums_wsus_ws_l(:,i)                              &
     1917                              * momentumflux_output_conversion ! w*u*
     1918             sums_l(:,15,i) = sums_wsvs_ws_l(:,i)                              &
     1919                              * momentumflux_output_conversion ! w*v*
    18641920             sums_l(:,30,i) = sums_us2_ws_l(:,i)        ! u*2
    18651921             sums_l(:,31,i) = sums_vs2_ws_l(:,i)        ! v*2
     
    18811937
    18821938          DO  i = 0, threads_per_task-1
    1883              sums_l(:,17,i) = sums_wspts_ws_l(:,i)      ! w*pt* from advec_s_ws
     1939             sums_l(:,17,i) = sums_wspts_ws_l(:,i)                             &
     1940                              * heatflux_output_conversion        ! w*pt* from advec_s_ws
    18841941             IF ( ocean ) sums_l(:,66,i) = sums_wssas_ws_l(:,i) ! w*sa*
    1885              IF ( humidity       )  sums_l(:,49,i)  = sums_wsqs_ws_l(:,i) !w*q*
     1942             IF ( humidity       )  sums_l(:,49,i)  = sums_wsqs_ws_l(:,i)      &
     1943                                            * waterflux_output_conversion !w*q*
    18861944             IF ( passive_scalar )  sums_l(:,116,i) = sums_wsss_ws_l(:,i) !w*s*
    18871945          ENDDO
     
    23782436                                 + ( w(k,j,i)   - w(k,j,i-1) ) * ddx           &
    23792437                                                               )               &
    2380                                * rmask(j,i,sr) * rflags_invers(j,i,k+1)
     2438                               * rmask(j,i,sr) * rflags_invers(j,i,k+1)        &
     2439                               * rho_air_zw(k)                                 &
     2440                               * momentumflux_output_conversion(k)
    23812441!
    23822442!--             Momentum flux w"v"
     
    23872447                                 + ( w(k,j,i)   - w(k,j-1,i) ) * ddy           &
    23882448                                                               )               &
    2389                                * rmask(j,i,sr) * rflags_invers(j,i,k+1)
     2449                               * rmask(j,i,sr) * rflags_invers(j,i,k+1)        &
     2450                               * rho_air_zw(k)                                 &
     2451                               * momentumflux_output_conversion(k)
    23902452!
    23912453!--             Heat flux w"pt"
    23922454                s3 = s3 - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )                 &
    23932455                                 * ( pt(k+1,j,i) - pt(k,j,i) )                 &
     2456                                 * rho_air_zw(k)                               &
     2457                                 * heatflux_output_conversion(k)               &
    23942458                                 * ddzu(k+1) * rmask(j,i,sr)                   &
    23952459                                 * rflags_invers(j,i,k+1)
     
    24352499                   s1 = s1 - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )              &
    24362500                                    * ( vpt(k+1,j,i) - vpt(k,j,i) )            &
     2501                                    * rho_air_zw(k)                            &
     2502                                    * heatflux_output_conversion(k)            &
    24372503                                    * ddzu(k+1) * rmask(j,i,sr)                &
    24382504                                    * rflags_invers(j,i,k+1)
    24392505                   s2 = s2 - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )              &
    24402506                                    * ( q(k+1,j,i) - q(k,j,i) )                &
     2507                                    * rho_air_zw(k)                            &
     2508                                    * waterflux_output_conversion(k)           &
    24412509                                    * ddzu(k+1) * rmask(j,i,sr)                &
    24422510                                    * rflags_invers(j,i,k+1)
     
    24592527                                       *  ( ( q(k+1,j,i) - ql(k+1,j,i) )       &
    24602528                                          - ( q(k,j,i) - ql(k,j,i) ) )         &
     2529                                       * rho_air_zw(k)                         &
     2530                                       * waterflux_output_conversion(k)        &
    24612531                                       * ddzu(k+1) * rmask(j,i,sr)             &
    24622532                                       * rflags_invers(j,i,k+1)
     
    25062576!
    25072577!--             Subgridscale fluxes in the Prandtl layer
    2508                 s1 = s1 + usws(j,i) * rmask(j,i,sr)     ! w"u"
    2509                 s2 = s2 + vsws(j,i) * rmask(j,i,sr)     ! w"v"
    2510                 s3 = s3 + shf(j,i)  * rmask(j,i,sr)     ! w"pt"
     2578                s1 = s1 + usws(j,i) * momentumflux_output_conversion(nzb)      &
     2579                                    * rmask(j,i,sr) ! w"u"
     2580                s2 = s2 + vsws(j,i) * momentumflux_output_conversion(nzb)      &
     2581                                    * rmask(j,i,sr) ! w"v"
     2582                s3 = s3 + shf(j,i)  * heatflux_output_conversion(nzb)          &
     2583                                    * rmask(j,i,sr) ! w"pt"
    25112584                s4 = s4 + 0.0_wp * rmask(j,i,sr)        ! u"pt"
    25122585                s5 = s5 + 0.0_wp * rmask(j,i,sr)        ! v"pt"
     
    25452618             DO  i = nxl, nxr
    25462619                DO  j =  nys, nyn
    2547                    s1 = s1 + qsws(j,i) * rmask(j,i,sr)  ! w"q" (w"qv")
     2620                   s1 = s1 + qsws(j,i) * waterflux_output_conversion(nzb)      &
     2621                                       * rmask(j,i,sr) ! w"q" (w"qv")
    25482622                   s2 = s2 + ( ( 1.0_wp + 0.61_wp * q(nzb,j,i) ) * shf(j,i)    &
    2549                                + 0.61_wp * pt(nzb,j,i) * qsws(j,i) )
     2623                               + 0.61_wp * pt(nzb,j,i) * qsws(j,i) )           &
     2624                             * heatflux_output_conversion(nzb)
    25502625                ENDDO
    25512626             ENDDO
     
    25642639                      s1 = s1 + ( ( 1.0_wp +                                   &
    25652640                                    0.61_wp * q(nzb,j,i) - ql(nzb,j,i) ) *     &
    2566                                   shf(j,i) + 0.61_wp * pt(nzb,j,i) * qsws(j,i) )
     2641                                 shf(j,i) + 0.61_wp * pt(nzb,j,i) * qsws(j,i) )&
     2642                                * heatflux_output_conversion(nzb)
    25672643                   ENDDO
    25682644                ENDDO
     
    25822658!
    25832659!--                   Formula does not work if ql(nzb) /= 0.0
    2584                       s1 = s1 + qsws(j,i) * rmask(j,i,sr)   ! w"q" (w"qv")
     2660                      s1 = s1 + qsws(j,i) * waterflux_output_conversion(nzb)   &
     2661                                          * rmask(j,i,sr)   ! w"q" (w"qv")
    25852662                   ENDDO
    25862663                ENDDO
     
    26242701          DO  i = nxl, nxr
    26252702             DO  j =  nys, nyn
    2626                 s1 = s1 + uswst(j,i) * rmask(j,i,sr)    ! w"u"
    2627                 s2 = s2 + vswst(j,i) * rmask(j,i,sr)    ! w"v"
    2628                 s3 = s3 + tswst(j,i)  * rmask(j,i,sr)   ! w"pt"
     2703                s1 = s1 + uswst(j,i) * momentumflux_output_conversion(nzt:nzt+1) &
     2704                                     * rmask(j,i,sr)    ! w"u"
     2705                s2 = s2 + vswst(j,i) * momentumflux_output_conversion(nzt:nzt+1) &
     2706                                     * rmask(j,i,sr)    ! w"v"
     2707                s3 = s3 + tswst(j,i) * heatflux_output_conversion(nzt:nzt+1)   &
     2708                                     * rmask(j,i,sr)    ! w"pt"
    26292709                s4 = s4 + 0.0_wp * rmask(j,i,sr)        ! u"pt"
    26302710                s5 = s5 + 0.0_wp * rmask(j,i,sr)        ! v"pt"
     
    26632743             DO  i = nxl, nxr
    26642744                DO  j =  nys, nyn
    2665                    s1 = s1 + qswst(j,i) * rmask(j,i,sr) ! w"q" (w"qv")
     2745                   s1 = s1 + qswst(j,i) * waterflux_output_conversion(nzt)     &
     2746                                        * rmask(j,i,sr) ! w"q" (w"qv")
    26662747                   s2 = s2 + ( ( 1.0_wp + 0.61_wp * q(nzt,j,i) ) * tswst(j,i) +&
    2667                                  0.61_wp * pt(nzt,j,i) * qswst(j,i) )
     2748                                 0.61_wp * pt(nzt,j,i) * qswst(j,i) )          &
     2749                             * heatflux_output_conversion(nzt)
    26682750                ENDDO
    26692751             ENDDO
     
    26832765                                    0.61_wp * q(nzt,j,i) - ql(nzt,j,i) ) *     &
    26842766                                  tswst(j,i) +                                 &
    2685                                   0.61_wp * pt(nzt,j,i) * qswst(j,i) )
     2767                                  0.61_wp * pt(nzt,j,i) * qswst(j,i) )         &
     2768                                * heatflux_output_conversion(nzt)
    26862769                   ENDDO
    26872770                ENDDO
     
    27012784!
    27022785!--                   Formula does not work if ql(nzb) /= 0.0
    2703                       s1 = s1 + qswst(j,i) * rmask(j,i,sr)  ! w"q" (w"qv")
     2786                      s1 = s1 + qswst(j,i) * waterflux_output_conversion(nzt)  &
     2787                                           * rmask(j,i,sr)  ! w"q" (w"qv")
    27042788                   ENDDO
    27052789                ENDDO
     
    27552839!--             Energy flux w*e* (has to be adjusted?)
    27562840                s3 = s3 + w(k,j,i) * 0.5_wp * ( ust**2 + vst**2 + w(k,j,i)**2 )&
    2757                                    * rmask(j,i,sr) * rflags_invers(j,i,k+1)
     2841                                   * rmask(j,i,sr) * rflags_invers(j,i,k+1)    &
     2842                                   * momentumflux_output_conversion(k)
    27582843             ENDDO
    27592844          ENDDO
     
    28222907                      s1 = s1 + 0.5_wp * ( vpt(k,j,i)   - hom(k,1,44,sr) +     &
    28232908                                           vpt(k+1,j,i) - hom(k+1,1,44,sr) ) * &
     2909                                         heatflux_output_conversion(k) *       &
    28242910                                         w(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
    28252911                   ENDDO
     
    28392925                         s1 = s1 + 0.5_wp * ( ( q(k,j,i)   - ql(k,j,i)   ) - hom(k,1,42,sr) +   &
    28402926                                              ( q(k+1,j,i) - ql(k+1,j,i) ) - hom(k+1,1,42,sr) ) &
     2927                                          * waterflux_output_conversion(k)                      &
    28412928                                          * w(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
    28422929                      ENDDO
     
    29283015                         s1 = s1 + 0.5_wp * ( vpt(k,j,i)   - hom(k,1,44,sr) +   &
    29293016                                              vpt(k+1,j,i) - hom(k+1,1,44,sr) ) &
     3017                                          * heatflux_output_conversion(k)       &
    29303018                                          * w(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
    29313019                      ENDDO
     
    29393027                !$acc parallel loop present( hom, sums_l )
    29403028                DO  k = nzb, nzt_diff
    2941                    sums_l(k,46,tn) = ( 1.0_wp + 0.61_wp * hom(k,1,41,sr) ) * sums_l(k,17,tn) + &
    2942                                                 0.61_wp * hom(k,1,4,sr) * sums_l(k,49,tn)
     3029                   sums_l(k,46,tn) = ( ( 1.0_wp + 0.61_wp * hom(k,1,41,sr) ) * &
     3030                                       sums_l(k,17,tn) + 0.61_wp *             &
     3031                                       hom(k,1,4,sr) * sums_l(k,49,tn)         &
     3032                                     ) * heatflux_output_conversion(k)
    29433033                ENDDO
    29443034                !$acc end parallel loop
     
    29933083                   s1 = s1 + 0.5_wp * ( w(k,j,i-1) + w(k,j,i) )                &
    29943084                                    * ust * rmask(j,i,sr)                      &
     3085                                    * momentumflux_output_conversion(k)        &
    29953086                                    * rflags_invers(j,i,k+1)
    29963087!
     
    29983089                   s2 = s2 + 0.5_wp * ( w(k,j-1,i) + w(k,j,i) )                &
    29993090                                    * vst * rmask(j,i,sr)                      &
     3091                                    * momentumflux_output_conversion(k)        &
    30003092                                    * rflags_invers(j,i,k+1)
    30013093                ENDDO
    30023094             ENDDO
    30033095             sums_l(k,13,tn) = s1
    3004              sums_l(k,15,tn) = s1
     3096             sums_l(k,15,tn) = s2
    30053097          ENDDO
    30063098          !$acc end parallel loop
     
    30213113                   s1 = s1 + 0.5_wp * ( pt(k,j,i)   - hom(k,1,4,sr) +          &
    30223114                                        pt(k+1,j,i) - hom(k+1,1,4,sr) )        &
     3115                                    * heatflux_output_conversion(k)            &
    30233116                                    * w(k,j,i) * rmask(j,i,sr)                 &
    30243117                                    * rflags_invers(j,i,k+1)
     
    30393132                      s1 = s1 + 0.5_wp * ( q(k,j,i)   - hom(k,1,41,sr) +       &
    30403133                                           q(k+1,j,i) - hom(k+1,1,41,sr) )     &
     3134                                       * waterflux_output_conversion(k)        &
    30413135                                       * w(k,j,i) * rmask(j,i,sr)              &
    30423136                                       * rflags_invers(j,i,k+1)
     
    31673261                                                   ( kh(k,j,i) + kh(k,j,i-1) ) &
    31683262                                                 * ( pt(k,j,i-1) - pt(k,j,i) ) &
     3263                                               * rho_air_zw(k)                 &
     3264                                               * heatflux_output_conversion(k) &
    31693265                                                 * ddx * rmask(j,i,sr)
    31703266                   sums_l(k,61,tn) = sums_l(k,61,tn) - 0.5_wp *                &
    31713267                                                   ( kh(k,j,i) + kh(k,j-1,i) ) &
    31723268                                                 * ( pt(k,j-1,i) - pt(k,j,i) ) &
     3269                                               * rho_air_zw(k)                 &
     3270                                               * heatflux_output_conversion(k) &
    31733271                                                 * ddy * rmask(j,i,sr)
    31743272!
     
    31773275                                     ( u(k,j,i) - hom(k,1,1,sr) ) * 0.5_wp *   &
    31783276                                     ( pt(k,j,i-1) - hom(k,1,4,sr) +           &
    3179                                        pt(k,j,i)   - hom(k,1,4,sr) )
     3277                                       pt(k,j,i)   - hom(k,1,4,sr) )           &
     3278                                     * heatflux_output_conversion(k)
    31803279                   pts = 0.5_wp * ( pt(k,j-1,i) - hom(k,1,4,sr) +              &
    31813280                                    pt(k,j,i)   - hom(k,1,4,sr) )
     
    31833282                                     ( v(k,j,i) - hom(k,1,2,sr) ) * 0.5_wp *   &
    31843283                                     ( pt(k,j-1,i) - hom(k,1,4,sr) +           &
    3185                                        pt(k,j,i)   - hom(k,1,4,sr) )
     3284                                       pt(k,j,i)   - hom(k,1,4,sr) )           &
     3285                                     * heatflux_output_conversion(k)
    31863286                ENDDO
    31873287             ENDDO
     
    34683568       END IF
    34693569
     3570       hom(:,1,121,sr) = rho_air       ! rho_air in Kg/m^3
     3571       hom(:,1,122,sr) = rho_air_zw    ! rho_air_zw in Kg/m^3
     3572
    34703573       hom(:,1,pr_palm,sr) =   sums(:,pr_palm)
    34713574                                       ! u*, w'u', w'v', t* (in last profile)
     
    35763679       IF ( hom(nzb,1,18,sr) > 1.0E-8_wp  .AND.  z_i(1) /= 0.0_wp )  THEN
    35773680          hom(nzb+8,1,pr_palm,sr) = &
    3578              ( g / hom(k_surface_level+1,1,4,sr) * hom(k_surface_level,1,18,sr)&
     3681             ( g / hom(k_surface_level+1,1,4,sr) *                             &
     3682             ( hom(k_surface_level,1,18,sr) / heatflux_output_conversion(nzb) )&
    35793683             * ABS( z_i(1) - mean_surface_level_height(sr) ) )**0.333333333_wp
    35803684       ELSE
  • palm/trunk/SOURCE/header.f90

    r2001 r2037  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! Anelastic approximation implemented
    2323!
    2424! Former revisions:
     
    552552!-- Numerical schemes
    553553    WRITE ( io, 110 )
     554    WRITE ( io, 121 )  TRIM( approximation )
    554555    IF ( psolver(1:7) == 'poisfft' )  THEN
    555556       WRITE ( io, 111 )  TRIM( fft_method )
     
    19061907110 FORMAT (/' Numerical Schemes:'/ &
    19071908             ' -----------------'/)
     1909121 FORMAT (' --> Use the ',A,' approximation for the model equations.')
    19081910111 FORMAT (' --> Solve perturbation pressure via FFT using ',A,' routines')
    19091911112 FORMAT (' --> Solve perturbation pressure via SOR-Red/Black-Schema'/ &
  • palm/trunk/SOURCE/init_3d_model.f90

    r2032 r2037  
    2020! Current revisions:
    2121! ------------------
    22 !
     22! Anelastic approximation implemented
    2323!
    2424! Former revisions:
     
    315315    USE arrays_3d
    316316
     317    USE cloud_parameters,                                                      &
     318        ONLY:  cp, l_v, r_d
     319
    317320    USE constants,                                                             &
    318321        ONLY:  pi
     
    324327   
    325328    USE grid_variables,                                                        &
    326         ONLY:  dx, dy
     329        ONLY:  dx, dy, ddx2_mg, ddy2_mg
    327330   
    328331    USE indices
     
    393396    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  ngp_2dh_outer_l    !<
    394397    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  ngp_2dh_s_inner_l  !<
     398
     399    REAL(wp)     ::  t_surface !< air temperature at the surface
     400
     401    REAL(wp), DIMENSION(:), ALLOCATABLE ::  p_hydrostatic !< hydrostatic pressure
     402
     403    INTEGER(iwp) ::  l       !< loop variable
     404    INTEGER(iwp) ::  nzt_l   !< index of top PE boundary for multigrid level
     405    REAL(wp) ::  dx_l !< grid spacing along x on different multigrid level
     406    REAL(wp) ::  dy_l !< grid spacing along y on different multigrid level
    395407
    396408    REAL(wp), DIMENSION(1:3) ::  volume_flow_area_l     !<
     
    647659
    648660!
     661!-- Allocation of anelastic and Boussinesq approximation specific arrays
     662    ALLOCATE( p_hydrostatic(nzb:nzt+1) )
     663    ALLOCATE( rho_air(nzb:nzt+1) )
     664    ALLOCATE( rho_air_zw(nzb:nzt+1) )
     665    ALLOCATE( drho_air(nzb:nzt+1) )
     666    ALLOCATE( drho_air_zw(nzb:nzt+1) )
     667
     668!
     669!-- Density profile calculation for anelastic approximation
     670    IF ( TRIM( approximation ) == 'anelastic' ) THEN
     671       t_surface = pt_surface * ( surface_pressure / 1000.0_wp )**( r_d / cp )
     672       DO  k = nzb, nzt+1
     673          p_hydrostatic(k)    = surface_pressure * 100.0_wp *                  &
     674                                ( 1 - ( g * zu(k) ) / ( cp * t_surface )       &
     675                                )**( cp / r_d )
     676          rho_air(k)          = ( p_hydrostatic(k) *                           &
     677                                  ( 100000.0_wp / p_hydrostatic(k)             &
     678                                  )**( r_d / cp )                              &
     679                                ) / ( r_d * pt_init(k) )
     680       ENDDO
     681       DO  k = nzb, nzt
     682          rho_air_zw(k) = 0.5_wp * ( rho_air(k) + rho_air(k+1) )
     683       ENDDO
     684       rho_air_zw(nzt+1)  = rho_air_zw(nzt)                                    &
     685                            + 2.0_wp * ( rho_air(nzt+1) - rho_air_zw(nzt)  )
     686    ELSE
     687       rho_air     = 1.0_wp
     688       rho_air_zw  = 1.0_wp
     689    ENDIF
     690
     691!-- compute the inverse density array in order to avoid expencive divisions
     692    drho_air    = 1.0_wp / rho_air
     693    drho_air_zw = 1.0_wp / rho_air_zw
     694
     695!
     696!-- Allocation of flux conversion arrays
     697    ALLOCATE( heatflux_input_conversion(nzb:nzt+1) )
     698    ALLOCATE( waterflux_input_conversion(nzb:nzt+1) )
     699    ALLOCATE( momentumflux_input_conversion(nzb:nzt+1) )
     700    ALLOCATE( heatflux_output_conversion(nzb:nzt+1) )
     701    ALLOCATE( waterflux_output_conversion(nzb:nzt+1) )
     702    ALLOCATE( momentumflux_output_conversion(nzb:nzt+1) )
     703
     704!
     705!-- calculate flux conversion factors according to approximation and in-/output mode
     706    DO  k = nzb, nzt+1
     707
     708        IF ( TRIM( flux_input_mode ) == 'kinematic' )  THEN
     709            heatflux_input_conversion(k)      = rho_air_zw(k)
     710            waterflux_input_conversion(k)     = rho_air_zw(k)
     711            momentumflux_input_conversion(k)  = rho_air_zw(k)
     712        ELSEIF ( TRIM( flux_input_mode ) == 'dynamic' ) THEN
     713            heatflux_input_conversion(k)      = 1.0_wp / cp
     714            waterflux_input_conversion(k)     = 1.0_wp / l_v
     715            momentumflux_input_conversion(k)  = 1.0_wp
     716        ENDIF
     717
     718        IF ( TRIM( flux_output_mode ) == 'kinematic' )  THEN
     719            heatflux_output_conversion(k)     = drho_air_zw(k)
     720            waterflux_output_conversion(k)    = drho_air_zw(k)
     721            momentumflux_output_conversion(k) = drho_air_zw(k)
     722        ELSEIF ( TRIM( flux_output_mode ) == 'dynamic' ) THEN
     723            heatflux_output_conversion(k)     = cp
     724            waterflux_output_conversion(k)    = l_v
     725            momentumflux_output_conversion(k) = 1.0_wp
     726        ENDIF
     727
     728        IF ( .NOT. humidity ) THEN
     729            waterflux_input_conversion(k)  = 1.0_wp
     730            waterflux_output_conversion(k) = 1.0_wp
     731        ENDIF
     732
     733    ENDDO
     734
     735!
     736!-- In case of multigrid method, compute grid lengths and grid factors for the
     737!-- grid levels with respective density on each grid
     738    IF ( psolver(1:9) == 'multigrid' )  THEN
     739
     740       ALLOCATE( ddx2_mg(maximum_grid_level) )
     741       ALLOCATE( ddy2_mg(maximum_grid_level) )
     742       ALLOCATE( dzu_mg(nzb+1:nzt+1,maximum_grid_level) )
     743       ALLOCATE( dzw_mg(nzb+1:nzt+1,maximum_grid_level) )
     744       ALLOCATE( f1_mg(nzb+1:nzt,maximum_grid_level) )
     745       ALLOCATE( f2_mg(nzb+1:nzt,maximum_grid_level) )
     746       ALLOCATE( f3_mg(nzb+1:nzt,maximum_grid_level) )
     747       ALLOCATE( rho_air_mg(nzb:nzt+1,maximum_grid_level) )
     748       ALLOCATE( rho_air_zw_mg(nzb:nzt+1,maximum_grid_level) )
     749
     750       dzu_mg(:,maximum_grid_level) = dzu
     751       rho_air_mg(:,maximum_grid_level) = rho_air
     752!       
     753!--    Next line to ensure an equally spaced grid.
     754       dzu_mg(1,maximum_grid_level) = dzu(2)
     755       rho_air_mg(nzb,maximum_grid_level) = rho_air(nzb) +                     &
     756                                             (rho_air(nzb) - rho_air(nzb+1))
     757
     758       dzw_mg(:,maximum_grid_level) = dzw
     759       rho_air_zw_mg(:,maximum_grid_level) = rho_air_zw
     760       nzt_l = nzt
     761       DO  l = maximum_grid_level-1, 1, -1
     762           dzu_mg(nzb+1,l) = 2.0_wp * dzu_mg(nzb+1,l+1)
     763           dzw_mg(nzb+1,l) = 2.0_wp * dzw_mg(nzb+1,l+1)
     764           rho_air_mg(nzb,l)    = rho_air_mg(nzb,l+1) + (rho_air_mg(nzb,l+1) - rho_air_mg(nzb+1,l+1))
     765           rho_air_zw_mg(nzb,l) = rho_air_zw_mg(nzb,l+1) + (rho_air_zw_mg(nzb,l+1) - rho_air_zw_mg(nzb+1,l+1))
     766           rho_air_mg(nzb+1,l)    = rho_air_mg(nzb+1,l+1)
     767           rho_air_zw_mg(nzb+1,l) = rho_air_zw_mg(nzb+1,l+1)
     768           nzt_l = nzt_l / 2
     769           DO  k = 2, nzt_l+1
     770              dzu_mg(k,l) = dzu_mg(2*k-2,l+1) + dzu_mg(2*k-1,l+1)
     771              dzw_mg(k,l) = dzw_mg(2*k-2,l+1) + dzw_mg(2*k-1,l+1)
     772              rho_air_mg(k,l)    = rho_air_mg(2*k-1,l+1)
     773              rho_air_zw_mg(k,l) = rho_air_zw_mg(2*k-1,l+1)
     774           ENDDO
     775       ENDDO
     776
     777       nzt_l = nzt
     778       dx_l  = dx
     779       dy_l  = dy
     780       DO  l = maximum_grid_level, 1, -1
     781          ddx2_mg(l) = 1.0_wp / dx_l**2
     782          ddy2_mg(l) = 1.0_wp / dy_l**2
     783          DO  k = nzb+1, nzt_l
     784             f2_mg(k,l) = rho_air_zw_mg(k,l) / ( dzu_mg(k+1,l) * dzw_mg(k,l) )
     785             f3_mg(k,l) = rho_air_zw_mg(k-1,l) / ( dzu_mg(k,l)   * dzw_mg(k,l) )
     786             f1_mg(k,l) = 2.0_wp * ( ddx2_mg(l) + ddy2_mg(l) ) &
     787                          * rho_air_mg(k,l) + f2_mg(k,l) + f3_mg(k,l)
     788          ENDDO
     789          nzt_l = nzt_l / 2
     790          dx_l  = dx_l * 2.0_wp
     791          dy_l  = dy_l * 2.0_wp
     792       ENDDO
     793
     794    ENDIF
     795
     796!
    649797!-- 3D-array for storing the dissipation, needed for calculating the sgs
    650798!-- particle velocities
     
    8831031             vsws = 0.0_wp
    8841032          ENDIF
    885           uswst = top_momentumflux_u
    886           vswst = top_momentumflux_v
     1033          uswst = top_momentumflux_u * momentumflux_input_conversion(nzt+1)
     1034          vswst = top_momentumflux_v * momentumflux_input_conversion(nzt+1)
    8871035
    8881036!
     
    10431191          us    = 1E-30_wp
    10441192          usws  = 0.0_wp
    1045           uswst = top_momentumflux_u
     1193          uswst = top_momentumflux_u * momentumflux_input_conversion(nzt+1)
    10461194          vsws  = 0.0_wp
    1047           vswst = top_momentumflux_v
     1195          vswst = top_momentumflux_v * momentumflux_input_conversion(nzt+1)
    10481196          IF ( humidity )       qs = 0.0_wp
    10491197          IF ( passive_scalar ) ss = 0.0_wp
     
    11651313                CALL disturb_heatflux
    11661314             ELSE
    1167                 shf = surface_heatflux
     1315                shf = surface_heatflux * heatflux_input_conversion(nzb)
    11681316!
    11691317!--             Initialize shf with data from external file LSF_DATA
     
    11781326                      DO  j = nysg, nyng
    11791327                         IF ( nzb_s_inner(j,i) /= 0 )  THEN
    1180                             shf(j,i) = wall_heatflux(0)
     1328                            shf(j,i) = wall_heatflux(0)                        &
     1329                                  * heatflux_input_conversion(nzb_s_inner(j,i))
    11811330                         ENDIF
    11821331                      ENDDO
     
    11941343             ENDIF
    11951344             IF ( constant_waterflux )  THEN
    1196                 qsws   = surface_waterflux
     1345                qsws   = surface_waterflux * waterflux_input_conversion(nzb)
    11971346!
    11981347!--             Over topography surface_waterflux is replaced by
     
    12031352                      DO  j = nysg, nyng
    12041353                         IF ( nzb_s_inner(j,i) /= 0 )  THEN
    1205                             qsws(j,i) = wall_qflux(0)
     1354                            qsws(j,i) = wall_qflux(0)                          &
     1355                                 * waterflux_input_conversion(nzb_s_inner(j,i))
    12061356                         ENDIF
    12071357                      ENDDO
     
    12411391!
    12421392!--       Prescribe to heat flux
    1243           IF ( constant_top_heatflux )  tswst = top_heatflux
     1393          IF ( constant_top_heatflux )  tswst = top_heatflux                   &
     1394                                             * heatflux_input_conversion(nzt+1)
    12441395!
    12451396!--       Prescribe zero latent flux at the top     
  • palm/trunk/SOURCE/init_grid.f90

    r2022 r2037  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! Anelastic approximation implemented
    2323!
    2424! Former revisions:
     
    222222         
    223223    USE grid_variables,                                                        &
    224         ONLY:  ddx, ddx2, ddx2_mg, ddy, ddy2, ddy2_mg, dx, dx2, dy, dy2, fwxm, &
     224        ONLY:  ddx, ddx2, ddy, ddy2, dx, dx2, dy, dy2, fwxm,                  &
    225225               fwxp, fwym, fwyp, fxm, fxp, fym, fyp, wall_e_x, wall_e_y,       &
    226226               wall_u, wall_v, wall_w_x, wall_w_y, zu_s_inner, zw_w_inner
     
    293293
    294294    REAL(wp) ::  dum           !< dummy variable to skip columns while reading topography file   
    295     REAL(wp) ::  dx_l          !< grid spacing along x on different multigrid level
    296     REAL(wp) ::  dy_l          !< grid spacing along y on different multigrid level
    297295    REAL(wp) ::  dz_stretched  !< stretched vertical grid spacing
    298296
     
    431429       ddzu_pres = ddzu
    432430       ddzu_pres(1) = ddzu_pres(2)  ! change for lowest level
    433     ENDIF   
    434 
    435 !
    436 !-- In case of multigrid method, compute grid lengths and grid factors for the
    437 !-- grid levels
    438     IF ( psolver(1:9) == 'multigrid' )  THEN
    439 
    440        ALLOCATE( ddx2_mg(maximum_grid_level), ddy2_mg(maximum_grid_level), &
    441                  dzu_mg(nzb+1:nzt+1,maximum_grid_level),                   &
    442                  dzw_mg(nzb+1:nzt+1,maximum_grid_level),                   &
    443                  f1_mg(nzb+1:nzt,maximum_grid_level),                      &
    444                  f2_mg(nzb+1:nzt,maximum_grid_level),                      &
    445                  f3_mg(nzb+1:nzt,maximum_grid_level) )
    446 
    447        dzu_mg(:,maximum_grid_level) = dzu
    448 !       
    449 !--    Next line to ensure an equally spaced grid.
    450        dzu_mg(1,maximum_grid_level) = dzu(2)
    451 
    452        dzw_mg(:,maximum_grid_level) = dzw
    453        nzt_l = nzt
    454        DO  l = maximum_grid_level-1, 1, -1
    455            dzu_mg(nzb+1,l) = 2.0_wp * dzu_mg(nzb+1,l+1)
    456            dzw_mg(nzb+1,l) = 2.0_wp * dzw_mg(nzb+1,l+1)
    457            nzt_l = nzt_l / 2
    458            DO  k = 2, nzt_l+1
    459               dzu_mg(k,l) = dzu_mg(2*k-2,l+1) + dzu_mg(2*k-1,l+1)
    460               dzw_mg(k,l) = dzw_mg(2*k-2,l+1) + dzw_mg(2*k-1,l+1)
    461            ENDDO
    462        ENDDO
    463 
    464        nzt_l = nzt
    465        dx_l  = dx
    466        dy_l  = dy
    467        DO  l = maximum_grid_level, 1, -1
    468           ddx2_mg(l) = 1.0_wp / dx_l**2
    469           ddy2_mg(l) = 1.0_wp / dy_l**2
    470           DO  k = nzb+1, nzt_l
    471              f2_mg(k,l) = 1.0_wp / ( dzu_mg(k+1,l) * dzw_mg(k,l) )
    472              f3_mg(k,l) = 1.0_wp / ( dzu_mg(k,l)   * dzw_mg(k,l) )
    473              f1_mg(k,l) = 2.0_wp * ( ddx2_mg(l) + ddy2_mg(l) ) + &
    474                           f2_mg(k,l) + f3_mg(k,l)
    475           ENDDO
    476           nzt_l = nzt_l / 2
    477           dx_l  = dx_l * 2.0_wp
    478           dy_l  = dy_l * 2.0_wp
    479        ENDDO
    480 
    481431    ENDIF
    482432
  • palm/trunk/SOURCE/ls_forcing_mod.f90

    r2001 r2037  
    2020! Current revisions:
    2121! ------------------
    22 !
     22! Anelastic approximation implemented
    2323!
    2424! Former revisions:
     
    109109           ONLY:  p_surf, pt_surf, q_surf, qsws_surf, shf_surf, td_lsa_lpt,    &
    110110                  td_lsa_q, td_sub_lpt, td_sub_q, time_surf, time_vert,        &
     111                  heatflux_input_conversion, waterflux_input_conversion,       &
    111112                  ug_vert, vg_vert, wsubs_vert, zu
    112113
     
    213214       ENDDO
    214215
     216       shf_surf  = shf_surf  * heatflux_input_conversion(nzb)
     217       qsws_surf = qsws_surf * waterflux_input_conversion(nzb)
    215218
    216219       IF ( time_surf(1) > end_time )  THEN
  • palm/trunk/SOURCE/modules.f90

    r2032 r2037  
    2020! Current revisions:
    2121! ------------------
    22 !
     22! Anelastic approximation implemented
    2323!
    2424! Former revisions:
     
    493493    REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE ::  rif_wall, tri
    494494
     495    REAL(wp), DIMENSION(:), ALLOCATABLE :: rho_air     !< air density profile on the uv grid
     496    REAL(wp), DIMENSION(:), ALLOCATABLE :: rho_air_zw  !< air density profile on the w grid
     497    REAL(wp), DIMENSION(:), ALLOCATABLE :: drho_air    !< inverse air density profile on the uv grid
     498    REAL(wp), DIMENSION(:), ALLOCATABLE :: drho_air_zw !< inverse air density profile on the w grid
     499
     500    REAL(wp), DIMENSION(:,:), ALLOCATABLE :: rho_air_mg    !< air density profiles on the uv grid for multigrid
     501    REAL(wp), DIMENSION(:,:), ALLOCATABLE :: rho_air_zw_mg !< air density profiles on the w grid for multigrid
     502
     503    REAL(wp), DIMENSION(:), ALLOCATABLE :: heatflux_input_conversion      !< conversion factor array for heatflux input
     504    REAL(wp), DIMENSION(:), ALLOCATABLE :: waterflux_input_conversion     !< conversion factor array for waterflux input
     505    REAL(wp), DIMENSION(:), ALLOCATABLE :: momentumflux_input_conversion  !< conversion factor array for momentumflux input
     506    REAL(wp), DIMENSION(:), ALLOCATABLE :: heatflux_output_conversion     !< conversion factor array for heatflux output
     507    REAL(wp), DIMENSION(:), ALLOCATABLE :: waterflux_output_conversion    !< conversion factor array for waterflux output
     508    REAL(wp), DIMENSION(:), ALLOCATABLE :: momentumflux_output_conversion !< conversion factor array for momentumflux output
     509
    495510
    496511    SAVE
     
    616631                             psolver = 'poisfft', &
    617632                             scalar_advec = 'ws-scheme'
     633    CHARACTER (LEN=20)   ::  approximation = 'boussinesq'
     634    CHARACTER (LEN=40)   ::  flux_input_mode = 'approximation-specific'
     635    CHARACTER (LEN=40)   ::  flux_output_mode = 'approximation-specific'
    618636    CHARACTER (LEN=20)   ::  bc_e_b = 'neumann', bc_lr = 'cyclic', &
    619637                             bc_ns = 'cyclic', bc_p_b = 'neumann', &
  • palm/trunk/SOURCE/netcdf_interface_mod.f90

    r2032 r2037  
    2020! Current revisions:
    2121! ------------------
    22 !
     22! Anelastic approximation implemented
    2323!
    2424! Former revisions:
     
    252252             ( 'unknown      ', i9 = 1, dots_max-42 ) /)
    253253
     254    CHARACTER (LEN=16) :: heatflux_output_unit     !< unit for heatflux output
     255    CHARACTER (LEN=16) :: waterflux_output_unit    !< unit for waterflux output
     256    CHARACTER (LEN=16) :: momentumflux_output_unit !< unit for momentumflux output
     257
    254258    CHARACTER (LEN=9), DIMENSION(300) ::  dopr_unit = 'unknown'
    255259
     
    368372            id_var_x_fl, id_var_y_fl, id_var_z_fl,  nc_stat,                   &
    369373            netcdf_data_format, netcdf_data_format_string, netcdf_deflate,     &
    370             netcdf_precision, output_for_t0
    371            
    372            
     374            netcdf_precision, output_for_t0, heatflux_output_unit,             &
     375            waterflux_output_unit, momentumflux_output_unit
     376
    373377    SAVE
    374378
  • palm/trunk/SOURCE/parin.f90

    r2036 r2037  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! Anelastic approximation implemented
    2323!
    2424! Former revisions:
     
    313313
    314314
    315     NAMELIST /inipar/  alpha_surface, background_communication, bc_e_b, bc_lr, &
     315    NAMELIST /inipar/  alpha_surface, approximation,                           &
     316                       background_communication, bc_e_b, bc_lr,                &
    316317                       bc_ns, bc_p_b, bc_p_t, bc_pt_b, bc_pt_t, bc_q_b,        &
    317318             bc_q_t,bc_s_b, bc_s_t, bc_sa_t, bc_uv_b, bc_uv_t,                 &
     
    332333             dt, dt_pr_1d, dt_run_control_1d, dx, dy, dz, dz_max,              &
    333334             dz_stretch_factor, dz_stretch_level, end_time_1d,                 &
    334              ensemble_member_nr,                                               &
    335              e_init, e_min, fft_method, galilei_transformation, humidity,      &
     335             ensemble_member_nr, e_init, e_min, fft_method,                    &
     336             flux_input_mode, flux_output_mode,                                &
     337             galilei_transformation, humidity,                                 &
    336338             inflow_damping_height, inflow_damping_width,                      &
    337339             inflow_disturbance_begin, inflow_disturbance_end,                 &
  • palm/trunk/SOURCE/poismg_mod.f90

    r2022 r2037  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! Anelastic approximation implemented (stll error in optimized multigrid)
    2323!
    2424! Former revisions:
     
    262262
    263263       USE arrays_3d,                                                          &
    264            ONLY:  f1_mg, f2_mg, f3_mg
     264           ONLY:  f1_mg, f2_mg, f3_mg, rho_air_mg
    265265
    266266       USE control_parameters,                                                 &
     
    309309                kp1 = k-ind_even_odd
    310310                r(k,j,i) = f_mg(k,j,i)                                         &
    311                       - ddx2_mg(l) *                                           &
     311                      - rho_air_mg(k,l) * ddx2_mg(l) *                         &
    312312                      ( p_mg(k,j,i+1) +  p_mg(k,j,i-1)  )                      &
    313                       - ddy2_mg(l) *                                           &
     313                      - rho_air_mg(k,l) * ddy2_mg(l) *                         &
    314314                      ( p_mg(k,j+1,i) + p_mg(k,j-1,i)  )                       &
    315315                      - f2_mg_b(k,l) * p_mg(kp1,j,i)                           &
     
    322322                kp1 = k+ind_even_odd+1
    323323                r(k,j,i) = f_mg(k,j,i)                                         &
    324                       - ddx2_mg(l) *                                           &
     324                      - rho_air_mg(k,l) * ddx2_mg(l) *                         &
    325325                      ( p_mg(k,j,i+1) +  p_mg(k,j,i-1)  )                      &
    326                       - ddy2_mg(l) *                                           &
     326                      - rho_air_mg(k,l) * ddy2_mg(l) *                         &
    327327                      ( p_mg(k,j+1,i) + p_mg(k,j-1,i)  )                       &
    328328                      - f2_mg_b(k,l) * p_mg(kp1,j,i)                           &
     
    693693
    694694       USE arrays_3d,                                                          &
    695            ONLY:  f1_mg, f2_mg, f3_mg
     695           ONLY:  f1_mg, f2_mg, f3_mg, rho_air_mg
    696696
    697697       USE control_parameters,                                                 &
     
    755755                         kp1 = k-ind_even_odd
    756756                         p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * (               &
    757                                  ddx2_mg(l) *                                  &
     757                                 rho_air_mg(k,l) * ddx2_mg(l) *                &
    758758                               ( p_mg(k,j,i+1) + p_mg(k,j,i-1) )               &
    759                                + ddy2_mg(l) *                                  &
     759                               + rho_air_mg(k,l) * ddy2_mg(l) *                &
    760760                               ( p_mg(k,j+1,i) + p_mg(k,j-1,i) )               &
    761761                               + f2_mg_b(k,l) * p_mg(kp1,j,i)                  &
     
    774774                         kp1 = k-ind_even_odd
    775775                         p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * (               &
    776                                  ddx2_mg(l) *                                  &
     776                                 rho_air_mg(k,l) * ddx2_mg(l) *                &
    777777                               ( p_mg(k,j,i+1) + p_mg(k,j,i-1) )               &
    778                                + ddy2_mg(l) *                                  &
     778                               + rho_air_mg(k,l) * ddy2_mg(l) *                &
    779779                               ( p_mg(k,j+1,i) + p_mg(k,j-1,i) )               &
    780780                               + f2_mg_b(k,l) * p_mg(kp1,j,i)                  &
     
    793793                         kp1 = k+ind_even_odd+1
    794794                         p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * (               &
    795                                  ddx2_mg(l) *                                  &
     795                                 rho_air_mg(k,l) * ddx2_mg(l) *                &
    796796                               ( p_mg(k,j,i+1) + p_mg(k,j,i-1) )               &
    797                                + ddy2_mg(l) *                                  &
     797                               + rho_air_mg(k,l) * ddy2_mg(l) *                &
    798798                               ( p_mg(k,j+1,i) + p_mg(k,j-1,i) )               &
    799799                               + f2_mg_b(k,l) * p_mg(kp1,j,i)                  &
     
    812812                         kp1 = k+ind_even_odd+1
    813813                         p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * (               &
    814                                  ddx2_mg(l) *                                  &
     814                                 rho_air_mg(k,l) * ddx2_mg(l) *                &
    815815                               ( p_mg(k,j,i+1) + p_mg(k,j,i-1) )               &
    816                                + ddy2_mg(l) *                                  &
     816                               + rho_air_mg(k,l) * ddy2_mg(l) *                &
    817817                               ( p_mg(k,j+1,i) + p_mg(k,j-1,i) )               &
    818818                               + f2_mg_b(k,l) * p_mg(kp1,j,i)                  &
     
    843843                         j   = jj
    844844                         p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * (               &
    845                                  ddx2_mg(l) *                                  &
     845                                 rho_air_mg(k,l) * ddx2_mg(l) *                &
    846846                               ( p_mg(k,j,i+1) + p_mg(k,j,i-1) )               &
    847                                + ddy2_mg(l) *                                  &
     847                               + rho_air_mg(k,l) * ddy2_mg(l) *                &
    848848                               ( p_mg(k,j+1,i) + p_mg(k,j-1,i) )               &
    849849                               + f2_mg_b(k,l) * p_mg(kp1,j,i)                  &
     
    852852                         j = jj+2
    853853                         p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * (               &
    854                                  ddx2_mg(l) *                                  &
     854                                 rho_air_mg(k,l) * ddx2_mg(l) *                &
    855855                               ( p_mg(k,j,i+1) + p_mg(k,j,i-1) )               &
    856                                + ddy2_mg(l) *                                  &
     856                               + rho_air_mg(k,l) * ddy2_mg(l) *                &
    857857                               ( p_mg(k,j+1,i) + p_mg(k,j-1,i) )               &
    858858                               + f2_mg_b(k,l) * p_mg(kp1,j,i)                  &
     
    869869                         j   = jj
    870870                         p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * (               &
    871                                  ddx2_mg(l) *                                  &
     871                                 rho_air_mg(k,l) * ddx2_mg(l) *                &
    872872                               ( p_mg(k,j,i+1) + p_mg(k,j,i-1) )               &
    873                                + ddy2_mg(l) *                                  &
     873                               + rho_air_mg(k,l) * ddy2_mg(l) *                &
    874874                               ( p_mg(k,j+1,i) + p_mg(k,j-1,i) )               &
    875875                               + f2_mg_b(k,l) * p_mg(kp1,j,i)                  &
     
    878878                         j = jj+2
    879879                         p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * (               &
    880                                  ddx2_mg(l) *                                  &
     880                                 rho_air_mg(k,l) * ddx2_mg(l) *                &
    881881                               ( p_mg(k,j,i+1) + p_mg(k,j,i-1) )               &
    882                                + ddy2_mg(l) *                                  &
     882                               + rho_air_mg(k,l) * ddy2_mg(l) *                &
    883883                               ( p_mg(k,j+1,i) + p_mg(k,j-1,i) )               &
    884884                               + f2_mg_b(k,l) * p_mg(kp1,j,i)                  &
     
    895895                         j   = jj
    896896                         p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * (               &
    897                                  ddx2_mg(l) *                                  &
     897                                 rho_air_mg(k,l) * ddx2_mg(l) *                &
    898898                               ( p_mg(k,j,i+1) + p_mg(k,j,i-1) )               &
    899                                + ddy2_mg(l) *                                  &
     899                               + rho_air_mg(k,l) * ddy2_mg(l) *                &
    900900                               ( p_mg(k,j+1,i) + p_mg(k,j-1,i) )               &
    901901                               + f2_mg_b(k,l) * p_mg(kp1,j,i)                  &
     
    904904                         j = jj+2
    905905                         p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * (               &
    906                                  ddx2_mg(l) *                                  &
     906                                 rho_air_mg(k,l) * ddx2_mg(l) *                &
    907907                               ( p_mg(k,j,i+1) + p_mg(k,j,i-1) )               &
    908                                + ddy2_mg(l) *                                  &
     908                               + rho_air_mg(k,l) * ddy2_mg(l) *                &
    909909                               ( p_mg(k,j+1,i) + p_mg(k,j-1,i) )               &
    910910                               + f2_mg_b(k,l) * p_mg(kp1,j,i)                  &
     
    921921                         j   = jj
    922922                         p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * (               &
    923                                  ddx2_mg(l) *                                  &
     923                                 rho_air_mg(k,l) * ddx2_mg(l) *                &
    924924                               ( p_mg(k,j,i+1) + p_mg(k,j,i-1) )               &
    925                                + ddy2_mg(l) *                                  &
     925                               + rho_air_mg(k,l) * ddy2_mg(l) *                &
    926926                               ( p_mg(k,j+1,i) + p_mg(k,j-1,i) )               &
    927927                               + f2_mg_b(k,l) * p_mg(kp1,j,i)                  &
     
    930930                         j = jj+2
    931931                         p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * (               &
    932                                  ddx2_mg(l) *                                  &
     932                                 rho_air_mg(k,l) * ddx2_mg(l) *                &
    933933                               ( p_mg(k,j,i+1) + p_mg(k,j,i-1) )               &
    934                                + ddy2_mg(l) *                                  &
     934                               + rho_air_mg(k,l) * ddy2_mg(l) *                &
    935935                               ( p_mg(k,j+1,i) + p_mg(k,j-1,i) )               &
    936936                               + f2_mg_b(k,l) * p_mg(kp1,j,i)                  &
  • palm/trunk/SOURCE/poismg_noopt.f90

    r2022 r2037  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! Anelastic approximation implemented
    2323!
    2424! Former revisions:
     
    268268
    269269    USE arrays_3d,                                                             &
    270         ONLY:  f1_mg, f2_mg, f3_mg
     270        ONLY:  f1_mg, f2_mg, f3_mg, rho_air_mg
    271271
    272272    USE control_parameters,                                                    &
     
    339339          DO  k = nzb+1, nzt_mg(l)
    340340             r(k,j,i) = f_mg(k,j,i)                                         &
    341                         - ddx2_mg(l) *                                      &
     341                        - rho_air_mg(k,l) * ddx2_mg(l) *                    &
    342342                            ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
    343343                                          ( p_mg(k,j,i) - p_mg(k,j,i+1) ) + &
    344344                              p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * &
    345345                                          ( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) &
    346                         - ddy2_mg(l) *                                      &
     346                        - rho_air_mg(k,l) * ddy2_mg(l) *                    &
    347347                            ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
    348348                                          ( p_mg(k,j,i) - p_mg(k,j+1,i) ) + &
     
    745745
    746746    USE arrays_3d,                                                             &
    747         ONLY:  f1_mg, f2_mg, f3_mg
     747        ONLY:  f1_mg, f2_mg, f3_mg, rho_air_mg
    748748
    749749    USE control_parameters,                                                    &
     
    846846
    847847                      p_mg(k,j,i) = 1.0_wp / f1_mg(k,l) * (                    &
    848                              ddx2_mg(l) *                                      &
     848                             rho_air_mg(k,l) * ddx2_mg(l) *                    &
    849849                               ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
    850850                                             ( p_mg(k,j,i) - p_mg(k,j,i+1) ) + &
    851851                                 p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * &
    852852                                             ( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) &
    853                            + ddy2_mg(l) *                                      &
     853                           + rho_air_mg(k,l) * ddy2_mg(l) *                    &
    854854                               ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
    855855                                             ( p_mg(k,j,i) - p_mg(k,j+1,i) ) + &
     
    869869                   DO  k = nzb+1, nzt_mg(l), 2
    870870                      p_mg(k,j,i) = 1.0_wp / f1_mg(k,l) * (                    &
    871                              ddx2_mg(l) *                                      &
     871                             rho_air_mg(k,l) * ddx2_mg(l) *                    &
    872872                               ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
    873873                                             ( p_mg(k,j,i) - p_mg(k,j,i+1) ) + &
    874874                                 p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * &
    875875                                             ( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) &
    876                            + ddy2_mg(l) *                                      &
     876                           + rho_air_mg(k,l) * ddy2_mg(l) *                    &
    877877                               ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
    878878                                             ( p_mg(k,j,i) - p_mg(k,j+1,i) ) + &
     
    892892                   DO  k = nzb+2, nzt_mg(l), 2
    893893                      p_mg(k,j,i) = 1.0_wp / f1_mg(k,l) * (                    &
    894                              ddx2_mg(l) *                                      &
     894                             rho_air_mg(k,l) * ddx2_mg(l) *                    &
    895895                               ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
    896896                                             ( p_mg(k,j,i) - p_mg(k,j,i+1) ) + &
    897897                                 p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * &
    898898                                             ( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) &
    899                            + ddy2_mg(l) *                                      &
     899                           + rho_air_mg(k,l) * ddy2_mg(l) *                    &
    900900                               ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
    901901                                             ( p_mg(k,j,i) - p_mg(k,j+1,i) ) + &
     
    915915                   DO  k = nzb+2, nzt_mg(l), 2
    916916                      p_mg(k,j,i) = 1.0_wp / f1_mg(k,l) * (                    &
    917                              ddx2_mg(l) *                                      &
     917                             rho_air_mg(k,l) * ddx2_mg(l) *                    &
    918918                               ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
    919919                                             ( p_mg(k,j,i) - p_mg(k,j,i+1) ) + &
    920920                                 p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * &
    921921                                             ( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) &
    922                            + ddy2_mg(l) *                                      &
     922                           + rho_air_mg(k,l) * ddy2_mg(l) *                    &
    923923                               ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
    924924                                             ( p_mg(k,j,i) - p_mg(k,j+1,i) ) + &
     
    947947                      j = jj
    948948                      p_mg(k,j,i) = 1.0_wp / f1_mg(k,l) * (                    &
    949                              ddx2_mg(l) *                                      &
     949                             rho_air_mg(k,l) * ddx2_mg(l) *                    &
    950950                               ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
    951951                                             ( p_mg(k,j,i) - p_mg(k,j,i+1) ) + &
    952952                                 p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * &
    953953                                             ( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) &
    954                            + ddy2_mg(l) *                                      &
     954                           + rho_air_mg(k,l) * ddy2_mg(l) *                    &
    955955                               ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
    956956                                             ( p_mg(k,j,i) - p_mg(k,j+1,i) ) + &
     
    964964                      j = jj+2
    965965                      p_mg(k,j,i) = 1.0_wp / f1_mg(k,l) * (                    &
    966                              ddx2_mg(l) *                                      &
     966                             rho_air_mg(k,l) * ddx2_mg(l) *                    &
    967967                               ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
    968968                                             ( p_mg(k,j,i) - p_mg(k,j,i+1) ) + &
    969969                                 p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * &
    970970                                             ( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) &
    971                            + ddy2_mg(l) *                                      &
     971                           + rho_air_mg(k,l) * ddy2_mg(l) *                    &
    972972                               ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
    973973                                             ( p_mg(k,j,i) - p_mg(k,j+1,i) ) + &
     
    986986                      j =jj
    987987                      p_mg(k,j,i) = 1.0_wp / f1_mg(k,l) * (                    &
    988                              ddx2_mg(l) *                                      &
     988                             rho_air_mg(k,l) * ddx2_mg(l) *                    &
    989989                               ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
    990990                                             ( p_mg(k,j,i) - p_mg(k,j,i+1) ) + &
    991991                                 p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * &
    992992                                             ( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) &
    993                            + ddy2_mg(l) *                                      &
     993                           + rho_air_mg(k,l) * ddy2_mg(l) *                    &
    994994                               ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
    995995                                             ( p_mg(k,j,i) - p_mg(k,j+1,i) ) + &
     
    10031003                      j = jj+2
    10041004                      p_mg(k,j,i) = 1.0_wp / f1_mg(k,l) * (                    &
    1005                              ddx2_mg(l) *                                      &
     1005                             rho_air_mg(k,l) * ddx2_mg(l) *                    &
    10061006                               ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
    10071007                                             ( p_mg(k,j,i) - p_mg(k,j,i+1) ) + &
    10081008                                 p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * &
    10091009                                             ( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) &
    1010                            + ddy2_mg(l) *                                      &
     1010                           + rho_air_mg(k,l) * ddy2_mg(l) *                    &
    10111011                               ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
    10121012                                             ( p_mg(k,j,i) - p_mg(k,j+1,i) ) + &
     
    10251025                      j =jj
    10261026                      p_mg(k,j,i) = 1.0_wp / f1_mg(k,l) * (                    &
    1027                              ddx2_mg(l) *                                      &
     1027                             rho_air_mg(k,l) * ddx2_mg(l) *                    &
    10281028                               ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
    10291029                                             ( p_mg(k,j,i) - p_mg(k,j,i+1) ) + &
    10301030                                 p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * &
    10311031                                             ( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) &
    1032                            + ddy2_mg(l) *                                      &
     1032                           + rho_air_mg(k,l) * ddy2_mg(l) *                    &
    10331033                               ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
    10341034                                             ( p_mg(k,j,i) - p_mg(k,j+1,i) ) + &
     
    10421042                      j = jj+2
    10431043                      p_mg(k,j,i) = 1.0_wp / f1_mg(k,l) * (                    &
    1044                              ddx2_mg(l) *                                      &
     1044                             rho_air_mg(k,l) * ddx2_mg(l) *                    &
    10451045                               ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
    10461046                                             ( p_mg(k,j,i) - p_mg(k,j,i+1) ) + &
    10471047                                 p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * &
    10481048                                             ( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) &
    1049                            + ddy2_mg(l) *                                      &
     1049                           + rho_air_mg(k,l) * ddy2_mg(l) *                    &
    10501050                               ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
    10511051                                             ( p_mg(k,j,i) - p_mg(k,j+1,i) ) + &
     
    10641064                      j =jj
    10651065                      p_mg(k,j,i) = 1.0_wp / f1_mg(k,l) * (                    &
    1066                              ddx2_mg(l) *                                      &
     1066                             rho_air_mg(k,l) * ddx2_mg(l) *                    &
    10671067                               ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
    10681068                                             ( p_mg(k,j,i) - p_mg(k,j,i+1) ) + &
    10691069                                 p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * &
    10701070                                             ( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) &
    1071                            + ddy2_mg(l) *                                      &
     1071                           + rho_air_mg(k,l) * ddy2_mg(l) *                    &
    10721072                               ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
    10731073                                             ( p_mg(k,j,i) - p_mg(k,j+1,i) ) + &
     
    10811081                      j = jj+2
    10821082                      p_mg(k,j,i) = 1.0_wp / f1_mg(k,l) * (                    &
    1083                              ddx2_mg(l) *                                      &
     1083                             rho_air_mg(k,l) * ddx2_mg(l) *                    &
    10841084                               ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
    10851085                                             ( p_mg(k,j,i) - p_mg(k,j,i+1) ) + &
    10861086                                 p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * &
    10871087                                             ( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) &
    1088                            + ddy2_mg(l) *                                      &
     1088                           + rho_air_mg(k,l) * ddy2_mg(l) *                    &
    10891089                               ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
    10901090                                             ( p_mg(k,j,i) - p_mg(k,j+1,i) ) + &
  • palm/trunk/SOURCE/pres.f90

    r2001 r2037  
    2020! Current revisions:
    2121! ------------------
    22 !
     22! Anelastic approximation implemented
    2323!
    2424! Former revisions:
     
    128128
    129129    USE arrays_3d,                                                             &
    130         ONLY:  d, ddzu, ddzu_pres, ddzw, dzw, p, p_loc, tend, u, v, w
     130        ONLY:  d, ddzu, ddzu_pres, ddzw, dzw, p, p_loc, rho_air, rho_air_zw,   &
     131               tend, u, v, w
    131132
    132133    USE control_parameters,                                                    &
     
    399400       DO  j = nys, nyn
    400401          DO  k = nzb_s_inner(j,i)+1, nzt
    401              d(k,j,i) = ( ( u(k,j,i+1) - u(k,j,i) ) * ddx + &
    402                           ( v(k,j+1,i) - v(k,j,i) ) * ddy + &
    403                           ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) ) * ddt_3d      &
    404                                                                 * d_weight_pres
     402             d(k,j,i) = ( ( u(k,j,i+1) - u(k,j,i) ) * rho_air(k) * ddx +       &
     403                          ( v(k,j+1,i) - v(k,j,i) ) * rho_air(k) * ddy +       &
     404                          ( w(k,j,i)   * rho_air_zw(k) -                       &
     405                            w(k-1,j,i) * rho_air_zw(k-1) ) * ddzw(k)           &
     406                        ) * ddt_3d * d_weight_pres
    405407          ENDDO
    406408!
     
    427429       DO  j = nys, nyn
    428430          DO  k = 1, nzt
    429              d(k,j,i) = ( ( u(k,j,i+1) - u(k,j,i) ) * ddx +             &
    430                         ( v(k,j+1,i) - v(k,j,i) ) * ddy +               &
    431                         ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) ) * ddt_3d  &
    432                         * d_weight_pres * rflags_s_inner(k,j,i)
     431             d(k,j,i) = ( ( u(k,j,i+1) - u(k,j,i) ) * rho_air(k) * ddx +       &
     432                          ( v(k,j+1,i) - v(k,j,i) ) * rho_air(k) * ddy +       &
     433                          ( w(k,j,i)   * rho_air_zw(k) -                       &
     434                            w(k-1,j,i) * rho_air_zw(k-1) ) * ddzw(k)           &
     435                        ) * ddt_3d * d_weight_pres * rflags_s_inner(k,j,i)
    433436          ENDDO
    434437       ENDDO
     
    808811          DO  j = nys, nyn
    809812             DO  k = nzb_s_inner(j,i)+1, nzt
    810                 d(k,j,i) = ( u(k,j,i+1) - u(k,j,i) ) * ddx +                   &
    811                            ( v(k,j+1,i) - v(k,j,i) ) * ddy +                   &
    812                            ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
     813             d(k,j,i) = ( u(k,j,i+1) - u(k,j,i) ) * rho_air(k) * ddx +         &
     814                        ( v(k,j+1,i) - v(k,j,i) ) * rho_air(k) * ddy +         &
     815                        ( w(k,j,i)   * rho_air_zw(k) -                         &
     816                          w(k-1,j,i) * rho_air_zw(k-1) ) * ddzw(k)
    813817             ENDDO
    814818             DO  k = nzb+1, nzt
     
    823827          DO  j = nys, nyn
    824828             DO  k = 1, nzt
    825                    d(k,j,i) = ( ( u(k,j,i+1) - u(k,j,i) ) * ddx +              &
    826                                 ( v(k,j+1,i) - v(k,j,i) ) * ddy +              &
    827                                 ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)            &
    828                               ) * rflags_s_inner(k,j,i)
     829                d(k,j,i) = ( ( u(k,j,i+1) - u(k,j,i) ) * rho_air(k) * ddx +    &
     830                             ( v(k,j+1,i) - v(k,j,i) ) * rho_air(k) * ddy +    &
     831                             ( w(k,j,i)   * rho_air_zw(k) -                    &
     832                               w(k-1,j,i) * rho_air_zw(k-1) ) * ddzw(k)        &
     833                           ) * rflags_s_inner(k,j,i)
    829834             ENDDO
    830835          ENDDO
  • palm/trunk/SOURCE/sor.f90

    r2001 r2037  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! Anelastic approximation implemented
    2323!
    2424! Former revisions:
     
    5959!------------------------------------------------------------------------------!
    6060 SUBROUTINE sor( d, ddzu, ddzw, p )
    61  
     61
     62    USE arrays_3d,                                                             &
     63        ONLY:  rho_air, rho_air_zw
    6264
    6365    USE grid_variables,                                                        &
     
    101103!-- Compute pre-factors.
    102104    DO  k = 1, nz
    103          f2(k) = ddzu(k+1) * ddzw(k)
    104          f3(k) = ddzu(k)   * ddzw(k)
    105          f1(k) = 2.0_wp * ( ddx2 + ddy2 ) + f2(k) + f3(k)
     105         f2(k) = ddzu(k+1) * ddzw(k) * rho_air_zw(k)
     106         f3(k) = ddzu(k)   * ddzw(k) * rho_air_zw(k-1)
     107         f1(k) = 2.0_wp * ( ddx2 + ddy2 ) * rho_air(k) + f2(k) + f3(k)
    106108    ENDDO
    107109
     
    131133             DO  k = nzb+1, nzt
    132134                p(k,j,i) = p(k,j,i) + omega_sor / f1(k) * (            &
    133                                ddx2 * ( p(k,j,i+1) + p(k,j,i-1) ) +    &
    134                                ddy2 * ( p(k,j+1,i) + p(k,j-1,i) ) +    &
    135                                f2(k) * p(k+1,j,i)                 +    &
    136                                f3(k) * p(k-1,j,i)                 -    &
    137                                d(k,j,i)                           -    &
    138                                f1(k) * p(k,j,i)           )
     135                           rho_air(k) * ddx2 * ( p(k,j,i+1) + p(k,j,i-1) ) +   &
     136                           rho_air(k) * ddy2 * ( p(k,j+1,i) + p(k,j-1,i) ) +   &
     137                           f2(k) * p(k+1,j,i)                              +   &
     138                           f3(k) * p(k-1,j,i)                              -   &
     139                           d(k,j,i)                                        -   &
     140                           f1(k) * p(k,j,i)           )
    139141             ENDDO
    140142          ENDDO
     
    144146          DO  j = nys1, nyn, 2
    145147             DO  k = nzb+1, nzt
    146                 p(k,j,i) = p(k,j,i) + omega_sor / f1(k) * (            &
    147                                ddx2 * ( p(k,j,i+1) + p(k,j,i-1) ) +    &
    148                                ddy2 * ( p(k,j+1,i) + p(k,j-1,i) ) +    &
    149                                f2(k) * p(k+1,j,i)                 +    &
    150                                f3(k) * p(k-1,j,i)                 -    &
    151                                d(k,j,i)                           -    &
    152                                f1(k) * p(k,j,i)           )
     148                p(k,j,i) = p(k,j,i) + omega_sor / f1(k) * (                    &
     149                           rho_air(k) * ddx2 * ( p(k,j,i+1) + p(k,j,i-1) ) +   &
     150                           rho_air(k) * ddy2 * ( p(k,j+1,i) + p(k,j-1,i) ) +   &
     151                           f2(k) * p(k+1,j,i)                              +   &
     152                           f3(k) * p(k-1,j,i)                              -   &
     153                           d(k,j,i)                                        -   &
     154                           f1(k) * p(k,j,i)           )
    153155             ENDDO
    154156          ENDDO
     
    176178             DO  k = nzb+1, nzt
    177179                p(k,j,i) = p(k,j,i) + omega_sor / f1(k) * (            &
    178                                ddx2 * ( p(k,j,i+1) + p(k,j,i-1) ) +    &
    179                                ddy2 * ( p(k,j+1,i) + p(k,j-1,i) ) +    &
    180                                f2(k) * p(k+1,j,i)                 +    &
    181                                f3(k) * p(k-1,j,i)                 -    &
    182                                d(k,j,i)                           -    &
    183                                f1(k) * p(k,j,i)           )
     180                           rho_air(k) * ddx2 * ( p(k,j,i+1) + p(k,j,i-1) ) +   &
     181                           rho_air(k) * ddy2 * ( p(k,j+1,i) + p(k,j-1,i) ) +   &
     182                           f2(k) * p(k+1,j,i)                              +   &
     183                           f3(k) * p(k-1,j,i)                              -   &
     184                           d(k,j,i)                                        -   &
     185                           f1(k) * p(k,j,i)           )
    184186             ENDDO
    185187          ENDDO
     
    190192             DO  k = nzb+1, nzt
    191193                p(k,j,i) = p(k,j,i) + omega_sor / f1(k) * (            &
    192                                ddx2 * ( p(k,j,i+1) + p(k,j,i-1) ) +    &
    193                                ddy2 * ( p(k,j+1,i) + p(k,j-1,i) ) +    &
    194                                f2(k) * p(k+1,j,i)                 +    &
    195                                f3(k) * p(k-1,j,i)                 -    &
    196                                d(k,j,i)                           -    &
    197                                f1(k) * p(k,j,i)           )
     194                           rho_air(k) * ddx2 * ( p(k,j,i+1) + p(k,j,i-1) ) +   &
     195                           rho_air(k) * ddy2 * ( p(k,j+1,i) + p(k,j-1,i) ) +   &
     196                           f2(k) * p(k+1,j,i)                              +   &
     197                           f3(k) * p(k-1,j,i)                              -   &
     198                           d(k,j,i)                                        -   &
     199                           f1(k) * p(k,j,i)           )
    198200             ENDDO
    199201          ENDDO
  • palm/trunk/SOURCE/surface_layer_fluxes_mod.f90

    r2012 r2037  
    2020! Current revisions:
    2121! ------------------
    22 !
     22! Anelastic approximation implemented
    2323!
    2424! Former revisions:
     
    158158        ONLY:  e, kh, nr, nrs, nrsws, ol, pt, q, ql, qr, qrs, qrsws, qs, qsws, &
    159159               s, shf, ss, ssws, ts, u, us, usws, v, vpt, vsws, zu, zw, z0,    &
    160                z0h, z0q
     160               z0h, z0q, drho_air_zw, rho_air_zw
    161161
    162162    USE cloud_parameters,                                                      &
     
    540540                      rib(j,i) = - g * z_mo * ( ( 1.0_wp + 0.61_wp             &
    541541                                 * q(k+1,j,i) ) * shf(j,i) + 0.61_wp           &
    542                                  * pt(k+1,j,i) * qsws(j,i) )   &
     542                                 * pt(k+1,j,i) * qsws(j,i) ) * drho_air_zw(k)  &
    543543                                 / ( uv_total(j,i)**3 * vpt(k+1,j,i) * kappa**2&
    544544                                 + 1.0E-20_wp)
    545545                   ELSE
    546                       rib(j,i) = - g * z_mo * shf(j,i)                         &
     546                      rib(j,i) = - g * z_mo * shf(j,i) * drho_air_zw(k)        &
    547547                           / ( uv_total(j,i)**3 * pt(k+1,j,i) * kappa**2       &
    548548                           + 1.0E-20_wp )
     
    828828!--       For a given heat flux in the surface layer:
    829829          !$OMP PARALLEL DO
    830           !$acc kernels loop private( j )
     830          !$acc kernels loop private( j, k )
    831831          DO  i = nxlg, nxrg
    832832             DO  j = nysg, nyng
    833                 ts(j,i) = -shf(j,i) / ( us(j,i) + 1E-30_wp )
     833                k   = nzb_s_inner(j,i)
     834                ts(j,i) = -shf(j,i) * drho_air_zw(k) / ( us(j,i) + 1E-30_wp )
    834835!
    835836!--             ts must be limited, because otherwise overflow may occur in case
     
    888889             DO  i = nxlg, nxrg
    889890                DO  j = nysg, nyng
    890                    qs(j,i) = -qsws(j,i) / ( us(j,i) + 1E-30_wp )
     891                   k   = nzb_s_inner(j,i)
     892                   qs(j,i) = -qsws(j,i) * drho_air_zw(k) / ( us(j,i) + 1E-30_wp )
    891893                ENDDO
    892894             ENDDO
     
    10251027                                     + psi_m( z0(j,i) / ol_mid ) )
    10261028
    1027              usws(j,i) = -usws(j,i) * 0.5_wp * ( us(j,i-1) + us(j,i) )
     1029             usws(j,i) = -usws(j,i) * 0.5_wp * ( us(j,i-1) + us(j,i) )         &
     1030                                    * rho_air_zw(k)
    10281031          ENDDO
    10291032       ENDDO
     
    10521055                                     + psi_m( z0(j,i) / ol_mid ) )
    10531056
    1054              vsws(j,i) = -vsws(j,i) * 0.5_wp * ( us(j,i-1) + us(j,i) )
     1057             vsws(j,i) = -vsws(j,i) * 0.5_wp * ( us(j,i-1) + us(j,i) )         &
     1058                                    * rho_air_zw(k)
    10551059
    10561060          ENDDO
     
    10741078             !$acc loop independent
    10751079             DO  j = nysg, nyng
    1076                 shf(j,i) = -ts(j,i) * us(j,i)
     1080                k   = nzb_s_inner(j,i)
     1081                shf(j,i) = -ts(j,i) * us(j,i) * rho_air_zw(k)
    10771082             ENDDO
    10781083          ENDDO
     
    10901095             !$acc loop independent
    10911096             DO  j = nysg, nyng
    1092                 qsws(j,i) = -qs(j,i) * us(j,i)
     1097                k   = nzb_s_inner(j,i)
     1098                qsws(j,i) = -qs(j,i) * us(j,i) * rho_air_zw(k)
    10931099             ENDDO
    10941100          ENDDO
  • palm/trunk/SOURCE/tridia_solver_mod.f90

    r2001 r2037  
    2020! Current revisions:
    2121! ------------------
    22 !
     22! Anelastic approximation implemented
    2323!
    2424! Former revisions:
     
    129129
    130130       USE arrays_3d,                                                          &
    131            ONLY:  ddzu_pres, ddzw
     131           ONLY:  ddzu_pres, ddzw, rho_air_zw
    132132
    133133       USE kinds
     
    140140
    141141       DO  k = 0, nz-1
    142           ddzuw(k,1) = ddzu_pres(k+1) * ddzw(k+1)
    143           ddzuw(k,2) = ddzu_pres(k+2) * ddzw(k+1)
     142          ddzuw(k,1) = ddzu_pres(k+1) * ddzw(k+1) * rho_air_zw(k)
     143          ddzuw(k,2) = ddzu_pres(k+2) * ddzw(k+1) * rho_air_zw(k+1)
    144144          ddzuw(k,3) = -1.0_wp * &
    145                        ( ddzu_pres(k+2) * ddzw(k+1) + ddzu_pres(k+1) * ddzw(k+1) )
     145                       ( ddzu_pres(k+2) * ddzw(k+1) * rho_air_zw(k+1) +        &
     146                         ddzu_pres(k+1) * ddzw(k+1) * rho_air_zw(k) )
    146147       ENDDO
    147148!
     
    168169
    169170          USE arrays_3d,                                                       &
    170               ONLY:  tric
     171              ONLY:  tric, rho_air
    171172
    172173          USE constants,                                                       &
     
    231232             DO  j = nys_z, nyn_z
    232233                DO  i = nxl_z, nxr_z
    233                    tric(i,j,k) = ddzuw(k,3) - ll(i,j)
     234                   tric(i,j,k) = ddzuw(k,3) - ll(i,j) * rho_air(k+1)
    234235                ENDDO
    235236             ENDDO
     
    491492
    492493       USE arrays_3d,                                                          &
    493            ONLY:  ddzu_pres, ddzw
     494           ONLY:  ddzu_pres, ddzw, rho_air, rho_air_zw
    494495
    495496       USE control_parameters,                                                 &
     
    526527       DO  k = 0, nz-1
    527528          DO  i = 0,nx
    528              tri_for_1d(2,i,k) = ddzu_pres(k+1) * ddzw(k+1)
    529              tri_for_1d(3,i,k) = ddzu_pres(k+2) * ddzw(k+1)
     529             tri_for_1d(2,i,k) = ddzu_pres(k+1) * ddzw(k+1) * rho_air_zw(k)
     530             tri_for_1d(3,i,k) = ddzu_pres(k+2) * ddzw(k+1) * rho_air_zw(k+1)
    530531          ENDDO
    531532       ENDDO
     
    591592          DO  k = 0, nz-1
    592593             DO  i = 0, nx
    593                 a = -1.0_wp * ddzu_pres(k+2) * ddzw(k+1)
    594                 c = -1.0_wp * ddzu_pres(k+1) * ddzw(k+1)
    595                 tri_for_1d(1,i,k) = a + c - l(i)
     594                a = -1.0_wp * ddzu_pres(k+2) * ddzw(k+1) * rho_air_zw(k+1)
     595                c = -1.0_wp * ddzu_pres(k+1) * ddzw(k+1) * rho_air_zw(k)
     596                tri_for_1d(1,i,k) = a + c - l(i) * rho_air(k+1)
    596597             ENDDO
    597598          ENDDO
Note: See TracChangeset for help on using the changeset viewer.