diff --git a/toyctm/toyctm.py b/toyctm/toyctm.py
index 18fde736f63d7ad1d5722baeea5d7036934ff030..0cccca590e0b167e051f9be33fa738528f144cae 100644
--- a/toyctm/toyctm.py
+++ b/toyctm/toyctm.py
@@ -322,7 +322,7 @@ class ctm:
           self.vdconc,vcfl=self.tdgrid.verticalTransportTrend(self.conc,intermdensity[0,:,:,:],self.MassFlux,self.nspec,self.vts,self.dt)
           self.maxcfl=max(self.maxcfl,vcfl)
           self.conc=self.conc+self.dt*self.vdconc
-
+          
   def strangtransport(self):
        """
        routine ctm.strangtransport performs 3d transport with Strang splitting
@@ -357,12 +357,13 @@ class ctm:
        if(self.vts.name=='null'):
        # We have no vertical tranport: things get simple
        # ENTIRE zonal time step, NO vertical transport  
+               self.tdgrid.fillMassFlux(self.MassFlux,self.t,self.dt,self.MassFluxFunc)
                self.mdconc,mcfl=self.tdgrid.meridionalTransportTrend(self.conc    ,intermdensity[0,:,:,:],self.MassFlux,self.nspec,self.ts,self.dt)
                mddensity,mcfl  =self.tdgrid.meridionalTransportTrend(intermdensity,intermdensity[0,:,:,:],self.MassFlux,1         ,self.ts,self.dt)
                self.maxcfl=max(self.maxcfl,mcfl)
                self.conc     = self.conc    +self.dt*self.mdconc
                intermdensity = intermdensity+self.dt*mddensity
-               
+               self.tdgrid.fillMassFlux(self.MassFlux,self.t+self.dt/2.,self.dt/2.,self.MassFluxFunc)   # Update mass flux for second time-step       
        else:  
        # HALF zonal time step and update concentratipn
                self.mdconc,mcfl=self.tdgrid.meridionalTransportTrend(self.conc    ,intermdensity[0,:,:,:],self.MassFlux,self.nspec,self.ts,self.dt/2.)
@@ -625,7 +626,7 @@ class ctm:
   def activeConc_ppb(self,spec):
      """ Returns numpy array of field shape self.tdgrid.active_field_shape with all the useful concentrations of species spec expressed in ppb"""
      return np.copy(self.tdgrid.extract_active(self.conc[self.specdict[spec],:,:,:]/self.density*1.e+9))
- 
+
   def returnslice(self,spec,unit='number_conc'):
      """plotslice(self,spec,fname) gives a vertical slice of species spec and the corner coordinates for later use with pcolormesh
      This routine is not portable, no design effort has been made, Quick and dirty