!                         WD23JP
!                                                         *****************
!                                                         *   SOCT4.FOR   *
!                                                         *  PURSER 1996  *
!                                                         *****************
!
!------------------------------------------------------------------------------C
!   R.J.Purser, NCEP, Washington D.C.   1996                                   C
!                   SUBROUTINE OTOC1                                           C
!   Transform to earth-centered cartesians from octagon-1 (OTOC1) or from      C
!   octagon-2 (OTOC2).                                                         C
!                                                                              C
! --> XM   x and y map-coordinate of point in octagon-1 (OTOC1) or             C
!           octagon-2 (OTOC2)                                                  C
! <-- XE    3-vector of earth-centered cartesian coordinates corresponding     C
!           to map location (X,Y).                                             C
!------------------------------------------------------------------------------C
      SUBROUTINE OTOC1(XM,XE)
		use proc1
      COMMON/CSTOCT/MAG,M1G,ROTM(3,3),SC,SCI,A
      DIMENSION XE(3),XM(2),DXCDXM(3,2)
      CALL XMTOXC(XM,XE,DXCDXM,1)
      CALL SCHMIDT(XE,XE,SC)
      CALL NMAPT(ROTM,XE,XE)
      RETURN

      ENTRY OTOC2(XM,XE)
      CALL XMTOXC(XM,XE,DXCDXM,2)
      CALL SCHMIDT(XE,XE,SC)
      CALL NMAPT(ROTM,XE,XE)
      RETURN
      END

!------------------------------------------------------------------------------C
!   R.J.Purser, National Meteorological Center, Washington D.C.  1995          C
!                   SUBROUTINE CTOO                                            C
!   Transform from earth-centered cartesians to octagon-1 (KMAP=1) or to       C
!   octagon-2 (KMAP=-1). The method uses Newton iterations                     C
!                                                                              C
! --> XC    earth-centered coordinates (3-vector) of a point                   C
! <-- XM   map-coordinates of this point                                       C
! <-- KMAP  map indicator (=1 for octagon-1, =-1 for octagon-2)                C
!------------------------------------------------------------------------------C
      SUBROUTINE CTOO(XC,XM,KMAP)
		use proc1
		use mulmm_m
      COMMON/CSTOCT/MAG,M1G,ROTM(3,3),SC,SCI,A
      DIMENSION XC(3),XM(2),XCT(3),DXCDXM(3,2),DXMDXC(2,3),G2(2,2)
      DATA NIT/10/,SMAX/1.E-10/
!  RE-ORIENT THE CARTESIAN POSITION VECTOR TO THE FRAME IN WHICH OCTAGON
!  CENTERS ARE AT (0,0,1) AND (0,0,-1):
      CALL NMAP(ROTM,XC,XC)
!  PERFORM THE SCHMIDT TRANSFORMATION:
      CALL SCHMIDT(XC,XC,SCI)
      IF(XC(3).GT.0.)THEN
       KMAP=1
      ELSE
       KMAP=-1
       XC(1)=-XC(1)
       XC(3)=-XC(3)
      ENDIF
      XM(1)=XC(1)
      XM(2)=XC(2)
      CALL XMTOXC(XM,XCT,DXCDXM,1)
      DO IT=1,NIT
       CALL MULTM(DXCDXM,DXCDXM,G2,2,3,2,3,3,2)
       DO J=1,3
       DO I=1,2
        DXMDXC(I,J)=DXCDXM(J,I)
       ENDDO
       ENDDO
       DET=G2(1,1)*G2(2,2)-G2(1,2)*G2(2,1)
       IF(ABS(DET).LE.1.E-14)RETURN
       CALL LINVMM(G2,DXMDXC,2,3,2,2)
       S=0.
       DO I=1,3
        XCT(I)=XCT(I)-XC(I)
        S=S+XCT(I)**2
       ENDDO
       CALL MSBMM(DXMDXC,XCT,XM,2,3,1,2,3,2)
       IF(S.LE.SMAX)RETURN
       CALL XMTOXC(XM,XCT,DXCDXM,1)
      ENDDO
      PRINT'('' ITERATIONS IN XCTOXM INSUFFICIENT FOR CONVERGENCE'')'
      PRINT'('' ARE YOU SURE THE INPUT VECTOR IS VALID?'')'
      PRINT'('' OR YOU MIGHT LOOSEN THE CONVERGENCE CRITERION SMAX'')'
      PRINT'('' XC='',3(1X,E12.6))',(XC(I),I=1,3)
      PRINT'('' S='',E12.6)',S
      PRINT'('' XM='',2(1X,E12.6))',(XM(I),I=1,2)
      PRINT'('' XCT='',3(1X,E12.6))',(XCT(I),I=1,3)
      STOP
      END

!------------------------------------------------------------------------------C
!   R.J.Purser, National Meteorological Center, Washington D.C.  1994          C
!                   SUBROUTINE STOE                                            C
!   Transform latitude and longitude (degrees) to earth-centered cartesian     C
!   coordinates.                                                               C
!  --> DLAT     latitude                                                       C
!  --> DLON     longitude                                                      C
!  <-- XE       three cartesian components.                                    C
!------------------------------------------------------------------------------C
      SUBROUTINE STOE(DLAT,DLON,XE)
      DIMENSION XE(3)
      DATA DTOR/1.745329251994E-2/
      RLAT=DTOR*DLAT
      RLON=DTOR*DLON
      SLA=SIN(RLAT)
      CLA=COS(RLAT)
      SLO=SIN(RLON)
      CLO=COS(RLON)
      XE(1)=CLA*CLO
      XE(2)=CLA*SLO
      XE(3)=SLA
      RETURN
      END

      SUBROUTINE SCHMIDT(XC1,XC2,S)
!  EVALUATE BASIC SCHMIDT TRANSFORMATION
      DIMENSION XC1(3),XC2(3)
      SS=S*S
      SSP=1.+SS
      SSC=1.-SS
      S2=S*2
      DI=1./(SSP+SSC*XC1(3))
      DIS2=DI*S2
      XC2(1)=DIS2*XC1(1)
      XC2(2)=DIS2*XC1(2)
      XC2(3)=DI*(SSC+SSP*XC1(3))
      RETURN
      END

      SUBROUTINE DSCHMIDT(XC1,XC2,D2D1,S)
!  EVALUATE BASIC SCHMIDT TRANSFORMATION TOGETHER WITH ITS DERIVATIVE
      DIMENSION XC1(3),XC2(3),D2D1(3,3),T(3)
! COPY XC1 TO T AS A PRECAUTION IN CASE XC1 AND XC2 OCCUPY SAME MEMORY:
      DO I=1,3
       T(I)=XC1(I)
      ENDDO
      SS=S*S
      SSP=1.+SS
      SSC=1.-SS
      S2=S*2
      DI=1./(SSP+SSC*XC1(3))
      DIS2=DI*S2
      SSCDI=SSC*DI
      XC2(1)=DIS2*XC1(1)
      XC2(2)=DIS2*XC1(2)
      XC2(3)=DI*(SSC+SSP*XC1(3))
      D2D1(1,1)=DIS2
      D2D1(2,1)=0.
      D2D1(3,1)=0.
      D2D1(1,2)=0.
      D2D1(2,2)=DIS2
      D2D1(3,2)=0.
      D2D1(1,3)=-XC2(1)*SSCDI
      D2D1(2,3)=-XC2(2)*SSCDI
      D2D1(3,3)=-XC2(3)*SSCDI+SSP*DI
      DO I=1,3
       V=XC2(I)
       DO J=1,3
        V=V-D2D1(I,J)*T(J)
       ENDDO
       DO J=1,3
        D2D1(I,J)=D2D1(I,J)+V*T(J)
       ENDDO
      ENDDO
      RETURN
      END