! version 0.2 16 June 2014
! ifort -o cubit cubit.f where.f90
PROGRAM cubit
INCLUDE 'box.size'
parameter ia=129,ja=129,ka=kn ! for the first x-y interpolation
parameter iz=129,jz=129,kz=129 ! for the final z interpolation
character fname*15,case*4,aname*15,numstr*10
DIMENSION AP(100),U(IN,JN,kn),V(IN,JN,kn),w(in,jn,kn),
# T(IN,JN,kn),P(IN,JN,kn),X(IN,JN,kn),Y(IN,JN,KN),z(in,jn,kn)
dimension pa(ia,ja,ka),ua(ia,ja,ka),va(ia,ja,ka),wa(ia,ja,ka),
# ta(ia,ja,ka),xa(ia,ja,ka),ya(ia,ja,ka),za(ia,ja,ka)
dimension pq(iz,jz,kz),uq(iz,jz,kz),vq(iz,jz,kz),wq(iz,jz,kz),
# tq(iz,jz,kz),xq(iz,jz,kz),yq(iz,jz,kz),zq(iz,jz,kz)
common/par/
# rei,reit,swirl,dtmax,alim,diff,pit,rlx,gamma,slip,
# ex1,ex2,ex3,ex4,
# time,cput,steps,dt,dt1,dt2
equivalence (ap(1),rei)
real xs(in,jn),ys(in,jn)
dimension ps(in,jn),us(in,jn),vs(in,jn),ws(in,jn),ts(in,jn)
real yas(ja,ka),zas(ja,ka)
dimension pas(ja,ka),uas(ja,ka),vas(ja,ka),was(ja,ka),tas(ja,ka)
logical inside
real xlim,ylim,zlim,wx,wy,wxg,yxg,val
integer nc
integer iargc
xlim=.25
ylim=.25
zlim=.5
fname='cased000000.sav'
aname='casea000000.sav'
if (iargc() .ne. 2) then
write (*,*) ' enter ''case'' ,num'
read (*,*) case,num
else
call getarg(1, case)
call getarg(2, numstr)
read(numstr,'(i)') num
endif
write(fname(1:4),'(a4)') case
write(fname(6:11),'(i6.6)')num
write(aname(1:4),'(a4)') case
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,t,u,v,w,p
print *, ap(98),ap(99),ap(100)
c ic=(in-1)/2+1
c jc=(jn-1)/2+1
c ib=(in-1)/2-(ia-1)/2
c jb=(jn-1)/2-(ja-1)/2
c write (*,*) ia,ja,ka
c print *, ib,jb
c print *, ic,jc
do k=1,ka
xs=x(:,:,k)
ys=y(:,:,k)
ps=p(:,:,k)
us=u(:,:,k)
vs=v(:,:,k)
ws=w(:,:,k)
ts=t(:,:,k)
print *, 'x-y interpolating, level:',k,' at z=',z(1,1,k)
do i=1,ia
do j=1,ja
xp=-xlim+2*xlim*(i-1.)/(ia-1.)
yp=-ylim+2*ylim*(j-1.)/(ja-1.)
xa(i,j,k)=xp
ya(i,j,k)=yp
za(i,j,k)=z(1,1,k)
wxg=1.
wyg=1.
call locate_in_grid(xs,ys,xp,yp,wxg,wyg,
# 10,in,jn,wx,wy,inside,nc)
if (i==10.and.j==12) print *, xp,yp
if (i==10.and.j==12) print *, inside,nc,wx,wy
call interpolate_from_grid(ps,wx,wy,in,jn,val)
pa(i,j,k)=val
call interpolate_from_grid(us,wx,wy,in,jn,val)
ua(i,j,k)=val
call interpolate_from_grid(vs,wx,wy,in,jn,val)
va(i,j,k)=val
call interpolate_from_grid(ws,wx,wy,in,jn,val)
wa(i,j,k)=val
call interpolate_from_grid(ts,wx,wy,in,jn,val)
ta(i,j,k)=val
end do
end do
end do
do i=1,iz
yas=ya(i,:,:)
zas=za(i,:,:)
pas=pa(i,:,:)
uas=ua(i,:,:)
vas=va(i,:,:)
was=wa(i,:,:)
tas=ta(i,:,:)
print *, 'y-z interpolating, level:',i,' at x=',xa(i,1,1)
do j=1,jz
do k=1,kz
yp = -ylim+2*ylim*(j-1.)/(jz-1.)
zp = zlim*(k-1.)/(kz-1.)
xq(i,j,k) = xa(i,1,1)
yq(i,j,k) = yp
zq(i,j,k) = zp
wxg=1.
wyg=1.
call locate_in_grid(yas,zas,yp,zp,wxg,wyg,
# 10,ja,ka,wx,wy,inside,nc)
if (j==10.and.k==12) print *, xp,yp
if (j==10.and.k==12) print *, inside,nc,wx,wy
call interpolate_from_grid(pas,wx,wy,ia,ka,val)
pq(i,j,k)=val
call interpolate_from_grid(uas,wx,wy,ia,ka,val)
uq(i,j,k)=val
call interpolate_from_grid(vas,wx,wy,ia,ka,val)
vq(i,j,k)=val
call interpolate_from_grid(was,wx,wy,ia,ka,val)
wq(i,j,k)=val
call interpolate_from_grid(tas,wx,wy,ia,ka,val)
tq(i,j,k)=val
end do
end do
end do
ap(98)=iz
ap(99)=jz
ap(100)=kz
rewind(31)
write (31) ap,xq,yq,zq,tq,uq,vq,wq,pq
close (unit=31)
print *, ' interpolated grid size=',xq(iz/2,1,1)-xq(iz/2-1,1,1)
print *, 'minimum model grid size=',x(in/2,1,1)-x(in/2-1,1,1)
print *,fname,' ',aname
end