subroutine smt(nlortz,wlortz,efc,tol,rmass,wmin,wmax,wmean &,wdel,wav,str,eng,nlines,tmp,prs,clm,ntau,nmax,wpp,opp,ibinmax) c ******************************************************************* c * * c * Purpose: * c * * c * To group frequencies with in a molecular band with similar * c * absorption coefficients at all levels in a model atmosphere * c * into bins thereby reduce the number of frequency integration * c * points used iin calculations of radiation fields and heating * c * rates. This subroutine is based upon the spectral mapping * c * algorithm described by West et al., JQSRT, vol xx, pp xx-yy, * c * 1990. * c * * c * Input Parameters: * c * * c * nlortz - maximum extent of a spectral line Lorentz line * c * widths * c * wlortz - Lorentz line width coefficient in cm-1/ubar * c * efc - exponentail factor in line strength T dependence * c * tol - Tolerance on the tranmission between layers of * c * the atmosphere used in the smt grouping algorithm * c * rmass - The mass of the molecule (gm) * c * wmin - Lower bound on the frequency interval considered * c * wmax - Upper bound on the frequency interval considered * c * wmean - Mean value of the frequency interval, used to * c * calculate Doppler widths * c * wdel - step size for sampling the absorption coefficient * c * wav - array of wavenumbers for spectral lines * c * str - array of spectral line strengths (cm-1/cm2) * c * eng - array of lower state energies for spectral lines * c * nlines - number of spectral lines * c * tmp - array of atmospheric temperatures * c * prs - array of atmospheric pressures * c * clm - array of overhead column abundance * c * ntau - number of atmospheric layers * c * nmax - maximum number of spectral groups * c * * c * * c * Output Variables: * c * * c * wpp - array of weights for spectral groups * c * opp - array of absorption coefficents for spectal groups * c * ibinmax - number of spectral groups * c * * c * * c * 17 August 1990 * c * Roger Yelle * c * * c ******************************************************************* implicit double precision (o-z, a-h) logical itest, istop parameter (zero = 0.d0, half = 5.d-1, one = 1.d0, two = 2.d0 & ,five=5.d0, rkb=1.3807d-16, cs = 2.9979d10 & , sqpi = 1.772453851d0, pi = 3.141592654d0 & , hck = 1.439565d0, nd1 = 101, nd2 = 3000) c ------------------------------------------------------------- c Input Arrays c ------------------------------------------------------------- dimension tmp(1), prs(1), clm(1), wav(1), str(1), eng(1) c ------------------------------------------------------------- c Output Arrays c ------------------------------------------------------------- dimension opp(nd1,1), wpp(1) c -------------------------------------------------------------- c Internal Arrays c -------------------------------------------------------------- dimension op(nd1),wdop(nd1),wlor(nd1),wreach(nd1),avgt(nd1) & ,xvgt(2000,nd1),yvgt(2000,nd1),dxvgt(nd1),nvgt(nd1) & ,hckt(nd1),qrat(nd1),xmax1(nd1),xmax2(nd1) & ,ymax(nd1),lw1(nd1),lw2(nd1),tpt(nd1),tpb(nd1,nd2) tref = 296.d0 c ---------------------------------------------------- c Calculate Voigt functions, and other arrays c ---------------------------------------------------- do nt = 1, ntau wdop(nt) = wmean*dsqrt(two*rkb*tmp(nt)/rmass)/cs wlor(nt) = wlortz * prs(nt) * dsqrt(tref / tmp(nt)) avgt(nt) = wlor(nt) / wdop(nt) wreach(nt) = dmax1(nlortz*wlor(nt),five*wdop(nt)) call vgrid(avgt(nt),xvgt(1,nt),yvgt(1,nt),dxvgt(nt) & ,nvgt(nt)) xmax2(nt) = xvgt(nvgt(nt),nt) xmax1(nt) = xmax2(nt) - dxvgt(nt) ymax(nt) = yvgt(nvgt(nt),nt) hckt(nt) = hck*(one/tmp(nt) - one/tref) qrat(nt) = (tref/tmp(nt))**efc 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 -------------------------------------------------------------- wlast = wmin ibinmax = 0 do while(wlast .lt. wmax) wnow = dmin1(wlast + wdel, wmax) wmid = half*(wnow+wlast) delw = wnow - wlast opmax = zero 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 do while((nt .ge. 1).and.(lw1(nt).le.lw2(nt))) c -------------------------------------------------------- c Find spectral lines which contribute to absorption c coefficient at wmid c -------------------------------------------------------- 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 c -------------------------------------------------------- c Calculate contribution of each line to absorption c coefficient c --------------------------------------------------------- sum = zero lw = lw1(nt) do while(lw .le. lw2(nt)) xw = dabs(wmid - wav(lw)) / wdop(nt) if (xw .ge. xmax1(nt)) then rw = (xmax2(nt)/xw) ** 2 vg = ymax(nt) * rw else i1 = xw/dxvgt(nt) i1 = i1 + 1 i2 = i1 + 1 rr = (xw-xvgt(i1,nt))/(xvgt(i2,nt)-xvgt(i1,nt)) vg = yvgt(i1,nt)+rr*(yvgt(i2,nt)-yvgt(i1,nt)) end if sum=sum+vg*str(lw)*dexp(-hckt(nt)*eng(lw)) lw = lw + 1 end do op(nt) = sum*qrat(nt)/wdop(nt) nt = nt - 1 end do tpt(1) = dexp(-two*op(1)*(clm(2)-clm(1))) do nt = 2, ntau-1 tpt(nt) = dexp(-two*op(nt)*(clm(nt+1)-clm(nt-1))) end do tpt(ntau) = dexp(-two*op(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.tol) 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(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(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 return end