!--------------------------------------------------------------------------      
      subroutine xyz_glb

!  XM = northward magnetic perturbation, nT
!  YM = eastward magnetic perturbation, nT
!  ZM = downward magnetic perturbation, nT
     
      use coef_module,only: fcmp,tmpo,nmp,mmp,nm,mm,
     |   c_mptextgrd,  ! magnetic potential external part @ ground
     |   c_mptintgrd,  ! magnetic potential induced part @ ground
     |   p,dp,
     |   r,pmopmmo,	! common/eqcr/
     |   f		! common/far/
      use param_module,only: rad,rmlon,rmlat,igeoplot
      use eqvcur_module,only: 
     |   reglat,    ! mag. dipole latitude from south to north pole
     |   reglon,    ! mag. dipole longitude
     |   regmlt,    ! reg grid mag. local time
     |   dB_phi,    ! magnetic perturbation deltaB_phi
     |   dB_lam,    ! magnetic perturbation deltaB_lam
     |   dB_dwn ,   ! magnetic perturbation deltaB_downward
     |   dB_phi_geo,    ! magnetic perturbation deltaB_phi (geo.direction)
     |   dB_lam_geo,    ! magnetic perturbation deltaB_lam (geo.direction)
     |   dB_dwn_geo     ! magnetic perturbation deltaB_downward (geo.direction)
      
      implicit none

      integer :: jmx,j,l,mp,m,n,np,mc,mq
      real :: th,ct,st,sts,ph,post(0:nm,0:mm),pm2,xm,ym,zm,cz,
     |   xmr,ymr,xmg,ymg,c_tot,fr,f1(2),f2(2)
     
      real :: xmlon
      common/ssolar/ xmlon
!      
      dB_phi = 0.
      dB_lam = 0.
      dB_dwn = 0.
      dB_phi_geo = 0.
      dB_lam_geo = 0.
      dB_dwn_geo = 0.      
!            
      DO 100 J=1,rmlat
        TH = (90. - reglat(J))/RAD
        CT = COS(TH)   
        ST = SIN(TH)                                                              
C Avoid poles, since there is division by ST
        ST = AMAX1(ST,1.E-7)
        STS = ST*ST                                                             
        DO 100 L=1,rmlon               
!         PH = regmlt(l)*15./RAD   ! use magn. local time          
         PH = (reglon(l)+180.-xmlon)/RAD        ! use magn. longitude
         CALL FCMP(PH,MM)                                                        
         P(0,0) = 1.
         POST(0,0) = 1./ST
         DO 80 MP=1,MMP    
           M = MP - 1                         
           IF (MP.NE.1) THEN 
	     POST(M,M) = PMOPMMO(M)*P(M-1,M-1)
	     P(M,M)    = POST(M,M)*ST 
           ENDIF
           DP(M,M) = float(M)*CT*POST(M,M)
           PM2 = 0.   
           DO 80 N=MP,NM             
             NP = N + 1                       
             POST(N,M) = (CT*POST(N-1,M) - R(N-1,M)*PM2)/R(N,M)
             P(N,M) = POST(N,M)*ST
             PM2 = POST(N-1,M)
             DP(N,M) = float(N)*CT*POST(N,M)-(2.*float(N)+1.)*
     |               R(N,M)*PM2
   80    continue
         XM = 0.
         YM = 0.
         ZM = 0.
         DO 120 m = -mm,mm		!Loop over order M  
           MP = 1 + IABS(M)                                                   
           DO 120 NP=MP,NMP 		!Loop over degree N   
             N = NP - 1 
	     c_tot =  c_mptextgrd(n,m) + c_mptintgrd(n,m)   ! v_n^m(ext) + v_n^m(int)      
             cz = float(n)*c_mptextgrd(n,m) - float(N+1)*   ! n  v_n^m(ext) - (n-1)v_n^m(int)    
     |          c_mptintgrd(n,m)
             XM = XM + C_TOT*DP(N,MP-1)*F(M)	! v_n^m(tot) d P_n^m/ d th * f(m)
             YM = YM - C_TOT*float(M)*POST(N,MP-1)*F(-m) ! - v_n^m(tot) m P_n^m* f(-m)
             ZM = ZM + CZ*P(N,MP-1)*F(M)	! cz P_n^m* f(m)
  120    CONTINUE 
  
! calculate magnetic perturbation delta  (in magnetic direction)
	 
         dB_phi(l,j) = dB_phi(l,j) + YM*1.E9
         dB_lam(l,j) = dB_lam(l,j) + XM*1.E9
         dB_dwn(l,j) = ZM*1.E9
	 
! calculate magnetic perturbation delta B_a (in geographic direction)
! dB = dB_qphi * vec(f2) x vec(k) + dB_qlam * vec(k) x vec(f1)
! dB_qphi = ym	dB_qlam = xm	
         if (igeoplot == 1) then
	    call cal_f1f2(l,j,fr,f1,f2)
            xmg =  f1(1)*XM-f2(1)*YM
            ymg = -f1(2)*XM+f2(2)*YM
                                                            
C Convert from T to nT 
	  dB_phi_geo(l,j) = dB_phi_geo(l,j) + YMG*1.E9
	  dB_lam_geo(l,j) = dB_lam_geo(l,j) + XMG*1.E9
          dB_dwn_geo(l,j) = ZM*sqrt(fr)*1.E9
	 endif

  100 CONTINUE 
           
      return
      end subroutine xyz_glb                                                                     
!---------------------------------------------------------------------------
      subroutine cal_f1f2(l,j,f,f1,f2)

!  am 2003-1-31
!   
!    calculate the base vector F and its components
!    uses the file "/tmp/maute/apex.d/apxntrp.grid.alt90.f90.1995" for f90
!    uses the file "/tmp/maute/apex.d/apxntrp.grid.alt90.f77.1995" for f77
!    files generated with program myglob.f / myglob.F

      use param_module,only: sglon,sglat,  ! geo.longitude & latitude
     |    fileapex,	! apex filename
     |    lwk,wk,apx_msgun,apx_iun, 
     |    hr, 	        ! reference height,altitude [km] used in TIEGCM
     |    alt    	! altitude [km] for apex
      use eqvcur_module,only: 
     |   bhat       ! unit vector along geomagnetic field direction
                    ! east north up
      
      implicit none

      integer,intent(in) :: l,j
      real,intent(out) :: f1(2),f2(2)

      integer :: ist
      real :: b(3),si,xlatm,vmp,w,d,be3,sim,f,
     |   d1(3),d2(3),d3(3),e1(3),e2(3),e3(3),
     |   bnrth,beast,bdown,babs,xlatqd,xlonm
      real :: xmlon,date
      common/ssolar/ xmlon,date
!
      if(l ==1.and.j==1) then
!        CALL COFRM (date)  ! not necessary since read from file
        CALL APXRDA (apx_msgun,fileapex,apx_iun,date,wk,lwk,ist)
      endif
!
!  get modified-apex and quasi-dipole coordinates and 
!  associated base vectors:
!  input  : GLAT,GLO ,ALT,HR, WK
!  output : B,BHAT,BABS,SI, 			       !Mag Fld
!           XLONM,            			       !Apx Lon
!           XLATM,VMP,W,D,BE3,SIM,D1,D2,D3,E1,E2,E3,   !Mod Apx
!           XLATQD,F,F1,F2 , IST 		       !Qsi-Dpl
      
      CALL APXMALL (SGLAT,SGLON,ALT,HR, WK,
     +             B,BHAT(l,j,:),BABS,SI,XLONM,
     +             XLATM,VMP,W,D,BE3,SIM,D1,D2,D3,E1,E2,E3,
     +             XLATQD,F,F1,F2 , IST)

      return
      end subroutine cal_f1f2
!----------------------------------------------------------------------------
      subroutine xyz_stat
C Input:
!   sglon             geograpgic longitude of point
!   sglat             geograpgic latitude of point
C calculated:
!   xx                magnetic local time of point
!   yy                latitude[degrees] from south to north -90 to 90
!   smlon             magnetic longitude of point
C  XM = northward magnetic perturbation, nT
C  YM = eastward magnetic perturbation, nT
C  ZM = downward magnetic perturbation, nT
     
      use coef_module,only: fcmp,tmpo,nmp,mmp,nm,mm,
     |   c_mptextgrd,  ! magnetic potential external part @ ground
     |   c_mptintgrd,  ! magnetic potential induced part @ ground
     |   p,dp,
     |   r,pmopmmo,	! common/eqcr/
     |   f		! common/far/
      use param_module,only: rad,sglon,sglat,fileapex,mtime
      
      implicit none

      integer :: mp,m,n,np,mc,mq
      real :: th,ct,st,sts,ph,post(0:nm,0:mm),pm2,xm,ym,zm,cz,c_tot,
     |   h,d
      real :: f1(2),f2(2),fr
      real :: xx,yy,smlon,smlat,d0
      real :: xmlon, date, colat,elon
      common/ssolar/ xmlon,date,colat,elon
!                                                                   
! calculate the f-base vectors therefore geo-coordinates needed
! input: sglon,sglat geo.lon&mag.lat
! output: f1,f2,f
! smlon,smlat:  mag.longitude, quasi dipole latitude [deg]

      call fvecstat(f1,f2,fr,smlon,smlat,d0)
      yy = smlat
      xx = (smlon-xmlon)*12./180.+12.    ! convert from mag.long to mag.loc.time
      if (xx .GE. 24.0) xx = xx - 24.0
      if (xx .LT.   0.) xx = xx + 24.0
!      write(55,*) xx,smlon,xmlon,mtime(2:3)

! yy: latitude southpole -90 to north pole +90 degree 
! th: colatitude southpole 180 to north pole 0 degree     
      TH = (90. - YY)/RAD
      CT = COS(TH)   
      ST = SIN(TH)                                                              
C Avoid poles, since there is division by ST
      ST = AMAX1(ST,1.E-7)
      STS = ST*ST 
! xx: magnetic local time not used here             
      PH = XX*15./RAD  
      CALL FCMP(PH,MMP)                                                        
      P(0,0) = 1.
      POST(0,0) = 1./ST
!   
      DO 80 MP=1,MMP	
        M = MP - 1			   
        IF (MP.NE.1) THEN 
          POST(M,M) = PMOPMMO(M)*P(M-1,M-1)
          P(M,M)    = POST(M,M)*ST 
        ENDIF
        DP(M,M) = float(M)*CT*POST(M,M)
        PM2 = 0.   
        DO 80 N=MP,NM		  
          NP = N + 1			   
          POST(N,M) = (CT*POST(N-1,M) - R(N-1,M)*PM2)/R(N,M)
          P(N,M) = POST(N,M)*ST
          PM2 = POST(N-1,M)
          DP(N,M) = float(N)*CT*POST(N,M)-(2.*float(N)+1.)*
     |  	  R(N,M)*PM2
   80 continue
!   
      XM = 0.
      YM = 0.
      ZM = 0. 
!   
      DO 120 m = -mm,mm 	     !Loop over order M  
        MP = 1 + IABS(M)						   
        DO 120 NP=MP,NMP	     !Loop over degree N   
          N = NP - 1 
	  c_tot =  c_mptextgrd(n,m) + c_mptintgrd(n,m)        
	  cz = float(n)*c_mptextgrd(n,m) - float(N+1)*
     |       c_mptintgrd(n,m)
!
!          c_tot =  c_mptextgrd(n,m)    
!          cz = float(n)*c_mptextgrd(n,m) 
          XM = XM + C_TOT*DP(N,MP-1)*F(M)
          YM = YM - C_TOT*float(M)*POST(N,MP-1)*F(-m)
          ZM = ZM + CZ*P(N,MP-1)*F(M)
  120 CONTINUE  					 
!
      XM = XM*1.E9   ! convert from T to nT 
      YM = YM*1.E9
      ZM = ZM*1.E9
!     
!      write(10,'(5(a,x,f10.5,x))') 'mlt',xx,'mlat',yy,'bp= ',xm,
!     |                           'bl= ',ym,'bz= ',zm

! take the base vectors into account
      XM =  XM*f1(1) - YM*f2(1)
      YM = -XM*f1(2) + YM*f2(2)
      ZM =  ZM*sqrt(fr)
!
      write(10,'(5(a,x,f10.5,x))') 'mlt',xx,'mlat',yy,'bx= ',xm,
     |                           'by= ',ym,'bz= ',zm
! convert to H and D
      h = xm*cos(d0) + ym*sin(d0)
      d = ym*cos(d0) - xm*sin(d0)
      
      write(10,'(5(a,x,f10.5,x))') 'mlt',xx,'mlat',yy,'h=  ',h,
     |                           'd=  ',d,'bz= ',zm
!     
      return                                                                  
      end subroutine xyz_stat                                                                      
!---------------------------------------------------------------------------
      subroutine fvecstat(f1,f2,f,xlonm,xlatqd,dec)

!  am 2001-9-6
!   
!    calculate the base vector F and its components
!    uses the file "/tmp/maute/apex.d/apxntrp.grid.alt90.f90.1995" for f90
!    uses the file "/tmp/maute/apex.d/apxntrp.grid.alt90.f77.1995" for f77
!    files generated with program myglob.f / myglob.F

      use param_module,only: sglon,sglat,  ! geo.longitude & latitude
     |    fileapex,	! apex filename
     |    lwk,wk,apx_msgun,apx_iun, 
     |    hr, 	        ! reference height,altitude [km] used in TIEGCM
     |    alt    	! altitude [km] for apex
      
      implicit none

      real,intent(out) :: f1(2),f2(2),f,xlonm,xlatqd,dec

      integer :: ist
      real :: b(3),bhat(3),si,xlatm,vmp,w,d,be3,sim,
     |   d1(3),d2(3),d3(3),e1(3),e2(3),e3(3),
     |   bnrth,beast,bdown,babs
      real :: xmlon,date
      common/ssolar/ xmlon,date
!
!      CALL COFRM (date) ! not necessary since read from file
      CALL APXRDA (apx_msgun,fileapex,apx_iun,date,wk,lwk,ist)
!
!  get modified-apex and quasi-dipole coordinates and 
!  associated base vectors:
!  input  : GLAT,GLO ,ALT,HR, WK
!  output : B,BHAT,BABS,SI, 			       !Mag Fld
!           XLONM,            			       !Apx Lon
!           XLATM,VMP,W,D,BE3,SIM,D1,D2,D3,E1,E2,E3,   !Mod Apx
!           XLATQD,F,F1,F2 , IST 		       !Qsi-Dpl
!        
      CALL APXMALL (SGLAT,SGLON,ALT,HR, WK,
     +             B,BHAT,BABS,SI,XLONM,
     +             XLATM,VMP,W,D,BE3,SIM,D1,D2,D3,E1,E2,E3,
     +             XLATQD,F,F1,F2 , IST)
!      write(55,*) 'geog. ',SGLON,SGLAT,XLONM,XLATQD,XLATM

!  Compute the DGRF/IGRF field components at the point GLAT,GLON,ALT.
!  COFRM must be called to establish coefficients for correct date
!  prior to calling FELDG.
!  IENTY = 1
!    input:
!      GLAT = Latitude of point (deg)
!      GLON = Longitude (east=+) of point (deg)
!      ALT  = Ht of point (km)
!    output:
!      BNRTH  north component of field vector (Gauss)
!      BEAST  east component of field vector  (Gauss)
!      BDOWN  downward component of field vector (Gauss)
!      BABS   magnitude of field vector (Gauss)
      call FELDG (1,SGLAT,SGLON,ALT, BNRTH,BEAST,BDOWN,BABS)

! calculate declination angle
      dec = atan2(beast,bnrth)
!      write(6,*) 'dec= ',dec/180.*3.1415926535898 ,BNRTH,BEAST

      return
      end subroutine fvecstat
!--------------------------------------------------------------------------     
      subroutine cal_bmagnitude

      use param_module,only: rad,rmlon,rmlat
      use eqvcur_module,only:
     |   vallat,    ! TIEGCM mag. quasi-dipole latitude from south to north pole      
     |   dB_phi,    ! magnetic perturbation deltaB_phi
     |   dB_lam,    ! magnetic perturbation deltaB_lam
     |   dB_dwn ,   ! magnetic perturbation deltaB_downward
     |   dB_phi_geo,! magnetic perturbation deltaB_phi (geo.direction)
     |   dB_lam_geo,! magnetic perturbation deltaB_lam (geo.direction)
     |   dB_dwn_geo,! magnetic perturbation deltaB_downward (geo.direction)
     |   bhat,      ! unit vector along geomagnetic field direction
                    ! east north up
     |   babs_mag,   ! total B_0 + dB
     |   babs_geo   ! total B_0 + dB
      
      implicit none

      integer :: i,j,k
      real :: coslm,sinlm,fac,cosIq,sinIq

! loop over all quasi-dipole coordinates
      
      do j = 1,rmlat
        coslm = cos(vallat(j)/rad)
        sinlm = sin(vallat(j)/rad)
	fac = sqrt(4.-3.*coslm**2)
	fac = 1/fac
        cosIq = coslm*fac
	sinIq = 2.*sinlm*fac
	
        do i = 1,rmlon 
	   babs_geo(i,j) =  bhat(i,j,1)*dB_phi_geo(i,j) +
     |       bhat(i,j,2)*dB_lam_geo(i,j) -  bhat(i,j,3)*dB_dwn_geo(i,j)
	   babs_mag(i,j) = cosIq*dB_lam(i,j) + sinIq*dB_dwn(i,j)
        enddo
      enddo
         
      return
      end subroutine cal_bmagnitude
!-------------------------------------------------------------------------- 
