! module qrj_module ! ! This software is part of the NCAR TIE-GCM. Use is governed by the ! Open Source Academic Research License Agreement contained in the file ! tiegcmlicense.txt. ! use params_module,only: nlevp1,nlonp4 use addfld_module,only: addfld implicit none ! integer,parameter :: lmax=37 ! number of bins from 0.5A-1750A integer,parameter :: l1=15 ! number of bins in (1050-1750A) real :: | euveff(nlevp1), | sigeuv(3,lmax), | rlmeuv(lmax), | feuv(lmax), | fsrc(l1), | sigsrc(l1), | rlmsrc(l1), | sigin4s(lmax), | quench(4), | wave1(lmax), ! short bound of wave bins | wave2(lmax), ! long bound of wave bins | sfmin(lmax), ! reference solar minimum flux of EUVAC model | afac(lmax), ! The A factor of EUVAC model | sflux(lmax) ! Solar flux for each time step ! ! Branching ratios for photon are branching ratios from photon absorption rate. ! Branching ratios for photolectron are branching ratios from photoionization rate. real:: | BPhotonI(3,lmax), ! Photoionization branching ratio for major species | BElectronI(3,lmax), ! Photoelectron ionization branching ratio ! for major species | brop2pPh(lmax), ! photoionization branching ratio for O+(2p) | brop2dPh(lmax), ! photoionization branching ratio for O+(2d) | brop4sPh(lmax), ! photoionization branching ratio for O+(4s) | bro2DPh(lmax), ! photodissociation branching ratio for O2 | brn2DPh(lmax), ! photodissociation branching ratio for N2 | bro2DIPh(lmax), ! Photon dissociative ionization branching ratio for O2 | brn2DIPh(lmax), ! Photon dissociative ionization branching ratio for N2 | brop2pEl(lmax), ! electron impact ionization branching ratio for O+(2p) | brop2dEl(lmax), ! electron impact ionization branching ratio for O+(2d) | brop4sEl(lmax), ! electron impact ionization branching ratio for O+(4s) | bro2DIEl(lmax), ! Photoelectron dissociative ionization branching ratio for O2 | brn2DIEl(lmax), ! Photoelectron dissociative ionization branching ratio for N2 | brn2DEl(lmax), ! Photoelectron dissociation branching ratio for N2 | bro2DEl(lmax) ! Photoelectron dissociation branching ratio for O2 ! ! Heating and ionization terms set by qrj, and used by other routines. ! These are allocated for task subdomains by alloc_q (called from allocdata) ! real,dimension(:,:,:),allocatable :: | rj, ! total o2 dissociation frequency (s^-1) | qtef, ! total N dissociative production rate (cm^-3 s^-1) | qtotal, ! total heating rate | qop2p, ! o+(2p) production rate | qop2d, ! o+(2d) production rate | qo2p, ! o2+ production rate | qop, ! o+ production rate | qn2p, ! n2+ production rate | qnp, ! n+ production rate | qnop ! no+ production rate ! contains !----------------------------------------------------------------------- subroutine qrj(sco2,sco1,scn2,tn,no,o2,o1,n4s,xnmbari, | lev0,lev1,lon0,lon1,lat) ! ! Calculate heating and dissociation rates. ! use input_module,only: f107,f107a use init_module,only: sfeps use cons_module,only: t0,expz,avo,rmassinv_n4s,rmassinv_no, | rmassinv_o2,rmassinv_o1,rmassinv_n2,expzmid,gask, | expzmid_inv,p0,check_exp use lbc,only: fb,b use chemrates_module,only: beta9 use fields_module,only: tlbc ! ! Args: integer,intent(in) :: lev0,lev1,lon0,lon1,lat real,dimension(lev0:lev1,lon0-2:lon1+2),intent(in) :: | sco2,sco1,scn2, ! chapman integrals | tn,no,o2,o1,n4s, ! tn and species mass mixing ratios | xnmbari ! p0*e(-z)*barm/kT at interfaces ! ! VT vampir tracing: ! #ifdef VT #include #endif ! ! Local: integer :: k,i,l,ier real,parameter :: | do2=8.203E-12 , | do22=1.1407E-11 , | aband=0.143 , ! shumann-runge | bband=9.64E8 , ! shumann-runge | cband=9.03E-19 , ! shumann-runge | e3=0.33 , | hc = 1.9845E-16 ! C(60) real :: rlmeuvinv(lmax),rlmsrcinv(l1),factor,fmin,fmax real,dimension(lev0:lev1,lon0:lon1) :: | o2i, ! o2 at interfaces (s1) | o1i, ! o at interfaces (s2) | n2i, ! n2 at interfaces (s3) | n4si, ! n4s at interfaces (s4) | tni ! tn at interfaces (s6) real,dimension(lev0:lev1,lon0:lon1) :: | quenchfac, ! (s8) | sigchap, ! (s9) | p3f ! (s7) ! temporary loop variables real:: | absorp_o, ! photoabsorption frequency of O | absorp_o2, ! photoabsoption frequency of o2 | absorp_n2, ! photoabsoption frequency of n2 | ioniz_o, ! photoionization frequency of o | ioniz_o2, ! photoionization frequency of o2 | ioniz_n2, ! photoionization frequency of n2 | htfac ! hc/wavelength, for calculation of heating real,dimension(lev0:lev1,lon0:lon1) :: | di_o2, ! total dissociative ionization frequency of o2 | di_n2, ! total dissociative ionization frequency of n2 | mn_o2, ! transfer mass density to number density (O2) | mn_o1, ! transfer mass density to number density (O) | mn_n2, ! transfer mass density to number density (N2) | mn_n ! transfer mass density to number density (N) real,dimension(lev0:lev1,lon0:lon1) :: | sum1, ! sum(o2,o,n2)(sigma*chapman) (s5) | sum2, ! sum(o2,o,n2)(sigma*psi/rmass) (s6) | sum3 ! sum(o2,o,n2,n4s)(sigmas) (s7) logical,parameter :: | debug=.false. ! insert print statements ! ! expo() (util.F) is used only if check_exp is true. This will avoid ! NaNS fpe, but will degrade performance. Check_exp is in cons.F. ! real,external :: expo ! used only when check_exp is set (util.F) ! if (debug) write(6,"('Enter qrj: lat=',i3,' lon0,1=',2i3)") | lat,lon0,lon1 ! #ifdef VT call vtbegin(118,ier) #endif ! ! Exec: ! call addfld('XNMBARI',' ',' ',xnmbari(:,lon0:lon1), ! | 'lev',lev0,lev1,'lon',lon0,lon1,lat) ! call addfsech('XNMBARI',' ',' ',xnmbari(:,lon0:lon1), ! | lon0,lon1,nlevs,nlevs,lat) ! ! calculate inverse of wave length ! do i=1,lmax rlmeuvinv(i) = 1./rlmeuv(i) enddo do i=1,l1 rlmsrcinv(i) = 1./rlmsrc(i) enddo ! ! O2,O,N4S at interface levels: do i=lon0,lon1 do k=lev0,lev1-1 o2i (k+1,i) = 0.5*(o2(k,i)+o2(k+1,i)) o1i (k+1,i) = 0.5*(o1(k,i)+o1(k+1,i)) n4si(k+1,i) = 0. enddo enddo ! ! Bottom boundary: do i=lon0,lon1 o2i(1,i) = .5*((b(i,1,1)+1.)*o2(1,i)+b(i,1,2)*o1(1,i)+ | fb(i,1)) o1i(1,i) = .5*(b(i,2,1)*o2(1,i)+(b(i,2,2)+1.)*o1(1,i)+ | fb(i,2)) n4si(1,i) = 0. enddo ! ! N2 at interfaces: do k=lev0,lev1 n2i(k,:) = 1.-o2i(k,:)-o1i(k,:) do i=lon0,lon1 !!! temporal fix CISM if(n2i(k,i)<=0.)n2i(k,i)=1.e-5 !!! temporal fix CISM enddo !!! temporal fix CISM enddo ! ! calculate variables for transferring mass density to number density do i=lon0,lon1 do k=lev0,lev1 mn_o2(k,i)=o2i(k,i)*rmassinv_o2 mn_o1(k,i)=o1i(k,i)*rmassinv_o1 mn_n2(k,i)=n2i(k,i)*rmassinv_n2 mn_n(k,i)=n4si(k,i)*rmassinv_n4s enddo enddo ! ! call addfld('O2I' ,' ',' ',o2i, ! | 'lev',lev0,lev1,'lon',lon0,lon1,lat) ! call addfld('O1I' ,' ',' ',o1i, ! | 'lev',lev0,lev1,'lon',lon0,lon1,lat) ! call addfld('N2I' ,' ',' ',n2i, ! | 'lev',lev0,lev1,'lon',lon0,lon1,lat) ! ! Initialize ionization/dissociation arrays on current processor domain: ! (global module data for use by other routines) rj (lev0:lev1,lon0:lon1,lat) = 0. qtef (lev0:lev1,lon0:lon1,lat) = 0. qtotal(lev0:lev1,lon0:lon1,lat) = 0. qop2p (lev0:lev1,lon0:lon1,lat) = 0. qop2d (lev0:lev1,lon0:lon1,lat) = 0. qo2p (lev0:lev1,lon0:lon1,lat) = 0. qop (lev0:lev1,lon0:lon1,lat) = 0. qn2p (lev0:lev1,lon0:lon1,lat) = 0. qnp (lev0:lev1,lon0:lon1,lat) = 0. qnop (lev0:lev1,lon0:lon1,lat) = 0. ! ! Initialize local arrays di_o2=0. di_n2=0. ! ! Summation over wavelength: do l=l1+1,lmax ! from 0.5A to 1050A sum1(:,:) = 0. ! sum(o2,o,n2)(sigma*chapman) sum2(:,:) = 0. ! sum(o2,o,n2)(sigma*psi/rmass) sum3(:,:) = 0. ! sum(o2,o,n2,n4s)(sigmas) htfac=hc*rlmeuvinv(l) do i=lon0,lon1 do k=lev0,lev1 sum1(k,i) = sum1(k,i)+sigeuv(1,l)*sco2(k,i)+ | sigeuv(2,l)*sco1(k,i)+ | sigeuv(3,l)*scn2(k,i) ! if (.not.check_exp) then sum1(k,i) = feuv(l)*exp(-sum1(k,i)) else sum1(k,i) = feuv(l)*expo(-sum1(k,i),0) endif sum2(k,i) = sum2(k,i)+sigeuv(1,l)*mn_o2(k,i)+ | sigeuv(2,l)*mn_o1(k,i)+ | sigeuv(3,l)*mn_n2(k,i) enddo enddo ! ! Longitude and column domain of current process: do i=lon0,lon1 do k=lev0,lev1 ! absorption/ionization frequency for the three major species (O2, O, and N2) absorp_o2= sum1(k,i)*sigeuv(1,l) absorp_o= sum1(k,i)*sigeuv(2,l) absorp_n2= sum1(k,i)*sigeuv(3,l) ioniz_o2=absorp_o2*BPhotonI(1,l) ioniz_o=absorp_o*BPhotonI(2,l) ioniz_n2=absorp_n2*BPhotonI(3,l) ! ! ionization/dissociative ionization frequency (s^-1) di_o2(k,i)=di_o2(k,i)+absorp_o2*bro2DIPh(l)+ | ioniz_o2*bro2DIEl(l) di_n2(k,i)=di_n2(k,i)+absorp_n2*brn2DIPh(l)+ | ioniz_n2*brn2DIEl(l) qnp(k,i,lat) = qnp(k,i,lat)+sigin4s(l)*sum1(k,i) qn2p(k,i,lat) = qn2p(k,i,lat)+ioniz_n2+ | ioniz_n2*BElectronI(3,l) qo2p(k,i,lat) = qo2p(k,i,lat)+ioniz_o2+ | ioniz_o2*BElectronI(1,l) qop2p(k,i,lat) = qop2p(k,i,lat)+absorp_o*brop2pPh(l)+ | ioniz_o*brop2pEl(l) qop2d(k,i,lat) = qop2d(k,i,lat)+absorp_o*brop2dPh(l)+ | ioniz_o*brop2dEl(l) qop(k,i,lat) = qop(k,i,lat)+absorp_o*brop4sPh(l)+ | ioniz_o*brop4sEl(l) ! ! total dissociation and EUV heating rj(k,i,lat) = rj(k,i,lat)+(absorp_o2*bro2DPh(l) | +ioniz_o2*bro2DEl(l)) qtef(k,i,lat) = qtef(k,i,lat)+absorp_n2*brn2DPh(l) | +ioniz_n2*brn2DEl(l) qtotal(k,i,lat) = qtotal(k,i,lat)+htfac* | sum1(k,i)*sum2(k,i) enddo enddo enddo ! l=l1+1,lmax ! ! call addfld('QOP2P' ,' ',' ', qop2p(lev0:lev1,lon0:lon1,lat), ! | 'lev',lev0,lev1,'lon',lon0,lon1,lat) ! call addfld('QOP2D' ,' ',' ', qop2d(lev0:lev1,lon0:lon1,lat), ! | 'lev',lev0,lev1,'lon',lon0,lon1,lat) ! call addfld('SUM1' ,' ',' ', sum1, ! | 'lev',lev0,lev1,'lon',lon0,lon1,lat) ! call addfld('Q0' ,' ',' ',qtotal(:,:,lat), ! | 'lev',lev0,lev1,'lon',lon0,lon1,lat) ! call addfld('RJ_QRJ',' ',' ',rj(:,:,lat), ! | 'lev',lev0,lev1,'lon',lon0,lon1,lat) ! ! Multiply Q by efficiency factor: do i=lon0,lon1 do k=lev0,lev1 qtotal(k,i,lat) = qtotal(k,i,lat)*euveff(k)*avo enddo enddo ! call addfld('Q1' ,' ',' ',qtotal(:,:,lat), ! | 'lev',lev0,lev1,'lon',lon0,lon1,lat) ! ! transfer frequency (s^-1) to production rates (cm^-3 s^-1) do i=lon0,lon1 do k=lev0,lev1 mn_o2(k,i)=mn_o2(k,i)*xnmbari(k,i) mn_n2(k,i)=mn_n2(k,i)*xnmbari(k,i) mn_o1(k,i)=mn_o1(k,i)*xnmbari(k,i) mn_n(k,i)=mn_n(k,i)*xnmbari(k,i) di_o2(k,i)=di_o2(k,i)*mn_o2(k,i) di_n2(k,i)=di_n2(k,i)*mn_n2(k,i) qo2p(k,i,lat) = qo2p(k,i,lat)*mn_o2(k,i)-di_o2(k,i) qn2p(k,i,lat) = qn2p(k,i,lat)*mn_n2(k,i)-di_n2(k,i) qnp(k,i,lat) = qnp(k,i,lat)*mn_n(k,i)+di_n2(k,i) qop(k,i,lat) = qop(k,i,lat)*mn_o1(k,i)+0.54*di_o2(k,i) qop2p(k,i,lat) = qop2p(k,i,lat)*mn_o1(k,i)+0.22*di_o2(k,i) qop2d(k,i,lat) = qop2d(k,i,lat)*mn_o1(k,i)+0.24*di_o2(k,i) qtef(k,i,lat) = 2.*qtef(k,i,lat)*mn_n2(k,i)+di_n2(k,i) enddo enddo ! ! call addfld('QOP',' ',' ',qop(:,lon0:lon1,lat), ! | 'lev',lev0,lev1,'lon',lon0,lon1,lat) ! ! Add no ionization to qnop: do i=lon0,lon1 qnop(lev0,i,lat) = qnop(lev0,i,lat)+ | beta9(lev0,i,lat)*no(lev0,i)*xnmbari(lev0,i)*rmassinv_no enddo do i=lon0,lon1 do k=lev0+1,lev1 qnop(k,i,lat) = qnop(k,i,lat)+beta9(k,i,lat)*.5*(no(k,i)+ | no(k-1,i))*xnmbari(k,i)*rmassinv_no enddo enddo ! ! tn at interfaces: do i=lon0,lon1 tni(lev0,i) = tlbc(i,lat) do k=lev0+1,lev1-1 tni(k,i) = .5*(tn(k-1,i)+tn(k,i)) enddo tni(lev1,i) = tn(lev1-1,i) ! nlevp1 <- nlev enddo ! ! Quench: do i=lon0,lon1 do k=lev0,lev1 factor = avo*p0/gask*expz(1)*expzmid**(2*k-3) quenchfac(k,i) = factor/((o2i(k,i)*rmassinv_o2+ | o1i(k,i)*rmassinv_o1+ | n2i(k,i)*rmassinv_n2)*tni(k,i)) quenchfac(k,i) = quenchfac(k,i)* | (quench(1)*n2i(k,i)*rmassinv_n2+ | quench(2)*o2i(k,i)*rmassinv_o2) quenchfac(k,i) = quench(3)*quenchfac(k,i)/ | (quench(4)+quenchfac(k,i)) enddo enddo ! ! Summation over wave length: sum1(:,:) = 0. do l=1,l1 ! from 1050A to 1750A do i=lon0,lon1 do k=lev0,lev1 ! ! Note check_exp should be true for debug runs only: if (.not.check_exp) then sigchap(k,i) = sigsrc(l)*fsrc(l)*exp(-sigsrc(l)*sco2(k,i)) else sigchap(k,i) = sigsrc(l)*fsrc(l)* | expo(-sigsrc(l)*sco2(k,i),0) endif ! sum1(k,i) = sum1(k,i)+sigchap(k,i)* | (hc*rlmsrcinv(l)-do22+quenchfac(k,i)) rj(k,i,lat) = rj(k,i,lat)+sigchap(k,i) enddo enddo enddo ! ! Update q: do i=lon0,lon1 do k=lev0,lev1 qtotal(k,i,lat) = qtotal(k,i,lat)+sum1(k,i)*avo* | o2i(k,i)*rmassinv_o2 enddo enddo ! ! call addfld('Q2' ,' ',' ',qtotal(:,:,lat), ! | 'lev',lev0,lev1,'lon',lon0,lon1,lat) ! ! Contributions from schumann runge bands: do i=lon0,lon1 do k=lev0,lev1 if (sco2(k,i) >= 1.e18) then p3f(k,i) = (1./(aband*sco2(k,i)+bband*sqrt(sco2(k,i))))* | (1.+0.11*(f107-65.)/165.)*sfeps else p3f(k,i) = cband*(1.+0.11*(f107-65.)/165.)*sfeps endif qtotal(k,i,lat) = qtotal(k,i,lat)+p3f(k,i)*avo* | o2i(k,i)*rmassinv_o2*e3 rj(k,i,lat) = rj(k,i,lat)+p3f(k,i)/do2 enddo enddo ! call addfld('Q' ,' ',' ',qtotal(lev0:lev1,lon0:lon1,lat), ! | 'lev',lev0,lev1,'lon',lon0,lon1,lat) ! call addfld('RJ' ,' ',' ',rj(lev0:lev1,lon0:lon1,lat), ! | 'lev',lev0,lev1,'lon',lon0,lon1,lat) ! call addfld('QOP2Pa' ,' ',' ',qop2p(lev0:lev1,lon0:lon1,lat), ! | 'lev',lev0,lev1,'lon',lon0,lon1,lat) ! call addfld('QOP2Da' ,' ',' ', qop2d(lev0:lev1,lon0:lon1,lat), ! | 'lev',lev0,lev1,'lon',lon0,lon1,lat) ! call addfld('QTEF' ,' ',' ', qtef(lev0:lev1,lon0:lon1,lat), ! | 'lev',lev0,lev1,'lon',lon0,lon1,lat) ! call addfld('QO2P' ,' ',' ', qo2p(lev0:lev1,lon0:lon1,lat), ! | 'lev',lev0,lev1,'lon',lon0,lon1,lat) ! call addfld('QN2P' ,' ',' ', qn2p(lev0:lev1,lon0:lon1,lat), ! | 'lev',lev0,lev1,'lon',lon0,lon1,lat) ! call addfld('QNP' ,' ',' ', qnp(lev0:lev1,lon0:lon1,lat), ! | 'lev',lev0,lev1,'lon',lon0,lon1,lat) ! call addfld('QOPa' ,' ',' ', qop(lev0:lev1,lon0:lon1,lat), ! | 'lev',lev0,lev1,'lon',lon0,lon1,lat) ! call addfld('QNOP' ,' ',' ', qnop(lev0:lev1,lon0:lon1,lat), ! | 'lev',lev0,lev1,'lon',lon0,lon1,lat) ! #ifdef VT ! code = 118 ; state = 'qrj' ; activity='ModelCode' call vtend(118,ier) #endif ! end subroutine qrj !--------------------------------------------------------------- subroutine init_sflux use input_module,only: f107,f107a,see_ncfile use init_module,only: sfeps,iday,iyear,secs,istep use soldata_module,only: get_soldata,soldata,nwave ! ! Flux initialization once per time step, called from advance. ! ! Local: integer :: n ! sflu: scaled solar flux returned by subroutine ssflux() (photons cm-2 s-1) ! if (len_trim(see_ncfile) == 0) then call ssflux(f107,f107a) do n = l1+1,lmax feuv(n) = sflux(n)*sfeps enddo do n = 1,l1 fsrc(n) = sflux(n)*sfeps enddo else call get_soldata(iyear,iday,int(secs),istep) if (nwave /= lmax) then write(6,"('init_sflux(): wave bins mismatch: nwave=',i4, | ' lmax=',i4)") nwave,lmax call shutdown('init_sflux') endif do n=1,l1 fsrc(n)=soldata(n) enddo do n=l1+1,lmax feuv(n)=soldata(n) enddo endif end subroutine init_sflux !----------------------------------------------------------------------- subroutine init_qrj integer :: m,n ! ! Called once per run, from tgcm.F ! ! Initialize bins (37 bins) ! wave1 = (/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, | 1027.00, 987.00, 975.00, 913.00, 913.00, | 913.00, 798.00, 798.00, 798.00, 650.00, | 650.00, 540.00, 320.00, 290.00, 224.00, | 155.00, 70.00, 32.00, 18.00, 8.00, | 4.00, 0.50/) wave2 = (/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, 1027.00, 987.00, 975.00, 975.00, | 975.00, 913.00, 913.00, 913.00, 798.00, | 798.00, 650.00, 540.00, 320.00, 290.00, | 224.00, 155.00, 70.00, 32.00, 18.00, | 8.00, 4.00/) ! ! Solar spectrum based on EUVAC and glow for wave length less than 1050 A ! and Woods for wavelength greater than 1050 A ! ! solar minimum flux (when P_index=80, unit:photon cm^-2 S^-1) ! sfmin=(/3.397e+11, 1.998e+11, 1.055e+11, 7.260e+10, | 5.080e+10, 2.802e+10, 1.824e+10, 1.387e+10, | 2.659e+10, 7.790e+09, 1.509e+10, 3.940e+11, | 8.399e+09, 3.200e+09, 3.298e+09, 4.235e+09, | 4.419e+09, 4.482e+09, 7.156e+08, 1.028e+09, | 3.818e+08, 8.448e+08, 3.655e+09, 2.364e+09, | 1.142e+09, 1.459e+09, 4.830e+09, 2.861e+09, | 8.380e+09, 4.342e+09, 5.612e+09, 1.270e+09, | 5.326e+08, 2.850e+07, 2.000e+06, 1.000e+04, | 5.010e+01/) ! ! scaling factor A as defined in EUVAC model ! afac=(/5.937e-04, 6.089e-04, 1.043e-03, 1.125e-03, | 1.531e-03, 1.202e-03, 1.873e-03, 2.632e-03, | 2.877e-03, 2.610e-03, 3.739e-03, 4.230e-03, | 2.541e-03, 2.099e-03, 3.007e-03, 4.825e-03, | 5.021e-03, 3.950e-03, 4.422e-03, 4.955e-03, | 4.915e-03, 5.437e-03, 5.261e-03, 5.310e-03, | 3.680e-03, 5.719e-03, 5.857e-03, 1.458e-02, | 7.059e-03, 2.575e-02, 1.433e-02, 9.182e-03, | 1.343e-02, 6.247e-02, 2.000e-01, 3.710e-01, | 6.240e-01/) ! ! 0.5*(wave1+wave2) in centimeter ! rlmeuv=(/1.725e-05, 1.675e-05, 1.625e-05, 1.575e-05, | 1.525e-05, 1.475e-05, 1.425e-05, 1.375e-05, | 1.325e-05, 1.275e-05, 1.225e-05, 1.216e-05, | 1.175e-05, 1.125e-05, 1.075e-05, 1.038e-05, | 1.007e-05, 9.810e-06, 9.440e-06, 9.440e-06, | 9.440e-06, 8.555e-06, 8.555e-06, 8.555e-06, | 7.240e-06, 7.240e-06, 5.950e-06, 4.300e-06, | 3.050e-06, 2.570e-06, 1.895e-06, 1.125e-06, | 5.100e-07, 2.500e-07, 1.300e-07, 6.000e-08, | 2.250e-08/) ! ! O2 absorption coefficient: ! sigeuv(1,:) = (/ | 5.00e-01, 1.50e+00, 3.40e+00, 6.00e+00, 1.00e+01, | 1.30e+01, 1.50e+01, 1.20e+01, 2.20e+00, 4.00e-01, | 1.30e+01, 1.00e-02, 1.40e+00, 4.00e-01, 1.00e+00, | 1.15e+00, 1.63e+00, 1.87e+01, 3.25e+01, 1.44e+01, | 1.34e+01, 1.33e+01, 1.09e+01, 1.05e+01, 2.49e+01, | 2.36e+01, 2.70e+01, 2.03e+01, 1.68e+01, 1.32e+01, | 7.63e+00, 2.63e+00, 6.46e-01, 2.10e-01, 2.25e-01, | 3.40e-02, 4.54e-03/) ! ! O absorption coefficient: ! sigeuv(2,:) = (/ | 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, 3.79e+00, 4.10e+00, 3.00e+00, 4.79e+00, | 8.52e+00, 1.31e+01, 1.07e+01, 7.72e+00, 6.02e+00, | 3.78e+00, 1.32e+00, 3.25e-01, 1.05e-01, 1.13e-01, | 1.70e-02, 2.27e-03/) ! ! N2 absorption coefficient: ! sigeuv(3,:) = (/ | 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, 2.55e+00, 1.15e+02, 1.44e+01, | 2.18e+00, 7.17e+01, 1.31e+01, 2.14e+00, 5.45e+01, | 2.30e+01, 2.31e+01, 1.97e+01, 1.17e+01, 9.94e+00, | 5.09e+00, 1.53e+00, 3.46e-01, 1.14e+00, 1.41e-01, | 2.01e-02, 2.53e-03/) ! ! The three major species' ionization branching ratio (off absorption): ! ! O2 ! BPhotonI(1,:) = (/ | 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, 6.13e-01, 8.30e-01, 6.20e-01, 7.86e-01, | 7.56e-01, 5.34e-01, 5.74e-01, 5.49e-01, 4.76e-01, | 6.73e-01, 9.83e-01, 1.00e+00, 1.00e+00, 1.00e+00, | 1.00e+00, 1.00e+00, 1.00e+00, 1.00e+00, 1.00e+00, | 1.00e+00, 1.00e+00/) ! ! O ! BPhotonI(2,:) = (/ | 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.00e+00, 1.00e+00, 1.00e+00, 1.00e+00, | 1.00e+00, 1.00e+00, 1.00e+00, 1.00e+00, 1.00e+00, | 1.00e+00, 1.00e+00, 1.00e+00, 1.00e+00, 1.00e+00, | 1.00e+00, 1.00e+00/) ! ! N2 ! BPhotonI(3,:) = (/ | 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, 0.00e+00, 0.00e+00, 0.00e+00, 4.29e-01, | 6.80e-01, 1.00e+00, 1.00e+00, 1.00e+00, 1.00e+00, | 1.00e+00, 1.00e+00, 1.00e+00, 1.00e+00, 1.00e+00, | 1.00e+00, 1.00e+00/) ! ! photon ionization branching ratio for O+(2p),o+(2d),O+(4s) ! (off O photon ionization) ! ! O+(2p) ! brop2pPh(:) = (/ | 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, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, | 8.56e-03, 2.52e-01, 2.60e-01, 2.46e-01, 2.41e-01, | 2.33e-01, 2.27e-01, 2.26e-01, 2.24e-01, 2.24e-01, | 2.24e-01, 2.24e-01/) ! ! O+(2d) ! brop2dPh(:) = (/ | 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, 0.00e+00, 0.00e+00, 0.00e+00, 6.98e-02, | 3.37e-01, 4.51e-01, 4.24e-01, 4.03e-01, 4.02e-01, | 3.92e-01, 3.77e-01, 3.74e-01, 3.78e-01, 3.78e-01, | 3.78e-01, 3.78e-01/) ! ! O+(4s) ! brop4sPh(:) = (/ | 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.00e+00, 1.00e+00, 1.00e+00, 9.30e-01, | 6.55e-01, 2.98e-01, 3.17e-01, 3.46e-01, 3.50e-01, | 3.67e-01, 3.89e-01, 3.93e-01, 3.90e-01, 3.90e-01, | 3.90e-01, 3.90e-01/) ! ! O2 photon dissociation braching ratio ! bro2DPh(:) = (/ | 1.00e+00, 1.00e+00, 1.00e+00, 1.00e+00, 1.00e+00, | 1.00e+00, 1.00e+00, 1.00e+00, 1.00e+00, 1.00e+00, | 1.00e+00, 1.00e+00, 1.00e+00, 1.00e+00, 1.00e+00, | 1.00e+00, 3.87e-01, 1.70e-01, 3.80e-01, 2.14e-01, | 2.44e-01, 4.66e-01, 4.26e-01, 4.51e-01, 5.24e-01, | 3.27e-01, 1.74e-02, 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/) ! ! n2 photon dissociation braching ratio for n(2d) ! brn2DPh(:) = (/ | 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.00e+00, 1.00e+00, 1.00e+00, | 1.00e+00, 1.00e+00, 1.00e+00, 1.00e+00, 5.71e-01, | 3.20e-01, 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/) ! ! O2 photon dissociative ionization braching ratio ! bro2DIPh(:) = (/ | 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, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, | 5.35e-04, 1.08e-01, 2.40e-01, 3.51e-01, 3.76e-01, | 4.47e-01, 6.53e-01, 8.92e-01, 1.00e+00, 1.00e+00, | 1.00e+00, 1.00e+00/) ! ! n2 photon dissociative ionization braching ratio ! brn2DIPh(:) = (/ | 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, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, | 0.00e+00, 0.00e+00, 4.72e-03, 9.27e-02, 2.46e-01, | 2.53e-01, 2.49e-01, 2.82e-01, 9.60e-01, 9.60e-01, | 9.60e-01, 9.60e-01/) ! ! n(4s) photoionization cross section ! sigin4s(:) = (/ | 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, 3.24e+00, 2.48e+00, 2.11e+00, 1.12e+01, | 1.18e+01, 1.16e+01, 9.35e+00, 6.43e+00, 4.80e+00, | 2.55e+00, 6.81e-01, 1.66e-01, 5.68e-01, 7.05e-02, | 1.00e-02, 1.27e-03/) ! ! The three major species' electron impact ionization branching ratio ! off its photon ionization rate ! ! O2 ! BElectronI(1,:) = (/ | 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, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, | 0.00e+00, 0.00e+00, 2.38e-02, 1.05e-01, 2.42e-01, | 5.79e-01, 1.61e+00, 4.27e+00, 6.00e+01, 2.03e+01, | 5.02e+01, 2.11e+02/) ! ! O ! BElectronI(2,:) = (/ | 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, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, | 0.00e+00, 0.00e+00, 1.27e-01, 4.18e-01, 6.94e-01, | 1.09e+00, 2.19e+00, 4.99e+00, 7.14e+01, 2.36e+01, | 5.06e+01, 2.17e+02/) ! ! N2 ! BElectronI(3,:) = (/ | 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, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, | 0.00e+00, 0.00e+00, 3.07e-02, 1.78e-01, 3.61e-01, | 9.33e-01, 2.86e+00, 7.79e+00, 1.08e+01, 3.22e+01, | 8.09e+01, 3.43e+02/) ! ! electron impact ionization branching ratio for O+(2p),o+(2d),O+(4s) ! (off O photon ionization) ! ! O+(2p) ! brop2pEl(:) = (/ | 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, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, | 0.00e+00, 0.00e+00, 9.15e-03, 6.12e-02, 1.16e-01, | 2.03e-01, 4.36e-01, 1.01e+00, 1.46e+01, 4.77e+00, | 1.10e+01, 4.74e+01/) ! ! O+(2d) ! brop2dEl(:) = (/ | 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, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, | 0.00e+00, 0.00e+00, 3.39e-02, 1.48e-01, 2.53e-01, | 4.18e-01, 8.53e-01, 1.96e+00, 2.82e+01, 9.36e+00, | 2.07e+01, 8.85e+01/) ! ! O+(4s) ! brop4sEl(:) = (/ | 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, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, | 0.00e+00, 0.00e+00, 8.36e-02, 2.09e-01, 3.25e-01, | 4.70e-01, 9.02e-01, 2.02e+00, 2.86e+01, 9.42e+00, | 1.89e+01, 8.12e+01/) ! ! photoelectron dissociative ionization branching ratio of O2 ! bro2DIEl(:) = (/ | 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, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, | 0.00e+00, 0.00e+00, 9.60e-04, 1.43e-02, 5.23e-02, | 1.63e-01, 5.21e-01, 1.44e+00, 2.03e+01, 6.98e+00, | 1.79e+01, 7.61e+01/) ! ! photoelectron dissociative ionization branching ratio of N2 ! brn2DIEl(:) = (/ | 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, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, | 0.00e+00, 0.00e+00, 1.84e-04, 8.49e-03, 3.66e-02, | 1.46e-01, 5.71e-01, 1.65e+00, 2.29e+00, 6.95e+00, | 1.83e+01, 7.87e+01/) ! ! photoelectron dissociation branching ratio of N2 ! brn2DEl(:) = (/ | 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, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, | 0.00e+00, 0.00e+00, 1.57e-01, 5.15e-01, 7.64e-01, | 1.37e+00, 2.91e+00, 6.53e+00, 9.05e+00, 2.53e+01, | 5.21e+01, 2.45e+02/) ! ! photoelectron dissociation branching ratio of o2 ! bro2DEl(:) = (/ | 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, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, | 0.00e+00, 1.10e-02, 6.53e-01, 7.62e-01, 9.96e-01, | 1.27e+00, 2.04e+00, 4.11e+00, 5.70e+01, 1.78e+01, | 2.03e+01, 8.79e+01/) ! ! transfer units of cross section to cm^2 do m=1,3 do n=1,lmax sigeuv(m,n)=sigeuv(m,n)*1.E-18 enddo enddo do n = 1,lmax sigin4s(n) = sigin4s(n)*1.E-18 enddo ! do n = 1,l1 rlmsrc(n) = rlmeuv(n) sigsrc(n) = sigeuv(1,n) enddo ! euveff(:) = 0.05 quench = (/7.E-11,5.E-11,3.1401E-12,9.1E-3/) end subroutine init_qrj !--------------------------------------------------------------- subroutine ssflux (f107, f107a) ! ! Args: real,intent(in) :: f107,f107a ! ! Local: real :: pind integer :: l ! ! solar model is the same as EUVAC, i.e., solar flux at ! (f107,f107a) is scale as: sflux=sfmin(1+afac(P-80.)) ! where P=0.5(f107+f107a) pind=0.5*(f107+f107a) do l=1,lmax sflux(l)=sfmin(l)*(1+afac(l)*(pind-80.)) ! set solar flux to be 80% of the value when pind=80 ! if it becomes negative if (sflux(l) .le. 0.8*sfmin(l)) sflux(l) = 0.8*sfmin(l) enddo ! end subroutine ssflux !----------------------------------------------------------------------- subroutine alloc_q(lon0,lon1,lat0,lat1) ! ! Args: integer,intent(in) :: lon0,lon1,lat0,lat1 ! ! Local: integer :: istat ! allocate(rj(nlevp1,lon0:lon1,lat0:lat1),stat=istat) if (istat /= 0) write(6,"('>>> alloc_q: error allocating', | ' rj: stat=',i3)") istat allocate(qtef(nlevp1,lon0:lon1,lat0:lat1),stat=istat) if (istat /= 0) write(6,"('>>> alloc_q: error allocating', | ' qtef: stat=',i3)") istat allocate(qtotal(nlevp1,lon0:lon1,lat0:lat1),stat=istat) if (istat /= 0) write(6,"('>>> alloc_q: error allocating', | ' qtotal: stat=',i3)") istat allocate(qop2p(nlevp1,lon0:lon1,lat0:lat1),stat=istat) if (istat /= 0) write(6,"('>>> alloc_q: error allocating', | ' qop2p: stat=',i3)") istat allocate(qop2d(nlevp1,lon0:lon1,lat0:lat1),stat=istat) if (istat /= 0) write(6,"('>>> alloc_q: error allocating', | ' qop2d: stat=',i3)") istat allocate(qo2p(nlevp1,lon0:lon1,lat0:lat1),stat=istat) if (istat /= 0) write(6,"('>>> alloc_q: error allocating', | ' qo2p: stat=',i3)") istat allocate(qop(nlevp1,lon0:lon1,lat0:lat1),stat=istat) if (istat /= 0) write(6,"('>>> alloc_q: error allocating', | ' qop: stat=',i3)") istat allocate(qn2p(nlevp1,lon0:lon1,lat0:lat1),stat=istat) if (istat /= 0) write(6,"('>>> alloc_q: error allocating', | ' qn2p: stat=',i3)") istat allocate(qnp(nlevp1,lon0:lon1,lat0:lat1),stat=istat) if (istat /= 0) write(6,"('>>> alloc_q: error allocating', | ' qnp: stat=',i3)") istat allocate(qnop(nlevp1,lon0:lon1,lat0:lat1),stat=istat) if (istat /= 0) write(6,"('>>> alloc_q: error allocating', | ' qnop: stat=',i3)") istat end subroutine alloc_q !----------------------------------------------------------------------- end module qrj_module