//a41456k job  region=6144k,class=b-z]' '[a-z]'
//fort77  exec  fortclg,parm.fort='ae,opt(3)',
//     parm.go='kami=016',xynlp=on
//fort.sysin    dd  dsn=a41456a.suba.fort(sub1ag),disp=shr
//              dd   *
c      jupit.cntl(ogncc96)  from  ognc.cntl(ogncc92)
       parameter (x0=1.0,y0=1.0,xl0=2.6,yl=2.6)
       parameter (dxl=0.2,dyl=0.3,gth=-60.0,th0=60.0,dth=10.0,jyj=2)
       parameter (nx=160,ny=50,nz=100,n1=nx+2,n2=n1*(ny+2),n3=n2*(nz+2))
       parameter (nb=8,nbb=11,n4=n3*nb,n5=n3*nbb,thh0=70.0)
       parameter (nxg=30,nzg=20,icu=1,mdfu=2)
       parameter (mx=2*nx+3,jyg2=nz,jyg=jyg2+2,my=2*jyg2+3)
       parameter (kk=mx*my,mx2=mx,kp=0,lank=20,kkk=n1*jyg*2,n6=kkk*4)
       parameter (last=2,iiq0=1,xxl=80.50,yyl=25.5,zzl=50.5,mm=2)
       parameter (nxp=30,aru=10.0,ar1=3.5,lan1=30,lan2=40,mod=2)
       parameter (ipen=1,iar=0,urmin=0.01,bis=-0.0e-4,ep1=1.0e-2)
       parameter (ro01=5.0e-4,pr01=3.56e-8,vsw=0.044,vsww=0.2*vsw)
       parameter (jxg=n1*jyg,arb=2.0,arbi=3.5,n7=n3*3,nx1=nx+1)
       parameter (ny1=ny+1,nz1=nz+1,nx2=nx+2,ny2=ny+2,nz2=nz+2)
       parameter (n11=n1+1,n12=n2+1,n21=n1+n2,n22=n21+1,nxz=nx*nz)
       parameter (np1=nx2/2+nxp,np1b=nx2/2-nxp)
       parameter (np2=np1*ny2,np3=np2*nz2,np4=np3*3)
c      parameter (ibb=20,ibl=80,ibd=20,kbb=2,kbl=18,kbd=4)
       parameter (ibb=20,ibl=80,ibd=20,kbb=2,kbl=14,kbd=4)
       parameter (mm1=n2*2,mm2=51,mm3=51*nb)
c      character tl*50
c
       dimension u(kk),f(n3),p(n6),ff(n6),pp(np4)
       dimension izm(9),zmin(9),zmax(9),ia(10),aa(22)
       data izm/2,2,2,2,2,2,2,2,2/
       data zmin/-3.0e-2,    0.0,    0.0,-0.5e-3,  -0.02,
     1           -0.5e-3,-0.4e-4,-0.5e-3,    0.0/
       data zmax/ 3.0e-2, 0.5e-6, 2.0e-3, 0.5e-3,   0.02,
     1            0.5e-3, 0.4e-4, 0.5e-3, 0.5e-5/
       rewind 12
c
       ia(1)=nx
       ia(2)=jyg2
       ia(3)=nx
       ia(4)=ny
       ia(5)=nz
       ia(6)=nxp
       ia(7)=6
       ia(8)=2
       aa(1)=th0
       aa(2)=aru
       aa(3)=ar1
       aa(4)=arb
       aa(5)=xxl
       aa(6)=yyl
       aa(7)=zzl
       aa(8)=1.0
       aa(9)=0.2
       aa(11)=-1.0
       aa(12)=0.0
       aa(15)=xl0
       aa(16)=yl
       aa(17)=gth
       aa(18)=-20.0
       aa(19)=80.0
       aa(21)=30.0
       aa(20)=ep1
       xep=0.5*float(nx1-2*nxp)/float(nx1)
c
c      graphic open
c      call gopen( 1, '**  o g i n o  **' )
       call xint(  '**  o g i n o  **' )
       xll1=3.0*xl+2.0*dxl+x0
       ar2=ar1*ar1
       pi=3.1415926
       nb1=nb+1
       nxp2=nxp*2
       mxa=mx-2
       mya=my-2
       hx=xxl/float(nx1)
       hy=yyl/float(ny1)
       hz=zzl/float(nz1)
       hx2=0.5*hx
       hy2=0.5*hy
       hz2=0.5*hz
       hxg=float(nx1)/float(nxg+1)
       hzg=float(jyg2+1)/float(nzg+1)
       hxx=xxl/float(nxg+1)
       hzz=(yyl+zzl-hy2-hz2)/float(nzg+1)
       nxzg=nxg*nzg
       iiq=iiq0-1
       jj=0
c
       do 100 ii=1,last
       iiq=iiq+1
c
       jxn=nx2
       ip=0
       lb=1
       ll=8
       do 20 m=1,ll
       do 170 j1=1,mm2
       i1=1+mm1*(j1-1)
       i2=mm1*j1
       read(12,end=9) (f(i),i=i1,i2)
  170  continue
       do 20 i=1,nx2
       do 20 k=1,nz2
       i1=i+n2*(k-1)
       i2=i+n1*(k-1)+jxg*(m-1)
       p(i2)=0.5*(f(i1)+f(i1+n1))
   20  continue
       if(iiq.ne.iiq0) go to 100
       iiq=0
c
       do 34 i=1,jxg
       i1=i
       i2=i1+jxg
       i3=i2+jxg
       i4=i3+jxg
       i5=i4+jxg
       i6=i5+jxg
       i7=i6+jxg
       i8=i7+jxg
       f(i1)=p(i1)
       f(i2)=p(i2)
       f(i3)=p(i3)
       f(i4)=p(i4)
       f(i5)=p(i5)
       f(i6)=p(i6)
       f(i7)=p(i7)
       f(i8)=p(i8)
   34  continue
c
       do 70 j=1,nz2
       do 70 i=1,n1
       x=float(2*i-1-n1+2*nxp)
       z=float(2*j-1-nz2)
       r=sqrt(x*x+z*z)
       i1=i+n1*(j-1)
       i2=i1+jxg
       i3=i2+jxg
       i4=i3+jxg
       i5=i4+jxg
       i6=i5+jxg
       i7=i6+jxg
       i8=i7+jxg
       i9=i8+jxg
c      f(i9)=sqrt(f(i6)*f(i6)+f(i7)*f(i7)+f(i8)*f(i8))
       f(i9)=(x*f(i6)+z*f(i8))/r
   70  continue
c
       call xviewp(0.0,0.0,32.0,25.0)
       call plot(0.0,0.0,-3)
       call factor(2.40)
       call newpen(ipen)
       xl=1.5*xl0
       i1=ii
       vo=0.0
       call data(1.0,0.5,last,i1,nxp,nx)
c
c       1 pressure-p
       im=2
       kp1=izm(im)
       vmin=zmin(im)
       vmax=zmax(im)
       xb=x0
       yb=y0
       iii=5
       call phas2(iii,nx,jyg2,f,u)
       call diff2(mm,mx,my,u,p)
       call grap4m(mxa,mya,nxp2,xb,yb,xl,yl,hx2,hz2,ar2,u,vmin,vmax,
     1             kp1,lan1,vo)
c       2 density-ro
       im=3
       kp1=izm(im)
       vmin=zmin(im)
       vmax=zmax(im)
       xb=x0+1.0*(xl+dxl)
       iii=1
       call phas2(iii,nx,jyg2,f,u)
       call diff2(mm,mx,my,u,p)
       call grap4m(mxa,mya,nxp2,xb,yb,xl,yl,hx2,hz2,ar2,u,vmin,vmax,
     1             kp1,lan1,vo)
c       3 mod-b
       im=4
       kp1=izm(im)
       vmin=zmin(im)
       vmax=zmax(im)
       xb=x0+2.0*(xl+dxl)
       iii=9
       call phas2(iii,nx,jyg2,f,u)
       call diff2(mod,mx,my,u,p)
c      call grap4m(mxa,mya,nxp2,xb,yb,xl,yl,hx2,hz2,ar2,u,vmin,vmax,
c    1             kp1,lan1,vo)
       call grap4o(mxa,mya,nxp2,xb,yb,xl,yl,hx2,hz2,ar2,u,vmin,vmax,
     1             kp1,lank,vo,th0,dth,ipen)
c       4 vr,vz
       im=5
       kp1=izm(im)
       vmin=zmin(im)
       vmax=zmax(im)
       yb=y0+yl+dyl
       xb=x0
       do 150 j=1,nzg
       do 150 i=1,nxg
       r1=float(i)*hxg
       z1=float(j)*hzg
       ir1=ifix(r1)+1
       iz1=ifix(z1)+1
       r1=r1-float(ir1-1)
       z1=z1-float(iz1-1)
       i1=ir1+n1*(iz1-1)+jxg
       i2=i1+jxg*2
       j1=i+nxg*(j-1)
       j2=j1+nxzg
       u(j1)=(1.0-r1)*(1.0-z1)*f(i1)+r1*z1*f(i1+n1+1)
     1       +(1.0-r1)*z1*f(i1+n1)+r1*(1.0-z1)*f(i1+1)
       u(j2)=(1.0-r1)*(1.0-z1)*f(i2)+r1*z1*f(i2+n1+1)
     1       +(1.0-r1)*z1*f(i2+n1)+r1*(1.0-z1)*f(i2+1)
  150  continue
       call grap1m(nxg,nzg,nxp,xb,yb,xl,yl,hxx,hzz,ar2,u,vmin,vmax,
     1             urmin,kp1,iar,xep)
c       5  br,bz
       im=8
       kp1=izm(im)
       vmin=zmin(im)
       vmax=zmax(im)
       xb=x0+xl+dxl
       do 180 j=1,nzg
       do 180 i=1,nxg
       r1=float(i)*hxg
       z1=float(j)*hzg
       ir1=ifix(r1)+1
       iz1=ifix(z1)+1
       r1=r1-float(ir1-1)
       z1=z1-float(iz1-1)
       i1=ir1+n1*(iz1-1)+jxg*5
       i2=i1+jxg*2
       j1=i+nxg*(j-1)
       j2=j1+nxzg
       u(j1)=(1.0-r1)*(1.0-z1)*f(i1)+r1*z1*f(i1+n1+1)
     1       +(1.0-r1)*z1*f(i1+n1)+r1*(1.0-z1)*f(i1+1)
       u(j2)=(1.0-r1)*(1.0-z1)*f(i2)+r1*z1*f(i2+n1+1)
     1       +(1.0-r1)*z1*f(i2+n1)+r1*(1.0-z1)*f(i2+1)
  180  continue
       call grap1m(nxg,nzg,nxp,xb,yb,xl,yl,hxx,hzz,ar2,u,vmin,vmax,
     1             urmin,kp1,iar,xep)
c
       do 194 j=1,jyg
       do 194 i=1,nx2
       i2=i+n1*(j-1)+jxg
       i3=i2+jxg
       i4=i3+jxg
       i5=i4+jxg
       i6=i5+jxg
       i7=i6+jxg
       i8=i7+jxg
       x=sqrt(f(i6)*f(i6)+f(i7)*f(i7)+f(i8)*f(i8))
       if(x.lt.1.0e-20) x=1.0e-20
       f(i5)=(f(i2)*f(i6)+f(i3)*f(i7)+f(i4)*f(i8))/x
  194  continue
c
       do 196 j=1,ny1
       do 196 i=1,nx2
       i2=i+n1*(j-1)+jxg
       i3=i2+jxg
       i4=i3+jxg
       i5=i4+jxg
       i6=i5+jxg
       i7=i6+jxg
       i8=i7+jxg
       f(i4)=f(i3)
       f(i6)=f(i7)
  196  continue
c       6  vx
       im=5
       kp1=izm(im)
       vmin=zmin(im)
       vmax=zmax(im)
       xb=x0+2.0*(xl+dxl)
       iii=2
       call phas2(iii,nx,jyg2,f,u)
       call diff2(mm,mx,my,u,p)
       call grap4o(mxa,mya,nxp2,xb,yb,xl,yl,hx2,hz2,ar2,u,vmin,vmax,
     1             kp1,lank,vo,th0,dth,ipen)
c       7  vy, vz
       im=5
       kp1=izm(im)
       vmin=zmin(im)
       vmax=zmax(im)
       xb=x0
       yb=y0+2.0*(yl+dyl)
       iii=4
       call phas2(iii,nx,jyg2,f,u)
       call diff2(mm,mx,my,u,p)
       call grap4o(mxa,mya,nxp2,xb,yb,xl,yl,hx2,hz2,ar2,u,vmin,vmax,
     1             kp1,lank,vo,th0,dth,ipen)
c       8  bx, by
       im=8
       kp1=izm(im)
       vmin=zmin(im)
       vmax=zmax(im)
       xb=x0+xl+dxl
       iii=6
       call phas2(iii,nx,jyg2,f,u)
       call diff2(mm,mx,my,u,p)
       call grap4o(mxa,mya,nxp2,xb,yb,xl,yl,hx2,hz2,ar2,u,vmin,vmax,
     1             kp1,lank,vo,th0,dth,ipen)
c       9  bz
       im=8
       kp1=izm(im)
       vmin=zmin(im)
       vmax=zmax(im)
       xb=x0+2.0*(xl+dxl)
       iii=8
       call phas2(iii,nx,jyg2,f,u)
       call diff2(mm,mx,my,u,p)
       call grap4o(mxa,mya,nxp2,xb,yb,xl,yl,hx2,hz2,ar2,u,vmin,vmax,
     1             kp1,lank,vo,th0,dth,ipen)
c      call chart
c
  100  continue
    9  continue
c      call gclose
       call xend
c      graphic closed
       stop
       end
//go.xyout     dd   space=(cyl,(16,5))
//go.ft12f001  dd   disp=shr,dsn=a41456a.ogc7103
//