!-----------------------------------------------------------------------
      subroutine printdiags(iyear,iday,ndays,ut,f107,f107a,ap,rzur)
!
! Report end-of-run diagnostics to stdout, similiar to old code.
!
      use params,only: nlev
      use flds_modelz,only: zp,zpht,gz
      use flds_atmos,only: tn,xno2,xnn2,xno,xno21d,xno21s,xnno,xn4s,
     |  xn2d,xno2,xno3,xno1d,xnco,xnco2,xnoh,xnho2,xnch4,xnh2,xnh2o2,
     |  xnhtot,xnh,xnh2o,xnhox,xnhe,xnarg,xnno,xnno2,ch3,ch3o,ch2o,
     |  ch3o2,ch3ooh,hco,sht,rho
      use flds_msisatm,only: tnms,zpms,zpmht,xo2m,xn2m,xom,xn4sm,xnhm,
     |  xnhem,xnargm,shtms,rhoms
      use composition,only: phij,phip,phie,fluxhox,ps_o2,ps_ox,ps_noz,
     |  ps_n4s,ps_ch4,ps_h2,ps_co2,ps_co,ps_hox,ps_h2o,ps_ar,ps_he,
     |  sumnap
      use airglow_module,only: zhoh,fohb,nohvalt,mxohvlev,ilow,ihigh,
     |  te5577,dung,dvng,dtng,teo200,duno2,dvno2,dtno2, 
     |  te6300,dunr,dvnr,dtnr
      use flds_airglow,only: e5577,eo200,e6300,eco215,eno53
      use flds_sodium,only: xnas,xnasp,xnan2p,xnaoh,xnao2,xnao,
     |  xnaco3,xnahco3,xnat,pnas
      use flds_fields,only: edy
      use flds_dissoc,only: do3chap,po3chap,do3hug,po3hug,do3har,
     |  po3har,do3hz,pdo3hz,do2hz,pdo2hz,do3srb,do2srb,pdo2srb,
     |  dh2osrb,dco2srb,do3src,do2src,pdo2src,dh2osrc,
     |  dco2src,dnchm,dichm,do2lya,do2t,do2euv,do3lya,do3euv,do3t,
     |  dco2t,dco2euv,dco2lya,dh2ot,dh2oeuv,dh2olya,dno2t,pnoeuv,
     |  dnoeuv,dnosrb,dh2o2t,dcoeuv,dch4lya,dch4t,dch4euv,xjdis
      use flds_heat,only: qchap,qpchap,qhug,qphug,qhar,qphar,qherz,
     |  qpherz,qsrb,qpsrb,qpsrc,qsrc,qntot,qnn,qnc,qn1d,qnlya,qjoul,
     |  gwheat,qo3pr,xlen,xnoc,xo3pc,xircool,xco2c,coolt,hce,hcm,
     |  xmetac,xmetah,xno21sc,xno21dc,qnic,qnnc,qntot,qn,qproh,
     |  xircool,xco2c,xnoc,xo3pc,ck,cp,cmu
      use flds_col,only: yeesrcj,yeesrch,cno,cno2,cnn2,clno3,
     |  txno,txno2,txnn2,txno3
      use flds_ionatm,only: xne,xneee,xiop4s,xio2p,xinop,xin2p,xinp,
     |  xihp,te,ti,sped,shal,parcd,alamxx,alamxy
      use flds_ionzrt,only: qian,qpro,qixray,qnp,qiaur,qcr,qn2p,
     |  qop,qtin,qnop,qo2p,dinolya
      use restart,only: qcool
      implicit none
!
! Args:
      integer,intent(in) :: iyear,iday,ndays
      real,intent(in) :: ut,f107,f107a,ap,rzur
!
! Local:
      integer :: ix,k
      character(len=120) :: char120
      real,dimension(nlev) :: dtot,dtot_ms,edyk,cmuu,cmuh
      real,dimension(nlev) :: qntot_dpd,qnn_dpd,qnc_dpd,qn1d_dpd,
     |  qnlya_dpd,qsrc_dpd,qsrb_dpd,qherz_dpd,qian_dpd,qjoul_dpd,
     |  qproh_dpd,xlen_dpd,qn_dpd,qhug_dpd,qhar_dpd,qchap_dpd,
     |  gwheat_dpd,qo3pr_dpd
      real,dimension(nlev) :: xircool_dpd,xco2c_dpd,xnoc_dpd,
     |  xo3pc_dpd,hcm_dpd,hce_dpd,coolt_dpd
      real,dimension(nlev) :: xno21dc_dpd,xno21sc_dpd,xmetac_dpd,
     |  xmetah_dpd

      real :: dum,sum_nas,sum_nat,fluxnat
      real :: pn2,po,po2,cf,cpp
      real,parameter :: ur = 8.3144e+7
!
      write(6,"(/,72('-'),/,'GLBMEAN post-run diagnostics:')")
!
! Current state:
      write(char120,"('iyear=',i4,' iday=',i4,' ndays=',i4,' ut=',
     |  f5.2,' F107=',f8.2,' F107a=',f7.2,' AP=',f5.2,' RZUR=',
     |  f5.2)") iyear,iday,ndays,ut,f107,f107a,ap,rzur
!
! Geopotential and neutral temperature (glbmean and msis):
      call fprint8(char120,nlev,
     |  zp , zpms , zpht , zpmht , tn , tnms , dum,dum,
     | 'ZP','ZPMS','ZPHT','ZPMHT','TN','TNMS', ' ',' ')
!
! Major species (glbmean and msis):
      char120 = ' '
      call fprint8(char120,nlev,
     |  zp, zpht, xno, xom  , xnn2 , xn2m  , xno2 , xo2m  ,
     | 'ZP','ZPHT','O' ,'O_MS','N2'  ,'N2_MS','O2'  ,'O2_MS')
!
! Number densities:
      call fprint8(char120,nlev,
     |  zp,  zpht  , xno21d , xno21s, dum,dum,dum,dum,
     | 'ZP','ZPHT' ,'XNO21D','XNO21S',' ',' ',' ',' ')

      char120 = 'Number Densities:'
      call fprint8(char120,nlev,
     |  zp,  zpht  ,xnno,xn4s , xn4sm  ,xn2d ,xnno2 ,dum,
     | 'ZP','ZPHT' ,'NO','N4S','N4S_MS','N2D','NO2',' ')
      call fprint8(char120,nlev,
     |  zp,  zpht ,xno3,xno1d,xno,xnco,xnco2,dum,
     | 'ZP','ZPHT','O3','O1D','O','CO','CO2',' ')
!
! HOX (see sub comphox in comp.F):
      write(6,"(/,2x,'FHOX        ','PHIJ        ','PHIE        ',
     |          'PHIP        ')")
      write(6,"(4(1pe12.4))") fluxhox,phij,phie,phip
!
! More number densities:
      char120 = ' '
      call fprint8(char120,nlev,
     |  zp,  zpht ,xnh2o,xnhox,xnh,xnhm  ,dum,dum,
     | 'ZP','ZPHT','H2O','HOX','H','H_MS',' ',' ')
!
      call fprint8(char120,nlev,
     |  zp,  zpht ,xnoh,xnho2,xnch4,xnh2,xnh2o2,xnhtot,
     | 'ZP','ZPHT','OH','HO2','CH4','H2','H2O2','HTOT')
!
      call fprint8(char120,nlev,
     |  zp,  zpht ,xnhe, xnhem ,xnarg, xnargm ,dum,dum,
     | 'ZP','ZPHT','HE','HE_MS','ARG','ARG_MS',' ',' ')
      call fprint8(char120,nlev,
     |  zp,  zpht , ch3 , ch3o2 , ch3ooh , ch3o , ch2o , hco,
     | 'ZP','ZPHT','CH3','CH3O2','CH3OOH','CH3O','CH2O','HCO')
!
! Volume mixing ratios:
! Total number densities: dtot(nlev), dtot_ms(nlev)
      dtot    = xno+xno2+xnn2  ! total number density (whole-array op (nlev))
      dtot_ms = xom+xo2m+xn2m  ! msis total number density
      char120 = 'Volume Mixing ratios'
!
! Major species, and msis:
      call fprint8(char120,nlev,
     |  zp, zpht ,xno/dtot,xom/dtot_ms,xnn2/dtot,xn2m/dtot_ms,xno2/dtot,
     |  xo2m/dtot_ms,
     | 'ZP','ZPHT',  'O'   ,  'O_MS'   ,  'N2'  ,   'N2_MS'  ,   'O2'  ,
     |   'O2_MS')
      call fprint8(char120,nlev,
     |  zp,  zpht , dtot, dtot_ms,xno3/dtot, dum, dum, dum,
     | 'ZP','ZPHT','XMT','XMTS'  ,  'O3'   , ' ', ' ', ' ')
!
! Nitrogen volume mixing ratios:
      call fprint8(char120,nlev,
     |  zp, zpht,xnno/dtot,xn4s/dtot,xn4sm/dtot_ms,xn2d/dtot,xnno2/dtot,
     |  dum,
     |'ZP','ZPHT','NO'    ,'N4S'    ,'N4S_MS'    ,'N2D'    ,'NO2',' ')
!
! O and C species volume mixing ratios:
      call fprint8(char120,nlev,
     |  zp,zpht ,xno3/dtot,xno1d/dtot,xno/dtot,xnco/dtot,xnco2/dtot,dum,
     | 'ZP','ZPHT','O3'   ,'O1D'     ,'O'     ,'CO'     ,'CO2'     ,' ')
!
! Hydrogen species volume mixing ratios:
      call fprint8(char120,nlev,
     |  zp ,zpht ,xnh2o/dtot,xnhox/dtot,xnh/dtot,xnhm/dtot_ms,xnoh/dtot,
     |                                                       xnho2/dtot,
     | 'ZP','ZPHT','H2O'    ,'HOX'     ,'H'     ,'H_MS'     ,'OH','HO2')

      call fprint8(char120,nlev,
     |  zp,  zpht ,xnch4/dtot,xnh2/dtot,xnh2o2/dtot,xnhtot/dtot,dum,dum,
     | 'ZP','ZPHT','CH4'     ,'H2'     ,'H2O2'     ,'HTOT'     ,' ',' ')
!
! More volume mixing ratios:
      call fprint8(char120,nlev,
     |  zp,  zpht ,xnhe/dtot,xnhem/dtot_ms,xnarg/dtot,xnargm/dtot_ms,
     |  dum,dum,
     | 'ZP','ZPHT','HE'     ,'HE_MS'      ,'ARG'     ,'ARG_MS',' ',' ')
!
! Chlorine sp:
      call fprint8(char120,nlev,
     |  zp,  zpht ,ch3/dtot,ch3o2/dtot,ch3ooh/dtot,ch3o/dtot,ch2o/dtot,
     |  hco/dtot,
     | 'ZP','ZPHT','CH3'   ,'CH3O2'   ,'CH3OOH'   ,'CH3O'   ,'CH2O'   ,
     | 'HCO')
!
! Mass mixing ratios:
      char120 = 'Mass mixing ratios of major and minor species'
      call fprint8(char120,nlev,
     | ps_o2,ps_ox,ps_noz,ps_n4s ,ps_ch4,ps_h2,ps_co2,dum,
     | 'O2' ,'O'  ,'NO'  ,'N4S'  ,'CH4' ,'H2' ,'CO2',' ')

      call fprint8(char120,nlev,
     | ps_co,ps_hox,ps_h2o,ps_ar,ps_he,dum,dum,dum,
     | 'CO' ,'HOX' ,'H2O' ,'AR' ,'HE' ,' ',' ',' ')
!
! OH emission (s.a., sub airglow in airglow.F): 
! write(6,"(1x,i3,2x,12(1pe10.2))") i,zhoh(i),fohb(ix,i)

      char120 = 'OH 62-BAND EMISSION'
      ix = (ihigh-ilow-1)*mxohvlev+ihigh+1 ! index to 1st dim of fohb
      call fprint8(char120,nohvalt,
     |  zhoh, fohb(ix,:), dum,dum,dum,dum,dum,dum,
     | 'ZPHT','OHEMISS' , ' ',' ',' ',' ',' ',' ')
!
! Other emissions (see airglow.F and mkemiss.F):
      char120 = 'Airglow emissions'
      call fprint8(char120,nlev,
     |  zp , zpht , e5577 , eo200 , e6300 , eco215 , eno53, dum,
     | 'ZP','ZPHT','E5577','EO200','E6300','ECO215','ENO53',' ') 
!
! Diffs in te6300 are due to diffs in xne in sub mke6300 (mkemiss.F):
      write(6,"(/,2X, 'TE5577=',E12.3,2X,'TEO200=',E12.3,2X,
     |  'TE6300=',E12.3,/)") TE5577,TEO200,TE6300
      write(6,"(2X, 'DTNG=',E12.3,2X,'DUNG=',E12.3,2X,'DVNG=',
     |  E12.3,/)") DTNG,DUNG,DVNG
      write(6,"(2X, 'DTNR=',E12.3,2X,'DUNR=',E12.3,2X,'DVNR=',
     |  E12.3,/)") DTNR,DUNR,DVNR
      write(6,"(2X,'DTNO2=',E12.3,2X,'DUNO2=',E12.3,2X,'DVNO2=',
     |  E12.3,/)") DTNO2,DUNO2,DVNO2
!
! Sodium, et.al:
      char120 = 'Sodium densities '
      call fprint8(char120,nlev,
     |  zp , zpht ,xnas ,xnasp ,xnan2p ,xnaoh ,xnao2,dum,
     | 'ZP','ZPHT','NAS','NAS+','NA+N2','NAOH','NAO2',' ')

      char120 = 'Sodium densities '
      call fprint8(char120,nlev,
     |  zp , zpht ,xnao ,xnaco3 ,xnahco3 ,xnat ,pnas,dum,
     | 'ZP','ZPHT','NAO','NACO3','NAHCO3','NAT','PNAS',' ')
!
      sum_nat = 0.
      sum_nas = 0.
      do k=2,nlev
        sum_nas = sum_nas+(xnas(k)+xnas(k-1))*
     |    0.5*(zpht(k)-zpht(k-1))*1.e+5 
        sum_nat = sum_nat+(xnat(k)+xnat(k-1))*
     |    0.5*(zpht(k)-zpht(k-1))*1.e+5 
        edyk(k)=edy(k)*sht(k)*sht(k)
      enddo ! k=1,nlev
      fluxnat=-edyk(2)*((xnat(3)-xnat(1))/((zpht(3)-zpht(1))*1.e+5)
     |         +xnat(2)/sht(2)+xnat(2)*(tn(3)-tn(1))/((zpht(3)
     |         -zpht(1))*1.e+5)/tn(2))
!
! 7/25/05 btf: sumnap diffs w/ old code:
      write(6,"(/,'COL NA = ',1pe12.4,' COL NAT = ',1pe12.4,
     |  ' COL NAP = ',1pe12.4,' FLUXNAT = ',1pe12.4)")
     |  sum_nas,sum_nat,sumnap,fluxnat
!
! Volume mixing ratios of sodium sp:
      char120 = 'Volume mixing ratios'
      call fprint8(char120,nlev,
     |  zp,zpht,xnas/dtot,xnasp/dtot,xnan2p/dtot,xnaoh/dtot,xnao2/dtot,
     |  dum,
     | 'ZP','ZPHT','NAS' ,'NAS+'    ,'NA+Nw'    ,'NAOH'    ,'NAO2',' ')

      call fprint8(char120,nlev,
     |  zp,zpht,xnao/dtot,xnaco3/dtot,xnahco3/dtot,xnat/dtot,dum,dum,
     | 'ZP','ZPHT','NAO' ,'NACO3'    ,'NAHCO3'    ,'NAT',' ',' ')
!
! Dissociation rates:
      char120 = 'Dissociation and heating rates'
      call fprint8(char120,nlev,
     |  zp,  zpht , do3chap  ,po3chap , qchap , qpchap ,dum,dum,
     | 'ZP','ZPHT','DO3CHAP','PO3CHAP','QCHAP','QPCHAP',' ',' ')
!
      call fprint8(char120,nlev,
     |  zp,  zpht , do3hug , po3hug , qhug , qphug ,dum,dum,
     | 'ZP','ZPHT','DO3HUG','PO3HUG','QHUG','QPHUG',' ',' ')
!
      call fprint8(char120,nlev,
     |  zp,  zpht , do3har , po3har , qhar , qphar ,dum,dum,
     | 'ZP','ZPHT','DO3HAR','PO3HAR','QHAR','QPHAR',' ',' ')
!
      call fprint8(char120,nlev,
     |  zp,  zpht , do3hz , pdo3hz , do2hz , pdo2hz , qherz , qpherz,
     | 'ZP','ZPHT','DO3HZ','PDO3HZ','DO2HZ','PDO2HZ','QHERZ','QPHERZ')
!
      call fprint8(char120,nlev,
     |  zp,  zpht , do3srb , do2srb , pdo2srb , qsrb , qpsrb, dum,
     | 'ZP','ZPHT','DO3SRB','DO2SRB','PDO2SRB','QSRB','QPSRB',' ')
!
      call fprint8(char120,nlev,
     |  zp,  zpht , dh2osrb , dco2srb , do3src ,dum,dum,dum,
     | 'ZP','ZPHT','DH2OSRB','DCO2SRB','DO3SRC',' ',' ',' ')
!
      call fprint8(char120,nlev,
     |  zp, zpht , do2src , pdo2src , qsrc , qpsrc , dh2osrc , dco2src,
     |'ZP','ZPHT','DO2SRC','PDO2SRC','QSRC','QPSRC','DH2OSRC','DCO2SRC')
!
      call fprint8(char120,nlev,
     |  zp, zpht , do2src , yeesrcj , qsrc , yeesrch ,dum,dum,
     |'ZP','ZPHT','DO2SRC','YEESRCJ','QSRC','YEESRCH',' ',' ')
!
! Ions and temps:
      char120 = 'Ions and temperatures'
      call fprint8(char120,nlev,
     |  zp, zpht , xne , xneee ,xiop4s,xio2p,xinop,xin2p,
     |'ZP','ZPHT','XNE','XNEEE','O+'  ,'O2+','NO+','N2+')

      call fprint8(char120,nlev,
     |  zp, zpht ,xinp,xihp, te , ti , tn ,dum,
     |'ZP','ZPHT','N+','H+','TE','TI','TN',' ')
!
! Conductivities and lamdas:
! 7/26/05 btf: Most of these have diffs w/ old code:
      char120 = 'Conductivities and lamdas'
      call fprint8(char120,nlev,
     |  zp, zpht  , xne , sped , shal , parcd , alamxx , alamxy,
     | 'ZP','ZPHT','XNE','SPED','SHAL','PARCD','ALAMXX','ALAMXY')
!
! Heating rates (ergs/gm-s)
!
      char120 = 'Heating rates ergs/(gm-s)'
!
! Diffs w/ old code in qnc.
      call fprint8(char120,nlev,
     |  zp, zpht  , qntot , qnn  , qnc  , qn1d , qnlya ,dum,
     | 'ZP','ZPHT','QNTOT','QNNC','QNIC','QN1D','QNLYA',' ')
!
! Diffs w/ old code in qjoul.
! (old code had QNAUR instead of QIAN) 
      call fprint8(char120,nlev,
     |  zp, zpht  , qsrc , qsrb , qherz , qian  , qjoul ,dum,
     | 'ZP','ZPHT','QSRC','QSRB','QHERZ','QIAN','QNJOUL',' ')

      call fprint8(char120,nlev,
     |  zp, zpht  , qhar , qhug , qchap , qo3pr , gwheat ,dum,
     | 'ZP','ZPHT','QHAR','QHUG','QCHAP','QO3PR','GWHEAT',' ')

      call fprint8(char120,nlev,
     |  zp, zpht  , qpro , xlen , qntot ,dum,dum,dum,
     | 'ZP','ZPHT','QPRO','XLEN','QNTOT',' ',' ',' ')

      call fprint8(char120,nlev,
     |  zp, zpht  , qntot , xircool , xco2c , xnoc , xo3pc ,dum,
     | 'ZP','ZPHT','QNTOT','XIRCOOL','XCO2C','XNOC','XO3PC',' ')

      coolt = xco2c+xnoc+xo3pc+xmetac-hcm-hce ! whole-array ops
      call fprint8(char120,nlev,
     |  zp, zpht  , hcm , hce , qntot , coolt ,dum,dum,
     | 'ZP','ZPHT','HCM','HCE','QNTOT','COOLT',' ',' ')

      call fprint8(char120,nlev,
     |  zp, zpht  , xno21dc , xno21sc , xmetac , xmetah ,dum,dum,
     | 'ZP','ZPHT','XNO21DC','XNO21SC','XMETAC','XMETAH',' ',' ')
!
      char120 = ' '
      call tqc
!
! Heating rates k/day:
!
      do k=1,nlev 
        pn2=xnn2(k)/dtot(k)
        po=xno(k)/dtot(k)
        po2=1.-pn2-po
        cf=0.5*ur*(po2*7./32.+pn2*7./28.+po*5./16.)
        cpp=86400./cf
        qntot_dpd(k) = qntot(k)*cpp
        qnn_dpd  (k) = qnn  (k)*cpp
        qnc_dpd  (k) = qnc  (k)*cpp
        qn1d_dpd (k) = qn1d (k)*cpp
        qnlya_dpd(k) = qnlya(k)*cpp
        qsrc_dpd (k) = qsrc (k)*cpp
        qsrb_dpd (k) = qsrb (k)*cpp
        qherz_dpd(k) = qherz(k)*cpp
        qian_dpd (k) = qian (k)*cpp
        qjoul_dpd(k) = qjoul(k)*cpp

        qhar_dpd  (k) = qhar  (k)*cpp
        qhug_dpd  (k) = qhug  (k)*cpp
        qchap_dpd (k) = qchap (k)*cpp
        qo3pr_dpd (k) = qo3pr (k)*cpp
        gwheat_dpd(k) = gwheat(k)*cpp
        qproh_dpd (k) = qproh (k)*cpp
        xlen_dpd  (k) = xlen  (k)*cpp
        qn_dpd    (k) = qn    (k)*cpp
      enddo ! k=1,nlev 

      char120 = 'Heating rates K/day'
!
! Diffs w/ old code in qnnc (also qntot_dpd == qnn_dpd?)
      call fprint8(char120,nlev,
     |  zp,  zpht ,qntot_dpd ,qnn_dpd,qnc_dpd,qn1d_dpd, qnlya_dpd, dum,
     | 'ZP','ZPHT','QNT'     ,'QNNC' ,'QNIC' ,'QN1D'  ,'QNLYA',' ')

! qnjoul diffs:
      call fprint8(char120,nlev,
     |  zp,  zpht , qsrc_dpd,qsrb_dpd,qherz_dpd,qian_dpd,qjoul_dpd,dum,
     | 'ZP','ZPHT','QSRC'   ,'QSRB'  ,'QHERZ'  ,'QIAN'  ,'QNJOUL',' ')

      call fprint8(char120,nlev,
     |  zp,  zpht , qhar_dpd,qhug_dpd,qchap_dpd,qo3pr_dpd,dum,dum,
     | 'ZP','ZPHT','QHAR'   ,'QHUG'  ,'QCHAP'  ,'QO3PR'  ,' ',' ')

      call fprint8(char120,nlev,
     |  zp,  zpht , gwheat_dpd,qproh_dpd,xlen_dpd,qn_dpd,dum,dum,
     | 'ZP','ZPHT','GWHEAT'   ,'QPROH'  ,'XLEN'  ,'QNTOT',' ',' ')

      do k=1,nlev 
        pn2=xnn2(k)/dtot(k)
        po=xno(k)/dtot(k)
        po2=1.-pn2-po
        cf=0.5*ur*(po2*7./32.+pn2*7./28.+po*5./16.)
        cpp=86400./cf
        xircool_dpd(k) = xircool(k)*cpp
        xco2c_dpd  (k) = xco2c(k)  *cpp
        xnoc_dpd   (k) = xnoc(k)   *cpp
        xo3pc_dpd  (k) = xo3pc(k)  *cpp
        hcm_dpd  (k)   = hcm(k)    *cpp
        hce_dpd  (k)   = hce(k)    *cpp
        coolt_dpd(k)   = coolt(k)  *cpp

        xno21dc_dpd(k) = xno21dc(k)*cpp
        xno21sc_dpd(k) = xno21sc(k)*cpp
        xmetac_dpd (k) = xmetac (k)*cpp
        xmetah_dpd (k) = xmetah (k)*cpp
      enddo ! k=1,nlev 

      call fprint8(char120,nlev,
     |  zp,zpht,qntot_dpd,xircool_dpd,xco2c_dpd,xnoc_dpd,xo3pc_dpd,dum, 
     |'ZP','ZPHT','QNTOT','XIRCOOL'  ,'XCO2C'  ,'XNOC'  ,'XO3PC',' ')

      call fprint8(char120,nlev,
     |  zp,  zpht , hcm_dpd, hce_dpd, qntot_dpd,coolt_dpd,dum,dum,
     | 'ZP','ZPHT','HCM'   ,'HCE'   ,'QNTOT'   ,'COOLT'  ,' ',' ')

      call fprint8(char120,nlev,
     |  zp,  zpht , qcool , xco2c_dpd, dum,dum,dum,dum,
     | 'ZP','ZPHT','QCOOL','XCO2C'   ,' ',' ',' ',' ')

      call fprint8(char120,nlev,
     |  zp,zpht ,xno21dc_dpd,xno21sc_dpd,xmetac_dpd,xmetah_dpd,dum,dum,
     | 'ZP','ZPHT','XNO21DC','XNO21SC'  ,'XMETAC'  ,'XMETAH',' ',' ')

      char120 = ' '
! diffs in edyk (maybe lbc only):
      call fprint8(char120,nlev,
     |  zp , zpht , edy , edyk , gz , sht , shtms ,dum,
     | 'ZP','ZPHT','EDY','EDYK','GZ','SHT','SHTMS',' ')

      call fprint8(char120,nlev,
     |  zp , zpht , rho , rhoms , ck , cmu , cp ,dum,
     | 'ZP','ZPHT','RHO','RHOMS','CK','CMU','CP',' ')

      do k=1,nlev
        cmuu(k) = cmu(k)/(rho(k)*sht(k)**2)
        cmuh(k) = cmu(k)/rho(k)
      enddo ! k=1,nlev

      call fprint8(char120,nlev,
     |  zp , zpht , edy , edyk , cmuu , cmuh , cmu ,dum,
     | 'ZP','ZPHT','EDY','EDYK','CMUU','CMUH','CMU',' ')

      call fprint8(char120,nlev,
     |  zp , zpht , cno , cno2 , cnn2 , clno3 ,dum,dum,
     | 'ZP','ZPHT','CNO','CNO2','CNN2','CLNO3',' ',' ')

      call fprint8(char120,nlev,
     |  zp , zpht , txno , txno2 , txnn2 , txno3 ,dum,dum,
     | 'ZP','ZPHT','TXNO','TXNO2','TXNN2','TXNO3',' ',' ')

      call fprint8(char120,nlev,
     |  zp , zpht , qtin , qop , qo2p , qnop , qn2p ,dum,
     | 'ZP','ZPHT','QTIN','QOP','QO2P','QNOP','QN2P',' ')

      call fprint8(char120,nlev,
     |  zp , zpht , qnp , qixray , qcr , qiaur , qpro ,dum,
     | 'ZP','ZPHT','QNP','QIXRAY','QCR','QIAUR','QPRO',' ')

      call fprint8(char120,nlev,
     |  zp , zpht , do2t , do2euv , do2lya , do2src , do2srb ,dum,
     | 'ZP','ZPHT','DO2T','DO2EUV','DO2LYA','DO2SRC','DO2SRB',' ')

      call fprint8(char120,nlev,
     |  zp , zpht , do2hz , pdo2hz , dnchm , dichm ,dum,dum,
     | 'ZP','ZPHT','DO2HZ','PDO2HZ','DNCHM','DICHM',' ',' ')

      call fprint8(char120,nlev,
     |  zp , zpht , do2srb , pdo2srb , do3har , po3har ,dum,dum,
     | 'ZP','ZPHT','DO2SRB','PDO2SRB','DO3HAR','PO3HAR',' ',' ')

      call fprint8(char120,nlev,
     |  zp , zpht , do3hug , po3hug , do3chap , po3chap ,dum,dum,
     | 'ZP','ZPHT','DO3HUG','PO3HUG','DO3CHAP','PO3CHAP',' ',' ')

      call fprint8(char120,nlev,
     |  zp , zpht , qhar , qphar , qhug , qphug , qchap , qpchap,
     | 'ZP','ZPHT','QHAR','QPHAR','QHUG','QPHUG','QCHAP','QPCHAP')

      call fprint8(char120,nlev, ! 9308
     |  zp , zpht , qsrb , qpsrb , qsrc , qpsrc ,dum,dum,
     | 'ZP','ZPHT','QSRB','QPSRB','QSRC','QPSRC',' ',' ')

      char120 = 'Dissociation rates'
      call fprint8(char120,nlev, ! 306
     |  zp , zpht , do3t , do3euv , do3lya , do3har ,dum,dum,
     | 'ZP','ZPHT','DO3T','DO3EUV','DO3LYA','DO3HAR',' ',' ')

      call fprint8(char120,nlev, ! 306
     |  zp , zpht , do3hug , do3chap , qherz , qpherz ,dum,dum,
     | 'ZP','ZPHT','DO3HUG','DO3CHAP','QHERZ','QPHERZ',' ',' ')

      call fprint8(char120,nlev, ! 324
     |  zp , zpht , dh2ot , dh2oeuv , dh2olya , dh2osrc , dh2osrb ,dum,
     | 'ZP','ZPHT','DH2OT','DH2OEUV','DH2OLYA','DH2OSRC','DH2OSRB',' ')

      call fprint8(char120,nlev, ! 324
     |  zp , zpht , dco2t , dco2euv , dco2lya , dco2src , dco2srb ,dum,
     | 'ZP','ZPHT','DCO2T','DCO2EUV','DCO2LYA','DCO2SRC','DCO2SRB',' ')

      pnoeuv = pnoeuv/xnn2 ! whole array op
      call fprint8(char120,nlev, ! 326
     |  zp , zpht , dnoeuv , dinolya , dnosrb , dno2t , pnoeuv ,dum,
     | 'ZP','ZPHT','DNOEUV','DINOLYA','DNOSRB','DNO2T','PNOEUV',' ')

      call fprint8(char120,nlev, ! 326
     |  zp , zpht , dch4t , dch4euv , dch4lya , dh2o2t , dcoeuv ,dum,
     | 'ZP','ZPHT','DCH4T','DCH4EUV','DCH4LYA','DH2O2T','DCOEUV',' ')
!
! xjdis is in flds_dissoc:
! real :: xjdis(nlev,44) ! see sub uars_solaruv in solheat.F
! xjdis(1)  == dbo2
! xjdis(2)  == dbh2o
! xjdis(3)  == dbco2
! xjdis(4)  == dbch4
! xjdis(5)  == dbho2
! xjdis(17) == dbch3ooh
!
      call fprint8(char120,nlev, ! 328
     |  zp , zpht , xjdis(:,1),xjdis(:,2),xjdis(:,3),xjdis(:,4),
     |              xjdis(:,5),xjdis(:,17),
     | 'ZP','ZPHT','DBO2'     ,'DBH2O'   ,'DBCO2'   ,'DBCH4'   ,
     |             'DBHO2'    ,'DBCH3OOH')
!
! xjdis(25)  == dbo3
! xjdis(26)  == dbn2o
! xjdis(27)  == dbh2o2
! xjdis(28)  == dbno2
!
      call fprint8(char120,nlev, ! 330
     |  zp,zpht,xjdis(:,25),xjdis(:,26),xjdis(:,27),xjdis(:,28),dum,dum, 
     | 'ZP','ZPHT','DBO3'  ,'DBN2O'    ,'DBH2O2'   ,'DBNO2',' ',' ')
!
! xjdis(29)  == dbh2o
! xjdis(39)  == dbo3-1d
! xjdis(42)  == dbch2o-2h
! xjdis(43)  == dbch2o-h2
!
      call fprint8(char120,nlev, ! 330
     |  zp,zpht,xjdis(:,29),xjdis(:,39),xjdis(:,42),xjdis(:,43),dum,dum, 
     | 'ZP','ZPHT','DBCH2O','DBO3-1D'  ,'DBCH2O-2H','DBCH2O-H2',' ',' ')


      end subroutine printdiags
!-----------------------------------------------------------------------
      subroutine tqc
      use params,only: nlev
      use flds_atmos,only: rho
      use flds_heat,only: qntot,qproh,xlen,gwheat,qchap,qo3pr,xircool,
     |  hcm,hce,xo3pc,xco2c,xnoc,qnlya,qsrc,qn1d,qnn,qnc,qsrb,qhar,
     |  qhug,qjoul,qherz
      use flds_ionzrt,only: qian
      use flds_modelz,only: zpht
      implicit none
!
! Local:
      integer :: k,klb
      real :: tqn,tqnn,tqnc,tqn1d,tqlya,tqsrc,tqsrb,tqherz,tqaur,tqjoul,
     |  tqhar,tqhug,tqchap,tco3pr,tgwh,tqpro,tqxlen,tcool,tcco2,tcnoc,
     |  tqo3pr,thcm,thce
      real :: rz,rzm,rw,dum
      character(len=120) :: char120
!
      tqn=0.   ;tqnn=0. ;tqnc=0.  ;tqn1d=0. ;tqlya=0.;tqsrc=0.;tqsrb=0.
      tqherz=0.;tqaur=0.;tqjoul=0.;tqhar=0. ;tqhug=0.;tqchap=0.
      tco3pr=0.;tgwh=0. ;tqpro=0. ;tqxlen=0.;tcool=0.;tcco2=0.;tcnoc=0.
      tqo3pr=0.;thcm=0. ;thce=0.
!
      klb = 18
      do k=klb,nlev
        rz = rho(k)
        rzm = rho(k-1)
        rw=0.5*(zpht(k)-zpht(k-1))*1.e+5
        tqn   =tqn    +(qntot(k)  *rz+qntot(k-1)  *rzm)*rw
        tqnn  =tqnn   +(qnn(k)    *rz+qnn(k-1)    *rzm)*rw
        tqnc  =tqnc   +(qnc(k)    *rz+qnc(k-1)    *rzm)*rw
        tqn1d =tqn1d  +(qn1d(k)   *rz+qn1d(k-1)   *rzm)*rw
        tqlya =tqlya  +(qnlya(k)  *rz+qnlya(k-1)  *rzm)*rw
        tqsrc =tqsrc  +(qsrc(k)   *rz+qsrc(k-1)   *rzm)*rw
        tqsrb =tqsrb  +(qsrb(k)   *rz+qsrb(k-1)   *rzm)*rw
        tqherz=tqherz +(qherz(k)  *rz+qherz(k-1)  *rzm)*rw
        tqaur =tqaur  +(qian(k)   *rz+qian(k-1)   *rzm)*rw
        tqjoul=tqjoul +(qjoul(k)  *rz+qjoul(k-1)  *rzm)*rw
        tqhar =tqhar  +(qhar(k)   *rz+qhar(k-1)   *rzm)*rw
        tqhug =tqhug  +(qhug(k)   *rz+qhug(k-1)   *rzm)*rw
        tqchap=tqchap +(qchap(k)  *rz+qchap(k-1)  *rzm)*rw
        tqo3pr=tqo3pr +(qo3pr(k)  *rz+qo3pr(k-1)  *rzm)*rw
        tgwh  =tgwh   +(gwheat(k) *rz+gwheat(k-1) *rzm)*rw
        tqpro =tqpro  +(qproh(k)  *rz+qproh(k-1)  *rzm)*rw
        tqxlen=tqxlen +(xlen(k)   *rz+xlen(k-1)   *rzm)*rw
        tcool =tcool  +(xircool(k)*rz+xircool(k-1)*rzm)*rw
        tcco2 =tcco2  +(xco2c(k)  *rz+xco2c(k-1)  *rzm)*rw
        tcnoc =tcnoc  +(xnoc(k)   *rz+xnoc(k-1)   *rzm)*rw
        tco3pr=tco3pr +(xo3pc(k)  *rz+xo3pc(k-1)  *rzm)*rw
        thcm  =thcm   +(hcm(k)    *rz+hcm(k-1)    *rzm)*rw
        thce  =thce   +(hce(k)    *rz+hce(k-1)    *rzm)*rw

      enddo ! k=klb,nlev

      rz=rho(klb+1)
      rzm=rho(klb)
      rw=0.5*(zpht(klb+1)-zpht(klb))*1.e+5
      thcm=(hcm(klb+1)*rz+hcm(klb)*rzm)*rw
      thce=(hce(klb+1)*rz+hce(klb)*rzm)*rw

      char120 = 'TQxxx (See sub tqc in advance.F)'
!
! tqnn diffs from old code (is same as tqn?).
      call fprint8(char120,1,
     |  tqn , tqnn , tqnc , tqn1d , tqlya ,dum,dum,dum, 
     | 'TQN','TQNN','TQNC','TQN1D','TQLYA',' ',' ',' ')

      call fprint8(char120,1,
     | tqsrc   , tqsrb , tqherz,  tqaur , tqjoul, dum,dum,dum,
     | 'HTQSRC','TQSRB','TQHERZ','TQAUR','TQJOUL',' ',' ',' ')

      call fprint8(char120,1,
     |  tqhar , tqhug , tqchap , tqo3pr , tgwh , tqpro , tqxlen , tqn, 
     | 'TQHAR','TQHUG','TQCHAP','TQO3PR','TGWH','TQPRO','TQXLEN','TQN')

      call fprint8(char120,1,
     |  tqn , tcool , tcco2 , tcnoc , tco3pr , thcm , thce, dum,
     | 'TQN','TCOOL','TCCO2','TCNOC','TCO3PR','THCM','THCE',' ')

      end subroutine tqc
