PROGRAM IRDUMPC
! This program reads the files irtemp1.dat and irtemp2.dat in McIDAS AREA
! format and dumps them in an ASCII format to files irdump1.dat irdump2.dat.
! These files have a header containing date/time and array size info, followed
! by three blocks of data (IR brightness temp, lon and lat). These are written
! in rows of 20 f9.3 values. The lon is in 0 to 360 deg E convention and lat is deg N.
! The lon and lat are written for every BT point, even when the lon/lat grid is 
! rectilinear or evenly spaced, for generality. 
!
! This version also reads an input file (irdumpc.inp) containing a latitude and
! longitude (0 to 360 deg convention) and "crops" the image by only printing
! out values near that point.
!
! Ver 2.0.0
!
! Modified 5/2/2016 to process three IR files instead of two for 2016 
!                   SHIPS/LGEM/RII implementation. (MD)
  USE merc_4km
  IMPLICIT NONE

  INTEGER:: i,j,k
  INTEGER:: i1,i2,j1,j2
  INTEGER:: icen,jcen,ibuf,jbuf
  INTEGER:: incxy,ludat,luinp
  INTEGER:: nx,ny

  REAL:: slon,slat,dll,dllt

  CHARACTER (LEN=80) :: FN

  ludat = 51
  luinp = 52

! Specify the number of pixels to include surrounding the
! center location
  ibuf = 200
  jbuf = 200

! Open and read the input file
  open(unit=luinp,file='irdumpc.inp',form='formatted',status='old')
  read(luinp,*) slat,slon
  close(luinp)

  FN='irtemp1.dat'

  IF (irar_get(FN)) THEN

!    Find the line and element closest to slat,slon
     dll = 10000.0
     icen = 0
     do i=1,nelex
        dllt = abs(slon-lon(i))
        if (dllt .lt. dll) then
           icen = i
           dll  = dllt
        endif
     enddo
     i1 = icen - ibuf
     i2 = icen + ibuf
     if (i1 .lt.     1) i1 = 1
     if (i2 .gt. nelex) i2 = nelex
 
     dll = 10000.0
     jcen = 0
     do j=1,nliny
        dllt = abs(slat-lat(j))
        if (dllt .lt. dll) then
           jcen = j
           dll  = dllt
        endif
     enddo
     j1 = jcen - jbuf
     j2 = jcen + jbuf
     if (j1 .lt.     1) j1 = 1
     if (j2 .gt. nliny) j2 = nliny
 
     IF (icen .gt. 0 .and. jcen .gt. 0) THEN
        open(unit=ludat,file='irdump1.dat',form='formatted',status='replace')
 
        nx = 1 + (i2-i1)
        ny = 1 + (j2-j1)

        write(ludat,100) dtgg,time,nx,ny
 100    format(i8,1x,i6.6,1x,i6,1x,i6)
     
        write(ludat,110) ((bt(i,j),i=i1,i2),j=j1,j2)
 110    format(20(f9.3,1x))

        write(ludat,110) ((lon(i),i=i1,i2),j=j1,j2)
        write(ludat,110) ((lat(j),i=i1,i2),j=j1,j2)

        close(ludat)
     END IF
        
  END IF

  FN='irtemp2.dat'

  IF (irar_get(FN)) THEN
!    Find the line and element closest to slat,slon
     dll = 10000.0
     icen = 0
     do i=1,nelex
        dllt = abs(slon-lon(i))
        if (dllt .lt. dll) then
           icen = i
           dll  = dllt
        endif
     enddo
     i1 = icen - ibuf
     i2 = icen + ibuf
     if (i1 .lt.     1) i1 = 1
     if (i2 .gt. nelex) i2 = nelex
 
     dll = 10000.0
     jcen = 0

     do j=1,nliny
        dllt = abs(slat-lat(j))
        if (dllt .lt. dll) then
           jcen = j
           dll  = dllt
        endif
     enddo
     j1 = jcen - jbuf
     j2 = jcen + jbuf
     if (j1 .lt.     1) j1 = 1
     if (j2 .gt. nliny) j2 = nliny
 
     IF (icen .gt. 0 .and. jcen .gt. 0) THEN
        open(unit=ludat,file='irdump2.dat',form='formatted',status='replace')

        nx = 1 + (i2-i1)
        ny = 1 + (j2-j1)

        write(ludat,100) dtgg,time,nx,ny
     
        write(ludat,110) ((bt(i,j),i=i1,i2),j=j1,j2)
        write(ludat,110) ((lon(i),i=i1,i2),j=j1,j2)
        write(ludat,110) ((lat(j),i=i1,i2),j=j1,j2)

        close(ludat)
     END IF

  END IF

  FN='irtemp3.dat'

  IF (irar_get(FN)) THEN
!    Find the line and element closest to slat,slon
     dll = 10000.0
     icen = 0
     do i=1,nelex
        dllt = abs(slon-lon(i))
        if (dllt .lt. dll) then
           icen = i
           dll  = dllt
        endif
     enddo
     i1 = icen - ibuf
     i2 = icen + ibuf
     if (i1 .lt.     1) i1 = 1
     if (i2 .gt. nelex) i2 = nelex
 
     dll = 10000.0
     jcen = 0

     do j=1,nliny
        dllt = abs(slat-lat(j))
        if (dllt .lt. dll) then
           jcen = j
           dll  = dllt
        endif
     enddo
     j1 = jcen - jbuf
     j2 = jcen + jbuf
     if (j1 .lt.     1) j1 = 1
     if (j2 .gt. nliny) j2 = nliny
 
     IF (icen .gt. 0 .and. jcen .gt. 0) THEN
        open(unit=ludat,file='irdump3.dat',form='formatted',status='replace')

        nx = 1 + (i2-i1)
        ny = 1 + (j2-j1)

        write(ludat,100) dtgg,time,nx,ny
     
        write(ludat,110) ((bt(i,j),i=i1,i2),j=j1,j2)
        write(ludat,110) ((lon(i),i=i1,i2),j=j1,j2)
        write(ludat,110) ((lat(j),i=i1,i2),j=j1,j2)

        close(ludat)
     END IF

  END IF

  STOP
END PROGRAM IRDUMPC