!=========================================================================
!=========================================================================
!=========================================================================
!
SUBROUTINE Move_plasma_grid_NEW (dt, ie_ask)
  USE Rcm_mod_subs
  IMPLICIT NONE
  REAL (rprec), INTENT (IN) :: dt
  INTEGER (iprec), INTENT (IN) :: ie_ask
!_____________________________________________________________________________
!   Subroutine to advance eta distribution for a time step
!   by using new CLAWPACK advection routines
!                                                                       
!   Created:     12-05-00
!_____________________________________________________________________________
!
!
  REAL (rprec) :: eeta2 (isize,jsize), veff  (isize,jsize),  &
                  dvefdi(isize,jsize), dvefdj(isize,jsize), &
                  didt, djdt, biold, bjold, rate, &
                  mass_factor
  INTEGER (iprec) :: i, j, kc, ie
!\\\
  integer :: CLAWiter, joff, icut
  real, dimension(-1:isize+2,-1:jsize-1):: loc_didt, loc_djdt
  real, dimension(-1:isize+2,-1:jsize-1):: loc_Eta, loc_rate
  double precision xlower,xupper,ylower,yupper, T1,T2
  logical, save :: FirstTime=.true.
  REAL (rprec) :: r_dist, eta_scale
  
  
  joff=jwrap-1
  
  if (FirstTime) then
     T1=0.
  else
     T1=T2
  end if
  T2=T1+dt
!///
!
!
!
  DO kc = 1, kcsize
!
!    If oxygen is to be added, must change this!
!
     ie = ikflavc (kc)
!
     if (ie /= ie_ask) CYCLE
!
     IF (kc > 1) THEN
        eta_scale = MAXVAL( eeta(:,:,kc))
     ELSE
        eta_scale = 1.0
     END IF
     eta_scale = 1.0
     
!    IF (eta_scale == 0.0) CYCLE
!
     mass_factor = SQRT (xmass(1)/xmass(ie))
!
!    1. Compute the effective potential for the kc energy channel:
!
     veff = v +vcorot - vpar + vm*alamc(kc)
!%   CALL V_eff_polar_cap (veff) !12/20/05, Sazykin: New B.C. on V specify it above the boundary
!
!    2. Differentiate Veff with respect to I and J:
!
     CALL Deriv_i_new (veff, isize, jsize, j1, j2, imin_J, dvefdi)
     CALL Deriv_j_new (veff, isize, jsize, j1, j2, imin_J, dvefdj)
!
!
     loc_Eta  = zero
     loc_didt = zero
     loc_djdt = zero
     loc_rate = zero
!
!
      icut=0
      do j=j1,j2
         icut=max(icut,imin_j(j))
         do i=imin_j(j),i2
            if (eeta(i,j,kc) > 1.0/eta_scale) icut=max(icut,i)
         end do
      end do
      icut=icut+5
!
      fac = 1.0E-3*signbe*bir*alpha*beta*dlam*dpsi*ri**2
!
     DO j = j1, j2
        DO i = 2, isize-1
           loc_didt (i,j-joff) = dvefdj (i-1,j) / fac(i-1,j)
           loc_djdt (i,j-joff) = - dvefdi (i,j-1) / fac(i-1,j)
           IF (i > icut) THEN
              loc_didt(i,j-joff) = 0.0
              loc_djdt(i,j-joff) = 0.0
           END IF
!
!          Determine the rate of loss for a given species:
!
           IF (ie == 1) THEN
              IF (kc == 1) THEN
!                Cold Electrons (plasmasphere). No pitch-angle scattering,
!                but account for the refilling from the ionosphere 
!                rate is > 0, so don't use is as a loss term!!!
                 if (i >= imin_j(j)) then
                 loc_rate(i,j-joff) = &
                    Plasmasphere_refill_rate (rmin(i,j), doy, sunspot_number,&
                                              eeta(i,j,kc), vm(i,j))
!                eeta (i,j,kc) = eeta(i,j,kc)+(loc_rate(i,j-joff)/eta_scale)*dt
                 eeta (i,j,kc) = eeta(i,j,kc)+(loc_rate(i,j-joff))*dt
                 end if
                 loc_rate(i,j-joff) = 0.0
              ELSE
!                Plasmasheet electrons, strong pitch-angle scattering loss:
                 loc_rate(i,j-joff) = 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 if it is on:
              IF (L_dktime .AND. i >= imin_j(j)) THEN
              r_dist = SQRT(xmin(i,j)**2+ymin(i,j)**2)
              loc_rate(i,j-joff) = Cexrat (ie, ABS(alamc(kc))*vm(i,j), &
                                           R_dist, &
                                           sunspot_number, dktime, &
                                           irdk,inrgdk,isodk,iondk)
              ELSE
                 loc_rate(i,j-joff) = 0.0
              END IF
           END IF
        END DO
!       eeta (1:imin_j(j)-1,j,kc) = etac (kc)
        loc_didt(isize,j-joff) = loc_didt(isize-1,j-joff)
        loc_djdt(isize,j-joff) = loc_djdt(isize-1,j-joff)
        loc_rate(isize,j-joff) = loc_rate(isize-1,j-joff)
     END DO
!
!
!Copy to local variables
     loc_Eta (1:isize, 1:jsize-jwrap) = eeta (1:isize, jwrap:jsize-1, kc)
!    !Set ghost cell values for clawpack solver
!    !  Pole
!    do i=1-2, 1-1
!       loc_Eta (i,j1-joff:j2-joff) = loc_Eta (1,j1-joff:j2-joff)
!       loc_didt(i,j1-joff:j2-joff) = loc_didt(1,j1-joff:j2-joff)
!       loc_djdt(i,j1-joff:j2-joff) = loc_djdt(1,j1-joff:j2-joff)
!    end do
!    !  Equator
!    do i=isize+1,isize+2
!       loc_Eta (i,j1-joff:j2-joff) = loc_Eta (isize,j1-joff:j2-joff)
!       loc_didt(i,j1-joff:j2-joff) = loc_didt(isize,j1-joff:j2-joff)
!       loc_djdt(i,j1-joff:j2-joff) = loc_djdt(isize,j1-joff:j2-joff)
!    end do
!    !  Periodic
!    loc_Eta (-1:isize+1,-1:0) = loc_Eta (-1:isize+1,jsize-4:jsize-3)
!    loc_didt(-1:isize+1,-1:0) = loc_didt(-1:isize+1,jsize-4:jsize-3)
!    loc_djdt(-1:isize+1,-1:0) = loc_djdt(-1:isize+1,jsize-4:jsize-3)
!    loc_Eta (-1:isize+1,jsize-joff:jsize-joff+1) = loc_Eta (-1:isize+1,1:2)
!    loc_didt(-1:isize+1,jsize-joff:jsize-joff+1) = loc_didt(-1:isize+1,1:2)
!    loc_djdt(-1:isize+1,jsize-joff:jsize-joff+1) = loc_djdt(-1:isize+1,1:2)
!    
     !Call clawpack
     xlower = 1
     xupper = isize
     ylower = zero
     yupper = jsize-3
     FirstTime=.true.
     CALL Claw2ez (FirstTime, T1,T2, xlower,xupper, ylower,yupper, &
                   CLAWiter, 2,isize-1+1,jsize-3, &
                   loc_Eta, loc_didt, loc_djdt, loc_rate)
     FirstTime=.false.
 
     !Copy out
     DO j = j1, j2
     DO i = imin_j(j)+1, isize - 1
!       eeta (i, j, kc) = loc_Eta (i, j-joff) * eta_scale
        eeta (i, j, kc) = loc_Eta (i, j-joff)
     END DO
     END DO
!
     eeta(:,:,kc) = eeta(:,:,kc) * eta_scale
     CALL Circle (eeta(:,:,kc))
!
  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
  
  SUBROUTINE Deriv_i_NEW (array, isize, jsize, j1, j2, imin_j, derivi)
    USE Rcm_mod_subs, ONLY : iprec, rprec
    IMPLICIT NONE
    INTEGER (iprec), INTENT (IN) :: isize, jsize, j1, j2, imin_j(jsize)
    REAL (rprec), INTENT (IN) :: array (isize,jsize)
    REAL (rprec), INTENT (OUT) :: Derivi (isize,jsize)
!
    INTEGER (iprec) :: i, j, idim, jdim
!
    DO j = 1, jsize
    DO i = 1, isize
       IF (i == 1) THEN
          Derivi(i,j) = -1.5*array(i,j) + 2.0*array(i+1,j) - 0.5*array(i+2,j)
       ELSE IF (i == isize) THEN
          Derivi (i,j) =  +1.5*array(i,j) - 2.0*array(i-1,j) + 0.5*array(i-2,j)
       ELSE
          Derivi(i,j) = 0.5*(array(i+1,j)-array(i-1,j))
       END IF
    END DO
    END DO
    RETURN
  END SUBROUTINE Deriv_i_NEW
  
  SUBROUTINE Deriv_j_NEW (array, isize, jsize, j1, j2, imin_j, derivJ)
    USE rcm_mod_subs, ONLY : iprec, rprec
    IMPLICIT NONE
    INTEGER (iprec), INTENT (IN) :: isize, jsize, j1, j2, imin_j(jsize)
    REAL (rprec), INTENT (IN) :: array (isize,jsize)
    REAL (rprec), INTENT (OUT) :: Derivj (isize,jsize)
!
    INTEGER (iprec) :: i, j
!
    DO j = j1, j2
    DO i = 1, isize
       Derivj (i,j) = (array(i,j+1)-array(i,j-1))*0.5
    END DO
    END DO
    CALL Circle (Derivj)
    RETURN
  END SUBROUTINE Deriv_j_NEW
  
END SUBROUTINE Move_plasma_grid_NEW