#include "dims.h"
!
      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"
#ifdef OSF1
!
! Floating point exception flags for osf compaq:
! If DIV0+INV, will abort on div by zero or invalid, and will flush
!   overflows to "huge" and underflows to zero.
!
#include "for_fpe_flags.f"
      integer :: new_fpe_flags = FPE_M_TRAP_DIV0 + FPE_M_TRAP_INV
      integer :: old_fpe_flags
      integer,external :: for_set_fpe 
#endif
!
! 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
      character (len=24) :: host
      integer :: val1(8),i,ier
!
! Externals:
#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)
      call gethostsname(host)
      write(6,"('Host = ',a)") trim(host)
!
#ifdef OSF1
!
! Set fpe trap flags on compaq:
!     old_fpe_flags = for_set_fpe(new_fpe_flags)
!     write(6,"('Set OSF fpe flags.')")
#endif
!
! 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 and derived constants:
      call init
!
! 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
!
! Set up magnetic coordinate grid and related parameters:
! 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 gravity wave parameterization:
! (gravity waves are not in tiegcm)
      call mgwinti
!
! 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) (ibm rtc function) = ',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
      write(6,"('Average wall clock seconds per time step =',f7.3)")
     |  wall_time/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) :: nsec
      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
        nsec = dv(7)+60*dv(6)+60*60*dv(5) ! secs
      end if
      t = real(nsec)
      end subroutine wtimer