C     STORM SURGE MODEL FOR SRILANKA REGION
C     WITH CURRENTS  menu driven model
C     IT IS A STAIR STEP MODEL
C     LAST UPDATED 31.5.95
      PROGRAM SIPALA
      DIMENSION VEL(53,72),FLUX1(53,72),FLUX2(53,72),ZETA1(53,72),
     *ZETA2(53,72),F(53,72),MAP(53,72),H(53,72),pz3(53,72)
     +TX(53,72),TY(53,72),XVEL(53,72),YVEL(53,72)
      DIMENSION PZ1(53,72),KO1(53,72),MIN1(53,72)
      DIMENSION PZ2(53,72),KO2(53,72),MIN2(53,72)
      DIMENSION ALAT1(20),ALONG1(20),DIFKOUNT(20),DTOTAL(20),AHOUR(20)
      DIMENSION TRX(53,72),TRY(53,72)
C	,UC(53,72),VC(53,72),CUR(53,72),
c     1DRN(53,72),XUC(53,72),YVC(53,72),WTX(53,72),WTY(53,72),
c    2WCUR(53,72),WDRN(53,72)
      CHARACTER *30FIL1,FIL2,FIL3,FIL4,FIL5,FIL6,FIL7,FIL8,NAME 
      COMMON/A/AL,B
      COMMON/D/ALAT1,ALONG1,DIFKOUNT,DTOTAL
      COMMON/C/VMAX,R,NAME
      
C      
C*************************READ INPUT*****************************
C      
c      CALL INPUT(KS,VMAX,NAME)  
      OPEN(4,FILE='RESDEP')
      OPEN(5,FILE='SRIMENU.DAT')
            
      DO 577 J=1,72
      READ(4,*)(H(I,J),I=1,53)
 577  CONTINUE
      READ(5,101)M,N,K1,INITIA
101   FORMAT(I4,3(/I4))
      READ(5,102)AL,PHI,DT,B
102   FORMAT(F10.0,3(/F10.0))
C
C******************************************************************    
C*********************DIFFERENT SCENARIOS**************************  
C
c      CALL CS(23)
c      WRITE(*,*)'A. SCENARIO(1):TO RUN THE PROGRAM FOR PREDICTED' 
c      WRITE(*,*)'   TRACK'
c      CALL CS(1)
c      WRITE(*,*)'B. SCENARIO(2):TO RUN THE PROGRAM FOR THE TRACK' 
c      WRITE(*,*)'   100 KM LEFT OF PREDICTED TRACK'
c      CALL CS(1)
c      WRITE(*,*)'C. SCENARIO(3):TO RUN THE PROGRAM FOR THE TRACK'
c      WRITE(*,*)'   100 KM RIGHT OF PREDICTED TRACK'
c      CALL CS(1)
c      WRITE(*,*)'SELECT REQUIRED SCENARIO'
c      WRITE(*,*)'ENTER (1) FOR SCENARIO(1)'
c      WRITE(*,*)'AND ENTER (2) FOR ALL ABOVE SCENARIOS:>'
c      READ(*,*)ANS      
c      IF(ANS.GT.1)GO TO 2222 
c2222  IF(ANS.EQ.2)GO TO 5656 
c5656   IIGO=1
c4444  IF(ANS.EQ.2.AND.IIGO.EQ.1)GO TO 3333
c      IF(ANS.EQ.2.AND.IIGO.EQ.2)GO TO 3334
c      IF(ANS.EQ.2.AND.IIGO.EQ.3)GO TO 3335
c3333    WRITE(*,*)'ENTER DATA FOR THE PREDICTED TRACK:>'
c      READ(*,*)IC
c      GO TO 3336        
c3334  WRITE(*,*)'ENTER DATA FOR THE TRACK 100 KM LEFT OF'
c      WRITE(*,*)'PREDICTED TRACK:>'
c      READ(*,*)IC
C      GO TO 3336
c3335  WRITE(*,*)'ENTER DATA FOR THE TRACK 100 KM RIGHT OF'
c      WRITE(*,*)'PREDICTED TRACK:>'   
cC      READ(*,*)IC
c3336  CALL CS(23)
c      WRITE(*,*)'ENTER NAME OF THE FILE FOR RESULT:>'
c      READ(*,1)FIL1
c      WRITE(*,*)'ENTER NAME OF THE FILE FOR ZETA2:>'
c      READ(*,1)FIL2
c      WRITE(*,*)'ENTER NAME OF THE FILE FOR PZ1:>'
c      READ(*,1)FIL3
c      WRITE(*,*)'ENTER NAME OF THE FILE FOR CURRENTS:>'
c     READ(*,1)FIL4
c      WRITE(*,*)'ENTER NAME OF THE FILE FOR CURRENTS IN SEA:>'
c      READ(*,1)FIL5
c      WRITE(*,*)'ENTER NAME OF THE FILE FOR WINDS :>'
c      READ(*,1)FIL6
c      WRITE(*,*)'ENTER NAME OF THE FILE FOR WINDS IN SEA:>'
c      READ(*,1)FIL7
c      WRITE(*,*)'ENTER NAME OF THE FILE FOR PEAK SURGE CONTOURS:>'
c      READ(*,1)FIL8      
c1     FORMAT(A30)      
      OPEN(6,FILE='a.res')
      OPEN(11,FILE='az2.dat')
C      OPEN(12,FILE='apz1.dat')      
C      OPEN(13,FILE='ac.dat')
C      OPEN(15,FILE='acs.dat')
C      OPEN(16,FILE='aw.dat')
C      OPEN(17,FILE='aws.dat')      
      OPEN(18,FILE='pz3.dat')
      OPEN(8,FILE='inp.dat')

c      WRITE(*,*)'ENTER RADIUS OF MAXIMUM WINDS IN METERS:>'
c      READ(*,*)R
c      CALL CS(23)
c      WRITE(*,*)'ENTER NO OF TRACK POSITIONS:>'
c      READ(*,*)ND
c      CALL CS(23)
c      WRITE(*,*)'ENTER LATITUDE, LONGITUDE AND TIME IN HOURS:>'
c      CALL CS(1)
       read(8,*)name
       read(8,*)VMAX
       read(8,*)KS
       read(8,*)R
       read(8,*)ND
      DO 9019 K=1,ND
      READ(8,*)ALAT1(K),ALONG1(K),AHOUR(K)
c1312  FORMAT(3F6.2)      
9019  CONTINUE
      WRITE(6,8676)
8676  FORMAT(54H******************************************************,
     *26H**************************,/,
     *24X,27HIIT DELHI STORM SURGE MODEL,/,
     *22X,31HCENTER FOR ATMOSPHERIC SCIENCES,/,
     *18X,41HINDIAN INSTITUTE OF TECHNOLOGY, NEW DELHI,/,
     *62H**************************************************************,
     *18H******************)      
      WRITE(*,*)  
      WRITE(6,*)M,N,ks,K1,INITIA,ND     
      WRITE(*,*)     
      WRITE(6,205)AL,B,PHI,DT,R,VMAX
205   FORMAT(1X,'NORTH SOUTH EXTENT OF ANALYSIS AREA IN METERS',
     *3X,F10.0/
     *' EAST WEST EXTENT OF ANALYSIS AREA IN METRES',5X,F10.0/
     *' LATITUDE OF OPEN SEA BOUNDARY IN RADIANS',7X,F10.4/
     *' TIME STEP IN SEC ',31X,F10.0/
     +' RADIUS OF MAXIMUM WIND',31X,F10.0/
     +' MAXIMUM WIND IN MPS ',31X,F10.0/)       
      WRITE(6,*)
      WRITE(6,*)'TRACK OF THE CYCLONE'
      WRITE(6,*)'LATITUDE      LONGITUDE'
      DO 3776 I=1,ND
      WRITE(6,*)ALAT1(I),ALONG1(I)
3776  CONTINUE     
C
C    ******************************************************************
C    ******************************************************************
C                        MAIN PROGRAM BEGINS           
C    ******************************************************************     
     
      M1=M-1
      M2=M-2
      N1=N-1
      N2=N-2
      CD=2.8E-3
      CF=2.6E-3
      RHO=1.026E3
      RHO1=1.0/RHO
      RHOAIR=1.226
      DX=B/FLOAT(M1)
      DX2=1.0/(2.0*DX)
      DX4=1.0/(4.0*DX)
      DY=AL/FLOAT(N1)
      DY2=1.0/(2.0*DY)
      DY4=1.0/(4.0*DY)
      DO  50 JJ =1,N
      J=N+1-JJ
      READ(5,103)(MAP(I,J),I=1,M)
C     WRITE(6,222)J,(MAP(I,J),I=1,M)
222   FORMAT(1X,I2,5X,53I2)
50    CONTINUE
103   FORMAT(53I1)
      DO 3 J=1,N
      S=2.0*7.292E-5*SIN(PHI+FLOAT(J-1)*DY/6371220)
      DO 3 I=1,M
      TX(I,J)=0.0
      TY(I,J)=0.0
      VEL(I,J)=0.0
      FLUX1(I,J)=0.0
      FLUX2(I,J)=0.0
      ZETA1(I,J)=0.0
      ZETA2(I,J)=0.0
      F(I,J)=S
      PZ1(I,J)=0.0
      KO1(I,J)=0
      MIN1(I,J)=0
      PZ2(I,J)=0.0
      pz3(i,j)=0.0
      KO2(I,J)=0
      MIN2(I,J)=0
      IF(MAP(I,J).EQ.0)GOTO 3
3     CONTINUE
      DO 777 I=29,38
      DO 777 J=26,42
      IF(MAP(I,J).EQ.0)GOTO 777
      H(I,J)=0.4*H(I,J)
 777  CONTINUE
      DO 888 I=22,27
      DO 888 J=15,29
      IF(MAP(I,J).EQ.0)GOTO 888
      H(I,J)=0.4*H(I,J)
 888  CONTINUE
      DO 77 I=22,27
      DO 77 J=15,25
      IF(MAP(I,J).EQ.0)GOTO 77
      H(I,J)=0.5*H(I,J)
  77  CONTINUE
      DO 88 I=30,34
      DO 88 J=31,35
      IF(MAP(I,J).EQ.0)GOTO 88
      H(I,J)=0.7*H(I,J)
  88  CONTINUE
      DO 66 I=36,38
      DO 66 J=24,25
      IF(MAP(I,J).EQ.0)GOTO 66
      H(I,J)=0.5*H(I,J)
  66  CONTINUE
      DO 43 I=31,33
      DO 43 J=33,42
      H(I,J)=0.7*H(I,J)
  43  CONTINUE
      H(36,26)=6.0*H(36,26)
      H(26,33)=2.0*H(26,33)
      H(27,33)=3.0+H(27,33)
      H(26,32)=3.0+H(26,32)
      H(22,29)=2.0*H(22,29)
      H(23,29)=3.0+H(23,29)
      H(32,31)=2.0*H(32,31)
      H(32,32)=2.0*H(32,32)
      H(27,32)=2.0+H(27,32)
      
      G=9.81
      G1=1.0/G
      KOUNT=0
      SUM=0.0
      DO 905 K=1,ND
      DIFKOUNT(K)=AHOUR(K)*20.0
      DTOTAL(K)=SUM+DIFKOUNT(K)
      SUM=SUM+DIFKOUNT(K)
  905 CONTINUE
      K=1
c        CALL CS(23)
c        WRITE(*,*)'   INDIAN INSTITUTE OF TECHNOLOGY, NEW DELHI'
c        WRITE(*,*)'        STORM SURGE MODEL IS RUNNING'
9999  CONTINUE

C     UPDATE ELEVATION
      DO 4 I=2,M1,2
      DO 4 J=3,N1,2
      IF(MAP(I,J).EQ.0)GOTO4
      X=(FLUX1(I+1,J)-FLUX1(I-1,J))*DX2
      X=X+(FLUX1(I,J+1)-FLUX1(I,J-1))*DY2
      ZETA2(I,J)=ZETA1(I,J)-X*DT
4     CONTINUE
C     RADIATION CONDITION ON SOUTHERN OPEN BOUNDARY
      DO 530 I=2,M1,2
      IF(MAP(I,1).EQ.0)GOTO 530
      ZETA2(I,1)=-ZETA2(I,3)-2.0*SQRT(H(I,2)*G1)*VEL(I,2)
 530  CONTINUE
C     RADIATION CONDITON ON EASTERN OPEN BOUNDARY
      DO 540 J=1,N1,2
      IF(MAP(M,J).EQ.0)GOTO 540
      ZETA2(M,J)=-ZETA2(M2,J)+2.0*SQRT(H(M1,J)*G1)*VEL(M1,J)
 540  CONTINUE
C     RADIATION CONDITION ON WESTERN OPEN BOUNDARY
      DO 550 J=1,N1,2
      IF(MAP(1,J).EQ.0)GOTO 550
      VEL(1,J)=-VEL(3,J)-2.0*SQRT(G/H(2,J))*ZETA2(2,J)
      FLUX2(1,J)=VEL(1,J)*(ZETA2(1,J)+H(1,J))
 550  CONTINUE
      DO  5 I=3,M2,2
      DO  5 J=1,N1,2
      IF(MAP(I,J).EQ.0)GOTO5
      ZETA2(I,J)=0.5*(ZETA2(I+1,J)+ZETA2(I-1,J))
 5    CONTINUE
      DO  6 J=2,N2,2
      DO  6 I=2,M
      IF(MAP(I,J).EQ.0)GOTO6
      ZETA2(I,J)=0.5*(ZETA2(I,J+1)+ZETA2(I,J-1))
 6    CONTINUE
      DO 7 J=1,N1
      IF(MAP(1,J).EQ.0)GOTO 7
      ZETA2(1,J)=2.0*ZETA2(2,J)-ZETA2(3,J)
 7    CONTINUE
      DO 9 J=1,N
      DO 9 I=4,38
      IF(MAP(I,J).EQ.0.AND.MAP(I+1,J).EQ.1)GOTO70
      IF(MAP(I,J).EQ.0.AND.MAP(I-1,J).EQ.1)GOTO71
      GOTO9
  70  ZETA2(I,J)=0.5*(3.0*ZETA2(I+1,J)-ZETA2(I+3,J))
      GOTO9
  71  ZETA2(I,J)=0.5*(3.0*ZETA2(I-1,J)-ZETA2(I-3,J))
   9  CONTINUE
      DO 14 J=10,40
      DO 14 I=4,38
      IF(MAP(I,J).EQ.0.AND.MAP(I,J+1).EQ.1)GOTO22
      IF(MAP(I,J).EQ.0.AND.MAP(I,J-1).EQ.1)GOTO23
      GO TO 14
  22  ZETA2(I,J)=0.5*(3.0*ZETA2(I,J+1)-ZETA2(I,J+3))
      GO TO 14
  23  ZETA2(I,J)=0.5*(3.0*ZETA2(I,J-1)-ZETA2(I,J-3))
  14  CONTINUE
C     ELEVATION UPDATED
      IF(FLOAT(KOUNT).EQ.DTOTAL(K+1))K=K+1
      XC=X0(FLOAT(KOUNT),K)
      YC=Y0(FLOAT(KOUNT),K)
      DO 51 I=3,M2,2
      X=FLOAT(I-1)*DX
      DO 51 J=3,N1,2
      Y=FLOAT(J-1)*DY
c      IF(MAP(I,J).EQ.0)GO TO 51
      RR=SQRT((X-XC)**2+(Y-YC)**2)
      IF(RR.LT.10.0)GO TO 54
      VTHETA=2.0*VMAX*R*RR/(R*R+RR*RR)
      TTHETA=2.8E-3*RHOAIR*(VTHETA**2)
      TX(I,J)=-TTHETA*(Y-YC)/RR
   54 CONTINUE
   51 CONTINUE
      DO 52 I=2,M1,2
      X=FLOAT(I-1)*DX
      DO 52 J=2,N2,2
      Y=FLOAT(J-1)*DY
c      IF(MAP(I,J).EQ.0)GO TO 52
      RR=SQRT((X-XC)**2+(Y-YC)**2)
      IF(RR.LT.10.0)GO TO 53
      VTHETA=2.0*VMAX*R*RR/(R*R+RR*RR)
      TTHETA=2.8E-3*RHOAIR*(VTHETA**2)
      TY(I,J)=TTHETA*(X-XC)/RR
   53 CONTINUE
   52 CONTINUE
C     UPDATE X-VELOCITY
      DO  8 I=3,M2,2
      DO  8 J=3,N1,2
      IF(MAP(I,J).EQ.0)GOTO8
      X1=VEL(I-2,J)*FLUX1(I-2,J)
      IF(I.EQ.M1)X2=(2.0*VEL(M1,J)-VEL(M-3,J))*(2.0*FLUX1(M1,J)
     *-FLUX1(M-3,J))
      IF(I.NE.M1)X2=VEL(I+2,J)*FLUX1(I+2,J)
      X=(X2-X1)*DX4
      X1=0.5*(VEL(I-1,J-1)+VEL(I+1,J-1))*0.5*(FLUX1(I,J)+FLUX1(I,J-2))
      X2=0.0
      IF(J.NE.N1)X2=0.5*(VEL(I-1,J+1)+VEL(I+1,J+1))*0.5*(FLUX1(I,J)+
     *FLUX1(I,J+2))
      X=X+(X2-X1)*DY2
      X=X-F(I,J)*0.25*(FLUX1(I+1,J+1)+FLUX1(I+1,J-1)
     2+FLUX1(I-1,J+1)+FLUX1(I-1,J-1))
      X=-X-G*(ZETA2(I,J)+H(I,J))*(ZETA2(I+1,J)-ZETA2(I-1,J))*DX2
      X=X+TX(I,J)*RHO1
      X1=0.25*(FLUX1(I+1,J+1)+FLUX1(I+1,J-1)+
     2FLUX1(I-1,J+1)+FLUX1(I-1,J-1))
      X3=SQRT(FLUX1(I,J)**2+X1**2)*CF/(ZETA2(I,J)+H(I,J))**2
      FLUX2(I,J)=(FLUX1(I,J)+X*DT)/(1.0+X3*DT)
      VEL(I,J)=FLUX2(I,J)/(ZETA2(I,J)+H(I,J))
8     CONTINUE
      DO 10 I=1,M2,2
      FLUX2(I,1)=2.0*FLUX2(I,3)-FLUX2(I,5)
      VEL(I,1)=FLUX2(I,1)/(ZETA2(I,1)+H(I,1))
 10   CONTINUE
C     X-VELOCITY UPDAEED
C     UPDATE Y-VELOCITY
      DO  12 I=2,M1,2
      DO  12 J=2,N2,2
      IF(MAP(I,J).EQ.0)GOTO12
      X1=0.0
      X2=0.0
      IF(I.NE.2)X1=0.5*(VEL(I-1,J+1)+VEL(I-1,J-1))*0.5*(FLUX1(I,J)+
     2FLUX1(I-2,J))
      X2=0.5*(VEL(I+1,J+1)+VEL(I+1,J-1))*0.5*(FLUX1(I,J)+FLUX1(I+2,J))
      X=(X2-X1)*DX2
      IF(J.EQ.2)X1=(2.0*VEL(I,2)-VEL(I,4))*(2.0*FLUX1(I,2)-FLUX1(I,4))
      IF(J.NE.2)X1=VEL(I,J-2)*FLUX1(I,J-2)
      X=X+(VEL(I,J+2)*FLUX1(I,J+2)-X1)*DY4
      X=X+F(I,J)*0.25*(FLUX2(I+1,J+1)+FLUX2(I+1,J-1)+
     2FLUX2(I-1,J+1)+FLUX2(I-1,J-1))
      X=-X-G*(ZETA2(I,J)+H(I,J))*(ZETA2(I,J+1)-ZETA2(I,J-1))*DY2
      X=X+TY(I,J)*RHO1
      X1=0.25*(FLUX1(I+1,J+1)+FLUX1(I+1,J-1)+FLUX1(I-1,J+1)+
     2FLUX1(I-1,J-1))
      X3=SQRT(X1**2+FLUX1(I,J)**2)*CF/(ZETA2(I,J)+H(I,J))**2
      FLUX2(I,J)=(FLUX1(I,J)+X*DT)/(1.0+X3*DT)
C     IF(FLUX2(I,J).LT.RV)FLUX2(I,J)=0.0
      VEL(I,J)=FLUX2(I,J)/(ZETA2(I,J)+H(I,J))
12    CONTINUE
      DO 156 J=2,N2,2
      FLUX2(M,J)=2.0*FLUX2(M2,J)-FLUX2(M-4,J)
      VEL(M,J)=FLUX2(M,J)/(ZETA2(M,J)+H(M,J))
 156  CONTINUE
      KOUNT=KOUNT+1
C     Y-VELOCITY UPDATED
C        DO 216 J=1,N
C        DO 216 I=1,M
C        XUC(I,J)=0.0
C        YVC(I,J)=0.0
C        UC(I,J)=0.0
C        VC(I,J)=0.0
C        CUR(I,J)=0.0
C        DRN(I,J)=0.0
C        WTX(I,J)=0.0
C        WTY(I,J)=0.0
C216     CONTINUE
C        DO 499 I=1,M,2
C        DO 499 J=1,N1,2
C        XUC(I,J)=VEL(I,J)
C        WTX(I,J)=TX(I,J)
C499     UC(I,J)=VEL(I,J)
C        DO 498 I=2,M1,2
C        DO 498 J=2,N,2
C        YVC(I,J)=VEL(I,J)
C        WTY(I,J)=TY(I,J)
C498     VC(I,J)=VEL(I,J)
C        DO 800 I=2,M1,2
C        DO 800 J=1,N1,2
C        WTX(I,J)=0.5*(WTX(I+1,J)+WTX(I-1,J))
C800     UC(I,J)=0.5*(UC(I+1,J)+UC(I-1,J))
C        DO 801 I=1,M,2
C        DO 801  J=2,N2,2
C        WTX(I,J)=0.5*(WTX(I,J+1)+WTX(I,J-1))
C801    UC(I,J)=0.5*(UC(I,J+1)+UC(I,J-1))
C        DO 802 I=2,M1,2
C        DO 802 J=2,N2,2
C        WTX(I,J)=0.25*(WTX(I+1,J+1)+WTX(I+1,J-1)+WTX(I-1,J+1)+
C     *  WTX(I-1,J-1))
C        UC(I,J)=0.25*(UC(I+1,J+1)+UC(I+1,J-1)+UC(I-1,J+1)+UC(I-1,J-1))
C802     CONTINUE
C        DO 803 I=3,M2,2
C        DO 803 J=2,N2,2
C        WTY(I,J)=0.5*(WTY(I+1,J)+WTY(I-1,J))
C803     VC(I,J)=0.5*(VC(I+1,J)+VC(I-1,J))
C        DO 504 I=2,M1,2
C        DO 504  J=3,N1,2
C        WTY(I,J)=0.5*(WTY(I,J+1)+WTY(I,J-1))
C504     VC(I,J)=0.5*(VC(I,J+1)+VC(I,J-1))
C        DO 505 I=3,M2,2
C        DO 505 J=3,N1,2
C        WTY(I,J)=0.25*(WTY(I+1,J+1)+WTY(I+1,J-1)+WTY(I-1,J+1)+
C     *  WTY(I-1,J-1))
C        VC(I,J)=0.25*(VC(I+1,J+1)+VC(I+1,J-1)+VC(I-1,J+1)+VC(I-1,J-1))
C505     CONTINUE
C        DO 506 I=1,M
C        DO 506 J=1,N 
C        WCUR(I,J)=SQRT(WTX(I,J)**2+WTY(I,J)**2)
C        CUR(I,J)=SQRT(UC(I,J)**2+VC(I,J)**2)
C        IF(UC(I,J).EQ.0.0.AND.VC(I,J).EQ.0.0) GO TO 506
C        IF(WTX(I,J).EQ.0.0.AND.WTY(I,J).EQ.0.0) GO TO 506
C        DRN(I,J)=ATAN2(VC(I,J),UC(I,J))*57.29578
C        WDRN(I,J)=ATAN2(WTY(I,J),WTX(I,J))*57.29578
c       DRN(I,J)=DRN(I,J)+180.0
506     CONTINUE



      do 4786 j=1,n1
      do 4786 i=1,m1
      if (zeta2(i,j).gt.pz3(i,j))pz3(i,j)=zeta2(i,j)
4786  continue
C     PEAK SURGE VALUES AT FIRST OFF-SHORE POINTS
      DO 3734 J=1,N1
      DO 3734 I=1,M1
      IF(MAP(I-1,J).EQ.0.AND.MAP(I,J).EQ.1)GOTO4040
      IF(MAP(I,J).EQ.1.AND.MAP(I+1,J).EQ.0)GOTO4040
      IF(MAP(I,J).EQ.1.AND.MAP(I,J+1).EQ.0)GOTO4040
      GOTO3734
 4040 IF(ZETA2(I,J).GT.PZ1(I,J))GOTO3735
      GOTO3734
 3735 PZ1(I,J)=ZETA2(I,J)
      KO1(I,J)=KOUNT*DT/3600
      MIN1(I,J)=(KOUNT*DT-KO1(I,J)*3600)/60
 3734 CONTINUE
C     PEAK SURGE VALUES AT BOUNDARY
      DO 5734 J=1,N
      DO 5734 I=2,M1
      IF(MAP(I,J).EQ.0.AND.MAP(I-1,J).EQ.1)GOTO6734
      IF(MAP(I,J).EQ.0.AND.MAP(I+1,J).EQ.1)GOTO6734
      GOTO5734
 6734 IF(ZETA2(I,J).GT.PZ2(I,J))GOTO6735
      GOTO5734
 6735 PZ2(I,J)=ZETA2(I,J)
      KO2(I,J)=KOUNT*DT/3600
      MIN2(I,J)=(KOUNT*DT-KO2(I,J)*3600)/60
 5734 CONTINUE
C      WRITE(*,899)KOUNT
C899   FORMAT(5X,'KOUNT=',I5)
      DO  15 I=1,M
      DO  15 J=1,N
      FLUX1(I,J)=FLUX2(I,J)
      ZETA1(I,J)=ZETA2(I,J)
15    CONTINUE
      IF(KOUNT.NE.KOUNT/K1*K1.AND.KOUNT.NE.KS)GOTO 758
      TIME=FLOAT(KOUNT)*DT/3600.0
202   FORMAT(1H2,'KOUNT=',I5,2x,'XC=',2x,F10.5,2X,'YC=',2X,F10.5,/)
C220   FORMAT(1H ,'TIME IN HOURS',F10.2)
C      WRITE(6,220)TIME
C     WRITE(6,215)
C215  FORMAT(' ','SURFACE ELEVATION IN METRES')

758   CONTINUE
c      IF(KOUNT.EQ.50)WRITE(*,*)'NUMBER OF ITERATIONS 
c     * COMPLETED = ',KOUNT
c      IF(KOUNT.EQ.250)WRITE(*,*)'NUMBER OF ITERATIONS 
c     * COMPLETED = ',KOUNT
c      IF(KOUNT.EQ.500)WRITE(*,*)'NUMBER OF ITERATIONS 
c     * COMPLETED= ',KOUNT
c      IF(KOUNT.EQ.750)WRITE(*,*)'NUMBER OF ITERATIONS 
c     * COMPLETED = ',KOUNT
c      IF(KOUNT.EQ.1000)WRITE(*,*)'NUMBER OF ITERATIONS 
c     * COMPLETED = ',KOUNT
c      IF(KOUNT.EQ.1250)WRITE(*,*)'NUMBER OF ITERATIONS 
c     * COMPLETED = ',KOUNT
c      IF(KOUNT.EQ.1750)WRITE(*,*)'NUMBER OF ITERATIONS 
c     * COMPLETED = ',KOUNT
c      IF(KOUNT.EQ.KS)WRITE(*,*)'NUMBER OF ITERATIONS 
c     * COMPLETED = ',KOUNT
      IF(KOUNT.LT.KS)GOTO 9999
          DO 320 J=1,72
320       TRX(1,J)=72.0
         DO 321 J=1,72
         DO 321 I=2,53
321      TRX(I,J)=TRX(I-1,J)+14.0/FLOAT(52)
         DO 322 I=1,53
 322     TRY(I,1)=20.0
         DO 323 I=1,53
         DO 323 J=2,72
 323     TRY(I,J)=TRY(1,J-1)-18.0/FLOAT(71)

      DO 2322 J=1,N
      DO 2322 I=1,M
      JJ=N+1-J
      WRITE(11,*)TRX(I,J),TRY(I,J),ZETA2(I,JJ)
C      WRITE(12,2394)TRX(I,J),TRY(I,J),PZ1(I,JJ)
C      WRITE(13,2395)TRX(I,J),TRY(I,J),CUR(I,JJ),DRN(I,JJ)
C      WRITE(15,2395)TRX(I,J),TRY(I,J),WCUR(I,JJ),WDRN(I,JJ)
      write(18,2394)trx(i,j),try(i,j),pz3(i,jj)
2394  FORMAT(2F5.1,F10.5)        
2395  FORMAT(2F5.1,2F15.5)
2322  CONTINUE      
C      DO 2387 J=2,N-1
C      DO 2387 I=2,M-1
C      JJ=N+1-J
C      IF(MAP(I,JJ).EQ.0) GO TO 2387
C      WRITE(17,2395)TRX(I,J),TRY(I,J),WCUR(I,JJ),WDRN(I,JJ)
C      WRITE(16,2395)TRX(I,J),TRY(I,J),CUR(I,JJ),DRN(I,JJ)
C2387  CONTINUE      
      WRITE(6,3738)
      DO 3750 J=1,N1
      DO 3740 I=1,M1
      IF(MAP(I-1,J).EQ.0.AND.MAP(I,J).EQ.1)GOTO4141
      IF(MAP(I,J).EQ.1.AND.MAP(I+1,J).EQ.0)GOTO4141
      IF(MAP(I,J).EQ.1.AND.MAP(I,J+1).EQ.0)GOTO4141
      GOTO3740
 4141 WRITE(6,3739)J,I,KO1(I,J),MIN1(I,J),PZ1(I,J)
 3740 CONTINUE
      WRITE(6,3770)
 3750 CONTINUE
 3738 FORMAT(12X,'PEAK SURGE ENVELOPE AT FIRST OFF-SHORE POINT'//
     +11X,'J-GRID',3X,'I-GRID',3X,'HOUR',3X,'MIN',5X,
     +'ELEVATIONS'//)
      WRITE(6,6738)
      DO 6750 J=1,N
      DO 6740 I=2,M1
      IF(MAP(I,J).EQ.0.AND.MAP(I-1,J).EQ.1)GOTO6161
      IF(MAP(I,J).EQ.0.AND.MAP(I+1,J).EQ.1)GOTO6161
      GOTO6740
 6161 WRITE(6,3739)J,I,KO2(I,J),MIN2(I,J),PZ2(I,J)
 6740 CONTINUE
      WRITE(6,3770)
 6750 CONTINUE
 6738 FORMAT(////12X,'PEAK SURGE ENVELOPE AT BOUNDARY'//
     +11X,'J-GRID',3X,'I-GRID',3X,'HOUR',3X,'MIN',5X,
     +'ELEVATIONS'//)
 3739 FORMAT(12X,I2,7X,I2,7X,I2,4X,I2,8X,F6.3)
 3770 FORMAT(2X,'***********************************************')
      CLOSE(4)
      CLOSE(5)
      CLOSE(6)
c        IIGO=IIGO+1
c        IF(ANS.EQ.2.AND.IIGO.LE.3)GO TO 4444        
c        CALL CS(15)
c        WRITE(*,*)'           MODEL RUN SUCCESSFULLY'        
      STOP
      END
C     *************************************************************
C                         MAIN PROGRAM ENDS      
C     *************************************************************

      FUNCTION X0(T,K)
      DIMENSION ALAT1(20),ALONG1(20),DIFKOUNT(20),DTOTAL(20)
      COMMON /D/ALAT1,ALONG1,DIFKOUNT,DTOTAL
      ALAT0=2.0*0.0174
      ALONG0=72.0*0.0174
      T1=T-DTOTAL(K)
      ALAT=(ALAT1(K)+(ALAT1(K+1)-ALAT1(K))/DIFKOUNT(K+1)*T1)*0.0174
      ALONG=(ALONG1(K)+(ALONG1(K+1)-ALONG1(K))/DIFKOUNT(K+1)*T1)*0.0174
C     WRITE(1,751)ALAT,ALONG,KOUNT
C751  FORMAT(3X,E15.8,3X,E15.8,3X,I5)
      X0=6371000.0*(ALONG-ALONG0)*COS(ALAT)
      RETURN
      END
      FUNCTION Y0(T,K)
      DIMENSION ALAT1(20),ALONG1(20),DIFKOUNT(20),DTOTAL(20)
      COMMON /D/ALAT1,ALONG1,DIFKOUNT,DTOTAL
      ALAT0=2.0*0.0174
      T1=T-DTOTAL(K)
      ALAT=(ALAT1(K)+(ALAT1(K+1)-ALAT1(K))/DIFKOUNT(K+1)*T1)*0.0174
      Y0=6371000.0*(ALAT-ALAT0)
      RETURN
      END
      SUBROUTINE CS(N)
      INTEGER I1
      INTEGER N
      DO 10 I1=1,N
         WRITE(*,*)
   10 CONTINUE
      RETURN  
      END