C
C
	SUBROUTINE IPOLSTER(ALON,ALAT,X,Y,XLAM,SLAT)
C
C	COMPUTES THE INVERSE POLAR STEROGRAPHIC TRANSFORMATION FOR (X,Y)
C	GIVEN IN KM WITH REFERENCES LON,LAT=(XLAM,SLAT).
C	OUTPUT LON,LAT=ALON,ALAT
C
C	ALGORITHM IS THE SAME AS USED FOR PROCESSING ERS-1 SAR IMAGES
C	AS RECEIVED FROM M. DRINKWATER (1994).  UPDATED BY D. LONG TO
C	IMPROVE ACCURACY USING ITERATION WITH FORWARD TRANSFORM.
C
C	WRITTEN BY: DGL OCT 1994
C	MODIFIED BY: DGL MAR 1995
C	+ INCORPORATED ITERATIVE UPDATE TO IMPROVE ACCURACY
C
	IMPLICIT NONE
	REAL ALON, ALAT, X, Y, XLAM, SLAT
C
	real sind,cosd,tand,atan2d  ! for some stupid compilers
C
	INTEGER ICNT
	REAL XX,YY,R,RR,RERR,A,AA,AERR,SN1
	DOUBLE PRECISION E2,RE,RHO,E,E22,E23
        DOUBLE PRECISION CM,CHI,T,PI2,X1,Y1,SN,SLAT1
C
        REAL ARCTAND
        EXTERNAL ARCTAND
C
	DATA E2/0.006693883D0/
	DATA RE/6378.273D0/
	DATA PI2/1.570796327D0/
C
C	FIRST USE APPROXIMATE INVERSE CALCULATION
C
	E=SQRT(E2)
	E22=E2*E2
	E23=E2*E2*E2
	X1=X
	Y1=Y
	RHO=X1*X1+Y1*Y1
	IF (RHO.GT.0.0) RHO=SQRT(RHO)
	IF (RHO.LT.0.05) THEN
		ALON=XLAM
		ALAT=SIGN(90.0,SLAT)
	ELSE
		SN=1.0D0
		SLAT1=SLAT
		IF (SLAT.LT.0.0) THEN
			SN=-1.0D0
			SLAT1=-SLAT
		ENDIF
		CM=COSD(SLAT1)/SQRT(1.0D0-E2*SIND(SLAT1)**2)
		T=TAND(45.0-0.5*SLAT1)/
     *		  ((1.D0-E*SIND(SLAT1))/(1.D0+E*SIND(SLAT1)))**(E*0.5D0)
		T=RHO*T/(RE*CM)
		CHI=PI2-2.0D0*ATAN(T)
		T=CHI+
     *		  (0.5D0*E2+5.0D0*E22/24.0D0+E23/12.0D0)*SIN(2.0D0*CHI)+
     *		  (7.0D0*E22/48.0D0+29.0D0*E23/240.0D0)*SIN(4.0D0*CHI)+
     *		  (7.0D0*E23/120.0D0)*SIN(6.0D0*CHI)
		ALAT=SN*(T*90.0D0/PI2)
		ALON=SN*ATAN2D(SN*X1,-SN*Y1)+XLAM
		IF (ALON.LT.-180.0) ALON=ALON+360.0
		IF (ALON.GT.180.0) ALON=ALON-360.0
	ENDIF
C
C	USING THE APPROXIMATE RESULT AS A STARTING POINT, ITERATE TO IMPROVE
C	THE ACCURACY OF THE INVERSE SOLUTION
C
	SN1=1.0
	IF (SLAT.LT.0.0) SN1=-1.0
	A=ARCTAND(Y,X)
	R=SQRT(X*X+Y*Y)
	ICNT=0
10	CONTINUE
	CALL POLSTER(ALON,ALAT,XX,YY,XLAM,SLAT)
	RR=SQRT(XX*XX+YY*YY)
	RERR=SN1*(RR-R)/180.0
	AA=ARCTAND(YY,XX)
	AERR=AA-A
	IF (ABS(AERR).GT.180.0) AERR=360.0-AERR
C
C	CHECK FOR CONVERGENCE
C
	IF ((ABS(RERR).LT.0.001.AND.ABS(AERR).LT.0.001)
     *       .OR.ICNT.GT.9) GOTO 40
C
C	CONSTRAIN UPDATES
C
	ALON=ALON+AERR
	IF (ABS(ALON).GT.360.0) ALON=MOD(ALON,360.0)
	IF (ALAT*SLAT.LT.0.0) THEN
		RERR=RERR*(1.0-SIND(ABS(ALAT)))
		IF (ABS(RERR).GT.2.0) RERR=SIGN(2.0,RERR)/FLOAT(ICNT)
	ENDIF
	ALAT=ALAT+RERR
	IF (ABS(ALAT).GT.90.0) ALAT=SIGN(90.0,ALAT)
	ICNT=ICNT+1
	GOTO 10
40	CONTINUE
	IF (ABS(ALON).GT.360.0) ALON=MOD(ALON,360.0)
	RETURN	
	END