c c----------------------------------------------------------------------- c subroutine mke6300(tn,xo2,xo,xn2,ht,te,xne,xo2p,id,fout,ie6300, + ut,glat,glon,f107d,f107a) ! ! 3/10/94: ! Given tn,te,o2,o,n2,o2+ (all dimensioned id), return volume ! emission rate at 6300A in fout. (from doptuv.f in this lib) ! 3/15/05: adapted for new glbmean (added implicit none, and ! use-associationg from new modules). ! use flds_airglow,only: sr63 implicit none ! ! Args: integer,intent(in) :: id,ie6300 real,intent(in) :: ut,glat,glon,f107d,f107a real,dimension(id),intent(in) :: tn,xo2,xo,xn2,ht,te,xne,xo2p real,dimension(id),intent(out) :: fout real,parameter :: b1d=1.2, rk3=8.e-12, a1d=0.0093, a6300=0.0071 ! ! Local: integer :: k real :: rk1,rk2,rkdr do k=1,id rk1 = 2.e-11 * exp(107.8/tn(k)) rk2 = 2.9e-11 * exp(67.5/tn(k)) if (te(k).ge.1200.) rkdr = 1.6e-7 * (te(k)/300.)**(-0.55) if (te(k).lt.1200.) rkdr = 1.95e-7 * (te(k)/300.)**(-0.7) fout(k) = (b1d*rkdr*xo2p(k)*xne(k)) / + (a1d+rk1*xn2(k)+rk2*xo2(k)+rk3*xo(k))*a6300 if (ie6300.gt.0) fout(k) = fout(k) + sr63(k) ! ! Diffs w/ old code in xne: ! write(6,"('mke6300: k=',i3,' rk1,2=',2e12.4,' te=',e12.4, ! | ' xne=',e12.4,' fout=',e12.4)") k,rk1,rk2,te(k),xne(k), ! | fout(k) enddo end subroutine mke6300 c c----------------------------------------------------------------------- c subroutine mke5577(tn,xo2,xo,xn2,ht,te,xne,xo2p,xo21d, + id,fout,ie5577,ut,glat,glon,f107d,f107a) c c Calculate greenline emission: c c 6/95: Up to five components may be requested: c ie5577(1) > 0 -> o1 recombination (original) c (need tn,xo2,xo,xn2) c ie5577(2) > 0 -> dissociative recombination of o2+ c (need te,xo2p,xne,xo21d) c ie5577(3) > 0 -> photoelectron impact c (need tn,ht,xo2,xo,xn2,xo21d and solred) c ie5577(4) > 0 -> airglow c (need tn,ht,xo2,xo,xn2,xo21d and solred) c ie5577(5) > 0 -> photo dissoc of o2 by solar lyman-beta c (qlyb returned by solred) c If solred is called, also need iyd, ut, glat, glon, f107d, f107a c Solred returns sr63 (see e6300), qn2ps, and qops c use flds_ionzrt,only: qops,qn2ps,qlyb implicit none ! ! Args: integer,intent(in) :: id,ie5577(5) real,intent(in) :: ut,glat,glon,f107d,f107a real,dimension(id),intent(in) :: tn,xo2,xo,xn2,ht,te,xne, | xo2p,xo21d real,dimension(id),intent(out) :: fout ! ! Local: integer :: k,kk real :: po1s1,po1s2,po1s3,po1s4,po1s5,o1oprat,rkdr,rk1,rk5,rk12, | qyld,xn2p real,parameter :: a5=1.18, a6=1.35, cp=15., cpp=211. c c O(1s)/O+ ratios and corresponding hts, obtained from Marianna Shepard c 6/95. This is used in po1s4 (component 4) calculation below. c integer,parameter :: no1op=48 real :: o1op(no1op) = +(/0.1377984, 0.1344115, 0.1306053, 0.1256085, 0.1185497, + 0.1175393, 0.1165371, 0.1155426, 0.1145558, 0.1133328, + 0.1121214, 0.1109214, 0.1092618, 0.1076284, 0.1057982, + 0.1037934, 0.1016451, 0.09900434, 0.09598991, 0.09234733, + 0.08779908, 0.08215878, 0.07523289, 0.07433964, 0.07339706, + 0.07248987, 0.07148170, 0.07047267, 0.06947681, 0.0683660, + 0.06724124, 0.06598753, 0.06464776, 0.06320941, 0.0616136, + 0.05986323, 0.05803178, 0.05598638, 0.05314200, 0.05069219, + 0.04738312, 0.04701184, 0.04661626, 0.04618508, 0.04574585, + 0.04538802, 0.04496423, 0.04443988/) real :: o1opht(no1op) = +(/100.2, 101.3, 102.6, 104.4, 107.1, 107.5, 107.9, 108.3, + 108.7, 109.2, 109.7, 110.2, 110.9, 111.6, 112.4, 113.3, + 114.3, 115.6, 117.2, 119.3, 122.2, 126.6, 134.1, 135.2, + 136.4, 137.6, 139.0, 140.5, 142.1, 144.0, 146.0, 148.3, + 150.9, 153.9, 157.4, 161.6, 166.6, 173.0, 181.2, 192.6, + 210.1, 212.4, 214.9, 217.5, 220.4, 223.4, 226.8, 230.4/) c c Column loop: do k=1,id c c po1s1: o1 recombination (need tn,xo2,xo,xn2) c The greenline volume emission rates due to the recombination source, which c dominates in the 85 to 110 km region at night. c (from Ian McDade, (mcdade@windic.yorku.ca) Mar 7, 1994 mail c communication to Roble) c c A5 k1 [O]^3 {[N2]+[O2]} c V1S = ----------------------------- (photons cm-3 sec-1) c {A6 + k5[O2]} {C'[O2] + C"[O]} c c where [O], [O2] and [N2] are number densities cm-3 c A5 = 1.18 sec-1 c A6 = 1.35 sec-1 c k1 = 4.7E-33 (300/T)^2 cm6 sec-1 c k5 = 4.0E-12 exp(-865/T) cm3 sec-1 c C' = 15 dimensionless c C" = 211 dimensionless c c po1s1: o1 recombination c (need tn,o,o2,n2) c po1s1 = 0. if (ie5577(1).gt.0) then rk1 = 4.7e-33*(300./tn(k))**2 rk5 = 4.0e-12*exp(-865./tn(k)) po1s1 = (a5 * rk1 * xo(k)**3 * (xn2(k)+xo2(k))) / + ((a6 + rk5*xo2(k)) * (cp*xo2(k) + cpp*xo(k))) ! write(6,"('mke5577: k=',i3,' rk1,rk5,po1s1=',3e12.4)") ! | k,rk1,rk5,po1s1 ! write(6,"('mke5577: k=',i3,' a5,6=',2e12.4,' cp,cpp=', ! | 2e12.4)") k,a5,a6,cp,cpp ! write(6,"('mke5577: k=',i3,' xo2,xo,xn2=',3e12.4)") ! | k,xo2(k),xo(k),xn2(k) endif c c po1s2: dissociative recombination of o2+ c (need te,xo2p,xne,xo21d) c po1s2 = 0. if (ie5577(2).gt.0) then if (te(k).ge.1200.) rkdr = 1.6e-7 * (te(k)/300.)**(-0.55) if (te(k).lt.1200.) rkdr = 1.95e-7 * (te(k)/300.)**(-0.7) c qyld = 0.12 ! average qyld = 0.04 ! 6/28: experimental lowest value po1s2 = qyld*rkdr*xo2p(k)*xne(k) endif c c po1s3: photoelectron impact c (need tn,ht,xo2,xo,xn2,xo21d and qn2ps from solred) c po1s3 = 0. if (ie5577(3).gt.0) then xn2p = 0.3896*qn2ps(k) / + (2.8e-11*xo(k)+2.3e-11*xo2(k)+0.77) po1s3 = 0.75*2.8e-11*xo(k)*xn2p endif c c po1s4: airglow (from Marianna Shepard) c (need tn,ht,xo2,xo,xn2,xo21d and qops from solred) c o1oprat is O(1s)/O+ ratio corresponding to height, as given by c Marianna Shepard (she ran msis90 to get these numbers): c po1s4 = 0. if (ie5577(4).gt.0) then if (ht(k).lt.o1opht(1)) then o1oprat = 0. elseif (ht(k).gt.o1opht(no1op)) then o1oprat = o1op(no1op) else do kk=1,no1op-1 if (ht(k).ge.o1opht(kk).and.ht(k).le.o1opht(kk+1)) then o1oprat = o1op(kk) if (o1opht(kk+1)-ht(k).lt.ht(k)-o1opht(kk)) + o1oprat = o1op(kk+1) endif enddo endif po1s4 = o1oprat*qops(k) endif c c po1s5: photo dissociation of o2 by solar lyman-beta: c (glby returned by solred) c po1s5 = 0. if (ie5577(5).gt.0) po1s5 = qlyb(k) c c Sum the components: c if (ie5577(2).gt.0.or.ie5577(3).gt.0.or.ie5577(4).gt.0.or. + ie5577(5).gt.0) then rk12 = 4.e-12*exp(-865.*tn(k)) fout(k) = (po1s2 + po1s3 + po1s4 + po1s5)*1.18 / + (1.35+rk12*xo2(k)+2.e-14*xo(k)+2.6e-10*xo21d(k)) + po1s1 elseif (ie5577(1).gt.0) then fout(k) = po1s1 endif enddo end subroutine mke5577 c c----------------------------------------------------------------------- c subroutine mkeo200(tn,xo2,xo,xn2,id,fout) c c 3/9/94: c Given tn, xo2,xo,xn2 (species in cm-3), dimensioned id, return c O2 (0-0) band emission: c (from Ian McDade, (mcdade@windic.yorku.ca) Mar 7, 1994 mail c communication to Roble) c c The nighttime O2 Atmospheric (0-0) band emission rates in the region may be c similarly modelled using: c c A1 k1 [O]^2 {[N2]+[O2]} [O2] c Vat = --------------------------------------------------(photons cm-3 sec-1) c {A2 + k2a[O2] + k2b[N2] + k2c[O]} {D'[O2] + D"[O]} c c where [O], [O2] and [N2] are number densities cm-3 c A1 = 0.079 sec-1 c A2 = 0.083 sec-1 c k1 = 4.7E-33 (300/T)^2 cm6 sec-1 c k2a = 4.0E-17 cm3 sec-1 c k2b = 2.2E-15 cm3 sec-1 c k2c = 8.0E-14 cm3 sec-1 c D' = 6.6 dimensionless c D" = 19 dimensionless c implicit none ! ! Args: integer,intent(in) :: id real,dimension(id),intent(in) :: tn,xo2,xo,xn2 real,dimension(id),intent(out) :: fout ! ! Local: integer :: k real :: rk1 real,parameter :: a1=.079, a2=.083, rk2a=4.e-17, rk2b=2.2e-15, + rk2c=8.e-14, dp=6.6, dpp=19. c do k=1,id rk1 = 4.7e-33*(300./tn(k))**2 fout(k) = (a1*rk1*xo(k)**2 * (xn2(k)+xo2(k)) * xo2(k)) / + ((a2 + rk2a*xo2(k) + rk2b*xn2(k) + rk2c*xo(k)) * + (dp*xo2(k)+dpp*xo(k))) enddo end subroutine mkeo200 c c----------------------------------------------------------------------- c subroutine mkeco215(tn,xo,xco2,id,fout) c c Return 15 micron co2 emission (from John Wise): c c N (15 mic) = (5.94E-26) Sqrt[Tk] exp[-960/Tk] [O] [CO2] c ----------------------------------------- c (4 Pi) (1.28 + 3.5E-13 Sqrt[Tk] [O]) c c The 15 micron term is only the O-CO2 collisional component c but it accounts for at least 80% of the radiance above 110 km c implicit none ! ! Args: integer,intent(in) :: id real,dimension(id),intent(in) :: tn,xo,xco2 real,dimension(id),intent(out) :: fout ! ! Local: integer :: k real,parameter :: pi=3.14156 c do k=1,id fout(k) = (5.94e-26*sqrt(tn(k))*exp(-960./tn(k))*xo(k)*xco2(k))/ + (4.*pi*(1.28+3.5e-13*sqrt(tn(k))*xo(k))) enddo return end c c----------------------------------------------------------------------- c subroutine mkeno53(tn,xo,xno,id,fout) c c Return 5.3 micron NO emission (from John Wise): c c N (5.3 mic) = (2.63E-22) exp[-2715/Tk] [O] [NO] c --------------------------------- c (4 Pi) (10.78 + 6.5E-11 [O]) c implicit none ! ! Args: integer,intent(in) :: id real,dimension(id),intent(in) :: tn,xo,xno real,dimension(id),intent(out) :: fout real,parameter :: pi=3.14156 ! ! Local: integer :: k c do k=1,id fout(k) = (2.63e-22*exp(-2715./tn(k))*xo(k)*xno(k)) / + (4.*pi*(10.78+6.5e-11*xo(k))) enddo return end c c------------------------------------------------------------------ c subroutine getohv(fohv,fohb,integ,fohv_integ,fohb_integ) c c Return fohv(nohvalt) as oh at mxohvlev (10) vibrational levels c (first level is ground state) c use params,only: nlev use flds_atmos,only: tn,xno2,xno,xnn2,xnh,xno3,xnho2,xnoh, | rho use flds_modelz,only: zpht use cons,only: pi implicit none integer,parameter :: mxohvlev=10,nohvalt=91 integer,parameter :: mxohvband=60 ! ! Args: integer,intent(in) :: integ real,intent(out) :: fohv(mxohvlev,nohvalt) ! ohrad output (cm3) real,intent(out) :: fohb(mxohvband,nohvalt) ! ohrad output (kilo-Rayleighs) real,intent(out) :: fohv_integ(mxohvlev) ! ht-integ output of ohrad real,intent(out) :: fohb_integ(mxohvband) ! ht-integ output of ohrad ! ! Local: integer :: i,ier,ip real,parameter :: spval=1.e36 real,parameter :: ohv_kmbot=30.,ohv_kmtop=120. real :: fohv1(nohvalt),ohv_alt(nohvalt) real :: + tn_ohv(nohvalt), xo2_ohv(nohvalt), xo_ohv(nohvalt), + xn2_ohv(nohvalt), xh_ohv(nohvalt), xo3_ohv(nohvalt), + xho2_ohv(nohvalt), xoh_ohv(nohvalt), rho_ohv(nohvalt) integer :: ideltav=2 integer :: ibohv_watts=0 ! ! External: real,external :: quadrat ! (util.F) c do i=1,nohvalt ohv_alt(i) = ohv_kmbot + float(i-1) enddo c c Interpolate required fields to the ohrad altitude grid: c subroutine intloc(gcmin,gcmht,nlev,hts,nhts,loght, c + gcmout,idim1,ndim1,idim2,ier,spval,iprnt) c call intloc(tn,zpht,nlev,ohv_alt,nohvalt,0,tn_ohv, + 1,1,nohvalt,ier,spval,0) call intloc(xno2,zpht,nlev,ohv_alt,nohvalt,1,xo2_ohv, + 1,1,nohvalt,ier,spval,0) call intloc(xno,zpht,nlev,ohv_alt,nohvalt,1,xo_ohv, + 1,1,nohvalt,ier,spval,0) call intloc(xnn2,zpht,nlev,ohv_alt,nohvalt,1,xn2_ohv, + 1,1,nohvalt,ier,spval,0) call intloc(xnh,zpht,nlev,ohv_alt,nohvalt,1,xh_ohv, + 1,1,nohvalt,ier,spval,0) call intloc(xno3,zpht,nlev,ohv_alt,nohvalt,1,xo3_ohv, + 1,1,nohvalt,ier,spval,0) call intloc(xnho2,zpht,nlev,ohv_alt,nohvalt,1,xho2_ohv, + 1,1,nohvalt,ier,spval,0) call intloc(xnoh,zpht,nlev,ohv_alt,nohvalt,1,xoh_ohv, + 1,1,nohvalt,ier,spval,0) call intloc(rho,zpht,nlev,ohv_alt,nohvalt,1,rho_ohv, + 1,1,nohvalt,ier,spval,0) c c Call Makhlouf's code for oh: c All 10 levels are returned in fohv(mxohvlev,nohvalt), including c ground state, which is fohv(1,:). c c subroutine ohrad(kmbot,kmtop,mvl,nkm,TA,SO2,SO,SN2,SH,SO3,SHO2, c + SOH0,SM,ideltav,dohv,bohv,iwatts,spval,iprint) c call ohrad(ifix(ohv_kmbot),ifix(ohv_kmtop),mxohvlev,nohvalt, + tn_ohv,xo2_ohv,xo_ohv,xn2_ohv,xh_ohv,xo3_ohv,xho2_ohv,xoh_ohv, + rho_ohv,ideltav,fohv,fohb,ibohv_watts,spval,0) c c Save ht-integrated ohv: c real function quadrat(parm,hts,npts,spval) c real hts(npts),parm(npts),htcm(npts) c if (integ.gt.0) then do ip=1,mxohvlev fohv1(:) = fohv(ip,:) fohv_integ(ip) = quadrat(fohv1,ohv_alt,nohvalt,spval) enddo do ip=1,mxohvband fohv1(:) = fohb(ip,:) fohb_integ(ip) = quadrat(fohv1,ohv_alt,nohvalt,spval) enddo c c 9/11/96: fohb_integ units are kilo-rayleighs if ibohv_watts=0, or c watts/cm2-str if ibohv_watts=1 (see also ohrad.f in lib/util) c if (ibohv_watts.le.0) then do ip=1,mxohvband fohb_integ(ip) = fohb_integ(ip) * 1.e-9 ! KR enddo else do ip=1,mxohvband fohb_integ(ip) = fohb_integ(ip) / (4.*pi) ! watts/cm2-str enddo endif endif return end !----------------------------------------------------------------------- subroutine intloc(gcmin,gcmht,nlev,hts,nhts,loght, + gcmout,idim1,ndim1,idim2,ier,spval,iprnt) c c Interpolate idim1 values at nhts hts (all values assumed c to be at some tgcm grid lat/lon) c c On input: c gcmin(idim1,nlev) = Tgcm parameter to interpolate c gcmht(idim1,nlev) = Tgcm heights c hts(nhts) = heights at which to do interpolation c loght = 0 or 1 for linear or log interpolation c spval = special value is returned in gcmout if c hts(ih) > max height in gcmht or < min height in gcmht c iprnt = 0 or 1 for print out if spval is used (height out of range) c c On output: c gcmout(idim1,idim2) is returned with interpolated data c ier = 0 if not error has occurred, > 0 if error has occurred c implicit none ! ! Args: integer,intent(in) :: nlev,nhts,loght,idim1,ndim1,idim2,iprnt integer,intent(out) :: ier real,intent(in) :: spval real,intent(in) :: gcmin(idim1,nlev),gcmht(idim1,nlev),hts(nhts) real,intent(out) :: gcmout(idim1,idim2) ! ! Local: integer :: k,id,ih,ialt1,ialt2 real :: htmax,htmin,exparg,f1 c c Check idim2: ier = 0 if (idim2.lt.nhts) then write(6,"('>>> intloc: bad idim2=',i3,' (is < nhts=',i3,')')") + idim2,nhts ier = 1 return endif c c First dimension loop: do 100 id = 1,ndim1 c c Find max height: htmax = -1.e36 htmin = 1.e36 do 110 k=1,nlev htmax = amax1(htmax,gcmht(id,k)) htmin = amin1(htmin,gcmht(id,k)) 110 continue c c Height loop: do 120 ih=1,nhts if (hts(ih).gt.htmax.or.hts(ih).lt.htmin) then if (iprnt.gt.0) then write(6,"('>>> intloc: id ih=',2i4,' height=',e12.4, + ' outside range: htmin max=',2e12.4,' returning ',e12.4)") + id,ih,hts(ih),htmin,htmax,spval endif gcmout(id,ih) = spval ier = 2 goto 120 c stop endif c c Bracket desired height: do 130 k=1,nlev-1 if (hts(ih).ge.gcmht(id,k).and.hts(ih).le.gcmht(id,k+1))then ialt1 = k ialt2 = k+1 goto 131 endif 130 continue if (iprnt.gt.0) + write(6,"('>>> intloc: could not find index to height=', + e12.4,' returning spval')") hts(ih) gcmout(id,ih) = spval goto 120 131 continue c c Do log interpolation: if (loght.eq.1) then if (gcmin(id,ialt1).eq.spval.or.gcmin(id,ialt2).eq.spval.or. + gcmin(id,ialt2).le.1.e-20.or.gcmin(id,ialt1).le.1.e-20 + .or.gcmht(id,ialt2)-gcmht(id,ialt1).le.1.e-20) then gcmout(id,ih) = spval goto 120 endif if (gcmin(id,ialt2).gt.0.) then exparg = (alog(gcmin(id,ialt2) / gcmin(id,ialt1)) / + (gcmht(id,ialt2)-gcmht(id,ialt1))) * + (hts(ih)-gcmht(id,ialt1)) gcmout(id,ih) = gcmin(id,ialt1)*exp(exparg) else gcmout(id,ih) = spval endif c c Do linear interpolation: else if (gcmht(id,ialt2).eq.spval.or.gcmht(id,ialt1).eq.spval.or. + gcmht(id,ialt2)-gcmht(id,ialt1).le.1.e-20) then gcmout(id,ih) = spval goto 120 endif f1 = (hts(ih)-gcmht(id,ialt1)) / + (gcmht(id,ialt2)-gcmht(id,ialt1)) gcmout(id,ih) = f1*gcmin(id,ialt2)+(1.-f1)*gcmin(id,ialt1) endif 120 continue c c End first dimension loop: 100 continue return end