C *********************************************************************
C *                                                                   *
C *     RANDOM calculate the molecular opacity in a specified         *
C *     frequency range.  Required input consists of an array of      *
C *     line positions, widths, strengths, and quantum numbers.       *
C *     The S-1 random model with a Curtis-Godson approximation       *
C *     is used.  The relevant equations are described by Wallace,    *
C *     Prather, and Belton, Ap. J., 1974, Vol 193, pp 481-493.       *
C *                                                                   *
C *********************************************************************

       subroutine random(uu,prs,tmp,col,ntau,str,wid,eng,nlines,efc
     &  ,rmass,wmid,delnu,a)
       implicit double precision (a-h,o-z)
       include "parameters.dat"

C      -----------------------------------------------------------------
C                               VARIABLES
C      -----------------------------------------------------------------
C      
C       INPUT
C
C            uu   - cosine of solar zenith  angle
C            prs  - pressure grid (atm)
C            tmp  - temperature grid (Kelvins)
C            col  - column density grid (cm-2)
C            str  - list of line strengths (cm-1)
C            wid  - Lorentz widths of lines at a ref. temperature (cm-1)
C            eng  - energy of lower state (cm-1)
C            nlines - number of lines in this frequency interval (cm-1)
C            efc - rotational partition function exponent
C            wmid  - frequency in center of interval (cm-1)
C            delnu - range of frequency covered by these lines
C
C       OUTPUT
C
C            a - average absorption in this frequency interval
C
C     -------------------------------------------------------------------

       dimension prs(1),tmp(1),col(1),str(1),wid(1),eng(1),a(1)

c      ------------------------------------------------------------------
c                         Internal arrays
c      ------------------------------------------------------------------
       dimension  wdop(nd1),wlor(nd1),hckt(nd1),qrat(nd1),sumstr(nd1)
     &  ,sumsqr(nd1),sumlor(nd1),u(nd1),y(nd1),f(nd1),alor(nd1),w(nd1)
     &  ,v(nd1),g(nd1),adop(nd1),st(nd1,500),dcol(nd1)

       tmp0 = 296.d0
       twopi = two*pi

       if(nlines .eq. 0) then
          do n = 1, ntau
             a(n) = zero 
          end do
          return
       end if

       do n = 1, ntau
           wdop(n) = wmid*dsqrt(two*rkb*tmp(n)/rmass)/cs
           wlor(n) = 1.d-6*prs(n) * dsqrt(tmp0 / tmp(n))
           hckt(n) = hck*(one/tmp(n) - one/tmp0)
           qrat(n) = ((tmp0/tmp(n))**efc)
       end do

C      --------------------------------------------------------------
C          Employ the Curtis-Godson approximation to determine 
C          appropriate averages for the pressure and column
C          abundance.  Use the local temperature as the reference.
C      --------------------------------------------------------------
       do lw = 1, nlines
          do n = 1, ntau
             st(n,lw)=str(lw)*qrat(n)*dexp(-hckt(n)*eng(lw))
          end do
       end do

       do n = 2, ntau
          dcol(n) = col(n)-col(n-1)
       end do

       do n = 1, ntau

c         -----------------------------------------------------------
c                    Calculate sum of line strengths 
c         -----------------------------------------------------------
          suml = zero
          sumw = zero
          do lw = 1, nlines
             s1=st(1,lw)
             sum = s1*col(1)
             do m = 2, n
                s2=st(m,lw)
                sum = sum + half*(s1+s2)*dcol(m)
                s1 = s2
             end do
             suml = suml + sum/uu
             sumw = sumw + dsqrt(sum/uu)
          enddo
          sumstr(n) = suml
          sumsqr(n) = sumw

c         -----------------------------------------------------------
c          Calculate sum of line strengths weighted by Lorentz width
c         -----------------------------------------------------------
          suml = zero
          do lw = 1, nlines
             s1=st(1,lw)*wid(lw)*wlor(1)
             sum = s1*col(1)
             do m = 2, n
                s2=st(m,lw)*wid(lw)*wlor(m)
                sum = sum + half*(s1+s2)*dcol(m)
                s1 = s2
             end do
             suml = suml + dsqrt(sum/uu)
          enddo
          sumlor(n) = suml

       end do

C      --------------------------------------------------------------
C                     Caclulate Lorentz opacity
C      --------------------------------------------------------------
       fac = two/pi
       do n = 1, ntau
          u(n) = fac*(sumstr(n)/sumlor(n))**2
       end do

       fac = one/four/delnu
       do n = 1, ntau
          y(n) = fac*(sumlor(n)**2)/sumstr(n)
       end do

       do n = 1, ntau       
          if(u(n) .le. four) then
             f(n) = u(n)/dsqrt(one+0.515*u(n))
          else
             f(n)=(half+four*u(n)-0.03125/u(n))/dsqrt(twopi*u(n))-one
          endif
       end do

       do n = 1, ntau
          alor(n) = twopi*y(n)*f(n)
       end do
           
C      --------------------------------------------------------------
C                     Cacluale Doppler opacity
C      --------------------------------------------------------------
       fac = four/sqpi
       do n = 1, ntau
          w(n) = fac*(sumstr(n)/sumsqr(n))**2/wdop(n)
       end do

       fac = one/four/delnu
       do n = 1, ntau
          v(n) = fac*wdop(n)*sumsqr(n)**2/sumstr(n)
       end do
     
       do n = 1, ntau
       if(w(n) .le. 1.d-3) then
          g(n) = w(n)
       else if(w(n) .le. four) then
          g(n)=3.8035d0*(dlog((w(n)+3.4097d0)/(w(n)+16.5903d0))
     &        +1.5822d0)
       else
          dum=dlog(sqpi*w(n))
          g(n)=0.1462d0+0.66666667d0*dum**1.5d0+0.3703d0/dum
       endif
       end do

       do n = 1, ntau
          adop(n) = two*v(n)*g(n)
       end do

       do n = 1, ntau
          a(n) = dsqrt(adop(n)**2 + alor(n)**2)
       end do

c       write(unit=*,fmt='(1x,i6,1p4e12.4)')
c     &  nlines,wmid,sumstr(1),sumsqr(1),sumlor(1)

       return
       end