subroutine bwpunch(iou,ntim,ntimr,tim,dt1,time,dt)
c
c Write/read results to/from the *.res file.
c The order of the *.res file is essentially given by the
c read/writes in this subroutine, and must be compatible between
c output subroutine wpunch.f and the input subroutine wread.f
c     
c     (c) Copyright 1991 to 1995 by Michael A. Beer, William D. Dorland, 
c     and Gregory W. Hammett. ALL RIGHTS RESERVED.
c
c CAUTION: DO NOT EDIT wread.f DIRECTLY!
c instead, it is generated from wpunch.f by the following command
c (which appears in its proper form only in wpunch.f, not in wread.f):
c
c   sed -e "s/write(/read(/" -e "s/wpunch(/wread(/" wpunch.f > wread.f
c

      implicit none
      include 'itg.par'
      include 'itg.cmn'
      include 'itg_ri.cmn'

c subroutine arguments 
      integer iou  ! on entry, fortran i/o unit number of the *.res file
      integer ntim ! on entry, the index of last time point previously used.
                   ! on exit, the index of the last time point now used.
      integer ntimr ! the same as for ntim, but for a copy of phi(t) at
                    ! one point with potentially higher time resolution
                    ! for frequency calculation.
c output subroutine arguments:
      real tim   ! final time when dt is the same for all ky's.
      real dt1   ! time step when dt is the same for all ky's.
      real time(mz)  ! final time for various ky's.
      real dt(mz)  ! time step for various ky's.
      real ve0,ve1
      complex zve0,zve1,up0,up1,ut0,ut1
      
c local variables:
      integer i,id,l,m,n,n1,nt1,nt2,ntf1,ntf2
      real xxx

c For backwards compatibility with old cosine/sine coding
c (where the amplitude of the sine part is the negative of the
c amplitude of the imaginary part), reverse the sign of the imaginary 
c parts of some of the fields:

      do n=1,nz
         do i=1,nspecz
            do m=1,mz
               do l=1,lz
                  density_ri(2,l,m,n,i)=-density_ri(2,l,m,n,i)
                  u_par_ri(2,l,m,n,i)=-u_par_ri(2,l,m,n,i)
                  t_par_ri(2,l,m,n,i)=-t_par_ri(2,l,m,n,i)
                  t_perp_ri(2,l,m,n,i)=-t_perp_ri(2,l,m,n,i)
                  q_par_ri(2,l,m,n,i)=-q_par_ri(2,l,m,n,i)
                  q_perp_ri(2,l,m,n,i)=-q_perp_ri(2,l,m,n,i)
               enddo
            enddo
         enddo
      enddo
      do n=1,nz
         do m=1,mz
            do l=1,lz
               potential_ri(2,l,m,n)=-potential_ri(2,l,m,n)
            enddo
         enddo
      enddo

c defaults values for backwards compatibility:
      xxx=0.0
      tim=time(md)
      dt1=dt(md)

c start reading from or writing to the *.res file:

      write(iou)ld,ldb,kd,nmin,nmax,nd,lin,kdpass
      write(iou)nstp,ne,nfreq,ntotal,md,nspecies,ikx,iperiod
      write(iou)nparmom,nperpmom,nemom,iodd,iflr,iphi00,ifilter,iexp,
     *  igradon
      write(iou)(mmin(n),n=1,nd)
      write(iou)(mmax(n),n=1,nd)
      write(iou)tim,shr,qsf,epsn,eps,epse,nueeff,alpha,
     * etai,nuii,rmu1,xxx,xxx,tiovte,vy,vz,rmime,etae
      write(iou)dt0,dt1,x0,y0,z0,xp

      n1=nspecies

      nt1=ntim+1
      nt2=ntim+ne
      ntim=ntim+ne

      ntf1=ntimr+1
      ntf2=ntimr+ntotal
      ntimr=ntimr+ntotal

c eventually, put the following arrays in a more logical order:
c For now, keep backward compability.  Put new arrays at end:
      write(iou)((((density_ri(2,l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n1)
      write(iou)((((density_ri(1,l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n1)

      write(iou)(((potential_ri(2,l,m,n),l=1,ld),m=1,md),n=1,nd)
      write(iou)(((potential_ri(1,l,m,n),l=1,ld),m=1,md),n=1,nd)

      write(iou)((((t_par_ri(2,l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n1)
      write(iou)((((t_par_ri(1,l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n1)

      write(iou)((((u_par_ri(2,l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n1)
      write(iou)((((u_par_ri(1,l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n1)

      write(iou)((((t_perp_ri(2,l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n1)
      write(iou)((((t_perp_ri(1,l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n1)

      write(iou)((((q_perp_ri(2,l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n1)
      write(iou)((((q_perp_ri(1,l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n1)

      write(iou)((((q_par_ri(2,l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n1)
      write(iou)((((q_par_ri(1,l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n1)

c The above replace 3rd and 4th perp moments, not yet used toroidally:
c     write(iou)((((rps(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n1)
c     write(iou)((((rpc(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n1)
c     write(iou)((((sps(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n1)
c     write(iou)((((spc(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n1)

      write(iou)(((e_density_ri(1,l,m,n),l=1,kdpass),m=1,md),n=1,nd)
      write(iou)(((e_density_ri(2,l,m,n),l=1,kdpass),m=1,md),n=1,nd)

      write(iou)(((e_p_ri(1,l,m,n),l=1,kdpass),m=1,md),n=1,nd)
      write(iou)(((e_p_ri(2,l,m,n),l=1,kdpass),m=1,md),n=1,nd)

      write(iou)(((e_r_ri(1,l,m,n),l=1,kdpass),m=1,md),n=1,nd)
      write(iou)(((e_r_ri(2,l,m,n),l=1,kdpass),m=1,md),n=1,nd)

      write(iou)(((e_t_ri(1,l,m,n),l=1,kdpass),m=1,md),n=1,nd)
      write(iou)(((e_t_ri(2,l,m,n),l=1,kdpass),m=1,md),n=1,nd)

      write(iou)(((phi_ba_ri(1,l,m,n),l=1,kdpass),m=1,md),n=1,nd)
      write(iou)(((phi_ba_ri(2,l,m,n),l=1,kdpass),m=1,md),n=1,nd)

      write(iou)(((phi_bk_ri(1,l,m,n),l=1,ld),m=1,md),n=1,nd)
      write(iou)(((phi_bk_ri(2,l,m,n),l=1,ld),m=1,md),n=1,nd)

      write(iou)(((e_denk_ri(1,l,m,n),l=1,ld),m=1,md),n=1,nd)
      write(iou)(((e_denk_ri(2,l,m,n),l=1,ld),m=1,md),n=1,nd)

      write(iou)((phiav_old_ri(1,m,n),m=1,md),n=1,nd)
      write(iou)((phiav_old_ri(2,m,n),m=1,md),n=1,nd)

      write(iou)(timo(i),i=nt1,nt2)     !time history plot array 
      write(iou)((wpfx(i,n),i=nt1,nt2),n=1,n1)     !ion thermal flux
      write(iou)((fluxi(i,n),i=nt1,nt2),n=1,n1)    !ion particle flux
      write(iou)(fluxe(i),i=nt1,nt2)             !electron particle flux
      write(iou)(qfluxe(i),i=nt1,nt2)            !electron thermal flux
      write(iou)(gamx(i),i=nt1,nt2)     !total growth rate?
      write(iou)(((mgamx(m,n,i),m=1,md),n=1,nd),i=nt1,nt2)  !gamma vs. mode
c      write(iou)(wrhx(i),i=nt1,nt2)
c      write(iou)(wrhy(i),i=nt1,nt2)
c      write(iou)(wdux(i),i=nt1,nt2)
c      write(iou)(wu2x(i),i=nt1,nt2)
c      write(iou)(wp2x(i),i=nt1,nt2)
c      write(iou)(wv2x(i),i=nt1,nt2)
c      write(iou)(wlux(i),i=nt1,nt2)
c      write(iou)(wdpx(i),i=nt1,nt2)
c      write(iou)(wdvx(i),i=nt1,nt2)
      write(iou)(wenx(i),i=nt1,nt2)     ! total energy
c      write(iou)(wdfx(i),i=nt1,nt2)
      write(iou)(pcerr(i),i=nt1,nt2)    ! energy conservation
      write(iou)(((wtif(i,m,n),i=nt1,nt2),m=1,md),n=1,nd)  !energy vs. mode
      write(iou)(((wkif(i,m,n),i=nt1,nt2),m=1,md),n=1,nd)  !kin.en. vs. mode
      write(iou)(((wpif(i,m,n),i=nt1,nt2),m=1,md),n=1,nd)  !p^2 vs mode
      write(iou)(((psp(i,m,n),i=nt1,nt2),m=1,md),n=1,nd) !phi*phi vs mode
      write(iou)(((wenrk(i,m,n),i=nt1,nt2),m=1,md),n=1,nd) !energy vs. mode
      write(iou)(((dketa(i,m,n),i=nt1,nt2),m=1,md),n=1,nd) !eta drive v mode
      write(iou)(((dkpar(i,m,n),i=nt1,nt2),m=1,md),n=1,nd) !par damp. v mode
      write(iou)(((dktor(i,m,n),i=nt1,nt2),m=1,md),n=1,nd) !tor drive v mode
      write(iou)(((dktdp(i,m,n),i=nt1,nt2),m=1,md),n=1,nd) !tor damp. v mode
      write(iou)(((wkups(i,m,n),i=nt1,nt2),m=1,md),n=1,nd) !wupsi v mode
      write(iou)(((dktot(i,m,n),i=nt1,nt2),m=1,md),n=1,nd) !tot. drive v mode
      write(iou)((cx(i,n),i=nt1,nt2),n=1,4*nd)  !phi(x) not aved. in y
      write(iou)((cy(i,m),i=nt1,nt2),m=1,4*md)  !phi(y) not aved. in x
      write(iou)((cz(i,n),i=nt1,nt2),n=1,ld)    !phi(z) not aved. in x,y
      write(iou)((czc(i,n),i=nt1,nt2),n=1,ld)   !phi(0)phi(z) aved. in x,y
      write(iou)((czcn(i,n),i=nt1,nt2),n=1,ld)  !phi(z)^2 aved. in x,y
      write(iou)((czn(i,n),i=nt1,nt2),n=1,ld)   !ne(0)ne(z) aved. in x,y
      write(iou)((cznn(i,n),i=nt1,nt2),n=1,ld)  !ne(z)^2 aved. in x,y
      write(iou)(ct(i),i=nt1,nt2)               !phi(0)phi(t) aved. in x,y
      write(iou) ninterv
      write(iou)((phi00(i,n),i=nt1,nt2),n=1,nd)   !flux surf. ave. potential
      write(iou)((den00(i,n),i=nt1,nt2),n=1,nd)   !flux surf. ave. density
      write(iou)((upar00(i,n),i=nt1,nt2),n=1,nd)  !flux surf. ave. velocity
      write(iou)((tpar00(i,n),i=nt1,nt2),n=1,nd)  !flux surf. ave. t_par
      write(iou)((tperp00(i,n),i=nt1,nt2),n=1,nd) !flux surf. ave. t_perp
      write(iou)((drt(id,i),id=1,14),i=nt1,nt2)   !drive terms
      write(iou)(phirms(i),i=nt1,nt2)   ! rms phi
c      write(iou)(((wvif(i,m,n),i=nt1,nt2),m=1,md),n=1,nd)
c      write(iou)(((wplx(i,m,n),i=nt1,nt2),m=1,md),n=1,nd)
      write(iou)(wakx(i),i=nt1,nt2)     ! mode width in x
      write(iou)(waky(i),i=nt1,nt2)     ! mode width in y
c      write(iou)(therm(i),i=1,ld)
      write(iou)(wxsp(i),i=nt1,nt2)     ! mode spread
      write(iou)(utor(i),i=nt1,nt2)     ! toroidal flow
      write(iou)(upol(i),i=nt1,nt2)     ! poloidal flow
      write(iou)(((grtmx(i,m,n),i=nt1,nt2),m=1,md),n=1,nd) ! amplitude 
      write(iou)(((utim_ri(1,i,m,n),i=ntf1,ntf2),m=1,md),n=1,nd) ! Phi(t)
      write(iou)(((utim_ri(2,i,m,n),i=ntf1,ntf2),m=1,md),n=1,nd)
      write(iou)(timf(i),i=ntf1,ntf2)   ! time base for freq(t) and phi(t)
      write(iou)avgflag
      write(iou)date
      write(iou)note

c
c Write multi-species stuff at the end for backwards compatibility:
c
      write(iou)(rmass(n),n=1,nspecies)
      write(iou)(charge(n),n=1,nspecies)
      write(iou)(n_I(n),n=1,nspecies)
      write(iou)(Ln(n),n=1,nspecies)
      write(iou)(eta(n),n=1,nspecies)
      write(iou)(eta_par(n),n=1,nspecies)
      write(iou)(tau(n),n=1,nspecies)
      write(iou)(((phisq(i,m,n),i=nt1,nt2),m=1,md),n=1,nd)
      write(iou)((((fluximn(i,m,n,l),i=nt1,nt2),m=1,md),n=1,nd),
     .     l=1,nspecies)
      write(iou)((((qfluximn(i,m,n,l),i=nt1,nt2),m=1,md),n=1,nd),
     .     l=1,nspecies)
      write(iou)(((fluxemn(i,m,n),i=nt1,nt2),m=1,md),n=1,nd)
      write(iou)(((qfluxemn(i,m,n),i=nt1,nt2),m=1,md),n=1,nd)

c more new variables tacked on at the end for upward compabitility:
      write(iou) icrit,ntimf
      write(iou)(time(m),m=1,md)
      write(iou)(dt(m),m=1,md)
      write(iou)(((rkperp2(l,m,n),l=1,ld),m=1,md),n=1,nd)
      write(iou) ((rdiff(m,n),m=1,md),n=1,nd)
      write(iou) chi_int
      

c For backwards compatibility with old cosine/sine coding
c (where the amplitude of the sine part is the negative of the
c amplitude of the imaginary part), reverse the sign of the imaginary 
c parts of some of the fields:

c Change the signs back:

      do n=1,nz
         do i=1,nspecz
            do m=1,mz
               do l=1,lz
                  density_ri(2,l,m,n,i)=-density_ri(2,l,m,n,i)
                  u_par_ri(2,l,m,n,i)=-u_par_ri(2,l,m,n,i)
                  t_par_ri(2,l,m,n,i)=-t_par_ri(2,l,m,n,i)
                  t_perp_ri(2,l,m,n,i)=-t_perp_ri(2,l,m,n,i)
                  q_par_ri(2,l,m,n,i)=-q_par_ri(2,l,m,n,i)
                  q_perp_ri(2,l,m,n,i)=-q_perp_ri(2,l,m,n,i)
               enddo
            enddo
         enddo
      enddo
      do n=1,nz
         do m=1,mz
            do l=1,lz
               potential_ri(2,l,m,n)=-potential_ri(2,l,m,n)
            enddo
         enddo
      enddo

      return
      end