integer function opengl_proc(); ! include ,nolist opengl_proc=2;end integer opengl_proc,window;external opengl_proc ! include ,nolist include ,nolist parameter(id=400, jd=100, jd8=jd+8) real f(0:id,0:jd),om(0:id,0:jd),u(0:id,0:jd),v(0:id,0:jd) real Re, h, Reh, om0, Dt, cal, sigma, sc integer mg(0:id,0:jd8),idm,jdm,iter,itr,Kshow Logical IO DATA Re/1000./,om0/0.3/,cal/.5/,sigma/0.3/,iter/1000/,Kshow/10/,sc/5/,IO/.true./ ier=winio@('%sp%ww[no_border]%pv%^og[double]%lw',0,0,id,jd8,opengl_proc,window) jdm=jd-1; idm=id-1 h = 1./jd; Reh=Re*h 10 continue u=1.;v=0.;om=0.;b=0.;e=0.;mg=0;Ud=.3 Do j=0,jd;f(0:id,j)=float(j-jd/2);EndDo om(0:20,jd/2+10:jd/2+20)=om0 om(0:20,jd/2-20:jd/2-10)=-om0 itr=0;call image;call image 1 CONTINUE;WRITE(*,*) ' 0-EXIT/1-EXE/2-Re/3-cal/4-sigma/7-ITER/8-SCALE/81-IO/9-Restart' READ (*,*) key;SELECT CASE(key) CASE(1) ! MAIN PROGRAM XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX! mg=0;DO itr=1,iter call motion_visco call psiPTM if(mod(itr,Kshow)==0)call image !read(*,*) ENDDO !------------------------------------------------------------------- CASE(2); WRITE(*,*) Re,' - Reynolds Number'; READ (*,*) Re; Reh=Re*h CASE(21);WRITE(*,*) om0,' - Enter Vortex'; READ (*,*) om0 om(0:2,jd/2+10:jd/2+20)=om0;om(0:2,jd/2-20:jd/2-10)=-om0 CASE(22);WRITE(*,*) om0,' - Enter Vortex'; READ (*,*) om0;om(0:2,jd/2-5:jd/2+5)=om0 CASE(3); WRITE(*,*) cal,' - calcul parameter'; READ (*,*) cal CASE(4); WRITE(*,*) sigma,' - calcul parameter'; READ (*,*) sigma CASE(7); WRITE(*,*) iter; READ (*,*) iter CASE(77);WRITE(*,*) Kshow; READ (*,*) Kshow CASE(8); WRITE(*,*) sc,' - SCALE'; READ (*,*) sc;call image CASE(81);if(IO)then;IO=.false.;else;IO=.true.;endif;call image;call image CASE(9); GO TO 10 CASE(0); GO TO 100; END SELECT GO TO 1 100 CONTINUE;!XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX CONTAINs subroutine motion_visco !Dt=dt/h 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 Vmax=max(maxval(abs(u)),maxval(abs(v))) Dt=cal/Vmax;sDt=Dt*.5;Ret=Dt/Reh Do j = 1, jdm; om1m=om(1,j)-om(0,j);u1m=u(1,j)-u(0,j); do i = 1, idm om1p=om(i+1,j)-om(i,j);om11=om1p+om1m;om2=om1p-om1m;u1p=u(i+1,j)-u(i,j) dx=-u(i,j)*Dt;if(dx>0.)then;dx=dx/(u1p*sDt+1.);else;dx=dx/(u1m*sDt+1.);endif If(abs(om2)0.)then;om(i,j)=om1p*dx+om(i,j);else;om(i,j)=om1m*dx+om(i,j);endif;EndIf om(i,j)=om2*Ret+om(i,j);om1m=om1p;u1m=u1p; enddo; u(id,j)=max(u(idm,j),Ud);v(id,j)=v(idm,j);EndDo do i = 1, idm; om1m=om(i,1)-om(i,0);v1m=v(i,1)-v(i,0);Do j = 1, jdm om1p=om(i,j+1)-om(i,j);om11=om1p+om1m;om2=om1p-om1m;v1p=v(i,j+1)-v(i,j) dy=-v(i,j)*Dt;if(dy>0.)then;dy=dy/(v1p*sDt+1.);else;dy=dy/(v1m*sDt+1.);endif If(abs(om2)0.)then;om(i,j)=om1p*dy+om(i,j);else;om(i,j)=om1m*dy+om(i,j);endif;EndIf om(i,j)=om2*Ret+om(i,j);om1m=om1p;v1m=v1p; EndDo; enddo Do j=1,jdm;u(id,j)=max(u(idm,j),Ud);v(id,j)=v(idm,j);EndDo endsubroutine subroutine psiPTM;real ff !Delta(f)=om, sigma=2*h*h/dt ff(i,j)=(f(i,j)*sm2+f(i+1,j)+f(i,j+1)+f(i-1,j)+f(i,j-1)-om(i,j))*osp2 sm2=sigma-2.;osp2=1./(sigma+2.) Do j = 1, jdm; do i = 1, idm; f(i,j)=ff(i,j); enddo;f(id,j)=f(idm,j); EndDo Do j=jdm,1,-1; do i=idm,1,-1; f(i,j)=ff(i,j); enddo; EndDo Do j=jdm,1,-1; do i = 1, idm; f(i,j)=ff(i,j); enddo;f(id,j)=f(idm,j); EndDo Do j = 1, jdm; do i=idm,1,-1; f(i,j)=ff(i,j); enddo; EndDo endsubroutine subroutine image integer rg,gb,wh;rg=257;gb=ishft(1,16)+256;wh=(gb+1)*127 Do j = 0, jd; do i = 0, id; if(IO)then;kol=mod(nint(om(i,j)*sc*100),128);else;kol=mod(nint(f(i,j)*sc),128);endif !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)=wh+rg*kol;else;mg(i,j)=wh-gb*kol;endif !White Phone enddo; EndDo i=itr*id/iter;mg(i,jd+1:jd+8)=ishft(127,8) call glDrawPixels(id+1,jd8+1,GL_rgba,GL_byte,mg);call swap_opengl_buffers() endsubroutine END !XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX