! module solarheat use params,only: nlev use input,only: f107,f107a use cons,only: cgm,eff,re,dtor,pi,erg implicit none integer,parameter :: lmax=59 integer,parameter :: ltorr=37 ! ! sigabs and sigion were formerly in crostrf.h: real :: | sigabs(9,lmax), ! o,o2,n2,co,co2,o3,co2,h2o,no | sigion(13,lmax) ! o,o2,n2,co,co2,o4s,o2d,o2p,o4p,o2q,he,n,h ! ! Fields formerly in common/euv/: real,dimension(ltorr) :: | wleuv1 ,wleuv2 ,sigao ,sigao2 ,sigan2 ,sigio ,sigio2 , | sigin2 ,brop4s ,brop2d ,brop2p ,sigop2p ,sigop2d ,sigop4s , | sigin ,brn2np ,bro2op ! real :: xj762 ! formerly in /solarj/ (solar.h) Used by comp.F. ! ! Function expo checks exp() arg and flushes the result to 0 if arg ! is too small, or to largest value allowed if arg is too big, avoiding ! an fpe underflow or overflow. This function is used only if EXPO ! is set (i.e., when fpe traps are set, see Makefile). Use of expo() ! does slow the code down considerably. ! real,external :: expo ! util.F ! contains !----------------------------------------------------------------------- subroutine solheat(ut) use flds_ionzrt,only: | qop,qo2p,qn2p,qnp,qop2p,qop2d,qop4s,qn2ps,qops,qlyb,qtin, | dinolya,qian,qpro,qixray,qnspe,qiaur,qia,qti,qnop,qcr use flds_dissoc,only: | diop ,dio2p ,din2p ,dinp ,diop2p ,diop2d ,diop4s, | do2euv ,do2lya ,do2src ,do2srb ,do2hz ,do2t, | dh2oeuv ,dh2olya ,dh2osrb ,dh2osrc ,dh2ot , | dco2euv ,dco2lya ,dco2src ,dco2srb ,dco2t , | dch4lya ,dch4euv ,dch4b ,dch4t ,dcoeuv , | dnoeuv ,pnoeuv ,dnosrb , | do3euv ,do3lya ,do3src ,do3srb ,do3hz ,do3har ,do3hug, | do3chap ,do3t ,diop , | sjmlya ,sjo2lya ,sjnmlya use flds_heat,only: | qsrc,qsrb,qnlya,qhar,qhug,qchap,qherz, | qproh,qaurh,qnight,qn,qnpe use flds_airglow,only: | sr63 use flds_col,only: | cnn2,cno2,cno,clno3,schap,txno,txno2,txnn2,txno3, | yeesrcj,yeesrch use flds_atmos,only: | xnn2,xno2,xno,xnno,xno3,rho,sht,xnh2o,xnch4,xnco2,xn4s,xno21d,tn use flds_modelz,only: zp,zpht use aurora,only: cosray,qaurora,protonp,proton use input,only: iaurora,isolproton,istep,nstep ! ! Args: real,intent(in) :: ut ! ! Local: integer :: k,l,ll,itimet,iscale real :: hlybr,fexvir,hlya,heiew,xuvfac,exray,fduxlya,xmlw, | a1d,quench,a6300,tauvn,zi,alt,xp,y,tt,yerf,chaps,scxno2, | xj,eps,sfac,tau,ctau,taur,taup,taur1,taur2,taur3,ctaup, | ctaur,ctaur1,ctaur2,ctaur3,rsp,rsq,cspi,cspj,wl1,sgch4, | alamav,essr,taubt,dum,qo2s1,qo2s2,aureff,alp,eaur,fluxar real :: thng,thng1,thng2,thng3,thng4,thng5,thng6,thng7,thnh real,dimension(lmax) :: wave1,wave2,sflux real,dimension(lmax) :: brn2,bro2,brn2s,bro2s real,dimension(nlev) :: xo21dl,qo2s,xno2s,xxn,col real,dimension(3) :: sigxo2,sigxn2,sigxo,exflxr,effx,xlam,xrflx real,dimension(13) :: sgh2o,sgco2 real :: fnuxlya = 2.e+9 real :: sa(3,3),si(3,3),al(3) real :: effpe = 0.02 real :: euvflx(ltorr) real,dimension(nlev) :: als,tdena,tu,ab,zprs ! for solar protons character(len=80) :: char80 ! sgh2o = (/5.0 ,5.0 ,5.0 ,3.0 ,1.5 ,0.8 ,0.8 ,1.1 ,0.,0.,0.,0.,0./) sgco2 = (/0.5 ,1.0 ,1.5 ,3.0 ,4.0 ,5.5 ,5.0 ,5.0 ,0.,0.,0.,0.,0./) al = (/1.E+7, 5.E+7, 5.E+7/) si(:,1) = (/0., 1.0e-18, 0./) si(:,2) = (/23.11e-18, 22.e-18, 10.24e-18/) si(:,3) = (/11.61e-18, 16.e-18, 8.40e-18/) sa(:,1) = (/0., 1.6e-18, 0./) sa(:,2) = (/23.11e-18, 22.e-18, 10.24e-18/) sa(:,3) = (/11.61e-18, 16.e-18, 8.4e-18/) sigxo2 = (/9.28, 4.4, 1.64/) sigxn2 = (/5.45, 2.35, 0.967/) sigxo = (/4.64, 2.0, 0.87/) exflxr = (/0.50, 0.4, 0.1/) effx = (/3., 5., 7./) xlam = (/70., 50., 35./) brn2s = (/ | 0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0., | 0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0., | 0.01, 0.04, 0.04, 0.03, 0.05, 0.05, | 0.15, 0.20, 0.20, 0.25, 0.32, 0.34, | 0.36, 0.36, 0.36, 0.36, 0.36, 0.36, 0.36,0.36,0.36,0.36,0.36/) bro2s = (/ | 0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0., | 0.,0.,0.,0.,0.,0.,0.,0.,0.,0., | .025, .036, .045, .120, .155, .189, .230, 0.20, 0.20, 0.20, | 0.23, 0.25, 0.29, | 0.33, 0.33, 0.33, 0.33, 0.33, 0.33, 0.33,0.33,0.33,0.33,0.33, | 0.33, 0.33, 0.33, 0.33, 0.33/) ! itimet = 1 ! itimet may not be necessary (see sub ssflux) iscale = 0 hlybr = 0. fexvir = 0. hlya = 1.e+11*(0.5839+0.3554*sqrt(f107a)+0.1730*(sqrt(f107) | -sqrt(f107a))) heiew = 0. xuvfac = 2.0 - (f107-68.0) / (243.0-68.0) if (xuvfac < 1.0) xuvfac = 1.0 ! ! Define fields formerly in common/euv/: call seteuv ! ! THE FLAG ISCALE CHOOSES THE FLUX MODELING METHOD ! FROM THE SUBROUTINE SSFLUX. ISCALE=0 DOES SOLOMON'S, ! ISCALE=1 DOES A LIN. INTERP.,AND ISCALE=2 DOES TOBISKA'S ! ! ISCALE =0 for Hinteregger contrast ratio method ! =1 for Hinteregger linear interpolation ! =2 for Tobiska EUV91 model ! =3 for Woods & Rottman 10 Nov. 1988 measurement ! =4 for Woods & Rottman 20 Jun. 1989 measurement ! F107 daily 10.7 cm flux (1.E-22 W m-2 Hz-1) ! F107A 81-day centered average 10.7 cm flux ! HLYBR ratio of H Ly-b 1026A flux to solar minimum value (optional) ! FEXVIR ratio of Fe XVI 335A flux to solar minimum value (optional) ! HLYA H Lyman-alpha flux (photons cm-2 s-1) (optional) ! HEIEW He I 10830A equivalent width, (milliAngstroms) (optional) ! XUVFAC factor for scaling flux 16-250A (optional) ! WAVE1 longwave bound of spectral intervals (Angstroms) ! WAVE2 shortwave bound of intervals (= WAVE1 for indiv. lines) ! SFLUX scaled solar flux returned by subroutine (photons cm-2 s-1) ! call ssflux(iscale,f107,f107a,hlybr,fexvir,hlya, | heiew,xuvfac,wave1,wave2,sflux,itimet) ! ! 11/17/04 btf: no diffs with old code: ! call fprint8('After ssflux:',lmax, ! | wave1,wave2,sflux,dum,dum,dum,dum,dum, ! | 'wave1','wave2','sflux',' ',' ',' ',' ',' ') ! do l=1,lmax brn2(l)=brn2s(l) bro2(l)=bro2s(l) enddo ! l=1,lmax ! call euvac(euvflx) ! write(6,"('after euvac: ltorr=',i4,' euvflx=',/,(6e12.4))") ! | ltorr,euvflx ! do l=1,ltorr ! 1,37 ll = l+15 ! 16,52 sflux(ll) = euvflx(l) wave1(ll) = wleuv2(l) wave2(ll) = wleuv1(l) ! sigabs(1,ll) = sigao(l) sigabs(2,ll) = sigao2(l) sigabs(3,ll) = sigan2(l) ! sigion(1,ll) = sigio(l) sigion(2,ll) = sigio2(l) sigion(3,ll) = sigin2(l) sigion(6,ll) = sigop4s(l) sigion(7,ll) = sigop2d(l) sigion(8,ll) = sigop2p(l) ! brn2(ll) = brn2np(l) bro2(ll) = bro2op(l) enddo ! ! call fprint8('After euvac:',lmax, ! | sflux,wave1,wave2,sigabs(1,:),sigabs(2,:),sigabs(3,:),dum,dum, ! | 'sflux','wave1','wave2','sigabs(1)','sigabs(2)','sigabs(3)', ! | ' ',' ') ! ! COMPLETELY QUIET SUN 0.1 ERGS/CM/S ! QUIET SUN LOW 0.2 ! QUIET SUN HIGH 0.4 ! ACTIVE 0.6 ! DISTURBED FLARES 1.0 ! EXCEPTIONAL FLARE > 1.0 ! ! EXRAY=0.5+1.*(F107-67.)/176. ! BASED ON SNOE STUDY OF XRAYS ! exray = 0.4+0.6*(f107-70.)/240. ! old model (commented) exray = 0.6+0.6*(f107-70.)/240. ! old model (used) ! exray = 0.6+2.0*(f107-70.)/240. ! Roble tuned ! exray = 1.5 ! Roble tuned ! fduxlya=sflux(12) ! fduxlya=sminflx(8) ! fduxlya=smaxflx(8) ! xmlw = 48. xj762 = 5.35e-9*cgm quench = 2.3e-11 a1d = 6.81e-3 ! local (not the one in flds_ratecoef) a6300 = 5.15e-3 ! ! Init ion production (module flds_ionzrt): qop = 0. qo2p = 0. qn2p = 0. qnp = 0. qop2p = 1.e-20 qop2d = 1.e-20 qop4s = 1.e-20 qn2ps = 1.e-20 qops = 1.e-20 qlyb = 1.e-20 qtin = 1.e-20 ! ! Init ion dissociation (module flds_dissoc): diop = 1.e-20 dio2p = 1.e-20 din2p = 1.e-20 dinp = 1.e-20 diop2p = 1.e-20 diop2d = 1.e-20 diop4s = 1.e-20 ! ! Init O2: qsrc = 1.e-20 qsrb = 1.e-20 qnlya = 1.e-20 qhar = 1.e-20 qhug = 1.e-20 qchap = 1.e-20 qherz = 1.e-20 ! do2euv = 1.e-20 do2lya = 1.e-20 do2src = 1.e-20 do2srb = 1.e-20 do2hz = 1.e-20 do2t = 1.e-20 ! ! Init H20: dh2oeuv = 1.e-20 dh2olya = 1.e-20 dh2osrb = 1.e-20 dh2osrc = 1.e-20 dh2ot = 1.e-20 ! ! Init CO2: dco2euv = 1.e-20 dco2lya = 1.e-20 dco2src = 1.e-20 dco2srb = 1.e-20 dco2t = 1.e-20 ! ! Init CH4: dch4lya = 1.e-20 dch4euv = 1.e-20 dch4b = 1.e-20 dch4t = 1.e-20 ! ! Init CO: dcoeuv = 1.e-20 ! ! Init NO: dnoeuv = 1.e-20 pnoeuv = 1.e-20 dinolya = 1.e-20 dnosrb = 1.e-20 ! ! Init O3: do3euv = 1.e-20 do3lya = 1.e-20 do3src = 1.e-20 do3srb = 1.e-20 do3hz = 1.e-20 do3har = 1.e-20 do3hug = 1.e-20 do3chap = 1.e-20 do3t = 1.e-20 ! ! Init ionization rates: qian = 1.e-20 qiaur = 1.e-20 qpro = 1.e-20 qaurh = 1.e-20 qproh = 1.e-20 qixray = 1.e-20 qnspe = 1.e-20 xo21dl = 1.e-20 diop = 1.e-20 dio2p = 1.e-20 ! ! Airglow: sr63 = 1.e-20 ! ! Aurora: qia(:,:) = 0. qti = 0. ! ! Nighttime ionization source do k=1,nlev do l=1,3 ! was DO 2 tauvn=exp(-(sa(1,l)*cnn2(k)+sa(2,l)*cno2(k)+sa(3,l)*cno(k))) thng1=al(l)*si(1,l)*xnn2(k)*tauvn qn2p(k)=qn2p(k)+thng1*0.86 ! qn2p def 1 qnp(k)=qnp(k)+thng1*0.14 thng2=al(l)*si(2,l)*xno2(k)*tauvn qo2p(k)=qo2p(k)+thng2*0.67 thng3=al(l)*si(3,l)*xno(k)*tauvn qop4s(k)=qop4s(k)+thng3+thng2*0.33 qnspe(k)=qnspe(k)+thng1*1.57 enddo ! l=1,3 ! write(6,"('qn2p def 1: k=',i4,' qn2p(k)=',e12.4)") k,qn2p(k) ! ! Nighttime ionization source for no: qnop(k)=fnuxlya*2.e-18*xnno(k)*exp(-8.e-21*cno2(k)) ! ! Total nighttime ionization source: qtin(k)=qnop(k)+qn2p(k)+qnp(k)+qo2p(k)+qop4s(k) qnight(k)=qtin(k)*35.*eff*1.602e-12/rho(k) qn(k)=qnight(k) xxn(k)=xno3(k) enddo ! k=1,nlev ! ! All apparently ok here (8/22/05) ! call fprint8('Nighttime ionization:',nlev, ! | qn2p,qnp,qo2p,qop4s,qnspe,dum,dum,dum, ! | 'qn2p','qnp','qo2p','qop4s','qnspe',' ',' ',' ') ! ! Column number density for O3: call colum(xmlw,xxn,col) ! ! Solar zenith angle: zi = 60. ! ! Main levels loop: do k=1,nlev alt = zpht(k) xp = (re+alt*1.e+5)/sht(k) y=0.5*xp*(cos(zi*dtor))**2 tt=sqrt(y) if (tt > 8.) then yerf=0.56498823/(0.06651874+tt) else yerf=(1.0606963+0.55643831*tt)/(1.0619896+1.7245609*tt+tt*tt) endif clno3(k) = col(k) chaps=sqrt(0.5*pi*xp)*yerf schap(k)=chaps ! ! Slant column number densities: txno(k) = cno(k) *chaps txno2(k) = cno2(k) *chaps txnn2(k) = cnn2(k) *chaps txno3(k) = clno3(k)*chaps ! ! SRC photodissociation (yee): scxno2 = txno2(k) call srcj(scxno2,xj) call srceps(scxno2,eps) yeesrcj(k) = xj*cgm yeesrch(k) = eps*xno2(k)/rho(k)*cgm ! ! SRB photodissociation: sfac=1.+0.11*(f107-65.)/165. dh2osrb(k)=1.2e-6 *exp(-1.e-7*txno2(k)**0.35)*cgm*sfac dco2srb(k)=1.90e-9*exp(-1.e-7*txno2(k)**0.35)*cgm*sfac ! ! Kockarts LY-A parameterization: ! #ifdef EXPO sjmlya(k) = (0.68431 *expo(-8.22114e-21*txno2(k),0)+ | 0.229841 *expo(-1.77556e-20*txno2(k),0)+ | 0.0865412*expo(-8.22112e-21*txno2(k),0))* | sflux(12)*cgm sjo2lya(k) = (6.0073e-21 *expo(-8.2166e-21 *txno2(k),0)+ | 4.28569e-21*expo(-1.63296e-20*txno2(k),0)+ | 1.28059e-20*expo(-4.85121e-17*txno2(k),0))* | sflux(12)*cgm qnop(k) = qnop(k)+2.02e-18*sjmlya(k)*xnno(k) sjnmlya(k) = (0.68431 *expo(-8.22114e-21*cno2(k),0)+ | 0.229841 *expo(-1.77556e-20*cno2(k),0)+ | 0.0865412*expo(-8.22112e-21*cno2(k),0)) #else sjmlya(k) = (0.68431 *exp(-8.22114e-21*txno2(k))+ | 0.229841 *exp(-1.77556e-20*txno2(k))+ | 0.0865412*exp(-8.22112e-21*txno2(k)))* | sflux(12)*cgm sjo2lya(k) = (6.0073e-21 *exp(-8.2166e-21 *txno2(k))+ | 4.28569e-21*exp(-1.63296e-20*txno2(k))+ | 1.28059e-20*exp(-4.85121e-17*txno2(k)))* | sflux(12)*cgm qnop(k) = qnop(k)+2.02e-18*sjmlya(k)*xnno(k) sjnmlya(k) = (0.68431 *exp(-8.22114e-21*cno2(k))+ | 0.229841 *exp(-1.77556e-20*cno2(k))+ | 0.0865412*exp(-8.22112e-21*cno2(k))) #endif qnop(k) = 2.02e-18*sflux(12)/100.*sjnmlya(k)*xnno(k) ! dinolya(k) = 2.02e-18*sjmlya(k) dh2olya(k) = 1.4e-17 *sjmlya(k) do3lya(k) = 2.27e-17*sjmlya(k) do2lya(k) = sjo2lya(k) dch4lya(k) = 1.85e-17*sjmlya(k) dco2lya(k) = 8.14e-20*sjmlya(k) qnlya(k)=(do2lya(k)*xno2(k)+dh2olya(k)*xnh2o(k)+dch4lya(k)* | xnch4(k)+dco2lya(k)*xnco2(k))*1.602e-12*(12400./1216.)/rho(k) ! ! Solar x-ray production: do l=1,3 ! was DO 20 xrflx(l)=exray*exflxr(l)*xlam(l)/(12400.*erg) tau=(sigxo2(l)*cno2(k)+sigxn2(l)*cnn2(k)+sigxo(l)*cno(k))* | 1.e-19 ctau=0.5*(exp(-tau)-tau*e1(tau)) ! func e1 is in solvers.F thng1=xno2(k)*sigxo2(l)*1.e-19*ctau*xrflx(l)*effx(l) thng2=xnn2(k)*sigxn2(l)*1.e-19*ctau*xrflx(l)*effx(l) thng3=xno(k) *sigxo(l) *1.e-19*ctau*xrflx(l)*effx(l) qo2p(k)=qo2p(k)+thng1*0.67 qn2p(k)=qn2p(k)+thng2*0.64 ! qn2p def 2 ! write(6,"('qn2p def 2: l=',i2,' k=',i4,' sigxn2(l)=',e12.4, ! | ' ctau=',e12.4,' xrflx(l)=',e12.4,' effx(l)=',e12.4)") ! | l,k,sigxn2(l),ctau,xrflx(l),effx(l) qnp(k) =qnp(k) +thng2*0.36 qop4s(k)=qop4s(k)+(thng3+0.33*thng1) qixray(k)=qixray(k)+thng1+thng2+thng3 qnspe(k)=qnspe(k)+1.57*thng2 qn(k)=qn(k)+qixray(k)*35.*eff*1.602e-12/rho(k) enddo ! write(6,"('qn2p def 2: k=',i4,' qn2p(k)=',e12.4)") k,qn2p(k) ! ! Solar EUV secondary electrons: taur=sigabs(1,49)*txno(k)+sigabs(3,49)*txnn2(k)+ | sigabs(2,49)*txno2(k) taup=sigabs(1,20)*txno(k)+sigabs(3,20)*txnn2(k)+ | sigabs(2,20)*txno2(k) taur1=1.3*taur taur2=2.0*taur taur3=2.5*taur if (taur > 9.) taur = 9. if (taup > 9.) taup = 9. if (taur1 > 9.) taur1 = 9. if (taur2 > 9.) taur2 = 9. if (taur3 > 9.) taur3 = 9. ! ! e1 is a function in solvers.F ctaur = 0.5*(exp(-taur) -taur *e1(taur)) ctaur1 = 0.5*(exp(-taur1)-taur1*e1(taur1)) ctaur2 = 0.5*(exp(-taur2)-taur2*e1(taur2)) ctaur3 = 0.5*(exp(-taur3)-taur3*e1(taur3)) ctaup = 0.5*(exp(-taup) -taup *e1(taup)) rsp=2.4*ctaur/(ctaur+2.*(ctaur1+ctaur2+ctaur3)) rsq=1.5*ctaur/(ctaur+2.*(ctaur1+ctaur2+ctaur3)+taup/taur*ctaup) cspi=1.0+rsp cspj=1.0+rsq ! ! Loop over solar euv and uv wavelength: do l=1,lmax tau=sigabs(1,l)*cno(k)+sigabs(2,l)*cno2(k)+sigabs(3,l)*cnn2(k) #ifdef EXPO ctau=0.5*(expo(-tau,0)-tau*e1(tau)) #else ctau=0.5*(exp(-tau)-tau*e1(tau)) #endif thng1=xno2(k)*sigion(2 ,l)*ctau*sflux(l) thng2=xnn2(k)*sigion(3 ,l)*ctau*sflux(l) thng3=xno(k) *sigion(1 ,l)*ctau*sflux(l) thng4=xn4s(k)*sigion(12,l)*ctau*sflux(l) thng5=xno(k) *sigion(8 ,l)*ctau*sflux(l) thng6=xno(k) *sigion(7 ,l)*ctau*sflux(l) thng7=xno(k) *sigion(6 ,l)*ctau*sflux(l) qo2p(k)=qo2p(k)+thng1*(1.-bro2(l)) qn2p(k)=qn2p(k)+thng2*(1.-brn2(l)) ! qn2p def 3 qnp(k) =qnp(k) +thng2*brn2(l)+thng4 qop2p(k)=qop2p(k)+thng5 qop2d(k)=qop2d(k)+thng6 qop4s(k)=qop4s(k)+thng7+thng1*bro2(l) diop(k) =diop(k) +sigion(1,l)*ctau*sflux(l) dio2p(k)=dio2p(k)+sigion(2,l)*ctau*sflux(l) din2p(k)=din2p(k)+sigion(3,l)*ctau*sflux(l) dinp(k) =dinp(k)+sigion(12,l)*ctau*sflux(l) diop2p(k)=diop2p(k)+sigion(8,l)*ctau*sflux(l) diop2d(k)=diop2p(k)+sigion(7,l)*ctau*sflux(l) diop4s(k)=diop4s(k)+sigion(6,l)*ctau*sflux(l) pnoeuv(k)=pnoeuv(k)+xnn2(k)*(sigabs(3,l)-sigion(3,l))*ctau* | sflux(l) dnoeuv(k)=dnoeuv(k)+sigabs(9,l)*ctau*sflux(l) qnspe(k)=qnspe(k)+1.57*thng2*rsp if (l >= 13) | do2euv(k)=do2euv(k)+(sigabs(2,l)-sigion(2,l))*ctau*sflux(l) ! ! EUV dissociation ! if (l >= 16) then dh2oeuv(k)=dh2oeuv(k)+sigabs(8,l)*ctau*sflux(l) dco2euv(k)=dco2euv(k)+sigabs(7,l)*ctau*sflux(l) dcoeuv(k) =dcoeuv(k) +1.7e-17*ctau*sflux(l) do3euv(k) =do3euv(k) +sigabs(6,l)*ctau*sflux(l) wl1=wave1(l) sgch4=1.e-30 if (wl1 < 300.) sgch4=1.e-17 if (wl1 >= 300. .and. wl1 < 900.) | sgch4=1.e-17+(4.e-17/600.)*(wl1-300.) if(wl1 >= 900. .and. wl1 < 1050.) | sgch4=5.e-17-(2.e-17/130.)*(wl1-900.) dch4euv(k)=dch4euv(k)+sgch4*ctau*sflux(l) endif ! l >= 16 ! ! Schumann-runge continuum ! if (l <= 11) then tau=sigabs(2,l)*cno2(k) #ifdef EXPO thng=sflux(l)*sigabs(2,l)*cgm*expo(-sigabs(2,l)*txno2(k),0) #else thng=sflux(l)*sigabs(2,l)*cgm*exp(-sigabs(2,l)*txno2(k)) #endif thnh=thng/sigabs(2,l) do2src(k)=yeesrcj(k) alamav=(wave1(l)+wave2(l))/2. essr= 12400./alamav-7.12 if (essr.lt.0.) essr=0. qsrc(k)=qsrc(k)+thng*xno2(k)*1.602e-12*essr/rho(k) dh2osrc(k)=dh2osrc(k)+thnh*sigabs(8,l) dco2src(k)=dco2src(k)+thnh*sigabs(7,l) sr63(k)=sr63(k)+xno2(k)*thng*a6300/ | (a1d*(1.+quench*xnn2(k)/a1d)) endif ! ! End wavelength loop: enddo ! l=1,lmax ! write(6,"('qn2p def 3: k=',i4,' qn2p(k)=',e12.4)") k,qn2p(k) ! qops(k)=qop2p(k)+qop2d(k)+qop4s(k) qn2p(k)=qn2p(k)*(1.+rsp) ! qn2p def 4 ! write(6,"('qn2p def 4: k=',i4,' qn2p(k)=',e12.4)") k,qn2p(k) qn2ps(k)=qn2p(k) qnp(k)=qnp(k)*(1.+rsp) qo2p(k)=qo2p(k)*(1.+rsq) qop2p(k)=qop2p(k)+qops(k)*rsp*0.20 qop2d(k)=qop2d(k)+qops(k)*rsp*0.24 qop4s(k)=qop4s(k)+qops(k)*rsp*0.56 qnop(k)=qnop(k)+dnoeuv(k)*xnno(k) qtin(k)=qo2p(k)+qn2p(k)+qnp(k)+qnop(k)+qop2p(k)+ | qop2d(k)+qop4s(k) qn2ps(k)=qn2p(k) qops(k)=qop2p(k)+qop2d(k)+qop4s(k) ! ! Photoelectron heating of neutral gas qnpe(k)=qtin(k) qn(k)=qn(k)+qnpe(k)*35.*effpe*1.602e-12/rho(k) ! ! qlyb = photo dissociation of o2 by solar lyman-beta, as 5th ! source of greenline e5577 emission: ! taubt = sigabs(1,17)*txno(k)+sigabs(2,17)*txno2(k)+ | sigabs(3,17)*txnn2(k) qlyb(k) = 0.1*sflux(17)*xno2(k)*sigabs(2,17)*exp(-taubt)*cgm ! ! End levels loop: enddo ! k=1,nlev ! ! 11/17/04 btf: Only diffs are due to machine range of 8-byte reals. ! IBM Dave provides 1.e-320, whereas lightning linux cluster ! flushes to zero. This occurs in yeesrj, qnop, sjnmlya, qnop. ! ! if (istep==nstep) then ! write(char80,"('solheat 1, step=',i5)") istep ! call fprint8(char80,nlev, ! | txno , txno2 , txnn2 , txno3 , yeesrcj , yeesrch,dum,dum, ! | 'txno','txno2','txnn2','txno3','yeesrcj','yeesrch',' ',' ') ! write(char80,"('solheat 2, step=',i5)") istep ! call fprint8(char80,nlev, ! | dh2osrb , dco2srb , sjmlya , sjo2lya , qnop , sjnmlya , qnop , ! | dum, ! | 'dh2osrb','dco2srb','sjmlya','sjo2lya','qnop','sjnmlya','qnop', ! | ' ') ! write(char80,"('solheat 3, step=',i5)") istep ! call fprint8(char80,nlev, ! | dinolya , dh2olya , do3lya , do2lya , dch4lya , dco2lya , ! | qnlya, dum, ! | 'dinolya','dh2olya','do3lya','do2lya','dch4lya','dco2lya', ! | 'qnlya',' ') ! write(char80,"('solheat 4, step=',i5)") istep ! call fprint8(char80,nlev, ! | qo2p , qn2p , qnp , qop4s , qixray , qnspe , qn , dum, ! | 'qo2p','qn2p','qnp','qop4s','qixray','qnspe','qn',' ') ! write(char80,"('solheat 5, step=',i5)") istep ! call fprint8(char80,nlev, ! | qop2p , qop2d , diop , dio2p , din2p , dinp , dum,dum, ! | 'qop2p','qop2d','diop','dio2p','din2p','dinp',' ',' ') ! write(char80,"('solheat 6, step=',i5)") istep ! call fprint8(char80,nlev, ! | diop2p , diop2d , diop4s , pnoeuv , dnoeuv , qnspe , dum,dum, ! | 'diop2p','diop2d','diop4s','pnoeuv','dnoeuv','qnspe',' ',' ') ! write(char80,"('solheat 7, step=',i5)") istep ! call fprint8(char80,nlev, ! | qops , qtin , qn2ps , qlyb , qnpe , dh2osrc , dco2src , dum, ! | 'qops','qtin','qn2ps','qlyb','qnpe','dh2osrc','dco2src', ' ') ! endif ! last timestep ! ! Solar uv UARS from LY-A to 800 NM: call uars_solaruv ! do k=1,nlev ! was DO 222 do3t(k)=do3har(k)+do3hug(k)+do3chap(k)+do3euv(k)+do3lya(k) | +do3hz(k)+do3srb(k)+do3src(k) do2t(k)=yeesrcj(k)+do2srb(k)+do2lya(k)+do2hz(k)+do2euv(k) qn(k)=qn(k)+yeesrch(k)+qsrb(k)+qnlya(k)+qhar(k)+qhug(k) | +qchap(k)+qherz(k) dch4t(k)=dch4lya(k)+dch4euv(k)+dch4b(k) dco2t(k)=dco2lya(k)+dco2srb(k)+dco2src(k)+dco2euv(k) dh2ot(k)=dh2olya(k)+dh2osrc(k)+dh2osrb(k)+dh2oeuv(k) enddo ! ! Cosmic ray ion production and O2(1SG) ! call cosray(45.) do k=1,nlev ! was DO 16 ! ! O2(1SG) xno2s(k)=xno21d(k) qo2s1=0.549e-09*exp(-2.406e-20*txno2(k))*cgm qo2s2=2.614e-09*exp(-8.508e-20*txno2(k))*cgm qo2s(k)=xno2s(k)*(qo2s1+qo2s2) ! ! Cosmic ray ion production qn2p(k) =qn2p (k)+0.585*qcr(k) qnp(k) =qnp (k)+0.185*qcr(k) qo2p(k) =qo2p (k)+0.154*qcr(k)+qo2s(k) qop4s(k)=qop4s(k)+0.076*qcr(k) qnspe(k)=qnspe(k)+0.585*qcr(k)*1.57 enddo ! k=1,nlev ! was DO 16 ! call fprint8('solheat',nlev, ! | qo2s , qn2p , qnp , qo2p , qop4s , qnspe,dum,dum, ! | 'qo2s','qn2p','qnp','qo2p','qop4s','qnspe',' ',' ') ! ! Aurora production: aureff=0.1 alp=0.75 eaur=0.025 fluxar=eaur/(2.*alp*1.602e-09) ! ! Qaurora (aurora.F) defines qia(1:5,nlev) (module flds_ionzrt) if (iaurora > 0) then call qaurora(fluxar,alp) ! do k=1,nlev qn2p(k) =qn2p (k)+qia(1,k) qo2p(k) =qo2p (k)+qia(2,k) qop4s(k)=qop4s(k)+qia(3,k) qnp(k) =qnp (k)+qia(5,k) qtin(k)=qop4s(k)+qo2p(k)+qn2p(k)+qnp(k)+qnop(k)+qop2p(k)+ | qop2d(k) qiaur(k)=qia(1,k)+qia(2,k)+qia(3,k)+qia(5,k) qian(k)=35.*qti(k)*1.602e-12*aureff/rho(k) qaurh(k)=qian(k) qn(k)=qn(k)+qian(k) qnspe(k)=qnspe(k)+1.57*qia(1,k) enddo ! k=1,nlev endif ! iaurora ! ! Solar proton enhancement over polar cap (reid) ! Sub protonp (aurora.F) does not appear to return anything ! that is used later (qpe and qpi are calculated locally) ! call protonp ! ! isolproton is in namelist (input.F), with default == 0. if (isolproton > 0) then do k=1,nlev ! was DO 31 zprs(k)=zp(k) als(k)=zpht(k) tu(k)=tn(k) tdena(k)=(xno(k)*16.+xno2(k)*32.+xnn2(k)*28.)*1.66e-24 ab(k)=xno(k)+xno2(k)+xnn2(k) enddo ! k=1,nlev ! call fprint6('solheat before proton:',nlev, ! | zprs , als , tu , tdena , ab , dum, ! | 'zprs','als','tu','tdena','ab', ' ') ! ! Sub proton (aurora.F) calculates total ion production in qpro ! (module flds_ionzrt, formerly ionzrt.h) ! 11/19/04 btf: sub proton is unfinished.. ! call proton(ut) ! do k=1,nlev ! was DO 32 qo2p(k)=qo2p(k)+qpro(k) qtin(k)=qtin(k)+qpro(k) qproh(k)=35.*qpro(k)*1.602e-12/rho(k) qian(k)=qian(k)+qproh(k) qn(k)=qn(k)+qproh(k) qnspe(k)=qnspe(k)+0.585*qpro(k)*1.57 enddo ! k=1,nlev endif ! isolproton ! do k=1,nlev ! was DO 337 qop(k)=qop4s(k)+qop2d(k)+qop2p(k) enddo ! k=1,nlev end subroutine solheat !----------------------------------------------------------------------- subroutine seteuv ! ! Init module data (formerly block data blk_euv). ! ! lambdas: wleuv1 = | (/1000.00, 1031.91, 1025.72, 950.00, 977.02, 900.00, | 850.00, 800.00, 750.00, 789.36, 770.41, 765.15, | 700.00, 703.36, 650.00, 600.00, 629.73, 609.76, | 550.00, 584.33, 554.31, 500.00, 450.00, 465.22, | 400.00, 350.00, 368.07, 300.00, 303.78, 303.31, | 250.00, 284.15, 256.30, 200.00, 150.00, 100.00, | 50.00/) wleuv2 = | (/1050.00, 1031.91, 1025.72, 1000.00, 977.02, 950.00, | 900.00, 850.00, 800.00, 789.36, 770.41, 765.15, | 750.00, 703.36, 700.00, 650.00, 629.73, 609.76, | 600.00, 584.33, 554.31, 550.00, 500.00, 465.22, | 450.00, 400.00, 368.07, 350.00, 303.78, 303.31, | 300.00, 284.15, 256.30, 250.00, 200.00, 150.00, | 100.00/) c c sigao sigao = | (/0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, 1.32e-18, | 4.55e-18, 3.50e-18, 5.09e-18, 3.75e-18, 3.89e-18, 4.00e-18, | 1.07e-17, 1.15e-17, 1.72e-17, 1.34e-17, 1.34e-17, 1.34e-17, | 1.30e-17, 1.31e-17, 1.26e-17, 1.21e-17, 1.21e-17, 1.19e-17, | 1.15e-17, 9.69e-18, 9.84e-18, 8.69e-18, 7.70e-18, 7.68e-18, | 6.46e-18, 7.08e-18, 6.05e-18, 5.20e-18, 3.73e-18, 1.84e-18, | 7.30e-19/) ! ! sigao2 sigao2 = | (/1.35e-18, 1.00e-18, 1.63e-18, 2.11e-17, 1.87e-17, 1.28e-17, | 8.56e-18, 1.66e-17, 2.21e-17, 2.67e-17, 1.89e-17, 2.08e-17, | 2.85e-17, 2.74e-17, 2.19e-17, 2.60e-17, 3.21e-17, 2.81e-17, | 2.66e-17, 2.28e-17, 2.60e-17, 2.46e-17, 2.31e-17, 2.19e-17, | 2.03e-17, 1.81e-17, 1.83e-17, 1.74e-17, 1.68e-17, 1.68e-17, | 1.44e-17, 1.58e-17, 1.34e-17, 1.09e-17, 7.51e-18, 3.81e-18, | 1.32e-18/) ! ! sigan2 sigan2 = | (/0.00e+00, 0.00e+00, 0.00e+00, 5.10e-17, 2.24e-18, 9.68e-18, | 2.02e-17, 1.70e-17, 3.36e-17, 1.65e-17, 1.42e-17, 1.20e-16, | 2.47e-17, 2.65e-17, 3.18e-17, 2.33e-17, 2.34e-17, 2.28e-17, | 2.28e-17, 2.24e-17, 2.41e-17, 2.45e-17, 2.35e-17, 2.32e-17, | 2.17e-17, 1.64e-17, 1.69e-17, 1.39e-17, 1.17e-17, 1.17e-17, | 1.05e-17, 1.09e-17, 1.02e-17, 8.39e-18, 4.96e-18, 2.26e-18, | 7.20e-19/) ! ! sigio sigio = | (/0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, 1.32e-18, | 4.55e-18, 3.50e-18, 5.09e-18, 3.75e-18, 3.89e-18, 4.00e-18, | 1.07e-17, 1.15e-17, 1.72e-17, 1.34e-17, 1.34e-17, 1.34e-17, | 1.30e-17, 1.31e-17, 1.26e-17, 1.21e-17, 1.21e-17, 1.19e-17, | 1.15e-17, 9.69e-18, 9.84e-18, 8.69e-18, 7.70e-18, 7.68e-18, | 6.46e-18, 7.08e-18, 6.05e-18, 5.20e-18, 3.73e-18, 1.84e-18, | 7.30e-19/) ! ! sigio2 sigio2 = | (/2.59e-19, 0.00e+00, 1.05e-18, 1.39e-17, 1.55e-17, 9.37e-18, | 5.49e-18, 6.41e-18, 1.06e-17, 1.02e-17, 8.47e-18, 1.17e-17, | 2.38e-17, 2.38e-17, 2.13e-17, 2.49e-17, 3.11e-17, 2.64e-17, | 2.66e-17, 2.28e-17, 2.60e-17, 2.46e-17, 2.31e-17, 2.19e-17, | 2.03e-17, 1.81e-17, 1.83e-17, 1.74e-17, 1.68e-17, 1.68e-17, | 1.44e-17, 1.58e-17, 1.34e-17, 1.09e-17, 7.51e-18, 3.81e-18, | 1.32e-18/) ! ! sigin2 sigin2 = | (/0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, | 0.00e+00, 0.00e+00, 1.43e-17, 8.86e-18, 8.50e-18, 6.58e-17, | 1.51e-17, 2.55e-17, 2.92e-17, 2.33e-17, 2.34e-17, 2.28e-17, | 2.28e-17, 2.24e-17, 2.41e-17, 2.45e-17, 2.35e-17, 2.32e-17, | 2.17e-17, 1.64e-17, 1.69e-17, 1.39e-17, 1.17e-17, 1.17e-17, | 1.05e-17, 1.09e-17, 1.02e-17, 8.39e-18, 4.96e-18, 2.26e-18, | 7.20e-19/) ! ! sigin sigin = | (/0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, | 0.21e-18,10.29e-18,11.71e-18,10.96e-18,11.24e-18,11.32e-18, | 12.10e-18,13.26e-18,12.42e-18,11.95e-18,11.21e-18,11.80e-18, | 11.76e-18,11.78e-18,11.77e-18,11.50e-18,11.02e-18,10.58e-18, | 9.56e-18, 8.15e-18, 8.30e-18, 7.30e-18, 6.41e-18, 6.40e-18, | 5.24e-18, 5.73e-18, 4.87e-18, 3.95e-18, 2.49e-18, 0.99e-18, | 0.33e-18/) ! ! O 4S brop4s = | (/0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, 1.00e+00, | 1.00e+00, 1.00e+00, 1.00e+00, 1.00e+00, 1.00e+00, 1.00e+00, | 6.30e-01, 4.30e-01, 3.00e-01, 3.20e-01, 2.90e-01, 2.70e-01, | 2.80e-01, 3.00e-01, 2.90e-01, 2.80e-01, 2.80e-01, 2.80e-01, | 2.70e-01, 2.60e-01, 2.60e-01, 2.60e-01, 2.50e-01, 2.50e-01, | 2.50e-01, 2.50e-01, 2.50e-01, 2.60e-01, 2.70e-01, 2.90e-01, | 3.00e-01/) ! ! O 2D brop2d = | (/0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, | 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, | 3.70e-01, 5.70e-01, 6.60e-01, 4.60e-01, 4.70e-01, 4.60e-01, | 4.50e-01, 4.60e-01, 4.50e-01, 4.50e-01, 4.50e-01, 4.50e-01, | 4.30e-01, 4.00e-01, 4.00e-01, 4.00e-01, 3.70e-01, 3.70e-01, | 3.60e-01, 3.70e-01, 3.50e-01, 3.50e-01, 3.30e-01, 3.20e-01, | 3.20e-01/) ! ! O 2P brop2p = | (/0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, | 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, | 0.00e+00, 0.00e+00, 4.00e-02, 2.20e-01, 2.40e-01, 2.70e-01, | 2.70e-01, 2.40e-01, 2.60e-01, 2.70e-01, 2.70e-01, 2.70e-01, | 2.60e-01, 2.50e-01, 2.60e-01, 2.50e-01, 2.50e-01, 2.50e-01, | 2.30e-01, 2.30e-01, 2.30e-01, 2.20e-01, 2.20e-01, 2.10e-01, | 2.10e-01/) ! ! N2 -> N+ brn2np = | (/0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, | 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, | 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, | 0.00e+00, 0.00e+00, 0.00e+00, 1.30e-03, 3.00e-02, 4.60e-02, | 4.50e-02, 1.05e-01, 9.00e-02, 1.63e-01, 2.13e-01, 2.13e-01, | 3.00e-01, 2.57e-01, 3.35e-01, 3.77e-01, 3.64e-01, 3.46e-01, | 3.85e-01/) ! ! O2 -> O+ bro2op = | (/0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, | 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, | 0.00e+00, 0.00e+00, 8.91e-03, 3.86e-02, 3.31e-02, 7.01e-02, | 1.44e-01, 1.71e-01, 1.71e-01, 2.08e-01, 2.31e-01, 2.25e-01, | 2.37e-01, 2.59e-01, 2.49e-01, 2.99e-01, 3.53e-01, 3.53e-01, | 3.71e-01, 3.73e-01, 3.66e-01, 3.93e-01, 4.48e-01, 3.84e-01, | 0.00e+00/) end subroutine seteuv !----------------------------------------------------------------------- subroutine ssflux(iscale,f107,f107a,hlybr,fexvir,hlya, | heiew,xuvfac,wave1,wave2,sflux,itimet) ! ! Args: integer,intent(in) :: iscale real,intent(in) :: hlybr,fexvir,hlya,heiew,xuvfac,f107,f107a real,dimension(lmax),intent(out) :: wave1,wave2,sflux integer,intent(inout) :: itimet ! ! Local: integer :: l,m real,dimension(lmax) :: wavel,waves,rflux,xflux,scale1,scale2, | tchr0,tchr1,tchr2,tcor0,tcor1,tcor2,war1,war2 real :: b1(3),b2(3),frat,r1,r2,hlymod,heimod,xuvf,dum real,dimension(lmax) :: | ssigao ,ssigao2 ,ssigan2 ,ssigio ,ssigin2 ,ssigio2 ,ssigaco , | ssigico ,ssigaco2,ssigico2,ssigio4s,ssigio2d,ssigio2p,ssigio4p , | ssigio2q,ssigihe ,ssigin ,ssigih ,ssigao3 ,ssigbco2,ssigah2o , | ssigano ! ! regression coefficients which reduce to solar min. spectrum: b1 = (/1.0, 0.0138 , 0.005/) b2 = (/1.0, 0.59425, 0.3811/) ! wavel = (/ 1750.00, 1700.00, 1650.00, 1600.00, 1550.00, 1500.00, | 1450.00, 1400.00, 1350.00, 1300.00, 1250.00, 1215.67, | 1200.00, 1150.00, 1100.00, 1050.00, 1031.91, 1025.72, | 1000.00, 977.02, 950.00, 900.00, 850.00, 800.00, | 789.36, 770.41, 765.15, 750.00, 703.31, 700.00, | 650.00, 629.73, 609.76, 600.00, 584.33, 554.37, | 550.00, 500.00, 465.22, 450.00, 400.00, 368.07, | 350.00, 303.78, 303.31, 300.00, 284.15, 256.30, | 250.00, 200.00, 150.00, 100.00, 50.00, 32.00, | 23.00, 16.00, 8.00, 4.00, 2.00/) waves = (/ 1700.00, 1650.00, 1600.00, 1550.00, 1500.00, 1450.00, | 1400.00, 1350.00, 1300.00, 1250.00, 1200.00, 1215.67, | 1150.00, 1100.00, 1050.00, 1000.00, 1031.91, 1025.72, | 950.00, 977.02, 900.00, 850.00, 800.00, 750.00, | 789.36, 770.41, 765.15, 700.00, 703.31, 650.00, | 600.00, 629.73, 609.76, 550.00, 584.33, 554.37, | 500.00, 450.00, 465.22, 400.00, 350.00, 368.07, | 300.00, 303.78, 303.31, 250.00, 284.15, 256.30, | 200.00, 150.00, 100.00, 50.00, 32.00, 23.00, | 16.00, 8.00, 4.00, 2.00, 1.00/) rflux = (/ 322.00, 168.00, 95.00, 62.00, 44.00, 25.00, | 16.90, 11.80, 19.50, 4.10, 11.10, 249.00, | 2.78, 0.70, 3.07, 3.64, 3.18, 4.38, | 1.78, 5.96, 4.22, 4.43, 1.93, 0.87, | 0.79, 0.24, 0.20, 0.17, 0.39, 0.22, | 0.17, 1.50, 0.45, 0.48, 1.58, 0.80, | 0.51, 0.31, 0.18, 0.39, 0.21, 0.74, | 0.87, 6.00, 0.24, 0.84, 0.10, 0.27, | 0.92, 1.84, 0.13, 0.38, 0.0215, 0.0067, | 1.E-3, 2.E-3, 1.E-5, 5.E-8, 1.E-10/) xflux = (/ 354.00, 191.00, 110.00, 76.00, 55.00, 28.00, | 19.60, 14.30, 25.30, 5.00, 17.20, 401.00, | 6.26, 1.51, 6.11, 8.66, 9.04, 13.12, | 4.42, 13.18, 12.03, 13.29, 5.01, 2.18, | 1.59, 0.67, 0.43, 0.43, 0.72, 0.46, | 0.48, 3.02, 1.46, 1.02, 4.86, 1.59, | 1.57, 1.67, 0.36, 0.99, 2.20, 1.39, | 5.63, 11.28, 2.50, 4.14, 3.16, 0.59, | 3.70, 4.85, 0.34, 1.15, 0.18, 0.08, | 2.5E-2, 5.E-2, 8.E-4, 3.E-5, 5.E-7/) scale1 = (/ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, | 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, | 1692.09, 405.95, 1516.20, 2731.70, 3314.57, 4375.00, | 1316.91, 3621.91, 3908.56, 4432.54, 1541.21, 531.73, | 364.83, 0.00, 116.00, 129.41, 162.48, 94.07, | 41.29, 709.50, 0.00, 268.47, 1561.05, 367.64, | 290.06, 184.36, 0.00, 86.15, 7.50, 0.00, | 0.00, 2220.00, 0.00, 61.00, 0.00, 86.95, | 206.00, 135.89, 60.35, 157.12, 7.06, 0.75, | 0.00, 0.00, 0.00, 0.00, 0.00/) scale2 = (/ 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, | 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, | 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, | 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, | 0.00, 5.34, 0.00, 0.00, 0.00, 0.54, | 3.30, 0.00, 12.60, 0.00, 0.00, 0.00, | 5.34, 11.63, 2.28, 5.56, 24.93, 8.16, | 60.69, 0.00, 28.20, 45.90, 40.80, 1.27, | 35.47, 42.80, 1.12, 6.19, 1.26, 0.69, | 0.23, 0.46, 7.6E-3, 2.9E-4, 4.8E-6/) tchr0 = (/ | 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, | 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, | 0.0, 0.0, 0.0,-4.290E+00,-5.709E+00,-8.493E+00, |-1.161E+00,-3.429E+00,-5.464E+00,-6.502E+00,-1.912E+00,-4.034E-01, |-1.448E-01, 0.000E+00,-9.702E-02,-6.591E-02,-2.338E-02,-1.273E-01, |-2.406E-01,-3.351E-01, 0.000E+00,-1.465E+00,-2.405E+00,-7.975E-02, |-4.197E-01,-1.971E-01, 0.000E+00,-5.895E-02,-5.815E-03, 0.000E+00, | 0.000E+00, 2.138E-01, 0.000E+00,-7.713E-02, 0.000E+00,-3.035E-02, |-2.039E-01,-1.749E-01,-1.041E-01,-2.638E-01,-1.094E-02, 0.000E+00, | 0.0, 0.0, 0.0, 0.0, 0.0/) tchr1 = (/ | 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, | 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, | 0.0, 0.0, 0.0,-3.023E-13,-3.745E-13,-5.385E-13, |-1.211E-13,-3.868E-13,-3.646E-13,-4.125E-13,-1.527E-13,-4.753E-14, |-3.411E-14, 0.000E+00,-1.190E-14,-1.034E-14,-1.343E-14,-1.539E-14, |-5.174E-14,-6.934E-14, 0.000E+00,-1.215E-13,-1.537E-13,-2.024E-14, |-4.596E-14,-1.562E-14, 0.000E+00,-1.221E-14,-1.123E-15, 0.000E+00, | 0.000E+00,-2.263E-13, 0.000E+00,-1.508E-14, 0.000E+00,-1.744E-14, |-2.100E-14,-1.805E-14,-8.224E-15,-1.919E-14,-7.944E-16, 0.000E+00, | 0.0, 0.0, 0.0, 0.0, 0.0/) tchr2 = (/ | 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, | 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, | 0.0, 0.0, 0.0, 3.275E-11, 4.057E-11, 6.160E-11, | 1.312E-11, 4.189E-11, 4.167E-11, 4.716E-11, 1.654E-11, 5.150E-12, | 3.901E-12, 0.000E+00, 1.289E-12, 1.120E-12, 1.455E-12, 1.667E-12, | 5.604E-12, 7.931E-12, 0.000E+00, 1.317E-11, 1.757E-11, 2.194E-12, | 4.978E-12, 1.693E-12, 0.000E+00, 1.324E-12, 1.285E-13, 0.000E+00, | 0.000E+00, 2.586E-11, 0.000E+00, 1.724E-12, 0.000E+00, 1.889E-12, | 2.400E-12, 2.063E-12, 8.911E-13, 2.193E-12, 9.090E-14, 0.000E+00, | 0.0, 0.0, 0.0, 0.0, 0.0/) tcor0 = (/ | 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, | 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, | 0.0, 0.0, 0.0, 0.000E+00, 0.000E+00, 0.000E+00, | 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00,-6.060E-02, | 0.000E+00,-3.399E-02, 0.000E+00, 0.000E+00, 0.000E+00, 4.866E-02, |-1.762E-01, 0.000E+00,-2.412E-01, 0.000E+00, 0.000E+00, 0.000E+00, |-4.743E-01,-9.713E-01, 5.891E-02,-1.263E-01,-1.246E+00, 2.870E-01, |-4.659E+00, 0.000E+00,-1.058E+00,-3.821E+00,-1.874E+00, 0.000E+00, |-1.896E+00,-8.505E-01,-2.101E-04,-2.012E-01,-6.097E-02,-2.925E-02, |-4.875E-03, 0.0, 0.0, 0.0, 0.0/) tcor1 = (/ | 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, | 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, | 0.0, 0.0, 0.0, 0.000E+00, 0.000E+00, 0.000E+00, | 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 2.877E-03, | 0.000E+00, 1.760E-03, 0.000E+00, 0.000E+00, 0.000E+00, 3.313E-04, | 3.643E-03, 0.000E+00, 5.225E-03, 0.000E+00, 0.000E+00, 0.000E+00, | 4.085E-03, 1.088E-02, 8.447E-04, 3.237E-03, 1.907E-02, 2.796E-03, | 4.460E-02, 0.000E+00, 1.007E-02, 3.481E-02, 1.604E-02, 0.000E+00, | 2.029E-02, 2.160E-02, 6.342E-04, 3.594E-03, 5.503E-04, 2.687E-04, | 4.479E-05, 0.0, 0.0, 0.0, 0.0/) tcor2 = (/ | 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, | 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, | 0.0, 0.0, 0.0, 0.000E+00, 0.000E+00, 0.000E+00, | 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 1.846E-03, | 0.000E+00, 1.127E-03, 0.000E+00, 0.000E+00, 0.000E+00, 1.891E-04, | 2.326E-03, 0.000E+00, 2.801E-03, 0.000E+00, 0.000E+00, 0.000E+00, | 2.446E-03, 7.121E-03, 5.204E-04, 1.983E-03, 1.204E-02, 1.721E-03, | 2.911E-02, 0.000E+00, 7.177E-03, 2.272E-02, 9.436E-03, 0.000E+00, | 1.316E-02, 1.398E-02, 4.098E-04, 2.328E-03, 3.574E-04, 1.745E-04, | 2.909E-05, 0.0, 0.0, 0.0, 0.0/) war1 = (/ 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, | 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, | 0.00, 0.00, 3.80, 6.25, 4.93, 6.06, | 2.70, 7.07, 8.62, 9.60, 4.54, 2.37, | 0.82, 0.33, 0.24, 0.67, 0.28, 0.55, | 1.56, 1.11, 0.77, 1.32, 1.71, 0.44, | 1.11, 0.95, 0.39, 0.81, 2.00, 1.49, | 6.81, 5.07, 1.63, 5.62, 2.08, 0.59, | 3.89, 5.19, 0.35, 1.18, 0.099, 0.04, | 0.007, 0.00, 0.00, 0.00, 0.00/) war2 = (/ 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, | 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, | 0.00, 0.00, 20.80, 17.90, 9.30, 14.30, | 6.90, 12.00, 15.60, 18.60, 10.10, 4.30, | 12.40, 8.00, 3.60, 1.80, 0.50, 1.40, | 3.90, 2.60, 1.60, 3.40, 4.10, 0.70, | 4.30, 4.30, 3.80, 2.60, 6.08, 1.35, | 12.60, 9.78, 2.96, 10.20, 4.11, 6.68, | 6.62, 8.07, 0.47, 1.73, 0.17, 0.075, | 0.012, 0.00, 0.00, 0.00, 0.00/) ! ssigao = (/ 0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0., | 0.00, 0.00, 2.12, 4.18, 4.38, 4.23, | 4.28, 4.18, 4.18, 8.00,11.35,10.04, | 12.21,12.22,12.23,11.90,12.17,12.13, | 11.91,11.64,11.25,11.21, 9.64, 9.95, | 8.67, 7.70, 7.68, 6.61, 7.13, 6.05, | 5.30, 2.90, 1.60, 0.59, 0.16, 0.05, | 0.51, 0.07, .012, .002, .0002/) ! ssigao2 = (/ 0.50, 1.50, 3.40, 6.00,10.00,13.00, | 15.00,12.00, 2.20, 0.30, 3.00, 0.01, | 0.30, 0.10, 1.00, 1.10, 1.00, 1.60, | 16.53, 4.00,15.54, 9.85,20.87,27.09, | 26.66,25.18,21.96,29.05,25.00,26.27, | 26.02,25.80,26.10,25.04,22.00,25.59, | 24.06,21.59,20.40,19.39,18.17,18.40, | 17.19,16.80,16.80,15.10,15.70,13.20, | 10.60, 7.10, 4.00, 1.18, 0.32, 0.10, | 1.02, 0.14, .024, .004, .0004/) ! ssigan2 = (/0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0., | 36.16, 0.70,16.99,46.63,15.05,30.71, | 19.26,26.88,35.46,30.94,26.30,29.75, | 23.22,23.20,23.10,22.38,23.20,24.69, | 24.53,21.85,21.80,21.07,17.51,18.00, | 13.00,11.60,11.60,10.30,10.60, 9.70, | 8.00, 4.40, 1.90, 0.60, 0.24, 1.16, | 0.48, 0.09, .015, .003, .0003/) ! ssigio = (/0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0., | 0.00, 0.00, 2.12, 4.18, 4.38, 4.23, | 4.28, 4.18, 4.18, 8.00,11.35,10.04, | 12.21,12.22,12.23,11.90,12.17,12.13, | 11.91,11.64,11.25,11.21, 9.64, 9.95, | 8.67, 7.70, 7.68, 6.61, 7.13, 6.05, | 5.30, 2.90, 1.60, 0.59, 0.16, 0.05, | 0.51, 0.07, .012, .002, .0002/) ! ssigio2 = (/ 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, | 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, | 0.00, 0.00, 0.00, 0.27, 0.00, 1.00, | 12.22, 2.50, 9.34, 4.69, 6.12, 9.39, | 11.05, 9.69, 8.59,23.81,23.00,22.05, | 25.94,25.80,26.10,25.04,22.00,25.59, | 24.06,21.59,20.40,19.39,18.17,18.40, | 17.19,16.80,16.80,15.10,15.70,13.20, | 10.60, 7.10, 4.00, 1.18, 0.32, 0.10, | 1.02, 0.14, .024, .004, .0004/) ! ssigin2 = (/0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0., | 0.00, 0.00, 0.00, 0.00, 0.00,16.75, | 10.18,18.39,23.77,23.20,23.00,25.06, | 23.22,23.20,23.10,22.38,23.20,24.69, | 24.53,21.85,21.80,21.07,17.51,18.00, | 13.00,11.60,11.60,10.30,10.60, 9.70, | 8.00, 4.40, 1.90, 0.60, 0.24, 1.16, | 0.48, 0.09, .015, .003, .0003/) ! ssigao3 = (/0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0., | 1.25, 0.92, 0.92, 0.92, 0.92, 0.80, | 0.95, 0.95, 1.30, 1.86, 2.96, 2.96, | 3.48, 3.48, 3.64, 0.44, 0.41, 3.62, | 3.46, 3.51, 3.51, 3.20, 3.06, 3.06, | 3.20, 0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0., | 0.,0.,0.,0.,0.,0.,0.,0.,0.,0./) ! ssigbco2 = (/0.05, 0.10, 0.15, 0.30, 0.40, 0.55, | 0.50, 0.50, 0.80, 0.50, 0.20, 0.00, | 0.00,18.52,14.80,18.50,14.20,15.10, | 29.60,42.90,74.10,18.50,14.80,22.20, | 31.90,37.00,93.90,22.20,35.20,18.50, | 29.60,34.30,35.30,31.50,34.20,30.40, | 33.30,25.20,27.60,20.40,27.80,27.80, | 18.50,23.40,23.40,22.90,25.90,28.20, | 22.20,11.10,0.,0.,0.,0.,0.,0.,0.,0.,0./) ! ssigah2o = (/5.00, 5.00, 5.00, 3.00, 1.50, 0.80, | 0.80, 1.10, 5.00, 8.00, 8.00, 0.00, | 4.44, 4.44, 4.44,18.52,14.10,29.60, | 15.56, 9.63,18.52,18.52,35.19,37.04, | 48.90,35.19,33.33,25.93,14.81,25.93, | 24.07,24.44,31.10,22.20,17.04,26.67, | 24.07,29.60,26.67,26.67,24.07,23.70, | 22.20,22.20,22.20,22.20,32.60,18.52, | 14.81,11.11,0.,0.,0.,0.,0.,0.,0.,0.,0./) ! ssigano = (/ 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, | 0.00, 0.00, 1.85, 2.04, 2.04, 0.00, | 2.41, 3.70, 6.48, 7.41, 3.70, 8.15, | 16.67,20.74,24.07,16.67,12.96,12.96, | 14.44,12.96, 7.41,14.44,14.44,18.52, | 20.00,20.00,20.00,20.00,20.00,20.00, | 18.52,16.67,24.81,22.22,22.22,21.85, | 18.52,22.96,22.96,25.93,19.26,25.93, | 22.22,22.22,0.,0.,0.,0.,0.,0.,0.,0.,0./) ! ! CO ssigaco = (/0.,0.,0.,0.,0.,0., | 1.92,3.53,5.48,8.02,10.02,11.70,11.01, | 12.52,12.47,13.61,15.43,15.69,18.01,19.92,20.09,21.61,22.28, | 22.52,22.41,18.42,18.60,19.78,25.59,24.45,25.98,26.28,15.26, | 33.22,21.35,22.59,37.64,49.44,28.50,52.9,0.,0.,0., | 0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0./) ! ssigico = (/0.,0.,0.,0.,0.,0., | 1.92,3.53,5.48,8.02,10.02,11.70,11.01, | 12.52,12.47,13.61,15.43,15.69,18.01,19.92,20.09,21.44,22.31, | 21.38,21.62,16.93,16.75,17.01,17.04,16.70,17.02,12.17,9.20, | 15.44,11.38,17.13,11.70,0.,0.,0.,0.,0.,0., | 0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0./) ! ! CO2 ssigaco2 = (/0.,0.,0.,0.,0.,0., | 4.42,7.51,11.03,14.98,17.88,21.21,20.00, | 23.44,23.44,23.88,25.70,25.81,27.52,28.48,29.27,31.61,33.20, | 34.21,34.00,25.31,25.86,25.88,25.96,21.76,22.48,53.96,26.48, | 21.79,31.83,12.84,49.06,70.89,29.91,34.41,0.,0.,0., | 0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0./) ! ssigico2 = (/0.,0.,0.,0.,0.,0., | 4.42,7.51,11.03,14.98,17.88,21.21,20.00, | 23.44,23.44,23.88,25.70,25.81,27.52,28.48,29.27,31.61,33.20, | 34.21,34.00,20.16,21.27,21.14,21.72,17.71,17.02,50.39,20.00, | 17.07,21.53,10.67,19.66,0.,0.,0.,0.,0.,0., | 0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0./) ! ! O+(4S) ssigio4s = (/0.,0.,0.,0.,0.,0., | .32,1.03,1.62,1.95,2.15,2.33,2.23,2.23, | 2.45,2.61,2.81,2.77,2.99,3.15,3.28,3.39,3.50,3.58,3.46,3.67, | 3.74,3.73,4.04,4.91,4.20,4.18,4.18,4.28,4.23,4.38,4.18,2.12, | 0.,0.,0.,0.,0., | 0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0./) ! ! O+(2D) ssigio2d = (/0.,0.,0.,0.,0.,0., | .34,1.14,2.00,2.62,3.02,3.39,3.18,3.62, | 3.63,3.98,4.37,4.31,4.75,5.04,5.23,5.36,5.47,5.49,5.30,5.51, | 5.50,5.50,5.52,6.44,3.80,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0., | 0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0./) ! ! O+(2P) ssigio2p = (/0.,0.,0.,0.,0.,0., | .22,.75,1.30,1.70,1.95,2.17,2.04,2.32, | 2.32,2.52,2.74,2.70,2.93,3.06,3.13,3.15,3.16,3.10,3.02,3.05, | 2.98,2.97,.47,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0., | 0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0./) ! ! O+(4P) ssigio4p = (/0.,0.,0.,0.,0.,0., | .10,.34,.58,.73,.82,.89,.85,.91,.91,.93, | .92,.92,.55,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0., | 0.,0.,0.,0.,0.,0.,0.,0., | 0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0./) ! ! O+(2P*) ssigio2q = (/0.,0.,0.,0.,0.,0., | .03,.27,.46,.54,.56,.49,.52,.41,.41, | 0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0., | 0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0., | 0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0./) ! ! HE+ ssigihe = (/0.,0.,0.,0.,0.,0., | .21,.53,1.02,1.71,2.16,2.67,2.38, | 3.05,3.05,3.65,4.35,4.25,5.51,6.53,7.09,.72, | 0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0., | 0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0./) ! ! N+ BANKS AND KOCKARTS X1.E-17 ssigin = (/0.,0.,0.,0.,0.,0., | 0.1,0.2,0.25,0.35,0.4,0.5,0.5,0.6,0.6,0.65,0.8, | 0.7,1.,1.,1.,1.1,1.15,1.2,1.1,1.2,1.2,1.2,1.2,1.2,1.1,1.1,1.1, | 1.0,1.0,1.0,1.0,1.0,0.,0.,0.,0.,0., | 0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0./) ! ! H+ BANKS AND KOCKARTS X1.E-18 ssigih = (/0.,0.,0.,0.,0.,0., | 0.05,0.02,0.05,0.12,0.16,0.20,0.23,0.27, | 0.27,0.36,0.44,0.53,0.8,0.9,1.0,1.4,1.6,1.8,1.8,2.,2.2,2.3,2.8, | 3.1,3.5,3.8,4.,4.,4.1,4.8,5.8,6.12,0.,0.,0.,0.,0., | 0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0./) ! ! Linear Interpolation between SC#21REFW and F79050: ! frat = (f107-68.) / (243.-68.) do l=1,lmax sflux(l) = rflux(l) + (xflux(l)-rflux(l)) * frat enddo ! ! Hinteregger contrast ratio method: ! if (iscale == 0) then if (hlybr .gt. 0.001) then r1 = hlybr else r1 = b1(1) + b1(2)*(f107a-71.5) + b1(3)*(f107-f107a+3.9) endif if (fexvir > 0.001) then r2 = fexvir else r2 = b2(1) + b2(2)*(f107a-71.5) + b2(3)*(f107-f107a+3.9) endif do l=13,lmax sflux(l) = (rflux(l) + ((r1-1.)*scale1(l) | + (r2-1.)*scale2(l)) / 1000.) enddo endif ! ! Tobiska EUV91 Method: ! if (iscale == 2) then if (hlya > 0.001) then hlymod = hlya else if (heiew > 0.001) then hlymod = heiew * 3.77847e9 + 8.40317e10 else hlymod = 8.70e8 * f107 + 1.90e11 endif endif if (heiew > 0.001) then heimod = heiew * 3.77847e9 + 8.40317e10 else heimod = hlymod endif do l=16,55 sflux(l) = tchr0(l) + tchr1(l)*hlymod + tchr2(l)*heimod | + tcor0(l) + tcor1(l)*f107 + tcor2(l)*f107a enddo endif ! ! Woods and Rottman (10 Nov. 1988) spectrum: ! if (iscale == 3) then do l=15,55 sflux(l) = war1(l) enddo endif ! ! Woods and Rottman (20 June 1989) spectrum: ! if (iscale == 4) then do l=15,55 sflux(l) = war2(l) enddo endif ! ! Substitute in H Lyman-alpha and XUVFAC if provided: ! if (hlya > 0.001) sflux(12) = hlya / 1.e9 if (xuvfac > 0.001) then xuvf = xuvfac else xuvf = 1.0 endif ! ! Convert from gigaphotons to photons, etc.: ! do l=1,lmax wave1(l) = wavel(l) wave2(l) = waves(l) if (sflux(l) < 0.0) sflux(l) = 0.0 if (wavel(l) < 251.0 .and. waves(l) > 15.0) | sflux(l)=sflux(l)*xuvf ! old model ! | sflux(l)=sflux(l)*4. ! Roble tune sflux(l) = sflux(l) * 1.e9 enddo ! ! write(6,"(/,' wave1 wave2 sflux')") ! do l=1,lmax ! write(6,"(3e12.4)") wave1(l),wave2(l),sflux(l) ! enddo ! ! Define module data sigabs and sigion from local ssigao, ssigao2, etc: ! sigabs(9,lmax), ! o,o2,n2,co,co2,o3,bco2,h2o,no ! sigion(13,lmax) ! o,o2,n2,co,co2,o4s,o2d,o2p,o4p,o2q,he,n,h ! if (itimet > 1) return ! (I don't think this ever happens) do l=1,lmax m=lmax-l+1 sigabs(1,l) = ssigao(l) * 1.E-18 sigabs(2,l) = ssigao2(l) * 1.E-18 sigabs(3,l) = ssigan2(l) * 1.E-18 sigabs(4,l) = ssigaco(m) * 1.E-18 sigabs(5,l) = ssigaco2(m)* 1.E-18 sigabs(6,l) = ssigao3(l) * 1.E-17 sigabs(7,l) = ssigbco2(l)* 1.E-18 sigabs(8,l) = ssigah2o(l)* 1.E-18 sigabs(9,l) = ssigano(l) * 1.E-18 ! sigion(1,l) = ssigio(l) * 1.E-18 sigion(2,l) = ssigio2(l) * 1.E-18 sigion(3,l) = ssigin2(l) * 1.E-18 sigion(4,l) = ssigico(m) * 1.E-18 sigion(5,l) = ssigico2(m)* 1.E-18 sigion(6,l) = ssigio4s(m)* 1.E-18 sigion(7,l) = ssigio2d(m)* 1.E-18 sigion(8,l) = ssigio2p(m)* 1.E-18 sigion(9,l) = ssigio4p(m)* 1.E-18 sigion(10,l)= ssigio2q(m)* 1.E-18 sigion(11,l)= ssigihe(m) * 1.E-18 sigion(12,l)= ssigin(m) * 1.E-17 sigion(13,l)= ssigih(m) * 1.E-18 enddo ! l=1,lmax itimet = itimet+1 ! ! 11/17/04 btf: no diffs with old code: ! ! call fprint8('From ssflux:',lmax, ! | sigabs(1,:),sigabs(2,:),sigabs(3,:),sigabs(4,:),sigabs(5,:), ! | sigabs(6,:),dum,dum, ! | 'sigabs(1)','sigabs(2)','sigabs(3)','sigabs(4)','sigabs(5)', ! | 'sigabs(6)',' ',' ') ! end subroutine ssflux !----------------------------------------------------------------------- subroutine euvac(euvflx) ! ! Args: real,intent(out) :: euvflx(ltorr) ! ! Local: integer :: i,ii real,dimension(ltorr) :: afac,f74113 real :: flxfac ! ! F74113 reference spectrum (doubled below 150-250 A, tripled <150) ! Will be multiplied by 1.0E9 later f74113 = (/ | 1.20 ,0.450 ,4.800 ,3.100 ,0.460 ,0.210 ,1.679 ,0.8, | 6.900 ,0.965 ,0.650 ,0.314 ,0.383 ,0.290 ,0.285 ,0.452 ,0.720, | 1.270 ,0.357 ,0.530 ,1.590 ,0.342 ,0.230 ,0.360 ,0.141 ,0.170, | 0.260 ,0.702 ,0.758 ,1.625 ,3.537 ,3.000 ,4.400 ,1.475 ,3.500, | 2.100 ,2.467/) ! ! Scaling factors(Ai) for the EUV flux afac = (/ | 1.0017E-02 ,7.1250E-03 ,1.3375E-02 ,1.9450E-02 ,2.7750E-03, | 1.3768E-01 ,2.6467E-02 ,2.5000E-02 ,3.3333E-03 ,2.2450E-02, | 6.5917E-03 ,3.6542E-02 ,7.4083E-03 ,7.4917E-03 ,2.0225E-02, | 8.7583E-03 ,3.2667E-03 ,5.1583E-03 ,3.6583E-03 ,1.6175E-02, | 3.3250E-03 ,1.1800E-02 ,4.2667E-03 ,3.0417E-03 ,4.7500E-03, | 3.8500E-03 ,1.2808E-02 ,3.2750E-03 ,4.7667E-03 ,4.8167E-03, | 5.6750E-03 ,4.9833E-03 ,3.9417E-03 ,4.4167E-03 ,5.1833E-03, | 5.2833E-03 ,4.3750E-03/) ! ! Loop through the wavelengths calculating the scaling factors and ! the resulting solar flux. The scaling factors are restricted to ! be greater than 0.8 ! do i=1,ltorr ii=i flxfac = (1.0 + afac(ii) * (0.5*(f107+f107a) - 80.0)) if(flxfac < 0.8) flxfac=0.8 euvflx(i) = f74113(ii) * flxfac * 1.0e9 enddo ! ! sigop2p, sigop2d, sigop4s are module data above ! (sigio and brxxxx are also module data) do ii=1,ltorr i=(ltorr+1)-ii sigop2p(i)=sigio(i)*brop2p(i) sigop2d(i)=sigio(i)*brop2d(i) sigop4s(i)=sigio(i)*brop4s(i) enddo ! end subroutine euvac !----------------------------------------------------------------------- subroutine srcj(fcol,jout) ! ! Formerly in yeesrc.F. ! ! Schumann Runge continuum parameterization from DeMajistre-Yee (2000) ! parameterization of the photolysis rate (J) and the energy deposition ! rate (eps) using UARS solstice solar flux. ! Inputs - ! fcol - column density of O2 (cm-2) ! f107 - F10.7 value (solar activity proxy) ! Output - ! jout - photolysis rate (s-1) ! ! Args: real,intent(in) :: fcol real,intent(out) :: jout ! ! Local: real :: a=1.031e-17 real :: ji(3) = (/2.118e-6, 5.658e-9, -9.179e-12/) real :: alpha(4,3) integer :: ic,jc real :: jinf,w ! alpha(:,1) = (/ 6.255e-1, 2.577e-1, 1.082e-1, 4.206e-3/) alpha(:,2) = (/ 2.488e-4, -1.659e-4, -7.647e-5, 7.783e-6/) alpha(:,3) = (/-4.871e-7, 3.577e-7, 1.281e-7, 1.752e-8/) ! jinf=0. do ic=1,3 jinf=jinf+ji(ic)*(f107**(ic-1)) enddo jout=0. do ic=1,4 w=0. do jc=1,3 w=w+alpha(ic,jc)*(f107**(jc-1)) enddo jout=jout+w*exp(-fcol*a*(4.**(1-ic))) enddo jout=jout*jinf end subroutine srcj !----------------------------------------------------------------------- subroutine srceps(fcol,eps) ! ! Formerly in yeesrc.F ! ! Inputs - ! fcol - column density of O2 (cm-2) ! f107 - F10.7 value (solar activity proxy) ! Output - ! eps - energy deposition rate (erg cm-3 s-1) ! ! Args: real,intent(in) :: fcol real,intent(out) :: eps ! ! Local: integer ic,jc real :: b=1.094e-17,einf,u real :: ei(3) = (/3.179e-18, 9.979e-21, -1.609e-23/) real :: beta(4,3) ! beta(:,1) = (/ 7.932e-1, 1.702e-1, 2.467e-1, 1.440e-3/) beta(:,2) = (/ 1.291e-4, -1.540e-4, 3.988e-5, -8.463e-6/) beta(:,3) = (/-3.110e-7, 3.815e-7, -1.092e-7, 2.474e-8/) ! einf=0. do ic=1,3 einf=einf+ei(ic)*(f107**(ic-1)) enddo ! eps=0. do ic=1,4 u=0. do jc=1,3 u=u+beta(ic,jc)*(f107**(jc-1)) enddo #ifdef EXPO eps=eps+u*expo(-fcol*b*(4.**(1-ic)),0) #else eps=eps+u*exp(-fcol*b*(4.**(1-ic))) #endif enddo eps=eps*einf end subroutine srceps !----------------------------------------------------------------------- subroutine uars_solaruv use restart,only: wv_low,wv_high,sminflx,smedflx,smaxflx ! uars data use guydat,only: nsigbin,nfsig,nfsigt,fsig,interp_sigt use flds_dissoc,only: dno2t,pdno2t,pdh2o2t,dh2o2t,pdh2o2t,dch4b, | dho2b,do3chap,po3chap,do3har,po3har,do3hz,do2hz,pdo2hz,pdo3hz, | do3srb,pdo2srb,do2srb,do3src,do2src,pdo2src,dnosrb use flds_dissoc,only: xjdis, ! (nlev,44), see fmodules.F | do3hug,po3hug use flds_heat,only: qchap,qpchap,qhar,qphar,qherz,qpherz2, | qpherz3,qpherz,qpsrb,qsrb,qsrc,qpsrc,qanosrb,qphug,qhug use flds_atmos,only: tn,rho,xno3,xno2,xnno use flds_ionzrt,only: xjch3oo,xjch2oa,xjch2ob use flds_col,only: txno2,txno3 use flds_modelz,only: zp,zpht ! for print only ! ! Local: integer :: l,k,j real,dimension(nlev) :: xxdis,qsheat,xxdis1,xxdis2,xxdis3, | qsheat1,qsheat2,qsheat3 real :: usflx(nsigbin) real :: fsigtn(nsigbin,nfsigt) ! interpolated Guy Brassier data integer :: ifsigtn(nfsigt) real :: tau,vla,tranhz1,tranhz2,alamsr,dum ! ! SOLAR MINIMUM UARS DAY 1417 - 29 JULY 1995 (95210) F107=72, F107A=77.59 ! SOLAR MAXIMUM UARS DAY 144 - 2 FEB. 1992 (92033) F107=280.1, F107A=202.1 ! do l=1,nsigbin ! was DO 100 ! usflx(l) = sminflx(l) usflx(l) = smedflx(l) ! usflx(l) = smaxflx(l) enddo ! l=1,nsigbin do k=1,nlev ! was DO 1 do j=1,44 xjdis(k,j) = 0. enddo enddo ifsigtn(:) = 1 ! ! Main k-loop: do k=1,nlev ! was DO 2 call interp_sigt(fsigtn,ifsigtn,tn(k),0,0.) ! fsigtn is returned do j=1,nfsig ! was DO 3 do l=1,nsigbin tau=fsig(l,1)*txno2(k)+fsigtn(l,1)*txno3(k) if (l /= 8) #ifdef EXPO | xjdis(k,j)=xjdis(k,j)+fsig(l,j)*usflx(l)*expo(-tau,0)*cgm #else | xjdis(k,j)=xjdis(k,j)+fsig(l,j)*usflx(l)*exp(-tau)*cgm #endif enddo enddo do j=nfsig+1,44 ! was DO 4 do l=1,nsigbin tau=fsig(l,1)*txno2(k)+fsigtn(l,1)*txno3(k) if (l /= 8) #ifdef EXPO | xjdis(k,j)=xjdis(k,j)+fsigtn(l,j-24)*usflx(l)* | expo(-tau,0)*cgm #else | xjdis(k,j)=xjdis(k,j)+fsigtn(l,j-24)*usflx(l)* | exp(-tau)*cgm #endif enddo enddo dno2t(k) = xjdis(k,28) pdno2t(k)=8.5e-3*exp(-1.0e-20*txno3(k)) pdh2o2t(k)=1.13e-4*exp(-5.0e-18*txno3(k))*cgm dh2o2t(k) = xjdis(k,27) pdh2o2t(k)=1.13e-4*exp(-5.0e-18*txno3(k))*cgm dch4b(k) = 1.e-20 dho2b(k) = xjdis(k,5) xjch3oo(k) = xjdis(k,17) xjch2oa(k) = xjdis(k,42) xjch2ob(k) = xjdis(k,43) xxdis(k)=0. qsheat(k)=0. do l=114,nsigbin ! was DO 5 vla=(wv_high(l)+wv_low(l))*0.5*10. tau=fsig(l,1)*txno2(k)+fsigtn(l,1)*txno3(k) xxdis(k)=xxdis(k)+fsigtn(l,1)*usflx(l)*exp(-tau)*cgm qsheat(k)=qsheat(k)+fsigtn(l,1)*usflx(l)*xno3(k)*exp(-tau)* | cgm*(12400./vla-1.09)*erg/rho(k) enddo do3chap(k)=xxdis(k) po3chap(k)=3.34e-4*exp(-3.16e-21*txno3(k))*cgm qchap(k)=qsheat(k) qpchap(k)=5.512e-16*exp(-3.16e-21*txno3(k))*xno3(k)/rho(k)*cgm xxdis(k)=0. qsheat(k)=0. do l=95,113 ! was DO 6 vla=(wv_high(l)+wv_low(l))*0.5*10. tau=fsig(l,1)*txno2(k)+fsigtn(l,1)*txno3(k) xxdis(k)=xxdis(k)+fsigtn(l,1)*usflx(l)*exp(-tau)*cgm qsheat(k)=qsheat(k)+fsigtn(l,1)*usflx(l)*xno3(k)*exp(-tau)* | cgm*(12400./vla-1.09)*erg/rho(k) enddo do3hug(k) = xxdis(k) po3hug(k) = (7.9444e-5*exp(-8.94e-19*txno3(k))+ | 2.1012e-5*exp(-1.09e-19*txno3(k))+ | 1.0188e-4*exp(-3.07e-18*txno3(k)))*cgm qhug(k)=qsheat(k) qphug(k)=(41.21e-16*exp(-8.94e-19*txno3(k))+ | 10.9 e-16*exp(-1.09e-19*txno3(k))+ | 52.80e-16*exp(-3.07e-18*txno3(k)))*cgm*xno3(k)/rho(k) xxdis(k)=0. qsheat(k)=0. do l=76,94 ! was DO 7 vla=(wv_high(l)+wv_low(l))*0.5*10. tau=fsig(l,1)*txno2(k)+fsigtn(l,1)*txno3(k) xxdis(k)=xxdis(k)+fsigtn(l,15)*usflx(l)*exp(-tau)*cgm qsheat(k)=qsheat(k)+fsigtn(l,15)*usflx(l)*xno3(k)*exp(-tau)* | cgm*(12400./vla-1.09-1.96-0.98)*erg/rho(k) enddo do3har(k)=xxdis(k) po3har(k)=(4.697e-3*exp(-9.65e-18*txno3(k))+6.33e-4* | exp(-4.79e-18*txno3(k)))*cgm qhar(k)=qsheat(k) qphar(k)=(69.67e-16*exp(-9.65e-18*txno3(k))+9.39e-16* | exp(-4.79e-18*txno3(k)))*cgm*xno3(k)/rho(k) xxdis(k)=0. xxdis1(k)=0. qsheat(k)=0. qsheat1(k)=0. do l=59,75 ! was DO 8 vla=(wv_high(l)+wv_low(l))*0.5*10. tau=fsig(l,1)*txno2(k)+fsigtn(l,1)*txno3(k) xxdis(k) =xxdis(k)+fsigtn(l,15)*usflx(l)*exp(-tau)*cgm xxdis1(k)=xxdis1(k)+fsig(l,1)*usflx(l)*exp(-tau)*cgm qsheat(k)=qsheat(k)+fsigtn(l,15)*usflx(l)*xno3(k)*exp(-tau)* | cgm*(12400./vla-1.09-1.96-0.98)*erg/rho(k) qsheat1(k)=qsheat1(k)+fsig(l,1)*usflx(l)*xno2(k)*exp(-tau)* | cgm*(12400./vla-5.12)*erg/rho(k) enddo do3hz(k)=xxdis(k) do2hz(k)=xxdis1(k) tranhz1=exp(-5.50e-24*txno2(k)-6.22e-18*txno3(k)) tranhz2=exp(-1.34e-26*txno2(k)-1.66e-18*txno3(k)) pdo2hz(k)=(6.545e-10*tranhz1+1.101e-12*tranhz2)* | cgm*0.8536 pdo3hz(k)=(6.593e-4*tranhz1+1.216e-4*tranhz2)*cgm qherz(k)=qsheat(k)+qsheat1(k) qpherz2(k)=(2.74e-22*tranhz1+4.61e-25*tranhz2)* | cgm*xno2(k)/rho(k)*1.494 qpherz3(k)=(14.94e-16*tranhz1+2.76e-16*tranhz2)* | cgm*xno3(k)/rho(k) qpherz(k) = qpherz2(k)+qpherz3(k) xxdis(k)=0. xxdis1(k)=0. xxdis2(k)=0. xxdis3(k)=0. qsheat(k)=0. qsheat1(k)=0. qsheat2(k)=0. qsheat3(k)=0. do l=45,58 ! was DO 9 vla=(wv_high(l)+wv_low(l))*0.5*10. tau=fsig(l,1)*txno2(k)+fsigtn(l,1)*txno3(k) xxdis(k)=xxdis(k)+fsigtn(l,15)*usflx(l)*exp(-tau)*cgm xxdis1(k)=xxdis1(k)+fsig(l,1)*usflx(l)*exp(-tau)*cgm xxdis2(k)=xxdis2(k)+fsig(l,2)*usflx(l)*exp(-tau)*cgm xxdis3(k)=xxdis3(k)+fsig(l,3)*usflx(l)*exp(-tau)*cgm qsheat(k)=qsheat(k)+fsigtn(l,15)*usflx(l)*xno3(k)*exp(-tau)* | cgm*(12400./vla-1.09-1.96-0.98)*erg/rho(k) qsheat1(k)=qsheat1(k)+fsig(l,1)*usflx(l)*xno2(k)*exp(-tau)* | cgm*(12400./vla-5.12)*erg/rho(k) enddo do3srb(k)=xxdis(k) pdo2srb(k)=1.1e-7*exp(-1.97e-10*txno2(k)**0.522)*cgm if(txno2(k).gt.1.e+19) pdo2srb(k)=1.45e+8/(txno2(k)**0.83)*cgm qpsrb(k)=2.833e-12*xno2(k)*pdo2srb(k)/rho(k) do2srb(k)=pdo2srb(k) qsrb(k)=qpsrb(k) xxdis(k)=0. xxdis1(k)=0. xxdis2(k)=0. xxdis3(k)=0. qsheat(k)=0. qsheat1(k)=0. qsheat2(k)=0. qsheat3(k)=0. do l=1,44 ! was DO 10 if (l /= 8) then vla=(wv_high(l)+wv_low(l))*0.5*10. tau=fsig(l,1)*txno2(k)+fsigtn(l,1)*txno3(k) #ifdef EXPO xxdis(k)=xxdis(k)+fsigtn(l,15)*usflx(l)*expo(-tau,0)*cgm xxdis1(k)=xxdis1(k)+fsig(l,1)*usflx(l)*expo(-tau,0)*cgm xxdis2(k)=xxdis2(k)+fsig(l,2)*usflx(l)*expo(-tau,0)*cgm xxdis3(k)=xxdis3(k)+fsig(l,3)*usflx(l)*expo(-tau,0)*cgm qsheat(k)=qsheat(k)+fsigtn(l,15)*usflx(l)*xno3(k)* | expo(-tau,0)*cgm*(12400./vla-1.09-1.96-0.98)*erg/rho(k) qsheat1(k)=qsheat1(k)+fsig(l,1)*usflx(l)*xno2(k)* | expo(-tau,0)*cgm*(12400./vla-7.12)*erg/rho(k) #else xxdis(k)=xxdis(k)+fsigtn(l,15)*usflx(l)*exp(-tau)*cgm xxdis1(k)=xxdis1(k)+fsig(l,1)*usflx(l)*exp(-tau)*cgm xxdis2(k)=xxdis2(k)+fsig(l,2)*usflx(l)*exp(-tau)*cgm xxdis3(k)=xxdis3(k)+fsig(l,3)*usflx(l)*exp(-tau)*cgm qsheat(k)=qsheat(k)+fsigtn(l,15)*usflx(l)*xno3(k)*exp(-tau)* | cgm*(12400./vla-1.09-1.96-0.98)*erg/rho(k) qsheat1(k)=qsheat1(k)+fsig(l,1)*usflx(l)*xno2(k)*exp(-tau)* | cgm*(12400./vla-7.12)*erg/rho(k) #endif endif enddo do3src(k)=xxdis(k) do2src(k)=xxdis1(k) qsrc(k)=qsheat(k)+qsheat1(k) #ifdef EXPO pdo2src(k)=(0.98*exp(-2.9e-19*txno2(k))-0.55*expo(-1.7e-18* | txno2(k),0)-0.43*expo(-1.15e-17*txno2(k),0))*cgm/txno2(k)* | 4.4286e+11 #else pdo2src(k)=(0.98*exp(-2.9e-19*txno2(k))-0.55*exp(-1.7e-18* | txno2(k))-0.43*exp(-1.15e-17*txno2(k)))*cgm/txno2(k)* | 4.4286e+11 #endif alamsr=1576. qpsrc(k)=pdo2src(k)*xno2(k)*1.602e-12*(12400./alamsr-7.12)/ | rho(k) #ifdef EXPO dnosrb(k)=7.0e-6*exp(-1.e-8*txno2(k)**0.38)*expo(-5.e-19* | txno3(k),0)*cgm #else dnosrb(k)=7.0e-6*exp(-1.e-8*txno2(k)**0.38)*exp(-5.e-19* | txno3(k))*cgm #endif qanosrb(k)=dnosrb(k)*xnno(k)*(12400./1900.-6.51)*erg/rho(k) enddo ! k=1,nlev ! ! 11/17/04 btf: no diffs with old code: ! call fprint8('uars_solaruv:',nlev, ! | zp , zpht , do3chap , po3chap , qchap , qpchap,dum,dum, ! | 'zp','zpht','do3chap','po3chap','qchap','qpchap',' ',' ') ! call fprint8('uars_solaruv:',nlev, ! | do3src , do2src , qsrc , pdo2src , qpsrc , dnosrb , qanosrb, dum, ! |'do3src','do2src','qsrc','pdo2src','qpsrc','dnosrb','qanosrb',' ') end subroutine uars_solaruv !----------------------------------------------------------------------- real function e1(x) ! ! 2/14/05 btf: this function moved from old solvers.F. ! implicit none ! ! Args: real,intent(in) :: x ! ! Local: integer :: ir real :: sum real :: tr(17),trr(12) real :: ar(17) = | (/3.9368857696, -8.0314874286, 3.8797325768, -1.6042971072, | 0.5630905453, -0.1704423017, 0.0452099390, -0.0106538986, | 0.0022562638, -0.0004335700, 0.0000762166, -0.0000123417, | 0.0000018519, -0.0000002588, 0.0000000338, -0.0000000041, | 0.0000000004/) real :: br(12) = | (/ 0.1077641888, 0.1028106215, -0.0045526707, 0.0003571613, | -0.0000379341, 0.0000049143, -0.0000007355, 0.0000001230, | -0.0000000225, 0.0000000044, -0.0000000009, 0.0000000002/) ! print *,'E1: x=',x IF (X.GT.4.) GO TO 11 TR(1)=1. TR(2)=X/4. DO 1 IR=3,17 TR(IR)=2.*X/4*TR(IR-1)-TR(IR-2) 1 CONTINUE SUM=0. DO 2 IR=1,17 2 SUM=SUM+AR(IR)*TR(IR) E1=-ALOG(ABS(X))-SUM RETURN 11 CONTINUE TRR(1)=1. TRR(2)=2.*4./X-1. DO 3 IR=3,12 3 TRR(IR)=(4.*4./X-2.)*TRR(IR-1)-TRR(IR-2) SUM=0. DO 4 IR=1,12 4 SUM=SUM+BR(IR)*TRR(IR) ! print *,'E1: x=',x,' sum=',sum,' exp(-x)=',exp(-x) ! print *,'E1: x=',x,' sum=',sum #ifdef EXPO e1=expo(-x,0)*sum #else e1=exp(-x)*sum #endif ! print *,' e1 returning: e1=',e1 end function e1 !----------------------------------------------------------------------- end module solarheat