!
      module map_module
!      
      use param_module,only: rglon,rglat
      
      implicit none
      
      integer :: im(rglon,rglat),jm(rglon,rglat)
      real :: dim(rglon,rglat),djm(rglon,rglat)
!           
      contains
!---------------------------------------------------------------------------
! set up regular geographic grid and prepare weighting factors
! for mapping from the magnetic to the geographic grid
!
      subroutine prep_mag2geo
      
      use param_module,only: alt, 	! satellite height         
     |   hr,  		! reference height 90 km from TIEGCM [km]  
     |   fileapex,	! apex filename
     |   rmlon,rmlat,	! dimension of regular mag.grid
     |   rad,	        ! 180./pi
     |   apx_msgun,apx_iun,apx_mlat,apx_mlon,apx_malt,lwk
     
      use eqvcur_module,only: 
     |   reglat,    ! mag. dipole latitude from south to north pole
     |   reglon,    ! mag. dipole longitude 
     |   reg_glat,  ! reg. grid geographic latitude from south to north pole
     |   reg_glon   ! reg grid geographic longitude
      
      implicit none

      real :: xmlon,date,bnrth,beast,bdown,b0,
     |   f,xlonm,xlatm,xlatqd,f1_tmp(2),f2_tmp(2),
     |   gdlat,gdlon,f1(2),f2(2),bhat(3),d,be3,si,sim,
     |   w,vmp,colat,elon,sbsllat,sbsllon
      common/ssolar/ xmlon,date,colat,elon,sbsllat,sbsllon
      
      integer :: i,j,jjm,ist
      
      real :: alonm(rglon,rglat),alatm(rglon,rglat)
      real :: dlatm,dlonm,xlonmi
      
      real :: b(3),d1(3),d2(3),d3(3),e1(3),e2(3),e3(3)
      real :: wk(lwk)

      dlonm = 360./(rmlon-1)
      dlatm = 180./(rmlat-1)
      
      call cofrm(date)
      write(6,'(a5,f12.6)') 'date ',date
      write(6,*) fileapex
      call apxrda (apx_msgun,fileapex,apx_iun,date,wk,lwk,ist)

! loop over all geographic coordinates
      
      do i = 1,rglon
        do j = 1,rglat 
	     
      
!  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 (reg_glat(j),reg_glon(i),alt,hr, wk,
     +             b,bhat,b0,si,xlonm,
     +             xlatm,vmp,w,d,be3,sim,d1,d2,d3,e1,e2,e3,
     +             xlatqd,f,f1_tmp,f2_tmp,ist)
          
          if (ist /= 0) stop 'prep_mag2geo: apxmall'
	  
          alonm(i,j)    =  xlonm
          alatm(i,j)    =  xlatqd
! Set up parameters for magnetic to geographic interpolation 
	  xlonmi = (alonm(i,j) - reglon(1))/dlonm
     
	  if (xlonmi < 0.) xlonmi = xlonmi + float(rmlon)
	  im(i,j) = xlonmi
	  dim(i,j) = xlonmi - float(im(i,j))
	  im(i,j) = im(i,j) + 1
	  if (im(i,j) >= rmlon) im(i,j) = im(i,j) - float(rmlon)
	  alatm(i,j) = amin1(alatm(i,j),reglat(rmlat))
          do jjm=2,rmlat
            if (alatm(i,j) > reglat(jjm)) cycle
            jm(i,j) = jjm - 1
	    djm(i,j) = (alatm(i,j) - reglat(jm(i,j)))/
     |                 (reglat(jjm) - reglat(jm(i,j)))
            exit
          enddo
	  
	enddo   ! j/latitude-loop
      enddo	! i/longitude loop
      
      return
      end subroutine prep_mag2geo
!-----------------------------------------------------------------------
! mapping from magnetic to geographic grid
!
      subroutine mag2geo(am,ag,im,jm,dim,djm,lg,lm,nlong,nlatg,nlonm,
     |  nlatm)
!
! Args:
      integer,intent(in) :: lg,lm,nlong,nlatg,nlonm,nlatm
      integer,intent(in) :: im(lg,*),jm(lg,*)
      real,intent(in) :: am(lm,*),dim(lg,*),djm(lg,*)
      real,intent(out) :: ag(lg,*)
!
! Local:
      integer :: ig,jg
!
      do jg=1,nlatg 
        do ig=1,nlong
          ag(ig,jg) =
     |    am(im(ig,jg)  ,jm(ig,jg))  *(1.-dim(ig,jg))*(1.-djm(ig,jg))+
     |    am(im(ig,jg)+1,jm(ig,jg))  *    dim(ig,jg) *(1.-djm(ig,jg))+
     |    am(im(ig,jg)  ,jm(ig,jg)+1)*(1.-dim(ig,jg))*djm(ig,jg)+
     |    am(im(ig,jg)+1,jm(ig,jg)+1)*    dim(ig,jg) *djm(ig,jg)
        enddo ! ig=1,nlong
      enddo ! jg=1,nlatg 

      return
      end subroutine mag2geo
!-----------------------------------------------------------------------
      
      end module map_module