Index: LMDZ6/branches/LMDZ_ECRad/libf/phylmd/ecrad/radiation_flux.F90
===================================================================
--- LMDZ6/branches/LMDZ_ECRad/libf/phylmd/ecrad/radiation_flux.F90	(revision 4175)
+++ LMDZ6/branches/LMDZ_ECRad/libf/phylmd/ecrad/radiation_flux.F90	(revision 4444)
@@ -100,4 +100,11 @@
   end type flux_type
 
+! Added for DWD (2020)
+#ifdef __SX__
+      logical, parameter :: use_indexed_sum_vec = .true.
+#else
+      logical, parameter :: use_indexed_sum_vec = .false.
+#endif
+
 contains
 
@@ -132,5 +139,5 @@
         if (config%n_spec_lw == 0) then
           write(nulerr,'(a)') '*** Error: number of LW spectral points to save not yet defined ' &
-               & // 'so cannot allocated spectral flux arrays'
+               & // 'so cannot allocate spectral flux arrays'
           call radiation_abort()
         end if
@@ -321,4 +328,11 @@
     end if
 
+    if (allocated(this%lw_dn_surf_g))               deallocate(this%lw_dn_surf_g)
+    if (allocated(this%lw_dn_surf_clear_g))         deallocate(this%lw_dn_surf_clear_g)
+    if (allocated(this%sw_dn_diffuse_surf_g))       deallocate(this%sw_dn_diffuse_surf_g)
+    if (allocated(this%sw_dn_direct_surf_g))        deallocate(this%sw_dn_direct_surf_g)
+    if (allocated(this%sw_dn_diffuse_surf_clear_g)) deallocate(this%sw_dn_diffuse_surf_clear_g)
+    if (allocated(this%sw_dn_direct_surf_clear_g))  deallocate(this%sw_dn_direct_surf_clear_g)
+
     if (lhook) call dr_hook('radiation_flux:deallocate',1,hook_handle)
 
@@ -349,28 +363,56 @@
     if (config%do_sw .and. config%do_surface_sw_spectral_flux) then
 
-      do jcol = istartcol,iendcol
-        call indexed_sum(this%sw_dn_direct_surf_g(:,jcol), &
-             &           config%i_band_from_reordered_g_sw, &
-             &           this%sw_dn_direct_surf_band(:,jcol))
-        call indexed_sum(this%sw_dn_diffuse_surf_g(:,jcol), &
-             &           config%i_band_from_reordered_g_sw, &
-             &           this%sw_dn_surf_band(:,jcol))
-        this%sw_dn_surf_band(:,jcol) &
-             &  = this%sw_dn_surf_band(:,jcol) &
-             &  + this%sw_dn_direct_surf_band(:,jcol)
-      end do
+      if (use_indexed_sum_vec) then
+        call indexed_sum_vec(this%sw_dn_direct_surf_g, &
+             &               config%i_band_from_reordered_g_sw, &
+             &               this%sw_dn_direct_surf_band, istartcol, iendcol)
+        call indexed_sum_vec(this%sw_dn_diffuse_surf_g, &
+             &               config%i_band_from_reordered_g_sw, &
+             &               this%sw_dn_surf_band, istartcol, iendcol)
+        do jcol = istartcol,iendcol
+          this%sw_dn_surf_band(:,jcol) &
+               &  = this%sw_dn_surf_band(:,jcol) &
+               &  + this%sw_dn_direct_surf_band(:,jcol)
+        end do
+      else
+        do jcol = istartcol,iendcol
+          call indexed_sum(this%sw_dn_direct_surf_g(:,jcol), &
+               &           config%i_band_from_reordered_g_sw, &
+               &           this%sw_dn_direct_surf_band(:,jcol))
+          call indexed_sum(this%sw_dn_diffuse_surf_g(:,jcol), &
+               &           config%i_band_from_reordered_g_sw, &
+               &           this%sw_dn_surf_band(:,jcol))
+          this%sw_dn_surf_band(:,jcol) &
+               &  = this%sw_dn_surf_band(:,jcol) &
+               &  + this%sw_dn_direct_surf_band(:,jcol)
+        end do
+      end if
 
       if (config%do_clear) then
-        do jcol = istartcol,iendcol
-          call indexed_sum(this%sw_dn_direct_surf_clear_g(:,jcol), &
-               &           config%i_band_from_reordered_g_sw, &
-               &           this%sw_dn_direct_surf_clear_band(:,jcol))
-          call indexed_sum(this%sw_dn_diffuse_surf_clear_g(:,jcol), &
-               &           config%i_band_from_reordered_g_sw, &
-               &           this%sw_dn_surf_clear_band(:,jcol))
-          this%sw_dn_surf_clear_band(:,jcol) &
-               &  = this%sw_dn_surf_clear_band(:,jcol) &
-               &  + this%sw_dn_direct_surf_clear_band(:,jcol)
-        end do
+        if (use_indexed_sum_vec) then
+          call indexed_sum_vec(this%sw_dn_direct_surf_clear_g, &
+               &               config%i_band_from_reordered_g_sw, &
+               &               this%sw_dn_direct_surf_clear_band, istartcol, iendcol)
+          call indexed_sum_vec(this%sw_dn_diffuse_surf_clear_g, &
+               &               config%i_band_from_reordered_g_sw, &
+               &               this%sw_dn_surf_clear_band, istartcol, iendcol)
+          do jcol = istartcol,iendcol
+            this%sw_dn_surf_clear_band(:,jcol) &
+                 &  = this%sw_dn_surf_clear_band(:,jcol) &
+                 &  + this%sw_dn_direct_surf_clear_band(:,jcol)
+          end do
+        else
+          do jcol = istartcol,iendcol
+            call indexed_sum(this%sw_dn_direct_surf_clear_g(:,jcol), &
+                 &           config%i_band_from_reordered_g_sw, &
+                 &           this%sw_dn_direct_surf_clear_band(:,jcol))
+            call indexed_sum(this%sw_dn_diffuse_surf_clear_g(:,jcol), &
+                 &           config%i_band_from_reordered_g_sw, &
+                 &           this%sw_dn_surf_clear_band(:,jcol))
+            this%sw_dn_surf_clear_band(:,jcol) &
+                 &  = this%sw_dn_surf_clear_band(:,jcol) &
+                 &  + this%sw_dn_direct_surf_clear_band(:,jcol)
+          end do
+        end if
       end if
 
@@ -383,12 +425,21 @@
         this%sw_dn_direct_surf_canopy (:,istartcol:iendcol) = this%sw_dn_direct_surf_g (:,istartcol:iendcol)
       else if (config%do_nearest_spectral_sw_albedo) then
-        do jcol = istartcol,iendcol
-          call indexed_sum(this%sw_dn_direct_surf_g(:,jcol), &
-               &           config%i_albedo_from_band_sw(config%i_band_from_reordered_g_sw), &
-               &           this%sw_dn_direct_surf_canopy(:,jcol))
-          call indexed_sum(this%sw_dn_diffuse_surf_g(:,jcol), &
-               &           config%i_albedo_from_band_sw(config%i_band_from_reordered_g_sw), &
-               &           this%sw_dn_diffuse_surf_canopy(:,jcol))
-        end do
+        if (use_indexed_sum_vec) then
+          call indexed_sum_vec(this%sw_dn_direct_surf_g, &
+               &               config%i_albedo_from_band_sw(config%i_band_from_reordered_g_sw), &
+               &               this%sw_dn_direct_surf_canopy, istartcol, iendcol)
+          call indexed_sum_vec(this%sw_dn_diffuse_surf_g, &
+               &               config%i_albedo_from_band_sw(config%i_band_from_reordered_g_sw), &
+               &               this%sw_dn_diffuse_surf_canopy, istartcol, iendcol)
+        else
+          do jcol = istartcol,iendcol
+            call indexed_sum(this%sw_dn_direct_surf_g(:,jcol), &
+                 &           config%i_albedo_from_band_sw(config%i_band_from_reordered_g_sw), &
+                 &           this%sw_dn_direct_surf_canopy(:,jcol))
+            call indexed_sum(this%sw_dn_diffuse_surf_g(:,jcol), &
+                 &           config%i_albedo_from_band_sw(config%i_band_from_reordered_g_sw), &
+                 &           this%sw_dn_diffuse_surf_canopy(:,jcol))
+          end do
+        end if
       else
         ! More accurate calculations using weights, but requires
@@ -425,17 +476,29 @@
         this%lw_dn_surf_canopy(:,istartcol:iendcol) = this%lw_dn_surf_g(:,istartcol:iendcol)
       else if (config%do_nearest_spectral_lw_emiss) then
-        do jcol = istartcol,iendcol
-          call indexed_sum(this%lw_dn_surf_g(:,jcol), &
-               &           config%i_emiss_from_band_lw(config%i_band_from_reordered_g_lw), &
-               &           this%lw_dn_surf_canopy(:,jcol))
-        end do
+        if (use_indexed_sum_vec) then
+          call indexed_sum_vec(this%lw_dn_surf_g, &
+               &               config%i_emiss_from_band_lw(config%i_band_from_reordered_g_lw), &
+               &               this%lw_dn_surf_canopy, istartcol, iendcol)
+        else
+          do jcol = istartcol,iendcol
+            call indexed_sum(this%lw_dn_surf_g(:,jcol), &
+                 &           config%i_emiss_from_band_lw(config%i_band_from_reordered_g_lw), &
+                 &           this%lw_dn_surf_canopy(:,jcol))
+          end do
+        end if
       else
         ! Compute fluxes in each longwave emissivity interval using
         ! weights; first sum over g points to get the values in bands
-        do jcol = istartcol,iendcol
-          call indexed_sum(this%lw_dn_surf_g(:,jcol), &
-               &           config%i_band_from_reordered_g_lw, &
-               &           lw_dn_surf_band(:,jcol))
-        end do
+        if (use_indexed_sum_vec) then
+          call indexed_sum_vec(this%lw_dn_surf_g, &
+               &               config%i_band_from_reordered_g_lw, &
+               &               lw_dn_surf_band, istartcol, iendcol)
+        else
+          do jcol = istartcol,iendcol
+            call indexed_sum(this%lw_dn_surf_g(:,jcol), &
+                 &           config%i_band_from_reordered_g_lw, &
+                 &           lw_dn_surf_band(:,jcol))
+          end do
+        end if
         nalbedoband = size(config%lw_emiss_weights,1)
         this%lw_dn_surf_canopy(:,istartcol:iendcol) = 0.0_jprb
@@ -465,5 +528,5 @@
 
     use yomhook,          only : lhook, dr_hook
-    use radiation_config, only : out_of_bounds_2d
+    use radiation_check,  only : out_of_bounds_2d
 
     class(flux_type), intent(inout) :: this
@@ -497,5 +560,5 @@
   function heating_rate_out_of_physical_bounds(this, nlev, istartcol, iendcol, pressure_hl) result(is_bad)
     
-    use radiation_config, only : out_of_bounds_2d
+    use radiation_check, only : out_of_bounds_2d
     use radiation_constants, only : AccelDueToGravity
 
@@ -581,4 +644,25 @@
   end subroutine indexed_sum
 
+  !---------------------------------------------------------------------
+  ! Vectorized version of "add_indexed_sum"
+  subroutine indexed_sum_vec(source, ind, dest, ist, iend)
+
+    real(jprb), intent(in)  :: source(:,:)
+    integer,    intent(in)  :: ind(:)
+    real(jprb), intent(out) :: dest(:,:)
+    integer,    intent(in)  :: ist, iend
+
+    integer :: ig, jg, jc
+
+    dest = 0.0
+
+    do jg = lbound(source,1), ubound(source,1)
+      ig = ind(jg)
+      do jc = ist, iend
+        dest(ig,jc) = dest(ig,jc) + source(jg,jc)
+      end do
+    end do
+
+  end subroutine indexed_sum_vec
 
   !---------------------------------------------------------------------
