Index: palm/trunk/SOURCE/land_surface_model_mod.f90
===================================================================
--- palm/trunk/SOURCE/land_surface_model_mod.f90	(revision 3145)
+++ palm/trunk/SOURCE/land_surface_model_mod.f90	(revision 3146)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! Adjustments for new surface structure
 ! 
 ! Former revisions:
@@ -2320,7 +2320,9 @@
 !
 !--       Update virtual potential temperature
-          vpt(k+k_off,j+j_off,i+i_off) = pt(k+k_off,j+j_off,i+i_off) *         &
+          surf%vpt_surface(m) = surf%pt_surface(m) *         &
                      ( 1.0_wp + 0.61_wp * q(k+k_off,j+j_off,i+i_off) )
 
+        
+                     
        ENDDO
 
Index: palm/trunk/SOURCE/surface_layer_fluxes_mod.f90
===================================================================
--- palm/trunk/SOURCE/surface_layer_fluxes_mod.f90	(revision 3145)
+++ 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*)
Index: palm/trunk/SOURCE/surface_mod.f90
===================================================================
--- palm/trunk/SOURCE/surface_mod.f90	(revision 3145)
+++ palm/trunk/SOURCE/surface_mod.f90	(revision 3146)
@@ -21,5 +21,5 @@
 ! Current revisions:
 ! ------------------
-! 
+! Added vpt_surface and vpt1
 ! 
 ! Former revisions:
@@ -251,6 +251,8 @@
        REAL(wp), DIMENSION(:), ALLOCATABLE ::  z0q       !< roughness length for humidity
 
-       REAL(wp), DIMENSION(:), ALLOCATABLE ::  pt1       !< Potential temperature at first grid level
+       REAL(wp), DIMENSION(:), ALLOCATABLE ::  pt1       !< potential temperature at first grid level
        REAL(wp), DIMENSION(:), ALLOCATABLE ::  qv1       !< mixing ratio at first grid level
+       REAL(wp), DIMENSION(:), ALLOCATABLE ::  vpt1      !< virtual potential temperature at first grid level
+       
        REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  css     !< scaling parameter chemical species
 !
@@ -302,4 +304,5 @@
 
        REAL(wp), DIMENSION(:), ALLOCATABLE   ::  pt_surface        !< skin-surface temperature
+       REAL(wp), DIMENSION(:), ALLOCATABLE   ::  vpt_surface       !< skin-surface virtual temperature
        REAL(wp), DIMENSION(:), ALLOCATABLE   ::  rad_net           !< net radiation 
        REAL(wp), DIMENSION(:), ALLOCATABLE   ::  rad_net_l         !< net radiation, used in USM
@@ -1047,4 +1050,6 @@
        DEALLOCATE ( surfaces%pt1 )
        DEALLOCATE ( surfaces%qv1 )
+       DEALLOCATE ( surfaces%vpt1 )
+       
 !
 !--       
@@ -1149,8 +1154,9 @@
        ALLOCATE ( surfaces%pt_surface(1:surfaces%ns) ) 
 !
-!--    Characteristic humidity and surface flux of latent heat
+!--    Characteristic humidity, surface flux of latent heat, and surface virtual potential temperature
        IF ( humidity )  THEN
           ALLOCATE ( surfaces%qs(1:surfaces%ns)   ) 
           ALLOCATE ( surfaces%qsws(1:surfaces%ns) )      
+          ALLOCATE ( surfaces%vpt_surface(1:surfaces%ns) ) 
        ENDIF 
 !
@@ -1171,4 +1177,5 @@
        ALLOCATE ( surfaces%pt1(1:surfaces%ns) )
        ALLOCATE ( surfaces%qv1(1:surfaces%ns) )
+       ALLOCATE ( surfaces%vpt1(1:surfaces%ns) )
 !
 !--       
@@ -1406,5 +1413,6 @@
        IF ( humidity )  THEN
           DEALLOCATE ( surfaces%qs ) 
-          DEALLOCATE ( surfaces%qsws )      
+          DEALLOCATE ( surfaces%qsws )  
+          DEALLOCATE ( surfaces%vpt_surface )
        ENDIF 
 !
@@ -1425,4 +1433,5 @@
        DEALLOCATE ( surfaces%pt1 )
        DEALLOCATE ( surfaces%qv1 )
+       DEALLOCATE ( surfaces%vpt1 )
 
        IF ( cloud_physics .AND. microphysics_seifert)  THEN
@@ -1527,6 +1536,7 @@
 !--    Characteristic humidity and surface flux of latent heat
        IF ( humidity )  THEN
-          ALLOCATE ( surfaces%qs(1:surfaces%ns)   ) 
-          ALLOCATE ( surfaces%qsws(1:surfaces%ns) )      
+          ALLOCATE ( surfaces%qs(1:surfaces%ns)          ) 
+          ALLOCATE ( surfaces%qsws(1:surfaces%ns)        )   
+          ALLOCATE ( surfaces%vpt_surface(1:surfaces%ns) )            
        ENDIF 
 !
@@ -1547,4 +1557,5 @@
        ALLOCATE ( surfaces%pt1(1:surfaces%ns) )
        ALLOCATE ( surfaces%qv1(1:surfaces%ns) )
+       ALLOCATE ( surfaces%vpt1(1:surfaces%ns) )
 
        IF ( cloud_physics .AND. microphysics_seifert)  THEN
@@ -2097,7 +2108,8 @@
                    surf%nrsws(num_h) = 0.0_wp
 
-                   surf%pt1(num_h) = 0.0_wp
-                   surf%qv1(num_h) = 0.0_wp
-
+                   surf%pt1(num_h)  = 0.0_wp
+                   surf%qv1(num_h)  = 0.0_wp
+                   surf%vpt1(num_h) = 0.0_wp
+                   
                 ENDIF
              ENDIF
@@ -2604,4 +2616,6 @@
                    IF ( ALLOCATED( surf_def_h(l)%pt_surface ) )                &
                       surf_h(l)%pt_surface(mm(l)) = surf_def_h(l)%pt_surface(m)
+                   IF ( ALLOCATED( surf_def_h(l)%vpt_surface ) )               &
+                      surf_h(l)%vpt_surface(mm(l)) = surf_def_h(l)%vpt_surface(m)                      
                    IF ( ALLOCATED( surf_def_h(l)%usws ) )                      &
                       surf_h(l)%usws(mm(l))    = surf_def_h(l)%usws(m)
@@ -2659,4 +2673,6 @@
                       IF ( ALLOCATED( surf_lsm_h%pt_surface ) )                &
                          surf_h(l)%pt_surface(mm(l)) = surf_lsm_h%pt_surface(m)
+                      IF ( ALLOCATED( surf_def_h(l)%vpt_surface ) )            &
+                      surf_h(l)%vpt_surface(mm(l)) = surf_def_h(l)%vpt_surface(m)
                       IF ( ALLOCATED( surf_lsm_h%usws ) )                      &
                          surf_h(0)%usws(mm(0))    = surf_lsm_h%usws(m)
@@ -2714,4 +2730,6 @@
                       IF ( ALLOCATED( surf_usm_h%pt_surface ) )                &
                          surf_h(l)%pt_surface(mm(l)) = surf_usm_h%pt_surface(m)
+                      IF ( ALLOCATED( surf_usm_h%vpt_surface ) )               &
+                         surf_h(l)%vpt_surface(mm(l)) = surf_usm_h%vpt_surface(m)
                       IF ( ALLOCATED( surf_usm_h%usws ) )                      &
                          surf_h(0)%usws(mm(0))    = surf_usm_h%usws(m)
@@ -2810,4 +2828,6 @@
                    IF ( ALLOCATED( surf_def_v(l)%pt_surface ) )                &
                       surf_v(l)%pt_surface(mm(l)) = surf_def_v(l)%pt_surface(m)
+                   IF ( ALLOCATED( surf_def_v(l)%vpt_surface ) )               &
+                      surf_v(l)%vpt_surface(mm(l)) = surf_def_v(l)%vpt_surface(m)
                    IF ( ALLOCATED( surf_def_v(l)%shf ) )                       &
                       surf_v(l)%shf(mm(l))     = surf_def_v(l)%shf(m)
@@ -2866,4 +2886,6 @@
                    IF ( ALLOCATED( surf_lsm_v(l)%pt_surface ) )                &
                       surf_v(l)%pt_surface(mm(l)) = surf_lsm_v(l)%pt_surface(m)
+                   IF ( ALLOCATED( surf_lsm_v(l)%vpt_surface ) )               &
+                      surf_v(l)%vpt_surface(mm(l)) = surf_lsm_v(l)%vpt_surface(m)
                    IF ( ALLOCATED( surf_lsm_v(l)%usws ) )                      &
                       surf_v(l)%usws(mm(l))    = surf_lsm_v(l)%usws(m)
@@ -2926,4 +2948,6 @@
                    IF ( ALLOCATED( surf_usm_v(l)%pt_surface ) )                &
                       surf_v(l)%pt_surface(mm(l)) = surf_usm_v(l)%pt_surface(m)
+                   IF ( ALLOCATED( surf_usm_v(l)%vpt_surface ) )               &
+                      surf_v(l)%vpt_surface(mm(l)) = surf_usm_v(l)%vpt_surface(m)
                    IF ( ALLOCATED( surf_usm_v(l)%usws ) )                      &
                       surf_v(l)%usws(mm(l))    = surf_usm_v(l)%usws(m)
@@ -3066,4 +3090,9 @@
           ENDIF
 
+          IF ( ALLOCATED ( surf_h(l)%vpt_surface ) )  THEN
+             CALL wrd_write_string( 'surf_h(' // dum // ')%vpt_surface' ) 
+             WRITE ( 14 )  surf_h(l)%vpt_surface
+          ENDIF
+
           IF ( ALLOCATED ( surf_h(l)%usws ) )  THEN
              CALL wrd_write_string( 'surf_h(' // dum // ')%usws' ) 
@@ -3191,4 +3220,9 @@
              CALL wrd_write_string( 'surf_v(' // dum // ')%pt_surface' )
              WRITE ( 14 )  surf_v(l)%pt_surface
+          ENDIF 
+
+          IF ( ALLOCATED ( surf_v(l)%vpt_surface ) )  THEN
+             CALL wrd_write_string( 'surf_v(' // dum // ')%vpt_surface' )
+             WRITE ( 14 )  surf_v(l)%vpt_surface
           ENDIF 
 
@@ -3430,4 +3464,7 @@
              IF ( ALLOCATED( surf_h(0)%pt_surface )  .AND.  kk == 1 )          &
                 READ ( 13 )  surf_h(0)%pt_surface
+          CASE ( 'surf_h(0)%vpt_surface' )         
+             IF ( ALLOCATED( surf_h(0)%vpt_surface )  .AND.  kk == 1 )         &
+                READ ( 13 )  surf_h(0)%vpt_surface
           CASE ( 'surf_h(0)%usws' )         
              IF ( ALLOCATED( surf_h(0)%usws )  .AND.  kk == 1 )                &
@@ -3507,4 +3544,7 @@
              IF ( ALLOCATED( surf_h(1)%pt_surface )  .AND.  kk == 1 )          &
                 READ ( 13 )  surf_h(1)%pt_surface
+          CASE ( 'surf_h(1)%vpt_surface' )         
+             IF ( ALLOCATED( surf_h(1)%vpt_surface )  .AND.  kk == 1 )          &
+                READ ( 13 )  surf_h(1)%vpt_surface
           CASE ( 'surf_h(1)%usws' )         
              IF ( ALLOCATED( surf_h(1)%usws )  .AND.  kk == 1 )                &
@@ -3584,4 +3624,7 @@
              IF ( ALLOCATED( surf_h(2)%pt_surface )  .AND.  kk == 1 )          &
                 READ ( 13 )  surf_h(2)%pt_surface
+          CASE ( 'surf_h(2)%vpt_surface' )         
+             IF ( ALLOCATED( surf_h(2)%vpt_surface )  .AND.  kk == 1 )          &
+                READ ( 13 )  surf_h(2)%vpt_surface
           CASE ( 'surf_h(2)%usws' )         
              IF ( ALLOCATED( surf_h(2)%usws )  .AND.  kk == 1 )                &
@@ -3663,4 +3706,7 @@
              IF ( ALLOCATED( surf_v(0)%pt_surface )  .AND.  kk == 1 )          &
                 READ ( 13 )  surf_v(0)%pt_surface
+          CASE ( 'surf_v(0)%vpt_surface' )         
+             IF ( ALLOCATED( surf_v(0)%vpt_surface )  .AND.  kk == 1 )          &
+                READ ( 13 )  surf_v(0)%vpt_surface
           CASE ( 'surf_v(0)%shf' )         
              IF ( ALLOCATED( surf_v(0)%shf )  .AND.  kk == 1 )                 &
@@ -3743,4 +3789,7 @@
              IF ( ALLOCATED( surf_v(1)%pt_surface )  .AND.  kk == 1 )          &
                 READ ( 13 )  surf_v(1)%pt_surface
+          CASE ( 'surf_v(1)%vpt_surface' )         
+             IF ( ALLOCATED( surf_v(1)%vpt_surface )  .AND.  kk == 1 )         &
+                READ ( 13 )  surf_v(1)%vpt_surface
           CASE ( 'surf_v(1)%shf' )         
              IF ( ALLOCATED( surf_v(1)%shf )  .AND.  kk == 1 )                 &
@@ -3823,4 +3872,7 @@
              IF ( ALLOCATED( surf_v(2)%pt_surface )  .AND.  kk == 1 )          &
                 READ ( 13 )  surf_v(2)%pt_surface
+          CASE ( 'surf_v(2)%vpt_surface' )         
+             IF ( ALLOCATED( surf_v(2)%vpt_surface )  .AND.  kk == 1 )          &
+                READ ( 13 )  surf_v(2)%vpt_surface
           CASE ( 'surf_v(2)%shf' )         
              IF ( ALLOCATED( surf_v(2)%shf )  .AND.  kk == 1 )                 &
@@ -3903,4 +3955,7 @@
              IF ( ALLOCATED( surf_v(3)%pt_surface )  .AND.  kk == 1 )          &
                 READ ( 13 )  surf_v(3)%pt_surface
+          CASE ( 'surf_v(3)%vpt_surface' )         
+             IF ( ALLOCATED( surf_v(3)%vpt_surface )  .AND.  kk == 1 )          &
+                READ ( 13 )  surf_v(3)%vpt_surface
           CASE ( 'surf_v(3)%shf' )         
              IF ( ALLOCATED( surf_v(3)%shf )  .AND.  kk == 1 )                 &
@@ -4097,4 +4152,10 @@
              ENDIF
 
+             IF ( INDEX( restart_string(1:length), '%vpt_surface' ) /= 0 )  THEN 
+                IF ( ALLOCATED( surf_target%vpt_surface )  .AND.               &
+                     ALLOCATED( surf_file%vpt_surface   ) )                    & 
+                   surf_target%vpt_surface(m_target) = surf_file%vpt_surface(m_file)
+             ENDIF
+             
              IF ( INDEX( restart_string(1:length), '%usws' ) /= 0 )  THEN 
                 IF ( ALLOCATED( surf_target%usws )  .AND.                      &
