core
/
eccodes
mirror of https://github.com/ecmwf/eccodes.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

ECC-109

This commit is contained in:
Enrico Fucile 2015-11-18 14:40:34 +00:00
parent cd5ad41d46
commit dfb1863d63
4 changed files with 229 additions and 98 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
examples/F90/CMakeLists.txt
View File

@ -39,6 +39,7 @@ list( APPEND tests
bufr_read_scatterometer
bufr_read_synop
bufr_read_temp
bufr_read_tropical_cyclone
bufr_set_keys
bufr_subset
get_product_kind

5
examples/F90/Makefile.am
View File

@ -6,7 +6,8 @@ TESTS = copy_message.sh grib_get_keys.sh get_data.sh get_pl.sh get_pv.sh grib_ke
grib_set_bitmap.sh set_missing.sh grib_set_pv.sh samples.sh count_messages.sh read_message.sh \
read_from_file.sh grib_index.sh get_set_uuid.sh bufr_attributes.sh grib_clone.sh bufr_clone.sh \
bufr_expanded.sh bufr_get_keys.sh bufr_read_header.sh bufr_read_synop.sh \
bufr_set_keys.sh bufr_keys_iterator.sh bufr_subset.sh get_product_kind.sh bufr_read_temp.sh bufr_read_scatterometer.sh
bufr_set_keys.sh bufr_keys_iterator.sh bufr_subset.sh get_product_kind.sh bufr_read_temp.sh \
bufr_tropical_cyclone.sh bufr_read_scatterometer.sh
noinst_PROGRAMS = eccodes_f_grib_index \
eccodes_f_copy_message \
@ -44,6 +45,7 @@ noinst_PROGRAMS = eccodes_f_grib_index \
eccodes_f_bufr_attributes \
eccodes_f_get_product_kind \
eccodes_f_bufr_read_temp \
eccodes_f_bufr_tropical_cyclone \
eccodes_f_bufr_read_scatterometer
eccodes_f_grib_index_SOURCES=grib_index.f90
@ -83,6 +85,7 @@ eccodes_f_bufr_read_temp_SOURCES=bufr_read_temp.f90
eccodes_f_bufr_set_keys_SOURCES=bufr_set_keys.f90
eccodes_f_bufr_subset_SOURCES=bufr_subset.f90
eccodes_f_get_product_kind_SOURCES=get_product_kind.f90
eccodes_f_bufr_tropical_cyclone_SOURCES=bufr_tropical_cyclone.f90
INCLUDES = -I$(top_builddir)/src

319
examples/F90/bufr_read_tropical_cyclone.f90
View File

@ -1,4 +1,4 @@
!
!Copyright 2005-2015 ECMWF.
!
! This software is licensed under the terms of the Apache Licence Version 2.0
@ -18,132 +18,259 @@
program bufr_read_tropical_cyclone
use eccodes
implicit none
integer :: ifile
integer :: iret
integer :: ibufr
integer :: metAttSignificance,metAttSignificance3
integer :: year,month,day
integer :: i, ii,count=0
integer(kind=4) :: numObs,nTimeSteps,rank1,rank3
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), dimension(:), allocatable :: latVal,lonVal,pressure,wind,latVal3,lonVal3,time_steps
!integer, allocatable ::time_steps(:)
integer(kind=4), dimension(:), allocatable :: ensemble_members
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=32) :: key_name_rank1,key_name_rank3,keyname
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
! 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)') 'message: ',count
write(*,'(A,I3,A)') '**************** MESSAGE: ',count,' *****************'
! we need to instruct ecCodes to expand all the descriptors
! i.e. unpack the data values
! we need to instruct ecCodes to unpack the data values
call codes_set(ibufr,"unpack",1);
! The BUFR file contains a single message with 2016 subsets in a compressed form.
! It means each subset has exactly the same structure: they store one location with
! several beams and one backscatter value in each beam.
!
! To print the backScatter values for beamIdentifier=2 from all the subsets
! we will simply access the key by condition (see below).
! Read the total number of subsets.
call codes_get(ibufr,'numberOfSubsets',numObs)
write(*,'(A,I5)') "Number of Subsets:",numObs
call codes_get(ibufr,'year',year);
call codes_get(ibufr,'month',month);
call codes_get(ibufr,'day',day);
!allocate (time_steps(50))
call codes_get(ibufr,'timePeriod',time_steps)
call codes_get(ibufr,'ensembleMemberNumber',ensemble_members)
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
nTimeSteps=size(time_steps)
print*, 'time', time_steps,nTimeSteps
call codes_get(ibufr,'stormIdentifier',stormIdentifier)
write(*,'(A,A)')'Storm identifier: ',stormIdentifier
do ii=1,nTimeSteps
!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)
enddo
!the numberOfPeriods includes the analysis (period=0)
numberOfPeriods=rankPeriod
rank1=(ii)*2+2
rank3=(ii)*2+3
call codes_get(ibufr,'ensembleMemberNumber',memberNumber)
write (key_name_rank1,'(I0)')rank1
write (key_name_rank3,'(I0)')rank3
write (keyname,'(I0)')ii+1
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
!Get latitude
call codes_get(ibufr,'#'//trim(key_name_rank1)//'#latitude',latVal);
!Get longitude
call codes_get(ibufr,'#'//trim(key_name_rank1)//'#longitude',lonVal);
!Get pressure
call codes_get(ibufr,'#'//trim(keyname)//'#pressureReducedToMeanSeaLevel',pressure);
call codes_get(ibufr,'#'//trim(key_name_rank1)//'#meteorologicalAttributeSignificance',metAttSignificance);
call codes_get(ibufr,'#'//trim(key_name_rank3)//'#meteorologicalAttributeSignificance',metAttSignificance3);
!Get latitude
call codes_get(ibufr,'#'//trim(key_name_rank3)//'#latitude',latVal3);
!Get longitude
call codes_get(ibufr,'#'//trim(key_name_rank3)//'#longitude',lonVal3);
!Get pressure
call codes_get(ibufr,'#'//trim(keyname)//'#windSpeedAt10M',wind);
! ---- Check that all arrays are same size
if (size(latVal)/= numObs .or. size(lonVal)/= numObs) then
print *,'inconsistent array dimension'
exit
! 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
endif
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
endif
! 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
endif
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)
endif
! 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
endif
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)
endif
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
endif
enddo
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
endif
enddo
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
endif
!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
endif
enddo
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
endif
enddo
! ---- Print the values --------------------------------
write(*,*) "--------------------------------------"
write(*,'(A,I2,A,I2,A,I4,A,F3.0)')'Date:',day,'.',month,'.',year,' Time step:', time_steps(ii)
write(*,'(2(A,I1) )') 'metAttSignificance:', metAttSignificance,' metAttSignificance:', metAttSignificance3
write(*,*) 'ens lat lon pressure lat lon wind10m'
do i=1,numObs
write(*,'( I4,2(F6.2,A),F8.0,A,2(F6.2,A),F4.1 )') ensemble_members(i) ,latVal(i),' ',lonVal(i),' ',pressure(i),' '&
& ,latVal3(i),' ',lonVal3(i),' ',wind(i)
end do
! free arrays
deallocate(latVal)
deallocate(lonVal)
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
endif
enddo
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)
endif
enddo
endif
enddo
! free 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(latVal3)
deallocate(lonVal3)
deallocate(latitudeWind)
deallocate(longitudeWind)
deallocate(wind)
! if(allocated(time_steps))deallocate(time_steps)
! deallocate(time_steps)
!do i=1,15
!time_steps(i)=-9.
!end do
end do
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

2
examples/F90/bufr_read_tropical_cyclone.sh
View File

@ -18,7 +18,7 @@ label="bufr_read_tropical_cyclone_f"
fTmp=${label}".tmp.txt"
rm -f $fTmp | true
#We check "syno_multi.bufr". The path is
#We check "tropical_cyclone.bufr". The path is
#hardcoded in the example
REDIRECT=/dev/null