subroutine initialize_amie

  use ModAMIE
  use ModFiles
  use ModConstants

  implicit none
  integer :: ierr

  integer, dimension(7) :: itime
  real                  :: rtime
  integer, external     :: jday

  real    :: blatfit, j10, blatf10, tdeg
  integer :: ith, i, j, k, l, m, n, lon, LunHeader

  integer :: in_ithmx, in_mrow, in_kmx, in_mmx, in_nq, in_NQS, in_MXNMX
  integer :: in_mcoef, in_ithtrns, ir

  character (len=100) :: temp, StrData, StrPE

  iteration_number = 0

  !\
  ! Get location of magnetic pole :
  !/

  call time_real_to_int(currenttime, itime)

  rtime = float(itime(1)) + &
       float(jday(itime(1),itime(2),itime(3))) / &
       float(jday(itime(1),      12,      31))

  call cofrm(rtime)
  call dypol(pole_colat, pole_east_lon, vp)

  call apxrda(6, apex_file, lun_apex, rtime, wkap, apex_parameter, ierr)

  write(temp,"(I5)") mcoef*100 + ithtrns
  call merge_str(qset_file, temp)
  qset_file = temp

  if (DebugLevel > 2) &
       write(*,*) "  ===> read in qset file, me: ",qset_file, me_world

  open (lun_qset,file=qset_file,status='old', iostat = ierr)

  if (ierr /= 0) then
     write(6,*) 'Error opening file :',qset_file
     stop
  endif

  read (lun_qset,"(e14.4,2i4)") dth,in_mcoef,in_ithtrns

  read (lun_qset,"(8i4)") in_ithmx,in_mrow,in_kmx,in_mmx,in_nq, &
       in_NQS ,in_MXNMX ,in_ITHTRNS

  do i=0,ithmx
     read (lun_qset,"(i3,3f8.3)") ir,st(i),ct(i),tdeg
     st(i) = sin(tdeg*dtor)
     ct(i) = cos(tdeg*dtor)
     if (i .eq. 0) st(0) = 1.e-10
     read (lun_qset,"(10e12.3)") (q(j,i),j=1,nq)
     read (lun_qset,"(10e12.3)") (dq(j,i),j=1,nq)
  enddo

  do i=0,mmx
     read (lun_qset,"(4i6)") ir,ns(i),nmx(i),nss(i)
  enddo

  do i=1,nq
     read (lun_qset,"(i3,3e14.4)") ir,wnt2(i),wnp2(i),qnorm(i)
  enddo

  do i=1,mxnmx
     read (lun_qset,"(10e12.3)") (vfmps(j,i),j=1,nqs)
     read (lun_qset,"(10e12.3)") (zfmps(j,i),j=1,nqs)
  enddo

  blatfit = 90. * (1. - float(ithtrns)/float(ithmx))
  j10 = (blatfit + 4.9999) / 10.
  blatf10 = j10 * 10

  wnt2s(1) = 0.
  wnp2s(1) = 0.

!  qsclr is defined to ithmx, but is only reasonable to ~ithtrns+2

  do ith=0,ithmx
     qsclr(1,ith) = 1.
  enddo
  
  do n=2,nmx(0)
     wnt2s(n) = wnt2(2*n-3)
     wnp2s(n) = 0.
     do ith=0,ithmx
        qsclr(n,ith) = q(2*n-3,ith)*qnorm(2*n-3)
     enddo
  enddo

  do m=1,mmx
     i = nss(m) - 1
     do n=1,nmx(m)
        i = i + 1
        j = ns(m) + 2*(n-1)
        wnt2s(i) = wnt2(j)
        wnp2s(i) = wnp2(j)
        do ith=0,ithmx
           qsclr(i,ith) = q(j,ith)*qnorm(j)
        enddo
     enddo
  enddo

  xlonmx=lonmx
  dlon = twopi/float(lonmx)
  dth  = pi/2.0/float(ithmx)
  mlt(0) = 0.

  do lon=1,lonmx
     mlt(lon) = lon*24/float(lonmx)
     sl(lon) = sin(lon*dlon)
     cl(lon) = cos(lon*dlon)
  enddo

  do ith=0,ithmx
     clatd(ith) = ith*dth*rtod
  enddo

  ibm(-mmx) = 1
  i = 0

  do m=-mmx,mmx
     do n=1,nmx(iabs(m))
        i = i + 1
        mi(i) = m
     enddo
     iem(m) = ibm(m) - 1 + nmx(iabs(m))
     if (m.ne.mmx) ibm(m+1) = nmx(iabs(m)) + ibm(m)
  enddo

  imx = iem(mmx)
  call fcmp(mmx,lonmx,cl,sl,amie_f)

  !\
  ! Open output files
  !/

  if (numprocs >= 2) then 
     if (me_world < 10) then
        write(StrPE,'(a3,I1)') '_00',me_world
     elseif (me_world < 100) then
        write(StrPE,'(a2,I2)') '_0',me_world
     elseif (me_world < 1000) then
        write(StrPE,'(a1,I3)') '_',me_world
     else
        write(*,*) 'Yikes. You have a lot of processors.'
        write(*,*) 'Why are you running AMIE on so many?'
        call stop_mpi('reduce the number of processors to < 1000')
     endif
  else
     StrPE = ' '
  endif
  temp = StrPE
  call merge_str(output_file,temp)

  if (IsOutputBinary) then
     open(LunMainOutput,file=temp,status='unknown',form='UNFORMATTED', &
          iostat = ierr)
  else
     open(LunMainOutput,file=temp,status='unknown',iostat = ierr)
  endif

  if (ierr.lt.0) then
     write(6,*) 'Error opening file :',temp
     call stop_mpi('')
  endif

  if (IsOutputBinary) then

     write(LunMainOutput) ithtrns+1,lonmx+1, ntimes
     write(LunMainOutput) (clatd(i),i=0,ithtrns)
     write(LunMainOutput) mlt
     write(LunMainOutput) 18
     write(LunMainOutput) 'Electric Potential            '
     write(LunMainOutput) 'Pedersen Conductance (Total)  '
     write(LunMainOutput) 'Pedersen Conductance (Aurora) '
     write(LunMainOutput) 'Hall Conductance (Total)      '
     write(LunMainOutput) 'Hall Conductance (Aurora)     '
     write(LunMainOutput) 'Auroral Mean Energy (keV)     '
     write(LunMainOutput) 'Auroral Energy Flux (W/m2)    '
     write(LunMainOutput) 'Electric Field (East)         '
     write(LunMainOutput) 'Electric Field (North)        '
     write(LunMainOutput) 'Horizontal Current (East)     '
     write(LunMainOutput) 'Horizontal Current (North)    '
     write(LunMainOutput) 'Horizontal Current (Magnitude)'
     write(LunMainOutput) 'Simple Joule Heating          '
     write(LunMainOutput) 'Auroral Mean Energy (model)   '
     write(LunMainOutput) 'Auroral Energy Flux (model)   '
     !  write(LunMainOutput) 'Field Aligned Current         '
     !  write(LunMainOutput) 'Electric Field Error          '
     write(LunMainOutput) 'Pedersen Conductance Model    '
     write(LunMainOutput) 'Hall Conductance Model        '
     write(LunMainOutput) 'Potential Model               '

  else

     write(LunMainOutput,'(3I4)') ithtrns+1,lonmx+1, ntimes
     write(LunMainOutput,*) (clatd(i),i=0,ithtrns)
     write(LunMainOutput,*) mlt
     write(LunMainOutput,*) 18
     write(LunMainOutput,'(a)') 'Electric Potential            '
     write(LunMainOutput,'(a)') 'Pedersen Conductance (Total)  '
     write(LunMainOutput,'(a)') 'Pedersen Conductance (Aurora) '
     write(LunMainOutput,'(a)') 'Hall Conductance (Total)      '
     write(LunMainOutput,'(a)') 'Hall Conductance (Aurora)     '
     write(LunMainOutput,'(a)') 'Auroral Mean Energy (keV)     '
     write(LunMainOutput,'(a)') 'Auroral Energy Flux (W/m2)    '
     write(LunMainOutput,'(a)') 'Electric Field (East)         '
     write(LunMainOutput,'(a)') 'Electric Field (North)        '
     write(LunMainOutput,'(a)') 'Horizontal Current (East)     '
     write(LunMainOutput,'(a)') 'Horizontal Current (North)    '
     write(LunMainOutput,'(a)') 'Horizontal Current (Magnitude)'
     write(LunMainOutput,'(a)') 'Simple Joule Heating          '
     !  write(LunMainOutput,'(a)') 'Field Aligned Current         '
     !  write(LunMainOutput,'(a)') 'Electric Field Error          '
     write(LunMainOutput,'(a)') 'Auroral Mean Energy (model)   '
     write(LunMainOutput,'(a)') 'Auroral Energy Flux (model)   '
     write(LunMainOutput,'(a)') 'Pedersen Conductance Model    '
     write(LunMainOutput,'(a)') 'Hall Conductance Model        '
     write(LunMainOutput,'(a)') 'Potential Model               '

  endif

  if (me_world == 0) then

     LunHeader = LunMainOutput + 1
     temp = '.h'
     call merge_str(output_file,temp)
     open(LunHeader,file=temp,status='unknown', iostat = ierr)
     write(LunHeader,*) numprocs

     do i=0,numprocs-1

        if (numprocs >= 2) then
           if (i < 10) then
              write(temp,'(a3,I1)') '_00',i
           elseif (i < 100) then
              write(temp,'(a2,I2)') '_0',i
           elseif (i < 1000) then
              write(temp,'(a1,I3)') '_',i
           endif
        else
           temp = ' '
        endif

        call merge_str(output_file,temp)

        call strlen(temp,j)
        write(LunHeader,'(a)') temp(1:j)

     enddo

     temp = ' '
     call merge_str(output_file,temp)

     call strlen(temp,j)
     write(LunHeader,'(a)') temp(1:j)

     close(LunHeader)

  endif

  if (OutputData) then

     StrData = "_data"
     temp = StrPE
     call merge_str(StrData,temp)
     call merge_str(output_file,temp)

     open(lun_dataoutput,file=temp,status='unknown', iostat = ierr)

     if (ierr.lt.0) then
        write(6,*) 'Error opening file :',temp
        call stop_mpi('')
     endif

  endif

end subroutine initialize_amie