Index: palm/trunk/SOURCE/surface_layer_fluxes_mod.f90
===================================================================
--- palm/trunk/SOURCE/surface_layer_fluxes_mod.f90	(revision 3130)
+++ palm/trunk/SOURCE/surface_layer_fluxes_mod.f90	(revision 3146)
@@ -21,5 +21,5 @@
 ! Current revisions:
 ! ------------------
-! 
+! Major bugfix in calculation of Obukhov length
 ! 
 ! Former revisions:
@@ -329,5 +329,10 @@
              surf => surf_def_h(l)
              CALL calc_pt_q
-             IF ( .NOT. neutral )  CALL calc_pt_surface
+             IF ( .NOT. neutral )  THEN
+                CALL calc_pt_surface
+                IF ( humidity )  THEN
+                   CALL calc_vpt_surface
+                ENDIF
+             ENDIF
           ENDIF
        ENDDO
@@ -354,6 +359,5 @@
           ENDIF
 
-!
-!--    Call for urban-type vertical surfaces
+!--       Call for urban-type vertical surfaces
           IF ( surf_usm_v(l)%ns >= 1  )  THEN
              surf => surf_usm_v(l)
@@ -1184,39 +1188,31 @@
           IF ( ibc_pt_b /= 1 )  THEN
              IF ( humidity )  THEN
-                !$OMP PARALLEL DO PRIVATE( i, j, k, z_mo )
+                !$OMP PARALLEL DO PRIVATE( z_mo )
                 DO  m = 1, surf%ns
 
-                   i   = surf%i(m)            
-                   j   = surf%j(m)
-                   k   = surf%k(m)
-
                    z_mo = surf%z_mo(m)
 
-                   surf%rib(m) = g * z_mo *                                    &
-                                         ( vpt(k,j,i) - vpt(k-1,j,i) ) /       &
-                      ( surf%uvw_abs(m)**2 * vpt(k,j,i) + 1.0E-20_wp )
+                   surf%rib(m) = g * z_mo                                      &
+                                 * ( surf%vpt1(m) - surf%vpt_surface(m) )      &
+                                 / ( surf%uvw_abs(m)**2 * surf%vpt1(m)         &
+                                 + 1.0E-20_wp )
                 ENDDO
              ELSE
-                !$OMP PARALLEL DO PRIVATE( i, j, k, z_mo )
+                !$OMP PARALLEL DO PRIVATE( z_mo )
                 DO  m = 1, surf%ns
 
-                   i   = surf%i(m)            
-                   j   = surf%j(m)
-                   k   = surf%k(m)
-
                    z_mo = surf%z_mo(m)
 
-                   surf%rib(m) = g * z_mo *                                    &
-                                         ( pt(k,j,i) - pt(k-1,j,i)   ) /       &
-                      ( surf%uvw_abs(m)**2 * pt(k,j,i)  + 1.0E-20_wp )
+                   surf%rib(m) = g * z_mo                                      &
+                                 * ( surf%pt1(m) - surf%pt_surface(m) )        &
+                                 / ( surf%uvw_abs(m)**2 * surf%pt1(m)          &
+                                 + 1.0E-20_wp )
                 ENDDO
              ENDIF
           ELSE
              IF ( humidity )  THEN
-                !$OMP PARALLEL DO PRIVATE( i, j, k, z_mo )
+                !$OMP PARALLEL DO PRIVATE( k, z_mo )
                 DO  m = 1, surf%ns
 
-                   i   = surf%i(m)            
-                   j   = surf%j(m)
                    k   = surf%k(m)
 
@@ -1224,16 +1220,14 @@
 
                    surf%rib(m) = - g * z_mo * ( ( 1.0_wp + 0.61_wp             &
-                           * q(k,j,i) ) * surf%shf(m) + 0.61_wp                &
-                           * pt(k,j,i) * surf%qsws(m) ) *                      &
+                           * surf%qv1(m) ) * surf%shf(m) + 0.61_wp             &
+                           * surf%pt1(m) * surf%qsws(m) ) *                    &
                              drho_air_zw(k-1)                /                 &
-                         ( surf%uvw_abs(m)**3 * vpt(k,j,i) * kappa**2          &
+                         ( surf%uvw_abs(m)**3 * surf%vpt1(m) * kappa**2        &
                            + 1.0E-20_wp )
                 ENDDO
              ELSE
-                !$OMP PARALLEL DO PRIVATE( i, j, k, z_mo )
+                !$OMP PARALLEL DO PRIVATE( k, z_mo )
                 DO  m = 1, surf%ns
 
-                   i   = surf%i(m)            
-                   j   = surf%j(m)
                    k   = surf%k(m)
 
@@ -1242,5 +1236,5 @@
                    surf%rib(m) = - g * z_mo * surf%shf(m) *                    &
                                         drho_air_zw(k-1)            /          &
-                        ( surf%uvw_abs(m)**3 * pt(k,j,i) * kappa**2            &
+                        ( surf%uvw_abs(m)**3 * surf%pt1(m) * kappa**2          &
                            + 1.0E-20_wp )
                 ENDDO
@@ -1418,14 +1412,10 @@
 
           IF ( .NOT. humidity )  THEN
-             !$OMP PARALLEL DO PRIVATE( i, j, k, z_mo )
+             !$OMP PARALLEL DO PRIVATE( z_mo )
              DO  m = 1, surf%ns
 
-                i   = surf%i(m)            
-                j   = surf%j(m)
-                k   = surf%k(m)
-
                 z_mo = surf%z_mo(m)
 
-                surf%ol(m) =  ( pt(k,j,i) *  surf%us(m)**2 ) /                 &
+                surf%ol(m) =  ( surf%pt1(m) *  surf%us(m)**2 ) /               &
                                        ( kappa * g *                           &
                                          surf%ts(m) + 1E-30_wp )
@@ -1453,5 +1443,5 @@
 
 
-                surf%ol(m) =  ( vpt(k,j,i) * surf%us(m)**2 ) /                 &
+                surf%ol(m) =  ( surf%vpt1(m) * surf%us(m)**2 ) /               &
                     ( kappa * g * ( surf%ts(m) +                               &
                          0.61_wp * surf%pt1(m) * surf%us(m)                    &
@@ -1472,16 +1462,12 @@
           ELSE
 
-             !$OMP PARALLEL DO PRIVATE( i, j, k, z_mo )
+             !$OMP PARALLEL DO PRIVATE( z_mo )
              DO  m = 1, surf%ns
 
-                i   = surf%i(m)            
-                j   = surf%j(m)
-                k   = surf%k(m)
-
                 z_mo = surf%z_mo(m)
 
-                surf%ol(m) =  ( vpt(k,j,i) *  surf%us(m)**2 ) /                &
-                    ( kappa * g * ( surf%ts(m) + 0.61_wp * pt(k,j,i) *         &
-                                    surf%qs(m) + 0.61_wp * q(k,j,i)  *         &
+                surf%ol(m) =  ( surf%vpt1(m) *  surf%us(m)**2 ) /              &
+                    ( kappa * g * ( surf%ts(m) + 0.61_wp * surf%pt1(m) *       &
+                                    surf%qs(m) + 0.61_wp * surf%qv1(m) *       &
                                     surf%ts(m) ) + 1E-30_wp )
 
@@ -1560,5 +1546,6 @@
 
 !
-!-- Calculate potential temperature and specific humidity at first grid level
+!-- Calculate potential temperature, specific humidity, and virtual potential
+!-- temperature at first grid level
 !-- ( only for upward-facing surfs )
     SUBROUTINE calc_pt_q
@@ -1590,4 +1577,8 @@
           ENDIF 
 
+          IF ( humidity )  THEN
+             surf%vpt1(m) = pt(k,j,i) * ( 1.0_wp + 0.61_wp * q(k,j,i) )
+          ENDIF
+          
        ENDDO
 
@@ -1596,5 +1587,5 @@
 
 !
-!-- Calculate potential temperature and specific humidity at first grid level
+!-- Calculate potential temperature at first grid level
 !-- ( only for upward-facing surfs )
     SUBROUTINE calc_pt_surface
@@ -1602,8 +1593,8 @@
        IMPLICIT NONE
 
-       INTEGER(iwp) ::  koff    !< index offset between surface and atmosphere grid point (-1 for upward-, +1 for downward-facing walls) 
+       INTEGER(iwp) ::  k_off    !< index offset between surface and atmosphere grid point (-1 for upward-, +1 for downward-facing walls) 
        INTEGER(iwp) ::  m       !< loop variable over all horizontal surf elements 
        
-       koff = surf%koff
+       k_off = surf%koff
        !$OMP PARALLEL DO PRIVATE( i, j, k )
        DO  m = 1, surf%ns 
@@ -1613,5 +1604,5 @@
           k   = surf%k(m)
 
-          surf%pt_surface(m) = pt(k+koff,j,i)
+          surf%pt_surface(m) = pt(k+k_off,j,i)
 
        ENDDO
@@ -1619,4 +1610,28 @@
     END SUBROUTINE calc_pt_surface
 
+!
+!-- Calculate virtual potential temperature at first grid level
+!-- ( only for upward-facing surfs )
+    SUBROUTINE calc_vpt_surface
+
+       IMPLICIT NONE
+
+       INTEGER(iwp) ::  k_off    !< index offset between surface and atmosphere grid point (-1 for upward-, +1 for downward-facing walls) 
+       INTEGER(iwp) ::  m       !< loop variable over all horizontal surf elements 
+       
+       k_off = surf%koff
+       !$OMP PARALLEL DO PRIVATE( i, j, k )
+       DO  m = 1, surf%ns 
+
+          i   = surf%i(m)            
+          j   = surf%j(m)
+          k   = surf%k(m)
+
+          surf%vpt_surface(m) = vpt(k+k_off,j,i)
+
+       ENDDO
+
+    END SUBROUTINE calc_vpt_surface
+    
 !
 !-- Calculate the other MOST scaling parameters theta*, q*, (qc*, qr*, nc*, nr*)
