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) C --- LAY if ((wmid .gt. 2809.030 .and. & wmid .lt. 2818.894) .and. & (uu .gt. 0.1 .and. & uu .lt. 0.12)) then idebug=1 write (*,*) 'DEBUG random' else idebug=0 end if if (idebug .eq. 1) then write (*,*) 'str' write(unit=*,fmt='(3e11.4)') str(1),str(2),str(3) write(unit=*,fmt='(A11,2e11.4)') & '..',str(nlines-1),str(nlines) write (*,*) 'wid' write(unit=*,fmt='(3e11.4)') wid(1),wid(2),wid(3) write(unit=*,fmt='(A11,2e11.4)') & '..',wid(nlines-1),wid(nlines) write (*,*) 'eng' write(unit=*,fmt='(3e11.4)') eng(1),eng(2),eng(3) write(unit=*,fmt='(A11,2e11.4)') & '..',eng(nlines-1),eng(nlines) end if C --- LAY 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) C --- LAY if (idebug .eq. 1) then write (*,*) 'tmp' write(unit=*,fmt='(3e11.4)') & tmp(1),tmp(2),tmp(3) write(unit=*,fmt='(A11,2e11.4)') & '..',tmp(ntau-1),tmp(ntau) write (*,*) 'sumstr' write(unit=*,fmt='(3e11.4)') & sumstr(1),sumstr(2),sumstr(3) write(unit=*,fmt='(A11,2e11.4)') & '..',sumstr(ntau-1),sumstr(ntau) write (*,*) 'sumsqr' write(unit=*,fmt='(3e11.4)') & sumsqr(1),sumsqr(2),sumsqr(3) write(unit=*,fmt='(A11,2e11.4)') & '..',sumsqr(ntau-1),sumsqr(ntau) write (*,*) 'wdop' write(unit=*,fmt='(3e11.4)') wdop(1),wdop(2),wdop(3) write(unit=*,fmt='(A11,2e11.4)') & '..',wdop(ntau-1),wdop(ntau) write (*,*) 'w' write(unit=*,fmt='(3e11.4)') w(1),w(2),w(3) write(unit=*,fmt='(A11,2e11.4)') & '..',w(ntau-1),w(ntau) write (*,*) 'g' write(unit=*,fmt='(3e11.4)') g(1),g(2),g(3) write(unit=*,fmt='(A11,2e11.4)') & '..',g(ntau-1),g(ntau) write (*,*) 'v' write(unit=*,fmt='(3e11.4)') v(1),v(2),v(3) write(unit=*,fmt='(A11,2e11.4)') & '..',v(ntau-1),v(ntau) write (*,*) 'adop' write(unit=*,fmt='(3e11.4)') adop(1),adop(2),adop(3) write(unit=*,fmt='(A11,2e11.4)') & '..',adop(ntau-1),adop(ntau) write (*,*) 'alor' write(unit=*,fmt='(3e11.4)') alor(1),alor(2),alor(3) write(unit=*,fmt='(A11,2e11.4)') & '..',alor(ntau-1),alor(ntau) write (*,*) 'a' write(unit=*,fmt='(3e11.4)') a(1),a(2),a(3) write(unit=*,fmt='(A11,2e11.4)') & '..',a(ntau-1),a(ntau) end if C --- LAY return end