c     @(#)curl.f	1.3
c
      subroutine curl
c
c
c
c   routine for calculating the surface velocity from the wind stress.
c
c
c_______________________________________________________________________________
c
      include'mladv.inc'
c
c_______________________________________________________________________________
c
      common/parm/kz,m,mm1,mm2,mp1,msq,tim,tmax,tstart,hy,rea,dt,
     +   xbasin,yy(mx0),icalc,kzp1,kzm1,iftop,ifbot,iftvv,ifrel,theta,
     +   ifsbl,ifeva

c
      common /cnn/n,nm1,nm2,np1

c
      common /extpar/ rlat0, rlng0, t0, v0, dhor, ht, time0,
     +                r0,f0star

      common/threed/hz(k0),scsig,sigz(k01),eigval(k0),ainv(k0,k0),
     +      amat(k0,k0)

      common/dynam/alpha,beta,rf

      common/fields/vort(mn0,k0),stream(mn0,k0),prvjac(mn0,k0),
     +     topden(mn0),botden(mn0),relief(mn0),prvtop(mn0),
     +     prvbot(mn0),topfct(k0),botfct(k0),pvort(mn0)

      dimension zet0(mn0,1),psi0(mn0),dx2psi(mn0),dy2psi(mn0),
     .          dxzet0(mn0),dyzet0(mn0),
     .          dxtaux(mn0),dytaux(mn0),
     .          dxtauy(mn0),dytauy(mn0)

c
c_______________________________________________________________________
c
c  now determine the surface velocity wsurf....
c
c
      if (iftvv.eq.0) then
         do 4 ip=1,msq
            wsurf(ip) = 0.0
  4       continue
         return
       endif

      fact= 1.0 / (dhor*f0star)
      fact2 = - (v0/dhor) / f0star
      bstar = beta / (dhor*t0)
      fact3 = fact2 * bstar / (v0/dhor)
      fact4 = 0.5 * hz(1) / (scsig*sigz(1))

      do 1 ip=1,msq
         psi0(ip) = stream(ip,1) + fact4 * topden(ip)
  1    continue

      call agrad2(psi0,dx2psi,dy2psi,m,n,hy)

      do 2 ip=1,msq
         zet0(ip,1) = dx2psi(ip) + dy2psi(ip)
  2    continue

      call agrad(zet0,dxzet0,dyzet0,m,n,hy)
      call agrad(taux,dxtaux,dytaux,m,n,hy)
      call agrad(tauy,dxtauy,dytauy,m,n,hy)

      do 3 ip=1,msq
         wsurf(ip) = fact * ( dxtauy(ip) - dytaux(ip) 
     .            + fact2 * ( tauy(ip) * dxzet0(ip) 
     .                      - taux(ip) * dyzet0(ip) 
     .                      - fact3 * taux(ip) ) ) 
  3    continue

      return
      end
