**[Abstract:]{.underline}** The new version (02/04/2008) of the GEOPACK
library includes 20 FORTRAN subroutines, to be used in various studies
that involve calculations of the geomagnetic field in the Earth's
magnetosphere using empirical models and/or spacecraft observations. The
library includes subroutines for the current (IGRF) and past (DGRF)
internal geomagnetic field models, a group of subroutines for
transformations between various coordinate systems, a field line tracer,
and two magnetopause model codes.

**PREFATORY NOTES TO THIS RELEASE (FEBRUARY 08, 2008)**

To avoid inadvertent use of obsolete subroutines from earlier versions,
a suffix 08 was added to the end of each subroutine's name. With respect
to their content, three essential changes/additions were made in this
version, as follows below.

The first modification is the possibility to calculate components of a
vector in the "Geocentric Solar-Wind" (GSW) coordinate system \[*Hones
et al.*, 1986\], also known as GSWM \[*Tsyganenko et al.*, 1998;
*Tsyganenko and Fairfield*, 2004\]. That system is analogous to the
standard geocentric solar-magnetospheric (GSM), except that X-axis in
the GSW coordinates is anti-parallel to the actual observed direction of
the solar wind flow, which not only aberrates by \~4^o^ from the
strictly radial line due to Earth's orbital motion around Sun, but also
often significantly fluctuates around its average direction. Orientation
of axes in the GSW system can be uniquely defined by specifying three
Cartesian components of the solar wind velocity vector (V~X~, V~Y~,
V~Z~), included as additional input parameters in the new version of the
subroutine RECALC_08. In the absence of reliable data on the solar wind
direction, one can either set V~Y~=V~Z~=0 (in which case the GSW system
becomes identical to the standard GSM) or set V~Y~=29.8 km/s, V~Z~=0,
thus taking into account the aberration effect. For more details, see
the documentation file GEOPACK-2008.DOC.

The second modification allows users to have more control over the
procedure of field line mapping using the subroutine TRACE_08. To that
end, three new input parameters were added in that subroutine, making it
possible to set (i) an upper limit, DSMAX, on the automatically adjusted
step size, (ii) a permissible step error, ERR, and (iii) maximal length,
LMAX, of arrays XX,YY,ZZ, where field line point coordinates are stored.
Minor changes in the tracing algorithm were also introduced, in order to
(i) make it more compact and easier to

understand, and (ii) prevent the algorithm from making uncontrollable
large number of multiple loops in some cases with plasmoid-like field
structures.

Finally, one more subroutine named GEODGEO_08 was added to the package,
making it possible to convert geodetic coordinates of a point in space
(altitude above the Earth\'s WGS84 ellipsoid and geodetic latitude) to
geocentric radial distance and colatitude, and vice versa.

The users should be aware of the following changes and caveats:

1\. The length of the common-block /GEOPACK1/ has been reduced in this
version from 35 to 34.

2\. Subroutine RECALC_08: (a) three GSE components of the solar wind
velocity have been added to the list of input formal arguments, (b)
common block /GEOPACK1/ has been reduced to 34 elements (see #1 above),
and (c) extensive changes were made in the body of the subroutine to
include the calculation of elements of GSW-GSE transformation matrix. As
already noted above, to convert vectors to and from the standard GSM
system, it suffices to set VGSEY=VGSEZ=0, with any non-zero value for
VGSEX.

3\. Subroutine STEP_08: (a) one more formal argument, DSMAX, was added,
that sets an upper limit on the step size.

4\. Subroutine TRACE_08: (a) three more formal arguments were added,
making it possible to explicitly specify a value of permissible step
error ERR and set any desired upper limits on the step size along the
field line, DSMAX, and on the maximal length LMAX of output arrays XX,
YY, ZZ, where coordinates of field line points are stored; (b) as in
many other subroutines, all vectors are now in the Geocentric Solar-Wind
GSW system (reduces to standard GSM by setting VGSEY=VGSEZ=0 in
RECALC_08, as detailed above).

5\. Subroutines SHUETAL_MGNP_08 and T96_MGNP_08: no changes, but the
position vectors (XGSM, YGSM, ZGSM) and (XMGNP,YMGNP,ZMGNP) can be
interpreted as referring to the GSW system (assuming that the
magnetopause orients itself following the direction of the solar wind
flow).

6\. Information on the solar wind direction (VGSEX, VGSEY, and VGSEZ)
can be obtained from available online resources of interplanetary data,
e.g., <http://cdaweb.gsfc.nasa.gov/cdaweb/sp_phys/> ,

or <http://omniweb.gsfc.nasa.gov/form/omni_min.html>. An important thing
to keep in mind is that, in most cases, the observed values of VGSEY in
those data are already reduced by 29.78 km/s, in order to take into
account the aberration effect due to Earth's orbital motion around Sun.
In order to restore the actual orientation of the observed solar wind
vector in the GSE coordinate system, the offset 29.78 km/s must be added
back to VGSEY, before using it as an input argument of the subroutine
RECALC_08. Whether or not to do that, should be determined by consulting
with the data originator, and/or by calculating average values of VGSEY
in the data set over a sufficiently long (e.g., one-year) time period.
If those averages are close to zero, then the 29.78 km/s correction must
be added back to individual values of VGSEY; otherwise (if close to \~30
km/s) they should be left intact.

**PREFATORY NOTES TO THE RELEASE OF MAY 4, 2005**

In this version, the IGRF coefficients were updated according to the
recently published table of IGRF-10 coefficients, so that the main field
model now extends through 2010 (a linear extrapolation is used for
2005 - 2010, based on the table of secular velocities). For more
details, see

<http://www.ngdc.noaa.gov/IAGA/vmod/igrf.html> (revision of 03/22/2005).

**FOREWORD TO THE RELEASE OF APRIL 22, 2003**

This collection of subroutines is a result of several upgrades of the
original package, developed almost a quarter century ago \[*Tsyganenko*,
1979\]. It represents an in-depth revision of the previous release
(January 5, 2001), with significant changes in the format of calling
statements. Users should familiarize themselves with the new formats and
rules, and accordingly adjust their source codes, as specified below.

The following changes were made to the previous release of GEOPACK:

\(1\) Subroutine IGRF, calculating the Earth\'s main field:

\(a\) Two versions of this subroutine are provided here. In the first
one (IGRF_GSM), both input (position) and output (field components) are
in the Geocentric Solar-Magnetospheric Cartesian coordinates, while the
second one (IGRF_GEO) uses spherical geographical (geocentric)
coordinates, as in the older releases.

\(b\) updating of all expansion coefficients is now made separately in
the s/r RECALC, which also takes into account the secular change of the
coefficients within a given year (at the Earth\'s surface, the rate of
the change can reach 7 nT/month).

\(c\) the optimal length of the spherical harmonic expansions is now
automatically set inside the code as a function of the radial distance,
so that the deviation from the full-length approximation does not exceed
0.01 nT. (In all previous versions, the upper limit on the order of
spherical harmonics had to be explicitly specified by users),

\(2\) Subroutine DIP, calculating the Earth\'s field in the dipole
approximation:

\(a\) no longer accepts the tilt angle via the list of formal
parameters. Instead, the sine SPS and cosine CPS of that angle are now
implicitly forwarded into DIP via the first common block /GEOPACK1/.
Accordingly, there are two options: (i) to implicitly calculate SPS and
CPS by invoking RECALC before calling DIP, or (ii) to specify them
explicitly. In the last case, the common block must be present in the
calling module **and** SPS and CPS should be specified **after** the
invocation of RECALC (otherwise they will be overridden by those
returned by RECALC).

\(b\) the Earth\'s dipole moment is now calculated by RECALC, based on
the table of the IGRF coef-ficients and their secular variation rates,
for a given year and day of the year, and the obtained value of the
moment is forwarded into DIP via the second common block /GEOPACK2/. (In
all previous versions, only a single fixed value was provided for the
geodipole moment, corresponding to the most

recent epoch).

\(3\) Subroutine RECALC: now consolidates in one module all calculations
needed to initialize and up-date the values of coefficients and
quantities that vary in time, either due to secular changes of the main
geomagnetic field or as a result of Earth\'s diurnal rotation and
orbital motion around Sun. That allowed us to simplify the codes and
make them more compiler-independent.

\(4\) Subroutine GEOMAG: is now identical in its structure to other
coordinate transformation subrou-tines. It no longer invokes RECALC from
within GEOMAG, but uses pre-calculated values of the ro-tation matrix
elements, obtained by a separate external invocation of RECALC. This
eliminates possi-ble interference of the two subroutines in the old
version of the package.

\(5\) Subroutine TRACE (and the subsidiary modules STEP and RHAND):

\(a\) no longer needs to specify the highest order of spherical
harmonics in the main geomagnetic field expansion. Instead, it is now

calculated automatically inside the IGRF_GSM (or IGRF_GEO) subroutine.

\(b\) the internal field model can now be explicitly chosen by
specifying the parameter INNAME (either as IGRF_GSM or DIP).

\(6\) A new subroutine BCARSP was added. It converts Cartesian field
components into spherical ones (an operation inverse to that performed
by the subroutine BSPCAR).

\(7\) Two new subroutines were added, SHUETAL_MGNP and T96_MGNP, to
calculate the position of the magnetopause, according to the model of
*Shue et al.* \[1998\] and the one used in the T96 magnetospheric
magnetic field model \[*Tsyganenko*, 1995, 1996\]. The model of *Shue et
al.* provides better accuracy, since it takes into account both the
solar wind ram pressure P and the IMF Bz component, while the T96 model
magnetopause is driven only by P (it provides here a starting
approximation boundary for the subroutine SHUETAL_MGNP).

**GEOPACK: A SET OF FORTRAN SUBROUTINES FOR COMPUTATIONS OF THE**

**GEOMAGNETIC FIELD IN THE EARTH\'S MAGNETOSPHERE**

Version of February 08, 2008

(Previous releases: October 1979, March 1987, April 16, 1996, January 5,
2001, April 22, 2003, May 4, 2005)

N. A. TSYGANENKO

Institute of Physics, University of St.-Petersburg, Petrodvoretz,
Saint-Petersburg 198504, Russian Federation.

I. **INTRODUCTION**

Recent studies in the solar-terrestrial physics led to recognizing the
role of the geomagnetic field as one of the most important
characteristics of human environment. The Earth\'s magnetic field links
the interplanetary medium with the upper atmosphere and ionosphere,
guides energetic charged particles ejected during solar flares, channels
the low-frequency electromagnetic waves and heat flux, confines the
radiation belt and auroral plasma, and serves as a giant accumulator of
the solar wind energy that eventually dissipates during the magnetic
storms. Understanding these phenomena is crucial for forecasting
conditions in the near-Earth geospace ("space weather") that strongly
affect modern space technologies.

In many applications one often needs numerical tools for evaluating the
components of the geomagnetic field in a wide range of distances and
tracing its force lines far away from the Earth\'s surface, calculate
geomagnetically conjugate points, and map a spacecraft\'s position with
respect to characteristic magnetospheric/ionospheric boundaries. This
requires using quantitative models of the Earth\'s magnetic field,
including its internal part due to the dynamo currents inside Earth, and
the external part, produced by the magnetospheric and ionospheric
electric current systems.

The present set of subroutines was conceived as a subsidiary software
package for calculating the geomagnetic field components at any point of
space from the Earth\'s surface up to the Moon\'s orbit. Upon specifying
year, day of year, and universal time as input parameters, it calculates
elements of the matrices of transformations between several most
frequently used coordinate systems. It also updates coefficients of
spherical harmonic expansions, approximating the Earth\'s internal
magnetic field (the IGRF/DGRF models). That field can be computed either
in a dipole approximation or by using a full-scale IGRF/DGRF model, in
which the length of expansions is automatically controlled to maintain
the needed precision. It also contains a field line mapping subroutine,
tracing the geomagnetic field lines from any point of space, based on
appropriate internal and external field models for a specified date/time
and/or the geodipole tilt angle. Like in the previous versions, we did
not include subroutines for calculating the external magnetic field, but
chose to restrict this package to only a general-use set of codes, which
is unlikely to significantly change in the future. In contrast, the
external field models rapidly improve and supersede older versions; for
that reason, it was considered reasonable to provide them separately.

A convenient way of using the GEOPACK-2008 subroutines is to separately
compile them and consolidate the corresponding object modules into a
dedicated library.

The GEOPACK subroutines were originally written in 1978; the present
version emerged as a result of their many upgrades, various tests, and
numerous helpful comments received from users since the first release of
the package. Although it now appears as a quite robust tool, there
should certainly exist room for further improvement. The author would
greatly appreciate any comments on the performance of the package,
reports on related problems, and any suggestions on how to make the
GEOPACK subroutines simpler, faster, more versatile, and easier to
understand.

Below is a sort of a manual guide for using the subroutines, including a
bibliography and two typical samples of a main program for tracing model
field lines, aimed at helping users to compile and debug their own
codes, avoiding common mistakes. The FORTRAN listings of the package
subroutines are placed in a separate file GEOPACK-2008.FOR.

II\. **DESCRIPTIONS OF INDIVIDUAL SUBROUTINES**

**1.** SUBROUTINE: **IGRF_GSW_08**

FUNCTION: Calculates three components of the main (internal) geomagnetic
field in the

Geocentric Solar-Wind (GSW, also known as GSWM) coordinate system.

FORTRAN STATEMENT:

CALL IGRF_GSW_08 (XGSW,YGSW,ZGSW,BXGSW,BYGSW,BZGSW)

INPUT PARAMETERS: XGSW,YGSW,ZGSW - position of the observation point in
Cartesian GSW

coordinates (in Earth\'s radii, R~E~=6371.2 km), at which the field
vector is to

be evaluated. XGSW axis is directed antiparallel to the observed solar
wind

flow vector, ZGSW axis points northward and lies in the plane defined by

XGSW and the Earth's magnetic dipole axis, YGSW axis completes the
right-

handed system. The GSW system reduces to the standard GSM in the case of

strictly radial solar wind flow.

OUTPUT PARAMETERS: BXGSW, BYGSW, BZGSW - Cartesian GSW components of the
main

geomagnetic field in nanotesla.

COMMON BLOCKS: /GEOPACK2/.

OTHER SUBROUTINES INVOKED: GEOGSW_08.

COMMENTS: The internal sources of the geomagnetic field are rigidly tied
to rotating Earth, and hence

their position with respect to Sun varies with universal time and day of
year. In addition,

the geomagnetic field slowly changes with time in the Earth\'s frame of
reference (secular

variation). Because of that, when calculating the geomagnetic field in
the GSW system,

one needs first to define a transformation matrix between the geographic
coordinates and

GSW and initialize/update the IGRF model coefficients, by specifying the
universal time

and date, as well as three components of the solar wind velocity in GSE
coordinates before

the first invocation of this subroutine, or if the time/date were
changed. This should be done

by calling the subroutine RECALC_08, which calculates the position of
Sun and takes into

account secular variations of the internal field. It includes the
harmonic coefficients for 9

epochs: 1965.0, 1970.0, 1975.0, 1980.0, 1985.0, 1990.0, 1995.0, 2000.0,
and 2005.0 so that

their values for any date are calculated by the linear interpolation
between the nearest

epochs. For the dates in the interval 2005\<IYEAR\<2010, the subroutine
uses extrapolated

coefficients, based on the secular velocities through the order 8. If
IYEAR\<1965 or

IYEAR\>2010, the coefficients are assumed equal to those for 1965 or
2005, respectively.

Subroutine RECALC_08 is regularly updated, as coefficients for the next
epoch become

available. When calculating the IGRF field at large distances, there is
no need to retain all

terms in the expansion, since the relative contribution from
higher-order multipoles rapidly

decreases with the altitude. To save the calculation, the subroutine
IGRF_GSW_08

automatically truncates the summation, so that at any radial distance
the error does not

exceed 0.01 nT. For more details on the main geomagnetic field modeling,
see, e.g., *Langel*

\[1987\] and *Tsyganenko* \[1990\].

**2.** SUBROUTINE: **IGRF_GEO_08**

FUNCTION: Calculates three components of the main (internal) geomagnetic
field in

spherical geocentric geographic (GEO) coordinate system.

FORTRAN STATEMENT: CALL IGRF_GEO_08 (R, THETA, PHI, BR, BTHETA, BPHI)

INPUT PARAMETERS: R, THETA, PHI - position in spherical geocentric GEO
coordinates:

R is the radial distance (in Earth\'s radii, R~E~ =6371.2 km), THETA and
PHI are the

colatitude and longitude (in radians), respectively.

OUTPUT PARAMETERS: BR, BTHETA, BPHI - spherical components of the main
geomagnetic

field (in nanotesla; positive BR outward, BTHETA southward, BPHI
eastward).

COMMON BLOCKS: /GEOPACK2/.

OTHER SUBROUTINES INVOKED: None.

COMMENTS: In this case the Earth\'s internal field is calculated in the
geographic (geocentric)

coordinate system, rigidly fixed with respect to Earth. In this system,
there are no

diurnal/seasonal variations of the internal field, and the position of
Sun is irrele-

vant. The only effect to be taken into account here is a slow (secular)
variation of

the internal field. As in the case of IGRF_GSW_08, this should be done
by calling the

subroutine RECALC_08, but only year and day of the year are important
here, while

the universal time of the day and the solar wind velocity do not matter
and can be

set arbitrarily.

**3.** SUBROUTINE: **DIP_08**

FUNCTION: Calculates Cartesian Geocentric Solar-Wind (GSW) components of
the Earth\'s magnetic

field, corresponding to the first (dipolar) term in its spherical
harmonic expansion for a

specified epoch.

FORTRAN STATEMENT: CALL DIP_08 (XGSW,YGSW,ZGSW,BXGSW,BYGSW,BZGSM)

INPUT PARAMETERS: XGSW,YGSW, ZGSW - position in the Cartesian GSW
coordinates (in

Earth\'s radii, R~E~ =6371.2 km).

OUTPUT PARAMETERS: BXGSW, BYGSW, BZGSW - geodipole field components, in
nanotesla.

COMMON BLOCKS: /GEOPACK1/, /GEOPACK2/.

OTHER SUBROUTINES INVOKED: None.

COMMENTS: (1) The Earth\'s dipole is assumed to be centered at the
origin.

\(2\) The geodipole axis is rigidly tied to the rotating Earth, and
hence its position with

respect to Sun varies with universal time and day of year. In addition,
the solar wind

direction fluctuates around the (aberrated) Sun-Earth line, and both the
magnitude and

orientation of the geodipole moment slowly change with time (secular
variation). Because

of that, before calculating the dipole field in the GSW system one
should first determine the

transformation between the geographic and GSW coordinates and
initialize/update the

value of the geodipole moment by specifying the universal time and date.
This should be

done before the first invocation of DIP_08 (or if the time/date/solar
wind direction have

changed) by calling the subroutine RECALC_08, which calculates the
current position of

the solar wind source with respect to the geodipole, taking into account
the secular variation

of the internal field. See also comments for IGRF_GSW_08 and RECALC_08.

**4.** SUBROUTINE: **SUN_08**

FUNCTION: This is a subsidiary subroutine, which is usually called from
the subroutine

RECALC_08. It calculates the orientation of the Earth-Sun line in the
geocentric

inertial coordinate system (GEI), and the Greenwich mean sidereal time.

FORTRAN STATEMENT:

CALL SUN_08 (IYEAR,IDAY,IHOUR,MIN,ISEC,GST,SLONG,SRASN,SDEC)

INPUT PARAMETERS: IYEAR,IDAY,IHOUR,MIN,ISEC - year (four digits), day of
year, hour,

minute, and second, respectively;

OUTPUT PARAMETERS:

GST - Greenwich mean sidereal time,

SLONG - ecliptic longitude of Sun

SRASN - right ascension of Sun

SDEC - declination of Sun

Note: all the above output parameters are **in radians**.

COMMON BLOCKS: None.

OTHER SUBROUTINES INVOKED: None.

COMMENTS: (1) 1901\<IYR\<2099

\(2\) IDAY=1 is January 1

\(3\) This subroutine, authored by G. D. Mead, was compiled with minor

changes from the paper by *Russell* \[1971\].

**5.** SUBROUTINE: **SPHCAR_08**

FUNCTION: Computes components of a position vector in spherical
coordinates from its

Cartesian components or vice versa.

FORTRAN STATEMENT: CALL SPHCAR_08 (R, THETA, PHI, X, Y, Z, J)

INPUT PARAMETERS: (a) J - integer switch parameter

\(b\) if J\>0, then spherical coordinates R,THETA,PHI (colatitude THETA

and longitude PHI are **in radians**)

if J\<0, then Cartesian coordinates X,Y,Z

OUTPUT PARAMETERS: if J\>0 then Cartesian coordinates X,Y,Z

if J\<0 then spherical coordinates R,THETA,PHI.

COMMON BLOCKS: None.

OTHER SUBROUTINES INVOKED: None.

COMMENTS: If X=0 and Y=0, then PHI is set equal to 0 (J\<0).

**6.** SUBROUTINE: **BSPCAR_08**

FUNCTION: Calculates Cartesian components of a field vector from its
known local spherical

components (e.g., returned by IGRF_GEO_08) and spherical angles THETA

and PHI of the observation point.

FORTRAN STATEMENT: CALL BSPCAR_08 (THETA,PHI,BR,BTHETA,BPHI,BX,BY,BZ)

INPUT PARAMETERS: THETA, PHI - colatitude and longitude of the
observation point

(**radians**)

BR, BTHETA, BPHI - spherical components of the field vector in

a local orthogonal coordinate system with its origin at the observation

point.

OUTPUT PARAMETERS: BX,BY,BZ - Cartesian components of the vector.

COMMON BLOCKS: None.

OTHER SUBROUTINES INVOKED: None.

**7.** SUBROUTINE: **BCARSP_08**

FUNCTION: Calculates spherical components of a field vector in a local
orthogonal coordinate

system, based on the known position of the observation point and the
field vector

components in the corresponding Cartesian coordinate system.

FORTRAN STATEMENT: CALL BCARSP_08 (X,Y,Z,BX,BY,BZ,BR,BTHETA,BPHI)

INPUT PARAMETERS: X, Y, Z - position of the observation point in the
Cartesian coordinate system,

BX, BY, BZ - Cartesian field components at that location

OUTPUT PARAMETERS: BR,BTHETA,BPHI - spherical components of the vector
in the **local**

orthogonal coordinate system with its origin at the observation point.

COMMON BLOCKS: None.

OTHER SUBROUTINES INVOKED: None.

**8.** SUBROUTINE: **RECALC_08**

FUNCTIONS: (1) Calculates the angles defining the geodipole axis
orientation for a given year, day of

year, and universal time of day, as well as elements of rotation
matrices, needed for

transformations between the following Cartesian geocentric coordinate
systems:

geographic geocentric (GEO), geomagnetic (MAG), solar-magnetic (SM),
geocentric

solar-wind (GSW), geocentric solar-ecliptic GSE), and geocentric
equatorial inertial (GEI).

\(2\) Initializes or updates the values of the spherical harmonic
coefficients for the model

(IGRF/DGRF) main geomagnetic field expansions, and coefficients of the
recursion

relations, entering in those expansions.

FORTRAN STATEMENT:

CALL RECALC_08 (IYEAR,IDAY,IHOUR,MIN,ISEC,VGSEX,VGSEY,VGSEZ)

INPUT PARAMETERS: IYEAR, IDAY, IHOUR, MIN, ISEC, VGSEX, VGSEY, VGSEZ -
year

(four digits), day, UT hour, minute, second, and 3 GSE components of the
solar wind velocity.

OUTPUT PARAMETERS: None.

OTHER SUBROUTINES INVOKED: SUN_08.

COMMON BLOCKS: Output parameters are placed into two named common
blocks:

/GEOPACK1/ (containing 34 four-byte variables) and

/GEOPACK2/ (containing 198 four-byte variables).

COMMENTS: (1) This subroutine should be called at least once before
using the following

subroutines from this package: IGRF_GEO_08, IGRF_GSW_08, DIP_08,

GEOMAG_08, GEOGSW_08, MAGSM_08, SMGSW_08, GSWGSE_08,

GEIGEO_08, TRACE_08.

\(2\) If only questionable information (or no information at all) is
available on the solar

wind flow velocity, or, if the standard GSM and/or SM coordinates are to
be used,

then set VGSEX=-400.0, VGSEY=0.0, and VGSEZ=0.0. In that case, the
Solar-

Wind (GSW) coordinate system becomes identical to the standard GSM.

**9.** SUBROUTINE: **GEOMAG_08**

FUNCTION: Transformation of components of a vector from the geographic
(geocentric)

Cartesian coordinate system (GEO) to the dipolar geomagnetic (MAG) or
vice

versa.

FORTRAN STATEMENT: CALL GEOMAG_08 (XGEO,YGEO,ZGEO,XMAG,YMAG,ZMAG,J)

INPUT PARAMETERS: (a) J - integer switch parameter,

\(b\) if J\>0, then geographic coordinates XGEO, YGEO, ZGEO

if J\<0, then geomagnetic coordinates XMAG, YMAG, ZMAG

OUTPUT PARAMETERS: If J\>0 then XMAG, YMAG, ZMAG; if J\<0 then XGEO,
YGEO, ZGEO.

COMMON BLOCKS: /GEOPACK1/

OTHER SUBROUTINES INVOKED: None.

COMMENTS: Subroutine RECALC_08 must be called before GEOMAG_08 in two
cases:

\(a\) before the first use, or

\(b\) if the current values of IYEAR and/or IDAY were changed

No information is required in this case on the solar wind velocity, so
one can set

VGSEX=-400.0, VGSEY=0.0, VGSEZ=0.0 in RECALC_08.

**10.** SUBROUTINE: **GEIGEO_08**

FUNCTION: Transformation of components of a vector from the Geocentric
Equatorial

Inertial (GEI) coordinate system into Geographic geocentric (GEO) or

vice versa.

FORTRAN STATEMENT: CALL GEIGEO_08 (XGEI,YGEI,ZGEI,XGEO,YGEO,ZGEO,J)

INPUT PARAMETERS: (a) J - integer switch parameter:

\(b\) if J\>0, then GEI coordinates XGEI, YGEI, ZGEI

if J\<0, then the Geographic coordinates XGEO, YGEO, ZGEO

OUTPUT PARAMETERS: If J\>0 then XGEO, YGEO, ZGEO; if J\<0 then XGEI,
YGEI, ZGEI.

COMMON BLOCKS: /GEOPACK1/

OTHER SUBROUTINES INVOKED: None.

COMMENTS: Subroutine RECALC_08 must be called before GEIGEO_08 in two
cases:

\(a\) before the first use, or

\(b\) if the current values of IYEAR, IDAY, IHOUR, MIN, ISEC have
changed

No information is required here on the solar wind velocity, so one can
set in RECALC_08

VGSEX=-400.0, VGSEY=0.0, VGSEZ=0.0.

**11.** SUBROUTINE: **MAGSM_08**

FUNCTION: Converts components of a vector from the Dipolar Geomagnetic
(MAG)

coordinate system into the Solar Magnetic (SM) system or vice versa.

FORTRAN STATEMENT: CALL MAGSM_08 (XMAG,YMAG,ZMAG,XSM,YSM,ZSM,J)

INPUT PARAMETERS: (a) J - integer switch parameter:

\(b\) if J\>0, then MAG coordinates XMAG, YMAG, ZMAG

if J\<0, then SM coordinates XSM, YSM, ZSM

OUTPUT PARAMETERS: If J\>0 then XSM, YSM, ZSM; if J\<0 then XMAG, YMAG,
ZMAG.

COMMON BLOCKS: /GEOPACK1/

OTHER SUBROUTINES INVOKED: None.

COMMENTS: 1. Subroutine RECALC_08 must be called before MAGSM_08 in
three cases:

\(a\) before the first use, or

\(b\) if the current values of IYEAR, IDAY, IHOUR, MIN, ISEC have
changed, and/or

\(c\) if the solar wind velocity components VGSEX, VGSEY, VGSEZ have
changed

2\. A general implicit assumption here is that the position of the noon
meridian plane in the

SM system is defined by the observed direction of the incoming solar
wind flow, rather

than by the Earth-Sun line. To convert vectors to/from the standard SM
system, set in

RECALC_08 VGSEX=-400.0, VGSEY=0.0, and VGSEZ=0.

**12.** SUBROUTINE: **GSWGSE_08**

FUNCTION: Transformation of components of a vector from the Geocentric
Solar-Wind (GSW)

coordinate system into the Geocentric Solar Ecliptic (GSE) system or
vice versa.

FORTRAN STATEMENT: CALL GSWGSE_08 (XGSW,YGSW,ZGSW,XGSE,YGSE,ZGSE,J)

INPUT PARAMETERS: (a) J - integer switch parameter:

\(b\) if J\>0 then GSW coordinates XGSW, YGSW, ZGSW

if J\<0 then GSE coordinates XGSE, YGSE, ZGSE

OUTPUT PARAMETERS: If J\>0 then XGSE, YGSE, ZGSE; if J\<0 then XGSW,
YGSW, ZGSW.

COMMON BLOCKS: /GEOPACK1/

OTHER SUBROUTINES INVOKED: None.

COMMENTS: 1. Subroutine RECALC_08 must be called before GSWGSE_08 in
three cases:

\(a\) before the first use, or

\(b\) if the current values of IYEAR, IDAY, IHOUR, MIN, ISEC have
changed, and/or

\(c\) if the solar wind velocity components VGSEX, VGSEY, VGSEZ have
changed

2\. Setting VGSEX=-400.0, VGSEY=0.0, and VGSEZ=0.0 in RECALC_08 results
in

reduction of all GSW vector components to those in the standard GSM; in
that case, this

subroutine becomes fully equivalent to GSMGSE from previous releases of
this package.

**13.** SUBROUTINE: **SMGSW_08**

FUNCTION: Transformation of components of a vector from the Solar
Magnetic (SM) coordinate

system into the Geocentric Solar-Wind (GSW) system and vice versa.

FORTRAN STATEMENT: CALL SMGSW_08 (XSM,YSM,ZSM,XGSW,YGSW,ZGSW,J)

INPUT PARAMETERS: (a) J - integer switch parameter:

\(b\) if J\>0 then SM coordinates XSM, YSM, ZSM

if J\<0 then GSW coordinates XGSW, YGSW, ZGSW

OUTPUT PARAMETERS: If J\>0 then XGSW, YGSW, ZGSW; if J\<0 then XSM, YSM,
ZSM.

COMMON BLOCKS: /GEOPACK1/

OTHER SUBROUTINES INVOKED: None.

COMMENTS: 1. Subroutine RECALC_08 must be called before SMGSW_08 in
three cases:

\(a\) before the first use, or

\(b\) if the current values of IYEAR, IDAY, IHOUR, MIN, ISEC have
changed, and/or

\(c\) if the solar wind velocity components VGSEX, VGSEY, VGSEZ have
changed

2\. A general implicit assumption here is that the position of the noon
meridian plane in the

SM system (as well as the orientation of X axis in GSW) are defined by
the observed

direction of the incoming solar wind flow, rather than by the Earth-Sun
line. To convert

vectors between standard SM and GSM systems, set in RECALC_08:
VGSEX=-400.0

(or any other non-zero value), VGSEY=0.0, and VGSEZ=0.0.

**14.** SUBROUTINE: **GEOGSW_08**

FUNCTION: Transformation of components of a vector from the Geographic
Geocentric (GEO)

coordinates into the Geocentric Solar-Wind (GSW) system and vice versa.

FORTRAN STATEMENT: CALL GEOGSW_08 (XGEO,YGEO,ZGEO,XGSW,YGSW,ZGSW,J)

INPUT PARAMETERS: (a) J - integer switch parameter:

\(b\) if J\>0 then GEO coordinates XGEO, YGEO, ZGEO

if J\<0 then GSW coordinates XGSW, YGSW, ZGSW

OUTPUT PARAMETERS: If J\>0 then XGSW, YGSW, ZGSW; if J\<0 then XGEO,
YGEO, ZGEO.

COMMON BLOCKS: /GEOPACK1/

OTHER SUBROUTINES INVOKED: None.

COMMENTS: 1. Subroutine RECALC_08 must be called before GEOGSW_08 in
three cases:

\(a\) before the first use, or

\(b\) if the current values of IYEAR, IDAY, IHOUR, MIN, ISEC have
changed, and/or

\(c\) if the solar wind velocity components VGSEX, VGSEY, VGSEZ have
changed

2\. Setting VGSEX=-400.0, VGSEY=0.0, and VGSEZ=0.0 in RECALC_08 results
in

reduction of all GSW vector components to those in the standard GSM; in
that case, this

subroutine becomes fully equivalent to GEOGSM from previous releases of
this package.

**15.** SUBROUTINE: **GEODGEO_08**

FUNCTION: Converts geodetic coordinates {H,****} of a point in space
(local vertical height, and geodetic latitude, respectively) to
geocentric coordinates {R,} (radial distance and spherical
co-latitude, respectively) and vice versa.

FORTRAN STATEMENT: CALL GEODGEO_08 (H,XMU,R,THETA,J)

INPUT PARAMETERS: (a) J - integer switch parameter:

\(b\) if J\>0 then geodetic coordinates H (km), XMU (radians)

if J\<0 then geocentric coordinates R (km), THETA (radians)

OUTPUT PARAMETERS: If J\>0 then R (km), THETA (radians); if J\<0 then H
(km), XMU (radians).

COMMON BLOCKS: None.

OTHER SUBROUTINES INVOKED: None.

**16.** SUBROUTINE: **RHAND_08**

FUNCTION: Subsidiary subroutine, calculating the right-hand side vector
in the field line

equation (that is, components of the unit vector along the local
direction of **B**)

FORTRAN STATEMENT:

CALL RHAND_08 (X,Y,Z,R1,R2,R3,IOPT,PARMOD,EXNAME,INNAME)

INPUT PARAMETERS:

X, Y, Z - GSW coordinates of the current point on a field line;

IOPT - integer index, reserved for specifying a version of the external
magnetic field

model. In the case of T87 or T89 model \[*Tsyganenko*, 1987; 1989\],
this parameter can

be used for specifying the range of the Kp-index;

PARMOD - a 10-element array, reserved for other input parameters of the
external

field model, e.g., solar wind pressure, interplanetary magnetic field
components, Dst-

index, as in the T96 model \[*Tsyganenko*, 1995; 1996\];

EXNAME - name of an external field model subroutine (see below the
description of

TRACE_08 for more details).

INNAME - name of an internal field model subroutine (see below
description of

TRACE_08 for more details).

OUTPUT PARAMETERS: R1, R2, R3 -- right hand side quantities (components
of the unit vector

**B**/B).

COMMON BLOCKS: /GEOPACK1/

OTHER SUBROUTINES INVOKED: EXNAME and INNAME (replaced by actual names
of

external and internal model field subroutines).

**17.** SUBROUTINE: **STEP_08**

FUNCTION: Makes one step along a model magnetic field line. The model
field is a sum of

contributions from internal and external sources, provided by the
subroutines

IGRF_GSW_08 (or DIP_08) and EXNAME, respectively.

FORTRAN STATEMENT:

CALL STEP_08 (X, Y, Z, DS, DSMAX, ERRIN, IOPT, PARMOD, EXNAME, INNAME)

INPUT PARAMETERS:

X, Y, Z - GSW position of the initial point on a field line;

DS - tentative step size along the line;

DSMAX -- upper limit on the step size

ERRIN - upper limit of the permissible error;

IOPT - integer index, reserved for specifying a version of the external
magnetic field

model (in the case of T87 or T89 model this parameter can be used for
speci-

fying the range of the Kp-index)

PARMOD - a 10-element array, reserved for other input parameters of the
external

field model (e.g., solar wind pressure, interplanetary magnetic field
compo-

nents, and Dst, as in the T96 model)

EXNAME - name of an external field model subroutine (see below
description of the

subroutine TRACE_08 for more details)

INNAME - name of an internal field model subroutine (see below
description of

TRACE_08 for more details).

OUTPUT PARAMETERS: X, Y, Z - GSW coordinates of the next point on the
field line: **note that the**

**input values of X,Y,Z are replaced by their output values**,
corresponding to the

next point on the line.

COMMON BLOCKS: /GEOPACK1/

OTHER SUBROUTINES INVOKED: RHAND_08.

COMMENTS: (1) Do not forget to specify the actual names of the external
and internal model field

subroutines in the EXTERNAL statement of your own subroutine or main
program, from

which STEP_08 is invoked (see also comments 1-3 in the description of
TRACE_08).

\(2\) The parameter DS can change its value as a result of having
invoked STEP_08,

because (a) if the actual error estimate exceeds ERRIN, the step is
repeated with

a halved value of the stepsize DS, (b) if the actual error is too small
(less than

0.04\*ERRIN), the current step is **not** repeated, and the value of DS
for the next

step is doubled (but not to exceed DSMAX), in anticipation of a further
decrease of the

field line curvature.

**18.** SUBROUTINE: **TRACE_08**

FUNCTION: Calculates GSW positions of a sequence of points lying on a
field line. The tracing starts at

a given initial point and terminates (i) upon crossing inner spherical
boundary of the tracing

region (Earth\'s surface or the ionosphere), (ii) upon exiting from the
modeling region, (iii) if

the field line makes more than 2 loops inside the tracing region around
an O-type quasi-

neutral line, or (iv) if a prescribed upper limit on the number of field
line points is exceeded.

FORTRAN STATEMENT:

CALL TRACE_08 (XI,YI,ZI,DIR,DSMAX,ERR,RLIM,R0,IOPT,

\* PARMOD,EXNAME,INNAME,XF,YF,ZF,XX,YY,ZZ,L,LMAX)

INPUT PARAMETERS:

XI, YI, ZI, - GSW position of the initial point (in Earth\'s radii,
R~E~);

DIR - defines the tracing direction: assign DIR=1.0 if the line should
be traced oppositely to

the **B** vector (e.g., from the northern polar ionosphere down the
northern tail lobe), or

DIR=-1.0 in the opposite case.

DSMAX -- a desired value (in R~E~) of the upper limit on the step size
along the field line.

ERR -- permissible error at each step. A reasonable estimate that
ensures a sufficient accuracy for

most applications is ERR=0.0001. Smaller/larger values will result in
larger/smaller number

of steps (hence, more field line points). Note that using much smaller
values of this parameter

(such as 10^-10^) may not necessarily result in a tangible improvement
of the accuracy, unless

the entire source code (including this package) is converted to double
precision.

RLIM - geocentric radius of the outer boundary of the tracing region (in
R~E~); tracing is

terminated once R\>RLIM;

R0 - radius of the inner boundary of the tracing region (usually R0=1.0,
corresponding to

the Earth\'s surface);

IOPT - an integer index, reserved for specifying a version of the
external magnetic field

model (in the case of T87 or T89 model, this parameter can be used for
specifying

the range of the Kp-index);

PARMOD - a 10-element array, reserved for other input parameters of the
external field

model (e.g., solar wind pressure, interplanetary magnetic field
components, and

Dst, as in the T96 model);

EXNAME - name of an external field model subroutine. The corresponding
actual argu-

ment should be specified in the EXTERNAL statement of the main module.
The

list of formal parameters of the subroutine EXNAME must be as follows:

(IOPT, PARMOD, PS, X, Y, Z, BX, BY, BZ), where IOPT is an integer,
PARMOD

is a 10-element array of REAL\*4 type, and the remaining parameters are
real

variables, standing for the dipole tilt angle (in radians), GSW position
(in Re) and

the external field components.

INNAME - name of a subroutine providing components of the internal
magnetic field

(either IGRF_GSW_08 or DIP_08).

LMAX -- maximal length of the arrays XX, YY, ZZ in which coordinates of
the field line

points are stored. LMAX should be set equal to the actual length of the
arrays,

defined in the main program as actual arguments of this subroutine

OUTPUT PARAMETERS:

XX, YY, ZZ -- arrays containing GSW coordinates (in R~E~) of the field
line points;

XF, YF, ZF -- coordinates of the field line ending point;

L -- the number of calculated field line points.

COMMON BLOCKS: /GEOPACK1/

OTHER SUBROUTINES INVOKED: STEP_08, RHAND_08.

COMMENTS:

\(1\) This subroutine allows two options:

\(a\) If INNAME=IGRF_GSW_08, then the IGRF model will be used for
calculating the

contribution from earth\'s internal sources. In this case, subroutine
RECALC_08 must be

called before using TRACE_08, with properly specified date, universal
time, and solar

wind velocity components, to calculate in advance all quantities needed
for the main

field model and for transformations between involved coordinate systems.

\(b\) If INNAME=DIP_08, then a pure dipole field will be used instead of
the IGRF model.

As in the case (a), the subroutine RECALC_08 must also be called before
TRACE_08.

Here we can choose either to (i) calculate dipole tilt angle based on
date, time, and

solar wind direction, or (ii) explicitly specify that angle, without any
reference to

date/time/solar wind. In the last case, the sine (SPS) and cosine (CPS)
of the

dipole tilt angle must be specified in advance (**but after** **the
invocation of** RECALC_08)

and forwarded in the common block /GEOPACK1/ (in its 11th and 12th
elements,

respectively). In this case the role of the subroutine RECALC_08 is
reduced to only

initializing the components of the Earth\'s dipole moment.

\(2\) Any external field model subroutine EXNAME can be used, provided
it has the same struc-

ture of the list of formal arguments and its actual name is included in
the EXTERNAL

statement in the user\'s module that invokes TRACE_08.

\(3\) If the field line is closed (i.e., returns back to Earth), then
the last (L-th) element in the arrays

XX, YY, ZZ (and the triplet XF, YF, ZF) will correspond to the footpoint
of the line at R=R0.

\(4\) If in the process of the tracing the number of field line points
exceeds LMAX, computation

is aborted and a warning message displayed.

For more details, see examples of using this subroutine in the end of
this file.

**19.** SUBROUTINE: **SHUETAL_MGNP_08**

FUNCTION: For a given point of space (either inside or outside of the
magnetosphere) finds a point on

the model magnetopause by *Shue et al*. \[1998\], lying in the plane
containing XGSM-axis

and the original point, and nearly closest to that point. The model
magnetopause is specified

by setting parameters of the incoming solar wind and Bz component of
IMF.

FORTRAN STATEMENT:

CALL SHUETAL_MGNP (XN_PD, VEL, BZIMF, XGSM, YGSM, ZGSM,

\*XMGNP, YMGNP, ZMGNP, DIST, ID)

INPUT PARAMETERS: XN_PD - either the solar wind proton number density in
particles per cubic

centimeter (if VEL\>0), or the solar wind ram pressure in nanoPascals
(if VEL\<0)

VEL - either the solar wind speed (km/s) or any negative value, which in

that case indicates that XN_PD should be interpreted as the solar wind

ram pressure in nPa.

BZIMF - Bz component of the interplanetary magnetic field.

XGSM, YGSM, ZGSM - GSM position of the input observation point (in
R~E~).

OUTPUT PARAMETERS:

XMGNP, YMGNP, ZMGNP - GSM position of the output boundary point

DIST - distance (in R~E~) between the input point (XGSM, YGSM, ZGSM) and
the model

magnetopause

ID - position flag: ID=+1 (-1) means that the input point lies inside
(outside) of the model

magnetopause.

COMMON BLOCKS: None.

OTHER SUBROUTINES INVOKED: T96_MGNP_08.

COMMENTS: Position vectors (XGSM, YGSM, ZGSM) and (XMGNP,YMGNP,ZMGNP)
can be

interpreted as referring to the GSW system (assuming that the
magnetopause orients itself

following the direction of the solar wind flow).

**20.** SUBROUTINE: **T96_MGNP_08**

FUNCTION: For a given point of space (either inside or outside of the
magnetosphere) finds a point on

the model magnetopause used in the T96 model \[*Tsyganenko*, 1996\],
lying on the plane

containing XGSM axis and the original point, and nearly closest to that
point. The model

magnetopause is specified by setting parameters of the incoming solar
wind.

FORTRAN STATEMENT:

CALL T96_MGNP_08 (XN_PD, VEL, XGSM, YGSM, ZGSM, XMGNP, YMGNP, ZMGNP,
DIST, ID)

INPUT PARAMETERS:

XN_PD - either the solar wind proton number density in particles per
cubic cm (if VEL\>0),

or the solar wind ram pressure in nPa (if VEL\<0).

VEL - either the solar wind speed (km/s) or any negative value, which in
this case

indicates that XN_PD should be interpreted as the solar wind ram
pressure.

BZIMF - Bz component of the interplanetary magnetic field.

XGSM, YGSM, ZGSM - GSM position of the input observation point (in Re).

OUTPUT PARAMETERS:

XMGNP, YMGNP, ZMGNP - GSM position of the output point on the boundary.

DIST - distance (in Re) between the observation point (XGSM, YGSM, ZGSM)
and the

model magnetopause

ID - position flag: ID=+1 (-1) means that the input point lies inside
(outside) of the model

magnetopause.

COMMON BLOCKS: None.

OTHER SUBROUTINES INVOKED: None.

COMMENTS: Position vectors (XGSM, YGSM, ZGSM) and (XMGNP,YMGNP,ZMGNP)
can be

interpreted as referring to the GSW system (assuming that the
magnetopause orients itself

following the direction of the solar wind flow).

III\. **REFERENCES**

Hones, E. W., Jr., R. D. Zwickl, T. A. Fritz, and S. J. Bame (1986),
Structural and dynamical

aspects of the distant magnetotail determined from ISEE-3 plasma
measurements, Planet. Space

Sci., v.34, 889, 1986.

Langel, R. A., Main field, Ch.IV, in Geomagnetism, edited by J. A.
Jacobs, Academic

Press, London, 1987.

Russell, C. T., Geophysical coordinate transformations. Cosmic
Electrodynamics, v.2,

p.184, 1971.

Shue, J.-H. et al., Magnetopause location under extreme solar wind
conditions,

J. Geophys. Res., v.103, p. 17,691, 1998.

Tsyganenko, N. A., Global quantitative models of the geomagnetic field
in the cislu-

nar magnetosphere for different disturbance levels, Planet.Space Sci.,
v.35,

p.1347, 1987.

Tsyganenko, N. A., A magnetospheric magnetic field model with a warped
tail plasma

sheet. Planet.Space Sci., v.37, p.5, 1989.

Tsyganenko, N. A., Quantitative models of the magnetospheric magnetic
field: Methods

and results. Space Sci. Rev., v.54, pp.75-186, 1990.

Tsyganenko, N. A., Modeling the Earth\'s magnetospheric magnetic field
confined within

a realistic magnetopause, J. Geophys.Res., v.100, p.5599, 1995.

Tsyganenko, N. A., Effects of the solar wind conditions on the global
magnetospheric

configuration as deduced from data-based field models, in European Space
Agency

Publication ESA SP-389, p.181, 1996.

Tsyganenko, N. A., S. B. P. Karlsson, S. Kokubun, T. Yamamoto, A. J.
Lazarus, K. W. Ogilvie,

and C. T. Russell (1998), Global configuration of the magnetotail
current sheet as derived from

Geotail, Wind, IMP 8, and ISEE 1/2 data, J. Geophys. Res., v.103, 6827.

Tsyganenko, N. A., and D. H. Fairfield, Global shape of the magnetotail
current sheet as derived from

Geotail and Polar data, J. Geophys. Res., v.109(A3), A03218, doi:
10.1029/2003JA010062, 2004.

IV\. **ACKNOWLEDGEMENT:** The author thanks Victor Sergeev for his
useful suggestions, in particular, to include in the current version
transformations that take into account the observed deflection of the
solar wind flow from strictly anti-sunward direction. Thanks are also
due to many other users who provided their helpful comments on the
previous releases of the package.

C

C
\############################################################################

C \# THE MAIN PROGRAMS BELOW GIVE TWO EXAMPLES OF TRACING FIELD LINES \#

C \# USING THE GEOPACK-2008 SOFTWARE (release of Feb 08, 2008) \#

C
\############################################################################

C

PROGRAM EXAMPLE1

C

C IN THIS EXAMPLE IT IS ASSUMED THAT WE KNOW GEOGRAPHIC COORDINATES OF A
FOOTPOINT

C OF A FIELD LINE AT THE EARTH\'S SURFACE AND TRACE THAT LINE FOR A
SPECIFIED

C MOMENT OF UNIVERSAL TIME, USING A FULL IGRF EXPANSION FOR THE INTERNAL
FIELD

C

PARAMETER (LMAX=500)

C

C LMAX IS THE UPPER LIMIT ON THE NUMBER OF FIELD LINE POINTS RETURNED BY
THE TRACER.

C IT CAN BE SET ARBITRARILY LARGE, DEPENDING ON THE SPECIFICS OF A
PROBLEM UNDER STUDY.

C IN THIS EXAMPLE, LMAX IS TENTATIVELY SET EQUAL TO 500.

C

DIMENSION XX(LMAX),YY(LMAX),ZZ(LMAX), PARMOD(10)

c

c Be sure to include an EXTERNAL statement in your codes, specifying the
names

c of external and internal field model subroutines in the package, as
shown below.

c In this example, the external and internal field models are T96_01 and
IGRF_GSW_08,

c respectively. Any other models can be used, provided they have the
same format

c and the same meaning of the input/output parameters.

c

EXTERNAL T96_01,IGRF_GSW_08

C

C DEFINE THE UNIVERSAL TIME AND PREPARE THE COORDINATE TRANSFORMATION
PARAMETERS

C BY INVOKING THE SUBROUTINE RECALC_08: IN THIS PARTICULAR CASE WE TRACE
A LINE

C FOR YEAR=1997, IDAY=350, UT HOUR = 21, MIN = SEC = 0

C

IYEAR=1997

IDAY=350

IHOUR=21

MIN=0

ISEC=0

C

C AT THAT TIME, ACCORDING TO THE OMNI DATABASE, THE SOLAR WIND VELOCITY
IN GSE HAD THE

C COMPONENTS

C

VGSEX=-304.0

VGSEY= -16.0

VGSEZ= 4.0

C

C NOTE, HOWEVER, THAT THE ABERRATION CORRECTION WAS ALREADY MADE IN THE
SOLAR WIND DATA.

C THEREFORE, TO CORRECTLY TRANSFORM THE DATA TO GSM COORDINATE SYSTEM,
WE HAVE TO RESTORE C VGSEY TO ITS ORIGINAL OBSERVED VALUE:

C

VGSEY=VGSEY+29.78

C

CALL RECALC_08 (IYEAR,IDAY,IHOUR,MIN,ISEC,VGSEX,VGSEY,VGSEZ)

C

OPEN(UNIT=1,FILE=\'LINTEST1.DAT\')

C

WRITE (1,100) IYEAR, IDAY, IHOUR, MIN

100 FORMAT (\' IYEAR=\',I4,\' IDAY=\',I3,\' IHOUR=\',I2,\' MIN=\',I2,/)

C

PARMOD(1)=3.

PARMOD(2)=-20.

PARMOD(3)=3.

PARMOD(4)=-5.

WRITE (1,110) PARMOD(1)

110 FORMAT(\' SOLAR WIND RAM PRESSURE (NANOPASCALS):\',F6.1,/)

C

WRITE (1,120) PARMOD(2)

120 FORMAT (\' DST-INDEX: \',F6.0,/)

C

WRITE (1,130) PARMOD(3),PARMOD(4)

130 FORMAT (\' IMF By and Bz: \',2F6.1,/)

C

C THE LINE WILL BE TRACED FROM A GROUND (RE=1.0) FOOTPOINT WITH
GEOGRAPHICAL

C LONGITUDE 45 DEGREES AND LATITUDE 75 DEGREES:

C

GEOLAT=75.

GEOLON=45.

RE=1.0

PRINT \*, \' GEOGRAPHIC (GEOCENTRIC) LATITUDE (degs): \',GEOLAT

PRINT \*, \' GEOGRAPHIC (GEOCENTRIC) LONGITUDE (degs): \',GEOLON

COLAT=(90.-GEOLAT)\*.01745329

XLON=GEOLON\*.01745329

C

C CONVERT SPHERICAL COORDS INTO CARTESIAN :

C

CALL SPHCAR_08 (RE,COLAT,XLON,XGEO,YGEO,ZGEO,1)

C

C TRANSFORM GEOGRAPHICAL GEOCENTRIC COORDS INTO SOLAR WIND
MAGNETOSPHERIC ONES:

C

CALL GEOGSW_08 (XGEO,YGEO,ZGEO,XGSW,YGSW,ZGSW,1)

C

c SPECIFY TRACING PARAMETERS:

C

DIR=1.

C (TRACE THE LINE WITH A FOOTPOINT IN THE NORTHERN HEMISPHERE, THAT IS,

C ANTIPARALLEL TO THE MAGNETIC FIELD)

C

DSMAX=1.0

C (MAXIMAL SPACING BETWEEN THE FIELD LINE POINTS SET EQUAL TO 1 RE)

C

ERR=0.0001

C (PERMISSIBLE STEP ERROR SET AT ERR=0.0001)

RLIM=60.

C (LIMIT THE TRACING REGION WITHIN R=60 Re)

C

R0=1.

C (LANDING POINT WILL BE CALCULATED ON THE SPHERE R=1,

C I.E. ON THE EARTH\'S SURFACE)

IOPT=0

C (IN THIS EXAMPLE IOPT IS JUST A DUMMY PARAMETER,

C WHOSE VALUE DOES NOT MATTER)

C

C TRACE THE FIELD LINE:

C

CALL TRACE_08 (XGSW,YGSW,ZGSW,DIR,DSMAX,ERR,RLIM,R0,IOPT,

\* PARMOD,T96_01,IGRF_GSW_08,XF,YF,ZF,XX,YY,ZZ,M,LMAX)

C

C WRITE THE RESULTS IN THE DATAFILE \'LINTEST1.DAT\':

C

WRITE (1,20)

20 FORMAT(\' THE LINE IN GSW COORDS:\',/)

WRITE (1,21) (XX(L),YY(L),ZZ(L),L=1,M)

21 FORMAT ((2X,3F6.2))

CLOSE(UNIT=1)

END

C\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*

C

PROGRAM EXAMPLE2

C

C THIS IS ANOTHER EXAMPLE OF USING THE GEOPACK SUBROUTINE \"TRACE_08\".
UNLIKE IN

C THE EXAMPLE1, HERE WE ASSUME A PURELY DIPOLAR APPROXIMATION FOR THE
EARTH\'S

C INTERNAL FIELD.

C IN THIS CASE WE ALSO EXPLICITLY SPECIFY THE TILT ANGLE OF THE
GEODIPOLE,

C INSTEAD OF CALCULATING IT FROM THE DATE/TIME.

C

PARAMETER (LMAX=500)

C

C LMAX IS THE UPPER LIMIT ON THE NUMBER OF FIELD LINE POINTS RETURNED BY
THE TRACER.

C IT CAN BE SET ARBITRARILY LARGE, DEPENDING ON THE SPECIFICS OF A
PROBLEM UNDER STUDY.

C IN THIS EXAMPLE, LMAX IS TENTATIVELY SET EQUAL TO 500.

C

DIMENSION XX(LMAX),YY(LMAX),ZZ(LMAX), PARMOD(10)

c

C Unlike in the EXAMPLE1, here we \"manually\" specify the tilt angle
and its sine/cosine.

c To forward them to the coordinate transformation subroutines, we need
to explicitly

c include the common block /GEOPACK1/:

C

COMMON /GEOPACK1/ AA(10),SPS,CPS,BB(3),PSI,CC(18)

C

c be sure to include an EXTERNAL statement with the names of (i) a
magnetospheric

c external field model and (ii) Earth\'s internal field model.

c

EXTERNAL T96_01, DIP_08

C

C First, call RECALC_08, to define the main field coefficients and,
hence, the magnetic

C moment of the geodipole for IYEAR=1997 and IDAY=350.

C The universal time and solar wind direction does not matter in this
example,

C because here we explicitly specify the tilt angle (hence, the
orientation of

C dipole in the GSW coordinates), so we arbitrarily set IHOUR=MIN=ISEC=0
and

C VGSEX=-400.0, VGSEY=VGSEZ=0 (any other values would be equally OK):

C

CALL RECALC_08 (1997,350,0,0,0,-400.0,0.0,0.0)

c

c Enter input parameters for T96_01:

c

PRINT \*, \' ENTER SOLAR WIND RAM PRESSURE IN NANOPASCALS\'

READ \*, PARMOD(1)

C

PRINT \*, \' ENTER DST \'

READ \*, PARMOD(2)

C

PRINT \*, \' ENTER IMF BY AND BZ\'

READ \*, PARMOD(3),PARMOD(4)

C

c Define the latitude (XLAT) and longitude (XLON) of the field line
footpoint

c in the GSW coordinate system:

c

XLAT=75.

XLON=180.

C

C Specify the dipole tilt angle PS, its sine SPS and cosine CPS,
entering

c in the common block /GEOPACK1/:

C

PSI=0.

SPS=SIN(PSI)

CPS=COS(PSI)

c

c Calculate Cartesian coordinates of the starting footpoint:

c

T=(90.-XLAT)\*.01745329

XL=XLON\*.01745329

XGSW=SIN(T)\*COS(XL)

YGSW=SIN(T)\*SIN(XL)

ZGSW=COS(T)

C

c SPECIFY TRACING PARAMETERS:

C

DIR=1.

C (TRACE THE LINE WITH A FOOTPOINT IN THE NORTHERN HEMISPHERE, THAT IS,

C ANTIPARALLEL TO THE MAGNETIC FIELD)

C

DSMAX=1.0

C (SETS THE MAXIMAL SPACING BETWEEN CONSECUTIVE POINTS ON THE LINE)

ERR=0.0001

C (PERMISSIBLE STEP ERROR SET AT ERR=0.0001)

RLIM=60.

C (LIMIT THE TRACING REGION WITHIN R=60 Re)

C

R0=1.

C (LANDING POINT WILL BE CALCULATED ON THE SPHERE R=1,

C I.E. ON THE EARTH\'S SURFACE)

c

IOPT=0

C (IN THIS EXAMPLE IOPT IS JUST A DUMMY PARAMETER,

C WHOSE VALUE DOES NOT MATTER)

c

c Trace the field line:

c

CALL TRACE_08 (XGSW,YGSW,ZGSW,DIR,DSMAX,ERR,RLIM,R0,IOPT,

\* PARMOD,T96_01,DIP_08,XF,YF,ZF,XX,YY,ZZ,M,LMAX)

C

C Write the result in the output file \'LINTEST2.DAT\':

C

OPEN(UNIT=1, FILE=\'LINTEST2.DAT\')

1 WRITE (1,20) (XX(L),YY(L),ZZ(L),L=1,M)

20 FORMAT((2X,3F6.2))

CLOSE(UNIT=1)

END