c***********************************************************************
c  NCB subroutine for new 2D model                                     *
c  David B. Considine, 3/11/94                                         *
c  This version of the ncb subroutine is the result of my modification *
c  of Anne Douglass' ncb_ratios version over a long period of time,    *
c  adapted for use in the new version of the 2D model                  *
c***********************************************************************

       subroutine noxclx(IJ,IK)
       
       include 'com2d.h'
       save


c common blocks for aircraft nox perturbations (dbc 1/26/91)

       common/subsonic/subnox(l$,z$)
       common/supersonic/supnox(l$,z$)

       DIMENSION bcclo24(l$,z$),bcn2o5(l$,z$),cno2(l$,z$),
     *           bcbro24(l$,z$)
       DIMENSION RHOD(l$,z$),THED(l$,z$),
     *   OMEGAD(l$,z$),OMEGAN(l$,z$)
       dimension gamloc(l$,z$),ano2init(l$,z$),acloinit(l$,z$)
       DIMENSION zz1(10,l$,z$),Z2(z$),Z3(z$)
       dimension yy(5,l$,z$),xx(4,l$,z$),rr(7,l$,z$)
       dimension alpha(l$,z$),beta(l$,z$)
       dimension isw2(l$,z$),isw(l$,z$),icount(l$,z$),
     c	isw3(l$,z$),isw4(l$,z$)


c **********************************************************************
c                                                                      
c   this part of the program is for part sunlight, part dark           
c                                                                      
c **********************************************************************
       entry ncdn(ij,ik)
        if(ij.eq.13.and.ik.eq.13)then
          iicount=iicount+1
cccef          write(6,*)iicount
        endif

       instep=24
       deltatl=(1.-tfd(ij))*dt/instep

       if(isw2(ij,ik) .eq. 1) then
              c(64,ij,ik)=0.
              cn(64,ij,ik)=0.
              isw(ij,ik)=0
              c(13,ij,ik)=cn(13,ij,ik)
              c(14,ij,ik)=cn(14,ij,ik)
              c(12,ij,ik)=cn(12,ij,ik)
c switch is only set for ij=1,2 and 17,18 so you don't need a second if
c may not need this part anymore
              if(ij .le. 3) then
                     c(6,ij,ik)=c(6,ij+2,ik)
                     c(28,ij,ik)=c(28,ij+2,ik)
              else
                     c(6,ij,ik)=c(6,ij-2,ik)
                     c(28,ij,ik)=c(28,ij-2,ik)
              endif
              isw2(ij,ik)=0
              isw3(ij,ik)=1
              isw4(ij,ik)=1
       endif
c sum aircraft perturbations and get in correct units (dbc 1/26/91)
       noxjet=(subnox(ij,ik)+supnox(ij,ik))/(deltaz(ij,ik)*1.e5)

c heterogeneous reactions are defined in a separate subroutine
c
c beta is the total "interchangeable odd chlorine"; 
c this does not contain
c HCl which will be transported separately
c
       beta(ij,ik)=c(33,ij,ik)-c(25,ij,ik)-c(29,ij,ik)
c

c Bry calculation: define "interchangeable odd bromine," alpha.  Essentially,
c you are dividing Bry into 3 boxes: HBr, which is assumed to remain
c constant over the course of the day, Bromine which is available
c for conversion to BrONO2, and Bromine which is not converted.

      if(isw4(ij,ik).eq.1)then

c This is the first time back after polar night.  HOBr and BrCl
c have built up over the night, so they should be included in alpha.
c also adjust alpha for any changes in Bry that may have occurred
c in the dynamics module:
	
         isw4(ij,ik)=0
         alpha(ij,ik)=(c(48,ij,ik)-c(46,ij,ik))*cn(48,ij,ik)/c(48,ij,ik)

      else

c On any other day, assume that today's daytime av. HOBr and BrCl
c will be similar to yesterday's value.  So, they should not be
c included in alpha:
         
         alpha(ij,ik)=(c(48,ij,ik)-c(46,ij,ik)-c(68,ij,ik)-c(60,ij,ik))
     c                *cn(48,ij,ik)/c(48,ij,ik)

	 if(alpha(ij,ik).le.0.)alpha(ij,ik)=(c(44,ij,ik)+c(47,ij,ik)+
     c                         c(45,ij,ik))*cn(48,ij,ik)/c(48,ij,ik)


      endif


c
c details for boundary conditions
c need N2O5 and ClONO2 at the start of the nighttime calculation

c *** N2O5 dawn
       alamda=(j(9,ij,ik)+k(31,ij,ik)*m(ij,ik))*tfd(ij)*dt
       alamda=alamda+(kh(3,ij,ik)*
     c (c(15,ij,ik))
     c +kh(4,ij,ik)*c(29,ij,ik))*tfd(ij)*dt
       q1=exp(-1.*alamda)
       bcn2o5(ij,ik)=c(8,ij,ik)*alamda/(1.-q1)

cdbc 10/26/95: when alamda -> 0, the quantity alamda/(1-q1) should -> 1. 
c              this doesn't happen because denominator -> 0 numerically.
c              Trap for this situation:

       if(q1.eq.1)bcn2o5(ij,ik)=c(8,ij,ik)


c *** CLONO2 dawn
       O4J=(J(15,ij,Ik)+j(58,ij,ik)+K(13,ij,Ik)*C(13,ij,Ik))
     c              +C(1,ij,Ik)*K(35,ij,Ik)
     c              +c(29,ij,ik)*kh(1,ij,ik)
     c         +(c(15,ij,ik))*kh(2,ij,ik)
       O4GD= K(32,ij,Ik)*M(ij,Ik)*c(6,ij,Ik)
       o4jgd= O4J+O4GD
       y=o4gd*(beta(ij,ik)-2.*c(61,ij,ik))

cdbc 10/2/92:
cdbc The production rate of ClONO2 can't be so large as to exhaust
cdbc all the NO2 available for conversion: Therefore, the rate
cdbc should be the minimum of y and [NO2]/(daytime):

c      maxprod=c(6,ij,ik)/(dt*tfd(ij))
c      if(y.gt.maxprod)then
c         write(6,*)ij,ik,y*dt*tfd(ij)/m(ij,ik),c(6,ij,ik)/m(ij,ik)
c         y=maxprod
c      endif
cdbc back to your regularly scheduled algorithm ...

       Q2=1.-exp(-o4jgd*tfd(ij)*dt)
       bcclono2(ij,ik)=(c(30,ij,ik)-y/o4jgd)*(o4jgd*tfd(ij)*dt)
     c              /q2 + y/o4jgd

cdbc  10/26/95 This can occasionally predict a dawn value of clono2 that is
c     less than zero.  This occurs after polar night when ClONO2
c     is small and the environment changes suddenly such that the
c     equilibrium P/L value is large.  This equation assumes that
c     the previous day's environment is nearly the same as today's,
c     so c(30) and y/o4jgd (P/L) are nearly the same.  I think just
c     setting the dawn value to zero in this case will be fine.

       if(bcclono2(ij,ik).lt.0.)bcclono2(ij,ik)=0.

cdbc 10/26/95 The term o4jgd*tfd*dt/q2 has functional form of x/(1-e^-x)
c    the limit as x -> 0 of this is 1 - but I am getting infinities
c    out of the code here.  so, make following adjustment for values
c    of x small enough to blow up code. Note other adjustments for
c    dusk values of clono2, n2o5, and brono2 as well

       if(q2.eq.0.)bcclono2(ij,ik)=c(30,ij,ik)

c Get dusk values of Brx species from daytime average values by
c assuming constant production and loss of Brx during the day:
c First get dawn value, then get dusk from dawn value.

       brj=j(28,ij,ik)+j(64,ij,ik)+(c(15,ij,ik))*kh(6,ij,ik)

       brgd=k(96,ij,ik)*m(ij,ik)*c(6,ij,ik)
       brjgd=brj+brgd
       z=brgd*alpha(ij,ik)
       q3=1.-exp(-brjgd*tfd(ij)*dt)

c dawn value:

       bcbrono2(ij,ik)=(c(47,ij,ik)-z/brjgd)*(brjgd*tfd(ij)*dt)
     c            /q3 + z/brjgd

c dusk value:

       abr11=brgd*alpha(ij,ik)
       bbr11=brjgd
       brono2d=abr11/bbr11-(abr11/bbr11-bcbrono2(ij,ik))
     c              *exp(-bbr11*tfd(ij)*dt)

c assume the rest of alpha is bro to get dusk value of bro:

      brod=alpha(ij,ik)-brono2d

c if daytime production of brono2 has been overestimated by the
c above approximation, brono2d may be larger than alpha and bro
c will be negative with the above approx.  So, trap for this
c situation:

      if(brod.lt.0.)then
         brod=0.
         brono2d=alpha(ij,ik)
      endif


c need boundaries for the start of the nighttime calculation
c the "d" stands for "dusk"
c clono2
c       a11=o4gd*(beta(ij,ik)-2.*c(61,ij,ik))
       a11=y
       b11=o4jgd
       clono2d=a11/b11-(a11/b11-bcclono2(ij,ik))
     c              *exp(-b11*tfd(ij)*dt)

c n2o5
       n2o5d=bcn2o5(ij,ik)*q1


c no2
       no2d=c(5,ij,ik)+c(6,ij,ik)+2.*(c(8,ij ,ik)-n2o5d)
     c	+(c(30,ij,ik)-clono2d)+(c(47,ij,ik)-brono2d)

       if(no2d.lt.0)no2d=0.
c no3

cdbc I am following Anne's lead here with this boundary condition

       no3=c(7,ij,ik)

c boundary for clo at the start of night
c need to approximate cl2o2; implicitly assumed that any cl2
c goes immediately to cl, clo

cdbc 10/26/95: changing initialization of clo, cl2o2:

c       cl2o2=k(102,ij,ik)*c(28,ij,ik)*c(28,ij,ik)*m(ij,ik)/
c     c       (j(36,ij,ik)+k(106,ij,ik)*m(ij,ik))
c       clod=beta(ij,ik)-clono2d-2.*cl2o2

        quada=m(ij,ik)*k(102,ij,ik)/(j(36,ij,ik)+k(106,ij,ik)*m(ij,ik))
        quadb=0.5
        quadc=0.5*(beta(ij,ik)-clono2d)

      if(j(36,ij,ik).eq.0)then
        cl2o2=c(61,ij,ik)
        clod=beta(ij,ik)-clono2d-2.*cl2o2
      else
        clod=(-1.*quadb+sqrt(quadb*quadb+4.*quada*quadc))/(2.*quada)
        cl2o2=quadc-quadb*clod
        if(cl2o2.lt.0)cl2o2=quada*clod*clod
      endif
        if(clod.lt.0.)clod=0.

cdbc end change in clo initialization 

c need approximate value for HCl just to get the nighttime loss of 
c ClONO2 correctly

       hcld=c(29,ij,ik)

c values of HOCl and HOBr at start of night.  Do same as with
c HCl:

       hocld=c(25,ij,ik)
       hobrd=c(68,ij,ik)

c do nighttime calculation
c these are the quantities I need for 24 hour averaged productions and
c losses for transported quantities

       pn2o5sum=0.
       ln2o5sum=0.
       pclono2sum=0.
       lclono2sum=0.
       lhclsum=0.
       phno3sum=0.
       poddnsum=0.
       loddnsum=0.
       pclxsum=0.
       lclxsum=0.

c need average nighttime values of ClONO2 and N2O5 to get the day
c to night ratios

       aclono2=0.
       an2o5=0.
       ano3=0.

c You want to calculate the average nighttime BrONO2 to get
c day/night ratios:

      abrono2=0.
      abro=0.


       do 100 it1=1,instep

	clod=clod/(1.+(k(32,ij,ik)*no2d+2.*k(102,ij,ik)*clod)
     c		*m(ij,ik)*deltatl)              
	pclono2=k(32,ij,ik)*no2d*clod*m(ij,ik)	
	lclono2=kh(1,ij,ik)*hcld
     c     +kh(2,ij,ik)*(c(15,ij,ik))
	clono2d=(clono2d+pclono2*deltatl)/(1.+lclono2*deltatl)
	pclono2sum=pclono2sum+pclono2/instep
	lclono2sum=lclono2sum+lclono2*clono2d/instep

	aclono2=aclono2+clono2d/instep

        hocld=(hocld+kh(2,ij,ik)*c(15,ij,ik)*clono2d*deltatl)/
     c           (1.+kh(5,ij,ik)*hcld*deltatl)

c no3
	pno3=k(10,ij,ik)*c(4,ij,ik)*no2d+k(31,ij,ik)*n2o5d*m(ij,ik)
	lno3=k(46,ij,ik)*no2d*m(ij,ik)
	no3=(no3+pno3*deltatl)/(1.+lno3*deltatl)
	ano3=ano3+no3/instep		

cdbc	no2d=no2d/(1.+(2.*k(10,ij,ik)*(c(4,ij,ik)+c(1,ij,ik))
cdbc     c		+k(32,ij,ik)*clod*m(ij,ik)
cdbc     c		+k(96,ij,ik)*brod*m(ij,ik))*deltatl)

c no2
	pno2=k(31,ij,ik)*n2o5d*m(ij,ik)
	lno2=(k(46,ij,ik)*no3+k(32,ij,ik)*clod+k(96,ij,ik)*brod)
     c  *m(ij,ik)
	no2d=(no2d+pno2*deltatl)/(1.+lno2*deltatl)

c Bro:  Br is assumed to be zero at night.  This rests on Br+O3
c being much faster than any of the BrO+ClO or BrO + BrO reactions
c converting BrO to Br.  Therefore, no loss of BrO to Br should be
c included in the integration of nighttime BrO, because no production
c from Br is included, either. BrO is therefore simply lost to BrONO2
c and BrCl:


         brod=brod/(1.+(k(105,ij,ik)*clod
     c                +k(96,ij,ik)*no2d*m(ij,ik))*deltatl)

         brono2d=(brono2d+k(96,ij,ik)*no2d*brod*m(ij,ik)*deltatl)/
     c           (1.+kh(6,ij,ik)*(c(15,ij,ik))*deltatl)

	 abrono2=abrono2+brono2d/instep
         abro=abro+brod/instep

         hobrd=(hobrd+kh(6,ij,ik)*c(15,ij,ik)*brono2d*deltatl)/
     c           (1.+kh(7,ij,ik)*hcld*deltatl)


cdbc	pn2o5=k(10,ij,ik)*(c(4,ij,ik)+c(1,ij,ik))*no2d
	pn2o5=k(46,ij,ik)*no3*no2d*m(ij,ik)
	ln2o5=kh(3,ij,ik)*(c(15,ij,ik))
     c        +kh(4,ij,ik)*hcld+k(31,ij,ik)*m(ij,ik)
	n2o5d=(n2o5d+pn2o5*deltatl)/(1.+ln2o5*deltatl)
	pn2o5sum=pn2o5sum+pn2o5/instep
	ln2o5sum=ln2o5sum+ln2o5*n2o5d/instep

	an2o5=an2o5+n2o5d/instep

	hcld=hcld/(1.+(kh(1,ij,ik)*clono2d+kh(4,ij,ik)*n2o5d
     c                 +kh(5,ij,ik)*hocld+kh(7,ij,ik)*hobrd)*deltatl)

	lhclsum=lhclsum+(kh(1,ij,ik)*clono2d+kh(4,ij,ik)*n2o5d
     c                 +kh(5,ij,ik)*hocld+kh(7,ij,ik)*hobrd)
     c			/instep

	phno3sum=phno3sum+(kh(1,ij,ik)*clono2d*hcld
     c			        +kh(4,ij,ik)*n2o5d*hcld
     c	+kh(3,ij,ik)*n2o5d*(c(15,ij,ik))*2.
     c	+kh(2,ij,ik)*clono2d*(c(15,ij,ik))
     c  +kh(6,ij,ik)*brono2d*c(15,ij,ik) )
     c				/instep
	loddnsum=loddnsum+(kh(1,ij,ik)*clono2d*hcld
     c		+kh(4,ij,ik)*n2o5d*hcld
     c		+kh(3,ij,ik)*n2o5d*
     c          (c(15,ij,ik))*2.
     c		+kh(2,ij,ik)*clono2d*
     c          (c(15,ij,ik)))/instep

	poddnsum=poddnsum+kh(4,ij,ik)*n2o5d*hcld/instep

100        continue

	bcbrono2(ij,ik)=brono2d

c approximate night/day ratios
              n2o5day=n2o5d*(1.-q1)/alamda
              rnday(8,ij,ik)=an2o5/n2o5day

              clono2day=y/o4jgd+(clono2d-y/o4jgd)
     c              *q2/(dt*tfd(ij)*o4jgd)
              rnday(30,ij,ik)=aclono2/clono2day

c Daytime average value of BrONO2 from new dawn value and night/day
c ratio:

	cn(47,ij,ik)=z/brjgd+(bcbrono2(ij,ik)-z/brjgd)
     c		*q3/(dt*tfd(ij)*brjgd)

	rnday(47,ij,ik)=abrono2/cn(47,ij,ik)

c transport ClONO2 and calculate the daytime value

c this is the 24 hour ClONO2
	prod=tfd(ij)*(k(32,ij,ik)*m(ij,ik))*c(6,ij,ik)
     c       *(beta(ij,ik)-2.*c(61,ij,ik))
	loss=tfd(ij)*(o4jgd+kh(1,ij,ik)*c(29,ij,ik)+kh(2,ij,ik)*
     c		(c(15,ij,ik)))
                                                  ! This may not be the complete ClONO2 prod/loss - EF, 2/97
        CTPROD(26,IJ,IK) = prod
        CTLOSS(26,IJ,IK) = loss

	cn(30,ij,ik)=(cn(30,ij,ik)+(prod+(1.-tfd(ij))
     c		*(pclono2sum-lclono2sum))*dt)/
     c		(1.+loss*dt/(tfd(ij)+(1.-tfd(ij))*rnday(30,ij,ik)))

c	if(ik .eq. 1) cn(30,ij,ik)=bval(30,ij)*m(ij,1)

c correct back to daytime average
	cn(30,ij,ik)=cn(30,ij,ik)/(tfd(ij)+(1.-tfd(ij))*rnday(30,ij,ik))

        if(cn(30,ij,ik).gt.0.9*cn(31,ij,ik))then
           cn(30,ij,ik)=0.9*cn(31,ij,ik)
c           write(6,*)'clono2 adjustment at',iicount,ij,ik
        endif

c transport N2O5 and calculate the daytime value

c this is the 24 hour N2O5
	prod=tfd(ij)*(k(46,ij,ik)*m(ij,ik))*c(6,ij,ik)*c(7,ij,ik)
	loss=(k(31,ij,ik)*m(ij,ik)+j(9,ij,ik)+kh(3,ij,ik)*c(15,ij,ik)
     c        +kh(4,ij,ik)*c(29,ij,ik))*tfd(ij)
                                                  ! This may not be the complete N2O5 prod/loss - EF, 2/97
        CTPROD(27,IJ,IK) = prod
        CTLOSS(27,IJ,IK) = loss

	cn(8,ij,ik)=(cn(8,ij,ik)+(prod+(1.-tfd(ij))
     c		*(pn2o5sum-ln2o5sum))*dt)/
     c		(1.+loss*dt/(tfd(ij)+(1.-tfd(ij))*rnday(8,ij,ik)))

c	if(ik .eq. 1) cn(8,ij,ik)=bval(8,ij)*m(ij,1)
c correct back to daytime average

	cn(8,ij,ik)=cn(8,ij,ik)/(tfd(ij)+(1.-tfd(ij))*rnday(8,ij,ik))

c transport HNO3
        hno3p=(k(18,ij,ik)*m(ij,ik))*c(13,ij,ik)
     c   *c(6,ij,ik)*tfd(ij)
c add production from heterogeneous processes
        hno3p=hno3p+kh(1,ij,ik)*cn(30,ij,ik)*c(29,ij,ik)*tfd(ij)
     c   +kh(2,ij,ik)*cn(30,ij,ik)
     c   *(c(15,ij,ik))*tfd(ij)
     c   +kh(3,ij,ik)*cn(8,ij,ik)
     c   *(c(15,ij,ik))*tfd(ij)*2.
     c   +kh(4,ij,ik)*cn(8,ij,ik)*c(29,ij,ik)*tfd(ij)
     c   +kh(6,ij,ik)*cn(47,ij,ik)
     c   *(c(15,ij,ik))*tfd(ij)
     c   +phno3sum*(1.-tfd(ij))
        CTPROD(14,IJ,IK) = hno3p       

        hno3l=(k(37,ij,ik)*c(13,ij,ik)+
     c      j(6,ij,ik))*tfd(ij)
     c      +2.4E-7*c(38,ij,ik)
        CTLOSS(14,IJ,IK) = hno3l

        cn(10,ij,ik)=(cn(10,ij,ik) + hno3p*dt)/(1. + hno3l*dt)

        if(ik .eq. 1) cn(10,ij,ik)=bval(10,ij)*m(ij,1)

c ClONO2/NO2
      zz1(2,ij,ik)=cn(30,ij,ik)/c(6,ij,ik)

c N2O5/NO2
      zz1(3,ij,ik)=cn(8,ij,ik)/c(6,ij,ik)

c HNO3/NO2
      zz1(6,ij,ik)=cn(10,ij,ik)/c(6,ij,ik)
      zz1(6,ij,ik)=amin1(cn(10,ij,ik)/c(6,ij,ik),1.e12)
c	if(zz1(6,ij,ik) .eq. 1.e12)type *,ij,ik,cn(10,ij,ik),c(6,ij,ik)

C  HO2NO2(day)/NO2(DAY)=ho2no2(24)/no2(day)
       zz1(1,ij,Ik)=K(34,ij,Ik)*C(14,ij,Ik)*M(ij,Ik)
     c       /(2.4E-7*C(38,IJ,IK)+(J(24,ij,Ik)+K(55,ij,Ik)*
     c       C(13,ij,Ik)+j(42,ij,ik))+K(70,ij,Ik)*M(ij,Ik))

c need to approximate BrO for first time back after polar night
	c(44,ij,ik)=alpha(ij,ik)-cn(47,ij,ik)

C  NO3(day)/NO2(DAY)
       T1=K(44,ij,Ik)*M(ij,Ik)*C(1,ij,Ik)+K(10,ij,Ik)*
     1   C(4,ij,Ik)+(J(24,ij,Ik)*zz1(1,ij,Ik)
     2 +J(15,ij,Ik)*zz1(2,ij,Ik))
       T2=K(56,ij,Ik)*C(21,ij,Ik)
       T3=K(35,ij,Ik)*C(1,ij,Ik)*zz1(2,ij,Ik)+
     1       K(13,ij,Ik)*C(13,ij,Ik)*zz1(2,ij,Ik)
       T4=k(31,ij,ik)*m(ij,ik)*zz1(9,ij,ik)+j(9,ij,ik)*zz1(9,ij,ik)
       B11=(J(5,ij,Ik)+J(17,ij,Ik))+k(74,ij,ik)*c(5,ij,ik)
       B22=k(46,ij,ik)*m(ij,ik)*c(6,ij,ik)
       zz1(5,ij,Ik)=(T1+T2+T3+T4)/(B11+B22)

C  NO/NO2(DAY)
       T1=J(7,ij,Ik)+J(17,ij,Ik)*zz1(5,ij,Ik)+K(43,ij,Ik)*C(1,ij,Ik)
       B11=K(9,ij,Ik)*C(4,ij,Ik)+K(28,ij,ik)*C(28,ij,Ik)
       B22=K(38,ij,Ik)*C(14,ij,Ik)+K(15,ij,Ik)*C(21,ij,Ik)
     c       +k(74,ij,ik)*zz1(5,ij,ik)*c(6,ij,ik)
       zz1(4,ij,Ik)=T1/(B11+B22)

    
C  NO2(DAY)/NOX

C  NOX=NO+NO2+NO3+CLONO2+HO2NO2+2*N2O5


cdbc 10/7/92 for ncb_ratios10.f, I am partitioning NOx', where
cdbc NOx'=NOx-ClONO2.  Accordingly, I want zz7=NO2/NOx', so I
cdbc am taking zz2 out of the calculation.

cdbc old way:

cdbc       zz1(7,ij,Ik)=1.0/(1.+zz1(4,ij,Ik)+zz1(5,ij,ik)+zz1(1,ij,Ik)
cdbc     c	+cn(47,ij,ik)/c(6,ij,ik)
cdbc     c	+zz1(2,ij,ik)+2.*zz1(3,ij,ik)+zz1(6,ij,ik))

cdbc new way:

       zz1(7,ij,Ik)=1.0/(1.+zz1(4,ij,Ik)+zz1(5,ij,ik)+zz1(1,ij,Ik)
     c	+cn(47,ij,ik)/c(6,ij,ik)
     c	+2.*zz1(3,ij,ik)+zz1(6,ij,ik))

c ***** transport odd n

cdbc 10/7/92:  With the redefinition of zz7 from zz7=NO2/NOx to 
cdbc zz7=NO2/NOx', where NOx' = NOx - ClONO2, rnox has to be
cdbc modified: rnox = (NO/NO2)*(NO2/NOx')*(NOx'/NOx).  So instead of:

cdbc       rnox=zz1(4,ij,ik)*zz1(7,ij,ik)

cdbc first approximate NOx'/NOx:

       noxpnox=1.- cn(30,ij,ik)/c(31,ij,ik)

cdbc then incorporate this in rnox:

       rnox=zz1(4,ij,ik)*zz1(7,ij,ik)*noxpnox

       gno=rnox*c(31,ij,ik)
       anum=2.*k(47,ij,ik)*j(16,ij,ik)*tfd(ij)*rnox**2
       bnum=k(47,ij,ik)*rnox
       cnum=k(20,ij,ik)*c(3,ij,ik)
              
       lossn=(2.*anum*c(31,ij,ik)/(bnum*c(31,ij,ik)
     c       +cnum)-anum*(c(31,ij,ik))**2*bnum/
     c       (bnum*c(31,ij,ik)+cnum)**2
     c       +2.4e-7*c(38,ij,ik)*(zz1(1,ij,ik)+zz1(6,ij,ik))
     c      *zz1(7,ij,ik)*noxpnox) 
        CTLOSS(2,IJ,IK) = lossn       
 
       prodn=(2.*k(45,ij,ik)*c(2,ij,ik)*c(11,ij,ik)*tfd(ij)
     c        +anum*(c(31,ij,ik))**2/
     c        (bnum*c(31,ij,ik)+cnum))
     c        -anum*(c(31,ij,ik))**2*bnum*c(31,ij,ik)
     c        /(bnum*c(31,ij,ik)+cnum)**2
        CTPROD(2,IJ,IK) = prodn       


       if(ik .eq. 1) then
              cn(31,ij,ik)=(bval(31,ij)+bval(10,ij))*m(ij,1)

       else
C ADD LIGHTNING PRODUCTION OF NOX  12/15/88  H2 SOURCE OF KO ET AL. 1986
              PLIGHTN=0.0
              IF(IK.GE.2 .AND. IK.LE.8)THEN
                     IF(IJ.GE.7 .AND. IJ.LE.12)THEN
                         PLIGHTN=1.E3
                     ENDIF
              ENDIF
C ADD GCR PRODUCTION OF NOX  12/13/88
              cn(31,ij,ik) = (cn(31,ij,ik)+(prodn+1.25*GCR(IJ,IK)+
     C                        PLIGHTN+noxjet)*dt)/(1.+lossn*dt)
       endif
c ***** calculate easy daytime species concentrations from ratios
c don't forget about BrONO2

cdbc 10/7/92: subtract cn(30,ij,ik) from cn(31,ij,ik) to remove ClONO2
cdbc from partitioning of odd nitrogen:
c       if(ij.eq.15.and.ik.eq.11)write(6,*)ij,ik,zz1(1,ij,ik),
c     c zz1(2,ij,ik),
c     c zz1(3,ij,ik),zz1(4,ij,ik),zz1(5,ij,ik),zz1(6,ij,ik),zz1(7,ij,ik),
c     c cn(31,ij,ik),cn(30,ij,ik),cn(10,ij,ik),cn(8,ij,ik),
c     c rnday(8,ij,ik),rnday(30,ij,ik),prodn,lossn,anum,bnum,
c     c cnum,aclono2,clono2day,c(6,ij,ik),cn(6,ij,ik),y,o4j,o4gd,o4jgd,
c     c y/o4jgd,q2,clono2d,beta(ij,ik),c(33,ij,ik)

       cn(6,ij,Ik)=zz1(7,ij,Ik)*(Cn(31,Ij,Ik)-cn(30,ij,ik))
       cn(32,ij,Ik)=zz1(1,ij,Ik)*cn(6,ij,Ik)
       cn(5,ij,Ik)=zz1(4,ij,Ik)*cn(6,ij,Ik)
       cn(7,ij,Ik)=zz1(5,ij,Ik)*cn(6,ij,Ik)
       cn(8,ij,ik)=zz1(3,ij,ik)*cn(6,ij,ik)
c       cn(30,ij,ik)=zz1(2,ij,ik)*cn(6,ij,ik)
       cn(10,ij,ik)=zz1(6,ij,ik)*cn(6,ij,ik)


c calculate new value for night-to-day ratio for odd N


c as below, this value for bcn2o5 wijl only be used for the next
c timestep puts this latitude into polar night
       bcn2o5(ij,ik)=bcn2o5(ij,ik)*cn(31,ij,ik)/c(31,ij,ik)
       rnday(5,ij,ik)=0.0
C  N(DAY)
       T1=J(16,ij,Ik)*Cn(5,ij,Ik)
       B11=K(47,ij,Ik)*Cn(5,ij,Ik)+K(20,ij,Ik)*C(3,ij,Ik)+
     c       k(68,ij,ik)*c(13,ij,ik)
       Cn(9,ij,Ik)=(T1/B11)


c *** calculate easy chlorine ratios
c $$ influence of BrO on ClO ratios?


c HCl/Cl

      top=(k(26,ij,ik)*c(18,ij,ik)+k(23,ij,ik)
     c      *c(17,ij,ik)+k(53,ij,ik)*c(14,ij,ik)
     c      +k(63,ij,ik)*c(23,ij,ik)
     c      +k(7,ij,ik)*c(16,ij,ik))*tfd(ij)

      bot=(k(27,ij,ik)*c(13,ij,ik)+j(18,ij,ik))*tfd(ij)
     c     +2.4E-7*c(38,ij,ik)
     c  +(1.-tfd(ij))*lhclsum
     c     +tfd(ij)*(kh(1,ij,ik)*cn(30,ij,ik)
     c     +kh(4,ij,ik)*cn(8,ij,ik)
     c     +kh(5,ij,ik)*c(25,ij,ik)
     c     +kh(7,ij,ik)*c(68,ij,ik))

c      yy(5,ij,ik)=top/bot

      yy(5,ij,ik)=top/bot*(1.-exp(-1.*bot*dt))+yy(5,ij,ik)
     c*exp(-1.*bot*dt)
      
c add Cl2O2

      yy(4,ij,ik)=k(102,ij,ik)*c(28,ij,ik)*m(ij,ik)/
     c     (j(36,ij,ik)+k(106,ij,ik)*m(ij,ik))

cdbc  10/26/95: This is a kluge to deal with the possibility
c               that ncdn is called even though j(36) = 0.
c               I am assuming no loss of cl2o2 and production
c               over time dt from ClO + ClO

      if(j(36,ij,ik).eq.0.)
     c yy(4,ij,ik)=c(61,ij,ik)/c(28,ij,ik)+
     c             k(102,ij,ik)*c(28,ij,ik)*m(ij,ik)*dt

      if(c(28,ij,ik).eq.0)yy(4,ij,ik)=0.

cdbc end kluge

      if(isw3(ij,ik) .eq. 0)then

C (HOCL/CLO)  DAY

         yy(1,ij,Ik)=K(6,ij,Ik)*C(14,ij,Ik)/(J(23,ij,Ik)
     c        +K(62,ij,Ik)*C(13,ij,Ik)+kh(5,ij,ik)*c(29,ij,ik))

C (CL/CLO)  DAY

cdbc ratio modified 4/2/92 to include ClO+BrO -> Cl + Br + O2
cdbc as part of the partitioning.

         t1=K(25,ij,Ik)*C(1,ij,ik)+K(28,ij,Ik)*C(5,ij,Ik)
     c     +K(6,ij,Ik)*C(14,ij,Ik)+(k(93,ij,ik)+k(103,ij,ik)
     c     +k(105,ij,ik))*c(44,ij,ik)
cdbc     c    +K(6,ij,Ik)*C(14,ij,Ik)
         t2=(K(60,ij,Ik)*C(13,ij,Ik)+K(61,ij,Ik)*C(13,ij,Ik))
         b11=K(24,ij,Ik)*C(4,ij,Ik)
         yy(2,ij,Ik)=(T1+t2)/B11

      else
C (HOCL/CLO)  DAY

         yy(1,ij,Ik)=K(6,ij,Ik)*Cn(14,ij,Ik)/(J(23,ij,Ik)
     c        +K(62,ij,Ik)*Cn(13,ij,Ik))

C (CL/CLO)  DAY
         t1=K(25,ij,Ik)*Cn(1,ij,ik)+K(28,ij,Ik)*Cn(5,ij,Ik)
     c        +K(6,ij,Ik)*Cn(14,ij,Ik)+(k(93,ij,ik)
     c        +k(103,ij,ik)+k(105,ij,ik))*c(44,ij,ik)
cdbc     c     +K(6,ij,Ik)*Cn(14,ij,Ik)

cdbc c(44 and cn(44 should be the same at this point in the code
cdbc I think, so there is no particular reason to use one or the other.

         t2=(K(60,ij,Ik)*Cn(13,ij,Ik)+K(61,ij,Ik)*Cn(13,ij,Ik))
         b11=K(24,ij,Ik)*C(4,ij,Ik)
         yy(2,ij,Ik)=(T1+t2)/B11
         isw3(ij,ik)=0.

      endif


c calculate clx production and loss
c include production from heterogeneous processes

cdbc 10/8/92: to be consistent with the way the odd nitrogen
cdbc partitioning is done, yy(3)=ClO/(Cl+HCl+ClO+HOCL+2xCl2O2)
cdbc so ClONO2 is not included.  However, the production and
cdbc loss of Clx includes ClONO2, so a correction has to be
cdbc applied to the loss term.

      clxpclx=1.-cn(30,ij,ik)/c(33,ij,ik)

       clxl=0.0
       clxprd1= c(36,ij,ik)*j(19,ij,ik)*tfd(ij)*4.0+
     c 	(J(20,ij,Ik)+K(16,ij,Ik)*c(13,ij,ik))*tfd(ij)*c(37,ij,ik)+
     c	J(21,ij,Ik)*c(34,ij,ik)*tfd(ij)*2.+
     c	c(35,ij,ik)*j(22,ij,ik)*2.*tfd(ij)
      clxprd2=(K(75,IJ,IK)*c(13,ij,ik)+J(26,IJ,IK))*TFD(IJ)
     c  *C(40,IJ,IK)*3.+c(51,ij,ik)*j(31,ij,ik)*tfd(ij)
     c  +c(52,ij,ik)*(j(32,ij,ik)+k(98,ij,ik)*c(13,ij,ik))*tfd(ij)
       clxprd3=c(53,ij,ik)*(j(33,ij,ik)+K(99,IJ,IK)*C(2,IJ,IK))
     c  *tfd(ij)*3.
     c  +c(54,ij,ik)*(j(34,ij,ik)+K(100,IJ,IK)*C(2,IJ,IK))*tfd(ij)*2.
     c  +c(55,ij,ik)*(j(35,ij,ik)+K(101,IJ,IK)*C(2,IJ,IK))*tfd(ij)

c Include proper Clx production from O(1D) reactions with CFC's and HALONs (EF, May 1997)

       clxprd4 = (c(34,ij,ik)*c(2,ij,ik)*k(83,IJ,IK)*2. +
     c            c(35,ij,ik)*c(2,ij,ik)*k(80,IJ,IK)*2. +
     c            c(36,ij,ik)*c(2,ij,ik)*k(54,IJ,IK)*4. +
     c            c(51,ij,ik)*c(2,ij,ik)*k(112,IJ,IK) +
     c            c(57,ij,ik)*c(2,ij,ik)*k(65,IJ,IK) +
     c            c(57,ij,ik)*j(43,ij,ik))*tfd(ij)

	yy(3,ij,ik)=1./(1.+yy(1,ij,ik)+yy(2,ij,ik)
     c		+2.*yy(4,ij,ik)+yy(5,ij,ik)*yy(2,ij,ik))

        clxl=2.4e-7*c(38,ij,ik)*yy(5,ij,ik)*yy(3,ij,ik)*yy(2,ij,ik)
     c  *clxpclx

        CTLOSS(3,IJ,IK) = clxl

        clxprd=clxprd1+clxprd2+clxprd3+clxprd4
        CTPROD(3,IJ,IK) = clxprd


       cn(33,ij,ik)=(Cn(33,IJ,IK) + clxprd*dt)/(1. + clxl*dt)

c       if(ik .eq. 1) cn(33,ij,ik)=bval(33,ij)*m(ij,1)


        cn(28,ij,ik)=yy(3,ij,ik)*(cn(33,ij,ik)-cn(30,ij,ik))


c calculate easy odd chlorine species from easy ratios
       cn(29,ij,ik)=cn(28,ij,ik)*yy(5,ij,ik)*yy(2,ij,ik)
       cn(25,ij,ik)=cn(28,ij,ik)*yy(1,ij,ik)
       cn(27,ij,ik)=cn(28,ij,ik)*yy(2,ij,ik)
       cn(61,ij,ik)=cn(28,ij,ik)*yy(4,ij,ik)
c       if(ik .eq. 1)cn(29,ij,ik)=bval(29,ij)*m(ij,1)

       bcclo24(ij,ik)=cclo*cn(33,ij,ik)/c(33,ij,ik)
       bcclono2(ij,ik)=clono2*cn(33,ij,ik)/c(33,ij,ik)

c DAY/NIGHT ODD BROMINE:
      
c approximate initial Brx

      brx=cn(48,ij,ik)
      brxcount=0
      brxmax=30
      getout=1.e-3

c Proper Brx production from O(1D)-HALON reactions included (EF, May 1997)

      brxprd=c(49,ij,ik)*(j(29,ij,ik)+k(97,ij,ik)*c(13,ij,ik))*tfd(ij)
     c       + c(50,ij,ik)*(j(30,ij,ik))*tfd(ij)
     c       + c(51,ij,ik)*(j(31,ij,ik))*tfd(ij)
     c       + c(51,ij,ik)*c(2,ij,ik)*k(112,IJ,IK)*tfd(ij)
     c       + c(50,ij,ik)*c(2,ij,ik)*k(113,IJ,IK)*tfd(ij)

      CTPROD(15,IJ,IK) = brxprd 

      
c begin repeat ... until loop:

 1065 continue

c BrONO2/BrO:

      if(brxcount.eq.0)then

         rr(1,ij,ik)=k(96,ij,ik)*cn(6,ij,ik)*m(ij,ik)/
     c               (j(28,ij,ik)+j(64,ij,ik)+kh(6,ij,ik)*
     c                cn(15,ij,ik))

      else

         rr(1,ij,ik)=cn(47,ij,ik)/bro

      endif

c HOBr/BrO:

      rr(2,ij,ik)=(k(108,ij,ik)*cn(14,ij,ik)+kh(6,ij,ik)*
     c            cn(15,ij,ik)*rr(1,ij,ik))/(j(62,ij,ik)+kh(7,ij,ik)*
     c            cn(29,ij,ik))

cccef      if(ij.eq.12.and.ik.eq.10)write(6,*)rr(2,ij,ik),k(109,ij,ik),
cccef     ccn(14,ij,ik),kh(6,ij,ik),cn(15,ij,ik),rr(1,ij,ik),j(62,ij,ik),
cccef     ckh(7,ij,ik),cn(29,ij,ik)

c BrCl/BrO:

      if(j(40,ij,ik).gt.0.)then
      rr(3,ij,ik)=(k(105,ij,ik)*cn(28,ij,ik)+kh(7,ij,ik)*cn(29,ij,ik)*
     c             rr(2,ij,ik))/j(40,ij,ik)
      else
      rr(3,ij,ik)=0.
      endif

c HBr/Br:

      rr(4,ij,ik)=(k(92,ij,ik)*cn(14,ij,ik)+k(109,ij,ik)*cn(23,ij,ik))/
     c            (k(95,ij,ik)*cn(13,ij,ik))

c BrO: Now solve a quadratic equation to get BrO:

      aa1=k(107,ij,ik)*cn(1,ij,ik)+j(27,ij,ik)+k(69,ij,ik)*cn(5,ij,ik)
     c    +(k(93,ij,ik)+k(103,ij,ik))*cn(28,ij,ik)
     c    +j(28,ij,ik)*rr(1,ij,ik)+j(62,ij,ik)*rr(2,ij,ik)
     c    +j(40,ij,ik)*rr(3,ij,ik)

      bb1=2.*k(94,ij,ik)*(1.+rr(4,ij,ik))/(k(91,ij,ik)*cn(4,ij,ik))

      bb2=1.+rr(1,ij,ik)+rr(2,ij,ik)+rr(3,ij,ik)+aa1*(1.+rr(4,ij,ik))/
     c    (k(91,ij,ik)*cn(4,ij,ik))

      bro=(-dfloat(bb2)+
     c dsqrt(dfloat(bb2)*dfloat(bb2)
     c +4.*dfloat(bb1)*dfloat(brx)))/(2.*bb1)

c Br/BrO:

      rr(5,ij,ik)=(brx/bro-(1.+rr(1,ij,ik)+
     c            rr(2,ij,ik)+rr(3,ij,ik)))/(1.+rr(4,ij,ik))



c BrO/Brx:

      rr(6,ij,ik)=bro/brx
      if(rr(6,ij,ik).gt.1.)rr(6,ij,ik)=1.

c HBr/Brx:

      rr(7,ij,ik)=rr(4,ij,ik)*rr(5,ij,ik)*rr(6,ij,ik)
      if(rr(7,ij,ik).gt.1.)rr(7,ij,ik)=1.

      brxl=2.4e-7*c(38,ij,ik)*rr(7,ij,ik)*tfd(ij)*dt

      CTLOSS(15,IJ,IK) = brxl/dt

c  Note CTLOSS for Bry is in proper units   

      brxnew=(cn(48,ij,ik) + brxprd*dt)/(1. + brxl)

cccef      if(ij.eq.12.and.ik.eq.10)
cccef     c  write(6,*)ij,ik,brxcount,rr(1,ij,ik),rr(2,ij,ik),rr(3,ij,ik),
cccef     crr(4,ij,ik),rr(5,ij,ik),rr(6,ij,ik),rr(7,ij,ik),aa1,bb1,
cccef     cbb2,bro,brx,brxnew,brxprd,j(28,ij,ik),j(64,ij,ik),j(40,ij,ik),
cccef     cj(62,ij,ik)

      brxcount=brxcount+1

      if(brxcount.gt.brxmax)go to 1066

      if(abs(brxnew-brx).gt.getout)then
         brx=brxnew  
         go to 1065
      else
         go to 1067
      endif
    
 1066 continue
cccef      write(6,*)'Failed convergence in Brx day/night at: ',
cccef     c iicount,ij,ik
 1067 continue

      cn(48,ij,ik)=brxnew
      cn(44,ij,ik)=rr(6,ij,ik)*cn(48,ij,ik)
      cn(46,ij,ik)=rr(7,ij,ik)*cn(48,ij,ik)
      cn(45,ij,ik)=rr(5,ij,ik)*cn(44,ij,ik)
      cn(60,ij,ik)=rr(3,ij,ik)*cn(44,ij,ik)
      cn(68,ij,ik)=rr(2,ij,ik)*cn(44,ij,ik)


      rnday(44,ij,ik)=abro/cn(44,ij,ik)

c      do 328 i=1,s$
c      if(cn(i,ij,ik).ne.cn(i,ij,ik))then
c         write(6,*)'NaN trapped at day/night: ',i,ij,ik
c      endif 
c 328  continue

c      do 3281 i=1,s$
c      if(cn(i,ij,ik).lt.0.)then
c         write(6,*)'Neg value, day/night at: ',i,ij,ik,cn(i,ij,ik)
c      endif 
c 3281 continue


        return 

c ************
c 24 hour day
c note that this is set up to reach this point only after a day-night 
c calculation has been done, i.e., no zero initial values
c ************
c there is no need to make any changes here - except to zero cl2o2.
c there wijl be no heterogeneous processes because it is too warm
        entry ncday(ij,ik)
c sum aircraft perturbations and get in correct units (dbc 1/26/91)
       noxjet=(subnox(ij,ik)+supnox(ij,ik))/(deltaz(ij,ik)*1.e5)
c *** calculate N2O5 and ClONO2 

c calculate ClONO2
	beta(ij,ik)=c(30,ij,ik)+c(28,ij,ik)+c(27,ij,ik)

	O4J=(J(15,ij,Ik)+j(58,ij,ik)+K(13,ij,Ik)*C(13,ij,Ik))
     c           +C(1,ij,Ik)*K(35,ij,Ik)
	O4GD= K(32,ij,Ik)*M(ij,Ik)*c(6,ij,Ik)
	o4jgd= O4J+O4GD
	y=o4gd*beta(ij,ik)

	prod=k(32,ij,ik)*m(ij,ik)*c(6,ij,ik)*beta(ij,ik)
	loss=o4jgd+kh(1,ij,ik)*c(29,ij,ik)
     c       +kh(2,ij,ik)
     c       *(c(15,ij,ik))

        CTPROD(26,IJ,IK) = prod
        CTLOSS(26,IJ,IK) = loss

	cn(30,ij,ik) = (cn(30,ij,ik) + prod*dt)/(1. + loss*dt)

c	if(ik .eq. 1) cn(30,ij,ik)=bval(30,ij)*m(ij,1)

	rnday(30,ij,ik)=1.

c calculate N2O5
	pn2o5=k(46,ij,ik)*c(6,ij,ik)*c(7,ij,ik)*m(ij,ik)
	ln2o5=(j(9,ij,ik)+k(31,ij,ik)*m(ij,ik))
	ln2o5=ln2o5+(kh(3,ij,ik)*
     c              (c(15,ij,ik))
     c        +kh(4,ij,ik)*c(29,ij,ik))

        CTPROD(27,IJ,IK) = pn2o5
        CTLOSS(27,IJ,IK) = ln2o5

	cn(8,ij,ik) = (cn(8,ij,ik) + pn2o5*dt)/(1. + ln2o5*dt)

	rnday(8,ij,ik)=1.


c transport HNO3
        hno3p=(k(18,ij,ik)*m(ij,ik))*c(13,ij,ik)
     c   *c(6,ij,ik)*tfd(ij)
     c   +kh(3,ij,ik)*c(8,ij,ik)
     c   *(c(15,ij,ik))*2.
     c	+kh(1,ij,ik)*c(30,ij,ik)*c(29,ij,ik)
     c	+kh(2,ij,ik)*c(30,ij,ik)
     c  *(c(15,ij,ik))
     c	+kh(6,ij,ik)*c(47,ij,ik)
     c  *(c(15,ij,ik))
     c	+kh(4,ij,ik)*c(8,ij,ik)*c(29,ij,ik)
        CTPROD(14,IJ,IK) = hno3p       

        hno3l=(k(37,ij,ik)*c(13,ij,ik)+
     c      j(6,ij,ik))*tfd(ij)
     c      +2.4E-7*c(38,ij,ik)
        CTLOSS(14,IJ,IK) = hno3l

        cn(10,ij,ik) = (cn(10,ij,ik) + hno3p*dt)/(1. + hno3l*dt)
        if(ik .eq. 1) cn(10,ij,ik)=bval(10,ij)*m(ij,1)

c calculate BrONO2 in an analagous manner to old ClONO2

 	 alpha(ij,ik)=c(44,ij,ik)+c(47,ij,ik)

         brj=j(28,ij,ik)+j(64,ij,ik)
         brgd=k(96,ij,ik)*m(ij,ik)*c(6,ij,ik)
         brjgd=brj+brgd

          z=brgd*alpha(ij,ik)
          cn(47,ij,ik)=z/brjgd+(c(47,ij,ik)-z/brjgd)*
     c                     exp(-brjgd*dt)

C  HO2NO2(day)/NO2(DAY)=ho2no2(24)/no2(day)
         zz1(1,ij,Ik)=K(34,ij,Ik)*C(14,ij,Ik)*M(ij,Ik)*tfd(ij)
     c       /(2.4E-7*C(38,IJ,IK)+(J(24,ij,Ik)+K(55,ij,Ik)*
     c       C(13,ij,Ik)+j(42,ij,ik))*tfd(ij)+K(70,ij,Ik)*M(ij,Ik))

C  NO3(day)/NO2(DAY)
         T1=K(44,ij,Ik)*M(ij,Ik)*C(1,ij,Ik)+K(10,ij,Ik)*
     1    C(4,ij,Ik)+(J(24,ij,Ik)*zz1(1,ij,Ik)+J(15,ij,Ik)
     2    *zz1(2,ij,Ik))
         T2=K(56,ij,Ik)*C(21,ij,Ik)
         T3=K(35,ij,Ik)*C(1,ij,Ik)*zz1(2,ij,Ik)+
     1       K(13,ij,Ik)*C(13,ij,Ik)*zz1(2,ij,Ik)
         T4=k(31,ij,ik)*m(ij,ik)*zz1(9,ij,ik)+j(9,ij,ik)*zz1(9,ij,ik)
         B11=(J(5,ij,ik)+J(17,ij,Ik))+k(74,ij,ik)*c(5,ij,ik)
         B22=k(46,ij,ik)*m(ij,ik)*c(6,ij,ik)
         zz1(5,ij,Ik)=(T1+T2+T3+T4)/(B11+B22)

C  NO/NO2(DAY)
         T1=J(7,ij,Ik)+J(17,ij,Ik)*zz1(5,ij,Ik)+K(43,ij,Ik)*C(1,ij,Ik)
         B11=K(9,ij,Ik)*C(4,ij,Ik)+K(28,ij,ik)*C(28,ij,Ik)
         B22=K(38,ij,Ik)*C(14,ij,Ik)+K(15,ij,Ik)*C(21,ij,Ik)
     c       +k(74,ij,ik)*zz1(5,ij,ik)*c(6,ij,ik)
         zz1(4,ij,Ik)=T1/(B11+B22)
c N2O5/NO2
	zz1(3,ij,ik)=cn(8,ij,ik)/c(6,ij,ik)
c ClONO2/NO2
	zz1(2,ij,ik)=cn(30,ij,ik)/c(6,ij,ik)
c BrONO2/NO2
	zz1(9,ij,ik)=cn(47,ij,ik)/c(6,ij,ik)
c HNO3/NO2

      zz1(6,ij,ik)=cn(10,ij,ik)/c(6,ij,ik)
      zz1(6,ij,ik)=amin1(cn(10,ij,ik)/c(6,ij,ik),1.e12)
c	if(zz1(6,ij,ik) .eq. 1.e12)type *,ij,ik,cn(10,ij,ik),c(6,ij,ik)

C  NO2(DAY)/NOX

C  NOX=NO+NO2+NO3+HNO3+ClONO2+HO2NO2+2*N2O5
cdbc 10/7/92 for ncb_ratios10.f, I am partitioning NOx', where
cdbc NOx'=NOx-ClONO2.  Accordingly, I want zz7=NO2/NOx', so I
cdbc am taking zz2 out of the calculation.

cdbc old way:

cdbc         zz1(7,ij,ik)=1.0/(1.+zz1(5,ij,ik)+zz1(1,ij,ik)+zz1(4,ij,ik)
cdbc     c	+zz1(2,ij,ik)+2.*zz1(3,ij,ik)+zz1(9,ij,ik)+zz1(6,ij,ik))


cdbc new way:

         zz1(7,ij,ik)=1.0/(1.+zz1(5,ij,ik)+zz1(1,ij,ik)+zz1(4,ij,ik)
     c	+2.*zz1(3,ij,ik)+zz1(9,ij,ik)+zz1(6,ij,ik))

cdbc 10/7/92:  With the redefinition of zz7 from zz7=NO2/NOx to 
cdbc zz7=NO2/NOx', where NOx' = NOx - ClONO2, rnox has to be
cdbc modified: rnox = (NO/NO2)*(NO2/NOx')*(NOx'/NOx).  So instead of:

cdbc       rnox=zz1(4,ij,ik)*zz1(7,ij,ik)

cdbc first approximate NOx'/NOx:

       noxpnox=1.- cn(30,ij,ik)/(c(31,ij,ik))

cdbc then incorporate this in rnox:

       rnox=zz1(4,ij,ik)*zz1(7,ij,ik)*noxpnox

         gno=rnox*c(31,ij,ik)
         anum=2.*k(47,ij,ik)*j(16,ij,ik)*tfd(ij)*rnox**2
         bnum=k(47,ij,ik)*rnox
         cnum=k(20,ij,ik)*c(3,ij,ik)
         lossn=(2.*anum*c(31,ij,ik)/(bnum*c(31,ij,ik)
     C        +cnum)-anum*(c(31,ij,ik))**2*bnum/
     c        (bnum*c(31,ij,ik)+cnum)**2
     c        +2.4e-7*c(38,ij,ik)*(zz1(1,ij,ik)+zz1(6,ij,ik))*
     c         zz1(7,ij,ik)*noxpnox)
        CTLOSS(2,IJ,IK) = lossn       

         prodn=(2.*k(45,ij,ik)*c(2,ij,ik)*c(11,ij,ik)*tfd(ij)
     c        +anum*(c(31,ij,ik))**2/
     c        (bnum*c(31,ij,ik)+cnum))
     c        -anum*(c(31,ij,ik))**2*bnum*c(31,ij,ik)
     c        /(bnum*c(31,ij,ik)+cnum)**2
        CTPROD(2,IJ,IK) = prodn       

       if(ik .eq. 1) then
              cn(31,ij,ik)=(bval(31,ij)+bval(10,ij))*m(ij,1)
       else
C ADD LIGHTNING PRODUCTION OF NOX  12/15/88  H2 SOURCE OF KO ET AL. 1986
              PLIGHTN=0.0
              IF(IK.GE.2 .AND. IK.LE.8)THEN
		                  IF(IJ.GE.7 .AND. IJ.LE.12)THEN
                    PLIGHTN=1.E3
                  ENDIF
              ENDIF
C ADD GCR PRODUCTION OF NOX  12/13/88
              cn(31,ij,ik)=(cn(31,ij,ik)+(prodn+1.25*GCR(IJ,IK)+
     C           PLIGHTN+noxjet)*dt)/(1.+lossn*dt)
       endif

	rnday(31,ij,ik)=1.
c ***** calculate easy daytime species concentrations from ratios
       cn(6,ij,ik)=zz1(7,ij,ik)*(cn(31,ij,ik)-cn(30,ij,ik))
       Cn(32,ij,ik)=zz1(1,ij,ik)*cn(6,ij,ik)
       Cn(5,ij,ik)=zz1(4,ij,ik)*cn(6,ij,ik)
       Cn(7,ij,ik)=zz1(5,ij,ik)*cn(6,ij,ik)
c       cn(30,ij,ik)=zz1(2,ij,ik)*cn(6,ij,ik)
       cn(8,ij,ik)=zz1(3,ij,ik)*cn(6,ij,ik)
       cn(10,ij,ik)=zz1(6,ij,ik)*cn(6,ij,ik)
c       cn(47,ij,ik)=zz1(9,ij,ik)*cn(6,ij,ik)
       rnday(5,ij,ik)=1.
C  N(DAY)
       T1=J(16,ij,Ik)*Cn(5,ij,Ik)
       B11=K(47,ij,Ik)*Cn(5,ij,Ik)+K(20,ij,Ik)*C(3,ij,Ik)+
     c       k(68,ij,ik)*c(13,ij,ik)
       Cn(9,ij,Ik)=(T1/B11)



C (HOCL/CLO)  DAY
       yy(1,ij,Ik)=K(6,ij,Ik)*C(14,ij,Ik)/(J(23,ij,Ik)
     c       +K(62,ij,Ik)*C(13,ij,Ik)+kh(5,ij,ik)*c(29,ij,ik))

c       bet3d(ij,ik)=bet3d(ij,ik)+yy(1,ij,ik)
C (CL/CLO)  DAY
       t1=K(25,ij,Ik)*C(1,ij,ik)+K(28,ij,Ik)*C(5,ij,Ik)
     c    +K(6,ij,Ik)*C(14,ij,Ik)+k(93,ij,ik)*c(44,ij,ik)
       t2=(K(60,ij,Ik)*C(13,ij,Ik)+K(61,ij,Ik)*C(13,ij,Ik))
       b11=K(24,ij,Ik)*C(4,ij,Ik)
       yy(2,ij,Ik)=(T1+t2)/B11

c HCl/Cl

       top=(k(26,ij,ik)*c(18,ij,ik)+k(23,ij,ik)
     c   *c(17,ij,ik)+k(53,ij,ik)*c(14,ij,ik)
     c   +k(63,ij,ik)*c(23,ij,ik)
     c   +k(7,ij,ik)*c(16,ij,ik))*tfd(ij)

      bot=(k(27,ij,ik)*c(13,ij,ik)+j(18,ij,ik))
     c   +2.4E-7*c(38,ij,ik)+kh(1,ij,ik)*c(30,ij,ik)
     c	+kh(4,ij,ik)*c(8,ij,ik)+kh(5,ij,ik)*c(25,ij,ik)
     c  +kh(7,ij,ik)*c(68,ij,ik)

      yy(5,ij,ik)=top/bot

c Cl2O2/ClO

      yy(4,ij,ik)=k(102,ij,ik)*c(28,ij,ik)*m(ij,ik)/
     c            (j(36,ij,ik)+k(106,ij,ik)*m(ij,ik))



c calculate clx production and loss

      clxpclx=1.-cn(30,ij,ik)/c(33,ij,ik)

       clxl=0.
       clxprd=(c(36,ij,ik)*j(19,ij,ik)*tfd(ij)*4.0+(J(20,ij,Ik)
     c       +K(16,ij,Ik)*C(13,ij,Ik))*tfd(ij)*c(37,ij,ik)+J(21,ij,Ik)
     c       *c(34,ij,ik)*tfd(ij)*2.+c(35,ij,ik)*j(22,ij,ik)*2.*tfd(ij)
     C       +(K(75,IJ,IK)*C(13,IJ,IK)+J(26,IJ,IK))*TFD(IJ)*
     C       C(40,IJ,IK)*3.)+c(51,ij,ik)*j(31,ij,ik)*tfd(ij)
     c  +c(52,ij,ik)*(j(32,ij,ik)+k(98,ij,ik)*C(13,IJ,IK))*tfd(ij)
     c  +c(53,ij,ik)*(j(33,ij,ik)+K(99,IJ,IK)*C(2,IJ,IK))*tfd(ij)*3.
     c  +c(54,ij,ik)*(j(34,ij,ik)+K(100,IJ,IK)*C(2,IJ,IK))*tfd(ij)*2.
     c  +c(55,ij,ik)*(j(35,ij,ik)+K(101,IJ,IK)*C(2,IJ,IK))*tfd(ij)

c Include proper Clx production from O(1D) reactions with CFC's and HALONs (EF, May 1997)

     c  + (c(34,ij,ik)*c(2,ij,ik)*k(83,IJ,IK)*2. + 
     c     c(35,ij,ik)*c(2,ij,ik)*k(80,IJ,IK)*2. +
     c     c(36,ij,ik)*c(2,ij,ik)*k(54,IJ,IK)*4. +
     c     c(51,ij,ik)*c(2,ij,ik)*k(112,IJ,IK) +
     c     c(57,ij,ik)*c(2,ij,ik)*k(65,IJ,IK) + 
     c     c(57,ij,ik)*j(43,ij,ik))*tfd(ij)

        CTPROD(3,IJ,IK) = clxprd

       yy(3,ij,ik)=1./(1.+yy(1,ij,ik)+yy(2,ij,ik)+2.*yy(4,ij,ik)
     c           +yy(5,ij,ik)*yy(2,ij,ik))
     
        clxl=2.4e-7*c(38,ij,ik)*yy(5,ij,ik)*yy(3,ij,ik)*yy(2,ij,ik)
     c  *clxpclx

        CTLOSS(3,IJ,IK) = clxl

       cn(33,ij,ik) = (cn(33,IJ,IK) + clxprd*dt)/(1. + clxl*dt)
     
c       if(ik .eq. 1)cn(33,ij,ik)=bval(33,ij)*m(ij,1)

       cn(28,ij,ik)=yy(3,ij,ik)*(cn(33,ij,ik)-cn(30,ij,ik))

c calculate easy odd chlorine species from easy ratios

       cn(25,ij,ik)=cn(28,ij,ik)*yy(1,ij,ik)
       cn(27,ij,ik)=cn(28,ij,ik)*yy(2,ij,ik)
       cn(61,ij,ik)=cn(28,ij,ik)*yy(4,ij,ik)
       cn(29,ij,ik)=cn(28,ij,ik)*yy(5,ij,ik)*yy(2,ij,ik)
       rnday(33,ij,ik)=1.

c DAYTIME ODD BROMINE:

c approximate initial Brx

      brx=cn(48,ij,ik)
      brxcount=0
      brxmax=30
      getout=1.e-3

c Proper Brx production from O(1D)-HALON reactions included (EF, May 1997)

      brxprd=c(49,ij,ik)*(j(29,ij,ik)+k(97,ij,ik)*c(13,ij,ik))
     c       +c(50,ij,ik)*(j(30,ij,ik))
     c       +c(51,ij,ik)*(j(31,ij,ik))
     c       + c(51,ij,ik)*c(2,ij,ik)*k(112,IJ,IK)
     c       + c(50,ij,ik)*c(2,ij,ik)*k(113,IJ,IK)

      CTPROD(15,IJ,IK) = brxprd 


c begin repeat ... until loop:

 1055 continue
c      write(6,*)ij,ik,brxcount,rr(1,ij,ik),rr(2,ij,ik),rr(3,ij,ik),
c     crr(4,ij,ik),rr(5,ij,ik),rr(6,ij,ik),rr(7,ij,ik),aa1,bb1,
c     cbb2,bro,brx,brxprd,j(28,ij,ik),j(64,ij,ik),j(40,ij,ik),
c     cj(62,ij,ik)
c BrONO2/BrO:

         rr(1,ij,ik)=k(96,ij,ik)*cn(6,ij,ik)*m(ij,ik)/
     c               (j(28,ij,ik)+j(64,ij,ik)+kh(6,ij,ik)*
     c                cn(15,ij,ik))

c HOBr/BrO:

      rr(2,ij,ik)=(k(108,ij,ik)*cn(14,ij,ik)+kh(6,ij,ik)*
     c            cn(15,ij,ik)*rr(1,ij,ik))/(j(62,ij,ik)+kh(7,ij,ik)*
     c            cn(29,ij,ik))


c BrCl/BrO:

      if(j(40,ij,ik).gt.0.)then
      rr(3,ij,ik)=(k(105,ij,ik)*cn(28,ij,ik)+kh(7,ij,ik)*cn(29,ij,ik)*
     c             rr(2,ij,ik))/j(40,ij,ik)
      else
      rr(3,ij,ik)=0.
      endif

c HBr/Br:

      rr(4,ij,ik)=(k(92,ij,ik)*cn(14,ij,ik)+k(109,ij,ik)*cn(23,ij,ik))/
     c            (k(95,ij,ik)*cn(13,ij,ik))

c BrO: Now solve a quadratic equation to get BrO:

      aa1=k(107,ij,ik)*cn(1,ij,ik)+j(27,ij,ik)+k(69,ij,ik)*cn(5,ij,ik)
     c    +(k(93,ij,ik)+k(103,ij,ik))*cn(28,ij,ik)
     c    +j(28,ij,ik)*rr(1,ij,ik)+j(62,ij,ik)*rr(2,ij,ik)
     c    +j(40,ij,ik)*rr(3,ij,ik)

      bb1=2.*k(94,ij,ik)*(1.+rr(4,ij,ik))/(k(91,ij,ik)*cn(4,ij,ik))

      bb2=1.+rr(1,ij,ik)+rr(2,ij,ik)+rr(3,ij,ik)+aa1*(1.+rr(4,ij,ik))/
     c    (k(91,ij,ik)*cn(4,ij,ik))

      bro=(-dfloat(bb2)+
     c dsqrt(dfloat(bb2)*dfloat(bb2)
     c +4.*dfloat(bb1)*dfloat(brx)))/(2.*bb1)

c Br/BrO:

      rr(5,ij,ik)=(brx/bro-(1.+rr(1,ij,ik)+
     c            rr(2,ij,ik)+rr(3,ij,ik)))/(1.+rr(4,ij,ik))


c BrO/Brx:

      rr(6,ij,ik)=bro/brx
      if(rr(6,ij,ik).gt.1.)rr(6,ij,ik)=1.

c HBr/Brx:

      rr(7,ij,ik)=rr(4,ij,ik)*rr(5,ij,ik)*rr(6,ij,ik)
      if(rr(7,ij,ik).gt.1.)rr(7,ij,ik)=1.

      brxl=2.4e-7*c(38,ij,ik)*rr(7,ij,ik)*dt
      CTLOSS(15,IJ,IK) = brxl/dt
c                                                             Note CTLOSS for Bry is in proper units   
   
      brxnew=(cn(48,ij,ik) + brxprd*dt)/(1. + brxl)

      brxcount=brxcount+1
      
      if(brxcount.gt.brxmax)go to 1056

      if(abs(brxnew-brx).gt.getout)then
         brx=brxnew  
         go to 1055
      else
         go to 1057
      endif
    
 1056 continue
c      write(6,*)'Failed convergence in Brx day at:  ',
c     c iicount,ij,ik
 1057 continue

      cn(48,ij,ik)=brxnew
      cn(44,ij,ik)=rr(6,ij,ik)*cn(48,ij,ik)
      cn(46,ij,ik)=rr(7,ij,ik)*cn(48,ij,ik)
      cn(45,ij,ik)=rr(5,ij,ik)*cn(44,ij,ik)
      cn(60,ij,ik)=rr(3,ij,ik)*cn(44,ij,ik)
      cn(68,ij,ik)=rr(2,ij,ik)*cn(44,ij,ik)
      cn(47,ij,ik)=rr(1,ij,ik)*cn(44,ij,ik)

c      do 329 i=1,s$
c      if(cn(i,ij,ik).ne.cn(i,ij,ik))then
c         write(6,*)'NaN trapped at day: ',i,ij,ik
c      endif 
c 329  continue

c      do 329 i=1,s$
c      if(cn(i,ij,ik).lt.0.)then
c         write(6,*)'Neg value day at: ',i,ij,ik
c      endif 
c 329  continue
c       if(ij.eq.2.and.ik.eq.4)write(6,*)'day'

       return
c *************
c 24 hour night
c note that this is set up to reach this point only after a day-night
c calculation has been done, i.e., no zero initial values
c *************
c there are lots of wintertime changes for heterogeneous chemistry
       entry ncnight(ij,ik)
       instep=96
       deltatl=dt/instep

       if(isw(ij,ik).eq.0)then
              isw(ij,ik)=1
              isw2(ij,ik)=1
c these save the model from destruction if not started at equinox!
              if(bcclo24(ij,ik) .lt. 1.)bcclo24(ij,ik)=c(28,ij,ik)
              if(bcn2o5(ij,ik) .lt. 1.) bcn2o5(ij,ik)=c(8,ij,ik)

cdbc cl2 has to be set = 0 on entry into subroutine.  cl2 is assumed
cdbc to photolyze immediately during daytime.

              c(64,ij,ik)=0.
              cn(64,ij,ik)=0.

        endif

c sum aircraft perturbations and get in correct units (dbc 1/26/91)
       noxjet=(subnox(ij,ik)+supnox(ij,ik))/(deltaz(ij,ik)*1.e5)

c set initial values for nighttime integration

	clod=(c(28,ij,ik)+c(27,ij,ik))
        hcld=c(29,ij,ik)
        clono2d=c(30,ij,ik)
        cl2o2d=c(61,ij,ik)
        hocld=c(25,ij,ik)
	no2d=(c(5,ij,ik)+c(6,ij,ik))
	no3=c(7,ij,ik)
	brono2d=c(47,ij,ik)
	brod=c(44,ij,ik)
	brd=c(45,ij,ik)
        brcld=c(60,ij,ik)
        hobrd=c(68,ij,ik)
        n2o5d=c(8,ij,ik)
        cl2d=c(64,ij,ik)

      if(cn(48,ij,ik).gt.0.and.c(48,ij,ik).gt.0)then
          brox=(brod+brd )*cn(48,ij,ik)/c(48,ij,ik)
       else
          brox=0.
       endif


c set night to day ratios

       rnday(8,ij,ik)=1.
       rnday(30,ij,ik)=1.
       rnday(33,ij,ik)=1.
       rnday(31,ij,ik)=1.


c set short lived radicals to zero
       cn(27,ij,ik)=0.0
       Cn(5,ij,ik)=0.0
       cn(9,ij,ik)=0.0


c HBr fraction stays fixed

      rr(7,ij,ik)=c(46,ij,ik)/c(48,ij,ik)

c transport Brx

       CTPROD(15,IJ,IK) = 0.0

       brxloss=2.4E-7*C(38,IJ,IK)*rr(7,ij,ik)
       CTLOSS(15,IJ,IK) = brxloss

       cn(48,ij,ik) = cn(48,ij,ik)/(1. + brxloss*dt)

       cn(46,ij,ik)=rr(7,ij,ik)*cn(48,ij,ik)



c don't change beta here - there is no sun to change cl2 i
c into cl, clo
c just let hcl go down, that wijl impact on some other things.
c also remember than HNO3 must go up

       pn2o5sum=0.
       ln2o5sum=0.
       pclono2sum=0.
       lclono2sum=0.
       lhclsum=0.
       phno3sum=0.
       phoclsum=0.
       poddnsum=0.
       loddnsum=0.
       pclxsum=0.
       lclxsum=0.

c need average nighttime values of ClONO2 and N2O5 to get the day
c to night ratios
       ahocl=0.
       ahcl=0.
       aclono2=0.
       acl2o2=0.
       acl2=0.
       an2o5=0.
       ano2=0.
       aclo=0.
       abrono2=0.
       abro=0.
       abrd=0.
       abrcl=0.
       ahobr=0.
       ano3=0.

       mnow=m(ij,ik)
c       if(ij.eq.1.and.ik.eq.11)
c     c write(6,*)kh(1,ij,ik),kh(2,ij,ik),
c     c (c(15,ij,ik))

       do 110 it1=1,instep

        cl2o2d=(cl2o2d+k(102,ij,ik)*mnow*clod*clod*deltatl)/
     c         (1.+k(106,ij,ik)*mnow*deltatl)         

        acl2o2=acl2o2+cl2o2d/instep

	clod=(clod+k(106,ij,ik)*mnow*cl2o2d*deltatl)
     c       /(1.+((k(32,ij,ik)*no2d+
     c		2.*k(102,ij,ik)*clod)*mnow+
     c          ((k(93,ij,ik)+k(103,ij,ik)+k(105,ij,ik))*brod))
     c           *deltatl)

	aclo=aclo+clod/instep

	pclono2=k(32,ij,ik)*no2d*clod*mnow	
	lclono2=kh(1,ij,ik)*hcld
     c          +kh(2,ij,ik)*(c(15,ij,ik))

	clono2d=(clono2d+pclono2*deltatl)/(1.+lclono2*deltatl)

	pclono2sum=pclono2sum+pclono2/instep
	lclono2sum=lclono2sum+lclono2/instep

	aclono2=aclono2+clono2d/instep

c NO3
	pno3=k(10,ij,ik)*c(4,ij,ik)*no2d+k(31,ij,ik)*n2o5d*mnow
	lno3=k(46,ij,ik)*no2d*mnow
	no3=(no3+pno3*deltatl)/(1.+lno3*deltatl)
	ano3=ano3+no3/instep		

cdbc	no2d=no2d/(1.+(2.*k(10,ij,ik)*c(4,ij,ik)+
cdbc     c       k(32,ij,ik)*clod*mnow
cdbc     c	     +k(96,ij,ik)*brod*mnow)*deltatl)

c NO2
	pno2=k(31,ij,ik)*n2o5d*mnow
	lno2=(k(46,ij,ik)*no3+k(32,ij,ik)*clod+k(96,ij,ik)*brod)
     c  *mnow+k(10,ij,ik)*c(4,ij,ik)
	no2d=(no2d+pno2*deltatl)/(1.+lno2*deltatl)
	ano2=ano2+no2d/instep

c Br/BrO
	rrr=((k(93,ij,ik)+k(103,ij,ik))*clod
     c      +2.*k(94,ij,ik)*brod)/
     c     (k(91,ij,ik)*c(4,ij,ik))

        hobrd=(hobrd+kh(6,ij,ik)*c(15,ij,ik)*brono2d*deltatl)/
     c           (1.+kh(7,ij,ik)*hcld*deltatl)

        ahobr=ahobr+hobrd/instep

        brcld=brcld+(k(105,ij,ik)*brod*clod
     c               +kh(7,ij,ik)*hcld*hobrd)*deltatl

        abrcl=abrcl+brcld/instep

	brono2d=(brono2d+k(96,ij,ik)*brod*no2d*mnow*deltatl)/
     c           (1.+kh(6,ij,ik)*(c(15,ij,ik))*deltatl)

	broxd=cn(48,ij,ik)-cn(46,ij,ik)-brono2d-hobrd-brcld

        if(broxd.lt.0.)then
           broxd=0.
           brono2d=cn(48,ij,ik)-cn(46,ij,ik)-hobrd-brcld
           if(brono2d.lt.0.)brono2d=0.
        endif

	brod=broxd/(1.+rrr)
	brd=brod*rrr	

	abrono2=abrono2+brono2d/instep
	abrod=abrod+brod/instep
	abrd=abrd+brd/instep

	pn2o5=k(46,ij,ik)*no2d*no3*mnow
cdbc	pn2o5=k(10,ij,ik)*c(4,ij,ik)*no2d
	ln2o5=kh(3,ij,ik)
     c  *(c(15,ij,ik))+kh(4,ij,ik)*hcld
     c  +k(31,ij,ik)*mnow
	n2o5d=(n2o5d+pn2o5*deltatl)/(1.+ln2o5*deltatl)
	
	pn2o5sum=pn2o5sum+pn2o5/instep
	ln2o5sum=ln2o5sum+ln2o5/instep
	an2o5=an2o5+n2o5d/instep

      hcld=hcld/(1.+(kh(1,ij,ik)*clono2d+kh(4,ij,ik)*n2o5d
     c +kh(5,ij,ik)*hocld+kh(7,ij,ik)*hobrd+2.4e-7*c(38,ij,ik))*deltatl)

	lhclsum=lhclsum+(kh(1,ij,ik)*clono2d+kh(4,ij,ik)*n2o5d
     c  +kh(5,ij,ik)*hocld+kh(7,ij,ik)*hobrd+2.4e-7*c(38,ij,ik))/instep

              ahcl=ahcl+hcld/instep

	phno3sum=phno3sum+(kh(1,ij,ik)*clono2d*hcld
     c		+kh(4,ij,ik)*n2o5d*hcld
     c	    	+kh(3,ij,ik)
     c          *(c(15,ij,ik))*2*n2o5d
     c		+kh(2,ij,ik)*clono2d
     c          *(c(15,ij,ik))
     c		+kh(6,ij,ik)*brono2d
     c          *c(15,ij,ik) )/instep

	loddnsum=loddnsum+(kh(1,ij,ik)*clono2d*hcld
     c		+kh(4,ij,ik)*n2o5d*hcld
     c		+kh(3,ij,ik)*n2o5d
     c          *(c(15,ij,ik))*2.
     c		+kh(2,ij,ik)*clono2d
     c          *(c(15,ij,ik)))/instep


	pclxsum=pclxsum+kh(1,ij,ik)*clono2d*hcld/instep+
     c		kh(4,ij,ik)*n2o5d*hcld/instep+
     c		kh(2,ij,ik)*clono2d
     c          *(c(15,ij,ik))/instep

c	lclxsum=lclxsum+k(32,ij,ik)*no2d*clod*mnow/instep

        phocl=kh(2,ij,ik)*clono2d
     c        *(c(15,ij,ik))
        lhocl=kh(5,ij,ik)*hcld
	phoclsum=phoclsum
     c           +kh(2,ij,ik)*clono2d
     c           *(c(15,ij,ik))/instep
        hocld=(hocld+phocl*deltatl)/(1.+lhocl*deltatl)
        ahocl=ahocl+hocld/instep

       cl2d=cl2d+deltatl*(kh(1,ij,ik)*clono2d
     c +kh(5,ij,ik)*hocld)*hcld
       acl2=acl2+cl2d/instep

 110        continue

c       if(ij.eq.1.and.ik.eq.11)write(6,3232)ij,cn(30,ij,ik)/mnow,
c     c c(30,ij,ik)/mnow,pclono2sum*dt/mnow,
c     c lclono2sum*dt/mnow,mnow,dt,aclono2/mnow

        CTPROD(26,IJ,IK) = pclono2sum
        CTLOSS(26,IJ,IK) = lclono2sum

      cn(30,ij,ik) = (cn(30,ij,ik) + pclono2sum*dt)/(1. + lclono2sum*dt)

cdbc 10/26/95: ClONO2 can't be larger than Cly.  So, if it is
c              calculated to exceed Cly, set it to Cly. (Because
c              ClONO2 is subtracted from Cly below to get other
c              family constituents, if this is not done then 
c              negative values for the other family members can
c              be obtained.)

       if(cn(30,ij,ik).ge.cn(33,ij,ik))cn(30,ij,ik)=cn(33,ij,ik)


c       if(ij.eq.1.and.ik.eq.11)write(6,3232)ij,cn(30,ij,ik)/mnow,
c     c c(30,ij,ik)/mnow,pclono2sum*dt/mnow,
c     c lclono2sum*dt/mnow,mnow,dt,difn(26,ij,ik)*dt/mnow



       sum=ahcl+aclo+ahocl+2.*acl2+2.*acl2o2
       hclrat=ahcl/sum
       clorat=aclo/sum
       hoclrat=ahocl/sum
       cl2o2rat=2.*acl2o2/sum
       cl2rat=2.*acl2/sum

       cn(44,ij,ik)=abro
       cn(45,ij,ik)=abrd 
       cn(47,ij,ik)=abrono2
       cn(60,ij,ik)=abrcl
       cn(68,ij,ik)=ahobr

c       if(ij.eq.1.and.ik.eq.3)write(6,*)abro,abrd,abrono2,
c     c abrcl,ahobr,cn(46,ij,ik),cn(48,ij,ik),rrr,brod,brd,
c     c broxd,brono2d,hobrd,brcld

        CTPROD(27,IJ,IK) = pn2o5sum
        CTLOSS(27,IJ,IK) = ln2o5sum

       cn(8,ij,ik) = (cn(8,ij,ik) + pn2o5sum*dt)/(1. + ln2o5sum*dt)


c	cn(8,ij,ik)=an2o5
c        write(6,*)ij,ik,aclo,aclono2,ahcl 	     

c transport HNO3
        hno3p=phno3sum
        hno3l=2.4E-7*c(38,ij,ik)

        CTPROD(14,IJ,IK) = hno3p       
        CTLOSS(14,IJ,IK) = hno3l       

       cn(10,ij,ik) = (cn(10,ij,ik) + hno3p*dt)/(1. + hno3l*dt)

        if(ik .eq. 1) cn(10,ij,ik)=bval(10,ij)*m(ij,1)


       noxprod=0.0
       noxl=0
        CTLOSS(2,IJ,IK) = noxl
        CTPROD(2,IJ,IK) = noxprod

       if(ik .eq. 1) then
             cn(31,ij,ik)=(bval(31,ij)+bval(10,ij))*m(ij,1)
       else
C ADD LIGHTNING PRODUCTION OF NOX  12/15/88  H4 SOURCE OF KO ET AL. 1986
             PLIGHTN=0.0
             IF(IK.GE.2 .AND. IK.LE.8)THEN
                  IF(IJ.GE.7 .AND. IJ.LE.12)THEN
                    PLIGHTN=1.E3
                  ENDIF
             ENDIF
C ADD GCR PRODUCTION OF NOX  12/13/88
             cn(31,ij,ik)=(cn(31,ij,ik)+(noxprod+1.25*GCR(IJ,IK)+
     C             PLIGHTN+noxjet)*dt)/(1.+noxl*dt)
       endif

       cn(32,ij,ik)=c(32,ij,ik)*cn(31,ij,ik)/c(31,ij,ik)
       cn(6,ij,ik)=ano2
	sum=2.*cn(8,ij,ik)+cn(32,ij,ik)+cn(6,ij,ik)+cn(47,ij,ik)
     c                  +cn(30,ij,ik)+cn(10,ij,ik)
	cn(6,ij,ik)=cn(6,ij,ik)*cn(31,ij,ik)/sum
	cn(32,ij,ik)=cn(32,ij,ik)*cn(31,ij,ik)/sum
	cn(8,ij,ik)=cn(8,ij,ik)*cn(31,ij,ik)/sum
c	cn(47,ij,ik)=cn(47,ij,ik)*cn(31,ij,ik)/sum
c 	cn(30,ij,ik)=cn(30,ij,ik)*cn(31,ij,ik)/sum
	cn(10,ij,ik)=cn(10,ij,ik)*cn(31,ij,ik)/sum
	sum2=cn(6,ij,ik)+cn(32,ij,ik)+2.*cn(8,ij,ik)+
     c		cn(47,ij,ik)+cn(30,ij,ik)

 3232   format(i3,7e11.3)
 3233   format(7e11.3)

       clxprd=0.
       clxl=2.4e-7*c(38,ij,ik)
       CTLOSS(3,IJ,IK) = clxl
       CTPROD(3,IJ,IK) = clxprd

       cn(33,IJ,IK) = (cn(33,IJ,IK) + clxprd*dt)/(1. + clxl*dt)

c       if(ik .eq. 1) cn(33,ij,ik)=bval(33,ij)*m(ij,1)

       cn(29,ij,ik)=hclrat*(cn(33,ij,ik)-cn(30,ij,ik))
       cn(28,ij,ik)=clorat*(cn(33,ij,ik)-cn(30,ij,ik))
       cn(25,ij,ik)=hoclrat*(cn(33,ij,ik)-cn(30,ij,ik))
       cn(61,ij,ik)=cl2o2rat/2.*(cn(33,ij,ik)-cn(30,ij,ik))
       cn(64,ij,ik)=cl2rat/2.*(cn(33,ij,ik)-cn(30,ij,ik))

  

c      do 331 i=1,s$
c      if(cn(i,ij,ik).ne.cn(i,ij,ik))then
c         write(6,*)'NaN trapped at night:  ',i,ij,ik
c      endif 
c 331  continue

c      do 331 i=1,s$
c      if(cn(i,ij,ik).lt.0.)then
c         write(6,*)'Neg value at night: ',i,ij,ik
c      endif 
c 331  continue

       return
       end