subroutine gsfc65 (year,rre,thet,phi,br,bt,bp,bb,ier)
c*
c***********************************************************************
c*
c*          "Copyright [c] CNES 98 - tous droits reserves"
c*          **********************************************
c*
c*PRO MAGLIB
c*
c*vER 99.03.31 - V 1.0
c*
c*AUT spec. GODDARD Space Flight Center
c*AUT adap. CNES - JC KOSIK - janvier 1991
c*AUT port. CISI
c*
c*ROL Theme : Modeles de champs magnetiques
c*ROL         Calcul du champ magnetique terrestre avec les
c*ROL         coefficients de 1965 correspondant au programme
c*ROL         original de Mac Ilwain.
c*
c*PAR year (I) : annee fractionnaire (>= 1965.)
c*
c*PAR rre  (I) : distance radiale geocentrique (rayons terrestres)
c*PAR thet (I) : colatitude geocentrique (radians)
c*PAR phi  (I) : longitude geocentrique (radians)
c*
c*PAR br   (O) : composante radiale du champ magnetique le long du
c*             : meridien positive vers l'exterieur (gauss)
c*PAR bt   (O) : composante tangentielle du champ magnetique le long
c*PAR          : du meridien positive vers le sud (gauss)
c*PAR bp   (O) : composante azimuthale du champ magnetique, positive
c*PAR          : vers l'est (gauss)
c*
c*PAR bb   (O) : module du champ magnetique (gauss)
c*
c*PAR ier  (O) : code de retour
c*
c*NOT ier      : sans objet
c*
c*NOT common   : util, util2
c*
c*INF utilise  : sans objet
c*
c*HST version 1.0 - 99.03.31 - creation de la maglib au CDPP
c*HST version 2.0 - 01.06.05 - correction de commentaires de code
c*HST version 2.1 - 03.01.06 - corrections en compilation avec g77
c*
c***********************************************************************
c*
      implicit none
c
c     ---------------------------------
c*FON Declaration identificateur rcs_id
c     ---------------------------------
c
      character rcs_id*100
c
c     --------------------------
c*FON Declaration des parametres
c     --------------------------
c
      double precision year, rre, thet, phi
      double precision br, bt, bp, bb
      integer ier
c
c     ----------------------------------
c*FON Declaration des variables communes
c     ----------------------------------
c
      double precision pi,dpi,rad,deg,pid,xmu,rayt
c
c*COM pi   : constante pi (obtenue a partir de acos(-1.))
c*COM dpi  : constante 2 * pi
c*COM pid  : constante pi / 2
c*COM rad  : facteur de conversion degres  ----> radians
c*COM deg  : facteur de conversion radians ----> degres
c*COM xmu  : constante de gravitation terrestre (km**3/sec**2)
c*COM rayt : rayon equatorial terrestre (km)
c
      common/util/pi,dpi,rad,deg,pid,xmu,rayt
c
      double precision rgmt
c
c*COM rgmt : rayon geocentrique terrestre (km)
c
      common/util2/rgmt
c
c     ---------------------------------
c*FON Declaration des variables locales
c     ---------------------------------
c
      integer lg(13,13),lgt(13,13)
c*LOC lg, lgt : coefficients 1965 correspondant au programme
c*LOC         : original de Mac Ilwain
c
      integer kmax
c*LOC kmax : nombre de pas de calcul (= 10)
c
      integer n,m,jj
c*LOC n,m,jj : indices de boucles et de tableaux
c
      double precision gg(13,13),ggt(13,13)
      double precision shmit(13,13)
      double precision g(13,13)
      double precision const(13,13),fn(13),fm(13)
      double precision p(13,13),dp(13,13),sp(13),cp(13)
      double precision aor,ar,cph,ct,f1,f2,par,sph,st,t
      double precision temp,tmold,tzero,rrmag
c*LOC Variables de travail intermediaires
c
      SAVE
c
c     ---------------------------------
c*FON Affectation identificateur rcs_id
c     ---------------------------------
c
      data rcs_id /"
     >$Id$"/
c
c     --------------------------
c*FON Affectation des constantes
c     --------------------------
c
      data lg/
     >   10,-304249,-15361, 13009, 9576,-2277,  498, 709,  48,  99,3*0,
     >57748, -21616, 30002,-19870, 8028, 3595,  607,-572,  67,  29,3*0,
     >-19498,  2043, 15853, 12904, 5026, 2313,   45,  56, -88,  74,3*0,
     > -4310,  2308, -1300,  8712,-3940, -312,-2417,  75,-138,-156,3*0,  
     >  1520, -2684,    29, -2505, 2714,-1573,  -12,-244, -33, 114,3*0, 
     >    86,  1212, -1160, -1104,  799, -652,    5,  -15,  71,111,3*0,
     >  -119,  1028,   609,  -272, -124, -116,-1091,  141, -56, 10,3*0,
     >  -540,  -244,   -91,    22,  276, -211, -201,   58, 117,  4*0,
     >    69,  -122,    58,   -170,  26,  236,  -25, -160,  64, 16,3*0,
     >  -220,   156,    51,    -35, -18,   96,   121,   2, -25, 15,3*0,
     > 39*0/
c
      data lgt/
     >   100,  2059, -2907,   266,  -86, 255, -70, 6*0,
     >  -394,   602,   121, -1003,  194,  -8,  99, 6*0,
     > -1369, -1578,   -70,   163, -117, 153,  85, 6*0,
     >   649,   293,  -924,  -130,  -54, -42, 211, 6*0,
     >  -177,  -154,   318,  -548, -417, -72, 157, 6*0,
     >   304,   288,  -186,   125,   80, 164,  -9, 6*0,
     >  -139,    12,   153,   -73,   -6,  45,   6, 6*0,
     > 78*0/
c
      data kmax /10/
c
      data p(1,1)  /1.0d0/
      data cp(1)   /1.0d0/
      data dp(1,1) /0.0d0/
      data sp(1)   /0.0d0/
c
      data shmit(1,1) /0.d0/
c
      data tmold /0.d0/
c
      data tzero /1960.d0/
c
c     ******************
c     Debut de programme
c     ******************
c
      ier = 0
c
      if (shmit(1,1) .eq. -1.d0) go to 40
c
c     --------------------------------------------
c*FON Initialisation des tableaux au premier appel
c     --------------------------------------------
c
      shmit(1,1) = -1.d0
c
      do 10 n = 1, 13
         fn(n) = dble(n)
         do 15 m = 1, 13
            fm(m)      = dble(m - 1)
            const(n,m) = dble((n - 2)**2 - (m - 1)**2) 
     >                   / ((2 * n - 3) * (2 * n - 5))
15       continue
10    continue
c
      do 20 n = 2, 13
         shmit(n,1) = dble(2 * n - 3) * shmit(n-1,1) / dble(n - 1)
         jj = 2
         do 25 m = 2, n
            shmit(n,m)   = shmit(n,m-1) * sqrt(dble((n - m + 1) * jj)
     >                     / (n + m - 2))
            shmit(m-1,n) = shmit(n,m)
           jj = 1
25       continue
20    continue
c
      f1 = dble(lg(1,1))
      f2 = dble(lgt(1,1))
c
      do 30 n = 1, kmax
         do 35 m = 1, kmax
            gg(n,m)  = dble(lg(n,m))  * shmit(n,m) / f1
            ggt(n,m) = dble(lgt(n,m)) * shmit(n,m) / f2
35       continue
30    continue
c
40    continue
c
      if (year .eq. tmold) go to 55
      tmold = year
      t     = year - tzero
c
      do 45 n = 1, kmax
         do 50 m = 1, kmax
            g(n,m) = gg(n,m) + t * ggt(n,m)
50       continue
45    continue
c
c     ---------------------------------------------------------------
c*FON Calcul des coordonnees du dipole boreal et du moment magnetique
c*FON (em) - les calculs commencent ici
c     ---------------------------------------------------------------
c
55    continue
c
      ct    = cos(thet)
      st    = sin(thet)
      cph   = cos(phi)
      sph   = sin(phi)
      sp(2) = sph
      cp(2) = cph
c
      do 60 m = 3, kmax
         sp(m) = sp(2) * cp(m-1) + cp(2) * sp(m-1)
         cp(m) = cp(2) * cp(m-1) - sp(2) * sp(m-1)
60    continue
c
      rrmag   = rre * rayt / rgmt
      aor     = 1.d0 / rrmag
      ar      = aor * aor * aor
      p(2,1)  = ct
      dp(2,1) = -st
      p(2,2)  = st
      dp(2,2) = ct
c
      br = - (ar + ar) * (g(2,1) * p(2,1) + p(2,2)
     >     * (g(2,2) * cp(2) + g(1,2) * sp(2)))
      bt = ar * (g(2,1) * dp(2,1) + dp(2,2) * 
     >     (g(2,2) * cp(2) + g(1,2) * sp(2)))
      bp = ar * (g(1,2) * cp(2) - g(2,2) * sp(2)) * p(2,2)
c
      do 65 n = 3, kmax
         ar = aor * ar
         do 64 m = 1, n
            if (m .eq. n) go to 75
            if (const(n,m) .eq. 0.0d0) go to 70
            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)
            go to 80
c
70          continue
c
            p(n,m)  = ct * p(n-1,m)
            dp(n,m) = ct * dp(n-1,m) - st * p(n-1,m)
            go to 80
c
75          continue
c
            p(n,n)  = st * p(n-1,n-1)
            dp(n,n) = st * dp(n-1,n-1) + ct * p(n-1,n-1)
c
80          continue
c
            par = p(n,m) * ar
            if (m .eq. 1) go to 85
            temp = g(n,m) * cp(m) + g(m-1,n) * sp(m)
            bp   = bp - (g(n,m) * sp(m) - g(m-1,n) * cp(m)) * fm(m)
     >             * par
            go to 90
c
85          continue
c
            temp = g(n,m)
c
90          continue
c
            br = br - temp * fn(n) * par
            bt = bt + temp * dp(n,m) * ar
c
64       continue
65    continue
c
c     ---------------------------------------
c*FON Expression du champ magnetique en gauss
c     ---------------------------------------
c
      bp = bp / st / 100000.d0
      br = br / 100000.d0
      bt = bt / 100000.d0
      bb = sqrt(br * br + bt * bt + bp * bp)
c
c     ****************
c     Fin de programme
c     ****************
c
      return
      end