#include "dims.h" SUBROUTINE ORORA use input_module,only: aurora use cons_module,only: kmax,kmaxp1,kmaxm1,jmax,len1,len2, | rmassinv,rmass,expz,grav,gask,dphi,pi,p0 implicit none C **** C **** CALCULATES AURORAL ADDITIONS TO IONIZATION RATES C **** ON ENTRY: C **** T3=DRIZL (PROFILE FOR AURORAL DRIZZLE) C **** T5=CUSP (PROFILE FOR AURORAL CUSP) C **** T6,T7=ALFA(1), ALFA(2) (CHARACTERISTIC PARTICLE ENERGIE C **** T8,T9=F(1),F(2) (PARTICLE FLUXES) C **** #include "params.h" #include "fcom.h" #include "vscr.h" #include "index.h" #include "buff.h" #include "phys.h" #include "ovalr.h" #include "aurlp.h" #include "crates_tdep.h" C **** C **** C **** ARRAYS AUREFF AND BDRIZ DEFINE VERTICAL PROFILES FOR C **** AURORAL HEATING EFFICIENCY AND BACKGROUND DRIZZLE C **** RESPECTIVELY C **** ! ! Local: real :: AUREFF(ZKMXP),BDRIZ(ZKMXP),userla,xn2 integer :: i,k,nmsk,nps1k,nps2k,ntk,m,nqpk,ishunk character(len=80) :: title ! ! Ion production from low energy proton source ! (see call low_proton below. Sub low_proton is in proton.f) real :: qia(zimxp,zkmxp,5) ! #if (NLEV==28) ! NLEV==28 for -7 to +7 by 0.50 AUREFF = |(/0.55,0.55,0.55,0.55,0.55,0.55,0.55,0.55,0.55,0.55,0.55,0.55, | 0.55,0.50,0.45,0.37,0.30,0.20,0.16,0.13,0.10,0.08,0.07,0.06, | 0.05,0.05,0.05,0.05,0.05/) BDRIZ = |(/2.17E-4,4.48E-3,3.78E-2,1.55E-1,3.26E-1,3.85E-1,3.50E-1,2.81E-1, | 2.04E-1,1.38E-1,8.74E-2,5.23E-2,3.01E-2,1.68E-2,9.27E-3,5.08E-3, | 2.79E-3,1.55E-3,8.62E-4,4.85E-4,2.76E-4,1.58E-4,9.11E-5,5.28E-5, | 3.08E-5,1.79E-5,1.04E-5,6.01E-6,3.49E-6/) #elif (NLEV==56) ! NLEV==56 for -7 to +7 by 0.25 AUREFF = |(/0.5500E+00, 0.5500E+00, 0.5500E+00, 0.5500E+00, 0.5500E+00, | 0.5500E+00, 0.5500E+00, 0.5500E+00, 0.5500E+00, 0.5500E+00, | 0.5500E+00, 0.5500E+00, 0.5500E+00, 0.5500E+00, 0.5500E+00, | 0.5500E+00, 0.5500E+00, 0.5500E+00, 0.5500E+00, 0.5500E+00, | 0.5500E+00, 0.5500E+00, 0.5500E+00, 0.5500E+00, 0.5500E+00, | 0.5250E+00, 0.5000E+00, 0.4750E+00, 0.4500E+00, 0.4100E+00, | 0.3700E+00, 0.3350E+00, 0.3000E+00, 0.2500E+00, 0.2000E+00, | 0.1800E+00, 0.1600E+00, 0.1450E+00, 0.1300E+00, 0.1150E+00, | 0.1000E+00, 0.9000E-01, 0.8000E-01, 0.7500E-01, 0.7000E-01, | 0.6500E-01, 0.6000E-01, 0.5500E-01, 0.5000E-01, 0.5000E-01, | 0.5000E-01, 0.5000E-01, 0.5000E-01, 0.5000E-01, 0.5000E-01, | 0.5000E-01, 0.5000E-01/) BDRIZ = |(/0.2170E-03, 0.2349E-02, 0.4480E-02, 0.2114E-01, 0.3780E-01, | 0.9640E-01, 0.1550E+00, 0.2405E+00, 0.3260E+00, 0.3555E+00, | 0.3850E+00, 0.3675E+00, 0.3500E+00, 0.3155E+00, 0.2810E+00, | 0.2425E+00, 0.2040E+00, 0.1710E+00, 0.1380E+00, 0.1127E+00, | 0.8740E-01, 0.6985E-01, 0.5230E-01, 0.4120E-01, 0.3010E-01, | 0.2345E-01, 0.1680E-01, 0.1303E-01, 0.9270E-02, 0.7175E-02, | 0.5080E-02, 0.3935E-02, 0.2790E-02, 0.2170E-02, 0.1550E-02, | 0.1206E-02, 0.8620E-03, 0.6735E-03, 0.4850E-03, 0.3805E-03, | 0.2760E-03, 0.2170E-03, 0.1580E-03, 0.1246E-03, 0.9110E-04, | 0.7195E-04, 0.5280E-04, 0.4180E-04, 0.3080E-04, 0.2435E-04, | 0.1790E-04, 0.1415E-04, 0.1040E-04, 0.8205E-05, 0.6010E-05, | 0.4750E-05, 0.3490E-05/) #endif ! call addfsech('QO2P_PRE',' ',' ',f(1,nqo2p),zimxp,zkmxp,zkmxp,j) ! call addfsech('QOP_PRE' ,' ',' ',f(1,nqop) ,zimxp,zkmxp,zkmxp,j) ! call addfsech('QN2P_PRE',' ',' ',f(1,nqn2p),zimxp,zkmxp,zkmxp,j) C **** C **** NO CALCULATION IF ISHUNK=0 C **** ! userla and ishunk are local (s.a. heelis): ! ishunk > 0 means j = 1-12 and 25-36 ! aion will be called 24x4 = 96 times per time step ! USERLA=(FLOAT(J-JMAX)-.5)*dphi+pi/2. ISHUNK=1 IF (aurora.EQ.0) ISHUNK=0 IF (ABS(USERLA).LT.pi/6.) ISHUNK=0 IF (ISHUNK.EQ.0) then qia = 1.e-20 qteaur = 1.e-20 ! ! Put qia on secondary histories (see also below, when ishunk > 0): ! ! call addfsech('QIA1_N2P', ! | 'N2+ production from low-energy proton source', ! | ' ',qia(:,:,1),zimxp,zkmxp,zkmx,j) ! call addfsech('QIA1_O2P', ! | 'O2+ production from low-energy proton source', ! | ' ',qia(:,:,2),zimxp,zkmxp,zkmx,j) ! call addfsech('QIA1_OP', ! | 'O+ production from low-energy proton source', ! | ' ',qia(:,:,3),zimxp,zkmxp,zkmx,j) ! call addfsech('QIA1_NOP', ! | 'NO+ production from low-energy proton source', ! | ' ',qia(:,:,4),zimxp,zkmxp,zkmx,j) ! call addfsech('QIA1_NP', ! | 'N+ production from low-energy proton source', ! | ' ',qia(:,:,5),zimxp,zkmxp,zkmx,j) return endif ! write(6,"('aurora_ions: lat=',i2,' alfa1=',/,(6e12.4))") j,t6 ! write(6,"('aurora_ions: lat=',i2,' alfa2=',/,(6e12.4))") j,t7 ! write(6,"('aurora_ions: lat=',i2,' cusp =',/,(6e12.4))") j,t5 ! write(6,"('aurora_ions: lat=',i2,' drizl=',/,(6e12.4))") j,t3 C **** C **** CALCULATE TOTAL IONIZATION, QI C **** DO 1 K=1,KMAX DO 1 I=1,LEN1 C **** C **** T3=DRIZL (PROFILE FOR AURORAL DRIZZLE) C **** T5=CUSP (PROFILE FOR AURORAL CUSP) C **** T6,T7=ALFA(1), ALFA(2) (CHARACTERISTIC PARTICLE ENERGIE C **** T8,T9=F(1),F(2) (PARTICLE FLUXES) C C **** S7,S8=X(1),X(2) (AURORA) C **** S9=X (CUSP) C **** S10=X (DRIZZLE) C **** WHERE X=(N/4.E-6)**0.606/ALFA C **** AND N=P0*EXP(-Z)/G C S6(I,K)=(p0*expz(K)/(grav*4.E-6))**0.606 S7(I,K)=S6(I,K)/T6(I) ! xalfa1 = p0zp/alfa1 S8(I,K)=S6(I,K)/T7(I) ! xalfa2 = p0zp/alfa2 S9(I,K)=S6(I,K)/ALFAC ! xcusp = p0zp/alfac S10(I,K)=S6(I,K)/ALFAD ! xdrizl = p0zp/alfad 1 CONTINUE C **** C **** S3,S4 = F(X(1)),F(X(2)) (AURORA) C **** S5 = F(X) (CUSP) C **** S6 = F(X) (DRIZZLE) C **** ! 9/99: aion is a cpu consumer because of **real ! and exp. It is located in inline.f. ! CALL AION(S7,S3,LEN2) ! (xalfa1,alfa1_ion) CALL AION(S8,S4,LEN2) ! (xalfa2,xalfa2_ion) CALL AION(S9,S5,LEN2) ! (cusp,cusp_ion) CALL AION(S10,S6,LEN2) ! (drizl,drizl_ion) call addfsech('ALFA1',' ',' ',s3 ,zimxp,zkmxp,zkmx,j) call addfsech('ALFA2',' ',' ',s4 ,zimxp,zkmxp,zkmx,j) call addfsech('CUSP ',' ',' ',s5 ,zimxp,zkmxp,zkmx,j) call addfsech('DRIZL',' ',' ',s6,zimxp,zkmxp,zkmx,j) ! ! Contribution of low energy protons: ! qia(zimxp,zkmxp,5) is output ion production for N2+, O2+, O+, NO+, N+ ! (top level zkmxp not defined) (s.a., aurht.f) ! ! call low_proton(alfa_lp,eflux_lp,qia,j) qia = 0. ! ! Put qia on secondary histories: ! subroutine addfsech(name,long_name,units,f2d,londim,levdim, ! | nlev,lat) ! ! call addfsech('QIA1_N2P', ! | 'N2+ production from low-energy proton source', ! | ' ',qia(:,:,1),zimxp,zkmxp,zkmx,j) ! call addfsech('QIA1_O2P', ! | 'O2+ production from low-energy proton source', ! | ' ',qia(:,:,2),zimxp,zkmxp,zkmx,j) ! call addfsech('QIA1_OP', ! | 'O+ production from low-energy proton source', ! | ' ',qia(:,:,3),zimxp,zkmxp,zkmx,j) ! call addfsech('QIA1_NOP', ! | 'NO+ production from low-energy proton source', ! | ' ',qia(:,:,4),zimxp,zkmxp,zkmx,j) ! call addfsech('QIA1_NP', ! | 'N+ production from low-energy proton source', ! | ' ',qia(:,:,5),zimxp,zkmxp,zkmx,j) ! DO 4 K=1,KMAX DO 4 I=1,LEN1 C **** C **** S7,S8=ALFA(1)*F(1),ALFA(2)*F(2) (AURORA) C **** S9=ALFA*F (CUSP) C **** S10=ALFA*F (DRIZZLE) C **** S7(I,K)=T6(I)*T8(I) S8(I,K)=T7(I)*T9(I) S9(I,K)=T5(I)*ALFAC*FC S10(I,K)=T3(I)*ALFAD*FD 4 CONTINUE NMSK=NJ+NMS NPS1K=NJ+NPS NPS2K=NJ+NPS2 NTK=NJ+NT DO 2 I=1,LEN2 C **** C **** S1=SUM(QI)=QT C **** =SUM(F*ALFA*FX(X)/(35.E-3*H)) C **** S11(I,1)=grav*.5*(F(I,NMSK)+F(I,NMSK+1))/(35.E-3*gask*F(I,NTK)) ! qsum = xalfa1 * alfa1_ion + xalfa2 * alfa2_ion S1(I,1) = S7(I,1) * S3(I,1) + S8(I,1) * S4(I,1) ! qsum = qsum + fcusp * cusp_ion S1(I,1) = S1(I,1) + S9(I,1) * S5(I,1) ! qsum = (qsum + fdrizl * drizl_ion) * barm_t S1(I,1) = (S1(I,1) + S10(I,1) * S6(I,1)) * S11(I,1) C **** C **** S3=DENOM=0.92*PS3/RMASS(3)+1.5*PS1/RMASS(1)+0.56*PS2/ C **** RMASS(2) C **** S3(I,1)=0.92*(1.-F(I,NPS1K)-F(I,NPS2K))*rmassinv(3)+1.5* | F(I,NPS1K)*rmassinv(1)+0.56*F(I,NPS2K)*rmassinv(2) 2 CONTINUE ! write(6,"('aurora: lat=',i3,' falfa1=',/,(6e12.4))") j,s7(:,1) ! write(6,"('aurora: lat=',i3,' falfa2=',/,(6e12.4))") j,s8(:,1) ! write(6,"('aurora: lat=',i3,' fcusp=',/,(6e12.4))") j,s9(:,1) ! write(6,"('aurora: lat=',i3,' fdrizl=',/,(6e12.4))") j,s10(:,1) C **** C **** ADD IONIZATION CONTRIBUTIONS TO NQO2P, NQOP, NQN2P, NQNP C **** C **** S13 = QO2P, S14 = QOP, S15 = QN2P (K+1/2) c c 6/98: Include production due to low energy protons. c qia(zimxp,zkmxp,5) is output ion production from low_proton c (called above) for N2+, O2+, O+, NO+, N+ respectively c (currently, NO+ is not defined) ! 2/28/00: added to tgcm14 c NPS1K=NJ+NPS NPS2K=NPS1K+KMAXP1 DO 40 I=1,LEN2 S13(I,1)=S1(I,1)*F(I,NPS1K)/(RMASS(1)*S3(I,1))+qia(i,1,2) S14(I,1)=S1(I,1)*(0.5*F(I,NPS1K)*rmassinv(1)+ 1 0.56*F(I,NPS2K)*rmassinv(2))/S3(I,1)+qia(i,1,3) xn2 = (1.-F(I,NPS1K)-F(I,NPS2K)) if (xn2 <= 1.e-8) xn2 = 1.e-8 S15(I,1)=S1(I,1)*0.7*xn2/(RMASS(3)*S3(I,1))+qia(i,1,1) 40 CONTINUE call addfsech('BARMT' ,' ',' ',s11,zimxp,zkmxp,zkmx,j) call addfsech('QSUM' ,' ',' ',s1 ,zimxp,zkmxp,zkmx,j) call addfsech('DENOM' ,' ',' ',s3 ,zimxp,zkmxp,zkmx,j) call addfsech('QO2P_AUR',' ',' ',s13,zimxp,zkmxp,zkmx,j) call addfsech('QOP_AUR' ,' ',' ',s14,zimxp,zkmxp,zkmx,j) call addfsech('QN2P_AUR',' ',' ',s15,zimxp,zkmxp,zkmx,j) call addfsech('QO2P_PRE',' ',' ',f(1,nqo2p),zimxp,zkmxp,zkmxp,j) call addfsech('QOP_PRE' ,' ',' ',f(1,nqop) ,zimxp,zkmxp,zkmxp,j) call addfsech('QN2P_PRE',' ',' ',f(1,nqn2p),zimxp,zkmxp,zkmxp,j) C **** SET NQO2P, NQOP, NQN2P DO 35 M=1,3 K=(3-M)*KMAXP1+1 NQPK=NQO2P+(M-1)*KMAXP1 C **** LEVELS 2 THRU KMAX DO 36 I=LEN1+1,LEN2 F(I,NQPK)=F(I,NQPK)+SQRT(S15(I,K)*S15(I,K-1)) 36 CONTINUE C **** LEVELS 1 AND KMAXP1 DO 37 I=1,LEN1 F(I,NQPK) = F(I,NQPK)+1.5*S15(I,K)-0.5*S15(I,K+1) F(I,NQPK+KMAX) = F(I,NQPK+KMAX)+1.5*S15(I,K+KMAXM1)-0.5* 1 S15(I,K+KMAX-2) 37 CONTINUE 35 CONTINUE C **** SET NQNP NQPK=NQNP C **** LEVELS 2 THRU KMAX DO 38 I=LEN1+1, LEN2 F(I,NQPK)=F(I,NQPK)+.22/.7*SQRT(S15(I,1)*S15(I,0)) 38 CONTINUE C **** LEVELS 1 AND KMAXP1 DO 39 I=1,LEN1 F(I,NQPK) = F(I,NQPK)+.22/.7*(1.5*S15(I,1)-0.5*S15(I,2)) F(I,NQPK+KMAX) = F(I,NQPK+KMAX)+.22/.7*(1.5*S15(I,KMAX) 1 -0.5*S15(I,KMAXM1)) 39 CONTINUE ! ! disn2p (see crates_tdep.h, qrj.f, qtieff.f, xray.f) DO 41 K=2,KMAX DO 41 I=1,LEN1 DISN2P(I,K)=DISN2P(I,K)+SQRT(S15(I,K-1)*S15(I,K)) 41 CONTINUE DO 42 I=1,LEN1 DISN2P(I,1) = DISN2P(I,1)+1.5*S15(I,1)-0.5*S15(I,2) DISN2P(I,KMAXP1) = DISN2P(I,KMAXP1)+1.5*S15(I,KMAX)-0.5* 1 S15(I,KMAX-1) 42 CONTINUE call addfsech('QO2P' ,' ',' ',f(1,nqo2p) ,zimxp,zkmxp,zkmxp,j) call addfsech('QOP' ,' ',' ',f(1,nqop) ,zimxp,zkmxp,zkmxp,j) call addfsech('QN2P' ,' ',' ',f(1,nqn2p) ,zimxp,zkmxp,zkmxp,j) call addfsech('QNP' ,' ',' ',f(1,nqnp) ,zimxp,zkmxp,zkmxp,j) call addfsech('DISN2P',' ',' ',disn2p,zimxp,zkmxp,zkmxp,j) RETURN END C