Index: palm/trunk/SOURCE/diffusion_v.f90
===================================================================
--- palm/trunk/SOURCE/diffusion_v.f90	(revision 77)
+++ palm/trunk/SOURCE/diffusion_v.f90	(revision 102)
@@ -4,5 +4,5 @@
 ! Actual revisions:
 ! -----------------
-! 
+! Momentumflux at top (vswst) included as boundary condition
 !
 ! Former revisions:
@@ -52,5 +52,6 @@
 ! Call for all grid points
 !------------------------------------------------------------------------------!
-    SUBROUTINE diffusion_v( ddzu, ddzw, km, km_damp_x, tend, u, v, vsws, w )
+    SUBROUTINE diffusion_v( ddzu, ddzw, km, km_damp_x, tend, u, v, vsws, &
+                            vswst, w )
 
        USE control_parameters
@@ -64,5 +65,5 @@
        REAL    ::  ddzu(1:nzt+1), ddzw(1:nzt+1), km_damp_x(nxl-1:nxr+1)
        REAL    ::  tend(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1)
-       REAL, DIMENSION(:,:),   POINTER ::  vsws
+       REAL, DIMENSION(:,:),   POINTER ::  vsws, vswst
        REAL, DIMENSION(:,:,:), POINTER ::  km, u, v, w
        REAL, DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) ::  vsus
@@ -153,5 +154,5 @@
 !--          Compute vertical diffusion. In case of simulating a Prandtl
 !--          layer, index k starts at nzb_v_inner+2.
-             DO  k = nzb_diff_v(j,i), nzt
+             DO  k = nzb_diff_v(j,i), nzt_diff
 !
 !--             Interpolate eddy diffusivities on staggered gridpoints
@@ -193,6 +194,26 @@
                       & + ( kmzp * ( w(k,j,i) - w(k,j-1,i)     ) * ddy       &
                       &   ) * ddzw(k)                                        &
-                      & + ( kmzp * ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1)     &
+                      & + ( kmzp * ( v(k+1,j,i) - v(k,j,i)     ) * ddzu(k+1) &
                       &   + vsws(j,i)                                        &
+                      &   ) * ddzw(k)
+             ENDIF
+
+!
+!--          Vertical diffusion at the first gridpoint below the top boundary,
+!--          if the momentum flux at the top is prescribed by the user
+             IF ( use_top_fluxes )  THEN
+                k = nzt
+!
+!--             Interpolate eddy diffusivities on staggered gridpoints
+                kmzp = 0.25 * &
+                       ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
+                kmzm = 0.25 * &
+                       ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) )
+
+                tend(k,j,i) = tend(k,j,i)                                    &
+                      & - ( kmzm *  ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy      &
+                      &   ) * ddzw(k)                                        &
+                      & + ( -vswst(j,i)                                      &
+                      &   - kmzm * ( v(k,j,i)   - v(k-1,j,i)    ) * ddzu(k)  &
                       &   ) * ddzw(k)
              ENDIF
@@ -208,5 +229,5 @@
 !------------------------------------------------------------------------------!
     SUBROUTINE diffusion_v_ij( i, j, ddzu, ddzw, km, km_damp_x, tend, u, v, &
-                               vsws, w )
+                               vsws, vswst, w )
 
        USE control_parameters
@@ -221,5 +242,5 @@
        REAL    ::  tend(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1)
        REAL, DIMENSION(nzb:nzt+1)      ::  vsus
-       REAL, DIMENSION(:,:),   POINTER ::  vsws
+       REAL, DIMENSION(:,:),   POINTER ::  vsws, vswst
        REAL, DIMENSION(:,:,:), POINTER ::  km, u, v, w
 
@@ -302,5 +323,5 @@
 !--    Compute vertical diffusion. In case of simulating a Prandtl layer,
 !--    index k starts at nzb_v_inner+2.
-       DO  k = nzb_diff_v(j,i), nzt
+       DO  k = nzb_diff_v(j,i), nzt_diff
 !
 !--       Interpolate eddy diffusivities on staggered gridpoints
@@ -338,9 +359,29 @@
                       & + ( kmzp * ( w(k,j,i) - w(k,j-1,i)     ) * ddy       &
                       &   ) * ddzw(k)                                        &
-                      & + ( kmzp * ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1)     &
+                      & + ( kmzp * ( v(k+1,j,i) - v(k,j,i)     ) * ddzu(k+1) &
                       &   + vsws(j,i)                                        &
                       &   ) * ddzw(k)
        ENDIF
 
+!
+!--    Vertical diffusion at the first gridpoint below the top boundary,
+!--    if the momentum flux at the top is prescribed by the user
+       IF ( use_top_fluxes )  THEN
+          k = nzt
+!
+!--       Interpolate eddy diffusivities on staggered gridpoints
+          kmzp = 0.25 * &
+                 ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
+          kmzm = 0.25 * &
+                 ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) )
+
+          tend(k,j,i) = tend(k,j,i)                                          &
+                      & - ( kmzm *  ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy      &
+                      &   ) * ddzw(k)                                        &
+                      & + ( -vswst(j,i)                                      &
+                      &   - kmzm * ( v(k,j,i)   - v(k-1,j,i)    ) * ddzu(k)  &
+                      &   ) * ddzw(k)
+       ENDIF
+
     END SUBROUTINE diffusion_v_ij
 
