##// END OF EJS Templates
DATABASES » madrigal
RSS
Add restart to nginx container
Add restart to nginx container
jespinoza - - Load All Authors

File last commit:

r0:b84e1135c2c4
r22:1c8c96a4f254 master
Show More
sfc_convert_geo_coord.f
902 lines | 33.2 KiB | text/x-fortran | FortranFixedLexer
/ source / madf / geolib / sfc_convert_geo_coord.f
History | Source | Raw |Copy content |Copy permalink
jespinoza
Initial
 r0 SUBROUTINE CONVERT_GEO_COORD(R_LAT_IN,R_LON_IN,R_HEIGHT_IN,
* R_LAT_OUT,R_LON_OUT,I_FLAG,I_ERROR)
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c $Revision: 1.2 $
c The initial version of this subroutine was written by RADEX, INC.
c for use on VAX/VMS systems.
c
c Subsequent revisions for UNIX systems have been made by KBB at
c The Johns Hopkins Univ. Applied Physics Laboratory. These revisions
c have been managed using the Revision Control System (RCS) and a
c log of the revisions will be found at the end of the comments.
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
C
C
C Purpose:
C
C This subroutine uses a set of spherical harmonic coefficients to
C perform a coordinate conversion between geographic and corrected
C geomagnetic coordinates (or vice versa).
C
C The spherical harmonics for the geographic to corrected geomagnetic
C coordinate system correspond to a conversion from geographic to
C the corrected geomagnetic coordinates (expressed in terms of
C centered dipole coordinates at the input altitude), and are then
C transformed to ground level.
C
C The spherical harmonic coefficients used for the inverse are
C computed relative to centered dipole coordinates at the input
C altitude. The input CGM coordinates are converted to equivalent
C values in this coordinate sytem using the inverse altitude
C algorithm before the evaluation of the geographic spherical
C harmonic expansion.
C
C Method:
C
C This subroutine uses a five-step process in converting a position
C in one coordinate system to a position in the other.
C The five steps are as follows:
C
C 1. The appropriate spherical harmonic coefficients for the
C coordinate conversion are computed for the input altitude.
C
C 2. The appropriate coordinates for use in the spherical harmonic
C expansion are computed. For the geographic ==> corrected
C geomagnetic coordinate version, these are geographic colatitude
C and longitude. For the inverse coordinate conversion, the
C input coordinates are first converted into equivalent dipole
C coordinates at the input altitude.
C
C 3. The Cartesian coordinates of the unit vector in the desired
C coordinate system are then computed using the appropriate
C spherical harmonic expansion.
C
C 4. For the geographic ==> corrected geomagnetic coordinate
C conversion, the dipole-equivalent coordinates at input
C altitude are converted to their cgm cartesian coordinates.
C
C 5. Standard trigonometric identities are used to compute the
C latitude and longitude in the desired output coordinate system.
C
C Input Arguments:
C
C R_LAT_IN - REAL*4 - The input latitude in degrees.
C This could be either geographic latitude
C or corrected geomagnetic latitude. The
C acceptable range of values is (-90 to +90
C degrees).
C
C R_LON_IN - REAL*4 - The input longitude in degrees east.
C This could be either geographic longitude
C or corrected geomagnetic longitude. The
C acceptable range of values is (0 to 360 degrees).
C
C R_HEIGHT_IN - REAL*4 - The height in kilometers for which the
C coordinate transformation will be accomplished.
C The acceptable range of values is (0 km to
C 7200 km)
C
C I_FLAG - INTEGER - The flag that indicates which way
C the conversion will proceed.
C
C = 1 convert geographic to corrected
C geomagnetic coordinates
C
C = 2 convert corrected geomagnetic
C to geographic coordinates
C
C Output Arguments:
C
C R_LAT_OUT - REAL*4 - The output latitude in degrees.
C This could be either the geographic latitude
C or corrected geomagnetic latitude. This
C value will be between -90 and +90 degrees.
C
C R_LON_OUT - REAL*4 - The output longitude in degrees.
C This could be either geographic longitude or
C corrected geomagnetic longitude. This
C value will be between 0 and 360 degrees.
C
C I_ERROR - INTEGER - The error flag
C
C = 0 normal processing
C = -2 R_HEIGHT_IN is outside the allowable range
C = -4 I_FLAG value is invalid
C = -8 R_LAT_IN is outside the allowable range
C = -16 R_LON_IN is outside the allowable range
C = -32 Magnitude of the "unit vector" of the
C target coordinate system deviates
C significantly (+/- 20% or more) from 1.
C = -64 For altitudes > 0, for the corrected
C geomagnetic to geographic coordinate
C conversion there is a range of latitudes
C for which the transformation is invalid.
C This flag is set when the requested input
C cgm latitude falls within the invalid
C region.
C
C Local Variables:
C
C C_CLASS - Character String, used to identify message
C class in error messages.
C
C C_FLAG - Character String used to store I_FLAG
C as a character string for use in error
C message.
C
C C_HEIGHT_IN - Character String used to store R_HEIGHT_IN
C as a character string for use in error
C message.
C
C C_LAT_IN - Character String used to store R_LAT_IN
C as a character string for use in error
C message.
C
C C_LON_IN - Character String used to store R_LON_IN
C as a character string for use in error
C message.
C
C C_R - Character String used to store D_R
C as a character string for use in error
C message.
C
C C_REP_PARAMS - Character String used to store error
C message.
C
C
C C_ROUTINE_NAME - Character String containing name of
C subroutine generating error message.
C
C C_UNDEFINED - Character String used to store input data
C (R_HEIGHT, R_LAT_IN, R_LON_IN) as a
C character string for use in error message.
C
C D_CINT(,,) - Double Precision Array (3-D) -
C Contains the spherical harmonic coefficients
C interpolated to the input height, R_HEIGHT_IN.
C
C D_COEF(,,,) - Double Precision Array (3-D) -
C Contains the spherical harmonic coefficients
C used to compute the Cartesian coordinates of
C unit vector in the target coordinate system.
C
C First index: sp. harm. coeff. index
C Second x, y, z components of unit vector
C Third altitude indices
C Fourth direction of conversion index
C
C coefficients for a given altitude have the form
C
C a0 + h a1 + h^2 * a2 + h^3 * a3
C where h = alt/7000 [km]
C
C D_COLAT_INPUT - Double Precision - colatitude (radians) in the
C input coordinate system
C
C D_COLAT_OUTPUT - Double Precision - colatitude (radians) in the
C output coordinate system
C
C D_LON_INPUT - Double Precision - longitude (radians) in the
C input input system
C
C D_LON_OUTPUT - Double Precision - longitude (radians) in the
C output coordinate system
C
C D_R - Double Precision - magnitude of the
C unit radius vector in the target coordinate
C system. The target coordinate system is
C the system to which this subroutine is
C converting the latitude and longitude.
C
C D_X - Double Precision - the X-component of the
C unit radius vector in the target coordinate
C system
C
C D_Y - Double Precision - the Y-component of the
C unit radius vector in the target coordinate
C system
C
C D_Z - Double Precision - the Z-component of the
C unit radius vector in the target coordinate
C system
C
C D_YLMVAL - Double Precision - the array of spherical
C harmonic basis functions evaluated at
C a particular colatitude and longitude.
C
C
C I_COND_VALUE - Integer - the condition value returned from the
C call to SFC$$PUT_USER_MSG
C
C I_ERR64 - Logical - error flag for CGM_TO_ALTITUDE
C Routine.
C
C
C I_SYSTEM_ERR - Integer - Set to the value, SFC$$_NORMAL in
C all calls to SFC$$PUT_USER_MSG
C
C K - Integer - first index of the D_CINT
C interpolated spherical harmonic coefficient
C array.
C
C L - Integer - order of each spherical harmonic
C coefficient and spherical harmonic function.
C
C M - Integer - zonal index of the spherical
C harmonic functions.
C
C R_HEIGHT_OLD() - Real*4 Array (1-D) - Variable containing
C previous height (km) used to determine whether
C the interpolation to compute D_CINT() from
C D_COEF() needs to be done.
C
C R_LAT_ADJ - Real*4 Variable used to store result of
C ALTITUDE_TO_CGM or CGM_TO_ALTITUDE
C conversion.
C
C R_LAT_ALT - Real*4 Variable used to store at altitude
C dipole latitude.
C
C R_R - Real*4 - magnitude of unit vector (= D_R)
C (D_X, D_Y, D_Z) in single precision
C
C Constants:
C
C
C I_MAX_LENGTH - Integer - the length of the D_COEF(,,,)
C and D_YLMVAL spherical harmonic function
C arrays, corresponding to the order of the
C spherical harmonic expansion used.
C I_MAX_LENGTH = (I_ORDER + 1) * (I_ORDER + 1)
C
C I_MSGID - Integer - contains a message shell number of
C 4014 which will be used in all calls to
C SFC$$PUT_USER_MSG
C
C I_NUM_AXES - Integer - the number of axes in a Cartesian
C coordinate system (3)
C
C I_NUM_FLAG - Integer - the number of coordinate trans-
C formation flags available (2)
C
C I_NUM_LEVEL - Integer - the number of grid levels
C available. This is the number of
C distinct heights for which there are
C spherical harmonic coefficients available
C
C I_ORDER - Integer - the order of the spherical harmonic
C expansion used.
C
C DEGRAD - Double Precision - the multiplicative
C conversion factor for converting degrees
C to radians
C
C PI - Double Precision - mathematical constant
C
C Subroutines Required:
C
C RYLM - This subroutine returns the spherical harmonic
C function value array for a given colatitude and
C longitude. The spherical harmonics are returned
C in the D_YLMVAL array.
C
C CG_ALT_DIP - Given altitude (km), latitude, and a direction
C flag, returns altitude corrected latitude, and
C an error flag. See comments in subroutine for
C additional information.
C
C
C SFC$$PUT_USER_MSG - This subroutine is used to send error messages
C to the Message Handling System.
C
C Files Used at Compile Time:
C
C SFCLIB$$DEFS$$:SFC$$SFCDEF/LIST
C
C - This include file contains useful system
C constants and definitions. (eg. this is where
C SFC$$_NORMAL is defined)
C
C
C Files Used at Run Time: None
C
C
C Databases Accessed: None
C
C
C Warnings:
C
C 1. For certain values of Corrected Geomagnetic Coordinates, the
C transformation to Geographic Coordinates (Geocentric) is
C undefined. When this situation occurs, an error flag is set
C I_ERROR = -64 and returned to the calling routine.
C
C
C 2. This subroutine contains the non-standard include statement
C
C
C
C
C Revision History:
C
C Original Routine, based upon the Kyle Baker routine was converted
C to a standard library subroutine for SFC$$ by MFS on 14 OCT 1992.
C
C 09/17/93 SPR1993-0203 (SWS) The return condition value needs to
C be initialized before coordinate conversion processing begins
C
C 08/18/93 SPR1993-0175 (SWS) Routine performed a comparison of
C Z with a -1.D0. The variable Z should be D_Z.
C
C The present routine represents the use of an improved algorithm for
C the computation of corrected geomagnetic coordinates and (where
C it exists) the inverse, together with the use of an updated
C magnetic field model (IGRF 1995).
C
C The present routine incorporates the functionality of Kyle Baker
C routine, and, uses the same subroutine calls. Relevent portions of
C the old FORTRAN code have been retained where possible.
C
C The new algorithm and the current code have been developed by
C Radex, Inc, 3 Preston Court, Bedford, MA 01730. A technical
C report describing the improved algorithm, and discussing its
C limitations is in preparation, and will be published in early
C 1997 as a Philips Laboratory technical report entitled:.
C
C An Expanded Algorithm for the Computation of Corrected Geomagnetic
C Coordinates, by C. Hein and K. Bhavnani.
C
C The principal changes implementing the new algorithm and code are
C as follows:
C
C (1) The set of spherical harmonic function values are computed
C recursively using a single call to a new subroutine RYLM
C rather than using repetitive function call. The Kyle Baker
C routine calculated each of the spherical harmonics separately
C as needed.
C
C Computationally the recursive caluclation is much faster.
C
C (2) The spherical harmonic coefficients are based upon
C the IGRF 95 magnetic field model. The order of the
C spherical harmonic expansion used here is 10 (previously
C was 4).
C
C (3) The coefficients represent a cubic fit to the spherical
C harmonic coefficients over a range of 24 altitudes from 0 to
C 7200 km. The previous (7/17/95) version used a quadratic fit
C at 0, 300 and 1200 km altitude, and had an allowable
C range of altitudes is 0 - 2000 km. The original Kyle Baker
C routine used polynomial interpolations of the coefficients
C at 0, 150, 300 and 450 km.
C
C (4) The spherical harmonic coefficients are computed in an
C auxilliary coordinate system, which is aligned with the
C desired target coordinate system. This results in a
C considerable simplification of the code, eliminating the
C need to perform multiple coordinate system rotations
C as was the case in the Kyle Baker routine.
C
C 07/17/95
C
C MODIFIED FOR THE IGRF 95 MAGNETIC FIELD MODEL.
C
C CHANGES: BLOCKDATA SECTION HAS BEEN REPLACED WITH COEFFICIENTS
C APPROPRIATE FOR THE IGRF 1995 MODEL
C
C REFERENCES TO IGRF 90 IN THE COMMENT LINES HAVE BEEN
C CHANGED TO IGRF 95.
C
C 07/01/96
C
C The BLOCKDATA section has been expanded to accomodate spherical
C harmonic coefficients for a cubic interpolation, in order to
C accomodate the expanded 0 - 7200 km altitude range. The
C spherical harmonic coefficients were caluclated on the basis
C of the IGRF 95 magnetic field model.
C
C **********************************************************************
C
C The version provided here does not include the BLOCK DATA source
C code. Instead individual BLOCK DATA files are provided as follows
C for the following Magnetic Field Models:
C
C Field Model File Name
c
C DGRF 1975 BLKDAT75.FOR
C DGRF 1980 BLKDAT80.FOR
C DGRF 1985 BLKDAT85.FOR
C DGRF 1990 BLKDAT90.FOR
C IGRF 1995 BLKDAT95.FOR
C
C Line 3 of the BLKDATxx.FOR files contains an identification field
C for the magnetic field model used to generate the coefficients.
C
C **********************************************************************
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c RCS revision log for POSIX systems
c
c $Log: not supported by cvs2svn $
c Revision 1.1 2004/07/08 18:01:56 brideout
c files added to support AACGM conversion
c
c Revision 1.8 1997/11/14 19:27:49 baker
c implementation of modifications to extend the coordinate
c system to 7200 km in altitude.
c
c Revision 1.7 1997/10/28 21:04:24 baker
c There was an error in the value of PI that has been fixed.
c
c Revision 1.6 1996/03/20 22:53:03 baker
c Explicitly declared the sizes of the arguments to the
c subroutine. This should help avoid compatibility
c problems when mixing C and Fortran with different
c machines and different compilers.
c
c Revision 1.5 1996/03/12 18:59:18 baker
c Added code to force a recomputation of the
c conversion coefficients at a given ghheight if the
c coordinates model has changed from the last call.
c
c Revision 1.4 1996/03/11 19:25:36 baker
c Modifications for the 1995 version of AACGM.
c
c
c Revision 1.3 94/10/17 12:35:32 12:35:32 baker (Kile Baker S1G)
cadded error code -64 to indicate invalid magnetic coordinates specified
c as input. This also requires a change in the call to cg_alt_dip
c
c Revision 1.2 94/10/14 10:53:36 10:53:36 baker (Kile Baker S1G)
c Added the SAVE instruction to make variables static
c
c Revision 1.1 94/10/12 15:28:38 15:28:38 baker (Kile Baker S1G)
c Initial revision
c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
C
C
C
C
C Set the order used in this algorithm to 10.
C
C
C Parameterize the maximum values of the indices of the
C spherical harmonic coefficient array D_COEF( , , , ) and the
C interpolated spherical harmonic coefficient array, D_CINT( , , ).
C
C
C
C
C
C
c
c The APL version of this code needs to keep track of
c the previous set of coefficients defining the conversion
c model (date). This is done by looking at the first coefficient
c in the model.
c
C
C CHARACTER*300 C_REP_PARAMS
C CHARACTER*8 C_LAT_IN
C CHARACTER*8 C_LON_IN
C CHARACTER*8 C_HEIGHT_IN
C CHARACTER*2 C_FLAG
C CHARACTER*6 C_R
C CHARACTER*70 C_UNDEFINED
C
C .. Parameters ..
DOUBLE PRECISION PI
PARAMETER (PI=3.141592653589793D0)
DOUBLE PRECISION DEGRAD
PARAMETER (DEGRAD=1.745329251994330D-2)
INTEGER I_ORDER
PARAMETER (I_ORDER=10)
INTEGER I_NUM_TERMS
PARAMETER (I_NUM_TERMS=121)
INTEGER I_NUM_AXES
PARAMETER (I_NUM_AXES=3)
INTEGER I_NUM_LEVEL
PARAMETER (I_NUM_LEVEL=5)
INTEGER I_NUM_FLAG
PARAMETER (I_NUM_FLAG=2)
C ..
C .. Scalar Arguments ..
DOUBLE PRECISION R_HEIGHT_IN,R_LAT_IN,R_LAT_OUT,R_LON_IN,R_LON_OUT
INTEGER I_ERROR,I_FLAG
C ..
C .. Local Scalars ..
DOUBLE PRECISION ALT_VAR,ALT_VAR_CU,ALT_VAR_QU,ALT_VAR_SQ,
* D_COLAT_INPUT,D_COLAT_OUTPUT,D_COLAT_TEMP,
* D_LON_INPUT,D_LON_OUTPUT,D_LON_TEMP,D_R,D_X,D_Y,
* D_Z,FIRST_COEFF_OLD,R_LAT_ADJ,R_LAT_ALT,R_R
INTEGER I,J,K,L,M
LOGICAL I_ERR64
C ..
C .. Local Arrays ..
DOUBLE PRECISION D_CINT(I_NUM_TERMS,I_NUM_AXES,I_NUM_FLAG),
* D_YLMVAL(I_NUM_TERMS),R_HEIGHT_OLD(2)
C ..
C .. External Subroutines ..
EXTERNAL ALTITUDE_TO_CGM,CGM_TO_ALTITUDE,RYLM
C ..
C .. Intrinsic Functions ..
INTRINSIC ABS,DACOS,DATAN2,DBLE,SQRT
C ..
C .. Common blocks ..
COMMON /SPH_HARM_MODEL/D_COEF
DOUBLE PRECISION D_COEF(I_NUM_TERMS,I_NUM_AXES,I_NUM_LEVEL,
* I_NUM_FLAG)
C ..
C .. Save statement ..
SAVE D_CINT,R_HEIGHT_OLD,FIRST_COEFF_OLD
C ..
C .. Data statements ..
C
C Initialize R_HEIGHT_OLD() to impossible values.
C
C
C
C
DATA R_HEIGHT_OLD(1)/-1.D0/
DATA R_HEIGHT_OLD(2)/-1.D0/
C ..
C **** DATA I_SYSTEM_ERR /SFC$$_NORMAL/
C
C
C
C The COMMON block / SPH_HARM_MODEL / array D_COEF contains
C the spherical harmonic coefficients used to generate the
C Cartesian coordinates of the unit vector in the target coordinate
C system. The target coordinate system is the system to which
C the algorithm is converting and is the coordinate system
C corresponding to the output variables, R_LAT_OUT and R_LON_OUT.
C
C The coefficient set is stored in a 4-dimensional array with
C indices defined as follows:
C
C First Index - Represents the number of terms used in the
C spherical harmonic expansion. Equal to
C (I_ORDER + 1) * (I_ORDER + 1)
C
C Second Index - Represents the Cartesian coordinate a particular
C coefficient will be used to generate. Indices are
C defined as follows:
C
C 1 - X-coordinate (the X-axis points from the
C center of the earth to where the prime
C meridian crosses the equator.
C
C 2 - Y-coordinate (the Y-axis points from the
C center of the earth to 90 degrees east
C of the X-axis)
C
C 3 - Z-coordinate (the Z-axis points from the
C center of the earth to the north pole)
C
C Third Index - Represents the terms of a quadratic fit to
C altitude (independent variable h =
C altitude [km]/1200) for a given spharical
C harmonic coefficient. The indices are defined
C as follows:
C
C 1 - Constant term: corresponds to the fit at
C 0 km altitude
C 2 - linear term
C 3 - quadratic term
C 4 - cubic term
C
C Fourth Index - Represents the direction of the coordinate
C transformation. Indices are defined as follows:
C
C 1 - Conversion of geographic coordinates to
C corrected geomagnetic coordinates.
C
C 2 - Conversion of corrected geomagnetic coordinates
C to geographic coordinates.
C
C
C The Data for D_COEF is provided in a BLOCK DATA file (see below)
C
C
C Initialize the error return condition indicator
C
I_ERROR = 0
C
C The APL version of this program allows the magnetic field
c model to change during the execution of the program.
c We therefore have to check to make sure that if the model
c has been changed that we recalculate the coordinate
c conversion coefficients.
C
C This IF statement checks to see if the magnetic
C coordinates model has been changed. If so, we
C have to force the recalculation of the conversion
C coefficients, by resetting the values of the old height
C
IF (FIRST_COEFF_OLD.NE.D_COEF(1,1,1,1)) THEN
R_HEIGHT_OLD(1) = -1
R_HEIGHT_OLD(2) = -1
END IF
FIRST_COEFF_OLD = D_COEF(1,1,1,1)
C
C
C The following IF-THEN-ELSE block checks the input arguments to
C ensure they are within allowable ranges. If any argument is
C outside the acceptable range, the error flag, I_ERROR will be
C set to some non-zero value, and control will be returned to the
C calling program.
C
IF ((R_HEIGHT_IN.LT.0.D0) .OR. (R_HEIGHT_IN.GT.7200.D0)) THEN
C
C The height in kilometers is outside the allowable range.
C
C Convert the height value R_HEIGHT_IN to a character string
C and send an error message to the user.
C
C WRITE(UNIT = C_HEIGHT_IN,
C $ FMT = '(F8.2)',
C $ IOSTAT = I_STATUS) R_HEIGHT_IN
C
C C_REP_PARAMS = C_ROUTINE_NAME//'\\HEIGHT\\'//C_HEIGHT_IN
C PRINT *, C_REP_PARAMS
C CALL PUT_USER_MSG(I_MSGID, C_CLASS, C_REP_PARAMS,
C $ I_SYSTEM_ERR, I_COND_VALUE)
C
C Set the error flag to the appropriate value and return
C
I_ERROR = -2
RETURN
C
ELSE IF ((I_FLAG.LT.1) .OR. (I_FLAG.GT.2)) THEN
C
C The conversion flag is neither 1 nor 2.
C
C Convert I_FLAG to a character string C_FLAG and send an error
C message to the user
C
C WRITE(UNIT = C_FLAG,
C $ FMT = '(I2)',
C $ IOSTAT = I_STATUS) I_FLAG
C
C C_REP_PARAMS = C_ROUTINE_NAME//'\\CONVERSION FLAG\\'//C_FLAG
C PRINT*,C_REP_PARAMS
C CALL PUT_USER_MSG(I_MSGID, C_CLASS, C_REP_PARAMS,
C $ I_SYSTEM_ERR, I_COND_VALUE)
C
C Set the error flag to the appropriate value and return
C
I_ERROR = -4
RETURN
C
ELSE IF (ABS(R_LAT_IN).GT.90.D0) THEN
C
C The latitude is outside the allowable range.
C
C Convert R_LAT_IN to a character string C_LAT_IN and send
C an error message to the user.
C
C WRITE(UNIT = C_LAT_IN,
C $ FMT = '(F8.2)',
C $ IOSTAT = I_STATUS) R_LAT_IN
C
C C_REP_PARAMS = C_ROUTINE_NAME//'\\LATITUDE\\'//C_LAT_IN
C PRINT*,C_REP_PARAMS
C CALL PUT_USER_MSG(I_MSGID, C_CLASS, C_REP_PARAMS,
C $ I_SYSTEM_ERR, I_COND_VALUE)
C
C Set the error flag to the appropriate value and return
C
I_ERROR = -8
RETURN
C
ELSE IF ((R_LON_IN.LT.0.D0) .OR. (R_LON_IN.GT.360.D0)) THEN
C
C The longitude is outside the allowable range.
C
C Convert R_LON_IN into a character string C_LON_IN and send
C an error message to the user.
C
C WRITE(UNIT = C_LON_IN,
C $ FMT = '(F8.2)',
C $ IOSTAT = I_STATUS) R_LON_IN
C
C C_REP_PARAMS = C_ROUTINE_NAME//'\\LONGITUDE\\'//C_LON_IN
C PRINT*,C_REP_PARAMS
C CALL PUT_USER_MSG(I_MSGID, C_CLASS, C_REP_PARAMS,
C $ I_SYSTEM_ERR, I_COND_VALUE)
C
C Set the error flag to the appropriate value and return
C
I_ERROR = -16
RETURN
C
END IF
C
C All input arguments are within allowable ranges.
C
C Compute Spherical Harmonic Coefficients for current
C altitude if required
C
IF (R_HEIGHT_IN.NE.R_HEIGHT_OLD(I_FLAG)) THEN
C
ALT_VAR = R_HEIGHT_IN/7200.0D0
ALT_VAR_SQ = ALT_VAR*ALT_VAR
ALT_VAR_CU = ALT_VAR*ALT_VAR_SQ
ALT_VAR_QU = ALT_VAR*ALT_VAR_CU
C
DO 20 I = 1,3
C
DO 10 J = 1,121
C
D_CINT(J,I,I_FLAG) = D_COEF(J,I,1,I_FLAG) +
* ALT_VAR*D_COEF(J,I,2,I_FLAG) +
* ALT_VAR_SQ*D_COEF(J,I,3,I_FLAG) +
* ALT_VAR_CU*D_COEF(J,I,4,I_FLAG) +
* ALT_VAR_QU*D_COEF(J,I,5,I_FLAG)
C
10 CONTINUE
C
20 CONTINUE
C
R_HEIGHT_OLD(I_FLAG) = R_HEIGHT_IN
C
END IF
C
C Zero Sums for Spherical Harmonic Expansion Computation
C
D_X = 0.0D0
D_Y = 0.0D0
D_Z = 0.0D0
C
C
C Prepare for Spherical Harmonic Expansion Computation
C
C
D_LON_INPUT = DBLE(R_LON_IN)*DEGRAD
C
IF (I_FLAG.EQ.1) THEN
C
C Computing CGM from Geographic Coordinates. No altitude
C correction required
C
D_COLAT_INPUT = (90.0D0-DBLE(R_LAT_IN))*DEGRAD
C
ELSE
C
C Computing Geographic Coordinates from CGM Coordinates.
C
C Convert CGM Latitude to Dipole Latitude at Altitude R_HEIGHT
C
CALL CGM_TO_ALTITUDE(R_HEIGHT_IN,R_LAT_IN,R_LAT_ADJ,I_ERR64)
C
IF (I_ERR64 .EQV. .TRUE.) THEN
C
C WRITE(UNIT = C_UNDEFINED,
C $ FMT = 90,
C $ IOSTAT = I_STATUS) R_LAT_IN, R_LON_IN, R_HEIGHT_IN
C
C C_REP_PARAMS = C_ROUTINE_NAME//C_UNDEFINED
C PRINT*,C_REP_PARAMS
C CALL PUT_USER_MSG(I_MSGID, C_CLASS, C_REP_PARAMS,
C $ I_SYSTEM_ERR, I_COND_VALUE)
C
I_ERROR = -64
RETURN
END IF
C
D_COLAT_INPUT = (90.0D0-DBLE(R_LAT_ADJ))*DEGRAD
C
END IF
C
C Generate the spherical harmonics at the coordinate point
C
CALL RYLM(D_COLAT_INPUT,D_LON_INPUT,I_ORDER,D_YLMVAL)
C
C Calculate the Cartesian coordinates of the unit vector
C in the target coordinate system.
C
DO 40 L = 0,I_ORDER
DO 30 M = -L,L
C
K = L*(L+1) + M + 1
C
C Add the contribution of each spherical harmonic to the
C Cartesian components of the unit vector in the
C appropriate coordinate system.
C
C
D_X = D_X + D_CINT(K,1,I_FLAG)*D_YLMVAL(K)
D_Y = D_Y + D_CINT(K,2,I_FLAG)*D_YLMVAL(K)
D_Z = D_Z + D_CINT(K,3,I_FLAG)*D_YLMVAL(K)
C
30 CONTINUE
40 CONTINUE
C
C Compute the magnitude of the Cartesian unit vector of the
C magnetic dipole coordinate system.
C
D_R = SQRT(D_X*D_X+D_Y*D_Y+D_Z*D_Z)
C
C If the magnitude of the unit vector differs significantly
C from 1, set the error flag and continue processing.
C
IF ((D_R.LT.0.9D0) .OR. (D_R.GT.1.1D0)) THEN
C
C The magnitude of the target unit vector differs significantly
C from 1. Convert the value D_R into the character string C_R,
C send a routine message to the user, and continue processing.
C
R_R = D_R
C
C WRITE(UNIT = C_R,
C $ FMT = '(F6.2)',
C $ IOSTAT = I_STATUS) R_R
C
C C_REP_PARAMS = C_ROUTINE_NAME//'\\UNIT VECTOR\\'//C_R
C PRINT*,C_REP_PARAMS
C CALL PUT_USER_MSG(I_MSGID, C_CLASS, C_REP_PARAMS,
C $ I_SYSTEM_ERR, I_COND_VALUE)
C
C Set the error flag to the appropriate value and return.
C
I_ERROR = -32
RETURN
END IF
C
C Adjust the components so they do represent the components of
C a unit vector. If D_R is equal to 1.0D0, this step will not
C change the values of D_X, D_Y, or D_Z.
C
D_Z = D_Z/D_R
D_X = D_X/D_R
D_Y = D_Y/D_R
C
C Obtain output co_latitude and longitude from the unit
C vector using standard formulas
C
IF (D_Z.GT.1.0D0) THEN
D_COLAT_TEMP = 0.0D0
ELSE IF (D_Z.LT.-1.0D0) THEN
D_COLAT_TEMP = PI
ELSE
D_COLAT_TEMP = DACOS(D_Z)
END IF
C
IF ((ABS(D_X).LT.1.0D-8) .AND. (ABS(D_Y).LT.1.0D-8)) THEN
D_LON_TEMP = 0.0D0
ELSE
D_LON_TEMP = DATAN2(D_Y,D_X)
END IF
C
C Prepare Latitude Data for Output
C
D_LON_OUTPUT = D_LON_TEMP
C
IF (I_FLAG.EQ.1) THEN
C
R_LAT_ALT = 90.0D0 - D_COLAT_TEMP/DEGRAD
C
CALL ALTITUDE_TO_CGM(R_HEIGHT_IN,R_LAT_ALT,R_LAT_ADJ)
C
D_COLAT_OUTPUT = (90.0D0-DBLE(R_LAT_ADJ))*DEGRAD
C
ELSE
C
D_COLAT_OUTPUT = D_COLAT_TEMP
C
END IF
C
C Convert colatitude into latitude and convert both coordinates
C from DOUBLE PRECISION radians to REAL degrees.
C
R_LAT_OUT = 90.0D0 - D_COLAT_OUTPUT/DEGRAD
R_LON_OUT = D_LON_OUTPUT/DEGRAD
C
IF (R_LON_OUT.LT.0.0D0) R_LON_OUT = R_LON_OUT + 360.D0
C
C
RETURN
9000 FORMAT (' INVERSE UNDEFINED AT CGM LAT ',F6.2,' CGM LON ',F7.2,
* ' ALTITUDE ',F8.2)
END