	SUBROUTINE DynamoSim(ed1,ed2,HB)
      use maggrid      ,only: nmlat,nmlon,nmlonp1
      INCLUDE 'EIG.h'
      COMMON/input/ F107,F107A,AP,Uth,TimeStep,Hindex(6),iyd,IYR,IMO,IDA
      REAL*8 ed1(nmlat,nmlonp1),ed2(nmlat,nmlonp1)

	  CALL MagPOL(IYR,iyd,dCLATP,dPOLON)
        CALL Dsubsol(IYR,iyd,Uth,dCLATP,dPOLON,GLATo,dSMLON)
        call EEmodel(dSMLON,HB)
	  CALL EfieldTrans(ed1,ed2,dSMLON,HB)
      END SUBROUTINE DynamoSim
C-------------------------------------------------------------------------
	SUBROUTINE EEmodel(dSMLON,HB)
      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)
	COMMON/EEmC/  EEm(NLi,NFi,2),EEH(NHN,2,2),dE(NLi,NFi,2)
      COMMON/input/ F107,F107A,AP,Uth,TimeStep,Hindex(6),iyd,IYR,IMO,IDA
	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 Vup(2),VWE(2),XMLA(NHN),XMLT(NHN),ET1(NHN),EP1(NHN),EPO1(
     &       NHN),ET2(NHN),EP2(NHN),EPO2(NHN),ET3(NFi),EP3(NFi),EPO3(NFi
     &       ),XMLAL(NFi),XMLTL(NFi)

        RX2=(RE+300.0D0)/(RE+HB)
	  RX1=RX2*Dsqrt(RX2)
        IHR=INT(Uth)
	  IMN=INT((Uth-IHR)*60)                  !High latitude Electric field
        XMLA(1:NHN)=MH300(1:NHN)               !North Hemisphere
        DO ij=1,NHN
          XMLT(ij)=(MlonH(ij)-dSMLON)/15.0D0+12.0D0
	    IF (XMLT(ij) .GE. 24.0D0) XMLT(ij) = XMLT(ij) - 24.0D0
	    IF (XMLT(ij) .LT.  0.0D0) XMLT(ij) = XMLT(ij) + 24.0D0
	  ENDDO
	  CALL Dweimer96(IYR,IMO,IDA,IHR,IMN,XMLA,XMLT,Hindex(1),Hindex(2),
     &               Hindex(3),NHN,EPT1,ET1,EP1,EPO1)
        XMLA(1:NHN)=0.0D0-XMLA(1:NHN)                  !South Hemisphere
	  CALL Dweimer96(IYR,IMO,IDA,IHR,IMN,XMLA,XMLT,Hindex(1),Hindex(2),
     &               Hindex(3),NHN,EPT2,ET2,EP2,EPO2)
	  XX=dmax1(EPT1,EPT2)
	  XX=Hindex(6)/XX
	  IF(XX.LE.0.0D0)XX=1.0D0
C     ET,EP= Etheta,Ephi(magnetic equatorward,eastward) E field component (V/m)
	  EEH(1:NHN,1,1)=EP1(1:NHN)*RX1*XX
	  EEH(1:NHN,2,1)=EP2(1:NHN)*RX1*XX
	  EEH(1:NHN,1,2)=ET1(1:NHN)*RX2/SIH(1:NHN)*XX
        EEH(1:NHN,2,2)=ET2(1:NHN)*RX2/SIH(1:NHN)*XX

        XMLTL(1:NFi)=XMLT(NHN-NFi+1:NHN) !Mid- and Low-latitude Electric field
	  DO k=1,NLi
	    IF(k.GT.9)Altk=alti(9)
		IF(k.LE.9)Altk=alti(k)
		Ml=MH300(NHN+k)
          DO j=1,NFi
            SLT=UTh+gloni(k,j,k)/15.0D0
		  if(SLT.GT.24.0D0)SLT=SLT-24.0D0
	      call dEdrift(MLONi(j),MLATi(lmi+1,k),iyd,SLT,Altk,Altk+1.0D0
     &                  ,Vup,VWE)      !Vup,VWE: upward, eastward drift in m/s
		  EEm(k,j,1)=Vup(1)*dE(k,j,1)
		  EEm(k,j,2)=VWE(1)*dE(k,j,2)
          ENDDO
          IF(Ml.GE.50.0D0)THEN !Electric field from Middle to high latitude
            XMLAL(1:NFi)=Ml
		  xy=(6.0D1-Ml)*0.1D0
	      IF(xy.LT.0.0D0)xy=0.0D0
		  xz=(1.0D0-xy)*0.5D0
	      CALL Dweimer96(IYR,IMO,IDA,IHR,IMN,XMLAL,XMLTL,Hindex(1),
     &                   Hindex(2),Hindex(3),NFi,EPT3,ET3,EP3,EPO3)
	      EEM(k,1:NFi,1)=EEM(k,1:NFi,1)*xy+xz*EP3(1:NFi)*RX1*XX
	      EEM(k,1:NFi,2)=EEM(k,1:NFi,2)*xy+xz*ET3(1:NFi)*RX2/SIi(k)*XX
            XMLAL(1:NFi)=0.0D0-Ml
	      CALL Dweimer96(IYR,IMO,IDA,IHR,IMN,XMLAL,XMLTL,Hindex(1),
     &                     Hindex(2),Hindex(3),NFi,EPT3,ET3,EP3,EPO3)
	      EEM(k,1:NFi,1)=EEM(k,1:NFi,1)+xz*EP3(1:NFi)*RX1*XX
	      EEM(k,1:NFi,2)=EEM(k,1:NFi,2)+xz*ET3(1:NFi)*RX2/SIi(k)*XX
	    ENDIF
	  ENDDO
      END SUBROUTINE EEmodel
c----------------------------------------------------------------------
	SUBROUTINE EfieldTrans(ed1,ed2,dSMLON,HB)
      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/input/ F107,F107A,AP,Uth,TimeStep,Hindex(6),iyd,IYR,IMO,IDA
	COMMON/EEmC/  EEm(NLi,NFi,2),EEH(NHN,2,2),dE(NLi,NFi,2)
      real*8 ed1(nmlat,nmlonp1),ed2(nmlat,nmlonp1),DyP(2),REP(2),XMLA(1)
     &       ,XMLT(1),ET1(1),EP1(1),EPO1(NHN)

	n1=2
      DO i=1,nmlat
	  DO j=1,nmlon
	    IF(ie(i,j,7).EQ.1)THEN
	      ix1=ie(i,j,1)
	      IF(ix1.GT.NLi)ix1=lmi+1-ix1
	      ix2=ie(i,j,2)
	      IF(ix2.GT.NLi)ix2=lmi+1-ix2
            DyP(1:n1)=EEm(ix1,ie(i,j,3),1:n1)*xie(i,j,1)
     &               +EEm(ix1,ie(i,j,4),1:n1)*xie(i,j,2)
     &               +EEm(ix2,ie(i,j,3),1:n1)*xie(i,j,3)
     &               +EEm(ix2,ie(i,j,4),1:n1)*xie(i,j,4)
          ELSE IF(ie(i,j,7).EQ.2)THEN
	      ix1=ie(i,j,1)
	      IF(ix1.GT.NLi)ix1=lmi+1-ix1
            DyP(1:n1)=EEm(ix1,ie(i,j,2),1:n1)*xie(i,j,1)
     &               +EEm(ix1,ie(i,j,3),1:n1)*xie(i,j,2)
     &               +EEH(ie(i,j,4),ie(i,j,6),1:n1)*xie(i,j,3)
     &               +EEH(ie(i,j,5),ie(i,j,6),1:n1)*xie(i,j,4)
          ELSE IF(ie(i,j,7).EQ.3)THEN
            DyP(1:n1)=EEH(ie(i,j,1),ie(i,j,6),1:n1)*xie(i,j,1)
     &               +EEH(ie(i,j,2),ie(i,j,6),1:n1)*xie(i,j,2)
     &               +EEH(ie(i,j,3),ie(i,j,6),1:n1)*xie(i,j,3)
     &               +EEH(ie(i,j,4),ie(i,j,6),1:n1)*xie(i,j,4)
          ELSE IF(ie(i,j,7).EQ.4)THEN
            REP(1:n1)=(EEH(1,ie(i,j,6),1:n1)+EEH(2,ie(i,j,6),1:n1)
     &               +EEH(3,ie(i,j,6),1:n1))/3.0D0
            DyP(1:n1)=EEH(ie(i,j,1),ie(i,j,6),1:n1)*xie(i,j,1)
     &               +EEH(ie(i,j,2),ie(i,j,6),1:n1)*xie(i,j,2)
     &               +REP(1:n1)*xie(i,j,3)
	    ENDIF
          ed1(i,j)=DyP(1)
          ed2(i,j)=DyP(2)
	  ENDDO
        ed1(i,nmlonp1)=ed1(i,1)
        ed2(i,nmlonp1)=ed2(i,1)
	ENDDO
	END SUBROUTINE EfieldTrans
c------------------------------------------------------------------------
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCCCCCCC               ISR ion drift model                CCCCCCCCCCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
	SUBROUTINE dEdrift(MLON,MLAT,iDAYNO,SLT,Hmin,Hmax,Vup,VWE)
	REAL*8 MLON,MLAT,SLT,Hmin,Hmax,Vup(2),VWE(2)
	REAL*4 rMLON,rMLAT,rDAYNO,rSLT,rHmin,rHmax,rVup(2),rVWE(2)
	  rMLON=real(MLON)
	  rMLAT=real(MLAT)
	  rDAYNO=real(iDAYNO)
	  rSLT=real(SLT)
	  rHmin=real(Hmin)
	  rHmax=real(Hmax)
        CALL Edrift(rMLON,rMLAT,rDAYNO,rSLT,rHmin,rHmax,rVup,rVWE)
	  Vup(1:2)=dble(rVup(1:2))
	  VWE(1:2)=dble(rVWE(1:2))
    	  RETURN
	END
C 
      SUBROUTINE Edrift(XMLON,GMLATs,DAYNO,XLT,Hmin,Hmax,Vup,VWE)
      PARAMETER(MAXLAT=2)
      DIMENSION ALT(MAXLAT),GMLAT(MAXLAT)
      DIMENSION Vup(MAXLAT),VWE(MAXLAT)
C     Hmin=100,Hmax=500
      ALT300=300.e3  !in Meters
         DO iMLAT=1,MAXLAT
            ALT(iMLAT)=Hmin+(Hmax-Hmin)*real(iMLAT-1)/real(MAXLAT-1)
            GMLAT(iMLAT)=GMLATs  !in degree
         END DO
      ISEASAV=0
      IUTAV=0
         UT=XLT-(XMLON-69.0)/15.0
100      DO iMLAT=1,MAXLAT
            XMLAT=GMLAT(iMLAT)
            CLAT=90.0-XMLAT
            CALL DipoleLine(CLAT,ALT(iMLAT)*1000.0,ALTeq,-1)
            IF(ALTeq.LT.ALT300) ALTeq=ALT300
            CALL DipoleLine(CLAT,-ALT300,ALTeq,1)
            XMLAT=90.0-CLAT
            CALL ISREFIELD(XMLAT,XMLON,DAYNO,UT,ISEASAV,IUTAV,POT,VU,VE)
c            CALL ToEquator(VUeq,VEeq,Vu,Ve,CLAT)
c            CALL FromEquator(VUeq,VEeq,Vu,Ve,CLAT)
            Vup(iMLAT)=Vu 
            VWE(iMLAT)=Ve
         ENDDO
      RETURN
	END
C
      SUBROUTINE DipoleLine(CLAT,Alt,ALTeq,MODE)
C  Gives Magnetic coordinats at a dipole field model;
C  if MODE>0, Convert coordinates from the equator to CLAT;
C             if ALT<0, then gives CLAT;
C             else gives ALT; 
C  else, Convert coordinates from CLAT to the equator;
C  ALT in meters, CLAT in degree;
C  r = Req * sin(CLAT)**2
      DegRad=4.0*atan(1.0)/180.0
      Re=6.37e6
      IF(MODE.GT.0) THEN
        Req=Re+ALTeq
        IF(Alt<0.0) THEN
          R=Re-Alt
          CLAT=asin(sqrt(R/Req))/DegRad
        ELSE
          Alt=Req*sin(CLAT*DegRad)**2-Re
        ENDIF
      ELSE
        R=Re+Alt
        ALTeq=R/sin(CLAT*DegRad)**2-Re
      ENDIF
      RETURN
      END
C
      SUBROUTINE FromEquator(VUeq,VEeq,Vu,Ve,CLAT)
C  Gives Vertical drift components at colatitude CLAT(degree)
C  from those at the magnetic equator;
      DegRad=4.0*atan(1.0)/180.0
      Vu=VUeq*sin(CLAT*DegRad)**3/sqrt(1.0+3.0*cos(CLAT*DegRad)**2)
      Ve=VEeq*sin(CLAT*DegRad)**3	      
      RETURN
      END
C
      SUBROUTINE ToEquator(VUeq,VEeq,Vu,Ve,CLAT)
C  Gives Vertical drift components at the magnetic equator
C  from those at colatitude CLAT(degree);
      DegRad=4.0*atan(1.0)/180.0
      VUeq=VU*sqrt(1.0+3.0*cos(CLAT*DegRad)**2)/sin(CLAT*DegRad)**3
      VEeq=VE/sin(CLAT*DegRad)**3	      
      RETURN
	END
C
      SUBROUTINE ISREFIELD(XMLAT,XMLON,DAYNO,UT,ISEASAV,IUTAV,POT,VU,VE)
C  GIVES QUIET-DAY IONOSPHERIC ELECTROSTATIC PSEUDO-POTENTIAL AND E X B
C  DRIFTS AT 300 KM FOR SOLAR MINIMUM CONDITIONS. SEE RICHMOND ET AL. (JGR,
C  1980, P.4658) FOR DEFINITIONS OF MAGNETIC COORDINATES, PSEUDO-POTENTIAL,
C  AND DRIFT COMPONENTS.
C
C  INPUT PARAMETERS -
C    XMLAT, XMLON ARE MAGNETIC LATITUDE AND EAST LONGITUDE IN DEGREES.
C    DAYNO IS DAY NUMBER OF THE YEAR FROM 1. TO 365.24, WITH 1. BEING JAN. 1.
C    UT IS UNIVERSAL TIME IN HOURS.
C    ISEASAV IS 0 IF NO SEASONAL AVERAGING IS DESIRED.
C            IS 1 FOR AVERAGE OVER NOV. - FEB.
C            IS 2 FOR AVERAGE OVER MAY -AUG.
C            IS 3 FOR AVERAGE OVER MAR., APR., SEPT., OCT.
C            IS 4 FOR AVERAGE OVER ENTIRE YEAR.
C      IF ISEASAV.NE.0, DAYNO IS IGNORED.
C    IUTAV IS 0 IF NO UT AVERAGING IS DESIRED.
C          IS 1 FOR AVERAGE OVER ALL UT AT THE FIXED LOCAL TIME GIVEN BY
C                   UT + (XMLON - 69.)/15.
C
C  OUTPUT PARAMETERS -
C    POT IS THE ELECTROSTATIC PSEUDO-POTENTIAL IN VOLTS.
C    VU IS THE DRIFT VELOCITY COMPONENT PERPENDICULAR TO THE GEOMAGNETIC FIELD
C      IN THE UPWARD/POLEWARD DIRECTION IN THE MAGNETIC MERIDIAN PLANE, IN M/S.
C    VE IS THE DRIFT VELOCITY COMPONENT IN THE MAGNETIC EASTWARD DIRECTION, IN
C      M/S.
C  THE OUTPUT VALUES ARE GEOPHYSICALLY MEANINGFUL ONLY FOR LATITUDES BETWEEN
C    ABOUT -65 AND +65 DEGREES.  IF ISEASAV OR IUTAV IS OUT OF RANGE POT, VU,
C    AND VE ARE SET TO -1/0.
      DIMENSION KF(128),LF(128),MF(128),NF(128),JF(128),Q(5),RS(16),
     1   FUT(5),RR(16),RSRS(16),P(16),VP(16),PA(5,9),VA(5,9),FS(3),
     2   FT(3,3),SML(4),CML(4),A(128),PB(9),VB(9)
      DATA IFIRST/1/,ISVP/5/,A/
     * -70.,-183.,  31.,-112.,  19., -39.,  -2.,   2., -33.,   2.,   2.,
     *-111.,  46.,  -4.,  -5.,   7.,   9., -17.,   2.,   9., -10.,   2.,
     *  -9.,  22., 145., -57., -42.,  -6.,   6.,  -5.,  -2.,  20.,  16.,
     *  16., -77., -18.,  13.,  -8.,  16., -52., -10.,   7.,   2.,  11.,
     * -28.,   2., -85., -82.,   3.,-281., -71., -25., -57., -50.,  21.,
     * -10.,  10., -81.,  24.,   7.,   5.,  30.,  32.,   5.,  -5.,  11.,
     * -31.,   8.,  10.,  20., -15., -42.,  32.,   7., -19.,   7.,  34.,
     * -11., -15.,  26.,  21.,   1.,  22.,  12.,  -2., 275., 777.,-318.,
     *-320.,-208.,  47., 429.,-523.,   8., -35.,-224.,-450., -66.,  -7.,
     *  -8.,-231.,  55.,   6., -28., -51., -81.,  48.,   9.,   2., -10.,
     *  54.,  16., 112.,  69., -33., 120., -47.,   5., -19., -17., -23.,
     * -40., -22., -21.,  -7., -30.,  15.,   3./
C***********************************************************************
C88/2/2 DEACTIVATE SHORTCUT, SINCE IT WON'T NECESSARILY WORK RIGHT IF
C  VARIABLES ARE NOT RETAINED WITH A SAVE STATEMENT.
C     IF (IFIRST.EQ.0) GO TO 510
C      IFIRST = 0
C***********************************************************************
C  (THROUGH STATEMENT 510) SET UP CONSTANT PARAMETERS IN FIRST CALL TO
C    SUBROUTINE.
C  SET UP VALUES OF XMLATP, TUP, TLP, AND DAYNOP WHICH ARE NOT EQUAL TO XMLAT,
C    UT, MAGNETIC LOCAL TIME, AND DAYNO, RESPECTIVELY.
      XMLATP = -361.
      IF (XMLATP.EQ.XMLAT) XMLATP = 0.
      TUP = -25.
      IF (TUP.EQ.UT) TUP = 0.
      TL = UT + (XMLON - 69.)/15.
      TLP = -25.
      IF (TLP.EQ.TL) TLP = 0.
      DAYNOP = -366.
      IF (DAYNOP.EQ.DAYNO) DAYNOP = 0.
C  (THROUGH STATEMENT 100) SELECT ONLY THOSE TERMS IN SERIES FOR WHICH
C    COEFFICIENTS A ARE DEFINED TO BE NON-ZERO.
      I = 0
      DO 100 KP=1,3
      LST = 4 - KP
      LEND = 2 + KP
      DO 100 LP = LST,LEND
      L = LP - 3
      IF (MOD(KP+LP,2) .NE.0) GO TO 100
      DO 90 NP=2,8
      N = NP - 1
      IF (KP.NE.3.AND.N.GT.6) GO TO 90
      IF (IABS(L).EQ.2.AND.N.GT.5) GO TO 90
      DO 80 MP = 1,9
      M = MP - 5
      IF (IABS(M).GT.N) GO TO 80
      IF (IABS(M).GT.3.AND.IABS(L).EQ.2) GO TO 80
      IF (MOD(N-M,2).NE.0) GO TO 80
      I = I + 1
      KF(I) = KP
      LF (I)= LP
      NF(I) = NP
      MF(I) = MP
   80 CONTINUE
   90 CONTINUE
  100 CONTINUE
      IMAX = I
      FT(1,1) = .75*SQRT(6.E0)/3.1415926535898
      FT(1,2) = 2.E0*FT(1,1)
      FT(1,3) = 1.E0
      FT(2,1) = FT(1,1)
      FT(2,2) = -FT(1,2)
      FT(2,3) = 1.E0
      FT(3,1) = FT(2,2)
      FT(3,2) = 0.
      FT (3,3) = 1.E0
      HRANG= 3.1415926535898/12.
      DANG = 3.1415926535898/182.62
      SQ2 = SQRT(2.E0)
      RAD = 180./3.1415926535898
      RB = -6.671E6*5.2E-5
C RB IS -(EARTH RADIUS + 3.E5 M) TIMES DIPOLE MAGNETIC FIELD AT POLE AT 300 KM.
      DO 400 I=1,IMAX
      MM = IABS(MF(I) - 5)
C  JF GIVES APPROPRIATE INDEX OF LEGENDRE POLYNOMIALS AS ORDERED BETWEEN
C    STATEMENTS 530 AND 595.
  400 JF(I) = (2*(NF(I) + 7*MM) - (MM - 1)**2 + 4)/4
C  (THROUGH STATEMENT 500) COMPUTE RR (DEFINED AS R(N,M)*R(N-1,M)) AND RSRS
C    (DEFINED AS R(N-1,M)**2 + R(N-2,M)**2) NEEDED FOR LEGENDRE POLYNOMIAL
C    GENERATING RECURSION RELATIONS, WHERE R(N,M) IS DEFINED AS
C    SQRT(N**2 - M**2)/SQRT(4*N**2 - 1).  ORDERING IS SAME AS FOR P(N,M).
      J = 0
      DO 500 MP=1,5
      M = MP - 1
      XM = M
      IF (M.NE.0) Q(MP) = SQRT((2.*XM + 1.)/(2.*XM))
      DO 500 NP = MP,8,2
      N = NP - 1
      XNS = N*N
      XNMS = (N-1)**2
      J = J + 1
      RS(J) = (XNS - XM*XM)/(4.*XNS - 1.)
      RSM    = (XNMS - XM*XM)/(4.*XNMS - 1.)
      RR(J) = SQRT(RS(J)*RSM)
      IF (NP.NE.MP) RSRS(J) = RSM    + RS(J - 1)
  500 CONTINUE
  510 CONTINUE
      IF (IUTAV.LT.0.OR.IUTAV.GT.1) GO TO 988
      IF (IABS(ISEASAV-2).GT.2) GO TO 988
C  (THROUGH STATEMENT 530) IF SEASONAL HARMONIC AMPLITUDES ARE DIFFERENT FROM
C    PREVIOUS CALL TO THIS SUBROUTINE, RECOMPUTE THEM.
C***********************************************************************
C88/2/2 DEACTIVATE SHORTCUT, SINCE IT WON'T NECESSARILY WORK RIGHT IF
C  VARIABLES ARE NOT RETAINED WITH A SAVE STATEMENT.
C     ICPT = 0
C     IF (ISEASAV.EQ.ISVP) GO TO 524
C***********************************************************************
      ICPT = 1
      ISVP = ISEASAV
C  (THROUGH STATEMENT 530) FS(1), FS(2), FS(3) ARE FACTORS FOR AMPLITUDE OF
C    SEMIANNUAL, ANNUAL, YEARLY AVERAGE COMPONENTS, RESPECTIVELY.
C  IF ISEASAV = 0, COMPUTE FS FOR GIVEN DAY OF THE YEAR.
      IF (ISEASAV.EQ.0)GO TO 525
C  IF ISEASAV = 4, USE ONLY YEARLY AVERAGE COMPONENT.
      IF (ISEASAV.EQ.4) GO TO 527
C  IF ISEASAV = 1 - 3, COMPUTE FS FOR APPROPRIATE SEASONAL AVERAGE.
      DO 521 K=1,3
  521 FS(K) = FT(ISEASAV,K)
      GO TO 530
  524 CONTINUE
      IF (ISEASAV.NE.0) GO TO 530
C***********************************************************************
C88/2/2 DEACTIVATE SHORTCUT, SINCE IT WON'T NECESSARILY WORK RIGHT IF
C  VARIABLES ARE NOT RETAINED WITH A SAVE STATEMENT.
C     IF (DAYNO.EQ.DAYNOP) GO TO 530
C***********************************************************************
  525 DAYNOP = DAYNO
      ICPT = 1
      ANG = (DAYNO + 9.)*DANG
      FS(1) = SQ2*COS(2.*ANG)
      FS(2) = SQ2*COS(ANG)
      FS(3) = 1.
      GO TO 530
  527 FS(1) = 0.
      FS(2) = 0.
      FS(3) = 1.
  530 CONTINUE
C  IF MAGNETIC LATITUDE IS SAME AS IN PREVIOUS CALL TO THIS SUBROUTINE, SKIP TO
C    STATEMENT 596.
C***********************************************************************
C88/2/2 DEACTIVATE SHORTCUT, SINCE IT WON'T NECESSARILY WORK RIGHT IF
C  VARIABLES ARE NOT RETAINED WITH A SAVE STATEMENT.
C     IF (XMLAT.EQ.XMLATP) GO TO 596
C***********************************************************************
      ICPT = 1
      XMLATP = XMLAT
      CT = SIN(XMLAT/RAD)
      CTS = CT*CT
      ST = SQRT(1. - CTS)
      RBT = RB*SQRT(.25 + .75*CTS)
C  (THROUGH STATEMENT 595) CALCULATE LEGENDRE POLYNOMIALS P(N,M) AS WELL AS VP,
C    DEFINED AS (DP(N,M)/D(COLATITUDE))/(RB*CT).
      J = 0
      DO 595 MP = 1,5
C  FOR MP=1, P=P(N,M).  FOR MP.GT.1, P=P(N,M)/ST.  THIS DIFFERENCE IS SO THAT
C    PROGRAM NEVER DIVIDES BY ST.
      J = J + 1
      XM = MP - 1
      MPP = MP + 2
      IF (MP.GT.1) GO TO 550
      X = 1.
      P(1) = 1.E0
      XZ = 2.*ST/RB
      VP(1) = 0.
      GO TO 560
  550 P(J) = Q(MP)*X
      X = P(J)*ST
      VP(J) = XM*P(J)/RB
      XZ = 2.*ST*ST/RB
  560 DO 590 NP = MPP,8,2
      J = J + 1
      Z = 0.
      ZP = 0.
      IF (NP.EQ.MPP) GO TO 570
      Z = RR(J-1)*P(J-2)
      ZP = RR(J-1) *VP(J-2)
  570 P(J) = ((CTS - RSRS(J))*P(J-1) - Z)/RR(J)
  590 VP(J) = ((CTS - RSRS(J))*VP(J-1) - ZP - XZ*P(J-1))/RR(J)
  595 CONTINUE
  596 CONTINUE
C  (THROUGH STATEMENT 600) CALCULATE ARRAYS OF FOURIER COEFFICIENTS, WITH LP
C    INDICATING HARMONIC OF UT AND MPF INDICATING HARMONIC OF MAGNETIC LOCAL
C    TIME.  IF ARRAYS ARE SAME AS IN PREVIOUS CALL TO THIS SUBROUTINE, SKIP TO
C    STATEMENT 601.
C***********************************************************************
C88/2/2 DEACTIVATE SHORTCUT, SINCE IT WON'T NECESSARILY WORK RIGHT IF
C  VARIABLES ARE NOT RETAINED WITH A SAVE STATEMENT.
C     IF (ICPT.EQ.0) GO TO 601
C***********************************************************************
      DO 597 MPF = 1,9
      DO 597 LP=1,5
      PA(LP,MPF) = 0.
  597 VA(LP,MPF) = 0.
      DO 600 I=1,IMAX
      LP = LF(I)
      MPF = MF(I)
      J = JF(I)
      X = A(I)*FS(KF(I))
      PA(LP,MPF) = PA(LP,MPF) + X*P(J)
      VA(LP,MPF) = VA(LP,MPF) + X*VP(J)
  600 CONTINUE
  601 CONTINUE
C  (THROUGH STATEMENT 607) CALCULATE FOURIER COEFFICIENTS PB AND VB AT GIVEN UT
C    FOR HARMONICS OF TL.
C  IF UT IS SAME AS IN PREVIOUS CALL TO THIS SUBROUTINE, SKIP TO STATEMENT 603.
C***********************************************************************
C88/2/2 DEACTIVATE SHORTCUT, SINCE IT WON'T NECESSARILY WORK RIGHT IF
C  VARIABLES ARE NOT RETAINED WITH A SAVE STATEMENT.
C     IF (UT.EQ.TUP) GO TO 603
C***********************************************************************
      TUP = UT
      ICPT = 1
      TUA = UT*HRANG
      SL = SIN(TUA)
      CL = COS(TUA)
      FUT(3) = 1.
      FUT(2) = SQ2*CL
      FUT(4) = SQ2*SL
      FUT(1) = CL*FUT(2) - SL*FUT(4)
      FUT(5) = CL*FUT(4) + SL*FUT(2)
  603 CONTINUE
      IF (ICPT.EQ.0) GO TO 608
      DO 607 MPF = 1,9
      PB(MPF) = 0.
      VB(MPF) = 0.
      DO 605 LP=1,5
      IF (IUTAV.NE.0.AND.LP.NE.3) GO TO 605
      PB(MPF) = PB(MPF) + FUT(LP)*PA(LP,MPF)
      VB(MPF) = VB(MPF) + FUT(LP)*VA(LP,MPF)
  605 CONTINUE
  607 CONTINUE
C  TL IS MAGNETIC LOCAL TIME.
  608 TL = UT + (XMLON - 69.)/15.
C  IF TL IS SAME AS IN PREVIOUS CALL TO THIS SUBROUTINE, SKIP TO STATEMENT 630.
C***********************************************************************
C88/2/2 DEACTIVATE SHORTCUT, SINCE IT WON'T NECESSARILY WORK RIGHT IF
C  VARIABLES ARE NOT RETAINED WITH A SAVE STATEMENT.
C     IF (TL.EQ.TLP) GO TO 630
C***********************************************************************
      TLP = TL
C  (THROUGH STATEMENT 610) CALCULATE SINES AND COSINES, TIMES SQ2, OF HARMONICS
C    OF MAGNETIC LOCAL TIME.
      TLA = TL*HRANG
      SL = SIN(TLA)
      CL = COS(TLA)
      SML(1) = SQ2*SL
      CML(1) = SQ2*CL
      DO 610 M=2,4
      SML(M) = CL*SML(M-1) + SL*CML(M-1)
  610 CML(M) = CL*CML(M-1) - SL*SML(M-1)
C  CALCULATE POT, VE, AND VU BY SUMMING FOURIER COEFFICIENTS MULTIPLIED BY
C    APPROPRIATE SINES AND COSINES.
  630 POT = 0.
      VE = 0.
      VU = 0.
      DO 640 M=1,4
      MPF = M+5
      MMF = 5 - M
      XM = M
      POT = POT + PB(MPF)*SML(M) + PB(MMF)*CML(M)
      VE = VE + VB(MPF)*SML(M) +   VB(MMF)*CML(M)
  640 VU = VU + XM*(PB(MPF)*CML(M) - PB(MMF)*SML(M))
      POT = POT*ST + PB(5)
      VE = VE + VB(5)
      VU = VU/RBT
      GO TO 999
  988 POT = -1/IFIRST
      VU = -1/IFIRST
      VE = -1/IFIRST
  999 RETURN
      END
C
      SUBROUTINE D1_Interpolate(X,Y,MAX,N,t,z)
C
      REAL X(MAX),Y(MAX)
C
      z=0.0
      IF(N.LT.1) RETURN
      IF(N.EQ.1) THEN
        z=Y(1)
        RETURN
      ENDIF
      IF(N.EQ.2) THEN
        z=(Y(1)*(t-X(2))-Y(2)*(t-X(1)))/(X(1)-X(2))
        RETURN
      ENDIF
      DO i=1,N-1
         IF(x(i).GE.t) GOTO 10
      ENDDO
10    k=i-4
      IF(k.LT.0) k=0
      m=i+3
      IF(m.GT.N-1) m=N-1
      DO i=k+1,m+1
         s=1.0
         DO j=k+1,m+1
            IF(i.NE.j) s=s*(t-X(j))/(X(i)-X(j))
         ENDDO
         z=z+s*Y(i)
      ENDDO
      RETURN
      END
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCCCC               THE END OF ISR drift model               CCCCCCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCCCC     Weimer96 polar cap ionospheric electric fields     CCCCCCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      SUBROUTINE Dweimer96(IYR,IMO,IDA,IHR,IMN,XMLA,XMLT,rBY,rBZ,rSWS,
     &	       NmN,Eps,ET,EP,EPOT)
c IYR,IMO,IDA,IHR,IMN means year,month,day,hour,minute; weimer96 was define in 
C magnetic coordinate, so RMLA(deg) means latitude, and MLT means magnetic 
C local time used to repersent longitude. By,Bz:IMF y,z component GSM coord
C (nT); plus toward dusk and Mag north.sws: solar wind speed in GSM (km/sec)
      INTEGER IYR,IMO,IDA,IHR,IMN
      REAL*8 XMLA(3360),XMLT(3360),rBY,rBZ,rSWS,ET(3360),EP(3360),EPOT
     &       (3360),Eps,XX,tmax,tmin
      REAL*4 RMLA,RMLT,BY,BZ,SWS
        BY=real(rBY)
	  BZ=real(rBZ)
	  SWS=real(rSWS)
	  IF (XMLA(1).LT.0) THEN
	    IHEM=-1
	  ELSE
	    IHEM=1
	  ENDIF
        Tmin=0.0D0
	  Tmax=0.0D0
	  R2D=57.2957795130823208767981548147
C ISETM=1 ! flag for WEIEPOT96 indicating new geophysical conditions        
	  BT=SQRT(BY*BY+BZ*BZ)
	  ANGL=ATAN2(BY,BZ)*R2D
	  HR=real(IHR)
	  TILT=GET_TILT(IYR,IMO,IDA,HR)
	  ANGL=IHEM*ANGL
	  TILT=IHEM*TILT
	  CALL SETMODEL96 (ANGL,BT,TILT,SWS)
        DO i=1,NmN
          RMLA=abs(real(XMLA(i)))
		RMLT=real(XMLT(i))
          EPOT4 =EPOTVAL96 (RMLA,RMLT)        ! potential
          CALL GECMP96 (RMLA,RMLT,ET4,EP4)    ! electric field
C     ET,EP= Etheta,Ephi(magnetic equatorward,eastward) E field component (V/m)
	    EP(i)=dble(EP4)
		ET(i)=dble(ET4)
		EPOT(i)=dble(EPOT4)*1.D3
	    if(EPOT(i).ge.tmax)tmax=EPOT(i)
		if(EPOT(i).le.tmin)tmin=EPOT(i)
	  ENDDO
	  Eps=Dabs(dble(tmax)-dble(tmin))*1e-3
      END SUBROUTINE Dweimer96
C If the next two lines begin with tilde (~), see the INSTALLATION notes
************************ Copyright 1996, Dan Weimer/MRC ***********************
* Subroutines to calculate the electric potentials from the Weimer '96 model of
* the polar cap ionospheric electric potentials.
*
* To use, first call subroutine ReadCoef once.
* Next, call SetModel with the specified input parameters.
* The function EpotVal(gLAT,gMLT) can then be used repeatively to get the
* electric potential at the desired location in geomagnetic coordinates.
*
* This code is protected by copyright and is
* distributed for research or educational use only.
* Commerical use without written permission from Dan Weimer/MRC is prohibited.
C NCAR      Revisions for use at NCAR:
C            (1) Change behavior at minimum magnetic latitude.  When approaching
C                the model minimum absolute latitude (45 degrees) the electric
C                potential returned used to go to zero discontinuously; although
C                intended as a flag, it created artificial gradients in the
C                electric field calculation.  Now the potential returned is
C                constant (that of the minimum latitude) for any latitude at
C                or equatorward of the minimum.
C            (2) Accommodate running sumultaneously 1996 and 2001 versions.  To
C                avoid name collisions this required: (i) revising names (e.g.,
C                adding '96') for differing subprograms, and (ii) relocating
C                common routines into another file (weicom.f).
C            (3) Pass the coefficients file name and unit number into RDCOEF96
C                rather than using hard coded values.
C            (4) Add wrapper subroutines for non-ANSI trig functions which
C                input angles in degrees.
C            (5) Add electric field routine (GECMP96) to deterine the electric
C                potential gradient.
C            (6) Add wrapper routine (WEIEPOT96) for use with AMIE; this is a
C                substitute for calling SETMODEL96 and EPOTVAL96.
C NCAR      NCAR changes are delimited by "CNCAR"
************************ Copyright 1996, Dan Weimer/MRC ***********************
CNCAR      Oct 01:  Distinguish from the 2001 version
	FUNCTION EPOTVAL96 (gLAT,gMLT)
C       FUNCTION EPOTVAL   (gLAT,gMLT)
CNCAR   Return the value of the electric potential in kV at
* corrected geomagnetic coordinates gLAT (degrees) and gMLT (hours).
*
* Must first call ReadCoef and SetModel to set up the model coeficients for
* the desired values of Bt, IMF clock angle, Dipole tilt angle, and SW Vel.
*
	REAL gLAT,gMLT
	Real Theta,Phi,Z,ct,Phim
	REAL Plm(0:20,0:20)
 
	REAL Coef(0:1,0:8,0:3),pi
	INTEGER ML,MM
	COMMON/SetCoef/ML,MM,Coef,pi
 
	r = 90.-gLAT
CNCAR      Sep 01: Limit the colatitude (r) such that an input absolute
C          magnetic latitude of 45 degrees or less always calculates the
C          potential at 45.
	R = AMIN1 (R,45.)
C       IF(r .LT. 45.)THEN
CNCAR
	  Theta=r*pi/45.
          Phi=gMLT*pi/12.
	  Z=Coef(0,0,0)
	  ct=COS(Theta)
	  CALL Legendre(ct,ML,MM,Plm)
	  DO l=1,ML
	    Z=Z + Coef(0,l,0)*Plm(l,0)
	    IF(l.LT.MM)THEN
	      limit=l
	    ELSE
	      limit=MM
	    ENDIF
	    DO m=1,limit
	      phim=phi*m
	      Z=Z + Coef(0,l,m)*Plm(l,m)*COS(phim) +
     &		   Coef(1,l,m)*Plm(l,m)*SIN(phim)
	    ENDDO
	  ENDDO
CNCAR
C       ELSE
C         Z=0.
C       ENDIF
	EPOTVAL96 = Z
C       EpotVal   = Z
CNCAR
	RETURN
	END 
************************ Copyright 1996, Dan Weimer/MRC *********************** 
CNCAR      Sep 01:  Pass in the unit no. and file name rather than hard coding
C          inside.  Also change name to distinguish 2001 version.
 	SUBROUTINE READCOEF96 (udat,cfile)
C       SUBROUTINE ReadCoef
CNCAR
* Read in the data file with the model coefficients
	INTEGER udat
CNCAR
	CHARACTER    cfile*(*), skip*15
C       CHARACTER*15 cfile,skip
CNCAR
	REAL C(0:3)
	REAL Cn( 0:3 , 0:1 , 0:4 , 0:1 , 0:8 , 0:3 )
	INTEGER MaxL,MaxM,MaxN
	COMMON /AllCoefs/MaxL,MaxM,MaxN,Cn
 
 
CNCAR      Jan 97:  Initialize constants used in GECMP96
C          Sep 01:  Omit unneeded min lat variables because of switch to
C                   constant potential at and below min lat (hardcoded
C                   in EPOTVAL96).
      COMMON /CECMP/ ALAMX,STPD,STP2,CSTP,SSTP
C            ALAMX = Absolute max latitude (deg) for normal gradient calc.
C            STPD  = Angular dist (deg) of step @ 300km above earth (r=6371km)
C            STP2  = Denominator in gradient calc
      PARAMETER ( D2R =  0.0174532925199432957692369076847 ,
     +            R2D = 57.2957795130823208767981548147)
      STEP = 10.
      STPR = STEP/6671.
      STPD = STPR*R2D
      STP2 = 2.*STEP
      CSTP = COS (STPR)
      SSTP = SQRT (1. - CSTP*CSTP)
      ALAMX = 90. - STPD
CNCAR
 
CNCAR      Sep 01:  udat,cfile are now a formal arguments
C       cfile='wei96.cofcnts'
C       udat=99
CNCAR
	OPEN(udat,FILE=cfile,STATUS='OLD')
  900 FORMAT(A15)
 1000	FORMAT(3I8)
 2000	FORMAT(3I2)
 3000	FORMAT(2I2,4E15.6)
 
	READ(udat,900) skip
	READ(udat,1000) MaxL,MaxM,MaxN
	DO l=0,MaxL
	  IF(l.LT.MaxM)THEN
	    mlimit=l
	  ELSE
	    mlimit=MaxM
	  ENDIF
	  DO m=0,mlimit
	    IF(m.LT.1)THEN
	      klimit=0
	    ELSE
	      klimit=1
	    ENDIF
	    DO k=0,klimit
	      READ(udat,2000) ll,mm,kk
	      IF(ll.NE.l .OR. mm.NE.m .OR. kk.NE.k)THEN
		PRINT *,'Data File Format Error'
		STOP
	      ENDIF
	      DO n=0,MaxN
	        IF(n.LT.1)THEN
	          ilimit=0
	        ELSE
	          ilimit=1
	        ENDIF
		DO i=0,ilimit
		  READ(udat,3000) nn,ii,C
	          IF(nn.NE.n .OR. ii.NE.i)THEN
		    PRINT *,'Data File Format Error'
		    STOP
		  ENDIF
		  Cn(0,i,n,k,l,m)=C(0)
		  Cn(1,i,n,k,l,m)=C(1)
		  Cn(2,i,n,k,l,m)=C(2)
		  Cn(3,i,n,k,l,m)=C(3)
		ENDDO
	      ENDDO
	    ENDDO
	  ENDDO
	ENDDO
 
	CLOSE(udat)
	RETURN
	END
************************ Copyright 1996, Dan Weimer/MRC ***********************
CNCAR      Oct 01:  Change name from SetModel to be distinct
 
	SUBROUTINE SETMODEL96 (angle,Bt,Tilt,SWVel)
C       SUBROUTINE SetModel   (angle,Bt,Tilt,SWVel)
CNCAR
*
* Calculate the complete set of spherical harmonic coefficients,
* given an arbitrary IMF angle (degrees from northward toward +Y),
* magnitude Bt (nT), dipole tilt angle (degrees),
* and solar wind velocity (km/sec).
* Returns the Coef in the common block SetCoef.
*
	REAL angle,Bt,Tilt,SWVel
	REAL FSC(0:1,0:4)
	REAL Cn( 0:3 , 0:1 , 0:4 , 0:1 , 0:8 , 0:3 )
	INTEGER MaxL,MaxM,MaxN
	COMMON /AllCoefs/MaxL,MaxM,MaxN,Cn
 
	REAL Coef(0:1,0:8,0:3),pi
	INTEGER ML,MM
	COMMON/SetCoef/ML,MM,Coef,pi
 
	pi=2.*ASIN(1.)
	ML=MaxL
	MM=MaxM
	SinTilt=SIND(Tilt)
 
	omega=angle*pi/180.
	DO l=0,MaxL
	  IF(l.LT.MaxM)THEN
	    mlimit=l
	  ELSE
	    mlimit=MaxM
	  ENDIF
	  DO m=0,mlimit
	    IF(m.LT.1)THEN
	      klimit=0
	    ELSE
	      klimit=1
	    ENDIF
	    DO k=0,klimit
* Retrieve the regression coefficients and evaluate the function
* as a function of Bt,Tilt,and SWVel to get each Fourier coefficient.
	      DO n=0,MaxN
	        IF(n.LT.1)THEN
	          ilimit=0
	        ELSE
	          ilimit=1
	        ENDIF
		DO i=0,ilimit
		  FSC(i,n)=Cn(0,i,n,k,l,m) + Bt*Cn(1,i,n,k,l,m) +
     &		   SinTilt*Cn(2,i,n,k,l,m) + SWVel*Cn(3,i,n,k,l,m)
		ENDDO
	      ENDDO
* Next evaluate the Fourier series as a function of angle.
      	      Coef(k,l,m)=FSVal(omega,MaxN,FSC)
	    ENDDO
	  ENDDO
	ENDDO
	RETURN
	END
************************ Copyright 1996, Dan Weimer/MRC ***********************
	SUBROUTINE LEGENDRE(x,lmax,mmax,Plm)
* compute Associate Legendre Function P_l^m(x)
* for all l up to lmax and all m up to mmax.
* returns results in array Plm
* if X is out of range ( abs(x)>1 ) then value is returned as if x=1.
	DIMENSION Plm(0:20,0:20)
	  DO l=0,20
	    DO m=0,20
		Plm(l,m)=0.
	    ENDDO
	  ENDDO
	xx=MIN(x,1.)
	xx=MAX(xx,-1.)
	IF(lmax .LT. 0 .OR. mmax .LT. 0 .OR. mmax .GT. lmax )THEN
	  Print *,'Bad arguments to Legendre'
	  RETURN
	ENDIF
* First calculate all Pl0 for l=0 to l
	Plm(0,0)=1.
	IF(lmax.GT.0)Plm(1,0)=xx
	IF (lmax .GT. 1 )THEN
	  DO L=2,lmax
	    Plm(L,0)=( (2.*L-1)*xx*Plm(L-1,0) - (L-1)*Plm(L-2,0) )/L
	  ENDDO
	ENDIF
	IF (mmax .EQ. 0 )RETURN
	fact=SQRT( (1.-xx)*(1.+xx) )
	DO M=1,mmax
	  DO L=m,lmax
	    lm2=MAX(L-2,0)
	    Plm(L,M)=Plm(lm2,M) - ( 2*L-1)*fact*Plm(L-1,M-1)
	  ENDDO
	ENDDO
	RETURN
	END
 
CNCAR
      SUBROUTINE GECMP96 (AMLA,RMLT,ET,EP)
C          Get Electric field components for the Weimer electrostatic
C          potential model.  Before use, first load coefficients (CALL
C          READCOEF) and initialize model conditions (CALL SETMODEL).
C          The electric field, the gradient of the electic potential,
C          is determined using finite differences.
C
C          INPUTS:
C            AMLA = Absolute value of magnetic latitude (deg)
C            RMLT = Magnetic local time (hours).
C          RETURNS:
C            ET = Etheta (magnetic equatorward*) E field component (V/m)
C            EP = Ephi   (magnetic eastward)     E field component (V/m)
C
C          * ET direction is along the magnetic meridian away from the
C            current hemisphere; i.e., when ET > 0, the direction is
C              southward when in northern magnetic hemisphere
C              northward when in southern magnetic hemisphere
C
C          HISTORY:
C          Jan 97:  Initial implementation at NCAR.  R.Barnes
C          Feb 01:  Error corrected in determining DPHI.  Old version
C          used wrong spherical right triangle formula (inadvertently
C          computing along a latitude line); more importantly, it
C          neglected to convert mlt from degrees to hour angle input
C          to EPOTVAL96.
C          Sep 01:  Revised equatorward boundary logic to the scheme Barbara
C          uses for AMIE:  As absolute magnetic latitude decreases to the
C          model minimum, the electric potential from EPOTVAL96 now goes to
C          a non-zero constant (rather than zero), thus obviating the need
C          for special handling here.  Therefore, special logic for determining
C          gradients at lower limit has been removed.
 
      PARAMETER ( D2R =  0.0174532925199432957692369076847 ,
     +            R2D = 57.2957795130823208767981548147)
 
C          CECMP contains constants initialized in READCOEF
      COMMON /CECMP/ ALAMX,STPD,STP2,CSTP,SSTP
 
      ET = -99999.
      EP = -99999.
      IF (AMLA .LT. 0.) GO TO 100
 
C          Calculate -(latitude gradient) by stepping along the magnetic
C          latitude line in each direction (flipping coordinates when
C          going over pole to keep lat <= 90).
      KPOL  = 0
      XMLT  = RMLT
   10 XMLT1 = XMLT
      AMLA1 = AMLA + STPD
      IF (AMLA1 .GT. 90.) THEN
	AMLA1 = 180. - AMLA1
	XMLT1 = XMLT1 + 12.
      ENDIF
      P1 = EPOTVAL96 (AMLA1    ,XMLT1)
      P2 = EPOTVAL96 (AMLA-STPD,XMLT )
      IF (KPOL .EQ. 1) GO TO 20
      ET = (P1 - P2) / STP2
 
C          Calculate -(lon gradient).  Step along the magnetic meridion
C          in both directions to obtain the electric potential
      IF (AMLA .LT. ALAMX) THEN
	AMLA1 = ASIN(SIN(AMLA*D2R)*CSTP)
 
C       DPHI  = ASIN (SSTP/SIN(AMLA1))*R2D          ! old and wrong! (changed Feb 01)
	DPHI  = ASIN (SSTP/COS(AMLA1))*R2D/15.      ! 15 converts from degrees to hours
 
	AMLA1 = AMLA1*R2D
	P1 = EPOTVAL96 (AMLA1,XMLT+DPHI)
	P2 = EPOTVAL96 (AMLA1,XMLT-DPHI)
      ELSE
C          At the pole.  Avoid a divide by zero by using Art's trick
C          where Ephi(90,lon) = Etheta(90,lon+90).
	AMLA = 90.
	XMLT = XMLT + 6.
	KPOL = 1
	GO TO 10
      ENDIF
   20 EP = (P2 - P1) / STP2
      IF (KPOL .EQ. 1) EP = -EP
 
  100 RETURN
      END
 
      SUBROUTINE WEIEPOT96 (IYR,IMO,IDA,IHR,IMN,SWS,BY,BZ,IHEM,ISETM,
     +                     RMLA,RMLT, ET,EP,EPOT)
C          Interface to Weimer-96 for AMIE's combined electrostatic
C          potential models.  This replaces two calls (SETMODEL96 and
C          EPOTVAL96) and automates sign changes for southern hemisphere.
C          INPUTS:
C            IYR   = UT year ((4-digit)
C            IMO   = month of IYR
C            IDA   = day of IMO
C            IHR   = hour of day
C            IMN   = min of hour
C            SWS   = Solar wind speed (km/s)
C            BY    = IMF By component in GSM coordinates (nt)
C            BZ    = IMF Bz component in GSM coordinates (nt)
C            IHEM  = Hemisphere flag: (-1) southern, (1) northern
C            ISETM = Model conditions change flag: (0) no-change
C                                                  (1) time, IMF or SW changed
C            RMLA  = Magnetic latitude (deg) of point to determine Potential.
C                    RMLA should be positive (NH) only, since SH values are
C                    obtained by changing sign of By (or ANGL) and TILT.
C            RMLT  = Magnetic longitude (hrs) of point to determine Potential
C          RETURNS:
C            ET    = Etheta (magnetic equatorward*) E field component (V/m)
C            EP    = Ephi   (magnetic eastward)     E field component (V/m)
C            EPOT  = Electric potential (kV)
C
C          * ET direction is along the magnetic meridian away from the
C            current hemisphere; i.e., when ET > 0, the direction is
C              southward when in northern magnetic hemisphere
C              northward when in southern magnetic hemisphere
C
C          Since the AMIE model assumes a NH solution although the
C          latitudes are negative for SH data, the sign of ET
C          for the SH should be changed outside of this routine.
 
C          HISTORY:
C          Jan 97:  Initial implementation. B. Emery.
C          Feb 01:  Remove bug (HR was not defined)
C          Sep 01:  Add common block s.t. computed values are available to
C          the calling routine without changing the argument list
C          Oct 01:  Revise routine names s.t. this can run with the 2001 version
      COMMON /WEIAT/ ANGL, TILT
      PARAMETER (R2D=57.2957795130823208767981548147)
 
      H = REAL (IHEM)
 
      IF (ISETM .EQ. 1) THEN
	ANGL = ATAN2 (BY,BZ)*R2D
	BT   = SQRT (BY*BY + BZ*BZ)
	HR   = FLOAT(IHR) + FLOAT(IMN)/60.       ! define HR (bug fix Feb 01)
	TILT = GET_TILT (IYR,IMO,IDA,HR)
 
	HANGL = H * ANGL
	HTILT = H * TILT
	CALL SETMODEL96 (HANGL,BT,HTILT,SWS)
C       WRITE (6,'('' WEIMER-1996: sw by bz bt angle tilt='',6f7.2)')
C    +                            SWS,BY,BZ,BT,ANGL,TILT
      ENDIF
 
C          NH assumed latitudes only
      AMLA = ABS (RMLA)
      EPOT = EPOTVAL96 (AMLA, RMLT)
      CALL GECMP96 (AMLA,RMLT,ET,EP)
 
      RETURN
      END
CNCAR
C          Routines common to both 1996 and 2001 versions of Weimer's
C          Ionospheric Electrostatic Potential Model
 
************************ Copyright 1996,2001 Dan Weimer/MRC ***********************
	FUNCTION FSVal(omega,MaxN,FSC)
* Fourier Series Value
* Return value of Sine/Cosine Fourier series for N terms up to MaxN
* at angle omega, given the F.S. coeficients in array FSC
	REAL omega,FSC(0:1,0:*)
	INTEGER MaxN,n
	REAL Y,theta
	Y=0.
	DO n=0,MaxN
	  theta=omega*n
	  Y=Y + FSC(0,n)*COS(theta) + FSC(1,n)*SIN(theta)
	ENDDO
	FSVal=Y
	RETURN
	END
************************ Copyright 1996,2001 Dan Weimer/MRC ***********************
* COORDINATE TRANSFORMATION UTILITIES
 
CNCAR      Feb 01:  Changed TRANS to GET_TILT s.t. the dipole tilt angle is
C          returned.
 
	FUNCTION GET_TILT (YEAR,MONTH,DAY,HOUR)
C       SUBROUTINE TRANS(YEAR,MONTH,DAY,HOUR,IDBUG)
CNCAR
 
	INTEGER YEAR,MONTH,DAY,IDBUG
	REAL HOUR
C
C      THIS SUBROUTINE DERIVES THE ROTATION MATRICES AM(I,J,K) FOR 11
C      TRANSFORMATIONS, IDENTIFIED BY K.
C          K=1 TRANSFORMS GSE to GEO
C          K=2     "      GEO to MAG
C          K=3     "      GSE to MAG
C          K=4     "      GSE to GSM
C          K=5     "      GEO to GSM
C          K=6     "      GSM to MAG
C          K=7     "      GSE to GEI
C          K=8     "      GEI to GEO
C          K=9     "      GSM to SM
C	   K=10    "      GEO to SM
C	   K=11    "      MAG to SM
C
C      IF IDBUG IS NOT 0, THEN OUTPUTS DIAGNOSTIC INFORMATION TO
C      FILE UNIT=IDBUG
C	
	INTEGER GSEGEO,GEOGSE,GEOMAG,MAGGEO
	INTEGER GSEMAG,MAGGSE,GSEGSM,GSMGSE
	INTEGER GEOGSM,GSMGEO,GSMMAG,MAGGSM
	INTEGER GSEGEI,GEIGSE,GEIGEO,GEOGEI
	INTEGER GSMSM,SMGSM,GEOSM,SMGEO,MAGSM,SMMAG
 
	PARAMETER (GSEGEO= 1,GEOGSE=-1,GEOMAG= 2,MAGGEO=-2)
	PARAMETER (GSEMAG= 3,MAGGSE=-3,GSEGSM= 4,GSMGSE=-4)
	PARAMETER (GEOGSM= 5,GSMGEO=-5,GSMMAG= 6,MAGGSM=-6)
	PARAMETER (GSEGEI= 7,GEIGSE=-7,GEIGEO= 8,GEOGEI=-8)
	PARAMETER (GSMSM = 9,SMGSM =-9,GEOSM =10,SMGEO=-10)
	PARAMETER (MAGSM =11,SMMAG =-11)
C
C      The formal names of the coordinate systems are:
C	GSE - Geocentric Solar Ecliptic
C	GEO - Geographic
C	MAG - Geomagnetic
C	GSM - Geocentric Solar Magnetospheric
C	SM  - Solar Magnetic
C	
C      THE ARRAY CX(I) ENCODES VARIOUS ANGLES, STORED IN DEGREES
C      ST(I) AND CT(I) ARE SINES & COSINES.
C
C      Program author:  D. R. Weimer
C
C      Some of this code has been copied from subroutines which had been
C      obtained from D. Stern, NASA/GSFC.  Other formulas are from "Space
C      Physics Coordinate Transformations: A User Guide" by M. Hapgood (1991).
C
C      The formulas for the calculation of Greenwich mean sidereal time (GMST)
C      and the sun's location are from "Almanac for Computers 1990",
C      U.S. Naval Observatory.
C
	REAL UT,T0,GMSTD,GMSTH,ECLIP,MA,LAMD,SUNLON
 
	COMMON/MFIELD/EPOCH,TH0,PH0,DIPOLE
	COMMON/TRANSDAT/CX(9),ST(6),CT(6),AM(3,3,11)
C
	DATA EPOCH,TH0,PH0,DIPOLE/1980.,11.19,-70.76,.30574/
 
CNCAR      Feb 01:  Prevent debug printing to a file
	IDBUG = 0
CNCAR
 
	IF(YEAR.LT.1900)THEN
	  IYR=1900+YEAR
	ELSE
	  IYR=YEAR
	ENDIF
	UT=HOUR
	JD=JULDAY(MONTH,DAY,IYR)
	MJD=JD-2400001
	T0=(FLOAT(MJD)-51544.5)/36525.0
	GMSTD=100.4606184 + 36000.770*T0 + 3.87933E-4*T0*T0 +
     &        15.0410686*UT
	CALL ADJUST(GMSTD)
	GMSTH=GMSTD*24./360.
	ECLIP=23.439 - 0.013*T0
	MA=357.528 + 35999.050*T0 + 0.041066678*UT
	CALL ADJUST(MA)
	LAMD=280.460 + 36000.772*T0 + 0.041068642*UT
	CALL ADJUST(LAMD)
	SUNLON=LAMD + (1.915-0.0048*T0)*SIND(MA) + 0.020*SIND(2.*MA)
	CALL ADJUST(SUNLON)
	IF(IDBUG.NE.0)THEN
	  WRITE(IDBUG,*) YEAR,MONTH,DAY,HOUR
	  WRITE(IDBUG,*) 'MJD=',MJD
	  WRITE(IDBUG,*) 'T0=',T0
	  WRITE(IDBUG,*) 'GMSTH=',GMSTH
	  WRITE(IDBUG,*) 'ECLIPTIC OBLIQUITY=',ECLIP
	  WRITE(IDBUG,*) 'MEAN ANOMALY=',MA
	  WRITE(IDBUG,*) 'MEAN LONGITUDE=',LAMD
	  WRITE(IDBUG,*) 'TRUE LONGITUDE=',SUNLON
	ENDIF
 
	CX(1)= GMSTD
	CX(2) = ECLIP
	CX(3) = SUNLON
	CX(4) = TH0
	CX(5) = PH0
c Derived later:
c       CX(6) = Dipole tilt angle
c       CX(7) = Angle between sun and magnetic pole
c       CX(8) = Subsolar point latitude
c       CX(9) = Subsolar point longitude
 
	DO I=1,5
	  ST(I) = SIND(CX(I))
	  CT(I) = COSD(CX(I))
	ENDDO
C
      AM(1,1,GSEGEI) = CT(3)
      AM(1,2,GSEGEI) = -ST(3)
      AM(1,3,GSEGEI) = 0.
      AM(2,1,GSEGEI) = ST(3)*CT(2)
      AM(2,2,GSEGEI) = CT(3)*CT(2)
      AM(2,3,GSEGEI) = -ST(2)
      AM(3,1,GSEGEI) = ST(3)*ST(2)
      AM(3,2,GSEGEI) = CT(3)*ST(2)
      AM(3,3,GSEGEI) = CT(2)
C
      AM(1,1,GEIGEO) = CT(1)
      AM(1,2,GEIGEO) = ST(1)
      AM(1,3,GEIGEO) = 0.
      AM(2,1,GEIGEO) = -ST(1)
      AM(2,2,GEIGEO) = CT(1)
      AM(2,3,GEIGEO) = 0.
      AM(3,1,GEIGEO) = 0.
      AM(3,2,GEIGEO) = 0.
      AM(3,3,GEIGEO) = 1.
C
      DO I=1,3
      DO J=1,3
        AM(I,J,GSEGEO) = AM(I,1,GEIGEO)*AM(1,J,GSEGEI) +
     &    AM(I,2,GEIGEO)*AM(2,J,GSEGEI) +  AM(I,3,GEIGEO)*AM(3,J,GSEGEI)
      ENDDO
      ENDDO
C
      AM(1,1,GEOMAG) = CT(4)*CT(5)
      AM(1,2,GEOMAG) = CT(4)*ST(5)
      AM(1,3,GEOMAG) =-ST(4)
      AM(2,1,GEOMAG) =-ST(5)
      AM(2,2,GEOMAG) = CT(5)
      AM(2,3,GEOMAG) = 0.
      AM(3,1,GEOMAG) = ST(4)*CT(5)
      AM(3,2,GEOMAG) = ST(4)*ST(5)
      AM(3,3,GEOMAG) = CT(4)
C
      DO I=1,3
      DO J=1,3
       AM(I,J,GSEMAG) = AM(I,1,GEOMAG)*AM(1,J,GSEGEO) +
     &   AM(I,2,GEOMAG)*AM(2,J,GSEGEO) +  AM(I,3,GEOMAG)*AM(3,J,GSEGEO)
      ENDDO
      ENDDO
C
      B32 = AM(3,2,GSEMAG)
      B33 = AM(3,3,GSEMAG)
      B3  = SQRT(B32*B32+B33*B33)
      IF (B33.LE.0.) B3 = -B3
C
      AM(2,2,GSEGSM) = B33/B3
      AM(3,3,GSEGSM) = AM(2,2,GSEGSM)
      AM(3,2,GSEGSM) = B32/B3
      AM(2,3,GSEGSM) =-AM(3,2,GSEGSM)
      AM(1,1,GSEGSM) = 1.
      AM(1,2,GSEGSM) = 0.
      AM(1,3,GSEGSM) = 0.
      AM(2,1,GSEGSM) = 0.
      AM(3,1,GSEGSM) = 0.
C
      DO I=1,3
      DO J=1,3
        AM(I,J,GEOGSM) = AM(I,1,GSEGSM)*AM(J,1,GSEGEO) +
     &    AM(I,2,GSEGSM)*AM(J,2,GSEGEO) + AM(I,3,GSEGSM)*AM(J,3,GSEGEO)
      ENDDO
      ENDDO
C
      DO I=1,3
      DO J=1,3
        AM(I,J,GSMMAG) = AM(I,1,GEOMAG)*AM(J,1,GEOGSM) +
     &   AM(I,2,GEOMAG)*AM(J,2,GEOGSM) + AM(I,3,GEOMAG)*AM(J,3,GEOGSM)
      ENDDO
      ENDDO
C
	ST(6) = AM(3,1,GSEMAG)
	CT(6) = SQRT(1.-ST(6)*ST(6))
	CX(6) = ASIND(ST(6))
 
        AM(1,1,GSMSM) = CT(6)
        AM(1,2,GSMSM) = 0.
        AM(1,3,GSMSM) = -ST(6)
        AM(2,1,GSMSM) = 0.
        AM(2,2,GSMSM) = 1.
        AM(2,3,GSMSM) = 0.
        AM(3,1,GSMSM) = ST(6)
        AM(3,2,GSMSM) = 0.
        AM(3,3,GSMSM) = CT(6)
C
      DO I=1,3
      DO J=1,3
        AM(I,J,GEOSM) = AM(I,1,GSMSM)*AM(1,J,GEOGSM) +
     &    AM(I,2,GSMSM)*AM(2,J,GEOGSM) +  AM(I,3,GSMSM)*AM(3,J,GEOGSM)
      ENDDO
      ENDDO
C
      DO I=1,3
      DO J=1,3
        AM(I,J,MAGSM) = AM(I,1,GSMSM)*AM(J,1,GSMMAG) +
     &   AM(I,2,GSMSM)*AM(J,2,GSMMAG) + AM(I,3,GSMSM)*AM(J,3,GSMMAG)
      ENDDO
      ENDDO
C
      CX(7)=ATAN2D( AM(2,1,11) , AM(1,1,11) )
      CX(8)=ASIND( AM(3,1,1) )
      CX(9)=ATAN2D( AM(2,1,1) , AM(1,1,1) )
 
      IF(IDBUG.NE.0)THEN
	  WRITE(IDBUG,*) 'Dipole tilt angle=',CX(6)
	  WRITE(IDBUG,*) 'Angle between sun and magnetic pole=',CX(7)
	  WRITE(IDBUG,*) 'Subsolar point latitude=',CX(8)
	  WRITE(IDBUG,*) 'Subsolar point longitude=',CX(9)
 
        DO K=1,11
         WRITE(IDBUG,1001) K
         DO I=1,3
           WRITE(IDBUG,1002) (AM(I,J,K),J=1,3)
         ENDDO
        ENDDO
 1001   FORMAT(' ROTATION MATRIX ',I2)
 1002   FORMAT(3F9.5)
      ENDIF
 
CNCAR      Mar 96: return the dipole tilt from this function call.
      GET_TILT = CX(6)
CNCAR
 
      RETURN
      END
******************************************************************************
CNCAR      Feb 01:  Eliminate unused routines from translib.for: ROTATE,
C          ROTATEV, FROMCART, TOCART, MLT, MAGLONG, SUNLOC.  Remaining
C          are ADJUST and JULDAY
CNCAR
	SUBROUTINE ADJUST(ANGLE)
C	ADJUST AN ANGLE IN DEGREES TO BE IN RANGE OF 0 TO 360.
 10	CONTINUE
	IF(ANGLE.LT.0.)THEN
	  ANGLE=ANGLE+360.
	  GOTO 10
	ENDIF
 20	CONTINUE
 	IF(ANGLE.GE.360.)THEN
	  ANGLE=ANGLE-360.
	  GOTO 20
	ENDIF
	RETURN
	END
******************************************************************************
      FUNCTION JULDAY(MM,ID,IYYY)
      PARAMETER (IGREG=15+31*(10+12*1582))
      IF (IYYY.EQ.0) PAUSE 'There is no Year Zero.'
      IF (IYYY.LT.0) IYYY=IYYY+1
      IF (MM.GT.2) THEN
        JY=IYYY
        JM=MM+1
      ELSE
        JY=IYYY-1
        JM=MM+13
      ENDIF
      JULDAY=INT(365.25*JY)+INT(30.6001*JM)+ID+1720995
      IF (ID+31*(MM+12*IYYY).GE.IGREG) THEN
        JA=INT(0.01*JY)
        JULDAY=JULDAY+2-JA+INT(0.25*JA)
      ENDIF
      RETURN
      END
 
CNCAR      Routines added to work around non-ANSI trig functions which
C          input degrees instead of radians:  SIND, COSD, ASIND, ATAN2D
 
      FUNCTION SIND (DEG)
      PARAMETER ( D2R =  0.0174532925199432957692369076847 ,
     +            R2D = 57.2957795130823208767981548147)
      SIND = SIN (DEG * D2R)
      RETURN
      END
 
      FUNCTION COSD (DEG)
      PARAMETER ( D2R =  0.0174532925199432957692369076847 ,
     +            R2D = 57.2957795130823208767981548147)
 
      COSD = COS (DEG * D2R)
      RETURN
      END
 
      FUNCTION ASIND (RNUM)
      PARAMETER ( D2R =  0.0174532925199432957692369076847 ,
     +            R2D = 57.2957795130823208767981548147)
      ASIND = R2D * ASIN (RNUM)
      RETURN
      END
 
      FUNCTION ATAN2D (RNUM1,RNUM2)
      PARAMETER ( D2R =  0.0174532925199432957692369076847 ,
     +            R2D = 57.2957795130823208767981548147)
      ATAN2D = R2D * ATAN2 (RNUM1,RNUM2)
      RETURN
      END
CNCAR
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCCCCCCC          THE END OF Weimer96 model          CCCCCCCCCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
