program spiegel c ******************************************************************* c * * c * Spiegel calculates scale dependent radiative damping * c * rates for homogeneous atmospheres. * c * * c * RVY 31 APR 91 * c * * c ******************************************************************* implicit double precision (o-z, a-h) include "parameters.dat" character fname*64 dimension arad(12,3),pcol(12,3),dgn(12,3),rdmol(3),rmcol(3) &,rmass(3),qrot(3),wave(12,3),wvert(50),rate(50,3),xcol(50,3) &,wav(nd4),str(nd4),wid(nd4),eng(nd4),ifund(4),fmix(3),wlortz(3) &,thin(3) c --------------------------------------------------------- c Data for Titan c --------------------------------------------------------- data rmn2 / 28.d0 / data rdn2 / 3.681d0 / data GMm / 4.16d-4 / data rplnt / 2575.d5 / data nlortz/300/ data wlortz/5.2d-2,1.d-1,1.d-1/ data wdel/10.d0/ data ifund/1,5,7,8/ data inmax/10/ data immax/3/ fourpi = four*pi twopi = two*pi do im = 1, immax do in = 1, inmax rate(in,im) = zero end do thin(im) = zero end do c --------------------------------------------------------------- c Vertical (spatial) wavenumber array c --------------------------------------------------------------- do in = 1, inmax wvert(in) = twopi/(2.d5*float(in)) end do c --------------------------------------------------------------- c Mixing ratios c --------------------------------------------------------------- write(unit=*,fmt='(22h Enter mixing ratios >,$)') read (unit=*,fmt=*) (fmix(im),im=1,immax) c --------------------------------------------------------------- c Temperatures, densities, and pressures c --------------------------------------------------------------- write(unit=*,fmt='(33h Enter temperature and pressure >,$)') read (unit=*,fmt=*) tmp,prs den0 = prs/rkb/tmp cp = 2.5d0*rkb*den0 do im = 1, immax do in = 1, inmax xcol(in,im) = fmix(im)*den0/wvert(in) end do end do c ---------------------------------------------------------- c Read in IR vibration transition parameters c ---------------------------------------------------------- call bandread(wave,arad,pcol,dgn,rmass,rdmol,qrot) efc = qrot(1) c ---------------------------------------------------------- c Calculate reduced mass and collision diameters c for hydrocarbon-N2 collisions c ---------------------------------------------------------- do im = 1,immax rdmol(im) = 1.d-8*half*(rdn2+rdmol(im)) rmcol(im) = amu*rmn2*rmass(im)/(rmn2+rmass(im)) end do c --------------------------------------------------------------- c Loop over molecular species c --------------------------------------------------------------- do im = 1, immax c -------------------------------------------------------------- c Read spectral line list c -------------------------------------------------------------- call sread(ifund(im),wav,str,wid,eng,nlines,stot,wmean & ,wstart,wend) write(unit=*,fmt=952) stot,nlines,wstart,wend c write(unit=*,fmt='(28h Enter wstart, wend, wdel > ,$)') c read (unit=*,fmt=*) wstart, wend, wdel c -------------------------------------------------------------- c Coarse frequency loop c -------------------------------------------------------------- op_sum = zero wmin = wstart do while(wmin .lt. wend) wmax = dmin1(wmin+wdel,wend) wavs = half*(wmin+wmax) write(unit=*,fmt=*) wavs c ----------------------------------------------------------- c Calculate Planck function and derivative c ----------------------------------------------------------- bp = bplanck(wavs,tmp) dbp = (bp/tmp)*dplanck(wavs,tmp) c ----------------------------------------------------------- c Prepare to calculate absorption coefficients c ----------------------------------------------------------- tref = 296.d0 wfc = wavs*dsqrt(two*rkb/rmass(im)/amu)/cs wstep = wfc/two wlc = 1.d-6/dsqrt(tref)/wfc wdop = wfc*dsqrt(tmp) wlor = wlc*prs wreach = dmax1(nlortz*wlortz(im)*wlor*wdop & ,five*wdop) hckt = hck*(one/tmp - one/tref) qw = (tref/tmp)**efc/wdop lwstart = min(locate(wav,nlines,wmin-wreach)+1,nlines) lwend = max(locate(wav,nlines,wmin+wreach),1) c -------------------------------------------------------------- c Fine frequency loop c -------------------------------------------------------------- wlast = wmin op_last = zero do while(wlast .lt. wmax) wnow = dmin1(wlast + wstep, wmax) wmid = half*(wnow+wlast) delw = wnow - wlast c ----------------------------------------------------- c Update lwstart and lwend c ----------------------------------------------------- do while(wmid-wav(lwstart) .gt. wreach) lwstart = lwstart + 1 end do do while(wav(lwend)-wmid .le. wreach) lwend = lwend + 1 end do c ------------------------------------------------------ c Calculate contribution of each line to absorption c coefficient c ------------------------------------------------------ suml = zero lw = lwstart do while(lw .le. min(lwend,nlines)) dw = dabs(wmid-wav(lw)) if(dw .le. wreach) then xw = dw / wdop avgt = wid(lw)*wlor call voigt(xw,avgt,vg) opv = vg*str(lw)*dexp(-hckt*eng(lw)) suml = suml + opv end if lw = lw + 1 end do op_now = qw*suml/sqpi c ----------------------------------------------------- c Loop over vertical (spatial) wavenumber c ----------------------------------------------------- op = half*(op_last+op_now) op_sum = op_sum + delw*op fac = delw*(fourpi*dbp/cp) thin(im) = thin(im) + fac*fmix(im)*op*den0 do in = 1, inmax tau = op*xcol(in,im) rate(in,im) = rate(in,im) + fac*wvert(in) & *tau*(one-tau*datan(one/tau)) end do op_last = op_now wlast = wnow c ----------------------------------------------------------- c End fine frequency loop c ----------------------------------------------------------- end do wmin = wmax c ------------------------------------------------------------- c End coarse frequency loop c ------------------------------------------------------------- end do write(unit=*,fmt='(28h Integrated band strength = ,e11.3)') & op_sum write(unit=*,fmt='(14h Thin limit = ,1pe11.3)') thin(im) do in = 1, inmax write(unit=*,fmt=*) 1.e-5/wvert(in), rate(in,im) end do c --------------------------------------------------------------- c End molecular species loop c --------------------------------------------------------------- end do c --------------------------------------------------------------- c Write output to file c --------------------------------------------------------------- write(unit=*,fmt='(25h Enter output file name >,$)') read (unit=*,fmt=900) fname open(unit=60,file=fname,status='new') do in = 1, inmax write(unit=60,fmt=901) wvert(in), (rate(in,im),im=1,immax) end do close(unit=60) stop c --------------------------------------------------------------- c Formats c --------------------------------------------------------------- 900 format(a) 901 format(1x,1p4e12.4) 952 format(' band strength = ',1pe12.4/ & ,' number of lines = ',i8/ & ,' wavenumber range = ',1p2e12.4) end