Index: /palm/trunk/SOURCE/land_surface_model_mod.f90
===================================================================
--- /palm/trunk/SOURCE/land_surface_model_mod.f90	(revision 4712)
+++ /palm/trunk/SOURCE/land_surface_model_mod.f90	(revision 4713)
@@ -25,4 +25,11 @@
 ! -----------------
 ! $Id$
+! - Optimize calculation in case of calc_soil_moisture=false
+! - Improve OMP parallelization
+! - Do not calculate surface prognostic temperature for water surfaces (avoid numerical
+!   instabilities)
+! Author: J. Resler
+!
+! 4694 2020-09-23 15:09:19Z pavelkrc
 ! Fix reading of surface data from MPI restart file
 !
@@ -1734,5 +1741,5 @@
     i_off = surf%ioff
 
-    !$OMP PARALLEL PRIVATE (m, i, j, k, lambda_h_sat, ke, lambda_soil, lambda_surface,             &
+    !$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,   &
@@ -2038,8 +2045,13 @@
 !--    Implicit solution when the surface layer has no heat capacity,
 !--    otherwise use RK3 scheme.
-       surf_t_surface_p%var_1d(m) = ( coef_1 * dt_3d * tsc(2) + c_surface_tmp *&
-                          surf_t_surface%var_1d(m) ) / ( c_surface_tmp + coef_2&
+       IF ( surf%water_surface(m) )  THEN
+!
+!--       water temperature is fixed (avoid numerical inaccuracies)
+          surf_t_surface_p%var_1d(m) = surf_t_surface%var_1d(m)
+       ELSE
+          surf_t_surface_p%var_1d(m) = ( coef_1 * dt_3d * tsc(2) + c_surface_tmp * &
+                          surf_t_surface%var_1d(m) ) / ( c_surface_tmp + coef_2    &
                                              * dt_3d * tsc(2) )
-
+       ENDIF
 !
 !--    Add RK3 term
@@ -5545,5 +5557,6 @@
        ENDIF
 
-       !$OMP PARALLEL PRIVATE (m, k, lambda_temp, lambda_h_sat, ke, tend, gamma_temp, h_vg, m_total)
+       !$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)
        DO  m = 1, surf%ns
@@ -5636,54 +5649,51 @@
 
 
-             DO  k = nzb_soil, nzt_soil
-
-!
-!--             In order to prevent water tranport through paved surfaces,
-!--             conductivity and diffusivity are set to zero
-                IF ( surf%pavement_surface(m)  .AND.                           &
-                     k <= surf%nzt_pavement(m) )  THEN
-                   lambda_temp(k) = 0.0_wp
-                   gamma_temp(k)  = 0.0_wp
-
-                ELSE
-
-!
-!--                Calculate soil diffusivity at the center of the soil layers
-                   lambda_temp(k) = (- b_ch * surf%gamma_w_sat(k,m) * psi_sat  &
-                                    / surf%m_sat(k,m) ) * (                    &
-                                    MAX( surf_m_soil%var_2d(k,m),              &
-                                    surf%m_wilt(k,m) ) / surf%m_sat(k,m) )**(  &
-                                    b_ch + 2.0_wp )
-
-!
-!--                Parametrization of Van Genuchten
-!--                Calculate the hydraulic conductivity after Van Genuchten (1980)
-                   h_vg = ( ( ( surf%m_res(k,m) - surf%m_sat(k,m) ) /          &
-                              ( surf%m_res(k,m) -                              &
-                                MAX( surf_m_soil%var_2d(k,m), surf%m_wilt(k,m) )&
-                              )                                                &
-                            )**(                                               &
-                          surf%n_vg(k,m) / ( surf%n_vg(k,m) - 1.0_wp )         &
-                               ) - 1.0_wp                                      &
-                          )**( 1.0_wp / surf%n_vg(k,m) ) / surf%alpha_vg(k,m)
-
-                   gamma_temp(k) = surf%gamma_w_sat(k,m) * ( ( ( 1.0_wp +      &
-                          ( surf%alpha_vg(k,m) * h_vg )**surf%n_vg(k,m)        &
-                                                                  )**(         &
-                              1.0_wp - 1.0_wp / surf%n_vg(k,m)) - (            &
-                          surf%alpha_vg(k,m) * h_vg )**( surf%n_vg(k,m)        &
-                              - 1.0_wp) )**2 )                                 &
-                              / ( ( 1.0_wp + ( surf%alpha_vg(k,m) * h_vg       &
-                              )**surf%n_vg(k,m) )**( ( 1.0_wp  - 1.0_wp        &
-                              / surf%n_vg(k,m) ) *                             &
-                              ( surf%l_vg(k,m) + 2.0_wp) ) )
-
-                ENDIF
-
-             ENDDO
-
-
              IF ( calc_soil_moisture )  THEN
 
+                DO  k = nzb_soil, nzt_soil
+   !
+   !--             In order to prevent water tranport through paved surfaces,
+   !--             conductivity and diffusivity are set to zero
+                   IF ( surf%pavement_surface(m)  .AND.                           &
+                        k <= surf%nzt_pavement(m) )  THEN
+                      lambda_temp(k) = 0.0_wp
+                      gamma_temp(k)  = 0.0_wp
+
+                   ELSE
+
+   !
+   !--                Calculate soil diffusivity at the center of the soil layers
+                      lambda_temp(k) = (- b_ch * surf%gamma_w_sat(k,m) * psi_sat  &
+                                       / surf%m_sat(k,m) ) * (                    &
+                                       MAX( surf_m_soil%var_2d(k,m),              &
+                                       surf%m_wilt(k,m) ) / surf%m_sat(k,m) )**(  &
+                                       b_ch + 2.0_wp )
+
+   !
+   !--                Parametrization of Van Genuchten
+   !--                Calculate the hydraulic conductivity after Van Genuchten (1980)
+                      h_vg = ( ( ( surf%m_res(k,m) - surf%m_sat(k,m) ) /          &
+                                 ( surf%m_res(k,m) -                              &
+                                   MAX( surf_m_soil%var_2d(k,m), surf%m_wilt(k,m) )&
+                                 )                                                &
+                               )**(                                               &
+                             surf%n_vg(k,m) / ( surf%n_vg(k,m) - 1.0_wp )         &
+                                  ) - 1.0_wp                                      &
+                             )**( 1.0_wp / surf%n_vg(k,m) ) / surf%alpha_vg(k,m)
+
+                      gamma_temp(k) = surf%gamma_w_sat(k,m) * ( ( ( 1.0_wp +      &
+                             ( surf%alpha_vg(k,m) * h_vg )**surf%n_vg(k,m)        &
+                                                                     )**(         &
+                                 1.0_wp - 1.0_wp / surf%n_vg(k,m)) - (            &
+                             surf%alpha_vg(k,m) * h_vg )**( surf%n_vg(k,m)        &
+                                 - 1.0_wp) )**2 )                                 &
+                                 / ( ( 1.0_wp + ( surf%alpha_vg(k,m) * h_vg       &
+                                 )**surf%n_vg(k,m) )**( ( 1.0_wp  - 1.0_wp        &
+                                 / surf%n_vg(k,m) ) *                             &
+                                 ( surf%l_vg(k,m) + 2.0_wp) ) )
+
+                   ENDIF
+
+                ENDDO
 !
 !--             Prognostic equation for soil moisture content. Only performed,
@@ -7157,5 +7167,4 @@
                 ALLOCATE( t_soil_h(1)%var_2d(nzb_soil:nzt_soil+1,              &
                                           1:surf_lsm_h(1)%ns) )
-             write(9,*) surf_lsm_h(1)%ns
              READ ( 13 )  tmp_walltype_h_2d(1)%var_2d
           ENDIF
@@ -7753,11 +7762,11 @@
 !       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)
        DO  m = 1, surf_lsm_h(0)%ns
 !--       only upward facin horizontal surfaces are considered for water surface processing
-          i   = surf_lsm_h(0)%i(m)
-          j   = surf_lsm_h(0)%j(m)
 
           IF ( surf_lsm_h(0)%water_surface(m) )  THEN
-
+             i   = surf_lsm_h(0)%i(m)
+             j   = surf_lsm_h(0)%j(m)
 !
 !--          Disabled: FLake parameterization. Ideally, the Charnock
@@ -7813,4 +7822,5 @@
           ENDIF
        ENDDO
+       !$OMP END DO
 #if defined( __parallel )
        CALL MPI_ALLREDUCE( MPI_IN_PLACE, flag_exceed_z0, 1, MPI_LOGICAL,       &
Index: /palm/trunk/SOURCE/radiation_model_mod.f90
===================================================================
--- /palm/trunk/SOURCE/radiation_model_mod.f90	(revision 4712)
+++ /palm/trunk/SOURCE/radiation_model_mod.f90	(revision 4713)
@@ -28,4 +28,7 @@
 ! -----------------
 ! $Id$
+! Correct OpenMP parallelization including cycles with cumulative variables (J. Resler)
+!
+! 4708 2020-09-28 17:42:58Z suehring
 ! - Bugfix, correct mapping of RRTMG heating rates, as well as incoming/outgoing
 !   radiation onto the topography-following grid
@@ -5967,5 +5970,5 @@
      REAL(wp)                          ::  pabs_pc_lwdifl     !< total absorbed LW radiation in plant canopy from sky in local processor (W)
      REAL(wp)                          ::  pabs_pc_lwdif      !< total absorbed LW radiation in plant canopy from sky in all processors (W)
-
+!-   rotation related variables
      REAL(wp)                          ::  sun_direct_factor  !< factor for direct normal radiation from direct horizontal
      REAL(wp)                          ::  sin_rot            !< sine of rotation_angle
@@ -6072,4 +6075,6 @@
 !--  Set up thermal radiation from surfaces
      mm = 1
+!--  following code depends on the order of the execution
+!--  do not parallelize by OpenMP
      DO  i = nxl, nxr
         DO  j = nys, nyn
@@ -6159,4 +6164,5 @@
 
      IF ( surface_reflections)  THEN
+        !$OMP DO PRIVATE (i, j, k, isvf, isurf, isurfsrc) SCHEDULE (STATIC)
         DO  isvf = 1, nsvfl
            isurf = svfsurf(1, isvf)
@@ -6173,7 +6179,9 @@
            ENDIF
         ENDDO
+        !$OMP END DO
      ENDIF
 !
 !--  diffuse radiation using sky view factor
+     !$OMP DO PRIVATE (i, j, d, isurf) REDUCTION(+:pinswl, pinlwl) SCHEDULE (STATIC)
      DO isurf = 1, nsurfl
         j = surfl(iy, isurf)
@@ -6193,6 +6201,8 @@
         ENDIF
      ENDDO
+     !$OMP END DO
 !
 !--  MRT diffuse irradiance
+     !$OMP DO PRIVATE (i, j, imrt) SCHEDULE (STATIC)
      DO  imrt = 1, nmrtbl
         j = mrtbl(iy, imrt)
@@ -6201,4 +6211,5 @@
         mrtinlw(imrt) = mrtsky(imrt) * rad_lw_in_diff(j,i)
      ENDDO
+     !$OMP END DO
 !
 !--  Direct radiation
@@ -6218,4 +6229,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)
         DO isurf = 1, nsurfl
            j = surfl(iy, isurf)
@@ -6227,6 +6239,8 @@
            pinswl = pinswl + surfinswdir(isurf) * facearea(d)
         ENDDO
+        !$OMP END DO
 !
 !--     MRT direct irradiance
+        !$OMP DO PRIVATE (i, j, imrt) SCHEDULE (STATIC)
         DO  imrt = 1, nmrtbl
            j = mrtbl(iy, imrt)
@@ -6235,7 +6249,9 @@
                                      * 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)
      DO  imrtf = 1, nmrtf
         imrt = mrtfsurf(1, imrtf)
@@ -6243,4 +6259,5 @@
         mrtinlw(imrt) = mrtinlw(imrt) + mrtf(imrtf) * surfoutl(isurfsrc)
      ENDDO
+     !$OMP END DO
 !
 !--  Absorption in each local plant canopy grid box from the first atmospheric
@@ -6252,4 +6269,6 @@
          pcbinlw(:) = 0.0_wp
 
+         !$OMP 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
              ipcgb = csfsurf(1, icsf)
@@ -6309,4 +6328,5 @@
              ENDIF
          ENDDO
+         !$OMP END DO
 
          pcbinsw(:) = pcbinswdir(:) + pcbinswdif(:)
@@ -6362,5 +6382,5 @@
      ENDIF
 
-!-- Next passes of radiation interactions:
+!--  Next passes of radiation interactions:
 !--  radiation reflections
 
@@ -6404,4 +6424,5 @@
 !
 !--      Reflected radiation
+         !$OMP DO PRIVATE (isvf, isurf, isurfsrc) SCHEDULE (STATIC)
          DO isvf = 1, nsvfl
              isurf = svfsurf(1, isvf)
@@ -6414,4 +6435,5 @@
              ENDIF
          ENDDO
+         !$OMP END DO
 !
 !--      NOTE: PC absorbtion and MRT from reflected can both be done at once
@@ -6421,4 +6443,5 @@
 !
 !--      Radiation absorbed by plant canopy
+         !$OMP DO PRIVATE (icsf, ipcgb, isurfsrc, asrc) SCHEDULE (STATIC)
          DO  icsf = 1, ncsfl
              ipcgb = csfsurf(1, icsf)
@@ -6435,6 +6458,8 @@
              ENDIF
          ENDDO
+         !$OMP END DO
 !
 !--      MRT reflected
+         !$OMP DO PRIVATE (imrtf, imrt, isurfsrc) SCHEDULE (STATIC)
          DO  imrtf = 1, nmrtf
             imrt = mrtfsurf(1, imrtf)
@@ -6443,4 +6468,5 @@
             mrtinlw(imrt) = mrtinlw(imrt) + mrtf(imrtf) * surfoutl(isurfsrc)
          ENDDO
+         !$OMP END DO
 
          IF ( trace_fluxes_above >= 0.0_wp )  THEN
@@ -6464,4 +6490,5 @@
      IF ( npcbl > 0 )  THEN
          pcm_heating_rate(:,:,:) = 0.0_wp
+         !$OMP DO PRIVATE (ipcgb, i, j, k, kk) REDUCTION(+:pabsswl) SCHEDULE (STATIC)
          DO ipcgb = 1, npcbl
              j = pcbl(iy, ipcgb)
@@ -6476,4 +6503,5 @@
              pabsswl = pabsswl + pcbinsw(ipcgb)
          ENDDO
+         !$OMP END DO
 
          IF ( humidity .AND. plant_canopy_transpiration ) THEN
@@ -6481,4 +6509,5 @@
              pcm_transpiration_rate(:,:,:) = 0.0_wp
              pcm_latent_rate(:,:,:) = 0.0_wp
+             !$OMP DO PRIVATE (ipcgb, i, j, k, kk) SCHEDULE (STATIC)
              DO ipcgb = 1, npcbl
                  i = pcbl(ix, ipcgb)
@@ -6490,5 +6519,6 @@
                                                    pcm_transpiration_rate(kk,j,i),  &
                                                    pcm_latent_rate(kk,j,i) )
-              ENDDO
+             ENDDO
+             !$OMP END DO
          ENDIF
      ENDIF
@@ -6506,4 +6536,6 @@
 !    and claculate relevant radiation model-RTM coupling terms
      mm = 1
+!--  following code depends on the order of the execution
+!--  do not parallelize by OpenMP
      DO  i = nxl, nxr
         DO  j = nys, nyn
@@ -6598,4 +6630,6 @@
      ENDDO
 
+     !$OMP DO PRIVATE (i, d) REDUCTION(+:pabsswl, pabslwl, pemitlwl, pabs_surf_lwdifl) &
+     !$OMP&   SCHEDULE (STATIC)
      DO  i = 1, nsurfl
         d  = surfl(id, i)
@@ -6612,6 +6646,8 @@
                            emiss_surf(i) * facearea(d) * surfinlwdif(i)
      ENDDO
+     !$OMP END DO
 
      DO l = 0, 1
+        !$OMP 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)  +          &
@@ -6620,4 +6656,6 @@
                                   surf_usm_h(l)%rad_lw_out(m)
         ENDDO
+        !$OMP END DO
+        !$OMP 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)  +          &
@@ -6626,6 +6664,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
         DO  m = 1, surf_usm_v(l)%ns
@@ -6635,5 +6676,7 @@
                                      surf_usm_v(l)%rad_lw_out(m)
         ENDDO
+        !$OMP END DO
 !--     land
+        !$OMP 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)  +          &
@@ -6643,4 +6686,5 @@
 
         ENDDO
+        !$OMP END DO
      ENDDO
 !
Index: /palm/trunk/SOURCE/urban_surface_mod.f90
===================================================================
--- /palm/trunk/SOURCE/urban_surface_mod.f90	(revision 4712)
+++ /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
 
