C C SUBROUTINE IEASEGRID(IOPT,ALON,ALAT,THELON,THELAT,ASCALE) C C COMPUTES THE INVERSE "EASE" GRID TRANSFORM C C GIVEN THE IMAGE TRANSFORMATION COORDINATES (THELON,THELAT) AND C THE SCALE (ASCALE) THE CORRESPONDING LON,LAT (ALON,ALAT) IS COMPUTED C USING THE "EASE GRID" (VERSION 1.0) TRANSFORMATION GIVEN IN FORTRAN C SOURCE CODE SUPPLIED BY NSIDC C C IOPT IS EASE TYPE: IOPT=11=NORTH, IOPT=12=SOUTH, IOPT=13=CYLINDRICAL C C THE RADIUS OF THE EARTH USED IN THIS PROJECTION IS IMBEDDED INTO C ASCALE WHILE THE PIXEL DIMENSION IN KM IS IMBEDDED IN BSCALE C THE BASE VALUES ARE: RADIUS EARTH= 6371.228 KM C PIXEL DIMEN =25.067525 KM C THEN, BSCALE = 2*(BASE_PIXEL_DIMEN) C ASCALE = 4*(BASE_RADIUS_EARTH) C C THESE VALUES ARE SET IN THE GETSIR AND PUTSIRROUTINES C C THE INVERSE ALGORITHM IS SPECIFIED BY THE EASE GRID DEFINITION. C IT IS ACTUALLY ONLY AN APPROXIMATION TO THE MODIFIED POLAR C STEREOGRAPHIC PROJECTION USED IN THE EASE GRID DEFINITION. C C CODE WRITTEN BY: DGL 4 MAR 1995 C IMPLICIT NONE INTEGER IOPT REAL ALON, ALAT, THELON, THELAT, ASCALE real sind,cosd,acosd,asind ! for some stupid compilers REAL ARCTAND ! EXTERNAL FUNCTION REAL PI2, X1, Y1,TEMP DATA PI2/1.57079633/ ! Pi/2 C X1=THELON Y1=THELAT IF (IOPT.EQ.11) THEN ! EASE GRID NORTH ALON=ARCTAND(X1,-Y1) IF (ABS(SIND(ALON)).GT.ABS(COSD(ALON))) THEN TEMP=(X1/SIND(ALON))/ASCALE ELSE TEMP=(-Y1/COSD(ALON))/ASCALE ENDIF IF (ABS(TEMP).LE.1.0) THEN ALAT=90.0-2.*ASIND(TEMP) ELSE ALAT=-90.0+2.*ASIND(1./TEMP) WRITE (*,*) '*** ERROR IN EASE INVERSE SINE ***',TEMP C ALAT=SIGN(90.0,TEMP) ENDIF ELSE IF (IOPT.EQ.12) THEN ! EASE GRID SOUTH ALON=ARCTAND(X1,Y1) IF (ABS(COSD(ALON)).GT.ABS(SIND(ALON))) THEN TEMP=(X1/COSD(ALON))/ASCALE ELSE TEMP=(Y1/SIND(ALON))/ASCALE ENDIF IF (ABS(TEMP).LE.1.0) THEN ALAT=90.0-2.*ACOSD(TEMP) ELSE WRITE (*,*) '*** ERROR IN EASE INVERSE COSINE ***',TEMP ALAT=0.0 ENDIF ELSE IF (IOPT.EQ.13) THEN ! EASE CYLINDRICAL ALON=((X1/ASCALE)/COSD(30.0))*90.0/PI2 TEMP=(Y1/COSD(30.0))/ASCALE IF (ABS(TEMP).LE.1.0) THEN ALAT=ASIND(TEMP) ELSE WRITE (*,*) '*** ERROR IN EASE INVERSE SINE ***',TEMP ALAT=SIGN(90.0,TEMP) ENDIF ENDIF RETURN END