
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
!      SUBROUTINE TRACE_PR(XI,YI,ZI,DIR,RLIM,R0,XX,YY,ZZ,L)
! RLIM fout la merde...
      SUBROUTINE TRACE_PR(XI,YI,ZI,DIR,R0,XX,YY,ZZ,L)
!
!     *****************************************************************
!     version modifie de trace de tsyganenko
!     *****************************************************************


      DIMENSION XX(1000),YY(1000),ZZ(1000)
     
!
      ERR=0.0001
      L=0
      DS=0.5*DIR
      X=XI
      Y=YI
      Z=ZI
      DD=DIR

      XR=0.
      YR=0.
      ZR=0.

      yzlim=20.
      xlim=20.

!
!  here we call RHAND_PRjust 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_PR(X,Y,Z,R1,R2,R3,iout)
      if(iout == 1) go to 11
      ! print*,'ro1'
      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):

!      IF (R.GT.RLIM.OR.RYZ.GT.yzlim.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_PR(X,Y,Z,DS,ERR,iout)
      if(iout == 1) go to 11
      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   continue
!     WRITE (*,10) XI,YI,ZI,X,Y,Z
      L=999
  8   XF=X
      YF=Y
      ZF=Z

!      write(*,'(a,3f9.1,a,3f9.1a,i4)') 'Start=',XI,YI,ZI,'   end=',XF,YF,ZF,'   N=',L
      RETURN
      
!10   FORMAT('start:',3f9.1,'    end:',3f9.1,'   more than 1000 pts')
 11   continue

 !     ! print*, 'iout=1'
 ! print*,'or'
      return
      END
!
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!=====================================================================================
!
      SUBROUTINE RHAND_PR(X,Y,Z,R1,R2,R3,iout)

      COMMON /GEOPACK1/ BID1,BID2,BID3,BID4,BID5,BID6,BID7,BID8,BID9, &
           BID10,BID11,BID12,BID13,BID14,BID15,PSI,BID17,BID18,BID19, &
           BID20,BID21,BID22,BID23,BID24,BID25,BID26,DS3,BID28, &
           BID29,BID30(6)
           
           ! print*,'rr'

 !     BX=BXGSM+HXGSM 
 !     BY=BYGSM+HYGSM
 !     BZ=BZGSM+HZGSM
  iout=0
     ! ------------------------------------
      call give_model_B(X,Y,Z,BX,BY,BZ,iout)
     ! ------------------------------------
     

      
  ! print*,'rhand: BX,BY,BZ=',BX,BY,BZ
  if(iout == 1) return

      B=DS3/SQRT(BX**2+BY**2+BZ**2)
      R1=BX*B
      R2=BY*B
      R3=BZ*B
      ! print*,'rgm'     
      RETURN
      END
!
!===================================================================================
!
      SUBROUTINE STEP_PR(X,Y,Z,DS,ERRIN,iout)
!
!   RE-CALCULATES {X,Y,Z}, MAKING A STEP ALONG A FIELD LINE.
!   DS IS THE STEP SIZE, ERRIN IS PERMISSIBLE ERROR VALUE, IOPT SPECIFIES THE EXTERNAL
!   MODEL VERSION, THE ARRAY PARMOD CONTAINS INPUT PARAMETERS FOR THAT MODEL
!   EXNAME IS THE NAME OF THE EXTERNAL FIELD SUBROUTINE
!   INNAME IS THE NAME OF THE INTERNAL FIELD SUBROUTINE (EITHER DIP OR IGRF)
!
!   ALL THE PARAMETERS ARE INPUT ONES; OUTPUT IS THE RENEWED TRIPLET X,Y,Z
!
!     LAST MODIFICATION:  MARCH 31, 2003
!
!     AUTHOR:  N. A. TSYGANENKO
!

      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,DS3,BID28,BID29, &
           BID30(6)
      EXTERNAL EXNAME,INNAME
      
      ! print*,'ss'

  1   DS3=-DS/3.
      CALL RHAND_PR(X,Y,Z,R11,R12,R13,iout)
      if(iout == 1) return
      CALL RHAND_PR(X+R11,Y+R12,Z+R13,R21,R22,R23,iout)
      if(iout == 1) return
      CALL RHAND_PR(X+.5*(R11+R21),Y+.5*(R12+R22),Z+.5*(R13+R23),R31,R32,R33,iout)
      if(iout == 1) return
      CALL RHAND_PR(X+.375*(R11+3.*R31),Y+.375*(R12+3.*R32),Z+.375*(R13+3.*R33),R41,R42,R43,iout)
      if(iout == 1) return
      CALL RHAND_PR(X+1.5*(R11-3.*R31+4.*R41),Y+1.5*(R12- &
      3.*R32+4.*R42),Z+1.5*(R13-3.*R33+4.*R43),R51,R52,R53,iout)
      if(iout == 1) return
      
      ERRCUR=ABS(R11-4.5*R31+4.*R41-.5*R51)+ABS(R12-4.5*R32+4.*R42-.5* &
      R52)+ABS(R13-4.5*R33+4.*R43-.5*R53)
      IF (ERRCUR.LT.ERRIN) GOTO 2
      DS=DS*.5
      GOTO 1
  2   X=X+.5*(R11+4.*R41+R51)
      Y=Y+.5*(R12+4.*R42+R52)
      Z=Z+.5*(R13+4.*R43+R53)
      IF(ERRCUR.LT.ERRIN*.04.AND.ABS(DS).LT.1.33) DS=DS*1.5
      RETURN
      END
!
!==============================================================================
!!
      SUBROUTINE RECALC (IYEAR,IDAY,IHOUR,MIN,ISEC)
!
!  1. PREPARES ELEMENTS OF ROTATION MATRICES FOR TRANSFORMATIONS OF VECTORS BETWEEN
!     SEVERAL COORDINATE SYSTEMS, MOST FREQUENTLY USED IN SPACE PHYSICS.
!
!  2. PREPARES COEFFICIENTS USED IN THE CALCULATION OF THE MAIN GEOMAGNETIC FIELD
!      (IGRF MODEL)
!
!  THIS SUBROUTINE SHOULD BE INVOKED BEFORE USING THE FOLLOWING SUBROUTINES:
!    IGRF_GEO, IGRF_GSM, DIP, GEOMAG, GEOGSM, MAGSM, SMGSM, GSMGSE, GEIGEO.
!
!  THERE IS NO NEED TO REPEATEDLY INVOKE RECALC, IF MULTIPLE CALCULATIONS ARE MADE
!    FOR THE SAME DATE AND TIME.
!
!-----INPUT PARAMETERS:
!
!     IYEAR   -  YEAR NUMBER (FOUR DIGITS)
!     IDAY  -  DAY OF YEAR (DAY 1 = JAN 1)
!     IHOUR -  HOUR OF DAY (00 TO 23)
!     MIN   -  MINUTE OF HOUR (00 TO 59)
!     ISEC  -  SECONDS OF MINUTE (00 TO 59)
!
!-----OUTPUT PARAMETERS:   NONE (ALL OUTPUT QUANTITIES ARE PLACED
!                         INTO THE COMMON BLOCKS /GEOPACK1/ AND /GEOPACK2/)
!
!    OTHER SUBROUTINES CALLED BY THIS ONE: SUN
!
!    AUTHOR:  N.A. TSYGANENKO
!    DATE:    DEC.1, 1991
!
!    LAST REVISION: APRIL 3, 2003
!      The code now includes preparation of the model coefficients for the subroutines
!      IGRF and GEOMAG. This eliminates the need for the SAVE statements, used in the
!      old versions, making the codes easier and more compiler-independent.
!
      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)
!
!  THE COMMON BLOCK /GEOPACK1/ CONTAINS ELEMENTS OF THE ROTATION MATRICES AND OTHER
!   PARAMETERS RELATED TO THE COORDINATE TRANSFORMATIONS PERFORMED BY THIS PACKAGE
!
      COMMON /GEOPACK2/ G(66),H(66),REC(66)
!
!  THE COMMON BLOCK /GEOPACK2/ CONTAINS COEFFICIENTS OF THE IGRF FIELD MODEL, CALCULATED
!    FOR A GIVEN YEAR AND DAY FROM THEIR STANDARD EPOCH VALUES. THE ARRAY REC CONTAINS
!    COEFFICIENTS USED IN THE RECURSION RELATIONS FOR LEGENDRE ASSOCIATE POLYNOMIALS.
!
      DIMENSION G65(66),H65(66),G70(66),H70(66),G75(66),H75(66),G80(66), &
       H80(66),G85(66),H85(66),G90(66),H90(66),G95(66),H95(66),G00(66), &
       H00(66),DG00(45),DH00(45)

      save G65,H65,G70,H70,G75,H75,G80,  &
           H80,G85,H85,G90,H90,G95,H95,G00, H00,DG00,DH00
!
      DATA G65/0.,-30334.,-2119.,-1662.,2997.,1594.,1297.,-2038.,1292., &
      856.,957.,804.,479.,-390.,252.,-219.,358.,254.,-31.,-157.,-62., &
      45.,61.,8.,-228.,4.,1.,-111.,75.,-57.,4.,13.,-26.,-6.,13.,1.,13., &
      5.,-4.,-14.,0.,8.,-1.,11.,4.,8.,10.,2.,-13.,10.,-1.,-1.,5.,1.,-2., &
      -2.,-3.,2.,-5.,-2.,4.,4.,0.,2.,2.,0./
      DATA H65/0.,0.,5776.,0.,-2016.,114.,0.,-404.,240.,-165.,0.,148., &
      -269.,13.,-269.,0.,19.,128.,-126.,-97.,81.,0.,-11.,100.,68.,-32., &
      -8.,-7.,0.,-61.,-27.,-2.,6.,26.,-23.,-12.,0.,7.,-12.,9.,-16.,4., &
      24.,-3.,-17.,0.,-22.,15.,7.,-4.,-5.,10.,10.,-4.,1.,0.,2.,1.,2., &
      6.,-4.,0.,-2.,3.,0.,-6./
!
      DATA G70/0.,-30220.,-2068.,-1781.,3000.,1611.,1287.,-2091.,1278., &
      838.,952.,800.,461.,-395.,234.,-216.,359.,262.,-42.,-160.,-56., &
      43.,64.,15.,-212.,2.,3.,-112.,72.,-57.,1.,14.,-22.,-2.,13.,-2., &
      14.,6.,-2.,-13.,-3.,5.,0.,11.,3.,8.,10.,2.,-12.,10.,-1.,0.,3., &
      1.,-1.,-3.,-3.,2.,-5.,-1.,6.,4.,1.,0.,3.,-1./
      DATA H70/0.,0.,5737.,0.,-2047.,25.,0.,-366.,251.,-196.,0.,167., &
      -266.,26.,-279.,0.,26.,139.,-139.,-91.,83.,0.,-12.,100.,72.,-37., &
      -6.,1.,0.,-70.,-27.,-4.,8.,23.,-23.,-11.,0.,7.,-15.,6.,-17.,6., &
      21.,-6.,-16.,0.,-21.,16.,6.,-4.,-5.,10.,11.,-2.,1.,0.,1.,1.,3., &
      4.,-4.,0.,-1.,3.,1.,-4./
!
      DATA G75/0.,-30100.,-2013.,-1902.,3010.,1632.,1276.,-2144.,1260., &
      830.,946.,791.,438.,-405.,216.,-218.,356.,264.,-59.,-159.,-49., &
      45.,66.,28.,-198.,1.,6.,-111.,71.,-56.,1.,16.,-14.,0.,12.,-5., &
      14.,6.,-1.,-12.,-8.,4.,0.,10.,1.,7.,10.,2.,-12.,10.,-1.,-1.,4., &
      1.,-2.,-3.,-3.,2.,-5.,-2.,5.,4.,1.,0.,3.,-1./
      DATA H75/0.,0.,5675.,0.,-2067.,-68.,0.,-333.,262.,-223.,0.,191., &
      -265.,39.,-288.,0.,31.,148.,-152.,-83.,88.,0.,-13.,99.,75.,-41., &
      -4.,11.,0.,-77.,-26.,-5.,10.,22.,-23.,-12.,0.,6.,-16.,4.,-19.,6., &
      18.,-10.,-17.,0.,-21.,16.,7.,-4.,-5.,10.,11.,-3.,1.,0.,1.,1.,3., &
      4.,-4.,-1.,-1.,3.,1.,-5./
!
      DATA G80/0.,-29992.,-1956.,-1997.,3027.,1663.,1281.,-2180.,1251., &
      833.,938.,782.,398.,-419.,199.,-218.,357.,261.,-74.,-162.,-48., &
      48.,66.,42.,-192.,4.,14.,-108.,72.,-59.,2.,21.,-12.,1.,11.,-2., &
      18.,6.,0.,-11.,-7.,4.,3.,6.,-1.,5.,10.,1.,-12.,9.,-3.,-1.,7.,2., &
      -5.,-4.,-4.,2.,-5.,-2.,5.,3.,1.,2.,3.,0./
      DATA H80/0.,0.,5604.,0.,-2129.,-200.,0.,-336.,271.,-252.,0.,212., &
      -257.,53.,-297.,0.,46.,150.,-151.,-78.,92.,0.,-15.,93.,71.,-43., &
      -2.,17.,0.,-82.,-27.,-5.,16.,18.,-23.,-10.,0.,7.,-18.,4.,-22.,9., &
      16.,-13.,-15.,0.,-21.,16.,9.,-5.,-6.,9.,10.,-6.,2.,0.,1.,0.,3., &
      6.,-4.,0.,-1.,4.,0.,-6./
!
      DATA G85/0.,-29873.,-1905.,-2072.,3044.,1687.,1296.,-2208.,1247., &
      829.,936.,780.,361.,-424.,170.,-214.,355.,253.,-93.,-164.,-46., &
      53.,65.,51.,-185.,4.,16.,-102.,74.,-62.,3.,24.,-6.,4.,10.,0.,21., &
      6.,0.,-11.,-9.,4.,4.,4.,-4.,5.,10.,1.,-12.,9.,-3.,-1.,7.,1.,-5., &
      -4.,-4.,3.,-5.,-2.,5.,3.,1.,2.,3.,0./
      DATA H85/0.,0.,5500.,0.,-2197.,-306.,0.,-310.,284.,-297.,0.,232., &
      -249.,69.,-297.,0.,47.,150.,-154.,-75.,95.,0.,-16.,88.,69.,-48., &
      -1.,21.,0.,-83.,-27.,-2.,20.,17.,-23.,-7.,0.,8.,-19.,5.,-23.,11., &
      14.,-15.,-11.,0.,-21.,15.,9.,-6.,-6.,9.,9.,-7.,2.,0.,1.,0.,3., &
      6.,-4.,0.,-1.,4.,0.,-6./
!
      DATA G90/0., -29775.,  -1848.,  -2131.,   3059.,   1686.,   1314., &
           -2239.,   1248.,    802.,    939.,    780.,    325.,   -423., &
             141.,   -214.,    353.,    245.,   -109.,   -165.,    -36., &
              61.,     65.,     59.,   -178.,      3.,     18.,    -96., &
              77.,    -64.,      2.,     26.,     -1.,      5.,      9., &
               0.,     23.,      5.,     -1.,    -10.,    -12.,      3., &
               4.,      2.,     -6.,      4.,      9.,      1.,    -12., &
               9.,     -4.,     -2.,      7.,      1.,     -6.,     -3., &
              -4.,      2.,     -5.,     -2.,      4.,      3.,      1., &
               3.,      3.,      0./

      DATA H90/0.,      0.,   5406.,      0.,  -2279.,   -373.,      0., &
            -284.,    293.,   -352.,      0.,    247.,   -240.,     84., &
            -299.,      0.,     46.,    154.,   -153.,    -69.,     97., &
               0.,    -16.,     82.,     69.,    -52.,      1.,     24., &
               0.,    -80.,    -26.,      0.,     21.,     17.,    -23., &
              -4.,      0.,     10.,    -19.,      6.,    -22.,     12., &
              12.,    -16.,    -10.,      0.,    -20.,     15.,     11., &
              -7.,     -7.,      9.,      8.,     -7.,      2.,      0., &
               2.,      1.,      3.,      6.,     -4.,      0.,     -2., &
               3.,     -1.,     -6./

      DATA G95/0., -29682.,  -1789.,  -2197.,   3074.,   1685.,   1329., &
           -2268.,   1249.,    769.,    941.,    782.,    291.,   -421., &
             116.,   -210.,    352.,    237.,   -122.,   -167.,    -26., &
              66.,     64.,     65.,   -172.,      2.,     17.,    -94., &
              78.,    -67.,      1.,     29.,      4.,      8.,     10., &
              -2.,     24.,      4.,     -1.,     -9.,    -14.,      4., &
               5.,      0.,     -7.,      4.,      9.,      1.,    -12., &
               9.,     -4.,     -2.,      7.,      0.,     -6.,     -3., &
              -4.,      2.,     -5.,     -2.,      4.,      3.,      1., &
               3.,      3.,      0./

      DATA H95/0.,      0.,   5318.,      0.,  -2356.,   -425.,      0., &
            -263.,    302.,   -406.,      0.,    262.,   -232.,     98., &
            -301.,      0.,     44.,    157.,   -152.,    -64.,     99., &
               0.,    -16.,     77.,     67.,    -57.,      4.,     28., &
               0.,    -77.,    -25.,      3.,     22.,     16.,    -23., &
              -3.,      0.,     12.,    -20.,      7.,    -21.,     12., &
              10.,    -17.,    -10.,      0.,    -19.,     15.,     11., &
              -7.,     -7.,      9.,      7.,     -8.,      1.,      0., &
               2.,      1.,      3.,      6.,     -4.,      0.,     -2., &
               3.,     -1.,     -6./

      DATA G00/0., -29615.,  -1728.,  -2267.,   3072.,   1672.,   1341., &
           -2290.,   1253.,    715.,    935.,    787.,    251.,   -405., &
             110.,   -217.,    351.,    222.,   -131.,   -169.,    -12., &
              72.,     68.,     74.,   -161.,     -5.,     17.,    -91., &
              79.,    -74.,      0.,     33.,      9.,      7.,      8., &
              -2.,     25.,      6.,     -9.,     -8.,    -17.,      9., &
               7.,     -8.,     -7.,      5.,      9.,      3.,    - 8., &
               6.,     -9.,     -2.,      9.,     -4.,     -8.,     -2., &
              -6.,      2.,     -3.,      0.,      4.,      1.,      2., &
               4.,      0.,     -1./


      DATA H00/0.,      0.,   5186.,      0.,  -2478.,   -458.,      0., &
            -227.,    296.,   -492.,      0.,    272.,   -232.,    119., &
            -304.,      0.,     44.,    172.,   -134.,    -40.,    107., &
               0.,    -17.,     64.,     65.,    -61.,      1.,     44., &
               0.,    -65.,    -24.,      6.,     24.,     15.,    -25., &
              -6.,      0.,     12.,    -22.,      8.,    -21.,     15., &
               9.,    -16.,     -3.,      0.,    -20.,     13.,     12., &
              -6.,     -8.,      9.,      4.,     -8.,      5.,      0., &
               1.,      0.,      4.,      5.,     -6.,     -1.,     -3., &
               0.,     -2.,     -8./


      DATA DG00/0.0,  14.6,    10.7,   -12.4,     1.1,    -1.1,     0.7, &
               -5.4,   0.9,    -7.7,    -1.3,     1.6,    -7.3,     2.9, &
               -3.2,   0.0,    -0.7,    -2.1,    -2.8,    -0.8,     2.5, &
                1.0,  -0.4,     0.9,     2.0,    -0.6,    -0.3,     1.2, &
               -0.4,  -0.4,    -0.3,     1.1,     1.1,    -0.2,     0.6, &
               -0.9,  -0.3,     0.2,    -0.3,     0.4,    -1.0,     0.3, &
               -0.5,  -0.7,    -0.4/

      DATA DH00/0.0,   0.0,   -22.5,     0.0,   -20.6,    -9.6,     0.0, &
                6.0,  -0.1,   -14.2,     0.0,     2.1,     1.3,     5.0, &
                0.3,   0.0,    -0.1,     0.6,     1.7,     1.9,     0.1, &
                0.0,  -0.2,    -1.4,     0.0,    -0.8,     0.0,     0.9, &
                0.0,   1.1,     0.0,     0.3,    -0.1,    -0.6,    -0.7, &
                0.2,   0.0,     0.1,     0.0,     0.0,     0.3,     0.6, &
               -0.4,   0.3,     0.7/
!
!
      IY=IYEAR
!
!  WE ARE RESTRICTED BY THE INTERVAL 1965-2005,
!  FOR WHICH THE IGRF COEFFICIENTS ARE KNOWN; IF IYEAR IS OUTSIDE THIS INTERVAL,
!  THE SUBROUTINE USES THE NEAREST LIMITING VALUE AND PRINTS A WARNING:
!
      IF(IY.LT.1965) THEN
       IY=1965
       WRITE (*,10) IYEAR,IY
      ENDIF

      IF(IY.GT.2005) THEN
       IY=2005
       WRITE (*,10) IYEAR,IY
      ENDIF

!
!  CALCULATE THE ARRAY REC, CONTAINING COEFFICIENTS FOR THE RECURSION RELATIONS,
!  USED IN THE IGRF SUBROUTINE FOR CALCULATING THE ASSOCIATE LEGENDRE POLYNOMIALS
!  AND THEIR DERIVATIVES:
!
      DO 20 N=1,11
         N2=2*N-1
         N2=N2*(N2-2)
         DO 20 M=1,N
            MN=N*(N-1)/2+M
20    REC(MN)=FLOAT((N-M)*(N+M-2))/FLOAT(N2)
!
      IF (IY.LT.1970) GOTO 50          !INTERPOLATE BETWEEN 1965 - 1970
      IF (IY.LT.1975) GOTO 60          !INTERPOLATE BETWEEN 1970 - 1975
      IF (IY.LT.1980) GOTO 70          !INTERPOLATE BETWEEN 1975 - 1980
      IF (IY.LT.1985) GOTO 80          !INTERPOLATE BETWEEN 1980 - 1985
      IF (IY.LT.1990) GOTO 90          !INTERPOLATE BETWEEN 1985 - 1990
      IF (IY.LT.1995) GOTO 100         !INTERPOLATE BETWEEN 1990 - 1995
      IF (IY.LT.2000) GOTO 110         !INTERPOLATE BETWEEN 1995 - 2000
!
!       EXTRAPOLATE BEYOND 2000:
!

      DT=FLOAT(IY)+FLOAT(IDAY-1)/365.25-2000.
      DO 40 N=1,66
         G(N)=G00(N)
         H(N)=H00(N)
         IF (N.GT.45) GOTO 40
         G(N)=G(N)+DG00(N)*DT
         H(N)=H(N)+DH00(N)*DT
40    CONTINUE
      GOTO 300
!
!       INTERPOLATE BETWEEEN 1965 - 1970:
!
50    F2=(FLOAT(IY)+FLOAT(IDAY-1)/365.25-1965.)/5.
      F1=1.-F2
      DO 55 N=1,66
         G(N)=G65(N)*F1+G70(N)*F2
55       H(N)=H65(N)*F1+H70(N)*F2
      GOTO 300
!
!       INTERPOLATE BETWEEN 1970 - 1975:
!
60    F2=(FLOAT(IY)+FLOAT(IDAY-1)/365.25-1970.)/5.
      F1=1.-F2
      DO 65 N=1,66
         G(N)=G70(N)*F1+G75(N)*F2
65       H(N)=H70(N)*F1+H75(N)*F2
      GOTO 300
!
!       INTERPOLATE BETWEEN 1975 - 1980:
!
70    F2=(FLOAT(IY)+FLOAT(IDAY-1)/365.25-1975.)/5.
      F1=1.-F2
      DO 75 N=1,66
         G(N)=G75(N)*F1+G80(N)*F2
75       H(N)=H75(N)*F1+H80(N)*F2
      GOTO 300
!
!       INTERPOLATE BETWEEN 1980 - 1985:
!
80    F2=(FLOAT(IY)+FLOAT(IDAY-1)/365.25-1980.)/5.
      F1=1.-F2
      DO 85 N=1,66
         G(N)=G80(N)*F1+G85(N)*F2
85       H(N)=H80(N)*F1+H85(N)*F2
      GOTO 300
!
!       INTERPOLATE BETWEEN 1985 - 1990:
!
90    F2=(FLOAT(IY)+FLOAT(IDAY-1)/365.25-1985.)/5.
      F1=1.-F2
      DO 95 N=1,66
         G(N)=G85(N)*F1+G90(N)*F2
95       H(N)=H85(N)*F1+H90(N)*F2
      GOTO 300
!
!       INTERPOLATE BETWEEN 1990 - 1995:
!
100   F2=(FLOAT(IY)+FLOAT(IDAY-1)/365.25-1990.)/5.
      F1=1.-F2
      DO 105 N=1,66
         G(N)=G90(N)*F1+G95(N)*F2
105      H(N)=H90(N)*F1+H95(N)*F2
      GOTO 300
!
!       INTERPOLATE BETWEEN 1995 - 2000:
!
110   F2=(FLOAT(IY)+FLOAT(IDAY-1)/365.25-1995.)/5.
      F1=1.-F2
      DO 115 N=1,66
         G(N)=G95(N)*F1+G00(N)*F2
115      H(N)=H95(N)*F1+H00(N)*F2
      GOTO 300
!
!   COEFFICIENTS FOR A GIVEN YEAR HAVE BEEN CALCULATED; NOW MULTIPLY
!   THEM BY SCHMIDT NORMALIZATION FACTORS:
!
300   S=1.
      DO 120 N=2,11
         MN=N*(N-1)/2+1
         S=S*FLOAT(2*N-3)/FLOAT(N-1)
         G(MN)=G(MN)*S
         H(MN)=H(MN)*S
         P=S
         DO 120 M=2,N
            AA=1.
            IF (M.EQ.2) AA=2.
            P=P*SQRT(AA*FLOAT(N-M+1)/FLOAT(N+M-2))
            MNN=MN+M-1
            G(MNN)=G(MNN)*P
120         H(MNN)=H(MNN)*P

           G10=-G(2)
           G11= G(3)
           H11= H(3)
!
!  NOW CALCULATE THE COMPONENTS OF THE UNIT VECTOR EzMAG IN GEO COORD.SYSTEM:
!   SIN(TETA0)*COS(LAMBDA0), SIN(TETA0)*SIN(LAMBDA0), AND COS(TETA0)
!         ST0 * CL0                ST0 * SL0                CT0
!
      SQ=G11**2+H11**2
      SQQ=SQRT(SQ)
      SQR=SQRT(G10**2+SQ)
      SL0=-H11/SQQ
      CL0=-G11/SQQ
      ST0=SQQ/SQR
      CT0=G10/SQR
      STCL=ST0*CL0
      STSL=ST0*SL0
      CTSL=CT0*SL0
      CTCL=CT0*CL0
!
      CALL SUN (IY,IDAY,IHOUR,MIN,ISEC,GST,SLONG,SRASN,SDEC)
!
!  S1,S2, AND S3 ARE THE COMPONENTS OF THE UNIT VECTOR EXGSM=EXGSE IN THE
!   SYSTEM GEI POINTING FROM THE EARTH'S CENTER TO THE SUN:
!
      S1=COS(SRASN)*COS(SDEC)
      S2=SIN(SRASN)*COS(SDEC)
      S3=SIN(SDEC)
      CGST=COS(GST)
      SGST=SIN(GST)
!
!  DIP1, DIP2, AND DIP3 ARE THE COMPONENTS OF THE UNIT VECTOR EZSM=EZMAG
!   IN THE SYSTEM GEI:
!
      DIP1=STCL*CGST-STSL*SGST
      DIP2=STCL*SGST+STSL*CGST
      DIP3=CT0
!
!  NOW CALCULATE THE COMPONENTS OF THE UNIT VECTOR EYGSM IN THE SYSTEM GEI
!   BY TAKING THE VECTOR PRODUCT D x S AND NORMALIZING IT TO UNIT LENGTH:
!
      Y1=DIP2*S3-DIP3*S2
      Y2=DIP3*S1-DIP1*S3
      Y3=DIP1*S2-DIP2*S1
      Y=SQRT(Y1*Y1+Y2*Y2+Y3*Y3)
      Y1=Y1/Y
      Y2=Y2/Y
      Y3=Y3/Y
!
!   THEN IN THE GEI SYSTEM THE UNIT VECTOR Z = EZGSM = EXGSM x EYGSM = S x Y
!    HAS THE COMPONENTS:
!
      Z1=S2*Y3-S3*Y2
      Z2=S3*Y1-S1*Y3
      Z3=S1*Y2-S2*Y1
!
!    THE VECTOR EZGSE (HERE DZ) IN GEI HAS THE COMPONENTS (0,-SIN(DELTA),
!     COS(DELTA)) = (0.,-0.397823,0.917462); HERE DELTA = 23.44214 DEG FOR
!   THE EPOCH 1978 (SEE THE BOOK BY GUREVICH OR OTHER ASTRONOMICAL HANDBOOKS).
!    HERE THE MOST ACCURATE TIME-DEPENDENT FORMULA IS USED:
!
      DJ=FLOAT(365*(IY-1900)+(IY-1901)/4 +IDAY) &
       -0.5+FLOAT(IHOUR*3600+MIN*60+ISEC)/86400.
      T=DJ/36525.
      OBLIQ=(23.45229-0.0130125*T)/57.2957795
      DZ1=0.
      DZ2=-SIN(OBLIQ)
      DZ3=COS(OBLIQ)
!
!  THEN THE UNIT VECTOR EYGSE IN GEI SYSTEM IS THE VECTOR PRODUCT DZ x S :
!
      DY1=DZ2*S3-DZ3*S2
      DY2=DZ3*S1-DZ1*S3
      DY3=DZ1*S2-DZ2*S1
!
!   THE ELEMENTS OF THE MATRIX GSE TO GSM ARE THE SCALAR PRODUCTS:
!   CHI=EM22=(EYGSM,EYGSE), SHI=EM23=(EYGSM,EZGSE), EM32=(EZGSM,EYGSE)=-EM23,
!     AND EM33=(EZGSM,EZGSE)=EM22
!
      CHI=Y1*DY1+Y2*DY2+Y3*DY3
      SHI=Y1*DZ1+Y2*DZ2+Y3*DZ3
      HI=ASIN(SHI)
!
!    TILT ANGLE: PSI=ARCSIN(DIP,EXGSM)
!
      SPS=DIP1*S1+DIP2*S2+DIP3*S3
      CPS=SQRT(1.-SPS**2)
      PSI=ASIN(SPS)
!
!    THE ELEMENTS OF THE MATRIX MAG TO SM ARE THE SCALAR PRODUCTS:
! CFI=GM22=(EYSM,EYMAG), SFI=GM23=(EYSM,EXMAG); THEY CAN BE DERIVED AS FOLLOWS:
!
! IN GEO THE VECTORS EXMAG AND EYMAG HAVE THE COMPONENTS (CT0*CL0,CT0*SL0,-ST0)
!  AND (-SL0,CL0,0), RESPECTIVELY.    HENCE, IN GEI THE COMPONENTS ARE:
!  EXMAG:    CT0*CL0*COS(GST)-CT0*SL0*SIN(GST)
!            CT0*CL0*SIN(GST)+CT0*SL0*COS(GST)
!            -ST0
!  EYMAG:    -SL0*COS(GST)-CL0*SIN(GST)
!            -SL0*SIN(GST)+CL0*COS(GST)
!             0
!  THE COMPONENTS OF EYSM IN GEI WERE FOUND ABOVE AS Y1, Y2, AND Y3;
!  NOW WE ONLY HAVE TO COMBINE THE QUANTITIES INTO SCALAR PRODUCTS:
!
      EXMAGX=CT0*(CL0*CGST-SL0*SGST)
      EXMAGY=CT0*(CL0*SGST+SL0*CGST)
      EXMAGZ=-ST0
      EYMAGX=-(SL0*CGST+CL0*SGST)
      EYMAGY=-(SL0*SGST-CL0*CGST)
      CFI=Y1*EYMAGX+Y2*EYMAGY
      SFI=Y1*EXMAGX+Y2*EXMAGY+Y3*EXMAGZ
!
      XMUT=(ATAN2(SFI,CFI)+3.1415926536)*3.8197186342
!
!  THE ELEMENTS OF THE MATRIX GEO TO GSM ARE THE SCALAR PRODUCTS:
!
!   A11=(EXGEO,EXGSM), A12=(EYGEO,EXGSM), A13=(EZGEO,EXGSM),
!   A21=(EXGEO,EYGSM), A22=(EYGEO,EYGSM), A23=(EZGEO,EYGSM),
!   A31=(EXGEO,EZGSM), A32=(EYGEO,EZGSM), A33=(EZGEO,EZGSM),
!
!   ALL THE UNIT VECTORS IN BRACKETS ARE ALREADY DEFINED IN GEI:
!
!  EXGEO=(CGST,SGST,0), EYGEO=(-SGST,CGST,0), EZGEO=(0,0,1)
!  EXGSM=(S1,S2,S3),  EYGSM=(Y1,Y2,Y3),   EZGSM=(Z1,Z2,Z3)
!                                                           AND  THEREFORE:
!
      A11=S1*CGST+S2*SGST
      A12=-S1*SGST+S2*CGST
      A13=S3
      A21=Y1*CGST+Y2*SGST
      A22=-Y1*SGST+Y2*CGST
      A23=Y3
      A31=Z1*CGST+Z2*SGST
      A32=-Z1*SGST+Z2*CGST
      A33=Z3
!
 10   FORMAT(//1X, &
      '**** RECALC WARNS: YEAR IS OUT OF INTERVAL 1965-2005: IYEAR=',I4, &
       /,6X,'CALCULATIONS WILL BE DONE FOR IYEAR=',I4,/)
      RETURN
      END
!
!====================================================================

!
!=====================================================================================
!*******************************************************************
!
      SUBROUTINE SUN (IYEAR,IDAY,IHOUR,MIN,ISEC,GST,SLONG,SRASN,SDEC)
!
!  CALCULATES FOUR QUANTITIES NECESSARY FOR COORDINATE TRANSFORMATIONS
!  WHICH DEPEND ON SUN POSITION (AND, HENCE, ON UNIVERSAL TIME AND SEASON)
!
!-------  INPUT PARAMETERS:
!  IYR,IDAY,IHOUR,MIN,ISEC -  YEAR, DAY, AND UNIVERSAL TIME IN HOURS, MINUTES,
!    AND SECONDS  (IDAY=1 CORRESPONDS TO JANUARY 1).
!
!-------  OUTPUT PARAMETERS:
!  GST - GREENWICH MEAN SIDEREAL TIME, SLONG - LONGITUDE ALONG ECLIPTIC
!  SRASN - RIGHT ASCENSION,  SDEC - DECLINATION  OF THE SUN (RADIANS)
!  ORIGINAL VERSION OF THIS SUBROUTINE HAS BEEN COMPILED FROM:
!  RUSSELL, C.T., COSMIC ELECTRODYNAMICS, 1971, V.2, PP.184-196.
!
!  LAST MODIFICATION:  MARCH 31, 2003 (ONLY SOME NOTATION CHANGES)
!
!     ORIGINAL VERSION WRITTEN BY:    Gilbert D. Mead
!
!!!   DOUBLE PRECISION DJ,FDAY
      DOUBLE PRECISION DJ,FDAY,RAD
      DATA RAD/57.295779513d0/
!
      IF(IYEAR.LT.1901.OR.IYEAR.GT.2099) RETURN
      FDAY=DFLOAT(IHOUR*3600+MIN*60+ISEC)/86400.D0
! jah, foresys : Integer division in a non-integer context
! jah, DJ=365*(IYEAR-1900)+(IYEAR-1901)/4+IDAY-0.5D0+FDAY
      ItempDJ=(IYEAR-1901)/4
      DJ=365.d0*dble(IYEAR-1900)+dble(ItempDJ)+dble(IDAY)-0.5D0+FDAY
! jah, foresys : Precision loss in assignment from double to real
! jah, T=DJ/36525.
      T=sngl(DJ)/36525.
! jah, foresys : Precision loss in assignment from real*8 to real
! jah, VL=DMOD(279.696678+0.9856473354*DJ,360.D0)
      VL=sngl(DMOD(279.696678d0+0.9856473354d0*DJ,360.D0))
! jah, foresys : Precision loss in assignment from real*8 to real
! jah, GST=DMOD(279.690983+.9856473354*DJ+360.*FDAY+180.,360.D0)/RAD
      GST=sngl(DMOD(279.690983d0+.9856473354d0*DJ+360.d0*FDAY+180.d0,360.D0)/ &
           RAD)
! jah, foresys : Precision loss in assignment from real*8 to real
! jah, G=DMOD(358.475845+0.985600267*DJ,360.D0)/RAD
      G=sngl(DMOD(358.475845d0+0.985600267d0*DJ,360.D0)/RAD)
      SLONG=(VL+(1.91946-0.004789*T)*SIN(G)+0.020094*SIN(2.*G))/sngl(RAD)
      IF(SLONG.GT.6.2831853) SLONG=SLONG-6.2831853
      IF (SLONG.LT.0.) SLONG=SLONG+6.2831853
      OBLIQ=(23.45229-0.0130125*T)/sngl(RAD)
      SOB=SIN(OBLIQ)
      SLP=SLONG-9.924E-5
!
!   THE LAST CONSTANT IS A CORRECTION FOR THE ANGULAR ABERRATION  DUE TO
!   THE ORBITAL MOTION OF THE EARTH
!
      SIND=SOB*SIN(SLP)
      COSD=SQRT(1.-SIND**2)
      SC=SIND/COSD
      SDEC=ATAN(SC)
      SRASN=3.141592654-ATAN2(COS(OBLIQ)/SOB*SC,-COS(SLP)/COSD)
      RETURN
      END
!
!================================================================================
