fluid_solve.f90

visual fortran 源码

!-------------块修正TDMA求解器----------------------------
subroutine solve()
use varible 
integer :: ISTF,JSTF,IT1,IT2,JT1,JT2
real   :: BL,BLP,BLM,BLC,DENOM
integer :: I,J,K,II,JJ
!write(*,*)"ist=",ist
!write(*,*)"jst=",jst
ISTF=IST-1
JSTF=JST-1
IT1=L2+IST
IT2=L3+IST
JT1=M2+JST
JT2=M3+JST
!-------------------------------------------------------------
do NT=1,NTIMES(NF)
 N=NF
!--------------------块修正器---------------------------------
 if(LBLK(NF)==.TRUE.)then
  PT(ISTF)=0
  QT(ISTF)=0
    do I=IST,L2
    BL=0
	BLP=0
	BLM=0
	BLC=0
       do J=JST,M2
       BL=BL+AP(I,J)
       if(J/=M2)then
	   BL=BL-AJP(I,J)
  	   end if
	   if(J/=JST)then
	   BL=BL-AJM(I,J)
	   end if
       BLP=BLP+AIP(I,J)
	   BLM=BLM+AIM(I,J)
       BLC=BLC+CON(I,J)+AIP(I,J)*F(I+1,J,N)+AIM(I,J)*F(I-1,J,N)+AJP(I,J)*F(I,J+1,N)+AJM(I,J)*F(I,J-1,N)-AP(I,J)*F(I,J,N)
	   continue
	   end do
    DENOM=BL-PT(I-1)*BLM
    if(ABS(DENOM/BL)<1E-10)then
    DENOM=1E35
    end if
    PT(I)=BLP/DENOM
    QT(I)=(BLC+BLM*QT(I-1))/DENOM
    continue
	end do
  !-----------------------------------------------------
  BL=0 
    do II=IST,L2
    I=IT1-II
    BL=BL*PT(I)+QT(I)
      do J=JST,M2
	  F(I,J,N)=F(I,J,N)+BL
	  end do
	end do
!-----------------------------------------------------
  PT(JSTF)=0
  QT(JSTF)=0
    do J=JST,M2
       BL=0
	   BLP=0
	   BLM=0
	   BLC=0
	   do I=IST,L2
	   BL=BL+AP(I,J)
	   if(I/=L2)then
	   BL=BL-AIP(I,J)
	   end if
	   if (I/=IST)then
	   BL=BL-AIM(I,J)
	   end if 
	   BLP=BLP+AJP(I,J)
	   BLM=BLM+AJM(I,J)
	   BLC=BLC+CON(I,J)+AIP(I,J)*F(I+1,J,N)+AIM(I,J)*F(I-1,J,N)+AJP(I,J)*F(I,J+1,N)+AJM(I,J)*F(I,J-1,N)-AP(I,J)*F(I,J,N)
	   continue
	   end do
	 DENOM=BL-PT(J-1)*BLM
	 if(ABS(DENOM/BL)<1E-10)then
	 DENOM=1E35
	 end if
	 PT(J)=BLP/DENOM
	 QT(J)=(BLC+BLM*QT(J-1))/DENOM
	 continue
	 end do

!--------------------------------------------------------
  BL=0
    do JJ=JST,M2
	J=JT1-JJ
	BL=BL*PT(J)+QT(J)
	  do I=IST,L2
	  F(I,J,N)=F(I,J,N)+BL
	  end do
	end do
 else
 continue
 end if
!----------------TDMA交替方向迭代-------------------
    do J=JST,M2
    PT(ISTF)=0
    QT(ISTF)=F(ISTF,J,N)
      do I=IST,L2
	  DENOM=AP(I,J)-PT(I-1)*AIM(I,J)
	  PT(I)=AIP(I,J)/DENOM
	  TEMP=CON(I,J)+AJP(I,J)*F(I,J+1,N)+AJM(I,J)*F(I,J-1,N)
	  QT(I)=(TEMP+AIM(I,J)*QT(I-1))/DENOM
	  continue
	  end do
	  do II=IST,L2
	  I=IT1-II
      F(I,J,N)=F(I+1,J,N)*PT(I)+QT(I)
	  end do
    continue
    end do
!---------------------------------------------------------
    do JJ=JST,M3
	J=JT2-JJ
	PT(ISTF)=0
	QT(ISTF)=F(ISTF,J,N)
	  do I=IST,L2
	  DENOM=AP(I,J)-PT(I-1)*AIM(I,J)
	  PT(I)=AIP(I,J)/DENOM
	  TEMP=CON(I,J)+AJP(I,J)*F(I,J+1,N)+AJM(I,J)*F(I,J-1,N)
	  QT(I)=(TEMP+AIM(I,J)*QT(I-1))/DENOM
      continue
	  end do
	  do II=IST,L2
	  I=IT1-II
	  F(I,J,N)=F(I+1,J,N)*PT(I)+QT(I)
      end do
	continue
	end do
!----------------------------------------------------------
    do I=IST,L2
	PT(JSTF)=0
	QT(JSTF)=F(I,JSTF,N)
	  do J=JST,M2
	  DENOM=AP(I,J)-PT(J-1)*AJM(I,J)
	  PT(J)=AJP(I,J)/DENOM
	  TEMP=CON(I,J)+AIP(I,J)*F(I+1,J,N)+AIM(I,J)*F(I-1,J,N)
	  QT(J)=(TEMP+AJM(I,J)*QT(J-1))/DENOM
	  continue
	  end do
	  do JJ=JST,M2
	  J=JT1-JJ
	  F(I,J,N)=F(I,J+1,N)*PT(J)+QT(J)
	  end do 
	continue
	end do
!----------------------------------------------------------
    do II=IST,L3
	I=IT2-II
	PT(JSTF)=0
	QT(JSTF)=F(I,JSTF,N)
	  do J=JST,M2
	  DENOM=AP(I,J)-PT(J-1)*AJM(I,J)
	  PT(J)=AJP(I,J)/DENOM
	  TEMP=CON(I,J)+AIP(I,J)*F(I+1,J,N)+AIM(I,J)*F(I-1,J,N)
	  QT(J)=(TEMP+AJM(I,J)*QT(J-1))/DENOM
	  continue
	  end do
	  do JJ=JST,M2
	  J=JT1-JJ
	  F(I,J,N)=F(I,J+1,N)*PT(J)+QT(J)
	  end do
	continue 
	end do
!------------------------------------------------------------
continue 
end do
do J=2,M2
   do I=2,L2
   CON(I,J)=0
   AP(I,J)=0
   continue
   end do
end do
return 
end subroutine solve