      subroutine rflsns(ichan,iofset,ilat,xlat,grref,clref,
     1 sza,pcloud,pteran,estozn,fteran,lprint,clfrac,ref,
     3 iloprf,ihiprf,dndr)
c
c***********************************************************************
c
c     rflsns
c
c     august, 2001 by charlie wellemeyer of ssai
c
c     purpose
c       compute reflectivty sensitivity given cloud fraction and refl 
c
c     method
c       establish assumed ground and cloud reflectivities.
c       select interpolation indices for table lookup.
c       perform table interpolations in pressure and ozone,
c       then calculate sensitivity for given refl or cloud fraction.
c
c     variables
c       name      type  i/o   description
c       ----      ----  ---   -----------
c      arguments
c       ichan     i*4    i    channel index for calculation
c       iofset    i*4    i    offset pointer.
c       ilat      i*4    i    latitude index (1 - low 2 - mid 3 - high)
c       xlat      r*8    i    latitude
c       grref     r*8    i    ground reflectivity
c       clref     r*8    i    cloud reflectivity
c       sza       r*8    i    solar zenith angle for channel ichan
c       pcloud    r*8    i    cloud top pressure
c       pteran    r*8    i    terrain pressure
c       estozn    r*8    i    current total ozone estimate
c       fteran    r*8    i    interpolation fraction for terrain pressure
c       lprint    l*4    i    if = .true., print
c       clfrac    r*8    o    cloud fraction
c       ref       r*8    o    reflectivity for wavelength ichan
c                             profiles (needed for clfrac < 0, > 1)
c       iloprf    i*4    i    profile index at low bound
c       ihiprf    i*4    i    profile index at high bound
c       dndr      r*8    o    refl sensitivity for wavelength ichan
c
c***********************************************************************
c
c     -- input parameters
      implicit none
      real xlat, pcloud, pteran, estozn, sza, ref, dndr
      real grref, clref, xnvalm(12)
      real ramcor(8,6)
      integer iofset(12), isnow, iloprf, ihiprf, ilat, ichan
      logical frstrf, lprint(30)
c     -- internal parameters
      real ezero, t, qgclc(4), qcclc(4), qsavlo, qsavhi
      real den, fteran, oznind, dlref, cnvrt, omeglo, omeghi,
     1     ozfrac, ezlggr, tloggr, sbgrnd, grad, ezlgcl, tlogcl,
     2     sbcld, crad, rad1, rad2, grad1, grad2, crad1, crad2
      real ezgr(2),ssgr(2),tgr(2),sbgr(2),ezcl(2),sscl(2),
     1  tcl(2),sbcl(2),radgr(2),radcl(2),qgcalc(2),qccalc(2),sb
      data cnvrt/0.017453/
      integer iprfl, ioz
c     -- output parameters
      real clfrac, pgrref, pclref, pcfrac, dnd331
c     -- common area and data statements
      include 'stndprof.com'
c
      if (lprint(11)) then
         write (6,1000) xlat, iofset(ichan), grref, clref, 
     1    ichan, pcloud, pteran, ilat, sza,
     2    estozn, fteran
      endif
c
c     -- use profile indeces determined by calling routine
c
      iprfl = iloprf - 2
c
c     -- do two calculations for bracketting ozone profiles
c
      do ioz=1,2
        iprfl = iprfl + 1
c
c     -- compute i0, tr, and sb values
c
        call iztrsb(ichan,xlat,iprfl,pteran,pcloud,
     1    iofset(ichan),lprint,ezgr(ioz),ssgr(ioz),tgr(ioz),
     2    sbgr(ioz),ezcl(ioz),sscl(ioz),tcl(ioz),sbcl(ioz))
c
      enddo
c
c     -- compute delta ozone including terrain height corrections
c
      omeglo = profl(iloprf) - terroz(iloprf)*fteran
      omeghi = profl(ihiprf) - terroz(ihiprf)*fteran
      ozfrac = (estozn-omeglo)/(omeghi-omeglo)
c
c     -- Interpolate log of table parameters to current ozone estimate
c     -- and convert to radiance for cloud fraction calculation
c
      ezlggr=alog10(ezgr(1))+alog10(ezgr(2)/ezgr(1))*ozfrac
      sbgrnd=sbgr(1)+(sbgr(2)-sbgr(1))*ozfrac
      if(tgr(2)/tgr(1).gt.0.0) then
        tloggr=alog10(tgr(1))+alog10(tgr(2)/tgr(1))*ozfrac
        grad=10.0**ezlggr+grref*10.0**tloggr/(1.0-grref*sbgrnd)
      else
        tloggr=-99999.0
        grad=10.0**ezlggr
      endif
c
      ezlgcl=alog10(ezcl(1))+alog10(ezcl(2)/ezcl(1))*ozfrac
      sbcld=sbcl(1)+(sbcl(2)-sbcl(1))*ozfrac
      if(tcl(2)/tcl(1).gt.0.0) then
        tlogcl=alog10(tcl(1))+alog10(tcl(2)/tcl(1))*ozfrac
        crad=10.0**ezlgcl+clref*10.0**tlogcl/(1.0-clref*sbcld)
      else
        tlogcl=-99999.0
        crad=10.0**ezlgcl
      endif
c
c     -- calculate reflectivity using version 6 method for cloud free
c
      if (clfrac.le.0.0d0) then
        rad1=10.0**ezlggr+grref*10.0**tloggr/(1.0-grref*sbgrnd)
        dlref=grref+0.02
        rad2=10.0**ezlggr+dlref*10.0**tloggr/(1.0-dlref*sbgrnd)
        dndr=-log10(rad2/rad1)/0.02
c
c     -- calculate reflectivity using version 6 method for cloudy case
c
      else if (clfrac.ge.1.0d0) then
        rad1=10.0**ezlgcl+clref*10.0**tlogcl/(1.0-clref*sbcld)
        dlref=clref+0.02
        rad2=10.0**ezlgcl+dlref*10.0**tlogcl/(1.0-dlref*sbcld)
        dndr=-log10(rad2/rad1)/0.02
c
c     -- calculate reflectivity using cloud fraction and nominal refls
c
      else
        dlref=grref+0.02
        grad1=10.0**ezlggr+grref*10.0**tloggr/(1.0-grref*sbgrnd)
        grad2=10.0**ezlggr+dlref*10.0**tloggr/(1.0-dlref*sbgrnd)
        dlref=clref+0.02
        crad1=10.0**ezlgcl+clref*10.0**tlogcl/(1.0-clref*sbcld)
        crad2=10.0**ezlgcl+dlref*10.0**tlogcl/(1.0-dlref*sbcld)
        rad1 = grad1 + clfrac*(crad1-grad1)
        rad2 = grad2 + clfrac*(crad2-grad2)
        dndr=-log10(rad2/rad1)/0.02
c
      endif
c
      if (lprint(11)) write (6,1200) dndr
c
      return
c
c     -- format statements
c
1000  format (/'Subroutine rflsns'/'Input:   ','xlat    = ',f8.3,
     1 ' iofset = ',i8/
     2 '         grref  = ',f8.4,' clref  = ',f8.4,' ichan  = ',i8/
     3 '         pcloud = ',f8.4,' pteran = ',f8.4,/,
     4 '         ilat =   ',i3,'      sza = ',f6.2,/,
     5 '         estozn = ',f7.1,'  fteran = ',f6.3)
1200  format ('Output:  rsns  = ',f8.4)
c
      end
