c     THIS IS THE LATEST VERSION OF IIT-D MENU MODE STORM SURGE MODEL
C     LAST UPDATED 16.6.94
      PROGRAM LOGCOASTAL
C     ANDHRA CYCLONE MAY 1990
C     COASTAL ZONE MODEL WITH NEW DEPTH FN. AND LOG IN TRANSFORM
C     WITH MOVING STORM MOTION
      CHARACTER *30OUTFILE
      CHARACTER *50NAME
      INTEGER ANS       
      DIMENSION VEL(20,49),V1(20,49),V2(20,49),ZETA1(20,49),
     1ZETA2(20,49),F(20,49),U(20,49),H(20,49),B1(100),B2(100),
     2B(100),TX(20,49),TY(20,49),XI(42)
      DIMENSION PZ1(20,49),PZ2(100),KO1(100),KO2(100),MIN1(100)
     *,MIN2(100),KK1(402),KK2(100),YY(20,49)
      DIMENSION A(20,49),C(20,49)
      DIMENSION ALAT1(20),ALONG1(20),DIFKOUNT(20),DTOTAL(20),AHOUR(20)
      DIMENSION D(200),FANG(100),US(20),DIR(20)
      DIMENSION DIST(49),PEAK1(402),XDIST1(402),DIST1(49)
      REAL MAP(20,49),MAP1(20,49)
      COMMON /A/AL
      COMMON /B/ALAT1,ALONG1,DIFKOUNT,DTOTAL
      COMMON /C/DP,R,VMAX,RHOAIR,
C    *******************************************************
C     THE FOLLOWING PARAMETERS ARE TO BE SUPPLIED FOR ANY CYCLONE
C
C       : OF THE CYCLONE;MAX. 30 CHARACTERS IN LENGTH
C     KS    :TOTAL NO. OF ITERATIONS FOR EXECUTION;
C            KS X DT GIVES TOTAL NO. OF HOURS OF INTEGRATION
C     K1    :GIVE K1=KS
C     ND    :TOTAL NO. OF OBSERVATIONS OF THE CYCLONE TRACK
C     IC    :LANDFALL POINT
C     R     :RADIOUS OF MAXIMUM WINDS OF THE CYCLONE
C     DP    :PRESURE DROP OF THE CYCLONE
C     ALAT1 :LATITUDE OF EACH OBSERVATION OF THE CYCLONE TRACK
C            IN DEGREES
C     ALONG1:LONGITUDE OF EACH OBSERVATION OF THE CYCLONE TRACK
C            IN DEGREES
C     AHOUR :TIME DURATION IN HOURS BETWEEN TWO OBSERVATIONS OF
C            THE CYCLONE TRACK
C
C     END OF DATA OF PARAMETERS SUPPLIED
C    ********************************************************
C
C    THE FOLLOWING PARAMETERS ARE EXPLAINED
C     M     :NO. OF GRID POINTS IN THE X-DIR.
C     N     :NO. OF GRID POINTS IN THE Y-DIR.
C     INITIA:PARAMETER USED FOR CONTINUATION OF INTEGRATION
C     B1    :DISTANCE MEASURED FROM ORIGIN TO THE WEST COAST
C            OF THE ANALYSIS AREA
C     B2    :DISTANCE MEASURED FROM ORIGIN TO THE EAST COAST
C            OF THE ANALYSIS AREA
C     H     :ACTUAL DEPTH PROVIDED FOR THE ANALYSIS AREA
C     DIST  :DISTANCE MEASURED BETWEEN TWO GRID POINTS IN THE
C            Y-DIR. ALONG THE COAST
C
      OPEN(5,FILE='DATA.DAT')
      OPEN(10,FILE='INP.DAT')
      OPEN(2,FILE='PEAK.DAT')
C      
C*************************READ INPUT*****************************
C      
C      CALL INPUT(KS,DP,NAME)
C      CALL CS(23)
C        WRITE(*,*)'ENTER NAME OF OUT PUT FILE FOR RESULT:>'
C        READ(*,1367)OUTFILE
1367    FORMAT(A30)        
        OPEN(6,FILE='RESULT')      
      READ(5,101)M,N,K1,INITIA
101   FORMAT(I4,3(/I4))
      READ(5,102)AL,STRESS,RHO,DT,CIO,EPS
102   FORMAT(F10.0,5(/F10.0))
C*************************READ B1 AND B2*************************
      DO  2  JJ=1,N
      J=N+1-JJ
      READ(5,104)B1(J),B2(J)
      B(J)=2.5
      B2(J)=B1(J)+B(J)
      B1(J)=1.0E5*B1(J)
      B2(J)=1.0E5*B2(J)
      B(J) =1.0E5*B(J)
2     CONTINUE
104   FORMAT(F4.2,F4.2)
C
C**************************READ DEPTH****************************
C
      DO 444 JJ=1,N
      J=N+1-JJ
      READ(5,5011)(H(I,J),I=1,10)
444   CONTINUE
5011  FORMAT(10F3.1)
      DO 5022 I=11,20
      DO 5022 J=1,N
5022  H(I,J)=50.0
      DO 5033 I=1,M
      DO 5033 J=1,N
5033  H(I,J)=H(I,J)*1.83
      DO  5034  J=1,N
      DO  5034  I=1,2
5034  H(I,J)=H(I,J)*0.8
      DO  5035  J=1,N
      DO  5035  I=3,7
5035  H(I,J)=H(I,J)*0.7

C
C*************************READ DISTANCES***************************
C
      II=1
      IM=7
      DO 905 J=1,7
      READ(5,906)(DIST(I),I=II,IM)
      II=IM+1
      IM=IM+7
 905  CONTINUE
 906  FORMAT(7F3.0)
      READ(10,*)NAME
      READ(10,*)DP
      READ(10,*)KS
      READ(10,*)IC
      READ(10,*)R
      READ(10,*)ND
      DO 9089 I=1,ND
      READ(10,*)ALAT1(I),ALONG1(I),AHOUR(I)
9089  CONTINUE 
C
C******************************************************************    
C*********************DIFFERENT SCENARIOS**************************  
C
C      CALL CS(23)
C      WRITE(*,*)'1. SCENARIO(1):TO RUN THE PROGRAM FOR PREDICTED' 
C      WRITE(*,*)'   LAND FALL POINT'
C      CALL CS(1)
C      WRITE(*,*)'2. SCENARIO(2):TO RUN THE PROGRAM FOR LANDFALL' 
C      WRITE(*,*)'   100 KM LEFT OF PREDICTED LAND FALL POINT'
C      CALL CS(1)
C      WRITE(*,*)'3. SCENARIO(3):TO RUN THE PROGRAM FOR LANDFALL '
C      WRITE(*,*)'   100 KM RIGHT OF PREDICTED LAND FALL POINT'
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 PREDICTED LAND FALL POINT:>'
C      READ(*,*)IC
C      GO TO 3336        
C3334  WRITE(*,*)'ENTER LAND FALL POINT 100 KM LEFT OF'
C      WRITE(*,*)'PREDICTED LAND FALL POINT:>'
C      READ(*,*)IC
C      GO TO 3336
C3335  WRITE(*,*)'ENTER LAND FALL POINT 100 KM RIGHT OF'
C      WRITE(*,*)'PREDICTED LAND FALL POINT:>'   
C      READ(*,*)IC
C3336  CALL CS(23)
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)
C      DO 9019 K=1,ND
C      READ(*,*)ALAT1(K),ALONG1(K),AHOUR(K)
Cc1312  FORMAT(3F6.2)      
C9019  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(6,8677)NAME
8677  FORMAT('NAME OF THE CYCLONE=',A30)
      WRITE(6,205)AL,STRESS,RHO,DT,R,DP
205   FORMAT('NORTH-SOUTH EXTENT OF ANALYSIS AREA IN METRES=',
     *3X,F10.0/'LATITUDE OF OPEN-SEA BOUNDARY IN RADIAN=',7X,F10.4/'DENS
     *ITY OF WATER IN GM/M3=',22X,F10.0/
     *'TIME-STEP IN SEC=',31X,F10.0/
     * 'E-FOLDING RADIUS OF CYCLONE IN METERS=',13X,F10.0/
     *'PRESSURE-DROP IN CYCLONE IN MB=',17X,F10.0)
      WRITE(6,8678)KS,IC,ND
8678  FORMAT('NO OF ITERAIONS= ',I4,/,
     *'LAND FALL POINT=',I4,/,
     *'NO OF CYCLONE TRACK POSITIONS=',I3,/,
     *'LATITUDE',3X,'LONGITUDE',/)
      WRITE(6,8679)(ALAT1(K),ALONG1(K),K=1,ND)      
8679  FORMAT(1X,F6.2,6X,F6.2)
C
C*******************************************************************
C***********************MAIN PROGRAM BEGINS*************************
      M1=M-1
      M2=M-2
      M3=M-3
      N1=N-1
      N2=N-2
      N3=N-3
      ND1=ND-1
      CF=2.6E-3
      RHOAIR=1.226E3
      ETAM=1.0+EPS*ALOG(1.0+1.0/CIO)
      DETA=ETAM/FLOAT(M1)
      DY=AL/FLOAT(N-1)
      CALL TRANS(XI,DETA,CIO,EPS,M)
      DO 3 J=1,N
      SAA=DY/6371220.0
      S=2.0*7.292E-5*SIN(STRESS+FLOAT(J-1)*SAA)
      DO 3I=1,M
      PZ1(I,J)=0.0
      PZ2(J)=0.0
      KO1(J)=0
      KO2(J)=0
      MIN1(J)=0
      MIN2(J)=0
      VEL(I,J)=0.0
      V1(I,J)=0.0
      V2(I,J)=0.0
      ZETA1(I,J)=0.0
      ZETA2(I,J)=0.0
      U(I,J)=0.0
      F(I,J)=S
3     CONTINUE

      DO  50 J=1,N
      DO  50 I=1,M
      MAP(I,J)=B1(J)+B(J)*XI(I)
50    MAP1(I,J)=MAP(I,J)
C     WRITE(6,225)
225   FORMAT(1H1,'COORDINATE POINTS IN METRES')
      SUM=0.0
      DO 907 K=1,ND
      DIFKOUNT(K)=AHOUR(K)*20.0
      DTOTAL(K)=SUM+DIFKOUNT(K)
      SUM=SUM+DIFKOUNT(K)
907   CONTINUE
      SUM=0.0
      DO 904 J=1,N
      SUM=DIST(J)+SUM
      DIST1(J)=SUM
904   CONTINUE
      RAD=1.745329E-2
      DO  99  K=1,ND1
      D1=ALAT1(K+1)-ALAT1(K)
      IF(D1.EQ.0.0)D1=0.1
      D2=ALONG1(K+1)-ALONG1(K)
      R1=SQRT(D1**2+D2**2)
      HR=(DTOTAL(K+1)-DTOTAL(K))/20
      US(K)=R1*100/HR*0.278
      DIR(K)=ASIN(D1/R1)/RAD
      IF(D2.LT.0.0)DIR(K)=180.0-DIR(K)
99    CONTINUE
      K=1
      G=9.78
      KOUNT=0
      VMAX=0.0
      CALL FORCE1(D,FANG)
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,M2,2
      DO  4 J=3,N2,2
      X1=0.0
      IF(I.NE.2)X1=(ZETA1(I-1,J)+H(I-1,J))*U(I-1,J)
      X1=B(J)*X1
      X2=(ZETA1(I+1,J)+H(I+1,J))*U(I+1,J)
      X2=B(J)*X2
      X=(X2-X1)/(2.0*DETA)
      X=X*(1.0+EPS/(XI(I)+CIO))
      X1=V1(I,J-1)
      X2=2.0*V1(I,N1)-V1(I,N3)
      IF(J.NE.N)X2=V1(I,J+1)
      X=X+(X2-X1)/(2.0*DY)
      X=-X
      ZETA2(I,J)=ZETA1(I,J)+DT*X/B(J)
4     CONTINUE
C     RADIATION CONDITION AT EASTERN BOUNDARY
      DO  155 J=3,N2,2
155   ZETA2(M,J)=-ZETA2(M2,J)+2.0*SQRT(H(M1,J)/G)*U(M1,J)*B(J)
C     RADIATION CONDITION AT SOUTHETN BOUNDARY
      DO 153 I=2,M,2
153   ZETA2(I,1)=-ZETA2(I,3)-2.0*SQRT(H(I,2)/G)*VEL(I,2)
C     RADIATION CONDITION AT NORTHERN BOUNDARY
      DO 1531 I=2,M,2
1531  ZETA2(I,N)=-ZETA2(I,N2)+2.0*SQRT(H(I,N1)/G)*VEL(I,N1)
      DO  5 I=3,M1,2
      DO  5 J=3,N,2
5     ZETA2(I,J)=0.5*(ZETA2(I-1,J)+ZETA2(I+1,J))
      DO  20 I=2,M,2
      DO  20 J=2,N1,2
20    ZETA2(I,J)=0.5*(ZETA2(I,J+1)+ZETA2(I,J-1))
      DO 6 J=1,N
      ZETA2(1,J)=0.5*(3.0*ZETA2(2,J)-ZETA2(4,J))
6     CONTINUE
C     ELEVATION UPDATED
C     UPDATE X-VELOCITY
      IF(FLOAT(KOUNT).EQ.DTOTAL(K+1))K=K+1
      XC=X0(FLOAT(KOUNT),K)
      YC=Y0(FLOAT(KOUNT),K)
      CALL STRES(D,FANG,MAP,TX,TY,XC,YC,US(K),DIR(K),M,N)
      DO  8 I=3,M1,2
      DO  8 J=3,N2,2
      X1=0.0
      IF(I.NE.3)X1=U(I-2,J)*V1(I-2,J)
      IF(I.EQ.M1)X2=(2.0*U(M1,J)-U(M3,J))*(2.0*V1(M1,J)-V1(M3,J))
      IF(I.NE.M1)X2=U(I+2,J)*V1(I+2,J)
      X=(X2-X1)/(4.0*DETA)
      X=X*(1.0+EPS/(XI(I)+CIO))
      X1=0.5*(VEL(I-1,J-1)+VEL(I+1,J-1))*0.5*(V1(I,J)+V1(I,J-2))
      X2=0.5*(VEL(I-1,J+1)+VEL(I+1,J+1))*0.5*(V1(I,J)+V1(I,J+2))
      X=X+(X2-X1)/(2.0*DY)
      X=X-F(I,J)*0.25*(V1(I-1,J-1)+V1(I-1,J+1)+V1(I+1,J-1)+V1(I+1,J+1))
      X1=(ZETA2(I+1,J)-ZETA2(I-1,J))/(2.0*DETA)
      X1=X1*(1.0+EPS/(XI(I)+CIO))
      X2=(ZETA2(I,J)+H(I,J))*G
      X=-X-X1*X2
      X1=0.25*(VEL(I-1,J-1)+VEL(I-1,J+1)+VEL(I+1,J-1)+VEL(I+1,J+1))
      X=X+B(J)*TX(I,J)/RHO
      X3=SQRT(VEL(I,J)**2+X1**2)*CF/(ZETA2(I,J)+H(I,J))
      V2(I,J)=(V1(I,J)+X*DT)/(1.0+X3*DT)
8     CONTINUE
      DO  10 J=3,N,2
      V2(1,J)=2.0*V2(3,J)-V2(5,J)
10    CONTINUE
      DO  9 I=1,M1,2
      V2(I,N)=2.0*V2(I,N2)-V2(I,N-4)
9     V2(I,1)=2.0*V2(I,3)-V2(I,5)
      DO  11 I=1,M1,2
      DO  11 J=1,N,2
11    VEL(I,J)=V2(I,J)/(B(J)*(ZETA2(I,J)+H(I,J)))
C     X-VELOCITY UPDATED
C     UPDATE Y-VELOCITY
      DO  12 I=2,M2,2
      DO  12 J=2,N1,2
      X1=0.0
      IF(I.NE.2)X1=0.5*(U(I-1,J+1)+U(I-1,J-1))*0.5*(V1(I,J)+V1(I-2,J))
      X2=0.5*(U(I+1,J+1)+U(I+1,J-1))*0.5*(V1(I,J)+V1(I+2,J))
      X=(X2-X1)/(2.0*DETA)
      X=X*(1.0+EPS/(XI(I)+CIO))
      IF(J.EQ.2)X1=(2.0*VEL(I,2)-VEL(I,4))*(2.0*V1(I,2)-V1(I,4))
      IF(J.NE.2)X1=VEL(I,J-2)*V1(I,J-2)
      IF(J.EQ.N1)X2=(2.0*VEL(I,N1)-VEL(I,N3))*(2.0*V1(I,N1)-V1(I,N3))
      IF(J.NE.N1)X2=VEL(I,J+2)*V1(I,J+2)
      X=X+(X2-X1)/(4.0*DY)
      X=X+0.25*(V2(I-1,J-1)+V2(I-1,J+1)+V2(I+1,J-1)+V2(I+1,J+1))*F(I,J)
      X=-X-G*B(J)*(ZETA2(I,J)+H(I,J))*(ZETA2(I,J+1)-ZETA2(I,J-1))/
     1(2.0*DY)
      X1=G*(ZETA2(I,J)+H(I,J))
      X2=(B1(J+1)-B1(J-1))/(2.0*DY)
      X3=  (B(J+1)-B(J-1))/(2.0*DY)
      X1=X1*(X2+X3*XI(I))
      X1=X1*(1.0+EPS/(XI(I)+CIO))
      IF(I.EQ.2)X4=(0.5*(ZETA2(4,J+1)+ZETA2(4,J-1))
     2-0.5*(ZETA2(2,J+1)+ZETA2(2,J-1)))/(2.0*DETA)
      IF(I.NE.2)X4=(0.5*(ZETA2(I+2,J+1)+ZETA2(I+2,J-1))-
     10.5*(ZETA2(I-2,J+1)+ZETA2(I-2,J-1)))/(4.0*DETA)
      X=X+X1*X4
      X=X+B(J)*TY(I,J)/RHO
      X1=VEL(I,J)
      X2=0.25*(VEL(I-1,J-1)+VEL(I-1,J+1)+VEL(I+1,J-1)+VEL(I+1,J+1))
      X3=SQRT(X1**2+X2**2)*CF/(ZETA2(I,J)+H(I,J))
      V2(I,J)=(V1(I,J)+X*DT)/(1.0+X3*DT)
12    CONTINUE
      DO  156 J=2,N1,2
156   V2(M,J)=2.0*V2(M2,J)-V2(M-4,J)
      DO  13 I=2,M,2
      DO  13 J=2,N1,2
13    VEL(I,J)=V2(I,J)/(B(J)*(ZETA2(I,J)+H(I,J)))
      DO  14 I=3,M1,2
      DO  14 J=3,N2,2
      X1=(B1(J+1)-B1(J-1))/(2.0*DY)+XI(I)*
     1(B(J+1)-B(J-1))/(2.0*DY)
      U(I,J)=(VEL(I,J)-X1*0.25*(VEL(I-1,J-1)+VEL(I-1,J+1)+VEL(I+1,J+1)
     1+VEL(I+1,J-1)))/B(J)
14    CONTINUE
      DO  141 I=3,M1,2
      U(I,N)=2.0*U(I,N2)-U(I,N-4)
141   CONTINUE
C     PEAK SURGE VALUES 
      DO  3734 J=1,N
      DO  3734 I=1,M
      IF(ZETA2(I,J).GT.PZ1(I,J))GOTO3735
      GOTO3736
3735  PZ1(I,J)=ZETA2(I,J)
      KO1(J)=KOUNT*DT/3600
      MIN1(J)=(KOUNT*DT-KO1(J)*3600)/60
3736  IF(ZETA2(2,J).GT.PZ2(J))GOTO3737
      GOTO3734
3737  PZ2(J)=ZETA2(2,J)
      KO2(J)=KOUNT*DT/3600
      MIN2(J)=(KOUNT*DT-KO2(J)*3600)/60
3734  CONTINUE
      KOUNT=KOUNT+1
      DO  15 I=1,M
      DO  15 J=1,N
      ZETA1(I,J)=ZETA2(I,J)
      V1(I,J)=V2(I,J)
15    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
CC      IF(KOUNT.EQ.606)CALL OUTPUT(ZETA2,M,N,CHNO)
      IF(KOUNT.LT.KS)GOTO9999
      WRITE(6,8680)
8680  FORMAT('PEAK SURGE ENVELOPE ALONG THE COAST',/,
     111X,'NO',3X,'HOUR',3X,'MINUTE',3X,'SURGE IN METERS')                 
      DO 2727 J=1,N
      WRITE(6,3739)J,KO1(J),MIN1(J),PZ1(1,J)
      KK2(J)=KO1(J)*60+MIN1(J)
2727  CONTINUE
3739  FORMAT(11X,I2,4X,I2,4X,I4,7X,F6.3)
C     TO CALCULATE PEAK SURGE FOR EVERY 5KM. DISTANCE ALONG
C     THE COAST
      DO  2017  I=1,402
      PEAK1(I)=0.0
2017  CONTINUE
      IM=2
      IN=0
      KK1(1)=KK2(1)
      PEAK1(1)=PZ1(1,1)
      DO 1000 J=1,N2,2
      II=(DIST1(J+2)-DIST1(J))/5.0
      II=II+IN
      I1=II-IN
      I2=II+1
      DO 1001 IJ=IM,I2
      PEAK1(IJ)=PZ1(1,J)*(I2-IJ)/I1+PZ1(1,J+2)*(I1-(I2-IJ))/I1
      KK1(IJ)=KK2(J)*(I2-IJ)/I1+KK2(J+2)*(I1-(I2-IJ))/I1
1001  CONTINUE
      IM=II+2
      IN=I1+IN
1000  CONTINUE
      DO 3344 JJ=1,402
      XDIST1(JJ)=JJ*5.0
      WRITE(2,*)XDIST1(JJ),PEAK1(JJ)
3344  CONTINUE
      CALL PLOT(PEAK1,XDIST1,IC,NAME)
C     *********INTERPOLATION OF PEAKSURGE VALUES************
        DO 9786 J=2,N1,2      
        DO 9786 I=3,M1,2
        PZ1(I,J)=0.25*(PZ1(I+1,J+1)+PZ1(I+1,J-1)+PZ1(I-1,J+1)
     *   +PZ1(I-1,J-1))
9786    CONTINUE      
        DO 9787 I=3,M1,2     
        PZ1(I,1)=0.5*(PZ1(I+1,1)+PZ1(I-1,1))      
9787    CONTINUE
C***************************************************************
        DELY1=AL/FLOAT(48)
        DO 37 I=1,20
37      YY(I,1)=0.0
        DO 5542 I=1,20
        DO 5542 J=2,49
5542    YY(I,J)=YY(I,J-1)+DELY1
C       DO 354 J=1,49
C       JJ=N+1-J
C       DO 354 I=1,20
C       MAP1(I,J)=MAP1(I,J)
C       WRITE(14,*)MAP1(I,JJ)/1000.0,YY(I,JJ),ZETA2(I,JJ)
C       WRITE(15,*)I,JJ,ZETA2(I,JJ)
C       WRITE(16,*)I,JJ,PZ1(I,JJ)
C      WRITE(17,*)MAP1(I,JJ)/1000.0,YY(I,JJ)/1000.0,PZ1(I,JJ)
354    CONTINUE       
        CALL CONTOUR(B1,B,XI,AL,PZ1)
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'        
C     TO TAKE PRINTING THE PEAK SURGE FOR EVERY 5KM. DISTANCE
C     ALONG THE COAST
C     CALL PLOT(402,PEAK,IC)
C     CALL COAST(B1,B2,N,XC,YC,DY,R)
3339    CONTINUE
      STOP
      END
      FUNCTION DEPTH (XI,Y,ALPHA,BETA,GAMMA,H0)
      COMMON /A/AL
      Y=Y/AL
      X=(XI**BETA)*((1.0-Y)**(Y**GAMMA))*ALPHA
      DEPTH=H0*EXP(X)
      RETURN
      END
      SUBROUTINE TRANS(XI,DETA,CIO,EPS,M)
      DIMENSION XI(42)
      XI(1)=0.0
      DO  2 I=2,M
      XI(I)=XI(I-1)
3     X=XI(I)
      Y=-FLOAT(I-1)*DETA+X+EPS*ALOG(1.0+X/CIO)
      Z=1.0+EPS/(X+CIO)
      XI(I)=X-Y/Z
      IF(ABS(X-XI(I)).GT.1.0E-3*ABS(X))GOTO3
2     CONTINUE
      RETURN
      END
      SUBROUTINE COAST(B1,B2,N,XC,YC,DY,RR)
      CHARACTER BB,P,Q,R,S
      DIMENSION B1(100),B2(100),BB(120)
      DATA P,Q,R,S/' ','*','C','$'/
      BMAX=0.0
      BMIN=0.0
      DO1J=1,N
      IF(B2(J).GT.BMAX)BMAX=B2(J)
      IF(B1(J).LT.BMIN)BMIN=B1(J)
1     CONTINUE
      B=BMAX-BMIN
      M3=1+INT(ABS(XC-BMIN)/B*119.0)
      II=INT(YC/DY)
      IK=0
      IR=INT(RR/DY)
      IPR=II+IR
      IMR=II-IR
      DO2JJ=1,N
      J=N+1-JJ
      M1=1+INT(ABS(B1(J)-BMIN)/B*119.0)
      M2=120-INT(ABS(B2(J)-BMAX)/B*119.0)
      DO3I=1,120
      BB(I)=P
      IF(I.EQ.M1.OR.I.EQ.M2)BB(I)=Q
      IF(J.GT.IPR.OR.J.LT.IMR)GOTO3
      IF(IPR.GT.N)IK=IPR-N
      J1=IR-IK
      T=FLOAT(J1)*DY/RR
      TH=ACOS(T)
      X1=XC+RR*SIN(TH)
      X2=XC-RR*SIN(TH)
      M4=1+INT(ABS(X1-BMIN)/B*119.0)
      IF(X2.LT.BMIN)GOTO6
      M5=1+INT(ABS(X2-BMIN)/B*119.0)
      BB(M5)=S
  6   BB(M4)=S
      IF(J.NE.II)GOTO3
      IF(XC.LT.BMIN.OR.II.GT.N)GOTO3
      BB(M3)=R
3     CONTINUE
      IF(J.LE.IPR)IK=IK+1
      IF(J.NE.N)GOTO4
      DO5I=M1,M2
5     BB(I)=Q
4     CONTINUE
      WRITE(6,201)J,(BB(I),I=1,120)
2     CONTINUE
201   FORMAT(1H ,I2,3X,120A1)
      RETURN
      END
      SUBROUTINE OUTPUT(Y,M,N,CHNO)
      INTEGER CHNO
      DIMENSION Y(20,49)
      MAX=126
      DO 1 J=1,N
      DO 1 I=1,M
      IF(Y(I,J).NE.0.0)GO TO 1000
 1    CONTINUE
      WRITE(6,100)
 100  FORMAT(/////'ALL ZERO',//)
      RETURN
 1000 BIGVAL=ABS(Y(1,1))
      DO 2 J=1,N
      DO 2 I=1,M
      BIGVAL=AMAX1(BIGVAL,ABS(Y(I,J)))
 2    CONTINUE
      RN=ALOG10(BIGVAL)
      IF(RN.GE.0)GO TO 9000
      IF(RN-INT(RN).EQ.0)GO TO 9000
      RN=RN-1
 9000 CONTINUE
      NN=INT(RN)
      WRITE(6,1004)NN
 1004 FORMAT(100X,'*10**',I3)
      SF=10.**NN
      IF(NN.EQ.0)GO TO 1001
      DO 3 J=1,N
      DO 3 I=1,M
      Y(I,J)=Y(I,J)/SF
 3    CONTINUE
 1001 CONTINUE
      IWDTH=(MAX-6)/6
      LOW=1
      IUPP=(LOW+IWDTH-1)
 1002 IF(IUPP.GT.M)IUPP=M
      WRITE(6,101)(I,I=LOW,IUPP)
 101  FORMAT(//,4X,20I6,//)
      DO 4 J=1,N
      JJ=N+1-J
      WRITE(6,102)JJ,(Y(I,JJ),I=LOW,IUPP)
 102  FORMAT(1H0,I3,2X,20F6.2)
 4    CONTINUE
      LOW=LOW+IWDTH
      IF(LOW.GT.M)GO TO 1003
      IUPP=IUPP+IWDTH
      WRITE(6,103)
103   FORMAT(1H2)
      GO TO 1002
 1003 IF(N.EQ.0)RETURN
      DO 5 J=1,N
      DO 5 I=1,M
      Y(I,J)=Y(I,J)*SF
 5    CONTINUE
      RETURN
      END
      FUNCTION X0(T,K)
      DIMENSION ALAT1(20),ALONG1(20),DIFKOUNT(20),DTOTAL(20)
      COMMON /B/ALAT1,ALONG1,DIFKOUNT,DTOTAL
      ALAT0=6.0*0.0174
      ALONG0=81.84*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
      X0=6371000.0*(ALONG-ALONG0)*COS(ALAT)
      RETURN
      END
      FUNCTION Y0(T,K)
      DIMENSION ALAT1(20),ALONG1(20),DIFKOUNT(20),DTOTAL(20)
      COMMON /B/ALAT1,ALONG1,DIFKOUNT,DTOTAL
      ALAT0=6.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 FORCE1(D,FANG)
C     JELESNIANSKI SEPT 1980TDL IBM 360/195
C     PURPOSE
C     THIS SUBROUTINE (FORCE1) COMPUTES,AT MILE INTERVALS FROM
C     THE STORM CENTER, FOR 800 MILES.
C     (1)SINE OF THE INFLOW ANGLE
C     (2)COSINE OF THE INFLOW ANGLE
C     (3)SURFACE PRESSURE
C     (4)STATIC HEIGHT
C     REFERENCE FOR THE ALGORITHM AND THEORY IS:
C     JELESNIANSKI AND TAYLOR, 1973, A PRELIMINARY
C     VIEW OF STORM SURGES BEFORE AND AFTER STORM
C     MODIFICATIONS, NOAA TECHNICAL MEMORANDUM
C     ERL WMPO-3,WEATHER MODIFICATION PROGRAM OFFICE,
C     ERL, NOAA, U.S.DEPARTMENT OF COMMERCE,WASHINGTON,
C     D.C.,MAY,33PP.
C     DATA SET USE
C     NONE
C     VARIABLES
C     VMAX=DUMMY VALUE FOR MAXIMUM WINDS
C     RMAX=DUMMY VALUE FOR STORM SIZE
C     XKS XKN=WIND FRICTION VALUES
C     CDRAG=DUMMY VALUE FOR WIND FRCTN CFFCNT,IN 'BLCKDT'
C     IC12=MAX DISTANCE FROM STRM CNTR,SET IN 'BLCKDT'
C     COR=CORIOLIS(1/SEC)
C     CORHR=CORIOLIS(1/HOUR)
C     C(800) S(800)=COS/SIN OF INFLOW ANGLE,MILE INTERVALS
C     P(800)=PRE GRADIENT,MILE INTERVALS
C     DELP(800)=STATIC HEIGHTS,MILE INTERVALS
C     RHOA=AIR DENSITY (MB HR**2 MILE **-2)
C     RHOW=WATER DENSITY (MB SEC**2 FT**-1 MILE**-1)
C     RHOWG=WATER DENSITY*GRAVITY(MB FT**-1)
C     PADD=PRE FROM A RADIUS TO INFINITY
C     PINFT=PRE DROP FROM ITERATION
C     PDROP=INPUT PRE DROP
C     X12(4)=STORM SIZE IN FEET
C     C24=DUMMY VALUE FOR STORM SIZE
C     X12(8)=MAX WINDS IN FT/SEC
C     X12(15)=SQUARE OF STORM SIZE
C     C19=DRAG COEFFICIENT
C     X12(17)=DRAG CFFCNT (3*10-6)TIMES SQUARE OF STRM SIZE
C     X12(18)=TWICE THE MAX WIND
C     ITIME=COUNTER FOR TIME STEPS
C     NOLTHR=NUMBER OF TIME STEPS IN ONE HOUR
C     IBGNT=SUBSCRIPT FOR STORM INITIALIZATION
C     JHR=SUBSCRIPT FOR NEAREST APPROACH OF STORM
C     GENERAL COMMENTS
C     THIS SUBROUTINE RESIDE IN OVERLAY `CMPUTE'.IT IS CALLED IN
C     BY TWO SUBROUTINES,`SETCMP'AND `STMVAL'.SUBROUTINE
C     'FORCE1'CALLS IN SUBROUTINE 'PRTSTM'.
      CHARACTER *30NAME
      COMMON /DUMMY4/ S(800),C(800),P(800),DELP(800)
      COMMON /C/PDROP,RMAX,VMAX,RHOAIR,NAME
      DIMENSION IBEG(3),IEND(3),IN(3)
      DIMENSION  VELR(20,49),VELRX(20,49),VELRY(20,49),THETA(20,49)
     1,MAP(20,49),FANG(100),D(200)
      EQUIVALENCE(C23,VMAX),(C24,RMAX)
      DATA IBEG/1,111,251/,IEND/101,231,751/,IN/2,10,50/
      DATA RAD/1.745329E-2/
      DATA IC12 /795/
C     INPUT LATITUDE OF THE STORM
      DATA ALTO/14./
C     INPUT PRE DROP AND STORM RADIUS
C
      DATA RHOA/2.298E-03/,RHOW/4905./,RHOWG/29.89/
      DATA EPSLN /1.E-5/
C     5    FUNCTIONAL FORMS FOLLOW
      VEEA(R,RMAX)=RMAX/(RMAX**2+R**2)
      VEE(R)=2.*VMAX*RMAX*R/(RMAX**2+R**2)
      DERLNV(R)=1./R-2.*R/(RMAX**2+R**2)
      PRGRAD(R,V)=RHOA*V*V*(XKS/SINPHI-DERLNV(R))
      ATANC(P1,Q1)=ATAN2(SIGN(AMAX1(1.0E-9,ABS(P1)),P1),Q1)
      PYTHAG(X)=SQRT((1.-X)*(1.+X))
      RMAX=RMAX*0.6214/1000.
      JSET=MIN0(IC12-1,MAX0(150,INT(15.*RMAX)))
      AL=ALTO*1.745529E-2
      COR=2.*(7.292116E-5)*SIN(AL)
      CORHR=3600.*COR
      CDRAG=SQRT(22./RMAX)
      CALL INTWND(PDROP,RMAX,VMAX,ALTO)
       VMAX=14.2*SQRT(PDROP)*1.152
C     BEGIN ITERATION ALGORITHM
 100  CONTINUE
C
C     FORM 2 FRICTION VALUES FOR A STORM
      XKS=.01*SQRT(RMAX/(.3*VMAX+60.))
      XKN=.01*SQRT(RMAX/(.4*VMAX+80.))
      XKS=XKS*CDRAG
      XKN=XKN*CDRAG
C     SET RANGES RELATIVE TO STORM CENTER,SCALED TO STORM SIZE
      J2=MIN1(3.*RMAX,50.)
      R0=JSET+50
      MMAX=JSET+9*J2+49
      DELR=1.
      VEER0=VEE(R0)*R0
      COTPHI=(CORHR*R0/VEER0+XKN)/XKS
      W0=VEER0*COTPHI/SQRT(1.+COTPHI*COTPHI)
C     DWDR=0.WITH THIS INITIAL CONDITION
      SINPHI=PYTHAG(W0/VEER0)
      PRESS=0.
      DPDR0=PRGRAD(R0,VEER0/R0)
              CALL INVALOP(L*4)
        CALL OVEFL(L*4)
        CALL DVCHK(L*4)      
C     RUNGE-KUTTA-HEUN-TAYLOR METHOD TO SOLVE THE EQUATION
C     DWDR=W*XKS/SIN(PHI)-CORHR*R-XKN*V*R=G(W,R)WHERE W=V*R*COS(PHI)
      DO 160 I=1,MMAX
      FACT=2./DELR
      SLP0=W0*XKS/SINPHI-CORHR*R0-XKN*VEER0
      R1=R0-DELR
      VEER1=VEE(R1)*R1
      W1=W0*(VEER1/VEER0)
C     NEWTON-RAPHSON METHOD TO FIND W1 SO THAT
C     .5*(G(W1,R1)+G(W0,R0))=(W0-W1)/(R0-R1)
      DO 120 K=1,75
      SINPHI=PYTHAG(W1/VEER1)
      SLP1=W1*XKS/SINPHI-CORHR*R1-XKN*VEER1
      PADD=FACT*(W0-W1)-SLP0
      DGDW=XKS/(SINPHI*SINPHI*SINPHI)
      DELW=(SLP1-PADD)/(DGDW+FACT)
      DELW=AMAX1(DELW,.9*(W1-VEER1))
      W1=W1-DELW
      IF(W1) 180,180,110
 110  IF(ABS(DELW/VEER1)-EPSLN) 140,140,120
 120  CONTINUE
      WRITE(6,130)
 130  FORMAT('0INSUFFICIENT STEPS IN NEWTON-RAPHSON LOOP')
      GO TO 999
 140  R0=R1
      K=R1+1.01
      VEER0=VEER1
      W0=W1
      C(K)=W0/VEER0
      SINPHI=PYTHAG(C(K))
      S(K)=SINPHI
      DPDR1=PRGRAD(R1,VEER1/R1)
      P(K)=DPDR1
      PRESS=PRESS-(DPDR0+DPDR1)/FACT
      DELP(K)=PRESS
      DPDR0=DPDR1
      IF(K-J2-1) 150,150,160
 150  DELR=.1
 160  CONTINUE
      DELP(1)=PRESS-DPDR0/FACT
      P(1)=0.
      S(1)=0.
      C(1)=1.
      R0=JSET-1
      TERM=R0/(RMAX*RMAX+R0*R0)
      PADD=2.*RHOA*((VMAX*RMAX)**2)*(XKS*(ATAN(RMAX/R0)/RMAX+TERM)/
     *S(JSET))+.5*RHOA*VEE(R0)**2
C     PADD=PRE IN MB FROM R0=JSET-1 TO R0=INFINITY,ASSUMING
C     SINPHI=S(JSET)=CONST. FROM R=JSET-1 TO  INFINITY
      PINFT=PADD+DELP(JSET)-DELP(1)
C     PINFT=PRE DROP IN MB FROM R=0 TO INFINITY,
C     THEN ON THIS ASSUMPTION
      PRS=PDROP
      GOTO200
180   WRITE(6,190)
190   FORMAT('0COS(PHI).LT.0')
      GOTO999
200   SINPHI=S(JSET-1)
      COSPHI=C(JSET-1)
      JSTP1=JSET+1
      DO 210 K=JSTP1,IC12
      R0=K+1
      S(K)=SINPHI
      C(K)=COSPHI
      P(K)=PRGRAD(R0,VEE(R0))
      DELP(K)=DELP(K-1)+.5*(P(K)+P(K-1))
 210  CONTINUE
      PRESS=PADD+DELP(JSET)
      DO 220 I=1,IC12
      DELP(I)=(PRESS-DELP(I))/RHOWG
      P(I)=P(I)/RHOW
 220  CONTINUE
      C233=VEE(RMAX)
      L=0
      DO  260 K=1,3
      J1=IBEG(K)
      J2=IEND(K)
      J3=IN(K)
      DO 260 I=J1,J2,J3
      L=L+1
      J=I-1
      PHI=ASIN(S(I))/RAD
      FANG(L)=PHI
      PRESS=PINFT-RHOWG*DELP(I)
      R0=J
      XX=VEE(R0)
      D(L)=J
 260  CONTINUE
      DO  104  KK=75,81
      D(KK+1)=D(KK)+50.0
      FANG(KK+1)=FANG(75)
104   CONTINUE
 999  RETURN
      END
      SUBROUTINE INTWND(PDROP,C24,C23,ALAT)
C      JELESNIANSKI AUG 1980TDL IBM 360/195
C      PURPOSE
C      THIS SUBROUTINE (INTWND) SETS UP AN INITIAL GUESS FOR THE
C      MAX WIND OF A STATIONARY STORM USING INITIAL INPUT STORM
C      PARAMETERS.THE INITIAL GUESS IS AN ENTRY VALUE TO
C      SUBROUTINE 'FORCE1'FOR ITERATION TO EXACT MAXIMUM WINDS.
C      DATA SET USE
C      NONE
C      VARIABLES
C      PDROP=PRESSURE DROP FOR 1ST HOUR
C      ALAT=LATITUDE OF BASIN CENTER
C      KPRES=DOCTORED PRESSURE DROP FOR LAT OTHER THAN 30 DEG
C      WM(150)=ARRAY FOR SLOPE OF LINEAR SLOPE FORMULA
C      WB(150)=ARRAY FOR INTERCEPT OF LINEAR SLOPE FORMULA
C      C24=RADIUS OF MAXIMUM WIND
C      C23=MAXIMUM WIND
C      GENERAL COMMENTS
C      THIS SUBROUTINE LIES IN OVERLAY STRMPR.IT IS CALLED IN BY
C      SUBROUTINE 'STRMPR'
C      FOR A GIVEN PRESSURE DROP 'P(100)',THE MAX WIND
C      VARIES ALMOST LINEARLY WITH STORM SIZE 'R(100)'.
C      SEE BELOW:
C
C             -       E
C             -      R
C             -     U
C             W    S  WM(KPRES)=SLOPE OF PRESSURE  LINE
C             I   S
C             N  E
C             D R
C             -P
C      WB(KPRES)=INTERCEPT-
C                         0---R=STORM SIZE----
C      THE SLOPE AND INTERCEPT OF THE PRESSURE LINE VARIES
C      WEAKLY WITH LATITUDE OR CORIOLIS,THE GIVEN WM,WB,ARE
C      FOR LATITUDE 30 DEG:THE ENTRY PRESSURE AT TIME 'IBGNT'
C      IS DOCTORED TO CORRECT FOR LATITUDE OTHER THAN 30 DEG.
      DIMENSION     WM(150),WB(150)
      DATA WM/.10,.14,.17,.20,.23,.25,.28,.30,.32,.33,.35,.37,.38,.39,
     1.41,.42,.43,.44,.45,.46,.46,.47,.48,.49,.49,.50,.51,.52,.52,.53,
     2.53,.54,.54,.55,.56,.56,.57,.57,.58,.58,.59,.59,.60,.60,.61,.62,
     3.62,.63,.63,.64,.65,.65,.66,.66,.67,.68,.69,.69,.70,.71,.71,.72,
     4.72,.73,.73,.74,.74,.75,.75,.76,.76,.76,.77,.77,.77,.77,.78,.78,
     5.78,.78,.79,.79,.79,.79,.80,.80,.80,.80,.80,.81,.81,.81,.81,.81,
     6.82,.82,.82,.82,.82,.82,.83,.83,.83,.83,.84,.84,.84,.84,.84,.84,
     7.84,.85,.85,.85,.85,.85,.85,.85,.86,.86,.86,.86,.86,.86,.86,.86,
     8.86,.86,.87,.87,.87,.87,.87,.87,.87,.87,.87,.87,.87,.87,.87,.87,
     9.87,.87,.87,.87,.87,.87,.87,.87/
      DATA WB/12.,23.,29.,33.,37.,39.,41.,43.,45.,47.,49.,51.,53.,55.,
     157.,58.,60.,62.,63.,65.,67.,68.,70.,71.,73.,74.,76.,77.,78.,80.,
     281.,82.,84.,85.,86.,88.,89.,90.,92.,93.,94.,95.,96.,98.,99.,100.,
     3101.,102.,104.,105.,106.,107.,108.,109.,110.,112.,113.,114.,115.,
     4116.,117.,118.,119.,120.,121.,122.,123.,123.,124.,125.,126.,127.,
     5128.,129.,130.,131.,132.,132.,133.,134.,135.,136.,137.,137.,138.,
     6139.,140.,141.,141.,142.,143.,144.,145.,145.,146.,147.,148.,148.,
     7149.,150.,151.,151.,152.,153.,154.,154.,155.,156.,157.,157.,158.,
     8158.,159.,160.,160.,161.,162.,162.,163.,164.,164.,165.,166.,166.,
     9167.,167.,168.,169.,169.,170.,171.,171.,172.,172.,173.,174.,174.,
     *175.,175.,176.,176.,177.,177.,178.,179.,179.,180.,180.,181.,182./ 254
C     DOCTOR ENTRY PRE FOR LATITUDE OTHER THAN 30 DEG
      PR=PDROP
      PF=((PR/67.7)**2)*.3*(30.-ALAT)
      KPRES=PR+PF+0.5
C     LINEAR SLOPE FORMULA FOR APPROXIMATE MAXIMUM WIND
      C23=-WM(KPRES)*C24+WB(KPRES)
      RETURN
      END
      SUBROUTINE STRES(D,FANG,MAP,TX,TY,XC,YC,USS,BETA,M,N)
      CHARACTER *30NAME
      DIMENSION D(200),FANG(100),MAP(20,49),TX(20,49),TY(20,49)
      REAL MAP
      COMMON /A/AL
      COMMON /C/DP,RMAX,VMAX,RHOAIR,NAME
      BETA1=BETA*0.0174
      USX=USS*COS(BETA1)
      USY=USS*SIN(BETA1)
      RS2=RMAX*RMAX
      VR2=2.0*VMAX*0.447
      DY=AL/(N-1)
      CD=2.8E-3
      L=82
      DO 100 J=1,N
      DO 100 I=1,M
      X=MAP(I,J)
      RR=SQRT((X-XC)**2+(DY*FLOAT(J-1)-YC)**2)
      RR=RR*0.6214/1000.
      RSQR=RR*RR
      CC=1.0/(RS2+RSQR)
      L=L-1
      DO 101 K=1,L
      IF(D(K).LT.RR.AND.D(K+1).GT.RR)GO TO 102
      GO TO 101
102   DIFF1=ABS(D(K)-RR)
      DIFF2=ABS(D(K+1)-RR)
      FANGLE=(DIFF1*FANG(K+1)+DIFF2*FANG(K))/(DIFF1+DIFF2)
      FANGLE=FANGLE*0.0174
      GO TO 103
101   CONTINUE
103   CONTINUE
      STHETA=(X-XC)*0.6214/(RR*1000.)
      CTHETA=(DY*FLOAT(J-1)-YC)*0.6214/(RR*1000.)
      SFAN=SIN(FANGLE)
      CFAN=COS(FANGLE)
      A=USX/VR2-(STHETA*CFAN+CTHETA*SFAN)
      B=USY/VR2+(CTHETA*CFAN-STHETA*SFAN)
      FCT=VR2*RMAX*RR/(RS2+RR*RR)
      FCT=FCT*FCT
      FCT1=CD*RHOAIR*FCT*SQRT(A*A+B*B)
      TX(I,J)=FCT1*A
      TY(I,J)=FCT1*B
100   CONTINUE
      RETURN
      END
      SUBROUTINE PLOT1(N,PEAK,IC)
      INTEGER IB(402)
      CHARACTER BB,P,Q,R
      CHARACTER *16STN
      CHARACTER *30NAME
      CHARACTER *4 AS(402)
      DIMENSION PEAK(402),STN(402),BB(402)
      COMMON /C/DP,RR,VMAX,RHOAIR,NAME
      DATA P,Q,R/' ','+','*'/
      VMAX=VMAX*0.447
      RR=RR/0.6214
      DO 800 I1=1,N
800   IB(I1)=0
      DO 900 I1=1,N
900   STN(I1)='                 '
      STN(107)='VEDARANNIYAM'
      STN(117)='NAGAPPATTINAM'
      STN(135)='CUDDALORE'
      STN(143)='PONDICHERRY'
      STN(160)='TBM'
      STN(167)='MADRAS'
      STN(199)='NELLORE'
      STN(207)='KAVALI'
      STN(225)='ONGOLE'
      STN(243)='DIVI'
      STN(247)='MACHILIPATNAM'
      STN(270)='KAKINADA'
      STN(289)='VISHAKHAPATNAM'
      STN(309)='KALINGAPATNAM'
      STN(328)='GOPALPUR'
      STN(345)='PURI'
      STN(362)='PARADEEP'
      STN(376)='CHANDABALI'
      STN(394)='BALASORE'
      STN(402)='CONTAI'
      STN(IC)='     LANDFALL-->'
      IB(IC)=0
      IN1=IC+5
      IN2=IC-5
      N1=(N-IC)/5
      N2=IC/5
      DO 60 I1=1,N1
      IB(IN1)=I1*25
      IN1=IN1+5
60    CONTINUE
      DO 70 I2=1,N2
      IB(IN2)=I2*25
      IN2=IN2-5
70    CONTINUE
C     ENCODE(620,15,AS)IB
C15    FORMAT(155I4)
      DO 50 I=1,N
      IF(I.EQ.IC)GO TO 50
      IF(IB(I).EQ.0) AS(I)='    '
50    CONTINUE
      BMAX=0.0
      DO 1 J=1,N
      IF(PEAK(J).GT.BMAX)BMAX=PEAK(J)
1     CONTINUE
      WRITE(6,52)NAME,DP,VMAX,RR
52    FORMAT(7X,'CYCLONE:',A30,'PRDROP(MB)=',F5.1,5X,'VMAX(M/S)=',
     *F5.1,5X,'RMAX(KM)=',F5.1//)
      WRITE(6,51)
51    FORMAT(7X,'STATION',9X,'DIST.',1X,'SURGE',10X,
     *'PEAK SURGE IN MTS------>'/,23X,'(KM)',2X,'(MTS)'/)
      IMIN=1
      IMAX=N
      IF(IMIN.LT.1)IMIN=1
      IF(IMAX.GT.N)IMAX=N
      DO 3 JJ=IMIN,IMAX
      J=IMAX+IMIN-JJ
      MAX=70
      M2=50-INT((BMAX-PEAK(J))/BMAX*49.0)
      DO 4 I=1,MAX
      BB(I)=P
      IF(I.EQ.M2)BB(I)=R
4     CONTINUE
      WRITE(6,101)J,STN(J),AS(J), PEAK(J),(BB(I),I=1,MAX)
3     CONTINUE
101   FORMAT(2X,I3,2X,A16,A4,F6.2,5X,70A1)
      RETURN
      END
        SUBROUTINE PLOT(PEAK,XDIST,IC,NAME) 
C       PROGRAM TO CREATE A PLOT FILE TO DRAW INDIA MAP
        DIMENSION PEAK(402),XDIST(402),XW(5),YW(5)
        CHARACTER *50INFILE,NAME,PLOTITLE
        DATA XW/2.0,12.0,12.0,2.0,2.0/
        DATA YW/1.0,1.0,8.0,8.0,1.0/
C        CALL CS(23)
C        WRITE(*,*)'ENTER NAME OF FILE FOR PLOTING PEAKSURGE' 
C        WRITE(*,*)'FILE NAME SHOULD BE WITH EXTENSION OF .PLT:>'
C        READ(*,199)INFILE
C        CALL CS(23)
C        WRITE(*,*)'ENTER NAME OF THE PLOT:>'
C        READ(*,199)PLOTITLE
199     FORMAT(A50)        
        OPEN(3,FILE='PS.PLT')
        X1=1.5
        Y1=4.0
        XL=11.0
        YL=6.0
        n=402
        XMIN=0.0
        XMAX=2010.0
        YMIN=0.0
        YMAX=8.0
        XS=XL/(XMAX-XMIN)
        YS=YL/(YMAX-YMIN)
        DO 5 I=1,N
        XDIST(I)=XS*XDIST(I)-XMIN*XS+X1
        PEAK(I)=YS*PEAK(I)-YMIN*YS+Y1
5       CONTINUE
C       DRAW THE WINDOW
        WRITE(3,96)
96      FORMAT('PI 0.000 0.000')
        WRITE(3,97)
97      FORMAT('TR 0.000 0.000')
        WRITE(3,98)
98      FORMAT('SC 1.000 1.000')
        WRITE(3,99)
99      FORMAT('SP1')
        WRITE(3,966)
966     FORMAT('PI 1.500 4.000')
        WRITE(3,967)
967     FORMAT('RO 0.000')
        WRITE(3,968)
968     FORMAT('MA 1.500 4.000')
        WRITE(3,969)
969     FORMAT('PA 12.500 4.000')
        WRITE(3,970)
970     FORMAT('RO 0.000')
        WRITE(3,971)
971     FORMAT('PI 1.500 4.000')
        WRITE(3,972)
972     FORMAT('RO 90.000 ')
        WRITE(3,973)
973     FORMAT('MA 1.500 4.000')
        WRITE(3,974)
974     FORMAT('PA 7.500 4.000')
        WRITE(3,981)
981     FORMAT('RO -90.000')
        WRITE(3,975)
975     FORMAT('MA 1.500 4.000')
        WRITE(3,976)
976     FORMAT('PA 12.500 4.000')
        WRITE(3,977)
977     FORMAT('MA 12.500 4.000')
        WRITE(3,978)
978     FORMAT('PA 12.500 10.000')
        WRITE(3,979)
979     FORMAT('MA 12.500 10.000')
        WRITE(3,980)
980     FORMAT('PA 1.500 10.000')
C       PLOTTING OF TICS AND LABELS ON X AXIS   
        PY=10.000
        P1Y=10.113
        P2Y=10.123
        P1X=0.0*XS+X1
        P2X=250.0*XS+X1
        P3X=500.0*XS+X1
        P4X=750.0*XS+X1
        P5X=1000.0*XS+X1
        P6X=1250.0*XS+X1
        P7X=1500.0*XS+X1
        P8X=1750.0*XS+X1
        P9X=2010.0*XS+X1
        CHT=0.1
        ANG=90.0
        WRITE(3,221)
221     FORMAT('PI 1.500 10.000')
        WRITE(3,222)
222     FORMAT('RO 0.000 ')
        WRITE(3,223)
223     FORMAT('SS "DEFAULT.SYM"')
        WRITE(3,224)P1X,PY
224     FORMAT('MA',1X,2F7.3)
        WRITE(3,225)P1X,P1Y
225     FORMAT('PA',1X,2F7.3)
        WRITE(3,226)P1X+0.05,P2Y,CHT,ANG
226     FORMAT('PS',1X,4F7.3,1X,'"0"')
        WRITE(3,224)P2X,PY
        WRITE(3,225)P2X,P1Y
        WRITE(3,227)P2X+0.05,P2Y,CHT,ANG
227     FORMAT('PS',1X,4F7.3,1X,'"250"')
        WRITE(3,224)P3X,PY
        WRITE(3,225)P3X,P1Y
        WRITE(3,228)P3X+0.05,P2Y,CHT,ANG
228     FORMAT('PS',1X,4F7.3,1X,'"500"')
        WRITE(3,224)P4X,PY
        WRITE(3,225)P4X,P1Y
        WRITE(3,229)P4X+0.05,P2Y,CHT,ANG
229     FORMAT('PS',1X,4F7.3,1X,'"750"')
        WRITE(3,224)P5X,PY
        WRITE(3,225)P5X,P1Y
        WRITE(3,230)P5X+0.05,P2Y,CHT,ANG
230     FORMAT('PS',1X,4F7.3,1X,'"1000"')
        WRITE(3,224)P6X,PY
        WRITE(3,225)P6X,P1Y
        WRITE(3,231)P6X+0.05,P2Y,CHT,ANG
231     FORMAT('PS',1X,4F7.3,1X,'"1250"')
        WRITE(3,224)P7X,PY
        WRITE(3,225)P7X,P1Y
        WRITE(3,232)P7X+0.05,P2Y,CHT,ANG
232     FORMAT('PS',1X,4F7.3,1X,'"1500"')
        WRITE(3,224)P8X,PY
        WRITE(3,225)P8X,P1Y
        WRITE(3,233)P8X+0.05,P2Y,CHT,ANG
233     FORMAT('PS',1X,4F7.3,1X,'"1750"')
        WRITE(3,224)P9X,PY
        WRITE(3,225)P9X,P1Y
        WRITE(3,234)P9X+0.05,P2Y,CHT,ANG
234     FORMAT('PS',1X,4F7.3,1X,'"2010"')
C       DRAW THE MAP
        WRITE(3,996)
996     FORMAT('SP 2')
        WRITE(3,966)
        WRITE(3,967)
        WRITE(3,102)XDIST(1),PEAK(1)
        DO 7 I=2,N
7       WRITE(3,103)XDIST(I),PEAK(I)
102     FORMAT('MA',2F6.3)
103     FORMAT('PA',2F6.3)
C       PLOTTING OF TICS AND LABELS ON Y AXIS
        PX=1.5
        P11X=4.113
        P22X=4.153
        P111Y=0.5*YS+X1
        P22Y=1.0*YS+X1
        P3Y=1.5*YS+X1
        P4Y=2.0*YS+X1
        P5Y=2.50*YS+X1
        P6Y=3.0*YS+X1
        P7Y=3.50*YS+X1
        P8Y=4.00*YS+X1
        P9Y=4.50*YS+X1
        P10Y=5.0*YS+X1
        P11Y=5.5*YS+X1
        P12Y=6.0*YS+X1
        p13y=6.5*ys+x1
        p14y=7.0*ys+x1
        p15y=7.5*ys+x1
        p16y=8.0*ys+x1
        ANG1=0.0
        WRITE(3,2214)
2214    FORMAT('SP4')
        WRITE(3,2211)
2211    FORMAT('PI 1.500 4.000')
        WRITE(3,1222)
1222    FORMAT('RO 90.000 ')
        WRITE(3,223)
        WRITE(3,224)PX,Y1
        WRITE(3,225)PX,P11X
        WRITE(3,226)PX-0.2,P22X,CHT,ANG1
2226    FORMAT('PS',1X,4F7.3,1X,'"0.0"')
        WRITE(3,224)P111Y,Y1
        WRITE(3,225)P111Y,P11X
        WRITE(3,1227)P111Y-0.2,P22X,CHT,ANG1
1227    FORMAT('PS',1X,4F7.3,1X,'"0.5"')
        WRITE(3,224)P22Y,Y1
        WRITE(3,225)P22Y,P11X
        WRITE(3,1228)P22Y-0.2,P22X,CHT,ANG1
1228    FORMAT('PS',1X,4F7.3,1X,'"1.0"')
        WRITE(3,224)P3Y,Y1
        WRITE(3,225)P3Y,P11X
        WRITE(3,1229)P3Y-0.2,P22X,CHT,ANG1
1229    FORMAT('PS',1X,4F7.3,1X,'"1.5"')
        WRITE(3,224)P4Y,Y1
        WRITE(3,225)P4Y,P11X
        WRITE(3,1230)P4Y-0.2,P22X,CHT,ANG1
1230    FORMAT('PS',1X,4F7.3,1X,'"2.0"')
        WRITE(3,224)P5Y,Y1
        WRITE(3,225)P5Y,P11X
        WRITE(3,1231)P5Y-0.2,P22X,CHT,ANG1
1231    FORMAT('PS',1X,4F7.3,1X,'"2.5"')
        WRITE(3,224)P6Y,Y1
        WRITE(3,225)P6Y,P11X
        WRITE(3,1232)P6Y-0.2,P22X,CHT,ANG1
1232    FORMAT('PS',1X,4F7.3,1X,'"3.0"')
        WRITE(3,224)P7Y,Y1
        WRITE(3,225)P7Y,P11X
        WRITE(3,1233)P7Y-0.2,P22X,CHT,ANG1
1233    FORMAT('PS',1X,4F7.3,1X,'"3.5"')
        WRITE(3,224)P8Y,Y1
        WRITE(3,225)P8Y,P11X
        WRITE(3,1234)P8Y-0.2,P22X,CHT,ANG1
1234    FORMAT('PS',1X,4F7.3,1X,'"4.0"')
        WRITE(3,224)P9Y,Y1
        WRITE(3,225)P9Y,P11X
        WRITE(3,1235)P9Y-0.2,P22X,CHT,ANG1
1235    FORMAT('PS',1X,4F7.3,1X,'"4.5"')
        WRITE(3,224)P10Y,Y1
        WRITE(3,225)P10Y,P11X
        WRITE(3,1236)P10Y-0.2,P22X,CHT,ANG1
1236    FORMAT('PS',1X,4F7.3,1X,'"5.0"')
        WRITE(3,224)P11Y,Y1
        WRITE(3,225)P11Y,P11X
        WRITE(3,1237)P11Y-0.2,P22X,CHT,ANG1
1237    FORMAT('PS',1X,4F7.3,1X,'"5.5"')
        WRITE(3,224)P12Y,Y1
        WRITE(3,225)P12Y,P11X
        WRITE(3,1238)P12Y-0.2,P22X,CHT,ANG1
1238    FORMAT('PS',1X,4F7.3,1X,'"6.0"')
        WRITE(3,224)P13Y,Y1
        WRITE(3,225)P13Y,P11X
        WRITE(3,1239)P13Y-0.2,P22X,CHT,ANG1
1239    FORMAT('PS',1X,4F7.3,1X,'"6.5"')
        WRITE(3,224)P14Y,Y1
        WRITE(3,225)P14Y,P11X
        WRITE(3,1240)P14Y-0.2,P22X,CHT,ANG1
1240    FORMAT('PS',1X,4F7.3,1X,'"7.0"')
        WRITE(3,224)P15Y,Y1
        WRITE(3,225)P15Y,P11X
        WRITE(3,1241)P15Y-0.2,P22X,CHT,ANG1
1241    FORMAT('PS',1X,4F7.3,1X,'"7.5"')
        WRITE(3,224)P16Y,Y1
        WRITE(3,225)P16Y,P11X
        WRITE(3,1242)P16Y-0.2,P22X,CHT,ANG1
1242    FORMAT('PS',1X,4F7.3,1X,'"8.0"')


        WRITE(3,3344)
3344    FORMAT('SP 5')
        WRITE(3,223)
        WRITE(3,3341)
3341    FORMAT('PS 3.863 4.353 0.100 0.000 "SURGE IN METERS"')
C       PLOTTING THE NAMES OF THE STATIONS
        TY=2.606
        TY1=3.1
        TY2=9.80
        TX1=107*5*XS+X1
        TX2=117*5*XS+X1
        TX3=135*5*XS+X1
        TX4=143*5*XS+X1
        TX5=167*5*XS+X1
        TX6=199*5*XS+X1
        TX7=207*5*XS+X1
        TX8=225*5*XS+X1
        TX9=243*5*XS+X1
        TX10=247*5*XS+X1
        TX11=270*5*XS+X1
        TX12=289*5*XS+X1
        TX13=309*5*XS+X1
        TX14=328*5*XS+X1
        TX15=345*5*XS+X1
        TX16=362*5*XS+X1
        TX17=376*5*XS+X1
        TX18=394*5*XS+X1
        TX19=402*5*XS+X1
        TX20=4.8
        TXLP=IC*5*XS+X1
        TX21=2.620
        TY21=9.148
        TX22=2.420
        TY22=8.748
        TX23=2.720
        TY23=1.740
        CHT2=0.220
        CHT3=0.150        
        WRITE(3,3310)
3310    FORMAT('SP6')
        WRITE(3,2211)
        WRITE(3,3311)
3311    FORMAT('RO -90.000')
        WRITE(3,223)
        WRITE(3,3312)TX1,TY,CHT,ANG
3312    FORMAT('PS',1X,4F7.3,1X,'"VEDANARAYANYAM"')
        WRITE(3,3313)TX2,TY,CHT,ANG
3313    FORMAT('PS',1X,4F7.3,1X,'"NAGAPATTINAM"')
        WRITE(3,3314)TX3,TY,CHT,ANG
3314    FORMAT('PS',1X,4F7.3,1X,'"CUDDALORE"')
        WRITE(3,3315)TX4,TY,CHT,ANG
3315    FORMAT('PS',1X,4F7.3,1X,'"PONDYCHERRY"')
        WRITE(3,3316)TX5,TY,CHT,ANG
3316    FORMAT('PS',1X,4F7.3,1X,'"MADRAS"')
        WRITE(3,3317)TX6,TY,CHT,ANG
3317    FORMAT('PS',1X,4F7.3,1X,'"NELLORE"')
        WRITE(3,3318)TX7,TY,CHT,ANG
3318    FORMAT('PS',1X,4F7.3,1X,'"KAVALI"')
        WRITE(3,3319)TX8,TY,CHT,ANG
3319    FORMAT('PS',1X,4F7.3,1X,'"ONGOLE"')
        WRITE(3,3320)TX9,TY,CHT,ANG
3320    FORMAT('PS',1X,4F7.3,1X,'"DIVI"')
        WRITE(3,3321)TX10,TY,CHT,ANG
3321    FORMAT('PS',1X,4F7.3,1X,'"MACHILIPATNAM"')
        WRITE(3,3322)TX11,TY,CHT,ANG
3322    FORMAT('PS',1X,4F7.3,1X,'"KAKINADA"')
        WRITE(3,3323)TX12,TY,CHT,ANG
3323    FORMAT('PS',1X,4F7.3,1X,'"VISAKHAPATNAM"')
        WRITE(3,3324)TX13,TY,CHT,ANG
3324    FORMAT('PS',1X,4F7.3,1X,'"KALINGAPATNAM"')
        WRITE(3,3325)TX14,TY,CHT,ANG
3325    FORMAT('PS',1X,4F7.3,1X,'"GOPALPUR"')
        WRITE(3,3326)TX15,TY,CHT,ANG
3326    FORMAT('PS',1X,4F7.3,1X,'"PURI"')
        WRITE(3,3327)TX16,TY,CHT,ANG
3327    FORMAT('PS',1X,4F7.3,1X,'"PARADIP"')
        WRITE(3,3328)TX17,TY,CHT,ANG
3328    FORMAT('PS',1X,4F7.3,1X,'"CHANDBALI"')
        WRITE(3,3329)TX18,TY,CHT,ANG
3329    FORMAT('PS',1X,4F7.3,1X,'"BALASORE"')
        WRITE(3,3330)TX19,TY,CHT,ANG
3330    FORMAT('PS',1X,4F7.3,1X,'"CONTAI"')
        WRITE(3,3332)TX20,TY2,CHT,ANG1
3332    FORMAT('PS',1X,4F7.3,1X,29H"DISTANCE ALONG EAST COAST OF, 
     #  1X,12HINDIA IN KM")
        WRITE(3,3333)
3333    FORMAT('SP 5')
        WRITE(3,223)
        WRITE(3,3331)TXLP,TY1,CHT,ANG
3331    FORMAT('PS',1X,4F7.3,1X,'"LANDFALL -->"')
        WRITE(3,3390)
3390    FORMAT('SP 6')
        WRITE(3,223)
        WRITE(3,3391)TX21,TY21,CHT2,ANG1
3391    FORMAT('PS',1X,4F7.3,1X,'"IIT-D SURGE MODEL"')
        WRITE(3,3392)
3392    FORMAT('SP 4')
        WRITE(3,223)
        WRITE(3,3393)TX22,TY22,CHT3,ANG1,NAME
3393    FORMAT('PS',1X,4F7.3,1X,'"',A50,'"')
        WRITE(3,3392)
3433    FORMAT('SP 4')
        WRITE(3,223)
        WRITE(3,3434)TX23,TY23,CHT3,ANG1
3434    FORMAT('PS',1X,4F7.3,1X,'"MAXIMUM PEAK SURGE ENVELOPE"')
        RETURN
        END
        SUBROUTINE INPUT(KS,DP,NAME)
C       THIS SUBROUTINE GIVES INPUT DATA FROM THE SCREEN
        CHARACTER NAME*50
        call cs(8)
      write(*,*)'               WELCOME TO IIT DELHI STORM SURGE MODEL'                        
      WRITE(*,*)'                   CENTER FOR ATMOSPHERIC SCIENCES'
      WRITE(*,*)'              INDIAN INSTITUTE OF TECNOLOGY, NEW DELHI'
        call cs(10)
        WRITE(*,*)'ENTER NAME OF THE CYCLONE:>'
        READ (*,199)NAME
199     FORMAT(A50)
C       call cs(23)
C       WRITE(*,*)'ENTER MAXIMUM RADIUS OF WINDS IN METERS:>'
C       call cs(23)
C       READ(*,*)R
        call cs(23) 
        WRITE(*,*)'ENTER PRESSURE DROP IN MILLI BARS:>'
        READ(*,*)DP
        call cs(23)
        WRITE(*,*)'ENTER TOTAL NUMBER OF ITERATIONS:>'
        call cs(23)
        READ(*,*)KS
C       call cs(23)
C       WRITE(*,*)'ENTER TOTAL NUMBER OF STORM TRACK POSITIONS:>'
C       READ(*,*)ND
        call cs(23)
C       WRITE(*,*)'ENTER LAND FALL POINT:>'
C       READ(*,*)IC
        END
      SUBROUTINE CS(N)
      INTEGER I1
      INTEGER N
      DO 10 I1=1,N
         WRITE(*,*)
   10 CONTINUE
      RETURN  
      END
      SUBROUTINE CONTOUR(B1,B,XI,AL,PZ1)
C     *******************************************************
C      CONVERSION TO LATITUDE AND LONGITUDE      
        DIMENSION B1(49),B(49),XI(42),PZ1(20,49),PZ3(40,49)
        DIMENSION XLAT(40,49),YLAT(20,49),YLONG(20,49),XLONG(40,49),
     1  BLONG(20,49),xxlong(20,49) 
        CHARACTER *30INFILE
        CALL CS(23)
        WRITE(*,*)'ENTER NAME OF FILE FOR PLOTING PEAK SURGE CONTOURS' 
        WRITE(*,*)'FILE NAME SHOULD BE WITH EXTENSION OF .DAT:>'
        READ(*,199)INFILE
199     FORMAT(A30)        
        OPEN(7,FILE=INFILE)     
      YLLL=1./(111.*1.0E3)
      DO 35 J=1,49   
      DO 35 I=1,20
      YLONG(I,J)=(81.84+B1(J)*YLLL)+YLLL*B(J)*XI(I)
c      ylong(i,j)=ylong(i,j)+0.3
35    CONTINUE      
      EEE=(250.*1.0e3)/FLOAT(19)
      EEE=EEE*YLL
      DO 315 J=1,49   
      DO 315 I=1,20
      XXlong(i,j)=(81.84+b1(j)*ylll)+ylll*b(j)*EEE
315    CONTINUE      
      DDD=(AL/FLOAT(48))*YLLL
      DO 36 I=1,20
36    YLAT(I,1)=6.0
      DO 38 I=1,20
      DO 38 J=2,49
38    YLAT(I,J)=YLAT(I,J-1)+DDD      
        DO 2434 J=1,49
2434    BLONG(1,J)=YLONG(1,J)-0.2
        DO 2435 J=1,49   
        DO 2435 I=2,20
2435    BLONG(I,J)=BLONG(I-1,J)-0.3
        DO 2436 J=1,49
        DO 2436 I=1,20
        XLAT(I,J)=YLAT(I,J)
        II=20+1-I
2436    XLONG(II,J)=BLONG(I,J)
        DO 2437 J=1,49
        DO 2437 I=1,20
        II=I+20 
        XLAT(II,J)=YLAT(I,J)
2437    XLONG(II,J)=YLONG(I,J)        
        DO 2441 J=1,49
        DO 2441 I=2,20
        IK=20+2-I
2441    PZ3(IK,J)=PZ1(I,J)    
        DO 2439 J=1,49
        DO 2439 I=1,20
        II=I+20
C        PZ3(IK,J)=PZ1(1,J)
2439    PZ3(II,J)=PZ1(I,J)
C        DO 2440 J=1,49
C        DO 2440 I=1,20
C        WRITE(18,*)YLONG(I,J),YLAT(I,J),PZ1(I,J)        
C2440    WRITE(21,*)XXLONG(I,J),YLAT(I,J),PZ1(I,J)        
        DO 2438 J=1,49
        DO 2438 I=1,40
C        XLAT(I,J)=XLAT(I,J)-0.5
2438    WRITE(7,*)XLONG(I,J),XLAT(I,J),PZ3(I,J)        
        RETURN
        END