c program ogncb7 c c file name /hearthb3/earthb14.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) parameter (nx=320,ny= 80,nz= 80,nxp= 80) c parameter (nx=480,ny=160,nz=160,nxp= 80) c parameter (nx=640,ny=160,nz=160,nxp=200) c parameter (nx=800,ny=200,nz=200,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) parameter (last=16384,nil=2,lll=18,ibig=0) c parameter (last= 4096,nil=2,lll=18,ibig=0) c parameter (last=12288,nil=2,lll=18,ibig=0) parameter (iiq0=8,iip0= 32,iis0= 8192,thx=2.3450007) c parameter (iiq0=8,iip0= 32,iis0= 4096,thx=1.0416666) c parameter (iiq0=8,iip0= 32,iis0=12288,thx=1.0416666) c parameter (iiq0=8,iip0= 32,iis0= 448,thx=0.9523810) c c parameter (rmuu=0.0020,eud0=0.0000,ro01=5.0e-4,pr01=3.56e-8) parameter (rmuu=0.0020,eud0=0.0020,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,xta=100.0,xtd=10.0) c parameter (eatt=0.0010,rrat=0.1,rra1=0.1) c parameter (eatt=0.0020,rrat=0.2,rra1=0.2) c parameter (eatt=0.0020,rrat=0.1,rra1=0.1) parameter (eatt=0.0020,rrat=0.2,rra1=0.2) c parameter (eatt=0.0020,rrat=0.3,rra1=0.3) 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=2) 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 (orid= 90.000,ths0=0.0,kk8=k3*8) c parameter (mfd=6,nfd=1441,iinob=last*140, ttno=409.6) parameter (mfd=6,nfd=1441,iinob=last*30, ttno=409.6) c parameter (mfd=6,nfd=1441,iinob=last*1410, ttno=448.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) dimension hhx(n1),fbb(16), 1 cj(10),cp(11),v(n2) 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$ 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) c 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 08 vsw ibr ibz data cp/128.4,32.4,32.4,3.5,0.0,2.68,1.35,7.0,.044,0.0,-0.0/ c c c rewind 12 c c open(10,file='earthbh0.data', c 1 access='sequential',form='unformatted') c open(11,file='earthbh1.data', 1 access='sequential',form='unformatted') c c bisz= 1.5 c bisz= 0.0 c bisz=-1.5 vsww=cp(9) c do 300 iii=1,1 c cp(9)=vsww cp(11)=bisz ntap=10+iii 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 bis,(cp(i),i=1,11),(cj(j),j=1,10) 12 format(1h ,5x,11i10/(1h ,5x,1p8e15.5)) 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(7+nb)=-1.0 c c initial condition c call equib7(f,pp,nxp,ro01,pr01,rra1,cj,cp) if(iii.eq.1) 1 call equib8(f,pp,nxp,ro01,pr01,rra1,cj,cp) 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(ntap) f if(nin.eq.0) go to 400 if(iii.ge.2) go to 400 c do 410 ii=1,itap c do 174 m=1,nb c do 174 k=1,nz2 c read(10) ((gf(i,j,k,m),i=1,nx2),j=1,ny2) c 174 continue c 410 continue c 400 continue c c!hpf$ asynchronous(id1),nobuffer begin c f(1:nx2,1:ny2,1:nz2,1:nb)=gf(1:nx2,1:ny2,1:nz2,1:nb) c!hpf$ end asynchronous c!hpf$ asyncwait(id1) c do 410 ii=1,itap do 174 m=1,nb do 174 k=1,nz2 c read(10) gf c!hpf$ asynchronous(id1),nobuffer begin c f(1:nx2,1:ny2,k,m) = gf(1:nx2,1:ny2) c!hpf$ end asynchronous c!hpf$ asyncwait(id1) cc read(10) f(1:nx2,1:ny2,k,m) c read(10,end=400) f(1:nx2,1:ny2,k,m) 174 continue 410 continue 400 continue 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 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 c !hpf$ independent,new(i,j,k,m) do 31 k=1,nz2 !hpf$ on home(f(:,:,k,:)),local(f,i,j,m,x,xx,xx3,xx4, !hpf$* hhx) begin do 30 m=1,nb if(k.eq.1) then do 311 j=2,ny1 do 311 i=2,nx1 f(i,j,1,m)=f(i,j,2,m)*fbb(m) 311 continue else if(k.eq.nz2) then do 312 j=2,ny1 f(2,j,nz2,m)=f(1,j,nz1,m) c f(2,j,k,m)=f(1,j,k-1,m) do 312 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,j,nz2,m)=x*f(i-1,j,nz1,m)+xx*f(i-2,j,nz,m) c f(i,j,k,m)=x*f(i-1,j,k-1,m)+xx*f(i-2,j,k-2,m) 312 continue end if 30 continue !hpf$ end on 31 continue c c !hpf$ independent,new(i,j,k,m) do 32 k=1,nz2 !hpf$ on home(f(:,:,k,:)),local(f,u,fbb,i,j,m,x,xx, !hpf$* xx3,xx4,hhx,pr02,ro02) begin do 3311 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) c 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) c f(i,1,k,m)=f(i,1,k,m)*fbb(m+nb) 3211 continue c c boundary condition at nx=nx2 do 3311 j=1,ny2 f(nx2,j,k,m)=f(nx1,j,k,m) 3311 continue c c do 3411 j=1,ny2 do 3411 i=1,nx2 f(i,j,k,5)=amax1(f(i,j,k,5),pr02) f(i,j,k,1)=amax1(f(i,j,k,1),ro02) u(i,j,k,5)=amax1(u(i,j,k,5),pr02) u(i,j,k,1)=amax1(u(i,j,k,1),ro02) 3411 continue !hpf$ end on 32 continue c c !hpf$ reflect f c if(ii.ne.1.and.iiq.ne.iiq0) go to 172 c do 372 m=1,nb !hpf$ independent,new(i,j,k,m) do 371 k=1,nz1 !hpf$ on home(u(:,:,k,:)),local(f,u,i,j,m) begin do 3711 m=1,nb do 3711 j=1,ny1 do 3711 i=1,nx1 u(i,j,k,m)=0.125*(f(i,j,k,m)+f(i+1,j,k,m) 1 +f(i,j+1,k,m)+f(i+1,j+1,k,m) 2 +f(i,j,k+1,m)+f(i+1,j,k+1,m) 3 +f(i,j+1,k+1,m)+f(i+1,j+1,k+1,m)) 3711 continue !hpf$ end on 371 continue c c 372 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 c if(nin.ne.0.or.ii.ne.1) go to 172 c do 171 m=1,nb c do 171 k=1,nz2 c write(ntap) ((gf(i,j,k,m),i=1,nx2),j=1,ny2) c 171 continue 172 continue c c c c first step c step of k=k !hpf$ independent,new(i,j,k,m) do 90 k=1,nz1 !hpf$ on home(u(:,:,k,:)),local(f,u,p,pp,ff,i,j,m, !hpf$* x,y,z,hx,hy,hz,xx4,hhx,ar2, !hpf$* ro2,ro3,dx1,dy1,dz1,gam,eud,gra,ar1, !hpf$* uxd,x1,x2,x3,x5,y1,y5,ro02,pr02,vmax) begin c c current do 38 j=1,ny1 do 38 i=1,nx1 p(i,j,k,11)=0.25*((f(i+1,j+1,k+1,7)+f(i+1,j,k+1,7) 1 +f(i+1,j+1,k,7)+f(i+1,j,k,7) 2 -f(i,j+1,k+1,7)-f(i,j,k+1,7)-f(i,j+1,k,7)-f(i,j,k,7))/hx 3 -(f(i+1,j+1,k+1,6)-f(i+1,j,k+1,6)+f(i+1,j+1,k,6)-f(i+1,j,k,6) 4 +f(i,j+1,k+1,6)-f(i,j,k+1,6)+f(i,j+1,k,6)-f(i,j,k,6))/hy) p(i,j,k,10)=0.25*((f(i+1,j+1,k+1,6)+f(i+1,j,k+1,6) 1 -f(i+1,j+1,k,6)-f(i+1,j,k,6) 2 +f(i,j+1,k+1,6)+f(i,j,k+1,6)-f(i,j+1,k,6)-f(i,j,k,6))/hz 3 -(f(i+1,j+1,k+1,8)+f(i+1,j,k+1,8)+f(i+1,j+1,k,8)+f(i+1,j,k,8) 4 -f(i,j+1,k+1,8)-f(i,j,k+1,8)-f(i,j+1,k,8)-f(i,j,k,8))/hx) p(i,j,k,9)=0.25*((f(i+1,j+1,k+1,8)-f(i+1,j,k+1,8) 1 +f(i+1,j+1,k,8)-f(i+1,j,k,8) 2 +f(i,j+1,k+1,8)-f(i,j,k+1,8)+f(i,j+1,k,8)-f(i,j,k,8))/hy 3 -(f(i+1,j+1,k+1,7)+f(i+1,j,k+1,7)-f(i+1,j+1,k,7)-f(i+1,j,k,7) 4 +f(i,j+1,k+1,7)+f(i,j,k+1,7)-f(i,j+1,k,7)-f(i,j,k,7))/hz) 38 continue do 382 j=1,ny1 do 382 i=1,nx1 p(i,j,k,11)=p(i,j,k,11)-pp(i,j,k,3) p(i,j,k,10)=p(i,j,k,10)-pp(i,j,k,2) p(i,j,k,9)=p(i,j,k,9)-pp(i,j,k,1) 382 continue c do 381 m=1,nb do 381 j=1,ny1 do 381 i=1,nx1 p(i,j,k,m)=0.125*(f(i,j,k,m)+f(i+1,j,k,m) 1 +f(i,j+1,k,m)+f(i+1,j+1,k,m) 2 +f(i,j,k+1,m)+f(i+1,j,k+1,m) 3 +f(i,j+1,k+1,m)+f(i+1,j+1,k+1,m)) c u(i,j,k,m)=p(i,j,k,m) 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 y=0.5*hy*float(2*j-2) do 502 m=1,nb do 502 i=1,nx1 ff(i,j,k,m)=u(i,j,k,m) 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) c u(i,j,k,8)=u(i,j,k,8) 1 +dx1*(f(i+1,j+1,k+1,4)*f(i+1,j+1,k+1,6) 1 -f(i+1,j+1,k+1,2)*f(i+1,j+1,k+1,8) 2 +f(i+1,j,k+1,4)*f(i+1,j,k+1,6)-f(i+1,j,k+1,2)*f(i+1,j,k+1,8) 3 +f(i+1,j+1,k,4)*f(i+1,j+1,k,6)-f(i+1,j+1,k,2)*f(i+1,j+1,k,8) 4 +f(i+1,j,k,4)*f(i+1,j,k,6)-f(i+1,j,k,2)*f(i+1,j,k,8) 5 -f(i,j+1,k+1,4)*f(i,j+1,k+1,6)+f(i,j+1,k+1,2)*f(i,j+1,k+1,8) 6 -f(i,j,k+1,4)*f(i,j,k+1,6)+f(i,j,k+1,2)*f(i,j,k+1,8) 7 -f(i,j+1,k,4)*f(i,j+1,k,6)+f(i,j+1,k,2)*f(i,j+1,k,8) 8 -f(i,j,k,4)*f(i,j,k,6)+f(i,j,k,2)*f(i,j,k,8)) 1 -dy1*(f(i+1,j+1,k+1,3)*f(i+1,j+1,k+1,8) 1 -f(i+1,j+1,k+1,4)*f(i+1,j+1,k+1,7) 2 -f(i+1,j,k+1,3)*f(i+1,j,k+1,8)+f(i+1,j,k+1,4)*f(i+1,j,k+1,7) 3 +f(i+1,j+1,k,3)*f(i+1,j+1,k,8)-f(i+1,j+1,k,4)*f(i+1,j+1,k,7) 4 -f(i+1,j,k,3)*f(i+1,j,k,8)+f(i+1,j,k,4)*f(i+1,j,k,7) 5 +f(i,j+1,k+1,3)*f(i,j+1,k+1,8)-f(i,j+1,k+1,4)*f(i,j+1,k+1,7) 6 -f(i,j,k+1,3)*f(i,j,k+1,8)+f(i,j,k+1,4)*f(i,j,k+1,7) 7 +f(i,j+1,k,3)*f(i,j+1,k,8)-f(i,j+1,k,4)*f(i,j+1,k,7) 8 -f(i,j,k,3)*f(i,j,k,8)+f(i,j,k,4)*f(i,j,k,7)) c u(i,j,k,7)=u(i,j,k,7) 1 +dz1*(f(i+1,j+1,k+1,3)*f(i+1,j+1,k+1,8) 1 -f(i+1,j+1,k+1,4)*f(i+1,j+1,k+1,7) 2 +f(i+1,j,k+1,3)*f(i+1,j,k+1,8)-f(i+1,j,k+1,4)*f(i+1,j,k+1,7) 3 -f(i+1,j+1,k,3)*f(i+1,j+1,k,8)+f(i+1,j+1,k,4)*f(i+1,j+1,k,7) 4 -f(i+1,j,k,3)*f(i+1,j,k,8)+f(i+1,j,k,4)*f(i+1,j,k,7) 5 +f(i,j+1,k+1,3)*f(i,j+1,k+1,8)-f(i,j+1,k+1,4)*f(i,j+1,k+1,7) 6 +f(i,j,k+1,3)*f(i,j,k+1,8)-f(i,j,k+1,4)*f(i,j,k+1,7) 7 -f(i,j+1,k,3)*f(i,j+1,k,8)+f(i,j+1,k,4)*f(i,j+1,k,7) 8 -f(i,j,k,3)*f(i,j,k,8)+f(i,j,k,4)*f(i,j,k,7)) 1 -dx1*(f(i+1,j+1,k+1,2)*f(i+1,j+1,k+1,7) 1 -f(i+1,j+1,k+1,3)*f(i+1,j+1,k+1,6) 2 +f(i+1,j,k+1,2)*f(i+1,j,k+1,7)-f(i+1,j,k+1,3)*f(i+1,j,k+1,6) 3 +f(i+1,j+1,k,2)*f(i+1,j+1,k,7)-f(i+1,j+1,k,3)*f(i+1,j+1,k,6) 4 +f(i+1,j,k,2)*f(i+1,j,k,7)-f(i+1,j,k,3)*f(i+1,j,k,6) 5 -f(i,j+1,k+1,2)*f(i,j+1,k+1,7)+f(i,j+1,k+1,3)*f(i,j+1,k+1,6) 6 -f(i,j,k+1,2)*f(i,j,k+1,7)+f(i,j,k+1,3)*f(i,j,k+1,6) 7 -f(i,j+1,k,2)*f(i,j+1,k,7)+f(i,j+1,k,3)*f(i,j+1,k,6) 8 -f(i,j,k,2)*f(i,j,k,7)+f(i,j,k,3)*f(i,j,k,6)) c u(i,j,k,6)=u(i,j,k,6) 1 +dy1*(f(i+1,j+1,k+1,2)*f(i+1,j+1,k+1,7) 1 -f(i+1,j+1,k+1,3)*f(i+1,j+1,k+1,6) 2 -f(i+1,j,k+1,2)*f(i+1,j,k+1,7)+f(i+1,j,k+1,3)*f(i+1,j,k+1,6) 3 +f(i+1,j+1,k,2)*f(i+1,j+1,k,7)-f(i+1,j+1,k,3)*f(i+1,j+1,k,6) 4 -f(i+1,j,k,2)*f(i+1,j,k,7)+f(i+1,j,k,3)*f(i+1,j,k,6) 5 +f(i,j+1,k+1,2)*f(i,j+1,k+1,7)-f(i,j+1,k+1,3)*f(i,j+1,k+1,6) 6 -f(i,j,k+1,2)*f(i,j,k+1,7)+f(i,j,k+1,3)*f(i,j,k+1,6) 7 +f(i,j+1,k,2)*f(i,j+1,k,7)-f(i,j+1,k,3)*f(i,j+1,k,6) 8 -f(i,j,k,2)*f(i,j,k,7)+f(i,j,k,3)*f(i,j,k,6)) 1 -dz1*(f(i+1,j+1,k+1,4)*f(i+1,j+1,k+1,6) 1 -f(i+1,j+1,k+1,2)*f(i+1,j+1,k+1,8) 2 +f(i+1,j,k+1,4)*f(i+1,j,k+1,6)-f(i+1,j,k+1,2)*f(i+1,j,k+1,8) 3 -f(i+1,j+1,k,4)*f(i+1,j+1,k,6)+f(i+1,j+1,k,2)*f(i+1,j+1,k,8) 4 -f(i+1,j,k,4)*f(i+1,j,k,6)+f(i+1,j,k,2)*f(i+1,j,k,8) 5 +f(i,j+1,k+1,4)*f(i,j+1,k+1,6)-f(i,j+1,k+1,2)*f(i,j+1,k+1,8) 6 +f(i,j,k+1,4)*f(i,j,k+1,6)-f(i,j,k+1,2)*f(i,j,k+1,8) 7 -f(i,j+1,k,4)*f(i,j+1,k,6)+f(i,j+1,k,2)*f(i,j+1,k,8) 8 -f(i,j,k,4)*f(i,j,k,6)+f(i,j,k,2)*f(i,j,k,8)) c x3=dx1*(f(i+1,j+1,k+1,2)+f(i+1,j,k+1,2) 1 +f(i+1,j+1,k,2)+f(i+1,j,k,2) 1 -f(i,j+1,k+1,2)-f(i,j,k+1,2)-f(i,j+1,k,2)-f(i,j,k,2)) 2 +dy1*(f(i+1,j+1,k+1,3)-f(i+1,j,k+1,3) 2 +f(i+1,j+1,k,3)-f(i+1,j,k,3) 3 +f(i,j+1,k+1,3)-f(i,j,k+1,3)+f(i,j+1,k,3)-f(i,j,k,3)) 4 +dz1*(f(i+1,j+1,k+1,4)+f(i+1,j,k+1,4) 4 -f(i+1,j+1,k,4)-f(i+1,j,k,4) 5 +f(i,j+1,k+1,4)+f(i,j,k+1,4)-f(i,j+1,k,4)-f(i,j,k,4)) c u(i,j,k,5)=u(i,j,k,5)-p(i,j,k,5)*gam*x3 1 -dx1*p(i,j,k,2)*(f(i+1,j+1,k+1,5)+f(i+1,j,k+1,5) 1 +f(i+1,j+1,k,5)+f(i+1,j,k,5) 2 -f(i,j+1,k+1,5)-f(i,j,k+1,5)-f(i,j+1,k,5)-f(i,j,k,5)) 3 -dy1*p(i,j,k,3)*(f(i+1,j+1,k+1,5)-f(i+1,j,k+1,5) 3 +f(i+1,j+1,k,5)-f(i+1,j,k,5) 4 +f(i,j+1,k+1,5)-f(i,j,k+1,5)+f(i,j+1,k,5)-f(i,j,k,5)) 5 -dz1*p(i,j,k,4)*(f(i+1,j+1,k+1,5)+f(i+1,j,k+1,5) 5 -f(i+1,j+1,k,5)-f(i+1,j,k,5) 6 +f(i,j+1,k+1,5)+f(i,j,k+1,5)-f(i,j+1,k,5)-f(i,j,k,5)) c c u(i,j,k,4)=u(i,j,k,4)+t1*((p(i,j,k,9)*p(i,j,k,7) 1 -p(i,j,k,10)*p(i,j,k,6))/p(i,j,k,1) 1 -0.0*eud*p(i,j,k,4)-gra*z/ro3) 2 -dx1*p(i,j,k,2)*(f(i+1,j+1,k+1,4)+f(i+1,j,k+1,4) 2 +f(i+1,j+1,k,4)+f(i+1,j,k,4) 3 -f(i,j+1,k+1,4)-f(i,j,k+1,4)-f(i,j+1,k,4)-f(i,j,k,4)) 4 -dy1*p(i,j,k,3)*(f(i+1,j+1,k+1,4)-f(i+1,j,k+1,4) 4 +f(i+1,j+1,k,4)-f(i+1,j,k,4) 5 +f(i,j+1,k+1,4)-f(i,j,k+1,4)+f(i,j+1,k,4)-f(i,j,k,4)) 6 -dz1*p(i,j,k,4)*(f(i+1,j+1,k+1,4)+f(i+1,j,k+1,4) 6 -f(i+1,j+1,k,4)-f(i+1,j,k,4) 7 +f(i,j+1,k+1,4)+f(i,j,k+1,4)-f(i,j+1,k,4)-f(i,j,k,4)) 8 -dz1*(f(i+1,j+1,k+1,5)+f(i+1,j,k+1,5) 8 -f(i+1,j+1,k,5)-f(i+1,j,k,5) 9 +f(i,j+1,k+1,5)+f(i,j,k+1,5)-f(i,j+1,k,5)-f(i,j,k,5)) 9 /p(i,j,k,1) c u(i,j,k,3)=u(i,j,k,3)+t1*((p(i,j,k,11)*p(i,j,k,6) 1 -p(i,j,k,9)*p(i,j,k,8))/p(i,j,k,1) 1 -0.0*eud*p(i,j,k,3)-gra*y/ro3) 2 -dx1*p(i,j,k,2)*(f(i+1,j+1,k+1,3)+f(i+1,j,k+1,3) 2 +f(i+1,j+1,k,3)+f(i+1,j,k,3) 3 -f(i,j+1,k+1,3)-f(i,j,k+1,3)-f(i,j+1,k,3)-f(i,j,k,3)) 4 -dy1*p(i,j,k,3)*(f(i+1,j+1,k+1,3)-f(i+1,j,k+1,3) 4 +f(i+1,j+1,k,3)-f(i+1,j,k,3) 5 +f(i,j+1,k+1,3)-f(i,j,k+1,3)+f(i,j+1,k,3)-f(i,j,k,3)) 6 -dz1*p(i,j,k,4)*(f(i+1,j+1,k+1,3)+f(i+1,j,k+1,3) 6 -f(i+1,j+1,k,3)-f(i+1,j,k,3) 7 +f(i,j+1,k+1,3)+f(i,j,k+1,3)-f(i,j+1,k,3)-f(i,j,k,3)) 8 -dy1*(f(i+1,j+1,k+1,5)-f(i+1,j,k+1,5) 8 +f(i+1,j+1,k,5)-f(i+1,j,k,5) 9 +f(i,j+1,k+1,5)-f(i,j,k+1,5)+f(i,j+1,k,5)-f(i,j,k,5)) 9 /p(i,j,k,1) c c uxd=amax1((x-4.0*ar1)/(4.0*ar1),0.0) uxd=amax1((x-xta)/xtd,0.0) uxd=amin1(uxd,1.0) uxd=amin1(uxd*p(i,j,k,2),0.0) c u(i,j,k,2)=u(i,j,k,2)+t1*((p(i,j,k,10)*p(i,j,k,8) 1 -p(i,j,k,11)*p(i,j,k,7))/p(i,j,k,1) 1 -eud*uxd-gra*x/ro3) 2 -dx1*p(i,j,k,2)*(f(i+1,j+1,k+1,2)+f(i+1,j,k+1,2) 2 +f(i+1,j+1,k,2)+f(i+1,j,k,2) 3 -f(i,j+1,k+1,2)-f(i,j,k+1,2)-f(i,j+1,k,2)-f(i,j,k,2)) 4 -dy1*p(i,j,k,3)*(f(i+1,j+1,k+1,2)-f(i+1,j,k+1,2) 4 +f(i+1,j+1,k,2)-f(i+1,j,k,2) 5 +f(i,j+1,k+1,2)-f(i,j,k+1,2)+f(i,j+1,k,2)-f(i,j,k,2)) 6 -dz1*p(i,j,k,4)*(f(i+1,j+1,k+1,2)+f(i+1,j,k+1,2) 6 -f(i+1,j+1,k,2)-f(i+1,j,k,2) 7 +f(i,j+1,k+1,2)+f(i,j,k+1,2)-f(i,j+1,k,2)-f(i,j,k,2)) 8 -dx1*(f(i+1,j+1,k+1,5)+f(i+1,j,k+1,5) 8 +f(i+1,j+1,k,5)+f(i+1,j,k,5) 9 -f(i,j+1,k+1,5)-f(i,j,k+1,5)-f(i,j+1,k,5)-f(i,j,k,5)) 9 /p(i,j,k,1) c u(i,j,k,1)=u(i,j,k,1) 1 -dx1*(f(i+1,j+1,k+1,2)*f(i+1,j+1,k+1,1) 1 +f(i+1,j,k+1,2)*f(i+1,j,k+1,1) 2 +f(i+1,j+1,k,2)*f(i+1,j+1,k,1)+f(i+1,j,k,2)*f(i+1,j,k,1) 3 -f(i,j+1,k+1,2)*f(i,j+1,k+1,1)-f(i,j,k+1,2)*f(i,j,k+1,1) 4 -f(i,j+1,k,2)*f(i,j+1,k,1)-f(i,j,k,2)*f(i,j,k,1)) 5 -dy1*(f(i+1,j+1,k+1,3)*f(i+1,j+1,k+1,1) 5 -f(i+1,j,k+1,3)*f(i+1,j,k+1,1) 6 +f(i+1,j+1,k,3)*f(i+1,j+1,k,1)-f(i+1,j,k,3)*f(i+1,j,k,1) 7 +f(i,j+1,k+1,3)*f(i,j+1,k+1,1)-f(i,j,k+1,3)*f(i,j,k+1,1) 8 +f(i,j+1,k,3)*f(i,j+1,k,1)-f(i,j,k,3)*f(i,j,k,1)) 1 -dz1*(f(i+1,j+1,k+1,4)*f(i+1,j+1,k+1,1) 1 +f(i+1,j,k+1,4)*f(i+1,j,k+1,1) 2 -f(i+1,j+1,k,4)*f(i+1,j+1,k,1)-f(i+1,j,k,4)*f(i+1,j,k,1) 3 +f(i,j+1,k+1,4)*f(i,j+1,k+1,1)+f(i,j,k+1,4)*f(i,j,k+1,1) 4 -f(i,j+1,k,4)*f(i,j+1,k,1)-f(i,j,k,4)*f(i,j,k,1)) 501 continue c do 520 i=1,nx1 x=hhx(i)+xx4 ro2=x*x+y*y+z*z x2=amax1(ro2/ar2-1.0,0.0) x2=aru*x2*x2 x2=x2/(x2+1.0) u(i,j,k,8)=u(i,j,k,8)*x2+ff(i,j,k,8)*(1.0-x2) u(i,j,k,7)=u(i,j,k,7)*x2+ff(i,j,k,7)*(1.0-x2) u(i,j,k,6)=u(i,j,k,6)*x2+ff(i,j,k,6)*(1.0-x2) u(i,j,k,5)=u(i,j,k,5)*x2+ff(i,j,k,5)*(1.0-x2) u(i,j,k,4)=u(i,j,k,4)*x2+ff(i,j,k,4)*(1.0-x2) u(i,j,k,3)=u(i,j,k,3)*x2+ff(i,j,k,3)*(1.0-x2) u(i,j,k,2)=u(i,j,k,2)*x2+ff(i,j,k,2)*(1.0-x2) u(i,j,k,1)=u(i,j,k,1)*x2+ff(i,j,k,1)*(1.0-x2) 520 continue c do 524 i=1,nx1 x5=u(i,j,k,5) x1=u(i,j,k,1) y5=amax1(x5,pr02) y1=amax1(x1,ro02) c x2=u(i,j,k,2)*u(i,j,k,2)+u(i,j,k,3)*u(i,j,k,3)+u(i,j,k,4)*u(i,j,k,4) c vmax=amax1(vmax,x2) x=amax1(x1,0.0)/y1 u(i,j,k,5)=y5 u(i,j,k,4)=x*u(i,j,k,4) u(i,j,k,3)=x*u(i,j,k,3) u(i,j,k,2)=x*u(i,j,k,2) u(i,j,k,1)=y1 524 continue 50 continue !hpf$ end on 90 continue c c c c do 70 m=1,nb !hpf$ independent,new(i,j,m) do 72 k=1,nz2 !hpf$ on home(ff(:,:,k,:)),local(f,ff,i,j,m) begin do 7211 m=1,nb do 7211 j=1,ny2 do 7211 i=1,nx2 ff(i,j,k,m)=f(i,j,k,m) 7211 continue !hpf$ end on 72 continue c 70 continue c !hpf$ reflect u !hpf$ reflect ff !hpf$ reflect pp c c second step c c !hpf$ independent,new(i,j,k,m),reduction(max:vmax) do 190 k=2,nz1 !hpf$ on home(f(:,:,k,:)),local(f,u,p,pp,ff,i,j,m, !hpf$* x,y,z,hx,hy,hz,gam,pmu, !hpf$* xx4,ar1,ar2,ro2,ro3,hhx,x1,x2,x3,x5,eat,eat0,po0, !hpf$* dx2,dy2,dz2,dx4,dy4,dz4,eud,gra,rmu,uxd, !hpf$* ro02,pr02,y1,y5,vmax,aru) begin c c c current do 68 j=2,ny1 do 68 i=2,nx1 p(i,j,k,11)=0.25*((u(i,j,k,7)+u(i,j-1,k,7) 1 +u(i,j,k-1,7)+u(i,j-1,k-1,7) 2 -u(i-1,j,k,7)-u(i-1,j-1,k,7) 3 -u(i-1,j,k-1,7)-u(i-1,j-1,k-1,7))/hx 4 -(u(i,j,k,6)-u(i,j-1,k,6) 5 +u(i,j,k-1,6)-u(i,j-1,k-1,6) 6 +u(i-1,j,k,6)-u(i-1,j-1,k,6) 7 +u(i-1,j,k-1,6)-u(i-1,j-1,k-1,6))/hy) p(i,j,k,10)=0.25*((u(i,j,k,6)+u(i,j-1,k,6) 1 -u(i,j,k-1,6)-u(i,j-1,k-1,6) 2 +u(i-1,j,k,6)+u(i-1,j-1,k,6) 3 -u(i-1,j,k-1,6)-u(i-1,j-1,k-1,6))/hz 4 -(u(i,j,k,8)+u(i,j-1,k,8) 5 +u(i,j,k-1,8)+u(i,j-1,k-1,8) 6 -u(i-1,j,k,8)-u(i-1,j-1,k,8) 7 -u(i-1,j,k-1,8)-u(i-1,j-1,k-1,8))/hx) p(i,j,k,9)=0.25*((u(i,j,k,8)-u(i,j-1,k,8) 1 +u(i,j,k-1,8)-u(i,j-1,k-1,8) 2 +u(i-1,j,k,8)-u(i-1,j-1,k,8) 3 +u(i-1,j,k-1,8)-u(i-1,j-1,k-1,8))/hy 4 -(u(i,j,k,7)+u(i,j-1,k,7) 5 -u(i,j,k-1,7)-u(i,j-1,k-1,7) 6 +u(i-1,j,k,7)+u(i-1,j-1,k,7) 7 -u(i-1,j,k-1,7)-u(i-1,j-1,k-1,7))/hz) 68 continue do 682 j=2,ny1 do 682 i=2,nx1 p(i,j,k,11)=p(i,j,k,11)-0.125*( 1 pp(i,j,k,3)+pp(i-1,j,k,3) 2 +pp(i,j-1,k,3)+pp(i-1,j-1,k,3) 3 +pp(i,j,k-1,3)+pp(i-1,j,k-1,3) 4 +pp(i,j-1,k-1,3)+pp(i-1,j-1,k-1,3)) p(i,j,k,10)=p(i,j,k,10)-0.125*( 1 pp(i,j,k,2)+pp(i-1,j,k,2) 2 +pp(i,j-1,k,2)+pp(i-1,j-1,k,2) 3 +pp(i,j,k-1,2)+pp(i-1,j,k-1,2) 4 +pp(i,j-1,k-1,2)+pp(i-1,j-1,k-1,2)) p(i,j,k,9)=p(i,j,k,9)-0.125*( 1 pp(i,j,k,1)+pp(i-1,j,k,1) 2 +pp(i,j-1,k,1)+pp(i-1,j-1,k,1) 3 +pp(i,j,k-1,1)+pp(i-1,j,k-1,1) 4 +pp(i,j-1,k-1,1)+pp(i-1,j-1,k-1,1)) 682 continue c do 681 m=1,nb do 681 j=2,ny1 do 681 i=2,nx1 p(i,j,k,m)=0.125*(u(i,j,k,m)+u(i-1,j,k,m) 1 +u(i,j-1,k,m)+u(i-1,j-1,k,m) 2 +u(i,j,k-1,m)+u(i-1,j,k-1,m) 3 +u(i,j-1,k-1,m)+u(i-1,j-1,k-1,m)) 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 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) c x1=ff(i,j,k,5)/(ff(i,j,k,1)*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(i+1,j,k,8)-2.0*ff(i,j,k,8)+ff(i-1,j,k,8)) 1 +dy4*(ff(i,j+1,k,8)-2.0*ff(i,j,k,8)+ff(i,j-1,k,8)) 2 +dz4*(ff(i,j,k+1,8)-2.0*ff(i,j,k,8)+ff(i,j,k-1,8))) f(i,j,k,8)=ff(i,j,k,8)+x1 1 +dx2*(u(i,j,k,4)*u(i,j,k,6)-u(i,j,k,2)*u(i,j,k,8) 2 +u(i,j-1,k,4)*u(i,j-1,k,6)-u(i,j-1,k,2)*u(i,j-1,k,8) 3 +u(i,j,k-1,4)*u(i,j,k-1,6)-u(i,j,k-1,2)*u(i,j,k-1,8) 4 +u(i,j-1,k-1,4)*u(i,j-1,k-1,6)-u(i,j-1,k-1,2)*u(i,j-1,k-1,8) 5 -u(i-1,j,k,4)*u(i-1,j,k,6)+u(i-1,j,k,2)*u(i-1,j,k,8) 6 -u(i-1,j-1,k,4)*u(i-1,j-1,k,6)+u(i-1,j-1,k,2)*u(i-1,j-1,k,8) 7 -u(i-1,j,k-1,4)*u(i-1,j,k-1,6)+u(i-1,j,k-1,2)*u(i-1,j,k-1,8) 8 -u(i-1,j-1,k-1,4)*u(i-1,j-1,k-1,6) 8 +u(i-1,j-1,k-1,2)*u(i-1,j-1,k-1,8)) 1 -dy2*(u(i,j,k,3)*u(i,j,k,8)-u(i,j,k,4)*u(i,j,k,7) 2 -u(i,j-1,k,3)*u(i,j-1,k,8)+u(i,j-1,k,4)*u(i,j-1,k,7) 3 +u(i,j,k-1,3)*u(i,j,k-1,8)-u(i,j,k-1,4)*u(i,j,k-1,7) 4 -u(i,j-1,k-1,3)*u(i,j-1,k-1,8)+u(i,j-1,k-1,4)*u(i,j-1,k-1,7) 5 +u(i-1,j,k,3)*u(i-1,j,k,8)-u(i-1,j,k,4)*u(i-1,j,k,7) 6 -u(i-1,j-1,k,3)*u(i-1,j-1,k,8)+u(i-1,j-1,k,4)*u(i-1,j-1,k,7) 7 +u(i-1,j,k-1,3)*u(i-1,j,k-1,8)-u(i-1,j,k-1,4)*u(i-1,j,k-1,7) 8 -u(i-1,j-1,k-1,3)*u(i-1,j-1,k-1,8) 8 +u(i-1,j-1,k-1,4)*u(i-1,j-1,k-1,7)) c x1=eat*(dx4*(ff(i+1,j,k,7)-2.0*ff(i,j,k,7)+ff(i-1,j,k,7)) 1 +dy4*(ff(i,j+1,k,7)-2.0*ff(i,j,k,7)+ff(i,j-1,k,7)) 2 +dz4*(ff(i,j,k+1,7)-2.0*ff(i,j,k,7)+ff(i,j,k-1,7))) f(i,j,k,7)=ff(i,j,k,7)+x1 1 +dz2*(u(i,j,k,3)*u(i,j,k,8)-u(i,j,k,4)*u(i,j,k,7) 2 +u(i,j-1,k,3)*u(i,j-1,k,8)-u(i,j-1,k,4)*u(i,j-1,k,7) 3 -u(i,j,k-1,3)*u(i,j,k-1,8)+u(i,j,k-1,4)*u(i,j,k-1,7) 4 -u(i,j-1,k-1,3)*u(i,j-1,k-1,8)+u(i,j-1,k-1,4)*u(i,j-1,k-1,7) 5 +u(i-1,j,k,3)*u(i-1,j,k,8)-u(i-1,j,k,4)*u(i-1,j,k,7) 6 +u(i-1,j-1,k,3)*u(i-1,j-1,k,8)-u(i-1,j-1,k,4)*u(i-1,j-1,k,7) 7 -u(i-1,j,k-1,3)*u(i-1,j,k-1,8)+u(i-1,j,k-1,4)*u(i-1,j,k-1,7) 8 -u(i-1,j-1,k-1,3)*u(i-1,j-1,k-1,8) 8 +u(i-1,j-1,k-1,4)*u(i-1,j-1,k-1,7)) 1 -dx2*(u(i,j,k,2)*u(i,j,k,7)-u(i,j,k,3)*u(i,j,k,6) 2 +u(i,j-1,k,2)*u(i,j-1,k,7)-u(i,j-1,k,3)*u(i,j-1,k,6) 3 +u(i,j,k-1,2)*u(i,j,k-1,7)-u(i,j,k-1,3)*u(i,j,k-1,6) 4 +u(i,j-1,k-1,2)*u(i,j-1,k-1,7)-u(i,j-1,k-1,3)*u(i,j-1,k-1,6) 5 -u(i-1,j,k,2)*u(i-1,j,k,7)+u(i-1,j,k,3)*u(i-1,j,k,6) 6 -u(i-1,j-1,k,2)*u(i-1,j-1,k,7)+u(i-1,j-1,k,3)*u(i-1,j-1,k,6) 7 -u(i-1,j,k-1,2)*u(i-1,j,k-1,7)+u(i-1,j,k-1,3)*u(i-1,j,k-1,6) 8 -u(i-1,j-1,k-1,2)*u(i-1,j-1,k-1,7) 8 +u(i-1,j-1,k-1,3)*u(i-1,j-1,k-1,6)) c x1=eat*(dx4*(ff(i+1,j,k,6)-2.0*ff(i,j,k,6)+ff(i-1,j,k,6)) 1 +dy4*(ff(i,j+1,k,6)-2.0*ff(i,j,k,6)+ff(i,j-1,k,6)) 2 +dz4*(ff(i,j,k+1,6)-2.0*ff(i,j,k,6)+ff(i,j,k-1,6))) f(i,j,k,6)=ff(i,j,k,6)+x1 1 +dy2*(u(i,j,k,2)*u(i,j,k,7)-u(i,j,k,3)*u(i,j,k,6) 2 -u(i,j-1,k,2)*u(i,j-1,k,7)+u(i,j-1,k,3)*u(i,j-1,k,6) 3 +u(i,j,k-1,2)*u(i,j,k-1,7)-u(i,j,k-1,3)*u(i,j,k-1,6) 4 -u(i,j-1,k-1,2)*u(i,j-1,k-1,7)+u(i,j-1,k-1,3)*u(i,j-1,k-1,6) 5 +u(i-1,j,k,2)*u(i-1,j,k,7)-u(i-1,j,k,3)*u(i-1,j,k,6) 6 -u(i-1,j-1,k,2)*u(i-1,j-1,k,7)+u(i-1,j-1,k,3)*u(i-1,j-1,k,6) 7 +u(i-1,j,k-1,2)*u(i-1,j,k-1,7)-u(i-1,j,k-1,3)*u(i-1,j,k-1,6) 8 -u(i-1,j-1,k-1,2)*u(i-1,j-1,k-1,7) 8 +u(i-1,j-1,k-1,3)*u(i-1,j-1,k-1,6)) 1 -dz2*(u(i,j,k,4)*u(i,j,k,6)-u(i,j,k,2)*u(i,j,k,8) 2 +u(i,j-1,k,4)*u(i,j-1,k,6)-u(i,j-1,k,2)*u(i,j-1,k,8) 3 -u(i,j,k-1,4)*u(i,j,k-1,6)+u(i,j,k-1,2)*u(i,j,k-1,8) 4 -u(i,j-1,k-1,4)*u(i,j-1,k-1,6)+u(i,j-1,k-1,2)*u(i,j-1,k-1,8) 5 +u(i-1,j,k,4)*u(i-1,j,k,6)-u(i-1,j,k,2)*u(i-1,j,k,8) 6 +u(i-1,j-1,k,4)*u(i-1,j-1,k,6)-u(i-1,j-1,k,2)*u(i-1,j-1,k,8) 7 -u(i-1,j,k-1,4)*u(i-1,j,k-1,6)+u(i-1,j,k-1,2)*u(i-1,j,k-1,8) 8 -u(i-1,j-1,k-1,4)*u(i-1,j-1,k-1,6) 8 +u(i-1,j-1,k-1,2)*u(i-1,j-1,k-1,8)) c x3=dx2*(u(i,j,k,2)+u(i,j-1,k,2)+u(i,j,k-1,2)+u(i,j-1,k-1,2) 1 -u(i-1,j,k,2)-u(i-1,j-1,k,2)-u(i-1,j,k-1,2)-u(i-1,j-1,k-1,2)) 2 +dy2*(u(i,j,k,3)-u(i,j-1,k,3)+u(i,j,k-1,3)-u(i,j-1,k-1,3) 3 +u(i-1,j,k,3)-u(i-1,j-1,k,3)+u(i-1,j,k-1,3)-u(i-1,j-1,k-1,3)) 4 +dz2*(u(i,j,k,4)+u(i,j-1,k,4)-u(i,j,k-1,4)-u(i,j-1,k-1,4) 5 +u(i-1,j,k,4)+u(i-1,j-1,k,4)-u(i-1,j,k-1,4)-u(i-1,j-1,k-1,4)) f(i,j,k,5)=ff(i,j,k,5)-p(i,j,k,5)*gam*x3 1 -dx2*p(i,j,k,2)* 1 (u(i,j,k,5)+u(i,j-1,k,5)+u(i,j,k-1,5)+u(i,j-1,k-1,5) 2 -u(i-1,j,k,5)-u(i-1,j-1,k,5)-u(i-1,j,k-1,5)-u(i-1,j-1,k-1,5)) 3 -dy2*p(i,j,k,3)* 3 (u(i,j,k,5)-u(i,j-1,k,5)+u(i,j,k-1,5)-u(i,j-1,k-1,5) 4 +u(i-1,j,k,5)-u(i-1,j-1,k,5)+u(i-1,j,k-1,5)-u(i-1,j-1,k-1,5)) 5 -dz2*p(i,j,k,4)* 5 (u(i,j,k,5)+u(i,j-1,k,5)-u(i,j,k-1,5)-u(i,j-1,k-1,5) 6 +u(i-1,j,k,5)+u(i-1,j-1,k,5)-u(i-1,j,k-1,5)-u(i-1,j-1,k-1,5)) 7 +pmu*(dx4*(ff(i+1,j,k,5)-2.0*ff(i,j,k,5)+ff(i-1,j,k,5)) 8 +dy4*(ff(i,j+1,k,5)-2.0*ff(i,j,k,5)+ff(i,j-1,k,5)) 9 +dz4*(ff(i,j,k+1,5)-2.0*ff(i,j,k,5)+ff(i,j,k-1,5))) c f(i,j,k,4)=ff(i,j,k,4)+t*((p(i,j,k,9)*p(i,j,k,7) 1 -p(i,j,k,10)*p(i,j,k,6))/p(i,j,k,1) 1 -0.0*eud*p(i,j,k,4)-gra*z/ro3) 2 -dx2*p(i,j,k,2)* 2 (u(i,j,k,4)+u(i,j-1,k,4)+u(i,j,k-1,4)+u(i,j-1,k-1,4) 3 -u(i-1,j,k,4)-u(i-1,j-1,k,4)-u(i-1,j,k-1,4)-u(i-1,j-1,k-1,4)) 4 -dy2*p(i,j,k,3)* 4 (u(i,j,k,4)-u(i,j-1,k,4)+u(i,j,k-1,4)-u(i,j-1,k-1,4) 5 +u(i-1,j,k,4)-u(i-1,j-1,k,4)+u(i-1,j,k-1,4)-u(i-1,j-1,k-1,4)) 6 -dz2*p(i,j,k,4)* 6 (u(i,j,k,4)+u(i,j-1,k,4)-u(i,j,k-1,4)-u(i,j-1,k-1,4) 7 +u(i-1,j,k,4)+u(i-1,j-1,k,4)-u(i-1,j,k-1,4)-u(i-1,j-1,k-1,4)) 8 -dz2*(u(i,j,k,5)+u(i,j-1,k,5)-u(i,j,k-1,5)-u(i,j-1,k-1,5) 9 +u(i-1,j,k,5)+u(i-1,j-1,k,5)-u(i-1,j,k-1,5)-u(i-1,j-1,k-1,5)) 9 /p(i,j,k,1) 1 +rmu*(dx4*(ff(i+1,j,k,4)-2.0*ff(i,j,k,4)+ff(i-1,j,k,4)) 2 +dy4*(ff(i,j+1,k,4)-2.0*ff(i,j,k,4)+ff(i,j-1,k,4)) 3 +dz4*(ff(i,j,k+1,4)-2.0*ff(i,j,k,4)+ff(i,j,k-1,4)))/ff(i,j,k,1) c f(i,j,k,3)=ff(i,j,k,3)+t*((p(i,j,k,11)*p(i,j,k,6) 1 -p(i,j,k,9)*p(i,j,k,8))/p(i,j,k,1) 1 -0.0*eud*p(i,j,k,3)-gra*y/ro3) 2 -dx2*p(i,j,k,2)* 2 (u(i,j,k,3)+u(i,j-1,k,3)+u(i,j,k-1,3)+u(i,j-1,k-1,3) 3 -u(i-1,j,k,3)-u(i-1,j-1,k,3)-u(i-1,j,k-1,3)-u(i-1,j-1,k-1,3)) 4 -dy2*p(i,j,k,3)* 4 (u(i,j,k,3)-u(i,j-1,k,3)+u(i,j,k-1,3)-u(i,j-1,k-1,3) 5 +u(i-1,j,k,3)-u(i-1,j-1,k,3)+u(i-1,j,k-1,3)-u(i-1,j-1,k-1,3)) 6 -dz2*p(i,j,k,4)* 6 (u(i,j,k,3)+u(i,j-1,k,3)-u(i,j,k-1,3)-u(i,j-1,k-1,3) 7 +u(i-1,j,k,3)+u(i-1,j-1,k,3)-u(i-1,j,k-1,3)-u(i-1,j-1,k-1,3)) 8 -dy2*(u(i,j,k,5)-u(i,j-1,k,5)+u(i,j,k-1,5)-u(i,j-1,k-1,5) 9 +u(i-1,j,k,5)-u(i-1,j-1,k,5)+u(i-1,j,k-1,5)-u(i-1,j-1,k-1,5)) 9 /p(i,j,k,1) 1 +rmu*(dx4*(ff(i+1,j,k,3)-2.0*ff(i,j,k,3)+ff(i-1,j,k,3)) 2 +dy4*(ff(i,j+1,k,3)-2.0*ff(i,j,k,3)+ff(i,j-1,k,3)) 3 +dz4*(ff(i,j,k+1,3)-2.0*ff(i,j,k,3)+ff(i,j,k-1,3)))/ff(i,j,k,1) c uxd=amax1((x-4.0*ar1)/(4.0*ar1),0.0) uxd=amax1((x-xta)/xtd,0.0) uxd=amin1(uxd,1.0) uxd=amin1(uxd*p(i,j,k,2),0.0) c f(i,j,k,2)=ff(i,j,k,2)+t*((p(i,j,k,10)*p(i,j,k,8) 1 -p(i,j,k,11)*p(i,j,k,7))/p(i,j,k,1) 1 -eud*uxd-gra*x/ro3) 2 -dx2*p(i,j,k,2)* 2 (u(i,j,k,2)+u(i,j-1,k,2)+u(i,j,k-1,2)+u(i,j-1,k-1,2) 3 -u(i-1,j,k,2)-u(i-1,j-1,k,2)-u(i-1,j,k-1,2)-u(i-1,j-1,k-1,2)) 4 -dy2*p(i,j,k,3)* 4 (u(i,j,k,2)-u(i,j-1,k,2)+u(i,j,k-1,2)-u(i,j-1,k-1,2) 5 +u(i-1,j,k,2)-u(i-1,j-1,k,2)+u(i-1,j,k-1,2)-u(i-1,j-1,k-1,2)) 6 -dz2*p(i,j,k,4)* 6 (u(i,j,k,2)+u(i,j-1,k,2)-u(i,j,k-1,2)-u(i,j-1,k-1,2) 7 +u(i-1,j,k,2)+u(i-1,j-1,k,2)-u(i-1,j,k-1,2)-u(i-1,j-1,k-1,2)) 8 -dx2*(u(i,j,k,5)+u(i,j-1,k,5)+u(i,j,k-1,5)+u(i,j-1,k-1,5) 9 -u(i-1,j,k,5)-u(i-1,j-1,k,5)-u(i-1,j,k-1,5)-u(i-1,j-1,k-1,5)) 9 /p(i,j,k,1) 1 +rmu*(dx4*(ff(i+1,j,k,2)-2.0*ff(i,j,k,2)+ff(i-1,j,k,2)) 2 +dy4*(ff(i,j+1,k,2)-2.0*ff(i,j,k,2)+ff(i,j-1,k,2)) 3 +dz4*(ff(i,j,k+1,2)-2.0*ff(i,j,k,2)+ff(i,j,k-1,2)))/ff(i,j,k,1) c f(i,j,k,1)=ff(i,j,k,1) 1 -dx2*(u(i,j,k,2)*u(i,j,k,1)+u(i,j-1,k,2)*u(i,j-1,k,1) 2 +u(i,j,k-1,2)*u(i,j,k-1,1)+u(i,j-1,k-1,2)*u(i,j-1,k-1,1) 3 -u(i-1,j,k,2)*u(i-1,j,k,1)-u(i-1,j-1,k,2)*u(i-1,j-1,k,1) 4 -u(i-1,j,k-1,2)*u(i-1,j,k-1,1) 4 -u(i-1,j-1,k-1,2)*u(i-1,j-1,k-1,1)) 5 -dy2*(u(i,j,k,3)*u(i,j,k,1)-u(i,j-1,k,3)*u(i,j-1,k,1) 6 +u(i,j,k-1,3)*u(i,j,k-1,1)-u(i,j-1,k-1,3)*u(i,j-1,k-1,1) 7 +u(i-1,j,k,3)*u(i-1,j,k,1)-u(i-1,j-1,k,3)*u(i-1,j-1,k,1) 8 +u(i-1,j,k-1,3)*u(i-1,j,k-1,1) 8 -u(i-1,j-1,k-1,3)*u(i-1,j-1,k-1,1)) 1 -dz2*(u(i,j,k,4)*u(i,j,k,1)+u(i,j-1,k,4)*u(i,j-1,k,1) 2 -u(i,j,k-1,4)*u(i,j,k-1,1)-u(i,j-1,k-1,4)*u(i,j-1,k-1,1) 3 +u(i-1,j,k,4)*u(i-1,j,k,1)+u(i-1,j-1,k,4)*u(i-1,j-1,k,1) 4 -u(i-1,j,k-1,4)*u(i-1,j,k-1,1) 4 -u(i-1,j-1,k-1,4)*u(i-1,j-1,k-1,1)) 5 +dfu*(dx4*(ff(i+1,j,k,1)-2.0*ff(i,j,k,1)+ff(i-1,j,k,1)) 6 +dy4*(ff(i,j+1,k,1)-2.0*ff(i,j,k,1)+ff(i,j-1,k,1)) 7 +dz4*(ff(i,j,k+1,1)-2.0*ff(i,j,k,1)+ff(i,j,k-1,1))) c 801 continue c do 820 i=2,nx1 x=hhx(i)+xx4 ro2=x*x+y*y+z*z x2=amax1(ro2/ar2-1.0,0.0) x2=aru*x2*x2 x2=x2/(x2+1.0) f(i,j,k,8)=f(i,j,k,8)*x2+ff(i,j,k,8)*(1.0-x2) f(i,j,k,7)=f(i,j,k,7)*x2+ff(i,j,k,7)*(1.0-x2) f(i,j,k,6)=f(i,j,k,6)*x2+ff(i,j,k,6)*(1.0-x2) f(i,j,k,5)=f(i,j,k,5)*x2+ff(i,j,k,5)*(1.0-x2) f(i,j,k,4)=f(i,j,k,4)*x2+ff(i,j,k,4)*(1.0-x2) f(i,j,k,3)=f(i,j,k,3)*x2+ff(i,j,k,3)*(1.0-x2) f(i,j,k,2)=f(i,j,k,2)*x2+ff(i,j,k,2)*(1.0-x2) f(i,j,k,1)=f(i,j,k,1)*x2+ff(i,j,k,1)*(1.0-x2) 820 continue c do 824 i=2,nx1 x5=f(i,j,k,5) x1=f(i,j,k,1) y5=amax1(x5,pr02) y1=amax1(x1,ro02) x2=f(i,j,k,2)*f(i,j,k,2)+f(i,j,k,3)*f(i,j,k,3) 1 +f(i,j,k,4)*f(i,j,k,4) vmax=amax1(vmax,x2) x=amax1(x1,0.0)/y1 f(i,j,k,5)=y5 f(i,j,k,4)=x*f(i,j,k,4) f(i,j,k,3)=x*f(i,j,k,3) f(i,j,k,2)=x*f(i,j,k,2) f(i,j,k,1)=y1 824 continue 80 continue c !hpf$ end on 190 continue c !xocl end spread max(vmax) 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 c !hpf$ independent,new(i,j,k) do 102 k=1,nz2 !hpf$ on home(f(:,:,k,:)),local(f,i,j,m,nx11,x,x1,hx,hhx, !hpf$* bis) begin do 1021 j=1,ny2 do 1021 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 f(i,j,k,8)=(1.0-x)*f(i,j,k,8)+bis*x f(i,j,k,7)=(1.0-x)*f(i,j,k,7) f(i,j,k,6)=(1.0-x)*f(i,j,k,6) f(i,j,k,5)=(1.0-x)*f(i,j,k,5)+pr01*x f(i,j,k,4)=(1.0-x)*f(i,j,k,4) f(i,j,k,3)=(1.0-x)*f(i,j,k,3) f(i,j,k,2)=(1.0-x)*f(i,j,k,2)+vsw*x f(i,j,k,1)=(1.0-x)*f(i,j,k,1)+ro01*x 1021 continue !hpf$ end on 102 continue c !hpf$ independent,new(i,j,k) do 104 k=1,nz1 !hpf$ on home(u(:,:,k,:)),local(u,i,j,m,nx11,x,x1,hx,hhx, !hpf$* bis) begin do 1041 j=1,ny1 do 1041 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 u(i,j,k,8)=(1.0-x)*u(i,j,k,8)+bis*x u(i,j,k,7)=(1.0-x)*u(i,j,k,7) u(i,j,k,6)=(1.0-x)*u(i,j,k,6) u(i,j,k,5)=(1.0-x)*u(i,j,k,5)+pr01*x u(i,j,k,4)=(1.0-x)*u(i,j,k,4) u(i,j,k,3)=(1.0-x)*u(i,j,k,3) u(i,j,k,2)=(1.0-x)*u(i,j,k,2)+vsw*x u(i,j,k,1)=(1.0-x)*u(i,j,k,1)+ro01*x 1041 continue !hpf$ end on 104 continue 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 c c!hpf$ asynchronous(id1),nobuffer begin c gf(1:nx2,1:ny2,1:nz2,1:nb)=f(1:nx2,1:ny2,1:nz2,1:nb) c!hpf$ end asynchronous c!hpf$ asyncwait(id1) c c write(ntap) f do 173 m=1,nb do 173 k=1,nz2 c!hpf$ asynchronous(id1),nobuffer begin c gf(1:nx2,1:ny2)=f(1:nx2,1:ny2,k,m) c!hpf$ end asynchronous c!hpf$ asyncwait(id1) c write(ntap) gf write(ntap) f(1:nx2,1:ny2,k,m) 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 !xocl end parallel c rewind 12 stop end subroutine equib7(f,pp,nxp,ro01,pr01,rrat,cj,cp) c implicit real*8 (a-h,o-z) parameter (npe=2) !hpf$ processors pe(npe) parameter (nx=320,ny= 80,nz= 80,nb=8) c parameter (nx=480,ny=160,nz=160,nb=8) c parameter (nx=640,ny=160,nz=160,nb=8) c parameter (nx=800,ny=200,nz=200,nb=8) parameter (nx1=nx+1,nx2=nx+2,ny1=ny+1,ny2=ny+2) parameter (nz1=nz+1,nz2=nz+2) c dimension f(nx2,ny2,nz2,nb),pp(nx2,ny2,nz2,3) c dimension gf(nx2,ny2,nz2,nb),gpp(nx2,ny2,nz2,3) c dimension gf(nx2,ny2,nz2,nb) dimension cj(1),cp(1) !hpf$ distribute f(*,*,block,*) onto pe !hpf$ distribute pp(*,*,block,*) onto pe !hpf$ shadow f(0,0,1:1,0) !hpf$ shadow pp(0,0,1:1,0) c !xocl global gf,gpp c equivalence (gf,f),(gpp,pp) c common /blk/f,pp 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 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 n1=nx2 n2=n1*ny2 n3=n2*nz2 hx=cp(1)/float(nx1) hy=cp(2)/float(ny1) hz=cp(3)/float(nz1) vsw=cp(9) c !hpf$ independent,new(i,j,k, !hpf$* xx,y1,px1,ro1,ro2,ro3,x,y,z,vr1) do 10 k=1,nz2 !hpf$ on home(f(:,:,k,:)),local(f,i,j,m,nxp,x,y,z, !hpf$* dn1,dv1,ax1,hx,hy,hz,dn,vr1,vsw,ro1, !hpf$* ro2,ro3,x1,ro02,cp,cj,p0,ro1,pr02,px1, !hpf$* y1,xx,bis,ro01) begin do 101 j=1,ny2 c vr1=vsw do 101 i=1,nx2 vr1=vsw 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 f(i,j,k,1)=1.0/ro3 c f(i,j,k,1)=ro1**(-x1) if(f(i,j,k,1).lt.ro02) f(i,j,k,1)=ro02 f(i,j,k,5)=1.0e-7*cp(6)/ro2 c f(i,j,k,5)=p0*ro1**(-x2) if(f(i,j,k,5).lt.pr02) f(i,j,k,5)=pr02 f(i,j,k,6)=-2.0*b0*x*z/(ro2*ro2) f(i,j,k,7)=-2.0*b0*y*z/(ro2*ro2) f(i,j,k,8)=b0*(x*x+y*y-z*z)/(ro2*ro2) f(i,j,k,6)=-3.0*b0*x*z/(ro2*ro3) f(i,j,k,7)=-3.0*b0*y*z/(ro2*ro3) f(i,j,k,8)=b0*(x*x+y*y-2.0*z*z)/(ro2*ro3) px1=sqrt(f(i,j,k,6)*f(i,j,k,6)+f(i,j,k,7)*f(i,j,k,7) 1 +f(i,j,k,8)*f(i,j,k,8)) if(abs(vsw).lt.1.0e-8) go to 101 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(i,j,k,2)=vr1 xx=vr1/vsw f(i,j,k,5)=f(i,j,k,5)+(pr01-pr02)*xx f(i,j,k,8)=f(i,j,k,8)+bis 101 continue !hpf$ end on 10 continue c c c current !hpf$ reflect f !hpf$ independent,new(i,j,k) do 38 k=1,nz1 !hpf$ on home(pp(:,:,k,:)),local(pp,f,i,j,m) begin do 381 j=1,ny1 do 381 i=1,nx1 pp(i,j,k,3)=0.25*((f(i+1,j+1,k+1,7)+f(i+1,j,k+1,7) 1 +f(i+1,j+1,k,7)+f(i+1,j,k,7) 2 -f(i,j+1,k+1,7)-f(i,j,k+1,7)-f(i,j+1,k,7)-f(i,j,k,7))/hx 3 -(f(i+1,j+1,k+1,6)-f(i+1,j,k+1,6)+f(i+1,j+1,k,6)-f(i+1,j,k,6) 4 +f(i,j+1,k+1,6)-f(i,j,k+1,6)+f(i,j+1,k,6)-f(i,j,k,6))/hy) pp(i,j,k,2)=0.25*((f(i+1,j+1,k+1,6)+f(i+1,j,k+1,6) 1 -f(i+1,j+1,k,6)-f(i+1,j,k,6) 2 +f(i,j+1,k+1,6)+f(i,j,k+1,6)-f(i,j+1,k,6)-f(i,j,k,6))/hz 3 -(f(i+1,j+1,k+1,8)+f(i+1,j,k+1,8)+f(i+1,j+1,k,8)+f(i+1,j,k,8) 4 -f(i,j+1,k+1,8)-f(i,j,k+1,8)-f(i,j+1,k,8)-f(i,j,k,8))/hx) pp(i,j,k,1)=0.25*((f(i+1,j+1,k+1,8)-f(i+1,j,k+1,8) 1 +f(i+1,j+1,k,8)-f(i+1,j,k,8) 2 +f(i,j+1,k+1,8)-f(i,j,k+1,8)+f(i,j+1,k,8)-f(i,j,k,8))/hy 3 -(f(i+1,j+1,k+1,7)+f(i+1,j,k+1,7)-f(i+1,j+1,k,7)-f(i+1,j,k,7) 4 +f(i,j+1,k+1,7)+f(i,j,k+1,7)-f(i,j+1,k,7)-f(i,j,k,7))/hz) 381 continue !hpf$ end on 38 continue return end subroutine equib8(f,pp,nxp,ro01,pr01,rrat,cj,cp) c implicit real*8 (a-h,o-z) parameter (npe=2) !hpf$ processors pe(npe) parameter (nx=320,ny= 80,nz= 80,nb=8) c parameter (nx=480,ny=160,nz=160,nb=8) c parameter (nx=640,ny=160,nz=160,nb=8) c parameter (nx=800,ny=200,nz=200,nb=8) parameter (nx1=nx+1,nx2=nx+2,ny1=ny+1,ny2=ny+2) parameter (nz1=nz+1,nz2=nz+2) c dimension f(nx2,ny2,nz2,nb),pp(nx2,ny2,nz2,3) c dimension gf(nx2,ny2,nz2,nb),gpp(nx2,ny2,nz2,3) c dimension gf(nx2,ny2,nz2,nb) dimension cj(1),cp(1) !hpf$ distribute f(*,*,block,*) onto pe !hpf$ distribute pp(*,*,block,*) onto pe !hpf$ shadow f(0,0,1:1,0) !hpf$ shadow pp(0,0,1:1,0) c !xocl global gf,gpp c equivalence (gf,f),(gpp,pp) c common /blk/f,pp 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 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 n1=nx2 n2=n1*ny2 n3=n2*nz2 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 !hpf$ independent,new(i,j,k, !hpf$* ro1m,ro2m,ro3m,ro1,ro2,ro3,x,y,z,vr1,x1m) do 10 k=1,nz2 !hpf$ on home(f(:,:,k,:)),local(f,i,j,m,nxp,x,y,z, !hpf$* dn1,dv1,ax1,vr1,vsw,x1m,xm,ro1,ro2,ro3, !hpf$* hx,hy,hz,dn,ro1m,ro2m,ro3m,x1,ro02,pr02, !hpf$* cp,cj,p0,b0) begin do 101 j=1,ny2 c vr1=vsw do 101 i=1,nx2 vr1=vsw 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 f(i,j,k,1)=1.0/ro3 c f(i,j,k,1)=ro1**(-x1) if(f(i,j,k,1).lt.ro02) f(i,j,k,1)=ro02 f(i,j,k,5)=1.0e-7*cp(6)/ro2 c f(i,j,k,5)=p0*ro1**(-x2) if(f(i,j,k,5).lt.pr02) f(i,j,k,5)=pr02 f(i,j,k,6)=-2.0*b0*x*z/(ro2*ro2) f(i,j,k,7)=-2.0*b0*y*z/(ro2*ro2) f(i,j,k,8)=b0*(x*x+y*y-z*z)/(ro2*ro2) f(i,j,k,6)=-3.0*b0*x*z/(ro2*ro3) f(i,j,k,7)=-3.0*b0*y*z/(ro2*ro3) f(i,j,k,8)=b0*(x*x+y*y-2.0*z*z)/(ro2*ro3) 101 continue !hpf$ end on 10 continue c c c current !hpf$ reflect f !hpf$ independent,new(i,j,k) do 38 k=1,nz1 !hpf$ on home(pp(:,:,k,:)),local(pp,f,i,j,m,hx,hy,hz) begin do 381 j=1,ny1 do 381 i=1,nx1 pp(i,j,k,3)=0.25*((f(i+1,j+1,k+1,7)+f(i+1,j,k+1,7) 1 +f(i+1,j+1,k,7)+f(i+1,j,k,7) 2 -f(i,j+1,k+1,7)-f(i,j,k+1,7)-f(i,j+1,k,7)-f(i,j,k,7))/hx 3 -(f(i+1,j+1,k+1,6)-f(i+1,j,k+1,6)+f(i+1,j+1,k,6)-f(i+1,j,k,6) 4 +f(i,j+1,k+1,6)-f(i,j,k+1,6)+f(i,j+1,k,6)-f(i,j,k,6))/hy) pp(i,j,k,2)=0.25*((f(i+1,j+1,k+1,6)+f(i+1,j,k+1,6) 1 -f(i+1,j+1,k,6)-f(i+1,j,k,6) 2 +f(i,j+1,k+1,6)+f(i,j,k+1,6)-f(i,j+1,k,6)-f(i,j,k,6))/hz 3 -(f(i+1,j+1,k+1,8)+f(i+1,j,k+1,8)+f(i+1,j+1,k,8)+f(i+1,j,k,8) 4 -f(i,j+1,k+1,8)-f(i,j,k+1,8)-f(i,j+1,k,8)-f(i,j,k,8))/hx) pp(i,j,k,1)=0.25*((f(i+1,j+1,k+1,8)-f(i+1,j,k+1,8) 1 +f(i+1,j+1,k,8)-f(i+1,j,k,8) 2 +f(i,j+1,k+1,8)-f(i,j,k+1,8)+f(i,j+1,k,8)-f(i,j,k,8))/hy 3 -(f(i+1,j+1,k+1,7)+f(i+1,j,k+1,7)-f(i+1,j+1,k,7)-f(i+1,j,k,7) 4 +f(i,j+1,k+1,7)+f(i,j,k+1,7)-f(i,j+1,k,7)-f(i,j,k,7))/hz) 381 continue !hpf$ end on 38 continue c c !hpf$ independent,new(i,j,k, !hpf$* xx,px1,vr1,y1,bx2,by2,bz2,bx1,by1,bz1, !hpf$* ro1m,ro2m,ro3m,ro1,ro2,ro3,x1m,x,y,z) do 40 k=1,nz2 !hpf$ on home(f(:,:,k,:)),local(f,i,j,m,nxp,x,y,z, !hpf$* dn1,dv1,ax1,vr1,vsw,x1m,xm,ro1,ro2,ro3, !hpf$* hx,hy,hz,dn,ro1m,ro2m,ro3m,x1,ro02,pr02, !hpf$* cp,cj,p0,b0,bx1,bx2,by1,by2,bz1,bz2,y1, !hpf$* alp1,bet,px1,ro01,pr01,xx,bis) begin do 401 j=1,ny2 c vr1=vsw do 401 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 f(i,j,k,1)=1.0/ro3 c f(i,j,k,1)=ro1**(-x1) if(f(i,j,k,1).lt.ro02) f(i,j,k,1)=ro02 f(i,j,k,5)=1.0e-7*cp(6)/ro2 c f(i,j,k,5)=p0*ro1**(-x2) if(f(i,j,k,5).lt.pr02) f(i,j,k,5)=pr02 f(i,j,k,6)=-2.0*b0*x*z/(ro2*ro2) f(i,j,k,7)=-2.0*b0*y*z/(ro2*ro2) f(i,j,k,8)=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(i,j,k,6)=bx1+bx2 f(i,j,k,7)=by1+by2 f(i,j,k,8)=bz1+bz2 if(x.lt.xm) f(i,j,k,6)=0.0 if(x.lt.xm) f(i,j,k,7)=0.0 if(x.lt.xm) f(i,j,k,8)=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(i,j,k,6)*f(i,j,k,6)+f(i,j,k,7)*f(i,j,k,7) 1 +f(i,j,k,8)*f(i,j,k,8)) if(abs(vsw).lt.1.0e-8) go to 401 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(i,j,k,2)=vr1 xx=vr1/vsw f(i,j,k,5)=f(i,j,k,5)+(pr01-pr02)*xx f(i,j,k,8)=f(i,j,k,8)+bis*y1 401 continue !hpf$ end on 40 continue c return end subroutine clock(ti) real*8 ti,ti1 c call gettod(ti1) call fjhpf_gettod(ti1) ti=1.0d-6*ti1 return end