C-----------------------------------------------------------------------------
      SUBROUTINE CFIT1
C-----------------------------------------------------------------------------
C       A horizontally uniform ocean sound speed model that fits
C       up to 100 measured C(Z) points with a smooth function
C       using the method of A.K.Paul, Radio Science .........
C       Reads C(Z) from a local file CZ.DAT .
C       The radius scaling factor G (.LE.1.) is provided by W(155).
C       The fitting error ERRMAX is provided by W(154).
C       For efficiency, keep number of data points as small as possible.

      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DOUBLE PRECISION MODC, X(100),D(100),RAD(100),A(100),B(100)
      LOGICAL THERE
      CHARACTER*9 FNAME
      COMMON//R,TH,PH,KR,KTH,KPH,RKVARS(14),TPULSE,CSTEP,DRDT(20)
      COMMON/CC/MODC(4),CS,PCST,PCSR,PCSTH,PCSPH
      COMMON/WW/ dum(10), MAXW,W(1000)
      EQUIVALENCE (EARTHR,W(1)),(W(154),ERRMAX),(W(155),G)                                         CWW30120
C
C-----------------------------------------------------------------------------
      ENTRY SETSPD
C-----------------------------------------------------------------------------
      CALL SETPSP
      RETURN
C
C-----------------------------------------------------------------------------
      ENTRY ISPEED
C-----------------------------------------------------------------------------
      CALL IPSPEED
      IF(W(150) .NE. 10.) THEN
           PRINT*,' Model check number mismatch: '
           STOP ' W(150) should be 10 if you are using model CFIT1.'
      ENDIF

      MODC(1)=5HCFIT1
      MODC(2)= w(152)

C       READ DATA FILE CZ.DAT AND COMPUTE MEAN VALUE, P0
C       NOTE THAT IN CZ.DAT, Z IS NEGATIVE FOR THE OCEAN AND THAT
C       Z MUST BE IN METERS AND C IN M/S, 10 COL. EACH

      SUM=0.
      FNAME='CZ.DAT'
6     INQUIRE (FILE=FNAME,EXIST=THERE)
      IF(THERE)THEN                       ! DOES CZ.DAT EXIST ?
           OPEN (2,FILE=FNAME)
      ELSE
           PRINT *,'I cannot find the file ',FNAME,' required by the CFI
     &T1 model'
           PRINT *,' Please enter the name of another data file:'
           READ (*,'(A9)') FNAME
           GOTO 6
      ENDIF
11    FORMAT (2F10.2)
      DO I=1,100
          READ (2,11,END=200) X(I),D(I)
          X(I)=X(I)/1000.         !CONVERT TO KM
          D(I)=D(I)/1000.         !CONVERT TO KM/S
          SUM = SUM + D(I)
          ND=I
      END DO
200   P0 = SUM / ND    ! MEAN OF DATA

C        COMPUTE A AND RAD FOR FIRST AND LAST DATA POINTS
      RAD(1)=1./(X(2)-X(1))/G  ! 1/ DISTANCE TO NEXT POINT * G
      A(1)=D(1)/P0 - 1.
      RAD(ND)=1./(X(ND)-X(ND-1))/G ! 1/ DISTANCE TO LAST POINT * G
      A(ND)=D(ND)/P0 -1.

C        COMPUTE A AND RAD FOR THE REST OF THE DATA POINTS
      ND_1=ND-1
      DO I = 2, ND_1
          A(I)=D(I)/P0-1.
          RAD(I)=1. / (MIN ( X(I+1) - X(I), X(I) - X(I-1) ))/G
      END DO

C       COMPUTE FUNCTION AT J-TH DATA POINT AND CORRECT A'S
      L=0       ! INITIALIZE ITERATION COUNT
50    ERR=0.D0  ! INITIALIZE ERROR CALCULATION
      DO 555 J=1,ND   ! LOOP ON NUMBER OF DATA POINTS
          PRODUCT=1.D0
          DO I=1,ND      ! COMPUTE FUNCTION AT J-TH POINT
              XRI = (X(J)-X(I)) * RAD(I)
              FI = 1. + A(I) / (1. + XRI * XRI)
              PRODUCT = PRODUCT * FI
          END DO
          C=P0 * PRODUCT        ! ESTIMATE OF C AT J-TH DATA POINT
          DA=(D(J)-C)/P0        ! ESTIMATE OF CORRECTION TO A(J)
          ERR=MAX(ERR,ABS(DA))  ! MAXIMUM OF ALL ERRORS
          IF ((ERR .LT. ERRMAX .AND. J.EQ.ND).OR.L.GT.99)GOTO 666
          B(J)=A(J)+.5*DA          ! B(J) = NEW A(J)
555   CONTINUE

      L=L+1        ! INCREMENT ITERATION COUNT
      DO K=1,ND    ! REPLACE OLD A'S WITH NEW ONES
          A(K)=B(K)
      END DO
      GOTO 50    ! CALCULATE FUNCTION AGAIN WITH NEW A'S
666   PRINT *,L,'  ITERATIONS TO COMPUTE CFIT1 FUNCTION'
      IF (L.GT.99)THEN
           PRINT *,'**** ITERATION NOT CONVERGING ****'
           STOP
      ENDIF

      CLOSE (2)
      RETURN

C-----------------------------------------------------------------------------
      ENTRY SPEED
C-----------------------------------------------------------------------------
      Z=R-EARTHR
      PRODUCT=1.D0
      SUM = 0.
      DO I=1,ND          ! EVALUATE FUNCTION AT POINT Z
          XRI = (Z-X(I)) * RAD(I)
          XRISQ=1./(1.+XRI*XRI)
          FI = 1. + A(I) * XRISQ
          PFIPX = -2.*XRI*RAD(I)*(FI-1.)*XRISQ
          PRODUCT = PRODUCT * FI
          SUM=SUM+PFIPX/FI
      END DO
      PPPX = PRODUCT*SUM
      C = P0 * PRODUCT
      CS=C*C           ! SOUND SPEED SQUARED

C       DERIVATIVES OF CS
      PCSR=2. * C * P0 * PPPX
      PCST=0.
      PCSTH=0.
      PCSPH=0.
      IF(W(999).EQ. 3.) PRINT *,Z,C,PCSR      ! DIAGNOSTIC PRINT
      CALL PSPEED
      RETURN
      END

