Index: palm/trunk/SOURCE/production_e.f90
===================================================================
--- palm/trunk/SOURCE/production_e.f90	(revision 4)
+++ 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
 
