Differences

This shows you the differences between two versions of the page.

--- bufr.pm:bufrdump [2012-01-27 15:41:35]
pals
+++ bufr.pm:bufrdump [2022-05-31 09:29:31] (current)
@@ Line 1: / Line 1: @@
-Use this in Makefile, with $(FC) set to your Fortran compiler (e.g. gfortran or g77), $(LDIR) set to the directory where libbufr.a is located, and $(FCFLAGS) set to -fbackslash for gfortran
+Use this in Makefile, with $(FC) set to your Fortran compiler (e.g. gfortran or g77), $(LDIR) set to the directory where libbufr.a is located, and $(FCFLAGS) set to -fbackslash for gfortran. Note that code for comfilter.f is found at the end of this page (below the code for bufrdump.F).
 <code>
@@ Line 5: / Line 5: @@
 	$(FC) $(FCFLAGS) -o bufrdump  $< -L $(LDIR) -lbufr
 </code>
 <code fortran bufrdump.F>
-C (C) Copyright 2010, met.no
+C (C) Copyright 2010-2016 MET Norway
 C
 C This program is free software; you can redistribute it and/or modify
@@ Line 24: / Line 25: @@
-C Extract BUFR messages from bufr file(s) using the Fortran program
+C Extract BUFR messages from bufr file(s) and print the data as
-C bufrdump and print the data as 'parameter=value' lines. See
+C 'parameter=value' lines. See usage_verbose for explanation of the
-C usage_verbose for explanation of the options allowed.
+C options allowed.
 C
 C Usage: bufrdump <bufr file>
 C        [--help]
 C        [--filter <filter file>]
-C        [--obstype] <amdar|ocea|surface|sounding>
+C        [--obstype] <amdar|ocea|surface|sounding|sounding->
 C        [--stop_on_error]
 C        [--rectangle x1 y1 x2 y2]
@@ Line 48: / Line 49: @@
 C export BUFR_TABLES=/usr/local/lib/bufrtables/
 C
-C Author: P.Sannes met.no 2010
+C Author: P.Sannes, MET Norway
@@ Line 160: / Line 161: @@
       INTEGER verbose
+      LOGICAL metar             ! Set to TRUE if metar (data subcategory 5)
       INTEGER kelem,kxelem      ! expected (max) number of expanded elements
       INTEGER kvals             ! expected (max) number of data values
@@ Line 207: / Line 209: @@
       kerr = 0
+      metar = .FALSE.
 C     Do a partial expansion without quality control
       kreq(1) = 1               ! All original elements without quality control
@@ Line 225: / Line 228: @@
 C
-C     Using parameter kelem in call to BUFREX might be too big for
+C     Using parameter kelem in call to BUFREX might be too big for some
-C     multisubset messages. Have copied the method used in decode_bufr.F
+C     multisubset messages (or too small for other messages). Have
-C     in libbufr, first calling BUS012 in order to get number of subsets
+C     copied the method used in decode_bufr.F in libbufr, first calling
-C     ksup(6)
+C     BUS012 in order to get number of subsets ksup(6)
       CALL BUS012(ilen,ibuff,ksup,ksec0,ksec1,ksec2,kerr)
       IF (kerr.NE.0) THEN
@@ Line 237: / Line 240: @@
       END IF
       kxelem = kvals/ksup(6)
-      IF (kxelem .GT. kelem) kxelem = kelem
+C     The second IF-condition is not in decode_bufr.F, but else I get 'KELEM
+C     ARGUMENT TOO SMALL' from BUUPWT when decoding large radiosondefiles
+      IF (kxelem.GT.kelem .AND. ksup(6).GT.10) kxelem = kelem
       CALL BUFREX (ilen,ibuff,ksup,ksec0,ksec1,ksec2,
@@ Line 252: / Line 257: @@
      +     'Number of subsets:',ksup(6)
-C     Convert messages with data category (BUFR table A) 0-2 and 4 only.
+C     Convert messages with data category (BUFR table A) 0-2,4,6 and 31 only.
 C     0 = Surface data - land, 1 = Surface data - sea, 2 = Vertical
 C     sounding (other than satellite) 4 = Single level upper-air data
-C     (other than satellite)
+C     (other than satellite) 6 = Radar data 31 = Oceanographic data
-      IF (ksec1(6).GT.2 .AND. ksec1(6).NE.4 .AND. ksec1(6).NE.31) RETURN
+      IF (ksec1(6).GT.2 .AND. ksec1(6).NE.4 .AND. ksec1(6).NE.6
+     +     .AND. ksec1(6).NE.31) RETURN
+C     MET (and perhaps also other ECMWF based software?) uses local
+C     subcategory 5 for metar (and BUFR edition 3)
+      IF (KSEC1(2).EQ.3 .AND. KSEC1(6).EQ.0 .AND. KSEC1(7).EQ.5) THEN
+         metar = .TRUE.
+      END IF
       IF (filter) THEN
@@ Line 292: / Line 304: @@
             IF (obstype.EQ.'surface') THEN
                CALL print_surface_values(ksub,kxelem,ktdexl,ktdexp,
-     +              values,cvals,rectangle,verbose)
+     +              values,cvals,rectangle,metar,verbose)
-            ELSE IF (obstype.EQ.'sounding') THEN
+            ELSE IF (obstype(1:8).EQ.'sounding') THEN
                CALL print_sounding_values(ksub,kxelem,ktdexl,ktdexp,
-     +              values,cvals,rectangle,verbose)
+     +              values,cvals,rectangle,obstype,verbose)
             ELSE IF (obstype.EQ.'amdar') THEN
                CALL print_amdar_values(ksub,kxelem,ktdexl,ktdexp,
@@ Line 304: / Line 316: @@
             ELSE IF (ksec1(6).LE.1) THEN
                CALL print_surface_values(ksub,kxelem,ktdexl,ktdexp,
-     +              values,cvals,rectangle,verbose)
+     +              values,cvals,rectangle,metar,verbose)
             ELSE IF (ksec1(6).EQ.2) THEN
                CALL print_sounding_values(ksub,kxelem,ktdexl,ktdexp,
-     +              values,cvals,rectangle,verbose)
+     +              values,cvals,rectangle,obstype,verbose)
             ELSE IF (ksec1(6).EQ.4) THEN
                CALL print_amdar_values(ksub,kxelem,ktdexl,ktdexp,
      +              values,cvals,rectangle,verbose)
+            ELSE IF (ksec1(6).EQ.6) THEN
+               CALL print_radar_profiler_values(ksub,kxelem,ktdexl,
+     +              ktdexp,values,cvals,rectangle,verbose)
             ELSE IF (ksec1(6).EQ.31) THEN
                CALL print_oceanographic_values(ksub,kxelem,ktdexl,
@@ Line 329: / Line 344: @@
             IF (obstype.EQ.'surface') THEN
                CALL print_surface_values(ksub,kxelem,ktdexl,ktdexp,
-     +              values,cvals,rectangle,verbose)
+     +              values,cvals,rectangle,metar,verbose)
-            ELSE IF (obstype.EQ.'sounding') THEN
+            ELSE IF (obstype(1:8).EQ.'sounding') THEN
                CALL print_sounding_values(ksub,kxelem,ktdexl,ktdexp,
-     +              values,cvals,rectangle,verbose)
+     +              values,cvals,rectangle,obstype,verbose)
             ELSE IF (obstype.EQ.'amdar') THEN
                CALL print_amdar_values(ksub,kxelem,ktdexl,ktdexp,
@@ Line 341: / Line 356: @@
             ELSE IF (ksec1(6).LE.1) THEN
                CALL print_surface_values(ksub,kxelem,ktdexl,ktdexp,
-     +              values,cvals,rectangle,verbose)
+     +              values,cvals,rectangle,metar,verbose)
             ELSE IF (ksec1(6).EQ.2) THEN
                CALL print_sounding_values(ksub,kxelem,ktdexl,ktdexp,
-     +              values,cvals,rectangle,verbose)
+     +              values,cvals,rectangle,obstype,verbose)
             ELSE IF (ksec1(6).EQ.4) THEN
                CALL print_amdar_values(ksub,kxelem,ktdexl,ktdexp,
      +              values,cvals,rectangle,verbose)
+            ELSE IF (ksec1(6).EQ.6) THEN
+               CALL print_radar_profiler_values(ksub,kxelem,ktdexl,
+     +              ktdexp,values,cvals,rectangle,verbose)
             ELSE IF (ksec1(6).EQ.31) THEN
 C     For ocea files at met.no data category is 31. But note that IOB
@@ Line 377: / Line 395: @@
       CHARACTER(LEN=200) argument ! Buffer for next argument
-      CHARACTER(LEN=1150) Usage
+      CHARACTER(LEN=1350) Usage
-      CHARACTER(LEN=900) Help
+      CHARACTER(LEN=950) Help
       INTEGER iargc,iarg
@@ Line 398: / Line 416: @@
      +     // '\t\t\tDecode observations meeting criteria in '
      +     // '<filter file> only\n'
-     +     // '\t--obstype <amdar|ocea|surface|sounding>\n'
+     +     // '\t--obstype <amdar|ocea|surface|sounding|sounding->\n'
      +     // '\t\t\tForce observation type. If this option '
      +     // 'is not set, bufrdump\n\t\t\twill make an '
      +     // 'educated guess of observation type based on\n'
-     +     // '\t\t\tmetadata in section 1 of each BUFR message\n'
+     +     // '\t\t\tmetadata in section 1 of each BUFR message. '
+     +     // 'The special value\n\t\t\t"sounding-" results in '
+     +     // 'levels with no vss being skipped\n'
      +     // '\t--stop_on_error\tIf a libbufr call returns an '
      +     // 'error, bufrdump will return\n\t\t\timmediately with '
@@ Line 439: / Line 459: @@
      +     // 'So far only filtering on exact matches on integer and\n'
      +     // 'character valued BUFR descriptors has been implemented.'
+     +     // ' But note that the\nclosely related program bufrdump.pl'
+     +     // ' has a lot more options for filtering.'
      +     // '\n\nIf an error occurs during decoding (typically '
      +     // 'because the required BUFR table\nis missing or message '
@@ Line 532: / Line 554: @@
       IF (obstype.NE.'' .AND. obstype.NE.'amdar' .AND. obstype.NE.'ocea'
-     +     .AND. obstype.NE.'surface' .AND. obstype.NE.'sounding') THEN
+     +     .AND. obstype.NE.'surface' .AND. obstype.NE.'sounding'
+     +     .AND. obstype.NE.'sounding-') THEN
          WRITE(*,*) 'Argument to --obstype must be one of amdar, ocea,',
-     +   ' surface or sounding'
+     +   ' surface, sounding or sounding-'
          CALL EXIT(1)
       END IF
@@ Line 544: / Line 567: @@
       SUBROUTINE err_msg(msg)
       IMPLICIT NONE
-      CHARACTER*80 msg
+      CHARACTER(LEN=*) msg
       WRITE(*,*) msg
       STOP
@@ Line 553: / Line 576: @@
       SUBROUTINE err_msg1(msg,kerr)
       IMPLICIT NONE
-      CHARACTER*80 msg
+      CHARACTER(LEN=*) msg
       INTEGER kerr
       WRITE(*,*) msg,kerr
@@ Line 655: / Line 678: @@
       CLOSE(11)
+      IF (verbose.GT.2) WRITE(*,*) 'nfidlines=',nfidlines
       IF (nfivlines.EQ.0) filter = .FALSE.
@@ Line 730: / Line 754: @@
       DO nsub = 1,num_subsets
 C     loop through all different conditions:
+         ifiv = 0
          DO i1 = 1,nfidlines
 C     loop through all filter value lines (for given) condition:
-            ifiv = 0
             DO ifvl = 1,nvl_fid(i1)
                ifiv = ifiv + 1
@@ Line 765: / Line 789: @@
                         END IF
                      END IF
-                   END DO        ! all data in subset
+                   END DO       ! all data in subset
                    IF (.NOT.match) ! there is one descriptor in condition which
-     +                  GO TO 300  ! subset does not match, so go to next value line
+     +                  EXIT       ! subset does not match, so go to next
-                END DO           ! all descriptors in condition
+                                   ! value line
+                END DO          ! all descriptors in condition
                 IF (match) THEN
                    num_found = num_found + 1
                    found_list(num_found) = nsub
                    GO TO 500    ! next subset
-               END IF
+                END IF
-           CONTINUE
             END DO              ! all filter value lines (for given) condition
         CONTINUE
@@ Line 786: / Line 810: @@
 C     Identify pressure, temperature etc and print parameter=value to
-C     screen. Note: have seen one example (high resolution data) where
+C     screen. Note: have seen high resolution data with 7002 HEIGHT OR
-C     7002 HEIGHT OR ALTITUDE was used as vertical coordinate instead of
+C     ALTITUDE and wind profiler data with 007009 GEOPOTENTIAL HEIGHT
-C     7004 PRESSURE (10004 was used for pressure). Not able to handle
+C     used as vertical coordinate instead of 7004 PRESSURE (10004 was
-C     that in present code.
+C     used for pressure in the first case). Not able to handle that in
+C     present code.
       SUBROUTINE print_sounding_values(ksub,kxelem,ktdexl,ktdexp,values,
-     +     cvals,rectangle,verbose)
+     +     cvals,rectangle,obstype,verbose)
       IMPLICIT NONE
@@ Line 801: / Line 826: @@
       CHARACTER*80 cvals(*)     ! Input: CCITTIA5 Bufr elements entries (one subset)
       LOGICAL rectangle         ! Input: TRUE if observations are wanted for a rectangle only
+      CHARACTER(LEN=*) obstype  ! Input: influence flow if set to 'sounding-'
       INTEGER verbose           ! Input: verbose level
@@ Line 807: / Line 833: @@
       REAL*8 II,iii,ix,year,month,day,hour,minute,
-     +     longitude,latitude,height
+     +     longitude,latitude,height,wigos_series
+      CHARACTER*5 wigos_issuer,wigos_issueno,missing5
       CHARACTER*9 call_sign,missing9
+      CHARACTER*8 flight_number,missing8
+      CHARACTER*16 wigos_localid,missing16
       CHARACTER one_bits
       REAL*8 value
-      INTEGER idx,cidx,desc,n,maxlevel,numlevels,ind
+      INTEGER idx,cidx,desc,n,maxlevel,numlevels,i,ind,ind2,ind3
-      PARAMETER(maxlevel=1000)
+      PARAMETER(maxlevel=100000)
       REAL*8 P(maxlevel),D(maxlevel),F(maxlevel),
@@ Line 819: / Line 849: @@
      +     wsb(maxlevel),wsa(maxlevel) ! absolute wind shear in 1 km layer below/above
       CHARACTER*8 vss(maxlevel) ! Vertical sounding significance
+      CHARACTER*8 vss_missing
+      LOGICAL reduce ! Set to TRUE if we should skip levels with
+                     ! vertical sounding significance missing or equal
+                     ! to 0
 C     Variables used for geographical filtering av observations
@@ Line 827: / Line 861: @@
       INTEGER lenstr
       CHARACTER*9 ctrim  ! length must be >= longest variable ctrim is used for
+      IF (obstype.EQ.'sounding-') THEN
+         reduce = .TRUE.
+      ELSE
+         reduce = .FALSE.
+      END IF
       one_bits = CHAR(255)
-      WRITE(missing9,'(9A)') one_bits,one_bits,one_bits,one_bits,
+      WRITE(missing5,'(5A)') one_bits,one_bits,one_bits,one_bits,
-     +     one_bits,one_bits,one_bits,one_bits,one_bits
+     +     one_bits
+      WRITE(missing8,'(8A)') one_bits,one_bits,one_bits,one_bits,
+     +     one_bits,one_bits,one_bits,one_bits
+      missing9 = missing8 // one_bits
+      missing16 = missing8 // missing8
+      vss_missing = '       '
 C     Initialize all parameters to missing values
       call_sign = missing9
+      flight_number = missing8
+      wigos_issuer = missing5
+      wigos_issueno = missing5
+      wigos_localid = missing16
+      wigos_series = rvind
       II = rvind
       iii= rvind
@@ Line 854: / Line 904: @@
          lo(n) = rvind
          tp(n) = rvind
-         vss(n) = '       '
+         vss(n) = vss_missing
          wsb(n) = rvind
          wsa(n) = rvind
@@ Line 875: / Line 925: @@
          ELSE IF (desc.EQ.1002) THEN ! WMO station number
             iii = value
-         ELSE IF (desc.EQ.1011) THEN  ! Ship or mobile land station identifier
+         ELSE IF (desc.EQ.1011.OR.   ! Ship or mobile land station identifier
+     +           desc.EQ.1195) THEN  ! Mobil land station identifier
                cidx = int(value/1000)
                call_sign = cvals(cidx) ! CCITTIA5 data
                call_sign = ctrim(call_sign,9,missing9)
+         ELSE IF (desc.EQ.1006) THEN ! Aircraft flight number, used for dropsondes
+            cidx = int(value/1000)
+            flight_number = cvals(cidx) ! CCITTIA5 data
+            flight_number = ctrim(flight_number,8,missing8)
+         ELSE IF (desc.EQ.1125) THEN ! WIGOS identifier series
+            wigos_series = value
+         ELSE IF (desc.EQ.1126) THEN ! WIGOS issuer of identifier
+            i = NINT(value)
+            WRITE(wigos_issuer,'(I5)') i
+            wigos_issuer = ADJUSTL(wigos_issuer)
+         ELSE IF (desc.EQ.1127) THEN ! WIGOS issue number
+            i = NINT(value)
+            WRITE(wigos_issueno,'(I5)') i
+            wigos_issueno = ADJUSTL(wigos_issueno)
+         ELSE IF (desc.EQ.1128) THEN  ! WIGOS local identifier
+            cidx = int(value/1000)
+            wigos_localid = cvals(cidx) ! CCITTIA5 data
+            wigos_localid = ctrim(wigos_localid,16,missing16)
          ELSE IF (desc.EQ.2001) THEN ! Type of station
             ix = value
@@ Line 900: / Line 969: @@
          ELSE IF (desc.EQ.6002 .AND. longitude.EQ.rvind) THEN ! Longitude (coarse accuracy)
             longitude = value
-         ELSE IF (desc.EQ.7001) THEN ! Height of station
+         ELSE IF (desc.EQ.7001.OR.  ! Height of station
-            height = value
+     +           desc.EQ.7030.OR.   ! Height of station ground above mean sea level
-         ELSE IF (desc.EQ.4086) THEN ! Long time period or displacement [second]
+     +           desc.EQ.7031.OR.   ! Height of barometer above mean sea level
+     +           desc.EQ.7007) THEN ! Height (i.e. of release of sonde above mean sea level)
+            IF (height.EQ.rvind) THEN
+               height = value
+            END IF
+         ELSE IF (desc.EQ.4086.AND.n.LT.maxlevel) THEN ! Long time period or displacement [second]
 C     In WMO template (309052) descriptors 004086 and 008042 comes BEFORE 7004 pressure
             tp(n+1) = value
-         ELSE IF (desc.EQ.8042) THEN ! Extended vertical sounding significance
+         ELSE IF (desc.EQ.8042.AND.n.LT.maxlevel) THEN ! Extended vertical sounding significance
             CALL vss_8042(NINT(value),vss(n+1))
          ELSE IF (desc.EQ.7004) THEN ! Pressure
             n = n + 1           ! new level
             IF (n.GT.maxlevel) THEN
+               n = maxlevel
                WRITE(*,*) 'Too many levels! Skipping rest of message'
                GOTO 110
@@ Line 915: / Line 990: @@
             P(n) = value
 C     All following descriptors come after pressure in 309052 (or Hirlam template)
-         ELSE IF (n.GT.0) THEN ! This condition secures against underflow
+         ELSE IF (n.GT.0 .AND. n.LE.maxlevel) THEN
             IF (desc.EQ.11001) THEN ! Wind direction
                D(n) = value
@@ Line 955: / Line 1030: @@
       END IF
-      WRITE(*,*)
       IF (II.NE.rvind .AND. iii.NE.rvind) THEN
+         WRITE(*,*)
          WRITE(*,'(A,I5.5)') 'wmonr=',NINT(II)*1000 + NINT(iii)
+      ELSE IF (wigos_series.NE.rvind .AND. wigos_issuer.NE.missing5
+     +        .AND. wigos_issueno.NE.missing5
+     +        .AND. wigos_localid.NE.missing16) THEN
+         ind = index(wigos_issuer,' ') - 1
+         IF (ind.EQ.-1) ind = 5
+         ind2 = index(wigos_issueno,' ') - 1
+         IF (ind2.EQ.-1) ind2 = 5
+         ind3 = index(wigos_localid,' ') - 1
+         IF (ind3.EQ.-1) ind3 = 16
+         WRITE(*,*)
+         WRITE(*,'(A,I1.1,A1,A,A1,A,A1,A)')
+     +        'wigosid=',NINT(wigos_series),
+     +        '-',wigos_issuer(1:ind),
+     +        '-',wigos_issueno(1:ind2),
+     +        '-',wigos_localid(1:ind3)
       ELSE IF (call_sign.NE.missing9) THEN
+         WRITE(*,*)
          WRITE(*,'(A,A)') 'call_sign=',call_sign(1:lenstr(call_sign,1))
+      ELSE IF (flight_number.NE.missing8) THEN
+         WRITE(*,*)
+         WRITE(*,'(A,A)') 'aircraft=',
+     +        flight_number(1:lenstr(flight_number,1))
       ELSE
          IF (verbose .GT. 1) THEN
-            WRITE(*,*) 'Both wmonr and call sign are missing!!!'
+            WRITE(*,*)
+            WRITE(*,*) 'Both wmonr, wigosid, call sign and aircraft'
+     +           // ' flight number are missing !!!'
          END IF
          RETURN
@@ Line 978: / Line 1076: @@
          END IF
       ENDIF
-      IF (NINT(ix).EQ.0) THEN
+      IF (ix.NE.rvind) THEN
-         WRITE(*,'(A,A)') 'type=Automatic'
+         IF (NINT(ix).EQ.0) THEN
-      ELSE IF (NINT(ix).EQ.1) THEN
+            WRITE(*,'(A,A)') 'type=Automatic'
-         WRITE(*,'(A,A)') 'type=Manned'
+         ELSE IF (NINT(ix).EQ.1) THEN
-      ELSE IF (NINT(ix).EQ.2) THEN
+            WRITE(*,'(A,A)') 'type=Manned'
-         WRITE(*,'(A,A)') 'type=Hybrid'
+         ELSE IF (NINT(ix).EQ.2) THEN
+            WRITE(*,'(A,A)') 'type=Hybrid'
+         END IF
       END IF
       IF (latitude.NE.rvind) THEN
@@ Line 996: / Line 1096: @@
       DO n=1,numlevels
+         IF (reduce .AND. vss(n).EQ.vss_missing) CYCLE
          WRITE(*,'(A,I12)'),'n=',n
 C     The following 3 parameters are not found in TAC TEMP
@@ Line 1036: / Line 1137: @@
       END SUBROUTINE print_sounding_values
+C     -----------------------------------------------------------------
+C     Identify height, wind and data quality parameters and print parameter=value to
+C     screen for the given heights.
+      SUBROUTINE print_radar_profiler_values(ksub,kxelem,ktdexl,ktdexp,
+     +     values,cvals,rectangle,verbose)
+      IMPLICIT NONE
+      INTEGER ksub              ! Input: number of subset currently processed
+      INTEGER kxelem            ! Input: expected (max) number of expanded elements
+      INTEGER ktdexl            ! Input: number of entries in list of expanded data descriptors
+      INTEGER ktdexp(*)         ! Input: array containing expanded data descriptors
+      REAL*8 values(*)          ! Input: expanded data values (one subset)
+      CHARACTER*80 cvals(*)     ! Input: CCITTIA5 Bufr elements entries (one subset)
+      LOGICAL rectangle         ! Input: TRUE if observations are wanted for a rectangle only
+      INTEGER verbose           ! Input: verbose level
+      REAL*8 rvind              ! missing value for real*8 data
+      PARAMETER (rvind=1.7E38)
+      REAL*8 II,iii,ix,year,month,day,hour,minute,
+     +     longitude,latitude,height
+      CHARACTER*9 call_sign,missing9
+      CHARACTER one_bits
+      REAL*8 value
+      INTEGER idx,cidx,desc,n,maxlevel,numlevels,ind
+      PARAMETER(maxlevel=1000)
+      REAL*8 HH(maxlevel),D(maxlevel),F(maxlevel),
+     +     Q1(maxlevel),WC(maxlevel),Q2(maxlevel)
+C     Variables used for geographical filtering av observations
+      REAL*8 x1,y1,x2,y2
+      COMMON /COM_RECTANGLE/  x1,y1,x2,y2
+C     Functions
+      INTEGER lenstr
+      CHARACTER*9 ctrim  ! length must be >= longest variable ctrim is used for
+      one_bits = CHAR(255)
+      WRITE(missing9,'(9A)') one_bits,one_bits,one_bits,one_bits,
+     +     one_bits,one_bits,one_bits,one_bits,one_bits
+C     Initialize all parameters to missing values
+      call_sign = missing9
+      II = rvind
+      iii= rvind
+      year = rvind
+      month = rvind
+      day = rvind
+      hour = rvind
+      minute = rvind
+      latitude = rvind
+      longitude = rvind
+      height = rvind
+      DO n=1,maxlevel
+         HH(n) = rvind
+         D(n) = rvind
+         F(n) = rvind
+         Q1(n) = rvind
+         WC(n) = rvind
+         Q2(n) = rvind
+      END DO
+C     Loop through all expanded descriptors
+      n = 0 ! Numbering the pressure levels
+      DO idx=1,ktdexl
+         desc = ktdexp(idx)
+         value = values(idx + (ksub-1)*kxelem)
+C     Nothing to do if value is missing (note that the missing value
+C     returned from libbufr might not be exactly equal to rvind)
+         IF (ABS(value - rvind)/rvind.LE.0.001) THEN
+            CYCLE
+         END IF
+         IF (desc.EQ.1001) THEN ! WMO block number
+            II = value
+         ELSE IF (desc.EQ.1002) THEN ! WMO station number
+            iii = value
+         ELSE IF (desc.EQ.1011.OR.   ! Ship or mobile land station identifier
+     +           desc.EQ.1195) THEN  ! Mobil land station identifier
+               cidx = int(value/1000)
+               call_sign = cvals(cidx) ! CCITTIA5 data
+               call_sign = ctrim(call_sign,9,missing9)
+         ELSE IF (desc.EQ.2001) THEN ! Type of station
+            ix = value
+         ELSE IF (desc.EQ.4001) THEN ! Year
+            year = value
+         ELSE IF (desc.EQ.4002) THEN ! Month
+            month = value
+         ELSE IF (desc.EQ.4003) THEN ! Day
+            day = value
+         ELSE IF (desc.EQ.4004) THEN ! Hour
+            hour = value
+         ELSE IF (desc.EQ.4005) THEN ! Minute
+            minute = value
+C     Seen 5002 and 6002 used for each level in high resolution data
+         ELSE IF (desc.EQ.5001 .AND. latitude.EQ.rvind) THEN ! Latitude (high accuracy)
+            latitude = value
+         ELSE IF (desc.EQ.5002 .AND. latitude.EQ.rvind) THEN ! Latitude (coarse accuracy)
+            latitude = value
+         ELSE IF (desc.EQ.6001 .AND. longitude.EQ.rvind) THEN ! Longitude (high accuracy)
+            longitude = value
+         ELSE IF (desc.EQ.6002 .AND. longitude.EQ.rvind) THEN ! Longitude (coarse accuracy)
+            longitude = value
+         ELSE IF (desc.EQ.7001) THEN  ! Height of station
+            IF (height.EQ.rvind) THEN
+               height = value
+            END IF
+         ELSE IF (desc.EQ.7007) THEN ! Height
+            n = n + 1           ! new level
+            IF (n.GT.maxlevel) THEN
+               n = maxlevel
+               WRITE(*,*) 'Too many levels! Skipping rest of message'
+               GOTO 110
+            END IF
+            HH(n) = value
+C     All following descriptors come after 7007 height
+         ELSE IF (n.GT.0 .AND. n.LE.maxlevel) THEN
+            IF (desc.EQ.11001) THEN ! Wind direction
+               D(n) = value
+            ELSE IF (desc.EQ.11002) THEN ! Wind speed
+               F(n) = value
+            ELSE IF (desc.EQ.11006) THEN ! W-component
+               WC(n) = value
+            ELSE IF (desc.EQ.33002) THEN ! Quality information, included 2 times per level
+               IF (ktdexp(idx-1).EQ.11002) THEN
+                  Q1(n) = value
+               ELSE
+                  Q2(n) = value
+               END IF
+            END IF
+         END IF
+      END DO
+  CONTINUE
+      numlevels = n
+      IF (rectangle) THEN
+         IF (longitude.EQ.rvind .OR. latitude.EQ.rvind
+     +        .OR. longitude.LT.x1 .OR. longitude.GT.x2
+     +        .OR. latitude.LT.y1 .OR. latitude.GT.y2) RETURN
+      END IF
+      IF (II.NE.rvind .AND. iii.NE.rvind) THEN
+         WRITE(*,*)
+         WRITE(*,'(A,I5.5)') 'wmonr=',NINT(II)*1000 + NINT(iii)
+      ELSE IF (call_sign.NE.missing9) THEN
+         WRITE(*,*)
+         WRITE(*,'(A,A)') 'call_sign=',call_sign(1:lenstr(call_sign,1))
+      ELSE
+         IF (verbose .GT. 1) THEN
+            WRITE(*,*)
+            WRITE(*,*) 'Both wmonr and call sign are missing!!!'
+         END IF
+         RETURN
+      END IF
+      IF (year.NE.rvind.AND.month.NE.rvind.AND.day.NE.rvind
+     +     .AND.hour.NE.rvind.AND.minute.NE.rvind) THEN
+         WRITE(*,900),NINT(year),NINT(month),NINT(day),
+     +        NINT(hour),NINT(minute)
+     FORMAT('obstime=',I4,'-',I2.2,'-',I2.2,' ',I2.2,':',I2.2,':00')
+      ELSE
+         IF (verbose .GT. 1) THEN
+            WRITE(*,*) 'obstime is missing!!!'
+            RETURN
+         END IF
+      ENDIF
+      IF (ix.NE.rvind) THEN
+         IF (NINT(ix).EQ.0) THEN
+            WRITE(*,'(A,A)') 'type=Automatic'
+         ELSE IF (NINT(ix).EQ.1) THEN
+            WRITE(*,'(A,A)') 'type=Manned'
+         ELSE IF (NINT(ix).EQ.2) THEN
+            WRITE(*,'(A,A)') 'type=Hybrid'
+         END IF
+      END IF
+      IF (latitude.NE.rvind) THEN
+         WRITE(*,'(A,F10.5)') 'lat=',latitude
+      END IF
+      IF (longitude.NE.rvind) THEN
+         WRITE(*,'(A,F10.5)') 'lon=',longitude
+      END IF
+      IF (height.NE.rvind) THEN
+         WRITE(*,'(A,I7)') 'height=',NINT(height)
+      END IF
+      DO n=1,numlevels
+         WRITE(*,'(A,I12)'),'n=',n
+         IF (HH(n).NE.rvind) THEN
+            WRITE(*,'(A,F11.1)') 'HH=',HH(n) ! m
+         END IF
+         IF (D(n).NE.rvind) THEN
+            WRITE(*,'(A,I11)') 'DD=',NINT(D(n))
+         END IF
+         IF (F(n).NE.rvind) THEN
+            WRITE(*,'(A,F11.1)') 'FF=',F(n)
+         END IF
+         IF (Q1(n).NE.rvind) THEN
+            WRITE(*,'(A,I11)') 'Q1=',NINT(Q1(n))
+         END IF
+         IF (WC(n).NE.rvind) THEN
+            WRITE(*,'(A,I11)') 'WC=',NINT(WC(n))
+         END IF
+         IF (Q2(n).NE.rvind) THEN
+            WRITE(*,'(A,I11)') 'Q2=',NINT(Q2(n))
+         END IF
+      END DO
+      END SUBROUTINE print_radar_profiler_values
 C     -----------------------------------------------------------------
       SUBROUTINE print_surface_values(ksub,kxelem,ktdexl,ktdexp,values,
-     +     cvals,rectangle,verbose)
+     +     cvals,rectangle,metar,verbose)
 C     Identify pressure, temperature etc and print parameter=value to screen
       IMPLICIT NONE
@@ Line 1051: / Line 1364: @@
       CHARACTER*80 cvals(*)     ! Input: CCITTIA5 Bufr elements entries (one subset)
       LOGICAL rectangle         ! Input: TRUE if observations are wanted for a rectangle only
+      LOGICAL metar             ! Input: TRUE if metar (data subcategory 5)
       INTEGER verbose           ! Input: verbose level
@@ Line 1056: / Line 1370: @@
       PARAMETER (rvind=1.7E38)
-      CHARACTER*8 icao_id,missing8,spc8
+      CHARACTER*5 wigos_issuer,wigos_issueno,missing5
-      CHARACTER*9 call_sign,missing9,spc9
+      CHARACTER*8 icao_id,missing8
-      CHARACTER*20 name,missing20,spc20
+      CHARACTER*9 call_sign,missing9
-      CHARACTER*32 long_name,missing32,spc32
+      CHARACTER*16 wigos_localid,missing16
+      CHARACTER*20 name,missing20
+      CHARACTER*32 long_name,missing32
 C     Parameters defined in Kvalobs
       REAL*8 AA,BI,CH,CI,CL,CM,DD,DG,DG_010,DG_1,DG_X,DI,DW1,DW2,
-     +     E,EM,ES,EV_1,EV_24,FF,FG,FG_010,FG_1,FG_X,FX,FX_1,FX_X,HL,
+     +     EE,ES,EV_1,EV_24,FF,FG,FG_010,FG_1,FG_X,FX,FX_1,FX_X,HL,
-     +     HW,HW1,HW2,NH,NN,OT_1,OT_24,PH,PO,PP,PR,PW,PW1,PW2,
+     +     HW,HW1,HW2,HWA,NH,NN,OT_1,OT_24,PH,PO,PP,PR,PW,PW1,PW2,PWA,
      +     QD,QE,QK,QL,QO,QS,QD_24,QE_24,QK_24,QL_24,QO_24,QS_24,
-     +     RR_1,RR_3,RR_6,RR_12,RR_24,RS,SA,SG,SI,
+     +     RR_1,RR_3,RR_6,RR_12,RR_24,RS,SA,SG,SI,SS_24,
      +     TA,TAN_12,TAX_12,TAN,TAX,TD,TGN_12,TW,UU,VV,W1,W2,WW,XIS,ZI,
 C     Other parameters
-     +     year,month,day,hour,minute,
+     +     year,month,day,hour,minute,a3,buoy_id5,buoy_id7,ds,
-     +     a3,buoy_id,ds,height,hhh,hour_p,II,iii,ix,latitude,longitude,
+     +     height,hhh,hp,hour_p,II,iii,ix,latitude,longitude,
      +     minute_p,TbTbTb,vert_sign_first,vs,wmo_region_number,
-     +     wmo_region_subarea,state_id,national_number
+     +     wmo_region_subarea,state_id,national_number,
+     +     wigos_series
       REAL*8 vert_sign(4),CC(4),HS(4),NS(4)
       INTEGER idx,cidx
       INTEGER cloud_type_count  ! Will be increased by one for each 020012
-                                ! (cloud type) encountered (0 initially)
+                                ! (cloud type) encountered (0 initially).
-      INTEGER num_cloud_layers  ! Number of individual cloud layers,
+                                ! Not used for metar
-                                ! set to value of 031001 (delayed
+      INTEGER num_cloud_layers  ! Number of individual cloud layers, set
-                                ! descriptor) if this is met immediately
+                                ! to value of 031001 (delayed descriptor
-                                ! after a 020012 descriptor (-1 initially)
+                                ! replication factor) if this is met immediately
+                                ! after a 020012 descriptor (-1 initially).
+                                ! For metar num_cloud_layers is increased
+                                ! by one for each new 020011
       LOGICAL bad_cloud_data    ! Set to true if something serious wrong is
                                 ! found in cloud data. No more cloud
@@ Line 1087: / Line 1407: @@
                                 ! surface') is encountered with a value different
                                 ! from 0
+      LOGICAL time_of_last_position ! Set to true if 008021 time significance is
+                                    ! included with value 26
       CHARACTER one_bits
       REAL*8 value
-      INTEGER desc,i,mm,hh,ind
+      INTEGER desc,i,mm,hh,ind,ind2,ind3
-      INTEGER degree2dir,HStoWMO_HSHS,NNtoWMO_N
+      INTEGER degree2dir,NNtoWMO_N
 C     Variables used for geographical filtering av observations
@@ Line 1102: / Line 1424: @@
       one_bits = CHAR(255)
+      WRITE(missing5,'(5A)') one_bits,one_bits,one_bits,one_bits,
+     +     one_bits
       WRITE(missing8,'(8A)') one_bits,one_bits,one_bits,one_bits,
      +     one_bits,one_bits,one_bits,one_bits
       missing9 = missing8 // one_bits
+      missing16 = missing8 // missing8
       missing20 = missing9 // missing9 // one_bits // one_bits
-      missing32 = missing20 // missing9 // one_bits // one_bits
+      missing32 = missing16 // missing16
-     +     // one_bits
-      spc8 = '        '
-      spc9 = '         '
-      spc20 = '                    '
-      spc32 = '                                '
       cloud_type_count = 0
@@ Line 1117: / Line 1437: @@
       num_cloud_layers = -1
       surface_data = .TRUE.
+      time_of_last_position = .FALSE.
 C     Initialize all parameters to missing values
@@ Line 1123: / Line 1444: @@
       icao_id = missing8
       call_sign = missing9
+      wigos_issuer = missing5
+      wigos_issueno = missing5
+      wigos_localid = missing16
       name = missing20
       long_name = missing32
@@ Line 1139: / Line 1463: @@
       DW1 = rvind
       DW2 = rvind
+      EE = rvind
       ES = rvind
-      E = rvind
-      EM = rvind
       EV_1 = rvind
       EV_24 = rvind
@@ Line 1156: / Line 1479: @@
       HW1 = rvind
       HW2 = rvind
+      HWA = rvind
       NH = rvind
       NN = rvind
@@ Line 1168: / Line 1492: @@
       PW1 = rvind
       PW2 = rvind
+      PWA = rvind
       QD = rvind
       QE = rvind
@@ Line 1189: / Line 1514: @@
       SG = rvind
       SI = rvind
+      SS_24 = rvind
       TA = rvind
       TAN_12 = rvind
@@ Line 1212: / Line 1538: @@
       longitude = rvind
       height = rvind
+      hp = rvind
       vert_sign_first = rvind
       II = rvind
@@ Line 1221: / Line 1548: @@
       vs = rvind
       TbTbTb = rvind
-      buoy_id = rvind
+      buoy_id5 = rvind
+      buoy_id7 = rvind
       wmo_region_number = rvind
       wmo_region_subarea = rvind
       state_id = rvind
       national_number = rvind
+      wigos_series = rvind
       DO i=1,4
@@ Line 1278: / Line 1607: @@
          ELSE IF (desc.EQ.4025) THEN ! Time period or displacement [minute]
             minute_p = value
+         ELSE IF (desc.EQ.8021) THEN ! Time significance
+            IF (NINT(value).EQ.26) THEN
+               time_of_last_position = .TRUE.
+            ELSE
+               time_of_last_position = .FALSE.
+            END IF
          ELSE IF (desc.EQ.1001 .AND. II.EQ.rvind) THEN ! WMO block number
-            II = value
+            IF (value.GE.0 .AND. value.LT.100) II = value
          ELSE IF (desc.EQ.1002 .AND. iii.EQ.rvind) THEN ! WMO station number
-            iii = value
+            IF (value.GE.0 .AND. value.LT.1000) iii = value
          ELSE IF (desc.EQ.1101 .AND. state_id.EQ.rvind) THEN ! WMO member state identifier
-            state_id = value
+            IF (value.GE.0 .AND. value.LT.1000) state_id = value
          ELSE IF (desc.EQ.1102 .AND. national_number.EQ.rvind) THEN ! National station number
-            national_number = value
+            IF (value.GE.0) national_number = value
          ELSE IF (desc.EQ.1015) THEN  ! Station or site name
             cidx = int(value/1000)
-            IF (cvals(cidx).NE.spc20) THEN
+            name = cvals(cidx)  ! CCITTIA5 data
-               name = cvals(cidx) ! CCITTIA5 data
+            name = ctrim(name,20,missing20)
-               name = ctrim(name,20,missing20)
-            END IF
          ELSE IF (desc.EQ.1019) THEN  ! Long station or site name
             cidx = int(value/1000)
-            IF (cvals(cidx).NE.spc32) THEN
+            long_name = cvals(cidx) ! CCITTIA5 data
-               long_name = cvals(cidx) ! CCITTIA5 data
+            long_name = ctrim(long_name,32,missing32)
-               long_name = ctrim(long_name,32,missing32)
+         ELSE IF (desc.EQ.1125) THEN ! WIGOS identifier series
-            END IF
+            wigos_series = value
+         ELSE IF (desc.EQ.1126) THEN ! WIGOS issuer of identifier
+            i = NINT(value)
+            WRITE(wigos_issuer,'(I5)') i
+            wigos_issuer = ADJUSTL(wigos_issuer)
+         ELSE IF (desc.EQ.1127) THEN ! WIGOS issue number
+            i = NINT(value)
+            WRITE(wigos_issueno,'(I5)') i
+            wigos_issueno = ADJUSTL(wigos_issueno)
+         ELSE IF (desc.EQ.1128) THEN  ! WIGOS local identifier
+            cidx = int(value/1000)
+            wigos_localid = cvals(cidx) ! CCITTIA5 data
+            wigos_localid = ctrim(wigos_localid,16,missing16)
          ELSE IF (desc.EQ.2001) THEN ! Type of station
             IF (ix.EQ.rvind) THEN
@@ Line 1303: / Line 1648: @@
             END IF
          ELSE IF (desc.EQ.4001) THEN ! Year
-            IF (year.EQ.rvind) THEN
+            IF (year.EQ.rvind.AND..NOT.time_of_last_position) THEN
                year = value
             END IF
          ELSE IF (desc.EQ.4002) THEN ! Month
-            IF (month.EQ.rvind) THEN
+            IF (month.EQ.rvind.AND..NOT.time_of_last_position) THEN
                month = value
             END IF
          ELSE IF (desc.EQ.4003) THEN ! Day
-            IF (day.EQ.rvind) THEN
+            IF (day.EQ.rvind.AND..NOT.time_of_last_position) THEN
                day = value
             END IF
          ELSE IF (desc.EQ.4004) THEN ! Hour
-            IF (hour.EQ.rvind) THEN
+            IF (hour.EQ.rvind.AND..NOT.time_of_last_position) THEN
                hour = value
             END IF
          ELSE IF (desc.EQ.4005) THEN ! Minute
-            IF (minute.EQ.rvind) THEN
+            IF (minute.EQ.rvind.AND..NOT.time_of_last_position) THEN
                minute = value
             END IF
-         ELSE IF (desc.EQ.5001) THEN ! Latitude (high accuracy)
+         ELSE IF (desc.EQ.5001.OR.  ! Latitude (high accuracy)
+     +           desc.EQ.5002) THEN ! Latitude (coarse accuracy)
             IF (latitude.EQ.rvind) THEN
                latitude = value
             END IF
-         ELSE IF (desc.EQ.5002) THEN ! Latitude (coarse accuracy)
+         ELSE IF (desc.EQ.6001.OR.  ! Longitude (high accuracy)
-            IF (latitude.EQ.rvind) THEN
+     +           desc.EQ.6002) THEN ! Longitude (coarse accuracy)
-               latitude = value
-            END IF
-         ELSE IF (desc.EQ.6001) THEN ! Longitude (high accuracy)
             IF (longitude.EQ.rvind) THEN
                longitude = value
             END IF
-         ELSE IF (desc.EQ.6002) THEN ! Longitude (coarse accuracy)
+         ELSE IF (desc.EQ.7001.OR.  ! Height of station
-            IF (longitude.EQ.rvind) THEN
+     +           desc.EQ.7030) THEN ! Height of station ground above mean sea level
-               longitude = value
-            END IF
-         ELSE IF (desc.EQ.7001) THEN ! Height of station
             IF (height.EQ.rvind) THEN
                height = value
             END IF
-         ELSE IF (desc.EQ.7030) THEN ! Height of station ground above mean sea level
+         ELSE IF (desc.EQ.7031) THEN ! Hp
-            IF (height.EQ.rvind) THEN
+            IF (hp.EQ.rvind) THEN
-               height = value
+               hp = value
             END IF
          ELSE IF (desc.EQ.10004) THEN ! Pressure
@@ Line 1362: / Line 1702: @@
                AA = value
             END IF
-         ELSE IF (desc.EQ.11011) THEN ! Wind direction at 10 m
+         ELSE IF (desc.EQ.11011.OR.  ! Wind direction at 10 m
+     +           desc.EQ.11001) THEN ! Wind direction
             IF (DD.EQ.rvind) THEN
                DD = value
             END IF
-         ELSE IF (desc.EQ.11001) THEN ! Wind direction
+         ELSE IF (desc.EQ.11012.OR.  ! Wind speed at 10 m
-            IF (DD.EQ.rvind) THEN
+     +           desc.EQ.11002) THEN ! Wind speed
-               DD = value
-            END IF
-         ELSE IF (desc.EQ.11012) THEN ! Wind speed at 10 m
-            IF (FF.EQ.rvind) THEN
-               FF = value
-            END IF
-         ELSE IF (desc.EQ.11002) THEN ! Wind speed
             IF (FF.EQ.rvind) THEN
                FF = value
@@ Line 1380: / Line 1714: @@
          ELSE IF (desc.EQ.11043) THEN ! Maximum wind gust direction
 C     DG is treated same way as FG
-            IF (value.NE.rvind
+            IF (minute_p.NE.rvind .AND. hour.NE.rvind) THEN
-     +           .AND. minute_p.NE.rvind .AND.hour.NE.rvind) THEN
                mm = NINT(minute_p)
                hh = NINT(hour)
@@ Line 1399: / Line 1732: @@
             END IF
          ELSE IF (desc.EQ.11041) THEN ! Maximum wind gust speed
-            IF (value.NE.rvind
+            IF (minute_p.NE.rvind .AND. hour.NE.rvind) THEN
-     +           .AND. minute_p.NE.rvind .AND. hour.NE.rvind) THEN
                mm = NINT(minute_p)
                hh = NINT(hour)
@@ Line 1423: / Line 1755: @@
             END IF
          ELSE IF (desc.EQ.11042) THEN ! Maximum wind speed (10-min mean wind)
-            IF (value.NE.rvind
+            IF (minute_p.NE.rvind .AND. hour.NE.rvind) THEN
-     +           .AND. minute_p.NE.rvind .AND. hour.NE.rvind) THEN
 C     FX is "Vindhastighet, maks. 10 minutt glidende middel siden
 C     forrige hovedobservasjon, m/s". Assume this is reported for
@@ Line 1441: / Line 1772: @@
                END IF
             END IF
-         ELSE IF (desc.EQ.12104) THEN ! Dry bulb temperature at 2m (data width 16 bits)
+         ELSE IF (desc.EQ.12104.OR. ! Dry bulb temperature at 2m (data width 16 bits)
+     +          desc.EQ.12004.OR.   ! Dry bulb temperature at 2m (12 bits)
+     +          desc.EQ.12101.OR.   ! Temperature/dry bulb temperature (16 bits)
+     +          desc.EQ.12001) THEN ! Temperature/dry bulb temperature (12 bits)
             IF (TA.EQ.rvind) THEN
                TA = value
             END IF
-         ELSE IF (desc.EQ.12004) THEN ! Dry bulb temperature at 2m (12 bits)
+         ELSE IF (desc.EQ.12106.OR.  ! Dew-point temperature at 2m (16 bits)
-            IF (TA.EQ.rvind) THEN
+     +           desc.EQ.12006.OR.   ! Dew-point temperature at 2m (12 bits)
-               TA = value
+     +           desc.EQ.12103.OR.   ! Dew-point temperature (16 bits)
-            END IF
+     +           desc.EQ.12003) THEN ! Dew-point temperature (12 bits)
-         ELSE IF (desc.EQ.12101) THEN ! Temperature/dry bulb temperature (16 bits)
-            IF (TA.EQ.rvind) THEN
-               TA = value
-            END IF
-         ELSE IF (desc.EQ.12001) THEN ! Temperature/dry bulb temperature (12 bits)
-            IF (TA.EQ.rvind) THEN
-               TA = value
-            END IF
-         ELSE IF (desc.EQ.12106) THEN ! Dew-point temperature at 2m (16 bits)
             IF (TD.EQ.rvind) THEN
                TD = value
             END IF
-         ELSE IF (desc.EQ.12006) THEN ! Dew-point temperature at 2m (12 bits)
+         ELSE IF (desc.EQ.12113.OR.  ! Ground minimum temperature at 2m (data width 16 bits)
-            IF (TD.EQ.rvind) THEN
+     +           desc.EQ.12013) THEN ! Ground minimum temperature at 2m (12 bits)
-               TD = value
-            END IF
-         ELSE IF (desc.EQ.12103) THEN ! Dew-point temperature (16 bits)
-            IF (TD.EQ.rvind) THEN
-               TD = value
-            END IF
-         ELSE IF (desc.EQ.12003) THEN ! Dew-point temperature (12 bits)
-            IF (TD.EQ.rvind) THEN
-               TD = value
-            END IF
-         ELSE IF (desc.EQ.12113) THEN ! Ground minimum temperature at 2m (data width 16 bits)
             IF (TGN_12.EQ.rvind) THEN
                TGN_12 = value
             END IF
-         ELSE IF (desc.EQ.12013) THEN ! Ground minimum temperature at 2m (12 bits)
+         ELSE IF (desc.EQ.12114.OR.  ! Maximum temperature at 2m, past 12 hours (16 bits)
-            IF (TGN_12.EQ.rvind) THEN
+     +           desc.EQ.12014) THEN ! Maximum temperature at 2m, past 12 hours (12 bits)
-               TGN_12 = value
-            END IF
-         ELSE IF (desc.EQ.12114) THEN ! Maximum temperature at 2m, past 12 hours (16 bits)
-            IF (TAX_12.EQ.rvind) THEN
-               TAX_12 = value
-            END IF
-         ELSE IF (desc.EQ.12014) THEN ! Maximum temperature at 2m, past 12 hours (12 bits)
             IF (TAX_12.EQ.rvind) THEN
                TAX_12 = value
@@ Line 1501: / Line 1808: @@
             END IF
 C     Do we also need to consider 12021 'Maximum temperature at 2m'?
-         ELSE IF (desc.EQ.12115) THEN ! Minimum temperature at 2m, past 12 hours (16 bits)
+         ELSE IF (desc.EQ.12115.OR.  ! Minimum temperature at 2m, past 12 hours (16 bits)
-            IF (TAN_12.EQ.rvind) THEN
+     +           desc.EQ.12015) THEN ! Minimum temperature at 2m, past 12 hours (12 bits)
-               TAN_12 = value
-            END IF
-         ELSE IF (desc.EQ.12015) THEN ! Minimum temperature at 2m, past 12 hours (12 bits)
             IF (TAN_12.EQ.rvind) THEN
                TAN_12 = value
@@ Line 1549: / Line 1853: @@
 C     number of cloud layers if previous descriptor is cloud
 C     type, according to all WMO recommended templates
+C     (but DNMI metar is an exception!)
             IF (ktdexp(idx - 1).EQ.20012) THEN
                IF (NINT(value).LE.4) THEN
@@ Line 1567: / Line 1872: @@
                cloud_layer = cloud_type_count - 2
                IF (cloud_layer.LE.num_cloud_layers) THEN
-                  IF (value.NE.rvind) THEN
+                  vert_sign(cloud_layer) = value
-                     vert_sign(cloud_layer) = value
-                  END IF
                END IF
             ELSE                ! rdb-files always have 0 or 4 cloud layers
                cloud_layer = cloud_type_count - 2
                IF (cloud_layer.LT.5) THEN
-                  IF (value.NE.rvind) THEN
+                  vert_sign(cloud_layer) = value
-                     vert_sign(cloud_layer) = value
-                  END IF
                END IF
             END IF
          ELSE IF (desc.EQ.20011 .AND. .NOT.bad_cloud_data) THEN ! Cloud amount
-            IF (cloud_type_count.EQ.0) THEN ! First occurrence
+            IF (metar) THEN
+               IF (num_cloud_layers.GT.-1) THEN
+                  num_cloud_layers = num_cloud_layers + 1
+               ELSE
+                  num_cloud_layers = 1
+               END IF
+               NS(num_cloud_layers) = value
+            ELSE IF (cloud_type_count.EQ.0) THEN ! First occurrence
                IF (NH.EQ.rvind) THEN
                   NH = value
@@ Line 1590: / Line 1898: @@
                cloud_layer = cloud_type_count - 2
                IF (cloud_layer.LE.num_cloud_layers) THEN
-                  IF (value.NE.rvind) THEN
+                  NS(cloud_layer) = value
-                     NS(cloud_layer) = value
-                  END IF
                END IF
             ELSE                ! rdb-files always have 0 or 4 cloud layers
                cloud_layer = cloud_type_count - 2
                IF (cloud_layer.LT.5) THEN
-                  IF (value.NE.rvind) THEN
+                  NS(cloud_layer) = value
-                     NS(cloud_layer) = value
-                  END IF
                END IF
             END IF
          ELSE IF (desc.EQ.20012 .AND. .NOT.bad_cloud_data) THEN ! Cloud type
-            cloud_type_count = cloud_type_count + 1
+            IF (metar) THEN
-            IF (cloud_type_count.GT.3) THEN
+               CC(num_cloud_layers) = value
-               cloud_layer = cloud_type_count - 3
+            ELSE
-               IF (num_cloud_layers .GT.-1) THEN
+               cloud_type_count = cloud_type_count + 1
-                  IF (cloud_layer.LE.num_cloud_layers) THEN
+               IF (cloud_type_count.GT.3) THEN
-                     IF (value.NE.rvind) THEN
+                  cloud_layer = cloud_type_count - 3
+                  IF (num_cloud_layers .GT.-1) THEN
+                     IF (value < 10.0 ! Accept one digit values only
+     +                    .AND. cloud_layer.LE.num_cloud_layers) THEN
                         CC(cloud_layer) = value
                      END IF
-                  END IF
+                  ELSE IF (cloud_layer.LT.5) THEN ! rdb-files always have 0 or 4 cloud layers
-               ELSE IF (cloud_layer.LT.5) THEN ! rdb-files always have 0 or 4 cloud layers
-                  IF (value.NE.rvind) THEN
                      CC(cloud_layer) = value
                   END IF
-               END IF
+               ELSE
-            ELSE
+                  IF (cloud_type_count.EQ.1) THEN
-               IF (cloud_type_count.EQ.1) THEN
+                     IF (CL.EQ.rvind) THEN
-                  IF (CL.EQ.rvind) THEN
+                        CL = value
-                     CL = value
+                     END IF
-                  END IF
+                  ELSE IF (cloud_type_count.EQ.2) THEN
-               ELSE IF (cloud_type_count.EQ.2) THEN
+                     IF (CM.EQ.rvind) THEN
-                  IF (CM.EQ.rvind) THEN
+                        CM = value
-                     CM = value
+                     END IF
-                  END IF
+                  ELSE IF (cloud_type_count.EQ.3) THEN
-               ELSE IF (cloud_type_count.EQ.3) THEN
+                     IF (CH.EQ.rvind) THEN
-                  IF (CH.EQ.rvind) THEN
+                        CH = value
-                     CH = value
+                     END IF
                   END IF
                END IF
             END IF
          ELSE IF (desc.EQ.20013 .AND. .NOT.bad_cloud_data) THEN ! Height of base of cloud
-            IF (cloud_type_count.EQ.0) THEN ! First occurrence
+            IF (metar) THEN
+               HS(num_cloud_layers) = value
+            ELSE IF (cloud_type_count.EQ.0) THEN ! First occurrence
                IF (HL.EQ.rvind) THEN
                   HL = value
@@ Line 1646: / Line 1953: @@
                cloud_layer = cloud_type_count - 3
                IF (cloud_layer.LE.num_cloud_layers) THEN
-                  IF (value.NE.rvind) THEN
+                  HS(cloud_layer) = value
-                     HS(cloud_layer) = value
-                  END IF
                END IF
             ELSE ! rdb-files always have 0 or 4 cloud layers
-               IF (value.NE.rvind) THEN
+               cloud_layer = cloud_type_count - 3
-                  cloud_layer = cloud_type_count - 3
+               IF (cloud_layer.LT.5) THEN
-                  IF (cloud_layer.LT.5) THEN
+                  HS(cloud_layer) = value
-                     IF (value.NE.rvind) THEN
-                        HS(cloud_layer) = value
-                     END IF
-                  END IF
                END IF
             END IF
@@ Line 1696: / Line 1997: @@
             END IF
          ELSE IF (desc.EQ.13013) THEN ! Total snow depth
-            IF (value.NE.rvind) THEN
 C     Don't check for SA.EQ.rvind, because SA might earlier have been set to 0 if
-C     E < 10, which probably means that 20062 has been wrongly encoded, not 13013
+C     EE < 10, which probably means that 20062 has been wrongly encoded, not 13013
-               SA = value
+            SA = value
+         ELSE IF (desc.EQ.20062) THEN ! State of the ground (with or without snow)
+            IF (EE.EQ.rvind) THEN
+               EE = value
             END IF
-         ELSE IF (desc.EQ.20062) THEN ! State of the ground (with or without snow)
+            IF (NINT(value).LE.10 .AND. SA.EQ.rvind) THEN
-            IF (value.NE.rvind) THEN
+               SA = 0
-               IF (NINT(value).LE.10) THEN ! E in TAC (3Ejjj)
+            END IF
-                  E = value
+C     Equating 013012 with SS_24 is dubious generally, but this is how
-                  IF (SA.EQ.rvind) THEN
+C     SS_24 is encoded in kvalobs for Norwegian avalanche stations (and
-                     SA = 0
+C     we have no better kvalobs parameter for depth of fresh snow anyway)
-                  END IF
+         ELSE IF (desc.EQ.13012) THEN ! Depth of fresh snow
-               ELSE IF (NINT(value).LE.20) THEN ! E' in TAC (4Esss)
+            IF (SS_24.EQ.rvind) THEN
-                  EM = value
+               SS_24 = value
-               END IF
             END IF
          ELSE IF (desc.EQ.14031) THEN ! Total sunshine
@@ Line 1792: / Line 2094: @@
             IF (hhh.EQ.rvind) THEN
                hhh = value
-            END IF
-         ELSE IF (desc.EQ.10008) THEN ! Geopotential (20 bits)
-            IF (hhh.EQ.rvind) THEN
-               hhh = value * 9.8
-            END IF
-         ELSE IF (desc.EQ.10003) THEN ! Geopotential (17 bits)
-            IF (hhh.EQ.rvind) THEN
-               hhh = value * 9.8
             END IF
 C     Special for ship or marine stations
          ELSE IF (desc.EQ.1011) THEN  ! Ship or mobile land station identifier
-            IF (value.NE.rvind) THEN
+            cidx = int(value/1000)
-               cidx = int(value/1000)
+            IF (cidx.GT.0) THEN
-               IF (cvals(cidx).NE.spc9) THEN
+               call_sign = cvals(cidx) ! CCITTIA5 data
-                  call_sign = cvals(cidx) ! CCITTIA5 data
+               call_sign = ctrim(call_sign,9,missing9)
-                  call_sign = ctrim(call_sign,9,missing9)
-               END IF
             END IF
          ELSE IF (desc.EQ.1012) THEN ! Direction of motion of moving observing platform
@@ Line 1819: / Line 2111: @@
             END IF
          ELSE IF (desc.EQ.7062) THEN ! Depth below sea/water surface
-            IF (value.NE.rvind) THEN
 C     Some buoy reports starts with depth 1.5 m, others starts with 0 m
 C     then 1.5 m and always have same sea/water temperature for these 2
 C     levels, so it seems like 0 m should be considered equivalent with 1.5 m
-               IF (value.LT.1.6) THEN
+            IF (value.LT.1.6) THEN
-                  surface_data = .TRUE.
+               surface_data = .TRUE.
-               ELSE
+            ELSE
-                  surface_data = .FALSE.
+               surface_data = .FALSE.
-               END IF
             END IF
-         ELSE IF (desc.EQ.22043) THEN ! Sea/water temperature (15 bits)
+         ELSE IF (desc.EQ.22043.OR.  ! Sea/water temperature (15 bits)
+     +           desc.EQ.22042.OR.   ! Sea/water temperature (12 bits)
+     +           desc.EQ.22049) THEN ! Sea-surface temperature (15 bits)
             IF (TW.EQ.rvind .AND. surface_data)THEN
                TW = value
             END IF
-         ELSE IF (desc.EQ.22042) THEN ! Sea/water temperature (12 bits)
+         ELSE IF (desc.EQ.12102.OR.  ! Wet-bulb temperature (16 bits)
-            IF (TW.EQ.rvind
+     +           desc.EQ.12005) THEN ! Wet-bulb temperature (12 bits)
-     +           .AND. surface_data) THEN
-               TW = value
-            END IF
-         ELSE IF (desc.EQ.12102) THEN ! Wet-bulb temperature (16 bits)
             IF (TbTbTb.EQ.rvind) THEN
                TbTbTb = value
             END IF
-         ELSE IF (desc.EQ.12005) THEN ! Wet-bulb temperature (12 bits)
+         ELSE IF (desc.EQ.22011) THEN !  Period of waves (instrumentally measured)
-            IF (TbTbTb.EQ.rvind) THEN
+            IF (PWA.EQ.rvind) THEN
-               TbTbTb = value
+               PWA = value
             END IF
-         ELSE IF (desc.EQ.22011) THEN !  Period of waves (BUFR doesn't distinguish between PW and PWA)
+         ELSE IF (desc.EQ.22012) THEN !  Period of wind waves (visually measured)
             IF (PW.EQ.rvind) THEN
                PW = value
             END IF
-         ELSE IF (desc.EQ.22021) THEN ! Heigth of waves (BUFR doesn't distinguish between HW and HWA)
+         ELSE IF (desc.EQ.22021) THEN ! Heigth of waves
+            IF (HWA.EQ.rvind) THEN
+               HWA = value
+            END IF
+         ELSE IF (desc.EQ.22022) THEN ! Heigth of wind waves
             IF (HW.EQ.rvind) THEN
                HW = value
             END IF
          ELSE IF (desc.EQ.22003) THEN ! Direction of swell waves
-            IF (value.NE.rvind) THEN
+            IF (DW1.EQ.rvind) THEN
-               IF (DW1.EQ.rvind) THEN
+               DW1 = value
-                  DW1 = value
+            ELSE
-               ELSE
+               DW2 = value
-                  DW2 = value
-               END IF
             END IF
          ELSE IF (desc.EQ.22013) THEN ! Period of swell waves
-            IF (value.NE.rvind) THEN
+            IF (PW1.EQ.rvind) THEN
-               IF (PW1.EQ.rvind) THEN
+               PW1 = value
-                  PW1 = value
+            ELSE
-               ELSE
+               PW2 = value
-                  PW2 = value
-               END IF
             END IF
          ELSE IF (desc.EQ.22023) THEN ! Height of swell waves
-            IF (value.NE.rvind) THEN
+            IF (HW1.EQ.rvind) THEN
-               IF (HW1.EQ.rvind) THEN
+               HW1 = value
-                  HW1 = value
+            ELSE
-               ELSE
+               HW2 = value
-                  HW2 = value
-               END IF
             END IF
          ELSE IF (desc.EQ.22061) THEN ! State of the sea
@@ Line 1915: / Line 2201: @@
             END IF
          ELSE IF (desc.EQ.1005) THEN ! Buoy/platform identifier
-            IF (buoy_id.EQ.rvind) THEN
+            IF (buoy_id5.EQ.rvind) THEN
-               buoy_id = value
+               buoy_id5 = value
             END IF
          ELSE IF (desc.EQ.1003) THEN ! WMO region number/geographical area
@@ Line 1925: / Line 2211: @@
             IF (wmo_region_subarea.EQ.rvind) THEN
                wmo_region_subarea = value
+            END IF
+         ELSE IF (desc.EQ.1087) THEN ! WMO Marine observing platform extended identifier
+            IF (buoy_id7.EQ.rvind) THEN
+               buoy_id7 = value
             END IF
 C     Special for metar
          ELSE IF (desc.EQ.1063) THEN  ! ICAO location indicator
-            IF (value.NE.rvind) THEN
+            cidx = int(value/1000)
-               cidx = int(value/1000)
+            icao_id = cvals(cidx) ! CCITTIA5 data
-               IF (cvals(cidx).NE.spc8) THEN
+            icao_id = ctrim(icao_id,8,missing8)
-                  icao_id = cvals(cidx) ! CCITTIA5 data
-                  icao_id = ctrim(icao_id,8,missing8)
-               END IF
-            END IF
          ELSE IF (desc.EQ.10052) THEN ! Altimeter setting (QNH)
             IF (PH.EQ.rvind) THEN
@@ Line 1948: / Line 2234: @@
       END IF
-      WRITE(*,*)
       IF (II.NE.rvind .AND. iii.NE.rvind) THEN
+         WRITE(*,*)
          WRITE(*,'(A,I5.5)') 'wmonr=',NINT(II)*1000 + NINT(iii)
+      ELSE IF (wigos_series.NE.rvind .AND. wigos_issuer.NE.missing5
+     +        .AND. wigos_issueno.NE.missing5
+     +        .AND. wigos_localid.NE.missing16) THEN
+         ind = index(wigos_issuer,' ') - 1
+         IF (ind.EQ.-1) ind = 5
+         ind2 = index(wigos_issueno,' ') - 1
+         IF (ind2.EQ.-1) ind2 = 5
+         ind3 = index(wigos_localid,' ') - 1
+         IF (ind3.EQ.-1) ind3 = 16
+         WRITE(*,*)
+         WRITE(*,'(A,I1.1,A1,A,A1,A,A1,A)')
+     +        'wigosid=',NINT(wigos_series),
+     +        '-',wigos_issuer(1:ind),
+     +        '-',wigos_issueno(1:ind2),
+     +        '-',wigos_localid(1:ind3)
       ELSE IF (state_id.NE.rvind .AND. national_number.NE.rvind) THEN
+         WRITE(*,*)
          WRITE(*,'(A,I3.3,A1,I10.10)') 'nationalnr=',NINT(state_id),
      +        '_',NINT(national_number)
@@ Line 1962: / Line 2264: @@
             ind = ind - 1
          END DO
+         WRITE(*,*)
          WRITE(*,'(A,A)') 'call_sign=',
      +        call_sign(1:ind)
-      ELSE IF (buoy_id.NE.rvind.AND.wmo_region_number.NE.rvind
+      ELSE IF (buoy_id7.NE.rvind) THEN
-     +        .AND.wmo_region_subarea.NE.rvind) THEN
+C     New templates introduced in 2014 for data category 1 use 001087
-         WRITE(*,'(A,I5)') 'buoy_id=',NINT(wmo_region_number)*10000
+C     WMO Marine observing platform extended identifier, 7 digits
-     +        + NINT(wmo_region_subarea)*1000 + NINT(buoy_id)
+         WRITE(*,*)
-      ELSE IF (buoy_id.NE.rvind.AND.buoy_id.GT.1000) THEN
+         WRITE(*,'(A,I7)') 'buoy_id=',NINT(buoy_id7)
+      ELSE IF (buoy_id5.NE.rvind) THEN
+         WRITE(*,*)
+         IF (wmo_region_number.EQ.rvind
+     +        .AND.wmo_region_subarea.EQ.rvind) THEN
 C     Old drau files (wrongly) includes wmo_region_number and
-C     wmo_region_subarea in buoy_id
+C     wmo_region_subarea in 001005 Buoy/platform identifier. Should we
-         WRITE(*,'(A,I5)') 'buoy_id=',NINT(buoy_id)
+C     expand this to 7 digits by inserting '00'?
+            WRITE(*,'(A,I5)') 'buoy_id=',NINT(buoy_id5)
+         ELSE IF (wmo_region_number.EQ.rvind
+     +        .AND.wmo_region_subarea.NE.rvind) THEN
+C     Some BUFR BUOYS on GTS have 'missing' value for wmo_region_number,
+C     but not for wmo_region_subarea
+            WRITE(*,'(A,I6)') 'buoy_id=',
+     +           NINT(wmo_region_subarea)*100000 + NINT(buoy_id5)
+         ELSE IF (wmo_region_number.NE.rvind
+     +        .AND.wmo_region_subarea.EQ.rvind) THEN
+C     Not easy to know how to display this case, but then I have never
+C     seen this in practice
+            WRITE(*,'(A,I5)') 'buoy_id=',NINT(buoy_id5)
+         ELSE IF (buoy_id5.GT.10000.AND.
+     +           (NINT(buoy_id5) - MOD(NINT(buoy_id5),1000)).EQ.
+     +           NINT(wmo_region_number)*10000
+     +           + NINT(wmo_region_subarea)*1000) THEN
+C     If first 2 digits of 5 digit buoy_id equals wmo_region_number and
+C     wmo_region_subarea respectively, this is almost certainly an
+C     encoding error and we choose to show last 3 digits of buoy_id5 only
+           WRITE(*,'(A,I7)') 'buoy_id=',NINT(wmo_region_number)*1000000
+     +           + NINT(wmo_region_subarea)*100000
+     +           + MOD(NINT(buoy_id5),1000)
+         ELSE
+           WRITE(*,'(A,I7)') 'buoy_id=',NINT(wmo_region_number)*1000000
+     +           + NINT(wmo_region_subarea)*100000 + NINT(buoy_id5)
+         END IF
       ELSE
          IF (verbose .GT. 1) THEN
-            WRITE(*,*) 'Both wmonr, nationalnr, call_sign and',
+            WRITE(*,*)
-     +           ' buoy_id are missing!!!'
+            WRITE(*,*) 'Both wmonr, wigosid, nationalnr, call_sign',
+     +           ' and buoy_id are missing!!!'
 C     Example: ISRZ47 EGRR has no proper station identification (except
 C     station name and position)
@@ Line 2002: / Line 2336: @@
          WRITE(*,'(A,A)') 'name=',long_name(1:lenstr(long_name,1))
       END IF
-      IF (NINT(ix).EQ.0) THEN
+      IF (ix.NE.rvind) THEN
-         WRITE(*,'(A,A)') 'type=Automatic'
+         IF (NINT(ix).EQ.0) THEN
-      ELSE IF (NINT(ix).EQ.1) THEN
+            WRITE(*,'(A,A)') 'type=Automatic'
-         WRITE(*,'(A,A)') 'type=Manned'
+         ELSE IF (NINT(ix).EQ.1) THEN
-      ELSE IF (NINT(ix).EQ.2) THEN
+            WRITE(*,'(A,A)') 'type=Manned'
-         WRITE(*,'(A,A)') 'type=Hybrid'
+         ELSE IF (NINT(ix).EQ.2) THEN
+            WRITE(*,'(A,A)') 'type=Hybrid'
+         END IF
       END IF
       IF (longitude.NE.rvind) THEN
@@ Line 2018: / Line 2354: @@
                                 ! so display as integer
          WRITE(*,'(A,I7)') 'height=',NINT(height)
+      END IF
+      IF (hp.NE.rvind) THEN
+         WRITE(*,'(A,F11.1)') 'hp=',hp
       END IF
       IF (vs.NE.rvind) THEN
@@ Line 2024: / Line 2363: @@
       IF (ds.NE.rvind) THEN
          WRITE(*,'(A,I11)') 'ds=',NINT(ds)
+      END IF
+      IF (a3.NE.rvind) THEN     ! Standard isobaric surface for which the
+                                ! geopotential is reported, no Kvalobs code exists
+         WRITE(*,'(A,I11)') 'a3=',NINT(a3/100) ! hPa
       END IF
       IF (hhh.NE.rvind) THEN    ! Geopotential height, no Kvalobs code exists
@@ Line 2112: / Line 2455: @@
          WRITE(*,'(A,F11.1)') 'PP=',PP/100 ! hPa
       END IF
-C     Precipitation = -0.1 means "trace" (less than 0.05 kg/m2), is
+C     Use the BUFR some special values for precipitation, *not* the
-C     represented as 0 in Kvalobs
+C     Kvalobs special values. So 'trace' (less than 0.05 kg/m2), is
+C     showed as -0.1 (not 0 as in Kvalobs)
       IF (RR_24.NE.rvind) THEN
-         IF (NINT(RR_24*10).EQ.-1) THEN
+         WRITE(*,'(A,F8.1)') 'RR_24=',RR_24
-            WRITE(*,'(A,F8.1)') 'RR_24=',0.0
-         ELSE
-            WRITE(*,'(A,F8.1)') 'RR_24=',RR_24
-         END IF
       END IF
       IF (RR_12.NE.rvind) THEN
-         IF (NINT(RR_12*10).EQ.-1) THEN
+         WRITE(*,'(A,F8.1)') 'RR_12=',RR_12
-            WRITE(*,'(A,F8.1)') 'RR_12=',0.0
-         ELSE
-            WRITE(*,'(A,F8.1)') 'RR_12=',RR_12
-         END IF
       END IF
       IF (RR_6.NE.rvind) THEN
-         IF (NINT(RR_6*10).EQ.-1) THEN
+         WRITE(*,'(A,F9.1)') 'RR_6=',RR_6
-            WRITE(*,'(A,F9.1)') 'RR_6=',0.0
-         ELSE
-            WRITE(*,'(A,F9.1)') 'RR_6=',RR_6
-         END IF
       END IF
       IF (RR_3.NE.rvind) THEN
-         IF (NINT(RR_3*10).EQ.-1) THEN
+         WRITE(*,'(A,F9.1)') 'RR_3=',RR_3
-            WRITE(*,'(A,F9.1)') 'RR_3=',0.0
-         ELSE
-            WRITE(*,'(A,F9.1)') 'RR_3=',RR_3
-         END IF
       END IF
       IF (RR_1.NE.rvind) THEN
-         IF (NINT(RR_1*10).EQ.-1) THEN
+         WRITE(*,'(A,F9.1)') 'RR_1=',RR_1
-            WRITE(*,'(A,F9.1)') 'RR_1=',0.0
-         ELSE
-            WRITE(*,'(A,F9.1)') 'RR_1=',RR_1
-         END IF
       END IF
       IF (WW.NE.rvind.AND.WW.LT.200) THEN ! 508-511 and w1w1 (in 333 9 group) ignored here
@@ Line 2170: / Line 2494: @@
          END IF
       END IF
-      IF (E.NE.rvind) THEN
+      IF (EE.NE.rvind) THEN
-         WRITE(*,'(A,I12)') 'E=',NINT(E)
+         WRITE(*,'(A,I11)') 'EE=',NINT(EE)
       END IF
-      IF (EM.NE.rvind) THEN
+C     Use the BUFR some special values for SA, *not* the Kvalobs special values
-         WRITE(*,'(A,I11)') 'EM=',NINT(EM)-10
+C     So 'trace' is showed as SA=-1, 'Snow cover not continuous' as SA=-2
-      END IF
-C     SA has some special values in BUFR as well as in Kvalobs
 C     Note that conversion from synop to BUFR normally will set SA=0 if E < 10
-c$$$      IF (E.NE.rvind .AND. NINT(E).LE.10
-c$$$     +     .AND. (SA.EQ.rvind .OR. NINT(SA).EQ.0) ) THEN
-c$$$         WRITE(*,'(A,I11)') 'SA=',-1 ! No snow
       IF (SA.NE.rvind) THEN
-         IF (NINT(SA*100).EQ.-1) THEN ! Trace: less than 0.5 cm snow
+         WRITE(*,'(A,I11)') 'SA=',NINT(SA*100)
-            WRITE(*,'(A,I11)') 'SA=',0
+      END IF
-         ELSE IF (NINT(SA*100).EQ.-2) THEN ! Snow cover not continuos
+      IF (SS_24.NE.rvind) THEN
-            WRITE(*,'(A,I11)') 'SA=',-1
+         WRITE(*,'(A,I8)') 'SS_24=',NINT(SS_24*100)
-         ELSE IF (NINT(SA*100).EQ.0) THEN ! 0 snow coded as -1 in Kvalobs. Stupid but true
-            WRITE(*,'(A,I11)') 'SA=',-1
-         ELSE
-            WRITE(*,'(A,I11)') 'SA=',NINT(SA*100)
-         END IF
       END IF
       IF (VV.NE.rvind) THEN
@@ Line 2224: / Line 2538: @@
       DO i=1,num_cloud_layers
          IF (NS(i).NE.rvind) THEN
-            WRITE(*,'(A,I1,A,I10)') 'NS',i,'=',NNtoWMO_N(NINT(NS(i)))
+            WRITE(*,'(A,I1,A,I10)') 'NS',i,'=',NINT(NS(i))
          END IF
          IF (CC(i).NE.rvind) THEN
@@ Line 2230: / Line 2544: @@
          END IF
          IF (HS(i).NE.rvind) THEN
-            WRITE(*,'(A,I1,A,I10)') 'HS',i,'=',HStoWMO_HSHS(HS(i))
+            WRITE(*,'(A,I1,A,I10)') 'HS',i,'=',NINT(HS(i))
          END IF
       END DO
       IF (SG.NE.rvind) THEN
          WRITE(*,'(A,I11)') 'SG=',NINT(SG)
+      END IF
+      IF (PWA.NE.rvind) THEN
+         WRITE(*,'(A,I10)') 'PWA=',NINT(PWA)
       END IF
       IF (PW.NE.rvind) THEN
          WRITE(*,'(A,I11)') 'PW=',NINT(PW)
+      END IF
+      IF (HWA.NE.rvind) THEN
+         WRITE(*,'(A,F10.1)') 'HWA=',HWA
       END IF
       IF (HW.NE.rvind) THEN
@@ Line 2366: / Line 2686: @@
       REAL*8 value
       REAL*8 d(maxlevel),s(maxlevel),T(maxlevel),v(maxlevel),z(maxlevel)
+      REAL*8 w(maxlevel)
 C     Variables used for geographical filtering av observations
@@ Line 2397: / Line 2718: @@
       DO n=1,maxlevel
          z(n) = rvind
+         w(n) = rvind
          T(n) = rvind
          s(n) = rvind
@@ Line 2435: / Line 2757: @@
                minute = value
             END IF
-         ELSE IF (desc.EQ.5001) THEN ! Latitude (high accuracy)
+         ELSE IF (desc.EQ.5001.OR.  ! Latitude (high accuracy)
+     +           desc.EQ.5002) THEN ! Latitude (coarse accuracy)
             IF (latitude.EQ.rvind) THEN
                latitude = value
             END IF
-         ELSE IF (desc.EQ.5002) THEN ! Latitude (coarse accuracy)
+         ELSE IF (desc.EQ.6001.OR.  ! Longitude (high accuracy)
-            IF (latitude.EQ.rvind) THEN
+     +           desc.EQ.6002) THEN ! Longitude (coarse accuracy)
-               latitude = value
-            END IF
-         ELSE IF (desc.EQ.6001) THEN ! Longitude (high accuracy)
-            IF (longitude.EQ.rvind) THEN
-               longitude = value
-            END IF
-         ELSE IF (desc.EQ.6002) THEN ! Longitude (coarse accuracy)
             IF (longitude.EQ.rvind) THEN
                longitude = value
@@ Line 2462: / Line 2778: @@
             n = n + 1           ! new level
             IF (n.GT.maxlevel) THEN
+               n = maxlevel
                WRITE(*,*) 'Too many levels! Skipping rest of message'
                GOTO 120
             END IF
             z(n) = value
-         ELSE IF (desc.EQ.22043) THEN ! Sea/water temperature (15 bits)
+         ELSE IF (desc.EQ.7065) THEN ! Water pressure
-            IF (n.GT.0) THEN
+            n = n + 1           ! new level
-               t(n) = value
+            IF (n.GT.maxlevel) THEN
-            ELSE IF (TW.EQ.rvind) THEN
+               n = maxlevel
-               TW = value
+               WRITE(*,*) 'Too many levels! Skipping rest of message'
+               GOTO 120
             END IF
-         ELSE IF (desc.EQ.22042) THEN ! Sea/water temperature (12 bits)
+            w(n) = value
-            IF (n.GT.0) THEN
+         ELSE IF (desc.EQ.22045.OR.  ! Sea/water temperature (19 bits)
+     +           desc.EQ.22043.OR.   ! Sea/water temperature (15 bits)
+     +           desc.EQ.22042) THEN ! Sea/water temperature (12 bits)
+            IF (n.GT.0 .AND. n.LE.maxlevel) THEN
                t(n) = value
-            ELSE IF (TW.EQ.rvind) THEN
+            ELSE IF (n.EQ.0 .AND. TW.EQ.rvind) THEN
                TW = value
             END IF
-         ELSE IF (desc.EQ.22062) THEN ! Salinity [part per thousand]
+         ELSE IF (desc.EQ.22064.OR.  ! Salinity [part per thousand] (17 bits)
-            IF (n.GT.0) THEN
+     +           desc.EQ.22062) THEN ! Salinity [part per thousand] (14 bits)
+            IF (n.GT.0 .AND. n.LE.maxlevel) THEN
                s(n) = value
-            ELSE IF (sal.EQ.rvind) THEN
+            ELSE IF (n.EQ.0 .AND. sal.EQ.rvind) THEN
                sal = value
             END IF
          ELSE IF (desc.EQ.22004) THEN ! Direction of current
-               IF (n.GT.0) THEN
+               IF (n.GT.0 .AND. n.LE.maxlevel) THEN
                   d(n) = value
-               ELSE IF (verbose.GT.1) THEN
+               ELSE IF (n.EQ.0 .AND. verbose.GT.1) THEN
                   WRITE(*,*) 'Find 22004 (dc) before first 7062 (zz)!'
                END IF
          ELSE IF (desc.EQ.22031) THEN ! Speed of current
-               IF (n.GT.0) THEN
+               IF (n.GT.0 .AND. n.LE.maxlevel) THEN
                   v(n) = value
-               ELSE IF (verbose.GT.1) THEN
+               ELSE IF (n.EQ.0 .AND. verbose.GT.1) THEN
                   WRITE(*,*) 'Find 22031 (vc) before first 7062 (zz)!'
                END IF
          ELSE IF (desc.EQ.22062) THEN ! Salinity [part per thousand]
-            IF (n.GT.0) THEN
+            IF (n.GT.0 .AND. n.LE.maxlevel) THEN
                s(n) = value
-            ELSE IF (sal.EQ.rvind) THEN
+            ELSE IF (n.EQ.0 .AND. sal.EQ.rvind) THEN
                sal = value
             END IF
@@ Line 2529: / Line 2851: @@
             END IF
          ELSE IF (desc.EQ.1087) THEN ! WMO marine observing platform extended identifier
-            IF (verbose .GT. 1) THEN
+            IF (buoy_id.EQ.rvind) THEN
-               WRITE(*,*) '001087 as stationid - unable to decode'
+               buoy_id = value
             END IF
-            RETURN
          END IF
       END DO
@@ Line 2544: / Line 2865: @@
       END IF
-      WRITE(*,*)
       IF (buoy_id.NE.rvind.AND.wmo_region_number.NE.rvind
      +     .AND.wmo_region_subarea.NE.rvind) THEN
-         WRITE(*,'(A,I5)') 'buoy_id=',NINT(wmo_region_number)*10000
+         WRITE(*,*)
-     +        + NINT(wmo_region_subarea)*1000 + NINT(buoy_id)
+         IF (buoy_id.LT.1000) THEN
+            WRITE(*,'(A,I5)') 'buoy_id=',NINT(wmo_region_number)*10000
+     +           + NINT(wmo_region_subarea)*1000 + NINT(buoy_id)
+         ELSE
+C     This is an error in encoding. We choose to show buoy_id only (and
+C     in all cases I have seen of this error, first 2 digits of buoy_id
+C     is wmo_region_number and wmo_region_subarea respectively).
+            IF (buoy_id.LT.10000) THEN
+               WRITE(*,'(A,I4)') 'buoy_id=',NINT(buoy_id)
+            ELSE IF (buoy_id.LT.100000) THEN
+               WRITE(*,'(A,I5)') 'buoy_id=',NINT(buoy_id)
+            ELSE IF (buoy_id.LT.1000000) THEN
+               WRITE(*,'(A,I6)') 'buoy_id=',NINT(buoy_id)
+            ELSE
+               WRITE(*,'(A,I7)') 'buoy_id=',NINT(buoy_id)
+            END IF
+         END IF
       ELSE IF (buoy_id.NE.rvind.AND.buoy_id.GT.1000) THEN
 C     Old drau files (wrongly) includes wmo_region_number and
 C     wmo_region_subarea in buoy_id
-         WRITE(*,'(A,I5)') 'buoy_id=',NINT(buoy_id)
+         WRITE(*,*)
+         IF (buoy_id.LT.10000) THEN ! 001005
+            WRITE(*,'(A,I5)') 'buoy_id=',NINT(buoy_id)
+         ELSE                       ! 001087
+            WRITE(*,'(A,I7)') 'buoy_id=',NINT(buoy_id)
+         END IF
       ELSE IF (call_sign.NE.missing9) THEN
          ind = index(call_sign,' ') - 1
@@ Line 2561: / Line 2902: @@
             ind = ind - 1
          END DO
+         WRITE(*,*)
          WRITE(*,'(A,A)') 'call_sign=',
      +        call_sign(1:ind)
       ELSE
          IF (verbose .GT. 1) THEN
+            WRITE(*,*)
             WRITE(*,*) 'Both buoy_id and call_sign are missing!!!'
          END IF
@@ Line 2606: / Line 2949: @@
 C     not quite uncommon)
          IF (z(n).EQ.rvind .AND. t(n).EQ.rvind .AND. s(n).EQ.rvind
-     +         .AND. d(n).EQ.rvind .AND. v(n).EQ.rvind) THEN
+     +        .AND. d(n).EQ.rvind .AND. v(n).EQ.rvind
+     +        .AND. w(n).EQ.rvind) THEN
             CYCLE
          END IF
@@ Line 2612: / Line 2956: @@
          IF (z(n).NE.rvind) THEN
             WRITE(*,'(A,I11)') 'zz=',NINT(z(n))
+         END IF
+         IF (w(n).NE.rvind) THEN
+            WRITE(*,'(A,I11)') 'wp=',NINT(w(n))
          END IF
          IF (t(n).NE.rvind) THEN
@@ Line 2648: / Line 2995: @@
       PARAMETER (rvind=1.7E38)
-      CHARACTER*8 aircraft,flight_number,missing8,spc8
+      CHARACTER*8 aircraft,flight_number,missing8
 C     Parameters defined in Kvalobs
       REAL*8 DD,FF,TT,PP,
 C     Other parameters
      +     year,month,day,hour,minute,second,latitude,longitude,
-     +     flight_level,phase,osn
+     +     flight_level,phase,osn,mixing_ratio
       CHARACTER*3 cphase,missing3
       INTEGER idx,cidx,ind
@@ Line 2674: / Line 3021: @@
       WRITE(missing8,'(8A)') one_bits,one_bits,one_bits,one_bits,
      +     one_bits,one_bits,one_bits,one_bits
-      spc8 = '        '
 C     Initialize all parameters to missing values
@@ Line 2695: / Line 3041: @@
       flight_level = rvind
       phase = rvind
+      mixing_ratio = rvind
 C     Loop through all expanded descriptors
@@ Line 2709: / Line 3056: @@
          IF (desc.EQ.1008) THEN ! Aircraft registration number or other identification
             cidx = int(value/1000)
-            IF (cvals(cidx).NE.spc8) THEN
+            aircraft = cvals(cidx) ! CCITTIA5 data
-               aircraft = cvals(cidx) ! CCITTIA5 data
+            aircraft = ctrim(aircraft,8,missing8)
-               aircraft = ctrim(aircraft,8,missing8)
-            END IF
          ELSE IF (desc.EQ.1006) THEN ! Aircraft flight number
             cidx = int(value/1000)
-            IF (cvals(cidx).NE.spc8) THEN
+            flight_number = cvals(cidx) ! CCITTIA5 data
-               flight_number = cvals(cidx) ! CCITTIA5 data
+            flight_number = ctrim(flight_number,8,missing8)
-               flight_number = ctrim(flight_number,8,missing8)
-            END IF
          ELSE IF (desc.EQ.1023) THEN ! Observation sequence number
             IF (osn.EQ.rvind) osn = value
@@ Line 2745: / Line 3088: @@
                second = value
             END IF
-         ELSE IF (desc.EQ.5001) THEN ! Latitude (high accuracy)
+         ELSE IF (desc.EQ.5001.OR.  ! Latitude (high accuracy)
+     +           desc.EQ.5002) THEN ! Latitude (coarse accuracy)
             IF (latitude.EQ.rvind) THEN
                latitude = value
             END IF
-         ELSE IF (desc.EQ.5002) THEN ! Latitude (coarse accuracy)
+         ELSE IF (desc.EQ.6001.OR.  ! Longitude (high accuracy)
-            IF (latitude.EQ.rvind) THEN
+     +           desc.EQ.6002) THEN ! Longitude (coarse accuracy)
-               latitude = value
-            END IF
-         ELSE IF (desc.EQ.6001) THEN ! Longitude (high accuracy)
             IF (longitude.EQ.rvind) THEN
                longitude = value
             END IF
-         ELSE IF (desc.EQ.6002) THEN ! Longitude (coarse accuracy)
+         ELSE IF (desc.EQ.7010.OR.   ! Flight level
-            IF (longitude.EQ.rvind) THEN
+     +           desc.EQ.7002.OR.    ! Height or altitude
-               longitude = value
+     +           desc.EQ.10070) THEN ! Indicated aircraft altitude
-            END IF
-         ELSE IF (desc.EQ.7010) THEN ! Flight level
             IF (flight_level.EQ.rvind) THEN
                flight_level = value
             END IF
-         ELSE IF (desc.EQ.7002) THEN ! Height or altitude
+         ELSE IF (desc.EQ.8009.OR.  ! Detailed phase of aircraft flight
-            IF (flight_level.EQ.rvind) THEN
+     +           desc.EQ.8004) THEN ! Phase of aircraft flight
-               flight_level = value
-            END IF
-         ELSE IF (desc.EQ.10070) THEN ! Indicated aircraft altitude
-            IF (flight_level.EQ.rvind) THEN
-               flight_level = value
-            END IF
-         ELSE IF (desc.EQ.8009) THEN ! Detailed phase of aircraft flight
             IF (phase.EQ.rvind) THEN
                phase = value
@@ Line 2783: / Line 3115: @@
                cphase = phase_8004(NINT(phase),missing3)
             END IF
-         ELSE IF (desc.EQ.12104) THEN ! Dry bulb temperature at 2m (data width 16 bits)
+         ELSE IF (desc.EQ.12104.OR. ! Dry bulb temperature at 2m (data width 16 bits)
+     +          desc.EQ.12004.OR.   ! Dry bulb temperature at 2m (12 bits)
+     +          desc.EQ.12101.OR.   ! Temperature/dry bulb temperature (16 bits)
+     +          desc.EQ.12001) THEN ! Temperature/dry bulb temperature (12 bits)
             IF (TT.EQ.rvind) THEN
                TT = value
             END IF
-         ELSE IF (desc.EQ.12004) THEN ! Dry bulb temperature at 2m (12 bits)
+         ELSE IF (desc.EQ.11011.OR.  ! Wind direction at 10 m
-            IF (TT.EQ.rvind) THEN
+     +           desc.EQ.11001) THEN ! Wind direction
-               TT = value
-            END IF
-         ELSE IF (desc.EQ.12101) THEN ! Temperature/dry bulb temperature (16 bits)
-            IF (TT.EQ.rvind) THEN
-               TT = value
-            END IF
-         ELSE IF (desc.EQ.12001) THEN ! Temperature/dry bulb temperature (12 bits)
-            IF (TT.EQ.rvind) THEN
-               TT = value
-            END IF
-         ELSE IF (desc.EQ.11011) THEN ! Wind direction at 10 m
             IF (DD.EQ.rvind) THEN
                DD = value
             END IF
-         ELSE IF (desc.EQ.11001) THEN ! Wind direction
+         ELSE IF (desc.EQ.11012.OR.  ! Wind speed at 10 m
-            IF (DD.EQ.rvind) THEN
+     +           desc.EQ.11002) THEN ! Wind speed
-               DD = value
-            END IF
-         ELSE IF (desc.EQ.11012) THEN ! Wind speed at 10 m
-            IF (FF.EQ.rvind) THEN
-               FF = value
-            END IF
-         ELSE IF (desc.EQ.11002) THEN ! Wind speed
             IF (FF.EQ.rvind) THEN
                FF = value
@@ Line 2818: / Line 3135: @@
             IF (PP.EQ.rvind) THEN
                PP = value
+            END IF
+         ELSE IF (desc.EQ.13002) THEN ! Mixing ratio
+            IF (mixing_ratio.EQ.rvind) THEN
+               mixing_ratio = value
             END IF
          END IF
@@ Line 2828: / Line 3149: @@
       END IF
-      WRITE(*,*)
       IF (aircraft.NE.missing8) THEN
+         WRITE(*,*)
          WRITE(*,'(A,A)') 'aircraft=',aircraft(1:lenstr(aircraft,1))
       ELSE IF (flight_number.NE.missing8) THEN
+         WRITE(*,*)
          WRITE(*,'(A,A)') 'aircraft=',
      +        flight_number(1:lenstr(flight_number,1))
       ELSE
          IF (verbose .GT. 1) THEN
+            WRITE(*,*)
             WRITE(*,*) 'Aircraft (001008/001006) is missing!!!'
          END IF
@@ Line 2882: / Line 3205: @@
       END IF
       IF (PP.NE.rvind) THEN
-         WRITE(*,'(A,I16)') 'PP=',NINT(PP)
+         WRITE(*,'(A,F16.1)') 'PP=',PP/100 ! hPa
+      END IF
+      IF (mixing_ratio.NE.rvind) THEN
+         WRITE(*,'(A,F16.2)') 'mr=',mixing_ratio*1000 ! g/kg
       END IF
@@ Line 2953: / Line 3279: @@
          phase_8009 = missing3
       ELSE
-         phase_8009 = indicator(phase)
+         phase_8009 = indicator(phase + 1)
       END IF
       RETURN
       END FUNCTION phase_8009
-C     -----------------------------------------------------------------
-C     Convert value of HS (in meter) into WMO code for hshs (table 1677)
-      INTEGER FUNCTION HStoWMO_HSHS(HS)
-      IMPLICIT NONE
-      REAL*8 HS
-      INTEGER HSHS
-      IF (NINT(HS).LT.1800) THEN
-         HSHS=NINT(HS/30)-1
-      ELSE IF (NINT(HS).LT.10500) THEN
-         HSHS=NINT((HS-1800)/300)+55
-      ELSE IF (NINT(HS).LE.21000) THEN
-         HSHS=NINT((HS-10500)/1500)+80
-      ELSE
-         HSHS=89
-      END IF
-      HStoWMO_HSHS = HSHS
-      RETURN
-      END FUNCTION HStoWMO_HSHS
 C     -----------------------------------------------------------------
@@ Line 3176: / Line 3480: @@
 c
       character*(*) text
-      integer       minlen
+      integer       k,l,lt,minlen
 c
       lt=len(text)
@@ Line 3190: / Line 3494: @@
       return
       end
 </code>
 <code fortran comfilter.f>
 C (C) Copyright 2010, met.no