diff --git a/Makefile b/Makefile
index 24f09335..4b4fb951 100644
--- a/Makefile
+++ b/Makefile
@@ -44,6 +44,31 @@ ifdef SINGLE_PRECISION
 CPPFLAGS += -DPARKIND1_SINGLE
 endif
 
+# Option to declare number of g-points (inner dimension) at compile time, beneficial when NG is small   
+# make PROFILE=... NG_SW=32 NG_LW=32 (when using 32-term ECCKD models)
+ifdef NG_SW
+CPPFLAGS += -DNG_SW=$(NG_SW)
+endif
+ifdef NG_LW
+CPPFLAGS += -DNG_LW=$(NG_LW)
+endif
+
+# Use the General Purpose Timing Library (GPTL)?
+ifeq ($(GPTL_TIMING),1)
+CPPFLAGS += -DUSE_TIMING
+TIMING_INCLUDE = -I$(TIME_DIR)/include
+LDFLAGS += -L$(TIME_DIR)/lib -Wl,-rpath=$(TIME_DIR)/lib 
+LIBS_TIMING = -lgptlf -lgptl -rdynamic -lunwind # -rdynamic  
+else ifeq ($(GPTL_TIMING),2) 
+# Use GPTL with PAPI to estimate FLOPS and other performance metrics with hardware counters
+# only tested to work with gcc, other compilers probably need different flags
+CPPFLAGS += -DUSE_TIMING  -DUSE_PAPI
+TIMING_INCLUDE = -I$(TIME_DIR)/include 
+LDFLAGS  += -L$(TIME_DIR)/lib -Wl,-rpath=$(TIME_DIR)/lib 
+LIBS_TIMING = -lgptlf -lgptl  -rdynamic -lpapi -lunwind
+# LIBS_TIMING      += -lgptl -lpapi
+endif
+
 # If PRINT_ENTRAPMENT_DATA=1 was given on the "make" command line
 # then the SPARTACUS shortwave solver will write data to fort.101 and
 # fort.102
@@ -64,9 +89,9 @@ endif
 # Consolidate flags
 export FC
 export FCFLAGS = $(WARNFLAGS) $(BASICFLAGS) $(CPPFLAGS) -I../include \
-	$(OPTFLAGS) $(DEBUGFLAGS) $(NETCDF_INCLUDE) $(OMPFLAG)
+	$(OPTFLAGS) $(DEBUGFLAGS) $(NETCDF_INCLUDE) $(TIMING_INCLUDE) $(OMPFLAG)
 export LIBS    = $(LDFLAGS) -L../lib -lradiation -lutilities \
-	-lifsrrtm -lifsaux $(FCLIBS) $(NETCDF_LIB) $(OMPFLAG)
+	-lifsrrtm -lifsaux $(FCLIBS) $(NETCDF_LIB) $(LIBS_TIMING) $(OMPFLAG)
 
 # Do we include Dr Hook from ECMWF's fiat library?
 ifdef FIATDIR
diff --git a/Makefile_include.intel_atos b/Makefile_include.intel_atos
index 967b8504..fa4b29bb 100644
--- a/Makefile_include.intel_atos
+++ b/Makefile_include.intel_atos
@@ -38,3 +38,6 @@ BASICFLAGS = -module ../mod -fpe0 -assume byterecl -align array64byte \
 	   -finline-functions -finline-limit=1500 -Winline -assume realloc_lhs
 DEBUGFLAGS = -g
 OMPFLAG = -qopenmp -qopenmp-threadprivate=compat
+
+# location of General Purpose Timing Library (compile with GPTL_TIMING=1 or GPTL_TIMING=2)
+TIME_DIR=/home/papu/gptl-ifort
\ No newline at end of file
diff --git a/driver/ecrad_ifs_driver_blocked.F90 b/driver/ecrad_ifs_driver_blocked.F90
index 79492de9..98fd11c8 100644
--- a/driver/ecrad_ifs_driver_blocked.F90
+++ b/driver/ecrad_ifs_driver_blocked.F90
@@ -65,7 +65,16 @@ program ecrad_ifs_driver
   use ecrad_driver_read_input,  only : read_input
   use easy_netcdf
   use ifs_blocking
-
+#ifdef USE_TIMING
+#ifndef USE_PAPI
+  ! Timing library
+ use gptl,                  only: gptlstart, gptlstop, gptlinitialize, gptlpr, gptlfinalize, gptlsetoption, &
+                                  gptlpercent, gptloverhead, gptlpr_file
+#else
+#include "f90papi.h"
+#include "gptl.inc"
+#endif    
+#endif  
   implicit none
 
 #include "radiation_scheme.intfb.h"
@@ -137,7 +146,29 @@ program ecrad_ifs_driver
 
 !  integer    :: iband(20), nweights
 !  real(jprb) :: weight(20)
-
+#ifdef USE_TIMING
+  integer :: ret
+  integer values(8)
+  character(len=100) :: timing_file_name, name_gas_model, name_solver
+  character(5) :: prefix, suffix
+  !
+  ! Initialize timers
+  !
+  ret = gptlsetoption (gptlpercent, 1)        ! Turn on "% of" print
+  ret = gptlsetoption (gptloverhead, 0)       ! Turn off overhead estimate
+
+#ifdef USE_PAPI  
+#ifdef PARKIND1_SINGLE
+  ret = GPTLsetoption (PAPI_SP_OPS, 1);
+#else
+  ret = GPTLsetoption (PAPI_DP_OPS, 1);
+#endif
+! ret = GPTLsetoption (GPTL_IPC, 1);
+ret = GPTLsetoption (PAPI_L1_DCM, 1);
+ret = GPTLsetoption (GPTL_L3MRT, 1);
+#endif  
+  ret = gptlinitialize()
+#endif
 
   ! --------------------------------------------------------
   ! Section 2: Configure
@@ -365,6 +396,10 @@ program ecrad_ifs_driver
     write(nulout,'(a)')  'Performing radiative transfer calculations'
   end if
 
+#ifdef USE_TIMING
+    ret =  gptlstart('radiation')
+#endif   
+
   ! Option of repeating calculation multiple time for more accurate
   ! profiling
 #ifndef NO_OPENMP
@@ -444,7 +479,9 @@ program ecrad_ifs_driver
 !    end if
 
   end do
-
+#ifdef USE_TIMING
+  ret =  gptlstop('radiation')
+#endif
 #ifndef NO_OPENMP
   tstop = omp_get_wtime()
   write(nulout, '(a,g12.5,a)') 'Time elapsed in radiative transfer: ', tstop-tstart, ' seconds'
@@ -497,5 +534,15 @@ program ecrad_ifs_driver
   if (driver_config%iverbose >= 2) then
     write(nulout,'(a)') '------------------------------------------------------------------------------------'
   end if
-
+#ifdef USE_TIMING
+  ! End timers
+  call yradiation%rad_config%get_gas_optics_name_sw(name_gas_model)
+  call yradiation%rad_config%get_solver_name_sw(name_solver)
+  call date_and_time(values=values)
+  write(timing_file_name,'(a,a,a,a,a,i0,a,i0,a,i4.4)') 'timing.', trim(name_gas_model), '_', trim(name_solver), &
+    & '_block', driver_config%nblocksize, '_nrep', driver_config%nrepeat, '_', values(8)
+  write(nulout,'(a,a)')  'Writing GPTL timing output to ', trim(timing_file_name)
+  ret = gptlpr_file(trim(timing_file_name))
+  ret = gptlfinalize()
+#endif
 end program ecrad_ifs_driver
diff --git a/radiation/ecrad_config.h b/radiation/ecrad_config.h
index 3e4d9c3a..816d1fb8 100644
--- a/radiation/ecrad_config.h
+++ b/radiation/ecrad_config.h
@@ -37,3 +37,15 @@
 ! tasks. The MPI version is not used for writing files.
 
 !#define EASY_NETCDF_READ_MPI 1
+
+! Allow size of inner dimension (number of g-points) to be known at compile time if NG_LW/NG_SW is defined
+#ifdef NG_LW
+  integer, parameter :: ng_lw = NG_LW
+#else
+#define ng_lw ng_lw_in
+#endif
+#ifdef NG_SW
+  integer, parameter :: ng_sw = NG_SW
+#else
+#define ng_sw ng_sw_in
+#endif
\ No newline at end of file
diff --git a/radiation/ecrad_definitions.F90 b/radiation/ecrad_definitions.F90
new file mode 100644
index 00000000..d7a58f30
--- /dev/null
+++ b/radiation/ecrad_definitions.F90
@@ -0,0 +1,22 @@
+module ecrad_definitions
+
+  use parkind1, only : jpim
+
+  implicit none
+
+  ! By making the number of regions a parameter NREGION, the compiler
+  ! can optimize better.  We also provide a preprocessor parameter as
+  ! in one or two cases the operations on 2x2 and 3x3 matrices
+  ! (e.g. matrix exponentials) are completely different algorithms.
+#define NUM_REGIONS 2
+  integer(jpim), parameter :: NREGION = NUM_REGIONS
+
+! Allow size of inner dimension (number of g-points) to be known at compile time if NG_LW is defined
+#ifdef NG_LW
+  integer(jpim), parameter :: ng_lw = NG_LW
+#else
+#define ng_lw ng_lw_in
+#endif
+
+
+end module ecrad_definitions
\ No newline at end of file
diff --git a/radiation/radiation_aerosol_optics.F90 b/radiation/radiation_aerosol_optics.F90
index a705fb37..0f5f6cde 100644
--- a/radiation/radiation_aerosol_optics.F90
+++ b/radiation/radiation_aerosol_optics.F90
@@ -18,12 +18,11 @@
 !   2022-03-27  R. Hogan  Add setup_general_aerosol_optics_legacy to use RRTM aerosol files with ecCKD
 !   2022-11-22  P. Ukkonen / R. Hogan  Optimizations to enhance vectorization
 
-#include "ecrad_config.h"
-
 module radiation_aerosol_optics
 
   implicit none
   public
+#include "ecrad_config.h"
 
 contains
 
diff --git a/radiation/radiation_aerosol_optics_data.F90 b/radiation/radiation_aerosol_optics_data.F90
index d4a372e9..b652a6c6 100644
--- a/radiation/radiation_aerosol_optics_data.F90
+++ b/radiation/radiation_aerosol_optics_data.F90
@@ -16,8 +16,6 @@
 !   2017-10-23  R. Hogan  Renamed single-character variables
 !   2018-04-20  A. Bozzo  Read optical properties at selected wavelengths
 
-#include "ecrad_config.h"
-
 module radiation_aerosol_optics_data
 
   use parkind1,      only : jprb
@@ -28,6 +26,8 @@ module radiation_aerosol_optics_data
 
   private :: get_line
 
+#include "ecrad_config.h"
+
   ! The following provide possible values for
   ! aerosol_optics_type%iclass, which is used to map the user's type
   ! index to the hydrophobic or hydrophilic types loaded from the
diff --git a/radiation/radiation_aerosol_optics_description.F90 b/radiation/radiation_aerosol_optics_description.F90
index 50c7ba81..b370fa94 100644
--- a/radiation/radiation_aerosol_optics_description.F90
+++ b/radiation/radiation_aerosol_optics_description.F90
@@ -13,8 +13,6 @@
 ! Email:   r.j.hogan@ecmwf.int
 !
 
-#include "ecrad_config.h"
-
 module radiation_aerosol_optics_description
 
   use parkind1,      only : jprb
@@ -22,6 +20,8 @@ module radiation_aerosol_optics_description
   implicit none
   public
 
+#include "ecrad_config.h"
+
   !---------------------------------------------------------------------
   ! This type holds the metadata from an aerosol optical property
   ! file, enabling the user to request the index to the aerosol type
diff --git a/radiation/radiation_cloud_cover.F90 b/radiation/radiation_cloud_cover.F90
index 57cb93ba..219b2880 100644
--- a/radiation/radiation_cloud_cover.F90
+++ b/radiation/radiation_cloud_cover.F90
@@ -18,14 +18,14 @@
 ! Modifications
 !   2020-10-07  R. Hogan  Ensure iobj1 initialized in case of alpha_obj==0
 
-#include "ecrad_config.h"
-
 module radiation_cloud_cover
 
   use parkind1, only           : jprb
 
   public
 
+#include "ecrad_config.h"
+
   ! Three overlap schemes.  Note that "Exponential" means that
   ! clear-sky regions have no special significance for computing the
   ! cumulative cloud cover: non-contiguous clouds are exponentially
diff --git a/radiation/radiation_cloud_optics_data.F90 b/radiation/radiation_cloud_optics_data.F90
index 46d92ca4..4bf5812a 100644
--- a/radiation/radiation_cloud_optics_data.F90
+++ b/radiation/radiation_cloud_optics_data.F90
@@ -13,7 +13,7 @@
 ! Email:   r.j.hogan@ecmwf.int
 !
 
-#include "ecrad_config.h"
+! #include "ecrad_config.h"
 
 module radiation_cloud_optics_data
 
@@ -22,6 +22,7 @@ module radiation_cloud_optics_data
   implicit none
   public
 
+#include "ecrad_config.h"
   !---------------------------------------------------------------------
   ! This type holds the configuration information to compute
   ! cloud optical properties
diff --git a/radiation/radiation_config.F90 b/radiation/radiation_config.F90
index 8f8bccb4..2308af91 100644
--- a/radiation/radiation_config.F90
+++ b/radiation/radiation_config.F90
@@ -33,8 +33,6 @@
 ! files in this directory, please inform Robin Hogan.
 !
 
-#include "ecrad_config.h"
-
 module radiation_config
 
   use parkind1,                      only : jprb
@@ -49,6 +47,7 @@ module radiation_config
 
   implicit none
   public
+#include "ecrad_config.h"
 
   ! Configuration codes: use C-style enumerators to avoid having to
   ! remember the numbers
@@ -641,6 +640,8 @@ module radiation_config
      procedure :: set_aerosol_wavelength_mono
      procedure :: consolidate_sw_albedo_intervals
      procedure :: consolidate_lw_emiss_intervals
+     procedure :: get_gas_optics_name_sw
+     procedure :: get_solver_name_sw
 
   end type config_type
 
@@ -2135,4 +2136,16 @@ subroutine print_enum(message_str, enum_str, name, val)
     write(nulout,'(a,a,a,a,i0,a)') str, ' (', name, '=', val,')'
   end subroutine print_enum
 
+  subroutine get_gas_optics_name_sw(this, out_str)
+    class(config_type),   intent(inout) :: this
+    character(len=*),   intent(out) :: out_str
+    write(out_str,'(a)') trim(GasModelName(this%i_gas_model_sw))
+  end subroutine get_gas_optics_name_sw
+
+  subroutine get_solver_name_sw(this, out_str)
+    class(config_type),   intent(inout) :: this
+    character(len=*),   intent(out) :: out_str
+    write(out_str,'(a)') trim(SolverName(this%i_solver_sw))
+  end subroutine get_solver_name_sw
+
 end module radiation_config
diff --git a/radiation/radiation_ecckd.F90 b/radiation/radiation_ecckd.F90
index 2089b8de..85c6631c 100644
--- a/radiation/radiation_ecckd.F90
+++ b/radiation/radiation_ecckd.F90
@@ -14,8 +14,6 @@
 ! License: see the COPYING file for details
 !
 
-#include "ecrad_config.h"
-
 module radiation_ecckd
 
   use parkind1, only : jprb
@@ -478,9 +476,15 @@ subroutine calc_optical_depth_ckd_model(this, ncol, nlev, istartcol, iendcol, nm
     ! Output variables
 
     ! Layer absorption optical depth for each g point
+#if defined NG_SW && defined NG_LW && NG_SW==NG_LW 
+    real(jprb),            intent(out) :: optical_depth_fl(NG_SW,nlev,istartcol:iendcol)
+    ! In the shortwave only, the Rayleigh scattering optical depth
+    real(jprb),  optional, intent(out) :: rayleigh_od_fl(NG_SW,nlev,istartcol:iendcol)
+#else 
     real(jprb),            intent(out) :: optical_depth_fl(this%ng,nlev,istartcol:iendcol)
     ! In the shortwave only, the Rayleigh scattering optical depth
     real(jprb),  optional, intent(out) :: rayleigh_od_fl(this%ng,nlev,istartcol:iendcol)
+#endif
 
     ! Local variables
 
diff --git a/radiation/radiation_ecckd_gas.F90 b/radiation/radiation_ecckd_gas.F90
index a1a67053..b2481501 100644
--- a/radiation/radiation_ecckd_gas.F90
+++ b/radiation/radiation_ecckd_gas.F90
@@ -14,8 +14,6 @@
 ! License: see the COPYING file for details
 !
 
-#include "ecrad_config.h"
-
 module radiation_ecckd_gas
 
   use parkind1, only : jprb
@@ -25,6 +23,8 @@ module radiation_ecckd_gas
 
   public
 
+#include "ecrad_config.h"
+
   ! Concentration dependence of individual gases
   enum, bind(c)
     enumerator :: IConcDependenceNone = 0, &
diff --git a/radiation/radiation_flux.F90 b/radiation/radiation_flux.F90
index 8f9c92bf..df0ed946 100644
--- a/radiation/radiation_flux.F90
+++ b/radiation/radiation_flux.F90
@@ -20,14 +20,13 @@
 !   2021-01-20  R. Hogan  Added heating_rate_out_of_physical_bounds function
 !   2022-12-07  R. Hogan  Added top-of-atmosphere spectral output
 
-#include "ecrad_config.h"
-
 module radiation_flux
 
   use parkind1, only : jprb
 
   implicit none
   public
+#include "ecrad_config.h"
 
   !---------------------------------------------------------------------
   ! This derived type contains the output from the radiation
diff --git a/radiation/radiation_general_cloud_optics_data.F90 b/radiation/radiation_general_cloud_optics_data.F90
index f62e1a11..d432bf1a 100644
--- a/radiation/radiation_general_cloud_optics_data.F90
+++ b/radiation/radiation_general_cloud_optics_data.F90
@@ -14,8 +14,6 @@
 ! License: see the COPYING file for details
 !
 
-#include "ecrad_config.h"
-
 module radiation_general_cloud_optics_data
 
   use parkind1, only : jprb
@@ -24,6 +22,8 @@ module radiation_general_cloud_optics_data
 
   public
 
+#include "ecrad_config.h"
+
   !---------------------------------------------------------------------
   ! This type holds the configuration information to compute optical
   ! properties for a particular type of cloud or hydrometeor in one of
diff --git a/radiation/radiation_interface.F90 b/radiation/radiation_interface.F90
index 2b413ee3..b628ddcc 100644
--- a/radiation/radiation_interface.F90
+++ b/radiation/radiation_interface.F90
@@ -202,7 +202,10 @@ subroutine radiation(ncol, nlev, istartcol, iendcol, config, &
 
     use parkind1,                 only : jprb
     use yomhook,                  only : lhook, dr_hook, jphook
-
+#ifdef USE_TIMING
+    ! Timing library
+    use gptl,                  only: gptlstart, gptlstop
+#endif
     use radiation_io,             only : nulout
     use radiation_config,         only : config_type, &
          &   IGasModelMonochromatic, IGasModelIFSRRTMG, IGasModelECCKD, &
@@ -305,6 +308,11 @@ subroutine radiation(ncol, nlev, istartcol, iendcol, config, &
 
     real(jphook) :: hook_handle
 
+#ifdef USE_TIMING
+    integer :: ret
+    ret =  gptlstart('radiation_interface:radiation')
+#endif
+
     if (lhook) call dr_hook('radiation_interface:radiation',0,hook_handle)
 
     if (thermodynamics%pressure_hl(istartcol,2) &
@@ -324,7 +332,9 @@ subroutine radiation(ncol, nlev, istartcol, iendcol, config, &
       call single_level%get_albedos(istartcol, iendcol, config, &
            &                        sw_albedo_direct, sw_albedo_diffuse, &
            &                        lw_albedo)
-
+#ifdef USE_TIMING
+    ret =  gptlstart('gas_optics')
+#endif     
       ! Compute gas absorption optical depth in shortwave and
       ! longwave, shortwave single scattering albedo (i.e. fraction of
       ! extinction due to Rayleigh scattering), Planck functions and
@@ -353,7 +363,10 @@ subroutine radiation(ncol, nlev, istartcol, iendcol, config, &
                &  incoming_sw=incoming_sw)
         end if
       end if
-
+#ifdef USE_TIMING
+    ret =  gptlstop('gas_optics')
+    ret =  gptlstart('cloud_optics')
+#endif  
       if (config%do_clouds) then
         ! Crop the cloud fraction to remove clouds that have too small
         ! a fraction or water content; after this, we can safely
@@ -381,7 +394,10 @@ subroutine radiation(ncol, nlev, istartcol, iendcol, config, &
                &  od_sw_cloud, ssa_sw_cloud, g_sw_cloud)
         end if
       end if ! do_clouds
-
+#ifdef USE_TIMING
+    ret =  gptlstop('cloud_optics')
+    ret =  gptlstart('aerosol_optics')
+#endif  
       if (config%use_aerosols) then
         if (config%i_gas_model_sw == IGasModelMonochromatic) then
 !          call add_aerosol_optics_mono(nlev,istartcol,iendcol, &
@@ -399,7 +415,9 @@ subroutine radiation(ncol, nlev, istartcol, iendcol, config, &
           g_lw(:,:,istartcol:iendcol)   = 0.0_jprb
         end if
       end if
-
+#ifdef USE_TIMING
+    ret =  gptlstop('aerosol_optics')
+#endif  
       ! For diagnostic purposes, save these intermediate variables to
       ! a NetCDF file
       if (config%do_save_radiative_properties) then
@@ -418,7 +436,9 @@ subroutine radiation(ncol, nlev, istartcol, iendcol, config, &
              &  od_lw_cloud, ssa_lw_cloud, g_lw_cloud, &
              &  od_sw_cloud, ssa_sw_cloud, g_sw_cloud)
       end if
-
+#ifdef USE_TIMING
+    ret =  gptlstart('solver_longwave')
+#endif  
       if (config%do_lw) then
         if (config%iverbose >= 2) then
           write(nulout,'(a)') 'Computing longwave fluxes'
@@ -438,7 +458,7 @@ subroutine radiation(ncol, nlev, istartcol, iendcol, config, &
                &  planck_hl, lw_emission, lw_albedo, flux)
         else if (config%i_solver_lw == ISolverTripleclouds) then
           ! Compute fluxes using the Tripleclouds longwave solver
-          call solver_tripleclouds_lw(nlev,istartcol,iendcol, &
+          call solver_tripleclouds_lw(config%n_g_lw,nlev,istartcol,iendcol, &
                &  config, cloud, & 
                &  od_lw, ssa_lw, g_lw, od_lw_cloud, ssa_lw_cloud, g_lw_cloud, &
                &  planck_hl, lw_emission, lw_albedo, flux)
@@ -455,7 +475,10 @@ subroutine radiation(ncol, nlev, istartcol, iendcol, config, &
                &  planck_hl, lw_emission, lw_albedo, flux)
         end if
       end if
-
+#ifdef USE_TIMING
+    ret =  gptlstop('solver_longwave')
+    ret =  gptlstart('solver_shortwave')
+#endif  
       if (config%do_sw) then
         if (config%iverbose >= 2) then
           write(nulout,'(a)') 'Computing shortwave fluxes'
@@ -477,7 +500,7 @@ subroutine radiation(ncol, nlev, istartcol, iendcol, config, &
                &  incoming_sw, flux)
         else if (config%i_solver_sw == ISolverTripleclouds) then
           ! Compute fluxes using the Tripleclouds shortwave solver
-          call solver_tripleclouds_sw(nlev,istartcol,iendcol, &
+          call solver_tripleclouds_sw(config%n_g_sw,nlev,istartcol,iendcol, &
                &  config, single_level, cloud, & 
                &  od_sw, ssa_sw, g_sw, od_sw_cloud, ssa_sw_cloud, &
                &  g_sw_cloud, sw_albedo_direct, sw_albedo_diffuse, &
@@ -497,7 +520,9 @@ subroutine radiation(ncol, nlev, istartcol, iendcol, config, &
                &  incoming_sw, flux)
         end if
       end if
-
+#ifdef USE_TIMING
+    ret =  gptlstop('solver_shortwave')
+#endif  
       ! Store surface downwelling, and TOA, fluxes in bands from
       ! fluxes in g points
       call flux%calc_surface_spectral(config, istartcol, iendcol)
@@ -506,7 +531,9 @@ subroutine radiation(ncol, nlev, istartcol, iendcol, config, &
     end if
     
     if (lhook) call dr_hook('radiation_interface:radiation',1,hook_handle)
-
+#ifdef USE_TIMING
+    ret =  gptlstop('radiation_interface:radiation')
+#endif
   end subroutine radiation
 
 
diff --git a/radiation/radiation_mcica_lw.F90 b/radiation/radiation_mcica_lw.F90
index f6bfcf98..51c6e348 100644
--- a/radiation/radiation_mcica_lw.F90
+++ b/radiation/radiation_mcica_lw.F90
@@ -18,11 +18,11 @@
 !   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
 
   public
+#include "ecrad_config.h"
 
 contains
 
diff --git a/radiation/radiation_mcica_sw.F90 b/radiation/radiation_mcica_sw.F90
index 38b58290..27b5ba43 100644
--- a/radiation/radiation_mcica_sw.F90
+++ b/radiation/radiation_mcica_sw.F90
@@ -17,11 +17,11 @@
 !   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
 
   public
+#include "ecrad_config.h"
 
 contains
 
diff --git a/radiation/radiation_spectral_definition.F90 b/radiation/radiation_spectral_definition.F90
index 2135e8ad..7a489568 100644
--- a/radiation/radiation_spectral_definition.F90
+++ b/radiation/radiation_spectral_definition.F90
@@ -14,8 +14,6 @@
 ! License: see the COPYING file for details
 !
 
-#include "ecrad_config.h"
-
 module radiation_spectral_definition
 
   use parkind1,    only : jprb
@@ -24,6 +22,8 @@ module radiation_spectral_definition
 
   public
 
+#include "ecrad_config.h"
+
   real(jprb), parameter :: SolarReferenceTemperature       = 5777.0_jprb ! K
   real(jprb), parameter :: TerrestrialReferenceTemperature = 273.15_jprb ! K
 
diff --git a/radiation/radiation_tripleclouds_lw.F90 b/radiation/radiation_tripleclouds_lw.F90
index e7c0f26a..ad2f4630 100644
--- a/radiation/radiation_tripleclouds_lw.F90
+++ b/radiation/radiation_tripleclouds_lw.F90
@@ -22,20 +22,23 @@
 
 module radiation_tripleclouds_lw
 
-  public
+  private
+  public :: solver_tripleclouds_lw
+#include "ecrad_config.h"
 
 contains
   ! Small routine for scaling cloud optical depth in the cloudy
   ! regions
 #include "radiation_optical_depth_scaling.h"
 
+
   !---------------------------------------------------------------------
   ! This module contains just one subroutine, the longwave
   ! "Tripleclouds" solver in which cloud inhomogeneity is treated by
   ! dividing each model level into three regions, one clear and two
   ! cloudy (with differing optical depth). This approach was described
   ! by Shonk and Hogan (2008).
-  subroutine solver_tripleclouds_lw(nlev,istartcol,iendcol, &
+  subroutine solver_tripleclouds_lw(ng_lw_in, nlev,istartcol,iendcol, &
        &  config, cloud, & 
        &  od, ssa, g, od_cloud, ssa_cloud, g_cloud, planck_hl, &
        &  emission, albedo, &
@@ -59,6 +62,8 @@ subroutine solver_tripleclouds_lw(nlev,istartcol,iendcol, &
     implicit none
 
     ! Inputs
+    integer, intent(in) :: ng_lw_in           ! number of g-points: leading dimension 
+    ! set to the procedure argument ng_lw_in if ng_lw not set at compile time (see ecrad_config.h)
     integer, intent(in) :: nlev               ! number of model levels
     integer, intent(in) :: istartcol, iendcol ! range of columns to process
     type(config_type),        intent(in) :: config
@@ -66,7 +71,7 @@ subroutine solver_tripleclouds_lw(nlev,istartcol,iendcol, &
 
     ! Gas and aerosol optical depth of each layer at each longwave
     ! g-point
-    real(jprb), intent(in), dimension(config%n_g_lw,nlev,istartcol:iendcol) :: od
+    real(jprb), intent(in), dimension(ng_lw,nlev,istartcol:iendcol) :: od
 
     ! Gas and aerosol single-scattering albedo and asymmetry factor,
     ! only if longwave scattering by aerosols is to be represented
@@ -85,19 +90,19 @@ subroutine solver_tripleclouds_lw(nlev,istartcol,iendcol, &
 
     ! Planck function (emitted flux from a black body) at half levels
     ! and at the surface at each longwave g-point
-    real(jprb), intent(in), dimension(config%n_g_lw,nlev+1,istartcol:iendcol) :: planck_hl
+    real(jprb), intent(in), dimension(ng_lw,nlev+1,istartcol:iendcol) :: planck_hl
 
     ! Emission (Planck*emissivity) and albedo (1-emissivity) at the
     ! surface at each longwave g-point
-    real(jprb), intent(in), dimension(config%n_g_lw, istartcol:iendcol) :: emission, albedo
+    real(jprb), intent(in), dimension(ng_lw, istartcol:iendcol) :: emission, albedo
 
     ! Optical depth, single scattering albedo and asymmetry factor in
     ! each g-point including gas, aerosol and clouds
-    real(jprb), dimension(config%n_g_lw) :: od_total, ssa_total, g_total
+    real(jprb), dimension(ng_lw) :: od_total, ssa_total, g_total
 
     ! Modified optical depth after Tripleclouds scaling to represent
     ! cloud inhomogeneity
-    real(jprb), dimension(config%n_g_lw) :: od_cloud_new
+    real(jprb), dimension(ng_lw) :: od_cloud_new
 
     ! Output
     type(flux_type), intent(inout):: flux
@@ -122,48 +127,48 @@ subroutine solver_tripleclouds_lw(nlev,istartcol,iendcol, &
          &                istartcol:iendcol) :: u_matrix, v_matrix
 
     ! Diffuse reflection and transmission matrices of each layer
-    real(jprb), dimension(config%n_g_lw, nregions, nlev) :: reflectance, transmittance
+    real(jprb), dimension(ng_lw, nregions, nlev) :: reflectance, transmittance
 
     ! Emission by a layer into the upwelling or downwelling diffuse
     ! streams
-    real(jprb), dimension(config%n_g_lw, nregions, nlev) &
+    real(jprb), dimension(ng_lw, nregions, nlev) &
          &  :: source_up, source_dn
 
     ! Clear-sky reflectance and transmittance
-    real(jprb), dimension(config%n_g_lw, nlev) &
+    real(jprb), dimension(ng_lw, nlev) &
          &  :: ref_clear, trans_clear
 
     ! ...clear-sky equivalent
-    real(jprb), dimension(config%n_g_lw, nlev) &
+    real(jprb), dimension(ng_lw, nlev) &
          &  :: source_up_clear, source_dn_clear
 
     ! Total albedo of the atmosphere/surface just above a layer
     ! interface with respect to downwelling diffuse radiation at that
     ! interface, where level index = 1 corresponds to the
     ! top-of-atmosphere
-    real(jprb), dimension(config%n_g_lw, nregions, nlev+1) :: total_albedo
+    real(jprb), dimension(ng_lw, nregions, nlev+1) :: total_albedo
 
     ! Upwelling radiation just above a layer interface due to emission
     ! below that interface, where level index = 1 corresponds to the
     ! top-of-atmosphere
-    real(jprb), dimension(config%n_g_lw, nregions, nlev+1) :: total_source
+    real(jprb), dimension(ng_lw, nregions, nlev+1) :: total_source
 
     ! Total albedo and source of the atmosphere just below a layer interface
-    real(jprb), dimension(config%n_g_lw, nregions) &
+    real(jprb), dimension(ng_lw, nregions) &
          &  :: total_albedo_below, total_source_below
 
     ! Downwelling flux below and above an interface between
     ! layers into a plane perpendicular to the direction of the sun
-    real(jprb), dimension(config%n_g_lw, nregions) &
+    real(jprb), dimension(ng_lw, nregions) &
          &  :: flux_dn, flux_dn_below, flux_up
 
     ! ...clear-sky equivalent (no distinction between "above/below")
-    real(jprb), dimension(config%n_g_lw, nlev+1) &
+    real(jprb), dimension(ng_lw, nlev+1) &
          &  :: flux_dn_clear, flux_up_clear
 
     ! Clear-sky equivalent, but actually its reciprocal to replace
     ! some divisions by multiplications
-    real(jprb), dimension(config%n_g_lw, nregions) :: inv_denom
+    real(jprb), dimension(ng_lw, nregions) :: inv_denom
 
     ! Identify clear-sky layers, with pseudo layers for outer space
     ! and below the ground, both treated as single-region clear skies
@@ -175,7 +180,7 @@ subroutine solver_tripleclouds_lw(nlev,istartcol,iendcol, &
     ! Index of the highest cloudy layer
     integer :: i_cloud_top
 
-    integer :: jcol, jlev, jg, jreg, jreg2, ng
+    integer :: jcol, jlev, jg, jreg, jreg2
 
     real(jphook) :: hook_handle
 
@@ -184,8 +189,6 @@ subroutine solver_tripleclouds_lw(nlev,istartcol,iendcol, &
     ! --------------------------------------------------------
     ! Section 1: Prepare general variables and arrays
     ! --------------------------------------------------------
-    ! Copy array dimensions to local variables for convenience
-    ng   = config%n_g_lw
 
     ! Compute the wavelength-independent region fractions and
     ! optical-depth scalings
@@ -236,26 +239,26 @@ subroutine solver_tripleclouds_lw(nlev,istartcol,iendcol, &
         ! No scattering in clear-sky flux calculation; note that here
         ! the first two dimensions of the input arrays are unpacked
         ! into vectors inside the routine        
-        call calc_no_scattering_transmittance_lw(ng*nlev, od(:,:,jcol), &
+        call calc_no_scattering_transmittance_lw(ng_lw*nlev, od(:,:,jcol), &
              &  planck_hl(:,1:nlev,jcol), planck_hl(:,2:nlev+1, jcol), &
              &  trans_clear, source_up_clear, source_dn_clear)
         ! Ensure that clear-sky reflectance is zero since it may be
         ! used in cloudy-sky case
         ref_clear = 0.0_jprb
         ! Simple down-then-up method to compute fluxes
-        call calc_fluxes_no_scattering_lw(ng, nlev, &
+        call calc_fluxes_no_scattering_lw(ng_lw, nlev, &
              &  trans_clear, source_up_clear, source_dn_clear, &
              &  emission(:,jcol), albedo(:,jcol), &
              &  flux_up_clear, flux_dn_clear)
       else
         ! Scattering in clear-sky flux calculation
-        call calc_ref_trans_lw(ng*nlev, &
+        call calc_ref_trans_lw(ng_lw*nlev, &
              &  od(:,:,jcol), ssa(:,:,jcol), g(:,:,jcol), &
              &  planck_hl(:,1:nlev,jcol), planck_hl(:,2:nlev+1,jcol), &
              &  ref_clear, trans_clear, &
              &  source_up_clear, source_dn_clear)
         ! Use adding method to compute fluxes
-        call adding_ica_lw(ng, nlev, &
+        call adding_ica_lw(ng_lw, nlev, &
              &  ref_clear, trans_clear, source_up_clear, source_dn_clear, &
              &  emission(:,jcol), albedo(:,jcol), &
              &  flux_up_clear, flux_dn_clear)
@@ -267,7 +270,7 @@ subroutine solver_tripleclouds_lw(nlev,istartcol,iendcol, &
           sum_up = 0.0_jprb
           sum_dn = 0.0_jprb
           !$omp simd reduction(+:sum_up, sum_dn)
-          do jg = 1,ng
+          do jg = 1,ng_lw
             sum_up = sum_up + flux_up_clear(jg,jlev)
             sum_dn = sum_dn + flux_dn_clear(jg,jlev)
           end do
@@ -276,7 +279,7 @@ subroutine solver_tripleclouds_lw(nlev,istartcol,iendcol, &
         end do
 
         ! Store surface spectral downwelling fluxes / TOA upwelling
-        do jg = 1,ng
+        do jg = 1,ng_lw
           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
@@ -357,7 +360,7 @@ subroutine solver_tripleclouds_lw(nlev,istartcol,iendcol, &
                        &     *  od_cloud_new / (ssa_total*od_total)
                 end where
               end if
-              call calc_ref_trans_lw(ng, &
+              call calc_ref_trans_lw(ng_lw, &
                    &  od_total, ssa_total, g_total, &
                    &  planck_hl(:,jlev,jcol), planck_hl(:,jlev+1,jcol), &
                    &  reflectance(:,jreg,jlev), transmittance(:,jreg,jlev), &
@@ -365,7 +368,7 @@ subroutine solver_tripleclouds_lw(nlev,istartcol,iendcol, &
             else
               ! No-scattering case: use simpler functions for
               ! transmission and emission
-              call calc_no_scattering_transmittance_lw(ng, od_total, &
+              call calc_no_scattering_transmittance_lw(ng_lw, od_total, &
                    &  planck_hl(:,jlev,jcol), planck_hl(:,jlev+1, jcol), &
                    &  transmittance(:,jreg,jlev), source_up(:,jreg,jlev), source_dn(:,jreg,jlev))
               reflectance(:,jreg,jlev) = 0.0_jprb
@@ -390,7 +393,7 @@ subroutine solver_tripleclouds_lw(nlev,istartcol,iendcol, &
       ! Calculate the upwelling radiation emitted by the surface, and
       ! copy the surface albedo into total_albedo 
       do jreg = 1,nregions
-        do jg = 1,ng
+        do jg = 1,ng_lw
           ! region_fracs(jreg,nlev,jcol) is the fraction of each region in the
           ! lowest model level
           total_source(jg,jreg,nlev+1) = region_fracs(jreg,nlev,jcol)*emission(jg,jcol)
@@ -408,7 +411,7 @@ subroutine solver_tripleclouds_lw(nlev,istartcol,iendcol, &
         if (is_clear_sky_layer(jlev)) then
           total_albedo_below = 0.0_jprb
           total_source_below = 0.0_jprb
-          do jg = 1,ng
+          do jg = 1,ng_lw
             inv_denom(jg,1) = 1.0_jprb &
                  &  / (1.0_jprb - total_albedo(jg,1,jlev+1)*reflectance(jg,1,jlev))
             total_albedo_below(jg,1) = reflectance(jg,1,jlev) &
@@ -436,7 +439,7 @@ subroutine solver_tripleclouds_lw(nlev,istartcol,iendcol, &
           total_albedo(:,:,jlev) = total_albedo_below(:,:)
           total_source(:,:,jlev) = total_source_below(:,:)
         else
-          total_source(:,:,jlev) = singlemat_x_vec(ng,ng,nregions,&
+          total_source(:,:,jlev) = singlemat_x_vec(ng_lw,ng_lw,nregions,&
                &  u_matrix(:,:,jlev,jcol), total_source_below)
           ! Use overlap matrix and exclude "anomalous" horizontal
           ! transport described by Shonk & Hogan (2008).  Therefore,
@@ -469,7 +472,7 @@ subroutine solver_tripleclouds_lw(nlev,istartcol,iendcol, &
         else
           sum_dn = 0.0_jprb
           !$omp simd reduction(+:sum_dn)
-          do jg = 1,ng
+          do jg = 1,ng_lw
             sum_dn = sum_dn + flux_dn_clear(jg,jlev)
           end do
           flux%lw_dn(jcol,jlev) = sum_dn
@@ -492,7 +495,7 @@ subroutine solver_tripleclouds_lw(nlev,istartcol,iendcol, &
 
       sum_up = 0.0_jprb
       !$omp simd reduction(+:sum_up)
-      do jg = 1,ng
+      do jg = 1,ng_lw
         sum_up = sum_up + flux_up(jg,1)
       end do
       flux%lw_up(jcol,i_cloud_top) = sum_up
@@ -506,7 +509,7 @@ subroutine solver_tripleclouds_lw(nlev,istartcol,iendcol, &
         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
+        do jg = 1,ng_lw
           sum_up = sum_up + flux_up(jg,1)
         end do
         flux%lw_up(jcol,jlev) = sum_up
@@ -533,7 +536,7 @@ subroutine solver_tripleclouds_lw(nlev,istartcol,iendcol, &
       do jlev = i_cloud_top,nlev
 
         if (is_clear_sky_layer(jlev)) then
-          do jg = 1,ng
+          do jg = 1,ng_lw
             flux_dn(jg,1) = (transmittance(jg,1,jlev)*flux_dn(jg,1) &
                  &  + reflectance(jg,1,jlev)*total_source(jg,1,jlev+1) + source_dn(jg,1,jlev) ) &
                  &  / (1.0_jprb - reflectance(jg,1,jlev)*total_albedo(jg,1,jlev+1))
@@ -552,7 +555,7 @@ subroutine solver_tripleclouds_lw(nlev,istartcol,iendcol, &
              &    .and. is_clear_sky_layer(jlev+1))) then
           ! Account for overlap rules in translating fluxes just above
           ! a layer interface to the values just below
-          flux_dn_below = singlemat_x_vec(ng,ng,nregions, &
+          flux_dn_below = singlemat_x_vec(ng_lw,ng_lw,nregions, &
                &  v_matrix(:,:,jlev+1,jcol), flux_dn)
           flux_dn = flux_dn_below
         end if ! Otherwise the fluxes in each region are the same so
@@ -563,7 +566,7 @@ subroutine solver_tripleclouds_lw(nlev,istartcol,iendcol, &
         sum_dn = 0.0_jprb
         do jreg = 1,nregions
           !$omp simd reduction(+:sum_up, sum_dn)
-          do jg = 1,ng
+          do jg = 1,ng_lw
             sum_up = sum_up + flux_up(jg,jreg)
             sum_dn = sum_dn + flux_dn(jg,jreg)
           end do
@@ -594,7 +597,7 @@ subroutine solver_tripleclouds_lw(nlev,istartcol,iendcol, &
         ! fluxes into the regions of the lowest layer; we sum over
         ! regions first to provide a simple spectral flux upwelling
         ! from the surface
-        call calc_lw_derivatives_region(ng, nlev, nregions, jcol, transmittance, &
+        call calc_lw_derivatives_region(ng_lw, nlev, nregions, jcol, transmittance, &
              &  u_matrix(:,:,:,jcol), sum(flux_up,2), flux%lw_derivatives)
       end if
       
diff --git a/radiation/radiation_tripleclouds_sw.F90 b/radiation/radiation_tripleclouds_sw.F90
index c2c0fb0a..c7a743c5 100644
--- a/radiation/radiation_tripleclouds_sw.F90
+++ b/radiation/radiation_tripleclouds_sw.F90
@@ -24,6 +24,7 @@
 module radiation_tripleclouds_sw
 
   public
+#include "ecrad_config.h"
 
 contains
   ! Provides elemental function "delta_eddington"
@@ -39,7 +40,7 @@ module radiation_tripleclouds_sw
   ! dividing each model level into three regions, one clear and two
   ! cloudy (with differing optical depth). This approach was described
   ! by Shonk and Hogan (2008).
-  subroutine solver_tripleclouds_sw(nlev,istartcol,iendcol, &
+  subroutine solver_tripleclouds_sw(ng_sw_in,nlev,istartcol,iendcol, &
        &  config, single_level, cloud, & 
        &  od, ssa, g, od_cloud, ssa_cloud, g_cloud, &
        &  albedo_direct, albedo_diffuse, incoming_sw, &
@@ -65,6 +66,8 @@ subroutine solver_tripleclouds_sw(nlev,istartcol,iendcol, &
     integer, parameter :: nregions = 3
     
     ! Inputs
+    integer, intent(in) :: ng_sw_in               ! number of g-points: leading dimension 
+    ! set to the procedure argument ng_sw_in if ng_sw not set at compile time (see ecrad_config.h)
     integer, intent(in) :: nlev               ! number of model levels
     integer, intent(in) :: istartcol, iendcol ! range of columns to process
     type(config_type),        intent(in) :: config
@@ -73,7 +76,7 @@ subroutine solver_tripleclouds_sw(nlev,istartcol,iendcol, &
 
     ! 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) &
+    real(jprb), intent(in), dimension(ng_sw,nlev,istartcol:iendcol) &
          &  :: od, ssa, g
 
     ! Cloud and precipitation optical depth, single-scattering albedo and
@@ -84,13 +87,13 @@ subroutine solver_tripleclouds_sw(nlev,istartcol,iendcol, &
     ! Optical depth, single scattering albedo and asymmetry factor in
     ! each g-point of each cloudy region including gas, aerosol and
     ! clouds
-    real(jprb), dimension(config%n_g_sw,2:nregions) &
+    real(jprb), dimension(ng_sw,2:nregions) &
          &  :: od_total, ssa_total, g_total
 
     ! Direct and diffuse surface albedos, and the incoming shortwave
     ! 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) &
+    real(jprb), intent(in), dimension(ng_sw,istartcol:iendcol) &
          &  :: albedo_direct, albedo_diffuse, incoming_sw
 
     ! Output
@@ -111,20 +114,20 @@ subroutine solver_tripleclouds_sw(nlev,istartcol,iendcol, &
 
     ! Diffuse reflection and transmission matrices in the cloudy
     ! regions of each layer
-    real(jprb), dimension(config%n_g_sw, 2:nregions, nlev) &
+    real(jprb), dimension(ng_sw, 2:nregions, nlev) &
          &  :: reflectance, transmittance
 
     ! Terms translating the direct flux entering the layer from above
     ! to the reflected radiation exiting upwards (ref_dir) and the
     ! scattered radiation exiting downwards (trans_dir_diff), along with the
     ! direct unscattered transmission matrix (trans_dir_dir).
-    real(jprb), dimension(config%n_g_sw, 2:nregions, nlev) &
+    real(jprb), dimension(ng_sw, 2:nregions, nlev) &
          &  :: ref_dir, trans_dir_diff, trans_dir_dir
 
     ! As above but for the clear regions; clear and cloudy layers are
     ! separated out so that calc_ref_trans_sw can be called on the
     ! entire clear-sky atmosphere at once
-    real(jprb), dimension(config%n_g_sw, nlev) &
+    real(jprb), dimension(ng_sw, nlev) &
     	 &  :: reflectance_clear, transmittance_clear, &
          &     ref_dir_clear, trans_dir_diff_clear, trans_dir_dir_clear
 
@@ -132,30 +135,30 @@ subroutine solver_tripleclouds_sw(nlev,istartcol,iendcol, &
     ! interface with respect to downwelling diffuse and direct
     ! (respectively) radiation at that interface, where level index =
     ! 1 corresponds to the top-of-atmosphere
-    real(jprb), dimension(config%n_g_sw, nregions, nlev+1) &
+    real(jprb), dimension(ng_sw, nregions, nlev+1) &
          &  :: total_albedo, total_albedo_direct
 
     ! ...equivalent values for clear-skies
-    real(jprb), dimension(config%n_g_sw, nlev+1) &
+    real(jprb), dimension(ng_sw, nlev+1) &
          &  :: total_albedo_clear, total_albedo_clear_direct
 
     ! Total albedo of the atmosphere just below a layer interface
-    real(jprb), dimension(config%n_g_sw, nregions) &
+    real(jprb), dimension(ng_sw, nregions) &
          &  :: total_albedo_below, total_albedo_below_direct
 
     ! Direct downwelling flux below and above an interface between
     ! layers into a plane perpendicular to the direction of the sun
-    real(jprb), dimension(config%n_g_sw, nregions) :: direct_dn
+    real(jprb), dimension(ng_sw, nregions) :: direct_dn
     ! Diffuse equivalents
-    real(jprb), dimension(config%n_g_sw, nregions) :: flux_dn, flux_up
+    real(jprb), dimension(ng_sw, nregions) :: flux_dn, flux_up
 
     ! ...clear-sky equivalent (no distinction between "above/below")
-    real(jprb), dimension(config%n_g_sw) &
+    real(jprb), dimension(ng_sw) &
          &  :: direct_dn_clear, flux_dn_clear, flux_up_clear
 
     ! Clear-sky equivalent, but actually its reciprocal to replace
     ! some divisions by multiplications
-    real(jprb), dimension(config%n_g_sw, nregions) :: inv_denom
+    real(jprb), dimension(ng_sw, nregions) :: inv_denom
 
     ! Identify clear-sky layers, with pseudo layers for outer space
     ! and below the ground, both treated as single-region clear skies
@@ -170,7 +173,7 @@ subroutine solver_tripleclouds_sw(nlev,istartcol,iendcol, &
     ! Local cosine of solar zenith angle
     real(jprb) :: mu0
 
-    integer :: jcol, jlev, jg, jreg, iband, jreg2, ng
+    integer :: jcol, jlev, jg, jreg, iband, jreg2
 
     real(jphook) :: hook_handle
 
@@ -179,8 +182,6 @@ subroutine solver_tripleclouds_sw(nlev,istartcol,iendcol, &
     ! --------------------------------------------------------
     ! Section 1: Prepare general variables and arrays
     ! --------------------------------------------------------
-    ! Copy array dimensions to local variables for convenience
-    ng   = config%n_g_sw
 
     ! Compute the wavelength-independent region fractions and
     ! optical-depth scalings
@@ -266,7 +267,7 @@ subroutine solver_tripleclouds_sw(nlev,istartcol,iendcol, &
       ! are computed for each layer
 
       ! Clear-sky quantities for all layers at once
-      call calc_ref_trans_sw(ng*nlev, &
+      call calc_ref_trans_sw(ng_sw*nlev, &
           &  mu0, od(:,:,jcol), ssa(:,:,jcol), g(:,:,jcol), &
           &  reflectance_clear, transmittance_clear, &
           &  ref_dir_clear, trans_dir_diff_clear, &
@@ -276,23 +277,42 @@ subroutine solver_tripleclouds_sw(nlev,istartcol,iendcol, &
       do jlev = 1,nlev ! Start at top-of-atmosphere
         if (.not. is_clear_sky_layer(jlev)) then
           do jreg = 2,nregions
-            do jg = 1,ng
-              ! Mapping from g-point to band
-              iband = config%i_band_from_reordered_g_sw(jg)
-              scat_od = od(jg,jlev,jcol)*ssa(jg,jlev,jcol)
-              scat_od_cloud = od_cloud(iband,jlev,jcol) &
-                   &  * ssa_cloud(iband,jlev,jcol) * od_scaling(jreg,jlev,jcol)
-              ! Add scaled cloud optical depth to clear-sky value
-              od_total(jg,jreg) = od(jg,jlev,jcol) &
-                   &  + od_cloud(iband,jlev,jcol)*od_scaling(jreg,jlev,jcol)
-              ! Compute single-scattering albedo and asymmetry
-              ! factor of gas-cloud combination
-              ssa_total(jg,jreg) = (scat_od+scat_od_cloud) &
-                   &  / od_total(jg,jreg)
-              g_total(jg,jreg) = (scat_od*g(jg,jlev,jcol) &
-                   &         + scat_od_cloud * g_cloud(iband,jlev,jcol)) &
-                   &      / (scat_od + scat_od_cloud)
-            end do
+            if (ng_sw== config%n_bands_sw) then 
+              ! ecCKD, simpler loop indices (can vectorize better)
+              do jg = 1,ng_sw
+                scat_od = od(jg,jlev,jcol)*ssa(jg,jlev,jcol)
+                scat_od_cloud = od_cloud(jg,jlev,jcol) &
+                    &  * ssa_cloud(jg,jlev,jcol) * od_scaling(jreg,jlev,jcol)
+                ! Add scaled cloud optical depth to clear-sky value
+                od_total(jg,jreg) = od(jg,jlev,jcol) &
+                    &  + od_cloud(jg,jlev,jcol)*od_scaling(jreg,jlev,jcol)
+                ! Compute single-scattering albedo and asymmetry
+                ! factor of gas-cloud combination
+                ssa_total(jg,jreg) = (scat_od+scat_od_cloud) &
+                    &  / od_total(jg,jreg)
+                g_total(jg,jreg) = (scat_od*g(jg,jlev,jcol) &
+                    &         + scat_od_cloud * g_cloud(jg,jlev,jcol)) &
+                    &      / (scat_od + scat_od_cloud)
+              end do  
+            else 
+              do jg = 1,ng_sw
+                ! Mapping from g-point to band
+                iband = config%i_band_from_reordered_g_sw(jg)
+                scat_od = od(jg,jlev,jcol)*ssa(jg,jlev,jcol)
+                scat_od_cloud = od_cloud(iband,jlev,jcol) &
+                    &  * ssa_cloud(iband,jlev,jcol) * od_scaling(jreg,jlev,jcol)
+                ! Add scaled cloud optical depth to clear-sky value
+                od_total(jg,jreg) = od(jg,jlev,jcol) &
+                    &  + od_cloud(iband,jlev,jcol)*od_scaling(jreg,jlev,jcol)
+                ! Compute single-scattering albedo and asymmetry
+                ! factor of gas-cloud combination
+                ssa_total(jg,jreg) = (scat_od+scat_od_cloud) &
+                    &  / od_total(jg,jreg)
+                g_total(jg,jreg) = (scat_od*g(jg,jlev,jcol) &
+                    &         + scat_od_cloud * g_cloud(iband,jlev,jcol)) &
+                    &      / (scat_od + scat_od_cloud)
+              end do
+            end if 
           end do
 
           if (config%do_sw_delta_scaling_with_gases) then
@@ -301,7 +321,7 @@ subroutine solver_tripleclouds_sw(nlev,istartcol,iendcol, &
             call delta_eddington(od_total, ssa_total, g_total)
           end if
           ! Both cloudy regions at once
-          call calc_ref_trans_sw(ng*(nregions-1), &
+          call calc_ref_trans_sw(ng_sw*(nregions-1), &
                &  mu0, od_total, ssa_total, g_total, &
                &  reflectance(:,:,jlev), transmittance(:,:,jlev), &
                &  ref_dir(:,:,jlev), trans_dir_diff(:,:,jlev), &
@@ -317,7 +337,7 @@ subroutine solver_tripleclouds_sw(nlev,istartcol,iendcol, &
       total_albedo_direct(:,:,:) = 0.0_jprb
 
       ! Copy surface albedos in clear-sky region
-      do jg = 1,ng
+      do jg = 1,ng_sw
         total_albedo(jg,1,nlev+1) = albedo_diffuse(jg,jcol)
         total_albedo_direct(jg,1,nlev+1) &
              &  = mu0 * albedo_direct(jg,jcol)
@@ -348,7 +368,7 @@ subroutine solver_tripleclouds_sw(nlev,istartcol,iendcol, &
           ! For clear-skies there is no need to consider "above" and
           ! "below" quantities since with no cloud overlap to worry
           ! about, these are the same
-          do jg = 1,ng
+          do jg = 1,ng_sw
             inv_denom(jg,1) = 1.0_jprb &
                  &  / (1.0_jprb - total_albedo_clear(jg,jlev+1)*reflectance_clear(jg,jlev))
             total_albedo_clear(jg,jlev) = reflectance_clear(jg,jlev) &
@@ -363,7 +383,7 @@ subroutine solver_tripleclouds_sw(nlev,istartcol,iendcol, &
         end if
 
         ! All-sky fluxes: first the clear region
-        do jg = 1,ng
+        do jg = 1,ng_sw
           inv_denom(jg,1) = 1.0_jprb &
                &  / (1.0_jprb - total_albedo(jg,1,jlev+1)*reflectance_clear(jg,jlev))
           total_albedo_below(jg,1) = reflectance_clear(jg,jlev) &
@@ -378,7 +398,7 @@ subroutine solver_tripleclouds_sw(nlev,istartcol,iendcol, &
         ! Then the cloudy regions if any cloud is present in this layer
         if (.not. is_clear_sky_layer(jlev)) then
           do jreg = 2,nregions
-            do jg = 1,ng
+            do jg = 1,ng_sw
               inv_denom(jg,jreg) = 1.0_jprb / (1.0_jprb &
                    &  - total_albedo(jg,jreg,jlev+1)*reflectance(jg,jreg,jlev))
               total_albedo_below(jg,jreg) = reflectance(jg,jreg,jlev) &
@@ -455,7 +475,7 @@ subroutine solver_tripleclouds_sw(nlev,istartcol,iendcol, &
       sum_dn_dir = 0.0_jprb
       do jreg = 1,nregions
         !$omp simd reduction(+:sum_up, sum_dn_dir)
-        do jg = 1,ng
+        do jg = 1,ng_sw
           sum_up     = sum_up     + flux_up(jg,jreg)
           sum_dn_dir = sum_dn_dir + direct_dn(jg,jreg)
         end do
@@ -469,7 +489,7 @@ subroutine solver_tripleclouds_sw(nlev,istartcol,iendcol, &
         sum_up     = 0.0_jprb
         sum_dn_dir = 0.0_jprb
         !$omp simd reduction(+:sum_up, sum_dn_dir)
-        do jg = 1,ng
+        do jg = 1,ng_sw
           sum_up     = sum_up     + flux_up_clear(jg)
           sum_dn_dir = sum_dn_dir + direct_dn_clear(jg)
         end do
@@ -518,7 +538,7 @@ subroutine solver_tripleclouds_sw(nlev,istartcol,iendcol, &
       ! Final loop back down through the atmosphere to compute fluxes
       do jlev = 1,nlev
         if (config%do_clear) then
-          do jg = 1,ng
+          do jg = 1,ng_sw
             flux_dn_clear(jg) = (transmittance_clear(jg,jlev)*flux_dn_clear(jg) + direct_dn_clear(jg) &
                &  * (trans_dir_dir_clear(jg,jlev)*total_albedo_clear_direct(jg,jlev+1)*reflectance_clear(jg,jlev) &
                &     + trans_dir_diff_clear(jg,jlev) )) &
@@ -530,7 +550,7 @@ subroutine solver_tripleclouds_sw(nlev,istartcol,iendcol, &
         end if
 
         ! All-sky fluxes: first the clear region...
-        do jg = 1,ng
+        do jg = 1,ng_sw
           flux_dn(jg,1) = (transmittance_clear(jg,jlev)*flux_dn(jg,1) + direct_dn(jg,1) &
                &  * (trans_dir_dir_clear(jg,jlev)*total_albedo_direct(jg,1,jlev+1)*reflectance_clear(jg,jlev) &
                &     + trans_dir_diff_clear(jg,jlev) )) &
@@ -547,7 +567,7 @@ subroutine solver_tripleclouds_sw(nlev,istartcol,iendcol, &
           direct_dn(:,2:nregions)= 0.0_jprb
         else
           do jreg = 2,nregions
-            do jg = 1,ng
+            do jg = 1,ng_sw
               flux_dn(jg,jreg) = (transmittance(jg,jreg,jlev)*flux_dn(jg,jreg) + direct_dn(jg,jreg) &
                    &  * (trans_dir_dir(jg,jreg,jlev)*total_albedo_direct(jg,jreg,jlev+1)*reflectance(jg,jreg,jlev) &
                    &     + trans_dir_diff(jg,jreg,jlev) )) &
@@ -563,9 +583,9 @@ subroutine solver_tripleclouds_sw(nlev,istartcol,iendcol, &
              &    .and. is_clear_sky_layer(jlev+1))) then
           ! Account for overlap rules in translating fluxes just above
           ! a layer interface to the values just below
-          flux_dn = singlemat_x_vec(ng,ng,nregions, &
+          flux_dn = singlemat_x_vec(ng_sw,ng_sw,nregions, &
                &  v_matrix(:,:,jlev+1,jcol), flux_dn)
-          direct_dn = singlemat_x_vec(ng,ng,nregions, &
+          direct_dn = singlemat_x_vec(ng_sw,ng_sw,nregions, &
                &  v_matrix(:,:,jlev+1,jcol), direct_dn)
         end if ! Otherwise the fluxes in each region are the same so
                ! nothing to do
@@ -579,7 +599,7 @@ subroutine solver_tripleclouds_sw(nlev,istartcol,iendcol, &
         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
+          do jg = 1,ng_sw
             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)
@@ -595,7 +615,7 @@ subroutine solver_tripleclouds_sw(nlev,istartcol,iendcol, &
           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
+          do jg = 1,ng_sw
             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)
diff --git a/radiation/radiation_two_stream.F90 b/radiation/radiation_two_stream.F90
index 794b92cd..1ca85638 100644
--- a/radiation/radiation_two_stream.F90
+++ b/radiation/radiation_two_stream.F90
@@ -21,14 +21,13 @@
 !   2022-11-22  P Ukkonen/R Hogan  Single precision uses no double precision
 !   2023-09-28  R Hogan  Increased security for single-precision SW "k"
 
-#include "ecrad_config.h"
-
 module radiation_two_stream
 
   use parkind1, only : jprb, jprd
 
   implicit none
   public
+#include "ecrad_config.h"
 
   ! Elsasser's factor: the effective factor by which the zenith
   ! optical depth needs to be multiplied to account for longwave
diff --git a/test/ifs/Makefile b/test/ifs/Makefile
index 485a2ab8..070e0b55 100644
--- a/test/ifs/Makefile
+++ b/test/ifs/Makefile
@@ -127,6 +127,10 @@ test_ecckd_noaer:
 		use_aerosols=false
 	$(DRIVER) config_ecckd_noaer.nam $(INPUT).nc $(INPUT)_ecckd_noaer_out.nc
 
+test_ifsdriver_blocked_ecckd_tc:
+	$(CHANGENAM) $(CONFIG_ECCKD) config_net.nam do_save_net_fluxes=true do_write_double_precision=true do_save_spectral_flux=false nblocksize=8 nrepeat=10
+	$(IFS_DRIVER_BLOCKED) config_net.nam $(INPUT).nc $(INPUT)_ifsdriver_blocked_ecckd_tc_out.nc | tee $(INPUT)_ifsdriver_blocked_ecckd_tc_out.log
+
 # Test the different ways that aerosol optical properties can be
 # averaged, outputing the aerosol properties in each gas-optics
 # spectral interval, producing the following: