Fortran语言--自由界面程序。

还是老规矩先宣传一下QQ群： 格子玻尔兹曼救星：293267908。

  PROGRAM dambreak
    USE IFPORT 
 
    implicit none
!	//LBM model
	Real*8,  parameter:: w(0:8) = (/4.0d0/9.0d0,1.0d0/9.0d0,1.0d0/9.0d0,1.0d0/9.0d0,1.0d0/9.0d0,1.0d0/36.0d0,1.0d0/36.0d0,1.0d0/36.0d0,1.0d0/36.0d0/)
    Real*8,  parameter:: rhoA = 1.d0
    Integer, parameter:: e(0:8,0:1)=(/(/0,1,0,-1, 0,1,-1,-1, 1/),(/0,0,1, 0,-1,1, 1,-1,-1/)/)
	Integer, parameter:: inv(0:8) = (/0,3,4,1,2,7,8,5,6/)
	Integer, parameter:: xDim=5900,yDim=920,xFluid=3400,yFluid=600
 
	
    integer, parameter:: itslip(0:8) = (/0,1,2,3,2,5,6,6,5/)
    integer, parameter:: ibslip(0:8) = (/0,1,4,3,4,8,7,7,8/)
    integer, parameter:: irslip(0:8) = (/0,1,2,1,4,5,5,8,8/)
    integer, parameter:: ilslip(0:8) = (/0,3,2,3,4,6,6,7,7/)
	
!   Real*8 :: w(0:8),e(0:8,0:1),inv(0:8)  !Real*8为双精度实型
    Real*8 :: rho(1:xDim,1:yDim),mass(1:xDim,1:yDim),u(1:xDim,1:yDim,0:1)
	Real*8 :: fEq(1:xDim,1:yDim,0:8),f(1:xDim,1:yDim,0:8),fpost(1:xDim,1:yDim,0:8)    
	Real*8 :: uxy,uSqr,F_eq,mass1,mass2,P_lt2ph
    Real*8 :: afb,aft,afr,afl,adb,adt,adr,adl
	Real*8 :: dx,dt,vis_phy,gra_phy,den_phy   
	Real*8 :: gravn,vis_lb,tau,omega,times,p_lb
	
    Integer:: iflag(0:xDim+1,0:yDim+1)
    Integer:: i,tStep,x,y,tMax,BCTYPE 
    Character ( len = 100 ) :: command_file_name = 'file_commands.txt'
    integer::csys
    
!   \\Deleting the former folder
    csys = SYSTEM('C:\Windows\System32\cmd.exe /E:ON /V:ON /K "phan.bat" ')
!    
    
	CALL Readparameters(tMax,BCTYPE,dx,dt,vis_phy,gra_phy,den_phy)
	print*,tMax,BCTYPE,dx,dt,vis_phy,gra_phy,den_phy
	
 
	if (BCTYPE.eq.1) then    !Non-slip BC
		adr = 1.0d0; afr=1.d0 - adr
		adb = 1.0d0; afb=1.d0 - adb
		adt = 1.0d0; aft=1.d0 - adt
		adl = 1.0d0; afl=1.d0 - adl
	elseif(BCTYPE.eq.2) then !Slip BC
		adr = 0.0d0; afr=1.d0 - adr
		adb = 0.0d0; afb=1.d0 - adb
		adt = 0.0d0; aft=1.d0 - adt
		adl = 0.0d0; afl=1.d0 - adl
	elseif(BCTYPE.eq.3) then !Hybrid BC
		adr = 0.0d0; afr=1.d0 - adr
		adb = 0.0d0; afb=1.d0 - adb
		adt = 1.0d0; aft=1.d0 - adt
		adl = 0.0d0; afl=1.d0 - adl
	endif
!   //Unit conversion
    P_lt2ph=dx*dx/(dt*dt)
!    print*,P_lt2ph
!    pause
!    stop
 
!	//Calculating relaxation time
	gravn=gra_phy*dt*dt/dx !重力的无量纲值
    vis_lb=vis_phy*dt/(dx*dx) !viscosity的无量纲值
    tau=3.d0*vis_lb+0.5d0
    omega=1.995d0
!    omega=1.d0/tau
    print*,omega,tau
!
    CALL iniCondi(rho,mass,u,iflag,f,fpost,xDim,yDim,xFluid,yFluid,w) 
!
!   goto 301
    open (11, file='pressure.txt',status='unknown')
    open (12, file='masstotal.txt',status='unknown')
	
	!循环开始--------------------------------------------------------------------------------------------------
    DO 300 tStep = 1, tMax
       CALL massUpdate(fpost,rho,mass,iflag,xDim,yDim,e,inv)
       CALL stream(f,fpost,mass,iflag,u,afb,aft,afr,afl,adb,adt,adr,adl,xDIm,yDim,rhoA,e,w,inv,ibslip,itslip,irslip,ilslip)
       CALL collide(f,fpost,iflag,rho,u,xDim,yDim,omega,gravn,w,e,p_lb,mass,P_lt2ph)
       CALL cellsUpdate(f,rho,u,iflag,mass,xDim,yDim,w,e)     
       
!     //Print out result
        IF(mod(tStep,100) .eq. 0)then
         print*,'Time Step = ',tStep      
!         CALL output(rho,u,iflag,mass,tStep,xDim, yDim)
        ENDIF
        
        times=dt*dble(tstep) ! dble 函数,把数据转换成双精度。
        write(11,*) times,p_lb       
!       
 
       do x=1,xDim
       do y=1,yDim
         do i=0,8
         fpost(x,y,i)=f(x,y,i) !每次循环完了就把f储存在fpost以待调用
         enddo
       enddo
       enddo
!   //Checking mass conservation
    mass1=0.d0
    mass2=0.d0
       do x=1,xDim
       do y=1,yDim
       mass1=mass1+mass(x,y) !总质量
       if (iflag(x,y).eq.1.or.iflag(x,y).eq.2) mass2=mass2+mass(x,y)  
       enddo
       enddo
       
      if (isnan(mass1)) then
        Print*,'THE NaN ERROR OCCURS'
        PAUSE
        STOP
      endif
       
   Write(12,116)tStep, mass1,mass2,mass1-mass2
116 format(i15,3f20.5)
!
300 CONTINUE
    close(11)
! 
301 continue   
!    call run_gnuplot ( command_file_name )
!    
    END PROGRAM dambreak
      
!-----------Equilibrium Function-------------------------------------------------------      
      Real*8 function F_eq(wi,rhoxy,uxy,uSqr)
      implicit none
      Real*8:: uxy,uSqr,wi,rhoxy
      
      F_eq = wi*rhoxy*(1.d0+3.0d0*uxy+4.5d0*uxy*uxy-1.5d0*uSqr)
      
      end function    
!-------------------------------------------------------------------------------------
!-----Reading parameters--------------------------------------------------------------
    SUBROUTINE Readparameters(tMax,BCTYPE,dx,dt,vis_phy,gra_phy,den_phy)
!
    implicit none
	Integer::tMax,BCTYPE
	Real*8:: dx,dt,vis_phy,gra_phy,den_phy
!
	open(1,file='PARAMETERS.txt')
	read(1,*)
	read(1,*)
	read(1,*)tMax
	read(1,*)
	read(1,*)dx,dt
	read(1,*)
	read(1,*)BCTYPE
	read(1,*)
	read(1,*)vis_phy
	read(1,*)
	read(1,*)gra_phy
	read(1,*)
	read(1,*)den_phy
	close(1)
    
	END SUBROUTINE
!--------------------------------------------------------------------------------------      
!------------------SUBROUTINE: iniCondi------------------------------------------------
!     Function: Generating Initial conditions
!      
      SUBROUTINE iniCondi(rho,mass,u,iflag,f,fpost,xDim,yDim,xFluid,yFluid,w)
      implicit none         
!      
      Integer:: xDim,yDim,xFluid,yFluid,x,y,i
      Integer:: iflag(0:xDim+1,0:yDim+1)
      Real*8 :: rho(1:xDim,1:yDim),mass(1:xDim,1:yDim),u(1:xDim,1:yDim,0:1)
	  Real*8 :: f(1:xDim,1:yDim,0:8),fpost(1:xDim,1:yDim,0:8)
      Real*8 :: w(0:8)    
!     //For fluid cells
      do x=1,xFluid
        do y=1,yFluid
          rho(x,y)  = 1.d0
          iflag(x,y)= 1
          mass(x,y) = rho(x,y)
          u(x,y,0)  = 0.d0
          u(x,y,1)  = 0.d0
          do i = 0, 8
            f(x,y,i)=w(i)
            fpost(x,y,i)=w(i)
          end do
        enddo
      enddo      
!     //For surface cells
      y=yFluid+1 !对于上表面，每个点的mass=流体的一半&#x