Index: palm/trunk/SOURCE/urban_surface_mod.f90
===================================================================
--- palm/trunk/SOURCE/urban_surface_mod.f90	(revision 4701)
+++ palm/trunk/SOURCE/urban_surface_mod.f90	(revision 4713)
@@ -27,4 +27,9 @@
 ! -----------------
 ! $Id$
+! - Do not change original fractions in USM energy balance
+! - Correct OpenMP parallelization
+! Author: J. Resler
+!
+! 4701 2020-09-27 11:02:15Z maronga
 ! Corrected parameter 33 for building_type 2 (ground floor window emissivity
 !
@@ -6763,10 +6768,10 @@
 !
 !-- First, treat horizontal surface elements
-    !$OMP PARALLEL PRIVATE (m, i, j, k, lambda_surface, lambda_surface_window,                      &
-    !$OMP&                  lambda_surface_green, qv1, rho_cp, rho_lv, drho_l_lv, f_shf,            &
-    !$OMP&                  f_shf_window, f_shf_green, m_total, f1, f2, e_s, e, f3, f_qsws_veg,     &
-    !$OMP&                  q_s, f_qsws_liq, f_qsws, e_s_dt, dq_s_dt, coef_1, coef_window_1,        &
-    !$OMP&                  coef_green_1, coef_2, coef_window_2, coef_green_2, stend_wall,          &
-    !$OMP&                  stend_window, stend_green, tend, m_liq_max)
+    !$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)
     DO  m = 1, surf%ns
@@ -6775,10 +6780,11 @@
 !--    Note, this is a temporary fix and needs to be removed later.
        IF ( during_spinup )  THEN
-          frac_win                         = surf%frac(m,ind_wat_win)
-          frac_wall                        = surf%frac(m,ind_veg_wall)
-          frac_green                       = surf%frac(m,ind_pav_green)
-          surf%frac(m,ind_wat_win)   = 0.0_wp
-          surf%frac(m,ind_veg_wall)  = 1.0_wp
-          surf%frac(m,ind_pav_green) = 0.0_wp
+          frac_win   = 0.0_wp
+          frac_wall  = 1.0_wp
+          frac_green = 0.0_wp
+       ELSE
+          frac_win   = surf%frac(m,ind_wat_win)
+          frac_wall  = surf%frac(m,ind_veg_wall)
+          frac_green = surf%frac(m,ind_pav_green)
        ENDIF
 !
@@ -6803,5 +6809,5 @@
        rho_cp  = c_p * hyp(k) / ( r_d * surf%pt1(m) * exner(k) )
 
-       IF ( surf%frac(m,ind_pav_green) > 0.0_wp )  THEN
+       IF ( frac_green > 0.0_wp )  THEN
 !
 !--       Calculate frequently used parameters
@@ -6843,13 +6849,13 @@
 !--           to calculation of roughness relative to concrete. Note, wind velocity is limited
 !--           to avoid division by zero. The nominator can become <= 0.0 for values z0 < 3*10E-4.
-              ueff        = MAX ( SQRT( ( ( u(k,j,i) + u(k,j,i+1) ) * 0.5_wp )**2 +                  &
-                                        ( ( v(k,j,i) + v(k,j+1,i) ) * 0.5_wp )**2 +                  &
-                                        ( ( w(k,j,i) + w(k-1,j,i) ) * 0.5_wp )**2 ),                 &
-                                  ( ( 4.0_wp + 0.1_wp )                                              &
-                                    / ( surf%z0(m) * d_roughness_concrete )                 &
-                                   - 11.8_wp ) / 4.2_wp                                              &
+              ueff        = MAX ( SQRT( ( ( u(k,j,i) + u(k,j,i+1) ) * 0.5_wp )**2 +               &
+                                        ( ( v(k,j,i) + v(k,j+1,i) ) * 0.5_wp )**2 +               &
+                                        ( ( w(k,j,i) + w(k-1,j,i) ) * 0.5_wp )**2 ),              &
+                                  ( ( 4.0_wp + 0.1_wp )                                           &
+                                    / ( surf%z0(m) * d_roughness_concrete )                       &
+                                   - 11.8_wp ) / 4.2_wp                                           &
                                  )
 
-              surf%r_a(m) = rho_cp / ( surf%z0(m) * d_roughness_concrete            &
+              surf%r_a(m) = rho_cp / ( surf%z0(m) * d_roughness_concrete                          &
                                      * ( 11.8_wp + 4.2_wp * ueff )  - 4.0_wp  )
        ENDIF
@@ -6869,5 +6875,5 @@
        f_shf_green   = rho_cp / surf%r_a_green(m)
 
-       IF ( surf%frac(m,ind_pav_green) > 0.0_wp )  THEN
+       IF ( frac_green > 0.0_wp )  THEN
 !
 !--       Adapted from LSM:
@@ -6876,5 +6882,5 @@
 
 !--       f1: Correction for incoming shortwave radiation (stomata close at night)
-          f1 = MIN( 1.0_wp, ( 0.004_wp * surf%rad_sw_in(m) + 0.05_wp ) /                     &
+          f1 = MIN( 1.0_wp, ( 0.004_wp * surf%rad_sw_in(m) + 0.05_wp ) /                          &
                     (0.81_wp * ( 0.004_wp * surf%rad_sw_in(m) + 1.0_wp ) ) )
 !
@@ -6914,5 +6920,5 @@
 !--       obsolete, as r_canopy is not used below.
 !--       To do: check for very dry soil -> r_canopy goes to infinity
-          surf%r_canopy(m) = surf%r_canopy_min(m) /                                    &
+          surf%r_canopy(m) = surf%r_canopy_min(m) /                                               &
                                   ( surf%lai(m) * f1 * f2 * f3 + 1.0E-20_wp )
 !
@@ -6968,5 +6974,5 @@
                + f_shf * surf%pt1(m) + lambda_surface * t_wall%val(nzb_wall,m)
 
-       IF ( ( .NOT. during_spinup ) .AND. (surf%frac(m,ind_wat_win) > 0.0_wp ) )  THEN
+       IF ( ( .NOT. during_spinup ) .AND. (frac_win > 0.0_wp ) )  THEN
           coef_window_1 = surf%rad_net_l(m) +  ( 3.0_wp + 1.0_wp )                                &
                           * surf%emissivity(m,ind_wat_win) * sigma_sb                             &
@@ -6974,5 +6980,5 @@
                           + lambda_surface_window * t_window%val(nzb_wall,m)
        ENDIF
-       IF ( ( humidity ) .AND. ( surf%frac(m,ind_pav_green) > 0.0_wp ) )  THEN
+       IF ( ( humidity ) .AND. ( frac_green > 0.0_wp ) )  THEN
           coef_green_1 = surf%rad_net_l(m) + ( 3.0_wp + 1.0_wp )                                  &
                          * surf%emissivity(m,ind_pav_green) * sigma_sb                            &
@@ -6990,9 +6996,9 @@
        coef_2 = 4.0_wp * surf%emissivity(m,ind_veg_wall) * sigma_sb * t_surf_wall%val(m)**3       &
                 + lambda_surface + f_shf / exner(k)
-       IF ( ( .NOT. during_spinup ) .AND. ( surf%frac(m,ind_wat_win) > 0.0_wp ) )  THEN
+       IF ( ( .NOT. during_spinup ) .AND. ( frac_win > 0.0_wp ) )  THEN
           coef_window_2 = 4.0_wp * surf%emissivity(m,ind_wat_win) * sigma_sb *                    &
                           t_surf_window%val(m)**3 + lambda_surface_window + f_shf_window / exner(k)
        ENDIF
-       IF ( ( humidity ) .AND. ( surf%frac(m,ind_pav_green) > 0.0_wp ) )  THEN
+       IF ( ( humidity ) .AND. ( frac_green > 0.0_wp ) )  THEN
           coef_green_2 = 4.0_wp * surf%emissivity(m,ind_pav_green) * sigma_sb *                   &
                          t_surf_green%val(m)**3 + f_qsws * dq_s_dt + lambda_surface_green         &
@@ -7007,5 +7013,5 @@
                             * t_surf_wall%val(m) )                                                &
                             / ( surf%c_surface(m) + coef_2 * dt_3d * tsc(2) )
-       IF ( ( .NOT. during_spinup ) .AND. (surf%frac(m,ind_wat_win) > 0.0_wp) )  THEN
+       IF ( ( .NOT. during_spinup ) .AND. (frac_win > 0.0_wp) )  THEN
           t_surf_window_p%val(m) = ( coef_window_1 * dt_3d * tsc(2) +                             &
                                    surf%c_surface_window(m) * t_surf_window%val(m) ) /            &
@@ -7027,7 +7033,7 @@
 !--    vpt_surface, which is, due to the lack of moisture on roofs, simply assumed to be the surface
 !--    temperature.
-       surf%pt_surface(m) = ( surf%frac(m,ind_veg_wall) * t_surf_wall_p%val(m)                    &
-                                 + surf%frac(m,ind_wat_win) * t_surf_window_p%val(m)              &
-                                 + surf%frac(m,ind_pav_green) * t_surf_green_p%val(m)             &
+       surf%pt_surface(m) = ( frac_wall * t_surf_wall_p%val(m)                                    &
+                                 + frac_win * t_surf_window_p%val(m)                              &
+                                 + frac_green * t_surf_green_p%val(m)                             &
                              ) / exner(k)
 
@@ -7069,11 +7075,11 @@
 !--    Calculate new fluxes
 !--    Rad_net_l is never used!
-       surf%rad_net_l(m) = surf%rad_net_l(m) + surf%frac(m,ind_veg_wall)                          &
+       surf%rad_net_l(m) = surf%rad_net_l(m) + frac_wall                                          &
                                  * sigma_sb * surf%emissivity(m,ind_veg_wall)                     &
                                  * ( t_surf_wall_p%val(m)**4 - t_surf_wall%val(m)**4 )            &
-                                 + surf%frac(m,ind_wat_win) * sigma_sb                            &
+                                 + frac_win * sigma_sb                                            &
                                  * surf%emissivity(m,ind_wat_win)                                 &
                                  * ( t_surf_window_p%val(m)**4 - t_surf_window%val(m)**4 )        &
-                                 + surf%frac(m,ind_pav_green) * sigma_sb                          &
+                                 + frac_green * sigma_sb                                          &
                                  * surf%emissivity(m,ind_pav_green)                               &
                                  * ( t_surf_green_p%val(m)**4 - t_surf_green%val(m)**4 )
@@ -7088,9 +7094,9 @@
 !--    Ground/wall/roof surface heat flux
        surf%wshf_eb(m)   = - f_shf  * ( surf%pt1(m) - t_surf_wall_p%val(m) / exner(k) )           &
-                                 * surf%frac(m,ind_veg_wall) - f_shf_window                       &
+                                 * frac_wall - f_shf_window                                       &
                                  * ( surf%pt1(m) - t_surf_window_p%val(m) / exner(k) )            &
-                                 * surf%frac(m,ind_wat_win) - f_shf_green                         &
+                                 * frac_win - f_shf_green                                         &
                                  * ( surf%pt1(m) - t_surf_green_p%val(m) / exner(k) )             &
-                                 * surf%frac(m,ind_pav_green)
+                                 * frac_green
 !
 !--    Store kinematic surface heat fluxes for utilization in other processes diffusion_s,
@@ -7103,5 +7109,5 @@
        ENDIF
 
-       IF ( humidity .AND.  surf%frac(m,ind_pav_green) > 0.0_wp )  THEN!
+       IF ( humidity .AND.  frac_green > 0.0_wp )  THEN!
 !--       Calculate true surface resistance
           IF ( upward ) THEN
@@ -7125,5 +7131,5 @@
                 IF ( m_liq_usm_h(l)%val(m) /= m_liq_max )  THEN
                    surf%qsws_liq(m) = surf%qsws_liq(m)                                            &
-                                            + surf%frac(m,ind_pav_green)                          &
+                                            + frac_green                                          &
                                             * prr(k+k_off,j+j_off,i+i_off) * hyrho(k+k_off)       &
                                             * 0.001_wp * rho_l * l_v
@@ -7191,12 +7197,4 @@
        ELSE
           surf%r_s(m) = 1.0E10_wp
-       ENDIF
-!
-!--    During spinup green and window fraction are set to zero. Here, the original values are
-!--    restored.
-       IF ( during_spinup )  THEN
-          surf%frac(m,ind_wat_win)   = frac_win
-          surf%frac(m,ind_veg_wall)  = frac_wall
-          surf%frac(m,ind_pav_green) = frac_green
        ENDIF
 
