**********************************************************************
	subroutine uvinitll(r,xjp,jdim,gf)
c
c initialize for grid characteristics
c
c input
c     r    - latitude of j = 1 (degrees)
c     xjp  - grid spacing in latitudinal direction (degrees)
c     jdim - dimension of grid in latitudinal direction
c
c output
c     gf - array of size jdim to contain calculation constants
c
	dimension gf(jdim)
c
c constants
c gw - constant equal to acc due to gravity (m/s) divided by 2 times
c       the rotational velocity of the earth
c
	parameter (g=9.80665,omega=7.292116e-5,gw=g/(2.*omega))
	parameter (pi=3.141592654,acon=pi/180.)
c
	do j=1,jdim
	  phi=r+xjp*(j-1)
	  sinphi=sin(phi*acon)
	  if (phi .ne. 0.) then
	    gf(j)=gw/sinphi
	  else
	    gf(j)=-999.9
	  endif
	enddo
c
	return
	end
**********************************************************************
	subroutine uvcalcll(a,u,v,gf,r,xip,xjp,idim,jdim,spval)
c
c computes geostrophic winds from a height grid onto a latitude/
c  longitude grid.  +u and +v are east and north, respectively.
c
c input
c     a - height array (idim by jdim)
c     gf - array of calculation constants (from uvinit)
c     r - latitude in degrees of j = 1
c     xip - grid spacing in degrees in longitudinal direction
c     xjp - grid spacing in degrees in latitudinal direction
c     idim,jdim - dimension of output wind component arrays
c     spval - special value used to indicate missing height on input
c             and missing winds on output
c
c output
c     u - u wind component (units determined by input heights -
c          i.e. - height in meters results in winds in meters/sec)
c     v - v wind component
c
c computes winds at all grid points with non-missing height values,
c  but values near the equator and near the pole are not reliable,
c  as well as values adjacent to missing height values.  missing
c  height values must be set to spval.
c
	parameter (pi=3.141592654,acon=pi/180.)
	dimension a(idim,jdim),u(idim,jdim),v(idim,jdim),gf(jdim)
	dimension b(360,360)
c
	gsi=111111.*xip
	gsj=111111.*xjp
c
	do j=1,jdim
	  do i=1,idim
	    u(i,j)=spval
	    v(i,j)=spval
	  enddo
	enddo
c
c extend height grid on both ends
c
	do j=1,jdim
	  do i=1,idim
	    b(i+2,j)=a(i,j)
	  enddo
	  b(1,j)=a(idim-1,j)
	  b(2,j)=a(idim,j)
	  b(idim+3,j)=a(1,j)
	  b(idim+4,j)=a(2,j)
	enddo
c
c compute v
c
	do j=1,jdim
	  if (b(1,j) .ne. spval .and. b(2,j) .ne. spval .and. b(3,j)
     +.ne. spval .and. b(4,j) .ne. spval) then
	    hi=(-b(1,j)+9.*b(2,j)+9.*b(3,j)-b(4,j))/16.
	    x1=2.5
	  else if (b(2,j) .ne. spval .and. b(3,j) .ne. spval) then
	    hi=0.5*(b(2,j)+b(3,j))
	    x1=2.5
	  else if (b(2,j) .ne. spval) then
	    hi=b(2,j)
	    x1=2.
	  else
	    hi=b(3,j)
	    x1=3.
	  endif
	  do i=1,idim
	    if (b(i+1,j) .ne. spval .and. b(i+2,j) .ne. spval .and.
     +b(i+3,j) .ne. spval .and. b(i+4,j) .ne. spval) then
		hj=(-b(i+1,j)+9.*b(i+2,j)+9.*b(i+3,j)-b(i+4,j))/16.
		x2=((i+2)+(i+3))*0.5
	    else if (b(i+2,j) .ne. spval .and. b(i+3,j) .ne. spval)
     +then
		hj=0.5*(b(i+2,j)+b(i+3,j))
		x2=((i+2)+(i+3))*0.5
	    else if (b(i+2,j) .ne. spval) then
		hj=b(i+2,j)
		x2=i+2
	    else
		hj=b(i+3,j)
		x2=i+3
	    endif
	    if (hi .ne. spval .and. hj .ne. spval .and. gf(j) .ne.
     +-999.9) then
		phi=(r+xjp*(j-1))*acon
		xx=(x2-x1)*gsi*cos(phi)
		v(i,j)=gf(j)*(hj-hi)/xx
	    endif
	    hi=hj
	    x1=x2
	  enddo
	enddo
c
c compute u
c
	do i=1,idim
	  hi=a(i,1)
	  x1=1.
	  if (a(i,1) .ne. spval .and. a(i,2) .ne. spval) then
	    hj=0.5*(a(i,1)+a(i,2))
	    x2=1.5
	  else if (a(i,1) .ne. spval) then
	    hi=spval
	    hj=a(i,1)
	    x2=1.
	  else
	    hj=a(i,2)
	    x2=2.
	  endif
	  if (gf(1) .ne. -999.9 .and. hi .ne. spval .and. hj .ne.
     +spval) then
	    xx=(x2-x1)*gsj
	    u(i,1)=gf(1)*(hi-hj)/xx
	  endif
	  hi=hj
	  x1=x2
	  do j=2,jdim-2
	    if (a(i,j-1) .ne. spval .and. a(i,j) .ne. spval .and.
     +a(i,j+1) .ne. spval .and. a(i,j+2) .ne. spval) then
		hj=(-a(i,j-1)+9.*a(i,j)+9.*a(i,j+1)-a(i,j+2))/16.
		x2=(j+(j+1))*0.5
	    else if (a(i,j) .ne. spval .and. a(i,j+1) .ne. spval) then
		hj=0.5*(a(i,j)+a(i,j+1))
		x2=(j+(j+1))*0.5
	    else if (a(i,j) .ne. spval) then
		hj=a(i,j)
		x2=j
	    else
		hj=a(i,j+1)
		x2=j+1
	    endif
	    if (gf(j) .ne. -999.9 .and. hi .ne. spval .and. hj .ne.
     +spval) then
		xx=(x2-x1)*gsj
		u(i,j)=gf(j)*(hi-hj)/xx
	    endif
	    hi=hj
	    x1=x2
	  enddo
	  if (a(i,jdim-1) .ne. spval .and. a(i,jdim) .ne. spval) then
	    hj=0.5*(a(i,jdim-1)+a(i,jdim))
	    x2=((jdim-1)+jdim)*0.5
	  else if (a(i,jdim-1) .ne. spval) then
	    hj=a(i,jdim-1)
	    x2=jdim-1
	  else
	    hi=spval
	    hj=a(i,jdim)
	    x2=jdim
	  endif
	  if (gf(jdim-1) .ne. -999.9 .and. hi .ne. spval .and. hj .ne.
     +spval) then
	    xx=(x2-x1)*gsj
	    u(i,jdim-1)=gf(jdim-1)*(hi-hj)/xx
	  endif
	  hi=hj
	  x1=x2
	  hj=a(i,jdim)
	  x2=jdim
	  if (gf(jdim) .ne. -999.9 .and. hi .ne. spval .and. hj .ne.
     +spval) then
	    xx=(x2-x1)*gsj
	    u(i,jdim)=gf(jdim)*(hi-hj)/xx
	  endif
	enddo
c
	return
	end