      subroutine fit(inop,imaxd,nskip,iter,errlim,azer,avetim,
     +epoch,nmax,nmaxt,nmaxtt,extfld,mono,prcoef) 
c ----------------------------------------------------------------------
c ** mag field fit  daniels version triangular
c ** mods by dave schmitz+j.cain to include xdot,ydot,zdot data 
c ** july 21, 1979
c ** reordered common/datar/ 25 mar 1980 jc 
c ** does not convert coefficients to gauss normalized
c ** changed id to character, etc  11/24/88 jc cain 
c ** mods by jc 4/18/99 to make some variables real*8
c-----------------------------------------------------------------------

      real*8 d,f,rdkk,fwt,sumd
      common/dd/d(1)
      integer extfld,prcoef,prp1
      character name*10,xid*72,type(11)*2
      common /datar/iskip,flatt,elong,alt,time,decl,declwt,dip,dipwt,hor
     +,horwt,x,xwt,y,ywt,z,zwt,b,bwt,xd,xdw,yd,ydw,zd,zdw,tfirst,tlast  
      common /coefs/g(61,61),h(61,61),gt(61,61),ht(61,61),gtt(61,61),
     +htt(61,61),e(3),azero,maxn
      equivalence (e(1),e1),(e(2),e2),(e(3),e3) 
      common /fld/cd,ci,ch,cf,cz,cx,cy,cddeg,cideg,wt,tm, 
     +cxd,cyd,czd,nor,nop,lat,lon,no,nonot,f(1) 
      dimension err(18,36),fno2(18,36),jerr(18),ierr(200),sig1(11), 
     +swt1(11),fno1(11) 
      parameter(itdata=100) 
      data type/'d','i','h','x','y','z','f','xd','yd','zd','*'/ 
      indx(i,j)=(i*(j-i))/2-1 
c 
c 
c      read initial coefficients
      read (5,46) xid 
      write (6,47) xid
      read(5,*) a,tzero,tf,tl 
      write(6,72) a,tzero,tf,tl 
 2    read (5,48) ln,lm,gnm,hnm,gtnm,htnm,gttnm,httnm 
      n=ln+1
      m=lm+1
      if (n.eq.0) go to 3 
      g(n,m)=gnm
      h(n,m)=hnm
      gt(n,m)=gtnm
      ht(n,m)=htnm
      gtt(n,m)=gttnm  
      htt(n,m)=httnm  
      maxn=max0(n,maxn) 
      go to 2 
c     read best set for external field
c     note that e1=ge10,e2=ge11 and e3=he11** 
    3 read (5,49) e1,e2,e3
      rewind 5
c     record starting parameters
      prp1=prcoef+1
      maxp=min0(maxn,prp1)
      write (6,48) ((n-1,m-1,g(n,m),h(n,m),gt(n,m),ht(n,m),gtt(n,m),
     +htt(n,m),m=1,n),n=mono,maxp)
      if(maxn.gt.prp1) print*,'print only to n= ',prcoef
      if(maxn.gt.prp1) print*,' maxn= ',maxn
      write (6,51) e1,e2,e3 
c      entry for successive calls using last coeffs from arrays g,h,etc 
c   note: if nmax decreased will use last higher terms as 
c   anomaly corrections unless zeroed out 
      entry fita(inop,imaxd,nskip,iter,errlim,azer) 
      azero=azer
      do 101 i=1,18 
      jerr(i)=0 
      do 101 j=1,36 
  101 err(i,j)=0. 
      tfirst=9999.
      tlast=-9999.
      j=1-mono
      nor=0 
      if(nmax.ne.0)   nor=nmax**2+j 
      no=nor
      if(nmaxt.ne.0)  nor=nor+nmaxt**2+j  
      nonot=nor 
      if(nmaxtt.ne.0) nor=nor+nmaxtt**2+j 
      if(extfld.ne.0) nor=nor+3 
      nop=nor+1 
      nopp=nop+1
      maxd=(nor*(nor+5))/2
      if(maxd.le.imaxd.and.nop.le.inop) go to 1 
      write (6,45) nop,inop,maxd,imaxd
c     call displa(' array exceeded,maxd=',maxd) 
c     call displa(4hnop=,nop) 
      stop'array exceeded'
    1 maxn=max0(nmax,maxn)
      write (6,50) nmax,nmaxt,nmaxtt,extfld,mono,nskip,iter,errlim,aveti
     +m,epoch,nop,maxd,azero
c     *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *
c     end initialization, begin data processing 
c     *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *
      do 39 itno=1,iter 
        ipr = 0
        ndata=0 
c do not convert coeffs to gauss normalization, only reduce gtt terms
        do 4 n=mono,maxn
          do 4 m=1,n 
            gtt(n,m)=0.5*gtt(n,m) 
            if(m.eq.1) go to 4  
            htt(n,m)=0.5*htt(n,m) 
    4   continue
        do 5 j=1,11 
      sig1(j)=0.
      fno1(j)=0.
   5  swt1(j)=0.
        do 6 i=1,200  
    6 ierr(i)=0 
        do 7 i=1,maxd 
    7 d(i)=0. 
      do 300 i=1,18
         do 299 j=1,36
            err(i,j) = 0.
            fno2(i,j) = 0.
 299     continue
 300  continue
        sumtm=0.0 
        line=0
        iskip=nskip 
c       read one data line
    8   call rdata(name)
c       iskip=0 ends iteration 
        if (iskip.le.0) goto23
        call fitfld (0,avetim,epoch,mono,nmax,nmaxt,nmaxtt,extfld)
c diagnostic patch 5/24/99 jc
c      if(nnn.lt.2)then
c        if(ipr.le.1)print*,'cd,ci,ch,cf,cz,cx,cy,cddeg,cideg,wt,tm'
c        print*, cd,ci,ch,cf,cz,cx,cy,cddeg,cideg,wt,tm
c        print *,'name = ',name
c        nnn=nnn+1
c        stop 'test point2'
c        endif
c       note:will not operate on data exactly zero  
        if (decl.eq.0.)go to   9
      itype=1 
      go to 15
    9   if (dip.eq.0.) go to  10
      itype=2 
      go to 15
   10   if (hor.eq.0.) go to  11
      itype=3 
      go to 15
   11   if (x.eq.0.)   go to  12
      itype=4 
      go to 15
   12   if (y.eq.0.)   go to  13
      itype=5 
      go to 15
   13   if (z.eq.0.)   go to  14
      itype=6 
      go to 15
   14   if (b.eq.0.)   go to 141
      itype=7 
      go to 15
  141   if (xd.eq.0.)  go to 142
      itype=8 
      go to 15
  142   if (yd.eq.0.)  go to 143
      itype=9 
      go to 15
  143   if (zd.eq.0.)  go to  21
      itype=10
   15 call compon(itype,avetim,epoch,mono,nmax,nmaxt,nmaxtt,extfld) 
  151 if (abs(f(nop)).lt.errlim) go to 17 
      ipr=ipr+1 
      if(ipr.gt.50) go to 17 
c*** note: may want to put in errlim for secular change 
  152 if(mod(line,29).ne.0) go to 16
        write(6,52)  
   16   write(6,530) flatt,elong,alt,time,decl,dip,hor,x,y,z,b,
     x   type(itype) 
        write(6,540) cddeg,cideg,ch,cx,cy,cz,cf,f(nop),name 
        line=line+1 
        go to (9,10,11,12,13,14,141,142,143,8) itype
c       sum deviations from computed values, weights, date
c       type counts, etc. for standard error estimates
   17   k=max(min(f(nop)/10.0+101.0,200.0),1.0) 
        ierr(k)=ierr(k)+1 
        ffw=wt*f(nop)**2
        sig1(itype)=sig1(itype)+ffw 
        fno1(itype)=fno1(itype)+1.
        swt1(itype)=swt1(itype)+wt
        sig1(11)=sig1(11)+ffw 
        fno1(11)=fno1(11)+1.
        swt1(11)=swt1(11)+wt
        sumtm=sumtm+tm
        if (itno.lt.iter) goto18
        fno2(lat,lon)=fno2(lat,lon)+wt
        err(lat,lon)=err(lat,lon)+f(nop)*wt 
c       compute squares and cross-products and
c       sum into triangular matrix d
c 
   18   k=1 
        do 20 i=1,nor 
          fwt=f(i)*wt 
          do 19 j=i,nop 
            d(k)=d(k)+fwt*f(j)  
   19     k=k+1 
   20   k=k+1 
c 
        go to (9,10,11,12,13,14,141,142,143,21),itype 
c       *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  
 21     ndata=ndata+1
c end processing one datum, go read another
        go to 8 
c       record data for stat subroutine on last iteration 
   23 print *,ndata,' observations, components rejected = ',ipr 
      if(itno.lt.iter) go to 240
        write (11) sig1,fno1,swt1,nor,nop,nmax,nmaxt,nmaxtt,extfld, 
     +  mono,maxd,maxn,avetim,epoch,azero,no
        write (11) g,h,gt,ht,gtt,htt,e
        write (11) (d(i),i=1,maxd)
        write (11) dummy
c       *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  
c       end data processing, begin solution 
c       of the least squares equations
  240   do 24 ij=1,11
          if(swt1(ij).eq.0.0) go to 24
          sig1(ij)=sqrt(sig1(ij)/swt1(ij))
   24   continue
c       compute sums for check column 
c       shirley's shorter version 
        nor3=nor+3
        k=nopp
        do 26 i=1,nor 
          ndx=i 
          sumd=d(i) 
          do 25 j=2,nor 
            if(j.le.i) ndx=ndx+nopp-j 
            ndx=ndx+1 
   25     sumd=sumd+d(ndx)
          d(k)=sumd 
   26   k=k+nopp-i
c       invert triangular matrix
c       kk picks up diagonal elements 
        kk=1
        nt1=2*nor+5 
        do 30 l=1,nor 
c         set up one complete column
          ndx=l 
          f(1)=d(l) 
          do 27 i=2,nor 
            if(i.le.l)ndx=ndx+nopp-i
            ndx=ndx+1 
            f(i)=d(ndx) 
   27     continue
          k=kk+1
          rdkk=1.0d0/d(kk)
          do 28 j=l,nop 
            d(k)=d(k)*rdkk
   28     k=k+1 
          kk=kk+nor3-l
          kjz=(l*(nt1-l))/2-nop 
          do 30 i=1,nor 
            if(i.eq.l) go to 30 
            kz=(i*(nt1-i))/2-nop
            do 29 j=l,nop 
              if(i.gt.j+1) go to 29 
              k=kz+j  
              d(k)=d(k)-f(i)*d(kjz+j) 
   29       continue  
   30   continue
c       *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  
c       end solution of the least squares equations,
c       begin estimation of the parameter corrections 
c       *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  
        write (6,55) (type(ij),sig1(ij),fno1(ij),ij=1,11) 
        write (6,56)  
        tfact=epoch-avetim
        i=1 
      if(nmax.eq.0) go to 330 
        do  33  n=mono,nmax 
          do  33  m=1,n 
c           compute time adjustment for the g(n,m) corrections
            k=indx(i,nt1) 
            kp=i+no 
            kp=indx(kp,nt1) 
            kpp=i+nonot 
            kpp=indx(kpp,nt1) 
            if (n.gt.nmaxt) go to 31
            d(k)=d(k)+d(kp)*tfact 
            if (n.gt.nmaxtt) go to 31 
            d(k)=d(k)+d(kpp)*tfact*tfact  
            d(kp)=d(kp)+2.0*d(kpp)*tfact  
c           correct g(n,m)
   31       g(n,m)=g(n,m)+d(k)  
c           record new g(n,m), correction, and the  
c           corresponding item in the check column  
            write (6,57) n,m,g(n,m),d(k),d(k+1) 
            if (m.eq.1) go to 33
            i=i+1 
c           compute time adjustment for the h(n,m) corrections
            k=indx(i,nt1) 
            kp=i+no 
            kp=indx(kp,nt1) 
            kpp=i+nonot 
            kpp=indx(kpp,nt1) 
            if (n.gt.nmaxt) go to 32
            d(k)=d(k)+d(kp)*tfact 
            if (n.gt.nmaxtt) go to 32 
            d(k)=d(k)+d(kpp)*tfact*tfact  
            d(kp)=d(kp)+2.0*d(kpp)*tfact  
c           correct h(n,m)
   32       h(n,m)=h(n,m)+d(k)  
c           record new h(n,m), correction, and the  
c           corresponding item in the check column  
            write (6,58) n,m,h(n,m),d(k),d(k+1)
   33   i=i+1
  330   if (nmaxt.eq.0) go to 36
        do 34 n=mono,nmaxt
          do 34 m=1,n 
c           correct gt(n,m) 
            k=indx(i,nt1) 
            gt(n,m)=gt(n,m)+d(k)
c           record new gt(n,m), correction, and the 
c           corresponding item in the check column  
            write (6,59) n,m,gt(n,m),d(k),d(k+1)
            if (m.eq.1) go to 34
            i=i+1 
c           correct ht(n,m) 
            k=indx(i,nt1) 
            ht(n,m)=ht(n,m)+d(k)
c           record new ht(n,m), correction, and the 
c           corresponding item in the check column  
            write (6,60) n,m,ht(n,m),d(k),d(k+1)
   34   i=i+1 
        if (nmaxtt.eq.0) go to 36 
        do 35 n=mono,nmaxtt 
          do 35 m=1,n 
c           correct gtt (n,m) 
            k=indx(i,nt1) 
            gtt(n,m)=gtt(n,m)+d(k)
c           record new gtt (n,m), correction, and the 
c           corresponding item in the check column  
            write (6,61) n,m,gtt(n,m),d(k),d(k+1) 
            if (m.eq.1) go to 35
            i=i+1 
c           correct htt(n,m)
            k=indx(i,nt1) 
            htt(n,m)=htt(n,m)+d(k)
c           record new htt(n,m), correction, and the
c           corresponding item in the check column  
            write (6,62) n,m,htt(n,m),d(k),d(k+1) 
   35   i=i+1 
c       record new e1,e2,e3 and the 
c       corresponding item in the check column
   36   if (extfld.eq.0) go to 37 
        k=indx(i,nt1) 
        e1=e1+d(k)
        write (6,63) e1,d(k),d(k+1) 
        i=i+1 
        k=indx(i,nt1) 
        e2=e2+d(k)
        write (6,64) e2,d(k),d(k+1) 
        i=i+1 
        k=indx(i,nt1) 
        e3=e3+d(k)
        write (6,65) e3,d(k),d(k+1) 
        i=i+1
c       changed 37 from avetim to avetime in order to avoid trying to
c       change the value of a parameter (otherwise run-time error) (dtm)
   37   avetime=sumtm/fno1(11)+epoch
c       convert back to true second derivatives 
        do 38 n=mono,maxn 
          do 38 m=1,n 
            gtt(n,m)=2.*gtt(n,m)
            if(m.eq.1) go to 38 
            htt(n,m)=2.*htt(n,m)
 38       continue
c       record entire array of g and h parameters 
        write (6,66) itno 
        write (6,48) ((n-1,m-1,g(n,m),h(n,m),gt(n,m),ht(n,m),gtt(n,m),  
     +  htt(n,m),m=1,n),n=mono,maxp)
      if(maxn.gt.prp1) print*,'print only to n= ',prcoef
      if(maxn.gt.prp1) print*,' maxn= ',maxn
        write (6,51) e1,e2,e3 
c     changed line to reflect change in line number 37 (dtm)
        write (6,67) avetime,epoch,tfirst,tlast
c       write file for input to next approximation  
      tf=tfirst 
      tl=tlast
      a=azero 
      tzero=epoch 
      rewind 12 
      write(12,46)xid 
      write(12,*) a,tzero,tf,tl 
        write(12,48)((n-1,m-1,g(n,m),h(n,m),gt(n,m),ht(n,m),gtt(n,m), 
     +  htt(n,m),m=1,n),n=mono,maxn)
        n=0 
        write(12,48) n-1
        write(12,49) e
   39 continue
c     record error distributions
      write (6,68) type(11) 
      do 41 jk=1,200,10 
        jl=jk+9 
        if (jk.ge.101) go to 40 
        jm=jk-101 
        go to 41
   40   jm=jk-100 
   41 write (6,69) jm,(ierr(ik),ik=jk,jl) 
c     record mean deviation for lat-long blocks 
      write (6,70) (l,l=10,90,10) 
      do 43 k=1,36
        do 42 j=1,18  
      jerr(j)=0 
      if(fno2(j,k).eq.0.) go to 42
       jerr(j)=(err(j,k)/fno2(j,k)) 
   42 continue
   43 write (6,71) k,(jerr(m),m=1,18) 
      return
   44 format(' processing ',f10.3,' components per second') 
   45 format('0array exceeded'/'  nop inop   maxd'  
     +,'  imaxd'/2i5,2i7) 
   46 format (a)
   47 format ('1',a)  
   48 format (2i3,2f16.8,4f10.5)
   49 format (6x,3f11.4)
   50 format ('0nmax  =',i5,' nmaxt =',i5,' nmaxtt=',i5,' extfld=', 
     +i5,' mono  =',i5/' nskip =',i5,' iter  =',i5,' errlim=',
     +f10.0/' avetim=',f10.0,'   epoch=',f8.0,'   nop=',i4, 
     +'   maxd=',i7,'   azero=',f7.1) 
   51 format ('0e1=',f13.4/' e2=',f13.4/' e3=',f13.4) 
   52 format ('1  lat   long    alt    time       decl    dip     h    
     + x       y        z      b    f(nop)  name'/)
c 53   format (1x,2f6.1,f7.1,f8.2,2f7.2,55x,a2/28x,2f7.2,5f7.0,f7.0, 
c     x f7.0,1x,a10)
 530  format (2f7.2,1x,f7.2,f8.2,1x,2f7.2,5f8.0,1x,a4)
 540  format(20x,11x, 2f7.2, 5f8.0, f7.1, 1x, a10) 
 54   format(' processing time is',f10.0,
     +' ns per darray element/component') 
   55 format (' 1sigma',4x,'components'/1x/(1x,a2,f5.0,4x,f7.0))
   56 format (5x,'n',2x,'m',15x,'p',14x,'dp',13x,'1.0') 
   57 format (' g ',2i3,4e20.8,f20.2) 
   58 format (' h ',2i3,4e20.8,f20.2) 
   59 format (' gt',2i3,4e20.8,f20.2) 
   60 format (' ht',2i3,4e20.8,f20.2) 
   61 format (' gtt',i2,i3,3e20.8)
   62 format (' htt',i2,i3,3e20.8)
   63 format (' e1',6x,3e20.8)  
   64 format (' e2',6x,3e20.8)  
   65 format (' e3',6x,3e20.8)  
   66 format ('1output coefficients',39x,'itno=',i4/'0')
   67 format ('0avetim=',f8.2,5x,'epoch=',f7.0,' data range=',
     +f7.1,' to',f7.1)
   68 format ('1error distribution for',3x,a2)
   69 format (i5,3x,10i6) 
   70 format ('1mean deviation for lat-long block    ',/'0',57x,9i6)
   71 format (1x,i2,3x,18i6)
   72 format('0a=',f7.1, ' tzero=',f7.1,' tf=',f7.2,
     x' tl=',f7.2/) 
      end 
