subroutine rovdh (year,rre,thetr,phir,rgvdh,rvdhg,ier)
c*
c***********************************************************************
c*
c*          "Copyright [c] CNES 98 - tous droits reserves"
c*          **********************************************
c*
c*PRO MAGLIB
c*
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.24 - V 5.0
c*VER 20.07.16 - V 6.0

c*
c*AUT spec. CNES - JC KOSIK - fevrier 2001
c*AUT port. CS-SI
c*AUT adapt. AKKA
c*
c*ROL Theme : Changements de repere
c*ROL         Calcul de la matrice de transformation du repere
c*ROL         geocentrique au repere Champ, defini avec zch 
c*ROL         suivant la verticale ascendante, ych est dans le plan zch,
c*ROL         B ou B est le champ xch complete le triedre.
c*
c*PAR year   (I) : annee fractionnaire
c*
c*PAR rre    (I) : distance radiale geocentrique (rayons terrestres)
c*PAR thetr  (I) : colatitude geocentrique (radians)
c*PAR phir   (I) : longitude geocentrique (radians)
c*
c*PAR rgvdh  (O) : matrice de rotation du repere geocentrique au 
c*ROL            : repere Champ
c*PAR rvdhg  (O) : matrice de rotation du repere Champ au repere
c*ROL            : geocentrique
c*
c*PAR ier    (O) : code de retour
c*
c*NOT ier        : sans objet
c*
c*NOT common     : util
c*

c*INF utilise    : angleg, igrf20, promat

c*
c*HST version 2.0 - 01.05.30 - Enrichissement de la maglib au CDPP
c*HST version 2.1 - 03.01.06 - corrections en compilation avec g77
c*HST version 3.0 - 01.01.07 - utilisation du champ IGRF 2005
c*HST version 4.0 - 13.10.10 - Introduction du champ IGRF10

c*HST version 5.0 - 17.02.24 - Introduction du champ IGRF15
c*HST version 6.0 - 20.07.16 - Introduction du champ IGRF20

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,thetr,phir
      double precision rgvdh(3,3),rvdhg(3,3)
      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
c     ----------------------------------
c*FON Declaration des fonctions externes
c     ----------------------------------
c
      external angleg
      double precision angleg
c
c     ---------------------------------
c*FON Declaration des variables locales
c     ---------------------------------
c
      double precision br,bt,bp
c*COM  br,bt,bp : composantes radiale, tangentielle et azimuthale du champ
c
      double precision bb
c*COM  bb : module du champ magnetique (gauss)
c
      double precision rgp(3,3),rot(3,3)
      double precision psi,theta
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 Calcul de la matrice de transformation geocentrique vers VDH
c     ------------------------------------------------------------
c
      rgp(1,1) =  dcos(thetr) * dcos(phir)
      rgp(1,2) =  dcos(thetr) * dsin(phir)
      rgp(1,3) = -dsin(thetr)
      rgp(2,1) = -dsin(phir)
      rgp(2,2) =  dcos(phir)
      rgp(2,3) =  0.d0
      rgp(3,1) =  dsin(thetr) * dcos(phir)
      rgp(3,2) =  dsin(thetr) * dsin(phir)
      rgp(3,3) =  dcos(thetr)
c
c     -----------------------------------------------------------
c*FON Calcul de l angle de rotation theta permettant d'annuler la
c*FON composante vers l'Est par rotation
c*FON la composante vers l'est Bp = D
c*FON la composante vers le Nord est Bt = - H
c*FON Recherche des composantes du champ dans tpar
c     -----------------------------------------------------------
c

      call igrf20(year,rre,thetr,phir,br,bt,bp,bb,ier)

c
      psi   = angleg(bt,bp,ier)
      theta = psi - pid
c
c     ---------------------------------------------------------
c*FON Calcul de la matrice de rotation passant de VDH au repere
c*FON *champ
c     ---------------------------------------------------------
c
      rot(1,1) =  dcos(theta)
      rot(1,2) =  dsin(theta)
      rot(1,3) =  0.d0
      rot(2,1) = -dsin(theta)
      rot(2,2) =  dcos(theta)
      rot(2,3) =  0.d0
      rot(3,1) =  0.d0
      rot(3,2) =  0.d0
      rot(3,3) =  1.d0
c
c     ------------------------------------------------------------
c*FON Calcul de la matrice de transformation pour passer du repere
c*FON geocentrique au repere champ et inverse
c     ------------------------------------------------------------
c
      call promat(rot,rgp,rgvdh,rvdhg,ier)
c
c     ****************
c     Fin de programme
c     ****************
c
      return
      end