C############################################################
      PROGRAM COMET
C############################################################
C
C     2D, Non-orthogonal grid, laminar, unsteady, multigrid
C     =====================================================
C
C                              M. Peric, IfS, Hamburg, 1995
C############################################################
      PARAMETER (ngit=3,nicv=21,njcv=8,NXA=NICV*(2**NGIT-1)+
     *           2*NGIT,NYA=NJCV*(2**NGIT-1)+2*NGIT,
     *           NX=NICV*2**(NGIT-1)+2,NY=NJCV*2**(NGIT-1)+2,
     *           NXYA=NICV*NJCV*(4**NGIT-1)/3+NICV*2*(2**NGIT-1)+
     *           NJCV*2*(2**NGIT-1)+4*NGIT,NXY=NX*NY,NPHI=5)
      COMMON/INDEX/  NI,NJ,NIM,NJM,NIJ,KGIT,IST,JST,IJST,NITP,
     *      LI(NX),LIG(NX),LSG(NGIT),LSR(NGIT),LSI(NGIT),MIT(NGIT),
     *      NIGIT(NGIT),NJGIT(NGIT),IGIT(NGIT),JGIT(NGIT),IJGIT(NGIT)
      COMMON/LOGIC/ WRIT3,READ3,TEST,AKSI,LOUTS,LOUTE,ISOL,LSOR,
     *      ICAL(NPHI)
      COMMON/CONTRO/ IMON,JMON,IJMON,IU,IV,IP,IEN,ISW,ISF,IPR,JPR,kgs,
     *      ITBS(NXA),ITBN(NXA),JTBW(NYA),JTBE(NYA),beta,gravx,gravy,
     *      NSW(NPHI),SOR(NPHI),RESOR(NPHI),PRTR(NPHI),URF(NPHI),
     *      GDS(NPHI),SNOR(NPHI),SORMAX,SLARGE,ALFA,GREAT,SMALL,TREF
      COMMON/VARF/ DENSIT,VISCOS,FLOREF,XMOMR,PRMR,PRR(NXYA),
     *      U(NXYA),V(NXYA),P(NXYA),PP(NXYA),T(NXYA),DEN(NXYA),
     *      VIS(NXYA),F1(NXYA),F2(NXYA),DP1(NXYA),DP2(NXYA)
      COMMON/VARG/ U1(NXYA),V1(NXYA),T1(NXYA),F1G(NXYA),F2G(NXYA),
     *      RESU(NXYA),RESV(NXYA),REST(NXYA)
      COMMON /GEO/ X(NXYA),Y(NXYA),XC(NXYA),YC(NXYA),FX(NXYA),FY(NXYA),
     *      VOL(NXYA),R(NXYA),RC(NXYA)
      COMMON/COEF/ AE(NXYA),AW(NXYA),AN(NXYA),AS(NXYA),AP(NXYA),
     *      SU(NXYA),SV(NXYA),SP(NXYA),D1(NXYA),D2(NXYA),E1(NXYA),
     *      E2(NXYA)
      COMMON/BCOEF/ BW(NXYA),BE(NXYA),BS(NXYA),BN(NXYA),BP(NXYA),
     *      RES(NXYA)
      CHARACTER*10 FILRES,FILGRD,FILIN,FILOUT
      LOGICAL WRIT3,READ3,TEST,AKSI,LOUTS,LOUTE,ISOL,LSOR,ICAL
C
      PRINT *, 'ENTER INPUT FILE NAME:  '
      READ(*,5) FILIN
      PRINT *,' ENTER OUTPUT FILE NAME:  '
      READ(*,5) FILOUT
      PRINT *,' ENTER GRID FILE NAME:  '
      READ(*,5) FILGRD
      PRINT *,' ENTER RESULTS FILE NAME:  '
      READ(*,5) FILRES
    5 FORMAT(A10)
      OPEN (UNIT=5,FILE=FILIN)
      OPEN (UNIT=6,FILE=FILOUT)
      OPEN (UNIT=4,FILE=FILGRD,FORM='UNFORMATTED')
      OPEN (UNIT=3,FILE=FILRES,FORM='UNFORMATTED')
      REWIND 3
      REWIND 4
      REWIND 5
      REWIND 6
C
c     initialize
c
      kgs=0
      CALL MODINP
      CALL MOOUT1
C
      kgs=kgs+1
      DO 100 KGIT=kgs,NGIT
      print *,'  Begining grid No.',kgit
      MAXIT=MIT(KGIT)
      LSGIT=LSG(KGIT)
      NI=NIGIT(KGIT)
      NJ=NJGIT(KGIT)
      IJST=IJGIT(KGIT)
      IST=IGIT(KGIT)
      JST=JGIT(KGIT)
      NIJ=NI*NJ
      NIM=NI-1
      NJM=NJ-1
      DO 110 I=1,NI
  110 LI(I)=(I-1)*NJ+IJST
C
       IF (KGIT.GT.1) THEN
C----- INTERPOLATE COARSE TO FINE
C
      IF(ICAL(IU))  CALL VINT(KGIT-1,U,U)
      IF(ICAL(IV))  CALL VINT(KGIT-1,V,V)
      IF(ICAL(IP))  CALL VINT(KGIT-1,P,P)
      IF(ICAL(IEN)) CALL VINT(KGIT-1,T,T)
      IF(ICAL(ISW)) CALL VINT(KGIT-1,WR,WR)
C
      DO 120 I=2,NI-2
      IJS=LI(I)+2
      IJE=LI(I)+NJM
      DO 120 IJ=IJS,IJE
      IJM=IJ-1
      IPJ=IJ+NJ
      FXE=FX(IJ)
      FXP=1.-FXE
      DXE=X(IJ)-X(IJM)
      DYE=Y(IJ)-Y(IJM)
      RE=0.5*(R(IJ)+R(IJM))
      UE=FXE*U(IPJ)+FXP*U(IJ)
      VE=FXE*V(IPJ)+FXP*V(IJ)
      DENE=FXE*DEN(IPJ)+FXP*DEN(IJ)
      F1(IJ)=(UE*DYE-VE*DXE)*DENE*RE
  120 CONTINUE
      DO 160 I=2,NIM
      IJS=LI(I)+2
      IJE=LI(I)+NJ-2
      DO 160 IJ=IJS,IJE
      IMJ=IJ-NJ
      IJP=IJ+1
      FYN=FY(IJ)
      FYP=1.-FYN
      DXN=X(IJ)-X(IMJ)
      DYN=Y(IJ)-Y(IMJ)
      RN=0.5*(R(IJ)+R(IMJ))
      UN=FYN*U(IJP)+FYP*U(IJ)
      VN=FYN*V(IJP)+FYP*V(IJ)
      DENN=FYN*DEN(IJP)+FYP*DEN(IJ)
      F2(IJ)=(VN*DXN-UN*DYN)*DENN*RN
  160 CONTINUE
      ENDIF
      DO 170 IJ=IJST+1,IJST+NIJ
  170 PP(IJ)=P(IJ)
      IIM=2**(KGIT-1)*(IMON-1)+1
      JJM=2**(KGIT-1)*(JMON-1)+1
      IJM=LI(IIM)+JJM
C
      CALL BCIN
C
      IF (LOUTS) CALL OUTPUT
      WRITE(6,6000) IMON,JMON
      DO 200 ITER=1,MAXIT
      print *, '    Begining cycle No.',iter
      ISOL=.TRUE.
      LSOR=.FALSE.
      DO 210 LS=1,LSGIT
      IF(ICAL(IU)) CALL CALCUV(KGIT)
      IF(ICAL(IU)) CALL OUTBC
      IF(ICAL(IP)) CALL CALCP(KGIT)
      IF(ICAL(IEN)) CALL CALCT (KGIT)
      DO 220 L=1,NPHI
  220 RESOR(L)=RESOR(L)*SNOR(L)
      WRITE(6,76) KGIT,ITER,LS,RESOR(IU),RESOR(IV),RESOR(IP),
     *            RESOR(IEN),RESOR(ISW),U(IJM),V(IJM),P(IJM),
     *            T(IJM)
      SOURCE=AMAX1(RESOR(IU),RESOR(IV),RESOR(IP),RESOR(IEN),
     *             RESOR(ISW))
      IF(SOURCE.GT.SLARGE) GO TO 500
      IF(SOURCE.LT.SORMAX) GO TO 530
  210 CONTINUE
      CALL PRESB(KGIT)
C
C----- START A  V-CYCLE
C
      IF(KGIT.GT.1) THEN
      WRITE(6,*) '    '
      DO 230 KK=KGIT,2,-1
  230 CALL RESTR(KK)
      WRITE(6,*) '    '
C
      DO 240 KK=1,KGIT-1
  240 CALL INTER(KK)
      WRITE(6,*) '    '
      ENDIF
      SOURCE=AMAX1(RESOR(IU),RESOR(IV),RESOR(IP),RESOR(IEN),
     *             RESOR(ISW))
      IF(SOURCE.GT.SLARGE) GO TO 500
      IF(SOURCE.LT.SORMAX) GO TO 530
C
  200 CONTINUE
C----- V-CYCLE FINISHED OR WAS NOT NECESSARY
  530 CONTINUE
C
      CALL PRESB(KGIT)
      IF (LOUTE) CALL OUTPUT
      CALL MOOUT2
      IF(WRIT3) WRITE(3) ITER,ALFA,NI,NJ,NIM,NJM,NIJ,
     *   (X(IJ),IJ=IJST+1,IJST+NIJ),(Y(IJ),IJ=IJST+1,IJST+NIJ),
     *   (XC(IJ),IJ=IJST+1,IJST+NIJ),(YC(IJ),IJ=IJST+1,IJST+NIJ),
     *   (F1(IJ),IJ=IJST+1,IJST+NIJ),(F2(IJ),IJ=IJST+1,IJST+NIJ),
     *   (U(IJ),IJ=IJST+1,IJST+NIJ),(V(IJ),IJ=IJST+1,IJST+NIJ),
     *   (P(IJ),IJ=IJST+1,IJST+NIJ),(T(IJ),IJ=IJST+1,IJST+NIJ)
  100 CONTINUE
      STOP
C
  500 PRINT *,'     *** TERMINATED - INNER ITERATIONS DIVERGE ***'
 6001 FORMAT(2X,I4,3X,1P3E10.3,31X,1P4E11.3)
 6007 FORMAT(//,10X,'TIME SPENT FOR GRID ',I2,': ',1PE12.4,/,
     *      10X,36(1H*),/)
   60 FORMAT(20X,'TIME STEP NO. ',I3,' ,  TIME =',E10.4,/,20X,
     *       36(1H*),/)
 6000 FORMAT(1X,'GRID',1X,'ITER',1X,'LSG',2X,
     *'I--------ABSOLUTE RESIDUAL SOURCE SUMS-----------I',3X,
     *'I----FIELD VALUES AT MONITORING LOCATION(',1X,I3,',',I3,')---I',
     */,1X,3('NO.',2X),5X,'UMOM',6X,'VMOM',6X,'MASS',6X,'ENER',7X,'WR',
     *       17X,'U',9X,'V',9X,'P',9X,'T',8X,'WR',/)
   76 FORMAT(1X,I2,2X,I4,1X,I3,6X,1P5E10.3,10X,1P5E10.3)
      IF(LOUTE) CALL OUTPUT
      STOP
      END
C
C***********************************************************
      SUBROUTINE RESTR(K)
C***********************************************************
      PARAMETER (ngit=3,nicv=21,njcv=8,NXA=NICV*(2**NGIT-1)+
     *           2*NGIT,NYA=NJCV*(2**NGIT-1)+2*NGIT,
     *           NX=NICV*2**(NGIT-1)+2,NY=NJCV*2**(NGIT-1)+2,
     *           NXYA=NICV*NJCV*(4**NGIT-1)/3+NICV*2*(2**NGIT-1)+
     *           NJCV*2*(2**NGIT-1)+4*NGIT,NXY=NX*NY,NPHI=5)
      COMMON/INDEX/  NI,NJ,NIM,NJM,NIJ,KGIT,IST,JST,IJST,NITP,
     *      LI(NX),LIG(NX),LSG(NGIT),LSR(NGIT),LSI(NGIT),MIT(NGIT),
     *      NIGIT(NGIT),NJGIT(NGIT),IGIT(NGIT),JGIT(NGIT),IJGIT(NGIT)
      COMMON/LOGIC/ WRIT3,READ3,TEST,AKSI,LOUTS,LOUTE,ISOL,LSOR,
     *      ICAL(NPHI)
      COMMON/CONTRO/ IMON,JMON,IJMON,IU,IV,IP,IEN,ISW,ISF,IPR,JPR,kgs,
     *      ITBS(NXA),ITBN(NXA),JTBW(NYA),JTBE(NYA),beta,gravx,gravy,
     *      NSW(NPHI),SOR(NPHI),RESOR(NPHI),PRTR(NPHI),URF(NPHI),
     *      GDS(NPHI),SNOR(NPHI),SORMAX,SLARGE,ALFA,GREAT,SMALL,TREF
      COMMON/VARF/ DENSIT,VISCOS,FLOREF,XMOMR,PRMR,PRR(NXYA),
     *      U(NXYA),V(NXYA),P(NXYA),PP(NXYA),T(NXYA),DEN(NXYA),
     *      VIS(NXYA),F1(NXYA),F2(NXYA),DP1(NXYA),DP2(NXYA)
      COMMON/VARG/ U1(NXYA),V1(NXYA),T1(NXYA),F1G(NXYA),F2G(NXYA),
     *      RESU(NXYA),RESV(NXYA),REST(NXYA)
      COMMON /GEO/ X(NXYA),Y(NXYA),XC(NXYA),YC(NXYA),FX(NXYA),FY(NXYA),
     *      VOL(NXYA),R(NXYA),RC(NXYA)
      COMMON/COEF/ AE(NXYA),AW(NXYA),AN(NXYA),AS(NXYA),AP(NXYA),
     *      SU(NXYA),SV(NXYA),SP(NXYA),D1(NXYA),D2(NXYA),E1(NXYA),
     *      E2(NXYA)
      COMMON/BCOEF/ BW(NXYA),BE(NXYA),BS(NXYA),BN(NXYA),BP(NXYA),
     *      RES(NXYA)
      COMMON /SYMB/ FIB,FIP,FIL,FIR,FAC,IJN
      LOGICAL WRIT3,READ3,TEST,AKSI,LOUTS,LOUTE,ISOL,LSOR,ICAL
C
      IJSTF=IJGIT(K)
      IJSTG=IJGIT(K-1)
      NIF=NIGIT(K)
      NIG=NIGIT(K-1)
      NJF=NJGIT(K)
      NJG=NJGIT(K-1)
      NIMG=NIG-1
      NJMG=NJG-1
      NIMF=NIF-1
      NJMF=NJF-1
      NI=NIF
      NJ=NJF
      NIJ=NI*NJ
      NIM=NI-1
      NJM=NJ-1
      IJST=IJSTF
      DO 5 I=1,NIF
    5 LI(I)=(I-1)*NJF+IJSTF
      DO 10 I=1,NIG
   10 LIG(I)=(I-1)*NJG+IJSTG
C----- Calculate only new coefficients
      ISOL=.FALSE.
      IF (ICAL(IU)) THEN
      CALL CALCUV(K)
C----- CALCULATE RESU,RESV
      DO 15 I=2,NIMF
      DO 15 J=2,NJMF
      IJ=LI(I)+J
      APTU=AE(IJ)+AW(IJ)+AN(IJ)+AS(IJ)+AP(IJ)
      APTV=AE(IJ)+AW(IJ)+AN(IJ)+AS(IJ)+SP(IJ)
      SU(IJ)=SU(IJ)+RESU(IJ)+AE(IJ)*U(IJ+NJF)+AW(IJ)*U(IJ-NJF)+
     *       AN(IJ)*U(IJ+1)+AS(IJ)*U(IJ-1)-APTU*U(IJ)
      SV(IJ)=SV(IJ)+RESV(IJ)+AS(IJ)*V(IJ-1)+AN(IJ)*V(IJ+1)+
     *       AE(IJ)*V(IJ+NJF)+AW(IJ)*V(IJ-NJF)-APTV*V(IJ)
   15 CONTINUE
C----- CALCULATIONS FOR COARSE GRID
      DO 20 I=2,NIMG
      IPF=I*2-1
      IMF=(I-1)*2
      DO 20 J=2,NJMG
      JPF=J*2-1
      JMF=(J-1)*2
      IJ=LIG(I)+J
      IPJPF=LI(IPF)+JPF
      IMJMF=LI(IMF)+JMF
      F1(IJ)=F1(IPJPF)+F1(IPJPF-1)
      F2(IJ)=F2(IPJPF)+F2(IMJMF+1)
      F1G(IJ)=F1(IJ)
      F2G(IJ)=F2(IJ)
      RESU(IJ)=SU(IPJPF)+SU(IPJPF-1)+SU(IMJMF)+SU(IMJMF+1)
      RESV(IJ)=SV(IPJPF)+SV(IPJPF-1)+SV(IMJMF)+SV(IMJMF+1)
   20 CONTINUE
C
      DO 25 J=2,NJMG
      JPF=J*2-1
      IJ=LIG(1)+J
      IJPF=LI(1)+JPF
      F1(IJ)=F1(IJPF)+F1(IJPF-1)
      F1G(IJ)=F1(IJ)
      IJ=LIG(NIMG)+J
      IJPF=LI(NIM)+JPF
      F1(IJ)=F1(IJPF)+F1(IJPF-1)
      F1G(IJ)=F1(IJ)
   25 CONTINUE
      ENDIF
C
C-----CALCULATE RESIDUALS FOR T AND W
C
      IF(ICAL(IEN)) THEN
      CALL CALCT(K)
      DO 30 I=2,NIMF
      IJS=LI(I)+2
      IJE=LI(I)+NJM
      DO 30 IJ=IJS,IJE
      APTS=AE(IJ)+AW(IJ)+AN(IJ)+AS(IJ)+SP(IJ)
      SU(IJ)=SU(IJ)+REST(IJ)+AE(IJ)*T(IJ+NJF)+AW(IJ)*T(IJ-NJF)+
     *       AN(IJ)*T(IJ+1)+AS(IJ)*T(IJ-1)-APTS*T(IJ)
   30 CONTINUE
C
C-----TRANSFER RESIDUALS TO COARSE GRID
C
      DO 35 I=2,NIMG
      IPF=I*2-1
      IMF=(I-1)*2
      DO 35 J=2,NJMG
      JPF=J*2-1
      JMF=(J-1)*2
      IJ=LIG(I)+J
      IPJPF=LI(IPF)+JPF
      IMJMF=LI(IMF)+JMF
      REST(IJ)=SU(IPJPF)+SU(IPJPF-1)+SU(IMJMF)+
     *         SU(IMJMF+1)
   35 CONTINUE
      ENDIF
C----- RESTRICT FINE GRID VALUES TO COARSE GRID
      IF (ICAL(IU)) CALL INTFG (K,U,U1)
      IF (ICAL(IV)) CALL INTFG (K,V,V1)
      IF(ICAL(IEN)) CALL INTFG(K,T,T1)
C
      NI=NIG
      NJ=NJG
      NIM=NI-1
      NJM=NJ-1
      NIJ=NI*NJ
      IJST=IJSTG
C----- COARSE GRID BECOMES FINE GRID
      DO 50 I=1,NI
   50 LI(I)=LIG(I)
C----- Symmetry boundary conditions for restricted variables
      IST=IGIT(K-1)
      JST=JGIT(K-1)
      DO 55 I=2,NIM
      ij=li(i)+2
      ibc=itbs(i+ist)
      IF(ibc.eq.3.or.ibc.eq.5) CALL BC1(IJ,IJ-1,IJ-NJ,IJ+NJ,IJ-1,ibc)
      ij=li(i)+njm
      ibc=itbn(i+ist)
      IF(ibc.eq.3.or.ibc.eq.5) CALL BC1(IJ,IJ+1,IJ+NJ,IJ-NJ,IJ-NJ,ibc)
   55 CONTINUE
      DO 57 J=2,NJM
      ij=li(2)+j
      ibc=jtbw(j+jst)
      IF(ibc.eq.3.or.ibc.eq.5) CALL BC1(IJ,IJ-NJ,IJ+1,IJ-1,IJ-NJ-1,ibc)
      ij=li(nim)+j
      ibc=jtbe(j+jst)
      IF(ibc.eq.3.or.ibc.eq.5) CALL BC1(IJ,IJ+NJ,IJ-1,IJ+1,IJ,ibc)
   57 CONTINUE
C
      DO 60 IJ=IJST+1,IJST+NIJ
      U (IJ)=U1(IJ)
      V (IJ)=V1(IJ)
      PP(IJ)=0.
      P (IJ)=0.
      T(IJ)=T1(IJ)
   60 CONTINUE
C
      ISOL=.TRUE.
      LSOR=.TRUE.
      NS=LSR(K-1)
      IIM=2**(KGIT-1)*(IMON-1)+1
      JJM=2**(KGIT-1)*(JMON-1)+1
      IJM=LI(IIM)+JJM
C
      DO 70 LS=1,NS
C----- START RESTRICTION-CYCLE
      IF (ICAL(IU)) CALL CALCUV (K-1)
      IF (ICAL(IP)) CALL CALCP (K-1)
      IF(ICAL(IEN)) CALL CALCT(K-1)
      DO 80 L=1,NPHI
   80 RESOR(L)=RESOR(L)*SNOR(L)
      WRITE(6,600) K-1,LS,RESOR(IU),RESOR(IV),RESOR(IP),RESOR(IEN),
     *             RESOR(ISW),U(IJM),V(IJM),P(IJM),T(IJM)
      LSOR=.FALSE.
   70 CONTINUE
      CALL PRESB(K-1)
  600 FORMAT(1X,I2,6X,I4,1X,5HRESTR,1P5E10.3,10X,1P5E10.3)
C
      RETURN
      END
C
C********************************************************
      SUBROUTINE INTFG(K,FI,FIG)
C********************************************************
      PARAMETER (ngit=3,nicv=21,njcv=8,NXA=NICV*(2**NGIT-1)+
     *           2*NGIT,NYA=NJCV*(2**NGIT-1)+2*NGIT,
     *           NX=NICV*2**(NGIT-1)+2,NY=NJCV*2**(NGIT-1)+2,
     *           NXYA=NICV*NJCV*(4**NGIT-1)/3+NICV*2*(2**NGIT-1)+
     *           NJCV*2*(2**NGIT-1)+4*NGIT,NXY=NX*NY,NPHI=5)
      COMMON/INDEX/  NI,NJ,NIM,NJM,NIJ,KGIT,IST,JST,IJST,NITP,
     *      LI(NX),LIG(NX),
     *      LSG(NGIT),LSR(NGIT),LSI(NGIT),MIT(NGIT),
     *      NIGIT(NGIT),NJGIT(NGIT),IGIT(NGIT),JGIT(NGIT),IJGIT(NGIT)
      COMMON/LOGIC/  WRIT3,READ3,TEST,AKSI,LOUTS,LOUTE,ISOL,LSOR,
     *      ICAL(NPHI)
      COMMON/CONTRO/ IMON,JMON,IJMON,IU,IV,IP,IEN,ISW,ISF,IPR,JPR,kgs,
     *      ITBS(NXA),ITBN(NXA),JTBW(NYA),JTBE(NYA),beta,gravx,gravy,
     *      NSW(NPHI),SOR(NPHI),RESOR(NPHI),PRTR(NPHI),URF(NPHI),
     *      GDS(NPHI),SNOR(NPHI),SORMAX,SLARGE,ALFA,GREAT,SMALL,TREF
      COMMON/VARF/ DENSIT,VISCOS,FLOREF,XMOMR,PRMR,PRR(NXYA),
     *      U(NXYA),V(NXYA),P(NXYA),PP(NXYA),T(NXYA),DEN(NXYA),
     *      VIS(NXYA),F1(NXYA),F2(NXYA),DP1(NXYA),DP2(NXYA)
      COMMON/VARG/ U1(NXYA),V1(NXYA),T1(NXYA),F1G(NXYA),F2G(NXYA),
     *      RESU(NXYA),RESV(NXYA),REST(NXYA)
      COMMON /GEO/ X(NXYA),Y(NXYA),XC(NXYA),YC(NXYA),FX(NXYA),FY(NXYA),
     *      VOL(NXYA),R(NXYA),RC(NXYA)
      COMMON/COEF/ AE(NXYA),AW(NXYA),AN(NXYA),AS(NXYA),AP(NXYA),
     *      SU(NXYA),SV(NXYA),SP(NXYA),D1(NXYA),D2(NXYA),E1(NXYA),
     *      E2(NXYA)
      COMMON/BCOEF/ BW(NXYA),BE(NXYA),BS(NXYA),BN(NXYA),BP(NXYA),
     *      RES(NXYA)
      LOGICAL WRIT3,READ3,TEST,AKSI,LOUTS,LOUTE,ISOL,LSOR,ICAL
      DIMENSION FI(NXYA),FIG(NXYA)
C
      IJSTF=IJGIT(K)
      IJSTG=IJGIT(K-1)
      NIF=NIGIT(K)
      NIG=NIGIT(K-1)
      NJF=NJGIT(K)
      NJG=NJGIT(K-1)
      NIMG=NIG-1
      NJMG=NJG-1
      DO 5 I=1,NIF
    5 LI(I)=(I-1)*NJF+IJSTF
      DO 10 I=1,NIG
   10 LIG(I)=(I-1)*NJG+IJSTG
C----- INDICES AND CONTROL PARAMETERS
      DO 100 IG=2,NIMG
      IIG=LIG(IG)
      I=2*IG-2
      II=LI(I)
      DO 100 JG=2,NJMG
      IJG=IIG+JG
      J=2*JG-2
      IJ=II+J
C----- INTERPOLATE FOR INTERIOR COARSE GRID
      IJE=IJ+NJF
      IJNE=IJE+1
      FIG(IJG)=FX(IJ)*( FY(IJ)*FI(IJ) + (1.-FY(IJ))*FI(IJ+1) )+
     *       (1.-FX(IJ))*( FY(IJ)*FI(IJE)+(1.-FY(IJ))*FI(IJNE) )
  100 CONTINUE
C----- SOUTH BOUNDARY
      DO 120 IG=2,NIMG
      I=2*IG-2
      IJG=LIG(IG)+1
      IJ=LI(I)+1
      FIG(IJG)=FX(IJ)*FI(IJ)+(1.-FX(IJ))*FI(IJ+NJF)
C----- NORTH BOUNDARY
      IJG=LIG(IG)+NJG
      IJ=LI(I)+NJF
      FIG(IJG)=FX(IJ)*FI(IJ)+(1.-FX(IJ))*FI(IJ+NJF)
  120 CONTINUE
C----- WEST BOUNDARY
      DO 130 JG=2,NJMG
      J=2*JG-2
      IJ=LI(1)+J
      IJG=LIG(1)+JG
      FIG(IJG)=FY(IJ)*FI(IJ)+(1.-FY(IJ))*FI(IJ+1)
C----- EAST BOUNDARY
      IJG=LIG(NIG)+JG
      IJ=LI(NIF)+J
      FIG(IJG)=FY(IJ)*FI(IJ)+(1.-FY(IJ))*FI(IJ+1)
  130 CONTINUE
C
C-----CORNER
      FIG(LIG(1)+1)=FI(LI(1)+1)
      FIG(LIG(1)+NJG)=FI(LI(1)+NJF)
      FIG(LIG(NIG)+1)=FI(LI(NIF)+1)
      FIG(LIG(NIG)+NJG)=FI(LI(NIF)+NJF)
      RETURN
      END
C
C*********************************************************
      SUBROUTINE INTER(K)
C*********************************************************
      PARAMETER (ngit=3,nicv=21,njcv=8,NXA=NICV*(2**NGIT-1)+
     *           2*NGIT,NYA=NJCV*(2**NGIT-1)+2*NGIT,
     *           NX=NICV*2**(NGIT-1)+2,NY=NJCV*2**(NGIT-1)+2,
     *           NXYA=NICV*NJCV*(4**NGIT-1)/3+NICV*2*(2**NGIT-1)+
     *           NJCV*2*(2**NGIT-1)+4*NGIT,NXY=NX*NY,NPHI=5)
      COMMON/INDEX/  NI,NJ,NIM,NJM,NIJ,KGIT,IST,JST,IJST,NITP,
     *      LI(NX),LIG(NX),
     *      LSG(NGIT),LSR(NGIT),LSI(NGIT),MIT(NGIT),
     *      NIGIT(NGIT),NJGIT(NGIT),IGIT(NGIT),JGIT(NGIT),IJGIT(NGIT)
      COMMON/LOGIC/  WRIT3,READ3,TEST,AKSI,LOUTS,LOUTE,ISOL,LSOR,
     *      ICAL(NPHI)
      COMMON/CONTRO/ IMON,JMON,IJMON,IU,IV,IP,IEN,ISW,ISF,IPR,JPR,kgs,
     *      ITBS(NXA),ITBN(NXA),JTBW(NYA),JTBE(NYA),beta,gravx,gravy,
     *      NSW(NPHI),SOR(NPHI),RESOR(NPHI),PRTR(NPHI),URF(NPHI),
     *      GDS(NPHI),SNOR(NPHI),SORMAX,SLARGE,ALFA,GREAT,SMALL,TREF
      COMMON/VARF/ DENSIT,VISCOS,FLOREF,XMOMR,PRMR,PRR(NXYA),
     *      U(NXYA),V(NXYA),P(NXYA),PP(NXYA),T(NXYA),DEN(NXYA),
     *      VIS(NXYA),F1(NXYA),F2(NXYA),DP1(NXYA),DP2(NXYA)
      COMMON/VARG/ U1(NXYA),V1(NXYA),T1(NXYA),F1G(NXYA),F2G(NXYA),
     *      RESU(NXYA),RESV(NXYA),REST(NXYA)
      COMMON /GEO/ X(NXYA),Y(NXYA),XC(NXYA),YC(NXYA),FX(NXYA),FY(NXYA),
     *      VOL(NXYA),R(NXYA),RC(NXYA)
      COMMON/COEF/ AE(NXYA),AW(NXYA),AN(NXYA),AS(NXYA),AP(NXYA),
     *      SU(NXYA),SV(NXYA),SP(NXYA),D1(NXYA),D2(NXYA),E1(NXYA),
     *      E2(NXYA)
      COMMON/BCOEF/ BW(NXYA),BE(NXYA),BS(NXYA),BN(NXYA),BP(NXYA),
     *      RES(NXYA)
      LOGICAL WRIT3,READ3,TEST,AKSI,LOUTS,LOUTE,ISOL,LSOR,ICAL
      DIMENSION VP(NXYA),UP(NXYA)
C
      IJSTG=IJGIT(K)
      IJSTF=IJGIT(K+1)
      NIG=NIGIT(K)
      NJG=NJGIT(K)
      NIF=NIGIT(K+1)
      NJF=NJGIT(K+1)
      IST=IGIT(K+1)
      JST=JGIT(K+1)
      NIMF=NIF-1
      NJMF=NJF-1
C
      NIJF=NIF*NJF
      NIJG=NIG*NJG
      NI=NIF
      NJ=NJF
      NIJ=NI*NJ
      NIM=NI-1
      NJM=NJ-1
      IJST=IJSTF
C
      DO 5 I=1,NIF
    5 LI(I)=(I-1)*NJF+IJSTF
      DO 10 I=1,NIG
   10 LIG(I)=(I-1)*NJG+IJSTG
C
C.....INTERPOLATE VELOCITY CORRECTIONS
C
      IF (ICAL(IU)) THEN
      DO 15 IJ=IJSTG+1,IJSTG+NIJG
   15 UP(IJ)=U(IJ)-U1(IJ)
      CALL VINT(K,UP,UP)
      DO 20 IJ=IJSTF+1,IJSTF+NIJF
   20 U(IJ)=U(IJ)+UP(IJ)
      ENDIF
      IF (ICAL(IV)) THEN
      DO 24 IJ=IJSTG+1,IJSTG+NIJG
   24 VP(IJ)=V(IJ)-V1(IJ)
      CALL VINT(K,VP,VP)
      DO 27 IJ=IJSTF+1,IJSTF+NIJF
   27 V(IJ)=V(IJ)+VP(IJ)
      ENDIF
C
C..... CALCULATE MASS FLUXES CORRECTIONS
C
      IF(ICAL(IU).OR.ICAL(IV)) THEN
      DO 30 I=2,NIF-2
      IJS=LI(I)+2
      IJE=LI(I)+NJMF
      DO 30 IJ=IJS,IJE
      IJM=IJ-1
      IPJ=IJ+NJ
      FXE=FX(IJ)
      FXP=1.-FXE
      DXE=X(IJ)-X(IJM)
      DYE=Y(IJ)-Y(IJM)
      RE=0.5*(R(IJ)+R(IJM))
      UE=FXE*UP(IPJ)+FXP*UP(IJ)
      VE=FXE*VP(IPJ)+FXP*VP(IJ)
      DENE=FXE*DEN(IPJ)+FXP*DEN(IJ)
      F1(IJ)=F1(IJ)+(UE*DYE-VE*DXE)*DENE*RE
   30 CONTINUE
      DO 35 I=2,NIMF
      IJS=LI(I)+2
      IJE=LI(I)+NJF-2
      DO 35 IJ=IJS,IJE
      IMJ=IJ-NJ
      IJP=IJ+1
      FYN=FY(IJ)
      FYP=1.-FYN
      DXN=X(IJ)-X(IMJ)
      DYN=Y(IJ)-Y(IMJ)
      RN=0.5*(R(IJ)+R(IMJ))
      UN=FYN*UP(IJP)+FYP*UP(IJ)
      VN=FYN*VP(IJP)+FYP*VP(IJ)
      DENN=FYN*DEN(IJP)+FYP*DEN(IJ)
      F2(IJ)=F2(IJ)+(VN*DXN-UN*DYN)*DENN*RN
   35 CONTINUE
      ENDIF
C-----T-CORRECTIONS
C
      IF(ICAL(IEN)) THEN
      DO 60 IJ=IJSTG+1,IJSTG+NIJG
   60 VP(IJ)=T(IJ)-T1(IJ)
      CALL VINT(K,VP,VP)
      DO 65 IJ=IJSTF+1,IJSTF+NIJF
   65 T(IJ)=T(IJ)+VP(IJ)
      ENDIF
C
C.....INTERPOLATE PRESSURE CORRECTIONS
C
      IF (ICAL(IP)) THEN
      CALL VINT(K,P,VP)
      DO 50 IJ=IJSTF+1,IJSTF+NIJF
      P(IJ)=P(IJ)+VP(IJ)
      PP(IJ)=P(IJ)
   50 CONTINUE
      ENDIF
C
C.....SOLVE ON FINE GRID (DO NS SWEEPS)
C
      NS=LSI(K+1)
      ISOL=.TRUE.
      IIM=2**(KGIT-1)*(IMON-1)+1
      JJM=2**(KGIT-1)*(JMON-1)+1
      IJM=LI(IIM)+JJM
C
      DO 80 LS=1,NS
      IF (ICAL(IU)) CALL CALCUV (K+1)
      IF (ICAL(IU).AND.K+1.EQ.KGIT) CALL OUTBC
      IF (ICAL(IP)) CALL CALCP (K+1)
      IF(ICAL(IEN)) CALL CALCT(K+1)
      DO 90 L=1,NPHI
   90 RESOR(L)=RESOR(L)*SNOR(L)
      WRITE(6,600) K+1,LS,RESOR(IU),RESOR(IV),RESOR(IP),RESOR(IEN),
     *             RESOR(ISW),U(IJM),V(IJM),P(IJM),T(IJM)
      SOURCE=AMAX1(RESOR(IU),RESOR(IV),RESOR(IP),RESOR(IEN),
     *             RESOR(ISW))
      IF(SOURCE.LT.SORMAX.and.k+1.eq.kgit) return
   80 CONTINUE
      CALL PRESB(K)
  600 FORMAT(1X,I2,6X,I4,1X,5HINTER,1P5E10.3,10X,1P5E10.3)
      RETURN
      END
C
C**********************************************************
      SUBROUTINE VINT(K,FIG,FI)
C**********************************************************
      PARAMETER (ngit=3,nicv=21,njcv=8,NXA=NICV*(2**NGIT-1)+
     *           2*NGIT,NYA=NJCV*(2**NGIT-1)+2*NGIT,
     *           NX=NICV*2**(NGIT-1)+2,NY=NJCV*2**(NGIT-1)+2,
     *           NXYA=NICV*NJCV*(4**NGIT-1)/3+NICV*2*(2**NGIT-1)+
     *           NJCV*2*(2**NGIT-1)+4*NGIT,NXY=NX*NY,NPHI=5)
      COMMON/INDEX/  NI,NJ,NIM,NJM,NIJ,KGIT,IST,JST,IJST,NITP,
     *      LI(NX),LIG(NX),
     *      LSG(NGIT),LSR(NGIT),LSI(NGIT),MIT(NGIT),
     *      NIGIT(NGIT),NJGIT(NGIT),IGIT(NGIT),JGIT(NGIT),IJGIT(NGIT)
      COMMON/LOGIC/  WRIT3,READ3,TEST,AKSI,LOUTS,LOUTE,ISOL,LSOR,
     *      ICAL(NPHI)
      COMMON/CONTRO/ IMON,JMON,IJMON,IU,IV,IP,IEN,ISW,ISF,IPR,JPR,kgs,
     *      ITBS(NXA),ITBN(NXA),JTBW(NYA),JTBE(NYA),beta,gravx,gravy,
     *      NSW(NPHI),SOR(NPHI),RESOR(NPHI),PRTR(NPHI),URF(NPHI),
     *      GDS(NPHI),SNOR(NPHI),SORMAX,SLARGE,ALFA,GREAT,SMALL,TREF
      COMMON/VARF/ DENSIT,VISCOS,FLOREF,XMOMR,PRMR,PRR(NXYA),
     *      U(NXYA),V(NXYA),P(NXYA),PP(NXYA),T(NXYA),DEN(NXYA),
     *      VIS(NXYA),F1(NXYA),F2(NXYA),DP1(NXYA),DP2(NXYA)
      COMMON/VARG/ U1(NXYA),V1(NXYA),T1(NXYA),F1G(NXYA),F2G(NXYA),
     *      RESU(NXYA),RESV(NXYA),REST(NXYA)
      COMMON /GEO/ X(NXYA),Y(NXYA),XC(NXYA),YC(NXYA),FX(NXYA),FY(NXYA),
     *      VOL(NXYA),R(NXYA),RC(NXYA)
      COMMON/COEF/ AE(NXYA),AW(NXYA),AN(NXYA),AS(NXYA),AP(NXYA),
     *      SU(NXYA),SV(NXYA),SP(NXYA),D1(NXYA),D2(NXYA),E1(NXYA),
     *      E2(NXYA)
      COMMON/BCOEF/ BW(NXYA),BE(NXYA),BS(NXYA),BN(NXYA),BP(NXYA),
     *      RES(NXYA)
      LOGICAL WRIT3,READ3,TEST,AKSI,LOUTS,LOUTE,ISOL,LSOR,ICAL
      DIMENSION FI(NXYA),FIG(NXYA)
C
      IJSTG=IJGIT(K)
      IJSTF=IJGIT(K+1)
      NIG=NIGIT(K)
      NIF=NIGIT(K+1)
      NJG=NJGIT(K)
      NJF=NJGIT(K+1)
      NIMG=NIG-1
      NJMG=NJG-1
C
      DO 5 I=1,NIF
    5 LI(I)=(I-1)*NJF+IJSTF
      DO 8 I=1,NIG
    8 LIG(I)=(I-1)*NJG+IJSTG
C----- MAIN LOOP
C----- START INTERPOLATION FOR FINE GRID NODES
      DO 200 IG=1,NIMG
      IIG=LIG(IG)
      I=2*IG-1
      II=LI(I)
      DO 200 JG=1,NJMG
      IJG=IIG+JG
      J=2*JG-1
      IJ=II+J
C----- PARAMETERS FOR POINTS BETWEEN FOUR COARSE GRID NODES
      IJP=IJ
      IJN=IJP+1
      IJE=IJP+NJF
      IJEN=IJE+1
      IJSW=IJG
      IJNW=IJSW+1
      IJSE=IJSW+NJG
      IJNE=IJSE+1
C----- INTERPOLATION FACTORS
      FY1L=FY(IJSW)*FY(IJP-1)
      FY1R=FY(IJSE)*FY(IJE-1)
      FY2L=(1.-FY(IJSW))*(1.-FY(IJN))
      FY2R=(1.-FY(IJSE))*(1.-FY(IJEN))
      FX1D=FX(IJSW)*FX(IJP-NJF)
      FX1U=FX(IJNW)*FX(IJN-NJF)
      FX2D=(1.-FX(IJSW))*(1.-FX(IJE))
      FX2U=(1.-FX(IJNW))*(1.-FX(IJEN))
C----- INTERPOLATE FOR VARIOUS VARIABLES
      FILD=FY1L*FIG(IJNW)+(1.-FY1L)*FIG(IJSW)
      FIRD=FY1R*FIG(IJNE)+(1.-FY1R)*FIG(IJSE)
      FILU=FY2L*FIG(IJSW)+(1.-FY2L)*FIG(IJNW)
      FIRU=FY2R*FIG(IJSE)+(1.-FY2R)*FIG(IJNE)
C
      FI(IJP)=FX1D*FIRD+(1.-FX1D)*FILD
      FI(IJN)=FX1U*FIRU+(1.-FX1U)*FILU
      FI(IJE)=FX2D*FILD+(1.-FX2D)*FIRD
      FI(IJEN)=FX2U*FILU+(1.-FX2U)*FIRU
  200 CONTINUE
      RETURN
      END
C
C***********************************************************
      SUBROUTINE STREAM
C***********************************************************
      PARAMETER (ngit=3,nicv=21,njcv=8,NXA=NICV*(2**NGIT-1)+
     *           2*NGIT,NYA=NJCV*(2**NGIT-1)+2*NGIT,
     *           NX=NICV*2**(NGIT-1)+2,NY=NJCV*2**(NGIT-1)+2,
     *           NXYA=NICV*NJCV*(4**NGIT-1)/3+NICV*2*(2**NGIT-1)+
     *           NJCV*2*(2**NGIT-1)+4*NGIT,NXY=NX*NY,NPHI=5)
      COMMON/INDEX/  NI,NJ,NIM,NJM,NIJ,KGIT,IST,JST,IJST,NITP,
     *      LI(NX),LIG(NX),
     *      LSG(NGIT),LSR(NGIT),LSI(NGIT),MIT(NGIT),
     *      NIGIT(NGIT),NJGIT(NGIT),IGIT(NGIT),JGIT(NGIT),IJGIT(NGIT)
      COMMON/VARF/ DENSIT,VISCOS,FLOREF,XMOMR,PRMR,PRR(NXYA),
     *      U(NXYA),V(NXYA),P(NXYA),PP(NXYA),T(NXYA),DEN(NXYA),
     *      VIS(NXYA),F1(NXYA),F2(NXYA),DP1(NXYA),DP2(NXYA)
      COMMON/COEF/ AE(NXYA),AW(NXYA),AN(NXYA),AS(NXYA),AP(NXYA),
     *      SU(NXYA),SV(NXYA),SP(NXYA),D1(NXYA),D2(NXYA),E1(NXYA),
     *      E2(NXYA)
C
      SP(LI(1)+1)=0.
      DO 10 I=1,NIM
      II=LI(I)
      IF(I.NE.1) SP(II+1)=SP(II-NJ+1)+F2(II+1)
      DO 10 J=2,NJ
      IJ=II+J
      SP(IJ)=SP(IJ-1)+F1(IJ)
   10 CONTINUE
C
      II=LI(NI)
      DO 20 J=1,NJ
      IJ=II+J
      SP(IJ)=SP(IJ-NJ)
   20 CONTINUE
      RETURN
      END
C
C*******************************************************
      SUBROUTINE CALCUV(KKK)
C*******************************************************
      PARAMETER (ngit=3,nicv=21,njcv=8,NXA=NICV*(2**NGIT-1)+
     *           2*NGIT,NYA=NJCV*(2**NGIT-1)+2*NGIT,
     *           NX=NICV*2**(NGIT-1)+2,NY=NJCV*2**(NGIT-1)+2,
     *           NXYA=NICV*NJCV*(4**NGIT-1)/3+NICV*2*(2**NGIT-1)+
     *           NJCV*2*(2**NGIT-1)+4*NGIT,NXY=NX*NY,NPHI=5)
      COMMON/INDEX/  NI,NJ,NIM,NJM,NIJ,KGIT,IST,JST,IJST,NITP,
     *      LI(NX),LIG(NX),
     *      LSG(NGIT),LSR(NGIT),LSI(NGIT),MIT(NGIT),
     *      NIGIT(NGIT),NJGIT(NGIT),IGIT(NGIT),JGIT(NGIT),IJGIT(NGIT)
      COMMON/LOGIC/  WRIT3,READ3,TEST,AKSI,LOUTS,LOUTE,ISOL,LSOR,
     *      ICAL(NPHI)
      COMMON/CONTRO/ IMON,JMON,IJMON,IU,IV,IP,IEN,ISW,ISF,IPR,JPR,kgs,
     *      ITBS(NXA),ITBN(NXA),JTBW(NYA),JTBE(NYA),beta,gravx,gravy,
     *      NSW(NPHI),SOR(NPHI),RESOR(NPHI),PRTR(NPHI),URF(NPHI),
     *      GDS(NPHI),SNOR(NPHI),SORMAX,SLARGE,ALFA,GREAT,SMALL,TREF
      COMMON/VARF/ DENSIT,VISCOS,FLOREF,XMOMR,PRMR,PRR(NXYA),
     *      U(NXYA),V(NXYA),P(NXYA),PP(NXYA),T(NXYA),DEN(NXYA),
     *      VIS(NXYA),F1(NXYA),F2(NXYA),DP1(NXYA),DP2(NXYA)
      COMMON/VARG/ U1(NXYA),V1(NXYA),T1(NXYA),F1G(NXYA),F2G(NXYA),
     *      RESU(NXYA),RESV(NXYA),REST(NXYA)
      COMMON /GEO/ X(NXYA),Y(NXYA),XC(NXYA),YC(NXYA),FX(NXYA),FY(NXYA),
     *      VOL(NXYA),R(NXYA),RC(NXYA)
      COMMON/COEF/ AE(NXYA),AW(NXYA),AN(NXYA),AS(NXYA),AP(NXYA),
     *      SU(NXYA),SV(NXYA),SP(NXYA),D1(NXYA),D2(NXYA),E1(NXYA),
     *      E2(NXYA)
      COMMON/BCOEF/ BW(NXYA),BE(NXYA),BS(NXYA),BN(NXYA),BP(NXYA),
     *      RES(NXYA)
      LOGICAL WRIT3,READ3,TEST,AKSI,LOUTS,LOUTE,ISOL,LSOR,ICAL
      DIMENSION SNSWV(NY),SNSW(NY),UNW(NY),VNW(NY),DW(NY)
C
C.....INITIALISATION OF TEMPORARILY STORED VARIABLES
C
      CALL CALCDP(KKK)
      URFU=1./URF(IU)
      URFV=1./URF(IV)
      GAM=GDS(IU)
      IJ=LI(1)+1
      UNE=U(IJ)
      VNE=V(IJ)
      VNW(1)=VNE
      UNW(1)=UNE
C-----WEST SIDE
      DO 10 J=2,NJM
      IJ=LI(1)+J
      IJM=IJ-1
      IPJ=IJ+NJ
      DXE=X(IJ)-X(IJM)
      DYE=Y(IJ)-Y(IJM)
      DXKS=XC(IPJ)-XC(IJ)
      DYKS=YC(IPJ)-YC(IJ)
      GAME=0.5*VIS(IJ)*(R(IJ)+R(IJM))/(DXKS*DYE-DXE*DYKS+SMALL)
      DW(J)=GAME*(DXE**2+DYE**2)
      USE=UNE
      UNE=U(IJ+1)*FY(IJ)+U(IJ)*(1.-FY(IJ))
      VSE=VNE
      VNE=V(IJ+1)*FY(IJ)+V(IJ)*(1.-FY(IJ))
      SNSW(J)=GAME*(DYE**2*(U(IPJ)-U(IJ))-(DXE*DXKS+2.*DYE*DYKS)*
     *         (UNE-USE)-DYE*DXE*(V(IPJ)-V(IJ))+DXE*DYKS*(VNE-VSE))
      SNSWV(J)=GAME*(DXE**2*(V(IPJ)-V(IJ))-(DYE*DYKS+2.*DXE*DXKS)*
     *         (VNE-VSE)-DXE*DYE*(U(IPJ)-U(IJ))+DXKS*DYE*(UNE-USE))
      UNW(J)=UNE
      VNW(J)=VNE
   10 CONTINUE
C-----SOUTH SIDE
      DO 100 I=2,NIM
      J=1
      IJ=LI(I)+J
      IMJ=IJ-NJ
      IJP=IJ+1
      FXE=FX(IJ)
      FXP=1.-FXE
      DXN=X(IJ)-X(IMJ)
      DYN=Y(IJ)-Y(IMJ)
      DXET=XC(IJP)-XC(IJ)
      DYET=YC(IJP)-YC(IJ)
      GAMN=0.5*VIS(IJ)*(R(IJ)+R(IMJ))/(DYET*DXN-DXET*DYN+SMALL)
      DN=GAMN*(DXN**2+DYN**2)
      UNE=U(IJ+NJ)*FXE+U(IJ)*FXP
      VNE=V(IJ+NJ)*FXE+V(IJ)*FXP
      SEWN=GAMN*(DYN**2*(U(IJP)-U(IJ))-(DXN*DXET+2.*DYN*DYET)*
     *     (UNE-UNW(J))-DXN*DYN*(V(IJP)-V(IJ))+DYET*DXN*(VNE-VNW(J)))
      SEWNV=GAMN*(DXN**2*(V(IJP)-V(IJ))-(DYN*DYET+2.*DXN*DXET)*
     *     (VNE-VNW(J))-DXN*DYN*(U(IJP)-U(IJ))+DXET*DYN*(UNE-UNW(J)))
      VNW(J)=VNE
      UNW(J)=UNE
C
C.....MAIN LOOP - ASSEMBLY OF COEFFICIENTS AND SOURCES
C
      IJS=LI(I)+2
      IJE=LI(I)+NJM
      DO 100 IJ=IJS,IJE
      J=IJ-IJS+2
      IPJ=IJ+NJ
      IMJ=IJ-NJ
      IJP=IJ+1
      IJM=IJ-1
      FXE=FX(IJ)
      FXP=1.-FXE
      FYN=FY(IJ)
      FYP=1.-FYN
C
      DXE=X(IJ)-X(IJM)
      DYE=Y(IJ)-Y(IJM)
      DXN=X(IJ)-X(IMJ)
      DYN=Y(IJ)-Y(IMJ)
      DXKS=XC(IPJ)-XC(IJ)
      DYKS=YC(IPJ)-YC(IJ)
      DXET=XC(IJP)-XC(IJ)
      DYET=YC(IJP)-YC(IJ)
      RE=0.5*(R(IJ)+R(IJM))
      RN=0.5*(R(IJ)+R(IMJ))
      GAME=(VIS(IJ)*FXP+VIS(IPJ)*FXE)*RE/(DXKS*DYE-DXE*DYKS+SMALL)
      GAMN=(VIS(IJ)*FYP+VIS(IJP)*FYN)*RN/(DYET*DXN-DYN*DXET+SMALL)
C
      DS=DN
      DE=GAME*(DXE**2+DYE**2)
      DN=GAMN*(DXN**2+DYN**2)
C
      CE=AMAX1(-F1(IJ),0.)
      CW=AMAX1(F1(IMJ),0.)
      CN=AMAX1(-F2(IJ),0.)
      CS=AMAX1(F2(IJM),0.)
      AE(IJ)=DE+CE
      AW(IJ)=DW(J)+CW
      AN(IJ)=DN+CN
      AS(IJ)=DS+CS
      AE1=-(CE+F1(IJ)*FXE)
      AW1=-(CW-F1(IMJ)*(1.-FX(IMJ)))
      AN1=-(CN+F2(IJ)*FYN)
      AS1=-(CS-F2(IJM)*(1.-FY(IJM)))
      APO=AE1+AW1+AN1+AS1
C
      DYEP=0.5*(Y(IJ)-Y(IJM)+Y(IMJ)-Y(IMJ-1))
      DYKP=0.5*(Y(IJ)-Y(IMJ)+Y(IJM)-Y(IMJ-1))
      DXKP=0.5*(X(IJ)-X(IMJ)+X(IJM)-X(IMJ-1))
      DXEP=0.5*(X(IJ)-X(IJM)+X(IMJ)-X(IMJ-1))
      RP=.25*(R(IJ)+R(IMJ)+R(IJM)+R(IMJ-1))
C
      USE=UNE
      UNE=(U(IJ)*FYP+U(IJP)*FYN)*FXP+(U(IPJ)*FYP+U(IPJ+1)*FYN)*FXE
      VSE=VNE
      VNE=(V(IJ)*FYP+V(IJP)*FYN)*FXP+(V(IPJ)*FYP+V(IPJ+1)*FYN)*FXE
      SEWS=SEWN
      SNSE=GAME*(DYE**2*(U(IPJ)-U(IJ))-(DXE*DXKS+2.*DYE*DYKS)*
     *     (UNE-USE)-DYE*DXE*(V(IPJ)-V(IJ))+DXE*DYKS*(VNE-VSE))
      SEWN=GAMN*(DYN**2*(U(IJP)-U(IJ))-(DXN*DXET+2.*DYN*DYET)*
     *     (UNE-UNW(J))-DXN*DYN*(V(IJP)-V(IJ))+DYET*DXN*(VNE-VNW(J)))
C
      SEWSV=SEWNV
      SNSEV=GAME*(DXE**2*(V(IPJ)-V(IJ))-(DYE*DYKS+2.*DXE*DXKS)*
     *        (VNE-VSE)-DXE*DYE*(U(IPJ)-U(IJ))+DXKS*DYE*(UNE-USE))
      SEWNV=GAMN*(DXN**2*(V(IJP)-V(IJ))-(DYN*DYET+2.*DXN*DXET)*
     *     (VNE-VNW(J))-DXN*DYN*(U(IJP)-U(IJ))+DXET*DYN*(UNE-UNW(J)))
C
      D1(IJ)=-DYEP*RP
      E1(IJ)= DYKP*RP
      D2(IJ)= DXEP*RP
      E2(IJ)=-DXKP*RP
C
      SP(IJ)=0.
      AP(IJ)=0.
      SUB=0.
      svb=0.
      IF(ICAL(IEN)) then
      heat=beta*den(ij)*vol(ij)*(t(ij)-tref)
      SUB=GRAVX*heat
      svb=gravy*heat
      endif
      SUP=D1(IJ)*DP1(IJ)+E1(IJ)*DP2(IJ)
      SVP=D2(IJ)*DP1(IJ)+E2(IJ)*DP2(IJ)
      SUC=GAM*(AE1*U(IPJ)+AW1*U(IMJ)+AN1*U(IJP)+AS1*U(IJM)-APO*U(IJ))
      SVC=GAM*(AE1*V(IPJ)+AW1*V(IMJ)+AN1*V(IJP)+AS1*V(IJM)-APO*V(IJ))
      IF(AKSI) THEN
      RPR=1./RP
      SP(IJ)=SP(IJ)+2.*VIS(IJ)*VOL(IJ)*RPR**2
      ENDIF
      SU(IJ)=SNSE-SNSW(J)+SEWN-SEWS+SUP+SUC+SUB
      SV(IJ)=SNSEV-SNSWV(J)+SEWNV-SEWSV+SVP+SVC+svb
C
      DW(J)=DE
      UNW(J)=UNE
      SNSW(J)=SNSE
      SNSWV(J)=SNSEV
      VNW(J)=VNE
  100 CONTINUE
C
C.....PROBLEM MODIFICATIONS - BOUNDARY CONDITIONS
C
      CALL MODUV (KKK)
C
      IF (.NOT.ISOL) RETURN
C----- ADDITIONAL SOURCES FOR FAS
      IF(LSOR) THEN
      DO 200 I=2,NIM
      IJS=LI(I)+2
      IJE=LI(I)+NJM
      DO 200 IJ=IJS,IJE
      APTU=AE(IJ)+AW(IJ)+AN(IJ)+AS(IJ)+AP(IJ)
      APTV=AE(IJ)+AW(IJ)+AN(IJ)+AS(IJ)+SP(IJ)
      SUC=AN(IJ)*U1(IJ+1)+AS(IJ)*U1(IJ-1)+AE(IJ)*U1(IJ+NJ)+
     *    AW(IJ)*U1(IJ-NJ)-APTU*U1(IJ)
      SVC=AN(IJ)*V1(IJ+1)+AS(IJ)*V1(IJ-1)+AE(IJ)*V1(IJ+NJ)+
     *    AW(IJ)*V1(IJ-NJ)-APTV*V1(IJ)
      RESU(IJ)=RESU(IJ)-SUC-SU(IJ)
      RESV(IJ)=RESV(IJ)-SVC-SV(IJ)
  200 CONTINUE
      ENDIF
C
      DO 300 I=2,NIM
      IJS=LI(I)+2
      IJE=LI(I)+NJM
      DO 300 IJ=IJS,IJE
      AP(IJ)=(AP(IJ)+AE(IJ)+AW(IJ)+AN(IJ)+AS(IJ))*URFU
      SU(IJ)=SU(IJ)+RESU(IJ)+(1.-URF(IU))*AP(IJ)*U(IJ)
      APU=1./AP(IJ)
      D1(IJ)=D1(IJ)*APU
      E1(IJ)=E1(IJ)*APU
  300 CONTINUE
C
C.....SOLVING F.D. EQUATION FOR U-VELOCITY
C
      CALL SIPSOL(U,IU)
C
      DO 400 I=2,NIM
      IJS=LI(I)+2
      IJE=LI(I)+NJM
      DO 400 IJ=IJS,IJE
      AP(IJ)=(AE(IJ)+AW(IJ)+AN(IJ)+AS(IJ)+SP(IJ))*URFV
      SU(IJ)=SV(IJ)+RESV(IJ)+(1.-URF(IV))*AP(IJ)*V(IJ)
      APV=1./AP(IJ)
      D2(IJ)=D2(IJ)*APV
      E2(IJ)=E2(IJ)*APV
  400 CONTINUE
C
C.....SOLVING F.D. EQUATION FOR V-VELOCITY
C
      CALL SIPSOL(V,IV)
C
      RETURN
      END
C
C**********************************************************
      SUBROUTINE CALCP(KKK)
C**********************************************************
      PARAMETER (ngit=3,nicv=21,njcv=8,NXA=NICV*(2**NGIT-1)+
     *           2*NGIT,NYA=NJCV*(2**NGIT-1)+2*NGIT,
     *           NX=NICV*2**(NGIT-1)+2,NY=NJCV*2**(NGIT-1)+2,
     *           NXYA=NICV*NJCV*(4**NGIT-1)/3+NICV*2*(2**NGIT-1)+
     *           NJCV*2*(2**NGIT-1)+4*NGIT,NXY=NX*NY,NPHI=5)
      COMMON/INDEX/  NI,NJ,NIM,NJM,NIJ,KGIT,IST,JST,IJST,NITP,
     *      LI(NX),LIG(NX),
     *      LSG(NGIT),LSR(NGIT),LSI(NGIT),MIT(NGIT),
     *      NIGIT(NGIT),NJGIT(NGIT),IGIT(NGIT),JGIT(NGIT),IJGIT(NGIT)
      COMMON/LOGIC/  WRIT3,READ3,TEST,AKSI,LOUTS,LOUTE,ISOL,LSOR,
     *      ICAL(NPHI)
      COMMON/CONTRO/ IMON,JMON,IJMON,IU,IV,IP,IEN,ISW,ISF,IPR,JPR,kgs,
     *      ITBS(NXA),ITBN(NXA),JTBW(NYA),JTBE(NYA),beta,gravx,gravy,
     *      NSW(NPHI),SOR(NPHI),RESOR(NPHI),PRTR(NPHI),URF(NPHI),
     *      GDS(NPHI),SNOR(NPHI),SORMAX,SLARGE,ALFA,GREAT,SMALL,TREF
      COMMON/VARF/ DENSIT,VISCOS,FLOREF,XMOMR,PRMR,PRR(NXYA),
     *      U(NXYA),V(NXYA),P(NXYA),PP(NXYA),T(NXYA),DEN(NXYA),
     *      VIS(NXYA),F1(NXYA),F2(NXYA),DP1(NXYA),DP2(NXYA)
      COMMON/VARG/ U1(NXYA),V1(NXYA),T1(NXYA),F1G(NXYA),F2G(NXYA),
     *      RESU(NXYA),RESV(NXYA),REST(NXYA)
      COMMON /GEO/ X(NXYA),Y(NXYA),XC(NXYA),YC(NXYA),FX(NXYA),FY(NXYA),
     *      VOL(NXYA),R(NXYA),RC(NXYA)
      COMMON/COEF/ AE(NXYA),AW(NXYA),AN(NXYA),AS(NXYA),AP(NXYA),
     *      SU(NXYA),SV(NXYA),SP(NXYA),D1(NXYA),D2(NXYA),E1(NXYA),
     *      E2(NXYA)
      COMMON/BCOEF/ BW(NXYA),BE(NXYA),BS(NXYA),BN(NXYA),BP(NXYA),
     *      RES(NXYA)
      LOGICAL WRIT3,READ3,TEST,AKSI,LOUTS,LOUTE,ISOL,LSOR,ICAL
C
C.....ASSEMBLY OF COEFFICIENTS AND SOURCES
C
      SUM=0.
      DO 100 I=2,NIM
      IJS=LI(I)+2
      IJE=LI(I)+NJM
      DO 100 IJ=IJS,IJE
      IPJ=IJ+NJ
      IJP=IJ+1
      IJM=IJ-1
      IMJ=IJ-NJ
      FXE=FX(IJ)
      FXP=1.-FXE
      FYN=FY(IJ)
      FYP=1.-FYN
C
      DXE=X(IJ)-X(IJM)
      DYE=Y(IJ)-Y(IJM)
      DXN=X(IJ)-X(IMJ)
      DYN=Y(IJ)-Y(IMJ)
      RE=0.5*(R(IJ)+R(IJM))
      RN=0.5*(R(IJ)+R(IMJ))
      D1E=D1(IJ)*FXP+D1(IPJ)*FXE
      D2E=D2(IJ)*FXP+D2(IPJ)*FXE
      E2N=E2(IJ)*FYP+E2(IJP)*FYN
      E1N=E1(IJ)*FYP+E1(IJP)*FYN
      DENE=DEN(IJ)*FXP+DEN(IPJ)*FXE
      DENN=DEN(IJ)*FYP+DEN(IJP)*FYN
C
      AE(IJ)=-DENE*(D1E*DYE-D2E*DXE)*RE
      AW(IJ)=AE(IMJ)
      AN(IJ)=-DENE*(E2N*DXN-E1N*DYN)*RN
      AS(IJ)=AN(IJM)
      DPE=P(IPJ)-P(IJ)
      DPN=P(IJP)-P(IJ)
      UE=(U(IPJ)-U1(IPJ)-D1(IPJ)*DP1(IPJ))*FXE+
     *   (U(IJ)- U1(IJ)- D1(IJ) *DP1(IJ)) *FXP+D1E*DPE
      VE=(V(IPJ)-V1(IPJ)-D2(IPJ)*DP1(IPJ))*FXE+
     *   (V(IJ)- V1(IJ)- D2(IJ) *DP1(IJ)) *FXP+D2E*DPE
      UN=(U(IJP)-U1(IJP)-E1(IJP)*DP2(IJP))*FYN+
     *   (U(IJ)- U1(IJ)- E1(IJ) *DP2(IJ)) *FYP+E1N*DPN
      VN=(V(IJP)-V1(IJP)-E2(IJP)*DP2(IJP))*FYN+
     *   (V(IJ)- V1(IJ)- E2(IJ) *DP2(IJ)) *FYP+E2N*DPN
C
      F1(IJ)=(UE*DYE-VE*DXE)*DENE*RE+F1G(IJ)
      F2(IJ)=(VN*DXN-UN*DYN)*DENN*RN+F2G(IJ)
C
      SU(IJ)=F1(IMJ)-F1(IJ)+F2(IJM)-F2(IJ)
      SUM=SUM+SU(IJ)
      SP(IJ)=0.0
  100 CONTINUE
      IF(TEST) WRITE(6,*) '       SUM = ',SUM
C
C.....PROBLEM MODIFICATIONS
C
      CALL MODPR (KKK)
      DO 200 I=2,NIM
      IJS=LI(I)+2
      IJE=LI(I)+NJM
      DO 200 IJ=IJS,IJE
      PP(IJ)=0.
      AP(IJ)=AE(IJ)+AW(IJ)+AN(IJ)+AS(IJ)+SP(IJ)
  200 CONTINUE
C
C.....CORRECTOR STEPS
C
      CALL SIPSOL(PP,IP)
C
C.....CORRECT FLUXES,VELOCITY COMPONENTS AND PRESSURE
C
      DO 500 L=1,NITP
C
      CALL CALCDP(KKK)
      IJPREF=LI(IPR)+JPR
      PPREF=PP(IJPREF)
C
      DO 300 I=2,NIM
      IJS=LI(I)+2
      IJE=LI(I)+NJM
      DO 300 IJ=IJS,IJE
      IPJ=IJ+NJ
      IJP=IJ+1
      F1(IJ)=F1(IJ)-AE(IJ)*(PP(IPJ)-PP(IJ))
      F2(IJ)=F2(IJ)-AN(IJ)*(PP(IJP)-PP(IJ))
      U(IJ)=U(IJ)+D1(IJ)*DP1(IJ)+E1(IJ)*DP2(IJ)
      V(IJ)=V(IJ)+D2(IJ)*DP1(IJ)+E2(IJ)*DP2(IJ)
      P(IJ)=P(IJ)+URF(IP)*(PP(IJ)-PPREF)
  300 CONTINUE
C
C.....SECOND FIRST CORRECTOR STEP
C
      IF(L.LT.NITP) THEN
      DO 400 I=2,NIM
      F2N=0.
      IJS=LI(I)+2
      IJE=LI(I)+NJM
      DO 400 IJ=IJS,IJE
      IPJ=IJ+NJ
      IJP=IJ+1
C-------EAST SIDE
      DYE=Y(IJ)-Y(IJ-1)
      DXE=X(IJ)-X(IJ-1)
      RE=0.5*(R(IJ)+R(IJ-1))
      DP2P=DP2(IJ)*(1.-FX(IJ))
      DP2E=DP2(IPJ)*FX(IJ)
      UE=E1(IJ)*DP2P+E1(IPJ)*DP2E
      VE=E2(IJ)*DP2P+E2(IPJ)*DP2E
      F1E=(UE*DYE-VE*DXE)*DENSIT*RE
      F1(IJ)=F1(IJ)+F1E
C-------NORTH SIDE
      DXN=X(IJ)-X(IJ-NJ)
      DYN=Y(IJ)-Y(IJ-NJ)
      RN=0.5*(R(IJ)+R(IJ-NJ))
      DP1P=DP1(IJ)*(1.-FY(IJ))
      DP1N=DP1(IJP)*FY(IJ)
      UN=D1(IJ)*DP1P+D1(IJP)*DP1N
      VN=D2(IJ)*DP1P+D2(IJP)*DP1N
      F2N=(VN*DXN-UN*DYN)*DENSIT*RN
      F2(IJ)=F2(IJ)+F2N
C-------SOURCE TERM
      SU(IJ)=F1(IJ-NJ)-F1(IJ)+F2(IJ-1)-F2(IJ)
      SUM=SUM+SU(IJ)
      PP(IJ)=0.
  400 CONTINUE
      IF(TEST) WRITE(6,*) '           SUMC =',SUM
C-------END
      CALL SOLVP
      ENDIF
  500 CONTINUE
      DO 550 IJ=IJST+1,IJST+NIJ
  550 PP(IJ)=P(IJ)
      RETURN
      END
C
C**********************************************************
      SUBROUTINE CALCT(KKK)
C**********************************************************
      PARAMETER (ngit=3,nicv=21,njcv=8,NXA=NICV*(2**NGIT-1)+
     *           2*NGIT,NYA=NJCV*(2**NGIT-1)+2*NGIT,
     *           NX=NICV*2**(NGIT-1)+2,NY=NJCV*2**(NGIT-1)+2,
     *           NXYA=NICV*NJCV*(4**NGIT-1)/3+NICV*2*(2**NGIT-1)+
     *           NJCV*2*(2**NGIT-1)+4*NGIT,NXY=NX*NY,NPHI=5)
      COMMON/INDEX/  NI,NJ,NIM,NJM,NIJ,KGIT,IST,JST,IJST,NITP,
     *      LI(NX),LIG(NX),
     *      LSG(NGIT),LSR(NGIT),LSI(NGIT),MIT(NGIT),
     *      NIGIT(NGIT),NJGIT(NGIT),IGIT(NGIT),JGIT(NGIT),IJGIT(NGIT)
      COMMON/LOGIC/  WRIT3,READ3,TEST,AKSI,LOUTS,LOUTE,ISOL,LSOR,
     *      ICAL(NPHI)
      COMMON/CONTRO/ IMON,JMON,IJMON,IU,IV,IP,IEN,ISW,ISF,IPR,JPR,kgs,
     *      ITBS(NXA),ITBN(NXA),JTBW(NYA),JTBE(NYA),beta,gravx,gravy,
     *      NSW(NPHI),SOR(NPHI),RESOR(NPHI),PRTR(NPHI),URF(NPHI),
     *      GDS(NPHI),SNOR(NPHI),SORMAX,SLARGE,ALFA,GREAT,SMALL,TREF
      COMMON/VARF/ DENSIT,VISCOS,FLOREF,XMOMR,PRMR,PRR(NXYA),
     *      U(NXYA),V(NXYA),P(NXYA),PP(NXYA),T(NXYA),DEN(NXYA),
     *      VIS(NXYA),F1(NXYA),F2(NXYA),DP1(NXYA),DP2(NXYA)
      COMMON/VARG/ U1(NXYA),V1(NXYA),T1(NXYA),F1G(NXYA),F2G(NXYA),
     *      RESU(NXYA),RESV(NXYA),REST(NXYA)
      COMMON /GEO/ X(NXYA),Y(NXYA),XC(NXYA),YC(NXYA),FX(NXYA),FY(NXYA),
     *      VOL(NXYA),R(NXYA),RC(NXYA)
      COMMON/COEF/ AE(NXYA),AW(NXYA),AN(NXYA),AS(NXYA),AP(NXYA),
     *      SU(NXYA),SV(NXYA),SP(NXYA),D1(NXYA),D2(NXYA),E1(NXYA),
     *      E2(NXYA)
      COMMON/BCOEF/ BW(NXYA),BE(NXYA),BS(NXYA),BN(NXYA),BP(NXYA),
     *      RES(NXYA)
      LOGICAL WRIT3,READ3,TEST,AKSI,LOUTS,LOUTE,ISOL,LSOR,ICAL
      DIMENSION SNSW(NY),DW(NY),FINW(NY)
C
C.....INITIALISATION OF TEMPORARILY STORED VARIABLES
C
      URFI=1./URF(IEN)
      GAM=GDS(IEN)
      IJ=LI(1)+1
      FINE=T(IJ)
      FINW(1)=FINE
C-----WEST SIDE
      DO 10 J=2,NJM
      IJ=LI(1)+J
      IJM=IJ-1
      IPJ=IJ+NJ
      DXE=X(IJ)-X(IJM)
      DYE=Y(IJ)-Y(IJM)
      DXKS=XC(IPJ)-XC(IJ)
      DYKS=YC(IPJ)-YC(IJ)
      RP=0.5*(R(IJ)+R(IJM))
      GAME=VIS(IJ)*PRR(IJ)*RP/(DXKS*DYE-DXE*DYKS+SMALL)
      DW(J)=GAME*(DXE**2+DYE**2)
      FISE=FINE
      FINE=T(IJ+1)*FY(IJ)+T(IJ)*(1.-FY(IJ))
      SNSW(J)=-GAME*(DXKS*DXE+DYKS*DYE)*(FINE-FISE)
      FINW(J)=FINE
   10 CONTINUE
C-----SOUTH SIDE
      DO 100 I=2,NIM
      J=1
      IJ=LI(I)+J
      IMJ=IJ-NJ
      IJP=IJ+1
      FXE=FX(IJ)
      FXP=1.-FXE
      DXN=X(IJ)-X(IMJ)
      DYN=Y(IJ)-Y(IMJ)
      DXET=XC(IJP)-XC(IJ)
      DYET=YC(IJP)-YC(IJ)
      RP=0.5*(R(IJ)+R(IMJ))
      GAMN=VIS(IJ)*PRR(IJ)*RP/(DYET*DXN-DXET*DYN+SMALL)
      DN=GAMN*(DXN**2+DYN**2)
      FINE=T(IJ+NJ)*FXE+T(IJ)*FXP
      SEWN=-GAMN*(DXN*DXET+DYN*DYET)*(FINE-FINW(J))
      FINW(J)=FINE
C
C.....MAIN LOOP - ASSEMBLY OF COEFFICIENTS AND SOURCES
C
      ijs=li(i)+2
      ije=li(i)+njm
      do 100 ij=ijs,ije
      j=ij-ijs+2
      IPJ=IJ+NJ
      IMJ=IJ-NJ
      IJP=IJ+1
      IJM=IJ-1
      FXE=FX(IJ)
      FXP=1.-FXE
      FYN=FY(IJ)
      FYP=1.-FYN
C
      DXE=X(IJ)-X(IJM)
      DYE=Y(IJ)-Y(IJM)
      DXN=X(IJ)-X(IMJ)
      DYN=Y(IJ)-Y(IMJ)
      DXKS=XC(IPJ)-XC(IJ)
      DYKS=YC(IPJ)-YC(IJ)
      DXET=XC(IJP)-XC(IJ)
      DYET=YC(IJP)-YC(IJ)
      RE=0.5*(R(IJ)+R(IJM))
      RN=0.5*(R(IJ)+R(IMJ))
      VISE=VIS(IJ)*PRR(IJ)*FXP+VIS(IPJ)*PRR(IJ)*FXE
      VISN=VIS(IJ)*PRR(IJ)*FYP+VIS(IJP)*PRR(IJ)*FYN
      GAME=VISE*RE/(DXKS*DYE-DXE*DYKS+SMALL)
      GAMN=VISN*RN/(DYET*DXN-DYN*DXET+SMALL)
C
      DS=DN
      DE=GAME*(DXE**2+DYE**2)
      DN=GAMN*(DXN**2+DYN**2)
C
      CE=AMAX1(-F1(IJ),0.)
      CW=AMAX1(F1(IMJ),0.)
      CN=AMAX1(-F2(IJ),0.)
      CS=AMAX1(F2(IJM),0.)
      AE(IJ)=DE+CE
      AW(IJ)=DW(J)+CW
      AN(IJ)=DN+CN
      AS(IJ)=DS+CS
      AE1=-(CE+F1(IJ)*FXE)
      AW1=-(CW-F1(IMJ)*(1.-FX(IMJ)))
      AN1=-(CN+F2(IJ)*FYN)
      AS1=-(CS-F2(IJM)*(1.-FY(IJM)))
      APO=AE1+AW1+AN1+AS1
C
      FISE=FINE
      FINE=(T(IJ)*FYP+T(IJP)*FYN)*FXP+
     *     (T(IPJ)*FYP+T(IPJ+1)*FYN)*FXE
      SEWS=SEWN
      SEWN=-GAMN*(DXN*DXET+DYN*DYET)*(FINE-FINW(J))
      SNSE=-GAME*(DXKS*DXE+DYKS*DYE)*(FINE-FISE)
      SUC=GAM*(AE1*T(IPJ)+AW1*T(IMJ)+AN1*T(IJP)+AS1*T(IJM)-
     *        APO*T(IJ))
      SU(IJ)=SNSE-SNSW(J)+SEWN-SEWS+SUC
      SP(IJ)=0.
      SNSW(J)=SNSE
      FINW(J)=FINE
      DW(J)=DE
  100 CONTINUE
C
      CALL MODT(KKK)
      if(.not.isol) return
C----- ADDITIONAL SOURCES FOR FAS
      IF(LSOR) THEN
      DO 200 I=2,NIM
      IJS=LI(I)+2
      IJE=LI(I)+NJM
      DO 200 IJ=IJS,IJE
      APT=AE(IJ)+AW(IJ)+AN(IJ)+AS(IJ)+SP(IJ)
      SUC=AN(IJ)*T1(IJ+1)+AS(IJ)*T1(IJ-1)+AE(IJ)*T1(IJ+NJ)+
     *    AW(IJ)*T1(IJ-NJ)-APT*T1(IJ)
      REST(IJ)=REST(IJ)-SUC-SU(IJ)
  200 CONTINUE
      ENDIF
C
C.....PROBLEM MODIFICATIONS - BOUNDARY CONDITIONS
C
      DO 250 I=2,NIM
      DO 250 J=2,NJM
      IJ=LI(I)+J
      AP(IJ)=(AW(IJ)+AE(IJ)+AN(IJ)+AS(IJ)+SP(IJ))*URFI
      SU(IJ)=SU(IJ)+REST(IJ)+(1.-URF(IEN))*AP(IJ)*T(IJ)
  250 CONTINUE
C
C.....SOLVING F.D. EQUATIONS
C
      CALL SIPSOL(T,IEN)
      RETURN
      END
C
C
C************************************************************
      SUBROUTINE PRINT(PHI,HEDPHI)
C************************************************************
      PARAMETER (ngit=3,nicv=21,njcv=8,NXA=NICV*(2**NGIT-1)+
     *           2*NGIT,NYA=NJCV*(2**NGIT-1)+2*NGIT,
     *           NX=NICV*2**(NGIT-1)+2,NY=NJCV*2**(NGIT-1)+2,
     *           NXYA=NICV*NJCV*(4**NGIT-1)/3+NICV*2*(2**NGIT-1)+
     *           NJCV*2*(2**NGIT-1)+4*NGIT,NXY=NX*NY,NPHI=5)
      COMMON/INDEX/  NI,NJ,NIM,NJM,NIJ,KGIT,IST,JST,IJST,NITP,
     *      LI(NX),LIG(NX),
     *      LSG(NGIT),LSR(NGIT),LSI(NGIT),MIT(NGIT),
     *      NIGIT(NGIT),NJGIT(NGIT),IGIT(NGIT),JGIT(NGIT),IJGIT(NGIT)
      DIMENSION PHI(NXYA)
      CHARACTER*6 HEDPHI
      WRITE(6,20) HEDPHI
      IS=-11
  100 IS=IS+12
      IE=IS+11
      IE=MIN0(NI,IE)
      WRITE(6,21) (I,I=IS,IE)
      WRITE(6,22)
      DO 101 JJ=1,NJ
      J=NJ+1-JJ
  101 WRITE(6,23) J,(PHI(LI(I)+J),I=IS,IE)
      IF(IE.LT.NI) GO TO 100
   20 FORMAT(2X,26('*-'),7X,A6,7X,26('-*'))
   21 FORMAT(3X,'I = ',I3,11I10)
   22 FORMAT(2X,'J')
   23 FORMAT(1X,I3,1P12E10.2)
      RETURN
      END
C
C********************************************************
      SUBROUTINE SIPSOL(FI,IFI)
C********************************************************
      PARAMETER (ngit=3,nicv=21,njcv=8,NXA=NICV*(2**NGIT-1)+
     *           2*NGIT,NYA=NJCV*(2**NGIT-1)+2*NGIT,
     *           NX=NICV*2**(NGIT-1)+2,NY=NJCV*2**(NGIT-1)+2,
     *           NXYA=NICV*NJCV*(4**NGIT-1)/3+NICV*2*(2**NGIT-1)+
     *           NJCV*2*(2**NGIT-1)+4*NGIT,NXY=NX*NY,NPHI=5)
      COMMON/INDEX/  NI,NJ,NIM,NJM,NIJ,KGIT,IST,JST,IJST,NITP,
     *      LI(NX),LIG(NX),
     *      LSG(NGIT),LSR(NGIT),LSI(NGIT),MIT(NGIT),
     *      NIGIT(NGIT),NJGIT(NGIT),IGIT(NGIT),JGIT(NGIT),IJGIT(NGIT)
      COMMON/LOGIC/  WRIT3,READ3,TEST,AKSI,LOUTS,LOUTE,ISOL,LSOR,
     *      ICAL(NPHI)
      COMMON/CONTRO/ IMON,JMON,IJMON,IU,IV,IP,IEN,ISW,ISF,IPR,JPR,kgs,
     *      ITBS(NXA),ITBN(NXA),JTBW(NYA),JTBE(NYA),beta,gravx,gravy,
     *      NSW(NPHI),SOR(NPHI),RESOR(NPHI),PRTR(NPHI),URF(NPHI),
     *      GDS(NPHI),SNOR(NPHI),SORMAX,SLARGE,ALFA,GREAT,SMALL,TREF
      COMMON/COEF/ AE(NXYA),AW(NXYA),AN(NXYA),AS(NXYA),AP(NXYA),
     *      SU(NXYA),SV(NXYA),SP(NXYA),D1(NXYA),D2(NXYA),E1(NXYA),
     *      E2(NXYA)
      COMMON/BCOEF/ BW(NXYA),BE(NXYA),BS(NXYA),BN(NXYA),BP(NXYA),
     *      RES(NXYA)
      LOGICAL WRIT3,READ3,TEST,AKSI,LOUTS,LOUTE,ISOL,LSOR,ICAL
      DIMENSION FI(NXYA)
C
      NS=NSW(IFI)
      DO 20 I=2,NIM
      IJS=LI(I)+2
      IJE=LI(I)+NJM
      DO 20 IJ=IJS,IJE
      IJM=IJ-1
      IMJ=IJ-NJ
      BW(IJ)=-AW(IJ)/(1.+ALFA*BN(IMJ))
      BS(IJ)=-AS(IJ)/(1.+ALFA*BE(IJM))
      P1=ALFA*BW(IJ)*BN(IMJ)
      P2=ALFA*BS(IJ)*BE(IJM)
      BP(IJ)=1./(AP(IJ)+P1+P2-BW(IJ)*BE(IMJ)-BS(IJ)*BN(IJM)+SMALL)
      BN(IJ)=(-AN(IJ)-P1)*BP(IJ)
      BE(IJ)=(-AE(IJ)-P2)*BP(IJ)
   20 CONTINUE
      DO 100 L=1,NS
      RESAB=0.
      DO 30 I=2,NIM
      IJS=LI(I)+2
      IJE=LI(I)+NJM
      DO 30 IJ=IJS,IJE
      RES(IJ)=AN(IJ)*FI(IJ+1)+AS(IJ)*FI(IJ-1)+AE(IJ)*FI(IJ+NJ)+
     *        AW(IJ)*FI(IJ-NJ)+SU(IJ)-AP(IJ)*FI(IJ)
      RESAB=RESAB+ABS(RES(IJ))
      RES(IJ)=(RES(IJ)-BS(IJ)*RES(IJ-1)-BW(IJ)*RES(IJ-NJ))*BP(IJ)
   30 CONTINUE
      IF(L.EQ.1) RESOR(IFI)=RESAB
      RSM=RESAB/(RESOR(IFI)+SMALL)
      DO 40 I=NIM,2,-1
      IJS=LI(I)+NJM
      IJE=LI(I)+2
      DO 40 IJ=IJS,IJE,-1
   40 RES(IJ)=RES(IJ)-BN(IJ)*RES(IJ+1)-BE(IJ)*RES(IJ+NJ)
      DO 50 I=2,NIM
      IJS=LI(I)+2
      IJE=LI(I)+NJM
      DO 50 IJ=IJS,IJE
   50 FI(IJ)=FI(IJ)+RES(IJ)
      IF(TEST) WRITE(6,*) '    ',L,'SWEEP, RESOR = ',RESAB
      IF(RSM.LT.SOR(IFI)) RETURN
  100 CONTINUE
      RETURN
      END
C
C******************************************************
      SUBROUTINE SOLVP
C******************************************************
      PARAMETER (ngit=3,nicv=21,njcv=8,NXA=NICV*(2**NGIT-1)+
     *           2*NGIT,NYA=NJCV*(2**NGIT-1)+2*NGIT,
     *           NX=NICV*2**(NGIT-1)+2,NY=NJCV*2**(NGIT-1)+2,
     *           NXYA=NICV*NJCV*(4**NGIT-1)/3+NICV*2*(2**NGIT-1)+
     *           NJCV*2*(2**NGIT-1)+4*NGIT,NXY=NX*NY,NPHI=5)
      COMMON/INDEX/  NI,NJ,NIM,NJM,NIJ,KGIT,IST,JST,IJST,NITP,
     *      LI(NX),LIG(NX),
     *      LSG(NGIT),LSR(NGIT),LSI(NGIT),MIT(NGIT),
     *      NIGIT(NGIT),NJGIT(NGIT),IGIT(NGIT),JGIT(NGIT),IJGIT(NGIT)
      COMMON/LOGIC/  WRIT3,READ3,TEST,AKSI,LOUTS,LOUTE,ISOL,LSOR,
     *      ICAL(NPHI)
      COMMON/CONTRO/ IMON,JMON,IJMON,IU,IV,IP,IEN,ISW,ISF,IPR,JPR,kgs,
     *      ITBS(NXA),ITBN(NXA),JTBW(NYA),JTBE(NYA),beta,gravx,gravy,
     *      NSW(NPHI),SOR(NPHI),RESOR(NPHI),PRTR(NPHI),URF(NPHI),
     *      GDS(NPHI),SNOR(NPHI),SORMAX,SLARGE,ALFA,GREAT,SMALL,TREF
      COMMON/VARF/ DENSIT,VISCOS,FLOREF,XMOMR,PRMR,PRR(NXYA),
     *      U(NXYA),V(NXYA),P(NXYA),PP(NXYA),T(NXYA),DEN(NXYA),
     *      VIS(NXYA),F1(NXYA),F2(NXYA),DP1(NXYA),DP2(NXYA)
      COMMON/COEF/ AE(NXYA),AW(NXYA),AN(NXYA),AS(NXYA),AP(NXYA),
     *      SU(NXYA),SV(NXYA),SP(NXYA),D1(NXYA),D2(NXYA),E1(NXYA),
     *      E2(NXYA)
      COMMON/BCOEF/ BW(NXYA),BE(NXYA),BS(NXYA),BN(NXYA),BP(NXYA),
     *      RES(NXYA)
      LOGICAL WRIT3,READ3,TEST,AKSI,LOUTS,LOUTE,ISOL,LSOR,ICAL
      NS=NSW(IP)
      DO 100 L=1,NS
      RESAB=0.
      DO 30 I=2,NIM
      IJS=LI(I)+2
      IJE=LI(I)+NJM
      DO 30 IJ=IJS,IJE
      RES(IJ)=AN(IJ)*PP(IJ+1)+AS(IJ)*PP(IJ-1)+AE(IJ)*PP(IJ+NJ)+
     *        AW(IJ)*PP(IJ-NJ)+SU(IJ)-AP(IJ)*PP(IJ)
      RESAB=RESAB+ABS(RES(IJ))
      RES(IJ)=(RES(IJ)-BS(IJ)*RES(IJ-1)-BW(IJ)*RES(IJ-NJ))*BP(IJ)
   30 CONTINUE
      RSM=RESAB/(RESOR(IP)+SMALL)
      DO 40 I=NIM,2,-1
      IJS=LI(I)+NJM
      IJE=LI(I)+2
      DO 40 IJ=IJS,IJE,-1
   40 RES(IJ)=RES(IJ)-BN(IJ)*RES(IJ+1)-BE(IJ)*RES(IJ+NJ)
      DO 50 I=2,NIM
      IJS=LI(I)+2
      IJE=LI(I)+NJM
      DO 50 IJ=IJS,IJE
   50 PP(IJ)=PP(IJ)+RES(IJ)
      IF(TEST) WRITE(6,*) '    ',L,' SWEEP, RESOR = ',RESAB
      IF(RSM.LT.SOR(IP)) RETURN
  100 CONTINUE
      RETURN
      END
C
C**********************************************************
      SUBROUTINE WALLBC(IJP,IJB,IJL,IJR,IJC,IFI)
C**********************************************************
      PARAMETER (ngit=3,nicv=21,njcv=8,NXA=NICV*(2**NGIT-1)+
     *           2*NGIT,NYA=NJCV*(2**NGIT-1)+2*NGIT,
     *           NX=NICV*2**(NGIT-1)+2,NY=NJCV*2**(NGIT-1)+2,
     *           NXYA=NICV*NJCV*(4**NGIT-1)/3+NICV*2*(2**NGIT-1)+
     *           NJCV*2*(2**NGIT-1)+4*NGIT,NXY=NX*NY,NPHI=5)
      COMMON /GEO/ X(NXYA),Y(NXYA),XC(NXYA),YC(NXYA),FX(NXYA),FY(NXYA),
     *      VOL(NXYA),R(NXYA),RC(NXYA)
      COMMON/VARF/ DENSIT,VISCOS,FLOREF,XMOMR,PRMR,PRR(NXYA),
     *      U(NXYA),V(NXYA),P(NXYA),PP(NXYA),T(NXYA),DEN(NXYA),
     *      VIS(NXYA),F1(NXYA),F2(NXYA),DP1(NXYA),DP2(NXYA)
      COMMON /WAAL/ TAU,SUU,SUP,SVU,SVP,QWAL
C
C.....Define unit vectors  p  and  n
C
      DXB=X(IJC)-XC(IJB)
      DYB=Y(IJC)-YC(IJB)
      RBC=1./SQRT(DXB**2+DYB**2)
      ALF=DXB*RBC
      BET=DYB*RBC
      ARB=2./RBC*RC(IJB)
C
C.....Determine neighbour and interp. factors, deln
C
      BPP=(XC(IJP)-XC(IJB))*ALF+(YC(IJP)-YC(IJB))*BET
      IJN=IJR
      IF(BPP.GT.0.) IJN=IJL
      FAC=BPP/((XC(IJP)-XC(IJN))*ALF+(YC(IJP)-YC(IJN))*BET)
      DNR=1./((YC(IJP)-YC(IJB))*ALF-(XC(IJP)-XC(IJB))*BET-
     *    FAC*((YC(IJP)-YC(IJN))*ALF-(XC(IJP)-XC(IJN))*BET))
      FACP=1.-FAC
C
C.....Calculate wall shear stress & force components
C
      IF (IFI.EQ.1) THEN
      COEF=VISCOS*ARB*DNR
      Q=(U(IJN)*ALF+V(IJN)*BET)*FAC-U(IJB)*ALF-V(IJB)*BET
      TAU=-VISCOS*DNR*((U(IJP)*ALF+V(IJP)*BET)*FACP+Q)
      SUU=-COEF*ALF*(V(IJP)*BET*FACP+Q)
      SUP=COEF*ALF**2*FACP
      SVU=-COEF*BET*(U(IJP)*ALF*FACP+Q)
      SVP=COEF*BET**2*FACP
      ELSEIF (IFI.EQ.4) THEN
      COEF=VISCOS*PRMR*ARB*DNR
      SUP=COEF*FACP
      SUU=COEF*(T(IJB)-T(IJN)*FAC)
      qwal=suu-sup*t(ijp)
      ENDIF
      RETURN
      END
C
C*************************************************************
      SUBROUTINE BC1(IJP,IJB,IJL,IJR,IJC,IBC)
C*************************************************************
      PARAMETER (ngit=3,nicv=21,njcv=8,NXA=NICV*(2**NGIT-1)+
     *           2*NGIT,NYA=NJCV*(2**NGIT-1)+2*NGIT,
     *           NX=NICV*2**(NGIT-1)+2,NY=NJCV*2**(NGIT-1)+2,
     *           NXYA=NICV*NJCV*(4**NGIT-1)/3+NICV*2*(2**NGIT-1)+
     *           NJCV*2*(2**NGIT-1)+4*NGIT,NXY=NX*NY,NPHI=5)
      COMMON/INDEX/  NI,NJ,NIM,NJM,NIJ,KGIT,IST,JST,IJST,NITP,
     *      LI(NX),LIG(NX),
     *      LSG(NGIT),LSR(NGIT),LSI(NGIT),MIT(NGIT),
     *      NIGIT(NGIT),NJGIT(NGIT),IGIT(NGIT),JGIT(NGIT),IJGIT(NGIT)
      COMMON/LOGIC/  WRIT3,READ3,TEST,AKSI,LOUTS,LOUTE,ISOL,LSOR,
     *      ICAL(NPHI)
      COMMON/VARF/ DENSIT,VISCOS,FLOREF,XMOMR,PRMR,PRR(NXYA),
     *      U(NXYA),V(NXYA),P(NXYA),PP(NXYA),T(NXYA),DEN(NXYA),
     *      VIS(NXYA),F1(NXYA),F2(NXYA),DP1(NXYA),DP2(NXYA)
      COMMON/VARG/ U1(NXYA),V1(NXYA),T1(NXYA),F1G(NXYA),F2G(NXYA),
     *      RESU(NXYA),RESV(NXYA),REST(NXYA)
      COMMON /GEO/ X(NXYA),Y(NXYA),XC(NXYA),YC(NXYA),FX(NXYA),FY(NXYA),
     *      VOL(NXYA),R(NXYA),RC(NXYA)
      COMMON /WAAL/ TAU,SUU,SUP,SVU,SVP,QWAL
      COMMON /SYMB/ FIB,FIP,FIL,FIR,FAC,IJN
      logical writ3,read3,test,aksi,louts,loute,isol,lsor,ical
C
      if(ibc.eq.3) then
      CALL SYMV1 (IJP,IJB,IJL,IJR,IJC)
      endif
      IF(ical(4)) then
      FIP=T1(IJP)
      FIL=T1(IJL)
      FIR=T1(IJR)
      CALL SYMBC (IJP,IJB,IJL,IJR,IJC)
      T1(IJB)=FIB
      ENDIF
      RETURN
      END
C
C**********************************************************
      SUBROUTINE SYMBC(IJP,IJB,IJL,IJR,IJC)
C**********************************************************
      PARAMETER (ngit=3,nicv=21,njcv=8,NXA=NICV*(2**NGIT-1)+
     *           2*NGIT,NYA=NJCV*(2**NGIT-1)+2*NGIT,
     *           NX=NICV*2**(NGIT-1)+2,NY=NJCV*2**(NGIT-1)+2,
     *           NXYA=NICV*NJCV*(4**NGIT-1)/3+NICV*2*(2**NGIT-1)+
     *           NJCV*2*(2**NGIT-1)+4*NGIT,NXY=NX*NY,NPHI=5)
      COMMON /GEO/ X(NXYA),Y(NXYA),XC(NXYA),YC(NXYA),FX(NXYA),FY(NXYA),
     *      VOL(NXYA),R(NXYA),RC(NXYA)
      COMMON /SYMB/ FIB,FIP,FIL,FIR,FAC,IJN
C
      DXB=X(IJC)-XC(IJB)
      DYB=Y(IJC)-YC(IJB)
      PBB=(XC(IJP)-XC(IJB))*DXB+(YC(IJP)-YC(IJB))*DYB
      IF(PBB.GT.0.) THEN
      FAC=PBB/((XC(IJP)-XC(IJL))*DXB+(YC(IJP)-YC(IJL))*DYB)
      FIB=FAC*FIL+(1.-FAC)*FIP
      IJN=IJL
      ELSE
      FAC=PBB/((XC(IJP)-XC(IJR))*DXB+(YC(IJP)-YC(IJR))*DYB)
      FIB=FAC*FIR+(1.-FAC)*FIP
      IJN=IJR
      ENDIF
      RETURN
      END
C
C*************************************************************
      SUBROUTINE SYMVEL(IJP,IJB,IJL,IJR,IJC)
C*************************************************************
      PARAMETER (ngit=3,nicv=21,njcv=8,NXA=NICV*(2**NGIT-1)+
     *           2*NGIT,NYA=NJCV*(2**NGIT-1)+2*NGIT,
     *           NX=NICV*2**(NGIT-1)+2,NY=NJCV*2**(NGIT-1)+2,
     *           NXYA=NICV*NJCV*(4**NGIT-1)/3+NICV*2*(2**NGIT-1)+
     *           NJCV*2*(2**NGIT-1)+4*NGIT,NXY=NX*NY,NPHI=5)
      COMMON /GEO/ X(NXYA),Y(NXYA),XC(NXYA),YC(NXYA),FX(NXYA),FY(NXYA),
     *      VOL(NXYA),R(NXYA),RC(NXYA)
      COMMON/VARF/ DENSIT,VISCOS,FLOREF,XMOMR,PRMR,PRR(NXYA),
     *      U(NXYA),V(NXYA),P(NXYA),PP(NXYA),T(NXYA),DEN(NXYA),
     *      VIS(NXYA),F1(NXYA),F2(NXYA),DP1(NXYA),DP2(NXYA)
      COMMON /WAAL/ TAU,SUU,SUP,SVU,SVP,QWAL
      COMMON /SYMB/ FIB,FIP,FIL,FIR,FAC,IJN
C
C.....Define unit vectors  p  and  n
C
      DXB=X(IJC)-XC(IJB)
      DYB=Y(IJC)-YC(IJB)
      RBC=1./SQRT(DXB**2+DYB**2)
      ALF=DXB*RBC
      BET=DYB*RBC
      ARB=2./RBC*RC(IJB)
C
C.....Calculate velocities at symmetry
C
      FIP=U(IJP)*ALF+V(IJP)*BET
      FIL=U(IJL)*ALF+V(IJL)*BET
      FIR=U(IJR)*ALF+V(IJR)*BET
      CALL SYMBC(IJP,IJB,IJL,IJR,IJC)
      U(IJB)=FIB*ALF
      V(IJB)=FIB*BET
C
C.....Calculate normal stress and force
C
      FACP=1.-FAC
      DNR=1./((YC(IJP)-YC(IJB))*ALF-(XC(IJP)-XC(IJB))*BET-
     *    FAC*((YC(IJP)-YC(IJN))*ALF-(XC(IJP)-XC(IJN))*BET))
      COEF=2.*VISCOS*ARB*DNR
      Q=(V(IJN)*ALF-U(IJN)*BET)*FAC
      SUU=COEF*BET*(V(IJP)*ALF*FACP+Q)
      SUP=COEF*BET**2*FACP
      SVU=COEF*ALF*(U(IJP)*BET*FACP-Q)
      SVP=COEF*ALF**2*FACP
      RETURN
      END
C
C*************************************************************
      SUBROUTINE SYMV1(IJP,IJB,IJL,IJR,IJC)
C*************************************************************
      PARAMETER (ngit=3,nicv=21,njcv=8,NXA=NICV*(2**NGIT-1)+
     *           2*NGIT,NYA=NJCV*(2**NGIT-1)+2*NGIT,
     *           NX=NICV*2**(NGIT-1)+2,NY=NJCV*2**(NGIT-1)+2,
     *           NXYA=NICV*NJCV*(4**NGIT-1)/3+NICV*2*(2**NGIT-1)+
     *           NJCV*2*(2**NGIT-1)+4*NGIT,NXY=NX*NY,NPHI=5)
      COMMON /GEO/ X(NXYA),Y(NXYA),XC(NXYA),YC(NXYA),FX(NXYA),FY(NXYA),
     *      VOL(NXYA),R(NXYA),RC(NXYA)
      COMMON/VARG/ U1(NXYA),V1(NXYA),T1(NXYA),F1G(NXYA),F2G(NXYA),
     *      RESU(NXYA),RESV(NXYA),REST(NXYA)
      COMMON /WAAL/ TAU,SUU,SUP,SVU,SVP,QWAL
      COMMON /SYMB/ FIB,FIP,FIL,FIR,FAC,IJN
C
C.....Define unit vectors  p  and  n
C
      DXB=X(IJC)-XC(IJB)
      DYB=Y(IJC)-YC(IJB)
      RBC=1./SQRT(DXB**2+DYB**2)
      ALF=DXB*RBC
      BET=DYB*RBC
C
C.....Calculate velocities at symmetry
C
      FIP=U1(IJP)*ALF+V1(IJP)*BET
      FIL=U1(IJL)*ALF+V1(IJL)*BET
      FIR=U1(IJR)*ALF+V1(IJR)*BET
      CALL SYMBC(IJP,IJB,IJL,IJR,IJC)
      U1(IJB)=FIB*ALF
      V1(IJB)=FIB*BET
      RETURN
      END
C
C**********************************************************
      SUBROUTINE USER
C**********************************************************
      PARAMETER (ngit=3,nicv=21,njcv=8,NXA=NICV*(2**NGIT-1)+
     *           2*NGIT,NYA=NJCV*(2**NGIT-1)+2*NGIT,
     *           NX=NICV*2**(NGIT-1)+2,NY=NJCV*2**(NGIT-1)+2,
     *           NXYA=NICV*NJCV*(4**NGIT-1)/3+NICV*2*(2**NGIT-1)+
     *           NJCV*2*(2**NGIT-1)+4*NGIT,NXY=NX*NY,NPHI=5)
      COMMON/INDEX/  NI,NJ,NIM,NJM,NIJ,KGIT,IST,JST,IJST,NITP,
     *      LI(NX),LIG(NX),
     *      LSG(NGIT),LSR(NGIT),LSI(NGIT),MIT(NGIT),
     *      NIGIT(NGIT),NJGIT(NGIT),IGIT(NGIT),JGIT(NGIT),IJGIT(NGIT)
      COMMON/LOGIC/  WRIT3,READ3,TEST,AKSI,LOUTS,LOUTE,ISOL,LSOR,
     *      ICAL(NPHI)
      COMMON/CONTRO/ IMON,JMON,IJMON,IU,IV,IP,IEN,ISW,ISF,IPR,JPR,kgs,
     *      ITBS(NXA),ITBN(NXA),JTBW(NYA),JTBE(NYA),beta,gravx,gravy,
     *      NSW(NPHI),SOR(NPHI),RESOR(NPHI),PRTR(NPHI),URF(NPHI),
     *      GDS(NPHI),SNOR(NPHI),SORMAX,SLARGE,ALFA,GREAT,SMALL,TREF
      COMMON/VARF/ DENSIT,VISCOS,FLOREF,XMOMR,PRMR,PRR(NXYA),
     *      U(NXYA),V(NXYA),P(NXYA),PP(NXYA),T(NXYA),DEN(NXYA),
     *      VIS(NXYA),F1(NXYA),F2(NXYA),DP1(NXYA),DP2(NXYA)
      COMMON/VARG/ U1(NXYA),V1(NXYA),T1(NXYA),F1G(NXYA),F2G(NXYA),
     *      RESU(NXYA),RESV(NXYA),REST(NXYA)
      COMMON /GEO/ X(NXYA),Y(NXYA),XC(NXYA),YC(NXYA),FX(NXYA),FY(NXYA),
     *      VOL(NXYA),R(NXYA),RC(NXYA)
      COMMON/COEF/ AE(NXYA),AW(NXYA),AN(NXYA),AS(NXYA),AP(NXYA),
     *      SU(NXYA),SV(NXYA),SP(NXYA),D1(NXYA),D2(NXYA),E1(NXYA),
     *      E2(NXYA)
      COMMON/BCOEF/ BW(NXYA),BE(NXYA),BS(NXYA),BN(NXYA),BP(NXYA),
     *      RES(NXYA)
      COMMON /WAAL/ TAU,SUU,SUP,SVU,SVP,QWAL
      COMMON /SYMB/ FIB,FIP,FIL,FIR,FAC,IJN
      LOGICAL WRIT3,READ3,TEST,AKSI,LOUTS,LOUTE,ISOL,LSOR,ICAL,
     *        lcont
      DIMENSION PW(NY),TAUS(NX),TAUN(NX),TAUW(NY),TAUE(NY),
     *          QENW(NY),QENE(NY),QENS(NX),QENN(NX)
      CHARACTER TITLE*80
C
C*********************
      ENTRY MODINP
C*********************
      READ(5,6) TITLE
    6 FORMAT(A80)
      READ(5,*) READ3,WRIT3,TEST,LOUTS,LOUTE,lcont
      READ(5,*) IMON,JMON,IPR,JPR,SORMAX,SLARGE,ALFA,NITP
      READ(5,*) DENSIT,VISCOS,PRM,GRAVX,gravy,BETA,
     *          TH,TC,TREF
      READ(5,*) UIN,VIN,TIN,kgs,ulid
      READ(5,*) (ICAL(I),I=1,NPHI)
      READ(5,*) (URF(I),I=1,NPHI)
      READ(5,*) (SOR(I),I=1,NPHI)
      READ(5,*) (NSW(I),I=1,NPHI)
      READ(5,*) (GDS(I),I=1,NPHI)
      READ(5,*) (LSG(I),I=1,NGIT)
      READ(5,*) (LSR(I),I=1,NGIT)
      READ(5,*) (LSI(I),I=1,NGIT)
      READ(5,*) (MIT(I),I=1,NGIT)
C
      IU=1
      IV=2
      IP=3
      IEN=4
      ISW=5
      ISF=6
      SMALL=1.E-20
      GREAT=1.E20
      PRMR=1./PRM
      DO 10 IJ=1,NXYA
      RC(IJ)=1.
      R(IJ)=1.
   10 DEN(IJ)=DENSIT
C
C.....READ GRID DATA GEOMETRY
C
      READ(4) AKSI,NIGIT,NJGIT,IJGIT,IGIT,JGIT,
     *        ITBS,ITBN,JTBW,JTBE,X,Y,XC,YC,VOL,FX,FY
      REWIND 4
      IF(AKSI) THEN
      DO 11 IJ=1,NXYA
      RC(IJ)=YC(IJ)
   11 R(IJ)=Y(IJ)
      ENDIF
C
C.....INITIALISE INTERIOR VARIABLE VALUES
C
      DO 13 LK=1,NGIT
      NI=NIGIT(LK)
      NJ=NJGIT(LK)
      IJST=IJGIT(LK)
      NIJ=NI*NJ
      DO 12 I=2,NI-1
      IJS=IJST+(I-1)*NJ+2
      IJE=IJST+(I-1)*NJ+NJ-1
      DO 12 IJ=IJS,IJE
      U(IJ)=UIN
      V(IJ)=VIN
      T(IJ)=TIN
      PRR(IJ)=PRMR
   12 VIS(IJ)=VISCOS
   13 CONTINUE
      DO 14 I=1,NIGIT(1)
   14 LI(I)=(I-1)*NJGIT(1)
      if (lcont) then
      do 15 l=1,kgs
      ijst=ijgit(l)
      ni=nigit(l)
      nj=njgit(l)
      nij=ni*nj 
      READ(3) ITER,DUMMY,NI,NJ,NIM,NJM,NIJ,
     *   (X(IJ),IJ=IJST+1,IJST+NIJ),(Y(IJ),IJ=IJST+1,IJST+NIJ),
     *   (XC(IJ),IJ=IJST+1,IJST+NIJ),(YC(IJ),IJ=IJST+1,IJST+NIJ),
     *   (F1(IJ),IJ=IJST+1,IJST+NIJ),(F2(IJ),IJ=IJST+1,IJST+NIJ),
     *   (U(IJ),IJ=IJST+1,IJST+NIJ),(V(IJ),IJ=IJST+1,IJST+NIJ),
     *   (P(IJ),IJ=IJST+1,IJST+NIJ),(T(IJ),IJ=IJST+1,IJST+NIJ)
   15 continue
      endif
      RETURN
C
C*************************
      ENTRY BCIN
C*************************
      XMOMR=0.
      FLOREF=0.
C.....READ INLET VALUES OF DEPENDENT VARIABLES
      DO 55 J=2,NJM
      IF(JTBW(J+JST).EQ.1) THEN
      IJ=LI(1)+J
C     READ(5,*) U(IJ),V(IJ),T(IJ)
c-------------------
      v(ij)=0.
      u(ij)=1.
c-------------------
      VIS(IJ)=VISCOS
      DYB=Y(IJ)-Y(IJ-1)
      DXB=X(IJ)-X(IJ-1)
      F1(IJ)=0.5*DEN(IJ)*(U(IJ)*DYB-V(IJ)*DXB)*(R(IJ)+R(IJ-1))
      XMOMR=XMOMR+F1(IJ)*SQRT(U(IJ)**2+V(IJ)**2)
      FLOREF=FLOREF+F1(IJ)
      ENDIF
   55 CONTINUE
      DO 56 J=2,NJM
      IF(JTBE(J+JST).EQ.1) THEN
      IJ=LI(NI)+J
      READ(5,*) U(IJ),V(IJ),T(IJ)
      VIS(IJ)=VISCOS
      DYB=Y(IJ)-Y(IJ-1)
      DXB=X(IJ)-X(IJ-1)
      F1(IJ-NJ)=0.5*DEN(IJ)*(U(IJ)*DYB-V(IJ)*DXB)*(R(IJ)+R(IJ-1))
      XMOMR=XMOMR-F1(IJ-NJ)*SQRT(U(IJ)**2+V(IJ)**2)
      FLOREF=FLOREF-F1(IJ-NJ)
      ENDIF
   56 CONTINUE
      DO 58 I=2,NIM
      IF(ITBS(I+IST).EQ.1) THEN
      IJ=LI(I)+1
      READ(5,*) U(IJ),V(IJ),T(IJ)
      VIS(IJ)=VISCOS
      DXB=X(IJ)-X(IJ-NJ)
      DYB=Y(IJ)-Y(IJ-NJ)
      F2(IJ)=0.5*DEN(IJ)*(V(IJ)*DXB-U(IJ)*DYB)*(R(IJ)+R(IJ-NJ))
      XMOMR=XMOMR+F2(IJ)*SQRT(U(IJ)**2+V(IJ)**2)
      FLOREF=FLOREF+F2(IJ)
      ENDIF
   58 CONTINUE
      DO 59 I=2,NIM
      IJ=LI(I)+NJ
      IF(ITBN(I+IST).EQ.1) THEN
      IJ=LI(I)+NJ
      READ(5,*) U(IJ),V(IJ),T(IJ)
      VIS(IJ)=VISCOS
      DXB=X(IJ)-X(IJ-NJ)
      DYB=Y(IJ)-Y(IJ-NJ)
      F2(IJ-1)=0.5*DEN(IJ)*(V(IJ)*DXB-U(IJ)*DYB)*(R(IJ)+R(IJ-NJ))
      XMOMR=XMOMR-F2(IJ-1)*SQRT(U(IJ)**2+V(IJ)**2)
      FLOREF=FLOREF-F2(IJ-1)
      ENDIF
   59 CONTINUE
c      U(LI(1)+1)=U(LI(1)+2)
C
C-----PRESET WALL velocity or TEMPERATURE
      DO 60 i=2,nim
      IJ=LI(I)+NJ
   60 u(ij)=ulid 
      do 70 j=1,nj
      ij=li(1)+j
      t(ij)=tc
      ij=li(ni)+j
   70 t(ij)=th
C.....SET NORMALISATION FACTORS
C
      DO 80 L=1,NPHI
      SNOR(L)=1.
   80 RESOR(L)=0.0
      IF(FLOREF.LT.1.E-20) FLOREF=1.
      IF(XMOMR.LT.1.E-20) XMOMR=1.
      SNOR(IU)=1./(XMOMR+SMALL)
      SNOR(IV)=SNOR(IU)
      SNOR(IP)=1./(FLOREF+SMALL)
      SNOR(IEN)=1./(FLOREF*TIN+SMALL)
      RETURN
C
C**********************
      ENTRY OUTBC
C**********************
      FLOW=0.
      DO 201 I=2,NIM
      IF(ITBS(I+IST).EQ.2) THEN
      IJ=LI(I)+1
      DXB=X(IJ)-X(IJ-NJ)
      DYB=Y(IJ)-Y(IJ-NJ)
      F2(IJ)=0.5*DEN(IJ)*(V(IJ+1)*DXB-U(IJ+1)*DYB)*(R(IJ)+R(IJ-NJ))
      FLOW=FLOW-F2(IJ)
      ENDIF
  201 CONTINUE
      DO 202 I=2,NIM
      IF(ITBN(I+IST).EQ.2) THEN
      IJ=LI(I)+NJM
      DXB=X(IJ)-X(IJ-NJ)
      DYB=Y(IJ)-Y(IJ-NJ)
      FLOW=FLOW+0.5*DEN(IJ+1)*(V(IJ)*DXB-U(IJ)*DYB)*(R(IJ)+R(IJ-NJ))
      ENDIF
  202 CONTINUE
      DO 203 J=2,NJM
      IF(JTBW(J+JST).EQ.2) THEN
      IJ=LI(1)+J
      DYB=Y(IJ)-Y(IJ-1)
      DXB=X(IJ)-X(IJ-1)
      F1(IJ)=0.5*DEN(IJ)*(U(IJ+NJ)*DYB-V(IJ+NJ)*DXB)*(R(IJ)+R(IJ-1))
      FLOW=FLOW-F1(IJ)
      ENDIF
  203 CONTINUE
      DO 204 J=2,NJM
      IF(JTBE(J+JST).EQ.2) THEN
      IJ=LI(NIM)+J
      DYB=Y(IJ)-Y(IJ-1)
      DXB=X(IJ)-X(IJ-1)
      F1(IJ)=0.5*DEN(IJ+NJ)*(U(IJ)*DYB-V(IJ)*DXB)*(R(IJ)+R(IJ-1))
      FLOW=FLOW+F1(IJ)
      ENDIF
  204 CONTINUE
C
      FAC=FLOREF/(FLOW+SMALL)
C
      DO 205 I=2,NIM
      IF(ITBS(I+IST).EQ.2) THEN
      IJ=LI(I)+1
      U(IJ)=U(IJ+1)*FAC
      V(IJ)=V(IJ+1)*FAC
      F2(IJ)=F2(IJ)*FAC
      T(IJ)=T(IJ+1)
      VIS(IJ)=VIS(IJ+1)
      ENDIF
  205 CONTINUE
      DO 206 I=2,NIM
      IF(ITBN(I+IST).EQ.2) THEN
      IJ=LI(I)+NJ
      U(IJ)=U(IJ-1)*FAC
      V(IJ)=V(IJ-1)*FAC
      T(IJ)=T(IJ-1)
      VIS(IJ)=VIS(IJ-1)
      ENDIF
  206 CONTINUE
      DO 207 J=2,NJM
      IF(JTBW(J+JST).EQ.2) THEN
      IJ=LI(1)+J
      U(IJ)=U(IJ+NJ)*FAC
      V(IJ)=V(IJ+NJ)*FAC
      F1(IJ)=F1(IJ)*FAC
      T(IJ)=T(IJ+NJ)
      VIS(IJ)=VIS(IJ+NJ)
      ENDIF
  207 CONTINUE
      DO 208 J=2,NJM
      IF(JTBE(J+JST).EQ.2) THEN
      IJ=LI(NIM)+J
      F1(IJ)=F1(IJ)*FAC
      U(IJ+NJ)=U(IJ)*FAC
      V(IJ+NJ)=V(IJ)*FAC
      T(IJ+NJ)=T(IJ)
      VIS(IJ+NJ)=VIS(IJ)
      ENDIF
  208 CONTINUE
C-----CORNER VALUES OUTFLOW/SYMMETRY
c      U(LI(NI)+1)=U(LI(NI)+2)
      RETURN
C
C********************
      ENTRY MODPR(KK)
C********************
      RETURN
C
C***********************
      ENTRY CALCDP(KK)
C***********************
      DO 260 I=2,NIM
      PN=0.
      DO 260 J=2,NJM
      IJ=LI(I)+J
      PS=PN
      PE=PP(IJ+NJ)*FX(IJ)+PP(IJ)*(1.-FX(IJ))
      PN=PP(IJ+ 1)*FY(IJ)+PP(IJ)*(1.-FY(IJ))
      DP1(IJ)=PE-PW(J)
      DP2(IJ)=PN-PS
      PW(J)=PE
  260 CONTINUE
      DO 261 I=2,NIM
      IJ=LI(I)+2
      DP2(IJ)=2.*(PP(IJ+ 1)-PP(IJ))*FY(IJ)
      IJ=LI(I)+NJM
  261 DP2(IJ)=2.*(PP(IJ)-PP(IJ-1))*(1.-FY(IJ-1))
      DO 262 J=2,NJM
      IJ=LI(2)+J
      DP1(IJ)=2.*(PP(IJ+NJ)-PP(IJ))*FX(IJ)
      IJ=LI(NIM)+J
  262 DP1(IJ)=2.*(PP(IJ)-PP(IJ-NJ))*(1.-FX(IJ-NJ))
      RETURN
C
C***********************
      ENTRY MODUV (KK)
C***********************
      IST=IGIT(KK)
      JST=JGIT(KK)
C-----SOUTH BOUNDARY
      DO 300 I=2,NIM
      IJ=LI(I)+2
      SUU=0.
      SUP=0.
      SVU=0.
      SVP=0.
      IF(ITBS(I+IST).EQ.1) THEN
      SUU=AS(IJ)*U(IJ-1)
      SUP=AS(IJ)
      SVU=AS(IJ)*V(IJ-1)
      SVP=AS(IJ)
      ELSEIF(ITBS(I+IST).EQ.3) THEN
      CALL SYMVEL(IJ,IJ-1,IJ-NJ,IJ+NJ,IJ-1)
      ELSEIF(ITBS(I+IST).EQ.4 .OR. ITBS(I+IST).EQ.5) THEN
      CALL WALLBC(IJ,IJ-1,IJ-NJ,IJ+NJ,IJ-1,IU)
      ENDIF
      SU(IJ)=SU(IJ)+SUU
      AP(IJ)=AP(IJ)+SUP
      SV(IJ)=SV(IJ)+SVU
      SP(IJ)=SP(IJ)+SVP
  300 AS(IJ)=0.0
C-----NORTH BOUNDARY
      DO 310 I=2,NIM
      IJ=LI(I)+NJM
      SUU=0.
      SUP=0.
      SVU=0.
      SVP=0.
      IF(ITBN(I+IST).EQ.1) THEN
      SUU=AN(IJ)*U(IJ+1)
      SUP=AN(IJ)
      SVU=AN(IJ)*V(IJ+1)
      SVP=AN(IJ)
      ELSEIF(ITBN(I+IST).EQ.3) THEN
      CALL SYMVEL(IJ,IJ+1,IJ+NJ,IJ-NJ,IJ-NJ)
      ELSEIF(ITBN(I+IST).EQ.4 .OR. ITBN(I+IST).EQ.5) THEN
      CALL WALLBC(IJ,IJ+1,IJ+NJ,IJ-NJ,IJ-NJ,IU)
      ENDIF
      SU(IJ)=SU(IJ)+SUU
      AP(IJ)=AP(IJ)+SUP
      SV(IJ)=SV(IJ)+SVU
      SP(IJ)=SP(IJ)+SVP
  310 AN(IJ)=0.0
C-----WEST BOUNDARY
      DO 320 J=2,NJM
      SUU=0.
      SUP=0.
      SVU=0.
      SVP=0.
      IJ=LI(2)+J
      IF(JTBW(J+JST).EQ.1) THEN
      SUU=AW(IJ)*U(IJ-NJ)
      SUP=AW(IJ)
      SVU=AW(IJ)*V(IJ-NJ)
      SVP=AW(IJ)
      ELSEIF(JTBW(J+JST).EQ.3) THEN
      CALL SYMVEL(IJ,IJ-NJ,IJ+1,IJ-1,IJ-NJ-1)
      ELSEIF(JTBW(J+JST).EQ.4 .OR. JTBW(J+JST).EQ.5) THEN
      CALL WALLBC(IJ,IJ-NJ,IJ+1,IJ-1,IJ-NJ-1,IU)
      ENDIF
      SU(IJ)=SU(IJ)+SUU
      AP(IJ)=AP(IJ)+SUP
      SV(IJ)=SV(IJ)+SVU
      SP(IJ)=SP(IJ)+SVP
  320 AW(IJ)=0.0
C-----EAST BOUNDARY
      DO 330 J=2,NJM
      IJ=LI(NIM)+J
      SUU=0.
      SUP=0.
      SVU=0.
      SVP=0.
      IF(JTBE(J+JST).EQ.1) THEN
      SUU=AE(IJ)*U(IJ+NJ)
      SUP=AE(IJ)
      SVU=AE(IJ)*V(IJ+NJ)
      SVP=AE(IJ)
      ELSEIF(JTBE(J+JST).EQ.3) THEN
      CALL SYMVEL(IJ,IJ+NJ,IJ-1,IJ+1,IJ)
      ELSEIF(JTBE(J+JST).EQ.4 .OR. JTBW(J+JST).EQ.5) THEN
      CALL WALLBC(IJ,IJ+NJ,IJ-1,IJ+1,IJ,IU)
      ENDIF
      SU(IJ)=SU(IJ)+SUU
      AP(IJ)=AP(IJ)+SUP
      SV(IJ)=SV(IJ)+SVU
      SP(IJ)=SP(IJ)+SVP
  330 AE(IJ)=0.0
      RETURN
C***********************
      ENTRY MODT (KK)
C***********************
      IST=IGIT(KK)
      JST=JGIT(KK)
C-----SOUTH BOUNDARY
      DO 400 I=2,NIM
      IJ=LI(I)+2
      SUU=0.
      SUP=0.
      IF(ITBS(I+IST).EQ.1) THEN
      SUU=AS(IJ)*T(IJ-1)
      SUP=AS(IJ)
      ELSEIF (ITBS(I+IST).EQ.3 .OR. ITBS(I+IST).EQ.5) THEN
      FIP=T(IJ)
      FIL=T(IJ-NJ)
      FIR=T(IJ+NJ)
      CALL SYMBC(IJ,IJ-1,IJ-NJ,IJ+NJ,IJ-1)
      T(IJ-1)=FIB
      ELSEIF (ITBS(I+IST).EQ.4) THEN
      CALL WALLBC(IJ,IJ-1,IJ-NJ,IJ+NJ,IJ-1,IEN)
      ENDIF
      SU(IJ)=SU(IJ)+SUU
      SP(IJ)=SP(IJ)+SUP
  400 AS(IJ)=0
C     NORTH BOUNDARY
      DO 410 I=2,NIM
      IJ=LI(I)+NJM
      SUU=0.
      SUP=0.
      IF(ITBN(I+IST).EQ.1) THEN
      SUU=AN(IJ)*T(IJ+1)
      SUP=AN(IJ)
      ELSEIF (ITBN(I+IST).EQ.3 .OR. ITBN(I+IST).EQ.5) THEN
      FIP=T(IJ)
      FIL=T(IJ+NJ)
      FIR=T(IJ-NJ)
      CALL SYMBC(IJ,IJ+1,IJ+NJ,IJ-NJ,IJ-NJ)
      T(IJ+1)=FIB
      ELSEIF (ITBN(I+IST).EQ.4) THEN
      CALL WALLBC(IJ,IJ+1,IJ+NJ,IJ-NJ,IJ-NJ,IEN)
      ENDIF
      SU(IJ)=SU(IJ)+SUU
      SP(IJ)=SP(IJ)+SUP
  410 AN(IJ)=0
C     WEST BOUNDARY
      DO 420 J=2,NJM
      IJ=LI(2)+J
      SUU=0.
      SUP=0.
      IF(JTBW(J+JST).EQ.1) THEN
      SUU=AW(IJ)*T(IJ-NJ)
      SUP=AW(IJ)
      ELSEIF (JTBW(J+JST).EQ.3 .OR. JTBW(J+JST).EQ.5) THEN
      FIP=T(IJ)
      FIL=T(IJ+1)
      FIR=T(IJ-1)
      CALL SYMBC(IJ,IJ-NJ,IJ+1,IJ-1,IJ-NJ-1)
      T(IJ-NJ)=FIB
      ELSEIF (JTBW(J+JST).EQ.4) THEN
      CALL WALLBC(IJ,IJ-NJ,IJ+1,IJ-1,IJ-NJ-1,IEN)
      ENDIF
      SU(IJ)=SU(IJ)+SUU
      SP(IJ)=SP(IJ)+SUP
  420 AW(IJ)=0
C     EAST BOUNDARY
      DO 430 J=2,NJM
      IJ=LI(NIM)+J
      SUU=0.
      SUP=0.
      IF(JTBE(J+JST).EQ.1) THEN
      SUU=AE(IJ)*T(IJ+NJ)
      SUP=AE(IJ)
      ELSEIF (JTBE(J+JST).EQ.3 .OR. JTBE(J+JST).EQ.5) THEN
      FIP=T(IJ)
      FIL=T(IJ-1)
      FIR=T(IJ+1)
      CALL SYMBC(IJ,IJ+NJ,IJ-1,IJ+1,IJ)
      T(IJ+nj)=FIB
      ELSEIF (JTBE(J+JST).EQ.4) THEN
      CALL WALLBC(IJ,IJ+NJ,IJ-1,IJ+1,IJ,IEN)
      ENDIF
      SU(IJ)=SU(IJ)+SUU
      SP(IJ)=SP(IJ)+SUP
  430 AE(IJ)=0
      RETURN
C
C***********************
      ENTRY PRESB(KK)
C***********************
      DO 700 I=2,NIM
      IJ=LI(I)+2
      P(IJ-1)=P(IJ)+(P(IJ)-P(IJ+1))*FY(IJ)
      IJ=LI(I)+NJM
  700 P(IJ+1)=P(IJ)+(P(IJ)-P(IJ-1))*(1.-FY(IJ-1))
      DO 710 J=1,NJ
      IJ=LI(2)+J
      P(IJ-NJ)=P(IJ)+(P(IJ)-P(IJ+NJ))*FX(IJ)
      IJ=LI(NIM)+J
  710 P(IJ+NJ)=P(IJ)+(P(IJ)-P(IJ-NJ))*(1.-FX(IJ-NJ))
      RETURN
C
C***********************
      ENTRY OUTPUT
C***********************
      IF(ICAL(IU)) CALL PRINT(U,'U  VEL')
      IF(ICAL(IV)) CALL PRINT(V,'V  VEL')
      IF(ICAL(IP)) CALL PRINT(P,'PRESS.')
      IF(ICAL(IEN)) CALL PRINT(T,'TEMP. ')
      IF(ICAL(ISW)) CALL PRINT(WR,'  WR  ')
      RETURN
C
C**********************
      ENTRY MOOUT1
C**********************
      WRITE(6,801) TITLE,DENSIT,VISCOS,ALFA
  801 FORMAT(1H1,//,50X,A80,/,50X,40('*'),/,50X,
     *        'FLUID DENSITY    :  ',E10.4,/,50X,
     *        'DYNAMIC VISCOSITY:  ',E10.4,/,50X,
     *        'ALFA  PARAMETER  :  ',F 6.2)
      IF(ICAL(IU)) WRITE(6,806) IU,URF(IU) ,IU,GDS(IU)
      IF(ICAL(IV)) WRITE(6,806) IV,URF(IV) ,IV,GDS(IV)
      IF(ICAL(IP)) WRITE(6,807) IP,URF(IP)
      IF(ICAL(IEN)) WRITE(6,806) IEN,URF(IEN) ,IEN,GDS(IEN)
      IF(ICAL(ISW)) WRITE(6,806) ISW,URF(ISW) ,ISW,GDS(ISW)
  806 FORMAT(50X,'URF(',I2,' ) = ',F4.2,' GDS(',I2,' ) = ',F4.2)
  807 FORMAT(50X,'URF(',I2,' ) = ',F4.2)
      WRITE(6,808)
  808 FORMAT(//)
  810 FORMAT(A40)
      RETURN
C
C************************
      ENTRY MOOUT2
C************************
      IST=IGIT(KGIT)
      JST=JGIT(KGIT)
      qws=0.
      DO 820 I=2,NIM
      IF(ITBS(I+IST).gt.3) THEN
      IJ=LI(I)+2
      CALL WALLBC(IJ,IJ-1,IJ-NJ,IJ+NJ,IJ-1,IU)
      CALL WALLBC(IJ,IJ-1,IJ-NJ,IJ+NJ,IJ-1,IEN)
      TAUS(I)=TAU
      qens(i)=qwal
      qws=qws+qwal
      ENDIF
  820 CONTINUE
      qwn=0.
      DO 830 I=2,NIM
      IF(ITBN(I+IST).gt.3) THEN
      IJ=LI(I)+NJM
      CALL WALLBC(IJ,IJ+1,IJ+NJ,IJ-NJ,IJ-NJ,IU)
      CALL WALLBC(IJ,IJ+1,IJ+NJ,IJ-NJ,IJ-NJ,IEN)
      TAUN(I)=TAU
      qenn(i)=qwal
      qwn=qwn+qwal
      ENDIF
  830 CONTINUE
      qww=0.
      DO 840 J=2,NJM
      IF(JTBW(J+JST).gt.3) THEN
      IJ=LI(2)+J
      CALL WALLBC(IJ,IJ-NJ,IJ+1,IJ-1,IJ-NJ-1,Iu)
      CALL WALLBC(IJ,IJ-NJ,IJ+1,IJ-1,IJ-NJ-1,IEN)
      TAUW(J)=TAU
      qenw(j)=qwal
      qww=qww+qwal
      ENDIF
  840 CONTINUE
      qwe=0.
      DO 850 J=2,NJM
      IF(JTBE(J+JST).gt.3) THEN
      IJ=LI(NIM)+J
      CALL WALLBC(IJ,IJ+NJ,IJ-1,IJ+1,IJ,IU)
      CALL WALLBC(IJ,IJ+NJ,IJ-1,IJ+1,IJ,IEN)
      TAUE(J)=TAU
      qene(j)=qwal
      qwe=qwe+qwal
      ENDIF
  850 CONTINUE
      WRITE(6,855) ' TAUS',(TAUS(I),I=1,NI)
      WRITE(6,855) ' TAUN',(TAUN(I),I=1,NI)
      WRITE(6,856) ' TAUW',(TAUW(J),J=1,NJ)
      WRITE(6,856) ' TAUE',(TAUE(J),J=1,NJ)
      IF(ICAL(IEN)) THEN
      WRITE(6,855) ' QS  ',(QENS(I),I=1,NI)
      WRITE(6,855) ' QN  ',(QENN(I),I=1,NI)
      WRITE(6,856) ' QW  ',(QENW(J),J=1,NJ)
      WRITE(6,856) ' QE  ',(QENE(J),J=1,NJ)
      write(6,*) '  qwall_south = ',qws
      write(6,*) '  qwall_north = ',qwn
      write(6,*) '  qwall_west = ',qww
      write(6,*) '  qwall_east = ',qwe
      ENDIF
  855 FORMAT(/,50X,A5,'(I), I=1,NI',/,(5X,1P10E12.4))
  856 FORMAT(/,50X,A5,'(J), J=1,NJ',/,(5X,1P10E12.4))
C
      CALL STREAM
      PSIMIN=1.E10
      psimax=-1.e10
      DO 860 I=2,NIM-1
      DO 860 J=2,NJM-1
      IJ=LI(I)+J
      IF(SP(IJ).LT.PSIMIN) THEN
      PSIMIN=SP(IJ)
      XPSIM=X(IJ)
      YPSIM=Y(IJ)
      ENDIF
      if(sp(ij).gt.psimax) then
      psimax=sp(ij)
      xpsix=x(ij)
      ypsix=y(ij)
      endif
  860 CONTINUE
      WRITE(6,*) ' PSIMIN = ',PSIMIN,' , XPSIM = ',XPSIM,
     *           ' , YPSIM = ',YPSIM
      WRITE(6,*) ' PSIMAX = ',PSIMAX,' , XPSIX = ',XPSIX,
     *           ' , YPSIX = ',YPSIX
      WRITE(6,*) '   '
C.....Velocity Profile along vertical centerline
      IM1=NI/2
      DO 910 J=2,NJ
      IJ=LI(IM1)+J
      qene(j)=t(ij+nj)*fx(ij)+t(ij)*(1.-fx(ij))
      TAUE(J)=u(IJ+NJ)*FX(IJ)+u(IJ)*(1.-FX(IJ))
  910 TAUW(J)=0.5*(Y(IJ)+Y(IJ-1))
C.....Velocity Profile along horizontal centerline
      JM1=NJ/2
      DO 920 I=2,NI
      IJ=LI(I)+JM1
      qenn(i)=t(ij+1)*fy(ij)+t(ij)*(1.-fy(ij))
      TAUN(I)=v(IJ+1)*FY(IJ)+v(IJ)*(1.-FY(IJ))
  920 TAUS(I)=0.5*(X(IJ)+X(IJ-NJ))
      WRITE(6,856) 'UPROF',(TAUE(J),J=1,NJ)
      WRITE(6,856) 'TPROF',(qene(J),J=1,NJ)
      WRITE(6,856) 'YPROF',(TAUW(J),J=1,NJ)
      WRITE(6,856) 'VPROF',(TAUN(I),I=1,NI)
      WRITE(6,856) 'TPROF',(qenn(I),I=1,NI)
      WRITE(6,856) 'XPROF',(TAUS(I),I=1,NI)
      RETURN
C
      END
