C $Header: /u/gcmpack/MITgcm/model/src/rotate_spherical_polar_grid.F,v 1.4 2011/12/25 22:24:35 jmc Exp $
C $Name:  $

#include "CPP_OPTIONS.h"

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
CBOP
C     !ROUTINE: ROTATE_SPHERICAL_POLAR_GRID
C     !INTERFACE:
      SUBROUTINE ROTATE_SPHERICAL_POLAR_GRID(
     U                  X, Y,
     I                  myThid )

C     !DESCRIPTION: \bv
C     *===================================================================*
C     | SUBROUTINE ROTATE_SPHERICAL_POLAR_GRID
C     | o rotate the model coordinates on the input arrays to
C     |   geographical coordinates
C     | o this is useful when a rotated spherical grid is used,
C     |   e.g., in order to avoid the pole singularity.
C     | The three Euler angles PhiEuler, ThetaEuler, and PsiEuler
C     | define the rotation about the original z-axis (of an sphere
C     | centered cartesian grid), the new x-axis, and the new z-axis,
C     | respectively.
C     | The input coordinates X, Y are assumed to be the model coordinates
C     | on a rotated grid defined by the Euler angles. In this S/R they
C     | are rotated *BACK* to the geographical coordinate; that is why
C     | all rotation matrices are the inverses of the original matrices.
C     | On exit X and Y are the geographical coordinates, that are
C     | used to compute the Coriolis parameter and also to interpolate
C     | forcing fields to as in pkg/exf/exf_interf.F
C     | Naturally, this feature does not work with all packages, so the
C     | some combinations are prohibited in config_summary (flt,
C     | flt_zonal, ecco, profiles), because there the coordinates are
C     | assumed to be regular spherical grid coordinates.
C     *===================================================================*
C     \ev

C     !USES:
      IMPLICIT NONE
C     === Global variables ===
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"

C     !INPUT/OUTPUT PARAMETERS:
C     == Routine arguments ==
C     X, Y   :: on entry: model coordinate location
C            :: on exit: geographical coordinate location
C     myThid :: my Thread Id Number
      _RS X(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RS Y(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      INTEGER myThid
CEOP

C     !LOCAL VARIABLES:
C     == Local variables ==
C     bi,bj  :: Tile indices
C     i, j   :: Loop counters
      INTEGER bi, bj
      INTEGER  I,  J, iA, jA, kA
C     Euler angles in radians
      _RL phiRad, thetaRad, psiRad
C     inverted rotation matrix
      _RL Ainv(3,3), Binv(3,3), Cinv(3,3), Dinv(3,3), CB(3,3)
C     cartesian coordinates
      _RL XYZgeo(3), XYZrot(3)
C     some auxilliary variables
      _RL hypotxy
CEOP

C     convert to radians
      phiRad    = phiEuler  *deg2rad
      thetaRad  = thetaEuler*deg2rad
      psiRad    = psiEuler  *deg2rad

C     create inverse of full rotation matrix
      Dinv(1,1) =  COS(phiRad)
      Dinv(1,2) = -SIN(phiRad)
      Dinv(1,3) =  0. _d 0
C
      Dinv(2,1) =  SIN(phiRad)
      Dinv(2,2) =  COS(phiRad)
      Dinv(2,3) =  0. _d 0
C
      Dinv(3,1) =  0. _d 0
      Dinv(3,2) =  0. _d 0
      Dinv(3,3) =  1. _d 0
C
      Cinv(1,1) =  1. _d 0
      Cinv(1,2) =  0. _d 0
      Cinv(1,3) =  0. _d 0
C
      Cinv(2,1) =  0. _d 0
      Cinv(2,2) =  COS(thetaRad)
      Cinv(2,3) = -SIN(thetaRad)
C
      Cinv(3,1) =  0. _d 0
      Cinv(3,2) =  SIN(thetaRad)
      Cinv(3,3) =  COS(thetaRad)
C
      Binv(1,1) =  COS(psiRad)
      Binv(1,2) = -SIN(psiRad)
      Binv(1,3) =  0. _d 0
C
      Binv(2,1) =  SIN(psiRad)
      Binv(2,2) =  COS(psiRad)
      Binv(2,3) =  0. _d 0
C
      Binv(3,1) =  0. _d 0
      Binv(3,2) =  0. _d 0
      Binv(3,3) =  1. _d 0
C     Ainv = Binv*Cinv*Dinv (matrix multiplications)
      DO jA=1,3
       DO iA=1,3
        Ainv(iA,jA) = 0. _d 0
        CB  (iA,jA) = 0. _d 0
       ENDDO
      ENDDO
      DO jA=1,3
       DO iA=1,3
        DO kA=1,3
         CB  (iA,jA) = CB  (iA,jA) + Cinv(iA,kA)*Binv(kA,jA)
        ENDDO
       ENDDO
      ENDDO
      DO jA=1,3
       DO iA=1,3
        DO kA=1,3
         Ainv(iA,jA) = Ainv(iA,jA) + Dinv(iA,kA)*CB(kA,jA)
        ENDDO
       ENDDO
      ENDDO
C

C     For each tile ...
      DO bj = myByLo(myThid), myByHi(myThid)
       DO bi = myBxLo(myThid), myBxHi(myThid)

C     loop over grid points
        DO J = 1-OLy,sNy+OLy
         DO I = 1-OLx,sNx+OLx
C     transform spherical coordinates with unit radius
C     to sphere centered cartesian coordinates
          XYZrot(1) =
     &         COS( Y(I,J,bi,bj)*deg2rad )*COS( X(I,J,bi,bj)*deg2rad )
          XYZrot(2) =
     &         COS( Y(I,J,bi,bj)*deg2rad )*SIN( X(I,J,bi,bj)*deg2rad )
          XYZrot(3) = SIN( Y(I,J,bi,bj)*deg2rad )
C     rotate cartesian coordinate (matrix multiplication)
          DO iA=1,3
           XYZgeo(iA) = 0. _d 0
          ENDDO
          DO iA=1,3
           DO jA=1,3
            XYZgeo(iA) = XYZgeo(iA) + Ainv(iA,jA)*XYZrot(jA)
           ENDDO
          ENDDO
C     tranform cartesian coordinates back to spherical coordinates
          hypotxy = SQRT( ABS(XYZgeo(1))*ABS(XYZgeo(1))
     &                  + ABS(XYZgeo(2))*ABS(XYZgeo(2)) )
          IF(XYZgeo(1) .EQ. 0. _d 0 .AND. XYZgeo(2) .EQ. 0. _d 0)THEN
C     happens exactly at the poles
           X(I,J,bi,bj) = 0. _d 0
          ELSE
           X(I,J,bi,bj) = ATAN2(XYZgeo(2),XYZgeo(1))/deg2rad
           IF ( X(I,J,bi,bj) .LT. 0. _d 0 )
     &          X(I,J,bi,bj) = X(I,J,bi,bj) + 360. _d 0
          ENDIF
          IF(hypotxy .EQ. 0. _d 0 .AND. XYZgeo(3) .EQ. 0. _d 0)THEN
C     this can not happen for a sphere with unit radius
           Y(I,J,bi,bj) = 0. _d 0
          ELSE
           Y(I,J,bi,bj) = ATAN2(XYZgeo(3),hypotxy)/deg2rad
          ENDIF
         ENDDO
        ENDDO
C     bi,bj-loops
       ENDDO
      ENDDO

      RETURN
      END