c #define f(a,b,c,d) F(d,a,b,c) #define ff(a,b,c,d) FF(d,a,b,c) #define p(a,b,c,d) P(d,a,b,c) #define u(a,b,c,d) U(d,a,b,c) c program ogndd7 c c MPI (Massage Passing Interface) c file name /zearthc2/earthc14.f c earthd14.f modified leap-frog scheme c ogtac.cntl(ogncb77) from ognc.cntl(ogncd75) c 3d mhd simulation of dipole tilt in 1/2 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 CC MPI START include 'mpif.h' integer istatus(mpi_status_size) common /para_info/ks,ks1,ke,ke1,kss,irank,isize c for mpi_gatherv parameter (npe=24) integer recvcount(npe),displs(npe) CC MPI END !hpf$ processors pe(npe) c parameter (nx=320,ny= 80,nz=160,nxp=100) c parameter (nx=500,ny=100,nz=200,nxp=150) parameter (nx=500,ny=318,nz=318,nxp=150) c parameter (nx=500,ny=300,nz=300,nxp=150) c parameter (nx=500,ny=200,nz=400,nxp=150) 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)) c parameter (nb=8,nbb=11,n4=n3*nb,n5=n3*nbb,n6=n3*18) parameter (nb=8,nbb=11) c parameter (last=2048,nil=2,lll=18,ibig=0) c parameter (last=8192,nil=2,lll=18,ibig=0) c parameter (last=4096,nil=2,lll=18,ibig=0) c parameter (last=1024,nil=2,lll=18,ibig=0) parameter (last= 8,nil=2,lll=18,ibig=0) c parameter (last=12288,nil=2,lll=18,ibig=0) c parameter (last= 448,nil=2,lll=18,ibig=0) c parameter (last= 32,nil=2,lll=18,ibig=0) c parameter (last= 32,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=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=2048,thx=2.0833333) c parameter (iiq0=8,iip0= 32,iis0=8192,thx=0.5208333) c parameter (iiq0=8,iip0= 32,iis0=4096,thx=3.1266676) c parameter (iiq0=8,iip0= 32,iis0=1024,thx=2.0844451) parameter (iiq0=1,iip0= 1,iis0= 2,thx=2.0844451) c parameter (iiq0=8,iip0= 32,iis0=8192,thx=1.5633338) c parameter (rmuu=0.0020,eud0=0.0020,ro01=5.0e-4,pr01=3.56e-8) 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) 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) 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=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) c parameter (orid= 45.000,ths0=0.0,kk8=k3*8) parameter (orid= 90.000,ths0=0.0,kk8=k3*8) c parameter (orid=315.000,ths0=0.0,kk8=k3*8) c parameter (orid= 0.000,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) c parameter (mfd=6,nfd=1441,iinob=last*780,ttno=512.0) c parameter (mfd=6,nfd=1441,iinob=last*0, ttno=102.4) c parameter (mfd=6,nfd=1441,iinob=last*0, ttno=204.8) c parameter (mfd=6,nfd=1441,iinob=last*25, ttno=819.2) 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) c parameter (mfd=6,nfd=1441,iinob=last*30, ttno=448.0) c parameter (mfd=6,nfd=1441,iinob=last*40,ttno=102.4) c parameter (mfd=6,nfd=1441,iinob=last*12*1,ttno=512.0) double precision zt0,zt1,zt2,zt3,zt integer (kind=8):: n6,n5,n4 c c dimension f(nx2,ny2,nz2,nb),u(nx2,ny2,nz2,nb), c 1 ff(nx2,ny2,nz2,nb),p(nx2,ny2,nz2,nbb), c 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) c dimension gf(nx2,ny2,nz2,nb),gfdd(mfd,nfd) c dimension gf(nx2,ny2),gfdd(mfd,nfd) CC MPI START parameter(nzz=(nz2-1)/npe+1) dimension f(nx2,ny2,0:nzz+1,nb),u(nx2,ny2,0:nzz+1,nb), 1 ff(nx2,ny2,0:nzz+1,nb),p(nx2,ny2,0:nzz+1,nbb), 2 pp(nx2,ny2,0:nzz+1,3) dimension ftemp1(nx2,ny2),ftemp2(nx2,ny2) dimension utemp1(nx2,ny2),utemp2(nx2,ny2) c for all_gather c dimension fg(nx2,ny2,nz2,nb) dimension fg(nx2,ny2,nz2) CC MPI END dimension gfdd(mfd,nfd) 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 data cj/10.0,0.3,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 til vsw ibr ibz data cp/150.3,90.3,90.3,3.5,0.00,2.68,1.35, 30.,.044,0.0,-0.0/ c data cp/150.3,45.3,90.3,3.5,0.00,2.68,1.35,-30.,.044,0.0,-0.0/ c c CC MPI START call clock(zt0) call mpi_init(ier) call mpi_comm_rank(mpi_comm_world,irank,ier) call mpi_comm_size(mpi_comm_world,isize,ier) c kk=nz2/isize kmod=mod(nz2,isize) c ks=1 kss=irank*kk+min(kmod,irank) c if (irank.lt.kmod) kk=kk+1 ke=ks+kk-1 kgs=ks+kss kge=ke+kss ks1=ks ke1=ke if (irank.eq.0) ks1=2 if (irank.eq.isize-1) ke1=ke-1 c c nword=(ke-ks+1)*nx2*ny2 nword=(ke-ks+1)*n2 call mpi_gather(nword,1,mpi_integer,recvcount, * 1,mpi_integer,0,mpi_comm_world,ier) displs(1)=0 do i=2,isize displs(i)=displs(i-1)+recvcount(i-1) end do c c for debug +++ c write(6,*) ' irank , ks , ke ',irank,ks,ke c if (irank.eq.0) write(6,*) ' value of recvcount ', c * (recvcount(ii),ii=1,isize) c if (irank.eq.0) write(6,*) ' value of displs ', c * (displs(ii),ii=1,isize) c CC MPI END c rewind 12 CC MPI START if(irank.eq.0) then c open(10,file='ecc1e030.data',access='sequential', c 1 form='unformatted') c open(11,file='edd1d005.data', c 1 access='sequential',form='unformatted') c open(9,file='wn991021.data',access='sequential', c 1 form='unformatted') c open(9,file='wn991022.data',access='sequential', c 1 form='unformatted') c open(10,file='./data1/earthd/edd1d420.data', c 1 access='sequential',form='unformatted') c c open(12,file='./data1/earthd/edd1d010.data', c 1 access='sequential',form='unformatted') c end if CC MPI END c c c read(9) gfdd do 660 i=1,nfd c read(9,662) (gfdd(j,i),j=1,mfd) do 664 j=1,mfd c fdd(j,i)=gfdd(j,i) 664 continue 662 format(1h ,6(1pe12.4)) 660 continue c c c til0=30.0 til0=0.0 c til0=90.0 tild=5.0 bis0=1.5 c c do 300 iii=1,5 c do 300 iii=1,10 do 300 iii=1,1 c do 300 iii=1,6 c ntap=10+iii c cp(11)=bisz*float(iii-2) cp(11)=bis0 c cp(8)=til0+5.0*float(iii) c cp(8)=til0+10.0*float(iii) cp(8)=til0 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 CC MPI START if(irank.eq.0) 1 write (6,12) iii,last,nx,ny,nz,n1,n2,n3,n4,n5,n6,eat0,rmu0,aru, 2 eud,rrat,hx,hy,hz,t,t1,ro01,pr01,gra,dx2,dy2,dz2,dx4,dy4,dz4, 3 bis,orid,th,bisy,bisz,(cp(i),i=1,11),(cj(j),j=1,10) CC MPI END 12 format(1h ,5x,11i10/(1h ,5x,1p8e15.5)) c if(iii.ge.2) go to 400 c c c do 20 k=1,nz2 CC MPI START do 20 k=ks,ke CC MPI END do 22 m=1,nb !hpf$ independent,new(i,j,k) !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 if(iii.eq.1) call equib7(f,pp,nxp,ro01,pr01,rra1,cj,cp) if(iii.eq.1) call equib8(f,pp,nxp,ro01,pr01,rra1,cj,cp) if(iii.ge.2) call equib8a(u,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 c do 410 ii=1,itap do 174 m=1,nb CC MPI START if(irank.eq.0) then c do 1742 k=1,nz2 c read(10) fg(1:nx2,1:ny2,k) 1742 continue end if c call mpi_scatterv(fg(1,1,1),recvcount,displs, c 1 mpi_real,f(1,1,1,m),nword,mpi_real,0, c 2 mpi_comm_world,ier) CC MPI END 174 continue 410 continue 400 continue c c c !hpf$ independent,new(i,j,k) c do 290 k=2,nz1 CC MPI START do 290 k=ks1,ke1 CC MPI END !hpf$ on home(f(:,:,k,:)),local(f,u,i,j, !hpf$* x,y,z,hx,hy,hz,xx4,ar2,ro2,hhx,x2, !hpf$* aru) begin c xx4=0.5*hx*float(2*nxp-nx1-2) z=0.5*hz*float(2*(kss+k)-nz2-1) do 280 j=2,ny1 c y=0.5*hy*float(2*j-3) y=0.5*hy*float(2*j-ny2-1) c do 2820 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+u(i,j,k,8)*(1.0-x2) f(i,j,k,7)=f(i,j,k,7)*x2+u(i,j,k,7)*(1.0-x2) f(i,j,k,6)=f(i,j,k,6)*x2+u(i,j,k,6)*(1.0-x2) f(i,j,k,5)=f(i,j,k,5)*x2+u(i,j,k,5)*(1.0-x2) f(i,j,k,4)=f(i,j,k,4)*x2+u(i,j,k,4)*(1.0-x2) f(i,j,k,3)=f(i,j,k,3)*x2+u(i,j,k,3)*(1.0-x2) f(i,j,k,2)=f(i,j,k,2)*x2+u(i,j,k,2)*(1.0-x2) f(i,j,k,1)=f(i,j,k,1)*x2+u(i,j,k,1)*(1.0-x2) 2820 continue 280 continue c !hpf$ end on 290 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 c do 100 ii=1,last c tilx=float(ii+last*(iii-1))/float(last) c cp(8)=til0+tild*tilx cp(8)=til0 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 c !hpf$ independent,new(i,j,k,m) c do 31 k=1,nz2 CC MPI START do 31 k=ks,ke CC MPI END !hpf$ on home(f(:,:,k,:)),local(f,i,j,m,x,xx,xx3,xx4, !hpf$* hhx) begin do 30 m=1,nb CC MPI START if(k.eq.1 .and. irank.eq.0) then CC MPI END do 311 j=2,ny1 f(2,j,1,m)=f(1,j,2,m) c f(2,j,k,m)=f(1,j,k+1,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 f(i,j,1,m)=x*f(i-1,j,2,m)+xx*f(i-2,j,3,m) c f(i,j,k,m)=x*f(i-1,j,k+1,m)+xx*f(i-2,j,k+2,m) 311 continue CC MPI START else if(k.eq.ke .and. irank.eq.isize-1) then CC MPI END do 312 j=2,ny1 c f(2,j,nz2,m)=f(1,j,nz1,m) 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 c f(i,j,nz2,m)=x*f(i-1,j,nz1,m)+xx*f(i-2,j,nz,m) 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) c do 32 k=1,nz2 CC MPI START do 32 k=ks,ke CC MPI END !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,1,k,m)=f(1,2,k,m) f(2,ny2,k,m)=f(1,ny1,k,m) c 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,1,k,m)=x*f(i-1,2,k,m)+xx*f(i-2,3,k,m) f(i,ny2,k,m)=x*f(i-1,ny1,k,m)+xx*f(i-2,ny,k,m) c 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 CC MPI START c do m=1,nb c len=nx2*ny2 c if (irank.gt.0) then c call mpi_send(f(1,1,ks,m),n2,mpi_real,irank-1, c & 100,mpi_comm_world,ier) c end if c if (irank.lt.isize-1) then c call mpi_recv(f(1,1,ke+1,m),n2,mpi_real,irank+1, c & 100,mpi_comm_world,istatus,ier) c end if c end do c call mpi_barrier(mpi_comm_world,ier) c c do m=1,nb c if (irank.gt.0) then c call mpi_sendrecv(f(1,1,ks,m),n2,mpi_real,irank-1,100, c & f(1,1,ke+1,m),n2,mpi_real,irank,100, c & mpi_comm_world,istatus,ier) c end if c end do c iright = irank + 1 ileft = irank - 1 if(irank.eq.0) then ileft = MPI_PROC_NULL else if(irank.eq.isize-1) then iright = MPI_PROC_NULL end if do m=1,nb ftemp1=f(:,:,ks,m) ftemp2=f(:,:,ke+1,m) call mpi_sendrecv(ftemp1,n2,mpi_real,ileft,100, & ftemp2,n2,mpi_real,iright,100, & mpi_comm_world,istatus,ier) f(:,:,ks,m)=ftemp1 f(:,:,ke+1,m)=ftemp2 end do c CC MPI END 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) c do 371 k=1,nz1 CC MPI START do 371 k=ks,ke1 CC MPI END !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) c do 90 k=1,nz1 CC MPI START do 90 k=ks,ke1 CC MPI END !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*(kss+k)-nz2) do 50 j=1,ny1 c y=0.5*hy*float(2*j-2) y=0.5*hy*float(2*j-ny2) 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) c do 72 k=1,nz2 CC MPI START do 72 k=ks,ke CC MPI END !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 CC MPI START c do m=1,nb c len=nx2*ny2 c if (irank.lt.isize-1) then c call mpi_send(u(1,1,ke,m),n2,mpi_real,irank+1, c & 110,mpi_comm_world,ier) c call mpi_send(ff(1,1,ke,m),n2,mpi_real,irank+1, c & 120,mpi_comm_world,ier) c end if c if (irank.gt.0) then c call mpi_recv(u(1,1,ks-1,m),n2,mpi_real,irank-1, c & 110,mpi_comm_world,istatus,ier) c call mpi_recv(ff(1,1,ks-1,m),n2,mpi_real,irank-1, c & 120,mpi_comm_world,istatus,ier) c end if c if (irank.gt.0) then c call mpi_send(ff(1,1,ks,m),n2,mpi_real,irank-1, c & 125,mpi_comm_world,ier) c end if c if (irank.lt.isize-1) then c call mpi_recv(ff(1,1,ke+1,m),n2,mpi_real,irank+1, c & 125,mpi_comm_world,istatus,ier) c end if c end do c c do m=1,nb c if (irank.lt.isize-1) then c call mpi_sendrecv(u(1,1,ke,m),n2,mpi_real,irank+1,110, c & u(1,1,ks-1,m),n2,mpi_real,irank,110, c & mpi_comm_world,istatus,ier) c end if c c if (irank.lt.isize-1) then c call mpi_sendrecv(ff(1,1,ke,m),n2,mpi_real,irank+1,120, c & ff(1,1,ks-1,m),n2,mpi_real,irank,120, c & mpi_comm_world,istatus,ier) c end if c if (irank.gt.0) then c call mpi_sendrecv(ff(1,1,ks,m),n2,mpi_real,irank-1,125, c & ff(1,1,ke+1,m),n2,mpi_real,irank,125, c & mpi_comm_world,istatus,ier) c end if c end do c c do m=1,3 c len=nx2*ny2 c if (irank.lt.isize-1) then c call mpi_send(pp(1,1,ke,m),n2,mpi_real,irank+1, c & 130,mpi_comm_world,ier) c end if c if (irank.gt.0) then c call mpi_recv(pp(1,1,ks-1,m),n2,mpi_real,irank-1, c & 130,mpi_comm_world,istatus,ier) c end if c end do c call mpi_barrier(mpi_comm_world,ier) c c do m=1,3 c if (irank.lt.isize-1) then c call mpi_sendrecv(pp(1,1,ke,m),n2,mpi_real,irank+1,130, c & pp(1,1,ks-1,m),n2,mpi_real,irank,130, c & mpi_comm_world,istatus,ier) c end if c end do c iright = irank + 1 ileft = irank - 1 if(irank.eq.0) then ileft = MPI_PROC_NULL else if(irank.eq.isize-1) then iright = MPI_PROC_NULL end if c do m=1,nb utemp1=u(:,:,ke,m) utemp2=u(:,:,ks-1,m) call mpi_sendrecv(utemp1,n2,mpi_real,iright,110, & utemp2,n2,mpi_real,ileft,110, & mpi_comm_world,istatus,ier) u(:,:,ke,m)=utemp1 u(:,:,ks-1,m)=utemp2 call mpi_sendrecv(ff(1,1,ke,m),n2,mpi_real,iright,120, & ff(1,1,ks-1,m),n2,mpi_real,ileft,120, & mpi_comm_world,istatus,ier) call mpi_sendrecv(ff(1,1,ks,m),n2,mpi_real,ileft,125, & ff(1,1,ke+1,m),n2,mpi_real,iright,125, & mpi_comm_world,istatus,ier) end do do m=1,3 call mpi_sendrecv(pp(1,1,ke,m),n2,mpi_real,iright,130, & pp(1,1,ks-1,m),n2,mpi_real,ileft,130, & mpi_comm_world,istatus,ier) end do c CC MPI END c c c second step c c !hpf$ independent,new(i,j,k,m),reduction(max:vmax) c do 190 k=2,nz1 CC MPI START do 190 k=ks1,ke1 CC MPI END !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*(kss+k)-nz2-1) do 80 j=2,ny1 c y=0.5*hy*float(2*j-3) y=0.5*hy*float(2*j-ny2-1) 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) CC MPI START call mpi_allreduce(vmax,vmax1,1,mpi_real,mpi_max, * mpi_comm_world,ier) vmax=vmax1 CC MPI END 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 c th=float(ii)+ths0*float(last) c th=orid*th*pi/180.0 c th=orid*pi/180.0 c ith=ii+last*(iii-1) c th1=180.0*float(ith)/float(last*4) c th2=orid+th1 c if(th2.gt.360.0) th2=th2-360.0 c th=th2*pi/180.0 c bisy=bis*cos(th) c bisz=bis*sin(th) ith=ii+last*(iii-1) th1=180.0*float(ith)/float(last*4) th2=orid+th1 c c iino=ii+last*(iii-1)+iinob c jjno=ifix(float(iino)/ttno)+1 x2=2.0 x1=2.0*x2 nx11=nx1-2*nxp c roo1=fdd(2,jjno) c vswx=-fdd(3,jjno) c pro1=fdd(4,jjno) c bisy=-fdd(5,jjno) c bisz=fdd(6,jjno) roo1=ro01 vswx=vsw pro1=pr01 c bisy=0.0 c bisz=bis c !hpf$ independent,new(i,j,k) c do 102 k=1,nz2 CC MPI START do 102 k=ks,ke CC MPI END !hpf$ on home(f(:,:,k,:)),local(f,i,j,m,nx11,x,x1,hx,hhx, !hpf$* bis,bisz,bisy,pro1,vswx,roo1) 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 c f(i,j,k,8)=(1.0-x)*f(i,j,k,8)+bis*x f(i,j,k,8)=(1.0-x)*f(i,j,k,8)+bisz*x f(i,j,k,7)=(1.0-x)*f(i,j,k,7)+bisy*x 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)+pro1*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)+vswx*x f(i,j,k,1)=(1.0-x)*f(i,j,k,1)+roo1*x 1021 continue !hpf$ end on 102 continue c !hpf$ independent,new(i,j,k) c do 104 k=1,nz1 CC MPI START do 104 k=ks,ke1 CC MPI END !hpf$ on home(u(:,:,k,:)),local(u,i,j,m,nx11,x,x1,hx,hhx, !hpf$* bis,bisz,bisy,pro1,vswx,roo1) 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 c u(i,j,k,8)=(1.0-x)*u(i,j,k,8)+bis*x u(i,j,k,8)=(1.0-x)*u(i,j,k,8)+bisz*x u(i,j,k,7)=(1.0-x)*u(i,j,k,7)+bisy*x 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)+pro1*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)+vswx*x u(i,j,k,1)=(1.0-x)*u(i,j,k,1)+roo1*x 1041 continue !hpf$ end on 104 continue c c c inner boundary condition call equib8a(ff,pp,nxp,ro01,pr01,rra1,cj,cp) c !hpf$ independent,new(i,j,k) c do 390 k=2,nz1 CC MPI START do 390 k=ks1,ke1 CC MPI END !hpf$ on home(f(:,:,k,:)),local(f,ff,i,j, !hpf$* x,y,z,hx,hy,hz,xx4,ar2,ro2,hhx,x2, !hpf$* aru) begin c xx4=0.5*hx*float(2*nxp-nx1-2) z=0.5*hz*float(2*(kss+k)-nz2-1) do 380 j=2,ny1 c y=0.5*hy*float(2*j-3) y=0.5*hy*float(2*j-ny2-1) c do 3820 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) 3820 continue 380 continue c !hpf$ end on 390 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 write(ntap) f do 173 m=1,nb CC MPI START call mpi_barrier(mpi_comm_world,ier) c call mpi_gatherv(f(1,1,1,m),nword,mpi_real,fg(1,1,1), c 1 recvcount,displs,mpi_real,0, c 2 mpi_comm_world,ier) c call mpi_barrier(mpi_comm_world,ier) if(irank.eq.0) then do 1732 k=1,nz2 c write(ntap) fg(1:nx2,1:ny2,k) 1732 continue end if CC MPI END 173 continue c c c CC MPI START if(irank.eq.0) * write(6,661) ii,last,iii,ith,nin,ntap,itap,th2,th,x2,vmax CC MPI END 661 format(1h ,1x,7i10/1x,1p4e15.7) c call clock(zt3) zt1=zt1-zt0 zt2=zt2-zt0 zt3=zt3-zt0 zt=zt2-zt1 CC MPI START if(irank.eq.0) * write(6,404) ii,zt0,zt1,zt2,zt3,zt CC MPI END 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 call mpi_finalize(ierr) stop end subroutine equib7(f,pp,nxp,ro01,pr01,rrat,cj,cp) c implicit real*8 (a-h,o-z) parameter (npe=24) !hpf$ processors pe(npe) c parameter (nx=320,ny= 80,nz=160,nb=8) c parameter (nx=500,ny=100,nz=200,nb=8) parameter (nx=500,ny=318,nz=318,nb=8) c parameter (nx=500,ny=300,nz=300,nb=8) c parameter (nx=500,ny=200,nz=400,nb=8) parameter (nx1=nx+1,nx2=nx+2,ny1=ny+1,ny2=ny+2) parameter (nz1=nz+1,nz2=nz+2) c 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) CC MPI START parameter(nzz=(nz2-1)/npe+1) dimension f(nx2,ny2,0:nzz+1,nb),pp(nx2,ny2,0:nzz+1,3) dimension ftemp1(nx2,ny2),ftemp2(nx2,ny2) dimension utemp1(nx2,ny2),utemp2(nx2,ny2) CC MPI END 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 CC MPI START common /para_info/ks,ks1,ke,ke1,kss,irank,isize CC MPI END 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) tl=pi*cp(8)/180.0 cotl=cos(tl) sitl=sin(tl) c !hpf$ independent,new(i,j,k, !hpf$* xx,y1,px1,bz,by,bx,ro3,ro2,ro1,zt,xt,x,vr1,y,z) c do 10 k=1,nz2 CC MPI START do 10 k=ks,ke CC MPI END !hpf$ on home(f(:,:,k,:)),local begin z=0.5*hz*float(2*(kss+k)-nz2-1) do 101 j=1,ny2 c vr1=vsw c y=0.5*hy*float(2*j-3) y=0.5*hy*float(2*j-ny2-1) do 101 i=1,nx2 vr1=vsw x=0.5*hx*float(2*i-nx2-1+2*nxp) c y=0.5*hy*float(2*j-3) c z=0.5*hz*float(2*k-nz2-1) xt=x*cotl+z*sitl zt=-x*sitl+z*cotl 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 c bx=-2.0*b0*xt*zt/(ro2*ro2) by=-2.0*b0*y*zt/(ro2*ro2) bz=b0*(xt*xt+y*y-zt*zt)/(ro2*ro2) bx=-3.0*b0*xt*zt/(ro2*ro3) by=-3.0*b0*y*zt/(ro2*ro3) bz=b0*(xt*xt+y*y-2.0*zt*zt)/(ro2*ro3) f(i,j,k,6)=bx*cotl-bz*sitl f(i,j,k,7)=by f(i,j,k,8)=bx*sitl+bz*cotl 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,1)=f(i,j,k,1)+(ro01-ro02)*xx 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 c CC MPI START c do m=1,nb c len=nx2*ny2 c if (irank.gt.0) then c call mpi_send(f(1,1,ks,m),n2,mpi_real,irank-1, c & 140,mpi_comm_world,ier) c end if c if (irank.lt.isize-1) then c call mpi_recv(f(1,1,ke+1,m),n2,mpi_real,irank+1, c & 140,mpi_comm_world,istatus,ier) c end if c end do c call mpi_barrier(mpi_comm_world,ier) c c do m=1,nb c if (irank.gt.0) then c call mpi_sendrecv(f(1,1,ks,m),n2,mpi_real,irank-1,140, c & f(1,1,ke+1,m),n2,mpi_real,irank,140, c & mpi_comm_world,istatus,ier) c end if c end do c iright = irank + 1 ileft = irank - 1 if(irank.eq.0) then ileft = MPI_PROC_NULL else if(irank.eq.isize-1) then iright = MPI_PROC_NULL end if do m=1,nb ftemp1=f(:,:,ks,m) ftemp2=f(:,:,ke+1,m) call mpi_sendrecv(ftemp1,n2,mpi_real,ileft,140, & ftemp2,n2,mpi_real,iright,140, & mpi_comm_world,istatus,ier) f(:,:,ks,m)=ftemp1 f(:,:,ke+1,m)=ftemp2 end do c CC MPI END c !hpf$ independent,new(i,j,k) c do 38 k=1,nz1 CC MPI START do 38 k=ks,ke1 CC MPI END !hpf$ on home(pp(:,:,k,:)),local 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 return end subroutine equib8(f,pp,nxp,ro01,pr01,rrat,cj,cp) c implicit real*8 (a-h,o-z) parameter (npe=24) !hpf$ processors pe(npe) c parameter (nx=320,ny= 80,nz=160,nb=8) c parameter (nx=500,ny=100,nz=200,nb=8) parameter (nx=500,ny=318,nz=318,nb=8) c parameter (nx=500,ny=300,nz=300,nb=8) c parameter (nx=500,ny=200,nz=400,nb=8) parameter (nx1=nx+1,nx2=nx+2,ny1=ny+1,ny2=ny+2) parameter (nz1=nz+1,nz2=nz+2) c 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) CC MPI START parameter(nzz=(nz2-1)/npe+1) dimension f(nx2,ny2,0:nzz+1,nb),pp(nx2,ny2,0:nzz+1,3) dimension ftemp1(nx2,ny2),ftemp2(nx2,ny2) dimension utemp1(nx2,ny2),utemp2(nx2,ny2) CC MPI END 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 CC MPI START include 'mpif.h' integer istatus(mpi_status_size) common /para_info/ks,ks1,ke,ke1,kss,irank,isize CC MPI END 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) tl=pi*cp(8)/180.0 cotl=cos(tl) sitl=sin(tl) 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$* bx,by,bz,ro1,ro2,ro3,xt,zt,x,y,z,vr1) c do 10 k=1,nz2 CC MPI START do 10 k=ks,ke CC MPI END !hpf$ on home(f(:,:,k,:)),local begin z=0.5*hz*float(2*(kss+k)-nz2-1) do 101 j=1,ny2 c vr1=vsw c y=0.5*hy*float(2*j-3) y=0.5*hy*float(2*j-ny2-1) do 101 i=1,nx2 vr1=vsw x=0.5*hx*float(2*i-nx2-1+2*nxp) c y=0.5*hy*float(2*j-3) c z=0.5*hz*float(2*k-nz2-1) xt=x*cotl+z*sitl zt=-x*sitl+z*cotl 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 c bx=-2.0*b0*xt*zt/(ro2*ro2) by=-2.0*b0*y*zt/(ro2*ro2) bz=b0*(xt*xt+y*y-zt*zt)/(ro2*ro2) bx=-3.0*b0*xt*zt/(ro2*ro3) by=-3.0*b0*y*zt/(ro2*ro3) bz=b0*(xt*xt+y*y-2.0*zt*zt)/(ro2*ro3) f(i,j,k,6)=bx*cotl-bz*sitl f(i,j,k,7)=by f(i,j,k,8)=bx*sitl+bz*cotl 101 continue !hpf$ end on 10 continue c c write(6,372) nx2,ny2,nz2,nxp, c 1 b0,cotl,sitl,ro02,pr02,rrat,ro01,pr01 c 372 format(1h ,4i4/8(1pe10.3)) c c current !hpf$ reflect f CC MPI START c do m=1,nb c len=nx2*ny2 c if (irank.gt.0) then c call mpi_send(f(1,1,ks,m),n2,mpi_real,irank-1, c & 140,mpi_comm_world,ier) c end if c if (irank.lt.isize-1) then c call mpi_recv(f(1,1,ke+1,m),n2,mpi_real,irank+1, c & 140,mpi_comm_world,istatus,ier) c end if c end do c call mpi_barrier(mpi_comm_world,ier) c c do m=1,nb c if (irank.gt.0) then c call mpi_sendrecv(f(1,1,ks,m),n2,mpi_real,irank-1,140, c & f(1,1,ke+1,m),n2,mpi_real,irank,140, c & mpi_comm_world,istatus,ier) c end if c end do c iright = irank + 1 ileft = irank - 1 if(irank.eq.0) then ileft = MPI_PROC_NULL else if(irank.eq.isize-1) then iright = MPI_PROC_NULL end if do m=1,nb ftemp1=f(:,:,ks,m) ftemp2=f(:,:,ke+1,m) call mpi_sendrecv(ftemp1,n2,mpi_real,ileft,140, & ftemp2,n2,mpi_real,iright,140, & mpi_comm_world,istatus,ier) f(:,:,ks,m)=ftemp1 f(:,:,ke+1,m)=ftemp2 end do c CC MPI END c !hpf$ independent,new(i,j,k) c do 38 k=1,nz1 CC MPI START do 38 k=ks,ke1 CC MPI END !hpf$ on home(pp(:,:,k,:)),local 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,bx22,by22,bz22,bx11,by11,bz11, !hpf$* bx2,by2,bz2,bx1,by1,bz1,ro1,ro2,ro3, !hpf$* ro1m,ro2m,ro3m,xtm,ztm,xt,zt,x1m,x,y,z) c do 40 k=1,nz2 CC MPI START do 40 k=ks,ke CC MPI END !hpf$ on home(f(:,:,k,:)),local begin z=0.5*hz*float(2*(kss+k)-nz2-1) do 401 j=1,ny2 c vr1=vsw c y=0.5*hy*float(2*j-3) y=0.5*hy*float(2*j-ny2-1) do 401 i=1,nx2 x=0.5*hx*float(2*i-nx2-1+2*nxp) c y=0.5*hy*float(2*j-3) c z=0.5*hz*float(2*k-nz2-1) x1m=x-2.0*xm c x1m=x+2.0*xm xt=x*cotl+z*sitl zt=-x*sitl+z*cotl xtm=x1m*cotl-z*sitl ztm=x1m*sitl+z*cotl 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 c bx1=-2.0*b0*xt*zt/(ro2*ro2) by1=-2.0*b0*y*zt/(ro2*ro2) bz1=b0*(xt*xt+y*y-zt*zt)/(ro2*ro2) bx1=-3.0*b0*xt*zt/(ro2*ro3) by1=-3.0*b0*y*zt/(ro2*ro3) bz1=b0*(xt*xt+y*y-2.0*zt*zt)/(ro2*ro3) bx2=-2.0*b0*xtm*ztm/(ro2m*ro2m) by2=-2.0*b0*y*ztm/(ro2m*ro2m) bz2=b0*(xtm*xtm+y*y-ztm*ztm)/(ro2m*ro2m) bx2=-3.0*b0*xtm*ztm/(ro2m*ro3m) by2=-3.0*b0*y*ztm/(ro2m*ro3m) bz2=b0*(xtm*xtm+y*y-2.0*ztm*ztm)/(ro2m*ro3m) c bx11=bx1*cotl-bz1*sitl by11=by1 bz11=bx1*sitl+bz1*cotl bx22=bx2*cotl+bz2*sitl by22=by2 bz22=-bx2*sitl+bz2*cotl c bx1=-3.0*b0*x*z/(ro2*ro3) c by1=-3.0*b0*y*z/(ro2*ro3) c bz1=b0*(x*x+y*y-2.0*z*z)/(ro2*ro3) c bx2=-3.0*b0*x1m*z/(ro2m*ro3m) c by2=-3.0*b0*y*z/(ro2m*ro3m) c bz2=b0*(x1m*x1m+y*y-2.0*z*z)/(ro2m*ro3m) f(i,j,k,6)=bx11+bx22 f(i,j,k,7)=by11+by22 f(i,j,k,8)=bz11+bz22 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 equib8a(u,pp,nxp,ro01,pr01,rrat,cj,cp) c implicit real*8 (a-h,o-z) parameter (npe=24) !hpf$ processors pe(npe) c parameter (nx=320,ny= 80,nz=160,nb=8) c parameter (nx=500,ny=100,nz=200,nb=8) parameter (nx=500,ny=318,nz=318,nb=8) c parameter (nx=500,ny=300,nz=300,nb=8) c parameter (nx=500,ny=200,nz=400,nb=8) parameter (nx1=nx+1,nx2=nx+2,ny1=ny+1,ny2=ny+2) parameter (nz1=nz+1,nz2=nz+2) c c dimension u(nx2,ny2,nz2,nb),pp(nx2,ny2,nz2,3) CC MPI START parameter(nzz=(nz2-1)/npe+1) dimension u(nx2,ny2,0:nzz+1,nb),pp(nx2,ny2,0:nzz+1,3) dimension ftemp1(nx2,ny2),ftemp2(nx2,ny2) dimension utemp1(nx2,ny2),utemp2(nx2,ny2) CC MPI END dimension cj(1),cp(1) !hpf$ distribute u(*,*,block,*) onto pe !hpf$ distribute pp(*,*,block,*) onto pe !hpf$ shadow u(0,0,1:1,0) !hpf$ shadow pp(0,0,1:1,0) 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 CC MPI START include 'mpif.h' integer istatus(mpi_status_size) common /para_info/ks,ks1,ke,ke1,kss,irank,isize CC MPI END 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) tl=pi*cp(8)/180.0 cotl=cos(tl) sitl=sin(tl) 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$* bx,by,bz,ro1,ro2,ro3,xt,zt,x,y,z,vr1) c do 10 k=1,nz2 CC MPI START do 10 k=ks,ke CC MPI END !hpf$ on home(u(:,:,k,:)),local begin z=0.5*hz*float(2*(kss+k)-nz2-1) do 101 j=1,ny2 c vr1=vsw c y=0.5*hy*float(2*j-3) y=0.5*hy*float(2*j-ny2-1) do 101 i=1,nx2 vr1=vsw x=0.5*hx*float(2*i-nx2-1+2*nxp) c y=0.5*hy*float(2*j-3) c z=0.5*hz*float(2*k-nz2-1) xt=x*cotl+z*sitl zt=-x*sitl+z*cotl ro2=x*x+y*y+z*z ro1=sqrt(ro2) ro3=ro1*ro2 u(i,j,k,1)=1.0/ro3 c u(i,j,k,1)=ro1**(-x1) if(u(i,j,k,1).lt.ro02) u(i,j,k,1)=ro02 u(i,j,k,5)=1.0e-7*cp(6)/ro2 c u(i,j,k,5)=p0*ro1**(-x2) if(u(i,j,k,5).lt.pr02) u(i,j,k,5)=pr02 c bx=-2.0*b0*xt*zt/(ro2*ro2) by=-2.0*b0*y*zt/(ro2*ro2) bz=b0*(xt*xt+y*y-zt*zt)/(ro2*ro2) bx=-3.0*b0*xt*zt/(ro2*ro3) by=-3.0*b0*y*zt/(ro2*ro3) bz=b0*(xt*xt+y*y-2.0*zt*zt)/(ro2*ro3) u(i,j,k,6)=bx*cotl-bz*sitl u(i,j,k,7)=by u(i,j,k,8)=bx*sitl+bz*cotl 101 continue !hpf$ end on 10 continue c c write(6,372) nx2,ny2,nz2,nxp, c 1 b0,cotl,sitl,ro02,pr02,rrat,ro01,pr01 c 372 format(1h ,4i4/8(1pe10.3)) c c current !hpf$ reflect u CC MPI START c do m=1,nb c len=nx2*ny2 c if (irank.gt.0) then c call mpi_send(u(1,1,ks,m),n2,mpi_real,irank-1, c & 140,mpi_comm_world,ier) c end if c if (irank.lt.isize-1) then c call mpi_recv(u(1,1,ke+1,m),n2,mpi_real,irank+1, c & 140,mpi_comm_world,istatus,ier) c end if c end do c call mpi_barrier(mpi_comm_world,ier) c c do m=1,nb c if (irank.gt.0) then c call mpi_sendrecv(u(1,1,ks,m),n2,mpi_real,irank-1,150, c & u(1,1,ke+1,m),n2,mpi_real,irank,150, c & mpi_comm_world,istatus,ier) c end if c end do c iright = irank + 1 ileft = irank - 1 if(irank.eq.0) then ileft = MPI_PROC_NULL else if(irank.eq.isize-1) then iright = MPI_PROC_NULL end if do m=1,nb utemp1=u(:,:,ks,m) utemp2=u(:,:,ke+1,m) call mpi_sendrecv(utemp1,n2,mpi_real,ileft,150, & utemp2,n2,mpi_real,iright,150, & mpi_comm_world,istatus,ier) u(:,:,ks,m)=utemp1 u(:,:,ke+1,m)=utemp2 end do c CC MPI END c !hpf$ independent,new(i,j,k) c do 38 k=1,nz1 CC MPI START do 38 k=ks,ke1 CC MPI END !hpf$ on home(pp(:,:,k,:)),local begin do 381 j=1,ny1 do 381 i=1,nx1 pp(i,j,k,3)=0.25*((u(i+1,j+1,k+1,7)+u(i+1,j,k+1,7) 1 +u(i+1,j+1,k,7)+u(i+1,j,k,7) 2 -u(i,j+1,k+1,7)-u(i,j,k+1,7)-u(i,j+1,k,7)-u(i,j,k,7))/hx 3 -(u(i+1,j+1,k+1,6)-u(i+1,j,k+1,6)+u(i+1,j+1,k,6)-u(i+1,j,k,6) 4 +u(i,j+1,k+1,6)-u(i,j,k+1,6)+u(i,j+1,k,6)-u(i,j,k,6))/hy) pp(i,j,k,2)=0.25*((u(i+1,j+1,k+1,6)+u(i+1,j,k+1,6) 1 -u(i+1,j+1,k,6)-u(i+1,j,k,6) 2 +u(i,j+1,k+1,6)+u(i,j,k+1,6)-u(i,j+1,k,6)-u(i,j,k,6))/hz 3 -(u(i+1,j+1,k+1,8)+u(i+1,j,k+1,8)+u(i+1,j+1,k,8)+u(i+1,j,k,8) 4 -u(i,j+1,k+1,8)-u(i,j,k+1,8)-u(i,j+1,k,8)-u(i,j,k,8))/hx) pp(i,j,k,1)=0.25*((u(i+1,j+1,k+1,8)-u(i+1,j,k+1,8) 1 +u(i+1,j+1,k,8)-u(i+1,j,k,8) 2 +u(i,j+1,k+1,8)-u(i,j,k+1,8)+u(i,j+1,k,8)-u(i,j,k,8))/hy 3 -(u(i+1,j+1,k+1,7)+u(i+1,j,k+1,7)-u(i+1,j+1,k,7)-u(i+1,j,k,7) 4 +u(i,j+1,k+1,7)+u(i,j,k+1,7)-u(i,j+1,k,7)-u(i,j,k,7))/hz) 381 continue !hpf$ end on 38 continue c return end subroutine clock(ti) double precision mpi_wtime real*8 ti,ti1 c call gettod(ti1) c call fjhpf_gettod(ti1) c call xclock(ti1,5) ti1=mpi_wtime() c ti=1.0d-6*ti1 ti=ti1 return end