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 = BASE_PIXEL_DIMEN C ASCALE = RADIUS_EARTH/BASE_PIXEL_DIMENSION C C CODE WRITTEN BY: DGL 4 MAR 1995 C MODIFIED BY: DGL 23 JUL 2005 C REVISED BY: DGL NOV 2005 C + eliminated COSD and SIND calls with DTR C IMPLICIT NONE INTEGER IOPT REAL ALON, ALAT, THELON, THELAT, ASCALE REAL ARCTAND ! EXTERNAL FUNCTION REAL PI2, X1, Y1,TEMP REAL DTR,RTD DATA DTR/0.01745329252D0/,RTD/57.2957795/ 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(SIN(DTR*ALON)).GT.ABS(COS(DTR*ALON))) THEN TEMP=(X1/SIN(DTR*ALON))/ASCALE ELSE TEMP=(-Y1/COS(DTR*ALON))/ASCALE ENDIF IF (ABS(TEMP).LE.1.0) THEN ALAT=90.0-2.*ASIN(TEMP)*RTD ELSE ALAT=-90.0+2.*ASIN(1./TEMP)*RTD 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(COS(DTR*ALON)).GT.ABS(SIN(DTR*ALON))) THEN TEMP=(Y1/COS(DTR*ALON))/ASCALE ELSE TEMP=(X1/SIN(DTR*ALON))/ASCALE ENDIF IF (ABS(TEMP).LE.1.0) THEN ALAT=90.0-2.*ACOS(TEMP)*RTD ELSE WRITE (*,*) '*** ERROR IN EASE INVERSE COSINE ***',TEMP ALAT=0.0 ENDIF ELSE IF (IOPT.EQ.13) THEN ! EASE CYLINDRICAL ALON=((X1/ASCALE)/COS(DTR*30.0))*90.0/PI2 TEMP=(Y1*COS(DTR*30.0))/ASCALE IF (ABS(TEMP).LE.1.0) THEN ALAT=ASIN(TEMP)*RTD ELSE WRITE (*,*) '*** ERROR IN EASE INVERSE SINE ***',TEMP ALAT=SIGN(90.0,TEMP) ENDIF ENDIF RETURN END