	SUBROUTINE IEASE2GRID(IOPT,ALON,ALAT,THELON,THELAT,ASCALE,BSCALE)
C
C	COMPUTE THE INVERSE "EASE2" GRID TRANSFORM
C
C	GIVEN THE IMAGE TRANSFORMATION COORDINATES (THELON,THELAT) AND
C	THE CORRESPONDING LON,LAT (ALON,ALAT) IS COMPUTED
C	USING THE "EASE GRID" (VERSION 2.0) TRANSFORMATION GIVEN IN IDL
C	SOURCE CODE SUPPLIED BY MJ BRODZIK
C
C	inputs:
C	  iopt: projection type 8=EASE2 N, 9-EASE2 S, 10=EASE2 T/M
C	  thelon: X coordinate in pixels (can be outside of image)
C	  thelat: Y coordinate in pixels (can be outside of image)
C         ascale and bscale should be integer valued)
C	  ascale: grid scale factor (0..5)  pixel size is (bscale/2^ascale)
C	  bscale: base grid scale index (ind=int(bscale))
C
C       for definitions of ind and isc, see EASE2_MAP_INFO
C
C	outputs:
C	  alon, alat: lon, lat location in deg  (can be outside of image)
C
C	WRITTEN BY: DGL 12 MAR 2014
C
	integer iopt
	real alat,alon,thelon,thelat,ascale,bscale

	integer isc,ind
	double precision map_equatorial_radius_m, map_eccentricity
	double precision e2, map_reference_latitude
	double precision map_reference_longitude
	double precision map_second_reference_latitude
	double precision sin_phi1, cos_phi1, kz
	double precision map_scale, r0, s0, epsilon 
	integer bcols, brows

	double precision easeconv_normalize_degrees

	double precision lam, arg, phi, beta, qp, rho2, x, y, e4, e6

	double precision RTD
	parameter (RTD=57.29577951308232)


	ind=int(bscale)
	isc=int(ascale)

	call EASE2_MAP_INFO(iopt,isc,ind,
     1    map_equatorial_radius_m, map_eccentricity,
     2     e2, map_reference_latitude, map_reference_longitude,
     3     map_second_reference_latitude, sin_phi1, cos_phi1, kz, 
     4     map_scale, bcols, brows, r0, s0, epsilon)

    
	e4 = e2 * e2
	e6 = e4 * e2

C       qp is the function q evaluated at phi = 90.0 deg
	qp = (1.0d0-e2) * ((1.0d0 /(1.0d0 - e2)) 
     1			 - (1.0d0 /(2.0d0*map_eccentricity)) 
     2			 * log((1.0d0-map_eccentricity) 
     3                         /(1.0d0+map_eccentricity)))

	x = (thelon - r0 - 0.5d0) * map_scale
	y = (thelat - 0.5d0 - s0) * map_scale ! 0 at bottom

	if (iopt==8) then          ! EASE2 grid north
	   rho2 = ( x * x ) + ( y * y )
	   arg=1.0d0- (rho2/(map_equatorial_radius_m * 
     1                       map_equatorial_radius_m * qp))
	   if (arg .gt.  1.0d0) arg=1.0      
	   if (arg .lt. -1.0d0) arg=-1.0
	   beta = asin( arg )
	   lam = atan2( x, -y )

	else if (iopt==9) then     ! EASE2 grid south
	   rho2 = ( x * x ) + ( y * y )
	   arg = 1.0d0-(rho2/(map_equatorial_radius_m * 
     1                        map_equatorial_radius_m * qp))
	   if (arg .gt.  1.0d0) arg=1.0d0      
	   if (arg .lt. -1.0d0) arg=-1.0d0
	   beta = -asin( arg )
	   lam = atan2( x, y )

	else if (iopt==10) then    ! EASE2 cylindrical
	   arg = 2.0d0*y*kz /(map_equatorial_radius_m * qp)
	   if (arg .gt.  1.0d0) arg=1.0d0    
	   if (arg .lt. -1.0d0) arg=-1.0d0
	   beta = asin( arg )
	   lam = x / (map_equatorial_radius_m * kz)
	else
	   beta = 0.0
	   write (*,*) '*** invalid EASE2 projection in iease2grid',iopt
	endif

	phi = beta 
     1     + ((( e2 / 3.0d0 ) + (( 31.0d0 / 180.0d0 ) * e4 ) 
     2	   + ((  517.0d0 / 5040.0) * e6 )) * sin( 2.0d0 * beta )) 
     3     + (((( 23.0d0 / 360.0d0) * e4) 
     4	   + ((  251.0d0 / 3780.0d0) * e6)) * sin( 4.0d0 * beta )) 
     5     + ((( 761.0d0 / 45360.0d0) * e6) * sin( 6.0d0 * beta ))
   
	alat = RTD * phi
	lam = map_reference_longitude + lam*RTD
	lam = easeconv_normalize_degrees( lam)
	alon = lam

	RETURN
	END
