model/src/rotate_spherical_polar_grid.F

7514c1bd55 Mart*0001 #include "CPP_OPTIONS.h"
                0002
4761358c97 Jean*0003 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
7514c1bd55 Mart*0004 CBOP
                0005 C     !ROUTINE: ROTATE_SPHERICAL_POLAR_GRID
                0006 C     !INTERFACE:
4761358c97 Jean*0007       SUBROUTINE ROTATE_SPHERICAL_POLAR_GRID(
                0008      U                  X, Y,
                0009      I                  myThid )
                0010
7514c1bd55 Mart*0011 C     !DESCRIPTION: \bv
4761358c97 Jean*0012 C     *===================================================================*
                0013 C     | SUBROUTINE ROTATE_SPHERICAL_POLAR_GRID
                0014 C     | o rotate the model coordinates on the input arrays to
                0015 C     |   geographical coordinates
                0016 C     | o this is useful when a rotated spherical grid is used,
                0017 C     |   e.g., in order to avoid the pole singularity.
                0018 C     | The three Euler angles PhiEuler, ThetaEuler, and PsiEuler
                0019 C     | define the rotation about the original z-axis (of an sphere
                0020 C     | centered cartesian grid), the new x-axis, and the new z-axis,
                0021 C     | respectively.
                0022 C     | The input coordinates X, Y are assumed to be the model coordinates
                0023 C     | on a rotated grid defined by the Euler angles. In this S/R they
                0024 C     | are rotated *BACK* to the geographical coordinate; that is why
                0025 C     | all rotation matrices are the inverses of the original matrices.
                0026 C     | On exit X and Y are the geographical coordinates, that are
                0027 C     | used to compute the Coriolis parameter and also to interpolate
                0028 C     | forcing fields to as in pkg/exf/exf_interf.F
                0029 C     | Naturally, this feature does not work with all packages, so the
                0030 C     | some combinations are prohibited in config_summary (flt,
                0031 C     | flt_zonal, ecco, profiles), because there the coordinates are
                0032 C     | assumed to be regular spherical grid coordinates.
                0033 C     *===================================================================*
7514c1bd55 Mart*0034 C     \ev
                0035
                0036 C     !USES:
                0037       IMPLICIT NONE
                0038 C     === Global variables ===
                0039 #include "SIZE.h"
                0040 #include "EEPARAMS.h"
                0041 #include "PARAMS.h"
                0042
                0043 C     !INPUT/OUTPUT PARAMETERS:
                0044 C     == Routine arguments ==
4761358c97 Jean*0045 C     X, Y   :: on entry: model coordinate location
                0046 C            :: on exit: geographical coordinate location
                0047 C     myThid :: my Thread Id Number
c679ff8f37 Jean*0048       _RS X(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0049       _RS Y(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
4761358c97 Jean*0050       INTEGER myThid
                0051 CEOP
7514c1bd55 Mart*0052
                0053 C     !LOCAL VARIABLES:
                0054 C     == Local variables ==
4761358c97 Jean*0055 C     bi,bj  :: Tile indices
                0056 C     i, j   :: Loop counters
7514c1bd55 Mart*0057       INTEGER bi, bj
                0058       INTEGER  I,  J, iA, jA, kA
                0059 C     Euler angles in radians
                0060       _RL phiRad, thetaRad, psiRad
                0061 C     inverted rotation matrix
                0062       _RL Ainv(3,3), Binv(3,3), Cinv(3,3), Dinv(3,3), CB(3,3)
                0063 C     cartesian coordinates
                0064       _RL XYZgeo(3), XYZrot(3)
                0065 C     some auxilliary variables
66314ec55a Mart*0066       _RL hypotxy
7514c1bd55 Mart*0067 CEOP
                0068
                0069 C     convert to radians
                0070       phiRad    = phiEuler  *deg2rad
                0071       thetaRad  = thetaEuler*deg2rad
                0072       psiRad    = psiEuler  *deg2rad
                0073
4761358c97 Jean*0074 C     create inverse of full rotation matrix
7514c1bd55 Mart*0075       Dinv(1,1) =  COS(phiRad)
                0076       Dinv(1,2) = -SIN(phiRad)
                0077       Dinv(1,3) =  0. _d 0
                0078 C
                0079       Dinv(2,1) =  SIN(phiRad)
                0080       Dinv(2,2) =  COS(phiRad)
                0081       Dinv(2,3) =  0. _d 0
                0082 C
                0083       Dinv(3,1) =  0. _d 0
                0084       Dinv(3,2) =  0. _d 0
                0085       Dinv(3,3) =  1. _d 0
                0086 C
                0087       Cinv(1,1) =  1. _d 0
                0088       Cinv(1,2) =  0. _d 0
                0089       Cinv(1,3) =  0. _d 0
                0090 C
                0091       Cinv(2,1) =  0. _d 0
                0092       Cinv(2,2) =  COS(thetaRad)
                0093       Cinv(2,3) = -SIN(thetaRad)
                0094 C
                0095       Cinv(3,1) =  0. _d 0
                0096       Cinv(3,2) =  SIN(thetaRad)
                0097       Cinv(3,3) =  COS(thetaRad)
                0098 C
                0099       Binv(1,1) =  COS(psiRad)
                0100       Binv(1,2) = -SIN(psiRad)
                0101       Binv(1,3) =  0. _d 0
                0102 C
                0103       Binv(2,1) =  SIN(psiRad)
                0104       Binv(2,2) =  COS(psiRad)
                0105       Binv(2,3) =  0. _d 0
                0106 C
                0107       Binv(3,1) =  0. _d 0
                0108       Binv(3,2) =  0. _d 0
                0109       Binv(3,3) =  1. _d 0
                0110 C     Ainv = Binv*Cinv*Dinv (matrix multiplications)
4761358c97 Jean*0111       DO jA=1,3
                0112        DO iA=1,3
                0113         Ainv(iA,jA) = 0. _d 0
                0114         CB  (iA,jA) = 0. _d 0
7514c1bd55 Mart*0115        ENDDO
                0116       ENDDO
4761358c97 Jean*0117       DO jA=1,3
                0118        DO iA=1,3
                0119         DO kA=1,3
                0120          CB  (iA,jA) = CB  (iA,jA) + Cinv(iA,kA)*Binv(kA,jA)
7514c1bd55 Mart*0121         ENDDO
                0122        ENDDO
                0123       ENDDO
4761358c97 Jean*0124       DO jA=1,3
                0125        DO iA=1,3
                0126         DO kA=1,3
                0127          Ainv(iA,jA) = Ainv(iA,jA) + Dinv(iA,kA)*CB(kA,jA)
7514c1bd55 Mart*0128         ENDDO
                0129        ENDDO
                0130       ENDDO
                0131 C
                0132
                0133 C     For each tile ...
                0134       DO bj = myByLo(myThid), myByHi(myThid)
                0135        DO bi = myBxLo(myThid), myBxHi(myThid)
                0136
                0137 C     loop over grid points
c679ff8f37 Jean*0138         DO J = 1-OLy,sNy+OLy
                0139          DO I = 1-OLx,sNx+OLx
4761358c97 Jean*0140 C     transform spherical coordinates with unit radius
7514c1bd55 Mart*0141 C     to sphere centered cartesian coordinates
4761358c97 Jean*0142           XYZrot(1) =
7514c1bd55 Mart*0143      &         COS( Y(I,J,bi,bj)*deg2rad )*COS( X(I,J,bi,bj)*deg2rad )
4761358c97 Jean*0144           XYZrot(2) =
7514c1bd55 Mart*0145      &         COS( Y(I,J,bi,bj)*deg2rad )*SIN( X(I,J,bi,bj)*deg2rad )
                0146           XYZrot(3) = SIN( Y(I,J,bi,bj)*deg2rad )
                0147 C     rotate cartesian coordinate (matrix multiplication)
                0148           DO iA=1,3
                0149            XYZgeo(iA) = 0. _d 0
                0150           ENDDO
                0151           DO iA=1,3
                0152            DO jA=1,3
                0153             XYZgeo(iA) = XYZgeo(iA) + Ainv(iA,jA)*XYZrot(jA)
                0154            ENDDO
                0155           ENDDO
4761358c97 Jean*0156 C     tranform cartesian coordinates back to spherical coordinates
7514c1bd55 Mart*0157           hypotxy = SQRT( ABS(XYZgeo(1))*ABS(XYZgeo(1))
                0158      &                  + ABS(XYZgeo(2))*ABS(XYZgeo(2)) )
                0159           IF(XYZgeo(1) .EQ. 0. _d 0 .AND. XYZgeo(2) .EQ. 0. _d 0)THEN
                0160 C     happens exactly at the poles
                0161            X(I,J,bi,bj) = 0. _d 0
                0162           ELSE
                0163            X(I,J,bi,bj) = ATAN2(XYZgeo(2),XYZgeo(1))/deg2rad
4761358c97 Jean*0164            IF ( X(I,J,bi,bj) .LT. 0. _d 0 )
7514c1bd55 Mart*0165      &          X(I,J,bi,bj) = X(I,J,bi,bj) + 360. _d 0
                0166           ENDIF
                0167           IF(hypotxy .EQ. 0. _d 0 .AND. XYZgeo(3) .EQ. 0. _d 0)THEN
                0168 C     this can not happen for a sphere with unit radius
                0169            Y(I,J,bi,bj) = 0. _d 0
                0170           ELSE
                0171            Y(I,J,bi,bj) = ATAN2(XYZgeo(3),hypotxy)/deg2rad
                0172           ENDIF
                0173          ENDDO
                0174         ENDDO
                0175 C     bi,bj-loops
                0176        ENDDO
                0177       ENDDO
                0178
                0179       RETURN
                0180       END