#include "dims.h"
!
      subroutine advnce
!
! Advance model from START to STOP:
!
      use input_module,only: calday,start,step,f107,f107a,power,ctpoten
      use init_module,only: secs,istep,iter,iyear,iday,uthr,igetgpi
      use hist_module,only: nstep,modeltime
      use gpi_module,only: rdgpi,getgpi
      use cons_module,only: jmax,dt,cs,re
      use timing_module,only: end_timing
#ifdef MPI
      use mpi_module,only: mytid,mp_close,mp_snddyn,mp_rcvdyn,
     |  mp_fgtomaster,mp_bndlats,ntask,lats_have,mp_updatephi,
     |  mp_sechtomaster
      implicit none
#include "mpif.h"
#else
      implicit none
#endif
#include "params.h"
#include "buff.h"
#include "phys.h"
#include "amie.h"
#include "fgcom.h"
!
! Local:
      integer :: m,k,i,jj,itemp,jx,iprint,ier,num_threads,
     |  iseries_p,iseries_s,jmaxv
      integer(kind=8) :: nsecs ! current model time in seconds,
                               ! including model day
      real :: sec
      real :: cputime0,cputime1,wall_time,secpstep
      logical :: time2write,wrprim,wrsech,save_prim,save_sech,
     |  newseries_sech,newseries_prim
!
! External:
      real,external :: cputime
      logical,external :: wrhist
      integer,external :: omp_get_num_threads
!
! Init timestep index, modeltime, and nsecs:
      istep = 0
      modeltime(1:3) = start(:,1)
      modeltime(4) = 0
      call modeltime_to_nsecs(modeltime,nsecs) ! sets nsecs, includes day
!
! Init buffer indices:
      modulo=8
      do m=1,modulo
        njtmp(m+1) = nbaddr(m)
      enddo
!
! Echo source history to primary history output if an initial run,
! and write initial secondary history if necessary (istep==0).
! (if mpi, only master writes histories)
!
#ifdef MPI
      if (mytid==0) call outhist(istep,modeltime)
#else
      call outhist(istep,modeltime)
#endif
!
      ixtimec = ixtimep+1 ! reset after reading/writing source hist
!
      cputime0 = cputime()
      write(6,"(/'Enter advance at iter ',i6,' total accum cpuTime=',
     |  f8.2)") iter,cputime0
!
! Main time loop:
  100 continue
      iter=iter+1
      istep = istep+1
!
! Increment current model time (day,hr,min,secs) by one time step:
!
      nsecs = nsecs+step
      call nsecs_to_modeltime(nsecs,modeltime) ! increments modeltime(4)
      secs = mod(iter*int(dt),86400)
      uthr = secs/3600.
!
! Update time, calendar day, sun's declination, amie, gpi:
      if (calday > 0) call advnceday
      sec = float(iter)*dt
      sectgcm = float(iter)*dt
      call sunloc
      if (iamie==1) call setamie(sec)
!
! Get gpi data. Data file is opened and read by rdgpi before the 
! first time step, defining module data in gpi_mod.f. This data is
! subsequently referenced at every time step by getgpi. If mpi job, 
! only the master task obtains the data file from mss.
!
      if (igetgpi > 0) then
        iprint = 0
        if (istep==1) then
#ifdef MPI
          call rdgpi(mytid)
#else
          call rdgpi(-1)
#endif
          iprint = 1
        endif
        call getgpi(iyear,iday,int(secs),iprint)
      endif
!
! Auroral parameters:
      iprint = 0
      if (istep==1) iprint = 1
      call tail(iprint)
!
! Report this time step to stdout:
! If doing <= 100 steps, print at every step, otherwise print every 
!   10th step. If getting GPI from database, print f107,power,ctpoten.
!
      call end_timing(wall_time)
      secpstep = wall_time/float(istep)
!     if (nstep <= 100) then
!       if (igetgpi <= 0) then
          write(6,"('Step ',i5,' of ',i5,' mtime ',i3,',',i2,',',i2,
     |      ',',i2,' wcsecs/step=',f6.2)") istep,nstep,modeltime,
     |      secpstep
!       else
!         write(6,"('Step ',i5,' of ',i5,' mtime ',i3,',',i2,',',i2,
!    |      ',',i2,' f107d=',f6.2,' power=',f6.2,' ctpoten=',
!    |      f6.2,' wcsecs/step=',f6.2)") istep,nstep,modeltime,f107,
!    |      power,ctpoten,secpstep
!       endif
!     else
!       if (mod(istep,10)==0) then
!         if (igetgpi <= 0) then
!           write(6,"('Step ',i5,' of ',i5,' mtime ',i3,',',i2,',',
!    |        i2,',',i2,' wcsecs/step=',f6.2)") istep,nstep,modeltime,
!    |        secpstep
!         else
!           write(6,"('Step ',i5,' of ',i5,' mtime ',i3,',',i2,',',
!    |        i2,',',i2,' f107d=',f6.2,' power=',f6.2,' ctpoten=',
!    |        f6.2,' wcsecs/step=',f6.2)") istep,nstep,modeltime,f107,
!    |        power,ctpoten,secpstep
!         endif
!       endif
!     endif
!
! init_sflux calls ssflux and euvac (see qrj.f):
      call init_sflux
!
! Calculate 3d electric field:
      call vdrift
!
! Calculate global NVC, NMS, and NZ at lats jj and store in fg:
! addiag: WRITE fg(nvc,nms,nz) <- READ fg(nv,nt,nps,nps2)
!
#ifdef MPI
!
      if (istep==1.or.ntask==1) then
        do jj=1,jmax
          call addiag(jj,ixtimep,0)
        enddo
      else
        do jj=lat0,lat1
          call addiag(jj,ixtimep,0)
        enddo
!
! Need to update bndry lats of nvc,nms, and nz
! (messages could be shortened to update only nvc,nmc,nz 
!  instead of all fields)
!
        call mp_bndlats(ixtimep)
      endif
!
#else
!
! Some compilers require that objects in openmp clauses must
! be variables (not parameter constants), hence use of local jmaxv.
      jmaxv = jmax
!
!$OMP PARALLEL DO PRIVATE(jj) SHARED(jmaxv,ixtimep)
!MIC$ DOALL PRIVATE(jj) SHARED(jmaxv,ixtimep)
      do jj=1,jmaxv
        call addiag(jj,ixtimep,0)
      enddo
#endif
!
! Make boundary latitudes:
! mklatbndfg: 
!   WRITE fg(j= 0) <- READ fg(j=1)
!   WRITE fg(j=-1) <- READ fg(j=2)
!   WRITE fg(j=jmax+1) <- READ fg(j=jmax)
!   WRITE fg(j=jmax+2) <- READ fg(j=jmax-1)
!
#ifdef MPI
      if (istep==1.or.ntask==1) then
        call mklatbndfg(0     ,1     ,ixtimep,.false.)
        call mklatbndfg(-1    ,2     ,ixtimep,.false.)
        call mklatbndfg(jmax+1,jmax  ,ixtimep,.false.)
        call mklatbndfg(jmax+2,jmax-1,ixtimep,.false.)
      else
!
! First task makes lats -1 and 0 from 1 and 2:
        if (mytid==0) then
          if (.not.any(lats_have(:,mytid)==1).or.
     |        .not.any(lats_have(:,mytid)==2)) then
            write(6,"('>>> advnce before mklatbndfg: task 0',
     |        ' must have lats 1 and 2 to make lats -1 and 0')")
            call mp_close
            stop 'mklatbndfg'
          endif
          call mklatbndfg(0     ,1     ,ixtimep,.false.)
          call mklatbndfg(-1    ,2     ,ixtimep,.false.)
!
! Last task makes lats jmax+1 and jmax+2 from jmax and jmax-1:
        elseif (mytid==ntask-1) then ! last task needs lats jmax+1,jmax+2
          if (.not.any(lats_have(:,mytid)==jmax).or.
     |        .not.any(lats_have(:,mytid)==jmax-1)) then
            write(6,"('>>> advnce before mklatbndfg: last task ',
     |        ' must have lats jmax and jmax-1 to make lats ',
     |        'jmax+1 and jmax+2')")
            call mp_close
            stop 'mklatbndfg'
          endif
          call mklatbndfg(jmax+1,jmax  ,ixtimep,.false.)
          call mklatbndfg(jmax+2,jmax-1,ixtimep,.false.)
        endif
      endif
#else
      call mklatbndfg(0     ,1     ,ixtimep,.false.)
      call mklatbndfg(-1    ,2     ,ixtimep,.false.)
      call mklatbndfg(jmax+1,jmax  ,ixtimep,.false.)
      call mklatbndfg(jmax+2,jmax-1,ixtimep,.false.)
#endif
!
! hdif1 saves nrh and kmh at jj in fnrh and fkmh (fgcom.h), using fg 
!   fields UNM and VNM at j+1 and j+2, and TNM and NMS at j+1.
! hdif1: WRITE fnrh,fkmh <- READ fg(nunm,nvnm,ntnm,nms)
!
#ifdef MPI
      if (istep==1.or.ntask==1) then
        do jj=-2,jmax
          call hdif1(jj,ixtimep)
        enddo
      else
        do jj=lat0-3,lat1
          call hdif1(jj,ixtimep)
        enddo
      endif
#else
      jmaxv = jmax
!
!$OMP PARALLEL DO PRIVATE(jj) SHARED(jmaxv,ixtimep)
!MIC$ DOALL PRIVATE(jj) SHARED(jmaxv,ixtimep)
      do jj=-2,jmaxv
        call hdif1(jj,ixtimep)
      enddo
#endif
!
! hdif2 saves kld fields in fg at j+1 (0->37)
! (hdif2 uses fg, and nrh, kmh at j-1 from hdif1)
! hdif2: WRITE fg(nkldu,v,t,ps,ps2) <- READ fg(nunm,nvnm,ntnm,npsnm,nps2nm)
!
#ifdef MPI
      if (istep==1.or.ntask==1) then
        do jj=-1,jmax
          call hdif2(jj,ixtimep)
        enddo
      else
        do jj=lat0-2,lat1
          call hdif2(jj,ixtimep)
        enddo
      endif
#else
      jmaxv = jmax
!
!$OMP PARALLEL DO PRIVATE(jj) SHARED(jmaxv,ixtimep)
!MIC$ DOALL PRIVATE(jj) SHARED(jmaxv,ixtimep)
      do jj=-1,jmaxv
        call hdif2(jj,ixtimep)
      enddo
#endif
!
! Report cpu time from top of time loop:
      cputime1 = cputime()
!     write(6,"('Time to main lat loop from beginning',
!    |  ' of iter: ',f8.2,' (secs)')") cputime1-cputime0
!
! Reinit cputime for lat loop timing:
      cputime0 = cputime()
!
! Main latitude loop:
!
! If this is *not* an MPI run (e.g. unicos or irix), lat0==1 and 
!   lat1==jmax, and OpenMP (irix) or MIC (unicos) distributes 
!   latitudes over number of processors.
! If this *is* an MPI run (e.g., aix), then lat0,lat1 are the first 
!   and last latitudes that each mpi task is responsible for. 
!   Data (e.g. updated boundary latitudes lat0-2,lat0-1,lat1+1,lat1+2
!   for 4th order finite differencing is passed via calls to mp_ 
!   routines in mpi_module.F.
!
!$OMP PARALLEL DO PRIVATE(jx,jj,i)
!
!MIC$  DOALL
!MIC$+   SHARED(lat0,lat1,c,cs,ixtimep,ixtimec,ndisk,njtmp,njnp)
!MIC$+   PRIVATE(j,jx,jj,i,racs)
!
      do 3000 jx=lat0,lat1
        j = jx                  ! j is in common in phys.h
!       write(6,"('advnce: jx=',i3,' lat0,1=',2i3)") jx,lat0,lat1
!
!       num_threads = omp_get_num_threads()
!       write(6,"('advnce begin main lat loop: jx=',i3,' istep=',
!    |    i6,' omp num_threads=',i3)") jx,istep,num_threads
!
        racs=1./(re*cs(jx))
!       write(6,"('jx=',i3,' racs=',e12.4)") jx,racs
!
! Each processor gets its own private f-array from the global fg-array,
! according to j (fg's 3rd dimension (latitude) is -1:jmax+2). Each
! f-array carries the j slice to be updated, and the 2 j slices on 
! either side, for finite differencing: j-2,j-1,j,j+1,j+2. These are
! from the ixtimep index of fg, i.e., the previous time step.
!
! j   njm2 njm1 nj njp1 njp2
! 1   -1    0    1    2   3
! 2    0    1    2    3   4
! 3    1    2    3    4   5
!  ...
! 34   32  33   34   35   36
! 35   33  34   35   36   37
! 36   34  35   36   37   38
!
        i = 0
        do jj=jx-2,jx+2
          i = i+1
          call putf(njtmp(i+1),jj,ixtimep,1,ndisk+2)
        enddo
!
! Dynamics driver (stores updated j in f(njnp)):
! (bulk of model code is called from sub dynamics)
        call dynamics(j,ixtimep)
!
! Save updated latitude slice in njnp to fg (at current time index):
        call putfg(njnp,jx,ixtimec,1,ndisk+2)
!
! End main parallel latitude loop:
!       write(6,"('advnce end main lat loop: jx=',i2)") jx
 3000 continue
!
! Report cpu time through main lat loop:
      cputime1 = cputime()
!     write(6,"('Time in main latitude loop = ',f8.2,' (secs)')")
!    |  cputime1-cputime0
!
#ifdef MPI
!
! Master mpi task must collect updated latitudes from
! slave tasks so it can calculate dynamo (dynamo is
! still serial code).
!
      if (ntask > 1) then
        if (mytid==0) then ! master: receive from slaves
          do jx=lat1+1,jmax
            call mp_rcvdyn(jx)
          enddo
        else               ! slaves: send to master
          do jx=lat0,lat1
            call mp_snddyn(jx,0)
          enddo
        endif
      endif
!
! Currently (2/00), only master mpi task can calculate dynamo:
      if (mytid==0) call dyn
!
! Distribute dynamo output back to slave tasks:
      if (ntask > 1) call mp_updatephi
!
! Write history output if necessary:
! All tasks use wrhist to determine if a history needs to be written.
!   If so, mp_fgtomaster transfers updated latitudes from slaves to 
!   the master, which can then write the history. 
! If its time to write a history, all tasks must call mp_fgtomaster, 
!   but only the master actually writes the history.
! If writing a secondary history, the master must receive lat slices
!   from the slaves for diagnostic fields (mp_sechtomaster).
!
      time2write = wrhist(istep,modeltime,
     |  wrprim,save_prim,newseries_prim, iseries_p,
     |  wrsech,save_sech,newseries_sech, iseries_s)
      if (time2write.and.ntask > 1) then
        call mp_fgtomaster
        if (wrsech) call mp_sechtomaster
      endif
      if (mytid==0) call outhist(istep,modeltime)
!
! Each mpi slave must receive the 2 lats before its first 
! updated lat (lat0-1,lat0-2), and the 2 lats after its last 
! updated lat (lat1+1,lat2+2). Each task must also send its 
! first 2 (lat0,lat0+1) and last 2 (lat1,lat1-1) updated lats
! to the appropriate "adjacent" tasks.
! 
      if (ntask > 1) call mp_bndlats(ixtimec)
!
! Non-mpi run: call dynamo and output:
#else
!
! Dynamo:
      call dyn
!
! Write history output if necessary:
      call outhist(istep,modeltime)
#endif
!
! Swap time indices of fg for next iteration.
! (current time becomes previous time for next iteration).
!
      itemp = ixtimec
      ixtimec = ixtimep
      ixtimep = itemp
!
! Return for next time iteration:
! If mpi job, set barrier.
!
      if (istep < nstep) then
#ifdef MPI
        if (ntask > 1) call mpi_barrier(MPI_COMM_WORLD,ier)
#endif
        goto 100
      endif
!
      write(6,"('Exit advance at iter  ',i6,' istep = ',i6)") iter,istep
      return
      end
!-----------------------------------------------------------------------
      subroutine advnceday
      use init_module,only: iter,iyear,iday,sfeps,sundec,sin_sundec,
     |  cos_sundec
      use cons_module,only: dt
      implicit none
!
! Local:
      integer :: idayit,idapr,iyrpr,iyr4,iyr100,lpyr,ienda
      real :: pi,delta,thet0
!
      IDAYIT = ITER*IFIX(dt)/86400
      IF (IDAYIT*86400 .EQ. ITER*IFIX(dt)) THEN
C       ****
C       ****     ADVANCE DAY
C       ****
        IDAPR = IDAY
        IYRPR = IYEAR
        IDAY = IDAY + 1
C       ****
C       ****     LPYR = 1(0)  IF IS (NOT) A LEAP YEAR
C       ****
        IYR4 = IYEAR/4
        IYR100 = IYEAR/100
        LPYR = 0
        IF (IYR4*4 .EQ. IYEAR .AND. IYR100*100 .NE. IYEAR) LPYR=1
        IENDA = 365 + LPYR
        IF (IDAY .GT. IENDA) THEN
         IYEAR = IYEAR + 1
         IDAY = IDAY - IENDA
        ENDIF                          !   FOR PAST YEAR'S END
C       ****
C       ****     RECALCULATE SUN'S DECLINATION
C       ****
        PI = 3.14159265358979
        DELTA = ATAN(TAN(23.5*PI/180.)*SIN(2.*PI*FLOAT(IDAY-80)/365.))
        sin_sundec = SIN(DELTA) ! C(95)
        cos_sundec = COS(DELTA) ! C(96)
!
! Update SFEPS:
C       ****  SFEPS IS 6% VARIATION IN SOLAR OUTPUT OVER A YEAR
C       ****  CAUSED BY THE ORBITAL ECCENTRICITY
        PI = 3.14159265358979
        THET0 = 2.*PI*FLOAT(iday)/365.
        SFEPS = 1.000110+0.034221*COS(THET0)+0.001280*SIN(THET0)
     1        +0.000719*COS(2.*THET0)+0.000077*SIN(2.*THET0)
!
        WRITE(6,"('ADVANCING DAY (PREVIOUS,PRESENT)=',4I5,' sfeps=',
     |   e12.4)") IDAPR,IYRPR,IDAY,IYEAR,sfeps
      ENDIF
      end subroutine advnceday
!-----------------------------------------------------------------------
      subroutine setamie(sec)
      implicit none
#include "params.h"
#include "amie.h"
!
! Args:
      real,intent(in) :: sec
!
! Local:
      integer :: ja,ia,ka
!
      IF (NTST.LT.NTIMS .AND. SEC.GT.SECUTA(2)) THEN
        NTST = NTST + 1
        SECUTA(1) = SECUTA(2)
        HPNHA(1) = HPNHA(2)
        HPSHA(1) = HPSHA(2)
        CPNHA(1) = CPNHA(2)
        CPSHA(1) = CPSHA(2)
        CUSPSTA(1) = CUSPSTA(2)
        CUSPSLA(1) = CUSPSLA(2)
        CUSPNTA(1) = CUSPNTA(2)
        CUSPNLA(1) = CUSPNLA(2)
        DO 120 JA=1,zjmx
          DO 120 IA=1,zimxp1
            POTKVA1(IA,JA) = POTKVA2(IA,JA)
            EKEVA1(IA,JA) = EKEVA2(IA,JA)
            EFLXA1(IA,JA) = EFLXA2(IA,JA)
  120   CONTINUE
        DO 121 KA=1,3
          DO 121 JA=1,zjmx
            DO 121 IA=1,zimxp1
              VIXYZA1(IA,JA,KA) = VIXYZA2(IA,JA,KA)
  121   CONTINUE
        READ (iuamie) SECUTA(2),HPSHA(2),HPNHA(2),CPSHA(2),CPNHA(2),
     1    CUSPSTA(2),CUSPNTA(2),CUSPSLA(2),CUSPNLA(2),POTKVA2,VIXYZA2,
     2    EKEVA2,EFLXA2
      ENDIF                    !  FOR READING AMIE
      end subroutine setamie
!-----------------------------------------------------------------------
      subroutine tstwrhist
!
! Routine to test wrhist (wrhist.f) and output logic (output.f). 
! Model is not run when this routine is called -- it only simulates 
! the time iteration loop for testing purposes. This routine is 
! optionally called from main tgcm.F instead of calling advnce.
! Can also set fakeflds in flds_mod.f for fast i/o.
!
      use input_module,only: start,step
      use hist_module,only: nstep,modeltime
!
! Local:
      integer(kind=8) :: nsecs
      logical,external :: wrhist
!
! Init for istep==0 at beginning of run:
      istep = 0
      modeltime(1:3) = start(:,1)
      modeltime(4) = 0
      call modeltime_to_nsecs(modeltime,nsecs) ! sets nsecs, includes day
      call outhist(istep,modeltime)
! 
! Time loop:
  100 continue
      istep = istep+1
      nsecs = nsecs+step
      call nsecs_to_modeltime(nsecs,modeltime) ! increments modeltime(4)
      call outhist(istep,modeltime)
!
! Next time-step, or done:
      if (istep < nstep) goto 100
      stop 'tstwrhist'
      end subroutine tstwrhist