subroutine plume2dyn(qplume,Nxplume,Lmplume,uref,vref,flag,
     . idim1,idim2,jdim1,jdim2,i1,i2,j1,j2,Nsx,Nsy,bi,bj,qdyn1,qdyn2)
C***********************************************************************
C Purpose:
C   To interpolate an arbitrary quantity from higher resolution plume
C         grid to the model's dynamics grid
C Algorithm:
C   Plumes -> Dynamics computes the plumes mean value, and in the case
C             of a vector field, preserves the direction of a vector 
C             given in (uref,vref)
C
C Input:
C   qplume... [idim2,jdim2,im,Lmplume,bi] Quantity on Input Grid
C   Nxplume . Longitude Dimension of Input
C   Lmplume.. Vertical  Dimension of Input
C   uref .... [im,jm,Lmplume,bi,bj] Reference u-component of velocity
C   vref .... [im,jm,plume,bi,bj] Reference v-component of velocity
C   flag .... Flag to indicate vector (1) or scalar (0) interpolation
C   idim1,2.. Beginning and ending dimension of output grid
C   jdim1,2.. Beginning and ending dimension of output grid
C   i1,2..... Beginning and ending x-direction span         
C   j1,2..... Beginning and ending y-direction span         
C   Nsx...... Number of processes in x-direction
C   Nsy...... Number of processes in y-direction
C   bi....... Index of process number in x-direction
C   bj....... Index of process number in x-direction
C
C Output:
C   qdyn1..... [im,jm,plume,bi,bj] Field at output grid (dynamics)
C   qdyn2..... [im,jm,plume,bi,bj] Field at output grid (dynamics)
C
C Notes:
C   1)  Assume (for now) that the number of vertical levels is the
C       same on both the input and output grids
C***********************************************************************
      implicit none
#include "CPP_OPTIONS.h"

      integer  Nxplume, Lmplume, Nsx, Nsy
      integer idim1, idim2, jdim1, jdim2, i1, i2, j1, j2
      integer bi, bj, flag
      _RL qplume(i2,j2,Nxplume,Lmplume,Nsx)
      _RL uref(idim1:idim2,jdim1:jdim2,Lmplume,Nsx,Nsy)
      _RL vref(idim1:idim2,jdim1:jdim2,Lmplume,Nsx,Nsy)
      _RL qdyn1(idim1:idim2,jdim1:jdim2,Lmplume,Nsx,Nsy)
      _RL qdyn2(idim1:idim2,jdim1:jdim2,Lmplume,Nsx,Nsy)

      integer i,j,L,iplume
      _RL qplumeav(i2,j2,Lmplume)
      _RL sqrtarg

C First step - compute the average of qplume over Nxplume
      do j = j1,j2
      do i = i1,i2
       do L = 1,Lmplume
        qplumeav(i,j,L) = 0.
        do iplume = 1,Nxplume
         qplumeav(i,j,L)=qplumeav(i,j,L)+qplume(i,j,iplume,L,bi)/Nxplume
        enddo
       enddo
      enddo
      enddo

C Now check the flag -- if a scalar, we are done - just assign
C the average to all the i and j points of the output grid.
C If a vector, there is some more work to do in order to preserve
C the angle given by uref and vref

      if (flag.eq.0) then
       do j = j1,j2
       do i = i1,i2
       do L = 1,Lmplume
        qdyn1(i,j,L,bi,bj) = qplumeav(i,j,L)
       enddo
       enddo
       enddo
      elseif (flag.eq.1) then
       do j = j1,j2
       do i = i1,i2
       do L = 1,Lmplume
        if(vref(i,j,L,bi,bj).ne.0.) then
         sqrtarg = (qplumeav(i,j,L)*qplumeav(i,j,L)) /
     .  ( ( (uref(i,j,L,bi,bj)*uref(i,j,L,bi,bj)) /
     .      (vref(i,j,L,bi,bj)*vref(i,j,L,bi,bj)) ) + 1. )
         qdyn2(i,j,L,bi,bj) = sqrt(sqrtarg)
         qdyn1(i,j,L,bi,bj) = qdyn2(i,j,L,bi,bj) *
     .                            (uref(i,j,L,bi,bj)/vref(i,j,L,bi,bj))
        else
         qdyn1(i,j,L,bi,bj) = qplumeav(i,j,L)
         qdyn2(i,j,L,bi,bj) = 0.
        endif
       enddo
       enddo
       enddo
      endif

      return
      end