Index: palm/trunk/SOURCE/surface_layer_fluxes_mod.f90
===================================================================
--- palm/trunk/SOURCE/surface_layer_fluxes_mod.f90	(revision 4360)
+++ palm/trunk/SOURCE/surface_layer_fluxes_mod.f90	(revision 4366)
@@ -26,4 +26,7 @@
 ! -----------------
 ! $Id$
+! vector version for calculation of Obukhov length via Newton iteration added
+! 
+! 4360 2020-01-07 11:25:50Z suehring
 ! Calculation of diagnostic-only 2-m potential temperature moved to 
 ! diagnostic_output_quantities.
@@ -112,5 +115,5 @@
                constant_waterflux, coupling_mode,                              &
                debug_output_timestep,                                          &
-               humidity,                                                       &
+               humidity, loop_optimization,                                    &
                ibc_e_b, ibc_pt_b, indoor_model,                                &
                land_surface, large_scale_forcing, lsf_surf, message_string,    &
@@ -859,4 +862,6 @@
        INTEGER(iwp) ::  m       !< loop variable over all horizontal wall elements 
 
+       LOGICAL, DIMENSION(surf%ns) ::  convergence_reached  !< convergence switch for vectorization
+
        REAL(wp)     :: f,      & !< Function for Newton iteration: f = Ri - [...]/[...]^2 = 0
                        f_d_ol, & !< Derivative of f
@@ -866,4 +871,7 @@
                        ol_u      !< Upper bound of L for Newton iteration
 
+       REAL(wp), DIMENSION(surf%ns) ::  ol_old_vec  !< temporary array required for vectorization
+       REAL(wp), DIMENSION(surf%ns) ::  z_mo_vec    !< temporary array required for vectorization
+
 !
 !--    Evaluate bulk Richardson number (calculation depends on
@@ -926,134 +934,265 @@
        ENDIF
 
-!
-!--    Calculate the Obukhov length using Newton iteration
-       !$OMP PARALLEL DO PRIVATE(i, j, z_mo) &
-       !$OMP PRIVATE(ol_old, ol_m, ol_l, ol_u, f, f_d_ol)
-       !$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
+       IF ( loop_optimization == 'cache' )  THEN
+!
+!--       Calculate the Obukhov length using Newton iteration
+          !$OMP PARALLEL DO PRIVATE(i, j, z_mo) &
+          !$OMP PRIVATE(ol_old, ol_m, ol_l, ol_u, f, f_d_ol)
+          !$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
+!
+!--                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
+!
+!--             Assure that Obukhov length does not become zero. If the limit is
+!--             reached, exit the loop.
+                IF ( ABS( surf%ol(m) ) < 1E-5_wp )  THEN
+                   surf%ol(m) = SIGN( 1E-5_wp, surf%ol(m) )
+                   EXIT
+                ENDIF
+!
+!--             Check for convergence
+                IF ( ABS( ( surf%ol(m) - ol_m ) /  surf%ol(m) ) < 1.0E-4_wp )  EXIT
+
+             ENDDO
+          ENDDO
+
+!
+!--    Vector Version
+       ELSE
+!
+!--       Calculate the Obukhov length using Newton iteration
+!--       First set arrays required for vectorization
+          DO  m = 1, surf%ns
+
+             z_mo_vec(m) = surf%z_mo(m)
+
+!
+!--          Store current value in case the Newton iteration fails
+             ol_old_vec(m) = 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
+
+          ENDDO
+
 !
 !--       Iteration to find Obukhov length
+          convergence_reached(:) = .FALSE.
           iter = 0
           DO
+
              iter = iter + 1
 !
-!--          In case of divergence, use the value of the previous time step
+!--          In case of divergence, use the value(s) of the previous time step
              IF ( iter > 1000 )  THEN
-                surf%ol(m) = ol_old
+                DO  m = 1, surf%ns
+                   IF ( .NOT. convergence_reached(m) )  surf%ol(1:surf%ns) = ol_old
+                ENDDO
                 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
-!
-!--          Assure that Obukhov length does not become zero. If the limit is
-!--          reached, exit the loop.
-             IF ( ABS( surf%ol(m) ) < 1E-5_wp )  THEN
-                surf%ol(m) = SIGN( 1E-5_wp, surf%ol(m) )
-                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
-       ENDDO
+
+             DO  m = 1, surf%ns
+
+                IF ( convergence_reached(m) )  CYCLE
+
+                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_vec(m) / ol_m ) * ( LOG( z_mo_vec(m) / surf%z0h(m) )   &
+                                                              - psi_h( z_mo_vec(m) / ol_m )        &
+                                                              + psi_h( surf%z0h(m) / ol_m )        &
+                                                              ) /                                  &
+                                                              ( LOG( z_mo_vec(m) / surf%z0(m) )    &
+                                                             - psi_m( z_mo_vec(m) / ol_m )         &
+                                                             + psi_m( surf%z0(m) /  ol_m )         &
+                                                              )**2
+
+!
+!--                Calculate df/dL
+                   f_d_ol = ( - ( z_mo_vec(m) / ol_u ) * ( LOG( z_mo_vec(m) / surf%z0h(m) )        &
+                                                         - psi_h( z_mo_vec(m) / ol_u )             &
+                                                         + psi_h( surf%z0h(m) / ol_u )             &
+                                                         ) /                                       &
+                                                         ( LOG( z_mo_vec(m) / surf%z0(m) )         &
+                                                         - psi_m( z_mo_vec(m) / ol_u )             &
+                                                         + psi_m( surf%z0(m) / ol_u )              &
+                                                         )**2                                      &
+                              + ( z_mo_vec(m) / ol_l ) * ( LOG( z_mo_vec(m) / surf%z0h(m) )        &
+                                                         - psi_h( z_mo_vec(m) / ol_l )             &
+                                                         + psi_h( surf%z0h(m) / ol_l )             &
+                                                         ) /                                       &
+                                                         ( LOG( z_mo_vec(m) / surf%z0(m) )         &
+                                                         - psi_m( z_mo_vec(m) / 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_vec(m) / ol_m ) / ( LOG( z_mo_vec(m) / surf%z0(m) )    &
+                                                              - psi_m( z_mo_vec(m) / ol_m )        &
+                                                              + psi_m( surf%z0(m) / ol_m )         &
+                                                              )**3
+
+!
+!--                Calculate df/dL
+                   f_d_ol = ( - ( z_mo_vec(m) / ol_u ) / ( LOG( z_mo_vec(m) / surf%z0(m) )         &
+                                                         - psi_m( z_mo_vec(m) / ol_u )             &
+                                                         + psi_m( surf%z0(m) / ol_u )              &
+                                                         )**3                                      &
+                              + ( z_mo_vec(m) / ol_l ) / ( LOG( z_mo_vec(m) / surf%z0(m) )         &
+                                                         - psi_m( z_mo_vec(m) / 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
+
+!
+!--             Check for convergence
+!--             This check does not modify surf%ol, therefore this is done first
+                IF ( ABS( ( surf%ol(m) - ol_m ) /  surf%ol(m) ) < 1.0E-4_wp )  THEN
+                   convergence_reached(m) = .TRUE.
+                ENDIF
+!
+!--             If unrealistic value occurs, set L to the maximum allowed value
+                IF ( ABS( surf%ol(m) ) > ol_max )  THEN
+                   surf%ol(m) = ol_max
+                   convergence_reached(m) = .TRUE.
+                ENDIF
+
+             ENDDO
+!
+!--          Assure that Obukhov length does not become zero
+             DO  m = 1, surf%ns
+                IF ( convergence_reached(m) )  CYCLE
+                IF ( ABS( surf%ol(m) ) < 1E-5_wp )  THEN
+                   surf%ol(m) = SIGN( 10E-6_wp, surf%ol(m) )
+                   convergence_reached(m) = .TRUE.
+                ENDIF
+             ENDDO
+
+             IF ( ALL( convergence_reached ) )  EXIT
+
+          ENDDO  ! end of iteration loop
+
+       ENDIF  ! end of vector branch
 
     END SUBROUTINE calc_ol
