Index: palm/trunk/DOC/palm2doxygen.config
===================================================================
--- palm/trunk/DOC/palm2doxygen.config	(revision 1708)
+++ palm/trunk/DOC/palm2doxygen.config	(revision 1709)
@@ -812,5 +812,5 @@
 # exclude all test directories for example use the pattern */test/*
 
-EXCLUDE_PATTERNS       = */message.f90 */cpulog.f90
+EXCLUDE_PATTERNS       = */message.f90 */cpulog.f90  */tridia_solver.f90
 
 # The EXCLUDE_SYMBOLS tag can be used to specify one or more symbol names
Index: palm/trunk/SOURCE/flow_statistics.f90
===================================================================
--- palm/trunk/SOURCE/flow_statistics.f90	(revision 1708)
+++ palm/trunk/SOURCE/flow_statistics.f90	(revision 1709)
@@ -19,5 +19,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! Updated output of Obukhov length
 ! 
 ! Former revisions:
@@ -1554,16 +1554,9 @@
        ts_value(20,sr) = hom(nzb+2,1,pr_palm,sr)    ! v'w'    at k=0
        ts_value(21,sr) = hom(nzb,1,48,sr)           ! w"q"    at k=0
- 
-!
-!--    Calculate obukhov length
-       IF ( ts_value(5,sr) /= 0.0_wp )  THEN
-!           IF ( most_method == 'circular' )  THEN
-!              ts_value(22,sr) = ts_value(4,sr)**2 /                             &
-!                           ( kappa * g * ts_value(5,sr) / ts_value(18,sr) ) 
-!           ELSE
-             ts_value(22,sr) = hom(nzb,1,114,sr)  
-!          ENDIF     
+
+       IF ( .NOT. neutral )  THEN
+          ts_value(22,sr) = hom(nzb,1,114,sr)          ! L
        ELSE
-          ts_value(22,sr) = 10000.0_wp
+          ts_value(22,sr) = 1.0E10_wp
        ENDIF
 
@@ -3503,8 +3496,8 @@
        ts_value(21,sr) = hom(nzb,1,48,sr)           ! w"q"    at k=0
 
-       IF ( ts_value(5,sr) /= 0.0_wp )  THEN
-          ts_value(22,sr) = hom(nzb,1,114,sr)
+       IF ( .NOT. neutral )  THEN
+          ts_value(22,sr) = hom(nzb,1,114,sr)          ! L
        ELSE
-          ts_value(22,sr) = 10000.0_wp
+          ts_value(22,sr) = 1.0E10_wp
        ENDIF
 
Index: palm/trunk/SOURCE/init_1d_model.f90
===================================================================
--- palm/trunk/SOURCE/init_1d_model.f90	(revision 1708)
+++ palm/trunk/SOURCE/init_1d_model.f90	(revision 1709)
@@ -19,5 +19,6 @@
 ! Current revisions:
 ! -----------------
-! 
+! Set initial time step to 10 s to avoid instability of the 1d model for small
+! grid spacings
 !
 ! Former revisions:
@@ -81,4 +82,5 @@
 !>
 !> @todo harmonize code with new surface_layer_fluxes module
+!> @bug 1D model crashes when using small grid spacings in the order of 1 m
 !------------------------------------------------------------------------------!
  SUBROUTINE init_1d_model
@@ -262,5 +264,5 @@
 !-- Determine the time step at the start of a 1D-simulation and
 !-- determine and printout quantities used for run control
-    CALL timestep_1d
+    dt_1d = 10.0_wp
     CALL run_control_1d
 
@@ -864,5 +866,5 @@
     USE pegrid
     
-    USE control_parameters,                                                              &
+    USE control_parameters,                                                    &
         ONLY:  message_string
 
Index: palm/trunk/SOURCE/land_surface_model.f90
===================================================================
--- palm/trunk/SOURCE/land_surface_model.f90	(revision 1708)
+++ palm/trunk/SOURCE/land_surface_model.f90	(revision 1709)
@@ -19,5 +19,9 @@
 ! Current revisions:
 ! -----------------
-! 
+! Renamed pt_1 and qv_1 to pt1 and qv1.
+! Bugfix: set initial values for t_surface_p in case of restart runs
+! Bugfix: zero resistance caused crash when using radiation_scheme = 'clear-sky'
+! Bugfix: calculation of rad_net when using radiation_scheme = 'clear-sky'
+! Added todo action
 !
 ! Former revisions:
@@ -87,5 +91,5 @@
 !>
 !> @todo Consider partial absorption of the net shortwave radiation by the 
-!>       surface layer.
+!>       skin layer.
 !> @todo Allow for water surfaces
 !> @todo Invert indices (running from -3 to 0. Currently: nzb_soil=0, 
@@ -93,7 +97,8 @@
 !> @todo Implement surface runoff model (required when performing long-term LES 
 !>       with considerable precipitation.
+!> @todo Fix crashes with radiation_scheme == 'constant'
 !>
-!> @note  No time step criterion is required as long as the soil layers do not 
-!>        become too thin.
+!> @note No time step criterion is required as long as the soil layers do not
+!>       become too thin.
 !------------------------------------------------------------------------------!
  MODULE land_surface_model_mod
@@ -124,9 +129,9 @@
     USE radiation_model_mod,                                                   &
         ONLY:  force_radiation_call, radiation_scheme, rad_net, rad_sw_in,     &
-               rad_lw_out, sigma_sb
+               rad_lw_out, rad_lw_out_change_0, sigma_sb
         
 #if defined ( __rrtmg )
     USE radiation_model_mod,                                                   &
-        ONLY:  rrtm_lwuflx_dt, rrtm_idrv
+        ONLY:  rrtm_idrv
 #endif
 
@@ -275,5 +280,5 @@
               rad_net_l,        & !< local copy of rad_net (net radiation at surface)
               r_a,              & !< aerodynamic resistance 
-              r_a_av,           & !< avergae of r_a
+              r_a_av,           & !< average of r_a
               r_canopy,         & !< canopy resistance
               r_soil,           & !< soil resistance
@@ -635,5 +640,5 @@
        INTEGER(iwp) ::  k !< running index
 
-       REAL(wp) :: pt_1   !< potential temperature at first grid level
+       REAL(wp) :: pt1   !< potential temperature at first grid level
 
 
@@ -788,5 +793,4 @@
 !--       Calculate surface temperature
           t_surface   = pt_surface * exn
-          t_surface_p = t_surface
 
 !
@@ -798,15 +802,15 @@
 
                 IF ( cloud_physics )  THEN
-                   pt_1 = pt(k+1,j,i) + l_d_cp * pt_d_t(k+1) * ql(k+1,j,i)
+                   pt1 = pt(k+1,j,i) + l_d_cp * pt_d_t(k+1) * ql(k+1,j,i)
                 ELSE
-                   pt_1 = pt(k+1,j,i)
+                   pt1 = pt(k+1,j,i)
                 ENDIF
 
 !
 !--             Assure that r_a cannot be zero at model start
-                IF ( pt_1 == pt(k,j,i) )  pt_1 = pt_1 + 1.0E-10_wp
+                IF ( pt1 == pt(k,j,i) )  pt1 = pt1 + 1.0E-10_wp
 
                 us(j,i)  = 0.1_wp
-                ts(j,i)  = (pt_1 - pt(k,j,i)) / r_a(j,i)
+                ts(j,i)  = (pt1 - pt(k,j,i)) / r_a(j,i)
                 shf(j,i) = - us(j,i) * ts(j,i)
              ENDDO
@@ -899,4 +903,7 @@
        m_soil_p    = m_soil
        m_liq_eb_p  = m_liq_eb
+       t_surface_p = t_surface
+
+
 
 !--    Store initial profiles of t_soil and m_soil (assuming they are 
@@ -926,5 +933,4 @@
        dots_unit(dots_soil+6:dots_soil+7) = "s/m"
 
-       RETURN
 
     END SUBROUTINE init_lsm
@@ -964,6 +970,6 @@
                    lambda_surface, & !< Current value of lambda_surface
                    m_liq_eb_max,   & !< maxmimum value of the liq. water reservoir
-                   pt_1,        & !< potential temperature at first grid level
-                   qv_1           !< specific humidity at first grid level
+                   pt1,         & !< potential temperature at first grid level
+                   qv1            !< specific humidity at first grid level
 
 !
@@ -989,12 +995,16 @@
 !--          equivalent to the ECMWF formulation using drag coefficients
              IF ( cloud_physics )  THEN
-                pt_1 = pt(k+1,j,i) + l_d_cp * pt_d_t(k+1) * ql(k+1,j,i)
-                qv_1 = q(k+1,j,i) - ql(k+1,j,i)
+                pt1 = pt(k+1,j,i) + l_d_cp * pt_d_t(k+1) * ql(k+1,j,i)
+                qv1 = q(k+1,j,i) - ql(k+1,j,i)
              ELSE
-                pt_1 = pt(k+1,j,i)
-                qv_1 = q(k+1,j,i)
-             ENDIF
-
-             r_a(j,i) = (pt_1 - pt(k,j,i)) / (ts(j,i) * us(j,i) + 1.0E-20_wp)
+                pt1 = pt(k+1,j,i)
+                qv1 = q(k+1,j,i)
+             ENDIF
+
+             r_a(j,i) = (pt1 - pt(k,j,i)) / (ts(j,i) * us(j,i) + 1.0E-20_wp)
+
+!
+!--          Make sure that the resistance does not drop to zero
+             IF ( ABS(r_a(j,i)) < 1.0E-10_wp )  r_a(j,i) = 1.0E-10_wp
 
 !
@@ -1005,10 +1015,11 @@
 !--          night)
              IF ( radiation_scheme /= 'constant' )  THEN
-                f1 = MIN( 1.0_wp, ( 0.004_wp * rad_sw_in(k,j,i) + 0.05_wp ) /&
-                              (0.81_wp * (0.004_wp * rad_sw_in(k,j,i)        &
+                f1 = MIN( 1.0_wp, ( 0.004_wp * rad_sw_in(k,j,i) + 0.05_wp ) /  &
+                              (0.81_wp * (0.004_wp * rad_sw_in(k,j,i)          &
                                + 1.0_wp)) )
              ELSE
                 f1 = 1.0_wp
              ENDIF
+
 
 !
@@ -1040,5 +1051,5 @@
 !
 !--             Calculate vapour pressure
-                e  = qv_1 * surface_pressure / 0.622_wp
+                e  = qv1 * surface_pressure / 0.622_wp
                 f3 = EXP ( -g_d(j,i) * (e_s - e) )
              ELSE
@@ -1082,5 +1093,5 @@
 !--          In case of dewfall, set evapotranspiration to zero
 !--          All super-saturated water is then removed from the air
-             IF ( humidity .AND. q_s <= qv_1 )  THEN
+             IF ( humidity .AND. q_s <= qv1 )  THEN
                 r_canopy(j,i) = 0.0_wp
                 r_soil(j,i)   = 0.0_wp
@@ -1108,12 +1119,8 @@
 !--          reservoir are not allowed to take up water from the super-saturated
 !--          air.
-             IF ( humidity )  THEN
-                IF ( q_s <= qv_1 )  THEN
-                   IF ( .NOT. dewfall )  THEN
-                      f_qsws_veg  = 0.0_wp
-                      f_qsws_soil = 0.0_wp
-                      f_qsws_liq  = 0.0_wp
-                   ENDIF
-                ENDIF
+             IF ( humidity .AND. q_s <= qv1 .AND. .NOT. dewfall )  THEN
+                f_qsws_veg  = 0.0_wp
+                f_qsws_soil = 0.0_wp
+                f_qsws_liq  = 0.0_wp
              ENDIF 
 
@@ -1133,12 +1140,11 @@
              rad_net_l(j,i) = rad_net(j,i) + rad_lw_out(nzb,j,i)
 
-
              IF ( humidity )  THEN
 #if defined ( __rrtmg )
 !
 !--             Numerator of the prognostic equation
-                coef_1 = rad_net_l(j,i) + rrtm_lwuflx_dt(0,nzb)                &
+                coef_1 = rad_net_l(j,i) + rad_lw_out_change_0(j,i)             &
                          * t_surface(j,i) - rad_lw_out(nzb,j,i)                &
-                         + f_shf * pt_1 + f_qsws * ( qv_1 - q_s                &
+                         + f_shf * pt1 + f_qsws * ( qv1 - q_s                  &
                          + dq_s_dt * t_surface(j,i) ) + lambda_surface         &
                          * t_soil(nzb_soil,j,i)
@@ -1146,5 +1152,5 @@
 !
 !--             Denominator of the prognostic equation
-                coef_2 = rrtm_lwuflx_dt(0,nzb) + f_qsws * dq_s_dt              &
+                coef_2 = rad_lw_out_change_0(j,i) + f_qsws * dq_s_dt           &
                          + lambda_surface + f_shf / exn
 #else
@@ -1154,12 +1160,12 @@
                 coef_1 = rad_net_l(j,i) + 3.0_wp * sigma_sb                    &
                          * t_surface(j,i) ** 4                                 &
-                         + f_shf * pt_1 + f_qsws * ( qv_1                      &
+                         + f_shf * pt1 + f_qsws * ( qv1                        &
                          - q_s + dq_s_dt * t_surface(j,i) )                    &
                          + lambda_surface * t_soil(nzb_soil,j,i)
 
 !
-!--                Denominator of the prognostic equation
-                   coef_2 = 4.0_wp * sigma_sb * t_surface(j,i) ** 3 + f_qsws   &
-                            * dq_s_dt + lambda_surface + f_shf / exn
+!--             Denominator of the prognostic equation
+                coef_2 = 4.0_wp * sigma_sb * t_surface(j,i) ** 3 + f_qsws      &
+                         * dq_s_dt + lambda_surface + f_shf / exn
  
 #endif
@@ -1169,12 +1175,12 @@
 !
 !--             Numerator of the prognostic equation
-                coef_1 = rad_net_l(j,i) + rrtm_lwuflx_dt(0,nzb)                &
+                coef_1 = rad_net_l(j,i) + rad_lw_out_change_0(j,i)             &
                          * t_surface(j,i) - rad_lw_out(nzb,j,i)                &
-                         + f_shf * pt_1  + lambda_surface                      &
+                         + f_shf * pt1  + lambda_surface                       &
                          * t_soil(nzb_soil,j,i)
 
 !
 !--             Denominator of the prognostic equation
-                coef_2 = rrtm_lwuflx_dt(0,nzb) + lambda_surface + f_shf / exn
+                coef_2 = rad_lw_out_change_0(j,i) + lambda_surface + f_shf / exn
 #else
 
@@ -1182,5 +1188,5 @@
 !--             Numerator of the prognostic equation
                 coef_1 = rad_net_l(j,i) + 3.0_wp * sigma_sb                    &
-                          * t_surface(j,i) ** 4 + f_shf * pt_1                 &
+                          * t_surface(j,i) ** 4 + f_shf * pt1                  &
                           + lambda_surface * t_soil(nzb_soil,j,i)
 
@@ -1191,5 +1197,4 @@
 #endif
              ENDIF
-
 
              tend = 0.0_wp
@@ -1238,12 +1243,12 @@
 !--          Calculate fluxes
 #if defined ( __rrtmg )
-             rad_net_l(j,i)   = rad_net_l(j,i) + rrtm_lwuflx_dt(0,nzb)         &
+             rad_net_l(j,i)   = rad_net_l(j,i) + rad_lw_out_change_0(j,i)      &
                                 * t_surface(j,i) - rad_lw_out(nzb,j,i)         &
-                                - rrtm_lwuflx_dt(0,nzb) * t_surface_p(j,i)
+                                - rad_lw_out_change_0(j,i) * t_surface_p(j,i)
 
              IF ( rrtm_idrv == 1 )  THEN
                 rad_net(j,i) = rad_net_l(j,i)
                 rad_lw_out(nzb,j,i) = rad_lw_out(nzb,j,i)                      &
-                                      + rrtm_lwuflx_dt(0,nzb)                  &
+                                      + rad_lw_out_change_0(j,i)               &
                                       * ( t_surface_p(j,i) - t_surface(j,i) )
              ENDIF
@@ -1252,29 +1257,29 @@
                                 * t_surface(j,i)**4 - 4.0_wp * sigma_sb        &
                                 * t_surface(j,i)**3 * t_surface_p(j,i)
-             rad_net(j,i) = rad_net_l(j,i)
 #endif
+
 
              ghf_eb(j,i)    = lambda_surface * (t_surface_p(j,i)               &
                               - t_soil(nzb_soil,j,i))
 
-             shf_eb(j,i)    = - f_shf * ( pt_1 - pt(k,j,i) )
+             shf_eb(j,i)    = - f_shf * ( pt1 - pt(k,j,i) )
 
              shf(j,i) = shf_eb(j,i) / rho_cp
 
              IF ( humidity )  THEN
-                qsws_eb(j,i)  = - f_qsws    * ( qv_1 - q_s + dq_s_dt           &
+                qsws_eb(j,i)  = - f_qsws    * ( qv1 - q_s + dq_s_dt            &
                                 * t_surface(j,i) - dq_s_dt * t_surface_p(j,i) )
 
                 qsws(j,i) = qsws_eb(j,i) / rho_lv
 
-                qsws_veg_eb(j,i)  = - f_qsws_veg  * ( qv_1 - q_s               &
+                qsws_veg_eb(j,i)  = - f_qsws_veg  * ( qv1 - q_s                &
                                     + dq_s_dt * t_surface(j,i) - dq_s_dt       &
                                     * t_surface_p(j,i) )
 
-                qsws_soil_eb(j,i) = - f_qsws_soil * ( qv_1 - q_s               &
+                qsws_soil_eb(j,i) = - f_qsws_soil * ( qv1 - q_s                &
                                     + dq_s_dt * t_surface(j,i) - dq_s_dt       &
                                     * t_surface_p(j,i) )
 
-                qsws_liq_eb(j,i)  = - f_qsws_liq  * ( qv_1 - q_s               &
+                qsws_liq_eb(j,i)  = - f_qsws_liq  * ( qv1 - q_s                &
                                     + dq_s_dt * t_surface(j,i) - dq_s_dt       &
                                     * t_surface_p(j,i) )
@@ -1285,5 +1290,5 @@
                 r_s(j,i) = 1.0E10_wp
              ELSE
-                r_s(j,i) = - rho_lv * ( qv_1 - q_s + dq_s_dt                   &
+                r_s(j,i) = - rho_lv * ( qv1 - q_s + dq_s_dt                   &
                            * t_surface(j,i) - dq_s_dt * t_surface_p(j,i) )     &
                            / qsws_eb(j,i) - r_a(j,i)
Index: palm/trunk/SOURCE/radiation_model.f90
===================================================================
--- palm/trunk/SOURCE/radiation_model.f90	(revision 1708)
+++ palm/trunk/SOURCE/radiation_model.f90	(revision 1709)
@@ -19,5 +19,6 @@
 ! Current revisions:
 ! -----------------
-! 
+! Bugfix: set initial value for rrtm_lwuflx_dt to zero, small formatting
+! corrections
 ! 
 ! Former revisions:
@@ -197,4 +198,5 @@
     REAL(wp), DIMENSION(:,:), ALLOCATABLE :: &
                 alpha,                       & !< surface broadband albedo (used for clear-sky scheme)
+                rad_lw_out_change_0,         & !< change in LW out due to change in surface temperature
                 rad_net,                     & !< net radiation at the surface
                 rad_net_av                     !< average of rad_net
@@ -363,13 +365,13 @@
            radiation_clearsky, radiation_rrtmg, radiation_scheme,              &
            radiation_tendency, rad_lw_in, rad_lw_in_av, rad_lw_out,            &
-           rad_lw_out_av, rad_lw_cs_hr, rad_lw_cs_hr_av, rad_lw_hr,            &
-           rad_lw_hr_av, rad_sw_in, rad_sw_in_av, rad_sw_out, rad_sw_out_av,   &
-           rad_sw_cs_hr, rad_sw_cs_hr_av, rad_sw_hr, rad_sw_hr_av, sigma_sb,   &
-           skip_time_do_radiation, sw_radiation, time_radiation, time_utc_init
+           rad_lw_out_av, rad_lw_out_change_0, rad_lw_cs_hr, rad_lw_cs_hr_av,  &
+           rad_lw_hr, rad_lw_hr_av, rad_sw_in, rad_sw_in_av, rad_sw_out,       &
+           rad_sw_out_av, rad_sw_cs_hr, rad_sw_cs_hr_av, rad_sw_hr,            &
+           rad_sw_hr_av, sigma_sb, skip_time_do_radiation, sw_radiation,       &
+           time_radiation, time_utc_init
 
 
 #if defined ( __rrtmg )
-    PUBLIC rrtm_aldif, rrtm_aldir, rrtm_asdif, rrtm_asdir, rrtm_idrv,          &
-           rrtm_lwuflx_dt
+    PUBLIC rrtm_aldif, rrtm_aldir, rrtm_asdif, rrtm_asdir, rrtm_idrv
 #endif
 
@@ -390,4 +392,11 @@
           ALLOCATE ( rad_net(nysg:nyng,nxlg:nxrg) )
           rad_net = 0.0_wp
+       ENDIF
+
+!
+!--    Allocate array for storing the surface net radiation
+       IF ( .NOT. ALLOCATED ( rad_lw_out_change_0 ) )  THEN
+          ALLOCATE ( rad_lw_out_change_0(nysg:nyng,nxlg:nxrg) )
+          rad_lw_out_change_0 = 0.0_wp
        ENDIF
 
@@ -654,5 +663,5 @@
        IF (  TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
           IF ( radiation_scheme == 'clear-sky' )  THEN
-              CALL radiation_clearsky
+             CALL radiation_clearsky
           ELSEIF ( radiation_scheme == 'rrtmg' )  THEN
              CALL radiation_rrtmg
@@ -679,6 +688,6 @@
        INTEGER(iwp) :: i, j, k   !< loop indices
        REAL(wp)     :: exn,   &  !< Exner functions at surface
-                       exn_1, &  !< Exner functions at first grid level
-                       pt_1      !< potential temperature at first grid level
+                       exn1,  &  !< Exner functions at first grid level
+                       pt1       !< potential temperature at first grid level
 
 !
@@ -696,9 +705,9 @@
 !--    Calculate radiation fluxes and net radiation (rad_net) for each grid 
 !--    point
-       DO i = nxl, nxr
-          DO j = nys, nyn
+       DO i = nxlg, nxrg
+          DO j = nysg, nyng
              k = nzb_s_inner(j,i)
 
-             exn_1 = (hyp(k+1) / 100000.0_wp )**0.286_wp
+             exn1 = (hyp(k+1) / 100000.0_wp )**0.286_wp
 
              rad_sw_in(0,j,i)  = solar_constant * sky_trans * zenith(0)
@@ -707,8 +716,8 @@
 
              IF ( cloud_physics )  THEN
-                pt_1 = pt(k+1,j,i) + l_d_cp / exn_1 * ql(k+1,j,i)
-                rad_lw_in(0,j,i)  = 0.8_wp * sigma_sb * (pt_1 * exn_1)**4
+                pt1 = pt(k+1,j,i) + l_d_cp / exn1 * ql(k+1,j,i)
+                rad_lw_in(0,j,i)  = 0.8_wp * sigma_sb * (pt1 * exn1)**4
              ELSE
-                rad_lw_in(0,j,i)  = 0.8_wp * sigma_sb * (pt(k+1,j,i) * exn_1)**4
+                rad_lw_in(0,j,i)  = 0.8_wp * sigma_sb * (pt(k+1,j,i) * exn1)**4
              ENDIF
 
@@ -718,12 +727,4 @@
           ENDDO
        ENDDO
-
-       CALL exchange_horiz_2d( rad_lw_in,  nbgp )
-       CALL exchange_horiz_2d( rad_lw_out, nbgp )
-       CALL exchange_horiz_2d( rad_sw_in,  nbgp )
-       CALL exchange_horiz_2d( rad_sw_out, nbgp )
-       CALL exchange_horiz_2d( rad_net,    nbgp )
-
-       RETURN
 
     END SUBROUTINE radiation_clearsky
@@ -901,4 +902,8 @@
                    rad_lw_cs_hr(k,j,i)  = rrtm_lwhrc(0,k) * d_hours_day
                 ENDDO
+
+!
+!--             Save change in LW heating rate
+                rad_lw_out_change_0(j,i) = rrtm_lwuflx_dt(0,nzb)
 
              ENDIF
@@ -953,4 +958,5 @@
 
        CALL exchange_horiz_2d( rad_net, nbgp )
+       CALL exchange_horiz_2d( rad_lw_out_change_0, nbgp )
 #endif
 
@@ -1347,5 +1353,5 @@
        ALLOCATE ( rrtm_lwuflxc_dt(0:0,nzb:nzt_rad+1) )
 
-       rrtm_lwuflxc_dt = 0.0_wp
+       rrtm_lwuflx_dt = 0.0_wp
        rrtm_lwuflxc_dt = 0.0_wp
 
Index: palm/trunk/SOURCE/read_3d_binary.f90
===================================================================
--- palm/trunk/SOURCE/read_3d_binary.f90	(revision 1708)
+++ palm/trunk/SOURCE/read_3d_binary.f90	(revision 1709)
@@ -19,5 +19,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! Added rad_lw_out_change_0, increased binary_version
 ! 
 ! Former revisions:
@@ -124,7 +124,7 @@
         ONLY: rad_net, rad_net_av, rad_lw_in, rad_lw_in_av, rad_lw_out,        &
               rad_lw_cs_hr, rad_lw_cs_hr_av, rad_lw_hr, rad_lw_hr_av,          &
-              rad_lw_out_av, rad_sw_in, rad_sw_in_av, rad_sw_out,              &
-              rad_sw_out_av, rad_sw_cs_hr, rad_sw_cs_hr_av, rad_sw_hr,         &
-              rad_sw_hr_av
+              rad_lw_out_av, rad_lw_out_change_0, rad_sw_in, rad_sw_in_av,     &
+              rad_sw_out, rad_sw_out_av, rad_sw_cs_hr, rad_sw_cs_hr_av,        &
+              rad_sw_hr, rad_sw_hr_av
 
     USE random_function_mod,                                                   &
@@ -348,5 +348,5 @@
 !--    First compare the version numbers
        READ ( 13 )  version_on_file
-       binary_version = '4.1'
+       binary_version = '4.2'
        IF ( TRIM( version_on_file ) /= TRIM( binary_version ) )  THEN
           WRITE( message_string, * ) 'version mismatch concerning data ',      &
@@ -877,4 +877,12 @@
                    rad_lw_out_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
                              tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+
+                CASE ( 'rad_lw_out_change_0' )
+                   IF ( .NOT. ALLOCATED( rad_lw_out_change_0 ) )  THEN
+                      ALLOCATE( rad_lw_out_change_0(nysg:nyng,nxlg:nxrg) )
+                   ENDIF
+                   IF ( k == 1 )  READ ( 13 )  tmp_2d
+                   rad_lw_out_change_0(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)&
+                              = tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
 
                 CASE ( 'rad_lw_cs_hr' )
Index: palm/trunk/SOURCE/surface_layer_fluxes.f90
===================================================================
--- palm/trunk/SOURCE/surface_layer_fluxes.f90	(revision 1708)
+++ 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  
Index: palm/trunk/SOURCE/write_3d_binary.f90
===================================================================
--- palm/trunk/SOURCE/write_3d_binary.f90	(revision 1708)
+++ palm/trunk/SOURCE/write_3d_binary.f90	(revision 1709)
@@ -19,5 +19,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! Added rad_lw_out_change_0, increased binary_version
 ! 
 ! Former revisions:
@@ -114,8 +114,8 @@
     USE radiation_model_mod,                                                   &
         ONLY: radiation, rad_net, rad_net_av, rad_lw_in, rad_lw_in_av,         &
-              rad_lw_out, rad_lw_out_av, rad_lw_cs_hr, rad_lw_cs_hr_av,        &
-              rad_lw_hr, rad_lw_hr_av, rad_sw_in, rad_sw_in_av, rad_sw_out,    &
-              rad_sw_out_av, rad_sw_cs_hr, rad_sw_cs_hr_av, rad_sw_hr,         &
-              rad_sw_hr_av
+              rad_lw_out, rad_lw_out_av, rad_lw_out_change_0, rad_lw_cs_hr,   &
+              rad_lw_cs_hr_av, rad_lw_hr, rad_lw_hr_av, rad_sw_in,             &
+              rad_sw_in_av, rad_sw_out, rad_sw_out_av, rad_sw_cs_hr,           &
+              rad_sw_cs_hr_av, rad_sw_hr, rad_sw_hr_av
 
     USE random_function_mod,                                                   &
@@ -139,5 +139,5 @@
 !
 !-- Write arrays.
-    binary_version = '4.1'
+    binary_version = '4.2'
 
     WRITE ( 14 )  binary_version
@@ -297,4 +297,8 @@
           WRITE ( 14 )  'rad_lw_out_av       ';  WRITE ( 14 )  rad_lw_out_av  
        ENDIF 
+       IF ( ALLOCATED( rad_lw_out_change_0 ) )  THEN
+          WRITE ( 14 )  'rad_lw_out_change_0 '
+          WRITE ( 14 )  rad_lw_out_change_0
+       ENDIF
        IF ( ALLOCATED( rad_lw_cs_hr ) )  THEN
           WRITE ( 14 )  'rad_lw_cs_hr        ';  WRITE ( 14 )  rad_lw_cs_hr
