c
       program ognca7
c
c      file name  /zeartha2/eartha14.f
c      eartha14.f  modified leap-frog scheme
c      ogtac.cntl(ogncb77) from ognc.cntl(ogncd75)
c      3d mhd simulation of 1/4 earth's magnetosphere
c      cartesian coordinate   finite resistivity  45 degree boudary
c      1990.10.22 modified  1986.06.16  by tatsuki ogino
c      implicit real*8 (a-h,o-z)
c
       parameter (npe=2)
!hpf$  processors pe(npe)
c      parameter (nx=320,ny= 80,nz=160,nxp=100)
       parameter (nx=500,ny=100,nz=200,nxp=190)
c      parameter (nx=640,ny=160,nz=320,nxp=200)
c      parameter (nx=800,ny=200,nz=400,nxp=250)
       parameter (n1=nx+2,n2=n1*(ny+2),n3=n2*(nz+2))
       parameter (nb=8,nbb=11,n4=n3*nb,n5=n3*nbb,n6=n3*18)
c      parameter (last=7680,nil=2,lll=18,ibig=0)
c      parameter (iiq0=8,iip0= 32,iis0=7680,thx=4.00)
c      parameter (last=1280,nil=2,lll=18,ibig=0)
c      parameter (last=2560,nil=2,lll=18,ibig=0)
c      parameter (last=5120,nil=2,lll=18,ibig=0)
c      parameter (last=512,nil=2,lll=18,ibig=0)
c      parameter (last= 32,nil=2,lll=18,ibig=0)
       parameter (last=  1,nil=2,lll=18,ibig=0)
c      parameter (last=1920,nil=2,lll=18,ibig=0)
c      parameter (last=3840,nil=2,lll=18,ibig=0)
c      parameter (last=7680,nil=2,lll=18,ibig=0)
c      parameter (iiq0=8,iip0= 32,iis0=3840,thx=4.00)
c      parameter (iiq0=8,iip0= 32,iis0=3840,thx=3.50)
c      parameter (iiq0=8,iip0= 32,iis0=3840,thx=3.00)
c      parameter (iiq0=8,iip0= 32,iis0=3840,thx=2.50)
c      parameter (iiq0=8,iip0= 32,iis0=3840,thx=2.00)
c      parameter (iiq0=8,iip0= 32,iis0=1280,thx=2.00)
c      parameter (iiq0=8,iip0= 32,iis0=2560,thx=1.00)
c      parameter (iiq0=8,iip0= 32,iis0=5120,thx=0.50)
c      parameter (iiq0=8,iip0= 32,iis0=512,thx=0.50)
c      parameter (iiq0=8,iip0= 32,iis0= 32,thx=0.50)
       parameter (iiq0=1,iip0=  1,iis0=  1,thx=0.50)
c      parameter (iiq0=8,iip0= 32,iis0=3840,thx=1.50)
c      parameter (iiq0=8,iip0= 32,iis0=1920,thx=2.00)
c      parameter (iiq0=8,iip0= 32,iis0=1920,thx=1.50)
c      parameter (iiq0=8,iip0= 32,iis0=1920,thx=1.00)
c      parameter (iiq0=8,iip0= 32,iis0=1920,thx=0.75)
c      parameter (iiq0=8,iip0= 32,iis0=1920,thx=0.50)
c      parameter (iiq0=8,iip0= 32,iis0=1920,thx=0.33333333)
c      parameter (iiq0=8,iip0= 32,iis0=1920,thx=0.25)
c      parameter (iiq0=8,iip0= 32,iis0=1920,thx=0.16666667)
c      parameter (rmuu=0.0020,eud0=0.0000,ro01=5.0e-4,pr01=3.56e-8)
       parameter (rmuu=0.0020,eud0=0.0000,ro01=5.0e-4,pr01=3.56e-8)
c      parameter (rmuu=0.0010,eud0=0.0000,ro01=5.0e-4,pr01=3.56e-8)
c      parameter (rmuu=0.0010,eud0=0.0010,ro01=5.0e-4,pr01=3.56e-8)
       parameter (aruu=30.0,pmu0=1.0,dfu0=1.0)
c      parameter (eatt=0.0010,rrat=0.1,rra1=0.1)
c      parameter (eatt=0.0020,rrat=0.2,rra1=0.2)
       parameter (eatt=0.0020,rrat=0.1,rra1=0.1)
c      parameter (eatt=0.0010,rrat=0.2,rra1=0.2)
       parameter (nx1=nx+1,nx2=nx+2,ny1=ny+1,ny2=ny+2)
       parameter (nz1=nz+1,nz2=nz+2,nz3=nz+3)
       parameter (nin=1,itap=1)
       parameter (k2=n2*2,k3=n2*3,k4=n2*4,k5=n2*5,k6=n2*6)
       parameter (k7=n2*7,k8=n2*8,k9=n2*9,k10=n2*10,k11=n2*11)
c      parameter (orid= 45.375,ths0=0.0,kk8=k3*8)
c      parameter (orid= 90.375,ths0=0.0,kk8=k3*8)
c      parameter (orid=270.375,ths0=0.0,kk8=k3*8)
       parameter (orid=270.750,ths0=0.0,kk8=k3*8)
c      parameter (orid=270.1875,ths0=0.0,kk8=k3*8)
c      parameter (mfd=6,nfd=1441,iinob=last*12*0,ttno=128.0)
c      parameter (mfd=6,nfd=1441,iinob=last*12*2,ttno=256.0)
c      parameter (mfd=6,nfd=1441,iinob=last*24*4,ttno=256.0)
c      parameter (mfd=6,nfd=1441,iinob=last*6*19,ttno=256.0)
c      parameter (mfd=6,nfd=1441,iinob=last*6*20,ttno=512.0)
       parameter (mfd=6,nfd=1441,iinob=last*300,ttno=512.0)
c      parameter (mfd=6,nfd=1441,iinob=last*12*1,ttno=512.0)
       double precision zt0,zt1,zt2,zt3,zt
c
       dimension f(nx2,ny2,nz2,nb),u(nx2,ny2,nz2,nb),
     1           ff(nx2,ny2,nz2,nb),p(nx2,ny2,nz2,nbb),
     2           pp(nx2,ny2,nz2,3),fdd(mfd,nfd)
c      dimension gf(nx2,ny2,nz2,nb),gu(nx2,ny2,nz2,nb),
c    1           gff(nx2,ny2,nz2,nb),gpp(nx2,ny2,nz2,3),
c    2           gfdd(mfd,nfd)
       dimension gf(nx2,ny2,nz2,nb),gfdd(mfd,nfd)
       dimension hhx(n1),fbb(16),
     1           cj(10),cp(11),v(n2)
       dimension fp(nx2,ny2,nz2,nb)
c
!hpf$  distribute f(*,*,block,*) onto pe
!hpf$  distribute u(*,*,block,*) onto pe
!hpf$  distribute pp(*,*,block,*) onto pe
!hpf$  distribute ff(*,*,block,*) onto pe
!hpf$  distribute p(*,*,block,*) onto pe
!hpf$  distribute fp(*,block,*,*) onto pe
!hpf$  shadow     f(0,0,1:1,0)
!hpf$  shadow     u(0,0,1:1,0)
!hpf$  shadow     pp(0,0,1:1,0)
!hpf$  shadow     ff(0,0,1:1,0)
!hpf$  shadow     p(0,0,1:1,0)
!hpf$  asyncid id1
c      !xocl global gf,gu,gff,gpp,gfdd
c      equivalence (gf,f),(gff,ff),(gu,u),(gpp,pp)
c      equivalence (gfdd,fdd)
       common /blk/f,pp
c      m  1ro 2vr 3vo 4vz 5pr 6br 7bo 8bz 9jr 10jo 11jz
c
c      cj      r1   dr  dz  b0  rm  a2    aid ar1 gam gra
       data cj/10.0,1.2,0.5,1.0,1.9,-0.01,0.7,0.8,0.9,1.0/
c
c      cp      xl   yl   zl   ra  vo0  p0   gra  07  vsw  ibr ibz
c      data cp/161.0,41.0,51.0,4.0,6.81,2.68,1.35,7.0,.044,9.0,-0.0/
c      data cp/128.4,32.4,64.4,4.0,0.00,2.68,1.35,7.0,.044,0.0,-0.0/
c      data cp/160.5,40.5,80.5,4.0,0.00,2.68,1.35,7.0,.044,0.0,-0.0/
       data cp/250.5,50.5,100.5,4.0,0.00,2.68,1.35,7.0,.044,0.0,-0.0/
c
c
c      do 300 iii=1,5
c      do 300 iii=1,30
c      do 300 iii=1,4
       do 300 iii=1,1
       bisz=1.5
c      bisz=3.0
c      bisz=0.0
       ntap=10+iii
c      cp(11)=bisz*float(iii-2)
       cp(11)=bisz
c
       eat0=eatt
       rmu0=rmuu
       rmu0=rmuu
       aru=aruu
       pi=3.1415926
       gam=5.0/3.0
       gm1=2.0/(gam-1.0)
       gm2=0.5*(gam-1.0)
       cp(6)=10.0*(gam-1.0)*cp(7)/gam
       gra=cp(7)*1.0e-6
       po0=cp(6)*1.0e-7
       bis=cp(11)*1.0e-4
c      th=0.0+ths0*float(last)
c      th=orid*th*pi/180.0
       th=orid*pi/180.0
c      th=(orid+15.0*float(iii))*pi/180.0
       bisy=bis*cos(th)
       bisz=bis*sin(th)
       cp(11)=0.0
c      cp(11)=cp(11)
       vsw=cp(9)
       ar1=cp(4)
       ar2=ar1*ar1
       hx=cp(1)/float(nx1)
       hy=cp(2)/float(ny1)
       hz=cp(3)/float(nz1)
       t=0.5*hx*thx
       t1=0.5*t
       cj(8)=ar1
       cj(9)=gam
       cj(10)=gra
c
       dx1=0.25*t1/hx
       dy1=0.25*t1/hy
       dz1=0.25*t1/hz
       dx2=0.25*t/hx
       dy2=0.25*t/hy
       dz2=0.25*t/hz
       dx3=t1/(hx*hx)
       dy3=t1/(hy*hy)
       dz3=t1/(hz*hz)
       dx4=t/(hx*hx)
       dy4=t/(hy*hy)
       dz4=t/(hz*hz)
c
       rmu=t*rmu0*ro01
       pmu=t*rmu0*pmu0
       dfu=t*rmu0*dfu0
       eud=eud0
       ro02=rrat*ro01
       pr02=rrat*pr01
c
       write (6,12) iii,last,nx,ny,nz,n1,n2,n3,n4,n5,n6,eat0,rmu0,aru,
     1 eud,rrat,hx,hy,hz,t,t1,ro01,pr01,gra,dx2,dy2,dz2,dx4,dy4,dz4,
     2 bis,orid,th,bisy,bisz,(cp(i),i=1,11),(cj(j),j=1,10)
   12  format(1h ,5x,11i10/(1h ,5x,1p8e15.5))
c
c      if(iii.ge.2) go to 400
c
       do 20 m=1,nb
c
!hpf$  independent,new(i,j,k)
       do 22 k=1,nz2
!hpf$  on home(u(:,:,k,:)),local(u,f,i,j,m) begin
       do 2211 j=1,ny2
       do 2211 i=1,nx2
       u(i,j,k,m)=0.0
       f(i,j,k,m)=0.0
 2211  continue
!hpf$  end on
   22  continue
   20  continue
c
       do 221 i=1,n1
  221  hhx(i)=hx*float(i)
       do 222 i=1,nb
       fbb(i+nb)=1.0
  222  fbb(i)=1.0
       fbb(4)=-1.0
       fbb(6)=-1.0
       fbb(7)=-1.0
       fbb(3+nb)=-1.0
       fbb(4+nb)=-1.0
       fbb(6+nb)=-1.0
c
c      initial condition
c
c      initial condition
c
c
c
c      three dimensional cartesian model
       time=0.0
       iiq=0
       iip=0
       iis=0
       vmax=0.0
c
       call clock(zt0)
       call clock(zt1)
c
       do 100 ii=1,last
c
c
c
       iiq=iiq+1
       iip=iip+1
       iis=iis+1
c      boundary condition at nz=1 and nz=nz2
c      boundary condition at ny=1 and ny=ny2
       xx4=0.5*hx*float(2*nxp-nx1-2)
       xx3=hhx(nx1)+xx4
c
!hpf$  asynchronous(id1),nobuffer begin
       fp(1:nx2,1:ny2,1:nz2,1:nb)=f(1:nx2,1:ny2,1:nz2,1:nb)
!hpf$  end asynchronous
!hpf$  asyncwait(id1)
c
!hpf$  independent,new(i,j,k)
       do 31 j=2,ny1
!hpf$  on home(fp(:,j,:,:)),local(fp,i,k,m,x,xx,xx3,xx4,
!hpf$*         hhx) begin
       do 311 m=1,nb
       fp(2,j,1,m)=fp(1,j,2,m)
       fp(2,j,nz2,m)=fp(1,j,nz1,m)
       do 311 i=3,nx1
       x=hhx(i)+xx4
       x=1.0-x/xx3
       x=amin1(x,1.0)
       x=1.0+x/3.0
       xx=1.0-x
       fp(i,j,1,m)=x*fp(i-1,j,2,m)+xx*fp(i-2,j,3,m)
       fp(i,j,nz2,m)=x*fp(i-1,j,nz1,m)+xx*fp(i-2,j,nz,m)
  311  continue
c
       if(j.eq.2) then
       do 312 m=1,nb
       do 312 k=1,nz2
       do 312 i=2,nx1
       fp(i,1,k,m)=fp(i,2,-k+nz3,m)
  312  continue
c
       end if
c
!hpf$  end on
   31  continue
c
!hpf$  asynchronous(id1),nobuffer begin
       f(1:nx2,1:ny2,1:nz2,1:nb)=fp(1:nx2,1:ny2,1:nz2,1:nb)
!hpf$  end asynchronous
!hpf$  asyncwait(id1)
c
c
c      !xocl spread move /indo,:
c      do 322 k=1,nz2
c      do 322 i=2,nx1
c      f(i,1,k,m)=gf(i,2,-k+nz3,m)
c 322  continue
c      !xocl end spread (id)
c      !xocl movewait (id)
c
c      !hpf$  asynchronous(id1),nobuffer begin
c      f(2:nx1,1,  1,m) = f(2:nx1,2,1+nz,m)
c      f(1:nx2,1:ny2,nz2,m) = f(1:nx2,1:ny2,   2,m)
c      f(2:nx1,1,1:nz2,1:nb) = f(2:nx1,2,nz2:1:-1,1:nb)
c      !hpf$  end asynchronous
c      !hpf$  asyncwait(id1)
c
c
!hpf$  independent,new(i,j,k)
       do 32 k=1,nz2
!hpf$  on home(f(:,:,k,:)),local(f,fbb,i,j,m,x,xx,
!hpf$*         xx3,xx4,hhx) begin
       do 3211 m=1,nb
       f(2,ny2,k,m)=f(1,ny1,k,m)
       f(2,1,k,m)=f(2,2,k,m)*fbb(m+nb)
       f(2,1,k,m)=f(2,1,k,m)*fbb(m+nb)
       do 3211 i=3,nx1
       x=hhx(i)+xx4
       x=1.0-x/xx3
       x=amin1(x,1.0)
       x=1.0+x/3.0
       xx=1.0-x
       f(i,ny2,k,m)=x*f(i-1,ny1,k,m)+xx*f(i-2,ny,k,m)
       f(i,1,k,m)=f(i,2,k,m)*fbb(m+nb)
       f(i,1,k,m)=f(i,1,k,m)*fbb(m+nb)
 3211  continue
!hpf$  end on
   32  continue
c
c      boundary condition at nx=nx2
!hpf$  independent,new(i,j,k)
       do 33 k=1,nz2
!hpf$  on home(f(:,:,k,:)),local(f,i,j,m) begin
       do 3311 m=1,nb
       do 3311 j=1,ny2
       f(nx2,j,k,m)=f(nx1,j,k,m)
 3311  continue
!hpf$  end on
   33  continue
c  30  continue
c
c
c
c
c
c
c
c
       if(vmax.gt.1.0) go to 402
c
       if(iiq.ne.iiq0) go to 401
       iiq=0
c
  401  if(iip.ne.iip0) go to 403
       iip=0
c
c
  403  if(iis.ne.iis0) go to 100
       iis=0
  402  continue
       call clock(zt2)
c      if(ii.lt.1024) go to 100
c
c      write(ntap) f
       do 173 m=1,nb
       do 173 k=1,nz2
c      write(ntap) ((gf(i,j,k,m),i=1,nx2),j=1,ny2)
  173  continue
       write(6,661) ii,last,iii,ith,nin,itap,th2,th,x2,vmax
  661  format(1h ,1x,6i10/1x,1p4e15.7)
c
       call clock(zt3)
       zt1=zt1-zt0
       zt2=zt2-zt0
       zt3=zt3-zt0
       zt=zt2-zt1
       write(6,404) ii,zt0,zt1,zt2,zt3,zt
  404  format(1h , i6,1pe12.3,4(0pf12.5))
       zt1=zt3+zt0
c
       if(vmax.gt.1.0) go to 9
c
  100  continue
  300  continue
    9  continue
c      rewind 12
       stop
       end
       subroutine clock(ti)
       real*8 ti,ti1
c      call gettod(ti1)
       call fjhpf_gettod(ti1)
       ti=1.0d-6*ti1
       return
       end