c*********************************************
c  zvrmpxa.f  
c      draw multi-pixel planes for type-a data 
c      with subroutine zvrsuba.f
c*********************************************
c
C  P V B  V x-y x-z heimen ni natteiru
C      earg2621.f  /EARTH/TEST261  OGNBB91  FILE NAME / MHD3CD / OGNBB91
C*     use subroutine sub3ahm.f (sub3ahm.o)
c*     long flow lines and contour lines
c
c      parameter (itap=20+1*4,iip00=8,iiq00=4,ICO=3)
       parameter (iir00=150,thx=1.0,itapp=11,ICO=3)
       parameter (ori1=315.0,tim1=960.0,ori2=360,last1=8)
       PARAMETER (X0=0.0,Y0=0.0,XL0=10.0,YL0=10.0)
       PARAMETER (DXL=1.8,DYL=1.5,GTH=-60.0,TH0=60.0,DTH=10.0,JYJ=2)
c
       PARAMETER (NX= 180,NY= 60,NZ= 120,LAS1=41,IXC=30)
       PARAMETER (N1=NX+2,N2=N1*(NY+2),N3=N2*(NZ+2),n32=n3*2)
       PARAMETER (NB=8,NBB=11,N4=N3*NB,N5=N3*NBB,THH0=70.0)
       PARAMETER (NHY=NY,NHY1=NHY+1,NHY2=NHY+2)
c      PARAMETER (NXG=20,NZG=20,ICU=1,MDFU=1,THBI=-60.)
       PARAMETER (NXG=40,NZG=20,ICU=1,MDFU=1,THBI=-60.)
       PARAMETER (MX=NX+2,JYG2=NZ,JYG=JYG2+2,MY=JYG2+2,MMY=MY)
       PARAMETER (KK=MX*MY,MX2=MX,KP=0,LANK=20,KKK=N1*JYG*2,N6=KKK*8)
c
       PARAMETER (LAST= 1,IIQ0=1,XXL= 90.5,YYL= 60.50,ZZL= 60.50,MM=2)
c      PARAMETER (NXP=100,ARU=10.0,AR1=3.0,LAN1=30,LAN2=40,MOD=2)
       PARAMETER (NXP= 30,ARU=10.0,AR1=3.0,LAN1=30,LAN2=40,MOD=2)
       PARAMETER (IPEN=2,IAR=1,URMIN=0.01,BIS=0.00E-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,NP2=NP1*NHY1,NP3=NP2*NZ2,NP4=NP3*3)
c      PARAMETER (IPE1=1,MM1=N2/4,MM2=4*NZ2,MM3=MM2*NB,KK2=KK*2)
       PARAMETER (IPE1=1,MM1=N2/1,MM2=1*NZ2,MM3=MM2*NB,KK2=KK*2)
       PARAMETER (IBB=12,IBL=18,IBD=2,KBB=3,KBL=10,KBD=2)
c
       parameter (idci=2,isz=1,fct=1.00,h1=0.32,ixcc=IXC)
       parameter (tra=0.5,iviewp=100)
c       *** tra:transparency of pixel planes
       parameter (ipn=11,icc=1)
c       *** ipn:number of planes 
c       *** icc:kind of planes(1:zx 2:yz 3:xy)
c       parameter (ori1=225,ori2=360,last1=8)
c      parameter (ity=1)
c       character chrt*24/"Bzo=0nT  Bz=-5nT  t=57m"/
c       character chrt*15/"Bz=-5nT  t=57m"/
       character chrt*24/"Rotation of Incoming IMF"/       
c      character chrt*24/"Incoming uniform IMF    "/       
c      character chrt*24/"Incoming southward IMF  "/       
c      character chrt*24/"No uniform IMF          "/       
       character chg(3)*46
c      data chg(1)/"MHD Simulation for 1996 November 17 Event     "/
c      data chg(1)/"MHD Simulation for 1999 March 19 Event        "/
       data chg(1)/"MHD Simulation for 1999 October 22 Event      "/
c      data chg(1)/"MHD Simulation for 1999 March 20 Event        "/
c      data chg(2)/"B= 2.5nT     t=960m (240m)                    "/
       data chg(3)/"Density and energy of cross section           "/
c
       character chr(6)*46
       data chr(1)/"T = 22-10-1999  06:30:00                      "/
       data chr(2)/"T = 19-03-1999  23:20:00                      "/
       data chr(3)/"T = 19-03-1999  23:30:00                      "/
       data chr(4)/"T = 19-03-1999  23:40:00                      "/
       data chr(5)/"T = 19-03-1999  23:50:00                      "/
       data chr(6)/"T = 19-03-1999  24:00:00                      "/
C
       DIMENSION U(KK2),F(N3),P(N6),PP(MM1),HHX(N1),FBB(16)
       DIMENSION IZM(9),ZMIN(9),ZMAX(9),IA(10),AA(23)
c      dimension SU(KK2),SP(N6)
       dimension p1(n32),p5(n32),p0(n32),ff(n32)
       integer g1,g2
       DATA IZM/2,2,2,2,2,2,2,2,2/
C      DATA ZMIN/-3.0E-2,     0.0,      0.0,    0.0, -0.0588,
C    1           -1.0E-3, -0.8E-4,-19.23E-4,    0.0/
C      DATA ZMAX/ 3.0E-2,2.064E-6,  10.0E-3, 1.0E-3,  0.0588,
C    1            1.0E-3,  0.8E-4, 19.23E-4, 1.0E-5/
       DATA ZMIN/-3.0E-2,     0.0,      0.0,    0.0, -0.0294,
     1           -1.0E-3, -0.8E-4,-3.000E-4,    0.0/
       DATA ZMAX/ 3.0E-2,1.032E-6,   1.0E-3, 1.0E-3,  0.0294,
     1            1.0E-3,  0.8E-4, 3.000E-4, 1.0E-5/
c       REWIND 12
       open(11,file='gg220410.data',
     1         access='sequential',form='unformatted')
c
       PI=3.1415926
       XL=1.000*XL0
       YL=1.000*YL0
c *** each bottom points of simulation box 
       xb=-xxl/2+xxl*nxp/(nx2-1)
       yb=-yyl/2
       zb=-zzl/2
c ***
       call initvrml
       call viewp(iviewp)
c
       IIQ=IIQ0-1
       IIQ=0
       JJ=0
C
       DO 100 II=1,LAST
       itap=ii+10
       IKP=MM1
       IIQ=IIQ+1
C
       IP=0
       LB=1
       LL=8
cc  kokokara do 20 no roop ha start
       DO 20 M=1,LL
C
       I1=0
       IF(IKP.EQ.MM1) GO TO 174
       I2=IKP+1
       DO 173 I=I2,MM1
       I1=I1+1
       F(I1)=PP(I)
  173  CONTINUE
  174  CONTINUE
       DO 170 J=1,MM2
       IF(M.LE.4) READ(itap,END=9) PP
       IF(M.GT.4) READ(itap,END=9) PP
       DO 170 I=1,MM1
       IKP=I
       I1=I1+1
       F(I1)=PP(I)
       IF(I1.EQ.N3) GO TO 172
  170  CONTINUE
  172  CONTINUE
c
       if (m.ne.1 .and. m.ne.5) goto 210
c
c      do 500 k=1,nz2
c      do 500 j=2,ny2
c      do 500 i=1,nx2
c      i1=i+n1*(j-1)+n2*(k-1)
c      i2=i+n1*(ny2-j)+nx2*(ny2*2-2)*(nz2-k)
c      i3=i+n1*(ny2-1+j-1-1)+nx2*(ny2*2-2)*(k-1)
       do 500 k=1,nz2
       do 500 j=1,ny2-1
       do 500 i=1,nx2
       i1=i+n1*(j-1)+n2*(k-1)
       i2=i+n1*(ny2-j-1)+nx2*(ny2*2-2)*(nz2-k)
       i3=i+n1*(ny2-1+j-1)+nx2*(ny2*2-2)*(k-1)
       ff(i2)=f(i1)
       ff(i3)=f(i1)
  500  continue
c
       if (icc.eq.1) then
         npt=nx2*nz2
         idd=(2*ny2-2)/(ipn-1)
         do 510 i=1,ipn
         ipt=idd*(i-1)+1
c        write(6,*)"ipt= ",ipt
c        call readpxz(nx2,ny2,nz2,ipt,f,p0) 
         call readpxz(nx2,2*ny2-2,nz2,ipt,ff,p0)
         do 510 j=1,npt
         i1=j+npt*(i-1)
         if (m.eq.1) then
           p1(i1)=p0(j)  
         else
           p5(i1)=p0(j)
         end if
  510    continue
       else if (icc.eq.2) then
         npt=(2*ny2-2)*nz2
         idd=nx2/(ipn-1)
         do 520 i=1,ipn
         ipt=idd*(i-1)+1
c        write(6,*)"ipt= ",ipt
         call readpyz(nx2,2*ny2-2,nz2,ipt,ff,p0)
         do 520 j=1,npt
         i1=j+npt*(i-1)
         if (m.eq.1) then
           p1(i1)=p0(j)
         else
           p5(i1)=p0(j)
         end if
  520    continue
       else if (icc.eq.3) then
         npt=nx2*(2*ny2-2)
         idd=nz2/(ipn-1)
         do 530 i=1,ipn
         ipt=idd*(i-1)+1
         call readpxy(nx2,2*ny2-2,nz2,ipt,ff,p0)
         do 530 j=1,npt
         i1=j+npt*(i-1)
         if (m.eq.1) then
           p1(i1)=p0(j)
         else
           p5(i1)=p0(j)
         end if
  530    continue
       end if

  210  continue
   20  CONTINUE
       IF(IIQ.NE.IIQ0) GO TO 100
       IIQ=0
C
       do 600 i=1,ipn
       do 600 j=1,npt
       i1=j+npt*(i-1)
       p1(i1)=p5(i1)/p1(i1)
  600  continue 
c
c      **************************c
       if (icc.eq.1) then
c      write xz planes
c      ************************** 
C       1 and 2 Temperature
c      ***************************
       IM=3
       KP1=IZM(IM)
       VMIN=ZMIN(IM)
       VMAX=1.0*ZMAX(IM)
       VMAX=0.5*ZMAX(IM)
c      III=1
       do 540 i=1,ipn
       III=i
c      CALL PHAS1(III,NX,JYG2,F,U)
       CALL PHAS1(III,nx,nz,p1,U)
       CALL DIFF2(MM,nx2,nz2,U,P)
c      xb=-9.0
c      yb=-6.0
c      xl=18.0
c      yl=12.0
       ipx0=4
c      icc=1
c      zcc=0.0
       zcc=yb+(i-1)*yyl/float(ipn)
       call pixel1(nx2,nz2,xb,zb,xxl,zzl,ipx0,ico,icc,zcc,vmin,vmax,
     1             tra,u)      
  540  continue
c      *****************************
       else if (icc.eq.3) then
c      write xy planes
c      ***************************** 
C       3 and 4 Temperature
c      *****************************
       IM=3
       KP1=IZM(IM)
       VMIN=ZMIN(IM)
       VMAX=1.0*ZMAX(IM)
       VMAX=0.5*ZMAX(IM)
c      III=5
       do 550 i=1,ipn
       III=i
c      CALL PHAS1(III,NX,JYG2,F,U)
       CALL PHAS1(III,nx,ny2*2-4,p1,U)
       CALL DIFF2(MM,nx2,ny2*2-2,U,P)
c      xb=-90.0
c      yb=-60.0
c      xl=180.0
c      yl=120.0
       ipx0=4
c      icc=3
c      zcc=0.0
c      zcc=-60.0+idd*(i-1)
       zcc=zb+(i-1)*zzl/float(ipn)
       call pixel1(nx2,ny2*2-2,xb,yb,xxl,yyl,ipx0,ico,icc,zcc,
     1             vmin,vmax,tra,u)
  550  continue
c
c      *****************************
       else if (icc.eq.2) then
c      write yz planes
c      *****************************
       IM=4
       KP1=IZM(IM)
       VMIN=ZMIN(IM)
       VMAX=1.0*ZMAX(IM)
       VMAX=0.5*ZMAX(IM)
c      III=1+8
       do 560 i=1,ipn
       III=i
c      CALL PHAS1(III,NX,JYG2,F,U)
       CALL PHAS1(III,ny2*2-2-2,nz,p1,U)
       CALL DIFF2(MM,ny2*2-2,nz2,U,P)
c      xb=-60.0
c      yb=-60.0
c      xl=120.0
c      yl=120.0
       ipx0=4
c      icc=2
c      ixcc1=idd*(i-1)
c      zcc=-zcc
c      zcc=-90+idd*(i-1)
       zcc=xb+(i-1)*xxl/float(ipn)
       call pixel1(ny2*2-2,nz2,yb,zb,yyl,zzl,ipx0,ico,icc,zcc,
     1             vmin,vmax,tra,u)
  560  continue
c
       end if
  100  CONTINUE
    9  CONTINUE
       STOP
       END
c      ***************************************
       SUBROUTINE PHAS1(III,NX,NZ,F,U)
       DIMENSION F(1),U(1)
       NX2=NX+2
       NZ2=NZ+2
       M1=NX2*NZ2
       DO 10 I=1,M1
       I1=I+M1*(III-1)
       U(I)=F(I1)
   10  CONTINUE
       RETURN
       END
c      ***************************************
       SUBROUTINE DIFF2(M,MX,MY,U,V)
       DIMENSION U(1),V(1)
C
       IF(M.EQ.99) GO TO 99
       ALP=1.0/8.0
       MX1=MX-1
       MY1=MY-1
       N1=MX*MY
       DO 100 II=1,M
       DO 12 I=1,N1
       V(I)=U(I)
   12  CONTINUE
       DO 14 J=2,MY1
       DO 14 I=2,MX1
       I1=I+MX*(J-1)
       U(I1)=(1.0-4.0*ALP)*V(I1)+ALP*(V(I1-1)+V(I1+1)+V(I1-MX)+V(I1+MX))
   14  CONTINUE
       MY2=MX*(MY-2)
C
       DO 16 I=2,MX1
       I1=I
       U(I1)=(1.0-4.0*ALP)*V(I1)+ALP*(V(I1-1)+V(I1+1)+2.0*V(I1+MX))
       I1=I+MX*(MY-1)
       U(I1)=(1.0-4.0*ALP)*V(I1)+ALP*(V(I1-1)+V(I1+1)+2.0*V(I1-MX))
   16  CONTINUE
       DO 18 J=2,MY1
       I1=1+MX*(J-1)
       U(I1)=(1.0-4.0*ALP)*V(I1)+ALP*(2.0*V(I1+1)+V(I1-MX)+V(I1+MX))
       I1=MX+MX*(J-1)
       U(I1)=(1.0-4.0*ALP)*V(I1)+ALP*(2.0*V(I1-1)+V(I1-MX)+V(I1+MX))
   18  CONTINUE
C
       DO 20 I=1,MX
       I1=I+MX*MY1/2
       U(I1)=V(I1)
   20  CONTINUE
C
  100  CONTINUE
   99  CONTINUE
       RETURN
       END
c      ***************************
       subroutine readpxz(nx,ny,nz,ipt,u,v)
c      for read data of zx plane  
c      ***************************
       dimension u(1),v(1)
c
c      write(6,*)"nx=",nx,"ny=",ny,"nz=",nz
       n1=nx
       n2=nx*ny
       do 10 k=1,nz
       do 10 i=1,nx
       i1=i+n1*(ipt-1)+n2*(k-1)
       i2=i+nx*(k-1)
       v(i2)=0.5*(u(i1)+u(i1+n1))
   10  continue
       return
       end
c      ***************************
       subroutine readpyz(nx,ny,nz,ipt,u,v)
c      for read data of yz plane
c      ***************************
       dimension u(1),v(1)
c
       n1=nx
       n2=nx*ny
       do 10 j=1,ny
       do 10 k=1,nz
       i1=ipt+n1*(j-1)+n2*(k-1)
       i2=k+nz*(j-1)
       v(i2)=0.5*(u(i1)+u(i1+1))
   10  continue
       return
       end
c      ***************************
       subroutine readpxy(nx,ny,nz,ipt,u,v)
c      for read data of xy plane
c      ***************************
       dimension u(1),v(1)
c
       n1=nx
       n2=nx*ny
       do 10 i=1,nx
       do 10 j=1,ny
       i1=i+n1*(j-1)+n2*(ipt-1)
c      i2=j+ny*(i-1)
       i2=i+nx*(j-1)
       v(i2)=0.5*(u(i1)+u(i1+n2))
   10  continue
       return
       end
 
