! ifx -o dublxxx dubl.f where.f90

      PROGRAM dubl ,13
      INCLUDE 'box.size'
!      parameter ia=65,ja=65,ka=33 !1998 qj
      parameter ia=181,ja=181,ka=91 !2005 resolution
!      parameter ia=361,ja=361,ka=181 !2007 resolution
      parameter ib=ia-1,jb=ja-1,kb=ka-1 !for ppa
      character fname*15,case*4,aname*15,numstr*10
      DIMENSION ap(100),u(in,jn,kn),v(in,jn,kn),w(in,jn,kn),
     # b(in,jn,kn),pp(in,jn,kn),x(in,jn,kn),y(in,jn,kn),z(in,jn,kn),
     # p(im,jm,km),vnu(kn),
     # dwdt(in,jn,kn,3),dudt(in,jn,kn,3),dvdt(in,jn,kn,3)
      dimension ppa(ia,ja,ka),ua(ia,ja,ka),va(ia,ja,ka),wa(ia,ja,ka),
     #          ba(ia,ja,ka),xa(ia,ja,ka),ya(ia,ja,ka),za(ia,ja,ka),
     # pa(ib,jb,kb),vnua(ka),dudta(ia,ja,ka,3),dvdta(ia,ja,ka,3),
     # dwdta(ia,ja,ka,3) 
      character casen*4,caser*4
      equivalence (ap(95),caser)
      integer nc
      integer iargc
c      stuffit=.true.
      stuffit=.false.
      fname='cased000000.sav'
      aname='caseE000000.sav'
      if (iargc() .ne. 3) then
          write (*,*) ' enter ''case'' ,num, ''newcase'' '
          read (*,*) case,num,casen
      else
         call getarg(1, case)
         call getarg(2, numstr)
         call getarg(3, casen )
         read(numstr,'(i)') num
      endif

      write(fname(1:4),'(a4)') case
      write(fname(6:11),'(i6.6)')num 
      write(aname(1:4),'(a4)') casen
      write(aname(6:11),'(i6.6)')num 
      open (unit=69,file=fname,STATUS= 'old',
     #form='unformatted')
      open (unit=31,file=aname,
     # status='unknown',form='unformatted')
      read(69) ap,x,y,z,b,u,v,w,pp,p,vnu,dudt,dvdt,dwdt
c  read (69) ap,x,y,z,t,u,v,w,p
      print *, ap(98),ap(99),ap(100)
      call interpolate(x,xa,in,jn,kn,ia,ja,ka)
      call interpolate(y,ya,in,jn,kn,ia,ja,ka)
      call interpolate(z,za,in,jn,kn,ia,ja,ka)
      call interpolate(b,ba,in,jn,kn,ia,ja,ka)
      call interpolate(u,ua,in,jn,kn,ia,ja,ka)
      call interpolate(v,va,in,jn,kn,ia,ja,ka)
      call interpolate(w,wa,in,jn,kn,ia,ja,ka)
      call interpolate(pp,ppa,in,jn,kn,ia,ja,ka)
      call interpolate(p,pa,im,jm,km,ib,jb,kb)
      do n=1,3
      call interpolate(dudt(1,1,1,n),dudta(1,1,1,n),in,jn,kn,ia,ja,ka)
      call interpolate(dvdt(1,1,1,n),dvdta(1,1,1,n),in,jn,kn,ia,ja,ka)
      call interpolate(dwdt(1,1,1,n),dwdta(1,1,1,n),in,jn,kn,ia,ja,ka)
      end do
      call interpolatz(vnu,vnua,kn,ka)
      caser=casen
      ap(98)=ia
      ap(99)=ja
      ap(100)=ka
      rewind(31)
      write(31) ap,xa,ya,za,ba,ua,va,wa,ppa,pa,vnua,
     # dudta,dvdta,dwdta
      close (unit=31)
      print *,fname,'  ',aname
      print *, ap(98),ap(99),ap(100)
      do i=1,ia
      print *, i,xa(i,j,2)
      enddo
      do k=1,ka
      print *,k,vnua(k)
      enddo

      end


      subroutine interpolatz(x,y,kn,ka) 1
      dimension y(ka),x(kn)
      do k=1,kn
      y(2*k-1)=x(k)
      enddo
      do k=1,kn-1
      kk=2*k
      y(kk)=.5*(y(kk-1)+y(kk+1))
      enddo
      end


      subroutine interpolate(x,y,in,jn,kn,ia,ja,ka) 12
      dimension x(in,jn,kn),y(ia,ja,ka)
      do i=1,in
      ii=2*i-1
      do j=1,jn
      jj=2*j-1
      do k=1,kn
      kk=2*k-1
      y(ii,jj,kk)=x(i,j,k)
      enddo
      enddo
      enddo

      do i=1,in-1
      ii=2*i
      do j=1,jn
      jj=2*j-1
      do k=1,kn
      kk=2*k-1
      y(ii,jj,kk)=.5*( y(ii-1,jj,kk)+ y(ii+1,jj,kk) )
      enddo
      enddo
      enddo

      do ii=1,ia
      do j=1,jn-1
      jj=2*j
      do k=1,kn
      kk=2*k-1
      y(ii,jj,kk)=.5*( y(ii,jj-1,kk)+ y(ii,jj+1,kk) )
      enddo
      enddo
      enddo

      do ii=1,ia
      do jj=1,ja
      do k=1,kn-1
      kk=2*k
      y(ii,jj,kk)=.5*( y(ii,jj,kk-1)+ y(ii,jj,kk+1) )
      enddo
      enddo
      enddo

      return
      end