!
!-----------------------------------------------------------------------
! Purpose:
! This is a test driver for the standalone single column WACCM-X 
! ionospheric subroutines (based on Ben's edynamo driver)
!
!-----------------------------------------------------------------------
!
program main 
  use shr_kind_mod,only: r8 => shr_kind_r8 ! 8-byte reals
  use params      ,only: set_cons,iulog
  use geogrid     ,only: nlon,nlat,jspole,jnpole,nlev,zlev,pMid,glat,glon
!  use maggrid     ,only: set_maggrid
!  use apex        ,only: get_apex   ! get mag apex coordinates
!  use my_mpi      ,only: mytid,mp_init,mp_decomp,mp_close,lev0,lev1,lon0,lon1,lat0,lat1
  use namelist    ,only: get_namelist,model_name,input_ncfile,output_ncfile,inLat,inLon
!  use heelis      ,only: heelis_model
!  use edynamo     ,only: dynamo
  use output      ,only: output_init,write_output
  use read_ncfile ,only:                  &
    get_geogrid,read_tgcm_wx,  & ! subroutines
    ntime,iday,year,ut,                   & ! time and date
    p0,                           & ! other input fields
    tn,un,vn,omega                  ! 3d input fields
!  use timing      ,only: report_timing,time_total
  use ionosphere,  only : ionos_intr

  implicit none
!#include <mpif.h>
  
  integer :: itime, ilev, oPBot, ionBot, iLt, iLn
  real(r8) :: secs
  real(r8) :: steptime=0._r8 ! wc time per step
  real(r8) :: &               ! for local timing
    begtotal,endtotal,begstep,endstep
!  real(r8) :: pres(97)
!  real(r8),allocatable,dimension(:),save :: pres
!
! Begin execution:
!
  call set_cons            ! set runtime constants

  write(*,*) 'Reading namelist'

  call get_namelist() ! get namelist input from user
!
! Read dimensions and coordinate vars from input file, 
! and set global geographic grid.
!
  write(*,*) 'Getting geographic grid'

  call get_geogrid(input_ncfile,model_name) 

! Initialize output file:
! (must be called once per run, before write_output)
!
  write(*,*) 'Initializing output'

  call output_init ! initialize output 
!
! Loop over number of histories on input file:
!
!  do itime=1,ntime   ! process all times on input file
 do itime=1,1       ! process first time only

! do iLt = 1, nlat
!   inLat = glat(iLt)
!   do iLn = 1, nlon
!     inLon = glon(iLn)
     
     write(*,*) 'Doing lat/lon ',inLat, inLon


!
  write(*,*) 'Reading TIME-GCM input file'

      call read_tgcm_wx(input_ncfile,itime)
!    elseif (trim(model_name) == 'WACCM') then
!      call read_waccm(input_ncfile,itime)
!    endif
! 
! Check pressure to find bottom index for O+ and Te calculations
! ionBot and oPBot are bottom vertical levels to calculate Te/Ti and O+/e respectively. 
! ionBot should correspond to about 60 km or 0.5 hPa. 
! oPBot should correspond to about 130 km or 0.00001 hPa.
! oPBot should correspond to about 100 km or 0.0001 hPa.
!
    do ilev = 1, nlev

!      if (pMid(1,ilev) <= 0.001) oPBot = ilev + 1
      if (pMid(1,ilev) <= 0.01) oPBot = ilev + 1
      if (pMid(1,ilev) <= 50.) ionBot = ilev + 1

    enddo
    
    write(*,*) 'oPBot, ionBot ', oPBot, ionBot

    write(*,*) 'Calling ionos_intr'
    call ionos_intr(1,oPBot,ionBot)

!enddo !lon
!enddo !lat
!
! Write to netcdf output file:
! (sub write_output is in output.F90)
!
   write(*,*) 'Writing output'

   call write_output(output_ncfile,itime)

!    endstep = mpi_wtime() ! end step-timing
!    steptime = endstep-begstep

!    write(iulog,"('Step ',i4,' of ',i4,' year=',i4,' day=',i4,' ut=',f6.2,' (WC=',f10.2,' seconds)')") &
!      itime,ntime,int(year),iday,ut,steptime

  enddo ! itime=1,ntime 

  write(iulog,"('NORMAL EXIT')")

end program main 
!-----------------------------------------------------------------------