C23456789012345678901234567890123456789012345678901234567890123456789012
C
      SUBROUTINE DIUR(NCAL,T,DTD,MONS,FIL,F8)
C
C
C     DIURNALLY AVERAGE PRODUCTION AND DESTRUCTION
C
C    *******************************************************************
C    * NR - TOTAL NUMBER OF CHEMICAL REACTIONS                         *
C    * N2 - BIMOLECULAR REACTIONS NUMBER                               *
C    * N3 - THREE MOLECULAR REACTIONS NUMBER                           *
C    * NMM- MONOMOLECULAR REACTIONS NUMBER                             *
C    * NHET-HETEROGENEOUS REACTIONS NUMBER                             *
C    * NCL- PSC REACTIONS NUMBER                                       *
C    * NAER-BACKGROUND AEROSOL REACTIONS NUMBER                        *
C    * L-NUMBER OF PRESSURE LEVELS                                     *
C    * NY-NUMBER OF CALCULATED GASES                                   *
C    *******************************************************************
      INCLUDE 'basp.h'
      INCLUDE 'din.h'     ! L,M,FI,PA,ZKM,AM,TR,KKAL,OPMD,NREL
      INCLUDE 'garr.h'    ! NC,NY,NYPR,PP,DD,PAV,DAV
      INCLUDE 'charr.h'   ! NR,N,PE,E,RO,R1,SN,SB,GAM,ST,RRA,AMU
      INCLUDE 'yarr.h'    ! Y,Y3,P,D
      INCLUDE 'pharr.h'   ! NSS,S,KS,VICTOR'S PARAMETERS
      INCLUDE 'hetarr.h'  ! AR,PSC1,PSC2
      INCLUDE 'wrm.h'     ! NGZ,NZAP
      INCLUDE 'trajarr.h' ! TRTIME,TRSZA,LEVP
      INCLUDE 'fam.h'     ! MS
C
      COMMON /JZIZIT/ JZ,IZ,JTZ
      COMMON /CHMBOX/ YBOX(L,NC),                 ! NLEV,...,1
     *                PBOX(L,NY),                 ! in SPS
     *                DBOX(L,NY)

      DIMENSION R(NR,L),SOL(NSS,L),F8(Len_tm)
      DATA TSUT/86400./,YMIN/1.E-20/,SCIC_MAX/157.5/,SCIC_MIN/12.3/
C
C     CHEMICAL AND HETEROGENEOUS REACTIONS CONSTANTS
C
      DO 614 I=1,L
         AM(I)=YBOX(I,NC)
         CALL SXR(R(1,I),NR,N2,N3,PE,E,TR(I),AM(I),RO,R1,SN,
     *   SB,NHET,YBOX(I,:),PA(I),ST,NMM,RRA,RA,N,NPSC1,NPSC2,
     *   KH2O,NHNO3,NHCL,GAM,AR(I),PSC1(I),PSC2(I),AMU,NST)
 614  CONTINUE
C
C     PHOTODISSOCIATION PARAMETERS FOR FULL DAY
C
      CALL PHOT_P(FLN8,ALT,TR,PA,L,L1,HH,AMDS,AM,O2S,
     *TAS,CONC,MNABS,WLT,Len_tm,Num_beg,Len_X,Ngas,n_abs,xs_all,
     *MNL,GAS,A_GAS,K_T,SIG)
C
C     DAILY STEPS
C
C/04/23/05/      T=0.
C     SS3/00 IT IS NECESSARY TO CALCULATE PHOTOLYSIS RATES AT
C            THE BEGINNING OF EACH DAY BECAUSE TEMPERATURE HAS CHANGED
C
      CALL PHTA(NSS,MNL,Len_tm,L,MNABS,N_O2,N_H2O,GAS,ALT,
     *WBD,we,wlt,UGG,ALB,EXT,SSAM,GM,SIG,CONC,AMDS,F8,HH,YCOL,E0,EPOT,
     *QFCL,EGR,S,xs_all,Num_beg,Len_X,TI,T,FIL,YBOX,KO3,UGMA,IP,NC)
C/04/23/05/      T=T+DTD                        
C
C     PHOTODISSOCIATION CONSTANTS CALCULATION
C
      DO 66 I=1,L
         DO 66 J=1,NSS
  66        SOL(J,I)=S(J,I)
      CALL PHTA(NSS,MNL,Len_tm,L,MNABS,N_O2,N_H2O,GAS,ALT,
     *WBD,we,wlt,UGG,ALB,EXT,SSAM,GM,SIG,CONC,AMDS,F8,HH,YCOL,E0,EPOT,
     *QFCL,EGR,S,xs_all,Num_beg,Len_X,TI,T+DTD,FIL,YBOX,KO3,UGMA,IP,NC)
!      if(FIL.eq.0.)print *,'dtd=',dtd,' t=',t/3600.,' UGG=',UGG,
!     $CONC(L,:)
      NSEP=MAX0(1,INT(DTD/DTF+0.9))
      TN=REAL(NSEP)/DTD
      TLN=1./(TN*TSUT)
      AMTSUT=DTD/TSUT
C      print *,' T=',T,' DTD=',DTD,' Nsep=',NSEP !' O3=',Y(KO3,1)
c      write(77,*)' T=',T/3600.,' O3=',Y(KO3,1),' NO2=',Y(34,1),
c     $' NO=',Y(33,1)
C     $ ,(S(J,1),J=2,3),(SOL(J,1),J=2,3)
      IPRI=0
C      if(FIL.gt.-5.and.FIL.lt.5.and.UGG.lt.60.)IPRI=1
      DO 617 I=1,L
C         write (*,*) 'beg SEM1=',I,L
         CALL SEM(YBOX(I,:),NY,MS,NBS,NFS,NUS)
C         write (*,*) 'end SEM1=',I,L
C
C        PHOTOCHEMICAL PRODUCTION AND DESTRUCTION
C
         CALL PHX(R(1,I),S(1,I),SOL(1,I),YBOX(I,:),NSS,N,NR,
     *   NSEP,TN,NCAL,ITM,DYM,KS,ISEM,I,IPRI)  !ISEM added by GM3/01
C         if(FIL.gt.-5.and.FIL.lt.5.and.UGG.lt.60.)print 2001,
C     $   FIL,UGG,ZKM(i),YBOX(I,30),PP(30),DD(30),S(2,i),YBOX(I,31),
C     $   YBOX(I,NC)
!       if(FIL.eq.0)then
!          print 2001,T/3600.,zkm(i),YBOX(I,2),PP(2),DD(2),
!     $    YBOX(I,30),YBOX(I,KH2O)
!       endif
 2001    format(2f6.1,7(1pe10.2))
C
C        DIURNALLY AVERAGE CONCENTRATION
C
!      if(S(1,I).ne.0..or.SOL(1,I).ne.0)then
!      print *,FIL,zkm(i),YBOX(I,NC-2),S(1,I),SOL(1,I),PP(1),dd(1)
!      pause
!      endif
         DO 617 J=1,NY
C/04/23/05/              D(I,J)=D(I,J)+DD(J)*AMTSUT
C/04/23/05/  617         P(I,J)=P(I,J)+PP(J)*AMTSUT
            PBOX(I,J)=PP(J)
 617        DBOX(I,J)=DD(J)
!      if(FIL.eq.0.)pause
   5  CONTINUE
      close(77)
 309  CONTINUE
   6  RETURN
      END
C
C----------------------------------------------------------
      SUBROUTINE SXR(R,N,M,K,P,S,T,A,B,C,D,E,NHT,X,Y,ST,L,RRA,RA,
     *NR,NPSC1,NPSC2,KH2O,NHNO3,NHCL,GAM,AR,PSC1,PSC2,AMU,NST)
C
C     CONSTANTS FOR CHEMICAL AND HETEROGENEOUS REACTIONS
C
      DIMENSION R(N),P(M),S(M),B(K),C(K),D(K),E(K),ST(3,L),RRA(N),
     *GAM(NHT),AMU(NHT),X(*),NR(2,N),NST(L)
      COMMON/CRNU/KCO,NHN,NHO2,NHO2H
      COMMON /JZIZIT/ JZ,IZ,JTZ
C
C     BIMOLECULAR REACTIONS
C
!      print *,' SXR',X(1:76),AR
      DO 1 JJ=1,M
   1     R(JJ)=P(JJ)*EXP(-S(JJ)/T)
C
C     TRIMOLECULAR REACTIONS
C
      DO 2 JJ=1,K
   2     R(M+JJ)=RTRZ(T,A,B(JJ),C(JJ),D(JJ),E(JJ))
C
C     SPECIAL CASES
C
C     1. OH+HNO3
C
!JPL05      AK3=1.9E-33*EXP(725./T)
!JPL05      AK2=4.1E-16*EXP(1440./T)
         AK3=6.5E-34*EXP(1335./T)
         AK2=2.7E-17*EXP(2199./T)
      IF(NHN.LE.N)R(NHN)=R(NHN)+AK3*A/(1.+AK3*A/AK2)
C
C     2. OH+CO
C
      IF(KCO.LE.N)R(KCO)=RTRZ(T,A,1.5E-13,2.1E9,-0.6,-6.1)/A+ ! JPL2010
     $   RTRZ(T,A,5.9E-33,1.1E-12,1.4,-1.3)
C
C     3. HO2+HO2
C
!jpl05      IF(NHO2.LE.N)R(NHO2)=R(NHO2)+1.7E-33*A*EXP(1000./T)
C
C     4. HO2+HO2+H2O
C
      IF(NHO2H.LE.N)R(NHO2H)=R(NHO2)*1.4E-21*EXP(2200./T)*X(KH2O)
C
C     MONOMOLECULAR REACTIONS
C
      NB=M-L
      IF(NB.GE.M)GOTO6
      DO 4 JJ=1,L
         J=NB+JJ
         RR=EXP(-0.5*S(J)/T)
!jpl05         REQ=RTRZ(T,A,P(J),ST(2,JJ),ST(1,JJ),ST(3,JJ))
!jpl05   4     R(J)=REQ*RR*RR/A
   4     R(J)=R(NST(JJ))*RR/(P(J)*A)*RR               !P(J)*RR*RR/A
C
C     HETEROGENEOUS REACTIONS CONSTANTS
C
   6  NN=M+K
      NSUL=NHT-NPSC1-NPSC2
      CALL RATE_HET(NHT,RA,T,X(KH2O),X(NHNO3),
     *X(NHCL),RRA(NN+1),GAM,AR,AMU)
      DO 3 JJ=1,NSUL
         IF(X(NR(2,JJ+NN)).GT.1.E-10)THEN
            R(JJ+NN)=GAM(JJ)*AR/X(NR(2,JJ+NN))
         ELSE
            R(JJ+NN)=0.
         ENDIF
!         print *,JJ,AR,GAM(JJ),X(NR(2,JJ+NN)),R(JJ+NN)
   3     CONTINUE
      NN=NN+NSUL
      DO 8 JJ=1,NPSC1
C         if (JZ.EQ.22.AND.IZ.EQ.19) write (*,*) ' beg RX8',JJ,JJ+NN
C      if (psc1.gt.1.e-8)then
C	   write (*,*) NR(2,JJ+NN),X(NR(2,JJ+NN)),psc1
C	   pause
C	endif
   8     R(NN+JJ)=GAM(JJ+NSUL)*PSC1/X(NR(2,JJ+NN))
      NN=NN+NPSC1
      DO 9 JJ=1,NPSC2
   9     R(NN+JJ)=GAM(JJ+NSUL+NPSC1)*PSC2/X(NR(2,JJ+NN))
C
C     RAINOUT CONSTANTS
C
      NNN=M+K+NHT+1
      IF(Y.GE.600.)THEN
         RMN=1.
      ELSE
         IF(Y.LE.200.)THEN
            RMN=0.
         ELSE
            RMN=EXP(2.-2.*600./Y)
         ENDIF
      ENDIF
      DO 7 JJ=NNN,N
   7     R(JJ)=RRA(JJ)*RMN
!      print *,Y,R(NNN:N)
!      pause
   5  CONTINUE
      RETURN
      END
C
      FUNCTION RTRZ(T,A,R,C,S,Y)
C
C     TRIMOLECULAR RECTIONS
C
      B=A*R*(300./T)**S
      IF(C.NE.0.)GOTO1
      RTRZ=B
      GOTO2
   1  D=B/(C*(T/300.)**(-Y))
      P=1./(1.+(ALOG10(D))**2)
      F=0.6**P
      RTRZ=B/(1.+D)*F
   2  RETURN
      END
C
C------------------------------------------------------------
C
      SUBROUTINE RATE_HET(NHET,RA,Temp,H2O,HNO3,HCl,RRA,GAM,AR,AMU)
C
C     HETEROGENEOUS REACTION RATES
C
      REAL Na,kBoltz
      DIMENSION GAM(NHET),AMU(NHET),RRA(NHET),WV(100)
      COMMON/NHETR/N1,N2,N3,N4,N5,N6,N7
      DATA
     &         kBoltz / 1.3806e-19 /,   ! Boltzman's const
     &         Na / 6.0225E23 /         ! Avogadro's const, mol^-1
!      print *,Temp,H2O,HnO3,HCL,AR
      DO 1 J=1,NHET
         WV(J)=3636.73*SQRT(TEMP/AMU(J))! Velocities of diff. molecules
   1     GAM(J)=RRA(J)
C
C     Special cases:
C
      PRM0=KBOLTZ*TEMP
      PRM=PRM0*H2O
      HNO3torr=PRM0*HNO3*0.76
      TICE=2663.5/(12.537-ALOG10(PRM*100.))
C
C     CLONO2+HCL (NAT)      HR,JGR,1993,22931
C
      GAM(N7)=AMAX1(1./(1./0.23+EXP(0.518*(TEMP-TICE))/0.7022),0.1)
      IF(AR.LT.1.E-10)GOTO2
C
C >> Aerosol reactions:
C    =================
C     N2O5+H2O     Zhang ET AL, 1995
C
C      IF(TEMP.GT.205.)THEN
C         GAM(N1) = 0.1
C      ELSE
C         GAM(N1)=-0.252-252.54*PRM+(0.001523+1.213*PRM)*TEMP
C      ENDIF
      con_ndpp = kBoltz*Temp/1013.25  ! convert. coef from number dens.
c                                       to partial pressure in atm.
      Va = ra/3.*1.E-7                ! *SA = specific aerosol vol.,cm^-3
C
C     CLONO2+H2O (Sulfate aerosol)
C     CLONO2+HCL (Silfate aerosol)
C
c ===============
c  >> R2 and R3 are estimated simultaneously
c     Hanson and Ravishankara, J.Phys.Chem., v.98, No.22, p.5734, 1994
      pH2O = H2O*kBoltz*Temp                       ! mbar
      pHCl = HCl*kBoltz*Temp/1013.25               ! atm
      psc = 10.**(9.217-2190./(Temp-12.7))
      aH2O = pH2O/psc
      g0 = 1.18e-4+9.1e-3*aH2O+.5*aH2O*aH2O
      HHCl = exp(6250./Temp-10.414)*aH2O**3.49
      Gs = 576.*aH2O*HHCl*pHCl
      P = 2.E3*HHCl*pHCl/aH2O
      Gclc = g0*sqrt(1.+P)
      adivl = ra/(1.4e-2*sqrt(1./aH2O))
      f = 1./tanh(adivl)-1./adivl
      gef = 1./(1./(Gs+f*Gclc)+1./0.3)
      GAM(N3) = gef*(Gs+f*Gclc*P/(1.+P))/(Gs+f*Gclc)
      GAM(N2) = gef - GAM(N3)
C
C     HOCL+HCL (Sulfate aerosol)
C
      z = alog(pH2O)
      w = ((-14.458+.62456*z)*Temp+3565.)/(44.777+1.3204*z-.19988*Temp)
      IF(w.LT.30.)w=30.
      HHCl = 10.**(15.514-0.1791*w)
      HCLstr = HHCl*pHCl
      x = 60.-w
      if ( x .ge. 0.) then
        HD = 15.+3.*x
c HOCl+HCl(a) is the bulk reaction, since l=1.3mkm>raer, so its rate
c is proportional to the volume of aerosol, not surface!
        volHOCl=1.E5*con_ndpp*HClstr*Va*Na
        R4 = 1.E4*HD*volHOCl
        Gef  = R4/WV(4)
        GAM(N4) = amin1(Gef,1.)
      else
c I do not allow gamma3=0. if wH2SO4>60% and fit it to Fig.5(HRS), but it
c seems unimportant for the model results
         GAM(N4) = 10.**(-0.408*w+21.735)
      endif
C
C     BRONO2+H2O   HANS&RAVI,JGR,1996
C
      GRXN=211.*(AH2O)**1.37
      GAM(N5)=1./(1./0.84+1./GRXN)
C
C     BRONO2+HCL   JPL97
C
      GAM(N6)=9./8.*GAM(N5)
C
C     HOBR+HCL
C
C      GAM(N6)=0.
C      H2Otorr=PRM*0.76
C      TMELT=3236./(11.502-ALOG10(PRM*0.76))
C      TSAT=9384./(39.-ALOG10(HNO3torr*H2Otorr))
C      IF(Temp.LT.TSAT.OR.TEMP.LE.TICE)THEN
C         GAM(N1)=0.006
C         GAM(N2)=AMIN1(0.002,10.**(24.659-0.1398*Temp))
C         GAM(N3)=EXP(-0.636-0.4802*(Temp-TICE))
C         GAM(N4)=AMIN1(0.1,GAM(N4))
C         GAM(N5)=0.006
C      ENDIF
      GAM(N3)=AMIN1(GAM(N3),0.23)
   2  DO 3 J=1,NHET
   3     GAM(J)=GAM(J)*WV(J)
      RETURN
      END
C
C-----------------------------------------------------------
      SUBROUTINE PHOT_P(FL,AL,TR,PA,L,L1,HH,AMDS,AM,O2S,TAS,CONC,MNABS,
     *WLT,Len_t,Num_beg,Len_X,Ng,n_abs,xs_all,NL,GAS,A_GAS,K_T,SIG)
C
C
C
      CHARACTER*32 GAS(Ng),A_GAS(MNabs)
      CHARACTER*60 FL
      DIMENSION AL(L),TR(L),PA(L),HH(2),AM(L),AMDS(L),TAS(L),CONC(L,2),
     *WLT(Len_t),xs_all(Len_t,L,Ng),K_T(Ng),SIG(Mnabs,L,NL),O2S(L)
      INTEGER Num_beg(Ng),Len_X(Ng),n_abs(MNabs)
      DATA RGAS/8.31436/, O2MIXR /0.209476/
C
C     GEOMETRIC ALTITUDE
C
      al(1)=0. !29.266e2*TR(1)*alog(1000./pa(1))
      do i=2,l
         al(i)=al(i-1)+29.266e2*0.5*(tr(i)+tr(i-1))*alog(pa(i-1)/pa(i))
      enddo
C
C     PHOTODISSOCIATION PARAMETERS FOR FULL DAY
C
      HH(1) = RGAS*TR(L)/9.8/32.*1.E5
      HH(2) = RGAS*TR(L)/9.8/48.*1.E5
      DO 92 II = 1,L1
         AMDS(II) = 0.5*(AM(II)+AM(II+1))
         O2S(II)=AMDS(II)*O2MIXR
         TAS(II) = 0.5*(TR(II)+TR(II+1))
         CONC(II,1)=O2S(II)
  92     CONTINUE
      CONC(L,1)=O2MIXR*AM(L)
      DO 91 IABS = 1, MNABS

Cgal         call vX_SEC(A_GAS(iabs),L,TAS,Len_t,wlt,Num_beg(n_abs(iabs)),
cgal     *   Len_X(n_abs(iabs)),xs_all(1,1,n_abs(iabs)),K_T(n_abs(iabs)),FL)

Cgal
      GO TO (1,2) , IABS
      write (*,*) 'error in PHOT_P, IABS = ', IABS
      STOP

    1 CONTINUE
Cgal       ( X. eq. 'O2      -->O       +O'' )
      call   vO2_X_97 ( L,Len_t,wlt,Num_beg(n_abs(iabs)),
     *   Len_X(n_abs(iabs)),xs_all(1,1,n_abs(iabs)),K_T(n_abs(iabs)))
      GO TO 3

    2 CONTINUE
Cgal       ( X .eq. 'O3tot' )
      call   vO3tot_X ( L,TAS,Len_t,wlt,Num_beg(n_abs(iabs)),
     *   Len_X(n_abs(iabs)),xs_all(1,1,n_abs(iabs)),K_T(n_abs(iabs)))
      GO TO 3

    3 CONTINUE
Cgal

         DO 91 il = 1, NL
            DO 91 ih = 1,L
  91           sig(iabs,ih,il)=xs_all(il,ih,n_abs(iabs))
      K_T (n_abs(MNabs)) = 1
      call vX_SEC_ALL(GAS,L,TR,Len_t,wlt,Num_beg,Len_X,XS_all,Ng,K_T,FL)
      RETURN
      END
C
C----------------------------------------------------------
      SUBROUTINE PHTA(NR,NL,Len_t,NH,NABS,NO2,NH2O,GAS,ALT,wb,
     *we,wlt,SZA,ALB,EXT,SSAM,GM,SIG,CONC,AIR,F8,HH,YCOL,E0,E,
     *QF,EGR,PJ,xs_all,Num_beg,Len_X,TI,T,FI,Y,KO3,UGMA,IP,NC)
      include 'jpl.par'
      INCLUDE 'trajarr.h'

      REAL  ALT(NH), WLT(Len_t), EXT(NH,NL), SSAM(NH,NL), GM(NH,NL),
     *      SIG(NABS,NH,NL), XS_all(Len_t,NH,NR), CONC(NH,NABS),
     *      HH(NABS), E0(NH,NL), E(NH,NL), QF(NL,NH), EGR(NH,NL),
     *      PJ(NR,NH), F8(Len_t), AIR(NH),
     *      YCOL(NH,NABS),Y(NH,NC),wb(Len_t), we(Len_t)
      integer   Num_beg(NR), Len_X(NR)
      character*32  GAS(NR)

      DO 625 J=1,NR
         DO 625 I=1,NH
 625        PJ(J,I)=0.
C
C     SOLAR ZENITH ANGLE CALCULATION
C
      CALL CSM0(TI,FI,SZA,T)
      IF(SZA.LT.0.OR.SZA.GT.UGMA)THEN
         IP=2
      ELSE
         IP=1
         DO 630 II=1,NH-1
 630        CONC(II,2)=0.5*(Y(II,KO3)+Y(II+1,KO3))
         CONC(NH,2)=Y(NH,KO3)
C234567890123456789012345678901234567890123456789012345678901234567890
C
C        PHOTODISSOCIATION CONSTANTS CALCULATION
C
         CALL PHOTS(NR,NL,Len_t,NH,NABS,NO2,NH2O,GAS,ALT,WB,we,WLT,
     *   SZA,ALB,EXT,SSAM,GM,SIG,CONC,       ! 'TR' NOT USED GM 1/00
     *   AIR,F8,HH,YCOL,E0,E,QF,EGR,PJ,xs_all,Num_beg,Len_X)
      ENDIF
      RETURN
      END
C
C-----------------------------------------------------------
      SUBROUTINE PHX(R,SR,SO,Y,NS,N,NR,NVA,TN,NCL,ITM,DYM,KS,ISM,IZ,
     $IPRI)
C
C     SHORT-LIVED GASES DIURNAL CALCULATIONS
C
      INCLUDE 'garr.h'
      INCLUDE 'chm.h'
      INCLUDE 'fam.h'
      DIMENSION Y(NC),R(*),SR(NS),KS(NS),SO(NS),S(100),N(2,NR),NCK(100)
      DATA YMINC/1.E-3/, YMIN/1.E-20/
      COMMON /JZIZIT/ JZ,IIZ,JTZ
      DO 1 J=1,NY
         PAV(J)=0.
   1     DAV(J)=0.
      DO 78 IVA=1,NVA
         DO 3 J=1,NY
            PP(J)=0.
   3        DD(J)=0.
         DO 77 JJ=1,NCL
            J=NDA(JJ)
  77        Y5(J)=Y(J)*TN
         AMM=(REAL(NVA-IVA))/REAL(NVA)
         AMP=1.-AMM
         SMAX=0.
         DO 66 J=1,NS
            S(J)=SO(J)*AMM+SR(J)*AMP
  66        SMAX=AMAX1(SMAX,S(J))
         IF(SMAX.GT.1.E-20)THEN
            IP=1
         ELSE
            IP=2
         ENDIF
         KCK=NCL
         DO 67 J=1,KCK
  67        NCK(J)=NDA(J)
         ITER=1
  45     KCKR=0
         AMM=0.3  ! 1./(3.+ITER)
         AMP=1.-AMM
         DO 24 JI=1,KCK
            J=NCK(JI)
            P=0.                                        !
            D=0.                                        !
            GOTO(26,27),IP
  26        DO 2 J1=KPR(J),LPR(J)
   2           P=P+R(NPR(J1))*Y(NPYJ1(J1))*Y(NPYJ2(J1))
            DO 5 J1=KPS(J),LPS(J)
C         if(ipri.eq.1.and.j.eq.30)print 121,' NPS ',j1,NPS(J1),
C     $   S(NPS(J1)),Y(KS(NPS(J1))),S(NPS(J1))*Y(KS(NPS(J1))),P
 121     format(a5,2i4,7(1pe10.2))
   5           P=P+S(NPS(J1))*Y(KS(NPS(J1)))
            IF(J-MS)46,46,47
  46        DO 18 J1=KDR(J),LDR(J)
  18           D=D+R(NDR(J1))*Y(NDYJ(J1))*Y(NDYJ2(J1))
            DO 16 J1=KDS(J),LDS(J)
  16           D=D+S(NDS(J1))*Y(KS(NDS(J1)))
            DO 17 J1=KDW(J),LW(J)
  17           D=D+R(NDW(J1))*Y(N(1,NDW(J1)))
            D=D/Y(J)
            GOTO25
  47        DO 8 J1=KDR(J),LDR(J)
C         if(ipri.eq.1.and.j.eq.30)print 121,' NDR ',j1,NDR(J1),
C     $   R(NDR(J1)),Y(NDYJ(J1)),R(NDR(J1))*Y(NDYJ(J1)),D
   8           D=D+R(NDR(J1))*Y(NDYJ(J1)) !*Y(NDYJ2(J1))
            DO 6 J1=KDS(J),LDS(J)
   6           D=D+S(NDS(J1)) !*Y(KS(NDS(J1)))
            GOTO35
  27        DO 28 J1=KPNR(J),LPNR(J)
  28           P=P+R(NPNR(J1))*Y(NPNYJ1(J1))*Y(NPNYJ2(J1))
            IF(J-MS)48,48,49
  48        DO 43 J1=KDNR(J),LDNR(J)
  43           D=D+R(NDNR(J1))*Y(NDNYJ(J1))*Y(NDNYJ2(J1))
            DO 37 J1=KDW(J),LW(J)
  37           D=D+R(NDW(J1))*Y(N(1,NDW(J1)))
            D=D/Y(J)
            GOTO25
  49        DO 33 J1=KDNR(J),LDNR(J)
  33           D=D+R(NDNR(J1))*Y(NDNYJ(J1)) !*Y(NDNYJ2(J1))
  35        DO 7 J1=KDW(J),LW(J)
   7           D=D+R(NDW(J1)) !*Y(N(1,NDW(J1)))
  25        YY=AMP*(P+Y5(J))/(D+TN)+AMM*Y(J)
C      if(D.gt.1.e10)then
C	   print *,j,ay(j),p,d,y(j),IP
C	   do j1=KDNR(J),LDNR(J)
C	      print *,NDNR(J1),R(NDNR(J1)),ay(NDNYJ(J1)),Y(NDNYJ(J1))
C	   enddo
C	   pause
C	endif
C            YY=(0.7*(P+Y5(J))/(D+TN*Y(J))+0.3)*Y(J)
            IF(Y(J).LT.YMINC)GOTO633
            DY=ABS(YY-Y(J))/Y(J)
            IF(DY.LE.DYM)GOTO633
               KCKR=KCKR+1
               NCK(KCKR)=J
 633        PP(J)=P
            DD(J)=D !/Y(J)         ! /05/12/05/  Units 1/s
  24        Y(J)=YY
         KCK=KCKR
         IF(KCK.LT.1)GOTO30
         IF(ITER.GE.ITM)GOTO30
         ITER=ITER+1
         GOTO45
  30     CONTINUE
         CALL SEM(Y,NY,MS,NBS,NFS,NUS)
         DO 44 JJ=1,NCL
            J=NDA(JJ)
  44        Y(J)=AMAX1(Y(J),YMIN)
  78     CONTINUE
 777  RETURN
      END
C
C----------------------------------------------------------
      SUBROUTINE SEM(Y,NC,NS,NB,NF,NU)
C
C     REDISTRIBUTION OF FAMILY GASES
C
      DIMENSION Y(NC),NB(NS),NF(NC,NS),NU(NC,NS)  !it was Y(*) GM3/01
C  NS is number of families, NB is number of species belonging to
C  each family, NF indicates which species belongs to which family, and
C  NU takes into account how many molecules of the total family value a
C  specific species makes up (e.g., N2O5=2, NO=1). GM3/01
C
      DO 1 IVA=1,10
         DO 1 J=1,NS
            YY=0.
            DO 2 JI=1,NB(J)
   2           YY=YY+Y(NF(JI,J))*NU(JI,J)
            IF(YY.NE.0.)THEN
               YYY=Y(J)/YY
            ELSE
               YYY=1.
            ENDIF
            DO 4 JI=1,NB(J)
C   4           Y(NF(JI,J))=(0.55*YYY+0.45)*Y(NF(JI,J))
   4           Y(NF(JI,J))=YYY*Y(NF(JI,J))         ! SPS /12/15/00
   1        CONTINUE
      RETURN
      END
