Index: palm/trunk/SOURCE/diffusion_e.f90
===================================================================
--- palm/trunk/SOURCE/diffusion_e.f90	(revision 64)
+++ palm/trunk/SOURCE/diffusion_e.f90	(revision 65)
@@ -64,58 +64,59 @@
  
 
-       DO  i = nxl, nxr
-          DO  j = nys, nyn
-!
-!--          First, calculate phi-function for eventually adjusting the &
-!--          mixing length to the prandtl mixing length
-             IF ( adjust_mixing_length  .AND.  prandtl_layer )  THEN
-                IF ( rif(j,i) >= 0.0 )  THEN
-                   phi_m = 1.0 + 5.0 * rif(j,i)
-                ELSE
-                   phi_m = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rif(j,i) ) )
+!
+!--    This if clause must be outside the k-loop because otherwise
+!--    runtime errors occur with -C hopt on NEC
+       IF ( use_pt_reference )  THEN
+
+          DO  i = nxl, nxr
+             DO  j = nys, nyn
+!
+!--             First, calculate phi-function for eventually adjusting the &
+!--             mixing length to the prandtl mixing length
+                IF ( adjust_mixing_length  .AND.  prandtl_layer )  THEN
+                   IF ( rif(j,i) >= 0.0 )  THEN
+                      phi_m = 1.0 + 5.0 * rif(j,i)
+                   ELSE
+                      phi_m = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rif(j,i) ) )
+                   ENDIF
                 ENDIF
-             ENDIF
-
-             DO  k = nzb_s_inner(j,i)+1, nzt
-!
-!--             Calculate the mixing length (for dissipation)
-                dpt_dz = ( theta(k+1,j,i) - theta(k-1,j,i) ) * dd2zu(k)
-                IF ( dpt_dz > 0.0 ) THEN
-                   IF ( use_pt_reference )  THEN
+
+                DO  k = nzb_s_inner(j,i)+1, nzt
+!
+!--                Calculate the mixing length (for dissipation)
+                   dpt_dz = ( theta(k+1,j,i) - theta(k-1,j,i) ) * dd2zu(k)
+                   IF ( dpt_dz > 0.0 ) THEN
                       l_stable = 0.76 * SQRT( e(k,j,i) ) / &
                                         SQRT( g / pt_reference * dpt_dz ) + 1E-5
                    ELSE
-                      l_stable = 0.76 * SQRT( e(k,j,i) ) / &
-                                        SQRT( g / theta(k,j,i) * dpt_dz ) + 1E-5
-                   ENDIF
-                ELSE
-                   l_stable = l_grid(k)
-                ENDIF
-!
-!--             Adjustment of the mixing length
-                IF ( wall_adjustment )  THEN
-                   l(k,j)  = MIN( wall_adjustment_factor * zu(k), l_grid(k), &
-                                  l_stable )
-                   ll(k,j) = MIN( wall_adjustment_factor * zu(k), l_grid(k) )
-                ELSE
-                   l(k,j)  = MIN( l_grid(k), l_stable )
-                   ll(k,j) = l_grid(k)
-                ENDIF
-                IF ( adjust_mixing_length  .AND.  prandtl_layer )  THEN
-                   l(k,j)  = MIN( l(k,j),  kappa * zu(k) / phi_m )
-                   ll(k,j) = MIN( ll(k,j), kappa * zu(k) / phi_m )
-                ENDIF
-
-             ENDDO
-          ENDDO
-!
-!--       Calculate the tendency terms
-          DO  j = nys, nyn
-             DO  k = nzb_s_inner(j,i)+1, nzt
-
-                 dissipation(k,j) = ( 0.19 + 0.74 * l(k,j) / ll(k,j) ) * &
-                                    e(k,j,i) * SQRT( e(k,j,i) ) / l(k,j)
-
-                 tend(k,j,i) = tend(k,j,i)                                     &
+                      l_stable = l_grid(k)
+                   ENDIF
+!
+!--                Adjustment of the mixing length
+                   IF ( wall_adjustment )  THEN
+                      l(k,j)  = MIN( wall_adjustment_factor * zu(k), l_grid(k),&
+                                     l_stable )
+                      ll(k,j) = MIN( wall_adjustment_factor * zu(k), l_grid(k) )
+                   ELSE
+                      l(k,j)  = MIN( l_grid(k), l_stable )
+                      ll(k,j) = l_grid(k)
+                   ENDIF
+                   IF ( adjust_mixing_length  .AND.  prandtl_layer )  THEN
+                      l(k,j)  = MIN( l(k,j),  kappa * zu(k) / phi_m )
+                      ll(k,j) = MIN( ll(k,j), kappa * zu(k) / phi_m )
+                   ENDIF
+
+                ENDDO
+             ENDDO
+
+!
+!--          Calculate the tendency terms
+             DO  j = nys, nyn
+                DO  k = nzb_s_inner(j,i)+1, nzt
+
+                    dissipation(k,j) = ( 0.19 + 0.74 * l(k,j) / ll(k,j) ) * &
+                                       e(k,j,i) * SQRT( e(k,j,i) ) / l(k,j)
+
+                    tend(k,j,i) = tend(k,j,i)                                  &
                                         + (                                    &
                           ( km(k,j,i)+km(k,j,i+1) ) * ( e(k,j,i+1)-e(k,j,i) )  &
@@ -132,19 +133,103 @@
                              - dissipation(k,j)
 
-             ENDDO
+                ENDDO
+             ENDDO
+
+!
+!--          Store dissipation if needed for calculating the sgs particle
+!--          velocities
+             IF ( use_sgs_for_particles )  THEN
+                DO  j = nys, nyn
+                   DO  k = nzb_s_inner(j,i)+1, nzt
+                      diss(k,j,i) = dissipation(k,j)
+                   ENDDO
+                ENDDO
+             ENDIF
+
           ENDDO
 
-!
-!--       Store dissipation if needed for calculating the sgs particle
-!--       velocities
-          IF ( use_sgs_for_particles )  THEN
-             DO  j = nys, nyn
+       ELSE
+
+          DO  i = nxl, nxr
+             DO  j = nys, nyn
+!
+!--             First, calculate phi-function for eventually adjusting the &
+!--             mixing length to the prandtl mixing length
+                IF ( adjust_mixing_length  .AND.  prandtl_layer )  THEN
+                   IF ( rif(j,i) >= 0.0 )  THEN
+                      phi_m = 1.0 + 5.0 * rif(j,i)
+                   ELSE
+                      phi_m = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rif(j,i) ) )
+                   ENDIF
+                ENDIF
+
                 DO  k = nzb_s_inner(j,i)+1, nzt
-                   diss(k,j,i) = dissipation(k,j)
-                ENDDO
-             ENDDO
-          ENDIF
-
-       ENDDO
+!
+!--                Calculate the mixing length (for dissipation)
+                   dpt_dz = ( theta(k+1,j,i) - theta(k-1,j,i) ) * dd2zu(k)
+                   IF ( dpt_dz > 0.0 ) THEN
+                      l_stable = 0.76 * SQRT( e(k,j,i) ) / &
+                                        SQRT( g / theta(k,j,i) * dpt_dz ) + 1E-5
+                   ELSE
+                      l_stable = l_grid(k)
+                   ENDIF
+!
+!--                Adjustment of the mixing length
+                   IF ( wall_adjustment )  THEN
+                      l(k,j)  = MIN( wall_adjustment_factor * zu(k), l_grid(k),&
+                                     l_stable )
+                      ll(k,j) = MIN( wall_adjustment_factor * zu(k), l_grid(k) )
+                   ELSE
+                      l(k,j)  = MIN( l_grid(k), l_stable )
+                      ll(k,j) = l_grid(k)
+                   ENDIF
+                   IF ( adjust_mixing_length  .AND.  prandtl_layer )  THEN
+                      l(k,j)  = MIN( l(k,j),  kappa * zu(k) / phi_m )
+                      ll(k,j) = MIN( ll(k,j), kappa * zu(k) / phi_m )
+                   ENDIF
+
+                ENDDO
+             ENDDO
+
+!
+!--          Calculate the tendency terms
+             DO  j = nys, nyn
+                DO  k = nzb_s_inner(j,i)+1, nzt
+
+                    dissipation(k,j) = ( 0.19 + 0.74 * l(k,j) / ll(k,j) ) * &
+                                       e(k,j,i) * SQRT( e(k,j,i) ) / l(k,j)
+
+                    tend(k,j,i) = tend(k,j,i)                                  &
+                                        + (                                    &
+                          ( km(k,j,i)+km(k,j,i+1) ) * ( e(k,j,i+1)-e(k,j,i) )  &
+                        - ( km(k,j,i)+km(k,j,i-1) ) * ( e(k,j,i)-e(k,j,i-1) )  &
+                                          ) * ddx2                             &
+                                        + (                                    &
+                          ( km(k,j,i)+km(k,j+1,i) ) * ( e(k,j+1,i)-e(k,j,i) )  &
+                        - ( km(k,j,i)+km(k,j-1,i) ) * ( e(k,j,i)-e(k,j-1,i) )  &
+                                          ) * ddy2                             &
+                                        + (                                    &
+               ( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1) &
+             - ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k)   &
+                                          ) * ddzw(k)                          &
+                             - dissipation(k,j)
+
+                ENDDO
+             ENDDO
+
+!
+!--          Store dissipation if needed for calculating the sgs particle
+!--          velocities
+             IF ( use_sgs_for_particles )  THEN
+                DO  j = nys, nyn
+                   DO  k = nzb_s_inner(j,i)+1, nzt
+                      diss(k,j,i) = dissipation(k,j)
+                   ENDDO
+                ENDDO
+             ENDIF
+
+          ENDDO
+
+       ENDIF
 
 !
