WAVZ2.FOR

C     INCLUDE (params - version 2 of WISWAV)
      PROGRAM WAVAD
      INCLUDE 'PARAMS'
      INCLUDE 'PARAMS2'

C  VERSION AS OF 9/30/93 - INCLUDES MODIFIED CURVATURE EFFECTS
C                          AND FETCH VARIATION WITH LATITUDE

C   NOTE THIS CODE USES A PARAMS2 INCLUDE FILE!!
C   THIS CONTAINS DIMENSIONS FOR THE "LARGE" MATRICES IN THE "SMOD"
C   COMMON BLOCK.
C   WHEN USING CODE AS WAVAD -	NO REFRACTION, SET THOSE DIMENSIONS
C   TO 1 IN PARAMS2
C   WHEN USING CODE AS FULL WISWAV -INCLUDING REFRACTION AND SHOALING,
C   SET THOSE DIMENSIONS TO THE SAME VALUES AS THEIR COUNTERPARTS IN
C   PARAMS

      COMMON /DEPN/ DEPNEW(IDMN,JDMN)
      COMMON /IFCST/NFCST
      COMMON /CCG/CA(IDMN,JDMN,IF),CGA(IDMN,JDMN,IF),DELX(JDMN)
      COMMON /OUTTTS/ IKOUTL,ISTWAV
      COMMON /GAMALP/ GAMOUT(IDMN,JDMN),ALPOUT(IDMN,JDMN)
      COMMON /MAXOUT/ WAVMAX(IDMN,JDMN),TPPMAX(IDMN,JDMN),
     1		      WDDMAX(IDMN,JDMN)
      COMMON /CDR/COEF1,COEF2,SPDFAC,XLAMZ,ALPSTR,BETA1
      COMMON /JAC/ ZJAC(IDMN,JDMN,IF)
      COMMON /DEPTH/ DEP(IDMN,JDMN)
      COMMON /IOP/ICURV,IDPRP,IUPDAT
      COMMON /SWL/ HSWL(IDMN,JDMN),TSWL(IDMN,JDMN),SWLDIR(IDMN,JDMN)
      COMMON /TM/ TIMTOT
      COMMON /PMOD/ XD(IF,JDMN,IA),XDU(IF,JDMN,IA),YD(IF,IA),
     1		    YDU(IF,IA)
      COMMON /DPR/ IUB(IA),JUA(IA),IAUA(IA,JDMN),AD(IF,JDMN,IA),
     1		    ADU(IF,JDMN,IA)
      COMMON /STREND/ IDFRST,IDLAST,IDBND1,IDBND2
c      integer*2 e
      COMMON /CEN/ E(IDMN,JDMN,IF,IA),ENXT(NBN,IF,IA),ELST(NBN,IF,IA)
      COMMON /O1/ FREQ(IF),SINR(IA),COSR(IA),ALFA0(IA),ALFAR0(IA),
     2		  DELF(IF),AINC,AINCR
      COMMON /O4/ EN(IDMN,JDMN,IA)
      COMMON /O5/ IBOUND(IDMN,JDMN)
      COMMON /A3/ G,PI,XN,DL,NANG,NFRQ,NX,NY,TWOPI,FAC,DT,RADDEG
      COMMON /OPT/ MSTA,NSTR,NORD,JPRC,NTMS,INDELF,NHR,
     2		  IOUT(70),JOUT(70),KTIMES,WL(725),WN(725),IG21,IR21
     3		  ,MXHR,IDEEP,DLAT,XLAT0,DZRO,NSTA(70),ALT(70),ALON(70)
      COMMON /WN/ IWSL(IDMN,JDMN),IWDL(IDMN,JDMN),IWS(IDMN,JDMN),
     2		  IWD(IDMN,JDMN),IWSN(IDMN,JDMN),IWDN(IDMN,JDMN)
      COMMON /S1/ PANG(IA,72,5)
      COMMON /S2/ FMA(IDMN,JDMN),
     2	     SWANG(IDMN,JDMN),KFRQA(IDMN,JDMN)
      COMMON /S3/ NF1,RADC,PI72,PI8
      COMMON /S4/ HSIG(IDMN,JDMN)
      COMMON /S5/ FMAN(IDMN,JDMN),E2(IF,IA),
     2	     EF(IF),EFN(IF),FKA(IF)
      COMMON /S6/ ES(IDMN,JDMN,IF),IWVDIR(IDMN,JDMN)
      COMMON /S7/ AVANG(IDMN,JDMN),TSS(IDMN,JDMN)
      COMMON /BD/ IBND,NBPO,IBPO(NOBPTS),JBPO(NOBPTS),NBPI,KNTBP,
     2	     IBPI(NBN),JBPI(NBN)
      COMMON /BDN/ FMANXT(NBN),ALPNXT(NBN),SWANXT(NBN),
     2	     GAMNXT(NBN),KFRNXT(NBN),WKMNXT(NBN),IWVNXT(NBN)
      COMMON /BDL/ FMALST(NBN),ALPLST(NBN),SWALST(NBN),
     2	     GAMLST(NBN),KFRLST(NBN),WKMLST(NBN),IWVLST(NBN)
      COMMON /MM/ INPLEV,IWRAP
      COMMON /HBND/  IBEGA(NBN),JBEGA(NBN),IENDA(NBN),JENDA(NBN),
     1 ITB(NBN,12),JTB(NBN,12),NLA(NBN),WTB(NBN,12),
     2	WTA(NBN,12),NSEG,NINT,ILS,JLS,ICBI(NBN),JCBI(NBN),
     3	IF1,JF1,IFN,JFN,NEXF1,NINPT,NINTR,NEXFN,NOF1,NOFN,
     4	IBOUN(IDMN,JDMN),ICLO
      COMMON /WW/IUWS,IUDIR
C
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C
      PRINT *,'OPEN FILES'
      OPEN(20,FILE='OPTIONS.DAT',STATUS='OLD')
      OPEN(10,FILE='TWOD.DAT',STATUS='UNKNOWN')
      OPEN(11,FILE='ONELNS.DAT',STATUS='UNKNOWN')
      OPEN(13,FILE='ONED.DAT',STATUS='UNKNOWN')
      OPEN(15,FILE='WAVES.DAT',STATUS='UNKNOWN')
      OPEN(16,FILE='NXTSTR.BN',STATUS='UNKNOWN',FORM='UNFORMATTED')
      OPEN(93,FILE='SUMM.DAT',STATUS='UNKNOWN')
      OPEN(23,FILE='WARM.BN',STATUS='UNKNOWN',FORM='UNFORMATTED')
      OPEN(41,FILE='BND.ASC',STATUS='UNKNOWN')
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

C

      PRINT *,'BEGIN EXECUTION '
      XLAT0 = 0.0
      DLAT = 0.0
      IENDWN=0
C	   GRAVITY, (M/S)
      G=9.8
      PI=3.1415927
      PI8=PI/4.
      TWOPI=2.*PI
      RADDEG=360./TWOPI
      RADC=TWOPI/360.
      PI72=TWOPI/72.
      PI2=PI/2.
      KKNT=0
C
C

C	   RDOPT IS THE SUBROUTINE WHICH READS INPUT OPTIONS

      CALL RDOPT

C     OPEN OPTIONAL FILES AS DIRECTED IN RDOPT

      IF (NORD.NE.1) OPEN(21,FILE='WINDS.DAT',STATUS='OLD')

      PRINT *,'OPTIONS READ COMPLETE'
      IF (IBND.GE.2) OPEN(25,FILE='BNDIN.BN',STATUS='OLD')
      IF (NX.LE.IDMN.AND.NY.LE.JDMN) GO TO 7000
      PRINT 7001,IDMN,NX,JDMN,NY
 7001 FORMAT ('LIMITS ON IDMN OR JDMN EXCEEDED ***********',/,
     2 'IDMN = ',I4,' NX = ',I4,' JDMN = ',I4,' NY = ',I4)
      PRINT *,'BECAUSE OF POTENTIAL TO WRITE OUT OF PROPER BOUNDS'
      PRINT *,'COMPUTER RUN IS DISCONTINUED !!!!!!!!!!!!!!!!!!!!!'
      STOP
 7000 AINCR=TWOPI/NANG
      AINC=360./NANG
      FAC=AINC
      DO 7200 I=1,NX
      DO 7200 J=1,NY
      DEPNEW(I,J)=DEP(I,J)
 7200 CONTINUE
C

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      NF1=NFRQ-4
C	   SET PROPAGATION ANGLE REFERENCE ARRAY
  453 CONTINUE
      TINC=DT
C
C
C *************************************************************
C     D E E P	W A T E R    P R O P A G A T I O N
C
C
      PRINT *,'DEEP OPTION COEFFICIENTS BEING CALCULATED'

C
C     DZRO IS THE CICUMFERENCE OF THE EARTH IN METRES
C
      DZRO=40000000.
C
C	   LOOP TO ADJUST GRID SPACING BY ANGLE AND LATITUDE
C	   CGG = GROUP VELOCITY OF WAVE TRAIN, (M/S)
C	   DS = DISTANCE TRAVELED BY ONE WAVE TRAIN
C		IN ONE TIME STEP, (M)
      DO 777 K=1,NFRQ
      PRINT *,'K,FREQ(K) ',K,FREQ(K)
      CGG=0.5*G/(TWOPI*FREQ(K))
      DS=CGG*TINC
      DO 777 J=2,NY-1
      DO 777 IIA=1,NANG
      XX=ABS(COSR(IIA))
      YY=ABS(SINR(IIA))
      ALAT=((J-1)*DLAT + XLAT0)*RADC
      XFAC=COS(ALAT)
      DELX(J)=DL*XFAC
C
C     XDU AND YDU ARE THE TRIANGULAR DISTANCES TRAVELED
C     THEY CANNOT BE GREATER THAN ONE BECAUSE THE ENERGY
C     CANNOT TRAVEL MORE THAN ONE GRID SPACING IN ONE TIME STEP

      ADU(K,J,IIA)=DS*TWOPI/(DZRO*AINCR)  *ABS(SIN(ALAT)) * XX

      IF (ADU(K,J,IIA).GT.1.) ADU(K,J,IIA)=1.
      XDU(K,J,IIA)=XX*DS/(DL*XFAC)
      IF (XDU(K,J,IIA).GT.1.) XDU(K,J,IIA)=1.0
      YDU(K,IIA)=YY*DS/DL
      IF (YDU(K,IIA).GT.1.0) YDU(K,IIA)=1.0
      XD(K,J,IIA)=1.-XDU(K,J,IIA)
      YD(K,IIA)=1.-YDU(K,IIA)
      AD(K,J,IIA)=1.-ADU(K,J,IIA)
  777 CONTINUE
C     IUB AND JUA ARE COEFFICIENTS USED TO DETERMINE
C     THE UPSTREAM POINT LOCATIONS
      DO 776 J=1,NY
      ISIGNH=-1
      ALAT=(J-1)*DLAT + XLAT0
      IF (ALAT.GE.0.) ISIGNH=1
      DO 776 IIA=1,NANG
      X=COSR(IIA)
      IF (ABS(X).LT.0.001) GO TO 7761
      IUB(IIA)=(ABS(X)+0.001)/X
      IAUA(IIA,J)=IIA+IUB(IIA)*ISIGNH
      IF (IAUA(IIA,J).GT.NANG) IAUA(IIA,J)=1
      GO TO 7762
 7761 IUB(IIA)=0
      IAUA(IIA,J)=IIA
 7762 Y=SINR(IIA)
      IF (ABS(Y).LT.0.001) GO TO 7763
      JUA(IIA)=(ABS(Y)+0.001)/Y
      GO TO 776
 7763 JUA(IIA)=0
  776 CONTINUE
      PRINT * ,'COMPLETED COEFFICIENT CALCULATION '
C
C ***************************************************************

      NXM1=NX-1
      NYM1=NY-1
C
C	   LOOP TO DETERMINE THE EQUILIBRIUM RANGE
C	   WKZO = DEEP WATER WAVE NUMBER
C	   WKIJ = ARBITRARY DEPTH WAVE NUMBER
C	   ZJAC FROM SHALLOW-WAVER WAVES I: THEORY BY RESIO
C	   PAGE 273, REF. KITAIGORODSKI(1983) K(-5/2) EQUILIBRIUM
      PRINT *,'BEGIN 2 LOOP'
      DO 2 K=1,NFRQ
C
      FKA(K)=G/(TWOPI**3*FREQ(K)**4)
      WKZ0=(TWOPI*FREQ(K))**2 / G
      FFF=FREQ(K)
      DO 9 I=1,NX
      DO 9 J=1,NY
      ZJAC(I,J,K)=1.
      DDD=DEP(I,J)
C     IF (IBOUND(I,J).NE.1) GO TO 9
      IF (DDD.LT.1.0) DDD=1.
      IF (DDD.GT.200.) DDD=200.
      WKIJ=WKFNC(FFF,DDD)
      ZJAC(I,J,K)=(WKZ0/WKIJ)**2.5
      CCC=TWOPI*FFF/WKIJ
      CA(I,J,K)=CCC
      CGA(I,J,K)=CGF(FFF,DDD,CCC)
    9 CONTINUE
C
C  MATCH FREQUENCY RELATED CONSTANTS ON BOUNDARY TO
C    VALUES AT NEXT INTERIOR POINT
	  DO 10 J=1,NY
	  ZJAC(1,J,K)=ZJAC(2,J,K)
	  CA(1,J,K)=CA(2,J,K)
	  CGA(1,J,K)=CGA(2,J,K)
	  ZJAC(NX,J,K)=ZJAC(NX-1,J,K)
	  CA(NX,J,K)=CA(NX-1,J,K)
	  CGA(NX,J,K)=CGA(NX-1,J,K)
  10	  CONTINUE
	  DO 11 I=1,NX
	  ZJAC(I,1,K)=ZJAC(I,2,K)
	  CA(I,1,K)=CA(I,2,K)
	  CGA(I,1,K)=CGA(I,2,K)
	  ZJAC(I,NY,K)=ZJAC(I,NY-1,K)
	  CA(I,NY,K)=CA(I,NY-1,K)
	  CGA(I,NY,K)=CGA(I,NY-1,K)
  11	  CONTINUE
    2 CONTINUE


      IF (IDPRP.NE.1) CALL INITSPR


C
C	   LOOP TO INTIALIZE VARIABLES AND PARAMETERS
      DO 250 I=1,NX
      DO 250 J=1,NY
      WAVMAX(I,J)=0.
      TPPMAX(I,J)=0.
      WDDMAX(I,J)=0.
      FMA(I,J)=1.0
      HSIG(I,J)=0.
      SWANG(I,J)=0.
      AVANG(I,J)=0.
      KFRQA(I,J)=NFRQ+1
      DO 255 IB=1,NANG
      EN(I,J,IB)=0
      DO 255 K=1,NFRQ
      E(I,J,K,IB)=0
  255 CONTINUE
  250 CONTINUE
C




      PRINT *,'INITIALIZATION COMPLETE'
C
C     DEFINE PANG MATRIX (ANGLE DIFFERENCE MATRIX)
C
C
C	   LOOP TO DEFINE DIRECTIONAL SPREADING FUNCTION
C	   CALCULATES 5 DIFFERENT VALUES OF COS FUNCTION
C	   COS(DIF)**4 WAS USED FOR GREAT LAKES STUDY
C	   IDIR=1,72 REPRESENTS 5 DEGREE DIR BAND INTERVALS
      DO 400 NNN=1,5
      NN=NNN*2
      DO 400 IIA=1,NANG
      ANG1=(IIA-1)*AINCR
      SUM=0.
      DO 401 IDIR=1,72
      PANG(IIA,IDIR,NNN)=0.
      ANG2=(IDIR-1)*PI72
      DIF=ABS(ANG1-ANG2)
      IF (DIF.GT.PI) DIF=TWOPI-DIF
      IF (DIF.GT.0.9*PI2) GO TO 401
      XX=COS(DIF)**NN
C
      SUM=SUM+XX
      PANG(IIA,IDIR,NNN)=XX
  401 CONTINUE
      SUM=SUM*PI72
      DO 402 IDIR=1,72
      PANG(IIA,IDIR,NNN)=PANG(IIA,IDIR,NNN)/SUM
  402 CONTINUE
  400 CONTINUE
C	   CKSTR SUBROUTINE READS & WRITES RESTART INFO
      PRINT *, 'CALLING CKSTR'
      CALL CKSTR(1)
      PRINT *, 'RETURNING FROM CALL TO CKSTR'

C	   ONLY ACCESSED IF READING IN BOUNDARY DATA
      IF(IBND.GE.2)THEN
C	   INITB SUBROUTINE INITIALIZES BOUNDARY CONDITIONS
      PRINT *, 'CALLING INITB'
      CALL INITB
      PRINT *,'BOUNDARY INITIALIZATIONS COMPLETE'

C	   RDBND SUBROUTINE READS IN THE NEXT SET OF
C	   BOUNDARY INFO
      CALL RDBND (IDM)
C	   LOOP TO EVALUATE THE NEXT FREQUENCY &
C	   ZERO THE ENERGY ARRAY
      DO 11008 K=1,KNTBP
      I  =ICBI(K)
      J  =JCBI(K)
      FMALST(K)=FMA(I,J)
      DO 11009 IFRC=1,NFRQ
      DO 11009 IANG=1,NANG
       ELST(K,IFRC,IANG)=0
11009 CONTINUE
11008 CONTINUE
      PRINT *,'INITIAL DATES BNDS = ',IDM
      ENDIF
C
C	  ONLY ACCESSED IF WIND INPUT IS READ IN
      IF (NORD.EQ.1) GO TO 1111
C	   RDWIN SUBROUTINE READS IN THE NEXT SET OF WINDS
      CALL RDWIN (ID,IENDWN)
      PRINT *,'INITIAL DATES WIND = ',ID
 1111 TIMTOT=0.
C
C
C     PRINT *,'BEGIN TIME STEPS'
C
C
C	   LOOP OVER THE MAXIMUM NUMBER OF WIND INPUTS TO PROCESS
      DO 1000 IHR=1,MXHR
      PRINT *,'INSIDE DO 1000 LOOP, IHR =',IHR
C	   ONLY ACCESSED IF WIND INPUT IS READ IN
      IF (NORD.EQ.1) GO TO 1100
      IDL=ID
      DO 1005 I=1,NX
      DO 1005 J=1,NY
      IWSL(I,J)=IWSN(I,J)
      IWDL(I,J)=IWDN(I,J)
 1005 CONTINUE

C     PRINT *,'BEFORE CALL TO RDWIN'
C	   SUBROUTINE RDWIN READS IN THE NEXT SET OF WINDS
C	   THIS SUBROUTINE IS CALLED TWICE SINCE TWO VALUES
C	   ARE NEEDED FOR THE INTERPOLATION SCHEME
      PRINT *,'CALLING RDWIN'
      CALL RDWIN(ID,IENDWN)
      PRINT *,'AFTER RDWIN'

C     INITIALIZING SHALLOW-WATER PROPAGATION MATRICES

      IF (IUPDAT.NE.1) GO TO 7300

      DO 7301 K=1,NFRQ
C
      FKA(K)=G/(TWOPI**3*FREQ(K)**4)
      WKZ0=(TWOPI*FREQ(K))**2 / G
      FFF=FREQ(K)
      DO 7309 I=1,NX
      DO 7309 J=1,NY
      ZJAC(I,J,K)=1.
      DEP(I,J)=DEPNEW(I,J)
      DDD=DEP(I,J)

C     IF (IBOUND(I,J).NE.1) GO TO 7309

      IF (DDD.LT.1.0) DDD=1.
      IF (DDD.GT.200.) DDD=200.
      WKIJ=WKFNC(FFF,DDD)
      ZJAC(I,J,K)=(WKZ0/WKIJ)**2.5
      CCC=TWOPI*FFF/WKIJ
      CA(I,J,K)=CCC
      CGA(I,J,K)=CGF(FFF,DDD,CCC)

 7309 CONTINUE
 7301 CONTINUE


      IF (IDPRP.NE.1) CALL INITSPR

 7300 CONTINUE


C	   IF END OF FILE WAS READ IN RDWIN, END RUN.
      PRINT *,'WINDS READ ID = ',ID

      IF (IENDWN.EQ.1) GO TO 9999

C	   THE FOLLOWING IS NOT PERFORMED FOR BOUNDARY DATA
C	   NOT READ IN
      IF (IBND.LT.2) GO TO 1152
C	   LOOP TO ZERO THE LAST FREQUENCY AND ENERGY VALUES
      DO 1008 K=1,KNTBP
      FMALST(K)=FMANXT(K)
      DO 1009 IFRC=1,NFRQ
      DO 1009 IANG=1,NANG
      ELST(K,IFRC,IANG)=ENXT(K,IFRC,IANG)
 1009 CONTINUE
 1008 CONTINUE
C	   SUBROUTINE RDBND READS IN THE NEXT SET OF BOUNDARY DATA
      CALL RDBND(IDM)
      PRINT *,'BOUNDARYS READ IDM = ',IDM
      IF (NORD.EQ.0) GO TO 1152
C	   BEGIN TIME STEP LOOP FOR WIND INPUT
C	   THE NEXT 10 LINES REFER TO TEST CASES ONLY
 1100 DO 1153 I=1,NX
      DO 1153 J=1,NY
      IWSL(I,J)=IUWS
      IWSN(I,J)=IUWS
      IWDL(I,J)=IUDIR
      IWDN(I,J)=IUDIR
 1153 CONTINUE
      IDN=IHR
      PRINT *,'WINDS SET IN TEST SITUATION AT SPEED,DIRECTION ',IUWS,
     2	IUDIR
 1152 CONTINUE
      PRINT *, 'ID= ', ID, ' IDFRST= ', IDFRST
      IF (ID.LT.IDFRST) GO TO 1000
      IF (ID.GT.IDLAST) GO TO 9999

C	   LOOP OVER THE NUMBER OF TIME STEPS BETWEEN WIND INPUT
      DO 1001 KTIME=1,KTIMES
      TIMTOT=TIMTOT+TINC
      PRINT *,'INSIDE DO 1001 LOOP, KTIME =',KTIME,' IHR = ',IHR
      KKNT=KKNT+1
C	   INTERPOLATE FOR CURRENT WINDS
C	   WHERE WL & WN ARE MULTIPLIERS FOR INTERPOLATION
C	   THEY ARE CALCULATED IN SUBROUTINE RDOPT
      PRINT *, 'GOING INTO 1020 LOOP'
      DO 1020 I=1,NX
      DO 1020 J=1,NY
C     IWS(I,J)=WL(KTIME)*IWSL(I,J)+WN(KTIME)*IWSN(I,J)
C     IF (ABS(IWDN(I,J)-IWDL(I,J)).LT.36) GO TO 1025
C     IF (IWDL(I,J).LT.36) IWDL(I,J)=IWDL(I,J)+72
C     IF (IWDN(I,J).LT.36) IWDN(I,J)=IWDN(I,J)+72
C1025 IWD(I,J)=0.5+WL(KTIME)*IWDL(I,J)+WN(KTIME)*IWDN(I,J)
C DO COMPONENT INTERPOLATION ON WINDS
	WDLL=((REAL(IWDL(I,J))-1.)*5.0)*RADC
	WXCL=REAL(IWSL(I,J))*COS(WDLL)
	WYCL=REAL(IWSL(I,J))*SIN(WDLL)
	WDNN=((REAL(IWDN(I,J))-1.)*5.0)*RADC
	WXCN=REAL(IWSN(I,J))*COS(WDNN)
	WYCN=REAL(IWSN(I,J))*SIN(WDNN)
	WXX=WL(KTIME)*WXCL+WN(KTIME)*WXCN
	WYY=WL(KTIME)*WYCL+WN(KTIME)*WYCN
	IWS(I,J)=SQRT(WXX**2+WYY**2)+0.5
	WDDD=ATAN2(WYY,WXX+0.0000001)
	IWD(I,J)=INT(WDDD*RADDEG/5.0+0.5)
      IF (IWD(I,J).LT.1) IWD(I,J)=IWD(I,J)+72
      IF (IWD(I,J).GT.72) IWD(I,J)=IWD(I,J)-72
 1020 CONTINUE

      PRINT *,'WINDS UPDATED'

      IF (IBND.LT.2) GO TO 1007
C	   INTERPOLATE FOR CURRENT BOUNDARY CONDITIONS
C      FMA AND ENERGY SPECTRA ONLY
C
      DO 1024 LNTBP=1,KNTBP
      IPT=ICBI(LNTBP)
      JPT=JCBI(LNTBP)
C
      FMA(IPT,JPT)=    WL(KTIME)*FMALST(LNTBP) + WN(KTIME)*FMANXT(LNTBP)
      DO 1125 IFRC=1,NFRQ
      DO 1125 IANG=1,NANG
      E(IPT,JPT,IFRC,IANG)=WL(KTIME)*ELST(LNTBP,IFRC,IANG)+
     2		       WN(KTIME)*ENXT(LNTBP,IFRC,IANG)
 1125 CONTINUE

 1024 CONTINUE
      PRINT *,'BOUNDARYS UPDATED'

 1007	   CONTINUE
C     PRINT *,'KTIME = ',KTIME
   93 FORMAT (21I6)
C
C
C     SET SIDE BOUNDARIES FOR 1-D TESTS
C

C	   SUBROUTINE SETSID SETS BOUNDARIES FOR A SPECIAL
C	   TEST CASE
      IF (NORD.EQ.1) CALL SETSID(NX,NY,NFRQ,NANG,ISPEC)
C
C
C	   BEGIN FREQUENCY LOOP
C
C
C     SET WRAP-AROUND GLOBAL BOUNDARIES IF IWRAP = 1
C
C	   FOR SMOOTH TRANSITION FOR GLOBAL WRAP-AROUND
      IF (IWRAP.NE.1) GO TO 2020
      DO 2019 J=2,NYM1
      DO 2018 K=1,NFRQ
      DO 2018 IANG=1,NANG
      E(1,J,K,IANG)=E(NXM1,J,K,IANG)
      E(NX,J,K,IANG)=E(2,J,K,IANG)
      FMA(1,J)=FMA(NXM1,J)
      FMA(NX,J)=FMA(2,J)
 2018 CONTINUE
 2019 CONTINUE
C
C
 2020 DO 2000 IFRC=1,NFRQ
C     PROPAGATION ROUTINE
C	   SUBROUTINE DPROP UPDATES ENERGY AT PREVIOUS TIME STEP

      IF (IDPRP.EQ.1) CALL DPROP(IFRC)
      IF (IDPRP.NE.1) CALL SPROP(IFRC)


 2000 CONTINUE
      PRINT *,'AFTER PROP IN MAIN'
C
C	   CALCULATE SOURCE TERMS
C
C	    SUBROUTINE FLXSTD CALLS SUBROUTINE SORPRP AND CALCULATES
C	    THE SOURCE TERMS AND UPDATES THE ENERGY VALUES
       CALL FLXSTD
       PRINT *,'AFTER FLXSTD IN MAIN'

      IF(KKNT.NE.1)GO TO 4451

      IF(IBND.EQ.1.OR.IBND.EQ.3)THEN
      PRINT *,' IBND=',IBND,' KKNT=',KKNT,' IHR=',IHR,' KTIME=',KTIME
C	   SUBROUTINE WBOUND WRITES THE OUTPUT BOUNDARY DATA
      CALL WBOUND(IDL)
      ENDIF

 4451 CONTINUE


C     WRITE OPTIONAL FORECAST RESTART INFORMATION
      IF (ABS(TIMTOT/3600.-NFCST).LT.0.5) CALL CKSTR(2)


C
C	   IF NO OUTPUT TO BE PRINTED, JUMP TO END OF TIME STEP LOOP.
      IF (MOD(KKNT,NTMS).NE.0) GO TO 1001

C	   UPDATE DATE-TIME CODE FOR OUTPUT.
C     PRINT *,' INSIDE OUTPUT LOOP'
      IF (NORD.NE.1) CALL NDATE(IDL,KTIME,IDN)
      IF (NORD.EQ.1) IDL=IHR
C	   CALL SUBROUTINE TO WRITE TWO-DIMENSIONAL SPECTRA FILE AND
C	     ONE-LINE SUMMARIES FILE
      IF (NORD.EQ.1) IDN=IHR
      PRINT *,' CALL OUTP'
C	   SUBROUTINE OUTP WRITES THE 1-LINE SUMMARY & 2D SPECTRA FILES
      CALL OUTP(IDN)
      PRINT *,'RETURNING FROM OUTP'
C	   END TIME STEPS
      IF (IDN.LT.IDBND1) GO TO 1000
      IF (IDN.GT.IDBND2) GO TO 1000
      IF (IBND.EQ.1.OR.IBND.EQ.3) CALL WBOUND(IDN)

 1001 CONTINUE
 1000 CONTINUE
 9999 CONTINUE
C	   SUBROUTINE MXXOUT WRITES THE MAXIMUM WAVE HEIGHT FILE
C     CALL MXXOUT
C	   SUBROUTINE CKSTR READS & WRITES THE RESTART DATA
      PRINT *, 'WRITING WARM START DATA'
      CALL CKSTR(2)

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      CLOSE(10)
      CLOSE(11)
      CLOSE(13)
      CLOSE(15)
      CLOSE(16)
      CLOSE(17)
      CLOSE(23)
      CLOSE(25)

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      STOP
      END
C
      SUBROUTINE MXXOUT
      INCLUDE 'PARAMS'
C      PARAMETER (IDMN=35,JDMN=51,IF=15,IA=16,NOBPTS=11,NBN=11)
C      PARAMETER(IDMN=78,JDMN=56)
      COMMON /MAXOUT/ WAVMAX(IDMN,JDMN),TPPMAX(IDMN,JDMN),
     1		      WDDMAX(IDMN,JDMN)
      OPEN(44,FILE='MAXFILE.DAT',STATUS='UNKNOWN')
      DO 10 J=JDMN,1,-1
   10 WRITE(44,20)(WAVMAX(I,J),I=1,IDMN)
   20 FORMAT(61F6.1)
      DO 30 J=JDMN,1,-1
   30 WRITE(44,20)(TPPMAX(I,J),I=1,IDMN)
      DO 40 J=JDMN,1,-1
   40 WRITE(44,50)(WDDMAX(I,J),I=1,IDMN)
   50 FORMAT(61F6.0)
      CLOSE(44)
      RETURN
      END
C
C
C
	SUBROUTINE RDOPT
C---------------------------------------------------------------------
C  THIS ROUTINE READS IN PROGRAM OPTION AND CONFIGURES MODEL RUN
C
C    PARAMETER STATEMENT
C	 IDMN	= MAXIMUM NUMBER OF GRID COLUMNS
C	 JDMN	= MAXIMUM NUMBER OF GRID ROWS
C	 IF	= MAXIMUM NUMBER OF FREQUENCY BANDS
C	 IA	= MAXIMUM NUMBER OF ANGLE BANDS
C	 NOBPTS = MAXIMUM NUMBER OF OUTPUT BOUNDARY POINTS
C	 NBN	= MAXIMUM NUMBER OF INPUT BOUNDARY POINTS
C    DATA READ IN
C----- OPTIONS RECORD 1
C	 NX	= NUMBER OF COLUMNS IN GRID
C	 NY	= NUMBER OF ROWS IN GRID
C	 NANG	= NUMBER OF DIRECTIONAL ANGLE BANDS
C	 NFRQ	= NUMBER OF FREQUENCY BANDS
C	 DL	= DISTANCE BETWEEN GRID POINTS (KM)
C	 DT	= TIME STEP (SEC)
C	 MSTA	= NUMBER OF SPECIAL OUTPUT LOCATIONS
C	 IFTYP	= DETERMINES HOW FREQUENCIES ARE OBTAINED
C		  0 - READ IN FREQUENCIES (SMALLEST TO LARGEST) (HZ)
C		  1 - READ IN FIRST FREQUENCY AND GEOMETRIC FACTOR FOR
C		      COMPUTED FREQUENCIES
C	 INPLEV = DETERMINES ELEVATION OF WINDS
C		  0 - WINDS AT 10 METRES
C		  1 - WINDS AT 20 METRES
C	 NSSPTS = NUMBER OF SUBSCALE FEATURES ( ALWAYS SET TO 0 )
C     NFCST = 0 NO FORCAST RESTART INFORMATION SAVED
C           = N>0  FORCAST RESTART INFORMATION SAVED AT HOUR N
C----- OPTIONS RECORD 2
C	 NSTR	= DETERMINES TYPE OF RUN MADE
C		  0 - COLD START, NO SAVE
C		  1 - COLD START, SAVE
C		  2 - WARM START, NO SAVE
C		  3 - WARM START, SAVE
C	 NORD	= DETERMINES TYPE OF INPUT FOR WINDS
C		  0 - WINDS READ IN
C		  1 - CONSTANT WINDS
C	 NTMS	= NUMBER OF TIME STEPS BETWEEN OUTPUT RESULTS
C	 NHR	= NUMBER OF HOURS BETWEEN INPUT WINDS
C	 MXHR	= MAXIMUM NUMBER OF WIND INPUTS
C	 IBND	= DETERMINES BOUNDARY DATA INPUT AND OUTPUT
C		  0 - NONE READ OR WRITTEN
C		  1 - DATA WRITTEN NOT READ
C		  2 - DATA READ NOT WRITTEN
C		  3 - DATA READ AND WRITTEN
C	 IWRAP	= DETERMINES IF DATA IS WRAPPED AROUND
C		  0 - NO WRAP AROUND
C		  1 - WRAP AROUND
C	 ISTDEP = DETERMINES HOW DEPTHS ARE INPUT
C		  0 - DEPTH FIELD READ IN
C		  1 - CONSTANT DEPTH OF 999 (M)
C----- OPTIONS RECORD 3
C	 DLAT	= LATITUDE INCREMENT OF GRID (DEG)
C	 XLAT0	= SOUTHERNMOST LATITUDE (DEG) ( - FOR SOUTHERN HEMISPHERE)
C	 ICURV	= CURVATURE EFFECTS INCLUDED IN PROPAGATION
C		  0 - DO NOT INCLUDE
C		  1 - INCLUDE
C	 IDPRP = DEEP OR SHALLOW WATER PROP. INDICATOR
C		 0  CALLS INITSPR FOR SHALLOW WATER
C		 1  ONLY DEEPWATER PROP.
C	 IUPDAT = INDICATOR FOR DEPTH CHANGES
C		 0 DOES NOT UPDATE CG'S AND K'S
C		 1 INCORPORATES DEPTH CHANGES
C----- OPTIONS RECORD 4
C---------- IF IFTYP = 0
C		FREQ(IFRC),IFRC =1,NFRQ -- FREQUENCIES FROM 1ST TO MAXIMUM
C			NFRQ = NUMBER OF FREQUENCY BANDS
C			FREQ = FREQUENCIES
C---------- IF IFTYP = 1
C		FZRO = INITIAL FREQUENCY
C		PWR  = FREQUENCY EXPONENTIAL
C----- OPTIONS RECORD 5
C---------- IF MSTA > 0
C		NSTA(ISTA) = STATION NUMBER OF SPECIAL OUTPUT LOCATIONS
C		IOUT(ISTA) = I LOCATION OF SPECIAL OUTPUT POINT
C		JOUT(ISTA) = J LOCATION OF SPECIAL OUTPUT POINT
C ISTA = 1,MSTA
C----- OPTIONS RECORD 6
C      DISCONTINUED IN WAVAD 10/01/93
C---------- IF NSSPTS NOT EQUAL TO 0 NEED SUBSCALE POINTS
C		INSSPS
C		JNSSPS
C		BLKX
C		BLKY
C----- OPTIONS RECORD 7
C---------- IF ISTDEP = 0
C		DEP(I,J),I=1,NX
C----- OPTIONS RECORD 8
C---------- IF ISTDEP = 0
C		IBOUND(I,J),I=1,NX
C			0 - LAND POINT
C			1 - WATER POINT
C			2 - EXTRAPOLATED BOUNDARY POINT AT LEADING EDGE
C			3 - INPUT BOUNDARY POINT
C			4 - INTERPOLATED BOUNDARY POINT
C			5 - EXTRAPOLATED BOUNDARY POINT AT TRAILING EDGE
C----- OPTIONS RECORD 9
C---------- IF IBND = 0
C		RETURN TO MAIN SUBROUTINE
C---------- IF IBND = 1
C		NBPO		 = NUMBER OF BOUNDARY POINTS TO BE OUTPUT
C		IBPO(K),K=1,NBPO = ORDERED I-LOCATIONS OF BOUNDARY OUTPUT POINTS
C		JBPO(K),K=1,NBPO = ORDERED J-LOCATIONS OF BOUNDARY OUTPUT POINTS
C---------- IF IBND = 2
C		ILS = STARTING I COORDINATE OF BOUNDARIES TO BE READ IN
C		JLS = STARTING J COORDINATE OF BOUNDARIES TO BE READ IN
C---------- IF IBND = 3
C		NBPO		 = NUMBER OF BOUNDARY POINTS TO BE OUTPUT
C		IBPO(K),K=1,NBPO = ORDERED I-LOCATIONS OF BOUNDARY OUTPUT POINTS
C		JBPO(K),K=1,NBPO = ORDERED J-LOCATIONS OF BOUNDARY OUTPUT POINTS
C		ILS = STARTING I COORDINATE OF BOUNDARIES TO BE READ IN
C		JLS = STARTING J COORDINATE OF BOUNDARIES TO BE READ IN
C----- OPTIONS RECORD 10
C---------- IF NORD.NE.1  NO READ EXECUTED (COMPLETE WIND FIELDS ARE INPUT)
C---------- IF NORD = 1
C      IUWS   CONSTANT WIND SPEED FOR TEST RUN  (M/SEC)
C      IUDIR  CONSTANT WIND DIRECTION DEGREES FROM AZIMUTH
C
C     COMMON STATEMENTS USED TO EXPORT DATA
C	IOP,OUTTTS,OPT
C
C
      INCLUDE 'PARAMS'
C     PARAMETER (IDMN=35,JDMN=51,IF=15,IA=16,NOBPTS=11,NBN=11)
      COMMON /IOP/ ICURV,IDPRP,IUPDAT
      COMMON /IFCST/NFCST
      COMMON /OUTTTS/ IKOUTL,ISTWAV
      COMMON /PMOD/ XD(IF,JDMN,IA),XDU(IF,JDMN,IA),YD(IF,IA),
     1		    YDU(IF,IA)
      COMMON /DPR/ IUB(IA),JUA(IA),IAUA(IA,JDMN),AD(IF,JDMN,IA),
     1		    ADU(IF,JDMN,IA)
      COMMON /STREND/ IDFRST,IDLAST,IDBND1,IDBND2
      COMMON /A3/ G,PI,XN,DL,NANG,NFRQ,NX,NY,TWOPI,FAC,DT,
     2		  RADDEG
      COMMON /OPT/ MSTA,NSTR,NORD,JPRC,NTMS,INDELF,NHR,
     2		   IOUT(70),JOUT(70),KTIMES,WL(725),WN(725),
     3 IG21,IR21,MXHR,IDEEP,DLAT,XLAT0,DZRO,NSTA(70),ALT(70),ALON(70)
      COMMON /WW/IUWS,IUDIR

C***  MOD 08/10/88
      COMMON /SSPT/ TRANSX(IDMN,JDMN),TRANSY(IDMN,JDMN)
C***

      COMMON /O1/ FREQ(IF),SINR(IA),COSR(IA),ALFA0(IA),ALFAR0(IA),
     2		   DELF(IF),AINC,AINCR
      COMMON /O5/ IBOUND(IDMN,JDMN)
      COMMON /BD/ IBND,NBPO,IBPO(NOBPTS),JBPO(NOBPTS),NBPI,KNTBP,
     2	     IBPI(NBN),JBPI(NBN)
      COMMON /MM/ INPLEV,IWRAP
      COMMON /HBND/  IBEGA(NBN),JBEGA(NBN),IENDA(NBN),JENDA(NBN),
     1 ITB(NBN,12),JTB(NBN,12),NLA(NBN),WTB(NBN,12),
     2	WTA(NBN,12),NSEG,NINT,ILS,JLS,ICBI(NBN),JCBI(NBN),
     3	IF1,JF1,IFN,JFN,NEXF1,NINPT,NINTR,NEXFN,NOF1,NOFN,
     4	IBOUN(IDMN,JDMN),ICLO
      COMMON /DEPTH/ DEP(IDMN,JDMN)
	  COMMON /DEPOR/  DEPORG(IDMN,JDMN)
      CHARACTER*4 TITL(2)
      DATA TITL /'COLD','WARM'/


      PRINT *,'READING OPTIONS RECORD 1'
      READ (20, * ) NX,NY,NANG,NFRQ,DL,DT,MSTA,IFTYP,NFCST
      PRINT *, NX,NY,NANG,NFRQ,DL,DT,MSTA,IFTYP,NFCST
      INPLEV=0
      NSSPTS=0

C	   CONVERT DL TO METERS
      DL=DL*1000.0
C  SET THE MAXIMUM NUMBER OF INPUT COLUMN FORMAT (28)
      IG21=NX/21
      IR21=MOD(NX,21)

      PRINT *,'READING OPTIONS RECORD 2'
      READ (20, * ) NSTR,NORD,NTMS,NHR,MXHR,IBND,IWRAP,ISTDEP
      PRINT *, NSTR,NORD,NTMS,NHR,MXHR,IBND,IWRAP,ISTDEP
      IKOUTL=0

      PRINT *,'READING OPTIONS RECORD 3'
      READ(20, * ) DLAT,XLAT0,ICURV,IDPRP
      PRINT *, DLAT,XLAT0,ICURV,IDPRP
      IUPDAT=0

      PRINT *,'READING OPTIONS RECORD 4'
C  IF IFTYP=1 GO TO 666 AND READ IN 1ST FREQUENCY AND THE
C  GEOMETRIC FACTOR TO COMPUTE ADDITIONAL FREQUENCIES
      IF (IFTYP.EQ.1) GO TO 666
C  READ IN FREQUENCIES FROM THE INPUT FILE
      READ (20, * ) (FREQ(IFRC),IFRC=1,NFRQ)
      PRINT * ,'FREQUENCIES BEING MODELED : '
      DO 7788 K=1,NFRQ
      PRINT *,'FREQUENCY NUMBER ',K,' AT ',FREQ(K),' HZ'
 7788 CONTINUE
C  SKIP OVER FREQUENCY COMPUTATIONS
      GO TO 667
  666 READ (20, * ) FZRO,PWR
C  COMPUTE FREQUENCIES
      DO 668 I=1,NFRQ
      FREQ(I)=FZRO*PWR**(I-1)
  668 CONTINUE
C  DETERMINE CENTER FREQUENCY AT EACH FREQUENCY INCREMENT
  667 F1=0.5*(FREQ(2)-FREQ(1))
      DO 120 IFRC=1,NFRQ-1
      F2=(FREQ(IFRC+1)-FREQ(IFRC))*0.5
      DELF(IFRC)=F1+F2
  120 F1=F2
C  SET DELF(J) TO FREQUENCY INCREMENT FOR EACH BAND
      DELF(NFRQ)=2.0*F2
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C  PRINT OUT INPUT INFORMATION
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      PRINT *
      PRINT 901,NY
      PRINT 902,NX
      PRINT 903,DL
      IF (NSTR.LE.1) PRINT 906,TITL(1)
      IF (NSTR.GE.2) PRINT 906,TITL(2)
      IF (NSTR.EQ.0.OR.NSTR.EQ.2) PRINT 909
      IF (NSTR.EQ.1.OR.NSTR.EQ.3) PRINT 910
      PRINT 911,NTMS
      PRINT 914,NFRQ
      PRINT 915,(FREQ(IFRC),IFRC=1,NFRQ)
      PRINT 916,NFRQ
      IF (INDELF.EQ.1) PRINT 912
      PRINT 915,(DELF(IFRC),IFRC=1,NFRQ)

      PRINT *,'READING OPTIONS RECORD 5'
C	   READ GRID COORDINATES FOR SPECIAL OUTPUT LOCATIONS
      IF(MSTA.GT.0)
     1READ(20,*)(NSTA(ISTA),IOUT(ISTA),JOUT(ISTA),ISTA=1,MSTA)
      PRINT 918,MSTA
      DO 190 ISTA=1,MSTA
      PRINT 919,ISTA,IOUT(ISTA),JOUT(ISTA)
  190 CONTINUE
C
C      PRINT *,'READING OPTIONS RECORD 6'
C      READ(20,*)IDFRST,IDLAST,IDBND1,IDBND2
C      PRINT *, IDFRST,IDLAST,IDBND1,IDBND2

      IDFRST=       0
      IDLAST=99999999
      IDBND1=       0
      IDBND2=99999999

C	READ SUB-SCALE POINTS
      DO 500 I=1,NX
      DO 500 J=1,NY
      TRANSX(I,J)=1.0
      TRANSY(I,J)=1.0
  500 CONTINUE
      IF(NSSPTS.NE.0)THEN
	DO 510 I=1,NSSPTS
	READ(20,*)INSSPS,JNSSPS,BLKX,BLKY
	PRINT *,INSSPS,JNSSPS,BLKX,BLKY
	II=INSSPS
	JJ=JNSSPS
	TRANSX(II,JJ)=1.0-BLKX
	TRANSY(II,JJ)=1.0-BLKY
  510	CONTINUE
      ENDIF

      PRINT *,'READING OPTIONS RECORD 7:  DEPTHS'

      IF(ISTDEP.EQ.0)THEN
C  READ IN WATER DEPTHS
	 DO 400 J=NY,1,-1
  400	 READ(20,401)(DEP(I,J),I=1,NX)
  401	 FORMAT(15F5.0)
	 PRINT *,'DEPTH READ COMPLETE'
	 DO 404 J=NY,1,-1
  404	 PRINT 403,(IFIX(DEP(I,J)+0.5),I=1,NX)
  403	 FORMAT(20I5)
      ELSE
C  DEPTHS NOT READ IN, HAVE CONSTANT DEPTH OF 999 METERS FOR ENTIRE GRID
	 DO 402 I=1,NX
	 DO 402 J=1,NY
  402	 DEP(I,J)=999.
	 PRINT *,'DEPTH SET TO 999. FOR ENTIRE GRID'
      ENDIF
C SAVE ORIGINAL DEPTH TO USE WITH SURGE IN OUTP
	DO 705 I=1,NX
	DO 705 J=1,NY
	DEPORG(I,J)=DEP(I,J)
 705	CONTINUE


C	   DEFINE NUMBER OF TIME STEPS BETWEEN WIND INPUTS
      SECBI=NHR*3600.
      KTIMES=SECBI/DT+0.5
      XTIMES=KTIMES
      PRINT 922,NHR,DT,KTIMES
      PRINT 956,MXHR
C	   CALCULATE MULTIPLIERS FOR INTERPOLATION THROUGH TIME
      DO 160 I=1,KTIMES
      WL(I)=(XTIMES-I+1)/XTIMES
      WN(I)=(I-1)/XTIMES
  160 CONTINUE


C	   DEFINE DIRECTIONS IN DEGREES AND RADIANS
      AINC=360./NANG
      AINCR=TWOPI/NANG
      DO 170 IANG=1,NANG
      ALFA0(IANG)=(IANG-1)*AINC
      ALFAR0(IANG)=(IANG-1)*AINCR
  170 CONTINUE


C	   DEFINE SINE AND COSINE OF EACH ANGLE
      DO 180 IANG=1,NANG
      SINR(IANG)=SIN(ALFAR0(IANG))
      COSR(IANG)=COS(ALFAR0(IANG))
  180 CONTINUE
      PRINT 9876,NX,NY
 9876 FORMAT (1X,'CHECK ON VALUES NX,NY ',2I4)

      PRINT *,'READING OPTIONS RECORD 8'


	 PRINT *,'READING LAND-SEA MATRIX'
C-------------------------------------------------
C  READ IN BOUNDARY AND BATHEMETRY
C	GRID ORIENTATION IS:
C	   1,1 IN LOWER CORNER
C	   1,NY IS UPPER LEFT CORNER
C	TO READ TOP-DOWN LOOP OVER J=NY,1,-1.
C-------------------------------------------------
	       DO 295 J=NY,1,-1
	       READ (20,957) (IBOUND(I,J),I=1,NX)
	       PRINT 959,(IBOUND(I,J),I=1,NX)
	       IF (ISTDEP.EQ.0) GO TO 295
	       DO 296 I=1,NX
	       DEP(I,J)=DEP(I,J)*IBOUND(I,J)
  296	       CONTINUE
 295	       CONTINUE

      PRINT *,'VERIFYING DEPTH AND LAND-SEA GRIDS'
C  TEST LAND-SEA MATRIX TO MAKE SURE IT MAKES SENSE
      DO 988 I=1,NX
      DO 988 J=1,NY
      IF(IBOUND(I,J).EQ.0)GO TO 988
      IF(DEP(I,J).GT.0.01)GO TO 988
      PRINT *,'RUN ABORT***************'
      PRINT *,' DEPTH EQUAL TO ZERO ON A COMPUTATIONAL GRID POINT'
      PRINT *,' I = ',I,'    J = ',J
      STOP
  988 CONTINUE


C     CHECKING SELECTED SPECIAL OUTPUT POINT(S) TO SEE IF ANY FALL ON LAND

      DO 333 M=1,MSTA
      I=IOUT(M)
      J=JOUT(M)
      IF (IBOUND(I,J).EQ.0) PRINT *,
     2'!!! WARNING SPECIAL OUTPUT POINT AT I= ',I,' AND J = ',J,
     3 'IS ON LAND'
  333 CONTINUE

      PRINT *,'READING OPTIONS RECORD 9'
  957 FORMAT (81I1)
  959 FORMAT (1X,80I1)
      PRINT *,'IBND = ',IBND
      IF (IBND.EQ.0) GO TO 1080
      IF (IBND.EQ.2) GO TO 895
      PRINT *,'TRYING TO READ NBPO'
      READ (20, * ) NBPO
      READ (20, * ) (IBPO(K),K=1,NBPO)
      READ (20, * ) (JBPO(K),K=1,NBPO)
  307 FORMAT (20I3)
      IF (IBND.EQ.1) RETURN
  895 PRINT *,'READ STARTING (I,J) FOR BOUNDARY INPUT'
      READ (20, * ) ILS,JLS
      PRINT *,'ILS = ',ILS,'    JLS = ',JLS



 1080 IF (NORD.EQ.0) GO TO 1347
C
C     IUWS = WINDSPEED FOR TEST CASE (KNOTS)
C     IUDIR = WIND DIRECTION FOR TEST CASE
      PRINT *,'READING OPTIONS RECORD 10:  TEST RUN '
      READ (20, * ) IUWS,IUDIR
      PRINT *,'************  IDEALIZED TEST CASE  ***************'
      PRINT 1386,IUWS,IUDIR
 1386 FORMAT ('WIND SPEED IN KNOTS = ',I4,/,'WIND DIRECTION IS = ',
     2	I5)
C     CONVERT TO INTERNAL VECTOR FORM
      IUDIR=270-IUDIR
      IF (IUDIR.LT.0) IUDIR=IUDIR+360
C     CONVERT DIRECTIONS TO INTERNAL 5-DEGREE DIRECTION CODE
      IUDIR=IUDIR/5 + 1.5
      IF (IUDIR.LT.0) IUDIR=IUDIR+72
      IF (IUDIR.GT.72) IUDIR=IUDIR-72
      PRINT 1387
 1387 FORMAT ('WIND CONDITIONS CONSTANT FOR TEST')
 1347 PRINT *,'END SUBROUTINE RDOPT'
      RETURN

  901 FORMAT (5X,'NUMBER OF ROWS IN GRID =',I5)
  902 FORMAT (5X,'NUMBER OF COLUMNS IN GRID =',I5)
  903 FORMAT (5X,'DISTANCE IN METERS BETWEEN GRID POINTS =',F15.5)
  906 FORMAT (5X,'THIS RUN HAS A ',A4,' START.')
  909 FORMAT (5X,'THIS RUN WILL NOT SAVE DATA FOR WARM START',
     2	       ' OF A CONSECUTIVE RUN.')
  910 FORMAT (5X,'DATA FOR WARM START OF NEXT RUN WILL BE SAVED.')
  911 FORMAT (5X,'DATA WILL BE WRITTEN FOR SPECIAL OUTPUT LOCATIONS',
     2	       ' EVERY',I3,' TIME STEPS.')
  912 FORMAT (5X,'FREQUENCY BAND WIDTHS WERE READ IN.')
  914 FORMAT (5X,'THERE ARE',I3,' FREQUENCIES:')
  915 FORMAT (5X,10F7.4)
  916 FORMAT (5X,'THE',I3,' FREQUENCY BAND WIDTHS ARE:')
  918 FORMAT (5X,'THERE ARE',I3,' SPECIAL OUTPUT LOCATIONS:')
  919 FORMAT (15X,'STATION',I3,'  AT I=',I2,', J=',I2)
  922 FORMAT (5X,'THERE ARE',I3,' HOURS BETWEEN WIND INPUTS.',
     2	      /,7X,'TIME STEP IN SECONDS =',F10.2,
     3	      /,7X,'THERE ARE',I4,' TIME STEPS BETWEEN WIND INPUTS.')
  923 FORMAT (7X,10F7.5)
  929 FORMAT (5X,'WEIGHTINGS ON LAST WINDS:')
  930 FORMAT (5X,'WEIGHTINGS ON NEXT WINDS:')
  925 FORMAT (5X,'THERE ARE',I4,' DIRECTION BANDS.',
     2	      /,7X,'EACH DIRECTION BAND WIDTH IS',F6.1,' DEGREES.')
  926 FORMAT (7X,10F7.1)
  927 FORMAT (10X,9F7.4)
  956 FORMAT (5X,'MAXIMUM NUMBER OF WIND INPUTS READ IN WILL BE ',I4)
      END
C
C
	SUBROUTINE RDWIN(ID,IENDWN)
C
C  READS IN SETS OF WIND INPUT.  WIND SPEED AND WIND DIRECTION.
C  WIND SPEEDS ARE IN UNITS OF KNOTS AND ARE ASSUMED TO BE
C  REPRESENTATIVE OF A 10-METER HEIGHT IF INPLEV = 0 AND A
C  20-METER HEIGHT IF INPLEV = 1.
C
C    PARAMETER STATEMENT
C	 IDMN	= MAXIMUM NUMBER OF GRID COLUMNS
C	 JDMN	= MAXIMUM NUMBER OF GRID ROWS
C	 IF	= MAXIMUM NUMBER OF FREQUENCY BANDS
C	 IA	= MAXIMUM NUMBER OF ANGLE BANDS
C	 NOBPTS = MAXIMUM NUMBER OF INPUT BOUNDARY POINTS
C	 NBN	= MAXIMUM NUMBER OF INPUT BOUNDARY POINTS
C    PARAMETERS PASSED
C	 ID	= DATE TIME CODE
C	 IENDWN = ERROR CHECK PARAMETER FOR ID
C    DATA READ IN
C-----
C	 ID		  = DATE TIME CODE
C	 WSN(I,J),I=IB,IE = WIND SPEED IN KNOTS
C	 WDN(I,J),I=IB,IE = WIND DIRECTION WIS FORMAT
C
      INCLUDE 'PARAMS'
C      PARAMETER (IDMN=35,JDMN=51,IF=15,IA=16,NOBPTS=11,NBN=11)
C	   READ IN ARRAYS OF WIND SPEEDS (IWSN) AND WIND DIRECTIONS (IWDN).
C	   CONVERT WIND DIRECTIONS TO RADIANS.
      COMMON /WN/ IWSL(IDMN,JDMN),IWDL(IDMN,JDMN),IWS(IDMN,JDMN),
     2		  IWD(IDMN,JDMN),IWSN(IDMN,JDMN),IWDN(IDMN,JDMN)
      COMMON /A3/ G,PI,XN,DL,NANG,NFRQ,NX,NY,TWOPI,FAC,DT,RADDEG
      COMMON /OPT/ MSTA,NSTR,NORD,JPRC,NTMS,INDELF,NHR,IOUT(70),
     2		   JOUT(70),KTIMES,WL(725),WN(725),IG21,IR21,MXHR,
     3		   IDEEP,DLAT,XLAT0,DZRO,NSTA(70),ALT(70),ALON(70)
	COMMON /DEPN/ DEPNEW(IDMN,JDMN)
      COMMON /IOP/ ICURV,IDPRP,IUPDAT
      DIMENSION WSN(IDMN,JDMN),WDN(IDMN,JDMN)

C     PRINT *,'INSIDE SUB RDWIN'
      READ (21,210,END=98) ID
c      write (15,210) id
  210 FORMAT (I10)
C     READ WIND SPEEDS (INPUT IS M/SEC)
      DO 211 J=NY,1,-1
      READ (21,*) (WSN(I,J),I=1,NX)
  211 CONTINUE
C     READ WIND DIRECTIONS (INPUT IS DEGREES CLOCKWISE FROM AZIMUTH)
C                          (FROM WHICH THE WIND IS COMING)
      DO 212 J=NY,1,-1
      READ (21,*) (WDN(I,J),I=1,NX)
  212 CONTINUE


C  CODE IS EXPECTING KNOTS AND CARTESIAN DIRECTION CONVENTION
C  THAT IS TOWARD EAST =0,TOWARD NORTH =90, ETC.
C	   CONVERT WIND DIRECTIONS TO CODES (DIRECTION BANDS)
C	   ROUND OFF WIND SPEEDS
   97 CONTINUE
      DO 450 I=1,NX
      DO 450 J=1,NY
      IWDN(I,J)=INT(WDN(I,J)/5.+1.5)
	 IWSN(I,J)=WSN(I,J)+0.5
  450 CONTINUE

C
      RETURN
   98 IENDWN=1

      RETURN
      END
C
C
C
      SUBROUTINE NDATE (IDL,KTIME,IDN)
C     SUBROUTINE NDATE (IDL,KTIME,TINC,IDN)
      COMMON /A3/ G,PI,XN,DL,NANG,NFRQ,NX,NY,TWOPI,FAC,TINC,RADDEG
C
C	    TINC*MAXKTIME SHOULD NEVER BE MORE THAN 1 DAY
C
      INC=TINC
      IDN=IDL+KTIME*INC/3600
      IDCK=IDN/100
      IHRCK=IDN-IDCK*100
C
C	    CHECK FOR MORE THAN 24 HOURS
C
      IF (IHRCK.LT.24) GO TO 1
      IDN=(IDCK+1)*100+(IHRCK-24)
      IDCK1=IDN/10000
      IDCK2=IDN/100
      IDAYCK=IDCK2-IDCK1*100
C
C	  CHECK FOR END OF MONTH, IDAYCK CONTAINS THE DAY
C
      IF (IDAYCK.LT.29) GO TO 1
      IYR=IDN/1000000
      IYRMN=IDN/10000
C
C	 IMNCK CONTAINS THE MONTH, IYR CONTAINS THE YEAR
C
      IMNCK=IYRMN-IYR*100
      GO TO (10,20,30,40) 33-IDAYCK
   97 PRINT 98, IDAYCK,IMNCK
   98 FORMAT (10X,'IDAYCK =',I12,5X,'IMNCK =',I6)
      RETURN
C
C	GET HERE IF IDAYCK = 32
C
   10 IF (IMNCK.EQ.2.OR.IMNCK.EQ.4.OR.IMNCK.EQ.6) GO TO 97
      IF (IMNCK.EQ.9.OR.IMNCK.EQ.11) GO TO 97
      GO TO 2
C
C	  GET HERE IF IDAYCK = 31
C
   20 IF (IMNCK.EQ.2) GO TO 97
      IF (IMNCK.NE.4.AND.IMNCK.NE.6.AND.IMNCK.NE.9.AND.
     * IMNCK.NE.11) GO TO 1
      GO TO 2
C
C	   GET HERE IF IDAYCK = 30
C
   30 IF ( IMNCK.EQ.2) GO TO 97
      GO TO 1
C
C	  GET HERE IF IDAYCK = 29
C
   40 IF (IMNCK.NE.2) GO TO 1
C
C	  MAKE ID - - MONTH+1,DAY=1,HOUR-24
C
    2 IDN=IYR*1000000+(IMNCK+1)*10000+100+(IHRCK-24)
      IF ((IMNCK+1).NE.13) GO TO 1
C
C	   CHANGE YEAR
C
      IDN=(IYR+1)*1000000+10000+100+(IHRCK-24)
   1	  CONTINUE
C   1 PRINT 91, IDL, KTIME,TINC,IDN
   91 FORMAT (1X,'IDL =',I10,5X,'KTIME =',I3,5X,'TINC =',F5.0,
     *	5X,'IDN =',I10)
      RETURN
      END
C
C
C
      SUBROUTINE OUTP(IDN)
      INCLUDE 'PARAMS'
C     PARAMETER (INMN=32, JDMN=51, IF=15, IA=16, NOBPTS=11, NBN=11)
      COMMON /MAXOUT/ WAVMAX(IDMN,JDMN),TPPMAX(IDMN,JDMN),
     1		      WDDMAX(IDMN,JDMN)
      COMMON /GAMALP/ GAMOUT(IDMN,JDMN),ALPOUT(IDMN,JDMN)
      COMMON /CDR/COEF1,COEF2,SPDFAC,XLAMZ,ALPSTR,BETA1
      COMMON /SWL/ HSWL(IDMN,JDMN),TSWL(IDMN,JDMN),SWLDIR(IDMN,JDMN)
      COMMON/CEN/E(IDMN,JDMN,IF,IA),ENXT(NBN,IF,IA),ELST(NBN,IF,IA)
      COMMON /O1/ FREQ(IF),SINR(IA),COSR(IA),ALFA0(IA),ALFAR0(IA),
     2		  DELF(IF),AINC,AINCR
      COMMON /A3/ G,PI,XN,DL,NANG,NFRQ,NX,NY,TWOPI,FAC,DT,RADDEG
      COMMON /O5/ IBOUND(IDMN,JDMN)
      COMMON /OPT/ MSTA,NSTR,NORD,JPRC,NTMS,INDELF,NHR,
     2	      IOUT(70),JOUT(70),KTIMES,WL(725),WN(725),IG21,IR21,
     3	      MXHR,IDEEP,DLAT,XLAT0,DZRO,NSTA(70),ALT(70),ALON(70)
      COMMON /WN/ IWSL(IDMN,JDMN),IWDL(IDMN,JDMN),IWS(IDMN,JDMN),
     2		  IWD(IDMN,JDMN),IWSN(IDMN,JDMN),IWDN(IDMN,JDMN)
      COMMON /S4/ HSIG(IDMN,JDMN)
      COMMON /S2/ FMA(IDMN,JDMN),
     2	     SWANG(IDMN,JDMN),KFRQA(IDMN,JDMN)
      COMMON /S7/ AVANG(IDMN,JDMN),TSS(IDMN,JDMN)
	COMMON /DEPTH/ DEP(IDMN,JDMN)
	COMMON /DEPN/ DEPNEW(IDMN,JDMN)
	  COMMON /DEPOR/ DEPORG(IDMN,JDMN)
      COMMON /TM/ TIMTOT
	COMMON /HSC/ HSCALE(IDMN,JDMN)
C	DIMENSION IE2(IDMN,JDMN,IF,IA),ENERGY(IDMN,JDMN,IF,IA)
      DIMENSION ETMP1(IF),iwh(idmn)

c      file 10 contains 2-d spectra
c      file 13 contains 1-d spectra
c      file 11 contains 1-line summary

      write (13,1312) idn
 1312 format (i8)
      write (10,1312) idn

      do 111 ista=1,msta
      i=iout(ista)
      j=jout(ista)
      ttt=1./(fma(i,j)+0.001)
      tswlt=tswl(i,j)
      wsp=iws(i,j)
      wdr=(iwd(i,j)-1)*5
      avangt=avang(i,j)*raddeg
      swangt=swang(i,j)*raddeg
      swldt=swldir(i,j)*raddeg

      write (10,1314) i,j
 1314 format (2i5)

      EMAX=0.
      FFMAX=1.0
      DO 1359 IFRC=1,NFRQ
      ETMP1(IFRC)=0.
      DO 1533 IANG=1,NANG
      ETMP1(IFRC)=ETMP1(IFRC)+E(I,J,IFRC,IANG)*AINCR
 1533 CONTINUE

      IF (EMAX.LT.ETMP1(IFRC)) THEN
         EMAX=ETMP1(IFRC)
         FFMAX=FREQ(IFRC)
         ENDIF

      write (10,1313) (e(i,j,ifrc,iang),iang=1,nang)
 1313 format (16f8.3)

 1359 CONTINUE
      TTM=1.0/FFMAX

      write (11,1111)idn,i,j,hsig(i,j),ttm,ttt,tswlt,avangt,swangt,
     2    swldt,wsp,wdr
 1111 format (i8,2i3,8f6.1,2f5.0)

      write (13,1314) i,j
      write (13,1313) (etmp1(k),k=1,nfrq)

  111 continue

       write (15,1111) idn

       do 1112 j=ny,1,-1
       do 113 i=1,nx
       iwh(i)=hsig(i,j)
       if (ibound(i,j).gt.1) iwh(i)=-1
  113 continue
      write (15,114) (iwh(i),i=1,nx)
  114 format (60i2)
 1112 continue

      do 115 j=ny,1,-1
      do 116 i=1,nx
      iwh(i)=1./(fma(i,j)+0.001)
      if (ibound(i,j).gt.1) iwh(i)=-1
      if (iwh(i).gt.99) iwh(i)=99
  116 continue
      write (15,114) (iwh(i),i=1,nx)
  115 continue

      do 1117 j=ny,1,-1
      do 1118 i=1,nx
      iwh(i)=avang(i,j)*raddeg /5.
      if (iwh(i).lt.0) iwh(i)=iwh(i)+72
 1118 continue
      write (15,114) (iwh(i),i=1,nx)
 1117 continue


       RETURN
       END
C
C
C
C
C
      SUBROUTINE CKSTR(ICALL)
      INCLUDE 'PARAMS'
C     PARAMETER (IDMN=35,JDMN=51,IF=15,IA=16,NOBPTS=11,NBN=11)
      COMMON /A3/ G,PI,XN,DL,NANG,NFRQ,NX,NY,TWOPI,FAC,DT,RADDEG
      COMMON /S2/ FMA(IDMN,JDMN),
     2	     SWANG(IDMN,JDMN),KFRQA(IDMN,JDMN)
      COMMON /S5/ FMAN(IDMN,JDMN),E2(IF,IA),
     2	     EF(IF),EFN(IF),FKA(IF)
      COMMON /S6/ ES(IDMN,JDMN,IF),IWVDIR(IDMN,JDMN)
      COMMON/CEN/E(IDMN,JDMN,IF,IA),ENXT(NBN,IF,IA),ELST(NBN,IF,IA)
      COMMON /OPT/ MSTA,NSTR,NORD,JPRC,NTMS,INDELF,NHR,IOUT(70),
     2		   JOUT(70),KTIMES,WL(725),WN(725),IG21,IR21,MXHR,
     3		   IDEEP,DLAT,XLAT0,DZRO,NSTA(70),ALT(70),ALON(70)
C
      GO TO (100,200) ICALL
      PRINT *,'INCORRECT ICALL IN CKSTR'
      STOP
  100 IF (NSTR.LE.1) RETURN
      READ (23) E,FMA
      READ (23) SWANG
      READ (23) KFRQA
      RETURN
C
  200 IF (NSTR.NE.1.AND.NSTR.NE.3) RETURN
      WRITE (16) E,FMA
      WRITE (16) SWANG
      WRITE (16) KFRQA
      RETURN
      END
C
C
C
      SUBROUTINE WBOUND(IDN)
      INCLUDE 'PARAMS'
C     PARAMETER (IDMN=35,JDMN=51,IF=15,IA=16,NOBPTS=11,NBN=11)
      COMMON /A3/ G,PI,XN,DL,NANG,NFRQ,NX,NY,TWOPI,FAC,DT,RADDEG
      COMMON /PMOD/ XD(IF,JDMN,IA),XDU(IF,JDMN,IA),YD(IF,IA),
     1		   YDU(IF,IA)
      COMMON /S2/ FMA(IDMN,JDMN),
     2	     SWANG(IDMN,JDMN),KFRQA(IDMN,JDMN)
      COMMON /BD/ IBND,NBPO,IBPO(NOBPTS),JBPO(NOBPTS),NBPI,KNTBP,
     2	     IBPI(NBN),JBPI(NBN)
      COMMON/CEN/E(IDMN,JDMN,IF,IA),ENXT(NBN,IF,IA),ELST(NBN,IF,IA)
C
C  NEXT TWO COMMON LINES ADDED TO WRITE BOUN OUTPUT TO SHALWV
      COMMON /S6/ ES(IDMN,JDMN,IF),IWVDIR(IDMN,JDMN)
      COMMON /GAMALP/ GAMOUT(IDMN,JDMN),ALPOUT(IDMN,JDMN)
      COMMON /DEPTH/ DEP(IDMN,JDMN)
C NBPO IS # OF BOUNDARY POINTS WRITTEN AS OUTPUT
C NOBPTS IS MAXIMUM NUMBER OF THESE POINTS DEFINED IN
C    THE PARAMETER STATEMENT
C DATA WRITTEN INCLUDES DATE, POINT COUNTER NUMBER,
C    (I,J) LOCATION OF BOUNDARY POINT, PEAK FREQUENCY,
C    AND THE ENERGY MATRIX
      DO 100 K=1,NBPO
      I=IBPO(K)
      J=JBPO(K)
      WRITE (41,*) IDN,K,I,J,FMA(I,J)
      WRITE (41,7171) ((E(I,J,IFRC,IANG),IANG=1,NANG),IFRC=1,NFRQ)
 7171 FORMAT (16F8.3)
C5050 CONTINUE
  100 CONTINUE
      PRINT *,'BOUNDARY DATA WRITTEN ',IDN
      RETURN
      END
C
C
C
C
      SUBROUTINE SETSID(NX,NY,NF,NA,ISPEC)
      INCLUDE 'PARAMS'
C     PARAMETER (IDMN=35,JDMN=51,IF=15,IA=16,NOBPTS=11,NBN=11)
      COMMON /O5/ IBOUND(IDMN,JDMN)
      COMMON /CEN/ E(IDMN,JDMN,IF,IA),ENXT(NBN,IF,IA),ELST(NBN,IF,IA)
      COMMON /S2/ FMA(IDMN,JDMN),
     2	     SWANG(IDMN,JDMN),KFRQA(IDMN,JDMN)
      IF (ISPEC.EQ.0) RETURN
      IF (ISPEC.NE.1) GO TO 10
      DO 1 I=2,NX-1
      IF (IBOUND(I,J).NE.1) GO TO 1
      KFRQA(I,1)=KFRQA(I,2)
      KFRQA(I,NY)=KFRQA(I,NY-1)
      FMA(I,1)=FMA(I,2)
      FMA(I,NY)=FMA(I,NY-1)
      DO 2 IAA=1,NA
      DO 2 ITF=1,NF
      E(I,1,ITF,IAA)=E(I,2,ITF,IAA)
      E(I,NY,ITF,IAA)=E(I,NY-1,ITF,IAA)
    2 CONTINUE
    1 CONTINUE
      RETURN
   10 IF (ISPEC.NE.2) GO TO 20
      DO 3 J=2,NY-1
      IF (IBOUND(I,J).NE.1) GO TO 3
      KFRQA(1,J)=KFRQA(2,J)
      KFRQA(NX,J)=KFRQA(NX-1,J)
      FMA(1,J)=FMA(2,J)
      FMA(NX,J)=FMA(NX-1,J)
      DO 4 IAA=1,NA
      DO 4 ITF=1,NF
      E(1,J,ITF,IAA)=E(2,J,ITF,IAA)
      E(NX,J,ITF,IAA)=E(NX-1,J,ITF,IAA)
    4 CONTINUE
    3 CONTINUE
      RETURN
   20 DO 5 I=2,NX-1
      IF (IBOUND(I,J).NE.1) GO TO 5
      KFRQA(I,1)=KFRQA(I,2)
      KFRQA(I,NY)=KFRQA(I,NY-1)
      FMA(I,1)=FMA(I,2)
      FMA(I,NY)=FMA(I,NY-1)
      DO 6 IAA=1,NA
      DO 6 ITF=1,NF
      E(I,1,ITF,IAA)=E(I,2,ITF,IAA)
      E(I,NY,ITF,IAA)=E(I,NY-1,ITF,IAA)
    6 CONTINUE
    5 CONTINUE
      DO 7 J=2,NY-1
      IF (IBOUND(I,J).NE.1) GO TO 7
      KFRQA(1,J)=KFRQA(2,J)
      KFRQA(NX,J)=KFRQA(NX-1,J)
      FMA(1,J)=FMA(2,J)
      FMA(NX,J)=FMA(NX-1,J)
      DO 8 IAA=1,NA
      DO 8 ITF=1,NF
      E(1,J,ITF,IAA)=E(2,J,ITF,IAA)
      E(NX,J,ITF,IAA)=E(NX-1,J,ITF,IAA)
    8 CONTINUE
    7 CONTINUE
      RETURN
      END
C
C
C
C
      SUBROUTINE INITB
C  THIS SUBROUTINE CHECKS THE LAND-SEA GRID TO
C  DEFINE THE LOCATION OF INPUT BOUNDARY ENERGY
C  LOCATIONS.  IT LOCATES BOUNDARY POINTS INTERPOLATED
C  BETWEEN THESE INPUT LOCATIONS AND SETS UP A
C  WEIGHTING SCHEME FOR THE INTERPOLATION OF THE BOUNDARY
C  INPUT ENERGY. IT ALSO IDENTIFIES POINTS EXTRAPOLATED
C  FROM THESE INPUT LOCATIONS.
C
C   NOTE THAT PARAMETER NEEDS VALUES FOR NOBPTS AND NBN.
C   NBN IS THE TOTAL NUMBER OF BOUNDARY POINTS TO BE INPUT
C     INCLUDING INTERPOLATED AND EXTRAPOLATED POINTS.
C   NOBPTS ARE THE NUMBER OF OUTPUT BOUNDARY LOCATIONS
C     WHERE DATA WILL BE WRITTEN.
      INCLUDE 'PARAMS'
C     PARAMETER (IDMN=35,JDMN=51,IF=15,IA=16,NOBPTS=11,NBN=11)
      COMMON /HBND/  IBEGA(NBN),JBEGA(NBN),IENDA(NBN),JENDA(NBN),
     1 ITB(NBN,12),JTB(NBN,12),NLA(NBN),WTB(NBN,12),
     2	WTA(NBN,12),NSEG,NINT,ILS,JLS,ICBI(NBN),JCBI(NBN),
     3	IF1,JF1,IFN,JFN,NEXF1,NINPT,NINTR,NEXFN,NOF1,NOFN,
     4	IBOUN(IDMN,JDMN),ICLO
      COMMON /BD/ IBND,NBPO,IBPO(NOBPTS),JBPO(NOBPTS),NBPI,KNTBP,
     2	     IBPI(NBN),JBPI(NBN)
      COMMON /O5/ IBOUND(IDMN,JDMN)
      COMMON /A3/ G,PI,XN,DL,NANG,NFRC,NX,NY,TOPI,FAC,DT,RADDEG
      DIMENSION IAA(8),JAA(8)
C THIS DATA INCREMENTS I AND J TO SEARCH FOR
C  A NEARBY POINT
      DATA IAA /1,0,-1,0,1,-1,-1,1/
      DATA JAA /0,-1,0,1,-1,-1,1,1/
C CALCULATE NUMBER OF POINTS IN THE LARGEST INTERIOR
C RECTANGLE WITHIN THE LAND-SEA GRID--USED TO SEE IF
C A CLOSED LOOP BOUNDARY SITUATION EXISTS
      PRINT *, 'IN INITB'
      NIMAX=(2*(NX-2))+(2*(NY-4))
C	      CHECK INPUT STARTING LOCATION
      IBZO=IBOUND(ILS,JLS)
C BOUNDARY POINT MUST BE EITHER A REAL INPUT BOUNDARY
C POINT OR AN EXTRAPOLATED BOUNDARY POINT
      IF(IBZO.LE.1.OR.IBZO.EQ.4)THEN
C THIS ERROR INDICATES STARTING LOCATION IS AN INTERPOLATED
C POINT OR A REGULAR WATER POINT
      PRINT *,'ERROR IN INITB SPECIFIED STARTING LOCATION NOT VALID'
      PRINT *,'IBOUND(ILS,JLS) = ',IBOUND(ILS,JLS)
      STOP 'ABORT IN INITB'
      ENDIF
      PRINT *
      PRINT *,'BEGIN BOUNDARY INITIALIZATIONS'
      PRINT *
C		  SCAN THE IBOUND MATRIX TO CLASSIFY AND COUNT
C		  BOUNDARY POINTS & INITIALISE TEMPORARY ARRAY
C		  IBOUN
C NEXF1 ARE # OF LEADING EDGE PTS(EXTRAPOLATED) IBOUND=2
C  NOTE LEADING EDGE POINTS ARE EXTRAPOLATED FROM
C  FIRST REAL BOUNDARY INPUT POINT
C NINPT ARE # OF INPUT BOUNDARY PTS IBOUND=3
C NINTR ARE # OF INTERPOLATED BOUNDARY PTS IBOUND=4
C NEXFN ARE # OF TRAILING EDGE PTS (EXTRAPOLATED) IBOUND=5
C NOTE TRAILING EDGE POINTS ARE EXTRAPOLATED FROM NTH
C   REAL BOUNDARY INPUT POINT
      NEXF1=0
      NINPT=0
      NINTR=0
      NEXFN=0
C THIS LOOP COUNTS # OF POINTS IN ABOVE CATEGORIES
      DO 10 I=1,NX
      DO 10 J=1,NY
      IB=IBOUND(I,J)
      IF(IB.EQ.2)NEXF1=NEXF1+1
      IF(IB.EQ.3)NINPT=NINPT+1
      IF(IB.EQ.4)NINTR=NINTR+1
      IF(IB.EQ.5)NEXFN=NEXFN+1
      IBOUN(I,J)=IB
   10 CONTINUE
C NTOT IS TOTAL # OF POINTS IN A BOUNDARY CAPACITY
C NREAL IS TOTAL # OF NON-EXTRAPOLATED BOUNDARY POINTS
C NBPI ARE # OF INPUT BOUNDARY POINTS WHERE REAL ENERGY
C   IS COMING IN(NOT EXTRAPOLATED OR INTERPOLATED)
      NTOT=NEXF1+NINPT+NINTR+NEXFN
      NREAL=NINPT+NINTR
      NBPI=NINPT
      PRINT 11,NEXF1,NINPT,NINTR,NEXFN
   11 FORMAT(I6,' PTS WILL BE EXTRAPOLATED FROM 1ST REAL POINT',/,
     1	     I6,' PTS WILL BE READ FROM INPUT',/,
     2	     I6,' PTS WILL BE INTERPOLATED',/,
     3	     I6,' PTS WILL BE EXTRAPOLATED FROM NTH REAL POINT')
      DO 12 J=NY,1,-1
   12 PRINT 9,(IBOUN(I,J),I=1,NX)
    9 FORMAT(80I1)
      PRINT 13,NTOT
   13 FORMAT('TOTAL NUMBER OF BOUNDARY POINTS = ',I6)
C CHECK TO SEE IF PARAMETER STATEMENT IS EXCEEDED
      IF(NTOT.GT.NBN)THEN
	  PRINT *,'NUMBER OF BOUNDARY SEGMENTS EXCEEDS MAX. LIMIT'
	  STOP 'ABORT IN INITB'
      ENDIF
C	  CHECK FOR CLOSED BOUNDARY LOOP
C  ICLO=0 IS OPEN LOOP, ICLO=1 IS CLOSED LOOP
      ICLO=0
      IF(NTOT.GE.NIMAX)THEN
	  ICLO=1
C IF INTERPOLATED POINTS EXIST
	  IF(NINTR.NE.0)THEN
C ADD ONE MORE INPUT BOUNDARY POINT TO RECOUNT FIRST
C  POINT AS LAST POINT FOR CLOSURE
C RESET NEW TOTALS BASED ON ONE ADDITIONAL POINT
	     NINPT=NINPT+1
	     NTOT=NEXF1+NINPT+NINTR+NEXFN
	     NREAL=NINPT+NINTR
	  ENDIF
      ENDIF
      PRINT *,'MMX',NINPT,NTOT,NIMAX,ICLO
      IFR=0
      K=0
      NINT=1
      NQI=0
C	STEP AROUND THE PERIMETER OF THE IBOUND MATRIX
C	TO DETERMINE I AND J COORDINATES OF "REAL" POINTS
C	("REAL" POINTS ARE EITHER READ OR INTERPOLATED).
C	POINTS REQUIRING EXTRAPOLATION WILL BE PROCESSED
C	IN SUBROUTINE RDBND
      IPKN=0
   14 CONTINUE
C LOOP BY DIRECTION IN 45 DEG. INCREMENTS
      DO 15 IK=1,8
C	CHECK FOR "REAL" STARTING POINT
C IF THIS IS FIRST TIME THROUGH
      IF(IFR.EQ.0)THEN
	IF(IBOUN(ILS,JLS).EQ.3)THEN
C GO THROUGH FOLLOWING SECTION IF POINT IS A REAL BOUNDARY
C  INPUT LOCATION(NOT INTERPOLATED)
C IPKN IS COUNTER NUMBER OF THE BOUNDARY POINT(REAL AND
C    INTERPOLATED)
C NQI IS COUNTING ONLY REAL INPUT BOUNDARY POINTS
C I AND J LOCATIONS ARE INDEXED BY COUNTER NUMBER OF POINT
C   _BPI ARE INPUT BOUNDARY LOCATIONS
C   _CBI ARE LOCATIONS USING TOTAL NUMBER OF POINTS
	  IPKN=IPKN+1
	  NQI=NQI+1
	  IBPI(NQI)=ILS
	  JBPI(NQI)=JLS
	  ICBI(IPKN)=ILS
	  JCBI(IPKN)=JLS
C IFR OVER ZERO INDICATES FIRST POINT HAS BEEN RECORDED
	  IFR=1
C RECORD BEGINNING VALUES OF REAL BOUNDARY POINT
	  IF1=ILS
	  JF1=JLS
	  NOF1=NQI
	  PRINT *,'INPUT PT # ',NQI,' I = ',IBPI(NQI),'  J = ',JBPI
     +	  (NQI),NINPT
	  PRINT *
	  PRINT *,'PT # ',IPKN,' I = ',ICBI(IPKN),'  J = ',JCBI(IPKN)
C RESET IBOUN
	  IF(ICLO.EQ.0)IBOUN(ILS,JLS)=1
	  GO TO 14
C GO THROUGH THIS SECTION IF POINT IS AN EXTRAPOLATED
C  BOUNDARY POINT RATHER THAN A REAL BOUNDARY INPUT POINT
	ELSEIF(IBOUN(ILS,JLS).EQ.2)THEN
	  IPKN=IPKN+1
	  ICBI(IPKN)=ILS
	  JCBI(IPKN)=JLS
	  PRINT *,'PT # ',IPKN,' I = ',ICBI(IPKN),'  J = ',JCBI(IPKN)
	  IBOUN(ILS,JLS)=1
	  IFR=1
	  GO TO 14
	ENDIF
      ENDIF
C	FIND NEXT POINT
C USE(I,J) INCREMENTATION DATA TO LOCATE THE NEXT BOUNDARY POINT
      IPT=ILS+IAA(IK)
      JPT=JLS+JAA(IK)
      IVAL1=IBOUN(IPT,JPT)
      IF(IVAL1.EQ.3)THEN
C THIS IS AN INPUT BOUNDARY POINT (IVAL1=3)
	  IPKN=IPKN+1
	  ICBI(IPKN)=IPT
	  JCBI(IPKN)=JPT
C	NEXT POINT IS AN INPUT POINT
C	ASSIGN IBPI AND JBPI
	  NQI=NQI+1
	  IBPI(NQI)=IPT
	  JBPI(NQI)=JPT
	  PRINT *
	  PRINT *,'INPUT PT # ',NQI,' I = ',IBPI(NQI),'  J = ',JBPI
     +	   (NQI),NINPT
	  PRINT *
	  PRINT *,'PT # ',IPKN,' I = ',ICBI(IPKN),'  J = ',JCBI(IPKN)

C			RESET IBOUN
	  IBOUN(IPT,JPT)=1
	  ILS=IPT
	  JLS=JPT
C		  DEFINE INTERPOLATION START & END
	  IF(NQI.EQ.1)THEN
C RECORD BOUNDARY POINT NUMBER AND (I,J) OF FIRST REAL BOUNDARY
C   POINT
	      IF1=IPT
	      JF1=JPT
	      NOF1=NQI
	  ENDIF
	  IF(NQI.EQ.NINPT)THEN
C RECORD BOUNDARY POINT NUMBER AND (I,J) OF LAST
C   REAL BOUNDARY POINT
	      IFN=IPT
	      JFN=JPT
	      NOFN=NQI
	  ENDIF
C	CHECK LOOP EXIT CONDITION
	  IF(IPKN.EQ.NTOT)GO TO 16
	  GO TO 14
      ELSEIF(IVAL1.GE.2)THEN
C		NEXT POINT MUST BE EXTRAPOLATED OR INTERPOLATED
	  IPKN=IPKN+1
	  ICBI(IPKN)=IPT
	  JCBI(IPKN)=JPT
	  PRINT *,'PT # ',IPKN,' I = ',ICBI(IPKN),'  J = ',JCBI(IPKN)
	  IBOUN(IPT,JPT)=1
	  ILS=IPT
	  JLS=JPT
	  IF(IPKN.EQ.NTOT)GO TO 16
	  GO TO 14
      ENDIF
   15 CONTINUE
      GO TO 14
   16 CONTINUE
C			END LOOP
C BOUNDARY POINTS HAVE BEEN NUMBERED AND LOCATIONS ARE
C  DEFINED BY NUMBER OF POINT
      PRINT *
      PRINT *,'FIRST REAL POINT IS ',IF1,JF1,NOF1
      PRINT *,'LAST REAL POINT IS ',IFN,JFN,NOFN
      IF(ICLO.EQ.1)PRINT *,'LOOP IS CLOSED'
      PRINT *
      IFR=0
      ILS=IF1
      JLS=JF1
C			RESET IBOUN
      DO 17 I=1,NX
      DO 17 J=1,NY
      IBOUN(I,J)=IBOUND(I,J)
	IF(IBOUN(I,J).EQ.2) IBOUN(I,J)=1
	IF(IBOUN(I,J).EQ.5) IBOUN(I,J)=1
 17	CONTINUE
C		BEGIN LOOP ON NUMBER OF LINE SEGMENTS WHICH CONTAIN
C		POINTS REQUIRING INTERPOLATION FOR CALCULATION OF
C		WEIGHTS.
C BOUNDARY SEGMENT CONSISTS OF A SET OF INTERPOLATED
C POINTS BETWEEN TWO REAL BOUNDARY POINTS.  K IS THE
C COUNTER FOR THE SEGMENT
    1 K=K+1
      KFR=0
C		SET COORDINATES FOR LINE SEGMENT BEGINNING
C DEFINE BEGINNING LOCATION FOR SEGMENT
      IBEGA(K)=ILS
      JBEGA(K)=JLS
      IBOUN(ILS,JLS)=1
C DO NOT RESET IBOUN IF FIRST POINT HAS TO BE
C  USED AS LAST POINT
      IF(ICLO.EQ.1.AND.IFR.NE.1)IBOUN(ILS,JLS)=3
      IFR=1
    2 CONTINUE
      IF(JLS.GE.NY.OR.ILS.GE.NX)THEN
	  PRINT *,'ERROR IN INITB...  ...POINT OUT OF BOUNDS'
	  STOP 'ABORT IN INITB'
      ENDIF
C		FIND NEXT POINT
      DO 100 IK=1,8
      IPT=ILS+IAA(IK)
      JPT=JLS+JAA(IK)
      IVAL1=IBOUN(IPT,JPT)
      IF(IVAL1.LE.2.OR.IVAL1.GE.5)GO TO 100
      IF(IVAL1.GE.3.OR.IVAL1.LE.4)GO TO 101
  100 CONTINUE
      PRINT *,'ERROR:  MALFORMED LAND SEA GRID'
      STOP 'ABORT IN INITB'
  101 CONTINUE
C		"REAL" POINT FOUND
      NINT=NINT+1
      IF(IVAL1.EQ.3)GO TO 102
C		INTERPOLATED POINT FOUND
C KFR IS COUNTER FOR INTERPOLATED POINT WITHIN SEGMENT
      KFR=KFR+1
C CHECK TO SEE IF NUMBER OF INTERPOLATED POINTS EXCEEDS
C  BOUNDARY CONDITIONS (CODE ADDED 12/89)
	IF(KFR.LE.12) GO TO 104
	PRINT *,'ERROR:INTERP PTS.GT.12'
	STOP 'ABORT IN INITB'
 104	CONTINUE
C (I,J) LOCATIONS DESCRIBED BY (SEG.#, PT. # WITHIN SEG.)
      ITB(K,KFR)=IPT
      JTB(K,KFR)=JPT
      JLS=JPT
      ILS=IPT
      IBOUN(ILS,JLS)=1
      GO TO 2
  102 CONTINUE
C		INPUT POINT FOUND:
C			ASSIGN COORDINATES OF SEGMENT ENDPOINT
      IENDA(K)=IPT
      JENDA(K)=JPT
C NUMBER OF LAST INTERPOLATED POINT WITHIN SEGMENT
      NLA(K)=KFR
C		COMPUTE WEIGHTS TO BE USED FOR INTERPOLATION
C		OF POINTS WITHIN THIS LINE SEGMENT
      JDIF=JENDA(K)-JBEGA(K)
      IDIF=IENDA(K)-IBEGA(K)
C DIST IS TOTAL DISTANCE BETWEEN TWO BOUNDARY PTS
      DIST=SQRT(FLOAT(JDIF*JDIF+IDIF*IDIF))
C WTB IS DISTANCE FROM BEGINNING OF SEGMENT TO
C  INTERPOLATED POINT DIVIDED BY TOTAL LENGTH
C  OF SEGMENT
C WTA IS DISTANCE FROM END OF SEGMENT TO INTERPOLATED
C  POINT DIVIDED BY TOTAL LENGTH OF SEGMENT
      DO 103 I=1,KFR
      JD=JTB(K,I)-JBEGA(K)
      ID=ITB(K,I)-IBEGA(K)
      DIS=SQRT(FLOAT(JD*JD+ID*ID))
      WTB(K,I)=DIS/DIST
      WTA(K,I)=1.0-WTB(K,I)
  103 CONTINUE
      ILS=IPT
      JLS=JPT
      IBOUN(ILS,JLS)=1
C		CHECK FOR LAST POINT (EXIT CONDITION)
      IF(ILS.EQ.IFN.AND.JLS.EQ.JFN)GO TO 500
      GO TO 1
  500 CONTINUE
      KNTBP=NTOT
      NSEG=K
      PRINT *
      PRINT *,'INITIALIZE BOUNDARY INTERPOLATION PARAMETERS'
      PRINT *
      PRINT *,'NREAL = ',NREAL
      PRINT *
      PRINT *,'NSEG = ',NSEG
      DO 650 K=1,NSEG
      IF(NLA(K).EQ.0)GO TO 650
      PRINT *
      PRINT *,IBEGA(K),JBEGA(K),IENDA(K),JENDA(K),NLA(K)
      PRINT *
      PRINT *,' NPSEG = ',NLA(K)
      DO 600 L=1,NLA(K)
      PRINT *,ITB(K,L),JTB(K,L),WTA(K,L),WTB(K,L)
  600 CONTINUE
  650 CONTINUE
      PRINT *
      DO 700 J=NY,1,-1
      PRINT 9,(IBOUN(I,J),I=1,NX)
  700 CONTINUE
      PRINT *
      PRINT *,'NBPI = ',NBPI
      PRINT *,'IBPI = ',(IBPI(I),I=1,NBPI)
      PRINT *,'JBPI = ',(JBPI(I),I=1,NBPI)
      PRINT *
      PRINT *,'NCBI = ',KNTBP
      PRINT *,'ICBI = ',(ICBI(I),I=1,KNTBP)
      PRINT *,'JCBI = ',(JCBI(I),I=1,KNTBP)
      PRINT *
      PRINT *,' END OF SUBROUTINE'
      RETURN
      END
C
C
C
C
      SUBROUTINE RDBND(ID)
      INCLUDE 'PARAMS'
C     PARAMETER (IDMN=35,JDMN=51,IF=15,IA=16,NOBPTS=11,NBN=11)
      COMMON /HBND/  IBEGA(NBN),JBEGA(NBN),IENDA(NBN),JENDA(NBN),
     1 ITB(NBN,12),JTB(NBN,12),NLA(NBN),WTB(NBN,12),
     2	WTA(NBN,12),NSEG,NINT,ILS,JLS,ICBI(NBN),JCBI(NBN),
     3	IF1,JF1,IFN,JFN,NEXF1,NINPT,NINTR,NEXFN,NOF1,NOFN,
     4	IBOUN(IDMN,JDMN),ICLO
      COMMON /BD/ IBND,NBPO,IBPO(NOBPTS),JBPO(NOBPTS),NBPI,KNTBP,
     2	     IBPI(NBN),JBPI(NBN)
      COMMON /O5/ IBOUND(IDMN,JDMN)
      COMMON /A3/ G,PI,XN,DL,NANG,NFRC,NX,NY,TOPI,FAC,DT,RADDEG
      COMMON /BDN/ FMANXT(NBN),ALPNXT(NBN),SWANXT(NBN),
     2	     GAMNXT(NBN),KFRNXT(NBN),WKMNXT(NBN),IWVNXT(NBN)
      COMMON /CEN/ E(IDMN,JDMN,IF,IA),ENXT(NBN,IF,IA),ELST(NBN,IF,IA)
      COMMON /O4/ EN(IDMN,JDMN,IA)
      DIMENSION IDUM(4),FMAT(IDMN,JDMN)
C
C	READ BOUNDARYS AND INTERPOLATE VALUES FOR ALL PERIMETER POINTS
C
C    FIRST SECTION OF CODE DEALS WITH READ(BINARY) OF THE
C      BOUNDARY INPUT INFORMATION FROM A PREVIOUS RUN
C    ENERGY IS PLACED IN PROPER BOUNDARY LOCATION (I,J) OF WAVE GRID
C    INTERPOLATED ENERGY IS PUT IN ITS PROPER LOCATION
C    LAST SECTION OF CODE DEALS WITH FILLING THE EXTRAPOLATED
C      BOUNDARY POINTS (LEADING AND TRAILING)
      PRINT *,'READING BOUNDARY INPUT'
      IF(NINT.NE.0)PRINT *,'I WILL PRINT DIAGNOSTICS ON FIRST CALL'
C NBPI IS # OF BOUNDARY POINTS WHERE REAL (NOT INTERPOLATED
C  OR EXTRAPOLATED) ENERGY IS COMING IN
C SEE SUB WBOUND FOR DEFINITION OF VARIABLES WRITTEN

      print *,'nfrc, nang ',nfrc,nang

      DO 30 K=1,NBPI
      IPT=IBPI(K)
      JPT=JBPI(K)
      READ (25,*,END=99999) (IDUM(L),L=1,4),FMAT(IPT,JPT)
      print *,(idum(l),l=1,4),fmat(ipt,jpt)
      READ (25,4141)((E(IPT,JPT,IFRC,IANG),IANG=1,NANG),IFRC=1,NFRC)
      PRINT *,'K = ',K,'  IPT = ',IPT,'  JPT = ',JPT,IDUM,FMAT(IPT,JPT)

c      write (41,4141) ((e(ipt,jpt,ifrc,iang),iang=1,16),ifrc=1,nfrc)
 4141 format (16f8.3)


   30 CONTINUE
      PRINT *,NBPI,' BOUNDARY POINTS READ'
      ID=IDUM(1)
C KNT IS THE NUMBER OF THE FIRST BOUNDARY POINT
C (ONE BEYOND THE LAST LEADING EDGE EXTRAPOLATED POINT)
C KNT WILL BE INCREMENTED IN FOLLOWING LOOPS AND WILL
C  FUNCTION AS THE BOUNDARY POINT COUNTER NUMBER
      KNT=NEXF1+1
C K IN FOLLOWING LOOP IS # OF SEGMENTS BETWEEN
C   INPUT LOCATIONS, NBPI
C NSEG IS MAX NUMBER OF SEGMENTS
C (I,J) LOCATIONS HAVE BEEN DEFINED IN INITB
C   FOR THE BEGINNING AND END OF EACH SEGMENT
C NPS IS # OF INTERPOLATED POINTS IN SEGMENT
      DO 50 K=1,NSEG
      IP1=IBEGA(K)
      JP1=JBEGA(K)
      IP2=IENDA(K)
      JP2=JENDA(K)
      NPS=NLA(K)
C FOLLOWING IF FOR NO INTERPOLATED POINTS
      IF(NPS.EQ.0)GO TO 54
      IF(NINT.NE.0)THEN
	  PRINT *,'BEGIN INTERPOLATION IN SEGMENT ',K
	  PRINT *,'THERE ARE ',NPS,' POINTS'
	  PRINT *,'COORDINATES OF SEGMENT END POINTS:'
	  PRINT *,'IP1 = ',IP1,'   JP1 = ',JP1
	  PRINT *,'IP2 = ',IP2,'   JP2 = ',JP2
      ENDIF
   54 CONTINUE
C THIS LOOP PUTS THE ENERGY IN THE MATRIX AT THE ACTUAL
C  BOUNDARY ENERGY LOCATION AT THE BEGINNING OF THE SEGMENT
C PEAK FREQUENCY MATRIX IS ALSO REFERENCED BY THIS
C   BOUNDARY POINT NUMBER
      DO 55 I=1,IF
      DO 55 J=1,IA
   55 ENXT(KNT,I,J)=E(IP1,JP1,I,J)
      FMANXT(KNT)=FMAT(IP1,JP1)
      IF(NINT.NE.0)PRINT *,'FMANXT(',KNT,') = ',FMANXT(KNT)
C INCREMENT BOUNDARY POINT COUNTER
      KNT=KNT+1
C INTERPOLATION LOOP
C INTERPOLATE ENERGY BETWEEN INPUT BOUNDARY LOCATIONS
C  USING WEIGHTING FACTORS FROM INITB
C PUT ENERGY AND PEAK FREQUENCY INTO MATRICES USING
C   BOUNDARY POINT COUNTER
      DO 60 L=1,NPS
      DO 70 I=1,IF
      DO 70 J=1,IA
      ENXT(KNT,I,J)=E(IP1,JP1,I,J)*WTA(K,L) + E(IP2,JP2,I,J)*WTB(K,L)
   70 CONTINUE
      FMANXT(KNT)=FMAT(IP1,JP1)*WTA(K,L) + FMAT(IP2,JP2)*WTB(K,L)
      IF(NINT.NE.0)THEN
	  PRINT *,'POINT # ',L,' IN SEGMENT ',K
	  PRINT *,'WEIGHT(1) = ',WTA(K,L),'  WEIGHT(2)= ',WTB(K,L)
	  PRINT *,'FMANXT(',KNT,') = ',FMANXT(KNT)
      ENDIF
C INCREMENT # OF BOUNDARY COUNTER
   60 KNT=KNT+1
   50 CONTINUE
C AT THE END OF THE 50 LOOP ALL BOUNDARY LOCATIONS
C  UP TO THE END POINT OF THE LAST SEGMENT AND
C  THE TRAILING EDGE POINTS HAVE HAD THEIR ENERGY
C  MATRIX AND PEAK FREQUENCY MATRIX FILLED
      IPT=IENDA(NSEG)
      JPT=JENDA(NSEG)
C THIS LOOP FILLS THE ENDPOINT OF THE LAST SEGMENT
      DO 80 I=1,IF
      DO 80 J=1,IA
   80 ENXT(KNT,I,J)=E(IPT,JPT,I,J)
      FMANXT(KNT)=FMAT(IPT,JPT)
      IF(NINT.NE.0)PRINT *,'FMANXT(',KNT,') = ',FMANXT(KNT)
C RETURN IF CLOSURE EXISTS (BOUNDARY COMPLETE)
      IF(ICLO.EQ.1)THEN
	  NINT=0
	  RETURN
      ENDIF
C		SET EXTRAPOLATED POINTS
C BEGIN WITH LEADING EDGE POINTS
C NEXF1 IS # OF LEADING EDGE BOUNDARY POINTS
C NXF IS # OF FIRST REAL BOUNDARY INPUT LOCATION
C  (NXF IS ONE BEYOND BEGINNING EXTRAPOLATED POINTS
C    AND CORRESPONDS TO COUNTER FOR FIRST REAL
C    BOUNDARY INPUT POINT)
C LOOP 100 FILLS LEADING EDGE BOUNDARY ENERGY
C   LOCATIONS WITH ENERGY AT FIRST REAL BOUNDARY
C   INPUT LOCATION
      NXF=NEXF1+1
      DO 100 II=1,NEXF1
      FMANXT(II)=FMANXT(NXF)
      IF(NINT.NE.0)PRINT *,'FMANXT(',II,') = ',FMANXT(II)
      DO 90 I=1,IF
      DO 90 J=1,IA
   90 ENXT(II,I,J)=ENXT(NXF,I,J)
  100 CONTINUE
C NXN IS # OF FIRST TRAILING EDGE EXTRAPOLATED
C  POINT.  KNT REFERENCES LAST REAL BOUNDARY
C  INPUT LOCATION.  ENERGY MATRICES ARE FILLED
C  WITH VALUES FROM THE LAST REAL BOUNDARY INPUT
C  LOCATION.
      NXN=NEXF1+NINPT+NINTR+1
      DO 120 II=NXN,KNTBP
      FMANXT(II)=FMANXT(KNT)
      IF(NINT.NE.0)PRINT *,'FMANXT(',II,') = ',FMANXT(II)
      DO 110 I=1,IF
      DO 110 J=1,IA
  110 ENXT(II,I,J)=ENXT(KNT,I,J)
  120 CONTINUE
      NINT=0
      RETURN
99999 PRINT *,'ERROR:  END OF FILE IN RDBND'
      STOP 'ABORT IN RDBND'
      END
C
C
C
C
      SUBROUTINE DPROP(K)
      INCLUDE 'PARAMS'
C     PARAMETER (IDMN=35,JDMN=51,IF=15,IA=16,NOBPTS=11,NBN=11)

C     SUBROUTINE FOR PURELY DEEP-WATER PROPAGATION
      COMMON /PMOD/XD(IF,JDMN,IA),XDU(IF,JDMN,IA),YD(IF,IA),YDU(IF,IA)
      COMMON /O5/ IBOUND(IDMN,JDMN)
      COMMON /A3/ G,PI,XN,DL,NANG,NFRC,NX,NY,TWOPI,FAC,DT,RADDEG
      COMMON /DPR/ IUB(IA),JUA(IA),IAUA(IA,JDMN),AD(IF,JDMN,IA),
     2		  ADU(IF,JDMN,IA)
      COMMON /O4/ EN(IDMN,JDMN,IA)
      COMMON /CEN/ E(IDMN,JDMN,IF,IA),ENXT(NBN,IF,IA),ELST(NBN,IF,IA)
      COMMON /IOP/ ICURV,IDPRP,IUPDAT
      COMMON /S6/ EONED(IDMN,JDMN,IF),IWVDIR(IDMN,JDMN)

C***  MOD 08/10/88
      COMMON /SSPT/ TRANSX(IDMN,JDMN),TRANSY(IDMN,JDMN)
C***

      NXM1=NX-1
      NYM1=NY-1
C
C     BEGIN PROPAGATE
C
C	   LOOPS FOR THE FIRST ORDER PROPAGATION
C	   THE GRID SPACING CHANGES ALONG J BUT
C	   REMAINS CONSTANT ALONG I
C	   THE 'U' IN XDU,YDU,IU & JU IS FOR UPSTREAM POINT
C	   TRANSX & TRANSY = 1; THIS ONLY CHANGES FOR A SUBGRID
      DO 100 I=2,NXM1
      DO 100 J=2,NYM1
C	   THIS ONLY PERFORMED ON WATER POINTS
      IF (IBOUND(I,J).NE.1) GO TO 100
      DO 101 IG=1,NANG
      IU=I-IUB(IG)
      EN(I,J,IG)=XD(K,J,IG)*E(I,J,K,IG)
     2		+XDU(K,J,IG)*E(IU,J,K,IG)*TRANSX(IU,J)
  101 CONTINUE
  100 CONTINUE

C     SET INFORMATION FROM BOUNDARIES INTO EN MATRIX

      DO 150 I=2,NXM1
      DO 150 J=2,NYM1
      IF (IBOUND(I,J).EQ.1) GO TO 150
      DO 151 IG=1,NANG
      EN(I,J,IG)=E(I,J,K,IG)
  151 CONTINUE
  150 CONTINUE

      DO 200 I=2,NXM1
      DO 200 J=2,NYM1
C	   THIS IS ONLY PERFORMED ON WATER POINTS
      IF (IBOUND(I,J).NE.1) GO TO 200
      DO 201 IG=1,NANG
      JU=J-JUA(IG)
      E(I,J,K,IG)=YD(K,IG)*EN(I,J,IG)
     2		 +YDU(K,IG)*EN(I,JU,IG)*TRANSY(I,JU)
  201 CONTINUE
  200 CONTINUE
      IF (ICURV.NE.1) RETURN
C
C     BEGIN CURVATURE ESTIMATION IN PROPAGATION
C
      DO 300 I=2,NXM1
      DO 300 J=2,NYM1
C      IF (IBOUND(I,J).EQ.0) GO TO 300
      DO 301 IG=1,NANG
      EN(I,J,IG)=E(I,J,K,IG)
  301 CONTINUE
  300 CONTINUE
      DO 400 I=2,NXM1
      DO 400 J=2,NYM1
      IF (IBOUND(I,J).NE.1) GO TO 400
      DO 401 IG=1,NANG
      IAU=IAUA(IG,J)
      E(I,J,K,IG)=AD(K,J,IG)*EN(I,J,IG)
     2		 +ADU(K,J,IG)*EN(I,J,IAU)
  401 CONTINUE
  400 CONTINUE
      RETURN
      END
C
C
C
C
      SUBROUTINE FLXSTD
C  THIS SUBROUTINE CALCULATES THE TWO DIMENSIONAL WAVE SPECTRUM
C  AT EACH GRID POINT FOR EACH TIME STEP. DO LOOPS IN THE FREQUENCY
C  DOMAIN ARE USED TO MODIFY ENERGY LEVELS BASED ON WIND INPUT AND
C  THE EFFECTS OF THE NON-LINEAR SOURCE TERM. REFER TO RESIO AND PERRIE
C  J. PHYSICAL OCN. VOL. 19, FEB.,1989, AND DOCUMENTATION REPORT FOR
C  TECHNICAL DETAILS.
C  THE FOLLOWING COMMENT LINES REFER TO THE EXECUTABLE LINES BELOW THE COMMENT.
      INCLUDE 'PARAMS'
C     PARAMETER (IDMN=35,JDMN=51,IF=15,IA=16,NOBPTS=11,NBN=11)
C  FIELDS OF GAMMA AND ALPHA
      COMMON /GAMALP/ GAMOUT(IDMN,JDMN),ALPOUT(IDMN,JDMN)
      COMMON /CDR/ COEF1,COEF2,SPDFAC,XLAM,ALPSTR,BETA1
C  DEPTHS OVER GRID
	COMMON /DEPN/ DEPNEW(IDMN,JDMN)
      COMMON /DEPTH/DEP(IDMN,JDMN)
      COMMON /OPT/ MSTA,NSTR,NORD,JPRC,NTMS,INDELF,NHR,
     2		  IOUT(70),JOUT(70),KTIMES,WL(725),WN(725),IG21,IR21
     3		  ,MXHR,IDEEP,DLAT,XLAT0,DZRO,NSTA(70),ALT(70),ALON(70)

C  ZJAC IS THE RATIO OF WAVE NUMBERS (DEEP/DEPTH DEPENDENT) TO THE 5/2 POWER.
C  IT REPRESENTS AN EQULIBRIUM OR LIMITATION ON ENERGY DUE TO LOCAL
C  STEEPNESS EXCEEDING A CRITICAL VALUE AND BREAKING OCCURRING.
C  SEE RESIO J.OF WATERWAY, PORT, COASTAL AND OCEAN ENGINEERING VOL.113, NO.3
C  MAY 1987, AND KITAIGORODSKII, JPO 13,1983.
      COMMON /JAC/ ZJAC(IDMN,JDMN,IF)
C  HEIGHT, PERIOD AND DIRECTION OF SWELL. THE SWELL PART OF THE SPECTRUM
C  IS FROM FREQUENCY BAND 1 TO 0.7 FM WHICH IS THE WIND SEA PEAK FREQUENCY.
      COMMON /SWL/ HSWL(IDMN,JDMN),TSWL(IDMN,JDMN),SWLDIR(IDMN,JDMN)
C  "WIND SEA" FREQUENCY SPECTRUM (M**2)
      DIMENSION FSEA(IF)
C  THE LAND - WATER MATRIX
      COMMON /O5/ IBOUND(IDMN,JDMN)
C  FREQUENCY ARRAY, SINES AND COSINES OF ANGLES, ANGLES IN DEG AND RADIANS
C  FREQUENCY INCREMENTS, ANGLE INCREMENTS IN DEG AND RADIANS.
      COMMON /O1/ FREQ(IF),SINR(IA),COSR(IA),ALFA0(IA),ALFAR0(IA),
     2	     DELF(IF),AINC,AINCR
      COMMON /A3/ G,PI,XN,XGRID,NANG,NFRQ,NI,NJ,TWOPI,XFAC,TINC,RADDEG
C LAST, PRESENT, AND NEXT WIND SPEED AND DIRECTION IN TIME.
      COMMON /WN/ IWSL(IDMN,JDMN),IWDL(IDMN,JDMN),IWND(IDMN,JDMN),
     2	     IWNDIR(IDMN,JDMN),IWSNXT(IDMN,JDMN),IWDNXT(IDMN,JDMN)
C  ENERGY AS A FUNCTION OF ANGLE
      COMMON /O4/ EN(IDMN,JDMN,IA)

C PRESENT, NEXT, AND LAST 2-D ENERGY SPECTRUM.
      COMMON/CEN/EL2(IDMN,JDMN,IF,IA),ENXT(NBN,IF,IA),ELST(NBN,IF,IA)
C  MAKE THE ABOVE ARRAYS INTEGERS
C  ARRAY OF COSINE ** N OF ANGLES EVERY 5 DEG, FOR N=2,4,6,8,10
C  ANGLES ARE EVERY 5 DEG ON EITHER SIDE OUT TO (+/- 40 DEG.) OF EACH ANGLE BAND
      COMMON /S1/ PANG(IA,72,5)
C  PEAK FREQUENCY, SEA WAVE ANGLE, AND BAND INDEX OF PEAK FREQ.
      COMMON /S2/ FMA(IDMN,JDMN),
     2	     SWANG(IDMN,JDMN),KFRQA(IDMN,JDMN)
      COMMON /S3/ NF1,RADC,PI72,PI8
C  ARRAY OF SIGNIFICANT WAVE HEIGHTS
      COMMON /S4/ HSIG(IDMN,JDMN)
C  PEAK FREQ., 2-D SPECTRUM, OLD + NEW 1-D SPECTRUM, G/(2*PI)**3*F**-4
      COMMON /S5/ FMAN(IDMN,JDMN),E2(IF,IA),
     2	     EF(IF),EFN(IF),FKA(IF)
C  1-D SPECTRUM, INDEX OF WAVE DIRECTION
      COMMON /S6/ E(IDMN,JDMN,IF),IWVDIR(IDMN,JDMN)
C  AVERAGE WAVE ANGLE, MEAN PERIOD
      COMMON /S7/ AVANG(IDMN,JDMN),TSS(IDMN,JDMN)
      COMMON /MM/ INPLEV,IWRAP
       COMMON/HSC/HSCALE(IDMN,JDMN)
C  GRID POINT LOCATIONS SURROUNDING POINT (I,J), SEE DATA STATEMENT BELOW.
      DIMENSION IPOS(8),JPOS(8)
C  PEAK FREQ. FROM PREVIOUS TIME STEP.
      DIMENSION FMA0(IDMN,JDMN)
C I AND J GRID POINT LOCATIONS TO DETERMINE UPSTREAM GRID POINT
      DATA IPOS /-1,-1,0,1,1,1,0,-1/
      DATA JPOS /0,-1,-1,-1,0,1,1,1/
C A SMALL VALUE OF ENERGY TO AVOID DIVIDING BY ZERO
      DATA EPSS /0.0001/

      PRINT *,'ENTERING FLXSTD'


C FREQUENCY USED IN DETERMINING PARAMETRIC ENERGY
      FMH=FREQ(NFRQ)+0.5*DELF(NFRQ)
C SET PEAK FREQUENCY TO VALUE AT PREVIOUS TIME STEP
      DO 4000 I=1,NI
      DO 4000 J=1,NJ
      FMA0(I,J)=FMA(I,J)
 4000 CONTINUE
C AN EMPIRICAL FACTOR USED IN FOLLOWING EQUATION
      EMPFAC=0.03
C AN ESTIMATE OF ENERGY LOSS FROM MID-RANGE FREQUENCIES TO
C THE SWELL PORTION OF THE SPECTRUM
      SWLCON=35.0*G**0.5*TINC*EMPFAC
C     DEFINITION OF 4 PRIMARY PARAMETERS AFFECTING WAVE GROWTH
C     ALPSTR IS THE 'UNIVERSAL' (NOT EVEN DEPENDENT ON ANYTHING)
C     EQUILIBRIUM CONSTANT
C     BETA1 IS THE RATIO OF THE ENERGY-WEIGHTED AVERAGE PHASE SPEED
C     TO THE PHASE SPEED OF THE SPECTRAL PEAK (ALONG THE
C     MEAN DIRECTION OF PROPAGATION)
C     XLAM IS THE FACTOR WHICH CONVERTS FROM FM PARAMETERS TO TOTAL
C     ENERGY. IT'S THE RATIO OF THE TOTAL ENERGY TO THE INTEGRAL OF THE
C     EQULIBRIUM RANGE FROM FM TO INFINITY.
C     CONJON IS THE JONSWAP CONSTANT FOR USTAR NORMALIZATION
C     REPRESENTED BY 'J' IN EQUATIONS OF DOC RPT. AND M1 IN RESIO & PERRIE.
C     CDFAC IS THE ASSUMED COEFFICIENT OF DRAG FOR THE JONSWAP DATA SET
C     USED IN CONVERTING FROM PARAMETERS NORMALIZED BY U TO USTAR
C     SIGB IS THE JONSWAP PARAMETER FOR F> FM
      ALPSTR=0.0395
      BETA1=0.73
      XLAM=1.76
      CDFAC=1.2E-3
      CONJON=1.7E-7/CDFAC
      SIGB=0.11
      C73=7./3.
      C37=3./7.
C     NOTE:76.49=(7/10) * 3 * (2*PI)**(10/3) / (4*PI), EQ. 2.37 DOC. RPT. Q3
      Q3=76.49*CONJON
      CONTT=Q3*BETA1/(ALPSTR*XLAM)
      ZXX= 0.199/(XLAM*ALPSTR)
C     ROFAC IS THE RATIO  OF AIR DENSITY TO WATER DENSITY
      ROFAC=.0012
C     COEF1 AND COEF2 ARE THE PARAMETERS FOR THE DRAG LAW. SET FOR 10 METERS
      COEF1=1.1
      COEF2=0.035
C SET COEFFICIENTS BASED ON INPUT ELEVATION LEVEL
C  INPLEV=1 IMPLIES 20 METER HEIGHT
      IF (INPLEV.EQ.1) COEF1=1.1
      IF (INPLEV.EQ.1) COEF2=0.027
C     ZCUT CONTROLS WHERE THE PEAK FREQ. IS FOR A FULLY-DEVELOPED PIERSON
C     MOSCOWITZ SPECTRUM
      ZCUT=0.90
      ZCHNG=0.90
      ZDIF=1.0/ZCUT - ZCHNG
C  ZDIF=0.21111
C     GAMMAX IS THE MAXIMUM ALLOWABLE GAMMA (PEAKEDNESS)
      GAMMAX=2.0
C     GAMMIN IS THE MINIMUM (FULLY-DEVELOPED) PEAKEDNESS VALUE
      GAMMIN=0.7
      GAMDIF=GAMMAX-GAMMIN
      ACONST=0.001
      WWCON=1.4E-4*BETA1/(2.*TWOPI*0.1388)


C THIS IS THE DO LOOP OVER I,J GRID POINTS AND SPANS THE
C THE REST OF THE SUBROUTINE

      DO 5000 I=1,NI
      DO 5000 J=1,NJ
C IF NOT A WATER POINT GO TO END OF MAIN DO LOOP
      IF (IBOUND(I,J).NE.1) GO TO 5000

C      PRINT *,'INSIDE 5000 LOOP I,J = ',I,J

      GAMMA=2.0
C  SUM OF ENERGY IN THE FOUR HIGHEST FREQ. BANDS. INITIALIZE TO 0
      SUMEH=0.
C CONVERT WIND SPEED FROM KNOTS TO METERS/SEC.
      U=IWND(I,J)
C IMPOSE A LOWER LIMIT OF 1 METER/SEC ON WIND SPEED
      IF(U.LT.1.0)U=1.0
C CALCULATE THE DRAG COEFFICIENT
C THIS FORMULATION IS SIMILAR TO LARGE AND POND JPO VOL.11 1981 324-336.
      CDRAG=(COEF1+COEF2*U)*0.001
C CONVERT WIND DIRECTION TO RADIANS
      WINDR=(IWNDIR(I,J)-1)*5.*RADC
C ABSOLUTE DIFFERENCE BETWEEN WIND DIRECTION AND SEA WAVE DIRECTION
      THDIF=ABS(WINDR-SWANG(I,J))
C KEEP DIFFERENCE BETWEEN 0 AND 180 DEGREES
      IF (THDIF.GT.PI) THDIF=TWOPI-THDIF
      AFACT=COS(THDIF/2.)
C SET AFACT EQUAL TO ZERO IF WIND CHANGES DIRECTION
      IF (AFACT.LT.0.1) AFACT=0.1
      USTSQ=CDRAG*U**2
      UST=SQRT(USTSQ)
C INDEX OF FREQUENCY BAND CONTAINING PEAK FREQUENCY AT I,J
      KFRQ=KFRQA(I,J)
C TOTAL ENERGY PROPAGATING EAST - WEST
      XSUM=0.
C TOTAL ENERGY PROPAGATING NORTH - SOUTH
      YSUM=0.
C     SET TEMPORARY WORK ARRAYS AND INTEGRATE TO GET EF(K)
C  NEXT TWO LOOPS CALCULATE FREQUENCY SPECTRUM FROM PREVIOUS 2-D SPECTRUM
      DO 101 K=1,NFRQ
C INITIALIZE ENERGY OF SEA PORTION OF SPECTRUM
      FSEA(K)=0.
C INITIALIZE 2-D SPECTRAL ENERGY
      SUM=0.
C  LOOP OVER ANGLE BANDS
      DO 102 IG=1,NANG
      IF (EL2(I,J,K,IG).LT.0 ) EL2(I,J,K,IG)=0
C SCALE 2-D INTEGER ENERGY VALUES
      E2(K,IG)=EL2(I,J,K,IG)
      SUM=SUM+E2(K,IG)
  102 CONTINUE
C FORM FREQUENCY SPECTRUM
      EF(K)=SUM*AINCR
      EFN(K)=EF(K)
      E(I,J,K)=EF(K)
  101 CONTINUE
C     ADD SEA SOURCE
      FMAIN = FMA(I,J)
C     COMPUTE PARAMETRIC CHANGE IN WAVENUMBER OF SPECTRAL PEAK
C PIERSON MOSKOWITZ PEAK FREQ.
      FMPM=ZCUT*G/(TWOPI*U+1.0E-10)
C  IF PEAK FREQ. HAS REACHED FULLY DEVELOPED VALUE, NO MORE DEVELOPMENT
C  NEEDED, SO GO TO COMPUTATION OF SPECTRUM.
      IF (FMAIN.LT.FMPM) GO TO 1100
C PROCEED WITH EVOLUTION OF PEAK FREQUENCY
      TM=1.0/(FMAIN+1.0E-10)
C EQ. 2.51 R&P CHANGE OF PEAK FREQUENCY OVER A TIME STEP
      TMN=(TM**C73 + CONTT*(UST/G)**1.333*TINC)**C37
      FMN=1.0/TMN
C     CHECK ON WIND FREQUENCY LIMIT FOR LOCAL SEA
C  IF PEAK FREQ. HAS NOT REACHED PM LIMIT GO ON TO FETCH LIMITATION
C  OTHERWISE SET PEAK FREQ. TO PM LIMIT.
      IF (FMN.GE.FMPM) GO TO 987
      FMN=FMPM
      TMN=1.0/FMN
C      FETCH-LIMITING ROUTINE
C INDEX OF WAVE ANGLE
  987 ITANG=SWANG(I,J)/PI8 + 1.5
      IF (ITANG.GT.8) ITANG=1
C UPSTREAM GRID POINT BASED ON WAVE DIRECTION
      IUP=I+IPOS(ITANG)
      JUP=J+JPOS(ITANG)

C GRID SPACING OR DIAGONAL OF A GRID CELL
      ALAT=((J-1)*DLAT + XLAT0)* RADC
      DXX=XGRID
      IF (ITANG.EQ.1.OR.ITANG.EQ.5) DXX=XGRID * COS(ALAT)
      IF (MOD(ITANG,2).EQ.0) DXX=SQRT(XGRID**2+(DXX*COS(ALAT))**2)

C     CALCULATE DEPTH-DEPENDENT FACTOR FOR FETCH LIMITED GROWTH
C   WAVENUMBER AS A FUNCTION OF PEAK FREQUENCY AND DEPTH
      WKMN=WKFNC(FMN,DEPNEW(I,J))
C  DEPTH DEPENDENT PHASE SPEED OF PEAK
      CMN=TWOPI*FMN/WKMN
C  PHASE SPEED OF PEAK IN DEEP WATER
      CMP=G*TMN/TWOPI
C  GROUP SPEED OF PEAK, DEPTH DEPENDENT
      CGMN=CGF(FMN,DEPNEW(I,J),CMN)
C  PERIOD WHICH WOULD RESULT FROM A GIVEN WIND SPEED AND FETCH
      TMLIM=(1.0/FMA0(IUP,JUP)**3.333 + ZXX*DXX*UST**1.333/G**2.333)
     2	 **0.3
C  IF PERIOD NOT FETCH LIMITED, CHECK ON DEPTH LIMITATION
C	PRINT *,'AFTER SORP'
      IF (TMN.LT.TMLIM) GO TO 167
C  IF FREQ. LIMITED, SET VALUES
      FMN=1.0/(TMLIM+0.00001)
      TMN=TMLIM
C  ENSEMBLE AVERAGE GROUP SPEED
      CGV=BETA1*CGMN
C  RATIO OF CELL SIZE TO DISTANCE ENERGY TRAVELS IN A TIME STEP
      XPT=DXX/(CGV*TINC)
C ENERGY MUST TRAVEL AT LEAST ONE DELTA X IN A TIME STEP OR IT WON'T
C BE RESOLVED BY THE GRID - COURANT CRITERIA
      IF (XPT.GT.1.) XPT=1.0
C WEIGHTED AVERAGE OF PEAK FREQ. AT (I,J) AND UPSTREAM POINT
      FMN2=FMA0(I,J)*(1.0-XPT) + FMA0(IUP,JUP)*XPT
C SET TO LOWER FREQ. VALUE
      IF (FMN.GT.FMN2) FMN=FMN2
C	PRINT *,'AFTER SORP'
      GO TO 167
C     COMPUTE EQUILIBRIUM SEA SPECTRUM
C SENT HERE IF PEAK HAS REACHED PM LIMIT
C SET PEAK FREQ., PERIOD, AND PHASE SPEED OF PEAK IN DEEP WATER
 1100 FMN=FMPM
      TMN=1.0/FMN
      CMP=G*TMN/TWOPI
C     CHECK ON GROWTH PAST KM*DEPTH = 0.6  LIMIT
  167 WKDLIM=0.63/(DEPNEW(I,J)+1.0E-10)
C  THIS IS THE DISPERSION RELATIONSHIP FOR ARBITRARY DEPTH
C  AND IS USED TO LIMIT THE MOVEMENT OF FM TO LOWER FREQUENCIES
C  DUE TO INCREASED WAVE WAVE INTERACTIONS IN SHALLOW WATER
C  SEE RESIO J WATERWAY, PORT, COASTAL, AND OCEAN ENG. VOL. 113, NO. 3
C  MAY 1987, PAGE 274.
      FDLIM=SQRT(G*WKDLIM*TANHH(WKDLIM*DEPNEW(I,J)))/TWOPI
C  IF THE PEAK FREQUENCY IS LOWER THAN THE LIMIT AND IT WAS GREATER
C  THAN THE LIMIT AT THE PREVIOUS TIME STEP, SET IT TO THE LIMIT.
      IF (FMN.LT.FDLIM.AND.FMAIN.GT.FDLIM) FMN=FDLIM
C  IF THE PEAK FREQ. IS GREATER THAN THE LARGEST SPECIFIED FREQ. BAND
C  GO TO THE PARAMETRIC GROWTH SECTION.
      IF (FMN.GT.FREQ(NFRQ)) GO TO 2000
      GAMMA=GAMMAX
C  DETERMINE THE VALUE OF GAMMA (JONSWAP PEAKEDNESS PARAMETER BASED ON
C  RATIO OF WAVE PHASE SPEED AT PEAK TO WIND SPEED. IF CP/U < 0.9,
C  GAMMA IS 2.0. IF CP/U > 1.25, GAMMA=0.8
       IF(G*TMN/(TWOPI*U).GT.ZCHNG)GAMMA=GAMMAX-(G*TMN/(TWOPI*U)
     2	-ZCHNG)*GAMDIF/ZDIF
      IF (GAMMA.LT.0.8) GAMMA=0.8
C  WIND DIRECTION
      IJWNDR = IWNDIR(I,J)
C  PEAK FREQ. BAND INDEX -5
      KB=KFRQ-5
C  BAND INDEX 2 IF KB < 2
      IF (KB.LT.2) KB=2
C  SOME OF THE TERMS IN EQ. 2.27 DOC.RPT.
      ECONST=ALPSTR*(UST**2*CMP)**0.333
      KK=0
C  DO LOOP TO CALCULATE THE PART OF THE SPECTRUM FROM PEAK TO LAST BAND
      DO 72 K=KB,NFRQ
      IF (FREQ(K).LT.FMN) GO TO 72
C  THE NUMBER OF FREQUENCY BANDS FROM FM TO NFREQ
      KK=KK+1
C  PART OF THE JONSWAP FORMULATION FOR F > FM
      XPON=-0.5*((FREQ(K)-FMN)/(SIGB*FMN))**2
      IF (XPON.LT.-25.) XPON=-25.
      XXX=GAMMA**EXP(XPON)
C  EQ. 2.27 DOC. RPT. WITH THE EXCLUSION OF ZJAC
      EFN(K)=ECONST*FKA(K)*XXX	  *  ZJAC(I,J,K)
      FSEA(K)=EFN(K)
C PHASE SPEED AT FREQ. K
      CCC=G/(TWOPI*FREQ(K))
C  THE WIND SOURCE TERM. SEE DISCUSSION IN R&P P.198
      ATERM=ACONST*EFN(K)*TINC * FREQ(K)  *  (U/CCC)
      BTERM=ROFAC*FREQ(K)*TINC * (U/CCC)
      SUM=0.
CALCULATE THE DIRECTIONAL SPECTRUM
      DO 73 IG=1,NANG
C  THE 2 BELOW IS NNN, SO THE COSINE OF THE ANGLE IS TO THE 2*NNN POWER
C  OR COSINE TO THE FOURTH IN THIS CASE
      ADIST=PANG(IG,IJWNDR,2)
      E2(K,IG)=E2(K,IG) + (BTERM*E2(K,IG)  +  ATERM) * ADIST
      SUM=SUM+E2(K,IG)
   73 CONTINUE
      IF (SUM.LT.1.0E-10) SUM=1.0E-10
C  RATIO OF ENERGY IN BAND K TO TOTAL
      RR=EFN(K)/(SUM*AINCR)
      IF (RR.GT.1.) RR=1.
C  CALCULATE THE EAST-WEST AND NORTH-SOUTH COMPONENTS OF ENERGY IN
C  EACH FREQ. DIRECTION BIN.
      DO 74 IG=1,NANG
      E2(K,IG)=E2(K,IG)*RR
      XSUM=XSUM+E2(K,IG)*COSR(IG)
      YSUM=YSUM+E2(K,IG)*SINR(IG)
   74 CONTINUE
   72 CONTINUE
C     REDEFINE SWANG AS MEAN SEA ANGLE
      SWANG0=SWANG(I,J)
      SWANG(I,J) = ATAN2(YSUM,XSUM+EPSS)
      IF (SWANG(I,J).LT.0.) SWANG(I,J)=SWANG(I,J)+TWOPI
      X=COS(SWANG0) + COS(SWANG(I,J))
      Y=SIN(SWANG0) + SIN(SWANG(I,J))
C  AVERAGE ANGLE OVER THIS AND LAST TIME STEP
      ANGAV=ATAN2(Y,X)
      IF (ANGAV.LT.0.) ANGAV=ANGAV+TWOPI
C     REDEFINE KFRQ AS NEW COUNTER FOR SPECTRAL PEAK
C  BAND NUMBER OF THE LARGEST FREQ. ON THE FRONT FACE
      KFRQ=NFRQ-KK+1
C     ADD FRONT FACE OF LOCAL SEA
C  EQ. 2.27 DOC. RPT. FOR F < FM
      EE2=GAMMA*AFACT*ALPSTR*(UST**2*CMP)**0.333*G/(TWOPI**3*FMN**4)
      FFCE=1.15
      KM1=KFRQ-1
      IF (KM1.GT.1) GO TO 7777
      PRINT *,'   W A R N I N G !!!!!! '
      PRINT *,'SPECTRAL PEAK HAS REACHED FIRST FREQUENCY *****'
      PRINT *,'GRID POINT I = ',I,' J = ',J
      WRITE (15,*) ' SPECTRAL PEAK IN 1ST FREQ AT I,J = ',I,J
 7777 IJWVDR = ANGAV/PI72 + 1.01
      IF (IJWVDR.GT.72) IJWVDR=1
C  PART OF THE ENERGY EQ. FOR THE CHANGE OF ENERGY WITH TIME
C  SEE R&P EQ. 2.15
      SCALE=WWCON*UST*UST/(FMN*FMN)  *	TINC
C  THIS LOOP ADDS THE FRONT FACE OF THE SPECTRUM FROM 1/2 FM TO THE PEAK
C  AND DETERMINES THE TIME RATE OF CHANGE OF ENERGY IN EACH BAND
      DO 81 K=1,KM1
      IF (FREQ(K).LT.0.5*FMN) GO TO 81
      XXX=FMN/FREQ(K)
      YYY= FFCE-FFCE*XXX**4
      IF (YYY.LT.-20.) YYY=-20.
      PSI=EXP(YYY)
      F2=EE2*EXP(YYY) *ZJAC(I,J,K)
      FSEA(K)=F2
      SUM=0.
C  CHANGE IN ENERGY BY FREQ. BAND FOR PRESENT TIME MINUS PREVIOUS TIME
      DFDT=F2-EF(K)
C  IF ENERGY DECREASING SET TO 0
      IF (DFDT.LT.0.) DFDT=0.
C  ADD INCREASE IN ENERGY TO SPECTRUM
      EFN(K)=EF(K)+DFDT
C  DETERMINE 2-D CHANGE IN ENERGY AND ADD TO 2-D SPECTRUM
      DFDT2=SCALE*PSI
      DO 82 IG=1,NANG
      E2(K,IG)=E2(K,IG)+DFDT2*PANG(IG,IJWVDR,2)
      SUM=SUM+E2(K,IG)
   82 CONTINUE
      IF (SUM.LT.1.0E-10) SUM=1.0E-10
      RR=EFN(K)/(SUM*AINCR)
      DO 83 IG=1,NANG
      E2(K,IG)=E2(K,IG)*RR
   83 CONTINUE
   81 CONTINUE
C     REDEFINE WAVE ANGLE COUNTER FOR NEXT TIME STEP
      IWVDIR(I,J)=SWANG(I,J)/PI72  +  1.5
      IF (IWVDIR(I,J).GT.72) IWVDIR(I,J)=1
C  BYPASS PARAMETRIC GROWTH SECTION
      GO TO 3000
C     PARAMETRIC GROWTH ONLY
C  ARRIVE HERE BECAUSE PEAK FREQ. IS LARGER THAN MAX DISCRETE FREQ.
 2000 CONTINUE
      GAMMA=GAMMAX
      IF (CMP/U.GT.ZCHNG) GAMMA=GAMMAX-(CMP/U-ZCHNG)*GAMDIF/ZDIF
      IF (GAMMA.LT.0.8) GAMMA=0.8
      SCALF=GAMMA*ALPSTR*(UST**2*CMP)**0.333*G/(TWOPI**3*FMN**4)
      KFRQ=NFRQ+1
      FFCE=1.00
      IJWNDR = IWNDIR(I,J)
      NF1=NFRQ-4
C  CALCULATE SPECTRUM FROM BAND NFRQ-4 TO NFRQ
C  CHECK TO SEE IF IT IS GROWING IN TIME. CALCULATE NEW 2-D SPECTRUM
      DO 85 K=NF1,NFRQ
      IF (FREQ(K).LT.0.5*FMN) GO TO 85
      XXX=FMN/FREQ(K)
      YYY= FFCE-FFCE*XXX**4
      IF (YYY.LT.-20.) YYY=-20.
      F2=SCALF*EXP(YYY)  *  ZJAC(I,J,K)
      FSEA(K)=F2
      DELFK=F2-EF(K)
      IF (DELFK.LT.0.) DELFK=0.
      EFN(K)=EF(K)+DELFK
      SUMEH=SUMEH+ F2*DELF(K)
      IF (DELFK.LT.0.) DELFK=0.0051*EF(K)*FMN*TINC
      DFDT2=DELFK/AINCR
      SUM=0.
      DO 86 IG=1,NANG
      E2(K,IG)=E2(K,IG)+DFDT2*PANG(IG,IJWNDR,2)
      SUM=SUM+E2(K,IG)
   86 CONTINUE
      IF (SUM.LT.1.0E-10) SUM=1.0E-10
      RR=EFN(K)/(SUM*AINCR)
      DO 87 IG=1,NANG
      E2(K,IG)=E2(K,IG)*RR
   87 CONTINUE
   85 CONTINUE
      SWANG(I,J)=RADC*(IWNDIR(I,J)-1.)*5.
C  SET WAVE DIRECTION TO WIND DIRECTION FOR PARAMETRIC GROWTH
      IWVDIR(I,J)=IWNDIR(I,J)
C     END OF LOCAL SEA SOURCE
C     SWELL SOURCE TREATMENT
C  ARRIVE HERE IF BYPASSING PARAMETRIC GROWTH
 3000 CONTINUE
      ESWL=0.
      SWMAX=0.001
      KMAX=1
C  SCAN SPECTRUM FROM BAND 1 TO 0.7*FM TO FIND BAND WITH LARGEST ENERGY
      DO 3001 K=1,NFRQ
      IF (FREQ(K).GT.0.7*FMN) GO TO 3001
      ESWL=ESWL+EFN(K)*DELF(K)
      IF (EFN(K).LT.SWMAX) GO TO 3001
      SWMAX=EFN(K)
      KMAX=K
 3001 CONTINUE
C  SKIP CALCULATING SWELL PARAMETERS IF BAND IS > NFRQ OR 1 OR ENERGY < 0.01
      IF (KMAX.GE.NFRQ) GO TO 3010
      IF (KMAX.LE.1) GO TO 3010
      IF (ESWL.LT.0.01) GO TO 3010
C  SWELL PEAK FREQUENCY
      SWLPK=FREQ(KMAX)
C  SWELL PEAK PERIOD
      TSWL(I,J)=1.0/SWLPK
C  SWELL WAVE NUMBER
      SWLKM=(TWOPI*SWLPK)**2/G
C  LOSS OF ENERGY FROM MID FREQS. TO SWELL PORTION OF SPECTRUM
C  SEE EQ.2.42 DOC.RPT.
      ELSS=SWLCON*SWLKM**4.5  * ESWL**3
C  DETERMINE ENERGY BETWEEN SWELL PEAK AND 0.7 * FM
      EBACK=0
      DO 3030 K=KMAX,NFRQ
      IF (FREQ(K).GT.0.7*FMN) GO TO 3030
      EBACK=EBACK+EFN(K) *DELF(K)
 3030 CONTINUE
C  LIMIT ENERGY LOSS TO 20% OF EBACK
      IF (ELSS.GT.0.2*EBACK) ELSS=0.2*EBACK
      RR=(EBACK-ELSS)/EBACK
      IF (RR.LT.0.8) RR=0.8
      ELSS=EBACK*(1.0-RR)
      EFF=0.2*ELSS
C  REDUCE FREQ. AND 2-D SPECTRUM FROM SWELL PEAK TO 0.7*FM BY ENERGY LOSS
      DO 3002 K=KMAX,NFRQ
      IF (FREQ(K).GT.0.7*FMN) GO TO 3002
      EFN(K)=EFN(K)*RR
      DO 3003 IG=1,NANG
      E2(K,IG)=E2(K,IG)*RR
 3003 CONTINUE
 3002 CONTINUE
C  TRANSFER ENERGY TO BANDS 1 TO SWELL PEAK MINUS 1
      EFRONT=0.
      DO 3040 K=1,KMAX-1
      EFRONT=EFRONT+EFN(K)*DELF(K)
 3040 CONTINUE
      RR=(EFRONT+EFF)/(EFRONT+1.0E-10)
      DO 3004 K=1,KMAX-1
      EFN(K)=EFN(K)*RR
      DO 3005 IG=1,NANG
      E2(K,IG)=E2(K,IG)*RR
 3005 CONTINUE
 3004 CONTINUE
C  DETERMINE SWELL HEIGHT
      ESWL=EFRONT+EBACK
      IF (ESWL.LT.0.) ESWL=0.
      HSWL(I,J)=4.0*SQRT(ESWL)
C  DETERMINE SWELL DIRECTION
      XCMP=0.
      YCMP=0.
      DO 3006 IG=1,NANG
      XCMP=XCMP+E2(KMAX,IG)*COSR(IG)
      YCMP=YCMP+E2(KMAX,IG)*SINR(IG)
 3006 CONTINUE
      SWDR=ATAN2(YCMP,XCMP+EPSS)
      IF (SWDR.LT.0.) SWDR=SWDR+TWOPI
      SWLDIR(I,J)=SWDR*RADDEG
C  FINISHED WITH SWELL CALCULATIONS
      GO TO 6000
C  ARRIVE HERE IF NO SIGNIFICANT SWELL
 3010 CONTINUE
      HSWL(I,J)=0.
      TSWL(I,J)=0.
      SWLDIR(I,J)=0.
C     RESET ALL VALUES INTO LARGE-PERMANENT ARRAYS
 6000 CONTINUE
      FSUM=0.
      SUM=0.
      SUMSW=0.
      YSW=0.
      XSW=0.
      KMAX=0
C  LOWER LIMIT ON SWELL ENERGY
      SWMAX=1.
C  LOOP OVER ENTIRE SPECTRUM
      DO 90 K=1,NFRQ
C  IF ENERGY IN BAND K IS < A SMALL VALUE GO TO END OF LOOP
      IF (EFN(K).LT.0.000001) GO TO 990
      E(I,J,K)=EFN(K)
      SUM=SUM+E(I,J,K)*DELF(K)
C  MEAN FREQUENCY
      FSUM=FSUM+E(I,J,K)*FREQ(K)*DELF(K)
C  BYPASS SWELL CALCULATIONS IF K > PEAK FREQ.
      IF (K.GE.KFRQ) GO TO 199
C  ESTIMATE SWELL ENERGY AS TOTAL - SEA FOR BAND K
      ENSWL=EFN(K)-FSEA(K)
C  DETERMINE PEAK FREQ. FOR SWELL
      IF (ENSWL.LE.SWMAX)GO TO 198
      KMAX=K
      SWMAX=ENSWL
  198 IF (ENSWL.LT.0) ENSWL=0.
      SUMSW=SUMSW+ENSWL*DELF(K)
C
  199 DO 190 IG=1,NANG
      EL2(I,J,K,IG)=E2(K,IG)
C  BYPASS IF K > FM
      IF (K.GE.KFRQ) GO TO 190
C  ADD FRONT FACE COMPONENT ENERGIES ONTO REAR FACE CALCULATED IN DO 72 LOOP
      XSUM=XSUM+E2(K,IG)*COSR(IG)
      YSUM=YSUM+E2(K,IG)*SINR(IG)
C  CALCULATE COMPONENTS FOR SWELL
      YSW=YSW+E2(K,IG)*SINR(IG)
      XSW=XSW+E2(K,IG)*COSR(IG)
  190 CONTINUE
      GO TO 90
C  ARRIVE HERE IF NEGLIGIBLE ENERGY AND SET 1-D AND 2-D SPECRTA TO 0
  990 E(I,J,K)=0.
      DO 991 IG=1,NANG
      EL2(I,J,K,IG)=0
  991 CONTINUE
   90 CONTINUE
C  SUM OF ENERGY FROM K=1,NFRQ MINUS ENERGY IN 4 HIGHEST BANDS
      SUM=SUM-SUMEH
      IF (SUM.LT.0.) SUM=0.
C  EQ. 2.26 IN DOC. RPT. TOTAL ENERGY IN SPECTRUM EXPRESSED IN TERMS OF FM
      EZRO =XLAM*ALPSTR*(UST**2*CMP)**0.333*G/(3.*TWOPI**3*FMN**3)
      EHIGH=EZRO
C  SAME AS ABOVE DIVIDED BY XLAM. TOTAL ENERGY FROM FM TO INFINITY
      IF (FMN.LT.FREQ(NFRQ)) EHIGH=ALPSTR*(UST**2*CMP)**0.333*G/
     2	 (3.*TWOPI**3*FMH**3)
      HSCALE(I,J)=EHIGH
      HSIG(I,J)=4.0*SQRT(SUM+EHIGH)
C  MEAN FREQ.  ENERGY WEIGHTED
      FSUM=FSUM/(SUM+1.0e-10)
      IF (FMN.GT.FREQ(NFRQ-1)) FSUM=1.4*FMN
C  MEAN PERIOD
      TSS(I,J)=1.0/(FSUM+1.0e-10)
C  SWELL HEIGHT AND PERIOD AND DIRECTION
      IF (KMAX.GT.0) TSWL(I,J)=1.0/FREQ(KMAX)
      IF (KMAX.EQ.0) TSWL(I,J)=0.
      HSWL(I,J)=4.0*SQRT(SUMSW)
      IF (HSWL(I,J).LT.0.01) HSWL(I,J)=0.
      SWLDIR(I,J)=ATAN2(YSW,XSW+EPSS)
      IF (SWLDIR(I,J).LT.0.) SWLDIR(I,J)=SWLDIR(I,J)+TWOPI
C  MEAN DIRECTION FOR K=1 TO NFRQ
      AVANG(I,J)=ATAN2(YSUM,XSUM+EPSS)
      IF (AVANG(I,J).LT.0.) AVANG(I,J)=AVANG(I,J)+TWOPI
C  IF ENERGY IS SMALL MEAN DIRECTION IS SEA WAVE ANGLE
      IF (XSUM+YSUM.LT.0.01) AVANG(I,J)=SWANG(I,J)
      IF (HSIG(I,J).LT.0.63*DEPNEW(I,J)) GO TO 678
C  LIMIT ENERGY VALUES AND HEIGHTS BY DEPTH
      PRINT *,'**** LIMITING HS FOR I,J,DEP ',I,J,DEPNEW(I,J)
      ERAT=(0.63*DEPNEW(I,J)/HSIG(I,J))**2
      DO 679 K=1,NFRQ
      DO 679 IG=1,NANG
      EL2(I,J,K,IG)=EL2(I,J,K,IG)*ERAT
  679 CONTINUE
      HSIG(I,J)=0.63*DEPNEW(I,J)
      HSWL(I,J)=SQRT(ERAT)*HSWL(I,J)
C  ARRIVE HERE IF NOT DEPTH LIMITED
C  FILL FIELDS OF PEAK FREQUENCY, PEAK FREQUENCY BAND INDEX, GAMMA, ALPHA
  678 FMA(I,J)=FMN
      KFRQA(I,J)=KFRQ
      GAMOUT(I,J)=GAMMA
	ALPHA=ALPSTR
      ALPOUT(I,J)=ALPHA
C  END OF I,J LOOP
 5000 CONTINUE
      RETURN
      END
C
C
C
       FUNCTION TANHH(X)
       IF(X.GT.10.)X=10.
       TANHH=TANH(X)
       RETURN
       END
C
C
C
      FUNCTION WKFNC(FM,D)
      WK0=(6.283185*FM)**2/9.8
      WK=WK0
      IF (WK*D.LT.0.1) GO TO 3
      IF (WK*D.GT.10.) GO TO 4
    1 WKN=WK0/TANH(WK*D)
      IF (ABS(WK-WKN).LT.0.000001) GO TO    2
      WK=WKN
      GO TO 1
    2 WKFNC=WKN
      RETURN
    3 WL=SQRT(9.8*D)/FM
      WKFNC=6.283185/WL
      RETURN
    4 WL=1.56/(FM*FM)
      WKFNC=6.283185/WL
      RETURN
      END
C
C
C
      FUNCTION CGF(F,D,C)
      WL=C/F
      WK=6.283/WL
      TKH=2.0*WK*D
      IF (TKH.GT.20.) TKH=20.
      CGF=0.5*C*(1.0+TKH/SINH(TKH))
      RETURN
      END



      SUBROUTINE INITSPR
      INCLUDE 'PARAMS'
      INCLUDE 'PARAMS2'
C     PARAMETER (IDMN=35,JDMN=51,IF=15,IA=16,NOBPTS=11,NBN=11)
      COMMON /SMOD/ DPX(IDMN2,JDMN2,IF2,IA2), DPY(IDMN2,JDMN2,IF2,IA2),
     2	   DTH(IDMN2,JDMN2,IF2,IA2),CCGRJX(IDMN2,JDMN2,IF2,IA2)
     +	   ,CCGRJY(IDMN2,
C    3	   JDMN,IF,IA)
     3	   JDMN2,IF2,IA2),LLX(IDMN2,JDMN2,IF2,IA2),LHX(IDMN2,JDMN2
     +	   ,IF2,IA2),
     4	   X1X(IDMN2,JDMN2,IF2,IA2),LLY(IDMN2,JDMN2,IF2,IA2),
     5	   LHY(IDMN2,JDMN2,IF2,IA2), X1Y(IDMN2,JDMN2,IF2,IA2)
      COMMON /CCG/ CA(IDMN,JDMN,IF),CGA(IDMN,JDMN,IF),DELX(JDMN)
      COMMON /A3/ G,PI,XN,DL,NANG,NFRQ,NX,NY,TWOPI,FAC,DT,RADDEG
      DIMENSION DCDX(IDMN2,JDMN2),DCDY(IDMN2,JDMN2)
      COMMON /O5/ IBOUND(IDMN,JDMN)
      COMMON /ANGS/ ANG(IA),SINA(IA),COSA(IA),ISLP(IDMN,JDMN)
      COMMON /DEPTH/ DEP(IDMN,JDMN)

      INTEGER*2 DPX,DPY,DTH,CCGRJX,CCGRJY,LLX,LHX,X1X,LLY,LHY,X1Y

      COMMON /BAND/ IUA(IA),JUA(IA)


	EP=0.0000001
      AINC=TWOPI/NANG
      ALIM=2*AINC
      DO 80 L=1,NANG
      ANG(L)=(L-1)*AINC
      SINA(L)=SIN(ANG(L))
      COSA(L)=COS(ANG(L))
   80 CONTINUE

C   CALCULATE DATA VALUES FOR IUA AND JUA
C   USING A VARIABLE VALUE FOR NANG
	DO 86 I=1,NANG
	X=COSA(I)
	Y=SINA(I)
	IUA(I)=-SIGN(1.0,X)
	JUA(I)=-SIGN(1.0,Y)
	IF (ABS(X).LT.0.0001) IUA(I)=0.0
	IF (ABS(Y).LT.0.0001) JUA(I)=0.0
  86	CONTINUE
	PRINT *, IUA
	PRINT *,JUA
      CURVLIM=1./(2.*DL)
C     PRINT *,DL,CURVLIM,PI,NANG,NX,NY,TWOPI,DT

      DO 50 K=1,NFRQ

      DO 60 I=1,NX
      DO 60 J=1,NY
      DO 60 L=1,NANG
      X1X(I,J,K,L)=0
      DPX(I,J,K,L)=0
      DPY(I,J,K,L)=0
      DTH(I,J,K,L)=0
      CCGRJX(I,J,K,L)=128
      CCGRJY(I,J,K,L)=128
   60 CONTINUE


      DO 1 I=2,NX-1
      DO 1 J=2,NY-1

      DCDX(I,J)=(CA(I+1,J,K)-CA(I-1,J,K))/(2.*DELX(J))
      DCDY(I,J)=(CA(I,J+1,K)-CA(I,J-1,K))/(2.*DL)

C     IF (J.EQ.9.AND.I.EQ. 4) PRINT*,'A ',I,J,DCGDX(I,J),DCGDY(I,J)
C    2			,DCDX(I,J),DCDY(I,J)

    1 CONTINUE

      DO 2 I=2,NX-1
      DO 2 J=2,NY-1

      IF (IBOUND(I,J).NE.1) GO TO 2
      DHX=DEP(I-1,J)-DEP(I+1,J)
      DHY=DEP(I,J-1)-DEP(I,J+1)

C     ANGSL IS THE SLOPE ANGLE (POSITIVE TOWARD SHALLOW WATER)

      ANGSL=ATAN2(DHY,DHX+0.0000001)
      ISLP(I,J)=0
      IF (ABS(COS(ANGSL)).GT.ABS(SIN(ANGSL))) ISLP(I,J)=1

      DO 3 L=1,NANG

      LL=L+NANG/2

C     NOTE LL REFERENCES THE "BACKTRACKING" RAY

      IF (LL.GT.NANG) LL=LL-NANG
C      CGAV=CGA(I,J,K)+ 0.5*DCGDX(I,J)*CGA(I,J,K)*COSA(LL)*DT
C     2		   + 0.5*DCGDY(I,J)*CGA(I,J,K)*SINA(LL)*DT
      DTHDS=(DCDX(I,J)*SINA(LL)-DCDY(I,J)*COSA(LL))/(CA(I,J,K)+EP)
      IF (ABS(DTHDS).GT.CURVLIM) DTHDS=CURVLIM*DTHDS/ABS(DTHDS)
C      TDPX=DELS*COS(THTAV)/PHI
C      TDPY=DELS*SIN(THTAV)/PHI
      IDPX=ABS(DPX(I,J,K,L))/DELX(J) * 255
      IDPY=ABS(DPY(I,J,K,L))/DL * 255
      IF (IDPX.GT.255) IDPX=255
      IF (IDPY.GT.255) IDPY=255
      DPX(I,J,K,L)=IDPX
      DPY(I,J,K,L)=IDPY

      II=I+IUA(L)
      CGUPX=CGA(II,J,K)
      CGRX=CGUPX/CGA(I,J,K)
      JJ=J+JUA(L)
      CGUPY=CGA(I,JJ,K)
      CGRY=CGUPY/CGA(I,J,K)

C     SOLVE FOR INTERSECTION WITH BOUNDARY

      C1=CA(I,J,K)
      C2=CA(II,J,K)
      RX=C1/C2
      C3=CA(I,JJ,K)
      RY=C1/C3

      ANGZ=ANG(L)-ANGSL
      IF (ANGZ.GT.PI) ANGZ=ANGZ-TWOPI
      IF (ANGZ.LT.-PI) ANGZ=ANGZ+TWOPI
      IF (ABS(ANGZ).GT.PI/2.) RX=1./RX
      IF (ABS(ANGZ).GT.PI/2.) RY=1./RY
      IF (ANGZ.GT.PI/2.) ANGZ=PI - ANGZ
      IF (ANGZ.LT.-PI/2.) ANGZ=-PI - ANGZ

      SINZ=SIN(ANGZ)

      IF (RX.GT.2.0) RX=2.0
      IF (RX.GT.1.) GO TO 76

      ANG0=ASIN(RX *SINZ)
      BX=ABS(COS(ANGZ)/(COS(ANG0)+EP))
C  SET BX AND BY FOR 90 OR 270 DEG INCIDENCE
	IF(ABS(SINZ).GT.0.99985) BX=1.0
      IF (ISLP(I,J).EQ.1) ZDTH=ANG0-ANGZ


      GO TO 61
   76 ANG0=ASIN(SINZ/RX)
      BX=ABS(COS(ANG0)/(COS(ANGZ)+EP))
      IF (ISLP(I,J).EQ.1) ZDTH=ANGZ-ANG0



   61 IF (RY.GT.2.0) RY=2.0
      IF (RY.GT.1.) GO TO 62
      ANG0=ASIN(RY*SINZ)
      BY=ABS(COS(ANGZ)/(COS(ANG0)+EP))
	IF(ABS(SINZ).GT.0.99985) BY=1.0
      IF (ISLP(I,J).NE.1) ZDTH=ANG0-ANGZ
      GO TO 63

   62 ANG0=ASIN(SINZ/RY)
      BY=ABS(COS(ANG0)/(COS(ANGZ)+EP))
	IF(ABS(SINZ).GT.0.99985) BY=1.0
      IF (ISLP(I,J).NE.1) ZDTH=ANGZ-ANG0


   63 IF (ZDTH.GT.ALIM) ZDTH=ALIM
      IF (ZDTH.LT.-ALIM) ZDTH=-ALIM
      IDTH=128+ZDTH*126/ALIM
      DTH(I,J,K,L)=IDTH

      IF (BX.GT.1.) BX = 1./CGRX
      IF (BY.GT.1.) BY = 1./CGRY

      ITMPX=CGRX*BX * 128
      IF (ITMPX.GT.255) ITMPX=255
      CCGRJX(I,J,K,L)=ITMPX
      ITMPY=CGRY*BY * 128
      IF (ITMPY.GT.255) ITMPY=255
      CCGRJY(I,J,K,L)=ITMPY

    3 CONTINUE

    2 CONTINUE

   50 CONTINUE

C  PRECALCULATE INDICIES FOR USE IN SPROP LOOP 7
C  ADDED BY JMH FEB.1991
       DO 4 K=1,NFRQ
       DO 4 I=2,NX-1
       DO 4 J=2,NY-1
       DO 4 L=1,NANG
       ANGZ=ANG(L)+(DTH(I,J,K,L)-128)*(ALIM/126.)*ISLP(I,J)
       IF(ABS(ANGZ).LT.0.02) ANGZ=0.0
       IF(ANGZ.GT.TWOPI) ANGZ=ANGZ-TWOPI
       IF(ANGZ.LT.0.0) ANGZ=ANGZ+TWOPI
       LLX(I,J,K,L)=ANGZ/AINC+1
       IF(LLX(I,J,K,L).LT.1) ANGZ=ANGZ+TWOPI
       IF(LLX(I,J,K,L).LT.1) LLX(I,J,K,L)=LLX(I,J,K,L)+NANG
       IF(LLX(I,J,K,L).GT.NANG) LLX(I,J,K,L)=LLX(I,J,K,L)-NANG
       LHX(I,J,K,L)=LLX(I,J,K,L)+1
       IF(LHX(I,J,K,L).GT.NANG) LHX(I,J,K,L)=LHX(I,J,K,L)-NANG
       IAX=(ANGZ-(LLX(I,J,K,L)-1)*AINC)/AINC * 128
       IF (IAX.GT.255) IAX=255
       IF (IAX.LT.1) IAX=1
       X1X(I,J,K,L)=IAX
    4  CONTINUE

C  PRECALCULATE INDICIES FOR USE IN SPROP LOOP 27
C  ADDED BY JMH FEB.1991
       DO 5 K=1,NFRQ
       DO 5 I=2,NX-1
       DO 5 J=2,NY-1
       DO 5 L=1,NANG
       ANGZ=ANG(L)+(DTH(I,J,K,L)-128)*(ALIM/126.)*(1-ISLP(I,J))
       IF(ABS(ANGZ).LT.0.02) ANGZ=0.0
       IF(ANGZ.GT.TWOPI) ANGZ=ANGZ-TWOPI
       IF(ANGZ.LT.0.0) ANGZ=ANGZ+TWOPI
       LLY(I,J,K,L)=ANGZ/AINC+1
       IF(LLY(I,J,K,L).LT.1) LLY(I,J,K,L)=LLY(I,J,K,L)+NANG
       IF(LLY(I,J,K,L).GT.NANG) LLY(I,J,K,L)=LLY(I,J,K,L)-NANG
       LHY(I,J,K,L)=LLY(I,J,K,L)+1
       IF(LHY(I,J,K,L).GT.NANG) LHY(I,J,K,L)=LHY(I,J,K,L)-NANG
       IAY=(ANGZ-(LLY(I,J,K,L)-1)*AINC)/AINC * 128
       IF (IAY.GT.255) IAY=255
       IF (IAY.LT.1) IAY=1
       X1Y(I,J,K,L)=IAY
    5  CONTINUE

C     PRINT *,'DPX SAMPLES'
      DO 500 J=NY,2,-1
C     PRINT *,'J = ',J
  180 FORMAT (8F7.4)
 500	 CONTINUE

C     PRINT *,'DPY SAMPLES'
      DO 501 J=NY,2,-1
C     PRINT *,'J = ',J
C     PRINT 180,(DPY( 8,J,1,L),L=1,NANG)
  501 CONTINUE

C     PRINT *,'DTH SAMPLES'
C     DO 502 J=NY,2,-1
C     PRINT *,'J = ',J
C     PRINT 180,(DTH( 8,J,1,L),L=1,NANG)

C     PRINT *,'CCGRJX  SAMPLES'
C     DO 503 J=NY,2,-1
C     PRINT *,'J = ',J
C     PRINT 180,(CCGRJX( 8,J,1,L),L=1,NANG)
  503 CONTINUE

C     PRINT *,'CCGRJY SAMPLES'
C     DO 505 J=NY,1,-1
C     PRINT *,'J = ',J
C     PRINT 180,(CCGRJY( 8,J,1,L),L=1,NANG)
  505 CONTINUE

      RETURN
      END




      SUBROUTINE SPROP(K)
C     SUBROUTINE ADAPTED TO ALLOW REFRACTION AND SHOALING ESTIMATES INTO
C     ADVECTION IN DWAVE

      INCLUDE 'PARAMS'
      INCLUDE 'PARAMS2'

C     PARAMETER (IDMN=105,JDMN=25,IF=15,IA=16,NOBPTS=5,NBN=5)
C    PARAMETER (IDMN2=15,JDMN2=51,IF2=15,IA2=16,NOBPTS2=11,NBN2=11)
      COMMON /SMOD/ DPX(IDMN2,JDMN2,IF2,IA2), DPY(IDMN2,JDMN2,IF2,IA2),
     2	   DTH(IDMN2,JDMN2,IF2,IA2),CCGRJX(IDMN2,JDMN2,IF2,IA2)
     +	   ,CCGRJY(IDMN2,
C    3	   JDMN,IF,IA)
     3	   JDMN2,IF2,IA2),LLX(IDMN2,JDMN2,IF2,IA2),LHX(IDMN2,
     +	   JDMN2,IF2,IA2),
     4 X1X(IDMN2,JDMN2,IF2,IA2),LLY(IDMN2,JDMN2,IF2,IA2),
     5LHY(IDMN2,JDMN2,IF2,IA2), X1Y(IDMN2,JDMN2,IF2,IA2)
      COMMON /CCG/ CA(IDMN,JDMN,IF),CGA(IDMN,JDMN,IF),DELX(JDMN)
      COMMON /A3/ G,PI,XN,DL,NANG,NFRQ,NX,NY,TWOPI,FAC,DT,RADDEG
      DIMENSION EOUTX(IDMN2,JDMN2,IA2),EOUTY(IDMN2,JDMN2,IA2),
     2	    EINX(IDMN2,JDMN2,IA2),
     2	EINY(IDMN,JDMN,IA),EINFX(IDMN,JDMN,IA),EINFY(IDMN,JDMN,IA)
      COMMON /ANGS/ ANG(IA),SINA(IA),COSA(IA),ISLP(IDMN,JDMN)
      COMMON /CEN/ E(IDMN,JDMN,IF,IA),ENXT(NBN,IF,IA),ELST(NBN,IF,IA)
      COMMON /O5/ IBOUND(IDMN,JDMN)
	COMMON /BAND/ IUA(IA),JUA(IA)

      INTEGER*2 DPX,DPY,DTH,CCGRJX,CCGRJY,LLX,LHX,X1X,LLY,LHY,X1Y


C     PRINT *,'ENTERING X-PROP'
C     X-PROPAGATION


      DO 1 I=1,NX
      DO 1 J=1,NY

C      IF (IBOUND(I,J).NE.1) GO TO 1

C      IF (IBOUND(I,J).EQ.0) GO TO 1

      DO 2 L=1,NANG
      EOUTX(I,J,L)=E(I,J,K,L)*DPX(I,J,K,L)/255.
      II=I+IUA(L)
      EINX(I,J,L)=E(II,J,K,L)

    2 CONTINUE


    1 CONTINUE


C     PRINT *,'ENTERING THETA-PROP - X'
C     THETA PROPAGATION

      DO 5 I=2,NX-1
      DO 5 J=2,NY-1

      DO 6 L=1,NANG
      EINFX(I,J,L)=0.
    6 CONTINUE

      IF (IBOUND(I,J).EQ.0) GO TO 5

C LOOP 7 MODIFIED BY JMH FEB.1991 TO VECTORIZE
      DO 7 L=1,NANG
      EINFX(I,J,LLX(I,J,K,L))=EINFX(I,J,LLX(I,J,K,L))+CCGRJX(I,J,K,L)/128.
     1 *EINX(I,J,L)*(1.-X1X(I,J,K,L)/255.)
     2 * DPX(I,J,K,LLX(I,J,K,L))/255.
      EINFX(I,J,LHX(I,J,K,L))=EINFX(I,J,LHX(I,J,K,L))+CCGRJX(I,J,K,L)/128.
     1 *EINX(I,J,L)*X1X(I,J,K,L)/255.
     2 * DPX(I,J,K,LHX(I,J,K,L))/255.
C     EINFX(I,J,LL)=EINFX(I,J,LL)+CCGRJX(I,J,K,L)/128.*EINX(I,J,L)*(1.-X1)
C    2	  * DPX(I,J,K,LL)/255.
C     EINFX(I,J,LH)=EINFX(I,J,LH)+CCGRJX(I,J,K,L)/128.*EINX(I,J,L)*X1
C    2	  * DPX(I,J,K,LH)/255.

    7 CONTINUE

      IF (IBOUND(I,J).NE.1) GO TO 5

      DO 8 L=1,NANG
      E(I,J,K,L)=E(I,J,K,L)-EOUTX(I,J,L)   +   EINFX(I,J,L)
      IF (E(I,J,K,L).LT.0.) E(I,J,K,L)=0.
    8 CONTINUE

    5 CONTINUE

C     Y-PROP
C     PRINT *,'ENTERING Y-PROP'

      DO 17 I=1,NX
      DO 17 J=1,NY

C      IF (IBOUND(I,J).NE.1) GO TO 17

      IF (IBOUND(I,J).EQ.0) GO TO 17


      DO 9 L=1,NANG
      EOUTY(I,J,L)=E(I,J,K,L)*DPY(I,J,K,L)/255.
      JJ=J+JUA(L)
      EINY(I,J,L)=E(I,JJ,K,L)

    9 CONTINUE


   17 CONTINUE


C     PRINT *,'ENTERING THETA-PROP - Y'
C     THETA PROPAGATION

      DO 3 I=2,NX-1
      DO 3 J=2,NY-1

      DO 16 L=1,NANG
      EINFY(I,J,L)=0.
   16 CONTINUE

      IF (IBOUND(I,J).EQ.0) GO TO 3

C  LOOP 27 MODIFIED BY JMH FEB.1991 TO VECTORIZE
      DO 27 L=1,NANG
      EINFY(I,J,LLY(I,J,K,L))=EINFY(I,J,LLY(I,J,K,L))+CCGRJY(I,J,K,L)/128.
     1	*EINY(I,J,L)*(1.-X1Y(I,J,K,L)/255.)
     2	* DPY(I,J,K,LLY(I,J,K,L))/255.
      EINFY(I,J,LHY(I,J,K,L))=EINFY(I,J,LHY(I,J,K,L))+CCGRJY(I,J,K,L)/128.
     1	*EINY(I,J,L)*X1Y(I,J,K,L)/255.
     2	* DPY(I,J,K,LHY(I,J,K,L))/255.
C     EINFY(I,J,LL)=EINFY(I,J,LL)+CCGRJY(I,J,K,L)/128.*EINY(I,J,L)*(1.-X1)
C    2	  * DPY(I,J,K,LL)/255.
C     EINFY(I,J,LH)=EINFY(I,J,LH)+CCGRJY(I,J,K,L)/128.*EINY(I,J,L)*X1
C    2	  * DPY(I,J,K,LH)/255.

   27 CONTINUE

      IF (IBOUND(I,J).NE.1) GO TO 3

      DO 18 L=1,NANG
      E(I,J,K,L)=E(I,J,K,L)-EOUTY(I,J,L)   +   EINFY(I,J,L)
      IF (E(I,J,K,L).LT.0.) E(I,J,K,L)=0.

   18 CONTINUE

    3 CONTINUE

      RETURN
      END