FUNCTION rk3, a, u, v, cr, WENO=weno

@c_block

  schemetitle = 'RK3-5TH WENO'

  if NOT keyword_set(WENO) then weno=0
  if NOT keyword_set(WENO) then schemetitle = 'RK3-5TH'

  dim = size(a)
  nx = dim(1)
  ny = dim(2)
  eps = 1.0e-10

; Set up dummy arrays.....

  b = a

; Flux arrays (domain is doubly periodic, so we dont need extra points)

  fx = fltarr(nx,ny)
  fy = fltarr(nx,ny)

; These are the standard weights that yield a 5th order scheme

  gi = fltarr(3)
  gi[0] = 1./10.
  gi[1] = 6./10.
  gi[2] = 3./10.

  FOR n = 0,2 DO BEGIN
 
    cra = cr / (3-n)

    cu   = cra * u[0:nx-1,0:ny-1]
    dir  = 0.5*sign(cu)
  
    ci   = b[0:nx-1,0:ny-1]
    cip3 = shift(ci,-3)
    cip2 = shift(ci,-2)
    cip1 = shift(ci,-1)
    cim1 = shift(ci, 1)
    cim2 = shift(ci, 2)
    cim3 = shift(ci, 3)

; This switches upwind directions.
     
    bip2 =  (0.5 + dir)*cip2 - (dir - 0.5)*cim1
    bip1 =  (0.5 + dir)*cip1 - (dir - 0.5)*ci
    bi   =  (0.5 + dir)*ci   - (dir - 0.5)*cip1
    bim1 =  (0.5 + dir)*cim1 - (dir - 0.5)*cip2
    bim2 =  (0.5 + dir)*cim2 - (dir - 0.5)*cip3

; Compute the three 3rd-order approximations

    f0 =  1./3.*bim2 - 7./6.*bim1 + 11./6.*bi
    f1 = -1./6.*bim1 + 5./6.*bi   + 1./3. *bip1
    f2 =  1./3.*bi   + 5./6.*bip1 - 1./6. *bip2

    IF ( weno ) THEN BEGIN

; Compute Taylor series approximations to check for smoothness

      beta0 = 13./12.*(bim2 - 2.*bim1 + bi  )^2 + 1./4.*(bim2 - 4.*bim1 + 3.*bi)^2
      beta1 = 13./12.*(bim1 - 2.*bi   + bip1)^2 + 1./4.*(bim1 - bip1)^2
      beta2 = 13./12.*(bi   - 2.*bip1 + bip2)^2 + 1./4.*(bip2 - 4.*bip1 + 3.*bi)^2
 
; Compute weights
 
      wk0 = gi[0] / (eps + beta0)^2
      wk1 = gi[1] / (eps + beta1)^2
      wk2 = gi[2] / (eps + beta2)^2
   
      sumwk = wk0 + wk1 + wk2
 
      wi0 = wk0 / sumwk
      wi1 = wk1 / sumwk
      wi2 = wk2 / sumwk

; Compute WENO flux

      fx[0:nx-1,0:ny-1] = cu * (wi0 * f0 + wi1*f1 + wi2*f2)

     ENDIF ELSE BEGIN

; Compute normal 5th-order flux

      fx[0:nx-1,0:ny-1] = cu * (gi[0]*f0 + gi[1]*f1 + gi[2]*f2)

     ENDELSE

; Y-direction
 
    cv   = cra * v[0:nx-1,0:ny-1]
    dir  = 0.5*sign(cv)

    ci   = b[0:nx-1,0:ny-1]
    cip3 = shift(ci,0,-3)
    cip2 = shift(ci,0,-2)
    cip1 = shift(ci,0,-1)
    cim1 = shift(ci,0, 1)
    cim2 = shift(ci,0, 2)
    cim3 = shift(ci,0, 3)
  
; This switches upwind directions.

    bip2 =  (0.5 + dir)*cip2 - (dir - 0.5)*cim1
    bip1 =  (0.5 + dir)*cip1 - (dir - 0.5)*ci
    bi   =  (0.5 + dir)*ci   - (dir - 0.5)*cip1
    bim1 =  (0.5 + dir)*cim1 - (dir - 0.5)*cip2
    bim2 =  (0.5 + dir)*cim2 - (dir - 0.5)*cip3
  
; Compute the three 3rd-order approximations

    f0 =  1./3.*bim2 - 7./6.*bim1 + 11./6.*bi
    f1 = -1./6.*bim1 + 5./6.*bi   + 1./3. *bip1
    f2 =  1./3.*bi   + 5./6.*bip1 - 1./6. *bip2

    IF ( weno ) THEN BEGIN

; Compute Taylor series approximations to check for smoothness

      beta0 = 13./12.*(bim2 - 2.*bim1 + bi  )^2 + 1./4.*(bim2 - 4.*bim1 + 3.*bi)^2
      beta1 = 13./12.*(bim1 - 2.*bi   + bip1)^2 + 1./4.*(bim1 - bip1)^2
      beta2 = 13./12.*(bi   - 2.*bip1 + bip2)^2 + 1./4.*(bip2 - 4.*bip1 + 3.*bi)^2

; Compute weights
 
      wk0 = gi[0] / (eps + beta0)^2
      wk1 = gi[1] / (eps + beta1)^2
      wk2 = gi[2] / (eps + beta2)^2
 
      sumwk = wk0 + wk1 + wk2
 
      wi0 = wk0 / sumwk
      wi1 = wk1 / sumwk
      wi2 = wk2 / sumwk

; Compute WENO flux

      fy[0:nx-1,0:ny-1] = cv * (wi0 * f0 + wi1*f1 + wi2*f2)

    ENDIF ELSE BEGIN

; Compute normal 5th-order flux

    fy[0:nx-1,0:ny-1] = cv * (gi[0]*f0 + gi[1]*f1 + gi[2]*f2)

    ENDELSE

;  Update field
 
    fim1 = shift(fx,1,0) 
    fjm1 = shift(fy,0,1) 

    b = a + (fim1 - fx) + (fjm1 - fy)

; End out RK loop

  ENDFOR

; RETURN NEW TIME LEVEL

return, b

END