Index: palm/trunk/SOURCE/radiation_model_mod.f90
===================================================================
--- palm/trunk/SOURCE/radiation_model_mod.f90	(revision 4708)
+++ 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
 !
