Differences

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

--- bufr.pm:bufrdump [2011-09-15 09:20:45]
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 26: / 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 50: / 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 162: / 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 209: / Line 209: @@
       kerr = 0
+      metar = .FALSE.
 C     Do a partial expansion without quality control
       kreq(1) = 1               ! All original elements without quality control
@@ Line 227: / 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 239: / 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 254: / 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 294: / 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 306: / 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 331: / 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 343: / 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 379: / 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 400: / 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 '
      +     // 'exit status 1\n'
      +     // '\t--rectangle x1 y1 x2 y2\n'
      +     // '\t\t\tDecode observations inside rectangle with lower '
@@ Line 441: / 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 533: / Line 553: @@
       END IF
       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 546: / Line 567: @@
       SUBROUTINE err_msg(msg)
       IMPLICIT NONE
-      CHARACTER*80 msg
+      CHARACTER(LEN=*) msg
       WRITE(*,*) msg
       STOP
@@ Line 555: / Line 576: @@
       SUBROUTINE err_msg1(msg,kerr)
       IMPLICIT NONE
-      CHARACTER*80 msg
+      CHARACTER(LEN=*) msg
       INTEGER kerr
       WRITE(*,*) msg,kerr
@@ Line 657: / Line 678: @@
       CLOSE(11)
+      IF (verbose.GT.2) WRITE(*,*) 'nfidlines=',nfidlines
       IF (nfivlines.EQ.0) filter = .FALSE.
@@ Line 729: / Line 751: @@
       num_found = 0
-      ifiv = 0
 C     loop through all subsets:
       DO nsub = 1,num_subsets
 C     loop through all different conditions:
+         ifiv = 0
          DO i1 = 1,nfidlines
-C     loop through all filter values lines (for given) condition:
+C     loop through all filter value lines (for given) condition:
             DO ifvl = 1,nvl_fid(i1)
                ifiv = ifiv + 1
@@ Line 753: / Line 775: @@
                         IF (values(iv).EQ.rvind)
      +                       GOTO 400 ! next condition, no point in checking
-                                      ! other filter values lines
+                                      ! other filter value lines
                         IF (int_fid) THEN
                            IF (nint(values(iv)).EQ.fivI(ifiv,i2)) THEN
@@ Line 767: / 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
-            END DO              ! all filter values lines (for given) condition
         CONTINUE
          END DO                 ! all different conditions
@@ Line 788: / 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 803: / 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 809: / 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),
      +     T(maxlevel),TD(maxlevel),h(maxlevel),
      +     la(maxlevel),lo(maxlevel),tp(maxlevel), ! displacements (lat, lon, time)
      +     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 826: / Line 858: @@
       COMMON /COM_RECTANGLE/  x1,y1,x2,y2
-      CHARACTER*9 ctrim ! Function
+C     Functions
+      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 865: / Line 915: @@
          value = values(idx + (ksub-1)*kxelem)
-C     The missing value returned from libbufr might not be exactly equal to rvind
+C     Nothing to do if value is missing (note that the missing value
-         IF (ABS(value - rvind)/rvind.LE.0.001) value = rvind
+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
@@ Line 872: / 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
-            IF (value.NE.rvind) THEN
+     +           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)
-            END IF
+         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 890: / Line 960: @@
          ELSE IF (desc.EQ.4005) THEN ! Minute
             minute = value
-         ELSE IF (desc.EQ.5001) THEN ! Latitude (high accuracy)
+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) THEN ! Latitude (coarse accuracy)
+         ELSE IF (desc.EQ.5002 .AND. latitude.EQ.rvind) THEN ! Latitude (coarse accuracy)
             latitude = value
-         ELSE IF (desc.EQ.6001) THEN ! Longitude (high accuracy)
+         ELSE IF (desc.EQ.6001 .AND. longitude.EQ.rvind) THEN ! Longitude (high accuracy)
             longitude = value
-         ELSE IF (desc.EQ.6002) THEN ! Longitude (coarse accuracy)
+         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
+     +           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.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
             END IF
             P(n) = value
-         ELSE IF (desc.EQ.11001) THEN ! Wind direction
+C     All following descriptors come after pressure in 309052 (or Hirlam template)
-            D(n) = value
+         ELSE IF (n.GT.0 .AND. n.LE.maxlevel) THEN
-         ELSE IF (desc.EQ.11002) THEN ! Wind speed
+            IF (desc.EQ.11001) THEN ! Wind direction
-            F(n) = value
+               D(n) = value
-         ELSE IF (desc.EQ.12101) THEN ! Temperature/dry bulb temperature (16 bits)
+            ELSE IF (desc.EQ.11002) THEN ! Wind speed
-            T(n) = value
+               F(n) = value
-         ELSE IF (desc.EQ.12001) THEN ! Temperature/dry bulb temperature (12 bits)
+            ELSE IF (desc.EQ.12101) THEN ! Temperature/dry bulb temperature (16 bits)
-            T(n) = value
+               T(n) = value
-         ELSE IF (desc.EQ.12103) THEN ! Dew-point temperature (16 bits)
+            ELSE IF (desc.EQ.12001) THEN ! Temperature/dry bulb temperature (12 bits)
-            TD(n) = value
+               T(n) = value
-         ELSE IF (desc.EQ.12003) THEN ! Dew-point temperature (12 bits)
+            ELSE IF (desc.EQ.12103) THEN ! Dew-point temperature (16 bits)
-            TD(n) = value
+               TD(n) = value
-         ELSE IF (desc.EQ.10003.AND.value.NE.rvind) THEN ! Geopotential
+            ELSE IF (desc.EQ.12003) THEN ! Dew-point temperature (12 bits)
-            h(n) = value/9.81
+               TD(n) = value
-         ELSE IF (desc.EQ.10009) THEN ! Geopotential height
+            ELSE IF (desc.EQ.10003.AND.value.NE.rvind) THEN ! Geopotential
-            h(n) = value
+               h(n) = value/9.81
-         ELSE IF (desc.EQ.5015) THEN ! Latitude displacement
+            ELSE IF (desc.EQ.10009) THEN ! Geopotential height
-            la(n) = value
+               h(n) = value
-         ELSE IF (desc.EQ.6015) THEN ! Longitude displacement
+            ELSE IF (desc.EQ.5015) THEN ! Latitude displacement
-            lo(n) = value
+               la(n) = value
-         ELSE IF (desc.EQ.4086) THEN ! Long time period or displacement [second]
+            ELSE IF (desc.EQ.6015) THEN ! Longitude displacement
-C     In WMO templates 004086 comes before 7004 pressure
+               lo(n) = value
-            tp(n+1) = value
+            ELSE IF (desc.EQ.8001 .AND. value.NE.rvind) THEN ! Vertical sounding significance
-         ELSE IF (desc.EQ.8001 .AND. value.NE.rvind) THEN ! Vertical sounding significance
+C     In HIRLAM template (309007) 008001 comes AFTER 7004 pressure
-C     In HIRLAM templates 008001 comes AFTER 7004 pressure
+               CALL vss_8001(NINT(value),vss(n))
-            CALL vss_8001(NINT(value),vss(n))
+            ELSE IF (desc.EQ.11061) THEN ! Absolute wind shear in 1 km layer below
-         ELSE IF (desc.EQ.8042 .AND. value.NE.rvind) THEN ! Extended vertical sounding significance
+               wsb(n) = value
-C     In WMO templates 008042 comes BEFORE 7004 pressure
+            ELSE IF (desc.EQ.11062) THEN ! Absolute wind shear in 1 km layer above
-            CALL vss_8042(NINT(value),vss(n+1))
+               wsa(n) = value
-         ELSE IF (desc.EQ.11061) THEN ! Absolute wind shear in 1 km layer below
+            END IF
-            wsb(n) = value
-         ELSE IF (desc.EQ.11062) THEN ! Absolute wind shear in 1 km layer above
-            wsa(n) = value
          END IF
       END DO
@@ Line 951: / 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
-         ind = index(call_sign,' ') - 1
+         WRITE(*,*)
-         IF (ind.EQ.-1) ind = 9
+         WRITE(*,'(A,A)') 'call_sign=',call_sign(1:lenstr(call_sign,1))
-C     Remove trailing NULL characters, which some centres erronously
+      ELSE IF (flight_number.NE.missing8) THEN
-C     insert
+         WRITE(*,*)
-         DO WHILE (IACHAR(call_sign(ind:ind)).EQ.0)
+         WRITE(*,'(A,A)') 'aircraft=',
-            ind = ind - 1
+     +        flight_number(1:lenstr(flight_number,1))
-         END DO
-         WRITE(*,'(A,A)') 'call_sign=',
-     +        call_sign(1:ind)
       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 982: / 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 1000: / 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 1040: / 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 1055: / 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 1060: / 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*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,
-     +     RR_1,RR_3,RR_6,RR_12,RR_24,RS,SA,SG,SI,
+     +     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,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
+     +     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 1089: / 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 1099: / Line 1419: @@
       COMMON /COM_RECTANGLE/  x1,y1,x2,y2
-      CHARACTER*20 ctrim ! Function
+C     Functions
+      CHARACTER*32 ctrim ! length must be >= longest variable ctrim is used for
+      INTEGER lenstr
       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
-      spc8 = '        '
+      missing32 = missing16 // missing16
-      spc9 = '         '
-      spc20 = '                    '
       cloud_type_count = 0
@@ Line 1114: / Line 1437: @@
       num_cloud_layers = -1
       surface_data = .TRUE.
+      time_of_last_position = .FALSE.
 C     Initialize all parameters to missing values
@@ Line 1120: / Line 1444: @@
       icao_id = missing8
       call_sign = missing9
+      wigos_issuer = missing5
+      wigos_issueno = missing5
+      wigos_localid = missing16
       name = missing20
+      long_name = missing32
       AA = rvind
       BI = rvind
@@ Line 1135: / Line 1463: @@
       DW1 = rvind
       DW2 = rvind
+      EE = rvind
       ES = rvind
-      E = rvind
-      EM = rvind
       EV_1 = rvind
       EV_24 = rvind
@@ Line 1152: / Line 1479: @@
       HW1 = rvind
       HW2 = rvind
+      HWA = rvind
       NH = rvind
       NN = rvind
@@ Line 1164: / Line 1492: @@
       PW1 = rvind
       PW2 = rvind
+      PWA = rvind
+      QD = rvind
+      QE = rvind
+      QK = rvind
+      QL = rvind
+      QO = rvind
+      QS = rvind
+      QD_24 = rvind
+      QE_24 = rvind
+      QK_24 = rvind
+      QL_24 = rvind
+      QO_24 = rvind
+      QS_24 = rvind
       RR_1 = rvind
       RR_3 = rvind
@@ Line 1173: / Line 1514: @@
       SG = rvind
       SI = rvind
+      SS_24 = rvind
       TA = rvind
       TAN_12 = rvind
@@ Line 1196: / Line 1538: @@
       longitude = rvind
       height = rvind
+      hp = rvind
       vert_sign_first = rvind
       II = rvind
@@ Line 1205: / 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 1221: / Line 1568: @@
          value = values(idx + (ksub-1)*kxelem)
-C     The missing value returned from libbufr might not be exactly equal to rvind
+         IF (ABS(value - rvind)/rvind.LE.0.001) THEN
-         IF (ABS(value - rvind)/rvind.LE.0.001) value = rvind
+C     Missing value. Only a few descriptors require a reset of
+C     corresponding parameter in this case. Note that the missing value
+C     returned from libbufr might not be exactly equal to rvind
+            IF (desc.EQ.4024) THEN ! Time period or displacement [hour]
+               hour_p = rvind
+            ELSE IF (desc.EQ.4025) THEN ! Time period or displacement [minute]
+               minute_p = rvind
+C     Delayed descriptor replication factor should never be missing
+            ELSE IF (desc.EQ.31001) THEN
+               IF (idx.GT.1 .AND. ktdexp(idx - 1).EQ.20012) THEN
+                  WRITE(*,*) 'WARNING: delayed descriptor replication'
+     +                 // ' factor after 020012 undefined!!!'
+                  bad_cloud_data = .TRUE.
+               ELSE
+                  WRITE(*,*) 'WARNING: delayed descriptor replication'
+     +                 // ' factor 31001 undefined!!!'
+               END IF
+C     Some counting needed for clouds even for missing values
+            ELSE IF (desc.EQ.20012 .AND. .NOT.bad_cloud_data) THEN ! Cloud type
+               cloud_type_count = cloud_type_count + 1
+               IF (cloud_type_count.GT.3) THEN
+                  cloud_layer = cloud_type_count - 3
+               END IF
+            END IF
+            CYCLE
+         END IF
+C     Continue the loop for non-missing value. For most variables we
+C     choose to not set the parameter if set before, because if a
+C     descriptor unexpectedly occurs more than once it is likely that
+C     the first occurrence is the 'standard' use of the descriptor,
+C     while the later occurrence(s) might for example be due to data
+C     required by regional or national reporting practices, added after
+C     a standard WMO template
          IF (desc.EQ.4024) THEN ! Time period or displacement [hour]
             hour_p = value
          ELSE IF (desc.EQ.4025) THEN ! Time period or displacement [minute]
             minute_p = value
-         ELSE IF (desc.EQ.1001) THEN ! WMO block number
+         ELSE IF (desc.EQ.8021) THEN ! Time significance
-            IF (value.NE.rvind .AND. II.EQ.rvind) THEN
+            IF (NINT(value).EQ.26) THEN
-               II = value
+               time_of_last_position = .TRUE.
-            END IF
+            ELSE
-         ELSE IF (desc.EQ.1002) THEN ! WMO station number
+               time_of_last_position = .FALSE.
-            IF (value.NE.rvind .AND. iii.EQ.rvind) THEN
-               iii = value
             END IF
+         ELSE IF (desc.EQ.1001 .AND. II.EQ.rvind) THEN ! WMO block number
+            IF (value.GE.0 .AND. value.LT.100) II = value
+         ELSE IF (desc.EQ.1002 .AND. iii.EQ.rvind) THEN ! WMO station number
+            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
+            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
+            IF (value.GE.0) national_number = value
          ELSE IF (desc.EQ.1015) THEN  ! Station or site name
-            IF (value.NE.rvind) THEN
+            cidx = int(value/1000)
-               cidx = int(value/1000)
+            name = cvals(cidx)  ! CCITTIA5 data
-               IF (cvals(cidx).NE.spc20) THEN
+            name = ctrim(name,20,missing20)
-                  name = cvals(cidx) ! CCITTIA5 data
+         ELSE IF (desc.EQ.1019) THEN  ! Long station or site name
-                  name = ctrim(name,20,missing20)
+            cidx = int(value/1000)
-               END IF
+            long_name = cvals(cidx) ! CCITTIA5 data
-            END IF
+            long_name = ctrim(long_name,32,missing32)
+         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
-            IF (value.NE.rvind .AND. ix.EQ.rvind) THEN
+            IF (ix.EQ.rvind) THEN
                ix = value
             END IF
          ELSE IF (desc.EQ.4001) THEN ! Year
-            IF (value.NE.rvind .AND. 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 (value.NE.rvind .AND. 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 (value.NE.rvind .AND. 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 (value.NE.rvind .AND. 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 (value.NE.rvind .AND. 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)
-            IF (value.NE.rvind .AND. latitude.EQ.rvind) THEN
+     +           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 (value.NE.rvind .AND. latitude.EQ.rvind) THEN
+     +           desc.EQ.6002) THEN ! Longitude (coarse accuracy)
-               latitude = value
+            IF (longitude.EQ.rvind) THEN
-            END IF
-         ELSE IF (desc.EQ.6001) THEN ! Longitude (high accuracy)
-            IF (value.NE.rvind .AND. longitude.EQ.rvind) THEN
-               longitude = value
-            END IF
-         ELSE IF (desc.EQ.6002) THEN ! Longitude (coarse accuracy)
-            IF (value.NE.rvind .AND. longitude.EQ.rvind) THEN
                longitude = value
             END IF
-         ELSE IF (desc.EQ.7001) THEN ! Height of station
+         ELSE IF (desc.EQ.7001.OR.  ! Height of station
-            IF (value.NE.rvind .AND. height.EQ.rvind) THEN
+     +           desc.EQ.7030) THEN ! Height of station ground above mean sea level
+            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 (value.NE.rvind .AND. height.EQ.rvind) THEN
+            IF (hp.EQ.rvind) THEN
-               height = value
+               hp = value
             END IF
          ELSE IF (desc.EQ.10004) THEN ! Pressure
-            IF (value.NE.rvind .AND. PO.EQ.rvind) THEN
+            IF (PO.EQ.rvind) THEN
                PO = value
             END IF
          ELSE IF (desc.EQ.10051) THEN ! Pressure reduced to mean sea level
-            IF (value.NE.rvind .AND. PR.EQ.rvind) THEN
+            IF (PR.EQ.rvind) THEN
                PR = value
             END IF
          ELSE IF (desc.EQ.10061) THEN ! 3-hour pressure change
-            IF (value.NE.rvind .AND. PP.EQ.rvind) THEN
+            IF (PP.EQ.rvind) THEN
                PP = value
             END IF
          ELSE IF (desc.EQ.10063) THEN ! Characteristic of pressure tendency
-            IF (value.NE.rvind .AND. AA.EQ.rvind) THEN
+            IF (AA.EQ.rvind) THEN
                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
-            IF (value.NE.rvind .AND. DD.EQ.rvind) THEN
+     +           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 (value.NE.rvind .AND. DD.EQ.rvind) THEN
+     +           desc.EQ.11002) THEN ! Wind speed
-               DD = value
+            IF (FF.EQ.rvind) THEN
-            END IF
-         ELSE IF (desc.EQ.11012) THEN ! Wind speed at 10 m
-            IF (value.NE.rvind .AND. FF.EQ.rvind) THEN
-               FF = value
-            END IF
-         ELSE IF (desc.EQ.11002) THEN ! Wind speed
-            IF (value.NE.rvind .AND. FF.EQ.rvind) THEN
                FF = value
             END IF
          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 1345: / 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 1369: / 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
 C     termins 0,3,6,9... only
                mm = NINT(minute_p)
                hh = NINT(hour)
@@ Line 1387: / 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)
-            IF (value.NE.rvind .AND. TA.EQ.rvind) THEN
+     +          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 (value.NE.rvind .AND. 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 (TD.EQ.rvind) THEN
-            IF (value.NE.rvind .AND. TA.EQ.rvind) THEN
-               TA = value
-            END IF
-         ELSE IF (desc.EQ.12001) THEN ! Temperature/dry bulb temperature (12 bits)
-            IF (value.NE.rvind .AND. TA.EQ.rvind) THEN
-               TA = value
-            END IF
-         ELSE IF (desc.EQ.12106) THEN ! Dew-point temperature at 2m (16 bits)
-            IF (value.NE.rvind .AND. TD.EQ.rvind) THEN
-               TD = value
-            END IF
-         ELSE IF (desc.EQ.12006) THEN ! Dew-point temperature at 2m (12 bits)
-            IF (value.NE.rvind .AND. TD.EQ.rvind) THEN
-               TD = value
-            END IF
-         ELSE IF (desc.EQ.12103) THEN ! Dew-point temperature (16 bits)
-            IF (value.NE.rvind .AND. TD.EQ.rvind) THEN
-               TD = value
-            END IF
-         ELSE IF (desc.EQ.12003) THEN ! Dew-point temperature (12 bits)
-            IF (value.NE.rvind .AND. TD.EQ.rvind) THEN
                TD = value
             END IF
-         ELSE IF (desc.EQ.12113) THEN ! Ground minimum temperature at 2m (data width 16 bits)
+         ELSE IF (desc.EQ.12113.OR.  ! Ground minimum temperature at 2m (data width 16 bits)
-            IF (value.NE.rvind .AND. TGN_12.EQ.rvind) THEN
+     +           desc.EQ.12013) THEN ! Ground minimum temperature at 2m (12 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 (value.NE.rvind .AND. TGN_12.EQ.rvind) THEN
+     +           desc.EQ.12014) THEN ! Maximum temperature at 2m, past 12 hours (12 bits)
-               TGN_12 = value
+            IF (TAX_12.EQ.rvind) THEN
-            END IF
-         ELSE IF (desc.EQ.12114) THEN ! Maximum temperature at 2m, past 12 hours (16 bits)
-            IF (value.NE.rvind .AND. 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 (value.NE.rvind .AND. TAX_12.EQ.rvind) THEN
                TAX_12 = value
             END IF
          ELSE IF (desc.EQ.12111) THEN ! Maximum temperature at height and over period specified
-            IF (value.NE.rvind .AND. TAX_12.EQ.rvind .AND. idx.GT.2
+            IF (TAX_12.EQ.rvind .AND. idx.GT.2
      +           .AND. ktdexp(idx-1).EQ.4024
      +           .AND. NINT(values(idx-1 + (ksub-1)*kxelem)).EQ.0
@@ Line 1447: / 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 (value.NE.rvind .AND. TAN_12.EQ.rvind) THEN
+     +           desc.EQ.12015) THEN ! Minimum temperature at 2m, past 12 hours (12 bits)
-               TAN_12 = value
+            IF (TAN_12.EQ.rvind) THEN
-            END IF
-         ELSE IF (desc.EQ.12015) THEN ! Minimum temperature at 2m, past 12 hours (12 bits)
-            IF (value.NE.rvind .AND. TAN_12.EQ.rvind) THEN
                TAN_12 = value
             END IF
          ELSE IF (desc.EQ.12112) THEN ! Minimum temperature at height and over period specified
-            IF (value.NE.rvind .AND. TAN_12.EQ.rvind .AND. idx.GT.2
+            IF (TAN_12.EQ.rvind .AND. idx.GT.2
      +           .AND. ktdexp(idx-1).EQ.4024
      +           .AND. values(idx-1 + (ksub-1)*kxelem).EQ.0
@@ Line 1468: / Line 1826: @@
 C     Do we also need to consider 12022 'Minimum temperature at 2m'?
          ELSE IF (desc.EQ.13003) THEN ! Relative humidity
-            IF (value.NE.rvind .AND. UU.EQ.rvind) THEN
+            IF (UU.EQ.rvind) THEN
                UU = value
             END IF
          ELSE IF (desc.EQ.20001) THEN ! Horizontal visibility
-            IF (value.NE.rvind .AND. VV.EQ.rvind) THEN
+            IF (VV.EQ.rvind) THEN
                VV = value
             END IF
          ELSE IF (desc.EQ.20003) THEN ! Present weather
-            IF (value.NE.rvind .AND. WW.EQ.rvind) THEN
+            IF (WW.EQ.rvind) THEN
                WW = value
             END IF
          ELSE IF (desc.EQ.20004) THEN ! Past weather (1)
-            IF (value.NE.rvind .AND. W1.EQ.rvind) THEN
+            IF (W1.EQ.rvind) THEN
                W1 = value
             END IF
          ELSE IF (desc.EQ.20005) THEN ! Past weather (2)
-            IF (value.NE.rvind .AND. W2.EQ.rvind) THEN
+            IF (W2.EQ.rvind) THEN
                W2 = value
             END IF
          ELSE IF (desc.EQ.20010) THEN ! Cloud cover (total)
-            IF (value.NE.rvind .AND. NN.EQ.rvind) THEN
+            IF (NN.EQ.rvind) THEN
                NN = value
             END IF
@@ Line 1495: / 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 (value.EQ.rvind) THEN
+               IF (NINT(value).LE.4) THEN
-                  WRITE(*,*) 'WARNING: delayed descriptor replication'
+                  num_cloud_layers = NINT(value)
-     +                 // ' factor after 020012 undefined!!!'
+               ELSE
                   bad_cloud_data = .TRUE.
-               ELSE
-                  IF (NINT(value).LE.4) THEN
-                     num_cloud_layers = NINT(value)
-                  ELSE
-                     bad_cloud_data = .TRUE.
-                  END IF
                END IF
             END IF
          ELSE IF (desc.EQ.8002 .AND. .NOT.bad_cloud_data) THEN ! Vertical significance (surface observations)
             IF (cloud_type_count.EQ.0) THEN ! First occurrence
-               IF (value.NE.rvind .AND. vert_sign_first.EQ.rvind) THEN
+               IF (vert_sign_first.EQ.rvind) THEN
                   vert_sign_first = value
                END IF
@@ Line 1519: / 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 (value.NE.rvind .AND. NH.EQ.rvind) 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
                END IF
@@ Line 1542: / 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 (value.NE.rvind .AND. 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 (value.NE.rvind .AND. 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 (value.NE.rvind .AND. 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
-               IF (value.NE.rvind .AND. HL.EQ.rvind) THEN
+               HS(num_cloud_layers) = value
+            ELSE IF (cloud_type_count.EQ.0) THEN ! First occurrence
+               IF (HL.EQ.rvind) THEN
                   HL = value
                END IF
@@ Line 1598: / 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
          ELSE IF (desc.EQ.13023) THEN ! Total precipitation past 24 hours
-            IF (value.NE.rvind .AND. RR_24.EQ.rvind) THEN
+            IF (RR_24.EQ.rvind) THEN
                RR_24 = value
             END IF
          ELSE IF (desc.EQ.13022) THEN ! Total precipitation past 12 hours
-            IF (value.NE.rvind .AND. RR_12.EQ.rvind) THEN
+            IF (RR_12.EQ.rvind) THEN
                RR_12 = value
             END IF
          ELSE IF (desc.EQ.13021) THEN ! Total precipitation past 6 hours
-            IF (value.NE.rvind .AND. RR_6.EQ.rvind) THEN
+            IF (RR_6.EQ.rvind) THEN
                RR_6 = value
             END IF
          ELSE IF (desc.EQ.13020) THEN ! Total precipitation past 3 hours
-            IF (value.NE.rvind .AND. RR_3.EQ.rvind) THEN
+            IF (RR_3.EQ.rvind) THEN
                RR_3 = value
             END IF
          ELSE IF (desc.EQ.13019) THEN ! Total precipitation past 1 hour
-            IF (value.NE.rvind .AND. RR_1.EQ.rvind) THEN
+            IF (RR_1.EQ.rvind) THEN
                RR_1 = value
             END IF
          ELSE IF (desc.EQ.13011) THEN ! Total precipitation/total water equivalent
-            IF (value.NE.rvind .AND. hour_p.NE.rvind) THEN
+            IF (hour_p.NE.rvind) THEN
                hh = NINT(hour_p)
                IF (hh.EQ.-24) THEN
@@ Line 1648: / 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
-            END IF
          ELSE IF (desc.EQ.20062) THEN ! State of the ground (with or without snow)
-            IF (value.NE.rvind) THEN
+            IF (EE.EQ.rvind) THEN
-               IF (NINT(value).LE.10) THEN ! E in TAC (3Ejjj)
+               EE = value
-                  E = value
+            END IF
-                  IF (SA.EQ.rvind) THEN
+            IF (NINT(value).LE.10 .AND. SA.EQ.rvind) THEN
-                     SA = 0
+               SA = 0
-                  END IF
+            END IF
-               ELSE IF (NINT(value).LE.20) THEN ! E' in TAC (4Esss)
+C     Equating 013012 with SS_24 is dubious generally, but this is how
-                  EM = value
+C     SS_24 is encoded in kvalobs for Norwegian avalanche stations (and
-               END IF
+C     we have no better kvalobs parameter for depth of fresh snow anyway)
+         ELSE IF (desc.EQ.13012) THEN ! Depth of fresh snow
+            IF (SS_24.EQ.rvind) THEN
+               SS_24 = value
             END IF
          ELSE IF (desc.EQ.14031) THEN ! Total sunshine
-            IF (value.NE.rvind .AND. hour_p.NE.rvind) THEN
+            IF (hour_p.NE.rvind) THEN
                hh = NINT(hour_p)
                IF (hh.EQ.-1) THEN
@@ Line 1674: / Line 2024: @@
             END IF
          ELSE IF (desc.EQ.13033) THEN ! Evaporation/evapotranspiration
-            IF (value.NE.rvind .AND. hour_p.NE.rvind) THEN
+            IF (hour_p.NE.rvind) THEN
                hh = NINT(hour_p)
                IF (hh.EQ.-1) THEN
@@ Line 1680: / Line 2030: @@
                ELSE IF (hh.EQ.-24) THEN
                   EV_24 = value
+               END IF
+            END IF
+         ELSE IF (desc.EQ.14002) THEN ! Long-wave radiation
+            IF (hour_p.NE.rvind) THEN
+               hh = NINT(hour_p)
+               IF (hh.EQ.-1) THEN
+                  QL = value
+               ELSE IF (hh.EQ.-24) THEN
+                  QL_24 = value
+               END IF
+            END IF
+         ELSE IF (desc.EQ.14004) THEN ! Short-wave radiation
+            IF (hour_p.NE.rvind) THEN
+               hh = NINT(hour_p)
+               IF (hh.EQ.-1) THEN
+                  QK = value
+               ELSE IF (hh.EQ.-24) THEN
+                  QK_24 = value
+               END IF
+            END IF
+         ELSE IF (desc.EQ.14016) THEN ! Net radiation
+            IF (hour_p.NE.rvind) THEN
+               hh = NINT(hour_p)
+               IF (hh.EQ.-1) THEN
+                  QE = value
+               ELSE IF (hh.EQ.-24) THEN
+                  QE_24 = value
+               END IF
+            END IF
+         ELSE IF (desc.EQ.14028) THEN ! Global solar radiation
+            IF (hour_p.NE.rvind) THEN
+               hh = NINT(hour_p)
+               IF (hh.EQ.-1) THEN
+                  QO = value
+               ELSE IF (hh.EQ.-24) THEN
+                  QO_24 = value
+               END IF
+            END IF
+         ELSE IF (desc.EQ.14029) THEN ! Diffuse solar radiation
+            IF (hour_p.NE.rvind) THEN
+               hh = NINT(hour_p)
+               IF (hh.EQ.-1) THEN
+                  QD = value
+               ELSE IF (hh.EQ.-24) THEN
+                  QD_24 = value
+               END IF
+            END IF
+         ELSE IF (desc.EQ.14030) THEN ! Direct solar radiation
+            IF (hour_p.NE.rvind) THEN
+               hh = NINT(hour_p)
+               IF (hh.EQ.-1) THEN
+                  QS = value
+               ELSE IF (hh.EQ.-24) THEN
+                  QS_24 = value
                END IF
             END IF
 C     Special for high altitude stations
          ELSE IF (desc.EQ.7004) THEN ! Pressure (location class)
-            IF (value.NE.rvind .AND. a3.EQ.rvind) THEN
+            IF (a3.EQ.rvind) THEN
                a3 = value
             END IF
          ELSE IF (desc.EQ.10009) THEN ! Geopotential height
-            IF (value.NE.rvind .AND. hhh.EQ.rvind) THEN
+            IF (hhh.EQ.rvind) THEN
                hhh = value
-            END IF
-         ELSE IF (desc.EQ.10008) THEN ! Geopotential (20 bits)
-            IF (value.NE.rvind .AND. hhh.EQ.rvind) THEN
-               hhh = value * 9.8
-            END IF
-         ELSE IF (desc.EQ.10003) THEN ! Geopotential (17 bits)
-            IF (value.NE.rvind .AND. 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
-            IF (value.NE.rvind .AND. ds.EQ.rvind) THEN
+            IF (ds.EQ.rvind) THEN
                ds = value
             END IF
          ELSE IF (desc.EQ.1013) THEN ! Speed of motion of moving observing platform
-            IF (value.NE.rvind .AND. vs.EQ.rvind) THEN
+            IF (vs.EQ.rvind) THEN
                vs = value
             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)
-            IF (value.NE.rvind .AND. TW.EQ.rvind .AND. surface_data)THEN
+     +           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 (value.NE.rvind .AND. TW.EQ.rvind
+     +           desc.EQ.12005) THEN ! Wet-bulb temperature (12 bits)
-     +           .AND. surface_data) THEN
+            IF (TbTbTb.EQ.rvind) THEN
-               TW = value
-            END IF
-         ELSE IF (desc.EQ.12102) THEN ! Wet-bulb temperature (16 bits)
-            IF (value.NE.rvind .AND. 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 (value.NE.rvind .AND. 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 (value.NE.rvind .AND. PW.EQ.rvind) THEN
+            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 (value.NE.rvind .AND. HW.EQ.rvind) THEN
+            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
-            IF (value.NE.rvind .AND. SG.EQ.rvind) THEN
+            IF (SG.EQ.rvind) THEN
                SG = value
             END IF
          ELSE IF (desc.EQ.20033) THEN ! Cause of ice accretion
-            IF (value.NE.rvind .AND. XIS.EQ.rvind) THEN
+            IF (XIS.EQ.rvind) THEN
                XIS = value
             END IF
          ELSE IF (desc.EQ.20031) THEN ! Ice deposit (thickness)
-            IF (value.NE.rvind .AND. ES.EQ.rvind) THEN
+            IF (ES.EQ.rvind) THEN
                ES = value
             END IF
          ELSE IF (desc.EQ.20032) THEN ! Rate of ice accretion
-            IF (value.NE.rvind .AND. RS.EQ.rvind) THEN
+            IF (RS.EQ.rvind) THEN
                RS = value
             END IF
          ELSE IF (desc.EQ.20034) THEN ! Sea ice concentration
-            IF (value.NE.rvind .AND. CI.EQ.rvind) THEN
+            IF (CI.EQ.rvind) THEN
                CI = value
             END IF
          ELSE IF (desc.EQ.20037) THEN ! Ice development
-            IF (value.NE.rvind .AND. SI.EQ.rvind) THEN
+            IF (SI.EQ.rvind) THEN
                SI = value
             END IF
          ELSE IF (desc.EQ.20035) THEN ! Amount and type of ice
-            IF (value.NE.rvind .AND. BI.EQ.rvind) THEN
+            IF (BI.EQ.rvind) THEN
                BI = value
             END IF
          ELSE IF (desc.EQ.20038) THEN ! Bearing of ice edge
-            IF (value.NE.rvind .AND. DI.EQ.rvind) THEN
+            IF (DI.EQ.rvind) THEN
                DI = value
             END IF
          ELSE IF (desc.EQ.20036) THEN ! Ice situation
-            IF (value.NE.rvind .AND. ZI.EQ.rvind) THEN
+            IF (ZI.EQ.rvind) THEN
                ZI = value
             END IF
          ELSE IF (desc.EQ.1005) THEN ! Buoy/platform identifier
-            IF (value.NE.rvind .AND. 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
-            IF (value.NE.rvind .AND. wmo_region_number.EQ.rvind) THEN
+            IF (wmo_region_number.EQ.rvind) THEN
                wmo_region_number = value
             END IF
          ELSE IF (desc.EQ.1020) THEN ! WMO region sub-area
-            IF (value.NE.rvind .AND. wmo_region_subarea.EQ.rvind) THEN
+            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 (value.NE.rvind .AND. PH.EQ.rvind) THEN
+            IF (PH.EQ.rvind) THEN
                PH = value
             END IF
@@ Line 1846: / 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)
       ELSE IF (call_sign.NE.missing9) THEN
          ind = index(call_sign,' ') - 1
@@ Line 1857: / 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, call_sign and buoy_id',
+            WRITE(*,*)
-     +           ' 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)
          END IF
          RETURN
@@ Line 1887: / Line 2328: @@
       ENDIF
       IF (icao_id.NE.missing8) THEN
-         WRITE(*,'(A,A)') 'icao_id=',icao_id
+         WRITE(*,'(A,A)') 'icao_id=',icao_id(1:lenstr(icao_id,1))
       END IF
       IF (name.NE.missing20) THEN
-         WRITE(*,'(A,A)') 'name=',name
+         WRITE(*,'(A,A)') 'name=',name(1:lenstr(name,1))
       END IF
-      IF (NINT(ix).EQ.0) THEN
+      IF (long_name.NE.missing32) THEN
-         WRITE(*,'(A,A)') 'type=Automatic'
+         WRITE(*,'(A,A)') 'name=',long_name(1:lenstr(long_name,1))
-      ELSE IF (NINT(ix).EQ.1) THEN
+      END IF
-         WRITE(*,'(A,A)') 'type=Manned'
+      IF (ix.NE.rvind) THEN
-      ELSE IF (NINT(ix).EQ.2) THEN
+         IF (NINT(ix).EQ.0) THEN
-         WRITE(*,'(A,A)') 'type=Hybrid'
+            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 (longitude.NE.rvind) THEN
@@ Line 1908: / 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 1914: / 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 2002: / 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 2060: / 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 2114: / 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 2120: / 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 2179: / Line 2609: @@
       IF (OT_24.NE.rvind) THEN
          WRITE(*,'(A,I8)') 'OT_24=',NINT(OT_24)
+      END IF
+      IF (QE.NE.rvind) THEN
+         WRITE(*,'(A,F11.2)') 'QE=',QE/3600 ! Wh/m2
+      END IF
+      IF (QO.NE.rvind) THEN
+         WRITE(*,'(A,F11.2)') 'QO=',QO/3600
+      END IF
+      IF (QL.NE.rvind) THEN
+         WRITE(*,'(A,F11.2)') 'QL=',QL/3600
+      END IF
+      IF (QK.NE.rvind) THEN
+         WRITE(*,'(A,F11.2)') 'QK=',QK/3600
+      END IF
+      IF (QD.NE.rvind) THEN
+         WRITE(*,'(A,F11.2)') 'QD=',QD/3600
+      END IF
+      IF (QS.NE.rvind) THEN
+         WRITE(*,'(A,F11.2)') 'QS=',QS/3600
+      END IF
+      IF (QE_24.NE.rvind) THEN
+         WRITE(*,'(A,F8.2)') 'QE_24=',QE_24/3600 ! Wh/m2
+      END IF
+      IF (QO_24.NE.rvind) THEN
+         WRITE(*,'(A,F8.2)') 'QO_24=',QO_24/3600
+      END IF
+      IF (QL_24.NE.rvind) THEN
+         WRITE(*,'(A,F8.2)') 'QL_24=',QL_24/3600
+      END IF
+      IF (QK_24.NE.rvind) THEN
+         WRITE(*,'(A,F8.2)') 'QK_24=',QK_24/3600
+      END IF
+      IF (QD_24.NE.rvind) THEN
+         WRITE(*,'(A,F8.2)') 'QD_24=',QD_24/3600
+      END IF
+      IF (QS_24.NE.rvind) THEN
+         WRITE(*,'(A,F8.2)') 'QS_24=',QS_24/3600
       END IF
       IF (EV_1.NE.rvind) THEN
@@ Line 2220: / 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 2226: / Line 2693: @@
       COMMON /COM_RECTANGLE/  x1,y1,x2,y2
-      CHARACTER*9 ctrim ! Function
+C     Functions
+      CHARACTER*9 ctrim ! length must be >= longest variable ctrim is used for
       one_bits = CHAR(255)
@@ Line 2250: / Line 2718: @@
       DO n=1,maxlevel
          z(n) = rvind
+         w(n) = rvind
          T(n) = rvind
          s(n) = rvind
@@ Line 2262: / Line 2731: @@
          value = values(idx + (ksub-1)*kxelem)
-C     The missing value returned from libbufr might not be exactly equal to rvind
+         IF (ABS(value - rvind)/rvind.LE.0.001) THEN
-         IF (ABS(value - rvind)/rvind.LE.0.001) value = rvind
+C     Missing value, nothing to do
+            CYCLE
+         END IF
+C     Continue the loop for non missing values
          IF (desc.EQ.4001) THEN ! Year
-            IF (value.NE.rvind .AND. year.EQ.rvind) THEN
+            IF (year.EQ.rvind) THEN
                year = value
-            END IF
+           END IF
          ELSE IF (desc.EQ.4002) THEN ! Month
-            IF (value.NE.rvind .AND. month.EQ.rvind) THEN
+            IF (month.EQ.rvind) THEN
                month = value
             END IF
          ELSE IF (desc.EQ.4003) THEN ! Day
-            IF (value.NE.rvind .AND. day.EQ.rvind) THEN
+            IF (day.EQ.rvind) THEN
                day = value
             END IF
          ELSE IF (desc.EQ.4004) THEN ! Hour
-            IF (value.NE.rvind .AND. hour.EQ.rvind) THEN
+            IF (hour.EQ.rvind) THEN
                hour = value
             END IF
          ELSE IF (desc.EQ.4005) THEN ! Minute
-            IF (value.NE.rvind .AND. minute.EQ.rvind) THEN
+            IF (minute.EQ.rvind) THEN
                minute = value
             END IF
-         ELSE IF (desc.EQ.5001) THEN ! Latitude (high accuracy)
+         ELSE IF (desc.EQ.5001.OR.  ! Latitude (high accuracy)
-            IF (value.NE.rvind .AND. latitude.EQ.rvind) THEN
+     +           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 (value.NE.rvind .AND. latitude.EQ.rvind) THEN
+     +           desc.EQ.6002) THEN ! Longitude (coarse accuracy)
-               latitude = value
+            IF (longitude.EQ.rvind) THEN
-            END IF
-         ELSE IF (desc.EQ.6001) THEN ! Longitude (high accuracy)
-            IF (value.NE.rvind .AND. longitude.EQ.rvind) THEN
-               longitude = value
-            END IF
-         ELSE IF (desc.EQ.6002) THEN ! Longitude (coarse accuracy)
-            IF (value.NE.rvind .AND. longitude.EQ.rvind) THEN
                longitude = value
             END IF
@@ Line 2312: / 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 (value.NE.rvind .AND. 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 (value.NE.rvind .AND. 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 (value.NE.rvind .AND. 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
-               d(n) = value
+               IF (n.GT.0 .AND. n.LE.maxlevel) THEN
+                  d(n) = value
+               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
-               v(n) = value
+               IF (n.GT.0 .AND. n.LE.maxlevel) THEN
+                  v(n) = value
+               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 (value.NE.rvind .AND. sal.EQ.rvind) THEN
+            ELSE IF (n.EQ.0 .AND. sal.EQ.rvind) THEN
                sal = value
             END IF
          ELSE IF (desc.EQ.22011) THEN !  Period of waves (BUFR doesn't distinguish between PW and PWA)
-            IF (value.NE.rvind .AND. PW.EQ.rvind) THEN
+            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)
-            IF (value.NE.rvind .AND. HW.EQ.rvind) THEN
+            IF (HW.EQ.rvind) THEN
                HW = value
             END IF
          ELSE IF (desc.EQ.1005) THEN ! Buoy/platform identifier
-            IF (value.NE.rvind .AND. buoy_id.EQ.rvind) THEN
+            IF (buoy_id.EQ.rvind) THEN
                buoy_id = value
             END IF
          ELSE IF (desc.EQ.1003) THEN ! WMO region number/geographical area
-            IF (value.NE.rvind .AND. wmo_region_number.EQ.rvind) THEN
+            IF (wmo_region_number.EQ.rvind) THEN
                wmo_region_number = value
             END IF
          ELSE IF (desc.EQ.1020) THEN ! WMO region sub-area
-            IF (value.NE.rvind .AND. wmo_region_subarea.EQ.rvind) THEN
+            IF (wmo_region_subarea.EQ.rvind) THEN
                wmo_region_subarea = value
             END IF
@@ Line 2369: / Line 2849: @@
                call_sign = cvals(cidx) ! CCITTIA5 data
                call_sign = ctrim(call_sign,9,missing9)
+            END IF
+         ELSE IF (desc.EQ.1087) THEN ! WMO marine observing platform extended identifier
+            IF (buoy_id.EQ.rvind) THEN
+               buoy_id = value
             END IF
          END IF
@@ Line 2381: / 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 2398: / 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 2443: / 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 2449: / 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 2485: / 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 2496: / Line 3006: @@
       REAL*8 value
       INTEGER desc
 C     Variables used for geographical filtering av observations
 C     (--rectangle option set)
@@ Line 2503: / Line 3013: @@
 C     Functions
-      CHARACTER*8 ctrim
+      INTEGER lenstr
+      CHARACTER*8 ctrim ! length must be >= longest variable ctrim is used for
       CHARACTER*3 phase_8004,phase_8009
@@ Line 2510: / 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 2531: / Line 3041: @@
       flight_level = rvind
       phase = rvind
+      mixing_ratio = rvind
 C     Loop through all expanded descriptors
       DO idx=1,ktdexl
@@ Line 2537: / Line 3048: @@
          value = values(idx + (ksub-1)*kxelem)
-C     The missing value returned from libbufr might not be exactly equal to rvind
+         IF (ABS(value - rvind)/rvind.LE.0.001) THEN
-         IF (ABS(value - rvind)/rvind.LE.0.001) value = rvind
+C     Missing value, nothing to do
+            CYCLE
+         END IF
+C     Continue the loop for non missing values
          IF (desc.EQ.1008) THEN ! Aircraft registration number or other identification
-            IF (value.NE.rvind) THEN
+            cidx = int(value/1000)
-               cidx = int(value/1000)
+            aircraft = cvals(cidx) ! CCITTIA5 data
-               IF (cvals(cidx).NE.spc8) THEN
+            aircraft = ctrim(aircraft,8,missing8)
-                  aircraft = cvals(cidx) ! CCITTIA5 data
-                  aircraft = ctrim(aircraft,8,missing8)
-               END IF
-            END IF
          ELSE IF (desc.EQ.1006) THEN ! Aircraft flight number
-            IF (value.NE.rvind) THEN
+            cidx = int(value/1000)
-               cidx = int(value/1000)
+            flight_number = cvals(cidx) ! CCITTIA5 data
-               IF (cvals(cidx).NE.spc8) THEN
+            flight_number = ctrim(flight_number,8,missing8)
-                  flight_number = cvals(cidx) ! CCITTIA5 data
-                  flight_number = ctrim(flight_number,8,missing8)
-               END IF
-            END IF
          ELSE IF (desc.EQ.1023) THEN ! Observation sequence number
-            IF (value.NE.rvind .AND. osn.EQ.rvind) THEN
+            IF (osn.EQ.rvind) osn = value
-               osn = value
-            END IF
          ELSE IF (desc.EQ.4001) THEN ! Year
-            IF (value.NE.rvind .AND. year.EQ.rvind) THEN
+            IF (year.EQ.rvind) THEN
                year = value
             END IF
          ELSE IF (desc.EQ.4002) THEN ! Month
-            IF (value.NE.rvind .AND. month.EQ.rvind) THEN
+            IF (month.EQ.rvind) THEN
                month = value
             END IF
          ELSE IF (desc.EQ.4003) THEN ! Day
-            IF (value.NE.rvind .AND. day.EQ.rvind) THEN
+            IF (day.EQ.rvind) THEN
                day = value
             END IF
          ELSE IF (desc.EQ.4004) THEN ! Hour
-            IF (value.NE.rvind .AND. hour.EQ.rvind) THEN
+            IF (hour.EQ.rvind) THEN
                hour = value
             END IF
          ELSE IF (desc.EQ.4005) THEN ! Minute
-            IF (value.NE.rvind .AND. minute.EQ.rvind) THEN
+            IF (minute.EQ.rvind) THEN
-               minute = value
-            END IF
-         ELSE IF (desc.EQ.4005) THEN ! Minute
-            IF (value.NE.rvind .AND. minute.EQ.rvind) THEN
                minute = value
             END IF
          ELSE IF (desc.EQ.4006) THEN ! Second
-            IF (value.NE.rvind .AND. second.EQ.rvind) THEN
+            IF (second.EQ.rvind) THEN
                second = value
             END IF
-         ELSE IF (desc.EQ.5001) THEN ! Latitude (high accuracy)
+         ELSE IF (desc.EQ.5001.OR.  ! Latitude (high accuracy)
-            IF (value.NE.rvind .AND. latitude.EQ.rvind) THEN
+     +           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 (value.NE.rvind .AND. latitude.EQ.rvind) THEN
+     +           desc.EQ.6002) THEN ! Longitude (coarse accuracy)
-               latitude = value
+            IF (longitude.EQ.rvind) THEN
-            END IF
-         ELSE IF (desc.EQ.6001) THEN ! Longitude (high accuracy)
-            IF (value.NE.rvind .AND. 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 (value.NE.rvind .AND. longitude.EQ.rvind) THEN
+     +           desc.EQ.7002.OR.    ! Height or altitude
-               longitude = value
+     +           desc.EQ.10070) THEN ! Indicated aircraft altitude
-            END IF
+            IF (flight_level.EQ.rvind) THEN
-         ELSE IF (desc.EQ.7010) THEN ! Flight level
-            IF (value.NE.rvind .AND. flight_level.EQ.rvind) THEN
-               flight_level = value
-            END IF
-         ELSE IF (desc.EQ.7002) THEN ! Height or altitude
-            IF (value.NE.rvind .AND. flight_level.EQ.rvind) THEN
-               flight_level = value
-            END IF
-         ELSE IF (desc.EQ.10070) THEN ! Indicated aircraft altitude
-            IF (value.NE.rvind .AND. flight_level.EQ.rvind) THEN
                flight_level = value
             END IF
-         ELSE IF (desc.EQ.8009) THEN ! Detailed phase of aircraft flight
+         ELSE IF (desc.EQ.8009.OR.  ! Detailed phase of aircraft flight
-            IF (value.NE.rvind .AND. phase.EQ.rvind) THEN
+     +           desc.EQ.8004) THEN ! Phase of aircraft flight
+            IF (phase.EQ.rvind) THEN
                phase = value
                cphase = phase_8009(NINT(phase),missing3)
             END IF
          ELSE IF (desc.EQ.8004) THEN ! Phase of aircraft flight
-            IF (value.NE.rvind .AND. phase.EQ.rvind) THEN
+            IF (phase.EQ.rvind) THEN
                phase = value
                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)
-            IF (value.NE.rvind .AND. TT.EQ.rvind) THEN
+     +          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 (value.NE.rvind .AND. TT.EQ.rvind) THEN
+     +           desc.EQ.11001) THEN ! Wind direction
-               TT = value
+            IF (DD.EQ.rvind) THEN
-            END IF
-         ELSE IF (desc.EQ.12101) THEN ! Temperature/dry bulb temperature (16 bits)
-            IF (value.NE.rvind .AND. TT.EQ.rvind) THEN
-               TT = value
-            END IF
-         ELSE IF (desc.EQ.12001) THEN ! Temperature/dry bulb temperature (12 bits)
-            IF (value.NE.rvind .AND. TT.EQ.rvind) THEN
-               TT = value
-            END IF
-         ELSE IF (desc.EQ.11011) THEN ! Wind direction at 10 m
-            IF (value.NE.rvind .AND. 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 (value.NE.rvind .AND. DD.EQ.rvind) THEN
+     +           desc.EQ.11002) THEN ! Wind speed
-               DD = value
+            IF (FF.EQ.rvind) THEN
-            END IF
-         ELSE IF (desc.EQ.11012) THEN ! Wind speed at 10 m
-            IF (value.NE.rvind .AND. FF.EQ.rvind) THEN
-               FF = value
-            END IF
-         ELSE IF (desc.EQ.11002) THEN ! Wind speed
-            IF (value.NE.rvind .AND. FF.EQ.rvind) THEN
                FF = value
             END IF
          ELSE IF (desc.EQ.7004) THEN ! Pressure
-            IF (value.NE.rvind .AND. PP.EQ.rvind) THEN
+            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 2671: / Line 3149: @@
       END IF
-      WRITE(*,*)
       IF (aircraft.NE.missing8) THEN
-         WRITE(*,'(A,A)') 'aircraft=',aircraft
+         WRITE(*,*)
+         WRITE(*,'(A,A)') 'aircraft=',aircraft(1:lenstr(aircraft,1))
       ELSE IF (flight_number.NE.missing8) THEN
-         WRITE(*,'(A,A)') 'aircraft=',flight_number
+         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 2724: / 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 2740: / Line 3224: @@
       ind = dim
       ascii = IACHAR(string(ind:ind))
-      DO WHILE (ascii.EQ.32 .OR. ascii.EQ.0)
+      DO WHILE (ind.GT.1 .AND. (ascii.EQ.32 .OR. ascii.EQ.0))
          ind = ind - 1
          ascii = IACHAR(string(ind:ind))
       END DO
-      IF (ind.GT.0) THEN
+      IF (ind.EQ.1) THEN
+         ascii = IACHAR(string(1:1))
+         IF (ascii.EQ.32 .OR. ascii.EQ.0) THEN
+            ctrim = missing
+         ELSE
+            ctrim = string(1:1)
+         END IF
+      ELSE IF (ind.GT.1) THEN
          ctrim = string(1:ind)
       ELSE
@@ Line 2753: / Line 3244: @@
       END FUNCTION ctrim
 C     -----------------------------------------------------------------
@@ Line 2760: / Line 3251: @@
       INTEGER phase
       CHARACTER*3 missing3,indicator(6)
       DATA indicator/'res','UNS','LVR','LVW','ASC','DES'/
 C                      1     2     3     4     5     6
@@ Line 2772: / Line 3263: @@
       RETURN
       END FUNCTION phase_8004
 C     -----------------------------------------------------------------
@@ Line 2784: / Line 3275: @@
      +     4*'ASC',4*'DES'/
 C            7-10    11-14
       IF (phase.LE.0 .OR. phase.GT.14) THEN
          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 2941: / Line 3410: @@
       significance = '        '
       jj = 8
       IF (btest(vss,18-15)) THEN
          significance(jj:jj) = 'r' ! level determined by regional decision
          jj = jj-1
@@ Line 2975: / Line 3444: @@
       END SUBROUTINE vss_8042
 C     -----------------------------------------------------------------
@@ Line 3011: / Line 3480: @@
 c
       character*(*) text
-      integer       minlen
+      integer       k,l,lt,minlen
 c
       lt=len(text)
@@ Line 3025: / Line 3494: @@
       return
       end
 </code>
 <code fortran comfilter.f>
 C (C) Copyright 2010, met.no