integer function opengl_proc();
!   include <opengl.ins>,nolist
   opengl_proc=2;end
   integer opengl_proc,window;external opengl_proc
!   include <clearwin.ins>,nolist
   include <opengl.ins>,nolist
parameter(id=200, jd=400, 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, Dt, cal, sigma, sc, time
integer mg(0:id,0:jd8),idm,jdm,iter,itr,Kshow
Logical IO
DATA Re/1000./,cal/.5/,sigma/0.3/,iter/1000/,Kshow/10/,sc/200/,IO/.false./

   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
f=0.;om=0.;u=0.;v=0.;time=0.
u(1:idm,jd-5:jd)=1.
om=0.003
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'
READ (*,*) key;SELECT CASE(key)
CASE(1) ! MAIN PROGRAM XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX!
mg=0;DO itr=1,iter
call motion_visco
!time=time+h*Dt
!call visco
call psiPTM
if(mod(itr,Kshow)==0)call image
!read(*,*)
     ENDDO
!-------------------------------------------------------------------
CASE(2); WRITE(*,*) Re,' - Reynolds Number'; READ (*,*) Re; Reh=Re*h
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(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)      !j=jd: 1.=f(i,jd+1)-f(i,jdm)
v(i,j)=f(i-1,j)-f(i+1,j)
enddo; EndDo
Vmax=max(maxval(abs(u)),maxval(abs(v)))
Dt=cal/Vmax;sDt=Dt*.5;Ret=Dt/Reh

Do j = 1, jdm; om(0,j)=f(1,j)*2; om(id,j)=f(idm,j)*2; EndDo
do i = 1, idm; om(i,0)=f(i,1)*2; om(i,jd)=f(i,jdm)*2+1.; enddo

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)<abs(om11))Then;om(i,j)=(om2*dx+om11)*dx*.5+om(i,j)
Else;if(dx>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; 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)<abs(om11))Then;om(i,j)=(om2*dy+om11)*dy*.5+om(i,j)
Else;if(dy>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
endsubroutine

subroutine visco;real oom !Jom/Jt=Delta(om)/Re, sigma=2Re*h*h/dt=2Reh/Dt
oom(i,j)=(om(i,j)*sm2+om(i+1,j)+om(i,j+1)+om(i-1,j)+om(i,j-1))*osp2
sm2=(Reh/Dt-1.)*2;osp2=1./(sm2+4.)
Do j = 1, jdm; do i = 1, idm; om(i,j)=oom(i,j); enddo; EndDo
Do j=jdm,1,-1; do i=idm,1,-1; om(i,j)=oom(i,j); enddo; EndDo 
endsubroutine

subroutine psiChess;real uu !Ju/Jt=Delta(u)+g, sigma=h*h/dt
uu(i,j)=(u(i,j)*sm2+u(i+1,j)+u(i,j+1)+u(i-1,j)+u(i,j-1)-om(i,j))*osp2
sm2=sigma-2.;osp2=1./(sigma+2.)
Do j = 1, jdm; do i = 2-mod(j,2), idm, 2; u(i,j)=uu(i,j); enddo; EndDo
Do j = 1, jdm; do i = mod(j,2)+1, idm, 2; u(i,j)=uu(i,j); enddo; 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; 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; 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