*-----------------------------------------------------------------------
*     COORDINATE ROTATION (3D SPHERICAL)
*                                              93/02/19   S.SAKAI
*                                              93/10/29   S.SAKAI
*                                            2000/04/28   S.SAKAI
*                                            2002/01/17   S.SAKAI
*-----------------------------------------------------------------------
*     Copyright (C) 2000-2016 GFD Dennou Club. All rights reserved.
*-----------------------------------------------------------------------
      SUBROUTINE CR3S(THETA, PHI, PSI, THETA0, PHI0, THETA1, PHI1)

      LOGICAL   LREQA, LFIRST

      EXTERNAL  RFPI, LREQA

      SAVE

      DATA      LFIRST /.TRUE./


      IF (LFIRST) THEN
        CALL GLRGET('REPSL',REPSL)
        PI = RFPI()
        LFIRST = .FALSE.
      END IF

      IF (LREQA(THETA, 0., REPSL)) THEN
        THETA1 = THETA0
        PHI1   = PHI0 - PHI - PSI
      ELSE
        CTE = COS(THETA)
        STE = SIN(THETA)

        CT0 = COS(THETA0)
        ST0 = SIN(THETA0)

        CP  = COS(PHI0 - PHI)
        SP  = SIN(PHI0 - PHI)

        A = CT0*CTE + ST0*STE*CP

        IF (ABS(A).LE. 0.8) THEN
          THETA1 = ACOS(A)
        ELSE
          C = SQRT((STE*SP)**2 + (CT0*STE*CP-ST0*CTE)**2)
          THETA1 = ASIN(C)
          IF(A.LE.0) THETA1 = PI - THETA1
        END IF

        A =    SP*ST0
        B = - CT0*STE + ST0*CTE*CP

        IF (A.EQ.0. .AND. B.EQ.0.) THEN
          PHI1 = PHI0 - PHI - PSI
        ELSE
          PHI1 = ATAN2(A, B) - PSI
        END IF
      END IF

      END
