diff --git a/fgspectra/frequency.py b/fgspectra/frequency.py
index 417120d..dca2a9a 100644
--- a/fgspectra/frequency.py
+++ b/fgspectra/frequency.py
@@ -12,7 +12,6 @@
 import numpy as np
 from scipy import constants
 from .model import Model
-from functools import wraps
 
 
 T_CMB = 2.72548
diff --git a/fgspectra/model.py b/fgspectra/model.py
index 4b50d3d..4bbddbc 100644
--- a/fgspectra/model.py
+++ b/fgspectra/model.py
@@ -3,6 +3,7 @@
 import inspect
 import yaml
 import numpy as np
+import importlib
 
 class Model(ABC):
     """ Abstract class for model definition
@@ -13,11 +14,28 @@ class Model(ABC):
 
     If the `eval` method of a hypotetical ``Child`` class calls the `eval`
     method of other `Model`s, it is likely that ``Child`` has also to override
-    `set_defaults`, `defaults` and `_get_repr`
+    `set_defaults`, `defaults` and `_get_repr`.
+    Note also that if you override `__init__`, you may want to allow to set the
+    defaults at construction time, including the module-level defaults
     """
 
     def __init__(self, **kwargs):
-        """ You can set defaults at construction time """
+        """You can set defaults at construction time using
+        1) the keyword arguments passed at construction time
+        2) if existing, the `defaults` dictionary defined globally in the module
+           of the class being constructed
+
+        1) overrides 2)
+        """
+        self.set_defaults_init(**kwargs)
+
+    def set_defaults_init(self, **kwargs):
+        # As set_defaults, but also considers module-level defaults, if existing
+        try:
+            module_defaults = importlib.import_module(self.__module__).defaults
+            kwargs = {**module_defaults, **kwargs}
+        except AttributeError:
+            pass
         self.set_defaults(**kwargs)
 
     def set_defaults(self, **kwargs):
diff --git a/fgspectra/power.py b/fgspectra/power.py
index 7e5e0da..90815b4 100644
--- a/fgspectra/power.py
+++ b/fgspectra/power.py
@@ -77,7 +77,7 @@ def __init__(self, filenames, **kwargs):
                                   mode='constant', constant_values=0)
 
             self._cl[i+(ell,)] = spec
-        self.set_defaults(**kwargs)
+        self.set_defaults_init(**kwargs)
 
     def eval(self, ell=None, ell_0=None, amp=1.0):
         """Compute the power spectrum with the given ell and parameters."""
@@ -192,6 +192,10 @@ def __init__(self, *power_spectra, **kwargs):
         self.n_comp = np.rint(-1 + np.sqrt(1 + 8 * len(power_spectra))) // 2
         self.n_comp = int(self.n_comp)
         assert (self.n_comp + 1) * self.n_comp // 2 == len(power_spectra)
+        try:
+            kwargs = {**defaults, **kwargs}
+        except NameError:
+            pass
         self.set_defaults(**kwargs)
 
     def set_defaults(self, **kwargs):
@@ -275,6 +279,10 @@ def __init__(self, *power_spectra, **kwargs):
         self.n_comp = np.rint(-1 + np.sqrt(1 + 8 * len(power_spectra))) // 2
         self.n_comp = int(self.n_comp)
         assert (self.n_comp + 1) * self.n_comp // 2 == len(power_spectra)
+        try:
+            kwargs = {**defaults, **kwargs}
+        except NameError:
+            pass
         self.set_defaults(**kwargs)
 
     def set_defaults(self, **kwargs):
diff --git a/tests/test_models.py b/tests/test_models.py
index ea54649..a0749d7 100755
--- a/tests/test_models.py
+++ b/tests/test_models.py
@@ -9,6 +9,7 @@ def test_ksz():
     assert fgp.kSZ_bat() is not None
 
 def test_ACT_models():
+    fgp.defaults = {'ell': np.array([2000])}
     # define the models from fgspectra
     ksz = fgc.FactorizedCrossSpectrum(fgf.ConstantSED(), fgp.kSZ_bat())
     cibp = fgc.FactorizedCrossSpectrum(fgf.ModifiedBlackBody(), fgp.PowerLaw())
@@ -25,8 +26,6 @@ def test_ACT_models():
     tsz = fgc.FactorizedCrossSpectrum(fgf.ThermalSZ(), fgp.tSZ_150_bat())
     cibc = fgc.FactorizedCrossSpectrum(fgf.CIB(), fgp.PowerLaw())
 
-    ells = np.array([2000])
-
 
     par = {
         'nu_0' : 150.0,
@@ -58,18 +57,18 @@ def test_ACT_models():
 
             par['a_tSZ'] * tsz(
                 {'nu':fsz, 'nu_0': par['nu_0']},
-                {'ell':ells, 'ell_0':par['ell_0']}) ,
+                {'ell_0':par['ell_0']}) ,
             par['a_kSZ'] * ksz(
                 {'nu':fsz},
-                {'ell':ells, 'ell_0':par['ell_0']}) ,
+                {'ell_0':par['ell_0']}) ,
 
             par['a_p'] * cibp(
                 {'nu': fdust, 'nu_0':par['nu_0'], 'temp':par['T_d'], 'beta':par['beta_p']},
-                {'ell':ells, 'ell_0':par['ell_0'], 'alpha':2}),
+                {'ell_0':par['ell_0'], 'alpha':2}),
 
             par['a_c'] * cibc(
                 {'nu': fdust, 'nu_0':par['nu_0'], 'temp':par['T_d'], 'beta':par['beta_c']},
-                {'ell':ells, 'ell_0':par['ell_0'], 'alpha':2 - par['n_CIBC']}),
+                {'ell_0':par['ell_0'], 'alpha':2 - par['n_CIBC']}),
 
             tSZ_and_CIB(
                 {'kwseq': (
@@ -77,20 +76,20 @@ def test_ACT_models():
                     {'nu': fdust, 'nu_0':par['nu_0'], 'temp':par['T_d'], 'beta':par['beta_c']} 
                     )},
                 {'kwseq': ( 
-                    {'ell':ells, 'ell_0':par['ell_0'], 
+                    {'ell_0':par['ell_0'], 
                      'amp':par['a_tSZ']},
-                    {'ell':ells, 'ell_0':par['ell_0'], 
+                    {'ell_0':par['ell_0'], 
                      'alpha':2-par['n_CIBC'], 'amp':par['a_c']},
-                    {'ell':ells, 'ell_0':par['ell_0'], 
+                    {'ell_0':par['ell_0'], 
                      'amp': -par['xi'] * np.sqrt(par['a_tSZ'] * par['a_c'])}
                     )}),
 
             par['a_s'] * radio(
                 {'nu': fsynch, 'nu_0':par['nu_0'], 'beta':-0.5 - 2},
-                {'ell':ells, 'ell_0':par['ell_0'], 'alpha':2}) ,
+                {'ell_0':par['ell_0'], 'alpha':2}) ,
             par['a_g'] * cirrus(
                 {'nu': fdust, 'nu_0':par['nu_0'], 'beta':3.8 - 2},
-                {'ell':ells, 'ell_0':par['ell_0'], 'alpha':-0.7})
+                {'ell_0':par['ell_0'], 'alpha':-0.7})
     )
 
     """