subroutine solarband (clm,den,ntau,dj0,xu,wu,nu,opp,dopp,wpp,npp
     & ,q,dq)

c *******************************************************************
c *                                                                 *
c *    solarband calculates the direct heating rate and the cascade *
c *    rate caused by absorption of solar photons in overtone and   *
c *    combination bands.                                           *
c *                                                                 *
c *******************************************************************

      parameter (zero=0.d0,half=5.d-1,one=1.d0,two=2.d0)

c     ---------------------------------------------------------------
c                    Input variables and arrays
c     ---------------------------------------------------------------
      real*8 clm(1),den(1),dj0,xu(1),wu(1),opp(1000,1),dopp(1000,1)
     & ,wpp(1)
      integer*4 ntau,nu,npp

c     ---------------------------------------------------------------
c                    Output variables and arrays
c     ---------------------------------------------------------------
      real*8 q(1), dq(65,1)

c     ---------------------------------------------------------------
c                     Internal variables and arrays
c     ---------------------------------------------------------------
      real*8 wjf,gf(65),dgf(65),tau(65),sum1,sum2,dj(65),dj2(65)
      integer*4 jt,jf,ju,j1,j2

c     ---------------------------------------------------------------
c                  Initialize variables and arrays
c     ---------------------------------------------------------------
      do j1 = 1, ntau
         do j2 = 1, ntau
            dq(j2,j1) = zero
         end do
         q(j1) = zero
      end do

c     ---------------------------------------------------------------
c                    Loop through frequencies 
c     ---------------------------------------------------------------
      do jf = 1, npp

         wjf =  wpp(jf)
         do jt = 1, ntau
            gf(jt) = opp(jf,jt)
            dgf(jt) = dopp(jf,jt)
         end do

c        ------------------------------------------------------------
c            Calculate overhead monochromatic optical depth
c        ------------------------------------------------------------
         tau(1) = gf(1)*clm(1)
         do jt = 2, ntau 
            tj = tau(jt-1)+half*(gf(jt)+gf(jt-1))
     &       *(clm(jt)-clm(jt-1))
            tau(jt) = dmin1(tj,33.d0)
         end do

c        ------------------------------------------------------------
c                Calculate diurnally averaged insolation
c        ------------------------------------------------------------
         do jt = 1, ntau
            sum1 = zero
            sum2 = zero
            do ju = 1, nu
               sum1 = sum1 + wu(ju)*dexp(-tau(jt)/xu(ju))
               sum2 = sum2 + wu(ju)*dexp(-tau(jt)/xu(ju))/xu(ju)
            end do
            dj(jt) = half*sum1*dj0
            dj2(jt) = half*sum2*dj0
         end do

c        ------------------------------------------------------------
c                Calculate solar energy absorption rate, q 
c        ------------------------------------------------------------
         do jt = 1, ntau
            q(jt) = q(jt) + wjf*gf(jt)*dj(jt)
         end do

c        ------------------------------------------------------------
c                Calculate temperature derivative of q, dq
c        ------------------------------------------------------------
         do jt = 1, ntau
            dq(jt,jt) = wjf*dgf(jt)*dj(jt)
         end do

         j1 = 1
         dq(j1,j1)=dq(j1,j1)+wjf*dgf(j1)*clm(j1)*dj2(j1)

         do j1 = 2, ntau 

            j2 = 1
            dq(j1,j2)=dq(j1,j2)+wjf*dgf(j2)*half
     &       *(clm(j2+1)-clm(j2))*dj2(j1)
            
            do j2 = 2, j1 - 1
               dq(j1,j2)=dq(j1,j2)+wjf*dgf(j2)*half
     &          *(clm(j2+1)-clm(j2-1))*dj2(j1)
            end do
              
            j2 = j1
            dq(j1,j2)=dq(j1,j2)+wjf*dgf(j2)*half
     &       *(clm(j2)-clm(j2-1))*dj2(j1)

         end do  

      end do

      do j1 = 1, ntau
         q(j1) = den(j1)*q(j1)
         do j2 = 1, ntau
            dq(j1,j2) = den(j1)*dq(j1,j2)
         end do
      end do 

      return
      end