Index: palm/trunk/SOURCE/production_e.f90
===================================================================
--- palm/trunk/SOURCE/production_e.f90	(revision 2478)
+++ palm/trunk/SOURCE/production_e.f90	(revision 2508)
@@ -25,4 +25,10 @@
 ! -----------------
 ! $Id$
+! - Bugfix in buoyancy production term, wrong base state was set.
+! - Consider use_single_reference_value case if humidity is used.
+! - In case of use_top_fluxes, use correct inverse density at model top
+! - Consider use_surface_fluxes and use_top_fluxes in ocean case  
+! 
+! 2478 2017-09-18 13:37:24Z suehring
 ! Bugfix, consider case where no constant-flux layer and no surfaces fluxes 
 ! are used
@@ -152,5 +158,5 @@
                   humidity, kappa, neutral, ocean, pt_reference,               &
                   rho_reference, use_single_reference_value,                   &
-                  use_surface_fluxes, use_top_fluxes
+                  use_surface_fluxes, use_top_fluxes, vpt_reference
 
        USE grid_variables,                                                     &
@@ -372,10 +378,5 @@
                 DO  m = surf_s, surf_e
                    k = surf_lsm_h%k(m)
-!
-!--                Please note, actually, an interpolation of u_0 and v_0 
-!--                onto the grid center would be required. However, this 
-!--                would require several data transfers between 2D-grid and
-!--                wall type. The effect of this missing interpolation is 
-!--                negligible. (See also production_e_init).
+
                    dudz(k,j) = ( u(k+1,j,i) - surf_lsm_h%u_0(m) ) * dd2zu(k)   
                    dvdz(k,j) = ( v(k+1,j,i) - surf_lsm_h%v_0(m) ) * dd2zu(k)
@@ -388,10 +389,5 @@
                 DO  m = surf_s, surf_e
                    k = surf_usm_h%k(m)
-!
-!--                Please note, actually, an interpolation of u_0 and v_0 
-!--                onto the grid center would be required. However, this 
-!--                would require several data transfers between 2D-grid and
-!--                wall type. The effect of this missing interpolation is 
-!--                negligible. (See also production_e_init).
+
                    dudz(k,j) = ( u(k+1,j,i) - surf_usm_h%u_0(m) ) * dd2zu(k)   
                    dvdz(k,j) = ( v(k+1,j,i) - surf_usm_h%v_0(m) ) * dd2zu(k)
@@ -405,15 +401,9 @@
                 DO  m = surf_s, surf_e
                    k = surf_def_h(1)%k(m)
-!
-!--                Please note, actually, an interpolation of u_0 and v_0 
-!--                onto the grid center would be required. However, this 
-!--                would require several data transfers between 2D-grid and
-!--                wall type. The effect of this missing interpolation is 
-!--                negligible. (See also production_e_init).
+
                    dudz(k,j) = ( surf_def_h(1)%u_0(m) - u(k-1,j,i) ) * dd2zu(k)   
                    dvdz(k,j) = ( surf_def_h(1)%v_0(m) - v(k-1,j,i) ) * dd2zu(k)
 
                 ENDDO
-
              ENDDO
 
@@ -493,175 +483,132 @@
              IF ( .NOT. humidity )  THEN
 
-                IF ( use_single_reference_value )  THEN
-
-                   IF ( ocean )  THEN
-!
-!--                   So far in the ocean no special treatment of density flux
-!--                   in the bottom and top surface layer
-                      DO  j = nys, nyn
-                         DO  k = nzb+1, nzt
-                            tend(k,j,i) = tend(k,j,i) +                        &
-                                          kh(k,j,i) * g / rho_reference *      &
-                                          ( prho(k+1,j,i) - prho(k-1,j,i) ) *  &
-                                          dd2zu(k) *                           &
-                                    MERGE( 1.0_wp, 0.0_wp,                     &
-                                             BTEST( wall_flags_0(k,j,i), 0 )   &
-                                          )
+                IF ( ocean )  THEN
+!
+!--                So far in the ocean no special treatment of density flux
+!--                in the bottom and top surface layer
+                   DO  j = nys, nyn
+                      DO  k = nzb+1, nzt
+                         tend(k,j,i) = tend(k,j,i) +                           &
+                                       kh(k,j,i) * g /                         &
+                              MERGE( rho_reference, prho(k,j,i),               &
+                                     use_single_reference_value ) *            &
+                                       ( prho(k+1,j,i) - prho(k-1,j,i) ) *     &
+                                       dd2zu(k) *                              &
+                                   MERGE( 1.0_wp, 0.0_wp,                      &
+                                          BTEST( wall_flags_0(k,j,i), 30 )     &
+                                        )                            *         &
+                                   MERGE( 1.0_wp, 0.0_wp,                      &
+                                          BTEST( wall_flags_0(k,j,i), 9 )      &
+                                        )  
+                      ENDDO
+!
+!--                   Treatment of near-surface grid points, at up- and down-
+!--                   ward facing surfaces
+                      IF ( use_surface_fluxes )  THEN
+                         DO  l = 0, 1
+                            surf_s = surf_def_h(l)%start_index(j,i)
+                            surf_e = surf_def_h(l)%end_index(j,i)
+                            DO  m = surf_s, surf_e
+                               k = surf_def_h(l)%k(m)
+                               tend(k,j,i) = tend(k,j,i) + g /                 &
+                                         MERGE( rho_reference, prho(k,j,i),    &
+                                                use_single_reference_value ) * &
+                                         drho_air_zw(k-1) *                    &
+                                         surf_def_h(l)%shf(m)
+                            ENDDO
                          ENDDO
+
+                      ENDIF
+
+                      IF ( use_top_fluxes )  THEN
+                         surf_s = surf_def_h(2)%start_index(j,i)
+                         surf_e = surf_def_h(2)%end_index(j,i)
+                         DO  m = surf_s, surf_e
+                            k = surf_def_h(2)%k(m)
+                            tend(k,j,i) = tend(k,j,i) + g /                    &
+                                         MERGE( rho_reference, prho(k,j,i),    &
+                                                use_single_reference_value ) * &
+                                         drho_air_zw(k) *                      &
+                                         surf_def_h(2)%shf(m) 
+                         ENDDO
+                      ENDIF
+
+                   ENDDO
+
+                ELSE
+
+                   DO  j = nys, nyn
+                      DO  k = nzb+1, nzt
+!
+!--                      Flag 9 is used to mask top fluxes, flag 30 to mask
+!--                      surface fluxes
+                         tend(k,j,i) = tend(k,j,i) -                           &
+                                       kh(k,j,i) * g /                         &
+                                   MERGE( pt_reference, pt(k,j,i),             &
+                                           use_single_reference_value ) *      &
+                                       ( pt(k+1,j,i) - pt(k-1,j,i) ) *         &
+                                       dd2zu(k)                      *         &
+                                   MERGE( 1.0_wp, 0.0_wp,                      &
+                                          BTEST( wall_flags_0(k,j,i), 30 )     &
+                                        )                            *         &
+                                   MERGE( 1.0_wp, 0.0_wp,                      &
+                                          BTEST( wall_flags_0(k,j,i), 9 )      &
+                                        )  
                       ENDDO
 
-                   ELSE
-
-                      DO  j = nys, nyn
-                         DO  k = nzb+1, nzt
-!
-!--                         Flag 9 is used to mask top fluxes, flag 30 to mask
-!--                         surface fluxes
-                            tend(k,j,i) = tend(k,j,i) -                        &
-                                          kh(k,j,i) * g / pt_reference  *      &
-                                          ( pt(k+1,j,i) - pt(k-1,j,i) ) *      &
-                                          dd2zu(k)                      *      &
-                                      MERGE( 1.0_wp, 0.0_wp,                   &
-                                             BTEST( wall_flags_0(k,j,i), 30 )  &
-                                           )                            *      &
-                                      MERGE( 1.0_wp, 0.0_wp,                   &
-                                             BTEST( wall_flags_0(k,j,i), 9 )   &
-                                           )  
+                      IF ( use_surface_fluxes )  THEN
+!
+!--                      Default surfaces, up- and downward-facing
+                         DO  l = 0, 1
+                            surf_s = surf_def_h(l)%start_index(j,i)
+                            surf_e = surf_def_h(l)%end_index(j,i)
+                            DO  m = surf_s, surf_e
+                               k = surf_def_h(l)%k(m)
+                               tend(k,j,i) = tend(k,j,i) + g /                 &
+                                    MERGE( pt_reference, pt(k,j,i),            &
+                                           use_single_reference_value )        &
+                                                      * drho_air_zw(k-1)       &
+                                                      * surf_def_h(l)%shf(m)    
+                            ENDDO      
                          ENDDO
-
-                         IF ( use_surface_fluxes )  THEN
-!
-!--                         Default surfaces, up- and downward-facing
-                            DO  l = 0, 1
-                               surf_s = surf_def_h(l)%start_index(j,i)
-                               surf_e = surf_def_h(l)%end_index(j,i)
-                               DO  m = surf_s, surf_e
-                                  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   
-                            ENDDO      
-!
-!--                         Natural surfaces
-                            surf_s = surf_lsm_h%start_index(j,i)
-                            surf_e = surf_lsm_h%end_index(j,i)
-                            DO  m = surf_s, surf_e
-                               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
-!
-!--                         Urban surfaces
-                            surf_s = surf_usm_h%start_index(j,i)
-                            surf_e = surf_usm_h%end_index(j,i)
-                            DO  m = surf_s, surf_e
-                               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                          
-                         ENDIF
-
-                         IF ( use_top_fluxes )  THEN
-                            surf_s = surf_def_h(2)%start_index(j,i)
-                            surf_e = surf_def_h(2)%end_index(j,i)
-                            DO  m = surf_s, surf_e
-                               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
-                         ENDIF
-                      ENDDO
-
-                   ENDIF
-
-                ELSE
-
-                   IF ( ocean )  THEN
-!
-!--                   So far in the ocean no special treatment of density flux
-!--                   in the bottom and top surface layer
-                      DO  j = nys, nyn
-                         DO  k = nzb+1, nzt
-                            tend(k,j,i) = tend(k,j,i) +                        &
-                                          kh(k,j,i) * g / prho(k,j,i) *        &
-                                          ( prho(k+1,j,i) - prho(k-1,j,i) ) *  &
-                                          dd2zu(k)                           * &
-                                     MERGE( 1.0_wp, 0.0_wp,                    &
-                                             BTEST( wall_flags_0(k,j,i), 0 )   &
-                                          )
+!
+!--                      Natural surfaces
+                         surf_s = surf_lsm_h%start_index(j,i)
+                         surf_e = surf_lsm_h%end_index(j,i)
+                         DO  m = surf_s, surf_e
+                            k = surf_lsm_h%k(m)
+                            tend(k,j,i) = tend(k,j,i) + g /                    &
+                                    MERGE( pt_reference, pt(k,j,i),            &
+                                           use_single_reference_value )        &
+                                                      * drho_air_zw(k-1)       &
+                                                      * surf_lsm_h%shf(m)   
                          ENDDO
-                      ENDDO
-
-                   ELSE
-
-                      DO  j = nys, nyn
-                         DO  k = nzb+1, nzt
-!
-!--                         Flag 9 is used to mask top fluxes, flag 30 to mask
-!--                         surface fluxes
-                            tend(k,j,i) = tend(k,j,i) -                        &
-                                          kh(k,j,i) * g / pt(k,j,i) *          &
-                                          ( pt(k+1,j,i) - pt(k-1,j,i) ) *      &
-                                          dd2zu(k)                      *      &
-                                      MERGE( 1.0_wp, 0.0_wp,                   &
-                                             BTEST( wall_flags_0(k,j,i), 30 )  &
-                                           )                            *      &
-                                      MERGE( 1.0_wp, 0.0_wp,                   &
-                                             BTEST( wall_flags_0(k,j,i), 9 )   &
-                                           ) 
+!
+!--                      Urban surfaces
+                         surf_s = surf_usm_h%start_index(j,i)
+                         surf_e = surf_usm_h%end_index(j,i)
+                         DO  m = surf_s, surf_e
+                            k = surf_usm_h%k(m)
+                            tend(k,j,i) = tend(k,j,i) + g /                    &
+                                    MERGE( pt_reference, pt(k,j,i),            &
+                                           use_single_reference_value )        &
+                                                      * drho_air_zw(k-1)       &
+                                                      * surf_usm_h%shf(m)   
+                         ENDDO                          
+                      ENDIF
+
+                      IF ( use_top_fluxes )  THEN
+                         surf_s = surf_def_h(2)%start_index(j,i)
+                         surf_e = surf_def_h(2)%end_index(j,i)
+                         DO  m = surf_s, surf_e
+                            k = surf_def_h(2)%k(m)
+                            tend(k,j,i) = tend(k,j,i) + g /                    &
+                                    MERGE( pt_reference, pt(k,j,i),            &
+                                           use_single_reference_value )        &
+                                                      * drho_air_zw(k)         &
+                                                      * surf_def_h(2)%shf(m) 
                          ENDDO
-
-                         IF ( use_surface_fluxes )  THEN
-!
-!--                         Default surfaces, up- and downwrd-facing
-                            DO  l = 0, 1
-                               surf_s = surf_def_h(l)%start_index(j,i)
-                               surf_e = surf_def_h(l)%end_index(j,i)
-                               DO  m = surf_s, surf_e
-                                  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  
-                            ENDDO
-!
-!--                         Natural surfaces
-                            surf_s = surf_lsm_h%start_index(j,i)
-                            surf_e = surf_lsm_h%end_index(j,i)
-                            DO  m = surf_s, surf_e
-                               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  
-!
-!--                         Urban surfaces
-                            surf_s = surf_usm_h%start_index(j,i)
-                            surf_e = surf_usm_h%end_index(j,i)
-                            DO  m = surf_s, surf_e
-                               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  
-                         ENDIF
-
-                         IF ( use_top_fluxes )  THEN
-                            surf_s = surf_def_h(2)%start_index(j,i)
-                            surf_e = surf_def_h(2)%end_index(j,i)
-                            DO  m = surf_s, surf_e
-                               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
-                         ENDIF
-                      ENDDO
-
-                   ENDIF
+                      ENDIF
+                   ENDDO
 
                 ENDIF
@@ -679,5 +626,7 @@
                          k2 = 0.61_wp * pt(k,j,i)
                          tend(k,j,i) = tend(k,j,i) - kh(k,j,i) *               &
-                                         g / vpt(k,j,i) *                      &
+                                         g /                                   &
+                                    MERGE( vpt_reference, vpt(k,j,i),          &
+                                           use_single_reference_value ) *      &
                                          ( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) +  &
                                            k2 * ( q(k+1,j,i) - q(k-1,j,i) )    &
@@ -704,5 +653,7 @@
                          ENDIF
                          tend(k,j,i) = tend(k,j,i) - kh(k,j,i) *               &
-                                         g / vpt(k,j,i) *                      &
+                                         g /                                   &
+                                    MERGE( vpt_reference, vpt(k,j,i),          &
+                                           use_single_reference_value ) *      &
                                          ( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) +  &
                                            k2 * ( q(k+1,j,i) - q(k-1,j,i) )    &
@@ -718,5 +669,7 @@
                          k2 = 0.61_wp * pt(k,j,i)
                          tend(k,j,i) = tend(k,j,i) -                           & 
-                                       kh(k,j,i) * g / vpt(k,j,i) *            &
+                                       kh(k,j,i) * g /                         &
+                                    MERGE( vpt_reference, vpt(k,j,i),          &
+                                           use_single_reference_value ) *      &
                                        ( k1 * ( pt(k+1,j,i)- pt(k-1,j,i) ) +   &
                                          k2 * ( q(k+1,j,i) -  q(k-1,j,i) ) -   &
@@ -739,5 +692,5 @@
                    DO  j = nys, nyn
 !
-!--                   Treat horizontal default surfaces, up- and downward-facing
+!--                   Treat horizontal default surfaces
                       DO  l = 0, 1
                          surf_s = surf_def_h(l)%start_index(j,i)
@@ -768,7 +721,9 @@
                             ENDIF
 
-                            tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) *       &
+                            tend(k,j,i) = tend(k,j,i) + g /                    &
+                                    MERGE( vpt_reference, vpt(k,j,i),          &
+                                           use_single_reference_value ) *      &
                                                ( 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
@@ -803,5 +758,7 @@
                          ENDIF
 
-                         tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) *          &
+                         tend(k,j,i) = tend(k,j,i) + g /                       &
+                                    MERGE( vpt_reference, vpt(k,j,i),          &
+                                           use_single_reference_value ) *      &
                                                ( k1 * surf_lsm_h%shf(m) +      &
                                                  k2 * surf_lsm_h%qsws(m)       &
@@ -837,5 +794,7 @@
                          ENDIF
 
-                         tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) *          &
+                         tend(k,j,i) = tend(k,j,i) + g /                       &
+                                    MERGE( vpt_reference, vpt(k,j,i),          &
+                                           use_single_reference_value ) *      &
                                                ( k1 * surf_usm_h%shf(m) +      &
                                                  k2 * surf_usm_h%qsws(m)       &
@@ -878,8 +837,10 @@
                          ENDIF
 
-                         tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) *          &
+                         tend(k,j,i) = tend(k,j,i) + g /                       &
+                                    MERGE( vpt_reference, vpt(k,j,i),          &
+                                           use_single_reference_value ) *      &
                                                ( k1 * surf_def_h(2)%shf(m) +   &
-                                                 k2 * surf_def_h(2)%qsws(m)       &
-                                               ) * drho_air_zw(k-1)
+                                                 k2 * surf_def_h(2)%qsws(m)    &
+                                               ) * drho_air_zw(k)
 
                       ENDDO
@@ -916,5 +877,5 @@
                   humidity, kappa, neutral, ocean, pt_reference,               &
                   rho_reference, use_single_reference_value,                   &
-                  use_surface_fluxes, use_top_fluxes
+                  use_surface_fluxes, use_top_fluxes, vpt_reference
 
        USE grid_variables,                                                     &
@@ -1065,5 +1026,5 @@
 !--             -1.0 for right-facing wall, 1.0 for left-facing wall
                 sign_dir = MERGE( 1.0_wp, -1.0_wp,                             &
-                                  BTEST( wall_flags_0(k,j,i-1), 0 ) ) 
+                                  BTEST( wall_flags_0(k,j,i-1), 0 ) )
                 dvdx(k) = sign_dir * vsus / ( km_neutral + 1E-10_wp )
                 dwdx(k) = sign_dir * wsus / ( km_neutral + 1E-10_wp )
@@ -1129,10 +1090,5 @@
           DO  m = surf_s, surf_e
              k = surf_lsm_h%k(m)
-!
-!--          Please note, actually, an interpolation of u_0 and v_0 
-!--          onto the grid center would be required. However, this 
-!--          would require several data transfers between 2D-grid and
-!--          wall type. The effect of this missing interpolation is 
-!--          negligible. (See also production_e_init).
+
              dudz(k)     = ( u(k+1,j,i) - surf_lsm_h%u_0(m) ) * dd2zu(k)   
              dvdz(k)     = ( v(k+1,j,i) - surf_lsm_h%v_0(m) ) * dd2zu(k)
@@ -1144,10 +1100,5 @@
           DO  m = surf_s, surf_e
              k = surf_usm_h%k(m)
-!
-!--          Please note, actually, an interpolation of u_0 and v_0 
-!--          onto the grid center would be required. However, this 
-!--          would require several data transfers between 2D-grid and
-!--          wall type. The effect of this missing interpolation is 
-!--          negligible. (See also production_e_init).
+
              dudz(k)     = ( u(k+1,j,i) - surf_usm_h%u_0(m) ) * dd2zu(k)   
              dvdz(k)     = ( v(k+1,j,i) - surf_usm_h%v_0(m) ) * dd2zu(k)
@@ -1160,10 +1111,5 @@
           DO  m = surf_s, surf_e
              k = surf_def_h(1)%k(m)
-!
-!--          Please note, actually, an interpolation of u_0 and v_0 
-!--          onto the grid center would be required. However, this 
-!--          would require several data transfers between 2D-grid and
-!--          wall type. The effect of this missing interpolation is 
-!--          negligible. (See also production_e_init).
+
              dudz(k)     = ( surf_def_h(1)%u_0(m) - u(k-1,j,i) ) * dd2zu(k)   
              dvdz(k)     = ( surf_def_h(1)%v_0(m) - v(k-1,j,i) ) * dd2zu(k) 
@@ -1232,166 +1178,127 @@
           IF ( .NOT. humidity )  THEN
 
-             IF ( use_single_reference_value )  THEN
-
-                IF ( ocean )  THEN
-!
-!--                So far in the ocean no special treatment of density flux in
-!--                the bottom and top surface layer
-                   DO  k = nzb+1, nzt
+             IF ( ocean )  THEN
+!
+!--             So far in the ocean no special treatment of density flux in
+!--             the bottom and top surface layer
+                DO  k = nzb+1, nzt
  
-                      tend(k,j,i) = tend(k,j,i) +                              &
-                                    kh(k,j,i) * g / rho_reference *            &
-                                    ( prho(k+1,j,i) - prho(k-1,j,i) ) *        &
-                                    dd2zu(k) *                                 &
-                                      MERGE( 1.0_wp, 0.0_wp,                   &
-                                             BTEST( wall_flags_0(k,j,i), 0 )   &
-                                           )                   
+                   tend(k,j,i) = tend(k,j,i) +                                 &
+                                 kh(k,j,i) * g /                               &
+                                 MERGE( rho_reference, prho(k,j,i),            &
+                                        use_single_reference_value ) *         &
+                                 ( prho(k+1,j,i) - prho(k-1,j,i) ) *           &
+                                 dd2zu(k) *                                    &
+                                   MERGE( 1.0_wp, 0.0_wp,                      &
+                                          BTEST( wall_flags_0(k,j,i), 30 )     &
+                                        ) *                                    &
+                                   MERGE( 1.0_wp, 0.0_wp,                      &
+                                          BTEST( wall_flags_0(k,j,i), 9 )      &
+                                        )           
+                ENDDO
+
+                IF ( use_surface_fluxes )  THEN
+!
+!--                Default surfaces, up- and downward-facing
+                   DO  l = 0, 1
+                      surf_s = surf_def_h(l)%start_index(j,i)
+                      surf_e = surf_def_h(l)%end_index(j,i)
+                      DO  m = surf_s, surf_e
+                         k = surf_def_h(l)%k(m)
+                         tend(k,j,i) = tend(k,j,i) + g /                       &
+                                   MERGE( rho_reference, prho(k,j,i),          &
+                                          use_single_reference_value ) *       &
+                                   drho_air_zw(k-1) *                          &
+                                   surf_def_h(l)%shf(m)
+                      ENDDO
                    ENDDO
 
-                ELSE
-
-                   DO  k = nzb+1, nzt
-!
-!--                   Flag 9 is used to mask top fluxes, flag 30 to mask
-!--                   surface fluxes
-                      tend(k,j,i) = tend(k,j,i) -                              &
-                                    kh(k,j,i) * g / pt_reference *             &
-                                    ( pt(k+1,j,i) - pt(k-1,j,i) ) * dd2zu(k) * &
-                                      MERGE( 1.0_wp, 0.0_wp,                   &
-                                             BTEST( wall_flags_0(k,j,i), 30 )  &
-                                           ) *                                 &
-                                      MERGE( 1.0_wp, 0.0_wp,                   &
-                                             BTEST( wall_flags_0(k,j,i), 9 )   &
-                                           )
-
+                ENDIF
+
+                IF ( use_top_fluxes )  THEN
+                   surf_s = surf_def_h(2)%start_index(j,i)
+                   surf_e = surf_def_h(2)%end_index(j,i)
+                   DO  m = surf_s, surf_e
+                      k = surf_def_h(2)%k(m)
+                      tend(k,j,i) = tend(k,j,i) + g /                          &
+                                   MERGE( rho_reference, prho(k,j,i),          &
+                                          use_single_reference_value ) *       &
+                                   drho_air_zw(k) *                            &
+                                   surf_def_h(2)%shf(m) 
                    ENDDO
-
-                   IF ( use_surface_fluxes )  THEN
-!
-!--                   Default surfaces, up- and downward-facing
-                      DO  l = 0, 1
-                         surf_s = surf_def_h(l)%start_index(j,i)
-                         surf_e = surf_def_h(l)%end_index(j,i)
-                         DO  m = surf_s, surf_e
-                            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
+                ENDIF
+
+             ELSE
+
+                DO  k = nzb+1, nzt
+!
+!--                Flag 9 is used to mask top fluxes, flag 30 to mask
+!--                surface fluxes
+                   tend(k,j,i) = tend(k,j,i) -                                 &
+                                 kh(k,j,i) * g /                               &
+                                   MERGE( pt_reference, pt(k,j,i),             &
+                                          use_single_reference_value ) *       &
+                                 ( pt(k+1,j,i) - pt(k-1,j,i) ) * dd2zu(k) *    &
+                                   MERGE( 1.0_wp, 0.0_wp,                      &
+                                          BTEST( wall_flags_0(k,j,i), 30 )     &
+                                        ) *                                    &
+                                   MERGE( 1.0_wp, 0.0_wp,                      &
+                                          BTEST( wall_flags_0(k,j,i), 9 )      &
+                                        )
+
+                ENDDO
+
+                IF ( use_surface_fluxes )  THEN
+!
+!--                Default surfaces, up- and downward-facing
+                   DO  l = 0, 1
+                      surf_s = surf_def_h(l)%start_index(j,i)
+                      surf_e = surf_def_h(l)%end_index(j,i)
+                      DO  m = surf_s, surf_e
+                         k = surf_def_h(l)%k(m)
+                         tend(k,j,i) = tend(k,j,i) + g /                       &
+                                   MERGE( pt_reference, pt(k,j,i),             &
+                                          use_single_reference_value ) *       &
+                                   drho_air_zw(k-1) *                          &
+                                   surf_def_h(l)%shf(m)
                       ENDDO
-!
-!--                   Natural surfaces
-                      surf_s = surf_lsm_h%start_index(j,i)
-                      surf_e = surf_lsm_h%end_index(j,i)
-                      DO  m = surf_s, surf_e
-                         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
-!
-!--                   Urban surfaces
-                      surf_s = surf_usm_h%start_index(j,i)
-                      surf_e = surf_usm_h%end_index(j,i)
-                      DO  m = surf_s, surf_e
-                         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
-                   ENDIF
-
-                   IF ( use_top_fluxes )  THEN
-                      surf_s = surf_def_h(2)%start_index(j,i)
-                      surf_e = surf_def_h(2)%end_index(j,i)
-                      DO  m = surf_s, surf_e
-                         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
-                   ENDIF
-
+                   ENDDO
+!
+!--                Natural surfaces
+                   surf_s = surf_lsm_h%start_index(j,i)
+                   surf_e = surf_lsm_h%end_index(j,i)
+                   DO  m = surf_s, surf_e
+                      k = surf_lsm_h%k(m)
+                      tend(k,j,i) = tend(k,j,i) + g /                          &
+                                   MERGE( pt_reference, pt(k,j,i),             &
+                                          use_single_reference_value ) *       &
+                                   drho_air_zw(k-1) *                          &
+                                   surf_lsm_h%shf(m) 
+                   ENDDO
+!
+!--                Urban surfaces
+                   surf_s = surf_usm_h%start_index(j,i)
+                   surf_e = surf_usm_h%end_index(j,i)
+                   DO  m = surf_s, surf_e
+                      k = surf_usm_h%k(m)
+                      tend(k,j,i) = tend(k,j,i) + g /                          &
+                                   MERGE( pt_reference, pt(k,j,i),             &
+                                          use_single_reference_value ) *       &
+                                   drho_air_zw(k-1) *                          &
+                                   surf_usm_h%shf(m) 
+                   ENDDO
                 ENDIF
 
-             ELSE
-
-                IF ( ocean )  THEN
-!
-!--                So far in the ocean no special treatment of density flux in
-!--                the bottom and top surface layer
-                   DO  k = nzb+1, nzt
-                      tend(k,j,i) = tend(k,j,i) +                              &
-                                    kh(k,j,i) * g / prho(k,j,i)       *        &
-                                    ( prho(k+1,j,i) - prho(k-1,j,i) ) *        &
-                                    dd2zu(k) *                                 &
-                                      MERGE( 1.0_wp, 0.0_wp,                   &
-                                             BTEST( wall_flags_0(k,j,i), 0 )   &
-                                           )           
+                IF ( use_top_fluxes )  THEN
+                   surf_s = surf_def_h(2)%start_index(j,i)
+                   surf_e = surf_def_h(2)%end_index(j,i)
+                   DO  m = surf_s, surf_e
+                      k = surf_def_h(2)%k(m)
+                      tend(k,j,i) = tend(k,j,i) + g /                          &
+                                   MERGE( pt_reference, pt(k,j,i),             &
+                                          use_single_reference_value ) *       &
+                                   drho_air_zw(k) *                            &
+                                   surf_def_h(2)%shf(m) 
                    ENDDO
-
-                ELSE
-
-                   DO  k = nzb+1, nzt
-!
-!--                   Flag 9 is used to mask top fluxes, flag 30 to mask
-!--                   surface fluxes
-                      tend(k,j,i) = tend(k,j,i) -                              &
-                                    kh(k,j,i) * g / pt(k,j,i) *                &
-                                    ( pt(k+1,j,i) - pt(k-1,j,i) ) * dd2zu(k) * &
-                                      MERGE( 1.0_wp, 0.0_wp,                   &
-                                             BTEST( wall_flags_0(k,j,i), 30 )  &
-                                           ) *                                 &
-                                      MERGE( 1.0_wp, 0.0_wp,                   &
-                                             BTEST( wall_flags_0(k,j,i), 9 )   &
-                                           )
-                   ENDDO
-
-                   IF ( use_surface_fluxes )  THEN
-!
-!--                   Default surfaces, up- and downward-facing
-                      DO  l = 0, 1
-                         surf_s = surf_def_h(l)%start_index(j,i)
-                         surf_e = surf_def_h(l)%end_index(j,i)
-                         DO  m = surf_s, surf_e
-                            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  
-                      ENDDO
-!
-!--                   Natural surfaces
-                      surf_s = surf_lsm_h%start_index(j,i)
-                      surf_e = surf_lsm_h%end_index(j,i)
-                      DO  m = surf_s, surf_e
-                         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  
-!
-!--                   Urban surfaces
-                      surf_s = surf_usm_h%start_index(j,i)
-                      surf_e = surf_usm_h%end_index(j,i)
-                      DO  m = surf_s, surf_e
-                         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
-                   ENDIF
-
-                   IF ( use_top_fluxes )  THEN
-                      surf_s = surf_def_h(2)%start_index(j,i)
-                      surf_e = surf_def_h(2)%end_index(j,i)
-                      DO  m = surf_s, surf_e
-                         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
-                   ENDIF
-
                 ENDIF
 
@@ -1406,5 +1313,7 @@
                    k1 = 1.0_wp + 0.61_wp * q(k,j,i)
                    k2 = 0.61_wp * pt(k,j,i)
-                   tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / vpt(k,j,i) *    &
+                   tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g /                 &
+                                   MERGE( vpt_reference, vpt(k,j,i),           &
+                                          use_single_reference_value ) *       &
                                          ( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) +  &
                                            k2 * ( q(k+1,j,i) - q(k-1,j,i) )    &
@@ -1431,5 +1340,7 @@
                       k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp )
                    ENDIF
-                   tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / vpt(k,j,i) *    &
+                   tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g /                 &
+                                   MERGE( vpt_reference, vpt(k,j,i),           &
+                                          use_single_reference_value ) *       &
                                          ( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) +  &
                                            k2 * ( q(k+1,j,i) - q(k-1,j,i) )    &
@@ -1444,5 +1355,7 @@
                    k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
                    k2 = 0.61_wp * pt(k,j,i)
-                   tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / vpt(k,j,i) *    &
+                   tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g /                 &
+                                   MERGE( vpt_reference, vpt(k,j,i),           &
+                                          use_single_reference_value ) *       &
                                      ( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) +      &
                                        k2 * ( q(k+1,j,i) - q(k-1,j,i) ) -      &
@@ -1489,5 +1402,7 @@
                       ENDIF
 
-                      tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) *             &
+                      tend(k,j,i) = tend(k,j,i) + g /                          &
+                                   MERGE( vpt_reference, vpt(k,j,i),           &
+                                          use_single_reference_value ) *       &
                                          ( k1 * surf_def_h(l)%shf(m) +         &
                                            k2 * surf_def_h(l)%qsws(m)          &
@@ -1524,5 +1439,7 @@
                    ENDIF
 
-                   tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) *                &
+                   tend(k,j,i) = tend(k,j,i) + g /                             &
+                                   MERGE( vpt_reference, vpt(k,j,i),           &
+                                          use_single_reference_value ) *       &
                                             ( k1 * surf_lsm_h%shf(m) +         &
                                               k2 * surf_lsm_h%qsws(m)          &
@@ -1558,5 +1475,7 @@
                    ENDIF
 
-                   tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) *                &
+                   tend(k,j,i) = tend(k,j,i) + g /                             &
+                                   MERGE( vpt_reference, vpt(k,j,i),           &
+                                          use_single_reference_value ) *       &
                                             ( k1 * surf_usm_h%shf(m) +         &
                                               k2 * surf_usm_h%qsws(m)          &
@@ -1596,8 +1515,10 @@
                    ENDIF
 
-                   tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) *                &
+                   tend(k,j,i) = tend(k,j,i) + g /                             &
+                                   MERGE( vpt_reference, vpt(k,j,i),           &
+                                          use_single_reference_value ) *       &
                                ( k1* surf_def_h(2)%shf(m) +                    &
                                  k2 * surf_def_h(2)%qsws(m)                    &
-                               ) * drho_air_zw(k-1)
+                               ) * drho_air_zw(k)
                 ENDDO
 
@@ -1682,5 +1603,5 @@
           ENDDO
 !
-!--       Default surfaces, downward-facing
+!--       Default surfaces, downward-facing surfaces
           !$OMP PARALLEL DO PRIVATE(i,j,k,m)
           DO  m = 1, surf_def_h(1)%ns
@@ -1689,26 +1610,20 @@
              j = surf_def_h(1)%j(m)
              k = surf_def_h(1)%k(m)
-!
-!--          Note, calculatione of u_0 and v_0 is not fully accurate, as u/v 
-!--          and km are not on the same grid. Actually, a further 
-!--          interpolation of km onto the u/v-grid is necessary. However, the
-!--          effect of this error is negligible. 
-!--          In case of downward-facing surfaces, gradient is calculated 
-!--          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)    )  /         &
+
+             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)    )  /         &
+             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 )  
 
-             IF ( ABS( surf_def_h(1)%u_0(m) - u(k-1,j,i) )  >                  &
-                  ABS( u(k+1,j,i)           - u(k-1,j,i) )                     &
+             IF ( ABS( surf_def_h(1)%u_0(m) - u(k-1,j,i) )  >               &
+                  ABS( u(k+1,j,i)           - u(k-1,j,i) )                  &
                 )  surf_def_h(1)%u_0(m) = u(k+1,j,i)
 
-             IF ( ABS( surf_def_h(1)%v_0(m) - v(k-1,j,i) )  >                  &
-                  ABS( v(k+1,j,i)           - v(k-1,j,i) )                     &
+             IF ( ABS( surf_def_h(1)%v_0(m) - v(k-1,j,i) )  >               &
+                  ABS( v(k+1,j,i)           - v(k-1,j,i) )                  &
                 )  surf_def_h(1)%v_0(m) = v(k+1,j,i)
 
