      SUBROUTINE DENMOD(UT,DT)                                          
      SAVE
      PARAMETER(KMX=96,KIY=40,MMX=46)
      include "blank.h"
      include "solars.h"
      COMMON/DENTRF/QNN(KMX),QNMN(KMX),QNC(KMX),QN1D(KMX),DSNCHM(KMX),  
     1DSICHM(KMX),QO3PR(KMX),E630Z(KMX),E63DR(KMX),E63N2D(KMX),         
     2E6300(KMX),QHOXO2(KMX),QHOXNO(KMX),QHOXT(KMX),PARCD(KMX),         
     3SPED(KMX),SHAL(KMX),ALAMXX(KMX),ALAMXY(KMX),XNOF(KMX),QJOUL(KMX)  
      COMMON/TRAVLH/RK(50),ALP(20),RB(20),RKM(60)                       
      COMMON/TRAVSL/AA,BB,CC,DD,EE,FF,GG,CI,HH,CJ,CK,RR                 
      COMMON/RGLWTR/QUENCH,A1D,A6300                                    
      COMMON/QIONDC/AR(54),ISWTC,ELEC                                   
      COMMON/PIONTR/PHOXO2,PHOXNO,PHOXI,PHOXT                           
      COMMON/IONTRFF/XNPI(KMX),XNNI(KMX),XNEEE(KMX)
      DIMENSION BCB(3),ROOT(6),B(3),R(6)                                
      DIMENSION ACOF(5),XMASS(40)                                       
      COMPLEX TOOT(4),RLAM                                              
      DIMENSION Y(4), ATOL(4), RWORK(74), IWORK(24)                     
      DATA XME/5.486E-4/,QCHR/1.602E-19/                                
      DATA XMASS/0.,0.,0.,0.,0.,30.,48.,66.,84.,58.,76.,94.,74.,92.,    
     1 110.,32.,64.,80.,50.,36.,19.,37.,55.,73.,91.,109.,127.,145.,     
     2 16.,32.,48.,64.,46.,62.,60.,76.,62.,60.,44.,34./                 
      ERG=1.602E-12                                                     
      FALP14S=0.2                                                       
      FALP12D=0.8                                                       
      FALP23P=0.15                                                      
      FALP21D=0.85                                                      
      FALP34S=0.5                                                       
      FALP32D=0.5                                                       
      A5577=1.28                                                        
      A1S=1.36                                                          
      Y1S=0.05                                                          
      FOE=0.85                                                          
      FN2D=1.                                                           
      QCHR2=QCHR*QCHR                                                   
      BGAUS=0.4                                                         
      BG=BGAUS*1.E-4                                                    
      CTTS=9.57946E+3*BGAUS                                             
      WE=CTTS/XME                                                       
C   ******************************                                      
      DO 55 K=1,KMX                                                     
      CALL RATECOE(TI(K),TE(K),TN(K))                                   
      XINP(K)=(QINP(K)+RK(22)*XIOP(K)*XN2D(K))                          
     1/((RK(6)+RK(7))*XNO2(K)+RK(8)*XNO(K))                             
      XIHP(K)=RK(12)*XNH(K)*XIOP(K)/(RK(11)*XNO(K))                     
  55  CONTINUE                                                          
  900 FORMAT(1X,*I AM HERE*)                                            
      M2=KMX-MMX+1                                                      
      DO 5 M=1,M2                                                       
      CALL RATECOE(TI(M),TE(M),TN(M))                                   
      XIOP(M)=(QOP(M)+RK(8)*XNO(M)*XINP(M)+RK(11)*XIHP(M)*XNO(M))/      
     1(RK(1)*XNO2(M)+RK(2)*XNN2(M)+RK(9)*XNH2(M)+RK(10)*XNH2O(M)+       
     2RK(12)*XNH(M)+RK(22)*XN2D(M))                                     
  5   CONTINUE                                                          
      CALL OPDIFSV(UT,DT)                                               
      DO 44 K=1,KMX                                                     
      CALL RATECOE(TI(K),TE(K),TN(K))                                   
      AA=QNOP(K)+RK(2)*XIOP(K)*XNN2(K)+RK(7)*XINP(K)*XNO2(K)            
      BB=QO2P(K)+RK(1)*XIOP(K)*XNO2(K)+RK(6)*XINP(K)*XNO2(K)            
C    1+RK(21)*XIN2P(K)*XNO2(K)
      CC=RK(4)*XN4S(K)+RK(5)*XNNO(K)                                    
      DD=QN2P(K)                                                        
      EE=RK(3)*XNO(K)+RK(21)*XNO2(K)                                    
      FF=XIOP(K)                                                        
      GG=XINP(K)                                                        
      A4=ALP(1)*ALP(2)*ALP(3)                                           
      A3=ALP(1)*(ALP(2)*EE+ALP(3)*CC)-ALP(1)*ALP(2)*ALP(3)*(FF+GG)      
      A2=ALP(1)*EE*CC-ALP(1)*(ALP(2)*EE+ALP(3)*CC)*(FF+GG)-ALP(1)*ALP(2)
     1*DD-ALP(1)*ALP(3)*BB-ALP(2)*ALP(3)*AA                             
      A1=-ALP(1)*(EE*CC*(FF+GG)+DD*CC+BB*EE)-ALP(2)*EE*(AA+DD)-ALP(3)*  
     1CC*(AA+BB)                                                        
      A0=-EE*CC*(AA+BB+DD)                                              
      ACOF(1)=A0                                                        
      ACOF(2)=A1                                                        
      ACOF(3)=A2                                                        
      ACOF(4)=A3                                                        
      ACOF(5)=A4                                                        
      CALL QUART(ACOF,SOOT,DELTA,G,H,RL,CDELTA,RLAM)                    
      XNE(K)=SOOT                                                       
      XIOP(K)=FF                                                        
      XIN2P(K)=DD/(EE+ALP(3)*XNE(K))                                    
      XINP(K)=GG                                                        
      XIO2P(K)=BB/(CC+ALP(2)*XNE(K))                                    
      XINOP(K)=(AA+EE*DD/(EE+ALP(3)*XNE(K))+CC*BB/(CC+ALP(2)*XNE(K      
     1)))/(ALP(1)*XNE(K))                                               
      XNPI(K)=XNE(K)
      XNEEE(K)=XNE(K)
      XNNI(K)=1.E-20
  44  CONTINUE                                                          
C******************************************************************     
C  *********************************************                        
C     CALL OZONE(UT)                                                    
      DO 4 K=1,KMX                                                      
      CALL RATECOE(TI(K),TE(K),TN(K))                                   
      CALL RATECOF(TI(K),TE(K),TN(K))                                   
      QNN(K)=(RB(1)*XN4S(K)*XNO2(K)*1.4+RB(2)*XN2D(K)*XNO2(K)*1.84+     
     1RB(3)*XN4S(K)*XNNO(K)*2.68+RB(4)*XN2D(K)*XNO(K)*2.38+RB(5)*XN2D(K)
     2*XNE(K)*2.38+RB(6)*XN2D(K)*XNNO(K)*5.63+RB(10)*XN4S(K)*XNOH(K)*2.1
     3+RB(11)*XNNO(K)*XNO3(K)*2.08+RB(12)*XNHO2(K)*XNNO(K)*0.35+RB(13)* 
     4XNNO2(K)*XNO(K)*1.98)*ERG/RHO(K)                                  
      XM=XNO(K)+XNO2(K)+XNN2(K)                                         
      QNMN(K)=(RKM(29)*XNH(K)*XNO3(K)*3.34+RKM(17)*XNO(K)*XNOH(K)*0.72  
     1+RKM(18)*XNO(K)*XNHO2(K)*2.39+RKM(28)*XNH(K)*XNO2(K)*XM*2.0       
     2+RKM(10)*XNO1D(K)*XNH2O(K)*1.23+RKM(11)*XNO1D(K)*XNH2(K)*1.88+    
     3RKM(19)*XNO(K)*XNH2O2(K)*3.44+RKM(20)*XNO(K)*XNH2(K)*0.08+        
     4RKM(21)*XNOH(K)*XNO3(K)*1.73+RKM(22)*XNOH(K)*XNOH(K)*.73+RKM(23)* 
     43.06*                                                             
     5XNOH(K)*XNHO2(K)+RKM(24)*XNOH(K)+XNH2O2(K)*1.35+RKM(25)*XNOH(K)*  
     6XNH2(K)*0.65+RKM(26)*XNHO2(K)*XNO3(K)*1.23+RKM(27)*XNHO2(K)*      
     7XNHO2(K)*1.71+RKM(30)*XNH(K)*XNHO2(K)*2.41+RKM(30)*XNH(K)*XNHO2(K)
     8*2.41+RKM(31)*XNH(K)*XNHO2(K)*1.61+RKM(32)*XNH(K)*XNHO2(K)*2.34+  
     9RKM(33)*XNH(K)*XNH(K)*XM*4.52+RKM(35)*XNCH4(K)*XNOH(K)*5.38+      
     1RKM(36)*XNO(K)*XNCH4(K)*4.65+RKM(37)*XNCH4(K)*XNO1D(K)            
     2+RKM(40)*XNO(K)*XM*XNCO(K)*5.51+RKM(41)*XNCO(K)*XNOH(K)*1.07      
     3)*ERG/RHO(K)                                                      
      QNN(K)=QNN(K)+QNMN(K)                                             
      QNC(K)=(RK(1)*XIOP(K)*XNO2(K)*1.555+RK(2)*XIOP(K)*XNN2(K)*1.088+  
     1RK(3)*XIN2P(K)*XNO(K)*0.70+RK(4)*XIO2P(K)*XN4S(K)*4.213+RK(5)*    
     2XIO2P(K)*XNNO(K)*2.813+RK(6)*XINP(K)*XNO2(K)*2.486+RK(7)*XINP(K)  
     3*XNO2(K)*6.699+ALP(1)*XINOP(K)*XNE(K)*(FALP12D*0.38+FALP14S*      
     42.75)+ALP(2)*XIO2P(K)*XNE(K)*(FALP23P*6.95+FALP21D*4.98)+         
     5ALP(3)*XIN2P(K)*XNE(K)*(FALP34S*5.82+FALP32D*3.44)+RK(8)*         
     6XINP(K)*XNO(K)*0.93+RK(21)*XIN2P(K)*XNO2(K)*3.52+                 
     7RK(22)*XIOP(K)*XN2D(K)*1.45)*ERG/RHO(K)                           
      QN1D(K)=(FALP21D*ALP(2)*XIO2P(K)*XNE(K)+RB(2)*XN2D(K)*XNO2(K))*   
     1 1.96*ERG*(2.3E-11*XNN2(K)/(9.E-3+2.3E-11*XNN2(K)))/RHO(K)        
      O1DEX(K)=FALP21D*ALP(2)*XIO2P(K)*XNE(K)+RB(2)*XN2D(K)*XNO2(K)     
      DSNCHM(K)=RB(1)*XN4S(K)+RB(2)*XN2D(K)                             
      DSICHM(K)=(RK(4)*XIO2P(K)*XN4S(K)+RK(7)*XINP(K)*XNO2(K)+ALP(2)*   
     1XIO2P(K)*XNE(K))/XNO2(K)                                          
      DSO2(K)=DSO2(K)+DSICHM(K)+DSNCHM(K)                               
      QO3PR(K)=(RKM(12)*XNO(K)*XM*XNO(K)*5.12+RKM(13)*XNO(K)*XNO(K)     
     1*XNO2(K)*1.05+RKM(14)*XNO(K)*XNO2(K)*XNO2(K)*1.05+RKM(15)*        
     2XNO(K)*XNO2(K)*XNN2(K)*1.05+RKM(16)*XNO(K)*XNO3(K)*4.06)*ERG/     
     3RHO(K)                                                            
      QN(K)=QNN(K)+QNC(K)+QN(K)+QN1D(K)+QO3PR(K)+QNO3(K)                
      DEN1=RKM(8)*XNO2(K)+RKM(9)*XNN2(K)+RKM(10)*XNH2O(K)+RKM(11)*      
     1XNH2(K)+RKM(37)*XNCH4(K)+A1D                                      
      SR63(K)=SR63(K)*A6300/DEN1                                        
      E630Z(K)=DO3HAR(K)*XNO3(K)*A6300/DEN1                             
      E63DR(K)=FOE*A6300*ALP(2)*XIO2P(K)*XNE(K)/DEN1                    
      E63N2D(K)=FN2D*A6300*RB(2)*XN2D(K)*XNO2(K)/DEN1                   
      E6300(K)=E63DR(K)+E63N2D(K)+SR63(K)+E630Z(K)                      
      XNT=XNO(K)+XNO2(K)+XNN2(K)                                        
      XMOLP=(XIN2P(K)*28.+XIO2P(K)*32.+XIOP(K)*16.+XINOP(K)*30.         
     1+XINP(K)*14.)/XNE(K)                                              
      BGAUS2=BGAUS*BGAUS                                                
      WPOSI=CTTS/XMOLP                                                  
      WNEGI=0.                                                          
      ALAM=0.                                                           
      VINN=2.6E-09*XNT/SQRT(XMOLP)                                      
      QHOXI(K)=0.                                                       
      QHOXO2(K)=0.                                                      
      QHOXNO(K)=0.                                                      
      QHOXT(K)=0.                                                       
      IF(K.GT.KIY) GO TO 10                                             
      HH=ZPHT(K)                                                        
      TTN=TN(K)                                                         
      AM=XNO(K)+XNO2(K)+XNN2(K)                                         
      AN2=XNN2(K)                                                       
      CO2=XNCO2(K)                                                      
      H2O=XNH2O(K)                                                      
      O2=XNO2(K)                                                        
      O1=XNO(K)                                                         
      O3=XNO3(K)                                                        
      ANO1=XNNO(K)                                                      
      ANO2=XNNO2(K)                                                     
C     O2S=XNO2S(K)                                                      
      O2S=1.E-20                                                        
      HNO3=1.E-20                                                       
      QIO2P=QO2P(K)+QN2P(K)+QOP(K)+QINP(K)                              
      QINOP=QNOP(K)                                                     
      O2COL=TXNO2(K)                                                    
      ALPH1=ALP(1)                                                      
      ALPH3=ALP(2)                                                      
      ISWTC=0                                                           
      CALL IONCOMP(HH,TTN,AM,AN2,CO2,H2O,O2,O3,O1,ANO1,ANO2,O2S,HNO3,   
     1QIO2P,QINOP,O2COL,ALPH1,ALPH3)                                    
      QHOXI(K)=PHOXT*QPRO(K)                                            
      QHOXO2(K)=PHOXO2                                                  
      QHOXNO(K)=PHOXNO                                                  
      QHOXT(K)=PHOXT                                                    
      XNPI(K)=AR(4)                                                     
      XNNI(K)=AR(5)                                                     
      XNEEE(K)=AR(2)                                                    
      XIO2P(K)=AR(16)                                                   
      XINOP(K)=AR(6)                                                    
      XNE(K)=AR(2)                                                      
      ALAM=AR(3)                                                        
      XPOSI=AR(4)                                                       
      XNEGI=AR(5)                                                       
      SUM=0.                                                            
      DO 46 J=6,28                                                      
  46  SUM=SUM+XMASS(J)*AR(J)                                            
      SUN=0.                                                            
      DO 47 J=29,40                                                     
  47  SUN=SUN+XMASS(J)*AR(J)                                            
      XMOLP=SUM/AR(4)                                                   
      XMOLN=SUN/AR(5)                                                   
      WNEGI=CTTS/XMOLN                                                  
      WPOSI=CTTS/XMOLP                                                  
      VINN=2.6E-09*XNT/SQRT(XMOLN)                                      
  10  CONTINUE                                                          
      TW=TE(K)                                                          
      TE12=SQRT(TW)                                                     
      VEN=2.33E-11*XNN2(K)*TW*(1.-1.21E-4*TW)+1.82E-10*XNO2(K)*TE12*(1.+
     13.6E-2*TE12)+8.2E-10*XNO(K)*TE12                                  
C     VINP=2.6E-09*XNT/SQRT(XMOLP)                                      
      COLFAC=1.5                                                        
      VINP=0.73E-9*SQRT(TN(K)/1000.)*XNO(K)*COLFAC+0.69E-9*XNN2(K)+     
     10.67E-9*XNO2(K)                                                   
      VEI=34.+4.18*ALOG10(TW**3/XNE(K))/TW**1.5                         
      CN=1./(BG*QCHR)                                                   
      CKOE=CN*WE/(VEI+VEN)                                              
      CKOP=CN*WPOSI/(VINP+VEI)                                          
      CKON=CN*WNEGI/(VINN+VEI)                                          
      PARCD(K)=XNE(K)*1.E+6*QCHR2*(CKOE+(1.+ALAM)*CKOP+ALAM*CKON)       
      VP=VINP                                                           
      VN=VINN                                                           
      VE=VEN+VEI                                                        
      CKPE=CN*WE*VE/(VE*VE+WE*WE)                                       
      CKPP=CN*WPOSI*VP/(VP*VP+WPOSI*WPOSI)                              
      CKPN=CN*WNEGI*VN/(VN*VN+WNEGI*WNEGI)                              
      SPED(K)=XNE(K)*1.E+6*QCHR2*(CKPE+(1.+ALAM)*CKPP+ALAM*CKPN)        
      CKHE=CN*WE*WE/(VE*VE+WE*WE)                                       
      CKHP=CN*WPOSI*WPOSI/(VP*VP+WPOSI*WPOSI)                           
      CKHN=CN*WNEGI*WNEGI/(VN*VN+WNEGI*WNEGI)                           
      SHAL(K)=XNE(K)*1.E+6*QCHR2*(CKHE-(1.+ALAM)*CKHP-ALAM*CKHN)        
      ALAM1=SPED(K)*BGAUS2/RHO(K)*1.E-11                                
      ALAM2=SHAL(K)*BGAUS2/RHO(K)*1.E-11                                
      ALAMXX(K)=ALAM1                                                   
      ALAMXY(K)=ALAM2                                                   
  202 FORMAT(1X,I5,2X,8E12.4)                                           
      XM=XNO(K)+XNO2(K)+XNN2(K)                                         
      XNOF(K)=XN2D(K)+XN4S(K)+XNNO(K)+XNNO2(K)                          
  4   CONTINUE                                                          
      TLAMXX=0.                                                         
      DO 988 K=2,KMX                                                    
      RZ=RHO(K)                                                         
      RZM=RHO(K-1)                                                      
      RW=0.5*(ZPHT(K)-ZPHT(K-1))*1.E+5                                  
      TLAMXX=TLAMXX+(ALAMXX(K)*RZ+ALAMXX(K-1)*RZM)*RW                   
  988 CONTINUE                                                          
      AX=5.163E+11                                                      
      TLAM=TLAMXX*AX                                                    
      AJOUL=7.5E+10                                                     
      EMVV=SQRT(AJOUL*(BGAUS*1.E-3)**2/TLAM)                            
      EMV=EMVV                                                          
      EFIEL=EMV/(BGAUS*1.E-3)
      DO 82 K=1,KMX                                                     
      QJOUL(K)=ALAMXX(K)*EFIEL*EFIEL
      QN(K)=QN(K)+QJOUL(K)                                              
  82  CONTINUE                                                          
      RETURN                                                            
      END                                                               
CDIR$ NOLIST                                                            
      SUBROUTINE QUART(A,ROOT,DELTA,G,H,RL,CDELTA,RLAM)                 
      SAVE
      COMPLEX CP                                                        
      DIMENSION A(5)                                                    
      DATA E/1.E-1200/                                                  
      A0=A(5)                                                           
      A1=A(4)/4.                                                        
      A2=A(3)/6.                                                        
      A3=A(2)/4.                                                        
      A4=A(1)                                                           
      G=A0**2*A3-3.*A0*A1*A2+2.*A1**3                                   
      H=A0*A2-A1**2                                                     
      RI=A0*A4-4.*A1*A3+3.*A2**2                                        
      RJ=A0*(A2*A4-A3**2)-A1*(A1*A4-A3*A2)+A2*(A1*A3-A2**2)             
      RK=A0**2*RI-3.*H**2                                               
      RL=12.*H**2-A0**2*RI                                              
      CH=-RI/12.                                                        
      CG=RJ/4.                                                          
      DELTA=RI**3-27.*RJ**2                                             
      CDELTA=CG**2+4.*CH**3                                             
      CP=CDELTA+E                                                       
      CP=(.5*(CG+CP**.5)+E)**(1./3.)                                    
      RLAM=-2.*REAL(CP)                                                 
      P=A0*RLAM+A1**2-A0*A2+E                                           
      P=SQRT(P)                                                         
      Q=(2.*RLAM+A2)**2-A0*A4+E                                         
      Q=SQRT(Q)                                                         
      PQ=2.*A1*RLAM+A1*A2-A0*A3+E                                       
      P=SIGN(P,Q*PQ)                                                    
      ROOT=P-A1                                                         
      ROOT=(ROOT+SQRT(ROOT**2-A0*(A2+2.*RLAM-Q)))/A0                    
      RETURN                                                            
      END                                                               
