      subroutine flux(indsol,latind,sza,satza,phi,psi,pteran,ozone,ref,
     1 rad,minref,res,algflg,snow,surf,cluvf,uvf,eryflux,eryclear)
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c  program flux
c
c   calculate downward UV flux for erythemal
c
c  author/date
c
c   colin seftor / Nov 99
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      include 'const.inc'
      include 'minref.inc'
      include 'uvflux.inc'
c     -- input arguments
      integer indsol, snow, surf, algflg, latind(2)
      real*8 sza, satza, phi, psi, pteran, ozone, ref, rad, minref, res
c     -- output arguments
      real*8 uvf(4), eryflux, cluvf(4), eryclear
c     -- internal parameters
      integer ihigh(3),ilow(3),lattab(3), iprofl(2)
      real*8 uvf_int(4), d_flux, df(2), df_clr
      real*8 uvf_pr(2)
      real*8 uvf_oz(2)
      real*8 nv(2) 
      real*8 profl(26), sza_l(4)
      real*8 daybnd(12), eck_ct
      real*8 ww, eflux(22), efclr(22)
      logical success
      real*8 terroz(26)
      real*8 lsza,lsatza,lphi,lpteran,lrs,lrad,tcl,lmref,taucloud,ct
      data ilow/1,7,17/, ihigh /6,16,26/, lattab/-3,5,15/
c     -- data statements
      data profl/225.,275.,325.,375.,425.,475.,125.,175.,225.,275.,325.,
     2 375.,425.,475.,525.,575.,125.,175.,225.,275.,325.,375.,425.,475.,
     3 525.,575./
      data terroz/18.5,18.5,18.4,18.4,18.4,18.4,11.1,11.2,19.7,20.0,
     1 20.0,20.0,19.8,20.2,20.5,21.8,11.1,11.2,12.0,17.4,17.9,18.3,
     2 18.7,19.8,21.2,21.8/
      data daybnd/15.,46.,74.,105.,135.,166.,196.,227.,258.,288.,
     1            319.,349./
c     ------
      if (ozone.lt.700. .and. sza.lt.90) then
        if (snow.gt.0 .or. (snow.eq.0 .and. algflg.ge.10)) then
          if (minref.lt.0.4) then
            if (ref.ge.0.4) then
              minref = .4
            else
              minref = ref
            endif
          else
            if (ref.lt.minref) minref = ref
          endif
        endif
        indref = minref*100. + 1.5
c       -- correct for cloud (using ed's tables and code)
        lsza    = sza
        lphi    = phi
        lsatza  = satza
        lpteran = pteran
        lmref   = minref
        lrad    = rad
        tcl=taucloud(lsza,lphi,lsatza,lpteran,lmref,lrad,success)
c       print *,'lsza,lphi,lsatza,lpteran,lmref,lrad,tcl'
c       print *,lsza,lphi,lsatza,lpteran,lmref,lrad,tcl
c       print *,' '
        do iwav = 1,26
          ww = wav(iwav)
c         print *,'iwav,wav,ww = ',iwav,wav(iwav),ww
c         -- determine cloud correction factor
          call ld_cloud_ccf(ww,tcl,lmref,lpteran,lsza,ct,success)
c         print *,'ww,tcl,lmref,lpteran,lsza,ct'
c         print *,ww,tcl,lmref,lpteran,lsza,ct
          if (ct.gt.1.d0) ct = 1.d0
          if (ct.lt.0.d0) ct = 0.d0
          if (minref.lt.1.) then
            eck_ct = (1.-ref/100.) / (1.-minref)
          else
            eck_ct = 0.
          endif
          if (eck_ct.lt.0.) then
            eck_ct = 0.
          else if (eck_ct.gt.1.) then
            eck_ct = 1.
          endif
          fteran = (1.0 - pteran)/0.6
          do ilat = 1,2
            iloprf = (ozone-25.)/50. + lattab(latind(ilat))
            if (iloprf.ge.ihigh(latind(ilat))) then
              iloprf = ihigh(latind(ilat)) - 1
            else if (iloprf.lt.ilow(latind(ilat))) then
              iloprf = ilow(latind(ilat))
            endif
            tot_oz = ozone + terroz(iloprf)*fteran
            if (tot_oz.gt.profl(iloprf+1) .and. 
     1       iloprf.ne.5 .and. iloprf.ne.9 .and. iloprf.ne.25) then
              iloprf = iloprf + 1
              tot_oz = ozone + terroz(iloprf)*fteran
            endif
            iprofl(1) = iloprf 
            iprofl(2) = iloprf + 1
            do iprf = 1,2
              do ip = 1,2
                do i = 1,4
                  uvf_int(i) = 
     1             uvflux(indsol+i-1,iprofl(iprf),indref,iwav,ip)
c                  print *,'indsol+i-1,iprofl(iprf),indref,iwav,ip'
c                  print *,indsol+i-1,iprofl(iprf),indref,iwav,ip
c    1             ,uvflux(indsol+i-1,iprofl(iprf),indref,iwav,ip)
                  sza_l(i)   = xzlog(indsol+i-1)
                enddo
c               -- use polynomial interpolation on log(sec(sza))
c               -- for angle
                call polint(sza_l,uvf_int,4,dlog(1./dcosd(sza)),
     1           uvf_pr(ip),dy)
                uvf_pr(ip) = exp(uvf_pr(ip))
              enddo
c             -- use linear interpolation on log pressure
              ptlog = dlog(pteran)
              fr_pr = (ptlog-p10log) / (p04log-p10log)
              uvf_oz(iprf) = (1.-fr_pr)*uvf_pr(1) + fr_pr*uvf_pr(2)
            enddo
c           -- calculate "nval"
            nv(1) = dlog(uvf_oz(1))
            nv(2) = dlog(uvf_oz(2))
c           -- use linear interpolation for ozone
            frac    = (tot_oz - profl(iprofl(1))) / 50.
            df(ilat)= exp((1.-frac)*nv(1) + frac*nv(2))
          enddo
c         -- interpolate on latitude profile bands
          flx = (1.-prfrac)*df(1) + prfrac*df(2)
c         -- correct for sun-earth distance
          d_flux = flx / r2
c         -- clear scene 
          df_clr = d_flux
c         -- correct for clouds
          if (ref.gt.minref .and. tcl.gt.0.) d_flux = d_flux * ct
c         -- correct for aerosol 
c         -- (use aerosol index as extinction coefficient)
          if (res.gt.0.5 .and. ref.lt.0.15)
     1     d_flux = df_clr * exp(-0.25*res)
          if (d_flux.lt.0.) d_flux = 0.
c         -- calculate erythemal flux
          if (iwav.le.4) then
            uvf(iwav)   = d_flux
            cluvf(iwav) = df_clr
          else 
            eflux(iwav-4) = d_flux
            efclr(iwav-4) = df_clr
          endif
        enddo
        eryflux = 0.
        eryclear = 0.
        do j = 1,21
          eryflux  = eryflux  + 
     1               5. * (eflux(j)*ery(j) + eflux(j+1)*ery(j+1)) / 2.
          eryclear = eryclear + 
     1               5. * (efclr(j)*ery(j) + efclr(j+1)*ery(j+1)) / 2.
        enddo
      endif
      return
      end
