MODULE production_e_mod

!------------------------------------------------------------------------------!
! Current revisions:
! -----------------
!
! Former revisions:
! -----------------
! $Id: production_e.f90 668 2010-12-23 13:22:58Z gryschka $
!
! 667 2010-12-23 12:06:00Z suehring/gryschka
! nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng
!
! 449 2010-02-02 11:23:59Z raasch
! test output from rev 410 removed
!
! 388 2009-09-23 09:40:33Z raasch
! Bugfix: wrong sign in buoyancy production of ocean part in case of not using
!         the reference density (only in 3D routine production_e)
! Bugfix to avoid zero division by km_neutral
!
! 208 2008-10-20 06:02:59Z raasch
! Bugfix concerning the calculation of velocity gradients at vertical walls 
! in case of diabatic conditions
!
! 187 2008-08-06 16:25:09Z letzel
! Change: add 'minus' sign to fluxes obtained from subroutine wall_fluxes_e for
! consistency with subroutine wall_fluxes
!
! 124 2007-10-19 15:47:46Z raasch
! Bugfix: calculation of density flux in the ocean now starts from nzb+1
!
! 108 2007-08-24 15:10:38Z letzel
! Bugfix: wrong sign removed from the buoyancy production term in the case
! use_reference = .T.,
! u_0 and v_0 are calculated for nxr+1, nyn+1 also (otherwise these values are
! not available in case of non-cyclic boundary conditions)
! Bugfix for ocean density flux at bottom
!
! 97 2007-06-21 08:23:15Z raasch
! Energy production by density flux (in ocean) added
! use_pt_reference renamed use_reference
!
! 75 2007-03-22 09:54:05Z raasch
! Wall functions now include diabatic conditions, call of routine wall_fluxes_e,
! reference temperature pt_reference can be used in buoyancy term,
! moisture renamed humidity
!
! 37 2007-03-01 08:33:54Z raasch
! Calculation extended for gridpoint nzt, extended for given temperature /
! 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)
!
! RCS Log replace by Id keyword, revision history cleaned up
!
! Revision 1.21  2006/04/26 12:45:35  raasch
! OpenMP parallelization of production_e_init
!
! Revision 1.1  1997/09/19 07:45:35  raasch
! Initial revision
!
!
! Description:
! ------------
! Production terms (shear + buoyancy) of the TKE
! WARNING: the case with prandtl_layer = F and use_surface_fluxes = T is
!          not considered well!
!------------------------------------------------------------------------------!

    USE wall_fluxes_mod

    PRIVATE
    PUBLIC production_e, production_e_init

    LOGICAL, SAVE ::  first_call = .TRUE.

    REAL, DIMENSION(:,:), ALLOCATABLE, SAVE ::  u_0, v_0

    INTERFACE production_e
       MODULE PROCEDURE production_e
       MODULE PROCEDURE production_e_ij
    END INTERFACE production_e
    
    INTERFACE production_e_init
       MODULE PROCEDURE production_e_init
    END INTERFACE production_e_init
 
 CONTAINS


!------------------------------------------------------------------------------!
! Call for all grid points
!------------------------------------------------------------------------------!
    SUBROUTINE production_e

       USE arrays_3d
       USE cloud_parameters
       USE control_parameters
       USE grid_variables
       USE indices
       USE statistics

       IMPLICIT NONE

       INTEGER ::  i, j, k

       REAL    ::  def, dudx, dudy, dudz, dvdx, dvdy, dvdz, dwdx, dwdy, dwdz, &
                   k1, k2, km_neutral, theta, temp

!       REAL, DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) ::  usvs, vsus, wsus, wsvs
       REAL, DIMENSION(nzb:nzt+1) ::   usvs, vsus, wsus, wsvs

!
!--    First calculate horizontal momentum flux u'v', w'v', v'u', w'u' at
!--    vertical walls, if neccessary
!--    So far, results are slightly different from the ij-Version.
!--    Therefore, ij-Version is called further below within the ij-loops.
!       IF ( topography /= 'flat' )  THEN
!          CALL wall_fluxes_e( usvs, 1.0, 0.0, 0.0, 0.0, wall_e_y )
!          CALL wall_fluxes_e( wsvs, 0.0, 0.0, 1.0, 0.0, wall_e_y )
!          CALL wall_fluxes_e( vsus, 0.0, 1.0, 0.0, 0.0, wall_e_x )
!          CALL wall_fluxes_e( wsus, 0.0, 0.0, 0.0, 1.0, wall_e_x )
!       ENDIF

!
!--    Calculate TKE production by shear
       DO  i = nxl, nxr

          DO  j = nys, nyn
             DO  k = nzb_diff_s_outer(j,i), nzt

                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
          ENDDO

          IF ( prandtl_layer )  THEN

!
!--          Position beneath wall
!--          (2) - Will allways be executed.
!--          'bottom and wall: use u_0,v_0 and wall functions'
             DO  j = nys, nyn

                IF ( ( wall_e_x(j,i) /= 0.0 ) .OR. ( wall_e_y(j,i) /= 0.0 ) ) &
                THEN

                   k = nzb_diff_s_inner(j,i) - 1
                   dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
                   dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &
                                  u_0(j,i)   - u_0(j,i+1)   ) * dd2zu(k)
                   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_0(j,i)   - v_0(j+1,i)   ) * dd2zu(k)
                   dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)

                   IF ( wall_e_y(j,i) /= 0.0 )  THEN
!                     
!--                   Inconsistency removed: as the thermal stratification is 
!--                   not taken into account for the evaluation of the wall 
!--                   fluxes at vertical walls, the eddy viscosity km must not 
!--                   be used for the evaluation of the velocity gradients dudy 
!--                   and dwdy
!--                   Note: The validity of the new method has not yet been 
!--                         shown, as so far no suitable data for a validation 
!--                         has been available
                      CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
                                          usvs, 1.0, 0.0, 0.0, 0.0 )
                      CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
                                          wsvs, 0.0, 0.0, 1.0, 0.0 )
                      km_neutral = kappa * ( usvs(k)**2 + wsvs(k)**2 )**0.25 * &
                                   0.5 * dy 
                      IF ( km_neutral > 0.0 )  THEN
                         dudy = - wall_e_y(j,i) * usvs(k) / km_neutral
                         dwdy = - wall_e_y(j,i) * wsvs(k) / km_neutral
                      ELSE
                         dudy = 0.0
                         dwdy = 0.0
                      ENDIF
                   ELSE
                      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
                      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
                   ENDIF

                   IF ( wall_e_x(j,i) /= 0.0 )  THEN
!                     
!--                   Inconsistency removed: as the thermal stratification is 
!--                   not taken into account for the evaluation of the wall 
!--                   fluxes at vertical walls, the eddy viscosity km must not 
!--                   be used for the evaluation of the velocity gradients dvdx 
!--                   and dwdx
!--                   Note: The validity of the new method has not yet been 
!--                         shown, as so far no suitable data for a validation 
!--                         has been available
                      CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
                                          vsus, 0.0, 1.0, 0.0, 0.0 )
                      CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
                                          wsus, 0.0, 0.0, 0.0, 1.0 )
                      km_neutral = kappa * ( vsus(k)**2 + wsus(k)**2 )**0.25 * &
                                   0.5 * dx
                      IF ( km_neutral > 0.0 )  THEN
                         dvdx = - wall_e_x(j,i) * vsus(k) / km_neutral
                         dwdx = - wall_e_x(j,i) * wsus(k) / km_neutral
                      ELSE
                         dvdx = 0.0
                         dwdx = 0.0
                      ENDIF
                   ELSE
                      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
                      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
                   ENDIF

                   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


!
!--                (3) - will be executed only, if there is at least one level
!--                between (2) and (4), i.e. the topography must have a
!--                minimum height of 2 dz. Wall fluxes for this case have
!--                already been calculated for (2).
!--                'wall only: use wall functions'

                   DO  k = nzb_diff_s_inner(j,i), nzb_diff_s_outer(j,i)-2

                      dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
                      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)
                      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)
                      dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)

                      IF ( wall_e_y(j,i) /= 0.0 )  THEN
!                     
!--                      Inconsistency removed: as the thermal stratification 
!--                      is not taken into account for the evaluation of the 
!--                      wall fluxes at vertical walls, the eddy viscosity km 
!--                      must not be used for the evaluation of the velocity 
!--                      gradients dudy and dwdy
!--                      Note: The validity of the new method has not yet 
!--                            been shown, as so far no suitable data for a 
!--                            validation has been available
                         km_neutral = kappa * ( usvs(k)**2 + & 
                                                wsvs(k)**2 )**0.25 * 0.5 * dy
                         IF ( km_neutral > 0.0 )  THEN
                            dudy = - wall_e_y(j,i) * usvs(k) / km_neutral
                            dwdy = - wall_e_y(j,i) * wsvs(k) / km_neutral
                         ELSE
                            dudy = 0.0
                            dwdy = 0.0
                         ENDIF
                      ELSE
                         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
                         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
                      ENDIF

                      IF ( wall_e_x(j,i) /= 0.0 )  THEN
!                     
!--                      Inconsistency removed: as the thermal stratification 
!--                      is not taken into account for the evaluation of the 
!--                      wall fluxes at vertical walls, the eddy viscosity km 
!--                      must not be used for the evaluation of the velocity 
!--                      gradients dvdx and dwdx
!--                      Note: The validity of the new method has not yet 
!--                            been shown, as so far no suitable data for a 
!--                            validation has been available
                         km_neutral = kappa * ( vsus(k)**2 + & 
                                                wsus(k)**2 )**0.25 * 0.5 * dx
                         IF ( km_neutral > 0.0 )  THEN
                            dvdx = - wall_e_x(j,i) * vsus(k) / km_neutral
                            dwdx = - wall_e_x(j,i) * wsus(k) / km_neutral
                         ELSE
                            dvdx = 0.0
                            dwdx = 0.0
                         ENDIF
                      ELSE
                         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
                         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
                      ENDIF

                      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

             ENDDO

!
!--          (4) - will allways be executed.
!--          'special case: free atmosphere' (as for case (0))
             DO  j = nys, nyn

                IF ( ( wall_e_x(j,i) /= 0.0 ) .OR. ( wall_e_y(j,i) /= 0.0 ) ) &
                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

             ENDDO

!
!--          Position without adjacent wall
!--          (1) - will allways be executed. 
!--          'bottom only: use u_0,v_0'
             DO  j = nys, nyn

                IF ( ( wall_e_x(j,i) == 0.0 ) .AND. ( wall_e_y(j,i) == 0.0 ) ) &
                THEN

                   k = nzb_diff_s_inner(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_0(j,i)   - u_0(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_0(j,i)   - v_0(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

             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

!
!--       Calculate TKE production by buoyancy
          IF ( .NOT. humidity )  THEN

             IF ( use_reference )  THEN

                IF ( ocean )  THEN
!
!--                So far in the ocean no special treatment of density flux in
!--                the bottom and top surface layer
                   DO  j = nys, nyn
                      DO  k = nzb_s_inner(j,i)+1, nzt
                         tend(k,j,i) = tend(k,j,i) +                   &
                                       kh(k,j,i) * g / rho_reference * &
                                       ( rho(k+1,j,i)-rho(k-1,j,i) ) * dd2zu(k)
                      ENDDO
                   ENDDO

                ELSE

                   DO  j = nys, nyn
                      DO  k = nzb_diff_s_inner(j,i), nzt_diff
                         tend(k,j,i) = tend(k,j,i) -                  &
                                       kh(k,j,i) * g / pt_reference * &
                                       ( 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
                         tend(k,j,i) = tend(k,j,i) + g / pt_reference * shf(j,i)
                      ENDIF

                      IF ( use_top_fluxes )  THEN
                         k = nzt
                         tend(k,j,i) = tend(k,j,i) + g / pt_reference * &
                                                     tswst(j,i)
                      ENDIF
                   ENDDO

                ENDIF

             ELSE

                IF ( ocean )  THEN
!
!--                So far in the ocean no special treatment of density flux in
!--                the bottom and top surface layer
                   DO  j = nys, nyn
                      DO  k = nzb_s_inner(j,i)+1, nzt
                         tend(k,j,i) = tend(k,j,i) +                &
                                       kh(k,j,i) * g / rho(k,j,i) * &
                                       ( rho(k+1,j,i)-rho(k-1,j,i) ) * dd2zu(k)
                      ENDDO
                   ENDDO

                ELSE

                   DO  j = nys, nyn
                      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
                         tend(k,j,i) = tend(k,j,i) + g / pt(k,j,i) * shf(j,i)
                      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

                ENDIF

             ENDIF

          ELSE

             DO  j = nys, nyn

                DO  k = nzb_diff_s_inner(j,i), nzt_diff

                   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) - kh(k,j,i) * g / vpt(k,j,i) * &
                                      ( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) +  &
                                        k2 * ( q(k+1,j,i) - q(k-1,j,i) )    &
                                      ) * dd2zu(k)
                ENDDO

             ENDDO

             IF ( use_surface_fluxes )  THEN

                DO  j = nys, nyn

                   k = nzb_diff_s_inner(j,i)-1

                   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* shf(j,i) + k2 * qsws(j,i) )
                ENDDO

             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

       ENDDO

    END SUBROUTINE production_e


!------------------------------------------------------------------------------!
! Call for grid point i,j
!------------------------------------------------------------------------------!
    SUBROUTINE production_e_ij( i, j )

       USE arrays_3d
       USE cloud_parameters
       USE control_parameters
       USE grid_variables
       USE indices
       USE statistics

       IMPLICIT NONE

       INTEGER ::  i, j, k

       REAL    ::  def, dudx, dudy, dudz, dvdx, dvdy, dvdz, dwdx, dwdy, dwdz, &
                   k1, k2, km_neutral, theta, temp

       REAL, DIMENSION(nzb:nzt+1) ::  usvs, vsus, wsus, wsvs

!
!--    Calculate TKE production by shear
       DO  k = nzb_diff_s_outer(j,i), nzt

          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

       IF ( prandtl_layer )  THEN

          IF ( ( wall_e_x(j,i) /= 0.0 ) .OR. ( wall_e_y(j,i) /= 0.0 ) )  THEN

!
!--          Position beneath wall
!--          (2) - Will allways be executed.
!--          'bottom and wall: use u_0,v_0 and wall functions'
             k = nzb_diff_s_inner(j,i)-1

             dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
             dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &
                            u_0(j,i)   - u_0(j,i+1)   ) * dd2zu(k)
             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_0(j,i)   - v_0(j+1,i)   ) * dd2zu(k)
             dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)

             IF ( wall_e_y(j,i) /= 0.0 )  THEN
!                     
!--             Inconsistency removed: as the thermal stratification 
!--             is not taken into account for the evaluation of the 
!--             wall fluxes at vertical walls, the eddy viscosity km 
!--             must not be used for the evaluation of the velocity 
!--             gradients dudy and dwdy
!--             Note: The validity of the new method has not yet 
!--                   been shown, as so far no suitable data for a 
!--                   validation has been available
                CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
                                    usvs, 1.0, 0.0, 0.0, 0.0 )
                CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
                                    wsvs, 0.0, 0.0, 1.0, 0.0 )
                km_neutral = kappa * ( usvs(k)**2 + wsvs(k)**2 )**0.25 * &
                             0.5 * dy
                IF ( km_neutral > 0.0 )  THEN
                   dudy = - wall_e_y(j,i) * usvs(k) / km_neutral
                   dwdy = - wall_e_y(j,i) * wsvs(k) / km_neutral
                ELSE
                   dudy = 0.0
                   dwdy = 0.0
                ENDIF
             ELSE
                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
                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
             ENDIF

             IF ( wall_e_x(j,i) /= 0.0 )  THEN
!                     
!--             Inconsistency removed: as the thermal stratification 
!--             is not taken into account for the evaluation of the 
!--             wall fluxes at vertical walls, the eddy viscosity km 
!--             must not be used for the evaluation of the velocity 
!--             gradients dvdx and dwdx
!--             Note: The validity of the new method has not yet 
!--                   been shown, as so far no suitable data for a 
!--                   validation has been available
                CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
                                    vsus, 0.0, 1.0, 0.0, 0.0 )
                CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
                                    wsus, 0.0, 0.0, 0.0, 1.0 )
                km_neutral = kappa * ( vsus(k)**2 + wsus(k)**2 )**0.25 * & 
                             0.5 * dx
                IF ( km_neutral > 0.0 )  THEN
                   dvdx = - wall_e_x(j,i) * vsus(k) / km_neutral
                   dwdx = - wall_e_x(j,i) * wsus(k) / km_neutral
                ELSE
                   dvdx = 0.0
                   dwdx = 0.0
                ENDIF
             ELSE
                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
                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
             ENDIF

             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

!
!--          (3) - will be executed only, if there is at least one level
!--          between (2) and (4), i.e. the topography must have a
!--          minimum height of 2 dz. Wall fluxes for this case have
!--          already been calculated for (2).
!--          'wall only: use wall functions'
             DO  k = nzb_diff_s_inner(j,i), nzb_diff_s_outer(j,i)-2

                dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
                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)
                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)
                dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)

                IF ( wall_e_y(j,i) /= 0.0 )  THEN
!                     
!--                Inconsistency removed: as the thermal stratification 
!--                is not taken into account for the evaluation of the 
!--                wall fluxes at vertical walls, the eddy viscosity km 
!--                must not be used for the evaluation of the velocity 
!--                gradients dudy and dwdy
!--                Note: The validity of the new method has not yet 
!--                      been shown, as so far no suitable data for a 
!--                      validation has been available
                   km_neutral = kappa * ( usvs(k)**2 + & 
                                          wsvs(k)**2 )**0.25 * 0.5 * dy
                   IF ( km_neutral > 0.0 )  THEN
                      dudy = - wall_e_y(j,i) * usvs(k) / km_neutral
                      dwdy = - wall_e_y(j,i) * wsvs(k) / km_neutral
                   ELSE
                      dudy = 0.0
                      dwdy = 0.0
                   ENDIF
                ELSE
                   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
                   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
                ENDIF

                IF ( wall_e_x(j,i) /= 0.0 )  THEN
!                     
!--                Inconsistency removed: as the thermal stratification 
!--                is not taken into account for the evaluation of the 
!--                wall fluxes at vertical walls, the eddy viscosity km 
!--                must not be used for the evaluation of the velocity 
!--                gradients dvdx and dwdx
!--                Note: The validity of the new method has not yet 
!--                      been shown, as so far no suitable data for a 
!--                      validation has been available
                   km_neutral = kappa * ( vsus(k)**2 + & 
                                          wsus(k)**2 )**0.25 * 0.5 * dx
                   IF ( km_neutral > 0.0 )  THEN
                      dvdx = - wall_e_x(j,i) * vsus(k) / km_neutral
                      dwdx = - wall_e_x(j,i) * wsus(k) / km_neutral
                   ELSE
                      dvdx = 0.0
                      dwdx = 0.0
                   ENDIF
                ELSE
                   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
                   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
                ENDIF

                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

!
!--          (4) - will allways be executed.
!--          'special case: free atmosphere' (as for case (0))
             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

          ELSE

!
!--          Position without adjacent wall
!--          (1) - will allways be executed. 
!--          'bottom only: use u_0,v_0'
             k = nzb_diff_s_inner(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_0(j,i)   - u_0(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_0(j,i)   - v_0(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

       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

!
!--    Calculate TKE production by buoyancy
       IF ( .NOT. humidity )  THEN

          IF ( use_reference )  THEN

             IF ( ocean )  THEN
!
!--             So far in the ocean no special treatment of density flux in the
!--             bottom and top surface layer
                DO  k = nzb_s_inner(j,i)+1, nzt
                   tend(k,j,i) = tend(k,j,i) + kh(k,j,i) * g / rho_reference * &
                                      ( rho(k+1,j,i) - rho(k-1,j,i) ) * dd2zu(k)
                ENDDO

             ELSE

                DO  k = nzb_diff_s_inner(j,i), nzt_diff
                   tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / pt_reference * &
                                       ( 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
                   tend(k,j,i) = tend(k,j,i) + g / pt_reference * shf(j,i)
                ENDIF

                IF ( use_top_fluxes )  THEN
                   k = nzt
                   tend(k,j,i) = tend(k,j,i) + g / pt_reference * tswst(j,i)
                ENDIF

             ENDIF

          ELSE

             IF ( ocean )  THEN
!
!--             So far in the ocean no special treatment of density flux in the
!--             bottom and top surface layer
                DO  k = nzb_s_inner(j,i)+1, nzt
                   tend(k,j,i) = tend(k,j,i) + kh(k,j,i) * g / rho(k,j,i) * &
                                      ( rho(k+1,j,i) - rho(k-1,j,i) ) * dd2zu(k)
                ENDDO

             ELSE

                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
                   tend(k,j,i) = tend(k,j,i) + g / pt(k,j,i) * shf(j,i)
                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

             ENDIF

          ENDIF

       ELSE

          DO  k = nzb_diff_s_inner(j,i), nzt_diff

             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) - kh(k,j,i) * g / vpt(k,j,i) *      &
                                      ( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) + &
                                        k2 * ( q(k+1,j,i) - q(k-1,j,i) )   &
                                      ) * dd2zu(k)
          ENDDO

          IF ( use_surface_fluxes )  THEN
             k = nzb_diff_s_inner(j,i)-1

             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* shf(j,i) + k2 * qsws(j,i) )
          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

    END SUBROUTINE production_e_ij


    SUBROUTINE production_e_init

       USE arrays_3d
       USE control_parameters
       USE grid_variables
       USE indices

       IMPLICIT NONE

       INTEGER ::  i, j, ku, kv

       IF ( prandtl_layer )  THEN

          IF ( first_call )  THEN
             ALLOCATE( u_0(nysg:nyng,nxlg:nxrg), &
                       v_0(nysg:nyng,nxlg:nxrg) )
             first_call = .FALSE.
          ENDIF

!
!--       Calculate a virtual velocity at the surface in a way that the
!--       vertical velocity gradient at k = 1 (u(k+1)-u_0) matches the
!--       Prandtl law (-w'u'/km). This gradient is used in the TKE shear
!--       production term at k=1 (see production_e_ij).
!--       The velocity gradient has to be limited in case of too small km
!--       (otherwise the timestep may be significantly reduced by large
!--       surface winds).
!--       Upper bounds are nxr+1 and nyn+1 because otherwise these values are
!--       not available in case of non-cyclic boundary conditions.
!--       WARNING: the exact analytical solution would require the determination
!--                of the eddy diffusivity by km = u* * kappa * zp / phi_m.
          !$OMP PARALLEL DO PRIVATE( ku, kv )
          DO  i = nxl, nxr+1
             DO  j = nys, nyn+1

                ku = nzb_u_inner(j,i)+1
                kv = nzb_v_inner(j,i)+1

                u_0(j,i) = u(ku+1,j,i) + usws(j,i) * ( zu(ku+1) - zu(ku-1) ) / &
                                 ( 0.5 * ( km(ku,j,i) + km(ku,j,i-1) ) +       &
                                   1.0E-20 )
!                                  ( us(j,i) * kappa * zu(1) )
                v_0(j,i) = v(kv+1,j,i) + vsws(j,i) * ( zu(kv+1) - zu(kv-1) ) / &
                                 ( 0.5 * ( km(kv,j,i) + km(kv,j-1,i) ) +       &
                                   1.0E-20 )
!                                  ( us(j,i) * kappa * zu(1) )

                IF ( ABS( u(ku+1,j,i) - u_0(j,i) )  > &
                     ABS( u(ku+1,j,i) - u(ku-1,j,i) ) )  u_0(j,i) = u(ku-1,j,i)
                IF ( ABS( v(kv+1,j,i) - v_0(j,i) )  > &
                     ABS( v(kv+1,j,i) - v(kv-1,j,i) ) )  v_0(j,i) = v(kv-1,j,i)

             ENDDO
          ENDDO

          CALL exchange_horiz_2d( u_0 )
          CALL exchange_horiz_2d( v_0 )

       ENDIF

    END SUBROUTINE production_e_init

 END MODULE production_e_mod