1! (C) Copyright 2005- ECMWF. 2! 3! This software is licensed under the terms of the Apache Licence Version 2.0 4! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. 5! 6! In applying this licence, ECMWF does not waive the privileges and immunities granted to it by 7! virtue of its status as an intergovernmental organisation nor does it submit to any jurisdiction. 8! 9! 10! Description: read and print radiosonde data from TEMP BUFR messages. 11! If available this version also lists the position information from the WMO list 12! (now OSCAR/Surface) appended to the reports by ECMWF 13! 14! Author: Bruce Ingleby 15! 16! Please note that TEMP reports can be encoded in various ways in BUFR. Therefore the code 17! below might not work directly for other types of TEMP messages than the one used in the 18! example. It is advised to use bufr_dump first to understand the structure of these messages. 19! 20program bufr_read_tempf 21 use eccodes 22 implicit none 23 integer :: ifile 24 integer :: iret 25 integer :: ibufr 26 integer :: i, count = 0 27 integer :: iflag 28 integer :: status_id, status_ht, status_time = 0, status_p 29 integer :: status_rsno, status_rssoft, status_balloonwt, statid_missing 30 integer(kind=4) :: sizews 31 integer(kind=4) :: blockNumber, stationNumber 32 integer(kind=4) :: ymd, hms 33 logical :: llstdonly = .True. ! Set True to list standard levels only 34 logical :: llskip 35 real(kind=8) :: year, month, day, hour, minute, second 36 real(kind=8) :: htg, htp, htec = 0, sondeType, balloonwt 37 ! real(kind=8), dimension(:), allocatable :: descriptors 38 real(kind=8), dimension(:), allocatable :: lat, lon 39 real(kind=8), dimension(:), allocatable :: timeVal, dlatVal, dlonVal, vssVal 40 real(kind=8), dimension(:), allocatable :: presVal, zVal, tVal, tdVal, wdirVal, wspVal 41 character(len=128) :: statid 42 character(len=16) :: rsnumber 43 character(len=16) :: rssoftware 44 45 call codes_open_file(ifile, '../../data/bufr/PraticaTemp.bufr', 'r') 46 47 ! the first BUFR message is loaded from file 48 ! ibufr is the BUFR id to be used in subsequent calls 49 call codes_bufr_new_from_file(ifile, ibufr, iret) 50 51 ! loop through all messages in the file 52 do while (iret /= CODES_END_OF_FILE .AND. status_time == CODES_SUCCESS) 53 54 ! Can check the template used 55 ! call codes_get(ibufr,'unexpandedDescriptors',descriptors) 56 ! write(0,*) 'Template: ', descriptors 57 ! IF (descriptors(1) /= 309056.0) GOTO 999 ! only list descent profiles 58 59 ! we need to instruct ecCodes to expand all the descriptors 60 ! i.e. unpack the data values 61 call codes_set(ibufr, "unpack", 1); 62 ! In our BUFR message verticalSoundingSignificance is always followed by 63 ! geopotential, airTemperature, dewpointTemperature, 64 ! windDirection, windSpeed and pressure. 65 66 count = count + 1 67 llskip = .False. 68 69 ! Metadata: 70 call codes_get(ibufr, 'shipOrMobileLandStationIdentifier', statid, status_id) 71 IF (status_id /= CODES_SUCCESS) statid = "UNKNOWN" 72 call codes_is_missing(ibufr, 'shipOrMobileLandStationIdentifier', statid_missing) 73 IF (statid_missing == 1) statid = "MISSING" 74 call codes_get(ibufr, 'blockNumber', blockNumber) 75 call codes_get(ibufr, 'stationNumber', stationNumber) 76 IF (blockNumber <= 99.0 .AND. stationNumber <= 1000) write (statid, '(I2.2,I3.3,3X)') blockNumber, stationNumber 77 78 call codes_get(ibufr, 'year', year) 79 call codes_get(ibufr, 'month', month) 80 call codes_get(ibufr, 'day', day) 81 call codes_get(ibufr, 'hour', hour) 82 call codes_get(ibufr, 'minute', minute) 83 call codes_get(ibufr, 'second', second, status_time) 84 IF (status_time /= CODES_SUCCESS) second = 0.0 85 call codes_get(ibufr, 'latitude', lat) 86 call codes_get(ibufr, 'longitude', lon) 87 call codes_get(ibufr, 'heightOfStationGroundAboveMeanSeaLevel', htg, status_ht) 88 IF (status_ht /= CODES_SUCCESS) htg = -999.0 89 call codes_get(ibufr, 'heightOfBarometerAboveMeanSeaLevel', htp, status_ht) 90 IF (status_ht /= CODES_SUCCESS) htp = -999.0 91 call codes_get(ibufr, 'radiosondeType', sondeType) 92 call codes_get(ibufr, 'heightOfStation', htec, status_ht) ! Height from WMO list (BUFR) 93 IF (status_ht == CODES_SUCCESS .AND. htg == -999.0) htg = htec 94 ymd = INT(year)*10000 + INT(month)*100 + INT(day) 95 hms = INT(hour)*10000 + INT(minute)*100 + INT(second) 96 call codes_get(ibufr, 'radiosondeSerialNumber', rsnumber, status_rsno) 97 call codes_get(ibufr, 'softwareVersionNumber', rssoftware, status_rssoft) 98 call codes_get(ibufr, 'weightOfBalloon', balloonwt, status_balloonwt) 99 IF (status_balloonwt /= CODES_SUCCESS) balloonwt = 0.0 100 101 ! Ascent (skip incomplete reports for now) 102 call codes_get(ibufr, 'timePeriod', timeVal, status_time) 103 IF (status_time /= CODES_SUCCESS) THEN 104 write (*, '(A,I7,A,A8,I9,I7.6,F9.3,F10.3,2F7.1,I4)') 'Ob: ', count, & 105 ' ', statid, ymd, hms, lat(1), lon(1), htg, htp, INT(sondeType) 106 write (*, '(A)') 'Missing times - skip' 107 llskip = .True. 108 END IF 109 call codes_get(ibufr, 'pressure', presVal, status_p) 110 IF (status_p /= CODES_SUCCESS) THEN 111 write (*, '(A,I7,A,A8,I9,I7.6,F9.3,F10.3,2F7.1,I4)') 'Ob: ', count, & 112 ' ', statid, ymd, hms, lat(1), lon(1), htg, htp, INT(sondeType) 113 write (*, '(A)') 'Missing pressures - skip' 114 llskip = .True. 115 END IF 116 call codes_get(ibufr, 'nonCoordinateGeopotentialHeight', zVal, status_ht) 117 IF (status_ht /= CODES_SUCCESS) THEN 118 write (*, '(A,I7,A,A8,I9,I7.6,F9.3,F10.3,2F7.1,I4)') 'Ob: ', count, & 119 ' ', statid, ymd, hms, lat(1), lon(1), htg, htp, INT(sondeType) 120 write (*, '(A)') 'Missing heights - skip' 121 llskip = .True. 122 END IF 123 ! IF (blockNumber /= 17 .OR. stationNumber /= 196) llskip=.True. ! FIX 124 ! IF (blockNumber /= 17.0) llskip=.True. ! FIX 125 126 IF (.NOT. llskip) THEN 127 call codes_get(ibufr, 'latitudeDisplacement', dlatVal) 128 call codes_get(ibufr, 'longitudeDisplacement', dlonVal) 129 call codes_get(ibufr, 'extendedVerticalSoundingSignificance', vssVal) 130 call codes_get(ibufr, 'airTemperature', tVal) 131 call codes_get(ibufr, 'dewpointTemperature', tdVal) 132 call codes_get(ibufr, 'windDirection', wdirVal) 133 call codes_get(ibufr, 'windSpeed', wspVal) 134 135 ! ---- Array sizes (pressure size can be larger - wind shear levels) 136 sizews = size(wspVal) 137 138 ! ---- Print the values -------------------------------- 139 write (*, '(A,I7,A,A8,I9,I7.6,F9.3,F10.3,2F7.1,I4,I5)') 'Ob: ', count, & 140 ' ', statid, ymd, hms, lat(1), lon(1), htg, htp, INT(sondeType), sizews 141 IF (status_rsno == CODES_SUCCESS) write (*, '(A,A,A,F7.3)') & 142 'RS number/software/balloonwt: ', rsnumber, rssoftware, balloonwt 143 IF (status_ht == CODES_SUCCESS .AND. SIZE(lat) > 1) write (*, '(A,F9.3,F10.3,F7.1)') & 144 'WMO list lat, lon, ht: ', lat(2), lon(2), htec 145 write (*, '(A)') 'level dtime dlat dlon pressure geopotH airTemp dewPtT windDir windSp signif' 146 do i = 1, sizews 147 iflag = vssVal(i) 148 IF (.NOT. llstdonly .OR. BTEST(iflag, 16)) & 149 write (*, '(I5,F7.1,2F7.3,F9.1,F8.1,4F8.2,I8)') i, timeVal(i), & 150 dlatVal(i), dlonVal(i), presVal(i), zVal(i), tVal(i), tdVal(i), & 151 wdirVal(i), wspVal(i), INT(vssVal(i)) 152 end do 153 154 END IF 155 ! free allocated arrays 156 IF (ALLOCATED(timeVal)) deallocate (timeVal) 157 IF (ALLOCATED(dlatVal)) deallocate (dlatVal) 158 IF (ALLOCATED(dlonVal)) deallocate (dlonVal) 159 IF (ALLOCATED(vssVal)) deallocate (vssVal) 160 IF (ALLOCATED(presVal)) deallocate (presVal) 161 IF (ALLOCATED(zVal)) deallocate (zVal) 162 IF (ALLOCATED(tVal)) deallocate (tVal) 163 IF (ALLOCATED(tdVal)) deallocate (tdVal) 164 IF (ALLOCATED(wdirVal)) deallocate (wdirVal) 165 IF (ALLOCATED(wspVal)) deallocate (wspVal) 166 IF (ALLOCATED(lat)) deallocate (lat) 167 IF (ALLOCATED(lon)) deallocate (lon) 168 169 ! 999 CONTINUE 170 ! release the BUFR message 171 call codes_release(ibufr) 172 173 ! load the next BUFR message 174 call codes_bufr_new_from_file(ifile, ibufr, iret) 175 176 end do 177 178 call codes_close_file(ifile) 179 180end program bufr_read_tempf 181