SUBROUTINE SETTEI CDIR$ VFUNCTION EXPHF implicit none C **** C **** CALCULATE ELECTRON AND ION TEMPERATURES C **** include "params.h" include "blnk.h" include "vscr.h" include "cons.h" include "index.h" include "strt.h" include "buff.h" include "phys.h" real :: rlatm,rlonm,dipmag,decmag,sndec,csdec,sn2dec,sncsdc, | dumdum COMMON/TRGM/RLATM(ZIMXP,ZJMX),RLONM(ZIMXP,ZJMX), 1 DIPMAG(ZIMXP,ZJMX),DECMAG(ZIMXP,ZJMX),SNDEC(ZIMXP,ZJMX), 2 CSDEC(ZIMXP,ZJMX),SN2DEC(ZIMXP,ZJMX),SNCSDC(ZIMXP,ZJMX), 3 DUMDUM(ZJMX,3) real :: alfa_lp,eflux_lp,qteaur common/aurlp/ alfa_lp(zimxp),eflux_lp(zimxp),qteaur(zimxp) ! ! Local: integer :: i,ntek,nek,nps1k,nps2k,ntk,k,nmsk,nqpk,ntik,nopk,no2pk, | nnopk,n,nqk,nuik,nvik,nuk,nvk,nlxxk,nlxyk,nlyxk,nlyyk real :: fpolar,del,alam,ad,sd,evterg,rtd,root2,f107te,rlat, | coslat,sinlat,expsk,exphf C C **** SET FPOLAR (POLAR TE FLUX) C DATA FPOLAR/-3.0E+9/ C C **** SMALL NUMBER DELTA C DATA DEL/1.E-6/ C C **** CORRECTION FACTOR FOR NEUTRAL HEATING DUE TO L(E,O1D) C DATA ALAM,AD,SD/0.0069,0.0091,2.3E-11/ C C **** CONVERSION FROM EV TO ERGS C DATA EVTERG/1.602E-12/ C C **** CONVERSION RADIANS TO DEGREES C RTD=180./C(110) C C **** SQRT(2.) C ROOT2 = SQRT(2.) C C **** CALCULATE HEAT FLUX, FE, AT UPPER BOUNDARY F107TE = C(61) IF (F107TE.GT.235.) F107TE = 235. C C **** CALCULATE SOLAR ZENITH ANGLE, CHI C RLAT=-.5*C(110)*FLOAT(1-NHEMI)+(FLOAT(J-1)+.5)*C(2) COSLAT=COS(RLAT) SINLAT=SIN(RLAT) DO 1 I = 1,LEN1 T2(I)=FLOAT(I-3) 1 CONTINUE DO 2 I=1,LEN1 C C **** T2 = LOCAL TIME C T2(I)=AMOD(SECS/3600.+(T2(I)*C(1)+C(110))*12./C(110),24.) C C **** T2 = CHI C T2(I)=ACOS(C(95)*SINLAT+C(96)*COSLAT*COS(C(110)*(T2(I)-12.) 1 /12.)) C C **** T3 = A = .5*(1.+SIN(PI*(ABS(RLATM)-PI/6.)/(PI/3.))) C **** FOR ABS(RLATM).LT.PI/3. C **** A = 1. FOR ABS(RLATM).GE.PI/3 C T3(I) = merge(1.,.5*(1.+SIN(C(110)*(ABS(RLATM(I,J))-C(110)/9.) 1 /(C(110)/4.5))),ABS(RLATM(I,J))-C(110)/4.5>=0.) C C **** T4 = FED T5 = FEN C C C **** INCREASED HEAT FLUX FOR TE FROM PROTONOSPHERE C C T4(I) =( -5.0E+7*F107TE*T3(I)-4.0E+7*F107TE)/2.5 T4(I) =( -5.0E+7*F107TE*T3(I)-4.0E+7*F107TE)*1.2 C T5(I) = T4(I)/3. T5(I) = T4(I)/2. T4(I) = T4(I)+QTEAUR(I) T5(I) = T5(I)+QTEAUR(I) C C **** T1 = FE C T1(I) = merge(T5(I),T4(I),T2(I)-.5*C(110)>=0.) T1(I) = merge(T1(I),.5*(T4(I)+T5(I))+.5*(T4(I)-T5(I)) 1 *COS(C(110)*(T2(I)*RTD-80.)/20.),(T2(I)*RTD-80.) 2 *(T2(I)*RTD-100.)>=0.)*EVTERG C C **** ADD FPOLAR IF MAGNETIC LATITUDE .GE. 60. DEG T1(I) = merge(T1(I)+FPOLAR*EVTERG,T1(I),ABS(RLATM(I,J)) 1 -C(110)/3.>=0.) 2 CONTINUE C C **** CALCULATE K0 = KE/TE**2.5 (K) C NTEK = NJ+NTE NEK = NJ+NE NPS1K = NJ+NPS NPS2K = NJ+NPS2 NTK = NJ+NT C C **** CALCULATE TE, PSI(O2), PSI(O), PSI(N2), T (K) C C **** LEVELS 2 THRU KMAX C DO 3 I = LEN1+1,LEN2 S14(I,1) = .5*(F(I,NTEK)+F(I,NTEK-1)) S13(I,1) = .5*(F(I,NPS1K)+F(I,NPS1K-1)) S12(I,1) = .5*(F(I,NPS2K)+F(I,NPS2K-1)) S10(I,1) = .5*(F(I,NTK)+F(I,NTK-1)) 3 CONTINUE C C **** LEVELS 1 AND KMAXP1 C DO 4 I = 1,LEN1 S14(I,1) = 1.5*F(I,NTEK)-.5*F(I,NTEK+1) S13(I,1) = 1.5*F(I,NPS1K)-.5*F(I,NPS1K+1) S12(I,1) = 1.5*F(I,NPS2K)-.5*F(I,NPS2K+1) S10(I,1) = 1.5*F(I,NTK)-.5*F(I,NTK+1) S14(I,KMAXP1) = 1.5*F(I,NTEK+KMAXM1)-.5*F(I,NTEK+KMAX-2) S13(I,KMAXP1) = 1.5*F(I,NPS1K+KMAXM1)-.5*F(I,NPS1K+KMAX-2) S12(I,KMAXP1) = 1.5*F(I,NPS2K+KMAXM1)-.5*F(I,NPS2K+KMAX-2) S10(I,KMAXP1) = 1.5*F(I,NTK+KMAXM1)-.5*F(I,NTK+KMAX-2) 4 CONTINUE C C **** S11 = PSI(N2) C DO 5 I = 1,LEN3 S11(I,1) = 1.-S12(I,1)-S13(I,1) S11(I,1) = merge(S11(I,1),0.,S11(I,1)>=0.) 5 CONTINUE C C **** S9 = EXP(-Z) (K) C C **** LEVELS 1 THRU KMAX C DO 6 K = 1,KMAX EXPSK = C(87)*EXPS(K) DO 6 I = 1,LEN1 S9(I,K) = EXPSK 6 CONTINUE C C **** LEVEL KMAXP1 C EXPSK = C(86)*EXPS(KMAX) DO 7 I = 1,LEN1 S9(I,KMAXP1) = EXPSK 7 CONTINUE C C **** C **** S13 = N(O2), S12 = N(O), S11 = N(N2) (K) C NMSK = NJ+NMS DO 8 I = 1,LEN3 C C **** S8 = P0*EXP(-Z)*MBAR/(K*T) (K) C S8(I,1) = C(81)*S9(I,1)*F(I,NMSK)/(C(84)*S10(I,1)) S13(I,1) = S8(I,1)*S13(I,1)/RMASS(1) S12(I,1) = S8(I,1)*S12(I,1)/RMASS(2) S11(I,1) = S8(I,1)*S11(I,1)/RMASS(3) C C **** S5 = SQRT(TE) C S14(I,1) = merge(S14(I,1),50.,S14(I,1)-50.>=0.) S5(I,1) = SQRT(S14(I,1)) C C **** S15 = K0 C S15(I,1) = 7.5E5/(1.+3.22E4*S14(I,1)**2/F(I,NEK)*(S5(I,1) 1 *(2.82E-17-3.41E-21*S14(I,1))*S11(I,1)+(2.2E-16+7.92E-18*S5(I,1) 2 )*S13(I,1)+1.1E-16*(1.+5.7E-4*S14(I,1))*S12(I,1)))*EVTERG 8 CONTINUE C C **** S7 = H (K+1/2) C **** S6 = H (K) C DO 9 I=1,LEN2 S7(I,1) = C(57)*F(I,NTK)/(.5*(F(I,NMSK)+F(I,NMSK+1))*C(54)) 9 CONTINUE DO 10 I = 1,LEN3 S6(I,1) = C(57)*S10(I,1)/(F(I,NMSK)*C(54)) 10 CONTINUE C C **** S1 = P, S3 = R, S2 = -(P+R), S4 = 0. (K+1/2) C C **** T2 = (SIN(DIPMAG))**2 C DO 11 I = 1,LEN1 T2(I) = (SIN(DIPMAG(I,J)))**2 C T2(I) = merge(T2(I),.01,T2(I)-.01>=0.) T2(I) = merge(T2(I),.10,T2(I)-.10>=0.) 11 CONTINUE C C **** S2 = T2 C DO 12 K = 1,KMAX DO 12 I = 1,LEN1 S2(I,K) = T2(I) 12 CONTINUE DO 13 I = 1,LEN2 C C **** S1 = 2./7.*(SIN(I))**2/(H*DZ**2) (K) C S1(I,1) = 2./7.*S2(I,1)/(S7(I,1)*C(3)**2) S3(I,1) = S1(I,1)*S15(I,2)/S6(I,2) S1(I,1) = S1(I,1)*S15(I,1)/S6(I,1) S2(I,1) = -(S1(I,1)+S3(I,1)) S4(I,1) = 0. 13 CONTINUE C C **** BOUNDARY CONDITIONS C DO 14 I = 1,LEN1 C C **** LOWER BOUNDARY C S2(I,1) = S2(I,1)-S1(I,1) S4(I,1) = S4(I,1)-2.*S1(I,1)*S10(I,1)**3.5 S1(I,1) = 0. C C **** UPPER BOUNDARY C S2(I,KMAX) = S2(I,KMAX)+S3(I,KMAX) S4(I,KMAX) = S4(I,KMAX)+S3(I,KMAX)*C(3)*3.5*S6(I,KMAXP1)*T1(I) 1 /S15(I,KMAXP1) S3(I,KMAX) = 0. 14 CONTINUE C C **** S15 = NE, S14 = N(O2), S13 = N(O), S12 = N(N2) (K+1/2) C **** S11 = QI(TOT) (K+1/2) C NEK = NJ+NE NPS1K = NJ+NPS NPS2K = NJ+NPS2 NQPK = NQO2P NTK = NJ+NT NMSK = NJ+NMS NTEK = NJ+NTE NTIK = NJ+NTI NOPK = NJ+NOP NO2PK = NJ+NO2P NNOPK = NNOP DO 15 I = 1,LEN3 S11(I,1)=0. 15 CONTINUE C C **** S11 = QI(TOT), TOTAL IONIZATION RATE (K+1/2) C DO 16 N = 1,7 DO 17 I = 1,LEN3 S11(I,1) = S11(I,1)+F(I,NQPK) S11(I,1) =merge(S11(I,1),1.E-20,S11(I,1)-1.E-20>=0.) 17 CONTINUE NQPK = NQPK+KMAXP1 16 CONTINUE DO 32 I = 1,LEN2 S11(I,1) = SQRT(S11(I,1)*S11(I,2)) 32 CONTINUE DO 31 I = 1,LEN1 S11(I,KMAXP1) = 0. 31 CONTINUE C C **** NOW SET NE (K+1/2) C DO 18 I = 1,LEN2 S10(I,1) = F(I,NEK)*F(I,NEK+1) S10(I,1) =merge(S10(I,1),1.E+4,S10(I,1)-1.E+4>=0.) S15(I,1) = SQRT(S10(I,1)) 18 CONTINUE C C **** S10 = EXP(-Z) (K+1/2) C DO 19 K = 1,KMAX EXPSK = EXPS(K) DO 19 I = 1,LEN1 S10(I,K) = EXPSK 19 CONTINUE DO 20 I = 1,LEN2 C C **** S6 = P0*EXP(-Z)*MBAR/(K*T) (K+1/2) C S6(I,1) = C(81)*S10(I,1)*.5*(F(I,NMSK)+F(I,NMSK+1))/ 1 (C(84)*F(I,NTK)) C C **** S14 = N(O2), S13 = N(O) (K+1/2) C S14(I,1) = S6(I,1)*F(I,NPS1K)/RMASS(1) S13(I,1) = S6(I,1)*F(I,NPS2K)/RMASS(2) C C **** S12 = N(N2) (K+1/2) C S12(I,1) = 1.-F(I,NPS1K)-F(I,NPS2K) S12(I,1) = merge(S12(I,1),0.,S12(I,1)>=0.) S12(I,1) = S6(I,1)*S12(I,1)/RMASS(3) C C **** CALCULATE SOURCE TERM, QE (K+1/2) C C **** S10 = R (K+1/2) C S10(I,1) = ALOG(S15(I,1)/(S14(I,1)+S12(I,1)+0.1*S13(I,1))) C C **** S10 = E (EFFICIENCY) (K+1/2) C S10(I,1) = exphf(-((((0.001996*S10(I,1)+0.08034)*S10(I,1)+1.166) C C CORRECTION FACTOR OF 2 INCREASE IN TE HEATING RATE 1 *S10(I,1)+6.941)*S10(I,1)+12.75))*1.0 C C **** S4 = S4 - QE (K+1/2) C S4(I,1) = S4(I,1)-S10(I,1)*S11(I,1)*EVTERG C C **** S11 = L0(E,N) = L(E,N)/(TE-TN) (K+1/2) C **** S8 = L0*(E,N) = L*(E,N)/(TE-TN) (K+1/2) C C **** S5 = SQRT(TE) (K+1/2) C S5(I,1) = SQRT(F(I,NTEK)) C C **** S9 = A(E,N2,VIB) (K+1/2) C S9(I,1) = merge(2.E-7*exphf(-4605.2/F(I,NTEK)),5.71E-8* | exphf(-3352.6/F(I,NTEK)),F(I,NTEK)-1000.>=0.) S9(I,1) = merge(S9(I,1),2.53E-6*S5(I,1)* | exphf(-17620./F(I,NTEK)),2000.-F(I,NTEK)>=0.) C C **** S10 = (1.-EXP(3200.*(1./TE-1./TN))/(TE-TN) (K+1/2) C **** (SAFEGUARDED FOR ZERO (TE-TN)) C S10(I,1) = 3200.*(1./F(I,NTEK)-1./F(I,NTK)) S10(I,1) = SIGN(ABS(S10(I,1))+DEL,S10(I,1)) S10(I,1) = -3200./(F(I,NTEK)*F(I,NTK))* | (1.-exphf(S10(I,1)))/S10(I,1) C C **** S10 = L0(E,N2,VIB)/(NE*N(N2)) (K+1/2) C S10(I,1) = 1.3E-4*S10(I,1)*S9(I,1) 20 CONTINUE DO 500 I=1,LEN2 C C **** S11 = (L0(E,N2)+L0(E,N2,ROT)+L0(E,N2,VIB))/NE C S11(I,1) = S12(I,1)*(1.77E-19*(1.-1.21E-4*F(I,NTEK))*F(I,NTEK) 1 +2.9E-14/S5(I,1)+S10(I,1)) C C **** S11 = S11+(L0(E,O2)+L0(E,O2,ROT)+L0(E,O2,VIB)/NE C S11(I,1) = S11(I,1)+S14(I,1)*(1.21E-18*(1.+3.6E-2*S5(I,1)) 1 *S5(I,1)+6.9E-14/S5(I,1)+3.125E-21*F(I,NTEK)**2) C C **** S10 = L0(E,O,1D)/(NE*N(O)) C S10(I,1) = 22713.*(1./F(I,NTEK)-1./F(I,NTK)) S10(I,1) = SIGN(ABS(S10(I,1))+DEL,S10(I,1)) S10(I,1) = 22713./(F(I,NTEK)*F(I,NTK))* | (1.-exphf(S10(I,1)))/S10(I,1) S10(I,1) = 1.57E-12*exphf((2.4E4+0.3*(F(I,NTEK)-1500.)-1.947E-5 1 *(F(I,NTEK)-1500.)*(F(I,NTEK)-4000.))*(F(I,NTEK)-3000.)/(3000. 2 *F(I,NTEK)))*S10(I,1) C C **** S11 = S11+(L0(E,O)+L0(E,O,F))/NE C S11(I,1) = S11(I,1)+S13(I,1)*(7.9E-19*(1.+5.7E-4*F(I,NTEK)) 1 *S5(I,1)+3.4E-12*(1.-7.E-5*F(I,NTEK))/F(I,NTK)*(150./F(I,NTEK) 2 +0.4)) C C **** S8 = L0*(E,N) C S8(I,1) = (S11(I,1)+S13(I,1)*(1.-ALAM/(AD+SD*S12(I,1))) 1 *S10(I,1))*S15(I,1)*EVTERG C C **** S11 = L0(E,N) C S11(I,1) = (S11(I,1)+S13(I,1)*S10(I,1))*S15(I,1)*EVTERG C C **** CALCULATE L0(E,I) = L(E,I)/(TE-TI), WHERE L(E,I) IS LOSS C **** DUE TO INTERACTIONS BETWEEN ELECTRONS AND IONS C S10(I,1) = 3.2E-8*S15(I,1)/(S5(I,1)*F(I,NTEK))*15. 1 *(F(I,NOPK)+0.5*F(I,NO2PK)+0.53*F(I,NNOPK))*EVTERG C C S5 = SQRT(TN) (K+1/2) C S5(I,1) = SQRT(F(I,NTK)) C C **** S9 = L0(I,N) =L(I,N)/(TI-TN) C S9(I,1) = ((6.6E-14*S12(I,1)+5.8E-14*S14(I,1)+0.21E-14*S13(I,1) 1 *ROOT2*S5(I,1))*F(I,NOPK)+(5.45E-14*S14(I,1)+5.9E-14*S12(I,1)+ 2 4.5E-14*S13(I,1))*F(I,NNOPK)+(5.8E-14*S12(I,1)+4.4E-14*S13(I,1) 3 +0.14E-14*S14(I,1)*S5(I,1))*F(I,NO2PK))*EVTERG C C **** COMPLETE TRIDIAGONAL MATRIX COEFFICIENTS AND RHS C C **** S2 = S2-(L0(E,N)+L0(E,I))/TE**2.5 = Q C S2(I,1) = S2(I,1)-(S11(I,1)+S10(I,1))/F(I,NTEK)**2.5 C C **** S4 = S4-L0(E,N)*TN-L0(E,I)*TI C S4(I,1) = S4(I,1)-S11(I,1)*F(I,NTK)-S10(I,1)*F(I,NTIK) 500 CONTINUE C C **** CALCULATE HEATING DUE TO ELECTRON/NEUTRAL AND C **** ELECTRON/ION COLLISIONS (K) C C **** UNITS ARE: ERGS/SEC/GM C C **** S7 = HEATING IN ERGS/SEC/GM (K+1/2) C DO 24 I=1,LEN2 S7(I,1)=merge(S10(I,1)*(F(I,NTEK)-F(I,NTIK)),0.,F(I,NTEK) 1 -F(I,NTIK)>=0.) S7(I,1) = (S8(I,1)*(F(I,NTEK)-F(I,NTK))+S7(I,1)) 1 *C(85)/S6(I,1) 24 CONTINUE C C **** ADD COLLISIONAL HEATING TO NQ FOR USE IN THERMODYNAMIC C **** EQUATION C NQK = NQ C C **** LEVELS 2 THRU KMAX C DO 27 I=1,LEN2-LEN1 F(I,NQK+1) = F(I,NQK+1)+.5*(S7(I,1)+S7(I,2)) 27 CONTINUE C C **** LEVELS 1 AND KMAXP1 C DO 28 I=1,LEN1 F(I,NQK) = F(I,NQK)+1.5*S7(I,1)-0.5*S7(I,2) F(I,NQK+KMAX) = F(I,NQK+KMAX)+1.5*S7(I,KMAX)-0.5*S7(I,KMAXM1) 28 CONTINUE C C **** S15 = RHO C DO 29 I = 1,LEN2 S15(I,1) = S6(I,1)/C(85) 29 CONTINUE C C **** SOLVE TRIDIAGONAL SYSTEM, S8 = G C NTEK = NJNP+NTE CALL TRSOLV(S1,S2,S3,S4,S7,F(1,NTEK),LEN1,3,LEN1-2,KMAXP1,1,KMAX 1,1) C C **** PERIODIC POINTS C DO 21 I = 1,2 DO 21 K = 0,KMAXM1 F(I,NTEK+K) = F(I+IMAX,NTEK+K) F(I+IMAXP2,NTEK+K) = F(I+2,NTEK+K) 21 CONTINUE NTK = NJ+NT NTEK =NJNP+NTE DO 33 I = 1,LEN2 F(I,NTEK) = merge(F(I,NTEK),F(I,NTK),F(I,NTEK)-F(I,NTK)>=0.) 33 CONTINUE C C **** TE = TE**(2./7.) C DO 23 K = 0,KMAXM1 DO 23 I = 1,LEN1 F(I,NTEK+K) = F(I,NTEK+K)**(2./7.) 23 CONTINUE C C **** S7 = QJI (ION JOULE HEATING) (K+1/2) C NUIK=NDJ+NUI NVIK=NDJ+NVI NUK=NJ+NU NVK=NJ+NV NLXXK=NLXX NLXYK=NLXY NLYXK=NLYX NLYYK=NLYY NTIK = NJNP+NTI NTEK = NJNP+NTE NTK = NJ+NT DO 30 I=1,LEN2 C C **** S6 = UI(K+1/2) C **** S5 = VI(K+1/2) C S6(I,1)=.5*(F(I,NUIK)+F(I,NUIK+1)) S5(I,1)=.5*(F(I,NVIK)+F(I,NVIK+1)) C C **** S7 = QJI C S7(I,1)=.5*((F(I,NLXXK)+F(I,NLXXK+1))*S6(I,1)*(S6(I,1)-F(I,NUK)) 1 +(F(I,NLXYK)+F(I,NLXYK+1))*S6(I,1)*(S5(I,1)-F(I,NVK))- 2 (F(I,NLYXK)+F(I,NLYXK+1))*S5(I,1)*(S6(I,1)-F(I,NUK))+(F(I,NLYYK) 3 +F(I,NLYYK+1))*S5(I,1)*(S5(I,1)-F(I,NVK))) C C **** SET TI C F(I,NTIK) = (S7(I,1)*S15(I,1)+S10(I,1)*F(I,NTEK)+S9(I,1)* 1 F(I,NTK))/(S10(I,1)+S9(I,1)) F(I,NTIK) = merge(F(I,NTIK),F(I,NTK),F(I,NTIK)-F(I,NTK)>=0.) 30 CONTINUE RETURN END C