!     FILE HARVARD.F
!     ******************************************************************
!
!                                                            ***********
!                                                            * UBCON3D *
!                                                            ***********
      SUBROUTINE UBCON3D
!     ******************************************************************
!     *** USER BOTTOM BOUNDARY CONDITION SUBROUTINE
!     ******************************************************************
!ljh  INCLUDE 'FOR3D.CMN'
      USE commons
      ZI=(N+1)*DZ
      DO 250 J=1,NSOL+2
        ANGR=-FLDW/2.0*PI/180.0+((J-2)*PHI)
        BOTY(J)=BOTX(J)
!       BOTX(J)=...
  250 CONTINUE
      RETURN
      END
!                                                            ***********
!                                                            * UBOTTOM *
!                                                            ***********
      SUBROUTINE UBOTTOM
!     ******************************************************************
!     ******************************************************************
!     * HARVARD BOTTOM                                                 *
!     * THIS SUBROUTINE IS DATA-DEPENDENT AND NOT GENERAL PURPOSE      *
!     * HARVARD FORMAT IS SUBJECT TO CHANGE                            *   
!     ******************************************************************
!     ******************************************************************
      USE commons
      PARAMETER (HU=3,MXLEV=11,MXY=7991)
!     *** HU IS HARVARD INPUT UNIT NUMBER
!     *** MXLEV IS MAXIMUM NUMBER OF LEVELS IN INPUT DATA SET
!     *** MXY IS NUMBER OF DATA POINTS IN EACH LEVEL
!ljh  INCLUDE 'FOR3D.CMN'
!ljh      DIMENSION SVP(MXY,MXLEV),ZLEV(MXLEV),SLAT(MXSOL),SLNG(MXSOL)
      DIMENSION ZLEV(MXLEV),SLAT(MXSOL),SLNG(MXSOL)
      DIMENSION BOT(MXY),DZOLD(MXSOL),DZNEW(MXSOL)
!ljh      CHARACTER*80 SPDFIL,BOTFIL
      DATA RAD/6378400/
!ljh
        SAVE
!ljh
!
!     *** INPUT PARAMETERS FOR HARVARD DATA 
!ljh      NDAY=U1      ! DAY OF DATA SET : 3 MEANS 3RD DAY
      NDAY=int(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
!
      DO 102 J=1,NSOL
        SLAT(J)=U2
        SLNG(J)=U3
        ITRK(J)=0
        DZOLD(J)=0
  102 CONTINUE
!
!     *** COMPUTE METERS PER DEGREE OF LATITUDE
      CONVLAT=PI*RAD/180.0
!
!     *** READ HARVARD DATA
!      OPEN(3)
!ljh   CALL ASSIGN(3,'E_LN3_F.SPD')
!ljh
      OPEN(UNIT=3,FILE='e_ln3_f.spd')
      READ(3,105)DLAT0,DLNG0,THTA,NX,NY,DX,DY,NLEV,    &
      ZLEV(1),(ZLEV(I),I=2,NLEV-1),ZLEV(NLEV)
  105 FORMAT(2X,2F7.3,F8.5,2I6,/,2(F7.1),I6,/,(10F7.1))
!
      IF(NLEV.GT.MXLEV) THEN
        WRITE(NPU)'ERROR. TOO MANY LEVELS. INCREASE MXLEV. RECOMPILE.'
        NSVP=0
        RETURN
      ELSE
      ENDIF
!
!     *** 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)
!     
      READ(3,107,END=192)NN,T,ITYP,K,ITMP1,ITMP2
  107 FORMAT(2X,I5,F10.5,4I5)
!
      IF(NN.NE.NXY)GO TO 190
!
!     *** READ BOTTOM DEPTHS
      READ(3,108)(BOT(I),I=1,NXY)
  108 FORMAT(1X,8F9.3)
!
      REWIND(3) ! END OF SOUND SPEED INPUT
!
      ZMAX=0.0
      CONVLNG0=RAD*COS(PI*SLAT0/180.)*PI/180.0
      DO 184 J=1,NSOL
        DZOLD(J)=0.0
  184 CONTINUE
!ljh  Can no longer have real number as do-loop counter.
!ljh      DO 185 R=RA,RMAX,DR
      R = RA - DR
      DO WHILE ( R .LT. RMAX )                                       ! Step out in range
        R = R + 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                                ! stheta is radial orienta wrt center DIR
          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
!
          RJS=(C1*SIN(PSI)/(DY))+1
          JS=INT(RJS)
          IF(JS.GT.NY) THEN
            WRITE(NPU,195)
            RETURN
          ELSE
          ENDIF
!
          IF(IS.LT.1) THEN
            WRITE(NPU,194)
            RETURN
          ELSE
          ENDIF
!
          IF(JS.LT.1) THEN
            WRITE(NPU,195)
            RETURN
          ELSE
          ENDIF
!
          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*    &
            (BOT(I3)-BOT(I1))/DY
          ELSE
            DEPTH=BOT(I4)+(DX-XOFF)*(BOT(I3)-BOT(I4))/DX+(DY-YOFF)*    &
            (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
!             *** 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
!ljh  185 CONTINUE
      END DO
!
!     *** 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)
!       *** INITIALIZE TRACK SEGMENT ARRAYS.
        R2(J)=TRACK(1,1,J)
        Z2(J)=TRACK(1,2,J)
  186 CONTINUE
!
!       *** 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
!
        IF(N.GT.MXNZ)THEN
          WRITE(NPU,*)'N TOO LARGE. INCREASE PARAMETER MXNZ. RECOMPILE.'
          STOP
        ELSE
        ENDIF
      RETURN
!
  190   CONTINUE
      WRITE(NPU,191)
  191   FORMAT(1X,'DATA MISMATCH. NX*NY DOES NOT EQUAL NN.')
      RETURN
!
  192   WRITE(NPU,193)
  193   FORMAT(1X,'READ ERROR. CHECK INPUT SOUND SPEED FILE.')
      NSVP=0
      RETURN
!
!  194   FORMAT(1X,'UBOTTOM: LONGITUDE OUTSIDE LIMITS OF SOUND SPEED    &
!       DOMAIN.')
!  195   FORMAT(1X,'UBOTTOM: LATITUDE OUTSIDE LIMITS OF SOUND SPEED    &
!       DOMAIN.')
  194   FORMAT(1X,'UBOTTOM: LONGITUDE OUTSIDE LIMITS OF SOUND SPEED DOMAIN.')
  195   FORMAT(1X,'UBOTTOM: LATITUDE OUTSIDE LIMITS OF SOUND SPEED DOMAIN.')
!ljh  196   FORMAT(1X,'LAT = ',F6.1,'  LON  = ',F6.2)
!
  200   CONTINUE
  300   CONTINUE
  400   CONTINUE
      RETURN
      END
!                                                            ***********
!                                                            * UEXACT  *
!                                                            ***********
      SUBROUTINE UEXACT
!     ******************************************************************
!     * EXACT SOLUTION  - USED FOR TEST PURPOSES
!     * USER COMPUTES EXACT SOLUTION UB AT PRESENT RANGE RA, DEPTH ZI,
!     * AND AZIMUTH ANG. MAIN PROGRAM COMPUTES UA. RELATIVE ERROR 
!     * MAY BE COMPUTED AS SHOWN BELOW. RA, ZI, ANG, UA, UB IN COMMON.
!     ******************************************************************
!ljh  INCLUDE 'FOR3D.CMN'
      USE commons
!      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
!                                                            ***********
!                                                            * UPORT3D *
!                                                            ***********
      SUBROUTINE UPORT3D
!     ******************************************************************
!     * PORT SIDEWALL BOUNDARY CONDITION
!     ******************************************************************
!ljh  INCLUDE 'FOR3D.CMN'
      USE commons
      J=0
      ANGR=-FLDW/2.0*PI/180.0+((J-1)*PHI)
      DO 100 I=1,N                                             
      ZI=I*DZ
      PORTY(I)=PORTX(I)
!      PORTX(I)=...
  100   CONTINUE
      RETURN                                                     
      END
!                                                            ***********
!                                                            * USCON3D *
!                                                            ***********
      SUBROUTINE USCON3D
!     *****************************************************************
!     *** USER SURFACE BOUNDARY CONDITION SUBROUTINE
!     *****************************************************************
!ljh  INCLUDE 'FOR3D.CMN'
      USE commons
      ZI=0.0
      DO 250 J=1,NSOL+2
      ANGR=-FLDW/2.0*PI/180.0+((J-2)*PHI)
      SURY(J)=SURX(J)
!      SURX(J)=....
  250   CONTINUE
      RETURN
      END
!                                                            ***********
!                                                            * USFLD3D *
!                                                            ***********
      SUBROUTINE USFLD3D
!     ******************************************************************
!     *** USER STARTING FIELD 
!     *** USER WRITES THIS SUBROUTINE IF GAUSSIAN FIELD NOT DESIRED
!     *** USFLD3D IS CALLED IF ISF = 1
!     ******************************************************************
!     *** UFIELD SUBROUTINE SUPPLIES:
!     *** U   - COMPLEX STARTING FIELD
!     ******************************************************************
!ljh  INCLUDE 'FOR3D.CMN'
      USE commons
      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
!      U(M+I)= ....
  100   CONTINUE
  200   CONTINUE
      RETURN                                                     
      END
!                                                            ***********
!                                                            * USTBD3D *
!                                                            ***********
      SUBROUTINE USTBD3D
!     ******************************************************************
!     * STARBOARD SIDEWALL BOUNDARY CONDITION
!     ******************************************************************
!ljh  INCLUDE 'FOR3D.CMN'
      USE commons
      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)
!      STBDX(I)=...
  100   CONTINUE
      RETURN                                                     
      END
!                                                            ***********
!                                                            * USVP3D  *
!                                                            ***********
      SUBROUTINE USVP3D
!     ******************************************************************
!     ******************************************************************
!     * HARVARD SOUND SPEED PROFILES                                   *
!     * THIS SUBROUTINE IS DATA-DEPENDENT AND NOT GENERAL PURPOSE      *
!     * HARVARD FORMAT IS SUBJECT TO CHANGE                            *   
!     ******************************************************************
!     ******************************************************************
      USE commons
      PARAMETER (HU=3,MXLEV=11,MXY=7991)
!     *** HU IS HARVARD INPUT UNIT NUMBER
!     *** MXLEV IS MAXIMUM NUMBER OF LEVELS IN INPUT DATA SET
!     *** MXY IS NUMBER OF DATA POINTS IN EACH LEVEL
!ljh  INCLUDE 'FOR3D.CMN'
      DIMENSION SVP(MXY,MXLEV),ZLEV(MXLEV),SLAT(MXSOL),SLNG(MXSOL)
      DIMENSION BOT(MXY)
!ljh      CHARACTER*80 SPDFIL,BOTFIL
      DATA IPASS/0/
      DATA RAD/6378400/
!ljh
        SAVE
!ljh
      GO TO (100,200,300,400) ,KSVP
      NSVP=0
      RETURN
  100 CONTINUE
!
!     *** IF THIS IS FIRST PASS, GET PARAMETERS AND SOUND SPEED DATA
      IF(IPASS.EQ.0) THEN
!       *** INPUT PARAMETERS FOR HARVARD DATA 
!ljh        NDAY=U1      ! DAY OF DATA SET : 3 MEANS 3RD DAY
        NDAY=int(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
!
        DO 102 J=1,NSOL
          SLAT(J)=U2
          SLNG(J)=U3
  102   CONTINUE
!
!       *** COMPUTE METERS PER DEGREE OF LATITUDE
        CONVLAT=PI*RAD/180.0
!
!       *** READ HARVARD DATA
!       OPEN(3)
        READ(3,105)DLAT0,DLNG0,THTA,NX,NY,DX,DY,NLEV,    &
        ZLEV(1),(ZLEV(I),I=2,NLEV-1),ZLEV(NLEV)
  105   FORMAT(2X,2F7.3,F8.5,2I6,/,2(F7.1),I6,/,(10F7.1))
!
        IF(NLEV.GT.MXLEV) THEN
          WRITE(NPU)'ERROR. TOO MANY LEVELS. INCREASE MXLEV. RECOMPILE.'
          NSVP=0
          RETURN
        ELSE
        ENDIF
!
!       *** 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)
!       
        READ(3,107,END=192)NN,T,ITYP,K,ITMP1,ITMP2
  107   FORMAT(2X,I5,F10.5,4I5)
!
        IF(NN.NE.NXY)GO TO 190
!
!       *** READ BOTTOM DEPTHS
        READ(3,108)(BOT(I),I=1,NXY)
  108   FORMAT(1X,8F9.3)
!
!       *** 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
        do i = 1, nxy
        write(8,*) SVP(i,8)
        end do
!
!       CALL CLOSE(3,STATUS='KEEP') ! END OF SOUND SPEED INPUT
!
        RAOLD=RA
        CONVLNG0=RAD*COS(PI*SLAT0/180.)*PI/180.0
        IPASS=1
      ELSE
      ENDIF
!
!     *** 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
!
        RJS=(C1*SIN(PSI)/(DY))+1
        JS=INT(RJS)
        IF(JS.GT.NY) THEN
          WRITE(NPU,195)
          KSVP=0
          RETURN
        ELSE
        ENDIF
!
        IF(IS.LT.1) THEN
          WRITE(NPU,194)
          KSVP=0
          RETURN
        ELSE
        ENDIF
!
        IF(JS.LT.1) THEN
          WRITE(NPU,195)
          KSVP=0
          RETURN
        ELSE
        ENDIF
!
        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*    &
                 (BOT(I3)-BOT(I1))/DY
        ELSE
          DEPTH=BOT(I4)+(DX-XOFF)*(BOT(I3)-BOT(I4))/DX+(DY-YOFF)*    &
                 (BOT(I2)-BOT(I4))/DY
        ENDIF
!
        IF(J.EQ.1) DEPTH1=DEPTH
        XOFF=XOFF/DX
        YOFF=YOFF/DY
!
!       *** 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
!
        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
!
        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))*(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
!
  140   ZSVP(NLEV+1,J)=DEPTH
        CSVP(NLEV+1,J)=CSVP(NLEV-1,J)+    &
        (CSVP(NLEV,J)-CSVP(NLEV-1,J))*(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
!
  141   IXSVP(ILYR,J)=NLEV
        ZLYR(ILYR,J)=ZSVP(NLEV,J)
        NSVP=NLEV
!
  145   CONTINUE
        ILYR=2
        RHO(ILYR,J)=BOTRHO
        RHOG(ILYR,J)=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
!
        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
!
  190 CONTINUE
      WRITE(NPU,191)
  191 FORMAT(1X,'DATA MISMATCH. NX*NY DOES NOT EQUAL NN.')
      NSVP=0
      RETURN
!
  192 WRITE(NPU,193)
  193 FORMAT(1X,'READ ERROR. CHECK INPUT SOUND SPEED FILE.')
      NSVP=0
      RETURN
!
  194 FORMAT(1X,'UBOTTOM: LONGITUDE OUTSIDE LIMITS OF SOUND SPEED DOMAIN.')
  195 FORMAT(1X,'UBOTTOM: LATITUDE OUTSIDE LIMITS OF SOUND SPEED DOMAIN.')
!ljh  196 FORMAT(1X,'LAT = ',F6.1,'  LON  = ',F6.2)
!
  200 CONTINUE
  300 CONTINUE
  400 CONTINUE
      NSVP=0
      RETURN
      END