! module aurora use params,only: nlev use input,only: f107 implicit none contains !----------------------------------------------------------------------- subroutine qaurora(fluxar,alp) use flds_atmos,only: rho,sht,xnn2,xno2,xno use flds_ionzrt,only: qti,qia ! ! Args: real,intent(in) :: fluxar,alp ! ! Local: integer :: k real :: delte = 35.e-3 real :: pi = 3.1415926 real :: theta,dencol,x,fy,denom,etan2,etao2,etao,dum ! theta=pi*0.5 do k=1,nlev dencol=rho(k)*sht(k) x=((dencol/4.e-6)**0.606)/alp fy=3.2333*x**2.56588*exp(-2.2541*x**0.7297198)+1.106907* | x**1.71349*exp(-1.8835444*x**0.86472135) qti(k)=fluxar*alp*fy/(delte*sht(k))+qti(k) denom=0.94*xnn2(k)+xno2(k)+0.55*xno(k) etan2=0.94*xnn2(k)*qti(k)/denom qia(1,k)=0.76*etan2+qia(1,k) qia(5,k)=0.24*etan2+qia(5,k) etao2=1.07*xno2(k)*etan2/xnn2(k) qia(2,k)= 0.67*etao2+qia(2,k) etao=0.59*xno(k)*etan2/xnn2(k) qia(3,k)=etao+0.33*etao2+qia(3,k) qia(4,k)=0. enddo ! k=1,nlev ! ! call fprint8('From qaurora:',nlev, ! | qia(1,:),qia(2,:),qia(3,:),qia(4,:),qia(5,:),dum,dum,dum, ! | 'qia(1)','qia(2)','qia(3)','qia(4)','qia(5)',' ',' ',' ') end subroutine qaurora !----------------------------------------------------------------------- subroutine cosray(phi) use flds_atmos,only: xno,xno2,xnn2 use flds_ionzrt,only: qcr ! ! Args: real,intent(in) :: phi ! ! Local: integer :: k real :: rad,x1,x2,x3,xm,a,b ! rad = 57.295 if (phi <= 53.) then x1=1.74e-18 x2=2.84e-17+(1.93e-17-2.84e-17)/135.*(f107-65.) x3=0.6+0.8*abs(cos(phi/rad)) do k=1,nlev ! was DO 4 xm=xno(k)+xno2(k)+xnn2(k) if (xm < 3.e+17) then qcr(k)=(x1+x2*(abs(sin(phi/rad)))**4.)*xm else qcr(k)=(x1+x2*(abs(sin(phi/rad)))**4.)*(3.e+17)**(1.-x3)* | (xm)**x3 endif enddo ! k=1,nlev else ! phi > 53 X1=1.44E-17 X2=4.92E-18 do k=1,nlev ! was DO 6 xm=xno(k)+xno2(k)+xnn2(k) a=x1*xm b=(x1+x2)*xm qcr(k)=b+(a-b)/135.*(f107-65.) enddo endif end subroutine cosray !----------------------------------------------------------------------- subroutine protonp use flds_atmos,only: rho,sht ! integer :: k real :: fperg,fpmev,emev,rang,rango,y,fy,qpe,qpi ! do k=1,nlev fperg=1. fpmev=fperg/1.602e-6 emev=10. rang=rho(k)*sht(k) rango=0.00271*(emev)**1.72 y=(rang/0.00271)**(1./1.72)/emev fy=0.12718*y**4.91199*exp(-1.8429*y**0.9936)+0.52472*y**1.556655 | *exp(-0.85732*y**1.4116) qpe=rho(k)*fpmev*fy/rang qpi=qpe*1.e+6/35. enddo end subroutine protonp !----------------------------------------------------------------------- subroutine proton(ut) C C THIS SUBROUTINE CALCULATES RATE OF PRODUCTION OF IONS BY AN INCOMING C FLUX OF PROTONS INCIDENT ISOTROPICALLY OVER THE UPPER HEMISPHERE. C INPUTS NEEDED ARE 6 GAUSSIAN QUADRATURE COEFFICIENTS, DETECTOR C ENERGY THRESHOLDS, AND A SERIES OF VALUES OF INTEGRAL DIRECTIONAL C FLUXES (PER CM2 SEC STERADIAN) OF PROTONS ABOVE EACH THRESHOLD. C THESE DATA ARE IN FILE "PRODATA". ENERGY THRESHOLDS IN THIS FILE C ARE THOSE OF THE PROTON DETECTORS ON THE NOAA SATELLITES. C THE CALCULATION IS CARRIED OUT FOR FIXED PRESSURE LEVELS, RATHER C THAN HEIGHTS, IN THE ATMOSPHERE. TO CONVERT TO HEIGHT CONSULT YOUR C FAVORITE MODEL ATMOSPHERE. THE TEMPERATURES CHOSEN ARE THOSE IN C HOUGHTON'S MODEL FOR SEPTEMBER AT 70 N LATITUDE. PRESSURE AND C TEMPERATURE ARE IN FILE "PRO2". C THE RATE OF ENERGY LOSS BY PROTONS IN MEV PER CM OF STANDARD C AIR IS APPROXIMATED BY A SERIES OF POWER-LAW RELATIONSHIPS OF C THE FORM DE/DX = K/E**N. THE ENERGY RANGE FROM 2150 MEV TO ZERO C IS BROKEN UP INTO 7 RANGES, EACH WITH DIFFERENT VALUES OF K AND N. C ! 11/04 btf: rewritten for new version of global mean model. ! ! Args: real,intent(in) :: ut ! ! Local: real :: qion(250),prion(250),flux(250) real :: press(nlev),temp(nlev) ! real :: x(6) = | (/.0337652,.1693953,.3806904,.6193096,.8306047,.9662348/) real :: a(6) = | (/.0856622,.1803808,.2339570,.2339570,.1803808,.0856622/) real :: amev(6) = (/.3,.8,2.5,16.,32.,80./) end subroutine proton !----------------------------------------------------------------------- end module aurora