C***************************************************************
C  THIS PROGRAM IS FOR TURBULENT FLOW OF SUDDEN EXPANSION.     *
C  THIS WAS ORIGINALLY USED AS 'TEACH' PROGRAM IN KARLSRUHE    *
C  UNIV., AND PROF.RODI PERMITTED ARAKAWA TO ADAPT IT FOR      *
C  ACADEMIC USE. ARAKAWA MODIFIED AND ARRANGED IT FOR LECTURES *
C  IN UNIVERSITY OF TOKYO AND THIS BOOK.                       *
C***************************************************************
C
      PROGRAM MAIN
      PARAMETER (NX=22, NY=22)
      COMMON
     &/UVEL/RESORU,NSWPU,URFU,DXEPU(NX),DXPWU(NX),SEWU(NX)
     &/VVEL/RESORV,NSWPV,URFV,DYNPV(NY),DYPSV(NY),SNSV(NY)
     &/PCOR/RESORM,NSWPP,URFP,DU(NX,NY),DV(NX,NY),IPREF,JPREF
     &/TEN/RESORK,NSWPK,URFK
     &/TDIS/RESORE,NSWPD,URFE
     &/VAR/U(NX,NY),V(NX,NY),P(NX,NY),PP(NX,NY),
     &     TE(NX,NY),ED(NX,NY)
     &/ALL/NI,NJ,NIM1,NJM1,GREAT
     &/GEOM/X(NX),Y(NY),DXEP(NX),DXPW(NX),DYNP(NY),DYPS(NY),
     &      SNS(NY),SEW(NX),XU(NX),YV(NY)
     &/FLUPR/URFVIS,VISCOS,DENSIT,PRANDT,VIS(NX,NY)
     &/KASE T1/UIN,TEIN,EDIN,FLOWIN,ALAMDA,YSMALL,YLARGE,
     &       AL1,AL2,JSTEP,JSTP1,JSTM1
     &/TURB/GEN(NX,NY),CD,CMU,C1,C2,CAPPA,ELOG,PRED,PRTE,PRM
     &/WALLF/YPLUSN(NX),XPLUSW(NY),TAUN(NX),TAUW(NY)
     &/COEF/AP(NX,NY),AN(NX,NY),AS(NX,NY),AE(NX,NY),AW(NX,NY),
     &      SU(NX,NY),SP(NX,NY)
      CHARACTER*12 HEDU,HEDV,HEDP,HEDK,HEDD,HEDM
      DATA HEDU/' U-VELOCITY '/
      DATA HEDV/' V-VELOCITY '/
      DATA HEDP/' PRESSURE   '/
      DATA HEDK/'TURB.-ENERGY'/
      DATA HEDD/' DISSIPATION'/
      DATA HEDM/' VISCOSITY  '/
C
CHAPTER  1  1  1  PARAMETERS AND CONTROL INDICES  1  1  1  1  1  1
        GREAT = 1.E30
        NITER = 0
        NSWPU = 2
        NSWPV = 2
        NSWPP = 3
        NSWPK = 2
        NSWPD = 2
C-----GRID
        NI  = 16
        NJ  = 16
        NIM1= NI-1
        NJM1= NJ-1
        CALL GRID
C-----FLUID PROPERTIES
        DENSIT = 1000.0
        VISCOS = 0.8E-3
C-----TURBULENCE CONSTANTS
        CMU  = 0.09
        CD   = 1.00
        C1   = 1.44
        C2   = 1.92
        CAPPA=.4187
        ELOG = 9.793
        PRED = 1.3
        PRTE = 1.0
C-----BOUNDARY VALUES
        UIN    = 1.0
        TURBIN = 0.03
        TEIN   = TURBIN*UIN**2
        ALAMDA = 0.005
        EDIN   = TEIN**1.5/(ALAMDA*YLARGE)
C-----PRESSURE CALCULATION
        IPREF = 2
        JPREF = 2
C-----PROGRAM CONTROL AND MONITOR
        MAXIT= 500
        IMON = 10
        JMON = 10
        URFU = 0.5
        URFV = 0.5
        URFP = 0.8
        URFE = 0.7
        URFK = 0.7
        URFVIS = 0.7
        SORMAX = 1.0E-3
C
CHAPTER  2  2  2  INITIAL OPERATIONS  2  2  2  2  2  2  2  2  2
C-----CALCULATE GEOMETRICAL QUANTITIES AND SET VARIABLES TO ZERO
        CALL INIT
C-----INITIALISE VARIABLE FIELDS
        DO  200  J=2,JSTEP
            U(2,J)  = UIN
            TE(1,J) = TEIN
            ED(1,J) = EDIN
  200 CONTINUE
        DO  202  I=2,NI
            FACTOR = YV(JSTP1)/YV(NJ)
            DO  202  J=2,NJM1
            TE(I,J) = TEIN
            ED(I,J) = EDIN
  202 IF(I.NE.2) U(I,J) = UIN*FACTOR
        DO  204  I=2,NIM1
  204 YPLUSN(I)=11.0
        DO  206  J=JSTEP,NJ
            XPLUSW(J) = 11.0
  206 IF(J.EQ.JSTEP) XPLUSW(J) = 0.0
      CALL PROPS
C-----INITIAL OUTPUT
        WRITE(6,210)
  210 FORMAT(1H ,5X,'LAMINAR FLOW IN SUDDEN EXPANSION CHANNEL')
        WRITE(6,220) UIN
  220 FORMAT(//1H 5X,'INLET FLUID VELOCITY   =   ',E11.3/)
        YSDRL  = YSMALL/YLARGE
        RE     = UIN*DENSIT*2.0*YSMALL/VISCOS
        WRITE(6,230) RE
  230 FORMAT(1H ,5X,  'REYNOLDS NUMBER       =   ',E11.3/)
        WRITE(6,240) YSDRL
  240 FORMAT(1H ,5X,  'EXPANSION RATIO       =   ',E11.3/)
        WRITE(6,250) DENSIT
  250 FORMAT(1H ,5X,   'FLUID DENSITY        =   ',E11.3/)
        WRITE(6,260) VISCOS
  260 FORMAT(1H ,5X,   'LAMINAR VISCOSITY    =   ',E11.3/)
C-----CALCULATE RESIDUAL SOURCES NORMALIZATION FACTORS
        FLOWIN = 0.0
        XMONIN = 0.0
        TEINS  = 0.0
        EDINS  = 0.0
        DO 270 J=2,NJM1
            ARDEN  = DENSIT*SNS(J)
            FLOWIN = FLOWIN+ARDEN*U(2,J)
            XMONIN = XMONIN+ARDEN*U(2,J)*U(2,J)
            TEINS  = TEINS +ARDEN*U(2,J)*U(2,J)*U(2,J)
            EDINS  = EDINS +ARDEN*U(2,J)**4/ALAMDA/YLARGE
  270 CONTINUE
        WRITE(6,280) IMON,JMON
  280 FORMAT(1H0,'ITER   ',
     &  'I---------------ABSOLUTE RESIDUAL SOURCE SUMS-------',
     &  '-------I '/1H ,7X,'I-------FIELD VALUES AT MONITORING',
     &  ' LOCATION  ', '(',I2,',',I2,') ------I' / 2X,'NO.',6X,
     &  'UMOM',8X,'VMOM',8X,'MAS',9X,'TKIN',8X,'DISP'/1H ,10X,
     &  'U',11X,'V',11X,'P',11X,'K',11X,'D'/)
C
CHAPTER  3  3  3  3  ITERATION LOOP  3  3  3  3  3  3  3  3  3
  300 NITER=NITER+1
        CALL CALCU
        CALL CALCV
        CALL CALCP
        CALL CALCTE
        CALL CALCED
C-----UPDATA FLUID PROPERTIES
        CALL PROPS
C-----INTERMEDIATE OUTPUT
        RESORM = RESORM/FLOWIN
        RESORU = RESORU/XMONIN
        RESORV = RESORV/XMONIN
        RESORK = RESORK/TEINS
        RESORE = RESORE/EDINS
        WRITE(6,310) NITER,RESORU,RESORV,RESORM,RESORK,RESORE,
     &                     U(IMON,JMON),V(IMON,JMON),P(IMON,JMON),
     &                     TE(IMON,NJM1),ED(IMON,NJM1)
  310 FORMAT(1H ,I3,2X,5E12.3/1H ,5X,5E12.3)
C-----TERMINATION TESTS
        SORCE  = AMAX1(RESORM,RESORU,RESORV)
        IF((NITER.EQ.20.AND.SORCE.GT.1.0E4*SORMAX)
     &                                 .OR.(NITER.EQ.MAXIT)) THEN
              WRITE(6,320)
  320       FORMAT(1H ,'SOLUTION DIVERGES OR NEEDS MORE ITERATIONS')
            GO TO 400
        ENDIF
        IF(SORCE.GT.SORMAX) GO TO 300
C
CHAPTER  4  4  4  FINAL OPERATIONS AND OUTPUT  4  4  4  4  4  4
  400 CONTINUE
        CALL PRINT(2,2,NI,NJ,XU,Y,U,HEDU)
        CALL PRINT(2,2,NI,NJ,X,YV,V,HEDV)
        CALL PRINT(2,2,NI,NJ,X,Y,P,HEDP)
        CALL PRINT(2,2,NI,NJ,X,Y,TE,HEDK)
        CALL PRINT(2,2,NI,NJ,X,Y,ED,HEDD)
C--OUTPUTS FOR GRAPHIC PROGRAMS
        OPEN (10,FILE='GRAPH.DAT',STATUS='UNKNOWN')
        WRITE(10,410)NI
        WRITE(10,410)NJ
        WRITE(10,420)X(NI)
        WRITE(10,420)Y(NJ)
        WRITE(10,430)
        DO 440 J=1,NJ
            DO 440 I=1,NI
              WRITE(10,420)X(I),Y(J),U(I,J),V(I,J),P(I,J),
     &                           TE(I,J),ED(I,J)
  440 CONTINUE
        CLOSE (10)
      
	OPEN (10,FILE='VECTOR.DAT',STATUS='UNKNOWN')
      DO 445 J=1,NJ
            DO 445 I=1,NI
              WRITE(10,420)X(I),Y(J),U(I,J)*0.1,V(I,J)*0.1
  445 CONTINUE
        CLOSE (10)
        STOP
  410 FORMAT(1H ,I5)
  420 FORMAT(1H ,7E11.3)
  430 FORMAT(1H ,1H )
        END
C
        SUBROUTINE GRID
        PARAMETER (NX=22, NY=22)
      COMMON
     &/ALL/NI,NJ,NIM1,NJM1,GREAT
     &/GEOM/X(NX),Y(NY),DXEP(NX),DYPW(NX),DYNP(NY),DYPS(NY),
     &            SNS(NY),SEW(NX),XU(NX),YV(NY)
     &/KASE T1/UIN,TEIN,EDIN,FLOWIN,ALAMDA,YSMALL,YLARGE,
     &             AL1,AL2,JSTEP,JSTP1,JSTM1
        JSTEP  = 8
        JSTP1  = JSTEP+1
        JSTM1  = JSTEP-1
        YLARGE = 0.1
        ALTOT  = 10.0*YLARGE
        EPSX   = 1.15
        SUMX   = 0.5*EPSX**(NI-4)+(EPSX**(NI-3)-1.)/(EPSX-1.)+0.5
        DX     = ALTOT/SUMX
        X(1)   =-0.5*DX
        X(2)   =-X(1)
        DO 100 I=3,NIM1
          X(I) = X(I-1)+DX
            DX   = EPSX*DX
  100 CONTINUE
        X(NI) = X(NIM1)-X(NI-2)+X(NIM1)
        DY    = YLARGE/FLOAT(NJ-2)
        Y(1)  =-0.5*DY
        DO  110  J=2,NJ
            Y(J) = Y(J-1)+DY
  110 CONTINUE
        YSMALL = 0.5*(Y(JSTEP)+Y(JSTP1))
        RETURN
        END
C
      SUBROUTINE INIT
      PARAMETER (NX=22, NY=22)
      COMMON
     &/UVEL/RESORU,NSWPU,URFU,DXEPU(NX),DXPWU(NX),SEWU(NX)
     &/VVEL/RESORV,NSWPV,URFV,DYNPV(NY),DYPSV(NY),SNSV(NY)
     &/PCOR/RESORM,NSWPP,URFP,DU(NX,NY),DV(NX,NY),IPREF,JPREF
     &/VAR/U(NX,NY),V(NX,NY),P(NX,NY),PP(NX,NY),
     &     TE(NX,NY),ED(NX,NY)
     &/ALL/NI,NJ,NIM1,NJM1,GREAT
     &/GEOM/X(NX),Y(NY),DXEP(NX),DXPW(NX),DYNP(NY),DYPS(NY),
     &      SNS(NY),SEW(NX),XU(NX),YV(NY)
     &/FLUPR/URFVIS,VISCOS,DENSIT,PRANDT,VIS(NX,NY)
     &/KASE T1/UIN,TEIN,EDIN,FLOWIN,ALAMDA,YSMALL,YLARGE,
     &       AL1,AL2,JSTEP,JSTP1,JSTM1
     &/TURB/GEN(NX,NY),CD,CMU,C1,C2,CAPPA,ELOG,PRED,PRTE,PRM
     &/COEF/AP(NX,NY),AN(NX,NY),AS(NX,NY),AE(NX,NY),AW(NX,NY),
     &      SU(NX,NY),SP(NX,NY)
C
CHAPTER  1  1  1  CALCULATE GEOMETRICAL QUANTITIES  1  1  1  1  1
        DXPW(1) = 0.0
        DXEP(NI)= 0.0
        DO 100 I=1,NIM1
            DXEP(I) = X(I+1)-X(I)
  100 DXPW(I+1)=DXEP(I)
        DYPS(1)  = 0.0
        DYNP(NJ) = 0.0
        DO 110 J=1,NJM1
            DYNP(J) = Y(J+1)-Y(J)
  110 DYPS(J+1) = DYNP(J)
        SEW(1) = 0.0
        SEW(NI)= 0.0
        DO 120 I=2,NIM1
  120 SEW(I) = 0.5*(DXEP(I)+DXPW(I))
        SNS(1) = 0.0
        SNS(NJ)= 0.0
        DO 130 J=2,NJM1
  130 SNS(J) = 0.5*(DYNP(J)+DYPS(J))
        XU(1)=0.0
        DO 140 I=2,NI
  140 XU(I) = 0.5*(X(I)+X(I-1))
        DXPWU(1) = 0.0
        DXPWU(2) = 0.0
        DXEPU(1) = 0.0
        DXEPU(NI)= 0.0
        DO 150 I=2,NIM1
            DXEPU(I) = XU(I+1)-XU(I)
  150 DXPWU(I+1) = DXEPU(I)
        SEWU(1) = 0.0
        SEWU(2) = 0.0
        DO 160 I=3,NIM1
  160 SEWU(I) = 0.5*(DXEPU(I)+DXPWU(I))
        YV(1) = 0.0
        DO 170 J=2,NJ
  170 YV(J) = 0.5*(Y(J)+Y(J-1))
        DYPSV(1) = 0.0
        DYPSV(2) = 0.0
        DYNPV(NJ)= 0.0
        DO 180 J=2,NJM1
            DYNPV(J) = YV(J+1)-YV(J)
  180 DYPSV(J+1) = DYNPV(J)
        SNSV(1) = 0.0
        SNSV(2) = 0.0
        SNSV(NJ)= 0.0
        DO 190 J=3,NJM1
  190 SNSV(J) = 0.5*(DYNPV(J)+DYPSV(J))
C
CHAPTER  2  2  2  2  SET VARIABLES TO ZERO  2  2  2  2  2  2
      DO 200 I=1,NI
      DO 200 J=1,NJ
          U(I,J) = 0.0
          V(I,J) = 0.0
          P(I,J) = 0.0
          PP(I,J)= 0.0
          TE(I,J)= 0.0
          ED(I,J)= 0.0
          VIS(I,J)= VISCOS
          DU(I,J) = 0.0
          DV(I,J) = 0.0
          SU(I,J) = 0.0
          SP(I,J) = 0.0
 200  CONTINUE
        RETURN
        END
C
      SUBROUTINE CALCU
      PARAMETER (NX=22, NY=22)
      COMMON
     &/UVEL/RESORU,NSWPU,URFU,DXEPU(NX),DXPWU(NX),SEWU(NX)
     &/PCOR/RESORM,NSWPP,URFP,DU(NX,NY),DV(NX,NY),IPREF,JPREF
     &/VAR/U(NX,NY),V(NX,NY),P(NX,NY),PP(NX,NY),
     &     TE(NX,NY),ED(NX,NY)
     &/ALL/NI,NJ,NIM1,NJM1,GREAT
     &/GEOM/X(NX),Y(NY),DXEP(NX),DXPW(NX),DYNP(NY),DYPS(NY),
     &      SNS(NY),SEW(NX),XU(NX),YV(NY)
     &/FLUPR/URFVIS,VISCOS,DENSIT,PRANDT,VIS(NX,NY)
     &/COEF/AP(NX,NY),AN(NX,NY),AS(NX,NY),AE(NX,NY),AW(NX,NY),
     &      SU(NX,NY),SP(NX,NY)
     &/KASE T1/UIN,TEIN,EDIN,FLOWIN,ALAMDA,YSMALL,YLARGE,
     &       AL1,AL2,JSTEP,JSTP1,JSTM1
     &/TURB/GEN(NX,NY),CD,CMU,C1,C2,CAPPA,ELOG,PRED,PRTE,PRM
     &/WALLF/YPLUSN(NX),XPLUSW(NY),TAUN(NX),TAUW(NY)
C
CHAPTER  1  1  1  1  ASSEMBLY OF COEFFICIENTS  1  1  1  1  1  1  1
        DO 100 I=3,NIM1
        DO 100 J=2,NJM1
C-----CALCULATE CONVECTION COEFFICIENTS
            GN  = DENSIT*V(I,J+1)
            GNW = DENSIT*V(I-1,J+1)
            GS  = DENSIT*V(I,J)
            GSW = DENSIT*V(I-1,J)
            GE  = DENSIT*U(I+1,J)
            GP  = DENSIT*U(I,J)
            GW  = DENSIT*U(I-1,J)
            CN  = 0.5*(GN+GNW)*SEWU(I)
            CS  = 0.5*(GS+GSW)*SEWU(I)
            CE  = 0.5*(GE+GP)*SNS(J)
            CW  = 0.5*(GP+GW)*SNS(J)
C-----CALCULATE DIFFUSION COEFFICIENTS
            VISN=0.25*(VIS(I,J)+VIS(I,J+1)+VIS(I-1,J)+VIS(I-1,J+1))
            VISS=0.25*(VIS(I,J)+VIS(I,J-1)+VIS(I-1,J)+VIS(I-1,J-1))
            DN  = VISN*SEWU(I)/DYNP(J)
            DS  = VISS*SEWU(I)/DYPS(J)
            DE  = VIS(I,J)*SNS(J)/DXEPU(I)
            DW  = VIS(I-1,J)*SNS(J)/DXPWU(I)
C-----ASSEMBLE MAIN COEFFICIENTS
            AN(I,J)=AMAX1(ABS(0.5*CN),DN)-0.5*CN
            AS(I,J)=AMAX1(ABS(0.5*CS),DS)+0.5*CS
            AE(I,J)=AMAX1(ABS(0.5*CE),DE)-0.5*CE
            AW(I,J)=AMAX1(ABS(0.5*CW),DW)+0.5*CW
            DU(I,J)=SNS(J)
            SU(I,J)=DU(I,J)*(P(I-1,J)-P(I,J))
            SP(I,J)=0.0
            DUDXP  =(U(I+1,J)-U(I,J))/SEW(I)
            DUDXM  =(U(I,J)-U(I-1,J))/SEW(I-1)
            SU(I,J)=SU(I,J)+(VIS(I,J)*DUDXP-VIS(I-1,J)*DUDXM)
     &                              		  /SEWU(I)*SEWU(I)*SNS(J)
            GAMP =0.25*(VIS(I,J)+VIS(I-1,J)+VIS(I,J+1)+VIS(I-1,J+1))
            DVDXP=(V(I,J+1)-V(I-1,J+1))/DXPW(I)
            GAMM =0.25*(VIS(I,J)+VIS(I-1,J)+VIS(I,J-1)+VIS(I-1,J-1))
            DVDXM=(V(I,J)-V(I-1,J))/DXPW(I)
            SU(I,J) =SU(I,J)+(GAMP*DVDXP-GAMM*DVDXM)
     &                                /SNS(J)*SEWU(I)*SNS(J)
  100 CONTINUE
C
CHAPTER  2  2  2  2  PROBLEM MODIFICATIONS  2  2  2  2  2  2  2
C-----TOP WALL
        CDTERM = CMU**0.25
        YP = YV(NJ)-Y(NJM1)
        J = NJM1
        DO 200  I=3,NIM1
          SQRTK = SQRT(0.5*(TE(I,J)+TE(I-1,J)))
            YPLUSA= 0.5*(YPLUSN(I)+YPLUSN(I-1))
            IF(YPLUSA.LE.11.63)          GO  TO  210
            TMULT = DENSIT*CDTERM*SQRTK*CAPPA/ALOG(ELOG*YPLUSA)
            GO  TO  220
  210      TMULT  = VISCOS/YP
  220      TAUN(I)=-TMULT*U(I,J)
            SP(I,J)= SP(I,J)-TMULT*SEWU(I)
  200 AN(I,J) = 0.0
          TAUN(2) = TAUN(3)
          TAUN(NI)= TAUN(NIM1)
C-----SIDE WALL
        DO  230  J=JSTP1,NJM1
  230 AW(3,J)=0.0
C-----SYMMETRY  AXIS
        DO  240  I=1,NI
            U(I,1) = U(I,2)
            AS(I,2) = 0.0
  240 CONTINUE
C-----OUTLET
        ARDENT = 0.0
        FLOW = 0.0
        DO  250  J=2,NJM1
            ARDEN  = DENSIT*SNS(J)
            ARDENT = ARDENT+ARDEN
            FLOW   = FLOW+ARDEN*U(NIM1,J)
  250 CONTINUE
        UINC = (FLOWIN-FLOW)/ARDENT
        DO  260  J=2,NJM1
  260 U(NI,J) = U(NIM1,J)+UINC
C
CHAPTER  3  FINAL COEFF. ASSEMBLY & RESIDUAL SOURCE CALCULATION
        RESORU=0.0
        DO 300 I=3,NIM1
        DO 300 J=2,NJM1
            AP(I,J) = AN(I,J)+AS(I,J)+AE(I,J)+AW(I,J)-SP(I,J)
            DU(I,J) = DU(I,J)/AP(I,J)
            RESOR=AN(I,J)*U(I,J+1)+AS(I,J)*U(I,J-1)+AE(I,J)*U(I+1,J)
     &             +AW(I,J)*U(I-1,J)-AP(I,J)*U(I,J)+SU(I,J)
            SORVOL=GREAT*SEWU(I)*SNS(J)
            IF(-SP(I,J).GT.0.5*SORVOL)  RESOR=RESOR/SORVOL
            RESORU = RESORU+ABS(RESOR)
C-------UNDER-RELAXATION
            AP(I,J) = AP(I,J)/URFU
            SU(I,J) = SU(I,J)+(1.-URFU)*AP(I,J)*U(I,J)
            DU(I,J) = DU(I,J)*URFU
  300 CONTINUE
C
CHAPTER  4  4  4  SOLUTION OF DIFFERENCE EQUATION  4  4  4  4
        DO 400 N=1,NSWPU
  400 CALL LISOLV(3,2,NI,NJ,U)
        RETURN
        END
C
      SUBROUTINE CALCV
      PARAMETER (NX=22, NY=22)
      COMMON
     &/VVEL/RESORV,NSWPV,URFV,DYNPV(NY),DYPSV(NY),SNSV(NY)
     &/PCOR/RESORM,NSWPP,URFP,DU(NX,NY),DV(NX,NY),IPREF,JPREF
     &/VAR/U(NX,NY),V(NX,NY),P(NX,NY),PP(NX,NY),
     &     TE(NX,NY),ED(NX,NY)
     &/ALL/NI,NJ,NIM1,NJM1,GREAT
     &/GEOM/X(NX),Y(NY),DXEP(NX),DXPW(NX),DYNP(NY),DYPS(NY),
     &      SNS(NY),SEW(NX),XU(NX),YV(NY)
     &/FLUPR/URFVIS,VISCOS,DENSIT,PRANDT,VIS(NX,NY)
     &/COEF/AP(NX,NY),AN(NX,NY),AS(NX,NY),AE(NX,NY),AW(NX,NY),
     &      SU(NX,NY),SP(NX,NY)
     &/KASE T1/UIN,TEIN,EDIN,FLOWIN,ALAMDA,YSMALL,YLARGE,
     &       AL1,AL2,JSTEP,JSTP1,JSTM1
     &/TURB/GEN(NX,NY),CD,CMU,C1,C2,CAPPA,ELOG,PRED,PRTE,PRM
     &/WALLF/YPLUSN(NX),XPLUSW(NY),TAUN(NX),TAUW(NY)
C
CHAPTER  1  1  1  1  ASSEMBLY OF COEFFICIENTS  1  1  1  1  1  1  1
        DO 100 I=2,NIM1
        DO 100 J=3,NJM1
C-----CALCULATE CONVECTION COEFFICIENTS
            GN  = DENSIT*V(I,J+1)
            GP  = DENSIT*V(I,J)
            GS  = DENSIT*V(I,J-1)
            GE  = DENSIT*U(I+1,J)
            GSE = DENSIT*U(I+1,J-1)
            GW  = DENSIT*U(I,J)
            GSW = DENSIT*U(I,J-1)
            CN  = 0.5*(GN+GP)*SEW(I)
            CS  = 0.5*(GP+GS)*SEW(I)
            CE  = 0.5*(GE+GSE)*SNSV(J)
            CW  = 0.5*(GW+GSW)*SNSV(J)
C-----CALCULATE DIFFUSION COEFFICIENTS
            VISE=0.25*(VIS(I,J)+VIS(I+1,J)+VIS(I,J-1)+VIS(I+1,J-1))
            VISW=0.25*(VIS(I,J)+VIS(I-1,J)+VIS(I,J-1)+VIS(I-1,J-1))
            DN  = VIS(I,J)*SEW(I)/DYNPV(J)
            DS  = VIS(I,J-1)*SEW(I)/DYPSV(J)
            DE  = VISE*SNSV(J)/DXEP(I)
            DW  = VISE*SNSV(J)/DXPW(I)
C-----ASSEMBLE MAIN COEFFICIENTS
            AN(I,J) = AMAX1(ABS(0.5*CN),DN)-0.5*CN
            AS(I,J) = AMAX1(ABS(0.5*CS),DS)+0.5*CS
            AE(I,J) = AMAX1(ABS(0.5*CE),DE)-0.5*CE
            AW(I,J) = AMAX1(ABS(0.5*CW),DW)+0.5*CW
            DV(I,J) = SEW(I)
            SU(I,J) = DV(I,J)*(P(I,J-1)-P(I,J))
            SP(I,J) = 0.0
            DUDYP  =(U(I+1,J)-U(I+1,J-1))/DYPS(J)
            GAMP =0.25*(VIS(I,J)+VIS(I+1,J)+VIS(I,J-1)+VIS(I+1,J-1))
            GAMM =0.25*(VIS(I,J)+VIS(I-1,J)+VIS(I,J-1)+VIS(I-1,J-1))
            DUDYM  = (U(I,J)-U(I,J-1))/DYPS(J)
            SU(I,J)= SU(I,J)+(GAMP*DUDYP-GAMM*DUDYM)
     &                                 /SEW(I)*SEW(I)*SNSV(J)
            DVDYP  = (V(I,J+1)-V(I,J))/SNS(J)
            RGAMP  = VIS(I,J)
            DVDYM  = (V(I,J)-V(I,J-1))/SNS(J-1)
            RGAMM  = VIS(I,J-1)
            SU(I,J) = SU(I,J)+(RGAMP*DVDYP-RGAMM*DVDYM)
     &                              	 /SNS(J)*SEW(I)*SNSV(J)
  100 CONTINUE
C
CHAPTER  2  2  2  2  2  PROBLEM MODIFICATIONS  2  2  2  2  2  2  2
C-----SIDE WALL
        CDTERM = CMU**0.25
        XP = X(2)-XU(2)
        I = 2
        DO  200  J=JSTP1,NJM1
          SQRTK = SQRT(0.5*(TE(I,J)+TE(I,J-1)))
            XPLUSA= 0.5*(XPLUSW(J)+XPLUSW(J-1))
            IF(XPLUSA.LE.11.63)          GO  TO  210
            TMULT = DENSIT*CDTERM*SQRTK*CAPPA/ALOG(ELOG*XPLUSA)
            GO  TO  220
  210      TMULT = VISCOS/XP
  220      TAUW(J) = -TMULT*V(I,J)
            SP(I,J) = SP(I,J)-TMULT*SNSV(J)
            AW(I,J) = 0.0
  200 CONTINUE
        TAUW(JSTEP) = TAUW(JSTP1)
          TAUW(NJ) = TAUW(NJM1)
C-----TOP WALL
        DO  230  I=2,NIM1
  230 AN(I,NJM1) = 0.0
C-----SYMMETRY  AXIS
        DO  240  I=2,NIM1
  240 AS(I,3) = 0.0
C
CHAPTER  3  FINAL COEFF. ASSEMBLY & RESIDUAL SOURCE CALCULATION
        RESORV = 0.0
        DO 300 I=2,NIM1
        DO 300 J=3,NJM1
            AP(I,J) = AN(I,J)+AS(I,J)+AE(I,J)+AW(I,J)-SP(I,J)
            DV(I,J) = DV(I,J)/AP(I,J)
            RESOR=AN(I,J)*V(I,J+1)+AS(I,J)*V(I,J-1)+AE(I,J)*V(I+1,J)
     &             +AW(I,J)*V(I-1,J)-AP(I,J)*V(I,J)+SU(I,J)
            SORVOL = GREAT*SEW(I)*SNSV(J)
            IF(-SP(I,J).GT.0.5*SORVOL)  RESOR=RESOR/SORVOL
            RESORV = RESORV+ABS(RESOR)
C-----UNDER-RELAXATION
            AP(I,J) = AP(I,J)/URFV
            SU(I,J) = SU(I,J)+(1.-URFV)*AP(I,J)*V(I,J)
            DV(I,J) = DV(I,J)*URFV
  300 CONTINUE
C
CHAPTER  4  4  4  SOLUTION OF DIFFERENCE EQUATION  4  4  4  4  4
        DO 400 N=1,NSWPV
  400 CALL LISOLV(2,3,NI,NJ,V)
        RETURN
        END
C
      SUBROUTINE CALCP
      PARAMETER (NX=22, NY=22)
      COMMON
     &/PCOR/RESORM,NSWPP,URFP,DU(NX,NY),DV(NX,NY),IPREF,JPREF
     &/VAR/U(NX,NY),V(NX,NY),P(NX,NY),PP(NX,NY),
     &     TE(NX,NY),ED(NX,NY)
     &/ALL/NI,NJ,NIM1,NJM1,GREAT
     &/GEOM/X(NX),Y(NY),DXEP(NX),DXPW(NX),DYNP(NY),DYPS(NY),
     &      SNS(NY),SEW(NX),XU(NX),YV(NY)
     &/FLUPR/URFVIS,VISCOS,DENSIT,PRANDT,VIS(NX,NY)
     &/COEF/AP(NX,NY),AN(NX,NY),AS(NX,NY),AE(NX,NY),AW(NX,NY),
     &      SU(NX,NY),SP(NX,NY)
        RESORM=0.0
C
CHAPTER  1  1  1  ASSEMBLY OF COEFFICIENTS  1  1  1  1  1  1  1
        DO 100 I=2,NIM1
        DO 100 J=2,NJM1
C-----CALCULATE COEFFICIENTS
            AN(I,J) = DENSIT*SEW(I)*DV(I,J+1)
            AS(I,J) = DENSIT*SEW(I)*DV(I,J)
            AE(I,J) = DENSIT*SNS(J)*DU(I+1,J)
            AW(I,J) = DENSIT*SNS(J)*DU(I,J)
C-----CALCULATE SOURCE TERMS
            CN = DENSIT*V(I,J+1)*SEW(I)
            CS = DENSIT*V(I,J)*SEW(I)
            CE = DENSIT*U(I+1,J)*SNS(J)
            CW = DENSIT*U(I,J)*SNS(J)
            SMP= CN-CS+CE-CW
            SP(I,J) = 0.0
            SU(I,J) =-SMP
C-----COMPUTE SUM OF ABSOLUTE MASS SOURCES
          RESORM = RESORM+ABS(SMP)
  100 CONTINUE
C
CHAPTER  3  3  3  3  FINAL COEFFICIENT ASSEMBLY   3  3  3  3
        DO 300 I=2,NIM1
        DO 300 J=2,NJM1
          AP(I,J)=AN(I,J)+AS(I,J)+AE(I,J)+AW(I,J)-SP(I,J)
  300 CONTINUE
C
CHAPTER  4  4  4  4  SOLUTION OF DIFFERENCE EQUATIONS  4  4  4
        DO 400 N=1,NSWPP
  400 CALL LISOLV(2,2,NI,NJ,PP)
C
CHAPTER  5  5  5  CORRECT VELOCITIES AND PRESSURE  5  5  5  5
C-----VELOCITIES
        DO 500 I=2,NIM1
        DO 500 J=2,NJM1
            IF(I.NE.2) U(I,J)=U(I,J)+DU(I,J)*(PP(I-1,J)-PP(I,J))
            IF(J.NE.2) V(I,J)=V(I,J)+DV(I,J)*(PP(I,J-1)-PP(I,J))
  500 CONTINUE
C-----PRESSURES (WITH PROVISION FOR UNDER-RELAXATION)
        PPREF=PP(IPREF,JPREF)
        DO 510 I=2,NIM1
        DO 510 J=2,NJM1
            P(I,J) = P(I,J)+URFP*(PP(I,J)-PPREF)
            PP(I,J)= 0.0
  510 CONTINUE
        RETURN
        END
C
      SUBROUTINE CALCTE
      PARAMETER (NX=22, NY=22)
      COMMON
     &/TEN/RESORK,NSWPK,URFK
     &/VAR/U(NX,NY),V(NX,NY),P(NX,NY),PP(NX,NY),
     &     TE(NX,NY),ED(NX,NY)
     &/ALL/NI,NJ,NIM1,NJM1,GREAT
     &/GEOM/X(NX),Y(NY),DXEP(NX),DXPW(NX),DYNP(NY),DYPS(NY),
     &      SNS(NY),SEW(NX),XU(NX),YV(NY)
     &/FLUPR/URFVIS,VISCOS,DENSIT,PRANDT,VIS(NX,NY)
     &/COEF/AP(NX,NY),AN(NX,NY),AS(NX,NY),AE(NX,NY),AW(NX,NY),
     &      SU(NX,NY),SP(NX,NY)
     &/TURB/GEN(NX,NY),CD,CMU,C1,C2,CAPPA,ELOG,PRED,PRTE,PRM
     &/WALLF/YPLUSN(NX),XPLUSW(NY),TAUN(NX),TAUW(NY)
     &/KASE T1/UIN,TEIN,EDIN,FLOWIN,ALAMDA,YSMALL,YLARGE,
     &       AL1,AL2,JSTEP,JSTP1,JSTM1
C
CHAPTER  1  1  1  1  ASSEMBLY OF COEFFICIENTS  1  1  1  1  1  1
        DO 100 I=2,NIM1
        DO 100 J=2,NJM1
C-----CALCULATE CONVECTION COEFFICIENTS
            GN  = DENSIT*V(I,J+1)
            GS  = DENSIT*V(I,J)
            GE  = DENSIT*U(I+1,J)
            GW  = DENSIT*U(I,J)
            CN  = GN*SEW(I)
            CS  = GS*SEW(I)
            CE  = GE*SNS(J)
            CW  = GW*SNS(J)
C-----CALCULATE DIFFUSION COEFFICIENTS
            GAMN = 0.5*(VIS(I,J)+VIS(I,J+1))/PRTE
            GAMS = 0.5*(VIS(I,J)+VIS(I,J-1))/PRTE
            GAME = 0.5*(VIS(I,J)+VIS(I+1,J))/PRTE
            GAMW = 0.5*(VIS(I,J)+VIS(I-1,J))/PRTE
            DN  = GAMN*SEW(I)/DYNP(J)
            DS  = GAMS*SEW(I)/DYPS(J)
            DE  = GAME*SNS(J)/DXEP(I)
            DW  = GAMW*SNS(J)/DXPW(I)
C-----SOURCE TERMS
            DUDX=(U(I+1,J)-U(I,J))/SEW(I)
            DVDY=(V(I,J+1)-V(I,J))/SNS(J)
            DUDY=((U(I,J)+U(I+1,J)+U(I,J+1)+U(I+1,J+1))/4.
     &            -(U(I,J)+U(I+1,J)+U(I,J-1)+U(I+1,J-1))/4.)/SNS(J)
            DVDX=((V(I,J)+V(I,J+1)+V(I+1,J)+V(I+1,J+1))/4.
     &            -(V(I,J)+V(I,J+1)+V(I-1,J)+V(I-1,J+1))/4.)/SEW(I)
            GEN(I,J)=(2.*(DUDX**2+DVDY**2)+(DUDY+DVDX)**2)*VIS(I,J)
C-----ASSEMBLE MAIN COEFFICIENTS -
            AN(I,J) = AMAX1(ABS(0.5*CN),DN)-0.5*CN
            AS(I,J) = AMAX1(ABS(0.5*CS),DS)+0.5*CS
            AE(I,J) = AMAX1(ABS(0.5*CE),DE)-0.5*CE
            AW(I,J) = AMAX1(ABS(0.5*CW),DW)+0.5*CW
            SU(I,J) = GEN(I,J)*SNS(J)*SEW(I)
            SP(I,J) =-CD*CMU*DENSIT**2*TE(I,J)
     &                     *SNS(J)*SEW(I)/VIS(I,J)
  100 CONTINUE
C
CHAPTER  2  2  2  2  PROBLEM MODIFICATIONS  2  2  2  2  2  2
C-----TOP WALL
        CDTERM = CMU**0.25
        YP = YV(NJ)-Y(NJM1)
        J  = NJM1
          DO  210  I=2,NIM1
                YPLUSN(I) = DENSIT*SQRT(TE(I,J))*CDTERM*YP/VISCOS
            TERM = ABS(TAUN(I+1)+TAUN(I))/2.
            TERM = TERM/(DENSIT*SQRT(CMU))
            SU(I,J) = GREAT*TERM
            SP(I,J) =-GREAT
  210      CONTINUE
C-----SIDE WALL
        XP = X(2)-XU(2)
        I = 2
        DO  220  J=JSTP1,NJM1
            XPLUSW(J) = DENSIT*SQRT(TE(I,J))*CDTERM*XP/VISCOS
            TERM = ABS(TAUW(J+1)+TAUW(J))/2.
            TERM = TERM/(DENSIT*SQRT(CMU))
            SU(I,J) = GREAT*TERM
            SP(I,J) =-GREAT
  220      CONTINUE
C-----SYMMETRY  AXIS
        DO  230  I=2,NIM1
  230 AS(I,2) = 0.0
C
CHAPTER  3  FINAL COEFF. ASSEMBLY & RESIDUAL SOURCE CALCULATION
        RESORK = 0.0
        DO 300 I=2,NIM1
        DO 300 J=2,NJM1
            AP(I,J) = AN(I,J)+AS(I,J)+AE(I,J)+AW(I,J)-SP(I,J)
            RESOR   = AN(I,J)*TE(I,J+1)+AS(I,J)*TE(I,J-1)+AE(I,J)
     &       *TE(I+1,J)+AW(I,J)*TE(I-1,J)-AP(I,J)*TE(I,J)+SU(I,J)
            SORVOL  = GREAT*SEW(I)*SNS(J)
            IF(-SP(I,J).GT.0.5*SORVOL)  RESOR=RESOR/SORVOL
            RESORK  = RESORK+ABS(RESOR)
C-----UNDER-RELAXATION
            AP(I,J) = AP(I,J)/URFK
            SU(I,J) = SU(I,J)+(1.-URFK)*AP(I,J)*TE(I,J)
  300 CONTINUE
C
CHAPTER  4  4  4  SOLUTION OF DIFFERENCE EQUATIONS  4  4  4  4
        DO 400 N=1,NSWPK
  400 CALL LISOLV(2,2,NI,NJ,TE)
        RETURN
        END
C
      SUBROUTINE CALCED
      PARAMETER (NX=22, NY=22)
      COMMON
     &/TDIS/RESORE,NSWPD,URFE
     &/ALL/NI,NJ,NIM1,NJM1,GREAT
     &/GEOM/X(NX),Y(NY),DXEP(NX),DXPW(NX),DYNP(NY),DYPS(NY),
     &      SNS(NY),SEW(NX),XU(NX),YV(NY)
     &/FLUPR/URFVIS,VISCOS,DENSIT,PRANDT,VIS(NX,NY)
     &/COEF/AP(NX,NY),AN(NX,NY),AS(NX,NY),AE(NX,NY),AW(NX,NY),
     &      SU(NX,NY),SP(NX,NY)
     &/TURB/GEN(NX,NY),CD,CMU,C1,C2,CAPPA,ELOG,PRED,PRTE,PRM
     &/WALLF/YPLUSN(NX),XPLUSW(NY),TAUN(NX),TAUW(NY)
     &/VAR/U(NX,NY),V(NX,NY),P(NX,NY),PP(NX,NY),
     &     TE(NX,NY),ED(NX,NY)
     &/KASE T1/UIN,TEIN,EDIN,FLOWIN,ALAMDA,YSMALL,YLARGE,
     &       AL1,AL2,JSTEP,JSTP1,JSTM1
C
CHAPTER  1  1  1  1  ASSEMBLY OF COEFFICIENTS  1  1  1  1  1  1
        DO 100 I=2,NIM1
        DO 100 J=2,NJM1
C-----CALCULATE CONVECTION COEFFICIENTS
            GN  = DENSIT*V(I,J+1)
            GS  = DENSIT*V(I,J)
            GE  = DENSIT*U(I+1,J)
            GW  = DENSIT*U(I,J)
            CN  = GN*SEW(I)
            CS  = GS*SEW(I)
            CE  = GE*SNS(J)
            CW  = GW*SNS(J)
C-----CALCULATE DIFFUSION COEFFICIENTS
            GAMN = 0.5*(VIS(I,J)+VIS(I,J+1))/PRED
            GAMS = 0.5*(VIS(I,J)+VIS(I,J-1))/PRED
            GAME = 0.5*(VIS(I,J)+VIS(I+1,J))/PRED
            GAMW = 0.5*(VIS(I,J)+VIS(I-1,J))/PRED
            DN  = GAMN*SEW(I)/DYNP(J)
            DS  = GAMS*SEW(I)/DYPS(J)
            DE  = GAME*SNS(J)/DXEP(I)
            DW  = GAMW*SNS(J)/DXPW(I)
C-----ASSEMBLE MAIN COEFFICIENTS -
            AN(I,J) = AMAX1(ABS(0.5*CN),DN)-0.5*CN
            AS(I,J) = AMAX1(ABS(0.5*CS),DS)+0.5*CS
            AE(I,J) = AMAX1(ABS(0.5*CE),DE)-0.5*CE
            AW(I,J) = AMAX1(ABS(0.5*CW),DW)+0.5*CW
            SU(I,J) = C1*GEN(I,J)*SNS(J)*SEW(I)*ED(I,J)/TE(I,J)
            SP(I,J) =-C2*DENSIT*ED(I,J)*SNS(J)*SEW(I)/TE(I,J)
  100 CONTINUE
C
CHAPTER  2  2  2  2  2  PROBLEM MODIFICATIONS  2  2  2  2  2  2
C-----TOP WALL
        YP = YV(NJ)-Y(NJM1)
        J = NJM1
        TERM = (CMU**.75)/(CAPPA*YP)
        DO  200  I=2,NIM1
            SU(I,J) = GREAT*TERM*TE(I,J)**1.5
            SP(I,J) =-GREAT
  200      CONTINUE
C-----SIDE WALL
        XP = X(2)-XU(2)
        I  = 2
        TERM = (CMU**.75)/(CAPPA*XP)
        NJM2 = NJ-2
        DO  210  J=JSTP1,NJM2
            SU(I,J) = GREAT*TERM*TE(I,J)**1.5
            SP(I,J) =-GREAT
  210      CONTINUE
C-----SYMMETRY  AXIS
        DO  220  I=2,NIM1
  220 AS(I,2) = 0.0
C
CHAPTER  3  FINAL COEFF. ASSEMBLY & RESIDUAL SOURCE CALCULATION
        RESORE = 0.0
        DO 300 I=2,NIM1
        DO 300 J=2,NJM1
            AP(I,J) = AN(I,J)+AS(I,J)+AE(I,J)+AW(I,J)-SP(I,J)
            RESOR   = AN(I,J)*ED(I,J+1)+AS(I,J)*ED(I,J-1)+AE(I,J)*
     &          ED(I+1,J)+AW(I,J)*ED(I-1,J)-AP(I,J)*ED(I,J)+SU(I,J)
            SORVOL  = GREAT*SNS(J)*SEW(I)
            IF(-SP(I,J).GT.0.5*SORVOL)  RESOR=RESOR/SORVOL
            RESORE  = RESORE+ABS(RESOR)
C-----UNDER-RELAXATION
            AP(I,J) = AP(I,J)/URFE
            SU(I,J) = SU(I,J)+(1.-URFE)*AP(I,J)*ED(I,J)
  300 CONTINUE
C
CHAPTER  4  4  4  SOLUTION OF DIFFERENCE EQUATIONS  4  4  4  4
        DO 400 N=1,NSWPD
  400 CALL LISOLV(2,2,NI,NJ,ED)
        RETURN
        END
C
      SUBROUTINE PROPS
      PARAMETER (NX=22, NY=22)
      COMMON
     &/FLUPR/URFVIS,VISCOS,DENSIT,PRANDT,VIS(NX,NY)
     &/VAR/U(NX,NY),V(NX,NY),P(NX,NY),PP(NX,NY),
     &     TE(NX,NY),ED(NX,NY)
     &/ALL/NI,NJ,NIM1,NJM1,GREAT
     &/TURB/GEN(NX,NY),CD,CMU,C1,C2,CAPPA,ELOG,PRED,PRTE,PRM
C-----TURBULENT VISCOSITY
        DO 100 I=2,NIM1
            VIS(I,1) = VIS(I,2)
        DO 100 J=2,NJM1
            VISOLD = VIS(I,J)
            IF(ED(I,J).EQ.0.) GO TO 120
            VIS(I,J) = DENSIT*TE(I,J)**2*CMU/ED(I,J)+VISCOS
            GO  TO  110
  120      VIS(I,J) = VISCOS
C-----UNDER-RELAX VISCOSITY
  110       VIS(I,J) = URFVIS*VIS(I,J)+(1.-URFVIS)*VISOLD
  100 CONTINUE
        RETURN
        END
C
      SUBROUTINE LISOLV(ISTART,JSTART,NI,NJ,PHI)
      PARAMETER (NX=22, NY=22)
      DIMENSION PHI(NX,NY),A(NY),B(NY),C(NY),D(NY)
      COMMON
     &/COEF/AP(NX,NY),AN(NX,NY),AS(NX,NY),AE(NX,NY),AW(NX,NY),
     &      SU(NX,NY),SP(NX,NY)
        NIM1 = NI-1
        NJM1 = NJ-1
        JSTM1 = JSTART - 1
        A(JSTM1)=0.0
C-----COMMENCE W-E SWEEP
        DO 100 I=ISTART,NIM1
            C(JSTM1) = PHI(I,JSTM1)
C-------COMMENCE S-N TRAVERSE
            DO 110 J = JSTART,NJM1
C---------ASSEMBLE TDMA COEFFICIENTS
              A(J) = AN(I,J)
              B(J) = AS(I,J)
              C(J) = AE(I,J)*PHI(I+1,J)+AW(I,J)*PHI(I-1,J)+SU(I,J)
              D(J) = AP(I,J)
C---------CALCULATE COEFFICIENTS OF RECURRENCE FORMULA
              TERM = 1./(D(J)-B(J)*A(J-1))
              A(J) = A(J)*TERM
              C(J) = (C(J)+B(J)*C(J-1))*TERM
  110       CONTINUE
C-------OBTAIN NEW PHI'S
            DO 120 JJ=JSTART,NJM1
              J = NJ+JSTM1-JJ
              PHI(I,J) = A(J)*PHI(I,J+1)+C(J)
  120      CONTINUE
  100 CONTINUE
        RETURN
        END
C
        SUBROUTINE PRINT(ISTART,JSTART,NI,NJ,X,Y,PHI,HEAD)
        PARAMETER  (NX=22, NY=22)
        DIMENSION PHI(NX,NY),X(NI),Y(NJ),STORE(NX)
        CHARACTER*12 HEAD
        ISKIP = 1
        JSKIP = 1
        WRITE(6,100)HEAD
  100 FORMAT(1H0,13(2H*-),4X,A12,4X,13(2H-*))
        ISTA = ISTART-6*ISKIP
  110 CONTINUE
          ISTA = ISTA+6*ISKIP
          IEND = ISTA+5*ISKIP
            IEND = MIN0(NI,IEND)
            WRITE(6,120)(I,I=ISTA,IEND,ISKIP)
  120       FORMAT(1H0,5H I = ,I3,5I10,7X,' Y = ')
            WRITE(6,130)
  130      FORMAT(3H  J)
            DO 140 JJ=JSTART,NJ,JSKIP
              J = JSTART+NJ-JJ
              DO 150 I=ISTA,IEND,ISKIP
                  A = PHI(I,J)
                  IF(ABS(A).LT.1.E-20) A=0.0
                  STORE(I) = A
  150        CONTINUE
              WRITE(6,160) J,(STORE(I),I=ISTA,IEND,ISKIP),Y(J)
  160        FORMAT(1H ,I2,7E10.3,E8.2)
  140      CONTINUE
            WRITE(6,170) (X(I),I=ISTA,IEND,ISKIP)
  170      FORMAT(3H0X=, 7E10.3)
        IF(IEND.LT.NI)GO TO 110
        RETURN
        END
