diff --git a/notebooks/Growth_with_w.ipynb b/notebooks/Growth_with_w.ipynb
index 793be030b..6ce21766f 100644
--- a/notebooks/Growth_with_w.ipynb
+++ b/notebooks/Growth_with_w.ipynb
@@ -58,7 +58,7 @@
     "def D_hypergeom(avec,csm):\n",
     "    bg = csm.get_background()\n",
     "    Om = csm.Omega0_m()\n",
-    "    if bg.has_key('(.)rho_lambda'):\n",
+    "    if '(.)rho_lambda' in bg:\n",
     "        Ol = bg['(.)rho_lambda'][-1]/bg['(.)rho_crit'][-1]\n",
     "    else:\n",
     "        Ol = bg['(.)rho_fld'][-1]/bg['(.)rho_crit'][-1]\n",
@@ -71,7 +71,7 @@
     "def f_hypergeom(avec,csm):\n",
     "    bg = csm.get_background()\n",
     "    Om = csm.Omega0_m()\n",
-    "    if bg.has_key('(.)rho_lambda'):\n",
+    "    if '(.)rho_lambda' in bg:\n",
     "        Ol = bg['(.)rho_lambda'][-1]/bg['(.)rho_crit'][-1]\n",
     "    else:\n",
     "        Ol = bg['(.)rho_fld'][-1]/bg['(.)rho_crit'][-1]\n",
@@ -84,7 +84,7 @@
     "def D_integral2(avec,csm):\n",
     "    bg = csm.get_background()\n",
     "    Om = csm.Omega0_m()\n",
-    "    if bg.has_key('(.)rho_lambda'):\n",
+    "    if '(.)rho_lambda' in bg:\n",
     "        Ol = bg['(.)rho_lambda'][-1]/bg['(.)rho_crit'][-1]\n",
     "        w0 = -1\n",
     "        wa = 0.0\n",
@@ -122,7 +122,7 @@
     "\n",
     "def D_linder(avec,csm):\n",
     "    bg = csm.get_background()\n",
-    "    if bg.has_key('(.)rho_lambda'):\n",
+    "    if '(.)rho_lambda' in bg:\n",
     "        Ol = bg['(.)rho_lambda'][-1]/bg['(.)rho_crit'][-1]\n",
     "        w0 = -1\n",
     "        wa = 0.0\n",
@@ -152,7 +152,7 @@
     "\n",
     "def D_linder2(avec,csm):\n",
     "    bg = csm.get_background()\n",
-    "    if bg.has_key('(.)rho_lambda'):\n",
+    "    if '(.)rho_lambda' in bg:\n",
     "        Ol = bg['(.)rho_lambda'][-1]/bg['(.)rho_crit'][-1]\n",
     "        w0 = -1\n",
     "        wa = 0.0\n",
@@ -314,7 +314,7 @@
     "fs=14\n",
     "fig, (ax1, ax2) = plt.subplots(2,1,figsize=(6,8),sharex=True,)\n",
     "#                              gridspec_kw = {'height_ratios':[2, 1]})\n",
-    "for M, csm in cosmo.iteritems():\n",
+    "for M, csm in cosmo.items():\n",
     "    if M!='LCDM':\n",
     "        w0, wa = M.strip('()').split(',')\n",
     "        if float(wa)!=0.0:\n",
diff --git a/notebooks/Thomas_PPF.ipynb b/notebooks/Thomas_PPF.ipynb
index 47850661e..3caa869e7 100644
--- a/notebooks/Thomas_PPF.ipynb
+++ b/notebooks/Thomas_PPF.ipynb
@@ -227,7 +227,7 @@
     "            phi_plus_phi_ppf = ptkppf['phi']+ptkppf['psi']\n",
     "            axes[1].semilogx(ptkppf['a'],phi_plus_phi_ppf,label=label+'_ppf')\n",
     "            axes[1].semilogx(ptkfld['a'],phi_plus_phi_fld,label=label+'_fld')\n",
-    "            print len(ptkppf['a']),len(ptkfld['a'])\n",
+    "            print(len(ptkppf['a']), len(ptkfld['a']))\n",
     "      \n",
     "        axes[0].legend(loc='lower left',ncol=2)\n",
     "        axes[0].set_xlim([2,300])\n",
@@ -394,8 +394,8 @@
     }
    ],
    "source": [
-    "print 0.31*0.64**2-0.022\n",
-    "print 0.26*0.74**2-0.022"
+    "print(0.31 * 0.64**2 - 0.022)\n",
+    "print(0.26 * 0.74**2 - 0.022)"
    ]
   },
   {
diff --git a/notebooks/distances.ipynb b/notebooks/distances.ipynb
index 84f660651..deec946ae 100644
--- a/notebooks/distances.ipynb
+++ b/notebooks/distances.ipynb
@@ -69,8 +69,8 @@
     "# Just to cross-check that Omega_Lambda is negligible \n",
     "# (but not exactly zero because we neglected radiation)\n",
     "derived = CDM.get_current_derived_parameters(['Omega0_lambda'])\n",
-    "print derived\n",
-    "print \"Omega_Lambda =\",derived['Omega0_lambda']"
+    "print(derived)\n",
+    "print(\"Omega_Lambda =\",derived['Omega0_lambda'])"
    ]
   },
   {
@@ -86,7 +86,7 @@
     "#Get background quantities and recover their names:\n",
     "baLCDM = LCDM.get_background()\n",
     "baCDM = CDM.get_background()\n",
-    "baCDM.viewkeys()"
+    "baCDM.keys()"
    ]
   },
   {
diff --git a/notebooks/many_times.ipynb b/notebooks/many_times.ipynb
index e8b616538..1149dc7ac 100644
--- a/notebooks/many_times.ipynb
+++ b/notebooks/many_times.ipynb
@@ -133,9 +133,9 @@
     "# (of Hubble crossing, sound horizon crossing, etc.) at different time\n",
     "#\n",
     "background = M.get_background() # load background table\n",
-    "#print background.viewkeys()\n",
+    "#print(background.keys())\n",
     "thermodynamics = M.get_thermodynamics() # load thermodynamics table\n",
-    "#print thermodynamics.viewkeys()    \n",
+    "#print(thermodynamics.keys())\n",
     "#\n",
     "background_tau = background['conf. time [Mpc]'] # read conformal times in background table\n",
     "background_z = background['z'] # read redshift\n",
@@ -173,10 +173,10 @@
     "# check and inform user whether intiial arbitrary choice of z_max_pk was OK\n",
     "max_z_needed = background_z_at_tau(tau[0])\n",
     "if max_z_needed > z_max_pk:\n",
-    "    print ('you must increase the value of z_max_pk to at least ',max_z_needed)\n",
-    "    () + 1  # this strange line is just a trick to stop the script execution there\n",
+    "    print('you must increase the value of z_max_pk to at least ', max_z_needed)\n",
+    "    raise ValueError\n",
     "else:\n",
-    "    print ('in a next run with the same values of tau, you may decrease z_max_pk from ',z_max_pk,' to ',max_z_needed)\n",
+    "    print('in a next run with the same values of tau, you may decrease z_max_pk from ', z_max_pk ,' to ', max_z_needed)\n",
     "#\n",
     "# get transfer functions at each time and build arrays Theta0(tau,k) and phi(tau,k)\n",
     "#\n",
@@ -201,7 +201,7 @@
     "#\n",
     "# inform user of the size of the grids (related to the figure resolution)\n",
     "#\n",
-    "print ('grid size:',len(k),len(tau),Theta0.shape)\n",
+    "print('grid size:', len(k), len(tau), Theta0.shape)\n",
     "#\n",
     "#################\n",
     "#\n",
@@ -214,9 +214,9 @@
     "# plot Theta0(k,tau)\n",
     "#\n",
     "ax_Theta = fig.add_subplot(121)\n",
-    "print ('> Plotting Theta_0')\n",
+    "print('> Plotting Theta_0')\n",
     "fig_Theta = ax_Theta.pcolormesh(K,T,Theta0,cmap='coolwarm',vmin=-Theta_amp, vmax=Theta_amp) #,shading='gouraud')\n",
-    "print ('> Done')\n",
+    "print('> Done')\n",
     "#\n",
     "# plot lines (characteristic times and scales)\n",
     "#\n",
@@ -264,9 +264,9 @@
     "ax_phi = fig.add_subplot(122)\n",
     "ax_phi.set_xlim(k[0],k[-1])\n",
     "#ax_phi.pcolor(K,T,phi,cmap='coolwarm')\n",
-    "print ('> Plotting phi')\n",
+    "print('> Plotting phi')\n",
     "fig_phi = ax_phi.pcolormesh(K,T,phi,cmap='coolwarm',vmin=-0., vmax=phi_amp)\n",
-    "print ('> Done')\n",
+    "print('> Done')\n",
     "#\n",
     "# plot lines (characteristic times and scales)\n",
     "#\n",
diff --git a/notebooks/neutrinohierarchy.ipynb b/notebooks/neutrinohierarchy.ipynb
index deaab14b8..e00cb9b3a 100644
--- a/notebooks/neutrinohierarchy.ipynb
+++ b/notebooks/neutrinohierarchy.ipynb
@@ -88,9 +88,9 @@
    "source": [
     "# test of this function, returning the 3 masses for total mass of 0.1eV\n",
     "m1,m2,m3 = get_masses(2.45e-3,7.50e-5,0.1,'NH')\n",
-    "print 'NH:',m1,m2,m3,m1+m2+m3\n",
+    "print('NH:',m1,m2,m3,m1+m2+m3)\n",
     "m1,m2,m3 = get_masses(2.45e-3,7.50e-5,0.1,'IH')\n",
-    "print 'IH:',m1,m2,m3,m1+m2+m3"
+    "print('IH:',m1,m2,m3,m1+m2+m3)"
    ]
   },
   {
diff --git a/notebooks/one_k.ipynb b/notebooks/one_k.ipynb
index ef81af017..bac55e245 100644
--- a/notebooks/one_k.ipynb
+++ b/notebooks/one_k.ipynb
@@ -86,7 +86,7 @@
     "# load perturbations\n",
     "#\n",
     "all_k = M.get_perturbations()  # this potentially constains scalars/tensors and all k values\n",
-    "print all_k['scalar'][0].viewkeys()\n",
+    "print(all_k['scalar'][0].keys())\n",
     "#    \n",
     "one_k = all_k['scalar'][0]     # this contains only the scalar perturbations for the requested k values\n",
     "#    \n",
@@ -107,14 +107,14 @@
     "# get the time of decoupling\n",
     "#\n",
     "quantities = M.get_current_derived_parameters(['tau_rec'])\n",
-    "# print times.viewkeys()\n",
+    "# print(times.keys())\n",
     "tau_rec = quantities['tau_rec']\n",
     "#\n",
     "# use table of background quantitites to find the time of\n",
     "# Hubble crossing (k / (aH)= 2 pi), sound horizon crossing (k * rs = 2pi)\n",
     "#\n",
     "background = M.get_background() # load background table\n",
-    "#print background.viewkeys()\n",
+    "#print(background.keys())\n",
     "#\n",
     "background_tau = background['conf. time [Mpc]'] # read confromal times in background table\n",
     "background_z = background['z'] # read redshift\n",
diff --git a/notebooks/one_time.ipynb b/notebooks/one_time.ipynb
index 20ae67415..b5cc7da6f 100644
--- a/notebooks/one_time.ipynb
+++ b/notebooks/one_time.ipynb
@@ -95,7 +95,7 @@
     "M.set(common_settings)\n",
     "M.compute()\n",
     "derived = M.get_current_derived_parameters(['z_rec','tau_rec','conformal_age'])\n",
-    "#print derived.viewkeys()\n",
+    "#print(derived.keys())\n",
     "z_rec = derived['z_rec']\n",
     "z_rec = int(1000.*z_rec)/1000. # round down at 4 digits after coma\n",
     "M.struct_cleanup()  # clean output\n",
@@ -111,7 +111,7 @@
     "# load transfer functions at recombination\n",
     "#\n",
     "one_time = M.get_transfer(z_rec)\n",
-    "print one_time.viewkeys()\n",
+    "print(one_time.keys())\n",
     "k = one_time['k (h/Mpc)']\n",
     "Theta0 = 0.25*one_time['d_g']\n",
     "phi = one_time['phi']\n",
@@ -129,7 +129,7 @@
     "# Hubble crossing (k = 2 pi a H), sound horizon crossing (k = 2pi / rs)\n",
     "#\n",
     "background = M.get_background() # load background table\n",
-    "#print background.viewkeys()\n",
+    "#print(background.keys())\n",
     "#\n",
     "background_tau = background['conf. time [Mpc]'] # read confromal times in background table\n",
     "background_z = background['z'] # read redshift\n",
diff --git a/notebooks/thermo.ipynb b/notebooks/thermo.ipynb
index a4cefd48f..0741a2109 100644
--- a/notebooks/thermo.ipynb
+++ b/notebooks/thermo.ipynb
@@ -72,7 +72,7 @@
     "M.compute()\n",
     "derived = M.get_current_derived_parameters(['tau_rec','conformal_age'])\n",
     "thermo = M.get_thermodynamics()\n",
-    "print thermo.viewkeys()"
+    "print(thermo.keys())"
    ]
   },
   {
diff --git a/notebooks/varying_neff.ipynb b/notebooks/varying_neff.ipynb
index 96790df9a..bd71ee52d 100644
--- a/notebooks/varying_neff.ipynb
+++ b/notebooks/varying_neff.ipynb
@@ -52,7 +52,7 @@
     "# i.e. such that omega_ur = omega_gamma * coeff * Neff:\n",
     "# coeff = omega_ur/omega_gamma/Neff_standard\n",
     "coeff = 1.710730e-05/2.472979e-05/3.046\n",
-    "print \"coeff=\",coeff\n",
+    "print(\"coeff=\", coeff)\n",
     "#\n",
     "#############################################\n",
     "#\n",
@@ -84,7 +84,7 @@
     "    alpha = (1.+coeff*N_ur)/(1.+coeff*3.046)\n",
     "    omega_cdm = (0.022032 + 0.12038)*alpha - 0.022032\n",
     "    h = 0.67556*math.sqrt(alpha)\n",
-    "    print ' * Compute with %s=%e, %s=%e, %s=%e'%('N_ur',N_ur,'omega_cdm',omega_cdm,'h',h)\n",
+    "    print(' * Compute with %s=%e, %s=%e, %s=%e' % ('N_ur', N_ur, 'omega_cdm', omega_cdm, 'h', h))\n",
     "    #\n",
     "    # call CLASS\n",
     "    #\n",
diff --git a/notebooks/varying_pann.ipynb b/notebooks/varying_pann.ipynb
index 25115d736..9a510e34e 100644
--- a/notebooks/varying_pann.ipynb
+++ b/notebooks/varying_pann.ipynb
@@ -96,7 +96,7 @@
     "#\n",
     "for i,var in enumerate(var_array):\n",
     "    #\n",
-    "    print ' * Compute with %s=%e'%(var_name,var)\n",
+    "    print(' * Compute with %s=%e'%(var_name,var))\n",
     "    #\n",
     "    # deal with colors and legends\n",
     "    #\n",
diff --git a/notebooks/warmup.ipynb b/notebooks/warmup.ipynb
index 30fb48ec7..8d4d2ad4b 100644
--- a/notebooks/warmup.ipynb
+++ b/notebooks/warmup.ipynb
@@ -36,7 +36,7 @@
     "# get all C_l output\n",
     "cls = LambdaCDM.lensed_cl(2500)\n",
     "# To check the format of cls\n",
-    "cls.viewkeys()"
+    "cls.keys()"
    ]
   },
   {