subroutine edotd (lon,lat, ubar,vbar,u,v, ED )
	dimension ubar(lon,lat),vbar(lon,lat), y(37),cy(37)
	dimension u(lon,lat), v(lon,lat),  ED(144,37)
	dimension dux(lon,lat),dvx(lon,lat),duy(lon,lat),dvy(lon,lat)

	pi = 4.*atan(1.)
	er = 637000.0
	detx=2.*pi/144.
	dety=pi/72.
	do j=1,37 
	   y(j)= float(j-1)*dety
	   cy(j)=cos(y(j))
	end do

	do j=5,lat-4
	do i=2,lon-1
	   dux(i,j)=(ubar(i+1,j)-ubar(i-1,j))/2./detx/er/cy(j)
	   dvx(i,j)=(vbar(i+1,j)-vbar(i-1,j))/2./detx/er/cy(j)
	   duy(i,j)=(ubar(i,j+1)-ubar(i,j-1))/2./dety/er
	   dvy(i,j)=(vbar(i,j+1)-vbar(i,j-1))/2./dety/er
	end do
	   i=1
	   dux(i,j)=(ubar(i+1,j)-ubar(144,j))/2./detx/er/cy(j)
	   dvx(i,j)=(vbar(i+1,j)-vbar(144,j))/2./detx/er/cy(j)
	   duy(i,j)=(ubar(i,j+1)-ubar(i,j-1))/2./dety/er
	   dvy(i,j)=(vbar(i,j+1)-vbar(i,j-1))/2./dety/er
	end do

	do j=5,lat-4
	do i=1,144
 	   ut = u(i,j)-ubar(i,j)
 	   vt = v(i,j)-vbar(i,j)

           ed1 = -ut*( ut*dux(i,j) + vt*duy(i,j) )
           ed2 = -vt*( ut*dvx(i,j) + vt*dvy(i,j) )
 	   ED(i,j) = ed1 + ed2
	end do
	end do

	return
	end