! (C) Copyright 2005- ECMWF.
!
! This software is licensed under the terms of the Apache Licence Version 2.0
! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
!
! In applying this licence, ECMWF does not waive the privileges and immunities granted to it by
! virtue of its status as an intergovernmental organisation nor does it submit to any jurisdiction.
!
!
!
! Description: How to read data for a tropical cyclone BUFR message.
!
! Please note that tropical cyclone tracks can be encoded in various ways in BUFR.
! Therefore the code below might not work directly for other types of messages
! than the one used in the example. It is advised to use bufr_dump to
! understand the structure of the messages.
!
program bufr_read_tropical_cyclone
  use eccodes
  implicit none
  integer            :: ifile
  integer            :: iret
  integer            :: ibufr, skipMember
  integer            :: significance
  integer            :: year, month, day, hour, minute
  integer            :: i, j, k, ierr, count = 1
  integer            :: rankPosition, rankSignificance, rankPressure, rankWind
  integer            :: rankPeriod, numberOfPeriods

  real(kind=8) :: latitudeCentre, longitudeCentre
  real(kind=8), dimension(:), allocatable :: latitudeMaxWind0, longitudeMaxWind0, windMaxWind0
  real(kind=8), dimension(:), allocatable :: latitudeAnalysis, longitudeAnalysis, pressureAnalysis
  real(kind=8), dimension(:, :), allocatable :: latitude, longitude, pressure
  real(kind=8), dimension(:, :), allocatable :: latitudeWind, longitudeWind, wind
  integer(kind=4), dimension(:), allocatable :: memberNumber, period
  real(kind=8), dimension(:), allocatable :: values
  integer(kind=4), dimension(:), allocatable :: ivalues

  character(len=8)   :: rankSignificanceStr, rankPositionStr, rankPressureStr, rankWindStr
  character(len=8)   :: stormIdentifier, rankPeriodStr

  call codes_open_file(ifile, '../../data/bufr/tropical_cyclone.bufr', 'r')

  ! The first BUFR message is loaded from file.
  ! ibufr is the BUFR id to be used in subsequent calls
  call codes_bufr_new_from_file(ifile, ibufr, iret)

  do while (iret /= CODES_END_OF_FILE)

    write (*, '(A,I3,A)') '**************** MESSAGE: ', count, '  *****************'

    ! We need to instruct ecCodes to unpack the data values
    call codes_set(ibufr, "unpack", 1);
    call codes_get(ibufr, 'year', year);
    call codes_get(ibufr, 'month', month);
    call codes_get(ibufr, 'day', day);
    call codes_get(ibufr, 'hour', hour);
    call codes_get(ibufr, 'minute', minute);
    write (*, '(A,I0,A,I0,A,I0,A,I0,A,I0,A,I0)') 'Date and time: ', day, '.', month, '.', year, '  ', hour, ':', minute

    call codes_get(ibufr, 'stormIdentifier', stormIdentifier)
    write (*, '(A,A)') 'Storm identifier: ', stormIdentifier

    ! How many different timePeriod in the data structure?
    rankPeriod = 0
    ierr = 0
    do while (ierr == 0)
      rankPeriod = rankPeriod + 1
      write (rankPeriodStr, '(I0)') rankPeriod
      call codes_get(ibufr, '#'//trim(rankPeriodStr)//'#timePeriod', period, ierr)
      if (allocated(period)) deallocate (period)
    end do
    ! The numberOfPeriods includes the analysis (period=0)
    numberOfPeriods = rankPeriod

    call codes_get(ibufr, 'ensembleMemberNumber', memberNumber)

    allocate (latitude(size(memberNumber), numberOfPeriods))
    allocate (longitude(size(memberNumber), numberOfPeriods))
    allocate (pressure(size(memberNumber), numberOfPeriods))
    allocate (latitudeWind(size(memberNumber), numberOfPeriods))
    allocate (longitudeWind(size(memberNumber), numberOfPeriods))
    allocate (wind(size(memberNumber), numberOfPeriods))
    allocate (values(size(memberNumber)))
    allocate (period(numberOfPeriods))
    period(1) = 0

    ! Observed Storm Centre
    call codes_get(ibufr, '#1#meteorologicalAttributeSignificance', significance);
    call codes_get(ibufr, '#1#latitude', latitudeCentre);
    call codes_get(ibufr, '#1#longitude', longitudeCentre);
    if (significance /= 1) then
      print *, 'ERROR: unexpected #1#meteorologicalAttributeSignificance'
      stop 1
    end if
    if (latitudeCentre == CODES_MISSING_DOUBLE .and. longitudeCentre == CODES_MISSING_DOUBLE) then
      write (*, '(a)') 'Observed storm centre position missing'
    else
      write (*, '(A,F8.2,A,F8.2)') 'Observed storm centre: latitude=', latitudeCentre, ' longitude=', longitudeCentre
    end if

    ! Location of storm in perturbed analysis
    call codes_get(ibufr, '#2#meteorologicalAttributeSignificance', significance);
    call codes_get(ibufr, '#2#latitude', latitudeAnalysis);
    call codes_get(ibufr, '#2#longitude', longitudeAnalysis);
    call codes_get(ibufr, '#1#pressureReducedToMeanSeaLevel', pressureAnalysis);
    if (significance /= 4) then
      print *, 'ERROR: unexpected #2#meteorologicalAttributeSignificance'
      stop 1
    end if
    if (size(latitudeAnalysis) == size(memberNumber)) then
      latitude(:, 1) = latitudeAnalysis
      longitude(:, 1) = longitudeAnalysis
      pressure(:, 1) = pressureAnalysis
    else
      latitude(:, 1) = latitudeAnalysis(1)
      longitude(:, 1) = longitudeAnalysis(1)
      pressure(:, 1) = pressureAnalysis(1)
    end if

    ! Location of Maximum Wind
    call codes_get(ibufr, '#3#meteorologicalAttributeSignificance', significance);
    call codes_get(ibufr, '#3#latitude', latitudeMaxWind0);
    call codes_get(ibufr, '#3#longitude', longitudeMaxWind0);
    if (significance /= 3) then
      print *, 'ERROR: unexpected #3#meteorologicalAttributeSignificance=', significance
      stop 1
    end if
    call codes_get(ibufr, '#1#windSpeedAt10M', windMaxWind0);
    if (size(latitudeMaxWind0) == size(memberNumber)) then
      latitudeWind(:, 1) = latitudeMaxWind0
      longitudeWind(:, 1) = longitudeMaxWind0
      wind(:, 1) = windMaxWind0
    else
      latitudeWind(:, 1) = latitudeMaxWind0(1)
      longitudeWind(:, 1) = longitudeMaxWind0(1)
      wind(:, 1) = windMaxWind0(1)
    end if

    rankSignificance = 3
    rankPosition = 3
    rankPressure = 1
    rankWind = 1
    rankPeriod = 0

    ! Loop on all periods excluding analysis period(1)=0
    do i = 2, numberOfPeriods

      rankPeriod = rankPeriod + 1
      write (rankPeriodStr, '(I0)') rankPeriod
      call codes_get(ibufr, '#'//trim(rankPeriodStr)//'#timePeriod', ivalues);
      do k = 1, size(ivalues)
        if (ivalues(k) /= CODES_MISSING_LONG) then
          period(i) = ivalues(k)
          exit
        end if
      end do
      deallocate (ivalues)

      ! Location of the storm
      rankSignificance = rankSignificance + 1
      write (rankSignificanceStr, '(I0)') rankSignificance
      call codes_get(ibufr, '#'//trim(rankSignificanceStr)//'#meteorologicalAttributeSignificance', ivalues);
      do k = 1, size(ivalues)
        if (ivalues(k) /= CODES_MISSING_LONG) then
          significance = ivalues(k)
          exit
        end if
      end do
      deallocate (ivalues)

      rankPosition = rankPosition + 1
      write (rankPositionStr, '(I0)') rankPosition
      call codes_get(ibufr, '#'//trim(rankPositionStr)//'#latitude', values);
      latitude(:, i) = values
      call codes_get(ibufr, '#'//trim(rankPositionStr)//'#longitude', values);
      longitude(:, i) = values

      if (significance == 1) then
        rankPressure = rankPressure + 1
        write (rankPressureStr, '(I0)') rankPressure
        call codes_get(ibufr, '#'//trim(rankPressureStr)//'#pressureReducedToMeanSeaLevel', values);
        pressure(:, i) = values
      else
        print *, 'ERROR: unexpected meteorologicalAttributeSignificance=', significance
        stop 1
      end if

      ! Location of maximum wind
      rankSignificance = rankSignificance + 1
      write (rankSignificanceStr, '(I0)') rankSignificance
      call codes_get(ibufr, '#'//trim(rankSignificanceStr)//'#meteorologicalAttributeSignificance', ivalues);
      do k = 1, size(ivalues)
        if (ivalues(k) /= CODES_MISSING_LONG) then
          significance = ivalues(k)
          exit
        end if
      end do
      deallocate (ivalues)

      rankPosition = rankPosition + 1
      write (rankPositionStr, '(I0)') rankPosition
      call codes_get(ibufr, '#'//trim(rankPositionStr)//'#latitude', values);
      latitudeWind(:, i) = values
      call codes_get(ibufr, '#'//trim(rankPositionStr)//'#longitude', values);
      longitudeWind(:, i) = values

      if (significance == 3) then
        rankWind = rankWind + 1
        write (rankWindStr, '(I0)') rankWind
        call codes_get(ibufr, '#'//trim(rankWindStr)//'#windSpeedAt10M', values);
        wind(:, i) = values
      else
        print *, 'ERROR: unexpected meteorologicalAttributeSignificance=,', significance
        stop 1
      end if

    end do

    ! Print the values
    do i = 1, size(memberNumber)
      skipMember = 1
      do j = 1, size(period)
        if (latitude(i, j) /= CODES_MISSING_DOUBLE .OR. latitudeWind(i, j) /= CODES_MISSING_DOUBLE) then
          skipMember = 0
          exit
        end if
      end do
      if (skipMember /= 1) then

        write (*, '(A,I3)') '== Member ', memberNumber(i)
        write (*, *) 'step  latitude   longitude   pressure    latitude   longitude   wind'
        do j = 1, size(period)
          if (latitude(i, j) /= CODES_MISSING_DOUBLE .OR. latitudeWind(i, j) /= CODES_MISSING_DOUBLE) then
            write (*, '( I4,2X,F8.2,4X,F8.2,3X,F9.1,2X,F8.2,4X,F8.2,2X,F8.2)') &
              period(j), latitude(i, j), longitude(i, j), pressure(i, j), &
              latitudeWind(i, j), longitudeWind(i, j), wind(i, j)
          end if
        end do

      end if
    end do

    ! deallocating the arrays is very important
    ! because the behaviour of the codes_get functions is as follows:
    ! if the array is not allocated then allocate
    ! if the array is already allocated only copy the values
    deallocate (values)
    deallocate (latitude)
    deallocate (longitude)
    deallocate (pressure)
    deallocate (latitudeWind)
    deallocate (longitudeWind)
    deallocate (wind)
    deallocate (period)
    deallocate (latitudeAnalysis)
    deallocate (longitudeAnalysis)
    deallocate (pressureAnalysis)
    deallocate (memberNumber)
    deallocate (latitudeMaxWind0)
    deallocate (longitudeMaxWind0)
    deallocate (windMaxWind0)

    ! Release the BUFR message
    call codes_release(ibufr)

    ! Load the next BUFR message
    call codes_bufr_new_from_file(ifile, ibufr, iret)

    count = count + 1

  end do

  ! Close file
  call codes_close_file(ifile)

end program bufr_read_tropical_cyclone