subroutine cread (in,ty,l,maxln,a)
   !  subroutine cread (in,out,ty,l,maxln,a)
! reads schmidt coeffs (with true second derivatives)
! and sets up for gfield and/or sphrc
! mod 5/15/83 to make less work of time update (j. cain)
! only updates after l=0 to maximum nonzero gt or ht value
! expanded coefficient format: 2i3,2f16.8,f,4f10.5
! also reads from unit 2, not 12
      implicit real (kind=8) (a-h,o-z) 
      
      parameter (md=101)
      common /coeff/tg(md,md),const(md,md),fm(md),fn(md),ge(8)
      dimension g(md,md),gt(md,md),gtt(md,md)
      character (len=72) :: ai
!     integer out
      
!      real (kind=8):: g, gt, gtt ,tg, const, fm, fn, ge
!      real (kind=8):: minn, maxn, maxnt, tf, tl, tzero
!      real (kind=8) l

      if (l /= 0) then
        minn = 2
        call con
        maxn = 0
        maxnt = 0
        read (in,9) ai 
        read (in,*) a,tzero,tf,tl
        do
          read (in,10) ln,lm,gnm,hnm,gtnm,htnm,gttnm,httnm
          n = ln + 1
          m = lm + 1
          if (lm*ln == 0) minn = 1
          if (ln < 0) exit
          maxn = max(n,maxn)
          if (gtnm /= 0. .or. htnm /= 0.) maxnt = max(n,maxnt)
          g(n,m)   = gnm
          gt(n,m)  = gtnm
          gtt(n,m) = gttnm
          if (lm == 0) cycle
          g(lm,n)   = hnm
          gt(lm,n)  = htnm 
          gtt(lm,n) = httnm
        enddo
        maxln = maxn - 1
        read (in,16) ge
        close(in)
!       write (out,11) ai,maxln,maxnt-1,a,tzero,tf,tl 
!       if (l <= 1) then
!         write (out,12)
!         do n = minn, maxn
!           ln = n - 1
!           do m = 1, n
!             lm = m - 1
!             if (m == 1) then
!               write (out,14) ln,lm,g(n,m),gt(n,m),gtt(n,m)
!             else
!               write (out,13) ln,lm,g(n,m),g(lm,n),gt(n,m),gt(lm,n),gtt(n,m),gtt(lm,n)
!             endif
!           enddo
!         enddo
!         write (out,44) ge
!       endif

        do n = 1, maxn 
          do m = 1, maxn 
            call cain_convrt(g(n,m),n,m,1)
!print*, 'ap. cain_convrt 1: '
            tg(n,m) = g(n,m)
            call cain_convrt(gt(n,m),n,m,1)
!print*, 'ap. cain_convrt 2: '
            gtt(n,m) = 0.5*gtt(n,m)
            call cain_convrt(gtt(n,m),n,m,1)
!print*, 'ap. cain_convrt 3: '
          enddo
        enddo
    
 !      if(ty > tl .or. ty < tf) write (out,15) ty
        l = 0 
      endif

      t = ty - tzero
    
      do n = 1, maxnt
        do m = 1, maxnt
          tg(n,m) = g(n,m) + t*(gt(n,m)+t*gtt(n,m))
        enddo
      enddo
      return

! two formats for coeffs
!  10 format (2i3,6f11.5)
   10 format (2i3,2f16.8,4f10.5)
    9 format (a)
!  11 format ('0',a,/2x,'mxln=',i3,2x,'mxlnt=',i3,2x,'a=',f6.1,2x,'tzero=',f7.1,2x,'tf=',f7.1,2x,'tl=',f7.1) 
!  12 format ('0 n  m',8x,'g',10x,'h',10x,'gt',9x,'ht',8x,'gtt',8x,'htt')
!  13 format (2i3,2f11.2,2f11.3,2f11.4) 
!  14 format (2i3,f11.2,11x,f11.3,11x,f11.4)
!  15 format ('0**warning,',f7.1,' is outside time limits') 
   16 format (6x,8f11.4)
!  44 format (' first externals',/6x,8f11.4)

      end 


      subroutine cain_convrt (g,i,l,k) 
! cain_convrt pour ne pas confondre avec convrt de la spicelib
! k=1 converts schmidt to gauss,k>1 converts gauss to schmidt
      implicit real (kind=8) (a-h,o-z) 

      parameter (md=101)
      dimension s(md,md)
      logical next
      save s
      data next/.false./

!print*,' cain_convrt g,i,l,k: ', g,i,l,k
      if (.not. next) then 
        next = .true.
        s(1,1) = -1.
        do n = 2, md
!         s(n,1) = s(n-1,1)*real(2*n-3)/real(n-1)   ! Cain 
          s(n,1) = s(n-1,1)*real(2*n-3, kind=8)/real(n-1, kind=8)  !dpo
          s(1,n) = 0.
          j = 2 
          do m = 2, n
!          s(n,m) = s(n,m-1)*sqrt(real((n-m+1)*j)/real(n+m-2))  !Cain
           s(n,m) = s(n,m-1)*dsqrt(real((n-m+1)*j, kind=8)/real(n+m-2, kind=8)) !dpo
           s(m-1,n) = s(n,m)
           j = 1
          enddo
        enddo 
      endif

      if (k > 1) then
        g = g/s(i,l)
      else
        g = g*s(i,l)
      endif
    
      return
      end 


      subroutine con
! sets up constants for use by sphrc (normally called by one
! of the subroutines that sets up the coefficients tg)

      implicit real (kind=8) (a-h,o-z) 

      parameter (md=101)
      common/coeff/tg(md,md),c(md,md),fm(md),fn(md),ge(8)
!     data c/3721*0./  ! assume computer core preset to zero or needed

      fm(1) = 0.
      do n = 2, md
        fm(n) = n-1
        fn(n) = n
        do m = 1, n
          c(n,m) = real( (n-2)**2 - (m-1)**2 ,kind=8)/real( (2*n-3)*(2*n-5), kind=8)
        enddo
      enddo
      
      return
      end 


      subroutine sphrc
! assume computer core is initiallized to zero
! otherwise should zero p,dp,sp arrays
! this version best if frequent dipole only computation.
! otherwise would be better to put n=2 in main loop
      implicit real (kind=8) (a-h,o-z) 

      parameter (md=101)
      common /gcom/ st,ct,sph,cph,aor,bt,bp,br,nmax,sind,cosd
      common /coeff/ g(md,md),const(md,md),fm(md),fn(md),ge(8)
      dimension p(md,md),dp(md,md),sp(md),cp(md)
      dimension pp(0:md,0:md),dpp(0:md,0:md),s(0:md,0:md),q(0:md,0:md)


      ar2 = aor**2
      ar = ar2*aor
      gc = g(2,2)*cph + g(1,2)*sph
      br = -ar2*g(1,1) - 2.*ar*(g(2,1)*ct+gc*st)
      bt = ar*(gc*ct-g(2,1)*st) 
      bp = ar*(g(1,2)*cph-g(2,2)*sph)
      if (nmax /= 2) then
        bp = bp*st
        p(2,1)  =  ct
        dp(2,1) = -st
        p(2,2)  =  st
        dp(2,2) =  ct
        sp(2) = sph
        cp(2) = cph
        p(1,1) = 1.
        do n = 3, nmax
          sp(n) = sph*cp(n-1) + cph*sp(n-1)
          cp(n) = cph*cp(n-1) - sph*sp(n-1)
          dp(n,1) = ct*dp(n-1,1) - st*p(n-1,1) - const(n,1)*dp(n-2,1)
          stm = g(n,1)*dp(n,1)
          spm = 0.
          p(n,1) = ct*p(n-1,1) - const(n,1)*p(n-2,1)
          srm = -g(n,1)*p(n,1)
          p(n,n) = st*p(n-1,n-1)
          dp(n,n) = fm(n)*ct*p(n-1,n-1)
          do m = 2, n
            if (n /= m) then
              p(n,m) = ct*p(n-1,m) - const(n,m)*p(n-2,m)
              dp(n,m) = ct*dp(n-1,m) - st*p(n-1,m) - const(n,m)*dp(n-2,m)
            endif
            gc = g(n,m)*cp(m) + g(m-1,n)*sp(m)
            stm = stm + gc*dp(n,m)
            spm = spm + (g(m-1,n)*cp(m)-g(n,m)*sp(m))*fm(m)*p(n,m)
            srm = srm - gc*p(n,m)
         enddo
          ar = aor*ar
          bt = bt + stm*ar
          bp = bp + spm*ar
          br = br + srm*fn(n)*ar
         enddo
          if(st == 0.) st = 1.0e-9   ! Cain ori.
! print*, 'sphcar >>> st: ', 
        bp = bp/st
      endif
 ! assume first order external usually non-zero
 
    roa = 1./aor
    c2ph = 2.*cph*cph - 1.
    s2ph = 2.*sph*cph
 
    q(1,0) = ge(1)
    q(1,1) = ge(2)
    s(1,1) = ge(3)
    q(2,0) = ge(4)
    q(2,1) = ge(5)
    s(2,1) = ge(6)
    q(2,2) = ge(7)
    s(2,2) = ge(8)
    
    pp(1,0) = p(2,1)
    pp(1,1) = p(2,2)
    pp(2,0) = p(3,1)
    pp(2,1) = p(3,2)
    pp(2,2) = p(3,3)

    dpp(1,0) = dp(2,1)
    dpp(1,1) = dp(2,2)
    dpp(2,0) = dp(3,1)
    dpp(2,1) = dp(3,2)
    dpp(2,2) = dp(3,3)
    

    e_r =   q(1,0)*pp(1,0)+(q(1,1)*cph+s(1,1)*sph)*pp(1,1) & 
         + (q(2,0)*pp(2,0)+(q(2,1)*cph+s(2,1)*sph)*pp(2,1)+(q(2,2)*c2ph+s(2,2)*s2ph)*pp(2,2))*2.0*roa
    e_t =   q(1,0)*dpp(1,0)+(q(1,1)*cph+s(1,1)*sph)*dpp(1,1) & 
         + (q(2,0)*dpp(2,0)+(q(2,1)*cph+s(2,1)*sph)*dpp(2,1)+(q(2,2)*c2ph+s(2,2)*s2ph)*dpp(2,2))*roa
    e_p =  (s(1,1)*cph-q(1,1)*sph) & 
         + ((s(2,1)*cph-q(2,1)*sph)*pp(2,1)+(s(2,2)*c2ph-q(2,2)*s2ph)*2.0*pp(2,2))*roa/st

    br = br - e_r
    bt = bt - e_t
    bp = bp - e_p
! print*, 'sphrc -> br, bt, bp: ', br, bt, bp
      return
      end 


      subroutine gfield(dlat,dlong,alt,mln,a,btt, bpp, brr, x,y,z,f)

      implicit real (kind=8) (a-h,o-z) 

      common/gcom/st,ct,sph,cph,aor,bt,bp,br,nmax,sind,cosd 
      save re,a2,a4,a2b2,b2,a4b4
   !  data re,tlast,rp/2*0.,3374.9/   Cain
      data re,tlast,rp/2*0.,3376.2/  !PCK00010.TPC

! print*, 'gfield -> dlat,dlong,alt: ', dlat,dlong,alt
      if (re /= 3396.19) then  !Cain: 3396.9, PCK00010.TPC: 3396.19
        re = 3396.19
        a2 = re**2
        a4 = a2**2
!        flat = 1.-1./298.25
        flat=rp/re
        b2 = (re*flat)**2
        a2b2 = a2*(1.-flat**2)
        a4b4 = a4*(1.-flat**4)
!       drad = atan(1.)/45.   ! Cain orig.
      end if
      drad = atan(1.)/45.    !dpo
      sinla = dsin(dlat*drad)
      rlong = dlong*drad
      cph = dcos(rlong)
      sph = dsin(rlong)
!print*, 'gfield drad rlong, cph, sph: ',drad,  rlong, cph, sph
      sinla2 = sinla**2
      cosla2 = 1. - sinla**2
      den2 = a2 - a2b2*sinla**2
      den = dsqrt(den2)
      fac = ((alt*den+a2)/(alt*den+b2))**2
      r = dsqrt(alt*(alt+2.*den)+(a4-a4b4*sinla2)/den2)
      ct = sinla/dsqrt(fac*cosla2+sinla2)
      st = dsqrt(1.-ct**2)
! print*, 'r, ct, st: ', r, ct, st
      nmax = mln + 1
      aor = a/r
      call sphrc
      f = dsqrt(bt**2+bp**2+br**2)

      btt=bt
      bpp=bp
      brr=br
! transforms to geodetic directions 
      sind = sinla*st-dsqrt(cosla2)*ct
      cosd = dsqrt(1.-sind**2)
      x = -bt*cosd-br*sind
      y = bp
      z = bt*sind-br*cosd
!print*,'### gfield dlat,dlong,alt, br, bt, bp : ', dlat,dlong,alt,br,bt,bp
      return
      end