C      stecpu1:/home/g3/ogino/earthbb/ognb4.f + tunesub3.f
C      FILE NAME OGINO/EARTHB/OGNBB93  FILE NAME /POSC/OGNBB94
       parameter (itapp=20+1*4,itappp=11)
       parameter (iip00=8,iiq00=4,iir00=150,thx=1.0)
       parameter (ori1=270.0,tim1=960.0,XLIM=0.0,XLI2=10.0)
c      PARAMETER (X0=2.0,Y0=5.0,XL0=7.5,YL=7.5,ori2=360,last1=8)
       PARAMETER (X0=2.5,Y0=3.0,XL0=10.0,YL=10.0,ori2=360,last1=8)
c      PARAMETER (X0=1.2,Y0=5.0,XL0=5.5,YL=5.5)
       PARAMETER (DXL=2.0,DYL=1.0,GTH=-45.0,TH0=60.0,DTH=10.0)
       PARAMETER (NX= 180,NY= 60,NZ= 60)
       PARAMETER (N1=NX+2,N2=N1*(NY+2),N3=N2*(NZ+2))
c      PARAMETER (NB=3,NBB=11,N4=N3*NB,N5=N3*NBB,THH0=70.0)
       PARAMETER (NB=3,NBB=11,N4=N3*NB,N5=N3*NBB,THH0=60.0)
c      PARAMETER (NXG=40,NZG=20,ICU=1,MDFU=2)
       PARAMETER (NXG=40,NZG=16,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)
c      PARAMETER (LAST=5,IIQ0=1,XXL=60.15,YYL=30.15,ZZL=30.15,MM=2)
       PARAMETER (LAST=1,IIQ0=1,XXL= 90.5,YYL=30.5,ZZL=30.5,MM=2)
       PARAMETER (NXP= 30,ARU=30.0,AR1=3.0,LAN1=30,LAN2=40,MOD=2)
       PARAMETER (IPEN=1,IAR=0,URMIN=0.01,BIS=-1.5E-4,EP1=1.0E-2)
       PARAMETER (RO01=5.0E-4,PR01=3.56E-8,VSW=0.044,VSWW=0.2*VSW)
c      PARAMETER (JXG=N1*JYG,ARB=2.0,ARBI=3.0,N7=N3*3,NX1=NX+1)
       PARAMETER (JXG=N1*JYG,ARB=2.0,ARBI=2.0,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=30,IBL=120,IBD=30,KBB=2,KBL=18,KBD=4)
       PARAMETER (IBB=30,IBL=120,IBD=30,KBB=2,KBL=14,KBD=4)
       PARAMETER (M0=1,M1=N2/M0,M2=NZ2*8*M0,M3=NZ2*5*M0)
       parameter (idci=2,ity=2,isz=1,fct=1.00,H1=0.45,H2=0.35)
c      character tl1*40/'Northward turning from southward IMF'/
c      character tl2*40/'Bz=0nT   -5nT   5nT  t=161.5m (104.5m)'/
c
       character chrt*47
c      data chrt/"Incoming Northward IMF                         "/
c      data chrt/"Incoming Southward IMF                         "/
       data chrt/"3D MHD Simulation of Earth's Magnetosphere     "/
c      data chrt/"Southward Turning from Northward IMF           "/
c      data chrt/"Northward Turning from Southward IMF           "/
c
       character chg(3)*47,chh(20)*49
       data chg(1)/"Incoming Southward IMF  Bz = -5nT              "/
       data chg(2)/"Bz= 5nT (-5nT)    t=450m (500m)                "/
c      data chg(3)/"B= 5nT      t=450m (500m)                      "/
c      data chg(1)/"B=-5nT      t=450m (500m)                      "/
c      data chg(1)/"Bz=-18.2nT  Nsw= 46/cc  Vsw=412km/s  t=240.00m "/
c      data chg(1)/"Bz=-18.2nT  Nsw= 92/cc  Vsw=412km/s  t=280.99m "/
c      data chg(1)/"Bz=-18.2nT  Nsw= 46/cc  Vsw=412km/s  t=322.00m "/
c      data chg(1)/"Bz=-18.2nT  Nsw=46.3/cc  Vsw=412km/s  t=153.75m"/
c      data chg(1)/"Bz=-18.2nT  Nsw=46.3/cc  Vsw=824km/s  t=146.25m"/
c      data chg(1)/"Bz=-18.2nT  Nsw= 46/cc  Vsw=412km/s  t=270.00m "/
c      data chg(1)/"B=-5nT      t=455m                             "/
c      data chh(1)/"Bz=-18.2nT  Nsw= 46/cc  Vsw=412km/s  t=302m ( 2m)"/
       data chh(1)/"Bz=  0.0nT  Nsw= 5/cc  Vsw=300km/s  t=180m       "/
       data chh(2)/"Bz= -5.0nT  Nsw= 5/cc  Vsw=300km/s  t=180m       "/
       data chh(3)/"Bz= -5.0nT  Nsw= 5/cc  Vsw=360km/s  t=210m       "/
       data chh(4)/"Bz= -5.0nT  Nsw= 5/cc  Vsw=420km/s  t=240m       "/
       data chh(5)/"Bz= -5.0nT  Nsw= 5/cc  Vsw=480km/s  t=270m       "/
       data chh(6)/"Bz= -5.0nT  Nsw= 5/cc  Vsw=540km/s  t=300m       "/
       data chh(7)/"Bz= -5.0nT  Nsw= 5/cc  Vsw=600km/s  t=330m       "/
       data chh(8)/"Bz= -5.0nT  Nsw= 5/cc  Vsw=660km/s  t=360m       "/
       data chh(9)/"Bz= -5.0nT  Nsw= 5/cc  Vsw=270km/s  t=210m       "/
       data chh(10)/"Bz= -5.0nT  Nsw= 5/cc  Vsw=240km/s  t=240m       "/
       data chh(11)/"Bz= -5.0nT  Nsw= 5/cc  Vsw=210km/s  t=270m       "/
       data chh(12)/"Bz= -5.0nT  Nsw= 5/cc  Vsw=180km/s  t=300m       "/
       data chh(13)/"Bz= -5.0nT  Nsw= 5/cc  Vsw=150km/s  t=330m       "/
       data chh(14)/"Bz= -5.0nT  Nsw= 5/cc  Vsw=120km/s  t=360m       "/
c
       DIMENSION F(N4),P(N6)
       DIMENSION IZM(9),ZMIN(9),ZMAX(9),IA(10),AA(24)
       DATA IZM/2,2,2,2,2,2,2,2,2/
       DATA ZMIN/-3.0E-2,    0.0,-0.0E-3,    0.0,  -0.02,
     &           -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,
     &            0.5E-3, 0.4E-4, 0.5E-3, 0.5E-5/
c
c
       open(11,file='earthb10.data',
     1         access='sequential',form='unformatted')
c
c
       call xyopen
       ino=0
       DO 300 JJJ=1,1
       ITAP=JJJ+11
       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
c      AA(18)=-20.0
c      AA(19)=120.0
c      AA(21)=25.0
       AA(18)=-20.0
       AA(19)=60.0
       AA(21)=25.0
       AA(20)=EP1
       AA(23)=XLIM
       AA(24)=XLI2
       XEP=0.5*FLOAT(NX1-2*NXP)/FLOAT(NX1)
C
C      GRAPHIC OPEN
c      CALL XYOPEN(0,0.0,0.0,33.0,33.0,0)
       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)
CC       HZZ=(ZZL-HY2-HZ2)/FLOAT(NZG+1)
       NXZG=NXG*NZG
C
       IA(1)=NXG
       IA(2)=NZG
       AA(1)=THH0
       AA(4)=ARBI
       AD=2.0
       TH1=GTH*PI/180.0
       XC1=COS(TH1)
       XS1=SIN(TH1)
       YLL1=YL+DYL
       XL=XL0*0.5*(AA(19)-AA(18))/AA(21)
       AA(15)=XL
       AA(16)=YL
       xla=xl*(-aa(18))/(aa(19)-aa(18))
       yla=yl*0.5
       xlb=yla*xc1
       ylb=-yla*xs1
       IIQ=IIQ0-1
       JJ=0
C
       DO 1003 II=1,LAST
c      ITAP1=10+II
       ITAP1=11
       IIQ=IIQ+1
       DO 724 I1=1,M2
       I2=M1*(I1-M3-1)+1
       I3=M1*(I1-M3)
       IF(I1.LE.M3) READ(ITAP1) (F(I),I=1,M1)
       IF(I1.GT.M3) READ(ITAP1) (F(I),I=I2,I3)
  724  CONTINUE
cc       do 555 i=m1*2*nz2*2+1,m1*2*nz2*2.5
cc       F(i)=-F(i)
cc  555  continue
       IF(IIQ.NE.IIQ0) GO TO 1003
       IIQ=0
       IF(II.LE.0) GO TO 1003
C
       X1=2.0
C      CALL BOUND3(NX,NY,NZ,NXP,MDFU,NB,BIS,RO01,PR01,VSW,
C    &             X1,HX,HY,HZ,ARU,AR2,F,P)
C
c      CALL PLOTS(NAME,16.2)
c      CALL FACTOR(1.50)
c      CALL PLOT(0.0,0.0,-3)
c
       ino=ino+1
       call plots(idci,ity,isz,ino)
       call factor(fct)
       CALL NEWPEN(IPEN)
       I1=II
c      CALL DATA(0.5,0.5,LAST,I1,NXP,NX)
C
C       3  B  BAGNETIC FIELD  BX,BY,BZ
       IA(1)=NXG/2
       IA(2)=NZG/2
       IA(1)=NZG
       IA(2)=NZG
       IM=3
       KP1=IZM(IM)
       VMIN=ZMIN(IM)
       VMAX=ZMAX(IM)
       X1=X0-0.5*YL*XC1
       AA(13)=AMAX1(X0,X1)
c      AA(14)=Y0-0.05*YL
       AA(14)=Y0+0.00*YL+0.00*DYL
       IA(7)=6
       write(6,*) ino,' No. 01'
cc   polar cap start B Vector under figure
c      CALL AINTE1(IA,AA,F,P)
       xb=x0
       yb=aa(14)
c      call symbol(xb-h2*2.0,yb+yla-h2*0.5,h2,'X',0.0,1)
c      call symbol(xb+xla+xlb-h2*0.5,yb+yla-ylb-h2*1.5,h2,'Y',0.0,1)
c      call symbol(xb+xla-h2*0.5,yb+yl+h2*0.5,h2,'Z',0.0,1)
C
       IA(1)=NXG
       IA(2)=NZG
       IA(2)=12
       IA(2)=16
       AA(13)=AMAX1(X0,X1)
c      AA(14)=Y0+1.30*YL+DYL
       AA(14)=Y0+0.00*YL+0.00*DYL
       IA(7)=6
       write(6,*) ino,' No. 02'
cc  from equator start B Vector upper figure
c      CALL AINTE2(IA,AA,F,P)
       CALL AINTE21(IA,AA,F,P)
       yb=aa(14)
       call symbol(xb-h2*2.0,yb+yla-h2*0.5,h2,'X',0.0,1)
       call symbol(xb+xla+xlb-h2*0.5,yb+yla-ylb-h2*1.5,h2,'Y',0.0,1)
       call symbol(xb+xla-h2*0.5,yb+yl+h2*0.5,h2,'Z',0.0,1)
       call symbol(xb-h2*4.0,yb+yla-h2*2.0,h2,'20Re',0.0,4)
       call symbol(xb+xl+h2*4.0,yb+yla-h2*2.0,h2,' -60Re',0.0,6)
       call symbol(xb+xla+xlb+h2*2.0,yb+yla-ylb-h2*1.5,h2,'25Re',0.0,4)
       call symbol(xb+xla+h2*2.0,yb+yl+h2*0.5,h2,'25Re',0.0,4)
c
c      xgb=x0+yl*0.2
       xgb=x0+yl*0.1
       xgg=xgb+h1*20.0
       ygg=y0+yl*1.10+dyl*0.5
ccccccccccccccccccccccccccccc
         ori=ori1
c 500  call kakudo(xgg,ygg+h1*3.5,ori,h1*1.0)
c      call symblc(xgg+5.7*h1,ygg+h1*3.5,h1*1.0,"\\260",0.0,5)
c      call symblc(xgg-h1*1.3,ygg+h1*3.5,h1*1.0,"\\161",0.0,5)
c
C
C       3  B  BAGNETIC FIELD  BX,BY,BZ
       IA(1)=NXG/2
       IA(2)=NZG/2
       IA(1)=NZG
       IA(2)=NZG
       IM=3
       KP1=IZM(IM)
       VMIN=ZMIN(IM)
       VMAX=ZMAX(IM)
       X1=X0-0.5*YL*XC1
       AA(13)=AMAX1(X0,X1)
c      AA(14)=Y0-0.05*YL
       AA(14)=Y0+0.00*YL+0.00*DYL
       IA(7)=6
       write(6,*) ino,' No. 03'
cc   polar cap start B Vector under figure
       CALL AINTE1(IA,AA,F,P)
       xb=x0
       yb=aa(14)
       call symbol(xb-h2*2.0,yb+yla-h2*0.5,h2,'X',0.0,1)
       call symbol(xb+xla+xlb-h2*0.5,yb+yla-ylb-h2*1.5,h2,'Y',0.0,1)
       call symbol(xb+xla-h2*0.5,yb+yl+h2*0.5,h2,'Z',0.0,1)
C
c
       hh2=h1*6
c
       tim=tim1
       xg1=xgb+hh2+h1*2
c      yg1=ygg+h1*3.5
       yg1=ygg+h1*3.0
c      call timesg(xg1,yg1,tim,h1)
c      call symbol(xg1-h1*2,yg1,h1*1.0," t=",0.0,3)
c      call symbol(xg1+h1*5,yg1,h1*1.0,"min",0.0,3)
c      call symbol(xgb,yg1,h1*1.0,chg,0.0,7)
c      call symbol(xgb,yg1+h1*2.5,h1*1.4,chrt,0.0,24)
c      call symbol(xgb,yg1,h1*1.0,chg,0.0,49)
       call symbol(xgb,yg1,h1*1.0,chh(ii),0.0,49)
       call symbol(xgb,yg1+h1*2.5,h1*1.4,chrt,0.0,47)
c      call symbol(xb,yg1,h1,tl2,0.0,40)
c      call symblc(xb+h1*5.4,yg1+h1*0.2,h1*0.7,'\256',0.0,4)
c      call symblc(xb+h1*10.2,yg1+h1*0.2,h1*0.7,'\256',0.0,4)
c      call symblc(xb+h1*09.7,yg1+h1*0.2,h1*0.7,'\256',0.0,4)
c      call symbol(xb,yg1+h1*1.0,h1,tl1,0.0,40)
       CALL PLOTE
 1003  CONTINUE
  300  CONTINUE
C 999  CONTINUE
c      CALL XYCLOS
       STOP
       END
       subroutine kakudo(xg,yg,ori,h1)
       character isymb*20
       write(isymb,40) ori
   40  format("=",f6.1)
       call symbol(xg,yg,h1,isymb,0.0,8)
c       write(6,*) "(",isymb,") s"
       return
       end
c
       subroutine timesg(xg,yg,ori,h1)
       character isymb*20
       write(isymb,40) ori
   40  format(" ",f6.1)
       call symbol(xg,yg,h1,isymb,0.0,7)
c       write(6,*) "(",isymb,") s"
       return
       end