fluid_setup.f90

visual fortran 源码

subroutine SETUP()
use varible
use non_common

!===========================================================================================
!==================================ENTRY-SETUP1=============================================
!===========================================================================================

!--------------------设置网格系统---------------------------------
ENTRY SETUP1()
 NP=NFMAX+1
 NRHO=NP+1
 NGAM=NRHO+1
!---------------------轴设置--------------------------------------
 X(1)=XU(2)
 do I=2,L2
 X(I)=0.5*(XU(I+1)+XU(I))
 end do
 X(L1)=XU(L1)
 Y(1)=YV(2)
 do J=2,M2
 Y(J)=0.5*(YV(J+1)+YV(J))
 end do
 Y(M1)=YV(M1)
!-----------------X方向控制体设置----------------------------------------
 do I=2,L1
 XDIF(I)=X(I)-X(I-1)
 end do
 do I=2,L2
 XCV(I)=XU(I+1)-XU(I)
 end do
 do I=3,L2
 XCVS(I)=XDIF(I)
 end do
 XCVS(3)=XCVS(3)+XDIF(2)
 XCVS(L2)=XCVS(L2)+XDIF(L1)
 do I=3,L3
 XCVI(I)=0.5*XCV(I)
 XCVIP(I)=XCVI(I)
 end do
 XCVIP(2)=XCV(2)
 XCVI(L2)=XCV(L2)
!-----------------Y方向控制体设置----------------------------------------
 do J=2,M1
 YDIF(J)=Y(J)-Y(J-1)
 end do 
 do J=2,M2
 YCV(J)=YV(J+1)-YV(J)
 end do
 do J=3,M2
 YCVS(J)=YDIF(J)
 end do
 YCVS(3)=YCVS(3)+YDIF(2)
 YCVS(M2)=YCVS(M2)+YDIF(M1)
!-----------------------------半径的设置--------------------------------------
 if(MODE==1)then
   do J=1,M1
   RMN(J)=1.0
   R(J)=1.0
   end do
 else
   do J=2,M1
   R(J)=R(J-1)+YDIF(J)
   end do
   RMN(2)=R(1)
   do J=3,M2
   RMN(J)=RMN(J-1)+YCV(J-1)
   end do
   RMN(M1)=R(M1)
 end if 
 continue
!----------------------------东西尺度系数的设置-------------------------------
 do J=1,M1
 SX(J)=1.0
 SXMN(J)=1.0
 if(MODE==3)then
   SX(J)=R(J)
   if(J/=1)then
   SXMN(J)=RMN(J)
   end if
 else
 continue
 end if 
 end do
!-------------------------
 do J=2,M2
 YCVR(J)=R(J)*YCV(J)
 ARX(J)=YCVR(J)
 if(MODE==3)then
 ARX(J)=YCV(J)
 else 
 continue
 end if
 end do
 do J=4,M3
 YCVRS(J)=0.5*(R(J)+R(J-1))*YDIF(J)
 end do
 YCVRS(3)=0.5*(R(3)+R(1))*YCVS(3)
 YCVRS(M2)=0.5*(R(M1)+R(M3))*YCVS(M2)
!-----------------------------------------------------------------------------------
 if (MODE==2)then
 do J=3,M3
 ARXJ(J)=0.25*(1.0+RMN(J)/R(J))*ARX(J)
 ARXJP(J)=ARX(J)-ARXJ(J)
 end do
 else
 do J=3,M3
 ARXJ(J)=0.5*ARX(J)
 ARXJP(J)=ARXJ(J)
 end do
 end if
 ARXJP(2)=ARX(2)
 ARXJ(M2)=ARX(M2)
!-------------------------------插值系数------------------------------------------------
 do J=3,M3
 FV(J)=ARXJP(J)/ARX(J)
 FVP(J)=1.0-FV(J)
 end do
 do I=3,L2
 FX(I)=0.5*XCV(I-1)/XDIF(I)
 FXM(I)=1.0-FX(I)
 end do
 FX(2)=0.0
 FXM(2)=1.0
 FX(L1)=1.0
 FXM(L1)=0.0
 do J=3,M2
 FY(J)=0.5*YCV(J-1)/YDIF(J)
 FYM(J)=1.0-FY(J)
 end do
 FY(2)=0.0
 FYM(2)=1.0
 FY(M1)=1.0
 FYM(M1)=0.0
!-----------------------------------INITIALIZE-------------------------------------------
 do J=1,M1
    do I=1,L1
    PC(I,J)=0.0
    U(I,J)=0.0
    V(I,J)=0.0
    CON(I,J)=0.0
    AP(I,J)=0.0
    RHO(I,J)=RHOCON
    P(I,J)=0.0
    continue
    end do 
 end do
return

!===========================================================================================
!===================================ENTRY-SETUP2============================================
!===========================================================================================
!--------------COEFFICIENTS FOR THE U EQUATION----------------------------------------------
ENTRY SETUP2()
  NF=1
if(LSOLVE(NF)==.TRUE.)then
  IST=3
  JST=2
  call GAMSOR()
  REL=1.0-RELAX(NF)
  do I=3,L2
  FL=XCVI(I)*V(I,2)*RHO(I,1)
  FLM=XCVIP(I-1)*V(I-1,2)*RHO(I-1,1)
  FLOW=R(1)*(FL+FLM)
  DIFF=R(1)*(XCVI(I)*GAM(I,1)+XCVIP(I-1)*GAM(I-1,1))/YDIF(2)
  call DIFLOW()
  AJM(I,2)=ACOF+AMAX1(0.0, FLOW)
  end do
  !-----------------------AIM(I,J),AIP(I,J)-------------------------------
  do J=2,M2
  FLOW=ARX(J)*U(2,J)*RHO(1,J)
  DIFF=ARX(J)*GAM(1,J)/(XCV(2)*SX(J))
  call DIFLOW()
  AIM(3,J)=ACOF+AMAX1(0.0,FLOW)
    do I=3,L2
	if(I==L2)then
    FLOW=ARX(J)*U(L1,J)*RHO(L1,J)
	DIFF=ARX(J)*GAM(L1,J)/(XCV(L2)*SX(J))
	else
	FL=U(I,J)*(FX(I)*RHO(I,J)+FXM(I)*RHO(I-1,J))
	FLP=U(I+1,J)*(FX(I+1)*RHO(I+1,J)+FXM(I+1)*RHO(I,J))
	FLOW=ARX(J)*0.5*(FL+FLP)
	DIFF=ARX(J)*GAM(I,J)/(XCV(I)*SX(J))
	end if
	call DIFLOW()
	AIM(I+1,J)=ACOF+AMAX1(0.0,FLOW)
	AIP(I,J)=AIM(I+1,J)-FLOW
	!-----------------------AJM(I,J),AJP(I,J)-------------------------------
	if(J==M2)then
	FL=XCVI(I)*V(I,M1)*RHO(I,M1)
	FLM=XCVIP(I-1)*V(I-1,M1)*RHO(I-1,M1)
	DIFF=R(M1)*(XCVI(I)*GAM(I,M1)+XCVIP(I-1)*GAM(I-1,M1))/YDIF(M1)
	else
	FL=XCVI(I)*V(I,J+1)*(FY(J+1)*RHO(I,J+1)+FYM(J+1)*RHO(I,J))
	FLM=XCVIP(I-1)*V(I-1,J+1)*(FY(J+1)*RHO(I-1,J+1)+FYM(J+1)*RHO(I-1,J))
	GM=GAM(I,J)*GAM(I,J+1)/(YCV(J)*GAM(I,J+1)+YCV(J+1)*GAM(I,J)+1.0E-30)*XCVI(I)
	GMM=GAM(I-1,J)*GAM(I-1,J+1)/(YCV(J)*GAM(I-1,J+1)+YCV(J+1)*GAM(I-1,J)+1.0E-30)*XCVIP(I-1)
	DIFF=RMN(J+1)*2*(GM+GMM)
	end if
	FLOW=RMN(J+1)*(FL+FLM)
	call DIFLOW()
	AJM(I,J+1)=ACOF+AMAX1(0.0,FLOW)
	AJP(I,J)=AJM(I,J+1)-FLOW
	!---------------------AP(I,J),CON(I,J),DU(I,J)--------------------------
	VOL=YCVR(J)*XCVS(I)
	APT=(RHO(I,J)*XCVI(I)+RHO(I-1,J)*XCVIP(I-1))/(XCVS(I)*DT)
	AP(I,J)=AP(I,J)-APT
	CON(I,J)=CON(I,J)+APT*U(I,J)
	AP(I,J)=(-AP(I,J)*VOL+AIP(I,J)+AIM(I,J)+AJP(I,J)+AJM(I,J))/RELAX(NF)
	CON(I,J)=CON(I,J)*VOL+REL*AP(I,J)*U(I,J)
	DU(I,J)=VOL/(XDIF(I)*SX(J))
	CON(I,J)=CON(I,J)+DU(I,J)*(P(I-1,J)-P(I,J))
	DU(I,J)=DU(I,J)/AP(I,J)
	end do
  end do
  call solve()

 ! do K=1,6
  !  do J=JST,M2
 !     do I=IST,L2
!	  COFU(I,J,K)=COF(I,J,K)
!	  end do
!	end do
 ! end do
else 
continue
end if


  !--------------------COEFFICIENTS FOR THE U EQUATION-----------------------
  NF=2
if(LSOLVE(NF)==.TRUE.)then
  IST=2
  JST=3
  call GAMSOR()
  REL=1.0-RELAX(NF)
  do I=2,L2
  AREA=R(1)*XCV(I)
  FLOW=AREA*V(I,2)*RHO(I,1)
  DIFF=AREA*GAM(I,1)/YCV(2)
  call DIFLOW()
  AJM(I,3)=ACOF+AMAX1(0.0,FLOW)
  end do
  !-----------------------AIM(I,J),AIP(I,J)-------------------------------
  do J=3,M2
  FL=ARXJ(J)*U(2,J)*RHO(1,J)
  FLM=ARXJP(J-1)*U(2,J-1)*RHO(1,J-1)
  FLOW=FL+FLM
  DIFF=(ARXJ(J)*GAM(1,J)+ARXJP(J-1)*GAM(1,J-1))/(XDIF(2)*SXMN(J))
  call DIFLOW()
  AIM(2,J)=ACOF+AMAX1(0.0,FLOW)
    do I=2,L2
	if(I==L2)then
	FL=ARXJ(J)*U(L1,J)*RHO(L1,J)
	FLM=ARXJP(J-1)*U(L1,J-1)*RHO(L1,J-1)
	DIFF=(ARXJ(J)*GAM(L1,J)+ARXJP(J-1)*GAM(L1,J-1))/(XDIF(L1)*SXMN(J))
	else
	FL=ARXJ(J)*U(I+1,J)*(FX(I+1)*RHO(I+1,J)+FXM(I+1)*RHO(I,J))
	FLM=ARXJP(J-1)*U(I+1,J-1)*(FX(I+1)*RHO(I+1,J-1)+FXM(I+1)*RHO(I,J-1))
	GM=GAM(I,J)*GAM(I+1,J)/(XCV(I)*GAM(I+1,J)+XCV(I+1)*GAM(I,J)+1.0E-30)*ARXJ(J)
	GMM=GAM(I,J-1)*GAM(I+1,J-1)/(XCV(I)*GAM(I+1,J-1)+XCV(I+1)*GAM(I,J-1)+1.0E-30)*ARXJP(J-1)
	DIFF=2.0*(GM+GMM)/SXMN(J)
	end if
	FLOW=FL+FLM
	call DIFLOW()
	AIM(I+1,J)=ACOF+AMAX1(0.0,FLOW)
	AIP(I,J)=AIM(I+1,J)-FLOW
	!-----------------------AJM(I,J),AJP(I,J)-------------------------------
	If(J==M2)then
	AREA=R(M1)*XCV(I)
	FLOW=AREA*V(I,M1)*RHO(I,M1)
	DIFF=AREA*GAM(I,M1)/YCV(M2)
	else
	AREA=R(J)*XCV(I)
	FL=V(I,J)*(FY(J)*RHO(I,J)+FYM(J)*RHO(I,J-1))*RMN(J)
	FLP=V(I,J+1)*(FY(J+1)*RHO(I,J+1)+FYM(J+1)*RHO(I,J))*RMN(J+1)
	FLOW=(FV(J)*FL+FVP(J)*FLP)*XCV(I)
	DIFF=AREA*GAM(I,J)/YCV(J)
	end if
	call DIFLOW()
	AJM(I,J+1)=ACOF+AMAX1(0.0,FLOW)
	AJP(I,J)=AJM(I,J+1)-FLOW
	!---------------------AP(I,J),CON(I,J),DU(I,J)--------------------------
	VOL=YCVRS(J)*XCV(I)
	SXT=SX(J)
	if(J==M2)then
	SXT=SX(M1)
	end if
	SXB=SX(J-1)
	if(J==3)then
	SXB=SX(1)
	end if
	APT=(ARXJ(J)*RHO(I,J)*0.5*(SXT+SXMN(J))+ARXJP(J-1)*RHO(I,J-1)*0.5*(SXB+SXMN(J)))/(YCVRS(J)*DT)

	AP(I,J)=AP(I,J)-APT
	CON(I,J)=CON(I,J)+APT*V(I,J)
	AP(I,J)=(-AP(I,J)*VOL+AIP(I,J)+AIM(I,J)+AJP(I,J)+AJM(I,J))/RELAX(NF)
	CON(I,J)=CON(I,J)*VOL+REL*AP(I,J)*V(I,J)
	DV(I,J)=VOL/YDIF(J)
    CON(I,J)=CON(I,J)+DV(I,J)*(P(I,J-1)-P(I,J))
	DV(I,J)=DV(I,J)/AP(I,J)
	continue
	end do
  end do
  call solve()
 ! do K=1,6
 !   do J=JST,M2
 !     do I=IST,L2
!	  COFV(I,J,K)=COF(I,J,K)
!	  end do
!	end do
  !end do
else
continue
end if
  !--------------------COEFFICIENTS FOR THE P EQUATION-----------------------
  NF=3
if(LSOLVE(NF)==.TRUE.)then
  IST=2
  JST=2
  call GAMSOR()
  SMAX=0.
  SSUM=0.
  do J=2,M2
     do I=2,L2
	 VOL=YCVR(J)*XCV(I)
	 CON(I,J)=CON(I,J)*VOL
	 end do
  end do
  do I=2,L2
  ARHO=R(1)*XCV(I)*RHO(I,1)
  CON(I,2)=CON(I,2)+ARHO*V(I,2)
  AJM(I,2)=0
  end do 

  do J=2,M2
  ARHO=ARX(J)*RHO(1,J)
  CON(2,J)=CON(2,J)+ARHO*U(2,J)
  AIM(2,J)=0
     do I=2,L2
	 if(I==L2)then
	 ARHO=ARX(J)*RHO(L1,J)
	 CON(I,J)=CON(I,J)-ARHO*U(L1,J)
	 AIP(I,J)=0
	 else
	 ARHO=ARX(J)*(FX(I+1)*RHO(I+1,J)+FXM(I+1)*RHO(I,J))
	 FLOW=ARHO*U(I+1,J)
	 CON(I,J)=CON(I,J)-FLOW
	 CON(I+1,J)=CON(I+1,J)+FLOW
	 AIP(I,J)=ARHO*DU(I+1,J)
	 AIM(I+1,J)=AIP(I,J)
	 end if

	 if(J==M2)then
	 ARHO=RMN(M1)*XCV(I)*RHO(I,M1)
	 CON(I,J)=CON(I,J)-ARHO*V(I,M1)
	 AJP(I,J)=0
	 else
	 AHRO=RMN(J+1)*XCV(I)*(FY(J+1)*RHO(I,J+1)+FYM(J+1)*RHO(I,J))
	 FLOW=ARHO*V(I,J+1)
	 CON(I,J)=CON(I,J)-FLOW
	 CON(I,J+1)=CON(I,J+1)+FLOW
	 AJP(I,J)=ARHO*DV(I,J+1)
	 AJM(I,J+1)=AJP(I,J)
	 end if

	 AP(I,J)=AIP(I,J)+AIM(I,J)+AJP(I,J)+AJM(I,J)
	 PC(I,J)=0.
     SMAX=AMAX1(SMAX,ABS(CON(I,J)))
     SSUM=SSUM+CON(I,J)
	 continue
	 end do
  end do
  call solve()

  !--------------------速度与压力修正值的亚松弛---------------------------
  do J=2,M2
    do I=2,L2
    P(I,J)=P(I,J)+PC(I,J)*RELAX(NP)
    IF(I .NE. 2) then
	U(I,J)=U(I,J)+DU(I,J)*(PC(I-1,J)-PC(I,J))
    end if
	IF(J .NE. 2) then
	V(I,J)=V(I,J)+DV(I,J)*(PC(I,J-1)-PC(I,J))
	end if
    continue
    end do
  end do
else
continue
end if

!----------------COEFFICIENTS FOR OTHER EQUATIONS----------------------
IST=2
JST=2
do N=4,NFMAX
NF=N
!--------------------------next row:lsove(nf)==.true.,right or wrong?
if(LSOLVE(NF)==.TRUE.)then
call GAMSOR()
   REL=1.0-RELAX(NF)
   do I=2,L2
   AREA=R(1)*XCV(I)
   FLOW=AREA*V(I,2)*RHO(I,1)
   DIFF=AREA*GAM(I,1)/YDIF(2)
   call DIFLOW()
   !write(*,*)"acof=",acof
   AJM(I,2)=ACOF+AMAX1(0.0,FLOW)
   end do 
   do J=2,M2
   FLOW=ARX(J)*U(2,J)*RHO(1,J)
   DIFF=ARX(J)*GAM(1,J)/(XDIF(2)*SX(J))
   call DIFLOW()
   AIM(2,J)=ACOF+AMAX1(0.0,FLOW)
      do I=2,L2
	  if(I==L2)then
	  FLOW=ARX(J)*U(L1,J)*RHO(L1,J)
	  DIFF=ARX(J)*GAM(L1,J)/(XDIF(L1)*SX(J))
	  else
	  FLOW=ARX(J)*U(I+1,J)*(FX(I+1)*RHO(I+1,J)+FXM(I+1)*RHO(I,J))
	  DIFF=ARX(J)*2*GAM(I,J)*GAM(I+1,J)/((XCV(I)*GAM(I+1,J)+XCV(I+1)*GAM(I,J)+1.0E-30)*SX(J))
	  end if 
	  call DIFLOW()
	  AIM(I+1,J)=ACOF+AMAX1(0.0,FLOW)
	  AIP(I,J)=AIM(I+1,J)-FLOW
	  AREA=RMN(J+1)*XCV(I)
	  if(J==M2)then
	  FLOW=AREA*V(I,M1)*RHO(I,M1)
	  DIFF=AREA*GAM(I,M1)/YDIF(M1)
	  else
	  FLOW=AREA*V(I,J+1)*(FY(J+1)*RHO(I,J+1)+FYM(J+1)*RHO(I,J))
	  DIFF=AREA*2*GAM(I,J)*GAM(I,J+1)/(YCV(J)*GAM(I,J+1)+YCV(J+1)*GAM(I,J)+1.0E-30)
	  end if
	  call DIFLOW()
	  AJM(I,J+1)=ACOF+AMAX1(0.0,FLOW)
	  AJP(I,J)=AJM(I,J+1)-FLOW
	  VOL=YCVR(J)*XCV(I)
	  APT=RHO(I,J)/DT
	  AP(I,J)=AP(I,J)-APT
	  CON(I,J)=CON(I,J)+APT*F(I,J,NF)
	  AP(I,J)=(-AP(I,J)*VOL+AIP(I,J)+AIM(I,J)+AJP(I,J)+AJM(I,J))/RELAX(NF)
	  CON(I,J)=CON(I,J)*VOL+REL*AP(I,J)*F(I,J,NF)
	  continue
	  end do
   end do
   call SOLVE()
   !write(101,"(f8.4,f8.4,f8.4,f8.4,f8.4,f8.4,f8.4,f8.4,f8.4,f8.4,f8.4,f8.4,f8.4,f8.4)")((t(i,j),i=1,ni),j=1,nj)
else
continue
end if
end do
TIME=TIME+DT
!WRITE(*,*)ITER
ITER=ITER+1
if(ITER<LAST)then
LSTOP=.FALSE.
else
LSTOP=.TRUE.
end if
return

end subroutine SETUP