Changeset 4713 for palm/trunk/SOURCE/radiation_model_mod.f90

Timestamp:

Sep 29, 2020 12:02:05 PM (4 years ago)

Author:

pavelkrc

Message:

Correct and improve OpenMP parallelization, fix numerical instabilities

File:

• palm/trunk/SOURCE/radiation_model_mod.f90 (modified) (30 diffs)

palm/trunk/SOURCE/radiation_model_mod.f90

@@ -28 +28 @@
 ! -----------------
 ! $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 +5970 @@
      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 +6075 @@
 !--  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 +6164 @@
 
      IF ( surface_reflections)  THEN
+        !$OMP DO PRIVATE (i, j, k, isvf, isurf, isurfsrc) SCHEDULE (STATIC)
         DO  isvf = 1, nsvfl
            isurf = svfsurf(1, isvf)
@@ -6173 +6179 @@
            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 +6201 @@
         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 +6211 @@
         mrtinlw(imrt) = mrtsky(imrt) * rad_lw_in_diff(j,i)
      ENDDO
+     !$OMP END DO
 !
 !--  Direct radiation
@@ -6218 +6229 @@
         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 +6239 @@
            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 +6249 @@
                                      * 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 +6259 @@
         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 +6269 @@
          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 +6328 @@
              ENDIF
          ENDDO
+         !$OMP END DO
 
          pcbinsw(:) = pcbinswdir(:) + pcbinswdif(:)
@@ -6362 +6382 @@
      ENDIF
 
-!-- Next passes of radiation interactions:
+!--  Next passes of radiation interactions:
 !--  radiation reflections
 
@@ -6404 +6424 @@
 !
 !--      Reflected radiation
+         !$OMP DO PRIVATE (isvf, isurf, isurfsrc) SCHEDULE (STATIC)
          DO isvf = 1, nsvfl
              isurf = svfsurf(1, isvf)
@@ -6414 +6435 @@
              ENDIF
          ENDDO
+         !$OMP END DO
 !
 !--      NOTE: PC absorbtion and MRT from reflected can both be done at once
@@ -6421 +6443 @@
 !
 !--      Radiation absorbed by plant canopy
+         !$OMP DO PRIVATE (icsf, ipcgb, isurfsrc, asrc) SCHEDULE (STATIC)
          DO  icsf = 1, ncsfl
              ipcgb = csfsurf(1, icsf)
@@ -6435 +6458 @@
              ENDIF
          ENDDO
+         !$OMP END DO
 !
 !--      MRT reflected
+         !$OMP DO PRIVATE (imrtf, imrt, isurfsrc) SCHEDULE (STATIC)
          DO  imrtf = 1, nmrtf
             imrt = mrtfsurf(1, imrtf)
@@ -6443 +6468 @@
             mrtinlw(imrt) = mrtinlw(imrt) + mrtf(imrtf) * surfoutl(isurfsrc)
          ENDDO
+         !$OMP END DO
 
          IF ( trace_fluxes_above >= 0.0_wp )  THEN
@@ -6464 +6490 @@
      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 +6503 @@
              pabsswl = pabsswl + pcbinsw(ipcgb)
          ENDDO
+         !$OMP END DO
 
          IF ( humidity .AND. plant_canopy_transpiration ) THEN
@@ -6481 +6509 @@
              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 +6519 @@
                                                    pcm_transpiration_rate(kk,j,i),  &
                                                    pcm_latent_rate(kk,j,i) )
-              ENDDO
+             ENDDO
+             !$OMP END DO
          ENDIF
      ENDIF
@@ -6506 +6536 @@
 !    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 +6630 @@
      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 +6646 @@
                            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 +6656 @@
                                   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 +6664 @@
                                   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 +6676 @@
                                      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 +6686 @@
 
         ENDDO
+        !$OMP END DO
      ENDDO
 !