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