!
!-----------------------------------------------------------------------
! 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
!  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
  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
!  call mp_init             ! initialize mpi library 
!  begtotal = mpi_wtime()   ! begin wallclock timing

  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) 
!
! Set parameter-based geomagnetic grid:
!
!  call set_maggrid         ! define parameter-based magnetic grid
!
! Decompose geographic and magnetic grids:
!
!  call mp_decomp(trim(model_name))  ! set up domain decomposition
!
! 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

!    begstep = mpi_wtime() ! start step-timing

!   write(iulog,"(/,'Main: begin step ',i4,' out of ',i4)") itime,ntime
!
! Get inputs from TIME-GCM:
!
!    if (trim(model_name) == 'TIMEGCM') then
!
  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.
!
    do ilev = 1, nlev

      if (pMid(1,ilev) >= 0.00001) oPBot = ilev - 1
      if (pMid(1,ilev) >= 0.5) ionBot = ilev - 1

    enddo
    
    write(*,*) 'oPBot, ionBot ', oPBot, ionBot
!
! Get apex coordinates and related parameters for current year:
! (sub get_apex is in apex.F90) (called once per run)
!
!    if (itime==1) call get_apex(year,nlon,nlat,jspole,jnpole)
!
! Get high-latitude potential from Heelis empirical model:
! (sub heelis_model is in heelis.F90)
!
!    secs = ut*3600._r8 ! ut in seconds
!    call heelis_model(iday,secs)
!
! Call main edynamo driver:
! (sub dynamo is in edynamo.F90)
!
!    call dynamo(tn,un,vn,omega,zpot,barm,ped,hall, &
!      lev0,lev1,lon0,lon1,lat0,lat1)
!
! Call WACCM-X extended ionosphere interface method:
! (sub dynamo is in edynamo.F90)
!
  write(*,*) 'Calling ionos_intr'
    call ionos_intr(1,oPBot,ionBot)
!
! 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 

!  endtotal = mpi_wtime()        ! end wctiming
!  time_total=endtotal-begtotal
!
!  call report_timing            ! report timing to stdout
!
!  call mp_close                 ! finalize MPI

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

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