! f90 easy BYU SIR file format interface code
!
! These routines define an easy-to-use fortran 90 programing interface
! for the BYU SIR file format.  The file SIREZ_f90.inc defines a
! custom data type for holding SIR file header information.  The routines
! call the standard (f77) routines for reading, writing, and geometric
! computation on the SIR files.
!
! written DGL 30 Jan. 2002
! (c) 2000 by BYU MERS 

! ***************************************************************

SUBROUTINE getsirhead(fname,lu,head,ierr)

!     read sir header information from file named fname into
!     data structure head using file unit lu.  Leave file open.

!     RETURNed error code ierr:
!     ierr   file read error code
!             set to 0 for successful read of header
!             set to -1 for open error
!             set to -2 for header read error

  INCLUDE 'SIREZ_f90.inc'

  CHARACTER(len=*), INTENT(IN) :: fname
  INTEGER, INTENT(inout) :: lu
  TYPE(sirhead), INTENT(inout) :: head
  INTEGER, INTENT(inout) :: ierr

  CALL readsirhead3(fname,lu,ierr, head%nhead, head%ndes, head%nhtype, &
       head%idatatype, head%nsx, head%nsy, head%xdeg, head%ydeg, &
       head%ascale, head%bscale, head%a0, head%b0, head%iopt, &
       head%ioff, head%iscale, &
       head%ixdeg_off, head%iydeg_off, head%ideg_sc, head%iscale_sc, &
       head%ia0_off, head%ib0_off, head%i0_sc, &
       head%iyear, head%isday, head%ismin, head%ieday, head%iemin, &
       head%iregion, head%itype, head%ipol, head%ifreqhm, head%ispare1, &
       head%anodata, head%vmin, head%vmax, &
       head%sensor, head%title, head%type, head%tag, head%crproc, &
       head%crtime, head%MAXDES, head%descrip, head%ldes, &
       head%MAXI, head%iaopt, head%nia)

  RETURN
END SUBROUTINE getsirhead


SUBROUTINE getsirdata(lu,head,stval,ierr,smin,smax)
!
!     read sir image data information from file attached to 
!     file unit lu into real array stval using sir header
!     data structure head%  Closes file on successful completion.
!     smin and smax RETURN the min,max of the data greater than the
!     nodata value.
!
!     RETURNed error code ierr:
!     ierr   file read error code
!             set to 0 for successful read of file
!             set to -3 for data read error
!

  INCLUDE 'SIREZ_f90.inc'

  INTEGER, INTENT(inout) :: lu
  TYPE(sirhead), INTENT(inout) :: head
  INTEGER, INTENT(inout) :: ierr
  REAL, INTENT(inout) :: smin, smax

  INTEGER :: ncnt

  CALL readsirf(lu, ierr, head%nhead, head%nhtype, head%idatatype, &
       stval, head%nsx, head%nsy, head%ioff, head%iscale, &
       smin, smax, ncnt, head%anodata, head%vmin, head%vmax)

  RETURN
END SUBROUTINE getsirdata


SUBROUTINE printsirhead(head)
!
!     print out summary contents of sir file header information
!     contained in sir header data structure head%  Output is to
!     standard out
!

  INCLUDE 'SIREZ_f90.inc'

  TYPE(sirhead), INTENT(in) :: head

  CALL printhead3(head%nhead, head%ndes, head%nhtype, head%idatatype, &
       head%nsx, head%nsy, head%xdeg, head%ydeg, &
       head%ascale, head%bscale, head%a0, head%b0, head%iopt, &
       head%ioff, head%iscale, &
       head%ixdeg_off, head%iydeg_off, head%ideg_sc, head%iscale_sc, &
       head%ia0_off, head%ib0_off, head%i0_sc, &
       head%iyear, head%isday, head%ismin, head%ieday, head%iemin, &
       head%iregion, head%itype, head%ipol, head%ifreqhm, head%ispare1, &
       head%anodata, head%vmin, head%vmax, &
       head%SENSOR, head%TITLE, head%TYPE, head%TAG, head%CRPROC, &
       head%CRTIME, head%DESCRIP, head%LDES, head%IAOPT, head%NIA)

  RETURN
END SUBROUTINE printsirhead


SUBROUTINE putsirfile(fname,lu,head,stval,ierr)
!
!     write sir image data and header to file
!     data structure head using file unit lu
!
!     RETURNed error code ierr:
!     ierr   file write error code
!             set to 0 for successful write
!             set to -1 for open error
!             set to -2 for write error
!             set to -3 for invalid image size

  INCLUDE 'SIREZ_f90.inc'

  CHARACTER(len=*), INTENT(IN) :: fname
  INTEGER, INTENT(inout) :: lu
  TYPE(sirhead), INTENT(inout) :: head
  INTEGER, INTENT(inout) :: ierr
  REAL, DIMENSION(*), INTENT(in) :: stval

  CALL writesir3(fname, lu, ierr, head%nhead, head%ndes, head%nhtype, &
       head%idatatype, head%nsx, head%nsy, head%xdeg, head%ydeg, &
       head%ascale, head%bscale, head%a0, head%b0, head%iopt, &
       head%ioff, head%iscale, &
       head%ixdeg_off, head%iydeg_off, head%ideg_sc, head%iscale_sc, &
       head%ia0_off, head%ib0_off, head%i0_sc, &
       head%iyear, head%isday, head%ismin, head%ieday, head%iemin, &
       head%iregion, head%itype, head%ipol, head%ifreqhm, head%ispare1, &
       head%anodata, head%vmin, head%vmax, &
       head%sensor, head%title, head%type, head%tag, head%crproc, &
       head%crtime,stval, head%descrip, head%ldes, &
       head%iaopt, head%nia)

  RETURN
END SUBROUTINE putsirfile


SUBROUTINE sirpix2latlon(x, y, alon, alat, head)
!
!     compute lat,lon corresponding to pixel x,y given info in sir
!     header data structure head%  x,y can be outside of the image
!     but this may cause an invalid lat,lon to be RETURNed or an algorithm
!     error
!
  REAL, INTENT(in)  :: x, y
  REAL, INTENT(out) :: alon, alat

  INCLUDE 'SIREZ_f90.inc'
  TYPE(sirhead), INTENT(in) :: head

  REAL :: thelon, thelat

  CALL pix2latlon(x, y, thelon, thelat, alon, alat, &
       head%iopt, head%xdeg, head%ydeg, head%ascale, head%bscale, &
       head%a0, head%b0)

  RETURN
END SUBROUTINE sirpix2latlon


INTEGER function isirlatlon2pix(alon, alat, x, y, head)
!
!     compute pixel coordinates corresponding to lat,lon given info in sir
!     header data structure head%  RETURNs lexicographic index to
!     corresponding pixel (only valid if > 0).  The image coordinate
!     range is defined as 1 <= x <= head%nsx and 0 <= y <= head%nsy
!     with the pixel (1,1) at index 1 and (nsx,nsy) at index nsx*nsy
!     if x or y out of image, an invalid index of 0 is RETURNed.
!
!
  REAL, INTENT(out) :: x, y
  REAL, INTENT(in)  :: alon, alat
  
  INCLUDE 'SIREZ_f90.inc'
  TYPE(sirhead), INTENT(in) :: head
  
  REAL :: thelon, thelat
  INTEGER :: ix, iy

  CALL latlon2pix(alon, alat, x, y, thelon, thelat, &
       head%iopt, head%xdeg, head%ydeg, head%ascale, head%bscale, &
       head%a0, head%b0)
  CALL f2ipix(x, y, ix, iy, head%nsx, head%nsy)
  IF (ix.NE.0.AND.iy.NE.0) THEN
     isirlatlon2pix=(iy-1) * head%nsx + ix
  ELSE
     isirlatlon2pix=0
  ENDIF
  RETURN
END FUNCTION isirlatlon2pix

INTEGER function isirpix(x, y, ix, iy, head)
!
!     compute pixel index from the coordinates given info in sir
!     header data structure head%  RETURNs lexicographic index to
!     corresponding pixel (only valid if > 0)  The image coordinate
!     range is defined as 1 <= x <= head%nsx and 0 <= y <= head%nsy
!     with the pixel (1,1) at index 1 and (nsx,nsy) at index nsx*nsy
!     if x or y out of image, an invalid index of 0 is RETURNed and
!     one of ix and iy are set to zero.
!

  REAL, INTENT(in)  :: x, y
  INTEGER, INTENT(out) :: ix, iy
  
  INCLUDE 'SIREZ_f90.inc'
  TYPE(sirhead), INTENT(in) :: head

  CALL f2ipix(x, y, ix, iy, head%nsx, head%nsy)
  IF (ix.NE.0.AND.iy.NE.0) THEN
     isirpix = (iy-1) * head%nsx + ix
  ELSE
     isirpix=0
  ENDIF
  RETURN
END FUNCTION isirpix


INTEGER function isirlex(ix, iy, head)
!
!     compute pixel coordinates corresponding to lat,lon given info in sir
!     header data structure head%  RETURNs lexicographic index to
!     corresponding pixel (only valid if > 0).  The image coordinate
!     range is defined as 1 <= x <= head%nsx and 0 <= y <= head%nsy
!     with the pixel (1,1) at index 1 and (nsx,nsy) at index nsx*nsy
!     if ix or iy out of image, an invalid index of 0 is RETURNed

  INTEGER, INTENT(in) :: ix, iy
  
  INCLUDE 'SIREZ_f90.inc'
  TYPE(sirhead), INTENT(in) :: head

  IF (ix.GT.0.AND.ix.LE.head%nsx.AND. &
      iy.GT.0.AND.iy.LE.head%nsy) THEN
     isirlex=(iy-1)*head%nsx+ix
  ELSE
     isirlex=0
  ENDIF
  RETURN
END FUNCTION isirlex