Previous: C API   Up: Madrigal developer's guide   Next: Tcl API

Madrigal Fortran API

Madrigal contains a comprehensive set of Fortran-language procedures for working with Madrigal Database files. This part of the Fortran API is simply a very thin wrapper around the C API. Like the C API, this Fortran API is split into two modules, a high level module Maddata, which hides the difference between derived and measured data from the user, and a low-level module Madrec, which allows users to work directly with Cedar files. Examples of using the Maddata and Madrec modules are also given.

The Fortran API also includes the geolib library, which includes methods dealing with geometry and time, which is written completely in Fortran.

The following are suggested lines to add to your Makefile when using the Madrigal Fortran API:


# Library directory
LIBDIR = $(MADROOT)/lib

LDLIBS = -L$(LIBDIR) -lmadrec -lgeo -lm -lnsl

    if solaris:

LDFLAGS = -R$(LIBDIR)

    if gnu:

LDFLAGS = -Xlinker -R$(LIBDIR)


Maddata module


/************************************************************************************************
      
      Fortran callable C routines for Fortran access to the C-language
      Madrigal API - Maddata module.
      
      These comments show how to call these methods from Fortran.
 
      To build, use -L$MADROOT/lib -lmadrec -lgeo
       See example program source/madc/testF77/tmaddataF77.f
      ________________________________________________________________________
      
      SUBROUTINE CRMADD(FILEC,
     *                  PARMS,
     *                  FLTSTR,
     *                  RFMADD,
     *                  NUMREC,
     *                  STATUS)
      CHARACTER FILEC*(*)
      CHARACTER PARMS*(*)
      CHARACTER FLTSTR*(*) 
      INTEGER RFMADD, NUMREC, STATUS
      (for 64 bit machines, use INTEGER*8 RFMADD)
C
      CRMADD  Creates Maddata given a file, a list of desired parameters,
      and filters.   
      CRMADD is meant to stand for "CReate MADData"
C            
      Input parameters: 
      FILEC - file name 
      PARMS - comma delimited list of desired parameters
              (not case-sensitive)  Example PARMS: "gdalt,Azm,F10.7,PH+,Fof2"
      FLTSTR - The filter string is the same string that is used in the new isprint
       command line.  Filters are separated by spaces.  The allowed filters
       are:  
C
          date1=mm/dd/yyyy  (starting date to be examined. If time1 not given, defaults to 0 UT.)
             Example: date1=01/20/1998 
C
          time1=hh:mm:ss (starting UT time to be examined. If date1 given, is applied to date1.
                          If not, applies on the first day of the experiment.)
             Example: time1=13:30:00
C
          date2=mm/dd/yyyy (ending date to be examined.  If time2 not given, defaults to 0 UT.)
             Example: date2=01/21/1998
C
          time2=hh:mm:ss (ending UT time to be examined - If date2 not given, ignored.)
             Example: time2=15:45:00
C
          In the follow arguments ranges are used.  If any range value is not given, it may be used to 
          indicate no lower or upper limit (but the comma is always required). Ranges are inclusive
          of the end points:
C
          z=lower alt limit1, upper alt limit1 [or lower alt limit2 , upper alt limit2 ...] (km)
             Example 1: z=100,500  (This would limit the geodetic altitude to 100 to 500 km.)
             Example 2: z=100,200or300,400  (This would limit the geodetic altitude to 100 to 200 km
                                             or 300 to 400 km.)
             Example 3: z=,200or300,400   (Since the lower limit of the first range is missing, this 
                                           would limit the geodetic altitude to anything below 200 km 
                                           or from 300 to 400 km.)
C
          az=lower az limit1, upper az limit1 [or lower az limit2 , upper az limit2 ...] (from -180 to 180 degrees)
             Example 1: az=100,120  (This would limit the azimuth to 100 to 120 degrees.)
             Example 2: az=-180,-90or90,180  (This would limit the azimuth to between -180 and -90 degrees or
                                             to between 90 and 180 degrees.  Note this allows a filter to go
                                             through 180 degrees.)
C
          el=lower el limit1, upper el limit1 [or lower el limit2 , upper el limit2 ...] (from 0 to 90) 
             Example 1: el=0,45  (This would limit the elevation from 0 to 45 degrees.) 
  
          plen=lower pl limit1, upper pl limit1 [or lower pl limit2 , upper pl limit2 ...] (pulse len in sec)
             Example 1: el=,5e-4  (This would limit the pulse length to 5e-4 seconds or less.)
C
C
          Free form filters using any mnemonic, or two mnemonics added, subtracted, multiplied, or divided.
          Any number of filters may be added: 
C
          filter=[mnemonic] or [mnemonic1,[+*-/]mnemonic2] , lower limit1 , upper limit1 [or lower limit2 , upper limit2 ...] 
             Example 1: filter=ti,500,1000or2000,3000   (Limits the data to points where Ti is between 500 and 1000 degrees
                                                       or between 2000 and 3000 degrees.  Note that the units are always
                                                       those of the Cedar standard.)
             Example 2: filter=gdalt,-,sdwht,0,    (This filter implies "gdalt - sdwht" must be greater than 0.0.  Since
                                                    sdwht is shadow height - the distance above any point on the earth 
                                                    where the sun is first visible - this filter implies that only data 
                                                    in direct sunlight will be displayed.)
             Example 3: filter=ti,/,Dti,100,   (Limits the data to points where the ratio Ti/dTi is more than 100.)
C
          So an full FLTSTR argument might be:
C 
             "date1=01/20/1998 time1=13:30:00 z=,200or300,400 filter=gdalt,-,sdwht,0, filter=ti,/,Dti,100," 
C   
      Output parameters:  
      RFMADD - a long integer reference to the Maddata created, 
      used by all other methods.  Will return 0 if error  
      occurs.  (for 64 bit machines, use INTEGER*8 RFMADD) 
      NUMREC - number of records in Maddata
      STATUS - 0 if success, -1 if failure  

      ***************************************************************
      
      SUBROUTINE CRNFMD(PARMS,
     *                  UT,
     *                  KINST,
     *                  ONED,
     *                  TWOD,
     *                  RFMADD,
     *                  NROW)
      CHARACTER PARMS*(*)
      DOUBLE PRECISION UT
      CHARACTER ONED*(*) 
      CHARACTER TWOD*(*)
      INTEGER KINST, RFMADD, NROW
C     (for 64 bit machines, use INTEGER*8 RFMADD)
C
      CRNFMD  Creates a Maddata record given some input data (no file needed) for a set time.   
      CRNFMD is meant to stand for "CReate NonFile MadData"
C            
      Input parameters: 
C
      PARMS - comma delimited list of desired parameters
              (not case-sensitive)  Example PARMS: "gdalt,Azm,F10.7,PH+,Fof2"
C
      UT - seconds since 1/1/1950 to calculate data at
C
      KINST - instrument id.  If not important, use 0    
C
      ONED - A string that sets one dimension data (That is, one value per parameter).  For example, 
             "gdalt=100.0 glon=45.0 gdlat=-20.0"   
C
      TWOD - A string that sets two dimension data (That is, NROW values per parameter, where every
             parameter must have the same number of values, if more than one).  For example, the
             follow TWOD string would produce 8 rows with every combination of gdlat = 45 or 50,
             glon = 20 or 30, and gdalt = 500 or 600:
C
             "gdlat=45,45,45,45,50,50,50,50 glon=20,20,30,30,20,20,30,30 gdalt=500,600,500,600,500,600,500,600"
C
      Output parameters:
C   
      RFMADD - a long integer reference to the Maddata created, 
      used by all other methods.  Will return 0 if error  
      occurs.  (for 64 bit machines, use INTEGER*8 RFMADD)
C
      NROW - number of rows in Record returned.  If TWOD string used, should be equal to the number of data points
             for each parameter.  If no TWOD data, should be 1.  If error, will be 0.
C
      ***************************************************************
 
 
      SUBROUTINE GTNROW(RFMADD, RECNUM, NROW)
      INTEGER RFMADD
      INTEGER RECNUM, NROW
      (for 64 bit machines, use INTEGER*8 RFMADD)
C 
      GTNROW  Gets the number of rows of data in record RECNUM.   
      GTNROW is meant to stand for "GeT Number of ROWs"
C              
      Input parameters: 
      RFMADD - reference to data created by CRMADD 
      RECNUM - record number of interest (starts at 1)
C 
      Output parameter:
      NROW - number of rows of data.  If only 1D data requested, will
      always be 1.

 
      ***************************************************************
 
      SUBROUTINE GTMADD(RFMADD, RECNUM, NROW, DATROW, STATUS)
      INTEGER RFMADD
      (for 64 bit machines, use INTEGER*8 RFMADD)
      INTEGER RECNUM
      INTEGER NROW, STATUS
      DOUBLE PRECISION DATROW(*)
C 
      GTMADD  Gets one row of data via the DATROW array.   
      GTMADD is meant to stand for "GeT MADData"
C 
      Input parameters: 
      RFMADD - reference to data created by CRMADD 
      RECNUM - record number of interest (starts at 1)
      NROW - row number of interest (starts at 1)
C 
      Output parameters:
      DATROW - Array of doubles to be filled with data.  Order of
      data is the same as the order of parameters in parms string
      passed into CRMADD.  User must be sure array size is equal to
      (or larger than) number of parameters requested.
C
      STATUS - 0 if success, -1 if failure

 
      ***************************************************************
 
 
      SUBROUTINE FRMADD(RFMADD)
      INTEGER RFMADD
      (for 64 bit machines, use INTEGER*8 RFMADD)
C 
      FRMADD  Frees the memory allocated by CRMADD.   
      FRMADD is meant to stand for "FRee MADData"
C 
      Input parameters: 
      RFMADD - reference to data created by CRMADD 
C 
      This method should be called if your program wants to call
      CRMADD more than once.  Call FRMADD just before CRMADD to
      avoid a memory leak.  Will set RFMADD = 0. (for 64 bit machines, use INTEGER*8 RFMADD)
      
      ***************************************************************
      
      
      SUBROUTINE STALOC(KINST,SLATGD,SLON,SALTGD,SLATGC,SR,ISTAT)
      INTEGER KINST,ISTAT
      DOUBLE PRECISION SLATGD,SLON,SALTGD,SLATGC,SR
      
C 
      STALOC  Returns location of a given instrument (kinst)   
      STALOC is meant to stand for "STAtion LOCation"
C 
      Input parameters: 
      KINST - instrument id 
C 
      Output parameters:
      SLATGD - instrument geodetic latitude
      SLON - instrument longitude (0 to 360)
      SALTGD - instrument geodetic altitude (km)
      SLATGC- instrument geocentric latitude
      SR - distance from center of earth (km)
      ISTAT - 0 if successful, -1 if not

*****************************************************************************/

Madrec Module


/************************************************************************************************

Fortran callable C routines for Fortran access to the C-language
Madrigal API - madrec module.
________________________________________________________________________

INTEGER FUNCTION MDOPEN(IOTYPE, FILEC)
   INTEGER IOTYPE
   CHARACTER FILEC(128)
C
   Opens CEDAR file for sequential reading or writing.
C
   IOTYPE -  file type as described below
       Open Cedar file for sequential reading:
             1 - Determine file type automatically
            10 - Madrigal file
            11 - Blocked Binary file
            12 - Cbf file
            13 - Unblocked Binary file
            14 - Ascii file
       Create Cedar file for update; discard previous contents
          if any:
             2 - Madrigal file
            20 - Madrigal file
            21 - Blocked Binary file
            22 - Cbf file
            23 - Unblocked Binary file
            24 - Ascii file
       Create Cedar file in memory for sequential and random read and write. 
            30 - Determine file type automatically
            40 - Madrigal file
            41 - Blocked Binary file
            42 - Cbf file
            43 - Unblocked Binary file
            44 - Ascii file 
       Fast create Cedar file in memory for sequential and random read and write.
          Does not calculate min and and max parameter values 
            50 - Determine file type automatically
            60 - Madrigal file
            61 - Blocked Binary file
            62 - Cbf file
            63 - Unblocked Binary file
            64 - Ascii file
            
   FILEC - file name
C
   Returns file handle (0-9) or negative value if file cannot be opened.
   Call MDGERR to get error description. At most 10 files can be open
   simultaneously.
________________________________________________________________________

INTEGER FUNCTION MDCLOS(MADFIL)
  INTEGER MADFIL

     Closes CEDAR file.

  MADFIL - File handle

  Returns: 0 - File closed successfully
           1 - Error closing file
           2 - File not open
           3 - Error flushing file
________________________________________________________________________

INTEGER FUNCTION MDREAD(MADFIL)
  INTEGER MADFIL

  Reads next CEDAR record

  MADFIL - File handle

  Returns: 0 - Record read successfully
           1 - Illegal file type
          -n - Error
________________________________________________________________________

INTEGER FUNCTION MDWRIT(MADFIL)
  INTEGER MADFIL

  Writes next Madrigal record

  MADFIL - File handle

  Returns: 0 - Record written successfully
           1 - Illegal file type
          -n - Error
          
________________________________________________________________________

INTEGER FUNCTION REWIND(MADFIL)
  INTEGER MADFIL

  Resets a file in memory to point to first record

  MADFIL - File handle

  Returns: 0 - Success
           1 - Illegal file type
          -n - Error
          
________________________________________________________________________

INTEGER FUNCTION GRECNO(MADFIL, RECNO)
  INTEGER MADFIL
  INTEGER RECNO

  Finds a given record number in a file (GRECNO stands for
  Get record by record number)

  MADFIL - File handle
  RECNO  - Record number (1 <= RECNO <= Total number of records)

  Returns: 0 - Success
          -1 - Specified record not in file

________________________________________________________________________

INTEGER FUNCTION GRCKEY(MADFIL, KEY)
  INTEGER MADFIL
  DOUBLE PRECISION KEY

  Finds a given record number in a file using the time key
  The record is the first data record for which key is
  greater than or equal to the start key of the
  record, and less than the start time of the
  following record. Thus, if the specified key
  corresponds to a time within a record, the
  first such record is returned. Header or catalog
  records are never returned. (GRCKEY stands for
  Get record by key)

  MADFIL - File handle
  KEY  - time in seconds since 1/1/1950

  Returns: 0 - Success
          -1 - Specified record not in file

_________________________________________________________________

INTEGER FUNCTION MDCOPY(MADFL1, MADFL2)
  INTEGER MADFL1
  INTEGER MADFL2

  Copies a record from one file to another

  MADFL1 - File handle of source record
  MADFL2 - File handle of destination record

  Returns: 0 - Success
           1 - Failure
________________________________________________________________________

INTEGER FUNCTION MDISDR(MADFIL)
  INTEGER MADFIL

  Identifies data records

  MADFIL - File handle

  Returns: 0 - Current record is not a data record
           1 - Current record is a data record
           
________________________________________________________________________

DOUBLE PECISION FUNCTION GTPMIN(MADFIL, PARCOD)
  INTEGER MADFIL
  INTEGER PARCOD

  Returns minimum value in file of given parcode

  MADFIL - File handle
  PARCOD - Parameter code

  Returns: Scaled minimum parameter value.  File must
  be opened in memory (types 30-44).  If not found,
  returns missing (1E-38).  Data rows with all error
  parameters invalid are not counted.
  
  
________________________________________________________________________

DOUBLE PECISION FUNCTION GTPMAX(MADFIL, PARCOD)
  INTEGER MADFIL
  INTEGER PARCOD

  Returns maximum value in file of given parcode

  MADFIL - File handle
  PARCOD - Parameter code

  Returns: Scaled maximum parameter value.  File must
  be opened in memory (types 30-44)  If not found,
  returns missing (1E-38).  Data rows with all error
  parameters invalid are not counted.

________________________________________________________________________

SUBROUTINE MDCREA(MADFIL,
            LPROL, JPAR, MPAR, NROW,
            KREC, KINST, KINDAT,
            YEAR1, MONTH1, DAY1,
            HOUR1, MIN1, SEC1, CSEC1,
            YEAR2, MONTH2, DAY2,
            HOUR2, MIN2, SEC2, CSEC2)    
  MADFIL - File handle
  INTEGER MADFIL,LPROL, JPAR, MPAR, NROW,KREC, KINST, KINDAT,
          YEAR1, MONTH1, DAY1, HOUR1, MIN1, SEC1, CSEC1,
          YEAR2, MONTH2, DAY2, HOUR2, MIN2, SEC2, CSEC2

    Creates a madrigal record with the specified size and prolog. The
    1d and 2d parameters must be inserted later by calls to
    MDS1DP and MDS2DP.

       MADFIL - File handle
       LPROL   - Length of prolog
       JPAR    - Number of parameters in 1d array
       MPAR    - Number of parameters in 2d array
       NROW    - Number of rows in 2d array
       KREC    - Kind of record
       KINST   - Instrument code
       KINDAT  - Kind-of-data code
       YEAR1-CSEC1 - Start date and time
       YEAR2-CSEC2 - End date and time
________________________________________________________________________


SUBROUTINE MDG1DP(MADFIL, PARCOD, PARM)
  INTEGER MADFIL, PARCOD
  DOUBLE PRECISION PARM

  Gets a 1d parameter from the current record.

  MADFIL - File handle
  PARCOD - Parameter code
  PARM    - Parameter value
________________________________________________________________________

INTEGER FUNCTION MDGNRW(MADFIL) {
  INTEGER MADFIL

  Gets number of rows in 2d array.

  MADFIL - File handle

  Returns: Number of rows in 2d array.
________________________________________________________________________

INTEGER FUNCTION MDGKST(MADFIL) {
  INTEGER MADFIL

  Gets Kind of instrument (Kinst) integer.
  MDGKST stands for MaDrec Get KinST

  MADFIL - File handle

  Returns: Kind of instrument (Kinst) integer.
  
________________________________________________________________________

INTEGER FUNCTION MDGKDT(MADFIL) {
  INTEGER MADFIL

  Gets Kind of data (Kindat) integer.
  MDGKDT stands for MaDrec Get Kind of DaTa

  MADFIL - File handle

  Returns: Kind of data (Kindat) integer.
________________________________________________________________________

INTEGER FUNCTION MDIBYR(MADFIL) {
  INTEGER MADFIL

  Gets beginning year integer.

  MADFIL - File handle

  Returns: beginning year integer.
________________________________________________________________________

INTEGER FUNCTION MDIBDT(MADFIL) {
  INTEGER MADFIL

  Gets beginning date (MMDD) integer.

  MADFIL - File handle

  Returns: beginning date (MMDD) integer.
________________________________________________________________________

INTEGER FUNCTION MDIBHM(MADFIL) {
  INTEGER MADFIL

  Gets beginning hour/minute IBHM (HHMM) integer.

  MADFIL - File handle

  Returns: beginning hour/minute (HHMM) integer.
________________________________________________________________________


INTEGER FUNCTION MDIBCS(MADFIL) {
  INTEGER MADFIL

  Gets beginning second/centisecond IBCS (SSCC) integer.

  MADFIL - File handle

  Returns: beginning second/centisecond (SSCC) integer.
________________________________________________________________________


INTEGER FUNCTION MDIEYR(MADFIL) {
  INTEGER MADFIL

  Gets ending year integer.

  MADFIL - File handle

  Returns: ending year integer.
________________________________________________________________________

INTEGER FUNCTION MDIEDT(MADFIL) {
  INTEGER MADFIL

  Gets ending date (MMDD) integer.

  MADFIL - File handle

  Returns: ending date (MMDD) integer.
________________________________________________________________________

INTEGER FUNCTION MDIEHM(MADFIL) {
  INTEGER MADFIL

  Gets ending hour/minute IEHM (HHMM) integer.

  MADFIL - File handle

  Returns: ending hour/minute (HHMM) integer.
________________________________________________________________________


INTEGER FUNCTION MDIECS(MADFIL) {
  INTEGER MADFIL

  Gets ending second/centisecond IESC (SSCC) integer.

  MADFIL - File handle

  Returns: ending second/centisecond (SSCC) integer.
________________________________________________________________________


SUBROUTINE MDG2DP(MADFIL, PARCOD, PARM)
  INTEGER MADFIL, PARCOD
  DOUBLE PRECISION PARM(NRANMX)

  Gets a 2d parameter from the current record.

  MADFIL - File handle
  PARCOD - Parameter code
  PARM    - Parameter value array
________________________________________________________________________

SUBROUTINE MDGFLT(MADFIL, PARCOD, PARM)
  INTEGER MADFIL, PARCOD
  DOUBLE PRECISION PARM(NRANMX)

  Gets a flattened parameter from the current record. If its
  a 1D parameter, its value is copied into the first NROW
  doubles in PARM array.  If its a 2D parameter, it acts just
  like MDG2DP.  With this method all parameters act like 2D
  parameters.
  
  Stands for MaDrigal Get FLaTtened parameter

  MADFIL - File handle
  PARCOD - Parameter code (can be 1D or 2D)
  PARM    - Parameter value array
________________________________________________________________________

SUBROUTINE MDS1DP(MADFIL, PARCOD, PARM, INDEX)
  INTEGER MADFIL, PARCOD, INDEX
  DOUBLE PRECISION PARM

  Sets 1d parameter

  MADFIL - File handle
  PARCOD - Parameter code
  PARM    - Parameter value
  INDEX   - Position of parameter in 1d array (1<=INDEX<=JPAR)
  
  Note: to set special value missing, use PARM=1.e-38.
      To set special value assumed (for error parms only), use PARM=2.e-38
      To set special value knownbad (for error parms only), use PARM=3.e-38
________________________________________________________________________

SUBROUTINE MDS2DP(MADFIL, PARCOD, PARM, INDEX)
  INTEGER MADFIL

  Sets 2d parameter array

  MADFIL - File handle
  PARCOD - Parameter code
  PARM    - Parameter value
  INDEX   - Position of parameter in 2d array (1<=INDEX<=MPAR)
  
    Note: to set special value missing, use PARM=1.e-38.
      To set special value assumed (for error parms only), use PARM=2.e-38
      To set special value knownbad (for error parms only), use PARM=3.e-38
________________________________________________________

DOUBLE PRECISION FUNCTION MDSSFA(PARCOD)
  INTEGER PARCOD

  Gets parameter scale factor

  PARCOD - Parameter code

  Returns: Parameter scale factor
________________________________________________________________________

SUBROUTINE MDGERR(ERROR)
  CHARACTER ERROR(128)

  Gets most recent error message

  ERROR  - Error message
________________________________________________________________________

DOUBLE PRECISION FUNCTION GETKEY(YEAR,MON,DAY,HOUR,MIN,SEC)
  INTEGER YEAR,MON,DAY,HOUR,MIN,SEC

  Gets key (number of seconds) since YEAR,MON,DAY,HOUR,MIN,SEC
  
  
________________________________________________________________________


SUBROUTINE MDGEOD(MADFIL, ROW, GDLAT, GLON, GDALT)
  INTEGER MADFIL
  DOUBLE PRECISION GDLAT(NRANMX)
  DOUBLE PRECISION GLON(NRANMX)
  DOUBLE PRECISION GDALT(NRANMX)  

  Gets GDLAT, GLON, GDALT from current record via cedarGetGeodetic.

  MADFIL - File handle, pointing to current record
  GDLAT - geodetic latitude array
  GLON - longitude array
  GDALT - geodetic altitude array
________________________________________________________________________

SUBROUTINE GTROOT(STROOT,LENGTH)
  CHARACTER STROOT(128)
  INTEGER LENGTH

  Gets MADROOT as set either in env variable or in cedar.h

  STROOT  - String holding MADROOT
  LENGTH - length of string copied
________________________________________________________________________


*****************************************************************************/

Example Fortran programs using maddataF77

The basic premise of the MaddataF77 module is to hide the difference between measured and derived parameters when working with Madrigal data. There are two ways to input measured data with the MaddataF77 module, with a file and without a file. In the first case, measured data comes from a particular file, and in the second the application supplies it directly through input parameters.

This example page includes simple examples of using measured data from a file, and using measured data supplied by the application.

Using the MaddataF77 module with a file

The following example prints both measured and derived data from a particular file, with various filters applied using a string very similar to the isprint command line. The key method is CRMADD.

C
C     This sample program prints out all requested parameters
C     from a madrigal file for given filters.  The requested
C     parameters can be either measured or derived - the API
C     hides these details from the user.
C
C     Note that RFMADD must be declared INTEGER or INTEGER*8, 
C     depending on whether 32-bit or 64-bit platform

C     
      PROGRAM TMADDATA
C
C     .. Local Scalars ..
      INTEGER RFMADD
C     ..If 64-bit, should be INTEGER*8 RFMADD
      INTEGER NUMREC
C     .. External Subroutines ..
      EXTERNAL CRMADD,GTNROW, GTMADD, FRMADD, STALOC
C     ..
C     .. Local Arrays ..
      CHARACTER*128 FILEC, MDROOT
      CHARACTER PARMS*100
      CHARACTER FLTSTR*1000
C     .. Local variables ..
      INTEGER STATUS
      DOUBLE PRECISION DATROW(100)
      DOUBLE PRECISION SLATGD,SLON,SALTGD,SLATGC,SR
      INTEGER REC, NROW, ROW, ISTAT, KINST
C     ..
C     ..
      RFMADD = 0
      NUMREC = 0
C     the requested madrigal file
      CALL GTROOT(MDROOT,LEN)
      FILEC=MDROOT(:LEN)//'/experiments/1998/mlh/20jan98/mil980120g.002'
C     the requested parms (measured or derived)
      PARMS = 'gdalt,ti,kp,nel'
C     the desired filters (exactly like isprint command line)
      FLTSTR = 'z=1000, filter=range,2500, filter=ti,1000,2000'
      CALL CRMADD(FILEC,PARMS,FLTSTR,RFMADD,NUMREC,STATUS)
      IF (STATUS .NE. 0) THEN
         PRINT*, "Create Maddata failed for ", FILEC
         STOP
      END IF
C     print all records (ignoring formatting)
      DO 30, REC = 1, NUMREC
         CALL GTNROW(RFMADD, REC, NROW) 
         PRINT*,PARMS
C        loop over all rows
         DO 20, ROW=1, NROW
            CALL GTMADD(RFMADD, REC, ROW, DATROW, STATUS)
C           Note that data is stored in datrow in the order of parms
            PRINT "(F10.5)",DATROW(1),DATROW(2),DATROW(3),DATROW(4)
20       CONTINUE
30    CONTINUE
C     free all data
      CALL FRMADD(RFMADD)
C
      KINST = 31
      CALL STALOC(KINST,SLATGD,SLON,SALTGD,SLATGC,SR,ISTAT)
      if (ISTAT .EQ. 0) THEN
         print*, 'Kinst 31:', SLATGD,SLON,SALTGD,SLATGC,SR,ISTAT
      else
         print *, 'now returned error code ', ISTAT, SLATGD
      END IF
      
      CALL STALOC(3331,SLATGD,SLON,SALTGD,SLATGC,SR,ISTAT)
      if (ISTAT .EQ. 0) THEN
         print*, 'Kinst 3331:', SLATGD,SLON,SALTGD,SLATGC,SR,ISTAT
      else
         print*, 'Kinst 3331 failed, as it should'
      END IF
      
      END
      

Using the MaddataF77 module without a file

The following example prints both measured and derived data based on user supplied data alone. The key method is CRNFMD.


C     $Id: dev_fortran.html 3304 2011-01-17 15:25:59Z brideout $
C
C     This sample program prints out all requested parameters
C     from a madrigal data passed in as arguments - no file needed.  The requested
C     parameters can be either measured or derived - the API
C     hides these details from the user.
C     
      PROGRAM TMADDATA_NOFILE
C
C     .. Local Scalars ..
      INTEGER RFMADD, KINST
      DOUBLE PRECISION UT
C     .. External Subroutines ..
      EXTERNAL CRNFMD,GTNROW,GTMADD,FRMADD
C     .. External Functions ..
      DOUBLE PRECISION GETKEY
C     ..
C     .. Local Arrays ..
      CHARACTER*300 ONED, TWOD
      CHARACTER PARMS*100
C     .. Local variables ..
      INTEGER STATUS
      DOUBLE PRECISION DROW(100)
      INTEGER NROW, ROW
C     ..
C     ..
      RFMADD = 0
      NUMREC = 0
      KINST=31
      NROW=0
      UT = GETKEY(1998,1,20,15,0,0)
C     the requested parms 
      PARMS = 'gdlat,glon,gdalt,bmag,kp,sdwht'
C     the one-D data (not used in this example, but can't hurt)
      ONED = 'pl=.001 sn=10'   
C     the two-D data
      TWOD="gdlat=45,45,50,50 glon=20,30,20,30 gdalt=500,500,500,500"    
      CALL CRNFMD(PARMS,UT,KINST,ONED,TWOD,RFMADD,NROW)
      IF (NROW .EQ. 0) THEN
         PRINT*, "Create non-file Maddata failed"
         STOP
      END IF
C     print the record (only 1 returned, so its always rec 1)
      PRINT*,PARMS
      CALL GTNROW(RFMADD, 1, NROW) 
C     loop over all rows
      DO 20, ROW=1, NROW
         CALL GTMADD(RFMADD, 1, ROW, DROW, STATUS)
C        Note that data is stored in DROW in the order of parms
         PRINT*,DROW(1),DROW(2),DROW(3),DROW(4),DROW(5),DROW(6)
20    CONTINUE
C     free all data
      CALL FRMADD(RFMADD)
C
      END
      

Example Fortran program using madrecF77

C
      PROGRAM TMADREC
C
C     This program is a simple example illustrating the use of the madrec
C     module from Fortran - see madrecF77.c. This module is appropriate for
C     dealing with Madrigal files - writing them or modifying them.  To deal
C     with Madrigal data at a higher level, where the differences between
C     derived and measured data can be ignored, use maddataF77.c methods.
C
C     This program contain 4 examples:
C        1. Reading an existing madrigal file in sequentially
C        2. Writing a new madrigal file sequentially
C        3. Appending to an existing madrigal file using in mem file
C        4. Searching and summarizing a file in memory
C
C     .. Local Scalars ..
      INTEGER STATUS,NUMREC,LEN
      CHARACTER*128 ERROR,FNAME,MDROOT
      DOUBLE PRECISION PARM, KEY
      INTEGER LPROL, JPAR, MPAR, NROW, KREC, KINST, KINDAT
      INTEGER YEAR1, MONTH1, DAY1, HOUR1, MIN1, SEC1, CSEC1
      INTEGER YEAR2, MONTH2, DAY2, HOUR2, MIN2, SEC2, CSEC2
      INTEGER MADFIL, MADFL1, MADFL2
C     ..
C     .. Local Arrays ..
      DOUBLE PRECISION PARM2D(500)
C     ..
C     .. External Functions ..
      INTEGER MDOPEN,MDCLOS
      INTEGER MDREAD,MDISDR
      INTEGER MDWRIT,REWIND
      INTEGER MDCOPY,GRCKEY
      DOUBLE PRECISION GETKEY,GTPMIN,GTPMAX
C     ..
C     .. External Subroutines ..
      EXTERNAL MDGERR,MDG1DP,MDG2DP,MDCREA
      EXTERNAL MDS1DP,MDS2DP,GTROOT
C
      NUMREC = 0      
      CALL GTROOT(MDROOT,LEN)
      FNAME=MDROOT(:LEN)//'/experiments/1998/mlh/20jan98/mil980120g.002'
      ERROR=''
C
C
C     .....Example 1 - read a madrigal file.....
C
C
      PRINT*, "\nExample 1 - read in a madrigal file"
      MADFIL = MDOPEN(1, FNAME)
      IF (MADFIL.LT.0) THEN
         CALL MDGERR(ERROR)
         PRINT*,ERROR
         STOP
      END IF
C     loop through all the records
10    STATUS = MDREAD(MADFIL)
      IF (STATUS.NE.0) THEN
          GOTO 100
      END IF
C     count data records, ignore header or catalog
      IF (MDISDR(MADFIL).EQ.1) THEN
          NUMREC = NUMREC + 1
      END IF
C     print some data for record 5
      IF (NUMREC.EQ.5) THEN
          PRINT*,"The following is data for the 5th data record:"
          CALL MDG1DP(MADFIL,402, PARM)
          PRINT*,"   Pulse length is "
          PRINT "(F10.5)",PARM
          CALL MDG2DP(MADFIL, 550, PARM2D)
          PRINT*,"   The first 4 Ti values are:"
          PRINT "(F10.5)",PARM2D(1),PARM2D(2),PARM2D(3),PARM2D(4)
      END IF
      GOTO 10
C
100   PRINT*,"The number of data records found in file:"
      PRINT"(I6)",NUMREC
C     ..Close the file..
      STATUS = MDCLOS(MADFIL)
C
C
C     ....Example 2: create a new madrigal file ....
C
C
      PRINT*,""
      PRINT*, "Example 2 - create new file tMadrecF77.out"
      MADFIL = MDOPEN(20, 'tMadrecF77.out')
      IF (MADFIL.LT.0) THEN
         CALL MDGERR(ERROR)
         PRINT*,ERROR
         STOP
      END IF
C     create a new record
      LPROL = 16
      JPAR  = 2
      MPAR  = 1
      NROW  = 3
      KREC  = 1002
      KINST = 32
      KINDAT = 3408
      YEAR1 = 2003
      MONTH1 = 3
      DAY1   = 19
      HOUR1  = 1
      MIN1   = 0
      SEC1   = 0
      CSEC1  = 0
      YEAR2 = 2003
      MONTH2 = 3
      DAY2   = 19
      HOUR2  = 1
      MIN2   = 2
      SEC2   = 59
      CSEC2  = 99
      CALL MDCREA(MADFIL,
     *      LPROL, JPAR, MPAR, NROW,
     *      KREC, KINST, KINDAT,
     *      YEAR1, MONTH1, DAY1,
     *      HOUR1, MIN1, SEC1, CSEC1,
     *      YEAR2, MONTH2, DAY2,
     *      HOUR2, MIN2, SEC2, CSEC2)
C     set the two 1D values (pl and systmp)
      PARM = 1.28000e-03
      CALL MDS1DP(MADFIL, 402, PARM, 1)
      PARM = 151.0
      CALL MDS1DP(MADFIL, 482, PARM, 2)
C     set the 1 2D parm (range) with 3 rows
      PARM2D(1)=100.0
      PARM2D(2)=150.0
      PARM2D(3)=200.0
      CALL MDS2DP(MADFIL, 120, PARM2D, 1)
C     write the new record to file
      STATUS = MDWRIT(MADFIL)
      IF (STATUS.NE.0) THEN
         CALL MDGERR(ERROR)
         PRINT*,ERROR
         STOP
      END IF
C     ..Close the file..
      STATUS = MDCLOS(MADFIL)
C
C
C     ....Example 3: Append to an existing madrigal file ....
C
C      
      PRINT*,""
      PRINT*, "Example 3 - append to existing file"
      PRINT*, "  and save as tMadrecF77_append.out"
C     read the old file into memory
      MADFL1 = MDOPEN(50, FNAME)
      IF (MADFL1.LT.0) THEN
         CALL MDGERR(ERROR)
         PRINT*,ERROR
         STOP
      END IF
C     create a new record to append to it
      CALL MDCREA(MADFL1,
     *      LPROL, JPAR, MPAR, NROW,
     *      KREC, KINST, KINDAT,
     *      YEAR1, MONTH1, DAY1,
     *      HOUR1, MIN1, SEC1, CSEC1,
     *      YEAR2, MONTH2, DAY2,
     *      HOUR2, MIN2, SEC2, CSEC2)
C     set the two 1D values (pl and systmp)
      PARM = 1.28000e-03
      CALL MDS1DP(MADFL1, 402, PARM, 1)
      PARM = 151.0
      CALL MDS1DP(MADFL1, 482, PARM, 2)
C     set the 1 2D parm (range) with 3 rows
      PARM2D(1)=100.0
      PARM2D(2)=150.0
      PARM2D(3)=200.0
      CALL MDS2DP(MADFL1, 120, PARM2D, 1)
C     append the new record to the in-memory file
      STATUS = MDWRIT(MADFL1)
      IF (STATUS.NE.0) THEN
         CALL MDGERR(ERROR)
         PRINT*,ERROR
         STOP
      END IF
C     next rewind the in-memory file
      STATUS = REWIND(MADFL1)    
      IF (STATUS.NE.0) THEN
         CALL MDGERR(ERROR)
         PRINT*,ERROR
         STOP
      END IF
C     now we open another file to write the appended file to
      MADFL2 = MDOPEN(20, 'tMadrecF77_append.out')
      IF (MADFL2.LT.0) THEN
         CALL MDGERR(ERROR)
         PRINT*,ERROR
         STOP
      END IF
C     loop through all the in memory records, copy to MADFL2
20    STATUS = MDREAD(MADFL1)
      IF (STATUS.NE.0) THEN
          GOTO 200
      END IF
      STATUS = MDCOPY(MADFL1, MADFL2)
C     write the copied record to file
      STATUS = MDWRIT(MADFL2)      
      GOTO 20  
C    
200   STATUS = MDCLOS(MADFL1)     
      STATUS = MDCLOS(MADFL2)
C
C
C     ....Example 4: Manipulate an in memory file ....
C
C      
      PRINT*,""
      PRINT*, "Example 4: Manipulate an in memory file"
C     read the file into memory, this time with summary info
      MADFIL = MDOPEN(30, FNAME)
      IF (MADFIL.LT.0) THEN
         CALL MDGERR(ERROR)
         PRINT*,ERROR
         STOP
      END IF
      KEY = GETKEY(1998, 1, 20, 15, 0, 0)
      STATUS = GRCKEY(MADFIL,KEY)
      IF (STATUS.NE.0) THEN
         CALL MDGERR(ERROR)
         PRINT*,ERROR
         STOP
      END IF
C     print some data from this record
      PRINT*,"The following are min and max values of Ti in file:"
      PRINT "(F10.5)",GTPMIN(MADFIL,550),GTPMAX(MADFIL,550)
      PRINT*,"The following is data for key 1/20/1998 15:00:"
      CALL MDG1DP(MADFIL,402, PARM)
      PRINT*,"   Pulse length is "
      PRINT "(F10.5)",PARM
      CALL MDG2DP(MADFIL, 550, PARM2D)
      PRINT*,"   The first 4 Ti values are:"
      PRINT "(F10.5)",PARM2D(1),PARM2D(2),PARM2D(3),PARM2D(4)
      STATUS = MDCLOS(MADFIL)
C
      PRINT*,"\nTest complete" 
C
      END

GEOLIB Procedure Synopsis

The geolib library (geolib.a) contains the following procedures focused on time and space.


***********************************************************************

      SUBROUTINE CARMEL(B,XI,VL)
C
C     Private/Internal subroutine. Part of Apex coordinate computation
C     package. See COORD for public API. Computes scaler VL as a
C     function of B and XI.
C
C       Input:
C           B - Scaler field strength value.
C          XI - integral invariant (see INTEG).
C
C      Output:
C          VL - McIlwain's L-shell parameter i.e. Invariant Latitude
C               = ACOS(DSQRT(1.0D0/VL))/DTR
C
C     .. Scalar Arguments ..
      DOUBLE PRECISION B,VL,XI
C     ..



***********************************************************************

      SUBROUTINE CONVRT(I,GDLAT,GDALT,GCLAT,RKM)
C
C     jmh - 11/79  ans fortran 66
C
C     Converts between geodetic and geocentric coordinates. the
C     reference geoid is that adopted by the iau in 1964. a=6378.16,
C     b=6356.7746, f=1/298.25. the equations for conversion from
C     geocentric to geodetic are from astron. j., vol 66, 1961, p. 15.
C
C       Input:
C             I - 1 geodetic to geocentric, 2 geocentric to geodetic
C       Input, Output:
C         GDLAT - geodetic latitude (degrees)
C         GDALT - altitude above geoid (km)
C         GCLAT - geocentric latitude (degrees)
C           RKM - geocentric radial distance (km)
C
C     .. Scalar Arguments ..
      DOUBLE PRECISION GCLAT,GDALT,GDLAT,RKM
      INTEGER I
C     ..



***********************************************************************

      SUBROUTINE COORD(SLATGD,SLON,SR,SLATGC,TM,AZ,EL,RANGE,GDLAT,GLON,
     *                 GDALT,RCOR)
C
C     Calculates the listed coordinates of a specified point. the
C     point may be specified either by slatgd, slon, sr, slatgc, tm,
C     az, el, range (range .ge. 0.) or by gdlat, glon, gdalt (range
C     .lt. 0.)
C
C     Input:
C        SLATGD - station geodetic latitude
C        SLON   - station longitude
C        SR     - radial distance of station from center of earth
C        SLATGC - station geocentric latitude
C        TM     - time in years (e.g. 1975.2)
C        AZ     - radar azimuth
C        EL     - radar elevation
C        RANGE  - range to observation point
C     Input, output:
C        GDLAT  - geodetic latitude of observation point
C        GLON   - longitude of observation point
C        GDALT  - altitude above spheroid of observation point
C     Output:
C        RCOR(7)  - b     - magnitude of geomagnetic field
C        RCOR(8)  - br    - radial component of geomagnetic field
C        RCOR(9)  - bt    - southward component of geomagnetic field
C        RCOR(10) - bp    - eastward component of geomagnetic field
C        RCOR(11) - rlatm - dip latitude
C        RCOR(12) - rlati - invariant latitude
C        RCOR(13) - rl    - magnetic l parameter
C        RCOR(14) - alat  - apex latitude
C        RCOR(15) - alon  - apex longitude
C
C        RCOR(16) - g(1,1) magnetic coordinate system metric tensor,
C                          upper half stored row-wise
C        RCOR(17) - g(2,1) "                                       "
C        RCOR(18) - g(2,1) "                                       "
C        RCOR(19) - g(2,1) "                                       "
C
C        RCOR(20) - south-direction cosine w/respect to geodetic coords.
C        RCOR(21) - east-direction cosine "                            "
C        RCOR(22) - upward-direction cosine "                          "
C
C        RCOR(23) - perpendicular to b in magnetic n - s plane
C                   (magnetic south)
C        RCOR(24) - perpendicular to b in horizontal plane
C                   (magnetic east)
C        RCOR(25) - upward along magnetic field
C
C        RCOR(26) - x-direction cosine of a vector perpendicular to
C                   l.o.s. w/respect to apex coords.
C        RCOR(27) - y-direction cosine "                          "
C        RCOR(28) - z-direction cosine "                          "
C
C        RCOR(29) - inclination of geomagnetic field
C        RCOR(30) - declination of geomagnetic field
C        RCOR(31) - gclat - geocentric latitude
C        RCOR(32) - aspct - aspect angle
C        RCOR(33) - conjugate geocentric latitude
C        RCOR(34) - conjugate geodetic latitude
C        RCOR(35) - conjugate longitude
C
C     .. Scalar Arguments ..
      DOUBLE PRECISION AZ,EL,GDALT,GDLAT,GLON,RANGE,SLATGC,SLATGD,SLON,
     *                 SR,TM
C     ..
C     .. Array Arguments ..
      DOUBLE PRECISION RCOR(38)
C     ..



***********************************************************************

      SUBROUTINE CSCONV(X,Y,Z,R,THETA,PHI,IMODE)
C
C     jmh - 11/79  ans fortran 66
C
C     Converts between cartesian coordinates x,y,z and spherical
C     coordinates r,theta,phi.  theta and phi are in degrees.
C
C       Input:
C         IMODE - 1 (x,y,z) -> (r,theta,phi)
C                 2 (r,theta,phi) -> (x,y,z)
C
C       Input, Output:
C              X, Y, Z - cartesian coordinates
C        R, THETA, PHI - spherical coordinates (degrees)
C
C     .. Scalar Arguments ..
      DOUBLE PRECISION PHI,R,THETA,X,Y,Z
      INTEGER IMODE
C     ..



***********************************************************************

      SUBROUTINE DIPLAT(TM,GDLAT,GLON,GDALT,AINC,DEC,RLATM)
C     .. Scalar Arguments ..
      DOUBLE PRECISION AINC,DEC,GDALT,GDLAT,GLON,RLATM,TM
C     ..



***********************************************************************

      SUBROUTINE DSF(GCLAT,GLON,RKM,ALT,HALT,ISTOP,DS)
C
C     Calculates an optimum integration step size for geomagnetic
C     field line tracing routine LINTRA, as an empirical function of
C     the geomagnetic dipole coordinates of the starting point. when
C     start and end points are in the same hemisphere and
C     abs(halt-alt)<10000, the empirical formula is modified so that
C     ds is no greater than 1/100 of the altitude difference between
C     start and end points.
C
C     input:
C         GCLAT - start point geocentric latitude (deg)
C          GLON - start point geocentric longitude (deg)
C           RKM - start point radial distance (km)
C           ALT - starting point geodetic altitude
C          HALT - tuning parameter for altitude test. 
C         ISTOP - set to -1 to enable abs(halt-alt)<10000 test
C     output:
C            DS - step size
C
C     .. Scalar Arguments ..
      DOUBLE PRECISION ALT,DS,GCLAT,GLON,HALT,RKM
      INTEGER ISTOP
C     ..



***********************************************************************

      SUBROUTINE GASPCT(SLATGD,SLON,SR,SLATGC,TM,AZ,EL,RANGE,GDLAT,GLON,
     *                  GDALT,B,CASPCT,ASPCT)
C
C     jmh - 3/88
C
C     *** warning *** this routing used to be called aspect. the name
C                     was changed on 3/2/88 to avoid conflict with a
C                     new routine of the same name in the geophysics
C                     library
C
C     Calculates the aspect angle between a radar beam and the
C     geomagnetic field at a specified point. the point may be
C     specified either by slatgd, slon, sr, slatgc, tm, az, el, range
C     (range .ge. 0.) or by gdlat, glon, gdalt (range .lt. 0.)
C
C     Input:
C        SLATGD - station geodetic latitude
C        SLON   - station longitude
C        SR     - radial distance of station from center of earth
C        SLATGC - station geocentric latitude
C        TM     - time in years (e.g. 1975.2)
C        AZ     - radar azimuth
C        EL     - radar elevation
C        RANGE  - range to observation point
C     Input, output:
C         GDLAT - geodetic latitude of observation point
C         GLON  - longitude of observation point
C         GDALT - altitude above spheroid of observation point
C     Output:
C        B      - geomagnetic field magnitude
C        CASPCT - cosine of aspect angle
C        ASPCT  - aspect angle (degrees)
C
C
C
C     .....subroutine parameter specifications.....
C
C     .....local specifications.....
C
C     .. Scalar Arguments ..
      DOUBLE PRECISION ASPCT,AZ,B,CASPCT,EL,GDALT,GDLAT,GLON,RANGE,
     *                 SLATGC,SLATGD,SLON,SR,TM
C     ..



***********************************************************************

      SUBROUTINE GDMAG(TM,GDLAT,GLON,GDALT,X,Y,Z,F,H,DEC,AINC)
C
C     jmh - 1/80  ans fortran 66
C
C     Evaluates the geomagnetic field at a point specified by its
C     geodetic coordinates. the reference geoid is that adopted by the
C     iau in 1964.
C
C     input:
C           TM - time in years for desired field (e.g. 1971.25)
C        GDLAT - geodetic latitude (degrees)
C         GLON - east longitude (degrees)
C        GDALT - altitude above geoid (km)
C
C     output:
C        X,Y,Z - geodetic field components (gauss)
C            F - magnitude of field (gauss)
C            H - horizontal intensity (gauss)
C          DEC - declination (degrees)
C         AINC - inclination (degrees)
C
C     .. Scalar Arguments ..
      DOUBLE PRECISION AINC,DEC,F,GDALT,GDLAT,GLON,H,TM,X,Y,Z
C     ..



***********************************************************************

      DOUBLE PRECISION FUNCTION GDRAN(SR,SLATGC,SLON,AZ,EL,ALT)
C
C     GDRAN uses a half-interval (binary) search technique to determine
C     the geodetic range to a point of observation given in terms of
C     azimuth, elevation, and geodetic altitude.
C
C       Input:
C            SR - radial distance of station from center of earth
C        SLATGC - station geocentric latitude
C          SLON - station longitude
C            AZ - radar azimuth
C            EL - radar elevation
C           ALT - geodetic altitude
C
C     harris fortran 77
C     rgm - 8/85
C
C     .. Scalar Arguments ..
      DOUBLE PRECISION ALT,AZ,EL,SLATGC,SLON,SR
C     ..



***********************************************************************

      SUBROUTINE GDV(GDLAT,GCLAT,FR,FT,FP,FX,FY,FZ)
C
C     jmh - 11/79  ans fortran 66
C
C     GDV converts a vector field f at geodetic latitude GDLAT and
C     geocentric latitude GCLAT from a geocentric based representation
C     to a geodetic based representation. the geocentric components
C     are FR (radial outward), FT (increasing geocentric colatitude,
C     e.g. southward) and FP (increasing east longitude). the
C     geodetic components are FX (northward, parallel to surface of
C     earth), FY (eastward, parallel to surface of earth) and FZ
C     (downward, perpendicular to surface of earth). FR,FT,FP thus
C     correspond to spherical coordinates r,theta,phi, with their
C     origin at the center of the earth. x,y,z are the coordinates
C     customarily used to describe the three components of the
C     geomagnetic field. FP and FY are the same.
C
C       Input:
C         GDLAT - geodetic latitude (degrees)
C         GCLAT - geocentric latitude (degrees)
C            FR - radial outward (geocentric).
C            FT - increasing geocentric colatitude (southward).
C            FP - increasing east longitude.
C
C      Output:
C            FX - northward, parallel to surface of earth (geodetic).
C            FY - eastward, parallel to surface of earth.
C            FZ - downward, perpendicular to surface of earth.
C
C     .. Scalar Arguments ..
      DOUBLE PRECISION FP,FR,FT,FX,FY,FZ,GCLAT,GDLAT
C     ..



***********************************************************************

      SUBROUTINE GMET(DXDQ,G)
C
C     jmh - 10/79  ans fortran 66
C
C     GMET calculates the metric tensor G of a coordinate system q
C     for which dx(i)/dq(j)=dxdq(i,j).
C
C       Input:
C         DXDQ - coordinate system array.
C
C      Output:
C            G - metric tensor.
C
C     .. Array Arguments ..
      DOUBLE PRECISION DXDQ(3,3),G(3,3)
C     ..



***********************************************************************

      SUBROUTINE GTS5(IYD,SEC,ALT,GLAT,GLONG,STL,F107A,F107,AP,MASS,D,T)
C        MSIS-86/CIRA 1986 Neutral Thermosphere Model
C         A.E.Hedin 3/15/85;2/26/87 (Variable Names Shortened)
C         9/21/87 M.E. Hagan (Non-Standard Statements Changed)
C     INPUT:
C        IYD - YEAR AND DAY AS YYDDD
C        SEC - UT(SEC)
C        ALT - ALTITUDE(KM) (GREATER THAN 85 KM)
C        GLAT - GEODETIC LATITUDE(DEG)
C        GLONG - GEODETIC LONGITUDE(DEG)
C        STL - LOCAL APPARENT SOLAR TIME(HRS)
C        F107A - 3 MONTH AVERAGE OF F10.7 FLUX
C        F107 - DAILY F10.7 FLUX FOR PREVIOUS DAY
C             UNITS:1.0e-22W/m2/Hz
C        AP - MAGNETIC INDEX(DAILY) OR WHEN SW(9)=-1. :
C           - ARRAY CONTAINING:
C             (1) DAILY AP
C             (2) 3 HR AP INDEX FOR CURRENT TIME
C             (3) 3 HR AP INDEX FOR 3 HRS BEFORE CURRENT TIME
C             (4) 3 HR AP INDEX FOR 6 HRS BEFORE CURRENT TIME
C             (5) 3 HR AP INDEX FOR 9 HRS BEFORE CURRENT TIME
C             (6) AVERAGE OF EIGHT 3 HR AP INDICIES FROM 12 TO 33 HRS
C                    PRIOR TO CURRENT TIME
C             (7) AVERAGE OF EIGHT 3 HR AP INDICIES FROM 36 TO 59 HRS
C                    PRIOR TO CURRENT TIME
C        MASS - MASS NUMBER (ONLY DENSITY FOR SELECTED GAS IS
C                 CALCULATED. MASS 0 IS TEMPERATURE. MASS 48 FOR ALL.
C
C     OUTPUT:
C        D(1) - HE NUMBER DENSITY(CM-3)
C        D(2) - O NUMBER DENSITY(CM-3)
C        D(3) - N2 NUMBER DENSITY(CM-3)
C        D(4) - O2 NUMBER DENSITY(CM-3)
C        D(5) - AR NUMBER DENSITY(CM-3)
C        D(6) - TOTAL MASS DENSITY(GM/CM3)
C        D(7) - H NUMBER DENSITY(CM-3)
C        D(8) - N NUMBER DENSITY(CM-3)
C        T(1) - EXOSPHERIC TEMPERATURE
C        T(2) - TEMPERATURE AT ALT
C
C      TO GET OUTPUT IN M-3 and KG/M3:   CALL METERS(.TRUE.)
C
C          ADDITIONAL COMMENTS
C           (1) LOWER BOUND QUANTITIES IN COMMON/GTS3C/
C           (2) TO TURN ON AND OFF PARTICULAR VARIATIONS CALL TSELEC(SW)
C               WHERE SW IS A 25 ELEMENT ARRAY CONTAINING 0. FOR OFF, 1.
C               FOR ON, OR 2. FOR MAIN EFFECTS OFF BUT CROSS TERMS ON
C               FOR THE FOLLOWING VARIATIONS
C               1 - F10.7 EFFECT ON MEAN  2 - TIME INDEPENDENT
C               3 - SYMMETRICAL ANNUAL    4 - SYMMETRICAL SEMIANNUAL
C               5 - ASYMMETRICAL ANNUAL   6 - ASYMMETRICAL SEMIANNUAL
C               7 - DIURNAL               8 - SEMIDIURNAL
C               9 - DAILY AP             10 - ALL UT/LONG EFFECTS
C              11 - LONGITUDINAL         12 - UT AND MIXED UT/LONG
C              13 - MIXED AP/UT/LONG     14 - TERDIURNAL
C              15 - DEPARTURES FROM DIFFUSIVE EQUILIBRIUM
C              16 - ALL TINF VAR         17 - ALL TLB VAR
C              18 - ALL T0 VAR           19 - ALL S VAR
C              20 - ALL Z0 VAR           21 - ALL NLB VAR
C              22 - ALL TR12 VAR         23 - TURBO SCALE HEIGHT VAR
C
C              To get current values of SW: CALL TRETRV(SW)
C
C     .. Scalar Arguments ..
      DOUBLE PRECISION ALT,F107,F107A,GLAT,GLONG,SEC,STL
      INTEGER IYD,MASS
      LOGICAL METER
C     ..
C     .. Array Arguments ..
      DOUBLE PRECISION AP(*),D(8),T(2)
C     ..
C     .. Scalars in Common ..
C     ..
C     .. Arrays in Common ..
C     ..



***********************************************************************

      DOUBLE PRECISION FUNCTION DENSS(ALT,DLB,TINF,TLB,XM,ALPHA,TZ,ZLB,
     *                                S2,T0,ZA,Z0,TR12)
C       Calculate Temperature and Density Profiles for MSIS models
C     .. Scalar Arguments ..
      DOUBLE PRECISION ALPHA,ALT,DLB,S2,T0,TINF,TLB,TR12,TZ,XM,Z0,ZA,ZLB
C     ..
C     .. Scalars in Common ..
C     ..



***********************************************************************

      DOUBLE PRECISION FUNCTION GLOBE5(YRD,SEC,LAT,LONG,TLOC,F107A,F107,
     *                 AP,P)
C       CALCULATE G(L) FUNCTION FOR MSIS-86/CIRA 1986
C       Upper Thermosphere Parameters
C     .. Scalar Arguments ..
      DOUBLE PRECISION F107,F107A,LAT,LONG,SEC,TLOC,YRD
C     ..
C     .. Array Arguments ..
      DOUBLE PRECISION AP(*),P(*)
C     ..
C     .. Scalars in Common ..
C     ..
C     .. Arrays in Common ..
C     ..



***********************************************************************

      SUBROUTINE TSELEC(SV)
C        SET SWITCHES
C     .. Array Arguments ..
      DOUBLE PRECISION SV(*)
C     ..
C     .. Scalars in Common ..
C     ..
C     .. Arrays in Common ..
C     ..



***********************************************************************

      DOUBLE PRECISION FUNCTION GLOB5L(P)
C      LIMITED PARAMETER VERSION OF GLOBE 9/2/82
C       CALCULATE G(L) FUNCTION FOR MSIS-86/CIRA 1986
C       Lower Thermosphere Parameters
C     .. Array Arguments ..
      DOUBLE PRECISION P(*)
C     ..
C     .. Scalars in Common ..
C     ..
C     .. Arrays in Common ..
C     ..



***********************************************************************

      DOUBLE PRECISION FUNCTION DNET(DD,DM,ZHM,XMM,XM)
C       8/20/80
C       TURBOPAUSE CORRECTION FOR MSIS MODELS
C       Eq. A12b
C     .. Scalar Arguments ..
      DOUBLE PRECISION DD,DM,XM,XMM,ZHM
C     ..



***********************************************************************

      DOUBLE PRECISION FUNCTION CCOR(ALT,R,H1,ZH)
C        CHEMISTRY/DISSOCIATION CORRECTION FOR MSIS MODELS
C     Eq. A20a or Eq. A21
C     .. Scalar Arguments ..
      DOUBLE PRECISION ALT,H1,R,ZH
C     ..



***********************************************************************

      SUBROUTINE PRMSG5
C          CIRA     11-FEB-86
C     .. Arrays in Common ..
C     ..



***********************************************************************

      DOUBLE PRECISION FUNCTION HFUN(SR,EL,RANGE)
C
C     jmh - 11/79  ans fortran 66
C
C     HFUN computes the height above a sphere (radius SR) of an
C     observation point at a specified elevation (EL) and range
C     (RANGE). SR and RANGE should be positive and EL should be in the
C     range 0.0 to 90.0.
C
C       Input:
C           SR - radial distance of station from center of earth
C           EL - radar elevation
C        RANGE - range to observation point
C
C     .. Scalar Arguments ..
      DOUBLE PRECISION EL,RANGE,SR
C     ..
C     .. Scalars in Common ..
C     ..



***********************************************************************

      SUBROUTINE INTEG(ARC,BEG,BEND,B,JEP,ECO,FI)
C
C     Private/Internal subroutine. Part of Apex coordinate computation
C     package. See COORD for public API. INTEG determines the value of
C     the integral invariant FI by numerically integrating along the
C     field line from the specified point of interest to its conjugate
C     point.
C
C       Input:
C         ARC - Altitudes in earth radii (array).
C         BEG - floating point array.
C        BEND - floating point array.
C           B - Magnitude of field (array)
C         ECO - floating point array.
C         JEP - floating point scaler.
C
C      Output:
C          FI - floating point scaler.
C
C     .. Scalar Arguments ..
      DOUBLE PRECISION FI
      INTEGER JEP
C     ..
C     .. Array Arguments ..
      DOUBLE PRECISION ARC(200),B(200),BEG(200),BEND(200),ECO(200)
C     ..



***********************************************************************

      SUBROUTINE INVAR(TM,FLAT,FLONG,ALT,ERR,BB,FL)
C
C     Private/Internal subroutine. Part of Apex coordinate computation
C     package. See COORD for public API. INVAR converts coordinates
C     TM, FLAT, FLON and ALT to L-shell coordinates FL and BB. The
C     uncertainty in FL is typically less than 10.*ERR*FL (percent)
C
C       Input:
C           TM - time in years for desired field (e.g. 1971.25)
C         FLAT - geocentric latitude (degrees)
C        FLONG - east longitude
C          ALT - altitude (km)
C          ERR - tolerance factor
C
C      Output:
C           BB - Magnetic Field strength at point.
C           FL - McIlwain's L-shell parameter i.e. 
C                Invariant Latitude = ACOS(DSQRT(1.0D0/FL))/DTR
C
C     .. Scalar Arguments ..
      DOUBLE PRECISION ALT,BB,ERR,FL,FLAT,FLONG,TM
C     ..



***********************************************************************

      SUBROUTINE INVLAT(TM,GDLAT,GLON,GDALT,RL,RLATI)
C     .. Scalar Arguments ..
      DOUBLE PRECISION GDALT,GDLAT,GLON,RL,RLATI,TM
C     ..



***********************************************************************

      SUBROUTINE ITERAT
C
C     Private/Internal subroutine. Part of Apex coordinate computation
C     package. See COORD for public API. ITERAT integrates magnetic
C     field line using a 4-point adams formula after initialization.
C     First 7 iterations advance point by 3*DS.
C
C       Input (via common block ITER):
C               L - step count. set l=1 first time thru,
C                   set l=l+1 thereafter.
C      B,BR,BT,BP - field + components at point y
C              ST - sine of geocentric colatitude
C             SGN - sgn=+1  traces in direction of field
C                   sgn=-1  traces in negative field direction
C              DS - integration stepsize (arc increment) in km
C
C      Input, Output (via common block ITER):
C               Y - R,DLAT,DLON: geocentric tracing point coordinates 
C                   (km,deg)
C
C      Output (via common block ITER):
C            YOLD - Y at iteration L-1
C
C     .. Scalars in Common ..
C     ..
C     .. Arrays in Common ..
C     ..



***********************************************************************

      SUBROUTINE LINES(R1,R2,R3,B,ARC,ERR,J,VP,VN,TM)
C
C     Private/Internal subroutine. Part of Apex coordinate computation
C     package. See COORD for public API. Makes repeated calls to the
C     IGRF, tracing magnetic field line to minimum B.
C
C       Input:
C          ERR - tolerance factor (see INVAR)
C           TM - Time in floating point years (e.g. 1995.7)
C
C       Input, Output:
C     R1,R2,R3 - field strength in geocentric directions.
C            B - Magnitude of field (array)
C          ARC - Altitudes in earth radii (array)
C           VP - Geocentric latitude (array)
C           VN - Geocentric longitude (array)
C
C      Output:
C            J - Number of points in trace.
C
C     .. Scalar Arguments ..
      DOUBLE PRECISION ERR,TM
      INTEGER J
C     ..
C     .. Array Arguments ..
      DOUBLE PRECISION ARC(200),B(200),R1(3),R2(3),R3(3),VN(3),VP(3)
C     ..



***********************************************************************

      SUBROUTINE LINTRA(TM,GCLAT,GLON,RKM,ALT,HALT,PLAT,PLON,PRKM,ARC,
     *                  ARAD,ALAT,ALON,ISTOP,NPR,INIT,IER)
C
C     jmh - 11/79  ans fortran 66
C
C     LINTRA traces the geomagnetic field line from a specified
C     starting point to either a specified altitude in either hemisphere
C     or to the apex of the field line. in either case, if the apex
C     is passed, the apex coordinates of the field line are calculated.
C     when points are found on either side of the end point or apex,
C     the final result is calculated by quadratic interpolation.
C
C       Input:
C            TM - time for desired field (years)
C         GCLAT - start point geocentric latitude (deg)
C         GCLON - start point geocentric longitude (deg)
C           RKM - start point radial distance (km)
C           ALT - start point geodetic altitude (km)
C         ISTOP = -1 - trace to altitude halt in same hemisphere
C         ISTOP = 0 - trace to apex of field line
C         ISTOP = +1 - trace to altitude halt in opposite hemisphere....
C         NPR=1 - return to calling program after each step
C
C       Input, Output:
C          HALT - end point geodetic altitude (km)
C        INIT=1 - set by calling program when returning to lintra after
C                 receiving intermediate results
C             2 - set by lintra when trace is complete
C
C      Output:
C          PLAT - end point geocentric latitude (deg)
C          PLON - end point geocentric longitude (deg)
C          PRKM - end point radial distance (km)
C           ARC - arc length of field line traced (km)
C          ARAD - apex radius of field line (earth radii)
C          ALAT - apex latitude of field line (deg)
C          ALON - apex longitude of field line (deg)
C           IER = 1 - error, number of steps exceeds maxs
C
C     .. Scalar Arguments ..
      DOUBLE PRECISION ALAT,ALON,ALT,ARAD,ARC,GCLAT,GLON,HALT,PLAT,PLON,
     *                 PRKM,RKM,TM
      INTEGER IER,INIT,ISTOP,NPR
C     ..
C     .. Scalars in Common ..
C     ..



***********************************************************************

      SUBROUTINE LOOK(SR,SLAT,SLON,PR,GLAT,GLON,AZ,EL,RANGE)
C
C     jmh - 1/80  ans fortran 66
C
C     LOOK calculates the azimuth, elevation and range from a radar
C     of a specified point.
C
C       Input:
C        SR    - distance of station from center of earth (km)
C        SLAT  - geocentric latitude of station (deg)
C        SLON  - longitude of station (deg)
C        PR    - distance from center of earth of observation point (km)
C        GLAT  - observation point geocentric latitude (deg)
C        GLON  - observation point longitude (deg)
C
C      Output:
C        AZ    - radar azimuth (deg)
C        EL    - radar elevation (deg)
C        RANGE - radar range (km)
C
C     ...calculate "observation-point" earth centered cartesian coords..
C     .. Scalar Arguments ..
      DOUBLE PRECISION AZ,EL,GLAT,GLON,PR,RANGE,SLAT,SLON,SR
C     ..



***********************************************************************

      SUBROUTINE MILMAG(TM,RKM,ST,CT,SPH,CPH,BR,BT,BP,B)
C
C     MILMAG evaluates the geomagnetic field at a point specified by
C     its geocentric coordinates. 
C
C     Modified by B. Rideout - Dec. 26, 2002
C     This method is now simply a thin wrapper around geo-cgm code,
C     method igrf.  See geo-cgm.f for details

C
C       Input:
C              TM - time in years for desired field (e.g. 1971.25)
C             RKM - geocentric distance (km)
C           ST,CT - sin and cos of geocentric colatitude
C         SPH,CPH - sin and cos of east longitude
C
C      Output:
C        BR,BT,BP - geocentric field components (gauss)
C               B - magnitude of field (gauss)
C     ..
C     .. Scalar Arguments ..
      DOUBLE PRECISION B,BP,BR,BT,CPH,CT,RKM,SPH,ST,TM
C     ..



***********************************************************************

      SUBROUTINE MINV(A,N,D,L,M)
C
C     Inverts general matrix A (overwrites A) using standard
C     gauss-jordan method. The determinant is also calculated. a
C     determinant of zero indicates that the matrix is singular.
C
C      Input:
C        n - order of matrix a
C        l - work vector of length n
C        m - work vector of length n
C
C      Input, Output:
C        a - input matrix, destroyed in computation and replaced by
C            resultant inverse.
C
C      Output:
C        d - resultant determinant
C
C     .. Scalar Arguments ..
      DOUBLE PRECISION D
      INTEGER N
C     ..
C     .. Array Arguments ..
      DOUBLE PRECISION A(*)
      INTEGER L(*),M(*)
C     ..



***********************************************************************

      SUBROUTINE MTRAN3(A)
C
C     jmh - 1/29/80  ans fortran 66
C
C     MTRAN3 calculates the transpose of a 3 x 3 matrix a
C
C       Input, Output:
C           A - input matrix, replaced with transpose.
C
C     .. Array Arguments ..
      DOUBLE PRECISION A(3,3)
C     ..




***********************************************************************

      SUBROUTINE POINT(SR,SLAT,SLON,AZ,EL,RANGE,PR,GLAT,GLON)
C
C     jmh - 1/80  ans fortran 66
C
C     POINT calculates the position of a point defined by the radar
C     line-of sight vector to that point.
C
C     Input:
C       SR    - distance of station from center of earth (km)
C       SLAT  - geocentric latitude of station (deg)
C       SLON  - longitude of station (deg)
C       AZ    - radar azimuth (deg)
C       EL    - radar elevation (deg)
C       RANGE - radar range (km)
C
C     Output:
C       PR    - distance from center of earth of observation point (km)
C       GLAT  - observation point geocentric latitude (deg)
C       GLON  - observation point longitude (deg)
C
C     ...calculate "line-of-sight" station centered cartesian coords...
C     .. Scalar Arguments ..
      DOUBLE PRECISION AZ,EL,GLAT,GLON,PR,RANGE,SLAT,SLON,SR
C     ..



***********************************************************************

      DOUBLE PRECISION FUNCTION RFUN(SR,EL,H)
C
C     jmh - 11/79  ans fortran 66
C
C     RFUN computes the range to an observation point at a specified
C     elevation (EL) and distance (H) above a sphere of radius SR.
C     SR and H should be positive and EL should be in the range
C     0.0 to 90.0.
C
C       Input:
C         SR - radius of sphere (km)
C         EL - elevation (deg)
C          H - distance above sphere (km)
C
C     .. Scalar Arguments ..
      DOUBLE PRECISION EL,H,SR
C     ..



***********************************************************************

      SUBROUTINE RPCART(SR,SLAT,SLON,AZ,EL,RANGE,RFX,RFY,RFZ,PFX,PFY,
     *                  PFZ)
C
C     jmh - 11/79  ans fortran 66
C
C     RPCART computes the components (RFX,RFY,RFZ) relative to an earth
C     centered cartesian coordinate system of the radar line of sight
C     vector from a radar with coordinates SR (distance from center
C     of earth), SLAT (geocentric latitude) and SLON (longitude). the
C     observation point is specified by AZ (azimuth), EL (elevation) and
C     RANGE (range). the cartesian coordinates of the observation
C     point are returned in (PFX,PFY,PFZ).
C
C
C     Input:
C       SR    - distance of station from center of earth (km)
C       SLAT  - geocentric latitude of station (deg)
C       SLON  - longitude of station (deg)
C       AZ    - radar azimuth (deg)
C       EL    - radar elevation (deg)
C       RANGE - radar range (km)
C
C     Output:
C       RFX,RFY,RFZ - earth centered cartesian coordinate components
C                     of radar line of sight.
C       PFX,PFY,PFZ - earth centered cartesian coordinate components
C                     of observation point.
C
C     .. Scalar Arguments ..
      DOUBLE PRECISION AZ,EL,PFX,PFY,PFZ,RANGE,RFX,RFY,RFZ,SLAT,SLON,SR
C     ..



***********************************************************************

      DOUBLE PRECISION FUNCTION SPROD(A,B)
C
C     jmh - 11/79  ans fortran 66
C
C     SPROD calculates the scalar product of two vectors A and B.
C
C     Input:
C        A - floating point vector of dimension 3
C        B - floating point vector of dimension 3
C
C     .. Array Arguments ..
      DOUBLE PRECISION A(3),B(3)
C     ..
      SPROD = A(1)*B(1) + A(2)*B(2) + A(3)*B(3)
      RETURN
C
      END


***********************************************************************

      SUBROUTINE STARTR(R1,R2,R3,B,ARC,V,TM)
C
C     Private/Internal subroutine. Part of Apex coordinate computation
C     package. See COORD for public API.
C
C       Input:
C           TM - time in years for desired field (e.g. 1971.25)
C
C       Input, Output:
C          ARC - Altitudes in earth radii (array).
C            V - floating point array.
C
C      Output:
C     R1,R2,R3 - field strength in geocentric directions.
C            B - Magnitude of field (array)
C
C     .. Scalar Arguments ..
      DOUBLE PRECISION TM
C     ..
C     .. Array Arguments ..
      DOUBLE PRECISION ARC(200),B(200),R1(3),R2(3),R3(3),V(3,3)
C     ..



***********************************************************************

      SUBROUTINE UTHM(UTM,IUH,IUM)
C     Converts UTM from the hour and fraction to HH:MM
C     .. Scalar Arguments ..
      DOUBLE PRECISION UTM
      INTEGER IUH,IUM
C     ..
C     .. Intrinsic Functions ..
      INTRINSIC INT,NINT
C     ..
      IUH = INT(UTM)
      IF (IUH.EQ.99) THEN
         IUM = 99
      ELSE
         IUM = NINT((UTM-IUH)*60)
      END IF
      RETURN
      END
C  *********************************************************************
C  =====================================================================
      SUBROUTINE GEOCGM01(ICOR,IYEAR,HI,DAT,PLA,PLO)
C   Version 2001 for GEO-CGM.FOR                              April 2001
C     Apr 11, 2001  GEOLOW is modified to account for interpolation of
C                CGM meridians near equator across the 360/0 boundary
C     AUTHORS:
C     Natalia E. Papitashvili (WDC-B2, Moscow, Russia, now at NSSDC,
C     NASA/Goddard Space Flight Center, Greenbelt, Maryland)
C     Vladimir O. Papitashvili (IZMIRAN, Moscow, Russia, now at SPRL,
C     The University of Michigan, Ann Arbor)
C     Conributions from Boris A. Belov and Vladimir A. Popov (both at
C     IZMIRAN), as well as from Therese Moretto (DMI, DSRI, now at
C     NASA/GSFC).
C     The original version of this code is described in the brochure by
C   N.A. Tsyganenko, A.V. Usmanov, V.O. Papitashvili, N.E. Papitashvili,
C     and V.A. Popov, Software for computations of geomagnetic field and
C  related coordinate systems, Soviet Geophys. Committ., Moscow, 58 pp.,
C     1987. A number of subroutines from the revised GEOPACK-96 software
C     package developed by Nikolai A. Tsyganenko and Mauricio Peredo are
C     utilized in this code with some modifications (see full version of
C    GEOPACK-96 on http://www-spof.gsfc.nasa.gov/Modeling/geopack.html).
C     This code consists of the main subroutine GEOCGM99, five functions
C  (OVL_ANG, CGMGLA, CGMGLO, DFRIDR, and AZM_ANG), eigth new and revised
C     subroutines from the above-mentioned brochure (MLTUT, MFC, FTPRNT,
C     GEOLOW, CORGEO, GEOCOR, SHAG, and RIGHT), and 9 subroutines from
C   GEOPACK-96 (IGRF, SPHCAR, BSPCAR, GEOMAG, MAGSM, SMGSM, RECALC, SUN)
C  =====================================================================
C     Input parameters:
C     icor = +1    geo to cgm
C            -1    cgm to geo
C     iyr  = year
C     hi   = altitude
C     slar = geocentric latitude
C     slor = geocentric longitude (east +)
C     These two pairs can be either input or output parameters
C     clar = cgm latitude
C     clor = cgm longitude (east +)
C     Output parameters:
C     Array DAT(11,4) consists of 11 parameters (slar, slor, clar, clor,
C     rbm, btr, bfr, brr, ovl, azm, utm) organized for the start point
C     (*,1), its conjugate point (*,2), then for the footprints at 1-Re
C     of the start (*,3) and conjugate (*,4) points
C     Description of parameters used in the subroutine:
C     slac = conjugate geocentric latitude
C     sloc = conjugate geocentric longitude
C     slaf = footprint geocentric latitude
C     slof = footprint geocentric longitude
C     rbm  = apex of the magnetic field line in Re (Re=6371.2 km)
C            (this parameter approximately equals the McIlwain L-value)
C     btr  = IGRF Magnetic field H (nT)
C     bfr  = IGRF Magnetic field D (deg)
C     brr  = IGRF Magnetic field Z (nT)
C     ovl  = oval_angle as the azimuth to "magnetic north":
C                + east in Northern Hemisphere
C                + west in Southern Hemisphere
C     azm  = meridian_angle as the azimuth to the CGM pole:
C                + east in Northern Hemisphere
C                + west in Southern Hemisphere
C     utm  = magnetic local time (MLT) midnight in UT hours
C     pla  = array of geocentric latitude and

C     plo  = array of geocentric longitudes for the CGM poles
C            in the Northern and Southern hemispheres at a given
C            altitude (indices 1 and 2) and then at the Earth's
C            surface - 1-Re or zero altitude - (indices 3 and 4)
C     dla  = dipole latitude
C     dlo  = dipole longitude
C  =====================================================================
C     Year (for example, as for Epoch 1995.0 - no fraction of the year)
C     .. Scalar Arguments ..
      DOUBLE PRECISION HI
      INTEGER ICOR,IYEAR
C     ..
C     .. Array Arguments ..
      DOUBLE PRECISION DAT(11,4),PLA(4),PLO(4)
C     ..



***********************************************************************

      DOUBLE PRECISION FUNCTION OVL_ANG(SLA,SLO,CLA,CLO,RR)
C   This function returns an estimate at the given location of the angle
C     (oval_angle) between the directions (tangents) along the constant
C     CGM and geographic latitudes by utilizing the function DFRIDR from
C     Numerical Recipes for FORTRAN.
C   This angle can be taken as the azimuth to the local "magnetic" north
C   (south) if the eastward (westward) tangent to the local CGM latitude
C     points south (north) from the local geographic latitude.
C  Written by Therese Moretto in August 1994 (revised by V. Papitashvili
C     in January 1999).
C   Ignore points which nearly coincide with the geographic or CGM poles
C    within 0.01 degree in latitudes; this also takes care if SLA or CLA
C     are dummy values (e.g., 999.99)
C     .. Scalar Arguments ..
      DOUBLE PRECISION CLA,CLO,RR,SLA,SLO
C     ..



***********************************************************************

      DOUBLE PRECISION FUNCTION CGMGLA(CLON)
C     This function returns the geocentric latitude as a function of CGM
C     longitude with the CGM latitude held in common block CGMGEO.
C     Essentially this function just calls the subroutine CORGEO.
C     .. Scalar Arguments ..
      DOUBLE PRECISION CLON
C     ..



***********************************************************************

      DOUBLE PRECISION FUNCTION CGMGLO(CLON)
C     Same as the function CGMGLA but this returns the geocentric
C    longitude. If cr360 is true, geolon+360 deg is returned when geolon
C     is less than 90 deg. If cr0 is true, geolon-360 deg is returned
C     when geolon is larger than 270 degrees.
C     .. Scalar Arguments ..
      DOUBLE PRECISION CLON
C     ..



***********************************************************************

      DOUBLE PRECISION FUNCTION AZM_ANG(SLA,SLO,CLA,PLA,PLO)
C     Computation of an angle between the north geographic meridian and
C     direction to the North (South) CGM pole: positive azimuth is
C     measured East (West) from geographic meridian, i.e., the angle is
C     measured between the great-circle arc directions to the geographic
C     and CGM poles. In this case the geomagnetic field components in
C     XYZ (NEV) system can be converted into the CGM system in both
C     hemispheres as:
C                           XM = X COS(alf) + Y SIN(alf)
C                           YM =-X SIN(alf) + Y COS(alf)
C   Written by V. O. Papitashvili in mid-1980s; revised in February 1999
C   Ignore points which nearly coincide with the geographic or CGM poles
C    within 0.01 degree in latitudes; this also takes care if SLA or CLA
C     are dummy values (e.g., 999.99)
C     .. Scalar Arguments ..
      DOUBLE PRECISION CLA,PLA,PLO,SLA,SLO
C     ..



***********************************************************************

      SUBROUTINE MLTUT(SLA,SLO,CLA,PLA,PLO,UT)
C     Calculates the MLT midnight in UT hours
C     Definition of the MLT midnight (MLTMN) here is different from the
C     approach described elsewhere. This definition does not take into
C    account the geomagnetic meridian of the subsolar point which causes
C    seasonal variations of the MLTMN in UT time. The latter approach is
C     perfectly applicable to the dipole or eccentric dipole magnetic
C    coordinates but it fails with the CGM coordinates because there are
C     forbidden areas near the geomagnetic equator where CGM coordinates
C     cannot be calculated by definition [e.g., Gustafsson et al., JATP,
C     54, 1609, 1992].
C  In this code the MLT midnight is defined as location of a given point
C    on (or above) the Earth's surface strictly behind the North (South)
C     CGM pole in such the Sun, the pole, and the point are lined up.
C     This approach was originally proposed and coded by Boris Belov
C     sometime in the beginning of 1980s; here it is slightly edited by
C     Vladimir Papitashvili in February 1999.
C   Ignore points which nearly coincide with the geographic or CGM poles
C    within 0.01 degree in latitudes; this also takes care if SLA or CLA
C     are dummy values (e.g., 999.99)
C     .. Scalar Arguments ..
      DOUBLE PRECISION CLA,PLA,PLO,SLA,SLO,UT
C     ..



***********************************************************************

      SUBROUTINE MFC(SLA,SLO,R,H,D,Z)
C     Computation of the IGRF magnetic field components
C     Extracted as a subroutine from the earlier version of GEO-CGM.FOR
C     V. Papitashvili, February 1999
C     This takes care if SLA or CLA are dummy values (e.g., 999.99)
C     .. Scalar Arguments ..
      DOUBLE PRECISION D,H,R,SLA,SLO,Z
C     ..



***********************************************************************

      SUBROUTINE FTPRNT(RH,SLA,SLO,CLA,CLO,ACLA,ACLO,SLAF,SLOF,RF)
C     Calculation of the magnetic field line footprint at the Earth's
C     (or any higher) surface.
C   Extracted as a subroutine from the earlier version of GEO-CGM.FOR by
C   V. Papitashvili in February 1999 but then the subroutine was revised
C    to obtain the Altitude Adjusted CGM coordinates. The AACGM approach
C     is proposed by Kile Baker of the JHU/APL, see their World Wide Web
C     site http://sd-www.jhuapl.edu/RADAR/AACGM/ for details.
C     If RF = 1-Re (i.e., at the Earth's surface), then the footprint
C     location is defined as the Altitude Adjusted (AA) CGM coordinates
C     for a given point (ACLA, ACLO).
C     If RF = 1.xx Re (i.e., at any altitude above or below the starting
C     point), then the conjunction between these two points can be found
C     along the field line.
C     This takes care if SLA or CLA are dummy values (e.g., 999.99)
C     .. Scalar Arguments ..
      DOUBLE PRECISION ACLA,ACLO,CLA,CLO,RF,RH,SLA,SLAF,SLO,SLOF
C     ..



***********************************************************************

      SUBROUTINE GEOLOW(SLAR,SLOR,RH,CLAR,CLOR,RBM,SLAC,SLOC)
C     Calculates CGM coordinates from geocentric ones at low latitudes
C    where the DGRF/IGRF magnetic field lines may never cross the dipole
C     equatorial plane and, therefore, the definition of CGM coordinates
C     becomes invalid.
C     The code is written by Natalia and Vladimir Papitashvili as a part
C   of the earlier versions of GEO-CGM.FOR; extracted as a subroutine by
C     V. Papitashvili in February 1999.
C     Apr 11, 2001  GEOLOW is modified to account for interpolation of
C                CGM meridians near equator across the 360/0 boundary
C     See the paper by  Gustafsson, G., N. E. Papitashvili, and V. O.
C    Papitashvili, A revised corrected geomagnetic coordinate system for
C     Epochs 1985 and 1990 [J. Atmos. Terr. Phys., 54, 1609-1631, 1992]
C     for detailed description of the B-min approach utilized here.
C     This takes care if SLA is a dummy value (e.g., 999.99)
C     .. Scalar Arguments ..
      DOUBLE PRECISION CLAR,CLOR,RBM,RH,SLAC,SLAR,SLOC,SLOR
C     ..



***********************************************************************

      SUBROUTINE CORGEO(SLA,SLO,RH,DLA,DLO,CLA,CLO,PMI)
C     Calculates geocentric coordinates from corrected geomagnetic ones.
C     The code is written by Vladimir Popov and Vladimir Papitashvili
C     in mid-1980s; revised by V. Papitashvili in February 1999
C     This takes care if CLA is a dummy value (e.g., 999.99)
C     .. Scalar Arguments ..
      DOUBLE PRECISION CLA,CLO,DLA,DLO,PMI,RH,SLA,SLO
C     ..



***********************************************************************

      SUBROUTINE GEOCOR(SLA,SLO,RH,DLA,DLO,CLA,CLO,PMI)
C     Calculates corrected geomagnetic coordinates from geocentric ones
C     The code is written by Vladimir Popov and Vladimir Papitashvili
C     in mid-1980s; revised by V. Papitashvili in February 1999
C     This takes care if SLA is a dummy value (e.g., 999.99)
C     .. Scalar Arguments ..
      DOUBLE PRECISION CLA,CLO,DLA,DLO,PMI,RH,SLA,SLO
C     ..



***********************************************************************

      SUBROUTINE SHAG(X,Y,Z,DS)
C     Similar to SUBR STEP from GEOPACK-1996 but SHAG takes into account
C     only internal sources
C     The code is re-written from Tsyganenko's subroutine STEP by
C     Natalia and Vladimir Papitashvili in mid-1980s
C     .. Scalar Arguments ..
      DOUBLE PRECISION DS,X,Y,Z
C     ..



***********************************************************************

      SUBROUTINE RIGHT(X,Y,Z,R1,R2,R3)
C   Similar to SUBR RHAND from GEOPACK-1996 but RIGHT takes into account
C     only internal sources
C     The code is re-written from Tsyganenko's subroutine RHAND
C     by Natalia and Vladimir Papitashvili in mid-1980s
C     .. Scalar Arguments ..
      DOUBLE PRECISION R1,R2,R3,X,Y,Z
C     ..



***********************************************************************

      SUBROUTINE IGRF(IY,NM,R,T,F,BR,BT,BF)
c     Jan 20, 2001: Subroutine IGRF is modified by V. Papitashvili - SHA
c     coefficients for IGRF-2000, and SV 2000-2005 are added (note that
c     IGRF-1995 has not been changed to DGRF-1995 this time
c     (see http://www.ngdc.noaa.gov/IAGA/wg8/igrf2000.html)
c     Aug 26, 1997: Subroutine IGRF is modified by V. Papitashvili - SHA
c     coefficients for DGRF-1990, IGRF-1995, and SV 1995-2000 are added
c     (EOS, v.77, No.16, p.153, April 16, 1996)
c   Feb 03, 1995: Modified by Vladimir Papitashvili (SPRL, University of
c     Michigan) to accept dates between 1945 and 2000
C  MODIFIED TO ACCEPT DATES BETWEEN 1965 AND 2000; COEFFICIENTS FOR IGRF
C     1985 HAVE BEEN REPLACED WITH DGRF1985 COEFFICIENTS [EOS TRANS. AGU
C     APRIL 21, 1992, C  P. 182]. ALSO, THE CODE IS MODIFIED TO ACCEPT
C    DATES BEYOND 1990, AND TO USE LINEAR EXTRAPOLATION BETWEEN 1990 AND
C     2000 BASED ON THE IGRF COEFFICIENTS FROM THE SAME EOS ARTICLE
C   Modified by Mauricio Peredo, Hughes STX at NASA/GSFC, September 1992
C     CALCULATES COMPONENTS OF MAIN GEOMAGNETIC FIELD IN SPHERICAL
C     GEOCENTRIC COORDINATE SYSTEM BY USING THIRD GENERATION IGRF MODEL
C     (J. GEOMAG. GEOELECTR. V.34, P.313-315, 1982; GEOMAGNETISM AND
C     AERONOMY V.26, P.523-525, 1986).
C     UPDATING OF COEFFICIENTS TO A GIVEN EPOCH IS MADE DURING THE FIRST
C     CALL AND AFTER EVERY CHANGE OF PARAMETER IY
C     ---INPUT PARAMETERS:
C     IY - YEAR NUMBER (FROM 1945 UP TO 1990)
C  NM - MAXIMAL ORDER OF HARMONICS TAKEN INTO ACCOUNT (NOT MORE THAN 10)
C   R,T,F - SPHERICAL COORDINATES OF THE POINT (R IN UNITS RE=6371.2 KM,
C     COLATITUDE T AND LONGITUDE F IN RADIANS)
C     ---OUTPUT PARAMETERS:
C     BR,BT,BF - SPHERICAL COMPONENTS OF MAIN GEOMAGNETIC FIELD (in nT)
C    AUTHOR: NIKOLAI A. TSYGANENKO, INSTITUTE OF PHYSICS, ST.-PETERSBURG
C      STATE UNIVERSITY, STARY PETERGOF 198904, ST.-PETERSBURG, RUSSIA
C      (now the NASA Goddard Space Fligth Center, Greenbelt, Maryland)
      IMPLICIT NONE
C     G0, G1, and H1 are used in SUBROUTINE DIP to calculate geodipole's
C     moment for a given year
C     .. Scalar Arguments ..
      DOUBLE PRECISION BF,BR,BT,F,R,T
      INTEGER IY,NM
C     ..



***********************************************************************

      SUBROUTINE RECALC(IYR,IDAY,IHOUR,MIN,ISEC)
C     THIS IS A MODIFIED VERSION OF THE SUBROUTINE RECOMP WRITTEN BY
C     N. A. TSYGANENKO. SINCE I WANT TO USE IT IN PLACE OF SUBROUTINE
C     RECALC, I HAVE RENAMED THIS ROUTINE RECALC AND ELIMINATED THE
C     ORIGINAL RECALC FROM THIS VERSION OF THE  PACKAGE.
C    THIS WAY ALL ORIGINAL CALLS TO RECALC WILL CONTINUE TO WORK WITHOUT
C     HAVING TO CHANGE THEM TO CALLS TO RECOMP.
C     AN ALTERNATIVE VERSION OF THE SUBROUTINE RECALC FROM THE GEOPACK
C     PACKAGE BASED ON A DIFFERENT APPROACH TO DERIVATION OF ROTATION
C     MATRIX ELEMENTS
C     THIS SUBROUTINE WORKS BY 20% FASTER THAN RECALC AND IS EASIER TO
C     UNDERSTAND
C     #####################################################
C     #  WRITTEN BY  N.A. TSYGANENKO ON DECEMBER 1, 1991  #
C     #####################################################
C     Modified by Mauricio Peredo, Hughes STX at NASA/GSFC Code 695,
C     September 1992
c    Modified to accept years up to 2005 (V. Papitashvili, January 2001)
c  Modified to accept dates up to year 2000 and updated IGRF coeficients
c     from 1945 (updated by V. Papitashvili, February 1995)
C     OTHER SUBROUTINES CALLED BY THIS ONE: SUN
C     IYR = YEAR NUMBER (FOUR DIGITS)
C     IDAY = DAY OF YEAR (DAY 1 = JAN 1)
C     IHOUR = HOUR OF DAY (00 TO 23)
C     MIN = MINUTE OF HOUR (00 TO 59)
C     ISEC = SECONDS OF DAY(00 TO 59)
      IMPLICIT NONE
C     .. Scalar Arguments ..
      INTEGER IDAY,IHOUR,ISEC,IYR,MIN
C     ..



***********************************************************************

      SUBROUTINE SUN(IYR,IDAY,IHOUR,MIN,ISEC,GST,SLONG,SRASN,SDEC)
C    CALCULATES FOUR QUANTITIES NECESSARY FOR COORDINATE TRANSFORMATIONS
C     WHICH DEPEND ON SUN POSITION (AND, HENCE, ON UNIVERSAL TIME AND
C     SEASON)
C     ---INPUT PARAMETERS:
C     IYR,IDAY,IHOUR,MIN,ISEC - YEAR, DAY, AND UNIVERSAL TIME IN HOURS,
C     MINUTES, AND SECONDS  (IDAY=1 CORRESPONDS TO JANUARY 1).
C     ---OUTPUT PARAMETERS:
C   GST - GREENWICH MEAN SIDEREAL TIME, SLONG - LONGITUDE ALONG ECLIPTIC
C     SRASN - RIGHT ASCENSION,  SDEC - DECLINATION  OF THE SUN (RADIANS)
C     THIS SUBROUTINE HAS BEEN COMPILED FROM:
C     RUSSELL C.T., COSM.ELECTRODYN., 1971, V.2,PP.184-196.
C     AUTHOR: Gilbert D. Mead
      IMPLICIT NONE
C     .. Scalar Arguments ..
      DOUBLE PRECISION GST,SDEC,SLONG,SRASN
      INTEGER IDAY,IHOUR,ISEC,IYR,MIN
C     ..



***********************************************************************

      SUBROUTINE SPHCAR(R,TETA,PHI,X,Y,Z,J)
C     CONVERTS SPHERICAL COORDS INTO CARTESIAN ONES AND VICA VERSA
C     (TETA AND PHI IN RADIANS).
C                  J>0            J<0
C     -----INPUT:   J,R,TETA,PHI     J,X,Y,Z
C     ----OUTPUT:      X,Y,Z        R,TETA,PHI
C    AUTHOR: NIKOLAI A. TSYGANENKO, INSTITUTE OF PHYSICS, ST.-PETERSBURG
C      STATE UNIVERSITY, STARY PETERGOF 198904, ST.-PETERSBURG, RUSSIA
C      (now the NASA Goddard Space Fligth Center, Greenbelt, Maryland)
      IMPLICIT NONE
C     .. Scalar Arguments ..
      DOUBLE PRECISION PHI,R,TETA,X,Y,Z
      INTEGER J
C     ..



***********************************************************************

      SUBROUTINE BSPCAR(TETA,PHI,BR,BTET,BPHI,BX,BY,BZ)
C     CALCULATES CARTESIAN FIELD COMPONENTS FROM SPHERICAL ONES
C     -----INPUT:   TETA,PHI - SPHERICAL ANGLES OF THE POINT IN RADIANS
C              BR,BTET,BPHI -  SPHERICAL COMPONENTS OF THE FIELD
C     -----OUTPUT:  BX,BY,BZ - CARTESIAN COMPONENTS OF THE FIELD
C    AUTHOR: NIKOLAI A. TSYGANENKO, INSTITUTE OF PHYSICS, ST.-PETERSBURG
C      STATE UNIVERSITY, STARY PETERGOF 198904, ST.-PETERSBURG, RUSSIA
C      (now the NASA Goddard Space Fligth Center, Greenbelt, Maryland)
      IMPLICIT NONE
C     .. Scalar Arguments ..
      DOUBLE PRECISION BPHI,BR,BTET,BX,BY,BZ,PHI,TETA
C     ..



***********************************************************************

      SUBROUTINE GEOMAG(XGEO,YGEO,ZGEO,XMAG,YMAG,ZMAG,J,IYR)
C   CONVERTS GEOCENTRIC (GEO) TO DIPOLE (MAG) COORDINATES OR VICA VERSA.
C     IYR IS YEAR NUMBER (FOUR DIGITS).
C                           J>0                J<0
C     -----INPUT:  J,XGEO,YGEO,ZGEO,IYR   J,XMAG,YMAG,ZMAG,IYR
C     -----OUTPUT:    XMAG,YMAG,ZMAG        XGEO,YGEO,ZGEO
C    AUTHOR: NIKOLAI A. TSYGANENKO, INSTITUTE OF PHYSICS, ST.-PETERSBURG
C      STATE UNIVERSITY, STARY PETERGOF 198904, ST.-PETERSBURG, RUSSIA
C      (now the NASA Goddard Space Fligth Center, Greenbelt, Maryland)
      IMPLICIT NONE
C     .. Scalar Arguments ..
      DOUBLE PRECISION XGEO,XMAG,YGEO,YMAG,ZGEO,ZMAG
      INTEGER IYR,J
C     ..



***********************************************************************

      SUBROUTINE MAGSM(XMAG,YMAG,ZMAG,XSM,YSM,ZSM,J)
C CONVERTS DIPOLE (MAG) TO SOLAR MAGNETIC (SM) COORDINATES OR VICA VERSA
C                    J>0              J<0
C     -----INPUT: J,XMAG,YMAG,ZMAG     J,XSM,YSM,ZSM
C     ----OUTPUT:    XSM,YSM,ZSM       XMAG,YMAG,ZMAG
C  ATTENTION: SUBROUTINE RECALC MUST BE CALLED BEFORE MAGSM IN TWO CASES
C     /A/  BEFORE THE FIRST USE OF MAGSM
C     /B/  IF THE CURRENT VALUES OF IYEAR,IDAY,IHOUR,MIN,ISEC ARE
C          DIFFERENT FROM THOSE IN THE PRECEDING CALL OF  MAGSM
C    AUTHOR: NIKOLAI A. TSYGANENKO, INSTITUTE OF PHYSICS, ST.-PETERSBURG
C      STATE UNIVERSITY, STARY PETERGOF 198904, ST.-PETERSBURG, RUSSIA
C      (now the NASA Goddard Space Fligth Center, Greenbelt, Maryland)
      IMPLICIT NONE
C     .. Scalar Arguments ..
      DOUBLE PRECISION XMAG,XSM,YMAG,YSM,ZMAG,ZSM
      INTEGER J
C     ..



***********************************************************************

      SUBROUTINE SMGSM(XSM,YSM,ZSM,XGSM,YGSM,ZGSM,J)
C CONVERTS SOLAR MAGNETIC (SM) TO SOLAR MAGNETOSPHERIC (GSM) COORDINATES
C     OR VICA VERSA.
C                  J>0                 J<0
C     -----INPUT: J,XSM,YSM,ZSM        J,XGSM,YGSM,ZGSM
C     ----OUTPUT:  XGSM,YGSM,ZGSM       XSM,YSM,ZSM
C  ATTENTION: SUBROUTINE RECALC MUST BE CALLED BEFORE SMGSM IN TWO CASES
C     /A/  BEFORE THE FIRST USE OF SMGSM
C     /B/  IF THE CURRENT VALUES OF IYEAR,IDAY,IHOUR,MIN,ISEC ARE
C          DIFFERENT FROM THOSE IN THE PRECEDING CALL OF SMGSM
C    AUTHOR: NIKOLAI A. TSYGANENKO, INSTITUTE OF PHYSICS, ST.-PETERSBURG
C      STATE UNIVERSITY, STARY PETERGOF 198904, ST.-PETERSBURG, RUSSIA
C      (now the NASA Goddard Space Fligth Center, Greenbelt, Maryland)
      IMPLICIT NONE
C     .. Scalar Arguments ..
      DOUBLE PRECISION XGSM,XSM,YGSM,YSM,ZGSM,ZSM
      INTEGER J
C     ..



***********************************************************************

      DOUBLE PRECISION FUNCTION TNF(TI,TE,NE,NHP,NO,NH,NN2,NO2,NHE,IER)
C
C     TNF calculates the ion temperature (tn) given the electron and
C     neutral temperatures (TE, TN) and the electron, o+, h+, o, h,
C     n2, o2 and he concentrations (NE, NOP, NHP, NO, NH, NN2, NO2,
C     NHE). o+ and h+ ions are assumed to be the only ions present.
C     only coulomb collisions and ion-neutral polarization and
C     charge-exchange interactions are considered in balancing
C     the ion gas energy source and sink terms.  the solution for
C     tn is an iterative procedure in which TI is the initial value
C     of tn for the first iteration. all concentrations are in
C     units of cm**-3.
C
C       Input:
C          TI - Ion Temperature (K)
C          TE - Electron Temperature
C          NE - Electron concentration (cm**-3)
C         NHP - H Ion concentration
C          NO - O Ion concentration
C         NN2 - N2 concentration
C         NO2 - O2 concentration
C         NHE - HE concentration
C
C      Output:
C         IER - If (IER.NE.0) an error has occurred.
C
C     .. Scalar Arguments ..
      DOUBLE PRECISION NE,NH,NHE,NHP,NN2,NO,NO2,TE,TI
      INTEGER IER
C     ..



***********************************************************************

      SUBROUTINE VADD(A,B,C)
C
C     jmh - 11/79  ans fortran 66
C
C     VADD calculates the sum of two vectors A and B, C = A + B.
C
C       Input:
C          A - floating point vector of dimension 3.
C          B - floating point vector of dimension 3.
C
C      Output:
C          C - floating point vector of dimension 3.
C
C     .. Array Arguments ..
      DOUBLE PRECISION A(3),B(3),C(3)
C     ..
      C(1) = A(1) + B(1)
      C(2) = A(2) + B(2)
      C(3) = A(3) + B(3)
      RETURN
C
      END


***********************************************************************

      SUBROUTINE VCTCNV(FX,FY,FZ,X,Y,Z,FR,FT,FP,R,THETA,PHI,IMODE)
C
C     jmh - 11/79  ans fortran 66
C
C     VCTCNV converts between the cartesian and spherical coordinate
C     representations of a vector field f. (FX,FY,FZ) are the
C     components of the field at (X,Y,Z). (FR,FT,FP) are the
C     components of the field at (R,THETA,PHI) in the directions of
C     increasing R, increasing THETA and increasing PHI. if IMODE=1,
C     (FX,FY,FZ,X,Y,Z) -> (FR,FT,FP,R,THETA,PHI). if IMODE=2,
C     (FR,FT,FP,R,THETA,PHI) -> (FX,FY,FZ,X,Y,Z). THETA and PHI are
C     in degrees.
C
C       Input:
C         IMODE - 1 cartesian to spherical
C                 2 spherical to cartesian
C
C     Input, output:
C      FX,FY,FZ - cartesian vector field components
C         X,Y,Z - cartesian vector field coordinates
C      FR,FT,FP - spherical vector field components
C   R,THETA,PHI - spherical vector field coordinates
C
C     .. Scalar Arguments ..
      DOUBLE PRECISION FP,FR,FT,FX,FY,FZ,PHI,R,THETA,X,Y,Z
      INTEGER IMODE
C     ..



***********************************************************************

      DOUBLE PRECISION FUNCTION VMAG(A)
C
C     jmh - 1/80  ans fortran 66
C
C     VMAG calculates the magnitude of a vector A
C
C     Input:
C        A - floating point vector of dimension 3
C
C     .. Array Arguments ..
      DOUBLE PRECISION A(3)
C     ..
C     .. Intrinsic Functions ..
      INTRINSIC DSQRT
C     ..
      VMAG = DSQRT(A(1)*A(1)+A(2)*A(2)+A(3)*A(3))
      RETURN
C
      END


***********************************************************************

      SUBROUTINE VSUB(A,B,C)
C
C     jmh - 11/79  ans fortran 66
C
C     VSUB calculates the difference of two vectors A and B, C = A - B.
C
C     Input:
C        A - floating point vector of dimension 3
C        B - floating point vector of dimension 3
C
C     Output:
C        C - floating point vector of dimension 3
C
C     .. Array Arguments ..
      DOUBLE PRECISION A(3),B(3),C(3)
C     ..
      C(1) = A(1) - B(1)
      C(2) = A(2) - B(2)
      C(3) = A(3) - B(3)
      RETURN
C
      END


***********************************************************************

      SUBROUTINE CALNDR(YEAR,DAYNO,DAY,MONTH,IER)
C
C     CALNDR returns DAY, MONTH and IER. THE FOLLOWING VARIABLES APPEAR
C     IN THE CALLING SEQUENCE (all INTEGERS).
C
C          YEAR - YEAR (1977)
C         DAYNO - DAY OF THE YEAR (60)
C           DAY - DAY OF THE MONTH (1)
C         MONTH - MONTH NUMBER (3)
C           IER - ERROR INDICATOR. IER IS RETURNED BY ALL ROUTINES IN
C                 THE PACKAGE. IER=0 IF NO ERRORS ARE DETECTED. IER=1
C                 IF AN ERROR IS DETECTED.
C
C     .. Scalar Arguments ..
      INTEGER DAY,DAYNO,IER,MONTH,YEAR
C     ..



***********************************************************************

      SUBROUTINE DATER(DATE,NCHAR,DAY,MONTH,YEAR,IER)
C
C     SUBROUTINES DATER, MONAME, MONUM, WKNAME, IDAY, CALNDR, JDAY,
C     JDATER, IZLR, IDMYCK AND DATES COMPRISE A COMPREHENSIVE DATE 
C     MANIPULATION
C     PACKAGE. THE FOLLOWING VARIABLES APPEAR IN THE CALLING SEQUENCE
C     OF ONE OR MORE SUBROUTINES. ALL ARE TYPED INTEGER. THE VALUE
C     CORRESPONDING TO MARCH 1, 1977 IS SHOWN IN PARENTHESES.
C
C           DAY - DAY OF THE MONTH (1)
C         MONTH - MONTH NUMBER (3)
C          YEAR - YEAR (1977)
C         DAYNO - DAY OF THE YEAR (60)
C        JDAYNO - JULIAN DAY NUMBER (2443204)
C          WDAY - WEEKDAY NUMBER (3)
C          DATE - DATE AS A STRING OF ALPHANUMERIC CHARS, 3 CHARS/WORD.
C                 WHEN AN INPUT VARIABLE, DATE MAY BE IN ANY REASONABLE
C                 FORMAT, E.G. MARCH 1 1977, 3/1/77, ETC. IF EXPRESSED
C                 AS THREE NUMERIC FIELDS, ORDER IS PRESUMED TO BE
C                 MONTH, DAY, YEAR. WHEN AN OUTPUT VARIABLE, THE FORMAT
C                 IS DETERMINED BY IOPT, AND SIX WORDS SHOULD BE
C                 RESERVED IN THE CALLING PROGRAM.
C
C          MSTR - MONTH, AS A STRING OF UPPER CASE ALPHABETIC
C                 CHARACTERS, 3 CHARACTERS/WORD. THREE WORDS SHOULD BE
C                 RESERVED IN THE CALLING PROGRAM. (MARCH)
C          WSTR - DAY OF THE WEEK, AS A STRING OF UPPER CASE ALPHABETIC
C                 CHARACTERS, 3 CHARACTERS/WORD. THREE WORDS SHOULD BE
C                 RESERVED IN THE CALLING PROGRAM. (TUESDAY)
C          IOPT - FORMAT INDICATOR WHEN DATE IS AN OUTPUT VARIABLE
C                  1 - 03/01/77
C                  3 - 01,MAR,1977
C                  4 - 1 MARCH, 1977
C                  5 - MARCH 1, 1977
C         NCHAR - NUMBER OF CHARACTERS TO BE SCANNED IN AN INPUT STRING.
C           IER - ERROR INDICATOR. IER IS RETURNED BY ALL ROUTINES IN
C                 THE PACKAGE. IER=0 IF NO ERRORS ARE DETECTED. IER=1
C                 IF AN ERROR IS DETECTED.
C
C     .. Scalar Arguments ..
      INTEGER DAY,IER,MONTH,NCHAR,YEAR
      CHARACTER*(*) DATE
C     ..



***********************************************************************

      SUBROUTINE DATES(DAY,MONTH,YEAR,IOPT,IER,DATE)
C
C     Returns DATE String in various formats given DAY, MONTH, YEAR
C     and IOPT (which specifies the desired format).
C
C     The following variables appear in the calling sequence
C
C       Input:
C           DAY - DAY OF THE MONTH (1)
C         MONTH - MONTH NUMBER (3)
C          YEAR - YEAR (1977)
C          IOPT - FORMAT INDICATOR WHEN DATE IS AN OUTPUT VARIABLE
C                  1 - 01/01/97
C                  2 - 01,JAN,1997
C                  3 - 1 JANUARY, 1997
C                  4 - JANUARY 1 1997
C                  5 - 01JAN97
C
C       Output:
C           IER - ERROR INDICATOR. IER IS RETURNED BY ALL ROUTINES IN
C                 THE PACKAGE. IER=0 IF NO ERRORS ARE DETECTED. IER=1
C                 IF AN ERROR IS DETECTED.
C          DATE - DATE AS A STRING OF ALPHANUMERIC CHARACTERS, 3
C                 CHARACTERS/WORD. WHEN AN INPUT VARIABLE, DATE MAY BE
C                 IN ANY REASONABLE FORMAT, E.G. MARCH 1 1977, 3/1/77,
C                 ETC. IF EXPRESSED AS THREE NUMERIC FIELDS, ORDER IS
C                 PRESUMED TO BE MONTH, DAY, YEAR. WHEN AN OUTPUT
C                 VARIABLE, THE FORMAT IS DETERMINED BY IOPT, AND SIX
C                 WORDS SHOULD BE RESERVED IN THE CALLING PROGRAM.
C
C     .. Scalar Arguments ..
      INTEGER DAY,IER,IOPT,MONTH,YEAR
      CHARACTER*(*) DATE
C     ..



***********************************************************************

      SUBROUTINE IDAY(DAY,MONTH,YEAR,DAYNO,IER)
C
C     Returns Day-of-year from DAY, MONTH, YEAR.
C
C      Input:
C        IDBFIL - Madrigal file number (see EXDCON)
C           DAY - Day of month (1-31)
C         MONTH - Month of year (1-12)
C          YEAR - Year (e.g. 1977)
C
C      Output:
C         DAYNO - Day-of-year (1-356)
C           IER - If (IER.NE.0) an error has occurred.
C
C     .. Scalar Arguments ..
      INTEGER DAY,DAYNO,IER,MONTH,YEAR
C     ..



***********************************************************************

      INTEGER FUNCTION IDMYCK(DAY,MONTH,YEAR)
C
C     Returns 1 if DAY, MONTH and YEAR are legal values, 0 otherwise.
C
C      Input:
C           DAY - Day of month (1-31)
C         MONTH - Month of year (1-12)
C          YEAR - Year (e.g. 1977)
C
C     .. Scalar Arguments ..
      INTEGER DAY,MONTH,YEAR
C     ..



***********************************************************************

      SUBROUTINE IZLR(DAY,MONTH,YEAR,WDAY,IER)
C
C     Returns day-of-the-week value from DAY, MONTH and YEAR.
C
C      Input:
C           DAY - Day of month (1-31)
C         MONTH - Month of year (1-12)
C          YEAR - Year (e.g. 1977)
C
C      Output:
C          WDAY - Day of week (1-7)
C           IER - If (IER.NE.0) an error has occurred.
C
C     .. Scalar Arguments ..
      INTEGER DAY,IER,MONTH,WDAY,YEAR
C     ..



***********************************************************************

      SUBROUTINE JDATER(JDAYNO,DAY,MONTH,YEAR,IER)
C
C     Returns DAY, MONTH, YEAR from Julian day (See JDAY for inverse).
C
C      Input:
C        JDAYNO - Julian day (e.g. 2447892)
C
C      Output:
C           DAY - Day of month (1-31)
C         MONTH - Month of year (1-12)
C          YEAR - Year (e.g. 1977)
C           IER - If (IER.NE.0) an error has occurred.
C
C     .. Scalar Arguments ..
      INTEGER DAY,IER,JDAYNO,MONTH,YEAR
C     ..



***********************************************************************

      SUBROUTINE JDAY(DAY,MONTH,YEAR,JDAYNO,IER)
C
C     Returns Julian day from DAY, MONTH, YEAR (See JDATER for
C     inverse).
C
C      Input:
C           DAY - Day of month (1-31)
C         MONTH - Month of year (1-12)
C          YEAR - Year (e.g. 1977)
C
C      Output:
C        JDAYNO - Julian day (e.g. 2447892)
C           IER - If (IER.NE.0) an error has occurred.
C
C     .. Scalar Arguments ..
      INTEGER DAY,IER,JDAYNO,MONTH,YEAR
C     ..



***********************************************************************

      SUBROUTINE MONAME(MONTH,MSTR,IER)
C
C     Returns Month string from Month integer.
C
C      Input:
C        MONTH - Month of year (1-12)
C
C      Output:
C         MSTR - Month string (e.g. 'JANUARY')
C          IER - If (IER.NE.0) an error has occurred.
C
C     .. Scalar Arguments ..
      INTEGER IER,MONTH
      CHARACTER*(*) MSTR
C     ..



***********************************************************************

      SUBROUTINE MONUM(MSTR,MONTH,IER)
C
C     Returns Month integer from Month string (case insignificant).
C
C      Input:
C         MSTR - Month string (e.g. 'January')
C
C      Output:
C        MONTH - Month of year (1-12)
C          IER - If (IER.NE.0) an error has occurred.
C
C     .. Scalar Arguments ..
      INTEGER IER,MONTH
      CHARACTER*(*) MSTR
C     ..



***********************************************************************

      SUBROUTINE WKNAME(WDAY,IER,WSTR)
C
C     Returns day of the week string from day of the week number.
C
C       Input:
C         WDAY - week day number (1-7)
C
C      Output:
C         IER  - If (IER.NE.0) an error has occurred.
C         WSTR - day of the week string (e.g. 'SUNDAY').
C
C     .. Scalar Arguments ..
      INTEGER IER,WDAY
      CHARACTER*(*) WSTR
C     ..

Previous: C API   Up: Madrigal developer's guide   Next: Tcl API