subroutine zetabar(i,j,delsqp)
c
c=======================================================================
c                                                                    ===
c  This routine computes barotropic vorticity from transport         ===
c  streamfunction at the point I,J. It is called by BOUNDARY         ===
c  to set the radiative CFVN open boundary conditions on the         ===
c  barotropic vorticity.                                             ===
c                                                                    ===
c=======================================================================
c
c-----------------------------------------------------------------------
c  Define global data.
c-----------------------------------------------------------------------
c
#include <cdefs.h>
#include <param.h>
#include <pconst.h>
#include <scalar.h>
#include <onedim.h>
#include <fullwd.h>
#include <fields.h>
#include <workspa.h>
c
c-----------------------------------------------------------------------
c  Define local data.
c-----------------------------------------------------------------------
c
      integer i,ii,j,jj
      FLOAT
     *      cfe,cfn,cfs,cfw,cpf,delsqp,fx,fxa,fxb,fxc
c
c-----------------------------------------------------------------------
c  Begin executable code.
c-----------------------------------------------------------------------
c
      jj=j
      ii=i
      if(i.eq.1) ii=2
      if(i.eq.imt) ii=imtm1
      if(j.eq.1) jj=2
      if(j.eq.jmt) jj=jmtm1
      fxa=c2*cstr(jj)*cstr(jj)
      fxb=c2*cs(jj  )*cstr(jj)*dytr(jj)*dyur(jj  )
      fxc=c2*cs(jj-1)*cstr(jj)*dytr(jj)*dyur(jj-1)
      fx=c1
      cfn=fxb/(c1/hv(ii-1,jj)+c1/hv(ii,jj))
      cfs=fxc/(c1/hv(ii-1,jj-1)+c1/hv(ii,jj-1))
      cfe=fxa*dxur(ii)*dxtr(ii)/(c1/hv(ii,jj)+c1/hv(ii,jj-1))
      cfw=fxa*dxur(ii-1)*dxtr(ii)/(c1/hv(ii-1,jj)+c1/hv(ii-1,jj-1))
      cpf=fx/(cfn+cfs+cfe+cfw)
      delsqp= p(ii+1,jj  )*cfe+p(ii-1,jj  )*cfw
     *       +p(ii  ,jj+1)*cfn+p(ii  ,jj-1)*cfs-p(ii,jj)/cpf
      return
      end