	Subroutine rdstab(lun4,lun3,nlam,nsza,no3,rlam,sza_tab,
     1		o3_tab,sdat,o2jdat)

C NAME: rdstab
C PURPOSE: read in radiation source function data table for subsequent
C		interpolation.  Also read table indexing arrays.
C CATEGORY:
C CALLING SEQUENCE: Call rdstab(lun4,lun3,nlam,nsza,no3,rlam,sza_tab,
C					o3_tab,sdat)
C INPUTS:
C   lun4 -- unit number of file containing table index values, ie., pressures,
C		wavelengths, sza, and o3 cols (no3,npr) used to generate stab
C		Values were written at same time as stab.
C   lun3 -- unit number of file containing s data table.  
C		Data must be on wavelength grid consistent with x-sections.
C		Table is written by wrt_stab.pro
C   nlam -- number of wavelength intervals used (checked against read-in)
C   nsza -- number of zenith angles in table
C   no3s -- number of column o3 values in table
C OPTIONAL INPUT PARAMETERS:
C KEYWORD PARAMETERS:
C OUTPUTS:
C   rlam -- array of wavelength values (nlam) used to generate s table
C   sza_tab -- array of sza values (nsza) used to generate s table
C   o3_tab -- array of col o3 values (no3,npr) used to generate s table
C   sdat -- array of s data (nzens,no3s,nprs,nwavelengths)
C   o2jdat -- array (nzens,no3s,nprs) of J(O2) values from srcfun output
C COMMON BLOCKS:
C SIDE EFFECTS:
C PROCEDURE:
C   Straightforward.
C RESTRICTIONS:
C   
C REQUIRED ROUTINES: Open (FORM='system') and Close units in calling
C			 program
C MODIFICATION HISTORY: 
C   Created 930824 - SR Kawa
C   Modified 960705 to read in J(O2) values from srcfun appended to s table,
C	and use 28-level table
C   $Header: /home/kawa/3d/rdstab.f $
C-

C40

	PARAMETER (npr=28)		!number of model press levels
	DIMENSION sdat(nsza,no3,npr,nlam)
	DIMENSION pr_tab(npr),rlam(nlam),sza_tab(nsza),o3_tab(no3,npr)
	DIMENSION o2jdat(nsza,no3,npr)

C Read in dimensions of index arrays and check for consistency
	READ (lun4,*) npr_in,nlam_in,nsza_in,no3_in
	PRINT *,'s data array dimensions:',nsza_in,no3_in,npr_in,nlam_in
	write (6,*) 
        IF((npr_in.NE.npr).OR.(nlam_in.NE.nlam).OR.(nsza_in.NE.nsza).OR.
     1		(no3_in.NE.no3)) THEN 
	   PRINT *,'Array dimensions of table indices do not match',
     1		' those expected',npr_in,npr,nlam_in,nlam,nsza_in,
     2		nsza,no3_in,no3
	   RETURN
	ENDIF

C  Read in indexing values
	READ (lun4,*) pr_tab
	READ (lun4,*) rlam
	READ (lun4,*) sza_tab
	READ (lun4,*) o3_tab

C Read in s function data in binary format.  
	READ (lun3) sdat
	READ (lun3) o2jdat

	RETURN
	END


       Subroutine jcalctraj(ilon,sza,colo3,t,xtab,sza_tab,
     1		o3_tab,sdat,aj,ipmax,ipmin,o2jdat)

C+
C NAME: jcalctraj
C PURPOSE:
C   calculate jrates for chemistry along trajectory model
C CATEGORY:
C CALLING SEQUENCE:
C   Call jcalc(ilon,sza,colo3,t,xtab,sza_tab,o3_tab,sdat,aj)
C INPUTS:
C   ilon -- longitude index for points at which Js are desired
C   sza -- array of zenith angle values (nlat)
C   colo3 -- array of overhead O3 values (nlat,nz)
C   t -- array of CTM temperature values
C   xtab -- table (nlam,nts,nx) of x-section data for range of Ts (150-349 K)
C   sza_tab -- array of sza values (nsza) used to generate s table
C   o3_tab -- array of col o3 values (no3,npr) used to generate s table
C   sdat -- array of s data (nsza,no3s,nprs,nwavelengths)
C   ipmax -- maximum pressure level (lowest) for which to calc Js
C   ipmin -- minimum pressure level (highest) for which to calc Js
C   o2jdat -- array (nzens,no3s,nprs) of J(O2) values from srcfun output
C OUTPUTS:
C   aj -- modified array of photolysis rates
C RESTRICTIONS:
C    currently set up for 22-J set (see var nx)
C REQUIRED ROUTINES:
C     interp_s, sortj
C MODIFICATION HISTORY: 
C   Created 930826 - SR Kawa
C   931123-SRK-remade xtab to do multiplication by solar flux beforehand and
C     removed inputs:
C     xdat -- cross-section, quantum yield, solar flux data read in once before
C     xnames -- character identifyers for xdat
C   940225 modified to add 3 additional Js, incl N2O
C   950322 modified to run over limited range of pressures for traj chem, SRK
C   951025 modified to include CH2O, CH3OOH, SRK
C   960705 modified to handle J(O2) separately and use 28 levels
C-

	PARAMETER (nlat=1,nlon=1,nz=28)
	PARAMETER (nx=22,nread=25,nlam=79,nsza=20,no3=12)
	PARAMETER (nts=200)

	DIMENSION sza(nlat),colo3(nlat,nz),t(nlon,nlat,nz)	!CTM inputs
	DIMENSION xtab(nlam,nts,nx)	!x-section data
	DIMENSION sdat(nsza,no3,nz,nlam)		!flux data
	DIMENSION o2jdat(nsza,no3,nz)			!jo2 data
	DIMENSION sza_tab(nsza),o3_tab(no3,nz)	!flux indexing data

	DIMENSION aj(37,nlat,nz)			!ctm j rates

	REAL s(nlat,nz,nlam),rjs(nlat,nz,nx)
	REAL jo2(nlat,nz)				!jo2 update
        real temp (nlat,nz,nx,nlam)
        integer ind(nlat,nz) 

        REAL a220(6),a300(6),prset(nz)		!for m-dep of J(CH2O)
        Data a220 /0.00000,5.38571e-20,1.28143e-19,2.02429e-19,
     1	    2.76714e-19,3.51000e-19/
 	Data a300 /1.42857e-21,3.21429e-20,6.28571e-20,9.35714e-20,
     1	    1.24286e-19,1.55000e-19/
 	Data prset /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 /
  
	DO iz=1,nz
	DO ilat=1,nlat
           ind(ilat,iz) = MIN(MAX(t(ilon,ilat,iz)-148.5,1.),200.) 
        end do
        end do       

c       interpolate radiative flux function values to model conds. 
        call interp_s_traj(ilon,nlam,sza,colo3,sdat,sza_tab,o3_tab,s,
     1	    ipmax,ipmin,o2jdat,jo2) 

c       calculate J values, xdat(nread) = solar flux

        do 31 iz=ipmax,ipmin
        do 32 ilat=1,nlat
          indt = ind(ilat,iz) 
          do 33 ix=1,nx 
           do 34 ilam=1,nlam
             temp(ilat,iz,ix,ilam) = s(ilat,iz,ilam)*
     1            xtab(ilam,indt,ix)
 34   	   continue
 33       continue
          DO ilam=55,59	   !CH2O quantum yield dependence on m, T, wavel
            ia=ilam-54
       	    atemp=AMIN1(a220(ia),AMAX1(0.,(t(ilon,ilat,iz)-220.)/80.
     1		*(a300(ia)-a220(ia))+a220(ia)))			!T-dep
            alpha=atemp*7.243e18*prset(iz)/t(ilon,ilat,iz)	!m-dep
            qy=1./(1.+alpha)			!total quantum yield at p & T
            qy300=1./(1.+a300(ia)*2.45e19)	!qy at ats press and room T
            rq1q2=xtab(ilam,indt,21)/xtab(ilam,indt,22)
            sigf=(xtab(ilam,indt,21)+xtab(ilam,indt,22))/qy300
            temp(ilat,iz,22,ilam) = s(ilat,iz,ilam)*sigf*qy/(rq1q2+1.)
            temp(ilat,iz,21,ilam) = temp(ilat,iz,22,ilam)*(rq1q2)
          ENDDO
 32   continue
 31   continue

        ilam=1
        do ix=1,nx 
          do iz=ipmax,ipmin
            do ilat=1,nlat
              rjs(ilat,iz,ix) = temp(ilat,iz,ix,ilam)
            end do 	
          end do
        end do
         
      do 51 ilam=2,nlam
         do 52 ix=1,nx 
          do 53 iz=ipmax,ipmin
            do 54 ilat=1,nlat
               rjs(ilat,iz,ix)=rjs(ilat,iz,ix)+temp(ilat,iz,ix,ilam)

CCCCC             rjs(ilat,iz,ix)=rjs(ilat,iz,ix)+s(ilat,iz,ilam)*
CCC  1            xtab(ilam,ind(ilat,iz),ix)*xdat(ilam,nread)
CCC  1            xtab(ilam,indt,ix)*xdat(ilam,nread)
 54   continue
 53   continue
 52   continue
 51   continue

C    Replace J(O2) values
	DO iz=1,nz					!jo2 update
	    DO ilat=1,nlat
	    	rjs(ilat,iz,15)=jo2(ilat,iz)
	    ENDDO
	ENDDO
        
CCCCC   write (*,*) '**CHECK VALUES** ',aj(17), aj(5), aj(2)
c       Order photolysis rates to match order in full chemistry model.  
        call sortj(rjs,aj,nx) 
       

        return
        end 



      Subroutine interp_s_traj(ilon,nlam,sza,ovho3,stab,sza_tab,o3_tab,
     1	    s,ipmax,ipmin,o2jdat,jo2)

C+
C NAME:
C   interp_s_traj
C PURPOSE:
C   Interpolate s values for each wavelength in table to specified O3
c   col and zenith angles
C CATEGORY:
C CALLING SEQUENCE:
C   Call interp_s(nlam,sza,ovho3,stab,sza_tab,o3_tab,s)
C INPUTS:
c   ilon -- longitude index 
C   nlam -- number of wavelength intervals used
C   sza --  input array (nlat) zenith angle values for which s
c           is desired
C   ovho3 --input array (nlat,nz) overhead o3 column value 
C   stab -- input array of values of radiative source function 
C           (nzens,no3,np_ctm,nlam) gridded to ctm p layers
C   sza_tab -- values of sza corresponding to stab table values
C   o3_tab -- array of overhead O3 values at each p-level (no3s,np_ctm)
C                used to index stab
C   ipmax -- maximum pressure level (lowest) for which to calc Js
C   ipmin -- minimum pressure level (highest) for which to calc Js
C   o2jdat -- array (nzens,no3s,nprs) of J(O2) values from srcfun output
C OPTIONAL INPUT PARAMETERS:
C OUTPUTS:
C   s -- array of s values (nlat,nz,nlam) for each wavelength 
C             at specified model p-levels interpolated to ovho3 and sza values
C   jo2 -- array (nlat,nz) of J(O2) values interped as above (REAL)
C COMMON BLOCKS:
C SIDE EFFECTS:
C PROCEDURE:
C   bi-linear interpolation, for sza>94 s=0, for o3 out of range use min/max
C RESTRICTIONS:
C REQUIRED ROUTINES:
C	stab and other tabs read in previously with rdstab.f
C MODIFICATION HISTORY: 
C   Created 930825 - SR Kawa
C   930825 - modified to operate over reduced range of pressures for traj chem
C   $Header$
C-
C
    	PARAMETER (nsza=20,no3=12,npr=28)
        parameter (nz=28,nlat=1) 

	real stab(nsza,no3,npr,nlam),sza_tab(nsza),o3_tab(no3,npr)
  	real o2jdat(nsza,no3,npr)			!jo2 update
	real s(nlat,nz,nlam),sza(nlat), ovho3(nlat,nz)  
	real jo2(nlat,nz)				!jo2 update

        integer locsza(nlat)
        real t(nlat)

c       For each input solar zenith angle, find the first element of
c       tabled sza_tab values that is greater than it.  Use this
c       table element and previous table element to determined 
c       interpolated value. 
        do ij=1,nlat
           do is=1,nsza
              if (sza_tab(is).gt.sza(ij)) go to 333 
           end do 
 333       locsza(ij) = is 
        end do
            
        do ij=1,nlat
           if (locsza(ij).eq.(nsza+1)) then 
c             Point is in darkness.  Set s's to zero. 
              do ik=1,nz 
              do il=1,nlam
                 s(ij,ik,il) = 0. 
                 jo2(ij,ik) = 0.
              end do 
              end do 
           else  
              ijj = locsza(ij)                       
              t(ij)=(sza(ij)-sza_tab(ijj-1))/
     1         (sza_tab(ijj)-sza_tab(ijj-1)) 
     
c             For each input overhead o3 column find the first element
c              of tabled o3_tab values that is > than it.  Use this
c              table element and previous table element to determine
c              interpolated value.             
Crk              do ik=1,nz 
              do ik=ipmax,ipmin 
                 do is=1,no3 
                    if (o3_tab(is,ik).gt.ovho3(ij,ik)) go to 334 
                 end do 
 334             ikk = is 
                 ikkm = ikk-1 
                 if ((ikk.gt.1).and.(ikk.le.no3)) then 
                    u=(ovho3(ij,ik)-o3_tab(ikkm,ik))/
     1               (o3_tab(ikk,ik)-o3_tab(ikkm,ik))          

c                   do bilinear interpolation at ik for each wavelength
c                   from numerical recipes, p.96

	            do il=1,nlam       
	               s(ij,ik,il)=(1.-t(ij))*(1.-u)*
     1                  stab(ijj-1,ikkm,ik,il)+t(ij)*(1.-u)*
     1		        stab(ijj,ikkm,ik,il)+t(ij)*u*
     1                  stab(ijj,ikk,ik,il)+
     1	         	(1.-t(ij))*u*stab(ijj-1,ikk,ik,il)
                    end do

	            jo2(ij,ik)=(1.-t(ij))*(1.-u)*o2jdat(ijj-1,ikkm,ik)
     1			+t(ij)*(1.-u)*o2jdat(ijj,ikkm,ik)
     2			+t(ij)*u*o2jdat(ijj,ikk,ik)
     3	         	+(1.-t(ij))*u*o2jdat(ijj-1,ikk,ik)
                
                 else if (ikk.eq.1) then 
c                    write (6,33) ilon,ij,ik,ovho3(ij,ik),o3_tab(1,ik)
c 33                 format(' Ovhd o3 col(',i4,i3,i3,') of ',e10.3,
c     1               ' too thin! Lowest tabled val=',e10.3)  
      
                    do il=1,nlam
	               s(ij,ik,il)=(1.-t(ij))*stab(ijj-1,1,ik,il)+
     1                  t(ij)*stab(ijj,1,ik,il)
	            end do        
	            jo2(ij,ik)=(1.-t(ij))*o2jdat(ijj-1,1,ik)+
     1                  t(ij)*o2jdat(ijj,1,ik)

                 else 
c                    write (6,34) ilon,ij,ik,ovho3(ij,ik),o3_tab(no3,ik)
c  34                 format(' Ovhd o3 col(',i4,i3,i3,') of ',e10.3,
c     1               ' too thick! Highest tabled val=',e10.3) 
                    do il=1,nlam
	               s(ij,ik,il)=(1.-t(ij))*stab(ijj-1,no3,ik,il)+
     1                  t(ij)*stab(ijj,no3,ik,il)
	            end do 
	            jo2(ij,ik)=(1.-t(ij))*o2jdat(ijj-1,no3,ik)+
     1                  t(ij)*o2jdat(ijj,no3,ik)
                 end if  
              end do 
	    end if
           end do 

           return
           end 


C ------------------------------------------------------------------------

	Subroutine sortj(rjs,aj,nx)

C sort rjs into CTM photolysis rate array, aj

	PARAMETER (nlat=1,nlon=1,nz=28)

	DIMENSION aj(37,nlat,nz)			!ctm j rates
	DIMENSION rjs(nlat,nz,nx)

C  Order of rjs:
C  1-J(BrONO2)
C  2-J(BrO)
C  3-J(Cl2O2)
C  4-J(ClONO2)
C  5-J(H2O2)
C  6-J(HCl)
C  7-J(HNO3)
C  8-J(HO2NO2)
C  9-J(HOCl)
C  10-J(N2O5)
C  11-J(NO2)
C  12-J(NO3_NO)
C  13-J(NO3_NO2)
C  14-J(OClO)
C  15-J(O2)
C  16-J(O3_O1D)
C  17-J(O3_3P)
C  18-J(HOBr)
C  19-J(CH3OOH)
C  20-J(N2O)
C  21-J(CH2O_HCO)
C  22-J(CH2O_CO)
        

        do iz=1,nz
        do ilat=1,nlat 
           aj(1,ilat,iz)=rjs(ilat,iz,15)
	   aj(2,ilat,iz)=rjs(ilat,iz,16)
	   aj(3,ilat,iz)=rjs(ilat,iz,17)
           aj(4,ilat,iz) = 0. 
   	   aj(5,ilat,iz)=rjs(ilat,iz,13)
	   aj(6,ilat,iz)=rjs(ilat,iz,7)
	   aj(7,ilat,iz)=rjs(ilat,iz,11)
	   aj(8,ilat,iz)=rjs(ilat,iz,5)
	   aj(9,ilat,iz)=rjs(ilat,iz,10)
           aj(10,ilat,iz) = rjs(ilat,iz,21)
           aj(11,ilat,iz) = rjs(ilat,iz,22) 
           aj(12,ilat,iz) = 0. 
           aj(13,ilat,iz) = rjs(ilat,iz,19)
	   aj(14,ilat,iz) = rjs(ilat,iz,20) 
           aj(15,ilat,iz)=rjs(ilat,iz,4)
           aj(16,ilat,iz) = 0. 
	   aj(17,ilat,iz)=rjs(ilat,iz,12)
	   aj(18,ilat,iz)=rjs(ilat,iz,6)
	   aj(19,ilat,iz) = 0. 
           aj(20,ilat,iz) = 0. 
           aj(21,ilat,iz) = 0. 
           aj(22,ilat,iz) = 0.
           aj(23,ilat,iz)=rjs(ilat,iz,9)
	   aj(24,ilat,iz)=rjs(ilat,iz,8)
	   aj(25,ilat,iz) = rjs(ilat,iz,18) 
           aj(26,ilat,iz) = 0.
           aj(27,ilat,iz)=rjs(ilat,iz,2)
	   aj(28,ilat,iz)=rjs(ilat,iz,1)
	   aj(29,ilat,iz) = 0. 
           aj(30,ilat,iz) = 0. 
           aj(31,ilat,iz) = 0. 
           aj(32,ilat,iz) = 0.
	   aj(33,ilat,iz) = 0. 
           aj(34,ilat,iz) = 0.
	   aj(35,ilat,iz) = 0.
           aj(36,ilat,iz)=rjs(ilat,iz,3)
	   aj(37,ilat,iz)=rjs(ilat,iz,14)
        end do
        end do 

	RETURN
	END
	
	
	Subroutine seto3(o3_tab,xpv,o3fit,o3prof)	
	
C    
C NAME: seto3
C PURPOSE: set O3 profile for photolysis look-up based on pv value 
C		according to input fit parameters
C CATEGORY:
C CALLING SEQUENCE: Call seto3(o3_tab,xpv,o3fit,o3prof)
C					
C INPUTS:
C   o3_tab -- array of col o3 values (no3,npr) used to generate s table,
C	read in calling routine using rdstab.
C   xpv -- scalar value of pv
C   o3fit -- array(8) of coefficients fitting total o3 (DU) to pv such that
C	total O3 = o3fit(1) + o3fit(2)*xpv + o3fit(3)*xpv**2 +... to a 
C	max of seventh order.
C
C OUTPUTS:
C   o3prof -- profile array (nz) of overhead ozone values (molec cm-2)
C
C COMMON BLOCKS:
C SIDE EFFECTS:
C PROCEDURE:
C   Straightforward.  If total O3 derived from pv and fit params is outside
C	table bounds then min/max table values are used with no notification.  
C RESTRICTIONS:   
C REQUIRED ROUTINES: 
C			
C MODIFICATION HISTORY: 
C   Created 950308 - SR Kawa
C   $Header$
C-

	PARAMETER (nz=28,no3=12)	!table dimensions
	REAL o3_tab(no3,nz)	
	REAL o3prof(nz)			!output profile (molec cm-2)
	DIMENSION o3fit(8)		!coefficients to fit total O3 to pv
	
c fit of total column O3 to pv	
	o3col=o3fit(1)						!in DU
	DO iorder=1,7
	    o3col=o3col+o3fit(iorder+1)*xpv**iorder
	ENDDO
	o3col=o3col*(4.464e-8*6.023e23)				!(molec cm-3)
 
	DO io3=2,no3
           IF (o3_tab(io3,1).GT.o3col) THEN
                frac=1.-(o3_tab(io3,1)-o3col)/
     1		    (o3_tab(io3,1)-o3_tab(io3-1,1)) 
                GOTO 334 
           ENDIF
        ENDDO 
        io3=no3				!o3col outside table max, use max prof
        frac=1.
        
 334  	CONTINUE
 	frac=AMAX1(frac,0.)		!o3col outside table min, use min prof
 
 	DO iz=1,nz
 	    diff=o3_tab(io3,iz)-o3_tab(io3-1,iz)
	    o3prof(iz)=o3_tab(io3-1,iz)+frac*diff
	ENDDO

	RETURN
	END

