subroutine corgm (magin,year,magout,indgm,indval,tilt,rb,
     >                  rggsm,rgsmg,rgdip,tetdip,phidip,
     >                  rre,theto,phio,rmax,np,
     >                  tr,tthet,tphi,tdr,tdtet,tdphi,ieq,dreq,
     >                  dteteq,dphieq,corthet,corlat,corlon,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 2017.02.23 - V 5.0
c*
c*AUT spec. CNES - JC KOSIK - janvier 1991
c*AUT port. CISI
c*AUT adapt. AKKA
c*
c*ROL Theme : Calculs de geophysique
c*ROL         Calcul des coordonnees geomagnetiques corrigees,  
c*ROL         latitude et longitude, tenant compte de champ magnetique
c*ROL         interne (de type DGRF 2000 ou DGRF 2005 ou DGRF 2010 ou IGRF 2015) 
c*ROL         et externe (de type Tsyganenko 1897 et Tsyganenko 1989 et Kosik 97).
c*
c*PAR magin   (I) : type de champ magnetique interne
c*
c*PAR year    (I) : annee fractionnaire
c*
c*PAR magout  (I) : type de champ magnetique externe
c*
c*PAR indgm   (I) : type d'indice geomagnetique Kp ou Ae
c*
c*PAR indval  (I) : indice geomagnetique : niveau d'amplitude du champ
c*
c*PAR tilt    (I) : angle de tilt (radians)
c*
c*PAR rb      (I) : distance subsolaire de la magnetopause
c*PAR             : (rayons terrestres) rb = 10. * re
c*
c*PAR rggsm   (I) : matrice (3,3) de passage du repere geographique au
c*PAR             : repere solaire magnetospherique
c*PAR rgsmg   (I) : matrice (3,3) de passage du repere solaire
c*PAR             : magnetospherique au repere geographique
c*
c*PAR rgdip   (I) : matrice (3,3) de passage du repere geographique
c*PAR             : au repere dipolaire
c*
c*PAR tetdip  (I) : colatitude geocentrique du dipole (radians)
c*PAR phidip  (I) : longitude geocentrique du dipole (radians)
c*
c*PAR dir     (I) : direction du trace des lignes de champ
c*
c*PAR rre     (I) : distance radiale geocentrique (rayons terrestres)
c*
c*PAR theto   (I) : colatitude geocentrique (radians)
c*PAR phio    (I) : longitude geocentrique (radians)
c*
c*PAR rmax    (I) : distance geocentrique maximum autorisee dans le trace
c*PAR             : (rayons terrestres)
c*
c*PAR np      (O) : nombre de points calcules
c*
c*PAR tr      (O) : tableau des distances geocentriques des point de la
c*PAR             : ligne de champ
c*
c*PAR tthet   (O) : tableau des colatitudes des points calcules de la
c*PAR             : ligne de champ
c*PAR tphi    (O) : tableau des longitudes des points calcules de la
c*PAR             : ligne de champ
c*
c*PAR tdr     (I) : tableau des distances dipolaires (r. terrestres)
c*
c*PAR tdtet   (I) : tableau des colatitudes dipolaires (radians)
c*PAR tdphi   (I) : tableau des longitudes dipolaires (radians)
c*
c*PAR ieq     (O) : indicateur de passage a l'equateur
c*
c*PAR dreq    (O) : distance radiale geocentrique du point equatorial
c*PAR             : de la ligne de champ (rayons terrestres)
c*PAR dteteq  (O) : colatitude geocentrique du point equatorial     
c*PAR             : de la ligne de champ (radians) 
c*PAR dphieq  (O) : longitude geocentrique du point equatorial
c*PAR             : de la ligne de champ (radians) 
c*
c*PAR corthet (O) : colatitude geomagnetique corrigee (radians)
c*PAR corlat  (O) : latitude geomagnetique corrigee (radians)
c*PAR corlon  (O) : longitude geomagnetique corrigee (radians)
c*
c*PAR ier     (O) : code de retour
c*
c*NOT magin      : 1 ---> champ dipolaire
c*NOT magin      : 2 ---> champ DGRF 2000
c*NOT magin      : 3 ---> champ DGRF 2005
c*NOT magin      : 4 ---> champ DGRF 2010
c*NOT magin      : 5 ---> champ IGRF 2015
c*
c*NOT year       : year doit etre >= 2000.
c*
c*NOT magout     : 0 = pas de champ externe
c*NOT magout     : 1 = Tsyganenko 87
c*NOT magout     : 2 = Tsyganenko 89
c*NOT magout     : 3 = Kosik 97
c*
c*NOT indgm      : 1 ---> indice geomagnetique Kp
c*NOT indval     : 1 ---> Kp =  0  , 0+
c*NOT indval     : 2 ---> Kp =  1- , 1 , 1+
c*NOT indval     : 3 ---> Kp =  2- , 2 , 2+
c*NOT indval     : 4 ---> Kp =  3- , 3 , 3+
c*NOT indval     : 5 ---> Kp =  4- , 4 , 4+
c*NOT indval     : 6 ---> Kp >  5-
c*
c*NOT indgm      : 2 ---> indice geomagnetique Ae
c*NOT indval     : 1 ---> Ae =    0 - 50
c*NOT indval     : 2 ---> Ae =   50 - 100
c*NOT indval     : 3 ---> Ae =  100 - 150
c*NOT indval     : 4 ---> Ae =  150 - 250
c*NOT indval     : 5 ---> Ae =  250 - 400
c*NOT indval     : 6 ---> Ae >= 400
c*
c*NOT dir        : +1 = trace vers les altitudes plus elevees
c*NOT            :      (depuis la surface vers l hemisphere oppose)
c*
c*NOT rfin       : point d'arret de la ligne de champ choisi
c*NOT            : pour la surface rfin=1
c*
c*NOT rmax       : point d arret des lignes de champ ouvertes
c*
c*NOT ieq        : 1 = passage a l'equateur
c*NOT            : 0 = non passage a l'equateur
c*
c*NOT ier        : 0 = OK
c*NOT            : 1 = la ligne de champ ne croise pas le
c*NOT            :     plan equatorial
c*
c*NOT common     : util, iter
c*
c*INF utilise    : magtot, iterat, mpause, equdip, angdip, corr
c*
c*HST version 1.0 - 99.03.31 - creation de la maglib au CDPP
c*HST version 2.0 - 01.05.30 - Ajout indice Kosik 97
c*HST version 2.1 - 03.01.06 - corrections en compilation avec g77
c*HST version 3.0 - 01.01.07 - adoption des coefficients definitifs
c*HST                          du champ 2000
c*HST version 4.0 - 13.10.10 - Adoption des coefficients definitifs
c*HST                          du champ 2005 et des coefficients provisoires
c*HST                          2010
c*HST version 5.0 - 2017.02.23 - Adoption des coefficients definitifs
c*HST                            du champ 2010 et des coefficients provisoires
c*HST                            2015
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
      integer magin
      double precision year
      integer magout, indgm, indval
      double precision tilt, rb
      double precision rggsm(3,3), rgsmg(3,3), rgdip(3,3)
      double precision tetdip, phidip, dir
      double precision rre, theto, phio
      double precision rmax
      integer np
      double precision tr(500), tthet(500), tphi(500)
      double precision tdr(500), tdtet(500), tdphi(500)
      integer ieq
      double precision dreq, dteteq, dphieq
      double precision corlat, corthet, corlon
      integer ier
c
c     ----------------------------------
c*FON Declaration des variables communes
c     ----------------------------------
c
      double precision y(3),yold(3)
      double precision br,bt,bp,b,st,sgn,ds
c*COM y        : composantes spheriques
c*COM yold     : composantes spheriques precedantes
c*COM br,bt,bp : composantes radiale, tangentielle et azimuthale du champ
c*COM b        : module du champ
c*COM st       : sinus theta
c*COM sgn      : sens du calcul (+1, -1)
c*COM ds       : intervalle d'integration le long de la ligne de champ
c
      common/iter/y,yold,br,bt,bp,b,st,sgn,ds
c
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 variables locales
c     ---------------------------------
c
      integer l
c*LOC l : pas d'integration
c
      integer npeq
c*LOC npeq : indice courant du point de la ligne de champ
c
      integer ipop
c*LOC ipop : indicateur de positionnement par rapport a la magnetosphere
c*
      integer ier1,ier2,ier3,ier4,ier5,ier6
c*LOC ier1,ier2,ier3,ier4,ier5,ier6 : codes retour des modules appeles
c
      double precision brm,btm,bpm,bm
c*LOC brm,btm,bpm,bm : coordonnees spheriques
      double precision cglond,singml
c*LOC cglond,singml : variables intermediaires pour la determination du ds
c
      double precision altgc
c*LOC altgc : altitude en km du 1er point de la ligne de champ
c
      double precision epsr
c*LOC epsr : epsilon pour validite de l'altitude
c
      double precision r,thet,phi,rp,thetp,phip
c*LOC r,thet,phi,rp,thetp,phip : coordonnees spheriques
c
      double precision tetdp,phidp,theti
c*LOC tetdp,phidp,theti : coordonnees spheriques
c
      double precision rdp,prod
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
c     -------------------------
c*FON Initialisations de depart
c     -------------------------
c
      ier = 0
      ier1 = 0
      ier2 = 0
      ier3 = 0
      ier4 = 0
      ier5 = 0
      ier6 = 0
      corlat = 999.d0 * rad
      corlon = 999.d0 * rad
c
      dir   = +1.d0
      sgn   = dir 
      ieq   = 0      
      altgc = (rre - 1.d0)
      epsr  = 0.01d0 
c 
c     -------------------
c*FON Determination du ds 
c     -------------------
c
      cglond = cos(phio - phidip)
      singml = cos(tetdip) * cos(theto) + sin(tetdip)
     >         * sin(theto) * cglond
      ds     = .06d0 * rre / (1.d0 - singml**2) - 0.058d0
c
      if (ds .gt. 1.d0) then
         ds = 0.3d0
      endif
c
      if (dir .lt. 0.d0) then
         if (altgc .gt. epsr .and. altgc .lt. 2.0d0) then
            ds = altgc / 20.d0
         else if (altgc .le. epsr)then
            ier = 1
         endif
      endif
c
      if (ier .eq. 1) goto 20
c
      r    = rre  
      thet = theto
      phi  = phio 
      l    = 0
      ieq  = 0 
c
c     --------------------------------------
c*FON Boucle de traitement tant que r > rfin
c*FON et que     = 0
c     --------------------------------------
c
10    continue 
c
c     ---------------------
c*FON Calcul du champ total
c     ---------------------
c
      st = sin(thet)
c
      call magtot(magin,year,magout,indgm,indval,tilt,rggsm,rgsmg,r,
     >            thet,phi,brm,btm,bpm,bm,ier1)
c
c     ----------------------------------------
c*FON Integration le long de la ligne de champ 
c     ----------------------------------------
c
      if (l .eq. 0) then
          yold(1) = 0.d0
          yold(2) = 0.d0
          yold(3) = 0.d0
          sgn = sign(1.d0, brm * dir)
      endif
      l = l + 1
c
      br = brm
      bt = btm
      bp = bpm
      b  = bm
      y(1) = r
      y(2) = thet
      y(3) = phi
      call iterat(l,ier2) 
      r    = y(1)
      thet = y(2)
      phi  = y(3)
      rp    = yold(1)
      thetp = yold(2)
      phip  = yold(3)
c
      if (l .gt. 500 .or. rp .gt. rmax) ier = 1
c
      if (ier .eq. 1) go to 20 
c
c     -----------------------------------------------
c*FON Calcul et test d'appartenance a la magnetopause
c*FON de type Shabansky
c     -----------------------------------------------
c
      call mpause(rb,rggsm,rp,thetp,phip,ipop,ier3) 
c
      if (ipop .eq. 2) then
         ier = 1
         go to 20 
      endif
c
      np = l
      tr(np) = rp
c
      if (thet .gt. pi) then 
         thet = dpi - thet
         phi  = pi  + phi
      endif
c
      if (thetp .gt. pi) then
         thetp = dpi - thetp
         phip  = pi  + phip
      endif
c
c     -----------------------------------------------------
c*FON Calcul de la colatitude et de la longitude dipolaires
c     -----------------------------------------------------
c
      call angdip(rgdip,rp,thetp,phip,rdp,tetdp,phidp,ier4)
c
      tdr(np)   = rdp 
      tthet(np) = thetp
      tdtet(np) = tetdp
      tphi(np)  = phip
      tdphi(np) = phidp
      prod      = (tdtet(np-1) - pid) * (tdtet(np) - pid)
c
      if (np .gt. 8 .and. prod .lt. 0.d0) then
         ieq  = 1
         npeq = np
c
c     --------------------------------------------------
c*FON Calcul de l'intersection de la ligne de champ avec
c*FON l'equateur geomagnetique dipolaire incline
c     --------------------------------------------------
c
         call equdip(npeq,tdtet,tdphi,tdr,dteteq,dphieq,dreq,ier5)
         theti = tdtet(1)
c
c        --------------------------------------------------------
c*FON    calcul de la colatitude geomagnetique et de la longitude
c*FON    geomagnetique corrigees
c        --------------------------------------------------------
c
         call corr(altgc,theti,dreq,dphieq,corlat,corlon,ier6)
c        
         corthet = pid - corlat
c
      endif
c
      if (ieq .eq. 0. and .ier .eq. 0) go to 10
c
c     -----------------------------------------------------
c*FON Calcul de la colatitude et de la longitude dipolaires
c     -----------------------------------------------------
c
      call angdip(rgdip,rp,thetp,phip,rdp,tetdp,phidp,ier4)
c
      tdtet(np) = tetdp
      tdphi(np) = phidp
      tdr(np)   = rdp
20    continue
c
c     ****************
c     Fin de programme
c     ****************
c
999    continue
      return
      end