      SUBROUTINE SETTEI
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"
      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)
      common/aurlp/ alfa_lp(zimxp),eflux_lp(zimxp),qteaur(zimxp)
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) = 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)*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*EXP(-4605.2/F(I,NTEK)),5.71E-8*EXP(-3352.6
     1  /F(I,NTEK)),F(I,NTEK)-1000.>=0.)
        S9(I,1) = merge(S9(I,1),2.53E-6*S5(I,1)*EXP(-17620./F(I,NTEK)),
     1  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.-EXP(S10(I,1)))
     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)))
     1  /S10(I,1)
        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)*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
