program opacity c ********************************************************************** c * * c * * c * * c * Roger Yelle * c * 01 Nov 1990 * c * * c ********************************************************************** implicit double precision (o-z, a-h) include "parameters.dat" character fname*64, iret*1 real*8 nld c --------------------------------------------------------- c Arrays : Model atmosphere c --------------------------------------------------------- dimension zalt(nd3),pres(nd3),temp(nd3),dens(nd3,0:4) & ,col(nd3,0:4),fmix(3) c --------------------------------------------------------- c Arrays : IR band parameters c --------------------------------------------------------- & , arad(12,3),pcol(12,3),dgn(12,3),rdmol(3),rmcol(3) & , rmass(3), qrot(3), wave(12,3), wlortz(3) c --------------------------------------------------------- c Arrays : line listing c --------------------------------------------------------- & , wav(nd4),str(nd4),wid(nd4),eng(nd4),nlines(4),nl(4) & , stot(4),opp0(4),wdopref(4),wmean(4),wstart(4),wend(4) & , wstep(4) c --------------------------------------------------------- c Arrays: absorption coefficients c --------------------------------------------------------- & , opf(nd1,nd2),wpf(nd2), wp(20000), op(20000,nd1), otmp(nd2) & ,nloc(nd2),wpec(20000),spec(20000,3),sres(200) c --------------------------------------------------------- c Arrays for R-T calculation c --------------------------------------------------------- & , z(nd1),tmp(nd1),prs(nd1),den(nd1,0:4),clm(nd1,0:4),dtdc(3) c --------------------------------------------------------- c Data for Titan c --------------------------------------------------------- data rmn2 / 28.d0 / data rdn2 / 3.d0 / data gm / 8.95d18 / data GMm / 4.16d-4 / data rplnt / 2575.d5 / data rh / 9.5d0 / data nu /3/ data nlortz/300/ data wlortz/5.2d-2,1.d-1,1.d-1/ data wdel/10.d0/ data pmin / 1.d-5 / data pmax / 1.d+2 / nld = 2.7d19 c write(unit=*, fmt=901) c read (unit=*, fmt=900) fname fname='titan.atm' c write(unit=*,fmt=902) c read (unit=*,fmt=*) cmin,pmax,ntau write(unit=*,fmt='(14h Enter ntau > ,$)') read (unit=*,fmt=*) ntau nmid = (ntau+1)/2 c write(unit=*,fmt=903) c read (unit=*,fmt=*) tref tref = 150.d0 write(unit=*,fmt=904) read (unit=*,fmt=*) otol c ---------------------------------------------------------- c Read in IR vibration transition parameters c ---------------------------------------------------------- call bandread(wave,arad,pcol,dgn,rmass,rdmol,qrot) c ---------------------------------------------------------- c Calculate reduced mass and collision diameters c for hydrocarbon-N2 collisions c ---------------------------------------------------------- do im = 1,3 rdmol(im) = 1.d-8*half*(rdn2+rdmol(im)) rmcol(im) = amu*rmn2*rmass(im)/(rmn2+rmass(im)) end do c ---------------------------------------------------------- c Read in line listing for fundamentals c ---------------------------------------------------------- write(unit=*,fmt='(14h Enter imol > ,$)') read (unit=*,fmt=*) imol if(imol .eq. 1) then jb = 1 else if(imol .eq. 2) then jb = 5 else if(imol .eq. 3) then jb = 7 else if(imol .eq. 4) then jb = 8 end if call sread(jb,wav,str,wid,eng,nlines(imol) & ,stot(imol),wmean(imol),wstart(imol),wend(imol)) wdopref(imol)=dsqrt(two*rkb*tref/rmass(imol)/amu) & *wmean(imol)/cs opp0(imol) = stot(imol)/wdopref(imol) wstep(imol) = wdopref(imol) write(unit=*,fmt=810) wstart(imol),wend(imol) 810 format(' wstart, wend = ',2f9.3) c ---------------------------------------------------------- c Read in Model Atmosphere c ---------------------------------------------------------- call atmread2(fname,pmax,natm,zalt,pres,temp,dens,col) call cgrid(gm,rplnt,pmin,pmax,ntau,zalt,temp,pres,dens & ,col,natm,z,tmp,prs,den,clm) c ---------------------------------------------------------- c Calculate absorption coefficients c ---------------------------------------------------------- write(unit=*,fmt='(20H Enter wmin, wmax > ,$)') read (unit=*,fmt=*) wmin,wmax wmid = half*(wmin+wmax) nmax = nd2 call smt(nlortz,wlortz(imol),qrot(imol),otol,amu*rmass(imol) & ,wmin,wmax,wmid,wstep(imol),wav,str,eng,wid & ,nlines(imol),tmp,prs,clm(1,imol),ntau,nmax,wpf,opf,npf,wp & ,op,np) write(unit=*,fmt=*) ' np = ',np,' npf = ',npf c ------------------------------------------------------------- c Sort absorption coefficients in decreasing order c ------------------------------------------------------------- do n = 1, npf otmp(n) = opf(nmid,n) nloc(n) = n end do call sort(npf,otmp,nloc) do n = 1, npf otmp(n) = opf(1,n) end do do n = 1, npf opf(1,n) = nld*otmp(nloc(n)) end do do n = 1, npf otmp(n) = opf(nmid,n) end do do n = 1, npf opf(nmid,n) = nld*otmp(nloc(n)) end do do n = 1, npf otmp(n) = opf(ntau,n) end do do n = 1, npf opf(ntau,n) = nld*otmp(nloc(n)) end do c ---------------------------------------------------------- c Smear to specified resolution for plotting c ---------------------------------------------------------- write(unit=*,fmt= & '(40H Enter resolution for output spectrum > ,$)') read (unit=*,fmt=*) wres npts = ifix(wres/(wp(2)-wp(1))) npts = 2*npts apts = npts da = one/apts da2 = da*da sum = zero do k = 1, npts/2 sres(k) = (k - half)*da2 sres(npts+1-k) = sres(k) sum = sum + sres(k) + sres(npts+1-k) end do do k = 1, npts sres(k) = sres(k)/sum write(unit=*,fmt=*) k, sres(k) end do do n = 1, np - npts sum0 = zero sum1 = zero sum2 = zero sum3 = zero do k = 1, npts l = k + n - 1 sum0 = sum0 + sres(k)*wp(l) sum1 = sum1 + sres(k)*op(l,1) sum2 = sum2 + sres(k)*op(l,nmid) sum3 = sum3 + sres(k)*op(l,ntau) end do wpec(n) = sum0 spec(n,1) = nld*sum1 spec(n,2) = nld*sum2 spec(n,3) = nld*sum3 end do c ---------------------------------------------------------- c Write output to Disc c ---------------------------------------------------------- write(unit=*,fmt='(28h Enter spectrum file name > ,$)') read (unit=*,fmt='(a)') fname open (unit=30, file=fname, status='new') do n = 1, np-npts write(unit=30,fmt='(1x,f10.4,2x,3e11.3)') & wpec(n),(spec(n,k),k=1,3) end do close(unit=30) write(unit=*,fmt='(23h Enter smt file name > ,$)') read (unit=*,fmt='(a)') fname open (unit=30, file=fname, status='new') wwv = zero do n = 1, npf write(unit=30,fmt='(1x,f10.4,2x,3e11.3)') & wwv,opf(1,n),opf(nmid,n),opf(ntau,n) wwv = wwv + wpf(nloc(n)) end do write(unit=30,fmt='(1x,f10.4,2x,3e11.3)') & wwv,opf(1,npf),opf(nmid,npf),opf(ntau,npf) close(unit=30) stop c ---------------------------------------------------------- c Formats c ---------------------------------------------------------- 900 format(a) 901 format(' Enter model atmosphere file name > ',$) 902 format(' Enter cmin, pmax, and ntau > ',$) 903 format(' Enter reference temperature > ',$) 904 format(' Enter otol > ',$) 905 format(' Enter filename for heating rates > ',$) 906 format(1x,0pf7.1,1x,0pf7.2,2x,1p9e11.3) 907 format(' Enter filename for source functions > ',$) 908 format(1x,0pf7.1,1x,1pe10.3,2x,1p3e11.3) 909 format(' Enter filename for time constants > ',$) 910 format(1x,0pf7.1,1x,0pf7.2,2x,1p5e10.3) end subroutine smt(nlortz,wlortz,efc,otol,rmass,wmin,wmax,wmean & ,wdel,wav,str,eng,wid,nlines,tmp,prs,clm,ntau,nmax,wpp,opp & ,ibinmax,wp,op,np) implicit double precision (o-z, a-h) include "parameters.dat" logical itest, istop c ------------------------------------------------------------- c Input Arrays c ------------------------------------------------------------- dimension tmp(nd1),prs(nd1),clm(nd1),wav(1),str(1),eng(1),wid(1) c ------------------------------------------------------------- c Output Arrays c ------------------------------------------------------------- dimension opp(nd1,nd2), wpp(nd2) c -------------------------------------------------------------- c Internal Arrays c -------------------------------------------------------------- dimension wp(20000),op(20000,nd1),wdop(nd1),wlor(nd1) & ,wreach(nd1),hckt(nd1),qw(nd1),lw1(nd1),lw2(nd1),tpt(nd1) & ,tpb(nd1,nd2) tref = 296.d0 do nt = 1, ntau wdop(nt) = wmean*dsqrt(two*rkb*tmp(nt)/rmass)/cs wlor(nt) = 1.d-6 * prs(nt) * dsqrt(tmp(nt)/tref)/wdop(nt) c wreach(nt) = dmax1(nlortz*wlortz*wlor(nt)*wdop(nt) c & ,five*wdop(nt)) wreach(nt) = 10.d0 hckt(nt) = hck*(one/tmp(nt) - one/tref) qw(nt) = (tref/tmp(nt))**efc/wdop(nt) end do c -------------------------------------------------------------- c Find the range of lines needed to calculate absorption c coefficient at the first frequency c -------------------------------------------------------------- do nt = 1, ntau w1 = wmin - wreach(nt) lw1(nt) = min(locate(wav,nlines,w1)+1,nlines) w2 = wmin + wreach(nt) lw2(nt) = max(locate(wav,nlines,w2),1) end do c -------------------------------------------------------------- c ***** Everything is ready, start frequency loop ***** c -------------------------------------------------------------- lw1nt = min(locate(wav,nlines,wmin)+1,nlines) lw2nt = max(locate(wav,nlines,wmax),1) wlast = wmin ibinmax = 0 iw = 0 do while(wlast .lt. wmax) wnow = dmin1(wlast + wdel, wmax) wmid = half*(wnow+wlast) delw = wnow - wlast opmax = zero iw = iw + 1 wp(iw) = wmid c -------------------------------------------------------- c Find spectral lines which contribute to absorption c coefficient at wmid c -------------------------------------------------------- do nt = 1, ntau do while((wmid-wav(lw1(nt)) .gt. wreach(nt)).and. & (lw1(nt) .lt. nlines)) lw1(nt) = lw1(nt) + 1 end do do while((wav(lw2(nt))-wmid .le. wreach(nt)).and. & (lw2(nt) .lt. nlines)) lw2(nt) = lw2(nt) + 1 end do end do c ----------------------------------------------------------- c Loop backwards through altitudes, stop loop once the c top of atmosphere is reached (nt = 1) or if there are c no spectral lines within wreach of wmid c ----------------------------------------------------------- nt = ntau c do while((nt .ge. 1).and.(lw1(nt).le.lw2(nt))) do while(nt .ge. 1) c --------------------------------------------------------- c Calculate contribution of each line to absorption c coefficient c --------------------------------------------------------- sumo = zero lw = lw1nt do while(lw .le. lw2nt) dw = dabs(wmid-wav(lw)) c if(dw .le. wreach(nt)) then xw = dw / wdop(nt) avgt = wid(lw)*wlor(nt) call voigt(xw,avgt,vg) if(vg .le. zero) then vg = avgt/xw/xw/pi end if opv = vg*str(lw)*dexp(-hckt(nt)*eng(lw)) sumo = sumo + opv c end if lw = lw + 1 end do op(iw,nt) = qw(nt)*sumo/sqpi nt = nt - 1 end do 910 format(' xw = ',1pe12.4,' avgt = ',1pe12.4,' vg = ',1pe12.4) tpt(1) = dexp(-two*op(iw,1)*(clm(2)-clm(1))) do nt = 2, ntau-1 tpt(nt) = dexp(-two*op(iw,nt)*(clm(nt+1)-clm(nt-1))) end do tpt(ntau) = dexp(-two*op(iw,ntau)*(clm(ntau)-clm(ntau-1))) c ----------------------------------------------------------- c Find a home for this frequency c Loop backwards through bins and stop when a home is found c ----------------------------------------------------------- istop = .true. ibin = ibinmax + 1 do while((ibin.gt.1).and. istop ) ibin = ibin - 1 c -------------------------------------------------------- c does this frequency fit in any pre-existing bins ? c -------------------------------------------------------- itest = .true. nt = 0 do while(itest .and. (nt .lt. ntau)) nt = nt + 1 if(dabs(tpt(nt)-tpb(nt,ibin)).gt.otol) itest= .false. end do c -------------------------------------------------------- c If a group is found, add into bin, redefine weight c and average c -------------------------------------------------------- if(itest) then wsum = wpp(ibin) + delw do nt = 1, ntau opp(nt,ibin)=(wpp(ibin)*opp(nt,ibin)+delw*op(iw,nt)) & /wsum end do wpp(ibin) = wsum tpb(1,ibin) = dexp(-two*opp(1,ibin)*(clm(2)-clm(1))) do nt = 2, ntau-1 tpb(nt,ibin) = dexp(-two*opp(nt,ibin) & *(clm(nt+1)-clm(nt-1))) end do tpb(ntau,ibin) = dexp(-two*opp(ntau,ibin) & *(clm(ntau)-clm(ntau-1))) istop = .false. end if end do c ----------------------------------------------------------- c If it didn't fit create a new bin c ----------------------------------------------------------- if(istop) then ibinmax = ibinmax + 1 do nt = 1, ntau opp(nt,ibinmax) = op(iw,nt) end do wpp(ibinmax) = delw do nt = 1, ntau tpb(nt,ibinmax) = tpt(nt) end do end if c ----------------------------------------------------------- c Done with this frequency, update wlast and continue c with loop c ----------------------------------------------------------- wlast = wnow end do np = iw return end