!
      subroutine dynamics(nstep,istep)
!
! Call dynamics layer of the model for current time step.
! Routines that do not have to do message-passing are generally called
!   from inside a latitude scan in this routine. Routines that do have
!   to make mpi calls are called once from this routine, and latitude
!   scan(s) are inside the called routine.
!
! Include the option of setting the lower boundary by using GSWM output
! this affects subroutine dt and duv
!
! am 6/06: change call to qjoule_ti to include ion and neutral density
!          factors, include vertical component
!          change call to qjoule_tn to include vertical velocity
!
      use fields_module
      use mpi_module,only: lon0,lon1,lat0,lat1,mytidi
#ifdef MPI
      use mpi_module,only: distribute_1d,mp_close
#endif
      use bndry_module,only: bndcmp
      use chemrates_module,only: chemrates_tdep
      use qrj_module,only: qrj,qtotal,qop,qo2p,qn2p,qnp,qnop
      use chapman_module,only: chapman
      use oplus_module,only: oplus
      use input_module,only: iaurora=>aurora
      use aurora_module,only: aurora
      use n4s_module,only: comp_n4s,minor_n4s
      use no_module,only: comp_no,minor_no
      use hdif_module,only: hdif3,hdif_bndlons,hdif_periodic
      use timing_module,only: timer,timing
      use dyndiag_module,only: dyndiag
!
      implicit none
!
! Args:
      integer,intent(in) :: nstep,istep
!
! VT vampir tracing:
!
#ifdef VT
#include <VT.inc>
#endif
!
! Local:
      integer :: i,j,n,lat,ier
      integer :: i0,i1,nk,nkm1,nlats
      logical,parameter :: debug=.false. ! add prints to stdout
      logical,parameter :: diags=.false. ! put fields on sec hist
      real :: time0,time1,time1_qrj,time1_cmpminor,time1_cmpmajor
!
! External:
      logical,external :: time2print
!
      if (debug) write(6,"(/,'Enter dynamics.')")
      i0=lon0 ; i1=lon1 ; nk=nlevp1 ; nkm1=nk-1 ! for addfsech
      nlats = lat1-lat0+1
!
      call bndcmp
      if (debug) write(6,"('dynamics after bndcmp')")
!
! Calculate specific heat and molecular viscosity.
! f4d itp fields are input, f3d fields are output.
!
      do lat=lat0,lat1
        call cpktkm(
     |    tn(levd0,lond0,lat,itp), ! 4d input
     |    o2(levd0,lond0,lat,itp), ! 4d input
     |    o1(levd0,lond0,lat,itp), ! 4d input
     |    cp(levd0,lond0,lat),     ! 3d output
     |    kt(levd0,lond0,lat),     ! 3d output
     |    km(levd0,lond0,lat),     ! 3d output
     |    1,nlevp1,lon0,lon1,lat)
      enddo ! lat=lat0,lat1
      if (debug) write(6,"('dynamics after cpktkm')")
      if (diags) then
        do lat=lat0,lat1
          call addfsech('CP' ,' ',' ',cp(:,i0:i1,lat),
     |      i0,i1,nk,nkm1,lat)
          call addfsech('KT' ,' ',' ',kt(:,i0:i1,lat),
     |      i0,i1,nk,nkm1,lat)
          call addfsech('KM' ,' ',' ',km(:,i0:i1,lat),
     |      i0,i1,nk,nkm1,lat)
        enddo
      endif
!
! Calculate omega for vertical velocity:
! (latitude scan is inside swdot)
!
      if (debug) write(6,"('dynamics before swdot')")
      call swdot(
     |  un(levd0,lond0,latd0,itp), ! un input
     |  vc(levd0,lond0,latd0,itp), ! vc input
     |  w (levd0,lond0,latd0,itc), ! omega output
     |  1,nlevp1,lon0,lon1,lat0,lat1,lat0,lat1)
      if (debug) write(6,"('dynamics after swdot')")
      if (diags) then
        do lat=lat0,lat1
          call addfsech('OMEGA' ,' ',' ',w(:,i0:i1,lat,itc),
     |      i0,i1,nk,nkm1,lat)
        enddo
      endif
!
! Calculate ExB ion drift velocities from electric field:
!
      call ionvel(
     |  z (levd0,lond0,latd0,itc), ! geopotential input
     |  ui(levd0,lond0,latd0),     ! U ion velocity output
     |  vi(levd0,lond0,latd0),     ! V ion velocity output
     |  wi(levd0,lond0,latd0),     ! W ion velocity output
     |  1,nlevp1,lon0,lon1,lat0,lat1)
      if (debug) write(6,"('dynamics after ionvel')")
      if (diags) then
        do lat=lat0,lat1
          call addfsech('UIVEL' ,' ',' ',ui(:,i0:i1,lat),
     |      i0,i1,nk,nkm1,lat)
          call addfsech('VIVEL' ,' ',' ',vi(:,i0:i1,lat),
     |      i0,i1,nk,nkm1,lat)
          call addfsech('WIVEL' ,' ',' ',wi(:,i0:i1,lat),
     |      i0,i1,nk,nkm1,lat)
        enddo
      endif
!
! Calculate column densities and line integrals of o2, o, and n2.
! Prognostics z,tn,o2,o1,barm are input, diagnostics vo2,vo1,vn2,
!   sco2,sco1,scn2 are output.
!
      do lat=lat0,lat1
        call chapman(
     |    z   (levd0,lond0,lat,itc), ! updated Z from addiag
     |    tn  (levd0,lond0,lat,itp), ! 4d in
     |    o2  (levd0,lond0,lat,itp), ! 4d in
     |    o1  (levd0,lond0,lat,itp), ! 4d in
     |    barm(levd0,lond0,lat,itp), ! 4d in
     |    vo2 (levd0,lond0,lat),     ! 3d out
     |    vo1 (levd0,lond0,lat),     ! 3d out
     |    vn2 (levd0,lond0,lat),     ! 3d out
     |    sco2(levd0,lond0,lat),     ! 3d out
     |    sco1(levd0,lond0,lat),     ! 3d out
     |    scn2(levd0,lond0,lat),     ! 3d out
     |    1,nlevp1,lon0,lon1,lat)
      enddo ! lat=lat0,lat1
      if (debug) write(6,"('dynamics after chapman')")
      if (diags) then
        do lat=lat0,lat1
          call addfsech('SCO2' ,' ',' ',sco2(:,i0:i1,lat),
     |      i0,i1,nk,nkm1,lat)
          call addfsech('SCO1' ,' ',' ',sco1(:,i0:i1,lat),
     |      i0,i1,nk,nkm1,lat)
          call addfsech('SCN2' ,' ',' ',scn2(:,i0:i1,lat),
     |      i0,i1,nk,nkm1,lat)
        enddo
      endif
!
! Calculate temperature dependent reaction rates (chemrates_module): 
!
      do lat=lat0,lat1
        call chemrates_tdep(
     |    tn  (levd0,lond0,lat,itp),
     |    te  (levd0,lond0,lat,itp),
     |    ti  (levd0,lond0,lat,itp),
     |    sco2(levd0,lond0,lat),
     |    vo2 (levd0,lond0,lat),
     |    1,nlevp1,lon0,lon1,lat)
      enddo
      if (debug) write(6,"('dynamics after chemrates_tdep')")
!
! Latitude scan for qrj, qinite, and aurora.
!
#ifdef VT
!     code = 122 ; state = 'ions' ; activity='ModelCode'
      call vtbegin(122,ier)
#endif
      do lat=lat0,lat1
!
! Calculate ionization and heating rates:
!
        if (debug) write(6,"('dynamics call qrj: lat=',i3)") lat
        call timer(time0,time1,'QRJ',0,0) ! start qrj timing for current lat
        call qrj(
     |    sco2(levd0,lond0,lat),
     |    sco1(levd0,lond0,lat),
     |    scn2(levd0,lond0,lat),
     |    tn  (levd0,lond0,lat,itp),
     |    no  (levd0,lond0,lat,itp),
     |    o2  (levd0,lond0,lat,itp),
     |    o1  (levd0,lond0,lat,itp),
     |    n4s (levd0,lond0,lat,itp),
     |    xnmbari(levd0,lond0,lat),
     |    1,nlevp1,lon0,lon1,lat)
        call timer(time0,time1,'QRJ',1,0) ! end qrj timing for current lat
        time1_qrj = time1_qrj+time1
        if (debug) write(6,"('dynamics after qrj: lat=',i3)") lat
          call addfsech('QRJ_QOP'    ,' ',' ',qop(:,i0:i1,lat),
     |      i0,i1,nk,nkm1,lat)
        if (diags) then
          call addfsech('QTOTAL' ,' ',' ',qtotal(:,i0:i1,lat),
     |      i0,i1,nk,nkm1,lat)
          call addfsech('QOP'    ,' ',' ',qop(:,i0:i1,lat),
     |      i0,i1,nk,nkm1,lat)
          call addfsech('QO2P'   ,' ',' ',qo2p(:,i0:i1,lat),
     |      i0,i1,nk,nkm1,lat)
          call addfsech('QN2P'   ,' ',' ',qn2p(:,i0:i1,lat),
     |      i0,i1,nk,nkm1,lat)
          call addfsech('QNP'    ,' ',' ',QNP(:,i0:i1,lat),
     |      i0,i1,nk,nkm1,lat)
          call addfsech('QNOP'   ,' ',' ',qnop(:,i0:i1,lat),
     |      i0,i1,nk,nkm1,lat)
        endif
!
! Calculate background (nightime) ionization:
        call qinite(
     |    o2     (levd0,lond0,lat,itp),
     |    o1     (levd0,lond0,lat,itp),
     |    no     (levd0,lond0,lat,itp),
     |    xnmbari(levd0,lond0,lat),
     |    vo2    (levd0,lond0,lat),
     |    vo1    (levd0,lond0,lat),
     |    vn2    (levd0,lond0,lat),
     |    1,nlevp1,lon0,lon1,lat)
        if (debug) write(6,"('dynamics after qinite: lat=',i3)") lat
!
! Add aurora particle precip and ionization:
        if (iaurora > 0) then
          call aurora(
     |      tn  (levd0,lond0,lat,itp),
     |      o2  (levd0,lond0,lat,itp),
     |      o1  (levd0,lond0,lat,itp),
     |      barm(levd0,lond0,lat,itp),
     |      1,nlevp1,lon0,lon1,lat)
          if (debug) write(6,"('dynamics after aurora: lat=',i3)") lat
        endif

        call addfsech('QAUR_QOP'    ,' ',' ',qop(:,i0:i1,lat),
     |    i0,i1,nk,nkm1,lat)

      enddo ! lat=lat0,lat1 for qrj, qinite, and aurora
      if (debug) write(6,"('dynamics after qrj,qinite,aurora')")
!
! Report timing from above latitude loop:
      if (timing%level >= 2.and.time2print(nstep,istep))
     |  write(6,"('Time in QRJ =          ',
     |    f12.2,' Dynamics: step ',i5)") time1_qrj,istep
!
#ifdef VT
!     code = 122 ; state = 'ions' ; activity='ModelCode'
      call vtend(122,ier)
#endif
!
! Advance O+:
! Oplus makes several calls to mp_bndlons_f3d, mp_periodic_f3d,
!   mp_gatherlons_f3d, and mp_scatterlons_f3.
! Full task subdomains (including ghost cells) are passed.
!
      call timer(time0,time1,'OPLUS',0,0) ! start oplus timing
      call oplus(
     |  tn    (levd0,lond0,latd0,itp),
     |  te    (levd0,lond0,latd0,itp),
     |  ti    (levd0,lond0,latd0,itp),
     |  o2    (levd0,lond0,latd0,itp),
     |  o1    (levd0,lond0,latd0,itp),
     |  n2d   (levd0,lond0,latd0,itp),
     |  ne    (levd0,lond0,latd0,itp),
     |  un    (levd0,lond0,latd0,itp),
     |  vn    (levd0,lond0,latd0,itp),
     |  w     (levd0,lond0,latd0,itp),
     |  barm  (levd0,lond0,latd0,itp),
     |  ui    ,
     |  vi    ,
     |  wi    ,
     |  xnmbar,
     |  op    (levd0,lond0,latd0,itp),
     |  op_nm (levd0,lond0,latd0,itp),
     |  op    (levd0,lond0,latd0,itc), ! out
     |  op_nm (levd0,lond0,latd0,itc), ! out
     |  xiop2p,                        ! out
     |  xiop2d,                        ! out
     |  1,nlevp1,lon0,lon1,lat0,lat1)
      if (debug) write(6,"('dynamics after oplus')")
      call timer(time0,time1,'OPLUS',1,0) ! end oplus timing
      if (timing%level >= 2.and.time2print(nstep,istep))
     |  write(6,"('Time in OPLUS =        ',
     |    f12.2,' Dynamics: step ',i5)") time1,istep
      if (diags) then
        do lat=lat0,lat1
          call addfsech('OPLUS' ,' ',' ',op(:,i0:i1,lat,itc),
     |      i0,i1,nk,nkm1,lat)
          call addfsech('XIOP2P',' ',' ',xiop2p(:,i0:i1,lat),
     |      i0,i1,nk,nkm1,lat)
          call addfsech('XIOP2D',' ',' ',xiop2d(:,i0:i1,lat),
     |      i0,i1,nk,nkm1,lat)
        enddo
      endif
!
! Electron density:
!
      do lat=lat0,lat1
        call elden(
     |    tn    (levd0,lond0,lat,itp),
     |    barm  (levd0,lond0,lat,itp),
     |    op    (levd0,lond0,lat,itp), ! O+ at current time step
     |    op    (levd0,lond0,lat,itc), ! O+ updated by sub oplus
     |    o2    (levd0,lond0,lat,itp),
     |    o1    (levd0,lond0,lat,itp),
     |    n2d   (levd0,lond0,lat,itp),
     |    no    (levd0,lond0,lat,itp),
     |    n4s   (levd0,lond0,lat,itp),
     |    xiop2p(levd0,lond0,lat),      ! from oplus
     |    xiop2d(levd0,lond0,lat),      ! from oplus
     |    z     (levd0,lond0,lat,itc),  ! updated Z from addiag
     |    nplus (levd0,lond0,lat),      ! N+  output
     |    n2p   (levd0,lond0,lat),      ! N2+ output
     |    nop   (levd0,lond0,lat),      ! NO+ output
     |    o2p   (levd0,lond0,lat,itc),  ! O2+ output
     |    ne    (levd0,lond0,lat,itc),  ! electron density output
     |    1,nlevp1,lon0,lon1,lat)
      enddo ! lat=lat0,lat1
      if (debug) write(6,"('dynamics after oplus')")
      do lat=lat0,lat1
        if (diags) then
          call addfsech('NPLUS'  ,' ',' ',nplus(:,i0:i1,lat),
     |      i0,i1,nk,nkm1,lat)
          call addfsech('ELD_N2P',' ',' ',n2p(:,i0:i1,lat),
     |      i0,i1,nk,nkm1,lat)
          call addfsech('ELD_NOP',' ',' ',nop(:,i0:i1,lat),
     |      i0,i1,nk,nkm1,lat)
          call addfsech('ELD_O2P',' ',' ',o2p(:,i0:i1,lat,itc),
     |      i0,i1,nk,nkm1,lat)
          call addfsech('ELD_NE' ,' ',' ',ne (:,i0:i1,lat,itc),
     |      i0,i1,nk,nkm1,lat)
        endif
      enddo ! lat=lat0,lat1

!
! Lamda ion drag coefficients:
!
      do lat=lat0,lat1
        call lamdas(
     |    tn  (levd0,lond0,lat,itp),
     |    barm(levd0,lond0,lat,itp),
     |    o2  (levd0,lond0,lat,itp),
     |    o1  (levd0,lond0,lat,itp),
     |    ti  (levd0,lond0,lat,itp),
     |    te  (levd0,lond0,lat,itp),
     |    o2p (levd0,lond0,lat,itp),
     |    op  (levd0,lond0,lat,itp),
     |    nop (levd0,lond0,lat),
     |    lxx (levd0,lond0,lat), ! LXX output
     |    lyy (levd0,lond0,lat), ! LYY output
     |    lxy (levd0,lond0,lat), ! LXY output
     |    lyx (levd0,lond0,lat), ! LYX output
     |    lam1(levd0,lond0,lat), ! LAM1 output
     |    ped (levd0,lond0,lat), ! pedersen conductivity output
     |    hall(levd0,lond0,lat), ! hall conductivity output
     |    1,nlevp1,lon0,lon1,lat)
        if (debug) write(6,"('dynamics after lamdas: lat=',i3)") lat
        if (diags) then
          call addfsech('LXX',' ',' ',lxx(:,i0:i1,lat),
     |      i0,i1,nk,nkm1,lat)
          call addfsech('LYY',' ',' ',lyy(:,i0:i1,lat),
     |      i0,i1,nk,nkm1,lat)
          call addfsech('LXY',' ',' ',lxy(:,i0:i1,lat),
     |      i0,i1,nk,nkm1,lat)
          call addfsech('LYX',' ',' ',lyx(:,i0:i1,lat),
     |      i0,i1,nk,nkm1,lat)
          call addfsech('LAM1','lam1','1/s',lam1(:,i0:i1,lat),
     |      i0,i1,nk,nkm1,lat)
          call addfsech('SIGMAPED',' ',' ',ped(:,i0:i1,lat),
     |      i0,i1,nk,nkm1,lat)
          call addfsech('SIGMAHAL',' ',' ',hall(:,i0:i1,lat),
     |      i0,i1,nk,nkm1,lat)
        endif
!
! Calculate n2d:
        call timer(time0,time1,'CMPMINOR',0,0) ! start comp minor timing
        call comp_n2d(
     |    o2  (levd0,lond0,lat,itp),
     |    o1  (levd0,lond0,lat,itp),
     |    no  (levd0,lond0,lat,itp),
     |    n4s (levd0,lond0,lat,itp),
     |    ne  (levd0,lond0,lat,itp),
     |    op  (levd0,lond0,lat,itp),
     |    n2p (levd0,lond0,lat),
     |    nop (levd0,lond0,lat),
     |    sco2(levd0,lond0,lat),
     |    xnmbarm(levd0,lond0,lat),
     |    n2d (levd0,lond0,lat,itc), ! N2D output
     |    1,nlevp1,lon0,lon1,lat)
        if (debug) write(6,"('dynamics after comp_n2d: lat=',i3)") lat
      enddo ! lat=lat0,lat1
      call timer(time0,time1,'CMPMINOR',1,0) ! suspend comp minor timing
      time1_cmpminor = time1_cmpminor+time1
      if (diags) then
        call addfsech('CMP_N2D',' ',' ',n2d(:,i0:i1,lat,itc),
     |    i0,i1,nk,nkm1,lat)
      endif
!
! Advance n4s:
!
      call timer(time0,time1,'CMPMINOR',0,0) ! resume comp minor timing
      do lat=lat0,lat1
        call comp_n4s(
     |    tn     (levd0,lond0,lat,itp),  ! TN (deg K)
     |    o2     (levd0,lond0,lat,itp),  ! O2 (mmr)
     |    o1     (levd0,lond0,lat,itp),  ! O (mmr)
     |    barm   (levd0,lond0,lat,itp),  ! mean mol weight
     |    xnmbarm(levd0,lond0,lat),      ! p0*e(-z)*barm
     |    no     (levd0,lond0,lat,itp),  ! NO
     |    n2d    (levd0,lond0,lat,itc),  ! N(2D) from comp_n2d
     |    ne     (levd0,lond0,lat,itp),  ! Ne
     |    o2p    (levd0,lond0,lat,itp),  ! O2+
     |    op     (levd0,lond0,lat,itp),  ! O+
     |    n2p    (levd0,lond0,lat),      ! N2+ from elden
     |    nplus  (levd0,lond0,lat),      ! N+  from elden
     |    nop    (levd0,lond0,lat),      ! NO+ from elden
     |    1,nlevp1,lon0,lon1,lat)
      enddo ! lat=lat0,lat1 
      if (debug) write(6,"('dynamics after comp_n4s')")
!
! Minor_n4s calls sub minor, which has 3d mpi calls, including gather/scatter
! for fft filtering. Full 3d task subdomains are passed.
!
      call minor_n4s(
     |    tn     (levd0,lond0,latd0,itp),  ! TN (deg K)
     |    o2     (levd0,lond0,latd0,itp),  ! O2 (mmr)
     |    o1     (levd0,lond0,latd0,itp),  ! O (mmr)
     |    n4s    (levd0,lond0,latd0,itp),  ! n4s from previous step
     |    n4s_nm (levd0,lond0,latd0,itp),  ! n4s at time n-1
     |    n4s    (levd0,lond0,latd0,itc),  ! output n4s
     |    n4s_nm (levd0,lond0,latd0,itc),  ! output n4s at time n-1
     |    1,nlevp1,lon0,lon1,lat0,lat1)
      if (debug) write(6,"('dynamics after minor_n4s')")
      call timer(time0,time1,'CMPMINOR',1,0) ! suspend comp minor timing
      time1_cmpminor = time1_cmpminor+time1
      if (diags) then
        call addfsech('CMP_N4S',' ',' ',n4s(:,i0:i1,lat,itc),
     |    i0,i1,nk,nkm1,lat)
      endif
! 
! Advance NO: 
!
      call timer(time0,time1,'CMPMINOR',0,0) ! resume comp minor timing
      do lat=lat0,lat1
        call comp_no(
     |    tn     (levd0,lond0,lat,itp),
     |    o2     (levd0,lond0,lat,itp),
     |    o1     (levd0,lond0,lat,itp),
     |    barm   (levd0,lond0,lat,itp),
     |    xnmbarm(levd0,lond0,lat),
     |    n4s    (levd0,lond0,lat,itp),
     |    n2d    (levd0,lond0,lat,itc),
     |    ne     (levd0,lond0,lat,itp),
     |    o2p    (levd0,lond0,lat,itp),
     |    1,nlevp1,lon0,lon1,lat)
      enddo ! lat=lat0,lat1 
      if (debug) write(6,"('dynamics after comp_no')")
!
! Minor_no calls minor, which has 3d mpi calls:
      call minor_no(
     |    tn   (levd0,lond0,latd0,itp),  ! TN (deg K)
     |    o2   (levd0,lond0,latd0,itp),  ! O2 (mmr)
     |    o1   (levd0,lond0,latd0,itp),  ! O (mmr)
     |    no   (levd0,lond0,latd0,itp),  ! current no
     |    no_nm(levd0,lond0,latd0,itp),  ! no at time n-1
     |    no   (levd0,lond0,latd0,itc),  ! output no
     |    no_nm(levd0,lond0,latd0,itc),  ! output no at time n-1
     |    1,nlevp1,lon0,lon1,lat0,lat1)
      if (debug) write(6,"('dynamics after minor_no')")
      call timer(time0,time1,'CMPMINOR',1,0) ! suspend comp minor timing
      time1_cmpminor = time1_cmpminor+time1
      if (timing%level >= 2.and.time2print(nstep,istep))
     |  write(6,"('Time in CMPMINOR =     ',f12.2,
     |    ' Dynamics: step ',i5)") time1_cmpminor,istep
      if (diags) then
        call addfsech('CMP_NO',' ',' ',no(:,i0:i1,lat,itc),
     |    i0,i1,nk,nkm1,lat)
      endif
!
! Calculate additions to neutral gas heating and O2 dissociation
! due to N, NO chemistry:
!
      do lat=lat0,lat1
        call qjnno(
     |    tn   (levd0,lond0,lat,itp),
     |    o2   (levd0,lond0,lat,itp),
     |    o1   (levd0,lond0,lat,itp),
     |    ne   (levd0,lond0,lat,itp),
     |    no   (levd0,lond0,lat,itp),
     |    n4s  (levd0,lond0,lat,itp),
     |    n2d  (levd0,lond0,lat,itc), ! from comp_n2d
     |    xnmbarm(levd0,lond0,lat),
     |    1,nlevp1,lon0,lon1,lat)
      enddo
      if (debug) write(6,"('dynamics after qjnno')")
!
! Calculate ion chemistry contribution to neutral gas heating
!   and O2 dissociation:
!
      do lat=lat0,lat1
        call qjion(
     |    tn    (levd0,lond0,lat,itp),
     |    o2    (levd0,lond0,lat,itp),
     |    o1    (levd0,lond0,lat,itp),
     |    o2p   (levd0,lond0,lat,itp),
     |    op    (levd0,lond0,lat,itp),
     |    n4s   (levd0,lond0,lat,itp),
     |    n2d   (levd0,lond0,lat,itp),
     |    no    (levd0,lond0,lat,itp),
     |    ne    (levd0,lond0,lat,itp),
     |    barm  (levd0,lond0,lat,itp),
     |    n2p   (levd0,lond0,lat),     ! from elden
     |    nplus (levd0,lond0,lat),     ! from elden
     |    nop   (levd0,lond0,lat),     ! from elden
     |    xiop2p(levd0,lond0,lat),     ! from oplus 
     |    xiop2d(levd0,lond0,lat),     ! from oplus
     |    1,nlevp1,lon0,lon1,lat)
      enddo ! lat=lat0,lat1
      if (debug) write(6,"('dynamics after qjion')")
!
! Electron and Ion temperatures:
!
      do lat=lat0,lat1
!
! Joule heating for ion temperature:
        call qjoule_ti(
     |    un (levd0,lond0,lat,itp), ! zonal neutral velocity input
     |    vn (levd0,lond0,lat,itp), ! meridional velocity input
     |    w  (levd0,lond0,lat,itp), ! vertical velocity input
     |    ui (levd0,lond0,lat),     ! zonal ion velocity input (ExB)
     |    vi (levd0,lond0,lat),     ! meridional ion velocity input (ExB)
     |    wi (levd0,lond0,lat),     ! vertical ion velocity input (ExB)
     |    lam1(levd0,lond0,lat),    ! lambda_1 ion drag coefficient [1/s]
     |    op (levd0,lond0,lat,itp), ! O+ 1/cm^3
     |    nplus(levd0,lond0,lat),   ! N+  1/cm^3
     |    n2p(levd0,lond0,lat),     ! N2+ 1/cm^3
     |    nop(levd0,lond0,lat),     ! NO+  1/cm^3
     |    o2p(levd0,lond0,lat,itp), ! O2+  1/cm^3
     |    barm(levd0,lond0,lat,itp),! mean molecular weight
     |    tn  (levd0,lond0,lat,itp),! neutral temperature
     |    qji_ti(levd0,lond0,lat),  ! ion joule heating output
     |    1,nlevp1,lon0,lon1,lat)
        if (debug) write(6,"('dynamics after qjoule_ti: lat=',i3)")lat
        if (diags) then
          call addfsech('QJI_TI',' ',' ',qji_ti(:,i0:i1,lat),
     |      i0,i1,nk,nkm1,lat)
        endif
!
        call settei(
     |    tn  (levd0,lond0,lat,itp),
     |    o2  (levd0,lond0,lat,itp),
     |    o1  (levd0,lond0,lat,itp),
     |    ne  (levd0,lond0,lat,itp),
     |    te  (levd0,lond0,lat,itp),
     |    ti  (levd0,lond0,lat,itp),
     |    op  (levd0,lond0,lat,itp),
     |    o2p (levd0,lond0,lat,itp),
     |    nop (levd0,lond0,lat),
     |    barm(levd0,lond0,lat,itp),
     |    qji_ti(levd0,lond0,lat),
     |    te  (levd0,lond0,lat,itc), ! output
     |    ti  (levd0,lond0,lat,itc), ! output
     |    1,nlevp1,lon0,lon1,lat)
        if (debug) write(6,"('dynamics after settei: lat=',i3)") lat
      enddo ! lat=lat0,lat1
      if (diags) then
        do lat=lat0,lat1
          call addfsech('TE_UPD',' ',' ',te(:,i0:i1,lat,itc),
     |      i0,i1,nk,nkm1,lat)
          call addfsech('TI_UPD',' ',' ',ti(:,i0:i1,lat,itc),
     |      i0,i1,nk,nkm1,lat)
        enddo
      endif
!
! Implicit and explicit cooling terms:
      do lat=lat0,lat1
        call newton(
     |    tn  (levd0,lond0,lat,itp),
     |    o2  (levd0,lond0,lat,itp),
     |    o1  (levd0,lond0,lat,itp),
     |    no  (levd0,lond0,lat,itp),
     |    barm(levd0,lond0,lat,itp),
     |    cp  (levd0,lond0,lat),
     |    cool_implicit(levd0,lond0,lat), ! output implicit cooling
     |    cool_explicit(levd0,lond0,lat), ! output explicit cooling
     |    1,nlevp1,lon0,lon1,lat)
        if (debug) write(6,"('dynamics after newton: lat=',i3)") lat
!
! Add O(3P) contribution to cooling terms:
        call newto3p(
     |    tn  (levd0,lond0,lat,itp),
     |    o1  (levd0,lond0,lat,itp),
     |    cp  (levd0,lond0,lat),
     |    cool_implicit(levd0,lond0,lat), ! output implicit cooling
     |    cool_explicit(levd0,lond0,lat), ! output explicit cooling
     |    1,nlevp1,lon0,lon1,lat)
        if (debug) write(6,"('dynamics after newto3p: lat=',i3)") lat
      enddo ! lat=lat0,lat1
      if (debug) write(6,"('dynamics after newton,newto3p')")
      if (diags) then
        do lat=lat0,lat1
          call addfsech('COOL_IMP',' ',' ',cool_implicit(:,i0:i1,lat),
     |      i0,i1,nk,nkm1,lat)
          call addfsech('COOL_EXP',' ',' ',cool_explicit(:,i0:i1,lat),
     |      i0,i1,nk,nkm1,lat)
        enddo
      endif
!
! Calculate horizontal diffusion terms for t,u,v,o2,o1, using
! coefficients from hdif2 (hdif2 was called from advance).
!
! Make boundary longitudes available for inputs to hdif3:
      call hdif_bndlons(kldt,kldu,kldv,kldo2,kldo1,
     |    1,nlevp1,lon0,lon1,lat0,lat1)
!
! Make hdt,u,v,o2,o1:
      do lat=lat0,lat1
        call hdif3(
     |    cp(levd0,lond0,lat),     ! specific heat input
     |    kldt ,                   ! 3d input from hdif2 for tn
     |    kldu ,                   ! 3d input from hdif2 for un
     |    kldv ,                   ! 3d input from hdif2 for vn
     |    kldo2,                   ! 3d input from hdif2 for o2
     |    kldo1,                   ! 3d input from hdif2 for o1
     |    hdt (levd0,lond0,lat),   ! 2d tn output at current lat
     |    hdu (levd0,lond0,lat),   ! 2d un output at current lat
     |    hdv (levd0,lond0,lat),   ! 2d vn output at current lat
     |    hdo2(levd0,lond0,lat),   ! 2d o2 output at current lat
     |    hdo1(levd0,lond0,lat),   ! 2d o1 output at current lat
     |    1,nlevp1,lon0,lon1,lat)
      enddo ! lat=lat0,lat1
      if (debug) write(6,"('dynamics after hdif3')")
!
! Periodic points for horizontal diffusion terms output from hdif3:
! This may not be necessary.
!     call hdif_periodic(hdt,hdu,hdv,hdo2,hdo1,
!    |    1,nlevp1,lon0,lon1,lat0,lat1)
!
! Calculate joule heating for tn (3d field qji_tn):
      do lat=lat0,lat1
        call qjoule_tn(
     |    un (levd0,lond0,lat,itp), ! zonal neutral velocity input
     |    vn (levd0,lond0,lat,itp), ! meridional velocity input
     |    w  (levd0,lond0,lat,itp), ! vertical velocity input
     |    ui (levd0,lond0,lat),     ! zonal ion velocity input
     |    vi (levd0,lond0,lat),     ! meridional ion velocity input
     |    wi (levd0,lond0,lat),     ! vertical ion velocity input
     |    lam1(levd0,lond0,lat),    ! lambda_1 ion drag coefficient
     |    barm(levd0,lond0,lat,itp),! mean molecular weight
     |    tn  (levd0,lond0,lat,itp),! neutral temperature
     |    qji_tn(levd0,lond0,lat),  ! ion joule heating output
     |    1,nlevp1,lon0,lon1,lat)
        if (diags) then
          call addfsech('QJI_TN',' ',' ',qji_tn(:,i0:i1,lat),
     |      i0,i1,nk,nkm1,lat)
        endif
!
! Calculate 2d (lon,lat) diagnostics for secondary histories:
! (see call dyndiag_sech below after lat loop)
!
        call dyndiag(
     |    ped (levd0,lond0,lat),     ! pedersen conductivity
     |    hall(levd0,lond0,lat),     ! hall conductivity
     |    qji_tn(levd0,lond0,lat),   ! ion Joule heating for tn
     |    z   (levd0,lond0,lat,itc), ! updated Z from addiag
     |    un  (levd0,lond0,lat,itp), ! neutral zonal velocity
     |    vn  (levd0,lond0,lat,itp), ! neutral meridional velocity
     |    ui  (levd0,lond0,lat),     ! ion zonal velocity
     |    vi  (levd0,lond0,lat),     ! ion meridional velocity
     |    1,nlevp1,lon0,lon1,lat)
      enddo ! lat=lat0,lat1
      if (debug) write(6,"('dynamics after qjoule_tn and dyndiag')")
!
! Advance neutral temperature:
!
      call timer(time0,time1,'DT',0,0) ! start dt timing
      call dt(
     |  tn	(levd0,lond0,latd0,itp),
     |  tn_nm (levd0,lond0,latd0,itp),
     |  un	(levd0,lond0,latd0,itp),
     |  vn	(levd0,lond0,latd0,itp),
     |  o2	(levd0,lond0,latd0,itp),
     |  o1	(levd0,lond0,latd0,itp),
     |  barm  (levd0,lond0,latd0,itp),    ! mean molecular weight
     |  cp	(levd0,lond0,latd0),	    ! specific heat	   (cpktkm.F)
     |  kt	(levd0,lond0,latd0),	    ! molecular diffusion  (cpktkm.F)
     |  km	(levd0,lond0,latd0),	    ! molecular viscosity  (cpktkm.F)
     |  hdt	(levd0,lond0,latd0),	    ! horizontal diffusion (hdif.F)
     |  qji_tn(levd0,lond0,latd0),	    ! joule heating (sub qjoule_tn)
     |  cool_implicit(levd0,lond0,latd0), ! implicit cooling (newton.F)
     |  cool_explicit(levd0,lond0,latd0), ! explicit cooling (newton.F)
     |  w	(levd0,lond0,latd0,itc),    ! updated W (swdot.F)
     |  tn	(levd0,lond0,latd0,itc),    ! output updated tn
     |  tn_nm (levd0,lond0,latd0,itc),    ! output updated tn at time n-1
     |  1,nlevp1,lon0,lon1,lat0,lat1)
      if (debug) write(6,"('dynamics after dt')")
      call timer(time0,time1,'DT',1,0) ! end dt timing
      if (timing%level >= 2.and.time2print(nstep,istep))
     |  write(6,"('Time in DT =           ',f12.2,
     |    ' Dynamics: step ',i5)") time1,istep
!
! Advance neutral velocities:
!
      call timer(time0,time1,'DUV',0,0) ! start duv timing
      call duv(
     |  tn	(levd0,lond0,latd0,itp), ! neutral temperature
     |  tn	(levd0,lond0,latd0,itc), ! updated neutral temperature (dt.F)
     |  tn_nm (levd0,lond0,latd0,itp), ! tn at time n-1
     |  un	(levd0,lond0,latd0,itp), ! zonal velocity
     |  vn	(levd0,lond0,latd0,itp), ! meridional velocity
     |  un_nm (levd0,lond0,latd0,itp), ! zonal velocity at time n-1
     |  vn_nm (levd0,lond0,latd0,itp), ! meridional velocity at time n-1
     |  w	(levd0,lond0,latd0,itc), ! updated vertical velocity (swdot.F)
     |  barm  (levd0,lond0,latd0,itp), ! mean molecular weight
     |  z	(levd0,lond0,latd0,itp), ! geopotential height
     |  hdu	(levd0,lond0,latd0),	 ! horizontal diffusion of U (hdif3.F)
     |  hdv	(levd0,lond0,latd0),	 ! horizontal diffusion of V (hdif3.F)
     |  ui	(levd0,lond0,latd0),	 ! zonal ion velocity input
     |  vi	(levd0,lond0,latd0),	 ! meridional ion velocity input
     |  lxx	(levd0,lond0,latd0),	 ! xx ion drag coefficient
     |  lyy	(levd0,lond0,latd0),	 ! yy ion drag coefficient
     |  lxy	(levd0,lond0,latd0),	 ! xy ion drag coefficient
     |  lyx	(levd0,lond0,latd0),	 ! yx ion drag coefficient
     |  km	(levd0,lond0,latd0),	 ! molecular viscosity  (cpktkm.F)
     |  un	(levd0,lond0,latd0,itc), ! output updated un
     |  un_nm (levd0,lond0,latd0,itc), ! output updated un at time n-1
     |  vn	(levd0,lond0,latd0,itc), ! output updated vn
     |  vn_nm (levd0,lond0,latd0,itc), ! output updated vn at time n-1
     |  1,nlevp1,lon0,lon1,lat0,lat1)
      if (debug) write(6,"('dynamics after duv')")
      call timer(time0,time1,'DUV',1,0) ! end duv timing
      if (timing%level >= 2.and.time2print(nstep,istep))
     |  write(6,"('Time in DUV =          ',f12.2,
     |    ' Dynamics: step ',i5)") time1,istep
!     
! Sources and sinks for major species composition O2, O:
!
      do lat=lat0,lat1
        call comp_o2o(
     |    tn   (levd0,lond0,lat,itp), ! neutral temperature
     |    o2   (levd0,lond0,lat,itp), ! molecular oxygen (mmr)
     |    o1   (levd0,lond0,lat,itp), ! atomic oxygen (mmr)
     |    barm (levd0,lond0,lat,itp), ! mean molecular weight
     |    op   (levd0,lond0,lat,itp), ! O+
     |    no   (levd0,lond0,lat,itp), ! NO
     |    n4s  (levd0,lond0,lat,itp), ! N(4S)
     |    n2d  (levd0,lond0,lat,itp), ! N(2D)
     |    o2p  (levd0,lond0,lat,itp), ! O2+
     |    ne   (levd0,lond0,lat,itp), ! electron density
     |    n2p  (levd0,lond0,lat),     ! N2+
     |    nplus(levd0,lond0,lat),     ! N+
     |    nop  (levd0,lond0,lat),     ! NO+
     |    1,nlevp1,lon0,lon1,lat)
      enddo ! lat=lat0,lat1
      if (debug) write(6,"('dynamics after comp_o2o')")
!
! Advance O2, O:
!
      time1_cmpmajor = 0.
      call timer(time0,time1,'CMPMAJOR',0,0) ! start comp major timing
      call comp(
     |  tn    (levd0,lond0,latd0,itp), ! neutral temperature
     |  o2    (levd0,lond0,latd0,itp), ! O2 (mmr)
     |  o2_nm (levd0,lond0,latd0,itp), ! O2 (mmr) at time n-1
     |  o1    (levd0,lond0,latd0,itp), ! O1 (mmr)
     |  o1_nm (levd0,lond0,latd0,itp), ! O1 (mmr) at time n-1
     |  un    (levd0,lond0,latd0,itp), ! zonal velocity
     |  vn    (levd0,lond0,latd0,itp), ! meridional velocity
     |  w     (levd0,lond0,latd0,itp), ! vertical velocity
     |  hdo2  (levd0,lond0,latd0),     ! O2 horizontal diffusion (hdif3)
     |  hdo1  (levd0,lond0,latd0),     ! O  horizontal diffusion (hdif3)
     |  o2    (levd0,lond0,latd0,itc), ! output: O2 updated for current step 
     |  o2_nm (levd0,lond0,latd0,itc), ! output: O2 updated for previous step
     |  o1    (levd0,lond0,latd0,itc), ! output: O1 updated for current step
     |  o1_nm (levd0,lond0,latd0,itc), ! output: O1 updated for previous step
     |  1,nlevp1,lon0,lon1,lat0,lat1)
      if (debug) write(6,"('dynamics after comp')")
      call timer(time0,time1,'CMPMAJOR',1,0) ! suspend comp major timing
      time1_cmpmajor = time1_cmpmajor+time1
      if (timing%level >= 2.and.time2print(nstep,istep))
     |  write(6,"('Time in COMP=          ',f12.2,
     |    ' Dynamics: step ',i5)") time1,istep
!
      if (timing%level >= 2.and.time2print(nstep,istep))
     |  write(6,"('Time in CMPMAJOR =     ',f12.2,
     |    ' Dynamics: step ',i5)") time1_cmpmajor,istep
      if (debug) write(6,"('dynamics returning')")
      end subroutine dynamics