C
C **************************************************************************
C
C       conventional F77 code used in the F90 interface
C
C **************************************************************************
C
        REAL FUNCTION ARCTAND(Y,X)
C       REVISED BY: DGL NOV 2005
C        + changed atan2d to atan2 with constant
        IF (Y.EQ.0.0.AND.X.EQ.0.0) THEN
                ARCTAND=0.0
        ELSE
                ARCTAND=ATAN2(Y,X)*57.2957795
        ENDIF
        RETURN
        END

        SUBROUTINE EASE2GRID(IOPT,ALAT,ALON,THELON,THELAT,ASCALE,BSCALE)
C
C        COMPUTE THE FORWARD "EASE2" GRID TRANSFORMATION
C
C        GIVEN THE IMAGE TRANSFORMATION COORDINATES (THELON,THELAT) AND
C        THE CORRESPONDING LON,LAT (ALON,ALAT) IS COMPUTED
C        USING THE "EASE2 GRID" (VERSION 2.0) TRANSFORMATION GIVEN IN IDL
C        SOURCE CODE SUPPLIED BY MJ BRODZIK
C        RADIUS EARTH=6378.137 KM (WGS 84)
C        MAP ECCENTRICITY=0.081819190843 (WGS84)
C
C        inputs:
C          iopt: projection type 8=EASE2 N, 9-EASE2 S, 10=EASE2 T/M
C          alon, alat: lon, lat (deg) to convert (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          NSIDC .grd file for isc=0
C           project type    ind=0     ind=1         ind=3
C              N         EASE2_N25km EASE2_N30km EASE2_N36km  
C             S         EASE2_S25km EASE2_S30km EASE2_S36km 
C             T/M       EASE2_T25km EASE2_M25km EASE2_M36km 
C
C          cell size (m) for isc=0 (scale is reduced by 2^isc)
C           project type    ind=0     ind=1            ind=3
C              N          25000.0     30000.0         36000.0
C             S          25000.0     30000.0         36000.0
C             T/M       T25025.26   M25025.2600081  M36032.220840584
C              
C          for a given base cell size isc is related to NSIDC .grd file names
C             isc        N .grd name   S .grd name   T .grd name
C              0              EASE2_N25km     EASE2_S25km     EASE2_T25km  
C              1              EASE2_N12.5km   EASE2_S12.5km   EASE2_T12.5km  
C              2              EASE2_N6.25km   EASE2_S6.25km   EASE2_T6.25km  
C              3              EASE2_N3.125km  EASE2_S3.125km  EASE2_T3.125km  
C              4              EASE2_N1.5625km EASE2_S1.5625km EASE2_T1.5625km  
C
C        outputs:
C          thelon: X coordinate in pixels (can be outside of image)
C          thelat: Y coordinate in pixels (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 dlon,phi,lam,sin_phi,q,qp,rho,x,y

        double precision easeconv_normalize_degrees

        double precision DTR
        parameter (DTR=0.01745329241994)

        ind=int(bscale)
        isc=int(ascale)
        dlon=alon
        phi= DTR * alat

        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)

        dlon = dlon - map_reference_longitude   
        dlon = easeconv_normalize_degrees( dlon )
        lam = DTR * dlon
  
        sin_phi=sin(phi)
        q=(1.0d0-e2)*((sin_phi/(1.0d0-e2*sin_phi*sin_phi)) 
     1        - (1.0d0/(2.0d0*map_eccentricity)) 
     2        * log((1.0d0-map_eccentricity*sin_phi) 
     3        / ( 1.0d0+map_eccentricity*sin_phi)))

        if (iopt==8) then          ! EASE2 grid north
           qp = 1.0d0-((1.0d0-e2)/(2.0d0*map_eccentricity) 
     1                   * log((1.0d0-map_eccentricity) 
     2             / (1.0d0 + map_eccentricity)))
           if (abs(qp-q).lt.epsilon) then 
              rho = 0.0
           else
              rho = map_equatorial_radius_m * sqrt( qp - q )
           endif
           x =  rho * sin( lam )
           y = -rho * cos( lam )

        else if (iopt==9) then     ! EASE2 grid south
           qp = 1.0d0-((1.0d0-e2)/( 2.0d0*map_eccentricity) 
     1                   * log((1.0d0 - map_eccentricity) 
     2                          / (1.0d0 + map_eccentricity)))
           if ( abs( qp + q ) < epsilon ) then
             rho = 0.0
           else
             rho = map_equatorial_radius_m * sqrt( qp + q )
           endif
           x = rho * sin( lam )
           y = rho * cos( lam )

        else if (iopt==10) then    ! EASE2 cylindrical
           x =   map_equatorial_radius_m * kz * lam
           y = ( map_equatorial_radius_m * q ) / ( 2.0D0 * kz )
        else
           write (*,*) '*** invalid EASE2 projection in ease2grid',iopt
           x = 0.0
           y = 0.0
        endif

        thelon = r0 + ( x / map_scale ) + 0.5
        thelat = s0 + ( y / map_scale ) + 0.5 ! 0 at bottom

        RETURN
        END

        SUBROUTINE 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)
C
C define EASE2 grid information 
C
C  inputs (integer)
C    iopt: projection type 8=EASE2 N, 9-EASE2 S, 10=EASE2 T/M
C    isc:  scale factor 0..5 grid size is (basesize(ind))/2^isc
C    ind:  base grid size index   (map units per cell in m
C 
C          NSIDC .grd file for isc=0
C           project type    ind=0     ind=1         ind=3
C               N         EASE2_N25km EASE2_N30km EASE2_N36km  
C              S         EASE2_S25km EASE2_S30km EASE2_S36km 
C              T/M       EASE2_T25km EASE2_M25km EASE2_M36km 
C
C          cell size (m) for isc=0 (scale is reduced by 2**isc)
C           project type    ind=0     ind=1            ind=3
C               N          25000.0     30000.0         36000.0
C              S          25000.0     30000.0         36000.0
C              T/M       T25025.26   M25025.26000  M36032.220840584
C              
C          for a given base cell size isc is related to NSIDC .grd file names
C             isc        N .grd name   S .grd name   T .grd name
C              0              EASE2_N25km     EASE2_S25km     EASE2_T25km  
C              1              EASE2_N12.5km   EASE2_S12.5km   EASE2_T12.5km  
C              2              EASE2_N6.25km   EASE2_S6.25km   EASE2_T6.25km  
C              3              EASE2_N3.125km  EASE2_S3.125km  EASE2_T3.125km  
C              4              EASE2_N1.5625km EASE2_S1.5625km EASE2_T1.5625km  
C
C  outputs (double precision)
C    map_equatorial_radius_m  EASE2 Earth equitorial radius (km) [WGS84]
C    map_eccentricity         EASE2 Earth eccentricity [WGS84]
C    map_reference_latitude   Reference latitude (deg) 
C    map_reference_longitude  Reference longitude (deg)
C    map_second_reference_latitude Secondary reference longitude* (deg)
C    sin_phi1, cos_phi1 kz    EASE2 Cylin parameters*
C    map_scale                EASE2 map projection pixel size (km)
C    bcols, brows,            EASE2 grid size in pixels
C    r0, s0                   EASE2 base projection size in pixels
C    epsilon                  EASE2 near-polar test factor
C
C    *these parameters only assigned values if projection is T
C
C   DGL Aug 2, 2014 + modified EASE2T vertical dimension from 538 to 540
C   DGL Aug 12, 2014 + modified EASE2 global M25km base scale to 25025.26000
C
        integer iopt,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 DTR
        parameter (DTR=0.01745329241994)

        double precision base
        integer m, nx, ny

        m=2**isc   ! compute power-law scale factor

        map_equatorial_radius_m = 6378137.0d0  ! WGS84
        map_eccentricity = 0.081819190843d0    ! WGS84
        e2 = map_eccentricity *  map_eccentricity
        map_reference_longitude = 0.0

        if (iopt==8) then          ! EASE2 grid north
           map_reference_latitude = 90.0d0
           if (ind==1) then           ! EASE2_N30km.gpd
              base=30000.0d0
              nx=600
              ny=600        
           else if (ind==2) then      ! EASE2_N36km.gpd
              base=36000.0d0    
              nx=500
              ny=500        
           else                       ! EASE2_N25km.gpd
              base=25000.0d0
              nx=720
              ny=720        
           endif
        else if (iopt==9) then     ! EASE2 grid south
           map_reference_latitude = -90.0d0
           if (ind==1) then           ! EASE2_S30km.gpd
              base=30000.0d0
              nx=600
              ny=600        
           else if (ind==2) then      ! EASE2_S36km.gpd
              base=36000.0d0    
              nx=500
              ny=500        
           else                       ! EASE2_S25km.gpd
              base=25000.0d0
              nx=720
              ny=720        
           endif
        else if (iopt==10) then    ! EASE2 cylindrical
           map_reference_latitude = 0.0d0
           map_second_reference_latitude = 30.0d0
           sin_phi1 = sin( DTR * map_second_reference_latitude)
           cos_phi1 = cos( DTR * map_second_reference_latitude)
           kz = cos_phi1 / sqrt( 1.0d0 - e2*sin_phi1*sin_phi1)
            if (ind==1) then           ! EASE2_M25km.gpd
              base=25025.26000d0 ! had trailing  81
              nx=1388
              ny=584        
           else if (ind==2) then      ! EASE2_M36km.gpd
              base=36032.220840584d0
              nx=964
              ny=406        
           else                       ! EASE2_T25km.gpd
              base=25025.26000d0
              nx=1388
              ny=540  ! was 538
           endif
        else
           base=25000.0d0
           nx=10
           ny=10
        endif

C       grid info
        map_scale = base / m
        !bcols = ceil(nx * m)
        !brows = ceil(ny * m)
        bcols = nx * m
        brows = ny * m
        r0 = (nx * m - 1) / 2.0d0
        s0 = (ny * m - 1) / 2.0d0

        RETURN
        END

        SUBROUTINE EASE2SF(iopt,ascale,bscale,fmap_scale,bcols,brows,fr0,fs0)
C
C return selected EASE2 grid information (external interface)
C
C  inputs (integer)
C    iopt: projection type 8=EASE2 N, 9-EASE2 S, 10=EASE2 T/M
C    (note: ascale, bscale should be integer valued)
C    ascale:  scale factor 0..5 grid size is (basesize(ind))/2^isc
C    bscale:  base grid size index   (map units per cell in m
C 
C          NSIDC .grd file for isc=0
C           project type    ind=0     ind=1         ind=3
C               N         EASE2_N25km EASE2_N30km EASE2_N36km  
C              S         EASE2_S25km EASE2_S30km EASE2_S36km 
C              T/M       EASE2_T25km EASE2_M25km EASE2_M36km 
C
C          cell size (m) for isc=0 (scale is reduced by 2**isc)
C           project type    ind=0     ind=1            ind=3
C               N          25000.0     30000.0         36000.0
C              S          25000.0     30000.0         36000.0
C              T/M       T25025.26   M25025.26000  M36032.220840584
C              
C          for a given base cell size isc is related to NSIDC .grd file names
C             isc        N .grd name   S .grd name   T .grd name
C              0              EASE2_N25km     EASE2_S25km     EASE2_T25km  
C              1              EASE2_N12.5km   EASE2_S12.5km   EASE2_T12.5km  
C              2              EASE2_N6.25km   EASE2_S6.25km   EASE2_T6.25km  
C              3              EASE2_N3.125km  EASE2_S3.125km  EASE2_T3.125km  
C              4              EASE2_N1.5625km EASE2_S1.5625km EASE2_T1.5625km  
C
C  outputs (real)
C    fmap_scale               EASE2 map projection pixel size (km)
C    bcols, brows,            EASE2 grid size in pixels
C    fr0, fs0                 EASE2 base projection size in pixels
C
      integer iopt, bcols, brows
      real ascale, bscale, fmap_scale, fr0, fs0
c
      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 

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

      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)

      fmap_scale=map_scale
      fr0=r0
      fs0=s0
      return
      end

        double precision function easeconv_normalize_degrees(dlon)
C
C       Return -180 <= dlon <= 180
C
        double precision dlon, out
        out=dlon
 10        continue
        if (out .lt. -180.0) then
           out=out+360.0
           goto 10
        endif
 20        continue
        if (out .gt. 180.0) then
           out=out-360.0
           goto 20
        endif
        easeconv_normalize_degrees=out
        return
        end

        SUBROUTINE EASEGRID(IOPT,ALAT,ALON,THELON,THELAT,ASCALE)
C
C        COMPUTES THE FORWARD "EASE1" GRID TRANSFORM
C
C        GIVEN A LAT,LON (ALAT,ALON) AND THE SCALE (ASCALE) THE IMAGE
C        TRANSFORMATION COORDINATES (THELON,THELAT) ARE COMUTED
C        USING THE "EASE GRID" (VERSION 1.0) TRANSFORMATION GIVEN IN FORTRAN
C        SOURCE CODE SUPPLIED BY NSIDC.
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        IOPT IS EASE1 TYPE: IOPT=11=NORTH, IOPT=12=SOUTH, IOPT=13=CYLINDRICAL
C
C        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 ALAT, ALON, THELON, THELAT, ASCALE
        REAL PI2
        REAL DTR
        DATA DTR/0.01745329252/
        DATA PI2/1.57079633/                ! Pi/2
C
        IF (IOPT.EQ.11) THEN                ! EASE GRID NORTH
                THELON= ASCALE*SIN(DTR*ALON)*SIN(DTR*(45.0-0.5*ALAT))
                THELAT=-ASCALE*COS(DTR*ALON)*SIN(DTR*(45.0-0.5*ALAT))
        ELSE IF (IOPT.EQ.12) THEN        ! EASE GRID SOUTH
                THELON=ASCALE*SIN(DTR*ALON)*COS(DTR*(45.0-0.5*ALAT))
                THELAT=ASCALE*COS(DTR*ALON)*COS(DTR*(45.0-0.5*ALAT))
        ELSE IF (IOPT.EQ.13) THEN        ! EASE CYLINDRICAL
                THELON=ASCALE*PI2*ALON*COS(DTR*30.0)/90.0
                THELAT=ASCALE*SIN(DTR*ALAT)/COS(DTR*30.0)
        ENDIF
        RETURN
        END
C
        SUBROUTINE F2IPIX(X,Y,IX,IY,NSX,NSY)
C
C     QUANTIZES F FLOATING POINT PIXEL VALUE TO THE "ACTUAL" INTEGER
C     PIXEL VALUE  (SEE LATLON2PIX)
C
        IF (X+0.0002.GE.1.0.AND.X+0.0002.LT.NSX+1) THEN
           IX=IFIX(X+0.0002)  ! ALLOW A SMALL AMOUNT OF ROUNDING
        ELSE
           IX=0
        ENDIF
        IF (Y+0.0002.GE.1.0.AND.Y+0.0002.LT.NSY+1) THEN
           IY=IFIX(Y+0.0002)  ! ALLOW A SMALL AMOUNT OF ROUNDING
        ELSE
           IY=0
        ENDIF
        RETURN
        END

        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          NSIDC .grd file for isc=0
C          project type    ind=0     ind=1         ind=3
C              N         EASE2_N25km EASE2_N30km EASE2_N36km  
C             S         EASE2_S25km EASE2_S30km EASE2_S36km 
C             T/M       EASE2_T25km EASE2_M25km EASE2_M36km 
C
C         cell size (m) for isc=0 (scale is reduced by 2^isc)
C          project type    ind=0     ind=1            ind=3
C              N          25000.0     30000.0         36000.0
C             S          25000.0     30000.0         36000.0
C             T/M       T25025.26   M25025.2600081  M36032.220840584
C              
C          for a given base cell size isc is related to NSIDC .grd file names
C             isc        N .grd name   S .grd name   T .grd name
C              0              EASE2_N25km     EASE2_S25km     EASE2_T25km  
C              1              EASE2_N12.5km   EASE2_S12.5km   EASE2_T12.5km  
C              2              EASE2_N6.25km   EASE2_S6.25km   EASE2_T6.25km  
C              3              EASE2_N3.125km  EASE2_S3.125km  EASE2_T3.125km  
C              4              EASE2_N1.5625km EASE2_S1.5625km EASE2_T1.5625km  
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

        SUBROUTINE IEASEGRID(IOPT,ALON,ALAT,THELON,THELAT,ASCALE)
C
C        COMPUTES THE INVERSE "EASE1" 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 EASE1 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
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
	        C=0.0
                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
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       REVISED BY: DGL NOV 2005
C        + eliminated COSD and SIND calls with DTR and RTD
C
        IMPLICIT NONE
        REAL ALON, ALAT, X, Y, XLAM, SLAT
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
        DOUBLE PRECISION DTR,RTD
        DATA DTR/0.01745329252D0/,RTD/57.2957795/
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=COS(DTR*SLAT1)/SQRT(1.0D0-E2*SIN(DTR*SLAT1)**2)
                T=TAN(DTR*(45.0-0.5*SLAT1))/
     *                  ((1.D0-E*SIN(DTR*SLAT1))/(1.D0+E*SIN(DTR*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*ATAN2(SN*X1,-SN*Y1)*RTD+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-SIN(DTR*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
C
C *********************************************************************
C
C        PROJECTION TRANSFORMATION ROUTINES
C
C *********************************************************************
C
        SUBROUTINE LAMBERT1(LAT,LON,X,Y,ORGLAT,ORGLON,IOPT)
C
C        COMPUTES THE 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
C
C        INPUTS:
C         LAT        (R): LATITUDE +90 TO -90 DEG WITH NORTH POSITIVE
C         LON        (R): LONGITUDE 0 TO +360 DEG WITH EAST POSITIVE
C                        OR -180 TO +180 WITH EAST MORE POSITIVE
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         X,Y        (R): RECTANGULAR COORDINATES IN KM
C
        IMPLICIT NONE
        REAL LAT, LON, X, Y, ORGLAT, ORGLON
        INTEGER  IOPT
        REAL LON1, ORGLON1, DENOM, AK
        REAL F, ERADEARTH, RADEARTH, ERA
        REAL DTR
        DATA DTR/0.01745329252/
C
        DATA RADEARTH/6378.135/, F/298.26/        ! EQUITORIAL EARTH RADIUS, 1/F
C                                                ! WGS 72 MODEL VALUES
        LON1=MOD(LON+720.0,360.0)
        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
        DENOM=1.0+SIN(DTR*ORGLAT)*SIN(DTR*LAT)+
     1                COS(DTR*ORGLAT)*COS(DTR*LAT)*COS(DTR*(LON1-ORGLON1))
        IF (DENOM.GT.0.0) THEN
                AK=SQRT(2.0/DENOM)
        ELSE
                WRITE(*,*) '*** DIVISION ERROR IN LAMBERT1 ROUTINE ***'
                AK=1.0
        ENDIF
        X=AK*COS(DTR*LAT)*SIN(DTR*(LON1-ORGLON1))
        Y=AK*(COS(DTR*ORGLAT)*SIN(DTR*LAT)-
     1                SIN(DTR*ORGLAT)*COS(DTR*LAT)*COS(DTR*(LON1-ORGLON1)))
        X=X*ERADEARTH
        Y=Y*ERADEARTH
        RETURN
        END
C
C
        SUBROUTINE LATLON2PIX(ALON,ALAT,X,Y,THELON,THELAT,
     $     IOPT,XDEG,YDEG,ASCALE,BSCALE,A0,B0)
C
C        CONVERT A LAT,LON COORDINATE (ALON,ALAT) TO AN IMAGE PIXEL LOCATION
C        (X,Y) (IN FLOATING POINT) AND IMAGE X,Y COORDINATE (THELON,THELAT).
C        TO COMPUTE INTEGER PIXEL INDICES (IX,IY): CHECK TO INSURE
C        1 <= X < NSX+1 AND 1 <= X < NSX+1 THEN IX=IFIX(X) IY=IFIX(Y)
C
C        WRITTEN BY:        DGL MARCH 1997
C        REVISED BY:        DGL 12 MAR 2014 + ADDED EASE2
C
        REAL X,Y,THELON,THELAT,ALON,ALAT,XDEG,YDEG,ASCALE,BSCALE,A0,B0
        INTEGER IOPT
C
        IF (IOPT.EQ.-1) THEN                ! IMAGE ONLY (CAN'T TRANSFORM!)
                X=ASCALE*(ALON-A0)+1.0
                Y=BSCALE*(ALAT-B0)+1.0
                THELON=ALON
                THELAT=ALAT
        ELSE IF (IOPT.EQ.0) THEN        ! LAT/LONG (ASCALE AND BSCALE > 0)
                THELON=ALON
                THELAT=ALAT
                X=ASCALE*(THELON-A0)+1.0
                Y=BSCALE*(THELAT-B0)+1.0
        ELSE IF (IOPT.EQ.1.OR.IOPT.EQ.2) THEN        ! LAMBERT (ASCALE AND BSCALE < 0)
                CALL LAMBERT1(ALAT,ALON,THELON,THELAT,YDEG,XDEG,IOPT)
                X=(THELON-A0)*ASCALE+1.0
                Y=(THELAT-B0)*BSCALE+1.0
        ELSE IF (IOPT.EQ.5) THEN        ! POLAR STEREOGRAPHIC
                CALL POLSTER(ALON,ALAT,THELON,THELAT,XDEG,YDEG)
                X=(THELON-A0)/ASCALE+1.0
                Y=(THELAT-B0)/BSCALE+1.0
        ELSE IF (IOPT.EQ.8.OR.IOPT.EQ.9.OR.IOPT.EQ.10) THEN        ! EASE2 GRID
                CALL EASE2GRID(IOPT,ALAT,ALON,THELON,THELAT,ASCALE,BSCALE)
                X=THELON+1.0-A0
                Y=THELAT+1.0-B0
        ELSE IF (IOPT.EQ.11.OR.IOPT.EQ.12.OR.IOPT.EQ.13) THEN        ! EASE1 GRID
                CALL EASEGRID(IOPT,ALAT,ALON,THELON,THELAT,ASCALE)
                THELON=THELON+XDEG
                THELAT=THELAT+YDEG
                X=THELON+1.0-(XDEG+A0)
                Y=THELAT+1.0-(YDEG+B0)
        ELSE                                ! UNKNOWN TRANSFORMATION TYPE
        ENDIF
C
        RETURN
        END
C
C **************************************************************************
C
        INTEGER FUNCTION LENGTH1(LINE)
C
C        COMPUTES THE LENGTH (NON-BLANK) OF A STRING
C
        CHARACTER*(*) LINE
        N=LEN(LINE)
        DO I=N,1,-1
                IF (LINE(I:I).NE.' '.AND.LINE(I:I).NE.CHAR(0)) THEN
                        LENGTH1=I
                        GOTO 10
                ENDIF
        END DO
        LENGTH1=1
10        RETURN
        END
C
C **************************************************************************
C
        INTEGER FUNCTION LENGTH(LINE)
C
C        COMPUTES THE LENGTH (NON-BLANK) OF A STRING
C
        CHARACTER*(*) LINE
        N=LEN(LINE)
        DO I=N,1,-1
                IF (LINE(I:I).NE.' '.AND.LINE(I:I).NE.CHAR(0)) THEN
                        LENGTH=I
                        GOTO 10
                ENDIF
        END DO
        LENGTH=1
10        RETURN
        END
C
C ***************************************************************************
C
C        IMAGE PIXEL TO LAT,LONG AND LAT,LONG TO IMAGE PIXEL
C       COORDINATE TRANSFORMATION ROUTINES.
C
C ***************************************************************************
C
        SUBROUTINE PIX2LATLON(X,Y,THELON,THELAT,ALON,ALAT,
     $     IOPT,XDEG,YDEG,ASCALE,BSCALE,A0,B0)
C
C        GIVEN AN IMAGE PIXEL LOCATION (X,Y) (AS FLOATING POINT NUMBERS)
C        COMPUTES THE X,Y TRANSFORMATION COORDINATES (THELON,THELAT) AND
C        THE LAT,LON COORDINATES (ALON,ALAT).   THE LAT,LON RETURNED 
C        CORRESPONDS TO THE LOWER-LEFT CORNER OF THE PIXEL.  IF LAT,LON
C        OF PIXEL CENTER IS DESIRED USE (X+0.5,Y+0.5) WHERE X,Y ARE INTEGER
C        VALUED.
C
C        NOTE: GIVEN INTEGER PIXEL INDICES, SET X=FLOAT(IX) Y=FLOAT(IY).
C        -- SEE LATLON2PIX.  WHILE ROUTINE WILL ATTEMPT TO CONVERT ANY (X,Y)
C        VALUES, ONLY (X,Y) VALUES WITH 1 <= X <= NSX+1 AND 1 <= Y <= NSY+1
C        ARE CONTAINED WITHIN IMAGE PIXELS.
C
C        WRITTEN BY: DGL MARCH 1997
C       REVISED BY: DGL SEPT 2005 + corrected EASE grid
C        REVISED BY: DGL 12 MAR 2014 + ADDED EASE2
C
        REAL X,Y,THELON,THELAT,ALON,ALAT,XDEG,YDEG,ASCALE,BSCALE,A0,B0
        INTEGER IOPT
C
        IF (IOPT.EQ.-1) THEN                ! IMAGE ONLY (CAN'T TRANSFORM!)
                THELON=(X-1.0)/ASCALE+A0
                THELAT=(Y-1.0)/BSCALE+B0
                ALON=THELON
                ALAT=THELAT
        ELSE IF (IOPT.EQ.0) THEN        ! LAT/LONG (ASCALE AND BSCALE > 0)
                THELON=(X-1.0)/ASCALE+A0
                THELAT=(Y-1.0)/BSCALE+B0
                ALON=THELON
                ALAT=THELAT
        ELSE IF (IOPT.EQ.1.OR.IOPT.EQ.2) THEN        ! LAMBERT (ASCALE AND BSCALE < 0)
                THELON=(X-1.0)/ASCALE+A0
                THELAT=(Y-1.0)/BSCALE+B0
                CALL ILAMBERT1(ALAT,ALON,THELON,THELAT,YDEG,XDEG,IOPT)
        ELSE IF (IOPT.EQ.5) THEN        ! POLAR STEREOGRAPHIC
                THELON=(X-1.0)*ASCALE+A0
                THELAT=(Y-1.0)*BSCALE+B0
                CALL IPOLSTER(ALON,ALAT,THELON,THELAT,XDEG,YDEG)
        ELSE IF (IOPT.EQ.8.OR.IOPT.EQ.9.OR.IOPT.EQ.10) THEN        ! EASE2 GRID
                THELON=X-1.0+A0
                THELAT=Y-1.0+B0
                CALL IEASE2GRID(IOPT,ALON,ALAT,THELON,THELAT,ASCALE,BSCALE)
        ELSE IF (IOPT.EQ.11.OR.IOPT.EQ.12.OR.IOPT.EQ.13) THEN        ! EASE1 GRID
                THELON=X-XDEG-1.0+(XDEG+A0)
                THELAT=Y-YDEG-1.0+(YDEG+B0)
                CALL IEASEGRID(IOPT,ALON,ALAT,THELON,THELAT,ASCALE)
        ELSE                                ! UNKNOWN TRANSFORM
        ENDIF
        RETURN
        END
C
C
        SUBROUTINE POLSTER(ALON,ALAT,X1,Y1,XLAM,SLAT)
C
C        COMPUTES THE POLAR STEROGRAPHIC TRANSFORMATION FOR A LON,LAT
C        INPUT OF (ALON,ALAT) WITH REFERENCE ORIGIN  LON,LAT=(XLAM,SLAT).
C        OUTPUT IS (X1,Y1) IN KM
C
C        ALGORITHM IS THE SAME AS USED FOR PROCESSING ERS-1 SAR IMAGES
C        AS RECEIVED FROM M. DRINKWATER (1994) 
C
C        WRITTEN BY: DGL OCT 1994
C       REVISED BY: DGL JUN 1998
C        + fixed sign error in computing CM
C       REVISED BY: DGL NOV 2005
C        + eliminated COSD, SIND, and TAND calls with DTR
C
        IMPLICIT NONE
        REAL ALON, ALAT, X1, Y1, XLAM, SLAT
C
        DOUBLE PRECISION E2,RE,E,T,TX,TY,CM,RHO,RLAT,SN
        DATA E2/0.006693883D0/
        DATA RE/6378.273D0/
        DOUBLE PRECISION DTR
        DATA DTR/0.01745329252D0/
C
        E=SQRT(E2)
        IF (SLAT.LT.0.0) THEN
                SN=-1.0
                RLAT=-ALAT
        ELSE
                SN=1.0
                RLAT=ALAT
        ENDIF
        T=((1.D0-E*SIN(DTR*RLAT))/(1.D0+E*SIN(DTR*RLAT)))**(E*0.5D0)
        TY=TAN(DTR*(45.0D0-0.5D0*RLAT))/T
        IF (SLAT.LT.0.0) THEN
                RLAT=-SLAT
        ELSE
                RLAT=SLAT
        ENDIF
        T=((1.D0-E*SIN(DTR*RLAT))/(1.D0+E*SIN(DTR*RLAT)))**(E*0.5D0)
        TX=TAN(DTR*(45.0D0-0.5D0*RLAT))/T
        CM=COS(DTR*RLAT)/SQRT(1.D0-E2*SIN(DTR*RLAT)**2)
        RHO=RE*CM*TY/TX
        X1= (SN*SIN(DTR*(SN*ALON-XLAM)))*RHO
        Y1=-(SN*COS(DTR*(SN*ALON-XLAM)))*RHO
        RETURN
        END
C
C **************************************************************************
C
        SUBROUTINE PRINTHEAD3(nhead,ndes,nhtype,idatatype,
     *     nsx,nsy,xdeg,ydeg,ascale,bscale,a0,b0,iopt,ioff,iscale,
     *     ixdeg_off,iydeg_off,ideg_sc,iscale_sc,ia0_off,ib0_off,i0_sc,
     *     iyear,isday,ismin,ieday,iemin,iregion,itype,
     *     ipol,ifreqhm,ispare1,anodata,vmin,vmax,
     *     SENSOR,TITLE,TYPE,TAG,CRPROC,CRTIME,
     *     DESCRIP,LDES,IAOPT,NIA)
C
C        PRINTS OUT SIR FILE HEADER (version 3)
C
C       Revised 12 Mar 2014 + added EASE2
C
        CHARACTER*(*) SENSOR,TITLE,TYPE,TAG,CRPROC,CRTIME,DESCRIP
        INTEGER IAOPT(1)
        INTEGER bcols, brows
        REAL ascale, bscale, fmap_scale, fr0, fs0
C
        IF (NHTYPE.LT.16) WRITE (*,*) 'Old style header',nhtype
C
         WRITE (*,*) 'Title:  "',TITLE(1:LENGTH1(TITLE)),'"'
        WRITE (*,*) 'Sensor: "',SENSOR(1:LENGTH1(SENSOR)),'"'
        IF (NHTYPE.GT.16) THEN
           WRITE (*,*) 'Type:   "',TYPE(1:LENGTH1(TYPE)),'"'
           WRITE (*,*) 'Tag:    "',TAG(1:LENGTH1(TAG)),'"'
           WRITE (*,*) 'Creator:"',CRPROC(1:LENGTH1(CRPROC)),'"'
           WRITE (*,*) 'Created:"',CRTIME(1:LENGTH1(CRTIME)),'"'
        ENDIF
        IF (IOPT.EQ.-1) THEN                ! IMAGE ONLY
         WRITE (*,*) 'Rectangular image-only form: ',ioff,iscale,iopt
         WRITE (*,*) ' x,y scale   ',ascale,bscale
         WRITE (*,*) ' x,y offsets:',a0,b0
         WRITE (*,*) ' x,y span:   ',xdeg,ydeg
        ELSE IF (IOPT.EQ.0) THEN        ! LAT/LONG
         WRITE (*,*) 'Lat/Lon Rectangular form: ',ioff,iscale,iopt
         WRITE (*,*) ' X and Y scales (in pixels/deg)',ascale,bscale
         WRITE (*,*) ' Offset values:',a0,b0
         WRITE (*,*) ' Degrees:      ',xdeg,ydeg
        ELSE IF (IOPT.EQ.1.OR.IOPT.EQ.2) THEN        ! LAMBERT
         WRITE (*,*) 'Lambert form: ',ioff,iscale,iopt
         WRITE (*,*) ' X and Y scales (in km/pixel)',1./ascale,1./bscale
         WRITE (*,*) ' Lower Left Corner:',a0,b0
         WRITE (*,*) ' Center Point:     ',xdeg,ydeg
         IF (IOPT.EQ.2) WRITE (*,*) ' Local radius used'
        ELSE IF (IOPT.EQ.5) THEN        ! POLAR STEREOGRAPHIC
         WRITE (*,*) 'Polar Sterographic form:     ',ioff,iscale,iopt
         WRITE (*,*) ' X and Y scales (in km/pixel)',ascale,bscale
         WRITE (*,*) ' Lower Left Corner: (km)     ',a0,b0
         WRITE (*,*) ' Center Lon, Lat:            ',xdeg,ydeg
        ELSE IF (IOPT.EQ.8.OR.IOPT.EQ.9) THEN ! EASE2 NORTH OR SOUTH 
         WRITE (*,*) 'EASE2 polar azimuthal form: ',ioff,iscale,iopt
         CALL EASE2SF(iopt,ascale,bscale,fmap_scale,bcols,brows,fr0,fs0)
         WRITE (*,*) ' Nscale, Etype:',ascale,bscale,'  ',fmap_scale*0.001,' km/pix'
         WRITE (*,*) ' Map origin (col/row):',a0,b0
         WRITE (*,*) ' Map center (col/row):',xdeg,ydeg
        ELSE IF (IOPT.EQ.10) THEN        ! EASE2 CYLINDRICAL
         WRITE (*,*) 'EASE2 cylindrical form: ',ioff,iscale,iopt
         CALL EASE2SF(iopt,ascale,bscale,fmap_scale,bcols,brows,fr0,fs0)
         WRITE (*,*) ' Nscale, Etype:',ascale,bscale,'  ',fmap_scale*0.001,' km/pix'
         WRITE (*,*) ' Map origin (col/row):',a0,b0
         WRITE (*,*) ' Map center (col/row):',xdeg,ydeg
        ELSE IF (IOPT.EQ.11.OR.IOPT.EQ.12) THEN        ! EASE1 NORTH OR SOUTH
         WRITE (*,*) 'EASE1 polar azimuthal form: ',ioff,iscale,iopt
         WRITE (*,*) ' A, B scales:',ascale,bscale
         WRITE (*,*) ' Map origin (col/row):',a0,b0
         WRITE (*,*) ' Map center (col/row):',xdeg,ydeg
        ELSE IF (IOPT.EQ.13) THEN        ! EASE1 CYLINDRICAL
         WRITE (*,*) 'EASE1 cylindrical form: ',ioff,iscale,iopt
         WRITE (*,*) ' A, B scales:',ascale,bscale
         WRITE (*,*) ' Map origin (col/row):',a0,b0
         WRITE (*,*) ' Map center (col/row):',xdeg,ydeg
        ELSE                                ! UNKNOWN TRANSFORMATION
         WRITE (*,*) '*** Unrecognized SIR file form ***',ioff,iscale,iopt
         WRITE (*,*) ' a,b scale   ',ascale,bscale
         WRITE (*,*) ' a,b offsets:',a0,b0
         WRITE (*,*) ' x,y span:   ',xdeg,ydeg
        ENDIF
        WRITE (*,*) ' Pixels:       ',nsx,' x ',nsy
        WRITE (*,*) ' Year:  ',iyear,'  Region: ',iregion
        WRITE (*,*) ' Start: ',isday,ismin,' End:   ',ieday,iemin
        WRITE (*,*) ' Type:  ',itype,'  Transform option: ',iopt
        IF (NHTYPE.GT.16) THEN
           IF (IPOL.EQ.1) THEN
              WRITE (*,*) ' Frequency:  ',ifreqhm*0.1,' GHz  H Pol '
           ELSE IF (IPOL.EQ.2) THEN
              WRITE (*,*) ' Frequency:  ',ifreqhm*0.1,' GHz  V Pol '
           ELSE
              WRITE (*,*) ' Frequency:  ',ifreqhm*0.1,' GHz  Pol: ',ipol
           ENDIF
           WRITE (*,*) ' Data:  ',idatatype,'  Header: ',nhead,nhtype,ndes
           WRITE (*,*) ' Offset:',ioff,'  Scale:',iscale
           WRITE (*,*) ' Head scale factors:',ixdeg_off,iydeg_off,ideg_sc,
     $                iscale_sc,ia0_off,ib0_off,i0_sc
           WRITE (*,*) ' Nodata:',anodata,'  Vn: ',vmin,'  Vx:',vmax
           IF (LDES.GT.0) WRITE (*,*) 'Description: "',DESCRIP(1:LDES),'"'
           IF (NIA.GT.0) THEN
              DO I=1,NIA
                 WRITE (*,*) ' Aopt: ',i,' ',iaopt(i)
              END DO
           ENDIF
        ENDIF
        RETURN
        END
C
C **************************************************************************
C
        SUBROUTINE PRINTHEAD(nhead,ndes,nhtype,idatatype,
     *     nsx,nsy,xdeg,ydeg,ascale,bscale,a0,b0,iopt,ioff,iscale,
     *     iyear,isday,ismin,ieday,iemin,iregion,itype,
     *     ipol,ifreqhm,ispare1,anodata,vmin,vmax,
     *     SENSOR,TITLE,TYPE,TAG,CRPROC,CRTIME,
     *     DESCRIP,LDES,IAOPT,NIA)
C
C        PRINTS OUT SIR FILE HEADER (version 3)
C
        CHARACTER*(*) SENSOR,TITLE,TYPE,TAG,CRPROC,CRTIME,DESCRIP
        INTEGER IAOPT(1)
c
        call PRINTHEAD3(nhead,ndes,nhtype,idatatype,
     *     nsx,nsy,xdeg,ydeg,ascale,bscale,a0,b0,iopt,ioff,iscale,
     *     ixdeg_off,iydeg_off,ideg_sc,iscale_sc,ia0_off,ib0_off,i0_sc,
     *     iyear,isday,ismin,ieday,iemin,iregion,itype,
     *     ipol,ifreqhm,ispare1,anodata,vmin,vmax,
     *     SENSOR,TITLE,TYPE,TAG,CRPROC,CRTIME,
     *     DESCRIP,LDES,IAOPT,NIA)
        RETURN
        END
c
c
        subroutine readsirf_block(iu,ierr,nhead,nhtype,idatatype,stval,nsx,nsy,
     *       ioff,iscale,smin,smax,ncnt,anodata,vmin,vmax,lx,ly,nx,ny,krec)
c
c        read part (a rectangulatr block of data) of a SIR file image
c
c inputs:
c     iu       fortran logical unit of previously opened file
c     ierr     error code from header read
c     nhead    number of 512 byte blocks in header
c     nhtype   header type
c     nsx,nsy  image dimensions
c     ioff,iscale image storage scale factors
c     anodata   specify nodata value
c     vmin,vmax specify min, max valid range
c     lx,ly    lower-left corner of image block (0<lx<=nsx, 0<ly<=nsy)
c     nx,ny    upper-right corner of image block (lx<=nx<=nsx,ly<=ny<=nsy)
c     krec     record counter for most recent read (on first call set to nhead)
c
c outputs:
c     ierr   file read error code
c             set to 0 for successful read of file
c             set to -3 for data read error
c     stval  output image array (nx-lx+1 by ny-ly+1)
c             indexed as ind=(y-1)*(nx-lx+1)+x where 0 < y <= (ny-ly+1) and
c              0 < x <= (nx-lx+1) and 1 <= ind<= (nx-lx+1) * (ny-ly+1).
c     smin,max min and max of image within valid range 
c              valid range = [vmin..vmax] and .ne. anodata
c     ncnt   count of pixels in valid range
c     krec   record counter for most recent read
c
        integer                iu,nhead,nhtype
        real                stval(*)
        integer         nsx,nsy,ncnt
        integer                ioff,iscale
        real            anodata,vmin,vmax,smin,smax
        integer         lx,ly,nx,ny,krec
c
c     working variables
c
        byte            tempb(512)  ! integer*1 (-127 .. 128)
        integer*2        temp(256)
        real            tempf(128)
        integer         ntot,jj,irec
        real            s,soff
c
c       platform-dependent parameters
c
        include 'SIR.inc'
c
        if (ierr.lt.0) return   ! return if file open previously failed
c
c       check block specification
c
        if (lx.lt.1 .or.lx.gt.nsx.or.ly.lt.1 .or.ly.gt.nsy .or.
     $      nx.lt.lx.or.nx.gt.nsx.or.ny.lt.ly.or.ny.gt.nsy) then
           print *,'*** readsirf_block error: invalid block',lx,ly,nx,ny,nsx,nsy
           ierr=-1
           return
        endif
c
        ntot =(nx-lx+1)*(ny-ly+1)  ! total number of pixels to be read
        if (iscale.eq.0) iscale=1
        s=1./float(iscale)
        soff=32767.0/float(iscale)
        nrec=256                   ! pixel values per 512 byte record
        if (idatatype.eq.1) then
           soff=128.0/float(iscale)
           nrec=512
        endif
        if (idatatype.eq.4) nrec=128
c
        jj=0                       ! output array index
        i=((ly-1)*nsx+lx-1)        ! index of (lx,ly) in file
        irec=i/nrec                ! initial 512 byte record in image
        ii=i-irec*nrec+1           ! initial value offset in record
        irec=nhead+irec            ! first 512 byte record to read
        smin=1.e25
        smax=-1.e25
        ncnt=0
200        continue
            irec=irec+1
            if (VMS) then        ! sequential i/o
               if (krec.ge.irec) then
                  rewind(iu)
                  krec=0
               endif
               if (krec.lt.irec-1) then
                  do i=krec,irec-1
                     read (iu,414,end=299) temp
                  end do
                  krec=irec-1
               endif
               krec=krec+1
               if (idatatype.eq.1) then
                  read (iu,413,end=299) tempb
 413                  format(512a1)
               else if (idatatype.eq.4) then
                  read (iu,412,end=299) tempf
 412                  format(128a4)
               else
                  read (iu,414,end=299) temp
 414                  format(256a2)
                  call swapbuf(temp,256)
               endif
            else                ! direct record-based i/o
               if (idatatype.eq.1) then
                  read(iu,rec=irec,err=298) tempb
               else if (idatatype.eq.4) then
                  read(iu,rec=irec,err=298) tempf
               else
                  read(iu,rec=irec,err=298) temp
                  if (BSWAP) call swapbuf(temp,256)
               endif
            endif
            do i=ii,nrec
               jj=jj+1
               if (idatatype.eq.1) then
                  k=tempb(i)
                  stval(jj)=s*k+soff+ioff
               else if (idatatype.eq.4) then
                  stval(jj)=tempf(i)
               else
                  stval(jj)=s*temp(i)+soff+ioff
               endif
               if (stval(jj).ne.anodata.and.
     *                   stval(jj).ge.vmin.and.
     *                   stval(jj).le.vmax) then
                  smin=min(smin,stval(jj))
                  smax=max(smax,stval(jj))
                  ncnt=ncnt+1
               endif
               if (jj.ge.ntot) goto 299
            end do
            ii=1
            goto 200
c
298     continue      ! error return
        write (*,98) irec
98        format(' *** ERROR reading input SIR file at ',I8)
c
c        close file on read error
c
        close(unit=iu)
        ierr=-3
        return
c
299        continue      ! successful completion
c
c        file is NOT closed on successful return
c
        krec=irec
        ierr=0
        return
        end




c
c standard version 3.0 header SIR read data routine
c
c
        subroutine readsirf(iu,ierr,nhead,nhtype,idatatype,stval,nsx,nsy,
     *       ioff,iscale,smin,smax,ncnt,anodata,vmin,vmax)
c
c        read SIR file image data
c
c inputs:
c     iu       fortran logical unit of previously opened file
c     ierr     error code from header read
c     nhead    number of 512 byte blocks in header
c     nhtype   header type
c     nsx,nsy  image dimensions
c     ioff,iscale image storage scale factors
c     anodata   specify nodata value
c     vmin,vmax specify min, max valid range
c
c outputs:
c     ierr   file read error code
c             set to 0 for successful read of file
c             set to -3 for data read error
c     stval  output image array
c     smin,max min and max of image within valid range 
c              valid range = [vmin..vmax] and .ne. anodata
c     ncnt   count of pixels in valid range
c
        integer                iu,nhead,nhtype
        real                stval(*)
        integer         nsx,nsy,ncnt
        integer                ioff,iscale
        real            anodata,vmin,vmax,smin,smax
c
c     working variables
c
c        byte            tempb(512)  ! integer*1 (-127 .. 128)
        integer*1       tempb(512)  ! integer*1 (-127 .. 128)
        integer*2        temp(256)
        real            tempf(128)
        integer         ntot,jj,irec
        real            s,soff
c
c       platform-dependent parameters
c
        include 'SIR.inc'
c
        if (ierr.lt.0) return   ! return if file open previously failed
c
        ntot = nsx*nsy
        if (iscale.eq.0) iscale=1
        s=1./float(iscale)
        soff=32767.0/float(iscale)
        nrec=256
        if (idatatype.eq.1) then
           soff=128.0/float(iscale)
           nrec=512
        endif
        if (idatatype.eq.4) nrec=128
c
        jj=0
        irec=nhead
        smin=1.e25
        smax=-1.e25
        ncnt=0
200        continue
                irec=irec+1
                if (VMS) then
                   if (idatatype.eq.1) then
                      read (iu,413,end=299) tempb
 413                      format(512a1)
                   else if (idatatype.eq.4) then
                      read (iu,412,end=299) tempf
 412                      format(128a4)
                   else
                      read (iu,414,end=299) temp
 414                      format(256a2)
                      call swapbuf(temp,256)
                   endif
                else
                   if (idatatype.eq.1) then
                      read(iu,rec=irec,err=298) tempb
                   else if (idatatype.eq.4) then
                      read(iu,rec=irec,err=298) tempf
                   else
                      read(iu,rec=irec,err=298) temp
                      if (BSWAP) call swapbuf(temp,256)
                   endif
                endif
                do i=1,nrec
                        jj=jj+1
                        if (idatatype.eq.1) then
                           k=tempb(i)
                           stval(jj)=s*k+soff+ioff
                        else if (idatatype.eq.4) then
                           stval(jj)=tempf(i)
                        else
                           stval(jj)=s*temp(i)+soff+ioff
                        endif
                        if (stval(jj).ne.anodata.and.
     *                             stval(jj).ge.vmin.and.
     *                             stval(jj).le.vmax) then
                           smin=min(smin,stval(jj))
                           smax=max(smax,stval(jj))
                           ncnt=ncnt+1
                        endif
                        if (jj.ge.ntot) goto 299
                end do
                goto 200
298     continue
        write (*,98)
98        format(' *** ERROR reading input SIR file ***')
c
c        close file on read error
c
        close(unit=iu)
        ierr=-3
        return
299        continue
c
c        close file on successful completion
c
        close(unit=iu)
        ierr=0
        return
        end

C
C **************************************************************************
C
C STANDARD SIR IMAGE FORMAT READ/WRITE ROUTINES
C
C **************************************************************************
C
      subroutine readsirhead3(fname,iu_in,ierr,nhead,ndes,nhtype,
     *     idatatype,
     *     nsx,nsy,xdeg,ydeg,ascale,bscale,a0,b0,iopt,ioff,iscale,
     *     ixdeg_off,iydeg_off,ideg_sc,iscale_sc,ia0_off,ib0_off,
     *     i0_sc,
     *     iyear,isday,ismin,ieday,iemin,iregion,itype,
     *     ipol,ifreqhm,ispare1,anodata,vmin,vmax,
     *     sensor,title,type,tag,crproc,crtime,
     *     maxdes,descrip,ldes,maxi,iaopt,nia)
c
c     Version 3.0 standard SIR image format header read routine
c
c     written by D.G. Long  28 Oct 2000
c     revised to initialize ierr
c
c inputs:
c     fname    file name to read from
c     iu_in    fortran logical unit to use
c     maxdes   maximum number of characters in descrip
c     maxi     maximum number of integers in iaopt
c
c outputs:
c     ierr   file read error code
c             set to 0 for successful read of header
c             set to -1 for open error
c             set to -2 for header read error
c     nhead  number of 512 byte blocks in header
c     nhtype file header type code from header
c     idatatype file data type code (0,2=i*2,1=i*1,4=i*4) def=i*2
c     nsx,nsy image dimensions
c     xdeg..b0 projection factors
c     iopt    projection type
c     ioff,iscale image storage scale factors
c     ixdeg_off..i0_sc projection offset and scale storage factors (see below)
c     iyear,isday,ismin  year,day,minute of start of data used in image
c     ieday,iemin day,minute of end of data used in image
c     iregion study region id number
c     itype   image type code (=0 none or unknown)
c             typically: 1's digit is polarization (0=H, 1=V)
c             10's digit and higher is frequency in hundreds of MHz
c     ispare* spare integers
c     sensor  20 character string: describes sensor
c     title   80 character string: image title
c     type    80 character string: description of image type
c     tag     creator byline
c     crproc  100 character string: name and version of creating program
c     crtime  40 character string: date and time of creation
c     descrip variable length character string: description and annontation
c     ldes    number of characters in descrip
c     iaopt   optional integers
c       by convention, first value iaopt is a code telling how to interpret
c       the rest of the array if nia>0
c     nia     number of integers
c
c     SIR files have headers which are multiples of 512 bytes followed by
c     image data stored as 2 byte integers zero padded to be a multiple
c     of 512 bytes.
c
        implicit none
c
        character*(*)         fname
        integer                iu_in,ierr,nhead,ndes,nhtype,idatatype
        integer         nsx,nsy,iopt
        real                xdeg,ydeg,ascale,bscale,a0,b0
        integer                ioff,iscale,iyear,isday,ismin,ieday,iemin
        integer         ixdeg_off,iydeg_off,ideg_sc,iscale_sc
        integer         ia0_off,ib0_off,i0_sc,maxi,maxdes
        integer                iregion,itype,ipol,ifreqhm,ispare1,ldes,nia
        real            anodata,vmin,vmax
        character*(*)        sensor,title,type,tag,crproc,crtime,descrip
        integer                iaopt(*)
c
c     working variables
c
        integer*2       i2(2)
        real            f2
        equivalence (i2,f2)      ! used in i*2 to real mapping
        integer*2        temp(256)
        character        sens1*41,tit1*81,type1*139
        character       tag1*41,crp1*101,crt1*29
        integer         iu,minv,i,j,icnt,irec,iv
        real            soff

        integer*2       i256  ! this awkward definition avoids some compiler
        parameter (i256=256)  ! variable typing problems for the mod function
c
c        header mapping (see code for scaling of <= variables)
c       note that character strings do not have to be null terminated
c
c        temp(1)        = nsx                        ! pixels in x direction
c        temp(2) = nsy                        ! pixels in y direction
c       temp(3) <= xdeg                        ! span of x
c        temp(4) <= ydeg                        ! span of y
c        temp(5) = nhtype                ! header type (old<15,20,30)
c        temp(6) <= ascale               ! x scaling
c        temp(7) <= bscale               ! y scaling
c        temp(8) <= a0                   ! x (or lon) origin
c         note: longitudes should be in the range -180 to 180
c        temp(9) <= b0                   ! y (or lat) origin
c         note: latitudes should be in the range -90 to 90
c
c     scaling for prorjection parameters are generally:
c
c        temp(3) = nint((xdeg + float(ixdeg_off)) * float(ideg_sc))
c        temp(4) = nint((ydeg + float(iydeg_off)) * float(ideg_sc))
c        temp(6) = nint(ascale * float(iscale_sc))
c        temp(7) = nint(bscale * float(iscale_sc))
c        temp(8) = nint((a0 + float(ia0_off)) * float(i0_sc))
c        temp(9) = nint((b0 + float(ib0_off)) * float(i0_sc))
c
c     with the following projection specific exceptions:
c
c        if (iopt.eq.1.or.iopt.eq.2) then                ! lambert
c           temp(6) = nint(float(iscale_sc)/ascale)
c           temp(7) = nint(float(iscale_sc)/bscale)
c        if (iopt.eq.11.or.iopt.eq.12.or.iopt.eq.13) then ! EASE1 grid
c           temp(6) = nint(float(iscale_sc)*anint(10.*ascale*
c     *                         25.067525/6371.228)*0.05)
c           temp(7) = nint(float(iscale_sc)*anint(10.*bscale/
c     *                         25.067525)*0.05)
c
c        temp(10) = ioff                        ! offset to be added to scale val
c        temp(11) = iscale                ! scale factor ival=(val-ioff)/iscale
c        temp(12) = iyear                ! year for data used
c        temp(13) = isday                ! starting JD
c        temp(14) = ismin                ! time of day for first data (in min)
c        temp(15) = ieday                ! ending JD
c        temp(16) = iemin                ! time of day for last data (in min)
c        temp(17) = iopt                        ! projection type
c                                        !  -1 = no projection, image only
c                                        !   0 = rectalinear lat/lon
c                                        !   1 = lambert equal area
c                                        !   2 = lambert equal area (local rad)
c                                        !   5 = polar stereographic
c                                        !   8 = EASE2 north equal area grid
c                                        !   9 = EASE2 south equal area grid
c                                        !  10 = EASE2 cylindrical grid
c                                        !  11 = EASE1 north equal area grid
c                                        !  12 = EASE1 south equal area grid
c                                        !  13 = EASE1 cylindrical grid
c        temp(18) = iregion                ! region id code
c        temp(19) = itype                ! image type code
c                                       ! some standard values: (0=unknown)
c                                       ! 1 = scatterometer A (dB)
c                                       ! 2 = scatterometer B (dB/deg)
c                                       ! 3 = radiometer Tb (K)
c                                       ! 9 = topography (m)
c        temp(20)-temp(39) 40 chars of sensor
c       temp(40) = iscale_sc            ! ascale/bscale scale factor
c       temp(41) = nhead                ! number of 512 byte header blocks
c       temp(42) = ndes                 ! number of 512 byte blocks description
c       temp(43) = ldes                 ! number of bytes of description
c       temp(44) = nia                  ! number of optional integers
c       temp(45) = ipol                 ! polarization (0=n/a,1=H,2=V)
c       temp(46) = ifreqhm              ! frequency in 100's MHz (0 if n/a)
c       temp(47) = ispare1              ! spare
c       temp(48) = idatatype            ! data type code 0,2=i*2,1=i*1,4=float
c
c       the value of idata type determines how data is stored and how
c       anodata, vmin, and vmax are stored.
c
c       if idatatype = 1 data is stored as bytes (minv=128)
c       if idatatype = 2 data is stored as 2 byte integers (minv=32766)
c       if idatatype = 4 data is stored as IEEE floating point
c
c       if idatatype = 1,2 anodata,vmin,vmax are stored as 2 byte integers
c         in temp(49)..temp(51)  minv, ioff and iscal used to convert
c         integers or bytes into floating point values
c         nodata, vmin, and vmax must be representable with ioff and iscale
c            temp(*) = (value-ioff)*iscale-minv
c            value = float(temp(*)+minv)/float(iscale)+ioff
c       idatatype=2 is considered the SIR standard format
c
c       if idatatype = f anodata,vmin,vmax are stored as floating points
c         in temp(42)..temp(57) and minv, ioff and iscale are ignored here
c         and when reading the file.
c         floating point numbers are NOT standard across platforms and
c         are therefore not recommended
c
c       temp(49) <= anodata           ! value representing no data
c       temp(50) <= vmin              ! minimum useful value from creator prg
c       temp(51) <= vmax              ! maximum useful value from creator prg
c       temp(52,53) = anodata         ! IEEE floating value of no data
c       temp(54,55) = vmin            ! IEEE floating minimum useful value
c       temp(56,57) = vmax            ! IEEE floating maximum useful value
c
c        temp(58)-temp(126) 150 chars of type
c       temp(127) = ixdeg_off         ! xdeg offset
c       temp(128) = iydeg_off         ! ydeg offset
c        temp(129)-temp(168) 80 chars of title
c       temp(169) = ideg_sc           ! xdeg,ydeg scale factor
c        temp(170)-temp(189) 40 chars of tag
c       temp(190) = ia0_off           ! b0 offset 
c        temp(191)-temp(240) 100 chars of crproc
c       temp(241) = ib0_off           ! b0 offset 
c        temp(242)-temp(255) 28 chars of crtime
c       temp(256) = i0_sc             ! a0,b0 scale factor
c
c     optional header blocks:
c
c        ndes header blocks of 512 bytes: chars of description
c        nhead-ndes-1 header blocks of 512 bytes: values of iaopt
c       by convention, first value iaopt is a code telling how to interpret
c       the rest of the array if nia>0.  Usage of additional blocks is
c       user dependent and non-standard.
c
c       remainder of file is image data in a multiple of 512 byte blocks
c
c       one,two byte integer scaling (idatatype=1,2) is
c          intval = (fvalue-ioff)*iscale-minv
c          fvalue = float(intval+minv)/float(iscale)+ioff
c       no scaling of float values for (idatatype=4)
c
c
c *******************************************************************
c
c     get platform-dependent parameters
c
        include 'SIR.inc'
c
c *******************************************************************
c
c     error check other input values
c
        iu=iu_in
        if (iu.lt.1) then
           iu=2
           write (*,*) '*** readsirhead3 error: invalid iu, using ',iu
        endif
        ierr=0
c
c        open input file
c
        if (VMS) then
          open(unit=iu,file=fname,status='old',form='formatted',
     *                recl=ISIRRECLENGTH,err=99) !,
c     *          readonly,recordtype='fixed',carriagecontrol='none')
        else
          open(unit=iu,file=fname,status='old',form='unformatted',
     *          recl=ISIRRECLENGTH,access='direct',err=99)
        endif
c
c        read first header block
c
        if (VMS) then
                read (iu,414) temp
 414                format(256a2)
                call swapbuf(temp,256)
        else
                read (iu,rec=1) temp
                if (BSWAP) call swapbuf(temp,256)
        endif
c
c        decode first header block
c
        nsx = temp(1)                        ! dimension
        nsy = temp(2)
        iopt = temp(17)                        ! transformation option
        nhtype = temp(5)                ! header type
        if (nhtype.lt.20) nhtype=1
c
        if (nhtype.lt.30) then
c
c     set version 3.0 parameters to header version 2.0 defaults
c     
           if (iopt.eq.-1) then ! image only
              ideg_sc=10
              iscale_sc=1000
              i0_sc=100
              ixdeg_off=0
              iydeg_off=0
              ia0_off=0
              ib0_off=0
           else if (iopt.eq.0) then ! rectalinear lat/lon
              ideg_sc=100
              iscale_sc=1000
              i0_sc=100
              ixdeg_off=-100
              iydeg_off=0
              ia0_off=0
              ib0_off=0
           else if (iopt.eq.1.or.iopt.eq.2) then ! lambert
              ideg_sc=100
              iscale_sc=1000
              i0_sc=1
              ixdeg_off=0
              iydeg_off=0
              ia0_off=0
              ib0_off=0
           else if (iopt.eq.5) then ! polar stereographic
              ideg_sc=100
              iscale_sc=100
              i0_sc=1
              ixdeg_off=-100
              iydeg_off=0
              ia0_off=0
              ib0_off=0
           else if (iopt.eq.8.or.iopt.eq.9.or.iopt.eq.10) then ! EASE2 grid
              ideg_sc=10
              iscale_sc=100
              i0_sc=10
              ixdeg_off=0
              iydeg_off=0
              ia0_off=0
              ib0_off=0
           else if (iopt.eq.11.or.iopt.eq.12.or.iopt.eq.13) then ! EASE1 grid
              ideg_sc=10
              iscale_sc=1000
              i0_sc=10
              ixdeg_off=0
              iydeg_off=0
              ia0_off=0
              ib0_off=0
           else                 ! unknown default scaling
              ideg_sc=100
              iscale_sc=1000
              i0_sc=100
              ixdeg_off=0
              iydeg_off=0
              ia0_off=0
              ib0_off=0
           endif
        else
c
c     get projection parameters offset and scale factors
c
           iscale_sc = temp( 40)
           ixdeg_off = temp(127)
           iydeg_off = temp(128)
           ideg_sc   = temp(169)
           ia0_off   = temp(190)
           ib0_off   = temp(241)
           i0_sc     = temp(256)
        endif
c
c     decode default projection parameters
c
        xdeg   = temp(3)/float(ideg_sc) - ixdeg_off
        ydeg   = temp(4)/float(ideg_sc) - iydeg_off
        ascale = temp(6)/float(iscale_sc)
        bscale = temp(7)/float(iscale_sc)
        a0     = temp(8)/float(i0_sc) - ia0_off
        b0     = temp(9)/float(i0_sc) - ib0_off
c
c        handle special cases which depend on transformation option
c
        if (iopt.eq.-1) then                                ! image only
        else if (iopt.eq.0) then                        ! rectalinear
        else if (iopt.eq.1.or.iopt.eq.2) then                ! lambert
                ascale=float(iscale_sc)/float(temp(6))
                bscale=float(iscale_sc)/float(temp(7))
        else if (iopt.eq.5) then                        ! polar stereographic
        else if (iopt.eq.8.or.iopt.eq.9.or.iopt.eq.10) then ! EASE2 grid
        else if (iopt.eq.11.or.iopt.eq.12.or.iopt.eq.13) then ! EASE1 grid
                ascale = 2.0d0*(float(temp(6))/float(iscale_sc))
     *               *6371.228d0/25.067525d0
                bscale = 2.0d0*(float(temp(7))/float(iscale_sc))
     *               *25.067525d0
        else
                write (*,*)'*** Unrecognized SIR option in READSIRHEAD ***'
        endif
c
c        decode scaling, data time, region, and type information
c
        ioff = temp(10)
        iscale = temp(11)
        if (iscale.eq.0) iscale=1
        iyear = temp(12)
        isday = temp(13)
        ismin = temp(14)
        ieday = temp(15)
        iemin = temp(16)
        iregion = temp(18)
        itype = temp(19)
        nhead = temp(41)        ! total number of header blocks
        if (nhead.eq.0) nhead=1
        if (nhtype.eq.1) nhead=1
        ndes = temp(42)         ! number of description blocks
        ldes = temp(43)
        nia = temp(44)
        ipol = temp(45)
        ifreqhm = temp(46)
        ispare1 = temp(47)
        idatatype = temp(48)    ! data type code
        if (idatatype.eq.0) idatatype=2
        if (nhtype.eq.1) then
           ndes=0
           ldes=0
           nia=0
           idatatype=2
        endif
c
        if (idatatype.eq.1) then
           minv=128
        else
           minv=32767
        endif
c
        soff=float(minv)/float(iscale)
        if (idatatype.eq.1) soff=float(minv)/float(iscale)
c
        anodata = temp(49)/float(iscale)+ioff+soff  ! int type values
        vmin = temp(50)/float(iscale)+ioff+soff
        vmax = temp(51)/float(iscale)+ioff+soff
c
        if (idatatype.eq.4) then  ! float type values
           i2(1)=temp(52)
           i2(2)=temp(53)
           anodata=f2
           i2(1)=temp(54)
           i2(2)=temp(55)
           vmin=f2
           i2(1)=temp(56)
           i2(2)=temp(57)
           vmax=f2
        endif
c
c        get character variables
c
        do i=1,20
                j=(i-1)*2+1
                sens1(j:j)=char(mod(temp(i+19),i256))
                sens1(j+1:j+1)=char(temp(i+19)/256)
        end do
        sens1(41:41)=char(0)  ! null terminate string
        sensor=sens1
        do i=1,40
                j=(i-1)*2+1
                tit1(j:j)=char(mod(temp(i+128),i256))
                tit1(j+1:j+1)=char(temp(i+128)/256)
        end do
        tit1(81:81)=char(0)  ! null terminate string
        title=tit1
        if (nhtype.gt.1) then
           do i=1,69
              j=(i-1)*2+1
              type1(j:j)=char(mod(temp(i+57),i256))
              type1(j+1:j+1)=char(temp(i+57)/256)
           end do
           type1(139:139)=char(0)  ! null terminate string
           do i=1,20
              j=(i-1)*2+1
              tag1(j:j)=char(mod(temp(i+169),i256))
              tag1(j+1:j+1)=char(temp(i+169)/256)
           end do
           tag1(41:41)=char(0)  ! null terminate string
           do i=1,50
              j=(i-1)*2+1
              crp1(j:j)=char(mod(temp(i+190),i256))
              crp1(j+1:j+1)=char(temp(i+190)/256)
           end do
           crp1(101:101)=char(0)  ! null terminate string
           do i=1,14
              j=(i-1)*2+1
              crt1(j:j)=char(mod(temp(i+241),i256))
              crt1(j+1:j+1)=char(temp(i+241)/256)
           end do
           crt1(29:29)=char(0)  ! null terminate string
        else
           type1=' '
           tag1=' '
           crp1=' '
           crt1=' '
           nhead=1
           nhtype=20  ! upgrade header version upon return
        endif
        type=type1
        tag=tag1
        crproc=crp1
        crtime=crt1
        descrip=' '
c
c        get extra storage information
c
        if (nhead.gt.1) then
           if (ndes.gt.0) then
c
c        read next header blocks
c
              icnt=0
              do irec=2,ndes+1
                 if (VMS) then
                    read (iu,414) temp
                    call swapbuf(temp,256)
                 else
                    read (iu,rec=irec) temp
                    if (BSWAP) call swapbuf(temp,256)
                 endif
                 do j=1,512
                    icnt=icnt+1
                    if (icnt.le.ldes.and.icnt.le.maxdes) then
                       i=(j-1)/2+1
                       iv=temp(i)
                       if (mod(j,2).eq.0) then
                          iv=iv/256
                       else
                          iv=mod(iv,256)
                       endif
                       descrip(icnt:icnt)=char(iv)
                    endif
                 end do                
              end do
              if (maxdes.lt.ldes) then
                 write (*,*) '*** readsirhead3 warning: extra char',
     *                ' head overflow ',nia,maxi
                 ldes=maxdes
              endif
           endif
           if (nhead-ndes-1.gt.0) then
              icnt=0
              do irec=ndes+2,nhead
                 if (VMS) then
                    read (iu,414) temp
                    call swapbuf(temp,256)
                 else
                    read (iu,rec=irec) temp
                    if (BSWAP) call swapbuf(temp,256)
                 endif
                 do i=1,256
                    icnt=icnt+1
                    if (icnt.le.nia.and.icnt.le.maxi) 
     *                 iaopt(icnt) = temp(i)
                 end do
              end do
              if (nia.gt.maxi) then
                 write (*,*) '*** readsirhead3 warning: extra ints',
     *                ' head overflow ',nia,maxi
                 nia=maxi
              endif
           endif
        endif
        return
c     
 300    continue
        write (*,398)
398        format(' *** readsirhead3 ERROR: read header of SIR file ***')
c
c        close file on header read error
c
        close(unit=iu)
c
        ierr=-2
        return

99        continue
        write (*,198)
198        format(' *** readsirhead3 ERROR: opening input SIR file ***')
        ierr=-1
        return
        end
C
C **************************************************************************
C
C header version 2.0 interface to version 3.0 standard sir read routine
C
C **************************************************************************
C
      subroutine readsirhead(fname,iu_in,ierr,nhead,ndes,nhtype,
     *     idatatype,
     *     nsx,nsy,xdeg,ydeg,ascale,bscale,a0,b0,iopt,ioff,iscale,
     *     iyear,isday,ismin,ieday,iemin,iregion,itype,
     *     ipol,ifreqhm,ispare1,anodata,vmin,vmax,
     *     sensor,title,type,tag,crproc,crtime,
     *     maxdes,descrip,ldes,maxi,iaopt,nia)
c
        implicit none
c
        character*(*)         fname
        integer                iu_in,ierr,nhead,ndes,nhtype,idatatype
        integer         nsx,nsy,iopt
        real                xdeg,ydeg,ascale,bscale,a0,b0
        integer                ioff,iscale,iyear,isday,ismin,ieday,iemin
        integer         ixdeg_off,iydeg_off,ideg_sc,iscale_sc
        integer         ia0_off,ib0_off,i0_sc,maxi,maxdes
        integer                iregion,itype,ipol,ifreqhm,ispare1,ldes,nia
        real            anodata,vmin,vmax
        character*(*)        sensor,title,type,tag,crproc,crtime,descrip
        integer                iaopt(*)
c
c     call version 3.0 read sir routine
c
        call readsirhead3(fname,iu_in,ierr,nhead,ndes,nhtype,
     *       idatatype,
     *       nsx,nsy,xdeg,ydeg,ascale,bscale,a0,b0,iopt,ioff,iscale,
     *       ixdeg_off,iydeg_off,ideg_sc,iscale_sc,ia0_off,ib0_off,
     *       i0_sc,
     *       iyear,isday,ismin,ieday,iemin,iregion,itype,
     *       ipol,ifreqhm,ispare1,anodata,vmin,vmax,
     *       sensor,title,type,tag,crproc,crtime,
     *       maxdes,descrip,ldes,maxi,iaopt,nia)
c
        if (nhtype.gt.30) then
           write (*,*) '*** readsirhead warning: file has version 3.0',
     *          ' head'
           write (*,*) '*** use readsirhead3 instead'
        endif
c
        return
        end

        SUBROUTINE SIRUPDATE(STVAL,NHTYPE,NSX,NSY,ANODATA,VMIN,VMAX)
	dimension STVAL(*)
C
C     UPDATE SIR HEADER INFO FOR OLD STYLE HEADERS
C
        IF (NHTYPE.EQ.1) THEN   ! OLD STYLE HEADER
C
C        DETERMINE MAX/MIN OF DATA ARRAY
C
           SMIN=1.E25
           SMAX=-1.E25
           DO NN1=0,NSY-1
              DO NN2=1,NSX
                 INDEX = NN1*NSX + NN2
                 IF (STVAL(INDEX).LT.SMIN) SMIN=STVAL(INDEX)
                 IF (STVAL(INDEX).GT.SMAX) SMAX=STVAL(INDEX)
              END DO
           END DO
           ANODATA=SMIN
           VMIN=SMIN
           VMAX=SMAX
           IF (ABS(SMIN+32.0).LT.2) THEN
              ANODATA=-32.0
              VMIN=-32.0
           ELSE IF (ABS(SMIN+3.2).LT.2) THEN
              ANODATA=-3.2
              VMIN=-3.0
           ENDIF
           IF (ABS(SMAX).LT.2) THEN
              VMAX=0.0
           ENDIF
        ENDIF
        RETURN
        END
c
c
        subroutine swapbuf(t,n)
c
c       swaps bytes for integer*2
c
        logical*1 t(1),s
        do i=1,n
                j=(i-1)*2+1
                s=t(j)
                t(j)=t(j+1)
                t(j+1)=s
        end do
        return
        end
c
c *************************************************************************
c
c  standard SIR file write routine
c
      subroutine writesir3(fname,iu_in,ierr,nhead,ndes,nhtype,idatatype,
     *     nsx,nsy,xdeg,ydeg,ascale,bscale,a0,b0,iopt,ioff,iscale,
     *     ixdeg_off,iydeg_off,ideg_sc,iscale_sc,ia0_off,ib0_off,i0_sc,
     *     iyear,isday,ismin,ieday,iemin,iregion,itype,
     *     ipol,ifreqhm,ispare1,anodata,vmin,vmax,
     *     sensor,title,type,tag,crproc,crtime,
     *     stval,descrip,ldes,iaopt,nia)
c
c     Version 3 standard SIR image format write routine
c
c     written by D.G. Long  28 Oct 2000
c     Revised by D.G. Long  12 Mar 2014 + added EASE2
c
c
c inputs:
c     fname    file name to write to
c              note: may overwrite an existing file of the same name
c     iu_in    fortran logical unit to use
c     idatatype file data type code (0,2=i*2,1=i*1,4=i*4) def=i*2
c     nsx,nsy image dimensions
c     xdeg..b0 projection factors
c     iopt    projection type
c     ioff,iscale image storage scale factors
c     ixdeg_off..i0_sc projection offset and scale storage factors (see below)
c     iyear,isday,ismin  year,day,minute of start of data used in image
c     ieday,iemin day,minute of end of data used in image
c     iregion study region id number
c     itype   image type code (=0 none or unknown)
c     ispare1 spare integer
c     sensor  20 character string: describes sensor
c     title   80 character string: image title
c     type    80 character string: description of image type
c     tag     creator byline
c     crproc  100 character string: name and version of creating program
c     crtime  40 character string: date and time of creation
c     stval   image array (floats)
c     descrip variable length character string: description and annontation
c     ldes    number of characters in descrip
c     iaopt   array of optional integers
c       by convention, first value iaopt is a code telling how to interpret
c       the rest of the array if nia>0
c     nia     number of integers in extra iaopt header
c
c outputs:
c     ierr   file write error code
c             set to 0 for successful write
c             set to -1 for open error
c             set to -2 for write error
c             set to -3 for invalid image size
c     nhead  number of 512 byte blocks in header
c     nhtype file header type code
c
c     SIR files have headers which are multiples of 512 bytes followed by
c     image data stored as 2 byte integers zero padded to be a multiple
c     of 512 bytes.
c
        implicit none
c
        character*(*)         fname
        integer                iu_in,ierr,nhead,ndes,nhtype,idatatype
        integer         nsx,nsy,iopt
        real                xdeg,ydeg,ascale,bscale,a0,b0
        integer         ixdeg_off,iydeg_off,ideg_sc,iscale_sc
        integer         ia0_off,ib0_off,i0_sc
        integer                ioff,iscale,iyear,isday,ismin,ieday,iemin
        integer                iregion,itype,ipol,ifreqhm,ispare1,ldes,nia
        real            anodata,vmin,vmax
        real            stval(*)
        character*(*)        sensor,title,type,tag,crproc,crtime,descrip
        integer                iaopt(*)
c
c     working variables
c
c        byte            tempb(512)    ! integer*1 (-127 .. 128)
        integer*1       tempb(512)    ! integer*1 (-127 .. 128)
        integer*2        temp(256)
        real            tempf(128)
        equivalence     (temp,tempb)  ! used only to save space
        equivalence     (temp,tempf)  ! used only to save space
        integer*2       i2(2)
        real            f2
        equivalence     (f2,i2)       ! used in i*2 to float mapping
        character        sens1*41,tit1*81,type1*139
        character       tag1*41,crp1*101,crt1*29
        integer         iu,minv,maxv,i,ii,j,k,ntot,nrec
        integer         irec,icnt,imax,iocnt,iucnt
c
c        header mapping (see code for scaling of <= variables)
c       note that character strings do not have to be null terminated
c
c        temp(1)        = nsx                        ! pixels in x direction
c        temp(2) = nsy                        ! pixels in y direction
c       temp(3) <= xdeg                        ! span of x
c        temp(4) <= ydeg                        ! span of y
c        temp(5) = nhtype                ! header type (old<15,20,30)
c        temp(6) <= ascale               ! x scaling
c        temp(7) <= bscale               ! y scaling
c        temp(8) <= a0                   ! x (or lon) origin
c         note: longitudes should be in the range -180 to 180
c        temp(9) <= b0                   ! y (or lat) origin
c         note: latitudes should be in the range -90 to 90
c
c     scaling for prorjection parameters are generally:
c
c        temp(3) = nint((xdeg + float(ixdeg_off)) * float(ideg_sc))
c        temp(4) = nint((ydeg + float(iydeg_off)) * float(ideg_sc))
c        temp(6) = nint(ascale * float(iscale_sc))
c        temp(7) = nint(bscale * float(iscale_sc))
c        temp(8) = nint((a0 + float(ia0_off)) * float(i0_sc))
c        temp(9) = nint((b0 + float(ib0_off)) * float(i0_sc))
c
c     with the following projection specific exceptions:
c
c        if (iopt.eq.1.or.iopt.eq.2) then                ! lambert
c           temp(6) = nint(float(iscale_sc)/ascale)
c           temp(7) = nint(float(iscale_sc)/bscale)
c        if (iopt.eq.11.or.iopt.eq.12.or.iopt.eq.13) then ! EASE1 grid
c           temp(6) = nint(float(iscale_sc)*anint(10.*ascale*
c     *                         25.067525/6371.228)*0.05)
c           temp(7) = nint(float(iscale_sc)*anint(10.*bscale/
c     *                         25.067525)*0.05)
c
c        temp(10) = ioff                        ! offset to be added to scale val
c        temp(11) = iscale                ! scale factor ival=(val-ioff)/iscale
c        temp(12) = iyear                ! year for data used
c        temp(13) = isday                ! starting JD
c        temp(14) = ismin                ! time of day for first data (in min)
c        temp(15) = ieday                ! ending JD
c        temp(16) = iemin                ! time of day for last data (in min)
c        temp(17) = iopt                        ! projection type
c                                        !  -1 = no projection, image only
c                                        !   0 = rectalinear lat/lon
c                                        !   1 = lambert equal area
c                                        !   2 = lambert equal area (local rad)
c                                        !   5 = polar stereographic
c                                        !   8 = EASE2 north equal area grid
c                                        !   9 = EASE2 south equal area grid
c                                        !  10 = EASE2 cylindrical grid
c                                        !  11 = EASE1 north equal area grid
c                                        !  12 = EASE1 south equal area grid
c                                        !  13 = EASE1 cylindrical grid
c        temp(18) = iregion                ! region id code
c        temp(19) = itype                ! image type code
c                                       ! some standard values: (0=unknown)
c                                       ! 1 = scatterometer A (dB)
c                                       ! 2 = scatterometer B (dB/deg)
c                                       ! 3 = radiometer Tb (K)
c                                       ! 9 = topography (m)
c        temp(20)-temp(39) 40 chars of sensor
c       temp(40) = iscale_sc            ! ascale/bscale scale factor
c       temp(41) = nhead                ! number of 512 byte header blocks
c       temp(42) = ndes                 ! number of 512 byte blocks description
c       temp(43) = ldes                 ! number of bytes of description
c       temp(44) = nia                  ! number of optional integers
c       temp(45) = ipol                 ! polarization (0=n/a,1=H,2=V)
c       temp(46) = ifreqhm              ! frequency in 100's MHz (0 if n/a)
c       temp(47) = ispare1              ! spare
c       temp(48) = idatatype            ! data type code 0,2=i*2,1=i*1,4=float
c
c       the value of idata type determines how data is stored and how
c       anodata, vmin, and vmax are stored.
c
c       if idatatype = 1 data is stored as bytes (minv=128)
c       if idatatype = 2 data is stored as 2 byte integers (minv=32766)
c       if idatatype = 4 data is stored as IEEE floating point
c
c       if idatatype = 1,2 anodata,vmin,vmax are stored as 2 byte integers
c         in temp(49)..temp(51)  minv, ioff and iscal used to convert
c         integers or bytes into floating point values
c         nodata, vmin, and vmax must be representable with ioff and iscale
c            temp(*) = (value-ioff)*iscale-minv
c            value = float(temp(*)+minv)/float(iscale)+ioff
c       idatatype=2 is considered the SIR standard format
c
c       if idatatype = f anodata,vmin,vmax are stored as floating points
c         in temp(42)..temp(57) and minv, ioff and iscale are ignored here
c         and when reading the file.
c         floating point numbers are NOT standard across platforms and
c         are therefore not recommended
c
c       temp(49) <= anodata           ! value representing no data
c       temp(50) <= vmin              ! minimum useful value from creator prg
c       temp(51) <= vmax              ! maximum useful value from creator prg
c       temp(52,53) = anodata         ! IEEE floating value of no data
c       temp(54,55) = vmin            ! IEEE floating minimum useful value
c       temp(56,57) = vmax            ! IEEE floating maximum useful value
c
c        temp(58)-temp(126) 150 chars of type
c       temp(127) = ixdeg_off         ! xdeg offset
c       temp(128) = iydeg_off         ! ydeg offset
c        temp(129)-temp(168) 80 chars of title
c       temp(169) = ideg_sc           ! xdeg,ydeg scale factor
c        temp(170)-temp(189) 40 chars of tag
c       temp(190) = ia0_off           ! b0 offset 
c        temp(191)-temp(240) 100 chars of crproc
c       temp(241) = ib0_off           ! b0 offset 
c        temp(242)-temp(255) 28 chars of crtime
c       temp(256) = i0_sc             ! a0,b0 scale factor
c
c     optional header blocks:
c
c        ndes header blocks of 512 bytes: chars of description
c        nhead-ndes-1 header blocks of 512 bytes: values of iaopt
c       by convention, first value iaopt is a code telling how to interpret
c       the rest of the array if nia>0.  Usage of additional blocks is
c       user dependent and non-standard.
c
c       remainder of file is image data in a multiple of 512 byte blocks
c
c       one,two byte integer scaling (idatatype=1,2) is
c          intval = (fvalue-ioff)*iscale-minv
c          fvalue = float(intval+minv)/float(iscale)+ioff
c       no scaling of float values for (idatatype=4)
c
c *******************************************************************
c
c     get platform-dependent parameters
c
        include 'SIR.inc'
c
c *******************************************************************
c
c     error check other input values
c
        iu=iu_in
        if (iu.lt.1) then
           iu=2
           write (*,*) '*** writesir3 error: invalid iu, using ',iu
        endif
c
        if (nsx.lt.1.or.nsy.lt.1.or.nsx.gt.32767
     *       .or.nsy.gt.32767) then
           write (*,*) '*** writesir3 error: invalid nsx,nsy ',nsx,nsy
           ierr=-3
           return
        endif
c
c     check projection scale and offset parameters
c
        if (nhtype.lt.30.or.
     *      ideg_sc.le.0.or.iscale_sc.le.0.or.i0_sc.le.0) then
           if (nhtype.ge.30.and.
     *          (ideg_sc.le.0.or.iscale_sc.le.0.or.i0_sc.le.0))
     *       write(*,*) '*** writsir3 warning: projection scale error'
     *          ,ideg_sc,iscale_sc,i0_sc,ixdeg_off,iydeg_off,
     *          ia0_off,ib0_off
c
c     set parameters to header version 2.0 defaults
c
           if (iopt.eq.-1) then ! image only
              ideg_sc=10
              iscale_sc=1000
              i0_sc=100
              ixdeg_off=0
              iydeg_off=0
              ia0_off=0
              ib0_off=0
           else if (iopt.eq.0) then ! rectalinear lat/lon
              ideg_sc=100
              iscale_sc=1000
              i0_sc=100
              ixdeg_off=-100
              iydeg_off=0
              ia0_off=0
              ib0_off=0
           else if (iopt.eq.1.or.iopt.eq.2) then ! lambert
              ideg_sc=100
              iscale_sc=1000
              i0_sc=1
              ixdeg_off=0
              iydeg_off=0
              ia0_off=0
              ib0_off=0
           else if (iopt.eq.5) then ! polar stereographic
              ideg_sc=100
              iscale_sc=100
              i0_sc=1
              ixdeg_off=-100
              iydeg_off=0
              ia0_off=0
              ib0_off=0
           else if (iopt.eq.8.or.iopt.eq.9.or.iopt.eq.10) then ! EASE2 grid
              ideg_sc=10
              iscale_sc=100
              i0_sc=10
              ixdeg_off=0
              iydeg_off=0
              ia0_off=0
              ib0_off=0
           else if (iopt.eq.11.or.iopt.eq.12.or.iopt.eq.13) then ! EASE1 grid
              ideg_sc=10
              iscale_sc=1000
              i0_sc=10
              ixdeg_off=0
              iydeg_off=0
              ia0_off=0
              ib0_off=0
           else ! unknown default scaling
              ideg_sc=100
              iscale_sc=1000
              i0_sc=100
              ixdeg_off=0
              iydeg_off=0
              ia0_off=0
              ib0_off=0
           endif
           if (nhtype.lt.30) nhtype=30  ! can only write new headers
c           if (nhtype.ge.30.and.
c     *          (ideg_sc.le.0.0.or.iscale_sc.le.0.0.or.i0_sc.le.0)
c     *       write(*,*) '*** values used:',ideg_sc,iscale_sc,i0_sc,
c     *          ixdeg_off,iydeg_off,ia0_off,ib0_off
        endif
c
c        make first header block
c
        ierr = 0
        temp(1)        = nsx                                ! image dimensions
        temp(2) = nsy
        if (nhtype.lt.30) nhtype = 31 ! nhtype=30 reserved for version 2.0
        temp(5) = nhtype                       ! version 3.0 header code
c
c       scale projection variables based on transform option
c
        temp(3) = nint((xdeg + float(ixdeg_off)) * float(ideg_sc))
        temp(4) = nint((ydeg + float(iydeg_off)) * float(ideg_sc))
        temp(6) = nint(ascale * float(iscale_sc))
        temp(7) = nint(bscale * float(iscale_sc))
        temp(8) = nint((a0 + float(ia0_off)) * float(i0_sc))
        temp(9) = nint((b0 + float(ib0_off)) * float(i0_sc))
c
        if (iopt.eq.-1)        then                        ! image only
        else if (iopt.eq.0) then                ! rectalinear lat/lon
        else if (iopt.eq.1.or.iopt.eq.2) then                ! lambert
           temp(6) = nint(float(iscale_sc)/ascale)
           temp(7) = nint(float(iscale_sc)/bscale)
        else if (iopt.eq.5) then                ! polar stereographic
        else if (iopt.eq.8.or.iopt.eq.9.or.iopt.eq.10) then ! EASE2 grid
        else if (iopt.eq.11.or.iopt.eq.12.or.iopt.eq.13) then ! EASE1 grid
           temp(6) = nint(float(iscale_sc)*anint(10.*ascale*
     *                         25.067525/6371.228)*0.05)
           temp(7) = nint(float(iscale_sc)*anint(10.*bscale/
     *                         25.067525)*0.05)
        else
           write (*,*) '*** Unrecognized SIR option in writesir3 ***'
        endif
c
c     store projection parameters offset and scale factors
c
        temp( 40) = iscale_sc
        temp(127) = ixdeg_off
        temp(128) = iydeg_off
        temp(169) = ideg_sc
        temp(190) = ia0_off
        temp(241) = ib0_off
        temp(256) = i0_sc
c
c        store scale factors, time, option, region, and image type
c
        temp(10) = ioff
        temp(11) = iscale
        if (iscale.eq.0) temp(11)=1
        temp(12) = iyear
        temp(13) = isday
        temp(14) = ismin
        temp(15) = ieday
        temp(16) = iemin
        temp(17) = iopt
        temp(18) = iregion
        temp(19) = itype
c
        nhead=1
        ndes=0
        if (ldes.gt.0) then
           ndes=ldes/512
           if (mod(ldes,512).gt.0) ndes=ndes+1
c           write (*,*) 'extra descriptor',ldes,ndes
           nhead=nhead+ndes
        endif
        if (nia.gt.0) then
c           write (*,*) 'extra integers',nia
           nhead=nhead+nia/256
           if (mod(nia,256).gt.0) nhead=nhead+1
        endif
c
        temp(41) = nhead
        temp(42) = ndes
        temp(43) = ldes
        temp(44) = nia
        temp(45) = ipol
        temp(46) = ifreqhm
        temp(47) = ispare1
        temp(48) = idatatype
c
        if (idatatype.eq.1) then
           maxv=127
           minv=128
        else
           maxv=32767
           minv=32767
        endif
c
        if (idatatype.ne.4) then
           k = nint(iscale*(anodata-ioff))-minv
           if (k.gt.maxv) then
              k=maxv
              write (*,*) '*** WRITESIR: overflow in "anodata"',
     *             anodata,ioff,iscale
           endif
           if (k.lt.-minv) then
              k=-minv
              write (*,*) '*** WRITESIR: underflow in "anodata"',
     *             anodata,ioff,iscale
           endif
           temp(49) = k
           k = nint(iscale*(vmin-ioff))-minv
           if (k.gt.maxv) then
              k=maxv
              write (*,*) '*** WRITESIR: overflow in "vmin"',
     *             vmin,ioff,iscale
           endif
           if (k.lt.-minv) then
              k=-minv
              write (*,*) '*** WRITESIR: underflow in "vmin"',
     *             vmin,ioff,iscale
           endif
           temp(50) = k
           k = nint(iscale*(vmax-ioff))-minv
           if (k.gt.maxv) then
              k=maxv
              write (*,*) '*** WRITESIR: overflow in "vmax"',
     *             vmax,ioff,iscale
           endif
           if (k.lt.-minv) then
              k=-minv
              write (*,*) '*** WRITESIR: underflow in "vmax"',
     *             vmax,ioff,iscale
           endif
           temp(51) = k
        endif
c
        f2=anodata
        temp(52)=i2(1)
        temp(53)=i2(2)
        f2=vmin
        temp(54)=i2(1)
        temp(55)=i2(2)
        f2=vmax
        temp(56)=i2(1)
        temp(57)=i2(2)
c
c        write character variables
c
        sens1=' '  ! initially make them blank
        tit1 =' '
        type1=' '
        tag1 =' '
        crp1 =' '
        crt1 =' '
        sens1=sensor
        do i=1,20
                j=(i-1)*2+1
                temp(i+19)=ichar(sens1(j:j))+256*ichar(sens1(j+1:j+1))
        end do

        type1=type
        do i=1,69
                j=(i-1)*2+1
                temp(i+57)=ichar(type1(j:j))+256*ichar(type1(j+1:j+1))
        end do

        tit1=title
        do i=1,40
                j=(i-1)*2+1
                temp(i+128)=ichar(tit1(j:j))+256*ichar(tit1(j+1:j+1))
        end do

        tag1=tag
        do i=1,20
                j=(i-1)*2+1
                temp(i+169)=ichar(tag1(j:j))+256*ichar(tag1(j+1:j+1))
        end do

        crp1=crproc
        do i=1,50
                j=(i-1)*2+1
                temp(i+190)=ichar(crp1(j:j))+256*ichar(crp1(j+1:j+1))
        end do

        crt1=crtime
        do i=1,14
                j=(i-1)*2+1
                temp(i+241)=ichar(crt1(j:j))+256*ichar(crt1(j+1:j+1))
        end do
c
c        open output file
c
        if (VMS) then
          open(unit=iu,file=fname,status='new',      ! for vax/vms
     *          form='formatted',err=999,
     *                recl=ISIRRECLENGTH)
c     *          ,recordtype='fixed',carriagecontrol='none')
        else
          open(unit=iu,file=fname,status='unknown',  ! generally portable
     *          form='unformatted',err=999,          ! unix file open
     *                recl=ISIRRECLENGTH,access='direct')
        endif
c
c        write first header block
c
        irec=1
        if (BSWAP) call swapbuf(temp,256)
        if (VMS) then
          write(iu,414,err=1999) temp
        else
          write(iu,rec=1,err=1999) temp
        endif
c
c       write additional header records if needed
c
        if (ndes.gt.0) then
           icnt=0
           do ii=1,ndes
              do i=1,256
                 k=0
                 icnt=icnt+1
                 if (icnt.le.ldes) then
                    k=ichar(descrip(icnt:icnt))
                    icnt=icnt+1
                    if (icnt.le.ldes) k=k+256*ichar(descrip(icnt:icnt))
                 endif
                 temp(i)=k
              end do
              if (BSWAP) call swapbuf(temp,256)
              irec=irec+1
              if (VMS) then
                 write(iu,414,err=1999) temp
              else
                 write(iu,rec=irec,err=1999) temp
              endif
           end do
        endif
c
c        write extra header parameters
c
        if (nhead-ndes-1.gt.0) then
           icnt=0
           do ii=nhead-ndes+1,nhead
              do i=1,256
                 icnt=icnt+1
                 if (icnt.le.nia) then
                    temp(i)=iaopt(icnt)
                 else
                    temp(i)=0
                 endif
              end do
              if (BSWAP) call swapbuf(temp,256)
              irec=irec+1
              if (VMS) then
                 write(iu,414,err=1999) temp
              else
                 write(iu,rec=irec,err=1999) temp
              endif
           end do
        endif
c
c        write image data
c
        iucnt=0
        iocnt=0
        ntot = nsx*nsy
        nrec=256                      ! output values/record (I*2)
        if (idatatype.eq.1) nrec=512  ! output values/record (byte)
        if (idatatype.eq.4) nrec=128  ! output values/record (float)
        do j = 0,ntot,nrec
           imax=min(nrec,ntot-j)
           if (imax.gt.0) then
              do i=1,imax
                 if (idatatype.eq.4) then
                    tempf(i)=stval(j+i)
                 else
                    k = nint(iscale*(stval(j+i)-ioff))-minv
                    if (k.gt.maxv) then
                       k=maxv
                       iocnt=iocnt+1
                    endif
                    if (k.lt.-minv) then
                       k=-minv
                       iucnt=iucnt+1
                    endif
                    if (idatatype.eq.1) then
                       tempb(i) = k
                    else
                       temp(i) = k
                    endif
                 endif
              end do
              if (imax.ne.nrec) then
                 do i=imax+1,nrec
                    if (idatatype.eq.1) then
                       tempb(i) = -minv
                    else if (idatatype.eq.4) then
                       tempf(i) = anodata
                    else
                       temp(i) = -minv
                    endif
                 end do
              endif
              irec=irec+1
              if (idatatype.eq.1) then ! byte size
                 if (VMS) then
                    write(iu,413,err=1999) tempb
 413                    format(512a1)
                 else
                    write(iu,rec=irec,err=1999) tempb
                 endif
              else if (idatatype.eq.4) then ! 4 byte float size
                 if (VMS) then
                    write(iu,412,err=1999) tempf
 412                    format(128a4)
                 else
                    write(iu,rec=irec,err=1999) tempf
                 endif
              else                ! integer*2 size
                 if (BSWAP) call swapbuf(temp,256)
                 if (VMS) then
                    write(iu,414,err=1999) temp
 414                    format(256a2)
                 else
                    write(iu,rec=irec,err=1999) temp
                 endif
              endif
           endif
        end do 
c
c        close file
c
        close(unit=iu)
c
c     summarize over/underflow stats
c
        if (iucnt.gt.0) write (*,50) iucnt,ntot
50        format(' *** WRITESIR3 Underflow count: ',I9,' of ',I9,' ***')
        if (iocnt.gt.0) write (*,60) iocnt,ntot
60        format(' *** WRITESIR3 Overflow count:  ',I9,' of ',I9,' ***')
        return
c
 999    continue
        write (*,5000)
 5000   format(' *** writesir3: ERROR OPENING OUTPUT FILE ***')
        ierr=-1
        return
 1999   continue
        write (*,5050)
 5050   format(' *** writesir3: ERROR WRITING OUTPUT FILE ***')
        ierr=-2
        return
        end
c
c *************************************************************************
c
c  header version 2.0 interface to version 3.0 standard sir write routine
c
      subroutine writesir(fname,iu_in,ierr,nhead,ndes,nhtype,idatatype,
     *     nsx,nsy,xdeg,ydeg,ascale,bscale,a0,b0,iopt,ioff,iscale,
     *     iyear,isday,ismin,ieday,iemin,iregion,itype,
     *     ipol,ifreqhm,ispare1,anodata,vmin,vmax,
     *     sensor,title,type,tag,crproc,crtime,
     *     stval,descrip,ldes,iaopt,nia)
c
        implicit none
c
        character*(*)         fname
        integer                iu_in,ierr,nhead,ndes,nhtype,idatatype
        integer         nsx,nsy,iopt
        real                xdeg,ydeg,ascale,bscale,a0,b0
        integer                ioff,iscale,iyear,isday,ismin,ieday,iemin
        integer                iregion,itype,ipol,ifreqhm,ispare1,ldes,nia
        real            anodata,vmin,vmax
        real            stval(*)
        character*(*)        sensor,title,type,tag,crproc,crtime,descrip
        integer                iaopt(*)
c
        integer         ixdeg_off,iydeg_off,ideg_sc,iscale_sc
        integer         ia0_off,ib0_off,i0_sc
c
c     set version 3.0 parameters to header version 2.0 defaults
c
        if (iopt.eq.-1) then    ! image only
           ideg_sc=10
           iscale_sc=1000
           i0_sc=100
           ixdeg_off=0
           iydeg_off=0
           ia0_off=0
           ib0_off=0
        else if (iopt.eq.0) then ! rectalinear lat/lon
           ideg_sc=100
           iscale_sc=1000
           i0_sc=100
           ixdeg_off=-100
           iydeg_off=0
           ia0_off=0
           ib0_off=0
        else if (iopt.eq.1.or.iopt.eq.2) then ! lambert
           ideg_sc=100
           iscale_sc=1000
           i0_sc=1
           ixdeg_off=0
           iydeg_off=0
           ia0_off=0
           ib0_off=0
        else if (iopt.eq.5) then ! polar stereographic
           ideg_sc=100
           iscale_sc=100
           i0_sc=1
           ixdeg_off=-100
           iydeg_off=0
           ia0_off=0
           ib0_off=0
        else if (iopt.eq.11.or.iopt.eq.12.or.iopt.eq.13) then ! EASE grid
           ideg_sc=10
           iscale_sc=1000
           i0_sc=10
           ixdeg_off=0
           iydeg_off=0
           ia0_off=0
           ib0_off=0
        else                    ! unknown default scaling
           ideg_sc=100
           iscale_sc=1000
           i0_sc=100
           ixdeg_off=0
           iydeg_off=0
           ia0_off=0
           ib0_off=0
        endif
c
c     call version 3.0 write sir routine
c
        call writesir3(fname,iu_in,ierr,nhead,ndes,nhtype,idatatype,
     *     nsx,nsy,xdeg,ydeg,ascale,bscale,a0,b0,iopt,ioff,iscale,
     *     ixdeg_off,iydeg_off,ideg_sc,iscale_sc,ia0_off,ib0_off,i0_sc,
     *     iyear,isday,ismin,ieday,iemin,iregion,itype,
     *     ipol,ifreqhm,ispare1,anodata,vmin,vmax,
     *     sensor,title,type,tag,crproc,crtime,
     *     stval,descrip,ldes,iaopt,nia)
c
        return
        end
