integer function opengl_proc();
   opengl_proc=2;end
   integer opengl_proc,window;external opengl_proc
   include <opengl.ins>,nolist
parameter(id=300, jd=192, 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, Dt, cal, sigma, scale
integer mg(0:id,0:jd8)
integer idm,jdm,iter,itr,Kshow

DATA Re/300./,cal/1./,sigma/0.1/,iter/2400/,Kshow/10/,scale/1.e4/
   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./float(jd)
f=0.; om=0.; u=0.; v=0.
call omega(id/4, jd/2,  10, 10.0)
call omega(id*3/4, jd/2, 10, 10.0)
!call omega(id/4, jd*6/16, 10, 10.0)
!call omega(id/4, jd*10/16, 10, -10.0)
call image
do it=1,100;call psiCal;enddo

1 CONTINUE;WRITE(*,*) ' 0-EXIT/1-EXE/2-Re/4-sigma/7-ITER/8-SCALE'
READ (*,*) key;SELECT CASE(key)
CASE(1) ! MAIN PROGRAM XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX!
mg=0
do itr=1,iter
call motion
call visco
call psiCal
if(mod(itr,Kshow)==0)call image
!read(*,*)
enddo
!-------------------------------------------------------------------
CASE(2); WRITE(*,*) Re,' - Reynolds Number'; READ (*,*) Re
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(*,*) scale,' - SCALE'; READ (*,*) scale;call image
CASE(0); GO TO 100;       
END SELECT
GO TO 1
100   CONTINUE;!XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
CONTAINs
subroutine omega(i0, j0, Lom, om0);LLom=Lom*Lom
Do jj=-Lom,Lom; jjjj=jj*jj; do ii=-Lom,Lom; rr=ii*ii+jjjj
if(rr<LLom)om(ii+i0,jj+j0)=om0*exp(-rr/Lom);enddo; EndDo
endsubroutine

subroutine motion
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;EndDo
Vmax=max(maxval(abs(u)),maxval(abs(v)));Dt=cal/Vmax

Do j=1,jdm;om1m=om(1,j)-om(0,j);do i=1,idm;om1p=om(i+1,j)-om(i,j)
s=-u(i,j)*Dt
if(s>0)then;om(i,j)=om1p*s+om(i,j);else;om(i,j)=om1m*s+om(i,j);endif
om1m=om1p;enddo;EndDo

do i=1,idm;om1m=om(i,1)-om(i,0);Do j=1,jdm;om1p=om(i,j+1)-om(i,j)
s=-v(i,j)*Dt
if(s>0)then;om(i,j)=om1p*s+om(i,j);else;om(i,j)=om1m*s+om(i,j);endif
om1m=om1p;enddo;EndDo
endsubroutine

subroutine visco
endsubroutine

subroutine psiCal;real ff !Delta(u)+g=0, sigma=2*h*h/dt
ff(i,j)=(f(i,j)*sm2+f(i,j-1)+f(i-1,j)+f(i+1,j)+f(i,j+1)-om(i,j))*osp2
sm2=sigma-2.;osp2=1./(sigma+2.)
DO I5=0,15; Do j = 1,jdm; do i =mod(j+I5,2)+1,idm,2; f(i,j)=ff(i,j); enddo; 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;
kol=mod(nint(om(i,j)*scale),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)=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,jd+9,GL_rgba,GL_byte,mg);call swap_opengl_buffers()
endsubroutine
END !XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX