diff --git a/python/magic/TOreaders.py b/python/magic/TOreaders.py
index c9e061c6..abf367b9 100644
--- a/python/magic/TOreaders.py
+++ b/python/magic/TOreaders.py
@@ -6,11 +6,12 @@
 import numpy as np
 import matplotlib.pyplot as plt
 from .npfile import npfile
+from .movie import Movie
 from .log import MagicSetup
 from .libmagic import scanDir
 
 
-class TOMovie:
+class TOMovie(Movie):
     """
     This class allows to read and display the :ref:`TO_mov.TAG <secTO_movieFile>`
     generated when :ref:`l_TOmovie=.true. <varl_TOmovie>` is True.
@@ -22,11 +23,17 @@ class TOMovie:
     >>> t = TOMovie(file='TO_mov.N0m2', avg=True, levels=65, cm='seismic')
     """
 
-    def __init__(self, file=None, iplot=True, cm='seismic',
-                 cut=0.5, levels=33, avg=True, precision=np.float32):
+    def __init__(self, file=None, iplot=True, cm='seismic', datadir='.',
+                 cut=0.5, levels=33, avg=True, precision=np.float32,
+                 lastvar=None, nvar='all'):
         """
         :param file: the filename of the TO_mov file
         :type file: str
+        :param nvar: the number of frames of the movie file we want to plot
+                     starting from the last line
+        :type nvar: int
+        :param lastvar: the index of the last timesteps to be read
+        :type lastvar: int
         :param cm: the name of the color map
         :type cm: str
         :param levels: the number of contour levels
@@ -41,6 +48,8 @@ def __init__(self, file=None, iplot=True, cm='seismic',
         :param precision: precision of the input file, np.float32 for single
                           precision, np.float64 for double precision
         :type precision: str
+        :param datadir: working directory
+        :type datadir: str
         """
 
         if file is None:
@@ -58,53 +67,26 @@ def __init__(self, file=None, iplot=True, cm='seismic',
         else:
             filename = file
 
-        mot = re.compile(r'TO_mov\.(.*)')
-        end = mot.findall(filename)[0]
-
-        # DETERMINE THE NUMBER OF LINES BY READING THE LOG FILE
-        logfile = open('log.{}'.format(end), 'r')
-        mot = re.compile(r' ! WRITING TO MOVIE FRAME NO\s*(\d*).*')
-        for line in logfile.readlines():
-            if mot.match(line):
-                nlines = int(mot.findall(line)[0])
-        logfile.close()
+        filename = os.path.join(datadir, filename)
 
-        self.nvar = nlines
+        # Get the version
+        self._get_version(filename)
 
-        # READ the movie file
+        # Read the movie file
         infile = npfile(filename, endian='B')
-        # HEADER
-        version = infile.fort_read('|S64')
-        n_type, n_surface, const, n_fields = infile.fort_read(precision)
-        movtype = infile.fort_read(precision)
-        n_fields = int(n_fields)
-        if n_fields == 8:
+
+        # Header
+        self._read_header(infile, 0, precision)
+        if self.n_fields == 8:
             self.l_phase = True
         else:
             self.l_phase = False
-        self.movtype = np.asarray(movtype)
-        n_surface = int(n_surface)
-
-        # RUN PARAMETERS
-        runid = infile.fort_read('|S64')
-        n_r_mov_tot, n_r_max, n_theta_max, n_phi_tot, self.minc, self.ra, \
-            self.ek, self.pr, self.prmag, \
-            self.radratio, self.tScale = infile.fort_read(precision)
-        n_r_mov_tot = int(n_r_mov_tot)
-        self.n_r_max = int(n_r_max)
-        self.n_theta_max = int(n_theta_max)
-        self.n_phi_tot = int(n_phi_tot)
-
-        # GRID
-        self.radius = infile.fort_read(precision)
-        self.radius = self.radius[:self.n_r_max]  # remove inner core
-        rout = 1./(1.-self.radratio)
-        # rin = self.radratio/(1.-self.radratio)
-        self.radius *= rout
-        self.theta = infile.fort_read(precision)
-        self.phi = infile.fort_read(precision)
 
-        shape = (self.n_theta_max, self.n_r_max)
+        # Get the number of lines
+        self._get_nlines(datadir, filename, nvar, lastvar)
+
+        # Shape
+        self._get_data_shape(precision, allocate=False)
 
         self.time = np.zeros(self.nvar, precision)
         self.asVphi = np.zeros((self.nvar, self.n_theta_max, self.n_r_max),
@@ -118,28 +100,41 @@ def __init__(self, file=None, iplot=True, cm='seismic',
         if self.l_phase:
             self.penalty = np.zeros_like(self.asVphi)
 
-        # READ the data
+        # Read the data
+        # If one skip the beginning, nevertheless read but do not store
+        for i in range(self.var2-self.nvar):
+            if self.version == 2:
+                n_frame, t_movieS, omega_ic, omega_ma, movieDipColat, \
+                    movieDipLon, movieDipStrength, movieDipStrengthGeo = \
+                    infile.fort_read(precision)
+            else:
+                t_movieS = infile.fort_read(precision)
+            for ll in range(self.n_fields):
+                dat = infile.fort_read(precision, shape=self.shape, order='F')
         for k in range(self.nvar):
-            n_frame, t_movieS, omega_ic, omega_ma, movieDipColat, \
-                movieDipLon, movieDipStrength, movieDipStrengthGeo \
-                = infile.fort_read(precision)
+            if self.version == 2:
+                n_frame, t_movieS, omega_ic, omega_ma, movieDipColat, \
+                    movieDipLon, movieDipStrength, movieDipStrengthGeo \
+                    = infile.fort_read(precision)
+            else:
+                t_movieS = infile.fort_read(precision)
             self.time[k] = t_movieS
-            self.asVphi[k, ...] = infile.fort_read(precision, shape=shape,
+            self.asVphi[k, ...] = infile.fort_read(precision, shape=self.shape,
                                                    order='F')
-            self.rey[k, ...] = infile.fort_read(precision, shape=shape,
+            self.rey[k, ...] = infile.fort_read(precision, shape=self.shape,
                                                 order='F')
-            self.adv[k, ...] = infile.fort_read(precision, shape=shape,
+            self.adv[k, ...] = infile.fort_read(precision, shape=self.shape,
                                                 order='F')
-            self.visc[k, ...] = infile.fort_read(precision, shape=shape,
+            self.visc[k, ...] = infile.fort_read(precision, shape=self.shape,
                                                  order='F')
-            self.lorentz[k, ...] = infile.fort_read(precision, shape=shape,
+            self.lorentz[k, ...] = infile.fort_read(precision, shape=self.shape,
                                                     order='F')
-            self.coriolis[k, ...] = infile.fort_read(precision, shape=shape,
+            self.coriolis[k, ...] = infile.fort_read(precision, shape=self.shape,
                                                      order='F')
-            self.dtVp[k, ...] = infile.fort_read(precision, shape=shape,
+            self.dtVp[k, ...] = infile.fort_read(precision, shape=self.shape,
                                                  order='F')
             if self.l_phase:
-                self.penalty[k, ...] = infile.fort_read(precision, shape=shape,
+                self.penalty[k, ...] = infile.fort_read(precision, shape=self.shape,
                                                         order='F')
 
         if iplot:
diff --git a/python/magic/__init__.py b/python/magic/__init__.py
index cca0e6b9..c98a3986 100644
--- a/python/magic/__init__.py
+++ b/python/magic/__init__.py
@@ -23,7 +23,6 @@
         from .potExtra import *
         from .graph2vtk import *
         #from .graph2rst import *
-        from .movie3D import *
         from .spectralTransforms import *
         from .potential import *
     except ImportError as e:
diff --git a/python/magic/averages.py b/python/magic/averages.py
index a899b4a5..54ecddfc 100644
--- a/python/magic/averages.py
+++ b/python/magic/averages.py
@@ -35,7 +35,7 @@ def to_json(o, level=0):
             ret += SPACE * INDENT * (level + 1)
             ret += '"' + str(k) + '":' + SPACE
             ret += to_json(v, level + 1)
-  
+
         ret += NEWLINE + SPACE * INDENT * level + "}"
     elif isinstance(o, str):
         ret += '"' + o + '"'
@@ -74,12 +74,12 @@ class AvgField:
     This class computes the time-average properties from time series, spectra
     and radial profiles. It will store the input starting time in a small file
     named ``tInitAvg``, such that the next time you use it you don't need to
-    provide ``tstart`` again. By default, the outputs are stored in a 
+    provide ``tstart`` again. By default, the outputs are stored in a
     fully documented JSON file named avg.json: this is split into several
     categories, namely numerical parameters, physical parameters, time averaged
     scalar quantities, time averaged spectra and time-averaged radial profiles.
     The quantities stored in the JSON file are entirely controlled by an input
-    model file wich enlists the different quantities of interest. 
+    model file wich enlists the different quantities of interest.
     A default example file named ``model.json`` is provided in
     ``$MAGIC_HOME/python/magic``, and an example of how to build a dedicated one
     is provided below.
@@ -97,7 +97,7 @@ class AvgField:
                        'spectra': {},
                        'radial_profiles': {'powerR': ['viscDiss', 'buoPower']}
                      }
-    >>> # Compute the selected averages in the dirctory mydir                 
+    >>> # Compute the selected averages in the dirctory mydir
     >>> a = AvgField(datadir='mydir', model=json_model)
     """
 
@@ -409,7 +409,7 @@ class AvgStack:
     avg files to compute a global output which summarises all the outputs in one
     single file
 
-    >>> # Simply go through the directories listed in "runs.txt" and produce a 
+    >>> # Simply go through the directories listed in "runs.txt" and produce a
     >>> # local file named "my_results.json"
     >>> st = AvgStack(dirList="runs.txt", filename="my_results.json")
     >>> # Now also recompute each individual avg.json file in each directory
@@ -423,8 +423,8 @@ def __init__(self, dirList='runs.txt', filename='avg.json', datadir='.',
         :param dirList: the filename of the list of directories, lines starting with
                         a hash are omitted by the reader
         :type dirList: str
-        :param filename: 
-        :param recompute: recompute each individual average file in each single 
+        :param filename:
+        :param recompute: recompute each individual average file in each single
                           directory when set to True
         :type recompute: bool
         :param datadir: working directory
diff --git a/python/magic/checker.py b/python/magic/checker.py
index 82f772a0..9f46ee0a 100644
--- a/python/magic/checker.py
+++ b/python/magic/checker.py
@@ -22,7 +22,7 @@ def MagicCheck(tstart=None):
     if the power.TAG and some spectra have been produced in the current directory.
     If in addition the tInitAvg file is also there in the directory it averages
     only from this starting time.
-    
+
     >>> MagicCheck(tstart=10.)
     """
 
@@ -48,7 +48,7 @@ def MagicCheck(tstart=None):
     ohmDiss_avg = avgField(ts.time[mask], ts.ohmDiss[mask])
     viscDiss_avg = avgField(ts.time[mask], ts.viscDiss[mask])
     ratio = 100*abs(buoPower_avg+buoPower_chem_avg+ohmDiss_avg+viscDiss_avg) / \
-            (buoPower_avg+buoPower_chem_avg) 
+            (buoPower_avg+buoPower_chem_avg)
 
     print(bcolors.BOLD + bcolors.UNDERLINE + 'Power balance:' + bcolors.ENDC)
     print('Power injected   : {:.5e}'.format(buoPower_avg+buoPower_chem_avg))
diff --git a/python/magic/checkpoint.py b/python/magic/checkpoint.py
index 294dd7a4..1d46d903 100644
--- a/python/magic/checkpoint.py
+++ b/python/magic/checkpoint.py
@@ -58,7 +58,7 @@ def get_map(lm_max, lmax, mmin, mmax, minc):
 def interp_one_field(field, rold, rnew, rfac=None):
     """
     This routine interpolates a complex input field from an old radial grid
-    to a new one. 
+    to a new one.
 
     :param field: the field to be interpolated
     :type field: numpy.ndarray
@@ -164,7 +164,7 @@ def __init__(self, l_read=True, filename=None, endian='l'):
 
     def read(self, filename, endian='l'):
         """
-        This routine is used to read a checkpoint file. 
+        This routine is used to read a checkpoint file.
 
         :param filename: name of the checkpoint file
         :type filename: str
@@ -179,7 +179,7 @@ def read(self, filename, endian='l'):
         self.version = np.fromfile(file, fmt, count=1)[0]
         fmt = '{}f8'.format(prefix)
         self.time = np.fromfile(file, fmt, count=1)[0]
-        
+
         # Time scheme
         self.tscheme_family = file.read(10).decode()
         nexp, nimp, nold = np.fromfile(file, dtype=np.int32, count=3)
@@ -574,7 +574,7 @@ def xshells2magic(self, xsh_trailing, n_r_max, rscheme='cheb',
 
         # When nalias=60 ntheta=lmax: trick to have lmax in MagIC's header
         self.nalias = 60
-        self.n_theta_max = self.l_max 
+        self.n_theta_max = self.l_max
         self.n_phi_tot = self.l_max
         self.n_r_ic_max = 1
 
@@ -727,7 +727,7 @@ def graph2rst(self, gr, filename='checkpoint_ave.from_chk'):
         s3D = rr3D*np.sin(th3D)
         omr = 1./s3D*(thetaderavg(np.sin(th3D)*gr.vphi, order=4) -
                       phideravg(gr.vtheta, minc=self.minc))
- 
+
         for i in range(self.n_r_max):
             om = sh.spat_spec(omr[:, :, i])
             self.ztor[i, 1:] = om[1:]/(sh.ell[1:]*(sh.ell[1:]+1)) * \
@@ -777,7 +777,7 @@ def graph2rst(self, gr, filename='checkpoint_ave.from_chk'):
                     rr = self.radius_ic[i-1]
                     rdep[1:] = (self.radius_ic[i-1]/ri)**(sh.ell[1:]+1)
                 self.bpol_ic[i, 1:] = vr[1:]/(sh.ell[1:]*(sh.ell[1:]+1))*rr**2
-                # Not stable: 
+                # Not stable:
                 #if self.radius_ic[i] >= 0.01:
                 #    mask = ( self.lm2l <= 2 ) * ( self.lm2m <= 2)
                 #    self.bpol_ic[i, mask] /= rdep[mask]
diff --git a/python/magic/cyl.py b/python/magic/cyl.py
index 70404c0b..73b088f4 100644
--- a/python/magic/cyl.py
+++ b/python/magic/cyl.py
@@ -96,7 +96,7 @@ def sph2cyl_plane(data, rad, ns):
     def zavg(input, radius, ns, minc, save=True, filename='vp.pickle',
              normed=True, colat=None):
         """
-        This function computes a z-integration of a list of input arrays 
+        This function computes a z-integration of a list of input arrays
         (on the spherical grid). This works well for 2-D (phi-slice) arrays.
         In case of 3-D arrays, only one element is allowed (too demanding
         otherwise).
@@ -114,7 +114,7 @@ def zavg(input, radius, ns, minc, save=True, filename='vp.pickle',
         :type save: bool
         :param filename: name of the output pickle when save=True
         :type filename: str
-        :param normed: a boolean to specify if ones wants to simply integrate 
+        :param normed: a boolean to specify if ones wants to simply integrate
                        over z or compute a z-average (default is True: average)
         :type normed: bool
         :param colat: an optional array containing the colatitudes
@@ -186,8 +186,8 @@ def zavg(input, radius, ns, minc, save=True, filename='vp.pickle',
     def zavg(input, radius, ns, minc, save=True, filename='vp.pickle',
              normed=True):
         """
-        This function computes a z-integration of a list of input arrays 
-        (on the spherical grid). This works well for 2-D (phi-slice) 
+        This function computes a z-integration of a list of input arrays
+        (on the spherical grid). This works well for 2-D (phi-slice)
         arrays. In case of 3-D arrays, only one element is allowed
         (too demanding otherwise).
 
diff --git a/python/magic/graph.py b/python/magic/graph.py
index 7c60aa9f..f6480f39 100644
--- a/python/magic/graph.py
+++ b/python/magic/graph.py
@@ -46,7 +46,7 @@ def getGraphEndianness(filename):
     except ValueError:
         f = open(filename, 'rb')
         # Little endian
-        version = np.fromfile(f, 'i4', count=1)[0]
+        version = np.fromfile(f, np.int32, count=1)[0]
         endian = 'l'
         if abs(version) > 100:
             f.close()
diff --git a/python/magic/libmelt.py b/python/magic/libmelt.py
index 81843eaf..b29bd35f 100644
--- a/python/magic/libmelt.py
+++ b/python/magic/libmelt.py
@@ -183,7 +183,7 @@ def stack_two(self, theta_new, time_new, rmelt_new):
         # If the number of latitudinal grid points has changed, one
         # needs to interpolate
         if len(self.colatitude) != len(theta_new):
-            it = interp1d(self.colatitude, self.rmelt, axis=-1, 
+            it = interp1d(self.colatitude, self.rmelt, axis=-1,
                           fill_value='extrapolate')
             self.rmelt = it(theta_new)
             self.colatitude = theta_new
@@ -208,7 +208,7 @@ def plot(self, levels=17, cm='viridis'):
         fig = plt.figure()
         ax = fig.add_subplot(111)
 
-        im = ax.contourf(self.time, self.colatitude, self.rmelt.T, levels, 
+        im = ax.contourf(self.time, self.colatitude, self.rmelt.T, levels,
                          cmap=plt.get_cmap(cm))
         ax.set_xlabel('Time')
         ax.set_ylabel('Colatitude')
diff --git a/python/magic/movie.py b/python/magic/movie.py
index 965eb880..32f3570e 100644
--- a/python/magic/movie.py
+++ b/python/magic/movie.py
@@ -71,14 +71,14 @@ def __init__(self, file=None, iplot=True, nstep=1, png=False,
                  std=False, dpi=80, normRad=False, precision=np.float32,
                  deminc=True, ifield=0, centeredCm=None, datadir='.'):
         """
-        :param nvar: the number of timesteps of the movie file we want to plot
+        :param nvar: the number of frames of the movie file we want to plot
                      starting from the last line
         :type nvar: int
         :param png: if png=True, write the png files instead of display
         :type png: bool
         :param iplot: if iplot=True, display otherwise just read
         :type iplot: bool
-        :param lastvar: the number of the last timesteps to be read
+        :param lastvar: the index of the last timesteps to be read
         :type lastvar: int
         :param nstep: the stepping between two timesteps
         :type nstep: int
@@ -160,196 +160,68 @@ def __init__(self, file=None, iplot=True, nstep=1, png=False,
             filename = file
 
         filename = os.path.join(datadir, filename)
-        mot = re.compile(r'.*[Mm]ov\.(.*)')
-        end = mot.findall(filename)[0]
+
+        # Get the version
+        self._get_version(filename)
 
         # Read the movie file
         infile = npfile(filename, endian='B')
-        # Header
-        version = infile.fort_read('|S64')[0].decode()
-        n_type, n_surface, const, n_fields = infile.fort_read(precision)
-        movtype = infile.fort_read(precision)
-        self.n_fields = int(n_fields)
-        if self.n_fields > 1:
-            print('!!! Warning: several fields in the movie file !!!')
-            print('!!! {} fields !!!'.format(self.n_fields))
-            print('!!! The one displayed is controlled by the    !!!')
-            print('!!! input parameter ifield (=0 by default)    !!!')
-        self.movtype = int(movtype[ifield])
-        n_surface = int(n_surface)
 
-        # Run parameters
-        runid = infile.fort_read('|S64')
-        n_r_mov_tot, n_r_max, n_theta_max, n_phi_tot, self.minc, self.ra, \
-            self.ek, self.pr, self.prmag, self.radratio, self.tScale =   \
-            infile.fort_read(precision)
-        self.minc = int(self.minc)
-        n_r_mov_tot = int(n_r_mov_tot)
-        self.n_r_max = int(n_r_max)
-        self.n_r_ic_max = n_r_mov_tot-self.n_r_max
-        self.n_theta_max = int(n_theta_max)
-        self.n_phi_tot = int(n_phi_tot)
+        # Read the movie header
+        self._read_header(infile, ifield, precision)
 
-        # Grid
-        if self.movtype in [100, 101, 102]:
-            self.cylRad = infile.fort_read(precision)
-            self.n_s_max = len(self.cylRad)
-        else:
-            self.n_s_max = 0
-        self.radius = infile.fort_read(precision)
-        self.radius_ic = np.zeros((self.n_r_ic_max+2), precision)
-        self.radius_ic[:-1] = self.radius[self.n_r_max-1:]
+        # Get the number of lines
+        self._get_nlines(datadir, filename, nvar, lastvar)
 
-        self.radius = self.radius[:self.n_r_max]  # remove inner core
-        # Overwrite radius to ensure double-precision of the
-        # grid (useful for Cheb der)
-        rout = 1./(1.-self.radratio)
-        # rin = self.radratio/(1.-self.radratio)
-        self.radius *= rout
-        self.radius_ic *= rout
-        # self.radius = chebgrid(self.n_r_max-1, rout, rin)
-        self.theta = infile.fort_read(precision)
-        self.phi = infile.fort_read(precision)
-
-        if filename.split('.')[-2].endswith('TO_mov'):
-            l_TOmov = True
-        else:
-            l_TOmov = False
-
-        # Determine the number of lines by reading the log.TAG file
-        logfile = open(os.path.join(datadir, 'log.{}'.format(end)), 'r')
-        mot = re.compile(r'  ! WRITING MOVIE FRAME NO\s*(\d*).*')
-        mot2 = re.compile(r' ! WRITING TO MOVIE FRAME NO\s*(\d*).*')
-        nlines = 0
-        for line in logfile.readlines():
-            if not l_TOmov and mot.match(line):
-                nlines = int(mot.findall(line)[0])
-            elif l_TOmov and mot2.match(line):
-                nlines = int(mot2.findall(line)[0])
-        logfile.close()
-
-        # In case no 'nlines' can be determined from the log file:
-        if nlines == 0:
-            nlines = getNlines(filename, endian='B', precision=precision)
-            nlines -= 8  # Remove 8 lines of header
-            nlines = nlines//(self.n_fields+1)
-
-        if lastvar is None:
-            self.var2 = nlines
-        else:
-            self.var2 = lastvar
-        if str(nvar) == 'all':
-            self.nvar = nlines
-            self.var2 = nlines
-        else:
-            self.nvar = nvar
-
-        if n_surface == 0:
-            self.surftype = '3d volume'
-            if self.movtype in [1, 2, 3, 31]:
-                shape = (self.n_phi_tot, self.n_theta_max, n_r_mov_tot+2)
-            else:
-                shape = (self.n_phi_tot, self.n_theta_max, self.n_r_max)
-            self.data = np.zeros((self.n_fields, self.nvar, self.n_phi_tot,
-                                  self.n_theta_max, self.n_r_max), precision)
-            self.data_ic = np.zeros((self.n_fields, self.nvar, self.n_phi_tot,
-                                    self.n_theta_max, self.n_r_ic_max+2),
-                                    precision)
-        elif n_surface == 1 or n_surface == -1:
-            self.surftype = 'r_constant'
-            shape = (self.n_phi_tot, self.n_theta_max)
-            self.data = np.zeros((self.n_fields, self.nvar, self.n_phi_tot,
-                                  self.n_theta_max), precision)
-        elif n_surface == 2:
-            self.surftype = 'theta_constant'
-            if self.movtype in [1, 2, 3, 14]:  # read inner core
-                shape = (self.n_phi_tot, n_r_mov_tot+2)
-            else:
-                shape = (self.n_phi_tot, self.n_r_max)
-            self.data = np.zeros((self.n_fields, self.nvar, self.n_phi_tot,
-                                  self.n_r_max), precision)
-            self.data_ic = np.zeros((self.n_fields, self.nvar, self.n_phi_tot,
-                                     self.n_r_ic_max+2), precision)
-        elif n_surface == -2:
-            self.surftype = 'theta_constant'
-            shape = (self.n_phi_tot, self.n_s_max)
-            self.data = np.zeros((self.n_fields, self.nvar, self.n_phi_tot,
-                                  self.n_s_max), precision)
-        elif n_surface >= 3:
-            self.surftype = 'phi_constant'
-            if self.movtype in [1, 2, 3, 14]:  # read inner core
-                shape = (self.n_theta_max, 2*(n_r_mov_tot+2))
-                self.n_theta_plot = 2*self.n_theta_max
-            elif self.movtype in [8, 9, 20, 21, 22, 23, 24, 25, 26]:
-                shape = (self.n_theta_max, n_r_mov_tot+2)
-                self.n_theta_plot = self.n_theta_max
-            elif self.movtype in [4, 5, 6, 7, 15, 16, 17, 18, 47, 54,
-                                  109, 112]:
-                shape = (self.n_theta_max, self.n_r_max, 2)
-                self.n_theta_plot = 2*self.n_theta_max
-            elif self.movtype in [10, 11, 12, 19, 92, 94, 95, 110, 111, 114,
-                                  115, 116]:
-                shape = (self.n_theta_max, self.n_r_max)
-                self.n_theta_plot = self.n_theta_max
-            # Inner core is not stored here
-            self.data = np.zeros((self.n_fields, self.nvar, self.n_theta_plot,
-                                  self.n_r_max), precision)
-            self.data_ic = np.zeros((self.n_fields, self.nvar,
-                                     self.n_theta_plot, self.n_r_ic_max+2),
-                                    precision)
+        # Get the shape of the data
+        self._get_data_shape(precision)
 
         self.time = np.zeros(self.nvar, precision)
 
         # Read the data
-
         # If one skip the beginning, nevertheless read but do not store
         for i in range(self.var2-self.nvar):
-            n_frame, t_movieS, omega_ic, omega_ma, movieDipColat, \
-                movieDipLon, movieDipStrength, movieDipStrengthGeo = \
-                infile.fort_read(precision)
+            if self.version == 2:
+                n_frame, t_movieS, omega_ic, omega_ma, movieDipColat, \
+                    movieDipLon, movieDipStrength, movieDipStrengthGeo = \
+                    infile.fort_read(precision)
+            else:
+                t_movieS = infile.fort_read(precision)
             for ll in range(self.n_fields):
-                dat = infile.fort_read(precision, shape=shape, order='F')
+                dat = infile.fort_read(precision, shape=self.shape, order='F')
         # then read the remaining requested nvar lines
         for k in range(self.nvar):
-            n_frame, t_movieS, omega_ic, omega_ma, movieDipColat, \
-                movieDipLon, movieDipStrength, movieDipStrengthGeo = \
-                infile.fort_read(precision)
+            if self.version == 2:
+                n_frame, t_movieS, omega_ic, omega_ma, movieDipColat, \
+                    movieDipLon, movieDipStrength, movieDipStrengthGeo = \
+                    infile.fort_read(precision)
+            else:
+                t_movieS = infile.fort_read(precision)
             self.time[k] = t_movieS
             for ll in range(self.n_fields):
-                dat = infile.fort_read(precision, shape=shape, order='F')
-                if n_surface == 0:
-                    if self.movtype in [1, 2, 3, 31]:
-                        self.data[ll, k, ...] = dat[:, :, :self.n_r_max]
+                dat = infile.fort_read(precision, shape=self.shape, order='F')
+                if self.n_surface == 0:
+                    self.data[ll, k, ...] = dat[:, :, :self.n_r_max]
+                    if self.lIC:
                         self.data_ic[ll, k, ...] = dat[:, :, self.n_r_max:]
-                    else:
-                        self.data[ll, k, ...] = dat
-                elif n_surface == 2:
-                    if self.movtype in [1, 2, 3, 14]:
-                        self.data[ll, k, ...] = dat[:, :self.n_r_max]
+                elif self.n_surface == 2:
+                    self.data[ll, k, ...] = dat[:, :self.n_r_max]
+                    if self.lIC:
                         self.data_ic[ll, k, ...] = dat[:, self.n_r_max:]
-                    else:
-                        self.data[ll, k, ...] = dat
-                elif n_surface == -2:
+                elif self.n_surface == -2:
                     self.data[ll, k, ...] = dat
-                elif n_surface >= 3:
-                    if self.movtype in [1, 2, 3, 14]:
-                        datoc0 = dat[:, :self.n_r_max]
-                        datoc1 = dat[:, self.n_r_max:2*self.n_r_max]
-                        datic0 = dat[:, 2*self.n_r_max:2*self.n_r_max+self.n_r_ic_max+2]
-                        datic1 = dat[:, 2*self.n_r_max+self.n_r_ic_max+2:]
-                        self.data[ll, k, ...] = np.vstack((datoc0, datoc1))
+                elif self.n_surface == 3:
+                    datoc0 = dat[:, :self.n_r_max]
+                    datoc1 = dat[:, self.n_r_max:2*self.n_r_max]
+                    self.data[ll, k, ...] = np.vstack((datoc0, datoc1))
+                    if self.lIC:
+                        datic0 = dat[:, 2*self.n_r_max:2*self.n_r_max+self.n_r_ic_max]
+                        datic1 = dat[:, 2*self.n_r_max+self.n_r_ic_max:]
                         self.data_ic[ll, k, ...] = np.vstack((datic0, datic1))
-                    elif self.movtype in [8, 9, 20, 21, 22, 23, 24, 25, 26]:
+                elif self.n_surface == 4:
+                    self.data[ll, k, ...] = dat[:, :self.n_r_max]
+                    if self.lIC:
                         self.data_ic[ll, k, ...] = dat[:, self.n_r_max:]
-                        self.data[ll, k, ...] = dat[:, :self.n_r_max]
-                    elif self.movtype in [4, 5, 6, 7, 15, 16, 17, 18, 47, 54, 91,
-                                          109, 112]:
-                        dat0 = dat[..., 0]
-                        dat1 = dat[..., 1]
-                        self.data[ll, k, ...] = np.vstack((dat0, dat1))
-                    elif self.movtype in [10, 11, 12, 19, 92, 94, 95, 110, 111,
-                                          114, 115, 116]:
-                        self.data[ll, k, ...] = dat
                 else:
                     self.data[ll, k, ...] = dat
                 if fluct:
@@ -379,6 +251,226 @@ def __init__(self, file=None, iplot=True, nstep=1, png=False,
             cmap = plt.get_cmap(cm)
             self.avgStd(ifield, std, cut, centeredCm, levels, cmap)
 
+    def _get_version(self, filename):
+        """
+        This routine determines the version of the movie files
+
+        :param filename: name of the movie file
+        :type filename: str
+        """
+        with open(filename, 'rb') as fi:
+            rm = np.fromfile(fi, np.float32, count=1) # record marker
+            self.version = np.fromfile(fi, '>i4', count=1)[0]
+            rm = np.fromfile(fi, np.float32, count=1)
+        if abs(self.version) > 100:
+            fi = npfile(filename, endian='B')
+            version = fi.fort_read('|S64')[0].decode().rstrip()
+            self.version = int(version[-1])
+            fi.close()
+
+    def _get_nlines(self, datadir, filename, nvar, lastvar):
+        """
+        This routine determines the number of frames stored in a movie file
+
+        :param filename: name of the movie file
+        :type filename: str
+        :param datadir: working directory
+        :type datadir: str
+        :param nvar: the number of frames of the movie file we want to plot
+                     starting from the last line
+        :type nvar: int
+        :param lastvar: the index of the last timesteps to be read
+        :type lastvar: int
+        """
+
+        # Get the number of lines
+        pattern = re.compile(r'.*[Mm]ov\.(.*)')
+        end = pattern.findall(filename)[0]
+        if filename.split('.')[-2].endswith('TO_mov'):
+            l_TOmov = True
+        else:
+            l_TOmov = False
+
+        # Determine the number of lines by reading the log.TAG file
+        logfile = open(os.path.join(datadir, 'log.{}'.format(end)), 'r')
+        pattern = re.compile(r'  ! WRITING MOVIE FRAME NO\s*(\d*).*')
+        mot2 = re.compile(r' ! WRITING TO MOVIE FRAME NO\s*(\d*).*')
+        nlines = 0
+        for line in logfile.readlines():
+            if not l_TOmov and pattern.match(line):
+                nlines = int(pattern.findall(line)[0])
+            elif l_TOmov and mot2.match(line):
+                nlines = int(mot2.findall(line)[0])
+        logfile.close()
+
+        # In case no 'nlines' can be determined from the log file:
+        if nlines == 0:
+            nlines = getNlines(filename, endian='B', precision=precision)
+            if self.version == 2:
+                nlines -= 8  # Remove 8 lines of header
+            else:
+                nlines -= 10 # Remove 10 lines of header
+            if self.movtype in [100, 101, 102]:
+                nlines -=1  # One additional line in the header for those movies
+            nlines = nlines//(self.n_fields+1)
+
+        if lastvar is None:
+            self.var2 = nlines
+        else:
+            self.var2 = lastvar
+        if str(nvar) == 'all':
+            self.nvar = nlines
+            self.var2 = nlines
+        else:
+            self.nvar = nvar
+
+    def _read_header(self, infile, ifield, precision):
+        """
+        This routine reads the header of the movie files
+
+        :param infile: the movie file loaded using npfile
+        :type infile: magic.npfile.npfile
+        :param ifield: the index of the considered field (in case the movie
+                       contains several)
+        :type ifield: int
+        :param precision: precision of the input file, np.float32 for single
+                          precision, np.float64 for double precision
+        :type precision: str
+        """
+        # Header
+        if self.version == 2:
+            version = infile.fort_read('|S64')[0].decode().rstrip()
+            n_type, n_surface, const, n_fields = infile.fort_read(precision)
+            self.n_surface = int(n_surface)
+            self.n_fields = int(n_fields)
+            movtype = infile.fort_read(precision).astype(np.int32)
+        else:
+            version = infile.fort_read(np.int32)[0]
+            self.n_surface, self.n_fields = infile.fort_read(np.int32)
+            const = infile.fort_read(precision)
+            movtype = infile.fort_read(np.int32)
+        self.movtype = movtype[ifield]
+
+        # Run parameters
+        runid = infile.fort_read('|S64')
+        if self.version == 2:
+            n_r_mov_tot, n_r_max, n_theta_max, n_phi_tot, self.minc, self.ra, \
+                self.ek, self.pr, self.prmag, self.radratio, self.tScale =   \
+                infile.fort_read(precision)
+            self.minc = int(self.minc)
+            self.n_r_mov_tot = int(n_r_mov_tot)
+            self.n_r_max = int(n_r_max)
+            self.n_r_ic_max = self.n_r_mov_tot-self.n_r_max
+            self.n_r_mov_tot += 2
+            self.n_r_ic_max += 2
+            self.raxi = 0.
+            self.sc = 0.
+            self.n_theta_max = int(n_theta_max)
+            self.n_phi_tot = int(n_phi_tot)
+            self.n_s_max = 0
+        else:
+            self.n_r_mov_tot, self.n_r_max, self.n_r_ic_max, self.n_theta_max, \
+                self.n_phi_tot, self.n_s_max, self.minc = infile.fort_read(np.int32)
+            self.ra, self.ek, self.pr, self.raxi, self.sc, self.prmag, \
+                self.radratio, self.tScale = infile.fort_read(precision)
+
+        # Grid
+        if self.movtype in [100, 101, 102]:
+            self.cylRad = infile.fort_read(precision)
+            if self.version == 2:
+                self.n_s_max = len(self.cylRad)
+        self.radius = infile.fort_read(precision)
+        if self.version == 2:
+            self.radius_ic = np.zeros((self.n_r_ic_max), precision)
+            self.radius_ic[:-1] = self.radius[self.n_r_max-1:]
+        else:
+            self.radius_ic = self.radius[self.n_r_max:]
+
+        self.radius = self.radius[:self.n_r_max]  # remove inner core
+        # Overwrite radius to ensure double-precision of the
+        # grid (useful for Cheb der)
+        rout = 1./(1.-self.radratio)
+        # rin = self.radratio/(1.-self.radratio)
+        self.radius *= rout
+        self.radius_ic *= rout
+        # self.radius = chebgrid(self.n_r_max-1, rout, rin)
+        self.theta = infile.fort_read(precision)
+        self.phi = infile.fort_read(precision)
+
+        if self.version == 2:
+            if self.movtype not in [1, 2, 3, 14, 8, 9, 20, 21, 22, 23, 24, 25, 26, 31]:
+                self.n_r_mov_tot = self.n_r_max
+                self.lIC = False
+            else:
+                self.lIC = True
+        else:
+            if self.n_r_mov_tot > self.n_r_max:
+                self.lIC = True
+            else:
+                self.lIC = False
+
+    def _get_data_shape(self, precision, allocate=True):
+        """
+        This determines the size of the movie frames depending on n_surface
+
+        :param precision: precision of the input file, np.float32 for single
+                          precision, np.float64 for double precision
+        :type precision: str
+        """
+
+        if self.n_surface == 0:
+            self.shape = (self.n_phi_tot, self.n_theta_max, self.n_r_mov_tot)
+            if allocate:
+                self.data = np.zeros((self.n_fields, self.nvar, self.n_phi_tot,
+                                      self.n_theta_max, self.n_r_max), precision)
+            if allocate and self.lIC:
+                self.data_ic = np.zeros((self.n_fields, self.nvar, self.n_phi_tot,
+                                        self.n_theta_max, self.n_r_ic_max),
+                                        precision)
+        elif self.n_surface == 1 or self.n_surface == -1:
+            self.shape = (self.n_phi_tot, self.n_theta_max)
+            if allocate:
+                self.data = np.zeros((self.n_fields, self.nvar, self.n_phi_tot,
+                                      self.n_theta_max), precision)
+        elif self.n_surface == 2:
+            self.shape = (self.n_phi_tot, self.n_r_mov_tot)
+            if allocate:
+                self.data = np.zeros((self.n_fields, self.nvar, self.n_phi_tot,
+                                      self.n_r_max), precision)
+            if allocate and self.lIC:
+                self.data_ic = np.zeros((self.n_fields, self.nvar, self.n_phi_tot,
+                                         self.n_r_ic_max), precision)
+        elif self.n_surface == -2:
+            self.shape = (self.n_phi_tot, self.n_s_max)
+            if allocate:
+                self.data = np.zeros((self.n_fields, self.nvar, self.n_phi_tot,
+                                      self.n_s_max), precision)
+        elif self.n_surface == 3:
+            if self.movtype in [8, 9, 10, 11, 12, 19, 20, 21, 22, 23, 24, 25, 26,
+                                92, 94, 95, 110, 111, 114, 115, 116]:
+                self.shape = (self.n_theta_max, self.n_r_mov_tot)
+                self.n_theta_plot = self.n_theta_max
+                self.n_surface == 4
+            else:
+                self.shape = (self.n_theta_max, 2*self.n_r_mov_tot)
+                self.n_theta_plot = 2*self.n_theta_max
+            if allocate:
+                self.data = np.zeros((self.n_fields, self.nvar, self.n_theta_plot,
+                                      self.n_r_max), precision)
+            if allocate and self.lIC:
+                self.data_ic = np.zeros((self.n_fields, self.nvar,
+                                         self.n_theta_plot, self.n_r_ic_max),
+                                        precision)
+        elif self.n_surface == 4:
+            self.shape = (self.n_theta_max, self.n_r_mov_tot)
+            if allocate:
+                self.data = np.zeros((self.n_fields, self.nvar, self.n_theta_max,
+                                      self.n_r_max), precision)
+            if allocate and self.lIC:
+                self.data_ic = np.zeros((self.n_fields, self.nvar,
+                                         self.n_theta_max, self.n_r_ic_max),
+                                        precision)
+
     def __add__(self, new):
         """
         Built-in function to sum two movies. In case, the spatial grid have been
@@ -393,7 +485,7 @@ def __add__(self, new):
         if self.data[0, 0, ...].shape != new.data[0, 0, ...].shape:
             new_shape = new.data[0, 0, ...].shape
             old_shape = self.data[0, 0, ...].shape
-            if self.surftype == 'r_constant':
+            if self.n_surface == 1:
                 if (new_shape[0] != old_shape[0]) and (new_shape[1] != old_shape[1]):
                     ip = interp1d(self.phi, self.data, axis=-2,
                                   fill_value='extrapolate')
@@ -409,7 +501,7 @@ def __add__(self, new):
                     it = interp1d(self.theta, self.data, axis=-1,
                                   fill_value='extrapolate')
                     self.data = it(new.theta)
-            elif self.surftype == 'theta_constant':
+            elif self.n_surface == 2:
                 if (new_shape[0] != old_shape[0]) and (new_shape[1] != old_shape[1]):
                     ip = interp1d(self.phi, self.data, axis=-2,
                                   fill_value='extrapolate')
@@ -425,8 +517,7 @@ def __add__(self, new):
                     ir = interp1d(self.radius[::-1], self.data[..., ::-1], axis=-1)
                     tmp = ir(new.radius[::-1])
                     self.data = self.data[..., ::-1]
-            elif self.surftype == 'phi_constant' and \
-                 self.movtype in [10, 11, 12, 19, 92, 94, 95, 110, 111]:
+            elif self.n_surface == 4:
                 if (new_shape[0] != old_shape[0]) and (new_shape[1] != old_shape[1]):
                     it = interp1d(self.theta, self.data, axis=-2,
                                   fill_value='extrapolate')
@@ -442,8 +533,7 @@ def __add__(self, new):
                     ir = interp1d(self.radius[::-1], self.data[..., ::-1], axis=-1)
                     tmp = ir(new.radius[::-1])
                     self.data = tmp[..., ::-1]
-            elif self.surftype == 'phi_constant' and \
-                 self.movtype not in [10, 11, 12, 19, 92, 94, 95, 110, 111]:
+            elif self.n_surface == 3:
                 if (new_shape[0] != old_shape[0]) and (new_shape[1] != old_shape[1]):
                     it = interp1d(self.theta, self.data[..., :self.n_theta_max, :],
                                   axis=-2, fill_value='extrapolate')
@@ -483,7 +573,7 @@ def __add__(self, new):
         return out
 
     def avgStd(self, ifield=0, std=False, cut=0.5, centeredCm=None,
-               levels=12, cmap='RdYlBu_r', ic=False):
+               levels=12, cmap='RdYlBu_r'):
         """
         Plot time-average or standard deviation
 
@@ -505,11 +595,11 @@ def avgStd(self, ifield=0, std=False, cut=0.5, centeredCm=None,
         """
         if std:
             avg = self.data[ifield, ...].std(axis=0)
-            if ic:
+            if self.lIC:
                 avg_ic = self.data_ic[ifield, ...].std(axis=0)
         else:
             avg = self.data[ifield, ...].mean(axis=0)
-            if ic:
+            if self.lIC:
                 avg_ic = self.data_ic[ifield, ...].mean(axis=0)
         if centeredCm is None:
             if avg.min() < 0 and avg.max() > 0:
@@ -525,12 +615,12 @@ def avgStd(self, ifield=0, std=False, cut=0.5, centeredCm=None,
             vmin = cut * avg.min()
         cs = np.linspace(vmin, vmax, levels)
 
-        if self.surftype == 'phi_constant':
-            if self.n_theta_plot == self.n_theta_max:
+        if self.n_surface in [3, 4]:
+            if self.n_surface == 4:
                 th = np.linspace(np.pi/2., -np.pi/2., self.n_theta_max)
                 fig = plt.figure(figsize=(4, 8))
                 th0 = th
-            else:
+            elif self.n_surface == 3:
                 th0 = np.linspace(np.pi/2., -np.pi/2., self.n_theta_max)
                 th1 = np.linspace(np.pi/2., 3.*np.pi/2., self.n_theta_max)
                 th = np.concatenate((th0, th1))
@@ -544,11 +634,11 @@ def avgStd(self, ifield=0, std=False, cut=0.5, centeredCm=None,
             yyout = rr.max() * np.sin(th0)
             xxin = rr.min() * np.cos(th0)
             yyin = rr.min() * np.sin(th0)
-            if ic:
+            if self.lIC:
                 rr, tth = np.meshgrid(self.radius_ic, th)
                 xx_ic = rr * np.cos(tth)
                 yy_ic = rr * np.sin(tth)
-        elif self.surftype == 'r_constant':
+        elif self.n_surface == 1:
             th = np.linspace(np.pi/2., -np.pi/2., self.n_theta_max)
             phi = np.linspace(-np.pi, np.pi, self.n_phi_tot)
             ttheta, pphi = np.meshgrid(th, phi)
@@ -556,7 +646,7 @@ def avgStd(self, ifield=0, std=False, cut=0.5, centeredCm=None,
             xxout, yyout = hammer2cart(th, -np.pi)
             xxin, yyin = hammer2cart(th, np.pi)
             fig = plt.figure(figsize=(8, 4))
-        elif self.surftype == 'theta_constant':
+        elif self.n_surface == 2:
             phi = np.linspace(0., 2.*np.pi, self.n_phi_tot)
             rr, pphi = np.meshgrid(self.radius, phi)
             xx = rr * np.cos(pphi)
@@ -565,8 +655,12 @@ def avgStd(self, ifield=0, std=False, cut=0.5, centeredCm=None,
             yyout = rr.max() * np.sin(pphi)
             xxin = rr.min() * np.cos(pphi)
             yyin = rr.min() * np.sin(pphi)
+            if self.lIC:
+                rr, pphi = np.meshgrid(self.radius_ic, phi)
+                xx_ic = rr * np.cos(pphi)
+                yy_ic = rr * np.sin(pphi)
             fig = plt.figure(figsize=(6, 6))
-        elif self.surftype == '3d volume':
+        elif n_surface == 0:
             self.data = self.data[ifield, ..., 0]
             th = np.linspace(np.pi/2., -np.pi/2., self.n_theta_max)
             phi = np.linspace(-np.pi, np.pi, self.n_phi_tot)
@@ -579,7 +673,7 @@ def avgStd(self, ifield=0, std=False, cut=0.5, centeredCm=None,
         fig.subplots_adjust(top=0.99, right=0.99, bottom=0.01, left=0.01)
         ax = fig.add_subplot(111, frameon=False)
         ax.contourf(xx, yy, avg, cs, cmap=cmap, extend='both')
-        if ic:
+        if self.lIC:
             ax.contourf(xx_ic, yy_ic, avg_ic, cs, cmap=cmap, extend='both')
         ax.plot(xxout, yyout, 'k-', lw=1.5)
         ax.plot(xxin, yyin, 'k-', lw=1.5)
@@ -646,12 +740,12 @@ def plot(self, ifield=0, cut=0.5, centeredCm=None, levels=12,
             # vmin, vmax = self.data.min(), self.data.max()
             cs = np.linspace(vmin, vmax, levels)
 
-        if self.surftype == 'phi_constant':
-            if self.n_theta_plot == self.n_theta_max:
+        if self.n_surface in [3, 4]:
+            if self.n_surface == 4:
                 th = np.linspace(np.pi/2., -np.pi/2., self.n_theta_max)
                 fig = plt.figure(figsize=(4, 8))
                 th0 = th
-            else:
+            elif self.n_surface == 3:
                 th0 = np.linspace(np.pi/2., -np.pi/2., self.n_theta_max)
                 th1 = np.linspace(np.pi/2., 3.*np.pi/2., self.n_theta_max)
                 th = np.concatenate((th0, th1))
@@ -670,7 +764,7 @@ def plot(self, ifield=0, cut=0.5, centeredCm=None, levels=12,
                 rr, tth = np.meshgrid(self.radius_ic, th)
                 xx_ic = rr * np.cos(tth)
                 yy_ic = rr * np.sin(tth)
-        elif self.surftype == 'r_constant':
+        elif self.n_surface == 1:
             th = np.linspace(np.pi/2., -np.pi/2., self.n_theta_max)
             if deminc:
                 phi = np.linspace(-np.pi, np.pi, self.n_phi_tot*self.minc+1)
@@ -684,7 +778,7 @@ def plot(self, ifield=0, cut=0.5, centeredCm=None, levels=12,
             ttheta, pphi = np.meshgrid(th, phi)
             xx, yy = hammer2cart(ttheta, pphi)
             fig = plt.figure(figsize=(8, 4))
-        elif self.surftype == 'theta_constant':
+        elif self.n_surface == 2:
             if deminc:
                 phi = np.linspace(0., 2.*np.pi, self.n_phi_tot*self.minc+1)
             else:
@@ -704,7 +798,7 @@ def plot(self, ifield=0, cut=0.5, centeredCm=None, levels=12,
                 xx_ic = rr * np.cos(pphi)
                 yy_ic = rr * np.sin(pphi)
             fig = plt.figure(figsize=(6, 6))
-        elif self.surftype == '3d volume':
+        elif self.n_surface == 0:
             self.data = self.data[ifield, ..., 0]
             th = np.linspace(np.pi/2., -np.pi/2., self.n_theta_max)
             phi = np.linspace(-np.pi, np.pi, self.n_phi_tot)
@@ -725,7 +819,7 @@ def plot(self, ifield=0, cut=0.5, centeredCm=None, levels=12,
                     vmin = cut * vmin
                     vmax = -vmin
                     cs = np.linspace(vmin, vmax, levels)
-                if self.surftype in ['r_constant', 'theta_constant']:
+                if self.n_surface in [1, 2]:
                     if deminc:
                         dat = symmetrize(self.data[ifield, k, ...], self.minc)
                         if ic:
@@ -765,7 +859,7 @@ def plot(self, ifield=0, cut=0.5, centeredCm=None, levels=12,
                         vmax = cut * self.data[ifield, k, ...].max()
                         vmin = cut * self.data[ifield, k, ...].min()
                     cs = np.linspace(vmin, vmax, levels)
-                if self.surftype in ['r_constant', 'theta_constant']:
+                if self.n_surface in [1, 2]:
                     if deminc:
                         dat = symmetrize(self.data[ifield, k, ...], self.minc)
                         if ic:
diff --git a/python/magic/movie2vtk.py b/python/magic/movie2vtk.py
index eccbec09..2644790c 100644
--- a/python/magic/movie2vtk.py
+++ b/python/magic/movie2vtk.py
@@ -20,7 +20,7 @@
     print("movie2vtk requires the use of evtk library!")
     print("You can get it from https://github.com/paulo-herrera/PyEVTK")
 
-class Movie2Vtk:
+class Movie2Vtk(Movie):
     """
     This class allows to transform an input Movie file
     :ref:`movie files <secMovieFile>` into a series of vts files that
@@ -125,110 +125,19 @@ def __init__(self, file=None, step=1, lastvar=None, nvar='all', fluct=False,
 
         # Read the movie file
         infile = npfile(os.path.join(datadir, filename), endian='B')
+
         # Header
-        version = infile.fort_read('|S64')
-        n_type, n_surface, const, n_fields = infile.fort_read(precision)
-        movtype = infile.fort_read(precision)
-        self.n_fields = int(n_fields)
-        if self.n_fields > 1:
-            print('!!! Warning: several fields in the movie file !!!')
-            print('!!! {} fields !!!'.format(self.n_fields))
-            print('!!! The one displayed is controlled by the    !!!')
-            print('!!! input parameter ifield (=0 by default)    !!!')
-        self.movtype = int(movtype[ifield])
-        n_surface = int(n_surface)
-
-        # Run parameters
-        runid = infile.fort_read('|S64')
-        n_r_mov_tot, n_r_max, n_theta_max, n_phi_tot, self.minc, self.ra, \
-            self.ek, self.pr, self.prmag, self.radratio, self.tScale =   \
-            infile.fort_read(precision)
-        self.minc = int(self.minc)
-        n_r_mov_tot = int(n_r_mov_tot)
-        self.n_r_max = int(n_r_max)
-        self.n_r_ic_max = n_r_mov_tot-self.n_r_max
-        self.n_theta_max = int(n_theta_max)
-        self.n_phi_tot = int(n_phi_tot)
-
-        # Grid
-        self.radius = infile.fort_read(precision)
-        self.radius_ic = np.zeros((self.n_r_ic_max+2), precision)
-        self.radius_ic[:-1] = self.radius[self.n_r_max-1:]
-
-        self.radius = self.radius[:self.n_r_max]  # remove inner core
-        # Overwrite radius to ensure double-precision of the
-        # grid (useful for Cheb der)
-        rout = 1./(1.-self.radratio)
-        self.radius *= rout
-        self.radius_ic *= rout
-        # rin = self.radratio/(1.-self.radratio)
-        # self.radius = chebgrid(self.n_r_max-1, rout, rin)
-        self.theta = infile.fort_read(precision)
+        self._read_header(infile, ifield, precision)
         if closePoles:
             self.theta = np.linspace(0., np.pi, self.n_theta_max)
-        self.phi = infile.fort_read(precision)
-
-        # Determine the number of lines by reading the log.TAG file
-        with open(os.path.join(datadir, 'log.{}'.format(end)), 'r') as logfile:
-            mot = re.compile(r'  ! WRITING MOVIE FRAME NO\s*(\d*).*')
-            mot2 = re.compile(r' ! WRITING TO MOVIE FRAME NO\s*(\d*).*')
-            nlines = 0
-            for line in logfile.readlines():
-                if mot.match(line):
-                    nlines = int(mot.findall(line)[0])
-                elif mot2.match(line):
-                    nlines = int(mot2.findall(line)[0])
-
-        # In case no 'nlines' can be determined from the log file:
-        if nlines == 0:
-            nlines = getNlines(os.path.join(datadir, filename),
-                               endian='B', precision=precision)
-            nlines -= 8  # Remove 8 lines of header
-            nlines = nlines // (self.n_fields+1)
-
-        if lastvar is None:
-            self.var2 = nlines
-        else:
-            self.var2 = lastvar
-        if str(nvar) == 'all':
-            self.nvar = nlines
-            self.var2 = nlines
-        else:
-            self.nvar = nvar
-        self.var2 = int(self.var2)
 
-        if n_surface == 0:
-            self.surftype = '3d volume'
-            if self.movtype in [1, 2, 3]:
-                shape = (self.n_phi_tot, self.n_theta_max, n_r_mov_tot+2)
-            else:
-                shape = (self.n_phi_tot, self.n_theta_max, self.n_r_max)
-        elif n_surface == 1:
-            self.surftype = 'r_constant'
-            shape = (self.n_phi_tot, self.n_theta_max)
-        elif n_surface == 2:
-            self.surftype = 'theta_constant'
-            if self.movtype in [1, 2, 3, 14]:  # read inner core
-                shape = (self.n_phi_tot, n_r_mov_tot+2)
-            else:
-                shape = (self.n_phi_tot, self.n_r_max)
-        elif n_surface >= 3:
-            self.surftype = 'phi_constant'
-            if self.movtype in [1, 2, 3, 14]:  # read inner core
-                shape = (self.n_theta_max, 2*(n_r_mov_tot+2))
-                self.n_theta_plot = 2*self.n_theta_max
-            elif self.movtype in [8, 9, 20, 21, 22, 23, 24, 25, 26]:
-                shape = (self.n_theta_max, n_r_mov_tot+2)
-                self.n_theta_plot = self.n_theta_max
-            elif self.movtype in [4, 5, 6, 7, 15, 16, 17, 18, 47, 54,
-                                  109, 112]:
-                shape = (self.n_theta_max, self.n_r_max, 2)
-                self.n_theta_plot = 2*self.n_theta_max
-            elif self.movtype in [10, 11, 12, 19, 92, 94, 95, 110, 111, 114,
-                                  115, 116]:
-                shape = (self.n_theta_max, self.n_r_max)
-                self.n_theta_plot = self.n_theta_max
+        # Get the number of lines
+        self._get_nlines(datadir, filename, nvar, lastar)
 
+        # Determine the shape of the data
+        self._get_data_shape(precision, allocate=False)
+
+        if self.n_surface == 3:
             if fluct or mean_field:
                 if 'B' in field:
                     prefix = 'AB_mov'
@@ -249,23 +158,28 @@ def __init__(self, file=None, step=1, lastvar=None, nvar='all', fluct=False,
 
         # If one skip the beginning, nevertheless read but do not store
         for i in range(self.var2-self.nvar):
-            n_frame, t_movieS, omega_ic, omega_ma, movieDipColat, \
-                movieDipLon, movieDipStrength, movieDipStrengthGeo = \
+            if self.version == 2:
+                n_frame, t_movieS, omega_ic, omega_ma, movieDipColat, \
+                    movieDipLon, movieDipStrength, movieDipStrengthGeo = \
                 infile.fort_read(precision)
+            else:
+                t_movieS = infile.fort_read(precision)
             for ll in range(self.n_fields):
-                dat = infile.fort_read(precision, shape=shape, order='F')
+                dat = infile.fort_read(precision, shape=self.shape, order='F')
         # then read the remaining requested nvar lines
         for k in range(self.nvar):
-            n_frame, t_movieS, omega_ic, omega_ma, movieDipColat, \
-                movieDipLon, movieDipStrength, movieDipStrengthGeo = \
-                infile.fort_read(precision)
+            if self.version == 2:
+                n_frame, t_movieS, omega_ic, omega_ma, movieDipColat, \
+                    movieDipLon, movieDipStrength, movieDipStrengthGeo = \
+                    infile.fort_read(precision)
+            else:
+                t_movieS = infile.fort_read(precision)
             for ll in range(self.n_fields):
-                dat = infile.fort_read(precision, shape=shape, order='F')
-                if n_surface == 0:
+                dat = infile.fort_read(precision, shape=self.shape, order='F')
+                if self.n_surface == 0:
                     fname = '{}{}{}_3D_{:05d}'.format(dir, os.sep, fieldName,
                                                       k+1+store_idx)
-                    if self.movtype in [1, 2, 3]:
-                        dat = dat[:, :, :self.n_r_max]
+                    dat = dat[:, :, :self.n_r_max]
                     if fluct:
                         dat -= dat.mean(axis=0)
                     elif mean_field:
@@ -274,11 +188,10 @@ def __init__(self, file=None, step=1, lastvar=None, nvar='all', fluct=False,
                             dat[ip, :, :] = tmp[:, :]
                     fname = os.path.join(datadir, fname)
                     self.scal3D2vtk(fname, dat, fieldName)
-                elif n_surface == 2:
+                elif self.n_surface == 2:
                     fname = '{}{}{}_eq_{:05d}'.format(dir, os.sep, fieldName,
                                                       k+1+store_idx)
-                    if self.movtype in [1, 2, 3, 14]:
-                        dat = dat[:, :self.n_r_max]
+                    dat = dat[:, :self.n_r_max]
                     if fluct:
                         dat -= dat.mean(axis=0)
                     elif mean_field:
@@ -288,61 +201,34 @@ def __init__(self, file=None, step=1, lastvar=None, nvar='all', fluct=False,
                         dat = tmp
                     fname = os.path.join(datadir, fname)
                     self.equat2vtk(fname, dat, fieldName)
-                elif n_surface >= 3:
+                elif self.n_surface == 3:
                     if fluct or mean_field:
                         field_m = mov_mean.data[0, k, ...]
-                    if self.movtype in [1, 2, 3, 14]:
-                        datoc0 = dat[:, :self.n_r_max]
-                        datoc1 = dat[:, self.n_r_max:2*self.n_r_max]
-                        if fluct:
-                            datoc0 -= field_m
-                            datoc1 -= field_m
-                        elif mean_field:
-                            datoc0 = field_m
-                            datoc1 = field_m
-
-                        fname = '{}{}{}_pcut{}_{:05d}'.format(dir, os.sep,
-                                                              fieldName,
-                                                              str(self.phiCut),
-                                                              k+1+store_idx)
-                        fname = os.path.join(datadir, fname)
-                        self.mer2vtk(fname, datoc0, self.phiCut, fieldName)
-                        name = str(self.phiCut+np.pi)
-                        if len(name) > 8:
-                            name = name[:8]
-                        fname = '{}{}{}_pcut{}_{:05d}'.format(dir, os.sep,
-                                                              fieldName,
-                                                              name, k+1,
-                                                              k+1+store_idx)
-                        fname = os.path.join(datadir, fname)
-                        self.mer2vtk(fname, datoc1, self.phiCut+np.pi,
-                                     fieldName)
-                    elif self.movtype in [4, 5, 6, 7, 15, 16, 17, 18, 47, 54,
-                                          91, 109, 112]:
-                        dat0 = dat[..., 0]
-                        dat1 = dat[..., 1]
-                        if fluct:
-                            dat0 -= field_m
-                            dat1 -= field_m
-                        elif mean_field:
-                            dat0 = field_m
-                            dat1 = field_m
-                        fname = '{}{}{}_pcut{}_{:05d}'.format(dir, os.sep,
-                                                              fieldName,
-                                                              str(self.phiCut),
-                                                              k+1+store_idx)
-                        fname = os.path.join(datadir, fname)
-                        self.mer2vtk(fname, dat0, self.phiCut, fieldName)
-                        name = str(self.phiCut+np.pi)
-                        if len(name) > 8:
-                            name = name[:8]
-                        fname = '{}{}{}_pcut{}_{:05d}'.format(dir, os.sep,
-                                                              fieldName,
-                                                              name,
-                                                              k+1+store_idx)
-                        fname = os.path.join(datadir, fname)
-                        self.mer2vtk(fname, dat1, self.phiCut+np.pi,
-                                     fieldName)
+                    datoc0 = dat[:, :self.n_r_max]
+                    datoc1 = dat[:, self.n_r_max:2*self.n_r_max]
+                    if fluct:
+                        datoc0 -= field_m
+                        datoc1 -= field_m
+                    elif mean_field:
+                        datoc0 = field_m
+                        datoc1 = field_m
+
+                    fname = '{}{}{}_pcut{}_{:05d}'.format(dir, os.sep,
+                                                          fieldName,
+                                                          str(self.phiCut),
+                                                          k+1+store_idx)
+                    fname = os.path.join(datadir, fname)
+                    self.mer2vtk(fname, datoc0, self.phiCut, fieldName)
+                    name = str(self.phiCut+np.pi)
+                    if len(name) > 8:
+                        name = name[:8]
+                    fname = '{}{}{}_pcut{}_{:05d}'.format(dir, os.sep,
+                                                          fieldName,
+                                                          name, k+1,
+                                                          k+1+store_idx)
+                    fname = os.path.join(datadir, fname)
+                    self.mer2vtk(fname, datoc1, self.phiCut+np.pi,
+                                 fieldName)
                 else:
                     fname = '{}{}{}_rcut{}_{:05d}'.format(dir, os.sep,
                                                           fieldName,
diff --git a/python/magic/movie3D.py b/python/magic/movie3D.py
deleted file mode 100644
index 1c14eb78..00000000
--- a/python/magic/movie3D.py
+++ /dev/null
@@ -1,215 +0,0 @@
-# -*- coding: utf-8 -*-
-import glob
-import re
-import os
-import numpy as np
-from magic.libmagic import symmetrize
-from magic.setup import buildSo
-from magic import ExtraPot
-from .npfile import *
-import sys
-if buildSo:
-    from .vtklib import *
-
-
-class Movie3D:
-    """
-    This class allows to read the 3D movie files :ref:`(B|V)_3D_.TAG<secMovieFile>`  and
-    transform them into a series of VTS files ``./vtsFiles/B3D_#.TAG`` that can be further
-    read using paraview.
-
-    >>> Movie3D(file='B_3D.TAG')
-    """
-
-    def __init__(self, file=None, step=1, lastvar=None, nvar='all', nrout=48,
-                 ratio_out=2., potExtra=False, precision=np.float32):
-        """
-        :param file: file name
-        :type file: str
-        :param nvar: the number of timesteps of the movie file we want to plot
-                     starting from the last line
-        :type nvar: int
-        :param lastvar: the number of the last timestep to be read
-        :type lastvar: int
-        :param step: the stepping between two timesteps
-        :type step: int
-        :param precision: precision of the input file, np.float32 for single
-                          precision, np.float64 for double precision
-        :type precision: str
-        :param potExtra: when set to True, potential extrapolation of the
-                         magnetic field outside the fluid domain is also
-                         computed
-        :type potExtra: bool
-        :param ratio_out: ratio of desired external radius to the CMB radius.
-                          This is is only used when potExtra=True
-        :type ratio_out: float
-        :param nrout: number of additional radial grid points to compute the
-                      potential extrapolation. This is only used when
-                      potExtra=True
-        :type nrout: int
-        """
-        if file is None:
-            dat = glob.glob('*_mov.*')
-            str1 = 'Which movie do you want ?\n'
-            for k, movie in enumerate(dat):
-                str1 += ' {}) {}\n'.format(k+1, movie)
-            index = int(input(str1))
-            try:
-                filename = dat[index-1]
-            except IndexError:
-                print('Non valid index: {} has been chosen instead'.format(dat[0]))
-                filename = dat[0]
-
-        else:
-            filename = file
-        mot = re.compile(r'.*_mov\.(.*)')
-        end = mot.findall(filename)[0]
-
-        # DETERMINE THE NUMBER OF LINES BY READING THE LOG FILE
-        logfile = open('log.{}'.format(end), 'r')
-        mot = re.compile(r'  ! WRITING MOVIE FRAME NO\s*(\d*).*')
-        for line in logfile.readlines():
-            if mot.match(line):
-                nlines = int(mot.findall(line)[0])
-        logfile.close()
-        if lastvar is None:
-            self.var2 = nlines
-        else:
-            self.var2 = lastvar
-        if str(nvar) == 'all':
-            self.nvar = nlines
-            self.var2 = nlines
-        else:
-            self.nvar = nvar
-
-        vecNames = np.r_[3]
-        scalNames = np.r_[-1]
-
-        # READ the movie file
-        infile = npfile(filename, endian='B')
-        # HEADER
-        version = infile.fort_read('|S64')
-        n_type, n_surface, const, n_fields = infile.fort_read(precision)
-        n_fields = int(n_fields)
-        n_surface = int(n_surface)
-        if n_fields == 1:
-            movtype = infile.fort_read(precision)
-            self.movtype = int(movtype)
-        else:
-            movtype = infile.fort_read(precision)
-
-        # RUN PARAMETERS
-        runid = infile.fort_read('|S64')
-        n_r_mov_tot, n_r_max, n_theta_max, n_phi_tot, minc, ra, \
-            ek, pr, prmag, radratio, tScale = infile.fort_read(precision)
-        minc = int(minc)
-        self.n_r_max = int(n_r_max)
-        self.n_theta_max = int(n_theta_max)
-        self.n_phi_tot = int(n_phi_tot)
-        n_r_mov_tot = int(n_r_mov_tot)
-
-        # GRID
-        if potExtra:
-            self.radius = np.zeros((n_r_mov_tot+1+nrout), np.float32)
-        else:
-            self.radius = np.zeros((n_r_mov_tot+2), np.float32)
-        tmp = infile.fort_read(precision)/(1.-radratio)
-        rcmb = tmp[0]
-        self.radius[2:n_r_mov_tot+2] = tmp[::-1]
-        self.theta = infile.fort_read(precision)
-        self.phi = infile.fort_read(precision)
-
-        shape = (n_r_mov_tot+2, self.n_theta_max, self.n_phi_tot)
-
-        self.time = np.zeros(self.nvar, precision)
-        if potExtra:
-            scals = np.zeros((n_r_mov_tot+1+nrout, self.n_theta_max,
-                              self.n_phi_tot*minc+1, 1), np.float32)
-        else:
-            scals = np.zeros((n_r_mov_tot+2, self.n_theta_max,
-                              self.n_phi_tot*minc+1, 1), np.float32)
-        vecr = np.zeros_like(scals)
-        vect = np.zeros_like(scals)
-        vecp = np.zeros_like(scals)
-
-        # Potential extrapolation
-        radii = self.radius
-        scals = scals.T
-        scals[0, ...] = radii
-
-        startdir = os.getcwd()
-        if not os.path.exists('vtsFiles'):
-            os.mkdir('vtsFiles')
-            os.chdir('vtsFiles')
-        else:
-            os.chdir('vtsFiles')
-        for i in range(self.var2-self.nvar):
-            n_frame, t_movieS, omega_ic, omega_ma, movieDipColat, \
-                 movieDipLon, movieDipStrength, \
-                 movieDipStrengthGeo = infile.fort_read(precision)
-            tmp = infile.fort_read(precision, shape=shape)
-            vecr[0:n_r_mov_tot+2, :, :, 0] = symmetrize(tmp, minc,
-                                                        reversed=True)
-            tmp = infile.fort_read(precision, shape=shape)
-            vect[0:n_r_mov_tot+2, :, :, 0] = symmetrize(tmp, minc,
-                                                        reversed=True)
-            tmp = infile.fort_read(precision, shape=shape)
-            vecp[0:n_r_mov_tot+2, :, :, 0] = symmetrize(tmp, minc,
-                                                        reversed=True)
-        for k in range(self.nvar):
-            n_frame, t_movieS, omega_ic, omega_ma, movieDipColat, \
-                 movieDipLon, movieDipStrength, \
-                 movieDipStrengthGeo = infile.fort_read(precision)
-            self.time[k] = t_movieS
-            if k % step == 0:
-                # print(k+self.var2-self.nvar)
-                tmp = infile.fort_read(precision, shape=shape)
-                brCMB = np.zeros((self.n_theta_max, self.n_phi_tot), np.float32)
-                brCMB = tmp[0, :, :]
-                brCMB = brCMB.T
-                if potExtra:
-                    vecr[nrout-1:nrout+n_r_mov_tot+2, :, :, 0] = \
-                        symmetrize(tmp, minc, reversed=True)
-                    tmp = infile.fort_read(precision, shape=shape)
-                    vect[nrout-1:nrout+n_r_mov_tot+2, :, :, 0] = \
-                        symmetrize(tmp, minc, reversed=True)
-                    tmp = infile.fort_read(precision, shape=shape)
-                    vecp[nrout-1:nrout+n_r_mov_tot+2, :, :, 0] = \
-                        symmetrize(tmp, minc, reversed=True)
-                else:
-                    vecr[0:n_r_mov_tot+2, :, :, 0] = \
-                        symmetrize(tmp, minc, reversed=True)
-                    tmp = infile.fort_read(precision, shape=shape)
-                    vect[0:n_r_mov_tot+2, :, :, 0] = \
-                        symmetrize(tmp, minc, reversed=True)
-                    tmp = infile.fort_read(precision, shape=shape)
-                    vecp[0:n_r_mov_tot+2, :, :, 0] = \
-                        symmetrize(tmp, minc, reversed=True)
-                filename = 'B3D_{:05d}'.format(k)
-                vecr = vecr[::-1, ...]
-                vect = vect[::-1, ...]
-                vecp = vecp[::-1, ...]
-                br = vecr.T
-                bt = vect.T
-                bp = vecp.T
-                if potExtra:
-                    pot = ExtraPot(rcmb, brCMB, minc, ratio_out=ratio_out,
-                                   nrout=nrout, cutCMB=True, deminc=True)
-                    br[0, ..., n_r_mov_tot+2:] = pot.brout
-                    bt[0, ..., n_r_mov_tot+2:] = pot.btout
-                    bp[0, ..., n_r_mov_tot+2:] = pot.bpout
-                    radii[n_r_mov_tot+2:] = pot.rout
-                else:
-                    radii = self.radius
-                vts(filename, radii, br, bt, bp, scals, scalNames, vecNames, 1)
-                print('write {}.vts'.format(filename))
-            else:  # Otherwise we read
-                vecr = infile.fort_read(precision, shape=shape)
-                vect = infile.fort_read(precision, shape=shape)
-                vecp = infile.fort_read(precision, shape=shape)
-
-        os.chdir(startdir)
-
-
-if __name__ == '__main__':
-    t1 = Movie3D(file='B_3D_mov.CJ2')
diff --git a/python/magic/plotlib.py b/python/magic/plotlib.py
index 61f6a2b5..df59b7b7 100644
--- a/python/magic/plotlib.py
+++ b/python/magic/plotlib.py
@@ -563,7 +563,7 @@ def radialContour(data, rad=0.85, label=None, proj='hammer', lon_0=0., vmax=None
     :type rasterized: bool
     :param gridColor: this is used to set the color of the grid
     :type gridColor: str
-    :param gridLineStyle: this allows to set the line style of the grid 
+    :param gridLineStyle: this allows to set the line style of the grid
                           (':', '-', '--')
     :type gridLineStyle: str
     :param gridLineWidth: this is used to tune the thickness of the lines used
diff --git a/python/magic/potential.py b/python/magic/potential.py
index 21c43443..cede6a86 100644
--- a/python/magic/potential.py
+++ b/python/magic/potential.py
@@ -31,7 +31,7 @@ def getPotEndianness(filename):
     try:
         f = npfile(filename, endian='l')
         try:  # This test is not 100% safe but in principle it should work
-            i = f.fort_read('i4')
+            i = f.fort_read(np.int32)
             if abs(i[0]) > 100:
                 raise TypeError
             endian = 'l'
@@ -39,7 +39,7 @@ def getPotEndianness(filename):
         except TypeError:
             f.close()
             f = npfile(filename, endian='B')
-            i = f.fort_read('i4')
+            i = f.fort_read(np.int32)
             if abs(i[0]) > 100:
                 raise TypeError
             endian = 'B'
@@ -47,7 +47,7 @@ def getPotEndianness(filename):
         record_marker = True
     except ValueError:
         f = open(filename, 'rb')
-        ver = np.fromfile(f, dtype='i4', count=1)[0]
+        ver = np.fromfile(f, dtype=np.int32, count=1)[0]
         if abs(ver) < 100:
             endian = 'l'
         else:
@@ -241,7 +241,7 @@ def read(self, filename, field, endian, record_marker, ic=False,
 
                 # Read header
                 self.l_max, self.n_r_max, self.n_r_ic_max, self.minc, \
-                    self.lm_max = infile.fort_read('i4')
+                    self.lm_max = infile.fort_read(np.int32)
                 self.m_max = int((self.l_max/self.minc)*self.minc)
                 self.n_m_max = int(self.m_max/self.minc+1)
                 self.ra, self.ek, self.pr, self.prmag, self.radratio, \
diff --git a/python/magic/series.py b/python/magic/series.py
index 6630967c..f83b4148 100644
--- a/python/magic/series.py
+++ b/python/magic/series.py
@@ -215,7 +215,7 @@ def plot(self):
             ax.legend(loc='best', frameon=False)
             ax.set_xlim(self.time[0], self.time[-1])
             fig.tight_layout()
-            
+
             fig1 = plt.figure()
             ax1 = fig1.add_subplot(111)
             if abs(self.lorentz_torque_ic).max() > 0.:
diff --git a/python/magic/setup.py b/python/magic/setup.py
index c08e2c35..8a9f14e5 100644
--- a/python/magic/setup.py
+++ b/python/magic/setup.py
@@ -188,7 +188,7 @@
                  '--fcompiler={}'.format(fcompiler),
                  '--compiler={}'.format(ccompiler),
                  '--opt={}'.format(f90options),
-                 omp_options, '-c', '-m', 
+                 omp_options, '-c', '-m',
                  'cylavg', omp_link,
                  'fortranLib/cyl.f90'], stderr=sp.PIPE, stdout=sp.PIPE)
 
diff --git a/python/magic/spectrum.py b/python/magic/spectrum.py
index 5aefa9fb..d1b2c171 100644
--- a/python/magic/spectrum.py
+++ b/python/magic/spectrum.py
@@ -1107,10 +1107,10 @@ def __init__(self, datadir='.', field='e_mag', iplot=False, ispec=None,
         print('Reading {}'.format(filename))
 
         if self.name == '2D_dtVrms_spec':
-            l_one = f.fort_read('i4')
+            l_one = f.fort_read(np.int32)
             if len(l_one) == 1:
                 self.version = l_one[0]
-                self.n_r_max, self.l_max = f.fort_read('i4')
+                self.n_r_max, self.l_max = f.fort_read(np.int32)
                 self.radius = f.fort_read(precision, shape=(self.n_r_max))
                 self.Cor_r_l = f.fort_read(precision,
                                            shape=(self.n_r_max, self.l_max+1))
@@ -1183,7 +1183,7 @@ def __init__(self, datadir='.', field='e_mag', iplot=False, ispec=None,
         else:
             if self.version == 'snap':
                 try:
-                    file_version = f.fort_read('i4')[0]
+                    file_version = f.fort_read(np.int32)[0]
                     if file_version > 20:
                         file_version = 0
                         # Close and reopen to rewind
@@ -1200,7 +1200,7 @@ def __init__(self, datadir='.', field='e_mag', iplot=False, ispec=None,
                 self.n_r_max, self.l_max, self.minc = out[1]
             elif self.version == 'ave':
                 try:
-                    file_version = f.fort_read('i4')[0]
+                    file_version = f.fort_read(np.int32)[0]
                     if file_version > 20:
                         file_version = 0
                         # Close and reopen to rewind