*-----------------------------------------------------------------------
*     MAP PROJECTION (LAMBERT CONICAL 2)               93/02/20 S.SAKAI
*-----------------------------------------------------------------------
*     Copyright (C) 2000-2016 GFD Dennou Club. All rights reserved.
*-----------------------------------------------------------------------
      SUBROUTINE MPFCOC(XLON, YLAT, X, Y)

      PARAMETER (EPSL = 1.E-5)

      EXTERNAL  RFPI, XMPLON

      SAVE


      TH  = PI/2 - S*YLAT
      IF (TH .GE. PI-EPSL) TH = PI - EPSL

      TT  = TAN(TH/2.)
      IF (ABS(TT) .LT. 1.E-6) TT = 0.
      R   = CK*TT**DK
      IF (R .GT. 10.) R = 10
      DLM = XMPLON(XLON)*DK

      X =    R*SIN(DLM)
      Y = -S*R*COS(DLM)

      RETURN
*-----------------------------------------------------------------------
      ENTRY MPICOC(X, Y, XLON, YLAT)

      R = SQRT(X*X + Y*Y)
      IF (R.EQ.0.) THEN
        XLON = 0.
        YLAT = PI/2.
        RETURN
      ELSE
        XLON = ATAN2(X, -S*Y)/DK
        IF (ABS(XLON) .LE. PI) THEN
          YLAT = S*(PI/2 - 2*ATAN((R/CK)**(1./DK)))
          RETURN
        END IF
      END IF
      XLON = RNA
      YLAT = RNA

      RETURN
*-----------------------------------------------------------------------
      ENTRY MPSCOC(YLAT1, YLAT2)

      PI = RFPI()
      CALL GLRGET('RUNDEF',RNA)

      S = SIGN(1., YLAT1+YLAT2)

      TH1 = PI/2 - MAX(ABS(YLAT1), ABS(YLAT2))
      TH2 = PI/2 - MIN(ABS(YLAT1), ABS(YLAT2))

      DK = (LOG(SIN(TH2  )) - LOG(SIN(TH1  )))
     #    /(LOG(TAN(TH2/2)) - LOG(TAN(TH1/2)))
      CK = SIN(TH1)/DK/TAN(TH1/2)**DK

      END
