3dmain.for

这是我下的fortran的文件

   95 CONTINUE
      WRITE(8,1)
      RETURN
*----------------------------------------------------------------------
      ENTRY SETUP2
*---COEFFICIENTS FOR THE U EQUATION----
      NF=1
      IF(.NOT. LSOLVE(NF)) GO TO 100
      IST=3
      JST=2
      KST=2
      CALL GAMSOR
      REL=1.-RELAX(NF)
      DO 102 I=3,L2
      DO 103 K=2,N2
      FL=XCVI(I)*V(I,2,K)*RHO(I,1,K)
      FLM=XCVIP(I-1)*V(I-1,2,K)*RHO(I-1,1,K)
      FLOW=ZCV(K)*(FL+FLM)
      DIFF=ZCV(K)*(XCVI(I)*GAM(I,1,K)+XCVIP(I-1)*GAM(I-1,1,K))/YDIF(2)
      CALL DIFLOW
  103 AJM(I,2,K)=ACOF+AMAX1(0.,FLOW)
      DO 104 J=2,M2
      FL=XCVI(I)*W(I,J,2)*RHO(I,J,1)
      FLM=XCVIP(I-1)*W(I-1,J,2)*RHO(I-1,J,1)
      FLOW=YCV(J)*(FL+FLM)
      DIFF=YCV(J)*(XCVI(I)*GAM(I,J,1)+XCVIP(I-1)*GAM(I-1,J,1))/ZDIF(2)
      CALL DIFLOW
  104 AKM(I,J,2)=ACOF+AMAX1(0.,FLOW)
  102 CONTINUE
      DO 105 K=2,N2
      DO 105 J=2,M2
      FLOW=AX(J,K)*U(2,J,K)*RHO(1,J,K)
      DIFF=AX(J,K)*GAM(1,J,K)/XCV(2)
      CALL DIFLOW
      AIM(3,J,K)=ACOF+AMAX1(0.,FLOW)
      DO 105 I=3,L2
      IF(I .EQ. L2) GO TO 106
      FL=U(I,J,K)*(FX(I)*RHO(I,J,K)+FXM(I)*RHO(I-1,J,K))
      FLP=U(I+1,J,K)*(FX(I+1)*RHO(I+1,J,K)+FXM(I+1)*RHO(I,J,K))
      FLOW=AX(J,K)*0.5*(FL+FLP)
      DIFF=AX(J,K)*GAM(I,J,K)/XCV(I)
      GO TO 107
  106 FLOW=AX(J,K)*U(L1,J,K)*RHO(L1,J,K)
      DIFF=AX(J,K)*GAM(L1,J,K)/XCV(L2)
  107 CALL DIFLOW
      AIM(I+1,J,K)=ACOF+AMAX1(0.,FLOW)
      AIP(I,J,K)=AIM(I+1,J,K)-FLOW
      IF (J .EQ. M2) GOTO 108
      FL=XCVI(I)*V(I,J+1,K)*(FY(J+1)*RHO(I,J+1,K)+FYM(J+1)*RHO(I,J,K))
      FLM=XCVIP(I-1)*V(I-1,J+1,K)*(FY(J+1)*RHO(I-1,J+1,K)+FYM(J+1)*
     &   RHO(I-1,J,K))
      GM=GAM(I,J,K)*GAM(I,J+1,K)/(YCV(J)*GAM(I,J+1,K)
     &+YCV(J+1)*GAM(I,J,K)+1.0E-30)*XCVI(I)
      GMM=GAM(I-1,J,K)*GAM(I-1,J+1,K)/(YCV(J)*GAM(I-1,J+1,K)+YCV(J+1)*
     &   GAM(I-1,J,K)+1.E-30)*XCVIP(I-1)
      DIFF=ZCV(K)*2.*(GM+GMM)
      GO TO 109
  108 FL=XCVI(I)*V(I,M1,K)*RHO(I,M1,K)
      FLM=XCVIP(I-1)*V(I-1,M1,K)*RHO(I-1,M1,K)
      DIFF=ZCV(K)*(XCVI(I)*GAM(I,M1,K)+XCVIP(I-1)*GAM(I-1,M1,K))/
     &  YDIF(M1)
  109 FLOW=ZCV(K)*(FL+FLM)
      CALL DIFLOW
      AJM(I,J+1,K)=ACOF+AMAX1(0.,FLOW)
      AJP(I,J,K)=AJM(I,J+1,K)-FLOW
      IF (K .EQ. N2) GOTO 110
      FL=XCVI(I)*W(I,J,K+1)*(FZ(K+1)*RHO(I,J,K+1)+FZM(K+1)*RHO(I,J,K))
      FLM=XCVIP(I-1)*W(I-1,J,K+1)*(FZ(K+1)*RHO(I-1,J,K+1)+FZM(K+1)*
     &   RHO(I-1,J,K))
      GM=GAM(I,J,K)*GAM(I,J,K+1)/(ZCV(K)*GAM(I,J,K+1)+
     &ZCV(K+1)*GAM(I,J,K)+1.0E-30)*XCVI(I)
      GMM=GAM(I-1,J,K)*GAM(I-1,J,K+1)/(ZCV(K)*GAM(I-1,J,K+1)+ZCV(K+1)*
     &   GAM(I-1,J,K)+1.E-30)*XCVIP(I-1)
      DIFF=YCV(J)*2.*(GM+GMM)
      GO TO 111
  110 FL=XCVI(I)*W(I,J,N1)*RHO(I,J,N1)
      FLM=XCVIP(I-1)*W(I-1,J,N1)*RHO(I-1,J,N1)
      DIFF=YCV(J)*(XCVI(I)*GAM(I,J,N1)+XCVIP(I-1)*GAM(I-1,J,N1))/
     & ZDIF(N1)
  111 FLOW=YCV(J)*(FL+FLM)
      CALL DIFLOW
      AKM(I,J,K+1)=ACOF+AMAX1(0.,FLOW)
      AKP(I,J,K)=AKM(I,J,K+1)-FLOW
      VOL=AX(J,K)*XCVS(I)
      APT=(RHO(I,J,K)*XCVI(I)+RHO(I-1,J,K)*XCVIP(I-1))
     &   /(XCVS(I)*DT)
      AP(I,J,K)=AP(I,J,K)-APT
      CON(I,J,K)=CON(I,J,K)+APT*U(I,J,K)
      AP(I,J,K)=(-AP(I,J,K)*VOL+AIP(I,J,K)+AIM(I,J,K)+AJP(I,J,K)+
     &AJM(I,J,K)+AKP(I,J,K)+AKM(I,J,K))/RELAX(NF)
      CON(I,J,K)=CON(I,J,K)*VOL+REL*AP(I,J,K)*U(I,J,K)
      DU(I,J,K)=VOL/XDIF(I)
      CON(I,J,K)=CON(I,J,K)+DU(I,J,K)*(P(I-1,J,K)-P(I,J,K))
      DU(I,J,K)=DU(I,J,K)/AP(I,J,K)
  105 CONTINUE
      CALL SOLVE
  100 CONTINUE
*---COEFFICIENTS FOR THE V EQUATION----
      NF=2
      IF(.NOT. LSOLVE(NF)) GO TO 200
      IST=2
      JST=3
      KST=2
      CALL GAMSOR
      REL=1.-RELAX(NF)
      DO 202 I=2,L2
      DO 203 K=2,N2
      FLOW=AY(I,K)*V(I,2,K)*RHO(I,1,K)
      DIFF=AY(I,K)*GAM(I,1,K)/YCV(2)
      CALL DIFLOW
  203 AJM(I,3,K)=ACOF+AMAX1(0.,FLOW)
      DO 204 J=3,M2
      FL=YCVJ(J)*W(I,J,2)*RHO(I,J,1)
      FLM=YCVJP(J-1)*W(I,J-1,2)*RHO(I,J-1,1)
      FLOW=XCV(I)*(FL+FLM)
      DIFF=XCV(I)*(YCVJ(J)*GAM(I,J,1)+YCVJP(J-1)*GAM(I,J-1,1))/ZDIF(2)
      CALL DIFLOW
  204 AKM(I,J,2)=ACOF+AMAX1(0.,FLOW)
  202 CONTINUE
      DO 205 K=2,N2
      DO 205 J=3,M2
      FL=YCVJ(J)*U(2,J,K)*RHO(1,J,K)
      FLM=YCVJP(J-1)*U(2,J-1,K)*RHO(1,J-1,K)
      FLOW=(FL+FLM)*ZCV(K)
      DIFF=ZCV(K)*(YCVJ(J)*GAM(1,J,K)+YCVJP(J-1)*GAM(1,J-1,K))/XDIF(2)
      CALL DIFLOW
      AIM(2,J,K)=ACOF+AMAX1(0.,FLOW)
      DO 205 I=2,L2
      IF(I .EQ. L2)GO TO 206
      FL=YCVJ(J)*U(I+1,J,K)*(FX(I+1)*RHO(I+1,J,K)+FXM(I+1)*RHO(I,J,K))
      FLM=YCVJP(J-1)*U(I+1,J-1,K)*(FX(I+1)*RHO(I+1,J-1,K)+FXM(I+1)*
     &   RHO(I,J-1,K))
      GM=GAM(I,J,K)*GAM(I+1,J,K)/(XCV(I)*GAM(I+1,J,K)+
     &XCV(I+1)*GAM(I,J,K)+1.E-30)*YCVJ(J)
      GMM=GAM(I,J-1,K)*GAM(I+1,J-1,K)/(XCV(I)*GAM(I+1,J-1,K)+XCV(I+1)*
     &   GAM(I,J-1,K)+1.0E-30)*YCVJP(J-1)
      DIFF=2.*(GM+GMM)*ZCV(K)
      GO TO 207
  206 FL=YCVJ(J)*U(L1,J,K)*RHO(L1,J,K)
      FLM=YCVJP(J-1)*U(L1,J-1,K)*RHO(L1,J-1,K)
      DIFF=ZCV(K)*(YCVJ(J)*GAM(L1,J,K)+YCVJP(J-1)*GAM(L1,J-1,K))/
     & XDIF(L1)
  207 FLOW=(FL+FLM)*ZCV(K)
      CALL DIFLOW
      AIM(I+1,J,K)=ACOF+AMAX1(0.,FLOW)
      AIP(I,J,K)=AIM(I+1,J,K)-FLOW
      IF(J .EQ. M2) GO TO 208
      FL=V(I,J,K)*(FY(J)*RHO(I,J,K)+FYM(J)*RHO(I,J-1,K))
      FLP=V(I,J+1,K)*(FY(J+1)*RHO(I,J+1,K)+FYM(J+1)*RHO(I,J,K))
      FLOW=0.5*AY(I,K)*(FL+FLP)
      DIFF=AY(I,K)*GAM(I,J,K)/YCV(J)
      GO TO 209
  208 FLOW=AY(I,K)*V(I,M1,K)*RHO(I,M1,K)
      DIFF=AY(I,K)*GAM(I,M1,K)/YCV(M2)
  209 CALL DIFLOW
      AJM(I,J+1,K)=ACOF+AMAX1(0.,FLOW)
      AJP(I,J,K)=AJM(I,J+1,K)-FLOW
      IF (K .EQ. N2) GOTO 210
      FL=YCVJ(J)*W(I,J,K+1)*(FZ(K+1)*RHO(I,J,K+1)+FZM(K+1)*RHO(I,J,K))
      FLM=YCVJP(J-1)*W(I,J-1,K+1)*(FZ(K+1)*RHO(I,J-1,K+1)+FZM(K+1)*
     &   RHO(I,J-1,K))
      GM=GAM(I,J,K)*GAM(I,J,K+1)/(ZCV(K)*GAM(I,J,K+1)+
     &ZCV(K+1)*GAM(I,J,K)+1.0E-30)*YCVJ(J)
      GMM=GAM(I,J-1,K)*GAM(I,J-1,K+1)/(ZCV(K)*GAM(I,J-1,K+1)+ZCV(K+1)*
     &   GAM(I,J-1,K)+1.E-30)*YCVJP(J-1)
      DIFF=XCV(I)*2.*(GM+GMM)
      GO TO 211
  210 FL=YCVJ(J)*W(I,J,N1)*RHO(I,J,N1)
      FLM=YCVJP(J-1)*W(I,J-1,N1)*RHO(I,J-1,N1)
      DIFF=XCV(I)*(YCVJ(J)*GAM(I,J,N1)+YCVJP(J-1)*GAM(I,J-1,N1))/
     & ZDIF(N1)
  211 FLOW=XCV(I)*(FL+FLM)
      CALL DIFLOW
      AKM(I,J,K+1)=ACOF+AMAX1(0.,FLOW)
      AKP(I,J,K)=AKM(I,J,K+1)-FLOW
      VOL=AY(I,K)*YCVS(J)
      APT=(RHO(I,J,K)*YCVJ(J)+RHO(I,J-1,K)*YCVJP(J-1))
     &   /(YCVS(J)*DT)
      AP(I,J,K)=AP(I,J,K)-APT
      CON(I,J,K)=CON(I,J,K)+APT*V(I,J,K)
      AP(I,J,K)=(-AP(I,J,K)*VOL+AIP(I,J,K)+AIM(I,J,K)+AJP(I,J,K)+
     &AJM(I,J,K)+AKP(I,J,K)+AKM(I,J,K))/RELAX(NF)
      CON(I,J,K)=CON(I,J,K)*VOL+REL*AP(I,J,K)*V(I,J,K)
      DV(I,J,K)=VOL/YDIF(J)
      CON(I,J,K)=CON(I,J,K)+DV(I,J,K)*(P(I,J-1,K)-P(I,J,K))
      DV(I,J,K)=DV(I,J,K)/AP(I,J,K)
  205 CONTINUE
      CALL SOLVE
  200 CONTINUE
*---COEFFICIENTS FOR THE W EQUATION----
      NF=3
      IF(.NOT. LSOLVE(NF)) GO TO 300
      IST=2
      JST=2
      KST=3
      CALL GAMSOR
      REL=1.-RELAX(NF)
      DO 302 I=2,L2
      DO 303 K=3,N2
      FL=ZCVK(K)*V(I,2,K)*RHO(I,1,K)
      FLM=ZCVKP(K-1)*V(I,2,K-1)*RHO(I,1,K-1)
      FLOW=XCV(I)*(FL+FLM)
      DIFF=XCV(I)*(ZCVK(K)*GAM(I,1,K)+ZCVKP(K-1)*GAM(I,1,K-1))/YDIF(2)
      CALL DIFLOW
  303 AJM(I,2,K)=ACOF+AMAX1(0.,FLOW)
      DO 304 J=2,M2
      FLOW=AZ(I,J)*W(I,J,2)*RHO(I,J,1)
      DIFF=AZ(I,J)*GAM(I,J,1)/ZCV(2)
      CALL DIFLOW
  304 AKM(I,J,3)=ACOF+AMAX1(0.,FLOW)
  302 CONTINUE
      DO 305 K=3,N2
      DO 305 J=2,M2
      FL=ZCVK(K)*U(2,J,K)*RHO(1,J,K)
      FLM=ZCVKP(K-1)*U(2,J,K-1)*RHO(1,J,K-1)
      FLOW=(FL+FLM)*YCV(J)
      DIFF=YCV(J)*(ZCVK(K)*GAM(1,J,K)+ZCVKP(K-1)*GAM(1,J,K-1))/XDIF(2)
      CALL DIFLOW
      AIM(2,J,K)=ACOF+AMAX1(0.,FLOW)
      DO 305 I=2,L2
      IF(I .EQ. L2)GO TO 306
      FL=ZCVK(K)*U(I+1,J,K)*(FX(I+1)*RHO(I+1,J,K)+FXM(I+1)*RHO(I,J,K))
      FLM=ZCVKP(K-1)*U(I+1,J,K-1)*(FX(I+1)*RHO(I+1,J,K-1)+FXM(I+1)*
     &   RHO(I,J,K-1))
      GM=GAM(I,J,K)*GAM(I+1,J,K)/(XCV(I)*GAM(I+1,J,K)+
     &XCV(I+1)*GAM(I,J,K)+1.E-30)*ZCVK(K)
      GMM=GAM(I,J,K-1)*GAM(I+1,J,K-1)/(XCV(I)*GAM(I+1,J,K-1)+XCV(I+1)*
     &   GAM(I,J,K-1)+1.0E-30)*ZCVKP(K-1)
      DIFF=2.*(GM+GMM)*YCV(J)
      GO TO 307
  306 FL=ZCVK(K)*U(L1,J,K)*RHO(L1,J,K)
      FLM=ZCVKP(K-1)*U(L1,J,K-1)*RHO(L1,J,K-1)
      DIFF=YCV(J)*(ZCVK(K)*GAM(L1,J,K)+ZCVKP(K-1)*GAM(L1,J,K-1))/
     & XDIF(L1)
  307 FLOW=(FL+FLM)*YCV(J)
      CALL DIFLOW
      AIM(I+1,J,K)=ACOF+AMAX1(0.,FLOW)
      AIP(I,J,K)=AIM(I+1,J,K)-FLOW
      IF(K .EQ. N2) GO TO 308
      FL=W(I,J,K)*(FZ(K)*RHO(I,J,K)+FZM(K)*RHO(I,J,K-1))
      FLP=W(I,J,K+1)*(FZ(K+1)*RHO(I,J,K+1)+FZM(K+1)*RHO(I,J,K))
      FLOW=0.5*AZ(I,J)*(FL+FLP)
      DIFF=AZ(I,J)*GAM(I,J,K)/ZCV(K)
      GO TO 309
  308 FLOW=AZ(I,J)*W(I,J,N1)*RHO(I,J,N1)
      DIFF=AZ(I,J)*GAM(I,J,N1)/ZCV(N2)
  309 CALL DIFLOW
      AKM(I,J,K+1)=ACOF+AMAX1(0.,FLOW)
      AKP(I,J,K)=AKM(I,J,K+1)-FLOW
      IF (J .EQ. M2) GOTO 310
      FL=ZCVK(K)*V(I,J+1,K)*(FY(J+1)*RHO(I,J+1,K)+FYM(J+1)*RHO(I,J,K))
      FLM=ZCVKP(K-1)*V(I,J+1,K-1)*(FY(J+1)*RHO(I,J+1,K-1)+FYM(J+1)*
     &   RHO(I,J,K-1))
      GM=GAM(I,J,K)*GAM(I,J+1,K)/(YCV(J)*GAM(I,J+1,K)+
     &YCV(J+1)*GAM(I,J,K)+1.0E-30)*ZCVK(K)
      GMM=GAM(I,J,K-1)*GAM(I,J+1,K-1)/(YCV(J)*GAM(I,J+1,K-1)+YCV(J+1)*
     &   GAM(I,J,K-1)+1.E-30)*ZCVKP(K-1)
      DIFF=XCV(I)*2.*(GM+GMM)
      GO TO 311
  310 FL=ZCVK(K)*V(I,M1,K)*RHO(I,M1,K)
      FLM=ZCVKP(K-1)*V(I,M1,K-1)*RHO(I,M1,K-1)
      DIFF=XCV(I)*(ZCVK(K)*GAM(I,M1,K)+ZCVKP(K-1)*GAM(I,M1,K-1))/
     & YDIF(M1)
  311 FLOW=XCV(I)*(FL+FLM)
      CALL DIFLOW
      AJM(I,J+1,K)=ACOF+AMAX1(0.,FLOW)
      AJP(I,J,K)=AJM(I,J+1,K)-FLOW
      VOL=AZ(I,J)*ZCVS(K)
      APT=(RHO(I,J,K)*ZCVK(K)+RHO(I,J,K-1)*ZCVKP(K-1))
     &   /(ZCVS(K)*DT)
      AP(I,J,K)=AP(I,J,K)-APT
      CON(I,J,K)=CON(I,J,K)+APT*W(I,J,K)
      AP(I,J,K)=(-AP(I,J,K)*VOL+AIP(I,J,K)+AIM(I,J,K)+AJP(I,J,K)+
     &AJM(I,J,K)+AKP(I,J,K)+AKM(I,J,K))/RELAX(NF)
      CON(I,J,K)=CON(I,J,K)*VOL+REL*AP(I,J,K)*W(I,J,K)
      DW(I,J,K)=VOL/ZDIF(K)
      CON(I,J,K)=CON(I,J,K)+DW(I,J,K)*(P(I,J,K-1)-P(I,J,K))
      DW(I,J,K)=DW(I,J,K)/AP(I,J,K)
  305 CONTINUE
      CALL SOLVE
  300 CONTINUE
*---COEFIICIENTS FOR THE PRESSURE CORRECTION EQUATION----
      NF=4
      IF(.NOT. LSOLVE(NF)) GO TO 500
      IST=2
      JST=2
      KST=2
      CALL GAMSOR
      SMAX=0.
      SSUM=0.
      DO 410 K=2,N2
      DO 410 J=2,M2
      DO 410 I=2,L2
      VOL=AX(J,K)*XCV(I)
  410 CON(I,J,K)=CON(I,J,K)*VOL
      DO 400 I=2,L2
      DO 401 K=2,N2
      ARHO=AY(I,K)*RHO(I,1,K)
      CON(I,2,K)=CON(I,2,K)+ARHO*V(I,2,K)
  401 AJM(I,2,K)=0.
      DO 402 J=2,M2
      ARHO=AZ(I,J)*RHO(I,J,1)
      CON(I,J,2)=CON(I,J,2)+ARHO*W(I,J,2)
  402 AKM(I,J,2)=0.
  400 CONTINUE
      DO 403 K=2,N2
      DO 403 J=2,M2
      ARHO=AX(J,K)*RHO(1,J,K)
      CON(2,J,K)=CON(2,J,K)+ARHO*U(2,J,K)
      AIM(2,J,K)=0.
      DO 403 I=2,L2
      IF(I .EQ. L2) GO TO 404
      ARHO=AX(J,K)*(FX(I+1)*RHO(I+1,J,K)+FXM(I+1)*RHO(I,J,K))
      FLOW=ARHO*U(I+1,J,K)
      CON(I,J,K)=CON(I,J,K)-FLOW
      CON(I+1,J,K)=CON(I+1,J,K)+FLOW
      AIP(I,J,K)=ARHO*DU(I+1,J,K)
      AIM(I+1,J,K)=AIP(I,J,K)
      GO TO 405
  404 ARHO=AX(J,K)*RHO(L1,J,K)
      CON(I,J,K)=CON(I,J,K)-ARHO*U(L1,J,K)
      AIP(I,J,K)=0.
  405 IF(J .EQ. M2) GO TO 406
      ARHO=AY(I,K)*(FY(J+1)*RHO(I,J+1,K)+FYM(J+1)*RHO(I,J,K))
      FLOW=ARHO*V(I,J+1,K)
      CON(I,J,K)=CON(I,J,K)-FLOW
      CON(I,J+1,K)=CON(I,J+1,K)+FLOW
      AJP(I,J,K)=ARHO*DV(I,J+1,K)
      AJM(I,J+1,K)=AJP(I,J,K)
      GO TO 407
  406 ARHO=AY(I,K)*RHO(I,M1,K)
      CON(I,J,K)=CON(I,J,K)-ARHO*V(I,M1,K)
      AJP(I,J,K)=0.
  407 IF(K .EQ. N2) GO TO 408
      ARHO=AZ(I,J)*(FZ(K+1)*RHO(I,J,K+1)+FZM(K+1)*RHO(I,J,K))