subroutine dipol (year,rr,thet,phi,brd,btd,bpd,bd,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*VER 01.05.30 - V 2.0
c*VER 03.01.06 - V 2.1
c*VER 01.01.07 - V 3.0
c*VER 13.10.10 - V 4.0

c*VER 17.02.23 - V 5.0
c*VER 20.07.15 - V 6.0

c*
c*AUT spec. CNES - JC KOSIK - janvier 1991
c*AUT port. CISI  
c*AUT adapt. AKKA 
c*
c*ROL Theme : Modeles de champs magnetiques
c*ROL         Calcul du champ dipolaire d'un dipole incline.
c*ROL         L'orientation du dipole est donnee par les termes g11 
c*ROL         et h11 du developpement en harmoniques spheriques du 

c*ROL         champ IGRF 2020 affectes de l'evolution seculaire.

c*

c*PAR year (I) : annee fractionnaire >= 2020.)

c*
c*PAR rr   (I) : distance radiale geocentrique (rayons terrestres)
c*PAR thet (I) : colatitude geocentrique (radians)
c*PAR phi  (I) : longitude geocentrique (radians)
c*
c*PAR brd  (O) : composante radiale du champ dipolaire le long du
c*             : meridien positive vers l'exterieur (gauss)
c*PAR btd  (O) : composante tangentielle du champ dipolaire 
c*             : le long du meridien, positive vers le sud (gauss)
c*PAR bpd  (O) : composante azimuthale du champ dipolaire positive
c*             : vers l'est (gauss)
c*
c*PAR bd   (O) : module du champ dipolaire (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.05.30 - correction de commentaires de code
c*HST version 2.1 - 03.01.06 - corrections en compilation avec g77
c*HST version 3.0 - 01.01.07 - adoption des coefficients de l'igrf 2005
c*HST version 4.0 - 13.10.10 - adoption des coefficients de l'igrf 2010

c*HST version 5.0 - 17.02.23 - adoption des coefficients de l'igrf 2015
c*HST version 6.0 - 20.07.15 - adoption des coefficients de l'igrf 2020

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
      double precision rr, thet, phi
      double precision brd, btd, bpd, bd
      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
      double precision r3,ar3,dt,g10,g11,h11,cph,sph,ct,st,rmag
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     Debut de programme
c     ******************
c
      ier = 0
c
c     -------------------------------------------
c*FON Initialisation a 0 des composantes du champ
c     -------------------------------------------
c
      brd = 0.0d0
      btd = 0.0d0
      bpd = 0.0d0
c
c     ----------------------
c*FON Calculs intermediaires
c     ----------------------
c
      rmag = rr * rayt/ rgmt
      r3   = rmag**3
      ar3  = 1.d0 / r3

      dt   = year - 2020.d0
      g10  = +29404.8d0 - 5.7d0  * dt
      g11  = +1450.9d0  - 7.4d0 * dt
      h11  = -4652.5d0  + 25.9d0 * dt

c
      cph =  cos(phi)
      sph =  sin(phi)
      ct  =  cos(thet)
      st  =  sin(thet)
c
c     ------------------------------------------
c*FON Calculs des composantes du champ dipolaire
c*FON    brd = composante radiale
c*FON    btd = composante tangentielle
c*FON    bpd = composante azimuthale
c     ------------------------------------------
c
      brd = -2.d0 * ar3 * (g10 * ct + (g11 * cph + h11 * sph) * st)
      btd = ar3 * (-g10 * st + (g11 * cph + h11 * sph) * ct)
      bpd = ar3 * (-g11 * sph + h11 * cph)
c
c     --------------------------------------------
c*FON Calcul du module du champ dipolaire en gauss
c     --------------------------------------------
c
      brd = brd / 100000.d0
      btd = btd / 100000.d0
      bpd = bpd / 100000.d0
      bd  = sqrt(brd * brd + btd * btd + bpd * bpd)
c
c     ****************
c     Fin de programme
c     ****************
c
      return
      end