Index: palm/trunk/SOURCE/diffusion_s.f90
===================================================================
--- palm/trunk/SOURCE/diffusion_s.f90	(revision 4)
+++ palm/trunk/SOURCE/diffusion_s.f90	(revision 19)
@@ -4,5 +4,6 @@
 ! Actual revisions:
 ! -----------------
-! 
+! Calculation extended for gridpoint nzt, fluxes can be given at top,
+! +s_flux_t in parameter list, s_flux renamed s_flux_b
 !
 ! Former revisions:
@@ -39,5 +40,5 @@
 ! Call for all grid points
 !------------------------------------------------------------------------------!
-    SUBROUTINE diffusion_s( ddzu, ddzw, kh, s, s_flux, tend )
+    SUBROUTINE diffusion_s( ddzu, ddzw, kh, s, s_flux_b, s_flux_t, tend )
 
        USE control_parameters
@@ -51,5 +52,5 @@
        REAL    ::  ddzu(1:nzt+1), ddzw(1:nzt)
        REAL    ::  tend(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1)
-       REAL, DIMENSION(:,:),   POINTER ::  s_flux
+       REAL, DIMENSION(:,:),   POINTER ::  s_flux_b, s_flux_t
        REAL, DIMENSION(:,:,:), POINTER ::  kh, s
 
@@ -58,5 +59,5 @@
 !
 !--          Compute horizontal diffusion
-             DO  k = nzb_s_outer(j,i)+1, nzt-1
+             DO  k = nzb_s_outer(j,i)+1, nzt
 
                 tend(k,j,i) = tend(k,j,i)                                     &
@@ -97,6 +98,7 @@
 !
 !--          Compute vertical diffusion. In case that surface fluxes have been
-!--          presribed or computed, index k starts at nzb+2.
-             DO  k = nzb_diff_s_inner(j,i), nzt-1
+!--          prescribed or computed at bottom and/or top, index k starts/ends at
+!--          nzb+2 or nzt-1, respectively.
+             DO  k = nzb_diff_s_inner(j,i), nzt_diff
 
                 tend(k,j,i) = tend(k,j,i)                                     &
@@ -108,5 +110,5 @@
 
 !
-!--          Vertical diffusion at the first computational gridpoint in &
+!--          Vertical diffusion at the first computational gridpoint along
 !--          z-direction
              IF ( use_surface_fluxes )  THEN
@@ -118,5 +120,21 @@
                                                * ( s(k+1,j,i)-s(k,j,i) )      &
                                                * ddzu(k+1)                    &
-                                           + s_flux(j,i)                      &
+                                           + s_flux_b(j,i)                    &
+                                         ) * ddzw(k)
+
+             ENDIF
+
+!
+!--          Vertical diffusion at the last computational gridpoint along
+!--          z-direction
+             IF ( use_top_fluxes )  THEN
+
+                k = nzt
+
+                tend(k,j,i) = tend(k,j,i)                                     &
+                                       + ( - s_flux_t(j,i)                    &
+                                         - 0.5 * ( kh(k-1,j,i)+kh(k,j,i) )    &
+                                               * ( s(k,j,i)-s(k-1,j,i) )      &
+                                               * ddzu(k)                      &
                                          ) * ddzw(k)
 
@@ -132,5 +150,6 @@
 ! Call for grid point i,j
 !------------------------------------------------------------------------------!
-    SUBROUTINE diffusion_s_ij( i, j, ddzu, ddzw, kh, s, s_flux, tend )
+    SUBROUTINE diffusion_s_ij( i, j, ddzu, ddzw, kh, s, s_flux_b, s_flux_t, &
+                               tend )
 
        USE control_parameters
@@ -144,10 +163,10 @@
        REAL    ::  ddzu(1:nzt+1), ddzw(1:nzt)
        REAL    ::  tend(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1)
-       REAL, DIMENSION(:,:),   POINTER ::  s_flux
+       REAL, DIMENSION(:,:),   POINTER ::  s_flux_b, s_flux_t
        REAL, DIMENSION(:,:,:), POINTER ::  kh, s
 
 !
 !--    Compute horizontal diffusion
-       DO  k = nzb_s_outer(j,i)+1, nzt-1
+       DO  k = nzb_s_outer(j,i)+1, nzt
 
           tend(k,j,i) = tend(k,j,i)                                           &
@@ -188,6 +207,7 @@
 !
 !--    Compute vertical diffusion. In case that surface fluxes have been
-!--    presribed or computed, index k starts at nzb+2.
-       DO  k = nzb_diff_s_inner(j,i), nzt-1
+!--    prescribed or computed at bottom and/or top, index k starts/ends at
+!--    nzb+2 or nzt-1, respectively.
+       DO  k = nzb_diff_s_inner(j,i), nzt_diff
 
           tend(k,j,i) = tend(k,j,i)                                           &
@@ -199,18 +219,31 @@
 
 !
-!--    Vertical diffusion at the first computational gridpoint in z-direction
+!--    Vertical diffusion at the first computational gridpoint along z-direction
        IF ( use_surface_fluxes )  THEN
 
           k = nzb_s_inner(j,i)+1
 
-          tend(k,j,i) = tend(k,j,i)                                           &
-                                       + ( 0.5 * ( kh(k,j,i)+kh(k+1,j,i) )    &
-                                               * ( s(k+1,j,i)-s(k,j,i) )      &
-                                               * ddzu(k+1)                    &
-                                           + s_flux(j,i)                      &
-                                         ) * ddzw(k)
+          tend(k,j,i) = tend(k,j,i) + ( 0.5 * ( kh(k,j,i)+kh(k+1,j,i) )  &
+                                            * ( s(k+1,j,i)-s(k,j,i) )    &
+                                            * ddzu(k+1)                  &
+                                        + s_flux_b(j,i)                  &
+                                      ) * ddzw(k)
 
        ENDIF
 
+!
+!--    Vertical diffusion at the last computational gridpoint along z-direction
+       IF ( use_top_fluxes )  THEN
+
+          k = nzt
+
+          tend(k,j,i) = tend(k,j,i) + ( - s_flux_t(j,i)                  &
+                                      - 0.5 * ( kh(k-1,j,i)+kh(k,j,i) )  &
+                                            * ( s(k,j,i)-s(k-1,j,i) )    &
+                                            * ddzu(k)                    &
+                                      ) * ddzw(k)
+
+       ENDIF
+
     END SUBROUTINE diffusion_s_ij
 
