PRO rdpgm,ht,plat,plong,glat,glong,dblsize=dblsize

; This routine reads in arrays of PACE geomagnetic coordinate and
; the corresponding geographic coordinates from a file.  The arrays
; are assumed to be square and all of the same size.  The altitude
; to which the coordinates pertain is given in the extension of the
; file.  The size of the arrays may vary over time so the routine
; determines the appropriate edge length of the square arrays from
; the total length of the file.

root='/jasper/cnsr3_data1/pace'
pf=FINDFILE(root+'/pgm2geo.'+STRING(FIX(ht),FORMAT='(I3.3)'),COUNT=npf)
IF npf GT 0 THEN BEGIN
    ; Read the coordinate arrays directly from the file
    pf=pf(0)
    OPENR,pfunit,pf,/GET_LUN
    pfstat=FSTAT(pfunit)
    edglen=SQRT(pfstat.SIZE/(4*4))	; pf contains 4 square arrays of 4-byte REALs
    plat=FLTARR(edglen,edglen)
    plong=plat
    glat=plat
    glong=plat
    READU,pfunit,plat,plong,glat,glong
    CLOSE,pfunit
    FREE_LUN,pfunit

    ; Rotate arrays so prime PGM meridian runs from center to bottom
    ; of plong array
    plat=ROTATE(plat,3)
    plong=ROTATE(plong,3)
    glat=ROTATE(glat,3)
    glong=ROTATE(glong,3)

    ENDIF ELSE BEGIN
    ; Interpolate in the coordinates for 150, 300, and 450 km
    ; Use three-point Lagrange interpolation
    ; First get the 150-km coordinates and initialize the interpolates
    h0=150.
    h1=300.
    h2=450.
    c0=((ht-h1)*(ht-h2))/((h0-h1)*(h0-h2))
    rdpgm,h0,plt,pln,glt,gln
    glat=c0*glt
    glong=c0*gln
    ; Next get the 300-km coordinates and update the interpolates
    c1=((ht-h0)*(ht-h2))/((h1-h0)*(h1-h2))
    rdpgm,h1,plt,pln,glt,gln
    plat=plt
    plong=pln
    glat=TEMPORARY(glat)+c1*glt
    glong=TEMPORARY(glong)+c1*gln
    ; Finally get the 450-km coordinates and complete the interpolation
    c2=((ht-h0)*(ht-h1))/((h2-h0)*(h2-h1))
    rdpgm,h2,plt,pln,glt,gln
    glat=TEMPORARY(glat)+c2*glt
    glong=TEMPORARY(glong)+c2*gln
    ENDELSE

IF KEYWORD_SET(dblsize) THEN BEGIN
; Interpolate the coordinate arrays to double their size
    latlon_interp,plat,plong,iplat,iplong,2,/degrees,/stereo
    plat=iplat
    plong=iplong
    latlon_interp,glat,glong,iglat,iglong,2,/degrees,/stereo
    glat=iglat
    glong=iglong
    ENDIF

RETURN
END