PROGRAM TEST_GEOPACK
    IMPLICIT NONE
    
!   PROGRAM INPUTS
    REAL                 ::  XGSWI,YGSWI,ZGSWI,BYIMF,BZIMF,PDYN,DST
    INTEGER              ::  IYEAR,IDAY,IHOUR,IMIN,ISEC
    REAL, DIMENSION(12)  ::  R_INPUTS

!   SUBROUTINE OUTPUTS (IF YOU WANT TO CHANGE XX/YY/ZZ TO LARGER ARRAYS YOU MUST
!   ALSO CHANGE THE LMAX INPUT IN THE CALL TO TRACE_08)
    REAL                   ::  XGEOF,YGEOF,ZGEOF,XGSWF,YGSWF,ZGSWF
    REAL, DIMENSION(5000)  ::  XX,YY,ZZ

!   RADIANS TO DEGREE VARIABLES
    REAL                   ::  PI,RADIUS,REI,TRADI,PRADI,REF,TRADF,PRADF

!   PROGRAM OUTPUTS
    REAL              ::  RF,THETAF,PHIF
    INTEGER           ::  L
    CHARACTER(LEN=8)  ::  HEMISPHERE
    
!   SUBROUTINE CONSTANTS
    REAL                 ::  VGSEX,VGSEY,VGSEZ,DSMAX,ERR,RMAX,RMIN
    REAL, DIMENSION(10)  ::  PARMOD
    EXTERNAL             ::  T96_01,IGRF_GSW_08
    
!   COMMAND LINE TEMP STRING
    CHARACTER(LEN=256)  ::  ARGUMENTS
    
!   LOOPING VARIABLES
    INTEGER    ::  II
    REAL       ::  MAPTO
    
!   READ THE COMMAND LINE
    DO II=1,12
        CALL GET_COMMAND_ARGUMENT(II,ARGUMENTS)
        READ(ARGUMENTS,*) R_INPUTS(II)
    END DO

!   ASSIGN INPUTS FROM THE COMMAND LINE
    IYEAR  = R_INPUTS(1)
    IDAY   = R_INPUTS(2)
    IHOUR  = R_INPUTS(3)
    IMIN   = R_INPUTS(4)
    ISEC   = R_INPUTS(5)
    XGSWI  = R_INPUTS(6)
    YGSWI  = R_INPUTS(7)
    ZGSWI  = R_INPUTS(8)
    BYIMF  = R_INPUTS(9)
    BZIMF  = R_INPUTS(10)
    PDYN   = R_INPUTS(11)
    DST    = R_INPUTS(12)
    
!   DEFINE CONSTANTS USING THE SAME VALUES AS THE PYTHON DEFAULTS
    VGSEX  = -400.0
    VGSEY  = 0.0
    VGSEZ  = 0.0
    DSMAX  = 0.01
    ERR    = 0.000001
    RMAX   = 60.0
    RMIN   = 1.0
    PARMOD = [PDYN,DST,BYIMF,BZIMF,0.0,0.0,0.0,0.0,0.0,0.0]
    
    CALL RECALC_08(IYEAR,IDAY,IHOUR,IMIN,ISEC,VGSEX,VGSEY,VGSEZ)
    
    DO MAPTO=-1,1,2
        CALL TRACE_08(XGSWI,YGSWI,ZGSWI,MAPTO,DSMAX,ERR,RMAX,RMIN,0,PARMOD,&
             T96_01,IGRF_GSW_08,XGSWF,YGSWF,ZGSWF,XX,YY,ZZ,L,5000)
        CALL GEOGSW_08(XGEOF,YGEOF,ZGEOF,XGSWF,YGSWF,ZGSWF,-1)
        CALL SPHCAR_08(REF,TRADF,PRADF,XGEOF,YGEOF,ZGEOF,-1)
        
!       CONVERT OUTPUTS TO DEGREES LATITUDE AND KILOMETRES
        PI = 3.1415926535
        RADIUS = 6371.2
        THETAF = 90.0 - (TRADF * 180.0 / PI)
        PHIF = PRADF * 180.0 / PI
        RF = REF * RADIUS

        IF(MAPTO.EQ.-1) HEMISPHERE = "NORTHERN"
        IF(MAPTO.EQ.1)  HEMISPHERE = "SOUTHERN"
        
        PRINT *,HEMISPHERE," HEMISPHERE:",THETAF,PHIF,RF
    
    END DO
    
END PROGRAM