SUBROUTINE LMLLIJ(KFILDO,ALAT,ALON,XMESHL,ORIENT,XLAT,
     1                  XLATLL,XLONLL,XI,YJ)
C
C        FEBRUARY 1995   GLAHN   TDL   MOS-2000
C        JULY     2002   GLAHN   ADDED DIAGNOSTIC PRINT 100
C        AUGUST   2015   GLAHN   ADDED DIAGNOSTIC PRINT 101
C
C        PURPOSE
C            TO CONVERT FROM LATITUDE, LONGITUDE TO GRID COORDINATES ON
C            A LAMBERT CONFORMAL PROJECTION WITH THE LATITUDE OF TANGENCY
C            IN THE NORTHERN HEMISPHERE.  ADAPTED FROM NMC'S W3FB11, WITH
C            THE OPTION FOR TANGENCY IN THE SOUTHERN HEMISPHERE REMOVED 
C            AND THE INPUT CHANGED FROM EAST TO WEST LONGITUDE.  ALTHOUGH
C            THE INPUT IS IN W LONGITUDE, THE COMPUTATIONS, BEING BASED
C            ON THOSE IN W3FB11, ARE IN TERMS OF E LONGITUDE.
C            FORMULAE AND NOTATION LOOSELY BASED ON HOKE, HAYES,
C            AND RENNINGER'S "MAP PROJECTIONS AND GRID SYSTEMS...",
C            MARCH 1981, AFGWC/TN-79/003.  ORIGINAL AUTHOR, STACKPOLE.
C
C        DATA SET USE
C            KFILDO - UNIT NUMBER OF OUTPUT (PRINT) FILE.  (OUTPUT)
C
C        VARIABLES
C            INPUT
C              KFILDO = UNIT NUMBER OF OUTPUT (PRINT) FILE.  (INPUT)
C                ALAT = NORTH LATITUDE IN DEGREES FOR WHICH THE GRID
C                       COORDINATES ARE WANTED.  NEGATIVE FOR 
C                       SOUTHERN HEMSIPHERE.
C                ALON = WEST LONGITUDE IN DEGREES FOR WHICH THE GRID
C                       COORDINATES ARE WANTED.  DON'T USE NEGATIVE. 
C              XMESHL = MESH LENGTH IN METERS AT XLAT DEGREES N LATITUDE.
C              ORIENT = ORIENTATION IN DEGREES WEST LONGITUDE.
C                XLAT = LATITUDE IN DEGREES AT WHICH XMESHL APPLIES.
C                       ALSO THE LATITUDE OF TANGENCY.  DON'T USE
C                       NEGATIVE. 
C              XLATLL = LATITUDE OF LOWER LEFT (1,1) CORNER POINT OF
C                       THE GRID.
C              XLONLL = LONGITUDE OF LOWER LEFT (1,1) CORNER POINT OF
C                       THE GRID.  DON'T USE NEGATIVE. 
C
C            OUTPUT
C               XI,YJ = IJ (LEFT TO RIGHT) AND JY (BOTTON TO TOP)
C                       GRIDPOINT NUMBERS OF THE POINT ALAT, ALON,
C                       CONSIDERING THE LOWER LEFT CORNER POINT TO BE (1,1).
C
C            INTERNAL
C              RADPDG = NUMBER OF RADIANS PER DEGREE.  SET BY PARAMETER.
C                  PI = PI.  SET BY PARAMETER.
C               RERTH = RADIUS OF THE EARTH IN METERS.  SET BY PARAMETER.
C                ELON = EAST LONGITUDE IN DEGREES.
C               ELON1 = E LONGITUDE OF LOWER LEFT (1,1) POINT OF GRID.
C               ELONV = THE ORIENTATION OF THE GRID.  I.E.,
C                       THE EAST LONGITUDE VALUE OF THE VERTICAL MERIDIAN
C                       WHICH IS PARALLEL TO THE Y-AXIS (OR COLUMNS OF
C                       OF THE GRID) ALON WHICH LATITUDE INCREASES AS
C                       THE Y-COORDINATE INCREASES.
C                       THIS IS ALSO THE MERIDIAN (ON THE BACK SIDE OF THE
C                       TANGENT CONE) ALON WHICH THE CUT IS MADE TO LAY
C                       THE CONE FLAT.
C
C        NONSYSTEM ROUTINES CALLED
C            NONE 
C
      PARAMETER (PI=3.14159,
     1           RERTH=6371200.,
     2           RADPDG=PI/180.)
C
      IF(XMESHL.LE.0..OR.ORIENT.LT.0..OR.XLAT.LT.0..OR.XLONLL.LT.0.)THEN
         WRITE(KFILDO,100)XMESHL,ORIENT,XLAT,XLONLL
 100     FORMAT(/' ****PROBLEM WITH EITHER'/
     1           '     XMESHL =',F12.4,','/
     2           '     ORIENT =',F12.4,','/
     3           '     XLAT   =',F12.4,', OR'/
     4           '     XLONLL =',F12.4,'.'/
     5           '     STOP IN LMLLIJ AT 100.')
         WRITE(KFILDO,101)XLATLL,ALAT,ALON
 101     FORMAT( '     XLATLL =',F12.4,','/
     1           '     ALAT   =',F12.4,','/
     2           '     ALON   =',F12.4,',')
         STOP 100
      ENDIF
C
      ELON=360.-ALON
      ELON1=360.-XLONLL
      ELONV=360.-ORIENT
C
C        PRELIMINARY VARIABLES AND REDIFINITIONS
C
      REBYXMESHL = RERTH/XMESHL
      ALATN1 = XLAT * RADPDG
      AN = SIN(ALATN1)
      COSLTN = COS(ALATN1)
C
C        MAKE SURE THAT INPUT LONGITUDES DO NOT PASS THROUGH
C        THE CUT ZONE (FORBIDDEN TERRITORY) OF THE FLAT MAP
C        AS MEASURED FROM THE VERTICAL (REFERENCE) LONGITUDE.
C
      ELON1L = ELON1
      IF((ELON1 - ELONV).GT. 180.)ELON1L = ELON1 - 360.
      IF((ELON1 - ELONV).LT.-180.)ELON1L = ELON1 + 360.
C
      ELONL = ELON
      IF((ELON - ELONV).GT. 180.)ELONL = ELON - 360.
      IF((ELON - ELONV).LT.-180.)ELONL = ELON + 360.
C
      ELONVR = ELONV *RADPDG
C
C        RADIUS TO LOWER LEFT HAND (LL) CORNER
C
      ALA1 =  XLATLL * RADPDG
      RMLL = REBYXMESHL * (((COSLTN)**(1.-AN))*(1.+AN)**AN) *
     1        (((COS(ALA1))/(1.+SIN(ALA1)))**AN)/AN
C
C        USE LL POINT INFO TO LOCATE POLE POINT.   
C
      ELO1 = ELON1L * RADPDG
      ARG = AN * (ELO1-ELONVR)
      POLEI = 1. - RMLL * SIN(ARG)
      POLEJ = 1. + RMLL * COS(ARG)
C
C        RADIUS TO DESIRED POINT AND THE I J TOO.
C
      ALA =  ALAT * RADPDG
      RM = REBYXMESHL * ((COSLTN**(1.-AN))*(1.+AN)**AN) *
     1        (((COS(ALA))/(1.+SIN(ALA)))**AN)/AN
C
      ELO = ELONL * RADPDG
      ARG = AN*(ELO-ELONVR)
      XI = POLEI + RM * SIN(ARG)
      YJ = POLEJ - RM * COS(ARG)
C
C        IF COORDINATE LESS THAN 1
C        COMPENSATE FOR ORIGIN AT (1,1)
C
C***      IF(XI.LT.1.)  XI = XI - 1.
C***      IF(YJ.LT.1.)  YJ = YJ - 1.
C***     ABOVE 2 STATEMENTS REMOVED 9/26/92 BY GLAHN
C***     THESE STATEMENTS WERE IN W3FB11.
D     WRITE(KFILDO,900)ALAT,ALON,XMESHL,ORIENT,XLAT,XLATLL,XLONLL,XI,YJ
D900  FORMAT(' LMLLIJ--ALAT,ALON,XMESHL,ORIENT,XLAT,XLATLL,XLONLL,',
D    1       'XI,YJ',9F10.2)
C
      RETURN
      END