C     FILE HARVARD.F
C     ******************************************************************
C
C                                                            ***********
C                                                            * UBCON3D *
C                                                            ***********
      SUBROUTINE UBCON3D
C     ******************************************************************
C     *** USER BOTTOM BOUNDARY CONDITION SUBROUTINE
C     ******************************************************************
      INCLUDE 'FOR3D.CMN'
      ZI=(N+1)*DZ
      DO 250 J=1,NSOL+2
      ANGR=-FLDW/2.0*PI/180.0+((J-2)*PHI)
      BOTY(J)=BOTX(J)
C      BOTX(J)=...
250   CONTINUE
      RETURN
      END
C                                                            ***********
C                                                            * UBOTTOM *
C                                                            ***********
      SUBROUTINE UBOTTOM
C     ******************************************************************
C     ******************************************************************
C     * HARVARD BOTTOM                                                 *
C     * THIS SUBROUTINE IS DATA DEPENDENT AND NOT GENERAL PURPOSE      *
C     * HARVARD FORMAT IS SUBJECT TO CHANGE                            *   
C     ******************************************************************
C     ******************************************************************
      PARAMETER (HU=3,MXLEV=11,MXY=7991)
C     *** HU IS HARVARD INPUT UNIT NUMBER
C     *** MXLEV IS MAXIMUM NUMBER OF LEVELS IN INPUT DATA SET
C     *** MXY IS NUMBER OF DATA POINTS IN EACH LEVEL
      INCLUDE 'FOR3D.CMN'
      DIMENSION SVP(MXY,MXLEV),ZLEV(MXLEV),SLAT(MXSOL),SLNG(MXSOL)
      DIMENSION BOT(MXY),DZOLD(MXSOL),DZNEW(MXSOL)
      CHARACTER*80 SPDFIL,BOTFIL
      DATA RAD/6378400/
C
C       *** INPUT PARAMETERS FOR HARVARD DATA 
        NDAY=U1      ! DAY OF DATA SET : 3 MEANS 3RD DAY
        SLAT0=U2     ! LATITUDE OF STARTING FIELD
        SLNG0=U3     ! LONGITUDE OF STARTING FIELD
        DIR=U4       ! DIRECTION OF PROPAGATION OF CENTER RAY
        BOTRHO=U5    ! DENSITY IN BOTTOM
        BOTRHOG=U6   ! DENSITY GRADIENT IN BOTTOM
        BOTBETA=U7   ! ATTENUATION IN BOTTOM
        BOTBETAG=U8  ! ATTENUATION GRADIENT IN BOTTOM
        CWCB=U9      ! SOUND SPEED RATIO AT BOTTOM INTERFACE
        CGRAD=U10    ! SOUND SPEED GRADIENT IN BOTTOM
        SEDZ=U11     ! SEDIMENT THICKNESS
	write(6,*)u1,u2,u3,u4,u5,u6,u7,u8,u9,u10,u11,u12
C
        DO 102 J=1,NSOL
        SLAT(J)=U2
        SLNG(J)=U3
        ITRK(J)=0
        DZOLD(J)=0
102     CONTINUE
C
C       *** COMPUTE METERS PER DEGREE OF LATITUDE
        CONVLAT=PI*RAD/180.0
C
C       *** READ HARVARD DATA
C        OPEN(3)
	CALL ASSIGN(3,'E_LN3_F.SPD')
        READ(3,105)DLAT0,DLNG0,THTA,NX,NY,DX,DY,NLEV,
     C  ZLEV(1),(ZLEV(I),I=2,NLEV-1),ZLEV(NLEV)
105     FORMAT(2X,2F7.3,F8.5,2I6,/,2(F7.1),I6,/,(10F7.1))
C
        IF(NLEV.GT.MXLEV) THEN
          WRITE(NPU)'ERROR. TOO MANY LEVELS. INCREASE MXLEV. RECOMPILE.'
          NSVP=0
          RETURN
        ELSE
        ENDIF
C
C       *** GRID SIZE IN METERS
        DX=DX*1000.0
        DY=DY*1000.0
        CONVLNG2=RAD*COS(PI*DLAT0/180.)*PI/180.0
        EPS=ATAN2(DY,DX)
        DXSEPS=DX*SIN(EPS)
        NXY=NX*NY
        XBASIN=(NX-1)*DX  
        YBASIN=(NY-1)*DY
        ANGLE=THTA+ATAN2(YBASIN,XBASIN)
        DCENTER=SQRT((.5*XBASIN)**2+(.5*YBASIN)**2)
C       
        READ(3,107,END=192)NN,T,ITYP,K,ITMP1,ITMP2
107     FORMAT(2X,I5,F10.5,4I5)
C
        IF(NN.NE.NXY)GO TO 190
C
C       *** READ BOTTOM DEPTHS
        READ(3,108)(BOT(I),I=1,NXY)
108     FORMAT(1X,8F9.3)
C
        REWIND(3) ! END OF SOUND SPEED INPUT
C
      ZMAX=0.0
      CONVLNG0=RAD*COS(PI*SLAT0/180.)*PI/180.0
	DO 184 J=1,NSOL
	DZOLD(J)=0.0
184	CONTINUE
      DO 185 R=RA,RMAX,DR
      DO 180 J=1,NSOL
      STHETA=DIR*PI/180.0-(FLDW/2.0-(J-1)*DTH)*PI/180.0
      STHETA=90.0*PI/180.0-STHETA
      SLAT(J)=SLAT0+R*SIN(STHETA)/CONVLAT
      CONVLNG=.5*(CONVLNG0+RAD*COS(PI*SLAT(J)/180.)*PI/180.0)
      SLNG(J)=SLNG0-R*COS(STHETA)/CONVLNG
      A1=((DLNG0+DCENTER*COS(ANGLE)/CONVLNG2)-SLNG(J))*CONVLNG2
      B1=(SLAT(J)-(DLAT0-DCENTER*SIN(ANGLE)/CONVLAT))*CONVLAT
      C1=SQRT(A1**2+B1**2)
      PSI=ATAN2(B1,A1)-THTA
      RIS=(C1*COS(PSI)/(DX))+1
      IS=INT(RIS)
      IF(IS.GT.NX) THEN
        WRITE(NPU,194)
        RETURN
      ELSE
      ENDIF
C
      RJS=(C1*SIN(PSI)/(DY))+1
      JS=INT(RJS)
      IF(JS.GT.NY) THEN
        WRITE(NPU,195)
        RETURN
      ELSE
      ENDIF
C
      IF(IS.LT.1) THEN
        WRITE(NPU,194)
        RETURN
      ELSE
      ENDIF
C
      IF(JS.LT.1) THEN
        WRITE(NPU,195)
        RETURN
      ELSE
      ENDIF
C
      XIP=C1*COS(PSI)
      YJP=C1*SIN(PSI)
      A2=DX*(IS-1)
      B2=DY*(JS-1)
      C2=SQRT(A2**2+B2**2)
      PSI2=ATAN2(B2,A2)
      XIS=C2*COS(PSI2)
      YJS=C2*SIN(PSI2)
      XOFF=(XIP-XIS)
      YOFF=(YJP-YJS)
      ZETA=ATAN2(YOFF,XOFF)
      GAMMA=PI-EPS-ZETA
      H1=SQRT(XOFF*XOFF+YOFF*YOFF)
      H=DXSEPS/SIN(GAMMA)
      I1=(JS-1)*NX+IS
      I2=I1+1
      I3=I1+NX
      I4=I1+NX+1
      IF(H1.LE.H) THEN
        DEPTH=BOT(I1)+XOFF*(BOT(I2)-BOT(I1))/DX+YOFF*
     C  (BOT(I3)-BOT(I1))/DY
      ELSE
        DEPTH=BOT(I4)+(DX-XOFF)*(BOT(I3)-BOT(I4))/DX+(DY-YOFF)*
     C  (BOT(I2)-BOT(I4))/DY
      ENDIF
      IF(ZMAX.LT.DEPTH)ZMAX=DEPTH
      ITRK(J)=ITRK(J)+1
      TRACK(ITRK(J),1,J)=R
      TRACK(ITRK(J),2,J)=DEPTH
      IF(R.GT.RA)THEN
	DZNEW(J)=DEPTH
	IF(ABS(DZOLD(J)-DZNEW(J)).LT..5*ZA/N.AND.ITRK(J).GT.2)THEN
C     *** LIMIT NUMBER OF TRACK POINTS
            ITRK(J)=ITRK(J)-1
            TRACK(ITRK(J),1,J)=R
            TRACK(ITRK(J),2,J)=DEPTH
          ELSE
	DZOLD(J)=DZNEW(J)
          ENDIF
      ELSE
      ENDIF
180   CONTINUE
185   CONTINUE
C
C     *** END TRACK.
      DO 186 J=1,NSOL
      TRACK(ITRK(J)+1,1,J)=1.0E+38
      TRACK(ITRK(J)+1,2,J)=TRACK(ITRK(J),2,J)
C     *** INITIALIZE TRACK SEGMENT ARRAYS.
      R2(J)=TRACK(1,1,J)
      Z2(J)=TRACK(1,2,J)
186   CONTINUE
C
C       *** IF DEPTH OF WATER PLUS DEPTH OF SEDIMENT GT ZA, ADJUST N.
        IF(ZMAX+SEDZ.GT.ZA) THEN
          DZ=ZA/N
          WRITE(6,*)'ZMAX+SEDZ GT ZA. ZA= ',ZA,'. N= ',N,'. DZ= ',DZ,'.'
          WRITE(6,*)'ZMAX = ',ZMAX,'.  SEDZ = ',SEDZ,'.'
          RN=(ZMAX+SEDZ)/DZ
          N=INT(RN)
          ZA=N*DZ
          WRITE(6,*)'N RESET TO ',N,'.   ZA RESET TO ',ZA,'.'
        ELSE
        ENDIF
C
        IF(N.GT.MXNZ)THEN
          WRITE(NPU,*)'N TOO LARGE. INCREASE PARAMETER MXNZ. RECOMPILE.'
          STOP
        ELSE
        ENDIF
      RETURN
C
190   CONTINUE
      WRITE(NPU,191)
191   FORMAT(1X,'DATA MISMATCH. NX*NY DOES NOT EQUAL NN.')
      RETURN
C
192   WRITE(NPU,193)
193   FORMAT(1X,'READ ERROR. CHECK INPUT SOUND SPEED FILE.')
      NSVP=0
      RETURN
C
194   FORMAT(1X,'LONGITUDE OUTSIDE LIMITS OF SOUND SPEED DOMAIN.')
195   FORMAT(1X,'LATITUDE OUTSIDE LIMITS OF SOUND SPEED DOMAIN.')
196   FORMAT(1X,'LAT = ',F6.1,'  LON  = ',F6.2)
C
200   CONTINUE
300   CONTINUE
400   CONTINUE
      RETURN
      END
C                                                            ***********
C                                                            * UEXACT  *
C                                                            ***********
      SUBROUTINE UEXACT
C     ******************************************************************
C     * EXACT SOLUTION  - USED FOR TEST PURPOSES
C     * USER COMPUTES EXACT SOLUTION UB AT PRESENT RANGE RA, DEPTH ZI,
C     * AND AZIMUTH ANG. MAIN PROGRAM COMPUTES UA. RELATIVE ERROR 
C     * MAY BE COMPUTED AS SHOWN BELOW. RA, ZI, ANG, UA, UB IN COMMON.
C     ******************************************************************
      INCLUDE 'FOR3D.CMN'
C      UB=... USER COMPUTES THIS
      PL=-20.0*ALOG10(CABS(UB))+10.0*ALOG10(RA)
      UC=(UA-UB)
      WRITE(NPU,120) PL,UB
120   FORMAT(20X,3X,F10.3,3X,'(',E12.5,2X,E12.5,' )')
      RELERR=CABS(UC/UB)
      WRITE(NPU,150) UC,RELERR
150   FORMAT(36X,'(',E12.5,2X,E12.5,' )',2X,E12.5)
      RETURN
      END
C                                                            ***********
C                                                            * UPORT3D *
C                                                            ***********
      SUBROUTINE UPORT3D
C     ******************************************************************
C     * PORT SIDEWALL BOUNDARY CONDITION
C     ******************************************************************
      INCLUDE 'FOR3D.CMN'
      J=0
      ANGR=-FLDW/2.0*PI/180.0+((J-1)*PHI)
      DO 100 I=1,N                                             
      ZI=I*DZ
      PORTY(I)=PORTX(I)
C      PORTX(I)=...
100   CONTINUE
      RETURN                                                     
      END
C                                                            ***********
C                                                            * USCON3D *
C                                                            ***********
      SUBROUTINE USCON3D
C     *****************************************************************
C     *** USER SURFACE BOUNDARY CONDITION SUBROUTINE
C     *****************************************************************
      INCLUDE 'FOR3D.CMN'
      ZI=0.0
      DO 250 J=1,NSOL+2
      ANGR=-FLDW/2.0*PI/180.0+((J-2)*PHI)
      SURY(J)=SURX(J)
C      SURX(J)=....
250   CONTINUE
      RETURN
      END
C                                                            ***********
C                                                            * USFLD3D *
C                                                            ***********
      SUBROUTINE USFLD3D
C     ******************************************************************
C     *** USER STARTING FIELD 
C     *** USER WRITES THIS SUBROUTINE IF GAUSSIAN FIELD NOT DESIRED
C     *** USFLD3D IS CALLED IF ISF = 1
C     ******************************************************************
C     *** UFIELD SUBROUTINE SUPPLIES:
C     *** U   - COMPLEX STARTING FIELD
C     ******************************************************************
      INCLUDE 'FOR3D.CMN'
      DO 200 J=1,NSOL
      M=(J-1)*N
      ANGR=-FLDW/2.0*PI/180.0+((J-1)*PHI)
      DO 100 I=1,N                                             
      ZI=I*DZ
C      U(M+I)= ....
100   CONTINUE
200   CONTINUE
      RETURN                                                     
      END
C                                                            ***********
C                                                            * USTBD3D *
C                                                            ***********
      SUBROUTINE USTBD3D
C     ******************************************************************
C     * STARBOARD SIDEWALL BOUNDARY CONDITION
C     ******************************************************************
      INCLUDE 'FOR3D.CMN'
      J=NSOL+1
      ANGR=-FLDW/2.0*PI/180.0+((J-1)*PHI)
      DO 100 I=1,N                                             
      ZI=I*DZ
      STBDY(I)=STBDX(I)
C      STBDX(I)=...
100   CONTINUE
      RETURN                                                     
      END
C                                                            ***********
C                                                            * USVP3D  *
C                                                            ***********
      SUBROUTINE USVP3D
C     ******************************************************************
C     ******************************************************************
C     * HARVARD SOUND SPEED PROFILES                                   *
C     * THIS SUBROUTINE IS DATA DEPENDENT AND NOT GENERAL PURPOSE      *
C     * HARVARD FORMAT IS SUBJECT TO CHANGE                            *   
C     ******************************************************************
C     ******************************************************************
      PARAMETER (HU=3,MXLEV=11,MXY=7991)
C     *** HU IS HARVARD INPUT UNIT NUMBER
C     *** MXLEV IS MAXIMUM NUMBER OF LEVELS IN INPUT DATA SET
C     *** MXY IS NUMBER OF DATA POINTS IN EACH LEVEL
      INCLUDE 'FOR3D.CMN'
      DIMENSION SVP(MXY,MXLEV),ZLEV(MXLEV),SLAT(MXSOL),SLNG(MXSOL)
      DIMENSION BOT(MXY)
      CHARACTER*80 SPDFIL,BOTFIL
      DATA IPASS/0/
      DATA RAD/6378400/
      GO TO (100,200,300,400) ,KSVP
      NSVP=0
      RETURN
100   CONTINUE
C
C     *** IF THIS IS FIRST PASS, GET PARAMETERS AND SOUND SPEED DATA
      IF(IPASS.EQ.0) THEN
C       *** INPUT PARAMETERS FOR HARVARD DATA 
        NDAY=U1      ! DAY OF DATA SET : 3 MEANS 3RD DAY
        SLAT0=U2     ! LATITUDE OF STARTING FIELD
        SLNG0=U3     ! LONGITUDE OF STARTING FIELD
        DIR=U4       ! DIRECTION OF PROPAGATION OF CENTER RAY
        BOTRHO=U5    ! DENSITY IN BOTTOM
        BOTRHOG=U6   ! DENSITY GRADIENT IN BOTTOM
        BOTBETA=U7   ! ATTENUATION IN BOTTOM
        BOTBETAG=U8  ! ATTENUATION GRADIENT IN BOTTOM
        CWCB=U9      ! SOUND SPEED RATIO AT BOTTOM INTERFACE
        CGRAD=U10    ! SOUND SPEED GRADIENT IN BOTTOM
        SEDZ=U11     ! SEDIMENT THICKNESS
C
        DO 102 J=1,NSOL
        SLAT(J)=U2
        SLNG(J)=U3
102     CONTINUE
C
C       *** COMPUTE METERS PER DEGREE OF LATITUDE
        CONVLAT=PI*RAD/180.0
C
C       *** READ HARVARD DATA
C        OPEN(3)
        READ(3,105)DLAT0,DLNG0,THTA,NX,NY,DX,DY,NLEV,
     C  ZLEV(1),(ZLEV(I),I=2,NLEV-1),ZLEV(NLEV)
105     FORMAT(2X,2F7.3,F8.5,2I6,/,2(F7.1),I6,/,(10F7.1))
C
        IF(NLEV.GT.MXLEV) THEN
          WRITE(NPU)'ERROR. TOO MANY LEVELS. INCREASE MXLEV. RECOMPILE.'
          NSVP=0
          RETURN
        ELSE
        ENDIF
C
C       *** GRID SIZE IN METERS
        DX=DX*1000.0
        DY=DY*1000.0
        CONVLNG2=RAD*COS(PI*DLAT0/180.)*PI/180.0
        EPS=ATAN2(DY,DX)
        DXSEPS=DX*SIN(EPS)
        NXY=NX*NY
        XBASIN=(NX-1)*DX  
        YBASIN=(NY-1)*DY
        ANGLE=THTA+ATAN2(YBASIN,XBASIN)
        DCENTER=SQRT((.5*XBASIN)**2+(.5*YBASIN)**2)
C       
        READ(3,107,END=192)NN,T,ITYP,K,ITMP1,ITMP2
107     FORMAT(2X,I5,F10.5,4I5)
C
        IF(NN.NE.NXY)GO TO 190
C
C       *** READ BOTTOM DEPTHS
        READ(3,108)(BOT(I),I=1,NXY)
108     FORMAT(1X,8F9.3)
C
C       *** READ NLEV LEVELS OF SPEED DATA FOR DAY NDAY
        DO 110 IDAY=1,NDAY
          DO 110 L=1,NLEV
          READ(3,107,END=192)NN,T,ITYP,K,ITMP1,ITMP2
          IF(NN.NE.NXY)GO TO 190
          READ(3,108)(SVP(I,L),I=1,NXY)
110     CONTINUE
C
C        CALL CLOSE(3,STATUS='KEEP') ! END OF SOUND SPEED INPUT
C
        RAOLD=RA
        CONVLNG0=RAD*COS(PI*SLAT0/180.)*PI/180.0
        IPASS=1
      ELSE
      ENDIF
C
C     *** COME HERE IMMEDIATELY IF IPASS NE 0
      DO 180 J=1,NSOL
      STHETA=DIR*PI/180.0-(FLDW/2.0-(J-1)*DTH)*PI/180.0
      STHETA=90.0*PI/180.0-STHETA
      SLAT(J)=SLAT0+RA*SIN(STHETA)/CONVLAT
      CONVLNG=.5*(CONVLNG0+RAD*COS(PI*SLAT(J)/180.)*PI/180.0)
      SLNG(J)=SLNG0-RA*COS(STHETA)/CONVLNG
      A1=((DLNG0+DCENTER*COS(ANGLE)/CONVLNG2)-SLNG(J))*CONVLNG2
      B1=(SLAT(J)-(DLAT0-DCENTER*SIN(ANGLE)/CONVLAT))*CONVLAT
      C1=SQRT(A1**2+B1**2)
      PSI=ATAN2(B1,A1)-THTA
      RIS=(C1*COS(PSI)/(DX))+1
      IS=INT(RIS)
      IF(IS.GT.NX) THEN
        WRITE(NPU,194)
        KSVP=0
        RETURN
      ELSE
      ENDIF
C
      RJS=(C1*SIN(PSI)/(DY))+1
      JS=INT(RJS)
      IF(JS.GT.NY) THEN
        WRITE(NPU,195)
        KSVP=0
        RETURN
      ELSE
      ENDIF
C
      IF(IS.LT.1) THEN
        WRITE(NPU,194)
        KSVP=0
        RETURN
      ELSE
      ENDIF
C
      IF(JS.LT.1) THEN
        WRITE(NPU,195)
        KSVP=0
        RETURN
      ELSE
      ENDIF
C
      XIP=C1*COS(PSI)
      YJP=C1*SIN(PSI)
      A2=DX*(IS-1)
      B2=DY*(JS-1)
      C2=SQRT(A2**2+B2**2)
      PSI2=ATAN2(B2,A2)
      XIS=C2*COS(PSI2)
      YJS=C2*SIN(PSI2)
      XOFF=(XIP-XIS)
      YOFF=(YJP-YJS)
      ZETA=ATAN2(YOFF,XOFF)
      GAMMA=PI-EPS-ZETA
      H1=SQRT(XOFF*XOFF+YOFF*YOFF)
      H=DXSEPS/SIN(GAMMA)
      I1=(JS-1)*NX+IS
      I2=I1+1
      I3=I1+NX
      I4=I1+NX+1
      IF(H1.LE.H) THEN
        DEPTH=BOT(I1)+XOFF*(BOT(I2)-BOT(I1))/DX+YOFF*
     C  (BOT(I3)-BOT(I1))/DY
      ELSE
        DEPTH=BOT(I4)+(DX-XOFF)*(BOT(I3)-BOT(I4))/DX+(DY-YOFF)*
     C  (BOT(I2)-BOT(I4))/DY
      ENDIF
C
      IF(J.EQ.1) DEPTH1=DEPTH
      XOFF=XOFF/DX
      YOFF=YOFF/DY
C
C     *** INTERPOLATE
      DO 120 LEV=1,NLEV
      P1=SVP(I1,LEV)
      P2=SVP(I2,LEV)
      P3=SVP(I3,LEV)
      P4=SVP(I4,LEV)
      A=(YOFF*P3+(1.0-YOFF)*P1)
      B=(YOFF*P4+(1.0-YOFF)*P2)
      CSVP(LEV,J)=(A*(1.0-XOFF)+B*XOFF)
      ZSVP(LEV,J)=ZLEV(LEV)
120   CONTINUE
C
      ILYR=1
      RHO(ILYR,J)=1.0
      RHOG(ILYR,J)=0.0
      BETA(ILYR,J)=0.0
      BETAG(ILYR,J)=0.0
      IF(DEPTH.GT.ZSVP(NLEV,J)) GO TO 140
      IF(DEPTH.EQ.ZSVP(NLEV,J)) GO TO 141
C
      DO 130 L=1,NLEV
      LEV=L
      IF(DEPTH.EQ.ZSVP(LEV,J)) GO TO 138
      IF(DEPTH.LT.ZSVP(LEV,J)) GO TO 135
130   CONTINUE
      WRITE(NPU,*)'ERROR IN LOGIC. ???'
135   CONTINUE
      CSVP(LEV,J)=CSVP(LEV-1,J)+(CSVP(LEV,J)-CSVP(LEV-1,J))*
     C(DEPTH-ZSVP(LEV-1,J))/(ZSVP(LEV,J)-ZSVP(LEV-1,J))
      ZSVP(LEV,J)=DEPTH
138   IXSVP(ILYR,J)=LEV
      ZLYR(ILYR,J)=DEPTH
      NSVP=LEV
      GO TO 145
C
140   ZSVP(NLEV+1,J)=DEPTH
      CSVP(NLEV+1,J)=CSVP(NLEV-1,J)+(CSVP(NLEV,J)-CSVP(NLEV-1,J))*
     C(ZSVP(NLEV+1,J)-ZSVP(NLEV-1,J))/(ZSVP(NLEV,J)-ZSVP(NLEV-1,J))
      IXSVP(ILYR,J)=NLEV+1
      ZLYR(ILYR,J)=DEPTH
      NSVP=NLEV+1
      GO TO 145
C
141   IXSVP(ILYR,J)=NLEV
      ZLYR(ILYR,J)=ZSVP(NLEV,J)
      NSVP=NLEV
C
145   CONTINUE
      ILYR=2
      RHO(ILYR,J)=BOTRHO
      RHOG(ILYR,J)=BOTRHOG
c	write(6,*)'usvp botrho,g ',botrho,botrhog
      ZG=ZLYR(ILYR,J)
      BETA(ILYR,J)=BOTBETA
      BETAG(ILYR,J)=BOTBETAG
      ZSVP(NSVP+1,J)=DEPTH
      CSVP(NSVP+1,J)=CSVP(NSVP,J)/CWCB
      IF(ZA-DEPTH.GE.SEDZ)THEN 
        ZZ=SEDZ
      ELSE
        ZZ=ZA-DEPTH
      ENDIF
C
      CSVP(NSVP+2,J)=CSVP(NSVP+1,J)+(ZZ)*CGRAD
      ZSVP(NSVP+2,J)=DEPTH+ZZ
      ZLYR(ILYR,J)=DEPTH+ZZ
      IXSVP(ILYR,J)=NSVP+2
      NLYRS(J)=ILYR
180   CONTINUE
      NSVP=NSVP+2
      RAOLD=RA
      RETURN
C
190   CONTINUE
      WRITE(NPU,191)
191   FORMAT(1X,'DATA MISMATCH. NX*NY DOES NOT EQUAL NN.')
      NSVP=0
      RETURN
C
192   WRITE(NPU,193)
193   FORMAT(1X,'READ ERROR. CHECK INPUT SOUND SPEED FILE.')
      NSVP=0
      RETURN
C
194   FORMAT(1X,'LONGITUDE OUTSIDE LIMITS OF SOUND SPEED DOMAIN.')
195   FORMAT(1X,'LATITUDE OUTSIDE LIMITS OF SOUND SPEED DOMAIN.')
196   FORMAT(1X,'LAT = ',F6.1,'  LON  = ',F6.2)
C
200   CONTINUE
300   CONTINUE
400   CONTINUE
      NSVP=0
      RETURN
      END