!-----------------------------------------------------------------------
      subroutine init_geo2mag
!
! Determine points to be passed between tasks for parallel geo to mag 
! transformation.
!
! Module data from mpi.F:
!     integer,allocatable,dimension(:,:,:) :: ! (mlon0:mlon1,mlat0:mlat1,4)
!    |  igeo,  ! geo lon indices at 4 pts boxing each mag pt in my subdomain
!    |  jgeo,  ! geo lat indices at 4 pts boxing each mag pt in my subdomain
!    |  idgeo  ! tid of mpi task with each of the 4 pts (owner)
!     integer,allocatable :: nrcvgeo(:) ! # geo pts to be recvd from each task
!     integer,allocatable :: nsndgeo(:) ! # geo pts to be sent to each ntask
!     integer,allocatable,dimension(:,:) :: ! (mxrcv[mxsnd],ntask) 
!    |  ircvgeo, ! lon indices to be received from each task
!    |  jrcvgeo, ! lat indices to be received from each task
!    |  isndgeo, ! lon indices to be sent from each task
!    |  jsndgeo  ! lat indices to be sent from each task
!
      use mpi_module,only: ntask,mlon0,mlon1,mlat0,mlat1,mymlons,
     |  mymlats,igeo,jgeo,idgeo,ircvgeo,jrcvgeo,isndgeo,jsndgeo,
     |  nsndgeo,nrcvgeo,tasks,ntaskj,mytid,lat0,lat1,lon0,lon1
      use init_module,only: glat,glon,gmlat,gmlon
!
! Local:
      integer :: imag,jmag,i,ii,j,n,nn,nmlats,nmlons,istat,nsend,lonend
      integer :: mxim,mxjm,mxig,mxjg
      integer :: mxsndgeo,mxrcvgeo
      logical,allocatable,dimension(:,:,:) :: 
     |  counted ! (nlonp4,nlat,ntask)
!
! Set magnetic grid in degrees:
!
! Allocate arrays (module data from mpi.F):
      allocate(igeo(mlon0:mlon1,mlat0:mlat1,4),stat=istat)
      if (istat /= 0) write(6,"('>>> init_geo2mag: error allocating',
     |  ' igeo: stat=',i3)") istat
      allocate(jgeo(mlon0:mlon1,mlat0:mlat1,4),stat=istat)
      if (istat /= 0) write(6,"('>>> init_geo2mag: error allocating',
     |  ' jgeo: stat=',i3)") istat
      allocate(idgeo(mlon0:mlon1,mlat0:mlat1,4),stat=istat)
      if (istat /= 0) write(6,"('>>> init_geo2mag: error allocating',
     |  ' idgeo: stat=',i3)") istat
      allocate(nrcvgeo(0:ntask-1),stat=istat)
      if (istat /= 0) write(6,"('>>> init_geo2mag: error allocating',
     |  ' nrcvgeo: stat=',i3)") istat
      allocate(nsndgeo(0:ntask-1),stat=istat)
      if (istat /= 0) write(6,"('>>> init_geo2mag: error allocating',
     |  ' nsndgeo: stat=',i3)") istat
!
! Assign igeo and jgeo from ig,jg (which were read by magfield):
! Numbering convention for the 4 geo pts surrounding a single mag pt
!   is 1 to 4, counterclockwise starting with 1 in lower left. 
!
      do jmag=mlat0,mlat1
        do imag=mlon0,mlon1
          igeo(imag,jmag,1) = ig(imag,jmag)  
          igeo(imag,jmag,2) = ig(imag,jmag)+1
          igeo(imag,jmag,3) = ig(imag,jmag)+1
          igeo(imag,jmag,4) = ig(imag,jmag)
!
          jgeo(imag,jmag,1) = jg(imag,jmag)  
          jgeo(imag,jmag,2) = jg(imag,jmag)
          jgeo(imag,jmag,3) = jg(imag,jmag)+1
          jgeo(imag,jmag,4) = jg(imag,jmag)+1
        enddo ! imag=mlon0,mlon1
      enddo ! jmag=mlat0,mlat1
!
! Out of curiosity:
      mxim = maxval(im) ; mxjm = maxval(jm) 
      mxig = maxval(ig) ; mxjg = maxval(jg)
      write(6,"('init_geo: mxim=',i3,' mxjm=',i3,' mxig=',i3,
     |  ' mxjg=',i3)") mxim,mxjm,mxig,mxjg
!
      nmlats = mlat1-mlat0+1
      nmlons = mlon1-mlon0+1
!
! Determine owners (tasks) of each mag grid pt in my subdomain:
      idgeo(:,:,:) = -1
      do jmag=mlat0,mlat1
        do imag=mlon0,mlon1
          do i=1,4
            do n=0,ntask-1
              if (any(tasks(n)%mylons==igeo(imag,jmag,i)).and.
     |            any(tasks(n)%mylats==jgeo(imag,jmag,i))) then
                idgeo(imag,jmag,i) = n ! this needed pt is owned by task n
              elseif (any(tasks(n)%mylons==igeo(imag,jmag,i)).and.
     |          jgeo(imag,jmag,i) > nlat .and.
     |          tasks(n)%mytidj==ntaskj-1) then ! northernmost tasks
                idgeo(imag,jmag,i) = n ! owner of nlat+1,nlat+2
              endif
            enddo ! n=0,ntask-1
            if (idgeo(imag,jmag,i)==-1) then
              write(6,"('>>> init_geo2mag note: no owner: ',
     |          ' jmag=',i3,' imag=',i3,' i=',i3,' igeo=',i3,' jgeo=',
     |          i3)") jmag,imag,i,igeo(imag,jmag,i),jgeo(imag,jmag,i)
              stop 'init_geo2mag'
            endif
          enddo ! i=1,4
        enddo ! imag=mlon0,mlon1
      enddo ! jmag=mlat0,mlat1
      if (any(idgeo==-1)) then
        write(6,"('>>> WARNING init_geo2mag: some idgeo == -1')")
        stop 'idgeo'
      endif
!
! Determine how many points I must receive from each task, for
!   allocation purposes below.
! Use logical counted to remove redundancy -- only need to send/recv 
!   each point once.
!
      allocate(counted(nlonp4,nlat,0:ntask-1),stat=istat)
      if (istat /= 0) write(6,"('>>> init_geo2mag: error allocating',
     |  ' counted: stat=',i3)") istat
      counted(:,:,:) = .false. ! (nlonp4,nlat,ntask)
      nrcvgeo(:) = 0
      do jmag=mlat0,mlat1
        do imag=mlon0,mlon1
          do i=1,4
            n = idgeo(imag,jmag,i)
            if (.not.counted(igeo(imag,jmag,i),jgeo(imag,jmag,i),n).and.
     |        n /= mytid) then
              nrcvgeo(n) = nrcvgeo(n)+1
              counted(igeo(imag,jmag,i),jgeo(imag,jmag,i),n) = .true.
            endif
          enddo ! i=1,4
        enddo ! imag=mlon0,mlon1
      enddo ! jmag=mlat0,mlat1
      do n=0,ntask-1
        write(6,"('init_geo2mag: n=',i3,' nrcvgeo=',i4)") n,nrcvgeo(n)
      enddo
!
! Allocate ircvgeo,jrcvgeo:
      mxrcvgeo = 0
      do n=0,ntask-1
        if (nrcvgeo(n) > mxrcvgeo) mxrcvgeo = nrcvgeo(n)
      enddo
      allocate(ircvgeo(mxrcvgeo,0:ntask-1),stat=istat)
      if (istat /= 0) write(6,"('>>> init_geo2mag: error allocating',
     |  ' ircvgeo: stat=',i3)") istat
      allocate(jrcvgeo(mxrcvgeo,0:ntask-1),stat=istat)
      if (istat /= 0) write(6,"('>>> init_geo2mag: error allocating',
     |  ' jrcvgeo: stat=',i3)") istat
!
! Determine lon,lat indices to be received from each task.
! 
      counted(:,:,:) = .false. ! (nlonp4,nlat,ntask)
      nrcvgeo(:) = 0

      write(6,"('sndtid',2x,'nrcv',2x,'  geo lon',8x,
     |  '  geo lat',8x,'  mag lon',2x,'  mag lat')") 

      do jmag=mlat0,mlat1
        do imag=mlon0,mlon1
          do i=1,4
            n = idgeo(imag,jmag,i) ! owner of current geo pt
            if (.not.counted(igeo(imag,jmag,i),jgeo(imag,jmag,i),n).and.
     |        n /= mytid) then
              nrcvgeo(n) = nrcvgeo(n)+1
              ircvgeo(nrcvgeo(n),n) = igeo(imag,jmag,i)
              jrcvgeo(nrcvgeo(n),n) = jgeo(imag,jmag,i)
              counted(igeo(imag,jmag,i),jgeo(imag,jmag,i),n) = .true.

              write(6,"(i6,2x,i4,2x,f9.2,' i=',i3,2x,f9.2,' j=',i3,2x,
     |          f9.2,2x,f9.2)") n,nrcvgeo(n),glon(ircvgeo(nrcvgeo(n),n))
     |          ,igeo(imag,jmag,i),glat(jrcvgeo(nrcvgeo(n),n)),
     |          jgeo(imag,jmag,i),gmlon(imag),gmlat(jmag)

            endif
          enddo ! i=1,4
        enddo ! imag=mlon0,mlon1
      enddo ! jmag=mlat0,mlat1
!
! Now determine how many of my points need to be sent to each 
! of the other tasks, so isndgeo and ircvgeo can be allocated.
!
      lonend = lon1
      if (lon1==nlonp4) lonend = nlon
      counted(:,:,:) = .false. ! (nlonp4,nlat,ntask)
      nsndgeo(:) = 0
      do j=lat0,lat1
        do i=lon0,lonend
          do n=0,ntask-1
            if (n /= mytid) then
              if (any(tasks(n)%mymlats==jm(i,j)).and.
     |            any(tasks(n)%mymlons==im(i,j)).and.
     |            .not.counted(i,j,n)) then
                nsndgeo(n) = nsndgeo(n)+1
                counted(i,j,n) = .true.
              endif
              if (any(tasks(n)%mymlats==jm(i+1,j)).and.
     |            any(tasks(n)%mymlons==im(i+1,j)).and.
     |            .not.counted(i+1,j,n)) then
                nsndgeo(n) = nsndgeo(n)+1
                counted(i+1,j,n) = .true.
              endif
              if (any(tasks(n)%mymlats==jm(i+1,j+1)).and.
     |            any(tasks(n)%mymlons==im(i+1,j+1)).and.
     |            .not.counted(i+1,j+1,n)) then
                nsndgeo(n) = nsndgeo(n)+1
                counted(i+1,j+1,n) = .true.
              endif
              if (any(tasks(n)%mymlats==jm(i,j+1)).and.
     |             any(tasks(n)%mymlons==im(i,j+1)).and.
     |            .not.counted(i,j+1,n)) then
                nsndgeo(n) = nsndgeo(n)+1
                counted(i,j+1,n) = .true.
              endif
            endif
          enddo ! n=0,ntask-1
        enddo ! i=lon0,lon1
      enddo ! j=lat0,lat1
      write(6,"(' ')")
      do n=0,ntask-1
        write(6,"('init_geo2mag: n=',i3,' nsndgeo=',i4)") n,nsndgeo(n)
      enddo ! n=0,ntask-1
!
! Allocate isndgeo and jsndgeo:
      mxsndgeo = 0
      do n=0,ntask-1
        if (nsndgeo(n) > mxsndgeo) mxsndgeo = nsndgeo(n)
      enddo
      allocate(isndgeo(mxsndgeo,0:ntask-1),stat=istat)
      if (istat /= 0) write(6,"('>>> init_geo2mag: error allocating',
     |  ' isndgeo: stat=',i3)") istat
      allocate(jsndgeo(mxsndgeo,0:ntask-1),stat=istat)
      if (istat /= 0) write(6,"('>>> init_geo2mag: error allocating',
     |  ' jsndgeo: stat=',i3)") istat
!
! Set isndgeo and jsndgeo, i.e., lat,lon indices of my geo points
!   that need to be sent to each of the other tasks (also reset nsndgeo).
!

      write(6,"('rcvtid',2x,'nsnd',2x,'  geo lon',8x,
     |  '  geo lat',8x,'  mag lon',2x,'  mag lat')") 

      lonend = lon1
      if (lon1==nlonp4) lonend = nlon
      counted(:,:,:) = .false. ! (nlonp4,nlat,ntask)
      nsndgeo(:) = 0
      isndgeo(:,:) = 0
      jsndgeo(:,:) = 0
      do j=lat0,lat1
        do i=lon0,lonend
          do n=0,ntask-1
            if (n /= mytid) then
! Geo point 1 (lower left of mag pt):
              if (any(tasks(n)%mymlats==jm(i,j)).and.
     |            any(tasks(n)%mymlons==im(i,j)).and.
     |            .not.counted(i,j,n)) then
                nsndgeo(n) = nsndgeo(n)+1
                isndgeo(nsndgeo(n),n) = i
                jsndgeo(nsndgeo(n),n) = j
                counted(i,j,n) = .true.

                write(6,"(i6,2x,i4,2x,f9.2,' i=',i3,2x,f9.2,' j=',i3,2x,
     |            f9.2,2x,f9.2)") n,nsndgeo(n),
     |            glon(isndgeo(nsndgeo(n),n)),i,
     |            glat(jsndgeo(nsndgeo(n),n)),j,gmlon(im(i,j)),
     |            gmlat(jm(i,j))

              endif

! Geo point 2 (lower right of mag pt):
              if (any(tasks(n)%mymlats==jm(i+1,j)).and.
     |            any(tasks(n)%mymlons==im(i+1,j)).and.
     |            .not.counted(i+1,j,n)) then
                nsndgeo(n) = nsndgeo(n)+1
                isndgeo(nsndgeo(n),n) = i+1
                jsndgeo(nsndgeo(n),n) = j
                counted(i+1,j,n) = .true.

                write(6,"(i6,2x,i4,2x,f9.2,' i=',i3,2x,f9.2,' j=',i3,2x,
     |            f9.2,2x,f9.2)") n,nsndgeo(n),
     |            glon(isndgeo(nsndgeo(n),n)),i,
     |            glat(jsndgeo(nsndgeo(n),n)),j,gmlon(im(i+1,j)),
     |            gmlat(jm(i+1,j))

              endif

! Geo point 3 (upper right of mag pt):
              if (any(tasks(n)%mymlats==jm(i+1,j+1)).and.
     |            any(tasks(n)%mymlons==im(i+1,j+1)).and.
     |            .not.counted(i+1,j+1,n)) then
                nsndgeo(n) = nsndgeo(n)+1
                isndgeo(nsndgeo(n),n) = i+1
                jsndgeo(nsndgeo(n),n) = j+1
                counted(i+1,j+1,n) = .true.

                write(6,"(i6,2x,i4,2x,f9.2,' i=',i3,2x,f9.2,' j=',i3,2x,
     |            f9.2,2x,f9.2)") n,nsndgeo(n),
     |            glon(isndgeo(nsndgeo(n),n)),i,
     |            glat(jsndgeo(nsndgeo(n),n)),j,gmlon(im(i+1,j+1)),
     |            gmlat(jm(i+1,j+1))

              endif

! Geo point 4 (upper left of mag pt):
              if (any(tasks(n)%mymlats==jm(i,j+1)).and.
     |            any(tasks(n)%mymlons==im(i,j+1)).and.
     |            .not.counted(i,j+1,n)) then
                nsndgeo(n) = nsndgeo(n)+1
                isndgeo(nsndgeo(n),n) = i
                jsndgeo(nsndgeo(n),n) = j+1
                counted(i,j+1,n) = .true.

                write(6,"(i6,2x,i4,2x,f9.2,' i=',i3,2x,f9.2,' j=',i3,2x,
     |            f9.2,2x,f9.2)") n,nsndgeo(n),
     |            glon(isndgeo(nsndgeo(n),n)),i,
     |            glat(jsndgeo(nsndgeo(n),n)),j,gmlon(im(i,j+1)),
     |            gmlat(jm(i,j+1))

              endif
            endif
          enddo ! n=0,ntask-1
        enddo ! i=lon0,lon1
      enddo ! j=lat0,lat1
!
      deallocate(counted)
      end subroutine init_geo2mag
!-----------------------------------------------------------------------
      end module magfield_module

