integer function opengl_proc(); include ,nolist opengl_proc=2;end integer opengl_proc,window;external opengl_proc !!!parameter(id=800, jd=400, jd8=jd+8) parameter(nr=48, jd=256, id=jd*2, jd8=jd+8) include ,nolist include ,nolist real u(0:id,0:jd),v(0:id,0:jd),f(0:id,0:jd),om(0:id,0:jd),fo(0:jd),alf(0:jd),ema(jd-1) real b(0:id,0:jd),e(0:id,0:jd) real Re, Dt, h, cal, sc,sc2, Vmax,Ud,oReh !Dt = dt/h integer mg(0:id,0:jd+jd8),Ib(-nr:nr),Ld,idm,jdm,io,jo,iter,itr,Kshow DATA Re/1000./,Ld/2/,iter/2500/,cal/0.9/,sc/1000./,Kshow/10/ ier=winio@('%sp%ww[no_border]%pv%^og[double]%lw',0,0,id,jd+jd8,opengl_proc,window) idm=id-1; jdm=jd-1; io=id/3;jo=jd/2 +1; h = 1.0/jd; oReh=1./(Re*h) u=1.;v=0.;om=0.;b=0.;e=0.;mg=0;Ud=.2 Do j=0,jd;f(0:id,j)=float(j-jd/2);fo(j)=float(j-jd/2);EndDo Do jj=0,jd/2;j=jj+jd/2;alf(j)=1.-float(jj)/(jd/2);alf(jd-j)=alf(j);EndDo Do j=1,jdm;ema(j)=(1.-alf(j))/alf(j);EndDo itr=0;sc2=1.e3;call Ib_Jb; call body_bnd; call image 1 CONTINUE;WRITE(*,*) '|0-EXIT|1-EXE|2-Re|3,31-cal,Ld|7-ITER|8-SCALE' READ (*,*) key;SELECT CASE(key) CASE(1) ! MAIN PROGRAM XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX! mg=127+ishft(127,8)+ishft(127,16) do itr=1,iter call TranVis !call conv;call visChess !call psiMgr call psi if(mod(itr,Kshow)==0)then;sc2=sc/maxval(abs(f));call image;endif !read(*,*) enddo !------------------------------------------------------------------- CASE(2); WRITE(*,*) Re; READ (*,*) Re; oReh=1./(Re*h) CASE(3); WRITE(*,*) cal,' - calcul parameter'; READ (*,*) cal CASE(30);WRITE(*,*) Ud,' - min Exit Velocity'; READ (*,*) Ud CASE(31);WRITE(*,*) Ld,' - calcul Levels'; READ (*,*) Ld CASE(32);WRITE(*,*) jo,' - jo'; READ (*,*) jo CASE(7); WRITE(*,*) iter; READ (*,*) iter CASE(77);WRITE(*,*) Kshow; READ (*,*) Kshow CASE(8); WRITE(*,*) sc,' - SC'; READ (*,*) sc;call image CASE(0); GO TO 100; END SELECT GO TO 1 100 CONTINUE;!XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX CONTAINs subroutine Ib_Jb;Ib=0 DO jj=0,nr-1;Ib(jj)=nint(sqrt(float(nr*nr-jj*jj)));if(Ib(jj)==nr)Ib(nr)=jj;Ib(-jj)=Ib(jj);ENDDO;Ib(-nr)=Ib(nr) endsubroutine subroutine body_bnd Do jj=1-nr,nr-1;j=jj+jo;do i=io-Ib(jj),io+Ib(jj);f(i,j)=0.;enddo;EndDo Do j=1,jdm;f(id,j)=(fo(j)*ema(j)+f(idm,j))*alf(j);om(id,j)=om(idm,j)*.5;EndDo !Do j=1,jdm;f(id,j)=(f(idm,j)+fo(j))*.5;om(id,j)=om(idm,j)*.5;EndDo !Do j=1,jdm;f(id,j)=(f(idm,j)*3.+fo(j))*.25;om(id,j)=om(idm,j)*.75;EndDo !Do j=1,jdm;u(id,j)=max(u(idm,j),Ud);f(id,j)=f(idm,j);om(id,j)=om(idm,j)*.5;EndDo endsubroutine subroutine TranVis Do j=1,jdm;do i=1,idm;u(i,j)=f(i,j+1)-f(i,j-1);v(i,j)=f(i-1,j)-f(i+1,j);enddo u(id,j)=max(u(idm,j),Ud);v(id,j)=v(idm,j);EndDo Vm=max(maxval(abs(u)),maxval(abs(v)));Dt=cal/Vm;sDt=Dt/2; Ret=Dt*oReh !Omega on Body's Boundary Do jj=1-nr,nr-1;j=jj+jo; i=io-Ib(jj);om(i,j)=(om(i,j)+f(i,j-1)+f(i-1,j)+f(i+1,j)+f(i,j+1))*.5 i=io+Ib(jj);om(i,j)=(om(i,j)+f(i,j-1)+f(i-1,j)+f(i+1,j)+f(i,j+1))*.5;EndDo Do j=1,jdm;u1m=0.;o1m=0.; do i=1,idm;u1p=u(i+1,j)-u(i,j);o1p=om(i+1,j)-om(i,j);o11=o1p+o1m;o2=o1p-o1m ds=-u(i,j)*Dt;if(ds>0)then;ds=ds/(u1p*sDt+1.);else;ds=ds/(u1m*sDt+1.);endif !Formula: ds=-u*Dt/(u1*sDt+1)! If(abs(o2)0)then;om(i,j)=o1p*ds+om(i,j);else;om(i,j)=o1m*ds+om(i,j);endif;Endif om(i,j)=o2*Ret+om(i,j);u1m=u1p;o1m=o1p;enddo;om(id,j)=om(idm,j);EndDo Do jj=2-nr,nr-2;j=jj+jo;do i=io+1-Ib(jj),io-1+Ib(jj);om(i,j)=0.;enddo;EndDo do i=1,idm;v1m=v(i,1);o1m=om(i,1) Do j=1,jdm;v1p=v(i,j+1)-v(i,j);o1p=om(i,j+1)-om(i,j);o11=o1p+o1m;o2=o1p-o1m ds=-v(i,j)*Dt;if(ds>0)then;ds=ds/(v1p*sDt+1.);else;ds=ds/(v1m*sDt+1.);endif If(abs(o2)0)then;om(i,j)=o1p*ds+om(i,j);else;om(i,j)=o1m*ds+om(i,j);endif;Endif om(i,j)=o2*Ret+om(i,j);v1m=v1p;o1m=o1p;EndDo;enddo ;Do j=1,jdm;om(id,j)=om(idm,j);EndDo Do jj=2-nr,nr-2;j=jj+jo;do i=io+1-Ib(jj),io-1+Ib(jj);om(i,j)=0.;enddo;EndDo endsubroutine subroutine conv !real b(0:id,0:jd),d(0:id,0:jd);b=0.;d=0. Vmax=0.1; Do j=1, jdm; do i=1, idm; u(i,j) = (f(i, j+1) - f(i, j-1))*0.5; v(i,j) = (f(i-1, j) - f(i+1, j))*0.5;enddo;EndDo Vmax=max(maxval(abs(u)),maxval(abs(v)));Dt = cal/Vmax; !Ret=Dt/(Re*h) !Omega on Body's Boundary Do jj=1-nr,nr-1;j=jj+jo; i=io-Ib(jj);om(i,j)=(om(i,j)+f(i,j-1)+f(i-1,j)+f(i+1,j)+f(i,j+1))*.5 i=io+Ib(jj);om(i,j)=(om(i,j)+f(i,j-1)+f(i-1,j)+f(i+1,j)+f(i,j+1))*.5;EndDo Do j=1, jdm; do i = 1,idm dx = -u(i,j)*Dt if(dx>0)then;bb=(om(i+1,j)-om(i,j))*dx;else;bb=(om(i,j)-om(i-1,j))*dx;endif dy = -v(i,j)*Dt if(dy>0)then;b(i,j)=(om(i,j+1)-om(i,j))*dy+bb;else;b(i,j)=(om(i,j)-om(i,j-1))*dy+bb;endif enddo;EndDo om=om+b endsubroutine subroutine visChess;real omega omega(i,j)=(om(i,j)*sm2+om(i+1,j)+om(i,j+1)+om(i-1,j)+om(i,j-1))*osp2 sgm =Re*h/Dt; sm2 = sgm-2.0; osp2 = 1.0/(sgm+2.0) DO k2=0,1;Do j=1,jdm;do i=mod(j+K2,2)+1,idm,2;om(i,j)=omega(i,j);enddo;EndDo ;ENDDO Do jj=2-nr,nr-2;j=jj+jo;do i=io+1-Ib(jj),io-1+Ib(jj);om(i,j)=0.;enddo;EndDo endsubroutine subroutine psi;real ff ff(i,j)=(-f(i,j)*1.5+f(i,j-1)+f(i-1,j)+f(i+1,j)+f(i,j+1)-om(i,j))*0.4 DO k8=0,7;Do j=1,jdm; do i=mod(j+K8,2)+1,idm,2; f(i,j)=ff(i,j); enddo; EndDo call body_bnd; ENDDO endsubroutine subroutine psiMgr !real e(0:id,0:jd);e=0. Do j=1,jdm; do i=1,idm; e(i,j)=-f(i,j)*4+f(i,j-1)+f(i-1,j)+f(i+1,j)+f(i,j+1)-om(i,j); enddo; EndDo DO L=0,Ld;m=ishft(1,L);m2=ishft(m,1) Do j=m,jdm,m; do i=m2,idm,m2; e(i,j)=(e(i,j)*4+e(i,j-m)+e(i-m,j)+e(i+m,j)+e(i,j+m))*.25; enddo; EndDo Do j=m2,jdm,m2; do i=m2,idm,m2; e(i,j)=(e(i,j)*4+e(i-m,j-m)+e(i+m,j-m)+e(i-m,j+m)+e(i+m,j+m))*.25; enddo; EndDo ENDDO Do j=m2,jdm,m2; do i=m2,idm,m2; e(i,j)=e(i,j)*.25; enddo; EndDo DO L=Ld,0,-1;m=ishft(1,L);m2=ishft(m,1) Do j=m,jdm,m2; do i=m2,idm,m2; e(i,j)=(e(i-m,j-m)+e(i+m,j-m)+e(i-m,j+m)+e(i+m,j+m)+e(i,j))*.25; enddo; EndDo Do j=m,jdm,m; do i=m,idm,m2; e(i,j)=(e(i,j-m)+e(i-m,j)+e(i+m,j)+e(i,j+m)+e(i,j))*.25; enddo; EndDo ENDDO !f=f+e;call body_bnd Do j=1,jdm; do i=1,idm; f(i,j)=f(i,j)+e(i,j); enddo; EndDo endsubroutine subroutine image integer rg,gb,red,green,blue rg=257;gb=ishft(1,16)+256 red=127;green=ishft(red,8);blue=ishft(red,16);white=red+green+blue !rg=red+green;gb=green+blue Do j=0,jd;do i=0,id; kol=mod(nint(f(i,j)*sc2),128); !if(kol<0)then;mg(i,j)=wh+rg*kol;else;mg(i,j)=wh-gb*kol;endif !White Phone !if(kol<0)then;mg(i,j)=white+rg*kol;else;mg(i,j)=white-gb*kol;endif !White Phone if(kol>0)then;mg(i,j+jd8)=kol;else;mg(i,j+jd8)=ishft(-kol,16);endif !Black Phone kol=mod(nint(om(i,j)*sc),128) !if(kol>0)then;mg(i,j)=kol;else;mg(i,j)=ishft(-kol,16);endif !Black Phone !if(kol<0)then;mg(i,j+jd8)=wh+rg*kol;else;mg(i,j+jd8)=wh-gb*kol;endif !White Phone if(kol<0)then;mg(i,j)=white+rg*kol;else;mg(i,j)=white-gb*kol;endif !White Phone enddo;EndDo i=itr*id/iter;mg(i,jd+1:jd+7)=0; !ishft(255,8) ! call glDrawPixels(id+1,jd+jd8+1,GL_rgba,GL_byte,mg);call swap_opengl_buffers() endsubroutine END !XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX