      subroutine kphoto(iozon,alb,terrain,hght,zenith,cldtrns,ldark,
     &                  temp,pres)
      use chmstry
      use o3colmap
      implicit none
c 
c----CAMx v7Beta6 190902
c 
c     KPHOTO adjusts all photolysis rate constants for height, ozone
c     column, surface albedo, terrain ht, and zenith angle
c                           
c     Copyright 1996 - 2018
c     Ramboll
c           
c     Modifications: 
c        4/2/03        Removed option for UAM-V type cloud adjustment
c       10/14/04       Modifications to handle Mechanism 10 (no photo rates)
c        7/14/10       Added code for in-line TUV cloud adjustment
c        04/02/12      Removed RADM cloud adjustment option, cloud/aerosol
c                      adjustments now always done with in-line TUV; AHO
c                      file is now just ozone column; replaced haze
c                      dimension with terrain height in photo file
c        04/04/12      Introduced 4-dim interpolation of J values over
c                      zenith, height, albedo, terrain
c        04/12/12      Added T and P adjustments to photolysis rates
c 
c     Input arguments: 
c        iozon               ozone column index
c        alb                 surface albedo
c        terrain             terrain height (km)
c        hght                height AGL (km)
c        zenith              solar zenith angle (deg)
c        cldtrns             enery transmission coefficient
c        ldark               darkness flag (T=dark)
c        temp                temperature (K)
c        pres                pressure (mb)
c 
c     Output arguments: 
c        none
c 
c     Routines Called: 
c        none
c 
c     Called by: 
c        CHEMDRIV
c        PIGDRIVE
c
      include "camx.prm"
c
c-----Arguments
c
      integer iozon
      real alb,terrain,hght,zenith,cldtrns
      real temp,pres
      logical ldark
c
c-----Local
c
      integer nphot,irxn,ihght,izen,ialb,itrn
      integer i
      real tmp,w1,w2,w3,w4
      real prkn1111,prkn1112,prkn1121,prkn1122,prkn1211,prkn1212,
     &     prkn1221,prkn1222,prkn2111,prkn2112,prkn2121,prkn2122,
     &     prkn2211,prkn2212,prkn2221,prkn2222,tmp11,tmp12,tmp13,
     &     tmp14,tmp15,tmp16,tmp17,tmp18,tmp21,tmp22,tmp23,tmp24,
     &     tmp31,tmp32
      real tref,pref,adj1,adj2
      real no2tadj(NZEN),o3tadj(NZEN),hchotadj1(NZEN),hchotadj2(NZEN),
     &     hchopadj(NZEN),ch3chopadj(NZEN)
c
      data tref,pref /298.,1013./
c
c-----Photolysis T and P factors (%/K, %/mb) for 10 zenith angles from TUV
c
      data no2tadj
     &     /0.071,0.072,0.072,0.072,0.072,0.073,0.074,0.076,0.080,0.081/
      data o3tadj
     &     /0.265,0.268,0.275,0.288,0.308,0.339,0.384,0.447,0.485,0.454/
      data hchotadj1
     &     /0.024,0.024,0.024,0.025,0.026,0.026,0.028,0.030,0.031,0.032/
      data hchotadj2
     &     /0.033,0.033,0.033,0.033,0.033,0.034,0.034,0.036,0.036,0.038/
      data hchopadj
     &     /-.017,-.017,-.017,-.017,-.018,-.018,-.019,-.020,-.021,-.023/
      data ch3chopadj
     &     /-1.393E-03,-1.396E-03,-1.402E-03,-1.412E-03,-1.427E-03,
     &      -1.446E-03,-1.473E-03,-1.504E-03,-1.528E-03,-1.536E-03/
c
c-----Entry point
c
c-----If nphot is zero, nothing to do
c
      nphot = nphot1 + nphot2
      if (nphot.eq.0) return
c
c-----If dark, zero all the photolysis reaction rate constants
c
      if (ldark) then
        do irxn = 1,nphot1
          rk(idphot1(irxn)) = 0.
        enddo
        do irxn = 1,nphot2
          rk(idphot2(irxn)) = 0.
        enddo
        goto 999
      endif       
c
c-----Locate the point in height axis
c
      ihght = 1
      do i = 1,NHGHT-1
        if (hght.ge.htint(i)) ihght = i
      enddo
      tmp = (hght - htint(ihght))/(htint(ihght+1) - htint(ihght))
      w1 = 1. - tmp
c
      if (tmp.gt.1.) then ! If height exceeds upper limit
        ihght = NHGHT - 1
        w1 = 0.
      endif
c
c-----Locate the point in zenith angle axis
c
      izen = 1
      do i = 1,NZEN-1
        if (zenith.ge.zenint(i)) izen = i
      enddo
      tmp = (zenith - zenint(izen))/(zenint(izen+1) - zenint(izen))
      w2 = 1. - tmp
c
      if (tmp.gt.1.) then ! If zenith exceeds upper limit
        izen = NZEN-1
        w2 = 0.
      endif
c
c-----Locate the point in albedo axis
c
      ialb = 1
      do i = 1,NALB-1
        if (alb.ge.albcl(i)) ialb = i
      enddo
      tmp = (alb - albcl(ialb))/(albcl(ialb+1) - albcl(ialb))
      w3 = 1. - tmp
c
      if (tmp.lt.0.) then ! If albedo exceeds lower limit
        w3 = 1.
      endif
      if (tmp.gt.1.) then ! If albedo exceeds upper limit
        ialb = NALB-1
        w3 = 0.
      endif
c
c-----Locate the point in terrain axis
c
      itrn = 1
      do i = 1,NTRN-1
        if (terrain.ge.trncl(i)) itrn = i
      enddo
      tmp = (terrain - trncl(itrn))/(trncl(itrn+1) - trncl(itrn))
      w4 = 1. - tmp
c
      if (tmp.lt.0.) then ! If terrain exceeds lower limit
        w4 = 1.
      endif
      if (tmp.gt.1.) then ! If terrain exceeds upper limit
        itrn = NTRN-1
        w4 = 0.
      endif
c
c-----Set the primary photolysis rates using interpolation
c     in height, zenith angle, albedo, and terrain elevation
c
      do irxn = 1,nphot1
        prkn1111 = prkn(izen,  irxn,ihght,  itrn,  ialb,  iozon)
        prkn1112 = prkn(izen,  irxn,ihght+1,itrn,  ialb,  iozon)
        prkn1121 = prkn(izen+1,irxn,ihght,  itrn,  ialb,  iozon)
        prkn1122 = prkn(izen+1,irxn,ihght+1,itrn,  ialb,  iozon)
        prkn1211 = prkn(izen,  irxn,ihght,  itrn,  ialb+1,iozon)
        prkn1212 = prkn(izen,  irxn,ihght+1,itrn,  ialb+1,iozon)
        prkn1221 = prkn(izen+1,irxn,ihght,  itrn,  ialb+1,iozon)
        prkn1222 = prkn(izen+1,irxn,ihght+1,itrn,  ialb+1,iozon)
        prkn2111 = prkn(izen,  irxn,ihght,  itrn+1,ialb,  iozon)
        prkn2112 = prkn(izen,  irxn,ihght+1,itrn+1,ialb,  iozon)
        prkn2121 = prkn(izen+1,irxn,ihght,  itrn+1,ialb,  iozon)
        prkn2122 = prkn(izen+1,irxn,ihght+1,itrn+1,ialb,  iozon)
        prkn2211 = prkn(izen,  irxn,ihght,  itrn+1,ialb+1,iozon)
        prkn2212 = prkn(izen,  irxn,ihght+1,itrn+1,ialb+1,iozon)
        prkn2221 = prkn(izen+1,irxn,ihght,  itrn+1,ialb+1,iozon)
        prkn2222 = prkn(izen+1,irxn,ihght+1,itrn+1,ialb+1,iozon)
        tmp11 = w1*prkn1111 + (1.-w1)*prkn1112
        tmp12 = w1*prkn1121 + (1.-w1)*prkn1122
        tmp13 = w1*prkn1211 + (1.-w1)*prkn1212
        tmp14 = w1*prkn1221 + (1.-w1)*prkn1222
        tmp15 = w1*prkn2111 + (1.-w1)*prkn2112
        tmp16 = w1*prkn2121 + (1.-w1)*prkn2122
        tmp17 = w1*prkn2211 + (1.-w1)*prkn2212
        tmp18 = w1*prkn2221 + (1.-w1)*prkn2222
        tmp21 = w2*tmp11 + (1.-w2)*tmp12
        tmp22 = w2*tmp13 + (1.-w2)*tmp14
        tmp23 = w2*tmp15 + (1.-w2)*tmp16
        tmp24 = w2*tmp17 + (1.-w2)*tmp18
        tmp31 = w3*tmp21 + (1.-w3)*tmp22
        tmp32 = w3*tmp23 + (1.-w3)*tmp24
        rk(idphot1(irxn)) = w4*tmp31 + (1.-w4)*tmp32
      enddo
c
c-----Cloud/aerosol adjustment
c
      do irxn = 1,nphot1
        rk(idphot1(irxn)) = rk(idphot1(irxn))*cldtrns
      enddo
c
c-----Set secondary photolysis rates as ratios to primary rates
c
      do irxn = 1,nphot2
        rk(idphot2(irxn)) = rk(idphot3(irxn))*phtscl(irxn)
      enddo
c
c-----Apply T and P adjustments for current conditions
c
      adj1 = 1. + no2tadj(izen)*(temp - tref)/100.
      adj2 = 1. + no2tadj(izen+1)*(temp - tref)/100.
      rk(jno2rxn) = rk(jno2rxn)*(w2*adj1 + (1.-w2)*adj2)

      adj1 = 1. + o3tadj(izen)*(temp - tref)/100.
      adj2 = 1. + o3tadj(izen+1)*(temp - tref)/100.
      rk(jo3rxn) = rk(jo3rxn)*(w2*adj1 + (1.-w2)*adj2)

      adj1 = 1. + hchotadj1(izen)*(temp - tref)/100.
      adj2 = 1. + hchotadj1(izen+1)*(temp - tref)/100.
      rk(jhcho1rxn) = rk(jhcho1rxn)*(w2*adj1 + (1.-w2)*adj2)

      adj1 = 1. + hchotadj2(izen)*(temp - tref)/100.
     &           *hchopadj(izen)*(pres*tref/temp  - pref)/100.
      adj2 = 1. + hchotadj2(izen+1)*(temp - tref)/100.
     &           *hchopadj(izen+1) *(pres*tref/temp  - pref)/100.
      rk(jhcho2rxn) = rk(jhcho2rxn)*(w2*adj1 + (1.-w2)*adj2)

      adj1 = 1. + ch3chopadj(izen)*(pres*tref/temp  - pref)/100.
      adj2 = 1. + ch3chopadj(izen+1)*(pres*tref/temp  - pref)/100.
      rk(jch3chorxn) = rk(jch3chorxn)*(w2*adj1 + (1.-w2)*adj2)
c
 999  return
      end
