Index: palm/trunk/SOURCE/advec_ws.f90
===================================================================
--- palm/trunk/SOURCE/advec_ws.f90	(revision 2328)
+++ 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
 
Index: palm/trunk/SOURCE/check_parameters.f90
===================================================================
--- palm/trunk/SOURCE/check_parameters.f90	(revision 2328)
+++ palm/trunk/SOURCE/check_parameters.f90	(revision 2329)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Bugfix: index corrected for rho_air and rho_air_zw output
+! 
+! 2320 2017-07-21 12:47:43Z suehring
 ! Modularize large-scale forcing and nudging
 !
@@ -2748,12 +2751,12 @@
 
           CASE ( 'rho_air' )
-             dopr_index(i)  = 121
+             dopr_index(i)  = 119
              dopr_unit(i)   = 'kg/m3'
-             hom(:,2,121,:) = SPREAD( zu, 2, statistic_regions+1 )
+             hom(:,2,119,:) = SPREAD( zu, 2, statistic_regions+1 )
 
           CASE ( 'rho_air_zw' )
-             dopr_index(i)  = 122
+             dopr_index(i)  = 120
              dopr_unit(i)   = 'kg/m3'
-             hom(:,2,122,:) = SPREAD( zw, 2, statistic_regions+1 )
+             hom(:,2,120,:) = SPREAD( zw, 2, statistic_regions+1 )
 
           CASE ( 'nc' )
Index: palm/trunk/SOURCE/init_3d_model.f90
===================================================================
--- palm/trunk/SOURCE/init_3d_model.f90	(revision 2328)
+++ palm/trunk/SOURCE/init_3d_model.f90	(revision 2329)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Removed temporary bugfix (r2327) as bug is properly resolved by this revision
+! 
+! 2327 2017-08-02 07:40:57Z maronga
 ! Temporary bugfix
 ! 
@@ -346,5 +349,4 @@
 !> or
 !> c) read values of a previous run
-!> @todo fix crashes for runs with topography and density /= 1
 !------------------------------------------------------------------------------!
  SUBROUTINE init_3d_model
@@ -737,14 +739,5 @@
        rho_air_zw(nzt+1)  = rho_air_zw(nzt)                                    &
                             + 2.0_wp * ( rho_air(nzt+1) - rho_air_zw(nzt)  )
-
-!
-!--    Temporary workaround as runs with topography do not work with density /= 1
-       rho_air    = 1.0_wp
-       rho_air_zw = 1.0_wp
-
-    ENDIF
-
-
-
+    ENDIF
 
 !-- compute the inverse density array in order to avoid expencive divisions
Index: palm/trunk/SOURCE/production_e.f90
===================================================================
--- palm/trunk/SOURCE/production_e.f90	(revision 2328)
+++ palm/trunk/SOURCE/production_e.f90	(revision 2329)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Bugfix: added division by density as kinematic fluxes are needed
+! 
+! 2233 2017-05-30 18:08:54Z suehring
 !
 ! 2232 2017-05-30 17:47:52Z suehring
@@ -135,5 +138,6 @@
 
        USE arrays_3d,                                                          &
-           ONLY:  ddzw, dd2zu, kh, km, prho, pt, q, ql, tend, u, v, vpt, w
+           ONLY:  ddzw, dd2zu, drho_air_zw, kh, km, prho, pt, q, ql, tend, u,  &
+                  v, vpt, w
 
        USE cloud_parameters,                                                   &
@@ -531,4 +535,5 @@
                                   k = surf_def_h(l)%k(m)
                                   tend(k,j,i) = tend(k,j,i) + g / pt_reference &
+                                                         * drho_air_zw(k-1)    &
                                                          * surf_def_h(l)%shf(m)   
                                ENDDO   
@@ -541,4 +546,5 @@
                                k = surf_lsm_h%k(m)
                                tend(k,j,i) = tend(k,j,i) + g / pt_reference    &
+                                                         * drho_air_zw(k-1)    &
                                                          * surf_lsm_h%shf(m)   
                             ENDDO
@@ -550,4 +556,5 @@
                                k = surf_usm_h%k(m)
                                tend(k,j,i) = tend(k,j,i) + g / pt_reference    &
+                                                         * drho_air_zw(k-1)    &
                                                          * surf_usm_h%shf(m)   
                             ENDDO                          
@@ -559,6 +566,7 @@
                             DO  m = surf_s, surf_e
                                k = surf_def_h(2)%k(m)
-                               tend(k,j,i) = tend(k,j,i) + g / pt_reference *     &
-                                                           surf_def_h(2)%shf(m) 
+                               tend(k,j,i) = tend(k,j,i) + g / pt_reference    &
+                                                         * drho_air_zw(k-1)    &
+                                                         * surf_def_h(2)%shf(m) 
                             ENDDO
                          ENDIF
@@ -613,4 +621,5 @@
                                   k = surf_def_h(l)%k(m)
                                   tend(k,j,i) = tend(k,j,i) + g / pt_reference &
+                                                         * drho_air_zw(k-1)    &
                                                          * surf_def_h(l)%shf(m)   
                                ENDDO  
@@ -623,4 +632,5 @@
                                k = surf_lsm_h%k(m)
                                tend(k,j,i) = tend(k,j,i) + g / pt_reference    &
+                                                         * drho_air_zw(k-1)    &
                                                          * surf_lsm_h%shf(m)   
                             ENDDO  
@@ -632,4 +642,5 @@
                                k = surf_usm_h%k(m)
                                tend(k,j,i) = tend(k,j,i) + g / pt_reference    &
+                                                         * drho_air_zw(k-1)    &
                                                          * surf_usm_h%shf(m)   
                             ENDDO  
@@ -641,6 +652,7 @@
                             DO  m = surf_s, surf_e
                                k = surf_def_h(2)%k(m)
-                               tend(k,j,i) = tend(k,j,i) + g / pt_reference *     &
-                                                           surf_def_h(2)%shf(m) 
+                               tend(k,j,i) = tend(k,j,i) + g / pt_reference    &
+                                                         * drho_air_zw(k-1)    &
+                                                         * surf_def_h(2)%shf(m) 
                             ENDDO
                          ENDIF
@@ -754,5 +766,6 @@
                             tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) *       &
                                                ( k1 * surf_def_h(l)%shf(m) +   &
-                                                 k2 * surf_def_h(l)%qsws(m) )
+                                                 k2 * surf_def_h(l)%qsws(m)       &
+                                               ) * drho_air_zw(k-1)
                          ENDDO
                       ENDDO
@@ -788,5 +801,6 @@
                          tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) *          &
                                                ( k1 * surf_lsm_h%shf(m) +      &
-                                                 k2 * surf_lsm_h%qsws(m) )
+                                                 k2 * surf_lsm_h%qsws(m)       &
+                                               ) * drho_air_zw(k-1)
                       ENDDO
 !
@@ -821,5 +835,6 @@
                          tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) *          &
                                                ( k1 * surf_usm_h%shf(m) +      &
-                                                 k2 * surf_usm_h%qsws(m) )
+                                                 k2 * surf_usm_h%qsws(m)       &
+                                               ) * drho_air_zw(k-1)
                       ENDDO
 
@@ -861,5 +876,6 @@
                          tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) *          &
                                                ( k1 * surf_def_h(2)%shf(m) +   &
-                                                 k2 * surf_def_h(2)%qsws(m) )
+                                                 k2 * surf_def_h(2)%qsws(m)       &
+                                               ) * drho_air_zw(k-1)
 
                       ENDDO
@@ -886,5 +902,6 @@
 
        USE arrays_3d,                                                          &
-           ONLY:  ddzw, dd2zu, kh, km, prho, pt, q, ql, tend, u, v, vpt, w
+           ONLY:  ddzw, dd2zu, drho_air_zw, kh, km, prho, pt, q, ql, tend, u,  &
+                  v, vpt, w
 
        USE cloud_parameters,                                                   &
@@ -1254,5 +1271,6 @@
                             k = surf_def_h(l)%k(m)
                             tend(k,j,i) = tend(k,j,i) + g / pt_reference *     &
-                                                           surf_def_h(l)%shf(m) 
+                                                           drho_air_zw(k-1) *  &
+                                                           surf_def_h(l)%shf(m)
                          ENDDO
                       ENDDO
@@ -1264,4 +1282,5 @@
                          k = surf_lsm_h%k(m)
                          tend(k,j,i) = tend(k,j,i) + g / pt_reference *        &
+                                                           drho_air_zw(k-1) *  &
                                                            surf_lsm_h%shf(m) 
                       ENDDO
@@ -1273,4 +1292,5 @@
                          k = surf_usm_h%k(m)
                          tend(k,j,i) = tend(k,j,i) + g / pt_reference *        &
+                                                           drho_air_zw(k-1) *  &
                                                            surf_usm_h%shf(m) 
                       ENDDO
@@ -1283,4 +1303,5 @@
                          k = surf_def_h(2)%k(m)
                          tend(k,j,i) = tend(k,j,i) + g / pt_reference *        &
+                                                           drho_air_zw(k-1) *  &
                                                            surf_def_h(2)%shf(m) 
                       ENDDO
@@ -1331,5 +1352,6 @@
                             k = surf_def_h(l)%k(m)
                             tend(k,j,i) = tend(k,j,i) + g / pt_reference       &
-                                                   * surf_def_h(l)%shf(m)   
+                                                   * drho_air_zw(k-1)          &
+                                                   * surf_def_h(l)%shf(m)
                          ENDDO  
                       ENDDO
@@ -1341,5 +1363,6 @@
                          k = surf_lsm_h%k(m)
                          tend(k,j,i) = tend(k,j,i) + g / pt_reference          &
-                                                     * surf_lsm_h%shf(m)   
+                                                     * drho_air_zw(k-1)        &
+                                                     * surf_lsm_h%shf(m)
                       ENDDO  
 !
@@ -1350,5 +1373,6 @@
                          k = surf_usm_h%k(m)
                          tend(k,j,i) = tend(k,j,i) + g / pt_reference          &
-                                                     * surf_usm_h%shf(m)   
+                                                     * drho_air_zw(k-1)        &
+                                                     * surf_usm_h%shf(m)
                       ENDDO
                    ENDIF
@@ -1360,5 +1384,6 @@
                          k = surf_def_h(2)%k(m)
                          tend(k,j,i) = tend(k,j,i) + g / pt_reference *        &
-                                                           surf_def_h(2)%shf(m) 
+                                                           drho_air_zw(k-1) *  &
+                                                           surf_def_h(2)%shf(m)
                       ENDDO
                    ENDIF
@@ -1461,5 +1486,6 @@
                       tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) *             &
                                          ( k1 * surf_def_h(l)%shf(m) +         &
-                                           k2 * surf_def_h(l)%qsws(m) )
+                                           k2 * surf_def_h(l)%qsws(m)          &
+                                         ) * drho_air_zw(k-1)
                    ENDDO
                 ENDDO
@@ -1495,5 +1521,6 @@
                    tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) *                &
                                             ( k1 * surf_lsm_h%shf(m) +         &
-                                              k2 * surf_lsm_h%qsws(m) )
+                                              k2 * surf_lsm_h%qsws(m)          &
+                                            ) * drho_air_zw(k-1)
                 ENDDO
 !
@@ -1528,5 +1555,6 @@
                    tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) *                &
                                             ( k1 * surf_usm_h%shf(m) +         &
-                                              k2 * surf_usm_h%qsws(m) )
+                                              k2 * surf_usm_h%qsws(m)          &
+                                            ) * drho_air_zw(k-1)
                 ENDDO
 
@@ -1565,5 +1593,6 @@
                    tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) *                &
                                ( k1* surf_def_h(2)%shf(m) +                    &
-                                 k2 * surf_def_h(2)%qsws(m) )
+                                 k2 * surf_def_h(2)%qsws(m)                    &
+                               ) * drho_air_zw(k-1)
                 ENDDO
 
@@ -1585,5 +1614,5 @@
 
        USE arrays_3d,                                                          &
-           ONLY:  kh, km, u, v, zu
+           ONLY:  kh, km, drho_air_zw, u, v, zu
 
        USE control_parameters,                                                 &
@@ -1630,7 +1659,9 @@
 !--          effect of this error is negligible. 
              surf_def_h(0)%u_0(m) = u(k+1,j,i) + surf_def_h(0)%usws(m) *       &
+                                        drho_air_zw(k-1) *                     &
                                         ( zu(k+1)    - zu(k-1)    )  /         &
                                         ( km(k,j,i)  + 1.0E-20_wp )  
              surf_def_h(0)%v_0(m) = v(k+1,j,i) + surf_def_h(0)%vsws(m) *       &
+                                        drho_air_zw(k-1) *                     &
                                         ( zu(k+1)    - zu(k-1)    )  /         &
                                         ( km(k,j,i)  + 1.0E-20_wp )  
@@ -1661,7 +1692,9 @@
 !--          between u_0 and u(k-1).
              surf_def_h(1)%u_0(m) = u(k-1,j,i) - surf_def_h(1)%usws(m) *       &
+                                        drho_air_zw(k-1) *                     &
                                         ( zu(k+1)    - zu(k-1)    )  /         &
                                         ( km(k,j,i)  + 1.0E-20_wp )  
              surf_def_h(1)%v_0(m) = v(k-1,j,i) - surf_def_h(1)%vsws(m) *       &
+                                        drho_air_zw(k-1) *                     &
                                         ( zu(k+1)    - zu(k-1)    )  /         &
                                         ( km(k,j,i)  + 1.0E-20_wp )  
@@ -1690,7 +1723,9 @@
 !--          effect of this error is negligible. 
              surf_lsm_h%u_0(m) = u(k+1,j,i) + surf_lsm_h%usws(m)      *        &
+                                           drho_air_zw(k-1) *                  &
                                            ( zu(k+1)   - zu(k-1)    )  /       &
                                            ( km(k,j,i) + 1.0E-20_wp )  
              surf_lsm_h%v_0(m) = v(k+1,j,i) + surf_lsm_h%vsws(m)      *        &
+                                           drho_air_zw(k-1) *                  &
                                            ( zu(k+1)   - zu(k-1)    )  /       &
                                            ( km(k,j,i) + 1.0E-20_wp )
@@ -1719,7 +1754,9 @@
 !--          effect of this error is negligible. 
              surf_usm_h%u_0(m) = u(k+1,j,i) + surf_usm_h%usws(m)      *        &
+                                           drho_air_zw(k-1) *                  &
                                            ( zu(k+1)   - zu(k-1)    )  /       &
                                            ( km(k,j,i) + 1.0E-20_wp )  
              surf_usm_h%v_0(m) = v(k+1,j,i) + surf_usm_h%vsws(m)      *        &
+                                           drho_air_zw(k-1) *                  &
                                            ( zu(k+1)   - zu(k-1)    )  /       &
                                            ( km(k,j,i) + 1.0E-20_wp )
