      SUBROUTINE OUTPUT(N,NSTEPS,TIME,IZPRNTFLG12,EMAX,EMAXK)

CZ hydro64 added printout of phi at j+/- 1/2
CZ hydro62 removed most commented out code; replaced CR1 w/ RHOR2; removed CR2;
CZ hydro60-61 multiple changes to debug conservation of mass
CZ hydro59 replaced remaining references to TF's FLUX with FLUXKZ; added CON column;
CZ hydro58 replaced WTM with WTP; added variables to bdry/midpoint table;
CZ hydro57 removed unneeded info in diffusion table; changed table to include J=1 and J=NZ
CZ3 MEZ changed to use DHUA to calc flux (removes -u from calc)
CZ Added flags, changed outputs
      use parameters,only: NQ,NQT,NR,NSP,NSP1,NMAX,protonM,SMALL,MSKIP
      use flds_atmosphere
      use flds_chemistry
      use flds_transport
      use flds_species
      use flds_heating,only: q,q_cond,q_adv,q_grav
      implicit none

      integer :: I,J,K,N,NSTEPS,ISKIP,JSKIP,IROW,K1,K2,L,NQ1,NPE,LR
     | ,LDISPLAY
     | ,LDISPLAY1,LDISPLAY2

CZ Added flags; set = 1 to print 
      INTEGER :: IZPRNTFLG12
      INTEGER :: IZFLUXTRANS = 0 !Print flux intermediate products TRANS_XX for H at levels near exobase
      INTEGER :: DRCNTRL = 1

C      
      real :: TIME,TIMEY,DUSOL,RL,DIF,u_h,f_mass,FLUXH2,CHEMH2
     | ,ZFL,ZFT,H2ESC,H2BAL,sfac_tf
C
      double precision :: FUP(NQT),FLOW(NQT),CON(NQT),R12

     1 ,ZF(NZM),TAUAER(5),ZKM(NZM)
     3 ,CONKZ(NQT),CONIN(NQT)     
     
CZ      
      REAL :: EMAX,EMAXK
CZ hydro58 dimensioned these to save for printout      
      DOUBLE PRECISION,DIMENSION(NSP,NZM) :: FT1
     | ,FT2,FT3,FT3A,FT3B,FT3C,CP,CL,CONFRAC
CZ hydro64 added for flux at top boundary
     | ,FLUXKZP,FT1P,FT3P,FT3AP,FT3BP,FT3CP
     
      DOUBLE PRECISION,DIMENSION(NZM) :: RHOR2,RHOR2_M

      LDISPLAY = LO
CZ      
      LDISPLAY1 = LHE
      LDISPLAY2 = LH      
  
      PRINT*,"Exobase at level NZ =",NZ

      ISKIP = 1
      JSKIP = ISKIP
      IF(ISKIP.GT.1 .AND. N.EQ.NSTEPS) ISKIP = 2
C      
      TIMEY = TIME/3600./24./365.
      PRINT 100, TIME,TIMEY
 100  FORMAT(/1X,"TIME =",1PE9.2,5X,"TIMEY =",E9.2,1X,"YEARS")
C      PRINT 101,NPHOT
 101  FORMAT(/1X,"NPHOT =",I5)
C


      R12 = R(1)**2

      DO J=1,NZ
       RHOR2(J) = RHO(J)*(R(J)**2)
cz       CRM(J) = DR(J)*RHO_M(J)*(R_M(J)**2)
       RHOR2_M(J) = RHO_M(J)*(R_M(J))**2
       ZF(J) = zk(J) + 0.5*DRC(J)*1E-5       
      ENDDO

CZ *** START SPECIES LOOP K FOR FLUX CALCULATION AND PRINTOUT ***

 1917 FORMAT(//1x,"The following are the fluxes of species:",I4,2x,A8)
C
      DO 3 K=1,NQ
C
      IF (N .gt. 1) THEN
 
         IF (K .eq. LDISPLAY .OR. K .eq. LDISPLAY1 .OR. 
     |   K .EQ. LDISPLAY2) then
         PRINT 1917, K, ISPEC(K)
         PRINT 1918
         ENDIF
      ENDIF

CZ ** Start Altitude Loop for Flux calc and printout **

      DO 4 J=2,NZ-1 ! (can insert ,MSKIP)  hydro60: changed J+1 to J AND J to J-1     

cz       HD(K,J) = HDA(K,J)-U_ADV(J) !HD now is defined in main, now included in fmodules

cz      FT1(K,J) = RHO(J)/(WGT(K)*(DR(J)+DR(J-1)))  ! this factor gives dimensional flux

      FT1(K,J) = 1/(DR(J)+DR(J-1))  ! this factor gives actual flux (above gives dimensional flux)
      FT1P(K,J) = 1/(DR(J)+DR(J+1))  ! this factor gives actual flux (above gives dimensional flux)


cz      FT2(K,J) = (R(J)/R(J-1))**2 ! Removed to keep fluxes separate from r-squared

CZ If we use the equation for the lower boundary, j-1/2, then the math works and Cj = USOL(I,J)
      FT3A(K,J) = 2*(USOL(K,J) - USOL(K,J-1))*DK(K,J)
      FT3B(K,J) = DTK(K,J) + HD(K,J) 
      FT3C(K,J) = DR(J)*USOL(K,J-1)+DR(J-1)*USOL(K,J)

      FT3(K,J) = FT3A(K,J) + FT3B(K,J)*FT3C(K,J)

      FT3AP(K,J) = 2*(USOL(K,J+1)-USOL(K,J))*DK(K,J+1)
      FT3BP(K,J) = DTK(K,J+1) + HD(K,J+1)
      FT3CP(K,J) = DR(J+1)*USOL(K,J) + DR(J)*USOL(K,J+1)

      FT3P(K,J) = FT3AP(K,J) + FT3BP(K,J)*FT3CP(K,J)

CZ hydro63 changed sign per JFK: error in derivation
CZ hydro64 removed FT2, was ratio of r squared	    

      FLUXKZ(K,J)= -FT1(K,J)*FT3(K,J)

      FLUXKZP(K,J) = -FT1P(K,J)*FT3P(K,J)

      FT3A(K,J) = FT3A(K,J)*FT1(K,J)
      FT3B(K,J) = FT3B(K,J)*FT1(K,J)
      FT3C(K,J) = FT3C(K,J)*FT1(K,J)

      FT3AP(K,J) = FT3AP(K,J)*FT1P(K,J)
      FT3BP(K,J) = FT3BP(K,J)*FT1P(K,J)
      FT3CP(K,J) = FT3CP(K,J)*FT1P(K,J)
      
            DUSOL = USOL(K,J+1) - USOL(K,J)

       CONKZ(K) = 0

CZ These factors used for Dimensional Flux
cz        CP(K,J) = WGT(K)*YP(K,J)/RHO_M(J)
cz        CL(K,J) = YL(K,J)*USOL(K,J)

CZ These factors used for Actual Flux
        CP(K,J) = YP(K,J)
        CL(K,J) = YL(K,J)*UNUM(K,J)

cz      CONKZ(K) = CONKZ(K) + DR(J)*RHOR2_M(J)*(CP(K,J)-CL(K,J))
cz     |      +RHOR2(J)*FLUXKZ(K,J)-RHOR2(J+1)*FLUXKZ(K,J+1)

cz      CONLAYER(K,J-1) = DR(J-1)*RHOR2_M(J-1)*(CP(K,J-1)-CL(K,J-1))
cz     |      +RHOR2(J-1)*FLUXKZ(K,J-1)-RHOR2(J)*FLUXKZ(K,J)

cz      CONLAYER(K,J-1) = DR(J-1)*RHOR2_M(J-1)*(CP(K,J-1)-CL(K,J-1))
cz     |      +R2(J-1)*FLUXKZ(K,J-1)-R2(J)*FLUXKZ(K,J)
     
cz      CONLAYER(K,J) = DR(J)*RHOR2_M(J)*(CP(K,J)-CL(K,J))
cz     |      +FLUXKZ(K,J)-FLUXKZ(K,J+1)     


      IF (N .gt. 1) THEN
      IF (K .eq. LDISPLAY .OR. K .eq. LDISPLAY1 .OR. 
     |   K .EQ. LDISPLAY2) THEN
      PRINT 1919, J, R(J), DR(J)
     | , UNUM(K,J)
     | , RHOR2(J)*FLUXKZ(K,J)/(WGT(K)*(R(1)**2))
     | , FLUXKZ(K,J)/UNUM(K,J),U_ADV(J),
     | DU(K,J), DHUA(K,J), DTU(K,J),
     | USOL(K,J), DUSOL, DEN(J), YP(K,J), YL(K,J)*UNUM(K,J)
     | ,DR(J)*(YP(K,J)-YL(K,J)*UNUM(K,J))
     | ,DK(K,J+1),DTK(K,J+1),HD(K,J+1)
     | ,DMU(J),WT(J)
     | ,FT1(K,J),FT3(K,J),FT3A(K,J),FT3B(K,J),FT3C(K,J) 
      ENDIF
      endif
   
   4  CONTINUE  ! end J altitude loop for flux and conservation calcs

cz?? add negative sign
CZ Calc flux for top and bottom boundaries J = 1 and J = NZ
      FLUXKZ(K,NZ) = (FLUXKZ(K,NZ-1)+DR(NZ-1)*(YP(K,NZ-1)
     |      -YL(K,NZ-1)*UNUM(K,NZ-1))) !*(R(NZ)/R(NZ-1))**2

      FLUXKZ(K,1) = (FLUXKZ(K,2)-DR(1)*(YP(K,1)-YL(K,1)*UNUM(K,1)))

      FLUXKZP(K,1) = 0
      FLUXKZP(K,NZ) = 0


CZ These factors used for Dimensional Flux
cz        CP(K,J) = WGT(K)*YP(K,1)/RHO_M(1)
cz        CL(K,J) = YL(K,1)*USOL(K,1)

CZ These factors used for Actual Flux
        CP(K,1) = YP(K,1)
        CL(K,1) = YL(K,1)*UNUM(K,1)


CZ Now calculate Conservation for all layers - do this later because need flux of layer above and NZ
      do J=1,NZ-1
            CONLAYER(K,J) = DR(J)*(CP(K,J)-CL(K,J))
     |      + (RHOR2(J)*FLUXKZ(K,J)
     |      - RHOR2(J+1)*FLUXKZ(K,J+1))
     |      /(WGT(K)*(R_M(J)**2))
         if(K .NE. LHE) THEN
cz          PRINT*,"NOT HELIUM, CALC CONFRAC = CONLAYER/(P-L)"
         CONFRAC(K,J)=CONLAYER(K,J)/(DR(J)*(CP(K,J)-CL(K,J)))
        ELSE
cz	  PRINT*,"HELIUM, CALC CONFRAC = CONLAYER/RHO R2 FLUX"
         CONFRAC(K,J)=CONLAYER(K,J)*WGT(K)*(R_M(J)**2)
     |  	 /(RHOR2(J+1)*FLUXKZ(K,J+1))
        ENDIF
      
      ENDDO      
    

CZ *** Print out table of values for 3 selected species *********

      IF (N .gt. 1) THEN
      IF (K .eq. LDISPLAY .OR. K .eq. LDISPLAY1 .OR. 
     |   K .EQ. LDISPLAY2) THEN
       DO J=NZ-1,NZ
       
      PRINT 1919, J, R(J), DR(J)
     | , UNUM(K,J), RHOR2(J)*FLUXKZ(K,J)/(WGT(K)*(R(1)**2))
     | , FLUXKZ(K,J)/UNUM(K,J),U_ADV(J),
     | DU(K,J), DHUA(K,J), DTU(K,J),
     | USOL(K,J), DUSOL, DEN(J), YP(K,J), YL(K,J)*UNUM(K,J)
     | ,DR(J)*(YP(K,J)-YL(K,J)*UNUM(K,J))
     | ,DK(K,J+1),DTK(K,J+1),HD(K,J+1)
     | ,DMU(J),WT(J)
     | ,FT1(K,J),FT3(K,J),FT3A(K,J),FT3B(K,J),FT3C(K,J)
      
       ENDDO
      PRINT 1918
      endif
      ENDIF

 1918 FORMAT(3x,"J",7x,"R",8x,"DR",8x,"UNUM",6x,"FLUX"
     | ,7x,"ui",6x,"U_ADV",7x,"DU"
     | ,7x,"DHUA",6x,"DTU",7x,"USOL",6x,"DUSOL",5x,"DEN"
     | ,8x,"YP",7x,"YL*Ni",3x,"DR*(P-L)"
     | ,5X,"DK",7X,"DTK",7X,"HD",8X,"DMU",8X,"WT"
     | ,7X,"FT1",7X,"FT3"
     | ,7X ,"FT3A",6X,"FT3B",6X,"FT3C")

 1919 FORMAT(1I4,1x,1P1E11.4,1P25E10.2)

CZ *** Print out table of values for 3 selected sp on BOUNDARIES/MIDPOINTS *********
      IF (N .gt. 1) THEN
      IF (K .eq. LDISPLAY .OR. K .eq. LDISPLAY1 .OR. 
     |   K .EQ. LDISPLAY2) THEN
       print*," "
       PRINT 2022,K,ISPEC(K)
 2022  FORMAT("Number Density and Flux of species no.",I3
     |  ,", species ",A)
     
       print*," "
       print 2019
       print 2018
       DO J=1,NZ-1
       
CZ Boundaries       
       print 2020, J,(R(J)-R_PLT)/1.E5,T(J),U_ADV(J),RHO(J)
     | ,RHOR2(J)*FLUXKZ(K,J)/(WGT(K)*(R(1)**2))
     |  ,FT1(K,J),FT3(K,J)
     |  ,FT3A(K,J),FT3B(K,J),FT3C(K,J)
     | ,RHOR2(J)*FLUXKZ(K,J)/(WGT(K)*(R_M(J)**2))
     | ,FLUXKZ(K,J),FLUXKZP(K,J)
     | ,FT1P(K,J),FT3P(K,J),FT3AP(K,J),FT3BP(K,J)
     | ,FT3CP(K,J)

CZ Midpoints
       print 2021, J,(R_M(J)-R_PLT)/1.E5,RHO_M(J),UNUM(K,J),USOL(K,J)
     |  ,DR(J)*(YP(K,J)-YL(K,J)*UNUM(K,J))
     |  ,CONFRAC(K,J)
     | ,CONLAYER(K,J)
     | ,DR(J)*(CP(K,J)-CL(K,J))

     
cz     | ,RHOR2(J)*FLUXKZ(K,J)/(WGT(K)*(R_M(J)**2))
cz     | ,RHOR2(J+1)*FLUXKZ(K,J+1)/(WGT(K)*(R_M(J)**2))
cz     | ,TA1(K,J),TB1(K,J),TC1(K,J)     

       enddo ! end print out table through NZ-1

cz        CP(K,J) = WGT(K)*YP(K,J)/RHO_M(J)
cz        CL(K,J) = YL(K,J)*USOL(K,J)

CZ Upper Boundary       
       J=NZ
       print 2020, J,(R(J)-R_PLT)/1.E5,T(J),U_ADV(J),RHO(J),FLUXKZ(K,J)
     |  ,FT1(K,J),FT3(K,J)
     |  ,FT3A(K,J),FT3B(K,J),FT3C(K,J)       
     | ,RHOR2(J)*FLUXKZ(K,J)/(WGT(K)*(R_M(J)**2))
     | ,FLUXKZ(K,J),FLUXKZP(K,J)
     | ,FT1P(K,J),FT3P(K,J),FT3AP(K,J),FT3BP(K,J)
     | ,FT3CP(K,J)
       
       endif
       endif
       
      if(N .GT. 1) THEN
        print*," "
        print*,"CONLAYER(K,5)=",CONLAYER(K,5)
        print*,"K=",K,"  species =",ISPEC(K)
        PRINT*,"CP(K,5)-CL(K,5)=",CP(K,5)-CL(K,5)
        PRINT*,"RHOR2(5)=",RHOR2(5)      
        PRINT*,"RHOR2(6)=",RHOR2(6)
        print*,"FLUX(K,5)=",FLUXKZ(K,5)
        print*,"FLUX(K,6)=",FLUXKZ(K,6)
cz        print*,"R2(5)*FLUX(K,5)=",R2(5)*FLUXKZ(K,5)
cz        print*,"R2(6)*FLUX(K,6)=",R2(6)*FLUXKZ(K,6)	
        print*,"RHOR2(5)*FLUX(K,5)=",RHOR2(5)*FLUXKZ(K,5)
        print*,"RHOR2(6)*FLUX(K,6)=",RHOR2(6)*FLUXKZ(K,6)
	
      ENDIF

   3  CONTINUE  ! End Species Loop for Flux calcs and printout

CZ Heading
 2019 FORMAT(8X,"J",4X,"Alt",4x,"Temp",4X,"U_ADV",3X,"RHO/RHO_M"
     | ,3X,"UNUM",6X,"USOL",4X,"NORMFLUX",2x,"DR*(P-L)"
     | ,6X,"FT1",7X,"FT3"
     | ,7X ,"FT3A",6X,"FT3B",6X,"FT3C"
     | ,4X,"CONFRAC",4x,"CONLAYER",3X,"DR(P-L)",1X,"RHOR2FLUX/"
     | ,3X,"FLUX-",6X,"FLUX+",6x,"FT1P",6X,"FT3P"
     | ,5X,"FT3AP",5X,"FT3BP",5X,"FT3CP")

cz     | ,2X,"TA",6X,"TB",7X,"TC")
     
CZ Line     
 2018 format("-----------------------------------------------------------
     | -----------------------------------------------------------
     |------------------------------------")

CZ Boundaries
 2020 FORMAT("BDRY---",1I3,F7.1,1X,F6.0,1P2E10.2,20x,1P1E10.2
     | ,12x,1P5E10.2,31x,1P8E10.2)
     
CZ Midpoints     
 2021 FORMAT("MID    ",1I3,F7.1,17X,1P3E10.2,10x,1P1E10.2,52X,1P1E10.2
     | ,1X,1P2E10.2,10x,1P1E10.2)

CZ ********* END OF MASS CONSERVATION TABLE PRINTOUT


      print*," "
      PRINT*, "VEFF and UNUM:"
      PRINT 120, (VEFF(K),K=1,10)
      PRINT 120, (UNUM(K,NZ),K=1,10)
C
CZ Write H2 level at 100km, initial and final, to screen
      PRINT*, "H2 @ 100km = "
      PRINT*,USOL(LH2,1)
CZ Write H2, NSP, and NZ to idlpass.out 

      OPEN(UNIT=6903,FILE='idlpass.out')
    
      WRITE(6903,*)USOL(LH2,1)
      WRITE(6903,*)NSP
      WRITE(6903,*)NZ
      WRITE(6903,*)EMAX
      WRITE(6903,*)EMAXK
      
      CLOSE(6903)
      
CZ      WRITE(6903,*)EMAXK 


CZ Write same to idlpass.out      
CZ      WRITE(6903,*)USOL(LH2,1)
CZ      WRITE(6903,*)NSP      

CZ
 222  PRINT 105,N
 105  FORMAT(/1X,"MASS MIXING RATIOS OF LONG-LIVED SPECIES",
     | 1x,"N=",I5/)
      IROW = 10
      LR = NQ/IROW + 1
      RL = FLOAT(NQ)/IROW + 1
      DIF = RL - LR
      IF (DIF.LT.0.001) LR = LR - 1
C
      DO 8 L=1,LR
      K1 = 1 + (L-1)*IROW
      K2 = K1 + IROW - 1
      IF (L.EQ.LR) K2 = NQ
      
      PRINT 110, (ISPEC(K),K=K1,K2)
 110  FORMAT(/"Level",5x,"Z",7X,13(A8,2X)/)
cz      DO 20 I=1,1
cz  20  PRINT 120, zk(I),(USOL(K,I),K=K1,K2)
cz      DO 21 I=2,NZ,ISKIP
      DO 21 I=1,NZ,MSKIP      
  21  PRINT 129, I,zk(I),(USOL(K,I),K=K1,K2)
      PRINT 110, (ISPEC(K),K=K1,K2)
 120  FORMAT(1X,1P13E11.2) ! hydro58 was 10.2
 129  FORMAT(I3,1X,1P1E11.3,1P12E10.2)
C
      IF (N.EQ.0) GO TO 8
      PRINT 140
 140  FORMAT(/1X,"TP, TL")
      PRINT 145, (TP(K),K=K1,K2)
      PRINT 145, (TL(K),K=K1,K2)
 145  FORMAT(10X,1P12E10.2)

 
   8  CONTINUE
C
C
      IF(N.ge.1 .and. N.LT.NSTEPS) RETURN
C      
C      IF (N.EQ.0) RETURN
C      
C  Print full altitude profiles for plotting
C      WRITE(120,34) ISPEC(LH2),ISPEC(LCO),ISPEC(LH),
C     2  ISPEC(LO),ISPEC(LOH),ISPEC(LCO2),ISPEC(LN2)
  34  FORMAT("#",4X,"Z",7X,7(A8,1X))
C      DO 33 I=1,NZ
C  33  WRITE(120,32) zk(I),USOL(LH2,I),USOL(LCO,I),
C     2  USOL(LH,I),USOL(LO,I),UNUM(LCO2,I)/DEN(I),
C     3  UNUM(LN2,I)/DEN(I)
  32  FORMAT(1X,1P7E9.2)
C
C      PRINT 151,O3COL
 151  FORMAT(//1X,"OZONE COLUMN DEPTH = ",1PE11.4)
C
C ***** PRINT ON LAST ITERATION ONLY *****
C
C      IF (N.EQ.0) RETURN
C      
      DO 115 K=1,NQ


 115  CONTINUE
 
      FLUXH2 = FLUXKZ(LH2,NZ1)
      CHEMH2 = (YP(LH2,NZ) - YL(LH2,NZ)*UNUM(LH2,NZ))*DR(NZ)

C                                                         
      PRINT 125, N
 125  FORMAT(/1X,"NUMBER DENSITIES OF LONG-LIVED SPECIES",
     | 1x,"N=",I5/)               
      IROW = 10
      LR = NQ/IROW + 1
      RL = FLOAT(NQ)/IROW + 1
      DIF = RL - LR
      IF (DIF.LT.0.001) LR = LR - 1
C      
 1200  FORMAT(I3,1X,1P13E10.2)
 1100  FORMAT(/"Level",5x,"Z",8X,13(A8,2X)/)
      DO 9 L=1,LR                                         
      K1 = 1 + (L-1)*IROW                                 
      K2 = K1 + IROW - 1                                  
      IF (L.EQ.LR) K2 = NQ
      PRINT 1100, (ISPEC(K),K=K1,K2)                       
cz      DO 22 I=1,NZ,ISKIP
      DO 22 I=1,NZ,MSKIP       
  22  PRINT 1200, I,zk(I),(UNUM(K,I),K=K1,K2)                   
      PRINT 1100, (ISPEC(K),K=K1,K2)                       
   9  CONTINUE
C
CZ Moved short lived species printout here
      PRINT 160,N
 160  FORMAT(/1X,"NUMBER DENSITIES OF SHORT-LIVED SPECIES",
     | 1x,"N=",I5)
      NPE = NSP - NQ
      IROW = 10
      LR = NPE/IROW + 1
      RL = FLOAT(NPE)/IROW + 1
      DIF = RL - LR
      IF (DIF.LT.0.001) LR = LR - 1
C
      NQ1 = NQ + 1
      DO 12 L=1,LR
      K1 = NQ1 + (L-1)*IROW
      K2 = K1 + IROW - 1
      IF (L.EQ.LR) K2 = NSP
      PRINT 1102, (ISPEC(K),K=K1,K2)
      DO 24 I=1,NZ,MSKIP
  24  PRINT 129, I,zk(I),(UNUM(K,I),K=K1,K2)
      PRINT 1102, (ISPEC(K),K=K1,K2)
  12  CONTINUE
1102  FORMAT(/"Level",5x,"Z",8X,10(A8,2X)/)


      IF (N.EQ.0) RETURN
      
C      PRINT 1251                                           
C 1251  FORMAT(/1X,"NUMBER DENSITIES OF LONG LIVED SPECIES"/)          
C      PRINT 1101, (ISPEC(K),K=1,NQ)                       
C 1101  FORMAT(/5X,"Z",8X,13(A8,1X))
C 1201  FORMAT(1X,1P10E10.2)

      IROW = 10
      LR = NQ/IROW + 1
      RL = FLOAT(NQ)/IROW + 1
      DIF = RL - LR
      IF (DIF.LT.0.001) LR = LR - 1
C
CZ hydro57 updated below to use JK grid levels
cz      ISKIP = 1
      PRINT 155                                           
 155  FORMAT(/1X,"FLUXES OF LONG-LIVED SPECIES"/)         
cd      ZFL = zk(1)-0.5*DRC(1)*1E-5
cd      ZFT = ZF(NZ)                                        
      DO 10 L=1,LR                                        
      K1 = 1 + (L-1)*IROW                                 
      K2 = K1 + IROW - 1                                  
      IF (L.EQ.LR) K2 = NQ                               
      PRINT 110, (ISPEC(K),K=K1,K2)                       
cd      PRINT 120,ZFL,(FLOW(K),K=K1,K2)                    
      DO 23 I=1,NZ,MSKIP                                  
  23  PRINT 129, I,zk(I),(FLUXKZ(K,I),K=K1,K2)     !hydro57 now prints updated fluxes           
cd      PRINT 120,ZFT,(FUP(K),K=K1,K2)                     
  10  CONTINUE                                            
C                                                         
      PRINT*, "VEFF and UNUM after FLUX printed:"
      PRINT 120, (VEFF(K),K=1,10)
      PRINT 120, (UNUM(K,NZ),K=1,10)
      PRINT 120, (VEFF(K)*UNUM(K,NZ),K=1,10)
      PRINT 120, (FLUXKZ(K,NZ1),K=1,10)
      PRINT 120, (FLUXKZ(K,NZ),K=1,10)
C
C-PK
cz      PRINT 176                                           
cz 176  FORMAT(//1X,"FENG'S LB, MB, TP, TL, FUP, FLOW, CON"/)

       PRINT 178
  178  FORMAT(//1X,"JFK/MEZ  LB, MB, TP, TL, FLUX(NZ), FLUX(1), 
     | CON, CONIN"/)   
      
      DO 13 L=1,LR                                        
      K1 = 1 + (L-1)*IROW                                 
      K2 = K1 + IROW - 1                                  
      IF (L.EQ.LR) K2 = NQ
      PRINT 110, (ISPEC(K),K=K1,K2)                       
C      PRINT 145, (SR(K),K=K1,K2)                          
C      PRINT 145, (PHIDEP(K),K=K1,K2)                       
C      PRINT 145, (TLOSS(K),K=K1,K2)                        
      PRINT 146, (LBOUND1(K),K=K1,K2)
      PRINT 146, (MBOUND(K),K=K1,K2)
 146  FORMAT(13X,13(I1,9X))                               
      PRINT 145                                           
      PRINT 145, (TP(K),K=K1,K2)                          
      PRINT 145, (TL(K),K=K1,K2)                          
      PRINT 145, (FLUXKZ(K,NZ),K=K1,K2)                         
      PRINT 145, (FLUXKZ(K,1),K=K1,K2)  
      PRINT 145, (CONKZ(K),K=K1,K2)      
      PRINT 145, (CONIN(K),K=K1,K2) 

  13  CONTINUE                                            
C
C   CHECK THAT H2 BUDGET IS BALANCED (H2 PRODUCTION = H2 LOSS)
C-PK      H2PROD = 0.
C-PK      DO 70 I=1,NQT
C-PK  70  H2PROD = H2PROD + FLOW(I)*HCOEFF(I)
C-PK      H2ESC = FUP(LH2)
C-PK      HESC = FUP(LH)
C-PK      H2LOSS = 0.5*HESC + H2ESC + H2CHEM
C-PK      H2BAL = H2PROD - H2LOSS
      H2ESC = FUP(LH2) + 0.5*FUP(LH)
      H2BAL = H2VOLC - H2ESC
C-PK Revised hydrogen budget
C      PRINT 310
 310  FORMAT(/1X,'H2 BUDGET BALANCE (PHOTO MODEL)'/)
C      PRINT 311, H2SURF,H2VOLC,H2CHEM,H2ESC,H2BAL
 311  FORMAT(5X,'H2SURF =',1PE10.3,3X,'H2VOLC =',E10.3,3X,'H2CHEM =',
     2  E10.3,4X,'H2 ESCAPE =',E10.3,4X,'H2BAL = ',E10.3/)
C-PK H2 vs. CH4 numbers
C      VPHOTO_H2 = ABS(FLOW(LH2)/UNUM(LH2,1))
C      PRINT 312, FLOW(LH2),VPHOTO_H2
C 312  FORMAT(5X,'PHI(H2) = ',1PE11.4,3X,'VPHOTO(H2) = ',1PE10.3,3X/)
C-PK CO numbers
C      VPHOTO_CO = ABS(FLOW(LCO)/UNUM(LCO,1))
C      PRINT 320, FLOW(LCO),VPHOTO_CO
C 320  FORMAT(5X,'PHI(CO) = ',1PE11.4,3X,'VPHOTO(CO) = ',1PE10.3/)
C
      PRINT 179
      print*, "NR=", NR
  179 FORMAT(/1X,"INTEGRATED REACTION RATES"/)
      PRINT 181
      IROW = 10
      LR = NR/IROW + 1
      RL = FLOAT(NR)/IROW + 1
      DIF = RL - LR
      IF (DIF.LT.0.001) LR = LR - 1
C
      DO 17 L=1,LR
      K1 = 1 + (L-1)*IROW
      K2 = K1 + IROW - 1
      IF (L.EQ.LR) THEN
        K2 = NR
        PRINT 186, K1,(RAT(K),K=K1,K2)
C-TF the following format statement needs to be changed based on
C-TF the number of reactions included in the last line of reaction list
C-TF For example, if there is only one reaction in the last line, then
C-TF "1P1E10.3" should be used. If there are 3, then "1P3E10.3" should 
C-TF be used
  186   FORMAT(I3,2X,1P9E10.3,12X)
        GO TO 17
      ENDIF
      PRINT 180, K1,(RAT(K),K=K1,K2),K2
  180 FORMAT(I3,2X,1P10E10.3,2X,I3)
   17 CONTINUE
      PRINT 181
  181 FORMAT(9X,"1",9X,"2",9X,"3",9X,"4",9X,"5",9X,"6",9X,"7",9X,
     2    "8",9X,"9",8X,"10")
C
      WRITE(140,900)
 900  FORMAT(1X,"Reaction rate profiles"/)
      WRITE(140,903),LR,IROW,NR,NZ
 903  FORMAT(1x,4I6)
      DO 90 L=1,LR
      K1 = 1 + (L-1)*IROW
      K2 = K1 + IROW - 1
      IF (L.EQ.LR) K2 = NR
      WRITE(140,901) (K,K=K1,K2)
 901  FORMAT(/5X,"Z  Rx#",4X,10(I4,7X))
      DO 91 J=1,NZ
  91  WRITE(140,902) zk(J),(RXTOT(K,J),K=K1,K2)
 902  FORMAT(1P11E11.3)
  90  CONTINUE
C

  
C
CZ Added flag and if statement
      IF (IZPRNTFLG12.eq.1) THEN
      PRINT 190
 190  FORMAT(/1X,"ATMOSPHERIC PARAMETERS")
      PRINT 195
 195  FORMAT(1x,"Level",3X,"Z",8X,"T",9X,"EDD",7X,"DEN",8x,"WT",8x,"RHO"
     | ,8x,'D(H)',5x,'D(H2)',5x,'D(O)',5x,'Ha')
      PRINT 200, (I,zk(I),T(I),EDD(I),DEN(I),WTP(I),RHO(I)
     | ,DI(LH,I),DI(LH2,I),DI(LO,I),HSCALE(I)
     | ,I=1,NZ,MSKIP)
 200  FORMAT(1X,I3,1x,1P10E10.2)
      PRINT 195
C
C
C
  224  PRINT 223,N
  223  FORMAT(/1x,"OUTPUT.f: Hydrogen Escape Fluxes",1x,"N=",I5/)
       PRINT 111
  111  FORMAT(/5X,"J",6x,"Z",8X,"R",8x,"RHO",7x,"WT"
     | ,7x,"C(H)",7x,"n(H)",7x,"U",7x,"u(H)",7x,"F(H)"
     | ,7x,"F",7x,"Fmass",5x,"q_adv",5x,"q",9x,"q_cond")
CZ
      ENDIF
CZ     
C-TF   NOTE that the unit of q_v (heat flux density or heat flux per volume
C-TF   and W_grav (gravity work density) are both g cm-1 s-3 (~ergs/s/cm3)
C-TF   q here is the heat flux per unit mass. To convert to heat
C-TF   flux per volume, it needs to be multiplied by RHO.
C
       DO 201 J=1,NZ
       u_h = FLUXKZ(LH,J)/UNUM(LH,J)*(R(NZ)/R(J))**2
       f_mass = RHO(J)*U_ADV(J)*(R(J)/R(NZ))**2
C
  201  PRINT 121, J,zk(J),R(J),RHO(J),WTP(J),USOL(LH,J),UNUM(LH,J)
     | , U_ADV(J),u_h,FLUXKZ(LH,J),f_mass/WT(J),f_mass
C     | , RHO(J)*U_ADV(J)*GZ(J), q(J)*RHO(J)
     | , q_adv(J), q(J), q_cond(J)
  121  FORMAT(1X,I5,1P15E10.2)
       PRINT 111
C
CZ Output column number density of H2 vector to the file idlpass_sp.out
CZ      WRITE(6904,*)USOL !(LH2,1:NZ)

CZ      DO 6905 J=1,NZ
CZ Write species levels to idlpass_sp.out     
CZ      WRITE(6904,3904),J,(USOL(I,J),I=1,NSP)
CZ 3904  FORMAT (I3,1P33E10.3)      
CZ 6905  CONTINUE
CZ

     

C
C 195  FORMAT(/4X,"Z",9X,"T",9X,"EDD",7X,"DEN",8X,"P",8X,"H2OSAT",
C     2  5X,"H2O",7X,"RELH",5X,"CONDEN",4x,"WT")
C      PRINT 200, (zk(I),T(I),EDD(I),DEN(I),P(I),H2OSAT(I),UH2O(I),
C     2  RELH(I),CONDEN(I),WT(I),I=1,NZ,ISKIP)

C      PRINT 200, (zk(NZ),T(NZ),EDD(NZ),DEN(NZ),P(NZ),H2OSAT(NZ),UH2O(NZ),
C     2  RELH(NZ),CONDEN(NZ),WT(I),I=NZ,NZ,1)
C 200  FORMAT(1X,1P10E10.3)
C
      RETURN
      END
C-PK **********************************
