diff --git a/docs/src/optics/LookUpTables.md b/docs/src/optics/LookUpTables.md
index 02669cf58..1f9cad9a5 100644
--- a/docs/src/optics/LookUpTables.md
+++ b/docs/src/optics/LookUpTables.md
@@ -8,5 +8,6 @@ CurrentModule = RRTMGP.LookUpTables
 LookUpLW
 LookUpSW
 LookUpCld
+PadeCld
 AbstractLookUp
 ```
diff --git a/src/optics/CloudOptics.jl b/src/optics/CloudOptics.jl
index 338b35c79..f21bf0a9e 100644
--- a/src/optics/CloudOptics.jl
+++ b/src/optics/CloudOptics.jl
@@ -65,8 +65,7 @@ end
 This function computed the TwoSteam cloud optics properties using either the 
 lookup table method or pade method.
 """
-function compute_cld_props(lkp_cld, as, cld_mask, glay, gcol, ibnd, igpt)
-    use_lut = lkp_cld.use_lut
+function compute_cld_props(lkp_cld::LookUpCld, as, cld_mask, glay, gcol, ibnd, igpt)
     re_liq = as.cld_r_eff_liq[glay, gcol]
     re_ice = as.cld_r_eff_ice[glay, gcol]
     ice_rgh = as.ice_rgh
@@ -75,83 +74,99 @@ function compute_cld_props(lkp_cld, as, cld_mask, glay, gcol, ibnd, igpt)
     FT = eltype(re_liq)
     τl, τl_ssa, τl_ssag = FT(0), FT(0), FT(0)
     τi, τi_ssa, τi_ssag = FT(0), FT(0), FT(0)
-    if use_lut # use LUT interpolation
-        (;
-            lut_extliq,
-            lut_ssaliq,
-            lut_asyliq,
-            lut_extice,
-            lut_ssaice,
-            lut_asyice,
-            radliq_lwr,
-            radliq_upr,
-            radice_lwr,
-            radice_upr,
-            nsize_liq,
-            nsize_ice,
-        ) = lkp_cld
-        Δr_liq = (radliq_upr - radliq_lwr) / FT(nsize_liq - 1)
-        Δr_ice = (radice_upr - radice_lwr) / FT(nsize_ice - 1)
-        # cloud liquid particles
-        if cld_path_liq > eps(FT)
-            loc = Int(max(min(unsafe_trunc(Int, (re_liq - radliq_lwr) / Δr_liq) + 1, nsize_liq - 1), 1))
-            fac = (re_liq - radliq_lwr - (loc - 1) * Δr_liq) / Δr_liq
-            fc1 = FT(1) - fac
-            τl = (fc1 * lut_extliq[loc, ibnd] + fac * lut_extliq[loc + 1, ibnd]) * cld_path_liq
-            τl_ssa = (fc1 * lut_ssaliq[loc, ibnd] + fac * lut_ssaliq[loc + 1, ibnd]) * τl
-            τl_ssag = (fc1 * lut_asyliq[loc, ibnd] + fac * lut_asyliq[loc + 1, ibnd]) * τl_ssa
-        end
-        # cloud ice particles
-        if cld_path_ice > eps(FT)
-            loc = Int(max(min(unsafe_trunc(Int, (re_ice - radice_lwr) / Δr_ice) + 1, nsize_ice - 1), 1))
-            fac = (re_ice - radice_lwr - (loc - 1) * Δr_ice) / Δr_ice
-            fc1 = FT(1) - fac
-            τi = (fc1 * lut_extice[loc, ibnd, ice_rgh] + fac * lut_extice[loc + 1, ibnd, ice_rgh]) * cld_path_ice
-            τi_ssa = (fc1 * lut_ssaice[loc, ibnd, ice_rgh] + fac * lut_ssaice[loc + 1, ibnd, ice_rgh]) * τi
-            τi_ssag = (fc1 * lut_asyice[loc, ibnd, ice_rgh] + fac * lut_asyice[loc + 1, ibnd, ice_rgh]) * τi_ssa
-        end
-    else # use pade interpolation
-        (;
-            pade_extliq,
-            pade_ssaliq,
-            pade_asyliq,
-            pade_extice,
-            pade_ssaice,
-            pade_asyice,
-            pade_sizreg_extliq,
-            pade_sizreg_ssaliq,
-            pade_sizreg_asyliq,
-            pade_sizreg_extice,
-            pade_sizreg_ssaice,
-            pade_sizreg_asyice,
-        ) = lkp_cld
-        m_ext, m_ssa_g = 3, 3
-        n_ext, n_ssa_g = 3, 2
-        # Finds index into size regime table
-        # This works only if there are precisely three size regimes (four bounds) and it's
-        # previously guaranteed that size_bounds(1) <= size <= size_bounds(4)
-        if cld_path_liq > eps(FT)
-            irad = Int(min(floor((re_liq - pade_sizreg_extliq[2]) / pade_sizreg_extliq[3]) + 2, 3))
-            τl = pade_eval(ibnd, re_liq, irad, m_ext, n_ext, pade_extliq) * cld_path_liq
+    (;
+        lut_extliq,
+        lut_ssaliq,
+        lut_asyliq,
+        lut_extice,
+        lut_ssaice,
+        lut_asyice,
+        radliq_lwr,
+        radliq_upr,
+        radice_lwr,
+        radice_upr,
+    ) = lkp_cld
+    nsize_liq = LookUpTables.get_nsize_liq(lkp_cld)
+    nsize_ice = LookUpTables.get_nsize_ice(lkp_cld)
+    Δr_liq = (radliq_upr - radliq_lwr) / FT(nsize_liq - 1)
+    Δr_ice = (radice_upr - radice_lwr) / FT(nsize_ice - 1)
+    # cloud liquid particles
+    if cld_path_liq > eps(FT)
+        loc = Int(max(min(unsafe_trunc(Int, (re_liq - radliq_lwr) / Δr_liq) + 1, nsize_liq - 1), 1))
+        fac = (re_liq - radliq_lwr - (loc - 1) * Δr_liq) / Δr_liq
+        fc1 = FT(1) - fac
+        τl = (fc1 * lut_extliq[loc, ibnd] + fac * lut_extliq[loc + 1, ibnd]) * cld_path_liq
+        τl_ssa = (fc1 * lut_ssaliq[loc, ibnd] + fac * lut_ssaliq[loc + 1, ibnd]) * τl
+        τl_ssag = (fc1 * lut_asyliq[loc, ibnd] + fac * lut_asyliq[loc + 1, ibnd]) * τl_ssa
+    end
+    # cloud ice particles
+    if cld_path_ice > eps(FT)
+        loc = Int(max(min(unsafe_trunc(Int, (re_ice - radice_lwr) / Δr_ice) + 1, nsize_ice - 1), 1))
+        fac = (re_ice - radice_lwr - (loc - 1) * Δr_ice) / Δr_ice
+        fc1 = FT(1) - fac
+        τi = (fc1 * lut_extice[loc, ibnd, ice_rgh] + fac * lut_extice[loc + 1, ibnd, ice_rgh]) * cld_path_ice
+        τi_ssa = (fc1 * lut_ssaice[loc, ibnd, ice_rgh] + fac * lut_ssaice[loc + 1, ibnd, ice_rgh]) * τi
+        τi_ssag = (fc1 * lut_asyice[loc, ibnd, ice_rgh] + fac * lut_asyice[loc + 1, ibnd, ice_rgh]) * τi_ssa
+    end
 
-            irad = Int(min(floor((re_liq - pade_sizreg_ssaliq[2]) / pade_sizreg_ssaliq[3]) + 2, 3))
-            τl_ssa = (FT(1) - max(FT(0), pade_eval(ibnd, re_liq, irad, m_ssa_g, n_ssa_g, pade_ssaliq))) * τl
+    τ = τl + τi
+    τ_ssa = τl_ssa + τi_ssa
+    τ_ssag = (τl_ssag + τi_ssag) / max(eps(FT), τ_ssa)
+    τ_ssa /= max(eps(FT), τ)
 
-            irad = Int(min(floor((re_liq - pade_sizreg_asyliq[2]) / pade_sizreg_asyliq[3]) + 2, 3))
-            τl_ssag = pade_eval(ibnd, re_liq, irad, m_ssa_g, n_ssa_g, pade_asyliq) * τl_ssa
-        end
+    return (τ, τ_ssa, τ_ssag)
+end
+
+function compute_cld_props(lkp_cld::PadeCld, as, cld_mask, glay, gcol, ibnd, igpt)
+    re_liq = as.cld_r_eff_liq[glay, gcol]
+    re_ice = as.cld_r_eff_ice[glay, gcol]
+    ice_rgh = as.ice_rgh
+    cld_path_liq = as.cld_path_liq[glay, gcol]
+    cld_path_ice = as.cld_path_ice[glay, gcol]
+    FT = eltype(re_liq)
+    τl, τl_ssa, τl_ssag = FT(0), FT(0), FT(0)
+    τi, τi_ssa, τi_ssag = FT(0), FT(0), FT(0)
 
-        if cld_path_ice > eps(FT)
-            irad = Int(min(floor((re_ice - pade_sizreg_extice[2]) / pade_sizreg_extice[3]) + 2, 3))
+    (;
+        pade_extliq,
+        pade_ssaliq,
+        pade_asyliq,
+        pade_extice,
+        pade_ssaice,
+        pade_asyice,
+        pade_sizreg_extliq,
+        pade_sizreg_ssaliq,
+        pade_sizreg_asyliq,
+        pade_sizreg_extice,
+        pade_sizreg_ssaice,
+        pade_sizreg_asyice,
+    ) = lkp_cld
+    m_ext, m_ssa_g = 3, 3
+    n_ext, n_ssa_g = 3, 2
+    # Finds index into size regime table
+    # This works only if there are precisely three size regimes (four bounds) and it's
+    # previously guaranteed that size_bounds(1) <= size <= size_bounds(4)
+    if cld_path_liq > eps(FT)
+        irad = Int(min(floor((re_liq - pade_sizreg_extliq[2]) / pade_sizreg_extliq[3]) + 2, 3))
+        τl = pade_eval(ibnd, re_liq, irad, m_ext, n_ext, pade_extliq) * cld_path_liq
 
-            τi = pade_eval(ibnd, re_ice, irad, m_ext, n_ext, pade_extice, ice_rgh) * cld_path_ice
+        irad = Int(min(floor((re_liq - pade_sizreg_ssaliq[2]) / pade_sizreg_ssaliq[3]) + 2, 3))
+        τl_ssa = (FT(1) - max(FT(0), pade_eval(ibnd, re_liq, irad, m_ssa_g, n_ssa_g, pade_ssaliq))) * τl
 
-            irad = Int(min(floor((re_ice - pade_sizreg_ssaice[2]) / pade_sizreg_ssaice[3]) + 2, 3))
-            τi_ssa = (FT(1) - max(FT(0), pade_eval(ibnd, re_ice, irad, m_ssa_g, n_ssa_g, pade_ssaice, ice_rgh))) * τi
+        irad = Int(min(floor((re_liq - pade_sizreg_asyliq[2]) / pade_sizreg_asyliq[3]) + 2, 3))
+        τl_ssag = pade_eval(ibnd, re_liq, irad, m_ssa_g, n_ssa_g, pade_asyliq) * τl_ssa
+    end
 
-            irad = Int(min(floor((re_ice - pade_sizreg_asyice[2]) / pade_sizreg_asyice[3]) + 2, 3))
-            τi_ssag = pade_eval(ibnd, re_ice, irad, m_ssa_g, n_ssa_g, pade_asyice, ice_rgh) * τi_ssa
-        end
+    if cld_path_ice > eps(FT)
+        irad = Int(min(floor((re_ice - pade_sizreg_extice[2]) / pade_sizreg_extice[3]) + 2, 3))
+
+        τi = pade_eval(ibnd, re_ice, irad, m_ext, n_ext, pade_extice, ice_rgh) * cld_path_ice
+
+        irad = Int(min(floor((re_ice - pade_sizreg_ssaice[2]) / pade_sizreg_ssaice[3]) + 2, 3))
+        τi_ssa = (FT(1) - max(FT(0), pade_eval(ibnd, re_ice, irad, m_ssa_g, n_ssa_g, pade_ssaice, ice_rgh))) * τi
+
+        irad = Int(min(floor((re_ice - pade_sizreg_asyice[2]) / pade_sizreg_asyice[3]) + 2, 3))
+        τi_ssag = pade_eval(ibnd, re_ice, irad, m_ssa_g, n_ssa_g, pade_asyice, ice_rgh) * τi_ssa
     end
 
     τ = τl + τi
@@ -162,6 +177,7 @@ function compute_cld_props(lkp_cld, as, cld_mask, glay, gcol, ibnd, igpt)
     return (τ, τ_ssa, τ_ssag)
 end
 
+
 """
     pade_eval(
         ibnd,
diff --git a/src/optics/GasOptics.jl b/src/optics/GasOptics.jl
index 1dcaa7caf..289f7e657 100644
--- a/src/optics/GasOptics.jl
+++ b/src/optics/GasOptics.jl
@@ -39,7 +39,7 @@ function compute_col_gas_kernel!(
 end
 """
     compute_interp_fractions(
-        lkp::AbstractLookUp{FT},
+        lkp::AbstractLookUp,
         vmr,
         p_lay,
         t_lay,
@@ -50,16 +50,7 @@ end
 
 compute interpolation fractions for binary species parameter, pressure and temperature.
 """
-@inline function compute_interp_fractions(
-    lkp::AbstractLookUp{FT},
-    vmr,
-    p_lay,
-    t_lay,
-    tropo,
-    ibnd,
-    glay,
-    gcol,
-) where {FT <: AbstractFloat}
+@inline function compute_interp_fractions(lkp::AbstractLookUp, vmr, p_lay, t_lay, tropo, ibnd, glay, gcol)
     jftemp = compute_interp_frac_temp(lkp, t_lay, glay, gcol)
     jfpress = compute_interp_frac_press(lkp, p_lay, tropo, glay, gcol)
     jfη, col_mix = compute_interp_frac_η(lkp, vmr, tropo, jftemp[1], ibnd, glay, gcol)
@@ -68,7 +59,7 @@ end
 
 """
     compute_interp_frac_temp(
-        lkp::AbstractLookUp{FT},
+        lkp::AbstractLookUp,
         t_lay,
         glay,
         gcol,
@@ -76,7 +67,7 @@ end
 
 compute interpolation fraction for temperature.
 """
-@inline function compute_interp_frac_temp(lkp::AbstractLookUp{FT}, t_lay, glay, gcol) where {FT <: AbstractFloat}
+@inline function compute_interp_frac_temp(lkp::AbstractLookUp, t_lay, glay, gcol)
     (; Δ_t_ref, n_t_ref, t_ref) = lkp
 
     @inbounds jtemp = loc_lower(t_lay, Δ_t_ref, n_t_ref, t_ref)
@@ -110,32 +101,24 @@ end
 
 """
     compute_interp_frac_η(
-        lkp::AbstractLookUp{FT},
+        lkp::AbstractLookUp,
         vmr,
         tropo,
         jtemp,
         ibnd,
         glay,
         gcol,
-    ) where {FT<:AbstractFloat}
+    )
 
 Compute interpolation fraction for binary species parameter.
 """
-@inline function compute_interp_frac_η(
-    lkp::AbstractLookUp{FT},
-    vmr,
-    tropo,
-    jtemp,
-    ibnd,
-    glay,
-    gcol,
-) where {FT <: AbstractFloat}
+@inline function compute_interp_frac_η(lkp::AbstractLookUp, vmr, tropo, jtemp, ibnd, glay, gcol)
     (; n_η, key_species, vmr_ref) = lkp
     ig = view(key_species, :, tropo, ibnd)
 
     vmr1 = get_vmr(vmr, ig[1], glay, gcol)
     vmr2 = get_vmr(vmr, ig[2], glay, gcol)
-
+    FT = eltype(vmr1)
     itemp = 1
     @inbounds η_half = vmr_ref[tropo, ig[1] + 1, jtemp + itemp - 1] / vmr_ref[tropo, ig[2] + 1, jtemp + itemp - 1]
     col_mix1 = vmr1 + η_half * vmr2
@@ -161,7 +144,7 @@ end
 
 """
     compute_τ_ssa_lw_src!(
-        lkp::AbstractLookUp{FT},
+        lkp::AbstractLookUp,
         vmr,
         col_dry,
         igpt,
@@ -176,7 +159,7 @@ Compute optical thickness, single scattering albedo, asymmetry parameter
 and longwave sources whenever applicable.
 """
 @inline function compute_τ_ssa_lw_src!(
-    lkp::AbstractLookUp{FT},
+    lkp::AbstractLookUp,
     vmr,
     col_dry,
     igpt,
diff --git a/src/optics/LookUpTables.jl b/src/optics/LookUpTables.jl
index 5af652e54..5bab1cdcf 100644
--- a/src/optics/LookUpTables.jl
+++ b/src/optics/LookUpTables.jl
@@ -3,18 +3,18 @@ module LookUpTables
 using DocStringExtensions
 using Adapt
 
-export AbstractLookUp, LookUpLW, LookUpSW, LookUpCld
+export AbstractLookUp, LookUpLW, LookUpSW, LookUpCld, PadeCld
 
 """
-    AbstractLookUp{FT}
+    AbstractLookUp
 
 Abstract lookup table for longwave and shortwave problems.
 """
-abstract type AbstractLookUp{FT} end
+abstract type AbstractLookUp end
 
 """
     LookUpLW{FT,UI8A1D,IA1D,IA2D,IA3D,FTA1D,FTA2D,FTA3D,FTA4D} <: 
-        AbstractLookUp{FT}
+        AbstractLookUp
 
 Longwave lookup tables, used to compute optical properties. 
 
@@ -31,7 +31,7 @@ struct LookUpLW{
     FTA2D <: AbstractArray{FT, 2},
     FTA3D <: AbstractArray{FT, 3},
     FTA4D <: AbstractArray{FT, 4},
-} <: AbstractLookUp{FT}
+} <: AbstractLookUp
     "number of gases used in the lookup table"
     n_gases::Int
     "number of longwave bands"
@@ -401,7 +401,7 @@ end
 
 """
     LookUpSW{FT,UI8A1D,IA1D,IA2D,IA3D,FTA1D,FTA2D,FTA3D,FTA4D} <:
-        AbstractLookUp{FT}
+        AbstractLookUp
 
 Shortwave lookup tables, used to compute optical properties. 
 
@@ -418,7 +418,7 @@ struct LookUpSW{
     FTA2D <: AbstractArray{FT, 2},
     FTA3D <: AbstractArray{FT, 3},
     FTA4D <: AbstractArray{FT, 4},
-} <: AbstractLookUp{FT}
+} <: AbstractLookUp
     "number of gases used in the lookup table"
     n_gases::Int
     "number of shortwave bands"
@@ -785,11 +785,11 @@ function LookUpSW(ds, ::Type{FT}, ::Type{DA}) where {FT <: AbstractFloat, DA}
 end
 
 """
-    LookUpCld{FT,FTA1D,FTA2D,FTA3D,FTA4D}
+    LookUpCld{FT,FTA2D,FTA3D}
 
 Lookup table for cloud optics.
 
-This struct stores the tables and Pade coefficients for determing extinction coeffient, 
+This struct stores the lookup tables for determing extinction coeffient, 
 single-scattering albedo, and asymmetry parameter g as a function of effective radius.
 We compute the optical depth tau (=exintinction coeff * condensed water path)
 and the products tau*ssa and tau*ssa*g for liquid and ice cloud separately.
@@ -798,25 +798,7 @@ These are used to determine the optical properties of ice and water cloud togeth
 # Fields
 $(DocStringExtensions.FIELDS)
 """
-struct LookUpCld{B, FT, FTA1D, FTA2D, FTA3D, FTA4D} <: AbstractLookUp{FT}
-    "number of bands"
-    nband::Int
-    "number of ice roughness types"
-    nrghice::Int
-    "number of liquid particle sizes"
-    nsize_liq::Int
-    "number of ice particle sizes"
-    nsize_ice::Int
-    "number of size regimes"
-    nsizereg::Int
-    "number of extinction coefficients for pade approximation"
-    ncoeff_ext::Int
-    "number of ssa/g coefficients for pade approximation"
-    ncoeff_ssa_g::Int
-    "number of size regime boundaries for pade interpolation"
-    nbound::Int
-    "pair = 2"
-    pair::Int
+struct LookUpCld{FT, FTA2D, FTA3D} <: AbstractLookUp
     "liquid particle size lower bound for LUT interpolation"
     radliq_lwr::FT
     "liquid particle size upper bound for LUT interpolation"
@@ -841,6 +823,82 @@ struct LookUpCld{B, FT, FTA1D, FTA2D, FTA3D, FTA4D} <: AbstractLookUp{FT}
     lut_ssaice::FTA3D
     "LUT ice asymmetry parameter (`nsize_ice, nband, nrghice`)"
     lut_asyice::FTA3D
+    "beginning and ending wavenumber for each band (`2, nband`) cm⁻¹"
+    bnd_lims_wn::FTA2D
+end
+Adapt.@adapt_structure LookUpCld
+
+# number of bands (cloud lookup table dimension)
+get_nband(lkp::LookUpCld) = size(lkp.lut_extliq, 2)
+# number of ice roughness types (cloud lookup table dimension)
+get_nrghice(lkp::LookUpCld) = size(lkp.lut_extice, 3)
+# number of liquid particle sizes (cloud lookup table dimension)
+get_nsize_liq(lkp::LookUpCld) = size(lkp.lut_extliq, 1)
+# number of ice particle sizes (cloud lookup table dimension)
+get_nsize_ice(lkp::LookUpCld) = size(lkp.lut_extice, 1)
+# pair = 2 (cloud lookup table dimension)
+get_pair(::LookUpCld) = 2
+
+function LookUpCld(ds, ::Type{FT}, ::Type{DA}) where {FT <: AbstractFloat, DA}
+    FTA2D = DA{FT, 2}
+    FTA3D = DA{FT, 3}
+
+    nband = Int(ds.dim["nband"])
+    nrghice = Int(ds.dim["nrghice"])
+    nsize_liq = Int(ds.dim["nsize_liq"])
+    nsize_ice = Int(ds.dim["nsize_ice"])
+    pair = Int(ds.dim["pair"])
+
+    radliq_lwr = FT(Array(ds["radliq_lwr"]))
+    radliq_upr = FT(Array(ds["radliq_upr"]))
+    radliq_fac = FT(Array(ds["radliq_fac"]))
+
+    radice_lwr = FT(Array(ds["radice_lwr"]))
+    radice_upr = FT(Array(ds["radice_upr"]))
+    radice_fac = FT(Array(ds["radice_fac"]))
+
+    lut_extliq = FTA2D(Array(ds["lut_extliq"]))
+    lut_ssaliq = FTA2D(Array(ds["lut_ssaliq"]))
+    lut_asyliq = FTA2D(Array(ds["lut_asyliq"]))
+
+    lut_extice = FTA3D(Array(ds["lut_extice"]))
+    lut_ssaice = FTA3D(Array(ds["lut_ssaice"]))
+    lut_asyice = FTA3D(Array(ds["lut_asyice"]))
+
+    bnd_lims_wn = FTA2D(Array(ds["bnd_limits_wavenumber"]))
+
+    return (LookUpCld{FT, FTA2D, FTA3D}(
+        radliq_lwr,
+        radliq_upr,
+        radliq_fac,
+        radice_lwr,
+        radice_upr,
+        radice_fac,
+        lut_extliq,
+        lut_ssaliq,
+        lut_asyliq,
+        lut_extice,
+        lut_ssaice,
+        lut_asyice,
+        bnd_lims_wn,
+    ))
+end
+
+"""
+    PadeCld{FTA1D,FTA2D,FTA3D,FTA4D}
+
+Pade coefficients for cloud optics.
+
+This struct stores the Pade coefficients for determing extinction coeffient, 
+single-scattering albedo, and asymmetry parameter g as a function of effective radius.
+We compute the optical depth tau (=exintinction coeff * condensed water path)
+and the products tau*ssa and tau*ssa*g for liquid and ice cloud separately.
+These are used to determine the optical properties of ice and water cloud together.
+
+# Fields
+$(DocStringExtensions.FIELDS)
+"""
+struct PadeCld{FTA1D, FTA2D, FTA3D, FTA4D} <: AbstractLookUp
     "pade coefficients for liquid extinction (`nband, nsizereg, ncoeff_ext`)"
     pade_extliq::FTA3D
     "pade coefficients for liquid single scattening albedo (`nband, nsizereg, ncoeff_ssa_g`)"
@@ -867,46 +925,29 @@ struct LookUpCld{B, FT, FTA1D, FTA2D, FTA3D, FTA4D} <: AbstractLookUp{FT}
     pade_sizreg_asyice::FTA1D
     "beginning and ending wavenumber for each band (`2, nband`) cm⁻¹"
     bnd_lims_wn::FTA2D
-    "use LUT (default) or pade interpolation"
-    use_lut::B
 end
-Adapt.@adapt_structure LookUpCld
-
-function LookUpCld(ds, ::Type{FT}, ::Type{DA}, use_lut::Bool = true) where {FT <: AbstractFloat, DA}
-
+Adapt.@adapt_structure PadeCld
+
+# number of bands for Pade approximation
+get_nband(lkp::PadeCld) = size(lkp.pade_extliq, 1)
+# number of size regimes for Pade approximation
+get_nsizereg(lkp::PadeCld) = size(lkp.pade_extliq, 2)
+# number of extinction coefficients for pade approximation
+get_ncoeff_ext(lkp::PadeCld) = size(lkp.pade_extliq, 3)
+# number of ssa/g coefficients for pade approximation
+get_ncoeff_ssa_g(lkp::PadeCld) = size(lkp.pade_ssaliq, 3)
+# number of size regime boundaries for pade interpolation
+get_nbound(lkp::PadeCld) = size(lkp.pade_sizreg_extliq, 1)
+
+function PadeCld(ds, ::Type{FT}, ::Type{DA}) where {FT <: AbstractFloat, DA}
     FTA1D = DA{FT, 1}
     FTA2D = DA{FT, 2}
     FTA3D = DA{FT, 3}
     FTA4D = DA{FT, 4}
 
-    nband = Int(ds.dim["nband"])
-    nrghice = Int(ds.dim["nrghice"])
-    nsize_liq = Int(ds.dim["nsize_liq"])
-    nsize_ice = Int(ds.dim["nsize_ice"])
     nsizereg = Int(ds.dim["nsizereg"])
-    ncoeff_ext = Int(ds.dim["ncoeff_ext"])
-    ncoeff_ssa_g = Int(ds.dim["ncoeff_ssa_g"])
-    nbound = Int(ds.dim["nbound"])
-    pair = Int(ds.dim["pair"])
-
     @assert nsizereg == 3 # RRTMGP pade approximation assumes exactly (3) size regimes
 
-    radliq_lwr = FT(Array(ds["radliq_lwr"]))
-    radliq_upr = FT(Array(ds["radliq_upr"]))
-    radliq_fac = FT(Array(ds["radliq_fac"]))
-
-    radice_lwr = FT(Array(ds["radice_lwr"]))
-    radice_upr = FT(Array(ds["radice_upr"]))
-    radice_fac = FT(Array(ds["radice_fac"]))
-
-    lut_extliq = FTA2D(Array(ds["lut_extliq"]))
-    lut_ssaliq = FTA2D(Array(ds["lut_ssaliq"]))
-    lut_asyliq = FTA2D(Array(ds["lut_asyliq"]))
-
-    lut_extice = FTA3D(Array(ds["lut_extice"]))
-    lut_ssaice = FTA3D(Array(ds["lut_ssaice"]))
-    lut_asyice = FTA3D(Array(ds["lut_asyice"]))
-
     pade_extliq = FTA3D(Array(ds["pade_extliq"]))
     pade_ssaliq = FTA3D(Array(ds["pade_ssaliq"]))
     pade_asyliq = FTA3D(Array(ds["pade_asyliq"]))
@@ -925,28 +966,7 @@ function LookUpCld(ds, ::Type{FT}, ::Type{DA}, use_lut::Bool = true) where {FT <
 
     bnd_lims_wn = FTA2D(Array(ds["bnd_limits_wavenumber"]))
 
-    return (LookUpCld{Bool, FT, FTA1D, FTA2D, FTA3D, FTA4D}(
-        nband,
-        nrghice,
-        nsize_liq,
-        nsize_ice,
-        nsizereg,
-        ncoeff_ext,
-        ncoeff_ssa_g,
-        nbound,
-        pair,
-        radliq_lwr,
-        radliq_upr,
-        radliq_fac,
-        radice_lwr,
-        radice_upr,
-        radice_fac,
-        lut_extliq,
-        lut_ssaliq,
-        lut_asyliq,
-        lut_extice,
-        lut_ssaice,
-        lut_asyice,
+    return (PadeCld{FTA1D, FTA2D, FTA3D, FTA4D}(
         pade_extliq,
         pade_ssaliq,
         pade_asyliq,
@@ -960,9 +980,7 @@ function LookUpCld(ds, ::Type{FT}, ::Type{DA}, use_lut::Bool = true) where {FT <
         pade_sizreg_ssaice,
         pade_sizreg_asyice,
         bnd_lims_wn,
-        use_lut,
     ))
-
 end
 
 end
diff --git a/src/optics/Optics.jl b/src/optics/Optics.jl
index 77744c46d..539301df2 100644
--- a/src/optics/Optics.jl
+++ b/src/optics/Optics.jl
@@ -140,7 +140,7 @@ end
         gcol::Int,
         igpt::Int,
         lkp::LookUpLW{FT},
-        lkp_cld::Union{LookUpCld,Nothing} = nothing,
+        lkp_cld::Union{LookUpCld,PadeCld,Nothing} = nothing,
     ) where {FT<:AbstractFloat}
 
 Computes optical properties for the longwave problem.
@@ -172,7 +172,7 @@ end
         gcol::Int,
         igpt::Int,
         lkp::LookUpSW{FT},
-        lkp_cld::Union{LookUpCld,Nothing} = nothing,
+        lkp_cld::Union{LookUpCld,PadeCld,Nothing} = nothing,
     ) where {FT<:AbstractFloat}
 
 Computes optical properties for the shortwave problem.
diff --git a/src/optics/OpticsKernels.jl b/src/optics/OpticsKernels.jl
index b400562b5..7c1bd2bef 100644
--- a/src/optics/OpticsKernels.jl
+++ b/src/optics/OpticsKernels.jl
@@ -12,7 +12,7 @@ function compute_optical_props_kernel!(
     sf::AbstractSourceLW{FT},
     igpt,
     lkp::LookUpLW{FT},
-    lkp_cld::LookUpCld,
+    lkp_cld::Union{LookUpCld, PadeCld},
 ) where {FT <: AbstractFloat}
     ibnd = lkp.major_gpt2bnd[igpt]
 
@@ -43,7 +43,7 @@ function compute_optical_props_kernel!(
     gcol,
     igpt,
     lkp::LookUpSW{FT},
-    lkp_cld::LookUpCld,
+    lkp_cld::Union{LookUpCld, PadeCld},
 ) where {FT <: AbstractFloat}
     ibnd = lkp.major_gpt2bnd[igpt]
 
diff --git a/src/rte/RTESolver.jl b/src/rte/RTESolver.jl
index ef3275da5..8494c1db1 100644
--- a/src/rte/RTESolver.jl
+++ b/src/rte/RTESolver.jl
@@ -24,7 +24,7 @@ include("RTESolverKernels.jl")
         slv::Solver,
         max_threads::Int,
         lookup_lw::Union{LookUpLW, Nothing} = nothing,
-        lookup_lw_cld::Union{LookUpCld, Nothing} = nothing,
+        lookup_lw_cld::Union{LookUpCld, PadeCld, Nothing} = nothing,
     )
 
 Solver for the longwave radiation problem
@@ -33,7 +33,7 @@ function solve_lw!(
     slv::Solver,
     max_threads::Int,
     lookup_lw::Union{LookUpLW, Nothing} = nothing,
-    lookup_lw_cld::Union{LookUpCld, Nothing} = nothing,
+    lookup_lw_cld::Union{LookUpCld, PadeCld, Nothing} = nothing,
 )
     (; as, op, bcs_lw, src_lw, flux_lw, fluxb_lw) = slv
     context = RTE.context(slv)
@@ -63,7 +63,7 @@ end
         max_threads,
         as::AbstractAtmosphericState{FT},
         lookup_lw::Union{LookUpLW, Nothing} = nothing,
-        lookup_lw_cld::Union{LookUpCld, Nothing} = nothing,
+        lookup_lw_cld::Union{LookUpCld, PadeCld, Nothing} = nothing,
     ) where {FT<:AbstractFloat}
 
 No scattering solver for the longwave problem.
@@ -82,7 +82,7 @@ function rte_lw_solve!(
     max_threads,
     as::AbstractAtmosphericState{FT},
     lookup_lw::Union{LookUpLW, Nothing} = nothing,
-    lookup_lw_cld::Union{LookUpCld, Nothing} = nothing,
+    lookup_lw_cld::Union{LookUpCld, PadeCld, Nothing} = nothing,
 ) where {FT <: AbstractFloat}
     (; nlay, ncol) = as
     nlev = nlay + 1
@@ -121,7 +121,7 @@ function rte_lw_noscat_solve_CUDA!(
     major_gpt2bnd,
     as::AbstractAtmosphericState{FT},
     lookup_lw::Union{LookUpLW, Nothing} = nothing,
-    lookup_lw_cld::Union{LookUpCld, Nothing} = nothing,
+    lookup_lw_cld::Union{LookUpCld, PadeCld, Nothing} = nothing,
 ) where {FT <: AbstractFloat}
     gcol = threadIdx().x + (blockIdx().x - 1) * blockDim().x # global id
     nlev = nlay + 1
@@ -174,7 +174,7 @@ function rte_lw_solve!(
     max_threads,
     as::AbstractAtmosphericState{FT},
     lookup_lw::Union{LookUpLW, Nothing} = nothing,
-    lookup_lw_cld::Union{LookUpCld, Nothing} = nothing,
+    lookup_lw_cld::Union{LookUpCld, PadeCld, Nothing} = nothing,
 ) where {FT <: AbstractFloat}
     (; nlay, ncol) = as
     nlev = nlay + 1
@@ -222,7 +222,7 @@ function rte_lw_2stream_solve_CUDA!(
     major_gpt2bnd::AbstractArray{UInt8, 1},
     as::AbstractAtmosphericState{FT},
     lookup_lw::Union{LookUpLW, Nothing} = nothing,
-    lookup_lw_cld::Union{LookUpCld, Nothing} = nothing,
+    lookup_lw_cld::Union{LookUpCld, PadeCld, Nothing} = nothing,
 ) where {FT <: AbstractFloat}
     gcol = threadIdx().x + (blockIdx().x - 1) * blockDim().x # global id
     nlev = nlay + 1
@@ -248,7 +248,7 @@ end
         slv::Solver,
         max_threads::Int,
         lookup_lw::Union{LookUpLW, Nothing} = nothing,
-        lookup_lw_cld::Union{LookUpCld, Nothing} = nothing,
+        lookup_lw_cld::Union{LookUpCld, PadeCld, Nothing} = nothing,
     )
 
 Solver for the shortwave radiation problem
@@ -257,7 +257,7 @@ function solve_sw!(
     slv::Solver,
     max_threads::Int,
     lookup_sw::Union{LookUpSW, Nothing} = nothing,
-    lookup_sw_cld::Union{LookUpCld, Nothing} = nothing,
+    lookup_sw_cld::Union{LookUpCld, PadeCld, Nothing} = nothing,
 )
     (; as, op, bcs_sw, src_sw, flux_sw, fluxb_sw) = slv
     context = RTE.context(slv)
@@ -318,7 +318,7 @@ end
         max_threads,
         as::AbstractAtmosphericState{FT},
         lookup_sw::Union{LookUpSW, Nothing} = nothing,
-        lookup_sw_cld::Union{LookUpCld, Nothing} = nothing,
+        lookup_sw_cld::Union{LookUpCld, PadeCld, Nothing} = nothing,
     ) where {FT<:AbstractFloat}
 
 No-scattering solver for the shortwave problem.
@@ -334,7 +334,7 @@ function rte_sw_noscat_solve!(
     max_threads,
     as::AbstractAtmosphericState{FT},
     lookup_sw::Union{LookUpSW, Nothing} = nothing,
-    lookup_sw_cld::Union{LookUpCld, Nothing} = nothing,
+    lookup_sw_cld::Union{LookUpCld, PadeCld, Nothing} = nothing,
 ) where {FT <: AbstractFloat}
     (; nlay, ncol) = as
     nlev = nlay + 1
@@ -372,7 +372,7 @@ function rte_sw_noscat_solve_CUDA!(
     solar_src_scaled::AbstractArray{FT, 1},
     as::AbstractAtmosphericState{FT},
     lookup_sw::Union{LookUpSW, Nothing} = nothing,
-    lookup_sw_cld::Union{LookUpCld, Nothing} = nothing,
+    lookup_sw_cld::Union{LookUpCld, PadeCld, Nothing} = nothing,
 ) where {FT <: AbstractFloat}
     gcol = threadIdx().x + (blockIdx().x - 1) * blockDim().x # global id
     nlev = nlay + 1
@@ -402,7 +402,7 @@ end
         max_threads,
         as::AbstractAtmosphericState{FT},
         lookup_sw::Union{LookUpSW, Nothing} = nothing,
-        lookup_sw_cld::Union{LookUpCld, Nothing} = nothing,
+        lookup_sw_cld::Union{LookUpCld, PadeCld, Nothing} = nothing,
     ) where {FT<:AbstractFloat}
 
 Two stream solver for the shortwave problem.
@@ -419,7 +419,7 @@ function rte_sw_2stream_solve!(
     max_threads,
     as::AbstractAtmosphericState{FT},
     lookup_sw::Union{LookUpSW, Nothing} = nothing,
-    lookup_sw_cld::Union{LookUpCld, Nothing} = nothing,
+    lookup_sw_cld::Union{LookUpCld, PadeCld, Nothing} = nothing,
 ) where {FT <: AbstractFloat}
     (; nlay, ncol) = as
     nlev = nlay + 1
@@ -482,7 +482,7 @@ function rte_sw_2stream_solve_CUDA!(
     solar_src_scaled::AbstractArray{FT, 1},
     as::AbstractAtmosphericState{FT},
     lookup_sw::Union{LookUpSW, Nothing} = nothing,
-    lookup_sw_cld::Union{LookUpCld, Nothing} = nothing,
+    lookup_sw_cld::Union{LookUpCld, PadeCld, Nothing} = nothing,
 ) where {FT <: AbstractFloat}
     gcol = threadIdx().x + (blockIdx().x - 1) * blockDim().x # global id
     nlev = nlay + 1
diff --git a/test/all_sky_utils.jl b/test/all_sky_utils.jl
index 081c51f99..6c3a82b23 100644
--- a/test/all_sky_utils.jl
+++ b/test/all_sky_utils.jl
@@ -57,7 +57,7 @@ function all_sky(
     close(ds_lw)
     # reading longwave cloud lookup data
     ds_lw_cld = Dataset(lw_cld_file, "r")
-    lookup_lw_cld = LookUpCld(ds_lw_cld, FT, DA, use_lut)
+    lookup_lw_cld = use_lut ? LookUpCld(ds_lw_cld, FT, DA) : PadeCld(ds_lw_cld, FT, DA)
     close(ds_lw_cld)
     #reading shortwave gas optics lookup data
     ds_sw = Dataset(sw_file, "r")
@@ -65,7 +65,7 @@ function all_sky(
     close(ds_sw)
     # reading longwave cloud lookup data
     ds_sw_cld = Dataset(sw_cld_file, "r")
-    lookup_sw_cld = LookUpCld(ds_sw_cld, FT, DA, use_lut)
+    lookup_sw_cld = use_lut ? LookUpCld(ds_sw_cld, FT, DA) : PadeCld(ds_sw_cld, FT, DA)
     close(ds_sw_cld)
     # reading input file 
     ds_in = Dataset(input_file, "r")