Index: /palm/trunk/SOURCE/modules.f90
===================================================================
--- /palm/trunk/SOURCE/modules.f90	(revision 36)
+++ /palm/trunk/SOURCE/modules.f90	(revision 37)
@@ -199,5 +199,5 @@
     CHARACTER (LEN=8)   ::  run_date, run_time
     CHARACTER (LEN=9)   ::  simulated_time_chr
-    CHARACTER (LEN=12)  ::  version = 'PALM   3.1b'
+    CHARACTER (LEN=12)  ::  version = 'PALM   3.1c'
     CHARACTER (LEN=16)  ::  momentum_advec = 'pw-scheme', &
                             psolver = 'poisfft', &
Index: /palm/trunk/SOURCE/production_e.f90
===================================================================
--- /palm/trunk/SOURCE/production_e.f90	(revision 36)
+++ /palm/trunk/SOURCE/production_e.f90	(revision 37)
@@ -5,5 +5,6 @@
 ! -----------------
 ! Calculation extended for gridpoint nzt, extended for given temperature /
-! humidity fluxes at the top
+! humidity fluxes at the top, wall-part is now executed in case that a
+! Prandtl-layer is switched on (instead of surfaces fluxes switched on)
 !
 ! Former revisions:
@@ -22,4 +23,6 @@
 ! ------------
 ! Production terms (shear + buoyancy) of the TKE
+! WARNING: the case with prandtl_layer = F and use_surface_fluxes = T is
+!          not considered well!
 !------------------------------------------------------------------------------!
 
@@ -99,5 +102,5 @@
           ENDDO
 
-          IF ( use_surface_fluxes )  THEN
+          IF ( prandtl_layer )  THEN
 
 !
@@ -320,4 +323,38 @@
              ENDDO
 
+          ELSEIF ( use_surface_fluxes )  THEN
+
+             DO  j = nys, nyn
+
+                k = nzb_diff_s_outer(j,i)-1
+
+                dudx  =        ( u(k,j,i+1) - u(k,j,i)     ) * ddx
+                dudy  = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
+                                 u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
+                dudz  = 0.5  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
+                                 u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
+
+                dvdx  = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
+                                 v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
+                dvdy  =        ( v(k,j+1,i) - v(k,j,i)     ) * ddy
+                dvdz  = 0.5  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
+                                 v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
+
+                dwdx  = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
+                                 w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
+                dwdy  = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
+                                 w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
+                dwdz  =        ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
+
+                def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
+                      dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + &
+                      dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
+
+                IF ( def < 0.0 )  def = 0.0
+
+                tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
+
+             ENDDO
+
           ENDIF
 
@@ -498,5 +535,5 @@
        ENDDO
 
-       IF ( use_surface_fluxes )  THEN
+       IF ( prandtl_layer )  THEN
 
           IF ( ( wall_e_x(j,i) /= 0.0 ) .OR. ( wall_e_y(j,i) /= 0.0 ) )  THEN
@@ -689,4 +726,34 @@
           ENDIF
 
+       ELSEIF ( use_surface_fluxes )  THEN
+
+          k = nzb_diff_s_outer(j,i)-1
+
+          dudx  =        ( u(k,j,i+1) - u(k,j,i)     ) * ddx
+          dudy  = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
+                           u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
+          dudz  = 0.5  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
+                           u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
+
+          dvdx  = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
+                           v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
+          dvdy  =        ( v(k,j+1,i) - v(k,j,i)     ) * ddy
+          dvdz  = 0.5  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
+                           v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
+
+          dwdx  = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
+                           w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
+          dwdy  = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
+                           w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
+          dwdz  =        ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
+
+          def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
+                dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + &
+                dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
+
+          IF ( def < 0.0 )  def = 0.0
+
+          tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
+
        ENDIF
 
@@ -807,5 +874,5 @@
        INTEGER ::  i, j, ku, kv
 
-       IF ( use_surface_fluxes )  THEN
+       IF ( prandtl_layer )  THEN
 
           IF ( first_call )  THEN
