C$$$  MAIN PROGRAM DOCUMENTATION BLOCK
C
C MAIN PROGRAM: EXT_6H
C   PRGMMR: CHEN             ORG: OST25       DATE: 2003-03-17
C
C ABSTRACT: 2-DIMENSIONAL MODEL FOR EXTRATROPICAL STORM SURGES
C           MODIFICATION: 1. CALIBRATED DRAG COEFFICIENT BASED DATA
C                         2. 6 HR'LY ANALYSIS FROM 4-CYCLE AVN RUNS
C
C PROGRAM HISTORY LOG:
C   94-10-25  CHEN
C   94-10-25  KIM C90 IMPLEMENTATION
C   95-12-22  MO PUT OLD FLE10GEN.F INTO THE END OF THE MAIN PROGRAM
C   01-01-01  CHEN MODIFY DRAG COEFFICIENT AND USE 6 HR'LY INCREMENT FOR
C             SPIN-UP PERIOD
C             GENERATE 10-ANALYSIS FILE EVERY CYCLE (2.5 DAYS), DO 
C             WITHOUT UP-KEEPING ROTATING FILE.
C
C USAGE:
C   FOR BAINS 'X' AND CYCLE T'HH'Z
C
C   INPUT FILES:
C     FORT.19  - ETBSN.X  EAST COAST BASIN DATA (FORMER FORT.9)
C     FORT.14  - ETTGP.X  TIDE GAGE STATIONS TO BE MONITORED (FORT.4)
C     FORT.11  - FT11.X   CONTROL WINDOW INFORMATION
C     FORT.21  - FT01.EGA   MODEL RUN TIME CONTROL FILE (FORMER FORT.1)
C     FORT.25  - ETLTLG.X LATITUDE-LONGITUDE OF MOMENTUM POINTS
C     FORT.33  - CYLF10.X  10 CYCLES AVN ANALYSES
C     FORT.34  - AVNPUV.HHX  AVN FORECAST
C     FORT.96  - SDS.HH LOG FILE OF AVN_PUV PROGRAM
C 
C   TEMPORARY FILES:
C     FORT.81  - SURGE HEIGHT   (FORMER FLE10)
C     FORT.82  - WIND SPEED     (FORMER FLE20)
C     FORT.83  - WIND DIRECTION (FORMER FLE30)
C
C   OUTPUTFILES:
C     FORT.52   - SSHISTORY.HHX    SURGE ARCHIVE
C
C   SUBROUTINES ASSOCIATED WITH THE ROOT OVERLAY ARE
C         1.  BLCKDT
C         2.  MAIN PROGRAM
C
C ATTRIBUTES:
C   LANGUAGE: FORTRAN77
C   MACHINE:  C90
C$$$
      COMMON /CESAV/  IPN(120),IPL(120),KHSPT
      COMMON /CESAV1/ STATS(120)
      COMMON /IFT11 / LATIN1,LATIN2,LONIN1,LONIN2
      COMMON /INZSTM/ IBGS,ITMS
      COMMON /GEO/  NOGEO
      COMMON /OPT/  NOPT
      COMMON /INTV/ INTV(24),KINPUT
      COMMON /KEYS/ KEY1
      COMMON /FLES/ FLE5,FLE9,FLE8,FLE91,FLE99,FLE10,FLE20,FLE30,FLE1
      COMMON /FLE/  FLE92
      CHARACTER*30  FLE5,FLE9,FLE8,FLE91,FLE99,FLE10,FLE20,FLE30,FLE1
      CHARACTER*30  FLE92
      CHARACTER*1   AAA,ACHR
      CHARACTER*30 STATS 
      CHARACTER*72 FIL_96,FIL_11,FIL_14,FIL_19,FIL_21,FIL_25,FIL_33,
     1 FIL_34,FIL_81,FIL_82,FIL_83,FIL_52

      CALL W3TAGB('EXT_6H',2007,0130,0000,'OST25')
C
C      CALL GETENV('FLE_96',VALUE=FIL_96)
C      OPEN (96,FILE=FIL_96)
C      CALL GETENV('FLE_11',VALUE=FIL_11)
C      OPEN (11,FILE=FIL_11)
C      CALL GETENV('FLE_14',VALUE=FIL_14)
C      OPEN (14,FILE=FIL_14)
C      CALL GETENV('FLE_19',VALUE=FIL_19)
C      OPEN (19,FILE=FIL_19)
C      CALL GETENV('FLE_21',VALUE=FIL_21)
C      OPEN (21,FILE=FIL_21)
C      CALL GETENV('FLE_25',VALUE=FIL_25)
C      OPEN (25,FILE=FIL_25,FORM='UNFORMATTED')
C      CALL GETENV('FLE_33',VALUE=FIL_33)
C      OPEN (33,FILE=FIL_33,FORM='UNFORMATTED')
C      CALL GETENV('FLE_34',VALUE=FIL_34)
C      OPEN (34,FILE=FIL_34,FORM='UNFORMATTED')
C      CALL GETENV('FLE_81',VALUE=FIL_81)
C      OPEN (81,FILE=FIL_81,FORM='UNFORMATTED')
C      CALL GETENV('FLE_82',VALUE=FIL_82)
C      OPEN (82,FILE=FIL_82,FORM='UNFORMATTED')
C      CALL GETENV('FLE_83',VALUE=FIL_83)
C      OPEN (83,FILE=FIL_83,FORM='UNFORMATTED')
C      CALL GETENV('FLE_52',VALUE=FIL_52)
C      OPEN (52,FILE=FIL_52,FORM='UNFORMATTED')
C
      READ (96,'(A1)') AAA
      IF (AAA.EQ.'N') THEN
      CALL W3TAGE('EXT_6H')
      STOP
      ENDIF

       READ (11,*) LATIN1,LATIN2
       READ (11,'(A1)') ACHR
       READ (11,*) LONIN1,LONIN2
       READ (11,*)

C ---------------- GAGE STATIONS USED FOR OUTPUT ------------
      READ (14,'(I3)') KHSPT
      DO 2001 N=1,KHSPT
      READ (14,'(A30,18X,2I4)') STATS(N),IPN(N),IPL(N)
2001  CONTINUE
      KEY1=1

      CALL INITLZ
      CALL CMPUTE
      REWIND (81)
      REWIND (82)
      REWIND (83)
      CALL FLEGEN

C     CALL PLOTDP
      CALL W3TAGE('EXT_6H')
      STOP
      END
C
         BLOCK DATA
C        J. CHEN     JULY 1988    MDL   FORTRAN 77 VERSION
C
C        PURPOSE
C        DATA INITIALIZATION
C           NONE
C
C        VARIABLES
C       IP(4) JP(4) = SHIFT SUBSCRIPTS T 4 ADJACENT MOMENTUM POINTS
C       IH(4) JH(4) = SHIFT SUBSCRIPTS TO 4 ADJACENT SURGE POINTS ABOUT
C                     A MOMUMTUM POINT
C     IIH(4) JJH(4) = SHIFT SUBSCRIPTS TO 4 ADJACENT SURGE POINTS ABOUT
C                     A SURGE POINT
C            IOPERL = 2  AS DEFAULT IN OPERATIONAL MODE
C                 G = GRAVITY
C                C7 = SLIP COEFFICIENT
C               C17 = AIR DENSITY
C               C19 = FRICTIONAL COEFFICIENT
C               C25 = EDDY VISCOSITY
C         MONTH(12) = ACCUMULATED DAYS AT MONTHLY INCREMENTS
C          AMON(12) = FIRST 3 CHARACTERS FOR EACH MONTH
C
C        GENERAL COMMENTS
C           THIS MEMBER  RESIDES IN THE ROOT OVERLAY
C
      COMMON /DUMB8/  IP(4),JP(4),IH(4),JH(4)
      COMMON /DUM88/  IIH(4),JJH(4),IHH(4),JHH(4)
      COMMON /SWTCH/  IOPERL(5)
      COMMON /DELH11/ MAPIN(10),NUMBER
      COMMON /EGTH/   DELS,DELT,G,COR
      COMMON /THRD/   C7,C17,C19,C25
      COMMON /MTH/    MONTH(12),AMON(12)
      CHARACTER*4     AMON
      COMMON /INTV/   INTV(24),KINPUT
C
      DATA IP/0,1,1,0/,JP/0,0,1,1/,IH/-1,0,-1,0/,JH/-1,-1,0,0/
      DATA IIH/0,1,0,-1/,JJH/-1,0,1,0/,IHH/-1,1,0,0/,JHH/0,0,-1,1/
      DATA IOPERL/2,2,2,2,2/
      DATA MAPIN/-9,-6,-3,0,2,4,6,9,2*100/
      DATA DELS,G/5280.,32.2/
      DATA C7,C17,C19,C25 / .006 , .00115 , 3.E-6,.25 /
      DATA MONTH /-1,30,58,89,119,150,180,211,242,272,303,333/
      DATA AMON/'JAN','FEB','MAR','APR','MAY','JUN','JUL','AUG','SEP',
     1     'OCT','NOV','DEC' /
      DATA INTV/1,1,1,1,1,1,1,1,1,1,1,13*2/
      END
C
      SUBROUTINE INTDLT
C        J. CHEN      AUGUST 1980 MDL   IBM 360/195
C
C        PURPOSE
C           THIS SUBROUTINE (INTDLT) ASCERTAINS TIME STEP 'DELT' FOR
C           COMPUTATIONS, INITIALIZES STORM AT SUBSCRIPT 'IBGNT' AND
C           SETS 'ITEND' FOR ENDING COMPUTATIONS AFTER NEAREST STORM
C           APPROACH.
C
C        DATA SET USE
C           NONE
C
C        VARIABLES
C              DELT = TIME INTERVAL FOR EXPLICIT FINITE DIFFERENCING
C              DLTB = THE SMALLEST REDUCTION OF DELT IN THE RUN
C            RLNGTH = DISTANCE OF STORM FROM NEAREST APPROACH
C               JHR = SUBSCRIPT FOR NEAREST APPROACH TO BASIN CENTER
C            P(100) = HOURLY PRESSURE DROPS
C     X(100) Y(100) = HOURLY COMPONENTS FOR STORM POSITIONS
C            R(100) = HOURLY STORM SIZES
C            XMOUTH = X-DISTANCE FROM POLE TO BASIN CENTER
C            YMOUTH = Y-DISTANCE (USUALLY EQUAL TO 0)
C             IBGNT = TABULATED SUBSCRIPT FOR STORM INITIALIZATION
C             ITEND = TABULATED SUBSCRIPT FOR STORM TERMINATION
C              IBGS . = NUMBER OF HOURS FOR COMPUTATIONS BEFORE
C                     NEAREST STORM APPROACH
C              ITMS = NUMBER OF REMAINING HOURS, AFTER NEAREST APPROACH,
C          DLTIN(3) = TIME STEPS FOR 3 CATAGORIES OF STORMS FROM DATA
C
C        GENERAL COMMENTS
C           THIS SUBROUTINE LIES IN OVERLAY 'STRMPR'. IT IS CALLED IN
C           BY SUBROUTINE 'STRMPR'.
C
CC      PARAMETER (IBGS=48, ITMS=48)
      COMMON /EGTH/   DELS,DELT,G,COR
      COMMON /DLTDTA/ DLTIN(3),DLTB
      COMMON /BSN/    PHI,ALTO,ALNO,PHI1,ALT1,ALN1,ALT1C
      COMMON /SPLN/   ALT(15),ALN(15),AX(15),AY(15),RL(15),ANGD(15),
     1                PT(15),RT(15),XLAT(100),YLONG(100),RLNGTH
      COMMON /STRMPS/ X(100),Y(100),P(100),R(100),DIR(100),SP(100)
      COMMON /SSTME/  ISTM,JHR,ITMADV,NHRAD,IBGNT,ITEND
      COMMON /INTV/   INTV(24),KINPUT
      COMMON /ORIGN/  XMOUTH,YMOUTH
      COMMON /INZSTM/ IBGS,ITMS
      COMMON /OPT/ NOPT
      COMMON /FLES/ FLE5,FLE9,FLE8,FLE91,FLE99,FLE10,FLE20,FLE30,FLE1
      CHARACTER*30  FLE5,FLE9,FLE8,FLE91,FLE99,FLE10,FLE20,FLE30,FLE1
C
C        SET INITIAL DELT.
      DELT=DLTIN(1)
C
C        COMPUTE TABULATED SUBSCRIPT FOR STORM INITIALIZATION
       IBGNT=MAX0(JHR-IBGS,1)
      RETURN
      END
      SUBROUTINE INTVAL
C        J. CHEN    AUGUST 1980   MDL   IBM 360/195
C        PURPOSE
C           THIS SUBROUTINE (INTVAL) SETS INITIAL STORM PARAMETERS,
C           SOME FIXED CONSTANTS FOR INITIALIZATION, AND OUTPUT TIMES.
C
C        DATA SET USE
C           NONE
C
C        VARIABLES
C              DELT = TIME INTERVAL FOR EXPLICIT FINITE DIFFERENCING
C     NDLTHR NDLTH2 = NUMBER OF TIME STEPS IN ONE AND 1/2 HOUR
C            IPRTSV = AT THIS TIME PRINT HISTORICAL SURGES (10 POINTS)
C            IPRTAD = ADVANCE 'IPRTSV' BY 'IPRTAD'
C             IBGNT = SUBSCRIPT FOR STORM INITIALIZATION
C     X(100) Y(100) = HOURLY COMPONENTS FOR STORM POSITIONS
C           SP(100) = HOURLY FORWARD STORM SPEEDS
C          DIR(100) = HOURLY STORM AZIMUTHS
C            R(100) = HOURLY STORM SIZES
C              IBGS = NO OF HRS FOR COMPUTATIONS BEFORE NEAREST APPROACH
C              ITMS = NO OF HRS FOR COMPUTATIONS AFTER NEAREST APPROACH
C             MHALT = NUMBER OF TIME STEPS FOR TOTAL COMPUATIONS
C             ITEND = TABULATED SUBSCRIPT FOR STORM TERMINATION
C              IC19 = GROWTH TIME TO MATURE STORM
C              ALTO = LATITUDE OF BASIN CENTER
C            IOPERL = 2 IF IN OPERATIONAL MODE
C            NUMBER = NEAREST HOUR TO NEAREST APPROACH
C
C        GENERAL COMMENTS
C           USE ONLY TIME STEPS DIVISIBLE BY 600 S(10 MIN),OR 900 S(
C           15 MIN).
C           THIS SUBROUTINE LIES IN OVERLAY STRMPR. IT IS CALLED IN BY
C           SUBROUTINE STRMPR.
C
      COMMON /EGTH/   DELS,DELT,G,COR
      COMMON /DELTH/  NDLTHR,NDLTH2,DLTHR,DLTH2
      COMMON /FRTH/   IC12,IC19,BT
      COMMON /PRTSV/  IPRTSV,IPRTAD
      COMMON /SSTME/  ISTM,JHR,ITMADV,NHRAD,IBGNT,ITEND
      COMMON /STRMPS/ X(100),Y(100),P(100),R(100),DIR(100),SP(100)
      COMMON /INZSTM/ IBGS,ITMS
      COMMON /FFTH/   ITIME,MHALT
      COMMON /BSN/    PHI,ALTO,ALNO,PHI1,ALT1,ALN1,ALT1C
      COMMON /DELH11/ MAPIN(10),NUMBER
      COMMON /SWTCH/  IOPERL(5)
      COMMON /MAPTME/ ISTMAP,IADMAP
C
C        NUMBER OF TIME STEPS IN AN HOUR/HALF-HOUR
      NDLTHR=3600./(DELT-.001)
      NDLTH2=NDLTHR/2
      DLTHR =NDLTHR
      DLTH2 =NDLTH2
C
C        NUMBER OF TIME STEPS BETWEEN TIME-HISTORY PRINT-OUT.
      IPRTAD=NDLTHR/4
      IF(MOD(NDLTHR,6).EQ.0) IPRTAD=NDLTHR/6
      IPRTSV=IPRTAD
C
C
C
C        TOTAL RUN TIME
      T=ITEND-IBGNT
      MHALT=(T*3600.+1.)/DELT
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C       INCREASE WIND STRESS FROM ZERO TO FULL STRENGTH IN 10 HRS.
      IC19=10*NDLTHR
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      RETURN
      END
       SUBROUTINE CNSTNT
C        J. CHEN    SEPTEMBER 1980 MDL   IBM 360/195
C
C        PURPOSE
C           THIS SUBROUTINE (CNSTNT) FIXES SOME PROGRAM CONSTANTS.
C
C        DATA SET USE
C           NONE
C
C        VARIABLES
C            X1(50) = STORED ARRAY OF CONSTANTS
C              DELT = TIME INCRMNT FOR FINITE DIFFRNCNG, SET IN 'INTDLT'
C               COR = CORIOLIS (1/SEC), (SET IN 'INTVAL')
C         IMXB JMXB = NO. OF GRID PTS ALONG A RAY/CIRCLE, IN 'CRDRD1'
C       IMXB1 JMXB1 = IMXB LESS 1, JMXB LESS 1
C       IMXB2 JMXB2 = IMXB LESS 2, JMXB LESS 2
C            DEGREE = ANGLE BETWEEN TWO RAYS, READ IN 'CRDRRD1'
C              DELA = VALUE OF 'DEGREE' IN RADIANS
C               RAD = RADIAN MEASURE
C           X1B(10) = STORED ARRAY OF CONSTANTS
C                 G = GRAVITY = 32.2 FT/SEC-2
C       AMAN->AGB13 = FIXED CONSTANT FOR SUBROUTINE 'CHANNL'
C
C        GENERAL COMMENTS
C           THIS SUBROUTINE IS CALLED BY SUBROUTINES INITLZ AND STMVAL.
C
      COMMON /FRST/   X1(50),X12(50)
      COMMON /EGTH/   DELS,DELT,G,COR
      COMMON /DUMB3/  IMXB,JMXB,IMXB1,JMXB1,IMXB2,JMXB2
      COMMON /POLAR/  DEGREE,RMOUTH,AZMTH,DELA
      COMMON /GAMAF/  GAMA,GAMA1,GAMFP(4)
      COMMON /DUMB4/  X1B(10)
      COMMON /DUMB4I/ CORL(90)
      X1B(1)=G*DELT/DELA
      X1B(2)=DELT/DELA
      X1B(4)=DELT*COR
      DO 10 L=1,90
      CORL(L)=2.*(7.292116E-5)*SIN((L-1)*1.74532925E-2)/COR
 10   CONTINUE
      X1B(5)=.5*X1B(2)
      X1B(6)=.5*X1B(1)
C        SET SMOOTHING CONSTSNTS FOR POLAR COORDINATES
      EXPDLA=EXP(DELA)
      COSHLA=.5*(EXPDLA+1./EXPDLA)
      GAMA =EXPDLA/COSHLA-1.
      GAMA1=-GAMA
      GAMFP(1)=0.
      GAMFP(3)=0.
      GAMFP(2)=GAMA
      GAMFP(4)=GAMA1
      RETURN
       END
      SUBROUTINE BTMSTR(ALTO)
C        J. CHEN    SEPTEMBER 1980 MDL   IBM 360/195
C
C        PURPOSE
C           THIS SUBROUTINE (BTMSTR) COMPUTES BOTTOM STRESS COEFFCNTS
C           AT 1 FOOT INTERVALS, FROM 1 TO 300 FEET.
C           THE EQUATIONS ARE IN REFERENCE: JELESNIANSKI,
C           "NUMERICAL COMPUTATIONS OF STORM SURGES WITH BOTTOM
C           STRESS", MONTHLY WEATHER REVIEW, VOL 95, NO. 11,
C           NOV 1967, PP 740-756.
C
C        DATA SET USE
C           NONE
C
C        VARIABLES
C               COR = CORIOLIS PARAMTER
C                 E = EKMAN PARAMETER, DEPTH*SQRT(COR/2*C25)
C               C25 = EDDY VISCOSITY COEFFICIENT, .25 FT**2/SEC
C                C7 = SLIP COEFFICIENT, .006 FT/SEC
C   AR(300) AI(300) = BOTTOM FRICTION COEFFS FOR CORIOLIS TERMS
C   BR(300) BI(300) = BOTTOM FRICTION COEFFS FOR SFC GRAD TERMS
C   CR(300) CI(300) = BOTTOM FRICTION COEFFS FOR SURFACE STRESS TERMS
C
C        GENERAL COMMENTS
C           THIS MEMBER 'BTMSTR' RESIDES IN OVERLAY 'INITLZ'. IT IS
C           CALLED IN BY SUBROUTINE 'INITLZ'.
C           PROGRAM REWRITTEN (1980) BY REPLACING SINH AND COSH
C           WITH SIN AND COS.
C
C
      COMMON /SCND/ AR(900),AI(900),BR(900),BI(900),CR(900),CI(900)
      COMPLEX SIGMAI,CO1,CO2,CO3,CO4,CO5,CO7,CO17
      COMPLEX CO16,CO10,CP
C      REAL*8 E
         DATA C25,C7/.25,.006/
C           CORIOLIS PARAMETER
      COR=2.*(7.292116E-5)*SIN(ALTO*1.74532925E-2)
C      PRINT   5
   5  FORMAT(//'  CONSTANTS FOR SLIP CONDITION AT BOTTOM'//
     1'     E       AR       AI       BR         BI       CR       CI
     2        DEPTH')
      NN=1
      NND=1
C
      DO 100 N=1,900
      A=N
      E=A*SQRT(COR/(2.*C25))
      COR1=COR*A/C7
C         SIGMA=E*(1+I), I*SIGMA=-E+I*E, I=IMAGINARY
      SIGMAI=CMPLX(-E,E)
      CO1=CCOS(SIGMAI)
      CO2=SIGMAI/CSIN(SIGMAI)
      CO1=CO1*CO2
C        CO1=SIGMA/TANH(SIGMA) , CO2=SIGMA/SINH(SIGMA)
C        CO1=SIGMAI/TAN(SIGMAI), CO2=SIGMAI/SIN(SIGMAI), SIGMAI=I*SIGMA
      CO3=1./(CO1+CMPLX(-1.,COR1))
      CO4=CO3*CO3
      CO5=CO2*CO2+CO1
      CP=1.-CO2
      CO7=(.5*CO5-1.)*CO4
      CO17=1./(1.+CO7)
      CO16=(1.+CO3)*CO17
      AR(N)=REAL(CO16)
      AI(N)=AIMAG(CO16)
      BR(N)=REAL(CO17)
      BI(N)=AIMAG(CO17)
      CO10=(1.+CP*CO3)*CO17
      CR(N)=REAL(CO10)
      CI(N)=AIMAG(CO10)
      IF (N.NE.NN) GOTO 100
      IF (N.GT.10) NND=15
      NN=NN+NND
C      PRINT 120,E,AR(N),AI(N),BR(N),BI(N),CR(N),CI(N),N
 120  FORMAT(F6.2,6F10.5,I10)
  100 CONTINUE
      RETURN
       END
      SUBROUTINE CNSDLT
C        J   CHEN  AUG 1982   IBM  360/195
C
C        PURPOSE
C           THIS SUBROUTINE (CNSDLT) FIXES COEFFICIENTS DEPENDING ON
C           TIME STEP FOR VARIATIONS IN THE COURSE OF COMPUTATION..
C
C        DATA SET USE
C           NONE
C
C        VARIABLES
C            X1(50) = STORED ARRAY OF CONSTANTS
C              DELT = TIME INCRMNT FOR FINITE DIFFRNCNG, SET IN 'INTDLT'
C               COR = CORIOLIS (1/SEC), (SET IN 'INTVAL')
C            DEGREE = ANGLE BETWEEN TWO RAYS, READ IN 'CRDRD1'
C              DELA = VALUE OF 'DEGREE' IN RADIANS
C               RAD = RADIAN MEASURE
C           X1B(10) = STORED ARRAY OF CONSTANTS
C                 G = GRAVITY = 32.2 FT/SEC-2
C
C        GENERAL COMMENTS
C           THIS SUBROUTINE IS CALLED BY SUBROUTINES INITLZ AND STMVAL.
C
      COMMON /EGTH/   DELS,DELT,G,COR
      COMMON /FRST/   X1(50),X12(50)
      COMMON /POLAR/  DEGREE,RMOUTH,AZMTH,DELA
      COMMON /DUMB4/  X1B(10)
C
      X1B(1)=G*DELT/DELA
      X1B(2)=DELT/DELA
      X1B(4)=DELT*COR
      X1B(5)=.5*X1B(2)
      X1B(6)=.5*X1B(1)
      RETURN
       END
      SUBROUTINE CRDRD1
C        JYE CHEN            MAY  1990 MDL   IBM 360/195
C        PURPOSE
C           THIS SUBROUTINE (CRDRD1) READS IN BASIN PROJECTION DATA,
C           AS INPUT DATA TO THE SLOSH PROGRAM. EACH BASIN DATA IS ON
C           A SEPARATE FILE, AND CALLED IN BY NAME.
C
C        DATA SET USE
C           FT09F001
C
C        VARIABLES
C               STA = ALPHAMERIC NAME OF A BASIN (10 LETTERS)
C              EBSN =  '$' INDICATES TYPE I ELLIPTIC COORDINATES,
C                      '+' INDICATES TYPE II ELLIPTIC COORDINATES,
C                          OTHERWISE, POLAR COORDINATES
C            DOLLAR = '$' USED FOR ISLAND, OTHERWISE FOR REGULAR BASIN
C             AZMTH = SLANT OF X-AXIS FROM NORTH
C         ALTO ALNO = LATITUDE/LONGITUDE OF BASIN CENTER OR ORIGIN
C            DEGREE = ANGLE BETWEEN TWO RAYS OF A POLAR GRID
C              DELA = RADIAN MEASURE OF DEGREE
C           IRB IRE = TOTAL GRID POINTS IN -+P DIRECTION FROM TANGENT PT
C           NDB NDE = TOTAL GRID POINTS IN -+Q DIRECTION FROM TANGENT PT
C            XMOUTH = X-COORDINATE FROM TANGENT POINT TO (X,Y)=(P,Q)=(0,
C            YMOUTH = Y-COORDINATE FROM TANGENT POINT TO (X,Y)=(P,Q)=(0,
C               PHI = SLANT FROM X-AXIS TO N/S-AXIS
C         IMXB JMXB = TOTAL GRID SPACINGS IN P/Q DIRECTIONS
C       IMXB1 JMXB1 = IMXB LESS 1, JMXB LESS 1
C
C        GENERAL COMMENTS
C
      COMMON /GPRT/   STA
      COMMON /GPRT1/  DOLLAR,EBSN
      COMMON /GPRT2/  EBSN1,EBSN2
      COMMON /POLAR/  DEGREE,RMOUTH,AZMTH,DELA
      COMMON /INDX/   IRB,IRE,NDB,NDE
      COMMON /BSN/    PHI,ALTO,ALNO,PHI1,ALT1,ALN1,ALT1C
      COMMON /DLTDTA/ DLTIN(3),DLTB
      COMMON /ORIGN/  XMOUTH,YMOUTH
      COMMON /ELLIP/  AAXIS,BAXIS,ABQAB
      COMMON /DUMB3/  IMXB,JMXB,IMXB1,JMXB1,IMXB2,JMXB2
      COMMON /EGTH/   DELS,DELT,G,COR
C
      CHARACTER*16  STA
      CHARACTER*2   DOLLAR
      CHARACTER*1   EBSN,EBSN1,EBSN2
      CHARACTER*30  FLE9
C ----------------------------------------------------------------------
C      WRITE(*,*)'  ENTER BASIN FILE '
C      READ (*,'(A)') FLE9
C      OPEN (19,FILE=FLE9)
C-----------------------------------------------------------------------
C        AT COLUMN 12, '$' INDICATES TYPE I ELLIPTIC COORDINATES,
C                      '+' INDICATES TYPE II ELLIPTIC COORDINATES,
C                          OTHERWISE, POLAR COORDINATES
C        AT COLUMN 13, '$ ' INDICATES FOR CLOSED ISLAND (PERIODIC B.C.)
C                      '2$' INDICATES FOR MSY BASIN
C        AT COLUMN 15, '&' INDICATES NO CORNER SMOOTHING
C        AT COLUMN 16, '+' INDICATES OUTPUT OPTION (J-I). DEFAULT (I-J)
C-----------------------------------------------------------------------
      READ (19,'(A10,1X,A1,A2,A1,A1)') STA,EBSN,DOLLAR,EBSN1,EBSN2
C
      READ (19,580) AZMTH
C       TLAT,TLONG NOT BEING USED. USE XMOUTH,YMOUTH, INSTEAD.
      IF (EBSN.EQ.'$'.OR.EBSN.EQ.'+') READ (19,580) TLAT,TLONG
C        LAT/LONG OF BASIN CENTER AND TANGENT POINT ON EARTH
      READ (19,580) ALTO,ALNO
C
C ------    FOR POLAR GRID SYSTEM, READ IN DEGREE OF SPREAD OF RAYS ---
      IF (EBSN.NE.'$'.AND.EBSN.NE.'+') THEN
         READ (19,580) DEGREE
         DELA=DEGREE*1.74532925E-2
       ENDIF
C ------   READ IN DISTANCE FROM BASIN CENTER TO 1ST HYPERBOLA,
C        ANGLE OF 1ST ASYMPTOTE FOR ELLIPTIC GRIDS I  ------
      IF (EBSN.EQ.'$') READ (19,580) DSTNT,ASMPT
C ------   READ IN AAXIS,BAXIS FOR ELLIPTIC GRIDS II  ---
      IF (EBSN.EQ.'+') READ (19,580) AAXIS,BAXIS,GSIZE
C
      READ (19,'(5I5)') IRB,IRE,NDB,NDE,NDB1
C        DEFINE DELTA THETA FOR ELLIPTIC GRIDS I AND II.
      IF (EBSN.EQ.'$') THEN
        IF (NDB.NE.0) THEN
         NDBZ=NDB
         ELSE
         NDBZ=NDE
         ENDIF
C     SPECIAL USE OF ELLIPTIC BASIN I (CARTESIAN APPROX.)
         IF (EBSN1.EQ.'$') NDBZ=NDB1
        DEGREE=(90.-ASMPT)/NDBZ
        DELA=DEGREE*1.74532925E-2
        ASMPT=ASMPT*1.74532925E-2
        AAXIS=0.
        BAXIS=DSTNT/COS(ASMPT)
        ENDIF
      IF (EBSN.EQ.'+') THEN
C        DEGREE=360./(NDB+NDE)
C        DELA=DEGREE*1.74532925E-2
        DELA=2.*GSIZE/(AAXIS+BAXIS)
        DEGREE=DELA/1.74532925E-2
        XMOUTH=AAXIS
        YMOUTH=0
        ENDIF
C
      PHI   =270.+AZMTH
      PHI1=1.74532925E-2*PHI
      ALT1=1.74532925E-2*ALTO
C
C        FORMULATE CORIOLIS, SEC-1
      COR=2.*(7.292116E-5)*SIN(ALT1)
      ALT1C=1.74532925E-2*(90.-ALTO)
      ALN1=1.74532925E-2*ALNO
C
C        READ IN DELTS FOR 3 CATAGORIES OF STORM
      READ (19,100) (DLTIN(I),I=1,3)
 100  FORMAT(3F10.6)
C
      IF (EBSN.NE.'+') READ (19,580) XMOUTH,YMOUTH
C         DEFINE AAXIS=BAXIS= DISATANCE FROM POLE TO TANGENT POINT
C         FOR POLAR GRIDS
      IF (EBSN.NE.'$'.AND.EBSN.NE.'+') THEN
         AAXIS=XMOUTH
         BAXIS=AAXIS
         ENDIF
C
C         DEFINE CONSTANTS FOR GENERAL CODES.
      RMOUTH=.5*(AAXIS+BAXIS)
      ABQAB=(AAXIS-BAXIS)/(AAXIS+BAXIS)
C
      IMXB  =IRB+IRE+1
      JMXB  =NDB+NDE+1
C
C     WRITE (*,*) '     ---- GRID TRANSFORMATION DATA ----'
C     WRITE (*,'(A)') '   AAXIS      BAXIS   DELA(DEG)  IMXB  JMXB  (A-B
C    1)/(A+B)'
C     WRITE (*,111) AAXIS,BAXIS,DELA/1.74532925E-2,IMXB,JMXB,ABQAB
 111  FORMAT(2F10.4,F10.6,2I6,F10.5)
C
      IMXB1 =IMXB-1
      JMXB1 =JMXB-1
      IMXB2=IMXB1-1
      JMXB2=JMXB1-1
C
 580  FORMAT (3F10.6)
      RETURN
       END
      SUBROUTINE TTIMER
C        J. CHEN    JULY 1989     MDL   IBM 360/195
C
C        PURPOSE
C           IT KEEPS TRACK OF TIME AND ACCUMULATES TIME STEPS.
C           IT UPDATES THE STORM POSITION AT EACH TIME STEP ON A
C           CARTESIAN GRID. IT COMPUTES A GROWTH FACTOR TO UPDATE
C           A STORM TO MATURITY.
C
C        DATA SET USE
C           NONE
C
C        VARIABLES
C              IC19 = STORM GROWTH TIME TO MATURITY, SET IN 'INTVAL'
C
C                           Y        (
C                           Y         *(X,Y)
C                           Y          )
C                           Y
C                     PXXXXX0XXXXXXXXXXXXXXXXXXXXXXXX
C                     0   ORIGIN=TANGENT POINT
C                     L     Y    ON EARTH
C                     E     Y
C                           Y
C
C             ITIME = INCREASES BY ONE FOR EACH TIME STEP
C                BT = STORM GROWTH COEFFICIENT, 0<BT<=1
C
C        GENERAL COMMENTS
C           THIS SUBROUTINE RESIDES IN THE ROOT OVERLAY
C
      COMMON /FFTH/   ITIME,MHALT
      COMMON /FRTH/   IC12,IC19,BT
      ITIME=ITIME+1
C      GROWTH FACTOR OF THE STORM FROM ZERO STRENGTH TO MATURITY
      BT=1.
      IF(ITIME.LT.IC19) THEN
      TIME=ITIME
      ANG=TIME/IC19
      BT=.5*(1.-COS(ANG*3.1415927))
         ENDIF
      RETURN
         END

       SUBROUTINE INITLZ
C        J. CHEN    SEPTEMBER 1980MDL   IBM 360/195
C
C        PURPOSE
C           THIS SECOND OVERLAY 'INITLZ' USES SUBROUTINE 'INITLZ' AS A
C           CALLING SEQUENCE FOR 5 SUBROUTINES.
C           THE PURPOSE OF THE OVERLAY IS TO FIX SOME PROGRAM CONSTANTS,
C           SOME POLAR GRID CONSTANTS; TO SET BATHY/TOPO AND
C           BARRIER HEIGHTS, AND TO INITIALIZE A TIMING SEQUENCE.
C
C        DATA SET USE
C           NONE
C
C        VARIABLES
C         ALTO (DEG)= LATITUDE OF BASIN CENTER OR ORIGIN
C
C        GENERAL COMMENTS
C     THIS SUBROUTINE RESIDES IN OVERLAY 'INITLZ'. IT IS CALLED
C           BY THE ROOT OVERLAY AND NOT BY ANY OTHER SUBROUTINE.
C
CC      PARAMETER (IBGS=48, ITMS=48)
      COMMON /BSN/    PHI,ALTO,ALNO,PHI1,ALT1,ALN1,ALT1C
      COMMON /SSTME/  ISTM,JHR,ITMADV,NHRAD,IBGNT,ITEND
      COMMON /SWTCH/  IOPERL(5)
      COMMON /GPRT1/  DOLLAR,EBSN
      COMMON /INZSTM/ IBGS,ITMS
      COMMON /DATUM/  SEADTM,DTMLAK
      COMMON /OPT/  NOPT
      COMMON /CESAV/  IPN(120),IPL(120),KHSPT
      COMMON /CESAV1/ STATS(120)
      COMMON /FLES/ FLE5,FLE9,FLE8,FLE91,FLE99,FLE10,FLE20,FLE30,FLE1
      CHARACTER*30  FLE5,FLE9,FLE8,FLE91,FLE99,FLE10,FLE20,FLE30,FLE1
      CHARACTER*2     DOLLAR
      CHARACTER*1     EBSN
      CHARACTER*30    STATS,FLESTA
C
C          READ BASIN PROJECTION DATA
      CALL CRDRD1
C
C        JHR=61, MAX 60 HOURS  SPIN-UP + HINCAST  TIME.
C        48 HR ANALYSES FIELDS EVERY 6 HR. (DEFAULT)
C        36 HR FORECAST FIELDS EVERY 3 HR.
      IBGS=48
      ITMS=96
      JHR=61

C----------- OPTION TO OVERRIDE THE SPIN-UP TIME (# HOURS)  -------
C 
      READ(21,'(2I3,4X,2F5.2)',END=2000) IBGS,ITMS,SEADTM,DTMLAK
C      READ(21,'(2I3)',END=2000) IBGS,ITMS
 2000 CONTINUE
C       CLOSE (1)
       DTMLAK=SEADTM
       IBGNT=MAX0(JHR-IBGS,1)
      ITEND=JHR+ITMS

C ---------------- GAGE STATIONS USED FOR OUTPUT ------------
CC      READ (14,'(I3)') KHSPT
C     WRITE (*,'(A,I3)') ' NO. OF STA',KHSPT
CC      DO 2001 N=1,KHSPT
CC      READ (14,'(A30,18X,2I4)') STATS(N),IPN(N),IPL(N)
C     WRITE (*,'(A28,20X,2I4)') STATS(N),IPN(N),IPL(N)
 2001 CONTINUE
C      CLOSE (4)
C
       CALL CRDRD2
       CALL DEPSFC
C
       CALL BTMSTR(ALTO)
       CALL SETDST
       CALL INTDLT
       CALL INTVAL
       CALL CNSTNT
      RETURN
       END
      SUBROUTINE SETDST
C        J. CHEN    SEPTEMBER 1980 MDL   IBM 360/195
C
C        PURPOSE
C           THIS SUBROUTINE (SETDST) COMPUTES SIZE SPECIFICS OF THE
C           POLAR GRID, AS WELL AS THE ORIENTATION OF THE (P,Q) POLAR
C           SYSTEM VERSUS THE (X,Y) SYSTEM.  THE END RESULTS ARE USED
C           FOR DRIVING FORCES TRANSFORMED FROM THE (X,Y) TO (P,Q)
C           SYSTEMS.  PRINT OUT OF POLAR SPECIFICS ARE EFFECTED.
C
C        DATA SET USE
C           NONE
C
C        VARIABLES
C              DELA = RADIAN MEASURE BETWEEN CONTIGUOUS RAYS
C         IMXB JMXB = TOTAL GRID SPACINGS, SET IN 'CRDRD1'
C       IMXB1 JMXB1 = IMXB LESS 1, JMXB LESS 1
C           IRB JRB = TOTAL GRID POINTS IN +-DIRCTN FROM ORIGIN
C           NDB NDE = TOTAL GRID POINTS IN +-Q DIRECTION FROM ORIGIN
C            RMOUTH = DISTANCE FROM POLE TO BASIN ORIGIN IN MILES
C        SIDEHZ(  ) = ARC LENGTH BETWEEN CONTIGUOUS GRID POINTS IN FT
C        ARESQ(  ) = AREA OF SQUARES, BETWEN CONCENTRIC CIRCLES IN FT-2
C        SIDEVT(  ) = DISTANCE BETWEEN CONCENTRIC CIRCLES IN MILES IN FE
C            DEGREE = ANGLE BETWEEN 2 CONTIGUOUS RAYS
C          ANGD(  ) = ANGLE BETWEEN P-AXIS AND ANY RAY, IN DEGREES
C               RAD = RADIAN MEASURE
C COST(  ) SINT(  ) = RAY PRJCTNS THRU MOMENTUM POINTS, ON (X,Y) AXIS
C COSL(  ) SINL(  ) = RAY PRJCTNS THRU SURGE POINTS, ON (X,Y) AXIS
C
C        GENERAL COMMENTS
C           THIS SUBROUTINE RESIDES IN OVERLAY 'INITLZ'. IT IS
C           CALLED IN BY SUBROUTINE 'INITLZ'.
C
C
      PARAMETER (M$=169,N$=239,L$=1700,NCT$=1700,LC$=L$)
      COMMON /DUMB5/  ZB(M$,N$),ZBM(M$,N$)
      COMMON /DUMB7/  UB(M$,N$),VB(M$,N$),HB(M$,N$)
      COMMON /SCRTCH/ HSUB(M$,N$)
      COMMON /DUMB10/ HMX(M$,N$)
      COMMON /ELPDST/ ELPCT(M$),ELPDT(M$),ELPCL(M$),ELPDL(M$)
      COMMON /ELPDS2/ ELPCL2(M$),ELPDL2(M$),SINL2(N$),SINT2(N$)
      COMMON /PLRPST/ SINT(N$),COST(N$),SINL(N$),COSL(N$)
      COMMON /SETMIX/ IS(N$),MS(N$),MF(N$)
      COMMON /SETMX1/ IE(N$),ME(N$)
      COMMON /ITREE/  ITREE(M$,N$)
      COMMON /KTREE/  KTREE(L$)
      CHARACTER*1     ITREE,KTREE
      COMMON /KSPCON/ KSKP(M$,N$)
      CHARACTER*1     KSKP
      COMMON /DUMB18/ ISQR(L$),JSQR(L$),ISIDE(L$),HWEIR(L$),
     1                ZBMIN(L$)
      COMMON /ACT1D/  FLWSQR(L$),F1DACT(L$)
      CHARACTER*1     F1DACT
      COMMON /DUMCUT/ CUTL(NCT$),CUTLI(NCT$),CUTLE(NCT$)
C
      COMMON /POLAR/  DEGREE,RMOUTH,AZMTH,DELA
      COMMON /ELLIP/  AAXIS,BAXIS,ABQAB
      COMMON /DUMB3/  IMXB,JMXB,IMXB1,JMXB1,IMXB2,JMXB2
      COMMON /INDX/   IRB,IRE,NDB,NDE
      COMMON /EGTH/   DELS,DELT,G,COR
      COMMON /FLWCPT/ NSQRS,NSQRW,NSQRWC,NPSS,NCUT
      COMMON /CHNL/   IPT0(2,5),JPT0(2,5),IPTL(2,5),JPTL(2,5),ENTEXT(5)
      COMMON /CHNARA/ ARESQ0(5),ARESQL(5)
      DIMENSION DUMMY(180)
C        COMPUTING RADII OF CONCENTRIC CIRCLES IN ST MILES
      DO 100 I=1,IMXB
      X=I-1-IRB
      RFEET1=EXP(X*DELA)
      GRID=RMOUTH*RFEET1*5280.
      Z1=1./(RFEET1*RFEET1)
      ELPCT(I)=(1.+ABQAB*Z1)*GRID
      ELPDT(I)=(1.-ABQAB*Z1)*GRID
 100  CONTINUE
C
C        COMPUTING AZIMUTHS FOR EACH RAY FROM STORM CENTER
      DO 130 J=1,JMXB
      X=J-1-NDB
      ANGL=X*DELA
C        AZIMUTHS ON MOMENTUM POINTS
      COST(J)=COS(ANGL)
      SINT(J)=SIN(ANGL)
      SINT2(J)=SINT(J)*SINT(J)
C        AZIMUTHS ON HEIGHT POINTS
      ANGH=ANGL+.5*DELA
      COSL(J)=COS(ANGH)
      SINL(J)=SIN(ANGH)
      SINL2(J)=SINL(J)*SINL(J)
 130  CONTINUE
C
C        COMPUTING RADIAL DISTANCES TO HT PTS, IN FEET
C        DETERMINING HEIGHT POINT JACOBIANS FOR CONTTY EQUATION IN FEET
      DO 150 I=1,IMXB1
      X=I-1-IRB
      X=X+0.5
      RFEET1=EXP(X*DELA)
      GRID=RMOUTH*RFEET1*5280.
C      GRID2=1./(RMOUTH*RFEET1*5280.)**2
      Z1=1./(RFEET1*RFEET1)
      ELPCL(I)=(1.+ABQAB*Z1)*GRID
      ELPDL(I)=(1.-ABQAB*Z1)*GRID
      ELPCL2(I)=ELPCL(I)*ELPCL(I)
      ELPDL2(I)=ELPDL(I)*ELPDL(I)
 150  CONTINUE
C        APPROX. THE AREAS OF THE ENTRANCE AND EXIT SQUARES OF CHANNELS.
      IF (NPSS.NE.0) THEN
      DO 310 NPASS=1,NPSS
      Z1=0.
      Z2=0.
      DO 300 K=1,2
      I=IPT0(K,NPASS)
      J=JPT0(K,NPASS)
      II=IPTL(K,NPASS)
      JJ=JPTL(K,NPASS)
      Z1=Z1+ELPDL2(I)+(ELPCL2(I)-ELPDL2(I))*SINL2(J)
      Z2=Z2+ELPDL2(II)+(ELPCL2(II)-ELPDL2(II))*SINL2(JJ)
 300  CONTINUE
      ARESQ0(NPASS)=.5*Z1*DELA*DELA
      ARESQL(NPASS)=.5*Z2*DELA*DELA
 310  CONTINUE
        ENDIF
C
C        PRINT OUT X-DISTANCES ALONG Q=0, IN MILES
C      WRITE (*,190)
      DO 140 I=1,IMXB
      DUMMY(I)=ELPCT(I)/5280.
 140  CONTINUE
C      WRITE (*,170)(DUMMY(I),I=1,IMXB)
C        PRINT OUT Y-DISTANCES ALONG P=0, IN MILES
C      WRITE (*,210)
      DO 135 J=1,JMXB1
      X=J-1-NDB
      ANGL=ELPDT(1)*SIN(X*DELA)/5280.
 135  DUMMY(J)=ANGL
C      WRITE (*,170)(DUMMY(J),J=1,JMXB1)
 160  FORMAT (10X,' RADIUS OF CIRCLES FROM ORIGIN, IN MILES')
 170  FORMAT(F7.1,19F6.1)
 180  FORMAT(3X,20I6)
 190  FORMAT(10X,' GRID DISTANCES ALONG X-AXIS (ALONG Q=0)')
 200  FORMAT(20F6.2)
 210  FORMAT(10X,' GRID DISTANCES ALONG Y-AXIS (ALONG P=0)')
      RETURN
      END
      SUBROUTINE CRDRD2
      PARAMETER (NBK$=1500)
      PARAMETER (M$=169,N$=239,L$=1700,NCT$=1700,LC$=L$)
      COMMON /DUMB5/  ZB(M$,N$),ZBM(M$,N$)
      COMMON /DUMB7/  UB(M$,N$),VB(M$,N$),HB(M$,N$)
      COMMON /SCRTCH/ HSUB(M$,N$)
      COMMON /DUMB10/ HMX(M$,N$)
      COMMON /ELPDST/ ELPCT(M$),ELPDT(M$),ELPCL(M$),ELPDL(M$)
      COMMON /ELPDS2/ ELPCL2(M$),ELPDL2(M$),SINL2(N$),SINT2(N$)
      COMMON /PLRPST/ SINT(N$),COST(N$),SINL(N$),COSL(N$)
      COMMON /SETMIX/ IS(N$),MS(N$),MF(N$)
      COMMON /SETMX1/ IE(N$),ME(N$)
      COMMON /ITREE/  ITREE(M$,N$)
      COMMON /KTREE/  KTREE(L$)
      CHARACTER*1     ITREE,KTREE
      COMMON /KSPCON/ KSKP(M$,N$)
      CHARACTER*1     KSKP
      COMMON /DUMB18/ ISQR(L$),JSQR(L$),ISIDE(L$),HWEIR(L$),
     1                ZBMIN(L$)
      COMMON /ACT1D/  FLWSQR(L$),F1DACT(L$)
      CHARACTER*1     F1DACT
      COMMON /DUMCUT/ CUTL(NCT$),CUTLI(NCT$),CUTLE(NCT$)
C
      COMMON /CS1DA/  COS1D(LC$),SIN1D(LC$)
      COMMON /BANK2/  DELCUT(LC$),BANK(LC$,2)
C
      COMMON /DUMB3/  IMXB,JMXB,IMXB1,JMXB1,IMXB2,JMXB2
      COMMON /BSN/    PHI,ALTO,ALNO,PHI1,ALT1,ALN1,ALT1C
      COMMON /SPLN/   ALT(15),ALN(15),AX(15),AY(15),RL(15),ANGD(15),
     1                PT(15),RT(15),XLAT(100),YLONG(100),RLNGTH
      COMMON /DUMBSV/ SAVEH(10),IPN(10),JPN(10),STATS(10)
      COMMON /PRTOVL/ NSHEET,JPGS(6),JPGF(6),IPGS(6),IPGF(6)
      COMMON /FLWCPT/ NSQRS,NSQRW,NSQRWC,NPSS,NCUT
      COMMON /DLTDTA/ DLTIN(3),DLTB
      COMMON /CHNLD/  CHLNH(5),CHWTH(5),CHDPH(5)
      COMMON /CHNL/   IPT0(2,5),JPT0(2,5),IPTL(2,5),JPTL(2,5),ENTEXT(5)
      COMMON /IPRX/   IX60
      COMMON /GPRT1/  DOLLAR,EBSN
      CHARACTER*1  EBSN
C
      CHARACTER*16 STATS
      DIMENSION   IANG(L$)
      CHARACTER*16 AAA
      DIMENSION   INW(7),JNW(7),HW(7)
      DIMENSION   IBNK(NBK$),JBNK(NBK$),ISBNK(NBK$),BKI(NBK$),
     1            BKE(NBK$),WCHN(NBK$)
      CHARACTER*2   AC,BC,IX60,XA1,YA1,ZA1,DOLLAR,SIG$
      DATA XA1/'X'/,YA1/'Y'/,ZA1/'Z'/,SIG$/'$'/
C
C
C        INITIALIZE ZBM(I,J) FOR ENTIRE GRIDS
      DO 120 J=1,JMXB
      DO 120 I=1,IMXB
 120  ZBM(I,J)=-300.
C
C         READING IN BARRIER HEIGHTS AT SELECTED GRID POINTS
      READ (19,520) INWP,AC
       IF(AC.EQ.'X') INWP=INWP+1000
       IF(AC.EQ.'Y') INWP=INWP+2000
C     WRITE (*,520) INWP
      IF (INWP.EQ.0) GO TO 220
      INWX=(INWP-1)/7+1
      I7X=INWP-(INWX-1)*7
      I7=7
      DO 210 IP=1,INWX
      IF(IP.EQ.INWX) I7=I7X
      READ (19,550) (INW(I),JNW(I),HW(I),I=1,I7)
      DO 200 N=1,I7
      I=INW(N)
      J=JNW(N)
C     TEST FOR I,J BUNDS
C      IF (I.GE.IMXB.OR.J.GE.JMXB) WRITE (*,190) I,J
 190  FORMAT(' MIS-PUNCHED BARRIER LOCATIONS',2I5)
 200  ZBM(I,J)=HW(N)
 210  CONTINUE
C
C        COMPLETE DEFINING ZBM( , ) IN SUBROUTINE 'DEPSFC'
C
C         READ IN DATA FOR 1-DIM FLOW, (I,J) POINTS, SIDES 1-4
 220  CONTINUE
C
C       INITIALIZE DATA FOR REGULAR 1D FLOW
      DO 230 I=1,L$
      BANK(I,1)=0.
      BANK(I,2)=0.
      HWEIR(I)=0.
      DELCUT(I)=1.
 230  CONTINUE
C       AC='X' IF 1D-FLOW POINTS EXCEED 1000 .
      READ (19,520) NSQRS,AC,BC
C     WRITE (*,520) NSQRS,AC,BC
      IF (AC.EQ.'X') NSQRS=1000+NSQRS
      IF (AC.EQ.'Y') NSQRS=2000+NSQRS
C
      IF (NSQRS.GT.L$) THEN
C      WRITE (*,*)' NO. OF 1D FLOW POINTS ',NSQRS
C      WRITE (*,*)' EXCEEDS THE MAXIMUM ALLOWED. MAX= ',L$
      STOP
      ENDIF
      IF (NSQRS.EQ.0) GO TO 300
      READ (19,540) (ISQR(I),JSQR(I),KTREE(I),ISIDE(I),IANG(I),
     1      I=1,NSQRS)
C
C        FOR BC='X' ENTER SPECIAL 1D FLOW SQUARES,BE SURE OF THE SAME
C        ORDER AS IN 1D FLOW
      IF (BC.NE.'X') GO TO 300
C
C       READ IN BANK DATA
      READ (19,520) NBNKS,AC
      IF (AC.EQ.'X') NBNKS=1000+NBNKS
C     WRITE (*,520) NBNKS
      READ (19,250)(IBNK(I),JBNK(I),AC,ISBNK(I),BKI(I),BKE(I),
     1 WCHN(I),I=1,NBNKS)
 250  FORMAT (4(2I3,A1,I1,2F3.0,F4.2))
C
C       REPOSITION BANKS WITH RESPECT TO 1D FLOW SQUARES
C       LIST OF BANK POINTS MUST BE IN THE SAME ORDER AS 1D FLOW PT.
      LLL=1
      DO 290 I=1,NBNKS
      II=IBNK(I)
      JJ=JBNK(I)
      K=ISBNK(I)
      LL=LLL
      DO 280 L=LL,NSQRS
      ISET=ISQR(L)
      IF (ISET.NE.II) GO TO 260
      JSET=JSQR(L)
      IF (JSET.NE.JJ) GO TO 260
      KSET=ISIDE(L)
      IF (KSET.NE.K) GO TO 260
      LLL=L+1
      BANK(L,1)=BKI(I)
      BANK(L,2)=BKE(I)
      DELCUT(L)=WCHN(I)
      GO TO 290
 260  IF (L.NE.NSQRS) GO TO 280
C     WRITE (*,270) IBNK(I),JBNK(I),ISBNK(I)
 270  FORMAT('   THE SPECIAL FLOW POINT FOR BANKS AT I=',I3,' J=',
     1I3,' SIDE=',I3,' DOES NOT EXIST IN THE 1D FLOW REGION')
 280  CONTINUE
 290  CONTINUE
C
C         READ IN DATA FOR RAISED WEIRS FOR 1-DIM FLOW
 300  READ (19,520) NWEIRS,AC
      IF (AC.EQ.'X') NWEIRS=NWEIRS+1000
      IF (NWEIRS.EQ.0) GO TO 320
      READ (19,530)(ISQR(NSQRS+I),JSQR(NSQRS+I),KTREE(NSQRS+I),
     1 ISIDE(NSQRS+I),HWEIR(NSQRS+I),I=1,NWEIRS)
      DO 310 I=1,NWEIRS
 310  IANG(I+NSQRS)=0
 320  NSQRW=NSQRS+NWEIRS
C
C
C         READ IN DATA FOR 10 SELECTED POINTS (SURGE HISTORY)
      READ (19,630) (STATS(I),I=1,10)
C     WRITE (*,630) (STATS(I),I=1,10)
      READ (19,610) (IPN(I) ,I=1,10)
      READ (19,610) (JPN(I) ,I=1,10)
C
C         READ IN DATA FOR CHANNELS
      READ (19,520) NPSS
      IF (NPSS.EQ.0) GO TO 340
      READ (19,650) (CHLNH(I),I=1,NPSS)
      READ (19,680) (CHWTH(I),I=1,NPSS)
      READ (19,690) (CHDPH(I),I=1,NPSS)
      READ (19,700) ((IPT0(I,J),I=1,2),J=1,NPSS)
      READ (19,700) ((JPT0(I,J),I=1,2),J=1,NPSS)
      READ (19,700) ((IPTL(I,J),I=1,2),J=1,NPSS)
      READ (19,700) ((JPTL(I,J),I=1,2),J=1,NPSS)
 340  CONTINUE
C
C         READ IN DATA FOR FLOW ACROSS CUTS
      READ (19,520) NCUT,AC
      NSQRWC=NSQRW+NCUT
      IF (NSQRWC.GT.L$) THEN
C      WRITE (*,*)' NO. OF 1D FLOW POINTS + CUTS',NSQRWC
C      WRITE (*,*)' EXCEEDS THE MAXIMUM ALLOWED. MAX= ',L$
      STOP
        ENDIF
C
      IF (NCUT.EQ.0) GO TO 390
C
      DO 380 I=1,NCUT
      M=I+NSQRW
      READ (19,370)ISQR(M),JSQR(M),KTREE(M),ISIDE(M),BANK(M,1),
     1 BANK(M,2),CUTL(I),CUTLI(I),CUTLE(I),IANG(M),HWEIR(M)
 370  FORMAT(2I3,A1,I1,2F4.0,3F5.2,I3,F5.1)
      IF (CUTL(I).EQ.0.) CUTL(I)=1.
      IF (CUTLI(I).EQ.0.) CUTLI(I)=1.
      IF (CUTLE(I).EQ.0.) CUTLE(I)=1.
 380  CONTINUE
 390  CONTINUE
C
C        COMPUTE SINE AND COSINE OF THE ANGLE THAT IS USED TO
C        ROTATE WIND STRESS DIRECTION FROM ACTUAL TO MODEL'S 1D FLOW
C        DIRECTIONS.
      DO 410 N=1,NSQRWC
      IF (IANG(N).EQ.0) GO TO 400
      ANG=IANG(N)*1.74532925E-2
      COS1D(N)=COS(ANG)
      SIN1D(N)=SIN(ANG)
      GO TO 410
 400  COS1D(N)=1.
      SIN1D(N)=0.
 410  CONTINUE
 420  CONTINUE
C        READ IN DATA FOR MANGROVE POINTS OR TREE LOCATIONS (MGPT=-1)
      READ (19,520) MGPT
      IF (MGPT.EQ.0) GOTO 470
C        READ IN LOCATIONS OF TREE POINTS (0 OR 1)
      DO 460 J=1,JMXB
      READ (19,450) (ITREE(I,J),I=1,IMXB)
 450  FORMAT(60A1)
 460  CONTINUE
C        SET UP BOUNDARY VALUES TO ZBM( , ) TO BE CHECKED FOR OPEN
C        BOUNDARY CALCULATIONS.
       IF (DOLLAR.NE.SIG$) THEN
      DO 160 J=1,JMXB,JMXB1
      DO 160 I=1,IMXB
      ZBM(I,J)=100.
 160  CONTINUE
          ENDIF
      DO 180 I=1,IMXB,IMXB1
      DO 180 J=1,JMXB
      ZBM(I,J)=100.
 180  CONTINUE
C
C         READ IN DATA FOR PRINT OUT OF SHEETS OF MAX SURGES
C        OPTIONAL FORMATS TO ACCOMMODATE DIFFERENT BASINS DATA
C        FORMAT-- I3,A1
C        IX60--' ' FOR DISPLAYING I=1,IXB,WHERE IXB < 60
C              'X' FOR DISPLAYING I=1,IXB,WHERE IXB >=60
C              'Y' FOR I=1,IPGF(I)
C              'Z' FOR I=IPGS(I),IPGF(I)  MOST GENERAL
 470  READ (19,520) NSHEET,IX60
C     WRITE (*,520) NSHEET,IX60
      READ (19,610) (JPGS(I),I=1,NSHEET)
C     WRITE (*,610) (JPGS(I),I=1,NSHEET)
      READ (19,610) (JPGF(I),I=1,NSHEET)
C     WRITE (*,610) (JPGF(I),I=1,NSHEET)
      DO 480 N=1,NSHEET
      IPGS(N)=1
      IPGF(N)=MIN0(59,IMXB1)
      IF (IX60.EQ.'X') IPGF(N)=IMXB1
 480  CONTINUE
      IF (IX60.EQ.'Z') THEN
       READ (19,610) (IPGS(I),I=1,NSHEET)
C      WRITE (*,610) (IPGS(I),I=1,NSHEET)
        ENDIF
      IF (IX60.EQ.YA1.OR.IX60.EQ.ZA1) THEN
        READ (19,610) (IPGF(I),I=1,NSHEET)
C       WRITE (*,610) (IPGF(I),I=1,NSHEET)
        ENDIF
C
C         READING IN DEPTH VALUES
      IMX=(IMXB-1)/15+1
      DO 500 K=1,IMX
      IST=(K-1)*15+1
      IFN=MIN0(IST+14,IMXB)
      DO 490 J=1,JMXB
      READ (19,710) (ZB(I,J),I=IST,IFN)
C     IF (J.EQ.JMXB) WRITE (*,710) (ZB(I,J),I=IST,IFN),AAA
C      DO 459 I=IST,IFN
C      IF (ZB(I,J).GE.900.) ZB(I,J)=2500.
C 459  CONTINUE
 490  CONTINUE
 500  CONTINUE
      IF (DOLLAR.EQ.'$') THEN
      DO 510 I=1,IMXB
 510  ZB(I,JMXB)=ZB(I,1)
      ELSE
      DO 512 I=1,IMXB
 512  ZB(I,JMXB)=ZB(I,JMXB1)
      ENDIF
C
 520  FORMAT (I3,2A1)
 530  FORMAT (5(2I3,A1,I1,F5.1))
 540  FORMAT (6(2I3,A1,I1,I3))
 550  FORMAT (7(I3,I3,F3.0))
 560  FORMAT (A10,2X,A1)
 570  FORMAT (F10.6)
 580  FORMAT (2F10.6)
 590  FORMAT (4I5)
 600  FORMAT (22I3)
 610  FORMAT (15I3)
 620  FORMAT (45I1)
 630  FORMAT (6A10)
 640  FORMAT (15F3.0)
 650  FORMAT (5F5.1)
 660  FORMAT (13F5.0)
 670  FORMAT (A4,3I3,1X,3I3)
 680  FORMAT (5F6.0)
 690  FORMAT (5F5.0)
 700  FORMAT (10I3)
 710  FORMAT (15F4.0,A16)
      RETURN
       END
      SUBROUTINE DEPSFC
C        J. CHEN    SEPTEMBER 1980 MDL   IBM 360/195
C
C        PURPOSE
C           THIS SUBROUTINE (DEPSFC) INITIALIZES INLAND SURFACE WATERS
C           TO LAKE DATUM (SEA WATER TO SEA DATUM IS DONE IN SUBROU-
C           TINE 'INILHT'). BARRIERS AND THEIR HEIGHTS ARE PLACED AT
C           SELECTED MOMENTUM GRID POINTS. FOR 1-DIM FLOW, DEPTHS AND
C           LATERAL BOUNDARIES ARE SET. BOUNDARIES ARE TEMPORARILY
C           SET AS A 100 FOOT WALL. PRINT OUTS ARE MADE OF THE DEPTH
C           FIELD, AND HIGHEST LAND VALUE AT MOMENTUM GRID POINTS
C           ARE SET.
C
C        DATA SET USE
C           NONE
C
C        VARIABLES
C         IMXB JMXB = TOTAL GRID SPACINGS, SET IN 'CRDRD1'
C       IMXB1 JMXB1 = IMXB LESS 1, JMXB LESS 1
C            DOLLAR = '$' FOR ISLAND, OTHERWISE FOR REGULAR BASIN
C         ZB(  ,  ) = BATHY/TOPO HEIGHTS,    RELATIVE TO NGVD
C        ZBM(  ,  ) = HIGHEST ELEVATION AT A MOMENTUM POINT
C                     OR, VALUES INDICATING BOUNDARY TYPES
C             NSQRW = NUMBER OF SQUARES WITH 1-DIM FLOW
C    ISQR JSQR(   ) = SUBSCRIPTS FOR SQUARES WITH 1-DIM FLOW
C        ISIDE(   ) = THE SIDE OF A SQUARE TO PASS 1-DIM FLOW
C     IIH(4) JJH(4) = SHIFT SUBSCRIPTS TO 4 MOMENTUM POINTS ON A SQUARE
C        HWEIR(   ) = LOWESTT BARRIER ELEVATION AT ISIDE(   )
C     IZ1(4) IZ2(4) = SET P SUBSCRIPTS ON 4 SIDES OF ISIDE(   )
C     JZ1(4) JZ2(4) = SET Q SUBSCRIPTS ON 4 SIDES OF ISIDE(   )
C        ZBMIN(   ) = LOWEST DEPTH OF TWO SQUARES BOUNDING ISIDE(30
C              MGPT = TOTAL NUMBER OF MANGROVE PTS
C     IMG( ) JMG( ) = I/J COORDINATES OF MANGROVE PTS
C      ITREE(  ,  ) = 0 FOR NO TREE, 1 FOR TREE ON A GRID POINT
C
C        GENERAL COMMENTS
C           THIS SUBROUTINE RESIDES IN OVERLAY 'INITLZ'. IT IS
C           CALLED IN SUBROUTINE 'INITLZ'.
C
      PARAMETER (M$=169,N$=239,L$=1700,NCT$=1700,LC$=L$)
      COMMON /DUMB5/  ZB(M$,N$),ZBM(M$,N$)
      COMMON /DUMB7/  UB(M$,N$),VB(M$,N$),HB(M$,N$)
      COMMON /SCRTCH/ HSUB(M$,N$)
      COMMON /DUMB10/ HMX(M$,N$)
      COMMON /ELPDST/ ELPCT(M$),ELPDT(M$),ELPCL(M$),ELPDL(M$)
      COMMON /ELPDS2/ ELPCL2(M$),ELPDL2(M$),SINL2(N$),SINT2(N$)
      COMMON /PLRPST/ SINT(N$),COST(N$),SINL(N$),COSL(N$)
      COMMON /SETMIX/ IS(N$),MS(N$),MF(N$)
      COMMON /SETMX1/ IE(N$),ME(N$)
      COMMON /ITREE/  ITREE(M$,N$)
      COMMON /KTREE/  KTREE(L$)
      CHARACTER*1     ITREE,KTREE
      COMMON /KSPCON/ KSKP(M$,N$)
      CHARACTER*1     KSKP
      COMMON /DUMB18/ ISQR(L$),JSQR(L$),ISIDE(L$),HWEIR(L$),
     1                ZBMIN(L$)
      COMMON /ACT1D/  FLWSQR(L$),F1DACT(L$)
      CHARACTER*1     F1DACT
      COMMON /DUMCUT/ CUTL(NCT$),CUTLI(NCT$),CUTLE(NCT$)
C
      COMMON /CS1DA/  COS1D(LC$),SIN1D(LC$)
      COMMON /BANK2/  DELCUT(LC$),BANK(LC$,2)
C
      COMMON /DUMB3/  IMXB,JMXB,IMXB1,JMXB1,IMXB2,JMXB2
      COMMON /DUMB8/  IP(4),JP(4),IH(4),JH(4)
      COMMON /DUM88/  IIH(4),JJH(4),IHH(4),JHH(4)
      COMMON /BSN/    PHI,ALTO,ALNO,PHI1,ALT1,ALN1,ALT1C
      COMMON /SPLN/   ALT(15),ALN(15),AX(15),AY(15),RL(15),ANGD(15),
     1                PT(15),RT(15),XLAT(100),YLONG(100),RLNGTH
      COMMON /BCPTS/  NBCPT,ISH(500),JSH(500)
      COMMON /FLWCPT/ NSQRS,NSQRW,NSQRWC,NPSS,NCUT
      COMMON /SWTCH/  IOPERL(5)
      COMMON /GPRT1/  DOLLAR,EBSN
      CHARACTER*2     DOLLAR
      CHARACTER*1     EBSN
      DIMENSION IZ1(4),IZ2(4),JZ1(4),JZ2(4)
      DATA IZ1/0,1,0,0/,  IZ2/0,1,1,1/, JZ1/1,1,0,1/, JZ2/0,0,0,1/
C       ZB=1500 FT STATIC HEIGHT BOUNDARY CONDITIONS.
       DO 401 J=2,JMXB1
       DO 401 I=2,IMXB1
       ZMX=AMAX1(-ZB(I-1,J-1),-ZB(I,J-1),-ZB(I-1,J),-ZB(I,J))
C -------- DEEP OCEAN BOUNDARY IS SET AT 1500 FT  ------------
C -------- OR, STATIC HEIGHT BUNDARY CONDITION IS SPECIFIED THERE ------
       IF (ZMX.LT.-1499.) ITREE(I,J)='6'
C     NO FLOODING, OVERTOPPING IS ALLOWED
      IF (ZMX.GT.0.) ITREE(I,J)='4'
C
 401   CONTINUE
C      DETERMINE THE ACTIVE/NONACTIVE SQUARES FOR HEIGHT COMPUTATIONS
       DO 400 J=1,JMXB1
       DO 400 I=1,IMXB1
       IF (ZB(I,J).LE.-35.) GOTO 1417
       KSKP(I,J)='1'
       KCT=1
C       SPECIAL TRUNCATED SQUARES IN DEEP OCEAN
       IF (ITREE(I,J).EQ.'6') KCT=KCT+1
       IF (ITREE(I+1,J).EQ.'6') KCT=KCT+1
       IF (ITREE(I,J+1).EQ.'6') KCT=KCT+1
       IF (ITREE(I+1,J+1).EQ.'6') KCT=KCT+1
       GOTO (400,415,415,415,417),KCT
 415   KSKP(I,J)='2'
       GOTO 400
 417   ZB(I,J)=-67.
 1417  KSKP(I,J)='0'
 400   CONTINUE
       DO 421 J=1,JMXB
 421   KSKP(IMXB,J)='0'
       DO 420 I=1,IMXB
 420   KSKP(I,JMXB)=KSKP(I,1)
C
C
C        ISOLATE HIGHEST LAND SQUARE ABOUT A TRANSPORT POINT
C        FOR ISLAND ZB( ,JMXB)=ZB( ,1)
      J2=JMXB1
      IF (DOLLAR.EQ.'$')  J2=JMXB
      DO 130 I=2,IMXB
      DO 130 J=2,J2
      ZMX=AMAX1(-ZB(I-1,J-1),-ZB(I,J-1),-ZB(I-1,J),-ZB(I,J))
      IF (ZBM(I,J).GE.ZMX) GO TO 130
      IF (ZBM(I,J).NE.-300.) GO TO 110
      ZBM(I,J)=ZMX
      GO TO 130
 110  WRITE (*,120)I,J,ZBM(I,J)
 120  FORMAT ('   ERROR IN BARRIER HEIGHT AT POSITION I=',I5,' J=',I5,
     1 '  HW=',F5.0)
      ZBM(I,J)=ZMX
 130  CONTINUE
      IF (DOLLAR.EQ.'$') THEN
        DO 132 I=1,IMXB
        ZBM(I,1)=ZBM(I,JMXB)
 132    CONTINUE
        ENDIF
C        IF ZBM GREATER THAN 35 FT IN INTERIOR POINTS, MOMENTURM PTS
C        ARE NONACTIVE.
       DO 40 J=2,J2
       DO 40 I=2,IMXB1
       IF (ZBM(I,J).GE.35.) ITREE(I,J)='4'
       IF (ITREE(I,J).EQ.'5') ITREE(I,J)='3'
 40    CONTINUE
C
C       REDEFINE MS(J). FOR EACH J, FIRST I THE COMPUTATION
C       STARTS IN MOMNTM (ZBM LESS THAN 35 FT).
      DO 50 J=2,JMXB
      MS(J)=1
      DO 60 I=1,IMXB
      IF (ZBM(I,J).LT.35.) GOTO 55
 60   CONTINUE
      I=I-1
 55   MS(J)=I
 50   CONTINUE
      IF (DOLLAR.EQ.'$') THEN
         MS(1)=MS(JMXB)
         ELSE
         MS(1)=MS(2)
         ENDIF
C       DEFINE IS(J). FOR EACH J, FIRST I THE COMPUTATION
C       STARTS IN CONTINUITY.
      DO 10 J=1,JMXB1
      IS(J)=1
      DO 20 I=1,IMXB1
      IF (KSKP(I,J).EQ.'1') GOTO 15
 20   CONTINUE
      I=I-1
 15   IS(J)=I
 10   CONTINUE
      IF (DOLLAR.EQ.'$')THEN
         IS(JMXB)=IS(1)
         ELSE
         IS(JMXB)=IS(JMXB1)
         ENDIF
C
C        DEFINE LOWER LIMIT OF I FOR MOMN. POINTS
      DO 450 J=2,JMXB1
      ME(J)=IMXB
      DO 460 II=2,IMXB1
      I=IMXB-II+1
      IF (ZBM(I,J).GE.35.) GOTO 460
      IF (ITREE(I,J).NE.'6') GOTO 455
 460  CONTINUE
      I=I+1
 455  ME(J)=I
 450  CONTINUE
      ME(1)=IMXB1
      ME(JMXB)=IMXB1
C        DEFINE LOWER LIMIT OF I FOR CONTTY POINTS
C
      DO 610 J=1,JMXB1
      IE(J)=IMXB1
      DO 620 II=1,IMXB1
      I=IMXB1-II+1
      IF (KSKP(I,J).EQ.'1') GOTO 615
 620  CONTINUE
      I=I+1
 615  IE(J)=I
 610  CONTINUE
      IE(JMXB)=IE(1)
C     PRINT 1111,(IS(J),J=1,JMXB)
C     PRINT 1111,(MS(J),J=1,JMXB)
C 1111 FORMAT(1X,15I3)
C     PRINT 1111,(IE(J),J=1,JMXB)
C     PRINT 1111,(ME(J),J=1,JMXB)
C
C
      DO 42 J=2,JMXB
      MF(J)=IMXB
      DO 45 II=1,IMXB1
      I=IMXB-II
      IF (ITREE(I,J).EQ.'1'.OR.ITREE(I,J).EQ.'3') GOTO 48
 45   CONTINUE
 48   MF(J)=I
 42   CONTINUE
      IF (DOLLAR.EQ.'$') THEN
        MF(1)=MF(JMXB)
        ELSE
        MF(1)=MF(2)
        ENDIF
C     IF NO BANKS, THEN SET BANK HEIGHTS TO BARRIER HEIGHT
      IF (NSQRWC.EQ.0) GO TO 250
      DO 200 L=1,NSQRWC
      IF (BANK(L,1).NE.0.) GO TO 190
      I=ISQR(L)
      J=JSQR(L)
      ISS=1
 180  BANK(L,ISS)=AMIN1(ZBM(I,J),ZBM(I,J+1),ZBM(I+1,J),ZBM(I+1,J+1))
       GO TO (190,200),ISS
 190  IF (BANK(L,2).NE.0.) GO TO 200
      K=ISIDE(L)
      I=ISQR(L)+IHH(K)
      J=JSQR(L)+JHH(K)
      ISS=2
      GO TO 180
 200  CONTINUE
C
C        SET UP DEPTHS AND SIDE HEIGHTS FOR 1-DIM FLOW AND RIVER BANKS
C        SET UP SPECIAL RAISED WEIRS OR SILLS FOR 1-DIM FLOW
C        AND ZBMIN FOR FLOW CUTS IN SUBROUTINE CRDRD2.
      IF (NSQRW.EQ.NSQRS) GOTO 2333
      N1=NSQRS+1
      DO 232 L=N1,NSQRW
      I=ISQR(L)
      J=JSQR(L)
      K=ISIDE(L)
      I1=I+IZ1(K)
      I2=I+IZ2(K)
      J1=J+JZ1(K)
      J2=J+JZ2(K)
      IF (HWEIR(L).EQ.0.) HWEIR(L)=AMIN1(ZBM(I1,J1),ZBM(I2,J2))
 232  CONTINUE
      DO 1232 N=N1,NSQRW
      I=ISQR(N)
      J=JSQR(N)
      K=ISIDE(N)
      I1=I+IZ1(K)
      I2=I+IZ2(K)
      J1=J+JZ1(K)
      J2=J+JZ2(K)
      ITREE(I1,J1)='Q'
      ITREE(I2,J2)='Q'
 1232 CONTINUE
 2333 CONTINUE
      DO 3222 J=2,JMXB1
      DO 3222 I=2,IMXB1
      IF (ITREE(I,J).EQ.'Q') THEN
      ZBM(I,J)=36.
      ITREE(I,J)='4'
      ENDIF
 3222 CONTINUE
      DO 240 L=1,NSQRWC
      I=ISQR(L)
      J=JSQR(L)
      ZINT=-ZB(I,J)
      K=ISIDE(L)
      II=I+IHH(K)
      JJ=J+JHH(K)
      ZEXT=-ZB(II,JJ)
C      FOR WEIR HEIGHTS NOT A SILL
      IF (HWEIR(L).NE.0.) GO TO 230
 220  HWEIR(L)=AMAX1(ZINT,ZEXT)
 230  I1=I+IZ1(K)
      I2=I+IZ2(K)
      J1=J+JZ1(K)
      J2=J+JZ2(K)
      IZBX=AMIN1(ZBM(I1,J1),ZBM(I2,J2))
      ZBMIN(L)=IZBX
 240  CONTINUE
C
      IF (NSQRW.EQ.NSQRS) GOTO 250
      DO 2221 N=N1,NSQRW
      I=ISQR(N)
      J=JSQR(N)
      K=ISIDE(N)
      I1=I+IZ1(K)
      I2=I+IZ2(K)
      J1=J+JZ1(K)
      J2=J+JZ2(K)
C     WRITE (*,2222) I,J,K,HWEIR(N),ZBMIN(N),ZBM(I1,J1),ZBM(I2,J2)
C 2222 FORMAT(3I5,4F6.1)
 2221 CONTINUE
C
 250  CONTINUE
C      DEFINE GRID COORDINATES FOR STATIC HEIGHT BOUNDARY CONDITION.
      N=0
      DO 1200 I=1,IMXB1
      DO 1200 J=1,JMXB1
      IF (KSKP(I,J).NE.'2') GOTO 1200
      N=N+1
      ISH(N)=I
      JSH(N)=J
      KSKP(I,J)='0'
 1200 CONTINUE
      NBCPT=N
C
      RETURN
       END
      SUBROUTINE FLW1DM
C        J. CHEN   SEPTEMBER 1980 MDL   IBM 360/195
C
C        PURPOSE
C           THIS SUBROUTINE (FLW1DM) COMPUTES 1-DIM FLOW ACROSS CHANNEL
C           SIDES OF INTERLOCKING SQUARES. WATER PASSES ONLY THRU 'OPEN'
C           SIDES OF A 1-DIM FLOW CHANNEL, SEE DIAGRAM BELOW:
C
C        DATA SET USE
C           NONE
C
C        VARIABLES
C           AI(800) = BOTTOM FRICTION, AT ONE FOOT INTERVALS
C             NSQRW = NO OF SQUARES WITH 1-DIM FLOW POSSIBLE
C  ISQR(L) JSQRL(L) = I/J-SUBSCRIPT FOR SQUARE WITH 1-DIM FLOW
C        HB(  ,   ) = SURGE FIELD
C        ISIDE(300) = SIDE OF FLOW RELATIVE TO INTERIOR 1-DIM FLOW SQRE
C     IHH(4) JHH(4) = SHIFT I/J-SUBSCRIPT FOR CONTIGUOUS SQUARE
C         HINT HEXT = INTERIOR/EXTERIOR SURGE FOR SFC GRADIENT
C       HINT1 HEXT1 = INTERIOR/EXTERIOR SURGE HEIGHTS ON SQUARES
C        ZBMIN(300) = LOWEST BARRIER VALUE AT END OF AN 'OPEN' SIDE
C                ZL = DUMMY VALUE FOR ZBMIN(300)
C         ZB(  ,  ) = DEPTH FIELD
C        HWEIR(300) = HIGHEST DEPTH  OF INT/EXT SQUARES, OR WEIR HEIGHT
C              HOVT = DUMMY VALUE FOR HWEIR(300)
C     FLWSQR(300  ) = TRANSPORT ACROSS AN 'OPEN' SIDE, TWO TIME LEVELS
C        GRIDR2(  ) = JACOBEAN FOR POLAR COORDINATES
C           SIGA(4) = SIGN FOR FLOW ACROSS SIDES OF SQUARES
C          DPH1D(2) = DEPTH+SURGE AT INT/EXT SQUARES, PRESENT TIME
C     CHNGHI CHNGHE = SURGE INCREASE ON INT/EXT SQRES, FROM 'OPEN' SIDE
C       DEPLI DEPLE = TOTAL DEPTH OF INT/EXT SQRES, FUTURE TIME
C               CII = ACCUMULATED WATER FROM 'OPEN' SIDES, INTERIOR SQR
C              SUBE = SURGE FROM EXT SQR TO PREVENT WATER BLW INT SQR
C         )    CEE = ACCUMULATED WATER FROM 'OPEN' SIDES, EXTERIOR SQR
C              SUBI = SURGE FROM INT SQR TO PREVENT WATER BLW EXT SQR
C            ID DPH = MEAN TOTAL DEPTH OF INT/EXT SQUARES, PRESENT TIME
C               DZL = HEIGHT OF HWEIR ABOVE MEAN DEPTHS
C           FXB FYB = COMPONENTS OF SURFACE STRESS
C            X1B(N) = FIXED COEFFICIENTS FOR FINITE DIFFERENCING
C   BR(300) BI(300) = FRICTION COEFFICIENTS AT 1-FOOT INTERVALS
C              JUMP = 2 FOR STRESS COMPUTATIONS, 1 FOR NO STRESS
C
C        GENERAL COMMENTS
C           THIS IS A VERSION FOR 2-STEP SCHEME (JULY 1981).
C           THIS SUBROUTINE RESIDES IN OVERLAY 'CMPUTE'. IT IS CALLED
C           IN BY SUBROUTINE 'CMPUTE'. SUBROUTINE 'FLW1DM' CALLS IN
C
      PARAMETER (M$=169,N$=239,L$=1700,NCT$=1700,LC$=L$)
      COMMON /DUMB5/  ZB(M$,N$),ZBM(M$,N$)
      COMMON /DUMB7/  UB(M$,N$),VB(M$,N$),HB(M$,N$)
      COMMON /SCRTCH/ HSUB(M$,N$)
      COMMON /DUMB10/ HMX(M$,N$)
      COMMON /ELPDST/ ELPCT(M$),ELPDT(M$),ELPCL(M$),ELPDL(M$)
      COMMON /ELPDS2/ ELPCL2(M$),ELPDL2(M$),SINL2(N$),SINT2(N$)
      COMMON /PLRPST/ SINT(N$),COST(N$),SINL(N$),COSL(N$)
      COMMON /SETMIX/ IS(N$),MS(N$),MF(N$)
      COMMON /SETMX1/ IE(N$),ME(N$)
      COMMON /ITREE/  ITREE(M$,N$)
      COMMON /KTREE/  KTREE(L$)
      CHARACTER*1     ITREE,KTREE
      COMMON /KSPCON/ KSKP(M$,N$)
      CHARACTER*1     KSKP
      COMMON /DUMB18/ ISQR(L$),JSQR(L$),ISIDE(L$),HWEIR(L$),
     1                ZBMIN(L$)
      COMMON /ACT1D/  FLWSQR(L$),F1DACT(L$)
      CHARACTER*1     F1DACT
      COMMON /DUMCUT/ CUTL(NCT$),CUTLI(NCT$),CUTLE(NCT$)
C
      COMMON /CS1DA/  COS1D(LC$),SIN1D(LC$)
      COMMON /BANK2/  DELCUT(LC$),BANK(LC$,2)
C
      COMMON /EGTH/   DELS,DELT,G,COR
      COMMON /SCND/   AR(900),AI(900),BR(900),BI(900),CR(900),CI(900)
      COMMON /DUMB4/  X1B(10)
      COMMON /DUMB8/  IP(4),JP(4),IH(4),JH(4)
      COMMON /DUM88/  IIH(4),JJH(4),IHH(4),JHH(4)
      COMMON /MF/     DPH,HPD(4),ID
      COMMON /FLWCPT/ NSQRS,NSQRW,NSQRWC,NPSS,NCUT
C
      DIMENSION       SIGA(4),DPH1D(2)
C
      DATA SIGA/-1.,1.,-1.,1./
C
      DO 100 L=1,NSQRWC
      I=ISQR(L)
      J=JSQR(L)
      K=ISIDE(L)
      II=I+IHH(K)
      JJ=J+JHH(K)
      HSUB(I,J)=0.
      HSUB(II,JJ)=0.
 100  F1DACT(L)='F'
C
CXXXXXXXXX BEGIN 1ST 'LOOP'. TREAT 2-DIM FLOW (FROM 'CONTTY') ONTO
CXXXXXXXXX A 1-DIM SQUARE AS A SOURCE.
C
      DO 510 L=1,NSQRWC
C
      IF (FLWSQR(L).EQ.0.) GO TO 510
C        ISOLATE A PARTICULAR SQUARE AND ITS SURFACE HEIGHT
      I=ISQR(L)
      J=JSQR(L)
C
C        DETERMINE WHICH OF 4 SIDES HAS FLOW AND ISOLATE ADJACENT SQUARE
C        AND ITS SURFACE HEIGHT
      K=ISIDE(L)
      II=I+IHH(K)
      JJ=J+JHH(K)
C
C        START CONTINUITY EQUATION FOR ONE SIDE OF A SQUARE
C        COMPUTE TOTAL FLOW IN/OUT OF INTERIOR SQUARE
 110  Z=ELPDL2(I)+(ELPCL2(I)-ELPDL2(I))*SINL2(J)
      Z1=ELPDL2(II)+(ELPCL2(II)-ELPDL2(II))*SINL2(JJ)
      CHNGHI=-SIGA(K)*X1B(5)*FLWSQR(L)/Z
      CHNGHE=-CHNGHI*Z/Z1
C
C     DETERMINE 1-DIM, OR EXPANSION, FLOW
      KCUT=1
      IF (L.GT.NSQRW) KCUT=2
      GO TO (120,130), KCUT
C
C     CHECK FOR EXISTENCE OF BANKS
 120  IF (DELCUT(L).GE.1.) GO TO 470
      GO TO 140
 130  LL=L-NSQRW
      IF (CUTLI(LL).GE.1..AND.CUTLE(LL).GE.1.) GO TO 420
C
C      TWO CASES EXIST: OUTFLOW VERSUS INFLOW
C      CASE 1: CHNGHI NEG (OUTFLOW), CHNGHE POS (INFLOW),  LEAP=1
C      CASE 2: CHNBHI POS (INFLOW),  CHNGHE NEG (OUTFLOW), LEAP=2
 140   LEAP=1
       IF (CHNGHI.GT.0.) LEAP=2
       GO TO (192,194),LEAP
 192   CHNG=CHNGHI
       CSTRT=CHNGHI
        HT=HB(I,J)
        BK=BANK(L,1)
        IF(KCUT.EQ.2) CUTW=CUTLI(LL)
        GO TO (200,300),LEAP
 194   CHNG=CHNGHE
       CSTRT=CHNGHE
       HT=HB(II,JJ)
       BK=BANK(L,2)
       IF(KCUT.EQ.2) CUTW=CUTLE(LL)
       GO TO (300,200),LEAP
C******LOGIC FOR NEGATIVE WATER CHANGES (OUTFLOW)
 200   IF (HT.GT.BK) GO TO 230
        GO TO (210,220),KCUT
C       INITIAL WATER BELOW BANK AND DESCENDING
 210   CHNG=CHNG
       GO TO 280
 220   CHNG=CHNG*CUTL(LL)/CUTW
       GO TO 272
C      INITIAL WATER ABOVE BANK AND DESCENDING
 230   DELB=HT-BK
       GO TO (240,260),KCUT
C      WATER DOES NOT DESCEND BELOW BANK, PURE 1-D FLOW
 240   IF(-CHNG*DELCUT(L).GT.DELB) GO TO 250
       CHNG=CHNG*DELCUT(L)
       GO TO 272
C      WATER DESCENDS BELOW BANK, PURE 1-D FLOW
 250   CHNG=CHNG-DELB+DELB/DELCUT(L)
       GO TO 272
C      WATER DOES NOT DESCEND BELOW BANK, CUT-TYPE FLOW
 260   IF (-CHNG*CUTL(LL).GT.DELB) GO TO 270
       CHNG=CHNG*CUTL(LL)
       GO TO 272
C      WATER DESCENDS BELOW BANK, CUT-TYPE FLOW
 270   CHNG=(CHNG*CUTL(LL)+DELB)/CUTW-DELB
 272   CONTINUE
 280   GO TO (290,295),LEAP
 290   CHNGHI=CHNG
       HSUB(I,J)=HSUB(I,J)+ABS(CSTRT-CHNG)
       GO TO 194
 295   CHNGHE=CHNG
       HSUB(II,JJ)=HSUB(II,JJ)+ABS(CSTRT-CHNG)
       GO TO 192
C$$$$$$LOGIC FOR POSITIVE WATER CHANGES (INFLOW)
 300   IF (HT.LT.BK) GO TO 330
       GO TO (310,320),KCUT
C      INITIAL WATER ABOVE BANK AND ASCENDING
 310   CHNG=CHNG*DELCUT(L)
       GO TO 372
 320   CHNG=CHNG*CUTL(LL)
       GO TO 372
C      INITIAL WATER BELOW BANK AND ASCENDING
 330   DELB=BK-HT
       GO TO (340,360),KCUT
C      WATER DOES NOT ASCEND ABOVE BANK, PURE 1-D FLOW
 340   IF (CHNG.GT.DELB) GO TO 350
       CHNG=CHNG
       GO TO 380
C      WATER ASCENDS ABOVE BANK, PURE 1-D FLOW
 350   CHNG=DELB+(CHNG-DELB)*DELCUT(L)
       GO TO 372
C      WATER DOES NOT ASCEND ABOVE BANK, CUT-TYPE FLOW
 360   IF (CHNG*CUTL(LL).GT.DELB*CUTW) GO TO 370
       CHNG=CHNG*CUTL(LL)/CUTW
       GO TO 372
C      WATER ASCENDS ABOVE BANK, CUT-TYPE FLOW
 370   CHNG=CHNG*CUTL(LL)-DELB*CUTW+DELB
 372   CONTINUE
 380   GO TO (390,395),LEAP
 390   CHNGHE=CHNG
       HSUB(II,JJ)=HSUB(II,JJ)+ABS(CSTRT-CHNG)
       GO TO 470
 395   CHNGHI=CHNG
       HSUB(I,J)=HSUB(I,J)+ABS(CSTRT-CHNG)
       GO TO 470
 420   CHNGHI=CHNGHI*CUTL(LL)/CUTLI(LL)
       CHNGHE=CHNGHE*CUTL(LL)/CUTLE(LL)
 470   HBJ=HB(I,J)+CHNGHI
      HBJJ=HB(II,JJ)+CHNGHE
C        FINISH CONTINUITY EQUATION FOR ONE SIDE OF A SQUARE BY
C        TESTING IF WATER LIES BELOW SQUARES
 480  DEPLI=HBJ+ZB(I,J)
      DEPLE=HBJJ+ZB(II,JJ)
      IF (DEPLI.GE.0.) GO TO 490
C        WATER LIES BELOW INTERIOR SQUARE (DEPLI IS NEGATIVE)
      CII=CHNGHI-DEPLI
      SUBE=-CII*Z/Z1
      HB(I,J)=-ZB(I,J)
      HB(II,JJ)=AMAX1(HB(II,JJ)-CHNGHE+SUBE,-ZB(II,JJ))
      FLWSQR(L)=FLWSQR(L)*CII/CHNGHI
      GO TO 510
C        WATER LIES BELOW EXTERIOR SQUARE (DEPLE IS NEGATIVE)
 490  IF(DEPLE.GE.0.) GO TO 500
      CEE=CHNGHE-DEPLE
      SUBI=-CEE*Z1/Z
      HB(I,J)=AMAX1(HB(I,J)-CHNGHI+SUBI,-ZB(I,J))
      HB(II,JJ)=-ZB(II,JJ)
      FLWSQR(L)=FLWSQR(L)*CEE/CHNGHE
      GO TO 510
C
C        BOTH INTERIOR AND EXTERIOR SQUARES HAVE WATER ABOVE TERRAIN
C        AFTER UPDATES.
 500  HB(I,J)=HBJ
      HB(II,JJ)=HBJJ
C
      F1DACT(L)='T'
 510  CONTINUE
C***********************************************************************
C        END OF HEIGHT UPDATES BY COMPUTED FLOW ACCROSS SIDES
C***********************************************************************
C
C        USE UPDATED HEIGHTS COMPUTE FLOW ACCROSS SIDES.
C
C        FRICTION CHANGE BETWEEN 1-2 FEET OF TOTAL DEPTH, USED
C        FOR EXTRAPOLATION FROM 1 TO 0 FT.
C        SET PARAMETERS FOR FRCPNT ARGUMENTS, LAKE WINDS, WITH STATIC
C        PRESSURE FORCE.
      SLPA=AI(2)-AI(1)
      NCATG=1
      NPLS=1
C
      DO 630 L=1,NSQRWC
C        ISOLATE A PARTICULAR SQUARE AND ITS SURFACE HEIGHT
      I=ISQR(L)
      J=JSQR(L)
      HINT=HB(I,J)
C        DETERMINE WHICH OF 4 SIDES HAS FLOW AND ISOLATE ADJACENT SQUARE
C        AND ITS SURFACE HEIGHT
      K=ISIDE(L)
      II=I+IHH(K)
      JJ=J+JHH(K)
      HEXT=HB(II,JJ)
C        ISOLATE MINIMUM ZBM ON SIDE CONNECTING TWO SQUARES
      ZL=ZBMIN(L)
C        TEST IF SFC LIES ABOVE ZL, E.G., 2-DIM FLOW IN BOTH SQUARES
      IF (AMIN1(HINT,HEXT).GE.ZL) GO TO 520
C        TEST IF BOTH SQUARES ARE DRY
      IF (HINT.EQ.-ZB(I,J).AND.HEXT.EQ.-ZB(II,JJ)) GO TO 520
C        DETERMINE HT OF WEIR OR HIGHEST DEPTH OF TWO SQUARES
      HOVT=HWEIR(L)
C        TEST IF ANY SQUARE(S) OVERTOP THE WEIR AT PRESENT TIME
      IF(AMAX1(HINT,HEXT).GT.HOVT) GO TO 530
 520  FLWSQR(L)=0.
      HSUB(I,J)=0.
      HSUB(II,JJ)=0.
      GO TO 630
 530  CONTINUE
C        COMPUTE MOMENTUM ON ONE SIDE OF A SQUARE
      HINT1=HB(I,J)
      HEXT1=HB(II,JJ)
      HINT=AMIN1(AMAX1(HINT1,HOVT),ZL)
      HEXT=AMIN1(AMAX1(HEXT1,HOVT),ZL)
      JUMP=2
      IF (AMAX1(HINT,HEXT).GE.ZL) JUMP=1
      HEAD=HEXT-HINT
      HEAD=SIGA(K)*HEAD
      DPH1D(1)=HINT+ZB(I,J)
      DPH1D(2)=HEXT+ZB(II,JJ)
      DPH=.5*(DPH1D(1)+DPH1D(2))
      ID=DPH
      ID=MIN(899,MAX0(ID,1))
      BIOBR=BI(ID)/BR(ID)
      IF (NSQRW.EQ.NSQRS) GOTO 1492
          IF (L.GT.NSQRS.AND.L.LE.NSQRW) BIOBR=0.
 1492 CONTINUE
      HCOFF=BR(ID)+BI(ID)*BIOBR
      FLWT=-X1B(1)*HCOFF*DPH*HEAD
         IF (JUMP.EQ.1)  GOTO 589
C
C        AVERAGE SHELTERING COEFF IF WATER BECOME LESS THAN 1 FT ABOVE
C        BED,OR 1 FT BELOW CHANNEL BANK,OR BARRIER HEIGHT.
      IF (KTREE(L).EQ.'T') THEN
         CSHLTR=(DPH/50.)**2
       ELSE
         CSHINT=1.-DIM(1.,DIM(HINT1,HOVT))
         CSHEXT=1.-DIM(1.,DIM(HEXT1,HOVT))
         CSH2I=1.-DIM(1.,DIM(ZL,HINT1))
         CSH2E=1.-DIM(1.,DIM(ZL,HEXT1))
         CSHLTR=.5*AMIN1(CSHINT+CSHEXT,CSH2I+CSH2E)
       ENDIF
C
      III=I
      JJJ=J
      GO TO (116,112,113,114),K
C        COMPUTE (X,Y) POSITION OF FORCE POINTS ON POLAR GRID
 112  III=I+1
 116  XR0=COSL(JJJ)
      YR0=SINL(JJJ)
      ELPCZ=ELPCT(III)
      ELPDZ=ELPDT(III)
      GO TO 115
 114  JJJ=J+1
 113  XR0=COST(JJJ)
      YR0=SINT(JJJ)
      ELPCZ=ELPCL(III)
      ELPDZ=ELPDL(III)
 115  CONTINUE
      XR=XR0*ELPCZ
      YR=YR0*ELPDZ
C
C      CALL FRCPNT(XR,YR,FXB,FYB,NCATG,CSHLTR,NPLS)
C
      GO TO (516,512,513,514),K
 516  III1=I
      JJJ1=J+1
      GOTO 1121
 512  III=III+1
      III1=III+1
      JJJ1=JJJ+1
      GOTO 1121
 513  III1=III+1
      JJJ1=JJJ
      GOTO 1121
 514  JJJ=JJJ+1
      III1=III+1
      JJJ1=JJJ+1
 1121 CONTINUE
      CALL FRCPNT(III,JJJ,FXB,FYB,NCATG,CSHLTR,NPLS)
      CALL FRCPNT(III1,JJJ1,FXB1,FYB1,NCATG,CSHLTR,NPLS)
       FXB=.5*(FXB+FXB1)
       FYB=.5*(FYB+FYB1)
C        CHANGE FORCING FROM (X,Y) SYSTEM TO IMAGE PLANE
      XR=XR0*ELPDZ
      YR=YR0*ELPCZ
      FXBP=XR*FXB+YR*FYB
      FYBP=XR*FYB-YR*FXB
C        CORRECT STRESS DIRECTION ACCORDING TO ACTURAL FLOW
      IF(COS1D(L).NE.1.) THEN
      FXB=COS1D(L)*FXBP-SIN1D(L)*FYBP
      FYB=SIN1D(L)*FXBP+COS1D(L)*FYBP
      ELSE
      FXB=FXBP
      FYB=FYBP
      ENDIF
C
      FRC1=.75*FXB
      FRC2=.75*FYB
      IF (K.EQ.3.OR.K.EQ.4) THEN
      Z=FRC2
      FRC2=-FRC1
      FRC1=Z
      ENDIF
C
      FRC12T= DELT*(FRC1+BIOBR*FRC2)
C           ENDED STRESS COMPUTATIONS.
      GOTO 590
C
C        INITIALIZE FLOW WITH GRAVITY AND WITHOUT SURFACE STRESS
 589  FRC12T=0.
C
 590   CONTINUE
C        TAKE AVERAGE OF 'AI' ON TWO SQUARES
      AINTRP=0.
      DO 620 KK=1,2
      IDD=AMIN1(299.,DPH1D(KK))
      IF (IDD.GT.1) GO TO 600
      AITR=AI(1)+SLPA*(DPH1D(KK)-1.)
      GO TO 610
 600  SLPAI=AI(IDD+1)-AI(IDD)
      DPHI=IDD
      AITR=AI(IDD)+SLPAI*(DPH1D(KK)-DPHI)
 610  AINTRP=AINTRP+AITR
 620  CONTINUE
      AINTRP=.5*AINTRP
      FLWCOF=X1B(4)*(AINTRP-AR(ID)*BIOBR)+1.
      FLW=FLWT+FLWCOF*FLWSQR(L)+FRC12T
      FWD=HSUB(II,JJ)-HSUB(I,J)
      IF (FWD.EQ.0.) GO TO 628
      FD2=-SIGA(K)*FWD*X1B(2)/(DPH*DPH)
      FD2=ABS(FD2)
      Z=ELPDL2(I)+(ELPCL2(I)-ELPDL2(I))*SINL2(J)
      Z1=ELPDL2(II)+(ELPCL2(II)-ELPDL2(II))*SINL2(JJ)
      XY=.5*(1./Z+1./Z1)
      FD2=FD2*XY
      FLW1=FLW
      DO 627 ITR=1,20
      FLW1=FLW/(1.+FD2*ABS(FLW1))
 627  CONTINUE
      FLWSQR(L)=FLW1
      GO TO 630
 628  FLWSQR(L)=FLW
 630  CONTINUE
C++++++++END MOMENTUM COMPUTATIONS FOR 1-DIM FLOW+++++++++++++++++++++++
C
C--------END MAIN 'LOOP'------------------------------------------------
      RETURN
       END
      SUBROUTINE GEOWND
C        J. CHEN               MAY 1992 MDL   IBM 360/195
C
C        PURPOSE
C           THIS SUBROUTINE CALULATES THE PRESSURE GRADIENT FORCINGS, AN
C           CONVERTS ALL FORCINGS IN PROPER UNITS USED IN SLSOH EQUATION
C           WIND SPEEDS ARE CONVERTED FROM METER/SEC TO FT/SEC, PRESSURE
C           FROM MB TO HEIGHT OF WATER COLUMN, PRESSURE GRADIENTS FROM
C           MB/DEG TO GRADIENT CURRENT UNIT, FT/SEC INCORPERATING CORIOL
C           PARAMETERS AT DIFFERENT LATITUDES.
C
C        DATA SET USE
C           NONE
C
C        VARIABLES
C     IU,IV,XIU,XIV = U/V SURFACE WIND COMPONENTS IN LAT/LON GRID
C            IP,XIP = PRESSURE IN TERMS OF WEIGHT OF A WATER COLUMN
C IUG,IVG,XIUG,XIVG = CURRENT UNIT DUE TO ATM. PRESSURE GRADIENTS
C        YLT,YLG(,) = LAT/LON OF THE SLOSH GRID OF A BASIN
C             RHOWG = WATER DENSITY*GRAVITY (MB FT**-1)
C
C        GENERAL COMMENTS
C          IU,IV,IP,IUG,IVG ARE INPUT DATA IN A SPECIAL UNITS, VALUES
C          TIME 10 STORED AS IN INTEGERS. IT NEED NOT BE SO IN FINAL
C          OPERATIONAL PROGRAM.
C
      PARAMETER (IMX$=36,JMX$=50)
      COMMON /UVP/ XIU(IMX$,JMX$,2),XIV(IMX$,JMX$,2),XIP(IMX$,JMX$,2),
     1               XIUG(IMX$,JMX$,2),XIVG(IMX$,JMX$,2)
      COMMON /IUV/   IU(IMX$,JMX$),IV(IMX$,JMX$),IP(IMX$,JMX$),
     1               IUG(IMX$,JMX$),IVG(IMX$,JMX$)
      COMMON /DATINC/ XLAT1,XLON1,DELLT,DELLN,IMXC,JMXC
      INTEGER*2    IU,IV,IP,IUG,IVG
      DIMENSION    ZP(IMX$,JMX$)
C
C     CNSTN REDUCES PRESSURE GRADIENT/CORIOLIS FROM MB/DEG*SEC TO
C         M/S(WIND SPEED UNIT ), GEOSTROPHIC WINDS.
C -------------------------------------------------------------------
      DATA CNSTN/4.769134/
C
C         CALCULATE THE GEOSTROPHIC WIND FROM FINITE DIFFERENCE
C         PRESSURE GRADIENT IN LAT/LONG COORDINATES.
C
        DO 220 J=1,JMXC
        DO 220 I=IMXC,1,-1
        ZP(I,J)=IP(I,J)
 220    CONTINUE
C
C
        DO 201 I=2,IMXC-1
        ZZ=SIN((XLAT1-I+1)*1.74532925E-2)
        Z1=COS((XLAT1-I+1)*1.74532925E-2)
        DO 201 J=2,JMXC-1
        XVG=.25*(ZP(I+1,J+1)-ZP(I+1,J-1)+ZP(I-1,J+1)-ZP(I-1,J-1))
        XUG=.25*(ZP(I+1,J+1)-ZP(I-1,J+1)+ZP(I+1,J-1)-ZP(I-1,J-1))
        IUG(I,J)=CNSTN*XUG/(ZZ*Z1)
        IVG(I,J)=CNSTN*XVG/ZZ
  201  CONTINUE
C
C     EXTENDED TO BOUNDARY VALURES
        K=-1
        DO 203 J=1,JMXC,JMXC-1
        K=-K
        DO 203 I=1,IMXC
        IUG(I,J)=IUG(I,J+K)
        IVG(I,J)=IVG(I,J+K)
 203    CONTINUE
        K=-1
        DO 204 I=1,IMXC,IMXC-1
        K=-K
        DO 204 J=1,JMXC
        IUG(I,J)=IUG(I+K,J)
        IVG(I,J)=IVG(I+K,J)
 204    CONTINUE
C
        RETURN
        END
      SUBROUTINE GEO2BL
      PARAMETER (IMX$=36,JMX$=50)
      COMMON /DATINC/ XLAT1,XLON1,DELLT,DELLN,IMXC,JMXC
      COMMON /IUV/   IU(IMX$,JMX$),IV(IMX$,JMX$),IP(IMX$,JMX$),
     1               IUG(IMX$,JMX$),IVG(IMX$,JMX$)
      INTEGER*2      IU,IV,IP,IUG,IVG
      COMMON /LANDS/ LANDS(IMX$,JMX$)
      CHARACTER*1    LANDS
      DATA CNSTN/4.769134/
      DATA XKS/5.E-1/,ROT/7.292/,XKSL/10.E-1/
C
       CNST=XKS/(2.*ROT)
       CNSTL=XKSL/(2.*ROT)
C        SIN AND COS OF CROSS-ISOBARIC ANGLE (SPECIFIED CONSTANT)
C       ANGC=COS(ANG*1.74532925E-2)
C       ANGS=SIN(ANG*1.74532925E-2)
C
C      REDUCE THE GEOSTROPHIC WIND TO SEA LEVEL. ACCOUNTING THE
C      BOUNDARY LAYER FRICTION BY ROTATING THE VECTOR TOWARD LOW
C      CENTER, WITH A SPECIFIED CROSS-ISOBARIC ANGLE AND REDUCING THE
C      SPEED BY A FACTOR.
C
C
CIGEO=1 MODE --- REPLACE INPUT WINDS IU,IV.
C
        DO 202 I=1,IMXC
        CORF=CNST/SIN((XLAT1-I+1)*1.74532925E-2)
        CORFL=CNSTL/SIN((XLAT1-I+1)*1.74532925E-2)
        DO 202 J=1,JMXC
C          ZIUG,ZIVG IN METER/SEC * 10.
       ZIUG=IUG(I,J)
       ZIVG=IVG(I,J)
       SPD=SQRT(ZIUG**2+ZIVG**2)/10.
       IF (LANDS(I,J).EQ.'0') THEN
        CORFZ=CORFL
        ELSE
        CORFZ=CORF
        ENDIF
C
       IF (SPD.GT.2.5) THEN
C
       ALPHA=SPD*CORFZ
C       COSINE AND SINE  ANG.
       ANGC=SQRT(-.5+.5*SQRT(1.+4.*ALPHA**2))/ALPHA
       ANGS=SQRT(1.-ANGC*ANGC)
C
       IU(I,J)=ANGC*(ANGC*ZIUG-ANGS*ZIVG)
       IV(I,J)=ANGC*(ANGS*ZIUG+ANGC*ZIVG)
       ENDIF
 202    CONTINUE
        RETURN
          END
      SUBROUTINE SETCP1(IGEO)
C        J. CHEN    MAY 1992 MDL   IBM 360/195
C
C        PURPOSE
C           THIS SUBROUTINE  READS IN SURFACE WIND COMPONENTS AND
C           PRESSURE FIELD, FOR THE FIRST TWO TIME PROJECTIONS.
C           IT CALLS IN THE INTERPOLATION PROGRAM TO SET UP FORCINGS
C           ON SLOSH GRIDS.
C
C        DATA SET USE
C
C        VARIABLES
C             ITIME = NUMBER OF TIME STEPS; INITIALIZED TO ZERO
C            NDLTHR = NUMBER OF TIME STEPS IN ONE HOUR
C             HHRAD = ACCUMULATIVE TIME STEPS ON THE HOUR
C             IBGNT = TABULATED SUBSCRIPT FOR STORM INITIALIZATION
C                     C21 INITIALLY SET IN 'INTVAL'
C              DELT = UNIT TIME STEP
C               C22 = DRCTN OF STORM MOTION, 1ST HR, SET IN 'INTVAL'
C            ITMADV = COUNTER FOR HOURLY TIME, 1ST HOUR
C
C        GENERAL COMMENTS
C
      PARAMETER (M$=169,N$=239,L$=1700,NCT$=1700,LC$=L$)
      COMMON /DUMB5/  ZB(M$,N$),ZBM(M$,N$)
      COMMON /DUMB7/  UB(M$,N$),VB(M$,N$),HB(M$,N$)
      COMMON /SCRTCH/ HSUB(M$,N$)
      COMMON /DUMB10/ HMX(M$,N$)
      COMMON /ELPDST/ ELPCT(M$),ELPDT(M$),ELPCL(M$),ELPDL(M$)
      COMMON /ELPDS2/ ELPCL2(M$),ELPDL2(M$),SINL2(N$),SINT2(N$)
      COMMON /PLRPST/ SINT(N$),COST(N$),SINL(N$),COSL(N$)
      COMMON /SETMIX/ IS(N$),MS(N$),MF(N$)
      COMMON /SETMX1/ IE(N$),ME(N$)
      COMMON /ITREE/  ITREE(M$,N$)
      COMMON /KTREE/  KTREE(L$)
      CHARACTER*1     ITREE,KTREE
      COMMON /KSPCON/ KSKP(M$,N$)
      CHARACTER*1     KSKP
      COMMON /DUMB18/ ISQR(L$),JSQR(L$),ISIDE(L$),HWEIR(L$),
     1                ZBMIN(L$)
      COMMON /ACT1D/  FLWSQR(L$),F1DACT(L$)
      CHARACTER*1     F1DACT
      COMMON /DUMCUT/ CUTL(NCT$),CUTLI(NCT$),CUTLE(NCT$)
C
      COMMON /CS1DA/  COS1D(LC$),SIN1D(LC$)
      COMMON /BANK2/  DELCUT(LC$),BANK(LC$,2)
      COMMON /LLXFLE/YLT(M$,N$),YLG(M$,N$)
      PARAMETER (IMX$=36,JMX$=50)
      COMMON /DATINC/ XLAT1,XLON1,DELLT,DELLN,IMXC,JMXC
      COMMON /UVP/ XIU(IMX$,JMX$,2),XIV(IMX$,JMX$,2),XIP(IMX$,JMX$,2),
     1               XIUG(IMX$,JMX$,2),XIVG(IMX$,JMX$,2)
      COMMON /IUV/   IU(IMX$,JMX$),IV(IMX$,JMX$),IP(IMX$,JMX$),
     1               IUG(IMX$,JMX$),IVG(IMX$,JMX$)
      INTEGER*2    IU,IV,IP,IUG,IVG
C
      COMMON /FFTH/   ITIME,MHALT
      COMMON /DUMB3/  IMXB,JMXB,IMXB1,JMXB1,IMXB2,JMXB2
      COMMON /DELTH/  NDLTHR,NDLTH2,DLTHR,DLTH2
      COMMON /SSTME/  ISTM,JHR,ITMADV,NHRAD,IBGNT,ITEND
      COMMON /WTIME/  NHRW,INCW,CLAT,CLON
      COMMON /STRMPS/ X(100),Y(100),PT(100),R(100),DIR(100),SP(100)
      COMMON /SPLN/   ALT(15),ALN(15),AX(15),AY(15),RL(15),ANGD(15),
     1                PZ(15),RT(15),XLAT(100),YLONG(100),RLNGTH
      COMMON /FRST/   X1(50),X12(50)
      COMMON /EGTH/   DELS,DELT,G,COR
      COMMON /BSN/    PHI,ALTO,ALNO,PHI1,ALT1,ALN1,ALT1C
      COMMON /THRD/   C7,C17,C19,C25
      COMMON /INTV/   INTV(24),KINPUT
      COMMON /DUMM/   DUM1,DUM2,DUM3
      COMMON /MTH/    MONTH(12),AMON(12)
      CHARACTER*40    FILNAM
      CHARACTER*80    DUM,DUM1,DUM2,DUM3
      CHARACTER*4     AMON,AMON1
C --------------------------------------------------------------------
      DATA RHOA/2.298E-3/,RHOW/4905./,RHOWG/29.89/
      REAL MTERFT
      DATA MTERFT/3.28084/
C-------------  TIME INTERVALS OF WIND INPUT ------
      INCW=6
      NHRW=0
      NHRAD=NDLTHR
      ITMADV=IBGNT+1
C ------------------------------------------------------
C          READ IN LAT/LONG OF THE SLOSH GRIDS IN USE
C      WRITE (*,*)' READ IN LAT/LONG OF THE GRIDS'
C      READ (*,'(A)') FILNAM
C      OPEN (25,FILE=FILNAM)
      READ (25) IDUM,JDUM
      READ (25) DUM,DUM
      READ (25) ((YLT(K,L),K=1,IMXB),L=1,JMXB)
      READ (25) ((YLG(K,L),K=1,IMXB),L=1,JMXB)
CC       CLOSE (25)
C------------  USE HOURLY STORM POSITIONS IN LAT/LONG -----
C
C      CLAT= XLAT(ITMADV+INCW-1)-XLAT(ITMADV-1)
C      CLON=YLONG(ITMADV+INCW-1)-YLONG(ITMADV-1)
       CALT=0.
       CLON=0.
C
C        INITIAL INPUT OF 2D WIND OF TWO TIME INTERVALS
C      JHR=61, IF IBGNT=JHR-48, KINPUT=3
      KINPUT=(IBGNT-1)/6+1
C      WRITE (*,*)'   JHR IBGNT KINPUT ', JHR,IBGNT,KINPUT
C      WRITE (*,*) 'ENTER FILENAME 1'
C      READ (*,'(A)') FILNAM
C      OPEN (33,FILE=FILNAM)
      KCNT=1
C     WRITE (*,*)'   KINPUT = ',KINPUT
 255  IF (KCNT.GT.KINPUT) GOTO 256
C     WRITE (*,*) ' KCNT',KCNT
C      JMXC MUST BE LESS OR EQUAL TO JMXC$.
      READ (33) IMXC,JMXC
      READ (33) DUM1,DUM2
C     WRITE (*,'(A)') DUM1
      READ (33) ((IU(I,J),I=1,IMXC),J=1,JMXC)
      READ (33) ((IV(I,J),I=1,IMXC),J=1,JMXC)
      READ (33) ((IP(I,J),I=1,IMXC),J=1,JMXC)
      KCNT=KCNT+1
      GOTO 255
 256  CONTINUE
C
       CALL GEOWND
       IF (IGEO.EQ.1) CALL GEO2BL
C
C
      DO 200 J=1,JMXC
      DO 200 I=1,IMXC
      XIU(I,J,1)=IU(I,J)*.1*MTERFT
      XIV(I,J,1)=IV(I,J)*.1*MTERFT
      XIP(I,J,1)=-(IP(I,J)*.1-12.)/RHOWG
C        DENSITY RATIO BETWEEN WATER AND AIR IS ASUMMED TO BE 1000.
      XIUG(I,J,1)=IUG(I,J)*.1*MTERFT
      XIVG(I,J,1)=IVG(I,J)*.1*MTERFT
 200  CONTINUE
C
C
      KINPUT=KINPUT+1
C     WRITE (*,*)'   KINPUT = ',KINPUT
      READ (33) IDUM,JMXC
      READ (33) DUM1,DUM2
C     WRITE (*,'(A)') DUM1
      READ (33) ((IU(I,J),I=1,IMXC),J=1,JMXC)
      READ (33) ((IV(I,J),I=1,IMXC),J=1,JMXC)
      READ (33) ((IP(I,J),I=1,IMXC),J=1,JMXC)
      CALL GEOWND
C
C       IF IGEO=1, FRICTIONALLY BALANCEED GEOSTRAPHIC WIND METHOD
C       IS USED.
         IF (IGEO.EQ.1) CALL GEO2BL
C
      DO 210 J=1,JMXC
      DO 210 I=1,IMXC
      XIU(I,J,2)=IU(I,J)*.1*MTERFT
      XIV(I,J,2)=IV(I,J)*.1*MTERFT
      XIP(I,J,2)=-(IP(I,J)*.1-12.)/RHOWG
      XIUG(I,J,2)=IUG(I,J)*.1*MTERFT
      XIVG(I,J,2)=IVG(I,J)*.1*MTERFT
 210  CONTINUE
C
C       LINEAR WEIGHTING IN TIME FOR TWO-PROJECTIONS, INCREMENTED BY
C       BY INCW HOURS. (12 OR 6 HRS)
C
      ZZ=NDLTHR*INCW
      WGT=(ITIME+1-NHRW)/ZZ
C      WRITE (*,*)'   ITIME NHRW', ITIME,NHRW,INCW,WGT,CLAT,CLON
C
      NBNDRY=0
      CALL WDPINT(WGT,CLAT,CLON,NBNDRY,IGEO)
C
      RETURN
       END
      SUBROUTINE STMVL1(IGEO)
C        J. CHEN                  MAY 1992 MDL   IBM 360/195
C
C        PURPOSE
C           THIS SUBROUTINE (STMVAL) TESTS FOR INCREMENTAL HOURS AND
C           SETS HOURLY TRACK AND METEOROLOGY VALUES. IT TESTS FOR 1-MB
C           CHANGE IN PRESSURE DROP, OR 1 MILE CHANGE IN STORM SIZE,
C           FOR EACH ADVANCE IN TIME. PREPARATORY VALUES ARE SET FOR
C           LATER ADVANCING IN TIME.
C
C        DATA SET USE
C           NONE
C
C        VARIABLES
C             ITIME = CUMULATIVE TIME COUNTER
C            NDLTHR = NUMBER OF TIME STEPS IN ONE HOUR
C            ITMADV = COUNTER FOR HOURLY TIME, SET IN 'SETCMP'
C
C        GENERAL COMMENTS
C
      PARAMETER (M$=169,N$=239,L$=1700,NCT$=1700,LC$=L$)
      COMMON /DUMB5/  ZB(M$,N$),ZBM(M$,N$)
      COMMON /DUMB7/  UB(M$,N$),VB(M$,N$),HB(M$,N$)
      COMMON /SCRTCH/ HSUB(M$,N$)
      COMMON /DUMB10/ HMX(M$,N$)
      COMMON /ELPDST/ ELPCT(M$),ELPDT(M$),ELPCL(M$),ELPDL(M$)
      COMMON /ELPDS2/ ELPCL2(M$),ELPDL2(M$),SINL2(N$),SINT2(N$)
      COMMON /PLRPST/ SINT(N$),COST(N$),SINL(N$),COSL(N$)
      COMMON /SETMIX/ IS(N$),MS(N$),MF(N$)
      COMMON /SETMX1/ IE(N$),ME(N$)
      COMMON /ITREE/  ITREE(M$,N$)
      COMMON /KTREE/  KTREE(L$)
      CHARACTER*1     ITREE,KTREE
      COMMON /KSPCON/ KSKP(M$,N$)
      CHARACTER*1     KSKP
      COMMON /DUMB18/ ISQR(L$),JSQR(L$),ISIDE(L$),HWEIR(L$),
     1                ZBMIN(L$)
      COMMON /ACT1D/  FLWSQR(L$),F1DACT(L$)
      CHARACTER*1     F1DACT
      COMMON /DUMCUT/ CUTL(NCT$),CUTLI(NCT$),CUTLE(NCT$)
C
      COMMON /CS1DA/  COS1D(LC$),SIN1D(LC$)
      COMMON /BANK2/  DELCUT(LC$),BANK(LC$,2)
      PARAMETER (IMX$=36,JMX$=50)
      COMMON /DATINC/ XLAT1,XLON1,DELLT,DELLN,IMXC,JMXC
      COMMON /UVP/ XIU(IMX$,JMX$,2),XIV(IMX$,JMX$,2),XIP(IMX$,JMX$,2),
     1               XIUG(IMX$,JMX$,2),XIVG(IMX$,JMX$,2)
      COMMON /IUV/   IU(IMX$,JMX$),IV(IMX$,JMX$),IP(IMX$,JMX$),
     1               IUG(IMX$,JMX$),IVG(IMX$,JMX$)
      INTEGER*2    IU,IV,IP,IUG,IVG
C
      COMMON /KEYS/   KEY1
      COMMON /FFTH/   ITIME,MHALT
      COMMON /BSN/    PHI,ALTO,ALNO,PHI1,ALT1,ALN1,ALT1C
      COMMON /DELTH/  NDLTHR,NDLTH2,DLTHR,DLTH2
      COMMON /EGTH/   DELS,DELT,G,COR
      COMMON /THRD/   C7,C17,C19,C25
      COMMON /STRMPS/ X(100),Y(100),PT(100),R(100),DIR(100),SP(100)
      COMMON /SPLN/   ALT(15),ALN(15),AX(15),AY(15),RL(15),ANGD(15),
     1                PZ(15),RT(15),XLAT(100),YLONG(100),RLNGTH
      COMMON /SSTME/  ISTM,JHR,ITMADV,NHRAD,IBGNT,ITEND
      COMMON /WTIME/  NHRW,INCW,CLAT,CLON
      COMMON /PRHSY/  IPRHR,JUMPR,IPRT
      COMMON /MAPTME/ ISTMAP,IADMAP
      COMMON /PRTSV/  IPRTSV,IPRTAD
      COMMON /SWTCH/  IOPERL(5)
      COMMON /DATE7/  IDATE(4)
      COMMON /INTV/   INTV(24),KINPUT
      COMMON /MTH/    MONTH(12),AMON(12)
      COMMON /DUMM/   DUM1,DUM2,DUM3
      COMMON /IDENT/  AIDENT(40),DACLOK(7)
      CHARACTER*4 AIDENT
      CHARACTER*40    FILNAM
      CHARACTER*80    DUM,DUM1,DUM2,DUM3
      CHARACTER*4     AMON,AMON1
C
C --------------------------------------------------------------------
      DATA RHOA/2.298E-3/,RHOW/4905./,RHOWG/29.89/
      REAL MTERFT
      DATA MTERFT/3.28084/
C
      ISKIP=0
      IF (ITIME.EQ.NHRAD) THEN
        NHRAD=NHRAD+NDLTHR
        ITMADV=ITMADV+1
        ENDIF
C        TEST FOR INCREMENTAL HOURLY CHANGES IN STORM TRACK
      IF (ITIME.GT.MHALT-2) RETURN
      IF(ITIME+1.NE.NHRW+NDLTHR*INCW) GO TO 100
C       WRITE (*,*)'  ITIME MHALT NHRW+',ITIME,MHALT,NHRW+NDLTHR*INCW
      NHRW=NHRW+NDLTHR*INCW
      NBNDRY=0
C
      DO 500 J=1,JMXC
      DO 500 I=1,IMXC
      XIU(I,J,1)=IU(I,J)*.1*MTERFT
      XIV(I,J,1)=IV(I,J)*.1*MTERFT
      XIP(I,J,1)=-(IP(I,J)*.1-12.)/RHOWG
      XIUG(I,J,1)=IUG(I,J)*.1*MTERFT
      XIVG(I,J,1)=IVG(I,J)*.1*MTERFT
 500  CONTINUE
C
C------------  LAT/LONG HOURLY STORM POSITIONS OF CENTER  -----
C              COMPSITE-FIELD INTERPOLATIONS ARE USED.
C      CLAT= XLAT(ITMADV+INCW-1)-XLAT(ITMADV-1)
C      CLON=YLONG(ITMADV+INCW-1)-YLONG(ITMADV-1)
C
C           LOCAL INTERPOLATIONS IN TIME IS USED.
       CLAT=0.
       CLON=0.
C
C      WRITE (*,*)' ITMADV-1',ITMADV-1,XLAT(ITMADV+INCW-1),
C     1     YLONG(ITMADV+INCW-1),CLAT,CLON
C
      KINPUT=KINPUT+1
C     WRITE (*,*)'   KINPUT = ',KINPUT
C     WRITE (*,*)' KINPUT INTV() ITIME',KINPUT,INTV(KINPUT),ITIME
C
C        KINPUT=11, FINISHED READING ANALYZED DATA
        IF (KINPUT.GE.11) THEN
C -------------------------------------------------------------------
C       WRITE (*,*) ' ENTER FILENAME 2 (UNIT=34)',KINPUT
C        WRITE (*,*) ' FORECAST WILL START 12 HOURS LATER.'
C        READ (*,'(A)') FILNAM
C        OPEN (34,FILE=FILNAM)
      READ (34) IDUM,JDUM
      READ (34) DUM1,DUM2
      READ (34) ((IU(I,J),I=1,IMXC),J=1,JMXC)
      READ (34) ((IV(I,J),I=1,IMXC),J=1,JMXC)
      READ (34) ((IP(I,J),I=1,IMXC),J=1,JMXC)
       IF (KINPUT.EQ.11) THEN
      DUM3=DUM1
      DO 344 K=1,20
      AIDENT(K)(1:4)=DUM3((K-1)*4+1:(K-1)*4+4)
 344  CONTINUE
C     WRITE (*,'(A)') DUM3
      DO 1001 L=1,4
      IDATE(L)=0
 1001  CONTINUE
       LB=18
      DO 1002 L=1,4
      LL=LB+L-1
      III=AMAX0(1,10**(4-L))
        IF (DUM3(LL:LL).EQ.' ') DUM3(LL:LL)='0'
C------------------------------------------------------------------
C        MAINFRAME   '0' = 240;   PC     '0' = 48.
C------------------------------------------------------------------
C         DECODE YEAR FROM ALPHAMERIC TO NUMERIC
      IDATE(1)=IDATE(1)+III*(ICHAR(DUM3(LL:LL))-240)
 1002  CONTINUE
C         DECODE DAY FROM ALPHAMERIC TO NUMERIC
      DO 1003 L=1,2
      LL=LB+5+L-1
      III=AMAX0(1,10*(2-L))
      IF (DUM3(LL:LL).EQ.' ') DUM3(LL:LL)='0'
      IDATE(3)=IDATE(3)+III*(ICHAR(DUM3(LL:LL))-240)
 1003  CONTINUE
C         DECODE HOUR FROM ALPHAMERIC TO NUMERIC
      DO 1004 L=1,2
      LL=LB+13+L-1
      III=AMAX0(1,10*(2-L))
      IF (DUM3(LL:LL).EQ.' ') DUM3(LL:LL)='0'
      IDATE(4)=IDATE(4)+III*(ICHAR(DUM3(LL:LL))-240)
 1004  CONTINUE
C       DECODE MONTH TO INTEGER
      AMON1=DUM3(LB+8:LB+10)
      DO 206 N=1,12
      IF (AMON1.EQ.AMON(N)) GOTO 207
 206  CONTINUE
 207  IDATE(2)=N
C      WRITE (*,*) IDATE
         ENDIF
C
C      FORECAST FIELDS, KINPUT > 11, ARE 6 HOURS APART.
         IF (KINPUT.GE.12) INCW=12.01/4.
C-----------------------------------------------------------------------
        ELSE
C ----------------------------------------------------------------------
      INCW=6.01
      READ (33) IDUM,JMXC
      DUM=DUM1
      READ (33) DUM1,DUM2
C     WRITE (*,'(A)') DUM1
      READ (33) ((IU(I,J),I=1,IMXC),J=1,JMXC)
      READ (33) ((IV(I,J),I=1,IMXC),J=1,JMXC)
      READ (33) ((IP(I,J),I=1,IMXC),J=1,JMXC)
        ENDIF
C-----------------------------------------------------------------------
C
C
      IF (DUM.EQ.DUM1) ISKIP=1
C
       CALL GEOWND
       IF (IGEO.EQ.1) CALL GEO2BL
C
      DO 200 J=1,JMXC
      DO 200 I=1,IMXC
      XIU(I,J,2)=IU(I,J)*.1*MTERFT
      XIV(I,J,2)=IV(I,J)*.1*MTERFT
      XIP(I,J,2)=-(IP(I,J)*.1-12.)/RHOWG
      XIUG(I,J,2)=IUG(I,J)*.1*MTERFT
      XIVG(I,J,2)=IVG(I,J)*.1*MTERFT
 200  CONTINUE
C
 100  CONTINUE
C
C
      IF (ISKIP.EQ.0) THEN
C       UPDATE THE WIND FILED ONCE AN HOUR FOR SAVING COMPUTATIONS.
      IF (MOD(ITIME+1,NDLTHR).EQ.0) THEN
      ZZ=NDLTHR*INCW
      WGT=(ITIME+1-NHRW)/ZZ
C     WRITE (*,*)'   ITIME NHRW', ITIME,NHRW,INCW,WGT,CLAT,CLON
C
      CALL WDPINT(WGT,CLAT,CLON,NBNDRY,IGEO)
      ENDIF
      ENDIF
C
      RETURN
       END
      SUBROUTINE WDPINT(WGT,CLAT,CLON,NBNDRY,IGEO)
C        J. CHEN               MAY 1992 MDL   IBM 360/195
C
C        PURPOSE
C           THIS SUBROUTINE INTERPOLATES DATA FIELDS FROM LAT/LON GRIDS
C           TO SLOSH GRIDS IN SPACE. FOR EACH TIME STEP BETWEEN HOURS
C           TWO DATA FIELDS ARE WEIGHTED THROUGH TRANSLATIONS AND
C           COMPSITING, BASING ON THE MOVEMENT VECTOR.
C
C        DATA SET USE
C           NONE
C
C        VARIABLES
C        YLT,YLG(,) = LAT/LON OF THE SLOSH GRID OF A BASIN
C         CLAT,CLON = TRANSLATION VECTOR FROM FIRST TO SECOND DATA FLD
C               WGT = THE WEIGHTING FACTOR FOR DATA FIELD AT TWO TIME
C                     LEVEL
C
      PARAMETER (M$=169,N$=239,L$=1700,NCT$=1700,LC$=L$)
      COMMON /DUMB5/  ZB(M$,N$),ZBM(M$,N$)
      COMMON /DUMB7/  UB(M$,N$),VB(M$,N$),HB(M$,N$)
      COMMON /SCRTCH/ HSUB(M$,N$)
      COMMON /DUMB10/ HMX(M$,N$)
      COMMON /ELPDST/ ELPCT(M$),ELPDT(M$),ELPCL(M$),ELPDL(M$)
      COMMON /ELPDS2/ ELPCL2(M$),ELPDL2(M$),SINL2(N$),SINT2(N$)
      COMMON /PLRPST/ SINT(N$),COST(N$),SINL(N$),COSL(N$)
      COMMON /SETMIX/ IS(N$),MS(N$),MF(N$)
      COMMON /SETMX1/ IE(N$),ME(N$)
      COMMON /ITREE/  ITREE(M$,N$)
      COMMON /KTREE/  KTREE(L$)
      CHARACTER*1     ITREE,KTREE
      COMMON /KSPCON/ KSKP(M$,N$)
      CHARACTER*1     KSKP
      COMMON /DUMB18/ ISQR(L$),JSQR(L$),ISIDE(L$),HWEIR(L$),
     1                ZBMIN(L$)
      COMMON /ACT1D/  FLWSQR(L$),F1DACT(L$)
      CHARACTER*1     F1DACT
      COMMON /DUMCUT/ CUTL(NCT$),CUTLI(NCT$),CUTLE(NCT$)
C
      COMMON /CS1DA/  COS1D(LC$),SIN1D(LC$)
      COMMON /BANK2/  DELCUT(LC$),BANK(LC$,2)
C
      PARAMETER (IMX$=36,JMX$=50)
      COMMON /DATINC/ XLAT1,XLON1,DELLT,DELLN,IMXC,JMXC
      COMMON /LLXFLE/YLT(M$,N$),YLG(M$,N$)
      COMMON /FORCE/ WX(M$,N$),WY(M$,N$),PP(M$,N$),PPX(M$,N$),PPY(M$,N$)
      COMMON /UVP/ XIU(IMX$,JMX$,2),XIV(IMX$,JMX$,2),XIP(IMX$,JMX$,2),
     1               XIUG(IMX$,JMX$,2),XIVG(IMX$,JMX$,2)
C
      COMMON /DUMB8/ IP(4),JP(4),IH(4),JH(4)
      COMMON /BSN/   PHI,ALTO,ALNO,PHI1,ALT1,ALN1,ALT1C
      COMMON /DUMB3/ IMXB,JMXB,IMXB1,JMXB1,IMXB2,JMXB2
      COMMON /DELTH/ NDLTHR,NDLTH2,DLTHR,DLTH2
      COMMON /BCPTS/ NBCPT,ISH(500),JSH(500)
C
C
      IF (NBNDRY.EQ.1) THEN
      DO 355 N=1,NBCPT
      DO 355 K=1,4
      I=ISH(N)+IP(K)
      J=JSH(N)+JP(K)
      XLT=YLT(I,J)
      XLN=YLG(I,J)
      CALL GRDINT(XLT,XLN,CLAT,CLON,WGT,OUT,XIP)
      PP(I,J)=OUT
 355  CONTINUE
      ELSE
C
      DO 350 I=1,IMXB
      DO 350 J=1,JMXB
      XLT=YLT(I,J)
      XLN=YLG(I,J)
      CALL GRDINT(XLT,XLN,CLAT,CLON,WGT,OUT,XIP)
      PP(I,J)=OUT
 350  CONTINUE
      ENDIF
C
      DO 150 J=1,JMXB
      DO 150 I=1,IMXB
      IF (ITREE(I,J).EQ.'6') GOTO 150
      XLT=YLT(I,J)
      XLN=YLG(I,J)
      CALL GRDINT(XLT,XLN,CLAT,CLON,WGT,OUT,XIU)
      IF (IGEO.EQ.1.AND.(ITREE(I,J).EQ.'3'.OR.ITREE(I,J).EQ.'1'))
     1       OUT=OUT*.8
      WX(I,J)=OUT
 150  CONTINUE
C
C
      DO 250 J=1,JMXB
      DO 250 I=1,IMXB
      IF (ITREE(I,J).EQ.'6') GOTO 250
      XLT=YLT(I,J)
      XLN=YLG(I,J)
      CALL GRDINT(XLT,XLN,CLAT,CLON,WGT,OUT,XIV)
      IF (IGEO.EQ.1.AND.(ITREE(I,J).EQ.'3'.OR.ITREE(I,J).EQ.'1'))
     1       OUT=OUT*.8
      WY(I,J)=OUT
 250  CONTINUE
C
C
      DO 450 J=1,JMXB
      DO 450 I=1,IMXB
      IF (ITREE(I,J).EQ.'6') GOTO 450
      XLT=YLT(I,J)
      CORR=2.*(7.292116E-5)*SIN(XLT*1.74532925E-2)
      XLN=YLG(I,J)
      CALL GRDINT(XLT,XLN,CLAT,CLON,WGT,OUT,XIUG)
      PPY(I,J)=-OUT*CORR/1000.
 450  CONTINUE
C
      DO 550 J=1,JMXB
      DO 550 I=1,IMXB
      IF (ITREE(I,J).EQ.'6') GOTO 550
      XLT=YLT(I,J)
      CORR=2.*(7.292116E-5)*SIN(XLT*1.74532925E-2)
      XLN=YLG(I,J)
      CALL GRDINT(XLT,XLN,CLAT,CLON,WGT,OUT,XIVG)
      PPX(I,J)=OUT*CORR/1000.
 550  CONTINUE
C
C        CHANGE NORTH/SOUTH COORDINATES INTO X/Y COORDINATE BY
C        ROTATING THE GRIDS
      CS=COS(PHI1)
      SS=SIN(PHI1)
C
      DO 600 J=1,JMXB
      DO 600 I=1,IMXB
      IF (ITREE(I,J).EQ.'6') GOTO 600
      ZA= CS*WX(I,J)-SS*WY(I,J)
      ZAA=SS*WX(I,J)+CS*WY(I,J)
      WX(I,J)=ZA
      WY(I,J)=ZAA
      ZA= CS*PPX(I,J)-SS*PPY(I,J)
      ZAA=SS*PPX(I,J)+CS*PPY(I,J)
      PPX(I,J)=ZA
      PPY(I,J)=ZAA
 600  CONTINUE
C
      RETURN
      END
       SUBROUTINE GRDINT(XLATIN,YLONIN,CLAT,CLON,WGT,FLDOUT,FLD)
C        J. CHEN                 MAY 1992 MDL   IBM 360/195
C
C        PURPOSE
C        THIS SUBROUTINE INTERPOLATES THE DATA BETWEEN TWO TIME LEVELS.
C        IT USES BI-LINEAR INTERPOLATIONS IN SPACE. IN TIME-INTERPOLA
C        TIONS, LINEAR WEIGHTING WITH THE COMPOSITE OF TWO FIELDS
C        FOLLOWING THE TRANSLATION VECTOR OF THE STORM CENTER AT TWO
C        TIME LEVELS. THE INPUT FIELDS ARE IN LAT/LONG GRID (1 DEGREE
C        RESOLUTION). TRANSLATES FORWARD FOR 1ST-FIELD, BACKWARD FOR
C        2ND-FIELD, THEN FORMS THE COMPOSITE WITH APPROPRIATE
C        WEIGHTINGS.
C
C        DATA SET USE
C           NONE
C
C        VARIABLES
C         IMXC,JMXC = DIMENSIONS OF INPUT FIELDS (1 DEG. INTERVAL)
C       XLAT1,XLON1 = LAT/LON OF NORTH/WEST CORNER OF THE DATA FIELD
C        FLD( , ,2) = DATA FIELDS OF TWO TIME LEVELS
C     XLATIN,YLONIN = A GRID LAT/LON ON WHICH INTERPOLATED VALUE IS TO
C                     BE ASSIGNED
C            FLDOUT = INTERPOLATED FIELD VALUE AT XLATON,YLONIN
C               WGT = WEIGHTING FACTOR FOR TWO TIME LEVEL, = 0 REDUCED
C                     TO 1ST TIME LEVEL, = 1 2ND TIME LEVEL
C         CLAT,CLON = TRANSLATION VECTOR BETWEEN TWO TIME LEVEL OF THE
C                     STORM CENTER (DEEP LOW PRESSURE CENTER)
C
      PARAMETER (IMX$=36,JMX$=50)
      COMMON /DATINC/ XLAT1,XLON1,DELLT,DELLN,IMXC,JMXC
       DIMENSION FLD(IMX$,JMX$,2)
       DIMENSION ZFLD(2)
C
C       XLATIN,YLONIN(-) = TARGET LAT/LON
C         FLD( ,) = SOURCE DATA FIELD
C       CLAT,CLON = TRANSLATION VECTOR FROM FIRST TO SECOND DATA FLD
C       WGT       = DISTANCE ALONG TRANSLATION VECTOR (NORMALIZE TO 1)
C                   ALSO, THE WEIGHTING FACTOR FOR TWO DATA FIELD
C
      DO 100 K=1,2
      IF (K.EQ.1) THEN
          FCT=WGT
          ELSE
          FCT=-1.+WGT
      ENDIF
      XLAT=XLATIN-CLAT*FCT
      YLON=YLONIN+CLON*FCT
      PA=XLAT1-XLAT+1.
      QA=YLON+XLON1+1.
      M=PA
      N=QA
C
      IF (M.LT.1.OR.N.LT.1.OR.M.GT.IMXC.OR.N.GT.JMXC) THEN
C      WRITE (*,*)'  LAT/LON',XLAT,YLON,XLATIN,YLONIN
C      WRITE (*,*) '  GRID POINT OUT OF BOUND ',PA,QA,M,N
       ELSE
      W1=PA-M
      W2=QA-N
      ZFLD(K)=(1.-W1)*(1.-W2)*FLD(M,N,K)+W1*(1.-W2)*FLD(M+1,N,K)
     1  +(1.-W1)*W2*FLD(M,N+1,K)+W1*W2*FLD(M+1,N+1,K)
C
      ENDIF
 100  CONTINUE
C
      FLDOUT=ZFLD(1)*(1.-WGT)+ZFLD(2)*WGT
C
       RETURN
       END
      SUBROUTINE FRCPNT(I,J,FXB,FYB,NCATG,CSHLTR,NPLS)
C        J. CHEN               MAY 1992 MDL   IBM 360/195
C
C        PURPOSE
C           THIS SUBROUTINE (FRCPNT) COMPUTES COMPONENTS OF SURFACE
C           DRIVING FORCES AT A MOMENTUM GRID POINT. COMPUTATIONS ARE
C        DATA SET USE
C           NONE
C
C        VARIABLES
C            CSHLTR = EXTINCTION COEF OF STRESS FOR WATER LESS THAN 1 FT
C            X12(4) = RADIUS OR MAX WIND IN FEET
C           X12(15) = SQUARE OF X12(4)
C    X12(9) X12(10) = COMPONENTS OF FORWARD SPEED, 1ST HOUR
C               C19 = FRIC CFFCNT (3X10-6)
C           FXX FYY = COMPONENTS OF STRESS
C           PXX PYY = COMPONENTS OF SURFACE PRESSURE GRADIENT FORCE
C         FXBP FYBP = COMPONENTS OF SURFACE DRIVING FORCES IN (X,Y)
C           FXB FYB = COMPONENTS OF SURFACE FORCES ON IMAGE PLANE
C        WX,WY( , ) = SURFACE WIND COMPONENTS IN X/Y COORDINATE
C      PPX,PPY( , ) = PRESSURE GRADIENT FORCING IN X/Y COORDINATE
C
C        GENERAL COMMENTS
C
      PARAMETER (M$=169,N$=239,L$=1700,NCT$=1700,LC$=L$)
      COMMON /DUMB5/  ZB(M$,N$),ZBM(M$,N$)
      COMMON /DUMB7/  UB(M$,N$),VB(M$,N$),HB(M$,N$)
      COMMON /SCRTCH/ HSUB(M$,N$)
      COMMON /DUMB10/ HMX(M$,N$)
      COMMON /ELPDST/ ELPCT(M$),ELPDT(M$),ELPCL(M$),ELPDL(M$)
      COMMON /ELPDS2/ ELPCL2(M$),ELPDL2(M$),SINL2(N$),SINT2(N$)
      COMMON /PLRPST/ SINT(N$),COST(N$),SINL(N$),COSL(N$)
      COMMON /SETMIX/ IS(N$),MS(N$),MF(N$)
      COMMON /SETMX1/ IE(N$),ME(N$)
      COMMON /ITREE/  ITREE(M$,N$)
      COMMON /KTREE/  KTREE(L$)
      CHARACTER*1     ITREE,KTREE
      COMMON /KSPCON/ KSKP(M$,N$)
      CHARACTER*1     KSKP
      COMMON /DUMB18/ ISQR(L$),JSQR(L$),ISIDE(L$),HWEIR(L$),
     1                ZBMIN(L$)
      COMMON /ACT1D/  FLWSQR(L$),F1DACT(L$)
      CHARACTER*1     F1DACT
      COMMON /DUMCUT/ CUTL(NCT$),CUTLI(NCT$),CUTLE(NCT$)
C
      COMMON /CS1DA/  COS1D(LC$),SIN1D(LC$)
      COMMON /BANK2/  DELCUT(LC$),BANK(LC$,2)
C
      COMMON /FORCE/ WX(M$,N$),WY(M$,N$),PP(M$,N$),PPX(M$,N$),PPY(M$,N$)
      COMMON /FRST/   X1(50),X12(50)
      COMMON /FRTH/   IC12,IC19,BT
      COMMON /THRD/   C7,C17,C19,C25
      COMMON /STRMSB/ C1,C2,C21,C22,AX,AY,PTENCY,RTENCY
      COMMON /SCND/   AR(900),AI(900),BR(900),BI(900),CR(900),CI(900)
      COMMON /MF/     DPH,HPD(4),ID
C
C        COMPUTING FORCE AT A GRID POINT
C       STRESS =DRAG COEFF* SPEED * WIND VECTOR
C       STRESS REDUCTION FOR WATER LESS THAN 1 FT, AND SPIN-UP PERIOD
      A=WX(I,J)
      B=WY(I,J)
C      FCTT=C19*SQRT(A*A+B*B)
C         DRAG COEFFICIENT DEPENDING ON WIND SPEEDS
      FCTW=SQRT(A*A+B*B)
C     WU'S FORMULA
C      C19X=1.0236E-6+5.366E-8*FCTW/3.28084
C
C      DRAG COEFF. TO BEST-FIT ET-SURGE HINDCASTS (JAN.-FEB., 1998)
C
      ZZZ=1.38+0.081*FCTW/3.28084
      ZZZ=MIN(3.,ZZZ)
      C19X=ZZZ*1.E-6
      FCTT=C19X*FCTW
C
      FCTT=FCTT*CSHLTR*BT
      FXX=FCTT*A
      FYY=FCTT*B
C
      FXB=CR(ID)*FXX-CI(ID)*FYY
      FYB=CR(ID)*FYY+CI(ID)*FXX
C        PRESSURE TERM IN THE FORCING FUNCTIONS
C
         IF (NPLS.NE.0) THEN
C
      PXX=-DPH*PPX(I,J)
      PYY=-DPH*PPY(I,J)
      FXB=FXB+BR(ID)*PXX-BI(ID)*PYY
      FYB=FYB+BR(ID)*PYY+BI(ID)*PXX
         ENDIF
C
      RETURN
       END
      SUBROUTINE BDRYHT
C        J. CHEN               MAY 1992 MDL   IBM 360/195
C
C        PURPOSE
C           THIS SUBROUTINE (BDRYHT) COMPUTES STATIC HEIGHTS ON DEEP
C           WATER BOUNDARIES ONLY, SEE DIAGRAM BELOW:
C
C        DATA SET USE
C           NONE
C
C        VARIABLES
C             AX AY = COMPS OF TOTAL STRM MOTION, ADVANCED IN 'STMVAL'
C             C1 C2 = INITIAL COMPS OF STORM, SET IN 'INTVAL'
C             IMXB1 = MAX I-SUBSCRIPT FOR HEIGHT POINTS
C            SEADTM = INITIAL HEIGHT OF THE SEA (NO STATIC HEIGHTS)
C         DELP(800) = STATIC HEIGHTS AT MILE INTERVALS FROM STORM CENTER
C        HB(  ,   ) = SURGE HEIGHTS
C        ITREE( , ) = A,   ON AT LEAST ONE BOUNDARY CORNER AS STATIC
C                     HEIGHT IS USED
C
C        GENERAL COMMENTS
C           THIS SUBROUTINE RESIDES IN OVERLAY 'CMPUTE'. IT IS CALLED
C           IN BY SUBROUTINE 'CMPUTE'.
C      THIS SUBROUTINE IS CALLED IN BY SUBROUTINE 'CONTTY'
C      WHICH RESIDES IN OVERLAY 'CMPUTE'
C
      PARAMETER (M$=169,N$=239,L$=1700,NCT$=1700,LC$=L$)
      COMMON /DUMB5/  ZB(M$,N$),ZBM(M$,N$)
      COMMON /DUMB7/  UB(M$,N$),VB(M$,N$),HB(M$,N$)
      COMMON /SCRTCH/ HSUB(M$,N$)
      COMMON /DUMB10/ HMX(M$,N$)
      COMMON /ELPDST/ ELPCT(M$),ELPDT(M$),ELPCL(M$),ELPDL(M$)
      COMMON /ELPDS2/ ELPCL2(M$),ELPDL2(M$),SINL2(N$),SINT2(N$)
      COMMON /PLRPST/ SINT(N$),COST(N$),SINL(N$),COSL(N$)
      COMMON /SETMIX/ IS(N$),MS(N$),MF(N$)
      COMMON /SETMX1/ IE(N$),ME(N$)
      COMMON /ITREE/  ITREE(M$,N$)
      COMMON /KTREE/  KTREE(L$)
      CHARACTER*1     ITREE,KTREE
      COMMON /KSPCON/ KSKP(M$,N$)
      CHARACTER*1     KSKP
      COMMON /DUMB18/ ISQR(L$),JSQR(L$),ISIDE(L$),HWEIR(L$),
     1                ZBMIN(L$)
      COMMON /ACT1D/  FLWSQR(L$),F1DACT(L$)
      CHARACTER*1     F1DACT
      COMMON /DUMCUT/ CUTL(NCT$),CUTLI(NCT$),CUTLE(NCT$)
C
      COMMON /CS1DA/  COS1D(LC$),SIN1D(LC$)
      COMMON /BANK2/  DELCUT(LC$),BANK(LC$,2)
C
      COMMON /FORCE/ WX(M$,N$),WY(M$,N$),PP(M$,N$),PPX(M$,N$),PPY(M$,N$)
      COMMON /BCPTS/  NBCPT,ISH(500),JSH(500)
      COMMON /DUMMY4/ S(800),C(800),P(800),DELP(800)
      COMMON /STRMSB/ C1,C2,C21,C22,AX,AY,PTENCY,RTENCY
      COMMON /SSTME/  ISTM,JHR,ITMADV,NHRAD,IBGNT,ITEND
      COMMON /DUMB3/  IMXB,JMXB,IMXB1,JMXB1,IMXB2,JMXB2
      COMMON /FRST/   X1(50),X12(50)
      COMMON /DATUM/  SEADTM,DTMLAK
C
C       STATIC HEIGHT BOUNDARY SQUARES DETERMINED IN SUBROUTINE DEPSFC
C
      DO 100 N=1,NBCPT
      I=ISH(N)
      J=JSH(N)
      HB(I,J)=.25*(PP(I,J)+PP(I+1,J)+PP(I,J+1)+PP(I+1,J+1))
     1          +SEADTM
C     HB(I,J)=MAX(0.,HB(I,J)-SEADTM)+SEADTM
  100 CONTINUE
C
      RETURN
       END
      SUBROUTINE INTLHT
C        J. CHEN   MAY 1992 MDL   IBM 360/195
C
C        PURPOSE
C           THIS SUBROUTINE (INTLHT) INITIALIZES THE SEA WITH QUIESCENT
C           WATER LEVELS PLUS STATIC HEIGHTS WITH THE STORM IN DEEP
C           WATER. STATIC HEIGHTS ARE NOT ADDED TO INLAND WATER BODIES.
C           INLAND WATER HEIGHTS ARE INITIALIZED IN SUBROUTINE 'DEPTHB'.
C
C        DATA SET USE
C           NONE
C
C        VARIABLES
C       IMXB1 JMXB1 = MAX SUBSCRIPTS FOR HT POINTS
C            MF(  ) = MOMENTUM P SUBSCRIPT, FOR END OF LAKE WINDS,
C             AX AY = COMPONENTS OF STORM TRAVERSE FOR 1ST UNIT OF TIME
C             C1 C2 = COMPONENTS OF STORM PSTN, TIME=0, SET IN 'INTVAL'
C COSL(  ) SINL(  ) = COMPONENTS OF AZIMUTH ON HEIGHT POINTS
C        GRIDRH(  ) = RADII OF CONCENTRIC CIRCLES THRU HT POINTS
C            SEADTM = INITTIAL, QUIESCENT, SEA SURFACE HEIGHT
C            DTMLAK = INITIAL WATER HEIGHT (INTERIOR WATER BODIES)
C         DELP(800) = STATIC HEIGTHS AT MILE INTERVALS FROM STRM CENTER
C        HB(  ,   ) = SURGE HEIGHTS
C      UB VB(  ,  ) = TRANSPORT COMPONENTS
C        HMX(  ,  ) = MAX SURGE HEIGHTS
C         ZB(  ,  ) = DEPTH VALUES ON TW0-DIMENSIONAL STAIR STEPS
C
C        GENERAL COMMENTS
C           THIS SUBROUTINE RESIDES IN OVERLAY 'CMPUTE'. IT IS
C           CALLED IN ONLY ONCE BY SUBROUTINE 'SETCMP'.
C
      PARAMETER (M$=169,N$=239,L$=1700,NCT$=1700,LC$=L$)
      COMMON /DUMB5/  ZB(M$,N$),ZBM(M$,N$)
      COMMON /DUMB7/  UB(M$,N$),VB(M$,N$),HB(M$,N$)
      COMMON /SCRTCH/ HSUB(M$,N$)
      COMMON /DUMB10/ HMX(M$,N$)
      COMMON /ELPDST/ ELPCT(M$),ELPDT(M$),ELPCL(M$),ELPDL(M$)
      COMMON /ELPDS2/ ELPCL2(M$),ELPDL2(M$),SINL2(N$),SINT2(N$)
      COMMON /PLRPST/ SINT(N$),COST(N$),SINL(N$),COSL(N$)
      COMMON /SETMIX/ IS(N$),MS(N$),MF(N$)
      COMMON /SETMX1/ IE(N$),ME(N$)
      COMMON /ITREE/  ITREE(M$,N$)
      COMMON /KTREE/  KTREE(L$)
      CHARACTER*1     ITREE,KTREE
      COMMON /KSPCON/ KSKP(M$,N$)
      CHARACTER*1     KSKP
      COMMON /DUMB18/ ISQR(L$),JSQR(L$),ISIDE(L$),HWEIR(L$),
     1                ZBMIN(L$)
      COMMON /ACT1D/  FLWSQR(L$),F1DACT(L$)
      CHARACTER*1     F1DACT
      COMMON /DUMCUT/ CUTL(NCT$),CUTLI(NCT$),CUTLE(NCT$)
C
      COMMON /CS1DA/  COS1D(LC$),SIN1D(LC$)
      COMMON /BANK2/  DELCUT(LC$),BANK(LC$,2)
      COMMON /FORCE/ WX(M$,N$),WY(M$,N$),PP(M$,N$),PPX(M$,N$),PPY(M$,N$)
C
      COMMON /DUMB3/  IMXB,JMXB,IMXB1,JMXB1,IMXB2,JMXB2
      COMMON /STRMSB/ C1,C2,C21,C22,AX,AY,PTENCY,RTENCY
      COMMON /SSTME/  ISTM,JHR,ITMADV,NHRAD,IBGNT,ITEND
      COMMON /FLWCPT/ NSQRS,NSQRW,NSQRWC,NPSS,NCUT
      COMMON /DATUM/  SEADTM,DTMLAK
      COMMON /DATUM1/ DTMCHN
      COMMON /FRST/   X1(50),X12(50)
      COMMON /DUMMY4/ S(800),C(800),P(800),DELP(800)
      COMMON /SWTCH/  IOPERL(5)
      COMMON /XOKE/   XOKE
      CHARACTER*1     XOKE
C
      COMMON /GPRT1/  DOLLAR,EBSN
      CHARACTER*2     DOLLAR
      CHARACTER*1     EBSN
C
C
C        INITIALIZE WATER HEIGHTS ACCORDING TO LAKE DATUM
      DO 90 J=1,JMXB1
      DO 90 I=1,IMXB1
      IF (ITREE(I,J).EQ.'2'.OR.ITREE(I,J).EQ.'5') GOTO 95
      IF (ITREE(I+1,J).EQ.'2'.OR.ITREE(I+1,J).EQ.'5') GOTO 95
      IF (ITREE(I,J+1).EQ.'2'.OR.ITREE(I,J+1).EQ.'5') GOTO 95
      IF (ITREE(I+1,J+1).EQ.'2'.OR.ITREE(I+1,J+1).EQ.'5') GOTO 95
      ZZ=.25*(PP(I,J)+PP(I+1,J)+PP(I,J+1)+PP(I+1,J+1)) +DTMLAK
C     ZZ=MAX(0.,ZZ-DTMLAK)+DTMLAK
      HB(I,J)=AMAX1(-ZB(I,J),AMAX1(ZZ,-ZB(I,J)))
      GOTO 90
C        STATIC HEIGHTS ON OCEAN,SEA OR GULF; NOT ON INLAND WATER BODIES
 95   CONTINUE
      HB(I,J)=.25*(PP(I,J)+PP(I+1,J)+PP(I,J+1)+PP(I+1,J+1))+SEADTM
C     HB(I,J)=MAX(0.,HB(I,J)-SEADTM)+SEADTM
      HB(I,J)=AMAX1(HB(I,J),-ZB(I,J))
 90   CONTINUE
 100  CONTINUE
 110  CONTINUE
C
C        INITIALIZE U/V/HMX
C
      DO 320 I=1,IMXB1
 320  HB(I,JMXB)=HB(I,1)
      DO 200 J=1,JMXB
      DO 200 I=1,IMXB
      UB(I,J)=0.
  200 VB(I,J)=0.
      DO 210 J=1,JMXB1
      DO 210 I=1,IMXB1
      HMX(I,J)=HB(I,J)
 210  CONTINUE
C
C        INITLIZE FLWSQR TO BE ZERO.
C
      DO 220 L=1,NSQRWC
 220  FLWSQR(L)=0.
C
      RETURN
       END
C
C-----------------------------------------------------------------------
       SUBROUTINE CONTTY
C        J. CHEN     SEPTEMBER 1992  MDL   IBM 360/195
C
C        PURPOSE
C           THIS SUBROUTINE (CONTTY) IS AN ALGORITHM FOR THE CONTINUITY
C           EQUATION, FOR AN IMAGE PLANE; POLAR COORDINATES ARE
C           TRANSFORMED ONTO THE IMAGE PLANE. THE SYSTEM AXIS ARE
C           (P,Q), BUT THE SUBROUTINE USES (I,J) FOR INTEGER SBSCRPTNG
C           TWO LEVELS IN TIME (FORWARD IN CONTTY PLUS BACKWARD IN MOMN)
C           ARE USED. THE GRID (SEE DIAGRAM BELOW) IS STAGGERED IN SPACE
C           A 5-POINT SCHEME (4 CORNERS OF A SQUARE) ARE USED FOR
C           THE MOMENTUM SPACE GRADIENT.
C
C
C        DATA SET USE
C
C        VARIABLES
C       IMXB1 JMXB1 = MAX SUBSCRIPTS FOR SURGE SQUARES
C        HB(  ,   ) = SURGE HEIGHTS
C       HSUB(  ,  ) = STORAGE FOR A SECOND COPY OF HB
C            IS(  ) = START LOOP, EXCLUDE HI TERRAIN, SET IN 'CRDRD2'
C            IB IFN = BEGINNING/ENDING A 'DO LOOP' IN I-DIRECTION
C             ISKIP = 1,FOR J=1,CHECK ZBM AT 2 BNDRY CORNERS FOR SKIPPIN
C                       STAATIC CONDITION.
C                     2,AND I=IMXB1 FOR BOTTOM BOUNDARY
C                     3,FOR J=JMXB1 FOR RIGHT BOUNDARY
C              JUMP = CONDITIONAL INTEGER FOR STANDARD LOOP,
C                     OR RETRIEVAL OF PREVIOUSLY COMPUTED SURGES,
C                     SEE DIAGRAM ABOVE FOR AFFECTED SQUARES
C              LEAP = CONDITIONAL 'GO TO' TO DRY AN (I,J) SQUARE,
C                     SEE DIAGRAM ABOVE FOR AN (I,J) SQUARE
C         IFST ISND = BEGIN/END I SUBSCRIPT ON UPDATED (J-1) SQUARES
C             IPASS = CONDITIONAL INTEGER, TO SET A SQUARE DRY
C         ZB(  ,  ) = DEPTHS AT CENTER OF SURGE SQUARES
C       IP(4) JP(4) = SHIFT SUBSCRIPTS TO 4 MOMENTUM POINTS
C     UB VB(  ,   ) = COMPONENTS OF MOMENTUM FIELDS
C              HPST = DUMMY VALUE FOR SUEGE VALUEC
C           DUX DVY = COMPONENTS OF MOMENTUM GRADIENT
C               FCT = FACTOR TO MULTIPLY MOMENTUM GRADIENT
C              CHNG = SURGE PLUS DEPTHS (TO TEST FOR DRY SQUARE)
C              DENM = CHANGE IN SURGE, ONE TIME ITERATION
C              ANUM = CHANGE IN SURGE HT, TO TO DRY A SQUARE
C              RFCT = RATIO, <1, OF ANUM TO DENM, REDUCTION FACTOR
C              ISET = I-SUBSCRIPT FOR DRIED SQUARE
C
C        GENERAL COMMENTS
C           THIS SUBROUTINE RESIDES IN OVERLAY 'CMPUTE'. IT IS CALLED
C           IN BY SUBROUTINE 'CMPUTE'
C           THE SYSTEM ROUTINE 'ERRSET' IS USED TO ISOLATE VARIABLES
C           IF AN UNDER/OVERFLOW OCCURS; VARIABLE PRINT OUT IS
C           AIDED VIA 'HELP' AND THE 'COMMON/CHELP/'.
C
C           THE SUBSCRIPT FOR COMPUTATIONS OF SURGES ARE OFFSET (.5,.5)
C           FROM THE MOMENTUM POINTS. SEE BELOW:
C
C             (I,J)       (I,J+1)
C              .         .
C                                  .=MOMENTUM POINTS (UB,VB)
C                                  +=SURGE    POINTS (HB)
C                 (I,J)
C                   +
C                                  .=BARRIER  POINTS (ZBM)
C                                  +=DEPTH    POINTS (ZB)
C
C
C             (I+1,J)   (I+1,J+1)
C
      PARAMETER (M$=169,N$=239,L$=1700,NCT$=1700,LC$=L$)
      COMMON /DUMB5/  ZB(M$,N$),ZBM(M$,N$)
      COMMON /DUMB7/  UB(M$,N$),VB(M$,N$),HB(M$,N$)
      COMMON /SCRTCH/ HSUB(M$,N$)
      COMMON /DUMB10/ HMX(M$,N$)
      COMMON /ELPDST/ ELPCT(M$),ELPDT(M$),ELPCL(M$),ELPDL(M$)
      COMMON /ELPDS2/ ELPCL2(M$),ELPDL2(M$),SINL2(N$),SINT2(N$)
      COMMON /PLRPST/ SINT(N$),COST(N$),SINL(N$),COSL(N$)
      COMMON /SETMIX/ IS(N$),MS(N$),MF(N$)
      COMMON /SETMX1/ IE(N$),ME(N$)
      COMMON /ITREE/  ITREE(M$,N$)
      COMMON /KTREE/  KTREE(L$)
      CHARACTER*1     ITREE,KTREE
      COMMON /KSPCON/ KSKP(M$,N$)
      CHARACTER*1     KSKP
      COMMON /DUMB18/ ISQR(L$),JSQR(L$),ISIDE(L$),HWEIR(L$),
     1                ZBMIN(L$)
      COMMON /ACT1D/  FLWSQR(L$),F1DACT(L$)
      CHARACTER*1     F1DACT
      COMMON /DUMCUT/ CUTL(NCT$),CUTLI(NCT$),CUTLE(NCT$)
C
      COMMON /CS1DA/  COS1D(LC$),SIN1D(LC$)
      COMMON /BANK2/  DELCUT(LC$),BANK(LC$,2)
C
      COMMON /DUMB3/  IMXB,JMXB,IMXB1,JMXB1,IMXB2,JMXB2
      COMMON /DUMB8/  IP(4),JP(4),IH(4),JH(4)
      COMMON /DUM88/  IIH(4),JJH(4),IHH(4),JHH(4)
      COMMON /DUMB4/  X1B(10)
      COMMON /CHELP/  ANUM,DENM,RFCT,I,J,HBJ,HPST,UP(4),VP(4)
C
C
C         SET PARAMETERS FOR J=1
      M=0
C
      DO 320 J=1,JMXB1
      IB=IS(J)
      IFN=IE(J)
C
      DO 320 I=IB,IFN
C
      IF (KSKP(I,J).EQ.'0') GOTO 320
C
      DUX=-UB(I,J)+UB(I+1,J)+UB(I+1,J+1)-UB(I,J+1)
C
C        TEST FOR DEPTH > 30 FT
      IF(ZB(I,J).GT.30.) GO TO 210
C
C        TEST FOR DRY SQUARES, IF U'S AT 4 CORNERS ARE 0., OR DUX=0.
C        IN EQUIVALENCE.
      IF (DUX.EQ.0.) GO TO 320
C
C        CONTINUITY EQUATION.
 210  HPST=HB(I,J)
      Z=ELPDL2(I)+(ELPCL2(I)-ELPDL2(I))*SINL2(J)
      FCT=X1B(5)/Z
      DVY=-VB(I,J)-VB(I+1,J)+VB(I+1,J+1)+VB(I,J+1)
      HB(I,J)=HPST-FCT*(DUX+DVY)
C
C**********END OF LOOP FOR REGULAR RETURN*******************************
C          ADJUSTMENT BELOW TERRAIN HAS BEEN REMOVED.
 320  CONTINUE
C
C***********  RETURN WITH A NEW J-LINE  **********************
C
      DO 360 I=1,IMXB1
 360  HB(I,JMXB)=HB(I,1)
C
      RETURN
       END
       SUBROUTINE MOMNTM
C        J. CHEN   SEPTEMBER 1980 MDL   IBM 360/195
C
C        PURPOSE
C           THIS SUBROUTINE (MOMNTM) COMPUTES FUTURE MOMENTUM COMPONENTS
C           ON ALL INTERIOR GRID POINTS, AFTER FIRST COMPUTING FUTURE
C           SURGE BY SUBROUTINE 'CONTTY'. IT CHECKS IF MOMNTM CAN EXIST
C           AT A GRID POINT; E.G.,  OVERTOPPING THE GRID POINT .
C           FUTURE SURGES UPDATED FROM SUB. 'CONTTY' ARE USED TO COMPUTE
C           HEIGHT GRADIENTS AND TO CHECK FOR ELIGIBLE MOMENTUM POINTS--
C           BACKWARD SCHEME---. CHECKS ARE MADE IF MOMENTUM IS TO BE
C           COMPUTED AT OPEN BOUNDARY POINTS.
C
C        DATA SET USE
C           NONE
C
C        VARIABLES
C             NCATG = 1 FOR LAKE WINDS, 2 FOR OCEAN WINDS
C              NPLS = 0 REMOVE STATIC PRESSURE TERM IN FORCING(INT.B.C.)
C                   = 1  OTHERWISE
C             JMXB1 = MAX NO. OF J-POINTS, LESS A BOUNDARY POINT
C              JUMP = 1, FULLY FLOODED; 2,PARTIALLY FLOODED (NO STRESS)
C            MS(  ) = BEGIN I-SUBSCRIPT ON A J-LINE
C       ITREE( ,  ) = INDICATOR ARRAY FOR LAKE/OCEAN, TREE/NO TREE
C                     A,B,C,D FOR BNOUDARY TYPES
C        ZBM(  ,  ) = MAXIMUM BARRIER HT AT A SURGE POINT
C       IH(4) JH(4) = SHIFT I/J-SUBSCRIPTS TO 4 SURRNDNG SRGE PTS
C       IP(4) JP(4) = SHIFTS I/J-SUBPTS FROM SURGE PTS TO 4 MOM. CORNERS
C        HB(  ,   ) = SURGE HEIGHTS
C         ZB(  ,  ) = DEPTHS
C             HF(4) = STORAGE FOR 4 SURRNDNG SURGE PTS,
C            HPD(4) = STORAGE FOR 4 SURRNDNG TOTAL DEPTH POINTS
C     UB VB(  ,   ) = COMPONENTS OF TRANSPORT
C           DHX DHY = COMPONENTS OF SURFACE GRADIENT
C            ID DPH = MEAN TOTAL DEPTH
C   AI(300) AR(300) = FRICTION COEFFICIENTS ON TRANSPORTS
C           FXB FYB = COMPONENTS OF DRIVING FORCES
C           X1B(4)  = DEL*CORIOLIS PARM
C            CSHLTR = EXTINCTION COEF OF STRESS FOR WATER LESS THAN 1 FT
C             ISKIP = 1, COMPUTE MOMENTUM ON J=1    BOUNDARY
C                     2, COMPUTE MOMENTUM ON I=2 BOUNDARY
C                     3, COMPUTE MOMENTUM ON J=JMXB BOUNDARY
C                     4, COMPUTE MOMENTUM ON I=IMXB BOUNDARY
C             IJUMP = 1, SHALLOW WATER BOUNDARY COMPUTATIONS
C                     2, INTER DEPTH   BOUNDARY COMPUTATIONS
C    ISBS JSBS(4) = SHIFT SUBSCRIPTS FROM INTERIOR TO BOUNDARY PT
C         ISHF(3,4) = SUBSCIPTS SHIFTS FROM INTERIOR TO 3 CORNER
C                     MOM-POINTS FOR ALL 4 CORNERS
C
C        GENERAL COMMENTS
C           THIS SUBROUTINE RESIDES IN OVERLAY 'CMPUTE'. IT IS CALLED
C           IN BY SUBROUTINE 'CMPUTE'. THE PRESENT SUBROUTINE CALLS
C           IN 2 SUBROUTINES, 'FRCPNT' AND 'MNTMBD'.
C
C
      PARAMETER (M$=169,N$=239,L$=1700,NCT$=1700,LC$=L$)
      COMMON /LLXFLE/ YLT(M$,N$),YLG(M$,N$)
      COMMON /DUMB5/  ZB(M$,N$),ZBM(M$,N$)
      COMMON /DUMB7/  UB(M$,N$),VB(M$,N$),HB(M$,N$)
      COMMON /SCRTCH/ HSUB(M$,N$)
      COMMON /DUMB10/ HMX(M$,N$)
      COMMON /ELPDST/ ELPCT(M$),ELPDT(M$),ELPCL(M$),ELPDL(M$)
      COMMON /ELPDS2/ ELPCL2(M$),ELPDL2(M$),SINL2(N$),SINT2(N$)
      COMMON /PLRPST/ SINT(N$),COST(N$),SINL(N$),COSL(N$)
      COMMON /SETMIX/ IS(N$),MS(N$),MF(N$)
      COMMON /SETMX1/ IE(N$),ME(N$)
      COMMON /ITREE/  ITREE(M$,N$)
      COMMON /KTREE/  KTREE(L$)
      CHARACTER*1     ITREE,KTREE
      COMMON /KSPCON/ KSKP(M$,N$)
      CHARACTER*1     KSKP
      COMMON /DUMB18/ ISQR(L$),JSQR(L$),ISIDE(L$),HWEIR(L$),
     1                ZBMIN(L$)
      COMMON /ACT1D/  FLWSQR(L$),F1DACT(L$)
      CHARACTER*1     F1DACT
      COMMON /DUMCUT/ CUTL(NCT$),CUTLI(NCT$),CUTLE(NCT$)
C
      COMMON /CS1DA/  COS1D(LC$),SIN1D(LC$)
      COMMON /BANK2/  DELCUT(LC$),BANK(LC$,2)
C
      COMMON /MF/     DPH,HPD(4),ID
      COMMON /DUMB3/  IMXB,JMXB,IMXB1,JMXB1,IMXB2,JMXB2
      COMMON /EGTH/   DELS,DELT,G,COR
      COMMON /SCND/   AR(900),AI(900),BR(900),BI(900),CR(900),CI(900)
      COMMON /DUMB8/  IP(4),JP(4),IH(4),JH(4)
      COMMON /DUM88/  IIH(4),JJH(4),IHH(4),JHH(4)
      COMMON /DUMB4/  X1B(10)
      COMMON /DUMB4I/ CORL(90)
      COMMON /GPRT1/  DOLLAR,EBSN
      DIMENSION     SLPAI(900),HFX(4),ISBS(4),JSBS(4),ISHF(3,4)
      CHARACTER*2   DOLLAR
      CHARACTER*1   EBSN
      DATA ISBS/0,-1,0,1/,  JSBS/-1,0,1,0/
      DATA SHGT/50./
      DATA ISHF/1,2,4, 1,2,3, 1,3,4, 2,3,4/
C
C
      DO 10 I=1,899
 10   SLPAI(I)=AI(I+1)-AI(I)
      SLPAI(900)=SLPAI(899)
C
      JEND=JMXB1
      IF (DOLLAR.EQ.'$') JEND=JMXB
C
      DO 570 J=2,JEND
C
C         ISKIP SET 3 RANGES FOR J
C
      ISKIP=2
C        FOR CILCULAR ISLAND, NO SIDE-BOUNDARIES (1 OR 3) TO BE TESTED.
C        FOR J=1, U/V  REPEAT THOSE AT J=JMXB (JEND).
      IF (DOLLAR.NE.'$') THEN
        IF(J.EQ.2) ISKIP=1
        IF(J.EQ.JEND) ISKIP=3
        ENDIF
C
      IST=MS(J)
      IFN=ME(J)
C        SKIP A J-COLUMN IF MS(J) IS SET TO IMXB OR GREATER.
      IF ((IFN-IST).LT.0) GO TO 570
C
      DO 560 I=IST,IFN
      NCATG=2
      NPLS=1
      CSHLTR=1.
      JUMP=1
C        NO TERRAIN HIGHER THAN 35 FT CAN BE FLOODED.
      IF (ITREE(I,J).EQ.'4'.OR.ITREE(I,J).EQ.'6') GO TO 560
      IF (ZBM(I,J).GT.0.) GOTO 560
C
C       TEMPORARILY STORE PRESENT U/V TO BE USED FOR BOUNDARY POINTS
      UBJ=UB(I,J)
      VBJ=VB(I,J)
C
C
C        ALL SURROUNDING SQUARES ARE FLOODED
      JUMP=1
      DHX=-HB(I-1,J-1)+HB(I,J-1)-HB(I-1,J)+HB(I,J)
      DHY=-HB(I-1,J-1)-HB(I,J-1)+HB(I-1,J)+HB(I,J)
      DPH=.25*(HB(I-1,J-1)+HB(I,J-1)+HB(I-1,J)+HB(I,J)+
     1         ZB(I-1,J-1)+ZB(I,J-1)+ZB(I-1,J)+ZB(I,J))
      DPH=MAX(1.,DPH)
      DPHZ=AMIN1(899.,DPH)
      ID=DPHZ
      DPHI=ID
      Z=DPHZ-DPHI
      AINTRP=AI(ID)+SLPAI(ID)*Z
C
 450  CONTINUE
      AII=1.+X1B(4)*AINTRP
C
C      CORIOLIS PARAMETER DEPENDS ON THE THE LATITUDE OF THE GRID.
      LLAT=YLT(I,J)
      ZCOR=CORL(LLAT+1)+(YLT(I,J)-LLAT)*(CORL(LLAT+2)-CORL(LLAT+1))
      ARJ=AR(ID)*X1B(4)*ZCOR
      FCT=X1B(6)*DPH
C
C        LOCATE GRID (I,J) ON PHYSICAT PLANE  Z=(XR,YR)=ZETA
      XR=ELPCT(I)*COST(J)
      YR=ELPDT(I)*SINT(J)
C
C      CALL FRCPNT(XR,YR,FXBP,FYBP,NCATG,CSHLTR,NPLS)
C
C     TWO DIMENSIONAL FORCING FUNCTIONS ARE USED. INTERPOLATIONS FROM
C     SOURCE DATA TO MODEL GRIDS LOCATIONS ARE NEEDED.
C
      CALL FRCPNT(I,J,FXBP,FYBP,NCATG,CSHLTR,NPLS)
C        CHANGE FORCING FROM (X,Y) SYSTEM BACK TO IMAGE PLANE,
C        CONJ(DZ/DZETA)*(FXBP,FYBP), WHERE DZ/DZETA=(XR,YR)
      XR=ELPDT(I)*COST(J)
      YR=ELPCT(I)*SINT(J)
      FXB=(XR*FXBP+YR*FYBP)
      FYB=(XR*FYBP-YR*FXBP)
C
      XTTX=DELT*FXB
      XTTY=DELT*FYB
C
C      AII=1.-3.E-3*DELT*SQRT((UBJ**2+VBJ**2)/(XR**2+YR**2))
C     1                   /DPH**2
C
      AIU=AII*UB(I,J)
      AIV=AII*VB(I,J)
C        COMPUTE MOMENTUM ON INTERIOR POINT, SIELECKI'S SCHEME FOR
C        CORIOLIS TERM IS USED.
      UB(I,J)=AIU-FCT*(BR(ID)*DHX-BI(ID)*DHY)+ARJ*VB(I,J)+XTTX
      VB(I,J)=AIV-FCT*(BR(ID)*DHY+BI(ID)*DHX)-ARJ*UB(I,J)+XTTY
C
C
C ------------------------------------------------------------------
C        CHECK IF MOMENTUM COMPUTATIONS REQUIRED ON BOUNDARY
C        ISKIP=1, LEFT BOUNDARY(J=1)
C             =2, AND I=2, NEXT TO TOP BOUNDARY
C             =3, RIGHT BOUNDARY
C             =4, AND I=IMXB1,NEXT TO BOTTOM BOUNDARY
C        IJUMP=1, INTERMEDIATE WATER; =2, SHALLOW WATER
C
 500  IJUMP=1
      GO TO (530,520,540),ISKIP
C
C        CHECK FOR 4 CORNER POINTS AND SET ICOR.
C        ICOR=1,LEFT-TOP; 2,LEFT-BOTTOM;3,RIGHT-TOP;4,RIGHT-BOTTOM
 530  IF(I.NE.ME(J)) GO TO 512
      ICOR=2
 512  IF(I.NE.MS(J)) GO TO 510
      ICOR=1
      GO TO 600
 540  IF(I.NE.ME(J)) GO TO 514
      ICOR=4
 514  IF(I.NE.MS(J)) GO TO 510
      ICOR=3
      GO TO 600
C
C        TEST FOR BOTTOM BOUNDARY
 520  IF(I.EQ.MS(J)) GO TO 510
      IF(I.EQ.ME(J)) GO TO 516
      GO TO 560
 516  ISKIP=4
 510  CONTINUE
      IA=ISBS(ISKIP)
      JA=JSBS(ISKIP)
      ISB=I+IA
      JSB=J+JA
C        IF AT (I,J),TRANSPORTS ARE ZERO,SET TRANSPORTS AT BOUNDY ZEROS.
      IF(UB(I,J).EQ.0.) THEN
      UB(ISB,JSB)=0.
      VB(ISB,JSB)=0.
        GOTO 560
        ENDIF
C        PRE-SET ITREEALONG BOUNDARIES; =4  FOR LAND (SKIP COMP.)
C                                       =9  FOR SHALLOW WATER
C                                       =8  FOR INTERMEDIATE DEPTHS
C                                       =6  FOR DEEP WATER (SKIP COMP.)
      IF(ITREE(ISB,JSB).EQ.'6'.OR.ITREE(ISB,JSB).EQ.'4') GO TO 560
      IF(ITREE(ISB,JSB).EQ.'8' ) IJUMP=2
C      WRITE (*,*) '  MNTMBD ISB,JSB',ISB,JSB,IJUMP
      CALL MNTMBD(JUMP,IJUMP,ISB,JSB,ARJ,AII,I,J,UBJ,VBJ,
     1   FXB,FYB,NCATG,CSHLTR,NPLS)
      GO TO 560
C        THREE MOMENTUM POINTS MUST BE CHECKED FOR COMPUTATIONS FOR
C
 600  DO 610 KK=1,3
      K=ISHF(KK,ICOR)
      ISB=I+IH(ICOR)+IP(K)
      JSB=J+JH(ICOR)+JP(K)
C        IF AT (I,J),TRANSPORTS ARE ZERO,SET TRANSPORTS AT BOUNDY ZEROS.
      IF(UB(I,J).EQ.0.) GOTO 630
C
      IF(ITREE(ISB,JSB).EQ.'6'.OR.ITREE(ISB,JSB).EQ.'4') GO TO 560
      IF(ITREE(ISB,JSB).EQ.'8' ) IJUMP=2
      CALL MNTMBD(JUMP,IJUMP,ISB,JSB,ARJ,AII,I,J,UBJ,VBJ,
     1   FXB,FYB,NCATG,CSHLTR,NPLS)
      GOTO 610
 630  UB(ISB,JSB)=0.
      VB(ISB,JSB)=0.
 610  CONTINUE
C -------------------------------------------------------------------
 560  CONTINUE
 570  CONTINUE
 580  CONTINUE
C       FOR CIRCULAR ISLAND, PERIODIC BOUNDARY CONDITIONS ARE USED.
      IF (DOLLAR.EQ.'$') THEN
        DO 700 I=1,IMXB
        UB(I,1)=UB(I,JMXB)
 700    VB(I,1)=VB(I,JMXB)
       ENDIF
      RETURN
       END
      SUBROUTINE MNTMBD(JUMP,IJUMP,ISB,JSB,ARJ,AII,I,J,UBJ,VBJ,
     1   FXB,FYB,NCATG,CSHLTR,NPLS)
C        J. CHEN   OCTOBER 1980   MDL   IBM 360/195
C
C        PURPOSE
C           THIS SUBROUTINE (MNTMBD) COMPUTES TRANSPORTS ON MOMENTUM
C           BOUNDARY POINTS, BETWEEN LAND AND DEEP WATERS. IN SHALLOW
C           WATERS THE BOUNDARY CONDITION EQUATES SURFACE GRADIENT AT
C           BOUNDARY AND NEAREST INTERIOR POINT.  IN INTERMEDIATE
C           WATERS THE BOUNDARY CONDITION IGNORES TIME VARIATIONS AND
C           EQUATES SURGE GRADIENT AS STATIC HEIGHT GRADIENT AT A
C           BOUNDARY POINT.
C           BELOW IS ARRANGEMENT OF MOMENTUM BOUNDARY LINE
C           AND NEAREST INTERIOR LINE.
C
C          .    .(I,J)          .    .(I,J)      .    . INTERIOR LINE
C                                                *    *
C            +                    +              * +  *
C                                                *    *
C          .    .<---IJUMP=2--->.    .< IJUMP=1 >.****. BOUNDARY LINE
C               (ISB,JSB)          (ISB,JSB)
C  DEEP--->      <INTERMEDIATE >      <SHALLOW  >     <---LAND-----
C
C        DATA SET USE
C           NONE
C
C        VARIABLES
C             ISKIP = 1 TO COMPUTE MOMENTUM ON J=1    BOUNDARY
C                     2 TO COMPUTE MOMENTUM ON I=IMXB BOUNDARY
C                     3 TO COMPUTE MOMENTUM ON J=JMXB BOUNDARY
C         ISB JSB = REVISED SUBSCRIPTS FOR A BOUNDARY POINT
C     UB VB(  ,   ) = TRANSPORTS
C              JUMP = 1 FOR SURFACE FORCING,2 FOR NO SURFACE FORCING
C             IJUMP = 1 FOR SHALLOW WATER BOUNDARY COMPUTATIONS
C                     2 FOR INTER DEPTH   BOUNDARY COMPUTATIONS
C           FXB FYB = COMPONENTS OF DRIVING FORCES
C       AIT AII AIU = PRE-COMPUTED VALUES IN 'MOMNTM'
C       AIV ARJ XTT = DITTO
C            AINTRP = DITTO
C           X1B(10) = FIXED CONSTANTS
C           UBJ VBJ = TRANSPORTS AT (I,J) AT PRESENT TIME,SAVED BEFORE
C                     UPDATED
C
C        GENERAL COMMENTS
C           THIS SUBROUTINE RESIDES IN OVERLAY 'CMPUTE'. IT IS CALLED
C           IN BY SUBROUTINE 'MOMNTM'. SUBROUTINE 'MOMNTM' IS CALLED
C          IN BY SUBROUTINE 'CMPUTE'. THIS SUBROUTINE ALSO CALLS IN
C           SUBROUTINE 'STRESS'.
C
      PARAMETER (M$=169,N$=239,L$=1700,NCT$=1700,LC$=L$)
      COMMON /DUMB5/  ZB(M$,N$),ZBM(M$,N$)
      COMMON /DUMB7/  UB(M$,N$),VB(M$,N$),HB(M$,N$)
      COMMON /SCRTCH/ HSUB(M$,N$)
      COMMON /DUMB10/ HMX(M$,N$)
      COMMON /ELPDST/ ELPCT(M$),ELPDT(M$),ELPCL(M$),ELPDL(M$)
      COMMON /ELPDS2/ ELPCL2(M$),ELPDL2(M$),SINL2(N$),SINT2(N$)
      COMMON /PLRPST/ SINT(N$),COST(N$),SINL(N$),COSL(N$)
      COMMON /SETMIX/ IS(N$),MS(N$),MF(N$)
      COMMON /SETMX1/ IE(N$),ME(N$)
      COMMON /ITREE/  ITREE(M$,N$)
      COMMON /KTREE/  KTREE(L$)
      CHARACTER*1     ITREE,KTREE
      COMMON /KSPCON/ KSKP(M$,N$)
      CHARACTER*1     KSKP
      COMMON /DUMB18/ ISQR(L$),JSQR(L$),ISIDE(L$),HWEIR(L$),
     1                ZBMIN(L$)
      COMMON /ACT1D/  FLWSQR(L$),F1DACT(L$)
      CHARACTER*1     F1DACT
      COMMON /DUMCUT/ CUTL(NCT$),CUTLI(NCT$),CUTLE(NCT$)
C
      COMMON /CS1DA/  COS1D(LC$),SIN1D(LC$)
      COMMON /BANK2/  DELCUT(LC$),BANK(LC$,2)
      COMMON /FORCE/ WX(M$,N$),WY(M$,N$),PP(M$,N$),PPX(M$,N$),PPY(M$,N$)
C
      COMMON /FRTH/   IC12,IC19,BT
      COMMON /EGTH/   DELS,DELT,G,COR
      COMMON /DUMB4/  X1B(10)
      COMMON /SCND/   AR(900),AI(900),BR(900),BI(900),CR(900),CI(900)
C
      AIUL=AII*(UB(ISB,JSB)-UBJ)
      AIVB=AII*(VB(ISB,JSB)-VBJ)
C       TEST FOR PARTIAL FLOODING (NO FORCING).
         IF (JUMP.EQ.2) THEN
      XTTX=0.
      XTTY=0.
      GO TO 140
         ENDIF
C           TEST FOR SHALLOW WATER BOUNDARY AND INTERMEDIATE WATER BRY.
      IF (IJUMP.EQ.2) GOTO 160
C
C        SHALLOW WATER COMPUTATIONS, FORCING COMPUTED ON BOUNDRY PTS
C
      FXD=FXB
      FYD=FYB
C        GRID TRANSFORMATION FROM (I,J) TO (X,Y)
      XR=ELPCT(ISB)*COST(JSB)
      YR=ELPDT(ISB)*SINT(JSB)
      NPLS=1
C      CALL FRCPNT(XR,YR,FXBP,FYBP,NCATG,CSHLTR,NPLS)
C
      CALL FRCPNT(ISB,JSB,FXBP,FYBP,NCATG,CSHLTR,NPLS)
C        CHANGE FORCING FROM (X,Y) SYSTEM BACK TO IMAGE PLANE
      XR=ELPDT(ISB)*COST(JSB)
      YR=ELPCT(ISB)*SINT(JSB)
      FXB=XR*FXBP+YR*FYBP
      FYB=XR*FYBP-YR*FXBP
C        STRESS DIFFERENCE BETWEEN BOUNDARY POINT AND INTERIOR POINT
      XTTX=DELT*(FXB-FXD)
      XTTY=DELT*(FYB-FYD)
C
C        COMPUTE MOMENTUM ON SHALLOW DEPTH BOUNDARIES. SURFACE GRADIENT
C        TERMS CANCELED BY ASSUMED B.C.
C        NOTE: USE SIELECKI'S SCHEME FOR CORIOLIS TERM.
C              VBJ = V-COMPONENT AT PRESENT TIME AT I,J
C              UB(I,J),VB(I,J) = U,V IN FUTURE TIME
  140 UB(ISB,JSB)=UB(I,J)+AIUL+ARJ*(VB(ISB,JSB)-VBJ)    +XTTX
      VB(ISB,JSB)=VB(I,J)+AIVB-ARJ*(UB(ISB,JSB)-UB(I,J))+XTTY
      GO TO 170
C
C        NO FORCIING WHEN JUMP=2
C
C        COMPUTE MOMENTUM ON INTERMEDIATE DEPTH BOUNDARIES
C        ELIMINATE STATIC PRESSURE TERM IN FRCPNT(NPLS=0).
  160 NPLS=0
      CSHLTR=1.
C        GRID TRANSFORMATION FROM (I,J) TO (X,Y)
      XR=ELPCT(ISB)*COST(JSB)
      YR=ELPDT(ISB)*SINT(JSB)
C
C      CALL FRCPNT(XR,YR,FXBP,FYBP,NCATG,CSHLTR,NPLS)
C
      CALL FRCPNT(ISB,JSB,FXBP,FYBP,NCATG,CSHLTR,NPLS)
C
C        CHANGE FROM (X,Y) SYSTEM TO IMAGE PLANE
      XR=ELPDT(ISB)*COST(JSB)
      YR=ELPCT(ISB)*SINT(JSB)
      FXB=XR*FXBP+YR*FYBP
      FYB=XR*FYBP-YR*FXBP
C
  166  CONTINUE
C      IF (UB(I,J).GT.0.) THEN
C            TEST CODES FOR RADIATION B.C.
C      I1=ISB-1
C      J1=JSB-1
C      HBB=.5*(HB(I1,J1)+HB(I1,JSB))-.5*(PP(ISB,JSB)+PP(I1,JSB))-SEADTM
C      UB(ISB,JSB)=5.65*SQRT(.5*(ZB(I1,J1)+ZB(I1,JSB)))*HBB
C
C       ELSE
      UB(ISB,JSB)=AII*UB(ISB,JSB)+ARJ*VB(ISB,JSB)+ DELT*FXB
C      ENDIF
      VB(ISB,JSB)=AII*VB(ISB,JSB)-ARJ*UB(ISB,JSB)+ DELT*FYB
  170 RETURN
       END
      SUBROUTINE HISTRY
C        J. CHEN   OCTOBER        MDL   IBM 360/195
C
C        PURPOSE
C           THIS SUBROUTINE (HISTRY) SMOOTHS THE SURGE AT 10
C           PRE-SELECTED GRID POINTS AND SAVES THEM ON TAPE FOR FUTURE
C           PRINT OUT. ALSO, WIND SPEED AND DIRECTION ARE ALSO SAVED AT
C           THE SAME 10 POINTS.
C        DATA SET USE
C           FT 10, TEMPORARY SCRATCH TAPE FOR SURGE STORAGE
C           FT 20, TEMPORARY SCRATCH TAPE FOR WIND SPEED STORAGE
C           FT 30, TEMPORARY SCRATCH TAPE FOR STORAGE OF WIND DIRECTION
C
C        VARIABLES
C   IPN(10) IPL(10) = SUBSCRIPTS FOR 10 SELECTED POINTS
C        HB(  ,  ) = SURGE FIELD
C        ZB(  ,   ) = DEPTH FIELD
C       IP(4) JP(4) = SHIFT SUBSCRIPTS TO 4 ADJACENT MOMENTUM POINTS
C       ZBM(  ,   ) = MAXIMUM BARRIER HEIGHT AT A MOMENTUM POINT
C COST(  ) SINT(  ) = RAY PROJECTIONS THRU MOMNTM POINTS, ON (X,Y) AXIS
C        GRIDRF(  ) = RADIAL DISTANCES TOMOMENTUM POINTS
C             AX AY = CMPNENTS OF TOTAL STRM TRVRSE, ADVNCD IN 'STMVAL'
C             C1 C2 = COMPONENTS OF INITIAL STORM POSITION
C             NCATG = 1 FOR LAKE WINDS, 2 FOR OCEAN WINDS
C          W1218(2) = TWICE THE MAXIMUM WINDS, LAKE/OCEAN WINDS
C   CW(800) SW(800) = INFLOW ANGLE FOR LAKE WINDS
C     C(800) S(800) = INFLOW ANGLE FOR OCEAN WINDS
C           X12(50) = FIXED COEFFICIENTS
C    X12(9) X12(10) = COMPONENTS OF FORWARD SPEED, 1ST HOUR
C          WIND(10) = WIND SPEED AT 10 SELECTED POINTS
C          TDIR(10) = WIND DIRECTION AT 10 SELECTED POINTS
C               PHI = SLANT, Y-AXIS TO N/S DIRECTION, AZMTH+270
C               RAD = RADIAN MEASURE
C             IPRHR = 1,WHEN ON THE HOUR; 2, WHEN HALF HOUR IN BETWEEN;
C                     3, OTHERWISE
C
C        GENERAL COMMENTS
C           THIS SUBROUTINE RESIDES IN OVERLAY 'CMPUTE'. IT IS CALLED
C           IN BY SUBROUTINE 'CMPUTE' AT SELECTED TIME INTERVALS.
C
      PARAMETER (M$=169,N$=239,L$=1700,NCT$=1700,LC$=L$)
      COMMON /DUMB5/  ZB(M$,N$),ZBM(M$,N$)
      COMMON /DUMB7/  UB(M$,N$),VB(M$,N$),HB(M$,N$)
      COMMON /SCRTCH/ HSUB(M$,N$)
      COMMON /DUMB10/ HMX(M$,N$)
      COMMON /ELPDST/ ELPCT(M$),ELPDT(M$),ELPCL(M$),ELPDL(M$)
      COMMON /ELPDS2/ ELPCL2(M$),ELPDL2(M$),SINL2(N$),SINT2(N$)
      COMMON /PLRPST/ SINT(N$),COST(N$),SINL(N$),COSL(N$)
      COMMON /SETMIX/ IS(N$),MS(N$),MF(N$)
      COMMON /SETMX1/ IE(N$),ME(N$)
      COMMON /ITREE/  ITREE(M$,N$)
      COMMON /KTREE/  KTREE(L$)
      CHARACTER*1     ITREE,KTREE
      COMMON /KSPCON/ KSKP(M$,N$)
      CHARACTER*1     KSKP
      COMMON /DUMB18/ ISQR(L$),JSQR(L$),ISIDE(L$),HWEIR(L$),
     1                ZBMIN(L$)
      COMMON /ACT1D/  FLWSQR(L$),F1DACT(L$)
      CHARACTER*1     F1DACT
      COMMON /DUMCUT/ CUTL(NCT$),CUTLI(NCT$),CUTLE(NCT$)
C
      COMMON /CS1DA/  COS1D(LC$),SIN1D(LC$)
      COMMON /BANK2/  DELCUT(LC$),BANK(LC$,2)
C
      COMMON /FORCE/ WX(M$,N$),WY(M$,N$),PP(M$,N$),PPX(M$,N$),PPY(M$,N$)
      COMMON /FRTH/   IC12,IC19,BT
      COMMON /PRHSY/  IPRHR,JUMPR,IPRT
      COMMON /CESAV/  IPN(120),IPL(120),KHSPT
      COMMON /CESAV1/ STATS(120)
      COMMON /STRMSB/ C1,C2,C21,C22,AX,AY,PTENCY,RTENCY
      COMMON /SSTME/  ISTM,JHR,ITMADV,NHRAD,IBGNT,ITEND
      COMMON /FRST/   X1(50),X12(50)
      COMMON /BSN/    PHI,ALTO,ALNO,PHI1,ALT1,ALN1,ALT1C
C---------------------------------------------------------
      COMMON /SUBDCE/ ZSUBCE,JSUB,JSUB1,ZSUB(N$)
      COMMON /GPRT1/  DOLLAR,EBSN
      DIMENSION       WIND(120),TDIR(120)
      DIMENSION       SAVEH(120)
      CHARACTER*2     DOLLAR
      CHARACTER*1     EBSN
C---------------------------------------------------------
      CHARACTER*30    STATS
      DATA KHR,KHR2,KHR3/1,0,0/
C
C
C       KEEP TRACK TIME, ON THE HOUR,HALF HOUR BETWEEN, OR ELSE
C       FROM MOD(IME,100)=1,OR 2,OR 3
C
       GO TO (10,20,30),IPRHR
   10  KHR=KHR+1
       IME=KHR*100+1
       GO TO 40
   20  KHR2=KHR2+1
       IME=KHR2*100+2
       GO TO 40
   30  KHR3=KHR3+1
      IME=KHR3*100+3
   40  CONTINUE
C********SMOOTHING HISTORICAL SURGES AT 10 SELECTED GRID POINTS*********
      DO 1190 L=1,KHSPT
      I=IPN(L)
      J=IPL(L)
C
      SAVEH(L)=HB(I,J)
C
      A=WX(I,J)
      B=WY(I,J)
      WIND(L)=SQRT(A*A+B*B)/3.28084
      IF (WIND(L).LT.1.E-5) THEN
       Z=180.
       ELSE
       Z =ATAN2(A,B)/1.74532925E-2+180.
       ENDIF
      TDIR(L)=AMOD(Z+PHI,360.)
 1190 CONTINUE
C        WRITING HISTORICAL SURGES AT SELECTED GRID PTS ONTO TAPE
      WRITE (81)   IME,(SAVEH(N),N=1,KHSPT)
C      WRITE(*,'(I5,10F6.2,F10.2)')   IME,(SAVEH(N),N=1,10),BT
      WRITE (82)   IME,(WIND(N),N=1,KHSPT)
      WRITE (83)   IME,(TDIR(N),N=1,KHSPT)
C       WRITE(*,'(5X,10F6.2,F10.2)')   (WIND(N),N=1,10)
C      WRITE(*,'(5X,10F6.2,F10.2)')   (TDIR(N),N=1,10)
      RETURN
      END
      SUBROUTINE FILTER
C        SEPTEMBER 1980    JYE CHEN    MDL   IBM 360/195
C        PURPOSE
C           TO SMOOTH THE SURGE HEIGHTS ONCE AN HOUR TO ELIMINATE
C           2-INTERVAL NOISES IN THE FIELDS.
C
C        DATA SET USE
C           NONE
C
C        VARIABLES
C       HSUB(  ,  ) = SCRATCH SPACES
C       IP(4) JP(4) = SHIFTS OF I J FROM HEIGHT PT TO 4 MOMNTN CORNERS
C          IBB(6,M) = RANGES AND INCREMENTS FOR 3 DO-LOOPS OF I,J
C                     M=1 INTERIOR,M=2 HORZNTL BNDY,M=3 VERTICAL BNDY
C                     RANGES ON THREE BASIN SEGMENTS ARE:
C
C                              222222222........2222222
C                              311111111........1111113
C                              311111111........1111113
C                              .                      .    J
C                              .                      .    J
C                              311111111........1111113    J
C                              311111111........1111113    J
C                              222222222........2222222    ---->I
C
C          J2(J) JC = SHIFT J-SUBSCRIPT TO PRESENT TIME OF SURGE FIELD
C        HB(  ,   ) = SURGE FIELD
C         ZB(  ,  ) = DEPTH FIELD
C     IIH(4) JJH(4) = I/J-SHIFT SUBSCRIPTS FOR SURGE POINTS
C       IP(4) JP(4) = I/J-SHIFT SUBSCRPTS FOR MOMNTM PTS, IN 'BLCKDT'
C
C           THE (I,J) SHIFTS TO MOMENTUM POINTS, VIA IP(4) AND JP(4),
C           ON CORNER POINTS K=1 TO 4 ARE:
C
C                        I+0,J+1      I+1,J+1
C                           K=4.-----.K=3
C                              I     I
C                              I I*J I
C                              I     I
C                           K=1.-----.K=2
C                        I+0,J+0      I+1,J+0
C
C        ZBM(  ,  ) = MAX BARRIER HEIGHTS AT MOMENTUM POINTS
C        GENERAL COMMENTS
C           THIS SUBROUTINE RESIDES IN OVERLAY 'CMPUTE'.
C           A 5-POINT SMOOTHING OPERATOR IS USED, NEIGHBORING GRIDS
C           WEIGHED CONDITIONALLY IF THEY ARE DRY OR BLOCKED
C           BY BARRIERS. SHIFTS OF (I,J) TO ADJACENT SQUARES ARE SET
C           VIA IIH(4) AND JJH(4). IIH/0,1,0,-1/,JJH/-1,0,1,0/
C
C
      PARAMETER (M$=169,N$=239,L$=1700,NCT$=1700,LC$=L$)
      COMMON /DUMB5/  ZB(M$,N$),ZBM(M$,N$)
      COMMON /DUMB7/  UB(M$,N$),VB(M$,N$),HB(M$,N$)
      COMMON /SCRTCH/ HSUB(M$,N$)
      COMMON /DUMB10/ HMX(M$,N$)
      COMMON /ELPDST/ ELPCT(M$),ELPDT(M$),ELPCL(M$),ELPDL(M$)
      COMMON /ELPDS2/ ELPCL2(M$),ELPDL2(M$),SINL2(N$),SINT2(N$)
      COMMON /PLRPST/ SINT(N$),COST(N$),SINL(N$),COSL(N$)
      COMMON /SETMIX/ IS(N$),MS(N$),MF(N$)
      COMMON /SETMX1/ IE(N$),ME(N$)
      COMMON /ITREE/  ITREE(M$,N$)
      COMMON /KTREE/  KTREE(L$)
      CHARACTER*1     ITREE,KTREE
      COMMON /KSPCON/ KSKP(M$,N$)
      CHARACTER*1     KSKP
      COMMON /DUMB18/ ISQR(L$),JSQR(L$),ISIDE(L$),HWEIR(L$),
     1                ZBMIN(L$)
      COMMON /ACT1D/  FLWSQR(L$),F1DACT(L$)
      CHARACTER*1     F1DACT
      COMMON /DUMCUT/ CUTL(NCT$),CUTLI(NCT$),CUTLE(NCT$)
C
      COMMON /CS1DA/  COS1D(LC$),SIN1D(LC$)
      COMMON /BANK2/  DELCUT(LC$),BANK(LC$,2)
C
      COMMON /FFTH/   ITIME,MHALT
      COMMON /DUMB3/  IMXB,JMXB,IMXB1,JMXB1,IMXB2,JMXB2
      COMMON /EGTH/   DELS,DELT,G,COR
      COMMON /DUMB8/  IP(4),JP(4),IH(4),JH(4)
      COMMON /DUM88/  IIH(4),JJH(4),IHH(4),JHH(4)
      COMMON /GAMAF/  GAMA,GAMA1,GAMFP(4)
      COMMON /GPRT1/  DOLLAR,EBSN
      DIMENSION       GAMF(4)
      DIMENSION       IBB(6,3)
      CHARACTER*2     DOLLAR
      CHARACTER*1     EBSN
      DATA HCRT/1./
      DATA IBB(1,1),IBB(1,2),IBB(1,3)/2,1,2/,
     1     IBB(3,1),         IBB(3,3)/1,  1/,
     2     IBB(4,1),IBB(4,2),IBB(4,3)/2,1,1/,
     3     IBB(6,1),IBB(6,2)         /1,1  /
      IBB(2,1)=JMXB2
      IBB(2,3)=JMXB2
      IBB(3,2)=JMXB2
      IBB(2,2)=JMXB1
      IBB(5,1)=IMXB2
      IBB(6,3)=IMXB2
      IBB(5,2)=IMXB1
      IBB(5,3)=IMXB1
      DO 100 J=1,JMXB1
      DO 100 I=1,IMXB1
 100  HSUB(I,J)=HB(I,J)
      DO 210 M=1,3
      JL  =IBB(1,M)
      JM  =IBB(2,M)
      JNCR=IBB(3,M)
      IL  =IBB(4,M)
      IM  =IBB(5,M)
      INCR=IBB(6,M)
      DO 200 J=JL,JM,JNCR
      DO 190 I=IL,IM,INCR
C        SKIP COMPUTATIONS FOR TERRAIN HIGHER THAN 35 FT.
      IF(ZB(I,J).LE.-35.) GO TO 190
      HST=HSUB(I,J)+ZB(I,J)
C        TEST FOR LAND WETTED LESS THAN 1 FT.
      IF(HST.LT.HCRT) GO TO 190
      HSM=0.
      HDIF=0.
C
C******** 'LOOP' FOR SMOOTHING******************************************
      GAM0=ELPDL2(I)+(ELPCL2(I)-ELPDL2(I))*SINL2(J)
      DO 170 K=1,4
      II=I+IIH(K)
      JJ=J+JJH(K)
      IF (M.EQ.1) GOTO 158
        IF (II.EQ.0.OR.JJ.EQ.0) THEN
        GAMF(K)=0.
        GO TO 160
        ENDIF
 158  CONTINUE
      IF (EBSN.EQ.'$'.OR.EBSN.EQ.'+') THEN
      GAMFK=ELPDL2(II)+(ELPCL2(II)-ELPDL2(II))*SINL2(JJ)
      Z=GAMFK/GAM0
      GAMF(K)=2.*Z/(1.+Z)-1.
      ELSE
      GAMF(K)=GAMFP(K)
      ENDIF
      IF (II.EQ.IMXB.OR.JJ.EQ.JMXB) GOTO 160
C        THE NEIGHBORING SQUARE IS MERELY WET(LESS THAN 1 FT).
      IF (HSUB(II,JJ)+ZB(II,JJ).LT.HCRT) GO TO 160
      KK=MOD(K,4)+1
      IA=I+IP(K)
      IB=I+IP(KK)
      JA=J+JP(K)
      JB=J+JP(KK)
      Z =ZBM(IA,JA)
      ZZ=ZBM(IB,JB)
      ZZZ=AMIN1(Z,ZZ)+1.
C        WATER ON EITHER SIDE MUST EXCEED THE BARRIER(LOWER) BY AT
C        LEAST 1 FT.
      IF( HSUB(I,J).LT.ZZZ.OR.HSUB(II,JJ).LT.ZZZ) THEN
          HSM=HSM+HSUB(I,J)
          HZZ=HSUB(I,J)
          ELSE
          HSM=HSM+HSUB(II,JJ)
          HZZ=HSUB(II,JJ)
          ENDIF
      GOTO 212
 160  CONTINUE
          HSM=HSM+HSUB(I,J)
          HZZ=HSUB(I,J)
 212  CONTINUE
C       WEIGHTING ADJUSTMENTS IN I-DIRECTION DUE TO POLAR GRIDS.
       HDIF=HDIF+HZZ*GAMF(K)
C
 170  CONTINUE
      HTEMP=.5*HSUB(I,J)+(HSM+HDIF)/8.
C******* END 'LOOP' FOR SMOOTHING **************************************
C
 180  HB(I,J)=AMAX1(-ZB(I,J),HTEMP)
 190  CONTINUE
 200  CONTINUE
 210  CONTINUE
      RETURN
       END
      SUBROUTINE PLOTDP
C        J. CHEN   OCTOBER 1980   MDL   IBM 360/195
C           DIRECTION AT 10 MINUTE INTERVALS FOR SELECTED STATIONS, AS
C           PREVIOUSLY STORED ON TAPE BY SUBROUTINE 'HISTRY'. HEADINGS
C           FOR THE HOUR, 30 MINUTES, AND 10 MINUTES INTERVALS ARE
C           COMPUTED SEPARATELY.
C
C        DATA SET USE
C           FT10, TEMPORARY SCRATCH FILE FOR SURGE STORAGE
C           FT20, TEMPORARY SCRATCH FILE FOR WIND SPEED STORAGE
C           FT30, TEMPORARY SCRATCH FILE FOR WIND DIRECTION STORAGE
C
C        VARIABLES
C          NDATA(3) = UNIT NUMBERS FOR 3 FILES
C         STATS(10) = ALPHAMERIC NAMES FOR 10 STATIONS
C             IBGNT = HOUR TO BEGIN COMPUTAIONS
C      IDAY<-->ISTM = TIME ENTRIES FOR 'JULHR'
C       IPN(10) IPL = GRID POSITIONS OF 10 STATIONS
C        IST<-->NHR = TIME OUTPUTS FROM 'DATEX'
C            KDAYWK = 3-LETTER ALPHAMERIC NAME FOR DAY OF WEEK
C          AMON(12) = 3-LETTER ALPHAMERIC NAME FOR MONTH OF THE YEA
C            IPRTAD = TIME BETWEEN PRINT OUTS
C         IME ITIME = TIME OF A PRINT OUT
C         ALINE(10) = PRINT OUT VALUES ON TAPE FOR THE 10 STATIONS
C            NDLTHR = NUMBER OF TIME STEPS IN ONE HOUR
C
C        GENERAL COMMENTS
C           THIS SUBROUTINE RESIDES IN OVERLAY 'PLOTS'. IT IS CALLED
C           IN BY SUBROUTINE 'PLOTS'.
C
      COMMON /IDENT/  AIDENT(40),DACLOK(7)
      CHARACTER*4   AIDENT
      CHARACTER*18 AAA
      COMMON /DMN/    IDAY,NMNTH,IYEAR,IMTH
      COMMON /CESAV/  IPN(120),IPL(120),KHSPT
      COMMON /CESAV1/ STATS(120)
      CHARACTER*30    STATS
      COMMON /SSTME/  ISTM,JHR,ITMADV,NHRAD,IBGNT,ITEND
      COMMON /MTH/    MONTH(12),AMON(12)
      COMMON /PRTSV/  IPRTSV,IPRTAD
      COMMON /FFTH/   ITIME,MHALT
      COMMON /DELTH/  NDLTHR,NDLTH2,DLTHR,DLTH2
      COMMON /DATE7/  IDATE(4)
      CHARACTER*4    DACLOK
      DIMENSION       NDATA(3)
      DIMENSION       ALNE(1000,120),IME(1000)
      CHARACTER*4 AMON,AM,KDAYWK
      DATA NDATA/51,52,53/
C
C      CALL DATIM2(DACLOK)
C**  DATIM2 IS STANDARD NOAA FORMATTED DAY-TIME ROUTINE. IT RETURNS THE
C**  CALENDAR DAY AND TIME OF RUN IN FORMAT 'BMMMBDD,BYYYYB-BHH:MM:SS.SS
C**  B' WHERE B IS BLANK, MMM IS 3-LETTER MONTH NAME, DD IS 2-DIGIT DAY,
C**  HH IS HOUR, MM IS MINUTE, AND SS.SS IS SECONDS TO HUNDREDTHS.
C
C     PRINT 101,DACLOK
 101  FORMAT(1H1,1H0,10X,' CALCULATIONS PERFORMED AT DAY - TIME ',7A4//)
C
       DO 170 NNN=1,1
      LUNIT=NDATA(NNN)
      REWIND LUNIT
      L=1
 100  READ (51,END=102) IME(L),(ALNE(L,I),I=1,KHSPT)
      L=L+1
      GOTO 100
 102  LMX=L
      MMM=(KHSPT-1)/10.+1
      DO 172 M=1,MMM
        M1=1+(M-1)*10
        M2=MIN(M1+9,KHSPT)
C        SKIP ONE PAGE, PRINT HEADINGS FOR SURGE, WIND SPEED/DIRECTION
       WRITE (*,'(1H1)')
       WRITE (*,'(A)') ' EXTRA-TROPICAL STORM SURGE FORECASTS '
       AAA(1:3)=AIDENT(5)(2:4)
       AAA(4:7)=AIDENT(6)(1:)
       AAA(8:11)=AIDENT(7)(1:)
       AAA(12:15)=AIDENT(8)(1:)
       AAA(16:17)=AIDENT(9)(1:2)
       WRITE (*,'(A,A,A)') ' COMPUTATION STARTING FROM ',AAA,
     1 '  --- OUTPUT IN EASTERN STANDARD TIME(-5 HOUR FROM UTC).'
         IF (NNN.EQ.1) WRITE(*,*)'           TIME HISTORIES OF SURGES',
     1 ' AT SELECTED GAGES '
         IF (NNN.EQ.2) WRITE (*,*)'           COMPUTED WIND SPEEDS',
     1 ' AT SELECTED GAGES OR GRID POINTS'
         IF (NNN.EQ.3) WRITE(*,*)'           COMPUTED WIND DIRECTIONS',
     1 ' AT SELECTED GAGES OR GRID POINTS'
      WRITE (*,220)(STATS(I),I=M1,M2)
C
C        PRINT STATION GRID POSITIONS
      WRITE (*,230)(IPN(I),IPL(I),I=M1,M2)
      WRITE (*,'(1H0)')
C
      JST=IBGNT-JHR
      REWIND LUNIT
      DO 160 L=1,LMX
      KPR=MOD(IME(L),100)
C     WRITE (*,*) (IDATE(I),I=1,4),JHR,IBGNT
      KHR=(IME(L)-KPR)/100
C-----------------------------------------------------------------------
C       OUTPUT ONLY FROM ANALYSES TO 36TH-HR FORECAST.
C       THIS IS FOR OPERATIONAL OUTPUT, IN HALF HOUR INTERVAL.
C-----------------------------------------------------------------------
      IF (KHR+JST.LT.1.AND.(ITEND-JHR).GT.0) GOTO 160
C-----------------------------------------------------------------------
      GO TO (920,155,155),KPR
 155  CONTINUE
      WRITE (*,910)(ALNE(L,I),I=M1,M2)
 910  FORMAT(19X,10F10.2)
      GO TO 160
 920  J=KHR+JST-6
C        PRINT OUT THE HOURLY VALUES
      CALL SSDATE(NHR,NDA,NMO,NYR,IDATE,J)
      AM=AMON(NMO)
      NHR00=NHR
      IMP=NHR00/10
      NHR00=NHR00-IMP*10
      IF (MOD(NHR,12).EQ.0) THEN
      WRITE (*,260) NDA,AM,NYR,IMP,NHR00,(ALNE(L,I),I=M1,M2)
  260 FORMAT(1X,I3,2X,A3,I5,',',1X,I1,I1,'E',10('----',F6.2))
      ELSE
      WRITE (*,261) IMP,NHR00,(ALNE(L,I),I=M1,M2)
  261 FORMAT(15X,1X,I1,I1,1X,10('----',F6.2))
      ENDIF
 160  CONTINUE
 172   CONTINUE
C      CLOSE (LUNIT)
  170 CONTINUE
  300 FORMAT(1X,20A4,/ 1X,20A4)
  220 FORMAT(/23X,10A10)
  230 FORMAT(//'  DAY MTH YEAR HOUR',3X,10(I4,',',I3,2X))
  250 FORMAT(19X,10F10.2)
  270 FORMAT(2X,I1,I1,'30',13X,        10F10.2)
      RETURN
      END
      SUBROUTINE SSDATE(ITHR,ITDY,ITMO,ITYR,IDATE,JINC)
C        PURPOSE - TO CALCULATE THE CALENDAR DATES/TIMES FROM
C                  INPUT IDATE(4) ADVANCED BY 'JINC' HOURS
      DIMENSION  JDATE(4),NDAYS(12),IDATE(4)
      DATA NDAYS/31,28,31,30,31,30,31,31,30,31,30,31/
C       RESET NDAYS(2)=29 IF IT STARTS AS A LEAP YEAR
      IF(IDATE(1).EQ.4*(IDATE(1)/4))NDAYS(2)=29
      IF(IDATE(1).EQ.100*INT(IDATE(1)/100)) NDAYS(2)=28
      IF(IDATE(1).EQ.400*INT(IDATE(1)/400)) NDAYS(2)=29
C        RECORD THE STARTING TIME
      DO 100 J=1,3
 100  JDATE(J)=IDATE(J)
      NDDD=(JINC+IDATE(4))/24
      ITHR=JINC+IDATE(4)-NDDD*24
      IF (ITHR.LT.0) THEN
        NDDD=NDDD-1
        ITHR=ITHR+24
        ENDIF
      JDATE(3)=JDATE(3)+NDDD
C
 300  IF(JDATE(3).LE.NDAYS(JDATE(2)))GOTO 500
      JDATE(3)=JDATE(3)-NDAYS(JDATE(2))
      JDATE(2)=JDATE(2)+1
C        IS THE MONTH DEC. OR ERLIER?
      IF(JDATE(2).LE.12) GOTO 300
      JDATE(2)=JDATE(2)-12
      JDATE(1)=JDATE(1)+1
      NDAYS(2)=28
C        IS IT A LEAP YEAR?
      IF(JDATE(1).EQ.4*(JDATE(1)/4))NDAYS(2)=29
      IF(JDATE(1).EQ.100*INT(JDATE(1)/100)) NDAYS(2)=28
      IF(JDATE(1).EQ.400*INT(JDATE(1)/400)) NDAYS(2)=29
 400  GOTO 300
 500  CONTINUE
C
      ITDY=JDATE(3)
      ITMO=JDATE(2)
      ITYR=JDATE(1)
      RETURN
      END
      SUBROUTINE CMPUTE
C
C        J. CHEN   SEPTEMBER 1980 MDL   IBM 360/195
C
C        PURPOSE
C           THIS 2ND OVERLAY 'CMPUTE' USES SUBROUTINE 'CMPUTE' AS
C           A DIRECT CALLING SEQUENCE FOR 12 SUBROUTINES.
C           THE PURPOSE OF THE OVERLAY IS TO COMPUTE SURGE HEIGHTS
C           USING ALGORITHMS FOR PARTIAL DIFFERENTIAL EQUATIONS OF
C           MOTION. SURGES ARE STORED FOR LATER PRINT OUT.
C
C        DATA SET USE
C           NONE
C
C        VARIABLES
C             ITIME = CUMULATIVE TIME COUNTER
C             MHALT = HALT COMPUTATIONS AT THIS TIME STEP
C              NPSS = NUMBER OF CHANNELS
C              NCUT = NUMBER OF CUTS ALONG BARRIORS
C            NSQRWC = NUMBER OF SQUARES FOR 1D FOLW PLUS FLOW CUTS
C            NDLTHR = NUMBER OF TIME STEPS IN ONE HOUR
C            IPRTSV = TIME TO SAVE HISTRCL SURGES, 10 PTS, SET 'INTVAL'
C            IPRTAD = UPDATE TIME OF IPRTSV, SET IN 'INTVAL'
C
C        GENERAL COMMENTS
C
      PARAMETER (IMX$=36,JMX$=50)
      COMMON /DATINC/ XLAT1,XLON1,DELLT,DELLN,IMXC,JMXC
      COMMON /DUMB3/  IMXB,JMXB,IMXB1,JMXB1,IMXB2,JMXB2
      COMMON /DELTH/  NDLTHR,NDLTH2,DLTHR,DLTH2
      COMMON /FFTH/   ITIME,MHALT
      COMMON /SSTME/  ISTM,JHR,ITMADV,NHRAD,IBGNT,ITEND
      COMMON /PRTSV/  IPRTSV,IPRTAD
      COMMON /FLWCPT/ NSQRS,NSQRW,NSQRWC,NPSS,NCUT
      COMMON /PRHSY/  IPRHR,JUMPR,IPRT
      COMMON /GPRT1/  DOLLAR,EBSN
      COMMON /GPRT2/  EBSN1,EBSN2
      COMMON /KEYS/   KEY1
      COMMON /SWTCH/  IOPERL(5)
      COMMON /POLAR/  DEGREE,RMOUTH,AZMTH,DELA
      COMMON /ELLIP/  AAXIS,BAXIS,ABQAB
      COMMON /INDX/   IRB,IRE,NDB,NDE
      COMMON /STRMSB/ C1,C2,C21,C22,AX,AY,PTENCY,RTENCY
      COMMON /LANDS/ LANDS(IMX$,JMX$)
      COMMON /IFT11 / LATIN1,LATIN2,LONIN1,LONIN2

      CHARACTER*1    LANDS
      CHARACTER*2    DOLLAR
      CHARACTER*1    EBSN,EBSN1,EBSN2,AIGEO
      CHARACTER*80 TITLE,FT11NM
C       SET ITEST =1 FOR TESTING WIND/PRESSURE INPUT ONLY.
       ITEST=0
C
       IGEO=0
C
C      READ U/V/P DATA DOMAIN INFORMATION, UPPER-LEFT LAT/LON.
C	WRITE (*,*) ' FT11 = '
C	READ(*,'(A)') FT11NM
C      OPEN (11,FILE=FT11NM)
C
CC       READ (11,*) LATIN1,LATIN2
CC       READ (11,'(A1)') ACHR
CC       READ (11,*) LONIN1,LONIN2
       XLAT1=LATIN1
       XLON1=LONIN1
C
       CALL SETCP1(IGEO)
       CALL INTLHT
C      OPEN (51,FILE='FLE10')
C      OPEN (52,FILE='FLE20')
C      OPEN (53,FILE='FLE30')
C
  100 CALL TTIMER
C
      IF (ITEST.NE.1) CALL CONTTY
      IF (ITEST.NE.1) CALL BDRYHT
      IF (ITEST.NE.1) THEN
      IF (NSQRWC.NE.0) CALL FLW1DM
      ENDIF
C
C      CALL HMXSV
C
      IPRHR=3
      IF (ITIME.NE.IPRTSV) GO TO 180
      IF(ITIME.EQ.NHRAD) IPRHR=1
      IF(ITIME.EQ.NHRAD-NDLTHR/2) IPRHR=2
      IPRTSV=IPRTSV+IPRTAD
C
C      CALL HISTRY
C      CHANGE TIME HISTRY OUTPUT IN EVERY HALF HOUR INTERVAL.
      IF (IPRHR.NE.3) CALL HISTRY
  180 CONTINUE
C
C      IF(MOD(ITIME+1,30).EQ.0)  CALL FILTER
C      SMOOTHING IS APPLIED EVERY TIME STEP.
       IF (ITEST.NE.1) CALL FILTER
       IF (ITEST.NE.1) CALL MOMNTM
C
       CALL STMVL1(IGEO)
C
      IF (ITIME.LT.MHALT) GO TO 100
C
       RETURN
       END

       SUBROUTINE FLEGEN
C   INPUT FILES:
C     FORT.11  - ETTGP  TIDE GAGE STATIONS TO BE MONITORED
C     FORT.81  - FLE21   SURGE HEIGHT (FLE10)
C     FORT.82  - FLE22   WIND SPEED   (FLE20)
C     FORT.83  - FLE23   WIND DIRECTION (FLE30)
C   OUTPUT FILE:
C     FORT.52   - SURGE.HH    SURGE ARCHIVE
C
CC      PARAMETER (IBGS=48, ITMS=48)
      COMMON /CESAV/  IPN(120),IPL(120),KHSPT
      COMMON /CESAV1/ STATS(120)
      COMMON /INZSTM/ IBGS,ITMS
      DIMENSION A(120)
      CHARACTER*30 STATS

       MHOUR=IBGS+ITMS
C       READ (11,'(I3)') NSTA
       WRITE (52) KHSPT,MHOUR
       DO 10 N=1,KHSPT
C       READ (11,'(A30,18X,2I4)') STA(N),IN(N),JN(N)
       WRITE (52) STATS(N),IPN(N),IPL(N)
  10   CONTINUE
       NH2=MHOUR*2
       DO 20 N=1,NH2
       READ (81,END=999) IME,(A(L),L=1,KHSPT)
       WRITE (52) IME,(A(L),L=1,KHSPT)
 20    CONTINUE
 999   CONTINUE
       DO 30 N=1,NH2
       READ (82) IME,(A(L),L=1,KHSPT)
       WRITE (52) IME,(A(L),L=1,KHSPT)
 30    CONTINUE
       DO 40 N=1,NH2
       READ (83) IME,(A(L),L=1,KHSPT)
       WRITE (52) IME,(A(L),L=1,KHSPT)
 40    CONTINUE
      RETURN
      END
