c
c
c

	program main


	parameter (Re = 6.49E6 )                ! 地球半径 (m)+(130km)
	parameter (PII = 3.141592 )             ! 円周率
	real*4  EPHIY (2)
	real*4  EPHIZ (2)

	INTEGER           EGNR,AGNR,OGNR
        COMMON/GENER/     UMR,ERA,AQUAD,BQUAD


      open(99,file='hmcoef.dat',status='old')


	RAD = 180./PII
	NMLT= 24
	DLMLT = 2.*PII/FLOAT(NMLT)
	DLILAT= FLOAT(1)/RAD
	DLY2D = 2.*DLMLT*Re
	DLZ2M = 2.*DLILAT*Re


C------------------------------------------------------------
C   MODLE  - INDEX FOR SPECIFYING POTENTIAL MODEL                      *
C            1 = HEPPNER-MAYNARD MODEL A                               *
C            2 = HEPPNER-MAYNARD MODEL BC                              *
C            3 = HEPPNER-MAYNARD MODEL DE                              *
C            4 = HEPPNER-MAYNARD MODEL BCP                             *
C            5 = HEPPNER-MAYNARD MODEL BCPP                            *
C            6 = HEPPNER-MAYNARD MODEL DEP                             *
C            7 = HEPPNER-MAYNARD MODEL DEPP                            *
C            8 = HEELIS MODEL                                          *

	MODLE= 1

	! ICHGHM
	! 0:No Change
	! 1:Change The HM Pattern.
	ICHGHM = 0 

	!
	! IMLT : Magnetic Local Time (0-24)
	!

	MLT = 18
	MLT = 18
	MLT = 6



	TLOC = MLT 
	GMLONG = TLOC*15.



	!
	! GMLAT > 60
	! 
        GMLAT =  65.4	! (deg)
        GMLAT =  70.0	! (deg)

	GMLAT = 60.0

	dt=50.e0


	do 2000  iday =1,10
	do 1000  iutsec =1, int(86400./dt)


	UTsec = real(iutsec)

	 ! 2-1) 磁気座標, 地理座標, UTから MLT(0-24) を計算
	  call get_MLT(
     I                  0,UTsec,Gmlat,Gmlong,
     O                  dMLT )

       TLOC = dMLT  !  - 12.       !TLOC (-12 : 12)

c	TLOC = GMLONG/15.



!
! Calculate Electric Potential (V)
!
	!
	!1) 経度方向の電位を求める
	!
	do 10 ILOC=1,2	

	 DLOC  = (real(ILOC)-1.5)*2.	! ILOC =1 : DLOC=-1
	 XLONG =  (TLOC + DLOC)*15.	! ILOC =2 : DLOC= 1

	 if(XLONG .lt. 0.) XLONG= 360.+XLONG

	 XMLAT =   GMLAT

		!
		! Heppner-Maynard-Rich Electric Field Model
		!	サブルーチンの呼出
		!
         CALL EPOT(XMLAT,XLONG,PHI,99,MODLE,16,ACHG,BCHG,DXCHG,DYCHG,
     X            AHM,BHM,DXHM,DYHM,ICHGHM)
         IF( ICHGHM .EQ. 1 ) PHI=PHI*CHGV

         EPHIY(ILOC) = PHI*1.E3		! (kV) => (V) 

   10 	continue

	!
	!2) 緯度方向の電位を求める
	!
	do 20 IMLAT=1,2	

	 XLONG =  GMLONG
	 XLONG =  TLOC *15.	


	 DMLAT =(real(IMLAT)-1.5)*2.	! IMLAT =1 : DDMLAT=-1
					! IMLAT =2 : DDMLAT= 1
		!
		! 緯度 88 度以上の処理
		!
	 if(GMLAT.le.87.)then
	  XMLAT =   GMLAT + DMLAT
	 else
	  XMLAT =   87. + DMLAT
	 endif

		!
		! Heppner-Maynard-Rich Electric Field Model
		!	サブルーチンの呼出
		!
         CALL EPOT(XMLAT,XLONG,PHI,99,MODLE,16,ACHG,BCHG,DXCHG,DYCHG,
     X            AHM,BHM,DXHM,DYHM,ICHGHM)

         IF( ICHGHM .EQ. 1 ) PHI=PHI*CHGV
         EPHIZ(IMLAT) = PHI*1.E3	! (kV) => (V) 
   20 	continue

	!
	! Calculate Electric Field (mV/M)
	!3) 電位差を求める
	!

	COLAT = (90.-GMLAT)/RAD
	DLY2M =DLY2D*SIN(COLAT)
	
	EY= 1.E3 * (EPHIY(2) - EPHIY(1))/DLY2M ! MLT(高地方時-低地方時)(mV/M)
	EX= 1.E3 * (EPHIZ(2) - EPHIZ(1))/DLZ2M ! MLAT(高緯度-低緯度)(mV/M)


!
! IGRF の呼出
!
! X: 北向き ; Y:東向き ; Z:下向き
	
	!初期化
	call INITIZE	


	!
	call FELDCOF(1999.,DIMO)
	call FELDG(GMLAT,TLOC*15.,130.,BNORTH,BEAST,BDOWN,BABS)
	CALL SHELLG(LATI,LONGI,HEIGHT,DIMO,XL,ICODE,BAB1)

	DIP=ASIN(BDOWN/BABS)/UMR
	DEC=ASIN(BEAST/SQRT(BEAST*BEAST+BNORTH*BNORTH))/UMR
			   
	! write(6,*) 'B(N,E,D) ,Value',BNORTH,BEAST,BDOWN,BABS
	! write(6,*) 'Ve:',sqrt(EY**2+EX**2)/1.E3 /BABS *1.E4 ! m/s
	! メモ sqrt(BNORTH**2+BEAST**2+BDOWN**2) = BABS (GAUSS)


	B_2 = BNORTH**2 + BEAST**2 + BDOWN**2

!
! ドリフト速度を求める
!
	VX =  EY*BDOWN/B_2 *1.e1	! EY(mV/m)*Gauss/Gauss**2 *1.e1 =>(m/s)
	VY =  -EX*BDOWN/B_2 *1.e1


!
! 次ステップでの  LT, MLAT を決定
!
	dMLAT  = dt*VX/Re*180./pii
	dMLONG = dt*VY/Re/cos(GMLAT*PII/180.)*180./pii

	GMLAT = GMLAT + dMLAT
	GMLONG = GMLONG + dMLONG


c	TLOC = GMLONG/15.

 	RR=Re*cos(GMLAT/180.*PII)

	 XLONG =  TLOC *15.	
	x=RR*cos(XLONG/180.*PII) 
	y=RR*sin(XLONG/180.*PII) 

c	if(iutsec .eq. int(iutsec/100)*100)then
	 write(6,200) x,y, real(iutsec),GMLAT,GMLONG,vy,dt*vy,B_2
c	endif
  200	format(6(f13.3,' '), 4(e13.5,' '))


 1000	continue
 2000	continue



	stop
	end


c
c get_MLT  磁気地方時を与える
c
	subroutine get_MLT(
     O		ideg,UTsec,dmlat,dmlong,
     I		dmlt)

c	INPUT parameter
c		ideb  : 入出力(deg/rad)のスイッチ
c			ideg : 0      入力が deg
c			ideg : 1      入力が rad
c		UTsec : 世界標準時 (sec)
c		dmlat : 磁気緯度 (NP : EQ : SP),(90 :0: -90),(pi/2:0:-pi/2)
c		dmlong: 磁気経度 (-180:180),(-pi:pi)
c
c	OUTPUT parameter
c		dmlt  :  磁気地方時 (0:24)
c
	PARAMETER (pi = 3.1415926535d0)
	PARAMETER (dlambda0 = 78.8d0/90 * pi/2 )  !  78.8 deg N
	PARAMETER (dphi0    = -70.9d0/180 * pi)   ! 289.1 deg E => 70.9 W



	! 入力が deg ならば一度 rad に直す
	if(ideg .eq. 0) then
		amlat = dmlat/90. * pi/2.
                amlong= dmlong/180. * pi
	else
		amlat = dmlat
		amlong= dmlong
	endif


	dmlt_gp = UTsec/3600. + dphi0/pi*12. ! 磁極のMLT

	dmlt = dmlt_gp  + amlong/pi*12.	     ! 求める MLT

	

 	if (dmlt .lt. 0.)  dmlt = 24. + dmlt
	if (dmlt .gt. 24.) dmlt = dmlt -24.
	

	return
	end
