#include "dims.h" SUBROUTINE ORORA use cons_module,only: len1,len2,kmax,kmaxp1,kmaxm1,jmax,pi,dphi, | p0,expz,boltz,rmass,avo,gask,grav,rmassinv_o2,rmassinv_o, | rmassinv_n2 use input_module,only: aurora 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 "crates_tdep.h" #include "ovalr.h" real :: alfa3,flux3,xn2 COMMON /AURXTRA/ ALFA3(ZIMXP),FLUX3(ZIMXP) !$OMP THREADPRIVATE (/aurxtra/) !DIR$ TASKCOMMON AURXTRA ! ! Local: integer :: i,k,ishunk,lem,m real :: userla,zlem integer :: ntk,nps1k,nps2k,nmsk,nqpk C **** C **** NO CALCULATION IF ISHUNK=0 C **** 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) return C **** C **** CALCULATE TOTAL IONIZATION, QI C **** NTK = NJ+NT-1 NPS1K = NJ+NPS-1 NPS2K = NJ+NPS2-1 NMSK = NJ+NMS-1 DO 1 K=1,KMAX NTK = NTK+1 NPS1K = NPS1K+1 NPS2K = NPS2K+1 NMSK = NMSK+1 DO I=1,LEN1 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 **** S1(I,K) = p0*expz(K)/(boltz*F(I,NTK)*(F(I,NPS1K)*rmassinv_o2 1 +F(I,NPS2K)*rmassinv_o+(1.-F(I,NPS1K)-F(I,NPS2K))* 2 rmassinv_n2))/avo S2(I,K) = gask*F(I,NTK)/(grav*F(I,NMSK))*S1(I,K) S12(I,K) = ((S2(I,K)/0.00271)**0.58140)/ALFAD2 S6(I,K)=(p0*expz(K)/(grav*4.E-6))**0.606 S7(I,K)=S6(I,K)/T6(I) S8(I,K)=S6(I,K)/T7(I) S9(I,K)=S6(I,K)/ALFAC S13(I,K)=S6(I,K)/ALFA3(I) S10(I,K)=S6(I,K)/ALFAD enddo 1 CONTINUE C **** C **** S3,S4 = F(X(1)),F(X(2)) (AURORA) C **** S5 = F(X) (CUSP) C **** S6 = F(X) (DRIZZLE) C **** DO 901 I=1,LEN2 S3(I,1) = 1.E-20 S4(I,1) = 1.E-20 S5(I,1) = 1.E-20 S14(I,1) = 1.E-20 S15(I,1) = 1.E-20 S6(I,1) = 1.E-20 901 CONTINUE ! ! AION calculates auroral electrons (ECR and Ed Su) ! BION calculates high-energy protons (D. Lummerzeim) ! Low-energy protons will be added in a later version (Marina Galand) ! Use zsb and dz to calculate index to zlem for any vertical resolution. ! zlem = -7.5 lem = int((zlem-zsb)/dz)+1 ! write(6,"('orora: j=',i2,' zlem=',f6.2,' lem=',i2)") j,zlem,lem CALL AION(S7,S3,LEM,LEN2) zlem = -5.0 lem = int((zlem-zsb)/dz)+1 ! write(6,"('orora: j=',i2,' zlem=',f6.2,' lem=',i2)") j,zlem,lem CALL AION(S9,S5,LEM,LEN2) zlem = -7.5 lem = int((zlem-zsb)/dz)+1 CALL AION(S10,S6,LEM,LEN2) LEM = 1 IF(FD2.GT.1.E-19) CALL BION(S12,S14,LEM,LEN2) IF(E30.GT.1.E-19) CALL AION(S13,S15,LEM,LEN2) DO 4 K=1,KMAX DO 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 S12(I,K) = T3(I)*S1(I,K)*FD2*1.E+6/(S2(I,K)*35.) S13(I,K) = ALFA3(I)*FLUX3(I) S10(I,K)=T3(I)*ALFAD*FD enddo 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)) S1(I,1)=S7(I,1)*S3(I,1)+S8(I,1)*S4(I,1) !falfa1*alfa1_ion+falfa2*alfa2_ion S1(I,1)=S1(I,1)+S13(I,1)*S15(I,1) !falfa3*alfa3_ion S1(I,1)=S1(I,1)+S9(I,1)*S5(I,1) !fcusp*cusp_ion ! qsum = (qsum +fdrizl *drizl_ion)*barm_t+xalfa_sp*alfasp_bion S1(I,1)=(S1(I,1)+S10(I,1)*S6(I,1))*S11(I,1)+S12(I,1)*S14(I,1) C **** C **** S3=DENOM=0.92*PS3/RMASS(3)+1.5*PS1/RMASS(1)+0.56*PS2/ C **** RMASS(2) S3(I,1)=0.92*(1.-F(I,NPS1K)-F(I,NPS2K))*rmassinv_n2+ | 1.5*F(I,NPS1K)*rmassinv_o2+0.56*F(I,NPS2K)*rmassinv_o 2 CONTINUE C **** C **** ADD IONIZATION CONTRIBUTIONS TO NQO2P, NQOP, NQN2P, NQNP C **** C **** S13 = QO2P, S14 = QOP, S15 = QN2P (K+1/2) ! Phys fields: ! 37 NQO2P 25981 KMAXP1 S13 ! 38 NQOP 26026 KMAXP1 S14 ! 39 NQN2P 26071 KMAXP1 S15 ! 40 NQNOP 26116 KMAXP1 (not defined here) ! 41 NQNP 26161 KMAXP1 ! 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)) S14(I,1)=S1(I,1)*(0.5*F(I,NPS1K)*rmassinv_o2+ 1 0.56*F(I,NPS2K)*rmassinv_o)/S3(I,1) ! S15(I,1)=S1(I,1)*0.7*merge(1.-F(I,NPS1K)-F(I,NPS2K),1.E-6, ! + (1.-F(I,NPS1K)-F(I,NPS2K))-1.E-6>=0.)/(RMASS(3)*S3(I,1)) ! xn2 = (1.-f(i,nps1k)-f(i,nps2k)) if (xn2 <= 1.e-6) xn2 = 1.e-6 s15(i,1) = s1(i,1)*0.7*xn2/(rmass(3)*s3(i,1)) 40 CONTINUE 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 ! ! 11/30/98: DISN2P added as per modsrc.kibo: 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 RETURN END