Index: palm/trunk/SOURCE/advec_ws.f90
===================================================================
--- palm/trunk/SOURCE/advec_ws.f90	(revision 2292)
+++ palm/trunk/SOURCE/advec_ws.f90	(revision 2329)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$ 
+! Bugfix concerning density in divergence correction close to buildings
+! 
+! 2292 2017-06-20 09:51:42Z schwenkel
 ! Implementation of new microphysic scheme: cloud_scheme = 'morrison' 
 ! includes two more prognostic equations for cloud drop concentration (nc)  
@@ -1094,5 +1097,5 @@
 
        USE arrays_3d,                                                          &
-           ONLY:  ddzw, drho_air, tend, u, v, w, rho_air, rho_air_zw
+           ONLY:  ddzw, drho_air, tend, u, v, w, rho_air_zw
 
        USE constants,                                                          &
@@ -1422,5 +1425,5 @@
                                          + IBITS(advc_flags_1(k,j,i-1),2,1)    &
                                          )                                     &
-                          ) * rho_air(k) * ddx                                 &
+                          ) * ddx                                              &
                         + ( v(k,j+1,i) * ( ibit3 + ibit4 + ibit5 )             &
                           - v(k,j,i)   * ( IBITS(advc_flags_1(k,j-1,i),3,1)    &
@@ -1428,5 +1431,5 @@
                                          + IBITS(advc_flags_1(k,j-1,i),5,1)    &
                                          )                                     &
-                          ) * rho_air(k) * ddy                                 &
+                          ) * ddy                                              &
                         + ( w(k,j,i) * rho_air_zw(k) *                         &
                                          ( ibit6 + ibit7 + ibit8 )             &
@@ -1436,5 +1439,5 @@
                                          + IBITS(advc_flags_1(k-1,j,i),8,1)    &
                                          )                                     &      
-                          ) * ddzw(k)
+                          ) * drho_air(k) * ddzw(k)
 
 
@@ -1525,8 +1528,9 @@
 !--       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)   ) * rho_air(k) * ddx      &
-                        + ( v(k,j+1,i) - v(k,j,i)   ) * rho_air(k) * ddy      &
+          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)   * rho_air_zw(k) -                      &
-                            w(k-1,j,i) * rho_air_zw(k-1) ) * ddzw(k)
+                            w(k-1,j,i) * rho_air_zw(k-1)                      &
+                          ) * drho_air(k) * ddzw(k)
 
           tend(k,j,i) = tend(k,j,i) - (                                       &
@@ -1659,5 +1663,5 @@
        USE arrays_3d,                                                         &
            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
+                  drho_air, rho_air_zw
 
        USE constants,                                                         &
@@ -1963,5 +1967,5 @@
                                       + IBITS(advc_flags_1(k,j,i-1),11,1)     &
                                       )                                       &
-                  ) * rho_air(k) * ddx                                        &
+                  ) * ddx                                                     &
                +  ( ( v_comp + gv ) * ( ibit12 + ibit13 + ibit14 )            &
                   - ( v(k,j,i)   + v(k,j,i-1 )  )                             &
@@ -1970,5 +1974,5 @@
                                       + IBITS(advc_flags_1(k,j-1,i),14,1)     &
                                       )                                       &
-                  ) * rho_air(k) * ddy                                        &
+                  ) * 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)             &
@@ -1977,5 +1981,5 @@
                                       + IBITS(advc_flags_1(k-1,j,i),17,1)     &
                                       )                                       &  
-                  ) * ddzw(k)   &
+                  ) * drho_air(k) * ddzw(k)                                   &
                 ) * 0.5_wp
 
@@ -2087,10 +2091,8 @@
 !--       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       &
-                                                                  * 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)                                                 &
+                  ) * drho_air(k) * ddzw(k)                                   &
                 ) * 0.5_wp
 
@@ -2151,5 +2153,5 @@
        USE arrays_3d,                                                          &
            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
+                  drho_air, rho_air_zw
 
        USE constants,                                                          &
@@ -2458,5 +2460,5 @@
                                          + IBITS(advc_flags_1(k,j,i-1),20,1)  &
                                          )                                    &
-                  ) * rho_air(k) * ddx                                        &
+                  ) * ddx                                                     &
                +  ( v_comp(k)                                                 &
                                        * ( ibit21 + ibit22 + ibit23 )         &
@@ -2466,5 +2468,5 @@
                                          + IBITS(advc_flags_1(k,j-1,i),23,1)  &
                                          )                                    &
-                  ) * rho_air(k) * ddy                                        &
+                  ) * 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)             &
@@ -2473,5 +2475,5 @@
                                          + IBITS(advc_flags_1(k-1,j,i),26,1)  &
                                          )                                    &
-                   ) * ddzw(k)   &
+                   ) * drho_air(k) * ddzw(k)                                  &
                 ) * 0.5_wp
 
@@ -2588,10 +2590,8 @@
 !--       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       &
-                                                                  * 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)                                                 &
+                  ) * drho_air(k) * ddzw(k)                                   &
                 ) * 0.5_wp
 
@@ -2652,5 +2652,5 @@
        USE arrays_3d,                                                         &
            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
+                  drho_air_zw, rho_air
 
        USE constants,                                                         &
@@ -2962,5 +2962,5 @@
                                       + IBITS(advc_flags_1(k,j,i-1),29,1)     &
                                       )                                       &
-                  ) * rho_air_zw(k) * ddx                                     &
+                  ) * ddx                                                     &
               +   ( ( v_comp + gv ) * ( ibit30 + ibit31 + ibit32 )            & 
                   - ( v(k+1,j,i) + v(k,j,i)   )                               &
@@ -2969,5 +2969,5 @@
                                       + IBITS(advc_flags_2(k,j-1,i),0,1)      &
                                       )                                       &
-                  ) * rho_air_zw(k) * ddy                                     &
+                  ) * ddy                                                     &
               +   ( w_comp * rho_air(k+1) * ( ibit33 + ibit34 + ibit35 )      &
                 - ( w(k,j,i)   + w(k-1,j,i)   ) * rho_air(k)                  &
@@ -2976,5 +2976,5 @@
                                       + IBITS(advc_flags_2(k-1,j,i),3,1)      &
                                       )                                       & 
-                  ) * ddzu(k+1)   &
+                  ) * drho_air_zw(k) * ddzu(k+1)                              &
                 ) * 0.5_wp
 
@@ -3078,10 +3078,8 @@
 !--       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         &
-                                                              * 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)                                               &
+                  ) * drho_air_zw(k) * ddzu(k+1)                              &
                 ) * 0.5_wp
 
@@ -3127,5 +3125,5 @@
 
        USE arrays_3d,                                                         &
-           ONLY:  ddzw, drho_air, tend, u, v, w, rho_air, rho_air_zw
+           ONLY:  ddzw, drho_air, tend, u, v, w, rho_air_zw
 
        USE constants,                                                         &
@@ -3448,5 +3446,5 @@
                                          + IBITS(advc_flags_1(k,j,i-1),2,1)    &
                                          )                                     &
-                          ) * rho_air(k) * ddx                                 &
+                          ) * ddx                                              &
                         + ( v(k,j+1,i) * ( ibit3 + ibit4 + ibit5 )             &
                           - v(k,j,i)   * ( IBITS(advc_flags_1(k,j-1,i),3,1)    &
@@ -3454,5 +3452,5 @@
                                          + IBITS(advc_flags_1(k,j-1,i),5,1)    &
                                          )                                     &
-                          ) * rho_air(k) * ddy                                 &
+                          ) * ddy                                              &
                         + ( w(k,j,i) * rho_air_zw(k) *                         &
                                          ( ibit6 + ibit7 + ibit8 )             &
@@ -3462,5 +3460,5 @@
                                          + IBITS(advc_flags_1(k-1,j,i),8,1)    &
                                          )                                     &      
-                          ) * ddzw(k)
+                          ) * drho_air(k) * ddzw(k)
 
 
@@ -3549,8 +3547,9 @@
 !--             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)   ) * rho_air(k) * ddx &
-                              + ( v(k,j+1,i) - v(k,j,i)   ) * rho_air(k) * ddy &
+                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)   * rho_air_zw(k) -                 &
-                                  w(k-1,j,i) * rho_air_zw(k-1) ) * ddzw(k)
+                                  w(k-1,j,i) * rho_air_zw(k-1)                 &
+                                ) * drho_air(k) * ddzw(k)
 
                 tend(k,j,i) = tend(k,j,i) - (                                 &
@@ -3682,5 +3681,5 @@
 
        USE arrays_3d,                                                          &
-           ONLY:  ddzw, drho_air, tend, u, v, w, rho_air, rho_air_zw
+           ONLY:  ddzw, drho_air, tend, u, v, w, rho_air_zw
 
        USE constants,                                                          &
@@ -3990,5 +3989,5 @@
                                       + IBITS(advc_flags_1(k,j,i-1),11,1)     &
                                       )                                       &
-                  ) * rho_air(k) * ddx                                        &
+                  ) * ddx                                                     &
                +  ( ( v_comp + gv ) * ( ibit12 + ibit13 + ibit14 )            &
                   - ( v(k,j,i)   + v(k,j,i-1 )  )                             &
@@ -3997,5 +3996,5 @@
                                       + IBITS(advc_flags_1(k,j-1,i),14,1)     &
                                       )                                       &
-                  ) * rho_air(k) * ddy                                        &
+                  ) * 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)             &
@@ -4004,5 +4003,5 @@
                                       + IBITS(advc_flags_1(k-1,j,i),17,1)     &
                                       )                                       &  
-                  ) * ddzw(k)   &
+                  ) * drho_air(k) * ddzw(k)                                   &
                 ) * 0.5_wp
 
@@ -4116,10 +4115,8 @@
 !--             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 &
-                                                                  * 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)                                           &
+                        ) * drho_air(k) * ddzw(k)                             &
                       ) * 0.5_wp
 
@@ -4178,5 +4175,5 @@
 
        USE arrays_3d,                                                          &
-           ONLY:  ddzw, drho_air, tend, u, v, w, rho_air, rho_air_zw
+           ONLY:  ddzw, drho_air, tend, u, v, w, rho_air_zw
 
        USE constants,                                                          &
@@ -4486,5 +4483,5 @@
                                          + IBITS(advc_flags_1(k,j,i-1),20,1)  &
                                          )                                    &
-                  ) * rho_air(k) * ddx                                        &
+                  ) * ddx                                                     &
                +  ( v_comp(k)                                                 &
                                        * ( ibit21 + ibit22 + ibit23 )         &
@@ -4494,5 +4491,5 @@
                                          + IBITS(advc_flags_1(k,j-1,i),23,1)  &
                                          )                                    &
-                  ) * rho_air(k) * ddy                                        &
+                  ) * ddy                                                     &
                +  ( w_comp * rho_air_zw(k)                                    &
                                        * ( ibit24 + ibit25 + ibit26 )         &
@@ -4502,5 +4499,5 @@
                                          + IBITS(advc_flags_1(k-1,j,i),26,1)  &
                                          )                                    &
-                   ) * ddzw(k)   &
+                   ) * drho_air(k) * ddzw(k)                                  &
                 ) * 0.5_wp
 
@@ -4619,10 +4616,8 @@
 !--             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 &
-                                                                  * 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)                                           &
+                        ) * drho_air(k) * ddzw(k)                             &
                       ) * 0.5_wp
  
@@ -4685,5 +4680,5 @@
 
        USE arrays_3d,                                                          &
-           ONLY:  ddzu, drho_air_zw, tend, u, v, w, rho_air, rho_air_zw
+           ONLY:  ddzu, drho_air_zw, tend, u, v, w, rho_air
 
        USE constants,                                                          &
@@ -4998,5 +4993,5 @@
                                       + IBITS(advc_flags_1(k,j,i-1),29,1)     &
                                       )                                       &
-                  ) * rho_air_zw(k) * ddx                                     &
+                  ) * ddx                                                     &
               +   ( ( v_comp + gv ) * ( ibit30 + ibit31 + ibit32 )            &
                   - ( v(k+1,j,i) + v(k,j,i)   )                               &
@@ -5005,5 +5000,5 @@
                                       + IBITS(advc_flags_2(k,j-1,i),0,1)      &
                                       )                                       &
-                  ) * rho_air_zw(k) * ddy                                     &
+                  ) * ddy                                                     &
               +   ( w_comp * rho_air(k+1) * ( ibit33 + ibit34 + ibit35 )      &
                 - ( w(k,j,i)   + w(k-1,j,i)   ) * rho_air(k)                  &
@@ -5012,5 +5007,5 @@
                                       + IBITS(advc_flags_2(k-1,j,i),3,1)      &
                                       )                                       & 
-                  ) * ddzu(k+1)   &
+                  ) * drho_air_zw(k) * ddzu(k+1)                              &
                 ) * 0.5_wp
 
@@ -5116,10 +5111,8 @@
 !--             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  &
-                                                             * 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)                                        &
+                        ) * drho_air_zw(k) * ddzu(k+1)                       &
                       ) * 0.5_wp
 
