subroutine rotate (glt, gln, phir, thetar, psir, rlt, rln) 
c 
c     double precision translator for drotate 
c 
c 
      real glt, gln, phir, thetar, psir, rlt, rln 
      real*8 dglt, dgln, dphir, dthetar, dpsir, drlt, drln 
c 
      dglt = glt 
      dgln = gln 
      dphir = phir 
      dthetar = thetar 
      dpsir = psir 
      drlt = rlt 
      drln = rln 
      call drotate (dglt, dgln, dphir, dthetar, dpsir, drlt, drln) 
      rlt = drlt 
      rln = drln 
c 
      return 
      end 
  

      subroutine drotate (glt, gln, phir, thetar, psir, rlt, rln) 
c 
c======================================================================= 
c     subroutine drotate takes a geographic latitude and longitude and 
c     finds the the equivalent latitude and longitude on a rotated grid. 
c     when going from a geographic grid to a rotated grid, all of the 
c     defined rotation angles given to drotate by the calling program 
c     are positive, but when going from a rotated grid back to the 
c     geographic grid, the calling program must reverse the angle order 
c     (phir and psir are switched) and all of the angles made negative. 
c 
c     the first rotation angle phir is defined as a rotation about the 
c     original z axis. the second rotation angle thetar is defined as a 
c     rotation about the new x axis. the final rotation angle psir is 
c     defined as a rotation about the new z axis. these rotation angles 
c     are just the Euler angles as defined in "classical mechanics" 
c     Goldstein (1951). 
c 
c     author:   M. Eby            e-mail [email protected] 
c======================================================================= 
c 
c 
      real*8 glt, gln, phir, thetar, psir, phis, thetas, rlt, rln 
     &,      gy, gx, gz, rx, ry, rz, d0, d1, d90, d180, d360, rad 
c 
c     g...  = geographic value 
c     r...  = rotated value 
c     ...lt = latitude (or equivalent spherical coordinate) 
c     ...ln = longitude (or equivalent spherical coordinate) 
c     ...x  = x coordinate 
c     ...y  = y coordinate 
c     ...z  = z coordinate 
c     psir, thetar, phir = Euler angles defining rotation 
c 
      d0   = 0.D+00 
      d1   = 1.D+00 
      d90  = 90.D+00 
      d180 = 180.D+00 
      d360 = 360.D+00 
c 
c     define rad for conversion to radians. 
      rad = acos(-d1)/d180 
c 
c     convert latitude and longitude to spherical coordinates 
      thetas = gln 
      if (thetas .gt. d180) thetas = thetas - d360 
      if (thetas .lt. -d180) thetas = thetas + d360 
      thetas = thetas*rad 
      phis = (d90 - glt)*rad 
c 
c     translate point into Cartesian coordinates for rotation. 
      gx = sin(phis)*cos(thetas) 
      gy = sin(phis)*sin(thetas) 
      gz = cos(phis) 
c 
c     rotate the point (gx, gy, gz) about the z axis by phir then the x 
c     axis by thetar and finally about the z axis by psir. 
c 
      rx = gx*(cos(psir)*cos(phir) - cos(thetar)*sin(phir)*sin(psir)) + 
     &     gy*(cos(psir)*sin(phir) + cos(thetar)*cos(phir)*sin(psir)) + 
     &     gz*sin(psir)*sin(thetar) 
c 
      ry = gx*(-sin(psir)*cos(phir) - cos(thetar)*sin(phir)*cos(psir)) + 
     &     gy*(-sin(psir)*sin(phir) + cos(thetar)*cos(phir)*cos(psir)) + 
     &     gz*(cos(psir)*sin(thetar)) 
c 
      rz = gx*(sin(thetar)*sin(phir)) + gy*(-sin(thetar)*cos(phir)) + 
     &     gz*(cos(thetar)) 
c 
c     convert rotated point back to spherical coordinates 
c 
c     check for rounding error (arccos(x): abs(x) must be .le. 1) 
      rz = min(rz, d1) 
      rz = max(rz, -d1) 
      rlt = acos(rz) 
c     if point is at a pole set rotated longitude equal to initial. 
      if (rlt .le. d0 .or. rlt .ge. d180*rad) then 
        rln = thetas 
      else 
c       if rln lies between -135 and -45 or between 45 and 135 degrees 
c       it is more accurate to use an arccos calculation. 
        if (abs(rx/sin(rlt)) .lt. cos(45.D+00*rad)) then 
          rln = acos(max(min(rx/sin(rlt), d1), -d1)) 
c         arccos will give rln between 0 and 180 degrees.  if the point 
c         is negative in y, rln must be equal to negative rln. 
          if (ry .lt. d0) rln = -rln 
        else 
c         if rln lies between -45 and 45 or between 135 and -135 degrees 
c         it is more accurate to use an arcsin calculation. 
          rln = asin(max(min(ry/sin(rlt), d1), -d1)) 
c         arcsin will give rln between -90 and 90 degrees. if the point 
c         is negative in x, rln must be equal to 180 degrees minus rln. 
          if (rx .lt. d0) rln = d180*rad - rln 
        endif 
      endif 
c 
c     convert back to degrees of latitude and longitude. 
      rlt = d90 - rlt/rad 
      rln = rln/rad 
      if (rln .gt. d180) rln = rln - d360 
      if (rln .le. -d180) rln = rln + d360 

      return 
      end