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(jd=256, id=jd*2, jd8=jd+8) include ,nolist include ,nolist real T(0:id,0:jd),u(0:id,0:jd),v(0:id,0:jd),f(0:id,0:jd),om(0:id,0:jd) real Gr, Pr, Dt, h, cal, sc,sc2,scNu, Vmax !Dt = dt/h integer mg(0:id,0:jd+jd8),Ld,idm,jdm,iter,itr,kNu DATA Gr/1.e3/,Pr/.7/,Ld/4/,iter/5000/,cal/0.1/,sc/1000./,scNu/500./,kNu/1/ ier=winio@('%sp%ww[no_border]%pv%^og[double]%lw',0,0,id,jd+jd8,opengl_proc,window) u=0.;v=0.;f=0.;om=0.;mg=0 idm=id-1; jdm=jd-1; h = 1.0/jd T=0.; Do j=1,jdm; do i=j+1, idm-j; T(i,j)=0.5; enddo; EndDo; T(0:id,0)=1. itr=0;call image 1 CONTINUE;WRITE(*,*) 'Gr=',Gr,'Pr=',Pr,'|0-EXIT|1-EXE|2,21-Gr,Pr|3,31-cal,Ld|4-Left|5-Right|7-ITER|8-SCALE' READ (*,*) key;SELECT CASE(key) CASE(1) ! MAIN PROGRAM XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX! !mg(0:id,jd:jd+jd8)=127+ishft(127,8)+ishft(127,16) mg=127+ishft(127,8)+ishft(127,16) do itr=1,iter call conv2 call vis !call visChess call psiMgr;call psiMgr! !call psi if(mod(itr,10)==0)then;sc2=sc/maxval(abs(f));call image;endif !read(*,*) enddo !------------------------------------------------------------------- CASE(2); WRITE(*,*) Gr; READ (*,*) Gr CASE(21);WRITE(*,*) Pr; READ (*,*) Pr CASE(3); WRITE(*,*) cal,' - calcul parameter'; READ (*,*) cal CASE(31);WRITE(*,*) Ld,' - calcul Levels'; READ (*,*) Ld CASE(4); WRITE(*,*) 'INPUT k'; READ (*,*) kk;do k=1,jd;T(k,k)=0.;if(k<=kk)T(k,k)=1.;enddo CASE(5); WRITE(*,*) 'INPUT k'; READ (*,*) kk;do k=1,jd;T(id-k,k)=0.;if(k<=kk)T(id-k,k)=1.;enddo CASE(6); open(4, file = 'Tn.dat'); Do j=1,jdm;write(4,'(2f7.3)')T(j+1,j),T(idm-j,j);EndDo;endfile 4;close(4) CASE(61); open(4, file = 'TnDown.dat'); do i=1,idm;write(4,'(f5.3)')T(i,1);enddo;endfile 4;close(4) CASE(7); WRITE(*,*) iter; READ (*,*) iter; CASE(8); WRITE(*,*) sc,' - SC'; READ (*,*) sc;call image CASE(88);WRITE(*,*) scNu,' - scNu'; READ (*,*) scNu;call image CASE(888);WRITE(*,*) kNu,' - kNu'; READ (*,*) kNu;call image CASE(9); call Heat_Flux CASE(0); GO TO 100; END SELECT GO TO 1 100 CONTINUE;!XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX CONTAINs subroutine conv2 Vmax=0.1; Do j=1, jdm; do i=j+1, idm-j; 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; Vmax=max(abs(u(i,j)),abs(v(i,j)),Vmax); enddo;EndDo Dt = cal/Vmax; sDt=Dt*.5; Grt=Gr*Dt*.5 !Dt = dt/h do i = 2, id - 2;om(i,0)=f(i,1)*2;enddo !Omega on Boundary Do j=1, jdm; om(j,j)=(f(j+1,j)+f(j,j-1))*2;om(id-j,j)=(f(idm-j,j)+f(id-j,j-1))*2; EndDo !u - Transit Do j=1,jdm;u1m=u(j+1,j)-u(j,j);T1m=T(j+1,j)-T(j,j);om1m=om(j+1,j)-om(j,j) do i=j+1,idm-j;u1p=u(i+1,j)-u(i,j);T1p=T(i+1,j)-T(i,j);om1p=om(i+1,j)-om(i,j) T11=T1p+T1m;T2=T1p-T1m;om11=om1p+om1m;om2=om1p-om1m dx=-u(i,j)*Dt;if(dx>0)then;dx=dx/(u1p*sDt+1.);else;dx=dx/(u1m*sDt+1.);endif !Formula: ds=-u*Dt/(u1*sDt+1)! If(abs(T2)0)then;T(i,j)=T1p*dx+T(i,j);else;T(i,j)=T1m*dx+T(i,j);endif;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 u1m=u1p;T1m=T1p;om1m=om1p;enddo;EndDo !v - Transit do i=2,id-2;v1m=v(i,1)-v(i,0);T1m=T(i,1)-T(i,0);om1m=om(i,1)-om(i,0) Jmin=min(i,id-i) do j=1,Jmin-1;v1p=v(i,j+1)-v(i,j);T1p=T(i,j+1)-T(i,j);om1p=om(i,j+1)-om(i,j) T11=T1p+T1m;T2=T1p-T1m;om11=om1p+om1m;om2=om1p-om1m dy=-v(i,j)*Dt;if(dy>0)then;dy=dy/(v1p*sDt+1.);else;dy=dy/(v1m*sDt+1.);endif !Formula: ds=-u*Dt/(u1*sDt+1)! If(abs(T2)0)then;T(i,j)=T1p*dy+T(i,j);else;T(i,j)=T1m*dy+T(i,j);endif;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)=(T(i-1,j)-T(i+1,j))*Grt+om(i,j); v1m=v1p;T1m=T1p;om1m=om1p;enddo;EndDo endsubroutine subroutine visChess;real omeg, Temper omeg(i,j)=(om(i,j)*sm2+om(i+1,j)+om(i,j+1)+om(i-1,j)+om(i,j-1))*osp2 Temper(i,j)=(T(i,j)*sTm2+T(i+1,j)+T(i,j+1)+T(i-1,j)+T(i,j-1))*osTp2 sgm =h/Dt; sm2 = sgm-2.0; osp2 = 1.0/(sgm+2.0) sgmT =Pr*h/Dt; sTm2 = sgmT-2.0; osTp2 = 1.0/(sgmT+2.0) DO k2=0,1;Do j=1,jdm; do i=mod(j+K2,2)+j+1,idm-j,2; om(i,j) = omeg(i,j); T(i,j)=Temper(i,j); enddo; EndDo; ENDDO endsubroutine subroutine vis;real omeg, Temper omeg(i,j)=(om(i,j)*sm2+om(i+1,j)+om(i,j+1)+om(i-1,j)+om(i,j-1))*osp2 Temper(i,j)=(T(i,j)*sTm2+T(i+1,j)+T(i,j+1)+T(i-1,j)+T(i,j-1))*osTp2 !sgm =2*h/Dt; sgm =4*h/Dt; sm2 = sgm-2.0; osp2 = 1.0/(sgm+2.0) sgmT =Pr*sgm; sTm2 = sgmT-2.0; osTp2 = 1.0/(sgmT+2.0) Do j=1,jdm; do i=j+1,idm-j; om(i,j) = omeg(i,j); T(i,j)=Temper(i,j); enddo; EndDo Do j=jdm,1,-1; do i=idm-j,j+1,-1; om(i,j) = omeg(i,j); T(i,j)=Temper(i,j); enddo; EndDo Do j=1,jdm; do i=idm-j,j+1,-1;om(i,j) = omeg(i,j); T(i,j)=Temper(i,j); enddo; EndDo Do j=jdm,1,-1;do i=j+1,idm-j; om(i,j) = omeg(i,j); T(i,j)=Temper(i,j); enddo; EndDo endsubroutine subroutine psi;real ff ff(i,j)=(-f(i,j)*1.5+f(i+1,j)+f(i,j+1)+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)+j+1,idm-j,2; f(i,j)=ff(i,j); enddo; EndDo; ENDDO endsubroutine subroutine psiMgr;real e(0:id,0:jd);e=0. Do j=1,jdm; do i=j+1,idm-j; 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=j+m2,idm-j,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=j+m2,idm-j,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=j+m2,idm-j,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=j+m2,idm-j,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=j+m,idm-j,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 Do j=1,jdm; do i=j+1,idm-j; f(i,j)=f(i,j)+e(i,j); enddo; EndDo endsubroutine subroutine Heat_Flux;c=1./sqrt(2.) !44 format('(f7.4)') open(4,file='HF.dat') do i=1,idm;j=i;ja=id-i; write(4,'(f7.4)')T(i,1)-T(i,0); if(i