C     $Id: carmel.f 3304 2011-01-17 15:25:59Z brideout $
C
      SUBROUTINE CARMEL(B,XI,VL)
C
C     Private/Internal subroutine. Part of Apex coordinate computation
C     package. See COORD for public API. Computes scaler VL as a
C     function of B and XI.
C
C       Input:
C           B - Scaler field strength value.
C          XI - integral invariant (see INTEG).
C
C      Output:
C          VL - McIlwain's L-shell parameter i.e. Invariant Latitude
C               = ACOS(DSQRT(1.0D0/VL))/DTR
C
C     .. Scalar Arguments ..
      DOUBLE PRECISION B,VL,XI
C     ..
C     .. Local Scalars ..
      DOUBLE PRECISION GG,XX
C     ..
C     .. Intrinsic Functions ..
      INTRINSIC DEXP,DLOG
C     ..
      IF (XI.GT.1.0D-36) THEN
         XX = 3.0D0*DLOG(XI)
         XX = XX + DLOG(B/0.311653D0)
         IF (XX.LE.-22.D0) THEN
            GG = .333338D0*XX + .30062102D0
         ELSE IF (XX.LE.-3.D0) THEN
            GG = ((((((((-8.1537735D-14*XX+8.3232531D-13)*XX+
     *           1.0066362D-9)*XX+8.1048663D-8)*XX+3.2916354D-6)*XX+
     *           8.2711096D-5)*XX+1.3714667D-3)*XX+.015017245D0)*XX+
     *           .43432642D0)*XX + .62337691D0
         ELSE IF (XX.LE.3.D0) THEN
            GG = ((((((((2.6047023D-10*XX+2.3028767D-9)*XX-
     *           2.1997983D-8)*XX-5.3977642D-7)*XX-3.3408822D-6)*XX+
     *           3.8379917D-5)*XX+1.1784234D-3)*XX+1.4492441D-2)*XX+
     *           .43352788D0)*XX + .6228644D0
         ELSE IF (XX.LE.11.7D0) THEN
            GG = ((((((((6.3271665D-10*XX-3.958306D-8)*XX+
     *           9.9766148D-07)*XX-1.2531932D-5)*XX+7.9451313D-5)*XX-
     *           3.2077032D-4)*XX+2.1680398D-3)*XX+1.2817956D-2)*XX+
     *           .43510529D0)*XX + .6222355D0
         ELSE IF (XX.GT.23.D0) THEN
            GG = XX - 3.0460681D0
         ELSE
            GG = (((((2.8212095D-8*XX-3.8049276D-6)*XX+2.170224D-4)*XX-
     *           6.7310339D-3)*XX+.12038224D0)*XX-.18461796D0)*XX +
     *           2.0007187D0
         END IF
         VL = (((1.0D0+DEXP(GG))*0.311653D0)/B)**(1.D0/3.D0)
C        end compute l
      ELSE
         VL = (0.311653D0/B)**(1.D0/3.D0)
      END IF
      RETURN
C
      END