CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CC             3D Theoretical Ionospheric Model of Earth in IGGCAS             CC
CC                                        Developed by Zhipeng Ren             CC
CC                                        Email: zpren@mail.iggcas.ac.cn       CC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
      SUBROUTINE TIME3Dmodule(TimeStep,Z,Tn,cO,cO2,cN2,cH,cHe,cN,Ws,We,
     &                        Wu,Qop,QHep,OpL,Qep,QHp,HepL,HpL,cO2p,cNOp
     &                        ,cN2p,ed1,ed2,Op,Te,Ti,Vi,Ped,Hall,zigm11,
     &                        zigm22,zigm2,zigmc,rim1,rim2,t)
      use geogrid_t3d  ,only: nlev,nlat,nlon
      use maggrid      ,only: nmlat,nmlon,nmlonp1
      INCLUDE 'EIG.h'                         !Hindex:,BY,BZ,SWS,HPN,HPS
      REAL*8 Tn(nlat,nlon,NLev),cO(nlat,nlon,NLev),cO2(nlat,nlon,NLev),
     &       cN2(nlat,nlon,NLev),cH(nlat,nlon,NLev),cHe(nlat,nlon,NLev),
     &       cN(nlat,nlon,NLev),Ws(nlat,nlon,NLev),We(nlat,nlon,NLev),Wu
     &       (nlat,nlon,NLev),Qop(nlat,nlon,NLev),QHep(nlat,nlon,NLev),
     &       OpL(nlat,nlon,NLev),Qep(nlat,nlon,NLev),cO2p(nlat,nlon,NLev
     &       ),cNOp(nlat,nlon,NLev),cN2p(nlat,nlon,NLev),Z(nlat,nlon,
     &       NLev),QHp(Nlat,Nlon,Nlev),HepL(Nlat,Nlon,Nlev),HpL(Nlat,
     &       Nlon,Nlev),ed1(nmlat,nmlonp1),ed2(nmlat,nmlonp1),
     &       Op(nlat,nlon,NLev),Te(nlat,nlon,NLev),Ti(nlat,nlon,NLev),Vi
     &       (nlat,nlon,NLev),zigm11(nmlat,nmlonp1),zigm22(nmlat,nmlonp1
     &       ),zigm2(nmlat,nmlonp1),zigmc(nmlat,nmlonp1),rim1(nmlat,
     &       nmlonp1),rim2(nmlat,nmlonp1),Ped(nlat,nlon,Nlev),Hall(nlat,
     &       nlon,Nlev)

        CALL GCMpre2alt(Z,Tn,cO,cO2,cN2,cH,cHe,cN,Ws,We,Wu,Qop,QHep,OpL,
     &                  Qep,cO2p,cNOp,cN2p,QHp,HepL,HpL)
        CALL GCM2TIME3D()
        CALL Dynamo2TIME3D(ed1,ed2)
  	CALL GridInterpolation(TimeStep)
        CALL ImodelDensityEquation(TimeStep)
        CALL ImodelEnergyEquation(TimeStep)
	CALL TIME3D2Dynamo(zigm11,zigm22,zigm2,zigmc,rim1,rim2)
        CALL TIME3D2GCM(Z,Op,Te,Ti,Vi,Ped,Hall)
      END SUBROUTINE TIME3Dmodule
c---------------------------------------------------------------------
      SUBROUTINE ImodelDensityEquation(TimeStep)
      INCLUDE 'EIG.h'
	COMMON/IGrid/ d1i(lmi,NFi,NLi,3),d2i(lmi,NFi,NLi,3),glati(lmi,NFi,
     &              NLi),gloni(lmi,NFi,NLi),Vsi(lmi,NFi,NLi),Wsi(lmi,NFi
     &              ,NLi),Bsi(lmi,NFi,NLi,4),gxi(lmi,NFi,NLi,3),dvv(lmi,
     &              NFi,NLi,2),mlov(NLi,NFi,2),d1ih(NHN,2,NLi,3),d2ih(
     &              NHN,2,NLi,3),glatih(NHN,2,NLi),glonih(NHN,2,NLi),
     &              Vsih(NHN,2,NLi),Wsih(NHN,2,NLi),Bsih(NHN,2,NLi,4),
     &              gxih(NHN,2,NLi,3),dvvh(NHN,2,NLi,2),mlovh(NHN,2,2),
     &              MH300(NHN+NLi)
	COMMON/EEmC/  EEm(NLi,NFi,2),EEH(NHN,2,2)
	COMMON/REIM3/ REIM(lmi,NFi,NLi,12),REIMold(lmi,NFi,NLi,8),REIH(NHN
     &             ,2,NLi,12),REIHold(NHN,2,NLi,8)
	COMMON/REIM3N/ PN(lmi,NFi,NLi,8),PLI(lmi,NFi,NLi,8),PNH(NHN,2,NLi,
     &               8),PLIH(NHN,2,NLi,8)
      COMMON/consti/ Bolt,Sim(3),Sif(3)
      REAL*8 DELTA(lmi),A(lmi),B(lmi),C(lmi),D(lmi),A1(lmi),A2(lmi),Ni(3
     &    ),Cij(3,3),Beta(3),Cin(5),DifC(lmi,3,9),DifWB(lmi,11),cHT(4,6)

c mid- and low-latitude
      DO k=1,NFi
        REIM(1,k,1,1)=dmax1(1.0D-25,(REIM(1,k,1,1)+PLI(1,k,1,1)*TimeStep
     &                  )/(PLI(1,k,1,4)*TimeStep+1.0D0))
        IF(iHe.eq.1)REIM(1,k,1,2)=dmax1(1.0D-25,(REIM(1,k,1,2)+PLI(1,k,1
     &                  ,2)*TimeStep)/(PLI(1,k,1,5)*TimeStep+1.0D0))
        REIM(1,k,1,3)=dmax1(1.0D-25,(REIM(1,k,1,3)+PLI(1,k,1,3)*TimeStep
     &                  )/(PLI(1,k,1,6)*TimeStep+1.0D0))
        REIM(1,k,1,7)=sum(REIM(1,k,1,1:6))
      ENDDO
      REIM(lmi,1:NFi,1,1:7)=REIM(1,1:NFi,1,1:7)
	isd=1
	IF(iHe.ne.1)isd=2
	A(1)=0.0D0
	B(1)=1.0D0
	C(1)=0.0D0
	xt1=1.0D0/TimeStep
	DO k=2,NLi
	  ik=2*k-1
	  l=lmi-k+1
	  A(ik)=0.0D0
	  B(ik)=1.0D0
	  C(ik)=0.0D0
	  DO j=1,NFi
	    DO i=1,k
	      Ni(1:3)=REIM(k,j,i,1:3)
            CALL Thdff(Ni,REIM(k,j,i,9),xt1,Cij,Beta) !Thdff(Ni,Ti,Cij,Betap)
c                CollisionFrequency(Ti,Tn,cO,cO2,cN2,cN,cH,cHe,Cin)
	      CALL CollisionFrequency(REIM(k,j,i,9),PN(k,j,i,1),PN(k,j,i,2
     &                              ),PN(k,j,i,3),PN(k,j,i,4),PN(k,j,i,7
     &                              ),PN(k,j,i,6),PN(k,j,i,5),Cin)
	      BsBo=Bsi(k,j,i,4)*bolt
		  xtte=REIM(k,j,i,8)/REIM(k,j,i,7)
	      gtg=gxi(k,j,i,3)+sum(EEm(k,j,1:2)*gxi(k,j,i,1:2))
	      DvvT=EEm(k,j,1)*dvv(k,j,i,1)+EEm(k,j,2)*dvv(k,j,i,2)
	      DifWB(i,1)=Bsi(k,j,i,4)*Wsi(k,j,i)
	      DO is=1,3,isd
	        Xi=sim(is)*(sum(Cij(is,1:3))+Cin(is)+0.004)                  !Xi
	        DifC(i,is,3)=BsBo/Xi
			DifC(i,is,8)=Xi/sim(is)              !Bs*k/Xi of Si_2
           DifC(i,is,1)=DifC(i,is,3)*(REIM(k,j,i,9)+xtte*REIM(k,j,i,is))    !Ri
	        DifC(i,is,2)=DifC(i,is,3)*xtte       !Bs*k*Te/Ne/Xi: part of Si_2
			DifC(i,is,4)=DifC(i,is,3)*Beta(is)   !Bs*k*Beta/Xi: part of Si_2
	        DifC(i,is,5)=sim(is)*(gtg+sum(REIM(k,j,i,10:12)*Cij(is,1:3
     &                     ))+PN(k,j,i,8)*(Cin(is)))/Xi                !Si_1/Xi
	        DifC(i,is,6)=REIMold(k,j,i,is)/timestep+PLI(k,j,i,is)
	        DifC(i,is,7)=PLI(k,j,i,is+3)
              IF(DvvT.GE.0.0D0)DifC(i,is,7)=DifC(i,is,7)+DvvT
              IF(DvvT.LT.0.0D0)DifC(i,is,6)=DifC(i,is,6)-DvvT*REIMold(k,
     &                                      j,i,is)
	        DifC(i,is,9)=REIMold(k,j,i,is+5)
	      ENDDO
	    if(i.ne.k)then
	      ij=ik-i+1
		  Ni(1:3)=REIM(l,j,i,1:3)
            CALL Thdff(Ni,REIM(l,j,i,9),xt1,Cij,Beta)
	      CALL CollisionFrequency(REIM(l,j,i,9),PN(l,j,i,1),PN(l,j,i,2
     &                              ),PN(l,j,i,3),PN(l,j,i,4),PN(l,j,i,7
     &                              ),PN(l,j,i,6),PN(l,j,i,5),Cin)
	      BsBo=Bsi(l,j,i,4)*bolt
		  xtte=REIM(l,j,i,8)/REIM(l,j,i,7)
	      gtg=gxi(l,j,i,3)+sum(EEm(k,j,1:2)*gxi(l,j,i,1:2))
	      DvvT=EEm(k,j,1)*dvv(l,j,i,1)+EEm(k,j,2)*dvv(l,j,i,2)
	      DifWB(ij,1)=Bsi(l,j,i,4)*Wsi(l,j,i)
	      DO is=1,3,isd
	        Xi=sim(is)*(sum(Cij(is,1:3))+Cin(is)+0.004)
	        DifC(ij,is,3)=BsBo/Xi
			DifC(ij,is,8)=Xi/sim(is)
              DifC(ij,is,1)=DifC(ij,is,3)*(REIM(l,j,i,9)+xtte*REIM(l,j,i
     &                      ,is))
	        DifC(ij,is,2)=DifC(ij,is,3)*xtte
			DifC(ij,is,4)=DifC(ij,is,3)*Beta(is)
	        DifC(ij,is,5)=sim(is)*(gtg+sum(REIM(l,j,i,10:12)*Cij(is,1:
     &                      3))+PN(l,j,i,8)*(Cin(is)))/Xi
	        DifC(ij,is,6)=REIMold(l,j,i,is)/timestep+PLI(l,j,i,is)
	        DifC(ij,is,7)=PLI(l,j,i,is+3)
              IF(DvvT.GE.0.0D0)DifC(ij,is,7)=DifC(ij,is,7)+DvvT
              IF(DvvT.LT.0.0D0)DifC(ij,is,6)=DifC(ij,is,6)-DvvT*REIMold(
     &                                       l,j,i,is)
	        DifC(ij,is,9)=REIMold(l,j,i,is+5)
	      ENDDO
	    endif
          ENDDO
	    DO i=1,k-1
	      DifWB(i,3)=1.0D0/(Vsi(k,j,i+1)-Vsi(k,j,i))
            DifWB(i,2)=(DifWB(i,1)+DifWB(i+1,1))*0.5D0
            DifWB(i,4:5)=(REIM(k,j,i+1,8:9)-REIM(k,j,i,8:9))*DifWB(i,3)
	      DifWB(i,6)=2.0D0/(Wsi(k,j,i+1)+Wsi(k,j,i))
	      ij=ik-i
		  DifWB(ij,3)=1.0D0/(Vsi(l,j,i)-Vsi(l,j,i+1))
            DifWB(ij,2)=(DifWB(ij,1)+DifWB(ij+1,1))*0.5D0
           DifWB(ij,4:5)=(REIM(l,j,i,8:9)-REIM(l,j,i+1,8:9))*DifWB(ij,3)
	      DifWB(ij,6)=2.0D0/(Wsi(l,j,i+1)+Wsi(l,j,i))
	    ENDDO
	    DO is=1,3,isd
	      DO i=1,k-1
              dtei=DifWB(i,4)+DifWB(i,5)
			dti=DifWB(i,5)
              dei=((REIM(k,j,i+1,7)-REIM(k,j,i+1,is))-(REIM(k,j,i,7)-
     &           REIM(k,j,i,is)))
	        xx=dabs(dei)/dmax1(1.0D-30,dabs(REIM(k,j,i,is)-REIM(k,j,
     &           i+1,is)))
	        dei=dei*DifWB(i,3)*dexp(Sif(is)*xx)       !if(xx.ge.1.0D0)
              DifWB(i,7)=DifWB(i,2)*(DifC(i,is,1)+DifC(i+1,is,1))*0.5D0
     &                   *DifWB(i,3)                          ! B*W*Ri/dv
!!!!! B*W*Si/2: Si= (Si_1 - Bs*k*Te/Ne*dei - Bs*k*dtei - Bs*k*Beta*dti)/Xi
              DifWB(i,8)=DifWB(i,2)*((DifC(i,is,5)+DifC(i+1,is,5))-(DifC
     &            (i,is,2)+DifC(i+1,is,2))*dei-(DifC(i,is,3)+DifC(i+1,is
     &            ,3))*dtei-(DifC(i,is,4)+DifC(i+1,is,4))*dti)*0.25D0
	        ij=ik-i
			dtei=DifWB(ij,4)+DifWB(ij,5)
			dti=DifWB(ij,5)
              dei=((REIM(l,j,i,7)-REIM(l,j,i,is))-(REIM(l,j,i+1,7)-
     &            REIM(l,j,i+1,is)))
	        xx=dabs(dei)/dmax1(1.0D-30,dabs(REIM(l,j,i,is)-REIM(l,j,
     &           i+1,is)))
	        dei=dei*DifWB(ij,3)*dexp(Sif(is)*xx)      !if(xx.ge.1.0D0)
              DifWB(ij,7)=DifWB(ij,2)*(DifC(ij,is,1)+DifC(ij+1,is,1))*
     &                    0.5D0*DifWB(ij,3)
              DifWB(ij,8)=DifWB(ij,2)*((DifC(ij,is,5)+DifC(ij+1,is,5))-(
     &         DifC(ij,is,2)+DifC(ij+1,is,2))*dei-(DifC(ij,is,3)+DifC(ij
     &        +1,is,3))*dtei-(DifC(ij,is,4)+DifC(ij+1,is,4))*dti)*0.25D0
		  ENDDO
            D(1)=dmax1(1.0D-25,(REIM(k,j,1,is)+PLI(k,j,1,is)*TimeStep)/(
     &            PLI(k,j,1,is+3)*TimeStep+1.0D0))
            D(ik)=dmax1(1.0D-25,(REIM(l,j,1,is)+PLI(l,j,1,is)*TimeStep)/
     &            (PLI(l,j,1,is+3)*TimeStep+1.0D0))
	      DifWB(1:ik-1,9)=1.0D0/DifWB(1:ik-1,3)    !dv
	      IF(kVC.eq.1)THEN
	        DifWB(1:ik-1,10)=1.0D0
	        DifWB(1:ik-1,11)=0.0D0
	      ELSE
	        DifWB(1:ik-1,10)=Dexp(0.0D0-(DifC(1:ik-1,is,8)+DifC(2:ik,
     &                 is,8))*0.5D0*TimeStep)
	        DifWB(1:ik-1,11)=0.25D0*(DifC(1:ik-1,is,9)+DifC(2:ik,is,9)
     &                 )*DifWB(1:ik-1,2)*DifWB(1:ik-1,10) ! B*W*exp(-vt)*V_old/2
             DifWB(1:ik-1,10)=1.0D0-DifWB(1:ik-1,10)
	      ENDIF
            DO i=2,ik-1
	        WDV=(DifWB(i-1,6)+DifWB(i,6))/(DifWB(i-1,9)+DifWB(i,9))
              X1=DifWB(i-1,7)
              X2=DifWB(i-1,8)
              X3=DifWB(i-1,10)
              X4=DifWB(i-1,11)
              X5=DifWB(i,7)
              X6=DifWB(i,8)
              X7=DifWB(i,10)
              X8=DifWB(i,11)
              xx=DifC(i,is,7)
			A(i)=(-X4-X3*(X2+X1))*WDV
              C(i)=(X8-X7*(X5-X6))*WDV
			B(i)=1.0D0/timestep+xx+(X3*(X1-X2)-X4+X8+X7*(X6+X5))*WDV
              D(i)=DifC(i,is,6)
	      ENDDO
            CALL TRIDAGI(DELTA,lmi,1,ik,A,B,C,D)
            DELTA(1:ik)=dmax1(1.0D-25,DELTA(1:ik))
            DO i=1,ik-1
	        XX=(DELTA(i)+DELTA(i+1))*0.5D0
c	        A1(i)=(DifWB(i,8)*2.0D0-DifWB(i,7)*(DELTA(i+1)-DELTA(i))/
c     &              XX)/DifWB(i,2)
	        A1(i)=(DifWB(i,11)*2.0D0+DifWB(i,10)*(DifWB(i,8)*2.0D0-
     &              DifWB(i,7)*(DELTA(i+1)-DELTA(i))/XX))/DifWB(i,2)
	      ENDDO
            A2(1)=A1(1)
		  A2(ik)=A1(ik-1)
		  A2(2:ik-1)=(A1(2:ik-1)+A1(1:ik-2))*0.5D0
	      DO i=1,k
              REIM(k,j,i,is)=DELTA(i)
			REIM(l,j,i,is)=DELTA(ik+1-i)
              REIM(k,j,i,is+9)=A2(i)
			REIM(l,j,i,is+9)=A2(ik+1-i)
	      ENDDO
	    ENDDO
	    DO i=1,k
            REIM(k,j,i,7)=sum(REIM(k,j,i,1:6))
	      REIM(l,j,i,7)=sum(REIM(l,j,i,1:6))
	    ENDDO
	  ENDDO
	ENDDO

c high-latitude
	A(lmi)=0.0D0
	B(lmi)=1.0D0
	C(lmi)=0.0D0
	DO k=1,NHN
	  cHT(1:4,1:6)=0.0D0
	  DO i=1,NLi
	    Ni(1:3)=REIH(k,1,i,1:3)
          CALL Thdff(Ni,REIH(k,1,i,9),xt1,Cij,Beta) !Thdff(Ni,Ti,Cij,Betap)
c              CollisionFrequency(Ti,Tn,cO,cO2,cN2,cN,cH,cHe,Cin)
	    CALL CollisionFrequency(REIH(k,1,i,9),PNH(k,1,i,1),PNH(k,1,i,2
     &                            ),PNH(k,1,i,3),PNH(k,1,i,4),PNH(k,1,i,
     &                            7),PNH(k,1,i,6),PNH(k,1,i,5),Cin)
	    BsBo=BsiH(k,1,i,4)*bolt
		xtte=REIH(k,1,i,8)/REIH(k,1,i,7)
	    gtg=gxiH(k,1,i,3)+sum(EEH(k,1,1:2)*gxiH(k,1,i,1:2))
	    DvvT=EEH(k,1,1)*dvvH(k,1,i,1)+EEH(k,1,2)*dvvH(k,1,i,2)
	    DifWB(i,1)=BsiH(k,1,i,4)*WsiH(k,1,i)
	    IF(i.GE.NLi-1)THEN       !Top edge
	      tie=REIH(k,1,i,9)+REIH(k,1,i,8)
		  mtt=REIH(k,1,i,8)/tie
            xx=sim(3)*gtg
		  cHT(1,1:2)=cHT(1,1:2)+sim(1:2)*gtg-mtt*xx  !mg: cHT(1,:)
	      cHT(1,3)=cHT(1,3)+xx
		  mti=1.0D0-mtt
		  cHT(2,1:2)=cHT(2,1:2)+BsBo*mti !Bsk:  cHT(2,:)
		  cHT(2,3)=cHT(2,3)+BsBo
		  xy=Beta(3)*BsBo
		  cHT(3,3)=cHT(3,3)+xy
            cHT(3,1:2)=cHT(3,1:2)+Beta(1:2)*BsBo-mtt*xy !Bs*k*beta:  cHT(3,:)
            cHT(4,3)=cHT(4,3)+BsBo*tie    !Bs*k*(Ti+Te):  cHT(4,:)
         cHT(4,1:2)=cHT(4,1:2)+BsBo*(REIH(k,1,i,9)+REIH(k,1,i,1:2)*xtte)
          ENDIF
	    IF(i.LT.NLi)THEN
	      DO is=1,3,isd
	        Xi=sim(is)*(sum(Cij(is,1:3))+Cin(is)+0.004)                  !Xi
	        DifC(i,is,3)=BsBo/Xi
			DifC(i,is,8)=Xi/sim(is)
           DifC(i,is,1)=DifC(i,is,3)*(REIH(k,1,i,9)+xtte*REIH(k,1,i,is)) !Ri
	        DifC(i,is,2)=DifC(i,is,3)*xtte
		    DifC(i,is,4)=DifC(i,is,3)*Beta(is)
	        DifC(i,is,5)=sim(is)*(gtg+sum(REIH(k,1,i,10:12)*Cij(is,1:3
     &                     ))+PNH(k,1,i,8)*(Cin(is)))/Xi
	        DifC(i,is,6)=REIHold(k,1,i,is)/timestep+PLIH(k,1,i,is)
	        DifC(i,is,7)=PLIH(k,1,i,is+3)
              IF(DvvT.GE.0.0D0)DifC(i,is,7)=DifC(i,is,7)+DvvT
              IF(DvvT.LT.0.0D0)DifC(i,is,6)=DifC(i,is,6)-DvvT*REIHold(k,
     &                                      1,i,is)
	        DifC(i,is,9)=REIHold(k,1,i,is+5)
	      ENDDO
          ENDIF
	    ij=lmi-i+1
		Ni(1:3)=REIH(k,2,i,1:3)
          CALL Thdff(Ni,REIH(k,2,i,9),xt1,Cij,Beta)
	    CALL CollisionFrequency(REIH(k,2,i,9),PNH(k,2,i,1),PNH(k,2,i,2
     &                            ),PNH(k,2,i,3),PNH(k,2,i,4),PNH(k,2,i,
     &                            7),PNH(k,2,i,6),PNH(k,2,i,5),Cin)
	    BsBo=BsiH(k,2,i,4)*bolt
		xtte=REIH(k,2,i,8)/REIH(k,2,i,7)
	    gtg=gxiH(k,2,i,3)+sum(EEH(k,2,1:2)*gxiH(k,2,i,1:2))
	    DvvT=EEH(k,2,1)*dvvH(k,2,i,1)+EEH(k,2,2)*dvvH(k,2,i,2)
	    DifWB(ij,1)=BsiH(k,2,i,4)*WsiH(k,2,i)
	    IF(i.GE.NLi-1)THEN       !Top edge
	      tie=REIH(k,2,i,9)+REIH(k,2,i,8)
		  mtt=REIH(k,2,i,8)/tie
            xx=sim(3)*gtg
		  cHT(1,4:5)=cHT(1,4:5)+sim(1:2)*gtg-mtt*xx  !mg: cHT(1,:)
	      cHT(1,6)=cHT(1,6)+xx
		  mti=1.0D0-mtt
		  cHT(2,4:5)=cHT(2,4:5)+BsBo*mti
		  cHT(2,6)=cHT(2,6)+BsBo !Bsk:  cHT(2,:)
		  xy=Beta(3)*BsBo
		  cHT(3,6)=cHT(3,6)+xy
            cHT(3,4:5)=cHT(3,4:5)+Beta(1:2)*BsBo-mtt*xy !Bs*k*beta:  cHT(3,:)
            cHT(4,6)=cHT(4,6)+BsBo*tie    !Bs*k*(Ti+Te):  cHT(4,:)
         cHT(4,4:5)=cHT(4,4:5)+BsBo*(REIH(k,2,i,9)+REIH(k,2,i,1:2)*xtte)
          ENDIF
	    IF(i.LT.NLi)THEN
	      DO is=1,3,isd
	        Xi=sim(is)*(sum(Cij(is,1:3))+Cin(is)+0.004)
	        DifC(ij,is,3)=BsBo/Xi
			DifC(ij,is,8)=Xi/sim(is)
         DifC(ij,is,1)=DifC(ij,is,3)*(REIH(k,2,i,9)+xtte*REIH(k,2,i,is))
	        DifC(ij,is,2)=DifC(ij,is,3)*xtte
		    DifC(ij,is,4)=DifC(ij,is,3)*Beta(is)
	        DifC(ij,is,5)=sim(is)*(gtg+sum(REIH(k,2,i,10:12)*Cij(is,1:
     &                      3))+PNH(k,2,i,8)*(Cin(is)))/Xi
	        DifC(ij,is,6)=REIHold(k,2,i,is)/timestep+PLIH(k,2,i,is)
	        DifC(ij,is,7)=PLIH(k,2,i,is+3)
              IF(DvvT.GE.0.0D0)DifC(ij,is,7)=DifC(ij,is,7)+DvvT
              IF(DvvT.LT.0.0D0)DifC(ij,is,6)=DifC(ij,is,6)-DvvT*REIHold(
     &                                       k,2,i,is)
	        DifC(ij,is,9)=REIHold(k,2,i,is+5)
	      ENDDO
          ENDIF
        ENDDO
	  DO i=1,NLi-1
	    DifWB(i,3)=1.0D0/(VsiH(k,1,i+1)-VsiH(k,1,i))
          DifWB(i,2)=(DifWB(i,1)+DifWB(i+1,1))*0.5D0
          DifWB(i,4:5)=(REIH(k,1,i+1,8:9)-REIH(k,1,i,8:9))*DifWB(i,3)
	    DifWB(i,6)=2.0D0/(WsiH(k,1,i+1)+WsiH(k,1,i))
	    ij=lmi-i
		DifWB(ij,3)=1.0D0/(VsiH(k,2,i)-VsiH(k,2,i+1))
          DifWB(ij,2)=(DifWB(ij,1)+DifWB(ij+1,1))*0.5D0
          DifWB(ij,4:5)=(REIH(k,2,i,8:9)-REIH(k,2,i+1,8:9))*DifWB(ij,3)
	    DifWB(ij,6)=2.0D0/(WsiH(k,2,i+1)+WsiH(k,2,i))
	  ENDDO
	  DO is=1,3,isd
	    DO i=1,NLi-2     !NLi-1
            dtei=DifWB(i,4)+DifWB(i,5)
		  dti=DifWB(i,5)
            dei=((REIH(k,1,i+1,7)-REIH(k,1,i+1,is))-(REIH(k,1,i,7)-
     &          REIH(k,1,i,is)))
	      xx=dabs(dei)/dmax1(1.0D-30,dabs(REIH(k,1,i,is)-REIH(k,1,i+
     &         1,is)))
	      dei=dei*DifWB(i,3)*dexp(Sif(is)*xx)         !if(xx.ge.1.0D0)
            DifWB(i,7)=DifWB(i,2)*(DifC(i,is,1)+DifC(i+1,is,1))*0.5D0*
     &                 DifWB(i,3)
            DifWB(i,8)=DifWB(i,2)*((DifC(i,is,5)+DifC(i+1,is,5))-(DifC(i
     &              ,is,2)+DifC(i+1,is,2))*dei-(DifC(i,is,3)+DifC(i+1,is
     &              ,3))*dtei-(DifC(i,is,4)+DifC(i+1,is,4))*dti)*0.25D0
	      ij=lmi-i
		  dtei=DifWB(ij,4)+DifWB(ij,5)
		  dti=DifWB(ij,5)
            dei=((REIH(k,2,i,7)-REIH(k,2,i,is))-(REIH(k,2,i+1,7)-REIH(
     &          k,2,i+1,is)))
	      xx=dabs(dei)/dmax1(1.0D-30,dabs(REIH(k,2,i,is)-REIH(k,2,i+
     &         1,is)))
	      dei=dei*DifWB(ij,3)*dexp(Sif(is)*xx)  !if(xx.ge.1.0D0)
            DifWB(ij,7)=DifWB(ij,2)*(DifC(ij,is,1)+DifC(ij+1,is,1))*
     &                  0.5D0*DifWB(ij,3)
            DifWB(ij,8)=DifWB(ij,2)*((DifC(ij,is,5)+DifC(ij+1,is,5))-(
     &        DifC(ij,is,2)+DifC(ij+1,is,2))*dei-(DifC(ij,is,3)+DifC(ij+
     &        1,is,3))*dtei-(DifC(ij,is,4)+DifC(ij+1,is,4))*dti)*0.25D0
     		ENDDO
	    i=NLi-1
		ij=NLi+1
		iv=is+3    !top edge
          dtei=DifWB(i,4)+DifWB(i,5)
		dti=DifWB(i,5)
          DifWB(i,7)=cHT(4,is)*DifWB(i,3)
          DifWB(i,8)=(cHT(1,is)-cHT(2,is)*dtei-cHT(3,is)*dti)*0.5D0
          dtei=DifWB(ij,4)+DifWB(ij,5)
		dti=DifWB(ij,5)
          DifWB(ij,7)=cHT(4,iv)*DifWB(ij,3)
          DifWB(ij,8)=(cHT(1,iv)-cHT(2,iv)*dtei-cHT(3,iv)*dti)*0.5D0
          D(1)=dmax1(1.0D-25,(REIH(k,1,1,is)+PLIH(k,1,1,is)*TimeStep)/(
     &            PLIH(k,1,1,is+3)*TimeStep+1.0D0))
	    DifWB(1:NLi-1,9)=1.0D0/DifWB(1:NLi-1,3)
	    DifWB(NLi+1:lmi-1,9)=1.0D0/DifWB(NLi+1:lmi-1,3)
	    IF(kVC.eq.1)THEN
	      DifWB(1:lmi-1,10)=1.0D0
	      DifWB(1:lmi-1,11)=0.0D0
	    ELSE
	      DifWB(1:lmi-1,10)=Dexp(0.0D0-(DifC(1:lmi-1,is,8)+DifC(2:lmi,
     &               is,8))*0.5D0*TimeStep)
	      DifWB(1:lmi-1,11)=0.25D0*(DifC(1:lmi-1,is,9)+DifC(2:lmi,is,9
     &               ))*DifWB(1:lmi-1,2)*DifWB(1:lmi-1,10) !B*W*exp(-vt)*V_old/2
            DifWB(1:lmi-1,10)=1.0D0-DifWB(1:lmi-1,10)
	    ENDIF
          DO i=2,NLi-1
	      WDV=(DifWB(i-1,6)+DifWB(i,6))/(DifWB(i-1,9)+DifWB(i,9))
            X1=DifWB(i-1,7)
	      X3=DifWB(i-1,10)
	      X7=DifWB(i,10)
            X2=DifWB(i-1,8)
	      X4=DifWB(i-1,11)
	      X8=DifWB(i,11)
		  X5=DifWB(i,7)
	      X6=DifWB(i,8)
	      xx=DifC(i,is,7)    
		  A(i)=(-X4-X3*(X2+X1))*WDV
	      C(i)=(X8-X7*(X5-X6))*WDV
		  B(i)=1.0D0/timestep+xx+(X3*(X1-X2)-X4+X8+X7*(X6+X5))*WDV
            D(i)=DifC(i,is,6)
	    ENDDO
		A(NLi)=-DifWB(NLi-1,7)-DifWB(NLi-1,8)
		B(NLi)=DifWB(NLi-1,7)-DifWB(NLi-1,8)
	    C(NLi)=0.0D0
          D(NLi)=0.0D0
		CALL TRIDAGI(DELTA,lmi,1,NLi,A,B,C,D)
          DELTA(1:NLi)=dmax1(1.0D-25,DELTA(1:NLi))
          DO i=1,NLi-1
	      XX=(DELTA(i)+DELTA(i+1))*0.5D0
	      A1(i)=(DifWB(i,11)*2.0D0+DifWB(i,10)*(DifWB(i,8)*2.0D0-DifWB
     &            (i,7)*(DELTA(i+1)-DELTA(i))/XX))/DifWB(i,2)
	    ENDDO
          A2(1)=A1(1)
		A2(NLi)=A1(NLi-1)   !REIH(k,1,1,is+9)=A1(1)
	    A2(2:NLi-1)=(A1(2:NLi-1)+A1(1:NLi-2))*0.5D0
          DO i=NLi+2,lmi-1
	      WDV=(DifWB(i-1,6)+DifWB(i,6))/(DifWB(i-1,9)+DifWB(i,9))
            X1=DifWB(i-1,7)
            X2=DifWB(i-1,8)
            X3=DifWB(i-1,10)
            X4=DifWB(i-1,11)
		  X5=DifWB(i,7)
            X6=DifWB(i,8)
            X7=DifWB(i,10)
            X8=DifWB(i,11)
            xx=DifC(i,is,7)
		  A(i)=(-X4-X3*(X2+X1))*WDV
            C(i)=(X8-X7*(X5-X6))*WDV
		  B(i)=1.0D0/timestep+xx+(X3*(X1-X2)-X4+X8+X7*(X6+X5))*WDV
            D(i)=DifC(i,is,6)
	    ENDDO
	    A(NLi+1)=0.0D0
		C(NLi+1)=DifWB(NLi+1,8)-DifWB(NLi+1,7)
          D(NLi+1)=0.0D0
		B(NLi+1)=DifWB(NLi+1,7)+DifWB(NLi+1,8)
          D(lmi)=dmax1(1.0D-25,(REIH(k,2,1,is)+PLIH(k,2,1,is)*TimeStep)/
     &            (PLIH(k,2,1,is+3)*TimeStep+1.0D0))
          CALL TRIDAGI(DELTA,lmi,NLi+1,lmi,A,B,C,D)
          DELTA(NLi+1:lmi)=dmax1(1.0D-25,DELTA(NLi+1:lmi))
          DO i=NLi+1,lmi-1
	      XX=(DELTA(i)+DELTA(i+1))*0.5D0
	      A1(i)=(DifWB(i,11)*2.0D0+DifWB(i,10)*(DifWB(i,8)*2.0D0-DifWB
     &            (i,7)*(DELTA(i+1)-DELTA(i))/XX))/DifWB(i,2)
	    ENDDO
          A2(lmi)=A1(lmi-1)
		A2(NLi+1)=A1(NLi+1)
	    A2(NLi+2:lmi-1)=(A1(NLi+2:lmi-1)+A1(NLi+1:lmi-2))*0.5D0
	    DO i=1,NLi                      !REIH(k,2,1,is+9)=A2(lmi)
            REIH(k,1,i,is)=DELTA(i)
		  REIH(k,2,i,is)=DELTA(lmi+1-i)
            REIH(k,1,i,is+9)=A2(i)
		  REIH(k,2,i,is+9)=A2(lmi+1-i)
	    ENDDO
	  ENDDO
	  DO i=1,NLi
          REIH(k,1,i,7)=sum(REIH(k,1,i,1:6))
	    REIH(k,2,i,7)=sum(REIH(k,2,i,1:6))
	  ENDDO
	ENDDO
      END SUBROUTINE ImodelDensityEquation
c---------------------------------------------------------------------
      SUBROUTINE ImodelEnergyEquation(TimeStep)
      INCLUDE 'EIG.h'
	COMMON/IGrid/ d1i(lmi,NFi,NLi,3),d2i(lmi,NFi,NLi,3),glati(lmi,NFi,
     &              NLi),gloni(lmi,NFi,NLi),Vsi(lmi,NFi,NLi),Wsi(lmi,NFi
     &              ,NLi),Bsi(lmi,NFi,NLi,4),gxi(lmi,NFi,NLi,3),dvv(lmi,
     &              NFi,NLi,2),mlov(NLi,NFi,2),d1ih(NHN,2,NLi,3),d2ih(
     &              NHN,2,NLi,3),glatih(NHN,2,NLi),glonih(NHN,2,NLi),
     &              Vsih(NHN,2,NLi),Wsih(NHN,2,NLi),Bsih(NHN,2,NLi,4),
     &              gxih(NHN,2,NLi,3),dvvh(NHN,2,NLi,2),mlovh(NHN,2,2),
     &              MH300(NHN+NLi),i300
	COMMON/EEmC/  EEm(NLi,NFi,2),EEH(NHN,2,2)
	COMMON/REIM3/ REIM(lmi,NFi,NLi,12),REIMold(lmi,NFi,NLi,8),REIH(NHN
     &             ,2,NLi,12),REIHold(NHN,2,NLi,8)
	COMMON/REIM3N/ PN(lmi,NFi,NLi,8),PLI(lmi,NFi,NLi,8),PNH(NHN,2,NLi,
     &               8),PLIH(NHN,2,NLi,8)
      COMMON/consti/ Bolt,Sim(3)
      REAL*8 DELTA(lmi),A(lmi),B(lmi),C(lmi),D(lmi),Ni(6),DWB(lmi,7),Xci
     &       (lmi,14),ne(lmi),co2(lmi),cn2(lmi),co(lmi),prodelec(lmi),
     &       length(lmi-1),yta(lmi),heatout(lmi),tei(2),tn(lmi),DEL(lmi)

c mid- and low-latitude
      REIM(1,1:NFi,1,8)=PN(1,1:NFi,1,1)
      REIM(1,1:NFi,1,9)=PN(1,1:NFi,1,1)
	REIM(lmi,1:NFi,1,8:9)=REIM(1,1:NFi,1,8:9)
	A(1)=0.0D0
	B(1)=1.0D0
	C(1)=0.0D0  !xt=timestep*0.0D-4
	DO k=2,NLi
	  ik=2*k-1
	  l=lmi-k+1
	  A(ik)=0.0D0
	  B(ik)=1.0D0
	  C(ik)=0.0D0
	  DO j=1,NFi
	    DO i=1,k
		  co(i)=PN(k,j,i,2)
	      co2(i)=PN(k,j,i,3)
	      cn2(i)=PN(k,j,i,4)
		  tn(i)=PN(k,j,i,1)
		  Ni(1:6)=REIM(k,j,i,1:6)
		  ne(i)=REIM(k,j,i,7)
		  prodelec(i)=PLI(k,j,i,7)
		  tei(1:2)=REIM(k,j,i,8:9)
	      Xci(i,1)=EThermCon(co(i),cn2(i),co2(i),ne(i),tei(1))
            Xci(i,2)=xIonThermConductivity(Ni(1),Ni(2),Ni(3),tei(2))
            CALL ElectronHeatingRate(co(i),cn2(i),co2(i),Ni,PN(k,j,i,6),
     &                      ne(i),tn(i),tei(1),tei(2),Xci(i,3),Xci(i,4))
            CALL IonHeatingRates(tei(2),tn(i),tei(1),co(i),cn2(i),co2(i)
     &                          ,PN(k,j,i,7),PN(k,j,i,6),PN(k,j,i,5),PLI
     &                          (k,j,i,8),ne(i),Ni,Xci(i,5),Xci(i,6))
	      Xci(i,14)=bolt*ne(i)*sum(EEm(k,j,1:2)*dvv(k,j,i,1:2))
		  Xci(i,4)=-Xci(i,4)
		  Xci(i,6)=-Xci(i,6)
	      Xci(i,7)=sum(REIM(k,j,i,1:3)*REIM(k,j,i,10:12))/ne(i)*Bsi(k,
     &               j,i,4)
	      Xci(i,8)=Bsi(k,j,i,4)
            Xci(i,9)=Wsi(k,j,i)
		  Xci(i,10)=1.5D0*ne(i)*bolt
		  Xci(i,11:12)=REIMold(k,j,i,4:5)
            Xci(i,13)=Xci(i,9)*Xci(i,8)*Xci(i,8)
	    if(i.ne.k)then
            ij=ik-i+1
		  tn(ij)=PN(l,j,i,1)
		  co(ij)=PN(l,j,i,2)
	      co2(ij)=PN(l,j,i,3)
		  cn2(ij)=PN(l,j,i,4)
		  ne(ij)=REIM(l,j,i,7)
		  Ni(1:6)=REIM(l,j,i,1:6)
	      tei(1:2)=REIM(l,j,i,8:9)
	      prodelec(ij)=PLI(l,j,i,7)
	      Xci(ij,1)=EThermCon(co(ij),cn2(ij),co2(ij),ne(ij),tei(1))
            Xci(ij,2)=xIonThermConductivity(Ni(1),Ni(2),Ni(3),tei(2))
            CALL ElectronHeatingRate(co(ij),cn2(ij),co2(ij),Ni,PN(l,j,i,
     &               6),ne(ij),tn(ij),tei(1),tei(2),Xci(ij,3),Xci(ij,4))
            CALL IonHeatingRates(tei(2),tn(ij),tei(1),co(ij),cn2(ij),co2
     &                         (ij),PN(l,j,i,7),PN(l,j,i,6),PN(l,j,i,5),
     &                       PLI(l,j,i,8),ne(ij),Ni,Xci(ij,5),Xci(ij,6))
            Xci(ij,14)=bolt*ne(ij)*sum(EEm(k,j,1:2)*dvv(l,j,i,1:2))
            Xci(ij,4)=-Xci(ij,4)
            Xci(ij,6)=-Xci(ij,6)
	      Xci(ij,7)=sum(REIM(l,j,i,1:3)*REIM(l,j,i,10:12))/ne(ij)*Bsi(
     &                l,j,i,4)
	      Xci(ij,8)=Bsi(l,j,i,4)
            Xci(ij,9)=Wsi(l,j,i)
		  Xci(ij,10)=1.5D0*ne(ij)*bolt
		  Xci(ij,11:12)=REIMold(l,j,i,4:5)
            Xci(ij,13)=Xci(ij,9)*Xci(ij,8)*Xci(ij,8)
	    endif
          ENDDO
          DWB(1:ik-1,2)=2.0D0/(Xci(2:ik,9)+Xci(1:ik-1,9))
	    DWB(1:ik-1,4)=Dsqrt(Xci(2:ik,13)*Xci(2:ik,1)*Xci(1:ik-1,13)*
     &                  Xci(1:ik-1,1))
	    DWB(1:ik-1,5)=Dsqrt(Xci(2:ik,13)*Xci(2:ik,2)*Xci(1:ik-1,13)*
     &                  Xci(1:ik-1,2))
          DWB(1:ik-1,6)=(Xci(2:ik,10)*Xci(2:ik,7)+Xci(1:ik-1,10)*Xci(1:
     &                  ik-1,7))*0.5D0
	    DO i=1,k-1
	      DWB(i,1)=Vsi(k,j,i+1)-Vsi(k,j,i)
	      DWB(i,3)=(Xci(i+1,9)*Xci(i+1,7)-Xci(i,9)*Xci(i,7))/DWB(i,1)*
     &               DWB(i,2)
	      ij=ik-i
		  DWB(ij,1)=Vsi(l,j,i)-Vsi(l,j,i+1)
	      DWB(ij,3)=(Xci(ij+1,9)*Xci(ij+1,7)-Xci(ij,9)*Xci(ij,7))/DWB(
     &                ij,1)*DWB(ij,2)
	    ENDDO
		yta(1:ik)=Xci(1:ik,8)
		heatout(1:lmi)=0.0D0
	    length(1:ik-1)=Dabs(DWB(1:ik-1,1)/dsqrt(Xci(1:ik-1,8)*Xci(2:ik
     &                   ,8)))
          call photoelectronheatUp(lmi,ik,i300+1,ik-i300,Ne,co2,cn2,co,
     &                               prodelec,length,yta,heatout)
	    Xci(1:ik,3)=Xci(1:ik,3)+heatout(1:ik)
	    Xci(1:ik,3:6)=Xci(1:ik,3:6)*1.6D-19
          D(1)=PN(k,j,1,1)
		D(ik)=PN(l,j,1,1)
		xxc=0.0D0
          DO i=2,ik-1
	      Wm=0.25D0*(Xci(i-1,9)+Xci(i+1,9))+0.5D0*Xci(i,9)
            XX1=2.0D0*DWB(i-1,4)/(DWB(i-1,1)*(DWB(i-1,1)+DWB(i,1))*Wm)
            XX2=2.0D0*DWB(i,4)/(DWB(i,1)*(DWB(i-1,1)+DWB(i,1))*Wm)
            B(i)=Xci(i,10)/timestep
		  D(i)=B(i)*Xci(i,11)+Xci(i,3)
            XX3=DWB(i-1,6)/DWB(i-1,1)*0.5D0
		  XX4=DWB(i,6)/DWB(i,1)*0.5D0
            xx=Xci(i,14)+(DWB(i-1,3)+DWB(i,3))*Xci(i,10)/3.0D0
	      IF(xx.LT.0.0D0)THEN
		    D(i)=D(i)-xx*Xci(i,11)
              xx=0.0D0
	      ENDIF
	      A(i)=-XX3-XX1
		  C(i)=XX4-XX2
		  B(i)=B(i)+(XX1+XX2)-XX4+XX3+xx+Xci(i,4)
	      if(A(i).GE.0.0D0.or.C(i).GE.0.0D0)xxc=1.0D0
          ENDDO
          CALL TRIDAGI(DELTA,lmi,1,ik,A,B,C,D)
          DELTA(1:ik)=dmax1(tn(1:ik),DELTA(1:ik))
          DELTA(1:ik)=dmin1(8.0D3,DELTA(1:ik))
	    DEL(1:ik)=DELTA(1:ik)
          DO i=2,ik-1
	      Wm=0.25D0*(Xci(i-1,9)+Xci(i+1,9))+0.5D0*Xci(i,9)
            XX1=2.0D0*DWB(i-1,5)/(DWB(i-1,1)*(DWB(i-1,1)+DWB(i,1))*Wm)
            XX2=2.0D0*DWB(i,5)/(DWB(i,1)*(DWB(i-1,1)+DWB(i,1))*Wm)
            XX3=DWB(i-1,6)/DWB(i-1,1)*0.5D0
		  XX4=DWB(i,6)/DWB(i,1)*0.5D0
            B(i)=Xci(i,10)/timestep
		  D(i)=B(i)*Xci(i,12)+Xci(i,5)
            xx=Xci(i,14)+(DWB(i-1,3)+DWB(i,3))*Xci(i,10)/3.0D0
	      IF(xx.LT.0.0D0)THEN
		    D(i)=D(i)-xx*Xci(i,12)
              xx=0.0D0
	      ENDIF
	      A(i)=-XX3-XX1
		  C(i)=XX4-XX2
		  B(i)=B(i)+(XX1+XX2)-XX4+XX3+xx+Xci(i,6)
	      if(dabs(A(i))+dabs(C(i)).GE.dabs(B(i)))xxc=1.0D0
          ENDDO
          CALL TRIDAGI(DELTA,lmi,1,ik,A,B,C,D)
          DELTA(1:ik)=dmax1(tn(1:ik),DELTA(1:ik))
          DELTA(1:ik)=dmin1(DEL(1:ik),DELTA(1:ik))
	    DO i=1,k
            REIM(k,j,i,8)=DEL(i)
		  REIM(l,j,i,8)=DEL(ik-i+1)
            REIM(k,j,i,9)=DELTA(i)
		  REIM(l,j,i,9)=DELTA(ik-i+1)
	    ENDDO
	  ENDDO
	ENDDO

c high-latitude
	A(lmi)=0.0D0
	B(lmi)=1.0D0
	C(lmi)=0.0D0
	DO k=1,NHN
        DO i=1,NLi
		co(i)=PNH(k,1,i,2)
          co2(i)=PNH(k,1,i,3)
	    cn2(i)=PNH(k,1,i,4)
		tn(i)=PNH(k,1,i,1)
		ne(i)=REIH(k,1,i,7)
		Ni(1:6)=REIH(k,1,i,1:6)
		tei(1:2)=REIH(k,1,i,8:9)
		prodelec(i)=PLIH(k,1,i,7)
	    Xci(i,1)=EThermCon(co(i),cn2(i),co2(i),ne(i),tei(1))
          Xci(i,2)=xIonThermConductivity(Ni(1),Ni(2),Ni(3),tei(2))
          CALL ElectronHeatingRate(co(i),cn2(i),co2(i),Ni,PNH(k,1,i,6),
     &                     ne(i),tn(i),tei(1),tei(2),Xci(i,3),Xci(i,4))
          CALL IonHeatingRates(tei(2),tn(i),tei(1),co(i),cn2(i),co2(i),
     &                       PNH(k,1,i,7),PNH(k,1,i,6),PNH(k,1,i,5),PLIH
     &                       (k,1,i,8),ne(i),Ni,Xci(i,5),Xci(i,6))
	    Xci(i,14)=bolt*ne(i)*sum(EEH(k,1,1:2)*dvvh(k,1,i,1:2))
	  	Xci(i,4)=-Xci(i,4)
		Xci(i,6)=-Xci(i,6)
	    Xci(i,7)=sum(REIH(k,1,i,1:3)*REIH(k,1,i,10:12))/ne(i)*BsiH(k,1
     &             ,i,4)
	    Xci(i,8)=BsiH(k,1,i,4)
          Xci(i,9)=WsiH(k,1,i)
		Xci(i,10)=1.5D0*ne(i)*bolt
		Xci(i,11:12)=REIHold(k,1,i,4:5)
          Xci(i,13)=Xci(i,9)*Xci(i,8)*Xci(i,8)
          ij=lmi-i+1
		co(ij)=PNH(k,2,i,2)
		co2(ij)=PNH(k,2,i,3)
		cn2(ij)=PNH(k,2,i,4)
		tn(ij)=PNH(k,2,i,1)
		Ni(1:6)=REIH(k,2,i,1:6)
		ne(ij)=REIH(k,2,i,7)
	    tei(1:2)=REIH(k,2,i,8:9)
	    prodelec(ij)=PLIH(k,2,i,7)
	    Xci(ij,1)=EThermCon(co(ij),cn2(ij),co2(ij),ne(ij),tei(1))
          Xci(ij,2)=xIonThermConductivity(Ni(1),Ni(2),Ni(3),tei(2))
          CALL ElectronHeatingRate(co(ij),cn2(ij),co2(ij),Ni,PNH(k,2,i,6
     &                ),ne(ij),tn(ij),tei(1),tei(2),Xci(ij,3),Xci(ij,4))
          CALL IonHeatingRates(tei(2),tn(ij),tei(1),co(ij),cn2(ij),co2(
     &                       ij),PNH(k,2,i,7),PNH(k,2,i,6),PNH(k,2,i,5),
     &                      PLIH(k,2,i,8),ne(ij),Ni,Xci(ij,5),Xci(ij,6))
          Xci(ij,14)=bolt*ne(ij)*sum(EEH(k,2,1:2)*dvvh(k,2,i,1:2))
          Xci(ij,4)=-Xci(ij,4)
	    Xci(ij,6)=-Xci(ij,6)
	    Xci(ij,7)=sum(REIH(k,2,i,1:3)*REIH(k,2,i,10:12))/ne(ij)*BsiH(k
     &              ,2,i,4)
	    Xci(ij,8)=BsiH(k,2,i,4)
          Xci(ij,9)=WsiH(k,2,i)
		Xci(ij,10)=1.5D0*ne(ij)*bolt
		Xci(ij,11:12)=REIHold(k,2,i,4:5)
          Xci(ij,13)=Xci(ij,9)*Xci(ij,8)*Xci(ij,8)
        ENDDO
        DWB(1:lmi-1,2)=2.0D0/(Xci(2:lmi,9)+Xci(1:lmi-1,9))
	  DWB(1:lmi-1,4)=Dsqrt(Xci(2:lmi,13)*Xci(2:lmi,1)*Xci(1:lmi-1,13)*
     &                 Xci(1:lmi-1,1))
	  DWB(1:lmi-1,5)=Dsqrt(Xci(2:lmi,13)*Xci(2:lmi,2)*Xci(1:lmi-1,13)*
     &                 Xci(1:lmi-1,2))
        DWB(1:lmi-1,6)=(Xci(2:lmi,10)*Xci(2:lmi,7)+Xci(1:lmi-1,10)*Xci(1
     &                 :lmi-1,7))*0.5D0
	  DO i=1,NLi-1
	    DWB(i,1)=VsiH(k,1,i+1)-VsiH(k,1,i)
	    DWB(i,3)=(Xci(i+1,9)*Xci(i+1,7)-Xci(i,9)*Xci(i,7))/DWB(i,1)*
     &             DWB(i,2)
	    ij=lmi-i
		DWB(ij,1)=VsiH(k,2,i)-VsiH(k,2,i+1)
	    DWB(ij,3)=(Xci(ij+1,9)*Xci(ij+1,7)-Xci(ij,9)*Xci(ij,7))/DWB(
     &              ij,1)*DWB(ij,2)
	  ENDDO
	  DWB(NLi,1)=VsiH(k,2,NLi)-VsiH(k,1,NLi)
	  DWB(NLi,3)=(Xci(NLi+1,9)*Xci(NLi+1,7)-Xci(NLi,9)*Xci(NLi,7))/DWB
     &             (NLi,1)*DWB(NLi,2)
	  yta(1:lmi)=Xci(1:lmi,8)
	  heatout(1:lmi)=0.0D0
	  length(1:lmi-1)=Dabs(DWB(1:lmi-1,1)/dsqrt(Xci(1:lmi-1,8)*Xci(2:
     &                  lmi,8)))
        call photoelectronheatUp(lmi,lmi,i300+1,lmi-i300,Ne,co2,cn2,co,
     &                           prodelec,length,yta,heatout)
	  Xci(1:lmi,3)=Xci(1:lmi,3)+heatout(1:lmi)
	  Xci(1:lmi,3:6)=Xci(1:lmi,3:6)*1.6D-19
        D(1)=PNH(k,1,1,1)
	  D(lmi)=PNH(k,2,1,1)
	  xxc=0.0D0
        DO i=2,lmi-1
	    Wm=0.25D0*(Xci(i-1,9)+Xci(i+1,9))+0.5D0*Xci(i,9)
          XX1=2.0D0*DWB(i-1,4)/(DWB(i-1,1)*(DWB(i-1,1)+DWB(i,1))*Wm)
          XX2=2.0D0*DWB(i,4)/(DWB(i,1)*(DWB(i-1,1)+DWB(i,1))*Wm)
          B(i)=Xci(i,10)/timestep
		D(i)=B(i)*Xci(i,11)+Xci(i,3)
          XX3=DWB(i-1,6)/DWB(i-1,1)*0.5D0
		XX4=DWB(i,6)/DWB(i,1)*0.5D0
          xx=Xci(i,14)+(DWB(i-1,3)+DWB(i,3))*Xci(i,10)/3.0D0
	    IF(xx.LT.0.0D0)THEN
		  D(i)=D(i)-xx*Xci(i,11)
            xx=0.0D0
	    ENDIF
	    A(i)=-XX3-XX1
		C(i)=XX4-XX2
		B(i)=B(i)+(XX1+XX2)-XX4+XX3+xx+Xci(i,4)
	    if(A(i).GE.0.0D0.or.C(i).GE.0.0D0)xxc=1.0D0
        ENDDO
        CALL TRIDAGI(DELTA,lmi,1,lmi,A,B,C,D)
        DELTA(1:lmi)=dmax1(tn(1:lmi),DELTA(1:lmi))
        DELTA(1:lmi)=dmin1(8.0D3,DELTA(1:lmi))
	  DEL(1:lmi)=DELTA(1:lmi)
        DO i=2,lmi-1
	    Wm=0.25D0*(Xci(i-1,9)+Xci(i+1,9))+0.5D0*Xci(i,9)
          XX1=2.0D0*DWB(i-1,5)/(DWB(i-1,1)*(DWB(i-1,1)+DWB(i,1))*Wm)
          XX2=2.0D0*DWB(i,5)/(DWB(i,1)*(DWB(i-1,1)+DWB(i,1))*Wm)
          XX3=DWB(i-1,6)/DWB(i-1,1)*0.5D0
		XX4=DWB(i,6)/DWB(i,1)*0.5D0
          B(i)=Xci(i,10)/timestep
		D(i)=B(i)*Xci(i,12)+Xci(i,5)
          xx=Xci(i,14)+(DWB(i-1,3)+DWB(i,3))*Xci(i,10)/3.0D0
	    IF(xx.LT.0.0D0)THEN
		  D(i)=D(i)-xx*Xci(i,12)
            xx=0.0D0
	    ENDIF
	    A(i)=-XX3-XX1
		C(i)=XX4-XX2
		B(i)=B(i)+(XX1+XX2)-XX4+XX3+xx+Xci(i,6)
	    if(dabs(A(i))+dabs(C(i)).GE.dabs(B(i)))xxc=1.0D0
        ENDDO
        CALL TRIDAGI(DELTA,lmi,1,lmi,A,B,C,D)
        DELTA(1:lmi)=dmax1(tn(1:lmi),DELTA(1:lmi))
        DELTA(1:lmi)=dmin1(DEL(1:lmi),DELTA(1:lmi))
	  DO i=1,NLi
          REIH(k,1,i,8)=DEL(i)
		REIH(k,2,i,8)=DEL(lmi-i+1)
          REIH(k,1,i,9)=DELTA(i)
		REIH(k,2,i,9)=DELTA(lmi-i+1)
	  ENDDO
	ENDDO
      END SUBROUTINE ImodelEnergyEquation
c------------------------------------------------------------------
      SUBROUTINE GridInterpolation(TimeStep)
      INCLUDE 'EIG.h'
	COMMON/IGrid/ d1i(lmi,NFi,NLi,3),d2i(lmi,NFi,NLi,3),glati(lmi,NFi,
     &              NLi),gloni(lmi,NFi,NLi),Vsi(lmi,NFi,NLi),Wsi(lmi,NFi
     &              ,NLi),Bsi(lmi,NFi,NLi,4),gxi(lmi,NFi,NLi,3),dvv(lmi,
     &              NFi,NLi,2),mlov(NLi,NFi,2),d1ih(NHN,2,NLi,3),d2ih(
     &              NHN,2,NLi,3),glatih(NHN,2,NLi),glonih(NHN,2,NLi),
     &              Vsih(NHN,2,NLi),Wsih(NHN,2,NLi),Bsih(NHN,2,NLi,4),
     &              gxih(NHN,2,NLi,3),dvvh(NHN,2,NLi,2),mlovh(NHN,2,2),
     &              MH300(NHN+NLi)
	COMMON/GRIDI/ alti(NLi),Cli(NLi),SIi(NLi),MLATi(lmi+1,NLi),Mloni(
     &              NFi),altH(NHN),ClH(NHN),SIH(NHN),MlatH(NHN,2,NLi),
     &              MlonH(NHN)
	COMMON/EEmC/  EEm(NLi,NFi,2),EEH(NHN,2,2)
	COMMON/REIM3/ REIM(lmi,NFi,NLi,12),REIMold(lmi,NFi,NLi,8),REIH(NHN
     &             ,2,NLi,12),REIHold(NHN,2,NLi,8)
	COMMON/constant/ PI,DtR,dF
      COMMON/consti/ Bolt,Sim(3)
	COMMON/Iinterp/ Sl(NLi+nk+1,2),dmlon(NLi+nk),Xli(NLi+nk+1),MLATT(
     &                NLi+nk+1,NLi),NHni(NHN,2),Nmlo(Nli+nk),NHX(Nli+1:
     &                Nli+nk,NFi)
      integer di,ic(2),jc(2),ih(2)
      REAL*8 Vd(NLi+nk,NFi,2,2),tre(2),MLATC(2),MLONC(2),V(2,2),tic(2),
     &       REX(lmi+1,9,2),Vx(2,2),REY(lmi+1,9,4),X(lmi+1),Y(lmi+1,8),
     &       T(lmi+1),Z(lmi+1,8),altv(NLi),Vhp(NLi+1:NLi+nk+1,NFi,2),Vh
     &       (NLi+1:NLi+nk+1,NFi,2,2),VHpx(2,2),REZ(2,NLi,9)

        Vd(1:NLi,1:NFi,1,1:2)=EEm(1:NLi,1:NFi,1:2)*mlov(1:NLi,1:NFi,1:2)
        Vd(1:NLi,1:NFi,2,1:2)=Vd(1:NLi,1:NFi,1,1:2)
	  n=NHX(NLi+nk-6,1)-1
	  XliH=Xli(NHni(n,1))
	  PI2=2.0D0*PI
	  P2I=pi*0.5D0
	DO k=1,NHN
        Vd(NHni(k,1),NHni(k,2),1:2,1:2)=EEH(k,1:2,1:2)*mlovH(k,1:2,1:2)
	  IF(k.LE.n)THEN       ! higher than 76 degree latitude
	    mlonx=dmlon(NHni(k,1))*DBLE(NHni(k,2)-1)
		mlatx=P2I-Xli(NHni(k,1))
		cmlon=dcos(mlonx)
		smlon=dsin(mlonx)
	    Vx(1:2,1:2)=Vd(NHni(k,1),NHni(k,2),1:2,1:2)
		Vx(1:2,2)=Vx(1:2,2)*dsin(mlatx)
	    Vh(NHni(k,1),NHni(k,2),1:2,1)=-Vx(1:2,1)*cmlon-Vx(1:2,2)*smlon
	    Vh(NHni(k,1),NHni(k,2),1:2,2)=Vx(1:2,2)*cmlon-Vx(1:2,1)*smlon
          Vhp(NHni(k,1),NHni(k,2),1)=mlatx*cmlon   
		Vhp(NHni(k,1),NHni(k,2),2)=mlatx*smlon
	  ENDIF
	ENDDO
c	 Vh(NLi+nk+1,1,1:2,1:2)=(Vh(NLi+nk,1,1:2,1:2)+Vh(NLi+nk,2,1:2,1:2)+Vh(NLi+nk,3,1:2,1:2))/3.0D0
c	 REZ(1:2,1:NLi,1:3)=Dlog((REIH(1,1:2,1:NLi,1:3)+REIH(2,1:2,1:NLi,1:3)+REIH(3,1:2,1:NLi,1:3))/3.0D0)
c       REZ(1:2,1:NLi,4:8)=(REIH(1,1:2,1:NLi,8:12)+REIH(2,1:2,1:NLi,8:12)+REIH(3,1:2,1:NLi,8:12))/3.0D0
c       REZ(1:2,1:NLi,9)=(VsiH(1,1:2,1:NLi)+VsiH(2,1:2,1:NLi)+VsiH(3,1:2,1:NLi))/3.0D0
	 Vh(NLi+nk+1,1,1:2,1:2)=0.0D0
	 REZ(1:2,1:NLi,1:9)=0.0D0
      DO i=1,Nmlo(NLi+nk)
	Vh(NLi+nk+1,1,1:2,1:2)=Vh(NLi+nk+1,1,1:2,1:2)+Vh(NLi+nk,i,1:2,1:2)
	  REZ(1:2,1:NLi,1:3)=REZ(1:2,1:NLi,1:3)+REIH(i,1:2,1:NLi,1:3)
        REZ(1:2,1:NLi,4:8)=REZ(1:2,1:NLi,4:8)+REIH(i,1:2,1:NLi,8:12)
        REZ(1:2,1:NLi,9)=REZ(1:2,1:NLi,9)+VsiH(i,1:2,1:NLi)
      ENDDO
	 Vh(NLi+nk+1,1,1:2,1:2)=Vh(NLi+nk+1,1,1:2,1:2)/dble(Nmlo(NLi+nk))
	 REZ(1:2,1:NLi,1:3)=Dlog(REZ(1:2,1:NLi,1:3)/dble(Nmlo(NLi+nk)))
       REZ(1:2,1:NLi,4:9)=REZ(1:2,1:NLi,4:9)/dble(Nmlo(NLi+nk))

	DO ix=1,NLi
        ji=NFi
	  IF(ix.eq.1)ji=NHN
	  DO jx=1,ji
		i=ix
		j=jx
		io=0
		jc(1:2)=0
          tre(1:2)=timestep
          if(ix.eq.1)THEN                                  !i=1ŒÆËãžßÎ³
		  i=NHni(jx,1)
		  j=NHni(jx,2)
		  io=1
          ENDIF
		V(1:2,1:2)=Vd(i,j,1:2,1:2)
		MLATC(1:2)=Xli(i)
		MLONC(1:2)=dmlon(i)*DBLE(j-1)
		ic(1:2)=i
	    IF(MLATC(1).GE.XliH)THEN
            jc(1:2)=1
		  V(1:2,1:2)=VH(i,j,1:2,1:2)
            VHpx(1,1:2)=Vhp(i,j,1:2)
		  VHpx(2,1:2)=Vhp(i,j,1:2)
	    ENDIF
      DO k=1,2
	  IF((k.eq.1).or.(io.eq.1))THEN
        DO while(tre(k).gt.0.0D0)           !Ê±ŒäÊ£ÓàÔòŒÌÐøŒÆËã
	    IF(MLATC(k).LT.XliH)THEN
	      DO iv=1,2
	        IF(DABS(V(k,iv))*tre(k).le.Sl(ic(k),iv))THEN
                tic(iv)=tre(k)*2.0D0
	        ELSE
                tic(iv)=Sl(ic(k),iv)/DABS(V(k,iv))
	        ENDIF
	      ENDDO
            dt=DMIN1(tre(k),tic(1),tic(2))
		  tre(k)=tre(k)-dt   !°ŽÊ£ÓàÊ±Œä¡¢ÒÔŒ°Áœžö·œÏò±ä»¯µÄÌØÕ÷Ê±ŒäÖÐ×îÐ¡µÄÒ»žöœøÐÐÊ±Œä²œœø
	      MLONC(k)=DMOD(MLONC(k)-V(k,2)*dt,PI2)
	      IF(MLONC(k).LT.0.0D0)MLONC(k)=MLONC(k)+PI2
		  MLATD=MLATC(k)
		  jc(k)=0
     		  MLATC(k)=MLATC(k)-V(k,1)*dt
		  IF(MLATC(k).GT.Xli(Nli))io=1
	    ELSE
		  DO iv=1,2
	        IF(DABS(V(k,iv))*tre(k).le.Sl(ic(k),1))THEN
                tic(iv)=tre(k)*2.0D0
	        ELSE
                tic(iv)=Sl(ic(k),1)/DABS(V(k,iv))
	        ENDIF
	      ENDDO
            dt=DMIN1(tre(k),tic(1),tic(2))
		  tre(k)=tre(k)-dt
	      IF(jc(k).eq.0)THEN
              mlatx=P2I-MLATC(k)
			VHpx(k,1)=mlatx*dcos(MLONC(k))
              VHpx(k,2)=mlatx*dsin(MLONC(k))
	      ENDIF
            VHpx(k,1:2)=VHpx(k,1:2)-V(k,1:2)*dt
		  MLATD=MLATC(k)
	      xm=Dsqrt(sum(VHpx(k,1:2)*VHpx(k,1:2)))
		  MLATC(k)=P2I-xm
		  io=1
		  jc(k)=1
            IF(xm.eq.0.0D0)THEN
              MLONC(k)=0.0D0
	      ELSE
              MLONC(k)=DACOS(VHpx(k,1)/xm)
	        IF(VHpx(k,2).LT.0.0D0)MLONC(k)=PI2-MLONC(k)
            ENDIF
	    ENDIF
          IF(MLATC(k).GE.XliH)THEN                             !Œ«ÇøŒÆËã
            IF(MLATC(k).GT.MLATD)THEN
              di=1
			IF(Xli(ic(k)).GT.MLATD)ic(k)=ic(k)-1
	      ELSE
              di=-1
			IF(Xli(ic(k)).LT.MLATD)ic(k)=ic(k)+1
            ENDIF
            i1=ic(k)+di
		  DO while((Xli(ic(k))-MLATC(k))*(Xli(i1)-MLATC(k)).GT.0.D0)
              ic(k)=i1
			i1=ic(k)+di
	      ENDDO
	      IF(ic(k).ne.NLi+nk+1)THEN
	        aj=MLONC(k)/dmlon(ic(k))+1.0D0
			j0=INT(aj)
			xx=aj-j0
			IF(j0.gt.Nmlo(ic(k)))j0=j0-Nmlo(ic(k))
		    j1=j0+1
			IF(j1.gt.Nmlo(ic(k)))j1=j1-Nmlo(ic(k))
              xx1=1.0D0-xx
	      ENDIF
            IF(tre(k).eq.0.0D0)THEN
	        IF(ic(k).eq.NLi+nk+1)THEN
 	          REY(1:NLi,1:9,1)=REZ(k,1:NLi,1:9)
	        ELSE
	          ix0=NHX(ic(k),j0)
			  ix1=NHX(ic(k),j1)
                REY(1:NLi,1:3,1)=Dlog(REIH(ix0,k,1:NLi,1:3)*xx1+REIH(ix1
     &                           ,k,1:NLi,1:3)*xx)   
                REY(1:NLi,4:8,1)=REIH(ix0,k,1:NLi,8:12)*xx1+REIH(ix1,k,1
     &                           :NLi,8:12)*xx
               REY(1:NLi,9,1)=VsiH(ix0,k,1:NLi)*xx1+VsiH(ix1,k,1:NLi)*xx
	        ENDIF
            ELSE
	        IF(ic(k).eq.NLi+nk+1)THEN
                Vx(1,1:2)=Vh(NLi+nk+1,1,k,1:2)
	        ELSE
                Vx(1,1:2)=Vh(ic(k),j0,k,1:2)*xx1+Vh(ic(k),j1,k,1:2)*xx
              ENDIF
	      ENDIF
	      IF(i1.ne.NLi+nk+1)THEN
	        aj=MLONC(k)/dmlon(i1)+1.0D0
			j0=INT(aj)
			xx=aj-j0
			IF(j0.gt.Nmlo(i1))j0=j0-Nmlo(i1)
			j1=j0+1
			IF(j1.gt.Nmlo(i1))j1=j1-Nmlo(i1)
			xx1=1.0D0-xx
	      ENDIF
            IF(tre(k).eq.0.0D0)THEN
	        IF(i1.eq.NLi+nk+1)THEN
 	          REY(1:NLi,1:9,2)=REZ(k,1:NLi,1:9)
	        ELSE
	          ix0=NHX(i1,j0)
			  ix1=NHX(i1,j1)
                REY(1:NLi,1:3,2)=Dlog(REIH(ix0,k,1:NLi,1:3)*xx1+REIH(ix1
     &                           ,k,1:NLi,1:3)*xx)  
                REY(1:NLi,4:8,2)=REIH(ix0,k,1:NLi,8:12)*xx1+REIH(ix1,k,1
     &                           :NLi,8:12)*xx
               REY(1:NLi,9,2)=VsiH(ix0,k,1:NLi)*xx1+VsiH(ix1,k,1:NLi)*xx
	        ENDIF
            ELSE
	        IF(i1.eq.NLi+nk+1)THEN
                Vx(2,1:2)=Vh(NLi+nk+1,1,k,1:2)
	        ELSE
                Vx(2,1:2)=Vh(i1,j0,k,1:2)*xx1+Vh(i1,j1,k,1:2)*xx
              ENDIF
	      ENDIF
            xx=(Xli(ic(k))-MLATC(k))/(Xli(ic(k))-Xli(i1))
		  xx1=1.0D0-xx
            IF(tre(k).eq.0.0D0)THEN
	        ih(k)=NLi     !i2=imin0(ic(k),i1,NLi)
              REX(1:NLi,1:9,k)=REY(1:NLi,1:9,1)*xx1+REY(1:NLi,1:9,2)*xx
            ELSE
              V(k,1:2)=Vx(1,1:2)*xx1+Vx(2,1:2)*xx
	      ENDIF
	    ELSE IF(MLATC(k).LE.Xli(1))THEN                          !Ž©³öÏÂ±ßœç
		  MLATC(k)=Xli(1)
		  ic(k)=1
		  MLONC(k)=DMOD(MLONC(k)-V(k,2)*tre(k),PI2)
		  tre(k)=0.0D0
	      IF(MLONC(k).LT.0.0D0)MLONC(k)=MLONC(k)+PI2
	      aj=MLONC(k)/dmlon(1)+1.0D0
		  j0=int(aj)
		  xx=aj-j0
		  if(j0.gt.Nmlo(1))j0=j0-Nmlo(1)
		  j1=j0+1
		  if(j1.gt.Nmlo(1))j1=j1-Nmlo(1)
		  xx1=1.0D0-xx
            REX(1,1:3,k)=Dlog(REIM(1,j0,1,1:3)*xx1+REIM(1,j1,1,1:3)*xx)
            REX(1,4:8,k)=REIM(1,j0,1,8:12)*xx1+REIM(1,j1,1,8:12)*xx
            REX(1,9,k)=Vsi(1,j0,1)*xx1+Vsi(1,j1,1)*xx
		  ih(1)=1
		  ih(2)=0
	    ELSE
c		  di=-1;     IF(V(k,1).LT.0.0D0)di=1;     i1=ic(k)+di;
            IF(MLATC(k).GT.MLATD)THEN
              di=1
			IF(Xli(ic(k)).GT.MLATD)ic(k)=ic(k)-1
	      ELSE
              di=-1
			IF(Xli(ic(k)).LT.MLATD)ic(k)=ic(k)+1
            ENDIF
            i1=ic(k)+di
		  DO while((Xli(ic(k))-MLATC(k))*(Xli(i1)-MLATC(k)).GT.0.D0)
              ic(k)=i1
			i1=ic(k)+di
	      ENDDO
	      aj=MLONC(k)/dmlon(ic(k))+1.0D0
		  j0=int(aj)
		  xx=aj-j0
		  IF(j0.gt.Nmlo(ic(k)))j0=j0-Nmlo(ic(k))
		  j1=j0+1
		  IF(j1.gt.Nmlo(ic(k)))j1=j1-Nmlo(ic(k))
            xx1=1.0D0-xx
            IF(tre(k).eq.0.0D0)THEN
              IF(ic(k).LE.NLi)THEN
	          ik=ic(k)
			  IF(k.eq.2)ik=lmi+1-ik
                REY(1:ic(k),1:3,1)=Dlog(REIM(ik,j0,1:ic(k),1:3)*xx1+REIM
     &                             (ik,j1,1:ic(k),1:3)*xx)
                REY(1:ic(k),4:8,1)=REIM(ik,j0,1:ic(k),8:12)*xx1+REIM(ik,
     &                             j1,1:ic(k),8:12)*xx
                REY(1:ic(k),9,1)=Vsi(ik,j0,1:ic(k))*xx1+Vsi(ik,j1,1:ic(k
     &                           ))*xx
                IF(io.eq.0)THEN
                  REY(1:ic(k),1:3,3)=Dlog(REIM(lmi+1-ic(k),j0,1:ic(k),1:
     &                       3)*xx1+REIM(lmi+1-ic(k),j1,1:ic(k),1:3)*xx)
                  REY(1:ic(k),4:8,3)=REIM(lmi+1-ic(k),j0,1:ic(k),8:12)
     &                       *xx1+REIM(lmi+1-ic(k),j1,1:ic(k),8:12)*xx
                  REY(1:ic(k),9,3)=Vsi(lmi+1-ic(k),j0,1:ic(k))*xx1+Vsi(
     &                       lmi+1-ic(k),j1,1:ic(k))*xx
	          ENDIF
	        ELSE
	          ix0=NHX(ic(k),j0)
			  ix1=NHX(ic(k),j1)
                REY(1:NLi,1:3,1)=Dlog(REIH(ix0,k,1:NLi,1:3)*xx1+REIH(ix1
     &                           ,k,1:NLi,1:3)*xx)   
                REY(1:NLi,4:8,1)=REIH(ix0,k,1:NLi,8:12)*xx1+REIH(ix1,k,1
     &                           :NLi,8:12)*xx
               REY(1:NLi,9,1)=VsiH(ix0,k,1:NLi)*xx1+VsiH(ix1,k,1:NLi)*xx
	        ENDIF
            ELSE
              Vx(1,1:2)=Vd(ic(k),j0,k,1:2)*xx1+Vd(ic(k),j1,k,1:2)*xx
	      ENDIF
	      IF((ic(k).GT.NLi+1).or.(i1.GT.NLi+1))THEN
	        aj=MLONC(k)/dmlon(i1)+1.0D0
			j0=INT(aj)
			xx=aj-j0
			IF(j0.gt.Nmlo(i1))j0=j0-Nmlo(i1)
			j1=j0+1
			IF(j1.gt.Nmlo(i1))j1=j1-Nmlo(i1)
			xx1=1.0D0-xx
	      ENDIF
            IF(tre(k).eq.0.0D0)THEN
              IF(i1.LE.NLi)THEN
	          ik=i1
			  IF(k.eq.2)ik=lmi+1-ik
                REY(1:i1,1:3,2)=Dlog(REIM(ik,j0,1:i1,1:3)*xx1+REIM(ik,j1
     &                          ,1:i1,1:3)*xx)   
                REY(1:i1,4:8,2)=REIM(ik,j0,1:i1,8:12)*xx1+REIM(ik,j1,1
     &                          :i1,8:12)*xx
                REY(1:i1,9,2)=Vsi(ik,j0,1:i1)*xx1+Vsi(ik,j1,1:i1)*xx
                IF(io.eq.0)THEN
                  REY(1:i1,1:3,4)=Dlog(REIM(lmi+1-i1,j0,1:i1,1:3)*xx1
     &                            +REIM(lmi+1-i1,j1,1:i1,1:3)*xx) 
                  REY(1:i1,4:8,4)=REIM(lmi+1-i1,j0,1:i1,8:12)*xx1+REIM(
     &                             lmi+1-i1,j1,1:i1,8:12)*xx
                  REY(1:i1,9,4)=Vsi(lmi+1-i1,j0,1:i1)*xx1+Vsi(lmi+1-i1,
     &                             j1,1:i1)*xx
	          ENDIF
	        ELSE
	          ix0=NHX(i1,j0)
			  ix1=NHX(i1,j1)
                REY(1:NLi,1:3,2)=Dlog(REIH(ix0,k,1:NLi,1:3)*xx1+REIH(ix1
     &                           ,k,1:NLi,1:3)*xx) 
                REY(1:NLi,4:8,2)=REIH(ix0,k,1:NLi,8:12)*xx1+REIH(ix1,k,1
     &                           :NLi,8:12)*xx
               REY(1:NLi,9,2)=VsiH(ix0,k,1:NLi)*xx1+VsiH(ix1,k,1:NLi)*xx
	        ENDIF
            ELSE
              Vx(2,1:2)=Vd(i1,j0,k,1:2)*xx1+Vd(i1,j1,k,1:2)*xx
	      ENDIF
            xx=(Xli(ic(k))-MLATC(k))/(Xli(ic(k))-Xli(i1))
		  xx1=1.0D0-xx
            IF(tre(k).eq.0.0D0)THEN
	        i2=imin0(ic(k),i1,NLi)
              REX(1:i2,1:9,k)=REY(1:i2,1:9,1)*xx1+REY(1:i2,1:9,2)*xx
              IF(io.eq.0)THEN
	          REX(1:i2,1:9,2)=REY(1:i2,1:9,3)*xx1+REY(1:i2,1:9,4)*xx
                ih(2)=i2
	        ENDIF
	        IF(i2.LT.NLi)THEN
                i2=i2+1
                REX(i2,1:9,k)=REY(ic(k),1:9,1)*xx1+REY(i1,1:9,2)*xx
		      IF(io.eq.1)REX(i2,9,k)=0.0D0
	        ENDIF
	        ih(k)=i2
            ELSE
              V(k,1:2)=Vx(1,1:2)*xx1+Vx(2,1:2)*xx
	      ENDIF
	    ENDIF
          IF(k.eq.1.and.io.eq.0)THEN
            tre(2)=tre(1)
		  V(2,1:2)=V(1,1:2)
		  ic(2)=ic(1)
		  MLATC(2)=MLATC(1)
		  MLONC(2)=MLONC(1)
          ENDIF
	  ENDDO
	  ENDIF
	ENDDO
          IF(ix.eq.1)THEN
	      DO k=1,2
	        X=0.0D0
			X(1:ih(k))=REX(1:ih(k),9,k)
			IF(k.eq.2)X=0.0D0-X
	        Y=0.0D0
			Y(1:ih(k),1:8)=REX(1:ih(k),1:8,k)
	        T(1:NLi)=VsiH(jx,k,1:NLi)
			IF(k.eq.2)T=0.0D0-T
              CALL DIGMAX(X,Y,lmi+1,8,ih(k),T,Z,lmi+1,8,NLi,8,ij,1)
	        IF(ij.LT.NLi)THEN
	          altv(ij+1:NLi)=(alti(ij+1:NLi)-alti(ij))/(alti(ij)-alti(
     &                         ij-1))
                DO ia=1,3
	            ax=dmin1(0.0D0,z(ij,ia)-z(ij-1,ia))
			    Z(ij+1:NLi,ia)=altv(ij+1:NLi)*ax+z(ij,ia)
                ENDDO
                Z(ij+1:NLi,4)=Z(ij,4)
			  Z(ij+1:NLi,5)=Z(ij,5)
                Z(ij+1:NLi,6:8)=0.0D0
	        ENDIF
              REIHold(jx,k,1:NLi,1)=DEXP(Z(1:NLi,1))
	        IF(iHe.eq.1)REIHold(jx,k,1:NLi,2)=DEXP(Z(1:NLi,2))
              REIHold(jx,k,1:NLi,3)=DEXP(Z(1:NLi,3))
              REIHold(jx,k,1:NLi,4:8)=Z(1:NLi,4:8)
	      ENDDO
	    ELSE
	      IF(io.ne.0)THEN
	        IF(DABS(REX(ih(1),9,1)).LE.1.0D-10.and.DABS(REX(ih(2),9,
     &           2)).LE.1.0D-10)THEN
              REX(ih(1),1:9,1)=(REX(ih(1),1:9,1)+REX(ih(2),1:9,2))*0.5D0
			  ih(2)=ih(2)-1
	        ENDIF
            ENDIF
	      X=0.0D0
		  X(1:ih(1))=REX(1:ih(1),9,1)
		  iht=ih(1)
	      Y=0.0D0
		  Y(1:ih(1),1:8)=REX(1:ih(1),1:8,1)
		  IF(ih(2).ne.0)THEN 
		    iht=ih(1)+ih(2)
			X(ih(1)+1:iht)=REX(ih(2):1:-1,9,2)
	        Y(ih(1)+1:iht,1:8)=REX(ih(2):1:-1,1:8,2)
            ENDIF
            iz=ix
		  T(1:ix)=Vsi(iz,jx,1:ix)
		  iz=lmi+1-iz
		  T(ix+1:2*ix-1)=Vsi(iz,jx,ix-1:1:-1)
            CALL DIGMAX(X,Y,lmi+1,8,iht,T,Z,lmi+1,8,2*ix-1,8,ij,0)
            REIMold(ix,jx,1:ix,1)=DEXP(Z(1:ix,1))
            REIMold(ix,jx,1:ix,3)=DEXP(Z(1:ix,3))
            REIMold(ix,jx,1:ix,4:8)=Z(1:ix,4:8)
            REIMold(iz,jx,1:ix,1)=DEXP(Z(2*ix-1:ix:-1,1))
            REIMold(iz,jx,1:ix,3)=DEXP(Z(2*ix-1:ix:-1,3))
            REIMold(iz,jx,1:ix,4:8)=Z(2*ix-1:ix:-1,4:8)
		  IF(iHe.eq.1)THEN
              REIMold(ix,jx,1:ix,2)=DEXP(Z(1:ix,2))
              REIMold(iz,jx,1:ix,2)=DEXP(Z(2*ix-1:ix:-1,2))
            ENDIF
	    ENDIF
	  ENDDO
	ENDDO
	END SUBROUTINE GridInterpolation
C---------------------------------------------------------------------
      SUBROUTINE DIGMAX(X,Y,M1,M2,N,T,Z,M3,M4,ii,jj,i,io)
C X Y old max M1*M2; T,Z new max M3*M4; use part: X Y N*jj, T Z ii*jj
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION X(M1),Y(M1,M2),T(M3),Z(M3,M4)

        i=1
	  DO WHILE((i.LE.ii).and.(T(i).LE.X(1)))
          Z(i,1:jj)=Y(1,1:jj)
	    i=i+1
	  ENDDO
	  j=2
	  j0=1
	  DO WHILE((i.LE.ii).and.(j.LE.N))
          DO WHILE((T(i).LE.X(j)).and.(i.LE.ii))
            s=(T(i)-X(j0))/(X(j)-X(j0))
		  s1=1.0D0-s
            Z(i,1:jj)=Y(j,1:jj)*s+Y(j0,1:jj)*s1
	      i=i+1
          ENDDO
          j=j+1
	    j0=j0+1
	  ENDDO
	  i=i-1
	  IF(io.eq.0.and.i.LT.ii)THEN
          DO k=i+1,ii
            Z(k,1:jj)=Y(N,1:jj)
          ENDDO
	    i=ii
        ENDIF
      END SUBROUTINE DIGMAX
c---------------------------------------------------------------------
      SUBROUTINE GCMpre2alt(Zg,Tng,cOg,cO2g,cN2g,cHg,cHeg,cNg,Wsg,Weg,
     &                      Wug,Qopg,QHepg,OpLg,Qepg,cO2pg,cNOpg,cN2pg
     &                      ,QHpg,HepLg,HpLg)
      use geogrid_t3d  ,only: nlev,nlat,nlon
      use TIME3Din     ,only: Tn,cO,cO2,cN2,cH,cHe,cN,Ws,We,Wu,Qop,QHep,
     &                        OpL,Qep,cO2p,cNOp,cN2p,QHp,HepL,HpL
      INCLUDE 'EIG.h'
	COMMON/GRIDI/ alti(NLi),Cli(NLi),SIi(NLi),MLATi(lmi+1,NLi),Mloni(
     &              NFi),altH(NHN),ClH(NHN),SIH(NHN),MlatH(NHN,2,NLi),
     &              MlonH(NHN)
      REAL*8 Tng(Nlat,Nlon,Nlev),cOg(Nlat,Nlon,Nlev),cO2g(Nlat,Nlon,Nlev
     &       ),cN2g(Nlat,Nlon,Nlev),cHg(Nlat,Nlon,Nlev),cHeg(Nlat,Nlon,
     &       Nlev),cNg(Nlat,Nlon,Nlev),Wsg(Nlat,Nlon,Nlev),Weg(Nlat,Nlon
     &       ,Nlev),Wug(Nlat,Nlon,Nlev),Qopg(Nlat,Nlon,Nlev),QHepg(Nlat,
     &       Nlon,Nlev),OpLg(Nlat,Nlon,Nlev),Qepg(Nlat,Nlon,Nlev),cO2pg(
     &       Nlat,Nlon,Nlev),cNOpg(Nlat,Nlon,Nlev),cN2pg(Nlat,Nlon,Nlev)
     &       ,Zg(Nlat,Nlon,Nlev),QHpg(Nlat,Nlon,Nlev),HepLg(Nlat,Nlon,
     &       Nlev),HpLg(Nlat,Nlon,Nlev)
      REAL*8 X(nlev),Y(nlev,17),Z(NLi,17),dhh(NLi)

      DO i=1,Nlat
        Do j=1,Nlon
	    X(1:nlev)=Zg(i,j,1:nlev)
          Y(1:nlev, 1)=Tng(i,j,1:nlev)
          Y(1:nlev, 2)=Wsg(i,j,1:nlev)
          Y(1:nlev, 3)=Weg(i,j,1:nlev)
          Y(1:nlev, 4)=Wug(i,j,1:nlev)
          Y(1:nlev, 5)=dlog( cOg(i,j,1:nlev))
          Y(1:nlev, 6)=dlog(cO2g(i,j,1:nlev))
          Y(1:nlev, 7)=dlog(cN2g(i,j,1:nlev))
          Y(1:nlev, 8)=dlog( cHg(i,j,1:nlev))
          Y(1:nlev, 9)=dlog( cNg(i,j,1:nlev))
          Y(1:nlev,10)=dlog(Qopg(i,j,1:nlev))
          Y(1:nlev,11)=dlog(OpLg(i,j,1:nlev))
          Y(1:nlev,12)=dlog(Qepg(i,j,1:nlev))
          Y(1:nlev,13)=dlog(cO2pg(i,j,1:nlev))
          Y(1:nlev,14)=dlog(cNOpg(i,j,1:nlev))
          Y(1:nlev,15)=dlog(cN2pg(i,j,1:nlev))
          Y(1:nlev,16)=dlog(QHpg(i,j,1:nlev))
          Y(1:nlev,17)=dlog(HpLg(i,j,1:nlev))
	    j0=17
	    jj=j0
          IF(iHe.eq.1)THEN
            Y(1:nlev,j0+1)=dlog(cHeg(i,j,1:lev))
            Y(1:nlev,j0+2)=dlog(QHepg(i,j,1:lev))
            Y(1:nlev,j0+3)=dlog(HepLg(i,j,1:nlev))
	      jj=j0+3
          ENDIF
          k=1
          DO WHILE((k.LE.NLi).and.(ALTi(k).LE.X(1)))
            Z(k,1:jj)=Y(1,1:jj)
	      k=k+1
	    ENDDO
	    l=2
	    l0=1
	    DO WHILE((k.LE.NLi).and.(l.LE.nlev))
            DO WHILE((ALTi(k).LE.X(l)).and.(k.LE.NLi))
              s=(ALTi(k)-X(l0))/(X(l)-X(l0))
		    s1=1.0D0-s
              Z(k,1:jj)=Y(l,1:jj)*s+Y(l0,1:jj)*s1
	        k=k+1
            ENDDO
            l=l+1
	      l0=l0+1
	    ENDDO
	    kh=k-1
	    IF(kh.LT.NLi)THEN
            Z(k:NLi, 1)=Tng(i,j,nlev)
            Z(k:NLi, 2)=Wsg(i,j,nlev)
            Z(k:NLi, 3)=Weg(i,j,nlev)
            Z(k:NLi, 4)=Wug(i,j,nlev)
            dco =(Y(nlev,5)-Y(nlev-5,5))/(X(nlev)-X(nlev-5))
            dco2=(Y(nlev,6)-Y(nlev-5,6))/(X(nlev)-X(nlev-5))
            dcn2=(Y(nlev,7)-Y(nlev-5,7))/(X(nlev)-X(nlev-5))
            dch =(Y(nlev,8)-Y(nlev-5,8))/(X(nlev)-X(nlev-5))
            dcn =(Y(nlev,9)-Y(nlev-5,9))/(X(nlev)-X(nlev-5))
            QopO=Qopg(i,j,nlev)/cOg(i,j,nlev)
            OpLN2=OpLg(i,j,nlev)/cN2g(i,j,nlev)
            QepO=Qepg(i,j,nlev)/cOg(i,j,nlev)
            cO2PO2=cO2pg(i,j,nlev)/cO2g(i,j,nlev)
            cNOPN2=cNOpg(i,j,nlev)/cN2g(i,j,nlev)
            cN2PN2=cN2pg(i,j,nlev)/cN2g(i,j,nlev)
            QHpH=QHpg(i,j,nlev)/cHg(i,j,nlev)
            HpLO=HpLg(i,j,nlev)/cOg(i,j,nlev)
            dhh(k:NLi)=Alti(k:NLi)-X(nlev)
            Z(k:NLi, 5)=Y(nlev, 5) + dco*dhh(k:NLi)
            Z(k:NLi, 6)=Y(nlev, 6) +dco2*dhh(k:NLi)
            Z(k:NLi, 7)=Y(nlev, 7) +dcn2*dhh(k:NLi)
            Z(k:NLi, 8)=Y(nlev, 8) + dch*dhh(k:NLi)
            Z(k:NLi, 9)=Y(nlev, 9) + dcn*dhh(k:NLi)
            IF(iHe.eq.1)THEN
              dcHe=(Y(nlev,j0+1)-Y(nlev-5,j0+1))/(X(nlev)-X(nlev-5))
              Z(k:NLi,j0+1)=Y(nlev,j0+1)+dcHe*dhh(k:NLi)
              QHepHe=QHepg(i,j,nlev)/cHeg(i,j,nlev)
              HepLN2=HepLg(i,j,nlev)/cN2g(i,j,nlev)
            ENDIF
          ENDIF
          Tn(i,j,1:NLi)=Z(1:NLi, 1)
          Ws(i,j,1:NLi)=Z(1:NLi, 2)
          We(i,j,1:NLi)=Z(1:NLi, 3)
          Wu(i,j,1:NLi)=Z(1:NLi, 4)
          cO(i,j,1:NLi)  =dmax1(1.0D-25,dexp(Z(1:NLi, 5)))
          cO2(i,j,1:NLi) =dmax1(1.0D-25,dexp(Z(1:NLi, 6)))
          cN2(i,j,1:NLi) =dmax1(1.0D-25,dexp(Z(1:NLi, 7)))
          cH(i,j,1:NLi)  =dmax1(1.0D-25,dexp(Z(1:NLi, 8)))
          cN(i,j,1:NLi)  =dmax1(1.0D-25,dexp(Z(1:NLi, 9)))
          Qop(i,j,1:kh)  =dexp(Z(1:kh, 10))
          OpL(i,j,1:kh)  =dexp(Z(1:kh, 11))
          Qep(i,j,1:kh)  =dexp(Z(1:kh, 12))
          cO2p(i,j,1:kh) =dmax1(1.0D-25,dexp(Z(1:kh, 13)))
          cNOp(i,j,1:kh) =dmax1(1.0D-25,dexp(Z(1:kh, 14)))
          cN2p(i,j,1:kh) =dmax1(1.0D-25,dexp(Z(1:kh, 15)))
          QHp(i,j,1:kh)  =dexp(Z(1:kh, 16))
          HpL(i,j,1:kh)  =dexp(Z(1:kh, 17))
	    IF(kh.LT.NLi)THEN
            Qop(i,j,k:NLi)=QopO*  cO(i,j,k:NLi)
            OpL(i,j,k:NLi)=OpLN2*cN2(i,j,k:NLi)
            Qep(i,j,k:NLi)=QepO*  cO(i,j,k:NLi)
            cO2p(i,j,k:NLi)=dmax1(1.0D-25,cO2PO2*cO2(i,j,k:NLi))
            cNOp(i,j,k:NLi)=dmax1(1.0D-25,cNOPN2*cN2(i,j,k:NLi))
            cN2p(i,j,k:NLi)=dmax1(1.0D-25,cN2PN2*cN2(i,j,k:NLi))
            QHp(i,j,k:NLi)=QHpH*  cH(i,j,k:NLi)
            HpL(i,j,k:NLi)=HpLO  *cO(i,j,k:NLi)
          ENDIF
          IF(iHe.eq.1)THEN
            cHe(i,j,1:NLi)=dmax1(1.0D-25,dexp(Z(1:NLi,j0+1)))
            QHep(i,j,1:kh)=dexp(Z(1:kh,j0+2))
            HepL(i,j,1:kh) =dexp(Z(1:kh, j0+3))
            IF(kh.LT.NLi)THEN
		    QHep(i,j,k:NLi)=QHepHe*cHe(i,j,k:NLi)
		    HepL(i,j,k:NLi)=HepLN2*cN2(i,j,k:NLi)
            ENDIF
          ENDIF
	  ENDDO
	ENDDO
      END SUBROUTINE GCMpre2alt
c---------------------------------------------------------------------
      SUBROUTINE GCM2TIME3D()    !QHp,HepL,HpL,
      use geogrid_t3d  ,only: nlev,nlat,nlon,ylatg,ylong
      use InterpTI     ,only: XMLAI,XMLOI,Cpole,xij,xijh,xix,ij,ijh,ix
      use TIME3Din     ,only: Tn,cO,cO2,cN2,cH,cHe,cN,Ws,We,Wu,Qop,QHep,
     &                        OpL,Qep,cO2p,cNOp,cN2p,QHp,HepL,HpL
      INCLUDE 'EIG.h'
	COMMON/REIM3/ REIM(lmi,NFi,NLi,12),REIMold(lmi,NFi,NLi,8),REIH(NHN
     &             ,2,NLi,12),REIHold(NHN,2,NLi,8)
	COMMON/REIM3N/ PN(lmi,NFi,NLi,8),PLI(lmi,NFi,NLi,8),PNH(NHN,2,NLi,
     &               8),PLIH(NHN,2,NLi,8)
	COMMON/DynamoWind/ Ud(lmi,NFi,NLi,2),Udh(NHN,2,NLi,2)
	COMMON/EEmC/  EEm(NLi,NFi,2),EEH(NHN,2,2)
	COMMON/IGrid/ d1i(lmi,NFi,NLi,3),d2i(lmi,NFi,NLi,3),glati(lmi,NFi,
     &              NLi),gloni(lmi,NFi,NLi),Vsi(lmi,NFi,NLi),Wsi(lmi,NFi
     &              ,NLi),Bsi(lmi,NFi,NLi,4),gxi(lmi,NFi,NLi,3),dvv(lmi,
     &              NFi,NLi,2),mlov(NLi,NFi,2),d1ih(NHN,2,NLi,3),d2ih(
     &              NHN,2,NLi,3),glatih(NHN,2,NLi),glonih(NHN,2,NLi),
     &              Vsih(NHN,2,NLi),Wsih(NHN,2,NLi),Bsih(NHN,2,NLi,4),
     &              gxih(NHN,2,NLi,3),dvvh(NHN,2,NLi,2),mlovh(NHN,2,2),
     &              MH300(NHN+NLi)
      REAL*8 BB(3),Wind(3),Pxy(NLi,4),Wpole(2)
      
	dtr=3.1415926D0/180.0D0
	ipole=0
c mid- and low-latitude
	DO j=1,NFi
	  DO k=1,NLi
	    DO i=k,lmi-k
          xi1=xij(i,j,k,1)
          xi2=xij(i,j,k,2)
          xi3=xij(i,j,k,3)
          xi4=xij(i,j,k,4)
          i1=ij(i,j,k,1)
          i2=ij(i,j,k,2)
          i3=ij(i,j,k,3)
          i4=ij(i,j,k,4)
	    IF(i2.ne.0)THEN
            PN(i,j,k,1)=Tn(i1,i3,k)*xi1  +Tn(i1,i4,k)*xi2
     &                 +Tn(i2,i3,k)*xi3  +Tn(i2,i4,k)*xi4
            PN(i,j,k,2)=cO(i1,i3,k)*xi1  +cO(i1,i4,k)*xi2 
     &                 +cO(i2,i3,k)*xi3  +cO(i2,i4,k)*xi4
            PN(i,j,k,3)=cO2(i1,i3,k)*xi1 +cO2(i1,i4,k)*xi2 
     &                 +cO2(i2,i3,k)*xi3 +cO2(i2,i4,k)*xi4
            PN(i,j,k,4)=cN2(i1,i3,k)*xi1 +cN2(i1,i4,k)*xi2 
     &                 +cN2(i2,i3,k)*xi3 +cN2(i2,i4,k)*xi4
            PN(i,j,k,5)=CHe(i1,i3,k)*xi1 +cHe(i1,i4,k)*xi2  
     &                 +cHe(i2,i3,k)*xi3 +cHe(i2,i4,k)*xi4
            PN(i,j,k,6)=cH(i1,i3,k)*xi1  +cH(i1,i4,k)*xi2  
     &                 +cH(i2,i3,k)*xi3  +cH(i2,i4,k)*xi4
            PN(i,j,k,7)=cN(i1,i3,k)*xi1  +cN(i1,i4,k)*xi2  
     &                 +cN(i2,i3,k)*xi3  +cN(i2,i4,k)*xi4
            Wind(1)=Ws(i1,i3,k)*xi1 +Ws(i1,i4,k)*xi2
     &             +Ws(i2,i3,k)*xi3 +Ws(i2,i4,k)*xi4
            Wind(2)=We(i1,i3,k)*xi1 +We(i1,i4,k)*xi2  
     &             +We(i2,i3,k)*xi3 +We(i2,i4,k)*xi4
            Wind(3)=Wu(i1,i3,k)*xi1 +Wu(i1,i4,k)*xi2 
     &             +Wu(i2,i3,k)*xi3 +Wu(i2,i4,k)*xi4
            PN(i,j,k,8)=-(Wind(1)*BSi(i,j,k,1)+Wind(2)*BSi(i,j,k,2)+
     &                  Wind(3)*BSi(i,j,k,3))/BSi(i,j,k,4)    !-B*W/|B|
            Ud(i,j,k,1)=Wind(1)*d1i(i,j,k,2)+Wind(2)*d1i(i,j,k,3)
     &                 +Wind(3)*d1i(i,j,k,1)
            Ud(i,j,k,2)=Wind(1)*d2i(i,j,k,2)+Wind(2)*d2i(i,j,k,3)
     &                 +Wind(3)*d2i(i,j,k,1)
            BB(1:3)=EEM(k,j,1)*d1i(i,j,k,1:3)+EEM(k,j,2)*d2i(i,j,k,1:3) !E-U¡ÁB
            BB(1)=BB(1)+Wind(2)*BSi(i,j,k,1)-Wind(1)*BSi(i,j,k,2)
	      BB(2)=BB(2)-Wind(2)*BSi(i,j,k,3)+Wind(3)*BSi(i,j,k,2)
		  BB(3)=BB(3)+Wind(1)*BSi(i,j,k,3)-Wind(3)*BSi(i,j,k,1)
           PLI(i,j,k,8)=sum(BB(1:3)*BB(1:3))/(BSi(i,j,k,4)*BSi(i,j,k,4))
            PLI(i,j,k,1)=QOp(i1,i3,k)*xi1 +QOp(i1,i4,k)*xi2   
     &                  +QOp(i2,i3,k)*xi3 +QOp(i2,i4,k)*xi4 
            PLI(i,j,k,3)=QHp(i1,i3,k)*xi1 +QHp(i1,i4,k)*xi2  
     &                  +QHp(i2,i3,k)*xi3 +QHp(i2,i4,k)*xi4 
            PLI(i,j,k,4)=OpL(i1,i3,k)*xi1 +OpL(i1,i4,k)*xi2  
     &                  +OpL(i2,i3,k)*xi3 +OpL(i2,i4,k)*xi4 
            PLI(i,j,k,6)=HpL(i1,i3,k)*xi1 +HpL(i1,i4,k)*xi2 
     &                  +HpL(i2,i3,k)*xi3 +HpL(i2,i4,k)*xi4
            PLI(i,j,k,7)=Qep(i1,i3,k)*xi1 +Qep(i1,i4,k)*xi2  
     &                  +Qep(i2,i3,k)*xi3 +Qep(i2,i4,k)*xi4 
	      PLI(i,j,k,2)=0.0D0
            PLI(i,j,k,5)=0.0D0
            IF(iHe.eq.1)THEN
              PLI(i,j,k,2)=QHep(i1,i3,k)*xi1 +QHep(i1,i4,k)*xi2 
     &                    +QHep(i2,i3,k)*xi3 +QHep(i2,i4,k)*xi4
		    PLI(i,j,k,5)=HepL(i1,i3,k)*xi1 +HepL(i1,i4,k)*xi2 
     &                    +HepL(i2,i3,k)*xi3 +HepL(i2,i4,k)*xi4
            ENDIF
            REIM(i,j,k,5)=cNOp(i1,i3,k)*xi1 +cNOp(i1,i4,k)*xi2  
     &                   +cNOp(i2,i3,k)*xi3 +cNOp(i2,i4,k)*xi4 
            REIM(i,j,k,4)=cO2p(i1,i3,k)*xi1 +cO2p(i1,i4,k)*xi2  
     &                   +cO2p(i2,i3,k)*xi3 +cO2p(i2,i4,k)*xi4 
            REIM(i,j,k,6)=cN2p(i1,i3,k)*xi1 +cN2p(i1,i4,k)*xi2  
     &                   +cN2p(i2,i3,k)*xi3 +cN2p(i2,i4,k)*xi4 
          ELSE
            PN(i,j,k,1)=Tn(i1,i3,k)*xi1  +Tn(i1,i4,k)*xi2  
     &                 +sum(Tn(i1,1:Nlon,k))/dble(Nlon)*xi3
            PN(i,j,k,2)=cO(i1,i3,k)*xi1  +cO(i1,i4,k)*xi2 
     &                 +sum(cO(i1,1:Nlon,k))/dble(Nlon)*xi3
            PN(i,j,k,3)=cO2(i1,i3,k)*xi1 +cO2(i1,i4,k)*xi2  
     &                 +sum(cO2(i1,1:Nlon,k))/dble(Nlon)*xi3
            PN(i,j,k,4)=cN2(i1,i3,k)*xi1 +cN2(i1,i4,k)*xi2 
     &                 +sum(cN2(i1,1:Nlon,k))/dble(Nlon)*xi3
            PN(i,j,k,5)=CHe(i1,i3,k)*xi1 +cHe(i1,i4,k)*xi2 
     &                 +sum(cHe(i1,1:Nlon,k))/dble(Nlon)*xi3
            PN(i,j,k,6)=cH(i1,i3,k)*xi1  +cH(i1,i4,k)*xi2 
     &                 +sum(cH(i1,1:Nlon,k))/dble(Nlon)*xi3
            PN(i,j,k,7)=cN(i1,i3,k)*xi1  +cN(i1,i4,k)*xi2 
     &                 +sum(cN(i1,1:Nlon,k))/dble(Nlon)*xi3
            IF(ipole.ne.1)THEN
	        ipole=1
	        DO ji=1,NLi
                Pxy(ji,1)=sum(Ws(1,1:NLon,k)*Cpole(1:NLon,1,1)+We(1,1:
     &                   NLon,k)*Cpole(1:NLon,1,2))/dble(NLon)
                Pxy(ji,2)=sum(Ws(1,1:NLon,k)*Cpole(1:NLon,1,3)+We(1,1:
     &                   NLon,k)*Cpole(1:NLon,1,4))/dble(NLon)
                Pxy(ji,3)=sum(Ws(Nlat,1:NLon,k)*Cpole(1:NLon,2,1)+We(
     &                   Nlat,1:NLon,k)*Cpole(1:NLon,2,2))/dble(NLon)
                Pxy(ji,4)=sum(Ws(Nlat,1:NLon,k)*Cpole(1:NLon,2,3)+We(
     &                   Nlat,1:NLon,k)*Cpole(1:NLon,2,4))/dble(NLon)
	        ENDDO
            ENDIF
            IF(i1.eq.1)THEN
	        Wpole(1)=Pxy(k,1)*dcos(gloni(i,j,k)*dtr-ylong(1))
     &                +Pxy(k,2)*dsin(gloni(i,j,k)*dtr-ylong(1))
	        Wpole(2)=Pxy(k,2)*dcos(gloni(i,j,k)*dtr-ylong(1))
     &	        	-Pxy(k,1)*dsin(gloni(i,j,k)*dtr-ylong(1))
		  ELSE
	        Wpole(1)=Pxy(k,3)*dcos(gloni(i,j,k)*dtr-ylong(1))
     &                -Pxy(k,4)*dsin(gloni(i,j,k)*dtr-ylong(1))
	        Wpole(2)=Pxy(k,4)*dcos(gloni(i,j,k)*dtr-ylong(1))
     &	        	+Pxy(k,3)*dsin(gloni(i,j,k)*dtr-ylong(1))
            ENDIF
            Wind(1)=Ws(i1,i3,k)*xi1 +Ws(i1,i4,k)*xi2 +Wpole(1)*xi3  
            Wind(2)=We(i1,i3,k)*xi1 +We(i1,i4,k)*xi2 +Wpole(2)*xi3  
            Wind(3)=Wu(i1,i3,k)*xi1 +Wu(i1,i4,k)*xi2  
     &             +sum(Wu(i1,1:Nlon,k))/dble(Nlon)*xi3
            PN(i,j,k,8)=-(Wind(1)*BSi(i,j,k,1)+Wind(2)*BSi(i,j,k,2)+
     &                   Wind(3)*BSi(i,j,k,3))/BSi(i,j,k,4)    !-B*W/|B|
            Ud(i,j,k,1)=Wind(1)*d1i(i,j,k,2)+Wind(2)*d1i(i,j,k,3)
     &                 +Wind(3)*d1i(i,j,k,1)
            Ud(i,j,k,2)=Wind(1)*d2i(i,j,k,2)+Wind(2)*d2i(i,j,k,3)
     &                 +Wind(3)*d2i(i,j,k,1)
            BB(1:3)=EEM(k,j,1)*d1i(i,j,k,1:3)+EEM(k,j,2)*d2i(i,j,k,1:3) !E-U¡ÁB
            BB(1)=BB(1)+Wind(2)*BSi(i,j,k,1)-Wind(1)*BSi(i,j,k,2)
	      BB(2)=BB(2)-Wind(2)*BSi(i,j,k,3)+Wind(3)*BSi(i,j,k,2)
		  BB(3)=BB(3)+Wind(1)*BSi(i,j,k,3)-Wind(3)*BSi(i,j,k,1)
           PLI(i,j,k,8)=sum(BB(1:3)*BB(1:3))/(BSi(i,j,k,4)*BSi(i,j,k,4))
            PLI(i,j,k,1)=QOp(i1,i3,k)*xi1 +QOp(i1,i4,k)*xi2  
     &                  +sum(QOp(i1,1:Nlon,k))/dble(Nlon)*xi3
            PLI(i,j,k,3)=QHp(i1,i3,k)*xi1 +QHp(i1,i4,k)*xi2 
     &                  +sum(QHp(i1,1:Nlon,k))/dble(Nlon)*xi3
            PLI(i,j,k,4)=OpL(i1,i3,k)*xi1 +OpL(i1,i4,k)*xi2  
     &                  +sum(OpL(i1,1:Nlon,k))/dble(Nlon)*xi3
            PLI(i,j,k,6)=HpL(i1,i3,k)*xi1 +HpL(i1,i4,k)*xi2  
     &                  +sum(HpL(i1,1:Nlon,k))/dble(Nlon)*xi3
            PLI(i,j,k,7)=Qep(i1,i3,k)*xi1 +Qep(i1,i4,k)*xi2  
     &                  +sum(Qep(i1,1:Nlon,k))/dble(Nlon)*xi3
	      PLI(i,j,k,2)=0.0D0
            PLI(i,j,k,5)=0.0D0
            IF(iHe.eq.1)THEN
              PLI(i,j,k,2)=QHep(i1,i3,k)*xi1 +QHep(i1,i4,k)*xi2  
     &                    +sum(QHep(i1,1:Nlon,k))/dble(Nlon)*xi3
		    PLI(i,j,k,5)=HepL(i1,i3,k)*xi1 +HepL(i1,i4,k)*xi2 
     &                    +sum(HepL(i1,1:Nlon,k))/dble(Nlon)*xi3
            ENDIF
            REIM(i,j,k,5)=cNOp(i1,i3,k)*xi1 +cNOp(i1,i4,k)*xi2 
     &                   +sum(cNOp(i1,1:Nlon,k))/dble(Nlon)*xi3 
            REIM(i,j,k,4)=cO2p(i1,i3,k)*xi1 +cO2p(i1,i4,k)*xi2 
     &                   +sum(cO2p(i1,1:Nlon,k))/dble(Nlon)*xi3
            REIM(i,j,k,6)=cN2p(i1,i3,k)*xi1 +cN2p(i1,i4,k)*xi2  
     &                   +sum(cN2p(i1,1:Nlon,k))/dble(Nlon)*xi3
	    ENDIF
            OpHr=2.5D-17*dsqrt(PN(i,j,k,1))*PN(i,j,k,6)             !O+ + H --> H+ +O
            HpOr=2.2D-17*dsqrt(REIM(i,j,k,9))*PN(i,j,k,2)           !H+ + O --> O+ +H
	      per=4.43D-18*(REIM(i,j,k,8)/3.0D2)**0.7D0*REIM(i,j,k,7) !O+ (He+, H+) +e --> O (He, H)
            PLI(i,j,k,4)=PLI(i,j,k,4)  +OpHr    +per
            PLI(i,j,k,6)=PLI(i,j,k,6)  +HpOr    +per
            PLI(i,j,k,1)=PLI(i,j,k,1)  +HpOr*REIM(i,j,k,3)
            PLI(i,j,k,3)=PLI(i,j,k,3)  +OpHr*REIM(i,j,k,1)
            IF(iHe.eq.1)THEN
              HepO2r=8.0d-16*PN(i,j,k,3)                            !He+ + O2 --> He +O + O+
              PLI(i,j,k,5)=PLI(i,j,k,5)  +HepO2r  +per
              PLI(i,j,k,1)=PLI(i,j,k,1)  +HepO2r*REIM(i,j,k,2)
            ENDIF
	    ENDDO
	    ki=lmi-k+1
	    PN(ki,j,k,1:8)=PN(k,j,k,1:8)
	    PLI(ki,j,k,1:8)=PLI(k,j,k,1:8)
	    REIM(ki,j,k,4:6)=REIM(k,j,k,4:6)
	    Ud(ki,j,k,1:2)=Ud(k,j,k,1:2)
	  ENDDO
	ENDDO
c high-latitude
	DO j=1,2
	  DO k=1,NLi
	    DO i=1,NHN
          xi1=xijh(i,j,k,1)
          xi2=xijh(i,j,k,2)
          xi3=xijh(i,j,k,3)
          xi4=xijh(i,j,k,4)
          i1=ijh(i,j,k,1)
          i2=ijh(i,j,k,2)
          i3=ijh(i,j,k,3)
          i4=ijh(i,j,k,4)
	    IF(i2.ne.0)THEN
            PNH(i,j,k,1)=Tn(i1,i3,k)*xi1  +Tn(i1,i4,k)*xi2  
     &                  +Tn(i2,i3,k)*xi3  +Tn(i2,i4,k)*xi4
            PNH(i,j,k,2)=cO(i1,i3,k)*xi1  +cO(i1,i4,k)*xi2 
     &                  +cO(i2,i3,k)*xi3  +cO(i2,i4,k)*xi4
            PNH(i,j,k,3)=cO2(i1,i3,k)*xi1 +cO2(i1,i4,k)*xi2 
     &                  +cO2(i2,i3,k)*xi3 +cO2(i2,i4,k)*xi4
            PNH(i,j,k,4)=cN2(i1,i3,k)*xi1 +cN2(i1,i4,k)*xi2 
     &                  +cN2(i2,i3,k)*xi3 +cN2(i2,i4,k)*xi4
            PNH(i,j,k,5)=CHe(i1,i3,k)*xi1 +cHe(i1,i4,k)*xi2
     &                  +cHe(i2,i3,k)*xi3 +cHe(i2,i4,k)*xi4   
            PNH(i,j,k,6)=cH(i1,i3,k)*xi1  +cH(i1,i4,k)*xi2 
     &                  +cH(i2,i3,k)*xi3  +cH(i2,i4,k)*xi4
            PNH(i,j,k,7)=cN(i1,i3,k)*xi1  +cN(i1,i4,k)*xi2 
     &                  +cN(i2,i3,k)*xi3  +cN(i2,i4,k)*xi4
            Wind(1)=Ws(i1,i3,k)*xi1 +Ws(i1,i4,k)*xi2 
     &             +Ws(i2,i3,k)*xi3 +Ws(i2,i4,k)*xi4
            Wind(2)=We(i1,i3,k)*xi1 +We(i1,i4,k)*xi2 
     &             +We(i2,i3,k)*xi3 +We(i2,i4,k)*xi4
            Wind(3)=Wu(i1,i3,k)*xi1 +Wu(i1,i4,k)*xi2 
     &             +Wu(i2,i3,k)*xi3 +Wu(i2,i4,k)*xi4
            PNH(i,j,k,8)=-(Wind(1)*BSiH(i,j,k,1)+Wind(2)*BSiH(i,j,k,2)+
     &                   Wind(3)*BSiH(i,j,k,3))/BSiH(i,j,k,4)    !-B*W/|B|
            Udh(i,j,k,1)=Wind(1)*d1ih(i,j,k,2)+Wind(2)*d1ih(i,j,k,3)
     &                  +Wind(3)*d1ih(i,j,k,1)
            Udh(i,j,k,2)=Wind(1)*d2ih(i,j,k,2)+Wind(2)*d2ih(i,j,k,3)
     &                  +Wind(3)*d2ih(i,j,k,1)
           BB(1:3)=EEH(k,j,1)*d1iH(i,j,k,1:3)+EEH(k,j,2)*d2iH(i,j,k,1:3) !E-U¡ÁB
            BB(1)=BB(1)+Wind(2)*BSiH(i,j,k,1)-Wind(1)*BSiH(i,j,k,2)
	      BB(2)=BB(2)-Wind(2)*BSiH(i,j,k,3)+Wind(3)*BSiH(i,j,k,2)
		  BB(3)=BB(3)+Wind(1)*BSiH(i,j,k,3)-Wind(3)*BSiH(i,j,k,1)
        PLIH(i,j,k,8)=sum(BB(1:3)*BB(1:3))/(BSiH(i,j,k,4)*BSiH(i,j,k,4))
            PLIH(i,j,k,1)=QOp(i1,i3,k)*xi1 +QOp(i1,i4,k)*xi2 
     &                   +QOp(i2,i3,k)*xi3 +QOp(i2,i4,k)*xi4
            PLIH(i,j,k,3)=QHp(i1,i3,k)*xi1 +QHp(i1,i4,k)*xi2 
     &                   +QHp(i2,i3,k)*xi3 +QHp(i2,i4,k)*xi4
            PLIH(i,j,k,4)=OpL(i1,i3,k)*xi1 +OpL(i1,i4,k)*xi2
     &                   +OpL(i2,i3,k)*xi3 +OpL(i2,i4,k)*xi4
            PLIH(i,j,k,6)=HpL(i1,i3,k)*xi1 +HpL(i1,i4,k)*xi2 
     &                   +HpL(i2,i3,k)*xi3 +HpL(i2,i4,k)*xi4
            PLIH(i,j,k,7)=Qep(i1,i3,k)*xi1 +Qep(i1,i4,k)*xi2 
     &                   +Qep(i2,i3,k)*xi3 +Qep(i2,i4,k)*xi4
	      PLIH(i,j,k,2)=0.0D0
            PLIH(i,j,k,5)=0.0D0
            IF(iHe.eq.1)THEN
              PLIH(i,j,k,2)=QHep(i1,i3,k)*xi1 +QHep(i1,i4,k)*xi2
     &                     +QHep(i2,i3,k)*xi3 +QHep(i2,i4,k)*xi4
		    PLIH(i,j,k,5)=HepL(i1,i3,k)*xi1 +HepL(i1,i4,k)*xi2
     &                     +HepL(i2,i3,k)*xi3 +HepL(i2,i4,k)*xi4
            ENDIF
            REIH(i,j,k,5)=cNOp(i1,i3,k)*xi1 +cNOp(i1,i4,k)*xi2 
     &                   +cNOp(i2,i3,k)*xi3 +cNOp(i2,i4,k)*xi4
            REIH(i,j,k,4)=cO2p(i1,i3,k)*xi1 +cO2p(i1,i4,k)*xi2 
     &                   +cO2p(i2,i3,k)*xi3 +cO2p(i2,i4,k)*xi4
            REIH(i,j,k,6)=cN2p(i1,i3,k)*xi1 +cN2p(i1,i4,k)*xi2 
     &                   +cN2p(i2,i3,k)*xi3 +cN2p(i2,i4,k)*xi4 
          ELSE
            PNH(i,j,k,1)=Tn(i1,i3,k)*xi1  +Tn(i1,i4,k)*xi2 
     &                  +sum(Tn(i1,1:Nlon,k))/dble(Nlon)*xi3
            PNH(i,j,k,2)=cO(i1,i3,k)*xi1  +cO(i1,i4,k)*xi2 
     &                  +sum(cO(i1,1:Nlon,k))/dble(Nlon)*xi3
            PNH(i,j,k,3)=cO2(i1,i3,k)*xi1 +cO2(i1,i4,k)*xi2 
     &                  +sum(cO2(i1,1:Nlon,k))/dble(Nlon)*xi3
            PNH(i,j,k,4)=cN2(i1,i3,k)*xi1 +cN2(i1,i4,k)*xi2 
     &                  +sum(cN2(i1,1:Nlon,k))/dble(Nlon)*xi3
            PNH(i,j,k,5)=CHe(i1,i3,k)*xi1 +cHe(i1,i4,k)*xi2  
     &                  +sum(cHe(i1,1:Nlon,k))/dble(Nlon)*xi3
            PNH(i,j,k,6)=cH(i1,i3,k)*xi1  +cH(i1,i4,k)*xi2 
     &                  +sum(cH(i1,1:Nlon,k))/dble(Nlon)*xi3
            PNH(i,j,k,7)=cN(i1,i3,k)*xi1  +cN(i1,i4,k)*xi2 
     &                  +sum(cN(i1,1:Nlon,k))/dble(Nlon)*xi3
            IF(ipole.ne.1)THEN
	        ipole=1
	        DO ji=1,NLi
                Pxy(ji,1)=sum(Ws(1,1:NLon,k)*Cpole(1:NLon,1,1)+We(1,1:
     &                   NLon,k)*Cpole(1:NLon,1,2))/dble(NLon)
                Pxy(ji,2)=sum(Ws(1,1:NLon,k)*Cpole(1:NLon,1,3)+We(1,1:
     &                   NLon,k)*Cpole(1:NLon,1,4))/dble(NLon)
                Pxy(ji,3)=sum(Ws(Nlat,1:NLon,k)*Cpole(1:NLon,2,1)+We(
     &                   Nlat,1:NLon,k)*Cpole(1:NLon,2,2))/dble(NLon)
                Pxy(ji,4)=sum(Ws(Nlat,1:NLon,k)*Cpole(1:NLon,2,3)+We(
     &                   Nlat,1:NLon,k)*Cpole(1:NLon,2,4))/dble(NLon)
	        ENDDO
            ENDIF
            IF(i1.eq.1)THEN
	        Wpole(1)=Pxy(k,1)*dcos(glonih(i,j,k)*dtr-ylong(1))
     &                +Pxy(k,2)*dsin(glonih(i,j,k)*dtr-ylong(1))
	        Wpole(2)=Pxy(k,2)*dcos(glonih(i,j,k)*dtr-ylong(1))
     &	        	-Pxy(k,1)*dsin(glonih(i,j,k)*dtr-ylong(1))
		  ELSE
	        Wpole(1)=Pxy(k,3)*dcos(glonih(i,j,k)*dtr-ylong(1))
     &                -Pxy(k,4)*dsin(glonih(i,j,k)*dtr-ylong(1))
	        Wpole(2)=Pxy(k,4)*dcos(glonih(i,j,k)*dtr-ylong(1))
     &	        	+Pxy(k,3)*dsin(glonih(i,j,k)*dtr-ylong(1))
            ENDIF
            Wind(1)=Ws(i1,i3,k)*xi1 +Ws(i1,i4,k)*xi2 +Wpole(1)*xi3
            Wind(2)=We(i1,i3,k)*xi1 +We(i1,i4,k)*xi2 +Wpole(2)*xi3
            Wind(3)=Wu(i1,i3,k)*xi1 +Wu(i1,i4,k)*xi2 
     &             +sum(Wu(i1,1:Nlon,k))/dble(Nlon)*xi3
            PNH(i,j,k,8)=-(Wind(1)*BSiH(i,j,k,1)+Wind(2)*BSiH(i,j,k,2)+
     &                   Wind(3)*BSiH(i,j,k,3))/BSiH(i,j,k,4)    !-B*W/|B|
            Udh(i,j,k,1)=Wind(1)*d1ih(i,j,k,2)+Wind(2)*d1ih(i,j,k,3)
     &                  +Wind(3)*d1ih(i,j,k,1)
            Udh(i,j,k,2)=Wind(1)*d2ih(i,j,k,2)+Wind(2)*d2ih(i,j,k,3)
     &                  +Wind(3)*d2ih(i,j,k,1)
           BB(1:3)=EEH(k,j,1)*d1iH(i,j,k,1:3)+EEH(k,j,2)*d2iH(i,j,k,1:3) !E-U¡ÁB
            BB(1)=BB(1)+Wind(2)*BSiH(i,j,k,1)-Wind(1)*BSiH(i,j,k,2)
	      BB(2)=BB(2)-Wind(2)*BSiH(i,j,k,3)+Wind(3)*BSiH(i,j,k,2)
		  BB(3)=BB(3)+Wind(1)*BSiH(i,j,k,3)-Wind(3)*BSiH(i,j,k,1)
            PLIH(i,j,k,8)=sum(BB(1:3)*BB(1:3))/(BSiH(i,j,k,4)*BSiH(i,j,k
     &                    ,4))
            PLIH(i,j,k,1)=QOp(i1,i3,k)*xi1 +QOp(i1,i4,k)*xi2 
     &                   +sum(QOp(i1,1:Nlon,k))/dble(Nlon)*xi3
            PLIH(i,j,k,3)=QHp(i1,i3,k)*xi1 +QHp(i1,i4,k)*xi2 
     &                   +sum(QHp(i1,1:Nlon,k))/dble(Nlon)*xi3
            PLIH(i,j,k,4)=OpL(i1,i3,k)*xi1 +OpL(i1,i4,k)*xi2 
     &                   +sum(OpL(i1,1:Nlon,k))/dble(Nlon)*xi3
            PLIH(i,j,k,6)=HpL(i1,i3,k)*xi1 +HpL(i1,i4,k)*xi2 
     &                   +sum(HpL(i1,1:Nlon,k))/dble(Nlon)*xi3
            PLIH(i,j,k,7)=Qep(i1,i3,k)*xi1 +Qep(i1,i4,k)*xi2 
     &                   +sum(Qep(i1,1:Nlon,k))/dble(Nlon)*xi3
	      PLIH(i,j,k,2)=0.0D0
            PLIH(i,j,k,5)=0.0D0
            IF(iHe.eq.1)THEN
              PLIH(i,j,k,2)=QHep(i1,i3,k)*xi1 +QHep(i1,i4,k)*xi2 
     &                     +sum(QHep(i1,1:Nlon,k))/dble(Nlon)*xi3
		    PLIH(i,j,k,5)=HepL(i1,i3,k)*xi1 +HepL(i1,i4,k)*xi2 
     &                     +sum(HepL(i1,1:Nlon,k))/dble(Nlon)*xi3
            ENDIF
            REIH(i,j,k,5)=cNOp(i1,i3,k)*xi1 +cNOp(i1,i4,k)*xi2  
     &                   +sum(cNOp(i1,1:Nlon,k))/dble(Nlon)*xi3
            REIH(i,j,k,4)=cO2p(i1,i3,k)*xi1 +cO2p(i1,i4,k)*xi2
     &                   +sum(cO2p(i1,1:Nlon,k))/dble(Nlon)*xi3
            REIH(i,j,k,6)=cN2p(i1,i3,k)*xi1 +cN2p(i1,i4,k)*xi2 
     &                   +sum(cN2p(i1,1:Nlon,k))/dble(Nlon)*xi3
	    ENDIF
            OpHr=2.5D-17*dsqrt(PNH(i,j,k,1))*PNH(i,j,k,6)           !O+ + H --> H+ +O
            HpOr=2.2D-17*dsqrt(REIH(i,j,k,9))*PNH(i,j,k,2)          !H+ + O --> O+ +H
	      per=4.43D-18*(REIH(i,j,k,8)/3.0D2)**0.7D0*REIH(i,j,k,7) !O+ (He+, H+) +e --> O (He, H)
            PLIH(i,j,k,4)=PLIH(i,j,k,4)  +OpHr    +per
            PLIH(i,j,k,6)=PLIH(i,j,k,6)  +HpOr    +per
            PLIH(i,j,k,1)=PLIH(i,j,k,1)  +HpOr*REIH(i,j,k,3)
            PLIH(i,j,k,3)=PLIH(i,j,k,3)  +OpHr*REIH(i,j,k,1)
            IF(iHe.eq.1)THEN
              HepO2r=8.0d-16*PNH(i,j,k,3)                           !He+ + O2 --> He +O + O+
              PLIH(i,j,k,5)=PLIH(i,j,k,5)  +HepO2r  +per
              PLIH(i,j,k,1)=PLIH(i,j,k,1)  +HepO2r*REIH(i,j,k,2)
            ENDIF
	    ENDDO
	  ENDDO
	ENDDO
      END SUBROUTINE GCM2TIME3D
c---------------------------------------------------------------------
      SUBROUTINE Dynamo2TIME3D(ed1,ed2)    !QHp,HepL,HpL,
      use maggrid      ,only: nmlat,nmlon,nmlonp1
      INCLUDE 'EIG.h'
	COMMON/EEmC/  EEm(NLi,NFi,2),EEH(NHN,2,2)
	COMMON/GRIDI/ alti(NLi),Cli(NLi),SIi(NLi),MLATi(lmi+1,NLi),Mloni(
     &              NFi),altH(NHN),ClH(NHN),SIH(NHN),MlatH(NHN,2,NLi),
     &              MlonH(NHN)
	COMMON/GridD/ Be3i(lmi,NFi),d11i(lmi,NFi,NLi),d12i(lmi,NFi,NLi),
     &              d22i(lmi,NFi,NLi),Dhi(lmi,NFi,NLi),Be3h(NHN,2),d11h(
     &              NHN,2,NLi),d12h(NHN,2,NLi),d22h(NHN,2,NLi),Dhh(NHN,2
     &              ,NLi),xije(NLi,NFi,4),xijeh(NHN,2,4),xie(nmlat,nmlon
     &              ,4),siy(nmlat+1),ije(NLi,NFi,4),ijeh(NHN,2,4),ie(
     &              nmlat,nmlon,7),idmax
      REAL*8 ed1(nmlat,nmlonp1),ed2(nmlat,nmlonp1),ed2x(nmlat,nmlonp1)

      ed2x=ed2
	DO i=1,nmlat
	  xy=siy(i)
        ed2x(i,1:nmlonp1)=ed2(i,1:nmlonp1)*xy
	ENDDO
	xy=siy(nmlat+1)

c mid- and low-latitude
	DO j=1,NFi
	  DO i=1,NLi
	    i1=ije(i,j,1)
	    i2=ije(i,j,2)
	    i3=ije(i,j,3)
	    i4=ije(i,j,4)
	    xi1=xije(i,j,1)
	    xi2=xije(i,j,2)
	    xi3=xije(i,j,3)
	    xi4=xije(i,j,4)
          EEm(i,j,1)=ed1(i1,i3)*xi1  +ed1(i1,i4)*xi2 
     &              +ed1(i2,i3)*xi3  +ed1(i2,i4)*xi4
	    EEm(i,j,1)=EEm(i,j,1)*xy
          EEm(i,j,2)=ed2x(i1,i3)*xi1 +ed2x(i1,i4)*xi2 
     &              +ed2x(i2,i3)*xi3 +ed2x(i2,i4)*xi4
	  ENDDO
	ENDDO
c high-latitude
	DO j=1,2
	  DO i=1,NHN
	    i1=ijeh(i,j,1)
	    i2=ijeh(i,j,2)
	    i3=ijeh(i,j,3)
	    i4=ijeh(i,j,4)
	    xi1=xijeh(i,j,1)
	    xi2=xijeh(i,j,2)
	    xi3=xijeh(i,j,3)
	    xi4=xijeh(i,j,4)
          EEH(i,j,1)=ed1(i1,i3)*xi1  +ed1(i1,i4)*xi2 
     &              +ed1(i2,i3)*xi3  +ed1(i2,i4)*xi4
	    EEH(i,j,1)=EEH(i,j,1)*xy
          EEH(i,j,2)=ed2x(i1,i3)*xi1 +ed2x(i1,i4)*xi2 
     &              +ed2x(i2,i3)*xi3 +ed2x(i2,i4)*xi4
	  ENDDO
	ENDDO
      END SUBROUTINE Dynamo2TIME3D
C----------------------------------------------------------------------
      REAL*8 FUNCTION EThermCon(cO,cN2,cO2,Ne,Te)
C Input: cO cN2 cO2 Ne:  Neutral O N2 O2 and Electron density of  in m-3 Te:  Electron 
C   temperature in K; Output: Electron Thermal Conductivity in J.m-1.s-1.k-1
C Refenrence: Schunk and Nagy, Rev. Geophys. Space Phys.,16, 355-399, 1978 
C    Millward G. H. & Baily G. J., Handbook of ionosphere model, 1996
      IMPLICIT REAL*8 (A-Z)
C Electron thermal conductivity
        dTe=dsqrt(Te)            
        qO=1.10D-16*(1.0D0+5.70D-4*Te)
	  qO2=2.20D-16*(1.0D0+3.60D-2*dTe)
        qN2=2.82D-17*dTe*(1.0D0-1.210D-4*Te)
	  ADD=qO*cO+qN2*cN2+qO2*cO2
        EThermCon=1.6D-19*7.7D7*Te*Te*dTe/(1.0D0+3.22D4*Te*Te*ADD/(Ne+
     &            1.0D-9))
      END FUNCTION EThermCon
C----------------------------------------------------------------------
      REAL*8 FUNCTION xIonThermConductivity(cOi,cHei,cHi,Ti)
C Input: cOi cHei cHi Ne:  O+ He+ H+ Electron density in m-3; Ti:  Ion temperature in K
C Output:  Ion thermal Conductivity in unit of J.m-1.s-1.k-1
C Refenrence: Rees and Roble, Rev. Geophys. Space Phys.,13(1), 201-242, 1975 
      IMPLICIT REAL*8 (A-Z)
        xIonThermConductivity=1.6D-19*1.15D6*Ti*Ti*DSQRT(Ti)*(cOi
     &                      +2.0D0*cHei+4.0D0*cHi)/(cOi+cHei+cHi)
	END FUNCTION xIonThermConductivity
C----------------------------------------------------------------------
      SUBROUTINE ElectronHeatingRate(cO,cN2,cO2,Ni,cH,Ne,Tn,Te,Ti,XT,YT)
C Input:  sdec: Solar zenith in degree   ALT:  Altitude in meters
C   cO cN2 cO2 cH Ne: Neutral O N2 H O2 Electron density of in m-3; Tn Te Ti: 
C   Neutral Electron Ion temperature in K;  F107: Solar 10.7 cm flux index
C Output:  XT YT: Coefficient XT YT in the equation HeatingRate=XT+YT*T
C Reference: Schunk and Nagy, Rev.Geophy. Space Phy.,16(3),355-399,1978
C   (Millward et.al, 239) and (Baily & Balan, 173,) handbook of ionospheric models, 1996
C NOTE:  Calculate heating rates(include positive and negative) for electron
      IMPLICIT  REAL*8 (A-Z)
	REAL*8 XN(3),DX(1:3),EX(1:3),Ni(6)
C Heating rate due to EUV radiation
C ---- This part is calculated out of this subroutine to improve the program speed Heating rate between electrons and neutrals
C  1---- Heating due to electron-neutral interactions
C   Reference: Schunk and Nagy, Rev.Geophy. Space Phy.,16(3),355-399,1978, EQ-43a - EQ-43e
       TDIF1=dsqrt(Te)
	 TDIF=Te-Tn
	 TDIFRT=TDIF/Tn/Te
	 LEN2=1.77D-25*Ne*cN2*(1.0D0-1.21D-4*Te)*Te
       LEO2=1.21D-24*Ne*cO2*(1.0D0+3.6D-2*TDIF1)*TDIF1
	 LEO=7.9D-25*Ne*cO*(1.0D0+5.7D-4*Te)*TDIF1
       LEH=9.63D-22*Ne*cH*(1.0D0-1.35D-4*Te)*dsqrt(Te)
C  2---- Heating due to rotational excitation of N2 and O2
C        Reference: Schunk and Nagy, Rev.Geophy. Space Phy.,16(3),355-399,1978, EQ-18 - EQ-19  
        LRN2=2.9D-20*Ne*cN2/TDIF1
	  LRO2=6.9D-20*Ne*cO2/TDIF1
C  3---- Heating due to vibrational excitation of N2 and O2
C     Reference: EQ-140 of Millward et.al, 239, in Handbook of ionospheric models
        LVN2b=6.5D-28*Ne*cN2*(Tn-3.1D2)*(Tn-3.1D2)*dexp(2.3D-3*TDIF)
        LVN2=(1.0D0+2.3D-3*TDIF)*LVN2b
        LVN2a=2.3D-3*TDIF*LVN2b*(Tn-Te)
C     Schunk and Nagy, Rev.Geophy. Space Phy.,16(3),355-399,1978, EQ-23,24
      HS=(3.3D3-8.39D2*dsin(1.91D-4*(Te-2.7D3)))*(Te-7.0D2)/(Te*7.0D2)
       LVO2=-5.196D-19*Ne*cO2*dexp(HS)*(dexp(-2.77D3*TDIFRT)-1.0D0)
C  4------ Heating due to fine structure excitation of O
c  from Dalgarno, A., Inelastic collisions at low energies, Can. J. Chem., 47, 1723-17291, 969.
        LFO=3.4D-18*Ne*cO*(1.0D0-7.D-5*Te)*(1.5D2/Te+0.4D0)/Tn
C  5-------- Heating due to electronic excitation of O
C            Reference: Schunk and Nagy, Rev.Geophy. Space Phy.,16(3),355-399,1978, EQ-34 - EQ-35
       DE=2.4D4+(Te-1.5D3)*(0.3D0-1.947D-5*(Te-4.0D3))
       EXH=3.3333D-4*DE*(Te-3.0D3)/Te
	 IF(EXH.GT.70.0D0) EXH=70.0D0
       EXH2=2.2713D4*TDIFRT
	 IF(EXH2.GT.70.0D0) EXH2=70.0D0
       LF1D=-1.57D-18*Ne*cO*dexp(EXH)*(dexp(-EXH2)-1.0D0)
       KEN=LFO+LEN2+LEO2+LEO+LEH+LRN2+LRO2+LVN2
       Qenb=-(LVO2+LF1D+LVN2a)+KEN*Tn
	 Qenk=-KEN
C Heating rate between electrons and ions
C   Reference: EQ-128 of Millward et.al, 239, in Handbook of ionospheric models, 1996
	 Qeik=7.7D-12*Ne*(Ni(1)/16.0d0+Ni(2)/4.0d0+Ni(3)+Ni(4)/32.0d0+Ni(5
     &      )/30.0d0+Ni(6)/28.0d0)/(Te*TDIF1)
       XT=Qenb+Qeik*Ti
	 YT=Qenk-Qeik
      RETURN
	END SUBROUTINE ElectronHeatingRate
C----------------------------------------------------------------------
      SUBROUTINE IonHeatingRates(Ti,Tn,Te,cO,cN2,cO2,cN,cH,cHe,dv2,Ne,
     &	                           Ni,XT,YT)
C Input: cO cN2 cO2 cN cH Ne: Neutral O N2 O2 N H Electron density of in m-3; Tn Te Ti: 
C    Ni: ion density in m-3; dv2: (V-U*B/Bs)^2 in (m/s)^2;
C Output: XT, YT:  Here HeatingRate: Qi + Fin => XT + YT * Ti.
C Reference: Millward et.al, handbook of ionospheric models, 239,1996
C    Baily & Balan, handbook of ionospheric models, 173,1996
C NOTE:  Heating rate of ions, in unit of eV m-3 s-1. Only consider O+
      IMPLICIT REAL*8 (A-Z)
      PARAMETER (ec=1.6D-19)        !unit charge
      REAL*8 Ni(6)
	 Atom=1.67D-27
C Heating rate due to ion-electron collisions
C   Reference: Equ.(128) of Millward et.al, 239,in Handbook of ionospheric models,1996
	 Qie=7.7D-12*Ne*(Ni(1)/16.0d0+Ni(2)/4.0d0+Ni(3)+Ni(4)/32.0d0+Ni(5)
     &     /30.0d0+Ni(6)/28.0d0)/(Te*dsqrt(Te))
C Heating rates due to ion-neutral elastic collisions
C   Reference: Table(5) of Baily & Balan, 173, in Handbook of ionospheric models,1996
       Tr=Ti+Tn
	 xTr=dsqrt(Tr)
       Qin=(2.1D-21*cO*xTr +6.60D-20*cN2 +5.80D-20*cO2 +3.30D-20*cH
     &+2.80D-20*cHe)*Ni(1) +(4.0D-21*cHe*xTr +5.70D-20*cO +5.3D-20*cN2 
     &+4.50D-20*cO2 +1.0D-19*cH)*Ni(2) +(1.4D-20*cH*xTr +3.50D-20*cO 
     &+3.1D-20*cN2 +2.80D-20*cO2 +5.5D-20*cHe)*Ni(3) +(1.4D-21*cO2*xTr 
     &+4.538D-20*cO +5.807D-20*cN2)*Ni(4) +(4.5D-20*cO+5.45D-20*cO2
     &+5.916D-20*cN2)*(Ni(5)+Ni(6))	 
C Heating rates due to ion-neutral inelastic collisions
       xTi=dsqrt(Ti)
	 Tr=Tr*0.5D0
	 iQin=3.8D-21*(8.0/9.0*cO*Ni(3)*Tn*xTi-cH*Ni(1)*Ti)
     &+3.4D-21*(9.0/8.0*cH*Ni(1)*Tn*dsqrt(Tn)-cO*Ni(3)*Ti*xTi)
C Frictional heating rate due to the relative motion between ions and neutral gas
       xtr=dsqrt(Tr)
	 ytr=dlog10(Tr)
       Cin=(3.67D-17*cO*xTr*(1.0-0.067*yTr)**2*8.0 +6.64D-16*cO2*10.667
     &+6.82D-16*cN2*10.18 +1.32D-16*cHe*3.2D0 +4.63D-18*cH*dsqrt(Tn+Ti
     &/16.0D0)*0.941 +4.62D-16*cN*7.467)*Ni(1) +(8.73D-17*cHe*xTr*(
     &1.0-0.093*yTr)**2*2.0D0  +10.1D-16*cO*3.2D0 +15.3D-16*cO2*3.556D0
     &+16.0D-16*cN2*3.5D0 +4.71D-16*cH*0.8D0)*Ni(2) +(2.65D-16*cH*xTr*(
     &1.0-0.083*yTr)**2*0.5D0 +6.61D-17*cO*xTr*(1.0-0.047*dlog10(Ti))**2
     &*0.9412 +32.0D-16*cO2*0.9697 +33.6D-16*cN2*0.966 +10.6D-16*cHe*0.8
     &)*Ni(3) +(2.59D-17*cO2*xTr*(1.0-0.073*yTr)**2*16.0 +2.31D-16*cO*
     &10.67 +4.13D-16*cN2*14.93 +0.7D-16*cHe*3.56 +0.7D-16*cH*0.97)*Ni(4
     &) +(2.44D-16*cO*10.43 +4.27D-16*cO2*15.48 +4.34D-16*cN2*14.48 
     &+0.74D-16*cHe*3.53 +0.69D-16*cH*0.9677)*Ni(5)
	 Fin=Atom*Cin*dv2
C HeatingRate: Qi + Fin => XT + YT * Ti.
       XT=Qie*Te +Qin*Tn +Fin/ec -iQin
	 YT=-Qie-Qin
	 END SUBROUTINE IonHeatingRates
C---------------------------------------------------------
      SUBROUTINE PreThdff() ! PreCollisions between ions coefficients
       IMPLICIT REAL*8 (A-H,O-Z)
	 COMMON/PreThdf/ A(3),Ast(3,3),cpre(3,3),Cst(3,3),cospp(3,3),
     &costpp(3,3)
       INTEGER s,t
       DATA A/16.0D0,4.0D0,1.0D0/

       DO s=1,3
         DO t=1,3
             Ast(s,t)=A(s)*A(t)/(A(s)+A(t))
             cpre(s,t)=1.27D-6*dsqrt(Ast(s,t))/A(s)
         ENDDO
       ENDDO
       DO s=1,3  ! Correction
         DO t=1,3
           Bst3=Ast(s,t)/A(s)*Ast(s,t)/A(s)
           Dst1=3.0D0*(Ast(s,t)/A(t))*(Ast(s,t)/A(t))-0.2D0*Bst3+
     &          Ast(s,t)/(A(s)+A(t))*0.1D0
           Dst4=3.0D0*(Ast(s,t)/A(s))*(Ast(s,t)/A(s))-0.2D0*Bst3-
     &          Ast(s,t)/(A(s)+A(t))*(1.6D0*A(t)/A(s)+1.5D0)
           Cst(s,t)=2.5D0-Ast(s,t)/A(t)*(2.5D0+A(t)/A(s))
           cospp(s,t)=1.25D0*(Dst1+1.5D0*Ast(s,t)/A(s))
           costpp(s,t)=1.25D0*(Dst4+1.5D0*Ast(s,t)/A(s))
          ENDDO
       ENDDO
      RETURN
      END SUBROUTINE PreThdff
C--------------------------------------------------------------------------
      SUBROUTINE Thdff(Ni,Ti,xt1,Cijx,Betap) !Collisions between ions and Thermal diffusion coef
      IMPLICIT REAL*8 (A-H,O-Z)
	COMMON/PreThdf/ A(3),Ast(3,3),cpre(3,3),Cst(3,3),cospp(3,3),
     &                 costpp(3,3)
      REAL*8 Cijx(3,3),Betap(3),Cij(3),eps(3,3),cost(3,3),Ni(3),Vstu(3,
     &       3),Wstu(3,3),cosp(3),costp(3,3),costx(3,3)
      INTEGER is,s,t,u,tx(2)
C Assume all ions with the same temperature Ti; Collisions frequency between ions; O+,He+,H+
        Tst=Ti
	  Cij(1:3)=0.0D0
	  Betap=0.0D0
	  Tst1=(Ti**1.5D0)  !1.0D0/
        DO s=1,3
	    cost(s,1:3)=cpre(s,1:3)*Ni(1:3)/Tst1  !
        ENDDO
        cosp(1)=cost(1,1)+cost(1,2)*cospp(1,2)+cost(1,3)*cospp(1,3)
        cosp(2)=cost(2,2)+cost(2,1)*cospp(2,1)+cost(2,3)*cospp(2,3)
        cosp(3)=cost(3,3)+cost(3,1)*cospp(3,1)+cost(3,2)*cospp(3,2)
        DO s=1,3  ! Correction
          DO t=1,3
            IF(s.EQ.t) CYCLE
              costp(s,t)=costpp(s,t)*cost(s,t)
	        costx(s,t)=costp(s,t)/cosp(s)
          ENDDO
        ENDDO
	  eps(1,2)=costx(1,2)*costx(2,1)
	  eps(2,1)=eps(1,2)
	  eps(2,3)=costx(2,3)*costx(3,2)
	  eps(3,2)=eps(2,3)
	  eps(3,1)=costx(3,1)*costx(1,3)
	  eps(1,3)=eps(3,1)
        DO s=1,3
          DO t=1,3
            IF(s.EQ.t) CYCLE
            u=6-s-t
	      Vstu(s,t)=costx(s,t)+costx(s,u)*costx(u,t) 
            Wstu(s,t)=1.0D0-eps(s,u)-Vstu(s,t)
          ENDDO
        ENDDO
	  xx1=(costp(1,2)*costp(2,3)*costp(3,1)+costp(2,1)*costp(3,2)*
     &      costp(1,3))/(cosp(1)*cosp(2)*cosp(3))
	  DO is=1,3
	    Cij(1:3)=cost(is,1:3)
          t=1
		IF(is.EQ.1)t=2
		u=6-is-t
          epsi=eps(is,t)+eps(t,u)+eps(u,is)+xx1
          faist=cost(is,t)*Ast(is,t)/Tst
          faisu=cost(is,u)*Ast(is,u)/Tst
	    x=15.0D0/(8.0D0*(1.0D0-epsi))*Ti
          beta=x/(cosp(is)*A(is))*(faist*Wstu(t,is)+faisu*Wstu(u,is))
          btast=x*A(is)/(A(t)*A(t)*cosp(t))
     &          *(faist*Wstu(is,t)-faisu*(Vstu(is,t)-Vstu(u,t)))
          btasu=x*A(is)/(A(u)*A(u)*cosp(u))
     &          *(faisu*Wstu(is,u)-faist*(Vstu(is,u)-Vstu(t,u)))
          dltst=0.4D0*(Cst(is,t)*beta+Ni(is)/Ni(t)*Cst(t,is)*btast)
          dltsu=0.4D0*(Cst(is,u)*beta+Ni(is)/Ni(u)*Cst(u,is)*btasu)
          dRstu=0.4D0/A(is)*(A(t)*cost(t,u)*Cst(t,u)*btast
     &          -A(u)*cost(u,t)*Cst(u,t)*btasu)
          Cijx(is,t)=Cij(t)*(1.0D0-dltst)-dRstu    !Modified Collision Frequency
          Cijx(is,u)=Cij(u)*(1.0D0-dltsu)+dRstu
		Cijx(is,is)=xt1    !0.0D0
          Betap(is)=beta-btast-btasu
	  ENDDO
      RETURN
      END SUBROUTINE Thdff
C-----------------TRIDAG---------------------------------------------
      SUBROUTINE TRIDAGI(DELTA,MS,NF,NL,A,B,C,D)
C Input: MS: The size of A,B,C,D,and DELTA; NF,NL: The index of the first and last data in A,B,C,D,and DELTA
C        A,B,C,D: Coefficients of LINEAR SIMULTANEOUS EQUATIONS;  OutPut: DELTA: SOLUTION VECTOR
C NOTE: Solving A SYSTEM OF LINEAR SIMULTANEOUS EQUATIONS WITH A TRIDIAGONAL COEFF MATRIX. THE EQNS ARE NUMBERED
C FROM NF TO NL, THEIR  SUB-DIAG. , DIAG. ,& SUPER-DIAG COEFFS. ARE STORED IN THE ARRAYS A, B, C. THE RIGHT HAND
C SIDE OF THE VECTOR IS STORED IN D. THE COMPUTED SOLUTION VECTOR IS STORED IN THE ARRAY DELTA. THIS ROUTINE COMES
C FROM CARNAHAN, LUTHER, AND WILKES, APPLIED NUMERICAL METHODS, WILEY, 1969,PAGE 446
      IMPLICIT REAL*8 (A-H,O-Z)
      INTEGER MS,NF,NL,NFP1,I,K,LAST
      REAL*8 A(MS),B(MS),C(MS),D(MS),DELTA(MS),alfa(MS),GAMMA(MS)
        alfa(NF)=B(NF)        !COMPUTE INTERMEDIATE ARRAYS alfa & GAMMA
        GAMMA(NF)=D(NF)/alfa(NF)
	  NFP1=NF+1
        DO I=NFP1,NL
          alfa(I)=B(I)-A(I)*C(I-1)/alfa(I-1)
          GAMMA(I)=(D(I)-A(I)*GAMMA(I-1))/alfa(I)
        ENDDO
        DELTA(NL)=GAMMA(NL)
	  LAST=NL-NF !COMPUTE FINAL SOLUTION VECTOR V
        DO K=1,LAST
          I=NL-K
		DELTA(I)=GAMMA(I)-C(I)*DELTA(I+1)/alfa(I)
        ENDDO
      END SUBROUTINE TRIDAGI
C--------------------------------------------------------------------
	SUBROUTINE CollisionFrequency(Tix,Tn,cO,cO2,cN2,cN,cH,cHe,Cin)
C Input: Ti,Tn: Ionic (O+,H+,He+), Neutral temperature in K; cO,cO2,cN2
C        ,cN,cH,cHe: Density of neutral O,O2,N2,N,H,He in m-3
C Output: Cin:  Collision frequency in s-1 for  O+,He+,H+,O2+,NO+
C Reference: Schunk, PAGEOPH, 127(2/3), 255-303,1988
      IMPLICIT REAL*8 (A-Z)
	real*8 Ti(3),Cin(5)
        Ti(1:3)=Tix
	  Tr=(Ti(1)+Tn)*0.5D0  !O+ - n
C     O+-O  Two choices ,one from Salah [1993], another from Pesnel
        Cin(1)=5.9D-17*cO*dsqrt(Tr)*(1.0D0-0.096d0*dlog10(Tr))**2
     &    +6.66D-16*cO2 +6.82D-16*cN2 +4.62D-16*cN
     &    +4.63D-18*cH*dsqrt(Tn+Ti(1)/16.0D0)+1.32D-16*cHe
        Cin(3)=6.61D-17*cO*dsqrt(Ti(2))*(1.0D0-4.7D-2*dlog10(Ti(2)))**2
     &    +3.36D-15*cN2 +3.20D-15*cO2 +1.06D-15*cHe +2.61D-15*cN+
     &    2.65D-16*cH*dsqrt(Tr)*(1.D0-8.2D-2*dlog10(Tr))**2
        Cin(2)=1.01D-15*cO +1.60D-15*cN2 +1.53D-15*cO2 +8.73D-17*cHe*
     & dsqrt(Tr)*(1.0D0-9.3D-2*dlog10(Tr))**2 +1.19D-15*cN +4.71D-16*cH
      END SUBROUTINE CollisionFrequency
c---------------------------------------------------------------
      subroutine photoelectronheatUp(NPt,NP,in300n,in300s,Ne,co2,cn2,co,
     &                               prodelec,length,yta,heatout)
      implicit real*8(a-z)
	integer in300n,in300s,i,j,NP,NPt
	dimension Ne(NPt),co2(NPt),cn2(NPt),co(NPt),prodelec(NPt),length(
     &          NPt),yta(NPt),heatout(NPt),xt(NPt),EEC(7),EEX(7)
	Data EEC /5.342D0,1.056D0,-4.392D-2,-5.900D-2,-9.346D-3,-5.755D-4,
     &         -1.249D-5/

	heatout=0.0
	EEX(1)=1.0D0
      if (in300n.LT.0.or.in300n.ge.in300s) then
	  do i=1,NP
          EEX(2)=Ne(i)/(co2(i)+cn2(i)+co(i))
	    if(EEX(2).lt.3.0D-8)EEX(2)=3.0D-8
		if(EEX(2).gt.10.0)EEX(2)=10.0
          EEX(2)=dlog(EEX(2))
		EEX(3)=EEX(2)*EEX(2)
		EEX(4)=EEX(3)*EEX(2)
          EEX(5)=EEX(4)*EEX(2)
		EEX(6)=EEX(5)*EEX(2)
		EEX(7)=EEX(6)*EEX(2)
          xtemp2=dexp(sum(EEX(1:7)*EEC(1:7)))
	    heatout(i)=xtemp2*prodelec(i)
	  enddo
	else
	  do i=1,in300n
          EEX(2)=Ne(i)/(co2(i)+cn2(i)+co(i))
	    if(EEX(2).lt.3.0D-8)EEX(2)=3.0D-8
		if(EEX(2).gt.10.0)EEX(2)=10.0
          EEX(2)=dlog(EEX(2))
		EEX(3)=EEX(2)*EEX(2)
		EEX(4)=EEX(3)*EEX(2)
          EEX(5)=EEX(4)*EEX(2)
		EEX(6)=EEX(5)*EEX(2)
		EEX(7)=EEX(6)*EEX(2)
          xtemp2=dexp(sum(EEX(1:7)*EEC(1:7)))
	    heatout(i)=xtemp2*prodelec(i)
	  enddo
	  do i=in300s,NP
          EEX(2)=Ne(i)/(co2(i)+cn2(i)+co(i))
	    if(EEX(2).lt.3.0D-8)EEX(2)=3.0D-8
		if(EEX(2).gt.10.0)EEX(2)=10.0
          EEX(2)=dlog(EEX(2))
		EEX(3)=EEX(2)*EEX(2)
		EEX(4)=EEX(3)*EEX(2)
          EEX(5)=EEX(4)*EEX(2)
		EEX(6)=EEX(5)*EEX(2)
		EEX(7)=EEX(6)*EEX(2)
          xtemp2=dexp(sum(EEX(1:7)*EEC(1:7)))
	    heatout(i)=xtemp2*prodelec(i)
	  enddo
	  xt(in300n)=0.0D0
	  q0s= heatout(in300s)/ne(in300s)
	  cqe= 3.0D-18
	  q0n= heatout(in300n)/ne(in300n)
        do i=in300n+1,in300s
          xt(i)=xt(i-1)+0.5*(ne(i)+ne(i-1))*length(i)
        enddo
	  do i=in300n+1,in300s-1
          xn=xt(i)
		xintn=dmin1(cqe*xn,30.0D0)
		xs=xt(in300s)-xn
		xints=dmin1(cqe*xs,30.0D0)
          xqs=ne(i)*q0s*yta(i)/yta(in300s)*dexp(-xints)
          xqn=ne(i)*q0n*yta(i)/yta(in300n)*dexp(-xintn)
          heatout(i)=xqs+xqn
        enddo
      endif
	return
      end subroutine photoelectronheatUp
c---------------------------------------------------------------------
	SUBROUTINE TIME3D2GCM(Z,Op,Te,Ti,Vi,Ped,Hall)
      use geogrid_t3d  ,only: nlev,nlat,nlon
      use InterpTI     ,only: XMLAI,XMLOI,Cpole,xij,xijh,xix,ij,ijh,ix
      INCLUDE 'EIG.h'
      PARAMETER (n0=8, n1=10, n2=n1-n0+2, n3=n2+1, n4=n2+2)
	COMMON/REIM3/ REIM(lmi,NFi,NLi,12),REIMold(lmi,NFi,NLi,8),REIH(NHN
     &             ,2,NLi,12),REIHold(NHN,2,NLi,8)
	COMMON/GRIDI/ alti(NLi),Cli(NLi),SIi(NLi),MLATi(lmi+1,NLi),Mloni(
     &              NFi),altH(NHN),ClH(NHN),SIH(NHN),MlatH(NHN,2,NLi),
     &              MlonH(NHN)
	COMMON/Iinterp/ Sl(NLi+nk+1,2),dmlon(NLi+nk),Xli(NLi+nk+1),MLATT(
     &                NLi+nk+1,NLi),NHni(NHN,2),Nmlo(Nli+nk),NHX(Nli+1:
     &                Nli+nk,NFi)
	COMMON/EConductance/ EconT(lmi,NFi,NLi,2),EconTH(NHN,2,NLi,2),imax
      real*8 Z(nlat,nlon,Nlev),Op(nlat,nlon,Nlev),Te(nlat,nlon,Nlev),Ti(
     &       nlat,nlon,Nlev),Vi(nlat,nlon,Nlev),Ped(nlat,nlon,Nlev),Hall
     &       (nlat,nlon,Nlev),REX(NLi,12),REP(12),Zg(Nlev+1),Xg(Nlev+1,
     &       12)

      DO i=1,Nlat
	  DO j=1,Nlon
	    DO k=1,NLi
	      i1=ix(i,j,k,1)
	      i2=ix(i,j,k,2)
	      i3=ix(i,j,k,3)
	      i4=ix(i,j,k,4)
	      i5=ix(i,j,k,5)
	      i6=ix(i,j,k,6)
	      i7=ix(i,j,k,7)
	      xi1=xix(i,j,k,1)
	      xi2=xix(i,j,k,2)
	      xi3=xix(i,j,k,3)
	      xi4=xix(i,j,k,4)
            IF(i7.EQ.1)THEN
              REX(k,1)=REIM(i1,i3,k,1)*xi1 +REIM(i1,i4,k,1)*xi2
     &                +REIM(i2,i3,k,1)*xi3 +REIM(i2,i4,k,1)*xi4
              REX(k,2:n2)=REIM(i1,i3,k,n0:n1)*xi1 +REIM(i1,i4,k,n0:n1)
     &           *xi2 +REIM(i2,i3,k,n0:n1)*xi3 +REIM(i2,i4,k,n0:n1)*xi4
            ELSE IF(i7.EQ.2)THEN
              REX(k,1)=REIM(i1,i2,k,1)*xi1 +REIM(i1,i3,k,1)*xi2
     &                +REIH(i4,i6,k,1)*xi3 +REIH(i5,i6,k,1)*xi4
              REX(k,2:n2)=REIM(i1,i2,k,n0:n1)*xi1 +REIM(i1,i3,k,n0:n1)
     &           *xi2 +REIH(i4,i6,k,n0:n1)*xi3 +REIH(i5,i6,k,n0:n1)*xi4
            ELSE IF(i7.EQ.3)THEN
              REX(k,1)=REIH(i1,i6,k,1)*xi1 +REIH(i2,i6,k,1)*xi2
     &                +REIH(i3,i6,k,1)*xi3 +REIH(i4,i6,k,1)*xi4
              REX(k,2:n2)=REIH(i1,i6,k,n0:n1)*xi1 +REIH(i2,i6,k,n0:n1)
     &           *xi2 +REIH(i3,i6,k,n0:n1)*xi3 +REIH(i4,i6,k,n0:n1)*xi4	
            ELSE IF(i7.EQ.4)THEN
c	        REP(1)=(REIH(1,i6,k,1)+REIH(2,i6,k,1) +REIH(3,i6,k,1))/3.0D0
c	        REP(2:n2)=(REIH(1,i6,k,n0:n1)+REIH(2,i6,k,n0:n1)+REIH(3,i6,k,n0:n1))/3.0D0  
              REP(1:n2)=0.0D0
              DO ik=1,Nmlo(NLi+nk)
	          REP(1)=REP(1)+REIH(ik,i6,k,1)
	          REP(2:n2)=REP(2:n2)+REIH(ik,i6,k,n0:n1)
              ENDDO
	        REP(1:n2)=REP(1:n2)/dble(Nmlo(NLi+nk))
             REX(k,1)=REIH(i1,i6,k,1)*xi1+REIH(i2,i6,k,1)*xi2+REP(1)*xi3
              REX(k,2:n2)=REIH(i1,i6,k,n0:n1)*xi1 +REIH(i2,i6,k,n0:n1)
     &                 *xi2 +REP(2:n2)*xi3
	      ENDIF	 
	    ENDDO
	    DO k=1,imax
            IF(i7.EQ.1)THEN
              REX(k,n3:n4)=EconT(i1,i3,k,1:2)*xi1+EconT(i1,i4,k,1:2)*xi2
     &                   +EconT(i2,i3,k,1:2)*xi3 +EconT(i2,i4,k,1:2)*xi4
            ELSE IF(i7.EQ.2)THEN
              REX(k,n3:n4)=EconT(i1,i2,k,1:2)*xi1+EconT(i1,i3,k,1:2)*xi2
     &                 +EconTH(i4,i6,k,1:2)*xi3 +EconTH(i5,i6,k,1:2)*xi4
            ELSE IF(i7.EQ.3)THEN
              REX(k,n3:n4)=EconTH(i1,i6,k,1:2)*xi1 +EconTH(i2,i6,k,1:2)
     &           *xi2 +EconTH(i3,i6,k,1:2)*xi3 +EconTH(i4,i6,k,1:2)*xi4	
            ELSE IF(i7.EQ.4)THEN
c	        REP(n3:n4)=(EconTH(1,i6,k,1:2)+EconTH(2,i6,k,1:2)+EconTH(3,i6,k,1:2))/3.0D0
              REP(n3:n4)=0.0D0
              DO ik=1,Nmlo(NLi+nk)
	          REP(n3:n4)=REP(n3:n4)+EconTH(ik,i6,k,1:2)
              ENDDO
	        REP(n3:n4)=REP(n3:n4)/dble(Nmlo(NLi+nk))
              REX(k,n3:n4)=EconTH(i1,i6,k,1:2)*xi1 +EconTH(i2,i6,k,1:2)
     &                    *xi2  +REP(n3:n4)*xi3
	      ENDIF	 
	    ENDDO
          xx=REX(imax,n3)/REX(imax,1)
          REX(imax+1:NLi,n3)=REX(imax+1:NLi,1)*xx
          xx=REX(imax,n4)/REX(imax,1)
          REX(imax+1:NLi,n4)=REX(imax+1:NLi,1)*xx
	    Zg(1:Nlev)=Z(i,j,1:Nlev)
          REX(1:NLi,1)=Dlog(REX(1:NLi,1))
          REX(1:NLi,n3:n4)=Dlog(REX(1:NLi,n3:n4))
          iy=0
          CALL DIGMAX(alti,REX,NLi,12,NLi,Zg,Xg,Nlev+1,12,Nlev,n4,iy,0)
          Op(i,j,1:Nlev)=Dexp(Xg(1:Nlev,1))
          Te(i,j,1:Nlev)=Xg(1:Nlev,2)
          Ti(i,j,1:Nlev)=Xg(1:Nlev,3)
          Vi(i,j,1:Nlev)=Xg(1:Nlev,4)
          Ped(i,j,1:Nlev)=Dexp(Xg(1:Nlev,5))
          Hall(i,j,1:Nlev)=Dexp(Xg(1:Nlev,6))
	  ENDDO	
	ENDDO
	END SUBROUTINE TIME3D2GCM
c---------------------------------------------------------------------
	SUBROUTINE TIME3D2Dynamo(zigm11,zigm22,zigm2,zigmc,rim1,rim2)
      use maggrid      ,only: nmlat,nmlon,nmlonp1
      INCLUDE 'EIG.h'
	COMMON/GridD/ Be3i(lmi,NFi),d11i(lmi,NFi,NLi),d12i(lmi,NFi,NLi),
     &              d22i(lmi,NFi,NLi),Dhi(lmi,NFi,NLi),Be3h(NHN,2),d11h(
     &              NHN,2,NLi),d12h(NHN,2,NLi),d22h(NHN,2,NLi),Dhh(NHN,2
     &              ,NLi),xije(NLi,NFi,4),xijeh(NHN,2,4),xie(nmlat,nmlon
     &              ,4),siy(nmlat+1),ije(NLi,NFi,4),ijeh(NHN,2,4),ie(
     &              nmlat,nmlon,7),idmax
	COMMON/GRIDI/ alti(NLi),Cli(NLi),SIi(NLi),MLATi(lmi+1,NLi),Mloni(
     &              NFi),altH(NHN),ClH(NHN),SIH(NHN),MlatH(NHN,2,NLi),
     &              MlonH(NHN)
	COMMON/IGrid/ d1i(lmi,NFi,NLi,3),d2i(lmi,NFi,NLi,3),glati(lmi,NFi,
     &              NLi),gloni(lmi,NFi,NLi),Vsi(lmi,NFi,NLi),Wsi(lmi,NFi
     &              ,NLi),Bsi(lmi,NFi,NLi,4),gxi(lmi,NFi,NLi,3),dvv(lmi,
     &              NFi,NLi,2),mlov(NLi,NFi,2),d1ih(NHN,2,NLi,3),d2ih(
     &              NHN,2,NLi,3),glatih(NHN,2,NLi),glonih(NHN,2,NLi),
     &              Vsih(NHN,2,NLi),Wsih(NHN,2,NLi),Bsih(NHN,2,NLi,4),
     &              gxih(NHN,2,NLi,3),dvvh(NHN,2,NLi,2),mlovh(NHN,2,2),
     &              MH300(NHN+NLi)
	COMMON/Iinterp/ Sl(NLi+nk+1,2),dmlon(NLi+nk),Xli(NLi+nk+1),MLATT(
     &                NLi+nk+1,NLi),NHni(NHN,2),Nmlo(Nli+nk),NHX(Nli+1:
     &                Nli+nk,NFi)
	COMMON/REIM3/ REIM(lmi,NFi,NLi,12),REIMold(lmi,NFi,NLi,8),REIH(NHN
     &             ,2,NLi,12),REIHold(NHN,2,NLi,8)
	COMMON/DynamoWind/ Ud(lmi,NFi,NLi,2),Udh(NHN,2,NLi,2)
	COMMON/REIM3N/ PN(lmi,NFi,NLi,8),PLI(lmi,NFi,NLi,8),PNH(NHN,2,NLi,
     &               8),PLIH(NHN,2,NLi,8)
	COMMON/EConductance/ EconT(lmi,NFi,NLi,2),EconTH(NHN,2,NLi,2),imax
      real*8 REP(6),DyP(6),DyPi(lmi,NFi,6),DyPH(NHN,2,6),IonD(4),Econd(
     &       3),zigm11(nmlat,nmlonp1),zigm22(nmlat,nmlonp1),zigm2(nmlat,
     &       nmlonp1),zigmc(nmlat,nmlonp1),rim1(nmlat,nmlonp1),rim2(
     &       nmlat,nmlonp1)

	DO j=1,NFi
	  DO i=1,lmi
	    k2=imin0(i,lmi+1-i,idmax)
	    DO k=1,k2
	      IonD(1)=REIM(i,j,k,1)
	      IonD(2:3)=REIM(i,j,k,4:5)
		  IonD(4)=sum(IonD(1:3))
            CALL ConductRate(IonD(4),REIM(i,j,k,8),REIM(i,j,k,9),PN(i,j,
     &                       k,1),PN(i,j,k,2),PN(i,j,k,3),PN(i,j,k,4),
     &                       IonD,Bsi(i,j,k,4),Econd)
            EconT(i,j,k,1:2)=Econd(2:3)
	    ENDDO
	  ENDDO
	ENDDO
	DO j=1,2
	  DO i=1,NHN
	    DO k=1,idmax
	      IonD(1)=REIH(i,j,k,1)
	      IonD(2:3)=REIH(i,j,k,4:5)
		  IonD(4)=sum(IonD(1:3))
            CALL ConductRate(IonD(4),REIH(i,j,k,8),REIH(i,j,k,9),PNH(i,j
     &                      ,k,1),PNH(i,j,k,2),PNH(i,j,k,3),PNH(i,j,k,4)
     &                      ,IonD,BsiH(i,j,k,4),Econd)
            EconTH(i,j,k,1:2)=Econd(2:3)
	    ENDDO
	  ENDDO
	ENDDO
      
	DO i=2,lmi-1
        k2=imin0(i,lmi+1-i,idmax)
	  DO j=1,NFi
	    Dy1=0.0D0
	    Dy2=0.0D0
	    Dy3=0.0D0
	    Dy4=0.0D0
	    Dy5=0.0D0
	    Dy6=0.0D0
	    DO k=1,k2
	      Ec1=EconT(i,j,k,1)
	      Ec2=EconT(i,j,k,2)
	      d11ix=d11i(i,j,k)
	      d22ix=d22i(i,j,k)
	      d12ix=d12i(i,j,k)
	      Dhix=Dhi(i,j,k)
	      U2x=Ud(i,j,k,2)
	      U1x=Ud(i,j,k,1)
	      x1=Ec1*d11ix*Dhix
	      x2=Ec1*d22ix*Dhix
	      x3=Ec2*Dhix
	      x4=Ec1*d12ix*Dhix
c	      Dy5=Dy5+ (Ec1*d11ix*U2x +(Ec2-Ec1*d12ix)*U1x)*Dhix 
c	      Dy6=Dy6+ ((Ec2-Ec1*d12ix)*U2x -Ec1*d22ix*U1x)*Dhix
	      Dy1=Dy1+ x1
	      Dy2=Dy2+ x2
	      Dy3=Dy3+ x3
	      Dy4=Dy4+ x4
	      Dy5=Dy5+ x1*U2x +(x3-x4)*U1x
	      Dy6=Dy6+ (x3+x4)*U2x -x2*U1x
	    ENDDO
	    Be3=Be3i(i,j)
	    DyPi(i,j,1)=Dy1
	    DyPi(i,j,2)=Dy2
	    DyPi(i,j,3)=Dy3
	    DyPi(i,j,4)=Dy4
          DyPi(i,j,5)=Dy5*Be3
          DyPi(i,j,6)=Dy6*Be3
	  ENDDO
	ENDDO

	DO j=1,2
	  DO i=1,NHN
	    Dy1=0.0D0
	    Dy2=0.0D0
	    Dy3=0.0D0
	    Dy4=0.0D0
	    Dy5=0.0D0
	    Dy6=0.0D0
	    DO k=1,idmax
	      Ec1=EconTH(i,j,k,1)
	      Ec2=EconTH(i,j,k,2)
	      d11ix=d11h(i,j,k)
	      d22ix=d22h(i,j,k)
	      d12ix=d12h(i,j,k)
	      Dhix=Dhh(i,j,k)
	      U2x=Udh(i,j,k,2)
	      U1x=Udh(i,j,k,1)
	      x1=Ec1*d11ix*Dhix
	      x2=Ec1*d22ix*Dhix
	      x3=Ec2*Dhix
	      x4=Ec1*d12ix*Dhix
	      Dy1=Dy1+ x1
	      Dy2=Dy2+ x2
	      Dy3=Dy3+ x3
	      Dy4=Dy4+ x4
	      Dy5=Dy5+ x1*U2x +(x3-x4)*U1x
	      Dy6=Dy6+ (x3+x4)*U2x -x2*U1x
	    ENDDO
	    Be3=Be3h(i,j)
	    DyPh(i,j,1)=Dy1
	    DyPh(i,j,2)=Dy2
	    DyPh(i,j,3)=Dy3
	    DyPh(i,j,4)=Dy4
          DyPh(i,j,5)=Dy5*Be3
          DyPh(i,j,6)=Dy6*Be3
	  ENDDO
	ENDDO

      DO i=1,nmlat
	  DO j=1,nmlon
	    i1=ie(i,j,1)
	    i2=ie(i,j,2)
	    i3=ie(i,j,3)
	    i4=ie(i,j,4)
	    i5=ie(i,j,5)
	    i6=ie(i,j,6)
	    i7=ie(i,j,7)
	    xi1=xie(i,j,1)
	    xi2=xie(i,j,2)
	    xi3=xie(i,j,3)
	    xi4=xie(i,j,4)
	    IF(i7.EQ.1)THEN
            DyP(1:6)=DyPi(i1,i3,1:6)*xi1 +DyPi(i1,i4,1:6)*xi2
     &              +DyPi(i2,i3,1:6)*xi3 +DyPi(i2,i4,1:6)*xi4
          ELSE IF(i7.EQ.2)THEN
            DyP(1:6)=DyPi(i1,i2,1:6)*xi1 +DyPi(i1,i3,1:6)*xi2
     &              +DyPH(i4,i6,1:6)*xi3 +DyPH(i5,i6,1:6)*xi4
          ELSE IF(i7.EQ.3)THEN
            DyP(1:6)=DyPH(i1,i6,1:6)*xi1 +DyPH(i2,i6,1:6)*xi2
     &              +DyPH(i3,i6,1:6)*xi3 +DyPH(i4,i6,1:6)*xi4	
          ELSE IF(i7.EQ.4)THEN
c           REP(1:6)=(DyPH(1,i6,1:6)+DyPH(2,i6,1:6)+DyPH(3,i6,1:6))/3.0D0
            REP(1:6)=0.0D0
            DO ij=1,Nmlo(NLi+nk)
	        REP(1:6)=REP(1:6)+DyPH(j,i6,1:6)
            ENDDO
	      REP(1:6)=REP(1:6)/dble(Nmlo(NLi+nk))
            DyP(1:6)=DyPH(i1,i6,1:6)*xi1 +DyPH(i2,i6,1:6)*xi2
     &               +REP(1:6)*xi3
	    ENDIF
          zigm11(i,j)=DyP(1)
          zigm22(i,j)=DyP(2)
          zigm2(i,j) =DyP(3)
          zigmc(i,j) =DyP(4)
          rim1(i,j)  =DyP(5)
          rim2(i,j)  =DyP(6)
	  ENDDO
        zigm11(i,nmlonp1)=zigm11(i,1)
        zigm22(i,nmlonp1)=zigm22(i,1)
        zigm2(i,nmlonp1) =zigm2(i,1)
        zigmc(i,nmlonp1) =zigmc(i,1)
        rim1(i,nmlonp1)  =rim1(i,1)
        rim2(i,nmlonp1)  =rim2(i,1)
	ENDDO
!
! Change sign of rim2 (whole-array operation):
        rim2 = -rim2
	END SUBROUTINE TIME3D2Dynamo
C--------------------------------------------------------------------
      SUBROUTINE ConductRate(Ne,Te,Ti,Tn,On,O2n,N2n,IonD,BMag,Econd)
      IMPLICIT REAL*8 (A-H,M,O-Z)
      DIMENSION Omega(4),Mass(4),Cin(4),IonD(4),ECond(3)
      REAL*8 IonD,Ne,Mass,N2n
	  ChargeUnit=1.602D-19
	  AtomUnit=1.66053D-27
	  Mass(1)=16.0D0*AtomUnit     !O+
	  Mass(2)=32.0D0*AtomUnit     !O2+
	  Mass(3)=30.0D0*AtomUnit     !NO+
	  Mass(4)=5.46448D-4*AtomUnit !e- 
	  Econd(1:3)=0.0D0
	  TTe=Tn
	  T300=Te/300.0D0  !T500=Tn/500.0D0
	  Tr=(Tn+Ti)*0.5D0
	  STr=dsqrt(Tr)
	  DlTr=dlog10(Tr)
	  XTr1=1.0D0-0.064*DlTr
	  XTr2=1.0D0-7.3D-2*DlTr
c       Taken from Schunk and Nagy, Rev. Geophys. Space Phys.,18,813,1980;
        Cin(1)= 6.82D-16*N2n +3.67D-17*STr*XTr1*XTr1*On +6.64D-16*O2n !O+ <-> n
        Cin(2)= 2.31D-16*On +2.59D-17*STr*XTr2*XTr2*O2n +4.13D-16*N2n !O2+ <-> n
        Cin(3)= 2.44D-16*On +4.34D-16*N2n +4.27D-16*O2n               !NO+ <-> n
	  Cin(4)=(4.11D0*N2n*(T300**0.95D0) +2.95D0*O2n*(T300**0.79D0)
     &         +1.09D0*On*(T300**0.85D0))*1.0D-26       !e- +n==>Cin(4)
	  T6=TTe*TTe
	  T6=T6*T6*T6
	  TTE15=1.0D0/(TTe*dsqrt(TTe))  !TTe**(-1.5D0)
	  Cie=Ne*(5.9D1+4.18D0*log10(T6/Ne))*TTE15*1.0D-6           !e+i>Cie
	  DO is=1,4
		Omega(is)=ChargeUnit*BMag/Mass(is)
          IF(is.eq.4)Cin(is)=Cin(is)*omega(is)/BMag+Cie
	  ENDDO
	  DO is=1,4
	    OC2=omega(is)*omega(is)+Cin(is)*Cin(is)
		al=IonD(is)*ChargeUnit/BMag*Cin(is)*omega(is)/OC2
		bl=IonD(is)*ChargeUnit/BMag*omega(is)*omega(is)/OC2
		IF(is.EQ.4)THEN
			Econd(2)=Econd(2)+al
			Econd(3)=Econd(3)+bl
		ELSE 
			Econd(2)=Econd(2)+al
			Econd(3)=Econd(3)-bl
		ENDIF
	  ENDDO
 	  IF(Econd(2).LT.1.0D-30)Econd(2)=1.0D-30
 	  IF(Econd(3).LT.1.0D-30)Econd(3)=1.0D-30
	END SUBROUTINE ConductRate
