!
!-----------------------------------------------------------------------
! Purpose:
! This is a test driver for the standalone edynamo module.
! Neutral atmosphere inputs are read from external model history files,
!   e.g., WACCM or TIMEGCM.
!
!-----------------------------------------------------------------------
!
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
  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
  use heelis      ,only: heelis_model
  use edynamo     ,only: dynamo
  use output      ,only: output_init,write_output
  use read_ncfile ,only:               &
    get_geogrid,read_waccm,read_tgcm,  & ! subroutines
    ntime,iday,year,ut,                & ! time and date
    tn,un,vn,omega,zpot,barm,ped,hall    ! 3d input fields
  use timing      ,only: report_timing,time_total
  implicit none
#include <mpif.h>
  
  integer :: itime
  real(r8) :: secs
  real(r8) :: steptime=0._r8 ! wc time per step
  real(r8) :: &               ! for local timing
    begtotal,endtotal,begstep,endstep
!
! Begin execution:
!
  call set_cons            ! set runtime constants
  call mp_init             ! initialize mpi library 
  begtotal = mpi_wtime()   ! begin wallclock timing
  call get_namelist(mytid) ! get namelist input from user
!
! Read dimensions and coordinate vars from input file, 
! and set global 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)
!
  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 edynamo inputs from neutral atmosphere model:
!
    if (trim(model_name) == 'TIMEGCM') then
      call read_tgcm(input_ncfile,itime)
    elseif (trim(model_name) == 'WACCM') then
      call read_waccm(input_ncfile,itime)
    endif
!
! 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)
!
! Write to netcdf output file:
! (sub write_output is in output.F90)
!
    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 
!-----------------------------------------------------------------------