mirror of https://github.com/ecmwf/eccodes.git
278 lines
10 KiB
Fortran
278 lines
10 KiB
Fortran
! (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
|