SUBROUTINE Move_plasma_grid (dt, ie_ask)
    IMPLICIT NONE
    REAL (rprec), INTENT (IN) :: dt
    INTEGER (iprec), INTENT (IN) :: ie_ask
!_____________________________________________________________________________
!   Subroutine to advance eta distribution for a time step
!   by a lifetime-based algorithm (raw doc dated 5/12/87)
!                                                                       
!   Last update: 05-11-88
!                01-29-96 ain ,min_j and calls to bndy added - frt
!   rws          06-05-97 etamov changed to reflect new use of
!                        eeta array in rcm697 version
!
!   CALLS:       BNDY, GNTRP, CIRCLE, RATEFN, BJMOD
!   CALLED FROM: TSTEP1
!_____________________________________________________________________________
!
!
    REAL (rprec) :: eeta2 (isize,jsize), veff  (isize,jsize), &
                    dvefdi(isize,jsize), dvefdj(isize,jsize), &
                    didt, djdt, biold, bjold, rate, mass_factor, r_dist
    INTEGER (iprec) :: i, j, kc, ie
!
!
    DO kc = 1, kcsize
!
       veff = v + vcorot - vpar + alamc(kc)*vm
       CALL V_eff_polar_cap (veff)
!
       ie = ikflavc(kc)
!
IF (ie /= ie_ask) CYCLE
!
       mass_factor = SQRT (xmass(1) / xmass(ie))
!
!
!      Compute spatial derivatives of effective potential for given energy:
!      (only inside the modeling region, boundaries possibly included).
!
       dvefdi = Deriv_i_of_Veff (veff)
       dvefdj = Deriv_j_of_Veff (veff)
!
!
!      3. Find the position from which particles of given energy have come
!         from to reach a given grid point (i,j) in one time step, and
!         interpolate EETA at that position. Then estimate that EETA at
!         (i,j) point is the same as that value, corrected for loss.
!
!         Notice: we are moving plasma only inside the modeling region:
!         (therefore, derivatives dvefdi and dvefdj are used only inside).
!
       DO j = j1, j2
          DO i = imin_j(j), i2
             fac (i,j) =  1.0E-3_rprec * bir(i,j) * alpha(i,j) * beta(i,j) *&
                          dlam * dpsi * ri**2 * signbe
             didt   =   dvefdj (i,j) / fac (i,j)
             djdt   = - dvefdi (i,j) / fac (i,j)
!            biold  = MAX (REAL(i,rprec) - didt * dt , Bndy (bndloc, REAL(j,rprec) ) )
             biold  = REAL(i,rprec) - didt*dt
             bjold  = Bjmod (REAL(j,rprec) - djdt * dt, jwrap, jsize )

             IF (ie == 1) THEN
                IF (kc ==1 ) THEN 
!                  Cold plasmaspheric electrons. Dont' do pitch angle scattering,
!                  but compute refilling rates for the plasmasphere
                   rate = - Plasmasphere_refill_rate (rmin(i,j), doy, &
                                        sunspot_number, eeta(i,j,kc), vm(i,j))
                   eeta(i,j,kc) = eeta(i,j,kc)+rate*dt
                   rate = 0.0
                 ELSE
!                   Plasmasheet electrons, pitch-angle scattering loss rate:
                    rate = Ratefn (fudgec(kc), alamc(kc), sini(i,j),&
                                           bir (i,j), vm(i,j), mass_factor)
                 END IF
             ELSE 
!              Positive ions, compute charge-exchange rate is it is on:
               IF (L_dktime .AND. i >= imin_j(j)) THEN
                r_dist = SQRT(xmin(i,j)**2+ymin(i,j)**2)
                rate = Cexrat (ie, ABS(alamc(kc))*vm(i,j), &
                                           R_dist, &
                                           sunspot_number, dktime, &
                                           irdk,inrgdk,isodk,iondk)
                ELSE
                     rate = 0.0
                END IF
             END IF
             eeta2(i,j) = Gntrp (eeta(:,:,kc), biold, bjold, 0)*EXP(-rate*dt)
          END DO
          eeta2 (isize,j) = eeta (isize,j,kc)
!
!         Points eeta(1:imin_j(j),:) are set in Get_eta_on_boundary, but:
!         eeta2 (1:imin_j(j),j) = etac(kc)
       END DO
!
       CALL Circle (eeta2)
       eeta (:, :, kc) = eeta2
!
    END DO
!
!
    RETURN
    CONTAINS
!
    FUNCTION Ratefn (fudgx, alamx, sinix, birx, vmx, xmfact)
    IMPLICIT NONE
    REAL (rprec), INTENT (IN) :: fudgx,alamx,sinix,birx,vmx,xmfact
    REAL (rprec)              :: Ratefn
!                                                                       
!   Function subprogram to compute precipitation rate
!   Last update:  04-04-88
!
    Ratefn = 0.0466_rprec*fudgx*SQRT(ABS(alamx))*(sinix/birx)*vmx**2
    Ratefn = xmfact * ratefn
    RETURN
    END FUNCTION Ratefn
    END SUBROUTINE Move_plasma_grid