program decode_eccodes

! Read a climate field from HIRLAM to determine the grid
! Then read all BUFR messages from the observation file,
! decode them, determine if they fall within the grid,
! and write the integer lat/lon index pair for its location,
! followed by the height.

use eccodes
implicit none

integer, parameter                                      :: max_strsize = 200
integer                                                 :: iret
integer                                                 :: ifile
integer                                                 :: ibufr
integer(kind=4)                                         :: iVal, iStat
real(kind=8)                                            :: rVal
character(len=max_strsize)                              :: sVal
integer(kind=4), dimension(:), allocatable              :: iValues
character(len=max_strsize) , dimension(:),allocatable   :: sValues
real(kind=8), dimension(:), allocatable                 :: rValues

character(len=max_strsize)                              :: infile_name

character*8 :: ident
real :: winddir, windspd, temp, dewp
real :: rh, vis, cc, rr, gust, snow, sst
integer :: iperiod

logical lopened(100:600)    ! Indicate which output files are opened
character*11  filetype(5) /'SYNOP','SHIP','TEMP','','AIREP-AMDAR'/

character*256 climate, obsfile, filename

! 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      orography(jplf)

integer   iobtype, iobsubtype, irec, nelem, ilen, ihour, iminute, iunit, inarea, nsynop, nship, ntemp, nairep
real      oblat, obslat, oblon, obslon, obheight

data      ibpw /0/, ibpwio /0/

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

! Get orography and grid definition of HIRLAM

call ASIMHM(ibpw, ibpwio, imdi)

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

ibs1 = imdi
ibs1(9:11) = [6, 105, 0] ! Parameter number, level, height of orography

ilf = jplf
call GETFD(1, ibs1, jplb1, ibs2, jplb2, zbs2, jplb2, zsup, jpsup, orography, ilf, 0, ierr)
if (ierr < 0) then
   write(0,*) 'DECODE-ECCODES: 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

lopened = .false. ! Initially, none of the output files is opened.

ibr2(6) = 0       ! GRIB Section 2 definition for "regular lat/lon grid"

nsynop = 0; nship = 0; ntemp = 0; nairep = 0
inarea = 0
irec   = 0

call codes_open_file(ifile, obsfile, 'r')
open(7, file='spanobs.txt')
do

   irec = irec + 1

   call codes_bufr_new_from_file(ifile, ibufr, iStat)
   if (iStat == CODES_END_OF_FILE) exit
   call codes_set(ibufr, 'unpack', 1)
   if(allocated(iValues)) deallocate(iValues)
   call codes_get(ibufr, 'delayedDescriptorReplicationFactor', iValues, iStat)
   call codes_get(ibufr, 'edition', iVal, iStat)
   call codes_get(ibufr, 'masterTableNumber', iVal, iStat)
   call codes_get(ibufr, 'bufrHeaderSubCentre', iVal, iStat)
   call codes_get(ibufr, 'bufrHeaderCentre', iVal, iStat)
   call codes_get(ibufr, 'updateSequenceNumber', iVal, iStat)
   call codes_get(ibufr, 'dataCategory', iVal, iStat)
   if (iStat == CODES_SUCCESS) then
      iobtype = iVal
   endif
   call codes_get(ibufr, 'dataSubCategory', iVal, iStat)
   if (iStat == CODES_SUCCESS) then
      iobsubtype = iVal
   endif
   call codes_get(ibufr, 'masterTablesVersionNumber', iVal, iStat)
   call codes_get(ibufr, 'localTablesVersionNumber', iVal, iStat)
   call codes_get(ibufr, 'typicalYearOfCentury', iVal, iStat)
   call codes_get(ibufr, 'typicalMonth', iVal, iStat)
   call codes_get(ibufr, 'typicalDay', iVal, iStat)
   call codes_get(ibufr, 'typicalHour', iVal, iStat)
   if (iStat == CODES_SUCCESS) then
      ihour   = iVal
   endif
   call codes_get(ibufr, 'typicalMinute', iVal, iStat)
   if (iStat == CODES_SUCCESS) then
      iminute = iVal
   endif
   call codes_get(ibufr, 'rdbType', iVal, iStat)
   call codes_get(ibufr, 'oldSubtype', iVal, iStat)
   call codes_get(ibufr, 'localYear', iVal, iStat)
   call codes_get(ibufr, 'localMonth', iVal, iStat)
   call codes_get(ibufr, 'localDay', iVal, iStat)
   call codes_get(ibufr, 'localHour', iVal, iStat)
   call codes_get(ibufr, 'localMinute', iVal, iStat)
   call codes_get(ibufr, 'localSecond', iVal, iStat)
   call codes_get(ibufr, 'rdbtimeDay', iVal, iStat)
   call codes_get(ibufr, 'rdbtimeHour', iVal, iStat)
   call codes_get(ibufr, 'rdbtimeMinute', iVal, iStat)
   call codes_get(ibufr, 'rdbtimeSecond', iVal, iStat)
   call codes_get(ibufr, 'rectimeDay', iVal, iStat)
   call codes_get(ibufr, 'rectimeHour', iVal, iStat)
   call codes_get(ibufr, 'rectimeMinute', iVal, iStat)
   call codes_get(ibufr, 'rectimeSecond', iVal, iStat)
   call codes_get(ibufr, 'restricted', iVal, iStat)
   call codes_get(ibufr, 'correction1', iVal, iStat)
   call codes_get(ibufr, 'correction1Part', iVal, iStat)
   call codes_get(ibufr, 'correction2', iVal, iStat)
   call codes_get(ibufr, 'correction2Part', iVal, iStat)
   call codes_get(ibufr, 'correction3', iVal, iStat)
   call codes_get(ibufr, 'correction3Part', iVal, iStat)
   call codes_get(ibufr, 'correction4', iVal, iStat)
   call codes_get(ibufr, 'correction4Part', iVal, iStat)
   call codes_get(ibufr, 'qualityControl', iVal, iStat)
   call codes_get(ibufr, 'newSubtype', iVal, iStat)
   call codes_get(ibufr, 'numberOfSubsets', iVal, iStat)
   call codes_get(ibufr,'localLatitude', rVal, iStat)
   if (iStat == CODES_SUCCESS) then
      oblat    = rVal
   endif
   call codes_get(ibufr,'localLongitude', rVal, iStat)
   if (iStat == CODES_SUCCESS) then
      oblon    = rVal
   endif
   call codes_get(ibufr, 'ident', sVal, iStat)
   if (iStat == CODES_SUCCESS) then
      ident    = sVal
   endif
   call codes_get(ibufr, 'observedData', iVal, iStat)
   call codes_get(ibufr, 'compressedData', iVal, iStat)
   if(allocated(iValues)) deallocate(iValues)
   call codes_get(ibufr, 'unexpandedDescriptors', iValues)
   call codes_get(ibufr, 'blockNumber', iVal, iStat)
   call codes_get(ibufr, 'stationNumber', iVal, iStat)
   call codes_get(ibufr, 'stationOrSiteName', sVal, iStat)
   call codes_get(ibufr, 'stationType', iVal, iStat)
   call codes_get(ibufr, 'year', iVal, iStat)
   call codes_get(ibufr, 'month', iVal, iStat)
   call codes_get(ibufr, 'day', iVal, iStat)
   call codes_get(ibufr, 'hour', iVal, iStat)
   call codes_get(ibufr, 'minute', iVal, iStat)
   call codes_get(ibufr, '#1#latitude', rVal, iStat)
   call codes_get(ibufr, '#1#longitude', rVal, iStat)
   call codes_get(ibufr, 'heightOfStationGroundAboveMeanSeaLevel', rVal, iStat)
   if (iobtype == 0) then ! Additional data to obtain from SYNOP observations
      winddir = 1.7e38
      windspd = 1.7e38
      temp = 1.7e38
      dewp = 1.7e38
      rh = 1.7e38
      vis = 1.7e38
      cc = 1.7e38
      rr = 1.7e38
      gust = 1.7e38
      snow = 1.7e38
      call codes_get(ibufr, 'heightOfBarometerAboveMeanSeaLevel', rVal, iStat)
      call codes_get(ibufr, 'nonCoordinatePressure', rVal, iStat)
      call codes_get(ibufr, 'pressureReducedToMeanSeaLevel', rVal, iStat)
      call codes_get(ibufr, '3HourPressureChange', rVal, iStat)
      call codes_get(ibufr, 'characteristicOfPressureTendency', iVal, iStat)
      call codes_get(ibufr, '#1#heightOfSensorAboveLocalGroundOrDeckOfMarinePlatform', rVal, iStat)
      call codes_get(ibufr, 'airTemperature', rVal, iStat)
      if (iStat == CODES_SUCCESS .and. rVal > 0 .and. rVal < 1.e6) then
         temp    = rVal
      endif
      call codes_get(ibufr, 'dewpointTemperature', rVal, iStat)
      if (iStat == CODES_SUCCESS .and. rVal > 0 .and. rVal < 1.e6) then
         dewp    = rVal
      endif
      call codes_get(ibufr, 'relativeHumidity', iVal, iStat)
      if (iStat == CODES_SUCCESS .and. iVal >= 0 .and. iVal < 1000000) then
         rh      = iVal
      endif
      call codes_get(ibufr, '#2#heightOfSensorAboveLocalGroundOrDeckOfMarinePlatform', rVal, iStat)
      call codes_get(ibufr, 'horizontalVisibility', rVal, iStat)
      if (iStat == CODES_SUCCESS .and. rVal >= 0 .and. rVal < 1.e6) then
         vis     = rVal
      endif
      call codes_get(ibufr, '#3#heightOfSensorAboveLocalGroundOrDeckOfMarinePlatform', rVal, iStat)
      call codes_get(ibufr, 'cloudCoverTotal', iVal, iStat)
      if (iStat == CODES_SUCCESS .and. iVal >= 0 .and. iVal < 1000000) then
         cc      = iVal
      endif
      call codes_get(ibufr, '#1#cloudAmount', iVal, iStat)
      call codes_get(ibufr, '#1#heightOfBaseOfCloud', rVal, iStat)
      call codes_get(ibufr, '#2#cloudAmount', iVal, iStat)
      call codes_get(ibufr, '#2#heightOfBaseOfCloud', rVal, iStat)
      call codes_get(ibufr, '#3#cloudAmount', iVal, iStat)
      call codes_get(ibufr, '#3#heightOfBaseOfCloud', rVal, iStat)
      call codes_get(ibufr, '#4#cloudAmount', iVal, iStat)
      call codes_get(ibufr, '#4#heightOfBaseOfCloud', rVal, iStat)
      call codes_get(ibufr, '#6#verticalSignificanceSurfaceObservations', iVal, iStat)
      call codes_get(ibufr, '#7#verticalSignificanceSurfaceObservations', iVal, iStat)
      call codes_get(ibufr, '#8#verticalSignificanceSurfaceObservations', iVal, iStat)
      call codes_get(ibufr, '#9#verticalSignificanceSurfaceObservations', iVal, iStat)
      call codes_get(ibufr, '#10#verticalSignificanceSurfaceObservations', iVal, iStat)
      call codes_get(ibufr, 'presentWeather', iVal, iStat)
      call codes_get(ibufr, '#1#timePeriod', iVal, iStat)
      call codes_get(ibufr, 'pastWeather1', iVal, iStat)
      call codes_get(ibufr, 'pastWeather2', iVal, iStat)
      call codes_get(ibufr, '#2#timePeriod', iVal, iStat)
      call codes_get(ibufr, '#3#timePeriod', iVal, iStat)
      call codes_get(ibufr, '#5#heightOfSensorAboveLocalGroundOrDeckOfMarinePlatform', rVal, iStat)
      call codes_get(ibufr, '#4#timePeriod', iVal, iStat)
      if (iStat == CODES_SUCCESS) then
         iperiod = iVal
      endif
      call codes_get(ibufr, '#1#totalPrecipitationOrTotalWaterEquivalent', rVal, iStat)
      if (iStat == CODES_SUCCESS .and. iperiod == -1 .and. rVal >= 0 .and. rVal < 1.e6) then
         rr     = rVal
      endif
      call codes_get(ibufr, '#5#timePeriod', iVal, iStat)
      if (iStat == CODES_SUCCESS) then
         iperiod = iVal
      endif
      call codes_get(ibufr, '#2#totalPrecipitationOrTotalWaterEquivalent', rVal, iStat)
      if (iStat == CODES_SUCCESS .and. iperiod == -1 .and. rVal >= 0 .and. rVal < 1.e6) then
         rr     = rVal
      endif
      call codes_get(ibufr, '#6#heightOfSensorAboveLocalGroundOrDeckOfMarinePlatform', rVal, iStat)
      call codes_get(ibufr, '#6#timePeriod', iVal, iStat)
      call codes_get(ibufr, '#7#timePeriod', iVal, iStat)
      call codes_get(ibufr, '#8#timePeriod', iVal, iStat)
      call codes_get(ibufr, '#9#timePeriod', iVal, iStat)
      call codes_get(ibufr, '#7#heightOfSensorAboveLocalGroundOrDeckOfMarinePlatform', rVal, iStat)
      call codes_get(ibufr, 'instrumentationForWindMeasurement', iVal, iStat)
      call codes_get(ibufr, '#1#timeSignificance', iVal, iStat)
      call codes_get(ibufr, '#10#timePeriod', iVal, iStat)
      call codes_get(ibufr, 'windDirection', iVal, iStat)
      if (iStat == CODES_SUCCESS .and. iVal >= 0 .and. iVal < 1000000) then
         winddir = iVal
      endif
      call codes_get(ibufr, 'windSpeed', rVal, iStat)
      if (iStat == CODES_SUCCESS .and. rVal >= 0 .and. rVal < 1.e6) then
         windspd = rVal
      endif
      call codes_get(ibufr, '#11#timePeriod', iVal, iStat)
      if (iStat == CODES_SUCCESS) then
         iperiod = iVal
      endif
      call codes_get(ibufr, '#1#maximumWindGustSpeed', rVal, iStat)
      if (iStat == CODES_SUCCESS .and. iperiod == -10 .and. rVal >= 0 .and. rVal < 1.e6) then
         gust = rVal
      endif
      call codes_get(ibufr, '#12#timePeriod', iVal, iStat)
      if (iStat == CODES_SUCCESS) then
         iperiod = iVal
      endif
      call codes_get(ibufr, '#2#maximumWindGustSpeed', rVal, iStat)
      if (iStat == CODES_SUCCESS .and. iperiod == -10 .and. rVal >= 0 .and. rVal < 1.e6) then
         gust = rVal
      endif
      call codes_get(ibufr, '#13#timePeriod', iVal, iStat)
      call codes_get(ibufr, '#14#timePeriod', iVal, iStat)
      call codes_get(ibufr, '#1#globalSolarRadiationIntegratedOverPeriodSpecified', rVal, iStat)
      call codes_get(ibufr, '#15#timePeriod', iVal, iStat)
      call codes_get(ibufr, '#2#globalSolarRadiationIntegratedOverPeriodSpecified', rVal, iStat)
      call codes_get(ibufr, '#16#timePeriod', iVal, iStat)
      call codes_get(ibufr, '#2#latitude', rVal, iStat)
      call codes_get(ibufr, '#2#longitude', rVal, iStat)
      call codes_get(ibufr, 'totalSnowDepth', rVal, iStat)
      if (iStat == CODES_SUCCESS .and. rVal >= 0 .and. rVal < 1.e6) then
         snow = rVal
      endif
   endif
   if (iobtype == 1) then ! Additional data from SHIP observations
      call codes_get(ibufr, 'seaSurfaceTemperature', rVal, iStat)
      if (iStat == CODES_SUCCESS .and. rVal >= 0 .and. rVal < 1.e6) then
         sst = rVal
      endif
   endif
   if (iobtype == 4) then
      call codes_get(ibufr, 'height', rVal, iStat)
      if (iStat == CODES_SUCCESS) then
         obheight = rVal
      else
         call codes_get(ibufr, 'flightLevel', rVal, iStat)
         if (iStat == CODES_SUCCESS) obheight = rVal
      endif
   else
      call codes_get(ibufr, 'heightOfStation', iVal, iStat)
      if (iStat == CODES_SUCCESS) obheight = iVal
   endif
   call codes_release(ibufr)

   if (iminute < 0 .or. iminute > 59) call abort
   if (iobtype /= 0 .and. iobtype /= 1 .and. iobtype /= 2 .and. iobtype /= 4) then
      write(0,*) 'DECODE-ECCODES: UNRECOGNISED OBSERVATION TYPE:', iobtype
      cycle
   endif

   if (obheight > 100000.) cycle               ! Missing data - cop out.

   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 = ihour + findloc([.false.,.true.],iminute>=30,dim=1) - 1
   ihour = mod(ihour, 24)

   if (iobtype == 0) nsynop = nsynop + 1
   if (iobtype == 1) nship  = nship  + 1
   if (iobtype == 2) ntemp  = ntemp  + 1
   if (iobtype == 4) nairep = nairep + 1
   inarea = inarea + 1
   iunit  = (iobtype + 1) * 100 + ihour
   if (.not. lopened(iunit)) then
      write(filename, '(a,i2.2,a)') trim(filetype(iobtype + 1)), ihour, '.txt'
      open(iunit, file = filename)
      lopened(iunit) = .true.
   endif
   if (iobtype == 0) then                  ! extra data obtained from SYNOP observations
      write(iunit, *) ilon, ilat, obheight, winddir, windspd, temp, dewp, rh, vis, cc, rr, gust, snow
   else
      write(iunit, *) ilon, ilat, obheight
   endif
   if      (iobtype == 0) then
      write(7,*) iobtype, iobsubtype, ident, oblat, oblon, obheight, ihour, iminute, temp, dewp, snow
   else if (iobtype == 1) then
      write(7,*) iobtype, iobsubtype, ident, oblat, oblon, obheight, ihour, iminute, sst, 1.7e38, 1.7e38
   endif

enddo
close(7)
call codes_close_file(ifile)

! Close output files

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

write(0,*)'DECODE-ECCODES: NUMBER OF OBSERVATIONS:',irec,', INSIDE AREA:',inarea, &
             ', SYNOPS:',nsynop,', SHIPS:',nship,', TEMPS:',ntemp,', AIREPS:',nairep

! The End

end program decode_eccodes