C
C
	SUBROUTINE ILAMBERT1(LAT,LON,X,Y,ORGLAT,ORGLON,IOPT)
C
C	COMPUTES THE INVERSE TRANSFORMATION FROM LAT/LON TO X/Y FOR THE
C	LAMBERT AZIMUTHAL EQUAL-AREA PROJECTION
C
C	SEE "MAP PROJECTIONS USED BY THE U.S. GEOLOGICAL SURVEY"
C	GEOLOGICAL SURVEY BULLETIN 1532, PGS 157-173
C
C	FOR THIS ROUTINE, A SPHERICAL EARTH IS ASSUMED FOR THE PROJECTION.
C	THE ERROR WILL BE SMALL FOR SMALL-SCALE MAPS.  
C	FOR IOPT=1 A FIXED, NOMINAL EARTH RADIUS IS USED.
C	FOR IOPT=2 THE LOCAL RADIUS OF THE EARTH IS USED BASED ON
C	THE 1972 WGS ELLIPSOID MODEL (BULLETIN PG 15).
C
C	WRITTEN BY: DGL 7 FEB 1992
C	MODIFIED:   DGL 4 MAR 1995
C	+ INCLUDED RADIUS OF THE EARTH COMPUTED AT LOCAL REFERENCE POINT
C       REVISED BY: DGL NOV 2005
C        + eliminated COSD and SIND calls with DTR and RTD
C
C	INPUTS:
C	 X,Y	(R): RECTANGULAR COORDINATES IN KM
C	 ORGLAT	(R): ORIGIN PARALLEL +90 TO -90 DEG WITH NORTH POSITIVE
C	 ORGLON	(R): CENTRAL MERIDIAN (LONGITUDE) 0 TO +360 DEG
C			OR -180 TO +180 WITH EAST MORE POSITIVE
C	 IOPT	(I): EARTH RADIUS OPTION
C
C	OUTPUTS:
C	 LAT	(R): LATITUDE +90 TO -90 DEG WITH NORTH POSITIVE
C	 LON	(R): LONGITUDE 0 TO +360 DEG WITH EAST POSITIVE
C
	IMPLICIT NONE
	REAL LAT, LON, X, Y, ORGLAT, ORGLON
	INTEGER  IOPT
	REAL ORGLON1, RHO, C, T1, T2, X1, Y1
	REAL F, ERADEARTH, RADEARTH, ERA
	REAL DTR,RTD
	DATA DTR/0.01745329252D0/,RTD/57.2957795/
C
	DATA RADEARTH/6378.135/, F/298.26/	! EQUITORIAL EARTH RADIUS, 1/F
C						! WGS 72 MODEL VALUES
	ORGLON1=MOD(ORGLON+720.0,360.0)
C
C	COMPUTE LOCAL RADIUS OF THE EARTH AT CENTER OF IMAGE
C
	ERADEARTH=6378.0			! USE FIXED NOMINAL VALUE
	IF (IOPT.EQ.2) THEN			! LOCAL RADIUS
	  ERA=(1.-1./F)
	  ERADEARTH=RADEARTH*ERA/SQRT(ERA*ERA*COS(DTR*ORGLAT)**2+
     *         SIN(DTR*ORGLAT)**2)
	ENDIF
C
	X1=X/ERADEARTH
	Y1=Y/ERADEARTH
	RHO=X1*X1+Y1*Y1
	IF (RHO.GT.0.0) THEN
		RHO=SQRT(RHO)
		C=2*ASIN(RHO*0.5)
		LAT=RTD*ASIN(COS(C)*SIN(DTR*ORGLAT)+Y1*SIN(C)*COS(DTR*ORGLAT)/RHO)
	ELSE
		LAT=ORGLAT
	ENDIF
	IF (ABS(ORGLAT).NE.90.0) THEN
		IF (RHO.EQ.0.0) THEN
			LON=ORGLON1
		ELSE
			T1=X1*SIN(C)
			T2=RHO*COS(DTR*ORGLAT)*COS(C)-Y1*SIN(DTR*ORGLAT)*SIN(C)
			LON=ORGLON1+ATAN2(T1,T2)*RTD
		ENDIF
	ELSE IF (ORGLAT.EQ.90.0) THEN
		LON=ORGLON1+ATAN2(X1,-Y1)*RTD
	ELSE
		LON=ORGLON1+ATAN2(X1,Y1)*RTD
	ENDIF
	LON=MOD(LON+720.0,360.0)
	IF (LON.GT.180.0) LON=LON-360.0
	RETURN
	END