Index: palm/trunk/SOURCE/check_parameters.f90
===================================================================
--- palm/trunk/SOURCE/check_parameters.f90	(revision 3667)
+++ palm/trunk/SOURCE/check_parameters.f90	(revision 3668)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! most_method removed
+! 
+! 3655 2019-01-07 16:51:22Z knoop
 ! Formatting
 ! 
@@ -3790,13 +3793,4 @@
     ENDIF
 
-!
-!-- Check for valid setting of most_method
-    IF ( TRIM( most_method ) /= 'circular'  .AND.                              &
-         TRIM( most_method ) /= 'newton'    .AND.                              &
-         TRIM( most_method ) /= 'lookup' )  THEN
-       message_string = 'most_method = "' // TRIM( most_method ) //            &
-                        '" is unknown'
-       CALL message( 'check_parameters', 'PA0416', 1, 2, 0, 6, 0 )
-    ENDIF
 
 !
Index: palm/trunk/SOURCE/land_surface_model_mod.f90
===================================================================
--- palm/trunk/SOURCE/land_surface_model_mod.f90	(revision 3667)
+++ palm/trunk/SOURCE/land_surface_model_mod.f90	(revision 3668)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Removed most_method
+! 
+! 3655 2019-01-07 16:51:22Z knoop
 ! nopointer option removed
 ! 
@@ -1324,6 +1327,5 @@
 
     USE control_parameters,                                                    &
-        ONLY:  bc_pt_b, bc_q_b, constant_flux_layer, message_string,           &
-               most_method
+        ONLY:  bc_pt_b, bc_q_b, constant_flux_layer, message_string
                      
     
@@ -1446,11 +1448,4 @@
     
     IF ( TRIM( surface_type ) == 'water' )  THEN
-
-       IF ( TRIM( most_method ) == 'lookup' )  THEN    
-          WRITE( message_string, * ) 'surface_type = ', surface_type,          &
-                                     ' is not allowed in combination with ',   &
-                                     'most_method = ', most_method
-          CALL message( 'lsm_check_parameters', 'PA0414', 1, 2, 0, 6, 0 )
-       ENDIF
 
        IF ( water_type == 0 )  THEN  
@@ -1531,11 +1526,4 @@
 
     IF ( TRIM( surface_type ) == 'netcdf' )  THEN
-       IF ( ANY( water_type_f%var /= water_type_f%fill )  .AND.                &
-            TRIM( most_method ) == 'lookup' )  THEN    
-          WRITE( message_string, * ) 'water-surfaces are not allowed in ' //   &
-                                     'combination with most_method = ',        &
-                                     TRIM( most_method )
-          CALL message( 'lsm_check_parameters', 'PA0999', 2, 2, 0, 6, 0 )
-       ENDIF
 !
 !--    MS: Some problme here, after calling message everythings stucks at
Index: palm/trunk/SOURCE/modules.f90
===================================================================
--- palm/trunk/SOURCE/modules.f90	(revision 3667)
+++ palm/trunk/SOURCE/modules.f90	(revision 3668)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Removed most_method
+! 
+! 3648 2019-01-02 16:35:46Z suehring
 ! -surface_data_output +surface_output
 ! 
@@ -1090,5 +1093,4 @@
     CHARACTER (LEN=8)    ::  coupling_char = ''                           !< appended to filenames in coupled or nested runs ('_O': ocean PE,
                                                                           !< '_NV': vertically nested atmosphere PE, '_N##': PE of nested domain ##
-    CHARACTER (LEN=8)    ::  most_method = 'newton'                       !< namelist parameter
     CHARACTER (LEN=10)   ::  run_date = ' '                               !< date of simulation run
     CHARACTER (LEN=8)    ::  run_time = ' '                               !< time of simulation run
Index: palm/trunk/SOURCE/parin.f90
===================================================================
--- palm/trunk/SOURCE/parin.f90	(revision 3667)
+++ palm/trunk/SOURCE/parin.f90	(revision 3668)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Removed most_method
+! 
+! 3649 2019-01-02 16:52:21Z suehring
 ! Delete debug-print statements
 ! 
@@ -573,5 +576,4 @@
              loop_optimization, lsf_exception, masking_method, mg_cycles,      &
              mg_switch_to_pe0_level, mixing_length_1d, momentum_advec,         &
-             most_method,                                                      &
              netcdf_precision, neutral, ngsrb,                                 &
              nsor, nsor_ini, nudging, nx, ny, nz, ocean_mode, omega,           &
@@ -644,5 +646,4 @@
              loop_optimization, lsf_exception, masking_method, mg_cycles,      &
              mg_switch_to_pe0_level, mixing_length_1d, momentum_advec,         &
-             most_method,                                                      &
              netcdf_precision, neutral, ngsrb,                                 &
              nsor, nsor_ini, nudging, nx, ny, nz, ocean_mode, omega,           &
Index: palm/trunk/SOURCE/read_restart_data_mod.f90
===================================================================
--- palm/trunk/SOURCE/read_restart_data_mod.f90	(revision 3667)
+++ palm/trunk/SOURCE/read_restart_data_mod.f90	(revision 3668)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Removed most_method and increased binary version
+! 
+! 3655 2019-01-07 16:51:22Z knoop
 ! Implementation of the PALM module interface
 ! 
@@ -203,5 +206,5 @@
        READ ( 13 )  version_on_file
 
-       binary_version_global = '4.7'
+       binary_version_global = '4.8'
        IF ( TRIM( version_on_file ) /= TRIM( binary_version_global ) )  THEN
           WRITE( message_string, * ) 'version mismatch concerning ',           &
@@ -488,6 +491,4 @@
              CASE ( 'momentum_advec' )
                 READ ( 13 )  momentum_advec
-             CASE ( 'most_method' )
-                READ ( 13 )  most_method
              CASE ( 'netcdf_precision' )
                 READ ( 13 )  netcdf_precision
Index: palm/trunk/SOURCE/surface_layer_fluxes_mod.f90
===================================================================
--- palm/trunk/SOURCE/surface_layer_fluxes_mod.f90	(revision 3667)
+++ palm/trunk/SOURCE/surface_layer_fluxes_mod.f90	(revision 3668)
@@ -26,4 +26,7 @@
 ! -----------------
 ! $Id$
+! Removed methods "circular" and "lookup"
+! 
+! 3655 2019-01-07 16:51:22Z knoop
 ! OpenACC port for SPEC
 ! 
@@ -228,11 +231,5 @@
 ! ------------
 !> Diagnostic computation of vertical fluxes in the constant flux layer from the
-!> values of the variables at grid point k=1. Three different methods are
-!> available:
-!> 1) the "old" version (most_method = 'circular') which is fast, but inaccurate
-!> 2) a Newton iteration method (most_method = 'newton'), which is accurate, but
-!>    slower
-!> 3) a method using a lookup table which is fast and accurate. Note, however,
-!>    that this method cannot be used in case of roughness heterogeneity
+!> values of the variables at grid point k=1 based on Newton iteration
 !>
 !> @todo (re)move large_scale_forcing actions
@@ -261,5 +258,5 @@
                intermediate_timestep_count, intermediate_timestep_count_max,   &
                land_surface, large_scale_forcing, lsf_surf, message_string,    &
-               most_method, neutral, passive_scalar, pt_surface, q_surface,    &
+               neutral, passive_scalar, pt_surface, q_surface,                 &
                run_coupled, surface_pressure, simulated_time, terminate_run,   &
                time_since_reference_point, urban_surface,                      &
@@ -400,106 +397,52 @@
 !--    First, calculate the new Obukhov length, then new friction velocity,
 !--    followed by the new scaling parameters (th*, q*, etc.), and the new
-!--    surface fluxes if required. The old routine ("circular") requires a 
-!--    different order of calls as the scaling parameters from the previous time
-!--    steps are used to calculate the Obukhov length
-
-!
-!--    Depending on setting of most_method use the "old" routine
-!--    Note, each routine is called for different surface types. 
+!--    surface fluxes if required. Note, each routine is called for different surface types. 
 !--    First call for default-type horizontal surfaces, for natural- and 
 !--    urban-type surfaces. Note, at this place only upward-facing horizontal
-!--    surfaces are treted. 
-       IF ( most_method == 'circular' )  THEN
-!
-!--       Default-type upward-facing horizontal surfaces
-          IF ( surf_def_h(0)%ns >= 1  )  THEN
-             surf => surf_def_h(0)
-             CALL calc_scaling_parameters
-             CALL calc_uvw_abs
-             IF ( .NOT. neutral )  CALL calc_ol
-             CALL calc_us
-             CALL calc_surface_fluxes
-             
-             IF ( do_output_at_2m )  THEN
-                CALL calc_pt_near_surface ( '2m' )
-             ENDIF
-          ENDIF
-!
-!--       Natural-type horizontal surfaces
-          IF ( surf_lsm_h%ns >= 1  )  THEN
-             surf => surf_lsm_h
-             CALL calc_scaling_parameters
-             CALL calc_uvw_abs
-             IF ( .NOT. neutral )  CALL calc_ol
-             CALL calc_us
-             CALL calc_surface_fluxes
-             IF ( do_output_at_2m )  THEN
-                CALL calc_pt_near_surface ( '2m' )
-             ENDIF
-          ENDIF
-!
-!--       Urban-type horizontal surfaces
-          IF ( surf_usm_h%ns >= 1  )  THEN
-             surf => surf_usm_h
-             CALL calc_scaling_parameters
-             CALL calc_uvw_abs
-             IF ( .NOT. neutral )  CALL calc_ol
-             CALL calc_us
-             CALL calc_surface_fluxes
-             IF ( do_output_at_2m )  THEN
-                CALL calc_pt_near_surface ( '2m' )
-             ENDIF
-             IF ( indoor_model )  THEN
-                CALL calc_pt_near_surface ( '10cm' )
-             ENDIF
-          ENDIF
-!
-!--    Use either Newton iteration or a lookup table for the bulk Richardson
-!--    number to calculate the Obukhov length 
-       ELSEIF ( most_method == 'newton'  .OR.  most_method == 'lookup' )  THEN
-!
-!--       Default-type upward-facing horizontal surfaces
-          IF ( surf_def_h(0)%ns >= 1  )  THEN
-             surf => surf_def_h(0)
-             CALL calc_uvw_abs
-             IF ( .NOT. neutral )  CALL calc_ol
-             CALL calc_us
-             CALL calc_scaling_parameters
-             CALL calc_surface_fluxes
-             IF ( do_output_at_2m )  THEN
-                CALL calc_pt_near_surface ( '2m' )
-             ENDIF
-          ENDIF
-!
-!--       Natural-type horizontal surfaces
-          IF ( surf_lsm_h%ns >= 1  )  THEN
-             surf => surf_lsm_h
-             CALL calc_uvw_abs
-             IF ( .NOT. neutral )  CALL calc_ol
-             CALL calc_us
-             CALL calc_scaling_parameters
-             CALL calc_surface_fluxes
-             IF ( do_output_at_2m )  THEN
-                CALL calc_pt_near_surface ( '2m' )
-             ENDIF
-          ENDIF
-!
-!--       Urban-type horizontal surfaces
-          IF ( surf_usm_h%ns >= 1  )  THEN
-             surf => surf_usm_h
-             CALL calc_uvw_abs
-             IF ( .NOT. neutral )  CALL calc_ol
-             CALL calc_us
-             CALL calc_scaling_parameters
-             CALL calc_surface_fluxes
-             IF ( do_output_at_2m )  THEN
-                CALL calc_pt_near_surface ( '2m' )
-             ENDIF
-             IF ( indoor_model )  THEN
-                CALL calc_pt_near_surface ( '10cm' )
-             ENDIF
-          ENDIF
-
+!--    surfaces are treated. 
+
+!
+!--    Default-type upward-facing horizontal surfaces
+       IF ( surf_def_h(0)%ns >= 1  )  THEN
+          surf => surf_def_h(0)
+          CALL calc_uvw_abs
+          IF ( .NOT. neutral )  CALL calc_ol
+          CALL calc_us
+          CALL calc_scaling_parameters
+          CALL calc_surface_fluxes
+          IF ( do_output_at_2m )  THEN
+             CALL calc_pt_near_surface ( '2m' )
+          ENDIF
        ENDIF
+!
+!--    Natural-type horizontal surfaces
+       IF ( surf_lsm_h%ns >= 1  )  THEN
+          surf => surf_lsm_h
+          CALL calc_uvw_abs
+          IF ( .NOT. neutral )  CALL calc_ol
+          CALL calc_us
+          CALL calc_scaling_parameters
+          CALL calc_surface_fluxes
+          IF ( do_output_at_2m )  THEN
+             CALL calc_pt_near_surface ( '2m' )
+          ENDIF
+       ENDIF
+!
+!--    Urban-type horizontal surfaces
+       IF ( surf_usm_h%ns >= 1  )  THEN
+          surf => surf_usm_h
+          CALL calc_uvw_abs
+          IF ( .NOT. neutral )  CALL calc_ol
+          CALL calc_us
+          CALL calc_scaling_parameters
+          CALL calc_surface_fluxes
+          IF ( do_output_at_2m )  THEN
+             CALL calc_pt_near_surface ( '2m' )
+          ENDIF
+          IF ( indoor_model )  THEN
+             CALL calc_pt_near_surface ( '10cm' )
+          ENDIF
+       ENDIF
+
 !
 !--    Treat downward-facing horizontal surfaces. Note, so far, these are 
@@ -549,47 +492,24 @@
 !--    flux in land-surface model in case of aero_resist_kray is not true.
        IF ( .NOT. aero_resist_kray )  THEN
-          IF ( most_method == 'circular' )  THEN
-             DO  l = 0, 1
-                IF ( surf_lsm_v(l)%ns >= 1  )  THEN
-                   surf => surf_lsm_v(l)
-!
-!--                Compute scaling parameters
-                   CALL calc_scaling_parameters
-!
-!--                Compute surface-parallel velocity
-                   CALL calc_uvw_abs_v_ugrid
-!
-!--                Compute Obukhov length
-                   IF ( .NOT. neutral )  CALL calc_ol
-!
-!--                Compute respective friction velocity on staggered grid
-                   CALL calc_us
-!
-!--                Compute respective surface fluxes for momentum and TKE
-                   CALL calc_surface_fluxes
-                ENDIF
-             ENDDO
-          ELSE
-             DO  l = 0, 1
-                IF ( surf_lsm_v(l)%ns >= 1  )  THEN
-                   surf => surf_lsm_v(l)
-!
-!--                Compute surface-parallel velocity
-                   CALL calc_uvw_abs_v_ugrid
-!
-!--                Compute Obukhov length
-                   IF ( .NOT. neutral )  CALL calc_ol
-!
-!--                Compute respective friction velocity on staggered grid
-                   CALL calc_us
-!
-!--                Compute scaling parameters
-                   CALL calc_scaling_parameters
-!
-!--                Compute respective surface fluxes for momentum and TKE
-                   CALL calc_surface_fluxes
-                ENDIF
-             ENDDO
-          ENDIF
+          DO  l = 0, 1
+             IF ( surf_lsm_v(l)%ns >= 1  )  THEN
+                surf => surf_lsm_v(l)
+!
+!--             Compute surface-parallel velocity
+                CALL calc_uvw_abs_v_ugrid
+!
+!--             Compute Obukhov length
+                IF ( .NOT. neutral )  CALL calc_ol
+!
+!--             Compute respective friction velocity on staggered grid
+                CALL calc_us
+!
+!--             Compute scaling parameters
+                CALL calc_scaling_parameters
+!
+!--             Compute respective surface fluxes for momentum and TKE
+                CALL calc_surface_fluxes
+             ENDIF
+          ENDDO
 !
 !--    No ts is required, so scaling parameters and Obukhov length do not need
@@ -652,47 +572,24 @@
 !--    flux in land-surface model in case of aero_resist_kray is not true.
        IF ( .NOT. aero_resist_kray )  THEN
-          IF ( most_method == 'circular' )  THEN
-             DO  l = 2, 3
-                IF ( surf_lsm_v(l)%ns >= 1  )  THEN
-                   surf => surf_lsm_v(l)
-!
-!--                Compute scaling parameters
-                   CALL calc_scaling_parameters
-!
-!--                Compute surface-parallel velocity
-                   CALL calc_uvw_abs_v_vgrid
-!
-!--                Compute Obukhov length
-                   IF ( .NOT. neutral )  CALL calc_ol
-!
-!--                Compute respective friction velocity on staggered grid
-                   CALL calc_us
-!
-!--                Compute respective surface fluxes for momentum and TKE
-                   CALL calc_surface_fluxes
-                ENDIF
-             ENDDO
-          ELSE
-             DO  l = 2, 3
-                IF ( surf_lsm_v(l)%ns >= 1  )  THEN
-                   surf => surf_lsm_v(l)
-!
-!--                Compute surface-parallel velocity
-                   CALL calc_uvw_abs_v_vgrid
-!
-!--                Compute Obukhov length
-                   IF ( .NOT. neutral )  CALL calc_ol
-!
-!--                Compute respective friction velocity on staggered grid
-                   CALL calc_us
-!
-!--                Compute scaling parameters
-                   CALL calc_scaling_parameters
-!
-!--                Compute respective surface fluxes for momentum and TKE
-                   CALL calc_surface_fluxes
-                ENDIF
-             ENDDO
-          ENDIF
+          DO  l = 2, 3
+             IF ( surf_lsm_v(l)%ns >= 1  )  THEN
+                surf => surf_lsm_v(l)
+!
+!--             Compute surface-parallel velocity
+                CALL calc_uvw_abs_v_vgrid
+!
+!--             Compute Obukhov length
+                IF ( .NOT. neutral )  CALL calc_ol
+!
+!--             Compute respective friction velocity on staggered grid
+                CALL calc_us
+!
+!--             Compute scaling parameters
+                CALL calc_scaling_parameters
+!
+!--             Compute respective surface fluxes for momentum and TKE
+                CALL calc_surface_fluxes
+             ENDIF
+          ENDDO
        ELSE
           DO  l = 2, 3
@@ -841,149 +738,4 @@
           IF ( surf_lsm_h%ns    >= 1 )  surf_lsm_h%ol    = 1.0E10_wp
           IF ( surf_usm_h%ns    >= 1 )  surf_usm_h%ol    = 1.0E10_wp
-       ENDIF
-
-       IF ( most_method == 'lookup' )  THEN
-
-!
-!--       Check for roughness heterogeneity. In that case terminate run and
-!--       inform user. Check for both, natural and non-natural walls.
-          IF ( surf_def_h(0)%ns >= 1 )  THEN
-             IF ( MINVAL( surf_def_h(0)%z0h ) /= MAXVAL( surf_def_h(0)%z0h )  .OR. &
-                  MINVAL( surf_def_h(0)%z0  ) /= MAXVAL( surf_def_h(0)%z0  ) )  THEN
-                terminate_run_l = .TRUE.
-             ENDIF
-          ENDIF
-          IF ( surf_lsm_h%ns >= 1 )  THEN
-             IF ( MINVAL( surf_lsm_h%z0h ) /= MAXVAL( surf_lsm_h%z0h )  .OR.       &
-                  MINVAL( surf_lsm_h%z0  ) /= MAXVAL( surf_lsm_h%z0  ) )  THEN
-                terminate_run_l = .TRUE.
-             ENDIF
-          ENDIF
-          IF ( surf_usm_h%ns >= 1 )  THEN
-             IF ( MINVAL( surf_usm_h%z0h ) /= MAXVAL( surf_usm_h%z0h )  .OR.       &
-                  MINVAL( surf_usm_h%z0  ) /= MAXVAL( surf_usm_h%z0  ) )  THEN
-                terminate_run_l = .TRUE.
-             ENDIF
-          ENDIF
-!
-!--       Check roughness homogeneity between differt surface types. 
-          IF ( surf_lsm_h%ns >= 1  .AND.  surf_def_h(0)%ns >= 1 )  THEN
-             IF ( MINVAL( surf_lsm_h%z0h ) /= MAXVAL( surf_def_h(0)%z0h )  .OR.    &
-                  MINVAL( surf_lsm_h%z0  ) /= MAXVAL( surf_def_h(0)%z0  ) )  THEN
-                terminate_run_l = .TRUE.
-             ENDIF
-          ENDIF
-          IF ( surf_usm_h%ns >= 1  .AND.  surf_def_h(0)%ns >= 1 )  THEN
-             IF ( MINVAL( surf_usm_h%z0h ) /= MAXVAL( surf_def_h(0)%z0h )  .OR.    &
-                  MINVAL( surf_usm_h%z0  ) /= MAXVAL( surf_def_h(0)%z0  ) )  THEN
-                terminate_run_l = .TRUE.
-             ENDIF
-          ENDIF
-          IF ( surf_usm_h%ns >= 1  .AND.  surf_lsm_h%ns >= 1 )  THEN
-             IF ( MINVAL( surf_usm_h%z0h ) /= MAXVAL( surf_lsm_h%z0h )  .OR.       &
-                  MINVAL( surf_usm_h%z0  ) /= MAXVAL( surf_lsm_h%z0  ) )  THEN
-                terminate_run_l = .TRUE.
-             ENDIF
-          ENDIF
-
-
-#if defined( __parallel )
-!
-!--       Make a logical OR for all processes. Force termiation of model if result
-!--       is TRUE
-          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
-          CALL MPI_ALLREDUCE( terminate_run_l, terminate_run, 1, MPI_LOGICAL,  &
-                              MPI_LOR, comm2d, ierr )
-#else
-          terminate_run = terminate_run_l
-#endif
-
-          IF ( terminate_run )  THEN
-             message_string = 'most_method = "lookup" cannot be used in ' //   &
-                              'combination with a prescribed roughness '  //   &
-                              'heterogeneity'
-             CALL message( 'surface_layer_fluxes', 'PA0116', 1, 2, 0, 6, 0 )
-          ENDIF
-
-          ALLOCATE(  zeta_tmp(0:num_steps-1) )
-
-!
-!--       Use the lowest possible value for z_mo
-          k    = nzb
-          z_mo = zu(k+1) - zw(k)
-
-!
-!--       Calculate z/L range from zeta_stretch to zeta_max using 90% of the
-!--       available steps (num_steps). The calculation is done with negative
-!--       values of zeta in order to simplify the stretching in the free
-!--       convection limit for the remaining 10% of steps.
-          zeta_tmp(0) = - zeta_max
-          num_steps_n = ( num_steps * 9 / 10 ) - 1
-          zeta_step   = (zeta_max - zeta_stretch) / REAL(num_steps_n)
-
-          DO li = 1, num_steps_n
-             zeta_tmp(li) = zeta_tmp(li-1) + zeta_step
-          ENDDO
-
-!
-!--       Calculate stretching factor for the free convection range
-          DO  WHILE ( ABS( (regr-regr_old) / regr_old ) > 1.0E-10_wp )
-             regr_old = regr
-             regr = ( 1.0_wp - ( -zeta_min / zeta_step ) * ( 1.0_wp - regr )   &
-                    )**( 10.0_wp / REAL(num_steps) )
-          ENDDO
-
-!
-!--       Calculate z/L range from zeta_min to zeta_stretch
-          DO li = num_steps_n+1, num_steps-1
-             zeta_tmp(li) = zeta_tmp(li-1) + zeta_step
-             zeta_step = zeta_step * regr
-          ENDDO
-
-!
-!--       Save roughness lengths to temporary variables
-          IF ( surf_def_h(0)%ns >= 1  )  THEN
-             z0h_min = surf_def_h(0)%z0h(1)
-             z0_min  = surf_def_h(0)%z0(1)
-          ELSEIF ( surf_lsm_h%ns >= 1 )  THEN
-             z0h_min = surf_lsm_h%z0h(1)
-             z0_min  = surf_lsm_h%z0(1)
-          ELSEIF ( surf_usm_h%ns >= 1 )  THEN
-             z0h_min = surf_usm_h%z0h(1)
-             z0_min  = surf_usm_h%z0(1)
-          ENDIF         
-!
-!--       Calculate lookup table for the Richardson number versus Obukhov length
-!--       The Richardson number (rib) is defined depending on the choice of 
-!--       boundary conditions for temperature
-          IF ( ibc_pt_b == 1 )  THEN
-             DO li = 0, num_steps-1
-                ol_tab(li)  = - z_mo / zeta_tmp(num_steps-1-li)
-                rib_tab(li) = z_mo / ol_tab(li)  / ( LOG( z_mo / z0_min )       &
-                                                - psi_m( z_mo / ol_tab(li) )    &
-                                                + psi_m( z0_min / ol_tab(li) )  &
-                                                  )**3
-             ENDDO  
-          ELSE
-             DO li = 0, num_steps-1
-                ol_tab(li)  = - z_mo / zeta_tmp(num_steps-1-li)
-                rib_tab(li) = z_mo / ol_tab(li)  * ( LOG( z_mo / z0h_min )     &
-                                              - psi_h( z_mo / ol_tab(li) )     &
-                                              + psi_h( z0h_min / ol_tab(li) )  &
-                                            )                                  &
-                                          / ( LOG( z_mo / z0_min )             &
-                                              - psi_m( z_mo / ol_tab(li) )     &
-                                              + psi_m( z0_min / ol_tab(li) )   &
-                                            )**2
-             ENDDO
-          ENDIF
-
-!
-!--       Determine minimum values of rib in the lookup table. Set upper limit 
-!--       to critical Richardson number (0.25)
-          rib_min  = MINVAL(rib_tab)
-          rib_max  = 0.25 !MAXVAL(rib_tab)
-
-          DEALLOCATE( zeta_tmp )
        ENDIF
 
@@ -1256,42 +1008,40 @@
                        ol_u      !< Upper bound of L for Newton iteration
 
-       IF ( TRIM( most_method ) /= 'circular' )  THEN
-!
-!--       Evaluate bulk Richardson number (calculation depends on
-!--       definition based on setting of boundary conditions
-          IF ( ibc_pt_b /= 1 )  THEN
-             IF ( humidity )  THEN
-                !$OMP PARALLEL DO PRIVATE( z_mo )
-                DO  m = 1, surf%ns
-
-                   z_mo = surf%z_mo(m)
-
-                   surf%rib(m) = g * z_mo                                      &
-                                 * ( surf%vpt1(m) - surf%vpt_surface(m) )      &
-                                 / ( surf%uvw_abs(m)**2 * surf%vpt1(m)         &
-                                 + 1.0E-20_wp )
-                ENDDO
-             ELSE
-                !$OMP PARALLEL DO PRIVATE( z_mo )
-                DO  m = 1, surf%ns
-
-                   z_mo = surf%z_mo(m)
-
-                   surf%rib(m) = g * z_mo                                      &
-                                 * ( surf%pt1(m) - surf%pt_surface(m) )        &
-                                 / ( surf%uvw_abs(m)**2 * surf%pt1(m)          &
-                                 + 1.0E-20_wp )
-                ENDDO
-             ENDIF
+!
+!--    Evaluate bulk Richardson number (calculation depends on
+!--    definition based on setting of boundary conditions
+       IF ( ibc_pt_b /= 1 )  THEN
+          IF ( humidity )  THEN
+             !$OMP PARALLEL DO PRIVATE( z_mo )
+             DO  m = 1, surf%ns
+
+                z_mo = surf%z_mo(m)
+
+                surf%rib(m) = g * z_mo                                         &
+                              * ( surf%vpt1(m) - surf%vpt_surface(m) )         &
+                              / ( surf%uvw_abs(m)**2 * surf%vpt1(m)            &
+                              + 1.0E-20_wp )
+             ENDDO
           ELSE
-             IF ( humidity )  THEN
-                !$OMP PARALLEL DO PRIVATE( k, z_mo )
-                DO  m = 1, surf%ns
-
-                   k   = surf%k(m)
-
-                   z_mo = surf%z_mo(m)
-
-                   surf%rib(m) = - g * z_mo * ( ( 1.0_wp + 0.61_wp             &
+             !$OMP PARALLEL DO PRIVATE( z_mo )
+             DO  m = 1, surf%ns
+
+                z_mo = surf%z_mo(m)
+
+                surf%rib(m) = g * z_mo                                         &
+                              * ( surf%pt1(m) - surf%pt_surface(m) )           &
+                              / ( surf%uvw_abs(m)**2 * surf%pt1(m) + 1.0E-20_wp )
+             ENDDO
+          ENDIF
+       ELSE
+          IF ( humidity )  THEN
+             !$OMP PARALLEL DO PRIVATE( k, z_mo )
+             DO  m = 1, surf%ns
+
+                k   = surf%k(m)
+
+                z_mo = surf%z_mo(m)
+
+                surf%rib(m) = - g * z_mo * ( ( 1.0_wp + 0.61_wp                &
                            * surf%qv1(m) ) * surf%shf(m) + 0.61_wp             &
                            * surf%pt1(m) * surf%qsws(m) ) *                    &
@@ -1299,272 +1049,144 @@
                          ( surf%uvw_abs(m)**3 * surf%vpt1(m) * kappa**2        &
                            + 1.0E-20_wp )
-                ENDDO
+             ENDDO
+          ELSE
+             !$OMP PARALLEL DO PRIVATE( k, z_mo )
+             !$ACC PARALLEL LOOP PRIVATE(k, z_mo) &
+             !$ACC PRESENT(surf, drho_air_zw)
+             DO  m = 1, surf%ns
+
+                k   = surf%k(m)
+
+                z_mo = surf%z_mo(m)
+
+                surf%rib(m) = - g * z_mo * surf%shf(m) *                    &
+                                     drho_air_zw(k-1)            /          &
+                        ( surf%uvw_abs(m)**3 * surf%pt1(m) * kappa**2       &
+                        + 1.0E-20_wp )
+             ENDDO
+          ENDIF
+       ENDIF
+
+!
+!--    Calculate the Obukhov length using Newton iteration
+       !$ACC PARALLEL LOOP PRIVATE(i, j, z_mo) &
+       !$ACC PRIVATE(ol_old, ol_m, ol_l, ol_u, f, f_d_ol) &
+       !$ACC PRESENT(surf)
+       DO  m = 1, surf%ns
+
+          i   = surf%i(m)            
+          j   = surf%j(m)
+
+          z_mo = surf%z_mo(m)
+
+!
+!--       Store current value in case the Newton iteration fails
+          ol_old = surf%ol(m)
+
+!
+!--       Ensure that the bulk Richardson number and the Obukhov 
+!--       length have the same sign
+          IF ( surf%rib(m) * surf%ol(m) < 0.0_wp  .OR.                      &
+               ABS( surf%ol(m) ) == ol_max )  THEN
+             IF ( surf%rib(m) > 1.0_wp ) surf%ol(m) =  0.01_wp
+             IF ( surf%rib(m) < 0.0_wp ) surf%ol(m) = -0.01_wp
+          ENDIF
+!
+!--       Iteration to find Obukhov length
+          iter = 0
+          DO
+             iter = iter + 1
+!
+!--          In case of divergence, use the value of the previous time step
+             IF ( iter > 1000 )  THEN
+                surf%ol(m) = ol_old
+                EXIT
+             ENDIF
+
+             ol_m = surf%ol(m)
+             ol_l = ol_m - 0.001_wp * ol_m
+             ol_u = ol_m + 0.001_wp * ol_m
+
+
+             IF ( ibc_pt_b /= 1 )  THEN
+!
+!--             Calculate f = Ri - [...]/[...]^2 = 0
+                f = surf%rib(m) - ( z_mo / ol_m ) * (                          &
+                                              LOG( z_mo / surf%z0h(m) )        &
+                                              - psi_h( z_mo / ol_m )           &
+                                              + psi_h( surf%z0h(m) /           &
+                                                       ol_m )                  &  
+                                                     )                         &
+                                           / ( LOG( z_mo / surf%z0(m) )        &
+                                              - psi_m( z_mo / ol_m )           &
+                                              + psi_m( surf%z0(m) /  ol_m )    &
+                                                 )**2
+
+!
+!--              Calculate df/dL
+                 f_d_ol = ( - ( z_mo / ol_u ) * ( LOG( z_mo /                  &
+                                                          surf%z0h(m) )        &
+                                         - psi_h( z_mo / ol_u )                &
+                                         + psi_h( surf%z0h(m) / ol_u )         &
+                                           )                                   &
+                                         / ( LOG( z_mo / surf%z0(m) )          &
+                                         - psi_m( z_mo / ol_u )                &
+                                         + psi_m( surf%z0(m) / ol_u )          &
+                                           )**2                                &
+                        + ( z_mo / ol_l ) * ( LOG( z_mo / surf%z0h(m) )        &
+                                         - psi_h( z_mo / ol_l )                &
+                                         + psi_h( surf%z0h(m) / ol_l )         &
+                                           )                                   &
+                                         / ( LOG( z_mo / surf%z0(m) )          &
+                                         - psi_m( z_mo / ol_l )                &
+                                         + psi_m( surf%z0(m) / ol_l )          &
+                                           )**2                                &
+                          ) / ( ol_u - ol_l )
              ELSE
-                !$OMP PARALLEL DO PRIVATE( k, z_mo )
-                !$ACC PARALLEL LOOP PRIVATE(k, z_mo) &
-                !$ACC PRESENT(surf, drho_air_zw)
-                DO  m = 1, surf%ns
-
-                   k   = surf%k(m)
-
-                   z_mo = surf%z_mo(m)
-
-                   surf%rib(m) = - g * z_mo * surf%shf(m) *                    &
-                                        drho_air_zw(k-1)            /          &
-                        ( surf%uvw_abs(m)**3 * surf%pt1(m) * kappa**2          &
-                           + 1.0E-20_wp )
-                ENDDO
+!
+!--             Calculate f = Ri - 1 /[...]^3 = 0
+                f = surf%rib(m) - ( z_mo / ol_m ) /                            &
+                                             ( LOG( z_mo / surf%z0(m) )        &
+                                         - psi_m( z_mo / ol_m )                &
+                                         + psi_m( surf%z0(m) / ol_m )          &
+                                             )**3
+
+!
+!--             Calculate df/dL
+                f_d_ol = ( - ( z_mo / ol_u ) / ( LOG( z_mo / surf%z0(m) )      &
+                                         - psi_m( z_mo / ol_u )                &
+                                         + psi_m( surf%z0(m) / ol_u )          &
+                                                  )**3                         &
+                           + ( z_mo / ol_l ) / ( LOG( z_mo / surf%z0(m) )      &
+                                         - psi_m( z_mo / ol_l )                &
+                                         + psi_m( surf%z0(m) / ol_l )          &
+                                            )**3                               &
+                          ) / ( ol_u - ol_l )
              ENDIF
-          ENDIF
-
-       ENDIF
-
-
-!
-!--    Calculate the Obukhov length either using a Newton iteration
-!--    method, via a lookup table, or using the old circular way
-       IF ( TRIM( most_method ) == 'newton' )  THEN
-
-          !$ACC PARALLEL LOOP PRIVATE(i, j, z_mo) &
-          !$ACC PRIVATE(ol_old, ol_m, ol_l, ol_u, f, f_d_ol) &
-          !$ACC PRESENT(surf)
-          DO  m = 1, surf%ns
-
-             i   = surf%i(m)            
-             j   = surf%j(m)
-
-             z_mo = surf%z_mo(m)
-
-!
-!--          Store current value in case the Newton iteration fails
-             ol_old = surf%ol(m)
+!
+!--          Calculate new L
+             surf%ol(m) = ol_m - f / f_d_ol
 
 !
 !--          Ensure that the bulk Richardson number and the Obukhov 
-!--          length have the same sign
-             IF ( surf%rib(m) * surf%ol(m) < 0.0_wp  .OR.                      &
-                  ABS( surf%ol(m) ) == ol_max )  THEN
-                IF ( surf%rib(m) > 1.0_wp ) surf%ol(m) =  0.01_wp
-                IF ( surf%rib(m) < 0.0_wp ) surf%ol(m) = -0.01_wp
+!--          length have the same sign and ensure convergence.
+             IF ( surf%ol(m) * ol_m < 0.0_wp )  surf%ol(m) = ol_m * 0.5_wp
+!
+!--          If unrealistic value occurs, set L to the maximum
+!--          value that is allowed
+             IF ( ABS( surf%ol(m) ) > ol_max )  THEN
+                surf%ol(m) = ol_max
+                EXIT
              ENDIF
 !
-!--          Iteration to find Obukhov length
-             iter = 0
-             DO
-                iter = iter + 1
-!
-!--             In case of divergence, use the value of the previous time step
-                IF ( iter > 1000 )  THEN
-                   surf%ol(m) = ol_old
-                   EXIT
-                ENDIF
-
-                ol_m = surf%ol(m)
-                ol_l = ol_m - 0.001_wp * ol_m
-                ol_u = ol_m + 0.001_wp * ol_m
-
-
-                IF ( ibc_pt_b /= 1 )  THEN
-!
-!--                Calculate f = Ri - [...]/[...]^2 = 0
-                   f = surf%rib(m) - ( z_mo / ol_m ) * (                       &
-                                                 LOG( z_mo / surf%z0h(m) )     &
-                                                 - psi_h( z_mo / ol_m )        &
-                                                 + psi_h( surf%z0h(m) /        &
-                                                          ol_m )               &
-                                                               )               &
-                                              / ( LOG( z_mo / surf%z0(m) )     &
-                                                 - psi_m( z_mo / ol_m )        &
-                                                 + psi_m( surf%z0(m) /         &
-                                                          ol_m )               &
-                                                 )**2
-
-!
-!--                 Calculate df/dL
-                    f_d_ol = ( - ( z_mo / ol_u ) * ( LOG( z_mo /               &
-                                                             surf%z0h(m) )     &
-                                            - psi_h( z_mo / ol_u )             &
-                                            + psi_h( surf%z0h(m) / ol_u )      &
-                                              )                                &
-                                            / ( LOG( z_mo / surf%z0(m) )       &
-                                            - psi_m( z_mo / ol_u )             &
-                                            + psi_m( surf%z0(m) / ol_u )       &
-                                              )**2                             &
-                           + ( z_mo / ol_l ) * ( LOG( z_mo / surf%z0h(m) )     &
-                                            - psi_h( z_mo / ol_l )             &
-                                            + psi_h( surf%z0h(m) / ol_l )      &
-                                              )                                &
-                                            / ( LOG( z_mo / surf%z0(m) )       &
-                                            - psi_m( z_mo / ol_l )             &
-                                            + psi_m( surf%z0(m) / ol_l )       &
-                                              )**2                             &
-                             ) / ( ol_u - ol_l )
-                ELSE
-!
-!--                Calculate f = Ri - 1 /[...]^3 = 0
-                   f = surf%rib(m) - ( z_mo / ol_m ) /                         &
-                                                ( LOG( z_mo / surf%z0(m) )     &
-                                            - psi_m( z_mo / ol_m )             &
-                                            + psi_m( surf%z0(m) / ol_m )       &
-                                                )**3
-
-!
-!--                Calculate df/dL
-                   f_d_ol = ( - ( z_mo / ol_u ) / ( LOG( z_mo /                &
-                                                            surf%z0(m) )       &
-                                            - psi_m( z_mo / ol_u )             &
-                                            + psi_m( surf%z0(m) / ol_u )       &
-                                                     )**3                      &
-                           + ( z_mo / ol_l ) / ( LOG( z_mo / surf%z0(m) )      &
-                                            - psi_m( z_mo / ol_l )             &
-                                            + psi_m( surf%z0(m) / ol_l )       &
-                                               )**3                            &
-                                  ) / ( ol_u - ol_l )
-                ENDIF
-!
-!--             Calculate new L
-                surf%ol(m) = ol_m - f / f_d_ol
-
-!
-!--             Ensure that the bulk Richardson number and the Obukhov 
-!--             length have the same sign and ensure convergence.
-                IF ( surf%ol(m) * ol_m < 0.0_wp )  surf%ol(m) = ol_m * 0.5_wp
-!
-!--             If unrealistic value occurs, set L to the maximum
-!--             value that is allowed
-                IF ( ABS( surf%ol(m) ) > ol_max )  THEN
-                   surf%ol(m) = ol_max
-                   EXIT
-                ENDIF
-!
-!--             Check for convergence
-                IF ( ABS( ( surf%ol(m) - ol_m ) /                              &
-                     surf%ol(m) ) < 1.0E-4_wp )  THEN
-                   EXIT
-                ELSE
-                   CYCLE
-                ENDIF
-
-             ENDDO
+!--          Check for convergence
+             IF ( ABS( ( surf%ol(m) - ol_m ) /  surf%ol(m) ) < 1.0E-4_wp )  THEN
+                EXIT
+             ELSE
+                CYCLE
+             ENDIF
+
           ENDDO
-
-       ELSEIF ( TRIM( most_method ) == 'lookup' )  THEN
-
-          !$OMP PARALLEL DO PRIVATE( i, j, z_mo, li ) FIRSTPRIVATE( li_bnd ) LASTPRIVATE( li_bnd )
-          DO  m = 1, surf%ns
-
-             i   = surf%i(m)            
-             j   = surf%j(m)
-!
-!--          If the bulk Richardson number is outside the range of the lookup
-!--          table, set it to the exceeding threshold value
-             IF ( surf%rib(m) < rib_min )  surf%rib(m) = rib_min
-             IF ( surf%rib(m) > rib_max )  surf%rib(m) = rib_max
-
-!
-!--          Find the correct index bounds for linear interpolation. As the
-!--          Richardson number will not differ very much from time step to
-!--          time step , use the index from the last step and search in the 
-!--          correct direction
-             li = li_bnd
-             IF ( rib_tab(li) - surf%rib(m) > 0.0_wp )  THEN
-                DO  WHILE ( rib_tab(li-1) - surf%rib(m) > 0.0_wp  .AND.  li > 0 )
-                   li = li-1
-                ENDDO
-             ELSE
-                DO  WHILE ( rib_tab(li) - surf%rib(m) < 0.0_wp                 &
-                           .AND.  li < num_steps-1 )
-                   li = li+1
-                ENDDO
-             ENDIF
-             li_bnd = li
-
-!
-!--          Linear interpolation to find the correct value of z/L
-             surf%ol(m) = ( ol_tab(li-1) + ( ol_tab(li) - ol_tab(li-1) )       &
-                         / (  rib_tab(li) - rib_tab(li-1) )                    &
-                         * ( surf%rib(m) - rib_tab(li-1) ) )
-
-          ENDDO
-
-       ELSEIF ( TRIM( most_method ) == 'circular' )  THEN
-
-          IF ( .NOT. humidity )  THEN
-             !$OMP PARALLEL DO PRIVATE( z_mo )
-             DO  m = 1, surf%ns
-
-                z_mo = surf%z_mo(m)
-
-                surf%ol(m) =  ( surf%pt1(m) *  surf%us(m)**2 ) /               &
-                                       ( kappa * g *                           &
-                                         surf%ts(m) + 1E-30_wp )
-!
-!--             Limit the value range of the Obukhov length.
-!--             This is necessary for very small velocities (u,v --> 1), because
-!--             the absolute value of ol can then become very small, which in
-!--             consequence would result in very large shear stresses and very
-!--             small momentum fluxes (both are generally unrealistic).
-                IF ( ( z_mo / ( surf%ol(m) + 1E-30_wp ) ) < zeta_min )         &
-                   surf%ol(m) = z_mo / zeta_min
-                IF ( ( z_mo / ( surf%ol(m) + 1E-30_wp ) ) > zeta_max )         &
-                   surf%ol(m) = z_mo / zeta_max
-
-             ENDDO
-          ELSEIF ( bulk_cloud_model  .OR.  cloud_droplets )  THEN
-             !$OMP PARALLEL DO PRIVATE( i, j, k, z_mo )
-             DO  m = 1, surf%ns
-
-                i   = surf%i(m)            
-                j   = surf%j(m)
-                k   = surf%k(m)
-
-                z_mo = surf%z_mo(m)
-
-
-                surf%ol(m) =  ( surf%vpt1(m) * surf%us(m)**2 ) /               &
-                    ( kappa * g * ( surf%ts(m) +                               &
-                         0.61_wp * surf%pt1(m) * surf%us(m)                    &
-                       + 0.61_wp * surf%qv1(m) * surf%ts(m) -                  &
-                                   surf%ts(m)  * ql(k,j,i) ) + 1E-30_wp )
-!
-!--             Limit the value range of the Obukhov length.
-!--             This is necessary for very small velocities (u,v --> 1), because
-!--             the absolute value of ol can then become very small, which in
-!--             consequence would result in very large shear stresses and very
-!--             small momentum fluxes (both are generally unrealistic).
-                IF ( ( z_mo / ( surf%ol(m) + 1E-30_wp ) ) < zeta_min )         &
-                   surf%ol(m) = z_mo / zeta_min
-                IF ( ( z_mo / ( surf%ol(m) + 1E-30_wp ) ) > zeta_max )         &
-                   surf%ol(m) = z_mo / zeta_max
-
-             ENDDO
-          ELSE
-
-             !$OMP PARALLEL DO PRIVATE( z_mo )
-             DO  m = 1, surf%ns
-
-                z_mo = surf%z_mo(m)
-
-                surf%ol(m) =  ( surf%vpt1(m) *  surf%us(m)**2 ) /              &
-                    ( kappa * g * ( surf%ts(m) + 0.61_wp * surf%pt1(m) *       &
-                                    surf%qs(m) + 0.61_wp * surf%qv1(m) *       &
-                                    surf%ts(m) ) + 1E-30_wp )
-
-!
-!--             Limit the value range of the Obukhov length.
-!--             This is necessary for very small velocities (u,v --> 1), because
-!--             the absolute value of ol can then become very small, which in
-!--             consequence would result in very large shear stresses and very
-!--             small momentum fluxes (both are generally unrealistic).
-                IF ( ( z_mo / ( surf%ol(m) + 1E-30_wp ) ) < zeta_min )         &
-                   surf%ol(m) = z_mo / zeta_min
-                IF ( ( z_mo / ( surf%ol(m) + 1E-30_wp ) ) > zeta_max )         &
-                   surf%ol(m) = z_mo / zeta_max
-
-             ENDDO
-
-          ENDIF
-
-       ENDIF
+       ENDDO
 
     END SUBROUTINE calc_ol
Index: palm/trunk/SOURCE/write_restart_data_mod.f90
===================================================================
--- palm/trunk/SOURCE/write_restart_data_mod.f90	(revision 3667)
+++ palm/trunk/SOURCE/write_restart_data_mod.f90	(revision 3668)
@@ -163,5 +163,5 @@
 
 
-       binary_version_global = '4.7'
+       binary_version_global = '4.8'
 
        CALL wrd_write_string( 'binary_version_global' )
@@ -466,7 +466,4 @@
        WRITE ( 14 )  momentum_advec
 
-       CALL wrd_write_string( 'most_method' ) 
-       WRITE ( 14 )  most_method
-
        CALL wrd_write_string( 'netcdf_precision' ) 
        WRITE ( 14 )  netcdf_precision
