*-----------------------------------------------------------------------
* INV MAP PROJECTION (NEWTON METHOD ON SPHERE)      2007-11-10 E. TOYODA
*-----------------------------------------------------------------------
*     Copyright (C) 2000-2016 GFD Dennou Club. All rights reserved.
*-----------------------------------------------------------------------
      SUBROUTINE MPNWTN(X, Y, XLON, YLAT, FUNC, XCYC)
        LOGICAL XCYC
        PARAMETER (EPSX = 1.0E-4)
        PARAMETER (EPS = 1.0E-4)
        PARAMETER (CONVT = 3.0E-5)
        PARAMETER (XMRGN = 0.01)
        PARAMETER (DETMIN = 1.0E-7)
        PARAMETER (IMAX = 40)
        PARAMETER (PI = ATAN(1.0) * 4.0)
        ACCEL = 1.0
        I = 0
  100   CONTINUE
          IF (XLON .GT. 0.0) THEN
            XLON2 = XLON - EPS
          ELSE
            XLON2 = XLON + EPS
          ENDIF
          IF (YLAT .GT. 0.0) THEN
            YLAT2 = YLAT - EPS
          ELSE
            YLAT2 = YLAT + EPS
          ENDIF
          CALL FUNC(XLON, YLAT, X0, Y0)
          CALL FUNC(XLON2, YLAT, X1, Y1)
          CALL FUNC(XLON, YLAT2, X2, Y2)
C         print "('foot:',4f13.6)", degree(xlon), degree(xlon2),
C    +     degree(ylat), degree(ylat2)
          XL = (X1 - X0) / (XLON2 - XLON)
          XB = (X2 - X0) / (YLAT2 - YLAT)
          YL = (Y1 - Y0) / (XLON2 - XLON)
          YB = (Y2 - Y0) / (YLAT2 - YLAT)
          DX = X - X0
          DY = Y - Y0
          DET = XL * YB - XB * YL
          IF (ABS(DET) .LT. DETMIN) THEN
C           print "('singularity')"
            GOTO 110
          ENDIF
          RDET = ACCEL / DET
C         print "('det=',f15.8)", det
          CX = RDET * (YB * DX - XB * DY)
          CY = RDET * (-YL * DX + XL * DY)
C         print "('corr=',2f13.6)", degree(cx), degree(cy)
          XLON = XLON + CX
          IF (XLON .GT. PI + EPSX) THEN
C           print "(' xlon > 180')"
            IF (XCYC) THEN
              XLON = XLON - 2.0 * PI
            ELSE
              XLON = PI - XMRGN
            ENDIF
            ACCEL = 0.5
          ELSE IF (XLON .LT. -PI - EPSX) THEN
C           print "(' xlon < -180')"
            IF (XCYC) THEN
              XLON = XLON + 2.0 * PI
            ELSE
              XLON = -PI + XMRGN
            ENDIF
            ACCEL = 0.5
          ENDIF
          YLAT = YLAT + CY
          IF (YLAT .GT. 0.5 * PI) THEN
C           print "(' ylat > 90')"
            YLAT = 0.5 * PI - XMRGN
            ACCEL = 0.5
          ELSE IF (YLAT .LT. -0.5 * PI) THEN
C           print "(' ylat < -90')"
            YLAT = -0.5 * PI + XMRGN
            ACCEL = 0.5
          ENDIF
          DIST = MAX(ABS(CX), ABS(CY))
          IF (DIST .LT. CONVT) GOTO 110
          I = I + 1
          IF (I .GT. IMAX) GOTO 911
        GOTO 100
  110   CONTINUE
        RETURN

  911   CONTINUE   
        CALL GLRGET('RUNDEF',RNA)
        XLON = RNA
        YLAT = RNA
      END
