PROGRAM curl

        integer nr,nc
        integer lshsec,ldwork,lvhaec

        open(unit=52,status='old')
        read(unit=52,fmt=*) nc,nr

        lwork=nr*nc*10
        lshsec=nr*nc*10
        ldwork=4*(nr+1)
        lvhaec=nr*nc*10

	call curl1(nc,nr,lwork,lshsec,ldwork,lvhaec)

	end

C -------------------------------------------------------------

	subroutine curl1(nc,nr,lwork,lshsec,ldwork,lvhaec)

	IMPLICIT NONE

	integer i,j,nr,nc,nlat,nlon,isym,nt,ids,jds,lwork,ierror
	integer lshsec,ldwork,lvhaec,ndab,mdab,ityp

	real u(nr,nc),v(nr,nc)
	real vort(nr,nc)
	real xlmbda,work(lwork),pertrb
	real wshsec(lshsec),wvhaec(lvhaec),br(nc,nc),bi(nc,nc),cr(nc,nc),ci(nc,nc)
	double precision dwork(ldwork)

	open(unit=52,status='old')
	do i=1,nr
	do j=1,nc
  	  read(unit=52,fmt=*) u(i,j),v(i,j)
	enddo
	enddo

	nlat=nr
	nlon=nc
	isym=0
	nt=1
	xlmbda=0
	ids=nlat
	jds=nlon
	ndab=nlat
	mdab=nlon
	ityp=0

	call vhaeci(nlat,nlon,wvhaec,lvhaec,dwork,ldwork,ierror)
	print *,'vhaeci err: ',ierror

	call shseci(nlat,nlon,wshsec,lshsec,dwork,ldwork,ierror)
	print *,'shseci err: ',ierror

	call vhaec(nlat,nlon,ityp,nt,v,u,nlat,nlon,br,bi,cr,ci,
     &             nlon,nlon,wvhaec,lvhaec,work,lwork,ierror)

	print *,'vhaec err: ',ierror

C nb: br,bi aren't used because (I think) the curl-bit is in cr,ci
C     and the div-bit in br,bi (see also div.f)

	call vrtec(nlat,nlon,isym,nt,vort,nlat,nlon,cr,ci,nlon,nlon,
     &             wshsec,lshsec,work,lwork,ierror)

	print *,'vrtec err code: ',ierror

	open(unit=51)
	do i=1,nr
	do j=1,nc
  	  write(unit=51,fmt=*) vort(i,j)
	enddo
	enddo
	
	end