
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      EXTERNAL FNA
      COMMON S(1001),P(1001),G(1001),T(1001),X(1001),VS(1001)
	COMMON DZ, SURF, NV
      DIMENSION RH(500)
	CHARACTER*15 FINP, FOT
C
C            MAIN PROGRAM TO CALL THE OTHER PIECES
C
C                REVISED APRIL,1985 (DCC)
C
c		revised to be DOS friendly 6/28/97 (KHB)
c		revised to update style 11/28/2005 (KHB)

c		Test with file "testcros.in"


	INPF = 10
	IOPF = 11


c		Set the horizontal grid resolution NN
c		To change the grid resolution, change NN here.
c			If you want NN > 500, change the dimension here and
c			elsewhere

	NN = 25

	WRITE(6,970)
	READ(5,*) FINP
	WRITE(6,971)
	READ(5,*) FOT

	OPEN (INPF,FILE=FINP)
	OPEN (IOPF,FILE=FOT)

c		nv = # of vertical grid points
c 		dd1 = 0 for rigid lid
c		dd1 = 1 for free surface

      READ(INPF,*)NV,DD1
      READ(INPF,*)F,XMAX

c 		F that you read is the true f/1e-5
c		xmax = max offshore distance in km

      XMAX=XMAX*1.E5
      F=F*1.E-5
      IF(DD1.EQ.0.)WRITE(IOPF,908)
      IF(DD1.EQ.1.)WRITE(IOPF,909)
      WRITE(IOPF,907)F,XMAX

     
      DX=XMAX/(NN-1)
      CALL DEP(DX,RH,NN)
c		DEP reads the depth profile as in BIGLOAD3
     
      DZ=RH(NN)/(NV-1)
      K1=1
      X(1)=0.0
      X1=0.0
      X2=0.0
      H1=0.0
      H2=RH(1)

      DO 100 I=2,NV
	      Z=FLOAT(I-1)*DZ
	      IF(Z.GE.RH(NN))Z=RH(NN)

160 	      IF(Z.LE.H2)GO TO 150
				K1=K1+1
				H1=RH(K1-1)
				H2=RH(K1)
				X1=DX*FLOAT(K1-2)
				X2=DX*FLOAT(K1-1)
		  GO TO 160

150 				DH=H2-H1
				SX=(X2-X1)/DH
				BX=-(X2*H1-X1*H2)/DH
				X(I)=SX*Z+BX
100   CONTINUE

      CALL NSQ(RH(NN),F,DZ,VS,NV)
c		NSQ reads stratification as in BIGLOAD3

      SURF=DD1*VS(1)*F*F/981.

      READ(INPF,*)NRS
c		NRS = number of modes to seek
      WRITE(IOPF,901)DZ
      WRITE(IOPF,902)(X(I),I=1,NV)
      WRITE(IOPF,903)
      WRITE(IOPF,902)(VS(I),I=1,NV)

      DO 200 IR=1,NRS
		READ(INPF,*)XMIN,XMAX,XINC
c			XMIN,XMAX, XINC = alongshore wavenumbers and increment*1e7
		XMIN=XMIN*1.E-7
		XMAX=XMAX*1.E-7
		XINC=XINC*1.E-7
		CALL ROOT(XMIN,XMAX,XINC,FNA,X3,IROOT)
		
		IF(IROOT.EQ.0)GO TO 200
		DO I=2,NV
     			P(I)=P(I)*DEXP(-X3*(X(I)+X(I-1))*.5)
		END DO	
		
		PSUM=0.0
		DO I=2,NV
  			PSUM=PSUM+P(I)**2
		END DO
		
		PSUM=DSQRT(DZ*PSUM)
		DO I=2,NV
  			P(I)=P(I)/PSUM
		END DO
		
		CF=F/X3
		WRITE(IOPF,904)X3,CF
		WRITE(IOPF,905)
		WRITE(IOPF,902)(P(I),I=2,NV)
200   CONTINUE
	
	CLOSE(IOPF)
	CLOSE(INPF)


901	FORMAT(//'  INVERSE TOPOGRAPHY, DZ=',D15.5/)
902   FORMAT(1X,10D12.5)
903   FORMAT(//'  NSQUARED/FSQUARED'/)
904   FORMAT(////'  DISPERSION CURVE CROSSES W=F AT L=',D15.5,5X,       &
     &  'PHASE SPEED=',D15.5/)
905   FORMAT(/'  PRESSURE FIELD (SURFACE TO BOTTOM)'/)
907   FORMAT('  F,XMAX=',2D15.5)
908   FORMAT(/'  RIGID LID'/)
909   FORMAT(//'  FREE SURFACE'/)
970	FORMAT(' input file?')
971	FORMAT(' output file?')
      
	STOP
      END



      SUBROUTINE DEP(DX,RH,NN)

      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION XL(1001),BX(1001),RH(500)
C
C          SUBROUTINE TO READ AND INTERPOLATE DEPTH PROFILE
C                  RH= DEPTH
C
	INPF = 10
	IOPF = 11

      READ(INPF,*)NRX
      DO I=1,NRX
		READ(INPF,*)XL(I),BX(I)
		BX(I)=BX(I)*100.
    		XL(I)=XL(I)*1.E5
      END DO

      IF(NRX.NE.1) THEN
		DO I=2,NRX
    			IF(BX(I).LE.BX(I-1))GO TO 100
		END DO
	END IF

      RHMA=BX(NRX)
      XMA=XL(NRX)
      RH(1)=BX(1)

      DO 20 N=2,25
		X=DX*FLOAT(N-1)
		IF(X.GT.XMA)GO TO 15

		IC=0
		DO 8 J=2,NRX
     			IF(IC.NE.0)GO TO 8
				I=J
			IF(X.GT.XL(I))GO TO 8
     				IC=I
8		CONTINUE

		IM=I-1
		RHX=(BX(I)-BX(IM))/(XL(I)-XL(IM))
		XX=X-XL(IM)
		RH(N)=BX(IM)+RHX*XX
		GO TO 20

15		RH(N)=RHMA
20    CONTINUE
      GO TO 200

100   WRITE(IOPF,901)

901   FORMAT(/'  TOPOGRAPHY IS NOT MONOTONIC'//)

200   RETURN
      END



      SUBROUTINE NSQ(RHMA,F,DZ,VS,NV)

      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION XL(10000),VS(1001)
C
C           SUBROUTINE TO READ AND INTERPOLATE BUOYANCY FREQUENCY
C
	INPF = 10
	IOPF = 11

      FSQ=F*F
      READ(INPF,*)NR,DZR,ALPH
      ALPH=ALPH*1.E5
      DZR=DZR*100.

      READ(INPF,*)CMLT
      DO I=1,NR
         READ(INPF,*)XL(I)
   	   XL(I)=XL(I)*CMLT
      END DO

      TF=XL(NR)
      ZZZ=DZR*FLOAT(NR-1)
      I=NR+1
50    Z=DZR*FLOAT(I-1)
      XL(I)=TF*DEXP((ZZZ-Z)/ALPH)

      IF(Z.GT.RHMA)GO TO 100
		I=I+1
      GO TO 50

100   K1=1
      VS(1)=XL(1)/FSQ
      H1=0.0
      H2=DZR

      DO 200 I=2,NV
		Z=FLOAT(I-1)*DZ
		IF(Z.GE.RHMA)Z=RHMA
160		IF(Z.LE.H2)GO TO 150
			K1=K1+1
			H1=FLOAT(K1-1)*DZR
			H2=FLOAT(K1)*DZR
		GO TO 160

150		DH=H2-H1
		SN=(XL(K1+1)-XL(K1))/DH
		BN=-(XL(K1+1)*H1-XL(K1)*H2)/DH
		VS(I)=(SN*Z+BN)/FSQ
200   CONTINUE

      RETURN
      END



      REAL*8 FUNCTION FNA(W)

      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      COMMON S(1001),P(1001),G(1001),T(1001),X(1001),VS(1001)
	COMMON DZ, SURF, NV
C
C          SUBFUNCTION TO INTEGRATE FROM Z=0 TO Z=-H
C
      NVM=NV-2
      W2=W*W
      TW=-2.*W

      DO I=1,NV
         EX=DEXP(TW*X(I))
         T(I)=W2*EX
	   S(I)=EX/VS(I)
      END DO

      P(NV)=1.0
      G(NV)=0.0
      
	DO I=1,NVM
	      J=NV-I
	      K=J+1
	      G(J)=G(K)-.5*DZ*(T(J)+T(K))*P(K)
	      P(J)=P(K)+DZ*G(J)/S(J)
      END DO

      G(1)=G(2)-.5*DZ*(T(1)+T(2))*P(2)
      FNA=G(1)+S(1)*SURF*.5*(3.*P(2)-P(3))

      RETURN
      END



      SUBROUTINE ROOT(X5,X6,X7,FNA,X3,IROOT)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
	EXTERNAL FNA
C
C          SUBROUTINE TO FIND ROOT OF FNA BETWEEN X5 AND X6
C
	INPF = 10
	IOPF = 11

      IROOT=1
220   IQ=INT((X6-X5)/X7)-1
      I1=0
240   X1=X5+I1*X7
      X2=X1+X7
      Y1=FNA(X1)
      Y2=FNA(X2)
      IF(Y1*Y2.LT.0.0)GO TO 340
		I1=I1+1
      IF(I1.NE.IQ)GO TO 240
		WRITE(IOPF,913)
		IROOT=0
      GO TO 400

340   IF(X1.NE.X5)GO TO 370
		X7=X7/2.
      GO TO 220

370   X2=X1
      X1=X1-0.1*X7
      Y1=FNA(X1)
      Y2=FNA(X2)
420   X3=(Y2*X1-Y1*X2)/(Y2-Y1)

      IF(X3.GT.X5.AND.X3.LT.X6)GO TO 200
		WRITE(IOPF,911)X3
200   Y3=FNA(X3)
     
      IF(DABS((X3-X2)/X2).LT.1.E-7)GO TO 400
		X1=X2
		X2=X3
		Y1=Y2
		Y2=Y3
      GO TO 420

911    FORMAT(/'  PROJECTION OUT OF INTERVAL:',D15.5//)
913    FORMAT(//'  NO ROOT FOUND'/)

400   RETURN
      END 