subroutine xtr_psi(lon,lat,map_p,imax)
c
c=======================================================================
c                                                                    ===
c  This subroutine extracts transport streamfunction and then        ===
c  interpolates it to the sub-domain grid (lon,lat).                 ===
c                                                                    ===
c  On Input:                                                         ===
c                                                                    ===
c     LON     sub-domain T-point longitude (degrees west, real array)===
c     LAT     sub-domain T-point latitude (degrees north, real array)===
c     IMAX    number of points in the x-direction to interpolate     ===
c             (integer)                                              ===
c                                                                    ===
c  On Output:                                                        ===
c                                                                    ===
c     MAP_P   sub-domain transport streamfunction (real array)       ===
c                                                                    ===
c  Calls:  LL2XY                                                     ===
c                                                                    ===
c=======================================================================
c
c-----------------------------------------------------------------------
c  Define global data.
c-----------------------------------------------------------------------
c
#include <cdefs.h>
#include <param.h>
#include <pconst.h>
#include <fields.h>
#include <moddat.h>
c
c-----------------------------------------------------------------------
c  Define local  data.
c-----------------------------------------------------------------------
c
      integer i,ic,imax,j,l
      integer icell(4),jcell(4)
      FLOAT
     *      x,xdis,y,ydis
      FLOAT
     *      lat(ximt),lon(ximt),map_p(ximt),p_cell(4)
c
c=======================================================================
c  Begin executable code.
c=======================================================================
c
c-----------------------------------------------------------------------
c  Extract and linearly interpolate transport streamfunction.
c-----------------------------------------------------------------------
c
      do 20 l=1,imax
c
c  Find indices for the current horizontal grid cell.
c
        call ll2xy (lon(l),lat(l),coord,imt,jmt,gridx,gridy,rlngd,rlatd,
     *              delx,dely,thetad,x,y)
        i=int(x)
        j=int(y)
        xdis=x-FLoaT(i)
        ydis=y-FLoaT(j)
c
        icell(1)=i
        icell(2)=i+1
        icell(3)=i+1
        icell(4)=i
        jcell(1)=j
        jcell(2)=j
        jcell(3)=j+1
        jcell(4)=j+1
c
c  Extract data for the current grid cell.
c
        do 10 ic=1,4
          p_cell(ic)=pb(icell(ic),jcell(ic))
  10    continue
c
c  Perform horizontal interpolation.
c
        map_p(l)=(c1-ydis)*((c1-xdis)*p_cell(1)+xdis*p_cell(2))+
     *           ydis*((c1-xdis)*p_cell(4)+xdis*p_cell(3))
  20  continue
      return
      end