!
! XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
! XX                                                                  XX
! XX  Complements a la tsygalib pour le calcul de lignes de forces    XX
! XX  Utilisee pour comagneto et tous les animagneto et anikp         XX
! XX  P. Robert, 1995 - derniere modif Avril 98                       XX
! XX                  - ajou de T89c et T96, Juin 2002                XX
! XX                  - ajou de T01 et T04, Novembre 2005             XX
! XX                  - ajou du  calcul des Kp  a partir de la base   XX
! XX                    des indices du CETP mise en format ascii,     XX
! XX                    Juillet 2003                                  XX
! XX                  - Portabilite gfortran, Juin 2010               XX
! XX                                                                  XX
! XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
!
!
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine calinca(ex,ikp,parmod,dir,xi,yi,zi,xl,yl,zl,nl)
!
!     *****************************************************************
!     calcule de la ligne de champ de Tsyganenko
!     point de depart en cartesien
!     *****************************************************************
!
      character*(*) ex
      real (kind=4), dimension(10) :: parmod
      real (kind=4), dimension(1000) :: xl,yl,zl

      external t87s,t87l,t89c,t96_01,t01_01,t04_s
      external dish,exzero,IGRF_GSM,DIP
!
!     verification des parametres en entres, ikp quelconque pour t96_01
!     t01_01 t04_s
!
      call check_ex(ex)
      call check_ikp(ex,ikp)
!
      rlim=80.
      r0=1.
      xlim=20.
      yzlim=40.
!
!
      if(ex(1:4).eq.'t89c' ) then
                       call trace2(xi,yi,zi,dir,rlim,xlim,yzlim, &
                                   r0,ikp,parmod,t89c,IGRF_GSM, &
                                   xf,yf,zf,xl,yl,zl,nl)
                       return
                       endif
!
      if(ex(1:4).eq.'t87s') then
                       call trace2(xi,yi,zi,dir,rlim,xlim,yzlim, &
                                   r0,ikp,parmod,t87s,IGRF_GSM, &
                                   xf,yf,zf,xl,yl,zl,nl)
                       return
                       endif
!
      if(ex(1:4).eq.'t87l') then
                       call trace2(xi,yi,zi,dir,rlim,xlim,yzlim, &
                                   r0,ikp,parmod,t87l,IGRF_GSM, &
                                   xf,yf,zf,xl,yl,zl,nl)
                       return
                       endif
!
      if(ex(1:6).eq.'t96_01') then
                       call trace2(xi,yi,zi,dir,rlim,xlim,yzlim, &
                                   r0,ikp,parmod,t96_01,IGRF_GSM, &
                                   xf,yf,zf,xl,yl,zl,nl)
                       return
                       endif
!
      if(ex(1:6).eq.'t01_01') then
                       call trace2(xi,yi,zi,dir,rlim,xlim,yzlim, &
                                   r0,ikp,parmod,t01_01,IGRF_GSM, &
                                   xf,yf,zf,xl,yl,zl,nl)
                       return
                       endif
!
      if(ex(1:5).eq.'t04_s') then
                       call trace2(xi,yi,zi,dir,rlim,xlim,yzlim, &
                                   r0,ikp,parmod,t04_s,IGRF_GSM, &
                                   xf,yf,zf,xl,yl,zl,nl)
                       return
                       endif
!
      xlim=100.
      yzlim=100.
!
      if(ex(1:4).eq.'dish') then
                       call trace2(xi,yi,zi,dir,rlim,xlim,yzlim, &
                                   r0,ikp,parmod,dish,IGRF_GSM, &
                                   xf,yf,zf,xl,yl,zl,nl)
                       return
                       endif
!
      if(ex(1:4).eq.'igrf') then
                       call trace2(xi,yi,zi,dir,rlim,xlim,yzlim, &
                                   r0,ikp,parmod,exzero,IGRF_GSM, &
                                   xf,yf,zf,xl,yl,zl,nl)
                       return
                       endif
!
      if(ex(1:4).eq.'dipo') then
                       call trace2(xi,yi,zi,dir,rlim,xlim,yzlim, &
                                   r0,ikp,parmod,exzero,DIP, &
                                   xf,yf,zf,xl,yl,zl,nl)
                       return
                       endif
!
      print*, '*** calinca: only t89c, t87s, t87l, t96_01, t01_01, t04_s, &
                            &dish, igrf, dipo'
!
      return
      end
!
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine calinpo(ex,ikp,parmod,ri,teti,phii,rf,tetf,phif,xl, &
           yl,zl,nl)
!
!     *****************************************************************
!     calcule de la ligne de champ de Tsyganenko
!     point de depart en polaire
!     *****************************************************************
!
      character*(*) ex
      real (kind=4), dimension(10) :: parmod
      real (kind=4), dimension(1000) :: xl,yl,zl
!
!
      pi=acos(-1.)
      pis2=pi/2.
!
      if(teti.lt.pis2) then
                          dir= 1.
                          else
                          dir=-1.
                          endif
!
      call sphcar(ri,teti,phii,xi,yi,zi,1)
!
      call calinca(ex,ikp,parmod,dir,xi,yi,zi,xl,yl,zl,nl)
!
      call sphcar(rf,tetf,phif,xl(nl),yl(nl),zl(nl),-1)
!
      return
      end
!
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine cokp(ifc,iyear,imon,iday,ih,rkp)
!
!     *****************************************************************
!     calcule le kp a partir de la base d'indices
!     *****************************************************************
!
      character*80 dirdat
      character*14 file
      integer kp8(8)
      save dirdat,file,kp8
!
      dirdat='./data/Kp_data/'
!
!     exemple de nom de fichier: kp2003_wdc.txt
!
      write (file, 100) iyear
!
      call strlen(dirdat,nc)
!
      print*, 'cokp: opening file ',dirdat(1:nc)//file, &
              ' with file code ',ifc
!
      close(ifc)
      open(ifc,file=dirdat(1:nc)//file,err=30)
!
  10  continue
      read(ifc,200,end=20) iyy,imm,idd,kp8
!
      if(iyy.ge.32) then
                    iyy=iyy+1900
                    else
                    iyy=iyy+2000
                    endif

      if(iyy.ne.iyear) then
                       print*, '*** cokp: year in file is ',iyy
      stop 'COKP           *** ABORTED ! Year not found in Kp data  ***'
                       endif
      if(imm.ne.imon) go to 10
      if(idd.ne.iday) go to 10
!
      if(ih.lt.0.or.ih.gt.23) &
      stop 'COKP           *** ABORTED ! hour must be >0 and <24    ***'
      rkp= float(kp8(ih/3 +1))/10.
      close(ifc)
      return
!
   20 continue
      print*, '*** cokp: end of Kp data file reached without good epoch'
      print*, '          set to not valib value = -1.'
      rkp= -1.
      return
!
   30 continue
      print*, '*** cokp: attempt to read a bad/inexistent Kp data file'
      print*, '          set to not valib value = -1.'
      rkp= -1.
      return
!
  100 format('kp',i4.4,'_wdc.txt')
  200 format(3i2,6x,8i2)
!
      end
!
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine convkptsy(ex,rkp,ikp)
!
!     *****************************************************************
!     calcule le ikp de Tsyganenko 87s, 87l, 89c a partir du rkp
!     *****************************************************************
!
      character*(*) ex
!
!
      ikp= -1
      if(rkp.lt.0.) return
!
! *** calcul de ikp pout T87s
!
      if(ex(1:4).ne.'t87s') go to 10
!
      ikp=8
!
      if(rkp.lt.5.33) ikp=7
      if(rkp.lt.4.66) ikp=6
      if(rkp.lt.3.66) ikp=5
      if(rkp.lt.2.66) ikp=4
      if(rkp.lt.2.00) ikp=3
      if(rkp.lt.1.33) ikp=2
      if(rkp.lt.0.66) ikp=1
!
   10 continue
!
! *** calcul de ikp pout T87l
!
      if(ex(1:4).ne.'t87l') go to 20
!
      ikp=6
!
      if(rkp.lt.4.66) ikp=5
      if(rkp.lt.3.66) ikp=4
      if(rkp.lt.2.66) ikp=3
      if(rkp.lt.1.66) ikp=2
      if(rkp.lt.0.66) ikp=1
!
   20 continue
!
! *** calcul de ikp pout T89c
!
      if(ex(1:4).ne.'t89c') go to 30
!
      ikp=7
!
      if(rkp.lt.5.66) ikp=6
      if(rkp.lt.4.66) ikp=5
      if(rkp.lt.3.66) ikp=4
      if(rkp.lt.2.66) ikp=3
      if(rkp.lt.1.66) ikp=2
      if(rkp.lt.0.66) ikp=1
!
   30 continue
!
      return
      end
!
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine convkptsy2(ex,ikp,rkp)
!
!     *****************************************************************
!     calcule le ikp de Tsyganenko 87s, 87l, 89c a partir du rkp
!     *****************************************************************
!
      character*(*) ex
!
!
      rkp= -1.
      if(ikp.lt.0) return
!
! *** calcul de rkp pout T87s
!
      if(ex(1:4).ne.'t87s') go to 10
!
      if(ikp.eq.8) rkp=5.66
      if(ikp.eq.7) rkp=5.33
      if(ikp.eq.6) rkp=4.66
      if(ikp.eq.5) rkp=3.66
      if(ikp.eq.4) rkp=2.66
      if(ikp.eq.3) rkp=2.00
      if(ikp.eq.2) rkp=1.33
      if(ikp.eq.1) rkp=0.66
!
   10 continue
!
! *** calcul de rkp pout T87l
!
      if(ex(1:4).ne.'t87l') go to 20
!
      ikp=6
!
      if(ikp.eq.6) rkp=5.00
      if(ikp.eq.5) rkp=4.66
      if(ikp.eq.4) rkp=3.66
      if(ikp.eq.3) rkp=2.66
      if(ikp.eq.2) rkp=1.66
      if(ikp.eq.1) rkp=0.66
!
   20 continue
!
! *** calcul de rkp pout T89c
!
      if(ex(1:4).ne.'t89c') go to 30
!
      if(ikp.eq.7) rkp=6.00
      if(ikp.eq.6) rkp=5.66
      if(ikp.eq.5) rkp=4.66
      if(ikp.eq.4) rkp=3.66
      if(ikp.eq.3) rkp=2.66
      if(ikp.eq.2) rkp=1.66
      if(ikp.eq.1) rkp=0.66
!
   30 continue
!
      return
      end
!
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine convrkp(rkp,ckp)
!
!     *****************************************************************
!     converti le rkp en ckp, ex: 2.66 -> 3-
!     *****************************************************************
!
      character*(*) ckp
!
      if(rkp.lt.0.) then
                    ckp='??'
                    return
                    endif
      ckp='8+'
!
!
      if(rkp.lt.8.66) ckp='8+'
      if(rkp.lt.8.33) ckp='8 '
      if(rkp.lt.8.00) ckp='8-'
!
      if(rkp.lt.7.66) ckp='7+'
      if(rkp.lt.7.33) ckp='7 '
      if(rkp.lt.7.00) ckp='7-'
!
      if(rkp.lt.6.66) ckp='6+'
      if(rkp.lt.6.33) ckp='6 '
      if(rkp.lt.6.00) ckp='6-'
!
      if(rkp.lt.5.66) ckp='5+'
      if(rkp.lt.5.33) ckp='5 '
      if(rkp.lt.5.00) ckp='5-'
!
      if(rkp.lt.4.66) ckp='4+'
      if(rkp.lt.4.33) ckp='4 '
      if(rkp.lt.4.00) ckp='4-'
!
      if(rkp.lt.3.66) ckp='3+'
      if(rkp.lt.3.33) ckp='3 '
      if(rkp.lt.3.00) ckp='3-'
!
      if(rkp.lt.2.66) ckp='2+'
      if(rkp.lt.2.33) ckp='2 '
      if(rkp.lt.2.00) ckp='2-'
!
      if(rkp.lt.1.66) ckp='1+'
      if(rkp.lt.1.33) ckp='1 '
      if(rkp.lt.1.00) ckp='1-'
!
      if(rkp.lt.0.66) ckp='0+'
      if(rkp.lt.0.33) ckp='0 '
      if(rkp.lt.0.00) ckp='0-'
!
!
      return
      end
!
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine dish(ikpbid,parmod,ps,x,y,z,bx,by,bz)
!
!     *****************************************************************
!     subroutine de transition entre igrf et exshor
!     ikpbid: pourcentage de champ "short"
!     ikpbid=0   -> 0
!     ikpbid=100 -> exshor avec ikp=1
!     Utilise avec trace2, donne la transition entre IGRF et IGRF+SHORT
!     *****************************************************************
!
      real (kind=4), dimension(10) :: parmod

      external t87s
!
      ikp=1
!
      call t87s(ikp,parmod,ps,x,y,z,bx,by,bz)
!
      bx=bx*float(ikpbid)/100.
      by=by*float(ikpbid)/100.
      bz=bz*float(ikpbid)/100.
!
      return
      end
!
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      SUBROUTINE TRACE2 (XI,YI,ZI,DIR,RLIM,xlim,yzlim,R0,IOPT,PARMOD, &
           EXNAME,INNAME,XF,YF,ZF,XX,YY,ZZ,L)
!
!     *****************************************************************
!     version modifie de trace de tsyganenko pour utilisation
!     avec option dipole et igrf en plus de t87s, t87l ou t89c
!     Doit etre appelle depuis calinca ou calinpo
!     *****************************************************************
!
!  TRACES A FIELD LINE FROM AN ARBITRARY POINT OF SPACE TO THE EARTH'S
!  SURFACE OR TO A MODEL LIMITING BOUNDARY.
!
!  THE HIGHEST ORDER OF SPHERICAL HARMONICS IN THE MAIN FIELD EXPANSION USED
!  IN THE MAPPING IS CALCULATED AUTOMATICALLY. IF INNAME=IGRF_GSM, THEN AN IGRF MODEL
!  FIELD WILL BE USED, AND IF INNAME=DIP, A PURE DIPOLE FIELD WILL BE USED.

!  IN ANY CASE, BEFORE CALLING TRACE, ONE SHOULD INVOKE RECALC, TO CALCULATE CORRECT
!  VALUES OF THE IGRF COEFFICIENTS AND ALL QUANTITIES NEEDED FOR TRANSFORMATIONS
!  BETWEEN COORDINATE SYSTEMS INVOLVED IN THIS CALCULATIONS.
!
!  ALTERNATIVELY, THE SUBROUTINE RECALC CAN BE INVOKED WITH THE DESIRED VALUES OF
!  IYEAR AND IDAY (TO SPECIFY THE DIPOLE MOMENT), WHILE THE VALUES OF THE DIPOLE
!  TILT ANGLE PSI (IN RADIANS) AND ITS SINE (SPS) AND COSINE (CPS) CAN BE EXPLICITLY
!  SPECIFIED AND FORWARDED TO THE COMMON BLOCK GEOPACK1 (11th, 12th, AND 16th ELEMENTS, RESP.)
!
!------------- INPUT PARAMETERS:
!
!   XI,YI,ZI - GSM COORDS OF INITIAL POINT (IN EARTH RADII, 1 RE = 6371.2 km),
!
!   DIR - SIGN OF THE TRACING DIRECTION: IF DIR=1.0 THEN WE MOVE ANTIPARALLEL TO THE
!     FIELD VECTOR (E.G. FROM NORTHERN TO SOUTHERN CONJUGATE POINT),
!     AND IF DIR=-1.0 THEN THE TRACING GOES IN THE OPPOSITE DIRECTION.
!
!   R0 -  RADIUS OF A SPHERE (IN RE) FOR WHICH THE FIELD LINE ENDPOINT COORDINATES
!     XF,YF,ZF  SHOULD BE CALCULATED.
!
!   RLIM - UPPER LIMIT OF THE GEOCENTRIC DISTANCE, WHERE THE TRACING IS TERMINATED.
!
!   yzlim : replace the value of 400 in original code
!   xlim  : replace the value of  20 in original code
!
!   IOPT - A MODEL INDEX; CAN BE USED FOR SPECIFYING AN OPTION OF THE EXTERNAL FIELD
!       MODEL (E.G., INTERVAL OF THE KP-INDEX). ALTERNATIVELY, ONE CAN USE THE ARRAY
!       PARMOD FOR THAT PURPOSE (SEE BELOW); IN THAT CASE IOPT IS JUST A DUMMY PARAMETER.
!
!   PARMOD -  A 10-ELEMENT ARRAY CONTAINING MODEL PARAMETERS, NEEDED FOR A UNIQUE
!      SPECIFICATION OF THE EXTERNAL FIELD. THE CONCRETE MEANING OF THE COMPONENTS
!      OF PARMOD DEPENDS ON A SPECIFIC VERSION OF THE EXTERNAL FIELD MODEL.
!
!   EXNAME - NAME OF A SUBROUTINE PROVIDING COMPONENTS OF THE EXTERNAL MAGNETIC FIELD
!    (E.G., T96_01).
!   INNAME - NAME OF A SUBROUTINE PROVIDING COMPONENTS OF THE INTERNAL MAGNETIC FIELD
!    (EITHER DIP OR IGRF_GSM).
!
!-------------- OUTPUT PARAMETERS:
!
!   XF,YF,ZF - GSM COORDS OF THE LAST CALCULATED POINT OF A FIELD LINE
!   XX,YY,ZZ - ARRAYS, CONTAINING COORDS OF FIELD LINE POINTS. HERE THEIR MAXIMAL LENGTH WAS
!      ASSUMED EQUAL TO 999.
!   L - ACTUAL NUMBER OF THE CALCULATED FIELD LINE POINTS. IF L EXCEEDS 999, TRACING
!     TERMINATES, AND A WARNING IS DISPLAYED.
!
!
!     LAST MODIFICATION:  MARCH 31, 2003.
!
!     AUTHOR:  N. A. TSYGANENKO
!
      DIMENSION XX(1000),YY(1000),ZZ(1000), PARMOD(10)

! jah, foresys :COMMON array size mismatch for AA in /geopack1/ from dip (is [26] expected [10])
! jah, COMMON /GEOPACK1/ AA(26),DD,BB(8)

      COMMON /GEOPACK1/ BID1,BID2,BID3,BID4,BID5,BID6,BID7,BID8,BID9, &
           BID10,BID11,BID12,BID13,BID14,BID15,BID16,BID17,BID18,BID19, &
           BID20,BID21,BID22,BID23,BID24,BID25,BID26,DD,BID28,BID29, &
           BID30(6)
      EXTERNAL EXNAME,INNAME
!
      ERR=0.0001
      L=0
      DS=0.5*DIR
      X=XI
      Y=YI
      Z=ZI
      DD=DIR
! jah, foresys : Result of assignment to AL is not used
! jah,       AL=0.

! jah, foresys : Variable XR,YR,ZR may not have been  initialized
      XR=0.
      YR=0.
      ZR=0.

!
!  here we call RHAND just to find out the sign of the radial component of the field
!   vector, and to determine the initial direction of the tracing (i.e., either away
!   or towards Earth):
!
      CALL RHAND (X,Y,Z,R1,R2,R3,IOPT,PARMOD,EXNAME,INNAME)
      AD=0.01
      IF (X*R1+Y*R2+Z*R3.LT.0.) AD=-0.01
!
!     |AD|=0.01 and its sign follows the rule:
! (1) if DIR=1 (tracing antiparallel to B vector) then the sign of AD is the same as of Br
! (2) if DIR=-1 (tracing parallel to B vector) then the sign of AD is opposite to that of Br
!     AD is defined in order to initialize the value of RR (radial distance at previous step):

      RR=SQRT(X**2+Y**2+Z**2)+AD
  1   L=L+1
      IF(L.GT.999) GOTO 7
      XX(L)=X
      YY(L)=Y
      ZZ(L)=Z
      RYZ=Y**2+Z**2
      R2=X**2+RYZ
      R=SQRT(R2)

!  check if the line hit the outer tracing boundary; if yes, then terminate
!   the tracing (label 8):

      yzlim2=yzlim**2
      IF (R.GT.RLIM.OR.RYZ.GT.yzlim2.OR.X.GT.xlim) GOTO 8
!
!  check whether or not the inner tracing boundary was crossed from outside,
!  if yes, then calculate the footpoint position by interpolation (go to label 6):
!
      IF (R.LT.R0.AND.RR.GT.R) GOTO 6

!  check if (i) we are moving outward, or (ii) we are still sufficiently
!    far from Earth (beyond R=5Re); if yes, proceed further:
!
      IF (R.GE.RR.OR.R.GT.5.) GOTO 5

!  now we moved closer inward (between R=3 and R=5); go to 3 and begin logging
!  previous values of X,Y,Z, to be used in the interpolation (after having
!  crossed the inner tracing boundary):

      IF (R.GE.3.) GOTO 3
!
!  we entered inside the sphere R=3: to avoid too large steps (and hence inaccurate
!  interpolated position of the footpoint), enforce the progressively smaller
!  stepsize values as we approach the inner boundary R=R0:
!
      FC=0.2
      IF(R-R0.LT.0.05) FC=0.05
      AL=FC*(R-R0+0.2)
      DS=DIR*AL
      GOTO 4
  3   DS=DIR
  4   XR=X
      YR=Y
      ZR=Z
  5   RR=R
      CALL STEP (X,Y,Z,DS,ERR,IOPT,PARMOD,EXNAME,INNAME)
      GOTO 1
!
!  find the footpoint position by interpolating between the current and previous
!   field line points:
!
  6   R1=(R0-R)/(RR-R)
      X=X-(X-XR)*R1
      Y=Y-(Y-YR)*R1
      Z=Z-(Z-ZR)*R1
      GOTO 8
  7   WRITE (*,10)
      L=999
  8   XF=X
      YF=Y
      ZF=Z
      RETURN
 10   FORMAT(//,1X,'**** COMPUTATIONS IN THE SUBROUTINE TRACE ARE', &
      ' TERMINATED: THE CURRENT NUMBER OF POINTS EXCEEDED 1000 ****'//)
      END
!
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine check_ex(ex)
!
      character*(*) ex
!
!     *****************************************************************
!     check validity of ex name
!     *****************************************************************
!
      if(ex(1:4).ne.'t89c'  .and. &
         ex(1:4).ne.'t87s'  .and. &
         ex(1:4).ne.'t87l'  .and. &
         ex(1:6).ne.'t96_01'.and. &
         ex(1:6).ne.'t01_01'.and. &
         ex(1:5).ne.'t04_s' .and. &
         ex(1:4).ne.'dish'  .and. &
         ex(1:4).ne.'igrf'  .and. &
         ex(1:4).ne.'dipo') then
         print*, 'ex=',ex, ' not allowed'
         
      print*, 'only t89c, t87s, t87l, t96_01, t01_01, t04_s,'
      print*, '     dish, igrf, dipo'
      stop 'COKP           *** ABORTED ! Model not allowed          ***'
      endif
!
      return
      end
!
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine check_ikp(ex,ikp)
!
      character*(*) ex
!
!     *****************************************************************
!     check validity of ikp value
!     *****************************************************************
!
      if(ex(1:4).eq.'t89c' ) then
                       if(ikp.lt.1.or.ikp.gt.7) then
                                                print*,'ikp: 1-7 only'
      stop 'COKP           *** ABORTED ! ikp must be >1 and <7      ***'
                                                else
                                                return
                                                endif
                       endif
!
      if(ex(1:4).eq.'t87s') then
                       if(ikp.lt.1.or.ikp.gt.8) then
                                                print*,'ikp: 1-8 only'
      stop 'COKP           *** ABORTED ! ikp must be >1 and <8      ***'
                                                else
                                                return
                                                endif
                       endif
!
      if(ex(1:4).eq.'t87l') then
                       if(ikp.lt.1.or.ikp.gt.6) then
                                                print*,'ikp: 1-6 only'
      stop 'COKP           *** ABORTED ! ikp must be >1 and <6      ***'
                                                else
                                                return
                                                endif
                        endif
!
! *** pour exdipo ou exigrf, ikp sans objets
!
      return
      end
!
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine cotsymag(ex,ikp,parmod,xgsm,ygsm,zgsm,bx,by,bz)
!
      character*(*) ex
      real (kind=4), dimension(10) :: parmod

      external t87s,t87l,t89c,t96_01,t01_01,t04_s
      external dish,exzero,IGRF_GSM,DIP
!
! ---------------------------------------------------------------------
!
!     compute_Tsyganenko_magnetic_field model
!
!     ex= only t89c, t87s, t87l, dish, igrf, dipo
!     ikp =1-7      for t89c
!          1-8      for t87s
!          1-6      for t87l
!          not used for igrf or dipole or t96_01 or t01_01 or t04_s
!
!     x,y,z   = point of space in GSM system and RE unit
!     by,by,bz= output field, in nT
!
!     recalc must be called before (init time parameters)
!
! ---------------------------------------------------------------------
!
!
!     verification des parametres en entres, ikp quelconque pour t96_01
!     t01_01 t04_s
!
      call check_ex(ex)
      call check_ikp(ex,ikp)
!
! *** test sur distance
!
      x=xgsm
      y=ygsm
      z=zgsm
!
      rr=sqrt(x**2 + y**2 + z**2)
      if(rr.lt.1.) then
                   print*, 'R=',rr,' taken to 1. Re in cotsymag'
                   call geogsm(0.,0.,1.001,x,y,z, 1)
                   endif
!
! *** calcul du dipole tilt
!
      call gdiptil(psi)
!
! *** calcul du champ
!
      if(ex(1:4).eq.'t89c' ) then
                            call IGRF_GSM(X,Y,Z,hx,hy,hz)
                            call t89c(ikp,parmod,psi,x,y,z,bx,by,bz)
                            bx=bx+hx
                            by=by+hy
                            bz=bz+hz
                            return
                            endif
!
      if(ex(1:4).eq.'t87s') then
                            call IGRF_GSM(X,Y,Z,hx,hy,hz)
                            call t87s(ikp,parmod,psi,x,y,z,bx,by,bz)
                            bx=bx+hx
                            by=by+hy
                            bz=bz+hz
                            return
                            endif
!
      if(ex(1:4).eq.'t87l') then
                            call IGRF_GSM(X,Y,Z,hx,hy,hz)
                            call t87l(ikp,parmod,psi,x,y,z,bx,by,bz)
                            bx=bx+hx
                            by=by+hy
                            bz=bz+hz
                            return
                            endif
!
      if(ex(1:6).eq.'t96_01') then
                            call IGRF_GSM(X,Y,Z,hx,hy,hz)
                            call t96_01(ikp,parmod,psi,x,y,z,bx,by,bz)
                            bx=bx+hx
                            by=by+hy
                            bz=bz+hz
                            return
                            endif
!
      if(ex(1:6).eq.'t01_01') then
                            call IGRF_GSM(X,Y,Z,hx,hy,hz)
                            call t01_01(ikp,parmod,psi,x,y,z,bx,by,bz)
                            bx=bx+hx
                            by=by+hy
                            bz=bz+hz
                            return
                            endif
!
      if(ex(1:5).eq.'t04_s') then
                            call IGRF_GSM(X,Y,Z,hx,hy,hz)
                            call t04_s(ikp,parmod,psi,x,y,z,bx,by,bz)
                            bx=bx+hx
                            by=by+hy
                            bz=bz+hz
                            return
                            endif
!
      if(ex(1:4).eq.'dish') then
                            call IGRF_GSM(X,Y,Z,hx,hy,hz)
                            call dish(ikp,parmod,psi,x,y,z,bx,by,bz)
                            bx=bx+hx
                            by=by+hy
                            bz=bz+hz
                            return
                            endif
!
      if(ex(1:4).eq.'igrf') then
                            call IGRF_GSM(X,Y,Z,bx,by,bz)
                            return
                            endif
!
      if(ex(1:4).eq.'dipo') then
                            call DIP(X,Y,Z,bx,by,bz)
                            return
                            endif
!
      print*
      print*, '*** costymag: ex=',ex
      print*, '***           only t89c, t87s, t87l, dish, igrf, dipo'
      stop 'COKP           *** ABORTED ! MODEL UNKNOWN              ***'
!
      end
!
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
       SUBROUTINE GDIPTIL(PS)
!
!-----INPUT PARAMETERS:  NONE
!-----OUTPUT PARAMETERS: PS - GEODIPOLE TILT ANGLE IN RADIANS
!
!     ATTENTION:  SUBROUTINE  RECALC  MUST BE INVOKED BEFORE GDIPTIL IN
!     TWO CASES:
!     /A/  BEFORE THE FIRST INVOCATION OF GDIPTIL
!     /B/  IF THE VALUES OF IYEAR,IDAY,IHOUR,MIN,ISEC HAVE BEEN CHANGED
!
!     P. Robert, June 2002.
!
!                      ----------------------
!
       COMMON /GEOPACK1/ ST0,CT0,SL0,CL0,CTCL,STCL,CTSL,STSL,SFI,CFI,  &
           SPS,CPS,SHI,CHI,HI,PSI,XMUT,A11,A21,A31,A12,A22,A32,A13,A23,&
            A33,DS3,CGST,SGST,BA(6)
!
      PS=PSI
!
      RETURN
      END
!
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine rap_exikp(ex,ikp)
!
      character*(*) ex
!
!     *****************************************************************
!
!     read and print ex et ikp for tsyganenko model
!
!     *****************************************************************
!
!
      print*, 'magnetic field model ? '
      print*, '    (only t89c, t87s, t87l, t96_01, t01_01, t04_s,'
      print*, '          dish, igrf, dipo)'
      read 100,  ex
      print * ,  ex
!
      call check_ex(ex)
!
      print*, 'ikp ? (ex:4)'
      read *, ikp
      print*, ikp
!
      call check_ikp(ex,ikp)
!
  100 format(a)
!
      return
      end
!
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine rap_parmod(parmod, idim)

      real (kind=4), dimension(idim) :: parmod

      do i=1,idim
!        init
         parmod(i)=0.0
!        print input variable
         write(*,*,advance='no') 'tsyganenko, parmod ',i,' ?'
!        read value
         read *, parmod(i)
!        print value of variable
         write(*,'(19x,i4)') parmod(i)
      enddo

      return
      end

!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX

      subroutine rap_unit(sc,unit)
!
! *** output: type (c or s) and unit (km or RE)
!
      character*(*) unit
      character*1 sc
!
      print*, 'input coordinates: spherical(s)  cartesian (c)'
      read(*,100) sc
      print*, sc
      if(sc.ne.'s'.and.sc.ne.'c') &
      stop 'COMAGLIB       *** ABORTED ! Coordinates must be s or c ***'
!
      print*, 'unit for distance ? (RE or km, 1 RE=6371.2 km)'
      read(*,100) unit
      print*, unit
      if(unit(1:2).ne.'RE'.and.unit(1:2).ne.'km') then
                    print*, 'unit must be only RE or km'
      stop 'COMAGLIB       *** ABORTED ! Units must be km or RE     ***'
                    endif
!
  100 format(a)
      return
      end
!
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine rap_reper(rep)
!
      character*(*) rep
!
      print*, 'Coordinate system ? gsm, gse, sma, mag,geo,gei'
      read *,rep
      print*,rep
!
      if(rep.eq.'gsm') return
      if(rep.eq.'gse') return
      if(rep.eq.'sma') return
      if(rep.eq.'mag') return
      if(rep.eq.'geo') return
      if(rep.eq.'gei') return
!
      print*, 'rep not allowed for this routine'
      stop 'COMAGLIB       *** ABORTED ! Coordinate system not valid***'
      end
!
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine rap_coord(sc,unit,x,y,z,ieof)
!
      character*(*) sc,unit
!
! *** output: cartesian in unit
!
      pisd=3.1415927/180.
      ieof=0
!
      if(unit(1:2).ne.'RE'.and.unit(1:2).ne.'km') then
                    print*, 'unit must be only RE or km'
      stop 'COMAGLIB       *** ABORTED ! Units must be km or RE     ***'
                    endif
!
      if(sc.ne.'s'.and.sc.ne.'c') &
      stop 'COMAGLIB       *** ABORTED ! Coordinates must be s or c ***'
!
! *** lecture des xyz suivant sc et unit
!
      if(sc.eq.'s') then
                    print*
                    print*, 'r,tet, phi (',unit(1:2),',deg.) ?'
                    read(5, *,end=20)  r,tet, phi
                    print*,  r,tet, phi, '  (deg.)'
                    tet=tet*pisd
                    phi=phi*pisd
                    call tsphcar(r,tet,phi,x,y,z)
                    print*, 'x,y,z (',unit(1:2),') '
                    print*,  x,y,z
!
                    else
                    print*
                    print*, 'x,y,z ? (',unit(1:2),') '
                    read(5, *,end=20)  x,y,z
                    print*,  x,y,z
                    call tcarsph(x,y,z,r,tet,phi)
                    print*, 'r,tet, phi  '
                    tet=tet/pisd
                    phi=phi/pisd
                    print*, r,tet, phi, '  (',unit(1:2),',deg.)'
                    endif
      return
!
   20 continue
      ieof=1
      end
!
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine strlen(str,nc)
!
!     ------------------------------------------------------------------
! *   Objet  : calcul de la longueur utile  d'une chaine de caracteres
! *   Classe : depouillement specifique Roproc
! *   Auteur : P. Robert, CETP, 1998
!     ------------------------------------------------------------------
!
      character*(*) str
!
!                    *********************
!
      nbc=len(str)
!
! *** nc= longueur au de la de laquelle il n'y a que des blancs
!         (mais la chaine 1:nc peu contenir des blancs)
!
      do 10 i=1,nbc
      ii=nbc-i+1
      if(str(ii:ii).ne.' ') go to 12
   10 continue
      ii=0
   12 continue
      nc=ii
!
      return
      end
!
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
!      subroutine longut(mot,nc)
!
!      character mot*(*)
!
!     *****************************************************************
!     calcul de la longueur utile
!     *****************************************************************
!
!      nbc=len(mot)
!
!      do 10 i=1,nbc
!      ii=nbc-i+1
!      if(mot(ii:ii).ne." ") go to 12
!   10 continue
!      ii=0
!   12 continue
!      nc=ii
!
!      return
!      end
!
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
