Index: palm/trunk/SOURCE/advec_s_pw.f90
===================================================================
--- palm/trunk/SOURCE/advec_s_pw.f90	(revision 18)
+++ palm/trunk/SOURCE/advec_s_pw.f90	(revision 19)
@@ -4,5 +4,5 @@
 ! Actual revisions:
 ! -----------------
-! 
+! Calculation extended for gridpoint nzt
 !
 ! Former revisions:
@@ -58,5 +58,5 @@
        DO  i = nxl, nxr
           DO  j = nys, nyn
-             DO  k = nzb_s_inner(j,i)+1, nzt-1
+             DO  k = nzb_s_inner(j,i)+1, nzt
                 tend(k,j,i) = tend(k,j,i)                                      &
               -0.5 * ( ( u(k,j,i+1) - u_gtrans ) * ( sk(k,j,i+1) - sk(k,j,i) ) &
@@ -93,5 +93,5 @@
 
 
-       DO  k = nzb_s_inner(j,i)+1, nzt-1
+       DO  k = nzb_s_inner(j,i)+1, nzt
           tend(k,j,i) = tend(k,j,i)                                            &
               -0.5 * ( ( u(k,j,i+1) - u_gtrans ) * ( sk(k,j,i+1) - sk(k,j,i) ) &
Index: palm/trunk/SOURCE/boundary_conds.f90
===================================================================
--- palm/trunk/SOURCE/boundary_conds.f90	(revision 18)
+++ palm/trunk/SOURCE/boundary_conds.f90	(revision 19)
@@ -4,5 +4,7 @@
 ! Actual revisions:
 ! -----------------
-! 
+! Boundary conditions for e(nzt), pt(nzt), and q(nzt) removed because these
+! values are now calculated by the prognostic equation,
+! Dirichlet and zero gradient condition for pt established at top boundary
 !
 ! Former revisions:
@@ -100,6 +102,13 @@
 !
 !--    Temperature at top boundary
-       IF ( ibc_pt_t == 1 )  THEN
-          pt(nzt,:,:)   = pt(nzt-1,:,:) + bc_pt_t_val * dzu(nzt)
+       IF ( ibc_pt_t == 0 )  THEN
+          IF ( timestep_scheme(1:5) /= 'runge' )  THEN
+             pt(nzt+1,:,:) = pt_m(nzt+1,:,:)
+          ELSE
+             pt(nzt+1,:,:) = pt_p(nzt+1,:,:)  ! pt_m not used for Runge-Kutta
+          ENDIF
+       ELSEIF ( ibc_pt_t == 1 )  THEN
+          pt(nzt+1,:,:) = pt(nzt,:,:)
+       ELSEIF ( ibc_pt_t == 2 )  THEN
           pt(nzt+1,:,:) = pt(nzt,:,:)   + bc_pt_t_val * dzu(nzt+1)
        ENDIF
@@ -114,5 +123,4 @@
              ENDDO
           ENDDO
-          e(nzt,:,:)   = e(nzt-1,:,:)
           e(nzt+1,:,:) = e(nzt,:,:)
        ENDIF
@@ -147,5 +155,4 @@
 !
 !--       Top boundary
-          q(nzt,:,:)   = q(nzt-1,:,:) + bc_q_t_val * dzu(nzt)
           q(nzt+1,:,:) = q(nzt,:,:)   + bc_q_t_val * dzu(nzt+1)
        ENDIF
Index: palm/trunk/SOURCE/calc_liquid_water_content.f90
===================================================================
--- palm/trunk/SOURCE/calc_liquid_water_content.f90	(revision 18)
+++ palm/trunk/SOURCE/calc_liquid_water_content.f90	(revision 19)
@@ -4,5 +4,5 @@
 ! Actual revisions:
 ! -----------------
-! 
+! Old comment removed
 !
 ! Former revisions:
@@ -79,7 +79,3 @@
     ENDDO
     
-!
-!-- Setting boundary values of ql
-!    CALL exchange_horiz( ql, 0, 0 )    
-
  END SUBROUTINE calc_liquid_water_content
Index: palm/trunk/SOURCE/calc_precipitation.f90
===================================================================
--- palm/trunk/SOURCE/calc_precipitation.f90	(revision 18)
+++ palm/trunk/SOURCE/calc_precipitation.f90	(revision 19)
@@ -4,5 +4,5 @@
 ! Actual revisions:
 ! -----------------
-! 
+! Calculation extended for gridpoint nzt
 !
 ! Former revisions:
@@ -54,5 +54,5 @@
        DO  i = nxl, nxr
           DO  j = nys, nyn
-             DO  k = nzb_2d(j,i)+1, nzt-1    
+             DO  k = nzb_2d(j,i)+1, nzt
 
                 IF ( ql(k,j,i) > ql_crit )  THEN
@@ -87,5 +87,5 @@
 
 
-       DO  k = nzb_2d(j,i)+1, nzt-1    
+       DO  k = nzb_2d(j,i)+1, nzt
 
           IF ( ql(k,j,i) > ql_crit )  THEN
Index: palm/trunk/SOURCE/check_parameters.f90
===================================================================
--- palm/trunk/SOURCE/check_parameters.f90	(revision 18)
+++ palm/trunk/SOURCE/check_parameters.f90	(revision 19)
@@ -4,5 +4,6 @@
 ! Actual revisions:
 ! -----------------
-! 
+! Temperature and humidity gradients at top are now calculated for nzt+1,
+! top_heatflux and respective boundary condition bc_pt_t is checked
 !
 ! Former revisions:
@@ -483,5 +484,5 @@
 !--    Store temperature gradient at the top boundary for possile Neumann
 !--    boundary condition
-       bc_pt_t_val = ( pt_init(nzt) - pt_init(nzt-1) ) / dzu(nzt)
+       bc_pt_t_val = ( pt_init(nzt+1) - pt_init(nzt) ) / dzu(nzt+1)
 
 !
@@ -539,5 +540,5 @@
 !--       Store humidity gradient at the top boundary for possile Neumann
 !--       boundary condition
-          bc_q_t_val = ( q_init(nzt) - q_init(nzt-1) ) / dzu(nzt)
+          bc_q_t_val = ( q_init(nzt+1) - q_init(nzt) ) / dzu(nzt+1)
 
        ENDIF
@@ -795,4 +796,6 @@
     ELSEIF ( bc_pt_t == 'neumann' )  THEN
        ibc_pt_t = 1
+    ELSEIF ( bc_pt_t == 'initial_gradient' )  THEN
+       ibc_pt_t = 2
     ELSE
        IF ( myid == 0 )  THEN
@@ -804,4 +807,9 @@
 
     IF ( surface_heatflux == 9999999.9 )  constant_heatflux = .FALSE.
+    IF ( top_heatflux     == 9999999.9 )  THEN
+       constant_top_heatflux = .FALSE.
+    ELSE
+       use_top_fluxes = .TRUE.    ! because this is currently the only choice
+    ENDIF
 
 !
@@ -823,4 +831,18 @@
           PRINT*, '    allowed with pt_surface_initial_change (/=0) = ', &
                   pt_surface_initial_change
+       ENDIF
+       CALL local_stop
+    ENDIF
+
+!
+!-- A given temperature at the top implies Dirichlet boundary condition for
+!-- temperature. In this case specification of a constant heat flux is
+!-- forbidden.
+    IF ( ibc_pt_t == 0  .AND.   constant_top_heatflux  .AND. &
+         top_heatflux /= 0.0 )  THEN
+       IF ( myid == 0 )  THEN
+          PRINT*, '+++ check_parameters:'
+          PRINT*, '    boundary_condition: bc_pt_t = ', bc_pt_t
+          PRINT*, '    is not allowed with constant_top_heatflux = .TRUE.'
        ENDIF
        CALL local_stop
Index: palm/trunk/SOURCE/diffusion_e.f90
===================================================================
--- palm/trunk/SOURCE/diffusion_e.f90	(revision 18)
+++ palm/trunk/SOURCE/diffusion_e.f90	(revision 19)
@@ -4,5 +4,5 @@
 ! Actual revisions:
 ! -----------------
-! 
+! Calculation extended for gridpoint nzt
 !
 ! Former revisions:
@@ -56,5 +56,5 @@
        REAL, DIMENSION(:,:), POINTER   ::  rif
        REAL, DIMENSION(:,:,:), POINTER ::  e, km, theta
-       REAL, DIMENSION(nzb+1:nzt-1,nys:nyn) ::  dissipation, l, ll
+       REAL, DIMENSION(nzb+1:nzt,nys:nyn) ::  dissipation, l, ll
  
 
@@ -72,5 +72,5 @@
              ENDIF
 
-             DO  k = nzb_s_inner(j,i)+1, nzt-1
+             DO  k = nzb_s_inner(j,i)+1, nzt
 !
 !--             Calculate the mixing length (for dissipation)
@@ -102,5 +102,5 @@
 !--       Calculate the tendency terms
           DO  j = nys, nyn
-             DO  k = nzb_s_inner(j,i)+1, nzt-1
+             DO  k = nzb_s_inner(j,i)+1, nzt
 
                  dissipation(k,j) = ( 0.19 + 0.74 * l(k,j) / ll(k,j) ) * &
@@ -130,5 +130,5 @@
           IF ( use_sgs_for_particles )  THEN
              DO  j = nys, nyn
-                DO  k = nzb_s_inner(j,i)+1, nzt-1
+                DO  k = nzb_s_inner(j,i)+1, nzt
                    diss(k,j,i) = dissipation(k,j)
                 ENDDO
@@ -171,5 +171,5 @@
        REAL, DIMENSION(:,:), POINTER   ::  rif
        REAL, DIMENSION(:,:,:), POINTER ::  e, km, theta
-       REAL, DIMENSION(nzb+1:nzt-1)    ::  dissipation, l, ll
+       REAL, DIMENSION(nzb+1:nzt)    ::  dissipation, l, ll
 
 
@@ -187,5 +187,5 @@
 !
 !--    Calculate the mixing length (for dissipation)
-       DO  k = nzb_s_inner(j,i)+1, nzt-1
+       DO  k = nzb_s_inner(j,i)+1, nzt
           dpt_dz = ( theta(k+1,j,i) - theta(k-1,j,i) ) * dd2zu(k)
           IF ( dpt_dz > 0.0 ) THEN
@@ -234,5 +234,5 @@
 !--    Store dissipation if needed for calculating the sgs particle velocities
        IF ( use_sgs_for_particles )  THEN
-          DO  k = nzb_s_inner(j,i)+1, nzt-1
+          DO  k = nzb_s_inner(j,i)+1, nzt
              diss(k,j,i) = dissipation(k)
           ENDDO
Index: palm/trunk/SOURCE/diffusion_s.f90
===================================================================
--- palm/trunk/SOURCE/diffusion_s.f90	(revision 18)
+++ 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
 
Index: palm/trunk/SOURCE/flow_statistics.f90
===================================================================
--- palm/trunk/SOURCE/flow_statistics.f90	(revision 18)
+++ palm/trunk/SOURCE/flow_statistics.f90	(revision 19)
@@ -4,5 +4,5 @@
 ! Actual revisions:
 ! -----------------
-! 
+! fluxes at top modified (tswst, qswst)
 !
 ! Former revisions:
@@ -347,4 +347,7 @@
                                  + ( w(k,j,i)   - w(k,j-1,i) ) * ddy           &
                                                                ) * rmask(j,i,sr)
+             ENDDO
+
+             DO  k = nzb_diff_s_outer(j,i)-1, nzt_diff
 !
 !--             Heat flux w"pt"
@@ -417,4 +420,34 @@
                    sums_l(nzb,48,tn) = sums_l(nzb,48,tn) + &
                                        qsws(j,i) * rmask(j,i,sr)  ! w"q" (w"qv")
+                ENDIF
+             ENDIF
+
+!
+!--          Subgridscale fluxes at the top surface
+             IF ( use_top_fluxes )  THEN
+                sums_l(nzt,16,tn) = sums_l(nzt,16,tn) + &
+                                    tswst(j,i)  * rmask(j,i,sr)   ! w"pt"
+                sums_l(nzt,58,tn) = sums_l(nzt,58,tn) + &
+                                    0.0 * rmask(j,i,sr)           ! u"pt"
+                sums_l(nzt,61,tn) = sums_l(nzt,61,tn) + &
+                                    0.0 * rmask(j,i,sr)           ! v"pt"
+                IF ( moisture )  THEN
+                   sums_l(nzt,48,tn) = sums_l(nzt,48,tn) + &
+                                       qswst(j,i) * rmask(j,i,sr)  ! w"q" (w"qv")
+                   IF ( cloud_physics )  THEN
+                      sums_l(nzt,45,tn) = sums_l(nzt,45,tn) + (           &
+                                          ( 1.0 + 0.61 * q(nzt,j,i) ) *   &
+                                          tswst(j,i) + 0.61 * pt(nzt,j,i) * &
+                                                     qsws(j,i)            &
+                                                              )
+!
+!--                   Formula does not work if ql(nzb) /= 0.0
+                      sums_l(nzt,51,tn) = sums_l(nzt,51,tn) + &   ! w"q" (w"qv")
+                                          qswst(j,i) * rmask(j,i,sr)
+                   ENDIF
+                ENDIF
+                IF ( passive_scalar )  THEN
+                   sums_l(nzt,48,tn) = sums_l(nzt,48,tn) + &
+                                       qswst(j,i) * rmask(j,i,sr)  ! w"q" (w"qv")
                 ENDIF
              ENDIF
Index: palm/trunk/SOURCE/header.f90
===================================================================
--- palm/trunk/SOURCE/header.f90	(revision 18)
+++ palm/trunk/SOURCE/header.f90	(revision 19)
@@ -345,7 +345,9 @@
     ENDIF
     IF ( ibc_pt_t == 0 )  THEN
-       roben  = TRIM( roben  ) // ' pt(nzt) = pt_top'
-    ELSE
-       roben  = TRIM( roben  ) // ' pt(nzt) = pt(nzt-1) + dpt/dz'
+       roben  = TRIM( roben  ) // ' pt(nzt+1) = pt_top'
+    ELSEIF( ibc_pt_t == 1 )  THEN
+       roben  = TRIM( roben  ) // ' pt(nzt+1) = pt(nzt)'
+    ELSEIF( ibc_pt_t == 2 )  THEN
+       roben  = TRIM( roben  ) // ' pt(nzt+1) = pt(nzt) + dpt/dz_ini'
     ENDIF
 
@@ -404,4 +406,14 @@
        IF ( passive_scalar  .AND.  constant_waterflux )  THEN
           WRITE ( io, 313 ) surface_waterflux
+       ENDIF
+    ENDIF
+
+    IF ( use_top_fluxes )  THEN
+       WRITE ( io, 304 )
+       IF ( constant_top_heatflux )  THEN
+          WRITE ( io, 306 )  top_heatflux
+       ENDIF
+       IF ( moisture  .OR.  passive_scalar )  THEN
+          WRITE ( io, 315 )
        ENDIF
     ENDIF
@@ -1220,7 +1232,8 @@
              ' B. bound.: ',A/ &
              ' T. bound.: ',A)
-303 FORMAT (/' Surface fluxes are used in diffusion terms at k=1')
-305 FORMAT (//'    Prandtl-Layer between surface and first computational ', &
-               'u,v-level:'// &
+303 FORMAT (/' Bottom surface fluxes are used in diffusion terms at k=1')
+304 FORMAT (/' Top surface fluxes are used in diffusion terms at k=nzt')
+305 FORMAT (//'    Prandtl-Layer between bottom surface and first ', &
+               'computational u,v-level:'// &
              '       zp = ',F6.2,' m   z0 = ',F6.4,' m   kappa = ',F4.2/ &
              '       Rif value range:   ',F6.2,' <= rif <=',F6.2)
@@ -1234,4 +1247,5 @@
 313 FORMAT ('       Predefined constant scalar flux: ',E10.3,' kg/(m**2 s)')
 314 FORMAT ('       Predefined scalar value at the surface')
+315 FORMAT ('       Humidity / scalar flux at top surface is 0.0')
 317 FORMAT (//' Lateral boundaries:'/ &
             '       left/right:  ',A/    &
Index: palm/trunk/SOURCE/impact_of_latent_heat.f90
===================================================================
--- palm/trunk/SOURCE/impact_of_latent_heat.f90	(revision 18)
+++ palm/trunk/SOURCE/impact_of_latent_heat.f90	(revision 19)
@@ -4,5 +4,5 @@
 ! Actual revisions:
 ! -----------------
-! 
+! Calculation extended for gridpoint nzt
 !
 ! Former revisions:
@@ -53,5 +53,5 @@
        DO  i = nxl, nxr
           DO  j = nys, nyn
-             DO  k = nzb_2d(j,i)+1, nzt-1    
+             DO  k = nzb_2d(j,i)+1, nzt
 
                 IF ( ql(k,j,i) > ql_crit )  THEN
@@ -87,5 +87,5 @@
 
 
-       DO  k = nzb_2d(j,i)+1, nzt-1    
+       DO  k = nzb_2d(j,i)+1, nzt
 
           IF ( ql(k,j,i) > ql_crit )  THEN
Index: palm/trunk/SOURCE/init_3d_model.f90
===================================================================
--- palm/trunk/SOURCE/init_3d_model.f90	(revision 18)
+++ palm/trunk/SOURCE/init_3d_model.f90	(revision 19)
@@ -7,5 +7,5 @@
 ! Actual revisions:
 ! -----------------
-! 
+! +handling of top fluxes
 !
 ! Former revisions:
@@ -77,15 +77,16 @@
     ALLOCATE( km_damp_x(nxl-1:nxr+1), km_damp_y(nys-1:nyn+1) )
 
-    ALLOCATE( rif_1(nys-1:nyn+1,nxl-1:nxr+1), shf_1(nys-1:nyn+1,nxl-1:nxr+1),  &
-              ts(nys-1:nyn+1,nxl-1:nxr+1), us(nys-1:nyn+1,nxl-1:nxr+1),        &
-              usws_1(nys-1:nyn+1,nxl-1:nxr+1), vsws_1(nys-1:nyn+1,nxl-1:nxr+1),&
-              z0(nys-1:nyn+1,nxl-1:nxr+1) )
+    ALLOCATE( rif_1(nys-1:nyn+1,nxl-1:nxr+1), shf_1(nys-1:nyn+1,nxl-1:nxr+1), &
+              ts(nys-1:nyn+1,nxl-1:nxr+1), tswst_1(nys-1:nyn+1,nxl-1:nxr+1),  &
+              us(nys-1:nyn+1,nxl-1:nxr+1), usws_1(nys-1:nyn+1,nxl-1:nxr+1),   &
+              vsws_1(nys-1:nyn+1,nxl-1:nxr+1), z0(nys-1:nyn+1,nxl-1:nxr+1) )
 
     IF ( timestep_scheme(1:5) /= 'runge' )  THEN
 !
 !--    Leapfrog scheme needs two timelevels of diffusion quantities
-       ALLOCATE( rif_2(nys-1:nyn+1,nxl-1:nxr+1),  &
-                 shf_2(nys-1:nyn+1,nxl-1:nxr+1),  &
-                 usws_2(nys-1:nyn+1,nxl-1:nxr+1), &
+       ALLOCATE( rif_2(nys-1:nyn+1,nxl-1:nxr+1),   &
+                 shf_2(nys-1:nyn+1,nxl-1:nxr+1),   &
+                 tswst_2(nys-1:nyn+1,nxl-1:nxr+1), &
+                 usws_2(nys-1:nyn+1,nxl-1:nxr+1),  &
                  vsws_2(nys-1:nyn+1,nxl-1:nxr+1) )
     ENDIF
@@ -121,8 +122,10 @@
 !--    2D-moisture/scalar arrays
        ALLOCATE ( qs(nys-1:nyn+1,nxl-1:nxr+1),     &
-                  qsws_1(nys-1:nyn+1,nxl-1:nxr+1) )
+                  qsws_1(nys-1:nyn+1,nxl-1:nxr+1), &
+                  qswst_1(nys-1:nyn+1,nxl-1:nxr+1) )
 
        IF ( timestep_scheme(1:5) /= 'runge' )  THEN
-          ALLOCATE( qsws_2(nys-1:nyn+1,nxl-1:nxr+1) )
+          ALLOCATE( qsws_2(nys-1:nyn+1,nxl-1:nxr+1), &
+                    qswst_2(nys-1:nyn+1,nxl-1:nxr+1) )
        ENDIF
 !
@@ -177,8 +180,9 @@
     IF ( timestep_scheme(1:5) /= 'runge' )  THEN
 
-       rif_m  => rif_1;   rif  => rif_2
-       shf_m  => shf_1;   shf  => shf_2
-       usws_m => usws_1;  usws => usws_2
-       vsws_m => vsws_1;  vsws => vsws_2
+       rif_m   => rif_1;    rif   => rif_2
+       shf_m   => shf_1;    shf   => shf_2
+       tswst_m => tswst_1;  tswst => tswst_2
+       usws_m  => usws_1;   usws  => usws_2
+       vsws_m  => vsws_1;   vsws  => vsws_2
        e_m  => e_1;   e  => e_2;   e_p  => e_3;   te_m  => e_3
        kh_m => kh_1;  kh => kh_2
@@ -190,5 +194,6 @@
 
        IF ( moisture  .OR.  passive_scalar )  THEN
-          qsws_m => qsws_1;  qsws => qsws_2
+          qsws_m  => qsws_1;   qsws  => qsws_2
+          qswst_m => qswst_1;  qswst => qswst_2
           q_m    => q_1;     q    => q_2;     q_p => q_3;     tq_m => q_3
           IF ( moisture )        vpt_m  => vpt_1;   vpt  => vpt_2
@@ -202,18 +207,20 @@
     ELSE
 
-       rif  => rif_1
-       shf  => shf_1
-       usws => usws_1
-       vsws => vsws_1
-       e    => e_1;   e_p  => e_2;   te_m  => e_3;   e_m  => e_3
-       kh   => kh_1
-       km   => km_1
-       pt   => pt_1;  pt_p => pt_2;  tpt_m => pt_3;  pt_m => pt_3
-       u    => u_1;   u_p  => u_2;   tu_m  => u_3;   u_m  => u_3
-       v    => v_1;   v_p  => v_2;   tv_m  => v_3;   v_m  => v_3
-       w    => w_1;   w_p  => w_2;   tw_m  => w_3;   w_m  => w_3
+       rif   => rif_1
+       shf   => shf_1
+       tswst => tswst_1
+       usws  => usws_1
+       vsws  => vsws_1
+       e     => e_1;   e_p  => e_2;   te_m  => e_3;   e_m  => e_3
+       kh    => kh_1
+       km    => km_1
+       pt    => pt_1;  pt_p => pt_2;  tpt_m => pt_3;  pt_m => pt_3
+       u     => u_1;   u_p  => u_2;   tu_m  => u_3;   u_m  => u_3
+       v     => v_1;   v_p  => v_2;   tv_m  => v_3;   v_m  => v_3
+       w     => w_1;   w_p  => w_2;   tw_m  => w_3;   w_m  => w_3
 
        IF ( moisture  .OR.  passive_scalar )  THEN
           qsws   => qsws_1
+          qswst  => qswst_1
           q      => q_1;     q_p  => q_2;     tq_m => q_3;    q_m => q_3
           IF ( moisture )        vpt  => vpt_1
@@ -453,5 +460,5 @@
 
 !
-!--    Initialize surface fluxes
+!--    Initialize fluxes at bottom surface
        IF ( use_surface_fluxes )  THEN
 
@@ -486,4 +493,24 @@
              ENDIF
           ENDIF
+
+       ENDIF
+
+!
+!--    Initialize fluxes at top surface
+!--    Currently, only the heatflux can be prescribed. The latent flux is
+!--    zeri in this case!
+       IF ( use_top_fluxes )  THEN
+
+          IF ( constant_top_heatflux )  THEN
+!
+!--          Heat flux is prescribed
+             tswst = top_heatflux
+             IF ( ASSOCIATED( tswst_m ) )  tswst_m = tswst
+
+             IF ( moisture  .OR.  passive_scalar )  THEN
+                qswst = 0.0
+                IF ( ASSOCIATED( qswst_m ) )  qswst_m = qswst
+             ENDIF
+         ENDIF
 
        ENDIF
Index: palm/trunk/SOURCE/init_grid.f90
===================================================================
--- palm/trunk/SOURCE/init_grid.f90	(revision 18)
+++ palm/trunk/SOURCE/init_grid.f90	(revision 19)
@@ -4,5 +4,5 @@
 ! Actual revisions:
 ! -----------------
-! 
+! Setting of nzt_diff
 !
 ! Former revisions:
@@ -202,5 +202,5 @@
 
 !
-!-- Define vertical gridpoint from which on the usual finite difference
+!-- Define vertical gridpoint from (or to) which on the usual finite difference
 !-- form (which does not use surface fluxes) is applied 
     IF ( prandtl_layer  .OR.  use_surface_fluxes )  THEN
@@ -208,4 +208,9 @@
     ELSE
        nzb_diff = nzb + 1
+    ENDIF
+    IF ( use_top_fluxes )  THEN
+       nzt_diff = nzt - 1
+    ELSE
+       nzt_diff = nzt
     ENDIF
 
Index: palm/trunk/SOURCE/init_pt_anomaly.f90
===================================================================
--- palm/trunk/SOURCE/init_pt_anomaly.f90	(revision 18)
+++ palm/trunk/SOURCE/init_pt_anomaly.f90	(revision 19)
@@ -4,5 +4,5 @@
 ! Actual revisions:
 ! -----------------
-! 
+! Calculation extended for gridpoint nzt
 !
 ! Former revisions:
@@ -47,5 +47,5 @@
     DO  i = nxl, nxr
        DO  j = nys, nyn
-          DO  k = nzb+1, nzt-1
+          DO  k = nzb+1, nzt
              x = ( i - ic ) * dx
              y = ( j - jc ) * dy
Index: palm/trunk/SOURCE/modules.f90
===================================================================
--- palm/trunk/SOURCE/modules.f90	(revision 18)
+++ palm/trunk/SOURCE/modules.f90	(revision 19)
@@ -5,5 +5,7 @@
 ! Actual revisions:
 ! -----------------
-! 
+! +constant_top_heatflux, top_heatflux, use_top_fluxes, +arrays for top fluxes,
+! +nzt_diff, default of bc_pt_t renamed "initial_gradient"
+! Bugfix: p is not a pointer
 !
 ! Former revisions:
@@ -75,9 +77,10 @@
 
     REAL, DIMENSION(:,:), ALLOCATABLE, TARGET ::                               &
-          qsws_1, qsws_2, rif_1, rif_2, shf_1, shf_2, usws_1, usws_2,          &
-          vsws_1, vsws_2
+          qsws_1, qsws_2, qswst_1, qswst_2, rif_1, rif_2, shf_1, shf_2,        &
+          tswst_1, tswst_2, usws_1, usws_2, vsws_1, vsws_2
 
     REAL, DIMENSION(:,:), POINTER ::                                           &
-          qsws, qsws_m, rif, rif_m, shf, shf_m, usws, usws_m, vsws, vsws_m
+          qsws, qsws_m, qswst, qswst_m, rif, rif_m, shf, shf_m, tswst,         &
+          tswst_m, usws, usws_m, vsws, vsws_m
 
     REAL, DIMENSION(:,:,:), ALLOCATABLE ::                                     &
@@ -88,10 +91,10 @@
 
     REAL, DIMENSION(:,:,:), ALLOCATABLE, TARGET ::                             &
-          e_1, e_2, e_3, kh_1, kh_2, km_1, km_2, pt_1, pt_2, pt_3, q_1, q_2,   &
-          q_3, ql_1, ql_2, u_1, u_2, u_3, v_1, v_2, v_3, vpt_1, vpt_2, w_1,    &
-          w_2, w_3
+          e_1, e_2, e_3, kh_1, kh_2, km_1, km_2, p, pt_1, pt_2, pt_3, q_1,     &
+          q_2, q_3, ql_1, ql_2, u_1, u_2, u_3, v_1, v_2, v_3, vpt_1, vpt_2,    &
+          w_1, w_2, w_3
 
     REAL, DIMENSION(:,:,:), POINTER ::                                         &
-          e, e_m, e_p, kh, kh_m, km, km_m, p, pt, pt_m, pt_p, q, q_m, q_p, ql, &
+          e, e_m, e_p, kh, kh_m, km, km_m, pt, pt_m, pt_p, q, q_m, q_p, ql,    &
           ql_c, te_m, tpt_m, tq_m, tu_m, tv_m, tw_m, u, u_m, u_p, v, v_m, v_p, &
           vpt, vpt_m, w, w_m, w_p
@@ -201,7 +204,7 @@
                             scalar_advec = 'pw-scheme'
     CHARACTER (LEN=20)  ::  bc_e_b = 'neumann', bc_lr = 'cyclic', &
-                            bc_ns = 'cyclic', &
-                            bc_p_b = 'neumann', bc_p_t = 'dirichlet', &
-                            bc_pt_b = 'dirichlet', bc_pt_t = 'neumann', &
+                            bc_ns = 'cyclic', bc_p_b = 'neumann', &
+                            bc_p_t = 'dirichlet', bc_pt_b = 'dirichlet', &
+                            bc_pt_t = 'initial_gradient', &
                             bc_q_b = 'dirichlet', bc_q_t = 'neumann', &
                             bc_s_b = 'dirichlet', bc_s_t = 'neumann', &
@@ -263,6 +266,7 @@
                 cloud_droplets = .FALSE., cloud_physics = .FALSE., &
                 conserve_volume_flow = .FALSE., constant_diffusion = .FALSE., &
-                constant_heatflux = .TRUE., constant_waterflux = .TRUE., &
-                create_disturbances = .TRUE., cut_spline_overshoot = .TRUE., &
+                constant_heatflux = .TRUE., constant_top_heatflux = .TRUE., &
+                constant_waterflux = .TRUE., create_disturbances = .TRUE., &
+                cut_spline_overshoot = .TRUE., &
                 data_output_2d_on_each_pe = .TRUE., do2d_at_begin = .FALSE., &
                 do3d_at_begin = .FALSE., do3d_compress = .FALSE., &
@@ -270,5 +274,5 @@
                 disturbance_created = .FALSE., &
                 first_call_advec_particles = .TRUE., &
-                force_print_header = .FALSE., galilei_transformation = .FALSE., &
+                force_print_header = .FALSE., galilei_transformation = .FALSE.,&
                 inflow_l = .FALSE., inflow_n = .FALSE., inflow_r = .FALSE., &
                 inflow_s = .FALSE., iso2d_output = .FALSE., &
@@ -281,6 +285,7 @@
                 random_heatflux = .FALSE., run_control_header = .FALSE., &
                 sloping_surface = .FALSE., stop_dt = .FALSE., &
-                terminate_run = .FALSE., use_prior_plot1d_parameters = .FALSE., &
-                use_surface_fluxes = .FALSE., use_ug_for_galilei_tr = .TRUE., &
+                terminate_run = .FALSE., use_prior_plot1d_parameters = .FALSE.,&
+                use_surface_fluxes = .FALSE., use_top_fluxes = .FALSE., &
+                use_ug_for_galilei_tr = .TRUE., &
                 use_upstream_for_tke = .FALSE., wall_adjustment = .TRUE.
 
@@ -336,8 +341,8 @@
              time_do_sla = 0.0, time_dvrp = 0.0, time_prel = 0.0, &
              time_restart = 9999999.9, time_run_control = 0.0, &
-             ug_surface = 0.0, u_gtrans = 0.0, ups_limit_e = 0.0, &
-             ups_limit_pt = 0.0, ups_limit_u = 0.0, ups_limit_v = 0.0, &
-             ups_limit_w = 0.0, vg_surface = 0.0, v_gtrans = 0.0, &
-             wall_adjustment_factor = 1.8, z_max_do1d = -1.0, &
+             top_heatflux = 9999999.9, ug_surface = 0.0, u_gtrans = 0.0, &
+             ups_limit_e = 0.0, ups_limit_pt = 0.0, ups_limit_u = 0.0, &
+             ups_limit_v = 0.0, ups_limit_w = 0.0, vg_surface = 0.0, &
+             v_gtrans = 0.0, wall_adjustment_factor = 1.8, z_max_do1d = -1.0, &
              z_max_do1d_normalized = -1.0, z_max_do2d = -1.0
 
@@ -500,5 +505,5 @@
     INTEGER ::  ngp_sums, nnx, nx = 0, nxa, nxl, nxr, nxra, nny, ny = 0, nya,  &
                 nyn, nyna, nys, nnz, nz = 0, nza, nzb, nzb_diff, nzt, nzta,    &
-                uxrp = 0, vynp = 0
+                nzt_diff, uxrp = 0, vynp = 0
 
     INTEGER, DIMENSION(:), ALLOCATABLE ::                                      &
Index: palm/trunk/SOURCE/parin.f90
===================================================================
--- palm/trunk/SOURCE/parin.f90	(revision 18)
+++ palm/trunk/SOURCE/parin.f90	(revision 19)
@@ -4,5 +4,5 @@
 ! Actual revisions:
 ! -----------------
-! 
+! +top_heatflux in inipar
 !
 ! Former revisions:
@@ -71,5 +71,5 @@
                        s_surface, s_surface_initial_change, &
                        s_vertical_gradient, s_vertical_gradient_level, &
-                       timestep_scheme, topography, ug_surface, &
+                       top_heatflux, timestep_scheme, topography, ug_surface, &
                        ug_vertical_gradient, ug_vertical_gradient_level, &
                        ups_limit_e, ups_limit_pt, ups_limit_u, ups_limit_v, &
Index: palm/trunk/SOURCE/production_e.f90
===================================================================
--- palm/trunk/SOURCE/production_e.f90	(revision 18)
+++ palm/trunk/SOURCE/production_e.f90	(revision 19)
@@ -4,5 +4,6 @@
 ! Actual revisions:
 ! -----------------
-! 
+! Calculation extended for gridpoint nzt, extended for given temperature /
+! humidity fluxes at the top
 !
 ! Former revisions:
@@ -67,5 +68,5 @@
 
           DO  j = nys, nyn
-             DO  k = nzb_diff_s_outer(j,i), nzt-1
+             DO  k = nzb_diff_s_outer(j,i), nzt
 
                 dudx  =        ( u(k,j,i+1) - u(k,j,i)     ) * ddx
@@ -327,8 +328,9 @@
              DO  j = nys, nyn
 
-                DO  k = nzb_diff_s_inner(j,i), nzt-1
+                DO  k = nzb_diff_s_inner(j,i), nzt_diff
                    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)
                 ENDDO
+
                 IF ( use_surface_fluxes )  THEN
                    k = nzb_diff_s_inner(j,i)-1
@@ -336,4 +338,9 @@
                 ENDIF
 
+                IF ( use_top_fluxes )  THEN
+                   k = nzt
+                   tend(k,j,i) = tend(k,j,i) + g / pt(k,j,i) * tswst(j,i)
+                ENDIF
+
              ENDDO
 
@@ -342,5 +349,5 @@
              DO  j = nys, nyn
 
-                DO  k = nzb_diff_s_inner(j,i), nzt-1
+                DO  k = nzb_diff_s_inner(j,i), nzt_diff
 
                    IF ( .NOT. cloud_physics )  THEN
@@ -375,5 +382,5 @@
                 DO  j = nys, nyn
 
-                   k = nzb_diff_s_inner(j,i)-1
+                   k = nzb_diff_s_inner(j,i)
 
                    IF ( .NOT. cloud_physics )  THEN
@@ -402,4 +409,35 @@
              ENDIF
 
+             IF ( use_top_fluxes )  THEN
+
+                DO  j = nys, nyn
+
+                   k = nzt
+
+                   IF ( .NOT. cloud_physics )  THEN
+                      k1 = 1.0 + 0.61 * q(k,j,i)
+                      k2 = 0.61 * pt(k,j,i)
+                   ELSE
+                      IF ( ql(k,j,i) == 0.0 )  THEN
+                         k1 = 1.0 + 0.61 * q(k,j,i)
+                         k2 = 0.61 * pt(k,j,i)
+                      ELSE
+                         theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
+                         temp  = theta * t_d_pt(k)
+                         k1 = ( 1.0 - q(k,j,i) + 1.61 *                 &
+                                    ( q(k,j,i) - ql(k,j,i) ) *          &
+                              ( 1.0 + 0.622 * l_d_r / temp ) ) /        &
+                              ( 1.0 + 0.622 * l_d_r * l_d_cp *          &
+                              ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
+                         k2 = theta * ( l_d_cp / temp * k1 - 1.0 )
+                      ENDIF
+                   ENDIF
+
+                   tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) * &
+                                         ( k1* tswst(j,i) + k2 * qswst(j,i) )
+                ENDDO
+
+             ENDIF
+
           ENDIF
 
@@ -430,5 +468,5 @@
 !
 !--    Calculate TKE production by shear
-       DO  k = nzb_diff_s_outer(j,i), nzt-1
+       DO  k = nzb_diff_s_outer(j,i), nzt
 
           dudx  =        ( u(k,j,i+1) - u(k,j,i)     ) * ddx
@@ -657,8 +695,9 @@
        IF ( .NOT. moisture )  THEN
 
-          DO  k = nzb_diff_s_inner(j,i), nzt-1
+          DO  k = nzb_diff_s_inner(j,i), nzt_diff
              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)
           ENDDO
+
           IF ( use_surface_fluxes )  THEN
              k = nzb_diff_s_inner(j,i)-1
@@ -666,7 +705,12 @@
           ENDIF
 
+          IF ( use_top_fluxes )  THEN
+             k = nzt
+             tend(k,j,i) = tend(k,j,i) + g / pt(k,j,i) * tswst(j,i)
+          ENDIF
+
        ELSE
 
-          DO  k = nzb_diff_s_inner(j,i), nzt-1
+          DO  k = nzb_diff_s_inner(j,i), nzt_diff
 
              IF ( .NOT. cloud_physics )  THEN
@@ -694,4 +738,5 @@
                                       ) * dd2zu(k)
           ENDDO
+
           IF ( use_surface_fluxes )  THEN
              k = nzb_diff_s_inner(j,i)-1
@@ -720,4 +765,30 @@
           ENDIF
 
+          IF ( use_top_fluxes )  THEN
+             k = nzt
+
+             IF ( .NOT. cloud_physics ) THEN
+                k1 = 1.0 + 0.61 * q(k,j,i)
+                k2 = 0.61 * pt(k,j,i)
+             ELSE
+                IF ( ql(k,j,i) == 0.0 )  THEN
+                   k1 = 1.0 + 0.61 * q(k,j,i)
+                   k2 = 0.61 * pt(k,j,i)
+                ELSE
+                   theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
+                   temp  = theta * t_d_pt(k)
+                   k1 = ( 1.0 - q(k,j,i) + 1.61 *                 &
+                              ( q(k,j,i) - ql(k,j,i) ) *          &
+                        ( 1.0 + 0.622 * l_d_r / temp ) ) /        &
+                        ( 1.0 + 0.622 * l_d_r * l_d_cp *          &
+                        ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
+                   k2 = theta * ( l_d_cp / temp * k1 - 1.0 )
+                ENDIF
+             ENDIF
+
+             tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) * &
+                                         ( k1* tswst(j,i) + k2 * qswst(j,i) )
+          ENDIF
+
        ENDIF
 
Index: palm/trunk/SOURCE/prognostic_equations.f90
===================================================================
--- palm/trunk/SOURCE/prognostic_equations.f90	(revision 18)
+++ palm/trunk/SOURCE/prognostic_equations.f90	(revision 19)
@@ -4,5 +4,6 @@
 ! Actual revisions:
 ! -----------------
-! 
+! Calculation of e, q, and pt extended for gridpoint nzt,
+! handling of given temperature/humidity/scalar fluxes at top surface
 !
 ! Former revisions:
@@ -22,5 +23,5 @@
 ! Description:
 ! ------------
-! Solving the prognostic equations and advecting particles.
+! Solving the prognostic equations.
 !------------------------------------------------------------------------------!
 
@@ -333,5 +334,5 @@
 !--       Tendency terms
           IF ( scalar_advec == 'bc-scheme' )  THEN
-             CALL diffusion_s( i, j, ddzu, ddzw, kh, pt, shf, tend )
+             CALL diffusion_s( i, j, ddzu, ddzw, kh, pt, shf, tswst, tend )
           ELSE
              IF ( tsc(2) == 2.0  .OR.  timestep_scheme(1:5) == 'runge' )  THEN
@@ -346,7 +347,8 @@
              IF ( tsc(2) == 2.0  .AND.  timestep_scheme(1:8) == 'leapfrog' ) &
              THEN
-                CALL diffusion_s( i, j, ddzu, ddzw, kh_m, pt_m, shf_m, tend )
+                CALL diffusion_s( i, j, ddzu, ddzw, kh_m, pt_m, shf_m, &
+                                  tswst_m, tend )
              ELSE
-                CALL diffusion_s( i, j, ddzu, ddzw, kh, pt, shf, tend )
+                CALL diffusion_s( i, j, ddzu, ddzw, kh, pt, shf, tswst, tend )
              ENDIF
           ENDIF
@@ -368,5 +370,5 @@
 !
 !--       Prognostic equation for potential temperature
-          DO  k = nzb_s_inner(j,i)+1, nzt-1
+          DO  k = nzb_s_inner(j,i)+1, nzt
              pt_p(k,j,i) = ( 1 - sat ) * pt_m(k,j,i) + sat * pt(k,j,i) + &
                            dt_3d * (                                           &
@@ -380,10 +382,10 @@
           IF ( timestep_scheme(1:5) == 'runge' )  THEN
              IF ( intermediate_timestep_count == 1 )  THEN
-                DO  k = nzb_s_inner(j,i)+1, nzt-1
+                DO  k = nzb_s_inner(j,i)+1, nzt
                    tpt_m(k,j,i) = tend(k,j,i)
                 ENDDO
              ELSEIF ( intermediate_timestep_count < &
                       intermediate_timestep_count_max )  THEN
-                DO  k = nzb_s_inner(j,i)+1, nzt-1
+                DO  k = nzb_s_inner(j,i)+1, nzt
                    tpt_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tpt_m(k,j,i)
                 ENDDO
@@ -429,5 +431,5 @@
 !--          Tendency-terms
              IF ( scalar_advec == 'bc-scheme' )  THEN
-                CALL diffusion_s( i, j, ddzu, ddzw, kh, q, qsws, tend )
+                CALL diffusion_s( i, j, ddzu, ddzw, kh, q, qsws, qswst, tend )
              ELSE
                 IF ( tsc(2) == 2.0  .OR.  timestep_scheme(1:5) == 'runge' ) THEN
@@ -443,7 +445,8 @@
                 THEN
                    CALL diffusion_s( i, j, ddzu, ddzw, kh_m, q_m, qsws_m, &
+                                     qswst_m, tend )
+                ELSE
+                   CALL diffusion_s( i, j, ddzu, ddzw, kh, q, qsws, qswst, &
                                      tend )
-                ELSE
-                   CALL diffusion_s( i, j, ddzu, ddzw, kh, q, qsws, tend )
                 ENDIF
              ENDIF
@@ -459,5 +462,5 @@
 !
 !--          Prognostic equation for total water content / scalar
-             DO  k = nzb_s_inner(j,i)+1, nzt-1
+             DO  k = nzb_s_inner(j,i)+1, nzt
                 q_p(k,j,i) = ( 1 - sat ) * q_m(k,j,i) + sat * q(k,j,i) +       &
                              dt_3d * (                                         &
@@ -471,10 +474,10 @@
              IF ( timestep_scheme(1:5) == 'runge' )  THEN
                 IF ( intermediate_timestep_count == 1 )  THEN
-                   DO  k = nzb_s_inner(j,i)+1, nzt-1
+                   DO  k = nzb_s_inner(j,i)+1, nzt
                       tq_m(k,j,i) = tend(k,j,i)
                    ENDDO
                 ELSEIF ( intermediate_timestep_count < &
                          intermediate_timestep_count_max )  THEN
-                   DO  k = nzb_s_inner(j,i)+1, nzt-1
+                   DO  k = nzb_s_inner(j,i)+1, nzt
                       tq_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tq_m(k,j,i)
                    ENDDO
@@ -577,5 +580,5 @@
 !--          reasons in the course of the integration. In such cases the old TKE
 !--          value is reduced by 90%.
-             DO  k = nzb_s_inner(j,i)+1, nzt-1
+             DO  k = nzb_s_inner(j,i)+1, nzt
                 e_p(k,j,i) = ( 1 - sat ) * e_m(k,j,i) + sat * e(k,j,i) +       &
                              dt_3d * (                                         &
@@ -589,10 +592,10 @@
              IF ( timestep_scheme(1:5) == 'runge' )  THEN
                 IF ( intermediate_timestep_count == 1 )  THEN
-                   DO  k = nzb_s_inner(j,i)+1, nzt-1
+                   DO  k = nzb_s_inner(j,i)+1, nzt
                       te_m(k,j,i) = tend(k,j,i)
                    ENDDO
                 ELSEIF ( intermediate_timestep_count < &
                          intermediate_timestep_count_max )  THEN
-                   DO  k = nzb_s_inner(j,i)+1, nzt-1
+                   DO  k = nzb_s_inner(j,i)+1, nzt
                       te_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * te_m(k,j,i)
                    ENDDO
@@ -809,7 +812,8 @@
              IF ( tsc(2) == 2.0  .AND.  timestep_scheme(1:8) == 'leapfrog' ) &
              THEN
-                CALL diffusion_s( i, j, ddzu, ddzw, kh_m, pt_m, shf_m, tend )
+                CALL diffusion_s( i, j, ddzu, ddzw, kh_m, pt_m, shf_m, &
+                                  tswst_m, tend )
              ELSE
-                CALL diffusion_s( i, j, ddzu, ddzw, kh, pt, shf, tend )
+                CALL diffusion_s( i, j, ddzu, ddzw, kh, pt, shf, tswst, tend )
              ENDIF
 
@@ -830,5 +834,5 @@
 !
 !--          Prognostic equation for potential temperature
-             DO  k = nzb_s_inner(j,i)+1, nzt-1
+             DO  k = nzb_s_inner(j,i)+1, nzt
                 pt_p(k,j,i) = ( 1.0-tsc(1) ) * pt_m(k,j,i) + tsc(1)*pt(k,j,i) +&
                               dt_3d * (                                        &
@@ -842,10 +846,10 @@
              IF ( timestep_scheme(1:5) == 'runge' )  THEN
                 IF ( intermediate_timestep_count == 1 )  THEN
-                   DO  k = nzb_s_inner(j,i)+1, nzt-1
+                   DO  k = nzb_s_inner(j,i)+1, nzt
                       tpt_m(k,j,i) = tend(k,j,i)
                    ENDDO
                 ELSEIF ( intermediate_timestep_count < &
                          intermediate_timestep_count_max )  THEN
-                   DO  k = nzb_s_inner(j,i)+1, nzt-1
+                   DO  k = nzb_s_inner(j,i)+1, nzt
                       tpt_m(k,j,i) = -9.5625 * tend(k,j,i) + &
                                       5.3125 * tpt_m(k,j,i)
@@ -870,7 +874,9 @@
                 IF ( tsc(2) == 2.0  .AND.  timestep_scheme(1:8) == 'leapfrog' )&
                 THEN
-                   CALL diffusion_s( i, j, ddzu, ddzw, kh_m, q_m, qsws_m, tend )
+                   CALL diffusion_s( i, j, ddzu, ddzw, kh_m, q_m, qsws_m, &
+                                     qswst_m, tend )
                 ELSE
-                   CALL diffusion_s( i, j, ddzu, ddzw, kh, q, qsws, tend )
+                   CALL diffusion_s( i, j, ddzu, ddzw, kh, q, qsws, qswst, &
+                                     tend )
                 ENDIF
        
@@ -885,5 +891,5 @@
 !
 !--             Prognostic equation for total water content / scalar
-                DO  k = nzb_s_inner(j,i)+1, nzt-1
+                DO  k = nzb_s_inner(j,i)+1, nzt
                    q_p(k,j,i) = ( 1.0-tsc(1) ) * q_m(k,j,i) + tsc(1)*q(k,j,i) +&
                                 dt_3d * (                                      &
@@ -897,10 +903,10 @@
                 IF ( timestep_scheme(1:5) == 'runge' )  THEN
                    IF ( intermediate_timestep_count == 1 )  THEN
-                      DO  k = nzb_s_inner(j,i)+1, nzt-1
+                      DO  k = nzb_s_inner(j,i)+1, nzt
                          tq_m(k,j,i) = tend(k,j,i)
                       ENDDO
                    ELSEIF ( intermediate_timestep_count < &
                             intermediate_timestep_count_max )  THEN
-                      DO  k = nzb_s_inner(j,i)+1, nzt-1
+                      DO  k = nzb_s_inner(j,i)+1, nzt
                          tq_m(k,j,i) = -9.5625 * tend(k,j,i) + &
                                         5.3125 * tq_m(k,j,i)
@@ -951,5 +957,5 @@
 !--             reasons in the course of the integration. In such cases the old
 !--             TKE value is reduced by 90%.
-                DO  k = nzb_s_inner(j,i)+1, nzt-1
+                DO  k = nzb_s_inner(j,i)+1, nzt
                    e_p(k,j,i) = ( 1.0-tsc(1) ) * e_m(k,j,i) + tsc(1)*e(k,j,i) +&
                                 dt_3d * (                                      &
@@ -963,10 +969,10 @@
                 IF ( timestep_scheme(1:5) == 'runge' )  THEN
                    IF ( intermediate_timestep_count == 1 )  THEN
-                      DO  k = nzb_s_inner(j,i)+1, nzt-1
+                      DO  k = nzb_s_inner(j,i)+1, nzt
                          te_m(k,j,i) = tend(k,j,i)
                       ENDDO
                    ELSEIF ( intermediate_timestep_count < &
                             intermediate_timestep_count_max )  THEN
-                      DO  k = nzb_s_inner(j,i)+1, nzt-1
+                      DO  k = nzb_s_inner(j,i)+1, nzt
                          te_m(k,j,i) = -9.5625 * tend(k,j,i) + &
                                         5.3125 * te_m(k,j,i)
@@ -1252,5 +1258,5 @@
 !-- pt-tendency terms with no communication
     IF ( scalar_advec == 'bc-scheme' )  THEN
-       CALL diffusion_s( ddzu, ddzw, kh, pt, shf, tend )
+       CALL diffusion_s( ddzu, ddzw, kh, pt, shf, tswst, tend )
     ELSE
        IF ( tsc(2) == 2.0  .OR.  timestep_scheme(1:5) == 'runge' )  THEN
@@ -1264,7 +1270,7 @@
        ENDIF
        IF ( tsc(2) == 2.0  .AND.  timestep_scheme(1:8) == 'leapfrog' )  THEN
-          CALL diffusion_s( ddzu, ddzw, kh_m, pt_m, shf_m, tend )
+          CALL diffusion_s( ddzu, ddzw, kh_m, pt_m, shf_m, tswst_m, tend )
        ELSE
-          CALL diffusion_s( ddzu, ddzw, kh, pt, shf, tend )
+          CALL diffusion_s( ddzu, ddzw, kh, pt, shf, tswst, tend )
        ENDIF
     ENDIF
@@ -1288,5 +1294,5 @@
     DO  i = nxl, nxr
        DO  j = nys, nyn
-          DO  k = nzb_s_inner(j,i)+1, nzt-1
+          DO  k = nzb_s_inner(j,i)+1, nzt
              pt_p(k,j,i) = ( 1 - sat ) * pt_m(k,j,i) + sat * pt(k,j,i) +       &
                            dt_3d * (                                           &
@@ -1304,5 +1310,5 @@
           DO  i = nxl, nxr
              DO  j = nys, nyn
-                DO  k = nzb_s_inner(j,i)+1, nzt-1
+                DO  k = nzb_s_inner(j,i)+1, nzt
                    tpt_m(k,j,i) = tend(k,j,i)
                 ENDDO
@@ -1313,5 +1319,5 @@
           DO  i = nxl, nxr
              DO  j = nys, nyn
-                DO  k = nzb_s_inner(j,i)+1, nzt-1
+                DO  k = nzb_s_inner(j,i)+1, nzt
                    tpt_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tpt_m(k,j,i)
                 ENDDO
@@ -1352,5 +1358,5 @@
 !--    Scalar/q-tendency terms with no communication
        IF ( scalar_advec == 'bc-scheme' )  THEN
-          CALL diffusion_s( ddzu, ddzw, kh, q, qsws, tend )
+          CALL diffusion_s( ddzu, ddzw, kh, q, qsws, qswst, tend )
        ELSE
           IF ( tsc(2) == 2.0  .OR.  timestep_scheme(1:5) == 'runge' )  THEN
@@ -1364,7 +1370,7 @@
           ENDIF
           IF ( tsc(2) == 2.0  .AND.  timestep_scheme(1:8) == 'leapfrog' )  THEN
-             CALL diffusion_s( ddzu, ddzw, kh_m, q_m, qsws_m, tend )
+             CALL diffusion_s( ddzu, ddzw, kh_m, q_m, qsws_m, qswst_m, tend )
           ELSE
-             CALL diffusion_s( ddzu, ddzw, kh, q, qsws, tend )
+             CALL diffusion_s( ddzu, ddzw, kh, q, qsws, qswst, tend )
           ENDIF
        ENDIF
@@ -1382,5 +1388,5 @@
        DO  i = nxl, nxr
           DO  j = nys, nyn
-             DO  k = nzb_s_inner(j,i)+1, nzt-1
+             DO  k = nzb_s_inner(j,i)+1, nzt
                 q_p(k,j,i) = ( 1 - sat ) * q_m(k,j,i) + sat * q(k,j,i) +       &
                              dt_3d * (                                         &
@@ -1398,5 +1404,5 @@
              DO  i = nxl, nxr
                 DO  j = nys, nyn
-                   DO  k = nzb_s_inner(j,i)+1, nzt-1
+                   DO  k = nzb_s_inner(j,i)+1, nzt
                       tq_m(k,j,i) = tend(k,j,i)
                    ENDDO
@@ -1407,5 +1413,5 @@
              DO  i = nxl, nxr
                 DO  j = nys, nyn
-                   DO  k = nzb_s_inner(j,i)+1, nzt-1
+                   DO  k = nzb_s_inner(j,i)+1, nzt
                       tq_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tq_m(k,j,i)
                    ENDDO
@@ -1503,5 +1509,5 @@
        DO  i = nxl, nxr
           DO  j = nys, nyn
-             DO  k = nzb_s_inner(j,i)+1, nzt-1
+             DO  k = nzb_s_inner(j,i)+1, nzt
                 e_p(k,j,i) = ( 1 - sat ) * e_m(k,j,i) + sat * e(k,j,i) +       &
                              dt_3d * (                                         &
@@ -1519,5 +1525,5 @@
              DO  i = nxl, nxr
                 DO  j = nys, nyn
-                   DO  k = nzb_s_inner(j,i)+1, nzt-1
+                   DO  k = nzb_s_inner(j,i)+1, nzt
                       te_m(k,j,i) = tend(k,j,i)
                    ENDDO
@@ -1528,5 +1534,5 @@
              DO  i = nxl, nxr
                 DO  j = nys, nyn
-                   DO  k = nzb_s_inner(j,i)+1, nzt-1
+                   DO  k = nzb_s_inner(j,i)+1, nzt
                       te_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * te_m(k,j,i)
                    ENDDO
Index: palm/trunk/SOURCE/read_3d_binary.f90
===================================================================
--- palm/trunk/SOURCE/read_3d_binary.f90	(revision 18)
+++ palm/trunk/SOURCE/read_3d_binary.f90	(revision 19)
@@ -4,5 +4,5 @@
 ! Actual revisions:
 ! -----------------
-! 
+! +qswst, qswst_m, tswst, tswst_m
 !
 ! Former revisions:
@@ -249,4 +249,8 @@
           CASE ( 'qsws_m' )
              READ ( 13 )  qsws_m
+          CASE ( 'qswst' )
+             READ ( 13 )  qswst
+          CASE ( 'qswst_m' )
+             READ ( 13 )  qswst_m
           CASE ( 'qv_av' )
              ALLOCATE( qv_av(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
@@ -266,4 +270,8 @@
           CASE ( 'shf_m' )
              READ ( 13 )  shf_m
+          CASE ( 'tswst' )
+             READ ( 13 )  tswst
+          CASE ( 'tswst_m' )
+             READ ( 13 )  tswst_m
           CASE ( 'spectrum_x' )
              READ ( 13 )  spectrum_x
Index: palm/trunk/SOURCE/read_var_list.f90
===================================================================
--- palm/trunk/SOURCE/read_var_list.f90	(revision 18)
+++ palm/trunk/SOURCE/read_var_list.f90	(revision 19)
@@ -4,5 +4,5 @@
 ! Actual revisions:
 ! -----------------
-! 
+! +top_heatflux
 !
 ! Former revisions:
@@ -340,4 +340,6 @@
           CASE ( 'topography' )
              READ ( 13 )  topography
+          CASE ( 'top_heatflux' )
+             READ ( 13 )  top_heatflux
           CASE ( 'tsc' )
              READ ( 13 )  tsc
Index: palm/trunk/SOURCE/spline_z.f90
===================================================================
--- palm/trunk/SOURCE/spline_z.f90	(revision 18)
+++ palm/trunk/SOURCE/spline_z.f90	(revision 19)
@@ -4,5 +4,5 @@
 ! Actual revisions:
 ! -----------------
-! 
+! Boundary condition for pt at top adjusted
 !
 ! Former revisions:
@@ -156,10 +156,10 @@
 !
 !--       Top boundary for temperature
-          IF ( ibc_pt_t == 1 )  THEN
-             vad(nzt,:)   = vad(nzt-1,:) + bc_pt_t_val * dz_spline(nzt)
+          IF ( ibc_pt_t == 0 )  THEN
+             vad(nzt+1,:) = pt(nzt+1,nys:nyn,i)
+          ELSEIF ( ibc_pt_t == 1 )  THEN
+             vad(nzt+1,:) = vad(nzt,:)
+          ELSEIF ( ibc_pt_t == 2 )  THEN
              vad(nzt+1,:) = vad(nzt,:)   + bc_pt_t_val * dz_spline(nzt+1)
-          ELSE
-             vad(nzt,:)   = pt(nzt,nys:nyn,i)
-             vad(nzt+1,:) = pt(nzt+1,nys:nyn,i)
           ENDIF
 
Index: palm/trunk/SOURCE/swap_timelevel.f90
===================================================================
--- palm/trunk/SOURCE/swap_timelevel.f90	(revision 18)
+++ palm/trunk/SOURCE/swap_timelevel.f90	(revision 19)
@@ -4,5 +4,5 @@
 ! Actual revisions:
 ! -----------------
-! 
+! Swaping of top fluxes
 !
 ! Former revisions:
@@ -122,4 +122,10 @@
                    rif_m  => rif_1;   rif  => rif_2
                 ENDIF
+                IF ( use_top_fluxes )  THEN
+                   tswst_m => tswst_1;  tswst => tswst_2
+                   IF ( moisture  .OR.  passive_scalar )  THEN
+                      qswst_m => qswst_1;  qswst => qswst_2
+                   ENDIF
+                ENDIF
              ENDIF
 
@@ -161,4 +167,10 @@
                    rif_m  => rif_2;   rif  => rif_1
                 ENDIF
+                IF ( use_top_fluxes )  THEN
+                   tswst_m  => tswst_2;  tswst => tswst_1
+                   IF ( moisture  .OR.  passive_scalar )  THEN
+                      qswst_m => qswst_2;  qswst => qswst_1
+                   ENDIF
+                ENDIF
              ENDIF
 
Index: palm/trunk/SOURCE/write_3d_binary.f90
===================================================================
--- palm/trunk/SOURCE/write_3d_binary.f90	(revision 18)
+++ palm/trunk/SOURCE/write_3d_binary.f90	(revision 19)
@@ -4,5 +4,5 @@
 ! Actual revisions:
 ! -----------------
-! 
+! +qswst, qswst_m, tswst, tswst_m
 !
 ! Former revisions:
@@ -111,4 +111,8 @@
           WRITE ( 14 )  'qsws_m              ';  WRITE ( 14 ) qsws_m
        ENDIF
+       WRITE ( 14 )  'qswst               ';  WRITE ( 14 ) qswst
+       IF ( timestep_scheme(1:5) /= 'runge' )  THEN
+          WRITE ( 14 )  'qswst_m             ';  WRITE ( 14 ) qswst_m
+       ENDIF
     ENDIF
     IF ( ALLOCATED( ql_c_av ) )  THEN
@@ -140,4 +144,8 @@
     IF ( ALLOCATED( ts_av ) )  THEN
        WRITE ( 14 )  'ts_av               ';  WRITE ( 14 )  ts_av
+    ENDIF
+    WRITE ( 14 )  'tswst               ';  WRITE ( 14 )  tswst
+    IF ( timestep_scheme(1:5) /= 'runge' )  THEN
+       WRITE ( 14 )  'tswst_m             ';  WRITE ( 14 )  tswst_m
     ENDIF
     WRITE ( 14 )  'u                   ';  WRITE ( 14 )  u
Index: palm/trunk/SOURCE/write_var_list.f90
===================================================================
--- palm/trunk/SOURCE/write_var_list.f90	(revision 18)
+++ palm/trunk/SOURCE/write_var_list.f90	(revision 19)
@@ -4,5 +4,5 @@
 ! Actual revisions:
 ! -----------------
-! 
+! +top_heatflux
 !
 ! Former revisions:
@@ -306,4 +306,6 @@
     WRITE ( 14 )  'topography                    '
     WRITE ( 14 )  topography
+    WRITE ( 14 )  'top_heatflux                  '
+    WRITE ( 14 )  top_heatflux
     WRITE ( 14 )  'tsc                           '
     WRITE ( 14 )  tsc
