**********************************************************************
	subroutine uvinitp(r,xip,xjp,idim,jdim,gf)
c
c initialize for grid characteristics
c
c input
c     r   - distance (on the projection) in grid units from pole to
c            equator
c     xip - i coord of pole point on the grid
c     xjp - j coord of pole point on the grid
c     idim,jdim - dimensions of grid
c
c output
c     gf - array of same size as height and wind arrays to contain
c           calculation constants
c
c initialization parameters for common polar-sterographic grids
c
c     grid                   r                xip,xjp      idim,jdim
c     nmc 47x51              31.204359052      24,26         47,51
c     nmc 65x65              31.204359052      33,33         65,65
c     nmc lfm (41x38 set)    62.408718104      19,42         41,38
c     navy 63x63             31.204359052      32,32         63,63
c     australian 47x47       23.403269292      24,24         47,47
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))
	dimension gf(idim,jdim)
c
	a=r*0.5
	a2=a*a
c
	do i=1,idim
	  xi=real(i)-xip
	  do j=1,jdim
	    xj=real(j)-xjp
	    xx=(xi*xi)+(xj*xj)
	    sinphi=(4*a2-xx)/(4*a2+xx)
	    gf(i,j)=gw/sinphi
c	write(*,*) i,j,gw,sinphi,gf(i,j)
	  enddo
	enddo
c
	return
	end
**********************************************************************
	subroutine uvcalcp(a,u,v,gf,r,idim,jdim,spval)
c
c computes geostrophic winds from a height grid on a
c  polar-stereographic grid.  +u and +v are +x and +y directions, not
c  east and north.
c
c input
c       a - height array
c       gf- array of calculation constants
c       r   - distance (on the projection) in grid units from pole to
c              equator (see uvinitp for r's for some common grids)
c       idim,jdim - dimension of all arrays
c       spval - special value used to indicate missing heights on
c                input and missing winds on output
c
c output
c       u - u wind component (units determined by input heights -
c            i.e. - heights in meters results in winds in meters/sec)
c       v - v wind component
c
c constants
c er - radius of the earth in meters
c
c computes winds at all grid points with non-missing height values,
c  but values at edges or adjacent to missing height values are not
c  reliable.  missing height values must be set to spval.
c
	parameter (er=6371.2213e03,pi=3.141592654)
	dimension a(idim,jdim),u(idim,jdim),v(idim,jdim),gf(idim,jdim)
c
	do j=1,jdim
	  do i=1,idim
	    u(i,j)=spval
	    v(i,j)=spval
	  enddo
	enddo
c
	gse=(er*pi/2.)/r
c
c compute v
c
	do j=1,jdim
	  hi=a(1,j)
	  x1=1.
	  if (a(1,j) .ne. spval .and. a(2,j) .ne. spval) then
	    hj=0.5*(a(1,j)+a(2,j))
	    x2=1.5
	  else if (a(1,j) .ne. spval) then
	    hi=spval
	    hj=a(1,j)
	    x2=1.
	  else
	    hj=a(2,j)
	    x2=2.
	  endif
	  if (hi .ne. spval .and. hj .ne. spval) then
	    xx=(x2-x1)*gse
	    v(1,j)=gf(1,j)*(hj-hi)/xx
	  endif
	  hi=hj
	  x1=x2
	  do i=2,idim-2
	    if (a(i-1,j) .ne. spval .and. a(i,j) .ne. spval .and.
     +a(i+1,j) .ne. spval .and. a(i+2,j) .ne. spval) then
		hj=(-a(i-1,j)+9.*a(i,j)+9.*a(i+1,j)-a(i+2,j))/16.
		x2=(i+(i+1))*0.5
	    else if (a(i,j) .ne. spval .and. a(i+1,j) .ne. spval) then
		hj=0.5*(a(i,j)+a(i+1,j))
		x2=(i+(i+1))*0.5
	    else if (a(i,j) .ne. spval) then
		hj=a(i,j)
		x2=i
	    else
		hj=a(i+1,j)
		x2=i+1
	    endif
	    if (hi .ne. spval .and. hj .ne. spval) then
		xx=(x2-x1)*gse
		v(i,j)=gf(i,j)*(hj-hi)/xx
	    endif
	    hi=hj
	    x1=x2
	  enddo
	  if (a(idim-1,j) .ne. spval .and. a(idim,j) .ne. spval) then
	    hj=0.5*(a(idim-1,j)+a(idim,j))
	    x2=((idim-1)+idim)*0.5
	  else if (a(idim-1,j) .ne. spval) then
	    hj=a(idim-1,j)
	    x2=idim-1
	  else
	    hi=spval
	    hj=a(idim,j)
	    x2=idim
	  endif
	  if (hi .ne. spval .and. hj .ne. spval) then
	    xx=(x2-x1)*gse
	    v(idim-1,j)=gf(idim-1,j)*(hj-hi)/xx
	  endif
	  hi=hj
	  x1=x2
	  hj=a(idim,j)
	  x2=idim
	  if (hi .ne. spval .and. hj .ne. spval) then
	    xx=(x2-x1)*gse
	    v(idim,j)=gf(idim,j)*(hj-hi)/xx
	  endif
	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 (hi .ne. spval .and. hj .ne. spval) then
	    xx=(x2-x1)*gse
	    u(i,1)=gf(i,1)*(hi-hj)/xx
c	write(*,*) i,u(i,1),gf(i,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 (hi .ne. spval .and. hj .ne. spval) then
		xx=(x2-x1)*gse
		u(i,j)=gf(i,j)*(hi-hj)/xx
c	write(*,*) i,j,u(i,j),gf(i,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 (hi .ne. spval .and. hj .ne. spval) then
	    xx=(x2-x1)*gse
	    u(i,jdim-1)=gf(i,jdim-1)*(hi-hj)/xx
c	write(*,*) i,jdim-1,u(i,jdim-1),gf(i,jdim-1),(hi-hj),xx
	  endif
	  hi=hj
	  x1=x2
	  hj=a(i,jdim)
	  x2=jdim
	  if (hi .ne. spval .and. hj .ne. spval) then
	    xx=(x2-x1)*gse
	    u(i,jdim)=gf(i,jdim)*(hi-hj)/xx
c	write(*,*) i,jdim,u(i,jdim),gf(i,jdim),(hi-hj),xx
	  endif
	enddo
c
	return
	end