!
      subroutine advance
!
! Advance model in time (this is called from main glbmean.F):
!
      use params,only: nlev
      use cons,only: boltz
      use input,only: 
     |  ndays,      ! number of days for time loop
     |  nutpday,    ! number of ut's for each day (inner time loop)
     |  hist,       ! ut interval to write histories
     |  start_day,start_year,
     |  istep,      ! current time step (1,nstep)
     |  nstep,      ! number of timesteps in this run
     |  f107,       ! solar flux
     |  f107a,      ! 80-day average solar flux
     |  ap,         ! solar index
     |  rzur        ! sunspot number
      use flds_modelz,only: zp, gz, zqht, zpht
      use flds_heat,only: qntot,qn,qnn,qno3,qchap,qhar,qhug
      use flds_fields,only: edy
      use flds_atmos,only: tn,rho,xno2,xnn2,xno,xnh,xnhe,sht,
     |  xnh2o,xno21d,xno21s

      use flds_atmos,only: xnhtot ! for debug

      use flds_msisatm,only: tnms,zpms,zpmht,xn2m,xom,xo2m
      use solarheat,only: solheat
      use denmodule,only: denmod
      use tetimodule,only: teti
      use composition,only: chemistry
      use output,only: writehist
      use temp,only: chemheat,glbmtemp,rhozmod
      use airglow_module,only: airglow
      implicit none
!
! Local:
      integer :: i,it,iu,k,iyd,iday,iyear
      real :: ut,dut,dt,xmb,press,h2omix,tnfp,dtfp,dum,
     |  dday,yrday
      real,dimension(nlev) :: edyk
      character(len=80) :: char80
      logical :: timehist
      logical :: debug = .false.
!
      dut = (24./float(nutpday)) ! timestep in hours
      dt = 3600.*dut             ! timestep in seconds
      iyd = (start_year-start_year/100*100)*1000+start_day
      iday = start_day
      timehist = .true.
      ut = 0.
      istep = 0
      nstep = ndays*nutpday
      iyear = start_year
!     write(6,"('advance: iyear = start_year = ',i4,' start_day=',i4)") 
!    |  start_year,start_day
!
! Main time loop (days):
      do it=1,ndays
        if (debug) write(6,"('advance: it=',i4,' ndays=',i4)")
     |    it,ndays
!
! Each day is broken into nutpday equal fractions:
        do iu = 1,nutpday
          istep = istep+1
          dday = float(istep)*dut/24.
          ut = ut+dut
          if (ut==24.) then
            ut = 0.
            iday = iday+1
            if (iday==366) then
              iday = 1
              iyd = iyd+1
              iyear = iyear+1
            endif
          endif
          yrday = float((iyear-1)*365+iday)+ut/24.
!
! Report step to stdout:
          write(6,"('Step ',i5,' (of ',i5,'), iyear=',i4,' iday=',i4,
     |      ' dday=',f10.2,' ut=',f6.2)") istep,nstep,
     |      iyear,iday,dday,ut

!         write(6,"('advance: begin istep=',i5,' xnhtot=',/,(6e12.4))") 
!    |      istep,xnhtot
!
! Update height dependent quantities:
          do k=1,nlev
            zqht(k) = zpht(k)*1.e5
            gz(k)=980.665*(1.-3.14466e-4*zpht(k))
            rho(k)=(32.*xno2(k)+28.*xnn2(k)+16.*xno(k)+1.*xnh(k)+4.*
     |        xnhe(k))*1.66e-24
            xmb=(28.*xnn2(k)+32.*xno2(k)+16.*xno(k))/(xnn2(k)+
     |        xno2(k)+xno(k))
            sht(k)=boltz*tn(k)/(xmb*gz(k)*1.66e-24)
            edyk(k)=edy(k)*sht(k)*sht(k)
            press = 5.e-7*exp(-zp(k))
            h2omix = xnh2o(k)/(xno2(k)+xnn2(k)+xno(k))
            tnfp=6077.4/(28.548-alog(h2omix)-alog(press*1.e+2))
            dtfp=tn(k)-tnfp
          enddo ! k=1,nlev
!
! Diffs with old code in sht after 1st timestep:
!         write(char80,"('advance: iday=',i3,' ut=',f5.2,': update',
!    |      ' height-dependent quantities')") iday,ut
!         call fprint8(char80,nlev,
!    |      zqht , gz , sht , edy , dum,dum,dum,dum,
!    |     'zqht','gz','sht','edy',' ',' ',' ',' ')
!
! Update physical processes:
          if (debug) 
     |      write(6,"('advance: istep=',i4,': calling ratecoef')") istep
          call ratecoef

          if (debug) 
     |      write(6,"('advance: istep=',i4,': calling column')") istep
          call column

          if (debug) 
     |      write(6,"('advance: istep=',i4,': calling solheat')") istep
          call solheat(ut)

          if (debug) 
     |      write(6,"('advance: istep=',i4,': calling denmod')") istep
          call denmod(iday,ut,dt)

          if (debug) 
     |      write(6,"('advance: istep=',i4,': calling teti')") istep
          call teti(iday,ut,dt)

          if (debug) 
     |      write(6,"('advance: istep=',i4,': calling chemistry')")istep
          call chemistry(iday,ut,dt)

          if (debug) 
     |      write(6,"('advance: istep=',i4,': calling chemheat')") istep
          call chemheat

          if (debug) 
     |      write(6,"('advance: istep=',i4,': calling glbmtemp')") istep
          call glbmtemp(ut,dt)

          if (debug) 
     |      write(6,"('advance: istep=',i4,': calling rhozmod')") istep
          call rhozmod          
!
! Define some diagnostics (whole array ops):
          qntot = qn
          qnn   = qn
          qno3  = qchap+qhar+qhug
!
! Write history every HIST time steps:
          if (mod(istep,hist)==0) then
            call writehist(yrday,iyear,iday,dday,ut,timehist)
          endif

!         write(6,"('advance: end of istep=',i5,' xnhtot=',/,(6e12.4))") 
!    |      istep,xnhtot
!
! End time loops:
        enddo ! iu=1,nutpday
      enddo ! it=1,ndays
!
! Airglow package:
      call airglow(ut)
!
! MSIS:
      call nrlmsis(iyd,ut)
!
! Write to stdout (print statements similiar to old code):
! (printdiags.F)
      call printdiags(iyear,iday,ndays,ut,f107,f107a,ap,rzur)
!
! Output time-independent variables to history file:
      timehist = .false.
      call writehist(yrday,iyear,iday,dday,ut,timehist)
      end subroutine advance
