0005c0005c0057c	This collection of subroutines was adapted from the0059c	GEOPACK.F set by N.A. Tsyganenko for use by M. Peredo0060c	at the Science Planning and Operations Facility (SPOF)0026c	of the ISTP program.0005c0065c	Please see the additional file geopack.doc for descriptions0044c	of what the individual subroutines do.0005c0060c	Major changes include: (i) modification of subroutine 0063c	IGRF to accept dates beyond 1990, (ii) adoption of RECOMP0059c	instead of RECALC (RECOMP has been renamed RECALC and0060c	the original RECALC has been removed from the packet),0067c	(iii) addition of new routine BTOT to compute the total field0063c	at an array of input points [X(count),Y(count),Z(count)],0066c	(iv) consolidation of old an new versions of the models into0063c	a single set, and modification of the models nomenclature0061c	to match the new nomenclature proposed in Appendix A of0065c	"Are Existing Magnetospheric Models Excessively Stretched?"0061c	by M. Peredo, D.P. Stern and N.A. Tsyganenko, JGR, 19930058c	(in press).  According to this new nomenclature, the0027c	available models are:0005c0065c	MODEL	  Based on    Activity levels   Range of   Subroutine0061c	          data from   and criterion     validity   name0005c0061c	T87a	  IMP-HEOS       6 by Kp        X = -70 Re  EXT87L0064c	T87b      IMP-HEOS       8 by Kp        X = -30 Re  EXT87T0064c	T87Wd     IMP-HEOS-ISEE  6 by Kp        X = -70 Re  EXT87W0064c	T87We     IMP-HEOS-ISEE  12 by Psw      X = -30 Re  EXT87W0020c				 and IMF-Bz0050c	T89a	  IMP-HEOS	 6 by Kp	X = -70 Re  EXT89KP0050c	T89b	  IMP-HEOS	 6 by AE	X = -70 Re  EXT89AE0055c	T89c	  IMP-HEOS-ISEE  7 by Kp 	X = -70 Re  EXT89M0005c0061c	NOTE that T87b and T87We are ONLY valid to X=-30 Re and0051c	SHOULD NOT be used tailward of that distance.0005c0005c0047c	Compiled and modified by: Mauricio Peredo0034c				  Hughes STX at NASA/GSFC0029c				  NASA/GSFC Code 6950030c				  Greenbelt, MD 207710014c				  USA0033c				  (301) 286-1526 (voice)0031c				  (301) 286-3271 (FAX)0018c				  E-mail:^^^^^^^^^^^^^^^^^0032c				  ISTP::PEREDO (DECnet)0048c				  peredo@istp1.gsfc.nasa.gov (Internet)0025c				  September 19920005c0005c0056c	February 19, 1993  Modified (M. Peredo) to add new0054c	warped T87 models (subroutine EXT87W) derived by0057c	M. Peredo for different solar wind dynamic pressure0057c	and IMF-Bz conditions (12 models) and for different0027c	Kp levels (6 models).0005c0062c	June 17, 1993  Modified (M. Peredo) to add new T89 model0062c	and revise introductory information to include reference0033c	to new nomenclature scheme.0005c0005c0005c0077C------------------------------------------------------------------------0061        SUBROUTINE BTOT(X,Y,Z,EXNAME,IOPT,BX,BY,BZ,count)0005C0005C0043C       Written by:     Mauricio Peredo0061C                       Hughes STX at NASA/GSFC, Code 6950042C                       September 19920005C00040066C       This routine computes the total magnetic field at each0067C       (input) point (Xi,Yi,Zi). It is assumed that subroutine0063C       RECALC has already been called to set quantities in0063C       the common block (in particular, the call to RECALC0065C       sets the dipole tilt angle, and all angles needed for0064C       conversions between various coordinate systems). The0063C       internal field is computed from the IGRF model (for0064C       appropriate year). In calling IGRF, we set the order0065C       of harmonics to use according to the following rules:0058C       if the radius is 10 Re or more we use a dipole0061C       approximation (IHARM=1), if the radius is between0063C       6 Re and 10 Re, we use dipole+quadrupole (IHARM=2),0063C       inside 6 Re we go for the full blown expansion with0062C       IHARM=10.  The external field is computed from the0062C       routine identified by the current value of EXNAME,0063C       and the "version" of that model selected depends on0062C       the value of the flag IOPT [For example if we have0067C       EXNAME=EX89KP and IOPT=1 the external field computation^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^0069C       is according to the T89 model with Kp = 0, 0+, and so on]0005C0005C0079C------ INPUT:  X,Y,Z = Arrays of Cartesian GSM coordinates with the points0069C                       where the field is to be computed (in Re)0074C               EXNAME = Name of the subroutine to use for computation0077C                        of the external field. Valid values are: EXT87L,0078C                        EXT87T, EX89KP, EX89AE (more will be added later)0074C               IOPT = Flag indication which "version" of the model to0077C                      use, the options correspond to different levels of0078C                      magnetospheric activity.  EXT87L, EX89KP and EX89AE0078C                      each allow IOPT = 1 to 6, EXT87T has IOPT = 1 to 8,0069c		       EXT87W has IOPT = 101 to 112 for models binned by solar0065c		       wind dynamic pressure and IMF-Bz, and IOPT = 1 to 60044c		       for models binned by Kp level.0041c		count =  Dimension of X,Y,Z arrays0005C0078C------ OUTPUT: BX,BY,BZ = Cartesian GSM components of the total (internal0078C                          plus external) magnetic field vector at (X,Y,Z)0005C0075C       ATTENTION: IT IS ASSUMED THAT SUBROUTINE RECALC HAS BEEN CALLED0062C                  PRIOR TO CALLING THIS SUBROUTINE!!!!!!!0005C0005C0055C       Duplicate the common block from RECOMP here0005C0018	IMPLICIT NONE00040036	INTEGER IY,K,count,i,iopt,IHARM00040020	EXTERNAL EXNAME00040062	REAL ST0,CT0,SL0,CL0,CTCL,STCL,CTSL,STSL,SFI,CFI,SPS,CPS,0069     1   SHI,CHI,HI,PSI,XMUT,A11,A21,A31,A12,A22,A32,A13,A23,A33,0064     2   DS3,X(count),Y(count),Z(count),BX(count),BY(count),0068     3   BZ(count),BXEXT,BYEXT,BZEXT,XGEO,YGEO,ZGEO,R,T,F,FX,FY,0057     4   FZ,BXINT,BYINT,BZINT,BR,BF,BT,BA(8),HX,HY,HZ00040074       COMMON/C1/ ST0,CT0,SL0,CL0,CTCL,STCL,CTSL,STSL,SFI,CFI,SPS,CPS,0071     * SHI,CHI,HI,PSI,XMUT,A11,A21,A31,A12,A22,A32,A13,A23,A33,DS3,0018     * K,IY,BA00040031        do 100 i = 1, count00040034C       Compute external field0004^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^0069           CALL EXNAME(IOPT,PSI,X(i),Y(i),Z(i),BXEXT,BYEXT,BZEXT)0005c0005C0078C       To compute internal field, we convert input point to GEO cartesian0078C       coordinates, then to GEO spherical coordinates, call IGRF with the0078C       appropriate value for IHARM (see comments above), then convert the0080C       resulting internal field to cartesian GEO coordinates and finally to0038C       cartesian GSM coordinates.0005C0060           CALL GEOGSM(XGEO,YGEO,ZGEO,X(i),Y(i),Z(i),-1)0051           CALL SPHCAR(R,T,F,XGEO,YGEO,ZGEO,-1)0038           IF (R.GE.9.999999) THEN0027              IHARM = 10042           ELSEIF (R.GE.5.999999) THEN0027              IHARM = 20019           ELSE0028              IHARM = 100020           ENDIF0049           CALL IGRF(IY,IHARM,R,T,F,BR,BT,BF)0058           CALL BSPCAR(T,F,BR,BT,BF,BXINT,BYINT,BZINT)0049	   call GEOGSM(BXINT,BYINT,BZINT,HX,HY,HZ,1)0005C0054C       Sum up internal and external contributions0005C0033           BX(i) = HX + BXEXT0033           BY(i) = HY + BYEXT0033           BZ(i) = HZ + BZEXT00040020 100    continue00040018        return0015        end0005C0005c0005c0047      SUBROUTINE IGRF(IY,NM,R,T,F,BR,BT,BF)0005C0005c0066C	MODIFIED TO ACCEPT DATES BETWEEN 1965 AND 2000; COEFFICIENTS0065C	FOR IGRF 1985 HAVE BEEN REPLACED WITH DGRF1985 COEFFICIENTS0067C	[EOS TRANS. AGU APRIL 21, 1992, P. 182].  ALSO, CODE MODIFIED0066C	TO ACCEPT DATES BEYOND 1990, AND TO USE LINEAR EXTRAPOLATION0067C	BETWEEN 1990 AND 2000 BASED ON THE IGRF COEFFICIENTS FROM THE0022C	SAME EOS ARTICLE0005c0005c0034c	Modified by:	Mauricio Peredo0031c			Hughes STX at NASA/GSFC0022c			September 19920005c0005c0067C  CALCULATES COMPONENTS OF MAIN GEOMAGNETIC FIELD IN SPHERICAL0069C  GEOGRAPHICAL COORD SYSTEM BY USING THIRD GENERATION IGRF MODEL0053C  (J. GEOMAG. GEOELECTR.(1982), V.34, P.313-315,0054C  GEOMAGN. AND AERONOMY (1986), V.26, P.523-525).0073C  UPDATING OF COEFFICIENTS TO A GIVEN EPOCH IS MADE DURING THE FIRST0051C  CALL AND AFTER EVERY CHANGE OF PARAMETER IY.^^^^^^^^^^^^^^^^^^^^^0028C------INPUT PARAMETERS:0046C  IY - YEAR NUMBER (FROM 1965 UP TO 1990)0076C  NM - MAXIMAL ORDER OF HARMONICS TAKEN INTO ACCOUNT (NOT MORE THAN 10)0075C  R,T,F - SPHERICAL COORDINATES OF THE POINT (R IN UNITS RE=6371.2 KM,0047C  COLATITUDE T AND LONGITUDE F IN RADIANS)0029C----- OUTPUT PARAMETERS:0073C  BR,BT,BF - SPHERICAL COMPONENTS OF MAIN GEOMAGN.FIELD IN NANOTESLA0005C0005C0053C                   AUTHOR: NIKOLAI A. TSYGANENKO0052C                           INSTITUTE OF PHYSICS0063C                           ST.-PETERSBURG STATE UNIVERSITY0053C                           STARY PETERGOF 1989040046C                           ST.-PETERSBURG0038C                           RUSSIA0005C0018	IMPLICIT NONE0005C0071      REAL A(11),B(11),G(66),H(66),REC(66),G65(66),H65(66),G70(66),0074     *H70(66),G75(66),H75(66),G80(66),H80(66),G85(66),H85(66),G90(66),0035     *H90(66),DG90(45),DH90(45)00040062	REAL R,T,F,BR,BT,BF,DT,F2,F1,S,P,AA,PP,D,BBR,BBF,U,CF,SF,0061     1       C,W,X,Y,Z,Q,BI,P2,D2,AN,E,HH,BBT,QQ,XK,DP,PM00040023      LOGICAL BK,BM00040052	INTEGER IY,NM,MA,IPR,IYR,KNM,N,N2,M,MNN,MN,K,MM0005c0075      DATA G65/0.,-30334.,-2119.,-1662.,2997.,1594.,1297.,-2038.,1292.,0073     *856.,957.,804.,479.,-390.,252.,-219.,358.,254.,-31.,-157.,-62.,0075     *45.,61.,8.,-228.,4.,1.,-111.,75.,-57.,4.,13.,-26.,-6.,13.,1.,13.,0076     *5.,-4.,-14.,0.,8.,-1.,11.,4.,8.,10.,2.,-13.,10.,-1.,-1.,5.,1.,-2.,0047     *-2.,-3.,2.,-5.,-2.,4.,4.,0.,2.,2.,0./0074      DATA H65/0.,0.,5776.,0.,-2016.,114.,0.,-404.,240.,-165.,0.,148.,0075     *-269.,13.,-269.,0.,19.,128.,-126.,-97.,81.,0.,-11.,100.,68.,-32.,0074     *-8.,-7.,0.,-61.,-27.,-2.,6.,26.,-23.,-12.,0.,7.,-12.,9.,-16.,4.,0073     *24.,-3.,-17.,0.,-22.,15.,7.,-4.,-5.,10.,10.,-4.,1.,0.,2.,1.,2.,0034     *6.,-4.,0.,-2.,3.,0.,-6./0005c0075      DATA G70/0.,-30220.,-2068.,-1781.,3000.,1611.,1287.,-2091.,1278.,0073     *838.,952.,800.,461.,-395.,234.,-216.,359.,262.,-42.,-160.,-56.,0073     *43.,64.,15.,-212.,2.,3.,-112.,72.,-57.,1.,14.,-22.,-2.,13.,-2.,^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^0072     *14.,6.,-2.,-13.,-3.,5.,0.,11.,3.,8.,10.,2.,-12.,10.,-1.,0.,3.,0055     *1.,-1.,-3.,-3.,2.,-5.,-1.,6.,4.,1.,0.,3.,-1./0073      DATA H70/0.,0.,5737.,0.,-2047.,25.,0.,-366.,251.,-196.,0.,167.,0075     *-266.,26.,-279.,0.,26.,139.,-139.,-91.,83.,0.,-12.,100.,72.,-37.,0073     *-6.,1.,0.,-70.,-27.,-4.,8.,23.,-23.,-11.,0.,7.,-15.,6.,-17.,6.,0073     *21.,-6.,-16.,0.,-21.,16.,6.,-4.,-5.,10.,11.,-2.,1.,0.,1.,1.,3.,0034     *4.,-4.,0.,-1.,3.,1.,-4./0005c0075      DATA G75/0.,-30100.,-2013.,-1902.,3010.,1632.,1276.,-2144.,1260.,0073     *830.,946.,791.,438.,-405.,216.,-218.,356.,264.,-59.,-159.,-49.,0072     *45.,66.,28.,-198.,1.,6.,-111.,71.,-56.,1.,16.,-14.,0.,12.,-5.,0073     *14.,6.,-1.,-12.,-8.,4.,0.,10.,1.,7.,10.,2.,-12.,10.,-1.,-1.,4.,0055     *1.,-2.,-3.,-3.,2.,-5.,-2.,5.,4.,1.,0.,3.,-1./0074      DATA H75/0.,0.,5675.,0.,-2067.,-68.,0.,-333.,262.,-223.,0.,191.,0074     *-265.,39.,-288.,0.,31.,148.,-152.,-83.,88.,0.,-13.,99.,75.,-41.,0075     *-4.,11.,0.,-77.,-26.,-5.,10.,22.,-23.,-12.,0.,6.,-16.,4.,-19.,6.,0074     *18.,-10.,-17.,0.,-21.,16.,7.,-4.,-5.,10.,11.,-3.,1.,0.,1.,1.,3.,0035     *4.,-4.,-1.,-1.,3.,1.,-5./0005c0075      DATA G80/0.,-29992.,-1956.,-1997.,3027.,1663.,1281.,-2180.,1251.,0073     *833.,938.,782.,398.,-419.,199.,-218.,357.,261.,-74.,-162.,-48.,0073     *48.,66.,42.,-192.,4.,14.,-108.,72.,-59.,2.,21.,-12.,1.,11.,-2.,0074     *18.,6.,0.,-11.,-7.,4.,3.,6.,-1.,5.,10.,1.,-12.,9.,-3.,-1.,7.,2.,0051     *-5.,-4.,-4.,2.,-5.,-2.,5.,3.,1.,2.,3.,0./0075      DATA H80/0.,0.,5604.,0.,-2129.,-200.,0.,-336.,271.,-252.,0.,212.,0074     *-257.,53.,-297.,0.,46.,150.,-151.,-78.,92.,0.,-15.,93.,71.,-43.,0075     *-2.,17.,0.,-82.,-27.,-5.,16.,18.,-23.,-10.,0.,7.,-18.,4.,-22.,9.,0073     *16.,-13.,-15.,0.,-21.,16.,9.,-5.,-6.,9.,10.,-6.,2.,0.,1.,0.,3.,0034     *6.,-4.,0.,-1.,4.,0.,-6./0005c0075      DATA G85/0.,-29873.,-1905.,-2072.,3044.,1687.,1296.,-2208.,1247.,0073     *829.,936.,780.,361.,-424.,170.,-214.,355.,253.,-93.,-164.,-46.,^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^0075     *53.,65.,51.,-185.,4.,16.,-102.,74.,-62.,3.,24.,-6.,4.,10.,0.,21.,0074     *6.,0.,-11.,-9.,4.,4.,4.,-4.,5.,10.,1.,-12.,9.,-3.,-1.,7.,1.,-5.,0047     *-4.,-4.,3.,-5.,-2.,5.,3.,1.,2.,3.,0./0075      DATA H85/0.,0.,5500.,0.,-2197.,-306.,0.,-310.,284.,-297.,0.,232.,0074     *-249.,69.,-297.,0.,47.,150.,-154.,-75.,95.,0.,-16.,88.,69.,-48.,0075     *-1.,21.,0.,-83.,-27.,-2.,20.,17.,-23.,-7.,0.,8.,-19.,5.,-23.,11.,0072     *14.,-15.,-11.,0.,-21.,15.,9.,-6.,-6.,9.,9.,-7.,2.,0.,1.,0.,3.,0034     *6.,-4.,0.,-1.,4.,0.,-6./0005c0075      DATA G90/0.,-29775.,-1851.,-2136.,3058.,1693.,1315.,-2240.,1246.,0074     *807.,939.,782.,324.,-423.,142.,-211.,353.,244.,-111.,-166.,-37.,0073     *61.,64.,60.,-178.,2.,17.,-96.,77.,-64.,4.,28.,1.,6.,10.,0.,22.,0076     *5.,-1.,-11.,-12.,4.,4.,3.,-6.,4.,10.,1.,-12.,9.,-4.,-1.,7.,2.,-6.,0047     *-4.,-4.,2.,-5.,-2.,4.,3.,1.,2.,3.,0./0075      DATA H90/0.,0.,5411.,0.,-2278.,-380.,0.,-287.,293.,-348.,0.,248.,0074     *-240.,87.,-299.,0.,47.,153.,-154.,-69.,98.,0.,-16.,83.,68.,-52.,0074     *2.,27.,0.,-81.,-27.,1.,20.,16.,-23.,-5.,0.,10.,-20.,7.,-22.,12.,0073     *11.,-16.,-11.,0.,-21.,15.,10.,-6.,-6.,9.,9.,-7.,2.,0.,1.,0.,3.,0034     *6.,-4.,0.,-1.,4.,0.,-6./0074      DATA DG90 /0.,18.0,10.6,-12.9,2.4,0.0,3.3,-6.7,0.1,-5.9,0.5,0.6,0074     *-7.0,0.5,-5.5,0.6,-0.1,-1.6,-3.1,-0.1,2.3,1.3,-0.2,1.8,1.3,-0.2,0075     *0.1,1.2,0.6,-0.5,-0.3,0.6,1.6,0.2,0.2,0.3,0.2,-0.7,-0.2,0.1,-1.1,0029     *0.0,-0.1,-0.5,-0.6/0076      DATA DH90 /0.,0.,-16.1,0.,-15.8,-13.8,0.,4.4,1.6,-10.6,0.,2.6,1.8,0072     *3.1,-1.4,0.,-0.1,0.5,0.4,1.7,0.4,0.,0.2,-1.3,0.0,-0.9,0.5,1.2,0072     *0.,0.6,0.2,0.8,-0.5,-0.2,0.0,0.0,0.,0.5,-0.2,0.3,0.3,0.4,-0.5,0019     *-0.3,0.6/0005c0032      DATA MA,IYR,IPR/0,0,0/0029      IF(MA.NE.1) GOTO 100031      IF(IY.NE.IYR) GOTO 300018      GOTO 130001410    MA=10016      KNM=150005C0022      DO 20 N=1,110021         N2=2*N-10025         N2=N2*(N2-2)0024         DO 20 M=1,N0030            MN=N*(N-1)/2+M004820    REC(MN)=FLOAT((N-M)*(N+M-2))/FLOAT(N2)0005C001630    IYR=IY^^^^^^^^^^0035      IF (IYR.LT.1965) IYR=19650035      IF (IYR.GT.2000) IYR=20000056      IF (IY.NE.IYR.AND.IPR.EQ.0) write(*,999)IY,IYR0030      IF (IYR.NE.IY) IPR=10068      IF (IYR.LT.1970) GOTO 50		!INTERPOLATE BETWEEN 1965 - 19700068      IF (IYR.LT.1975) GOTO 60		!INTERPOLATE BETWEEN 1970 - 19750068      IF (IYR.LT.1980) GOTO 70		!INTERPOLATE BETWEEN 1975 - 19800068      IF (IYR.LT.1985) GOTO 80		!INTERPOLATE BETWEEN 1980 - 19850068      IF (IYR.LT.1990) GOTO 90		!INTERPOLATE BETWEEN 1985 - 19900005C0005C0037C	EXTRAPOLATE BETWEEN 1990 - 20000005C0029      DT=FLOAT(IYR)-1990.0022      DO 40 N=1,660024         G(N)=G90(N)0024         H(N)=H90(N)0033         IF (N.GT.45) GOTO 400033         G(N)=G(N)+DG90(N)*DT0033         H(N)=H(N)+DH90(N)*DT001840    CONTINUE0018      GOTO 3000005C0005C0038C	INTERPOLATE BETWEEEN 1965 - 19700005C002650    F2=(IYR-1965)/5.0018      F1=1.-F20022      DO 55 N=1,660037         G(N)=G65(N)*F1+G70(N)*F2003755       H(N)=H65(N)*F1+H70(N)*F20018      GOTO 3000005C0005C0037C	INTERPOLATE BETWEEN 1970 - 19750005C002660    F2=(IYR-1970)/5.0018      F1=1.-F20022      DO 65 N=1,660037         G(N)=G70(N)*F1+G75(N)*F2003765       H(N)=H70(N)*F1+H75(N)*F20018      GOTO 3000005C0005C0037C	INTERPOLATE BETWEEN 1975 - 19800005C002670    F2=(IYR-1975)/5.0018      F1=1.-F20022      DO 75 N=1,660037         G(N)=G75(N)*F1+G80(N)*F2003775       H(N)=H75(N)*F1+H80(N)*F20018      GOTO 3000005C0005C0037C	INTERPOLATE BETWEEN 1980 - 19850005C002680    F2=(IYR-1980)/5.0018      F1=1.-F20022      DO 85 N=1,660037         G(N)=G80(N)*F1+G85(N)*F2003785       H(N)=H80(N)*F1+H85(N)*F20018      GOTO 3000005C0005C0037C	INTERPOLATE BETWEEN 1985 - 19900005C002690    F2=(IYR-1985)/5.0018      F1=1.-F20022      DO 95 N=1,660037         G(N)=G85(N)*F1+G90(N)*F2003795       H(N)=H85(N)*F1+H90(N)*F20018      GOTO 3000005C0005C0078C	GET HERE WHEN COEFFICIENTS FOR APPROPRIATE IGRF MODEL HAVE BEEN ASSIGNED0005C0014300   S=1.0023      DO 120 N=2,110027         MN=N*(N-1)/2+10040         S=S*FLOAT(2*N-3)/FLOAT(N-1)^^^^^^^^^^^^^^^^^^0026         G(MN)=G(MN)*S0026         H(MN)=H(MN)*S0016         P=S0025         DO 120 M=2,N0021            AA=1.0033            IF (M.EQ.2) AA=2.0054            P=P*SQRT(AA*FLOAT(N-M+1)/FLOAT(N+M-2))0026            MNN=MN+M-10031            G(MNN)=G(MNN)*P0031120         H(MNN)=H(MNN)*P0005C0033130   IF(KNM.EQ.NM) GO TO 1400016      KNM=NM0017      K=KNM+10017140   PP=1./R0014      P=PP0022      DO 150 N=1,K0019         P=P*PP0019         A(N)=P0021150      B(N)=P*N0014      P=1.0014      D=0.0016      BBR=0.0016      BBT=0.0016      BBF=0.0013      U=T0019      CF=COS(F)0019      SF=SIN(F)0018      C=COS(U)0018      S=SIN(U)0025      BK=(S.LT.1.E-5)0022      DO 200 M=1,K0024         BM=(M.EQ.1)0028         IF(BM) GOTO 1600019         MM=M-10016         W=X0024         X=W*CF+Y*SF0024         Y=Y*CF-W*SF0021         GOTO 1700017160      X=0.0017         Y=1.0016170      Q=P0016         Z=D0018         BI=0.0018         P2=0.0018         D2=0.0025         DO 190 N=M,K0023            AN=A(N)0030            MN=N*(N-1)/2+M0023            E=G(MN)0024            HH=H(MN)0026            W=E*Y+HH*X0032            BBR=BBR+B(N)*W*Q0030            BBT=BBT-AN*W*Z0031            IF(BM) GOTO 1800020            QQ=Q0027            IF(BK) QQ=Z0038            BI=BI+AN*(E*X-HH*Y)*QQ0026180         XK=REC(MN)0032            DP=C*Z-S*Q-XK*D20028            PM=C*Q-XK*P20020            D2=Z0020            P2=Q0020            Z=DP0019190        Q=PM0022         D=S*D+C*P0018         P=S*P0028         IF(BM) GOTO 2000021         BI=BI*MM0023         BBF=BBF+BI0018200   CONTINUE0005C0016      BR=BBR0016      BT=BBT0025      IF(BK) GOTO 2100018      BF=BBF/S0018      GOTO 2200030210   IF(C.LT.0.) BBF=-BBF0016      BF=BBF0018220   CONTINUE0016      RETURN0005C0022999   FORMAT(//1X,0074     * 'IGRF WARNS:**** YEAR IS OUT OF INTERVAL 1965-2000: IY =',I5,/,0070     *',        CALCULATIONS WILL BE DONE FOR IYR =',I5,' ****'//)0013      END0005C0005C0074c---------------------------------------------------------------------0005c^^^^^^^^^^^^^^^0044       SUBROUTINE DIP(PS,X,Y,Z,BX,BY,BZ)0005C0074C  CALCULATES GSM COMPONENTS OF GEODIPOLE FIELD WITH THE DIPOLE MOMENT0035C  CORRESPONDING TO EPOCH=1980.0034C------------INPUT PARAMETERS:0075C   PS - GEODIPOLE TILT ANGLE IN RADIANS, X,Y,Z - GSM COORDINATES IN RE0035C------------OUTPUT PARAMETERS:0064C   BX,BY,BZ - FIELD COMPONENTS IN GSM SYSTEM, IN NANOTESLA.0005C0005C0053C                   AUTHOR: NIKOLAI A. TSYGANENKO0052C                           INSTITUTE OF PHYSICS0063C                           ST.-PETERSBURG STATE UNIVERSITY0053C                           STARY PETERGOF 1989040046C                           ST.-PETERSBURG0038C                           RUSSIA0005C0018	IMPLICIT NONE0005C0049	REAL PS,X,Y,Z,BX,BY,BZ,PSI,SPS,CPS,P,U,V,T,Q0014	INTEGER M00040026      DATA M,PSI/0,5./0052      IF(M.EQ.1.AND.ABS(PS-PSI).LT.1.E-5) GOTO 10021      SPS=SIN(PS)0021      CPS=COS(PS)0016      PSI=PS0013      M=10016  1   P=X**20016      U=Z**20018      V=3.*Z*X0016      T=Y**20033      Q=30574./SQRT(P+T+U)**50037      BX=Q*((T+U-2.*P)*SPS-V*CPS)0034      BY=-3.*Y*Q*(X*SPS+Z*CPS)0037      BZ=Q*((P+T-2.*U)*CPS-V*SPS)0016      RETURN0013      END0077c------------------------------------------------------------------------0005c0005c0005C0051      SUBROUTINE EXT87L(IOPT,PS,X,Y,Z,BX,BY,BZ)0005C0068C      'LONG' VERSION OF THE MAGNETOSPHERIC MAGNETIC FIELD MODEL0072C  COMPUTES GSM COMPONENTS OF THE MAGNETIC FIELD OF EXTRATERRESTRIAL0070C  CURRENT SYSTEMS UP TO GEOCENTRIC RADIAL DISTANCES ABOUT 70 RE .0073C  CORRESPONDS TO MAGNETOSPHERIC FIELD MODEL (N.A.TSYGANENKO, PLANET.0038C  SPACE SCI., V.35, P.1347, 1987)0055C  BASED ON IMP-A,C,D,E,F,G,H,I,J (YEARS 1966-1974)0058C  AND HEOS-1,2 (1969-1974) SATELLITE MERGED DATA SET.0074C----- INPUT PARAMETERS:  IOPT - SPECIFIES MODEL VERSION, ACCORDING TO0056C    GROUND DISTURBANCE LEVEL, IN THE FOLLOWING WAY:0005C0064C  IOPT=   1       2         3          4         5        60028C         CORRESPOND TO:0066C   KP=   0,0+   1-,1,1+   2-,2,2+    3-,3,3+   4-,4,4+   >=5-0005C^^^^^^^0045C    PS - GEODIPOLE TILT ANGLE IN RADIANS0060C    X,Y,Z - GSM COORDINATES OF THE POINT IN EARTH RADII0074C------ OUTPUT PARAMETERS:   BX,BY,BZ - GSM COMPONENTS OF THE EXTERNAL0036C    MAGNETIC FIELD IN NANOTESLA0005C0018	IMPLICIT NONE00040064	REAL PS,X,Y,Z,BX,BY,BZ,PSI,HPI,FC,RT,X1,X2,C1,RRC2,DSTR,XN,0070     1     RH,DY,DELX,B0,B1,B2,XN21,XN2,XNR,XN22,ADLN,SPS,CPS,RPS,0070     2     ZS,ZP,ZM,FY,XNX,XNX2,XC1,XC2,XC22,XR2,XC12,B20,B2P,B2M,0070     3     B,BP,BM,XA1,XAP1,XAM1,XA2,XAP2,XAM2,XNA,XNAP,XNAM,F,FP,0070     4     FM,XLN1,XLNP1,XLNM1,XLN2,XLNP2,XLNM2,ALN,S0,S0P,S0M,S1,0070     5     S1P,S1M,S2,S2P,S2M,G1,G1P,G1M,G2,G2P,G2M,EX1,EX2,Y2,Z2,0054     6     YZ,XSM,ZSM,RR2,RR,ZN,BRSM,BZSM,BY1,BXSM00040022	INTEGER IOPT,IP,I00040069      REAL GA1(32),GA2(32),GA3(32),GA4(32),GA5(32),GA6(32),PA(32)0075      DATA GA1/-.09673,-10.63,1.21,34.57,-.04502,-.06553,-.02952,.3852,0075     *-.03665,-2.084,.001795,.00638,-23.49,.06082,.01642,-.02137,32.21,0071     *-.04373,-.02311,-.2832,-.002303,-.000631,-6.397,-967.,-8650.,0053     *-20.55,5.18,-2.796,2.715,13.58,8.038,29.21/0073      DATA GA2/-.485,-12.84,1.856,40.06,-.0294,-.09071,-.02993,.5465,0075     *-.04928,-2.453,.001587,.007402,-29.41,.08101,.02322,-.1091,40.75,0072     *-.07995,-.03859,-.2755,-.002759,-.000408,-6.189,-957.8,-7246.,0054     *-25.51,5.207,-4.184,2.641,16.56,7.795,29.36/0075      DATA GA3/-1.132,-18.05,2.625,48.55,-.004868,-.1087,-.03824,.8514,0075     *-.0522,-2.881,-.000295,.009055,-29.48,.06394,.03864,-.2288,41.77,0071     *-.05849,-.06443,-.4683,.001222,-.000519,-3.696,-991.1,-6955.,0054     *-31.43,4.878,-3.151,3.277,19.19,7.248,28.99/0071      DATA GA4/-1.003,-16.98,3.14,52.81,-.08625,-.1478,-.03501,.55,0073     *-.07778,-2.97,.002086,.01275,-26.79,.06328,.03622,.08345,39.72,0071     *-.06009,-.07825,-.9698,.000178,-.000573,-.9328,-872.5,-5851.,0053     *-39.68,4.902,-3.848,2.79,20.91,6.193,26.81/0075      DATA GA5/-1.539,-14.29,3.479,53.36,-.004201,-.2043,-.03932,.6409,0072     *-.1058,-3.221,-.00114,.02166,-30.43,.04049,.05464,.008884,42.,^0075     *-.01035,-.1053,-1.63,.003802,-.001029,4.204,-665.6,-1011.,-43.49,0044     *4.514,-2.948,2.99,21.59,6.005,22./0072      DATA GA6/-2.581,-7.726,5.045,53.31,.02262,-.1972,-.01981,.428,0073     *-.1055,-5.075,.002762,.03277,-27.35,.04986,.06119,-.1211,47.48,0070     *-.0502,-.1477,.838,-.01008,-.0057,9.231,-674.3,-900.,-74.43,0044     *4.658,-3.245,3.39,21.8,5.62,25.17/0069      DATA IP,PSI,HPI,FC,RT,X1,X2/100,10.,1.5707963,0.3183099031,0020     *30.,4.,5./0032      IF (IOPT.EQ.IP) GOTO 20017      IP=IOPT0021      DO 1 I=1,320035      IF (IP.EQ.1) PA(I)=GA1(I)0035      IF (IP.EQ.2) PA(I)=GA2(I)0035      IF (IP.EQ.3) PA(I)=GA3(I)0035      IF (IP.EQ.4) PA(I)=GA4(I)0035      IF (IP.EQ.5) PA(I)=GA5(I)0035      IF (IP.EQ.6) PA(I)=GA6(I)0018  1   CONTINUE0022      C1=PA(29)**20024      RRC2=PA(27)**20029      DSTR=PA(26)/RRC2*4.0019      XN=PA(28)0019      RH=PA(31)0019      DY=PA(30)0021      DELX=PA(32)0019      B0=PA(23)0019      B1=PA(24)0019      B2=PA(25)0025      XN21=(XN-X1)**20019      XN2=XN-X20020      XNR=1./XN20021      XN22=XN2**20030      ADLN=ALOG(XN22/XN21)0041  2   IF(ABS(PS-PSI).LT.1.E-6) GOTO 30016      PSI=PS0021      SPS=SIN(PS)0021      CPS=COS(PS)0020      RPS=RH*SPS0005C0076C   COMPUTATION BEGINS HERE IF PARAMETERS IOPT AND PS REMAINED UNCHANGED0051C   AFTER THE PRECEDING CALL OF THIS SUBROUTINE0005C0018  3   ZS=Z-RPS0017      ZP=Z-RT0017      ZM=Z+RT0030      FY=FC/(1.+(Y/DY)**2)0018      XNX=XN-X0021      XNX2=XNX**20018      XC1=X-X10018      XC2=X-X20021      XC22=XC2**20021      XR2=XC2*XNR0021      XC12=XC1**20022      B20=ZS**2+C10022      B2P=ZP**2+C10022      B2M=ZM**2+C10021      B=SQRT(B20)0022      BP=SQRT(B2P)0022      BM=SQRT(B2M)0022      XA1=XC12+B200023      XAP1=XC12+B2P0023      XAM1=XC12+B2M0027      XA2=1./(XC22+B20)0028      XAP2=1./(XC22+B2P)0028      XAM2=1./(XC22+B2M)0022      XNA=XNX2+B200023      XNAP=XNX2+B2P0023      XNAM=XNX2+B2M0020      F=B20-XC220021      FP=B2P-XC220021      FM=B2M-XC220029      XLN1=ALOG(XN21/XNA)0031      XLNP1=ALOG(XN21/XNAP)^^^^^^^^^^0031      XLNM1=ALOG(XN21/XNAM)0024      XLN2=XLN1+ADLN0026      XLNP2=XLNP1+ADLN0026      XLNM2=XLNM1+ADLN0040      ALN=0.25*(XLNP1+XLNM1-2.*XLN1)0032      S0=(ATAN(XNX/B)+HPI)/B0035      S0P=(ATAN(XNX/BP)+HPI)/BP0035      S0M=(ATAN(XNX/BM)+HPI)/BM0033      S1=(XLN1*.5+XC1*S0)/XA10037      S1P=(XLNP1*.5+XC1*S0P)/XAP10037      S1M=(XLNM1*.5+XC1*S0M)/XAM10044      S2=(XC2*XA2*XLN2-XNR-F*XA2*S0)*XA20051      S2P=(XC2*XAP2*XLNP2-XNR-FP*XAP2*S0P)*XAP20051      S2M=(XC2*XAM2*XLNM2-XNR-FM*XAM2*S0M)*XAM20038      G1=(B20*S0-0.5*XC1*XLN1)/XA10042      G1P=(B2P*S0P-0.5*XC1*XLNP1)/XAP10042      G1M=(B2M*S0M-0.5*XC1*XLNM1)/XAM10053      G2=((0.5*F*XLN2+2.*S0*B20*XC2)*XA2+XR2)*XA20059      G2P=((0.5*FP*XLNP2+2.*S0P*B2P*XC2)*XAP2+XR2)*XAP20059      G2M=((0.5*FM*XLNM2+2.*S0M*B2M*XC2)*XAM2+XR2)*XAM20047      BX=FY*(B0*(ZS*S0-0.5*(ZP*S0P+ZM*S0M))0074     * +B1*(ZS*S1-0.5*(ZP*S1P+ZM*S1M))+B2*(ZS*S2-0.5*(ZP*S2P+ZM*S2M)))0015      BY=0.0068      BZ=FY*(B0*ALN+B1*(G1-0.5*(G1P+G1M))+B2*(G2-0.5*(G2P+G2M)))0005C0072C    CALCULATION OF THE MAGNETOTAIL CURRENT CONTRIBUTION IS FINISHED0005C0025      EX1=EXP(X/DELX)0020      EX2=EX1**20017      Y2=Y**20017      Z2=Z**20016      YZ=Y*Z0072      BX=BX+(EX1*PA(1)+EX2*PA(3))*Z*CPS+(EX1*PA(2)+EX2*(PA(4)+PA(5)*0027     *Y2+PA(6)*Z2))*SPS0072      BY=(EX1*PA(7)+EX2*PA(9))*YZ*CPS+(EX1*PA(8)+EX2*(PA(10)+PA(11)*0030     *Y2+PA(12)*Z2))*Y*SPS0072      BZ=BZ+(EX1*(PA(13)+PA(14)*Y2+PA(15)*Z2)+EX2*(PA(17)+PA(18)*Y2+0072     *PA(19)*Z2))*CPS+(EX1*PA(16)+EX2*(PA(20)+PA(21)*Y2+PA(22)*Z2))*0015     *Z*SPS0005C0064C   DCF AND FAC CONTRIBUTION HAS BEEN ADDED TO BX,BY, AND BZ0005C0025      XSM=X*CPS-Z*SPS0025      ZSM=X*SPS+Z*CPS0019      Z2=ZSM**20025      RR2=XSM**2+Y**20022      RR=SQRT(RR2)0038      ZN=SQRT((RR2+Z2)/RRC2+4.)**50029      BRSM=DSTR*3.*ZSM/ZN0042      BZSM=DSTR*(2.*Z2-RR2+8.*RRC2)/ZN0020      BY1=BRSM*Y0023      BXSM=BRSM*XSM0033      BX=BX+BXSM*CPS+BZSM*SPS0033      BZ=BZ-BXSM*SPS+BZSM*CPS0019      BY=BY+BY10005C0054C   RING CURRENT FIELD HAS BEEN TAKEN INTO ACCOUNT0005C0016      RETURN^^^^^0013      END0077C------------------------------------------------------------------------0005C0051      SUBROUTINE EXT87T(IOPT,PS,X,Y,Z,BX,BY,BZ)0005C0072C  COMPUTES GSM COMPONENTS OF THE MAGNETIC FIELD OF EXTRATERRESTRIAL0070C  CURRENT SYSTEMS UP TO GEOCENTRIC RADIAL DISTANCES ABOUT 30 RE .0073C  CORRESPONDS TO MAGNETOSPHERIC FIELD MODEL (N.A.TSYGANENKO, Planet.0059C  Space Sci.,v.35, p.1347, 1987), 'TRUNCATED VERSION',0055C  BASED ON IMP-A,C,D,E,F,G,H,I,J (YEARS 1966-1980)0058C  AND HEOS-1,2 (1969-1974) SATELLITE MERGED DATA SET.0074C----- INPUT PARAMETERS:  IOPT - SPECIFIES MODEL VERSION, ACCORDING TO0056C    GROUND DISTURBANCY LEVEL, IN THE FOLLOWING WAY:0005C0074C  IOPT=   1       2       3       4        5        6        7      80041C                      CORRESPOND TO:0076C   KP=   0,0+    1-,1    1+,2-   2,2+   3-,3,3+   4-,4,4+   >=5-   >=5+0005C0045C    PS - GEODIPOLE TILT ANGLE IN RADIANS0060C    X,Y,Z - GSM COORDINATES OF THE POINT IN EARTH RADII0074C------ OUTPUT PARAMETERS:   BX,BY,BZ - GSM COMPONENTS OF THE EXTERNAL0036C    MAGNETIC FIELD IN NANOTESLA0005C0005C0053C                   AUTHOR: NIKOLAI A. TSYGANENKO0052C                           INSTITUTE OF PHYSICS0043C			    ST.-PETERSBURG STATE UNIVERSITY0033C			    STARY PETERGOF 1989040026C			    ST.-PETERSBURG0018C			    RUSSIA0005C0018	IMPLICIT NONE00040061	REAL PS,X,Y,Z,BX,BY,BZ,PSI,HPI,FC,RT,X1,C1,RRC2,DSTR,XN,0049     1     RH,DY,DELX,B0,B1,XN21,SPS,CPS,RPS,0057     2     ZS,ZP,ZM,FY,XNX,XNX2,XC1,XC12,B20,B2P,B2M,0051     3     B,BP,BM,XA1,XAP1,XAM1,XNA,XNAP,XNAM,0050     4     XLN1,XLNP1,XLNM1,ALN,S0,S0P,S0M,S1,0043     5     S1P,S1M,G1,G1P,G1M,EX,Y2,Z2,0054     6     YZ,XSM,ZSM,RR2,RR,ZN,BRSM,BZSM,BY1,BXSM00040022	INTEGER IOPT,IP,I00040071      REAL GA1(24),GA2(24),GA3(24),GA4(24),GA5(24),GA6(24),GA7(24),0026     *  GA8(24),PA(24)0075      DATA GA1/1.126,26.66,-.077,-.06102,-.06197,-2.048,.00327,.008473,0070     *12.72,-.00867,-.01953,-.3437,-.002903,-.000999,18.41,-270.3,0056     *-25.94,5.21,-6.2,2.29,11.96,8.315,44.22,11.15/^^^^^^^^^0075      DATA GA2/1.403,29.24,-.0693,-.0864,-.07202,-2.068,.00286,.007438,0073     *16.37,-.02705,-.0281,-.604,-.002256,.000152,20.2,-140.1,-29.65,0049     *5.62,-5.52,2.02,14.66,8.06,27.76,10.94/0070      DATA GA3/1.589,31.07,-.06527,-.07447,-.07632,-2.413,.002719,0075     *.01098,16.2,-.02355,-.03475,-.4377,-.002169,-.001383,18.7,-292.6,0055     *-35.25,5.29,-5.18,2.21,14.03,7.66,17.56,10.9/0069      DATA GA4/1.699,36.28,-.07514,-.1448,-.08049,-2.209,.000919,0074     *.01084,17.38,-.03516,-.03886,-1.169,.004239,.000881,21.79,-162.,0056     *-41.87,5.15,-3.62,2.35,17.26,7.61,17.99,10.74/0074      DATA GA5/2.141,41.51,-.1518,-.1857,-.1015,-2.929,.004584,.01589,0074     *18.29,-.02514,-.05927,-1.336,.00185,.001066,21.31,-358.8,-47.91,0048     *5.13,-3.74,2.07,17.23,6.33,32.51,9.73/0069      DATA GA6/2.252,39.35,-.04525,-.2062,-.1491,-3.059,-.000183,0074     *.02614,15.48,-.02144,-.06608,-1.855,.006199,-.00013,23.91,-161.,0056     *-51.48,4.61,-3.32,1.68,15.22,6.68,.6765,8.007/0073      DATA GA7/2.773,40.95,.00667,-.133,-.1304,-5.187,.004623,.03651,0073     *20.,-.03765,-.09066,.5838,-.01462,-.007189,24.87,-186.1,-74.81,0047     *4.57,-4.03,1.7,12.15,6.87,-1.746,8.9/0075      DATA GA8/2.919,34.96,2*0.,-.1609,-5.077,2*0.,22.1,-.05915,-.1051,0073     *.6321,2*0.,28.11,-330.1,-86.82,4.,-3.,1.73,12.56,5.11,4.,7.866/0067      DATA IP,PSI,HPI,FC,RT/100,10.,1.5707963,0.3183099031,30./0032      IF (IOPT.EQ.IP) GOTO 20017      IP=IOPT0021      DO 1 I=1,240035      IF (IP.EQ.1) PA(I)=GA1(I)0035      IF (IP.EQ.2) PA(I)=GA2(I)0035      IF (IP.EQ.3) PA(I)=GA3(I)0035      IF (IP.EQ.4) PA(I)=GA4(I)0035      IF (IP.EQ.5) PA(I)=GA5(I)0035      IF (IP.EQ.6) PA(I)=GA6(I)0035      IF (IP.EQ.7) PA(I)=GA7(I)0035      IF (IP.EQ.8) PA(I)=GA8(I)0018  1   CONTINUE0022      C1=PA(20)**20024      RRC2=PA(18)**20029      DSTR=PA(17)/RRC2*4.0019      XN=PA(19)0019      RH=PA(22)0019      X1=PA(23)0019      DY=PA(21)0021      DELX=PA(24)0019      B0=PA(15)0019      B1=PA(16)0025      XN21=(XN-X1)**20041  2   IF(ABS(PS-PSI).LT.1.E-6) GOTO 3^^^^^0016      PSI=PS0021      SPS=SIN(PS)0021      CPS=COS(PS)0020      RPS=RH*SPS0005C0076C   COMPUTATION BEGINS HERE IF PARAMETERS IOPT AND PS REMAINED UNCHANGED0051C   AFTER THE PRECEDING CALL OF THIS SUBROUTINE0005C0018  3   ZS=Z-RPS0017      ZP=Z-RT0017      ZM=Z+RT0030      FY=FC/(1.+(Y/DY)**2)0018      XNX=XN-X0021      XNX2=XNX**20018      XC1=X-X10021      XC12=XC1**20022      B20=ZS**2+C10022      B2P=ZP**2+C10022      B2M=ZM**2+C10021      B=SQRT(B20)0022      BP=SQRT(B2P)0022      BM=SQRT(B2M)0022      XA1=XC12+B200023      XAP1=XC12+B2P0023      XAM1=XC12+B2M0022      XNA=XNX2+B200023      XNAP=XNX2+B2P0023      XNAM=XNX2+B2M0029      XLN1=ALOG(XN21/XNA)0031      XLNP1=ALOG(XN21/XNAP)0031      XLNM1=ALOG(XN21/XNAM)0040      ALN=0.25*(XLNP1+XLNM1-2.*XLN1)0032      S0=(ATAN(XNX/B)+HPI)/B0035      S0P=(ATAN(XNX/BP)+HPI)/BP0035      S0M=(ATAN(XNX/BM)+HPI)/BM0033      S1=(XLN1*.5+XC1*S0)/XA10037      S1P=(XLNP1*.5+XC1*S0P)/XAP10037      S1M=(XLNM1*.5+XC1*S0M)/XAM10038      G1=(B20*S0-0.5*XC1*XLN1)/XA10042      G1P=(B2P*S0P-0.5*XC1*XLNP1)/XAP10042      G1M=(B2M*S0M-0.5*XC1*XLNM1)/XAM10047      BX=FY*(B0*(ZS*S0-0.5*(ZP*S0P+ZM*S0M))0043     * +B1*(ZS*S1-0.5*(ZP*S1P+ZM*S1M)))0015      BY=0.0046      BZ=FY*(B0*ALN+B1*(G1-0.5*(G1P+G1M)))0005C0072C    CALCULATION OF THE MAGNETOTAIL CURRENT CONTRIBUTION IS FINISHED0005C0024      EX=EXP(X/DELX)0017      Y2=Y**20017      Z2=Z**20016      YZ=Y*Z0062      BX=BX+EX*(CPS*PA(1)*Z+SPS*(PA(2)+PA(3)*Y2+PA(4)*Z2))0062      BY=EX*(CPS*PA(5)*YZ+SPS*Y*(PA(6)+PA(7)*Y2+PA(8)*Z2))0065      BZ=BZ+EX*(CPS*(PA(9)+PA(10)*Y2+PA(11)*Z2)+SPS*Z*(PA(12)0034     *  +PA(13)*Y2+PA(14)*Z2))0005C0064C   DCF AND FAC CONTRIBUTION HAS BEEN ADDED TO BX,BY, AND BZ0005C0025      XSM=X*CPS-Z*SPS0025      ZSM=X*SPS+Z*CPS0019      Z2=ZSM**20025      RR2=XSM**2+Y**20022      RR=SQRT(RR2)0038      ZN=SQRT((RR2+Z2)/RRC2+4.)**50029      BRSM=DSTR*3.*ZSM/ZN0042      BZSM=DSTR*(2.*Z2-RR2+8.*RRC2)/ZN0020      BY1=BRSM*Y0023      BXSM=BRSM*XSM0033      BX=BX+BXSM*CPS+BZSM*SPS0033      BZ=BZ-BXSM*SPS+BZSM*CPS^^^^^^0019      BY=BY+BY10005C0054C   RING CURRENT FIELD HAS BEEN TAKEN INTO ACCOUNT0005C0016      RETURN0013      END0005c0005c0005c0055 	      	SUBROUTINE EXT89KP(IOPT,PS,X,Y,Z,BX,BY,BZ)0072C   COMPUTES GSM COMPONENTS OF THE MAGNETIC FIELD PRODUCED BY EXTRA-0073C   TERRESTRIAL CURRENT SYSTEMS IN THE GEOMAGNETOSPHERE. THE MODEL IS0074C   VALID UP TO GEOCENTRIC DISTANCES OF 70 RE AND IS BASED ON THE MER-0075C   GED IMP-A,C,D,E,F,G,H,I,J (1966-1974) AND HEOS-1 AND -2 (1969-1974)0074C   SPACECRAFT DATA SET.  REFERENCE: N.A. TSYGANENKO, A MAGNETOSPHERIC0075C   MAGNETIC FIELD MODEL WITH A WARPED TAIL CURRENT SHEET: PLANET.SPACE0034C   SCI., V.37, PP.5-20, 1989.0005C0073C----INPUT PARAMETERS: IOPT - SPECIFIES THE GROUND DISTURBANCE LEVEL:0005C0070C   IOPT= 1         2         3          4           5           60037C                  CORRESPOND TO:0072C    KP= 0,0+    1-,1,1+    2-,2,2+    3-,3,3+     4-,4,4+      >=5-0005C0045C    PS - GEODIPOLE TILT ANGLE IN RADIANS0063C    X, Y, Z  - GSM COORDINATES OF THE POINT IN EARTH RADII0005C0075C----OUTPUT PARAMETERS: BX,BY,BZ - GSM COMPONENTS OF THE MODEL MAGNETIC0048C                        FIELD IN NANOTESLAS0005C0005C0052C              AUTHOR:     NIKOLAI A. TSYGANENKO0078C                          INSTITUTE OF PHYSICS, ST.-PETERSBURG UNIVERSITY0032C			   STARY PETERGOF 1989040025C			   ST.-PETERSBURG0017C			   RUSSIA0005C0018	IMPLICIT NONE00040063	REAL PS,X,Y,Z,BX,BY,BZ,PSI,DELX,ADR,D0,DD,RC,G,A,DY,SX,HA,0072     1       HA02,RDYC2,HLWC2M,DRDYCM,HLDXM,DDEL,RDY2,DRDY2M,HXLW2M,0069     2       DADD05,DADD18,SPS,CPS,HTPS,GSPM,XSM,ZSM,X2SM,Y2,RO2,0072     3       XXD,XXD2L2,RSQXDL,H,HSX,XSIXT,XSIXTD,SXSIX,DDOP,DDOPDX,0070     4       D,DDX,DDY,XRC,SXRC16,Y4,Y410,HY,HYS,ZS,DZSX,DZSY,XSX,0071     5       RQ,SRQ,FY,W,DWX,DWY,ZR,T,AT,S1,F5,F7,F3,F9,XWYW,FR,FS,0072     6       WT,WTFS,BRZR1,BRZR2,BXT,BYT,BZT,RX2A2,SRX2A2,FDR,DFDRX,0070     7       DDR,TDR,ADRT,ADRT2,ADT2R2,FK1,FK2,BXDR,EX,Z2,YZ,BXCF,0067     8       BYCF,BZCF,FCY,XSXC,RQC2,SRQC2,WC,DWCX,DWCY,XWCYWC,^^^^^^^^^^^^^^^^^^^0066     9       ZP,ZM,SP,SM,WCSP,WCSM,RSPRT,RSMRT,FXP,FXM,FYP,FYM00040065	REAL FZP,FZM,AA4SPS,BXC,BYC,BZC,A02,DYC,XLWC2,XLD2,DEL,XLW2,0038     1       DADD,XD,GAM,F1,SXC,RT00040022	INTEGER IOPT,IP,I00040069      REAL GA1(28),GA2(28),GA3(28),GA4(28),GA5(28),GA6(28),PA(28)0074      DATA GA1/-98.72,-10014.,15.03,76.62,-10237.,1.813,31.1,-0.07464,0071     *   -0.07764,0.003303,-1.129,0.001663,0.000988,18.21,-0.03018,0067     *   4*0.,24.74,8.16,2.08,-0.88,9.08,3.84,13.55,26.94,5.75/0074      DATA GA2/-35.65,-12800.,14.37,124.5,-13543.,2.316,35.64,-0.0741,0066     * -0.1081,0.003924,-1.451,0.00202,0.00111,21.37,-0.04567,0066     * 4*0.,22.33,8.12,1.664,0.932,9.24,2.43,13.81,28.83,6.05/0075      DATA GA3/-77.45,-14588.,64.85,123.9,-16229.,2.641,42.46,-0.07611,0067     * -0.1579,0.004078,-1.391,0.00153,0.000727,21.86,-0.04199,0063     * 4*0.,20.9,6.28,1.54,4.18,9.61,6.59,15.08,30.57,7.43/0073      DATA GA4/-70.12,-16125.,90.71,38.08,-19630.,3.181,47.5,-0.1327,0067     * -0.1864,0.01382,-1.488,0.002962,0.000897,22.74,-0.04095,0065     * 4*0.,18.64,6.27,0.935,5.39,8.57,5.94,15.63,31.47,8.10/0073      DATA GA5/-162.5,-15806.,160.6,5.888,-27534.,3.607,51.1,-0.1006,0067     * -0.1927,0.03353,-1.392,0.001594,0.002439,22.41,-0.04925,0064     * 4*0.,18.31,6.2,0.768,5.07,10.06,6.67,16.1,30.04,8.26/0073      DATA GA6/-128.4,-16184.,149.1,215.5,-36435.,4.09,49.09,-0.0231,0067     * -0.1359,0.01989,-2.298,0.004911,0.003421,21.79,-0.05447,0066     * 4*0.,19.48,5.83,0.332,6.47,10.47,9.08,15.85,25.27,7.98/0061      DATA DEL,GAM,DYC,XD,XLD2,XLW2,A02,RT,SXC,XLWC2,DADD0054     * /0.01,4.,20.,0.,40.,170.,25.,30.,4.,50.,1./0030      DATA IP,PSI/100,10./0027	IF (IOPT.EQ.IP) GOTO 20005C0077C   INITIAL CALCULATION OF CONSTANTS FROM A GIVEN SET OF MODEL PARAMETERS0045C   SPECIFIED BY THE OPTION NUMBER (IOPT)0005C0012	IP=IOPT0016	DO 1 I=1,280030	IF (IP.EQ.1) PA(I)=GA1(I)0037        IF (IP.EQ.2) PA(I)=GA2(I)0030	IF (IP.EQ.3) PA(I)=GA3(I)0030	IF (IP.EQ.4) PA(I)=GA4(I)0030	IF (IP.EQ.5) PA(I)=GA5(I)0030	IF (IP.EQ.6) PA(I)=GA6(I)^^^^^^^^^^^^^^^^^^^0020   1    CONTINUE0016	DELX=PA(20)0015	ADR=PA(21)0014	D0=PA(22)0014	DD=PA(23)0014	RC=PA(24)0013	G=PA(25)0013	A=PA(26)0014	DY=PA(27)0014	SX=PA(28)0013	HA=0.5*A0017	HA02=0.5*A020020	RDYC2=1./DYC**20029        HLWC2M=-0.5*XLWC20021	DRDYCM=-2.*RDYC20020	HLDXM=-0.5*XLD20016	DDEL=2.*DEL0018	RDY2=1./DY**20020	DRDY2M=-2.*RDY20021	HXLW2M=-0.5*XLW20020	DADD05=DADD*0.50021	DADD18=-18.*DADD0005C0062C   COEFFICIENTS PA(16)-PA(19) ARE FOUND FROM PA(6)-PA(13)0052C   SO THAT THE MAGNETIC FIELD BE DIVERGENCELESS0005C0036	PA(16)=-0.5*(PA(6)/DELX+PA(10))0032	PA(17)=-(PA(7)/DELX+PA(11))0035	PA(18)=-(PA(8)/DELX+3.*PA(12))0035	PA(19)=-(PA(9)/DELX+PA(13))/3.0005C0011	GOTO 30043   2    IF(ABS(PS-PSI).LT.1.E-6) GOTO 40018   3    PSI=PS0016	SPS=SIN(PS)0016	CPS=COS(PS)0022	HTPS=SPS/(2.*CPS)0016	GSPM=-G*SPS0005C0059C   COMPUTATION BEGINS HERE IF IOPT AND PS ARE THE SAME0005C0027   4    XSM=X*CPS-Z*SPS0020	ZSM=Z*CPS+X*SPS0016	X2SM=XSM**20013  	Y2=Y*Y0016	RO2=X2SM+Y20015	XXD=XSM-XD0028	XXD2L2=1./(XXD**2+XLD2)0024	RSQXDL=SQRT(XXD2L2)0026	H=0.5*(1.+XXD*RSQXDL)0029	HSX=-HLDXM*XXD2L2*RSQXDL0018	XSIXT=XSM+16.0029	XSIXTD=1./(XSIXT**2+36.)0023	SXSIX=SQRT(XSIXTD)0033	DDOP=DADD05*(1.-XSIXT*SXSIX)0031	DDOPDX=DADD18*XSIXTD*SXSIX0027	D=D0+DEL*Y2+GAM*H+DDOP0023	DDX=GAM*HSX+DDOPDX0015	DDY=DDEL*Y0015	XRC=XSM+RC0028	SXRC16=SQRT(XRC**2+16.)0013	Y4=Y2*Y20022	Y410=1./(Y4+1.E4)0017	HY=Y2*Y410*Y0021	HYS=HY*Y410*4.E40012	HY=HY*Y0033	ZS=HTPS*(XRC-SXRC16)+GSPM*HY0030	DZSX=HTPS*(1.-XRC/SXRC16)0018	DZSY=GSPM*HYS0005C0067C  ZS=ZS(XSM,YSM) DEFINES THE SHAPE OF THE WARPED CURRENT SHEET0005C0015	XSX=XSM-SX0024	RQ=1./(XSX**2+XLW2)0017	SRQ=SQRT(RQ)0023	FY=1./(1.+Y2*RDY2)0026	W=0.5*(1.-XSX*SRQ)*FY0025	DWX=HXLW2M*RQ*SRQ*FY0022	DWY=DRDY2M*W*Y*FY0014	ZR=ZSM-ZS0023	T=SQRT(ZR**2+D**2)0011	AT=A+T0023	S1=SQRT(AT**2+RO2)0013	F5=1./S10018	F7=1./(S1+AT)0013	F1=F5*F70013	F3=F5**30013	F9=AT*F30023	XWYW=XSM*DWX+Y*DWY0028	FR=ZR*(XSM*DZSX+Y*DZSY)0028	FS=FR-D*(XSM*DDX+Y*DDY)0011	WT=W/T0015	WTFS=WT*FS0016	BRZR1=WT*F10016	BRZR2=WT*F30037	BXT=(PA(1)*BRZR1+PA(2)*BRZR2)*ZR0014	BYT=BXT*Y^^^^^^^0016	BXT=BXT*XSM0066	BZT=PA(1)*(W*F5+XWYW*F7+WTFS*F1)+PA(2)*(W*F9+XWYW*F1+WTFS*F3)0005C0066C   CONTRIBUTION FROM CENTRAL TAIL CURRENT SHEET (BXT,BYT,BZT)0064C   IS FOUND.  NOW LET US PROCEED TO THE RING CURRENT FIELD:0005C0024	RX2A2=1./(X2SM+A02)0023	SRX2A2=SQRT(RX2A2)0028	FDR=0.5*(1.+XSM*SRX2A2)0028	DFDRX=HA02*RX2A2*SRX2A20023	DDR=D0+DD*FDR+DDOP0027	TDR=SQRT(ZR**2+DDR**2)0017	ADRT=ADR+TDR0018	ADRT2=ADRT**20026	ADT2R2=1./(ADRT2+RO2)0031	FK1=ADT2R2**2*SQRT(ADT2R2)0024	FK2=3.*ADRT*FK1/TDR0022	BXDR=PA(5)*ZR*FK20019	BYT=BYT+BXDR*Y0021	BXT=BXT+BXDR*XSM0060	BZT=BZT+PA(5)*((2.*ADRT2-RO2)*FK1+FK2*(FR-DDR*(DD*DFDRX0026     *  +DDOPDX)*XSM))0005C0043C   NOW CALCULATE CHAPMAN-FERRARO FIELD0055C   PLUS AVERAGE FIELD-ALIGNED CURRENT CONTRIBUTION0005C0019	EX=EXP(X/DELX)0011	Z2=Z*Z0011	YZ=Y*Z0056	BXCF=EX*(CPS*PA(6)*Z+SPS*(PA(7)+PA(8)*Y2+PA(9)*Z2))0063	BYCF=EX*(CPS*PA(10)*YZ+SPS*Y*(PA(11)+PA(12)*Y2+PA(13)*Z2))0068        BZCF=EX*(CPS*(PA(14)+PA(15)*Y2+PA(16)*Z2)+SPS*Z*(PA(17)+0035     *    PA(18)*Y2+PA(19)*Z2))0005C0066C   MAGNETOTAIL RETURN CURRENT FIELD COMPONENTS (BXC,BYC,BZC):0005C0025	FCY=1./(1.+Y2*RDYC2)0015	XSXC=X-SXC0028	RQC2=1./(XSXC**2+XLWC2)0021	SRQC2=SQRT(RQC2)0031	WC=0.5*(1.-XSXC*SRQC2)*FCY0031	DWCX=HLWC2M*RQC2*SRQC2*FCY0025	DWCY=DRDYCM*WC*Y*FCY0025	XWCYWC=X*DWCX+Y*DWCY0016	RO2=Y2+X**20012	ZP=Z+RT0012	ZM=Z-RT0023	SP=SQRT(ZP**2+RO2)0023	SM=SQRT(ZM**2+RO2)0015	WCSP=WC/SP0015	WCSM=WC/SM0021	RSPRT=1./(SP+ZP)0021	RSMRT=1./(SM-ZM)0019	FXP=WCSP*RSPRT0020	FXM=-WCSM*RSMRT0014	FYP=FXP*Y0014	FYM=FXM*Y0014	FXP=FXP*X0014	FXM=FXM*X0026	FZP=WCSP+XWCYWC*RSPRT0026	FZM=WCSM+XWCYWC*RSMRT0021	AA4SPS=PA(4)*SPS0041	BXC=PA(3)*(FXP+FXM)+(FXP-FXM)*AA4SPS0041	BYC=PA(3)*(FYP+FYM)+(FYP-FYM)*AA4SPS0041	BZC=PA(3)*(FZP+FZM)+(FZP-FZM)*AA4SPS0005C0062C   AT LAST, SUM UP THE FIELDS.  THE CENTRAL CURRENT SHEET0061C   FIELD COMPONENTS ARE TRANSFORMED INTO GSM COORDINATES0005C0032	BX=BXC+BXT*CPS+BZT*SPS+BXCF0020	BY=BYC+BYT+BYCF0032	BZ=BZC+BZT*CPS-BXT*SPS+BZCF0011	RETURN0008	END^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^0077c------------------------------------------------------------------------0005c0055 	      	SUBROUTINE EXT89AE(IOPT,PS,X,Y,Z,BX,BY,BZ)0005C0072C   COMPUTES GSM COMPONENTS OF THE MAGNETIC FIELD PRODUCED BY EXTRA-0073C   TERRESTRIAL CURRENT SYSTEMS IN THE GEOMAGNETOSPHERE. THE MODEL IS0074C   VALID UP TO GEOCENTRIC DISTANCES OF 70 RE AND IS BASED ON THE MER-0075C   GED IMP-A,C,D,E,F,G,H,I,J (1966-1974) AND HEOS-1 AND -2 (1969-1974)0074C   SPACECRAFT DATA SET.  REFERENCE: N.A. TSYGANENKO, A MAGNETOSPHERIC0075C   MAGNETIC FIELD MODEL WITH A WARPED TAIL CURRENT SHEET: PLANET.SPACE0034C   SCI., V.37, PP.5-20, 1989.0005C0073C----INPUT PARAMETERS: IOPT - SPECIFIES THE GROUND DISTURBANCE LEVEL:0005C0070C   IOPT= 1         2         3          4           5           60037C                  CORRESPOND TO:0072C    AE=0-50     50-100    100-150    150-250     250-400      >=4000005C0045C    PS - GEODIPOLE TILT ANGLE IN RADIANS0063C    X, Y, Z  - GSM COORDINATES OF THE POINT IN EARTH RADII0005C0075C----OUTPUT PARAMETERS: BX,BY,BZ - GSM COMPONENTS OF THE MODEL MAGNETIC0048C                        FIELD IN NANOTESLAS0005C0052C              AUTHOR:     NIKOLAI A. TSYGANENKO0078C                          INSTITUTE OF PHYSICS, ST.-PETERSBURG UNIVERSITY0032C			   STARY PETERGOF 1989040025C			   ST.-PETERSBURG0017C			   RUSSIA0005C0005C0018	IMPLICIT NONE00040063	REAL PS,X,Y,Z,BX,BY,BZ,PSI,DELX,ADR,D0,DD,RC,G,A,DY,SX,HA,0072     1       HA02,RDYC2,HLWC2M,DRDYCM,HLDXM,DDEL,RDY2,DRDY2M,HXLW2M,0069     2       DADD05,DADD18,SPS,CPS,HTPS,GSPM,XSM,ZSM,X2SM,Y2,RO2,0072     3       XXD,XXD2L2,RSQXDL,H,HSX,XSIXT,XSIXTD,SXSIX,DDOP,DDOPDX,0070     4       D,DDX,DDY,XRC,SXRC16,Y4,Y410,HY,HYS,ZS,DZSX,DZSY,XSX,0071     5       RQ,SRQ,FY,W,DWX,DWY,ZR,T,AT,S1,F5,F7,F3,F9,XWYW,FR,FS,0072     6       WT,WTFS,BRZR1,BRZR2,BXT,BYT,BZT,RX2A2,SRX2A2,FDR,DFDRX,0070     7       DDR,TDR,ADRT,ADRT2,ADT2R2,FK1,FK2,BXDR,EX,Z2,YZ,BXCF,0067     8       BYCF,BZCF,FCY,XSXC,RQC2,SRQC2,WC,DWCX,DWCY,XWCYWC,^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^0066     9       ZP,ZM,SP,SM,WCSP,WCSM,RSPRT,RSMRT,FXP,FXM,FYP,FYM00040065	REAL FZP,FZM,AA4SPS,BXC,BYC,BZC,A02,DYC,XLWC2,XLD2,DEL,XLW2,0038     1       DADD,XD,GAM,F1,SXC,RT00040022	INTEGER IOPT,IP,I00040069      REAL GA1(28),GA2(28),GA3(28),GA4(28),GA5(28),GA6(28),PA(28)0074      DATA GA1/-77.38,-10921.,-2.89,193.9,-10954.,2.048,28.83,-0.0367,0068     *   -0.0625,0.00104,-1.246,0.00166,0.000792,20.81,-0.03608,0072     *   4*0.,25.82,7.656,2.106,-0.30,8.753,2.733,13.17,27.09,5.184/0073      DATA GA2/-59.2,-13647.,11.6,237.2,-14465.,2.326,36.54,-0.07084,0068     * -0.1182,0.000146,-1.626,0.002466,0.001355,23.49,-0.04602,0065     * 4*0.,23.15,7.6,1.6,1.457,8.373,3.883,14.18,28.79,5.97/0075      DATA GA3/-65.91,-15267.,64.48,230.6,-14870.,2.712,43.28,-0.09687,0064     * -0.1671,0.009814,-1.835,0.0026,0.0017,22.73,-0.05049,0069     * 4*0.,20.76,6.314,1.42,4.141,9.436,6.266,15.07,30.35,6.852/0074      DATA GA4/-57.97,-16106.,86.28,199.7,-16796.,3.033,46.80,-0.1057,0068     * -0.1992,0.01509,-1.534,0.001783,-0.000542,22.29,-0.04664,0068     * 4*0.,19.79,5.788,1.177,5.007,9.60,7.32,15.78,33.53,7.401/0074      DATA GA5/-118.1,-16473.,136.3,158.7,-21134.,3.135,49.78,-0.1088,0068     * -0.2245,0.01759,-1.874,0.003243,-0.000268,21.97,-0.04913,0070     * 4*0.,18.48,6.25,0.8643,5.821,9.552,6.051,16.14,31.27,9.062/0074      DATA GA6/-223.7,-15054.,219.1,83.84,-31140.,3.777,51.08,-0.1261,0067     * -0.2393,0.02507,-1.426,0.001678,0.002039,19.10,-0.05344,0070     * 4*0.,16.73,6.532,0.382,5.807,8.099,6.726,16.21,25.63,9.431/0061      DATA DEL,GAM,DYC,XD,XLD2,XLW2,A02,RT,SXC,XLWC2,DADD0054     * /0.01,4.,20.,0.,40.,170.,25.,30.,4.,50.,1./0030      DATA IP,PSI/100,10./0027	IF (IOPT.EQ.IP) GOTO 20005C0077C   INITIAL CALCULATION OF CONSTANTS FROM A GIVEN SET OF MODEL PARAMETERS0045C   SPECIFIED BY THE OPTION NUMBER (IOPT)0005C0012	IP=IOPT0016	DO 1 I=1,280030	IF (IP.EQ.1) PA(I)=GA1(I)0037        IF (IP.EQ.2) PA(I)=GA2(I)0030	IF (IP.EQ.3) PA(I)=GA3(I)0030	IF (IP.EQ.4) PA(I)=GA4(I)0030	IF (IP.EQ.5) PA(I)=GA5(I)^^^^^^^^^^^^^^^^^^^^^^^^^^0030	IF (IP.EQ.6) PA(I)=GA6(I)0020   1    CONTINUE0016	DELX=PA(20)0015	ADR=PA(21)0014	D0=PA(22)0014	DD=PA(23)0014	RC=PA(24)0013	G=PA(25)0013	A=PA(26)0014	DY=PA(27)0014	SX=PA(28)0013	HA=0.5*A0017	HA02=0.5*A020020	RDYC2=1./DYC**20029        HLWC2M=-0.5*XLWC20021	DRDYCM=-2.*RDYC20020	HLDXM=-0.5*XLD20016	DDEL=2.*DEL0018	RDY2=1./DY**20020	DRDY2M=-2.*RDY20021	HXLW2M=-0.5*XLW20020	DADD05=DADD*0.50021	DADD18=-18.*DADD0005C0062C   COEFFICIENTS PA(16)-PA(19) ARE FOUND FROM PA(6)-PA(13)0052C   SO THAT THE MAGNETIC FIELD BE DIVERGENCELESS0005C0036	PA(16)=-0.5*(PA(6)/DELX+PA(10))0032	PA(17)=-(PA(7)/DELX+PA(11))0035	PA(18)=-(PA(8)/DELX+3.*PA(12))0035	PA(19)=-(PA(9)/DELX+PA(13))/3.0005C0011	GOTO 30043   2    IF(ABS(PS-PSI).LT.1.E-6) GOTO 40018   3    PSI=PS0016	SPS=SIN(PS)0016	CPS=COS(PS)0022	HTPS=SPS/(2.*CPS)0016	GSPM=-G*SPS0005C0059C   COMPUTATION BEGINS HERE IF IOPT AND PS ARE THE SAME0005C0027   4    XSM=X*CPS-Z*SPS0020	ZSM=Z*CPS+X*SPS0016	X2SM=XSM**20013  	Y2=Y*Y0016	RO2=X2SM+Y20015	XXD=XSM-XD0028	XXD2L2=1./(XXD**2+XLD2)0024	RSQXDL=SQRT(XXD2L2)0026	H=0.5*(1.+XXD*RSQXDL)0029	HSX=-HLDXM*XXD2L2*RSQXDL0018	XSIXT=XSM+16.0029	XSIXTD=1./(XSIXT**2+36.)0023	SXSIX=SQRT(XSIXTD)0033	DDOP=DADD05*(1.-XSIXT*SXSIX)0031	DDOPDX=DADD18*XSIXTD*SXSIX0027	D=D0+DEL*Y2+GAM*H+DDOP0023	DDX=GAM*HSX+DDOPDX0015	DDY=DDEL*Y0015	XRC=XSM+RC0028	SXRC16=SQRT(XRC**2+16.)0013	Y4=Y2*Y20022	Y410=1./(Y4+1.E4)0017	HY=Y2*Y410*Y0021	HYS=HY*Y410*4.E40012	HY=HY*Y0033	ZS=HTPS*(XRC-SXRC16)+GSPM*HY0030	DZSX=HTPS*(1.-XRC/SXRC16)0018	DZSY=GSPM*HYS0005C0067C  ZS=ZS(XSM,YSM) DEFINES THE SHAPE OF THE WARPED CURRENT SHEET0005C0015	XSX=XSM-SX0024	RQ=1./(XSX**2+XLW2)0017	SRQ=SQRT(RQ)0023	FY=1./(1.+Y2*RDY2)0026	W=0.5*(1.-XSX*SRQ)*FY0025	DWX=HXLW2M*RQ*SRQ*FY0022	DWY=DRDY2M*W*Y*FY0014	ZR=ZSM-ZS0023	T=SQRT(ZR**2+D**2)0011	AT=A+T0023	S1=SQRT(AT**2+RO2)0013	F5=1./S10018	F7=1./(S1+AT)0013	F1=F5*F70013	F3=F5**30013	F9=AT*F30023	XWYW=XSM*DWX+Y*DWY0028	FR=ZR*(XSM*DZSX+Y*DZSY)0028	FS=FR-D*(XSM*DDX+Y*DDY)0011	WT=W/T0015	WTFS=WT*FS0016	BRZR1=WT*F10016	BRZR2=WT*F3^^^^^^^^^^^^^^^^^^^^^^^^^^^^0037	BXT=(PA(1)*BRZR1+PA(2)*BRZR2)*ZR0014	BYT=BXT*Y0016	BXT=BXT*XSM0066	BZT=PA(1)*(W*F5+XWYW*F7+WTFS*F1)+PA(2)*(W*F9+XWYW*F1+WTFS*F3)0005C0066C   CONTRIBUTION FROM CENTRAL TAIL CURRENT SHEET (BXT,BYT,BZT)0064C   IS FOUND.  NOW LET US PROCEED TO THE RING CURRENT FIELD:0005C0024	RX2A2=1./(X2SM+A02)0023	SRX2A2=SQRT(RX2A2)0028	FDR=0.5*(1.+XSM*SRX2A2)0028	DFDRX=HA02*RX2A2*SRX2A20023	DDR=D0+DD*FDR+DDOP0027	TDR=SQRT(ZR**2+DDR**2)0017	ADRT=ADR+TDR0018	ADRT2=ADRT**20026	ADT2R2=1./(ADRT2+RO2)0031	FK1=ADT2R2**2*SQRT(ADT2R2)0024	FK2=3.*ADRT*FK1/TDR0022	BXDR=PA(5)*ZR*FK20019	BYT=BYT+BXDR*Y0021	BXT=BXT+BXDR*XSM0060	BZT=BZT+PA(5)*((2.*ADRT2-RO2)*FK1+FK2*(FR-DDR*(DD*DFDRX0026     *  +DDOPDX)*XSM))0005C0043C   NOW CALCULATE CHAPMAN-FERRARO FIELD0055C   PLUS AVERAGE FIELD-ALIGNED CURRENT CONTRIBUTION0005C0019	EX=EXP(X/DELX)0011	Z2=Z*Z0011	YZ=Y*Z0056	BXCF=EX*(CPS*PA(6)*Z+SPS*(PA(7)+PA(8)*Y2+PA(9)*Z2))0063	BYCF=EX*(CPS*PA(10)*YZ+SPS*Y*(PA(11)+PA(12)*Y2+PA(13)*Z2))0068        BZCF=EX*(CPS*(PA(14)+PA(15)*Y2+PA(16)*Z2)+SPS*Z*(PA(17)+0035     *    PA(18)*Y2+PA(19)*Z2))0005C0066C   MAGNETOTAIL RETURN CURRENT FIELD COMPONENTS (BXC,BYC,BZC):0005C0025	FCY=1./(1.+Y2*RDYC2)0015	XSXC=X-SXC0028	RQC2=1./(XSXC**2+XLWC2)0021	SRQC2=SQRT(RQC2)0031	WC=0.5*(1.-XSXC*SRQC2)*FCY0031	DWCX=HLWC2M*RQC2*SRQC2*FCY0025	DWCY=DRDYCM*WC*Y*FCY0025	XWCYWC=X*DWCX+Y*DWCY0016	RO2=Y2+X**20012	ZP=Z+RT0012	ZM=Z-RT0023	SP=SQRT(ZP**2+RO2)0023	SM=SQRT(ZM**2+RO2)0015	WCSP=WC/SP0015	WCSM=WC/SM0021	RSPRT=1./(SP+ZP)0021	RSMRT=1./(SM-ZM)0019	FXP=WCSP*RSPRT0020	FXM=-WCSM*RSMRT0014	FYP=FXP*Y0014	FYM=FXM*Y0014	FXP=FXP*X0014	FXM=FXM*X0026	FZP=WCSP+XWCYWC*RSPRT0026	FZM=WCSM+XWCYWC*RSMRT0021	AA4SPS=PA(4)*SPS0041	BXC=PA(3)*(FXP+FXM)+(FXP-FXM)*AA4SPS0041	BYC=PA(3)*(FYP+FYM)+(FYP-FYM)*AA4SPS0041	BZC=PA(3)*(FZP+FZM)+(FZP-FZM)*AA4SPS0005C0062C   AT LAST, SUM UP THE FIELDS.  THE CENTRAL CURRENT SHEET0061C   FIELD COMPONENTS ARE TRANSFORMED INTO GSM COORDINATES0005C0032	BX=BXC+BXT*CPS+BZT*SPS+BXCF0020	BY=BYC+BYT+BYCF0032	BZ=BZC+BZT*CPS-BXT*SPS+BZCF0011	RETURN0008	END^^^^^^^^^^^^^^^^^^^^^^^^^0073c--------------------------------------------------------------------0005c0070      SUBROUTINE SUN(IYR,IDAY,IHOUR,MIN,ISEC,GST,SLONG,SRASN,SDEC)0005C0074C  CALCULATES FOUR QUANTITIES NECESSARY FOR COORDINATE TRANSFORMATIONS0078C  WHICH DEPEND ON SUN POSITION (AND, HENCE, ON UNIVERSAL TIME AND SEASON)0005C0031C-------  INPUT PARAMETERS:0082C  IYR,IDAY,IHOUR,MIN,ISEC -  YEAR, DAY, AND UNIVERSAL TIME IN HOURS, MINUTES,0056C    AND SECONDS  (IDAY=1 CORRESPONDS TO JANUARY 1).0005C0032C-------  OUTPUT PARAMETERS:0075C  GST - GREENWICH MEAN SIDEREAL TIME, SLONG - LONGITUDE ALONG ECLIPTIC0073C  SRASN - RIGHT ASCENSION,  SDEC - DECLINATION  OF THE SUN (RADIANS)0074C  THIS SUBROUTINE HAS BEEN COMPILED FROM: RUSSELL C.T., COSM.ELECTRO-0034C  DYN., 1971, V.2,PP.184-196.0005C0005C0047C                   AUTHOR: Gilbert D. Mead0005C0005C0018	IMPLICIT NONE00040061	REAL GST,SLONG,SRASN,SDEC,RAD,T,VL,G,OBLIQ,SOB,SLP,SIND,0024     1       COSD,SC00040036	INTEGER IYR,IDAY,IHOUR,MIN,ISEC00040034      DOUBLE PRECISION DJ,FDAY0032      DATA RAD/57.295779513/0047      IF(IYR.LT.1901.OR.IYR.GT.2099) RETURN0054      FDAY=DFLOAT(IHOUR*3600+MIN*60+ISEC)/86400.D00056      DJ=365*(IYR-1900)+(IYR-1901)/4+IDAY-0.5D0+FDAY0021      T=DJ/36525.0052      VL=DMOD(279.696678+0.9856473354*DJ,360.D0)0071      GST=DMOD(279.690983+.9856473354*DJ+360.*FDAY+180.,360.D0)/RAD0054      G=DMOD(358.475845+0.985600267*DJ,360.D0)/RAD0071      SLONG=(VL+(1.91946-0.004789*T)*SIN(G)+0.020094*SIN(2.*G))/RAD0054      IF(SLONG.GT.6.2831853) SLONG=SLONG-6.28318530048      IF (SLONG.LT.0.) SLONG=SLONG+6.28318530042      OBLIQ=(23.45229-0.0130125*T)/RAD0024      SOB=SIN(OBLIQ)0028      SLP=SLONG-9.924E-50005C0076C   THE LAST CONSTANT IS A CORRECTION FOR THE ANGULAR ABERRATION  DUE TO0039C   THE ORBITAL MOTION OF THE EARTH0005C0027      SIND=SOB*SIN(SLP)0031      COSD=SQRT(1.-SIND**2)0022      SC=SIND/COSD0023      SDEC=ATAN(SC)0067      SRASN=3.141592654-ATAN2(COS(OBLIQ)/SOB*SC,-COS(SLP)/COSD)0016      RETURN0013      END^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^0076c-----------------------------------------------------------------------0005c0047      SUBROUTINE SPHCAR(R,TETA,PHI,X,Y,Z,J)0005C0068C   CONVERTS SPHERICAL COORDS INTO CARTESIAN ONES AND VICA VERSA0035C    (TETA AND PHI IN RADIANS).0005C0041C                  J>0            J<00043C-----INPUT:   J,R,TETA,PHI     J,X,Y,Z0045C----OUTPUT:      X,Y,Z        R,TETA,PHI0005C0005C0053C                   AUTHOR: NIKOLAI A. TSYGANENKO0052C                           INSTITUTE OF PHYSICS0063C                           ST.-PETERSBURG STATE UNIVERSITY0053C                           STARY PETERGOF 1989040046C                           ST.-PETERSBURG0038C                           RUSSIA0005C0018	IMPLICIT NONE00040029	REAL R,TETA,PHI,X,Y,Z,SQ00040014	INTEGER J00040027      IF(J.GT.0) GOTO 30022      SQ=X**2+Y**20025      R=SQRT(SQ+Z**2)0030      IF (SQ.NE.0.) GOTO 20016      PHI=0.0029      IF (Z.LT.0.) GOTO 10017      TETA=0.0016      RETURN0026  1   TETA=3.1415926540016      RETURN0021  2   SQ=SQRT(SQ)0024      PHI=ATAN2(Y,X)0026      TETA=ATAN2(SQ,Z)0043      IF (PHI.LT.0.) PHI=PHI+6.283185310016      RETURN0024  3   SQ=R*SIN(TETA)0023      X=SQ*COS(PHI)0023      Y=SQ*SIN(PHI)0023      Z=R*COS(TETA)0016      RETURN0013      END0073c--------------------------------------------------------------------0005c0059      SUBROUTINE BSPCAR(TETA,PHI,BR,BTET,BPHI,BX,BY,BZ)0065C   CALCULATES CARTESIAN FIELD COMPONENTS FROM SPHERICAL ONES0070C-----INPUT:   TETA,PHI - SPHERICAL ANGLES OF THE POINT IN RADIANS0068C              BR,BTET,BPHI -  SPHERICAL COMPONENTS OF THE FIELD0063C-----OUTPUT:  BX,BY,BZ - CARTESIAN COMPONENTS OF THE FIELD0005C0005C0053C                   AUTHOR: NIKOLAI A. TSYGANENKO0052C                           INSTITUTE OF PHYSICS0063C                           ST.-PETERSBURG STATE UNIVERSITY0053C                           STARY PETERGOF 1989040046C                           ST.-PETERSBURG0038C                           RUSSIA0005C0018	IMPLICIT NONE00040053	REAL TETA,PHI,BR,BTET,BPHI,BX,BY,BZ,S,C,SF,CF,BE0004^^^^^^^0021      S=SIN(TETA)0021      C=COS(TETA)0021      SF=SIN(PHI)0021      CF=COS(PHI)0024      BE=BR*S+BTET*C0026      BX=BE*CF-BPHI*SF0026      BY=BE*SF+BPHI*CF0024      BZ=BR*C-BTET*S0016      RETURN0013      END0005c0075C----------------------------------------------------------------------0005C0052      SUBROUTINE RECALC(IYR,IDAY,IHOUR,MIN,ISEC)0005C0005C0068C	THIS IS A MODIFIED VERSION OF THE SUBROUTINE RECOMP WRITTEN BY0069C	N.A. TSYGANENKO.  SINCE I WANT TO USE IT IN PLACE OF SUBROUTINE0066C	RECALC I HAVE RENAMED THIS ROUTINE RECALC AND ELIMINATED THE0069C	ORIGINAL RECALC FROM THIS VERSION OF THE GEOPACK.F PACKET. THIS0069C	WAY ALL ORIGINAL CALLS TO RECALC WILL CONTINUE TO WORK WITHOUT 0047C	HAVING TO CHANGE THEM TO CALLS TO RECOMP.0005C0005C0079C  AN ALTERNATIVE VERSION OF THE SUBROUTINE RECALC FROM THE GEOPACK PACKAGE0078C  BASED ON A DIFFERENT APPROACH TO DERIVATION OF ROTATION MATRIX ELEMENTS0005C0081C THIS SUBROUTINE WORKS BY 20% FASTER THAN RECALC AND IS EASIER TO UNDERSTAND0005C0056C   ################################################0056C   #  WRITTEN BY  N.A. TSYGANENKO ON DEC.1, 1991  #0056C   ################################################0005C0005c0034c	Modified by:	Mauricio Peredo0040c			Hughes STX at NASA/GSFC Code 6950022c			September 19920005c0005c0063c	Modified to accept dates up to year 2000 and updated IGRF0028C	coeeficients for 1985.0005c0049C   OTHER SUBROUTINES CALLED BY THIS ONE: SUN0005C0041C     IYR = YEAR NUMBER (FOUR DIGITS)0044C     IDAY = DAY OF YEAR (DAY 1 = JAN 1)0040C     IHOUR = HOUR OF DAY (00 TO 23)0041C     MIN = MINUTE OF HOUR (00 TO 59)0041C     ISEC = SECONDS OF DAY(00 TO 59)0005C0018	IMPLICIT NONE00040062	REAL ST0,CT0,SL0,CL0,CTCL,STCL,CTSL,STSL,SFI,CFI,SPS,CPS,0069     1       SHI,CHI,HI,PSI,XMUT,A11,A21,A31,A12,A22,A32,A13,A23,0063     2       A33,DS3,F2,F1,G10,G11,H11,DT,SQ,SQQ,SQR,S1,S2,0068     3       S3,CGST,SGST,DIP1,DIP2,DIP3,Y1,Y2,Y3,Y,Z1,Z2,Z3,DJ,0070     4       T,OBLIQ,DZ1,DZ2,DZ3,DY1,DY2,DY3,EXMAGX,EXMAGY,EXMAGZ,^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^0057     5       EYMAGX,EYMAGY,GST,SLONG,SRASN,SDEC,BA(8)00040053	INTEGER IYR,IDAY,IHOUR,MIN,ISEC,K,IY,IDE,IYE,IPR00040074       COMMON/C1/ ST0,CT0,SL0,CL0,CTCL,STCL,CTSL,STSL,SFI,CFI,SPS,CPS,0071     * SHI,CHI,HI,PSI,XMUT,A11,A21,A31,A12,A22,A32,A13,A23,A33,DS3,0018     * K,IY,BA0005c0031      DATA IYE,IDE,IPR/3*0/0048      IF (IYR.EQ.IYE.AND.IDAY.EQ.IDE) GOTO 50005C0064C  IYE AND IDE ARE THE CURRENT VALUES OF YEAR AND DAY NUMBER0005C0016      IY=IYR0018      IDE=IDAY0032      IF(IY.LT.1965) IY=19650032      IF(IY.GT.2000) IY=20000005C0051C  WE ARE RESTRICTED BY THE INTERVAL 1965-2000,0081C  FOR WHICH THE IGRF COEFFICIENTS ARE KNOWN; IF IYR IS OUTSIDE THIS INTERVAL0081C   THE SUBROUTINE GIVES A WARNING (BUT DOES NOT REPEAT IT AT THE NEXT CALLS)0005C0052      IF(IY.NE.IYR.AND.IPR.EQ.0) PRINT 10,IYR,IY0029      IF(IY.NE.IYR) IPR=10016      IYE=IY0005C0074C  LINEAR INTERPOLATION OF THE GEODIPOLE MOMENT COMPONENTS BETWEEN THE0037C  VALUES FOR THE NEAREST EPOCHS:0005C0039	IF (IY.LT.1970) THEN				!1965-19700048	   F2=(FLOAT(IY)+FLOAT(IDAY)/365.-1965.)/5.0018	   F1=1.D0-F20031	   G10=30334.*F1+30220.*F20030	   G11=-2119.*F1-2068.*F20029	   H11=5776.*F1+5737.*F20042	ELSEIF (IY.LT.1975) THEN			!1970-19750048	   F2=(FLOAT(IY)+FLOAT(IDAY)/365.-1970.)/5.0018	   F1=1.D0-F20031	   G10=30220.*F1+30100.*F20030	   G11=-2068.*F1-2013.*F20029	   H11=5737.*F1+5675.*F20042	ELSEIF (IY.LT.1980) THEN			!1975-19800050	   F2=(DFLOAT(IY)+DFLOAT(IDAY)/365.-1975.)/5.0018	   F1=1.D0-F20031	   G10=30100.*F1+29992.*F20030	   G11=-2013.*F1-1956.*F20029	   H11=5675.*F1+5604.*F20043 	ELSEIF (IY.LT.1985) THEN			!1980-19850048	   F2=(FLOAT(IY)+FLOAT(IDAY)/365.-1980.)/5.0018	   F1=1.D0-F20031	   G10=29992.*F1+29873.*F20030	   G11=-1956.*F1-1905.*F20029	   H11=5604.*F1+5500.*F20042	ELSEIF (IY.LT.1990) THEN			!1985-19900048	   F2=(FLOAT(IY)+FLOAT(IDAY)/365.-1985.)/5.0018	   F1=1.D0-F20031	   G10=29873.*F1+29775.*F20030	   G11=-1905.*F1-1851.*F20029	   H11=5500.*F1+5411.*F20025	ELSE						!1990-20000043	   DT=FLOAT(IY)+FLOAT(IDAY)/365.-1990.^^^^^^^^^^^^0026	   G10=29775.-18.0*DT0026	   G11=-1851.+10.6*DT0025	   H11=5411.-16.1*DT0010	ENDIF0005C0081C  NOW CALCULATE THE COMPONENTS OF THE UNIT VECTOR EzMAG IN GEO COORD.SYSTEM:0072C   SIN(TETA0)*COS(LAMBDA0), SIN(TETA0)*SIN(LAMBDA0), AND COS(TETA0)0067C         ST0 * CL0                ST0 * SL0                CT00005C0026      SQ=G11**2+H11**20022      SQQ=SQRT(SQ)0029      SQR=SQRT(G10**2+SQ)0022      SL0=-H11/SQQ0022      CL0=-G11/SQQ0021      ST0=SQQ/SQR0021      CT0=G10/SQR0022      STCL=ST0*CL00022      STSL=ST0*SL00022      CTSL=CT0*SL00022      CTCL=CT0*CL00005C0073C      THE CALCULATIONS ARE TERMINATED IF ONLY GEO-MAG TRANSFORMATION0076C       IS TO BE DONE  (IHOUR>24 IS THE AGREED CONDITION FOR THIS CASE):0005C0034   5   IF (IHOUR.GT.24) RETURN0005C0063      CALL SUN(IY,IDAY,IHOUR,MIN,ISEC,GST,SLONG,SRASN,SDEC)0005C0077C  S1,S2, AND S3 ARE THE COMPONENTS OF THE UNIT VECTOR EXGSM=EXGSE IN THE0063C   SYSTEM GEI POINTING FROM THE EARTH'S CENTER TO THE SUN:0005C0033      S1=COS(SRASN)*COS(SDEC)0033      S2=SIN(SRASN)*COS(SDEC)0022      S3=SIN(SDEC)0023      CGST=COS(GST)0023      SGST=SIN(GST)0005C0076C  DIP1, DIP2, AND DIP3 ARE THE COMPONENTS OF THE UNIT VECTOR EZSM=EZMAG0026C   IN THE SYSTEM GEI:0005C0034      DIP1=STCL*CGST-STSL*SGST0034      DIP2=STCL*SGST+STSL*CGST0018      DIP3=CT00005C0078C  NOW CALCULATE THE COMPONENTS OF THE UNIT VECTOR EYGSM IN THE SYSTEM GEI0077C   BY TAKING THE VECTOR PRODUCT D x S AND NORMALIZING IT TO UNIT LENGTH:0005C0028      Y1=DIP2*S3-DIP3*S20028      Y2=DIP3*S1-DIP1*S30028      Y3=DIP1*S2-DIP2*S10035      Y=SQRT(Y1*Y1+Y2*Y2+Y3*Y3)0017      Y1=Y1/Y0017      Y2=Y2/Y0017      Y3=Y3/Y0005C0080C   THEN IN THE GEI SYSTEM THE UNIT VECTOR Z = EZGSM = EXGSM x EYGSM = S x Y0028C    HAS THE COMPONENTS:0005C0024      Z1=S2*Y3-S3*Y20024      Z2=S3*Y1-S1*Y30024      Z3=S1*Y2-S2*Y10005C0077C    THE VECTOR EZGSE (HERE DZ) IN GEI HAS THE COMPONENTS (0,-SIN(DELTA),0078C     COS(DELTA)) = (0.,-0.397823,0.917462); HERE DELTA = 23.44214 DEG FOR^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^0082C   THE EPOCH 1978 (SEE THE BOOK BY GUREVICH OR OTHER ASTRONOMICAL HANDBOOKS).0063C    HERE THE MOST ACCURATE TIME-DEPENDENT FORMULA IS USED:0005C0074      DJ=FLOAT(365*(IY-1900)+(IY-1901)/4 +IDAY)-0.5+FLOAT(ISEC)/86400.0021      T=DJ/36525.0049      OBLIQ=(23.45229-0.0130125*T)/57.29577950016      DZ1=0.0025      DZ2=-SIN(OBLIQ)0024      DZ3=COS(OBLIQ)0005C0078C  THEN THE UNIT VECTOR EYGSE IN GEI SYSTEM IS THE VECTOR PRODUCT DZ x S :0005C0027      DY1=DZ2*S3-DZ3*S20027      DY2=DZ3*S1-DZ1*S30027      DY3=DZ1*S2-DZ2*S10005C0070C   THE ELEMENTS OF THE MATRIX GSE TO GSM ARE THE SCALAR PRODUCTS:0081C   CHI=EM22=(EYGSM,EYGSE), SHI=EM23=(EYGSM,EZGSE), EM32=(EZGSM,EYGSE)=-EM23,0037C     AND EM33=(EZGSM,EZGSE)=EM220005C0034      CHI=Y1*DY1+Y2*DY2+Y3*DY30034      SHI=Y1*DZ1+Y2*DZ2+Y3*DZ30022      HI=ASIN(SHI)0005C0042C    TILT ANGLE: PSI=ARCSIN(DIP,EXGSM)0005C0037      SPS=DIP1*S1+DIP2*S2+DIP3*S30029      CPS=SQRT(1.-SPS**2)0023      PSI=ASIN(SPS)0005C0070C    THE ELEMENTS OF THE MATRIX MAG TO SM ARE THE SCALAR PRODUCTS:0083C CFI=GM22=(EYSM,EYMAG), SFI=GM23=(EYSM,EXMAG); THEY CAN BE DERIVED AS FOLLOWS:0005C0083C IN GEO THE VECTORS EXMAG AND EYMAG HAVE THE COMPONENTS (CT0*CL0,CT0*SL0,-ST0)0075C  AND (-SL0,CL0,0), RESPECTIVELY.    HENCE, IN GEI THE COMPONENTS ARE:0050C  EXMAG:    CT0*CL0*COS(GST)-CT0*SL0*SIN(GST)0050C            CT0*CL0*SIN(GST)+CT0*SL0*COS(GST)0021C            -ST00043C  EYMAG:    -SL0*COS(GST)-CL0*SIN(GST)0043C            -SL0*SIN(GST)+CL0*COS(GST)0019C             00072C  THE COMPONENTS OF EYSM IN GEI WERE FOUND ABOVE AS Y1, Y2, AND Y3;0071C  NOW WE ONLY HAVE TO COMBINE THE QUANTITIES INTO SCALAR PRODUCTS:0005C0040      EXMAGX=CT0*(CL0*CGST-SL0*SGST)0040      EXMAGY=CT0*(CL0*SGST+SL0*CGST)0021      EXMAGZ=-ST00037      EYMAGX=-(SL0*CGST+CL0*SGST)0037      EYMAGY=-(SL0*SGST-CL0*CGST)0033      CFI=Y1*EYMAGX+Y2*EYMAGY0043      SFI=Y1*EXMAGX+Y2*EXMAGY+Y3*EXMAGZ0005C0057      XMUT=(ATAN2(SFI,CFI)+3.1415926536)*3.81971863420005C0069C  THE ELEMENTS OF THE MATRIX GEO TO GSM ARE THE SCALAR PRODUCTS:0005C^^^^0064C   A11=(EXGEO,EXGSM), A12=(EYGEO,EXGSM), A13=(EZGEO,EXGSM),0064C   A21=(EXGEO,EYGSM), A22=(EYGEO,EYGSM), A23=(EZGEO,EYGSM),0064C   A31=(EXGEO,EZGSM), A32=(EYGEO,EZGSM), A33=(EZGEO,EZGSM),0005C0068C   ALL THE UNIT VECTORS IN BRACKETS ARE ALREADY DEFINED IN GEI:0005C0063C  EXGEO=(CGST,SGST,0), EYGEO=(-SGST,CGST,0), EZGEO=(0,0,1)0062C  EXGSM=(S1,S2,S3),  EYGSM=(Y1,Y2,Y3),   EZGSM=(Z1,Z2,Z3)0079C                                                           AND  THEREFORE:0005C0029      A11=S1*CGST+S2*SGST0030      A12=-S1*SGST+S2*CGST0016      A13=S30029      A21=Y1*CGST+Y2*SGST0030      A22=-Y1*SGST+Y2*CGST0016      A23=Y30029      A31=Z1*CGST+Z2*SGST0030      A32=-Z1*SGST+Z2*CGST0016      A33=Z30005C0022 10   FORMAT(//1X,0075     * '****RECOMP WARNS:  YEAR IS OUT OF INTERVAL 1965-2000: IYR=',I4,0058     * /,6X,'CALCULATIONS WILL BE DONE FOR IYR=',I4,/)0016      RETURN0013      END0079C--------------------------------------------------------------------------0005C0064      SUBROUTINE GEOMAG(XGEO,YGEO,ZGEO,XMAG,YMAG,ZMAG,J,IYR)0005C0074C CONVERTS GEOGRAPHIC (GEO) TO DIPOLE (MAG) COORDINATES OR VICA VERSA.0039C IYR IS YEAR NUMBER (FOUR DIGITS).0005C0054C                           J>0                J<00061C-----INPUT:  J,XGEO,YGEO,ZGEO,IYR   J,XMAG,YMAG,ZMAG,IYR0057C-----OUTPUT:    XMAG,YMAG,ZMAG        XGEO,YGEO,ZGEO0005C0005C0053C                   AUTHOR: NIKOLAI A. TSYGANENKO0052C                           INSTITUTE OF PHYSICS0063C                           ST.-PETERSBURG STATE UNIVERSITY0053C                           STARY PETERGOF 1989040046C                           ST.-PETERSBURG0038C                           RUSSIA0005C0018	IMPLICIT NONE00040061	REAL XGEO,YGEO,ZGEO,XMAG,YMAG,ZMAG,ST0,CT0,SL0,CL0,CTCL,0044     1       STCL,CTSL,STSL,AB(19),BB(8)00040026	INTEGER J,IYR,K,IY,II00040067      COMMON/C1/ ST0,CT0,SL0,CL0,CTCL,STCL,CTSL,STSL,AB,K,IY,BB0020      DATA II/1/0030      IF(IYR.EQ.II) GOTO 10016      II=IYR0034      CALL RECALC(II,0,25,0,0)0018  1   CONTINUE0027      IF(J.LT.0) GOTO 2^^^^^^^^^^^^^^^^^^^0043      XMAG=XGEO*CTCL+YGEO*CTSL-ZGEO*ST00032      YMAG=YGEO*CL0-XGEO*SL00043      ZMAG=XGEO*STCL+YGEO*STSL+ZGEO*CT00016      RETURN0043  2   XGEO=XMAG*CTCL-YMAG*SL0+ZMAG*STCL0043      YGEO=XMAG*CTSL+YMAG*CL0+ZMAG*STSL0032      ZGEO=ZMAG*CT0-XMAG*ST00016      RETURN0013      END0077C------------------------------------------------------------------------0005C0056      SUBROUTINE MAGSM(XMAG,YMAG,ZMAG,XSM,YSM,ZSM,J)0005C0076C CONVERTS DIPOLE (MAG) TO SOLAR MAGNETIC (SM) COORDINATES OR VICA VERSA0005C0045C                    J>0              J<00051C-----INPUT: J,XMAG,YMAG,ZMAG     J,XSM,YSM,ZSM0052C----OUTPUT:    XSM,YSM,ZSM       XMAG,YMAG,ZMAG0080C  ATTENTION:  SUBROUTINE  RECALC  MUST BE CALLED BEFORE MAGSM IN TWO CASES:0044C     /A/  BEFORE THE FIRST USE OF MAGSM0079C     /B/  IF THE CURRENT VALUES OF IYEAR,IDAY,IHOUR,MIN,ISEC ARE DIFFERENT0057C          FROM THOSE IN THE PRECEDING CALL OF  MAGSM0005C0053C                   AUTHOR: NIKOLAI A. TSYGANENKO0052C                           INSTITUTE OF PHYSICS0063C                           ST.-PETERSBURG STATE UNIVERSITY0053C                           STARY PETERGOF 1989040046C                           ST.-PETERSBURG0038C                           RUSSIA0005C0018	IMPLICIT NONE00040055	REAL XMAG,YMAG,ZMAG,XSM,YSM,ZSM,SFI,CFI,A(8),B(7),0029     1       AB(10),BA(8)00040019	INTEGER J,K,IY00040043      COMMON/C1/ A,SFI,CFI,B,AB,K,IY,BA0028      IF (J.LT.0) GOTO 10031      XSM=XMAG*CFI-YMAG*SFI0031      YSM=XMAG*SFI+YMAG*CFI0018      ZSM=ZMAG0016      RETURN0030  1   XMAG=XSM*CFI+YSM*SFI0030      YMAG=YSM*CFI-XSM*SFI0018      ZMAG=ZSM0016      RETURN0013      END0078C-------------------------------------------------------------------------0005C0061       SUBROUTINE GSMGSE(XGSM,YGSM,ZGSM,XGSE,YGSE,ZGSE,J)0005C0074C CONVERTS SOLAR MAGNETOSPHERIC (GSM) TO SOLAR ECLIPTICAL (GSE) COORDS0022C   OR VICA VERSA.0047C                    J>0                J<00053C-----INPUT: J,XGSM,YGSM,ZGSM    J,XGSE,YGSE,ZGSE0053C----OUTPUT:   XGSE,YGSE,ZGSE      XGSM,YGSM,ZGSM^^^^^^^^^^^^^^^0081C  ATTENTION:  SUBROUTINE  RECALC  MUST BE CALLED BEFORE GSMGSE IN TWO CASES:0046C     /A/  BEFORE THE FIRST CALL OF GSMGSE0079C     /B/  IF THE CURRENT VALUES OF IYEAR,IDAY,IHOUR,MIN,ISEC ARE DIFFERENT0059C          FROM THOSE IN THE PRECEDING  CALL  OF GSMGSE0005C0005C0053C                   AUTHOR: NIKOLAI A. TSYGANENKO0052C                           INSTITUTE OF PHYSICS0063C                           ST.-PETERSBURG STATE UNIVERSITY0053C                           STARY PETERGOF 1989040046C                           ST.-PETERSBURG0038C                           RUSSIA0005C0018	IMPLICIT NONE00040048	REAL XGSM,YGSM,ZGSM,XGSE,YGSE,ZGSE,SHI,CHI,0035     1       A(12),AB(13),BA(8)0019	INTEGER J,K,IY00040041      COMMON/C1/ A,SHI,CHI,AB,K,IY,BA0027      IF(J.LT.0) GOTO 10019      XGSE=XGSM0032      YGSE=YGSM*CHI-ZGSM*SHI0032      ZGSE=YGSM*SHI+ZGSM*CHI0016      RETURN00191     XGSM=XGSE0032      YGSM=YGSE*CHI+ZGSE*SHI0032      ZGSM=ZGSE*CHI-YGSE*SHI0016      RETURN0013      END0079C--------------------------------------------------------------------------0005C0057       SUBROUTINE SMGSM(XSM,YSM,ZSM,XGSM,YGSM,ZGSM,J)0005C0076C CONVERTS SOLAR MAGNETIC (SM) TO SOLAR MAGNETOSPHERIC (GSM) COORDINATES0022C   OR VICA VERSA.0046C                  J>0                 J<00054C-----INPUT: J,XSM,YSM,ZSM        J,XGSM,YGSM,ZGSM0050C----OUTPUT:  XGSM,YGSM,ZGSM       XSM,YSM,ZSM0005C0078C  ATTENTION:  SUBROUTINE RECALC MUST BE CALLED BEFORE SMGSM IN TWO CASES:0044C     /A/  BEFORE THE FIRST USE OF SMGSM0079C     /B/  IF THE CURRENT VALUES OF IYEAR,IDAY,IHOUR,MIN,ISEC ARE DIFFERENT0056C          FROM THOSE IN THE PRECEDING CALL OF SMGSM0005C0005C0050C                AUTHOR: NIKOLAI A. TSYGANENKO0049C                        INSTITUTE OF PHYSICS0055C                        ST.-PETERSBURG  UNIVERSITY0050C                        STARY PETERGOF 1989040043C                        ST.-PETERSBURG0037C                        U.S.S.R.0005C0018	IMPLICIT NONE00040062	REAL XSM,YSM,ZSM,XGSM,YGSM,ZGSM,SPS,CPS,A(10),B(15),AB(8)^^^^^^^^^^^^^^^^^0019	INTEGER J,K,IY00040040      COMMON/C1/ A,SPS,CPS,B,K,IY,AB0028      IF (J.LT.0) GOTO 10030      XGSM=XSM*CPS+ZSM*SPS0018      YGSM=YSM0030      ZGSM=ZSM*CPS-XSM*SPS0016      RETURN0031  1   XSM=XGSM*CPS-ZGSM*SPS0018      YSM=YGSM0031      ZSM=XGSM*SPS+ZGSM*CPS0016      RETURN0013      END0080C---------------------------------------------------------------------------0005C0060      SUBROUTINE GEOGSM(XGEO,YGEO,ZGEO,XGSM,YGSM,ZGSM,J)0005C0076C CONVERTS GEOGRAPHIC TO SOLAR MAGNETOSPHERIC COORDINATES OR VICA VERSA.0005C0049C                   J>0                   J<00055C----- INPUT:  J,XGEO,YGEO,ZGEO    J,XGSM,YGSM,ZGSM0055C---- OUTPUT:    XGSM,YGSM,ZGSM      XGEO,YGEO,ZGEO0081C  ATTENTION:  SUBROUTINE  RECALC  MUST BE CALLED BEFORE GEOGSM IN TWO CASES:0045C     /A/  BEFORE THE FIRST USE OF GEOGSM0080C     /B/  IF THE CURRENT VALUES OF IYEAR,IDAY,IHOUR,MIN,ISEC  ARE DIFFERENT0067C            FROM THOSE IN THE PREVIOUS CALL OF THIS SUBROUTINE0005C0005C0053C                   AUTHOR: NIKOLAI A. TSYGANENKO0052C                           INSTITUTE OF PHYSICS0063C                           ST.-PETERSBURG STATE UNIVERSITY0053C                           STARY PETERGOF 1989040046C                           ST.-PETERSBURG0038C                           RUSSIA0005C0018	IMPLICIT NONE00040056	REAL XGEO,YGEO,ZGEO,XGSM,YGSM,ZGSM,A11,A21,A31,A12,0050     1       A22,A32,A13,A23,A33,D,AA(17),B(8)0019	INTEGER J,K,IY00040068      COMMON/C1/ AA,A11,A21,A31,A12,A22,A32,A13,A23,A33,D,K,IY,B0028      IF (J.LT.0) GOTO 10041      XGSM=A11*XGEO+A12*YGEO+A13*ZGEO0041      YGSM=A21*XGEO+A22*YGEO+A23*ZGEO0041      ZGSM=A31*XGEO+A32*YGEO+A33*ZGEO0016      RETURN0041  1   XGEO=A11*XGSM+A21*YGSM+A31*ZGSM0041      YGEO=A12*XGSM+A22*YGSM+A32*ZGSM0041      ZGEO=A13*XGSM+A23*YGSM+A33*ZGSM0016      RETURN0013      END0078C-------------------------------------------------------------------------0005C0054      SUBROUTINE RHAND(X,Y,Z,R1,R2,R3,IOPT,EXNAME)0005C0076C COMPUTES RIGHT HAND EXPRESSIONS IN THE GEOMAGNETIC FIELD LINE EQUATION^^^^^^^^^^^^^^^0055C      (a subsidiary subroutine for the subr. STEP)0005C0053C                   AUTHOR: NIKOLAI A. TSYGANENKO0052C                           INSTITUTE OF PHYSICS0063C                           ST.-PETERSBURG STATE UNIVERSITY0053C                           STARY PETERGOF 1989040046C                           ST.-PETERSBURG0038C                           RUSSIA0005C0018	IMPLICIT NONE00040061	REAL X,Y,Z,R1,R2,R3,PSI,DS3,BX,BY,BZ,A(15),AA(10),BB(8),0060     1       XG,YG,ZG,R,T,F,BR,BT,BF,FX,FY,FZ,HX,HY,HZ,B0022	INTEGER IOPT,K,IY0020	EXTERNAL EXNAME00040041      COMMON/C1/ A,PSI,AA,DS3,K,IY,BB0046      CALL EXNAME(IOPT,PSI,X,Y,Z,BX,BY,BZ)0028      IF (K.EQ.0) GOTO 10040      CALL GEOGSM(XG,YG,ZG,X,Y,Z,-1)0040      CALL SPHCAR(R,T,F,XG,YG,ZG,-1)0040      CALL IGRF(IY,K,R,T,F,BR,BT,BF)0044      CALL BSPCAR(T,F,BR,BT,BF,FX,FY,FZ)0042      CALL GEOGSM(FX,FY,FZ,HX,HY,HZ,1)0016      GOTO 20038  1   CALL DIP(PSI,X,Y,Z,HX,HY,HZ)0018  2   BX=BX+HX0018      BY=BY+HY0018      BZ=BZ+HZ0039      B=DS3/SQRT(BX**2+BY**2+BZ**2)0017      R1=BX*B0017      R2=BY*B0017      R3=BZ*B0016      RETURN0013      END0077C------------------------------------------------------------------------0005C0055      SUBROUTINE STEP(N,X,Y,Z,DS,ERRIN,IOPT,EXNAME)0005C0076C RE-CALCULATES COORDS X,Y,Z FOR ONE STEP ALONG FIELD LINE. N IS MAXIMUM0075C ORDER OF HARMONICS IN MAIN FIELD EXPANSION, DS IS STEP SIZE, ERRIN IS0081C PERMISSIBLE ERROR VALUE, IOPT AND EXNAME - SEE COMMENTS TO SUBROUTINE TRACE0077C  ALL THE PARAMETERS ARE INPUT ONES; OUTPUT IS THE RENEWED TRIPLET X,Y,Z0005C0053C                   AUTHOR: NIKOLAI A. TSYGANENKO0052C                           INSTITUTE OF PHYSICS0063C                           ST.-PETERSBURG STATE UNIVERSITY0053C                           STARY PETERGOF 1989040046C                           ST.-PETERSBURG0038C                           RUSSIA0005C0018	IMPLICIT NONE00040053	REAL X,Y,Z,DS,ERRIN,DS3,R11,R12,R13,R21,R22,R23,0070     1       R31,R32,R33,R41,R42,R43,R51,R52,R53,ERRCUR,A(26),B(8)0024	INTEGER N,IOPT,K,IY0004^^0033      COMMON/C1/ A,DS3,K,IY,B0025      EXTERNAL EXNAME0013      K=N0020  1   DS3=-DS/3.0051      CALL RHAND(X,Y,Z,R11,R12,R13,IOPT,EXNAME)0063      CALL RHAND(X+R11,Y+R12,Z+R13,R21,R22,R23,IOPT,EXNAME)0056      CALL RHAND(X+.5*(R11+R21),Y+.5*(R12+R22),Z+.5*0044     *(R13+R23),R31,R32,R33,IOPT,EXNAME)0059      CALL RHAND(X+.375*(R11+3.*R31),Y+.375*(R12+3.*R320056     *),Z+.375*(R13+3.*R33),R41,R42,R43,IOPT,EXNAME)0058      CALL RHAND(X+1.5*(R11-3.*R31+4.*R41),Y+1.5*(R12-0051     *3.*R32+4.*R42),Z+1.5*(R13-3.*R33+4.*R43),0034     *R51,R52,R53,IOPT,EXNAME)0074      ERRCUR=ABS(R11-4.5*R31+4.*R41-.5*R51)+ABS(R12-4.5*R32+4.*R42-.5*0045     *R52)+ABS(R13-4.5*R33+4.*R43-.5*R53)0037      IF (ERRCUR.LT.ERRIN) GOTO 20018      DS=DS*.50016      GOTO 10033  2   X=X+.5*(R11+4.*R41+R51)0033      Y=Y+.5*(R12+4.*R42+R52)0033      Z=Z+.5*(R13+4.*R43+R53)0063      IF(ERRCUR.LT.ERRIN*.04.AND.ABS(DS).LT.1.33) DS=DS*1.50016      RETURN0013      END0075C----------------------------------------------------------------------0005C0069      SUBROUTINE TRACE(XI,YI,ZI,DIR,RLIM,R0,IHARM,NP,IOPT,EXNAME,0030     *XF,YF,ZF,XX,YY,ZZ,L)0005C0071C   TRACES FIELD LINE FROM ARBITRARY POINT OF SPACE UP TO THE EARTH0049C   SURFACE OR UP TO MODEL LIMITING BOUNDARY.0037C-------------- INPUT PARAMETERS:0064C   XI,YI,ZI - GSM COORDS OF INITIAL POINT (IN EARTH RADII),0071C   DIR - SIGN OF TRACING DIRECTION: IF DIR=1. THEN ANTIPARALLEL TO0068C     B VECTOR (E.G. FROM NORTHERN TO SOUTHERN CONJUGATE POINT),0044C     AND IF DIR=-1. THEN PARALLEL TO B.0076C   R0 IS RADIUS OF SPHERE (IN RE) FOR WHICH 'LANDING POINT' COORDINATES0040C     XF,YF,ZF  SHOULD BE CALCULATED0069C   RLIM - UPPER GEOCENTRIC DISTANCE WHICH LIMITS TRACING REGION.0075C   IHARM - MAXIMAL ORDER OF SPH.HARMONICS IN THE MAIN FIELD EXPANSION,0073C     IT DEPENDS ON THE CURRENT VERSION OF INTERNAL FIELD MODEL. SUB-0073C     ROUTINE IGRF OF THIS ISSUE CORRESPONDS TO IHARM=10. IF THE MAIN0072C     FIELD SHOULD BE ASSUMED TO BE PURELY DIPOLAR, THEN PUT IHARM=0^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^0068C     AND SPECIFY VALUE OF GEODIPOLE TILT ANGLE PSI (IN RADIANS)0072C     IN THE 16-TH ELEMENT OF THE COMMON BLOCK BEFORE CALLING TRACE.0064C     OTHERWISE CALL SUBROUTINE RECALC BEFORE CALLING TRACE.0067C   NP - UPPER ESTIMATE OF NUMBER OF STEPS ALONG THE FIELD LINE0045C     (OF THE ORDER OF SEVERAL HUNDREDS).0070C   IOPT - SPECIFIES OPTION OF EXTERNAL FIELD MODEL (E.G. LEVEL OF0071C     DISTURBANCE) - SEE COMMENTS TO SUBROUTINES EXTERN AND EXLONG.0068C   EXNAME - NAME OF THE EXTRATERRESTRIAL FIELD MODEL SUBROUTINE0047C     (EXTERN OR EXLONG OR SOME OTHER ONE).0038C-------------- OUTPUT PARAMETERS:0044C   XF,YF,ZF - GSM COORDS OF FINAL POINT0076C   XX,YY,ZZ - ARRAYS (LENGTH NP) CONTAINING COORDS OF FIELD LINE POINTS0072C   L - ACTUAL NUMBER OF FIELD LINE POINTS. IF L EXCEEDS NP, TRACING0048C     TERMINATES, AND A WARNING IS DISPLAYED0005C0005C0053C                   AUTHOR: NIKOLAI A. TSYGANENKO0052C                           INSTITUTE OF PHYSICS0063C                           ST.-PETERSBURG STATE UNIVERSITY0053C                           STARY PETERGOF 1989040046C                           ST.-PETERSBURG0038C                           RUSSIA0005C0018	IMPLICIT NONE00040063	REAL XI,YI,ZI,DIR,RLIM,R0,XF,YF,ZF,DD,ERR,DS,AA(26),BB(8),0058     1       X,Y,Z,AL,R1,R2,R3,AD,RR,RYZ,R,FC,XR,YR,ZR0038	INTEGER IHARM,NP,IOPT,L,K1,K2,J,I00040035      REAL XX(NP),YY(NP),ZZ(NP)0035      COMMON/C1/ AA,DD,K1,K2,BB0025      EXTERNAL EXNAME0067 10   FORMAT(//,1X,'**** COMPUTATIONS IN THE SUBROUTINE TRACE',0052     *' ARE TERMINATED: NP IS TOO SMALL ****'//)0017      J=IHARM0020      ERR=0.00050013      L=00020      DS=0.5*DIR0014      X=XI0014      Y=YI0014      Z=ZI0016      DD=DIR0018      K1=IHARM0015      AL=0.0048      CALL RHAND(X,Y,Z,R1,R2,R3,IOPT,EXNAME)0038      AD=SIGN(0.01,X*R1+Y*R2+Z*R3)0036      RR=SQRT(X**2+Y**2+Z**2)+AD0015  1   L=L+10028      IF(L.GT.NP) GOTO 70017      XX(L)=X0017      YY(L)=Y0017      ZZ(L)=Z0023      RYZ=Y**2+Z**20021      R2=X**2+RYZ0020      R=SQRT(R2)^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^0057      IF(R.GT.RLIM.OR.RYZ.GT.900..OR.X.GT.15.) GOTO 80040      IF(R.LT.R0.AND.RR.GT.R) GOTO 60028      IF(R.GE.RR) GOTO 50028      IF(R.GT.5.) GOTO 50028      IF(R.GE.3.) GOTO 30016      FC=0.20034      IF(R-R0.LT.0.05) FC=0.050026      AL=FC*(R-R0+0.2)0019      DS=DIR*AL0016      GOTO 40016  3   DS=DIR0015      AL=1.0014  4   XR=X0014      YR=Y0014      ZR=Z0014  5   RR=R0039      IF(IHARM.NE.0) J=1+10./(R-AL)0032      IF(J.GT.IHARM) J=IHARM0047      CALL STEP(J,X,Y,Z,DS,ERR,IOPT,EXNAME)0016      GOTO 10026  6   R1=(R0-R)/(RR-R)0023      X=X-(X-XR)*R10023      Y=Y-(Y-YR)*R10023      Z=Z-(Z-ZR)*R10016      GOTO 80021  7   write(*,10)0014      L=NP0016      RETURN0014  8   XF=X0014      YF=Y0014      ZF=Z0021      DO 9 I=L,NP0018      XX(I)=XF0018      YY(I)=YF0018  9   ZZ(I)=ZF0016      RETURN0013      END0077C------------------------------------------------------------------------0005c0005c0005C0005c0054        SUBROUTINE EXT87W(IOPT,PSI,X,Y,Z,BX,BY,BZ)0005c0057c	This is a dummy subroutine to provide the interface0056c	for the routines that do the T87W computation into0057c	the GEOPACK format.  Different version of T87W have0053c	been derived by M. Peredo based on the combined0057c	IMP/HEOS/ISEE database (as of Feb 1993).  Version 10058c	includes 12 different models where the data has been0058c	binned simultaneously by solar wind dynamic pressure0057c	and IMF-Bz values.  Version 2 includes 6 different 0055c	models based on ranges of the Kp index.  To avoid0056c	having to include many coefficients files, all the0054c	coefficients have been put into DATA statements,0056c	and the user selects the appropriate model via the0059c	IOPT parameter (as in the original GEOPACK routines).0054c	This subroutine calls subroutine FUNCTS which is0056c	actually used in the derivation of the models, and0054c	thus not only computes the magnetic field B, but0063c	also the derivates with respect to the model coefficients0063c	which are not of interest in this GEOPACK, but have been ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^0063c	left in there to avoid having to write a separate version0059c	of the model computation.  If you are concerned with 0061c	cpu resources or space you will probably want to remove0033c	those sections of the code.0005c0005c0005c0068c	IOPT VALUES FOR VERSION WITH DATA BINNED BY Solar Wind and IMF0005c0062C----	INPUT PARAMETERS: IOPT - for Solar wind pressure and0061c	IMF-Bz version SPECIFIES solar wind and IMF conditions:0005C0076C   IOPT=          101       102       103     104        105        1060005c0023C    CORRESPOND TO:0005c0077c  Psw (nPa) =   (0,1.5)   (1.5,2.5)   >2.5   (0,1.5)   (1.5,2.5)    >2.50078c  IMF-Bz(nT) =    >4         >4        >4     (0,4)      (0,4)      (0,4)0005c0005c0076C   IOPT=          107       108       109     110        111        1120005c0023C    CORRESPOND TO:0005c0077c  Psw (nPa) =   (0,1.5)   (1.5,2.5)   >2.5   (0,1.5)   (1.5,2.5)    >2.50076c  IMF-Bz(nT) =   (-4,0)     (-4,0)   (-4,0)    <-4       <-4        <-40005c0005C0045C    PS - GEODIPOLE TILT ANGLE IN RADIANS0063C    X, Y, Z  - GSM COORDINATES OF THE POINT IN EARTH RADII0005C0005c0052c	IOPT VALUES FOR VERSION WITH DATA BINNED BY Kp0005c0050C----	INPUT PARAMETERS: IOPT - for Kp version 0045c	SPECIFIES THE GROUND DISTURBANCE LEVEL:0005C0070C   IOPT= 1        2          3          4           5           60005c0023C    CORRESPOND TO:0005c0072C    KP= 0,0+    1-,1,1+    2-,2,2+    3-,3,3+     4-,4,4+      >=5-0005C0045C    PS - GEODIPOLE TILT ANGLE IN RADIANS0063C    X, Y, Z  - GSM COORDINATES OF THE POINT IN EARTH RADII0005C0075C----OUTPUT PARAMETERS: BX,BY,BZ - GSM COMPONENTS OF THE MODEL MAGNETIC0048C                        FIELD IN NANOTESLAS0005C0005C0042C              AUTHOR:	MAURICIO PEREDO0043C			HUGHES STX CORPORATION AT NASA/GSFC0057C			ISTP/GGS SCIENCE PLANNING AND OPERATIONS FACILITY0040C			NASA/GODDARD SPACE FLIGHT CENTER0016C			CODE 6950027C			GREENBELT, MD 207710005C0005C0035        IMPLICIT REAL (A-H,O-Z)0005c0029        COMMON/WARP/IWARP0031        COMMON/MODPAR/A0,MA^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^0064        DIMENSION A0(50),BMODL(3),DBDA(50,3),XI(4),LISTA(50)0005c0062C	Define arrays to store model coefficients CC# for models0060c	binned by solar wind conditions with # = 1-12, DD# for0059c	models binned by Kp # = 1-6. Then define coefficients0064c	in data statements and use value of IOPT to assigne to A0.0005c0063	DIMENSION CC1(30),CC2(30),CC3(30),CC4(30),CC5(30),CC6(30),0072     1            CC7(30),CC8(30),CC9(30),CC10(30),CC11(30),CC12(30)0005C0062	DIMENSION DD1(30),DD2(30),DD3(30),DD4(30),DD5(30),DD6(30)0005c0060	DATA CC1/-10.995829,-933.756676,-5843.461608,-4.572607,0069     1    4.000000,5.000000,9.355965,30.000000,1.615685,9.653729,0064     2    -34.731545,4.557584,31.564412,0.716487,-38.516600,0071     3    0.656160,48.290421,0.019605,0.127894,-37.080917,0.155075,0072     4    -0.002958,-1.415130,47.345877,-0.145129,0.027851,1.363943,0043     5    -0.006602,-0.006339,8.783022/0005c0061	DATA CC2/-13.248172,-1156.877177,-5828.473359,-4.282023,0070     1    4.000000,5.000000,6.787052,30.000000,1.459454,-8.781281,0064     2    -35.240314,3.898964,28.291729,0.934324,-31.551167,0071     3    0.271897,46.554661,0.058721,0.187937,-44.192214,0.161475,0072     4    -0.023661,-2.294721,52.864011,-0.115014,0.077023,1.592545,0043     5    -0.001752,-0.001883,1.671365/0005c0061	DATA CC3/-22.316235,-1808.087785,-5778.839891,-3.909181,0071     1    4.000000,5.000000,5.095250,30.000000,3.403580,-12.400002,0065     2    -87.131244,2.190164,27.987300,-1.961441,-54.497060,0072     3    3.377586,75.278621,-0.075011,0.338224,-75.714712,0.311287,0072     4    0.078189,0.776392,94.799245,-0.268782,-0.064356,-1.785332,0044     5    -0.018772,-0.006775,14.287582/0005c0059	DATA CC4/-7.931272,-963.417318,-5844.081178,-4.415401,0070     1    4.000000,5.000000,8.336015,30.000000,2.401653,12.812140,0065     2    -34.674483,3.790038,32.424094,-0.058155,-30.771768,0071     3    0.925526,38.290867,0.052702,0.062453,-33.988194,0.091117,0072     4    0.031624,-0.855770,40.778195,-0.070502,-0.010291,0.837109,^^^^^^^^^^^^^^^^^^^^^^^^^0043     5    -0.003457,-0.002969,6.091908/0005c0061	DATA CC5/-11.063954,-1090.623119,-5834.133724,-4.852020,0070     1    4.000000,5.000000,8.413674,30.000000,1.883874,12.840780,0065     2    -35.632183,4.231252,30.397919,-0.248066,-41.489682,0071     3    1.669672,56.244204,0.047346,0.075572,-42.956890,0.119860,0072     4    0.049894,-1.552812,54.280951,-0.095582,-0.041235,1.092763,0043     5    -0.001148,-0.002929,9.761491/0005c0060	DATA CC6/-6.784540,-1167.241524,-3039.881883,-4.533466,0070     1    4.000000,5.000000,8.194062,30.000000,2.589804,12.709873,0065     2    -58.663401,3.276254,30.756523,-1.496764,-27.420540,0072     3    3.148070,57.932875,0.005259,-0.069929,-58.359064,0.200928,0073     4    0.046648,-0.326635,75.632561,-0.197418,-0.043888,-0.272762,0043     5    -0.006044,-0.000245,9.158945/0005c0059	DATA CC7/-6.809388,-986.547070,-5842.105739,-4.283139,0070     1    4.000000,5.000000,9.593448,30.000000,2.448913,15.147771,0065     2    -49.444601,3.763814,31.539745,-0.770789,-21.717920,0071     3    1.771927,37.431627,0.010152,0.008205,-35.308285,0.088541,0072     4    0.054370,-0.792636,42.588444,-0.070980,-0.043926,0.641295,0043     5    -0.002519,-0.001816,8.627818/0005c0060	DATA CC8/-8.325135,-1062.591858,-5836.150808,-4.178678,0070     1    4.000000,5.000000,8.090198,30.000000,1.686869,15.158005,0065     2    -60.253322,3.348521,27.098434,-0.492334,-25.318088,0071     3    2.137669,43.057453,0.016511,0.071893,-38.856806,0.095539,0072     4    0.071196,-0.401377,47.694962,-0.069481,-0.074697,0.250798,0043     5    -0.005439,-0.004443,8.569006/0005c0060	DATA CC9/-1.000573,-1046.591650,-3655.387277,-4.157887,0070     1    4.000000,5.000000,7.479074,30.000000,1.913858,14.589570,0065     2    -60.882753,3.236338,25.339141,-1.908383,-35.240215,0071     3    4.015674,62.079877,0.092327,0.002652,-46.990432,0.171173,0073     4    0.092966,-0.166584,64.926193,-0.176772,-0.131028,-0.280559,0043     5    -0.010686,-0.004174,8.179988/0005c0061	DATA CC10/-9.758139,-1133.614663,-5831.984090,-3.620778,^^^^^^0070     1    4.000000,5.000000,6.657985,30.000000,2.231157,16.048939,0065     2    -60.964410,3.513972,30.251830,-1.065994,-34.069729,0073     3    2.889400,53.636569,0.115655,-0.258396,-24.231761,-0.016044,0071     4    0.029005,-3.939334,23.904930,0.081847,-0.017146,3.213165,0041     5    0.007211,0.000902,1.830427/0005c0060	DATA CC11/2.370642,-1154.831109,-5827.805797,-1.940344,0070     1    4.000000,5.000000,8.486536,30.000000,2.677452,15.070445,0065     2    -57.089650,3.954353,24.128469,-3.235949,-21.150270,0073     3    4.771500,57.867921,-0.040256,-0.172150,-35.728413,0.089310,0072     4    0.164077,-3.547037,46.940782,-0.073820,-0.211344,0.734539,0041     5    0.027804,0.006272,8.723675/0005c0060	DATA CC12/1.471088,-1358.404720,-2997.728720,-3.004939,0070     1    4.000000,5.000000,8.600677,30.000000,1.953163,14.096985,0064     2    -79.085812,1.861273,22.408443,-5.206454,16.582682,0072     3    7.634370,31.973577,0.004823,-0.276016,-44.643930,0.006612,0071     4    0.233930,-4.473169,51.238758,0.142571,-0.302632,4.809631,0043     5    -0.000917,-0.006846,9.885635/0005c0005C0059	DATA DD1/-5.867103,-883.738293,-5850.097616,-5.188391,0070     1    4.000000,5.000000,8.460622,30.000000,2.616945,12.799227,0065     2    -29.295097,3.470999,30.012121,-0.165900,-17.637193,0071     3    1.081232,32.113211,0.007220,0.005960,-28.524636,0.088043,0071     4    0.020503,0.020194,36.238796,-0.066922,-0.016724,0.008272,0043     5    -0.006085,-0.001389,6.678375/0005c0059	DATA DD2/-6.594555,-933.159395,-5846.402067,-4.642404,0070     1    4.000000,5.000000,8.370964,30.000000,2.328611,14.082344,0074     2    -34.807249,3.714065,28.461083,-0.270863,-22.326947,1.457701,0072     3    38.819806,0.015287,-0.005407,-32.756903,0.106588,0.031201,0073     4    -0.284590,42.074933,-0.088637,-0.035114,0.052602,-0.004850,0033     5    -0.001412,6.264312/0005c0060	DATA DD3/-5.666452,-1027.528801,-5839.524379,-4.544371,0070     1    4.000000,5.000000,8.457785,30.000000,2.425008,15.078018,^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^0065     2    -47.079808,3.570725,28.926932,-1.060386,-28.460225,0073     3     2.418158,46.661731,0.020598,-0.008466,-36.450002,0.112016,0073     4     0.050126,-0.353297,47.387990,-0.094750,-0.062448,0.130280,0044     5     -0.003003,-0.002618,9.488941/0005c0060	DATA DD4/-5.099814,-1058.816704,-5836.919924,-4.075619,0070     1    4.000000,5.000000,7.691272,30.000000,2.133305,15.703607,0065     2    -47.666381,4.076097,26.609532,-1.278862,-32.633222,0072     3    3.050778,52.891613,0.019174,-0.014803,-37.618785,0.115264,0072     4    0.058449,0.349566,50.725635,-0.106374,-0.082590,-0.676248,0043     5    -0.006160,-0.004273,8.215376/0005c0060	DATA DD5/-6.778978,-1128.003458,-5832.090531,-3.575514,0070     1    4.000000,5.000000,7.776490,30.000000,1.587892,15.397206,0065     2    -68.493384,3.294595,23.353014,-0.770688,-35.364872,0072     3    3.046147,59.877481,0.074692,-0.119341,-37.080440,0.100374,0073     4    0.073908,-0.419508,46.000258,-0.068167,-0.103513,-0.730740,0043     5    -0.002787,-0.003340,8.249916/0005c0061	DATA DD6/-13.206068,-1709.633188,-9597.093218,-1.695995,0070     1    4.000000,5.000000,7.612087,30.000000,1.983091,15.804040,0066     2    -119.766571,2.834474,24.440544,-2.357557,-38.950308,0072     3    4.649310,75.589106,0.014460,-0.161612,-48.407884,0.084010,0072     4    0.118535,-1.248589,56.197803,-0.007591,-0.165901,0.760503,0043     5    -0.009544,-0.001897,9.579404/0005c0005C0067C	Assign coefficients corresponding to model identified by IOPT0005c0017	DO 10 I=1,300037	   IF (IOPT.EQ.101) A0(I)=CC1(I)0037	   IF (IOPT.EQ.102) A0(I)=CC2(I)0037	   IF (IOPT.EQ.103) A0(I)=CC3(I)0037	   IF (IOPT.EQ.104) A0(I)=CC4(I)0037	   IF (IOPT.EQ.105) A0(I)=CC5(I)0037	   IF (IOPT.EQ.106) A0(I)=CC6(I)0037	   IF (IOPT.EQ.107) A0(I)=CC7(I)0037	   IF (IOPT.EQ.108) A0(I)=CC8(I)0037	   IF (IOPT.EQ.109) A0(I)=CC9(I)0038	   IF (IOPT.EQ.110) A0(I)=CC10(I)0038	   IF (IOPT.EQ.111) A0(I)=CC11(I)0038	   IF (IOPT.EQ.112) A0(I)=CC12(I)0035	   IF (IOPT.EQ.1) A0(I)=DD1(I)0035	   IF (IOPT.EQ.2) A0(I)=DD2(I)^^^^^^^^^^^^^0035	   IF (IOPT.EQ.3) A0(I)=DD3(I)0035	   IF (IOPT.EQ.4) A0(I)=DD4(I)0035	   IF (IOPT.EQ.5) A0(I)=DD5(I)0035	   IF (IOPT.EQ.6) A0(I)=DD6(I)0048	   IF ( (IOPT.LE.0).OR.(IOPT.GT.112) ) THEN0050		WRITE(*,*) 'IMPROPER VALUE OF IOPT -- ABORT'0012		RETURN0053	   ELSEIF ( (IOPT.GT.6).AND.(IOPT.LT.101) ) THEN0050		WRITE(*,*) 'IMPROPER VALUE OF IOPT -- ABORT'0012		RETURN0013	   ENDIF001510	CONTINUE0005C0005C0040	IWARP =1	!Warp turned ON by default0043	MA = 30		!This model has 30 parameters0005c0043        TILT = PSI * (180./3.141592654)0052        CALL FUNCTS(X,Y,Z,TILT,A0,BMODL,DBDA,MA)0025        BX = BMODL(1)0025        BY = BMODL(2)0025        BZ = BMODL(3)0005C0018        RETURN0015        END0005c0005c0005C0060           SUBROUTINE FUNCTS(X,Y,Z,TILT,A,BMODL,DBDA,MA)0005C0067C       This subroutine computes the 'model' magnetic field and0068c       derivatives with respect to fitting coefficients for the0072c       nonlinear least square fitting program MRQMIN.  It is called0068c       from subroutine MRQCOF which itself is called by MRQMIN.0077c       The first four arguments are the X,Y,Z,and TILT values (input) at0077c       the Ith data point (ie. from POS(I,1-4). On input, A contains the0075c       current choice of parameters to compute the model with while MA0079c       is the total number of parameters. Outputs are  BMODL a 3-component0078c       vector containing the 3 cartesian components of the magnetic field0080c       at the Ith data point (output), and DBDA a 50 by 3 matrix containing0072c       the derivatives with respect to the fitting coefficients (up0067c       to 50) of the 3 components of the model magnetic field.0005C0005C0037          IMPLICIT REAL (A-H,O-Z)0005C0031        COMMON/PARAM/PR(40)0065        DIMENSION BMODL(3),DBDA(50,3),A(50),BTAIL(3),BRING(3)0073        DIMENSION DBTAIL(50,3),DBRING(50,3),DBMPAUSE(50,3),BMPAUSE(3)0005C0005C0068c       Assign parameters from input array A to array PR used in0066c       David's subroutines to evaluate the T87 magnetic field0005c^^^^^^^^^^^^^^^^^^^0039        PR(1) = A(1)            !B00039        PR(2) = A(2)            !B10039        PR(3) = A(3)            !B20039        PR(4) = A(4)            !XN0039        PR(5) = A(5)            !X10039        PR(6) = A(6)            !X20039        PR(7) = A(7)            !RH0039        PR(8) = A(8)            !RT0059        PR(9) = ABS( A(9) )     !D -- ABS to insure D>00049        PR(10) = A(10)          !DY (delta-Y)0040        PR(13) = A(11)          !Brc0063        PR(14) = ABS( A(12) )   !Rrc -- ABS to insure Rrc>00051        PR(15) = A(13)          !DX1 (delta-X1)0039        PR(16) = A(14)          !a10039        PR(17) = A(15)          !a20039        PR(18) = A(16)          !a30039        PR(19) = A(17)          !a40039        PR(20) = A(18)          !a50039        PR(21) = A(19)          !a60005c0073c       PR(22) - PR(27) correspond to b1 - b6 which are derived below0046c       input list continues with c1 - c100005c0039        PR(28) = A(20)          !c10039        PR(29) = A(21)          !c20039        PR(30) = A(22)          !c30039        PR(31) = A(23)          !c40039        PR(32) = A(24)          !c50039        PR(33) = A(25)          !c60039        PR(34) = A(26)          !c70039        PR(35) = A(27)          !c80039        PR(36) = A(28)          !c90040        PR(37) = A(29)          !c100068        PR(38) = A(30)          !GWARP KOLYA'S WARPING AMPLITUDE0005c0030c       Now derive b1 - b60005c0068        DX2 = 0.5* PR(15)               !delta-X2 = 0.5 delta-X10005c0063        PR(22) = -( 2.*PR(30) + (PR(16)/PR(15)) )       !b10063        PR(23) = -( PR(31) + (PR(17)/PR(15)) )          !b20063        PR(24) = -( 2.*PR(34) + (PR(18)/DX2) )          !b30063        PR(25) = -( PR(35) + (PR(19)/DX2) )             !b40063        PR(26) = -( PR(36) + (PR(20)/DX2) )/3.          !b50063        PR(27) = -( 3.*PR(37) + (PR(21)/DX2) )          !b60005c0068C       All appropriate parameters for the model are defined now0068c       we are ready to call clone routines based on MPAUSE.FOR,^^^^^^^^^^^^^^^^^^^^0067c       RING87.FOR and TAIL87.FOR to compute the field.  Recall0069c       that we also need derivatives with respect ot each of the0069c       coefficients in array A (i.e. don't need derivatives with0066c       respect to bi's).  Clone subroutines have been adapted0065c       from subroutines created by David for use with LSFIT.0005c0078        CALL NLTAIL(X,Y,Z,TILT,MA,BTAIL,DBTAIL)         !Tail contribution0081        CALL NLRING(X,Y,Z,TILT,MA,BRING,DBRING)    !Ring current contribution0083        CALL NLMPAUSE(X,Y,Z,TILT,MA,BMPAUSE,DBMPAUSE)   !Magnetopause/Birkeland0005c0005C0072        BMODL(1) = BTAIL(1) + BRING(1) + BMPAUSE(1)     !X-component0072        BMODL(2) = BTAIL(2) + BRING(2) + BMPAUSE(2)     !Y-component0072        BMODL(3) = BTAIL(3) + BRING(3) + BMPAUSE(3)     !Z-component0005c0024        DO 11 I=1,MA0081          DBDA(I,1) = DBTAIL(I,1) + DBRING(I,1) + DBMPAUSE(I,1)  !X-component0081          DBDA(I,2) = DBTAIL(I,2) + DBRING(I,2) + DBMPAUSE(I,2)  !Y-component0081          DBDA(I,3) = DBTAIL(I,3) + DBRING(I,3) + DBMPAUSE(I,3)  !Z-component002011      CONTINUE0005C0018        RETURN0015        END0005C0005C0005c0057        SUBROUTINE NLTAIL(X,Y,Z,TILT,MA,BTAIL,DBTAIL)0005C0068C       This subroutine was adapted from David's LSTAIL routine.0078c       Modifications include: (i) replacing position and tilt from common0077c       i/o into arguments in call, (ii) call also includes arguments for0079c       number of parameters and arrays for output field and derivatives of0074c       field components w.r.t. parameters, (iii) computation of field0077c       components was specifically added from terms used in derivatives.0073c       In addition, computation of sine of tilt angle is done inside0076c       this subroutine now.  Parameters are in array PR (instead of PAR0032c       in original LSTAIL).0005c0081c       Input are X,Y,Z,TILT specifying the 4-dim position at which the field0081c       is to be computed, and MA, the number of parameters in the prescribed^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^0050c       representation for the magnetic field.0080c       Output are BTAIL, a 3-component array containing the three cartesian0080c       components of the tail current contribution to the external magnetic0083c       field, and DBTAIL, an MA by 3 matrix containing the partial derivatives0076c       of each of the 3 field components wrt each of the MA parameters.0005c0005c0037          IMPLICIT REAL (A-H,O-Z)0005C0031        COMMON/PARAM/PR(40)0029        COMMON/WARP/IWARP0043        DIMENSION BTAIL(3),DBTAIL(50,3)0005C0060C       Set PR(39) according to value of warp flag IWARP0005C0039        IF (IWARP.EQ.0) PR(39) = 0.0039        IF (IWARP.EQ.1) PR(39) = 1.0005C0066c       Initialize all field and derivative components to zero0071        DO 10 J=1,3                     !Loop over X,Y,Z components0029           BTAIL(J) = 0.00072           DO 10 I=1,MA                 !Loop over fitted parameters0034              DBTAIL(I,J) = 0.002010      CONTINUE00040005c0028        PI = 3.1415926540052        PIHF = 1.5707963                ! 0.5 pi0034        S6 = SIN(TILT*PI/180.)0034        XI1 = PR(5) - X ! xi-10050        XI12= XI1*XI1           !  "   squared0034        XIN = PR(4) - X ! xi-N0050        XIN2= XIN*XIN           !  "   squared0005c0070c       Information on warping scale length and computation of y^40005c0024        Y4 = Y*Y*Y*Y0026        GWARP = PR(38)0020        LY = 10.0029        LY4 = LY*LY*LY*LY0005C0073      IF (PR(39).EQ.0.) THEN                    ! PR(39) gives signal0048        Z0 = PR(7)*S6           ! no warping0060      ELSE                      !Kolya's warping formula0048        Z0 = PR(7)*S6 - GWARP*S6*Y4/(Y4+LY4)0016      END IF0059        ZR  = Z - Z0            ! z-r, median z of tail0052        ZPL = Z - PR(8)                 ! z-plus0053        ZMI = Z + PR(8)                 ! z-minus0055        D2 = PR(9)*PR(9)                ! D-squared0005c0058        BTA2 = ZR*ZR + D2               ! beta-squared0063        BPL2 = ZPL*ZPL + D2             ! beta-plus squared^^^^^^^^^^^^^^0064        BMI2 = ZMI*ZMI + D2             ! beta-minus squared0028        BTA = SQRT(BTA2)0028        BPL = SQRT(BPL2)0028        BMI = SQRT(BMI2)0005c0061        GX1  = XI12 + BTA2                      ! gamma-10066        GPL1 = XI12 + BPL2                      ! gamma-1-plus0067        GMI1 = XI12 + BMI2                      ! gamma-1-minus0061        X1D  = (PR(4)-PR(5))                    ! XN - X10069        X1D2 = X1D*X1D                          !    "    squared0005c0064        PX1  = 0.5*(LOG(X1D2/(XIN2 + BTA2)))/GX1        ! P10069        PPL1 = 0.5*(LOG(X1D2/(XIN2 + BPL2)))/GPL1       ! P1-plus0070        PMI1 = 0.5*(LOG(X1D2/(XIN2 + BMI2)))/GMI1       ! P1-minus0005c0066        XI2  = PR(6) - X                                ! xi-20074        XI22 = XI2*XI2                                  !   "  squared0069        GX2  = XI22 + BTA2                              ! gamma-20074        GPL2 = XI22 + BPL2                              ! gamma-2-plus0075        GMI2 = XI22 + BMI2                              ! gamma-2-minus0034          G2IN  = 1./(GX2*GX2)0036          GP2IN = 1./(GPL2*GPL2)0036          GM2IN = 1./(GMI2*GMI2)0069        X2D  = PR(4) - PR(6)                            ! XN - X20077        X2D2 = X2D*X2D                                  !    "    squared0064          PX2  =  G2IN*LOG(X2D2/(XIN2 + BTA2))          ! P20069          PPL2 = GP2IN*LOG(X2D2/(XIN2 + BPL2))          ! P2-plus0070          PMI2 = GM2IN*LOG(X2D2/(XIN2 + BMI2))          ! P2-minus0005c0076        S0   = (PIHF + ATAN2(XIN,BTA))/BTA     ! S0 ; ATAN of neg. angle0077        S0PL = (PIHF + ATAN2(XIN,BPL))/BPL     ! S0-plus;   should be neg0061        S0MI = (PIHF + ATAN2(XIN,BMI))/BMI     ! S0-minus0005c0056        S1   = PX1  - (XI1*S0)/GX1              ! S10061        S1PL = PPL1 - (XI1*S0PL)/GPL1           ! S1-plus0062        S1MI = PMI1 - (XI1*S0MI)/GMI1           ! S1-minus0005c0073        G0 = 0.5*LOG((XIN2 + BTA2)/SQRT((XIN2 + BPL2)*(XIN2 + BMI2)))0005c0056        G1   = (BTA2*S0)/GX1 + XI1*PX1          ! G1^^^0061        G1PL = (BPL2*S0PL)/GPL1 + XI1*PPL1      ! G1-plus0062        G1MI = (BMI2*S0MI)/GMI1 + XI1*PMI1      ! G1-minus0005c0031        T2   = 1./(X2D*GX2)0032        T2PL = 1./(X2D*GPL2)0032        T2MI = 1./(X2D*GMI2)0070        S2   = ((XI22-BTA2)*S0*G2IN - XI2*PX2) - T2           ! S20075        S2PL = ((XI22-BPL2)*S0PL*GP2IN - XI2*PPL2) - T2PL     ! S2-plus0076        S2MI = ((XI22-BMI2)*S0MI*GM2IN - XI2*PMI2) - T2MI     ! S2-minus0005c0067        G2   = 0.5*(BTA2-XI22)*PX2 - (2.*G2IN*BTA2*S0 + T2)*XI20073        G2PL = 0.5*(BPL2-XI22)*PPL2 - (2.*GP2IN*BPL2*S0PL + T2PL)*XI20073        G2MI = 0.5*(BMI2-XI22)*PMI2 - (2.*GM2IN*BMI2*S0MI + T2MI)*XI20005C0005C0052c                       field vector calculation0005C0025        YR = Y/PR(10)0074        FY = 1./(PI*(1.+ (YR*YR)))                              ! f(y)0005c0005c0059        TERMX0 = FY*(ZR*S0 - 0.5*(ZPL*S0PL + ZMI*S0MI))0059        TERMX1 = FY*(ZR*S1 - 0.5*(ZPL*S1PL + ZMI*S1MI))0059        TERMX2 = FY*(ZR*S2 - 0.5*(ZPL*S2PL + ZMI*S2MI))0026        TERMZ0 = FY*G00048        TERMZ1 = FY*(G1 - 0.5*(G1PL + G1MI))0048        TERMZ2 = FY*(G2 - 0.5*(G2PL + G2MI))0005C0072c       Note that Y-component is zero, but write comment as reminder0080        BTAIL(1) = PR(1)*TERMX0 + PR(2)*TERMX1 + PR(3)*TERMX2   !X-component0080c       BTAIL(2) = 0.                                           !Y-component0080        BTAIL(3) = PR(1)*TERMZ0 + PR(2)*TERMZ1 + PR(3)*TERMZ2   !Z-component0005C0058C       Now compute partial derivatives wrt parameters0005c0066c       All Y-components are zero so will not write them again0005c0060        DBTAIL(1,1) = TERMX0                    !DBx/DB00060        DBTAIL(2,1) = TERMX1                    !DBx/DB10060        DBTAIL(3,1) = TERMX2                    !DBx/DB20005c0060        DBTAIL(1,3) = TERMZ0                    !DBz/DB00060        DBTAIL(2,3) = TERMZ1                    !DBz/DB10060        DBTAIL(3,3) = TERMZ2                    !DBz/DB20005c0079c       The remaining derivatives are a bit more complicated and we'll code^^^^^^0026c       them in steps.0005c0005C0017c       D/DXN0005c0062        TEMPA1 = (ZR/(BTA2+XIN2))-0.5*((ZPL/(BPL2+XIN2)) +0039     1           (ZMI/(BMI2+XIN2)))0031        TEMPA2 = TEMPA1/X1D0032        TEMPA3 = TEMPA1/X2D20005C0062        TEMPA4 = ((XIN2-BTA2)*(BPL2+BMI2) + 2.*BPL2*BMI2 -0073     1   2.*XIN2*BTA2)*XIN/( 2.*(BTA2+XIN2)*(BPL2+XIN2)*(BMI2+XIN2) )0061        TEMPA5 = (1./(BTA2+XIN2))-0.5*((1./(BPL2+XIN2)) +0038     1           (1./(BMI2+XIN2)))0039        TEMPA6 = XIN * TEMPA5 / X1D0040        TEMPA7 = XIN * TEMPA5 / X2D20005C0075        DBTAIL(4,1) = FY*( PR(1)*TEMPA1 + PR(2)*TEMPA2 + PR(3)*TEMPA3 )0075        DBTAIL(4,3) = FY*( PR(1)*TEMPA4 + PR(2)*TEMPA6 + PR(3)*TEMPA7 )0005C0005C0017C       D/DX10005C0069        TEMPB1 = ( 2.*XI1*PX1 - (1./X1D) - (BTA2-XI12)*S0/GX1 ) *0030     1            ZR / GX10074        TEMPB2 = -0.5*(2.*XI1*PPL1 - (1./X1D) - (BPL2-XI12)*S0PL/GPL1)0031     1            *ZPL/GPL10074        TEMPB3 = -0.5*(2.*XI1*PMI1 - (1./X1D) - (BMI2-XI12)*S0MI/GMI1)0031     1            *ZMI/GMI10005C0075        TEMPB4 = ( (BTA2-XI12)*PX1 - 2.*BTA2*XI1*S0/GX1 - XI1/X1D )/GX10074        TEMPB5 = -0.5*((BPL2-XI12)*PPL1-2.*BPL2*XI1*S0PL/GPL1-XI1/X1D)0025     1          /GPL10074        TEMPB6 = -0.5*((BMI2-XI12)*PMI1-2.*BMI2*XI1*S0MI/GMI1-XI1/X1D)0025     1          /GMI10005C0065        DBTAIL(5,1) = FY * PR(2) * (TEMPB1 + TEMPB2 + TEMPB3)0065        DBTAIL(5,3) = FY * PR(2) * (TEMPB4 + TEMPB5 + TEMPB6)0005C0005C0017C       D/DX20005C0065        TEMPC1 = 4.*XI2*G2IN/X2D+(4.*XI22/GX2-1.)*PX2-T2/X2D+0059     1           2.*XI2*S0*(1.-2.*(XI22-BTA2)/GX2)*G2IN0070        TEMPC2 = 4.*XI2*GP2IN/X2D+(4.*XI22/GPL2-1.)*PPL2-T2PL/X2D+0063     1           2.*XI2*S0PL*(1.-2.*(XI22-BPL2)/GPL2)*GP2IN0070        TEMPC3 = 4.*XI2*GM2IN/X2D+(4.*XI22/GMI2-1.)*PMI2-T2MI/X2D+0063     1           2.*XI2*S0MI*(1.-2.*(XI22-BMI2)/GMI2)*GM2IN0005C0064        TEMPC4 = -XI2 * ( 1. + 2.* (BTA2-XI22)/GX2 ) * PX2 +0069     1           2. * BTA2 * ( 4.* XI22 - GX2 ) * S0 * G2IN / GX2^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^0067     2          -2. * (BTA2-XI22) * G2IN / X2D - XI2 * T2 / X2D0066        TEMPC5 = -XI2 * ( 1. + 2.* (BPL2-XI22)/GPL2 ) * PPL2 +0072     1         2. * BPL2 * ( 4.* XI22 - GPL2 ) * S0PL * GP2IN / GPL20068     2        -2. * (BPL2-XI22) * GP2IN / X2D - XI2 * T2PL / X2D0066        TEMPC6 = -XI2 * ( 1. + 2.* (BMI2-XI22)/GMI2 ) * PMI2 +0072     1         2. * BMI2 * ( 4.* XI22 - GMI2 ) * S0MI * GM2IN / GMI20068     2        -2. * (BMI2-XI22) * GM2IN / X2D - XI2 * T2MI / X2D0005C0074        DBTAIL(6,1) = FY*PR(3)*(ZR*TEMPC1-0.5*(ZPL*TEMPC2+ZMI*TEMPC3))0069        DBTAIL(6,3) = FY*PR(3)*(TEMPC4 - 0.5 * (TEMPC5 + TEMPC6))0005C0005C0017C       D/DRH0005C0073        TEMPD1 = (ZR*ZR/BTA2 -1.)*S0 + ZR*ZR*XIN/( BTA2*(XIN2+BTA2) )0072        TEMPD2 = -ZR*ZR*XI1*S0/(GX1*BTA2) + (2.*ZR*ZR/GX1 - 1.)*S1 -0068     1           ZR*ZR*(XI1*XIN-BTA2) / ( GX1*BTA2*(XIN2+BTA2) )0067        TEMPD3 = ZR*ZR*S0/(GX2*BTA2) + (4.*ZR*ZR/GX2 - 1.)*S2 +0063     1           ZR*ZR*G2IN*( 2./X2D - 2.*XI2/(XIN2+BTA2) +0057     2           XIN*(XI22-BTA2)/(BTA2*(XIN2+BTA2)) )0005C0036        TEMPD4 = -1./(XIN2+BTA2)0065        TEMPD5 = ( 2.*G1 - S0 + (XIN+XI1)/(XIN2+BTA2) ) / GX10068        TEMPD6 = 4.*G2/GX2 + 2.*XI2*G2IN*S0 - PX2 + 2.*XI2*G2IN/0069     1           X2D + (BTA2 -XI22 - 2.*XI2*XIN)*G2IN/(XIN2+BTA2)0005C0074        DBTAIL(7,1) = FY*S6*( PR(1)*TEMPD1+PR(2)*TEMPD2+PR(3)*TEMPD3 )0064        DBTAIL(7,3) = FY*S6*ZR*( PR(1)*TEMPD4 + PR(2)*TEMPD50052     1                          + PR(3)*TEMPD6 )0005C0005C0020C       D/DGWARP0005C0066C       DUE TO SIMILARITY TO D/DRH, INSERT HERE TERMS D/DGWARP0074C       WHICH FOLLOW BY SUBSTITUTING IN D/DRH TERMS S6 (I.E. SIN(TILT)0062C       BY - S6*Y4/(Y4+LY4), I.E. SIN(TILT)*Y^4/(Y^4+LY^4)0005C0051        DBTAIL(30,1) = -DBTAIL(7,1)*Y4/(Y4+LY4)0051        DBTAIL(30,3) = -DBTAIL(7,3)*Y4/(Y4+LY4)0005C0005C0017C       D/DRT0005c0071        TEMPE1 = (ZPL*ZPL/BPL2) * ( S0PL + XIN/(XIN2+BPL2) ) - S0PL0071     1          -(ZMI*ZMI/BMI2) * ( S0MI + XIN/(XIN2+BMI2) ) + S0MI0005c^^^^^^^^^^^^^^^^^^^^^^^^^^^^0075        TEMPE2 = 2.*ZPL*ZPL*S1PL/GPL1-S1PL-ZPL*ZPL*XI1*S0PL/(GPL1*BPL2)0061     1     + ZPL*ZPL*(1.-XI1*XIN/BPL2)/(GPL1*(XIN2+BPL2))0074     2     -2.*ZMI*ZMI*S1MI/GMI1 + S1MI + ZMI*ZMI*XI1*S0MI/(GMI1*BMI2)0061     3     - ZMI*ZMI*(1.-XI1*XIN/BMI2)/(GMI1*(XIN2+BMI2))0005c0071        TEMPE3 = 4.*ZPL*ZPL*S2PL/GPL2 + 2.*ZPL*ZPL*GP2IN/X2D - S2PL0074     1    + ZPL*ZPL*S0PL/(GPL2*BPL2)-2.*ZPL*ZPL*XI2/(GPL2*(XIN2+BPL2))0064     2    + ZPL*ZPL*XIN*(XI22-BPL2)*GP2IN/(BPL2*(XIN2+BPL2))0065     3    -4.*ZMI*ZMI*S2MI/GMI2 - 2.*ZMI*ZMI*GM2IN/X2D + S2MI0074     4    - ZMI*ZMI*S0MI/(GMI2*BMI2)+2.*ZMI*ZMI*XI2/(GMI2*(XIN2+BMI2))0064     5    - ZMI*ZMI*XIN*(XI22-BMI2)*GM2IN/(BMI2*(XIN2+BMI2))0005C0054        TEMPE4 = ZPL/(XIN2+BPL2) - ZMI/(XIN2+BMI2)0074        TEMPE5 = ZPL*(BPL2-XI12)*S0PL/(GPL1*GPL1)+2.*ZPL*XI1*PPL1/GPL10053     1         + ZPL*(XIN+XI1)/(GPL1*(XIN2+BPL2))0074     2         - ZMI*(BMI2-XI12)*S0MI/(GMI1*GMI1)-2.*ZMI*XI1*PMI1/GMI10053     3         - ZMI*(XIN+XI1)/(GMI1*(XIN2+BMI2))0005c0072        TEMPE6 = ZPL*(BPL2-3.*XI22)*PPL2/GPL2 - 2.*ZPL*XI2*GP2IN/X2D0062     1         + 2.*ZPL*XI2*(XI22-3.*BPL2)*S0PL*GP2IN/GPL20054     2         + ZPL*(BPL2-XI22)*GP2IN/(XIN2+BPL2)0053     3         - 2.*ZPL*XI2*XIN*GP2IN/(XIN2+BPL2)0072     4         - ZMI*(BMI2-3.*XI22)*PMI2/GMI2 + 2.*ZMI*XI2*GM2IN/X2D0062     5         - 2.*ZMI*XI2*(XI22-3.*BMI2)*S0MI*GM2IN/GMI20054     6         - ZMI*(BMI2-XI22)*GM2IN/(XIN2+BMI2)0053     7         + 2.*ZMI*XI2*XIN*GM2IN/(XIN2+BMI2)0005C0074        DBTAIL(8,1) = -0.5*FY*(PR(1)*TEMPE1+PR(2)*TEMPE2+PR(3)*TEMPE3)0074        DBTAIL(8,3) =  0.5*FY*(PR(1)*TEMPE4-PR(2)*TEMPE5-PR(3)*TEMPE6)0005c0005c0016C       d/dD0005C0061        TEMPF1 = - ZR * ( XIN/(XIN2 + BTA2) + S0 ) / BTA20070     1           + 0.5 * ZPL * ( XIN/(XIN2 + BPL2) + S0PL ) / BPL20070     2           + 0.5 * ZMI * ( XIN/(XIN2 + BMI2) + S0MI ) / BMI20005c0064        TEMPF2 = ZR * ( -2.*PX1/GX1 + XI1*(XI12+3.*BTA2)*S0/0072     1    (GX1*GX1*BTA2) + (XI1*XIN - BTA2)/(GX1*BTA2*(XIN2+BTA2)) )^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^0065     2   -0.5*ZPL* ( -2.*PPL1/GPL1 + XI1*(XI12+3.*BPL2)*S0PL/0073     3    (GPL1*GPL1*BPL2) + (XI1*XIN-BPL2)/(GPL1*BPL2*(XIN2+BPL2)) )0065     4   -0.5*ZMI* ( -2.*PMI1/GMI1 + XI1*(XI12+3.*BMI2)*S0MI/0073     5    (GMI1*GMI1*BMI2) + (XI1*XIN-BMI2)/(GMI1*BMI2*(XIN2+BMI2)) )0005c0074        TEMPF3 = ZR*( -4.*S2/GX2-S0/(GX2*BTA2)+2.*XI2*G2IN/(XIN2+BTA2)0069     1   - XIN*(XI22-BTA2)*G2IN/(BTA2*(XIN2+BTA2)) -2.*G2IN/X2D )0071     2  -0.5*ZPL*( -4.*S2PL/GPL2 - S0PL/(GPL2*BPL2) + 2.*XI2*GP2IN/0059     3             (XIN2+BPL2) - XIN*(XI22-BPL2)*GP2IN/0057     4             (BPL2*(XIN2+BPL2)) -2.*GP2IN/X2D )0071     5  -0.5*ZMI*( -4.*S2MI/GMI2 - S0MI/(GMI2*BMI2) + 2.*XI2*GM2IN/0059     6             (XIN2+BMI2) - XIN*(XI22-BMI2)*GM2IN/0057     7             (BMI2*(XIN2+BMI2)) -2.*GM2IN/X2D )0005C0039        TEMPF4 = 1. / (XIN2 + BTA2)0071     1       - (2.*XIN2+BPL2+BMI2) / ( 2.*(XIN2+BPL2)*(XIN2+BMI2) )0005c0065        TEMPF5 = ( S0 - 2.*G1 - (XIN+XI1)/(XIN2+BTA2) ) / GX10071     1     -0.5 * ( S0PL - 2.*G1PL - (XIN+XI1)/(XIN2+BPL2) ) / GPL10071     2     -0.5 * ( S0MI - 2.*G1MI - (XIN+XI1)/(XIN2+BMI2) ) / GMI10005c0074        TEMPF6 =  ( PX2 - 4.*G2/GX2 - 2.*XI2*S0*G2IN - 2.*XI2*G2IN/X2D0062     1   + ( 2.*XIN*XI2 - BTA2 + XI22 )*G2IN/(XIN2+BTA2) )0074     2   -( PPL2 - 4.*G2PL/GPL2 - 2.*XI2*S0PL*GP2IN - 2.*XI2*GP2IN/X2D0066     3   + (2.*XIN*XI2 - BPL2 + XI22)*GP2IN/(XIN2+BPL2) ) / 2.0074     4   -( PMI2 - 4.*G2MI/GMI2 - 2.*XI2*S0MI*GM2IN - 2.*XI2*GM2IN/X2D0066     5   + (2.*XIN*XI2 - BMI2 + XI22)*GM2IN/(XIN2+BMI2) ) / 2.0005C0038        DBTAIL(9,1) = FY*SQRT(D2)*0066     1                (PR(1)*TEMPF1+PR(2)*TEMPF2+PR(3)*TEMPF3)0038        DBTAIL(9,3) = FY*SQRT(D2)*0066     1                (PR(1)*TEMPF4+PR(2)*TEMPF5+PR(3)*TEMPF6)0005C0005C0022C       d/dDELTA-Y0005C0057        DBTAIL(10,1) = 2.*PI*FY*YR*YR*BTAIL(1)/PR(10)0057        DBTAIL(10,3) = 2.*PI*FY*YR*YR*BTAIL(3)/PR(10)0005C0018        RETURN0015        END0005C0005C0005c0057        SUBROUTINE NLRING(X,Y,Z,TILT,MA,BRING,DBRING)0005C0005c^^^^^^^^^^^^^^^^^^0074c       After finding a "bug" in this subroutine (improper computation0074c       of d/dRrc terms), I decided to re-ro the whole subroutine from0066c       scratch (instead of modifying David's LSRING routine).00040078c       Notice that parameters are passed as arguments as in NLTAIL above.0077c       In addition, keep in mind that ring contribution must be computed0077c       in SM coordinates. Transformations have been incorporated here in0074c       the subroutine rather than with calls to ROTATE since they are0043c       relatively simple to implement.0005c0081c       Input are X,Y,Z,TILT specifying the 4-dim position at which the field0081c       is to be computed, and MA, the number of parameters in the prescribed0050c       representation for the magnetic field.0080c       Output are BRING, a 3-component array containing the three cartesian0080c       components of the ring current contribution to the external magnetic0083c       field, and DBRING, an MA by 3 matrix containing the partial derivatives0076c       of each of the 3 field components wrt each of the MA parameters.0005c0005c0037          IMPLICIT REAL (A-H,O-Z)0005C0031        COMMON/PARAM/PR(40)0043        DIMENSION BRING(3),DBRING(50,3)0005C0005c0005c0064c       Zero all derivatives so that only those compute here0029c       will be non-zero.0005c0024        DO 10 I=1,MA0026           DO 10 J=1,30034              DBRING(I,J) = 0.002010      CONTINUE0005C0005C0028        PI = 3.1415926540034        S6 = SIN(TILT*PI/180.)0034        C6 = COS(TILT*PI/180.)0005C0046C       CONVERT INPUT X,Y,Z FROM GSM TO SM0005C0029        XSM = X*C6 - Z*S60019        YSM = Y0029        ZSM = X*S6 + Z*C60005C0005c0040        BRC = PR(13)            !Brc0061        RRC = ABS(PR(14))       !FORCE Rrc TO BE POSITIVE0026        RRC2 = RRC*RRC0027        RRC3 = RRC2*RRC0005C0061        RHO2 = (XSM*XSM + YSM*YSM)      !Notice no /Rrc^20050        DEN1 = (RHO2 + ZSM*ZSM + 4.*RRC2)**2.50050        DEN2 = (RHO2 + ZSM*ZSM + 4.*RRC2)**3.50005C^^^^^^^^^^^^^^^^^^^^^^^^^^^^0042        TERM1 = (12.* ZSM * RRC3)/DEN10062        TERM2 = 4.*RRC3*(2.*ZSM*ZSM - RHO2 + 8.*RRC2)/DEN10005C0046c       Additional terms needed for D/DRrc0005c0066        TERM3 = 12.*ZSM*RRC2*(3.*RHO2+3.*ZSM*ZSM-8.*RRC2)/DEN20043        TERM4 = 2.*ZSM*ZSM-RHO2+8.*RRC20044        TERM5 = RHO2 + ZSM*ZSM + 4.*RRC20071        TERM6 = (3.*TERM4*TERM5+16.*RRC2*TERM5-20.*RRC2*TERM4)*RRC20005C0068C       Compute field components divided by Brc, that way we get0068c       derivative terms immediately, and for full field we just0028c       multiply by Brc.0005c0028        BXSM = TERM1*XSM0028        BYSM = TERM1*YSM0024        BZSM = TERM20005C0064C       Transform to GSM coordinates -- assignments first on0059c       derivatives, then multiply by Brc to get field.0005c0067        DBRING(11,1) = BXSM*C6 + BZSM*S6        !DBx/DBrc (GSM)0067        DBRING(11,2) = BYSM                     !DBy/DBrc (GSM)0067        DBRING(11,3) = -BXSM*S6 + BZSM*C6       !DBz/DBrc (GSM)0005C0045        DBRING(12,1) = BRC*XSM*TERM3*C6 +0067     1                 4.*BRC*TERM6*S6/DEN2     !DBx/DRrc (GSM)0067        DBRING(12,2) = BRC*YSM*TERM3            !DBy/DRrc (GSM)0046        DBRING(12,3) = -BRC*XSM*TERM3*S6 +0067     1                  4.*BRC*TERM6*C6/DEN2    !DBz/DRrc (GSM)0005c0039        BRING(1) = BRC*DBRING(11,1)0039        BRING(2) = BRC*DBRING(11,2)0039        BRING(3) = BRC*DBRING(11,3)0005C0005c0018        RETURN0015        END0005C0005C0005c0063        SUBROUTINE NLMPAUSE(X,Y,Z,TILT,MA,BMPAUSE,DBMPAUSE)0005C0060C       This subroutine was adapted from scratch.  Input0075c       are X,Y,Z,TILT specifying the 4-dim position at which the field0081c       is to be computed, and MA, the number of parameters in the prescribed0050c       representation for the magnetic field.0082c       Output are BMPAUSE, a 3-component array containing the three cartesian0080c       components of the magnetopause/Birkeland current contribution to the0083c       external magnetic field, and DBMPAUSE, an MA by 3 matrix containing the^^^^^^^^^^^^^^^^^^^^0081c       partial derivatives of each of the 3 field components wrt each of the0026c       MA parameters.0005c0005c0037          IMPLICIT REAL (A-H,O-Z)0005C0031        COMMON/PARAM/PR(40)0047        DIMENSION BMPAUSE(3),DBMPAUSE(50,3)0005C0077c       Initialize all derivatives to zero so that only the ones that are0063c       computed inside this routine may get nonzero values0005c0065        DO 10 I=1,MA                    !Loop over parameters0075           DO 10 J=1,3                  !Loop over components of vector0036              DBMPAUSE(I,J) = 0.002010      CONTINUE0005C0028        PI = 3.1415926540034        S6 = SIN(TILT*PI/180.)0034        C6 = COS(TILT*PI/180.)0005C0052        E1 = EXP(X/PR(15))      !EXP(X/DELTA-X1)0071        E2 = EXP(2.*X/PR(15))   !EXP(2X/DELTA-X1) = EXP(X/DELTA-X2)0005C0075c       Since I have carried out the derivation in paper by components,0073c       it is simplest to follow the same format here. So rather than0072c       compute all intermediate factors and define output variables0070c       at the end, the process will be carried out by components.0005c0043C       Evaluate Bx and its derivatives0005c0040        T1 = PR(16)*Z*C6 + PR(17)*S60072        T2 = PR(18)*Z*C6 + ( PR(19) + PR(20)*Y*Y + PR(21)*Z*Z ) * S60070        T3 = X / ( PR(15)*PR(15) )              ! x / (delta-x1)^20005C0064        BMPAUSE(1) = E1*T1 + E2*T2              !X-COMPONENT0005C0043C       Now for the nonzero derivatives0005c0074        DBMPAUSE(13,1) = - T3 * (E1*T1 + 2.*E2*T2)      !dBx/dDELTA-X10068        IF ( PR(16).NE.0.) DBMPAUSE(14,1) = E1*Z*C6     !dBx/da10068        IF ( PR(17).NE.0.) DBMPAUSE(15,1) = E1*S6       !dBx/da20068        DBMPAUSE(16,1) = E2*Z*C6                        !dBx/da30068        DBMPAUSE(17,1) = E2*S6                          !dBx/da40068        DBMPAUSE(18,1) = E2*Y*Y*S6                      !dBx/da50068        DBMPAUSE(19,1) = E2*Z*Z*S6                      !dBx/da60005C0005c0045c       Move on to By and its derivatives0005c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^0056        T4 = - ( PR(16)/PR(15)  + 2.*PR(30) )*Y*Z*C60050        T5 = - ( PR(17)/PR(15) + PR(31) )*Y*S60055        T6 = -2.* ( PR(18)/PR(15) + PR(34) )*Y*Z*C60053        T7 = - ( 2.*PR(19)/PR(15) + PR(35) )*Y*S60059        T8 = - (2.*PR(20)/PR(15) + PR(36) )*Y*Y*Y*S6/3.0059        T9 = - (2.*PR(21)/PR(15) + 3.*PR(37) )*Y*Z*Z*S60005c0030        S1 = X/PR(15) + 1.0033        S2 = 2.*X/PR(15) + 1.0051        S3 = ( PR(16)*S1 + 2.*PR(30)*X )*Y*Z*C60046        S4 = ( PR(17)*S1 + PR(31)*X )*Y*S60054        S5 = ( 2.*PR(18)*S2 + 4.*PR(34)*X )*Y*Z*C60005c0052c       Skip S6 since already used for sin(tilt)0005c0051        S7 = ( 2.*PR(19)*S2*Y + 2.*PR(35)*X*Y +0062     1      2.*PR(20)*S2*Y*Y*Y/3. + 2.*PR(36)*X*Y*Y*Y/3. +0061     2      2.*PR(21)*S2*Y*Z*Z + 6.*PR(37)*X*Y*Z*Z ) * S60005C0005C0072        BMPAUSE(2) = E1*(T4+T5) + E2*(T6+T7+T8+T9)      !Y-COMPONENT0005C0049C       Now computing the nonzero derivatives0005c0053        DBMPAUSE(13,2) = (E1*(S3+S4)+E2*(S5+S7))/0082     1                   (PR(15)*PR(15))                        !dBy/ddelta-x10076        IF ( PR(16).NE.0.) DBMPAUSE(14,2) = -E1*Y*Z*C6/PR(15)   !dBy/da10076        IF ( PR(17).NE.0.) DBMPAUSE(15,2) = -E1*Y*S6/PR(15)     !dBy/da20076        DBMPAUSE(16,2) = -2.*E2*Y*Z*C6/PR(15)                   !dBy/da30076        DBMPAUSE(17,2) = -2.*E2*Y*S6/PR(15)                     !dBy/da40076        DBMPAUSE(18,2) = -2.*E2*Y*Y*Y*S6/(3.*PR(15))            !dBy/da50076        DBMPAUSE(19,2) = -2.*E2*Y*Z*Z*S6/PR(15)                 !dBy/da60076        IF ( PR(30).NE.0.) DBMPAUSE(22,2) = -2.*E1*Y*Z*C6       !dBy/dc30076        IF ( PR(31).NE.0.) DBMPAUSE(23,2) = -E1*Y*S6            !dBy/dc40076        DBMPAUSE(26,2) = -2.*E2*Y*Z*C6                          !dBy/dc70076        DBMPAUSE(27,2) = -E2*Y*S6                               !dBy/dc80076        DBMPAUSE(28,2) = -E2*Y*Y*Y*S6/3.                        !dBy/dc90077        DBMPAUSE(29,2) = -3.*E2*Y*Z*Z*S6                        !dby/dc100005C0005C0042C       Next, proceed with Z-component0005c^^^^^^^^^^^^^^^0049        R1 = PR(28) + PR(29)*Y*Y + PR(30)*Z*Z0049        R2 = PR(32) + PR(33)*Y*Y + PR(34)*Z*Z0055        R3 = PR(35)*Z + PR(36)*Z*Y*Y + PR(37)*Z*Z*Z0005C0080        BMPAUSE(3) = E1*(R1*C6+PR(31)*Z*S6) + E2*(R2*C6+R3*S6)  !Z-component0005c0031C       NOW THE DERIVATIVES0005C0059        DBMPAUSE(13,3) = -T3*( E1*(R1*C6 + PR(31)*Z*S6)0074     1                        + 2.*E2*(R2*C6+R3*S6) )   !dBz/ddelta-x10068        IF ( PR(28).NE.0.) DBMPAUSE(20,3) = E1*C6       !dBz/dc10068        IF ( PR(29).NE.0.) DBMPAUSE(21,3) = E1*Y*Y*C6   !dBz/dc20068        IF ( PR(30).NE.0.) DBMPAUSE(22,3) = E1*Z*Z*C6   !dBz/dc30068        IF ( PR(31).NE.0.) DBMPAUSE(23,3) = E1*Z*S6     !dBz/dc40068        DBMPAUSE(24,3) = E2*C6                          !dBz/dc50068        DBMPAUSE(25,3) = E2*Y*Y*C6                      !dBz/dc60068        DBMPAUSE(26,3) = E2*Z*Z*C6                      !dBz/dc70068        DBMPAUSE(27,3) = E2*Z*S6                        !dBz/dc80068        DBMPAUSE(28,3) = E2*Z*Y*Y*S6                    !dBz/dc90069        DBMPAUSE(29,3) = E2*Z*Z*Z*S6                    !dBz/dc100005c0018        RETURN0015        END00040073c  ******************************************************************0005c0005c0005C0051      SUBROUTINE EXT89M(IOPT,PS,X,Y,Z,BX,BY,BZ)0005C0074C  CALCULATES GSM COMPONENTS OF THE EXTERNAL MAGNETIC FIELD BY CALLING0038C  THE SUBROUTINE T89M (see below)0005C0024C  INPUT PARAMETERS:0078c-------------------------------------------------------------------------0071c     IOPT   defines the Kp-index interval, according to the table:0005c0071c    IOPT= 1       2         3         4         5        6       70005c0073c    KP=  0,0+  1-,1,1+   2-,2,2+   3-,3,3+   4-,4,4+  5-,5,5+   >=6-0005c0079c--------------------------------------------------------------------------0048c     PS is the dipole tilt angle in radians0079c--------------------------------------------------------------------------0079c     X, Y, Z   are the GSM coordinates of the point (in the Earth's radii)^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^0079c--------------------------------------------------------------------------0005c0025c  OUTPUT PARAMETERS:0077C------------------------------------------------------------------------0069C    BX, BY, BZ   are the external field components, in nanotesla0079c--------------------------------------------------------------------------0005c0073c   THE COEFFICIENTS BELOW CORRESPOND TO THE NEW VERSION OF THE MODEL0075C    BASED ON THE EXTENDED DATASETS (INCLUDING 4 YEARS OF ISEE-1 AND -20016C    DATA). 0067C   NOTE ALSO THAT THE INITIAL DATA WERE ABERRATED BY 4 DEGREES0075C    IN THE GSE SYSTEM, BEFORE USING THEM FOR THE LEAST-SQUARES FITTING0075C    THE MODEL PARAMETERS.  THEREFORE, X, Y, AND Z SHOULD BE CONSIDERED0041C    AS THE ABERRATED GSM COORDINATES0075C======================================================================0005C0022C    REFERENCES:  0074C     1.  Tsyganenko, N.A., A magnetospheric magnetic field model with0078c         the warped tail current sheet. Planet.Space Sci., v.37, pp.5-20,0019c         1989.0075c     2.  Peredo, M., D.P. Stern, and N.A. Tsyganenko, Are the existing0074c         magnetospheric magnetic field models excessively stretched ?0058c         J.Geophys.Res., v.98, 1993, to be published.0005C0005C0053       DIMENSION XI(4),F(3),DER(3,30),PARAM(30,7)0029       COMMON /PAR/ A(30)0033       DOUBLE PRECISION F,DER0075        DATA PARAM/-116.53,-10719.,42.375,59.753,-11363.,1.7844,30.268,0069     * -0.35372E-01,-0.66832E-01,0.16456E-01,-1.3024,0.16529E-02,0072     * 0.20293E-02,20.289,-0.25203E-01,224.91,-9234.8,22.788,7.8813,0064     * 1.8362,-0.27228,8.8184,2.8714,14.468,32.177,0.01,0.0,0065     * 7.0459,4.0,20.0,-55.553,-13198.,60.647,61.072,-16064.,0071     * 2.2534,34.407,-0.38887E-01,-0.94571E-01,0.27154E-01,-1.3901,0070     * 0.13460E-02,0.13238E-02,23.005,-0.30565E-01,55.047,-3875.7,0068     * 20.178,7.9693,1.4575,0.89471,9.4039,3.5215,14.474,36.555,0074     * 0.01,0.0,7.0787,4.0,20.0,-101.34,-13480.,111.35,12.386,-24699.,^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^0067     * 2.6459,38.948,-0.34080E-01,-0.12404,0.29702E-01,-1.4052,0076     * 0.12103E-02,0.16381E-02,24.49,-0.37705E-01,-298.32,4400.9,18.692,0076     * 7.9064,1.3047,2.4541,9.7012,7.1624,14.288,33.822,0.01,0.0,6.7442,0073     * 4.0,20.0,-181.69,-12320.,173.79,-96.664,-39051.,3.2633,44.968,0074     * -0.46377E-01,-0.16686,0.048298,-1.5473,0.10277E-02,0.31632E-02,0074     * 27.341,-0.50655E-01,-514.10,12482.,16.257,8.5834,1.0194,3.6148,0072     * 8.6042,5.5057,13.778,32.373,0.01,0.0,7.3195,4.0,20.0,-436.54,0069     * -9001.0,323.66,-410.08,-50340.,3.9932,58.524,-0.38519E-01,0072     * -0.26822,0.74528E-01,-1.4268,-0.10985E-02,0.96613E-02,27.557,0075     * -0.56522E-01,-867.03,20652.,14.101,8.3501,0.72996,3.8149,9.2908,0074     *  6.4674,13.729,28.353,0.01,0.0,7.4237,4.0,20.0,-707.77,-4471.9,0070     * 432.81,-435.51,-60400.,4.6229,68.178,-0.88245E-01,-0.21002,0071     * 0.11846,-2.6711,0.22305E-02,0.10910E-01,27.547,-0.54080E-01,0075     * -424.23,1100.2,13.954,7.5337,0.89714,3.7813,8.2945,5.174,14.213,0073     * 25.237,0.01,0.0,7.0037,4.0,20.0,-1190.4,2749.9,742.56,-1110.3,0071     * -77193.,7.6727,102.05,-0.96015E-01,-0.74507,0.11214,-1.3614,0070     * 0.15157E-02,0.22283E-01,23.164,-0.74146E-01,-2219.1,48253.,0073     * 12.714,7.6777,0.57138,2.9633,9.3909,9.7263,11.123,21.558,0.01,0031     * 0.0,4.4518,4.0,20.0/0005c0023	DATA ID/1/,IOP/10/0005c0026	IF (IOP.NE.IOPT) then0005c0015	  IOP=IOPT0013	  ID = 10005c0026          DO 10 I=1,300025	  A(I)=PARAM(I,IOPT)001710	  continue0005C0010	ENDIF0005c0019        XI(1)=X0019        XI(2)=Y0019        XI(3)=Z0020        XI(4)=PS0037         CALL T89M(ID,A,XI,F,DER)0023	 IF (ID.EQ.1) ID=20019        BX=F(1)0019        BY=F(2)0019        BZ=F(3)0018        RETURN0014       END0005c0072C-------------------------------------------------------------------0005c0050          SUBROUTINE  T89M (ID, A, XI, F, DER)0005C0056C        ***  N.A. Tsyganenko ***  8-10.12.1991  ***0056C        ***     (REVISED APRIL 1992)            ***0005C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^0065C      Calculates dependent model variables and their deriva-0062C  tives for given independent variables and model parame-0066C  ters.  Specifies model functions with free parameters which0061C  must be determined by means of least squares fits (RMS0031C  minimization procedure).0005C0037C      Description of parameters:0005C0081C  ID  - number of the data point in a set (initial assignments are performed0049c        only for ID=1, saving thus CPU time)0054C  A   - input vector containing model parameters;0059C  XI  - input vector containing independent variables;0054C  F   - output double precision vector containing0054C        calculated values of dependent variables;0056C  DER   - output double precision vector containing0059C        calculated values for derivatives of dependent0055C        variables with respect to model parameters0005C0063C - - - - - - - - - - - - - - - - - - - - - - - - - - - - -0005C0065C      T89M represents external magnetospheric magnetic field0070C  in Cartesian SOLAR MAGNETOSPHERIC coordinates (Tsyganenko N.A.,0078C  Planet. Space Sci., 1989, v.37, p.5-20; the "T89 model" with the warped0075c  tail current sheet) + A MODIFICATION ADDED IN APRIL 1992 (SEE BELOW)0005C0076C      Model formulae for the magnetic field components contain in total0066c  30 free parameters (17 linear and 13 nonlinear parameters).0082C      First 2 independent linear parameters A(1)-A(2) correspond to contribu-0081c  tion from the tail current system, then follow A(3) and A(4) which are the0079c  amplitudes of symmetric and antisymmetric terms in the contribution from0080c  the closure currents; A(5) is the ring current amplitude. Then follow the0083c coefficients A(6)-A(15) which define Chapman-Ferraro+Birkeland current field.0074c    The coefficients c16-c19  (see Formula 20 in the original paper),0080c   due to DivB=0 condition, are expressed thru A(6)-A(15) and hence are not0026c    independent ones.0080c  A(16) AND A(17) CORRESPOND TO THE TERMS WHICH YIELD THE TILT ANGLE DEPEN-^^^^^^^^^^^^^^^^^^^^^^0082C    DENCE OF THE TAIL CURRENT INTENSITY (ADDED ON APRIL 9, 1992; SEE COMMENTS0038C   IN THE BROWN GODDARD NOTEBOOK)0005C0032C      Nonlinear parameters:0005C0081C    A(18) : DX - Characteristic scale of the Chapman-Ferraro field along the0019c        X-axis0070C    A(19) : ADR (aRC) - Characteristic radius of the ring current0068c    A(20) : D0 - Basic half-thickness of the tail current sheet0079C    A(21) : DD (GamRC)- defines rate of thickening of the ring current, as0046c             we go from night- to dayside0075C    A(22) : Rc - an analog of "hinging distance" entering formula (11)0073C    A(23) : G - amplitude of tail current warping in the Y-direction0063C    A(24) : aT - Characteristic radius of the tail current0079c    A(25) : Dy - characteristic scale distance in the Y direction entering0039c                 in W(x,y) in (13)0081c    A(26) : Delta - defines the rate of thickening of the tail current sheet0070c                 in the Y-direction (in T89 it was fixed at 0.01)0079c    A(27) : Q - this parameter was fixed at 0 in the final version of T89;0075c              initially it was introduced for making Dy to depend on X0054c    A(28) : Sx (Xo) - enters in W(x,y) ; see (13)0081c    A(29) : Gam (GamT) - enters in DT in (13) and defines rate of tail sheet0081c              thickening on going from night to dayside; in T89 fixed at 4.00080c    A(30) : Dyc - the Dy parameter for closure current system; in T89 fixed0027c               at 20.00064c  - - - - - - - - - - - - - - - - - - - - - - - - - - - - -0005C0040	 IMPLICIT  REAL * 8  (A - H, O - Z)0005C0030         REAL  A(1), XI(1)0005C0032	 DIMENSION  F(3), DER(3,30)0005C0031         INTEGER   ID, I, L0072      DATA A02,XLW2,YN,RPI,RT/25.D0,170.D0,30.D0,0.31830989D0,30.D0/0034      DATA XD,XLD2/0.D0,40.D0/0005C0080C   The last two quantities define variation of tail sheet thickness along X0005C0036      DATA SXC,XLWC2/4.D0,50.D0/0005C0080C   The two quantities belong to the function WC which confines tail closure^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^0039c    current in X- and Y- direction0005C0030      DATA IK,DXL/0,20.D0/0005C0005C0034         IF (ID.NE.1)  GOTO  30024	   DO  2  I = 1, 300025	     DO  1  L = 1, 30035  1            DER(L,I) = 0.0D00023  2        CONTINUE0005C0016        IK=10005C0020       DYC=A(30)0022       DYC2=DYC**20019       DX=A(18)0025       HA02=0.5D0*A020028       RDX2M=-1.D0/DX**20022       RDX2=-RDX2M0026       RDYC2=1.D0/DYC20030       HLWC2M=-0.5D0*XLWC20029       DRDYC2=-2.D0*RDYC20041       DRDYC3=2.D0*RDYC2*DSQRT(RDYC2)0029       HXLW2M=-0.5D0*XLW20020       ADR=A(19)0019       D0=A(20)0019       DD=A(21)0019       RC=A(22)0018       G=A(23)0019       AT=A(24)0016       DT=D00020       DEL=A(26)0018       P=A(25)0018       Q=A(27)0019       SX=A(28)0020       GAM=A(29)0029       HXLD2M=-0.5D0*XLD20020       ADSL=0.D00020       XGHS=0.D00017       H=0.D00018       HS=0.D00020       GAMH=0.D00023       W1=-0.5D0/DX0026       DBLDEL=2.D0*DEL0021       W2=W1*2.D00024       W4=-1.D0/3.D00019       W3=W4/DX0020       W5=-0.5D00019       W6=-3.D00019       AK1=A(1)0019       AK2=A(2)0019       AK3=A(3)0019       AK4=A(4)0019       AK5=A(5)0019       AK6=A(6)0019       AK7=A(7)0019       AK8=A(8)0019       AK9=A(9)0021       AK10=A(10)0021       AK11=A(11)0021       AK12=A(12)0021       AK13=A(13)0021       AK14=A(14)0021       AK15=A(15)0022        AK16=A(16)0022        AK17=A(17)0019       SXA=0.D00019       SYA=0.D00019       SZA=0.D00031       AK610=AK6*W1+AK10*W50028       AK711=AK7*W2-AK110031       AK812=AK8*W2+AK12*W60031       AK913=AK9*W3+AK13*W40024       RDXL=1.D0/DXL0027       HRDXL=0.5D0*RDXL0024       A6H=AK6*0.5D00023       A9T=AK9/3.D00027       YNP=RPI/YN*0.5D00022       YND=2.D0*YN0005C0021  3      CONTINUE0005C0025           X  = XI(1)0018	   Y  = XI(2)0018	   Z  = XI(3)0032               TLT2=XI(4)**20024	   SPS = sin(XI(4))0045           CPS = DSQRT (1.0D0 - SPS ** 2)0005C0017       X2=X*X0017       Y2=Y*Y0017       Z2=Z*Z0022       TPS=SPS/CPS0024       HTP=TPS*0.5D00020       GSP=G*SPS^^^^^^^^^^^^^^^^^^^^^^^^0026       XSM=X*CPS-Z*SPS0026       ZSM=X*SPS+Z*CPS0005C0078C   CALCULATE THE FUNCTION ZS DEFINING THE SHAPE OF THE TAIL CURRENT SHEET0037C    AND ITS SPATIAL DERIVATIVES:0005C0021       XRC=XSM+RC0029       XRC16=XRC**2+16.D00028       SXRC=DSQRT(XRC16)0019       Y4=Y2*Y20023       Y410=Y4+1.D40023       SY4=SPS/Y4100021       GSY4=G*SY40029       ZS1=HTP*(XRC-SXRC)0025       DZSX=-ZS1/SXRC0025       ZS=ZS1-GSY4*Y40037       D2ZSGY=-SY4/Y410*4.D4*Y2*Y0024       DZSY=G*D2ZSGY0005C0066C   CALCULATE THE COMPONENTS OF THE RING CURRENT CONTRIBUTION:0005C0022       XSM2=XSM**20031       DSQT=DSQRT(XSM2+A02)0036       FA0=0.5D0*(1.D0+XSM/DSQT)0024       DDR=D0+DD*FA00028       DFA0=HA02/DSQT**30020       ZR=ZSM-ZS0033       TR=DSQRT(ZR**2+DDR**2)0022       RTR=1.D0/TR0022       RO2=XSM2+Y20022       ADRT=ADR+TR0024       ADRT2=ADRT**20030       FK=1.D0/(ADRT2+RO2)0025       DSFC=DSQRT(FK)0024       FC=FK**2*DSFC0034       FACXY=3.0D0*ADRT*FC*RTR0021       XZR=XSM*ZR0019       YZR=Y*ZR0026       DBXDP=FACXY*XZR0029       DER(2,5)=FACXY*YZR0031       XZYZ=XSM*DZSX+Y*DZSY0038       FAQ=ZR*XZYZ-DDR*DD*DFA0*XSM0046       DBZDP=FC*(2.D0*ADRT2-RO2)+FACXY*FAQ0039       DER(1,5)=DBXDP*CPS+DBZDP*SPS0039       DER(3,5)=DBZDP*CPS-DBXDP*SPS0005C0053C  CALCULATE THE TAIL CURRENT SHEET CONTRIBUTION:0005C0023       DELY2=DEL*Y20021       D=DT+DELY20041       IF (DABS(GAM).LT.1.D-6) GOTO 80021       XXD=XSM-XD0033       RQD=1.D0/(XXD**2+XLD2)0026       RQDS=DSQRT(RQD)0034       H=0.5D0*(1.D0+XXD*RQDS)0030       HS=-HXLD2M*RQD*RQDS0021       GAMH=GAM*H0019       D=D+GAMH0026       XGHS=XSM*GAM*HS0023       ADSL=-D*XGHS0018   8   D2=D**20028       T=DSQRT(ZR**2+D2)0022       XSMX=XSM-SX0036       RDSQ2=1.D0/(XSMX**2+XLW2)0028       RDSQ=DSQRT(RDSQ2)0035       V=0.5D0*(1.D0-XSMX*RDSQ)0032       DVX=HXLW2M*RDSQ*RDSQ20042       OM=DSQRT(DSQRT(XSM2+16.D0)-XSM)0036       OMS=-OM/(OM*OM+XSM)*0.5D00028       RDY=1.D0/(P+Q*OM)0021       OMSV=OMS*V0022       RDY2=RDY**20033       FY=1.D0/(1.D0+Y2*RDY2)0017       W=V*FY^^^^^^^^^^^^^^^^^^^^^^^^^^^0031       YFY1=2.D0*FY*Y2*RDY20024       FYPR=YFY1*RDY0023       FYDY=FYPR*FY0033       DWX=DVX*FY+FYDY*Q*OMSV0026       YDWY=-V*YFY1*FY0023       DDY=DBLDEL*Y0019       ATT=AT+T0031       S1=DSQRT(ATT**2+RO2)0021       F5=1.D0/S10027       F7=1.D0/(S1+ATT)0019       F1=F5*F70019       F3=F5**30020       F9=ATT*F30034       FS=ZR*XZYZ-D*Y*DDY+ADSL0028       XDWX=XSM*DWX+YDWY0021       RTT=1.D0/T0019       WT=W*RTT0022       BRRZ1=WT*F10022       BRRZ2=WT*F30026       DBXC1=BRRZ1*XZR0026       DBXC2=BRRZ2*XZR0029       DER(2,1)=BRRZ1*YZR0029       DER(2,2)=BRRZ2*YZR0037          DER(2,16)=DER(2,1)*TLT20037          DER(2,17)=DER(2,2)*TLT20021       WTFS=WT*FS0037       DBZC1=W*F5+XDWX*F7+WTFS*F10037       DBZC2=W*F9+XDWX*F1+WTFS*F30039       DER(1,1)=DBXC1*CPS+DBZC1*SPS0039       DER(1,2)=DBXC2*CPS+DBZC2*SPS0039       DER(3,1)=DBZC1*CPS-DBXC1*SPS0039       DER(3,2)=DBZC2*CPS-DBXC2*SPS0037          DER(1,16)=DER(1,1)*TLT20037          DER(1,17)=DER(1,2)*TLT20037          DER(3,16)=DER(3,1)*TLT20037          DER(3,17)=DER(3,2)*TLT20005C0055C  CALCULATE CONTRIBUTION FROM THE CLOSURE CURRENTS0005C0019       ZPL=Z+RT0019       ZMN=Z-RT0023       ROGSM2=X2+Y20035       SPL=DSQRT(ZPL**2+ROGSM2)0035       SMN=DSQRT(ZMN**2+ROGSM2)0021       XSXC=X-SXC0036       RQC2=1.D0/(XSXC**2+XLWC2)0026       RQC=DSQRT(RQC2)0035       FYC=1.D0/(1.D0+Y2*RDYC2)0039       WC=0.5D0*(1.D0-XSXC*RQC)*FYC0035       DWCX=HLWC2M*RQC2*RQC*FYC0031       DWCY=DRDYC2*WC*FYC*Y0030       SZRP=1.D0/(SPL+ZPL)0030       SZRM=1.D0/(SMN-ZMN)0029       XYWC=X*DWCX+Y*DWCY0022       WCSP=WC/SPL0022       WCSM=WC/SMN0025       FXYP=WCSP*SZRP0025       FXYM=WCSM*SZRM0022       FXPL=X*FXYP0023       FXMN=-X*FXYM0022       FYPL=Y*FXYP0023       FYMN=-Y*FXYM0030       FZPL=WCSP+XYWC*SZRP0030       FZMN=WCSM+XYWC*SZRM0029       DER(1,3)=FXPL+FXMN0035       DER(1,4)=(FXPL-FXMN)*SPS0029       DER(2,3)=FYPL+FYMN0035       DER(2,4)=(FYPL-FYMN)*SPS0029       DER(3,3)=FZPL+FZMN0035       DER(3,4)=(FZPL-FZMN)*SPS0005C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^0075C   NOW CALCULATE CONTRIBUTION FROM CHAPMAN-FERRARO SOURCES + ALL OTHER0005C0028           EX=DEXP(X/DX)0024           EC=EX*CPS0024           ES=EX*SPS0023           ECZ=EC*Z0023           ESZ=ES*Z0020	   ESZY2=ESZ*Y20020	   ESZZ2=ESZ*Z20018	   ECZ2=ECZ*Z0016	   ESY=ES*Y0005C0020	   DER(1,6)=ECZ0019	   DER(1,7)=ES0022	   DER(1,8)=ESY*Y0022	   DER(1,9)=ESZ*Z0023	   DER(2,10)=ECZ*Y0021	   DER(2,11)=ESY0024	   DER(2,12)=ESY*Y20024	   DER(2,13)=ESY*Z20020	   DER(3,14)=EC0023	   DER(3,15)=EC*Y20024	   DER(3,6)=ECZ2*W10025	   DER(3,10)=ECZ2*W50023	   DER(3,7)=ESZ*W20022	   DER(3,11)=-ESZ0025	   DER(3,8)=ESZY2*W20026	   DER(3,12)=ESZY2*W60025	   DER(3,9)=ESZZ2*W30026	   DER(3,13)=ESZZ2*W40005C0065C  FINALLY, CALCULATE NET EXTERNAL MAGNETIC FIELD COMPONENTS,0048C    BUT FIRST OF ALL THOSE FOR C.-F. FIELD:0005C0065      SX1=AK6*DER(1,6)+AK7*DER(1,7)+AK8*DER(1,8)+AK9*DER(1,9)0073      SY1=AK10*DER(2,10)+AK11*DER(2,11)+AK12*DER(2,12)+AK13*DER(2,13)0070      SZ1=AK14*DER(3,14)+AK15*DER(3,15)+AK610*ECZ2+AK711*ESZ+AK8120029     * *ESZY2+AK913*ESZZ20041       BXCL=AK3*DER(1,3)+AK4*DER(1,4)0041       BYCL=AK3*DER(2,3)+AK4*DER(2,4)0041       BZCL=AK3*DER(3,3)+AK4*DER(3,4)0076       BXT=AK1*DER(1,1)+AK2*DER(1,2)+BXCL +AK16*DER(1,16)+AK17*DER(1,17)0076       BYT=AK1*DER(2,1)+AK2*DER(2,2)+BYCL +AK16*DER(2,16)+AK17*DER(2,17)0076       BZT=AK1*DER(3,1)+AK2*DER(3,2)+BZCL +AK16*DER(3,16)+AK17*DER(3,17)0040       F(1)=BXT+AK5*DER(1,5)+SX1+SXA0040       F(2)=BYT+AK5*DER(2,5)+SY1+SYA0040       F(3)=BZT+AK5*DER(3,5)+SZ1+SZA0005C0017       RETURN0014       END^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^