Program phottst 
		
      parameter (long=1,lat=0)		!fossils	
      parameter (levels=1,nlevphot=28)   
      parameter (nx=22,nlam=79,nsza=20,numo3=12)
      parameter (nts=200)

c     Arrays containing photolysis rates, overlying o3, and
c     solar zenith angle. 
      real  aj(37,0:lat,levels) 
      REAL o3prof(nlevphot)
      REAL aj25(37,0:lat,nlevphot),kel25(long,0:lat,nlevphot)
      real  colo3(0:lat,nlevphot),szan(0:lat)  

c     Define arrays needed for calculation of photolysis rates. 
      REAL o3fit(8)			!totalO3-pv fitting coefficients
      real xtab(nlam,nts,nx)
      real sdat(nsza,numo3,nlevphot,nlam)
      real rlam(nlam),sza_tab(nsza),o3_tab(numo3,nlevphot)
      REAL o2jdat(nsza,numo3,nlevphot)

	REAL plevs(nlevphot)
 	Data plevs /921.95, 771.36, 648.07, 547.72, 447.21,
     1	    346.41, 273.86, 223.60, 187.08, 162.02, 139.64,
     2	    118., 100., 88., 77., 68.1292, 57., 46.4159,
     3	    31.6228, 21.5443, 14.678, 10., 6.81292, 4.64159, 3.16228, 
     4	    2.15443, 1.4678, 0.681292 /

C   sample O3 profile (molec cm-2), #7 of the table, about 375 DU
        REAL colo3_7(nlevphot)
        DATA colo3_7/1.1516e+19,1.1447e+19,1.1355e+19,1.1248e+19,
     1	1.1098e+19,1.0851e+19,1.0529e+19,1.0153e+19,
     2	9.7430e+18,9.3623e+18,8.9262e+18,8.3872e+18,
     3	7.8259e+18,7.3824e+18,6.9183e+18,6.4981e+18,
     4	5.9010e+18,5.2369e+18,4.0343e+18,2.9843e+18,
     5	2.1290e+18,1.4368e+18,9.0456e+17,5.3860e+17,
     6	3.0354e+17,1.6312e+17,8.4567e+16,2.1794e+16/


c     Read in lookup table photolysis rate calc.  

c   Read in cross-section, quantum yield, and solar flux data.    
            open(unit=23,status='old',disp='keep',form='unformatted'
     1        ,   file='xsect_22_79.tab')
            lun5=23 
            read (lun5) xtab 
            CLOSE(23)  
C            open(unit=21,status='old',disp='keep',form='system')
            open(unit=21,file='stab_ctm28.xdr',
     1    	status='old',disp='keep',form='system')
            open(unit=22,file='indx_ctm28.dat',
     1    	status='old',disp='keep',form='formatted')
            lun3=21
            lun4=22 
            call rdstab(lun4,lun3,nlam,nsza,numo3,rlam,sza_tab,
     1	    	o3_tab,sdat,o2jdat)
	    CLOSE(21)
	    CLOSE(22)
	    
C	range of pressure set from input data, limiting speeds calculations
	pmax=500.
	pmin=40.

      OPEN(unit=15,name='phot1dtst.out',form='unformatted')

      DO ipt=1,100
        szan(0)=FLOAT(ipt)	!input SZA here
	temp=220.	!these could vary also as in flight track data
	press=50.

Crk  Figure out over what range of pressures to calc photolysis rates 
	DO ip=1,nlevphot-1
	      IF (pmax.GE.plevs(1)) THEN 
	        ipmax=1
	      ELSEIF (pmax.LE.plevs(ip).AND.pmax.GT.plevs(ip+1)) THEN
	    	ipmax=ip
	      ENDIF
	      IF (pmin.GE.plevs(1)) THEN 
	        ipmin=1
	      ELSEIF (pmin.LE.plevs(ip).AND.pmin.GT.plevs(ip+1)) THEN
	    	ipmin=ip+1
	      ENDIF
	      IF (pmin.LT.plevs(nlevphot)) THEN ipmin=nlevphot
	ENDDO

C  Calculate photolysis rates.  In this version, a partial profile of J rates
C    is calculated and then interpolated to the trajectory point in pressure.
c    SZA comes from traj file  
c  Calculate overlying ozone column from total O3 fit to pv. 
Cbells and whistles         Call seto3(o3_tab,xpv,o3fit,o3prof)
	

C  set profile of temperature and O3 for J profile calc
         do ik=1,nlevphot
         do ilat=0,lat
         do ilon=1,long
	    kel25(ilon,ilat,ik)=temp
	    colo3(ilat,ik)=colo3_7(ik)
         end do
         end do
         end do      
        
c        Calculate photolysis rates. 
         call jcalctraj(1,szan,colo3,kel25,xtab,sza_tab,o3_tab,sdat,
     1	    aj25,ipmax,ipmin,o2jdat)      

C  Interpolate J profiles onto pressure of traj point
         do ilat=0,lat
         do ij=1,37
           do ik=2,nlevphot
             IF (plevs(ik).LT.press) THEN
                 frac=(press-plevs(ik-1))/(plevs(ik)-plevs(ik-1))
		GOTO 334
	     ENDIF
	   ENDDO			!nlevphot
	   ik=nlevphot			!press over top of table
	   frac=1.
  334	   CONTINUE
 	   frac=AMAX1(frac,0.)		!press outside table max use bottom prof
  	   diff=aj25(ij,ilat,ik)-aj25(ij,ilat,ik-1)
  	   aj(ij,ilat,1)=aj25(ij,ilat,ik-1)+diff*frac
  	 ENDDO
  	 ENDDO  
  	 
  	 Write(15) szan
  	 WRITE(15) aj 		!write output js at sza
  	    	    
      ENDDO
      CLOSE(15)
      end