Pull request #6: Dwd explicit loops

Merge in ESCAPE/dwarf-p-radiation-ecrad from dwd-explicit-loops to master

* commit 'c09f783c90e004e3004303f6b78c50e674a7b7db':
  [dwd-explicit-loops] Added comments when and by whom compiler directives were introduced
  [dwd-explicit-loops] Removed accidentally committed renaming of yomhook
  [dwd-explicit-loops] Update radiation/radiation_mcica_lw.F90 and radiation/radiation_mcica_sw.F90
  [dwd-explicit-loops] Update radiation/radiation_liquid_optics_socrates.F90
  [dwd-explicit-loops] Update radiation/radiation_ifs_rrtm.F90
  [dwd-explicit-loops] Update radiation/radiation_ice_optics_fu.F90
  [dwd-explicit-loops] Update radiation/radiation_cloud_optics.F90
  [dwd-explicit-loops] Update for radiation/radiation_aerosol_optics.F90
  [dwd-explicit-loops] Update radiation/radiation_adding_ica_lw.F90

radiation/radiation_adding_ica_lw.F90

@@ -295,7 +295,7 @@ subroutine calc_fluxes_no_scattering_lw(ncol, nlev, &
     real(jprb), intent(out), dimension(ncol, nlev+1) :: flux_up, flux_dn
 
     ! Loop index for model level
-    integer :: jlev
+    integer :: jlev, jcol
 
     real(jprb) :: hook_handle
 
@@ -306,17 +306,25 @@ subroutine calc_fluxes_no_scattering_lw(ncol, nlev, &
 
     ! Work down through the atmosphere computing the downward fluxes
     ! at each half-level
+! Added for DWD (2020)
+!NEC$ outerloop_unroll(8)
     do jlev = 1,nlev
-      flux_dn(:,jlev+1) = transmittance(:,jlev)*flux_dn(:,jlev) + source_dn(:,jlev)
+      do jcol = 1,ncol
+        flux_dn(jcol,jlev+1) = transmittance(jcol,jlev)*flux_dn(jcol,jlev) + source_dn(jcol,jlev)
+      end do
     end do
 
     ! Surface reflection and emission
     flux_up(:,nlev+1) = emission_surf + albedo_surf * flux_dn(:,nlev+1)
 
     ! Work back up through the atmosphere computing the upward fluxes
     ! at each half-level
+! Added for DWD (2020)
+!NEC$ outerloop_unroll(8)
     do jlev = nlev,1,-1
-      flux_up(:,jlev) = transmittance(:,jlev)*flux_up(:,jlev+1) + source_up(:,jlev)
+      do jcol = 1,ncol
+        flux_up(jcol,jlev) = transmittance(jcol,jlev)*flux_up(jcol,jlev+1) + source_up(jcol,jlev)
+      end do
     end do
 
     if (lhook) call dr_hook('radiation_adding_ica_lw:calc_fluxes_no_scattering_lw',1,hook_handle)

radiation/radiation_aerosol_optics.F90

@@ -385,14 +385,14 @@ subroutine add_aerosol_optics_direct(nlev,istartcol,iendcol, &
     ! asymmetry-times-scattering-optical-depth for all the aerosols at
     ! a point in space for each spectral band of the shortwave and
     ! longwave spectrum
-    real(jprb), dimension(config%n_bands_sw) &
+    real(jprb), dimension(config%n_bands_sw,nlev) &
          & :: od_sw_aerosol, scat_sw_aerosol, scat_g_sw_aerosol
-    real(jprb), dimension(config%n_bands_lw) :: od_lw_aerosol
-    real(jprb), dimension(config%n_bands_lw_if_scattering) &
+    real(jprb), dimension(config%n_bands_lw,nlev) :: od_lw_aerosol
+    real(jprb), dimension(config%n_bands_lw_if_scattering,nlev) &
          & :: scat_lw_aerosol, scat_g_lw_aerosol
 
     ! Loop indices for column, level, g point and band
-    integer :: jcol, jlev, jg
+    integer :: jcol, jlev, jg, jb
 
     ! Range of levels over which aerosols are present
     integer :: istartlev, iendlev
@@ -421,33 +421,38 @@ subroutine add_aerosol_optics_direct(nlev,istartcol,iendcol, &
 
       ! Loop over position
       do jcol = istartcol,iendcol
+! Added for DWD (2020)
+!NEC$ forced_collapse
         do jlev = istartlev,iendlev
-          od_sw_aerosol = aerosol%od_sw(:,jlev,jcol)
-          scat_sw_aerosol = aerosol%ssa_sw(:,jlev,jcol) * od_sw_aerosol
-          scat_g_sw_aerosol = aerosol%g_sw(:,jlev,jcol) * scat_sw_aerosol
-
-          if (.not. config%do_sw_delta_scaling_with_gases) then
-            ! Delta-Eddington scaling on aerosol only.  Note that if
-            ! do_sw_delta_scaling_with_gases==.true. then the delta
-            ! scaling is done to the cloud-aerosol-gas mixture inside
-            ! the solver
-            call delta_eddington_extensive(od_sw_aerosol, scat_sw_aerosol, &
-                 &                         scat_g_sw_aerosol)
-          end if
-
-          ! Combine aerosol shortwave scattering properties with gas
-          ! properties (noting that any gas scattering will have an
-          ! asymmetry factor of zero)
-          if (od_sw_aerosol(1) > 0.0_jprb) then
+          do jb = 1,config%n_bands_sw
+            od_sw_aerosol(jb,jlev) = aerosol%od_sw(jb,jlev,jcol)
+            scat_sw_aerosol(jb,jlev) = aerosol%ssa_sw(jb,jlev,jcol) * od_sw_aerosol(jb,jlev)
+            scat_g_sw_aerosol(jb,jlev) = aerosol%g_sw(jb,jlev,jcol) * scat_sw_aerosol(jb,jlev)
+
+            if (.not. config%do_sw_delta_scaling_with_gases) then
+              ! Delta-Eddington scaling on aerosol only.  Note that if
+              ! do_sw_delta_scaling_with_gases==.true. then the delta
+              ! scaling is done to the cloud-aerosol-gas mixture
+              ! inside the solver
+              call delta_eddington_extensive(od_sw_aerosol(jb,jlev), scat_sw_aerosol(jb,jlev), &
+                   &                         scat_g_sw_aerosol(jb,jlev))
+            end if
+          end do
+        end do
+        ! Combine aerosol shortwave scattering properties with gas
+        ! properties (noting that any gas scattering will have an
+        ! asymmetry factor of zero)
+        do jlev = istartlev,iendlev
+          if (od_sw_aerosol(1,jlev) > 0.0_jprb) then
             do jg = 1,config%n_g_sw
               iband = config%i_band_from_reordered_g_sw(jg)
-              local_od = od_sw(jg,jlev,jcol) + od_sw_aerosol(iband)
+              local_od = od_sw(jg,jlev,jcol) + od_sw_aerosol(iband,jlev)
               local_scat = ssa_sw(jg,jlev,jcol) * od_sw(jg,jlev,jcol) &
-                   &  + scat_sw_aerosol(iband)
+                   &  + scat_sw_aerosol(iband,jlev)
               ! Note that asymmetry_sw of gases is zero so the following
               ! simply weights the aerosol asymmetry by the scattering
               ! optical depth
-              g_sw(jg,jlev,jcol) = scat_g_sw_aerosol(iband) / local_scat
+              g_sw(jg,jlev,jcol) = scat_g_sw_aerosol(iband,jlev) / local_scat
               ssa_sw(jg,jlev,jcol) = local_scat / local_od
               od_sw (jg,jlev,jcol) = local_od
             end do
@@ -475,27 +480,32 @@ subroutine add_aerosol_optics_direct(nlev,istartcol,iendcol, &
 
         ! Loop over position
         do jcol = istartcol,iendcol
+! Added for DWD (2020)
+!NEC$ forced_collapse
           do jlev = istartlev,iendlev
-            od_lw_aerosol = aerosol%od_lw(:,jlev,jcol)
-            scat_lw_aerosol = aerosol%ssa_lw(:,jlev,jcol) * od_lw_aerosol
-            scat_g_lw_aerosol = aerosol%g_lw(:,jlev,jcol) * scat_lw_aerosol
-
-            call delta_eddington_extensive(od_lw_aerosol, scat_lw_aerosol, &
-                 &                         scat_g_lw_aerosol)
+            do jb = 1,config%n_bands_lw
+              od_lw_aerosol(jb,jlev) = aerosol%od_lw(jb,jlev,jcol)
+              scat_lw_aerosol(jb,jlev) = aerosol%ssa_lw(jb,jlev,jcol) * od_lw_aerosol(jb,jlev)
+              scat_g_lw_aerosol(jb,jlev) = aerosol%g_lw(jb,jlev,jcol) * scat_lw_aerosol(jb,jlev)
 
+              call delta_eddington_extensive(od_lw_aerosol(jb,jlev), scat_lw_aerosol(jb,jlev), &
+                   &                         scat_g_lw_aerosol(jb,jlev))
+            end do
+          end do
+          do jlev = istartlev,iendlev
             do jg = 1,config%n_g_lw
               iband = config%i_band_from_reordered_g_lw(jg)
-              if (od_lw_aerosol(iband) > 0.0_jprb) then
+              if (od_lw_aerosol(iband,jlev) > 0.0_jprb) then
                 ! All scattering is due to aerosols, therefore the
                 ! asymmetry factor is equal to the value for aerosols
-                if (scat_lw_aerosol(iband) > 0.0_jprb) then
-                  g_lw(jg,jlev,jcol) = scat_g_lw_aerosol(iband) &
-                       &  / scat_lw_aerosol(iband)
+                if (scat_lw_aerosol(iband,jlev) > 0.0_jprb) then
+                  g_lw(jg,jlev,jcol) = scat_g_lw_aerosol(iband,jlev) &
+                       &  / scat_lw_aerosol(iband,jlev)
                 else
                   g_lw(jg,jlev,jcol) = 0.0_jprb
                 end if
-                local_od = od_lw(jg,jlev,jcol) + od_lw_aerosol(iband)
-                ssa_lw(jg,jlev,jcol) = scat_lw_aerosol(iband) / local_od
+                local_od = od_lw(jg,jlev,jcol) + od_lw_aerosol(iband,jlev)
+                ssa_lw(jg,jlev,jcol) = scat_lw_aerosol(iband,jlev) / local_od
                 od_lw (jg,jlev,jcol) = local_od
               end if
             end do
@@ -506,15 +516,21 @@ subroutine add_aerosol_optics_direct(nlev,istartcol,iendcol, &
 
         ! Loop over position
         do jcol = istartcol,iendcol
+! Added for DWD (2020)
+!NEC$ forced_collapse
           do jlev = istartlev,iendlev
             ! If aerosol longwave scattering is not included then we
             ! weight the optical depth by the single scattering
             ! co-albedo
-            od_lw_aerosol = aerosol%od_lw(:,jlev,jcol) &
-                 &  * (1.0_jprb - aerosol%ssa_lw(:,jlev,jcol))
+            do jb = 1, config%n_bands_lw
+              od_lw_aerosol(jb,jlev) = aerosol%od_lw(jb,jlev,jcol) &
+                 &  * (1.0_jprb - aerosol%ssa_lw(jb,jlev,jcol))
+            end do
+          end do
+          do jlev = istartlev,iendlev
             do jg = 1,config%n_g_lw
               od_lw(jg,jlev,jcol) = od_lw(jg,jlev,jcol) &
-                   &  + od_lw_aerosol(config%i_band_from_reordered_g_lw(jg))
+                   &  + od_lw_aerosol(config%i_band_from_reordered_g_lw(jg),jlev)
             end do
           end do
         end do

radiation/radiation_cloud_optics.F90

@@ -270,7 +270,7 @@ subroutine cloud_optics(nlev,istartcol,iendcol, &
     ! access
     type(cloud_optics_type), pointer :: ho
 
-    integer    :: jcol, jlev
+    integer    :: jcol, jlev, jb
 
     real(jprb) :: hook_handle
 
@@ -457,29 +457,39 @@ subroutine cloud_optics(nlev,istartcol,iendcol, &
 
           ! Combine liquid and ice 
           if (config%do_lw_cloud_scattering) then
-            od_lw_cloud(:,jlev,jcol) = od_lw_liq + od_lw_ice
-            where (scat_od_lw_liq+scat_od_lw_ice > 0.0_jprb)
-              g_lw_cloud(:,jlev,jcol) = (g_lw_liq * scat_od_lw_liq &
-                   &  + g_lw_ice * scat_od_lw_ice) &
-                   &  / (scat_od_lw_liq+scat_od_lw_ice)
-            elsewhere
-              g_lw_cloud(:,jlev,jcol) = 0.0_jprb
-            end where
-            ssa_lw_cloud(:,jlev,jcol) = (scat_od_lw_liq + scat_od_lw_ice) &
-                 &                    / (od_lw_liq + od_lw_ice)
+! Added for DWD (2020)
+!NEC$ shortloop
+            do jb = 1, config%n_bands_lw
+              od_lw_cloud(jb,jlev,jcol) = od_lw_liq(jb) + od_lw_ice(jb)
+              if (scat_od_lw_liq(jb)+scat_od_lw_ice(jb) > 0.0_jprb) then
+                g_lw_cloud(jb,jlev,jcol) = (g_lw_liq(jb) * scat_od_lw_liq(jb) &
+                   &  + g_lw_ice(jb) * scat_od_lw_ice(jb)) &
+                   &  / (scat_od_lw_liq(jb)+scat_od_lw_ice(jb))
+              else
+                g_lw_cloud(jb,jlev,jcol) = 0.0_jprb
+              end if
+              ssa_lw_cloud(jb,jlev,jcol) = (scat_od_lw_liq(jb) + scat_od_lw_ice(jb)) &
+                 &                    / (od_lw_liq(jb) + od_lw_ice(jb))
+            enddo
           else
             ! If longwave scattering is to be neglected then the
             ! best approximation is to set the optical depth equal
             ! to the absorption optical depth
+! Added for DWD (2020)
+!NEC$ shortloop
             od_lw_cloud(:,jlev,jcol) = od_lw_liq - scat_od_lw_liq &
                  &                   + od_lw_ice - scat_od_lw_ice
           end if
-          od_sw_cloud(:,jlev,jcol) = od_sw_liq + od_sw_ice
-          g_sw_cloud(:,jlev,jcol) = (g_sw_liq * scat_od_sw_liq &
-               &  + g_sw_ice * scat_od_sw_ice) &
-               &  / (scat_od_sw_liq + scat_od_sw_ice)
-          ssa_sw_cloud(:,jlev,jcol) &
-               &  = (scat_od_sw_liq + scat_od_sw_ice) / (od_sw_liq + od_sw_ice)
+! Added for DWD (2020)
+!NEC$ shortloop
+          do jb = 1, config%n_bands_sw
+            od_sw_cloud(jb,jlev,jcol) = od_sw_liq(jb) + od_sw_ice(jb)
+            g_sw_cloud(jb,jlev,jcol) = (g_sw_liq(jb) * scat_od_sw_liq(jb) &
+               &  + g_sw_ice(jb) * scat_od_sw_ice(jb)) &
+               &  / (scat_od_sw_liq(jb) + scat_od_sw_ice(jb))
+            ssa_sw_cloud(jb,jlev,jcol) &
+               &  = (scat_od_sw_liq(jb) + scat_od_sw_ice(jb)) / (od_sw_liq(jb) + od_sw_ice(jb))
+          enddo
         end if ! Cloud present
       end do ! Loop over column
     end do ! Loop over level

radiation/radiation_ice_optics_fu.F90

@@ -60,6 +60,7 @@ subroutine calc_ice_optics_fu_sw(nb, coeff, ice_wp, &
     ! Ice water path in g m-2
     real (jprb) :: iwp_gm_2
 
+    integer :: jb
     !real(jprb)  :: hook_handle
 
     !if (lhook) call dr_hook('radiation_ice_optics:calc_ice_optics_fu_sw',0,hook_handle)
@@ -69,12 +70,16 @@ subroutine calc_ice_optics_fu_sw(nb, coeff, ice_wp, &
     inv_de_um = 1.0_jprb / de_um
     iwp_gm_2  = ice_wp * 1000.0_jprb
 
-    od = iwp_gm_2 * (coeff(1:nb,1) + coeff(1:nb,2) * inv_de_um)
-    scat_od = od * (1.0_jprb - (coeff(1:nb,3) + de_um*(coeff(1:nb,4) &
-         &  + de_um*(coeff(1:nb,5) + de_um*coeff(1:nb,6)))))
-    g = min(coeff(1:nb,7) + de_um*(coeff(1:nb,8) &
-         &  + de_um*(coeff(1:nb,9) + de_um*coeff(1:nb,10))), &
+! Added for DWD (2020)
+!NEC$ shortloop
+    do jb = 1, nb
+      od(jb) = iwp_gm_2 * (coeff(jb,1) + coeff(jb,2) * inv_de_um)
+      scat_od(jb) = od(jb) * (1.0_jprb - (coeff(jb,3) + de_um*(coeff(jb,4) &
+         &  + de_um*(coeff(jb,5) + de_um*coeff(jb,6)))))
+      g(jb) = min(coeff(jb,7) + de_um*(coeff(jb,8) &
+         &  + de_um*(coeff(jb,9) + de_um*coeff(jb,10))), &
          &  MaxAsymmetryFactor)
+    end do
 
     !if (lhook) call dr_hook('radiation_ice_optics:calc_ice_optics_fu_sw',1,hook_handle)
 
@@ -105,6 +110,7 @@ subroutine calc_ice_optics_fu_lw(nb, coeff, ice_wp, &
     ! Ice water path in g m-2
     real (jprb) :: iwp_gm_2
 
+    integer :: jb
     !real(jprb)  :: hook_handle
 
     !if (lhook) call dr_hook('radiation_ice_optics:calc_ice_optics_fu_lw',0,hook_handle)
@@ -115,13 +121,17 @@ subroutine calc_ice_optics_fu_lw(nb, coeff, ice_wp, &
     inv_de_um = 1.0_jprb / de_um
     iwp_gm_2  = ice_wp * 1000.0_jprb
 
-    od = iwp_gm_2 * (coeff(1:nb,1) + inv_de_um*(coeff(1:nb,2) &
-         &  + inv_de_um*coeff(1:nb,3)))
-    scat_od = od - iwp_gm_2*inv_de_um*(coeff(1:nb,4) + de_um*(coeff(1:nb,5) &
-         &  + de_um*(coeff(1:nb,6) + de_um*coeff(1:nb,7))))
-    g = min(coeff(1:nb,8) + de_um*(coeff(1:nb,9) &
-         &  + de_um*(coeff(1:nb,10) + de_um*coeff(1:nb,11))), &
+! Added for DWD (2020)
+!NEC$ shortloop
+    do jb = 1, nb
+      od(jb) = iwp_gm_2 * (coeff(jb,1) + inv_de_um*(coeff(jb,2) &
+         &  + inv_de_um*coeff(jb,3)))
+      scat_od(jb) = od(jb) - iwp_gm_2*inv_de_um*(coeff(jb,4) + de_um*(coeff(jb,5) &
+         &  + de_um*(coeff(jb,6) + de_um*coeff(jb,7))))
+      g(jb) = min(coeff(jb,8) + de_um*(coeff(jb,9) &
+         &  + de_um*(coeff(jb,10) + de_um*coeff(jb,11))), &
          &  MaxAsymmetryFactor)
+    end do
 
     !if (lhook) call dr_hook('radiation_ice_optics:calc_ice_optics_fu_lw',1,hook_handle)
 

radiation/radiation_ifs_rrtm.F90

@@ -606,7 +606,7 @@ subroutine planck_function_atmos(nlev,istartcol,iendcol, &
     ! Temperature (K) of a half-level
     real(jprb) :: temperature
 
-    real(jprb) :: factor
+    real(jprb) :: factor, planck_tmp(istartcol:iendcol,config%n_g_lw)
     real(jprb) :: ZFLUXFAC
 
     integer :: jlev, jgreorder, jg, ig, iband, jband, jcol, ilevoffset
@@ -691,7 +691,10 @@ subroutine planck_function_atmos(nlev,istartcol,iendcol, &
         else
           do jg = 1,config%n_g_lw
             iband = config%i_band_from_g_lw(jg)
-            planck_hl(jg,jlev,:) = planck_store(:,iband) * PFRAC(:,jg,nlev+2-jlev)
+            planck_tmp(:,jg) = planck_store(:,iband) * PFRAC(:,jg,nlev+2-jlev)
+          enddo
+          do jcol = istartcol,iendcol
+            planck_hl(:,jlev,jcol) = planck_tmp(jcol,:)
           end do
         end if
       end if

radiation/radiation_liquid_optics_socrates.F90

@@ -51,6 +51,7 @@ subroutine calc_liq_optics_socrates(nb, coeff, lwp, re_in, od, scat_od, g)
     ! Total optical depth, scattering optical depth and asymmetry factor
     real(jprb), intent(out) :: od(nb), scat_od(nb), g(nb)
 
+    integer    :: jb
     ! Local effective radius (m), after applying bounds
     real(jprb) :: re
 
@@ -61,14 +62,18 @@ subroutine calc_liq_optics_socrates(nb, coeff, lwp, re_in, od, scat_od, g)
     ! Apply the bounds of validity to effective radius
     re = max(MinEffectiveRadius, min(re_in, MaxEffectiveRadius))
 
-    od = lwp * (coeff(1:nb,1) + re*(coeff(1:nb,2) + re*coeff(1:nb,3))) &
-         &  / (1.0_jprb + re*(coeff(1:nb,4) + re*(coeff(1:nb,5) &
-         &  + re*coeff(1:nb,6))))
-    scat_od = od * (1.0_jprb &
-         &  - (coeff(1:nb,7) + re*(coeff(1:nb,8) + re*coeff(1:nb,9))) &
-         &  / (1.0_jprb + re*(coeff(1:nb,10) + re*coeff(1:nb,11))))
-    g = (coeff(1:nb,12) + re*(coeff(1:nb,13) + re*coeff(1:nb,14))) &
-         &  / (1.0_jprb + re*(coeff(1:nb,15) + re*coeff(1:nb,16)))
+! Added for DWD (2020)
+!NEC$ shortloop
+    do jb = 1, nb
+      od(jb) = lwp * (coeff(jb,1) + re*(coeff(jb,2) + re*coeff(jb,3))) &
+         &  / (1.0_jprb + re*(coeff(jb,4) + re*(coeff(jb,5) &
+         &  + re*coeff(jb,6))))
+      scat_od(jb) = od(jb) * (1.0_jprb &
+         &  - (coeff(jb,7) + re*(coeff(jb,8) + re*coeff(jb,9))) &
+         &  / (1.0_jprb + re*(coeff(jb,10) + re*coeff(jb,11))))
+      g(jb) = (coeff(jb,12) + re*(coeff(jb,13) + re*coeff(jb,14))) &
+         &  / (1.0_jprb + re*(coeff(jb,15) + re*coeff(jb,16)))
+    enddo
 
     !if (lhook) call dr_hook('radiation_liquid_optics_socrates:calc_liq_optics_socrates',1,hook_handle)