Index: /palm/trunk/SOURCE/advec_ws.f90
===================================================================
--- /palm/trunk/SOURCE/advec_ws.f90	(revision 2036)
+++ /palm/trunk/SOURCE/advec_ws.f90	(revision 2037)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! ------------------
-! 
+! Anelastic approximation implemented
 ! 
 ! Former revisions:
@@ -902,5 +902,5 @@
 
        USE arrays_3d,                                                          &
-           ONLY:  ddzw, tend, u, v, w
+           ONLY:  ddzw, drho_air, tend, u, v, w, rho_air, rho_air_zw
 
        USE constants,                                                          &
@@ -1211,5 +1211,5 @@
 
 
-          flux_t(k) = w(k,j,i) * (                                            &
+          flux_t(k) = w(k,j,i) * rho_air_zw(k) * (                            &
                      ( 37.0_wp * ibit8 * adv_sca_5                            &
                   +     7.0_wp * ibit7 * adv_sca_3                            &
@@ -1225,5 +1225,5 @@
                                  )
 
-          diss_t(k) = -ABS( w(k,j,i) ) * (                                    &
+          diss_t(k) = -ABS( w(k,j,i) ) * rho_air_zw(k) * (                    &
                      ( 10.0_wp * ibit8 * adv_sca_5                            &
                   +     3.0_wp * ibit7 * adv_sca_3                            &
@@ -1266,9 +1266,9 @@
              fd_1   =  ( w(k-1,j,i)   * ( sk(k,j,i) + sk(k-1,j,i) )           &
                    -ABS( w(k-1,j,i) ) * ( sk(k,j,i) - sk(k-1,j,i) )           &
-                       ) * adv_sca_1 
+                       ) * adv_sca_1 * rho_air_zw(k)
 
              ft_1   =  ( w(k,j,i)   * ( sk(k+1,j,i) + sk(k,j,i) )             &
                    -ABS( w(k,j,i) ) * ( sk(k+1,j,i) - sk(k,j,i) )             &
-                      ) * adv_sca_1 
+                      ) * adv_sca_1 * rho_air_zw(k)
 !
 !--          Calculate ratio of upwind gradients. Note, Min/Max is just to 
@@ -1386,5 +1386,5 @@
                                          + IBITS(wall_flags_0(k,j,i-1),2,1)    &
                                          )                                     &
-                          ) * ddx                                              &
+                          ) * rho_air(k) * ddx                                 &
                         + ( v(k,j+1,i) * ( ibit3 + ibit4 + ibit5 )             &
                           - v(k,j,i)   * ( IBITS(wall_flags_0(k,j-1,i),3,1)    &
@@ -1392,7 +1392,9 @@
                                          + IBITS(wall_flags_0(k,j-1,i),5,1)    &
                                          )                                     &
-                          ) * ddy                                              &
-                        + ( w(k,j,i)   * ( ibit6 + ibit7 + ibit8 )             &
-                          - w(k-1,j,i) * ( IBITS(wall_flags_0(k-1,j,i),6,1)    &
+                          ) * rho_air(k) * ddy                                 &
+                        + ( w(k,j,i) * rho_air_zw(k) *                         &
+                                         ( ibit6 + ibit7 + ibit8 )             &
+                          - w(k-1,j,i) * rho_air_zw(k-1) *                     &
+                                         ( IBITS(wall_flags_0(k-1,j,i),6,1)    &
                                          + IBITS(wall_flags_0(k-1,j,i),7,1)    &
                                          + IBITS(wall_flags_0(k-1,j,i),8,1)    &
@@ -1406,6 +1408,7 @@
                       + ( flux_n(k) + diss_n(k) - swap_flux_y_local(k,tn)   - &
                           swap_diss_y_local(k,tn)              ) * ddy        &
-                      + ( flux_t(k) + diss_t(k) - flux_d - diss_d             &
-                                                               ) * ddzw(k)    &
+                      + ( ( flux_t(k) + diss_t(k) ) -                         &
+                          ( flux_d    + diss_d    )                           &
+                                                    ) * drho_air(k) * ddzw(k) &
                                       ) + sk(k,j,i) * div
 
@@ -1454,5 +1457,5 @@
           k_mm  = k - 2 * ibit8
 
-          flux_t(k) = w(k,j,i) * (                                            &
+          flux_t(k) = w(k,j,i) * rho_air_zw(k) * (                            &
                     ( 37.0_wp * ibit8 * adv_sca_5                             &
                  +     7.0_wp * ibit7 * adv_sca_3                             &
@@ -1468,5 +1471,5 @@
                                  )
 
-          diss_t(k) = -ABS( w(k,j,i) ) * (                                    &
+          diss_t(k) = -ABS( w(k,j,i) ) * rho_air_zw(k) * (                    &
                     ( 10.0_wp * ibit8 * adv_sca_5                             &
                  +     3.0_wp * ibit7 * adv_sca_3                             &
@@ -1511,9 +1514,9 @@
              fd_1   =  ( w(k-1,j,i)   * ( sk(k,j,i) + sk(k-1,j,i) )           &
                    -ABS( w(k-1,j,i) ) * ( sk(k,j,i) - sk(k-1,j,i) )           &
-                       ) * adv_sca_1 
+                       ) * adv_sca_1  * rho_air_zw(k)
 
              ft_1   =  ( w(k,j,i)   * ( sk(k+1,j,i) + sk(k,j,i) )             &
                    -ABS( w(k,j,i) ) * ( sk(k+1,j,i) - sk(k,j,i) )             &
-                       ) * adv_sca_1 
+                       ) * adv_sca_1  * rho_air_zw(k)
 !
 !--          Calculate ratio of upwind gradients. Note, Min/Max is just to 
@@ -1626,7 +1629,8 @@
 !--       correction is needed to overcome numerical instabilities introduced
 !--       by a not sufficient reduction of divergences near topography.
-          div         =   ( u(k,j,i+1) - u(k,j,i)   ) * ddx                   &
-                        + ( v(k,j+1,i) - v(k,j,i)   ) * ddy                   &
-                        + ( w(k,j,i)   - w(k-1,j,i) ) * ddzw(k)
+          div         =   ( u(k,j,i+1) - u(k,j,i)   ) * rho_air(k) * ddx      &
+                        + ( v(k,j+1,i) - v(k,j,i)   ) * rho_air(k) * ddy      &
+                        + ( w(k,j,i)   * rho_air_zw(k) -                      &
+                            w(k-1,j,i) * rho_air_zw(k-1) ) * ddzw(k)
 
           tend(k,j,i) = tend(k,j,i) - (                                       &
@@ -1635,6 +1639,7 @@
                       + ( flux_n(k) + diss_n(k) - swap_flux_y_local(k,tn)   - &
                           swap_diss_y_local(k,tn)              ) * ddy        &
-                      + ( flux_t(k) + diss_t(k) - flux_d - diss_d             &
-                                                               ) * ddzw(k)    &
+                      + ( ( flux_t(k) + diss_t(k) ) -                         &
+                          ( flux_d    + diss_d    )                           &
+                                                    ) * drho_air(k) * ddzw(k) &
                                       ) + sk(k,j,i) * div
 
@@ -1734,5 +1739,6 @@
 
        USE arrays_3d,                                                         &
-           ONLY:  ddzw, diss_l_u, diss_s_u, flux_l_u, flux_s_u, tend, u, v, w
+           ONLY:  ddzw, diss_l_u, diss_s_u, flux_l_u, flux_s_u, tend, u, v, w,&
+                  drho_air, rho_air, rho_air_zw
 
        USE constants,                                                         &
@@ -1999,5 +2005,5 @@
 
           w_comp    = w(k,j,i) + w(k,j,i-1)
-          flux_t(k) = w_comp  * (                                             &
+          flux_t(k) = w_comp * rho_air_zw(k) * (                              &
                      ( 37.0_wp * ibit17 * adv_mom_5                           &
                   +     7.0_wp * ibit16 * adv_mom_3                           &
@@ -2014,5 +2020,5 @@
                                  )
 
-          diss_t(k) = - ABS( w_comp ) * (                                     &
+          diss_t(k) = - ABS( w_comp ) * rho_air_zw(k) * (                     &
                      ( 10.0_wp * ibit17 * adv_mom_5                           &
                   +     3.0_wp * ibit16 * adv_mom_3                           &
@@ -2037,6 +2043,6 @@
                                       + IBITS(wall_flags_0(k,j,i-1),10,1)     &
                                       + IBITS(wall_flags_0(k,j,i-1),11,1)     &
-                                      )                                       &   
-                  ) * ddx                                                     &
+                                      )                                       &
+                  ) * rho_air(k) * ddx                                        &
                +  ( ( v_comp + gv ) * ( ibit12 + ibit13 + ibit14 )            &
                   - ( v(k,j,i)   + v(k,j,i-1 )  )                             &
@@ -2045,7 +2051,7 @@
                                       + IBITS(wall_flags_0(k,j-1,i),14,1)     &
                                       )                                       &
-                  ) * ddy                                                     &
-               +  ( w_comp          * ( ibit15 + ibit16 + ibit17 )            &
-                - ( w(k-1,j,i) + w(k-1,j,i-1) )                               &
+                  ) * rho_air(k) * ddy                                        &
+               +  ( w_comp * rho_air_zw(k) * ( ibit15 + ibit16 + ibit17 )     &
+                - ( w(k-1,j,i) + w(k-1,j,i-1) ) * rho_air_zw(k-1)             &
                                     * ( IBITS(wall_flags_0(k-1,j,i),15,1)     &
                                       + IBITS(wall_flags_0(k-1,j,i),16,1)     &
@@ -2061,6 +2067,7 @@
                           + ( flux_n(k) + diss_n(k)                           &
                           -   flux_s_u(k,tn) - diss_s_u(k,tn)     ) * ddy     &
-                          + ( flux_t(k) + diss_t(k)                           &
-                          -   flux_d    - diss_d                  ) * ddzw(k) &
+                          + ( ( flux_t(k) + diss_t(k) )                       &
+                          -   ( flux_d    + diss_d )                          &
+                                                    ) * drho_air(k) * ddzw(k) &
                                        ) + div * u(k,j,i)
 
@@ -2127,5 +2134,5 @@
 
           w_comp    = w(k,j,i) + w(k,j,i-1)
-          flux_t(k) = w_comp  * (                                             &
+          flux_t(k) = w_comp * rho_air_zw(k) * (                              &
                      ( 37.0_wp * ibit17 * adv_mom_5                           &
                   +     7.0_wp * ibit16 * adv_mom_3                           &
@@ -2142,5 +2149,5 @@
                                  )
 
-          diss_t(k) = - ABS( w_comp ) * (                                     &
+          diss_t(k) = - ABS( w_comp ) * rho_air_zw(k) * (                     &
                      ( 10.0_wp * ibit17 * adv_mom_5                           &
                   +     3.0_wp * ibit16 * adv_mom_3                           &
@@ -2161,6 +2168,10 @@
 !--       by a not sufficient reduction of divergences near topography.
           div = ( ( u_comp(k)   - ( u(k,j,i)   + u(k,j,i-1)   ) ) * ddx       &
+                                                                  * rho_air(k)&
                +  ( v_comp + gv - ( v(k,j,i)   + v(k,j,i-1 )  ) ) * ddy       &
-               +  ( w_comp      - ( w(k-1,j,i) + w(k-1,j,i-1) ) ) * ddzw(k)   &
+                                                                  * rho_air(k)&
+               +  (   w_comp                      * rho_air_zw(k)   -         &
+                    ( w(k-1,j,i) + w(k-1,j,i-1) ) * rho_air_zw(k-1)           &
+                  ) * ddzw(k)                                                 &
                 ) * 0.5_wp
 
@@ -2170,6 +2181,7 @@
                           + ( flux_n(k) + diss_n(k)                           &
                           -   flux_s_u(k,tn) - diss_s_u(k,tn)     ) * ddy     &
-                          + ( flux_t(k) + diss_t(k)                           &
-                          -   flux_d    - diss_d                  ) * ddzw(k) &
+                          + ( ( flux_t(k) + diss_t(k) )                       &
+                          -   ( flux_d    + diss_d    )                       &
+                                                    ) * drho_air(k) * ddzw(k) &
                                        ) + div * u(k,j,i)
 
@@ -2219,5 +2231,6 @@
 
        USE arrays_3d,                                                          &
-           ONLY:  ddzw, diss_l_v, diss_s_v, flux_l_v, flux_s_v, tend, u, v, w
+           ONLY:  ddzw, diss_l_v, diss_s_v, flux_l_v, flux_s_v, tend, u, v, w, &
+                  drho_air, rho_air, rho_air_zw
 
        USE constants,                                                          &
@@ -2486,5 +2499,5 @@
 
           w_comp    = w(k,j-1,i) + w(k,j,i)
-          flux_t(k) = w_comp  * (                                             &
+          flux_t(k) = w_comp * rho_air_zw(k) * (                              &
                      ( 37.0_wp * ibit26 * adv_mom_5                           &
                   +     7.0_wp * ibit25 * adv_mom_3                           &
@@ -2501,5 +2514,5 @@
                                  )
 
-          diss_t(k) = - ABS( w_comp ) * (                                     &
+          diss_t(k) = - ABS( w_comp ) * rho_air_zw(k) * (                     &
                      ( 10.0_wp * ibit26 * adv_mom_5                           &
                   +     3.0_wp * ibit25 * adv_mom_3                           &
@@ -2525,6 +2538,6 @@
                                          + IBITS(wall_flags_0(k,j,i-1),19,1)  &
                                          + IBITS(wall_flags_0(k,j,i-1),20,1)  &
-                                         )                                    &   
-                  ) * ddx                                                     &
+                                         )                                    &
+                  ) * rho_air(k) * ddx                                        &
                +  ( v_comp(k)                                                 &
                                        * ( ibit21 + ibit22 + ibit23 )         &
@@ -2533,9 +2546,8 @@
                                          + IBITS(wall_flags_0(k,j-1,i),22,1)  &
                                          + IBITS(wall_flags_0(k,j-1,i),23,1)  &
-                                         )                                    &   
-                  ) * ddy                                                     &
-               +  ( w_comp                                                    &
-                                       * ( ibit24 + ibit25 + ibit26 )         &
-                - ( w(k-1,j-1,i) + w(k-1,j,i) )                               &
+                                         )                                    &
+                  ) * rho_air(k) * ddy                                        &
+               +  ( w_comp * rho_air_zw(k) * ( ibit24 + ibit25 + ibit26 )     &
+                - ( w(k-1,j-1,i) + w(k-1,j,i) ) * rho_air_zw(k-1)             &
                                        * ( IBITS(wall_flags_0(k-1,j,i),24,1)  &
                                          + IBITS(wall_flags_0(k-1,j,i),25,1)  &
@@ -2551,6 +2563,7 @@
                        + ( flux_n(k) + diss_n(k)                              &
                        -   flux_s_v(k,tn) - diss_s_v(k,tn)       ) * ddy      &
-                       + ( flux_t(k) + diss_t(k)                              &
-                       -   flux_d    - diss_d                    ) * ddzw(k)  &
+                       + ( ( flux_t(k) + diss_t(k) )                          &
+                       -   ( flux_d    + diss_d    )                          &
+                                                   ) * drho_air(k) * ddzw(k)  &
                                       ) + v(k,j,i) * div
 
@@ -2622,5 +2635,5 @@
 
           w_comp    = w(k,j-1,i) + w(k,j,i)
-          flux_t(k) = w_comp  * (                                             &
+          flux_t(k) = w_comp * rho_air_zw(k) * (                              &
                      ( 37.0_wp * ibit26 * adv_mom_5                           &
                   +     7.0_wp * ibit25 * adv_mom_3                           &
@@ -2637,5 +2650,5 @@
                                  )
 
-          diss_t(k) = - ABS( w_comp ) * (                                     &
+          diss_t(k) = - ABS( w_comp ) * rho_air_zw(k) * (                     &
                      ( 10.0_wp * ibit26 * adv_mom_5                           &
                   +     3.0_wp * ibit25 * adv_mom_3                           &
@@ -2656,6 +2669,10 @@
 !--       by a not sufficient reduction of divergences near topography.
           div = ( ( u_comp + gu - ( u(k,j-1,i)   + u(k,j,i)   ) ) * ddx       &
+                                                                  * rho_air(k)&
                +  ( v_comp(k)   - ( v(k,j,i)     + v(k,j-1,i) ) ) * ddy       &
-               +  ( w_comp      - ( w(k-1,j-1,i) + w(k-1,j,i) ) ) * ddzw(k)   &
+                                                                  * rho_air(k)&
+               +  (   w_comp                      * rho_air_zw(k)   -         &
+                    ( w(k-1,j-1,i) + w(k-1,j,i) ) * rho_air_zw(k-1)           &
+                  ) * ddzw(k)                                                 &
                 ) * 0.5_wp
 
@@ -2665,6 +2682,7 @@
                        + ( flux_n(k) + diss_n(k)                              &
                        -   flux_s_v(k,tn) - diss_s_v(k,tn)       ) * ddy      &
-                       + ( flux_t(k) + diss_t(k)                              &
-                       -   flux_d    - diss_d                    ) * ddzw(k)  &
+                       + ( ( flux_t(k) + diss_t(k) )                          &
+                       -   ( flux_d    + diss_d    )                          &
+                                                   ) * drho_air(k) * ddzw(k)  &
                                       ) + v(k,j,i) * div
 
@@ -2714,5 +2732,6 @@
 
        USE arrays_3d,                                                         &
-           ONLY:  ddzu, diss_l_w, diss_s_w, flux_l_w, flux_s_w, tend, u, v, w
+           ONLY:  ddzu, diss_l_w, diss_s_w, flux_l_w, flux_s_w, tend, u, v, w,&
+                  drho_air_zw, rho_air, rho_air_zw
 
        USE constants,                                                         &
@@ -2984,5 +3003,5 @@
 
           w_comp    = w(k+1,j,i) + w(k,j,i)
-          flux_t(k) = w_comp  * (                                             &
+          flux_t(k) = w_comp * rho_air(k+1) * (                               &
                      ( 37.0_wp * ibit35 * adv_mom_5                           &
                   +     7.0_wp * ibit34 * adv_mom_3                           &
@@ -2999,5 +3018,5 @@
                                 )
 
-          diss_t(k) = - ABS( w_comp ) * (                                     &
+          diss_t(k) = - ABS( w_comp ) * rho_air(k+1) * (                      &
                      ( 10.0_wp * ibit35 * adv_mom_5                           &
                   +     3.0_wp * ibit34 * adv_mom_3                           &
@@ -3024,5 +3043,5 @@
                                       + IBITS(wall_flags_0(k,j,i-1),29,1)     &
                                       )                                       &
-                  ) * ddx                                                     &
+                  ) * rho_air_zw(k) * ddx                                     &
               +   ( ( v_comp + gv ) * ( ibit30 + ibit31 + ibit32 )            & 
                   - ( v(k+1,j,i) + v(k,j,i)   )                               &
@@ -3031,7 +3050,7 @@
                                       + IBITS(wall_flags_00(k,j-1,i),0,1)     &
                                       )                                       &
-                  ) * ddy                                                     &
-              +   ( w_comp          * ( ibit33 + ibit34 + ibit35 )            &
-                - ( w(k,j,i)   + w(k-1,j,i)   )                               &
+                  ) * rho_air_zw(k) * ddy                                     &
+              +   ( w_comp * rho_air(k+1) * ( ibit33 + ibit34 + ibit35 )      &
+                - ( w(k,j,i)   + w(k-1,j,i)   ) * rho_air(k)                  &
                                     * ( IBITS(wall_flags_00(k-1,j,i),1,1)     &
                                       + IBITS(wall_flags_00(k-1,j,i),2,1)     &
@@ -3047,6 +3066,7 @@
                     + ( flux_n(k) + diss_n(k)                                 &
                     -   flux_s_w(k,tn) - diss_s_w(k,tn)       ) * ddy         &
-                    + ( flux_t(k) + diss_t(k)                                 &
-                    -   flux_d    - diss_d                    ) * ddzu(k+1)   &
+                    + ( ( flux_t(k) + diss_t(k) )                             &
+                    -   ( flux_d    + diss_d    )                             &
+                                              ) * drho_air_zw(k) * ddzu(k+1)  &
                                       ) + div * w(k,j,i)
 
@@ -3105,5 +3125,5 @@
 
           w_comp    = w(k+1,j,i) + w(k,j,i)
-          flux_t(k) = w_comp  * (                                             &
+          flux_t(k) = w_comp * rho_air(k+1) * (                               &
                      ( 37.0_wp * ibit35 * adv_mom_5                           &
                   +     7.0_wp * ibit34 * adv_mom_3                           &
@@ -3120,5 +3140,5 @@
                                 )
 
-          diss_t(k) = - ABS( w_comp ) * (                                     &
+          diss_t(k) = - ABS( w_comp ) * rho_air(k+1) * (                      &
                      ( 10.0_wp * ibit35 * adv_mom_5                           &
                   +     3.0_wp * ibit34 * adv_mom_3                           &
@@ -3139,6 +3159,10 @@
 !--       by a not sufficient reduction of divergences near topography.
           div = ( ( u_comp + gu - ( u(k+1,j,i) + u(k,j,i)   ) ) * ddx         &
+                                                              * rho_air_zw(k) &
               +   ( v_comp + gv - ( v(k+1,j,i) + v(k,j,i)   ) ) * ddy         &
-              +   ( w_comp      - ( w(k,j,i)   + w(k-1,j,i) ) ) * ddzu(k+1)   &
+                                                              * rho_air_zw(k) &
+              +   (   w_comp                    * rho_air(k+1) -              &
+                    ( w(k,j,i)   + w(k-1,j,i) ) * rho_air(k)                  &
+                  ) * ddzu(k+1)                                               &
                 ) * 0.5_wp
 
@@ -3148,6 +3172,7 @@
                     + ( flux_n(k) + diss_n(k)                                 &
                     -   flux_s_w(k,tn) - diss_s_w(k,tn)       ) * ddy         &
-                    + ( flux_t(k) + diss_t(k)                                 &
-                    -   flux_d    - diss_d                    ) * ddzu(k+1)   &
+                    + ( ( flux_t(k) + diss_t(k) )                             &
+                    -   ( flux_d    + diss_d    )                             &
+                                              ) * drho_air_zw(k) * ddzu(k+1)  &
                                       ) + div * w(k,j,i)
 
@@ -3183,5 +3208,5 @@
 
        USE arrays_3d,                                                         &
-           ONLY:  ddzw, tend, u, v, w
+           ONLY:  ddzw, drho_air, tend, u, v, w, rho_air, rho_air_zw
 
        USE constants,                                                         &
@@ -3484,5 +3509,5 @@
 
 
-                flux_t(k) = w(k,j,i) * (                                      &
+                flux_t(k) = w(k,j,i) * rho_air_zw(k) * (                      &
                            ( 37.0_wp * ibit8 * adv_sca_5                      &
                         +     7.0_wp * ibit7 * adv_sca_3                      &
@@ -3498,5 +3523,5 @@
                                        )
 
-                diss_t(k) = -ABS( w(k,j,i) ) * (                              &
+                diss_t(k) = -ABS( w(k,j,i) ) * rho_air_zw(k) * (              &
                            ( 10.0_wp * ibit8 * adv_sca_5                      &
                         +     3.0_wp * ibit7 * adv_sca_3                      &
@@ -3539,9 +3564,9 @@
                    fd_1   = (  w(k-1,j,i)   * ( sk(k,j,i) + sk(k-1,j,i) )     &
                          -ABS( w(k-1,j,i) ) * ( sk(k,j,i) - sk(k-1,j,i) )     &
-                            ) * adv_sca_1 
+                            ) * adv_sca_1 * rho_air_zw(k)
 
                    ft_1   = (  w(k,j,i)   * ( sk(k+1,j,i) + sk(k,j,i) )       &
                          -ABS( w(k,j,i) ) * ( sk(k+1,j,i) - sk(k,j,i) )       &
-                            ) * adv_sca_1 
+                            ) * adv_sca_1 * rho_air_zw(k)
 !
 !--                Calculate ratio of upwind gradients. Note, Min/Max is just
@@ -3659,5 +3684,5 @@
                                          + IBITS(wall_flags_0(k,j,i-1),2,1)    &
                                          )                                     &
-                          ) * ddx                                              &
+                          ) * rho_air(k) * ddx                                 &
                         + ( v(k,j+1,i) * ( ibit3 + ibit4 + ibit5 )             &
                           - v(k,j,i)   * ( IBITS(wall_flags_0(k,j-1,i),3,1)    &
@@ -3665,7 +3690,9 @@
                                          + IBITS(wall_flags_0(k,j-1,i),5,1)    &
                                          )                                     &
-                          ) * ddy                                              &
-                        + ( w(k,j,i)   * ( ibit6 + ibit7 + ibit8 )             &
-                          - w(k-1,j,i) * ( IBITS(wall_flags_0(k-1,j,i),6,1)    &
+                          ) * rho_air(k) * ddy                                 &
+                        + ( w(k,j,i) * rho_air_zw(k) *                         &
+                                         ( ibit6 + ibit7 + ibit8 )             &
+                          - w(k-1,j,i) * rho_air_zw(k-1) *                     &
+                                         ( IBITS(wall_flags_0(k-1,j,i),6,1)    &
                                          + IBITS(wall_flags_0(k-1,j,i),7,1)    &
                                          + IBITS(wall_flags_0(k-1,j,i),8,1)    &
@@ -3679,6 +3706,7 @@
                       + ( flux_n(k) + diss_n(k) - swap_flux_y_local(k)   -    &
                           swap_diss_y_local(k)              ) * ddy           &
-                      + ( flux_t(k) + diss_t(k) - flux_d - diss_d             &
-                                                               ) * ddzw(k)    &
+                      + ( ( flux_t(k) + diss_t(k) ) -                         &
+                          ( flux_d    + diss_d    )                           &
+                                                    ) * drho_air(k) * ddzw(k) &
                                             ) + sk(k,j,i) * div
 
@@ -3725,5 +3753,5 @@
 
 
-                flux_t(k) = w(k,j,i) * (                                      &
+                flux_t(k) = w(k,j,i) * rho_air_zw(k) * (                      &
                            ( 37.0_wp * ibit8 * adv_sca_5                      &
                         +     7.0_wp * ibit7 * adv_sca_3                      &
@@ -3739,5 +3767,5 @@
                                        )
 
-                diss_t(k) = -ABS( w(k,j,i) ) * (                              &
+                diss_t(k) = -ABS( w(k,j,i) ) * rho_air_zw(k) * (              &
                            ( 10.0_wp * ibit8 * adv_sca_5                      &
                         +     3.0_wp * ibit7 * adv_sca_3                      &
@@ -3780,9 +3808,9 @@
                    fd_1   = (  w(k-1,j,i)   * ( sk(k,j,i) + sk(k-1,j,i) )     &
                          -ABS( w(k-1,j,i) ) * ( sk(k,j,i) - sk(k-1,j,i) )     &
-                            ) * adv_sca_1 
+                            ) * adv_sca_1 * rho_air_zw(k)
 
                    ft_1   = (  w(k,j,i)   * ( sk(k+1,j,i) + sk(k,j,i) )       &
                          -ABS( w(k,j,i) ) * ( sk(k+1,j,i) - sk(k,j,i) )       &
-                            ) * adv_sca_1 
+                            ) * adv_sca_1 * rho_air_zw(k)
 !
 !--                Calculate ratio of upwind gradients. Note, Min/Max is just
@@ -3895,7 +3923,8 @@
 !--             correction is needed to overcome numerical instabilities introduced
 !--             by a not sufficient reduction of divergences near topography.
-                div         =   ( u(k,j,i+1) - u(k,j,i)   ) * ddx             &
-                              + ( v(k,j+1,i) - v(k,j,i)   ) * ddy             &
-                              + ( w(k,j,i)   - w(k-1,j,i) ) * ddzw(k)
+                div         =   ( u(k,j,i+1) - u(k,j,i)   ) * rho_air(k) * ddx &
+                              + ( v(k,j+1,i) - v(k,j,i)   ) * rho_air(k) * ddy &
+                              + ( w(k,j,i)   * rho_air_zw(k) -                 &
+                                  w(k-1,j,i) * rho_air_zw(k-1) ) * ddzw(k)
 
                 tend(k,j,i) = tend(k,j,i) - (                                 &
@@ -3904,6 +3933,7 @@
                       + ( flux_n(k) + diss_n(k) - swap_flux_y_local(k)   -    &
                           swap_diss_y_local(k)              ) * ddy           &
-                      + ( flux_t(k) + diss_t(k) - flux_d - diss_d             &
-                                                               ) * ddzw(k)    &
+                      + ( ( flux_t(k) + diss_t(k) ) -                         &
+                          ( flux_d    + diss_d    )                           &
+                                                    ) * drho_air(k) * ddzw(k) &
                                             ) + sk(k,j,i) * div
 
@@ -4004,5 +4034,5 @@
 
        USE arrays_3d,                                                         &
-           ONLY:  ddzw, tend, u, v, w
+           ONLY:  ddzw, drho_air, tend, u, v, w, rho_air, rho_air_zw
 
        USE constants,                                                         &
@@ -4278,5 +4308,5 @@
                 k_mm  = k - 2 * ibit8
 
-                flux_t    = w(k,j,i) * (                                      &
+                flux_t    = w(k,j,i) * rho_air_zw(k) * (                      &
                            ( 37.0_wp * ibit8 * adv_sca_5                      &
                         +     7.0_wp * ibit7 * adv_sca_3                      &
@@ -4292,5 +4322,5 @@
                                        )
 
-                diss_t    = -ABS( w(k,j,i) ) * (                              &
+                diss_t    = -ABS( w(k,j,i) ) * rho_air_zw(k) * (              &
                            ( 10.0_wp * ibit8 * adv_sca_5                      &
                         +     3.0_wp * ibit7 * adv_sca_3                      &
@@ -4333,9 +4363,9 @@
                    fd_1   = (  w(k-1,j,i)   * ( sk(k,j,i) + sk(k-1,j,i) )     &
                         -ABS( w(k-1,j,i) ) * ( sk(k,j,i) - sk(k-1,j,i) )      &
-                            ) * adv_sca_1 
+                            ) * adv_sca_1 * rho_air_zw(k)
 
                    ft_1   = (  w(k,j,i)   * ( sk(k+1,j,i) + sk(k,j,i) )       &
                         -ABS( w(k,j,i) ) * ( sk(k+1,j,i) - sk(k,j,i) )        &
-                            ) * adv_sca_1 
+                            ) * adv_sca_1 * rho_air_zw(k)
 !
 !--                Calculate ratio of upwind gradients. Note, Min/Max is just
@@ -4447,5 +4477,5 @@
                                          + IBITS(wall_flags_0(k,j,i-1),2,1)    &
                                          )                                     &
-                          ) * ddx                                              &
+                          ) * rho_air(k) * ddx                                 &
                         + ( v(k,j+1,i) * ( ibit3 + ibit4 + ibit5 )             &
                           - v(k,j,i)   * ( IBITS(wall_flags_0(k,j-1,i),3,1)    &
@@ -4453,7 +4483,9 @@
                                          + IBITS(wall_flags_0(k,j-1,i),5,1)    &
                                          )                                     &
-                          ) * ddy                                              &
-                        + ( w(k,j,i)   * ( ibit6 + ibit7 + ibit8 )             &
-                          - w(k-1,j,i) * ( IBITS(wall_flags_0(k-1,j,i),6,1)    &
+                          ) * rho_air(k) * ddy                                 &
+                        + ( w(k,j,i) * rho_air_zw(k) *                         &
+                                         ( ibit6 + ibit7 + ibit8 )             &
+                          - w(k-1,j,i) * rho_air_zw(k-1) *                     &
+                                         ( IBITS(wall_flags_0(k-1,j,i),6,1)    &
                                          + IBITS(wall_flags_0(k-1,j,i),7,1)    &
                                          + IBITS(wall_flags_0(k-1,j,i),8,1)    &
@@ -4465,5 +4497,7 @@
                                ( flux_r + diss_r - flux_l - diss_l ) * ddx    &
                              + ( flux_n + diss_n - flux_s - diss_s ) * ddy    &
-                             + ( flux_t + diss_t - flux_d - diss_d ) * ddzw(k)&
+                             + ( ( flux_t + diss_t ) -                        &
+                                 ( flux_d + diss_d )                          &
+                                                    ) * drho_air(k) * ddzw(k) &
                                 ) + div * sk(k,j,i)
 
@@ -4511,5 +4545,5 @@
 
        USE arrays_3d,                                                          &
-           ONLY:  ddzw, tend, u, v, w
+           ONLY:  ddzw, drho_air, tend, u, v, w, rho_air, rho_air_zw
 
        USE constants,                                                          &
@@ -4780,5 +4814,5 @@
 
                 w_comp    = w(k,j,i) + w(k,j,i-1)
-                flux_t(k) = w_comp  * (                                      &
+                flux_t(k) = w_comp * rho_air_zw(k) * (                       &
                           ( 37.0_wp * ibit17 * adv_mom_5                        &
                        +     7.0_wp * ibit16 * adv_mom_3                        &
@@ -4795,5 +4829,5 @@
                                       )
 
-                diss_t(k) = - ABS( w_comp ) * (                              &
+                diss_t(k) = - ABS( w_comp ) * rho_air_zw(k) * (              &
                           ( 10.0_wp * ibit17 * adv_mom_5                        &
                        +     3.0_wp * ibit16 * adv_mom_3                        &
@@ -4818,6 +4852,6 @@
                                       + IBITS(wall_flags_0(k,j,i-1),10,1)     &
                                       + IBITS(wall_flags_0(k,j,i-1),11,1)     &
-                                      )                                       &   
-                  ) * ddx                                                     &
+                                      )                                       &
+                  ) * rho_air(k) * ddx                                        &
                +  ( ( v_comp + gv ) * ( ibit12 + ibit13 + ibit14 )            &
                   - ( v(k,j,i)   + v(k,j,i-1 )  )                             &
@@ -4826,7 +4860,7 @@
                                       + IBITS(wall_flags_0(k,j-1,i),14,1)     &
                                       )                                       &
-                  ) * ddy                                                     &
-               +  ( w_comp          * ( ibit15 + ibit16 + ibit17 )            &
-                - ( w(k-1,j,i) + w(k-1,j,i-1) )                               &
+                  ) * rho_air(k) * ddy                                        &
+               +  ( w_comp * rho_air_zw(k) * ( ibit15 + ibit16 + ibit17 )     &
+                - ( w(k-1,j,i) + w(k-1,j,i-1) ) * rho_air_zw(k-1)             &
                                     * ( IBITS(wall_flags_0(k-1,j,i),15,1)     &
                                       + IBITS(wall_flags_0(k-1,j,i),16,1)     &
@@ -4843,7 +4877,7 @@
                + ( flux_n(k) + diss_n(k)                                       &
                -   swap_flux_y_local_u(k)   - swap_diss_y_local_u(k)   ) * ddy &
-               + ( flux_t(k) + diss_t(k)                                       &
-               -   flux_d    - diss_d                                          &
-                                                                  ) * ddzw(k)  &
+               + ( ( flux_t(k) + diss_t(k) )                                   &
+               -   ( flux_d    + diss_d    )                                   &
+                                                    ) * drho_air(k) * ddzw(k)  &
                                            ) + div * u(k,j,i)
 
@@ -4911,5 +4945,5 @@
 
                 w_comp    = w(k,j,i) + w(k,j,i-1)
-                flux_t(k) = w_comp  * (                                      &
+                flux_t(k) = w_comp * rho_air_zw(k) * (                       &
                           ( 37.0_wp * ibit17 * adv_mom_5                        &
                        +     7.0_wp * ibit16 * adv_mom_3                        &
@@ -4926,5 +4960,5 @@
                                       )
 
-                diss_t(k) = - ABS( w_comp ) * (                              &
+                diss_t(k) = - ABS( w_comp ) * rho_air_zw(k) * (              &
                           ( 10.0_wp * ibit17 * adv_mom_5                        &
                        +     3.0_wp * ibit16 * adv_mom_3                        &
@@ -4945,7 +4979,10 @@
 !--             by a not sufficient reduction of divergences near topography.
                 div = ( ( u_comp(k)   - ( u(k,j,i)   + u(k,j,i-1)   ) ) * ddx &
+                                                                  * rho_air(k)&
                      +  ( v_comp + gv - ( v(k,j,i)   + v(k,j,i-1 )  ) ) * ddy &
-                     +  ( w_comp      - ( w(k-1,j,i) + w(k-1,j,i-1) ) )       &
-                                                                    * ddzw(k) &
+                                                                  * rho_air(k)&
+                     +  (   w_comp                      * rho_air_zw(k) -     &
+                          ( w(k-1,j,i) + w(k-1,j,i-1) ) * rho_air_zw(k-1)     &
+                        ) * ddzw(k)                                           &
                       ) * 0.5_wp
 
@@ -4955,7 +4992,7 @@
                + ( flux_n(k) + diss_n(k)                                       &
                -   swap_flux_y_local_u(k)   - swap_diss_y_local_u(k)   ) * ddy &
-               + ( flux_t(k) + diss_t(k)                                       &
-               -   flux_d    - diss_d                                          &
-                                                                  ) * ddzw(k)  &
+               + ( ( flux_t(k) + diss_t(k) )                                   &
+               -   ( flux_d    + diss_d    )                                   &
+                                                    ) * drho_air(k) * ddzw(k)  &
                                            ) + div * u(k,j,i)
 
@@ -5004,5 +5041,5 @@
 
        USE arrays_3d,                                                          &
-           ONLY:  ddzw, tend, u, v, w
+           ONLY:  ddzw, drho_air, tend, u, v, w, rho_air, rho_air_zw
 
        USE constants,                                                          &
@@ -5267,5 +5304,5 @@
 
                 w_comp    = w(k,j,i) + w(k,j,i-1)
-                flux_t    = w_comp  * (                                      &
+                flux_t    = w_comp * rho_air_zw(k) * (                       &
                           ( 37.0_wp * ibit17 * adv_mom_5                        &
                        +     7.0_wp * ibit16 * adv_mom_3                        &
@@ -5282,5 +5319,5 @@
                                       )
 
-                diss_t    = - ABS( w_comp ) * (                              &
+                diss_t    = - ABS( w_comp ) * rho_air_zw(k) * (              &
                           ( 10.0_wp * ibit17 * adv_mom_5                        &
                        +     3.0_wp * ibit16 * adv_mom_3                        &
@@ -5305,6 +5342,6 @@
                                       + IBITS(wall_flags_0(k,j,i-1),10,1)     &
                                       + IBITS(wall_flags_0(k,j,i-1),11,1)     &
-                                      )                                       &   
-                  ) * ddx                                                     &
+                                      )                                       &
+                  ) * rho_air(k) * ddx                                        &
                +  ( ( v_comp + gv ) * ( ibit12 + ibit13 + ibit14 )            &
                   - ( v(k,j,i)   + v(k,j,i-1 )  )                             &
@@ -5313,11 +5350,11 @@
                                       + IBITS(wall_flags_0(k,j-1,i),14,1)     &
                                       )                                       &
-                  ) * ddy                                                     &
-               +  ( w_comp          * ( ibit15 + ibit16 + ibit17 )            &
-                - ( w(k-1,j,i) + w(k-1,j,i-1) )                               &
+                  ) * rho_air(k) * ddy                                        &
+               +  ( w_comp * rho_air_zw(k) * ( ibit15 + ibit16 + ibit17 )     &
+                - ( w(k-1,j,i) + w(k-1,j,i-1) ) * rho_air_zw(k-1)             &
                                     * ( IBITS(wall_flags_0(k-1,j,i),15,1)     &
                                       + IBITS(wall_flags_0(k-1,j,i),16,1)     &
                                       + IBITS(wall_flags_0(k-1,j,i),17,1)     &
-                                      )                                       &  
+                                      )                                       &
                   ) * ddzw(k)   &
                 ) * 0.5_wp
@@ -5327,5 +5364,7 @@
                                ( flux_r + diss_r - flux_l - diss_l ) * ddx     &
                              + ( flux_n + diss_n - flux_s - diss_s ) * ddy     &
-                             + ( flux_t + diss_t - flux_d - diss_d ) * ddzw(k) &
+                             + ( ( flux_t + diss_t ) -                         &
+                                 ( flux_d + diss_d )                           &
+                                                     ) * drho_air(k) * ddzw(k) &
                                 ) + div * u(k,j,i)
 
@@ -5365,5 +5404,5 @@
 
        USE arrays_3d,                                                          &
-           ONLY:  ddzw, tend, u, v, w
+           ONLY:  ddzw, drho_air, tend, u, v, w, rho_air, rho_air_zw
 
        USE constants,                                                          &
@@ -5633,5 +5672,5 @@
 
                 w_comp    = w(k,j-1,i) + w(k,j,i)
-                flux_t(k) = w_comp  * (                                      &
+                flux_t(k) = w_comp * rho_air_zw(k) * (                       &
                           ( 37.0_wp * ibit26 * adv_mom_5                        &
                        +     7.0_wp * ibit25 * adv_mom_3                        &
@@ -5648,5 +5687,5 @@
                                       )
 
-                diss_t(k) = - ABS( w_comp ) * (                              &
+                diss_t(k) = - ABS( w_comp ) * rho_air_zw(k) * (              &
                           ( 10.0_wp * ibit26 * adv_mom_5                        &
                        +     3.0_wp * ibit25 * adv_mom_3                        &
@@ -5672,6 +5711,6 @@
                                          + IBITS(wall_flags_0(k,j,i-1),19,1)  &
                                          + IBITS(wall_flags_0(k,j,i-1),20,1)  &
-                                         )                                    &   
-                  ) * ddx                                                     &
+                                         )                                    &
+                  ) * rho_air(k) * ddx                                        &
                +  ( v_comp(k)                                                 &
                                        * ( ibit21 + ibit22 + ibit23 )         &
@@ -5680,9 +5719,9 @@
                                          + IBITS(wall_flags_0(k,j-1,i),22,1)  &
                                          + IBITS(wall_flags_0(k,j-1,i),23,1)  &
-                                         )                                    &   
-                  ) * ddy                                                     &
-               +  ( w_comp                                                    &
+                                         )                                    &
+                  ) * rho_air(k) * ddy                                        &
+               +  ( w_comp * rho_air_zw(k)                                    &
                                        * ( ibit24 + ibit25 + ibit26 )         &
-                - ( w(k-1,j-1,i) + w(k-1,j,i) )                               &
+                - ( w(k-1,j-1,i) + w(k-1,j,i) ) * rho_air_zw(k-1)             &
                                        * ( IBITS(wall_flags_0(k-1,j,i),24,1)  &
                                          + IBITS(wall_flags_0(k-1,j,i),25,1)  &
@@ -5700,7 +5739,7 @@
                      -   swap_flux_y_local_v(k) - swap_diss_y_local_v(k)      &
                        ) * ddy                                                &
-                     + ( flux_t(k) + diss_t(k)                                &
-                     -   flux_d    - diss_d                                   &
-                       ) * ddzw(k)                                            &
+                     + ( ( flux_t(k) + diss_t(k) )                            &
+                     -   ( flux_d    + diss_d    )                            &
+                       ) * drho_air(k) * ddzw(k)                              &
                                             )  + v(k,j,i) * div
 
@@ -5772,5 +5811,5 @@
 
                 w_comp    = w(k,j-1,i) + w(k,j,i)
-                flux_t(k) = w_comp  * (                                      &
+                flux_t(k) = w_comp * rho_air_zw(k) * (                       &
                           ( 37.0_wp * ibit26 * adv_mom_5                        &
                        +     7.0_wp * ibit25 * adv_mom_3                        &
@@ -5787,5 +5826,5 @@
                                       )
 
-                diss_t(k) = - ABS( w_comp ) * (                              &
+                diss_t(k) = - ABS( w_comp ) * rho_air_zw(k) * (              &
                           ( 10.0_wp * ibit26 * adv_mom_5                        &
                        +     3.0_wp * ibit25 * adv_mom_3                        &
@@ -5806,7 +5845,10 @@
 !--             by a not sufficient reduction of divergences near topography.
                 div = ( ( u_comp + gu - ( u(k,j-1,i)   + u(k,j,i)   ) ) * ddx &
+                                                                  * rho_air(k)&
                      +  ( v_comp(k)   - ( v(k,j,i)     + v(k,j-1,i) ) ) * ddy &
-                     +  ( w_comp      - ( w(k-1,j-1,i) + w(k-1,j,i) ) )       &
-                                                                    * ddzw(k) &
+                                                                  * rho_air(k)&
+                     +  (   w_comp                      * rho_air_zw(k) -     &
+                          ( w(k-1,j-1,i) + w(k-1,j,i) ) * rho_air_zw(k-1)     &
+                        ) * ddzw(k)                                           &
                       ) * 0.5_wp
  
@@ -5818,7 +5860,7 @@
                      -   swap_flux_y_local_v(k) - swap_diss_y_local_v(k)      &
                        ) * ddy                                                &
-                     + ( flux_t(k) + diss_t(k)                                &
-                     -   flux_d    - diss_d                                   &
-                       ) * ddzw(k)                                            &
+                     + ( ( flux_t(k) + diss_t(k) )                            &
+                     -   ( flux_d    + diss_d    )                            &
+                       ) * drho_air(k) * ddzw(k)                              &
                                             )  + v(k,j,i) * div
 
@@ -5869,5 +5911,5 @@
 
        USE arrays_3d,                                                          &
-           ONLY:  ddzw, tend, u, v, w
+           ONLY:  ddzw, drho_air, tend, u, v, w, rho_air, rho_air_zw
 
        USE constants,                                                          &
@@ -6132,5 +6174,5 @@
 
                 w_comp    = w(k,j-1,i) + w(k,j,i)
-                flux_t    = w_comp  * (                                      &
+                flux_t    = w_comp * rho_air_zw(k) * (                       &
                           ( 37.0_wp * ibit26 * adv_mom_5                        &
                        +     7.0_wp * ibit25 * adv_mom_3                        &
@@ -6147,5 +6189,5 @@
                                       )
 
-                diss_t    = - ABS( w_comp ) * (                              &
+                diss_t    = - ABS( w_comp ) * rho_air_zw(k) * (              &
                           ( 10.0_wp * ibit26 * adv_mom_5                        &
                        +     3.0_wp * ibit25 * adv_mom_3                        &
@@ -6171,6 +6213,6 @@
                                          + IBITS(wall_flags_0(k,j,i-1),19,1)  &
                                          + IBITS(wall_flags_0(k,j,i-1),20,1)  &
-                                         )                                    &   
-                  ) * ddx                                                     &
+                                         )                                    &
+                  ) * rho_air(k) * ddx                                        &
                +  ( v_comp                                                    &
                                        * ( ibit21 + ibit22 + ibit23 )         &
@@ -6179,9 +6221,9 @@
                                          + IBITS(wall_flags_0(k,j-1,i),22,1)  &
                                          + IBITS(wall_flags_0(k,j-1,i),23,1)  &
-                                         )                                    &   
-                  ) * ddy                                                     &
-               +  ( w_comp                                                    &
+                                         )                                    &
+                  ) * rho_air(k) * ddy                                        &
+               +  ( w_comp * rho_air_zw(k)                                    &
                                        * ( ibit24 + ibit25 + ibit26 )         &
-                - ( w(k-1,j-1,i) + w(k-1,j,i) )                               &
+                - ( w(k-1,j-1,i) + w(k-1,j,i) ) * rho_air_zw(k-1)             &
                                        * ( IBITS(wall_flags_0(k-1,j,i),24,1)  &
                                          + IBITS(wall_flags_0(k-1,j,i),25,1)  &
@@ -6195,5 +6237,7 @@
                                ( flux_r + diss_r - flux_l - diss_l ) * ddx     &
                              + ( flux_n + diss_n - flux_s - diss_s ) * ddy     &
-                             + ( flux_t + diss_t - flux_d - diss_d ) * ddzw(k) &
+                             + ( ( flux_t + diss_t ) -                         &
+                                 ( flux_d + diss_d )                           &
+                               ) * drho_air(k) * ddzw(k)                       &
                                 ) + div * v(k,j,i)
 
@@ -6235,5 +6279,5 @@
 
        USE arrays_3d,                                                          &
-           ONLY:  ddzu, tend, u, v, w
+           ONLY:  ddzu, drho_air_zw, tend, u, v, w, rho_air, rho_air_zw
 
        USE constants,                                                          &
@@ -6509,5 +6553,5 @@
 
                 w_comp    = w(k+1,j,i) + w(k,j,i)
-                flux_t(k) = w_comp  * (                                      &
+                flux_t(k) = w_comp * rho_air(k+1) * (                        &
                           ( 37.0_wp * ibit35 * adv_mom_5                        &
                        +     7.0_wp * ibit34 * adv_mom_3                        &
@@ -6524,5 +6568,5 @@
                                        )
 
-                diss_t(k) = - ABS( w_comp ) * (                              &
+                diss_t(k) = - ABS( w_comp ) * rho_air(k+1) * (               &
                           ( 10.0_wp * ibit35 * adv_mom_5                        &
                        +     3.0_wp * ibit34 * adv_mom_3                        &
@@ -6548,6 +6592,6 @@
                                       + IBITS(wall_flags_0(k,j,i-1),29,1)     &
                                       )                                       &
-                  ) * ddx                                                     &
-              +   ( ( v_comp + gv ) * ( ibit30 + ibit31 + ibit32 )            & 
+                  ) * rho_air_zw(k) * ddx                                     &
+              +   ( ( v_comp + gv ) * ( ibit30 + ibit31 + ibit32 )            &
                   - ( v(k+1,j,i) + v(k,j,i)   )                               &
                                     * ( IBITS(wall_flags_0(k,j-1,i),30,1)     &
@@ -6555,7 +6599,7 @@
                                       + IBITS(wall_flags_00(k,j-1,i),0,1)     &
                                       )                                       &
-                  ) * ddy                                                     &
-              +   ( w_comp          * ( ibit33 + ibit34 + ibit35 )            &
-                - ( w(k,j,i)   + w(k-1,j,i)   )                               &
+                  ) * rho_air_zw(k) * ddy                                     &
+              +   ( w_comp * rho_air(k+1) * ( ibit33 + ibit34 + ibit35 )      &
+                - ( w(k,j,i)   + w(k-1,j,i)   ) * rho_air(k)                  &
                                     * ( IBITS(wall_flags_00(k-1,j,i),1,1)     &
                                       + IBITS(wall_flags_00(k-1,j,i),2,1)     &
@@ -6574,7 +6618,7 @@
                     -   swap_flux_y_local_w(k)   - swap_diss_y_local_w(k)     &
                       ) * ddy                                                 &
-                    + ( flux_t(k) + diss_t(k)                                 &
-                    -   flux_d    - diss_d                                    &
-                      ) * ddzu(k+1)                                           &
+                    + ( ( flux_t(k) + diss_t(k) )                             &
+                    -   ( flux_d    + diss_d    )                             &
+                      ) * drho_air_zw(k) * ddzu(k+1)                          &
                                             )  + div * w(k,j,i)
 
@@ -6632,5 +6676,5 @@
 
                 w_comp    = w(k+1,j,i) + w(k,j,i)
-                flux_t(k) = w_comp  * (                                      &
+                flux_t(k) = w_comp * rho_air(k+1) * (                        &
                           ( 37.0_wp * ibit35 * adv_mom_5                        &
                        +     7.0_wp * ibit34 * adv_mom_3                        &
@@ -6647,5 +6691,5 @@
                                        )
 
-                diss_t(k) = - ABS( w_comp ) * (                              &
+                diss_t(k) = - ABS( w_comp ) * rho_air(k+1) * (               &
                           ( 10.0_wp * ibit35 * adv_mom_5                        &
                        +     3.0_wp * ibit34 * adv_mom_3                        &
@@ -6666,7 +6710,10 @@
 !--             by a not sufficient reduction of divergences near topography.
                 div = ( ( u_comp + gu - ( u(k+1,j,i) + u(k,j,i)   ) ) * ddx  &
+                                                             * rho_air_zw(k) &
                     +   ( v_comp + gv - ( v(k+1,j,i) + v(k,j,i)   ) ) * ddy  &
-                    +   ( w_comp      - ( w(k,j,i)   + w(k-1,j,i) ) )        &
-                                                                 * ddzu(k+1) &
+                                                             * rho_air_zw(k) &
+                    +   (   w_comp                    * rho_air(k+1) -       &
+                          ( w(k,j,i)   + w(k-1,j,i) ) * rho_air(k)           &
+                        ) * ddzu(k+1)                                        &
                       ) * 0.5_wp
 
@@ -6678,7 +6725,7 @@
                     -   swap_flux_y_local_w(k)   - swap_diss_y_local_w(k)     &
                       ) * ddy                                                 &
-                    + ( flux_t(k) + diss_t(k)                                 &
-                    -   flux_d    - diss_d                                    &
-                      ) * ddzu(k+1)                                           &
+                    + ( ( flux_t(k) + diss_t(k) )                             &
+                    -   ( flux_d    + diss_d    )                             &
+                      ) * drho_air_zw(k) * ddzu(k+1)                          &
                                             )  + div * w(k,j,i)
 
@@ -6714,5 +6761,5 @@
 
        USE arrays_3d,                                                          &
-           ONLY:  ddzu, tend, u, v, w
+           ONLY:  ddzu, drho_air_zw, tend, u, v, w, rho_air, rho_air_zw
 
        USE constants,                                                          &
@@ -6976,5 +7023,5 @@
 
                 w_comp    = w(k+1,j,i) + w(k,j,i)
-                flux_t    = w_comp  * (                                      &
+                flux_t    = w_comp * rho_air(k+1) * (                        &
                           ( 37.0_wp * ibit35 * adv_mom_5                        &
                        +     7.0_wp * ibit34 * adv_mom_3                        &
@@ -6991,5 +7038,5 @@
                                        )
 
-                diss_t    = - ABS( w_comp ) * (                              &
+                diss_t    = - ABS( w_comp ) * rho_air(k+1) * (               &
                           ( 10.0_wp * ibit35 * adv_mom_5                        &
                        +     3.0_wp * ibit34 * adv_mom_3                        &
@@ -7015,6 +7062,6 @@
                                       + IBITS(wall_flags_0(k,j,i-1),29,1)     &
                                       )                                       &
-                  ) * ddx                                                     &
-              +   ( ( v_comp + gv ) * ( ibit30 + ibit31 + ibit32 )            & 
+                  ) * rho_air_zw(k) * ddx                                     &
+              +   ( ( v_comp + gv ) * ( ibit30 + ibit31 + ibit32 )            &
                   - ( v(k+1,j,i) + v(k,j,i)   )                               &
                                     * ( IBITS(wall_flags_0(k,j-1,i),30,1)     &
@@ -7022,7 +7069,7 @@
                                       + IBITS(wall_flags_00(k,j-1,i),0,1)     &
                                       )                                       &
-                  ) * ddy                                                     &
-              +   ( w_comp          * ( ibit33 + ibit34 + ibit35 )            &
-                - ( w(k,j,i)   + w(k-1,j,i)   )                               &
+                  ) * rho_air_zw(k) * ddy                                     &
+              +   ( w_comp * rho_air(k+1) * ( ibit33 + ibit34 + ibit35 )      &
+                - ( w(k,j,i)   + w(k-1,j,i)   ) * rho_air(k)                  &
                                     * ( IBITS(wall_flags_00(k-1,j,i),1,1)     &
                                       + IBITS(wall_flags_00(k-1,j,i),2,1)     &
@@ -7036,5 +7083,7 @@
                                ( flux_r + diss_r - flux_l - diss_l ) * ddx       &
                              + ( flux_n + diss_n - flux_s - diss_s ) * ddy       &
-                             + ( flux_t + diss_t - flux_d - diss_d ) * ddzu(k+1) &
+                             + ( ( flux_t + diss_t ) -                           &
+                                 ( flux_d + diss_d ) * rho_air(k)                &
+                               ) * drho_air_zw(k) * ddzu(k+1)                    &
                                  ) + div * w(k,j,i)
 
Index: /palm/trunk/SOURCE/check_parameters.f90
===================================================================
--- /palm/trunk/SOURCE/check_parameters.f90	(revision 2036)
+++ /palm/trunk/SOURCE/check_parameters.f90	(revision 2037)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! Anelastic approximation implemented
 ! 
 ! Former revisions:
@@ -432,5 +432,6 @@
     USE netcdf_interface,                                                      &
         ONLY:  dopr_unit, do2d_unit, do3d_unit, netcdf_data_format,            &
-               netcdf_data_format_string
+               netcdf_data_format_string, dots_unit, heatflux_output_unit,     &
+               waterflux_output_unit, momentumflux_output_unit
 
     USE particle_attributes
@@ -837,4 +838,89 @@
 
     ENDIF
+
+!
+!-- Check approximation
+    IF ( TRIM( approximation ) /= 'boussinesq'   .AND.   &
+         TRIM( approximation ) /= 'anelastic' )  THEN
+       message_string = 'unknown approximation: approximation = "' //    &
+                        TRIM( approximation ) // '"'
+       CALL message( 'check_parameters', 'PA0000', 1, 2, 0, 6, 0 )
+    ENDIF
+
+!
+!-- Check approximation requirements
+    IF ( TRIM( approximation ) == 'anelastic'   .AND.   &
+         TRIM( momentum_advec ) /= 'ws-scheme' )  THEN
+       message_string = 'Anelastic approximation requires: ' //                &
+                        'momentum_advec = "ws-scheme"'
+       CALL message( 'check_parameters', 'PA0000', 1, 2, 0, 6, 0 )
+    ENDIF
+    IF ( TRIM( approximation ) == 'anelastic'   .AND.   &
+         TRIM( psolver ) == 'multigrid' )  THEN
+       message_string = 'Anelastic approximation currently only supports: ' // &
+                        'psolver = "poisfft", ' // &
+                        'psolver = "sor" and ' // &
+                        'psolver = "multigrid_noopt"'
+       CALL message( 'check_parameters', 'PA0000', 1, 2, 0, 6, 0 )
+    ENDIF
+    IF ( TRIM( approximation ) == 'anelastic'   .AND.   &
+         conserve_volume_flow )  THEN
+       message_string = 'Anelastic approximation is not allowed with:' //      &
+                        'conserve_volume_flow = .TRUE.'
+       CALL message( 'check_parameters', 'PA0000', 1, 2, 0, 6, 0 )
+    ENDIF
+
+!
+!-- Check flux input mode
+    IF ( TRIM( flux_input_mode ) /= 'dynamic'    .AND.   &
+         TRIM( flux_input_mode ) /= 'kinematic'  .AND.   &
+         TRIM( flux_input_mode ) /= 'approximation-specific' )  THEN
+       message_string = 'unknown flux input mode: flux_input_mode = "' //      &
+                        TRIM( flux_input_mode ) // '"'
+       CALL message( 'check_parameters', 'PA0000', 1, 2, 0, 6, 0 )
+    ENDIF
+!-- Set flux input mode according to approximation if applicable
+    IF ( TRIM( flux_input_mode ) == 'approximation-specific' )  THEN
+       IF ( TRIM( approximation ) == 'anelastic' )  THEN
+          flux_input_mode = 'dynamic'
+       ELSEIF ( TRIM( approximation ) == 'boussinesq' )  THEN
+          flux_input_mode = 'kinematic'
+       ENDIF
+    ENDIF
+
+!
+!-- Check flux output mode
+    IF ( TRIM( flux_output_mode ) /= 'dynamic'    .AND.   &
+         TRIM( flux_output_mode ) /= 'kinematic'  .AND.   &
+         TRIM( flux_output_mode ) /= 'approximation-specific' )  THEN
+       message_string = 'unknown flux output mode: flux_output_mode = "' //    &
+                        TRIM( flux_output_mode ) // '"'
+       CALL message( 'check_parameters', 'PA0000', 1, 2, 0, 6, 0 )
+    ENDIF
+!-- Set flux output mode according to approximation if applicable
+    IF ( TRIM( flux_output_mode ) == 'approximation-specific' )  THEN
+       IF ( TRIM( approximation ) == 'anelastic' )  THEN
+          flux_output_mode = 'dynamic'
+       ELSEIF ( TRIM( approximation ) == 'boussinesq' )  THEN
+          flux_output_mode = 'kinematic'
+       ENDIF
+    ENDIF
+
+!
+!-- set the flux output units according to flux_output_mode
+    IF ( TRIM( flux_output_mode ) == 'kinematic' ) THEN
+        heatflux_output_unit              = 'K m/s'
+        waterflux_output_unit             = 'kg/kg m/s'
+        momentumflux_output_unit          = 'm2/s2'
+    ELSEIF ( TRIM( flux_output_mode ) == 'dynamic' ) THEN
+        heatflux_output_unit              = 'W/m2'
+        waterflux_output_unit             = 'W/m2'
+        momentumflux_output_unit          = 'N/m2'
+    ENDIF
+
+!-- set time series output units for fluxes
+    dots_unit(14:16) = heatflux_output_unit
+    dots_unit(21)    = waterflux_output_unit
+    dots_unit(19:20) = momentumflux_output_unit
 
 !
@@ -2115,5 +2201,5 @@
           CASE ( 'w"u"' )
              dopr_index(i) = 12
-             dopr_unit(i)  = 'm2/s2'
+             dopr_unit(i)  = TRIM ( momentumflux_output_unit )
              hom(:,2,12,:) = SPREAD( zw, 2, statistic_regions+1 )
              IF ( constant_flux_layer )  hom(nzb,2,12,:) = zu(1)
@@ -2121,10 +2207,10 @@
           CASE ( 'w*u*' )
              dopr_index(i) = 13
-             dopr_unit(i)  = 'm2/s2'
+             dopr_unit(i)  = TRIM ( momentumflux_output_unit )
              hom(:,2,13,:) = SPREAD( zw, 2, statistic_regions+1 )
 
           CASE ( 'w"v"' )
              dopr_index(i) = 14
-             dopr_unit(i)  = 'm2/s2'
+             dopr_unit(i)  = TRIM ( momentumflux_output_unit )
              hom(:,2,14,:) = SPREAD( zw, 2, statistic_regions+1 )
              IF ( constant_flux_layer )  hom(nzb,2,14,:) = zu(1)
@@ -2132,25 +2218,25 @@
           CASE ( 'w*v*' )
              dopr_index(i) = 15
-             dopr_unit(i)  = 'm2/s2'
+             dopr_unit(i)  = TRIM ( momentumflux_output_unit )
              hom(:,2,15,:) = SPREAD( zw, 2, statistic_regions+1 )
 
           CASE ( 'w"pt"' )
              dopr_index(i) = 16
-             dopr_unit(i)  = 'K m/s'
+             dopr_unit(i)  = TRIM ( heatflux_output_unit )
              hom(:,2,16,:) = SPREAD( zw, 2, statistic_regions+1 )
 
           CASE ( 'w*pt*' )
              dopr_index(i) = 17
-             dopr_unit(i)  = 'K m/s'
+             dopr_unit(i)  = TRIM ( heatflux_output_unit )
              hom(:,2,17,:) = SPREAD( zw, 2, statistic_regions+1 )
 
           CASE ( 'wpt' )
              dopr_index(i) = 18
-             dopr_unit(i)  = 'K m/s'
+             dopr_unit(i)  = TRIM ( heatflux_output_unit )
              hom(:,2,18,:) = SPREAD( zw, 2, statistic_regions+1 )
 
           CASE ( 'wu' )
              dopr_index(i) = 19
-             dopr_unit(i)  = 'm2/s2'
+             dopr_unit(i)  = TRIM ( momentumflux_output_unit )
              hom(:,2,19,:) = SPREAD( zw, 2, statistic_regions+1 )
              IF ( constant_flux_layer )  hom(nzb,2,19,:) = zu(1)
@@ -2158,5 +2244,5 @@
           CASE ( 'wv' )
              dopr_index(i) = 20
-             dopr_unit(i)  = 'm2/s2'
+             dopr_unit(i)  = TRIM ( momentumflux_output_unit )
              hom(:,2,20,:) = SPREAD( zw, 2, statistic_regions+1 )
              IF ( constant_flux_layer )  hom(nzb,2,20,:) = zu(1)
@@ -2164,10 +2250,10 @@
           CASE ( 'w*pt*BC' )
              dopr_index(i) = 21
-             dopr_unit(i)  = 'K m/s'
+             dopr_unit(i)  = TRIM ( heatflux_output_unit )
              hom(:,2,21,:) = SPREAD( zw, 2, statistic_regions+1 )
 
           CASE ( 'wptBC' )
              dopr_index(i) = 22
-             dopr_unit(i)  = 'K m/s'
+             dopr_unit(i)  = TRIM ( heatflux_output_unit )
              hom(:,2,22,:) = SPREAD( zw, 2, statistic_regions+1 )
 
@@ -2335,15 +2421,15 @@
           CASE ( 'w"vpt"' )
              dopr_index(i) = 45
-             dopr_unit(i)  = 'K m/s'
+             dopr_unit(i)  = TRIM ( heatflux_output_unit )
              hom(:,2,45,:) = SPREAD( zw, 2, statistic_regions+1 )
 
           CASE ( 'w*vpt*' )
              dopr_index(i) = 46
-             dopr_unit(i)  = 'K m/s'
+             dopr_unit(i)  = TRIM ( heatflux_output_unit )
              hom(:,2,46,:) = SPREAD( zw, 2, statistic_regions+1 )
 
           CASE ( 'wvpt' )
              dopr_index(i) = 47
-             dopr_unit(i)  = 'K m/s'
+             dopr_unit(i)  = TRIM ( heatflux_output_unit )
              hom(:,2,47,:) = SPREAD( zw, 2, statistic_regions+1 )
 
@@ -2356,5 +2442,5 @@
              ELSE
                 dopr_index(i) = 48
-                dopr_unit(i)  = 'kg/kg m/s'
+                dopr_unit(i)  = TRIM ( waterflux_output_unit )
                 hom(:,2,48,:) = SPREAD( zw, 2, statistic_regions+1 )
              ENDIF
@@ -2368,5 +2454,5 @@
              ELSE
                 dopr_index(i) = 49
-                dopr_unit(i)  = 'kg/kg m/s'
+                dopr_unit(i)  = TRIM ( waterflux_output_unit )
                 hom(:,2,49,:) = SPREAD( zw, 2, statistic_regions+1 )
              ENDIF
@@ -2380,5 +2466,5 @@
              ELSE
                 dopr_index(i) = 50
-                dopr_unit(i)  = 'kg/kg m/s'
+                dopr_unit(i)  = TRIM ( waterflux_output_unit )
                 hom(:,2,50,:) = SPREAD( zw, 2, statistic_regions+1 )
              ENDIF
@@ -2423,9 +2509,9 @@
              IF ( humidity  .AND.  .NOT.  cloud_physics )  THEN
                 dopr_index(i) = 48
-                dopr_unit(i)  = 'kg/kg m/s'
+                dopr_unit(i)  = TRIM ( waterflux_output_unit )
                 hom(:,2,48,:) = SPREAD( zw, 2, statistic_regions+1 )
              ELSEIF ( humidity  .AND.  cloud_physics )  THEN
                 dopr_index(i) = 51
-                dopr_unit(i)  = 'kg/kg m/s'
+                dopr_unit(i)  = TRIM ( waterflux_output_unit )
                 hom(:,2,51,:) = SPREAD( zw, 2, statistic_regions+1 )
              ELSE
@@ -2441,9 +2527,9 @@
              THEN
                 dopr_index(i) = 49
-                dopr_unit(i)  = 'kg/kg m/s'
+                dopr_unit(i)  = TRIM ( waterflux_output_unit )
                 hom(:,2,49,:) = SPREAD( zw, 2, statistic_regions+1 )
              ELSEIF( humidity .AND. cloud_physics ) THEN
                 dopr_index(i) = 52
-                dopr_unit(i)  = 'kg/kg m/s'
+                dopr_unit(i)  = TRIM ( waterflux_output_unit )
                 hom(:,2,52,:) = SPREAD( zw, 2, statistic_regions+1 )
              ELSE
@@ -2458,9 +2544,9 @@
              IF ( humidity  .AND.  .NOT.  cloud_physics )  THEN
                 dopr_index(i) = 50
-                dopr_unit(i)  = 'kg/kg m/s'
+                dopr_unit(i)  = TRIM ( waterflux_output_unit )
                 hom(:,2,50,:) = SPREAD( zw, 2, statistic_regions+1 )
              ELSEIF ( humidity  .AND.  cloud_physics )  THEN
                 dopr_index(i) = 53
-                dopr_unit(i)  = 'kg/kg m/s'
+                dopr_unit(i)  = TRIM ( waterflux_output_unit )
                 hom(:,2,53,:) = SPREAD( zw, 2, statistic_regions+1 )
              ELSE
@@ -2503,30 +2589,30 @@
           CASE ( 'u"pt"' )
              dopr_index(i) = 58
-             dopr_unit(i)  = 'K m/s'
+             dopr_unit(i)  = TRIM ( heatflux_output_unit )
              hom(:,2,58,:) = SPREAD( zu, 2, statistic_regions+1 )
 
           CASE ( 'u*pt*' )
              dopr_index(i) = 59
-             dopr_unit(i)  = 'K m/s'
+             dopr_unit(i)  = TRIM ( heatflux_output_unit )
              hom(:,2,59,:) = SPREAD( zu, 2, statistic_regions+1 )
 
           CASE ( 'upt_t' )
              dopr_index(i) = 60
-             dopr_unit(i)  = 'K m/s'
+             dopr_unit(i)  = TRIM ( heatflux_output_unit )
              hom(:,2,60,:) = SPREAD( zu, 2, statistic_regions+1 )
 
           CASE ( 'v"pt"' )
              dopr_index(i) = 61
-             dopr_unit(i)  = 'K m/s'
+             dopr_unit(i)  = TRIM ( heatflux_output_unit )
              hom(:,2,61,:) = SPREAD( zu, 2, statistic_regions+1 )
              
           CASE ( 'v*pt*' )
              dopr_index(i) = 62
-             dopr_unit(i)  = 'K m/s'
+             dopr_unit(i)  = TRIM ( heatflux_output_unit )
              hom(:,2,62,:) = SPREAD( zu, 2, statistic_regions+1 )
 
           CASE ( 'vpt_t' )
              dopr_index(i) = 63
-             dopr_unit(i)  = 'K m/s'
+             dopr_unit(i)  = TRIM ( heatflux_output_unit )
              hom(:,2,63,:) = SPREAD( zu, 2, statistic_regions+1 )
 
@@ -2621,6 +2707,16 @@
           CASE ( 'hyp' )
              dopr_index(i) = 72
-             dopr_unit(i)  = 'dbar'
+             dopr_unit(i)  = 'hPa'
              hom(:,2,72,:) = SPREAD( zu, 2, statistic_regions+1 )
+
+          CASE ( 'rho_air' )
+             dopr_index(i)  = 121
+             dopr_unit(i)   = 'kg/m3'
+             hom(:,2,121,:) = SPREAD( zu, 2, statistic_regions+1 )
+
+          CASE ( 'rho_air_zw' )
+             dopr_index(i)  = 122
+             dopr_unit(i)   = 'kg/m3'
+             hom(:,2,122,:) = SPREAD( zw, 2, statistic_regions+1 )
 
           CASE ( 'nr' )
Index: /palm/trunk/SOURCE/diffusion_e.f90
===================================================================
--- /palm/trunk/SOURCE/diffusion_e.f90	(revision 2036)
+++ /palm/trunk/SOURCE/diffusion_e.f90	(revision 2037)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! Anelastic approximation implemented
 ! 
 ! Former revisions:
@@ -128,5 +128,6 @@
 
        USE arrays_3d,                                                          &
-           ONLY:  dd2zu, ddzu, ddzw, diss, e, km, l_grid, tend, zu, zw
+           ONLY:  dd2zu, ddzu, ddzw, diss, e, km, l_grid, tend, zu, zw,        &
+                  drho_air, rho_air_zw
            
        USE control_parameters,                                                 &
@@ -222,6 +223,8 @@
                                         + (                                    &
                ( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1) &
+                                                             * rho_air_zw(k)   &
              - ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k)   &
-                                          ) * ddzw(k)                          &
+                                                             * rho_air_zw(k-1) &
+                                          ) * ddzw(k) * drho_air(k)            &
                              - dissipation(k,j)
 
@@ -294,6 +297,8 @@
                                         + (                                    &
                ( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1) &
+                                                             * rho_air_zw(k)   &
              - ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k)   &
-                                          ) * ddzw(k)                          &
+                                                             * rho_air_zw(k-1) &
+                                          ) * ddzw(k) * drho_air(k)            &
                              - dissipation(k,j)
 
@@ -339,5 +344,6 @@
 
        USE arrays_3d,                                                          &
-           ONLY:  dd2zu, ddzu, ddzw, diss, e, km, l_grid, tend, zu, zw
+           ONLY:  dd2zu, ddzu, ddzw, diss, e, km, l_grid, tend, zu, zw,        &
+                  drho_air, rho_air_zw
           
        USE control_parameters,                                                 &
@@ -430,6 +436,8 @@
                                         + (                                    &
                ( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1) &
+                                                             * rho_air_zw(k)   &
              - ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k)   &
-                                          ) * ddzw(k)                          &
+                                                             * rho_air_zw(k-1) &
+                                          ) * ddzw(k) * drho_air(k)            &
                                   - dissipation
 
@@ -497,6 +505,8 @@
                                         + (                                    &
                ( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1) &
+                                                             * rho_air_zw(k)   &
              - ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k)   &
-                                          ) * ddzw(k)                          &
+                                                             * rho_air_zw(k-1) &
+                                          ) * ddzw(k) * drho_air(k)            &
                                   - dissipation
 
@@ -542,5 +552,6 @@
 
        USE arrays_3d,                                                          &
-           ONLY:  dd2zu, ddzu, ddzw, diss, e, km, l_grid, tend, zu, zw
+           ONLY:  dd2zu, ddzu, ddzw, diss, e, km, l_grid, tend, zu, zw,        &
+                  drho_air, rho_air_zw
           
        USE control_parameters,                                                 &
@@ -625,6 +636,8 @@
                                        + (                                    &
               ( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1) &
+                                                            * rho_air_zw(k)   &
             - ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k)   &
-                                         ) * ddzw(k)                          &
+                                                            * rho_air_zw(k-1) &
+                                         ) * ddzw(k) * drho_air(k)            &
                                        - dissipation(k)
 
Index: /palm/trunk/SOURCE/diffusion_s.f90
===================================================================
--- /palm/trunk/SOURCE/diffusion_s.f90	(revision 2036)
+++ /palm/trunk/SOURCE/diffusion_s.f90	(revision 2037)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! ------------------
-! 
+! Anelastic approximation implemented
 ! 
 ! Former revisions:
@@ -113,5 +113,5 @@
 
        USE arrays_3d,                                                          &
-           ONLY:  ddzu, ddzw, kh, tend
+           ONLY:  ddzu, ddzw, kh, tend, drho_air, rho_air_zw
        
        USE control_parameters,                                                 & 
@@ -190,6 +190,8 @@
                                        + 0.5_wp * (                            &
             ( kh(k,j,i) + kh(k+1,j,i) ) * ( s(k+1,j,i)-s(k,j,i) ) * ddzu(k+1)  &
+                                                            * rho_air_zw(k)    &
           - ( kh(k,j,i) + kh(k-1,j,i) ) * ( s(k,j,i)-s(k-1,j,i) ) * ddzu(k)    &
-                                                  ) * ddzw(k)
+                                                            * rho_air_zw(k-1)  &
+                                                  ) * ddzw(k) * drho_air(k)
              ENDDO
 
@@ -205,6 +207,7 @@
                                                   * ( s(k+1,j,i)-s(k,j,i) )    &
                                                   * ddzu(k+1)                  &
+                                                  * rho_air_zw(k)              &
                                            + s_flux_b(j,i)                     &
-                                         ) * ddzw(k)
+                                         ) * ddzw(k) * drho_air(k)
 
              ENDIF
@@ -222,5 +225,6 @@
                                                     * ( s(k,j,i)-s(k-1,j,i) )  &
                                                     * ddzu(k)                  &
-                                         ) * ddzw(k)
+                                                    * rho_air_zw(k-1)          &
+                                         ) * ddzw(k) * drho_air(k)
 
              ENDIF
@@ -240,5 +244,5 @@
 
        USE arrays_3d,                                                          &
-           ONLY:  ddzu, ddzw, kh, tend
+           ONLY:  ddzu, ddzw, kh, tend, drho_air, rho_air_zw
            
        USE control_parameters,                                                 & 
@@ -324,6 +328,8 @@
                                        + 0.5_wp * (                            &
             ( kh(k,j,i) + kh(k+1,j,i) ) * ( s(k+1,j,i)-s(k,j,i) ) * ddzu(k+1)  &
+                                                            * rho_air_zw(k)    &
           - ( kh(k,j,i) + kh(k-1,j,i) ) * ( s(k,j,i)-s(k-1,j,i) ) * ddzu(k)    &
-                                                  ) * ddzw(k)
+                                                            * rho_air_zw(k-1)  &
+                                                  ) * ddzw(k) * drho_air(k)
                 ENDIF
              ENDDO
@@ -338,6 +344,7 @@
                                                      * ( s(k+1,j,i)-s(k,j,i) ) &
                                                      * ddzu(k+1)               &
+                                                     * rho_air_zw(k)           &
                                               + s_flux_b(j,i)                  &
-                                            ) * ddzw(k)
+                                            ) * ddzw(k) * drho_air(k)
                 ENDIF
 
@@ -351,5 +358,6 @@
                                                        * ( s(k,j,i)-s(k-1,j,i) )  &
                                                        * ddzu(k)                  &
-                                            ) * ddzw(k)
+                                                       * rho_air_zw(k-1)          &
+                                            ) * ddzw(k) * drho_air(k)
                 ENDIF
              ENDDO
@@ -370,5 +378,5 @@
 
        USE arrays_3d,                                                          &
-           ONLY:  ddzu, ddzw, kh, tend
+           ONLY:  ddzu, ddzw, kh, tend, drho_air, rho_air_zw
            
        USE control_parameters,                                                 & 
@@ -446,6 +454,8 @@
                                        + 0.5_wp * (                            &
             ( kh(k,j,i) + kh(k+1,j,i) ) * ( s(k+1,j,i)-s(k,j,i) ) * ddzu(k+1)  &
+                                                            * rho_air_zw(k)    &
           - ( kh(k,j,i) + kh(k-1,j,i) ) * ( s(k,j,i)-s(k-1,j,i) ) * ddzu(k)    &
-                                                  ) * ddzw(k)
+                                                            * rho_air_zw(k-1)  &
+                                                  ) * ddzw(k) * drho_air(k)
        ENDDO
 
@@ -459,6 +469,7 @@
                                                * ( s(k+1,j,i)-s(k,j,i) )       &
                                                * ddzu(k+1)                     &
+                                               * rho_air_zw(k)                 &
                                         + s_flux_b(j,i)                        &
-                                      ) * ddzw(k)
+                                      ) * ddzw(k) * drho_air(k)
 
        ENDIF
@@ -474,5 +485,6 @@
                                                * ( s(k,j,i)-s(k-1,j,i) )       &
                                                * ddzu(k)                       &
-                                      ) * ddzw(k)
+                                               * rho_air_zw(k-1)               &
+                                      ) * ddzw(k) * drho_air(k)
 
        ENDIF
Index: /palm/trunk/SOURCE/diffusion_u.f90
===================================================================
--- /palm/trunk/SOURCE/diffusion_u.f90	(revision 2036)
+++ /palm/trunk/SOURCE/diffusion_u.f90	(revision 2037)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! Anelastic approximation implemented
 ! 
 ! Former revisions:
@@ -116,5 +116,6 @@
 
        USE arrays_3d,                                                          &
-           ONLY:  ddzu, ddzw, km, tend, u, usws, uswst, v, w
+           ONLY:  ddzu, ddzw, km, tend, u, usws, uswst, v, w,                  &
+                  drho_air, rho_air_zw
        
        USE control_parameters,                                                 &
@@ -217,9 +218,9 @@
                       & + ( kmzp * ( ( u(k+1,j,i) - u(k,j,i)   ) * ddzu(k+1)   &
                       &            + ( w(k,j,i)   - w(k,j,i-1) ) * ddx         &
-                      &            )                                           &
+                      &            ) * rho_air_zw(k)                           &
                       &   - kmzm * ( ( u(k,j,i)   - u(k-1,j,i)   ) * ddzu(k)   &
                       &            + ( w(k-1,j,i) - w(k-1,j,i-1) ) * ddx       &
-                      &            )                                           &
-                      &   ) * ddzw(k)
+                      &            ) * rho_air_zw(k-1)                         &
+                      &   ) * ddzw(k) * drho_air(k)
              ENDDO
 
@@ -242,9 +243,9 @@
 
                 tend(k,j,i) = tend(k,j,i)                                      &
-                      & + ( kmzp * ( w(k,j,i)   - w(k,j,i-1)   ) * ddx         &
-                      &   ) * ddzw(k)                                          &
-                      & + ( kmzp * ( u(k+1,j,i) - u(k,j,i)     ) * ddzu(k+1)   &
-                      &   + usws(j,i)                                          &
-                      &   ) * ddzw(k)
+                      & + ( kmzp * ( ( u(k+1,j,i) - u(k,j,i)   ) * ddzu(k+1)   &
+                      &            + ( w(k,j,i)   - w(k,j,i-1) ) * ddx         &
+                      &            ) * rho_air_zw(k)                           &
+                      &   - ( -usws(j,i) )                                     &
+                      &   ) * ddzw(k) * drho_air(k)
              ENDIF
 
@@ -260,9 +261,9 @@
 
                 tend(k,j,i) = tend(k,j,i)                                      &
-                      & - ( kmzm * ( w(k-1,j,i) - w(k-1,j,i-1) ) * ddx         &
-                      &   ) * ddzw(k)                                          &
-                      & + ( -uswst(j,i)                                        &
-                      &   - kmzm * ( u(k,j,i)   - u(k-1,j,i)   ) * ddzu(k)     &
-                      &   ) * ddzw(k)
+                      & + ( ( -uswst(j,i) )                                    &
+                      &   - kmzm * ( ( u(k,j,i)   - u(k-1,j,i)   ) * ddzu(k)   &
+                      &            + ( w(k-1,j,i) - w(k-1,j,i-1) ) * ddx       &
+                      &            ) * rho_air_zw(k-1)                         &
+                      &   ) * ddzw(k) * drho_air(k)
              ENDIF
 
@@ -281,5 +282,6 @@
 
        USE arrays_3d,                                                          &
-           ONLY:  ddzu, ddzw, km, tend, u, usws, uswst, v, w
+           ONLY:  ddzu, ddzw, km, tend, u, usws, uswst, v, w,                  &
+                  drho_air, rho_air_zw
        
        USE control_parameters,                                                 &
@@ -390,9 +392,9 @@
                          & + ( kmzp * ( ( u(k+1,j,i) - u(k,j,i)   ) * ddzu(k+1)&
                          &            + ( w(k,j,i)   - w(k,j,i-1) ) * ddx      &
-                         &            )                                        &
+                         &            ) * rho_air_zw(k)                        &
                          &   - kmzm * ( ( u(k,j,i)   - u(k-1,j,i)   ) * ddzu(k)&
                          &            + ( w(k-1,j,i) - w(k-1,j,i-1) ) * ddx    &
-                         &            )                                        &
-                         &   ) * ddzw(k)
+                         &            ) * rho_air_zw(k-1)                      &
+                         &   ) * ddzw(k) * drho_air(k)
                 ENDIF
              ENDDO
@@ -423,9 +425,9 @@
 
                 tend(k,j,i) = tend(k,j,i)                                      &
-                      & + ( kmzp * ( w(k,j,i)   - w(k,j,i-1)   ) * ddx         &
-                      &   ) * ddzw(k)                                          &
-                      & + ( kmzp * ( u(k+1,j,i) - u(k,j,i)     ) * ddzu(k+1)   &
-                      &   + usws(j,i)                                          &
-                      &   ) * ddzw(k)
+                      & + ( kmzp * ( ( u(k+1,j,i) - u(k,j,i)   ) * ddzu(k+1)   &
+                      &            + ( w(k,j,i)   - w(k,j,i-1) ) * ddx         &
+                      &            ) * rho_air_zw(k)                           &
+                      &   - ( -usws(j,i) )                                     &
+                      &   ) * ddzw(k) * drho_air(k)
              ENDDO
           ENDDO
@@ -449,9 +451,9 @@
 
                 tend(k,j,i) = tend(k,j,i)                                      &
-                      & - ( kmzm * ( w(k-1,j,i) - w(k-1,j,i-1) ) * ddx         &
-                      &   ) * ddzw(k)                                          &
-                      & + ( -uswst(j,i)                                        &
-                      &   - kmzm * ( u(k,j,i)   - u(k-1,j,i)   ) * ddzu(k)     &
-                      &   ) * ddzw(k)
+                      & + ( ( -uswst(j,i) )                                    &
+                      &   - kmzm * ( ( u(k,j,i)   - u(k-1,j,i)   ) * ddzu(k)   &
+                      &            + ( w(k-1,j,i) - w(k-1,j,i-1) ) * ddx       &
+                      &            ) * rho_air_zw(k-1)                         &
+                      &   ) * ddzw(k) * drho_air(k)
              ENDDO
           ENDDO
@@ -471,5 +473,6 @@
 
        USE arrays_3d,                                                          &
-           ONLY:  ddzu, ddzw, km, tend, u, usws, uswst, v, w
+           ONLY:  ddzu, ddzw, km, tend, u, usws, uswst, v, w,                  &
+                  drho_air, rho_air_zw
        
        USE control_parameters,                                                 &
@@ -559,9 +562,9 @@
                       & + ( kmzp * ( ( u(k+1,j,i) - u(k,j,i)   ) * ddzu(k+1)   &
                       &            + ( w(k,j,i)   - w(k,j,i-1) ) * ddx         &
-                      &            )                                           &
+                      &            ) * rho_air_zw(k)                           &
                       &   - kmzm * ( ( u(k,j,i)   - u(k-1,j,i)   ) * ddzu(k)   &
                       &            + ( w(k-1,j,i) - w(k-1,j,i-1) ) * ddx       &
-                      &            )                                           &
-                      &   ) * ddzw(k)
+                      &            ) * rho_air_zw(k-1)                         &
+                      &   ) * ddzw(k) * drho_air(k)
        ENDDO
 
@@ -583,9 +586,9 @@
 
           tend(k,j,i) = tend(k,j,i)                                            &
-                      & + ( kmzp * ( w(k,j,i)   - w(k,j,i-1)   ) * ddx         &
-                      &   ) * ddzw(k)                                          &
-                      & + ( kmzp * ( u(k+1,j,i) - u(k,j,i)     ) * ddzu(k+1)   &
-                      &   + usws(j,i)                                          &
-                      &   ) * ddzw(k)
+                      & + ( kmzp * ( ( u(k+1,j,i) - u(k,j,i)   ) * ddzu(k+1)   &
+                      &            + ( w(k,j,i)   - w(k,j,i-1) ) * ddx         &
+                      &            ) * rho_air_zw(k)                           &
+                      &   - ( -usws(j,i) )                                     &
+                      &   ) * ddzw(k) * drho_air(k)
        ENDIF
 
@@ -597,13 +600,12 @@
 !
 !--       Interpolate eddy diffusivities on staggered gridpoints
-          kmzm = 0.25_wp *                                                     &
-                 ( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
+          kmzm = 0.25_wp * ( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
 
           tend(k,j,i) = tend(k,j,i)                                            &
-                      & - ( kmzm * ( w(k-1,j,i) - w(k-1,j,i-1) ) * ddx         &
-                      &   ) * ddzw(k)                                          &
-                      & + ( -uswst(j,i)                                        &
-                      &   - kmzm * ( u(k,j,i)   - u(k-1,j,i)   ) * ddzu(k)     &
-                      &   ) * ddzw(k)
+                      & + ( ( -uswst(j,i) )                                    &
+                      &   - kmzm * ( ( u(k,j,i)   - u(k-1,j,i)   ) * ddzu(k)   &
+                      &            + ( w(k-1,j,i) - w(k-1,j,i-1) ) * ddx       &
+                      &            ) * rho_air_zw(k-1)                         &
+                      &   ) * ddzw(k) * drho_air(k)
        ENDIF
 
Index: /palm/trunk/SOURCE/diffusion_v.f90
===================================================================
--- /palm/trunk/SOURCE/diffusion_v.f90	(revision 2036)
+++ /palm/trunk/SOURCE/diffusion_v.f90	(revision 2037)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! Anelastic approximation implemented
 ! 
 ! Former revisions:
@@ -111,5 +111,6 @@
 
        USE arrays_3d,                                                          &
-           ONLY:  ddzu, ddzw, km, tend, u, v, vsws, vswst, w
+           ONLY:  ddzu, ddzw, km, tend, u, v, vsws, vswst, w,                  &
+                  drho_air, rho_air_zw
        
        USE control_parameters,                                                 &
@@ -212,9 +213,9 @@
                       & + ( kmzp * ( ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1)     &
                       &            + ( w(k,j,i) - w(k,j-1,i) ) * ddy           &
-                      &            )                                           &
+                      &            ) * rho_air_zw(k)                           &
                       &   - kmzm * ( ( v(k,j,i)   - v(k-1,j,i)   ) * ddzu(k)   &
                       &            + ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy       &
-                      &            )                                           &
-                      &   ) * ddzw(k)
+                      &            ) * rho_air_zw(k-1)                         &
+                      &   ) * ddzw(k) * drho_air(k)
              ENDDO
 
@@ -237,9 +238,9 @@
 
                 tend(k,j,i) = tend(k,j,i)                                      &
-                      & + ( kmzp * ( w(k,j,i) - w(k,j-1,i)     ) * ddy         &
-                      &   ) * ddzw(k)                                          &
-                      & + ( kmzp * ( v(k+1,j,i) - v(k,j,i)     ) * ddzu(k+1)   &
-                      &   + vsws(j,i)                                          &
-                      &   ) * ddzw(k)
+                      & + ( kmzp * ( ( v(k+1,j,i) - v(k,j,i)   ) * ddzu(k+1)   &
+                      &            + ( w(k,j,i)   - w(k,j-1,i) ) * ddy         &
+                      &            ) * rho_air_zw(k)                           &
+                      &   - ( -vsws(j,i) )                                     &
+                      &   ) * ddzw(k) * drho_air(k)
              ENDIF
 
@@ -255,9 +256,9 @@
 
                 tend(k,j,i) = tend(k,j,i)                                      &
-                      & - ( kmzm *  ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy        &
-                      &   ) * ddzw(k)                                          &
-                      & + ( -vswst(j,i)                                        &
-                      &   - kmzm * ( v(k,j,i)   - v(k-1,j,i)    ) * ddzu(k)    &
-                      &   ) * ddzw(k)
+                      & + ( ( -vswst(j,i) )                                    &
+                      &   - kmzm * ( ( v(k,j,i)   - v(k-1,j,i)   ) * ddzu(k)   &
+                      &            + ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy       &
+                      &            ) * rho_air_zw(k-1)                         &
+                      &   ) * ddzw(k) * drho_air(k)
              ENDIF
 
@@ -276,5 +277,6 @@
 
        USE arrays_3d,                                                          &
-           ONLY:  ddzu, ddzw, km, tend, u, v, vsws, vswst, w
+           ONLY:  ddzu, ddzw, km, tend, u, v, vsws, vswst, w,                  &
+                  drho_air, rho_air_zw
        
        USE control_parameters,                                                 &
@@ -385,9 +387,9 @@
                          & + ( kmzp * ( ( v(k+1,j,i) - v(k,j,i)   ) * ddzu(k+1)&
                          &            + ( w(k,j,i)   - w(k,j-1,i) ) * ddy      &
-                         &            )                                        &
+                         &            ) * rho_air_zw(k)                        &
                          &   - kmzm * ( ( v(k,j,i)   - v(k-1,j,i)   ) * ddzu(k)&
                          &            + ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy    &
-                         &            )                                        &
-                         &   ) * ddzw(k)
+                         &            ) * rho_air_zw(k-1)                      &
+                         &   ) * ddzw(k) * drho_air(k)
                 ENDIF
              ENDDO
@@ -418,9 +420,9 @@
 
                 tend(k,j,i) = tend(k,j,i)                                      &
-                      & + ( kmzp * ( w(k,j,i) - w(k,j-1,i)     ) * ddy         &
-                      &   ) * ddzw(k)                                          &
-                      & + ( kmzp * ( v(k+1,j,i) - v(k,j,i)     ) * ddzu(k+1)   &
-                      &   + vsws(j,i)                                          &
-                      &   ) * ddzw(k)
+                      & + ( kmzp * ( ( v(k+1,j,i) - v(k,j,i)   ) * ddzu(k+1)   &
+                      &            + ( w(k,j,i)   - w(k,j-1,i) ) * ddy         &
+                      &            ) * rho_air_zw(k)                           &
+                      &   - ( -vsws(j,i) )                                     &
+                      &   ) * ddzw(k) * drho_air(k)
              ENDDO
           ENDDO
@@ -444,9 +446,9 @@
 
                 tend(k,j,i) = tend(k,j,i)                                      &
-                      & - ( kmzm *  ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy        &
-                      &   ) * ddzw(k)                                          &
-                      & + ( -vswst(j,i)                                        &
-                      &   - kmzm * ( v(k,j,i)   - v(k-1,j,i)    ) * ddzu(k)    &
-                      &   ) * ddzw(k)
+                      & + ( ( -vswst(j,i) )                                    &
+                      &   - kmzm * ( ( v(k,j,i)   - v(k-1,j,i)   ) * ddzu(k)   &
+                      &            + ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy       &
+                      &            ) * rho_air_zw(k-1)                         &
+                      &   ) * ddzw(k) * drho_air(k)
              ENDDO
           ENDDO
@@ -466,5 +468,6 @@
 
        USE arrays_3d,                                                          &
-           ONLY:  ddzu, ddzw, km, tend, u, v, vsws, vswst, w
+           ONLY:  ddzu, ddzw, km, tend, u, v, vsws, vswst, w,                  &
+                  drho_air, rho_air_zw
        
        USE control_parameters,                                                 &
@@ -556,9 +559,9 @@
                       & + ( kmzp * ( ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1)     &
                       &            + ( w(k,j,i) - w(k,j-1,i) ) * ddy           &
-                      &            )                                           &
+                      &            ) * rho_air_zw(k)                           &
                       &   - kmzm * ( ( v(k,j,i)   - v(k-1,j,i)   ) * ddzu(k)   &
                       &            + ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy       &
-                      &            )                                           &
-                      &   ) * ddzw(k)
+                      &            ) * rho_air_zw(k-1)                         &
+                      &   ) * ddzw(k) * drho_air(k)
        ENDDO
 
@@ -580,9 +583,9 @@
 
           tend(k,j,i) = tend(k,j,i)                                            &
-                      & + ( kmzp * ( w(k,j,i) - w(k,j-1,i)     ) * ddy         &
-                      &   ) * ddzw(k)                                          &
-                      & + ( kmzp * ( v(k+1,j,i) - v(k,j,i)     ) * ddzu(k+1)   &
-                      &   + vsws(j,i)                                          &
-                      &   ) * ddzw(k)
+                      & + ( kmzp * ( ( v(k+1,j,i) - v(k,j,i)   ) * ddzu(k+1)   &
+                      &            + ( w(k,j,i)   - w(k,j-1,i) ) * ddy         &
+                      &            ) * rho_air_zw(k)                           &
+                      &   - ( -vsws(j,i) )                                     &
+                      &   ) * ddzw(k) * drho_air(k)
        ENDIF
 
@@ -594,13 +597,12 @@
 !
 !--       Interpolate eddy diffusivities on staggered gridpoints
-          kmzm = 0.25_wp * &
-                 ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) )
+          kmzm = 0.25_wp * ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) )
 
           tend(k,j,i) = tend(k,j,i)                                            &
-                      & - ( kmzm *  ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy        &
-                      &   ) * ddzw(k)                                          &
-                      & + ( -vswst(j,i)                                        &
-                      &   - kmzm * ( v(k,j,i)   - v(k-1,j,i)    ) * ddzu(k)    &
-                      &   ) * ddzw(k)
+                      & + ( ( -vswst(j,i) )                                    &
+                      &   - kmzm * ( ( v(k,j,i)   - v(k-1,j,i)   ) * ddzu(k)   &
+                      &            + ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy       &
+                      &            ) * rho_air_zw(k-1)                         &
+                      &   ) * ddzw(k) * drho_air(k)
        ENDIF
 
Index: /palm/trunk/SOURCE/diffusion_w.f90
===================================================================
--- /palm/trunk/SOURCE/diffusion_w.f90	(revision 2036)
+++ /palm/trunk/SOURCE/diffusion_w.f90	(revision 2037)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! Anelastic approximation implemented
 ! 
 ! Former revisions:
@@ -119,5 +119,5 @@
 
        USE arrays_3d,                                                          &          
-           ONLY :  ddzu, ddzw, km, tend, u, v, w
+           ONLY :  ddzu, ddzw, km, tend, u, v, w, drho_air_zw, rho_air
            
        USE control_parameters,                                                 & 
@@ -184,6 +184,8 @@
                       & + 2.0_wp * (                                           &
                       &   km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1)  &
+                      &               * rho_air(k+1)                           &
                       & - km(k,j,i)   * ( w(k,j,i)   - w(k-1,j,i) ) * ddzw(k)  &
-                      &            ) * ddzu(k+1)
+                      &               * rho_air(k)                             &
+                      &            ) * ddzu(k+1) * drho_air_zw(k)
              ENDDO
 
@@ -227,6 +229,8 @@
                                  + 2.0_wp * (                                  &
                            km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1) &
+                                       * rho_air(k+1)                          &
                          - km(k,j,i)   * ( w(k,j,i)   - w(k-1,j,i) ) * ddzw(k) &
-                                            ) * ddzu(k+1)
+                                       * rho_air(k)                            &
+                                            ) * ddzu(k+1) * drho_air_zw(k)
                 ENDDO
              ENDIF
@@ -246,5 +250,5 @@
 
        USE arrays_3d,                                                          &
-           ONLY :  ddzu, ddzw, km, tend, u, v, w
+           ONLY :  ddzu, ddzw, km, tend, u, v, w, drho_air_zw, rho_air
            
        USE control_parameters,                                                 &
@@ -316,6 +320,8 @@
                          & + 2.0_wp * (                                          &
                          &   km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1) &
+                         &               * rho_air(k+1)                          &
                          & - km(k,j,i)   * ( w(k,j,i)   - w(k-1,j,i) ) * ddzw(k) &
-                         &            ) * ddzu(k+1)
+                         &               * rho_air(k)                            &
+                         &            ) * ddzu(k+1) * drho_air_zw(k)
                 ENDIF
              ENDDO
@@ -361,6 +367,8 @@
                                  + 2.0_wp * (                                  &
                            km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1) &
+                                       * rho_air(k+1)                          &
                          - km(k,j,i)   * ( w(k,j,i)   - w(k-1,j,i) ) * ddzw(k) &
-                                            ) * ddzu(k+1)
+                                       * rho_air(k)                            &
+                                            ) * ddzu(k+1) * drho_air_zw(k)
                 ENDIF
              ENDDO
@@ -381,5 +389,5 @@
 
        USE arrays_3d,                                                          &          
-           ONLY :  ddzu, ddzw, km, tend, u, v, w
+           ONLY :  ddzu, ddzw, km, tend, u, v, w, drho_air_zw, rho_air
            
        USE control_parameters,                                                 & 
@@ -430,6 +438,8 @@
                       & + 2.0_wp * (                                           &
                       &   km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1)  &
+                      &               * rho_air(k+1)                           &
                       & - km(k,j,i)   * ( w(k,j,i)   - w(k-1,j,i) ) * ddzw(k)  &
-                      &            ) * ddzu(k+1)
+                      &               * rho_air(k)                             &
+                      &            ) * ddzu(k+1) * drho_air_zw(k)
        ENDDO
 
@@ -485,6 +495,8 @@
                                  + 2.0_wp * (                                  &
                            km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1) &
+                                       * rho_air(k+1)                          &
                          - km(k,j,i)   * ( w(k,j,i)   - w(k-1,j,i) ) * ddzw(k) &
-                                            ) * ddzu(k+1)
+                                       * rho_air(k)                            &
+                                            ) * ddzu(k+1) * drho_air_zw(k)
           ENDDO
        ENDIF
Index: /palm/trunk/SOURCE/disturb_heatflux.f90
===================================================================
--- /palm/trunk/SOURCE/disturb_heatflux.f90	(revision 2036)
+++ /palm/trunk/SOURCE/disturb_heatflux.f90	(revision 2037)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! ------------------
-! 
+! Anelastic approximation implemented
 ! 
 ! Former revisions:
@@ -65,5 +65,5 @@
 
     USE arrays_3d,                                                             &
-        ONLY:  shf
+        ONLY:  shf, heatflux_input_conversion
         
     USE control_parameters,                                                    &
@@ -76,5 +76,5 @@
     
     USE indices,                                                               &
-        ONLY:  nx, nxl, nxr, ny, nyn, nys, nzb_s_inner
+        ONLY:  nx, nxl, nxr, ny, nyn, nys, nzb, nzb_s_inner
 
     IMPLICIT NONE
@@ -97,9 +97,11 @@
           THEN
              IF ( nzb_s_inner(j,i) == 0 )  THEN
-                shf(j,i) = randomnumber * surface_heatflux
+                shf(j,i) = randomnumber * surface_heatflux                     &
+                           * heatflux_input_conversion(nzb)
              ELSE
 !
 !--             Over topography surface_heatflux is replaced by wall_heatflux(0)
-                shf(j,i) = randomnumber * wall_heatflux(0)
+                shf(j,i) = randomnumber * wall_heatflux(0)                     &
+                           * heatflux_input_conversion(nzb_s_inner(j,i))
              ENDIF
           ENDIF
Index: /palm/trunk/SOURCE/flow_statistics.f90
===================================================================
--- /palm/trunk/SOURCE/flow_statistics.f90	(revision 2036)
+++ /palm/trunk/SOURCE/flow_statistics.f90	(revision 2037)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! Anelastic approximation implemented
 ! 
 ! Former revisions:
@@ -214,9 +214,11 @@
 
     USE arrays_3d,                                                             &
-        ONLY:  ddzu, ddzw, e, hyp, km, kh, nr, ol, p, prho, prr, pt, q, qc, ql,&
-               qr, qs, qsws, qswst, rho_ocean, s, sa, ss, ssws, sswst, saswsb,       &
-               saswst, shf, td_lsa_lpt, td_lsa_q, td_sub_lpt, td_sub_q,        &
-               time_vert, ts, tswst, u, ug, us, usws, uswst, vsws, v, vg, vpt, &
-               vswst, w, w_subs, zw
+        ONLY:  ddzu, ddzw, e, heatflux_output_conversion, hyp, km, kh,         &
+               momentumflux_output_conversion, nr, ol, p, prho, prr, pt, q,    &
+               qc, ql, qr, qs, qsws, qswst, rho_air, rho_air_zw, rho_ocean, s, &
+               sa, ss, ssws, sswst, saswsb, saswst, shf, td_lsa_lpt, td_lsa_q, &
+               td_sub_lpt, td_sub_q, time_vert, ts, tswst, u, ug, us, usws,    &
+               uswst, vsws, v, vg, vpt, vswst, w, w_subs,                      &
+               waterflux_output_conversion, zw
         
     USE cloud_parameters,                                                      &
@@ -326,5 +328,6 @@
        sums_l(:,11,:) = sums_l_l(:,sr,:)      ! mixing length from diffusivities
 !--    WARNING: next line still has to be adjusted for OpenMP 
-       sums_l(:,21,0) = sums_wsts_bc_l(:,sr)  ! heat flux from advec_s_bc
+       sums_l(:,21,0) = sums_wsts_bc_l(:,sr) *                                 &
+                        heatflux_output_conversion  ! heat flux from advec_s_bc
        sums_l(nzb+9,pr_palm,0)  = sums_divold_l(sr)  ! old divergence from pres
        sums_l(nzb+10,pr_palm,0) = sums_divnew_l(sr)  ! new divergence from pres
@@ -358,6 +361,8 @@
 !
 !--          Swap the turbulent quantities evaluated in advec_ws.
-             sums_l(:,13,i) = sums_wsus_ws_l(:,i)       ! w*u* 
-             sums_l(:,15,i) = sums_wsvs_ws_l(:,i)       ! w*v* 
+             sums_l(:,13,i) = sums_wsus_ws_l(:,i)                              &
+                              * momentumflux_output_conversion ! w*u*
+             sums_l(:,15,i) = sums_wsvs_ws_l(:,i)                              &
+                              * momentumflux_output_conversion ! w*v*
              sums_l(:,30,i) = sums_us2_ws_l(:,i)        ! u*2 
              sums_l(:,31,i) = sums_vs2_ws_l(:,i)        ! v*2 
@@ -373,7 +378,9 @@
 
           DO  i = 0, threads_per_task-1
-             sums_l(:,17,i)                        = sums_wspts_ws_l(:,i) ! w*pt* 
+             sums_l(:,17,i)                        = sums_wspts_ws_l(:,i)      &
+                                           * heatflux_output_conversion  ! w*pt*
              IF ( ocean          ) sums_l(:,66,i)  = sums_wssas_ws_l(:,i) ! w*sa*
-             IF ( humidity       ) sums_l(:,49,i)  = sums_wsqs_ws_l(:,i)  ! w*q*
+             IF ( humidity       ) sums_l(:,49,i)  = sums_wsqs_ws_l(:,i)       &
+                                           * waterflux_output_conversion  ! w*q*
              IF ( passive_scalar ) sums_l(:,116,i) = sums_wsss_ws_l(:,i)  ! w*s*
           ENDDO
@@ -722,5 +729,7 @@
                                    ( u(k+1,j,i) - u(k,j,i)   ) * ddzu(k+1)     &
                                  + ( w(k,j,i)   - w(k,j,i-1) ) * ddx           &
-                                                               ) * rmask(j,i,sr)
+                                                           ) * rmask(j,i,sr)   &
+                                         * rho_air_zw(k)                       &
+                                         * momentumflux_output_conversion(k)
 !
 !--             Momentum flux w"v"
@@ -730,5 +739,7 @@
                                    ( v(k+1,j,i) - v(k,j,i)   ) * ddzu(k+1)     &
                                  + ( w(k,j,i)   - w(k,j-1,i) ) * ddy           &
-                                                               ) * rmask(j,i,sr)
+                                                           ) * rmask(j,i,sr)   &
+                                         * rho_air_zw(k)                       &
+                                         * momentumflux_output_conversion(k)
 !
 !--             Heat flux w"pt"
@@ -736,4 +747,6 @@
                                          - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )&
                                                * ( pt(k+1,j,i) - pt(k,j,i) )   &
+                                               * rho_air_zw(k)                 &
+                                               * heatflux_output_conversion(k) &
                                                * ddzu(k+1) * rmask(j,i,sr)
 
@@ -754,8 +767,12 @@
                                          - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )&
                                                * ( vpt(k+1,j,i) - vpt(k,j,i) ) &
+                                               * rho_air_zw(k)                 &
+                                               * heatflux_output_conversion(k) &
                                                * ddzu(k+1) * rmask(j,i,sr)
                    sums_l(k,48,tn) = sums_l(k,48,tn)                           &
                                          - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )&
                                                * ( q(k+1,j,i) - q(k,j,i) )     &
+                                               * rho_air_zw(k)                 &
+                                               * waterflux_output_conversion(k)&
                                                * ddzu(k+1) * rmask(j,i,sr)
 
@@ -765,4 +782,6 @@
                                                * ( ( q(k+1,j,i) - ql(k+1,j,i) )&
                                                 - ( q(k,j,i) - ql(k,j,i) ) )   &
+                                               * rho_air_zw(k)                 &
+                                               * waterflux_output_conversion(k)&
                                                * ddzu(k+1) * rmask(j,i,sr) 
                    ENDIF
@@ -784,8 +803,11 @@
              IF ( use_surface_fluxes )  THEN
                 sums_l(nzb,12,tn) = sums_l(nzb,12,tn) + &
+                                    momentumflux_output_conversion(nzb) * &
                                     usws(j,i) * rmask(j,i,sr)     ! w"u"
                 sums_l(nzb,14,tn) = sums_l(nzb,14,tn) + &
+                                    momentumflux_output_conversion(nzb) * &
                                     vsws(j,i) * rmask(j,i,sr)     ! w"v"
                 sums_l(nzb,16,tn) = sums_l(nzb,16,tn) + &
+                                    heatflux_output_conversion(nzb) * &
                                     shf(j,i)  * rmask(j,i,sr)     ! w"pt"
                 sums_l(nzb,58,tn) = sums_l(nzb,58,tn) + &
@@ -799,19 +821,23 @@
                 IF ( humidity )  THEN
                    sums_l(nzb,48,tn) = sums_l(nzb,48,tn) +                     &
+                                       waterflux_output_conversion(nzb) *      &
                                        qsws(j,i) * rmask(j,i,sr)  ! w"q" (w"qv")
                    sums_l(nzb,45,tn) = sums_l(nzb,45,tn) + (                   &
                                        ( 1.0_wp + 0.61_wp * q(nzb,j,i) ) *     &
                                        shf(j,i) + 0.61_wp * pt(nzb,j,i) *      &
-                                                  qsws(j,i) )
+                                                  qsws(j,i) )                  &
+                                       * heatflux_output_conversion(nzb)
                    IF ( cloud_droplets )  THEN
                       sums_l(nzb,45,tn) = sums_l(nzb,45,tn) + (                &
                                          ( 1.0_wp + 0.61_wp * q(nzb,j,i) -     &
                                            ql(nzb,j,i) ) * shf(j,i) +          &
-                                           0.61_wp * pt(nzb,j,i) * qsws(j,i) )
+                                           0.61_wp * pt(nzb,j,i) * qsws(j,i) ) &
+                                          * heatflux_output_conversion(nzb)
                    ENDIF
                    IF ( cloud_physics )  THEN
 !
 !--                   Formula does not work if ql(nzb) /= 0.0
-                      sums_l(nzb,51,tn) = sums_l(nzb,51,tn) +                  & 
+                      sums_l(nzb,51,tn) = sums_l(nzb,51,tn) +                  &
+                                          waterflux_output_conversion(nzb) *   &
                                           qsws(j,i) * rmask(j,i,sr) ! w"q" (w"qv")
                    ENDIF
@@ -860,8 +886,11 @@
              IF ( use_top_fluxes )  THEN
                 sums_l(nzt:nzt+1,12,tn) = sums_l(nzt:nzt+1,12,tn) + &
+                                    momentumflux_output_conversion(nzt:nzt+1) * &
                                     uswst(j,i) * rmask(j,i,sr)    ! w"u"
                 sums_l(nzt:nzt+1,14,tn) = sums_l(nzt:nzt+1,14,tn) + &
+                                    momentumflux_output_conversion(nzt:nzt+1) * &
                                     vswst(j,i) * rmask(j,i,sr)    ! w"v"
                 sums_l(nzt:nzt+1,16,tn) = sums_l(nzt:nzt+1,16,tn) + &
+                                    heatflux_output_conversion(nzt:nzt+1) * &
                                     tswst(j,i)  * rmask(j,i,sr)   ! w"pt"
                 sums_l(nzt:nzt+1,58,tn) = sums_l(nzt:nzt+1,58,tn) + &
@@ -876,14 +905,17 @@
                 IF ( humidity )  THEN
                    sums_l(nzt,48,tn) = sums_l(nzt,48,tn) +                     &
+                                       waterflux_output_conversion(nzt) *      &
                                        qswst(j,i) * rmask(j,i,sr) ! w"q" (w"qv")
                    sums_l(nzt,45,tn) = sums_l(nzt,45,tn) + (                   &
                                        ( 1.0_wp + 0.61_wp * q(nzt,j,i) ) *     &
                                        tswst(j,i) + 0.61_wp * pt(nzt,j,i) *    &
-                                                             qswst(j,i) )
+                                                             qswst(j,i) )      &
+                                       * heatflux_output_conversion(nzt)
                    IF ( cloud_droplets )  THEN
                       sums_l(nzt,45,tn) = sums_l(nzt,45,tn) + (                &
                                           ( 1.0_wp + 0.61_wp * q(nzt,j,i) -    &
                                             ql(nzt,j,i) ) * tswst(j,i) +       &
-                                            0.61_wp * pt(nzt,j,i) * qswst(j,i) )
+                                           0.61_wp * pt(nzt,j,i) * qswst(j,i) )&
+                                           * heatflux_output_conversion(nzt)
                    ENDIF
                    IF ( cloud_physics )  THEN
@@ -891,4 +923,5 @@
 !--                   Formula does not work if ql(nzb) /= 0.0
                       sums_l(nzt,51,tn) = sums_l(nzt,51,tn) + &   ! w"q" (w"qv")
+                                          waterflux_output_conversion(nzt) *   &
                                           qswst(j,i) * rmask(j,i,sr)
                    ENDIF
@@ -943,4 +976,5 @@
                                     vpt(k+1,j,i) - hom(k+1,1,44,sr) )
                       sums_l(k,46,tn) = sums_l(k,46,tn) + pts * w(k,j,i) *     &
+                                               heatflux_output_conversion(k) * &
                                                           rmask(j,i,sr)
                       sums_l(k,54,tn) = sums_l(k,54,tn) + ql(k,j,i) * rmask(j,i,sr)
@@ -953,4 +987,5 @@
                               hom(k+1,1,42,sr) )
                          sums_l(k,52,tn) = sums_l(k,52,tn) + pts * w(k,j,i) *  &
+                                             waterflux_output_conversion(k) *  &
                                                              rmask(j,i,sr)
                          sums_l(k,75,tn) = sums_l(k,75,tn) + qc(k,j,i) *       &
@@ -971,11 +1006,13 @@
                                           vpt(k+1,j,i) - hom(k+1,1,44,sr) )
                          sums_l(k,46,tn) = sums_l(k,46,tn) + pts * w(k,j,i) *  &
+                                              heatflux_output_conversion(k) *  &
                                                              rmask(j,i,sr)
                       ELSE IF ( ws_scheme_sca .AND. sr == 0 )  THEN
-                         sums_l(k,46,tn) = ( 1.0_wp + 0.61_wp *                & 
-                                             hom(k,1,41,sr) ) *                &
-                                           sums_l(k,17,tn) +                   &
-                                           0.61_wp * hom(k,1,4,sr) *           &
-                                           sums_l(k,49,tn)
+                         sums_l(k,46,tn) = ( ( 1.0_wp + 0.61_wp *              & 
+                                               hom(k,1,41,sr) ) *              &
+                                             sums_l(k,17,tn) +                 &
+                                             0.61_wp * hom(k,1,4,sr) *         &
+                                             sums_l(k,49,tn)                   &
+                                           ) * heatflux_output_conversion(k)
                       END IF
                    END IF
@@ -996,4 +1033,5 @@
                 sums_l(k,37,tn) = sums_l(k,37,tn) + w(k,j,i) * 0.5_wp *        &
                                              ( ust**2 + vst**2 + w(k,j,i)**2 ) &
+                                           * momentumflux_output_conversion(k) &
                                              * rmask(j,i,sr)
              ENDDO
@@ -1017,4 +1055,5 @@
                   sums_l(k,13,tn) = sums_l(k,13,tn) + 0.5_wp *                 &
                                                     ( w(k,j,i-1) + w(k,j,i) )  &
+                                          * momentumflux_output_conversion(k)  &
                                                     * ust * rmask(j,i,sr)
 !
@@ -1022,4 +1061,5 @@
                   sums_l(k,15,tn) = sums_l(k,15,tn) + 0.5_wp *                 &
                                                     ( w(k,j-1,i) + w(k,j,i) )  &
+                                          * momentumflux_output_conversion(k)  &
                                                     * vst * rmask(j,i,sr)
                ENDDO
@@ -1038,4 +1078,5 @@
                            ( pt(k,j,i)   - hom(k,1,4,sr) +                     &
                              pt(k+1,j,i) - hom(k+1,1,4,sr) )                   &
+                           * heatflux_output_conversion(k)                     &
                            * w(k,j,i) * rmask(j,i,sr)
                   IF ( humidity )  THEN
@@ -1043,4 +1084,5 @@
                                       q(k+1,j,i) - hom(k+1,1,41,sr) )
                      sums_l(k,49,tn) = sums_l(k,49,tn) + pts * w(k,j,i) *      &
+                                       waterflux_output_conversion(k) *        &
                                        rmask(j,i,sr)
                   ENDIF
@@ -1147,8 +1189,12 @@
                                                    ( kh(k,j,i) + kh(k,j,i-1) ) &
                                                  * ( pt(k,j,i-1) - pt(k,j,i) ) &
+                                               * rho_air_zw(k)                 &
+                                               * heatflux_output_conversion(k) &
                                                  * ddx * rmask(j,i,sr)
                    sums_l(k,61,tn) = sums_l(k,61,tn) - 0.5_wp *                &
                                                    ( kh(k,j,i) + kh(k,j-1,i) ) &
                                                  * ( pt(k,j-1,i) - pt(k,j,i) ) &
+                                               * rho_air_zw(k)                 &
+                                               * heatflux_output_conversion(k) &
                                                  * ddy * rmask(j,i,sr)
 !
@@ -1157,5 +1203,6 @@
                                                   ( u(k,j,i) - hom(k,1,1,sr) ) &
                                     * 0.5_wp * ( pt(k,j,i-1) - hom(k,1,4,sr) + &
-                                                 pt(k,j,i)   - hom(k,1,4,sr) )
+                                                 pt(k,j,i)   - hom(k,1,4,sr) ) &
+                                               * heatflux_output_conversion(k)
                    pts = 0.5_wp * ( pt(k,j-1,i) - hom(k,1,4,sr) +              &
                                     pt(k,j,i)   - hom(k,1,4,sr) )
@@ -1163,5 +1210,6 @@
                                                   ( v(k,j,i) - hom(k,1,2,sr) ) &
                                     * 0.5_wp * ( pt(k,j-1,i) - hom(k,1,4,sr) + &
-                                                 pt(k,j,i)   - hom(k,1,4,sr) )
+                                                 pt(k,j,i)   - hom(k,1,4,sr) ) &
+                                               * heatflux_output_conversion(k)
                 ENDDO
              ENDDO
@@ -1489,4 +1537,7 @@
        ENDIF
 
+       hom(:,1,121,sr) = rho_air       ! rho_air in Kg/m^3
+       hom(:,1,122,sr) = rho_air_zw    ! rho_air_zw in Kg/m^3
+
        hom(:,1,pr_palm,sr) =   sums(:,pr_palm)
                                        ! u*, w'u', w'v', t* (in last profile)
@@ -1597,5 +1648,6 @@
        THEN
           hom(nzb+8,1,pr_palm,sr) = &
-             ( g / hom(k_surface_level+1,1,4,sr) * hom(k_surface_level,1,18,sr)&
+             ( g / hom(k_surface_level+1,1,4,sr) *                             &
+             ( hom(k_surface_level,1,18,sr) / heatflux_output_conversion(nzb) )&
              * ABS( z_i(1) - mean_surface_level_height(sr) ) )**0.333333333_wp
        ELSE
@@ -1708,9 +1760,10 @@
 
     USE arrays_3d,                                                             &
-        ONLY:  ddzu, ddzw, e, hyp, km, kh, nr, p, prho, pt, q, qc, ql, qr, qs, &
-               qsws, qswst, rho_ocean, s, sa, saswsb, saswst, shf, ss, ssws, sswst,  &
-               td_lsa_lpt, td_lsa_q, td_sub_lpt, td_sub_q, time_vert, ts,      &
-               tswst, u, ug, us, usws, uswst, vsws, v, vg, vpt, vswst, w,      &
-               w_subs, zw
+        ONLY:  ddzu, ddzw, e, heatflux_output_conversion, hyp, km, kh,         &
+               momentumflux_output_conversion, nr, p, prho, pt, q, qc, ql, qr, &
+               qs, qsws, qswst, rho_air, rho_air_zw, rho_ocean, s, sa, saswsb, &
+               saswst, shf, ss, ssws, sswst, td_lsa_lpt, td_lsa_q, td_sub_lpt, &
+               td_sub_q, time_vert, ts, tswst, u, ug, us, usws, uswst, vsws,   &
+               v, vg, vpt, vswst, w, w_subs, waterflux_output_conversion, zw
                  
         
@@ -1828,5 +1881,6 @@
        sums_l(:,11,:) = sums_l_l(:,sr,:)      ! mixing length from diffusivities
 !--    WARNING: next line still has to be adjusted for OpenMP
-       sums_l(:,21,0) = sums_wsts_bc_l(:,sr)  ! heat flux from advec_s_bc
+       sums_l(:,21,0) = sums_wsts_bc_l(:,sr) *                                 &
+                        heatflux_output_conversion  ! heat flux from advec_s_bc
        sums_l(nzb+9,pr_palm,0)  = sums_divold_l(sr)  ! old divergence from pres
        sums_l(nzb+10,pr_palm,0) = sums_divnew_l(sr)  ! new divergence from pres
@@ -1860,6 +1914,8 @@
 !
 !--          Swap the turbulent quantities evaluated in advec_ws.
-             sums_l(:,13,i) = sums_wsus_ws_l(:,i)       ! w*u*
-             sums_l(:,15,i) = sums_wsvs_ws_l(:,i)       ! w*v*
+             sums_l(:,13,i) = sums_wsus_ws_l(:,i)                              &
+                              * momentumflux_output_conversion ! w*u*
+             sums_l(:,15,i) = sums_wsvs_ws_l(:,i)                              &
+                              * momentumflux_output_conversion ! w*v*
              sums_l(:,30,i) = sums_us2_ws_l(:,i)        ! u*2
              sums_l(:,31,i) = sums_vs2_ws_l(:,i)        ! v*2
@@ -1881,7 +1937,9 @@
 
           DO  i = 0, threads_per_task-1
-             sums_l(:,17,i) = sums_wspts_ws_l(:,i)      ! w*pt* from advec_s_ws
+             sums_l(:,17,i) = sums_wspts_ws_l(:,i)                             &
+                              * heatflux_output_conversion        ! w*pt* from advec_s_ws
              IF ( ocean ) sums_l(:,66,i) = sums_wssas_ws_l(:,i) ! w*sa*
-             IF ( humidity       )  sums_l(:,49,i)  = sums_wsqs_ws_l(:,i) !w*q*
+             IF ( humidity       )  sums_l(:,49,i)  = sums_wsqs_ws_l(:,i)      &
+                                            * waterflux_output_conversion !w*q*
              IF ( passive_scalar )  sums_l(:,116,i) = sums_wsss_ws_l(:,i) !w*s*
           ENDDO
@@ -2378,5 +2436,7 @@
                                  + ( w(k,j,i)   - w(k,j,i-1) ) * ddx           &
                                                                )               &
-                               * rmask(j,i,sr) * rflags_invers(j,i,k+1)
+                               * rmask(j,i,sr) * rflags_invers(j,i,k+1)        &
+                               * rho_air_zw(k)                                 &
+                               * momentumflux_output_conversion(k)
 !
 !--             Momentum flux w"v"
@@ -2387,9 +2447,13 @@
                                  + ( w(k,j,i)   - w(k,j-1,i) ) * ddy           &
                                                                )               &
-                               * rmask(j,i,sr) * rflags_invers(j,i,k+1)
+                               * rmask(j,i,sr) * rflags_invers(j,i,k+1)        &
+                               * rho_air_zw(k)                                 &
+                               * momentumflux_output_conversion(k)
 !
 !--             Heat flux w"pt"
                 s3 = s3 - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )                 &
                                  * ( pt(k+1,j,i) - pt(k,j,i) )                 &
+                                 * rho_air_zw(k)                               &
+                                 * heatflux_output_conversion(k)               &
                                  * ddzu(k+1) * rmask(j,i,sr)                   &
                                  * rflags_invers(j,i,k+1)
@@ -2435,8 +2499,12 @@
                    s1 = s1 - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )              &
                                     * ( vpt(k+1,j,i) - vpt(k,j,i) )            &
+                                    * rho_air_zw(k)                            &
+                                    * heatflux_output_conversion(k)            &
                                     * ddzu(k+1) * rmask(j,i,sr)                &
                                     * rflags_invers(j,i,k+1)
                    s2 = s2 - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )              &
                                     * ( q(k+1,j,i) - q(k,j,i) )                &
+                                    * rho_air_zw(k)                            &
+                                    * waterflux_output_conversion(k)           &
                                     * ddzu(k+1) * rmask(j,i,sr)                &
                                     * rflags_invers(j,i,k+1)
@@ -2459,4 +2527,6 @@
                                        *  ( ( q(k+1,j,i) - ql(k+1,j,i) )       &
                                           - ( q(k,j,i) - ql(k,j,i) ) )         &
+                                       * rho_air_zw(k)                         &
+                                       * waterflux_output_conversion(k)        &
                                        * ddzu(k+1) * rmask(j,i,sr)             & 
                                        * rflags_invers(j,i,k+1)
@@ -2506,7 +2576,10 @@
 !
 !--             Subgridscale fluxes in the Prandtl layer
-                s1 = s1 + usws(j,i) * rmask(j,i,sr)     ! w"u"
-                s2 = s2 + vsws(j,i) * rmask(j,i,sr)     ! w"v"
-                s3 = s3 + shf(j,i)  * rmask(j,i,sr)     ! w"pt"
+                s1 = s1 + usws(j,i) * momentumflux_output_conversion(nzb)      &
+                                    * rmask(j,i,sr) ! w"u"
+                s2 = s2 + vsws(j,i) * momentumflux_output_conversion(nzb)      &
+                                    * rmask(j,i,sr) ! w"v"
+                s3 = s3 + shf(j,i)  * heatflux_output_conversion(nzb)          &
+                                    * rmask(j,i,sr) ! w"pt"
                 s4 = s4 + 0.0_wp * rmask(j,i,sr)        ! u"pt"
                 s5 = s5 + 0.0_wp * rmask(j,i,sr)        ! v"pt"
@@ -2545,7 +2618,9 @@
              DO  i = nxl, nxr
                 DO  j =  nys, nyn
-                   s1 = s1 + qsws(j,i) * rmask(j,i,sr)  ! w"q" (w"qv")
+                   s1 = s1 + qsws(j,i) * waterflux_output_conversion(nzb)      &
+                                       * rmask(j,i,sr) ! w"q" (w"qv")
                    s2 = s2 + ( ( 1.0_wp + 0.61_wp * q(nzb,j,i) ) * shf(j,i)    &
-                               + 0.61_wp * pt(nzb,j,i) * qsws(j,i) )
+                               + 0.61_wp * pt(nzb,j,i) * qsws(j,i) )           &
+                             * heatflux_output_conversion(nzb)
                 ENDDO
              ENDDO
@@ -2564,5 +2639,6 @@
                       s1 = s1 + ( ( 1.0_wp +                                   &
                                     0.61_wp * q(nzb,j,i) - ql(nzb,j,i) ) *     &
-                                  shf(j,i) + 0.61_wp * pt(nzb,j,i) * qsws(j,i) )
+                                 shf(j,i) + 0.61_wp * pt(nzb,j,i) * qsws(j,i) )&
+                                * heatflux_output_conversion(nzb)
                    ENDDO
                 ENDDO
@@ -2582,5 +2658,6 @@
 !
 !--                   Formula does not work if ql(nzb) /= 0.0
-                      s1 = s1 + qsws(j,i) * rmask(j,i,sr)   ! w"q" (w"qv")
+                      s1 = s1 + qsws(j,i) * waterflux_output_conversion(nzb)   &
+                                          * rmask(j,i,sr)   ! w"q" (w"qv")
                    ENDDO
                 ENDDO
@@ -2624,7 +2701,10 @@
           DO  i = nxl, nxr
              DO  j =  nys, nyn
-                s1 = s1 + uswst(j,i) * rmask(j,i,sr)    ! w"u"
-                s2 = s2 + vswst(j,i) * rmask(j,i,sr)    ! w"v"
-                s3 = s3 + tswst(j,i)  * rmask(j,i,sr)   ! w"pt"
+                s1 = s1 + uswst(j,i) * momentumflux_output_conversion(nzt:nzt+1) &
+                                     * rmask(j,i,sr)    ! w"u"
+                s2 = s2 + vswst(j,i) * momentumflux_output_conversion(nzt:nzt+1) &
+                                     * rmask(j,i,sr)    ! w"v"
+                s3 = s3 + tswst(j,i) * heatflux_output_conversion(nzt:nzt+1)   &
+                                     * rmask(j,i,sr)    ! w"pt"
                 s4 = s4 + 0.0_wp * rmask(j,i,sr)        ! u"pt"
                 s5 = s5 + 0.0_wp * rmask(j,i,sr)        ! v"pt"
@@ -2663,7 +2743,9 @@
              DO  i = nxl, nxr
                 DO  j =  nys, nyn
-                   s1 = s1 + qswst(j,i) * rmask(j,i,sr) ! w"q" (w"qv")
+                   s1 = s1 + qswst(j,i) * waterflux_output_conversion(nzt)     &
+                                        * rmask(j,i,sr) ! w"q" (w"qv")
                    s2 = s2 + ( ( 1.0_wp + 0.61_wp * q(nzt,j,i) ) * tswst(j,i) +&
-                                 0.61_wp * pt(nzt,j,i) * qswst(j,i) )
+                                 0.61_wp * pt(nzt,j,i) * qswst(j,i) )          &
+                             * heatflux_output_conversion(nzt)
                 ENDDO
              ENDDO
@@ -2683,5 +2765,6 @@
                                     0.61_wp * q(nzt,j,i) - ql(nzt,j,i) ) *     &
                                   tswst(j,i) +                                 &
-                                  0.61_wp * pt(nzt,j,i) * qswst(j,i) )
+                                  0.61_wp * pt(nzt,j,i) * qswst(j,i) )         &
+                                * heatflux_output_conversion(nzt)
                    ENDDO
                 ENDDO
@@ -2701,5 +2784,6 @@
 !
 !--                   Formula does not work if ql(nzb) /= 0.0
-                      s1 = s1 + qswst(j,i) * rmask(j,i,sr)  ! w"q" (w"qv")
+                      s1 = s1 + qswst(j,i) * waterflux_output_conversion(nzt)  &
+                                           * rmask(j,i,sr)  ! w"q" (w"qv")
                    ENDDO
                 ENDDO
@@ -2755,5 +2839,6 @@
 !--             Energy flux w*e* (has to be adjusted?)
                 s3 = s3 + w(k,j,i) * 0.5_wp * ( ust**2 + vst**2 + w(k,j,i)**2 )&
-                                   * rmask(j,i,sr) * rflags_invers(j,i,k+1)
+                                   * rmask(j,i,sr) * rflags_invers(j,i,k+1)    &
+                                   * momentumflux_output_conversion(k)
              ENDDO
           ENDDO
@@ -2822,4 +2907,5 @@
                       s1 = s1 + 0.5_wp * ( vpt(k,j,i)   - hom(k,1,44,sr) +     &
                                            vpt(k+1,j,i) - hom(k+1,1,44,sr) ) * &
+                                         heatflux_output_conversion(k) *       &
                                          w(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
                    ENDDO
@@ -2839,4 +2925,5 @@
                          s1 = s1 + 0.5_wp * ( ( q(k,j,i)   - ql(k,j,i)   ) - hom(k,1,42,sr) +   &
                                               ( q(k+1,j,i) - ql(k+1,j,i) ) - hom(k+1,1,42,sr) ) &
+                                          * waterflux_output_conversion(k)                      &
                                           * w(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
                       ENDDO
@@ -2928,4 +3015,5 @@
                          s1 = s1 + 0.5_wp * ( vpt(k,j,i)   - hom(k,1,44,sr) +   &
                                               vpt(k+1,j,i) - hom(k+1,1,44,sr) ) &
+                                          * heatflux_output_conversion(k)       &
                                           * w(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
                       ENDDO
@@ -2939,6 +3027,8 @@
                 !$acc parallel loop present( hom, sums_l )
                 DO  k = nzb, nzt_diff
-                   sums_l(k,46,tn) = ( 1.0_wp + 0.61_wp * hom(k,1,41,sr) ) * sums_l(k,17,tn) + &
-                                                0.61_wp * hom(k,1,4,sr) * sums_l(k,49,tn)
+                   sums_l(k,46,tn) = ( ( 1.0_wp + 0.61_wp * hom(k,1,41,sr) ) * &
+                                       sums_l(k,17,tn) + 0.61_wp *             &
+                                       hom(k,1,4,sr) * sums_l(k,49,tn)         &
+                                     ) * heatflux_output_conversion(k)
                 ENDDO
                 !$acc end parallel loop
@@ -2993,4 +3083,5 @@
                    s1 = s1 + 0.5_wp * ( w(k,j,i-1) + w(k,j,i) )                &
                                     * ust * rmask(j,i,sr)                      &
+                                    * momentumflux_output_conversion(k)        &
                                     * rflags_invers(j,i,k+1)
 !
@@ -2998,9 +3089,10 @@
                    s2 = s2 + 0.5_wp * ( w(k,j-1,i) + w(k,j,i) )                &
                                     * vst * rmask(j,i,sr)                      & 
+                                    * momentumflux_output_conversion(k)        &
                                     * rflags_invers(j,i,k+1)
                 ENDDO
              ENDDO
              sums_l(k,13,tn) = s1
-             sums_l(k,15,tn) = s1
+             sums_l(k,15,tn) = s2
           ENDDO
           !$acc end parallel loop
@@ -3021,4 +3113,5 @@
                    s1 = s1 + 0.5_wp * ( pt(k,j,i)   - hom(k,1,4,sr) +          &
                                         pt(k+1,j,i) - hom(k+1,1,4,sr) )        &
+                                    * heatflux_output_conversion(k)            &
                                     * w(k,j,i) * rmask(j,i,sr)                 & 
                                     * rflags_invers(j,i,k+1)
@@ -3039,4 +3132,5 @@
                       s1 = s1 + 0.5_wp * ( q(k,j,i)   - hom(k,1,41,sr) +       &
                                            q(k+1,j,i) - hom(k+1,1,41,sr) )     &
+                                       * waterflux_output_conversion(k)        &
                                        * w(k,j,i) * rmask(j,i,sr)              &
                                        * rflags_invers(j,i,k+1)
@@ -3167,8 +3261,12 @@
                                                    ( kh(k,j,i) + kh(k,j,i-1) ) &
                                                  * ( pt(k,j,i-1) - pt(k,j,i) ) &
+                                               * rho_air_zw(k)                 &
+                                               * heatflux_output_conversion(k) &
                                                  * ddx * rmask(j,i,sr)
                    sums_l(k,61,tn) = sums_l(k,61,tn) - 0.5_wp *                &
                                                    ( kh(k,j,i) + kh(k,j-1,i) ) &
                                                  * ( pt(k,j-1,i) - pt(k,j,i) ) &
+                                               * rho_air_zw(k)                 &
+                                               * heatflux_output_conversion(k) &
                                                  * ddy * rmask(j,i,sr)
 !
@@ -3177,5 +3275,6 @@
                                      ( u(k,j,i) - hom(k,1,1,sr) ) * 0.5_wp *   &
                                      ( pt(k,j,i-1) - hom(k,1,4,sr) +           &
-                                       pt(k,j,i)   - hom(k,1,4,sr) )
+                                       pt(k,j,i)   - hom(k,1,4,sr) )           &
+                                     * heatflux_output_conversion(k)
                    pts = 0.5_wp * ( pt(k,j-1,i) - hom(k,1,4,sr) +              &
                                     pt(k,j,i)   - hom(k,1,4,sr) )
@@ -3183,5 +3282,6 @@
                                      ( v(k,j,i) - hom(k,1,2,sr) ) * 0.5_wp *   &
                                      ( pt(k,j-1,i) - hom(k,1,4,sr) +           &
-                                       pt(k,j,i)   - hom(k,1,4,sr) )
+                                       pt(k,j,i)   - hom(k,1,4,sr) )           &
+                                     * heatflux_output_conversion(k)
                 ENDDO
              ENDDO
@@ -3468,4 +3568,7 @@
        END IF
 
+       hom(:,1,121,sr) = rho_air       ! rho_air in Kg/m^3
+       hom(:,1,122,sr) = rho_air_zw    ! rho_air_zw in Kg/m^3
+
        hom(:,1,pr_palm,sr) =   sums(:,pr_palm)
                                        ! u*, w'u', w'v', t* (in last profile)
@@ -3576,5 +3679,6 @@
        IF ( hom(nzb,1,18,sr) > 1.0E-8_wp  .AND.  z_i(1) /= 0.0_wp )  THEN
           hom(nzb+8,1,pr_palm,sr) = &
-             ( g / hom(k_surface_level+1,1,4,sr) * hom(k_surface_level,1,18,sr)&
+             ( g / hom(k_surface_level+1,1,4,sr) *                             &
+             ( hom(k_surface_level,1,18,sr) / heatflux_output_conversion(nzb) )&
              * ABS( z_i(1) - mean_surface_level_height(sr) ) )**0.333333333_wp
        ELSE
Index: /palm/trunk/SOURCE/header.f90
===================================================================
--- /palm/trunk/SOURCE/header.f90	(revision 2036)
+++ /palm/trunk/SOURCE/header.f90	(revision 2037)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! Anelastic approximation implemented
 ! 
 ! Former revisions:
@@ -552,4 +552,5 @@
 !-- Numerical schemes
     WRITE ( io, 110 )
+    WRITE ( io, 121 )  TRIM( approximation )
     IF ( psolver(1:7) == 'poisfft' )  THEN
        WRITE ( io, 111 )  TRIM( fft_method )
@@ -1906,4 +1907,5 @@
 110 FORMAT (/' Numerical Schemes:'/ &
              ' -----------------'/)
+121 FORMAT (' --> Use the ',A,' approximation for the model equations.')
 111 FORMAT (' --> Solve perturbation pressure via FFT using ',A,' routines')
 112 FORMAT (' --> Solve perturbation pressure via SOR-Red/Black-Schema'/ &
Index: /palm/trunk/SOURCE/init_3d_model.f90
===================================================================
--- /palm/trunk/SOURCE/init_3d_model.f90	(revision 2036)
+++ /palm/trunk/SOURCE/init_3d_model.f90	(revision 2037)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! ------------------
-! 
+! Anelastic approximation implemented
 ! 
 ! Former revisions:
@@ -315,4 +315,7 @@
     USE arrays_3d
 
+    USE cloud_parameters,                                                      &
+        ONLY:  cp, l_v, r_d
+
     USE constants,                                                             &
         ONLY:  pi
@@ -324,5 +327,5 @@
     
     USE grid_variables,                                                        &
-        ONLY:  dx, dy
+        ONLY:  dx, dy, ddx2_mg, ddy2_mg
     
     USE indices
@@ -393,4 +396,13 @@
     INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  ngp_2dh_outer_l    !<
     INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  ngp_2dh_s_inner_l  !<
+
+    REAL(wp)     ::  t_surface !< air temperature at the surface
+
+    REAL(wp), DIMENSION(:), ALLOCATABLE ::  p_hydrostatic !< hydrostatic pressure
+
+    INTEGER(iwp) ::  l       !< loop variable
+    INTEGER(iwp) ::  nzt_l   !< index of top PE boundary for multigrid level
+    REAL(wp) ::  dx_l !< grid spacing along x on different multigrid level
+    REAL(wp) ::  dy_l !< grid spacing along y on different multigrid level
 
     REAL(wp), DIMENSION(1:3) ::  volume_flow_area_l     !<
@@ -647,4 +659,140 @@
 
 !
+!-- Allocation of anelastic and Boussinesq approximation specific arrays
+    ALLOCATE( p_hydrostatic(nzb:nzt+1) )
+    ALLOCATE( rho_air(nzb:nzt+1) )
+    ALLOCATE( rho_air_zw(nzb:nzt+1) )
+    ALLOCATE( drho_air(nzb:nzt+1) )
+    ALLOCATE( drho_air_zw(nzb:nzt+1) )
+
+!
+!-- Density profile calculation for anelastic approximation
+    IF ( TRIM( approximation ) == 'anelastic' ) THEN
+       t_surface = pt_surface * ( surface_pressure / 1000.0_wp )**( r_d / cp )
+       DO  k = nzb, nzt+1
+          p_hydrostatic(k)    = surface_pressure * 100.0_wp *                  &
+                                ( 1 - ( g * zu(k) ) / ( cp * t_surface )       &
+                                )**( cp / r_d )
+          rho_air(k)          = ( p_hydrostatic(k) *                           &
+                                  ( 100000.0_wp / p_hydrostatic(k)             &
+                                  )**( r_d / cp )                              &
+                                ) / ( r_d * pt_init(k) )
+       ENDDO
+       DO  k = nzb, nzt
+          rho_air_zw(k) = 0.5_wp * ( rho_air(k) + rho_air(k+1) )
+       ENDDO
+       rho_air_zw(nzt+1)  = rho_air_zw(nzt)                                    &
+                            + 2.0_wp * ( rho_air(nzt+1) - rho_air_zw(nzt)  )
+    ELSE
+       rho_air     = 1.0_wp
+       rho_air_zw  = 1.0_wp
+    ENDIF
+
+!-- compute the inverse density array in order to avoid expencive divisions
+    drho_air    = 1.0_wp / rho_air
+    drho_air_zw = 1.0_wp / rho_air_zw
+
+!
+!-- Allocation of flux conversion arrays
+    ALLOCATE( heatflux_input_conversion(nzb:nzt+1) )
+    ALLOCATE( waterflux_input_conversion(nzb:nzt+1) )
+    ALLOCATE( momentumflux_input_conversion(nzb:nzt+1) )
+    ALLOCATE( heatflux_output_conversion(nzb:nzt+1) )
+    ALLOCATE( waterflux_output_conversion(nzb:nzt+1) )
+    ALLOCATE( momentumflux_output_conversion(nzb:nzt+1) )
+
+!
+!-- calculate flux conversion factors according to approximation and in-/output mode
+    DO  k = nzb, nzt+1
+
+        IF ( TRIM( flux_input_mode ) == 'kinematic' )  THEN
+            heatflux_input_conversion(k)      = rho_air_zw(k)
+            waterflux_input_conversion(k)     = rho_air_zw(k)
+            momentumflux_input_conversion(k)  = rho_air_zw(k)
+        ELSEIF ( TRIM( flux_input_mode ) == 'dynamic' ) THEN
+            heatflux_input_conversion(k)      = 1.0_wp / cp
+            waterflux_input_conversion(k)     = 1.0_wp / l_v
+            momentumflux_input_conversion(k)  = 1.0_wp
+        ENDIF
+
+        IF ( TRIM( flux_output_mode ) == 'kinematic' )  THEN
+            heatflux_output_conversion(k)     = drho_air_zw(k)
+            waterflux_output_conversion(k)    = drho_air_zw(k)
+            momentumflux_output_conversion(k) = drho_air_zw(k)
+        ELSEIF ( TRIM( flux_output_mode ) == 'dynamic' ) THEN
+            heatflux_output_conversion(k)     = cp
+            waterflux_output_conversion(k)    = l_v
+            momentumflux_output_conversion(k) = 1.0_wp
+        ENDIF
+
+        IF ( .NOT. humidity ) THEN
+            waterflux_input_conversion(k)  = 1.0_wp
+            waterflux_output_conversion(k) = 1.0_wp
+        ENDIF
+
+    ENDDO
+
+!
+!-- In case of multigrid method, compute grid lengths and grid factors for the
+!-- grid levels with respective density on each grid
+    IF ( psolver(1:9) == 'multigrid' )  THEN
+
+       ALLOCATE( ddx2_mg(maximum_grid_level) )
+       ALLOCATE( ddy2_mg(maximum_grid_level) )
+       ALLOCATE( dzu_mg(nzb+1:nzt+1,maximum_grid_level) )
+       ALLOCATE( dzw_mg(nzb+1:nzt+1,maximum_grid_level) )
+       ALLOCATE( f1_mg(nzb+1:nzt,maximum_grid_level) )
+       ALLOCATE( f2_mg(nzb+1:nzt,maximum_grid_level) )
+       ALLOCATE( f3_mg(nzb+1:nzt,maximum_grid_level) )
+       ALLOCATE( rho_air_mg(nzb:nzt+1,maximum_grid_level) )
+       ALLOCATE( rho_air_zw_mg(nzb:nzt+1,maximum_grid_level) )
+
+       dzu_mg(:,maximum_grid_level) = dzu
+       rho_air_mg(:,maximum_grid_level) = rho_air
+!       
+!--    Next line to ensure an equally spaced grid. 
+       dzu_mg(1,maximum_grid_level) = dzu(2)
+       rho_air_mg(nzb,maximum_grid_level) = rho_air(nzb) +                     &
+                                             (rho_air(nzb) - rho_air(nzb+1))
+
+       dzw_mg(:,maximum_grid_level) = dzw
+       rho_air_zw_mg(:,maximum_grid_level) = rho_air_zw
+       nzt_l = nzt
+       DO  l = maximum_grid_level-1, 1, -1
+           dzu_mg(nzb+1,l) = 2.0_wp * dzu_mg(nzb+1,l+1)
+           dzw_mg(nzb+1,l) = 2.0_wp * dzw_mg(nzb+1,l+1)
+           rho_air_mg(nzb,l)    = rho_air_mg(nzb,l+1) + (rho_air_mg(nzb,l+1) - rho_air_mg(nzb+1,l+1))
+           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))
+           rho_air_mg(nzb+1,l)    = rho_air_mg(nzb+1,l+1)
+           rho_air_zw_mg(nzb+1,l) = rho_air_zw_mg(nzb+1,l+1)
+           nzt_l = nzt_l / 2
+           DO  k = 2, nzt_l+1
+              dzu_mg(k,l) = dzu_mg(2*k-2,l+1) + dzu_mg(2*k-1,l+1)
+              dzw_mg(k,l) = dzw_mg(2*k-2,l+1) + dzw_mg(2*k-1,l+1)
+              rho_air_mg(k,l)    = rho_air_mg(2*k-1,l+1)
+              rho_air_zw_mg(k,l) = rho_air_zw_mg(2*k-1,l+1)
+           ENDDO
+       ENDDO
+
+       nzt_l = nzt
+       dx_l  = dx
+       dy_l  = dy
+       DO  l = maximum_grid_level, 1, -1
+          ddx2_mg(l) = 1.0_wp / dx_l**2
+          ddy2_mg(l) = 1.0_wp / dy_l**2
+          DO  k = nzb+1, nzt_l
+             f2_mg(k,l) = rho_air_zw_mg(k,l) / ( dzu_mg(k+1,l) * dzw_mg(k,l) )
+             f3_mg(k,l) = rho_air_zw_mg(k-1,l) / ( dzu_mg(k,l)   * dzw_mg(k,l) )
+             f1_mg(k,l) = 2.0_wp * ( ddx2_mg(l) + ddy2_mg(l) ) &
+                          * rho_air_mg(k,l) + f2_mg(k,l) + f3_mg(k,l)
+          ENDDO
+          nzt_l = nzt_l / 2
+          dx_l  = dx_l * 2.0_wp
+          dy_l  = dy_l * 2.0_wp
+       ENDDO
+
+    ENDIF
+
+!
 !-- 3D-array for storing the dissipation, needed for calculating the sgs
 !-- particle velocities
@@ -883,6 +1031,6 @@
              vsws = 0.0_wp
           ENDIF
-          uswst = top_momentumflux_u
-          vswst = top_momentumflux_v
+          uswst = top_momentumflux_u * momentumflux_input_conversion(nzt+1)
+          vswst = top_momentumflux_v * momentumflux_input_conversion(nzt+1)
 
 !
@@ -1043,7 +1191,7 @@
           us    = 1E-30_wp 
           usws  = 0.0_wp
-          uswst = top_momentumflux_u
+          uswst = top_momentumflux_u * momentumflux_input_conversion(nzt+1)
           vsws  = 0.0_wp
-          vswst = top_momentumflux_v
+          vswst = top_momentumflux_v * momentumflux_input_conversion(nzt+1)
           IF ( humidity )       qs = 0.0_wp
           IF ( passive_scalar ) ss = 0.0_wp
@@ -1165,5 +1313,5 @@
                 CALL disturb_heatflux
              ELSE
-                shf = surface_heatflux
+                shf = surface_heatflux * heatflux_input_conversion(nzb)
 !
 !--             Initialize shf with data from external file LSF_DATA
@@ -1178,5 +1326,6 @@
                       DO  j = nysg, nyng
                          IF ( nzb_s_inner(j,i) /= 0 )  THEN
-                            shf(j,i) = wall_heatflux(0)
+                            shf(j,i) = wall_heatflux(0)                        &
+                                  * heatflux_input_conversion(nzb_s_inner(j,i))
                          ENDIF
                       ENDDO
@@ -1194,5 +1343,5 @@
              ENDIF
              IF ( constant_waterflux )  THEN
-                qsws   = surface_waterflux
+                qsws   = surface_waterflux * waterflux_input_conversion(nzb)
 !
 !--             Over topography surface_waterflux is replaced by 
@@ -1203,5 +1352,6 @@
                       DO  j = nysg, nyng
                          IF ( nzb_s_inner(j,i) /= 0 )  THEN
-                            qsws(j,i) = wall_qflux(0)
+                            qsws(j,i) = wall_qflux(0)                          &
+                                 * waterflux_input_conversion(nzb_s_inner(j,i))
                          ENDIF
                       ENDDO
@@ -1241,5 +1391,6 @@
 !
 !--       Prescribe to heat flux
-          IF ( constant_top_heatflux )  tswst = top_heatflux
+          IF ( constant_top_heatflux )  tswst = top_heatflux                   &
+                                             * heatflux_input_conversion(nzt+1)
 !
 !--       Prescribe zero latent flux at the top      
Index: /palm/trunk/SOURCE/init_grid.f90
===================================================================
--- /palm/trunk/SOURCE/init_grid.f90	(revision 2036)
+++ /palm/trunk/SOURCE/init_grid.f90	(revision 2037)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! Anelastic approximation implemented
 ! 
 ! Former revisions:
@@ -222,5 +222,5 @@
          
     USE grid_variables,                                                        &
-        ONLY:  ddx, ddx2, ddx2_mg, ddy, ddy2, ddy2_mg, dx, dx2, dy, dy2, fwxm, &
+        ONLY:  ddx, ddx2, ddy, ddy2, dx, dx2, dy, dy2, fwxm,                   &
                fwxp, fwym, fwyp, fxm, fxp, fym, fyp, wall_e_x, wall_e_y,       &
                wall_u, wall_v, wall_w_x, wall_w_y, zu_s_inner, zw_w_inner
@@ -293,6 +293,4 @@
 
     REAL(wp) ::  dum           !< dummy variable to skip columns while reading topography file    
-    REAL(wp) ::  dx_l          !< grid spacing along x on different multigrid level
-    REAL(wp) ::  dy_l          !< grid spacing along y on different multigrid level
     REAL(wp) ::  dz_stretched  !< stretched vertical grid spacing
 
@@ -431,52 +429,4 @@
        ddzu_pres = ddzu
        ddzu_pres(1) = ddzu_pres(2)  ! change for lowest level
-    ENDIF   
-
-!
-!-- In case of multigrid method, compute grid lengths and grid factors for the
-!-- grid levels
-    IF ( psolver(1:9) == 'multigrid' )  THEN
-
-       ALLOCATE( ddx2_mg(maximum_grid_level), ddy2_mg(maximum_grid_level), &
-                 dzu_mg(nzb+1:nzt+1,maximum_grid_level),                   &
-                 dzw_mg(nzb+1:nzt+1,maximum_grid_level),                   &
-                 f1_mg(nzb+1:nzt,maximum_grid_level),                      &
-                 f2_mg(nzb+1:nzt,maximum_grid_level),                      &
-                 f3_mg(nzb+1:nzt,maximum_grid_level) )
-
-       dzu_mg(:,maximum_grid_level) = dzu
-!       
-!--    Next line to ensure an equally spaced grid. 
-       dzu_mg(1,maximum_grid_level) = dzu(2)
-
-       dzw_mg(:,maximum_grid_level) = dzw
-       nzt_l = nzt
-       DO  l = maximum_grid_level-1, 1, -1
-           dzu_mg(nzb+1,l) = 2.0_wp * dzu_mg(nzb+1,l+1)
-           dzw_mg(nzb+1,l) = 2.0_wp * dzw_mg(nzb+1,l+1)
-           nzt_l = nzt_l / 2
-           DO  k = 2, nzt_l+1
-              dzu_mg(k,l) = dzu_mg(2*k-2,l+1) + dzu_mg(2*k-1,l+1)
-              dzw_mg(k,l) = dzw_mg(2*k-2,l+1) + dzw_mg(2*k-1,l+1)
-           ENDDO
-       ENDDO
-
-       nzt_l = nzt
-       dx_l  = dx
-       dy_l  = dy
-       DO  l = maximum_grid_level, 1, -1
-          ddx2_mg(l) = 1.0_wp / dx_l**2
-          ddy2_mg(l) = 1.0_wp / dy_l**2
-          DO  k = nzb+1, nzt_l
-             f2_mg(k,l) = 1.0_wp / ( dzu_mg(k+1,l) * dzw_mg(k,l) )
-             f3_mg(k,l) = 1.0_wp / ( dzu_mg(k,l)   * dzw_mg(k,l) )
-             f1_mg(k,l) = 2.0_wp * ( ddx2_mg(l) + ddy2_mg(l) ) + &
-                          f2_mg(k,l) + f3_mg(k,l)
-          ENDDO
-          nzt_l = nzt_l / 2
-          dx_l  = dx_l * 2.0_wp
-          dy_l  = dy_l * 2.0_wp
-       ENDDO
-
     ENDIF
 
Index: /palm/trunk/SOURCE/ls_forcing_mod.f90
===================================================================
--- /palm/trunk/SOURCE/ls_forcing_mod.f90	(revision 2036)
+++ /palm/trunk/SOURCE/ls_forcing_mod.f90	(revision 2037)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! ------------------
-! 
+! Anelastic approximation implemented
 ! 
 ! Former revisions:
@@ -109,4 +109,5 @@
            ONLY:  p_surf, pt_surf, q_surf, qsws_surf, shf_surf, td_lsa_lpt,    &
                   td_lsa_q, td_sub_lpt, td_sub_q, time_surf, time_vert,        &
+                  heatflux_input_conversion, waterflux_input_conversion,       &
                   ug_vert, vg_vert, wsubs_vert, zu
 
@@ -213,4 +214,6 @@
        ENDDO
 
+       shf_surf  = shf_surf  * heatflux_input_conversion(nzb)
+       qsws_surf = qsws_surf * waterflux_input_conversion(nzb)
 
        IF ( time_surf(1) > end_time )  THEN
Index: /palm/trunk/SOURCE/modules.f90
===================================================================
--- /palm/trunk/SOURCE/modules.f90	(revision 2036)
+++ /palm/trunk/SOURCE/modules.f90	(revision 2037)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! ------------------
-! 
+! Anelastic approximation implemented
 ! 
 ! Former revisions:
@@ -493,4 +493,19 @@
     REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE ::  rif_wall, tri
 
+    REAL(wp), DIMENSION(:), ALLOCATABLE :: rho_air     !< air density profile on the uv grid
+    REAL(wp), DIMENSION(:), ALLOCATABLE :: rho_air_zw  !< air density profile on the w grid
+    REAL(wp), DIMENSION(:), ALLOCATABLE :: drho_air    !< inverse air density profile on the uv grid
+    REAL(wp), DIMENSION(:), ALLOCATABLE :: drho_air_zw !< inverse air density profile on the w grid
+
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE :: rho_air_mg    !< air density profiles on the uv grid for multigrid
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE :: rho_air_zw_mg !< air density profiles on the w grid for multigrid
+
+    REAL(wp), DIMENSION(:), ALLOCATABLE :: heatflux_input_conversion      !< conversion factor array for heatflux input
+    REAL(wp), DIMENSION(:), ALLOCATABLE :: waterflux_input_conversion     !< conversion factor array for waterflux input
+    REAL(wp), DIMENSION(:), ALLOCATABLE :: momentumflux_input_conversion  !< conversion factor array for momentumflux input
+    REAL(wp), DIMENSION(:), ALLOCATABLE :: heatflux_output_conversion     !< conversion factor array for heatflux output
+    REAL(wp), DIMENSION(:), ALLOCATABLE :: waterflux_output_conversion    !< conversion factor array for waterflux output
+    REAL(wp), DIMENSION(:), ALLOCATABLE :: momentumflux_output_conversion !< conversion factor array for momentumflux output
+
 
     SAVE
@@ -616,4 +631,7 @@
                              psolver = 'poisfft', &
                              scalar_advec = 'ws-scheme'
+    CHARACTER (LEN=20)   ::  approximation = 'boussinesq'
+    CHARACTER (LEN=40)   ::  flux_input_mode = 'approximation-specific'
+    CHARACTER (LEN=40)   ::  flux_output_mode = 'approximation-specific'
     CHARACTER (LEN=20)   ::  bc_e_b = 'neumann', bc_lr = 'cyclic', &
                              bc_ns = 'cyclic', bc_p_b = 'neumann', &
Index: /palm/trunk/SOURCE/netcdf_interface_mod.f90
===================================================================
--- /palm/trunk/SOURCE/netcdf_interface_mod.f90	(revision 2036)
+++ /palm/trunk/SOURCE/netcdf_interface_mod.f90	(revision 2037)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! ------------------
-! 
+! Anelastic approximation implemented
 ! 
 ! Former revisions:
@@ -252,4 +252,8 @@
              ( 'unknown      ', i9 = 1, dots_max-42 ) /)
 
+    CHARACTER (LEN=16) :: heatflux_output_unit     !< unit for heatflux output
+    CHARACTER (LEN=16) :: waterflux_output_unit    !< unit for waterflux output
+    CHARACTER (LEN=16) :: momentumflux_output_unit !< unit for momentumflux output
+
     CHARACTER (LEN=9), DIMENSION(300) ::  dopr_unit = 'unknown'
 
@@ -368,7 +372,7 @@
             id_var_x_fl, id_var_y_fl, id_var_z_fl,  nc_stat,                   &
             netcdf_data_format, netcdf_data_format_string, netcdf_deflate,     &
-            netcdf_precision, output_for_t0
-            
-            
+            netcdf_precision, output_for_t0, heatflux_output_unit,             &
+            waterflux_output_unit, momentumflux_output_unit
+
     SAVE
 
Index: /palm/trunk/SOURCE/parin.f90
===================================================================
--- /palm/trunk/SOURCE/parin.f90	(revision 2036)
+++ /palm/trunk/SOURCE/parin.f90	(revision 2037)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! Anelastic approximation implemented
 ! 
 ! Former revisions:
@@ -313,5 +313,6 @@
 
 
-    NAMELIST /inipar/  alpha_surface, background_communication, bc_e_b, bc_lr, &
+    NAMELIST /inipar/  alpha_surface, approximation,                           &
+                       background_communication, bc_e_b, bc_lr,                &
                        bc_ns, bc_p_b, bc_p_t, bc_pt_b, bc_pt_t, bc_q_b,        &
              bc_q_t,bc_s_b, bc_s_t, bc_sa_t, bc_uv_b, bc_uv_t,                 &
@@ -332,6 +333,7 @@
              dt, dt_pr_1d, dt_run_control_1d, dx, dy, dz, dz_max,              &
              dz_stretch_factor, dz_stretch_level, end_time_1d,                 &
-             ensemble_member_nr,                                               &
-             e_init, e_min, fft_method, galilei_transformation, humidity,      &
+             ensemble_member_nr, e_init, e_min, fft_method,                    &
+             flux_input_mode, flux_output_mode,                                &
+             galilei_transformation, humidity,                                 &
              inflow_damping_height, inflow_damping_width,                      &
              inflow_disturbance_begin, inflow_disturbance_end,                 &
Index: /palm/trunk/SOURCE/poismg_mod.f90
===================================================================
--- /palm/trunk/SOURCE/poismg_mod.f90	(revision 2036)
+++ /palm/trunk/SOURCE/poismg_mod.f90	(revision 2037)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! Anelastic approximation implemented (stll error in optimized multigrid)
 ! 
 ! Former revisions:
@@ -262,5 +262,5 @@
 
        USE arrays_3d,                                                          &
-           ONLY:  f1_mg, f2_mg, f3_mg
+           ONLY:  f1_mg, f2_mg, f3_mg, rho_air_mg
 
        USE control_parameters,                                                 &
@@ -309,7 +309,7 @@
                 kp1 = k-ind_even_odd
                 r(k,j,i) = f_mg(k,j,i)                                         &
-                      - ddx2_mg(l) *                                           &
+                      - rho_air_mg(k,l) * ddx2_mg(l) *                         &
                       ( p_mg(k,j,i+1) +  p_mg(k,j,i-1)  )                      &
-                      - ddy2_mg(l) *                                           &
+                      - rho_air_mg(k,l) * ddy2_mg(l) *                         &
                       ( p_mg(k,j+1,i) + p_mg(k,j-1,i)  )                       &
                       - f2_mg_b(k,l) * p_mg(kp1,j,i)                           &
@@ -322,7 +322,7 @@
                 kp1 = k+ind_even_odd+1
                 r(k,j,i) = f_mg(k,j,i)                                         &
-                      - ddx2_mg(l) *                                           &
+                      - rho_air_mg(k,l) * ddx2_mg(l) *                         &
                       ( p_mg(k,j,i+1) +  p_mg(k,j,i-1)  )                      &
-                      - ddy2_mg(l) *                                           &
+                      - rho_air_mg(k,l) * ddy2_mg(l) *                         &
                       ( p_mg(k,j+1,i) + p_mg(k,j-1,i)  )                       &
                       - f2_mg_b(k,l) * p_mg(kp1,j,i)                           &
@@ -693,5 +693,5 @@
 
        USE arrays_3d,                                                          &
-           ONLY:  f1_mg, f2_mg, f3_mg
+           ONLY:  f1_mg, f2_mg, f3_mg, rho_air_mg
 
        USE control_parameters,                                                 &
@@ -755,7 +755,7 @@
                          kp1 = k-ind_even_odd
                          p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * (               &
-                                 ddx2_mg(l) *                                  &
+                                 rho_air_mg(k,l) * ddx2_mg(l) *                &
                                ( p_mg(k,j,i+1) + p_mg(k,j,i-1) )               &
-                               + ddy2_mg(l) *                                  &
+                               + rho_air_mg(k,l) * ddy2_mg(l) *                &
                                ( p_mg(k,j+1,i) + p_mg(k,j-1,i) )               &
                                + f2_mg_b(k,l) * p_mg(kp1,j,i)                  &
@@ -774,7 +774,7 @@
                          kp1 = k-ind_even_odd
                          p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * (               &
-                                 ddx2_mg(l) *                                  &
+                                 rho_air_mg(k,l) * ddx2_mg(l) *                &
                                ( p_mg(k,j,i+1) + p_mg(k,j,i-1) )               &
-                               + ddy2_mg(l) *                                  &
+                               + rho_air_mg(k,l) * ddy2_mg(l) *                &
                                ( p_mg(k,j+1,i) + p_mg(k,j-1,i) )               &
                                + f2_mg_b(k,l) * p_mg(kp1,j,i)                  &
@@ -793,7 +793,7 @@
                          kp1 = k+ind_even_odd+1
                          p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * (               &
-                                 ddx2_mg(l) *                                  &
+                                 rho_air_mg(k,l) * ddx2_mg(l) *                &
                                ( p_mg(k,j,i+1) + p_mg(k,j,i-1) )               &
-                               + ddy2_mg(l) *                                  &
+                               + rho_air_mg(k,l) * ddy2_mg(l) *                &
                                ( p_mg(k,j+1,i) + p_mg(k,j-1,i) )               &
                                + f2_mg_b(k,l) * p_mg(kp1,j,i)                  &
@@ -812,7 +812,7 @@
                          kp1 = k+ind_even_odd+1
                          p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * (               &
-                                 ddx2_mg(l) *                                  &
+                                 rho_air_mg(k,l) * ddx2_mg(l) *                &
                                ( p_mg(k,j,i+1) + p_mg(k,j,i-1) )               &
-                               + ddy2_mg(l) *                                  &
+                               + rho_air_mg(k,l) * ddy2_mg(l) *                &
                                ( p_mg(k,j+1,i) + p_mg(k,j-1,i) )               &
                                + f2_mg_b(k,l) * p_mg(kp1,j,i)                  &
@@ -843,7 +843,7 @@
                          j   = jj
                          p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * (               &
-                                 ddx2_mg(l) *                                  &
+                                 rho_air_mg(k,l) * ddx2_mg(l) *                &
                                ( p_mg(k,j,i+1) + p_mg(k,j,i-1) )               &
-                               + ddy2_mg(l) *                                  &
+                               + rho_air_mg(k,l) * ddy2_mg(l) *                &
                                ( p_mg(k,j+1,i) + p_mg(k,j-1,i) )               &
                                + f2_mg_b(k,l) * p_mg(kp1,j,i)                  &
@@ -852,7 +852,7 @@
                          j = jj+2
                          p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * (               &
-                                 ddx2_mg(l) *                                  &
+                                 rho_air_mg(k,l) * ddx2_mg(l) *                &
                                ( p_mg(k,j,i+1) + p_mg(k,j,i-1) )               &
-                               + ddy2_mg(l) *                                  &
+                               + rho_air_mg(k,l) * ddy2_mg(l) *                &
                                ( p_mg(k,j+1,i) + p_mg(k,j-1,i) )               &
                                + f2_mg_b(k,l) * p_mg(kp1,j,i)                  &
@@ -869,7 +869,7 @@
                          j   = jj
                          p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * (               &
-                                 ddx2_mg(l) *                                  &
+                                 rho_air_mg(k,l) * ddx2_mg(l) *                &
                                ( p_mg(k,j,i+1) + p_mg(k,j,i-1) )               &
-                               + ddy2_mg(l) *                                  &
+                               + rho_air_mg(k,l) * ddy2_mg(l) *                &
                                ( p_mg(k,j+1,i) + p_mg(k,j-1,i) )               &
                                + f2_mg_b(k,l) * p_mg(kp1,j,i)                  &
@@ -878,7 +878,7 @@
                          j = jj+2
                          p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * (               &
-                                 ddx2_mg(l) *                                  &
+                                 rho_air_mg(k,l) * ddx2_mg(l) *                &
                                ( p_mg(k,j,i+1) + p_mg(k,j,i-1) )               &
-                               + ddy2_mg(l) *                                  &
+                               + rho_air_mg(k,l) * ddy2_mg(l) *                &
                                ( p_mg(k,j+1,i) + p_mg(k,j-1,i) )               &
                                + f2_mg_b(k,l) * p_mg(kp1,j,i)                  &
@@ -895,7 +895,7 @@
                          j   = jj
                          p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * (               &
-                                 ddx2_mg(l) *                                  &
+                                 rho_air_mg(k,l) * ddx2_mg(l) *                &
                                ( p_mg(k,j,i+1) + p_mg(k,j,i-1) )               &
-                               + ddy2_mg(l) *                                  &
+                               + rho_air_mg(k,l) * ddy2_mg(l) *                &
                                ( p_mg(k,j+1,i) + p_mg(k,j-1,i) )               &
                                + f2_mg_b(k,l) * p_mg(kp1,j,i)                  &
@@ -904,7 +904,7 @@
                          j = jj+2
                          p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * (               &
-                                 ddx2_mg(l) *                                  &
+                                 rho_air_mg(k,l) * ddx2_mg(l) *                &
                                ( p_mg(k,j,i+1) + p_mg(k,j,i-1) )               &
-                               + ddy2_mg(l) *                                  &
+                               + rho_air_mg(k,l) * ddy2_mg(l) *                &
                                ( p_mg(k,j+1,i) + p_mg(k,j-1,i) )               &
                                + f2_mg_b(k,l) * p_mg(kp1,j,i)                  &
@@ -921,7 +921,7 @@
                          j   = jj
                          p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * (               &
-                                 ddx2_mg(l) *                                  &
+                                 rho_air_mg(k,l) * ddx2_mg(l) *                &
                                ( p_mg(k,j,i+1) + p_mg(k,j,i-1) )               &
-                               + ddy2_mg(l) *                                  &
+                               + rho_air_mg(k,l) * ddy2_mg(l) *                &
                                ( p_mg(k,j+1,i) + p_mg(k,j-1,i) )               &
                                + f2_mg_b(k,l) * p_mg(kp1,j,i)                  &
@@ -930,7 +930,7 @@
                          j = jj+2
                          p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * (               &
-                                 ddx2_mg(l) *                                  &
+                                 rho_air_mg(k,l) * ddx2_mg(l) *                &
                                ( p_mg(k,j,i+1) + p_mg(k,j,i-1) )               &
-                               + ddy2_mg(l) *                                  &
+                               + rho_air_mg(k,l) * ddy2_mg(l) *                &
                                ( p_mg(k,j+1,i) + p_mg(k,j-1,i) )               &
                                + f2_mg_b(k,l) * p_mg(kp1,j,i)                  &
Index: /palm/trunk/SOURCE/poismg_noopt.f90
===================================================================
--- /palm/trunk/SOURCE/poismg_noopt.f90	(revision 2036)
+++ /palm/trunk/SOURCE/poismg_noopt.f90	(revision 2037)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! Anelastic approximation implemented
 ! 
 ! Former revisions:
@@ -268,5 +268,5 @@
 
     USE arrays_3d,                                                             &
-        ONLY:  f1_mg, f2_mg, f3_mg
+        ONLY:  f1_mg, f2_mg, f3_mg, rho_air_mg
 
     USE control_parameters,                                                    &
@@ -339,10 +339,10 @@
           DO  k = nzb+1, nzt_mg(l)
              r(k,j,i) = f_mg(k,j,i)                                         &
-                        - ddx2_mg(l) *                                      &
+                        - rho_air_mg(k,l) * ddx2_mg(l) *                    &
                             ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
                                           ( p_mg(k,j,i) - p_mg(k,j,i+1) ) + &
                               p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * &
                                           ( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) &
-                        - ddy2_mg(l) *                                      &
+                        - rho_air_mg(k,l) * ddy2_mg(l) *                    &
                             ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
                                           ( p_mg(k,j,i) - p_mg(k,j+1,i) ) + &
@@ -745,5 +745,5 @@
 
     USE arrays_3d,                                                             &
-        ONLY:  f1_mg, f2_mg, f3_mg
+        ONLY:  f1_mg, f2_mg, f3_mg, rho_air_mg
 
     USE control_parameters,                                                    &
@@ -846,10 +846,10 @@
 
                       p_mg(k,j,i) = 1.0_wp / f1_mg(k,l) * (                    &
-                             ddx2_mg(l) *                                      &
+                             rho_air_mg(k,l) * ddx2_mg(l) *                    &
                                ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
                                              ( p_mg(k,j,i) - p_mg(k,j,i+1) ) + &
                                  p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * &
                                              ( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) &
-                           + ddy2_mg(l) *                                      &
+                           + rho_air_mg(k,l) * ddy2_mg(l) *                    &
                                ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
                                              ( p_mg(k,j,i) - p_mg(k,j+1,i) ) + &
@@ -869,10 +869,10 @@
                    DO  k = nzb+1, nzt_mg(l), 2 
                       p_mg(k,j,i) = 1.0_wp / f1_mg(k,l) * (                    &
-                             ddx2_mg(l) *                                      &
+                             rho_air_mg(k,l) * ddx2_mg(l) *                    &
                                ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
                                              ( p_mg(k,j,i) - p_mg(k,j,i+1) ) + &
                                  p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * &
                                              ( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) &
-                           + ddy2_mg(l) *                                      &
+                           + rho_air_mg(k,l) * ddy2_mg(l) *                    &
                                ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
                                              ( p_mg(k,j,i) - p_mg(k,j+1,i) ) + &
@@ -892,10 +892,10 @@
                    DO  k = nzb+2, nzt_mg(l), 2
                       p_mg(k,j,i) = 1.0_wp / f1_mg(k,l) * (                    &
-                             ddx2_mg(l) *                                      &
+                             rho_air_mg(k,l) * ddx2_mg(l) *                    &
                                ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
                                              ( p_mg(k,j,i) - p_mg(k,j,i+1) ) + &
                                  p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * &
                                              ( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) &
-                           + ddy2_mg(l) *                                      &
+                           + rho_air_mg(k,l) * ddy2_mg(l) *                    &
                                ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
                                              ( p_mg(k,j,i) - p_mg(k,j+1,i) ) + &
@@ -915,10 +915,10 @@
                    DO  k = nzb+2, nzt_mg(l), 2
                       p_mg(k,j,i) = 1.0_wp / f1_mg(k,l) * (                    &
-                             ddx2_mg(l) *                                      &
+                             rho_air_mg(k,l) * ddx2_mg(l) *                    &
                                ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
                                              ( p_mg(k,j,i) - p_mg(k,j,i+1) ) + &
                                  p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * &
                                              ( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) &
-                           + ddy2_mg(l) *                                      &
+                           + rho_air_mg(k,l) * ddy2_mg(l) *                    &
                                ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
                                              ( p_mg(k,j,i) - p_mg(k,j+1,i) ) + &
@@ -947,10 +947,10 @@
                       j = jj
                       p_mg(k,j,i) = 1.0_wp / f1_mg(k,l) * (                    &
-                             ddx2_mg(l) *                                      &
+                             rho_air_mg(k,l) * ddx2_mg(l) *                    &
                                ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
                                              ( p_mg(k,j,i) - p_mg(k,j,i+1) ) + &
                                  p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * &
                                              ( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) &
-                           + ddy2_mg(l) *                                      &
+                           + rho_air_mg(k,l) * ddy2_mg(l) *                    &
                                ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
                                              ( p_mg(k,j,i) - p_mg(k,j+1,i) ) + &
@@ -964,10 +964,10 @@
                       j = jj+2
                       p_mg(k,j,i) = 1.0_wp / f1_mg(k,l) * (                    &
-                             ddx2_mg(l) *                                      &
+                             rho_air_mg(k,l) * ddx2_mg(l) *                    &
                                ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
                                              ( p_mg(k,j,i) - p_mg(k,j,i+1) ) + &
                                  p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * &
                                              ( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) &
-                           + ddy2_mg(l) *                                      &
+                           + rho_air_mg(k,l) * ddy2_mg(l) *                    &
                                ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
                                              ( p_mg(k,j,i) - p_mg(k,j+1,i) ) + &
@@ -986,10 +986,10 @@
                       j =jj
                       p_mg(k,j,i) = 1.0_wp / f1_mg(k,l) * (                    &
-                             ddx2_mg(l) *                                      &
+                             rho_air_mg(k,l) * ddx2_mg(l) *                    &
                                ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
                                              ( p_mg(k,j,i) - p_mg(k,j,i+1) ) + &
                                  p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * &
                                              ( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) &
-                           + ddy2_mg(l) *                                      &
+                           + rho_air_mg(k,l) * ddy2_mg(l) *                    &
                                ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
                                              ( p_mg(k,j,i) - p_mg(k,j+1,i) ) + &
@@ -1003,10 +1003,10 @@
                       j = jj+2
                       p_mg(k,j,i) = 1.0_wp / f1_mg(k,l) * (                    &
-                             ddx2_mg(l) *                                      &
+                             rho_air_mg(k,l) * ddx2_mg(l) *                    &
                                ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
                                              ( p_mg(k,j,i) - p_mg(k,j,i+1) ) + &
                                  p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * &
                                              ( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) &
-                           + ddy2_mg(l) *                                      &
+                           + rho_air_mg(k,l) * ddy2_mg(l) *                    &
                                ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
                                              ( p_mg(k,j,i) - p_mg(k,j+1,i) ) + &
@@ -1025,10 +1025,10 @@
                       j =jj
                       p_mg(k,j,i) = 1.0_wp / f1_mg(k,l) * (                    &
-                             ddx2_mg(l) *                                      &
+                             rho_air_mg(k,l) * ddx2_mg(l) *                    &
                                ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
                                              ( p_mg(k,j,i) - p_mg(k,j,i+1) ) + &
                                  p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * &
                                              ( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) &
-                           + ddy2_mg(l) *                                      &
+                           + rho_air_mg(k,l) * ddy2_mg(l) *                    &
                                ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
                                              ( p_mg(k,j,i) - p_mg(k,j+1,i) ) + &
@@ -1042,10 +1042,10 @@
                       j = jj+2
                       p_mg(k,j,i) = 1.0_wp / f1_mg(k,l) * (                    &
-                             ddx2_mg(l) *                                      &
+                             rho_air_mg(k,l) * ddx2_mg(l) *                    &
                                ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
                                              ( p_mg(k,j,i) - p_mg(k,j,i+1) ) + &
                                  p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * &
                                              ( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) &
-                           + ddy2_mg(l) *                                      &
+                           + rho_air_mg(k,l) * ddy2_mg(l) *                    &
                                ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
                                              ( p_mg(k,j,i) - p_mg(k,j+1,i) ) + &
@@ -1064,10 +1064,10 @@
                       j =jj
                       p_mg(k,j,i) = 1.0_wp / f1_mg(k,l) * (                    &
-                             ddx2_mg(l) *                                      &
+                             rho_air_mg(k,l) * ddx2_mg(l) *                    &
                                ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
                                              ( p_mg(k,j,i) - p_mg(k,j,i+1) ) + &
                                  p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * &
                                              ( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) &
-                           + ddy2_mg(l) *                                      &
+                           + rho_air_mg(k,l) * ddy2_mg(l) *                    &
                                ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
                                              ( p_mg(k,j,i) - p_mg(k,j+1,i) ) + &
@@ -1081,10 +1081,10 @@
                       j = jj+2
                       p_mg(k,j,i) = 1.0_wp / f1_mg(k,l) * (                    &
-                             ddx2_mg(l) *                                      &
+                             rho_air_mg(k,l) * ddx2_mg(l) *                    &
                                ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
                                              ( p_mg(k,j,i) - p_mg(k,j,i+1) ) + &
                                  p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * &
                                              ( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) &
-                           + ddy2_mg(l) *                                      &
+                           + rho_air_mg(k,l) * ddy2_mg(l) *                    &
                                ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
                                              ( p_mg(k,j,i) - p_mg(k,j+1,i) ) + &
Index: /palm/trunk/SOURCE/pres.f90
===================================================================
--- /palm/trunk/SOURCE/pres.f90	(revision 2036)
+++ /palm/trunk/SOURCE/pres.f90	(revision 2037)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! ------------------
-! 
+! Anelastic approximation implemented
 ! 
 ! Former revisions:
@@ -128,5 +128,6 @@
 
     USE arrays_3d,                                                             &
-        ONLY:  d, ddzu, ddzu_pres, ddzw, dzw, p, p_loc, tend, u, v, w
+        ONLY:  d, ddzu, ddzu_pres, ddzw, dzw, p, p_loc, rho_air, rho_air_zw,   &
+               tend, u, v, w
 
     USE control_parameters,                                                    &
@@ -399,8 +400,9 @@
        DO  j = nys, nyn
           DO  k = nzb_s_inner(j,i)+1, nzt
-             d(k,j,i) = ( ( u(k,j,i+1) - u(k,j,i) ) * ddx + &
-                          ( v(k,j+1,i) - v(k,j,i) ) * ddy + &
-                          ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) ) * ddt_3d      &
-                                                                * d_weight_pres 
+             d(k,j,i) = ( ( u(k,j,i+1) - u(k,j,i) ) * rho_air(k) * ddx +       &
+                          ( v(k,j+1,i) - v(k,j,i) ) * rho_air(k) * ddy +       &
+                          ( w(k,j,i)   * rho_air_zw(k) -                       &
+                            w(k-1,j,i) * rho_air_zw(k-1) ) * ddzw(k)           &
+                        ) * ddt_3d * d_weight_pres 
           ENDDO
 !
@@ -427,8 +429,9 @@
        DO  j = nys, nyn
           DO  k = 1, nzt
-             d(k,j,i) = ( ( u(k,j,i+1) - u(k,j,i) ) * ddx +             &
-                        ( v(k,j+1,i) - v(k,j,i) ) * ddy +               &
-                        ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) ) * ddt_3d  &
-                        * d_weight_pres * rflags_s_inner(k,j,i)
+             d(k,j,i) = ( ( u(k,j,i+1) - u(k,j,i) ) * rho_air(k) * ddx +       &
+                          ( v(k,j+1,i) - v(k,j,i) ) * rho_air(k) * ddy +       &
+                          ( w(k,j,i)   * rho_air_zw(k) -                       &
+                            w(k-1,j,i) * rho_air_zw(k-1) ) * ddzw(k)           &
+                        ) * ddt_3d * d_weight_pres * rflags_s_inner(k,j,i)
           ENDDO
        ENDDO
@@ -808,7 +811,8 @@
           DO  j = nys, nyn
              DO  k = nzb_s_inner(j,i)+1, nzt
-                d(k,j,i) = ( u(k,j,i+1) - u(k,j,i) ) * ddx +                   &
-                           ( v(k,j+1,i) - v(k,j,i) ) * ddy +                   &
-                           ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
+             d(k,j,i) = ( u(k,j,i+1) - u(k,j,i) ) * rho_air(k) * ddx +         &
+                        ( v(k,j+1,i) - v(k,j,i) ) * rho_air(k) * ddy +         &
+                        ( w(k,j,i)   * rho_air_zw(k) -                         &
+                          w(k-1,j,i) * rho_air_zw(k-1) ) * ddzw(k)
              ENDDO
              DO  k = nzb+1, nzt
@@ -823,8 +827,9 @@
           DO  j = nys, nyn
              DO  k = 1, nzt
-                   d(k,j,i) = ( ( u(k,j,i+1) - u(k,j,i) ) * ddx +              &
-                                ( v(k,j+1,i) - v(k,j,i) ) * ddy +              &
-                                ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)            &
-                              ) * rflags_s_inner(k,j,i)
+                d(k,j,i) = ( ( u(k,j,i+1) - u(k,j,i) ) * rho_air(k) * ddx +    &
+                             ( v(k,j+1,i) - v(k,j,i) ) * rho_air(k) * ddy +    &
+                             ( w(k,j,i)   * rho_air_zw(k) -                    &
+                               w(k-1,j,i) * rho_air_zw(k-1) ) * ddzw(k)        &
+                           ) * rflags_s_inner(k,j,i)
              ENDDO
           ENDDO
Index: /palm/trunk/SOURCE/sor.f90
===================================================================
--- /palm/trunk/SOURCE/sor.f90	(revision 2036)
+++ /palm/trunk/SOURCE/sor.f90	(revision 2037)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! Anelastic approximation implemented
 ! 
 ! Former revisions:
@@ -59,5 +59,7 @@
 !------------------------------------------------------------------------------!
  SUBROUTINE sor( d, ddzu, ddzw, p )
- 
+
+    USE arrays_3d,                                                             &
+        ONLY:  rho_air, rho_air_zw
 
     USE grid_variables,                                                        &
@@ -101,7 +103,7 @@
 !-- Compute pre-factors.
     DO  k = 1, nz
-         f2(k) = ddzu(k+1) * ddzw(k)
-         f3(k) = ddzu(k)   * ddzw(k)
-         f1(k) = 2.0_wp * ( ddx2 + ddy2 ) + f2(k) + f3(k)
+         f2(k) = ddzu(k+1) * ddzw(k) * rho_air_zw(k)
+         f3(k) = ddzu(k)   * ddzw(k) * rho_air_zw(k-1)
+         f1(k) = 2.0_wp * ( ddx2 + ddy2 ) * rho_air(k) + f2(k) + f3(k)
     ENDDO
 
@@ -131,10 +133,10 @@
              DO  k = nzb+1, nzt
                 p(k,j,i) = p(k,j,i) + omega_sor / f1(k) * (            &
-                               ddx2 * ( p(k,j,i+1) + p(k,j,i-1) ) +    &
-                               ddy2 * ( p(k,j+1,i) + p(k,j-1,i) ) +    &
-                               f2(k) * p(k+1,j,i)                 +    &
-                               f3(k) * p(k-1,j,i)                 -    &
-                               d(k,j,i)                           -    &
-                               f1(k) * p(k,j,i)           )
+                           rho_air(k) * ddx2 * ( p(k,j,i+1) + p(k,j,i-1) ) +   &
+                           rho_air(k) * ddy2 * ( p(k,j+1,i) + p(k,j-1,i) ) +   &
+                           f2(k) * p(k+1,j,i)                              +   &
+                           f3(k) * p(k-1,j,i)                              -   &
+                           d(k,j,i)                                        -   &
+                           f1(k) * p(k,j,i)           )
              ENDDO
           ENDDO
@@ -144,11 +146,11 @@
           DO  j = nys1, nyn, 2
              DO  k = nzb+1, nzt
-                p(k,j,i) = p(k,j,i) + omega_sor / f1(k) * (            &
-                               ddx2 * ( p(k,j,i+1) + p(k,j,i-1) ) +    &
-                               ddy2 * ( p(k,j+1,i) + p(k,j-1,i) ) +    &
-                               f2(k) * p(k+1,j,i)                 +    &
-                               f3(k) * p(k-1,j,i)                 -    &
-                               d(k,j,i)                           -    &
-                               f1(k) * p(k,j,i)           )
+                p(k,j,i) = p(k,j,i) + omega_sor / f1(k) * (                    &
+                           rho_air(k) * ddx2 * ( p(k,j,i+1) + p(k,j,i-1) ) +   &
+                           rho_air(k) * ddy2 * ( p(k,j+1,i) + p(k,j-1,i) ) +   &
+                           f2(k) * p(k+1,j,i)                              +   &
+                           f3(k) * p(k-1,j,i)                              -   &
+                           d(k,j,i)                                        -   &
+                           f1(k) * p(k,j,i)           )
              ENDDO
           ENDDO
@@ -176,10 +178,10 @@
              DO  k = nzb+1, nzt
                 p(k,j,i) = p(k,j,i) + omega_sor / f1(k) * (            &
-                               ddx2 * ( p(k,j,i+1) + p(k,j,i-1) ) +    &
-                               ddy2 * ( p(k,j+1,i) + p(k,j-1,i) ) +    &
-                               f2(k) * p(k+1,j,i)                 +    &
-                               f3(k) * p(k-1,j,i)                 -    &
-                               d(k,j,i)                           -    &
-                               f1(k) * p(k,j,i)           )
+                           rho_air(k) * ddx2 * ( p(k,j,i+1) + p(k,j,i-1) ) +   &
+                           rho_air(k) * ddy2 * ( p(k,j+1,i) + p(k,j-1,i) ) +   &
+                           f2(k) * p(k+1,j,i)                              +   &
+                           f3(k) * p(k-1,j,i)                              -   &
+                           d(k,j,i)                                        -   &
+                           f1(k) * p(k,j,i)           )
              ENDDO
           ENDDO
@@ -190,10 +192,10 @@
              DO  k = nzb+1, nzt
                 p(k,j,i) = p(k,j,i) + omega_sor / f1(k) * (            &
-                               ddx2 * ( p(k,j,i+1) + p(k,j,i-1) ) +    &
-                               ddy2 * ( p(k,j+1,i) + p(k,j-1,i) ) +    &
-                               f2(k) * p(k+1,j,i)                 +    &
-                               f3(k) * p(k-1,j,i)                 -    &
-                               d(k,j,i)                           -    &
-                               f1(k) * p(k,j,i)           )
+                           rho_air(k) * ddx2 * ( p(k,j,i+1) + p(k,j,i-1) ) +   &
+                           rho_air(k) * ddy2 * ( p(k,j+1,i) + p(k,j-1,i) ) +   &
+                           f2(k) * p(k+1,j,i)                              +   &
+                           f3(k) * p(k-1,j,i)                              -   &
+                           d(k,j,i)                                        -   &
+                           f1(k) * p(k,j,i)           )
              ENDDO
           ENDDO
Index: /palm/trunk/SOURCE/surface_layer_fluxes_mod.f90
===================================================================
--- /palm/trunk/SOURCE/surface_layer_fluxes_mod.f90	(revision 2036)
+++ /palm/trunk/SOURCE/surface_layer_fluxes_mod.f90	(revision 2037)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! ------------------
-! 
+! Anelastic approximation implemented
 ! 
 ! Former revisions:
@@ -158,5 +158,5 @@
         ONLY:  e, kh, nr, nrs, nrsws, ol, pt, q, ql, qr, qrs, qrsws, qs, qsws, &
                s, shf, ss, ssws, ts, u, us, usws, v, vpt, vsws, zu, zw, z0,    &
-               z0h, z0q
+               z0h, z0q, drho_air_zw, rho_air_zw
 
     USE cloud_parameters,                                                      &
@@ -540,9 +540,9 @@
                       rib(j,i) = - g * z_mo * ( ( 1.0_wp + 0.61_wp             &
                                  * q(k+1,j,i) ) * shf(j,i) + 0.61_wp           &
-                                 * pt(k+1,j,i) * qsws(j,i) )   &
+                                 * pt(k+1,j,i) * qsws(j,i) ) * drho_air_zw(k)  &
                                  / ( uv_total(j,i)**3 * vpt(k+1,j,i) * kappa**2&
                                  + 1.0E-20_wp)
                    ELSE
-                      rib(j,i) = - g * z_mo * shf(j,i)                         &
+                      rib(j,i) = - g * z_mo * shf(j,i) * drho_air_zw(k)        &
                            / ( uv_total(j,i)**3 * pt(k+1,j,i) * kappa**2       &
                            + 1.0E-20_wp )
@@ -828,8 +828,9 @@
 !--       For a given heat flux in the surface layer:
           !$OMP PARALLEL DO
-          !$acc kernels loop private( j )
+          !$acc kernels loop private( j, k )
           DO  i = nxlg, nxrg
              DO  j = nysg, nyng
-                ts(j,i) = -shf(j,i) / ( us(j,i) + 1E-30_wp )
+                k   = nzb_s_inner(j,i)
+                ts(j,i) = -shf(j,i) * drho_air_zw(k) / ( us(j,i) + 1E-30_wp )
 !
 !--             ts must be limited, because otherwise overflow may occur in case
@@ -888,5 +889,6 @@
              DO  i = nxlg, nxrg
                 DO  j = nysg, nyng
-                   qs(j,i) = -qsws(j,i) / ( us(j,i) + 1E-30_wp )
+                   k   = nzb_s_inner(j,i)
+                   qs(j,i) = -qsws(j,i) * drho_air_zw(k) / ( us(j,i) + 1E-30_wp )
                 ENDDO
              ENDDO
@@ -1025,5 +1027,6 @@
                                      + psi_m( z0(j,i) / ol_mid ) )
 
-             usws(j,i) = -usws(j,i) * 0.5_wp * ( us(j,i-1) + us(j,i) )
+             usws(j,i) = -usws(j,i) * 0.5_wp * ( us(j,i-1) + us(j,i) )         &
+                                    * rho_air_zw(k)
           ENDDO
        ENDDO
@@ -1052,5 +1055,6 @@
                                      + psi_m( z0(j,i) / ol_mid ) )
 
-             vsws(j,i) = -vsws(j,i) * 0.5_wp * ( us(j,i-1) + us(j,i) )
+             vsws(j,i) = -vsws(j,i) * 0.5_wp * ( us(j,i-1) + us(j,i) )         &
+                                    * rho_air_zw(k)
 
           ENDDO
@@ -1074,5 +1078,6 @@
              !$acc loop independent
              DO  j = nysg, nyng
-                shf(j,i) = -ts(j,i) * us(j,i)
+                k   = nzb_s_inner(j,i)
+                shf(j,i) = -ts(j,i) * us(j,i) * rho_air_zw(k)
              ENDDO
           ENDDO
@@ -1090,5 +1095,6 @@
              !$acc loop independent
              DO  j = nysg, nyng
-                qsws(j,i) = -qs(j,i) * us(j,i)
+                k   = nzb_s_inner(j,i)
+                qsws(j,i) = -qs(j,i) * us(j,i) * rho_air_zw(k)
              ENDDO
           ENDDO
Index: /palm/trunk/SOURCE/tridia_solver_mod.f90
===================================================================
--- /palm/trunk/SOURCE/tridia_solver_mod.f90	(revision 2036)
+++ /palm/trunk/SOURCE/tridia_solver_mod.f90	(revision 2037)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! ------------------
-! 
+! Anelastic approximation implemented
 ! 
 ! Former revisions:
@@ -129,5 +129,5 @@
 
        USE arrays_3d,                                                          &
-           ONLY:  ddzu_pres, ddzw
+           ONLY:  ddzu_pres, ddzw, rho_air_zw
 
        USE kinds
@@ -140,8 +140,9 @@
 
        DO  k = 0, nz-1
-          ddzuw(k,1) = ddzu_pres(k+1) * ddzw(k+1)
-          ddzuw(k,2) = ddzu_pres(k+2) * ddzw(k+1)
+          ddzuw(k,1) = ddzu_pres(k+1) * ddzw(k+1) * rho_air_zw(k)
+          ddzuw(k,2) = ddzu_pres(k+2) * ddzw(k+1) * rho_air_zw(k+1)
           ddzuw(k,3) = -1.0_wp * &
-                       ( ddzu_pres(k+2) * ddzw(k+1) + ddzu_pres(k+1) * ddzw(k+1) )
+                       ( ddzu_pres(k+2) * ddzw(k+1) * rho_air_zw(k+1) +        &
+                         ddzu_pres(k+1) * ddzw(k+1) * rho_air_zw(k) )
        ENDDO
 !
@@ -168,5 +169,5 @@
 
           USE arrays_3d,                                                       &
-              ONLY:  tric
+              ONLY:  tric, rho_air
 
           USE constants,                                                       &
@@ -231,5 +232,5 @@
              DO  j = nys_z, nyn_z
                 DO  i = nxl_z, nxr_z
-                   tric(i,j,k) = ddzuw(k,3) - ll(i,j)
+                   tric(i,j,k) = ddzuw(k,3) - ll(i,j) * rho_air(k+1)
                 ENDDO
              ENDDO
@@ -491,5 +492,5 @@
 
        USE arrays_3d,                                                          &
-           ONLY:  ddzu_pres, ddzw
+           ONLY:  ddzu_pres, ddzw, rho_air, rho_air_zw
 
        USE control_parameters,                                                 &
@@ -526,6 +527,6 @@
        DO  k = 0, nz-1
           DO  i = 0,nx
-             tri_for_1d(2,i,k) = ddzu_pres(k+1) * ddzw(k+1)
-             tri_for_1d(3,i,k) = ddzu_pres(k+2) * ddzw(k+1)
+             tri_for_1d(2,i,k) = ddzu_pres(k+1) * ddzw(k+1) * rho_air_zw(k)
+             tri_for_1d(3,i,k) = ddzu_pres(k+2) * ddzw(k+1) * rho_air_zw(k+1)
           ENDDO
        ENDDO
@@ -591,7 +592,7 @@
           DO  k = 0, nz-1
              DO  i = 0, nx
-                a = -1.0_wp * ddzu_pres(k+2) * ddzw(k+1)
-                c = -1.0_wp * ddzu_pres(k+1) * ddzw(k+1)
-                tri_for_1d(1,i,k) = a + c - l(i)
+                a = -1.0_wp * ddzu_pres(k+2) * ddzw(k+1) * rho_air_zw(k+1)
+                c = -1.0_wp * ddzu_pres(k+1) * ddzw(k+1) * rho_air_zw(k)
+                tri_for_1d(1,i,k) = a + c - l(i) * rho_air(k+1)
              ENDDO
           ENDDO
