!
      module advance_module
      implicit none
!
! Advance the model in time.
!
      contains
!-----------------------------------------------------------------------
      subroutine advance
!
! Advance the model nstep time steps.
!
      use fields_module
      use params_module,only: nlon,nlat,nlev,nlevp1
      use hist_module,only: nstep,modeltime,nsource
      use init_module,only: istep,uthr,iter,secs,iday,iyear

      use input_module,only: start,step,iaurora=>aurora,idynamo=>dynamo,
     |  calendar_advance,gpi_ncfile,potential_model,meped_file,
     |  step_compqrj,current_pg,current_kq
      use input_module,only:
     |  power  ,power_time  ,ntimes_power, rd_power,
     |  ctpoten,ctpoten_time,ntimes_ctpoten, rd_ctpoten,
     |  bximf  ,bximf_time  ,ntimes_bximf,
     |  byimf  ,byimf_time  ,ntimes_byimf,
     |  bzimf  ,bzimf_time  ,ntimes_bzimf,
     |  swden  ,swden_time  ,ntimes_swden,
     |  swvel  ,swvel_time  ,ntimes_swvel,
     |  al     ,al_time     ,ntimes_al,
     |  kp     ,kp_time     ,ntimes_kp,
     |  f107   ,f107_time   ,ntimes_f107,
     |  f107a  ,f107a_time  ,ntimes_f107a
      use cons_module,only: dt
      use magfield_module,only: sunloc_apex
      use timing_module,only: timer,timing
      use mpitime_module,only: mpi_timer
      use gpi_module,only: rdgpi,getgpi
      use solgar_module,only: solgar_bndry,solgar_import
      use aurora_module,only: aurora_cons,echar_meped, epower_meped
      use qrj_module,only: init_sflux
      use pdynamo_module,only: pefield,prepare_phig3d
      use hdif_module,only: hdif1,hdif2
      use mgw_module,only: backgrnd
      use current,only: noso_crrt,noso_crdens
!
! Routines and timing for parallel dynamo (pdynamo.F):
      use pdynamo_module,only: dynamo_inputs,pdynamo
      use pdynamo_module,only: nmlat0,phihm ! high-lat heelis phi

      use heelis_module,only: heelis,colath
      use weimer_module,only: weimer01
      use output,only: outhist
      use meped,only: rdmeped
      use addfld_module,only: addfld
      use lbc,only: tuvz_lbc
!      
! current due to plasma pressure & gravity -> add to RHS of dynamo    
      use magpres_g, only: magpres_grav
!
#ifdef MPI
      use mpi
      use mpi_module,only: mytid,mp_gather2root,ntask,lat0,lat1,lon0,
     |  lon1,mp_polelats,mp_bndlats,mp_bndlons,mp_periodic_f4d
#else
      use mpi_module,only: lat0,lat1,lon0,lon1,mytid,ntask
#endif
      implicit none
!
! Local:
      integer(kind=8) :: nsecs ! total model time in seconds (includes day)
      integer :: itmp,ier,i,j,iprint
      real :: fmin,fmax,dday
      logical :: time2write,wrprim,wrsech,newseries_sech,
     |  newseries_prim,iseries_sech,iseries_prim
      integer :: lev0=1,lev1=nlevp1
!
! External:
      logical,external :: time2print
      real,external :: hp_from_kp,ctpoten_from_kp
!
      real :: secs_per_step,
     |  time0,          time1,
     |  time0_step,     time1_step,
     |  time0_dynamics, time1_dynamics,
     |  time0_dynamo,   time1_dynamo,
     |  time0_phist,    time1_phist,
     |  time0_shist,    time1_shist,
     |  time0_prep,     time1_prep,
     |  time0_init,     time1_init
      character(len=3) :: timing_type
!
      logical :: debug=.false.
!
! External:
      logical,external :: wrhist
!
      write(6,"(/,'Enter advance: iter=',i10,' nstep=',i10)") iter,nstep
!
      timing_type = 'sys'
      if (timing%rtc) timing_type = 'rtc'
      write(6,"('Advance: timing_type=',a)") timing_type
      timing_type = 'sys'
      if (timing%rtc) timing_type = 'rtc'
      call timer(time0_init,time1_init,'INIT',0,0)
!
! 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
!
! Source was read into data(:,:,:,itp), and output always writes from
! data(:,:,:,itc), so swap indices to echo itp source data:
      itmp = itp
      itp = itc
      itc = itmp
!
! Start phist timing (source hist may or may not be written):
      call timer(time0_phist,time1_phist,'PHIST',0,0) ! start phist
!
! Echo source history to primary history output if an initial run,
      time2write = wrhist(istep,modeltime,
     |  wrprim, newseries_prim, iseries_prim,
     |  wrsech, newseries_sech, iseries_sech)

      if (time2write) then
#ifdef MPI
!
! If an initial run, echo source history to output before first time step.
!   Because only subdomains are read by rdhist, and only the root task
!   outputs histories, the subdomain data must be gathered to the root.
!
        call mp_gather2root(itc,'prim')
        if (mytid==0) call outhist(istep,modeltime)
#else
        do i=1,nf4d
          foutput(:,lon0:lon1,lat0:lat1,i) =        ! foutput is in fields.F
     |      f4d(i)%data(:,lon0:lon1,lat0:lat1,itc)
        enddo
        tlbc_glb(lon0:lon1,lat0:lat1) = tlbc(lon0:lon1,lat0:lat1)
        ulbc_glb(lon0:lon1,lat0:lat1) = ulbc(lon0:lon1,lat0:lat1)
        vlbc_glb(lon0:lon1,lat0:lat1) = vlbc(lon0:lon1,lat0:lat1)

        call outhist(istep,modeltime)
#endif
!
! nstep may be zero if user only wants to copy source history:
        if (nstep==0) then
          write(6,"(/,'Advance: model is not advanced in time because ',
     |      ' start==stop.')")
          return
        endif
      endif ! time2write
      call timer(time0_phist,time1,'PHIST',1,0) ! suspend phist timing
      time1_phist = time1
!
! Reswap indices, so model reads itp data and updates itc data:
      itmp = itp
      itp = itc
      itc = itmp
!
      call timer(time0_init,time1_init,'INIT',1,0) ! end init timing
      if (timing%level >= 2)
     |  write(6,"('Time in INIT =          ',
     |    f12.3,' Advance: step ',i5)") time1_init,istep
!
! Main time loop:
!
  100 continue
      call timer(time0_step,time1_step,'STEP',0,0) ! start step timing
      call timer(time0_prep,time1_prep,'PREP',0,0) ! start prep timing
      iter = iter+1
      istep = istep+1
      iprint = 0
      if (istep==1) iprint = 1
      nsecs = nsecs+step
!
! Note model time is incremented even when calendar_advance==0
      call nsecs_to_modeltime(nsecs,modeltime) ! increments modeltime(4)
      secs = amod(float(iter)*dt,86400.)
      uthr = secs/3600.
!
! Report to stdout (data is actually from previous step):
      if (time2print(nstep,istep)) then
        if (istep > 1) then
          write(6,"('Step ',i6,' of ',i6,' mtime=',3i3,
     |      ' secs/step (',a,') =',f6.2)") istep,nstep,modeltime(1:3),
     |      timing_type,secs_per_step
        else
          write(6,"('Step ',i6,' of ',i6,' mtime=',3i3)")
     |      istep,nstep,modeltime(1:3)
        endif
      endif ! time2print
!
! Write output history if necessary (master task writes data%(:,:,:,itc)).
! This is called early in the timestep time2write is available to other
! modules who may need it.
!
      time2write = wrhist(istep,modeltime,
     |  wrprim, newseries_prim, iseries_prim,
     |  wrsech, newseries_sech, iseries_sech)
!
! Update calendar day and sun's location and declination:
! (sub advance_day also updates sfeps)
! (note iday is incremented only if calendar_advance > 0)
!
      if (calendar_advance > 0) call advance_day
      call sunloc_apex(iyear,iday,secs) ! locate sun's longitude
!
! Update lower boundaries of t,u,v,z.
!
      call tuvz_lbc
!
! Get Solomon-Garcia lower boundaries and imports:
!
      call solgar_bndry(iyear,iday,int(secs))
      call solgar_import(iyear,iday,int(secs))
!
! Read MEPED data if this is first timestep (see meped.F):
      if (istep==1.and.len_trim(meped_file) > 0) 
     |  call rdmeped(meped_file)
!
! Interpolate input parameters to current model time, if time-dependent
! values were read from input. If namelist read parameter indices_interp
! is zero, interpolation is not done.
!
      if (ntimes_ctpoten > 0)
     |  call set_index(ctpoten_time,ntimes_ctpoten,nsecs,ctpoten,
     |    'ctpoten')
      if (ntimes_power > 0)
     |  call set_index(power_time,ntimes_power,nsecs,power,'power')
      if (ntimes_bximf > 0)
     |  call set_index(bximf_time,ntimes_bximf,nsecs,bximf,'bximf')
      if (ntimes_byimf > 0)
     |  call set_index(byimf_time,ntimes_byimf,nsecs,byimf,'byimf')
      if (ntimes_bzimf > 0)
     |  call set_index(bzimf_time,ntimes_bzimf,nsecs,bzimf,'bzimf')
      if (ntimes_swden > 0)
     |  call set_index(swden_time,ntimes_swden,nsecs,swden,'swden')
      if (ntimes_swvel > 0)
     |  call set_index(swvel_time,ntimes_swvel,nsecs,swvel,'swvel')
      if (ntimes_al > 0)
     |  call set_index(al_time,ntimes_al,nsecs,al,'al')
      if (ntimes_kp > 0)
     |  call set_index(kp_time,ntimes_kp,nsecs,kp,'kp')
      if (ntimes_f107 > 0)
     |  call set_index(f107_time,ntimes_f107,nsecs,f107,'f107')
      if (ntimes_f107a > 0)
     |  call set_index(f107a_time,ntimes_f107a,nsecs,f107a,'f107a')
!
! Get gpi data if necessary:
!
      if (len_trim(gpi_ncfile) > 0) then
        if (istep==1) call rdgpi ! read gpi data
        call getgpi(iyear,iday,int(secs),iprint)
        if (time2print(nstep,istep)) then
          if (kp /= spval) then ! User-provided Kp will override GPI (see below):
            write(6,"('GPI run: istep=',i6,' f107=',f8.3,' f107a=',
     |        f8.3,' (User-provided Kp will be used to calculate ',
     |        'power and ctpoten')") istep,f107,f107a
          else 
            write(6,"('GPI run: istep=',i6,' power=',f8.3,' ctpoten=',
     |        f8.3,' f107=',f8.3,' f107a=',f8.3)") istep,power,ctpoten,
     |        f107,f107a
          endif
        endif
      endif
!
! If Kp was provided by the user, but power and ctpoten were not, 
! then calculate power and/or ctpoten from the given Kp. If both
! Kp and GPI file were provided (see above), then the power and
! ctpoten calculated from Kp will override GPI, and only f10.7
! will be used from the GPI data file. This is accomplished by
! checking rd_power and rd_ctpoten from namelist, then overwriting 
! power,ctpoten from GPI.
!
      if (kp /= spval) then ! Kp was provided in namelist input
        if (rd_power==spval) power = hp_from_kp(kp)
        if (rd_ctpoten==spval) ctpoten = ctpoten_from_kp(kp)
        if (len_trim(gpi_ncfile) > 0.and.time2print(nstep,istep))
     |    write(6,"('NOTE: User-provided Kp overrides GPI: ',
     |      'power_from_kp=',f8.3,' ctpoten_from_kp=',f8.3)") 
     |      power,ctpoten
      endif
!
! Update auroral parameters:
!
      call aurora_cons(iprint,iyear,iday,secs)
      if (len_trim(meped_file) > 0.and.time2print(nstep,istep))
     |  write(6,"('  MEPED epower=',f10.2,' echar=',f10.2)")
     |    epower_meped,echar_meped 
!
! Init gravity wave parameterization (mgw.F, see also mgwinit):
!
      call backgrnd
!
! init_sflux calls ssflux and euvac (in qrj_module, qrj.F):
!
      if (istep==1.or.mod(istep,step_compqrj/step)==0) 
     |  call init_sflux(iprint)
!
! Calculate electric field on magnetic subdomains, and regrid
! to geographic subdomains in ex,ey,ez.
!
      call pefield   ! parallel efield
!
! Calculate vc, barm, xnmbar, and z (tn,o2,ox,vn are input).
!
      call addiag(
     |  tn     (levd0,lond0,latd0,itp),     ! in
     |  o2     (levd0,lond0,latd0,itp),     ! in
     |  ox     (levd0,lond0,latd0,itp),     ! in
     |  he       (levd0,lond0,latd0,itp),   ! in
     |  vn     (levd0,lond0,latd0,itp),     ! in
     |  vc     (levd0,lond0,latd0,itp),     ! out
     |  barm   (levd0,lond0,latd0,itp),     ! out
     |  xnmbar (levd0,lond0,latd0),         ! out
     |  xnmbari(levd0,lond0,latd0),         ! out
     |  xnmbarm(levd0,lond0,latd0),         ! out
     |  z      (levd0,lond0,latd0,itp),     ! out (note z output is at itp)
     |  zg     (levd0,lond0,latd0),         ! out (z with varying gravity)
     |  n2     (levd0,lond0,latd0),         ! out (1.-o2-o1-he)
     |  lon0,lon1,1,nlevp1,lat0,lat1)
!
!     do j=lat0,lat1
!       call addfld('un',' ',' ',un(:,lon0:lon1,j,itp),
!    |    'lev',1,nlevp1,'lon',lon0,lon1,j)
!       call addfld('vc',' ',' ',vc(1:nlev,lon0:lon1,j,itp),
!    |    'lev',1,nlev,'lon',lon0,lon1,j)
!       call addfld('advbarm',' ',' ',barm(:,lon0:lon1,j,itp),
!    |    'lev',1,nlevp1,'lon',lon0,lon1,j)
!     enddo
!
#ifdef MPI
!
! Define boundary latitudes -1,0 across the south pole, and
!   nlat+1,nlat+2 across the north pole:
!
      call mp_polelats(itp)
!
! Update subdomain boundary latitudes, longitudes:
      f4d(:)%mpi = .true.
      call mp_bndlats(f4d,nf4d,itp)
      call mp_bndlons(f4d,nf4d,itp)
#else
      call mk_polelat( 0,1,itp)
      call mk_polelat(-1,2,itp)
      call mk_polelat(nlat+1,nlat  ,itp)
      call mk_polelat(nlat+2,nlat-1,itp)
#endif
!
! For Z, itc==itp (itp was set by addiag):
      z(:,:,:,itc) = z(:,:,:,itp)
!
! Horizontal diffusion, first step:
! hdif1 saves nrh and kmh at lats -2->nlat in nrh and kmh using
!   fields un_nm and vn_nm at j+1 and j+2, and tn_nm and barm at j+1.
! hdif1 calls mp_bnd[lat,lons]_kmh
!
      call hdif1(
     |  tn_nm(levd0,lond0,latd0,itp),
     |  un_nm(levd0,lond0,latd0,itp),
     |  vn_nm(levd0,lond0,latd0,itp),
     |  barm (levd0,lond0,latd0,itp),
     |  1,nlevp1,lon0,lon1,lat0,lat1)
!
! Horizontal diffusion, second step:
! hdif2 saves 3d diagnostics f3d (fkldt,u,v,o2,o) at j+1 (0->37)
! hdif2: WRITE (fkldt,u,v,o2,o1) <- READ (tn_nm,un_nm,vn_nm,o2_nm,ox_nm)
! The kld terms will be used later in hdif3.
!
      call hdif2(
     |  tn_nm (levd0,lond0,latd0,itp),  ! 4d input
     |  un_nm (levd0,lond0,latd0,itp),  ! 4d input
     |  vn_nm (levd0,lond0,latd0,itp),  ! 4d input
     |  o2_nm (levd0,lond0,latd0,itp),  ! 4d input
     |  ox_nm (levd0,lond0,latd0,itp),  ! 4d input
     |  n4s_nm(levd0,lond0,latd0,itp),  ! 4d input
     |  noz_nm(levd0,lond0,latd0,itp),  ! 4d input
     |  co_nm (levd0,lond0,latd0,itp),  ! 4d input
     |  co2_nm(levd0,lond0,latd0,itp),  ! 4d input
     |  h2o_nm(levd0,lond0,latd0,itp),  ! 4d input
     |  h2_nm (levd0,lond0,latd0,itp),  ! 4d input
     |  hox_nm(levd0,lond0,latd0,itp),  ! 4d input
     |  ch4_nm(levd0,lond0,latd0,itp),  ! 4d input
     |  ar_nm (levd0,lond0,latd0,itp),  ! 4d input
     |  he_nm (levd0,lond0,latd0,itp),  ! 4d input
     |  nat_nm(levd0,lond0,latd0,itp),  ! 4d input
     |  kldt ,      ! 3d output
     |  kldu ,      ! 3d output
     |  kldv ,      ! 3d output
     |  kldo2,      ! 3d output
     |  kldox,      ! 3d output
     |  kldn4s,     ! 3d output
     |  kldnoz,     ! 3d output
     |  kldco,      ! 3d output
     |  kldco2,     ! 3d output
     |  kldh2o,     ! 3d output
     |  kldh2,      ! 3d output
     |  kldhox,     ! 3d output
     |  kldch4,     ! 3d output
     |  kldar,      ! 3d output
     |  kldhe,      ! 3d output
     |  kldnat,     ! 3d output
     |  1,nlevp1,lon0,lon1,lat0,lat1)
!
! Dynamo calls Heelis (heelis.F), Weimer (wei01gcm.F), or neither
!   for high latitude electric potential, according to user-provided
!   namelist parameter "model_potential".
! Get high latitude (Heelis or other) colatitudes, NH pfrac, and poten phihm.
!  If Weimer is used, then theta0,phid etc is changed before use in aurora
!   in dynamics.
!
      if (potential_model == 'WEIMER') then
        if (debug) write(6,"('advance calling weimer: istep=',i5)") 
     |    istep
        call weimer01(iyear,iday,uthr)
        if (len_trim(gpi_ncfile) > 0.and.time2print(nstep,istep)) then
          write(6,"('NOTE: ctpoten from GPI is overwritten by the',
     |      ' Weimer potential model: ctpoten=',f9.3)") ctpoten
        endif
      elseif (potential_model == 'HEELIS') then
        if (debug) write(6,"('advance calling heelis: istep=',i5)")
     |    istep
        call heelis
      else  !  potential_model='NONE'
        do j=1,nmlat0
          do i=1,nmlonp1
            phihm(i,j) = 0.
          enddo ! i=1,nmlonp1
        enddo ! j=1,nmlat0
        call colath
      endif
      call timer(time0_prep,time1,'PREP',1,0) ! suspend step-prep timing
      time1_prep = time1
!
! Model dynamics. Sub dynamics calls physics, chemistry, and dynamics
! routines for current time step (dynamics.F).
!
      if (debug) write(6,"('advance calling dynamics: istep=',i5)") 
     |  istep
      call timer(time0_dynamics,time1_dynamics,'DYNAMICS',0,0)

      call dynamics(nstep,istep)

      call timer(time0_dynamics,time1_dynamics,'DYNAMICS',1,0)
      if (debug) write(6,"('advance after dynamics: istep=',i5)") 
     |  istep
      if (timing%level >= 2.and.time2print(nstep,istep))
     |  write(6,"('Time in DYNAMICS =     ',f12.3,
     |    ' Advance: step ',i5)") time1_dynamics,istep
!
! Start dynamo timing. This will include magpres_grav, prep_dynamo,
! and dynamo proper.
!
      call timer(time0_dynamo,time1_dynamo,'DYNAMO',0,0)
!
! Call dynamo if namelist idynamo > 0:
!
      if (idynamo > 0) then
!
! calculate addition to electrodynamo RHS (current due to plasma 
! pressure and gravity)
!
        if(current_pg > 0) then   
          call timer(time0,time1,"magpres_grav",0,0)
          call magpres_grav (
     |    z     (levd0,lond0,latd0,itp),   ! geopotential input
     |    te    (levd0,lond0,latd0,itp),   ! electron temperature K
     |    ti    (levd0,lond0,latd0,itp),   ! ion temperature K
     |    ne    (levd0,lond0,latd0,itp),   ! electron density 1/cm^3
     |    op    (levd0,lond0,latd0,itp),   ! O+ 1/cm^3
     |    nplus (levd0,lond0,latd0),       ! N+  1/cm^3
     |    n2p   (levd0,lond0,latd0),       ! N2+ 1/cm^3
     |    nop   (levd0,lond0,latd0),       ! NO+  1/cm^3
     |    o2p   (levd0,lond0,latd0,itp),   ! O2+  1/cm^3
     |    1,nlevp1,lon0,lon1,lat0,lat1)
          if (debug) write(6,"('advance after magpres_grav: istep=',
     |      i4)") istep
          call timer(time0,time1,"magpres_grav",1,0)
        endif 
!
! Prepare neutral inputs for parallel dynamo:
! Calculate vertical velocity wn:
        call calc_wn(
     |    tn   (lev0:lev1,lon0:lon1,lat0:lat1,itp),
     |    w    (lev0:lev1,lon0:lon1,lat0:lat1,itp),
     |    barm (lev0:lev1,lon0:lon1,lat0:lat1,itp),
     |    wn   (lev0:lev1,lon0:lon1,lat0:lat1),     ! output
     |    1,nlevp1,lon0,lon1,lat0,lat1)

        call mpi_timer('dynamo_inputs',0,0)
        call dynamo_inputs(
     |    un     (lev0:lev1,lon0:lon1,lat0:lat1,itp),
     |    vn     (lev0:lev1,lon0:lon1,lat0:lat1,itp),
     |    wn     (lev0:lev1,lon0:lon1,lat0:lat1),
     |    z      (lev0:lev1,lon0:lon1,lat0:lat1,itp),
     |    ped    (lev0:lev1,lon0:lon1,lat0:lat1),     ! from lamdas.F
     |    hall   (lev0:lev1,lon0:lon1,lat0:lat1),     ! from lamdas.F
     |    1,nlevp1,lon0,lon1,lat0,lat1)
        call mpi_timer('dynamo_inputs',1,0)
        if (debug) write(6,"('advance after dynamo_inputs: istep=',
     |    i5)") istep
!
! Call parallel dynamo (pdynamo.F has its own timing):
        call pdynamo
!
! Calculate current diagnostics (current.F90):
! (note sub nosocoef is called by the root task from complete_integrals 
!  in pdynamo. Sub nosocoef is a stub in current.F90 that calls sub noso_coef)
!
        if (current_kq > 0) then
          if (mytid==0) call noso_crrt
          call noso_crdens
        endif
!
! Prepare potential on geographic grid if writing history:
        call prepare_phig3d(wrprim,wrsech)
        if (debug) write(6,"('advance after pdynamo: istep=',i5)")
     |    istep
!
! If idynamo <= 0, dynamo is not called:
      else
        if (time2print(nstep,istep)) 
     |    write(6,"('Note: idynamo=',i3,' -> dynamo was not called')")
     |      idynamo
      endif ! idynamo
!
! End timing for dynamo (old or new)
      call timer(time0_dynamo,time1_dynamo,'DYNAMO',1,0)
      if (timing%level >= 2.and.time2print(nstep,istep))
     |  write(6,"('Time in DYNAMO =       ',
     |    f12.3,' Advance: step ',i5)") time1_dynamo,istep
#ifdef VT
!     code = 127 ; state = 'timestep' ; activity='ModelCode'
      call vtend(127,ier)
#endif
!
! If its time to write a history, root task must gather subdomain
! data from slave tasks before writing the history. This is done
! by sub mp_gather2root. Only root task writes to the history.
!
      if (time2write) then
!
! History i/o timing is not included in STEP segment:
        call timer(time0_step,time1,'STEP',1,0) ! suspend step timing for i/o
        time1_step = time1

        if (wrprim) call timer(time0_phist,time1_phist,'PHIST',0,0) ! resume phist
        if (wrsech) call timer(time0_shist,time1_shist,'SHIST',0,0)
#ifdef MPI
        if (ntask > 1) then
          call mp_gather2root(itc,'prim')
          if (wrsech) call mp_gather2root(itc,'sech')
        elseif (ntask==1) then
          do i=1,nf4d
            foutput(:,lon0:lon1,lat0:lat1,i) =
     |        f4d(i)%data(:,lon0:lon1,lat0:lat1,itc)
          enddo
          tlbc_glb(lon0:lon1,lat0:lat1) = tlbc(lon0:lon1,lat0:lat1)
          ulbc_glb(lon0:lon1,lat0:lat1) = ulbc(lon0:lon1,lat0:lat1)
          vlbc_glb(lon0:lon1,lat0:lat1) = vlbc(lon0:lon1,lat0:lat1)
          tlbc_nm_glb(lon0:lon1,lat0:lat1)=tlbc_nm(lon0:lon1,lat0:lat1)
          ulbc_nm_glb(lon0:lon1,lat0:lat1)=ulbc_nm(lon0:lon1,lat0:lat1)
          vlbc_nm_glb(lon0:lon1,lat0:lat1)=vlbc_nm(lon0:lon1,lat0:lat1)
        endif
#else
!
! Update foutput if serial non-MPI run:
        do i=1,nf4d
          foutput(:,lon0:lon1,lat0:lat1,i) =
     |      f4d(i)%data(:,lon0:lon1,lat0:lat1,itc)
        enddo
        tlbc_glb(lon0:lon1,lat0:lat1) = tlbc(lon0:lon1,lat0:lat1)
        ulbc_glb(lon0:lon1,lat0:lat1) = ulbc(lon0:lon1,lat0:lat1)
        vlbc_glb(lon0:lon1,lat0:lat1) = vlbc(lon0:lon1,lat0:lat1)
        tlbc_nm_glb(lon0:lon1,lat0:lat1) = tlbc_nm(lon0:lon1,lat0:lat1)
        ulbc_nm_glb(lon0:lon1,lat0:lat1) = ulbc_nm(lon0:lon1,lat0:lat1)
        vlbc_nm_glb(lon0:lon1,lat0:lat1) = vlbc_nm(lon0:lon1,lat0:lat1)
#endif
!
! Write to history files:
        if (mytid==0) call outhist(istep,modeltime)
!
! Update i/o timing:
        if (wrprim) then
          call timer(time0_phist,time1,'PHIST',1,0) ! end phist timing
          time1_phist = time1_phist+time1
        endif
        if (wrsech) call timer(time0_shist,time1_shist,'SHIST',1,0)
        call timer(time0_step,time1,'STEP',0,0) ! resume step timing
      endif ! time2write
#ifdef MPI
!
! Each mpi task 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.
!
      call timer(time0_prep,time1,'PREP',0,0) ! resume prep timing
      f4d(:)%mpi = .true.
      call mp_bndlats(f4d,nf4d,itc)
!
! Same for longitude:
      call mp_bndlons(f4d,nf4d,itc)
!
! Periodic points for all updated fields:
      call mp_periodic_f4d(itc)
#else
! Non-mpi serial run:
      call mk_polelat(0     ,1     ,itc)
      call mk_polelat(-1    ,2     ,itc)
      call mk_polelat(lat1+1,lat1  ,itc)
      call mk_polelat(lat1+2,lat1-1,itc)
      call set_periodic_f4d(itc)
#endif
!
! Swap field data time indices, so current updated data becomes previous
! step data for next step:
      itmp = itp
      itp = itc
      itc = itmp
!
      call timer(time0_prep,time1,'PREP',1,0) ! end prep timing
      time1_prep = time1_prep+time1
      call timer(time0_step,time1,'STEP',1,0) ! end step timing
      time1_step = time1_step+time1
      if (timing%level >= 2.and.time2print(nstep,istep)) then
        write(6,"('Time in PREP =         ',
     |    f12.3,' Advance: step ',i5)") time1_prep,istep
        write(6,"('Time in STEP =         ',
     |    f12.3,' Advance: step ',i5)") time1_step,istep
      endif
      secs_per_step = time1_step
!
! Return for next time step:
      if (istep < nstep) then
!
! Go back up for next timestep:
        goto 100
      endif
      end subroutine advance
!-----------------------------------------------------------------------
      subroutine advance_day
!
! Advance calendar day if needed. Also update sfeps.
!
      use init_module,only: iter,iyear,iday,sfeps,sundec,sin_sundec,
     |  cos_sundec
      use cons_module,only: dt,pi
      implicit none
!
! Local:
      integer :: idayit,idayprev,iyearprev,iyr4,iyr100,lpyr,ienda
      real :: theta0
!
      idayit = iter*int(dt)/86400
      if (idayit*86400 == iter*int(dt)) then
        idayprev = iday
        iyearprev = iyear
        iday = iday + 1
!
! lpyr = 1(0) if is (not) a leap year
        iyr4 = iyear/4
        iyr100 = iyear/100
        lpyr = 0
        if (iyr4*4 == iyear .and. iyr100*100 /= iyear) lpyr=1
        ienda = 365 + lpyr
        if (iday > ienda) then
          iyear = iyear + 1
          iday = iday - ienda
        endif                          !   for past year's end
!
! Recalculate sun's declination
!
! 1/27/05 btf: Changed local delta to sundec (sundec is in init_module, 
!  where it is initialized at beginning of run). This should not change 
!  any results, because only sin_sundec and cos_sundec are actually used 
!  (in chapman), but am updating sundec here so it can be used in future
!  code development.
!
        sundec = atan(tan(23.5*pi/180.)*sin(2.*pi*float(iday-80)/365.))
        sin_sundec = sin(sundec) ! C(95)
        cos_sundec = cos(sundec) ! C(96)
!
! Update sfeps:
! sfeps is 6% variation in solar output over a year
!   caused by the orbital eccentricity.
!
        theta0 = 2.*pi*float(iday)/365.
        sfeps = 1.000110+0.034221*cos(theta0)+0.001280*sin(theta0)
     1        +0.000719*cos(2.*theta0)+0.000077*sin(2.*theta0)
!
        write(6,"('Advancing day (previous,present)=',4i5,' sfeps=',
     |   e12.4,' sundec=',e12.4)") idayprev,iyearprev,iday,iyear,sfeps,
     |   sundec
      endif
      end subroutine advance_day
!-----------------------------------------------------------------------
      subroutine set_index(rindex,ntimes,msecs,outindex,name)
!
! User has provided time-dependent ctpoten (ctpoten_time) and/or
!   power (power_time) via namelist input. This routine interpolates
!   these inputs to the current model time msecs, returning outindex.
! This routine is called separately (from advance) for ctpoten and power,
!   (i.e., rindex will be either ctpoten_time or power_time from input).
! Note times in seconds are 8-byte integers.
!
      use params_module,only: 
     |  mxind_time ! max number of time-dependent solar index points
      implicit none
!
! Args:
      real,intent(in) :: 
     |  rindex(4,mxind_time) ! user input times and values (day,hour,min,value)
      integer,intent(in) :: 
     |  ntimes  ! number of valid time/values in rindex(:,1:ntimes)
      integer(kind=8),intent(in) ::
     |  msecs   ! current model time in seconds
      real,intent(out) :: outindex ! output interpolated value
      character(len=*),intent(in) :: name
!
! Local:
      integer :: i
      integer(kind=8) :: nsec0,nsec1
!
! External:
      integer(kind=8),external :: mtime_to_nsec
      real,external :: finterp_bigints
!
! If model time is beyond last rindex time, use last rindex time:
      nsec1 = mtime_to_nsec(int(rindex(1:3,ntimes)))
      if (msecs > nsec1) then
        outindex = rindex(4,ntimes)
        goto 100
      endif 
!
! Bracket model time:
      do i=1,ntimes-1
        nsec0 = mtime_to_nsec(int(rindex(1:3,i)))
        nsec1 = mtime_to_nsec(int(rindex(1:3,i+1)))
!
! If model time is at a provided time, interpolation is not necessary:
        if (nsec0 == msecs) then
          outindex = rindex(4,i)
          goto 100
        endif
        if (nsec1 == msecs) then
          outindex = rindex(4,i+1)
          goto 100
        endif
!
! Interpolate to model time msecs:
        if (msecs >= nsec0 .and. msecs <= nsec1) then
          outindex = finterp_bigints(rindex(4,i),rindex(4,i+1),nsec0,
     |      nsec1,msecs)
          goto 100
        endif
      enddo ! i=1,ntimes-1
!
! Error if model time could not be bracketed. This should not happen,
! but you never know...
      write(6,"(/,'>>> set_index: could not bracket model time ',
     |  i10)") msecs
      write(6,"('ntimes=',i3)") ntimes
      do i=1,ntimes
        write(6,"('  i=',i3,' ntimes=',i3,' day,hr,min=',3f7.2,
     |    ' value=',e12.4)") i,ntimes,rindex(1:3,i),rindex(4,i)
      enddo
      call shutdown('set_index')
!
! Report to stdout:
 100  continue
!     write(6,"('set_index: ',a,' = ',e12.4)") name,outindex
      end subroutine set_index
!-----------------------------------------------------------------------
      subroutine calc_wn(tn,omega,barm,wn,lev0,lev1,lon0,lon1,lat0,lat1)
      use cons_module,only: gask,grav
!
! Calculate vertical velocity wn for dynamo inputs:
!
! Args:
      integer,intent(in) :: lev0,lev1,lon0,lon1,lat0,lat1
      real,dimension(lev0:lev1,lon0:lon1,lat0:lat1),intent(in) :: 
     |  tn,omega,barm
      real,dimension(lev0:lev1,lon0:lon1,lat0:lat1),intent(out) :: wn
!
! Local:
      integer :: j,i,k
      real,dimension(lev0:lev1,lon0:lon1,lat0:lat1) :: scheight
!
      do j=lat0,lat1
        do i=lon0,lon1
          do k=lev0,lev1-1
            scheight(k,i,j) = gask*tn(k,i,j)/(.5*(barm(k,i,j)+
     |        barm(k+1,i,j))*grav)
            wn(k,i,j)=0.5*(omega(k,i,j)+omega(k+1,i,j))*scheight(k,i,j)
          enddo
          scheight(lev1,i,j) = 0.
          wn(lev1,i,j) = 0.
        enddo
      enddo
      end subroutine calc_wn
!-----------------------------------------------------------------------
      end module advance_module