Index: palm/trunk/SOURCE/surface_layer_fluxes_mod.f90
===================================================================
--- palm/trunk/SOURCE/surface_layer_fluxes_mod.f90	(revision 4592)
+++ palm/trunk/SOURCE/surface_layer_fluxes_mod.f90	(revision 4593)
@@ -20,9 +20,16 @@
 ! Current revisions:
 ! -----------------
-! 
-! 
+!
+!
 ! Former revisions:
 ! -----------------
 ! $Id$
+! - Pre-calculate ln(z/z0) at each timestep in order to reduce the number of log-calculations
+! - Bugfix - add missing density to fluxes of passive-scalars, chemistry and cloud-phyiscal
+!   quantities at upward-facing surfaces
+! - Move if-statement out of inner loop
+! - Remove unnecessary index referencing
+! 
+! 4562 2020-06-12 08:38:47Z raasch
 ! File re-formatted to follow the PALM coding standard
 !
@@ -268,6 +275,37 @@
     downward      = .FALSE. !< flag indicating downward-facing surface elements
 !
+!-- First, precalculate ln(z/z0) for all surfaces. This is done each timestep, in order to account
+!-- for time-dependent roughness or user-modifications.
+    DO  l = 0, 1
+       IF ( surf_def_h(l)%ns >= 1 )  THEN
+          surf => surf_def_h(l)
+          CALL calc_ln
+       ENDIF
+    ENDDO
+    IF ( surf_lsm_h%ns >= 1 )  THEN
+       surf => surf_lsm_h
+       CALL calc_ln
+    ENDIF
+    IF ( surf_usm_h%ns >= 1 )  THEN
+       surf => surf_usm_h
+       CALL calc_ln
+    ENDIF
+
+    DO  l = 0, 3
+       IF ( surf_def_v(l)%ns >= 1 )  THEN
+          surf => surf_def_v(l)
+          CALL calc_ln
+       ENDIF
+       IF ( surf_lsm_v(l)%ns >= 1 )  THEN
+          surf => surf_lsm_v(l)
+          CALL calc_ln
+       ENDIF
+       IF ( surf_usm_v(l)%ns >= 1 )  THEN
+          surf => surf_usm_v(l)
+          CALL calc_ln
+       ENDIF
+    ENDDO
+!
 !-- Derive potential temperature and specific humidity at first grid level from the fields pt and q
-!
 !-- First call for horizontal default-type surfaces (l=0 - upward facing, l=1 - downward facing)
     DO  l = 0, 1
@@ -314,11 +352,12 @@
 
 !
-!-- First, calculate the new Obukhov length, then new friction velocity, followed by the new scaling
+!-- First, calculate the new Obukhov length from precalculated values of log(z/z0) and wind speeds.
+!-- As a second step, then calculate new friction velocity, followed by the new scaling
 !-- parameters (th*, q*, etc.), and the new 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, here only upward-facing horizontal surfaces are treated.
-
-!
-!-- Default-type upward-facing horizontal surfaces
+!-- Note, calculation of log(z/z0) is redone each timestep, in order to account for time-dependent
+!-- values.
+!-- Start with default-type upward-facing horizontal surfaces
     IF ( surf_def_h(0)%ns >= 1 )  THEN
        surf => surf_def_h(0)
@@ -645,4 +684,29 @@
 ! Description:
 ! ------------
+!> Compute ln(z/z0).
+!--------------------------------------------------------------------------------------------------!
+ SUBROUTINE calc_ln
+
+    INTEGER(iwp) ::  m     !< running index surface elements
+
+!
+!-- Note, ln(z/z0h) and ln(z/z0q) is also calculated even if neural simulations are applied.
+!-- This is because the scalar coefficients are also used for other scalars such as passive scalars,
+!-- chemistry and aerosols.
+    !$OMP PARALLEL  DO PRIVATE( z_mo )
+    !$ACC PARALLEL LOOP PRIVATE(z_mo) &
+    !$ACC PRESENT(surf)
+    DO  m = 1, surf%ns
+       z_mo = surf%z_mo(m)
+       surf%ln_z_z0(m)  = LOG( z_mo / surf%z0(m) )
+       surf%ln_z_z0h(m) = LOG( z_mo / surf%z0h(m) )
+       surf%ln_z_z0q(m) = LOG( z_mo / surf%z0q(m) )
+    ENDDO
+
+ END SUBROUTINE calc_ln
+ 
+!--------------------------------------------------------------------------------------------------!
+! Description:
+! ------------
 !> Compute the absolute value of the horizontal velocity (relative to the surface) for horizontal
 !> surface elements. This is required by all methods.
@@ -663,37 +727,58 @@
     ibit = MERGE( 1, 0, .NOT. downward )
 
-    !$OMP PARALLEL DO PRIVATE(i, j, k, w_lfc)
-    !$ACC PARALLEL LOOP PRIVATE(i, j, k, w_lfc) &
-    !$ACC PRESENT(surf, u, v)
-    DO  m = 1, surf%ns
-       i   = surf%i(m)
-       j   = surf%j(m)
-       k   = surf%k(m)
-
-!
-!--    Calculate free convection velocity scale w_lfc is use_free_convection_scaling = .T.. This
-!--    will maintain a horizontal velocity even for very weak wind convective conditions. SIGN is
-!--    used to set w_lfc to zero under stable conditions.
-       IF ( use_free_convection_scaling )  THEN
+    IF ( use_free_convection_scaling )  THEN
+       !$OMP PARALLEL DO PRIVATE(i, j, k, w_lfc)
+       !$ACC PARALLEL LOOP PRIVATE(i, j, k, w_lfc) &
+       !$ACC PRESENT(surf, u, v)
+       DO  m = 1, surf%ns
+          i   = surf%i(m)
+          j   = surf%j(m)
+          k   = surf%k(m)
+
+!
+!--       Calculate free convection velocity scale w_lfc is use_free_convection_scaling = .T.. This
+!--       will maintain a horizontal velocity even for very weak wind convective conditions. SIGN is
+!--       used to set w_lfc to zero under stable conditions.
           w_lfc = ABS(g / surf%pt1(m) * surf%z_mo(m) * surf%shf(m))
           w_lfc = ( 0.5_wp * ( w_lfc + SIGN(w_lfc,surf%shf(m)) ) )**(0.33333_wp)
-       ELSE
-          w_lfc = 0.0_wp
-       ENDIF
-
-!
-!--    Compute the absolute value of the horizontal velocity. (relative to the surface in case the
-!--    lower surface is the ocean). Please note, in new surface modelling concept the index values
-!--    changed, i.e. the reference grid point is not the surface-grid point itself but the first
-!--    grid point outside of the topography. Note, in case of coupled ocean-atmosphere simulations
-!--    relative velocity with respect to the ocean surface is used, hence, (k-1,j,i) values are used
-!--    to calculate the absolute velocity. However, this does not apply for downward-facing walls.
-!--    To mask this, use ibit, which checks for upward/downward-facing surfaces.
-       surf%uvw_abs(m) = SQRT( ( 0.5_wp * ( u(k,j,i) + u(k,j,i+1) - ( u(k-1,j,i) + u(k-1,j,i+1) )  &
-                                            * ibit ) )**2                                          &
-                               + ( 0.5_wp * ( v(k,j,i) + v(k,j+1,i) - ( v(k-1,j,i) + v(k-1,j+1,i)  &
-                                                                        ) * ibit )                 &
+!
+!--       Compute the absolute value of the horizontal velocity. (relative to the surface in case the
+!--       lower surface is the ocean). Please note, in new surface modelling concept the index values
+!--       changed, i.e. the reference grid point is not the surface-grid point itself but the first
+!--       grid point outside of the topography. Note, in case of coupled ocean-atmosphere simulations
+!--       relative velocity with respect to the ocean surface is used, hence, (k-1,j,i) values are used
+!--       to calculate the absolute velocity. However, this does not apply for downward-facing walls.
+!--       To mask this, use ibit, which checks for upward/downward-facing surfaces.
+          surf%uvw_abs(m) = SQRT( ( 0.5_wp * ( u(k,j,i) + u(k,j,i+1) - ( u(k-1,j,i) + u(k-1,j,i+1)  &
+                                                                       ) * ibit )                   &
+                                  )**2                                                              &
+                                + ( 0.5_wp * ( v(k,j,i) + v(k,j+1,i) - ( v(k-1,j,i) + v(k-1,j+1,i)  &
+                                                                       ) * ibit )                   &
                                   )**2  + w_lfc**2 )
-    ENDDO
+       ENDDO
+    ELSE
+       !$OMP PARALLEL DO PRIVATE(i, j, k)
+       !$ACC PARALLEL LOOP PRIVATE(i, j, k) &
+       !$ACC PRESENT(surf, u, v)
+       DO  m = 1, surf%ns
+          i   = surf%i(m)
+          j   = surf%j(m)
+          k   = surf%k(m)
+!
+!--       Compute the absolute value of the horizontal velocity. (relative to the surface in case the
+!--       lower surface is the ocean). Please note, in new surface modelling concept the index values
+!--       changed, i.e. the reference grid point is not the surface-grid point itself but the first
+!--       grid point outside of the topography. Note, in case of coupled ocean-atmosphere simulations
+!--       relative velocity with respect to the ocean surface is used, hence, (k-1,j,i) values are used
+!--       to calculate the absolute velocity. However, this does not apply for downward-facing walls.
+!--       To mask this, use ibit, which checks for upward/downward-facing surfaces.
+          surf%uvw_abs(m) = SQRT( ( 0.5_wp * ( u(k,j,i) + u(k,j,i+1) - ( u(k-1,j,i) + u(k-1,j,i+1)  &
+                                                                       ) * ibit )                   &
+                                  )**2                                                              &
+                                + ( 0.5_wp * ( v(k,j,i) + v(k,j+1,i) - ( v(k-1,j,i) + v(k-1,j+1,i)  &
+                                                                       ) * ibit )                   &
+                                  )**2 )
+       ENDDO
+    ENDIF
 
  END SUBROUTINE calc_uvw_abs
@@ -980,35 +1065,34 @@
 !
 !--             Calculate f = Ri - [...]/[...]^2 = 0
-                f = surf%rib(m) - ( z_mo / ol_m ) * ( LOG( z_mo / surf%z0h(m) )                    &
+                f = surf%rib(m) - ( z_mo / ol_m ) * ( surf%ln_z_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
+                    ( surf%ln_z_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) )                         &
+                f_d_ol = ( - ( z_mo / ol_u ) * ( surf%ln_z_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 )
+                           ( surf%ln_z_z0(m) - psi_m( z_mo / ol_u )                                &
+                                             + psi_m( surf%z0(m) / ol_u ) )**2                     &
+                           + ( z_mo / ol_l ) * ( surf%ln_z_z0h(m) - psi_h( z_mo / ol_l )           &
+                                                                  + psi_h( surf%z0h(m) / ol_l ) ) /&
+                           ( surf%ln_z_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
+                    ( surf%ln_z_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) )                          &
+                f_d_ol = ( - ( z_mo / ol_u ) / ( surf%ln_z_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) )                          &
+                           + ( z_mo / ol_l ) / ( surf%ln_z_z0(m)                                   &
                                                  - psi_m( z_mo / ol_l )                            &
                                                  + psi_m( surf%z0(m) / ol_l )                      &
@@ -1095,9 +1179,9 @@
 !
 !--             Calculate f = Ri - [...]/[...]^2 = 0
-                f = surf%rib(m) - ( z_mo_vec(m) / ol_m ) * ( LOG( z_mo_vec(m) / surf%z0h(m) )      &
+                f = surf%rib(m) - ( z_mo_vec(m) / ol_m ) * ( surf%ln_z_z0(m)                       &
                                                            - psi_h( z_mo_vec(m) / ol_m )           &
                                                            + psi_h( surf%z0h(m) / ol_m )           &
                                                            ) /                                     &
-                                                           ( LOG( z_mo_vec(m) / surf%z0(m) )       &
+                                                           ( surf%ln_z_z0(m)                       &
                                                           - psi_m( z_mo_vec(m) / ol_m )            &
                                                           + psi_m( surf%z0(m) /  ol_m )            &
@@ -1106,17 +1190,17 @@
 !
 !--             Calculate df/dL
-                f_d_ol = ( - ( z_mo_vec(m) / ol_u ) * ( LOG( z_mo_vec(m) / surf%z0h(m) )           &
+                f_d_ol = ( - ( z_mo_vec(m) / ol_u ) * ( surf%ln_z_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) )            &
+                                                      ( surf%ln_z_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) )           &
+                           + ( z_mo_vec(m) / ol_l ) * ( surf%ln_z_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) )            &
+                                                      ( surf%ln_z_z0(m)                            &
                                                         - psi_m( z_mo_vec(m) / ol_l )              &
                                                         + psi_m( surf%z0(m) / ol_l )               &
@@ -1126,5 +1210,5 @@
 !
 !--             Calculate f = Ri - 1 /[...]^3 = 0
-                f = surf%rib(m) - ( z_mo_vec(m) / ol_m ) / ( LOG( z_mo_vec(m) / surf%z0(m) )       &
+                f = surf%rib(m) - ( z_mo_vec(m) / ol_m ) / ( surf%ln_z_z0(m)                       &
                                                              - psi_m( z_mo_vec(m) / ol_m )         &
                                                              + psi_m( surf%z0(m) / ol_m )          &
@@ -1133,9 +1217,9 @@
 !
 !--             Calculate df/dL
-                f_d_ol = ( - ( z_mo_vec(m) / ol_u ) / ( LOG( z_mo_vec(m) / surf%z0(m) )            &
+                f_d_ol = ( - ( z_mo_vec(m) / ol_u ) / ( surf%ln_z_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) )            &
+                           + ( z_mo_vec(m) / ol_l ) / ( surf%ln_z_z0(m)                            &
                                                         - psi_m( z_mo_vec(m) / ol_l )              &
                                                         + psi_m( surf%z0(m) / ol_l )               &
@@ -1210,5 +1294,5 @@
 !
 !--          Compute u* at the scalars' grid points
-             surf%us(m) = kappa * surf%uvw_abs(m) / ( LOG( z_mo / surf%z0(m) )                     &
+             surf%us(m) = kappa * surf%uvw_abs(m) / ( surf%ln_z_z0(m)                              &
                                                       - psi_m( z_mo / surf%ol(m) )                 &
                                                       + psi_m( surf%z0(m) / surf%ol(m) ) )
@@ -1217,12 +1301,11 @@
 !--    Compute u* at downward-facing surfaces. This case, do not consider any stability.
        ELSE
-          !$OMP PARALLEL  DO PRIVATE( z_mo )
-          !$ACC PARALLEL LOOP PRIVATE(z_mo) &
+          !$OMP PARALLEL  DO
+          !$ACC PARALLEL LOOP &
           !$ACC PRESENT(surf)
           DO  m = 1, surf%ns
-             z_mo = surf%z_mo(m)
 !
 !--          Compute u* at the scalars' grid points
-             surf%us(m) = kappa * surf%uvw_abs(m) / LOG( z_mo / surf%z0(m) )
+             surf%us(m) = kappa * surf%uvw_abs(m) / surf%ln_z_z0(m)
           ENDDO
        ENDIF
@@ -1231,10 +1314,9 @@
 !-- No stability is considered in this case.
     ELSE
-       !$OMP PARALLEL DO PRIVATE( z_mo )
-       !$ACC PARALLEL LOOP PRIVATE(z_mo) &
+       !$OMP PARALLEL DO
+       !$ACC PARALLEL LOOP &
        !$ACC PRESENT(surf)
        DO  m = 1, surf%ns
-          z_mo = surf%z_mo(m)
-          surf%us(m) = kappa * surf%uvw_abs(m) / LOG( z_mo / surf%z0(m) )
+          surf%us(m) = kappa * surf%uvw_abs(m) / surf%ln_z_z0(m)
        ENDDO
     ENDIF
@@ -1411,6 +1493,6 @@
           z_mo = surf%z_mo(m)
           surf%ts(m) = kappa * ( surf%pt1(m) - surf%pt_surface(m) )                                &
-                       / ( LOG( z_mo / surf%z0h(m) ) - psi_h( z_mo / surf%ol(m) )                  &
-                                                     + psi_h( surf%z0h(m) / surf%ol(m) ) )
+                       / ( surf%ln_z_z0h(m) - psi_h( z_mo / surf%ol(m) )                           &
+                                            + psi_h( surf%z0h(m) / surf%ol(m) ) )
        ENDDO
 
@@ -1482,6 +1564,6 @@
                 z_mo = surf%z_mo(m)
                 surf%qs(m) = kappa * ( surf%qv1(m) - surf%q_surface(m) )                           &
-                             / ( LOG( z_mo / surf%z0q(m) ) - psi_h( z_mo / surf%ol(m) )            &
-                                                           + psi_h( surf%z0q(m) / surf%ol(m) ) )
+                             / ( surf%ln_z_z0q(m) - psi_h( z_mo / surf%ol(m) )                     &
+                                                  + psi_h( surf%z0q(m) / surf%ol(m) ) )
              ENDDO
           ELSE
@@ -1493,6 +1575,6 @@
                 z_mo = surf%z_mo(m)
                 surf%qs(m) = kappa * ( q(k,j,i) - q(k-1,j,i) )                                     &
-                             / ( LOG( z_mo / surf%z0q(m) ) - psi_h( z_mo / surf%ol(m) )            &
-                                                           + psi_h( surf%z0q(m) / surf%ol(m) ) )
+                             / ( surf%ln_z_z0q(m) - psi_h( z_mo / surf%ol(m) )                     &
+                                                  + psi_h( surf%z0q(m) / surf%ol(m) ) )
              ENDDO
           ENDIF
@@ -1536,6 +1618,6 @@
 
              surf%ss(m) = kappa * ( s(k,j,i) - s(k-1,j,i) )                                        &
-                          / ( LOG( z_mo / surf%z0h(m) ) - psi_h( z_mo / surf%ol(m) )               &
-                                                        + psi_h( surf%z0h(m) / surf%ol(m) ) )
+                          / ( surf%ln_z_z0h(m) - psi_h( z_mo / surf%ol(m) )                        &
+                                               + psi_h( surf%z0h(m) / surf%ol(m) ) )
           ENDDO
        ENDIF
@@ -1594,10 +1676,10 @@
 
           surf%qcs(m) = kappa * ( qc(k,j,i) - qc(k-1,j,i) )                                        &
-                        / ( LOG( z_mo / surf%z0q(m) ) - psi_h( z_mo / surf%ol(m) )                 &
-                                                      + psi_h( surf%z0q(m) / surf%ol(m) ) )
+                        / ( surf%ln_z_z0q(m) - psi_h( z_mo / surf%ol(m) )                          &
+                                             + psi_h( surf%z0q(m) / surf%ol(m) ) )
 
           surf%ncs(m) = kappa * ( nc(k,j,i) - nc(k-1,j,i) )                                        &
-                        / ( LOG( z_mo / surf%z0q(m) ) - psi_h( z_mo / surf%ol(m) )                 &
-                                                      + psi_h( surf%z0q(m) / surf%ol(m) ) )
+                        / ( surf%ln_z_z0q(m) - psi_h( z_mo / surf%ol(m) )                          &
+                                             + psi_h( surf%z0q(m) / surf%ol(m) ) )
        ENDDO
 
@@ -1616,10 +1698,10 @@
 
           surf%qrs(m) = kappa * ( qr(k,j,i) - qr(k-1,j,i) )                                        &
-                        / ( LOG( z_mo / surf%z0q(m) ) - psi_h( z_mo / surf%ol(m) )                 &
-                                                      + psi_h( surf%z0q(m) / surf%ol(m) ) )
+                        / ( surf%ln_z_z0q(m) - psi_h( z_mo / surf%ol(m) )                          &
+                                             + psi_h( surf%z0q(m) / surf%ol(m) ) )
 
           surf%nrs(m) = kappa * ( nr(k,j,i) - nr(k-1,j,i) )                                        &
-                        / ( LOG( z_mo / surf%z0q(m) ) - psi_h( z_mo / surf%ol(m) )                 &
-                                                      + psi_h( surf%z0q(m) / surf%ol(m) ) )
+                        / ( surf%ln_z_z0q(m) - psi_h( z_mo / surf%ol(m) )                          &
+                                             + psi_h( surf%z0q(m) / surf%ol(m) ) )
        ENDDO
 
@@ -1671,6 +1753,6 @@
 
              surf%usws(m) = kappa * ( u(k,j,i) - u(k-1,j,i) )                                      &
-                            / ( LOG( z_mo / surf%z0(m) ) - psi_m( z_mo / surf%ol(m) )              &
-                                                         + psi_m( surf%z0(m) / surf%ol(m) ) )
+                            / ( surf%ln_z_z0(m) - psi_m( z_mo / surf%ol(m) )                       &
+                                                + psi_m( surf%z0(m) / surf%ol(m) ) )
 !
 !--          Please note, the computation of usws is not fully accurate. Actually a further
@@ -1684,5 +1766,5 @@
 !--    At downward-facing surfaces
        ELSE
-          !$OMP PARALLEL DO PRIVATE( i, j, k, z_mo )
+          !$OMP PARALLEL DO PRIVATE( i, j, k )
           DO  m = 1, surf%ns
              i = surf%i(m)
@@ -1690,7 +1772,5 @@
              k = surf%k(m)
 
-             z_mo = surf%z_mo(m)
-
-             surf%usws(m) = kappa * u(k,j,i) / LOG( z_mo / surf%z0(m) )
+             surf%usws(m) = kappa * u(k,j,i) / surf%ln_z_z0(m)
              surf%usws(m) = surf%usws(m) * surf%us(m) * rho_air_zw(k)
           ENDDO
@@ -1713,6 +1793,6 @@
 
              surf%vsws(m) = kappa * ( v(k,j,i) - v(k-1,j,i) )                                      &
-                            / ( LOG( z_mo / surf%z0(m) ) - psi_m( z_mo / surf%ol(m) )              &
-                                                         + psi_m( surf%z0(m) / surf%ol(m) ) )
+                            / ( surf%ln_z_z0(m) - psi_m( z_mo / surf%ol(m) )                       &
+                                                + psi_m( surf%z0(m) / surf%ol(m) ) )
 !
 !--          Please note, the computation of vsws is not fully accurate. Actually a further
@@ -1726,5 +1806,5 @@
 !--    Downward-facing surfaces
        ELSE
-          !$OMP PARALLEL DO PRIVATE( i, j, k, z_mo )
+          !$OMP PARALLEL DO PRIVATE( i, j, k )
           DO  m = 1, surf%ns
              i = surf%i(m)
@@ -1732,19 +1812,16 @@
              k = surf%k(m)
 
-             z_mo = surf%z_mo(m)
-
-             surf%vsws(m) = kappa * v(k,j,i) / LOG( z_mo / surf%z0(m) )
+             surf%vsws(m) = kappa * v(k,j,i) / surf%ln_z_z0(m)
              surf%vsws(m) = surf%vsws(m) * surf%us(m) * rho_air_zw(k)
           ENDDO
        ENDIF
 !
-!--    Compute the vertical kinematic heat flux
+!--    Compute the vertical kinematic heat flux. Note, only upward-facing surfaces are considered,
+!--    at downward-facing surfaces the flux is not parametrized with a scaling parameter.
        IF (  .NOT.  constant_heatflux  .AND.  ( ( time_since_reference_point <=  skip_time_do_lsm  &
              .AND.  simulated_time > 0.0_wp )  .OR.  .NOT.  land_surface )  .AND.                  &
              .NOT.  urban_surface  .AND.  .NOT. downward )  THEN
-          !$OMP PARALLEL DO PRIVATE( i, j, k )
-          DO  m = 1, surf%ns
-             i = surf%i(m)
-             j = surf%j(m)
+          !$OMP PARALLEL DO PRIVATE( k )
+          DO  m = 1, surf%ns
              k = surf%k(m)
              surf%shf(m) = -surf%ts(m) * surf%us(m) * rho_air_zw(k-1)
@@ -1757,8 +1834,6 @@
             .OR.  .NOT.  land_surface  )  .AND.  .NOT.  urban_surface  .AND.  .NOT.  downward )    &
             THEN
-          !$OMP PARALLEL DO PRIVATE( i, j, k )
-          DO  m = 1, surf%ns
-             i    = surf%i(m)
-             j    = surf%j(m)
+          !$OMP PARALLEL DO PRIVATE( k )
+          DO  m = 1, surf%ns
              k    = surf%k(m)
              surf%qsws(m) = -surf%qs(m) * surf%us(m) * rho_air_zw(k-1)
@@ -1768,9 +1843,7 @@
 !--    Compute the vertical scalar flux
        IF (  .NOT.  constant_scalarflux  .AND.  passive_scalar  .AND.  .NOT.  downward )  THEN
-          !$OMP PARALLEL DO PRIVATE( i, j )
-          DO  m = 1, surf%ns
-             i = surf%i(m)
-             j = surf%j(m)
-             surf%ssws(m) = -surf%ss(m) * surf%us(m)
+          !$OMP PARALLEL DO
+          DO  m = 1, surf%ns
+             surf%ssws(m) = -surf%ss(m) * surf%us(m) * rho_air_zw(k-1)
           ENDDO
        ENDIF
@@ -1779,9 +1852,7 @@
        DO  lsp = 1, nvar
           IF (  .NOT.  constant_csflux(lsp)  .AND.  air_chemistry  .AND.  .NOT.  downward )  THEN
-             !$OMP PARALLEL DO PRIVATE( i, j )
+             !$OMP PARALLEL DO
              DO  m = 1, surf%ns
-                i = surf%i(m)
-                j = surf%j(m)
-                surf%cssws(lsp,m) = -surf%css(lsp,m) * surf%us(m)
+                surf%cssws(lsp,m) = -surf%css(lsp,m) * surf%us(m) * rho_air_zw(k-1)
              ENDDO
           ENDIF
@@ -1791,11 +1862,8 @@
 !--    Compute (turbulent) fluxes of cloud water content and cloud drop conc.
        IF ( bulk_cloud_model  .AND.  microphysics_morrison  .AND.  .NOT. downward)  THEN
-          !$OMP PARALLEL DO PRIVATE( i, j )
-          DO  m = 1, surf%ns
-             i = surf%i(m)
-             j = surf%j(m)
-
-             surf%qcsws(m) = -surf%qcs(m) * surf%us(m)
-             surf%ncsws(m) = -surf%ncs(m) * surf%us(m)
+          !$OMP PARALLEL DO
+          DO  m = 1, surf%ns
+             surf%qcsws(m) = -surf%qcs(m) * surf%us(m) * rho_air_zw(k-1)
+             surf%ncsws(m) = -surf%ncs(m) * surf%us(m) * rho_air_zw(k-1)
           ENDDO
        ENDIF
@@ -1803,11 +1871,8 @@
 !--    Compute (turbulent) fluxes of rain water content and rain drop conc.
        IF ( bulk_cloud_model  .AND.  microphysics_seifert  .AND.  .NOT. downward)  THEN
-          !$OMP PARALLEL DO PRIVATE( i, j )
-          DO  m = 1, surf%ns
-             i = surf%i(m)
-             j = surf%j(m)
-
-             surf%qrsws(m) = -surf%qrs(m) * surf%us(m)
-             surf%nrsws(m) = -surf%nrs(m) * surf%us(m)
+          !$OMP PARALLEL DO
+          DO  m = 1, surf%ns
+             surf%qrsws(m) = -surf%qrs(m) * surf%us(m) * rho_air_zw(k-1)
+             surf%nrsws(m) = -surf%nrs(m) * surf%us(m) * rho_air_zw(k-1)
           ENDDO
        ENDIF
@@ -1832,4 +1897,5 @@
 !
 !-- Calcuate surface fluxes at vertical surfaces. No stability is considered.
+!-- Further, no density needs to be considered here.
     ELSE
 !
@@ -1841,5 +1907,5 @@
           flag_u = MERGE( 1.0_wp, 0.0_wp, l == 0  .OR.  l == 1 )
           flag_v = MERGE( 1.0_wp, 0.0_wp, l == 2  .OR.  l == 3 )
-          !$OMP PARALLEL  DO PRIVATE( i, j, k, z_mo )
+          !$OMP PARALLEL  DO PRIVATE( i, j, k )
           DO  m = 1, surf%ns
              i = surf%i(m)
@@ -1847,8 +1913,6 @@
              k = surf%k(m)
 
-             z_mo = surf%z_mo(m)
-
              surf%mom_flux_uv(m) = kappa * ( flag_u * u(k,j,i) + flag_v * v(k,j,i) )  /            &
-                                   LOG( z_mo / surf%z0(m) )
+                                   surf%ln_z_z0(m)
 
             surf%mom_flux_uv(m) = - surf%mom_flux_uv(m) * surf%us(m)
@@ -1858,5 +1922,5 @@
 !--    Compute wsus l={0,1} and wsvs l={2,3}
        IF ( mom_w )  THEN
-          !$OMP PARALLEL  DO PRIVATE( i, j, k, z_mo )
+          !$OMP PARALLEL  DO PRIVATE( i, j, k )
           DO  m = 1, surf%ns
              i = surf%i(m)
@@ -1864,7 +1928,5 @@
              k = surf%k(m)
 
-             z_mo = surf%z_mo(m)
-
-             surf%mom_flux_w(m) = kappa * w(k,j,i) / LOG( z_mo / surf%z0(m) )
+             surf%mom_flux_w(m) = kappa * w(k,j,i) / surf%ln_z_z0(m)
 
              surf%mom_flux_w(m) = - surf%mom_flux_w(m) * surf%us(m)
@@ -1907,12 +1969,10 @@
           ENDIF
 
-          !$OMP PARALLEL DO PRIVATE( i, j, dum, z_mo )
+          !$OMP PARALLEL DO PRIVATE( i, j, dum )
           DO  m = 1, surf%ns
              i = surf%i(m)
              j = surf%j(m)
 
-             z_mo = surf%z_mo(m)
-
-             dum = kappa / LOG( z_mo / surf%z0(m) )
+             dum = kappa / surf%ln_z_z0(m)
 !
 !--          usvs (l=0,1) and vsus (l=2,3)
@@ -1946,10 +2006,6 @@
     CHARACTER(LEN = *), INTENT(IN) ::  z_char  !< string identifier to identify z level
 
-    INTEGER(iwp) ::  i  !< grid index x-dimension
-    INTEGER(iwp) ::  j  !< grid index y-dimension
-    INTEGER(iwp) ::  k  !< grid index z-dimension
     INTEGER(iwp) ::  m  !< running index for surface elements
 
-
     SELECT CASE ( z_char)
 
@@ -1957,9 +2013,4 @@
 
           DO  m = 1, surf%ns
-
-             i = surf%i(m)
-             j = surf%j(m)
-             k = surf%k(m)
-
              surf%pt_10cm(m) = surf%pt_surface(m) + surf%ts(m) / kappa                             &
                                * ( LOG( 0.1_wp /  surf%z0h(m) ) - psi_h( 0.1_wp / surf%ol(m) )     &
Index: palm/trunk/SOURCE/surface_mod.f90
===================================================================
--- palm/trunk/SOURCE/surface_mod.f90	(revision 4592)
+++ palm/trunk/SOURCE/surface_mod.f90	(revision 4593)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Add arrays for pre-calculated ln(z/z0)
+!
+! 4586 2020-07-01 16:16:43Z gronemeier
 ! renamed Richardson flux number into gradient Richardson number (1D model)
 !
@@ -248,4 +251,9 @@
 
        REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  css    !< scaling parameter chemical species
+!
+!--    Pre-defined arrays for ln(z/z0)
+       REAL(wp), DIMENSION(:), ALLOCATABLE ::  ln_z_z0  !< ln(z/z0)
+       REAL(wp), DIMENSION(:), ALLOCATABLE ::  ln_z_z0h !< ln(z/z0h)
+       REAL(wp), DIMENSION(:), ALLOCATABLE ::  ln_z_z0q !< ln(z/z0q)
 !
 !--    Define arrays for surface fluxes
@@ -1281,4 +1289,9 @@
     DEALLOCATE ( surfaces%uvw_abs )
 !
+!-- Pre-calculated ln(z/z0)
+    DEALLOCATE ( surfaces%ln_z_z0  )
+    DEALLOCATE ( surfaces%ln_z_z0h )
+    DEALLOCATE ( surfaces%ln_z_z0q )
+!
 !-- Roughness
     DEALLOCATE ( surfaces%z0 )
@@ -1411,4 +1424,9 @@
     ALLOCATE ( surfaces%uvw_abs(1:surfaces%ns) )
 !
+!-- Precalculated ln(z/z0)
+    ALLOCATE( surfaces%ln_z_z0(1:surfaces%ns)  )
+    ALLOCATE( surfaces%ln_z_z0h(1:surfaces%ns) )
+    ALLOCATE( surfaces%ln_z_z0q(1:surfaces%ns) )
+!
 !-- Roughness
     ALLOCATE ( surfaces%z0(1:surfaces%ns)  )
@@ -1520,4 +1538,7 @@
     !$ACC DELETE(surfaces%z_mo(1:surfaces%ns)) &
     !$ACC DELETE(surfaces%uvw_abs(1:surfaces%ns)) &
+    !$ACC DELETE(surfaces%ln_z_z0(1:surfaces%ns)) &
+    !$ACC DELETE(surfaces%ln_z_z0h(1:surfaces%ns)) &
+    !$ACC DELETE(surfaces%ln_z_z0q(1:surfaces%ns)) &
     !$ACC DELETE(surfaces%z0(1:surfaces%ns)) &
     !$ACC COPYOUT(surfaces%us(1:surfaces%ns)) &
@@ -1561,4 +1582,7 @@
     !$ACC COPYIN(surfaces%z_mo(1:surfaces%ns)) &
     !$ACC COPYIN(surfaces%uvw_abs(1:surfaces%ns)) &
+    !$ACC COPYIN(surfaces%ln_z_z0(1:surfaces%ns)) &
+    !$ACC COPYIN(surfaces%ln_z_z0h(1:surfaces%ns)) &
+    !$ACC COPYIN(surfaces%ln_z_z0q(1:surfaces%ns)) &
     !$ACC COPYIN(surfaces%z0(1:surfaces%ns)) &
     !$ACC COPYIN(surfaces%us(1:surfaces%ns)) &
@@ -1830,9 +1854,13 @@
     DEALLOCATE ( surfaces%uvw_abs )
 !
+!-- Precalculated ln(z/z0)
+    DEALLOCATE ( surfaces%ln_z_z0  )
+    DEALLOCATE ( surfaces%ln_z_z0h )
+    DEALLOCATE ( surfaces%ln_z_z0q )
+!
 !-- Roughness
     DEALLOCATE ( surfaces%z0 )
     DEALLOCATE ( surfaces%z0h )
     DEALLOCATE ( surfaces%z0q )
-
 !
 !-- Friction velocity
@@ -1879,6 +1907,6 @@
 !-- Scaling parameter (cs*) and surface flux of chemical species
     IF ( air_chemistry )  THEN
-          DEALLOCATE ( surfaces%css )
-          DEALLOCATE ( surfaces%cssws )
+       DEALLOCATE ( surfaces%css )
+       DEALLOCATE ( surfaces%cssws )
     ENDIF
 !
@@ -1958,4 +1986,9 @@
     ALLOCATE ( surfaces%uvw_abs(1:surfaces%ns) )
 !
+!-- Precalculated ln(z/z0)
+    ALLOCATE( surfaces%ln_z_z0(1:surfaces%ns) )
+    ALLOCATE( surfaces%ln_z_z0h(1:surfaces%ns) )
+    ALLOCATE( surfaces%ln_z_z0q(1:surfaces%ns) )
+!
 !-- Roughness
     ALLOCATE ( surfaces%z0(1:surfaces%ns)  )
@@ -2007,6 +2040,6 @@
 !-- Scaling parameter (cs*) and surface flux of chemical species
     IF ( air_chemistry )  THEN
-          ALLOCATE ( surfaces%css(1:nvar,1:surfaces%ns) )
-          ALLOCATE ( surfaces%cssws(1:nvar,1:surfaces%ns) )
+       ALLOCATE ( surfaces%css(1:nvar,1:surfaces%ns) )
+       ALLOCATE ( surfaces%cssws(1:nvar,1:surfaces%ns) )
     ENDIF
 !
@@ -2062,4 +2095,7 @@
     !$ACC DELETE(surfaces%k(1:surfaces%ns)) &
     !$ACC DELETE(surfaces%uvw_abs(1:surfaces%ns)) &
+    !$ACC DELETE(surfaces%ln_z_z0(1:surfaces%ns) ) &
+    !$ACC DELETE(surfaces%ln_z_z0h(1:surfaces%ns) ) &
+    !$ACC DELETE(surfaces%ln_z_z0q(1:surfaces%ns) ) &
     !$ACC DELETE(surfaces%z0(1:surfaces%ns)) &
     !$ACC DELETE(surfaces%rib(1:surfaces%ns)) &
@@ -2094,4 +2130,7 @@
     !$ACC COPYIN(surfaces%k(1:surfaces%ns)) &
     !$ACC COPYIN(surfaces%uvw_abs(1:surfaces%ns)) &
+    !$ACC COPYIN(surfaces%ln_z_z0(1:surfaces%ns) ) &
+    !$ACC COPYIN(surfaces%ln_z_z0h(1:surfaces%ns) ) &
+    !$ACC COPYIN(surfaces%ln_z_z0q(1:surfaces%ns) ) &
     !$ACC COPYIN(surfaces%z0(1:surfaces%ns)) &
     !$ACC COPYIN(surfaces%rib(1:surfaces%ns)) &
@@ -2569,4 +2608,9 @@
     surf%rib(num_h)     = 0.0_wp
     surf%uvw_abs(num_h) = 0.0_wp
+!
+!-- Initialize ln(z/z0)
+    surf%ln_z_z0(num_h)  = LOG( surf%z_mo(num_h) / surf%z0(num_h) )
+    surf%ln_z_z0h(num_h) = LOG( surf%z_mo(num_h) / surf%z0h(num_h) )
+    surf%ln_z_z0q(num_h) = LOG( surf%z_mo(num_h) / surf%z0q(num_h) )
 
     IF ( .NOT. constant_diffusion )  THEN
@@ -2871,4 +2915,9 @@
 
     surf%us(num_v)  = 0.0_wp
+!
+!-- Initialize ln(z/z0)
+    surf%ln_z_z0(num_v)  = LOG( surf%z_mo(num_v) / surf%z0(num_v)  )
+    surf%ln_z_z0h(num_v) = LOG( surf%z_mo(num_v) / surf%z0h(num_v) )
+    surf%ln_z_z0q(num_v) = LOG( surf%z_mo(num_v) / surf%z0q(num_v) )
 !
 !-- If required, initialize Obukhov length
@@ -5106,5 +5155,5 @@
     ENDIF
     IF ( INDEX( restart_string(1:length), '%cssws' ) /= 0 )  THEN
-       IF ( ALLOCATED( surf_target%cssws )  .AND.  ALLOCATED( surf_file%cssws   ) )  THEN
+       IF ( ALLOCATED( surf_target%cssws )  .AND.  ALLOCATED( surf_file%cssws ) )  THEN
           DO  lsp = 1, nvar
              surf_target%cssws(lsp,m_target) = surf_file%cssws(lsp,m_file)
@@ -5507,10 +5556,9 @@
 !
 !-- Redistribute surface elements on its respective type.
-
     DO  l = 0, 2
        CALL restore_surface_elements( surf_def_h(l), surf_h(l) )
-       CALL restore_surface_elements( surf_lsm_h, surf_h(l) )
-       CALL restore_surface_elements( surf_usm_h, surf_h(l) )
     ENDDO
+    CALL restore_surface_elements( surf_lsm_h, surf_h(0) )
+    CALL restore_surface_elements( surf_usm_h, surf_h(0) )
 
     DO  l = 0, 3
