        program amie2aspen
        ! program takes AMIE ascii output files and extracts relevant
        ! data, writes output file for ASPEN runs

        ! original version Jul 2003, Chris Hackert
        ! version 1, May 2005, write POTEN in ASPEN magnetic coords

        use netcdf_module, only: nc_close
        implicit none
! includes
        integer :: ithmx, lonmx, mmx, nq, nqs, mxnmx, mcoef, ncoef,
     |     mrow,kmx,ithtrns
        include 'qset.h'

        include 'rdamie.h'
        
! variables
        real, parameter :: PI = 3.14159265358979323
        real, parameter :: dtor = PI/180.
        ! AMIE magnetic grid
        real, parameter :: dlatm = 90./real(ITHMX) * dtor
        real, parameter :: dlonm = 360./real(LONMX) * dtor

        ! source magnetic grid (AMIE * 2 for both hemispheres)
        integer, parameter :: imaxm = LONMX+1
        integer, parameter :: jmaxm = 2*(ITHTRNS+1) ! 2 hemispheres
        real, parameter :: dmlt = 24./LONMX  ! MLT interval in AMIEfile
        real :: ylatm(jmaxm)

        ! ASPEN geographic grid
        integer, parameter :: imaxg = 73
        integer, parameter :: jmaxg = 36
        real alatm(imaxg,jmaxg), alonm(imaxg,jmaxg)

        ! ASPEN magnetic grid
        integer, parameter :: zimxm = 80+1 
        integer, parameter :: zjmxm = 97
        real zlatm(zjmxm), zlonm(zimxm)
        
        !interpolation arrays
        integer im(imaxg,jmaxg), jm(imaxg,jmaxg)
        real    dim(imaxg,jmaxg),djm(imaxg, jmaxg)

        ! AMIE file info
        integer nfiles, index
        character*120 amiefileN, amiefileS
        integer iuan, iuas, numtot  ! file unit, number of times on file
        integer njh, nloop, mspall, iflgaur, ihppem, i0prsz, iblat,
     |     isolcn, ntm, iyear, imon, ida1, ida, ihr, imin, ndst, ierr
        real pppnt, pdifnt, oltxnt, omltxnt

        ! working vars
        integer i,j,nf, ncid
        real EKEVG(imaxg,jmaxg), EFLXG(imaxg,jmaxg), EPOTG(imaxg,jmaxg)
        real EKEVM(imaxm,jmaxm), EFLXM(imaxm,jmaxm), EPOTM(imaxm,jmaxm)
        real EPOTZ(zimxm,zjmxm)
        real xmlt, del, rot  ! MLT to Apex lon vars
        integer m, mp1       ! MLT to Apex lon vars
        ! NH and SH cusp loc in magnetic colat and long (radians)
        real*8 cuspnclt, cuspnlon, cuspsclt, cuspslon 
        real*8 ftimen, ftimes ! time of AMIE pattern in fractional day of year
        real temp

        ! functions
        integer nc_outfile
        external nc_outfile
        real*8 fday
        external fday

! begin execution

        ! read ASPEN magfield file for coordinates
        call rdfield(imaxg,jmaxg,alatm,alonm)

        ! setup ASPEN magnetic field coordinates
        call defmag(zimxm,zjmxm,zlatm,zlonm)

        ! Northern hemisphere AMIE magnetic latitudes
        do j=1,jmaxm/2
           ylatm(j) = PI/2. - dlatm * (j-1)
        enddo
        ! Southern hemisphere AMIE magnetic latitudes
        do j=1,jmaxm/2
           ylatm(j+jmaxm/2) = -PI/2. + dlatm * (jmaxm/2 - j)
        enddo

        ! setup geomagnetic to geographic transform using NCAR grdset
        call grdset(alatm,alonm,im,jm,dim,djm,imaxg,imaxg,jmaxg,
     |              imaxm,jmaxm,ylatm,-PI)

        ! prep output netcdf file
        ncid = nc_outfile('amievol.nc')

        ! open AMIE file, get header
        iuan = 50
        iuas = 51

        ! loop over AMIE file pairs (NH, SH)
        index = 1
        read(*,*) nfiles

        do nf = 1,nfiles
        
        read(*,*) amiefileN
        read(*,*) amiefileS

        !WARN: should have error checking that NH and SH files are
        !compatible (time increment, number of times, day of year, etc)
        open(unit=iuan, file=amiefileN, status='old')
        open(unit=iuas, file=amiefileS, status='old')
        call rdtamie(iuan,njh,nloop,numtot,mspall,iflgaur,ihppem,
     |               i0prsz,iblat,isolcn)
        call rdtamie(iuas,njh,nloop,numtot,mspall,iflgaur,ihppem,
     |               i0prsz,iblat,isolcn)

        ! loop over AMIE times
        do ntm=1,numtot

           ! read NH data into common blocks
           call rd1amie(iuan,njh,ntm ,iyear,imon,ida1,ida,ihr,imin,
     |                  pppnt, pdifnt, oltxnt, omltxnt, ndst, ierr)

           ftimen = fday(iyear,imon,ida,ihr,imin)
           ! Convert AMIE MLT output to APEX longitudes using
           !   AMIETERP method
           rot = 250./15. - ihr - imin/60.

           do i=1,imaxm
              xmlt = mod((i-1)*dmlt - rot + 24., 24.)
              m = ifix(xmlt/dmlt + 0.01)  ! m = 0 --> imaxm-1
              mp1 = m + 1
              if (mp1 .ge. imaxm) mp1 = 1 ! m=0 is same as m=imaxm-1
              del = xmlt - (m)*dmlt
              ! copy NH AMIE arrays into our arrays
              do j=1,jmaxm/2
                 EKEVM(i,j) = (1.-del)*EMKEV(m,j-1) + del*EMKEV(mp1,j-1)
                 EFLXM(i,j) = (1.-del)*WPM2(m,j-1) + del*WPM2(mp1,j-1)
                 EPOTM(i,j) = (1.-del)*EFPOT(m,j-1) + del*EFPOT(mp1,j-1)
              enddo
           enddo ! i

           ! NH cusp location (magnetic grid, radians)
           ! 12/1/03 Geoff says don't use AMIE oltxnt,omltxnt

           ! place cusp at 75 deg mag. lat. (15 deg colat)
           cuspnclt = 15. * PI/180.
           ! place cusp at min(abs(POTEN)) between 10 and 14 SLT
           ! find j index for CUSP
           j = int(cuspnclt/dlatm) + 1
           ! find index where POTEN changes sign (passes thru 0) on dayside
           i = int(6/dmlt)
           do while (i .lt. int(18/dmlt) 
     |                 .and. EFPOT(i,j)*EFPOT(i+1,j) .gt. 0)
             i=i+1
           enddo
           temp = (i+1) - efpot(i+1,j) / (efpot(i+1,j) - efpot(i,j))
           temp = temp*dmlt
           if (temp .lt. 10.) temp = 10.
           if (temp .gt. 14.) temp = 14.
           ! ASPEN wants cusp lon relative to subsolar longitude
           cuspnlon = PI/180. * ((temp - 12.) * 15.)
           if (cuspnlon .gt. PI) cuspnlon = cuspnlon - 2.*PI
           if (cuspnlon .lt. -PI) cuspnlon = cuspnlon + 2.*PI

           ! read SH data into common blocks
           call rd1amie(iuas,njh,ntm ,iyear,imon,ida1,ida,ihr,imin,
     |                  pppnt, pdifnt, oltxnt, omltxnt, ndst, ierr)

           ! check time
           ftimes = fday(iyear,imon,ida,ihr,imin)
           if (ifix(100*ftimen) .ne. ifix(100*ftimes)) then
                   write(*,*) 'NH ftime != SH ftime',ftimen, ftimes
                   STOP
           endif

           ! Convert AMIE MLT output to APEX longitudes using
           !   AMIETERP method
           rot = 250./15. - ihr - imin/60.

           do i=1,imaxm
              xmlt = mod((i-1)*dmlt - rot + 24., 24.)
              m = ifix(xmlt/dmlt + 0.01)  ! m = 0 --> imaxm-1
              mp1 = m + 1
              if (mp1 .ge. imaxm) mp1 = 1 ! m=0 is same as m=imaxm-1
              del = xmlt - (m)*dmlt
              ! copy SH AMIE arrays into our arrays
              do j=1,jmaxm/2
                 EKEVM(i,jmaxm-j+1) = (1.-del)*EMKEV(m,j-1) + 
     |                                del*EMKEV(mp1,j-1)
                 EFLXM(i,jmaxm-j+1) = (1.-del)*WPM2(m,j-1) + 
     |                                del*WPM2(mp1,j-1)
                 EPOTM(i,jmaxm-j+1) = (1.-del)*EFPOT(m,j-1) + 
     |                                del*EFPOT(mp1,j-1)
              enddo
           enddo ! i

           ! SH cusp location (magnetic grid, radians)
           ! NOTE: ASPEN wants SH colat from S pole
           ! 12/1/03 Geoff says don't use AMIE oltxnt,omltxnt

           ! place cusp at 75 deg mag. lat. (15 deg colat)
           cuspsclt = 15. * PI/180.
           ! place cusp at min(abs(POTEN)) between 10 and 14 SLT
           ! find j index for CUSP
           j = int(cuspsclt/dlatm) + 1
           ! find index where POTEN changes sign (passes thru 0) on dayside
           i = int(6/dmlt)
           do while (i .lt. int(18/dmlt) 
     |                 .and. EFPOT(i,j)*EFPOT(i+1,j) .gt. 0)
             i=i+1
           enddo
           temp = (i+1) - efpot(i+1,j) / (efpot(i+1,j) - efpot(i,j))
           temp = temp*dmlt
           if (temp .lt. 10.) temp = 10.
           if (temp .gt. 14.) temp = 14.
           ! ASPEN wants cusp lon relative to subsolar longitude
           cuspslon = PI/180. * ((temp - 12.) * 15.)
           if (cuspslon .gt. PI) cuspslon = cuspslon - 2.*PI
           if (cuspslon .lt. -PI) cuspslon = cuspslon + 2.*PI

           ! interpolate magnetic to geographic coordinates
           call grdtrp(EKEVM, EKEVG, im, jm, dim, djm, imaxg, imaxm,
     |                 imaxg, jmaxg, imaxm, jmaxm)
           call grdtrp(EFLXM, EFLXG, im, jm, dim, djm, imaxg, imaxm,
     |                 imaxg, jmaxg, imaxm, jmaxm)
! no longer converting POT to ASPEN geographic
!           call grdtrp(EPOTM, EPOTG, im, jm, dim, djm, imaxg, imaxm,
!     |                 imaxg, jmaxg, imaxm, jmaxm)

           ! do AMIE magnetic to ASPEN magnetic interpolation for POTEN
           call aterp(imaxm,jmaxm,ylatm,EPOTM,zimxm,zjmxm,zlatm,EPOTZ)

           ! output
           call nc_add(ncid,index, ftimen, cuspnclt, cuspnlon, cuspsclt,
     |                 cuspslon, ekevg, eflxg, epotz)

           index = index + 1
        enddo ! ntm
        enddo ! nf

        call nc_close(ncid)

        end

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

        integer function nc_outfile(filename)
        use netcdf_module, only: nc_open
        implicit none
        character(*) filename

        include "netcdf.inc"

        integer ncid, istat
        ! dimensions
        integer id_time, id_lon, id_lat, id_mln, id_mlt
        ! variables
        integer id_pot, id_eflx, id_ekev, id_ftim
        integer id_cuspnclt, id_cuspnlon, id_cuspsclt, id_cuspslon
        integer ids(3)
        character*80 char80

        common /netcdf/ id_time,id_lon,id_lat, id_ftim,id_pot,id_eflx,
     |     id_ekev, id_cuspnclt, id_cuspnlon, id_cuspsclt, id_cuspslon,
     |     id_mln, id_mlt

        ! ASPEN geographic grid
        integer, parameter :: imaxg = 73
        integer, parameter :: jmaxg = 36

        ! ASPEN magnetic grid
        integer, parameter :: zimxm = 80+1 
        integer, parameter :: zjmxm = 97


        ! create netcdf file
        call nc_open(ncid, filename, 'REPLACE','WRITE')

        ! define netcdf file
        istat = nf_def_dim(ncid,"time",NF_UNLIMITED,id_time)
        istat = nf_def_dim(ncid,"glon",imaxg,id_lon)
        istat = nf_def_dim(ncid,"glat",jmaxg,id_lat)
        istat = nf_def_dim(ncid,"mlon",zimxm,id_mln)
        istat = nf_def_dim(ncid,"mlat",zjmxm,id_mlt)

        ! define time varying scalars
        ids(1) = id_time

        istat = nf_def_var(ncid,"ftime",NF_DOUBLE,1,ids,id_ftim)
        write(char80,"('time in fraction day-of-year')")
        istat = nf_put_att_text(ncid,id_ftim,"long_name",
     |          len_trim(char80),trim(char80))

        istat = nf_def_var(ncid,"cuspnclt",NF_DOUBLE,1,ids,id_cuspnclt)
        write(char80,"('N cusp colatitude (radians)')")
        istat = nf_put_att_text(ncid,id_cuspnclt,"long_name",
     |          len_trim(char80),trim(char80))

        istat = nf_def_var(ncid,"cuspnlon",NF_DOUBLE,1,ids,id_cuspnlon)
        write(char80,"('N cusp longitude (radians)')")
        istat = nf_put_att_text(ncid,id_cuspnlon,"long_name",
     |          len_trim(char80),trim(char80))

        istat = nf_def_var(ncid,"cuspsclt",NF_DOUBLE,1,ids,id_cuspsclt)
        write(char80,"('S cusp colatitude (radians)')")
        istat = nf_put_att_text(ncid,id_cuspsclt,"long_name",
     |          len_trim(char80),trim(char80))

        istat = nf_def_var(ncid,"cuspslon",NF_DOUBLE,1,ids,id_cuspslon)
        write(char80,"('S cusp longitude (radians)')")
        istat = nf_put_att_text(ncid,id_cuspslon,"long_name",
     |          len_trim(char80),trim(char80))


        ! define time varying arrays
        ids(1) = id_lon
        ids(2) = id_lat
        ids(3) = id_time

        istat = nf_def_var(ncid,"eflx",NF_DOUBLE,3,ids,id_eflx)
        write(char80,"('AMIE electron flux (W/m^2)')")
        istat = nf_put_att_text(ncid,id_eflx,"long_name",
     |          len_trim(char80),trim(char80))

        istat = nf_def_var(ncid,"ekev",NF_DOUBLE,3,ids,id_ekev)
        write(char80,"('AMIE electron energy (keV)')")
        istat = nf_put_att_text(ncid,id_ekev,"long_name",
     |          len_trim(char80),trim(char80))


        ids(1) = id_mln
        ids(2) = id_mlt
        istat = nf_def_var(ncid,"poten",NF_DOUBLE,3,ids,id_pot)
        write(char80,"('AMIE potential (V)')")
        istat = nf_put_att_text(ncid,id_pot,"long_name",
     |          len_trim(char80),trim(char80))

        ! end define mode
        istat = nf_enddef(ncid)

        nc_outfile = ncid
        end


ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

        subroutine nc_add(ncid,index, ftime,cuspnclt,cuspnlon,cuspsclt,
     |             cuspslon, ekev, eflx, epot)
        implicit none
        ! ASPEN geographic grid
        integer, parameter :: imaxg = 73
        integer, parameter :: jmaxg = 36
        ! ASPEN magnetic grid
        integer, parameter :: zimxm = 80+1 
        integer, parameter :: zjmxm = 97
        ! arguments
        integer ncid, index
        real*8 cuspnclt, cuspnlon, cuspsclt, cuspslon, ftime
        real ekev(imaxg,jmaxg), eflx(imaxg,jmaxg), epot(zimxm,zjmxm)

        include "netcdf.inc"
        
        ! dimensions
        integer id_time, id_lon, id_lat
        ! variables
        integer id_pot, id_eflx, id_ekev, id_ftim
        integer id_cuspnclt, id_cuspnlon, id_cuspsclt, id_cuspslon
        common /netcdf/ id_time, id_lon, id_lat, id_ftim,id_pot,id_eflx,
     |     id_ekev, id_cuspnclt, id_cuspnlon, id_cuspsclt, id_cuspslon

        real*8 dekev(imaxg,jmaxg), deflx(imaxg,jmaxg),
     |         depot(zimxm,zjmxm)

        integer i,j,istat
        integer start3d(3), count3d(3)

        ! put time
        istat = nf_put_var1_double(ncid,id_ftim,index,ftime)

        ! put cusp info
        istat = nf_put_var1_double(ncid,id_cuspnclt,index,cuspnclt)
        istat = nf_put_var1_double(ncid,id_cuspnlon,index,cuspnlon)
        istat = nf_put_var1_double(ncid,id_cuspsclt,index,cuspsclt)
        istat = nf_put_var1_double(ncid,id_cuspslon,index,cuspslon)

        ! convert reals to doubles
        do i=1,imaxg
        do j=1,jmaxg
           dekev(i,j) = ekev(i,j)
           deflx(i,j) = eflx(i,j)
        enddo
        enddo
        do i=1,zimxm
        do j=1,zjmxm
           depot(i,j) = epot(i,j)
        enddo
        enddo

        ! put array info
        start3d(1) = 1
        start3d(2) = 1
        start3d(3) = index
        count3d(1) = imaxg
        count3d(2) = jmaxg
        count3d(3) = 1

        istat = nf_put_vara_double(ncid,id_ekev,start3d,count3d,dekev)
        istat = nf_put_vara_double(ncid,id_eflx,start3d,count3d,deflx)
        count3d(1) = zimxm
        count3d(2) = zjmxm
        istat = nf_put_vara_double(ncid,id_pot,start3d,count3d,depot)

        end

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

        real*8 function fday(iyear,imon,ida,ihr,imin)
        implicit none
        ! arguments
        integer iyear, imon, ida, ihr, imin
        ! local vars
        integer i
        integer days(12)
        data days /31,28,31,30,31,30,31,31,30,31,30,31/

        ! add months
        fday = 0.
        do i=1,imon-1
          fday=fday + days(i)
        enddo

        ! add leap year day (doesn't work for year 1000, 3000, etc.)
        if ((mod(iyear,4) .eq. 0) .and. (imon .gt. 2)) fday = fday+1

        ! add day, hr, min
        fday = fday + ida + (ihr + imin/60.d0)/24.d0

        return
        end
        
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
C
      SUBROUTINE defmag(imaxm,jmaxm,zlatm,zlonm)
      implicit none
C     ****
C     ****     PARAMETERS DEFINE ASPEN GEOMAGNETIC GRIDS
C     ****     adapted from ASPEN's CNSTNT subroutine
C     ****
      !arguments
      integer :: imaxm,jmaxm
      real :: zlatm(jmaxm),zlonm(imaxm)
! Local:
      real,parameter :: alfa=1.668
      real :: TANTH0(JMAXM),TANTHS(JMAXM),
     1  THETA0(JMAXM),HAMH0(JMAXM)
      real :: r1,pi,dlatm,dlonm
      integer :: i,j,n
!
! Init:
      r1 = 1.6387955E-2
      PI = 4.*ATAN(1.)
      DLATM = PI/FLOAT(JMAXM-1)
      DLONM = 2.*PI/FLOAT(IMAXM-1)
C     ****
C     *+**     FILL ARRAY YLATM (EQUALLY SPACED IN THETAO BUT HOLDS
C     ****       CORRESPONDING VALUE OF THETAS)
C     ****
      DO 3 J = 1,JMAXM
C       ****    THETA0 = EQUALLY SPACED GRID VALUES
        THETA0(J) = -PI/2.+FLOAT(J-1)*DLATM
    3 CONTINUE
      DO 4 J = 2,JMAXM-1
C       ****     TANTH0 = ABS(TAN(THETA0))
        TANTH0(J) = ABS(TAN(THETA0(J)))
C       ****     HAMH0 = HA-H0
        HAMH0(J) = R1*TANTH0(J)+TANTH0(J)**(2.+2.*ALFA)/
     1    (1.+TANTH0(J)**2)**ALFA
C       ****     TANTHS = ABS(TAN(THETAS))
        TANTHS(J) = SQRT(HAMH0(J))
C       ****     YLATM = TANTHS
        ZLATM(J) = SIGN(ATAN(TANTHS(J)),THETA0(J))
    4 CONTINUE
C     ****
C     ****     NOW DO POLES
C     ****
      ZLATM(1) = THETA0(1)
      ZLATM(JMAXM) = THETA0(JMAXM)
C     ****
C     ****     FILL YLONM
C     ****
      DO 5 I = 1,IMAXM
        ZLONM(I) = -PI+FLOAT(I-1)*DLONM
    5 CONTINUE
      END
C
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine aterp(imaxm,jmaxm,ylatm,EPOTM,zimxm,zjmxm,zlatm,EPOTZ)
      ! interpolation from one magnetic grid to another
      implicit none
      ! arguments
      integer :: imaxm,jmaxm,zimxm,zjmxm
      real :: ylatm(jmaxm),zlatm(zjmxm)
      real :: EPOTM(imaxm,jmaxm), EPOTZ(zimxm,zjmxm)
      ! local vars
      integer :: i,j,ii,jj
      integer :: i1(zimxm),i2(zimxm), j1(zjmxm),j2(zjmxm)
      real :: ylonm(imaxm), zlonm(zimxm)
      real :: fi(zimxm), fj(zjmxm)
      real :: pi

      ! begin execution

      ! set up longitudes
      pi = 4.* atan(1.)
      do i=1,imaxm
         ylonm(i) = -pi + float(i-1)*(2.*pi/(imaxm-1))
      enddo
      do i=1,zimxm
         zlonm(i) = -pi + float(i-1)*(2.*pi/(zimxm-1))
      enddo

      ! set up interpolation parameters - longitude
      i1(1) = 1
      i2(1) = 2
      fi(1) = 1.0
      i1(zimxm) = imaxm
      i2(zimxm) = 1
      fi(zimxm) = 1.0
      do i=2,zimxm-1
         ii=2
         do while (zlonm(i) > ylonm(ii))
            ii = ii+1
         end do
         i1(i) = ii-1
         i2(i) = ii
         fi(i) = (ylonm(ii) - zlonm(i))/(ylonm(ii) - ylonm(ii-1))
      enddo

      ! set up interpolation parameters - latitude
      ! ylatm array is decreasing, zlatm array is increasing
      j1(1) = jmaxm
      j2(1) = jmaxm-1
      fj(1) = 1.0
      j1(zjmxm) = 1
      j2(zjmxm) = 2
      fj(zjmxm) = 1.0
      do j=2,zjmxm-1
         jj=jmaxm
         do while (zlatm(j) > ylatm(jj))
            jj = jj-1
         end do
         j1(j) = jj+1
         j2(j) = jj
         fj(j) = (ylatm(jj) - zlatm(j))/(ylatm(jj) - ylatm(jj+1))
      enddo

      ! do interpolation
      do i=1,zimxm
      do j=1,zjmxm
         EPOTZ(i,j) = fi(i)*     fj(j) *EPOTM(i1(i),j1(j))
     |              + fi(i)*(1.0-fj(j))*EPOTM(i1(i),j2(j))
     |              + (1.0-fi(i))*fj(j)*EPOTM(i2(i),j1(j))
     |              + (1.0-fi(i))*(1.0-fj(j))*EPOTM(i2(i),j2(j))
      enddo
      enddo
      

      end subroutine
