INUND.FOR

C  Example 2-3, Inundation Analysis in Ogura Reclaimed Area
C  Revised on 26 Nov., 2001
C
      COMMON /CONST/IMAX,JMAX,DT,DX,DY,G,EPS,
     1        DT2DX,DT2DY,DTGDX,DTGDY,DT2,DXDY
      COMMON /GEODT/IP(40,50),ZB(40,50),RN(40,50),RNGX(40,50),
     1        RNGY(40,50)
      COMMON /SOLTN/SMO(40,50),SMN(40,50),SNO(40,50),SNN(40,50),
     1        HO(40,50),HN(40,50),ZS(40,50),HCV(40,50),SMXCV(40,50),
     2        SNYCV(40,50)
      COMMON /BRKPT/IBR,JBR,IBRD,JBRD,QBR
      COMMON /QHYDR/NHT,TRLX,QHYD(40)
      COMMON /OUTPT/IJOUTM,IOUT(1000),JOUT(1000)
      CHARACTER IP*1
C
C  NPRINT=INTERVAL for PRINT OUTPUT, NFILE=INTERVAL for FILE OUTPUT,
C  NFILE MUST BE INTEGER TIMES of NPRINT.
      DATA  NPRINT, NFILE/ 180, 360/
C
C  INPUT of DATA, and INITIALIZATION
      CALL START
      NFINAL=FLOAT(NHT-1)*3600.0/DT2
C
      NSTEP=0
      TIME=FLOAT(NSTEP)*DT2
C
C  DEFINITION of OUTPUT FILE
      OPEN(20,FILE='SOLTN.DAT',STATUS='UNKNOWN')
C
C  CHECK of WATER VOLUME
      VIN=0.0
      S0=0.0
      GO TO 400
C
  300 NSTEP=NSTEP+1
      TIME=FLOAT(NSTEP)*DT2
      CALL QBREAK( TIME, QBR )
      CALL INDFLW( NSTEP, TIME )
      VIN=VIN+QBR*DT2
C
C  VARIABLES for CONVECTIVE TERM COMPUTATION
      DO 21 I=1,IMAX
      DO 21 J=1,JMAX
      HCV(I,J)=(HO(I,J)+HN(I,J))*0.5
      SMXCV(I,J)=(SMO(I,J)+SMN(I,J))*0.5
      SNYCV(I,J)=(SNO(I,J)+SNN(I,J))*0.5
   21 CONTINUE
C
C  REPLACEMENT of VARIABLES for NEXT STEP
      DO 25 I=1,IMAX
      DO 25 J=1,JMAX
      SMO(I,J)=SMN(I,J)
      SNO(I,J)=SNN(I,J)
      HO(I,J)=HN(I,J)
      ZS(I,J)=ZB(I,J)+HO(I,J)
   25 CONTINUE
C
  400 CONTINUE
      N1=NSTEP/NPRINT
      IF(NSTEP.NE.N1*NPRINT) GO TO 450
C  CHECK of WATER VOLUME
      SN=0.0
      DO 10 I=1,IMAX
      DO 10 J=1,JMAX
      IF(IP(I,J).EQ.'M'.OR.IP(I,J).EQ.'B') GO TO 10
      SN=SN+HO(I,J)
   10 CONTINUE
      SN=SN*DXDY
      S0VIN=S0+VIN
C
      HOUR=TIME/3600.0
      IHR=HOUR
      FMIN=(HOUR-FLOAT(IHR))*60.0
      WRITE(*,260)NSTEP,TIME,HOUR,IHR,FMIN,QBR,S0,SN,VIN
  260 FORMAT(1H /1H ,'NSTEP=',I6,'  T(SEC)=',F10.1,'  HOUR=',F6.2,
     1  '     HR:MIN=',I3,':',F5.2/1H ,3X,'QBR(M3/S)=',F6.1,
     2  '  S0, SN, VIN=',1P,3E14.6,0P)
C
      N1=NSTEP/NFILE
      IF(NSTEP.NE.N1*NFILE) GO TO 450
      WRITE(20,200)NSTEP,TIME,HOUR,IHR,FMIN,NFINAL,IJOUTM,
     1      QBR,S0,SN,VIN
  200 FORMAT('NSTEP=',I6,'  TIME(SEC)=',F10.1,F6.2,I4,F6.2,2I6/
     1  10X,1P,4E15.6,0P)
      WRITE(20,201)(IOUT(I),JOUT(I),HO(IOUT(I),JOUT(I)),
     1              SMO(IOUT(I),JOUT(I)),SNO(IOUT(I),JOUT(I)),
     2              ZS(IOUT(I),JOUT(I)),I=1,IJOUTM)
  201 FORMAT(2I5,1P,4E15.6,0P)
C
  450 IF(NSTEP.LT.NFINAL) GO TO 300
      STOP 'SUCCESSFULLY ENDED'
      END
C
      BLOCK DATA
      COMMON /CONST/IMAX,JMAX,DT,DX,DY,G,EPS,
     1        DT2DX,DT2DY,DTGDX,DTGDY,DT2,DXDY
      COMMON /BRKPT/IBR,JBR,IBRD,JBRD,QBR
      DATA   IMAX, JMAX,  DX,      DY,     G,    EPS,    DT/
     1       36,   48,    285.44,  231.0,  9.8,  0.001,  2.5/
      DATA   IBR, JBR, IBRD, JBRD/ 29, 42, -1, 0/
      END
C
C  INPUT of DATA,  and INITIALIZATION
      SUBROUTINE START
      COMMON /CONST/IMAX,JMAX,DT,DX,DY,G,EPS,
     1        DT2DX,DT2DY,DTGDX,DTGDY,DT2,DXDY
      COMMON /GEODT/IP(40,50),ZB(40,50),RN(40,50),RNGX(40,50),
     1        RNGY(40,50)
      COMMON /SOLTN/SMO(40,50),SMN(40,50),SNO(40,50),SNN(40,50),
     1        HO(40,50),HN(40,50),ZS(40,50),HCV(40,50),SMXCV(40,50),
     2        SNYCV(40,50)
      COMMON /BRKPT/IBR,JBR,IBRD,JBRD,QBR
      COMMON /QHYDR/NHT,TRLX,QHYD(40)
      COMMON /CVTXY/CUM(40,50),CVM(40,50),CUN(40,50),CVN(40,50)
      COMMON /FROF/IFROF(40,50),JFROF(40,50)
      COMMON /OUTPT/IJOUTM,IOUT(1000),JOUT(1000)
      CHARACTER DUMMY*20,IP*1,IFROF*1,JFROF*1,FLDN*70
C
C  INPUT of GEOGRAPHICAL INFORMATION
      OPEN(10,FILE='GEO2D.DAT',STATUS='OLD')
      READ(10,100)IMAX,JMAX,DX,DY
  100 FORMAT(2I5,2F10.2)
C
      READ(10,101)DUMMY
  101 FORMAT(A20)
      DO 11 J=1,JMAX
      READ(10,102)(IP(I,J),I=1,IMAX)
  102 FORMAT(5X,40A1)
   11 CONTINUE
C
      READ(10,103)DUMMY
  103 FORMAT(A20)
      DO 13 J=1,JMAX
      READ(10,104)(ZB(I,J),I=1,IMAX)
  104 FORMAT(5X,8F8.2)
   13 CONTINUE
C
      READ(10,105)DUMMY
  105 FORMAT(A20)
      DO 15 J=1,JMAX
      READ(10,106)(RN(I,J),I=1,IMAX)
  106 FORMAT(5X,8F8.3)
   15 CONTINUE
      CLOSE(10)
C
C  CHANGE of IP for LEVEE-BREAK POINT
      IP(IBR,JBR)='B'
C
C  MESH WHERE HO, SMO, SNO ARE OUTPUTTED.
      II=0
      DO 16 I=1,IMAX
      DO 16 J=1,JMAX
      IF(IP(I,J).EQ.'M') GO TO 16
      II=II+1
      IOUT(II)=I
      JOUT(II)=J
   16 CONTINUE
      IJOUTM=II
C
      WRITE(*,260)IMAX,JMAX,DX,DY
  260 FORMAT(1H //1H ,'OGURA RECLAIMED AREA, IMAX, JMAX=',2I5,3X,
     1  'DX(m), DY(m)=',2F8.2)
C
C  INPUT of DISCHARGE HYDROGRAPH into INUNDATION AREA
      OPEN(11,FILE='FLOOD.DAT',STATUS='OLD')
      READ(11,110)FLDN,NHT,TRLX,(QHYD(I),I=1,NHT)
  110 FORMAT(A70/I5,F10.1/(8F8.0))
      CLOSE(11)
C
      WRITE(*,261)IBR,JBR,IBRD,JBRD
  261 FORMAT(1H //1H ,'LEVEE BREAKAGE, IBR, JBR, IBRD, JBRD=',4I5)
      WRITE(*,262)FLDN,NHT,TRLX,(QHYD(I),I=1,NHT)
  262 FORMAT(1H //1H ,'DISCHARGE HYDROGRAPH Q(m3/s)-T(hour)'//
     1  1H ,2X,'FLOOD= ',A70/1H ,5X,'NHT=',I5,'  T RELAX(sec)=',F10.1/
     2  (1H ,2X,8F8.1))
C
C  INITIALIZATION of VARIABLES
      DO 20 I=1,IMAX
      DO 20 J=1,JMAX
      SMO(I,J)=0.0
      SNO(I,J)=0.0
      HO(I,J)=0.0
      ZS(I,J)=ZB(I,J)
      SMN(I,J)=0.0
      SNN(I,J)=0.0
      HN(I,J)=0.0
      SMXCV(I,J)=0.0
      SNYCV(I,J)=0.0
      HCV(I,J)=0.0
      CUM(I,J)=0.0
      CVM(I,J)=0.0
      CUN(I,J)=0.0
      CVN(I,J)=0.0
      IFROF(I,J)='N'
      JFROF(I,J)='N'
      IF(I.EQ.1) GO TO 21
      IF(IP(I-1,J).NE.'I') GO TO 21
      RNGX(I,J)=((RN(I,J)+RN(I-1,J))/2.0)**2*G*DT
      GO TO 22
   21 RNGX(I,J)=RN(I,J)**2*G*DT
   22 IF(J.EQ.1) GO TO 23
      IF(IP(I,J-1).NE.'I') GO TO 23
      RNGY(I,J)=((RN(I,J)+RN(I,J-1))/2.0)**2*G*DT
      GO TO 20
   23 RNGY(I,J)=RN(I,J)**2*G*DT
   20 CONTINUE
C
C  CONSTANTS
      DT2=DT*2.0
      DT2DX=DT*2.0/DX
      DT2DY=DT*2.0/DY
      DTGDX=DT*G/DX*2.0
      DTGDY=DT*G/DY*2.0
      DXDY=DX*DY
      RETURN
      END
C
C  DISCHARGE into INUNDATION AREA through LEVEE-BREAKAGE
      SUBROUTINE QBREAK( TIME, QBR )
      COMMON /QHYDR/NHT,TRLX,QHYD(40)
      IF(TIME.LT.TRLX) GO TO 300
      THR=TIME/3600.0
      IT=THR+1.0
      IF(IT.GE.NHT) IT=NHT-1
      THR0=FLOAT(IT-1)
      QBR=(QHYD(IT+1)-QHYD(IT))*(THR-THR0)+QHYD(IT)
      RETURN
  300 QBR=QHYD(1)/TRLX*TIME
      RETURN
      END
C
C  2-D INUNDATION FLOW
      SUBROUTINE INDFLW( NSTEP, TIME )
      COMMON /CONST/IMAX,JMAX,DT,DX,DY,G,EPS,
     1        DT2DX,DT2DY,DTGDX,DTGDY,DT2,DXDY
      COMMON /GEODT/IP(40,50),ZB(40,50),RN(40,50),RNGX(40,50),
     1        RNGY(40,50)
      COMMON /SOLTN/SMO(40,50),SMN(40,50),SNO(40,50),SNN(40,50),
     1        HO(40,50),HN(40,50),ZS(40,50),HCV(40,50),SMXCV(40,50),
     2        SNYCV(40,50)
      COMMON /BRKPT/IBR,JBR,IBRD,JBRD,QBR
      COMMON /CVTXY/CUM(40,50),CVM(40,50),CUN(40,50),CVN(40,50)
      COMMON /FROF/IFROF(40,50),JFROF(40,50)
      CHARACTER IP*1,IFROF*1,JFROF*1
C
      CALL FROVF
      CALL CONVX
      CALL CONVY
C
C  DISCHARGE FLUX
      DO 10 I=1,IMAX
      DO 10 J=1,JMAX
      IF(I.EQ.1.OR.J.EQ.1) GO TO 10
      IF(IP(I,J).EQ.'M'.OR.IP(I,J).EQ.'B') GO TO 30
C
C  X-DIRECTION
      IF(IP(I-1,J).EQ.'M') GO TO 15
      IF(IFROF(I,J).EQ.'Y') GO TO 11
      IF((ZS(I,J).GT.ZS(I-1,J).AND.HO(I,J).LE.EPS).OR.
     1   (ZS(I,J).LT.ZS(I-1,J).AND.HO(I-1,J).LE.EPS)) GO TO 15
      HH=(HO(I,J)+HO(I-1,J))*0.5
      IF(HH.LE.EPS) GO TO 12
      UX=SMO(I,J)/HH
      VY=(SNO(I-1,J)+SNO(I,J)+SNO(I,J+1)+SNO(I-1,J+1))*0.25/HH
      SFM=RNGX(I,J)*SQRT(UX**2+VY**2)/HH**1.3333333
      GO TO 13
   12 SFM=0.0
   13 SMN(I,J)=(SMO(I,J)*(1.0-SFM)
     1         -(CUM(I,J)-CUM(I-1,J))*DT2DX-(CVM(I,J+1)-CVM(I,J))*DT2DY
     3         -(ZS(I,J)-ZS(I-1,J))*HH*DTGDX)/(1.0+SFM)
      IF((HO(I,J).LE.EPS.AND.SMN(I,J).LT.0.0).OR.
     1   (HO(I-1,J).LE.EPS.AND.SMN(I,J).GT.0.0)) GO TO 15
   11 IF(ABS(SMN(I,J)).LT.5.0E-5) SMN(I,J)=0.0
      GO TO 20
   15 SMN(I,J)=0.0
C
C  Y-DIRECTION
   20 CONTINUE
      IF(IP(I,J-1).EQ.'M') GO TO 25
      IF(JFROF(I,J).EQ.'Y') GO TO 21
      IF((ZS(I,J).GT.ZS(I,J-1).AND.HO(I,J).LE.EPS).OR.
     1   (ZS(I,J).LT.ZS(I,J-1).AND.HO(I,J-1).LE.EPS)) GO TO 25
      HH=(HO(I,J)+HO(I,J-1))*0.5
      IF(HH.LE.EPS) GO TO 22
      UX=(SMO(I,J)+SMO(I,J-1)+SMO(I+1,J-1)+SMO(I+1,J))*0.25/HH
      VY=SNO(I,J)/HH
      SFN=RNGY(I,J)*SQRT(UX**2+VY**2)/HH**1.3333333
      GO TO 23
   22 SFN=0.0
   23 SNN(I,J)=(SNO(I,J)*(1.0-SFN)
     1         -(CUN(I+1,J)-CUN(I,J))*DT2DX-(CVN(I,J)-CVN(I,J-1))*DT2DY
     3         -(ZS(I,J)-ZS(I,J-1))*HH*DTGDY)/(1.0+SFN)
      IF((HO(I,J).LE.EPS.AND.SNN(I,J).LT.0.0).OR.
     1   (HO(I,J-1).LE.EPS.AND.SNN(I,J).GT.0.0)) GO TO 25
   21 IF(ABS(SNN(I,J)).LT.5.0E-5) SNN(I,J)=0.0
      GO TO 10
   25 SNN(I,J)=0.0
      GO TO 10
   30 SMN(I,J)=0.0
      SNN(I,J)=0.0
   10 CONTINUE
C
C  POINT of LEVEE BREAK
      IF(IBRD.NE.0.AND.JBRD.NE.0) GO TO 31
      IF(IBRD.EQ.0) GO TO 32
      IF(JBRD.EQ.0) GO TO 33
      GO TO 35
   31 QQ=QBR/(DX+DY)
      SMN(IBR,JBR)=FLOAT(IBRD)*QQ
      SNN(IBR,JBR)=FLOAT(JBRD)*QQ
      GO TO 39
   32 QQ=QBR/DX
      SMN(IBR,JBR)=0.0
      SNN(IBR,JBR)=FLOAT(JBRD)*QQ
      GO TO 39
   33 QQ=QBR/DY
      SMN(IBR,JBR)=FLOAT(IBRD)*QQ
      SNN(IBR,JBR)=0.0
      GO TO 39
   35 SMN(IBR,JBR)=0.0
      SNN(IBR,JBR)=0.0
   39 CONTINUE
C
C  WATER-DEPTH
      DO 40 I=1,IMAX
      DO 40 J=1,JMAX
      IF(IP(I,J).EQ.'M'.OR.IP(I,J).EQ.'B') GO TO 41
      HN(I,J)=HO(I,J)-(SMN(I+1,J)-SMN(I,J))*DT2DX
     1               -(SNN(I,J+1)-SNN(I,J))*DT2DY
      IF(HN(I,J).GE.0.0) GO TO 40
   41 HN(I,J)=0.0
   40 CONTINUE
      RETURN
      END
C
C  CONVECTIVE TERM in X-WISE MONMENTUM EQUATION
      SUBROUTINE CONVX
      COMMON /CONST/IMAX,JMAX,DT,DX,DY,G,EPS,
     1        DT2DX,DT2DY,DTGDX,DTGDY,DT2,DXDY
      COMMON /GEODT/IP(40,50),ZB(40,50),RN(40,50),RNGX(40,50),
     1        RNGY(40,50)
      COMMON /SOLTN/SMO(40,50),SMN(40,50),SNO(40,50),SNN(40,50),
     1        HO(40,50),HN(40,50),ZS(40,50),HCV(40,50),SMXCV(40,50),
     2        SNYCV(40,50)
      COMMON /CVTXY/CUM(40,50),CVM(40,50),CUN(40,50),CVN(40,50)
      CHARACTER IP*1
      DATA EPSCV/ 0.100 /
C
      DO 10 I=1,IMAX
      DO 10 J=1,JMAX
      IF(I.EQ.1.OR.J.EQ.1) GO TO 10
C
C  CONVECTIVE TERM:  D(UM)/DX on (I+1/2,J) - (I+1/2,J+1)
      IF(IP(I,J).EQ.'M') GO TO 610
      HHE=(HCV(I+1,J)+HCV(I,J))*0.5
      HH=(HCV(I,J)+HCV(I-1,J))*0.5
      IF(HHE.LE.EPSCV) GO TO 601
      UEE=SMXCV(I+1,J)/HHE
      GO TO 602
  601 UEE=0.0
  602 IF(HH.LE.EPSCV) GO TO 603
      UU=SMXCV(I,J)/HH
      GO TO 604
  603 UU=0.0
  604 UE=(UEE+UU)*0.5
      IF(UE.GE.0.0) GO TO 605
      CUM(I,J)=UE*SMXCV(I+1,J)
      GO TO 615
  605 CUM(I,J)=UE*SMXCV(I,J)
      GO TO 615
  610 CUM(I,J)=0.0
  615 CONTINUE
C
C  CONVECTIVE TERM: D(VM)/DY on (I-1/2,J) - (I+1/2,J)
      HHSW=(HCV(I-1,J)+HCV(I-1,J-1))*0.5
      HHSE=(HCV(I,J-1)+HCV(I,J))*0.5
      IF(HHSW.LE.EPSCV) GO TO 621
      VSW=SNYCV(I-1,J)/HHSW
      GO TO 622
  621 VSW=0.0
  622 IF(HHSE.LE.EPSCV) GO TO 623
      VSE=SNYCV(I,J)/HHSE
      GO TO 624
  623 VSE=0.0
  624 VS=(VSE+VSW)*0.5
      IF(VS.GE.0.0) GO TO 625
      CVM(I,J)=VS*SMXCV(I,J)
      GO TO 10
  625 CVM(I,J)=VS*SMXCV(I,J-1)
   10 CONTINUE
      RETURN
      END
C
C  CONVECTIVE TERM in Y-WISE MOMENTUM EQUATION
      SUBROUTINE CONVY
      COMMON /CONST/IMAX,JMAX,DT,DX,DY,G,EPS,
     1        DT2DX,DT2DY,DTGDX,DTGDY,DT2,DXDY
      COMMON /GEODT/IP(40,50),ZB(40,50),RN(40,50),RNGX(40,50),
     1        RNGY(40,50)
      COMMON /SOLTN/SMO(40,50),SMN(40,50),SNO(40,50),SNN(40,50),
     1        HO(40,50),HN(40,50),ZS(40,50),HCV(40,50),SMXCV(40,50),
     2        SNYCV(40,50)
      COMMON /CVTXY/CUM(40,50),CVM(40,50),CUN(40,50),CVN(40,50)
      CHARACTER IP*1
      DATA EPSCV/ 0.100 /
C
      DO 10 I=1,IMAX
      DO 10 J=1,JMAX
      IF(I.EQ.1.OR.J.EQ.1) GO TO 10
C
C  CONVECTIVE TERM:  D(VN)/DY on (I,J+1/2) - (I+1,J+1/2)
      IF(IP(I,J).EQ.'M') GO TO 610
      HHN=(HCV(I,J+1)+HCV(I,J))*0.5
      HH=(HCV(I,J)+HCV(I,J-1))*0.5
      IF(HHN.LE.EPSCV) GO TO 601
      VNN=SNYCV(I,J+1)/HHN
      GO TO 602
  601 VNN=0.0
  602 IF(HH.LE.EPSCV) GO TO 603
      VV=SNYCV(I,J)/HH
      GO TO 604
  603 VV=0.0
  604 VN=(VNN+VV)*0.5
      IF(VN.GE.0.0) GO TO 605
      CVN(I,J)=VN*SNYCV(I,J+1)
      GO TO 615
  605 CVN(I,J)=VN*SNYCV(I,J)
      GO TO 615
  610 CVN(I,J)=0.0
  615 CONTINUE
C
C  CONVECTIVE TERM: D(UN)/DX on (I,J-1/2) - (I,J+1/2)
      HHNW=(HCV(I-1,J)+HCV(I,J))*0.5
      HHSW=(HCV(I,J-1)+HCV(I-1,J-1))*0.5
      IF(HHNW.LE.EPSCV) GO TO 621
      UNW=SMXCV(I,J)/HHNW
      GO TO 622
  621 UNW=0.0
  622 IF(HHSW.LE.EPSCV) GO TO 623
      USW=SMXCV(I,J-1)/HHSW
      GO TO 624
  623 USW=0.0
  624 UW=(UNW+USW)*0.5
      IF(UW.GE.0.0) GO TO 625
      CUN(I,J)=UW*SNYCV(I,J)
      GO TO 10
  625 CUN(I,J)=UW*SNYCV(I-1,J)
   10 CONTINUE
      RETURN
      END
C
C  FLUX of FREE OVERFALL
      SUBROUTINE FROVF
      COMMON /CONST/IMAX,JMAX,DT,DX,DY,G,EPS,
     1        DT2DX,DT2DY,DTGDX,DTGDY,DT2,DXDY
      COMMON /GEODT/IP(40,50),ZB(40,50),RN(40,50),RNGX(40,50),
     1        RNGY(40,50)
      COMMON /SOLTN/SMO(40,50),SMN(40,50),SNO(40,50),SNN(40,50),
     1        HO(40,50),HN(40,50),ZS(40,50),HCV(40,50),SMXCV(40,50),
     2        SNYCV(40,50)
      COMMON /FROF/IFROF(40,50),JFROF(40,50)
      CHARACTER IP*1,IFROF*1,JFROF*1
      DATA  CQ/ 0.544 /
C  CQ:DISCHARGE COEFFICIENT at FREE OVERFALL
C
      DO 10 I=1,IMAX
      DO 10 J=1,JMAX
      IF(I.EQ.1.OR.J.EQ.1) GO TO 10
      IF(IP(I,J).EQ.'M'.OR.IP(I,J).EQ.'B') GO TO 10
C
C  X-DIRECTION
      IF(IP(I-1,J).EQ.'M') GO TO 20
      IF(ZB(I-1,J).GT.ZS(I,J)) GO TO 11
      IF(ZB(I,J).GT.ZS(I-1,J)) GO TO 12
      IFROF(I,J)='N'
      GO TO 20
   11 ID=1
      HH=HO(I-1,J)
      GO TO 13
   12 ID=-1
      HH=HO(I,J)
   13 IFROF(I,J)='Y'
      IF(HH.LE.EPS) GO TO 14
      SMN(I,J)=FLOAT(ID)*CQ*HH*SQRT(G*HH)
      GO TO 20
   14 SMN(I,J)=0.0
C
C  Y-DIRECTION
   20 IF(IP(I,J-1).EQ.'M') GO TO 10
      IF(ZB(I,J-1).GT.ZS(I,J)) GO TO 21
      IF(ZB(I,J).GT.ZS(I,J-1)) GO TO 22
      JFROF(I,J)='N'
      GO TO 10
   21 ID=1
      HH=HO(I,J-1)
      GO TO 23
   22 ID=-1
      HH=HO(I,J)
   23 JFROF(I,J)='Y'
      IF(HH.LE.EPS) GO TO 24
      SNN(I,J)=FLOAT(ID)*CQ*HH*SQRT(G*HH)
      GO TO 10
   24 SNN(I,J)=0.0
   10 CONTINUE
      RETURN
      END