      subroutine regrot(pxreg,pyreg,pxrot,pyrot,kxdim,kydim,kx,ky,
     +                  pxcen,pycen,kcall)
c
      implicit none
c
c-------------------------------------------------------------------------------
c 
c*    conversion between regular and rotated spherical coordinates.
c*
C*    pxreg     longitudes of the regular coordinates
c*    pyreg     latitudes of the regular coordinates
c*    pxrot     longitudes of the rotated coordinates
c*    pyrot     latitudes of the rotated coordinates
c*              all coordinates given in degrees N (negative values for S)
c*              and degrees E (negative values for W)
c*    kxdim     dimension of the gridpoint fields in the x-direction
c*    kydim     dimension of the gridpoint fieLds in the y-direction
c*    kx        number of gridpoint in the x-direction
c*    ky        number of gridpoints in the y-direction
c*    pxcen     regular longitude of the south pole of the rotated grid
c*    pycen     regular latitude of the south pole of the rotated grid
c*
C*    kcall=-1: find regular as functions of rotated coordinates.
c*    kcall= 1: find rotated as functions of regular coordinates.
c*
c*    j.e. haugen   hirlam   june -92
c
c-------------------------------------------------------------------------------
c 
      integer kxdim,kydim,kx,ky,kcall
      real pxreg(kxdim,kydim),pyreg(kxdim,kydim),
     +     pxrot(kxdim,kydim),pyrot(kxdim,kydim),
     +     pxcen,pycen 
c
c-------------------------------------------------------------------------------
c
      real pi,zrad,zsycen,zcycen,zxmxc,zsxmxc,zcxmxc,zsyreg,zcyreg,
     x     zsyrot,zcyrot,zcxrot,zsxrot,zradi 
      integer jy,jx
c
c-------------------------------------------------------------------------------
c
      pi = 4.*atan(1.)
      zrad = pi/180. 
      zradi = 1./zrad 
      zsycen = sin(zrad*(pycen+90.))
      zcycen = cos(zrad*(pycen+90.))
c
      if (kcall.eq.1) then
c
      do jy = 1,ky
      do jx = 1,kx
c
      zxmxc  = zrad*(pxreg(jx,jy) - pxcen)
      zsxmxc = sin(zxmxc)
      zcxmxc = cos(zxmxc)
      zsyreg = sin(zrad*pyreg(jx,jy))
      zcyreg = cos(zrad*pyreg(jx,jy))
      zsyrot = zcycen*zsyreg - zsycen*zcyreg*zcxmxc
      zsyrot = amax1(zsyrot,-1.0)
      zsyrot = amin1(zsyrot,+1.0)
c
      pyrot(jx,jy) = asin(zsyrot)*zradi 
c
      zcyrot = cos(pyrot(jx,jy)*zrad)
      zcxrot = (zcycen*zcyreg*zcxmxc +
     +          zsycen*zsyreg)/zcyrot
      zcxrot = amax1(zcxrot,-1.0)
      zcxrot = amin1(zcxrot,+1.0)
      zsxrot = zcyreg*zsxmxc/zcyrot
c
      pxrot(jx,jy) = acos(zcxrot)*zradi 
c
      if (zsxrot.lt.0.0) pxrot(jx,jy) = -pxrot(jx,jy)
c
      enddo
      enddo
c
      elseif (kcall.eq.-1) then
c
      do jy = 1,ky
      do jx = 1,kx
c
      zsxrot = sin(zrad*pxrot(jx,jy))
      zcxrot = cos(zrad*pxrot(jx,jy))
      zsyrot = sin(zrad*pyrot(jx,jy))
      zcyrot = cos(zrad*pyrot(jx,jy))
      zsyreg = zcycen*zsyrot + zsycen*zcyrot*zcxrot
      zsyreg = amax1(zsyreg,-1.0)
      zsyreg = amin1(zsyreg,+1.0)
c
      pyreg(jx,jy) = asin(zsyreg)*zradi
c
      zcyreg = cos(pyreg(jx,jy)*zrad)
      zcxmxc = (zcycen*zcyrot*zcxrot -
     +          zsycen*zsyrot)/zcyreg
      zcxmxc = amax1(zcxmxc,-1.0)
      zcxmxc = amin1(zcxmxc,+1.0)
      zsxmxc = zcyrot*zsxrot/zcyreg
      zxmxc  = acos(zcxmxc)*zradi 
      if (zsxmxc.lt.0.0) zxmxc = -zxmxc
c
      pxreg(jx,jy) = zxmxc + pxcen
c
      enddo
      enddo
c
      else
      write(6,'(1x,''invalid kcall in regrot'')')
      stop
      endif
c
      return
      end



