!
      program tgcm
!
! Main tgcm model program:
!
      use input_module,only: input
      use init_module,only: init,iyear
      use bndry_module,only: lowbound
      use apex_module,only: apxparm
      use hist_module,only: nstep
#ifdef MPI
      use mpi_module,only: mp_init,mp_close,mytid,ntask
#endif
      implicit none
      include "params.h"
      include "fgcom.h" ! for lat0,lat1
!
! Local:
      real ::
     | total_cpu,    ! from f95 intrinsic cpu_time()
     | etime,        ! per processor elapsed time (AIX only)
     | total_etime,  ! total elapsed time for all processors (AIX only)
     | rtc_secs,     ! real time clock from aix extension rtc()
     | wall_time     ! wall clock time from date_and_time intrinsic
      character (len=8) :: date
      character (len=10) :: time
      character (len=5) :: zone
      integer :: val1(8),i,ier,istat
#ifdef AIX
      real,external :: timef,global_timef,rtc ! aix extensions
#endif
      call date_and_time(date,time,zone,val1)
      write(6,"('Begin execution of ',a,' at mm/dd/yy ',a,'/',a,'/',a,
     |  '  hh:mm:ss ',a,':',a,':',a)") tgcm_version,
     |  date(5:6),date(7:8),date(3:4),
     |  time(1:2),time(3:4),time(5:6)
!
! Initialize timing structures:
!
#ifdef AIX
      etime = timef()
      total_etime = global_timef()
      rtc_secs = rtc() 
#endif
      call wtimer(wall_time)
      call cpu_time(total_cpu) ! f95 intrinsic
!
! Set lat0,1 (fgcom.h) for non-mpi job. 
! For mpi job, these will be set in mp_init below.
      lat0 = 1
      lat1 = zjmx
!
! If MPI job, init mpi and distribute latitude indices 
! across tasks.
!
#ifdef MPI
      call mp_init
      call input(mytid) ! get user input (mpi)
#else
      call input(-1)    ! get user input (non-mpi)
#endif
!
! ncarg for debug (not on all machines)
!     call opngks
!
! Initialize variables:
      call init
!
! Set constants:
      call con
!
! Dynamo constants:
      call cnstnt
!
! Read source history. First attempt to read source as a
! netcdf file (readsource in rdsource.f). If this fails, 
! try reading as an old cray-blocked history (read_oldsrc.F),
! but only if running under UNICOS (i.e., cannot read old 
! cray-blocked histories on non-unicos platforms because we 
! do not want to be dependent on ncaru (crayopen, crayread, etc)).
! 
      call readsource(ier)
      if (ier > 0) then
#ifdef UNICOS
        write(6,"(/,'>>> tgcm: error return from readsource')")
        write(6,"(10x,'Will try reading as an old cray-blocked',
     |    ' file..')")
        call read_oldsrc
#else
        write(6,"(/,'>>> tgcm: error return from readsource')")
        write(6,"('  This may be an old cray-blocked history',
     |    ' file. If so, please try running on a unicos machine.')")
        stop 'source'
#endif
      endif
!
! Read magnetic data file:
#ifdef MPI
      call rdmag(mytid)
#else
      call rdmag(-1)
#endif
!
! apxparm obtained from Art 2/25/00. This is first step toward
!   new dynamo.
! 4/11/00: this is now working -- can run old dynamo with either
!   apxparm or rdmag.
!
!     write(6,"('tgcm call apxparm: iyear=',i5,'...')") iyear
!     call apxparm(real(iyear))
!     call magdyn
!     call dynpotmag
!
! Set lower boundaries:
      call lowbound
!
! Init sigmas, segeuv, etc (see qrj.f):
      call init_sigmas
!
! Advance the model in time:
      call advnce
!     call tstwrhist ! call instead of advnce for testing, debug
!
! ncarg for debug (not on all machines)
!     call clsgks
!
! If MPI job, finalize mpi and report timing:
!
#ifdef MPI
      call mp_close
      write(6,"('MPI run with ntask = ',i3)") ntask
#endif
!
! timef, global_timef, and rtc() are AIX extensions.
! All 3 appear to report the same time. Levesque says
! we should be using rtc. 
!
#ifdef AIX
!     etime = timef()               ! elapsed time for this proc (ms)
!     total_etime = global_timef()  ! total elapsed time all procs (ms)
!     write(6,"('Task  elapsed time (secs) = ',f15.3,' (',f5.2,
!    |  ' hours)')") etime/1000.,etime/(1000.*3600.)
!     write(6,"('Total elapsed time (secs) = ',f15.3,' (',f5.2,
!    |  ' hours)')") total_etime/1000.,total_etime/(1000.*3600.)
      rtc_secs = rtc()-rtc_secs
      write(6,"('Real clock time (secs) = ',f15.3)") rtc_secs
#endif
!
! Report total cpu time (f95 intrinsic cpu_time())
! (elapsed time is reported at the script level from timex)
!
      call wtimer(wall_time)
      write(6,"('Wall clock time (secs) = ',f15.3)") wall_time/1000.
      write(6,"('Average wall clock seconds per time step =',f7.3)")
     |  wall_time/1000./float(nstep)
      call cpu_time(total_cpu)
      write(6,"('Total cpu time (secs)  = ',f15.3)") total_cpu
      write(6,"('NORMAL EXIT')")
      end program tgcm
!-----------------------------------------------------------------------
      subroutine wtimer(t)
!
! Use fortran intrinsic date_and_time to return elapsed wall clock 
! time t since this routine was last called.
!
      implicit none
      real,intent(out) :: t
      character (len=8) :: date
      character (len=10) :: time
      character (len=5) :: zone
      integer :: val1(8), val2(8), dv(8), its
      integer(kind=8) :: ms
      data its / 0 /
      save
! 
      if ( its == 0 ) then
        call date_and_time(date,time,zone,val1)
        t = 0.0
        its = 1
      else
        call date_and_time(date,time,zone,val2)
        dv = val2-val1
        ms = dv(8)+1000*dv(7)+60*1000*dv(6)+60*60*1000*dv(5)
      end if
      t = real(ms)
      end subroutine wtimer