000400040074A SET OF FORTRAN SUBROUTINES FOR COMPUTATIONS OF THE GEOMAGNETIC FIELD0050                  IN THE EARTH'S MAGNETOSPHERE00040037                        (GEOPACK)00040040                     N.A. TSYGANENKO0045                     INSTITUTE OF PHYSICS0053                     UNIVERSITY OF ST.-PETERSBURG0046                     STARY PETERGOF 1989040039                     ST.-PETERSBURG0031                     RUSSIA000400040074* Note: Modified by M. Peredo to include new models, and subroutines *0074* used by the ISTP/Science Planning and Operations Facility.         *0074* Descriptions for the new routines are appended at the end of this  *0074* files                                                              *000400040043                           INTRODUCTION00040075     Recent studies in the solar-terrestrial physics led to recognizing0076the role of the geomagnetic field as one of the most important characte-0077ristics of human environment. The Earth's magnetic field links the inter-0075planetary medium  with the upper atmosphere and the ionosphere,  guides0075the energetic charged particles  ejected during solar flares,  channels0075the low-frequency electromagnetic waves and heat flux, confines the ra-0076diation belt and auroral plasma particles, and serves as a giant accumu-0073lator of the solar wind energy  that eventually dissipates during the0075magnetic storms.  Investigations of these phenomena are closely related0074to the problem of forecasting the "weather" in the near-Earth space,  0077conditions of the radiowave propagation,  safety of manned space flights,0052reliability of communication systems, and so on.00040075    For this reason, it is of great importance to have tools for compu-0075ting the structure of the geomagnetic field, capable of taking into ac-0076count the magnetospheric disturbance and other factors defining variati-0076ons of the near-Earth electric current systems. In many applications it ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^0075is often necessary to evaluate the components of the geomagnetic field 0076vector in a wide range of distances, trace the field lines far away from0074the Earth's surface,  calculate the geomagnetically conjugate points, 0073and map a spacecraft position with respect to characteristic magneto-0075spheric/ionospheric boundaries.   This requires to use quantitative mo-0073dels of the Earth's internal field as well as those for the external 0066field produced by the magnetospheric electric current systems.00040076    The present set of subroutines computes the geomagnetic field compo-0074nents at any point of space within the Earth's magnetosphere up to the0078Moon's orbit and traces the force lines starting at a given initial point.0078Upon specifying the universal time and day of year as input parameters, it0077automatically performs all the necessary rotations of coordinate axes and0074takes into account the tilt angle of the Earth's magnetic dipole axis.0076The Earth's internal magnetic field can be computed either in the dipole0077approximation or by taking into account any specified number of spherical0075harmonic terms in the DGRF/IGRF field model, updated for a given epoch.0073It also contains subsidiary codes for transformations between several0072coordinate systems most widely used in geophysics and space physics.00040076     The present version of GEOPACK does not include subroutines for an 0077external field model.  This area is now rapidly developing: there already0078exist several models, each of them has its own advantages and deficiencies0078and it is thus not easy to choose between them. Also, we expect better mo-0077dels to appear in the nearest future. Therefore, it was considered reaso-0077nable to divide the software into two files;  the first one (GEOPACK) in-0077cludes only the most "stable" general-use codes, which are unlikely to be0077significantly changed in future, while the external field models are pro-0040vided separately in the second file.^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^0077     A convenient way of using the GEOPACK subroutines is to compile them0077separately from the main user's code and include the corresponding object0043modules in the user's personal library.00040074     The set of subroutines was developed originally in 1978 and since0076then had been used by researchers at the Leningrad University and other 0076geophysical institutions in the USSR and abroad. The present version of 0075the GEOPACK is a result of amending and upgrading its earlier versions,0075based on the author's own experience and many critical comments of his 0078colleagues received during the last 15 years, and it proved to be a robust0077tool.   However, it does not mean that all the possibilities for further 0077improvement are already exhausted. The author would thus greatly appreci-0077ate any comments on the performance of the codes, possible problems, and 0077any suggestions on how to make the GEOPACK subroutines clearer, simpler, 0031faster, and more versatile.00040077      Below is a kind of "manual guide" for using the subroutines, inclu-0077ding the list of references.    The  FORTRAN texts are placed in the file0076GEOPACK.FOR, which also contains two examples of a typical main program 0078for tracing the model field lines, with the purpose to help users to debug0047their own codes and avoid typical mistakes.00040004000400040051                DESCRIPTIONS OF THE SUBROUTINES0004000400040030     1. SUBROUTINE:   IGRF00040073        FUNCTION:     Computes three spherical components of the main0077(internal) geomagnetic field in the geographical coordinate system (GEO).00040064        FORTRAN STATEMENT:   CALL IGRF(IY,NM,R,T,F,BR,BT,BF)00040073        INPUT PARAMETERS:    IY is year (four digits) which should be0073specified in the interval 1965 <= IY <= 1990;  NM  is the maximal or-0070der of spherical harmonics taken into account in the corresponding0074expansions (1<=NM<=10);  R,T,F  are spherical coordinates of the point^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^0073(radial geocentric distance in Earth radii Re=6371.2 km,  colatitude,0049and east longitude in radians, respectively).00040074        OUTPUT PARAMETERS:  BR, BT, BF are the spherical components of0071the main geomagnetic field in nanotesla: BR is the radial component0075(BR>0 corresponds to outward direction), BT is the southward component,0037and BP is the eastward component.00040033        COMMON BLOCKS:  none.00040050        REFERENCES TO OTHER SUBROUTINES: none.00040074        COMMENTS:   The subroutine uses harmonic coefficients for five0076epochs: 1965.0, 1970.0, 1975.0, 1980.0, and 1985.0 ; the appropriate va-0076lues of coefficients are automatically calculated by means of linear in-0075terpolation between the nearest epochs. If 1985<IY<1990, the subroutine0074extrapolates by using the secular coefficients. If IY<1965 or IY>1990,0076the coefficients are assumed equal to those for IY=1965 or IY=1990, res-0076pectively.   The subroutine can be easily upgraded, as soon as the coef-0060ficients for the next epochs become available in future.0008    0004000400040028     2. SUBROUTINE:  DIP00040073        FUNCTION:  Computes Cartesian solar-magnetospheric (GSM) com-0075ponents of the Earth's magnetic field corresponding to the first (dipo-0050lar) term in the spherical harmonic expansion.00040059        FORTRAN STATEMENT:  CALL DIP(PS,X,Y,Z,BX,BY,BZ)00040075        INPUT PARAMETERS:  PS is the angle in radians between the dipo-0073le moment axis and the ZGSM axis (PS<0 in winter and PS>0 in summer);0071X,Y,Z  are Cartesian GSM coordinates of the point in Earth's radii.00040073        OUTPUT PARAMETERS: BX,BY,BZ are the GSM magnetic field compo-0023nents in nanotesla.00040032        COMMON BLOCKS: none.00040051        REFERENCES TO OTHER SUBROUTINES:  none.00040069        COMMENTS: (1) The dipole is assumed to be centered at the0069origin; the value of its magnetic moment corresponds to the epoch00311980.0 (N.W. Peddie, 1982).0072                  (2) The angle PS can either be specified directly,^^^^^^^^^^^0070or calculated by using the subroutine RECOMP; in the last case the0071unlabelled COMMON block should be present in the user's module from0071which RECOMP and DIP are called. The value of the dipole tilt angle0069(in radians) will be placed by RECOMP in the 16-th element of the0070COMMON-block. The structure of the COMMON-block should be the same0072as in the subroutine RECOMP (see the examples in the end of the file0017GEOPACK.FOR).00040004000400040028      3. SUBROUTINE: SUN00040071         FUNCTION: This is a subsidiary subroutine which is usually0071called from the subroutine RECOMP and computes the spherical angles0072for the Earth-Sun line in the geocentric inertial coordinate system,0041and the Greenwich mean sidereal time.00040073      FORTRAN STATEMENT:  CALL SUN(IYR,IDAY,IHOUR,MIN,ISEC,GST,SLONG,0047                                SRASN,SDEC)00040073      INPUT PARAMETERS:  IYR,IDAY,IHOUR,MIN,ISEC  are the year (four-0067digit number), the day, hour, minute, and second, respectively;00040077      OUTPUT PARAMETERS:  GST is the Greenwich mean sidereal time, SLONG,0077SRASN,SDEC are the ecliptical longitude, right ascension, and declination0072of the Sun (all the output parameters are in radians), respectively.00040034         COMMON BLOCKS:  none.00040051         REFERENCES TO OTHER SUBROUTINES: none.00040041         COMMENTS:  (1) 1901<IYR<20990059                    (2) January 1 corresponds to IDAY=10076                    (3) The subroutine was compiled (with minor changes)0067from the paper by C.T. Russell (1971), its author is G.D. Mead.00040004000400040032      4. SUBROUTINE:  SPHCAR00040074         FUNCTION:  Computes spherical coordinates from Cartesian ones0018or vica versa.00040064         FORTRAN STATEMENT:  CALL SPHCAR(R,TETA,PHI,X,Y,Z,J)00040075         INPUT PARAMETERS:  J is the integer key parameter: if J>0 then0073spherical coordinates R,TETA,PHI are the input quantities (colatitude0077TETA and longitude PHI are in radians); if J<0 then Cartesian coordinates0025X,Y,Z are the inputs.^00040073         OUTPUT PARAMETERS:  If J>0 then Cartesian coordinates X,Y,Z;0027if J<0 then R,TETA,PHI.00040034         COMMON BLOCKS:  none.00040051         REFERENCES TO OTHER SUBROUTINES: none.00040075         COMMENTS:  If, for J<0, X=0.0 and Y=0.0, then PHI is set equal0009to 0.00040004000400040032      5. SUBROUTINE:  BSPCAR00040076         FUNCTION:  Calculation of Cartesian components of a vector from0067the spherical ones and the spherical coordinates theta and phi.00040076         FORTRAN STATEMENT:  CALL BSPCAR(TETA,PHI,BR,BTET,BPHI,BX,BY,BZ)00040077         INPUT PARAMETERS:  TETA and PHI are the colatitude and longitude0075of the point in radians; BR,BTET, and BPHI are the vector components in0032the local coordinate system.00040077         OUTPUT PARAMETERS: BX,BY,BZ are Cartesian components of the vec-0008tor.0033         COMMON BLOCKS: none.00040051         REFERENCES TO OTHER SUBROUTINES: none.000400040004000400040031      6. SUBROUTINE: RECALC00040077         FUNCTION:  Computes the angles defining the geodipole orientati-0078on for a given UT moment as well as elements of matrices for transformati-0070ons between the following Cartesian geocentric coordinate systems:0075geographic (GEO), dipole geomagnetic (MAG), solar-magnetic (SM), solar-0053magnetospheric (GSM), and solar-ecliptical (GSE).00040069         FORTRAN STATEMENT:  CALL RECALC(IYR,IDAY,IHOUR,MIN,ISEC)00040077         INPUT PARAMETERS:  IYR,IDAY,IHOUR,MIN, and ISEC  are the same as0032those in the subroutine SUN.00040037         OUTPUT PARAMETERS: none.00040065         REFERENCES TO OTHER SUBROUTINES:  SUN, MAGSM, SMGSM.00040075         COMMON BLOCKS:  the results are placed in the unlabeled common0030block containing 37 words.00040075         COMMENTS:  If it is necessary to make calculations only in the0072coordinate systems GEO and MAG (no UT dependence) then the parameter0077IHOUR can be assigned any value larger than 24; in this case other matri-0028ces will not be defined.^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^0075      The subroutine can be easily upgraded for 1990 and later epochs, 0078as soon as the corresponding internal field coefficients become available.0004000400040033      7. SUBROUTINE:  GEOMAG.00040076         FUNCTION:  Transformation of the geographic Cartesian coordina-0057tes into the dipolar geomagnetic ones and vica versa.00040075         FORTRAN STATEMENT:  CALL GEOMAG(XGEO,YGEO,ZGEO,XMAG,YMAG,ZMAG,0042                                J,IYR)00040075         INPUT PARAMETERS:  J is the integer key parameter: if J>0 then0071the geographic Cartesian coordinates XGEO,YGEO,ZGEO are the inputs;0073if J<0 then the geomagnetic dipole coordinates XMAG,YMAG,ZMAG are the0074inputs; IYR is year (the subroutine takes into account secular changes0034of the Earth's dipole moment).00040073          OUTPUT PARAMETERS:  If J>0 then XMAG,YMAG,ZMAG; if J<0 then0050                               XGEO,YGEO,ZGEO.00040052          COMMON BLOCKS:  The same as in RECOMP.00040055          REFERENCES TO OTHER SUBROUTINES:  RECOMP.0014          0075          COMMENTS: The subroutine can be easily upgraded for 1990 and 0075later epochs, as soon as the corresponding internal field coefficients 0021become available.0004000400040004000400040033       8. SUBROUTINE:  MAGSM.00040073          FUNCTION:  Transformation of the Cartesian dipolar magnetic0059coordinates into the solar magnetic ones or vica versa.00040073         FORTRAN STATEMENT:  CALL MAGSM(XMAG,YMAG,ZMAG,XSM,YSM,ZSM,J)00040075         INPUT PARAMETERS:  J is the integer key parameter: if J>0 then0067the dipolar magnetic coordinates XMAG,YMAG,ZMAG are the inputs;0066if J<0 then the solar magnetic coordinates XSM,YSM,ZSM are the0011inputs.00040070          OUTPUT PARAMETERS:  If J>0 then XSM,YSM,ZSM; if J<0 then0050                               XMAG,YMAG,ZMAG.00040052          COMMON BLOCKS:  The same as in RECOMP.00040053          REFERENCES TO OTHER SUBROUTINES:  none.000400040004000400040033      9. SUBROUTINE:  GSMGSE.0004^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^0076          FUNCTION:  Transformation of the geocentric solar magnetosphe-0077ric coordinates into the geocentric solar ecliptical ones and vica versa.00040077          FORTRAN STATEMENT: CALL GSMGSE(XGSM,YGSM,ZGSM,XGSE,YGSE,ZGSE,J)00040076          INPUT PARAMETERS:  J is the integer key parameter: if J>0 then0071the solar magnetospheric coordinates XGSM,YGSM,ZGSM are the inputs;0075if J<0 then the solar ecliptical coordinates XGSE,YGSE,ZGSE are the in-0009puts.00040073          OUTPUT PARAMETERS:  If J>0 then XGSE,YGSE,ZGSE; if J<0 then0050                               XGSM,YGSM,ZGSM.00040052          COMMON BLOCKS:  The same as in RECOMP.00040053          REFERENCES TO OTHER SUBROUTINES:  none.000400040004000400040033      10. SUBROUTINE:  SMGSM.00040073          FUNCTION:  Transformation of the solar magnetic coordinates0055 into the solar magnetospheric ones and vica versa.00040074          FORTRAN STATEMENT:  CALL SMGSM(XSM,YSM,ZSM,XGSM,YGSM,ZGSM,J)00040071          INPUT PARAMETERS:  J is the integer key parameter: if J>00067then the solar magnetic coordinates XSM,YSM,ZSM are the inputs;0071if J<0 then the solar magnetospheric coordinates XGSM,YGSM,ZGSM are0015the inputs.00040074           OUTPUT PARAMETERS:  If J>0 then XGSM,YGSM,ZGSM; if J<0 then0047                               XSM,YSM,ZSM.00040053           COMMON BLOCKS:  The same as in RECOMP.00040054           REFERENCES TO OTHER SUBROUTINES:  none.000400040004000400040035       11. SUBROUTINE:  GEOGSM.00040078           FUNCTION:  Transformation of the geographic Cartesian coordina-0058tes into the solar magnetospheric ones and vica versa.00040079           FORTRAN STATEMENT:  CALL GEOGSM(XGEO,YGEO,ZGEO,XGSM,YGSM,ZGSM,J)00040077           INPUT PARAMETERS:  J is the integer key parameter: if J>0 then0071the geographic Cartesian coordinates XGEO,YGEO,ZGEO are the inputs;0075if J<0 then the solar magnetospheric coordinates XGSM,YGSM,ZGSM are the0011inputs.00040074           OUTPUT PARAMETERS:  If J>0 then XGSM,YGSM,ZGSM; if J<0 then^^^^^^^^^^^^^^0050                               XGEO,YGEO,ZGEO.00040053           COMMON BLOCKS:  The same as in RECOMP.00040054           REFERENCES TO OTHER SUBROUTINES:  none.000400040004000400040034       12. SUBROUTINE:  RHAND.00040075           FUNCTION: This is a subsidiary subroutine which computes the0076right-hand sides of the field line equations (that is, the components of0074the unit vector defining the local direction of B) which are necessary0073for making one step along the field line made by the subroutine STEP.00040073           FORTRAN STATEMENT:  CALL RHAND(X,Y,Z,R1,R2,R3,IOPT,EXNAME)00040073           INPUT PARAMETERS:  X,Y,Z  are GSM Cartesian coordinates of0075the current point of space, IOPT is the number of the option of the ex-0072ternal magnetic field model; EXNAME is the name of an external field0072model subroutine  (see below the description of the subroutine TRACE0022for more details).00040074            OUTPUT PARAMETERS:  R1,R2,R3 are the right-hand-side quan-0040tities (the unit vector components).00040053            COMMON BLOCKS: the same as in RECOMP.00040075            REFERENCES TO OTHER SUBROUTINES: DIP, EXNAME, IGRF, GEOGSM,0063                                            SPHCAR, BSPCAR.000400040004000400040033        13. SUBROUTINE: STEP.00040075           FUNCTION:  Makes one step along the force line of the magne-0073tic field corresponding to the sum of contributions from internal and0074external sources computed by the subroutines IGRF (or DIP) and EXNAME,0017respectively.00040074           FORTRAN STATEMENT:  CALL STEP(N,X,Y,Z,DS,ERRIN,IOPT,EXNAME)00040076           INPUT PARAMETERS:  N is the maximal order of spherical harmo-0072nics in the main field expansions; X,Y,Z  are GSM coordinates of the0075initial point; DS is the step size; ERRIN is the estimate of the admis-0074sible error; IOPT is the number of the option of the external magnetic0073field model; EXNAME is the name of an external field model subroutine^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^0073(see below the description of the subroutine TRACE for more details).00040072           OUTPUT PARAMETERS:  X,Y,Z are GSM coordinates of the next0059point on the traced field line (after making one step).00040052           COMMON BLOCKS: the same as in RECOMP.00040055           REFERENCES TO OTHER SUBROUTINES:  RHAND.00040076           COMMENTS:  (1) If the subroutine is used separately (that is,0073called not via the subroutine TRACE), then it is necessary to specify0069the actual name of the subroutine EXNAME in the FORTRAN statement0075EXTERNAL in the module , from which STEP is to be called (see also com-0040ments 1-3 for the subroutine TRACE).0072                     (2) The parameter DS can change its value after0073having called STEP, since in case of too large actual error (that is,0074exceeding ERRIN) the step is repeated once again with the halved value0075of DS; if the actual error appears to be much smaller (less than 4 per-0072cent) than ERRIN, then the step is not repeated, but the value of DS0033for the next step is doubled.000400040004000400040034       14. SUBROUTINE:  TRACE.00040075           FUNCTION:  Computation of GSM coordinates of points lying on0075the geomagnetic field line, beginning from a given initial point, up to0073its intersection  with the inner or the outer spherical boundaries of0023the tracing region.00040076           FORTRAN STATEMENT:  CALL TRACE(XI,YI,ZI,DIR,RLIM,R0,IHARM,NP,0070                                  IOPT,EXNAME,XF,YF,ZF,XX,YY,ZZ,L)00040074           INPUT PARAMETERS:  XI,YI,ZI, are the GSM coordinates of the0072initial point (in Earth's radii); DIR defines the tracing direction:0072for DIR=1.0 it is opposite to that of the B vector, and for DIR=-1.00074it coincides with B direction. Parameter RLIM is the radius of spheri-0072cal outer boundary of the tracing region (in Earth's radii), so that0075the tracing is ceased once R>RLIM;  R0 is the radius of the inner sphe-0075rical boundary of the tracing region (usually R0=1.0, which corresponds^^^^^^^^^^^^^^^^^0074to the Earth's surface; see also the comment No.4); IHARM is the maxi-0072mal order of spherical harmonics taken into account in computing the0074main part of the B-field (see the comment No.2); NP is the upper esti-0071mate for the expected number of steps along the field line segment;0073IOPT is the number of the option of the external magnetic field model0072(see the description of the corresponding subroutine); EXNAME is any0074name of a subroutine that returns the components of the magnetic field0075contributed by the external sources. The corresponding actual parameter0075should have the type CHARACTER*n (n is the total number of symbols) and0071must be specified in the EXTERNAL statement of the main module. The0074list of formal parameters of the subroutine EXNAME must be as follows:0078 (IOPT, PS, X, Y, Z, BX, BY, BZ), where IOPT is an integer number reserved0075for specifying a concrete version of the external field model (e.g. the0078KP-index), the other 6 parameters are identical to those in the subroutine0008DIP.0069   See an example of using this subroutine in the end of the file0033GEOPACK.FOR for more details.00040053           COMMON BLOCKS:  the same as in RECALC.00040060           REFERENCES TO OTHER SUBROUTINES: STEP, RHAND.00040071           COMMENTS:  (1) Before making field line computations for0073concrete geophysical situations it is necessary to specify some quan-0073tities defining relative orientation of geographical and solar-magne-0074tospheric coordinate systems. This should be done by having called the0075subroutine  RECOMP with appropriate values of parameters defining year,0074day, and universal time of the day. In this case, it is recommended to0074take IHARM between 7 and 10, bearing in mind that the actual length of0076the corresponding expansions will be set automatically, depending on the0076current value of the geocentric distance, significantly saving thus com-0015puter time.0074                      (2) If the dipolar approximation for the main B-^^^^^^^^^^^^^^^^^^^^^^0073field is sufficient, then it is possible to set IHARM=0; in this case0073the subroutine DIP will be used instead of IGRF, and the computations0074will run much more rapidly. The corresponding value of the dipole tilt0075angle should be specified before tracing, by putting it (in radians) in0073the 16-th element of the unlabelled COMMON-block. Another possibility0070is to set IHARM=1 and call RECOMP before tracing; in this case the0076tilt angle and other quantities that are necessary for coordinate trans-0077formations will be automatically found from the year, day, and UT values.00040070                     (3) Any other external field model subroutine0075can be used, provided it has the same list of formal parameters and its0072name is included in the EXTERNAL statement in the user's module that0016calls TRACE.00040073                     (4) If the field line returns to the Earth, then0074the last (L-th) elements in the arrays XX,YY,ZZ as well as the triplet0071XF,YF,ZF  correspond to the intersection point of the line with the0024sphere of radius R0.0070                     (5) If L<NP, then the rest (NP-L) elements of0070XX,YY,ZZ are filled with the values of the final point coordinates0067(XF,YF,ZF). If L>NP (that is, not enough space reserved for the0070computation results) then the computation is aborted and a warning0071diagnostics is displayed. In most cases, it is enough to take NP of0021the order of 300.0004000400040075******************** Description of new subroutines *******************0075********************      added by M. Peredo        *******************0075********************             and                *******************0075******************** description of external field  *******************0075********************   subroutines old and new      *******************0004000400040033       14. SUBROUTINE:  BTOT.00040073           FUNCTION: This subroutine can be used to compute the total0070(ie. internal plus external) magnetic field at a series of points.0004^^^^^^^^^^^^^^^^^^^^^^0078           FORTRAN STATEMENT:  CALL BTOT(X,Y,Z,EXNAME,IOPT,BX,BY,BZ,count)00040076           INPUT PARAMETERS:  X,Y,Z  are arrays of size count containing0074the GSM Cartesian coordinates of the points at which the field must be0079determined, IOPT is the number of the option of the external magnetic field0078model; EXNAME is the name of an external field model subroutine, and count0050is the dimension of the input position arrays.00040080            OUTPUT PARAMETERS:  BX,BY,BZ are arrays of size count containing0080the GSM components of the magnetic field (in nT) at the points corresponding0047to the input positions in the arrays X,Y,Z.00040053            COMMON BLOCKS: the same as in RECALC.00040070            REFERENCES TO OTHER SUBROUTINES: EXNAME, IGRF, GEOGSM,0063                                            SPHCAR, BSPCAR.0004000400040035       15. SUBROUTINE:  EXT87W.00040069           FUNCTION: This subroutine provides an interface to the0075subroutine FUNCTS which computes the external magnetic field components0076according to the warped T87 model, i.e. T87W; depending on the value of 0077IOPT, versions of the model derived from data binned by Kp or by Psw and 0074IMF-Bz are selected. A dummy interface subroutine is used to match the0070calling sequence in the original external field subroutines in the0074geopack set; and thus allow consistent calling for any of the external0017field models.00040071           FORTRAN STATEMENT:  CALL EXT87W(IOPT,PSI,X,Y,Z,BX,BY,BZ)00040077           INPUT PARAMETERS:  IOPT is a parameter to select which version0076of the T87W model is used, X,Y,Z are the GSM coordinates of the point at0077which the field is computed, and PSI is the dipole tilt angle in radians.00040075            OUTPUT PARAMETERS:  BX,BY,BZ are the GSM components of the 0075magnetic field (in nT) at the point corresponding to the input position0013in X,Y,Z.00040044            COMMON BLOCKS: WARP, MODPAR.00040056            REFERENCES TO OTHER SUBROUTINES: FUNCTS.00040004^^^^^^^^^^^^^^^^^^^^^^^^^^^^^0035       16. SUBROUTINE:  FUNCTS.00040077           FUNCTION: This subroutine computes the external magnetic field0078components according to the warped T87 model, i.e. T87W; depending on the 0078value of IOPT, versions of the model derived from data binned by Kp or by 0032Psw and IMF-Bz are selected.00040074           FORTRAN STATEMENT:  CALL FUNCTS(X,Y,Z,TILT,A,BMODL,DBDA,MA)00040079           INPUT PARAMETERS:  X,Y,Z are the GSM coordinates of the point at0078which the field is computed, and TILT is the dipole tilt angle in degrees.0073A contains the model parameters while MA is the number of parameters.00040078            OUTPUT PARAMETERS:  BMODL a 3-component vector containing the 0080cartesian components of the magnetic field, DBDA a 50 by 3 matrix containing0081the derivatives with respect to the fitting coefficients of the 3 components 0032of the model magnetic field.00040038            COMMON BLOCKS: MODPAR.00040074            REFERENCES TO OTHER SUBROUTINES: NLTAIL, NLRING, NLMPAUSE.000400040054       17. SUBROUTINES:  NLTAIL, NLRING, NLMPAUSE.00040070           FUNCTION: These subroutines compute, respectively, the 0068contributions from the tail current system, the ring current and0069the combined magnetopause and Birkeland current systems (i.e. the0040polynomial terms) in the T87W model.00040074           FORTRAN STATEMENT:  CALL NLTAIL(X,Y,Z,TILT,MA,BTAIL,DBTAIL)0053			       CALL NLRING(X,Y,Z,TILT,MA,BRING,DBRING)0059			       CALL NLMPAUSE(X,Y,Z,TILT,MA,BMPAUSE,DBMPAUSE)00040079           INPUT PARAMETERS:  X,Y,Z are the GSM coordinates of the point at0078which the field is computed, and TILT is the dipole tilt angle in degrees.0073A contains the model parameters while MA is the number of parameters.00040076            OUTPUT PARAMETERS:  BTAIL, BRING or BMPAUSE are 3-component 0078vectors containing the cartesian components of the magnetic field, DBTAIL,0076DBRING or DBMPAUSE are 50 by 3 matrices containing the derivatives with ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^0073respect to the fitting coefficients of the 3 components of the model 0019magnetic field.00040062            COMMON BLOCKS: MODPAR, WARP (only for NLTAIL).00040074            REFERENCES TO OTHER SUBROUTINES: NLTAIL, NLRING, NLMPAUSE.0004000400040035       18. SUBROUTINE:  EXT87L.00040072           FUNCTION: This subroutine computes the external magnetic 0077field components according to the original T87 (long) model the value of 0070IOPT identifies the version of the model (by Kp value) to be used.00040071           FORTRAN STATEMENT:  CALL EXT87L(IOPT,PSI,X,Y,Z,BX,BY,BZ)00040077           INPUT PARAMETERS:  IOPT is a parameter to select which version0075of the T87 model is used, X,Y,Z are the GSM coordinates of the point at0077which the field is computed, and PSI is the dipole tilt angle in radians.00040075            OUTPUT PARAMETERS:  BX,BY,BZ are the GSM components of the 0075magnetic field (in nT) at the point corresponding to the input position0013in X,Y,Z.00040031            COMMON BLOCKS: 00040049            REFERENCES TO OTHER SUBROUTINES: 000400040035       19. SUBROUTINE:  EXT87T.00040072           FUNCTION: This subroutine computes the external magnetic 0078field components according to the original T87 (short) model the value of 0070IOPT identifies the version of the model (by Kp value) to be used.00040071           FORTRAN STATEMENT:  CALL EXT87T(IOPT,PSI,X,Y,Z,BX,BY,BZ)00040077           INPUT PARAMETERS:  IOPT is a parameter to select which version0075of the T87 model is used, X,Y,Z are the GSM coordinates of the point at0077which the field is computed, and PSI is the dipole tilt angle in radians.00040075            OUTPUT PARAMETERS:  BX,BY,BZ are the GSM components of the 0075magnetic field (in nT) at the point corresponding to the input position0013in X,Y,Z.00040031            COMMON BLOCKS: 00040049            REFERENCES TO OTHER SUBROUTINES: 000400040036       20. SUBROUTINE:  EXT89KP.00040072           FUNCTION: This subroutine computes the external magnetic ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^0070field components according to the original T89 model the value of 0070IOPT identifies the version of the model (by Kp value) to be used.00040072           FORTRAN STATEMENT:  CALL EXT89KP(IOPT,PSI,X,Y,Z,BX,BY,BZ)00040077           INPUT PARAMETERS:  IOPT is a parameter to select which version0075of the T87 model is used, X,Y,Z are the GSM coordinates of the point at0077which the field is computed, and PSI is the dipole tilt angle in radians.00040075            OUTPUT PARAMETERS:  BX,BY,BZ are the GSM components of the 0075magnetic field (in nT) at the point corresponding to the input position0013in X,Y,Z.00040031            COMMON BLOCKS: 00040049            REFERENCES TO OTHER SUBROUTINES: 0004000400040036       21. SUBROUTINE:  EXT89AE.00040072           FUNCTION: This subroutine computes the external magnetic 0061field components according to the T89 model the value of 0070IOPT identifies the version of the model (by AE value) to be used.00040072           FORTRAN STATEMENT:  CALL EXT89AE(IOPT,PSI,X,Y,Z,BX,BY,BZ)00040077           INPUT PARAMETERS:  IOPT is a parameter to select which version0075of the T87 model is used, X,Y,Z are the GSM coordinates of the point at0077which the field is computed, and PSI is the dipole tilt angle in radians.00040075            OUTPUT PARAMETERS:  BX,BY,BZ are the GSM components of the 0075magnetic field (in nT) at the point corresponding to the input position0013in X,Y,Z.00040031            COMMON BLOCKS: 00040049            REFERENCES TO OTHER SUBROUTINES: 000400040035       22. SUBROUTINE:  EXT89M.00040074           FUNCTION: This subroutine provides a dummy interface to the0080routine T89M which computes the external magnetic field components according0080to the new version of the T89 model based on an expanded data set containing0075ISEE observations in addition to the IMP-HEOS observations on which the0081original T89 was based; the value of IOPT identifies the version of the model0029(by Kp value) to be used.0004^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^0071           FORTRAN STATEMENT:  CALL EXT89M(IOPT,PSI,X,Y,Z,BX,BY,BZ)00040077           INPUT PARAMETERS:  IOPT is a parameter to select which version0075of the T87 model is used, X,Y,Z are the GSM coordinates of the point at0077which the field is computed, and PSI is the dipole tilt angle in radians.00040075            OUTPUT PARAMETERS:  BX,BY,BZ are the GSM components of the 0075magnetic field (in nT) at the point corresponding to the input position0013in X,Y,Z.00040034            COMMON BLOCKS: PAR00040053            REFERENCES TO OTHER SUBROUTINES: T89M0004000400040033       23. SUBROUTINE:  T89M.00040077           FUNCTION: This subroutine computes the external magnetic field0073components according to the new version of the T89 model based on an 0078expanded data set containing ISEE observations in addition to the IMP-HEOS0053observations on which the original T89 was based.00040064           FORTRAN STATEMENT:  CALL T89M (ID, A, XI, F, DER)00040078           INPUT PARAMETERS:  ID is the number of the data point in a set 0075(initial assignments are performed only for ID=1, saving thus CPU time)0073A is an input vector containing model parameters; XI is input vector 0036containing independent variables00040077            OUTPUT PARAMETERS:  F is a double precision vector containing0076the calculated values of dependent variables; DER is a double precision 0080vector containing the calculated values for the derivatives of the dependent0050variables with respect to the model parameters000400040030            COMMON BLOCKS:00040048            REFERENCES TO OTHER SUBROUTINES:0004000400040004000400040036                      REFERENCES00040074    Golovkov, V.P. and G.I. Kolomiytseva. The International Analytical0070Field and its secular trend for the 1980-1990 period. Geomagnetism0057and Aeronomy (Engl.Transl.), v.26, No.3, p.439, 1986.00040072    Peddie, N.W. A third generation of International Geomagnetic Re-0074ference Field  models. J. Geomagn.Geoelectr., v.34, No.6, p.310, 1982.0004^^^^^^^^^^^^^^^^^^^^^^0076Peredo, M., D.P. Stern and N.A. Tsyganenko, "Are Existing Magnetospheric0057Models Excessively Stretched?", JGR, 1993 (in press).00040074    Russell, C.T. Geophysical coordinate transformations. Cosmic Elec-0034trodynamics, v.2, p.184, 1971.00040071    Tsyganenko, N.A., A.V. Usmanov, V.O. Papitashvili, N.E. Papita-0073shvili, and V.A. Popov.  Software for computations of the geomagnetic0072field and related coordinate systems, Soviet Geophys.Committee, Mos-0014cow, 1987.00040070    Tsyganenko, N.A.  A magnetospheric magnetic field model with a0065warped tail plasma sheet. Planet.Space Sci., v.37, p.5, 1989.00040077    Tsyganenko, N.A.  Quantitative models of the magnetospheric magnetic 0071field: methods and results. Space Sci. Rev., v.54, pp.75-186, 1990.0005c0005c0005c000400040005C0077C------------------------------------------------------------------------0005C0075C######################################################################0075C     THE CODES BELOW GIVE TWO EXAMPLES OF TRACING THE FIELD LINES    #0075C                   BY USING THE GEOPACK SOFTWARE                     #0075C    (just for illustration; do not include them in your library)     #0075C######################################################################0005C0026      PROGRAM EXAMPLE10005C0074C    IN THIS EXAMPLE IT IS ASSUMED THAT WE KNOW GEOGRAPHIC COORDINATES0079C     OF A POINT AT THE EARTH'S SURFACE AND WANT TO TRACE THE CORRESPONDING0081C     FIELD LINE AT A SPECIFIED MOMENT OF UNIVERSAL TIME, TAKING INTO ACCOUNT0070C     THE HIGHER-ORDER MULTIPOLES IN THE INTERNAL FIELD EXPANSION.0005C0043      DIMENSION XX(500),YY(500),ZZ(500)0005c0080c  do not forget to include the following EXTERNAL statement in your codes !0005c0025      EXTERNAL EXTERN0005C0063c	Replace EXTERN in the last statement with the actual name0065c	of the external field routine you wish to use (e.g. EXT87L,0050c 	EXT87T, EXT87W, EXT89KP, EXT89AE or EXT89M)0005c0072C   FIRST OF ALL WE DEFINE THE TIME MOMENT AND COMPUTE THE NECESSARY^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^0073C   QUANTITIES FOR COORDINATE TRANSFORMATIONS BY USING THE SUBROUTINE0068C   RECOMP: IN THIS PARTICULAR CASE WE NEED TO TRACE FIELD LINES0045C   FOR YEAR=1967, IDAY=350, UT HOUR = 210005C0038      CALL RECALC(1967,350,21,0,0)0005C0071C  IOPT=1 CORRESPOND TO THE MOST QUIET EXTERNAL FIELD MODEL VERSION0005C0016      IOPT=10005C0067C    THE LINE WILL BE TRACED FROM THE POINT HAVING GEOGRAPHICAL0054C    LONGITUDE 45 DEGREES AND LATITUDE 75 DEGREES:0005C0020      GEOLAT=75.0020      GEOLON=45.0038      COLAT=(90.-GEOLAT)*.017453290031      XLON=GEOLON*.017453290005C0049C   CONVERT SPHERICAL COORDS INTO CARTESIAN :0005C0053      CALL SPHCAR(1.,COLAT,XLON,XGEO,YGEO,ZGEO,1)0005C0065C   TRANSFORM GEOGRAPHICAL COORDS INTO SOLAR MAGNETOSPHERIC :0005C0054      CALL GEOGSM(XGEO,YGEO,ZGEO,XGSM,YGSM,ZGSM,1)0005C0076C   TRACE THE FIELD LINE (EXTERN IS THE NAME OF THE EXTERNAL FIELD MODEL0078C     SUBROUTINE, THAT SHOULD BE COMPILED AND LINKED WITH THE GEOPACK SUB-0020C     ROUTINES):0005C0052      CALL TRACE(XGSM,YGSM,ZGSM,1.,60.,1.,7,500,0044     *  IOPT,EXTERN,XF,YF,ZF,XX,YY,ZZ,M)0005C0057C   WRITE THE RESULTS IN THE DATAFILE 'LINTEST1.DAT':0005C0057        OPEN(UNIT=1,FILE='LINTEST1.DAT',STATUS='NEW')0050        WRITE (1,20) (XX(L),YY(L),ZZ(L),L=1,M)0030 20     FORMAT((2X,3F6.2))0025        CLOSE(UNIT=1)0014      STOP0013      END0072C-------------------------------------------------------------------0005C0026      PROGRAM EXAMPLE20005C0054C   THIS IS AN EXAMPLE OF USING 'TRACE' SUBROUTINE0065C    IN CASE OF EXPLICIT SPECIFICATION OF DIPOLE TILT ANGLE .0073C    IT IS ASSUMED THAT WE KNOW THE GSM LATITUDE AND LONGITUDE OF THE0046C    STARTING POINT AT THE EARTH'S SURFACE0082C      ONE MORE ASSUMPTION IS THAT THE MAIN (INTERNAL) FIELD IS PURELY DIPOLAR0005C0059      COMMON AA(10),SPS,CPS,BB(3),PS,CC(11),KK(2),DD(8)0043      DIMENSION XX(500),YY(500),ZZ(500)0005c0080c  do not forget to include the following EXTERNAL statement in your codes !0005c0025      EXTERNAL EXTERN^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^0069c       Replace EXTERN in the last statement with the actual name0071c       of the external field routine you wish to use (e.g. EXT87L,0055c       EXT87T, EXT87W, EXT89KP, EXT89AE or EXT89M)0005c0005C0016      IOPT=10005c0018      XLAT=75.0019      XLON=180.0005C0081c  we have defined the latitude (XLAT) and longitude (XLON) in the GSM system0005c0016       PS=0.0022       SPS=SIN(PS)0022       CPS=COS(PS)0005C0083C WE HAVE SPECIFIED THE DIPOLE TILT ANGLE (PS), ITS SINE (SPS) AND COSINE (CPS)0005C0033       T=(90.-XLAT)*.017453290028       XL=XLON*.017453290029      XGSM=SIN(T)*COS(XL)0029      YGSM=SIN(T)*SIN(XL)0021      ZGSM=COS(T)0052      CALL TRACE(XGSM,YGSM,ZGSM,1.,60.,1.,0,500,0044     *  IOPT,EXTERN,XF,YF,ZF,XX,YY,ZZ,M)0005C0066C    THE MEANING OF INPUT PARAMETERS FOR TRACE IS EXPLAINED IN0031C     THE FILE  GEOPACK.TXT0005C0057C   WRITE THE RESULTS IN THE DATAFILE 'LINTEST2.DAT':0005C0057       OPEN(UNIT=1, FILE='LINTEST2.DAT',STATUS='NEW')0048  1   WRITE (1,20) (XX(L),YY(L),ZZ(L),L=1,M)0028 20   FORMAT((2X,3F6.2))0024       CLOSE(UNIT=1)0014      STOP0013      END0004^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^