Index: palm/trunk/SOURCE/CURRENT_MODIFICATIONS
===================================================================
--- palm/trunk/SOURCE/CURRENT_MODIFICATIONS	(revision 363)
+++ palm/trunk/SOURCE/CURRENT_MODIFICATIONS	(revision 364)
@@ -113,4 +113,6 @@
 Errors:
 ------
+Bugfix to avoid zero division by km_neutral (production_e)
+
 Bugfix for extended NetCDF files: In order to avoid 'data mode' errors if 
 updated attributes are larger than their original size, NF90_PUT_ATT is called
@@ -163,3 +165,3 @@
 Bugfix: initial setting of time_coupling in coupled restart runs (time_integration)
 
-advec_particles, check_parameters, cpu_log, data_output_2d, data_output_3d, header, init_3d_model, init_particles, modules, netcdf, prandtl_fluxes, read_var_list, time_integration, write_var_list
+advec_particles, check_parameters, cpu_log, data_output_2d, data_output_3d, header, init_3d_model, init_particles, modules, netcdf, prandtl_fluxes, production_e, read_var_list, time_integration, write_var_list
Index: palm/trunk/SOURCE/production_e.f90
===================================================================
--- palm/trunk/SOURCE/production_e.f90	(revision 363)
+++ palm/trunk/SOURCE/production_e.f90	(revision 364)
@@ -4,5 +4,5 @@
 ! Actual revisions:
 ! -----------------
-! 
+! Bugfix to avoid zero division by km_neutral
 !
 ! Former revisions:
@@ -185,6 +185,11 @@
                       km_neutral = kappa * ( usvs(k)**2 + wsvs(k)**2 )**0.25 * &
                                    0.5 * dy 
-                      dudy = - wall_e_y(j,i) * usvs(k) / km_neutral
-                      dwdy = - wall_e_y(j,i) * wsvs(k) / km_neutral
+                      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) - &
@@ -210,6 +215,11 @@
                       km_neutral = kappa * ( vsus(k)**2 + wsus(k)**2 )**0.25 * &
                                    0.5 * dx
-                      dvdx = - wall_e_x(j,i) * vsus(k) / km_neutral
-                      dwdx = - wall_e_x(j,i) * wsus(k) / km_neutral
+                      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) - &
@@ -257,6 +267,11 @@
                          km_neutral = kappa * ( usvs(k)**2 + & 
                                                 wsvs(k)**2 )**0.25 * 0.5 * dy
-                         dudy = - wall_e_y(j,i) * usvs(k) / km_neutral
-                         dwdy = - wall_e_y(j,i) * wsvs(k) / km_neutral
+                         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) - &
@@ -278,6 +293,11 @@
                          km_neutral = kappa * ( vsus(k)**2 + & 
                                                 wsus(k)**2 )**0.25 * 0.5 * dx
-                         dvdx = - wall_e_x(j,i) * vsus(k) / km_neutral
-                         dwdx = - wall_e_x(j,i) * wsus(k) / km_neutral
+                         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) - &
@@ -686,6 +706,11 @@
                 km_neutral = kappa * ( usvs(k)**2 + wsvs(k)**2 )**0.25 * &
                              0.5 * dy
-                dudy = - wall_e_y(j,i) * usvs(k) / km_neutral
-                dwdy = - wall_e_y(j,i) * wsvs(k) / km_neutral
+                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) - &
@@ -711,6 +736,11 @@
                 km_neutral = kappa * ( vsus(k)**2 + wsus(k)**2 )**0.25 * & 
                              0.5 * dx
-                dvdx = - wall_e_x(j,i) * vsus(k) / km_neutral
-                dwdx = - wall_e_x(j,i) * wsus(k) / km_neutral
+                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) - &
@@ -756,6 +786,11 @@
                    km_neutral = kappa * ( usvs(k)**2 + & 
                                           wsvs(k)**2 )**0.25 * 0.5 * dy
-                   dudy = - wall_e_y(j,i) * usvs(k) / km_neutral
-                   dwdy = - wall_e_y(j,i) * wsvs(k) / km_neutral
+                   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) - &
@@ -777,6 +812,11 @@
                    km_neutral = kappa * ( vsus(k)**2 + & 
                                           wsus(k)**2 )**0.25 * 0.5 * dx
-                   dvdx = - wall_e_x(j,i) * vsus(k) / km_neutral
-                   dwdx = - wall_e_x(j,i) * wsus(k) / km_neutral
+                   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) - &
