Index: /palm/trunk/SOURCE/land_surface_model_mod.f90
===================================================================
--- /palm/trunk/SOURCE/land_surface_model_mod.f90	(revision 4716)
+++ /palm/trunk/SOURCE/land_surface_model_mod.f90	(revision 4717)
@@ -25,4 +25,8 @@
 ! -----------------
 ! $Id$
+! Fixes and optimizations of OpenMP parallelization, formatting of OpenMP
+! directives (J. Resler)
+!
+! 4716 2020-09-30 22:06:37Z pavelkrc
 ! Revert change at water surfaces (previous)
 !
@@ -1744,9 +1748,8 @@
     i_off = surf%ioff
 
-    !$OMP PARALLEL PRIVATE (m, i, j, k, lambda_h_sat, ke, lambda_soil, lambda_surface, ueff,       &
-    !$OMP&                  c_surface_tmp, f1,m_total, f2, e_s, e, f3, m_min, m_liq_max, q_s,      &
-    !$OMP&                  f_qsws_veg, f_qsws_soil, f_qsws_liq, f_shf, f_qsws, e_s_dt, dq_s_dt,   &
-    !$OMP&                  coef_1, coef_2, tend)
-    !$OMP DO SCHEDULE (STATIC)
+    !$OMP PARALLEL DO PRIVATE (m, i, j, k, lambda_h_sat, ke, lambda_soil, lambda_surface, ueff,    &
+    !$OMP&                     c_surface_tmp, f1,m_total, f2, e_s, e, f3, m_min, m_liq_max, q_s,   &
+    !$OMP&                     f_qsws_veg, f_qsws_soil, f_qsws_liq, f_shf, f_qsws, e_s_dt, dq_s_dt,&
+    !$OMP&                     coef_1, coef_2, tend) SCHEDULE (STATIC)
     DO  m = 1, surf%ns
 
@@ -2109,5 +2112,5 @@
        surf%shf(m) = - f_shf * ( surf%pt1(m) - surf%pt_surface(m) ) / c_p
 !
-! update the 3d field of rad_lw_out array to have consistent output
+!--    update the 3d field of rad_lw_out array to have consistent output
        IF ( upward )  THEN
           IF ( radiation_scheme == 'rrtmg' ) THEN
@@ -2246,5 +2249,4 @@
 
     ENDDO
-    !$OMP END PARALLEL
 
 !
@@ -2294,7 +2296,5 @@
        REAL(wp) ::  resistance    !< aerodynamic and soil resistance term
 
-
-       !$OMP PARALLEL PRIVATE (m, i, j, k, e_s, q_s, resistance)
-       !$OMP DO SCHEDULE (STATIC)
+       !$OMP PARALLEL DO PRIVATE (m, i, j, k, e_s, q_s, resistance) SCHEDULE (STATIC)
        DO  m = 1, surf%ns
 
@@ -2330,8 +2330,5 @@
                                   ( 1.0_wp + 0.61_wp * surf%q_surface(m) )
 
-
-
        ENDDO
-       !$OMP END PARALLEL
 
     END SUBROUTINE calc_q_surface
@@ -5555,7 +5552,6 @@
        ENDIF
 
-       !$OMP PARALLEL PRIVATE (m, k, lambda_temp, lambda_h_sat, ke, tend, gamma_temp, h_vg,       &
-       !$OMP&                  m_total, root_extr)
-       !$OMP DO SCHEDULE (STATIC)
+       !$OMP PARALLEL DO PRIVATE (m, k, lambda_temp, lambda_h_sat, ke, tend, gamma_temp, h_vg,    &
+       !$OMP&                     m_total, root_extr) SCHEDULE (STATIC)
        DO  m = 1, surf%ns
 
@@ -5827,5 +5823,4 @@
 
        ENDDO
-       !$OMP END PARALLEL
 !
 !--    Debug location message
@@ -7760,5 +7755,6 @@
 !       REAL(wp) :: re_0 !< near-surface roughness Reynolds number
 
-       !$OMP DO PRIVATE (m, i, j) SCHEDULE (STATIC) REDUCTION(.OR.:flag_exceed_z0, flag_exceed_z0h)
+       !$OMP PARALLEL DO PRIVATE (m, i, j) SCHEDULE (STATIC)                                      &
+       !$OMP&            REDUCTION(.OR.:flag_exceed_z0, flag_exceed_z0h)
        DO  m = 1, surf_lsm_h(0)%ns
 !--       only upward facin horizontal surfaces are considered for water surface processing
@@ -7820,5 +7816,4 @@
           ENDIF
        ENDDO
-       !$OMP END DO
 #if defined( __parallel )
        CALL MPI_ALLREDUCE( MPI_IN_PLACE, flag_exceed_z0, 1, MPI_LOGICAL,       &
Index: /palm/trunk/SOURCE/poisfft_mod.f90
===================================================================
--- /palm/trunk/SOURCE/poisfft_mod.f90	(revision 4716)
+++ /palm/trunk/SOURCE/poisfft_mod.f90	(revision 4717)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Formatting of OpenMP directives (J. Resler)
+!
+! 4671 2020-09-09 20:27:58Z pavelkrc
 ! OMP bugfix
 !
@@ -1001,5 +1004,5 @@
 
     tn = 0              ! Default thread number in case of one thread
-!$OMP  PARALLEL DO PRIVATE ( i, j, k, m, n, tn, work_fftx, work_trix )
+    !$OMP  PARALLEL DO PRIVATE ( i, j, k, m, n, tn, work_fftx, work_trix )
     DO  j = nys_x, nyn_x
 
@@ -1142,6 +1145,6 @@
 !
 !--    Code for vector processors
-!$OMP  PARALLEL PRIVATE ( i, j, k )
-!$OMP  DO
+       !$OMP  PARALLEL PRIVATE ( i, j, k )
+       !$OMP  DO
        DO  i = 0, nx
 
@@ -1154,5 +1157,5 @@
        ENDDO
 
-!$OMP  DO
+       !$OMP  DO
        DO  j = nys, nyn
 
@@ -1166,5 +1169,5 @@
 
        ENDDO
-!$OMP  END PARALLEL
+       !$OMP  END PARALLEL
 
     ELSE
@@ -1172,6 +1175,6 @@
 !
 !--    Cache optimized code (there might still be a potential for better optimization).
-!$OMP  PARALLEL PRIVATE (i,j,k)
-!$OMP  DO
+       !$OMP  PARALLEL PRIVATE (i,j,k)
+       !$OMP  DO
        DO  i = 0, nx
 
@@ -1184,5 +1187,5 @@
        ENDDO
 
-!$OMP  DO
+       !$OMP  DO
        DO  j = nys, nyn
           DO  k = 1, nz
@@ -1196,5 +1199,5 @@
 
        ENDDO
-!$OMP  END PARALLEL
+       !$OMP  END PARALLEL
 
     ENDIF
Index: /palm/trunk/SOURCE/pres.f90
===================================================================
--- /palm/trunk/SOURCE/pres.f90	(revision 4716)
+++ /palm/trunk/SOURCE/pres.f90	(revision 4717)
@@ -25,9 +25,12 @@
 ! -----------------
 ! $Id$
+! Fixes and optimizations of OpenMP parallelization, formatting of OpenMP
+! directives (J. Resler)
+!
+! 4651 2020-08-27 07:17:45Z raasch
 ! preprocessor branch for ibm removed
 !
 ! 4649 2020-08-25 12:11:17Z raasch
 ! File re-formatted to follow the PALM coding standard
-!
 !
 ! 4457 2020-03-11 14:20:43Z raasch
@@ -171,5 +174,4 @@
     REAL(wp) ::  ddt_3d            !<
     REAL(wp) ::  d_weight_pres     !<
-    REAL(wp) ::  localsum          !<
     REAL(wp) ::  threadsum         !<
     REAL(wp) ::  weight_pres_l     !<
@@ -363,5 +365,5 @@
 
     IF ( psolver(1:9) == 'multigrid' )  THEN
-       !$OMP PARALLEL DO SCHEDULE( STATIC ) PRIVATE (i,j,k)
+       !$OMP PARALLEL DO PRIVATE (i,j,k) SCHEDULE( STATIC )
        DO  i = nxl-1, nxr+1
           DO  j = nys-1, nyn+1
@@ -372,5 +374,5 @@
        ENDDO
     ELSE
-       !$OMP PARALLEL DO SCHEDULE( STATIC ) PRIVATE (i,j,k)
+       !$OMP PARALLEL DO PRIVATE (i,j,k) SCHEDULE( STATIC )
        !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
        !$ACC PRESENT(d)
@@ -384,9 +386,5 @@
     ENDIF
 
-    localsum  = 0.0_wp
-    threadsum = 0.0_wp
-
-    !$OMP PARALLEL PRIVATE (i,j,k)
-    !$OMP DO SCHEDULE( STATIC )
+    !$OMP PARALLEL DO PRIVATE (i,j,k) SCHEDULE( STATIC )
     !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
     !$ACC PRESENT(u, v, w, rho_air, rho_air_zw, ddzw, wall_flags_total_0) &
@@ -403,5 +401,4 @@
        ENDDO
     ENDDO
-    !$OMP END PARALLEL
 
 !
@@ -410,6 +407,6 @@
     IF ( intermediate_timestep_count == intermediate_timestep_count_max  .OR.                      &
          intermediate_timestep_count == 0 )  THEN
-       !$OMP PARALLEL PRIVATE (i,j,k) FIRSTPRIVATE(threadsum) REDUCTION(+:localsum)
-       !$OMP DO SCHEDULE( STATIC )
+       threadsum = 0.0_wp
+       !$OMP PARALLEL DO PRIVATE (i,j,k) REDUCTION(+:threadsum) SCHEDULE( STATIC )
        !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i,j,k) &
        !$ACC REDUCTION(+:threadsum) COPY(threadsum) &
@@ -422,6 +419,5 @@
           ENDDO
        ENDDO
-       localsum = localsum + threadsum * dt_3d * weight_pres_l
-       !$OMP END PARALLEL
+       threadsum = threadsum + threadsum * dt_3d * weight_pres_l
     ENDIF
 
@@ -430,5 +426,5 @@
     IF ( intermediate_timestep_count == intermediate_timestep_count_max  .OR.                      &
          intermediate_timestep_count == 0 )  THEN
-       sums_divold_l(0:statistic_regions) = localsum
+       sums_divold_l(0:statistic_regions) = threadsum
     ENDIF
 
@@ -445,5 +441,5 @@
 !
 !--    Store computed perturbation pressure and set boundary condition in z-direction
-       !$OMP PARALLEL DO PRIVATE (i,j,k)
+       !$OMP PARALLEL DO PRIVATE (i,j,k) SCHEDULE( STATIC )
        !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
        !$ACC PRESENT(d, tend)
@@ -465,5 +461,5 @@
 !--       natural and urban surfaces
 !--       Upward facing
-          !$OMP PARALLEL DO PRIVATE( i, j, k )
+          !$OMP PARALLEL DO PRIVATE( m, i, j, k ) SCHEDULE( STATIC )
           !$ACC PARALLEL LOOP PRIVATE(i, j, k) &
           !$ACC PRESENT(bc_h, tend)
@@ -476,5 +472,5 @@
 !
 !--       Downward facing
-          !$OMP PARALLEL DO PRIVATE( i, j, k )
+          !$OMP PARALLEL DO PRIVATE( m, i, j, k ) SCHEDULE( STATIC )
           !$ACC PARALLEL LOOP PRIVATE(i, j, k) &
           !$ACC PRESENT(bc_h, tend)
@@ -491,5 +487,5 @@
 !--       urban surfaces
 !--       Upward facing
-          !$OMP PARALLEL DO PRIVATE( i, j, k )
+          !$OMP PARALLEL DO PRIVATE( m, i, j, k ) SCHEDULE( STATIC )
           DO  m = 1, bc_h(0)%ns
              i = bc_h(0)%i(m)
@@ -500,5 +496,5 @@
 !
 !--       Downward facing
-          !$OMP PARALLEL DO PRIVATE( i, j, k )
+          !$OMP PARALLEL DO PRIVATE( m, i, j, k ) SCHEDULE( STATIC )
           DO  m = 1, bc_h(1)%ns
              i = bc_h(1)%i(m)
@@ -515,5 +511,5 @@
 !
 !--       Neumann
-          !$OMP PARALLEL DO PRIVATE (i,j,k)
+          !$OMP PARALLEL DO PRIVATE (i,j) SCHEDULE( STATIC )
           DO  i = nxlg, nxrg
              DO  j = nysg, nyng
@@ -525,5 +521,5 @@
 !
 !--       Dirichlet
-          !$OMP PARALLEL DO PRIVATE (i,j,k)
+          !$OMP PARALLEL DO PRIVATE (i,j) SCHEDULE( STATIC )
           !$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i, j) &
           !$ACC PRESENT(tend)
@@ -590,6 +586,5 @@
 !-- Ghost layers are added in the output routines (except sor-method: see below)
     IF ( intermediate_timestep_count <= 1 )  THEN
-       !$OMP PARALLEL PRIVATE (i,j,k)
-       !$OMP DO
+       !$OMP PARALLEL DO PRIVATE (i,j,k) SCHEDULE( STATIC )
        !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
        !$ACC PRESENT(p, tend)
@@ -601,9 +596,7 @@
           ENDDO
        ENDDO
-       !$OMP END PARALLEL
 
     ELSEIF ( intermediate_timestep_count > 1 )  THEN
-       !$OMP PARALLEL PRIVATE (i,j,k)
-       !$OMP DO
+       !$OMP PARALLEL DO PRIVATE (i,j,k) SCHEDULE( STATIC )
        !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
        !$ACC PRESENT(p, tend)
@@ -615,5 +608,4 @@
           ENDDO
        ENDDO
-       !$OMP END PARALLEL
 
     ENDIF
@@ -635,6 +627,5 @@
 !-- the velocities, zero-gradient conditions for the pressure are set, so that no modification is
 !-- imposed at the boundaries.
-    !$OMP PARALLEL PRIVATE (i,j,k)
-    !$OMP DO
+    !$OMP PARALLEL DO PRIVATE (i,j,k) SCHEDULE( STATIC )
     !$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i, j, k) &
     !$ACC PRESENT(u, v, w, tend, ddzu, wall_flags_total_0)
@@ -660,5 +651,4 @@
        ENDDO
     ENDDO
-    !$OMP END PARALLEL
 
 !
@@ -675,15 +665,11 @@
     IF ( conserve_volume_flow  .AND.  bc_lr_cyc  .AND.  bc_ns_cyc  .AND.  nxr == nx )  THEN
 
-       !$OMP PARALLEL PRIVATE (j,k)
-       !$OMP DO
+       !$OMP PARALLEL DO PRIVATE (j,k) REDUCTION (volume_flow_l(1))
        DO  j = nys, nyn
-          !$OMP CRITICAL
           DO  k = nzb+1, nzt
              volume_flow_l(1) = volume_flow_l(1) + u(k,j,nxr) * dzw(k)                             &
                                 * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,nxr), 1 ) )
           ENDDO
-          !$OMP END CRITICAL
-       ENDDO
-       !$OMP END PARALLEL
+       ENDDO
 
     ENDIF
@@ -691,15 +677,11 @@
     IF ( conserve_volume_flow  .AND.  bc_ns_cyc  .AND.  bc_lr_cyc  .AND. nyn == ny )  THEN
 
-       !$OMP PARALLEL PRIVATE (i,k)
-       !$OMP DO
-       DO  i = nxl, nxr
-          !$OMP CRITICAL
+       !$OMP PARALLEL DO PRIVATE (i,k) REDUCTION (volume_flow_l(2))
+       DO  i = nxl, nxr
           DO  k = nzb+1, nzt
              volume_flow_l(2) = volume_flow_l(2) + v(k,nyn,i) * dzw(k)                             &
                                 * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,nyn,i), 2 ) )
            ENDDO
-          !$OMP END CRITICAL
-       ENDDO
-       !$OMP END PARALLEL
+       ENDDO
 
     ENDIF
@@ -719,6 +701,5 @@
                                  / volume_flow_area(1:2)
 
-       !$OMP PARALLEL PRIVATE (i,j,k)
-       !$OMP DO
+       !$OMP PARALLEL DO PRIVATE (i,j,k) SCHEDULE( STATIC )
        DO  i = nxl, nxr
           DO  j = nys, nyn
@@ -733,6 +714,4 @@
           ENDDO
        ENDDO
-
-       !$OMP END PARALLEL
 
     ENDIF
@@ -757,9 +736,5 @@
        IF ( topography /= 'flat' )  d = 0.0_wp
 
-       localsum  = 0.0_wp
-       threadsum = 0.0_wp
-
-       !$OMP PARALLEL PRIVATE (i,j,k) FIRSTPRIVATE(threadsum) REDUCTION(+:localsum)
-       !$OMP DO SCHEDULE( STATIC )
+       !$OMP PARALLEL DO PRIVATE (i,j,k) SCHEDULE( STATIC )
        !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
        !$ACC PRESENT(u, v, w, rho_air, rho_air_zw, ddzw, wall_flags_total_0) &
@@ -778,5 +753,6 @@
 !
 !--    Compute possible PE-sum of divergences for flow_statistics
-       !$OMP DO SCHEDULE( STATIC )
+       threadsum = 0.0_wp
+       !$OMP PARALLEL DO PRIVATE (i,j,k) REDUCTION(+:threadsum) SCHEDULE( STATIC )
        !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
        !$ACC REDUCTION(+:threadsum) COPY(threadsum) &
@@ -790,11 +766,8 @@
        ENDDO
 
-       localsum = localsum + threadsum
-       !$OMP END PARALLEL
-
 !
 !--    For completeness, set the divergence sum of all statistic regions to those of the total
 !--    domain
-       sums_divnew_l(0:statistic_regions) = localsum
+       sums_divnew_l(0:statistic_regions) = threadsum
 
        CALL cpu_log( log_point_s(1), 'divergence', 'stop' )
Index: /palm/trunk/SOURCE/prognostic_equations.f90
===================================================================
--- /palm/trunk/SOURCE/prognostic_equations.f90	(revision 4716)
+++ /palm/trunk/SOURCE/prognostic_equations.f90	(revision 4717)
@@ -25,9 +25,12 @@
 ! -----------------
 ! $Id$
+! Fixes and optimizations of OpenMP parallelization, formatting of OpenMP
+! directives (J. Resler)
+!
+! 4671 2020-09-09 20:27:58Z pavelkrc
 ! Implementation of downward facing USM and LSM surfaces
 !
 ! 4649 2020-08-25 12:11:17Z raasch
 ! File re-formatted to follow the PALM coding standard
-!
 !
 ! 4370 2020-01-10 14:00:44Z raasch
@@ -397,6 +400,5 @@
     CALL cpu_log( log_point(32), 'all progn.equations', 'start' )
 
-    !$OMP PARALLEL PRIVATE (i,j)
-    !$OMP DO
+    !$OMP PARALLEL DO PRIVATE (i,j) SCHEDULE( STATIC )
     DO  i = nxl, nxr
        DO  j = nys, nyn
@@ -409,5 +411,4 @@
 !-- Module Inferface for exchange horiz after non_advective_processes but before advection.
 !-- Therefore, non_advective_processes must not run for ghost points.
-    !$OMP END PARALLEL
     CALL module_interface_exchange_horiz()
 !
Index: /palm/trunk/SOURCE/radiation_model_mod.f90
===================================================================
--- /palm/trunk/SOURCE/radiation_model_mod.f90	(revision 4716)
+++ /palm/trunk/SOURCE/radiation_model_mod.f90	(revision 4717)
@@ -28,4 +28,8 @@
 ! -----------------
 ! $Id$
+! Fixes and optimizations of OpenMP parallelization, formatting of OpenMP
+! directives (J. Resler)
+!
+! 4713 2020-09-29 12:02:05Z pavelkrc
 ! Correct OpenMP parallelization including cycles with cumulative variables (J. Resler)
 !
@@ -5954,4 +5958,5 @@
      REAL(wp)                          ::  asrc               !< area of source face
      REAL(wp)                          ::  pcrad              !< irradiance from plant canopy
+     REAL(wp)                          ::  temp               !< temporary variable for calculation
 !-   variables for coupling the radiation modle (e.g. RRTMG) and RTM
      REAL(wp)                          ::  pabsswl            !< total absorbed SW radiation energy in local processor (W)
@@ -6164,5 +6169,5 @@
 
      IF ( surface_reflections)  THEN
-        !$OMP DO PRIVATE (i, j, k, isvf, isurf, isurfsrc) SCHEDULE (STATIC)
+        !$OMP PARALLEL DO PRIVATE (i, j, k, isvf, isurf, isurfsrc, temp) SCHEDULE (STATIC)
         DO  isvf = 1, nsvfl
            isurf = svfsurf(1, isvf)
@@ -6174,14 +6179,15 @@
 !--        For surface-to-surface factors we calculate thermal radiation in 1st pass
            IF ( plant_lw_interact )  THEN
-              surfinl(isurf) = surfinl(isurf) + svf(1,isvf) * svf(2,isvf) * surfoutl(isurfsrc)
+              temp = svf(1,isvf) * svf(2,isvf) * surfoutl(isurfsrc)
            ELSE
-              surfinl(isurf) = surfinl(isurf) + svf(1,isvf) * surfoutl(isurfsrc)
+              temp = svf(1,isvf) * surfoutl(isurfsrc)
            ENDIF
+           !$OMP ATOMIC
+           surfinl(isurf) = surfinl(isurf) + temp
         ENDDO
-        !$OMP END DO
      ENDIF
 !
 !--  diffuse radiation using sky view factor
-     !$OMP DO PRIVATE (i, j, d, isurf) REDUCTION(+:pinswl, pinlwl) SCHEDULE (STATIC)
+     !$OMP PARALLEL DO PRIVATE (i, j, d, isurf) REDUCTION(+:pinswl, pinlwl) SCHEDULE (STATIC)
      DO isurf = 1, nsurfl
         j = surfl(iy, isurf)
@@ -6193,16 +6199,13 @@
         IF ( plant_lw_interact )  THEN
            surfinlwdif(isurf) = rad_lw_in_diff(j,i) * skyvft(isurf)
-!-         update received LW energy for RTM coupling
-           pinlwl = pinlwl + surfinlwdif(isurf) * facearea(d)
         ELSE
            surfinlwdif(isurf) = rad_lw_in_diff(j,i) * skyvf(isurf)
-!-         update received LW energy for RTM coupling
-           pinlwl = pinlwl + surfinlwdif(isurf) * facearea(d)
         ENDIF
+!-      update received LW energy for RTM coupling
+        pinlwl = pinlwl + surfinlwdif(isurf) * facearea(d)
      ENDDO
-     !$OMP END DO
 !
 !--  MRT diffuse irradiance
-     !$OMP DO PRIVATE (i, j, imrt) SCHEDULE (STATIC)
+     !$OMP PARALLEL DO PRIVATE (i, j, imrt) SCHEDULE (STATIC)
      DO  imrt = 1, nmrtbl
         j = mrtbl(iy, imrt)
@@ -6211,5 +6214,4 @@
         mrtinlw(imrt) = mrtsky(imrt) * rad_lw_in_diff(j,i)
      ENDDO
-     !$OMP END DO
 !
 !--  Direct radiation
@@ -6229,5 +6231,5 @@
         isd = dsidir_rev(j, i)
 !-- TODO: check if isd = -1 to report that this solar position is not precalculated
-        !$OMP DO PRIVATE (i, j, d, isurf)  REDUCTION(+:pinswl) SCHEDULE (STATIC)
+        !$OMP PARALLEL DO PRIVATE (i, j, d, isurf)  REDUCTION(+:pinswl) SCHEDULE (STATIC)
         DO isurf = 1, nsurfl
            j = surfl(iy, isurf)
@@ -6236,11 +6238,10 @@
            surfinswdir(isurf) = rad_sw_in_dir(j,i) *                        &
                 costheta(surfl(id, isurf)) * dsitrans(isurf, isd) * sun_direct_factor
-!-        update received SW energy for RTM coupling
+!--        update received SW energy for RTM coupling
            pinswl = pinswl + surfinswdir(isurf) * facearea(d)
         ENDDO
-        !$OMP END DO
 !
 !--     MRT direct irradiance
-        !$OMP DO PRIVATE (i, j, imrt) SCHEDULE (STATIC)
+        !$OMP PARALLEL DO PRIVATE (i, j, imrt) SCHEDULE (STATIC)
         DO  imrt = 1, nmrtbl
            j = mrtbl(iy, imrt)
@@ -6249,15 +6250,15 @@
                                      * sun_direct_factor / 4.0_wp ! normal to sphere
         ENDDO
-        !$OMP END DO
      ENDIF
 !
 !--  MRT first pass thermal
-     !$OMP DO PRIVATE (imrtf, imrt, isurfsrc) SCHEDULE (STATIC)
+     !$OMP PARALLEL DO PRIVATE (imrtf, imrt, isurfsrc, temp) SCHEDULE (STATIC)
      DO  imrtf = 1, nmrtf
         imrt = mrtfsurf(1, imrtf)
         isurfsrc = mrtfsurf(2, imrtf)
-        mrtinlw(imrt) = mrtinlw(imrt) + mrtf(imrtf) * surfoutl(isurfsrc)
+        temp = mrtf(imrtf) * surfoutl(isurfsrc)
+        !$OMP ATOMIC
+        mrtinlw(imrt) = mrtinlw(imrt) + temp
      ENDDO
-     !$OMP END DO
 !
 !--  Absorption in each local plant canopy grid box from the first atmospheric
@@ -6269,5 +6270,5 @@
          pcbinlw(:) = 0.0_wp
 
-         !$OMP DO PRIVATE (icsf, ipcgb, i, j, k, kk, isurfsrc, pc_abs_frac, pcrad, asrc)          &
+         !$OMP PARALLEL DO PRIVATE (icsf, ipcgb, i, j, k, kk, isurfsrc, pc_abs_frac, pcrad, asrc) &
          !$OMP&   REDUCTION(+:pinswl, pinlwl, pabslwl, pemitlwl, pabs_pc_lwdifl) SCHEDULE (STATIC)
          DO icsf = 1, ncsfl
@@ -6328,5 +6329,4 @@
              ENDIF
          ENDDO
-         !$OMP END DO
 
          pcbinsw(:) = pcbinswdir(:) + pcbinswdif(:)
@@ -6424,5 +6424,5 @@
 !
 !--      Reflected radiation
-         !$OMP DO PRIVATE (isvf, isurf, isurfsrc) SCHEDULE (STATIC)
+         !$OMP PARALLEL DO PRIVATE (isvf, isurf, isurfsrc) SCHEDULE (STATIC)
          DO isvf = 1, nsvfl
              isurf = svfsurf(1, isvf)
@@ -6435,5 +6435,4 @@
              ENDIF
          ENDDO
-         !$OMP END DO
 !
 !--      NOTE: PC absorbtion and MRT from reflected can both be done at once
@@ -6443,5 +6442,5 @@
 !
 !--      Radiation absorbed by plant canopy
-         !$OMP DO PRIVATE (icsf, ipcgb, isurfsrc, asrc) SCHEDULE (STATIC)
+         !$OMP PARALLEL DO PRIVATE (icsf, ipcgb, isurfsrc, asrc, temp) SCHEDULE (STATIC)
          DO  icsf = 1, ncsfl
              ipcgb = csfsurf(1, icsf)
@@ -6453,20 +6452,26 @@
 !--          stored within `csf'
              asrc = facearea(surf(id, isurfsrc))
-             pcbinsw(ipcgb) = pcbinsw(ipcgb) + csf(1,icsf) * surfouts(isurfsrc) * asrc
+             temp = csf(1,icsf) * surfouts(isurfsrc) * asrc
+             !$OMP ATOMIC
+             pcbinsw(ipcgb) = pcbinsw(ipcgb) + temp
              IF ( plant_lw_interact )  THEN
-                pcbinlw(ipcgb) = pcbinlw(ipcgb) + csf(1,icsf) * surfoutl(isurfsrc) * asrc
+                temp = csf(1,icsf) * surfoutl(isurfsrc) * asrc
+                !$OMP ATOMIC
+                pcbinlw(ipcgb) = pcbinlw(ipcgb) + temp
              ENDIF
          ENDDO
-         !$OMP END DO
 !
 !--      MRT reflected
-         !$OMP DO PRIVATE (imrtf, imrt, isurfsrc) SCHEDULE (STATIC)
+         !$OMP PARALLEL DO PRIVATE (imrtf, imrt, isurfsrc, temp) SCHEDULE (STATIC)
          DO  imrtf = 1, nmrtf
             imrt = mrtfsurf(1, imrtf)
             isurfsrc = mrtfsurf(2, imrtf)
-            mrtinsw(imrt) = mrtinsw(imrt) + mrtft(imrtf) * surfouts(isurfsrc)
-            mrtinlw(imrt) = mrtinlw(imrt) + mrtf(imrtf) * surfoutl(isurfsrc)
+            temp = mrtft(imrtf) * surfouts(isurfsrc)
+            !$OMP ATOMIC
+            mrtinsw(imrt) = mrtinsw(imrt) + temp
+            temp = mrtf(imrtf) * surfoutl(isurfsrc)
+            !$OMP ATOMIC
+            mrtinlw(imrt) = mrtinlw(imrt) + temp
          ENDDO
-         !$OMP END DO
 
          IF ( trace_fluxes_above >= 0.0_wp )  THEN
@@ -6490,5 +6495,5 @@
      IF ( npcbl > 0 )  THEN
          pcm_heating_rate(:,:,:) = 0.0_wp
-         !$OMP DO PRIVATE (ipcgb, i, j, k, kk) REDUCTION(+:pabsswl) SCHEDULE (STATIC)
+         !$OMP PARALLEL DO PRIVATE (ipcgb, i, j, k, kk) REDUCTION(+:pabsswl) SCHEDULE (STATIC)
          DO ipcgb = 1, npcbl
              j = pcbl(iy, ipcgb)
@@ -6503,5 +6508,4 @@
              pabsswl = pabsswl + pcbinsw(ipcgb)
          ENDDO
-         !$OMP END DO
 
          IF ( humidity .AND. plant_canopy_transpiration ) THEN
@@ -6509,5 +6513,5 @@
              pcm_transpiration_rate(:,:,:) = 0.0_wp
              pcm_latent_rate(:,:,:) = 0.0_wp
-             !$OMP DO PRIVATE (ipcgb, i, j, k, kk) SCHEDULE (STATIC)
+             !$OMP PARALLEL DO PRIVATE (ipcgb, i, j, k, kk) SCHEDULE (STATIC)
              DO ipcgb = 1, npcbl
                  i = pcbl(ix, ipcgb)
@@ -6520,5 +6524,4 @@
                                                    pcm_latent_rate(kk,j,i) )
              ENDDO
-             !$OMP END DO
          ENDIF
      ENDIF
@@ -6630,6 +6633,6 @@
      ENDDO
 
-     !$OMP DO PRIVATE (i, d) REDUCTION(+:pabsswl, pabslwl, pemitlwl, pabs_surf_lwdifl) &
-     !$OMP&   SCHEDULE (STATIC)
+     !$OMP PARALLEL DO PRIVATE (i, d) REDUCTION(+:pabsswl, pabslwl, pemitlwl, pabs_surf_lwdifl) &
+     !$OMP&            SCHEDULE (STATIC)
      DO  i = 1, nsurfl
         d  = surfl(id, i)
@@ -6646,8 +6649,7 @@
                            emiss_surf(i) * facearea(d) * surfinlwdif(i)
      ENDDO
-     !$OMP END DO
 
      DO l = 0, 1
-        !$OMP DO PRIVATE (m) SCHEDULE (STATIC)
+        !$OMP PARALLEL DO PRIVATE (m) SCHEDULE (STATIC)
         DO  m = 1, surf_usm_h(l)%ns
            surf_usm_h(l)%surfhf(m) = surf_usm_h(l)%rad_sw_in(m)  +          &
@@ -6656,6 +6658,5 @@
                                   surf_usm_h(l)%rad_lw_out(m)
         ENDDO
-        !$OMP END DO
-        !$OMP DO PRIVATE (m) SCHEDULE (STATIC)
+        !$OMP PARALLEL DO PRIVATE (m) SCHEDULE (STATIC)
         DO  m = 1, surf_lsm_h(l)%ns
            surf_lsm_h(l)%surfhf(m) = surf_lsm_h(l)%rad_sw_in(m)  +          &
@@ -6664,10 +6665,9 @@
                                   surf_lsm_h(l)%rad_lw_out(m)
         ENDDO
-        !$OMP END DO
      ENDDO
 
      DO  l = 0, 3
-        !$OMP DO PRIVATE (m) SCHEDULE (STATIC)
 !--     urban
+        !$OMP PARALLEL DO PRIVATE (m) SCHEDULE (STATIC)
         DO  m = 1, surf_usm_v(l)%ns
            surf_usm_v(l)%surfhf(m) = surf_usm_v(l)%rad_sw_in(m)  +          &
@@ -6676,7 +6676,6 @@
                                      surf_usm_v(l)%rad_lw_out(m)
         ENDDO
-        !$OMP END DO
 !--     land
-        !$OMP DO PRIVATE (m) SCHEDULE (STATIC)
+        !$OMP PARALLEL DO PRIVATE (m) SCHEDULE (STATIC)
         DO  m = 1, surf_lsm_v(l)%ns
            surf_lsm_v(l)%surfhf(m) = surf_lsm_v(l)%rad_sw_in(m)  +          &
@@ -6686,5 +6685,4 @@
 
         ENDDO
-        !$OMP END DO
      ENDDO
 !
Index: /palm/trunk/SOURCE/surface_layer_fluxes_mod.f90
===================================================================
--- /palm/trunk/SOURCE/surface_layer_fluxes_mod.f90	(revision 4716)
+++ /palm/trunk/SOURCE/surface_layer_fluxes_mod.f90	(revision 4717)
@@ -25,4 +25,8 @@
 ! -----------------
 ! $Id$
+! Fixes and optimizations of OpenMP parallelization, formatting of OpenMP
+! directives (J. Resler)
+!
+! 4691 2020-09-22 14:38:38Z suehring
 ! Bugfix for commit 4593 in vector branch of calc_ol
 !
@@ -720,5 +724,5 @@
 !-- 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 )
+    !$OMP PARALLEL DO PRIVATE( z_mo )
     !$ACC PARALLEL LOOP PRIVATE(z_mo) &
     !$ACC PRESENT(surf)
@@ -1050,6 +1054,5 @@
 !
 !--    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)
+       !$OMP PARALLEL DO PRIVATE(i, j, z_mo, ol_old, iter, 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) &
@@ -1583,5 +1586,5 @@
 
           IF ( bulk_cloud_model  .OR.  cloud_droplets )  THEN
-            !$OMP PARALLEL DO PRIVATE( i, j, k, z_mo )
+             !$OMP PARALLEL DO PRIVATE( i, j, k, z_mo )
              DO  m = 1, surf%ns
                 i = surf%i(m)
@@ -1952,5 +1955,5 @@
 !--    Compute wsus l={0,1} and wsvs l={2,3}
        IF ( mom_w )  THEN
-          !$OMP PARALLEL  DO PRIVATE( i, j, k )
+          !$OMP PARALLEL DO PRIVATE( i, j, k )
           DO  m = 1, surf%ns
              i = surf%i(m)
Index: /palm/trunk/SOURCE/timestep.f90
===================================================================
--- /palm/trunk/SOURCE/timestep.f90	(revision 4716)
+++ /palm/trunk/SOURCE/timestep.f90	(revision 4717)
@@ -20,14 +20,17 @@
 ! Current revisions:
 ! -----------------
-! 
-! 
+!
+!
 ! Former revisions:
 ! -----------------
 ! $Id$
+! Fixes and optimizations of OpenMP parallelization, formatting of OpenMP
+! directives (J. Resler)
+! 
+! 4564 2020-06-12 14:03:36Z raasch
 ! Vertical nesting method of Huq et al. (2019) removed
 ! 
 ! 4540 2020-05-18 15:23:29Z raasch
 ! File re-formatted to follow the PALM coding standard
-!
 !
 ! 4444 2020-03-05 15:59:50Z raasch
@@ -308,6 +311,5 @@
        ENDDO
 
-       !$OMP PARALLEL private(i,j,k) reduction(MIN: dt_diff_l)
-       !$OMP DO
+       !$OMP PARALLEL DO private(i,j,k) reduction(MIN: dt_diff_l)
        !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i,j,k) &
        !$ACC COPY(dt_diff_l) REDUCTION(MIN: dt_diff_l) &
@@ -321,5 +323,4 @@
           ENDDO
        ENDDO
-       !$OMP END PARALLEL
 #if defined( __parallel )
        IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
Index: /palm/trunk/SOURCE/transpose.f90
===================================================================
--- /palm/trunk/SOURCE/transpose.f90	(revision 4716)
+++ /palm/trunk/SOURCE/transpose.f90	(revision 4717)
@@ -20,11 +20,13 @@
 ! Current revisions:
 ! -----------------
-! 
-! 
+!
+!
 ! Former revisions:
 ! -----------------
 ! $Id$
+! Formatting of OpenMP directives (J. Resler)
+! 
+! 4540 2020-05-18 15:23:29Z raasch
 ! File re-formatted to follow the PALM coding standard
-!
 !
 ! 4429 2020-02-27 15:24:30Z raasch
@@ -197,5 +199,5 @@
 !
 !--    Reorder transposed array
-!$OMP  PARALLEL PRIVATE ( i, j, k, l, ys )
+       !$OMP  PARALLEL PRIVATE ( i, j, k, l, ys )
        DO  l = 0, pdims(2) - 1
           ys = 0 + l * ( nyn_x - nys_x + 1 )
@@ -214,5 +216,5 @@
           !$OMP END DO NOWAIT
        ENDDO
-!$OMP  END PARALLEL
+       !$OMP  END PARALLEL
 #endif
 
@@ -221,6 +223,6 @@
 !
 !--    Reorder transposed array
-!$OMP  PARALLEL PRIVATE ( i, j, k )
-!$OMP  DO
+       !$OMP  PARALLEL PRIVATE ( i, j, k )
+       !$OMP  DO
 #if __acc_fft_device
        !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i,j,k) &
@@ -234,5 +236,5 @@
           ENDDO
        ENDDO
-!$OMP  END PARALLEL
+       !$OMP  END PARALLEL
 
     ENDIF
@@ -271,9 +273,9 @@
 !-- Rearrange indices of input array in order to make data to be send by MPI contiguous.
 !-- In case of parallel fft/transposition, scattered store is faster in backward direction!!!
-   !$OMP  PARALLEL PRIVATE ( i, j, k )
-   !$OMP  DO
-#if __acc_fft_device
-   !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i,j,k) &
-   !$ACC PRESENT(f_out, f_inv)
+    !$OMP  PARALLEL PRIVATE ( i, j, k )
+    !$OMP  DO
+#if __acc_fft_device
+    !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i,j,k) &
+    !$ACC PRESENT(f_out, f_inv)
 #endif
     DO  i = nxl, nxr
@@ -430,6 +432,6 @@
 !
 !--    Reorder the array in a way that the z index is in first position
-!$OMP  PARALLEL PRIVATE ( i, j, k )
-!$OMP  DO
+       !$OMP PARALLEL PRIVATE ( i, j, k )
+       !$OMP DO
 #if __acc_fft_device
        !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i,j,k) &
@@ -443,5 +445,5 @@
           ENDDO
        ENDDO
-!$OMP  END PARALLEL
+       !$OMP END PARALLEL
 
     ENDIF
@@ -481,6 +483,6 @@
 !
 !-- Rearrange indices of input array in order to make data to be send by MPI contiguous.
-   !$OMP  PARALLEL PRIVATE ( i, j, k )
-   !$OMP  DO
+   !$OMP PARALLEL PRIVATE ( i, j, k )
+   !$OMP DO
 #if __acc_fft_device
    !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i,j,k) &
@@ -494,5 +496,5 @@
         ENDDO
     ENDDO
-    !$OMP  END PARALLEL
+    !$OMP END PARALLEL
 
  END SUBROUTINE resort_for_yx
@@ -560,5 +562,5 @@
 !
 !--    Reorder input array for transposition
-!$OMP  PARALLEL PRIVATE ( i, j, k, l, ys )
+       !$OMP PARALLEL PRIVATE ( i, j, k, l, ys )
        DO  l = 0, pdims(2) - 1
           ys = 0 + l * ( nyn_x - nys_x + 1 )
@@ -577,5 +579,5 @@
           !$OMP END DO NOWAIT
        ENDDO
-!$OMP  END PARALLEL
+       !$OMP END PARALLEL
 
 !
@@ -610,6 +612,6 @@
 !
 !--    Reorder array f_in the same way as ALLTOALL did it.
-!$OMP  PARALLEL PRIVATE ( i, j, k )
-!$OMP  DO
+       !$OMP PARALLEL PRIVATE ( i, j, k )
+       !$OMP DO
 #if __acc_fft_device
        !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i,j,k) &
@@ -623,5 +625,5 @@
           ENDDO
        ENDDO
-!$OMP  END PARALLEL
+       !$OMP END PARALLEL
 
     ENDIF
@@ -748,6 +750,6 @@
 !
 !-- Rearrange indices of input array in order to make data to be send by MPI contiguous.
-   !$OMP  PARALLEL PRIVATE ( i, j, k )
-   !$OMP  DO
+    !$OMP PARALLEL PRIVATE ( i, j, k )
+    !$OMP DO
 #if __acc_fft_device
     !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i,j,k) &
@@ -761,5 +763,5 @@
         ENDDO
     ENDDO
-    !$OMP  END PARALLEL
+    !$OMP END PARALLEL
 
  END SUBROUTINE resort_for_yz
@@ -827,6 +829,6 @@
     IF ( pdims(1) == 1 )  THEN
 
-!$OMP  PARALLEL PRIVATE ( i, j, k )
-!$OMP  DO
+       !$OMP PARALLEL PRIVATE ( i, j, k )
+       !$OMP DO
 #if __acc_fft_device
        !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i,j,k) &
@@ -840,5 +842,5 @@
           ENDDO
        ENDDO
-!$OMP  END PARALLEL
+       !$OMP END PARALLEL
 
     ELSE
@@ -873,5 +875,5 @@
 !
 !--    Reorder transposed array
-!$OMP  PARALLEL PRIVATE ( i, j, k, l, zs )
+       !$OMP PARALLEL PRIVATE ( i, j, k, l, zs )
        DO  l = 0, pdims(1) - 1
           zs = 1 + l * ( nzt_y - nzb_y + 1 )
@@ -890,5 +892,5 @@
           !$OMP END DO NOWAIT
        ENDDO
-!$OMP  END PARALLEL
+       !$OMP END PARALLEL
 #endif
 
@@ -927,9 +929,9 @@
 !
 !-- Rearrange indices of input array in order to make data to be send by MPI contiguous.
-   !$OMP  PARALLEL PRIVATE ( i, j, k )
-   !$OMP  DO
-#if __acc_fft_device
-   !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i,j,k) &
-   !$ACC PRESENT(f_in, f_inv)
+    !$OMP PARALLEL PRIVATE ( i, j, k )
+    !$OMP DO
+#if __acc_fft_device
+    !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i,j,k) &
+    !$ACC PRESENT(f_in, f_inv)
 #endif
     DO  i = nxl, nxr
@@ -940,5 +942,5 @@
         ENDDO
     ENDDO
-    !$OMP  END PARALLEL
+    !$OMP END PARALLEL
 
  END SUBROUTINE resort_for_zx
@@ -1016,6 +1018,6 @@
     IF ( pdims(1) == 1 )  THEN
 
-!$OMP  PARALLEL PRIVATE ( i, j, k )
-!$OMP  DO
+       !$OMP PARALLEL PRIVATE ( i, j, k )
+       !$OMP DO
 #if __acc_fft_device
        !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i,j,k) &
@@ -1029,5 +1031,5 @@
           ENDDO
        ENDDO
-!$OMP  END PARALLEL
+       !$OMP END PARALLEL
 
     ELSE
@@ -1080,5 +1082,5 @@
        ELSE
 
-          !$OMP  PARALLEL PRIVATE ( i, j, k, l, xs )
+          !$OMP PARALLEL PRIVATE ( i, j, k, l, xs )
           DO  l = 0, pdims(1) - 1
              xs = 0 + l * nnx
@@ -1097,5 +1099,5 @@
              !$OMP END DO NOWAIT
           ENDDO
-          !$OMP  END PARALLEL
+          !$OMP END PARALLEL
 
        ENDIF
@@ -1139,6 +1141,6 @@
 !
 !-- Rearrange indices of input array in order to make data to be send by MPI contiguous.
-    !$OMP  PARALLEL PRIVATE ( i, j, k )
-    !$OMP  DO
+    !$OMP PARALLEL PRIVATE ( i, j, k )
+    !$OMP DO
 #if __acc_fft_device
     !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i,j,k) &
@@ -1152,5 +1154,5 @@
         ENDDO
     ENDDO
-    !$OMP  END PARALLEL
+    !$OMP END PARALLEL
 
  END SUBROUTINE resort_for_zy
@@ -1220,5 +1222,5 @@
 !
 !--    Reorder input array for transposition
-!$OMP  PARALLEL PRIVATE ( i, j, k, l, zs )
+       !$OMP PARALLEL PRIVATE ( i, j, k, l, zs )
        DO  l = 0, pdims(1) - 1
           zs = 1 + l * ( nzt_y - nzb_y + 1 )
@@ -1237,5 +1239,5 @@
           !$OMP END DO NOWAIT
        ENDDO
-!$OMP  END PARALLEL
+       !$OMP END PARALLEL
 
 !
@@ -1269,6 +1271,6 @@
 !
 !--    Reorder the array in the same way like ALLTOALL did it
-!$OMP  PARALLEL PRIVATE ( i, j, k )
-!$OMP  DO
+       !$OMP PARALLEL PRIVATE ( i, j, k )
+       !$OMP DO
 #if __acc_fft_device
        !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i,j,k) &
@@ -1282,5 +1284,5 @@
           ENDDO
        ENDDO
-!$OMP  END PARALLEL
+       !$OMP END PARALLEL
 
     ENDIF
Index: /palm/trunk/SOURCE/turbulence_closure_mod.f90
===================================================================
--- /palm/trunk/SOURCE/turbulence_closure_mod.f90	(revision 4716)
+++ /palm/trunk/SOURCE/turbulence_closure_mod.f90	(revision 4717)
@@ -25,4 +25,8 @@
 ! -----------------
 ! $Id$
+! Fixes and optimizations of OpenMP parallelization, formatting of OpenMP
+! directives (J. Resler)
+!
+! 4674 2020-09-10 10:36:55Z pavelkrc
 ! Update ACC directives for downward facing USM and LSM surfaces
 !
@@ -2103,5 +2107,5 @@
 !--    Not available in case of non-cyclic boundary conditions.
 !--    Default surfaces, upward-facing
-       !$OMP PARALLEL DO PRIVATE(i,j,k,m)
+       !$OMP PARALLEL DO PRIVATE(i, j, k, m, km_sfc)
        !$ACC PARALLEL LOOP PRIVATE(i, j, k, m, km_sfc) &
        !$ACC PRESENT(surf_def_h(0), u, v, drho_air_zw, zu)
@@ -2132,6 +2136,6 @@
 !
 !--    Default surfaces, downward-facing surfaces
-       !$OMP PARALLEL DO PRIVATE(i,j,k,m)
-       !$ACC PARALLEL LOOP PRIVATE(i, j, k, m, km_sfc) &
+       !$OMP PARALLEL DO PRIVATE(i, j, k, m)
+       !$ACC PARALLEL LOOP PRIVATE(i, j, k, m) &
        !$ACC PRESENT(surf_def_h(1), u, v, drho_air_zw, zu, km)
        DO  m = 1, surf_def_h(1)%ns
@@ -2158,5 +2162,5 @@
 !
 !--    Natural surfaces, upward- and downward facing
-       !$OMP PARALLEL DO PRIVATE(i,j,k,m)
+       !$OMP PARALLEL DO PRIVATE(i, j, k, m, km_sfc)
        !$ACC PARALLEL LOOP PRIVATE(i, j, k, m, km_sfc) &
        !$ACC PRESENT(surf_lsm_h(0), u, v, drho_air_zw, zu)
@@ -2187,6 +2191,6 @@
 !
 !--    Natural surfaces, downward-facing surfaces
-       !$OMP PARALLEL DO PRIVATE(i,j,k,m)
-       !$ACC PARALLEL LOOP PRIVATE(i, j, k, m, km_sfc) &
+       !$OMP PARALLEL DO PRIVATE(i, j, k, m)
+       !$ACC PARALLEL LOOP PRIVATE(i, j, k, m) &
        !$ACC PRESENT(surf_lsm_h(1), u, v, drho_air_zw, zu, km)
        DO  m = 1, surf_lsm_h(1)%ns
@@ -2213,5 +2217,5 @@
 !
 !--    Urban surfaces, upward-facing
-       !$OMP PARALLEL DO PRIVATE(i,j,k,m)
+       !$OMP PARALLEL DO PRIVATE(i, j, k, m, km_sfc)
        !$ACC PARALLEL LOOP PRIVATE(i, j, k, m, km_sfc) &
        !$ACC PRESENT(surf_usm_h(0), u, v, drho_air_zw, zu)
@@ -2242,6 +2246,6 @@
 !
 !--    Urban surfaces, downward-facing surfaces
-       !$OMP PARALLEL DO PRIVATE(i,j,k,m)
-       !$ACC PARALLEL LOOP PRIVATE(i, j, k, m, km_sfc) &
+       !$OMP PARALLEL DO PRIVATE(i, j, k, m)
+       !$ACC PARALLEL LOOP PRIVATE(i, j, k, m) &
        !$ACC PRESENT(surf_usm_h(1), u, v, drho_air_zw, zu, km)
        DO  m = 1, surf_usm_h(1)%ns
@@ -5192,5 +5196,5 @@
                                                 ! Data output
     !$ACC END KERNELS
-!$OMP END PARALLEL
+    !$OMP END PARALLEL
 
  END SUBROUTINE tcm_diffusivities_default
Index: /palm/trunk/SOURCE/urban_surface_mod.f90
===================================================================
--- /palm/trunk/SOURCE/urban_surface_mod.f90	(revision 4716)
+++ /palm/trunk/SOURCE/urban_surface_mod.f90	(revision 4717)
@@ -27,4 +27,8 @@
 ! -----------------
 ! $Id$
+! Fixes and optimizations of OpenMP parallelization, formatting of OpenMP
+! directives (J. Resler)
+!
+! 4713 2020-09-29 12:02:05Z pavelkrc
 ! - Do not change original fractions in USM energy balance
 ! - Correct OpenMP parallelization
@@ -5159,5 +5163,4 @@
     ENDIF
 
-    !$OMP PARALLEL PRIVATE (m, i, j, k, kw, wtend, wintend, win_absorp, wall_mod)
     wall_mod=1.0_wp
     IF ( usm_wall_mod  .AND.  during_spinup )  THEN
@@ -5184,5 +5187,5 @@
 !
 !-- Cycle for all surfaces in given direction
-    !$OMP DO SCHEDULE (STATIC)
+    !$OMP PARALLEL DO PRIVATE (m, i, j, k, kw, wtend, wintend, win_absorp, wall_mod) SCHEDULE (STATIC)
     DO  m = 1, surf%ns
 !
@@ -5194,24 +5197,24 @@
 !--    Prognostic equation for ground/roof temperature t_wall
        wtend(:) = 0.0_wp
-       wtend(nzb_wall) = ( 1.0_wp / surf%rho_c_wall(nzb_wall,m) )                     &
-                          * ( surf%lambda_h(nzb_wall,m) * wall_mod(nzb_wall)          &
-                              * ( t_wall%val(nzb_wall+1,m) - t_wall%val(nzb_wall,m) ) &
-                              * surf%ddz_wall(nzb_wall+1,m)                           &
-                              + surf%frac(m,ind_veg_wall)                             &
-                              / ( surf%frac(m,ind_veg_wall)                           &
-                                  + surf%frac(m,ind_pav_green) )                      &
-                              * surf%wghf_eb(m)                                       &
-                              - surf%frac(m,ind_pav_green)                            &
-                              / ( surf%frac(m,ind_veg_wall)                           &
-                                  + surf%frac(m,ind_pav_green) )                      &
-                              * ( surf%lambda_h_green(nzt_wall,m)                     &
-                              * wall_mod(nzt_wall)                                    &
-                              * surf%ddz_green(nzt_wall,m)                            &
-                              + surf%lambda_h(nzb_wall,m)                             &
-                              * wall_mod(nzb_wall)                                    &
-                              * surf%ddz_wall(nzb_wall,m) )                           &
-                              / ( surf%ddz_green(nzt_wall,m)                          &
-                              + surf%ddz_wall(nzb_wall,m) )                           &
-                              * ( t_wall%val(nzb_wall,m) - t_green%val(nzt_wall,m) )  &
+       wtend(nzb_wall) = ( 1.0_wp / surf%rho_c_wall(nzb_wall,m) )                                 &
+                          * ( surf%lambda_h(nzb_wall,m) * wall_mod(nzb_wall)                      &
+                              * ( t_wall%val(nzb_wall+1,m) - t_wall%val(nzb_wall,m) )             &
+                              * surf%ddz_wall(nzb_wall+1,m)                                       &
+                              + surf%frac(m,ind_veg_wall)                                         &
+                              / ( surf%frac(m,ind_veg_wall)                                       &
+                                  + surf%frac(m,ind_pav_green) )                                  &
+                              * surf%wghf_eb(m)                                                   &
+                              - surf%frac(m,ind_pav_green)                                        &
+                              / ( surf%frac(m,ind_veg_wall)                                       &
+                                  + surf%frac(m,ind_pav_green) )                                  &
+                              * ( surf%lambda_h_green(nzt_wall,m)                                 &
+                              * wall_mod(nzt_wall)                                                &
+                              * surf%ddz_green(nzt_wall,m)                                        &
+                              + surf%lambda_h(nzb_wall,m)                                         &
+                              * wall_mod(nzb_wall)                                                &
+                              * surf%ddz_wall(nzb_wall,m) )                                       &
+                              / ( surf%ddz_green(nzt_wall,m)                                      &
+                              + surf%ddz_wall(nzb_wall,m) )                                       &
+                              * ( t_wall%val(nzb_wall,m) - t_green%val(nzt_wall,m) )              &
                             ) * surf%ddz_wall_stag(nzb_wall,m)
 
@@ -5226,23 +5229,21 @@
 
        DO  kw = nzb_wall+1, nzt_wall-1
-          wtend(kw) = ( 1.0_wp / surf%rho_c_wall(kw,m) )                   &
-                         * ( surf%lambda_h(kw,m) * wall_mod(kw)            &
-                             * ( t_wall%val(kw+1,m) - t_wall%val(kw,m) )   &
-                             * surf%ddz_wall(kw+1,m)                       &
-                             - surf%lambda_h(kw-1,m) * wall_mod(kw-1)      &
-                             * ( t_wall%val(kw,m) - t_wall%val(kw-1,m) )   &
-                             * surf%ddz_wall(kw,m)                         &
-                           ) * surf%ddz_wall_stag(kw,m)
+          wtend(kw) = ( 1.0_wp / surf%rho_c_wall(kw,m) )                                          &
+                      * ( surf%lambda_h(kw,m) * wall_mod(kw)                                      &
+                      * ( t_wall%val(kw+1,m) - t_wall%val(kw,m) )                                 &
+                      * surf%ddz_wall(kw+1,m)                                                     &
+                      - surf%lambda_h(kw-1,m) * wall_mod(kw-1)                                    &
+                      * ( t_wall%val(kw,m) - t_wall%val(kw-1,m) )                                 &
+                      * surf%ddz_wall(kw,m)                                                       &
+                      ) * surf%ddz_wall_stag(kw,m)
        ENDDO
-
-       wtend(nzt_wall) = ( 1.0_wp / surf%rho_c_wall(nzt_wall,m) )                       &
-                            * ( -surf%lambda_h(nzt_wall-1,m) * wall_mod(nzt_wall-1)     &
-                                * ( t_wall%val(nzt_wall,m) - t_wall%val(nzt_wall-1,m) ) &
-                                * surf%ddz_wall(nzt_wall,m)                             &
-                                + surf%iwghf_eb(m)                                      &
-                              ) * surf%ddz_wall_stag(nzt_wall,m)
-
-       t_wall_p%val(nzb_wall:nzt_wall,m) = t_wall%val(nzb_wall:nzt_wall,m) + dt_3d      &
-                                         * ( tsc(2) * wtend(nzb_wall:nzt_wall) + tsc(3) &
+       wtend(nzt_wall) = ( 1.0_wp / surf%rho_c_wall(nzt_wall,m) )                                 &
+                         * ( -surf%lambda_h(nzt_wall-1,m) * wall_mod(nzt_wall-1)                  &
+                             * ( t_wall%val(nzt_wall,m) - t_wall%val(nzt_wall-1,m) )              &
+                              * surf%ddz_wall(nzt_wall,m) + surf%iwghf_eb(m)                      &
+                           ) * surf%ddz_wall_stag(nzt_wall,m)
+
+       t_wall_p%val(nzb_wall:nzt_wall,m) = t_wall%val(nzb_wall:nzt_wall,m) + dt_3d                &
+                                         * ( tsc(2) * wtend(nzb_wall:nzt_wall) + tsc(3)           &
                                              * surf%tt_wall_m(nzb_wall:nzt_wall,m) )
 
@@ -5255,12 +5256,12 @@
 !--       of shortwave radiation into account
           wintend(:) = 0.0_wp
-          wintend(nzb_wall) = ( 1.0_wp / surf%rho_c_window(nzb_wall,m) )                    &
-                                * ( surf%lambda_h_window(nzb_wall,m)                        &
-                                * ( t_window%val(nzb_wall+1,m) - t_window%val(nzb_wall,m) ) &
-                                * surf%ddz_window(nzb_wall+1,m)                             &
-                                + surf%wghf_eb_window(m)                                    &
-                                + surf%rad_sw_in(m)                                         &
-                                * ( 1.0_wp - exp( -win_absorp                               &
-                                    * surf%zw_window(nzb_wall,m) ) )                        &
+          wintend(nzb_wall) = ( 1.0_wp / surf%rho_c_window(nzb_wall,m) )                          &
+                                * ( surf%lambda_h_window(nzb_wall,m)                              &
+                                * ( t_window%val(nzb_wall+1,m) - t_window%val(nzb_wall,m) )       &
+                                * surf%ddz_window(nzb_wall+1,m)                                   &
+                                + surf%wghf_eb_window(m)                                          &
+                                + surf%rad_sw_in(m)                                               &
+                                * ( 1.0_wp - exp( -win_absorp                                     &
+                                    * surf%zw_window(nzb_wall,m) ) )                              &
                               ) * surf%ddz_window_stag(nzb_wall,m)
 
@@ -5272,39 +5273,39 @@
           ENDIF
 
+
           DO  kw = nzb_wall+1, nzt_wall-1
-             wintend(kw) = ( 1.0_wp / surf%rho_c_window(kw,m) )                        &
-                              * ( surf%lambda_h_window(kw,m)                           &
-                                  * ( t_window%val(kw+1,m) - t_window%val(kw,m) )      &
-                                  * surf%ddz_window(kw+1,m)                            &
-                                  - surf%lambda_h_window(kw-1,m)                       &
-                                 * ( t_window%val(kw,m) - t_window%val(kw-1,m) )       &
-                                 * surf%ddz_window(kw,m)                               &
-                                 + surf%rad_sw_in(m)                                   &
-                                 * ( exp( -win_absorp * surf%zw_window(kw-1,m) )       &
-                                     - exp(-win_absorp * surf%zw_window(kw,m) )        &
-                                   )                                                   &
-                                 ) * surf%ddz_window_stag(kw,m)
+             wintend(kw) = ( 1.0_wp / surf%rho_c_window(kw,m) )                                   &
+                           * ( surf%lambda_h_window(kw,m)                                         &
+                           * ( t_window%val(kw+1,m) - t_window%val(kw,m) )                        &
+                           * surf%ddz_window(kw+1,m)                                              &
+                           - surf%lambda_h_window(kw-1,m)                                         &
+                           * ( t_window%val(kw,m) - t_window%val(kw-1,m) )                        &
+                           * surf%ddz_window(kw,m)                                                &
+                           + surf%rad_sw_in(m)                                                    &
+                           * ( exp( -win_absorp * surf%zw_window(kw-1,m) )                        &
+                           - exp(-win_absorp * surf%zw_window(kw,m) )                             &
+                             )                                                                    &
+                             ) * surf%ddz_window_stag(kw,m)
 
           ENDDO
 
-          wintend(nzt_wall) = ( 1.0_wp / surf%rho_c_window(nzt_wall,m) )                     &
-                              * ( -surf%lambda_h_window(nzt_wall-1,m)                        &
-                                 * ( t_window%val(nzt_wall,m) - t_window%val(nzt_wall-1,m) ) &
-                                 * surf%ddz_window(nzt_wall,m)                               &
-                                 + surf%iwghf_eb_window(m)                                   &
-                                 + surf%rad_sw_in(m)                                         &
-                                 * ( exp( -win_absorp                                        &
-                                     * surf%zw_window(nzt_wall-1,m) )                        &
-                                   - exp( -win_absorp                                        &
-                                     * surf%zw_window(nzt_wall,m) )                          &
-                                   )                                                         &
-                                ) * surf%ddz_window_stag(nzt_wall,m)
-
-          t_window_p%val(nzb_wall:nzt_wall,m) = t_window%val(nzb_wall:nzt_wall,m) + dt_3d       &
-                                              * ( tsc(2) * wintend(nzb_wall:nzt_wall) + tsc(3)  &
+          wintend(nzt_wall) = ( 1.0_wp / surf%rho_c_window(nzt_wall,m) )                          &
+                                 * ( -surf%lambda_h_window(nzt_wall-1,m)                          &
+                                 * ( t_window%val(nzt_wall,m) - t_window%val(nzt_wall-1,m) )      &
+                                 * surf%ddz_window(nzt_wall,m)                                    &
+                                 + surf%iwghf_eb_window(m)                                        &
+                                 + surf%rad_sw_in(m)                                              &
+                                 * ( exp( -win_absorp                                             &
+                                 * surf%zw_window(nzt_wall-1,m) )                                 &
+                                 - exp( -win_absorp                                               &
+                                 * surf%zw_window(nzt_wall,m) )                                   &
+                                   )                                                              &
+                                   ) * surf%ddz_window_stag(nzt_wall,m)
+
+          t_window_p%val(nzb_wall:nzt_wall,m) = t_window%val(nzb_wall:nzt_wall,m) + dt_3d         &
+                                              * ( tsc(2) * wintend(nzb_wall:nzt_wall) + tsc(3)    &
                                               * surf%tt_window_m(nzb_wall:nzt_wall,m) )
 
        ENDIF
-
 !
 !--    Calculate t_wall tendencies for the next Runge-Kutta step
@@ -5316,5 +5317,5 @@
            ELSEIF ( intermediate_timestep_count < intermediate_timestep_count_max )  THEN
                DO  kw = nzb_wall, nzt_wall
-                  surf%tt_wall_m(kw,m) = -9.5625_wp * wtend(kw) +                            &
+                  surf%tt_wall_m(kw,m) = -9.5625_wp * wtend(kw) +                                 &
                                                5.3125_wp * surf%tt_wall_m(kw,m)
                ENDDO
@@ -5332,5 +5333,5 @@
               ELSEIF ( intermediate_timestep_count < intermediate_timestep_count_max )  THEN
                   DO  kw = nzb_wall, nzt_wall
-                     surf%tt_window_m(kw,m) = -9.5625_wp * wintend(kw) +                     &
+                     surf%tt_window_m(kw,m) = -9.5625_wp * wintend(kw) +                          &
                                                     5.3125_wp * surf%tt_window_m(kw,m)
                   ENDDO
@@ -5339,5 +5340,4 @@
        ENDIF
     ENDDO
-    !$OMP END PARALLEL
 
     IF ( debug_output_timestep )  THEN
@@ -5403,7 +5403,6 @@
        ENDIF
 
-       !$OMP PARALLEL PRIVATE (m, i, j, k, kw, lambda_h_green_sat, ke, lambda_green_temp, gtend,    &
-       !$OMP&                  tend, h_vg, gamma_green_temp, m_total, root_extr_green)
-       !$OMP DO SCHEDULE (STATIC)
+       !$OMP PARALLEL DO PRIVATE (m, i, j, k, kw, lambda_h_green_sat, ke, lambda_green_temp,      &
+       !$OMP&  gtend, tend, h_vg, gamma_green_temp, m_total, root_extr_green) SCHEDULE (STATIC)
        DO  m = 1, surf%ns
           IF (surf%frac(m,ind_pav_green) > 0.0_wp)  THEN
@@ -5417,16 +5416,16 @@
 !
 !--             Calculate volumetric heat capacity of the soil, taking into account water content
-                surf%rho_c_total_green(kw,m) = (surf%rho_c_green(kw,m)              &
-                                                     * (1.0_wp - swc_sat_h(l)%val(kw,m))                     &
+                surf%rho_c_total_green(kw,m) = (surf%rho_c_green(kw,m)                            &
+                                                     * (1.0_wp - swc_sat_h(l)%val(kw,m))          &
                                                      + rho_c_water * swc_h(l)%val(kw,m))
 
 !
 !--             Calculate soil heat conductivity at the center of the soil layers
-                lambda_h_green_sat = lambda_h_green_sm ** ( 1.0_wp - swc_sat_h(l)%val(kw,m) )                &
+                lambda_h_green_sat = lambda_h_green_sm ** ( 1.0_wp - swc_sat_h(l)%val(kw,m) )     &
                                      * lambda_h_water ** swc_h(l)%val(kw,m)
 
                 ke = 1.0_wp + LOG10( MAX( 0.1_wp,swc_h(l)%val(kw,m) / swc_sat_h(l)%val(kw,m) ) )
 
-                lambda_green_temp(kw) = ke * (lambda_h_green_sat - lambda_h_green_dry)                &
+                lambda_green_temp(kw) = ke * (lambda_h_green_sat - lambda_h_green_dry)            &
                                         + lambda_h_green_dry
 
@@ -5439,5 +5438,5 @@
 !--          For pavement surface, the true pavement depth is considered
              DO  kw = nzb_wall, nzt_wall
-                surf%lambda_h_green(kw,m) = ( lambda_green_temp(kw+1) + lambda_green_temp(kw) ) &
+                surf%lambda_h_green(kw,m) = ( lambda_green_temp(kw+1) + lambda_green_temp(kw) )   &
                                                   * 0.5_wp
              ENDDO
@@ -5447,25 +5446,25 @@
 !--          Prognostic equation for ground/roof temperature t_green_h
              gtend(:) = 0.0_wp
-             gtend(nzb_wall) = ( 1.0_wp / surf%rho_c_total_green(nzb_wall,m) )                &
-                               * ( surf%lambda_h_green(nzb_wall,m)                            &
-                                   * ( t_green_h(l)%val(nzb_wall+1,m)                                         &
-                                   - t_green_h(l)%val(nzb_wall,m) )                                           &
-                                   * surf%ddz_green(nzb_wall+1,m)                             &
-                                   + surf%wghf_eb_green(m)                                    &
+             gtend(nzb_wall) = ( 1.0_wp / surf%rho_c_total_green(nzb_wall,m) )                    &
+                               * ( surf%lambda_h_green(nzb_wall,m)                                &
+                                   * ( t_green_h(l)%val(nzb_wall+1,m)                             &
+                                   - t_green_h(l)%val(nzb_wall,m) )                               &
+                                   * surf%ddz_green(nzb_wall+1,m)                                 &
+                                   + surf%wghf_eb_green(m)                                        &
                                  ) * surf%ddz_green_stag(nzb_wall,m)
 
              DO  kw = nzb_wall+1, nzt_wall
-                gtend(kw) = ( 1.0_wp / surf%rho_c_total_green(kw,m) )                         &
-                            * ( surf%lambda_h_green(kw,m)                                     &
-                                * ( t_green_h(l)%val(kw+1,m) - t_green_h(l)%val(kw,m) )                              &
-                                * surf%ddz_green(kw+1,m)                                      &
-                                - surf%lambda_h_green(kw-1,m)                                 &
-                                * ( t_green_h(l)%val(kw,m) - t_green_h(l)%val(kw-1,m) )                              &
-                                * surf%ddz_green(kw,m)                                        &
+                gtend(kw) = ( 1.0_wp / surf%rho_c_total_green(kw,m) )                             &
+                            * ( surf%lambda_h_green(kw,m)                                         &
+                                * ( t_green_h(l)%val(kw+1,m) - t_green_h(l)%val(kw,m) )           &
+                                * surf%ddz_green(kw+1,m)                                          &
+                                - surf%lambda_h_green(kw-1,m)                                     &
+                                * ( t_green_h(l)%val(kw,m) - t_green_h(l)%val(kw-1,m) )           &
+                                * surf%ddz_green(kw,m)                                            &
                               ) * surf%ddz_green_stag(kw,m)
              ENDDO
 
-             t_green_h_p(l)%val(nzb_wall:nzt_wall,m) = t_green_h(l)%val(nzb_wall:nzt_wall,m)  &
-                            + dt_3d * ( tsc(2) * gtend(nzb_wall:nzt_wall) + tsc(3)            &
+             t_green_h_p(l)%val(nzb_wall:nzt_wall,m) = t_green_h(l)%val(nzb_wall:nzt_wall,m)      &
+                            + dt_3d * ( tsc(2) * gtend(nzb_wall:nzt_wall) + tsc(3)                &
                             * surf%tt_green_m(nzb_wall:nzt_wall,m) )
 
@@ -5480,5 +5479,5 @@
                  ELSEIF ( intermediate_timestep_count < intermediate_timestep_count_max )  THEN
                      DO  kw = nzb_wall, nzt_wall
-                        surf%tt_green_m(kw,m) = -9.5625_wp * gtend(kw) + 5.3125_wp             &
+                        surf%tt_green_m(kw,m) = -9.5625_wp * gtend(kw) + 5.3125_wp                 &
                                                       * surf%tt_green_m(kw,m)
                      ENDDO
@@ -5490,7 +5489,7 @@
 !
 !--             Calculate soil diffusivity at the center of the soil layers
-                lambda_green_temp(kw) = ( - b_ch * surf%gamma_w_green_sat(kw,m) * psi_sat      &
-                                          / swc_sat_h(l)%val(kw,m) )                           &
-                                        * ( MAX( swc_h(l)%val(kw,m), wilt_h(l)%val(kw,m) )     &
+                lambda_green_temp(kw) = ( - b_ch * surf%gamma_w_green_sat(kw,m) * psi_sat          &
+                                          / swc_sat_h(l)%val(kw,m) )                               &
+                                        * ( MAX( swc_h(l)%val(kw,m), wilt_h(l)%val(kw,m) )         &
                                           / swc_sat_h(l)%val(kw,m) )**( b_ch + 2.0_wp )
 
@@ -5500,22 +5499,22 @@
 !
 !--                Calculate the hydraulic conductivity after Van Genuchten (1980)
-                   h_vg = ( ( (swc_res_h(l)%val(kw,m) - swc_sat_h(l)%val(kw,m))                &
-                              / ( swc_res_h(l)%val(kw,m) -                 &
-                              MAX( swc_h(l)%val(kw,m), wilt_h(l)%val(kw,m) ) ) )**                                    &
-                              ( surf%n_vg_green(m) / (surf%n_vg_green(m) - 1.0_wp ) ) &
-                              - 1.0_wp                                                                  &
+                   h_vg = ( ( (swc_res_h(l)%val(kw,m) - swc_sat_h(l)%val(kw,m))                    &
+                              / ( swc_res_h(l)%val(kw,m) -                                         &
+                              MAX( swc_h(l)%val(kw,m), wilt_h(l)%val(kw,m) ) ) )**                 &
+                              ( surf%n_vg_green(m) / (surf%n_vg_green(m) - 1.0_wp ) )              &
+                              - 1.0_wp                                                             &
                           )** ( 1.0_wp / surf%n_vg_green(m) ) / surf%alpha_vg_green(m)
 
 
-                   gamma_green_temp(kw) = surf%gamma_w_green_sat(kw,m)                          &
-                                          * ( ( ( 1.0_wp + ( surf%alpha_vg_green(m) * h_vg )**  &
-                                                surf%n_vg_green(m) )**                          &
-                                                ( 1.0_wp - 1.0_wp / surf%n_vg_green(m) )        &
-                                                - ( surf%alpha_vg_green(m) * h_vg )**           &
-                                                ( surf%n_vg_green(m) - 1.0_wp) )**2             &
-                                            ) / ( ( 1.0_wp + ( surf%alpha_vg_green(m) * h_vg )**&
-                                                  surf%n_vg_green(m) )**                        &
-                                                  ( ( 1.0_wp  - 1.0_wp / surf%n_vg_green(m) )   &
-                                                    *( surf%l_vg_green(m) + 2.0_wp) )           &
+                   gamma_green_temp(kw) = surf%gamma_w_green_sat(kw,m)                             &
+                                          * ( ( ( 1.0_wp + ( surf%alpha_vg_green(m) * h_vg )**     &
+                                                surf%n_vg_green(m) )**                             &
+                                                ( 1.0_wp - 1.0_wp / surf%n_vg_green(m) )           &
+                                                - ( surf%alpha_vg_green(m) * h_vg )**              &
+                                                ( surf%n_vg_green(m) - 1.0_wp) )**2                &
+                                            ) / ( ( 1.0_wp + ( surf%alpha_vg_green(m) * h_vg )**   &
+                                                  surf%n_vg_green(m) )**                           &
+                                                  ( ( 1.0_wp  - 1.0_wp / surf%n_vg_green(m) )      &
+                                                    *( surf%l_vg_green(m) + 2.0_wp) )              &
                                                 )
 
@@ -5523,5 +5522,5 @@
 !--             Parametrization of Clapp & Hornberger
                 ELSE
-                   gamma_green_temp(kw) = surf%gamma_w_green_sat(kw,m) * ( swc_h(l)%val(kw,m)   &
+                   gamma_green_temp(kw) = surf%gamma_w_green_sat(kw,m) * ( swc_h(l)%val(kw,m)      &
                                           / swc_sat_h(l)%val(kw,m) )**( 2.0_wp * b_ch + 3.0_wp )
                 ENDIF
@@ -5538,9 +5537,9 @@
                 DO  kw = nzb_wall, nzt_wall-1
 
-                   surf%lambda_w_green(kw,m) = ( lambda_green_temp(kw+1)                     &
-                                                       + lambda_green_temp(kw) )                      &
+                   surf%lambda_w_green(kw,m) = ( lambda_green_temp(kw+1)                          &
+                                                       + lambda_green_temp(kw) )                  &
                                                      * 0.5_wp
-                   surf%gamma_w_green(kw,m)  = ( gamma_green_temp(kw+1)                      &
-                                                       + gamma_green_temp(kw) )                       &
+                   surf%gamma_w_green(kw,m)  = ( gamma_green_temp(kw+1)                           &
+                                                       + gamma_green_temp(kw) )                   &
                                                      * 0.5_wp
 
@@ -5588,42 +5587,42 @@
                 tend(:) = 0.0_wp
 
-                tend(nzb_wall) = ( surf_usm_h(l)%lambda_w_green(nzb_wall,m)                           &
-                                 * ( swc_h(l)%val(nzb_wall+1,m) - swc_h(l)%val(nzb_wall,m) )          &
-                                 * surf_usm_h(l)%ddz_green(nzb_wall+1,m)                              &
-                                 - surf_usm_h(l)%gamma_w_green(nzb_wall,m)                            &
-                                 - ( root_extr_green(nzb_wall) * surf_usm_h(l)%qsws_veg(m)            &
-!                                    + surf_usm_h(l)%qsws_soil_green(m)                               &
-                                   ) * drho_l_lv )                                                    &
+                tend(nzb_wall) = ( surf_usm_h(l)%lambda_w_green(nzb_wall,m)                       &
+                                 * ( swc_h(l)%val(nzb_wall+1,m) - swc_h(l)%val(nzb_wall,m) )      &
+                                 * surf_usm_h(l)%ddz_green(nzb_wall+1,m)                          &
+                                 - surf_usm_h(l)%gamma_w_green(nzb_wall,m)                        &
+                                 - ( root_extr_green(nzb_wall) * surf_usm_h(l)%qsws_veg(m)        &
+!                                    + surf_usm_h(l)%qsws_soil_green(m)                           &
+                                   ) * drho_l_lv )                                                &
                                  * surf_usm_h(l)%ddz_green_stag(nzb_wall,m)
 
                 DO  kw = nzb_wall+1, nzt_wall-1
-                   tend(kw) = ( surf_usm_h(l)%lambda_w_green(kw,m)                                    &
-                                * ( swc_h(l)%val(kw+1,m) - swc_h(l)%val(kw,m) )                       &
-                                * surf_usm_h(l)%ddz_green(kw+1,m)                                     &
-                                - surf_usm_h(l)%gamma_w_green(kw,m)                                   &
-                                - surf_usm_h(l)%lambda_w_green(kw-1,m)                                &
-                                * ( swc_h(l)%val(kw,m) - swc_h(l)%val(kw-1,m) )                       &
-                                * surf_usm_h(l)%ddz_green(kw,m)                                       &
-                                + surf_usm_h(l)%gamma_w_green(kw-1,m)                                 &
-                                - (root_extr_green(kw)                                                &
-                                * surf_usm_h(l)%qsws_veg(m)                                           &
-                                * drho_l_lv)                                                          &
+                   tend(kw) = ( surf_usm_h(l)%lambda_w_green(kw,m)                                &
+                                * ( swc_h(l)%val(kw+1,m) - swc_h(l)%val(kw,m) )                   &
+                                * surf_usm_h(l)%ddz_green(kw+1,m)                                 &
+                                - surf_usm_h(l)%gamma_w_green(kw,m)                               &
+                                - surf_usm_h(l)%lambda_w_green(kw-1,m)                            &
+                                * ( swc_h(l)%val(kw,m) - swc_h(l)%val(kw-1,m) )                   &
+                                * surf_usm_h(l)%ddz_green(kw,m)                                   &
+                                + surf_usm_h(l)%gamma_w_green(kw-1,m)                             &
+                                - (root_extr_green(kw)                                            &
+                                * surf_usm_h(l)%qsws_veg(m)                                       &
+                                * drho_l_lv)                                                      &
                              ) * surf_usm_h(l)%ddz_green_stag(kw,m)
 
                 ENDDO
-                tend(nzt_wall) = ( - surf_usm_h(l)%gamma_w_green(nzt_wall,m)                          &
-                                   - surf_usm_h(l)%lambda_w_green(nzt_wall-1,m)                       &
-                                   * (swc_h(l)%val(nzt_wall,m)                                        &
-                                   - swc_h(l)%val(nzt_wall-1,m))                                      &
-                                   * surf_usm_h(l)%ddz_green(nzt_wall,m)                              &
-                                   + surf_usm_h(l)%gamma_w_green(nzt_wall-1,m)                        &
-                                   - ( root_extr_green(nzt_wall)                                      &
-                                   * surf_usm_h(l)%qsws_veg(m)                                        &
-                                   * drho_l_lv  )                                                     &
+                tend(nzt_wall) = ( - surf_usm_h(l)%gamma_w_green(nzt_wall,m)                      &
+                                   - surf_usm_h(l)%lambda_w_green(nzt_wall-1,m)                   &
+                                   * (swc_h(l)%val(nzt_wall,m)                                    &
+                                   - swc_h(l)%val(nzt_wall-1,m))                                  &
+                                   * surf_usm_h(l)%ddz_green(nzt_wall,m)                          &
+                                   + surf_usm_h(l)%gamma_w_green(nzt_wall-1,m)                    &
+                                   - ( root_extr_green(nzt_wall)                                  &
+                                   * surf_usm_h(l)%qsws_veg(m)                                    &
+                                   * drho_l_lv  )                                                 &
                                  ) * surf_usm_h(l)%ddz_green_stag(nzt_wall,m)
 
-                swc_h_p(l)%val(nzb_wall:nzt_wall,m) = swc_h(l)%val(nzb_wall:nzt_wall,m) + dt_3d       &
-                                               * ( tsc(2) * tend(:) + tsc(3)                          &
-                                                   * surf_usm_h(l)%tswc_h_m(:,m)                      &
+                swc_h_p(l)%val(nzb_wall:nzt_wall,m) = swc_h(l)%val(nzb_wall:nzt_wall,m) + dt_3d   &
+                                               * ( tsc(2) * tend(:) + tsc(3)                      &
+                                                   * surf_usm_h(l)%tswc_h_m(:,m)                  &
                                                   )
 
@@ -5643,5 +5642,5 @@
                    ELSEIF ( intermediate_timestep_count < intermediate_timestep_count_max )  THEN
                       DO  kw = nzb_wall, nzt_wall
-                         surf_usm_h(l)%tswc_h_m(kw,m) = -9.5625_wp * tend(kw) + 5.3125_wp             &
+                         surf_usm_h(l)%tswc_h_m(kw,m) = -9.5625_wp * tend(kw) + 5.3125_wp         &
                                                      * surf_usm_h(l)%tswc_h_m(kw,m)
                       ENDDO
@@ -5652,5 +5651,4 @@
           ENDIF
        ENDDO
-       !$OMP END PARALLEL
     ELSE
        IF ( horizontal) THEN
@@ -5665,5 +5663,5 @@
           t_green      => t_green_v(l)
        ENDIF
-       !$OMP DO SCHEDULE (STATIC)
+       !$OMP PARALLEL DO PRIVATE (m, i, j, k, kw) SCHEDULE (STATIC)
        DO  m = 1, surf%ns
           IF (surf%frac(m,ind_pav_green) > 0.0_wp)  THEN
@@ -5684,27 +5682,27 @@
 ! !--          Prognostic equation for green temperature t_green_v
 !              gtend(:) = 0.0_wp
-!              gtend(nzb_wall) = (1.0_wp / surf%rho_c_green(nzb_wall,m)) * &
-!                                      ( surf%lambda_h_green(nzb_wall,m) * &
-!                                        ( t_green%val(nzb_wall+1,m)             &
-!                                        - t_green%val(nzb_wall,m) ) *           &
-!                                        surf%ddz_green(nzb_wall+1,m)      &
-!                                      + surf%wghf_eb(m) ) *               &
+!              gtend(nzb_wall) = (1.0_wp / surf%rho_c_green(nzb_wall,m)) *                        &
+!                                      ( surf%lambda_h_green(nzb_wall,m) *                        &
+!                                        ( t_green%val(nzb_wall+1,m)                              &
+!                                        - t_green%val(nzb_wall,m) ) *                            &
+!                                        surf%ddz_green(nzb_wall+1,m)                             &
+!                                      + surf%wghf_eb(m) ) *                                      &
 !                                        surf%ddz_green_stag(nzb_wall,m)
 !
 !              DO  kw = nzb_wall+1, nzt_wall
-!                 gtend(kw) = (1.0_wp / surf%rho_c_green(kw,m))          &
-!                           * (   surf%lambda_h_green(kw,m)              &
-!                             * ( t_green%val(kw+1,m) - t_green%val(kw,m) ) &
-!                             * surf%ddz_green(kw+1,m)                   &
-!                           - surf%lambda_h(kw-1,m)                      &
-!                             * ( t_green%val(kw,m) - t_green%val(kw-1,m) ) &
-!                             * surf%ddz_green(kw,m) )                   &
+!                 gtend(kw) = (1.0_wp / surf%rho_c_green(kw,m))                                   &
+!                           * (   surf%lambda_h_green(kw,m)                                       &
+!                             * ( t_green%val(kw+1,m) - t_green%val(kw,m) )                       &
+!                             * surf%ddz_green(kw+1,m)                                            &
+!                           - surf%lambda_h(kw-1,m)                                               &
+!                             * ( t_green%val(kw,m) - t_green%val(kw-1,m) )                       &
+!                             * surf%ddz_green(kw,m) )                                            &
 !                           * surf%ddz_green_stag(kw,m)
 !              ENDDO
 !
-!              t_green_v_p(l)%val(nzb_wall:nzt_wall,m) =                              &
-!                                   t_green%val(nzb_wall:nzt_wall,m)             &
-!                                 + dt_3d * ( tsc(2)                                &
-!                                 * gtend(nzb_wall:nzt_wall) + tsc(3)               &
+!              t_green_v_p(l)%val(nzb_wall:nzt_wall,m) =                                          &
+!                                   t_green%val(nzb_wall:nzt_wall,m)                              &
+!                                 + dt_3d * ( tsc(2)                                              &
+!                                 * gtend(nzb_wall:nzt_wall) + tsc(3)                             &
 !                                 * surf%tt_green_m(nzb_wall:nzt_wall,m) )
 !
@@ -5716,9 +5714,9 @@
 !                        surf%tt_green_m(kw,m) = gtend(kw)
 !                     ENDDO
-!                  ELSEIF ( intermediate_timestep_count <                           &
+!                  ELSEIF ( intermediate_timestep_count <                                         &
 !                           intermediate_timestep_count_max )  THEN
 !                      DO  kw = nzb_wall, nzt_wall
-!                         surf%tt_green_m(kw,m) =                          &
-!                                     - 9.5625_wp * gtend(kw) +                     &
+!                         surf%tt_green_m(kw,m) =                                                 &
+!                                     - 9.5625_wp * gtend(kw) +                                   &
 !                                       5.3125_wp * surf%tt_green_m(kw,m)
 !                      ENDDO
@@ -6768,11 +6766,10 @@
 !
 !-- First, treat horizontal surface elements
-    !$OMP PARALLEL PRIVATE (m, i, j, k, frac_win, frac_wall, frac_green, lambda_surface, ueff,    &
-    !$OMP&                  lambda_surface_window, lambda_surface_green, qv1, rho_cp, rho_lv,     &
-    !$OMP&                  drho_l_lv, f_shf, f_shf_window, f_shf_green, m_total, f1, f2, e_s, e, &
-    !$OMP&                  f3, f_qsws_veg, q_s, f_qsws_liq, f_qsws, e_s_dt, dq_s_dt, coef_1,     &
-    !$OMP&                  coef_window_1, coef_green_1, coef_2, coef_window_2, coef_green_2,     &
-    !$OMP&                  stend_wall, stend_window, stend_green, tend, m_liq_max)
-    !$OMP DO SCHEDULE (STATIC)
+    !$OMP PARALLEL DO PRIVATE (m, i, j, k, frac_win, frac_wall, frac_green, lambda_surface,       &
+    !$OMP&             lambda_surface_window, lambda_surface_green, ueff, qv1, rho_cp, rho_lv,    &
+    !$OMP&             drho_l_lv, f_shf, f_shf_window, f_shf_green, m_total, f1, f2, e_s, e,      &
+    !$OMP&             f3, f_qsws_veg, q_s, f_qsws_liq, f_qsws, e_s_dt, dq_s_dt, coef_1,          &
+    !$OMP&             coef_window_1, coef_green_1, coef_2, coef_window_2, coef_green_2,          &
+    !$OMP&             stend_wall, stend_window, stend_green, tend, m_liq_max) SCHEDULE (STATIC)
     DO  m = 1, surf%ns
 !
@@ -7200,5 +7197,4 @@
 
     ENDDO
-    !$OMP END PARALLEL
 
 !
