Index: LMDZ6/branches/cirrus/libf/phylmd/ecrad/lmdz/calcul_cloud_overlap_decorr_len.F90
===================================================================
--- LMDZ6/branches/cirrus/libf/phylmd/ecrad/lmdz/calcul_cloud_overlap_decorr_len.F90	(revision 4945)
+++ LMDZ6/branches/cirrus/libf/phylmd/ecrad/lmdz/calcul_cloud_overlap_decorr_len.F90	(revision 5202)
@@ -146,7 +146,7 @@
 !  ENDIF
 ENDIF
-CALL writefield_phy('latitude',latitude_deg,1)
-CALL writefield_phy('pressure_hl',pressure_hl,klev+1)
-CALL writefield_phy('Ldecorel',PDECORR_LEN_EDGES_M,klev)
+!CALL writefield_phy('latitude',latitude_deg,1)
+!CALL writefield_phy('pressure_hl',pressure_hl,klev+1)
+!CALL writefield_phy('Ldecorel',PDECORR_LEN_EDGES_M,klev)
 ! -------------------------------------------------------------------
 
Index: LMDZ6/branches/cirrus/libf/phylmd/ecrad/lmdz/radiation_setup.F90
===================================================================
--- LMDZ6/branches/cirrus/libf/phylmd/ecrad/lmdz/radiation_setup.F90	(revision 4945)
+++ LMDZ6/branches/cirrus/libf/phylmd/ecrad/lmdz/radiation_setup.F90	(revision 5202)
@@ -141,5 +141,4 @@
            &  -9, &
            &   4 /)
-!   rad_config%aerosol_optics_override_file_name = 'aerosol_optics_lmdz.nc'
 
 
Index: LMDZ6/branches/cirrus/libf/phylmd/ecrad/lmdz/readaerosol_optic_ecrad.F90
===================================================================
--- LMDZ6/branches/cirrus/libf/phylmd/ecrad/lmdz/readaerosol_optic_ecrad.F90	(revision 4945)
+++ LMDZ6/branches/cirrus/libf/phylmd/ecrad/lmdz/readaerosol_optic_ecrad.F90	(revision 5202)
@@ -4,5 +4,5 @@
      flag_aerosol, flag_bc_internal_mixture, itap, rjourvrai, &
      pdtphys, pplay, paprs, t_seri, rhcl, presnivs, &
-     tr_seri, mass_solu_aero, mass_solu_aero_pi) 
+     tr_seri, mass_solu_aero, mass_solu_aero_pi, m_allaer) 
 !     tau_aero, piz_aero, cg_aero, &
 !     tausum_aero, drytausum_aero, tau3d_aero )
@@ -18,5 +18,5 @@
        concso4,concno3,concoa,concbc,concss,concdust,loadso4,loadoa,loadbc,loadss,loaddust, &
        loadno3,load_tmp1,load_tmp2,load_tmp3,load_tmp4,load_tmp5,load_tmp6,load_tmp7, & 
-       load_tmp8,load_tmp9,load_tmp10,m_allaer
+       load_tmp8,load_tmp9,load_tmp10
 
   USE infotrac_phy, ONLY: tracers, nqtot, nbtr
@@ -49,4 +49,8 @@
   REAL, DIMENSION(klon,klev), INTENT(OUT)     :: mass_solu_aero    ! Total mass for all soluble aerosols
   REAL, DIMENSION(klon,klev), INTENT(OUT)     :: mass_solu_aero_pi !     -"-     preindustrial values
+  REAL, DIMENSION(klon,klev,naero_tot), INTENT(OUT) :: m_allaer
+  ! AI a passer par la suite en argument si besoin pour ecrad
+  !REAL, DIMENSION(klon,klev,naero_tot), INTENT(OUT) :: m_allaer_pi !RAF
+
 !  REAL, DIMENSION(klon,klev,2,NSW), INTENT(OUT) :: tau_aero    ! Aerosol optical thickness
 !  REAL, DIMENSION(klon,klev,2,NSW), INTENT(OUT) :: piz_aero    ! Single scattering albedo aerosol
@@ -86,6 +90,5 @@
   REAL, DIMENSION(klon,klev)   :: nitrinscoarse_pi
   REAL, DIMENSION(klon,klev)   :: pdel, zrho
-!  REAL, DIMENSION(klon,klev,naero_tot) :: m_allaer
-  REAL, DIMENSION(klon,klev,naero_tot) :: m_allaer_pi !RAF  
+  REAL, DIMENSION(klon,klev,naero_tot) :: m_allaer_pi
 
   integer :: id_ASBCM, id_ASPOMM, id_ASSO4M, id_ASMSAM, id_CSSO4M, id_CSMSAM, id_SSSSM
Index: LMDZ6/branches/cirrus/libf/phylmd/ecrad/radiation/radiation_mcica_lw.F90
===================================================================
--- LMDZ6/branches/cirrus/libf/phylmd/ecrad/radiation/radiation_mcica_lw.F90	(revision 4945)
+++ LMDZ6/branches/cirrus/libf/phylmd/ecrad/radiation/radiation_mcica_lw.F90	(revision 5202)
@@ -18,6 +18,4 @@
 !   2017-07-12  R. Hogan  Call fast adding method if only clouds scatter
 !   2017-10-23  R. Hogan  Renamed single-character variables
-
-#include "ecrad_config.h"
 
 module radiation_mcica_lw
@@ -126,11 +124,4 @@
     ! Identify clear-sky layers
     logical :: is_clear_sky_layer(nlev)
-
-    ! Temporary storage for more efficient summation
-#ifdef DWD_REDUCTION_OPTIMIZATIONS
-    real(jprb), dimension(nlev+1,2) :: sum_aux
-#else
-    real(jprb) :: sum_up, sum_dn
-#endif
 
     ! Index of the highest cloudy layer
@@ -188,28 +179,6 @@
 
       ! Sum over g-points to compute broadband fluxes
-#ifdef DWD_REDUCTION_OPTIMIZATIONS
-      sum_aux(:,:) = 0.0_jprb
-      do jg = 1,ng
-        do jlev = 1,nlev+1
-          sum_aux(jlev,1) = sum_aux(jlev,1) + flux_up_clear(jg,jlev)
-          sum_aux(jlev,2) = sum_aux(jlev,2) + flux_dn_clear(jg,jlev)
-        end do
-      end do
-      flux%lw_up_clear(jcol,:) = sum_aux(:,1)
-      flux%lw_dn_clear(jcol,:) = sum_aux(:,2)
-#else
-      do jlev = 1,nlev+1
-        sum_up = 0.0_jprb
-        sum_dn = 0.0_jprb
-        !$omp simd reduction(+:sum_up, sum_dn)
-        do jg = 1,ng
-          sum_up = sum_up + flux_up_clear(jg,jlev)
-          sum_dn = sum_dn + flux_dn_clear(jg,jlev)
-        end do
-        flux%lw_up_clear(jcol,jlev) = sum_up
-        flux%lw_dn_clear(jcol,jlev) = sum_dn
-      end do
-#endif
-
+      flux%lw_up_clear(jcol,:) = sum(flux_up_clear,1)
+      flux%lw_dn_clear(jcol,:) = sum(flux_dn_clear,1)
       ! Store surface spectral downwelling fluxes
       flux%lw_dn_surf_clear_g(:,jcol) = flux_dn_clear(:,nlev+1)
@@ -310,10 +279,8 @@
           else
             ! Clear-sky layer: copy over clear-sky values
-            do jg = 1,ng
-              reflectance(jg,jlev) = ref_clear(jg,jlev)
-              transmittance(jg,jlev) = trans_clear(jg,jlev)
-              source_up(jg,jlev) = source_up_clear(jg,jlev)
-              source_dn(jg,jlev) = source_dn_clear(jg,jlev)
-            end do
+            reflectance(:,jlev) = ref_clear(:,jlev)
+            transmittance(:,jlev) = trans_clear(:,jlev)
+            source_up(:,jlev) = source_up_clear(:,jlev)
+            source_dn(:,jlev) = source_dn_clear(:,jlev)
           end if
         end do
@@ -340,36 +307,13 @@
         
         ! Store overcast broadband fluxes
-#ifdef DWD_REDUCTION_OPTIMIZATIONS
-        sum_aux(:,:) = 0._jprb
-        do jg = 1, ng
-          do jlev = 1, nlev+1
-            sum_aux(jlev,1) = sum_aux(jlev,1) + flux_up(jg,jlev)
-            sum_aux(jlev,2) = sum_aux(jlev,2) + flux_dn(jg,jlev)
-          end do
-        end do
-        flux%lw_up(jcol,:) = sum_aux(:,1)
-        flux%lw_dn(jcol,:) = sum_aux(:,2)
-#else
-        do jlev = 1,nlev+1
-          sum_up = 0.0_jprb
-          sum_dn = 0.0_jprb
-          !$omp simd reduction(+:sum_up, sum_dn)
-          do jg = 1,ng
-            sum_up = sum_up + flux_up(jg,jlev)
-            sum_dn = sum_dn + flux_dn(jg,jlev)
-          end do
-          flux%lw_up(jcol,jlev) = sum_up
-          flux%lw_dn(jcol,jlev) = sum_dn
-        end do
-#endif
+        flux%lw_up(jcol,:) = sum(flux_up,1)
+        flux%lw_dn(jcol,:) = sum(flux_dn,1)
 
         ! Cloudy flux profiles currently assume completely overcast
         ! skies; perform weighted average with clear-sky profile
-        do jlev = 1,nlev+1
-          flux%lw_up(jcol,jlev) =  total_cloud_cover *flux%lw_up(jcol,jlev) &
-             &       + (1.0_jprb - total_cloud_cover)*flux%lw_up_clear(jcol,jlev)
-          flux%lw_dn(jcol,jlev) =  total_cloud_cover *flux%lw_dn(jcol,jlev) &
-             &       + (1.0_jprb - total_cloud_cover)*flux%lw_dn_clear(jcol,jlev)
-        end do
+        flux%lw_up(jcol,:) =  total_cloud_cover *flux%lw_up(jcol,:) &
+             &  + (1.0_jprb - total_cloud_cover)*flux%lw_up_clear(jcol,:)
+        flux%lw_dn(jcol,:) =  total_cloud_cover *flux%lw_dn(jcol,:) &
+             &  + (1.0_jprb - total_cloud_cover)*flux%lw_dn_clear(jcol,:)
         ! Store surface spectral downwelling fluxes
         flux%lw_dn_surf_g(:,jcol) = total_cloud_cover*flux_dn(:,nlev+1) &
@@ -391,8 +335,6 @@
         ! No cloud in profile and clear-sky fluxes already
         ! calculated: copy them over
-        do jlev = 1,nlev+1
-          flux%lw_up(jcol,jlev) = flux%lw_up_clear(jcol,jlev)
-          flux%lw_dn(jcol,jlev) = flux%lw_dn_clear(jcol,jlev)
-        end do
+        flux%lw_up(jcol,:) = flux%lw_up_clear(jcol,:)
+        flux%lw_dn(jcol,:) = flux%lw_dn_clear(jcol,:)
         flux%lw_dn_surf_g(:,jcol) = flux%lw_dn_surf_clear_g(:,jcol)
         if (config%do_lw_derivatives) then
Index: LMDZ6/branches/cirrus/libf/phylmd/ecrad/radiation/radiation_mcica_sw.F90
===================================================================
--- LMDZ6/branches/cirrus/libf/phylmd/ecrad/radiation/radiation_mcica_sw.F90	(revision 4945)
+++ LMDZ6/branches/cirrus/libf/phylmd/ecrad/radiation/radiation_mcica_sw.F90	(revision 5202)
@@ -17,6 +17,4 @@
 !   2017-04-22  R. Hogan  Store surface fluxes at all g-points
 !   2017-10-23  R. Hogan  Renamed single-character variables
-
-#include "ecrad_config.h"
 
 module radiation_mcica_sw
@@ -121,11 +119,4 @@
     ! Total cloud cover output from the cloud generator
     real(jprb) :: total_cloud_cover
-
-    ! Temporary storage for more efficient summation
-#ifdef DWD_REDUCTION_OPTIMIZATIONS
-    real(jprb), dimension(nlev+1,3) :: sum_aux
-#else
-    real(jprb) :: sum_up, sum_dn_diff, sum_dn_dir
-#endif
 
     ! Number of g points
@@ -184,47 +175,18 @@
         
         ! Sum over g-points to compute and save clear-sky broadband
-        ! fluxes. Note that the built-in "sum" function is very slow,
-        ! and before being replaced by the alternatives below
-        ! accounted for around 40% of the total cost of this routine.
-#ifdef DWD_REDUCTION_OPTIMIZATIONS
-        ! Optimized summation for the NEC architecture
-        sum_aux(:,:) = 0.0_jprb
-        do jg = 1,ng
-          do jlev = 1,nlev+1
-            sum_aux(jlev,1) = sum_aux(jlev,1) + flux_up(jg,jlev)
-            sum_aux(jlev,2) = sum_aux(jlev,2) + flux_dn_direct(jg,jlev)
-            sum_aux(jlev,3) = sum_aux(jlev,3) + flux_dn_diffuse(jg,jlev)
-          end do
-        end do
-        flux%sw_up_clear(jcol,:) = sum_aux(:,1)
-        flux%sw_dn_clear(jcol,:) = sum_aux(:,2) + sum_aux(:,3)
+        ! fluxes
+        flux%sw_up_clear(jcol,:) = sum(flux_up,1)
         if (allocated(flux%sw_dn_direct_clear)) then
-          flux%sw_dn_direct_clear(jcol,:) = sum_aux(:,2)
+          flux%sw_dn_direct_clear(jcol,:) &
+               &  = sum(flux_dn_direct,1)
+          flux%sw_dn_clear(jcol,:) = sum(flux_dn_diffuse,1) &
+               &  + flux%sw_dn_direct_clear(jcol,:)
+        else
+          flux%sw_dn_clear(jcol,:) = sum(flux_dn_diffuse,1) &
+               &  + sum(flux_dn_direct,1)
         end if
-#else
-        ! Optimized summation for the x86-64 architecture
-        do jlev = 1,nlev+1
-          sum_up      = 0.0_jprb
-          sum_dn_diff = 0.0_jprb
-          sum_dn_dir  = 0.0_jprb
-          !$omp simd reduction(+:sum_up, sum_dn_diff, sum_dn_dir)
-          do jg = 1,ng
-            sum_up      = sum_up      + flux_up(jg,jlev)
-            sum_dn_diff = sum_dn_diff + flux_dn_diffuse(jg,jlev)
-            sum_dn_dir  = sum_dn_dir  + flux_dn_direct(jg,jlev)
-          end do
-          flux%sw_up_clear(jcol,jlev) = sum_up
-          flux%sw_dn_clear(jcol,jlev) = sum_dn_diff + sum_dn_dir
-          if (allocated(flux%sw_dn_direct_clear)) then
-            flux%sw_dn_direct_clear(jcol,jlev) = sum_dn_dir
-          end if
-        end do
-#endif
-        
         ! Store spectral downwelling fluxes at surface
-        do jg = 1,ng
-          flux%sw_dn_diffuse_surf_clear_g(jg,jcol) = flux_dn_diffuse(jg,nlev+1)
-          flux%sw_dn_direct_surf_clear_g(jg,jcol)  = flux_dn_direct(jg,nlev+1)
-        end do
+        flux%sw_dn_diffuse_surf_clear_g(:,jcol) = flux_dn_diffuse(:,nlev+1)
+        flux%sw_dn_direct_surf_clear_g(:,jcol)  = flux_dn_direct(:,nlev+1)
 
         ! Do cloudy-sky calculation
@@ -287,11 +249,9 @@
             else
               ! Clear-sky layer: copy over clear-sky values
-              do jg = 1,ng
-                reflectance(jg,jlev) = ref_clear(jg,jlev)
-                transmittance(jg,jlev) = trans_clear(jg,jlev)
-                ref_dir(jg,jlev) = ref_dir_clear(jg,jlev)
-                trans_dir_diff(jg,jlev) = trans_dir_diff_clear(jg,jlev)
-                trans_dir_dir(jg,jlev) = trans_dir_dir_clear(jg,jlev)
-              end do
+              reflectance(:,jlev) = ref_clear(:,jlev)
+              transmittance(:,jlev) = trans_clear(:,jlev)
+              ref_dir(:,jlev) = ref_dir_clear(:,jlev)
+              trans_dir_diff(:,jlev) = trans_dir_diff_clear(:,jlev)
+              trans_dir_dir(:,jlev) = trans_dir_dir_clear(:,jlev)
             end if
           end do
@@ -304,73 +264,42 @@
           
           ! Store overcast broadband fluxes
-#ifdef DWD_REDUCTION_OPTIMIZATIONS
-          sum_aux(:,:) = 0.0_jprb
-          do jg = 1,ng
-            do jlev = 1,nlev+1
-              sum_aux(jlev,1) = sum_aux(jlev,1) + flux_up(jg,jlev)
-              sum_aux(jlev,2) = sum_aux(jlev,2) + flux_dn_direct(jg,jlev)
-              sum_aux(jlev,3) = sum_aux(jlev,3) + flux_dn_diffuse(jg,jlev)
-            end do
-          end do
-          flux%sw_up(jcol,:) = sum_aux(:,1)
-          flux%sw_dn(jcol,:) = sum_aux(:,2) + sum_aux(:,3)
+          flux%sw_up(jcol,:) = sum(flux_up,1)
           if (allocated(flux%sw_dn_direct)) then
-            flux%sw_dn_direct(jcol,:) = sum_aux(:,2)
+            flux%sw_dn_direct(jcol,:) = sum(flux_dn_direct,1)
+            flux%sw_dn(jcol,:) = sum(flux_dn_diffuse,1) &
+                 &  + flux%sw_dn_direct(jcol,:)
+          else
+            flux%sw_dn(jcol,:) = sum(flux_dn_diffuse,1) &
+                 &  + sum(flux_dn_direct,1)
           end if
-#else
-          do jlev = 1,nlev+1
-            sum_up      = 0.0_jprb
-            sum_dn_diff = 0.0_jprb
-            sum_dn_dir  = 0.0_jprb
-            !$omp simd reduction(+:sum_up, sum_dn_diff, sum_dn_dir)
-            do jg = 1,ng
-              sum_up      = sum_up      + flux_up(jg,jlev)
-              sum_dn_diff = sum_dn_diff + flux_dn_diffuse(jg,jlev)
-              sum_dn_dir  = sum_dn_dir  + flux_dn_direct(jg,jlev)
-            end do
-            flux%sw_up(jcol,jlev) = sum_up
-            flux%sw_dn(jcol,jlev) = sum_dn_diff + sum_dn_dir
-            if (allocated(flux%sw_dn_direct)) then
-              flux%sw_dn_direct(jcol,jlev) = sum_dn_dir
-            end if
-          end do
-#endif
-          
+
           ! Cloudy flux profiles currently assume completely overcast
           ! skies; perform weighted average with clear-sky profile
-          do jlev = 1, nlev+1
-            flux%sw_up(jcol,jlev) =  total_cloud_cover *flux%sw_up(jcol,jlev) &
-                 &     + (1.0_jprb - total_cloud_cover)*flux%sw_up_clear(jcol,jlev)
-            flux%sw_dn(jcol,jlev) =  total_cloud_cover *flux%sw_dn(jcol,jlev) &
-                 &     + (1.0_jprb - total_cloud_cover)*flux%sw_dn_clear(jcol,jlev)
-            if (allocated(flux%sw_dn_direct)) then
-              flux%sw_dn_direct(jcol,jlev) = total_cloud_cover *flux%sw_dn_direct(jcol,jlev) &
-                   &  + (1.0_jprb - total_cloud_cover)*flux%sw_dn_direct_clear(jcol,jlev)
-            end if
-          end do
+          flux%sw_up(jcol,:) =  total_cloud_cover *flux%sw_up(jcol,:) &
+               &  + (1.0_jprb - total_cloud_cover)*flux%sw_up_clear(jcol,:)
+          flux%sw_dn(jcol,:) =  total_cloud_cover *flux%sw_dn(jcol,:) &
+               &  + (1.0_jprb - total_cloud_cover)*flux%sw_dn_clear(jcol,:)
+          if (allocated(flux%sw_dn_direct)) then
+            flux%sw_dn_direct(jcol,:) = total_cloud_cover *flux%sw_dn_direct(jcol,:) &
+                 &  + (1.0_jprb - total_cloud_cover)*flux%sw_dn_direct_clear(jcol,:)
+          end if
           ! Likewise for surface spectral fluxes
-          do jg = 1,ng
-            flux%sw_dn_diffuse_surf_g(jg,jcol) = flux_dn_diffuse(jg,nlev+1)
-            flux%sw_dn_direct_surf_g(jg,jcol)  = flux_dn_direct(jg,nlev+1)
-            flux%sw_dn_diffuse_surf_g(jg,jcol) = total_cloud_cover *flux%sw_dn_diffuse_surf_g(jg,jcol) &
-                 &                 + (1.0_jprb - total_cloud_cover)*flux%sw_dn_diffuse_surf_clear_g(jg,jcol)
-            flux%sw_dn_direct_surf_g(jg,jcol)  = total_cloud_cover *flux%sw_dn_direct_surf_g(jg,jcol) &
-                 &                 + (1.0_jprb - total_cloud_cover)*flux%sw_dn_direct_surf_clear_g(jg,jcol)
-          end do
-
+          flux%sw_dn_diffuse_surf_g(:,jcol) = flux_dn_diffuse(:,nlev+1)
+          flux%sw_dn_direct_surf_g(:,jcol)  = flux_dn_direct(:,nlev+1)
+          flux%sw_dn_diffuse_surf_g(:,jcol) = total_cloud_cover *flux%sw_dn_diffuse_surf_g(:,jcol) &
+               &     + (1.0_jprb - total_cloud_cover)*flux%sw_dn_diffuse_surf_clear_g(:,jcol)
+          flux%sw_dn_direct_surf_g(:,jcol) = total_cloud_cover *flux%sw_dn_direct_surf_g(:,jcol) &
+               &     + (1.0_jprb - total_cloud_cover)*flux%sw_dn_direct_surf_clear_g(:,jcol)
+          
         else
           ! No cloud in profile and clear-sky fluxes already
           ! calculated: copy them over
-          do jlev = 1, nlev+1
-            flux%sw_up(jcol,jlev) = flux%sw_up_clear(jcol,jlev)
-            flux%sw_dn(jcol,jlev) = flux%sw_dn_clear(jcol,jlev)
-            if (allocated(flux%sw_dn_direct)) then
-              flux%sw_dn_direct(jcol,jlev) = flux%sw_dn_direct_clear(jcol,jlev)
-            end if
-          end do
-          do jg = 1,ng
-            flux%sw_dn_diffuse_surf_g(jg,jcol) = flux%sw_dn_diffuse_surf_clear_g(jg,jcol)
-            flux%sw_dn_direct_surf_g(jg,jcol)  = flux%sw_dn_direct_surf_clear_g(jg,jcol)
-          end do
+          flux%sw_up(jcol,:) = flux%sw_up_clear(jcol,:)
+          flux%sw_dn(jcol,:) = flux%sw_dn_clear(jcol,:)
+          if (allocated(flux%sw_dn_direct)) then
+            flux%sw_dn_direct(jcol,:) = flux%sw_dn_direct_clear(jcol,:)
+          end if
+          flux%sw_dn_diffuse_surf_g(:,jcol) = flux%sw_dn_diffuse_surf_clear_g(:,jcol)
+          flux%sw_dn_direct_surf_g(:,jcol)  = flux%sw_dn_direct_surf_clear_g(:,jcol)
 
         end if ! Cloud is present in profile
@@ -378,22 +307,18 @@
       else
         ! Set fluxes to zero if sun is below the horizon
-        do jlev = 1, nlev+1
-          flux%sw_up(jcol,jlev) = 0.0_jprb
-          flux%sw_dn(jcol,jlev) = 0.0_jprb
-          if (allocated(flux%sw_dn_direct)) then
-            flux%sw_dn_direct(jcol,jlev) = 0.0_jprb
-          end if
-          flux%sw_up_clear(jcol,jlev) = 0.0_jprb
-          flux%sw_dn_clear(jcol,jlev) = 0.0_jprb
-          if (allocated(flux%sw_dn_direct_clear)) then
-            flux%sw_dn_direct_clear(jcol,jlev) = 0.0_jprb
-          end if
-        end do
-        do jg = 1,ng
-          flux%sw_dn_diffuse_surf_g(jg,jcol) = 0.0_jprb
-          flux%sw_dn_direct_surf_g(jg,jcol)  = 0.0_jprb
-          flux%sw_dn_diffuse_surf_clear_g(jg,jcol) = 0.0_jprb
-          flux%sw_dn_direct_surf_clear_g(jg,jcol)  = 0.0_jprb
-        end do
+        flux%sw_up(jcol,:) = 0.0_jprb
+        flux%sw_dn(jcol,:) = 0.0_jprb
+        if (allocated(flux%sw_dn_direct)) then
+          flux%sw_dn_direct(jcol,:) = 0.0_jprb
+        end if
+        flux%sw_up_clear(jcol,:) = 0.0_jprb
+        flux%sw_dn_clear(jcol,:) = 0.0_jprb
+        if (allocated(flux%sw_dn_direct_clear)) then
+          flux%sw_dn_direct_clear(jcol,:) = 0.0_jprb
+        end if
+        flux%sw_dn_diffuse_surf_g(:,jcol) = 0.0_jprb
+        flux%sw_dn_direct_surf_g(:,jcol)  = 0.0_jprb
+        flux%sw_dn_diffuse_surf_clear_g(:,jcol) = 0.0_jprb
+        flux%sw_dn_direct_surf_clear_g(:,jcol)  = 0.0_jprb
       end if ! Sun above horizon
 
Index: LMDZ6/branches/cirrus/libf/phylmd/ecrad/radiation/radiation_tripleclouds_lw.F90
===================================================================
--- LMDZ6/branches/cirrus/libf/phylmd/ecrad/radiation/radiation_tripleclouds_lw.F90	(revision 4945)
+++ LMDZ6/branches/cirrus/libf/phylmd/ecrad/radiation/radiation_tripleclouds_lw.F90	(revision 5202)
@@ -170,7 +170,4 @@
     logical :: is_clear_sky_layer(0:nlev+1)
 
-    ! Temporaries to speed up summations
-    real(jprb) :: sum_dn, sum_up
-    
     ! Index of the highest cloudy layer
     integer :: i_cloud_top
@@ -264,21 +261,9 @@
       if (config%do_clear) then
         ! Sum over g-points to compute broadband fluxes
-        do jlev = 1,nlev+1
-          sum_up = 0.0_jprb
-          sum_dn = 0.0_jprb
-          !$omp simd reduction(+:sum_up, sum_dn)
-          do jg = 1,ng
-            sum_up = sum_up + flux_up_clear(jg,jlev)
-            sum_dn = sum_dn + flux_dn_clear(jg,jlev)
-          end do
-          flux%lw_up_clear(jcol,jlev) = sum_up
-          flux%lw_dn_clear(jcol,jlev) = sum_dn
-        end do
-
+        flux%lw_up_clear(jcol,:) = sum(flux_up_clear,1)
+        flux%lw_dn_clear(jcol,:) = sum(flux_dn_clear,1)
         ! Store surface spectral downwelling fluxes / TOA upwelling
-        do jg = 1,ng
-          flux%lw_dn_surf_clear_g(jg,jcol) = flux_dn_clear(jg,nlev+1)
-          flux%lw_up_toa_clear_g (jg,jcol) = flux_up_clear(jg,1)
-        end do
+        flux%lw_dn_surf_clear_g(:,jcol) = flux_dn_clear(:,nlev+1)
+        flux%lw_up_toa_clear_g (:,jcol) = flux_up_clear(:,1)
         ! Save the spectral fluxes if required
         if (config%do_save_spectral_flux) then
@@ -468,10 +453,5 @@
           end if
         else
-          sum_dn = 0.0_jprb
-          !$omp simd reduction(+:sum_dn)
-          do jg = 1,ng
-            sum_dn = sum_dn + flux_dn_clear(jg,jlev)
-          end do
-          flux%lw_dn(jcol,jlev) = sum_dn
+          flux%lw_dn(jcol,:) = sum(flux_dn_clear(:,jlev))
           if (config%do_save_spectral_flux) then
             call indexed_sum(flux_dn_clear(:,jlev), &
@@ -490,12 +470,5 @@
            &  + total_albedo(:,1,i_cloud_top)*flux_dn_clear(:,i_cloud_top)
       flux_up(:,2:) = 0.0_jprb
-
-      sum_up = 0.0_jprb
-      !$omp simd reduction(+:sum_up)
-      do jg = 1,ng
-        sum_up = sum_up + flux_up(jg,1)
-      end do
-      flux%lw_up(jcol,i_cloud_top) = sum_up
-
+      flux%lw_up(jcol,i_cloud_top) = sum(flux_up(:,1))
       if (config%do_save_spectral_flux) then
         call indexed_sum(flux_up(:,1), &
@@ -505,10 +478,5 @@
       do jlev = i_cloud_top-1,1,-1
         flux_up(:,1) = trans_clear(:,jlev)*flux_up(:,1) + source_up_clear(:,jlev)
-        sum_up = 0.0_jprb
-        !$omp simd reduction(+:sum_up)
-        do jg = 1,ng
-          sum_up = sum_up + flux_up(jg,1)
-        end do
-        flux%lw_up(jcol,jlev) = sum_up
+        flux%lw_up(jcol,jlev) = sum(flux_up(:,1))
         if (config%do_save_spectral_flux) then
           call indexed_sum(flux_up(:,1), &
@@ -560,15 +528,6 @@
 
         ! Store the broadband fluxes
-        sum_up = 0.0_jprb
-        sum_dn = 0.0_jprb
-        do jreg = 1,nregions
-          !$omp simd reduction(+:sum_up, sum_dn)
-          do jg = 1,ng
-            sum_up = sum_up + flux_up(jg,jreg)
-            sum_dn = sum_dn + flux_dn(jg,jreg)
-          end do
-        end do
-        flux%lw_up(jcol,jlev+1) = sum_up
-        flux%lw_dn(jcol,jlev+1) = sum_dn
+        flux%lw_up(jcol,jlev+1) = sum(sum(flux_up,1))
+        flux%lw_dn(jcol,jlev+1) = sum(sum(flux_dn,1))
 
         ! Save the spectral fluxes if required
Index: LMDZ6/branches/cirrus/libf/phylmd/ecrad/radiation/radiation_tripleclouds_lw.F90.or
===================================================================
--- LMDZ6/branches/cirrus/libf/phylmd/ecrad/radiation/radiation_tripleclouds_lw.F90.or	(revision 4945)
+++ LMDZ6/branches/cirrus/libf/phylmd/ecrad/radiation/radiation_tripleclouds_lw.F90.or	(revision 5202)
@@ -170,4 +170,7 @@
     logical :: is_clear_sky_layer(0:nlev+1)
 
+    ! Temporaries to speed up summations
+    real(jprb) :: sum_dn, sum_up
+    
     ! Index of the highest cloudy layer
     integer :: i_cloud_top
@@ -249,5 +252,5 @@
         call calc_ref_trans_lw(ng*nlev, &
              &  od(:,:,jcol), ssa(:,:,jcol), g(:,:,jcol), &
-             &  planck_hl(:,1:jlev,jcol), planck_hl(:,2:jlev+1,jcol), &
+             &  planck_hl(:,1:nlev,jcol), planck_hl(:,2:nlev+1,jcol), &
              &  ref_clear, trans_clear, &
              &  source_up_clear, source_dn_clear)
@@ -261,9 +264,21 @@
       if (config%do_clear) then
         ! Sum over g-points to compute broadband fluxes
-        flux%lw_up_clear(jcol,:) = sum(flux_up_clear,1)
-        flux%lw_dn_clear(jcol,:) = sum(flux_dn_clear,1)
+        do jlev = 1,nlev+1
+          sum_up = 0.0_jprb
+          sum_dn = 0.0_jprb
+          !$omp simd reduction(+:sum_up, sum_dn)
+          do jg = 1,ng
+            sum_up = sum_up + flux_up_clear(jg,jlev)
+            sum_dn = sum_dn + flux_dn_clear(jg,jlev)
+          end do
+          flux%lw_up_clear(jcol,jlev) = sum_up
+          flux%lw_dn_clear(jcol,jlev) = sum_dn
+        end do
+
         ! Store surface spectral downwelling fluxes / TOA upwelling
-        flux%lw_dn_surf_clear_g(:,jcol) = flux_dn_clear(:,nlev+1)
-        flux%lw_up_toa_clear_g (:,jcol) = flux_up_clear(:,1)
+        do jg = 1,ng
+          flux%lw_dn_surf_clear_g(jg,jcol) = flux_dn_clear(jg,nlev+1)
+          flux%lw_up_toa_clear_g (jg,jcol) = flux_up_clear(jg,1)
+        end do
         ! Save the spectral fluxes if required
         if (config%do_save_spectral_flux) then
@@ -453,5 +468,10 @@
           end if
         else
-          flux%lw_dn(jcol,:) = sum(flux_dn_clear(:,jlev))
+          sum_dn = 0.0_jprb
+          !$omp simd reduction(+:sum_dn)
+          do jg = 1,ng
+            sum_dn = sum_dn + flux_dn_clear(jg,jlev)
+          end do
+          flux%lw_dn(jcol,jlev) = sum_dn
           if (config%do_save_spectral_flux) then
             call indexed_sum(flux_dn_clear(:,jlev), &
@@ -470,5 +490,12 @@
            &  + total_albedo(:,1,i_cloud_top)*flux_dn_clear(:,i_cloud_top)
       flux_up(:,2:) = 0.0_jprb
-      flux%lw_up(jcol,i_cloud_top) = sum(flux_up(:,1))
+
+      sum_up = 0.0_jprb
+      !$omp simd reduction(+:sum_up)
+      do jg = 1,ng
+        sum_up = sum_up + flux_up(jg,1)
+      end do
+      flux%lw_up(jcol,i_cloud_top) = sum_up
+
       if (config%do_save_spectral_flux) then
         call indexed_sum(flux_up(:,1), &
@@ -478,5 +505,10 @@
       do jlev = i_cloud_top-1,1,-1
         flux_up(:,1) = trans_clear(:,jlev)*flux_up(:,1) + source_up_clear(:,jlev)
-        flux%lw_up(jcol,jlev) = sum(flux_up(:,1))
+        sum_up = 0.0_jprb
+        !$omp simd reduction(+:sum_up)
+        do jg = 1,ng
+          sum_up = sum_up + flux_up(jg,1)
+        end do
+        flux%lw_up(jcol,jlev) = sum_up
         if (config%do_save_spectral_flux) then
           call indexed_sum(flux_up(:,1), &
@@ -528,6 +560,15 @@
 
         ! Store the broadband fluxes
-        flux%lw_up(jcol,jlev+1) = sum(sum(flux_up,1))
-        flux%lw_dn(jcol,jlev+1) = sum(sum(flux_dn,1))
+        sum_up = 0.0_jprb
+        sum_dn = 0.0_jprb
+        do jreg = 1,nregions
+          !$omp simd reduction(+:sum_up, sum_dn)
+          do jg = 1,ng
+            sum_up = sum_up + flux_up(jg,jreg)
+            sum_dn = sum_dn + flux_dn(jg,jreg)
+          end do
+        end do
+        flux%lw_up(jcol,jlev+1) = sum_up
+        flux%lw_dn(jcol,jlev+1) = sum_dn
 
         ! Save the spectral fluxes if required
Index: LMDZ6/branches/cirrus/libf/phylmd/ecrad/radiation/radiation_tripleclouds_sw.F90
===================================================================
--- LMDZ6/branches/cirrus/libf/phylmd/ecrad/radiation/radiation_tripleclouds_sw.F90	(revision 4945)
+++ LMDZ6/branches/cirrus/libf/phylmd/ecrad/radiation/radiation_tripleclouds_sw.F90	(revision 5202)
@@ -74,11 +74,12 @@
     ! Gas and aerosol optical depth, single-scattering albedo and
     ! asymmetry factor at each shortwave g-point
-    real(jprb), intent(in), dimension(config%n_g_sw,nlev,istartcol:iendcol) &
-         &  :: od, ssa, g
+!    real(jprb), intent(in), dimension(istartcol:iendcol,nlev,config%n_g_sw) :: &
+    real(jprb), intent(in), dimension(config%n_g_sw,nlev,istartcol:iendcol) :: &
+         &  od, ssa, g
 
     ! Cloud and precipitation optical depth, single-scattering albedo and
     ! asymmetry factor in each shortwave band
-    real(jprb), intent(in), dimension(config%n_bands_sw,nlev,istartcol:iendcol) &
-         &  :: od_cloud, ssa_cloud, g_cloud
+    real(jprb), intent(in), dimension(config%n_bands_sw,nlev,istartcol:iendcol) :: &
+         &  od_cloud, ssa_cloud, g_cloud
 
     ! Optical depth, single scattering albedo and asymmetry factor in
@@ -91,6 +92,6 @@
     ! flux into a plane perpendicular to the incoming radiation at
     ! top-of-atmosphere in each of the shortwave g points
-    real(jprb), intent(in), dimension(config%n_g_sw,istartcol:iendcol) &
-         &  :: albedo_direct, albedo_diffuse, incoming_sw
+    real(jprb), intent(in), dimension(config%n_g_sw,istartcol:iendcol) :: &
+         &  albedo_direct, albedo_diffuse, incoming_sw
 
     ! Output
@@ -165,8 +166,4 @@
     real(jprb) :: scat_od, scat_od_cloud
 
-    ! Temporaries to speed up summations
-    real(jprb) :: sum_dn_diff, sum_dn_dir, sum_up
-
-    ! Local cosine of solar zenith angle
     real(jprb) :: mu0
 
@@ -447,33 +444,13 @@
       end if
       
-      ! Store the TOA broadband fluxes, noting that there is no
-      ! diffuse downwelling at TOA. The intrinsic "sum" command has
-      ! been found to be very slow; better performance is found on
-      ! x86-64 architecture with explicit loops and the "omp simd
-      ! reduction" directive.
-      sum_up     = 0.0_jprb
-      sum_dn_dir = 0.0_jprb
-      do jreg = 1,nregions
-        !$omp simd reduction(+:sum_up, sum_dn_dir)
-        do jg = 1,ng
-          sum_up     = sum_up     + flux_up(jg,jreg)
-          sum_dn_dir = sum_dn_dir + direct_dn(jg,jreg)
-        end do
-      end do
-      flux%sw_up(jcol,1) = sum_up
-      flux%sw_dn(jcol,1) = mu0 * sum_dn_dir
+      ! Store the TOA broadband fluxes
+      flux%sw_up(jcol,1) = sum(sum(flux_up,1))
+      flux%sw_dn(jcol,1) = mu0 * sum(sum(direct_dn,1))
       if (allocated(flux%sw_dn_direct)) then
         flux%sw_dn_direct(jcol,1) = flux%sw_dn(jcol,1)
       end if
       if (config%do_clear) then
-        sum_up     = 0.0_jprb
-        sum_dn_dir = 0.0_jprb
-        !$omp simd reduction(+:sum_up, sum_dn_dir)
-        do jg = 1,ng
-          sum_up     = sum_up     + flux_up_clear(jg)
-          sum_dn_dir = sum_dn_dir + direct_dn_clear(jg)
-        end do
-        flux%sw_up_clear(jcol,1) = sum_up
-        flux%sw_dn_clear(jcol,1) = mu0 * sum_dn_dir
+        flux%sw_up_clear(jcol,1) = sum(flux_up_clear)
+        flux%sw_dn_clear(jcol,1) = mu0 * sum(direct_dn_clear)
         if (allocated(flux%sw_dn_direct_clear)) then
           flux%sw_dn_direct_clear(jcol,1) = flux%sw_dn_clear(jcol,1)
@@ -490,5 +467,6 @@
              &           config%i_spec_from_reordered_g_sw, &
              &           flux%sw_dn_band(:,jcol,1))
-        flux%sw_dn_band(:,jcol,1) = mu0 * flux%sw_dn_band(:,jcol,1)
+        flux%sw_dn_band(:,jcol,1) = &
+             &  mu0 * flux%sw_dn_band(:,jcol,1)
         if (allocated(flux%sw_dn_direct_band)) then
           flux%sw_dn_direct_band(:,jcol,1) = flux%sw_dn_band(:,jcol,1)
@@ -571,38 +549,22 @@
                ! nothing to do
 
-        ! Store the broadband fluxes. The intrinsic "sum" command has
-        ! been found to be very slow; better performance is found on
-        ! x86-64 architecture with explicit loops and the "omp simd
-        ! reduction" directive.
-        sum_up      = 0.0_jprb
-        sum_dn_dir  = 0.0_jprb
-        sum_dn_diff = 0.0_jprb
-        do jreg = 1,nregions
-          !$omp simd reduction(+:sum_up, sum_dn_diff, sum_dn_dir)
-          do jg = 1,ng
-            sum_up      = sum_up      + flux_up(jg,jreg)
-            sum_dn_diff = sum_dn_diff + flux_dn(jg,jreg)
-            sum_dn_dir  = sum_dn_dir  + direct_dn(jg,jreg)
-          end do
-        end do
-        flux%sw_up(jcol,jlev+1) = sum_up
-        flux%sw_dn(jcol,jlev+1) = mu0 * sum_dn_dir + sum_dn_diff
+        ! Store the broadband fluxes
+        flux%sw_up(jcol,jlev+1) = sum(sum(flux_up,1))
         if (allocated(flux%sw_dn_direct)) then
-          flux%sw_dn_direct(jcol,jlev+1) = mu0 * sum_dn_dir
+          flux%sw_dn_direct(jcol,jlev+1) = mu0 * sum(sum(direct_dn,1))
+          flux%sw_dn(jcol,jlev+1) &
+               &  = flux%sw_dn_direct(jcol,jlev+1) + sum(sum(flux_dn,1))
+        else
+          flux%sw_dn(jcol,jlev+1) = mu0 * sum(sum(direct_dn,1)) + sum(sum(flux_dn,1))   
         end if
         if (config%do_clear) then
-          sum_up      = 0.0_jprb
-          sum_dn_dir  = 0.0_jprb
-          sum_dn_diff = 0.0_jprb
-          !$omp simd reduction(+:sum_up, sum_dn_diff, sum_dn_dir)
-          do jg = 1,ng
-            sum_up      = sum_up      + flux_up_clear(jg)
-            sum_dn_diff = sum_dn_diff + flux_dn_clear(jg)
-            sum_dn_dir  = sum_dn_dir  + direct_dn_clear(jg)
-          end do
-          flux%sw_up_clear(jcol,jlev+1) = sum_up
-          flux%sw_dn_clear(jcol,jlev+1) = mu0 * sum_dn_dir + sum_dn_diff
+          flux%sw_up_clear(jcol,jlev+1) = sum(flux_up_clear)
           if (allocated(flux%sw_dn_direct_clear)) then
-            flux%sw_dn_direct_clear(jcol,jlev+1) = mu0 * sum_dn_dir
+            flux%sw_dn_direct_clear(jcol,jlev+1) = mu0 * sum(direct_dn_clear)
+            flux%sw_dn_clear(jcol,jlev+1) &
+                 &  = flux%sw_dn_direct_clear(jcol,jlev+1) + sum(flux_dn_clear)
+          else
+            flux%sw_dn_clear(jcol,jlev+1) = mu0 * sum(direct_dn_clear) &
+                 &  + sum(flux_dn_clear)
           end if
         end if
@@ -643,4 +605,5 @@
           end if
         end if
+
       end do ! Final loop over levels
       
