#ifndef AsselinFilt
      subroutine okeep(iunit)
#else
      subroutine okeep(iunit,itau,itaum1)
#endif
C
C=======================================================================
C                                                                    ===
C  This routine writes to virtual disk, slab by slab, the volume     ===
c  data stored in common block  VOLDAT. The velocity and tracers     ===
c  fields are re-written into disk after  some additional volume     ===
c  operations had been carried out, like Shapiro filtering.          ===
c                                                                    ===
c=======================================================================
c
c-----------------------------------------------------------------------
c  Define global data.
c-----------------------------------------------------------------------
c
#include <cdefs.h>
#include <param.h>
#if defined secondmean | defined AsselinFilt
# include <pconst.h>
#endif
#ifdef secondmean
# include <vertslabs.h>
#endif
#include <onedim.h>
#include <fields.h>
#include <workspa.h>
#include <voldat.h>
#ifdef AsselinFilt
# include <workspr.h>
# include <fullwd.h>
# ifdef AsselinFilt_cod
#  include <cbiopnh.h>
# endif
#endif
c
c-----------------------------------------------------------------------
c  Define local data.
c-----------------------------------------------------------------------
c
      integer i,ip,iunit,j,k,m
#ifdef secondmean
      FLOAT
     *        umz,vmz,fx
#endif
#ifdef AsselinFilt
      integer itau,itaum1
#endif
c
c-----------------------------------------------------------------------
c  Begin executable code.
c-----------------------------------------------------------------------
c
c  Write out slab data
c
#ifdef secondmean
      do 130 j=1,jmt-1
       call setvert(j)
#else
      do 130 j=1,jmt-1     
#endif
        do 106 i=1,imt
        ip=(j-1)*imt+i
#ifdef secondmean
        umz=c0
        vmz=c0        
        fx=c0
#endif
        do 100 k=1,km
          ua(i,k)=xu(ip,k)
          va(i,k)=xv(ip,k)
#ifdef secondmean
          umz=umz+dzvqz(i,k,0)*ua(i,k)
          vmz=vmz+dzvqz(i,k,0)*va(i,k)
          fx=fx+dzvqz(i,k,0)
#endif
 100    continue
#ifdef secondmean
        umz=umz/fx
        vmz=vmz/fx
        do 105 k=1,km
           ua(i,k)=ua(i,k)-umz
           va(i,k)=va(i,k)-vmz
 105    continue
#endif
 106    continue
        do 110 i=1,imt
          ip=(j-1)*imt+i
          bcon(i,1)=bkeep(ip,1)
          bcon(i,2)=bkeep(ip,2)
 110    continue
        do 120 m=1,nt
        do 120 k=1,km
        do 120 i=1,imt
          ip=(j-1)*imt+i
          ta(i,k,m)=xt(ip,k,m)
 120    continue
#ifdef AsselinFilt
        call oget(itau  ,nslab,(j-1)*nslab+1,_t)
        call oget(itaum1,nslab,(j-1)*nslab+1,_tb)
# ifndef AsselinFilt_cod
        do 1000 i=1,imt
        do 1000 k=1,km
          _u(i,k)=_u(i,k)+_nu*(_ub(i,k)+ua(i,k)-c2*_u(i,k))
          _v(i,k)=_v(i,k)+_nu*(_vb(i,k)+va(i,k)-c2*_v(i,k))
 1000  continue

        do 1200 m=1,nt
        do 1200 k=1,km
        do 1200 i=1,imt
          _t(i,k,m)=_t(i,k,m)+_nu*(_tb(i,k,m)+ta(i,k,m)-
     $          c2*_t(i,k,m))
 1200  continue
# else
        m=icod
        do 1200 k=1,km
        do 1200 i=1,imt
          _t(i,k,m)=_t(i,k,m)+_nu*(_tb(i,k,m)+ta(i,k,m)-
     $          c2*_t(i,k,m))
 1200  continue
# endif
#endif
        call oput(iunit,nslab,(j-1)*nslab+1,ta)
#ifdef AsselinFilt
        call oput(itau  ,nslab,(j-1)*nslab+1,_t)
#endif
 130  continue
      return
      end