program use_of_span
use, intrinsic :: iso_fortran_env
implicit none

! Read a boundary field from HIRLAM to determine the grid.
! Then read the list of the result data from the surface data assimilation.
! Then, determine if they fall within the grid
! and write the integer lat/lon index pair for its location,
! followed by the flag for its use (1 or 4).

character*256 climate, spanfile, filename
integer :: iunit, istatus, ihour, inarea, irec, iused, irejected
logical lopened(10:113) ! Indicate which output files are opened
real obslat, obslon

! Contents of the surface data assimiliation result data file

character*8  :: statid1, statid2, name(10)
real         :: pobs(18+4*16)
real         :: oblat, oblon, obheight
integer      :: jobs, nobs
integer      :: idate, itime
integer      :: iobtype, iflag1, iflag

! Declarations for reading ASIMOF files

integer   jplb1, jplb2, jpsup, jplf
parameter (jplb1 = 30, jplb2 = 400, jpsup = 40, jplf = 1600000)

integer   ilf
integer   ibs1(jplb1)
integer   ibs2(jplb2)
integer   ibr2(jplb2)
integer   ibpw, ibpwio, imdi
integer   ierr
integer   ilon, ilat
integer   nlon, nlat
real      west, south
real      zdir
real      dlon, dlat
real      plon, plat
real      zbs2(jplb2)
real      zbr2(jplb2)
real      zsup(jpsup)
real      landseamask(jplf)

data      ibpw /0/, ibpwio /0/

! Get the command line arguments

call get_command_argument(1, climate)
call get_command_argument(2, spanfile)

! Get landseamask and grid definition of HIRLAM

call ASIMHM(ibpw, ibpwio, imdi)

call ASIMHO(1, climate, 0, ierr)
if (ierr /= 0) then
   write(0,*) 'USE-SPAN: CANNOT OPEN ',trim(climate),' ASIMHO RETURNS:', ierr
   call abort
endif

ibs1 = imdi
ibs1(9:11) = [81, 105, 0] ! Parameter number, level type, level number of land/sea mask

ilf = jplf
call GETFD(1, ibs1, jplb1, ibs2, jplb2, zbs2, jplb2, zsup, jpsup, landseamask, ilf, 0, ierr)
if (ierr < 0) then
   write(0,*) 'USE-SPAN: ERROR IN GETFD; IERR=',ierr
   write(0,'(30(1x,i3))') ibs1
   call abort
endif

call ASIMHC(1, ierr)

nlon = ibs2( 7)
nlat = ibs2( 9)
south= ibs2(11) / 1000.0
west = ibs2(14) / 1000.0
dlon = ibs2(24) / 1000.0
dlat = ibs2(26) / 1000.0
plat = ibs2(33) / 1000.0
plon = ibs2(36) / 1000.0

! Decode and write out

ibr2(6)   = 0       ! Defines regular lat/lon grid.
lopened   = .false.
irec      = 0
inarea    = 0
iused     = 0
irejected = 0

name(1) = 'SPAN'
name(9) = 'PSEUDO'

open(unit = 7, file = spanfile)

do

   READ(unit = 7,fmt = '(8f12.2,1x,a8,1x,a8/2X,8f10.4/(16f10.4))',iostat = istatus) &
     (pobs(jobs),jobs=1,8),statid1,statid2,    &
     (pobs(jobs),jobs=11,nint(pobs(1)))
   if (istatus == iostat_end) exit
   if (istatus /= 0) cycle

   iobtype = pobs(4)
   if (iobtype /= 1 .and. iobtype /= 9) then
      write(0,*) 'USE-SPAN: UNRECOGNISED OBSERVATION TYPE:', iobtype
      cycle
   endif
   ihour = nint(pobs(7) / 100.)
   oblat = pobs(11)
   oblon = pobs(12)
   nobs = (nint(pobs(1)) - 18) / 16
   iflag1 = 0
   do jobs = 1, nobs 
      irec = irec + 1

      iflag = pobs(18 + (jobs-1) * 16 + 3)
      if (iflag1 == 0) iflag1 = iflag       ! Keep only 1 flags,
      if (iflag  == 1) iflag1 = iflag       ! if any.

      call gb2gca(oblon, oblat, obslon, obslat, zdir, 1, ibr2, zbr2, ibs2, zbs2)

      ilat = int((obslat - south) / dlat + 1.0)
      if (ilat < 1 .or. ilat > nlat) cycle               ! Not in area - cop out.

      ilon = int((obslon - west ) / dlon + 1.0)
      if (ilon < 1 .or. ilon > nlon) cycle               ! Not in area - cop out.

      inarea = inarea + 1
      if (iflag1 == 1) iused = iused + 1
      if (iflag1 == 4) irejected = irejected + 1
      iunit  = 10 * iobtype + ihour
      if (.not. lopened(iunit)) then
         write(filename, '(a,i2.2,a)') trim(name(iobtype)), ihour, '.txt'
         open(iunit, file = filename)
         lopened(iunit) = .true.
      endif
      write(iunit, *) ilon, ilat, iflag1
   enddo

enddo

! Close input file

close(unit = 7)

! Close output files

do iunit = lbound(lopened, dim=1), ubound(lopened, dim=1)
   if (lopened(iunit)) close(iunit)
enddo

open(unit=8,file='SPAN_COUNT.txt'); write(8,*) iused, irejected; close(unit=8)
write(0,*)'USE-SPAN: NUMBER OF OBSERVATIONS:',irec, &
   ', INSIDE AREA:',inarea,', USED:',iused,', REJECTED:',irejected

! The End

end program use_of_span