#include "dims.h" SUBROUTINE SETTEI use input_module,only: f107 use init_module,only: sin_sundec,cos_sundec use cons_module,only: imax,imaxp2,kmax,kmaxm1,kmaxp1,expz,len1, | len2,len3,rmassinv,boltz,p0,gask,avo,grav,dphi,pi,dlamda, | expzmid,expzmid_inv,evergs,dipmin C **** C **** CALCULATE ELECTRON AND ION TEMPERATURES C **** use init_module,only: secs implicit none #include "params.h" #include "fcom.h" #include "vscr.h" #include "index.h" #include "buff.h" #include "phys.h" #include "aurlp.h" #include "trgm.h" ! ! 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,rtd,root2,f107te,rlat, | coslat,sinlat,expsk C C **** SET FPOLAR (POLAR TE FLUX) C C DATA FPOLAR/-3.0E+9/ DATA FPOLAR/-3.0E+7/ 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 RADIANS TO DEGREES C RTD=180./pi C C **** SQRT(2.) C ROOT2 = SQRT(2.) C C **** CALCULATE HEAT FLUX, FE, AT UPPER BOUNDARY F107TE = f107 IF (F107TE.GT.235.) F107TE = 235. C C **** CALCULATE SOLAR ZENITH ANGLE, CHI C RLAT=-.5*pi+(FLOAT(J-1)+.5)*dphi 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)*dlamda+pi)*12./pi,24.) C C **** T2 = CHI C T2(I)=ACOS(sin_sundec*SINLAT+cos_sundec*COSLAT* | COS(pi*(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 if (ABS(RLATM(I,J))-pi/4.5>=0.) then t3(i) = 1. else t3(i) = .5*(1.+SIN(pi*(ABS(RLATM(I,J))-pi/9.) | /(pi/4.5))) endif 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 t4(i) = 0. C t5(i) = 0. C C **** T1 = FE C if (T2(I)-.5*pi>=0.) then t1(i) = t5(i) else t1(i) = t4(i) endif if ((T2(I)*RTD-80.)*(T2(I)*RTD-100.)>=0.) then t1(i) = t1(i)*evergs else t1(i) = (.5*(T4(I)+T5(I))+.5*(T4(I)-T5(I))*COS(pi* | (T2(I)*RTD-80.)/20.))*evergs endif C C **** ADD FPOLAR IF MAGNETIC LATITUDE .GE. 60. DEG if (ABS(RLATM(I,J))-pi/3.>=0.) t1(i) = T1(I)+FPOLAR*evergs 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) if (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 = expzmid_inv*expz(K) DO 6 I = 1,LEN1 S9(I,K) = EXPSK 6 CONTINUE C C **** LEVEL KMAXP1 C EXPSK = expzmid*expz(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) = p0*S9(I,1)*F(I,NMSK)/(boltz*S10(I,1)) S13(I,1) = S8(I,1)*S13(I,1)*rmassinv(1) S12(I,1) = S8(I,1)*S12(I,1)*rmassinv(2) S11(I,1) = S8(I,1)*S11(I,1)*rmassinv(3) C C **** S5 = SQRT(TE) C if (s14(i,1)-50. < 0.) s14(i,1) = 50. 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)))*evergs 8 CONTINUE C C **** S7 = H (K+1/2) C **** S6 = H (K) C DO 9 I=1,LEN2 S7(I,1) = gask*F(I,NTK)/(.5*(F(I,NMSK)+F(I,NMSK+1))*grav) 9 CONTINUE DO 10 I = 1,LEN3 S6(I,1) = gask*S10(I,1)/(F(I,NMSK)*grav) 10 CONTINUE C C **** S1 = P, S3 = R, S2 = -(P+R), S4 = 0. (K+1/2) C C **** T2 = (SIN(DIPMAG))**2 ! dipmin is in trgm.h, initialized in types.f: C DO 11 I = 1,LEN1 if (abs(DIPMAG(I,J)) >= dipmin) then T2(I) = (SIN(DIPMAG(I,J)))**2 else T2(I) = (SIN(DIPMIN))**2 endif if (t2(i)-.10 < 0.) t2(i) = .10 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)*dz**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)*dz*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) if (s11(i,1)-1.e-20 < 0.) s11(i,1) = 1.e-20 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) if (s10(i,1)-1.e+4 < 0.) s10(i,1) = 1.e+4 S15(I,1) = SQRT(S10(I,1)) 18 CONTINUE C C **** S10 = EXP(-Z) (K+1/2) C DO 19 K = 1,KMAX EXPSK = expz(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) = p0*S10(I,1)*.5*(F(I,NMSK)+F(I,NMSK+1))/ 1 (boltz*F(I,NTK)) C C **** S14 = N(O2), S13 = N(O) (K+1/2) C S14(I,1) = S6(I,1)*F(I,NPS1K)*rmassinv(1) S13(I,1) = S6(I,1)*F(I,NPS2K)*rmassinv(2) C C **** S12 = N(N2) (K+1/2) C S12(I,1) = 1.-F(I,NPS1K)-F(I,NPS2K) if (s12(i,1) < 0.) s12(i,1) = 0. S12(I,1) = S6(I,1)*S12(I,1)*rmassinv(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) = exp(-((((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)*evergs 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 if (f(i,ntek)-1000. >= 0.) then S9(I,1) = 2.E-7*exp(-4605.2/F(I,NTEK)) else S9(I,1) = 5.71E-8*exp(-3352.6/F(I,NTEK)) endif if (2000.-f(i,ntek) < 0.) | s9(i,1) = 2.53E-6*S5(I,1)*exp(-17620./F(I,NTEK)) 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.-exp(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.-exp(S10(I,1)))/S10(I,1) ! ! Gpi runs may stop here with overflow in exp() when TE gets too hot. ! Sometimes S10 may already be -INF from above. ! S10(I,1) = 1.57E-12*exp((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)*evergs C C **** S11 = L0(E,N) C S11(I,1) = (S11(I,1)+S13(I,1)*S10(I,1))*S15(I,1)*evergs 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))*evergs 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))*evergs 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 if (f(i,ntek)-f(i,ntik) >= 0.) then S7(I,1)=S10(I,1)*(F(I,NTEK)-F(I,NTIK)) else s7(i,1) = 0. endif S7(I,1) = (S8(I,1)*(F(I,NTEK)-F(I,NTK))+S7(I,1)) 1 *avo/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)/avo 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 if (f(i,ntek)-f(i,ntk) < 0.) f(i,ntek) = f(i,ntk) ! if (f(i,ntek) > 5000.) f(i,ntek) = 5000. 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=NUI NVIK=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)) if (f(i,ntik)-f(i,ntk) < 0.) f(i,ntik) = f(i,ntk) ! if (f(i,ntik) > 5000.) f(i,ntik) = 5000. 30 CONTINUE RETURN END C