Index: palm/trunk/SOURCE/surface_layer_fluxes.f90
===================================================================
--- palm/trunk/SOURCE/surface_layer_fluxes.f90	(revision 1706)
+++ palm/trunk/SOURCE/surface_layer_fluxes.f90	(revision 1709)
@@ -19,5 +19,7 @@
 ! Current revisions:
 ! -----------------
-! 
+! Bugfix: division by zero could occur when calculating rib at low wind speeds
+! Bugfix: calculation of uv_total for neutral = .T., initial value for ol for
+! neutral = .T.
 ! 
 ! Former revisions:
@@ -96,5 +98,4 @@
 !>    that this method cannot be used in case of roughness heterogeneity
 !>
-!> @todo limiting of uv_total might be required in some cases
 !> @todo (re)move large_scale_forcing actions
 !> @todo check/optimize OpenMP and OpenACC directives
@@ -208,4 +209,6 @@
           CALL calc_scaling_parameters
 
+          CALL calc_uv_total
+
           IF ( .NOT. neutral )  THEN
              CALL calc_ol
@@ -220,4 +223,6 @@
 !--    number to calculate the Obukhov length 
        ELSEIF ( most_method == 'newton' .OR. most_method == 'lookup' )  THEN
+
+          CALL calc_uv_total
 
           IF ( .NOT. neutral )  THEN
@@ -273,4 +278,10 @@
 !--    Allocate field for storing the horizontal velocity
        ALLOCATE ( uv_total(nysg:nyng,nxlg:nxrg) )
+
+
+!
+!--    In case of runs with neutral statification, set Obukhov length to a
+!--    large value
+       IF ( neutral ) ol = 1.0E10_wp
 
        IF ( most_method == 'lookup' )  THEN
@@ -384,4 +395,43 @@
 ! Description:
 ! ------------
+!> Compute the absolute value of the horizontal velocity (relative to the
+!> surface). This is required by all methods
+!------------------------------------------------------------------------------!
+    SUBROUTINE calc_uv_total
+
+       IMPLICIT NONE
+
+
+       !$OMP PARALLEL DO PRIVATE( k, z_mo )
+       !$acc kernels loop
+       DO  i = nxl, nxr
+          DO  j = nys, nyn
+
+             k   = nzb_s_inner(j,i)
+             uv_total(j,i) = SQRT( ( 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1)      &
+                                         - u(k,j,i)   - u(k,j,i+1) ) )**2 +    &
+                              ( 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i)           &
+                                         - v(k,j,i)   - v(k,j+1,i) ) )**2 )
+
+!
+!--          For too small values of the local wind, MOST does not work. A
+!--          threshold value is thus set if required
+!            uv_total(j,i) = MAX(0.01_wp,uv_total(j,i))
+
+          ENDDO
+       ENDDO
+
+!
+!--    Values of uv_total need to be exchanged at the ghost boundaries
+       !$acc update host( uv_total )
+       CALL exchange_horiz_2d( uv_total )
+       !$acc update device( uv_total )
+
+    END SUBROUTINE calc_uv_total
+
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
 !> Calculate the Obukhov length (L) and Richardson flux number (z/L)
 !------------------------------------------------------------------------------!
@@ -394,40 +444,12 @@
 
        REAL(wp), DIMENSION(nysg:nyng,nxlg:nxrg) :: rib !< Bulk Richardson number
-                       
+
        REAL(wp)     :: f,      & !< Function for Newton iteration: f = Ri - [...]/[...]^2 = 0
                        f_d_ol, & !< Derivative of f
                        ol_l,   & !< Lower bound of L for Newton iteration
                        ol_m,   & !< Previous value of L for Newton iteration
-                       ol_old, & !< Previous time step value of L 
+                       ol_old, & !< Previous time step value of L
                        ol_u      !< Upper bound of L for Newton iteration
 
-!
-!--    Compute the absolute value of the horizontal velocity (relative to the 
-!--    surface). This is required by all methods
-       !$OMP PARALLEL DO PRIVATE( k, z_mo )
-       !$acc kernels loop
-       DO  i = nxl, nxr
-          DO  j = nys, nyn
-
-             k   = nzb_s_inner(j,i)
-
-             uv_total(j,i) = SQRT( ( 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1)      &
-                                         - u(k,j,i)   - u(k,j,i+1) ) )**2 +    &
-                              ( 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i)           &
-                                         - v(k,j,i)   - v(k,j+1,i) ) )**2 )
-
-!
-!--          For too small values of the local wind, MOST does not work. A
-!--          threshold value is thus set if required
-!            uv_total(j,i) = MAX(0.01_wp,uv_total(j,i))
-
-          ENDDO
-       ENDDO
-
-!
-!--    Values of uv_total need to be exchanged at the ghost boundaries
-       !$acc update host( uv_total )
-       CALL exchange_horiz_2d( uv_total )
-       !$acc update device( uv_total )
 
        IF ( TRIM( most_method ) /= 'circular' )  THEN
@@ -447,8 +469,8 @@
                    IF ( humidity )  THEN
                       rib(j,i) = g * z_mo * ( vpt(k+1,j,i) - vpt(k,j,i) )      &
-                           / ( uv_total(j,i)**2 * vpt(k+1,j,i) )
+                           / ( uv_total(j,i)**2 * vpt(k+1,j,i) + 1.0E-20_wp )
                    ELSE
                       rib(j,i) = g * z_mo * ( pt(k+1,j,i) - pt(k,j,i) )        &
-                           / ( uv_total(j,i)**2 * pt(k+1,j,i) )
+                           / ( uv_total(j,i)**2 * pt(k+1,j,i)  + 1.0E-20_wp )
                    ENDIF      
                 ELSE
@@ -462,8 +484,10 @@
                                  * q(k+1,j,i) ) * shf(j,i) + 0.61_wp           &
                                  * pt(k+1,j,i) * qsws(j,i) )   &
-                                 / ( uv_total(j,i)**3 * vpt(k+1,j,i) * kappa**2 )
+                                 / ( uv_total(j,i)**3 * vpt(k+1,j,i) * kappa**2&
+                                 + 1.0E-20_wp)
                    ELSE
                       rib(j,i) = - g * z_mo * shf(j,i)                         &
-                           / ( uv_total(j,i)**3 * pt(k+1,j,i) * kappa**2 )
+                           / ( uv_total(j,i)**3 * pt(k+1,j,i) * kappa**2       &
+                           + 1.0E-20_wp )
                    ENDIF
                 ENDIF  
