c program ogncb7 c c file name earthb10.f c earthb10.f modified leap-frog scheme 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 parameter (nx=180,ny=60,nz=60) 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) parameter (last=1024,nil=2,lll=18,ibig=0) c parameter (last= 64,nil=2,lll=18,ibig=0) c parameter (iiq0=8,iip0= 32,iis0= 512,thx=4.00) parameter (iiq0=8,iip0= 32,iis0=1024,thx=4.00) c parameter (iiq0=8,iip0= 32,iis0= 32,thx=4.00) parameter (rmuu=0.0010,eud0=0.00,ro01=5.0e-4,pr01=3.56e-8) parameter (nxp= 30,aruu=30.0,pmu0=1.0,dfu0=1.0) parameter (eatt=0.0020,rrat=0.1,rra1=0.1) parameter (n11=n1+1,n12=n2+1,n21=n1+n2,n22=n21+1) parameter (nx1=nx+1,nx2=nx+2,ny1=ny+1,ny2=ny+2) parameter (nz1=nz+1,nz2=nz+2,nxz=nx*nz) parameter (nin=1,itap=2,mn1=n1*ny1,mn2=n2*nz1) parameter (m2=n3*2,m3=n3*3,m4=n3*4,m5=n3*5,m6=n3*6,m7=n3*7) parameter (m8=n3*8,m9=n3*9,m10=n3*10,m11=n3*11,m12=n3*12) 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) parameter (kk2=k3*2,kk3=k3*3,kk4=k3*4,kk5=k3*5) parameter (kk6=k3*6,kk7=k3*7,kk8=k3*8) parameter (np1=nx2/2+nxp,np1b=nx2/2-nxp) parameter (np2=np1*ny1,np3=np2*nz1,np4=np3*3) parameter (np11=np1+1,np12=np2+1,np21=np1+np2,np22=np21+1) parameter (mm0=1,mm1=n2/mm0,mm2=mm0*nz2*nb) c dimension f(n4),ff(kk8),u(n4),p(k11),hhx(n1),fbb(16), 1 pp(np4),zmax(18),zmin(18),cj(10),cp(11),v(n2) 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/ data cp/90.50,30.5,30.5,3.5,0.00,2.68,1.35,7.0,.044,0.0,-0.0/ c c rewind 12 c open(10,file='earthb10.data', c 1 access='sequential',form='unformatted') open(11,file='earthb12.data', 1 access='sequential',form='unformatted') c bisz=1.5 do 300 iii=1,1 cp(11)=bisz*float(iii-2) 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 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 (cp(i),i=1,11),(cj(j),j=1,10) 12 format(1h ,5x,11i10/(1h ,5x,1p8e15.5)) c do 20 i=1,n4 u(i)=0.0 20 f(i)=0.0 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(7+nb)=-1.0 c c initial condition call equib8(nx,ny,nz,nxp,ro01,pr01,rra1,cj,cp,f,pp) c go to 882 c write(6,212) (f(i),i=1,n4) 212 format(1h ,2x,1p10e12.4) c 882 continue c c c initial condition c c read(10) f c if(i.eq.1) write(11) f if(nin.eq.0) go to 400 do 410 ii=1,itap do 174 j1=1,mm2 i1=1+mm1*(j1-1) i2=mm1*j1 c read(10) (f(i),i=i1,i2) 174 continue 410 continue 400 continue c c c three dimensional cartesian model time=0.0 iiq=0 iip=0 iis=0 vmax=0.0 c c do 100 ii=1,last 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 do 30 m=1,nb i3=n2*2 do 31 j=2,ny1 l1=n1*(j-1)+n3*(m-1) i1=2+l1 i2=i1+mn2 p(2)=f(i2-n2-1) v(2)=f(i1+n2)*fbb(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 i1=i+l1 i2=i1+mn2 p(i)=x*f(i2-n2-1)+xx*f(i2-i3-2) 311 v(i)=f(i1+n2)*fbb(m) do 31 i=2,nx1 i2=i+l1+mn2 f(i2)=p(i) 31 f(i2-mn2)=v(i) i3=n1*2 do 32 k=1,nz2 l1=n2*(k-1)+n3*(m-1) i1=2+l1 i2=i1+mn1 p(2)=f(i2-n1-1) v(2)=f(i1+n1)*fbb(m+nb) do 321 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 i1=i+l1 i2=i1+mn1 p(i)=x*f(i2-n1-1)+xx*f(i2-i3-2) 321 v(i)=f(i1+n1)*fbb(m+nb) do 32 i=2,nx1 i2=i+l1+mn1 f(i2)=p(i) 32 f(i2-mn1)=v(i) c boundary condition at nx=nx2 do 33 k=1,nz2 i1=n1+n2*(k-1)+n3*(m-1) do 331 j=1,ny2 p(j)=f(i1-1) 331 i1=i1+n1 i1=n1+n2*(k-1)+n3*(m-1) do 33 j=1,ny2 f(i1)=p(j) 33 i1=i1+n1 30 continue do 34 i=1,n3 f(i+m4)=amax1(f(i+m4),pr02) f(i)=amax1(f(i),ro02) u(i+m4)=amax1(u(i+m4),pr02) u(i)=amax1(u(i),ro02) 34 continue c if(ii.ne.1.and.iiq.ne.iiq0) go to 172 do 371 m=1,nb do 371 k=1,nz1 do 371 j=1,ny1 l1=n1*(j-1)+n2*(k-1)+n3*(m-1) do 371 i=1,nx1 i1=i+l1 u(i1)=0.125*(f(i1)+f(i1+1)+f(i1+n1)+f(i1+n11) 1 +f(i1+n2)+f(i1+n12)+f(i1+n21)+f(i1+n22)) 371 continue t1=0.5*t dx1=0.25*t1/hx dy1=0.25*t1/hy dz1=0.25*t1/hz dx3=t1/(hx*hx) dy3=t1/(hy*hy) dz3=t1/(hz*hz) c if(nin.ne.0.or.ii.ne.1) go to 172 do 171 j1=1,mm2 i1=1+mm1*(j1-1) i2=mm1*j1 c write(12) (f(i),i=i1,i2) 171 continue 172 continue c c step of k=1 k=1 c c current do 3810 j=1,ny1 l1=n1*(j-1) l2=n1*(j-1) do 3810 i=1,nx1 i6=i+l1+m5 i7=i6+n3 i8=i7+n3 j9=i+l2+k8 p(j9+k2)=0.25*((f(i7+n22)+f(i7+n12)+f(i7+n11)+f(i7+1) 1 -f(i7+n21)-f(i7+n2)-f(i7+n1)-f(i7))/hx 2 -(f(i6+n22)-f(i6+n12)+f(i6+n11)-f(i6+1) 3 +f(i6+n21)-f(i6+n2)+f(i6+n1)-f(i6))/hy) p(j9+n2)=0.25*((f(i6+n22)+f(i6+n12)-f(i6+n11)-f(i6+1) 1 +f(i6+n21)+f(i6+n2)-f(i6+n1)-f(i6))/hz 2 -(f(i8+n22)+f(i8+n12)+f(i8+n11)+f(i8+1) 3 -f(i8+n21)-f(i8+n2)-f(i8+n1)-f(i8))/hx) p(j9)=0.25*((f(i8+n22)-f(i8+n12)+f(i8+n11)-f(i8+1) 1 +f(i8+n21)-f(i8+n2)+f(i8+n1)-f(i8))/hy 2 -(f(i7+n22)+f(i7+n12)-f(i7+n11)-f(i7+1) 3 +f(i7+n21)+f(i7+n2)-f(i7+n1)-f(i7))/hz) 3810 continue do 3812 j=1,ny1 l1=np1*(j-1)+np2*(k-1) l2=n1*(j-1)+k8 do 3814 i=1,np1 j9=i+np1b-1+l2 i1=i+l1 p(j9+k2)=p(j9+k2)-pp(i1+2*np3) p(j9+n2)=p(j9+n2)-pp(i1+np3) p(j9)=p(j9)-pp(i1) 3814 continue do 3812 i=2,np1b j9=np1b-i+1+l2 i1=i+l1 p(j9+k2)=p(j9+k2)+pp(i1+2*np3) p(j9+n2)=p(j9+n2)+pp(i1+np3) p(j9)=p(j9)-pp(i1) 3812 continue c do 3811 m=1,nb do 3811 j=1,ny1 l1=n1*(j-1)+n3*(m-1) l2=n1*(j-1)+n2*(m-1) l3=n1*(j-1)+k3*(m-1) do 3811 i=1,nx1 i1=i+l1 j1=i+l2 j2=i+l3 p(j1)=0.125*(f(i1)+f(i1+1)+f(i1+n1)+f(i1+n11) 1 +f(i1+n2)+f(i1+n12)+f(i1+n21)+f(i1+n22)) c u(j2)=p(j1) 3811 continue c c zeroth step xx4=0.5*hx*float(2*nxp-nx1-1) z=0.5*hz*float(2*k-2) do 5010 j=1,ny1 l1=n1*(j-1) l2=n1*(j-1) l3=l2 y=0.5*hy*float(2*j-2) do 5012 m=1,nb l4=l1+n3*(m-1) l5=n1*(m-1) do 5012 i=1,nx1 i1=i+l4 i2=i+l5 v(i2)=u(i1) 5012 continue c do 5011 i=1,nx1 x=hhx(i)+xx4 ro2=x*x+y*y+z*z ro2=amax1(ro2,ar2) ro3=ro2*sqrt(ro2) i1=i+l1 i2=i1+n3 i3=i2+n3 i4=i3+n3 i5=i4+n3 i6=i5+n3 i7=i6+n3 i8=i7+n3 j1=i+l2 j2=j1+n2 j3=j2+n2 j4=j3+n2 j5=j4+n2 j6=j5+n2 j7=j6+n2 j8=j7+n2 j9=j8+n2 j10=j9+n2 j11=j10+n2 c jj1=i+l3 c u(i1+m7)=u(i1+m7) 1 +dx1*(f(i4+n22)*f(i6+n22)-f(i2+n22)*f(i8+n22) 2 +f(i4+n12)*f(i6+n12)-f(i2+n12)*f(i8+n12) 3 +f(i4+n11)*f(i6+n11)-f(i2+n11)*f(i8+n11) 4 +f(i4+1)*f(i6+1)-f(i2+1)*f(i8+1) 5 -f(i4+n21)*f(i6+n21)+f(i2+n21)*f(i8+n21) 6 -f(i4+n2)*f(i6+n2)+f(i2+n2)*f(i8+n2) 7 -f(i4+n1)*f(i6+n1)+f(i2+n1)*f(i8+n1) 8 -f(i4)*f(i6)+f(i2)*f(i8)) 1 -dy1*(f(i3+n22)*f(i8+n22)-f(i4+n22)*f(i7+n22) 2 -f(i3+n12)*f(i8+n12)+f(i4+n12)*f(i7+n12) 3 +f(i3+n11)*f(i8+n11)-f(i4+n11)*f(i7+n11) 4 -f(i3+1)*f(i8+1)+f(i4+1)*f(i7+1) 5 +f(i3+n21)*f(i8+n21)-f(i4+n21)*f(i7+n21) 6 -f(i3+n2)*f(i8+n2)+f(i4+n2)*f(i7+n2) 7 +f(i3+n1)*f(i8+n1)-f(i4+n1)*f(i7+n1) 8 -f(i3)*f(i8)+f(i4)*f(i7)) u(i1+m6)=u(i1+m6) 1 +dz1*(f(i3+n22)*f(i8+n22)-f(i4+n22)*f(i7+n22) 2 +f(i3+n12)*f(i8+n12)-f(i4+n12)*f(i7+n12) 3 -f(i3+n11)*f(i8+n11)+f(i4+n11)*f(i7+n11) 4 -f(i3+1)*f(i8+1)+f(i4+1)*f(i7+1) 5 +f(i3+n21)*f(i8+n21)-f(i4+n21)*f(i7+n21) 6 +f(i3+n2)*f(i8+n2)-f(i4+n2)*f(i7+n2) 7 -f(i3+n1)*f(i8+n1)+f(i4+n1)*f(i7+n1) 8 -f(i3)*f(i8)+f(i4)*f(i7)) 1 -dx1*(f(i2+n22)*f(i7+n22)-f(i3+n22)*f(i6+n22) 2 +f(i2+n12)*f(i7+n12)-f(i3+n12)*f(i6+n12) 3 +f(i2+n11)*f(i7+n11)-f(i3+n11)*f(i6+n11) 4 +f(i2+1)*f(i7+1)-f(i3+1)*f(i6+1) 5 -f(i2+n21)*f(i7+n21)+f(i3+n21)*f(i6+n21) 6 -f(i2+n2)*f(i7+n2)+f(i3+n2)*f(i6+n2) 7 -f(i2+n1)*f(i7+n1)+f(i3+n1)*f(i6+n1) 8 -f(i2)*f(i7)+f(i3)*f(i6)) u(i1+m5)=u(i1+m5) 1 +dy1*(f(i2+n22)*f(i7+n22)-f(i3+n22)*f(i6+n22) 2 -f(i2+n12)*f(i7+n12)+f(i3+n12)*f(i6+n12) 3 +f(i2+n11)*f(i7+n11)-f(i3+n11)*f(i6+n11) 4 -f(i2+1)*f(i7+1)+f(i3+1)*f(i6+1) 5 +f(i2+n21)*f(i7+n21)-f(i3+n21)*f(i6+n21) 6 -f(i2+n2)*f(i7+n2)+f(i3+n2)*f(i6+n2) 7 +f(i2+n1)*f(i7+n1)-f(i3+n1)*f(i6+n1) 8 -f(i2)*f(i7)+f(i3)*f(i6)) 1 -dz1*(f(i4+n22)*f(i6+n22)-f(i2+n22)*f(i8+n22) 2 +f(i4+n12)*f(i6+n12)-f(i2+n12)*f(i8+n12) 3 -f(i4+n11)*f(i6+n11)+f(i2+n11)*f(i8+n11) 4 -f(i4+1)*f(i6+1)+f(i2+1)*f(i8+1) 5 +f(i4+n21)*f(i6+n21)-f(i2+n21)*f(i8+n21) 6 +f(i4+n2)*f(i6+n2)-f(i2+n2)*f(i8+n2) 7 -f(i4+n1)*f(i6+n1)+f(i2+n1)*f(i8+n1) 8 -f(i4)*f(i6)+f(i2)*f(i8)) c x3=dx1*(f(i2+n22)+f(i2+n12)+f(i2+n11)+f(i2+1) 1 -f(i2+n21)-f(i2+n2)-f(i2+n1)-f(i2)) 2 +dy1*(f(i3+n22)-f(i3+n12)+f(i3+n11)-f(i3+1) 3 +f(i3+n21)-f(i3+n2)+f(i3+n1)-f(i3)) 4 +dz1*(f(i4+n22)+f(i4+n12)-f(i4+n11)-f(i4+1) 5 +f(i4+n21)+f(i4+n2)-f(i4+n1)-f(i4)) u(i1+m4)=u(i1+m4)-p(j5)*gam*x3 1 -dx1*p(j2)*(f(i5+n22)+f(i5+n12)+f(i5+n11)+f(i5+1) 2 -f(i5+n21)-f(i5+n2)-f(i5+n1)-f(i5)) 3 -dy1*p(j3)*(f(i5+n22)-f(i5+n12)+f(i5+n11)-f(i5+1) 4 +f(i5+n21)-f(i5+n2)+f(i5+n1)-f(i5)) 5 -dz1*p(j4)*(f(i5+n22)+f(i5+n12)-f(i5+n11)-f(i5+1) 6 +f(i5+n21)+f(i5+n2)-f(i5+n1)-f(i5)) u(i1+m3)=u(i1+m3)+t1*((p(j9)*p(j7)-p(j10)*p(j6))/p(j1) 1 -eud*p(j4)-gra*z/ro3) 2 -dx1*p(j2)*(f(i4+n22)+f(i4+n12)+f(i4+n11)+f(i4+1) 3 -f(i4+n21)-f(i4+n2)-f(i4+n1)-f(i4)) 4 -dy1*p(j3)*(f(i4+n22)-f(i4+n12)+f(i4+n11)-f(i4+1) 5 +f(i4+n21)-f(i4+n2)+f(i4+n1)-f(i4)) 6 -dz1*p(j4)*(f(i4+n22)+f(i4+n12)-f(i4+n11)-f(i4+1) 7 +f(i4+n21)+f(i4+n2)-f(i4+n1)-f(i4)) 8 -dz1*(f(i5+n22)+f(i5+n12)-f(i5+n11)-f(i5+1) 9 +f(i5+n21)+f(i5+n2)-f(i5+n1)-f(i5))/p(j1) u(i1+m2)=u(i1+m2)+t1*((p(j11)*p(j6)-p(j9)*p(j8))/p(j1) 1 -eud*p(j3)-gra*y/ro3) 2 -dx1*p(j2)*(f(i3+n22)+f(i3+n12)+f(i3+n11)+f(i3+1) 3 -f(i3+n21)-f(i3+n2)-f(i3+n1)-f(i3)) 4 -dy1*p(j3)*(f(i3+n22)-f(i3+n12)+f(i3+n11)-f(i3+1) 5 +f(i3+n21)-f(i3+n2)+f(i3+n1)-f(i3)) 6 -dz1*p(j4)*(f(i3+n22)+f(i3+n12)-f(i3+n11)-f(i3+1) 7 +f(i3+n21)+f(i3+n2)-f(i3+n1)-f(i3)) 8 -dy1*(f(i5+n22)-f(i5+n12)+f(i5+n11)-f(i5+1) 9 +f(i5+n21)-f(i5+n2)+f(i5+n1)-f(i5))/p(j1) u(i1+n3)=u(i1+n3)+t1*((p(j10)*p(j8)-p(j11)*p(j7))/p(j1) 1 -eud*p(j2)-gra*x/ro3) 2 -dx1*p(j2)*(f(i2+n22)+f(i2+n12)+f(i2+n11)+f(i2+1) 3 -f(i2+n21)-f(i2+n2)-f(i2+n1)-f(i2)) 4 -dy1*p(j3)*(f(i2+n22)-f(i2+n12)+f(i2+n11)-f(i2+1) 5 +f(i2+n21)-f(i2+n2)+f(i2+n1)-f(i2)) 6 -dz1*p(j4)*(f(i2+n22)+f(i2+n12)-f(i2+n11)-f(i2+1) 7 +f(i2+n21)+f(i2+n2)-f(i2+n1)-f(i2)) 8 -dx1*(f(i5+n22)+f(i5+n12)+f(i5+n11)+f(i5+1) 9 -f(i5+n21)-f(i5+n2)-f(i5+n1)-f(i5))/p(j1) u(i1)=u(i1) 1 -dx1*(f(i2+n22)*f(i1+n22)+f(i2+n12)*f(i1+n12) 2 +f(i2+n11)*f(i1+n11)+f(i2+1)*f(i1+1) 3 -f(i2+n21)*f(i1+n21)-f(i2+n2)*f(i1+n2) 4 -f(i2+n1)*f(i1+n1)-f(i2)*f(i1)) 5 -dy1*(f(i3+n22)*f(i1+n22)-f(i3+n12)*f(i1+n12) 6 +f(i3+n11)*f(i1+n11)-f(i3+1)*f(i1+1) 7 +f(i3+n21)*f(i1+n21)-f(i3+n2)*f(i1+n2) 8 +f(i3+n1)*f(i1+n1)-f(i3)*f(i1)) 1 -dz1*(f(i4+n22)*f(i1+n22)+f(i4+n12)*f(i1+n12) 2 -f(i4+n11)*f(i1+n11)-f(i4+1)*f(i1+1) 3 +f(i4+n21)*f(i1+n21)+f(i4+n2)*f(i1+n2) 4 -f(i4+n1)*f(i1+n1)-f(i4)*f(i1)) 5011 continue c do 5201 i=1,nx1 x=hhx(i)+xx4 ro2=x*x+y*y+z*z i1=i i2=i1+n1 i3=i2+n1 i4=i3+n1 i5=i4+n1 i6=i5+n1 i7=i6+n1 i8=i7+n1 j1=i+l1 x2=amax1(ro2/ar2-1.0,0.0) x2=aru*x2*x2 x2=x2/(x2+1.0) u(j1+m7)=u(j1+m7)*x2+v(i8)*(1.0-x2) u(j1+m6)=u(j1+m6)*x2+v(i7)*(1.0-x2) u(j1+m5)=u(j1+m5)*x2+v(i6)*(1.0-x2) u(j1+m4)=u(j1+m4)*x2+v(i5)*(1.0-x2) u(j1+m3)=u(j1+m3)*x2+v(i4)*(1.0-x2) u(j1+m2)=u(j1+m2)*x2+v(i3)*(1.0-x2) u(j1+n3)=u(j1+n3)*x2+v(i2)*(1.0-x2) u(j1)=u(j1)*x2+v(i1)*(1.0-x2) 5201 continue c do 5202 i=1,nx1 i1=i+l1 i2=i1+n3 i3=i2+n3 i4=i3+n3 i5=i4+n3 p(i+k4)=u(i5) p(i+k3)=u(i4) p(i+k2)=u(i3) p(i+n2)=u(i2) p(i)=u(i1) 5202 continue c do 5204 i=1,nx1 i1=i+l1 x5=p(i+k4) x1=p(i) y5=amax1(x5,pr02) y1=amax1(x1,ro02) c x2=p(i+n2)*p(i+n2)+p(i+k2)*p(i+k2)+p(i+k3)*p(i+k3) c vmax=amax1(vmax,x2) x=amax1(x1,0.0)/y1 u(i1+m4)=y5 u(i1+m3)=x*p(i+k3) u(i1+m2)=x*p(i+k2) u(i1+n3)=x*p(i+n2) u(i1)=y1 5204 continue 5010 continue c do 620 m=1,nb l1=n3*(m-1) l2=k3*(m-1) do 620 i=1,n2 i1=i+l1 j1=i+l2 ff(j1)=f(i1) 620 continue c c step of k=k do 90 k=2,nz1 c c current do 38 j=1,ny1 l1=n1*(j-1)+n2*(k-1) l2=n1*(j-1) do 38 i=1,nx1 i6=i+l1+m5 i7=i6+n3 i8=i7+n3 j9=i+l2+k8 p(j9+k2)=0.25*((f(i7+n22)+f(i7+n12)+f(i7+n11)+f(i7+1) 1 -f(i7+n21)-f(i7+n2)-f(i7+n1)-f(i7))/hx 2 -(f(i6+n22)-f(i6+n12)+f(i6+n11)-f(i6+1) 3 +f(i6+n21)-f(i6+n2)+f(i6+n1)-f(i6))/hy) p(j9+n2)=0.25*((f(i6+n22)+f(i6+n12)-f(i6+n11)-f(i6+1) 1 +f(i6+n21)+f(i6+n2)-f(i6+n1)-f(i6))/hz 2 -(f(i8+n22)+f(i8+n12)+f(i8+n11)+f(i8+1) 3 -f(i8+n21)-f(i8+n2)-f(i8+n1)-f(i8))/hx) p(j9)=0.25*((f(i8+n22)-f(i8+n12)+f(i8+n11)-f(i8+1) 1 +f(i8+n21)-f(i8+n2)+f(i8+n1)-f(i8))/hy 2 -(f(i7+n22)+f(i7+n12)-f(i7+n11)-f(i7+1) 3 +f(i7+n21)+f(i7+n2)-f(i7+n1)-f(i7))/hz) 38 continue do 382 j=1,ny1 l1=np1*(j-1)+np2*(k-1) l2=n1*(j-1)+k8 do 384 i=1,np1 j9=i+np1b-1+l2 i1=i+l1 p(j9+k2)=p(j9+k2)-pp(i1+2*np3) p(j9+n2)=p(j9+n2)-pp(i1+np3) p(j9)=p(j9)-pp(i1) 384 continue do 382 i=2,np1b j9=np1b-i+1+l2 i1=i+l1 p(j9+k2)=p(j9+k2)+pp(i1+2*np3) p(j9+n2)=p(j9+n2)+pp(i1+np3) p(j9)=p(j9)-pp(i1) 382 continue c do 381 m=1,nb do 381 j=1,ny1 l1=n1*(j-1)+n2*(k-1)+n3*(m-1) l2=n1*(j-1)+n2*(m-1) l3=n1*(j-1)+n2+k3*(m-1) do 381 i=1,nx1 i1=i+l1 j1=i+l2 j2=i+l3 p(j1)=0.125*(f(i1)+f(i1+1)+f(i1+n1)+f(i1+n11) 1 +f(i1+n2)+f(i1+n12)+f(i1+n21)+f(i1+n22)) c u(j2)=p(j1) 381 continue c c first step xx4=0.5*hx*float(2*nxp-nx1-1) z=0.5*hz*float(2*k-2) do 50 j=1,ny1 l1=n1*(j-1)+n2*(k-1) l2=n1*(j-1) l3=l2+n2 y=0.5*hy*float(2*j-2) do 502 m=1,nb l4=l1+n3*(m-1) l5=n1*(m-1) do 502 i=1,nx1 i1=i+l4 i2=i+l5 v(i2)=u(i1) 502 continue c do 501 i=1,nx1 x=hhx(i)+xx4 ro2=x*x+y*y+z*z ro2=amax1(ro2,ar2) ro3=ro2*sqrt(ro2) i1=i+l1 i2=i1+n3 i3=i2+n3 i4=i3+n3 i5=i4+n3 i6=i5+n3 i7=i6+n3 i8=i7+n3 j1=i+l2 j2=j1+n2 j3=j2+n2 j4=j3+n2 j5=j4+n2 j6=j5+n2 j7=j6+n2 j8=j7+n2 j9=j8+n2 j10=j9+n2 j11=j10+n2 c jj1=i+l3 c u(i1+m7)=u(i1+m7) 1 +dx1*(f(i4+n22)*f(i6+n22)-f(i2+n22)*f(i8+n22) 2 +f(i4+n12)*f(i6+n12)-f(i2+n12)*f(i8+n12) 3 +f(i4+n11)*f(i6+n11)-f(i2+n11)*f(i8+n11) 4 +f(i4+1)*f(i6+1)-f(i2+1)*f(i8+1) 5 -f(i4+n21)*f(i6+n21)+f(i2+n21)*f(i8+n21) 6 -f(i4+n2)*f(i6+n2)+f(i2+n2)*f(i8+n2) 7 -f(i4+n1)*f(i6+n1)+f(i2+n1)*f(i8+n1) 8 -f(i4)*f(i6)+f(i2)*f(i8)) 1 -dy1*(f(i3+n22)*f(i8+n22)-f(i4+n22)*f(i7+n22) 2 -f(i3+n12)*f(i8+n12)+f(i4+n12)*f(i7+n12) 3 +f(i3+n11)*f(i8+n11)-f(i4+n11)*f(i7+n11) 4 -f(i3+1)*f(i8+1)+f(i4+1)*f(i7+1) 5 +f(i3+n21)*f(i8+n21)-f(i4+n21)*f(i7+n21) 6 -f(i3+n2)*f(i8+n2)+f(i4+n2)*f(i7+n2) 7 +f(i3+n1)*f(i8+n1)-f(i4+n1)*f(i7+n1) 8 -f(i3)*f(i8)+f(i4)*f(i7)) u(i1+m6)=u(i1+m6) 1 +dz1*(f(i3+n22)*f(i8+n22)-f(i4+n22)*f(i7+n22) 2 +f(i3+n12)*f(i8+n12)-f(i4+n12)*f(i7+n12) 3 -f(i3+n11)*f(i8+n11)+f(i4+n11)*f(i7+n11) 4 -f(i3+1)*f(i8+1)+f(i4+1)*f(i7+1) 5 +f(i3+n21)*f(i8+n21)-f(i4+n21)*f(i7+n21) 6 +f(i3+n2)*f(i8+n2)-f(i4+n2)*f(i7+n2) 7 -f(i3+n1)*f(i8+n1)+f(i4+n1)*f(i7+n1) 8 -f(i3)*f(i8)+f(i4)*f(i7)) 1 -dx1*(f(i2+n22)*f(i7+n22)-f(i3+n22)*f(i6+n22) 2 +f(i2+n12)*f(i7+n12)-f(i3+n12)*f(i6+n12) 3 +f(i2+n11)*f(i7+n11)-f(i3+n11)*f(i6+n11) 4 +f(i2+1)*f(i7+1)-f(i3+1)*f(i6+1) 5 -f(i2+n21)*f(i7+n21)+f(i3+n21)*f(i6+n21) 6 -f(i2+n2)*f(i7+n2)+f(i3+n2)*f(i6+n2) 7 -f(i2+n1)*f(i7+n1)+f(i3+n1)*f(i6+n1) 8 -f(i2)*f(i7)+f(i3)*f(i6)) u(i1+m5)=u(i1+m5) 1 +dy1*(f(i2+n22)*f(i7+n22)-f(i3+n22)*f(i6+n22) 2 -f(i2+n12)*f(i7+n12)+f(i3+n12)*f(i6+n12) 3 +f(i2+n11)*f(i7+n11)-f(i3+n11)*f(i6+n11) 4 -f(i2+1)*f(i7+1)+f(i3+1)*f(i6+1) 5 +f(i2+n21)*f(i7+n21)-f(i3+n21)*f(i6+n21) 6 -f(i2+n2)*f(i7+n2)+f(i3+n2)*f(i6+n2) 7 +f(i2+n1)*f(i7+n1)-f(i3+n1)*f(i6+n1) 8 -f(i2)*f(i7)+f(i3)*f(i6)) 1 -dz1*(f(i4+n22)*f(i6+n22)-f(i2+n22)*f(i8+n22) 2 +f(i4+n12)*f(i6+n12)-f(i2+n12)*f(i8+n12) 3 -f(i4+n11)*f(i6+n11)+f(i2+n11)*f(i8+n11) 4 -f(i4+1)*f(i6+1)+f(i2+1)*f(i8+1) 5 +f(i4+n21)*f(i6+n21)-f(i2+n21)*f(i8+n21) 6 +f(i4+n2)*f(i6+n2)-f(i2+n2)*f(i8+n2) 7 -f(i4+n1)*f(i6+n1)+f(i2+n1)*f(i8+n1) 8 -f(i4)*f(i6)+f(i2)*f(i8)) c x3=dx1*(f(i2+n22)+f(i2+n12)+f(i2+n11)+f(i2+1) 1 -f(i2+n21)-f(i2+n2)-f(i2+n1)-f(i2)) 2 +dy1*(f(i3+n22)-f(i3+n12)+f(i3+n11)-f(i3+1) 3 +f(i3+n21)-f(i3+n2)+f(i3+n1)-f(i3)) 4 +dz1*(f(i4+n22)+f(i4+n12)-f(i4+n11)-f(i4+1) 5 +f(i4+n21)+f(i4+n2)-f(i4+n1)-f(i4)) u(i1+m4)=u(i1+m4)-p(j5)*gam*x3 1 -dx1*p(j2)*(f(i5+n22)+f(i5+n12)+f(i5+n11)+f(i5+1) 2 -f(i5+n21)-f(i5+n2)-f(i5+n1)-f(i5)) 3 -dy1*p(j3)*(f(i5+n22)-f(i5+n12)+f(i5+n11)-f(i5+1) 4 +f(i5+n21)-f(i5+n2)+f(i5+n1)-f(i5)) 5 -dz1*p(j4)*(f(i5+n22)+f(i5+n12)-f(i5+n11)-f(i5+1) 6 +f(i5+n21)+f(i5+n2)-f(i5+n1)-f(i5)) u(i1+m3)=u(i1+m3)+t1*((p(j9)*p(j7)-p(j10)*p(j6))/p(j1) 1 -eud*p(j4)-gra*z/ro3) 2 -dx1*p(j2)*(f(i4+n22)+f(i4+n12)+f(i4+n11)+f(i4+1) 3 -f(i4+n21)-f(i4+n2)-f(i4+n1)-f(i4)) 4 -dy1*p(j3)*(f(i4+n22)-f(i4+n12)+f(i4+n11)-f(i4+1) 5 +f(i4+n21)-f(i4+n2)+f(i4+n1)-f(i4)) 6 -dz1*p(j4)*(f(i4+n22)+f(i4+n12)-f(i4+n11)-f(i4+1) 7 +f(i4+n21)+f(i4+n2)-f(i4+n1)-f(i4)) 8 -dz1*(f(i5+n22)+f(i5+n12)-f(i5+n11)-f(i5+1) 9 +f(i5+n21)+f(i5+n2)-f(i5+n1)-f(i5))/p(j1) u(i1+m2)=u(i1+m2)+t1*((p(j11)*p(j6)-p(j9)*p(j8))/p(j1) 1 -eud*p(j3)-gra*y/ro3) 2 -dx1*p(j2)*(f(i3+n22)+f(i3+n12)+f(i3+n11)+f(i3+1) 3 -f(i3+n21)-f(i3+n2)-f(i3+n1)-f(i3)) 4 -dy1*p(j3)*(f(i3+n22)-f(i3+n12)+f(i3+n11)-f(i3+1) 5 +f(i3+n21)-f(i3+n2)+f(i3+n1)-f(i3)) 6 -dz1*p(j4)*(f(i3+n22)+f(i3+n12)-f(i3+n11)-f(i3+1) 7 +f(i3+n21)+f(i3+n2)-f(i3+n1)-f(i3)) 8 -dy1*(f(i5+n22)-f(i5+n12)+f(i5+n11)-f(i5+1) 9 +f(i5+n21)-f(i5+n2)+f(i5+n1)-f(i5))/p(j1) u(i1+n3)=u(i1+n3)+t1*((p(j10)*p(j8)-p(j11)*p(j7))/p(j1) 1 -eud*p(j2)-gra*x/ro3) 2 -dx1*p(j2)*(f(i2+n22)+f(i2+n12)+f(i2+n11)+f(i2+1) 3 -f(i2+n21)-f(i2+n2)-f(i2+n1)-f(i2)) 4 -dy1*p(j3)*(f(i2+n22)-f(i2+n12)+f(i2+n11)-f(i2+1) 5 +f(i2+n21)-f(i2+n2)+f(i2+n1)-f(i2)) 6 -dz1*p(j4)*(f(i2+n22)+f(i2+n12)-f(i2+n11)-f(i2+1) 7 +f(i2+n21)+f(i2+n2)-f(i2+n1)-f(i2)) 8 -dx1*(f(i5+n22)+f(i5+n12)+f(i5+n11)+f(i5+1) 9 -f(i5+n21)-f(i5+n2)-f(i5+n1)-f(i5))/p(j1) u(i1)=u(i1) 1 -dx1*(f(i2+n22)*f(i1+n22)+f(i2+n12)*f(i1+n12) 2 +f(i2+n11)*f(i1+n11)+f(i2+1)*f(i1+1) 3 -f(i2+n21)*f(i1+n21)-f(i2+n2)*f(i1+n2) 4 -f(i2+n1)*f(i1+n1)-f(i2)*f(i1)) 5 -dy1*(f(i3+n22)*f(i1+n22)-f(i3+n12)*f(i1+n12) 6 +f(i3+n11)*f(i1+n11)-f(i3+1)*f(i1+1) 7 +f(i3+n21)*f(i1+n21)-f(i3+n2)*f(i1+n2) 8 +f(i3+n1)*f(i1+n1)-f(i3)*f(i1)) 1 -dz1*(f(i4+n22)*f(i1+n22)+f(i4+n12)*f(i1+n12) 2 -f(i4+n11)*f(i1+n11)-f(i4+1)*f(i1+1) 3 +f(i4+n21)*f(i1+n21)+f(i4+n2)*f(i1+n2) 4 -f(i4+n1)*f(i1+n1)-f(i4)*f(i1)) 501 continue c do 520 i=1,nx1 x=hhx(i)+xx4 ro2=x*x+y*y+z*z i1=i i2=i1+n1 i3=i2+n1 i4=i3+n1 i5=i4+n1 i6=i5+n1 i7=i6+n1 i8=i7+n1 j1=i+l1 x2=amax1(ro2/ar2-1.0,0.0) x2=aru*x2*x2 x2=x2/(x2+1.0) u(j1+m7)=u(j1+m7)*x2+v(i8)*(1.0-x2) u(j1+m6)=u(j1+m6)*x2+v(i7)*(1.0-x2) u(j1+m5)=u(j1+m5)*x2+v(i6)*(1.0-x2) u(j1+m4)=u(j1+m4)*x2+v(i5)*(1.0-x2) u(j1+m3)=u(j1+m3)*x2+v(i4)*(1.0-x2) u(j1+m2)=u(j1+m2)*x2+v(i3)*(1.0-x2) u(j1+n3)=u(j1+n3)*x2+v(i2)*(1.0-x2) u(j1)=u(j1)*x2+v(i1)*(1.0-x2) 520 continue c do 522 i=1,nx1 i1=i+l1 i2=i1+n3 i3=i2+n3 i4=i3+n3 i5=i4+n3 p(i+k4)=u(i5) p(i+k3)=u(i4) p(i+k2)=u(i3) p(i+n2)=u(i2) p(i)=u(i1) 522 continue c do 524 i=1,nx1 i1=i+l1 x5=p(i+k4) x1=p(i) y5=amax1(x5,pr02) y1=amax1(x1,ro02) c x2=p(i+n2)*p(i+n2)+p(i+k2)*p(i+k2)+p(i+k3)*p(i+k3) c vmax=amax1(vmax,x2) x=amax1(x1,0.0)/y1 u(i1+m4)=y5 u(i1+m3)=x*p(i+k3) u(i1+m2)=x*p(i+k2) u(i1+n3)=x*p(i+n2) u(i1)=y1 524 continue 50 continue c c second step c do 62 m=1,nb l1=n2*(k-1)+n3*(m-1) l2=n2+k3*(m-1) do 62 i=1,n2 i1=i+l1 j1=i+l2 ff(j1+n2)=f(i1+n2) ff(j1)=f(i1) 62 continue c c current do 68 j=2,ny1 l1=n1*(j-1)+n2*(k-1) l2=n1*(j-1) l3=l2+n2 do 68 i=2,nx1 i6=i+l1+m5 i7=i6+n3 i8=i7+n3 j9=i+l2+k8 p(j9+k2)=0.25*((u(i7)+u(i7-n1)+u(i7-n2)+u(i7-n21) 1 -u(i7-1)-u(i7-n11)-u(i7-n12)-u(i7-n22))/hx 2 -(u(i6)-u(i6-n1)+u(i6-n2)-u(i6-n21) 3 +u(i6-1)-u(i6-n11)+u(i6-n12)-u(i6-n22))/hy) p(j9+n2)=0.25*((u(i6)+u(i6-n1)-u(i6-n2)-u(i6-n21) 1 +u(i6-1)+u(i6-n11)-u(i6-n12)-u(i6-n22))/hz 2 -(u(i8)+u(i8-n1)+u(i8-n2)+u(i8-n21) 3 -u(i8-1)-u(i8-n11)-u(i8-n12)-u(i8-n22))/hx) p(j9)=0.25*((u(i8)-u(i8-n1)+u(i8-n2)-u(i8-n21) 1 +u(i8-1)-u(i8-n11)+u(i8-n12)-u(i8-n22))/hy 2 -(u(i7)+u(i7-n1)-u(i7-n2)-u(i7-n21) 3 +u(i7-1)+u(i7-n11)-u(i7-n12)-u(i7-n22))/hz) 68 continue do 682 j=2,ny1 l1=np1*(j-1)+np2*(k-1) l2=n1*(j-1)+k8 do 684 i=2,np1 j9=i+np1b-1+l2 i1=i+l1 i2=i1+np3 i3=i2+np3 p(j9+k2)=p(j9+k2)-0.125*(pp(i3)+pp(i3-1)+pp(i3-np1)+pp(i3-np2) 1 +pp(i3-np11)+pp(i3-np12)+pp(i3-np21)+pp(i3-np22)) p(j9+n2)=p(j9+n2)-0.125*(pp(i2)+pp(i2-1)+pp(i2-np1)+pp(i2-np2) 1 +pp(i2-np11)+pp(i2-np12)+pp(i2-np21)+pp(i2-np22)) p(j9)=p(j9)-0.125*(pp(i1)+pp(i1-1)+pp(i1-np1)+pp(i1-np2) 1 +pp(i1-np11)+pp(i1-np12)+pp(i1-np21)+pp(i1-np22)) 684 continue do 682 i=2,np1b j9=np1b-i+2+l2 i1=i+l1 i2=i1+np3 i3=i2+np3 p(j9+k2)=p(j9+k2)+0.125*(pp(i3)+pp(i3-1)+pp(i3-np1)+pp(i3-np2) 1 +pp(i3-np11)+pp(i3-np12)+pp(i3-np21)+pp(i3-np22)) p(j9+n2)=p(j9+n2)+0.125*(pp(i2)+pp(i2-1)+pp(i2-np1)+pp(i2-np2) 1 +pp(i2-np11)+pp(i2-np12)+pp(i2-np21)+pp(i2-np22)) p(j9)=p(j9)-0.125*(pp(i1)+pp(i1-1)+pp(i1-np1)+pp(i1-np2) 1 +pp(i1-np11)+pp(i1-np12)+pp(i1-np21)+pp(i1-np22)) 682 continue c do 681 m=1,nb do 681 j=2,ny1 l2=n1*(j-1)+n2*(m-1) l1=n1*(j-1)+n2*(k-1)+n3*(m-1) do 681 i=2,nx1 i1=i+l1 p(i+l2)=0.125*(u(i1)+u(i1-1)+u(i1-n1)+u(i1-n11) 1 +u(i1-n2)+u(i1-n12)+u(i1-n21)+u(i1-n22)) 681 continue c c second step xx4=0.5*hx*float(2*nxp-nx1-2) z=0.5*hz*float(2*k-3) do 80 j=2,ny1 l1=n1*(j-1)+n2*(k-1) l2=n1*(j-1) l3=l2+n2 y=0.5*hy*float(2*j-3) c do 801 i=2,nx1 x=hhx(i)+xx4 ro2=x*x+y*y+z*z ro2=amax1(ro2,ar2) ro3=ro2*sqrt(ro2) i1=i+l1 i2=i1+n3 i3=i2+n3 i4=i3+n3 i5=i4+n3 i6=i5+n3 i7=i6+n3 i8=i7+n3 if1=i+l3 if2=if1+k3 if3=if2+k3 if4=if3+k3 if5=if4+k3 if6=if5+k3 if7=if6+k3 if8=if7+k3 j1=i+l2 j2=j1+n2 j3=j2+n2 j4=j3+n2 j5=j4+n2 j6=j5+n2 j7=j6+n2 j8=j7+n2 j9=j8+n2 j10=j9+n2 j11=j10+n2 jj1=i+l1 c x1=ff(if5)/(ff(if1)*po0) x2=ro2/ar2-1.0 x2=aru*x2*x2 x2=x2/(x2+1.0) eat=eat0*x2 c eat=eat0*x2/sqrt(x1*x1*x1) c x1=eat*(dx4*(ff(if8+1)-2.0*ff(if8)+ff(if8-1)) 1 +dy4*(ff(if8+n1)-2.0*ff(if8)+ff(if8-n1)) 2 +dz4*(ff(if8+n2)-2.0*ff(if8)+ff(if8-n2))) f(jj1+m7)=ff(if8)+x1 1 +dx2*(u(i4)*u(i6)-u(i2)*u(i8) 2 +u(i4-n1)*u(i6-n1)-u(i2-n1)*u(i8-n1) 3 +u(i4-n2)*u(i6-n2)-u(i2-n2)*u(i8-n2) 4 +u(i4-n21)*u(i6-n21)-u(i2-n21)*u(i8-n21) 5 -u(i4-1)*u(i6-1)+u(i2-1)*u(i8-1) 6 -u(i4-n11)*u(i6-n11)+u(i2-n11)*u(i8-n11) 7 -u(i4-n12)*u(i6-n12)+u(i2-n12)*u(i8-n12) 8 -u(i4-n22)*u(i6-n22)+u(i2-n22)*u(i8-n22)) 1 -dy2*(u(i3)*u(i8)-u(i4)*u(i7) 2 -u(i3-n1)*u(i8-n1)+u(i4-n1)*u(i7-n1) 3 +u(i3-n2)*u(i8-n2)-u(i4-n2)*u(i7-n2) 4 -u(i3-n21)*u(i8-n21)+u(i4-n21)*u(i7-n21) 5 +u(i3-1)*u(i8-1)-u(i4-1)*u(i7-1) 6 -u(i3-n11)*u(i8-n11)+u(i4-n11)*u(i7-n11) 7 +u(i3-n12)*u(i8-n12)-u(i4-n12)*u(i7-n12) 8 -u(i3-n22)*u(i8-n22)+u(i4-n22)*u(i7-n22)) x1=eat*(dx4*(ff(if7+1)-2.0*ff(if7)+ff(if7-1)) 1 +dy4*(ff(if7+n1)-2.0*ff(if7)+ff(if7-n1)) 2 +dz4*(ff(if7+n2)-2.0*ff(if7)+ff(if7-n2))) f(jj1+m6)=ff(if7)+x1 1 +dz2*(u(i3)*u(i8)-u(i4)*u(i7) 2 +u(i3-n1)*u(i8-n1)-u(i4-n1)*u(i7-n1) 3 -u(i3-n2)*u(i8-n2)+u(i4-n2)*u(i7-n2) 4 -u(i3-n21)*u(i8-n21)+u(i4-n21)*u(i7-n21) 5 +u(i3-1)*u(i8-1)-u(i4-1)*u(i7-1) 6 +u(i3-n11)*u(i8-n11)-u(i4-n11)*u(i7-n11) 7 -u(i3-n12)*u(i8-n12)+u(i4-n12)*u(i7-n12) 8 -u(i3-n22)*u(i8-n22)+u(i4-n22)*u(i7-n22)) 1 -dx2*(u(i2)*u(i7)-u(i3)*u(i6) 2 +u(i2-n1)*u(i7-n1)-u(i3-n1)*u(i6-n1) 3 +u(i2-n2)*u(i7-n2)-u(i3-n2)*u(i6-n2) 4 +u(i2-n21)*u(i7-n21)-u(i3-n21)*u(i6-n21) 5 -u(i2-1)*u(i7-1)+u(i3-1)*u(i6-1) 6 -u(i2-n11)*u(i7-n11)+u(i3-n11)*u(i6-n11) 7 -u(i2-n12)*u(i7-n12)+u(i3-n12)*u(i6-n12) 8 -u(i2-n22)*u(i7-n22)+u(i3-n22)*u(i6-n22)) x1=eat*(dx4*(ff(if6+1)-2.0*ff(if6)+ff(if6-1)) 1 +dy4*(ff(if6+n1)-2.0*ff(if6)+ff(if6-n1)) 2 +dz4*(ff(if6+n2)-2.0*ff(if6)+ff(if6-n2))) f(jj1+m5)=ff(if6)+x1 1 +dy2*(u(i2)*u(i7)-u(i3)*u(i6) 2 -u(i2-n1)*u(i7-n1)+u(i3-n1)*u(i6-n1) 3 +u(i2-n2)*u(i7-n2)-u(i3-n2)*u(i6-n2) 4 -u(i2-n21)*u(i7-n21)+u(i3-n21)*u(i6-n21) 5 +u(i2-1)*u(i7-1)-u(i3-1)*u(i6-1) 6 -u(i2-n11)*u(i7-n11)+u(i3-n11)*u(i6-n11) 7 +u(i2-n12)*u(i7-n12)-u(i3-n12)*u(i6-n12) 8 -u(i2-n22)*u(i7-n22)+u(i3-n22)*u(i6-n22)) 1 -dz2*(u(i4)*u(i6)-u(i2)*u(i8) 2 +u(i4-n1)*u(i6-n1)-u(i2-n1)*u(i8-n1) 3 -u(i4-n2)*u(i6-n2)+u(i2-n2)*u(i8-n2) 4 -u(i4-n21)*u(i6-n21)+u(i2-n21)*u(i8-n21) 5 +u(i4-1)*u(i6-1)-u(i2-1)*u(i8-1) 6 +u(i4-n11)*u(i6-n11)-u(i2-n11)*u(i8-n11) 7 -u(i4-n12)*u(i6-n12)+u(i2-n12)*u(i8-n12) 8 -u(i4-n22)*u(i6-n22)+u(i2-n22)*u(i8-n22)) c x3=dx2*(u(i2)+u(i2-n1)+u(i2-n2)+u(i2-n21) 1 -u(i2-1)-u(i2-n11)-u(i2-n12)-u(i2-n22)) 2 +dy2*(u(i3)-u(i3-n1)+u(i3-n2)-u(i3-n21) 3 +u(i3-1)-u(i3-n11)+u(i3-n12)-u(i3-n22)) 4 +dz2*(u(i4)+u(i4-n1)-u(i4-n2)-u(i4-n21) 5 +u(i4-1)+u(i4-n11)-u(i4-n12)-u(i4-n22)) f(jj1+m4)=ff(if5)-p(j5)*gam*x3 1 -dx2*p(j2)*(u(i5)+u(i5-n1)+u(i5-n2)+u(i5-n21) 2 -u(i5-1)-u(i5-n11)-u(i5-n12)-u(i5-n22)) 3 -dy2*p(j3)*(u(i5)-u(i5-n1)+u(i5-n2)-u(i5-n21) 4 +u(i5-1)-u(i5-n11)+u(i5-n12)-u(i5-n22)) 5 -dz2*p(j4)*(u(i5)+u(i5-n1)-u(i5-n2)-u(i5-n21) 6 +u(i5-1)+u(i5-n11)-u(i5-n12)-u(i5-n22)) 7 +pmu*(dx4*(ff(if5+1)-2.0*ff(if5)+ff(if5-1)) 8 +dy4*(ff(if5+n1)-2.0*ff(if5)+ff(if5-n1)) 9 +dz4*(ff(if5+n2)-2.0*ff(if5)+ff(if5-n2))) f(jj1+m3)=ff(if4)+t*((p(j9)*p(j7)-p(j10)*p(j6))/p(j1) 1 -eud*p(j4)-gra*z/ro3) 2 -dx2*p(j2)*(u(i4)+u(i4-n1)+u(i4-n2)+u(i4-n21) 3 -u(i4-1)-u(i4-n11)-u(i4-n12)-u(i4-n22)) 4 -dy2*p(j3)*(u(i4)-u(i4-n1)+u(i4-n2)-u(i4-n21) 5 +u(i4-1)-u(i4-n11)+u(i4-n12)-u(i4-n22)) 6 -dz2*p(j4)*(u(i4)+u(i4-n1)-u(i4-n2)-u(i4-n21) 7 +u(i4-1)+u(i4-n11)-u(i4-n12)-u(i4-n22)) 8 -dz2*(u(i5)+u(i5-n1)-u(i5-n2)-u(i5-n21) 9 +u(i5-1)+u(i5-n11)-u(i5-n12)-u(i5-n22))/p(j1) 1 +rmu*(dx4*(ff(if4+1)-2.0*ff(if4)+ff(if4-1)) 2 +dy4*(ff(if4+n1)-2.0*ff(if4)+ff(if4-n1)) 3 +dz4*(ff(if4+n2)-2.0*ff(if4)+ff(if4-n2)))/ff(if1) f(jj1+m2)=ff(if3)+t*((p(j11)*p(j6)-p(j9)*p(j8))/p(j1) 1 -eud*p(j3)-gra*y/ro3) 2 -dx2*p(j2)*(u(i3)+u(i3-n1)+u(i3-n2)+u(i3-n21) 3 -u(i3-1)-u(i3-n11)-u(i3-n12)-u(i3-n22)) 4 -dy2*p(j3)*(u(i3)-u(i3-n1)+u(i3-n2)-u(i3-n21) 5 +u(i3-1)-u(i3-n11)+u(i3-n12)-u(i3-n22)) 6 -dz2*p(j4)*(u(i3)+u(i3-n1)-u(i3-n2)-u(i3-n21) 7 +u(i3-1)+u(i3-n11)-u(i3-n12)-u(i3-n22)) 8 -dy2*(u(i5)-u(i5-n1)+u(i5-n2)-u(i5-n21) 9 +u(i5-1)-u(i5-n11)+u(i5-n12)-u(i5-n22))/p(j1) 1 +rmu*(dx4*(ff(if3+1)-2.0*ff(if3)+ff(if3-1)) 2 +dy4*(ff(if3+n1)-2.0*ff(if3)+ff(if3-n1)) 3 +dz4*(ff(if3+n2)-2.0*ff(if3)+ff(if3-n2)))/ff(if1) f(jj1+n3)=ff(if2)+t*((p(j10)*p(j8)-p(j11)*p(j7))/p(j1) 1 -eud*p(j2)-gra*x/ro3) 2 -dx2*p(j2)*(u(i2)+u(i2-n1)+u(i2-n2)+u(i2-n21) 3 -u(i2-1)-u(i2-n11)-u(i2-n12)-u(i2-n22)) 4 -dy2*p(j3)*(u(i2)-u(i2-n1)+u(i2-n2)-u(i2-n21) 5 +u(i2-1)-u(i2-n11)+u(i2-n12)-u(i2-n22)) 6 -dz2*p(j4)*(u(i2)+u(i2-n1)-u(i2-n2)-u(i2-n21) 7 +u(i2-1)+u(i2-n11)-u(i2-n12)-u(i2-n22)) 8 -dx2*(u(i5)+u(i5-n1)+u(i5-n2)+u(i5-n21) 9 -u(i5-1)-u(i5-n11)-u(i5-n12)-u(i5-n22))/p(j1) 1 +rmu*(dx4*(ff(if2+1)-2.0*ff(if2)+ff(if2-1)) 2 +dy4*(ff(if2+n1)-2.0*ff(if2)+ff(if2-n1)) 3 +dz4*(ff(if2+n2)-2.0*ff(if2)+ff(if2-n2)))/ff(if1) f(jj1)=ff(if1) 1 -dx2*(u(i2)*u(i1)+u(i2-n1)*u(i1-n1) 2 +u(i2-n2)*u(i1-n2)+u(i2-n21)*u(i1-n21) 3 -u(i2-1)*u(i1-1)-u(i2-n11)*u(i1-n11) 4 -u(i2-n12)*u(i1-n12)-u(i2-n22)*u(i1-n22)) 5 -dy2*(u(i3)*u(i1)-u(i3-n1)*u(i1-n1) 6 +u(i3-n2)*u(i1-n2)-u(i3-n21)*u(i1-n21) 7 +u(i3-1)*u(i1-1)-u(i3-n11)*u(i1-n11) 8 +u(i3-n12)*u(i1-n12)-u(i3-n22)*u(i1-n22)) 1 -dz2*(u(i4)*u(i1)+u(i4-n1)*u(i1-n1) 2 -u(i4-n2)*u(i1-n2)-u(i4-n21)*u(i1-n21) 3 +u(i4-1)*u(i1-1)+u(i4-n11)*u(i1-n11) 4 -u(i4-n12)*u(i1-n12)-u(i4-n22)*u(i1-n22)) 5 +dfu*(dx4*(ff(if1+1)-2.0*ff(if1)+ff(if1-1)) 6 +dy4*(ff(if1+n1)-2.0*ff(if1)+ff(if1-n1)) 7 +dz4*(ff(if1+n2)-2.0*ff(if1)+ff(if1-n2))) 801 continue c do 820 i=2,nx1 x=hhx(i)+xx4 ro2=x*x+y*y+z*z i1=i+l3 i2=i1+k3 i3=i2+k3 i4=i3+k3 i5=i4+k3 i6=i5+k3 i7=i6+k3 i8=i7+k3 j1=i+l1 x2=amax1(ro2/ar2-1.0,0.0) x2=aru*x2*x2 x2=x2/(x2+1.0) f(j1+m7)=f(j1+m7)*x2+ff(i8)*(1.0-x2) f(j1+m6)=f(j1+m6)*x2+ff(i7)*(1.0-x2) f(j1+m5)=f(j1+m5)*x2+ff(i6)*(1.0-x2) f(j1+m4)=f(j1+m4)*x2+ff(i5)*(1.0-x2) f(j1+m3)=f(j1+m3)*x2+ff(i4)*(1.0-x2) f(j1+m2)=f(j1+m2)*x2+ff(i3)*(1.0-x2) f(j1+n3)=f(j1+n3)*x2+ff(i2)*(1.0-x2) f(j1)=f(j1)*x2+ff(i1)*(1.0-x2) 820 continue c do 822 i=2,nx1 i1=i+l1 i2=i1+n3 i3=i2+n3 i4=i3+n3 i5=i4+n3 p(i+k4)=f(i5) p(i+k3)=f(i4) p(i+k2)=f(i3) p(i+n2)=f(i2) p(i)=f(i1) 822 continue c do 824 i=2,nx1 i1=i+l1 x5=p(i+k4) x1=p(i) y5=amax1(x5,pr02) y1=amax1(x1,ro02) x2=p(i+n2)*p(i+n2)+p(i+k2)*p(i+k2)+p(i+k3)*p(i+k3) vmax=amax1(vmax,x2) x=amax1(x1,0.0)/y1 f(i1+m4)=y5 f(i1+m3)=x*p(i+k3) f(i1+m2)=x*p(i+k2) f(i1+n3)=x*p(i+n2) f(i1)=y1 824 continue 80 continue c do 830 m=1,nb l1=k3*(m-1) do 830 i=1,n2 i1=i+l1 c u(i1)=u(i1+n2) ff(i1)=ff(i1+n2) 830 continue 90 continue c c time=time+t c if(ii.ne.1.and.iiq.ne.iiq0) go to 399 t1=t dx1=0.25*t1/hx dy1=0.25*t1/hy dz1=0.25*t1/hz dx3=t1/(hx*hx) dy3=t1/(hy*hy) dz3=t1/(hz*hz) 399 continue 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 x2=2.0 x1=2.0*x2 nx11=nx1-2*nxp do 102 k=1,nz2 do 102 j=1,ny2 l1=n1*(j-1)+n2*(k-1) do 102 i=1,nx2 x=(2.0*hhx(i)/hx-float(2+nx11))/float(nx11) x=(x1+1.0)*x+x1 x=-x x=amax1(x,0.0) x=amin1(x,1.0) x=x*x i1=i+l1 f(i1+m7)=(1.0-x)*f(i1+m7)+bis*x f(i1+m6)=(1.0-x)*f(i1+m6) f(i1+m5)=(1.0-x)*f(i1+m5) f(i1+m4)=(1.0-x)*f(i1+m4)+pr01*x f(i1+m3)=(1.0-x)*f(i1+m3) f(i1+m2)=(1.0-x)*f(i1+m2) f(i1+n3)=(1.0-x)*f(i1+n3)+vsw*x f(i1)=(1.0-x)*f(i1)+ro01*x 102 continue c do 104 k=1,nz1 do 104 j=1,ny1 l1=n1*(j-1)+n2*(k-1) do 104 i=1,nx1 x=(2.0*hhx(i)/hx-float(1+nx11))/float(nx11) x=(x1+1.0)*x+x1 x=-x x=amax1(x,0.0) x=amin1(x,1.0) x=x*x i1=i+l1 u(i1+m7)=(1.0-x)*u(i1+m7)+bis*x u(i1+m6)=(1.0-x)*u(i1+m6) u(i1+m5)=(1.0-x)*u(i1+m5) u(i1+m4)=(1.0-x)*u(i1+m4)+pr01*x u(i1+m3)=(1.0-x)*u(i1+m3) u(i1+m2)=(1.0-x)*u(i1+m2) u(i1+n3)=(1.0-x)*u(i1+n3)+vsw*x u(i1)=(1.0-x)*u(i1)+ro01*x 104 continue c 403 if(iis.ne.iis0) go to 100 iis=0 c if(ii.lt.1024) go to 100 c c write(12) f itapo=11 402 do 173 j1=1,mm2 i1=1+mm1*(j1-1) i2=mm1*j1 write(itapo) (f(i),i=i1,i2) 173 continue write(6,661) ii,last,nin,itap,x2,vmax 661 format(1h ,5x,4i10,1p2e15.7) c if(vmax.gt.1.0) go to 9 c jyg=ny1+nz1 jyg2=jyg-2 jxg=n1*jyg jxn=nx2 lank=36 ip=0 lb=1 ll=8 x=2.0*ar2 do 250 m=1,ll i3=m if(m.gt.nb) i3=1 do 250 i=1,nx2 do 252 k=2,nz2 i1=i+n2*(k-1)+n3*(i3-1) i2=i+n1*(ny1+k-2)+jxg*(m-1) ff(i2)=0.5*(f(i1)+f(i1+n1)) 252 continue do 250 j=2,ny2 i1=i+n1*(j-1)+n3*(i3-1) i2=i+n1*(ny1+1-j)+jxg*(m-1) ff(i2)=0.5*(f(i1)+f(i1+n2)) 250 continue i1=ii c call scalm(nx,jyg2,nxp,lb,ll,lll,hx,hz,x,zmin,zmax,ff) ll1=2 ll2=ll/ll1 do 150 j=1,ll2 lb=1+ll1*(j-1) c call grapd(jxg,jyg,lb,ll1,zmin,zmax,i1,lll,jxn,lank,ip,ff) c call grape(jxg,jyg,lb,ll1,zmin,zmax,i1,lll,jxn,lank,ip,ff) 150 continue c jyg=nz2 jyg2=jyg-2 jxg=ny2*jyg jxn=ny2 lank=36 arr1=0.1 ip=0 do 260 il=1,nil lb=1 ll=8 i=nx2/2+1+(nx2*il)/(2*(nil+1)) do 262 m=1,ll i3=m if(m.gt.nb) i3=1 do 262 k=1,nz2 do 262 j=1,ny2 i1=i+n1*(j-1)+n2*(k-1)+n3*(i3-1) i2=j+ny2*(k-1)+jxg*(m-1) ff(i2)=f(i1) 262 continue i1=i c call scalm(ny,jyg2,nxp,lb,ll,lll,hx,hz,arr1,zmin,zmax,ff) ll1=4 ll2=ll/ll1 do 260 j=1,ll2 lb=1+ll1*(j-1) c call grapd(jxg,jyg,lb,ll1,zmin,zmax,i1,lll,jxn,lank,ip,ff) c call grape(jxg,jyg,lb,ll1,zmin,zmax,i1,lll,jxn,lank,ip,ff) 260 continue if(vmax.gt.1.0) go to 9 c 100 continue 300 continue 9 continue c rewind 12 stop end subroutine equib7(nx,ny,nz,nxp,ro01,pr01,rrat,cj,cp,f,p) dimension p(1),f(1),cj(1),cp(1) c cj 1-r1 2-dr 3-dz 4-b0 5-rm 6-a2 7-aid 8-ar1 9-gam 10-gra c cp 1-r0 2-z0 3-ra 4-vo0 5-p0 6-gra 7-07 8-vsw 9-09 10-10 c aid=cj(7) ar1=cj(8) gam=cj(9) gra=cj(10) p0=(gam-1.0)*gra/gam bis=cp(11)*1.0e-4 x1=1.0/(gam-1.0) x2=gam*x1 b0=cj(4) ro02=rrat*ro01 pr02=pr01 c ic0=4 icd=6 pi=3.1415926 nx1=nx+1 nx2=nx+2 ny1=ny+1 ny2=ny+2 nz1=nz+1 nz2=nz+2 n1=nx2 n2=n1*ny2 n3=n2*nz2 n11=n1+1 n12=n2+1 n21=n1+n2 n22=n21+1 nb=8 hx=cp(1)/float(nx1) hy=cp(2)/float(ny1) hz=cp(3)/float(nz1) vsw=cp(9) c do 10 k=1,nz2 do 10 j=1,ny2 vr1=vsw do 10 i=1,nx2 x=0.5*hx*float(2*i-nx2-1+2*nxp) y=0.5*hy*float(2*j-3) z=0.5*hz*float(2*k-3) ro2=x*x+y*y+z*z ro1=sqrt(ro2) ro3=ro1*ro2 i1=i+n1*(j-1)+n2*(k-1) i2=i1+n3 i3=i2+n3 i4=i3+n3 i5=i4+n3 i6=i5+n3 i7=i6+n3 i8=i7+n3 f(i1)=1.0/ro3 c f(i1)=ro1**(-x1) if(f(i1).lt.ro02) f(i1)=ro02 f(i5)=1.0e-7*cp(6)/ro2 c f(i5)=p0*ro1**(-x2) if(f(i5).lt.pr02) f(i5)=pr02 f(i6)=-2.0*b0*x*z/(ro2*ro2) f(i7)=-2.0*b0*y*z/(ro2*ro2) f(i8)=b0*(x*x+y*y-z*z)/(ro2*ro2) f(i6)=-3.0*b0*x*z/(ro2*ro3) f(i7)=-3.0*b0*y*z/(ro2*ro3) f(i8)=b0*(x*x+y*y-2.0*z*z)/(ro2*ro3) px1=sqrt(f(i6)*f(i6)+f(i7)*f(i7)+f(i8)*f(i8)) if(abs(vsw).lt.1.0e-8) go to 10 y1=(ro01-cj(1)*(px1/vsw)**2)/ro01 if(y1.le.0.0) y1=0.0 y1=vsw*sqrt(y1) vr1=amin1(vr1,y1) f(i2)=vr1 xx=vr1/vsw f(i5)=f(i5)+(pr01-pr02)*xx f(i8)=f(i8)+bis 10 continue c np1=nx2/2+nxp np1b=nx2/2-nxp np2=np1*ny1 np3=np2*nz1 c current do 38 k=1,nz1 do 38 j=1,ny1 l1=n1*(j-1)+n2*(k-1)+n3*5 l2=np1*(j-1)+np2*(k-1) do 38 i=1,np1 i6=i+np1b-1+l1 i7=i6+n3 i8=i7+n3 j9=i+l2 p(j9+2*np3)=0.25*((f(i7+n22)+f(i7+n12)+f(i7+n11)+f(i7+1) 1 -f(i7+n21)-f(i7+n2)-f(i7+n1)-f(i7))/hx 2 -(f(i6+n22)-f(i6+n12)+f(i6+n11)-f(i6+1) 3 +f(i6+n21)-f(i6+n2)+f(i6+n1)-f(i6))/hy) p(j9+np3)=0.25*((f(i6+n22)+f(i6+n12)-f(i6+n11)-f(i6+1) 1 +f(i6+n21)+f(i6+n2)-f(i6+n1)-f(i6))/hz 2 -(f(i8+n22)+f(i8+n12)+f(i8+n11)+f(i8+1) 3 -f(i8+n21)-f(i8+n2)-f(i8+n1)-f(i8))/hx) p(j9)=0.25*((f(i8+n22)-f(i8+n12)+f(i8+n11)-f(i8+1) 1 +f(i8+n21)-f(i8+n2)+f(i8+n1)-f(i8))/hy 2 -(f(i7+n22)+f(i7+n12)-f(i7+n11)-f(i7+1) 3 +f(i7+n21)+f(i7+n2)-f(i7+n1)-f(i7))/hz) 38 continue return end subroutine equib8(nx,ny,nz,nxp,ro01,pr01,rrat,cj,cp,f,p) dimension p(1),f(1),cj(1),cp(1) c cj 1-r1 2-dr 3-dz 4-b0 5-rm 6-a2 7-aid 8-ar1 9-gam 10-gra c cp 1-r0 2-z0 3-ra 4-vo0 5-p0 6-gra 7-07 8-vsw 9-09 10-10 c aid=cj(7) ar1=cj(8) gam=cj(9) gra=cj(10) p0=(gam-1.0)*gra/gam bis=cp(11)*1.0e-4 x1=1.0/(gam-1.0) x2=gam*x1 b0=cj(4) ro02=rrat*ro01 pr02=pr01 c ic0=4 icd=6 pi=3.1415926 nx1=nx+1 nx2=nx+2 ny1=ny+1 ny2=ny+2 nz1=nz+1 nz2=nz+2 n1=nx2 n2=n1*ny2 n3=n2*nz2 n11=n1+1 n12=n2+1 n21=n1+n2 n22=n21+1 nb=8 hx=cp(1)/float(nx1) hy=cp(2)/float(ny1) hz=cp(3)/float(nz1) vsw=cp(9) alp=cj(1) bet=cj(2) alp1=1.0/sqrt(alp) xm=(2.0*alp*b0*b0/(ro01*vsw*vsw))**(1.0/6.0) xm=-xm c do 10 k=1,nz2 do 10 j=1,ny2 c vr1=vsw do 10 i=1,nx2 x=0.5*hx*float(2*i-nx2-1+2*nxp) y=0.5*hy*float(2*j-3) z=0.5*hz*float(2*k-3) x1m=x-2.0*xm ro2=x*x+y*y+z*z ro1=sqrt(ro2) ro3=ro1*ro2 ro2m=x1m*x1m+y*y+z*z ro1m=sqrt(ro2m) ro3m=ro1m*ro2m i1=i+n1*(j-1)+n2*(k-1) i2=i1+n3 i3=i2+n3 i4=i3+n3 i5=i4+n3 i6=i5+n3 i7=i6+n3 i8=i7+n3 f(i1)=1.0/ro3 c f(i1)=ro1**(-x1) if(f(i1).lt.ro02) f(i1)=ro02 f(i5)=1.0e-7*cp(6)/ro2 c f(i5)=p0*ro1**(-x2) if(f(i5).lt.pr02) f(i5)=pr02 f(i6)=-2.0*b0*x*z/(ro2*ro2) f(i7)=-2.0*b0*y*z/(ro2*ro2) f(i8)=b0*(x*x+y*y-z*z)/(ro2*ro2) f(i6)=-3.0*b0*x*z/(ro2*ro3) f(i7)=-3.0*b0*y*z/(ro2*ro3) f(i8)=b0*(x*x+y*y-2.0*z*z)/(ro2*ro3) 10 continue c np1=nx2/2+nxp np1b=nx2/2-nxp np2=np1*ny1 np3=np2*nz1 c current do 38 k=1,nz1 do 38 j=1,ny1 l1=n1*(j-1)+n2*(k-1)+n3*5 l2=np1*(j-1)+np2*(k-1) do 38 i=1,np1 i6=i+np1b-1+l1 i7=i6+n3 i8=i7+n3 j9=i+l2 p(j9+2*np3)=0.25*((f(i7+n22)+f(i7+n12)+f(i7+n11)+f(i7+1) 1 -f(i7+n21)-f(i7+n2)-f(i7+n1)-f(i7))/hx 2 -(f(i6+n22)-f(i6+n12)+f(i6+n11)-f(i6+1) 3 +f(i6+n21)-f(i6+n2)+f(i6+n1)-f(i6))/hy) p(j9+np3)=0.25*((f(i6+n22)+f(i6+n12)-f(i6+n11)-f(i6+1) 1 +f(i6+n21)+f(i6+n2)-f(i6+n1)-f(i6))/hz 2 -(f(i8+n22)+f(i8+n12)+f(i8+n11)+f(i8+1) 3 -f(i8+n21)-f(i8+n2)-f(i8+n1)-f(i8))/hx) p(j9)=0.25*((f(i8+n22)-f(i8+n12)+f(i8+n11)-f(i8+1) 1 +f(i8+n21)-f(i8+n2)+f(i8+n1)-f(i8))/hy 2 -(f(i7+n22)+f(i7+n12)-f(i7+n11)-f(i7+1) 3 +f(i7+n21)+f(i7+n2)-f(i7+n1)-f(i7))/hz) 38 continue c do 40 k=1,nz2 do 40 j=1,ny2 c vr1=vsw do 40 i=1,nx2 x=0.5*hx*float(2*i-nx2-1+2*nxp) y=0.5*hy*float(2*j-3) z=0.5*hz*float(2*k-3) x1m=x-2.0*xm ro2=x*x+y*y+z*z ro1=sqrt(ro2) ro3=ro1*ro2 ro2m=x1m*x1m+y*y+z*z ro1m=sqrt(ro2m) ro3m=ro1m*ro2m i1=i+n1*(j-1)+n2*(k-1) i2=i1+n3 i3=i2+n3 i4=i3+n3 i5=i4+n3 i6=i5+n3 i7=i6+n3 i8=i7+n3 f(i1)=1.0/ro3 c f(i1)=ro1**(-x1) if(f(i1).lt.ro02) f(i1)=ro02 f(i5)=1.0e-7*cp(6)/ro2 c f(i5)=p0*ro1**(-x2) if(f(i5).lt.pr02) f(i5)=pr02 f(i6)=-2.0*b0*x*z/(ro2*ro2) f(i7)=-2.0*b0*y*z/(ro2*ro2) f(i8)=b0*(x*x+y*y-z*z)/(ro2*ro2) bx1=-3.0*b0*x*z/(ro2*ro3) by1=-3.0*b0*y*z/(ro2*ro3) bz1=b0*(x*x+y*y-2.0*z*z)/(ro2*ro3) bx2=-3.0*b0*x1m*z/(ro2m*ro3m) by2=-3.0*b0*y*z/(ro2m*ro3m) bz2=b0*(x1m*x1m+y*y-2.0*z*z)/(ro2m*ro3m) f(i6)=bx1+bx2 f(i7)=by1+by2 f(i8)=bz1+bz2 if(x.lt.xm) f(i6)=0.0 if(x.lt.xm) f(i7)=0.0 if(x.lt.xm) f(i8)=0.0 y1=alp1+(1.0-alp1)*(x-xm)/(bet-1.0)/xm if(x.gt.xm) y1=0.0 y1=amax1(y1,0.0) y1=amin1(y1,1.0) vr1=vsw*y1 c px1=sqrt(f(i6)*f(i6)+f(i7)*f(i7)+f(i8)*f(i8)) if(abs(vsw).lt.1.0e-8) go to 40 c y1=(ro01-cj(1)*(px1/vsw)**2)/ro01 c if(y1.le.0.0) y1=0.0 c y1=vsw*sqrt(y1) c vr1=amin1(vr1,y1) f(i2)=vr1 xx=vr1/vsw f(i5)=f(i5)+(pr01-pr02)*xx f(i8)=f(i8)+bis*y1 40 continue c return end