program use_of_synop
use, intrinsic :: iso_fortran_env
implicit none

! Read a boundary field from HIRLAM to determine the grid.
! Then read the list of synop data processed by the 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, synopfile, filename
integer :: iunit, istatus, ihour, inarea, irec, iused, iredundant, irejected
logical lopened(10:33) ! Indicate which output files are opened
real obslat, obslon

! Contents of the synop use data file

character*10 :: idname
integer      :: idint1, idint2
real         :: oblat, oblon, obheight
integer      :: idate, itime
integer      :: iflag1, iflag2

! 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, synopfile)

! 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-SYNOP: 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-SYNOP: 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
iredundant= 0

open(unit = 7, file = synopfile)

do

   read(unit = 7, fmt = *, iostat = istatus) &
      idname, idint1, idint2,                &
      oblat, oblon, obheight,                &
      idate, itime,                          &
      iflag1, iflag2
   if (istatus == iostat_end) exit
   if (istatus /= 0) cycle

   irec = irec + 1

   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.

   ihour = itime / 10000
   ihour = ihour + findloc([.false.,.true.],mod(itime,10000)>=3000,dim=1) - 1
   ihour = mod(ihour, 24)

   inarea = inarea + 1
   if (iflag1 == 1) iused = iused + 1
   if (iflag1 == 4) then
      if (iflag2 == 2) then
         iflag1    = 3
         iredundant= iredundant + 1
      else
         irejected = irejected + 1
      endif
   endif
   iunit  = 10 + ihour
   if (.not. lopened(iunit)) then
      write(filename, '(a,i2.2,a)') 'SYNOP', ihour, '.txt'
      open(iunit, file = filename)
      lopened(iunit) = .true.
   endif
   write(iunit, *) ilon, ilat, iflag1

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='SYNOP_COUNT.txt'); write(8,*) iused, iredundant, irejected; close(unit=8)
write(0,*)'USE-SYNOP: NUMBER OF OBSERVATIONS:',irec, &
   ', INSIDE AREA:',inarea,', USED:',iused,', REDUNDANT:',iredundant,', REJECTED:',irejected

! The End

end program use_of_synop