subroutine step_plumes_corr (myTime, myIter, myThid)
c----------------------------------------------------------------------
c  Subroutine step_plumes_corr - 'Wrapper' routine to advance
c        the plumes state and make the new value. 
c        At this point, increment with the "correction term"
c        which includes the dynamics tendency and the integral
c        constraint to enforce agreement with the dynamics state
c        Also: Set up "bi, bj loop" and some timers and clocks here.
c
c Call:plume2dyn (4) (interpolate plumes state to dynamics grid 
c                           for use in the correction terms)
c      AtoC         (convert plumes state on dynamics grid to C-Grid)
c      CtoA         (convert correction term on dynamics grid to A-Grid)
c      dyn2plume (4) (interpolate A-Grid correction term to plumes grid)
c      step_plumes (advance plumes state by correction term)
c-----------------------------------------------------------------------
       implicit none
#include "CPP_OPTIONS.h"
#include "SIZE.h"
#include "GRID.h"
#include "plumes_SIZE.h"
#include "DYNVARS.h"
#include "plumes_coms.h"
#include "gridalt_mapping.h"
#include "EEPARAMS.h"

       integer myTime, myIter, myThid

       _RL udyntemp(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
       _RL vdyntemp(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
       _RL thdyntemp(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
       _RL sdyntemp(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
       _RL uplumetemp(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
       _RL vplumetemp(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
       _RL thplumetemp(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
       _RL splumetemp(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)

       integer i, j, L, Lbotij, bi, bj
       integer im1, im2, jm1, jm2, idim1, idim2, jdim1, jdim2
       _RL dt

       im1 = 1-OLx
       im2 = sNx+OLx
       jm1 = 1-OLy
       jm2 = sNy+OLy
       idim1 = 1
       idim2 = sNx
       jdim1 = 1
       jdim2 = sNy
       dt = 1.

       CALL TIMER_START('PLUME2DYN         [STEP_PLUME_CORR]',mythid)
       do bj = myByLo(myThid), myByHi(myThid)
       do bi = myBxLo(myThid), myBxHi(myThid)

c Compute correction term (new dyn state-plume state to dyn) on plumes grid:
c    First: interp plumes state to dynamics grid
        call plume2dyn(uplume,im1,im2,jm1,jm2,Nr,Nsx,Nsy,
     .        1,sNx,1,sNy,bi,bj,Nr,udyntemp)
        call plume2dyn(vplume,im1,im2,jm1,jm2,Nr,Nsx,Nsy,
     .        1,sNx,1,sNy,bi,bj,Nr,vdyntemp)
        call plume2dyn(thplume,im1,im2,jm1,jm2,Nr,Nsx,Nsy,
     .        1,sNx,1,sNy,bi,bj,Nr,thdyntemp)
        call plume2dyn(splume,im1,im2,jm1,jm2,Nr,Nsx,Nsy,
     .        1,sNx,1,sNy,bi,bj,Nr,sdyntemp)

       enddo
       enddo
       CALL TIMER_STOP('PLUME2DYN         [STEP_PLUME_CORR]',mythid)

c    Second: Convert plumes state on dynamics grid to C-Grid

       CALL TIMER_START('ATOC              [STEP_PLUME_CORR]',mythid)
       call AtoC(myThid,udyntemp,vdyntemp,maskC,im1,im2,jm1,jm2,Nr,
     .                     Nsx,Nsy,1,sNx,1,sNy,udyntemp,vdyntemp)
       CALL TIMER_STOP('ATOC              [STEP_PLUME_CORR]',mythid)

c    Third: Subtract plumes state on dyn. grid from new dynamics state
       do bj = myByLo(myThid), myByHi(myThid)
       do bi = myBxLo(myThid), myBxHi(myThid)

        do L = 1,Nr
        do j = jdim1,jdim2
        do i = idim1,idim2
        udyntemp(i,j,L,bi,bj)=uvel(i,j,L,bi,bj)-udyntemp(i,j,L,bi,bj)
        vdyntemp(i,j,L,bi,bj)=vvel(i,j,L,bi,bj)-vdyntemp(i,j,L,bi,bj)
        thdyntemp(i,j,L,bi,bj)=theta(i,j,L,bi,bj)-thdyntemp(i,j,L,bi,bj)
        sdyntemp(i,j,L,bi,bj)=salt(i,j,L,bi,bj)-sdyntemp(i,j,L,bi,bj)
        enddo
        enddo
        enddo

       enddo
       enddo

c    Fourth: Convert correction terms to A-Grid
       CALL TIMER_START('CTOA              [STEP_PLUME_CORR]',mythid)
        call CtoA(myThid,udyntemp,vdyntemp,maskW,maskS,im1,im2,jm1,jm2,
     .     Nr,Nsx,Nsy,1,sNx,1,sNy,udyntemp,vdyntemp)
       CALL TIMER_STOP('CTOA              [STEP_PLUME_CORR]',mythid)

c    Fifth: Interpolate correction terms to plumes grid
       CALL TIMER_START('DYN2PLUME         [STEP_PLUME_CORR]',mythid)
       do bj = myByLo(myThid), myByHi(myThid)
       do bi = myBxLo(myThid), myBxHi(myThid)

        call dyn2plume(udyntemp,pedyn,im1,im2,jm1,jm2,Nr,Nsx,Nsy,1,sNx,
     .      1,sNy,bi,bj,uplumetemp)
        call dyn2plume(vdyntemp,pedyn,im1,im2,jm1,jm2,Nr,Nsx,Nsy,1,sNx,
     .      1,sNy,bi,bj,vplumetemp)
        call dyn2plume(thdyntemp,pedyn,im1,im2,jm1,jm2,Nr,Nsx,Nsy,1,sNx,
     .     1,sNy,bi,bj,thplumetemp)
        call dyn2plume(sdyntemp,pedyn,im1,im2,jm1,jm2,Nr,Nsx,Nsy,1,sNx,
     .      1,sNy,bi,bj,splumetemp)
       enddo
       enddo
       CALL TIMER_STOP('DYN2PLUME         [STEP_PLUME_CORR]',mythid)

c    Last: Increment plumes state by the correction term
       do bj = myByLo(myThid), myByHi(myThid)
       do bi = myBxLo(myThid), myBxHi(myThid)
        call step_plumes(uplume,vplume,thplume,splume,dt,
     .   Nxplumes,Nyplumes,Nr,Nsx,Nsy,
     .                    uplumetemp,vplumetemp,thplumetemp,splumetemp)

       enddo
       enddo

       return
       end