User Tools

Site Tools


bufr.pm:bufrdump

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).

bufrdump: bufrdump.F comfilter.f
	$(FC) $(FCFLAGS) -o bufrdump  $< -L $(LDIR) -lbufr
bufrdump.F
C (C) Copyright 2010-2016 MET Norway
C
C This program is free software; you can redistribute it and/or modify
C it under the terms of the GNU General Public License as published by
C the Free Software Foundation; either version 2 of the License, or
C (at your option) any later version.
C
C This program is distributed in the hope that it will be useful, but
C WITHOUT ANY WARRANTY; without even the implied warranty of
C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
C General Public License for more details.
C
C You should have received a copy of the GNU General Public License
C along with this program; if not, write to the Free Software
C Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
C 02110-1301, USA.
 
 
C Extract BUFR messages from bufr file(s) and print the data as
C 'parameter=value' lines. See usage_verbose for explanation of the
C options allowed.
C
C Usage: bufrdump <bufr file>
C        [--help]
C        [--filter <filter file>]
C        [--obstype] <amdar|ocea|surface|sounding|sounding->
C        [--stop_on_error]
C        [--rectangle x1 y1 x2 y2]
C        or one or more of
C        [--lon1 x1]
C        [--lat1 y1]
C        [--lon2 x2]
C        [--lat2 y2]
C
C Extract BUFR messages from bufr file and print the data as
C 'parameter = value' lines. See subroutine get_arguments for
C explanation of the options and DokIT for more documentation.
C
C Probably you should set
C export PRINT_TABLE_NAMES=false
C export BUFR_TABLES=/usr/local/lib/bufrtables/
C
C Author: P.Sannes, MET Norway
 
 
      PROGRAM BUFRDUMP
 
      IMPLICIT NONE
 
      CHARACTER(LEN=200) bufr_file       ! Bufr file to read from
      LOGICAL data_only                  ! TRUE if only data section is to be printed
      LOGICAL filter                     ! TRUE if observations should be filtered
      LOGICAL rectangle                  ! TRUE if observations are wanted for a rectangle only
      LOGICAL stop_on_error              ! TRUE if program should stop if a libbufr call returns an error
      CHARACTER(LEN=10) obstype
      CHARACTER(LEN=200) filter_file     ! File containing the filter criteria
      CHARACTER(LEN=200) descriptor_file ! File containing the descriptors requested for partial expansion
      INTEGER verbose                    ! Level of verbose output: 0 - 3 (default 1)
      INTEGER nlibbufr_errors            ! Number of errors encountered in libbufr calls
 
C     Dimensions of arrays used by libbufr routines are the same as those used in
C     example program decode_bufr.F in bufr_000280/example/
 
      INTEGER ibuff(512000)               ! Buffer to hold BUFR message
      INTEGER ibflen                      ! Size in BYTES of the array ibuff
      INTEGER ilen                        ! Size in BYTES of the BUFR message read
      INTEGER wlen                        ! Size in WORDS of the BUFR message read
      INTEGER iunit,iret
      INTEGER iloop
      LOGICAL eof
 
      PARAMETER (ibflen=4*512000)
 
      INCLUDE 'comfilter.f'     ! contains variables used when --filter is set
 
      verbose = 1
 
      CALL get_arguments(bufr_file,filter,filter_file,rectangle,
     +     stop_on_error,obstype)
 
      IF (filter) CALL read_filter(filter,filter_file,verbose)
 
C     Open bufr file - for read
      CALL PBOPEN(iunit,bufr_file,'r',iret)
      IF (iret.NE.0) CALL err_msg('PBOPEN failed for ' // bufr_file)
 
      nlibbufr_errors = 0
      iloop = 0
      DO WHILE (.TRUE.)
         iloop = iloop + 1
 
C     Get the next BUFR message
         CALL get_next_bufr_message(iunit,ibuff,ibflen,ilen,
     +        iloop,verbose,eof)
         IF (eof) GOTO 900
 
C     Decode Bufr message
         wlen = ilen/4 + 1   ! shouldn't really add 1 when ilen is divisible with 4,
                             ! but copy how it is done in decode_bufr in libbufr
         CALL BUFRdecode(ibuff,wlen,filter,rectangle,verbose,
     +        stop_on_error,obstype,nlibbufr_errors)
      END DO
 900  CONTINUE
 
      IF (nlibbufr_errors.GT.0) THEN
         WRITE(*,*) '\n\n  NOTE: At least ',nlibbufr_errors,
     +        ' BUFR message(s) skipped due to errors in libbufr calls',
     +        ' Check output for details of the errors encountered!',
     +        ' (Search on "ERROR")'
      END IF
 
      END PROGRAM BUFRDUMP
 
 
C     -----------------------------------------------------------------
 
      SUBROUTINE get_next_bufr_message(iunit,ibuff,ibflen,ilen,
     +        iloop,verbose,eof)
      IMPLICIT none
      INTEGER ibuff(*)
      INTEGER iunit,ibflen,ilen,iloop,iret,verbose
      LOGICAL eof
 
      eof = .FALSE.
      iret = 0
      CALL PBBUFR(iunit,ibuff,ibflen,ilen,iret)
      IF (iret .EQ. -1 ) THEN
         eof = .TRUE.
         RETURN
      ELSEIF (iret .EQ. -2 ) THEN
         CALL err_msg('ERROR: File handling problem')
      ELSEIF (iret .EQ. -3 ) THEN
         CALL err_msg('ERROR: Array too small')
      END IF
      IF (verbose .GT. 2) THEN
         WRITE (*,'(1X,A,I10)') 'BUFR message number= ',iloop
         WRITE (*,'(1X,A,I10)') 'Length of BUFR message = ',ilen
      END IF
      END SUBROUTINE get_next_bufr_message
 
C     -----------------------------------------------------------------
 
      SUBROUTINE BUFRdecode(ibuff,ilen,filter,rectangle,
     +     verbose,stop_on_error,obstype,nlibbufr_errors)
      IMPLICIT NONE
      INTEGER ibuff(*)
      INTEGER ilen              ! size in WORDS (not bytes) of the BUFR message read
      LOGICAL filter
      LOGICAL rectangle
      LOGICAL stop_on_error
      CHARACTER(LEN=*) obstype
      INTEGER nlibbufr_errors   ! Number of errors encountered in libbufr calls
      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
 
C     Dimensions of arrays used by libbufr routines are the same as those used in
C     example program decode_bufr.F in bufr_000310/example/
      PARAMETER (kelem = 160000,kvals = 4096000)
 
C     BUFREX variables
      INTEGER ksup(9)           ! Integer array containing suplementary information
      INTEGER ksec0(3)          ! Bufr section 0 information
      INTEGER ksec1(40)
      INTEGER ksec2(4096)
      INTEGER ksec3(4)
      INTEGER ksec4(2)
      CHARACTER*64 cnames(kelem) ! Bufr Table B element names
      CHARACTER*24 cunits(kelem) ! Bufr Table B units
      CHARACTER*80 cvals(kelem) ! Bufr code table or CCITTIA5 Bufr elements entries
                                ! Using kelem as dimension seems somewhat arbitrary,
                                ! - just copying from decode_bufr.F in libbufr
      REAL*8 values(kvals)      ! expanded data values
 
C     BUSEL variables
      INTEGER ktdlst(2*kelem)     ! array containing data descriptors in section 3
      INTEGER ktdexp(2*kelem)     ! array containing expanded data descriptors
                                  ! Use 2*kelem, not kelem; copied from decode_bufr.F in libbufr
      INTEGER ktdlen            ! number of data descriptors in section 3
      INTEGER ktdexl            ! number of entries in list of expanded data descriptors
      INTEGER kerr              ! returned error code - but note that some libbufr
                                ! routines will return immediately if this is not 0 on input
      INTEGER ksub,ksub1,ksub2
 
C     BUSQR variables
      INTEGER kreq(2)
      INTEGER krql
      INTEGER krq(kelem)
      REAL*8 rqv(kelem)
      INTEGER rqd(101)           ! Requested descriptors
 
C     BUUKEY variables
      INTEGER key(46)
 
      INTEGER found_list(10000),num_found
 
      INTEGER i
 
      kerr = 0
      metar = .FALSE.
 
C     Do a partial expansion without quality control
      kreq(1) = 1               ! All original elements without quality control
      kreq(2) = 0
      krql = 0
 
C     Partial expansion
      CALL BUSRQ(kreq,krql,rqd,rqv,kerr)
      IF (kerr.NE.0) THEN
         IF (stop_on_error) STOP 1
         WRITE(*,*) 'ERROR IN BUSRQ: KERR= ',kerr
         nlibbufr_errors = nlibbufr_errors + 1
         RETURN
      END IF
C
C     Decode BUFR message into fully decoded form.
C
 
C     Using parameter kelem in call to BUFREX might be too big for some
C     multisubset messages (or too small for other messages). Have
C     copied the method used in decode_bufr.F in libbufr, first calling
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
         IF (stop_on_error) STOP 1
         WRITE(*,*) 'ERROR IN BUS012: KERR= ',kerr
         nlibbufr_errors = nlibbufr_errors + 1
         RETURN
      END IF
      kxelem = kvals/ksup(6)
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,
     +     ksec3,ksec4,kxelem,cnames,cunits,kvals,values,
     +     cvals,kerr)
      IF (kerr.EQ.45) RETURN      ! No requested elements in data
      IF (kerr.NE.0) THEN
         IF (stop_on_error) STOP 1
         WRITE(*,*) 'ERROR IN BUFREX: KERR= ',kerr
         nlibbufr_errors = nlibbufr_errors + 1
         RETURN
      END IF
      IF (verbose .GT. 2) WRITE (*,'(1X,A,I10)')
     +     'Number of subsets:',ksup(6)
 
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) 6 = Radar data 31 = Oceanographic data
      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
C     Return list of Data Descriptors from Section 3 of Bufr message,
C     and total/requested list of elements. BUFREX must have been called
C     before BUSEL. If filter is set, BUSEL2 will later be called to be
C     used for printing each subset, but for filter we need to know the
C     first members of ktdlst and ktdexl already now, to be able to call
C     find_stations.
         CALL BUSEL(ktdlen,ktdlst,ktdexl,ktdexp,kerr)
         IF (kerr .GT. 0) THEN
            IF (stop_on_error) STOP 1
            WRITE(*,*) 'ERROR IN BUSEL: KERR= ',kerr
            nlibbufr_errors = nlibbufr_errors + 1
            RETURN
         END IF
      END IF
 
      IF (filter) THEN
         CALL find_stations(ktdexp,ktdexl,values,cvals,ksup(6),
     +         kxelem,found_list,num_found)
      ENDIF
      IF (filter .AND. num_found.GT.0) THEN
         DO i=1,num_found
            ksub = found_list(i)
            CALL BUSEL2(ksub,kxelem,ktdlen,ktdlst,ktdexl,ktdexp,
     +           cnames,cunits,kerr)
            IF (kerr.NE.0) THEN
               IF (stop_on_error) STOP 1
               WRITE(*,*) 'ERROR IN BUSEL2: KERR= ',kerr
               nlibbufr_errors = nlibbufr_errors + 1
               RETURN
            END IF
 
            IF (obstype.EQ.'surface') THEN
               CALL print_surface_values(ksub,kxelem,ktdexl,ktdexp,
     +              values,cvals,rectangle,metar,verbose)
            ELSE IF (obstype(1:8).EQ.'sounding') THEN
               CALL print_sounding_values(ksub,kxelem,ktdexl,ktdexp,
     +              values,cvals,rectangle,obstype,verbose)
            ELSE IF (obstype.EQ.'amdar') THEN
               CALL print_amdar_values(ksub,kxelem,ktdexl,ktdexp,
     +              values,cvals,rectangle,verbose)
            ELSE IF (obstype.EQ.'ocea') THEN
               CALL print_oceanographic_values(ksub,kxelem,ktdexl,
     +              ktdexp,values,cvals,rectangle,verbose)
            ELSE IF (ksec1(6).LE.1) THEN
               CALL print_surface_values(ksub,kxelem,ktdexl,ktdexp,
     +              values,cvals,rectangle,metar,verbose)
            ELSE IF (ksec1(6).EQ.2) THEN
               CALL print_sounding_values(ksub,kxelem,ktdexl,ktdexp,
     +              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,
     +              ktdexp,values,cvals,rectangle,verbose)
            END IF
         ENDDO
      ELSE IF (.NOT. filter) THEN
         DO ksub=1,ksup(6)
            CALL BUSEL2(ksub,kxelem,ktdlen,ktdlst,ktdexl,ktdexp,
     +           cnames,cunits,kerr)
            IF (kerr.NE.0) THEN
               IF (stop_on_error) STOP 1
               WRITE(*,*) 'ERROR IN BUSEL2: KERR= ',kerr
               nlibbufr_errors = nlibbufr_errors + 1
               RETURN
            END IF
 
            IF (obstype.EQ.'surface') THEN
               CALL print_surface_values(ksub,kxelem,ktdexl,ktdexp,
     +              values,cvals,rectangle,metar,verbose)
            ELSE IF (obstype(1:8).EQ.'sounding') THEN
               CALL print_sounding_values(ksub,kxelem,ktdexl,ktdexp,
     +              values,cvals,rectangle,obstype,verbose)
            ELSE IF (obstype.EQ.'amdar') THEN
               CALL print_amdar_values(ksub,kxelem,ktdexl,ktdexp,
     +              values,cvals,rectangle,verbose)
            ELSE IF (obstype.EQ.'ocea') THEN
               CALL print_oceanographic_values(ksub,kxelem,ktdexl,
     +              ktdexp,values,cvals,rectangle,verbose)
            ELSE IF (ksec1(6).LE.1) THEN
               CALL print_surface_values(ksub,kxelem,ktdexl,ktdexp,
     +              values,cvals,rectangle,metar,verbose)
            ELSE IF (ksec1(6).EQ.2) THEN
               CALL print_sounding_values(ksub,kxelem,ktdexl,ktdexp,
     +              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
C     bulletins all have data category 1 (and sub-category 25 or 255),
C     so for these we end up calling print_surface_values (unless
C     --obstype ocea is used)
               CALL print_oceanographic_values(ksub,kxelem,ktdexl,
     +              ktdexp,values,cvals,rectangle,verbose)
            END IF
         END DO
      END IF
 
      END SUBROUTINE BUFRdecode
 
C     -----------------------------------------------------------------
 
      SUBROUTINE get_arguments(bufr_file,filter,filter_file,
     +     rectangle,stop_on_error,obstype)
 
      IMPLICIT NONE
 
      CHARACTER(LEN=*) bufr_file ! Bufr file to read from
      LOGICAL data_only          ! TRUE if only data section is to be printed
      LOGICAL filter             ! TRUE if observations should be filtered
      LOGICAL rectangle          ! TRUE if observations are wanted for a rectangle only
      LOGICAL stop_on_error      ! TRUE if program should stop if a libbufr call returns an error
      CHARACTER(LEN=*) filter_file ! File containing the filter criteria
      CHARACTER(LEN=*) obstype
 
      CHARACTER(LEN=200) argument ! Buffer for next argument
      CHARACTER(LEN=1350) Usage
      CHARACTER(LEN=950) Help
      INTEGER iargc,iarg
 
C     Variables used for geographical filtering av observations
      REAL*8 x1,y1,x2,y2
      COMMON /COM_RECTANGLE/  x1,y1,x2,y2
      x1 = -180.00001
      y1 = -90.00001
      x2 = 180.00001
      y2 = 90.00001
 
      Usage = '\nUsage: bufrdump <bufr file> [options]\n\n'
     +     // 'Will print section 4 in BUFR messages in <bufr file>'
     +     // ' as "parameter=value" lines.'
     +     // '\nOptions (may be abbreviated to one character)'
     +     // ' are:\n'
     +     // '\t--help\t\tPrint this Usage plus some more help\n'
     +     // '\t--filter <filter file>\n'
     +     // '\t\t\tDecode observations meeting criteria in '
     +     // '<filter file> only\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. '
     +     // '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 '
     +     // 'left coordinates\n\t\t\t(x1,y1) (i.e.lon,lat) '
     +     // 'and upper right coordinates (x2,y2) only\n'
     +     // 'As an alternative to --rectangle you can use one or '
     +     // 'more of\n'
     +     // '\t--lon1 x1'
     +     // '\tDecode observations with longitude >= x1 only\n'
     +     // '\t--lat1 y1'
     +     // '\tDecode observations with latitude >= y1 only\n'
     +     // '\t--lon2 x2'
     +     // '\tDecode observations with longitude <= x2 only\n'
     +     // '\t--lat2 y2'
     +     // '\tDecode observations with latitude <= y2 only\n'
     +     // '\t\t\tx1,y1,x2,y2 should be degrees or decimal degrees'
 
      Help = '\nYou should probably set\n'
     +     // '\texport BUFR_TABLES=/usr/local/lib/bufrtables/\n'
     +     // '\texport PRINT_TABLE_NAMES=false\n'
     +     // '\nUsing --filter will decode only those observations '
     +     // 'that meet the criteria in\n<filter file>. An example '
     +     // 'of a filter file is\n\n'
     +     // 'D: 001001 I2.2\n'
     +     // '01\n'
     +     // 'D: 001001 I2.2 001002 I3.3\n'
     +     // '03 895\n'
     +     // '06 252\n'
     +     // 'D: 001011 A9\n'
     +     // 'LMEL     \n\n'
     +     // 'which decodes all observations with block number 01, '
     +     // 'two other specific wmo\nstations and one ship. '
     +     // '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 '
     +     // 'is corrupt) the message is skipped, and the number '
     +     // 'of\nerrors is reported at end of output - unless '
     +     // 'option --stop_on_error has been\napplied (note that '
     +     // 'some error messages from libbufr might still have been '
     +     // 'sent\nto STDOUT).\n'
 
C     Default values
      filter = .FALSE.
      obstype = ''
      rectangle = .FALSE.
      stop_on_error = .FALSE.
      bufr_file = ''
      iarg = 0
 
      IF (iargc() .LT. 1) THEN
         WRITE(*,*) Usage
         CALL EXIT(0)
      END IF
 
C     Loop: reading of command options
 100  CONTINUE
      iarg = iarg + 1
      CALL getarg(iarg,argument)
      IF (argument .EQ. '--help' .OR.
     +     argument .EQ. '--h' .OR. argument .EQ. '-h') THEN
         WRITE(*,*) Usage
         WRITE(*,*) Help
         CALL EXIT(0)
      ELSE IF (argument .EQ. '--filter' .OR.
     +     argument .EQ. '--f' .OR. argument .EQ. '-f') THEN
         filter = .TRUE.
         iarg = iarg + 1
         CALL getarg(iarg,filter_file)
      ELSE IF (argument .EQ. '--obstype' .OR.
     +     argument .EQ. '--o' .OR. argument .EQ. '-o') THEN
         iarg = iarg + 1
         CALL getarg(iarg,obstype)
      ELSE IF (argument .EQ. '--stop_on_error') THEN
         stop_on_error = .TRUE.
      ELSE IF (argument .EQ. '--rectangle' .OR.
     +     argument .EQ. '--r' .OR. argument .EQ. '-r') THEN
         rectangle = .TRUE.
         iarg = iarg + 1
         CALL getarg(iarg,argument)
         READ (argument,*) x1
         iarg = iarg + 1
         CALL getarg(iarg,argument)
         READ (argument,*) y1
         iarg = iarg + 1
         CALL getarg(iarg,argument)
         READ (argument,*) x2
         iarg = iarg + 1
         CALL getarg(iarg,argument)
         READ (argument,*) y2
      ELSE IF (argument .EQ. '--lon1') THEN
         rectangle = .TRUE.
         iarg = iarg + 1
         CALL getarg(iarg,argument)
         READ (argument,*) x1
      ELSE IF (argument .EQ. '--lat1') THEN
         rectangle = .TRUE.
         iarg = iarg + 1
         CALL getarg(iarg,argument)
         READ (argument,*) y1
      ELSE IF (argument .EQ. '--lon2') THEN
         rectangle = .TRUE.
         iarg = iarg + 1
         CALL getarg(iarg,argument)
         READ (argument,*) x2
      ELSE IF (argument .EQ. '--lat2') THEN
         rectangle = .TRUE.
         iarg = iarg + 1
         CALL getarg(iarg,argument)
         READ (argument,*) y2
      ELSE IF (argument(1:2) .EQ. '--' .OR.
     +     argument .EQ. '--d' .OR. argument .EQ. '-d') THEN
         WRITE(*,*) 'Unknown option ',argument
         WRITE(*,*) Usage
         CALL EXIT(1)
      ELSE
C     Read in Bufr file to decode
         bufr_file = argument
      END IF
      IF (iarg .LT. iargc()) GOTO 100
 
      IF (bufr_file .EQ. '') THEN
         WRITE(*,*) 'No BUFR file given'
         CALL EXIT(1)
      END IF
 
      IF (obstype.NE.'' .AND. obstype.NE.'amdar' .AND. obstype.NE.'ocea'
     +     .AND. obstype.NE.'surface' .AND. obstype.NE.'sounding'
     +     .AND. obstype.NE.'sounding-') THEN
         WRITE(*,*) 'Argument to --obstype must be one of amdar, ocea,',
     +   ' surface, sounding or sounding-'
         CALL EXIT(1)
      END IF
 
      END SUBROUTINE get_arguments
 
C     -----------------------------------------------------------------
 
      SUBROUTINE err_msg(msg)
      IMPLICIT NONE
      CHARACTER(LEN=*) msg
      WRITE(*,*) msg
      STOP
      END SUBROUTINE err_msg
 
C     -----------------------------------------------------------------
 
      SUBROUTINE err_msg1(msg,kerr)
      IMPLICIT NONE
      CHARACTER(LEN=*) msg
      INTEGER kerr
      WRITE(*,*) msg,kerr
      STOP
      END SUBROUTINE err_msg1
 
C     -----------------------------------------------------------------
 
      SUBROUTINE read_filter(filter,filter_file,verbose)
C
C     Reads the content of filter file into matrices fid and fidformat,
C     arrays nd_fid and nvl_fid for the descriptors, into matrices fivI
C     and fivC for the values to filter on. Also sets nfidlines and
C     nfidvalues. See comfilter.f for explanation of variables. See Help
C     in get_arguments for an example of filter file.
C
      IMPLICIT NONE
      LOGICAL filter            ! TRUE on input. Might be changed to FALSE on output
                                ! if no filter condition is found in filter file
      CHARACTER(LEN=200) filter_file
      INTEGER verbose,ios,inplen,i1,i2,i3,ifid,jfid,lenstr,fmt_len
      CHARACTER*100 inpline,readformat
      CHARACTER*10 fmt
      LOGICAL new_fid
 
      INCLUDE 'comfilter.f'
 
      OPEN (UNIT = 11,FILE = filter_file,STATUS = 'OLD',
     +     IOSTAT = ios)
      IF (ios.NE.0) CALL err_msg('Cannot open ' // filter_file)
 
      nfidlines = 0
      nfivlines = 0
      ifid = 0
      filter_loop: DO WHILE (.TRUE.)
C     Note: a line consisting of blanks or the new line character only
C     marks end of filter conditions
         READ (11,'(A100)',END=11) inpline
         IF (verbose.GT.2) WRITE (*,*) inpline
         inplen = lenstr(inpline,1)
         IF (inplen.EQ.0) EXIT filter_loop
         i1 = 0
         i2 = 0
         CALL advance_word(inpline,inplen,i1,i2)
         IF (i1 .EQ. 0) EXIT filter_loop
         IF (inpline(i1:i1).EQ.'#') CYCLE filter_loop ! Comment line
 
         IF (inpline(i1:i1+1).EQ.'D:') THEN
C     Descriptor line (e.g. D: 001001 I2.2 001002 I3.3)
            ifid = ifid + 1
            IF (ifid.GT.dim1_fid) CALL err_msg(
     +      'bufrread: dimension 1 limit exeeded for array fid')
            jfid = 0
            DO WHILE (i1 .NE. 0)
               CALL advance_word(inpline,inplen,i1,i2)
               IF (i1 .NE. 0) THEN
                  jfid = jfid + 1
                  IF (jfid.GT.dim1_fid) CALL err_msg(
     +                 'bufrread: dimension 2 exeeded for array fid')
C     Read descriptor
                  READ (inpline(i1:i2),*) fid(ifid,jfid)
                  CALL advance_word(inpline,inplen,i1,i2)
                  IF (i1 .EQ. 0)
     +                 CALL err_msg('Wrong format in filter file')
C     Read format
                  READ (inpline(i1:i2),*) fidformat(ifid,jfid)
               END IF
            END DO
            nd_fid(ifid) = jfid
            nfidlines = ifid
         ELSE
 
C     Value line (e.g. 03 895 or: North Shields)
            nfivlines = nfivlines + 1
            IF (nfivlines.GT.dim_fiv) CALL err_msg(
     +      'bufrread: dimension limit exeeded for array fivI/C')
            nvl_fid(ifid) =  nvl_fid(ifid) + 1
C     Read in value field, descriptor for descriptor:
            DO jfid = 1,nd_fid(ifid)
               fmt = fidformat(ifid,jfid)
               IF (fmt(1:1).EQ.'I') THEN
                  READ (inpline(i1:i2),*) fivI(nfivlines,jfid)
                  IF (verbose.GT.2) WRITE(*,*) 'nfivlines,jfid,fivI=',
     +                 nfivlines,jfid,fivI(nfivlines,jfid)
               ELSE IF (fmt(1:1).EQ.'A') THEN
C                 String may include space, therefore i2 may need to be adjusted
                  READ (fmt(2:lenstr(fmt,0)),*) fmt_len
                  i2 = i1 + fmt_len
                  READ (inpline(i1:i2),*) fivC(nfivlines,jfid)
                  IF (verbose.GT.2) WRITE(*,*) 'nfivlines,jfid,fivC=',
     +                 nfivlines,jfid,' ',fivC(nfivlines,jfid)
               ELSE
                  CALL err_msg('Wrong value format in filter file')
               END IF
               CALL advance_word(inpline,inplen,i1,i2)
            END DO
         ENDIF
      END DO filter_loop
 
 11   CONTINUE
      CLOSE(11)
 
      IF (verbose.GT.2) WRITE(*,*) 'nfidlines=',nfidlines
      IF (nfivlines.EQ.0) filter = .FALSE.
 
      END SUBROUTINE read_filter
 
C     -----------------------------------------------------------------
 
      subroutine advance_word(str,length,beginp,endp)
C
C     The character variable 'str' contains words of non-blank
C     characters separated by blanks. 'beginp' and 'endp' points to the
C     beginning and end of a word (char indices within 'str'). This
C     subroutine advances these pointers to point to the next word
C     within 'str'. 'length' is the length of 'str'. Initially, 'beginp'
C     and 'endp' is set to 0 (outside this subroutine). The first call
C     to 'advance_word' will then set these pointers to the first word
C     within 'str'. When no more words are found, 'beginp' and 'endp'
C     are again set to 0.
C
      implicit none
      character*(*) str
      integer length,beginp,endp
      integer i1,i2
      i1 = endp+1
      do while (i1 .le. length)
         if (str(i1:i1) .ne. ' ') then
            beginp = i1
            i2 = index(str(i1:length),' ')
            if (i2 .eq. 0) then
               endp = length
               return
            else
               endp = i1 + i2 - 2
               return
            end if
         else
            i1 = i1 + 1
         end if
      end do
      beginp = 0
      endp = 0
      end subroutine advance_word
 
C     -----------------------------------------------------------------
 
      SUBROUTINE find_stations(ktdexp,ktdexl,values,cvals,num_subsets,
     +        kxelem,found_list,num_found)
C
C     Return the number of subsets in BUFR message which contain an
C     observation matching one of the criteria in filter file as
C     num_found, the subset numbers in found(list(i)), i=1,...,num_found
C
      IMPLICIT NONE
      INTEGER found_list(10000),num_found
      INTEGER ktdexp(*)         ! array containing expanded data descriptors
      INTEGER ktdexl            ! number of entries in list of expanded data descriptors
      INTEGER kxelem            ! number of data elements for each section 4 report
      CHARACTER*80 cvals(*)     ! Bufr code table or CCITTIA5 Bufr elements entries
      REAL*8 values(*)          ! expanded data values
      INTEGER num_subsets
 
      INTEGER i1,i2,nsub,ifvl,itd,iv,icv,ifiv
      LOGICAL match             ! Set to true if a value in data matches value in
                                ! filter, for a single descriptor
      LOGICAL int_fid           ! Set to true if filter descriptor is integer type
      LOGICAL char_fid          ! Set to true if filter descriptor is character type
 
      INCLUDE 'comfilter.f'
 
      REAL*8 rvind              ! missing value for real*8 data
      PARAMETER (rvind=1.7E38)
 
      num_found = 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 value lines (for given) condition:
            DO ifvl = 1,nvl_fid(i1)
               ifiv = ifiv + 1
C     loop through all descriptors in condition:
               DO i2 = 1,nd_fid(i1)
                  match = .FALSE.
                  int_fid = .FALSE.
                  char_fid = .FALSE.
                  IF (fidformat(i1,i2)(1:1).EQ.'I') THEN
                     int_fid = .TRUE.
                  ELSE IF (fidformat(i1,i2)(1:1).EQ.'A') THEN
                     char_fid = .TRUE.
                  END IF
C     loop through all data in subset:
                  DO itd = 1,ktdexl
                     IF (ktdexp(itd).EQ.fid(i1,i2)) THEN
                        iv = (nsub-1)*kxelem + itd
                        IF (values(iv).EQ.rvind)
     +                       GOTO 400 ! next condition, no point in checking
                                      ! other filter value lines
                        IF (int_fid) THEN
                           IF (nint(values(iv)).EQ.fivI(ifiv,i2)) THEN
                              match = .TRUE.
                              EXIT ! exit innermost loop
                           END IF
                        ELSE IF (char_fid) THEN
                           icv = nint(values(iv))/1000
                           IF (cvals(icv).EQ.fivC(ifiv,i2)) THEN
                              match = .TRUE.
                              EXIT ! exit innermost loop
                           END IF
                        END IF
                     END IF
                   END DO       ! all data in subset
                   IF (.NOT.match) ! there is one descriptor in condition which
     +                  EXIT       ! subset does not match, so go to next
                                   ! 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 DO              ! all filter value lines (for given) condition
 400        CONTINUE
         END DO                 ! all different conditions
 500     CONTINUE
      END DO                    ! all subsets
 
      END SUBROUTINE find_stations
 
C     -----------------------------------------------------------------
 
C     Identify pressure, temperature etc and print parameter=value to
C     screen. Note: have seen high resolution data with 7002 HEIGHT OR
C     ALTITUDE and wind profiler data with 007009 GEOPOTENTIAL HEIGHT
C     used as vertical coordinate instead of 7004 PRESSURE (10004 was
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,obstype,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
      CHARACTER(LEN=*) obstype  ! Input: influence flow if set to 'sounding-'
      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,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,i,ind,ind2,ind3
      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
      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
 
      IF (obstype.EQ.'sounding-') THEN
         reduce = .TRUE.
      ELSE
         reduce = .FALSE.
      END IF
 
      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
      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
      year = rvind
      month = rvind
      day = rvind
      hour = rvind
      minute = rvind
      latitude = rvind
      longitude = rvind
      height = rvind
      DO n=1,maxlevel
         P(n) = rvind
         D(n) = rvind
         F(n) = rvind
         T(n) = rvind
         TD(n) = rvind
         h(n) = rvind
         la(n) = rvind
         lo(n) = rvind
         tp(n) = rvind
         vss(n) = vss_missing
         wsb(n) = rvind
         wsa(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.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
         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.OR.  ! Height of station
     +           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
C     All following descriptors come after pressure in 309052 (or Hirlam template)
         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.12101) THEN ! Temperature/dry bulb temperature (16 bits)
               T(n) = value
            ELSE IF (desc.EQ.12001) THEN ! Temperature/dry bulb temperature (12 bits)
               T(n) = value
            ELSE IF (desc.EQ.12103) THEN ! Dew-point temperature (16 bits)
               TD(n) = value
            ELSE IF (desc.EQ.12003) THEN ! Dew-point temperature (12 bits)
               TD(n) = value
            ELSE IF (desc.EQ.10003.AND.value.NE.rvind) THEN ! Geopotential
               h(n) = value/9.81
            ELSE IF (desc.EQ.10009) THEN ! Geopotential height
               h(n) = value
            ELSE IF (desc.EQ.5015) THEN ! Latitude displacement
               la(n) = value
            ELSE IF (desc.EQ.6015) THEN ! Longitude displacement
               lo(n) = value
            ELSE IF (desc.EQ.8001 .AND. value.NE.rvind) THEN ! Vertical sounding significance
C     In HIRLAM template (309007) 008001 comes AFTER 7004 pressure
               CALL vss_8001(NINT(value),vss(n))
            ELSE IF (desc.EQ.11061) THEN ! Absolute wind shear in 1 km layer below
               wsb(n) = value
            ELSE IF (desc.EQ.11062) THEN ! Absolute wind shear in 1 km layer above
               wsa(n) = value
            END IF
         END IF
      END DO
 110  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 (wigos_series.NE.rvind .AND. wigos_issuer.NE.missing5
     +        .AND. wigos_issueno.NE.missing5
     +        .AND. wigos_localid.NE.missing16) THEN
         ind = index(wigos_issuer,' ') - 1
         IF (ind.EQ.-1) ind = 5
         ind2 = index(wigos_issueno,' ') - 1
         IF (ind2.EQ.-1) ind2 = 5
         ind3 = index(wigos_localid,' ') - 1
         IF (ind3.EQ.-1) ind3 = 16
         WRITE(*,*)
         WRITE(*,'(A,I1.1,A1,A,A1,A,A1,A)')
     +        'wigosid=',NINT(wigos_series),
     +        '-',wigos_issuer(1:ind),
     +        '-',wigos_issueno(1:ind2),
     +        '-',wigos_localid(1:ind3)
      ELSE IF (call_sign.NE.missing9) THEN
         WRITE(*,*)
         WRITE(*,'(A,A)') 'call_sign=',call_sign(1:lenstr(call_sign,1))
      ELSE IF (flight_number.NE.missing8) THEN
         WRITE(*,*)
         WRITE(*,'(A,A)') 'aircraft=',
     +        flight_number(1:lenstr(flight_number,1))
      ELSE
         IF (verbose .GT. 1) THEN
            WRITE(*,*)
            WRITE(*,*) 'Both wmonr, wigosid, call sign and aircraft'
     +           // ' flight number 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)
 900     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
         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
         IF (tp(n).NE.rvind) THEN
            WRITE(*,'(A,I12)') 't=',NINT(tp(n))
         END IF
         IF (la(n).NE.rvind) THEN
            WRITE(*,'(A,F11.5)') 'la=',la(n)
         END IF
         IF (lo(n).NE.rvind) THEN
            WRITE(*,'(A,F11.5)') 'lo=',lo(n)
         END IF
C     Choose to display vss even when it is empty
         WRITE(*,'(A,A10)') 'vss=',vss(n)
         IF (P(n).NE.rvind) THEN
            WRITE(*,'(A,F11.1)') 'PP=',P(n)/100 ! hPa
         END IF
         IF (h(n).NE.rvind) THEN
            WRITE(*,'(A,I11)') 'gh=',NINT(h(n))
         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 (T(n).NE.rvind) THEN
            WRITE(*,'(A,F11.1)') 'TT=',T(n)-273.15
         END IF
         IF (TD(n).NE.rvind) THEN
            WRITE(*,'(A,F11.1)') 'TD=',TD(n)-273.15
         END IF
         IF (wsb(n).NE.rvind) THEN
            WRITE(*,'(A,F10.1)') 'wsb=',wsb(n)
         END IF
         IF (wsa(n).NE.rvind) THEN
            WRITE(*,'(A,F10.1)') 'wsa=',wsa(n)
         END IF
      END DO
 
      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
 110  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)
 900     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,metar,verbose)
C     Identify pressure, temperature etc and print parameter=value to screen
      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
      LOGICAL metar             ! Input: TRUE if metar (data subcategory 5)
      INTEGER verbose           ! Input: verbose level
 
      REAL*8 rvind              ! missing value for real*8 data
      PARAMETER (rvind=1.7E38)
 
      CHARACTER*5 wigos_issuer,wigos_issueno,missing5
      CHARACTER*8 icao_id,missing8
      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,
     +     EE,ES,EV_1,EV_24,FF,FG,FG_010,FG_1,FG_X,FX,FX_1,FX_X,HL,
     +     HW,HW1,HW2,HWA,NH,NN,OT_1,OT_24,PH,PO,PP,PR,PW,PW1,PW2,PWA,
     +     QD,QE,QK,QL,QO,QS,QD_24,QE_24,QK_24,QL_24,QO_24,QS_24,
     +     RR_1,RR_3,RR_6,RR_12,RR_24,RS,SA,SG,SI,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,a3,buoy_id5,buoy_id7,ds,
     +     height,hhh,hp,hour_p,II,iii,ix,latitude,longitude,
     +     minute_p,TbTbTb,vert_sign_first,vs,wmo_region_number,
     +     wmo_region_subarea,state_id,national_number,
     +     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).
                                ! Not used for metar
      INTEGER num_cloud_layers  ! Number of individual cloud layers, set
                                ! to value of 031001 (delayed descriptor
                                ! 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
                                ! data will then be attempted decoded.
      INTEGER cloud_layer       ! Numbers the individual cloud layers
      LOGICAL surface_data      ! Set to false if 007062 ('Depth below sea/water
                                ! 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,ind2,ind3
      INTEGER degree2dir,NNtoWMO_N
 
C     Variables used for geographical filtering av observations
C     (--rectangle option set)
      REAL*8 x1,y1,x2,y2
      COMMON /COM_RECTANGLE/  x1,y1,x2,y2
 
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
      missing32 = missing16 // missing16
 
      cloud_type_count = 0
      bad_cloud_data = .FALSE.
      num_cloud_layers = -1
      surface_data = .TRUE.
      time_of_last_position = .FALSE.
 
C     Initialize all parameters to missing values
      hour_p = rvind
      minute_p = rvind
      icao_id = missing8
      call_sign = missing9
      wigos_issuer = missing5
      wigos_issueno = missing5
      wigos_localid = missing16
      name = missing20
      long_name = missing32
      AA = rvind
      BI = rvind
      CH = rvind
      CI = rvind
      CL = rvind
      CM = rvind
      DD = rvind
      DG = rvind
      DG_010 = rvind
      DG_1 = rvind
      DG_X = rvind
      DI = rvind
      DW1 = rvind
      DW2 = rvind
      EE = rvind
      ES = rvind
      EV_1 = rvind
      EV_24 = rvind
      FF = rvind
      FG = rvind
      FG_010 = rvind
      FG_1 = rvind
      FG_X = rvind
      FX = rvind
      FX_1 = rvind
      FX_X = rvind
      HL = rvind
      HW = rvind
      HW1 = rvind
      HW2 = rvind
      HWA = rvind
      NH = rvind
      NN = rvind
      PO = rvind
      OT_1 = rvind
      OT_24 = rvind
      PH = rvind
      PP = rvind
      PR = rvind
      PR = rvind
      PW = rvind
      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
      RR_6 = rvind
      RR_12 = rvind
      RR_24 = rvind
      RS = rvind
      SA = rvind
      SG = rvind
      SI = rvind
      SS_24 = rvind
      TA = rvind
      TAN_12 = rvind
      TAX_12 = rvind
      TAN = rvind
      TAX = rvind
      TD = rvind
      TGN_12 = rvind
      TW = rvind
      UU = rvind
      VV = rvind
      W1 = rvind
      W2 = rvind
      WW = rvind
      XIS = rvind
      ZI = rvind
      year = rvind
      month = rvind
      day = rvind
      hour = rvind
      minute = rvind
      latitude = rvind
      longitude = rvind
      height = rvind
      hp = rvind
      vert_sign_first = rvind
      II = rvind
      iii = rvind
      ix = rvind
      a3 = rvind
      hhh = rvind
      ds = rvind
      vs = rvind
      TbTbTb = 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
         vert_sign(i) = rvind
         CC(i) = rvind
         HS(i) = rvind
         NS(i) = rvind
      END DO
 
C     Loop through all expanded descriptors
      DO idx=1,ktdexl
         desc = ktdexp(idx)
         value = values(idx + (ksub-1)*kxelem)
 
         IF (ABS(value - rvind)/rvind.LE.0.001) THEN
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.8021) THEN ! Time significance
            IF (NINT(value).EQ.26) THEN
               time_of_last_position = .TRUE.
            ELSE
               time_of_last_position = .FALSE.
            END IF
         ELSE IF (desc.EQ.1001 .AND. II.EQ.rvind) THEN ! WMO block number
            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
            cidx = int(value/1000)
            name = cvals(cidx)  ! CCITTIA5 data
            name = ctrim(name,20,missing20)
         ELSE IF (desc.EQ.1019) THEN  ! Long station or site name
            cidx = int(value/1000)
            long_name = cvals(cidx) ! CCITTIA5 data
            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 (ix.EQ.rvind) THEN
               ix = value
            END IF
         ELSE IF (desc.EQ.4001) THEN ! Year
            IF (year.EQ.rvind.AND..NOT.time_of_last_position) THEN
               year = value
            END IF
         ELSE IF (desc.EQ.4002) THEN ! Month
            IF (month.EQ.rvind.AND..NOT.time_of_last_position) THEN
               month = value
            END IF
         ELSE IF (desc.EQ.4003) THEN ! Day
            IF (day.EQ.rvind.AND..NOT.time_of_last_position) THEN
               day = value
            END IF
         ELSE IF (desc.EQ.4004) THEN ! Hour
            IF (hour.EQ.rvind.AND..NOT.time_of_last_position) THEN
               hour = value
            END IF
         ELSE IF (desc.EQ.4005) THEN ! Minute
            IF (minute.EQ.rvind.AND..NOT.time_of_last_position) THEN
               minute = value
            END IF
         ELSE IF (desc.EQ.5001.OR.  ! Latitude (high accuracy)
     +           desc.EQ.5002) THEN ! Latitude (coarse accuracy)
            IF (latitude.EQ.rvind) THEN
               latitude = value
            END IF
         ELSE IF (desc.EQ.6001.OR.  ! Longitude (high accuracy)
     +           desc.EQ.6002) THEN ! Longitude (coarse accuracy)
            IF (longitude.EQ.rvind) THEN
               longitude = value
            END IF
         ELSE IF (desc.EQ.7001.OR.  ! Height of station
     +           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.7031) THEN ! Hp
            IF (hp.EQ.rvind) THEN
               hp = value
            END IF
         ELSE IF (desc.EQ.10004) THEN ! Pressure
            IF (PO.EQ.rvind) THEN
               PO = value
            END IF
         ELSE IF (desc.EQ.10051) THEN ! Pressure reduced to mean sea level
            IF (PR.EQ.rvind) THEN
               PR = value
            END IF
         ELSE IF (desc.EQ.10061) THEN ! 3-hour pressure change
            IF (PP.EQ.rvind) THEN
               PP = value
            END IF
         ELSE IF (desc.EQ.10063) THEN ! Characteristic of pressure tendency
            IF (AA.EQ.rvind) THEN
               AA = value
            END IF
         ELSE IF (desc.EQ.11011.OR.  ! Wind direction at 10 m
     +           desc.EQ.11001) THEN ! Wind direction
            IF (DD.EQ.rvind) THEN
               DD = value
            END IF
         ELSE IF (desc.EQ.11012.OR.  ! Wind speed at 10 m
     +           desc.EQ.11002) THEN ! Wind speed
            IF (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 (minute_p.NE.rvind .AND. hour.NE.rvind) THEN
               mm = NINT(minute_p)
               hh = NINT(hour)
               IF (mm.EQ.-10) THEN
                  IF (DG_010.EQ.rvind) DG_010 = value
               ELSE IF (mm.GE.-70 .AND. mm.LE.-50) THEN
                  IF (DG_1.EQ.rvind) DG_1 = value
               ELSE IF (mm.EQ.-360 .AND. MOD(hh,6).EQ.0) THEN
                  IF (DG.EQ.rvind) DG = value
               ELSE IF (mm.EQ.-180 .AND. MOD(hh,3).EQ.0
     +                 .AND. MOD(hh,6).NE.0) THEN
                  IF (DG.EQ.rvind) DG = value
               ELSE
C     Actually DG_X is not defined in Kvalobs (but FG_X is!)
                  IF (DG_X.EQ.rvind) DG_X = value
               END IF
            END IF
         ELSE IF (desc.EQ.11041) THEN ! Maximum wind gust speed
            IF (minute_p.NE.rvind .AND. hour.NE.rvind) THEN
               mm = NINT(minute_p)
               hh = NINT(hour)
               IF (mm.EQ.-10) THEN
                  IF (FG_010.EQ.rvind) FG_010 = value
               ELSE IF (mm.GE.-70 .AND. mm.LE.-50) THEN
                  IF (FG_1.EQ.rvind) FG_1 = value
C     For time periods > 1 hour, we choose to decode this as FG for
C     termins 0,6,12,18 if time period is 6 hours, and for termins
C     3,9,15,21 if time period is 3 hours. FG is defined as max wind
C     gust since last synoptic termin (0,6,12,18) and it is unlikely to
C     be reported for other termins than 0,3,6,9...
               ELSE IF (mm.EQ.-360 .AND. MOD(hh,6).EQ.0) THEN
                  IF (FG.EQ.rvind) FG = value
               ELSE IF (mm.EQ.-180 .AND. MOD(hh,3).EQ.0
     +                 .AND. MOD(hh,6).NE.0) THEN
                  IF (FG.EQ.rvind) FG = value
               ELSE
C     All other periods are put in FG_X (arbitrary period)
                  IF (FG_X.EQ.rvind) FG_X = value
               END IF
            END IF
         ELSE IF (desc.EQ.11042) THEN ! Maximum wind speed (10-min mean wind)
            IF (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)
               IF ((mm.EQ.-360 .AND. MOD(hh,6).EQ.0)
     +              .OR. (mm.EQ.-180 .AND. MOD(hh,3).EQ.0
     +              .AND. MOD(hh,6).NE.0) .AND. FX.EQ.rvind) THEN
                  IF (FX.EQ.rvind) FX = value
               ELSE IF (mm.GE.-70 .AND. mm.LE.-50) THEN
                  IF (FX_1.EQ.rvind) FX_1 = value
               ELSE
C     All other periods are put in FX_X (arbitrary period)
                  IF (FX_X.EQ.rvind) FX_X = value
               END IF
            END IF
         ELSE IF (desc.EQ.12104.OR. ! Dry bulb temperature at 2m (data width 16 bits)
     +          desc.EQ.12004.OR.   ! Dry bulb temperature at 2m (12 bits)
     +          desc.EQ.12101.OR.   ! Temperature/dry bulb temperature (16 bits)
     +          desc.EQ.12001) THEN ! Temperature/dry bulb temperature (12 bits)
            IF (TA.EQ.rvind) THEN
               TA = value
            END IF
         ELSE IF (desc.EQ.12106.OR.  ! Dew-point temperature at 2m (16 bits)
     +           desc.EQ.12006.OR.   ! Dew-point temperature at 2m (12 bits)
     +           desc.EQ.12103.OR.   ! Dew-point temperature (16 bits)
     +           desc.EQ.12003) THEN ! Dew-point temperature (12 bits)
            IF (TD.EQ.rvind) THEN
               TD = value
            END IF
         ELSE IF (desc.EQ.12113.OR.  ! Ground minimum temperature at 2m (data width 16 bits)
     +           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.12114.OR.  ! Maximum temperature at 2m, past 12 hours (16 bits)
     +           desc.EQ.12014) THEN ! Maximum temperature at 2m, past 12 hours (12 bits)
            IF (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 (TAX_12.EQ.rvind .AND. idx.GT.2
     +           .AND. ktdexp(idx-1).EQ.4024
     +           .AND. NINT(values(idx-1 + (ksub-1)*kxelem)).EQ.0
     +           .AND. ktdexp(idx-2).EQ.4024) THEN
               IF (NINT(values(idx-2 + (ksub-1)*kxelem)).EQ.-12) THEN
                  IF (TAX_12.EQ.rvind) TAX_12 = value
               ELSE IF (NINT(values(idx-2 + (ksub-1)*kxelem)).EQ.-1)THEN
                  IF (TAX.EQ.rvind) TAX = value
               END IF
            END IF
C     Do we also need to consider 12021 'Maximum temperature at 2m'?
         ELSE IF (desc.EQ.12115.OR.  ! Minimum temperature at 2m, past 12 hours (16 bits)
     +           desc.EQ.12015) THEN ! Minimum temperature at 2m, past 12 hours (12 bits)
            IF (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 (TAN_12.EQ.rvind .AND. idx.GT.2
     +           .AND. ktdexp(idx-1).EQ.4024
     +           .AND. values(idx-1 + (ksub-1)*kxelem).EQ.0
     +           .AND. ktdexp(idx-2).EQ.4024) THEN
               IF (NINT(values(idx-2 + (ksub-1)*kxelem)).EQ.-12) THEN
                  IF (TAN_12.EQ.rvind) TAN_12 = value
               ELSE IF (NINT(values(idx-2 + (ksub-1)*kxelem)).EQ.-1)THEN
                  IF (TAN.EQ.rvind) TAN = value
               END IF
            END IF
C     Do we also need to consider 12022 'Minimum temperature at 2m'?
         ELSE IF (desc.EQ.13003) THEN ! Relative humidity
            IF (UU.EQ.rvind) THEN
               UU = value
            END IF
         ELSE IF (desc.EQ.20001) THEN ! Horizontal visibility
            IF (VV.EQ.rvind) THEN
               VV = value
            END IF
         ELSE IF (desc.EQ.20003) THEN ! Present weather
            IF (WW.EQ.rvind) THEN
               WW = value
            END IF
         ELSE IF (desc.EQ.20004) THEN ! Past weather (1)
            IF (W1.EQ.rvind) THEN
               W1 = value
            END IF
         ELSE IF (desc.EQ.20005) THEN ! Past weather (2)
            IF (W2.EQ.rvind) THEN
               W2 = value
            END IF
         ELSE IF (desc.EQ.20010) THEN ! Cloud cover (total)
            IF (NN.EQ.rvind) THEN
               NN = value
            END IF
         ELSE IF (desc.EQ.31001) THEN
C     Delayed descriptor replication factor; this should be
C     number of cloud layers if previous descriptor is cloud
C     type, according to all WMO recommended templates
C     (but DNMI metar is an exception!)
            IF (ktdexp(idx - 1).EQ.20012) THEN
               IF (NINT(value).LE.4) THEN
                  num_cloud_layers = NINT(value)
               ELSE
                  bad_cloud_data = .TRUE.
               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 (vert_sign_first.EQ.rvind) THEN
                  vert_sign_first = value
               END IF
            ELSE IF (cloud_type_count.LT.3) THEN
               bad_cloud_data = .TRUE. ! There should always be three 020012
                                       ! following first 008002
            ELSE IF (num_cloud_layers.GT.-1) THEN
               cloud_layer = cloud_type_count - 2
               IF (cloud_layer.LE.num_cloud_layers) THEN
                  vert_sign(cloud_layer) = value
               END IF
            ELSE                ! rdb-files always have 0 or 4 cloud layers
               cloud_layer = cloud_type_count - 2
               IF (cloud_layer.LT.5) THEN
                  vert_sign(cloud_layer) = value
               END IF
            END IF
         ELSE IF (desc.EQ.20011 .AND. .NOT.bad_cloud_data) THEN ! Cloud amount
            IF (metar) THEN
               IF (num_cloud_layers.GT.-1) THEN
                  num_cloud_layers = num_cloud_layers + 1
               ELSE
                  num_cloud_layers = 1
               END IF
               NS(num_cloud_layers) = value
            ELSE IF (cloud_type_count.EQ.0) THEN ! First occurrence
               IF (NH.EQ.rvind) THEN
                  NH = value
               END IF
            ELSE IF (cloud_type_count.LT.3) THEN
               bad_cloud_data = .TRUE. ! There should always be three 020012
                                       ! following first 008002
            ELSE IF (num_cloud_layers.GT.-1) THEN
               cloud_layer = cloud_type_count - 2
               IF (cloud_layer.LE.num_cloud_layers) THEN
                  NS(cloud_layer) = value
               END IF
            ELSE                ! rdb-files always have 0 or 4 cloud layers
               cloud_layer = cloud_type_count - 2
               IF (cloud_layer.LT.5) THEN
                  NS(cloud_layer) = value
               END IF
            END IF
         ELSE IF (desc.EQ.20012 .AND. .NOT.bad_cloud_data) THEN ! Cloud type
            IF (metar) THEN
               CC(num_cloud_layers) = value
            ELSE   
               cloud_type_count = cloud_type_count + 1
               IF (cloud_type_count.GT.3) 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
                  ELSE IF (cloud_layer.LT.5) THEN ! rdb-files always have 0 or 4 cloud layers
                     CC(cloud_layer) = value
                  END IF
               ELSE
                  IF (cloud_type_count.EQ.1) THEN
                     IF (CL.EQ.rvind) THEN
                        CL = value
                     END IF
                  ELSE IF (cloud_type_count.EQ.2) THEN
                     IF (CM.EQ.rvind) THEN
                        CM = value
                     END IF
                  ELSE IF (cloud_type_count.EQ.3) THEN
                     IF (CH.EQ.rvind) THEN
                        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 (metar) 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
C     Note that 020013 in individual cloud layers comes
C     AFTER 020012 and therefore must be treated
C     differently than 008002 and 020011
            ELSE IF (cloud_type_count.LT.4) THEN
               bad_cloud_data = .TRUE. ! There should always be three 020012
                                       ! following first 008002
            ELSE IF (num_cloud_layers.GT.-1) THEN
               cloud_layer = cloud_type_count - 3
               IF (cloud_layer.LE.num_cloud_layers) THEN
                  HS(cloud_layer) = value
               END IF
            ELSE ! rdb-files always have 0 or 4 cloud layers
               cloud_layer = cloud_type_count - 3
               IF (cloud_layer.LT.5) THEN
                  HS(cloud_layer) = value
               END IF
            END IF
         ELSE IF (desc.EQ.13023) THEN ! Total precipitation past 24 hours
            IF (RR_24.EQ.rvind) THEN
               RR_24 = value
            END IF
         ELSE IF (desc.EQ.13022) THEN ! Total precipitation past 12 hours
            IF (RR_12.EQ.rvind) THEN
               RR_12 = value
            END IF
         ELSE IF (desc.EQ.13021) THEN ! Total precipitation past 6 hours
            IF (RR_6.EQ.rvind) THEN
               RR_6 = value
            END IF
         ELSE IF (desc.EQ.13020) THEN ! Total precipitation past 3 hours
            IF (RR_3.EQ.rvind) THEN
               RR_3 = value
            END IF
         ELSE IF (desc.EQ.13019) THEN ! Total precipitation past 1 hour
            IF (RR_1.EQ.rvind) THEN
               RR_1 = value
            END IF
         ELSE IF (desc.EQ.13011) THEN ! Total precipitation/total water equivalent
            IF (hour_p.NE.rvind) THEN
               hh = NINT(hour_p)
               IF (hh.EQ.-24) THEN
                  RR_24 = value
               ELSE IF (hh.EQ.-12) THEN
                  RR_12 = value
               ELSE IF (hh.EQ.-6) THEN
                  RR_6 = value
               ELSE IF (hh.EQ.-3) THEN
                  RR_3 = value
               ELSE IF (hh.EQ.-1) THEN
                  RR_1 = value
               END IF
            END IF
         ELSE IF (desc.EQ.13013) THEN ! Total snow depth
C     Don't check for SA.EQ.rvind, because SA might earlier have been set to 0 if
C     EE < 10, which probably means that 20062 has been wrongly encoded, not 13013
            SA = value
         ELSE IF (desc.EQ.20062) THEN ! State of the ground (with or without snow)
            IF (EE.EQ.rvind) THEN
               EE = value
            END IF
            IF (NINT(value).LE.10 .AND. SA.EQ.rvind) THEN
               SA = 0
            END IF
C     Equating 013012 with SS_24 is dubious generally, but this is how
C     SS_24 is encoded in kvalobs for Norwegian avalanche stations (and
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 (hour_p.NE.rvind) THEN
               hh = NINT(hour_p)
               IF (hh.EQ.-1) THEN
                  OT_1 = value
               ELSE IF (hh.EQ.-24) THEN
                  OT_24 = value
               END IF
            END IF
         ELSE IF (desc.EQ.13033) THEN ! Evaporation/evapotranspiration
            IF (hour_p.NE.rvind) THEN
               hh = NINT(hour_p)
               IF (hh.EQ.-1) THEN
                  EV_1 = value
               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 (a3.EQ.rvind) THEN
               a3 = value
            END IF
         ELSE IF (desc.EQ.10009) THEN ! Geopotential height
            IF (hhh.EQ.rvind) THEN
               hhh = value
            END IF
C     Special for ship or marine stations
         ELSE IF (desc.EQ.1011) THEN  ! Ship or mobile land station identifier
            cidx = int(value/1000)
            IF (cidx.GT.0) THEN
               call_sign = cvals(cidx) ! CCITTIA5 data
               call_sign = ctrim(call_sign,9,missing9)
            END IF
         ELSE IF (desc.EQ.1012) THEN ! Direction of motion of moving observing platform
            IF (ds.EQ.rvind) THEN
               ds = value
            END IF
         ELSE IF (desc.EQ.1013) THEN ! Speed of motion of moving observing platform
            IF (vs.EQ.rvind) THEN
               vs = value
            END IF
         ELSE IF (desc.EQ.7062) THEN ! Depth below sea/water surface
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
               surface_data = .TRUE.
            ELSE
               surface_data = .FALSE.
            END IF
         ELSE IF (desc.EQ.22043.OR.  ! Sea/water temperature (15 bits)
     +           desc.EQ.22042.OR.   ! Sea/water temperature (12 bits)
     +           desc.EQ.22049) THEN ! Sea-surface temperature (15 bits)
            IF (TW.EQ.rvind .AND. surface_data)THEN
               TW = value
            END IF
         ELSE IF (desc.EQ.12102.OR.  ! Wet-bulb temperature (16 bits)
     +           desc.EQ.12005) THEN ! Wet-bulb temperature (12 bits)
            IF (TbTbTb.EQ.rvind) THEN
               TbTbTb = value
            END IF
         ELSE IF (desc.EQ.22011) THEN !  Period of waves (instrumentally measured)
            IF (PWA.EQ.rvind) THEN
               PWA = value
            END IF
         ELSE IF (desc.EQ.22012) THEN !  Period of wind waves (visually measured)
            IF (PW.EQ.rvind) THEN
               PW = value
            END IF
         ELSE IF (desc.EQ.22021) THEN ! Heigth of waves
            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 (DW1.EQ.rvind) THEN
               DW1 = value
            ELSE
               DW2 = value
            END IF
         ELSE IF (desc.EQ.22013) THEN ! Period of swell waves
            IF (PW1.EQ.rvind) THEN
               PW1 = value
            ELSE
               PW2 = value
            END IF
         ELSE IF (desc.EQ.22023) THEN ! Height of swell waves
            IF (HW1.EQ.rvind) THEN
               HW1 = value
            ELSE
               HW2 = value
            END IF
         ELSE IF (desc.EQ.22061) THEN ! State of the sea
            IF (SG.EQ.rvind) THEN
               SG = value
            END IF
         ELSE IF (desc.EQ.20033) THEN ! Cause of ice accretion
            IF (XIS.EQ.rvind) THEN
               XIS = value
            END IF
         ELSE IF (desc.EQ.20031) THEN ! Ice deposit (thickness)
            IF (ES.EQ.rvind) THEN
               ES = value
            END IF
         ELSE IF (desc.EQ.20032) THEN ! Rate of ice accretion
            IF (RS.EQ.rvind) THEN
               RS = value
            END IF
         ELSE IF (desc.EQ.20034) THEN ! Sea ice concentration
            IF (CI.EQ.rvind) THEN
               CI = value
            END IF
         ELSE IF (desc.EQ.20037) THEN ! Ice development
            IF (SI.EQ.rvind) THEN
               SI = value
            END IF
         ELSE IF (desc.EQ.20035) THEN ! Amount and type of ice
            IF (BI.EQ.rvind) THEN
               BI = value
            END IF
         ELSE IF (desc.EQ.20038) THEN ! Bearing of ice edge
            IF (DI.EQ.rvind) THEN
               DI = value
            END IF
         ELSE IF (desc.EQ.20036) THEN ! Ice situation
            IF (ZI.EQ.rvind) THEN
               ZI = value
            END IF
         ELSE IF (desc.EQ.1005) THEN ! Buoy/platform identifier
            IF (buoy_id5.EQ.rvind) THEN
               buoy_id5 = value
            END IF
         ELSE IF (desc.EQ.1003) THEN ! WMO region number/geographical area
            IF (wmo_region_number.EQ.rvind) THEN
               wmo_region_number = value
            END IF
         ELSE IF (desc.EQ.1020) THEN ! WMO region sub-area
            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
            cidx = int(value/1000)
            icao_id = cvals(cidx) ! CCITTIA5 data
            icao_id = ctrim(icao_id,8,missing8)
         ELSE IF (desc.EQ.10052) THEN ! Altimeter setting (QNH)
            IF (PH.EQ.rvind) THEN
               PH = value
            END IF
         END IF
      END DO
 
      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 (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
         IF (ind.EQ.-1) ind = 9
C     Remove trailing NULL characters, which some centres erronously
C     insert
         DO WHILE (IACHAR(call_sign(ind:ind)).EQ.0)
            ind = ind - 1
         END DO
         WRITE(*,*)
         WRITE(*,'(A,A)') 'call_sign=',
     +        call_sign(1:ind)
      ELSE IF (buoy_id7.NE.rvind) THEN
C     New templates introduced in 2014 for data category 1 use 001087
C     WMO Marine observing platform extended identifier, 7 digits
         WRITE(*,*)
         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 001005 Buoy/platform identifier. Should we
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(*,*)
            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
      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)
 900     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 (icao_id.NE.missing8) THEN
         WRITE(*,'(A,A)') 'icao_id=',icao_id(1:lenstr(icao_id,1))
      END IF
      IF (name.NE.missing20) THEN
         WRITE(*,'(A,A)') 'name=',name(1:lenstr(name,1))
      END IF
      IF (long_name.NE.missing32) THEN
         WRITE(*,'(A,A)') 'name=',long_name(1:lenstr(long_name,1))
      END IF
      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 (longitude.NE.rvind) THEN
         WRITE(*,'(A,F10.5)') 'lon=',longitude
      END IF
      IF (latitude.NE.rvind) THEN
         WRITE(*,'(A,F10.5)') 'lat=',latitude
      END IF
      IF (height.NE.rvind) THEN ! 7030 has scale 1 but 7003 has scale 0 and is more common,
                                ! 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
         WRITE(*,'(A,F11.1)') 'vs=',vs
      END IF
      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
         WRITE(*,'(A,F10.1)') 'hhh=',hhh
      END IF
      IF (DD.NE.rvind) THEN
         WRITE(*,'(A,I11)') 'DD=',NINT(DD)
      END IF
      IF (FF.NE.rvind) THEN
         WRITE(*,'(A,F11.1)') 'FF=',FF
      END IF
      IF (FX.NE.rvind) THEN
         WRITE(*,'(A,F11.1)') 'FX=',FX
      END IF
      IF (FX_1.NE.rvind) THEN
         WRITE(*,'(A,F9.1)') 'FX_1=',FX_1
      END IF
      IF (FX_X.NE.rvind) THEN
         WRITE(*,'(A,F9.1)') 'FX_X=',FX_X
      END IF
      IF (DG.NE.rvind) THEN
         WRITE(*,'(A,F11.1)') 'DG=',DG
      END IF
      IF (FG.NE.rvind) THEN
         WRITE(*,'(A,F11.1)') 'FG=',FG
      END IF
      IF (DG_010.NE.rvind) THEN
         WRITE(*,'(A,F7.1)') 'DG_010=',DG_010
      END IF
      IF (FG_010.NE.rvind) THEN
         WRITE(*,'(A,F7.1)') 'FG_010=',FG_010
      END IF
      IF (DG_1.NE.rvind) THEN
         WRITE(*,'(A,F9.1)') 'DG_1=',DG_1
      END IF
      IF (FG_1.NE.rvind) THEN
         WRITE(*,'(A,F9.1)') 'FG_1=',FG_1
      END IF
      IF (DG_X.NE.rvind) THEN
         WRITE(*,'(A,F9.1)') 'DG_X=',DG_X
      END IF
      IF (FG_X.NE.rvind) THEN
         WRITE(*,'(A,F9.1)') 'FG_X=',FG_X
      END IF
      IF (TA.NE.rvind) THEN
         WRITE(*,'(A,F11.1)') 'TA=',TA-273.15
      END IF
      IF (TD.NE.rvind) THEN
         WRITE(*,'(A,F11.1)') 'TD=',TD-273.15
      END IF
      IF (TAX_12.NE.rvind) THEN
         WRITE(*,'(A,F7.1)') 'TAX_12=',TAX_12-273.15
      END IF
      IF (TAN_12.NE.rvind) THEN
         WRITE(*,'(A,F7.1)') 'TAN_12=',TAN_12-273.15
      END IF
      IF (TAX.NE.rvind) THEN
         WRITE(*,'(A,F10.1)') 'TAX=',TAX-273.15
      END IF
      IF (TAN.NE.rvind) THEN
         WRITE(*,'(A,F10.1)') 'TAN=',TAN-273.15
      END IF
      IF (TGN_12.NE.rvind) THEN
         WRITE(*,'(A,F7.1)') 'TGN_12=',TGN_12-273.15
      END IF
      IF (TW.NE.rvind) THEN
         WRITE(*,'(A,F11.1)') 'TW=',TW-273.15
      END IF
      IF (TbTbTb.NE.rvind) THEN
         WRITE(*,'(A,F11.1)') 'Tb=',TbTbTb-273.15
      END IF
      IF (UU.NE.rvind) THEN
         WRITE(*,'(A,I11)') 'UU=',NINT(UU)
      END IF
      IF (PO.NE.rvind) THEN
         WRITE(*,'(A,F11.1)') 'PO=',PO/100 ! hPa
      END IF
      IF (PR.NE.rvind) THEN
         WRITE(*,'(A,F11.1)') 'PR=',PR/100 ! hPa
      END IF
      IF (PH.NE.rvind) THEN
         WRITE(*,'(A,F11.1)') 'PH=',PH/100 ! hPa
      END IF
      IF (AA.NE.rvind) THEN
         WRITE(*,'(A,I11)') 'AA=',NINT(AA)
      END IF
      IF (PP.NE.rvind) THEN
         WRITE(*,'(A,F11.1)') 'PP=',PP/100 ! hPa
      END IF
C     Use the BUFR some special values for precipitation, *not* the
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
         WRITE(*,'(A,F8.1)') 'RR_24=',RR_24
      END IF
      IF (RR_12.NE.rvind) THEN
         WRITE(*,'(A,F8.1)') 'RR_12=',RR_12
      END IF
      IF (RR_6.NE.rvind) THEN
         WRITE(*,'(A,F9.1)') 'RR_6=',RR_6
      END IF
      IF (RR_3.NE.rvind) THEN
         WRITE(*,'(A,F9.1)') 'RR_3=',RR_3
      END IF
      IF (RR_1.NE.rvind) THEN
         WRITE(*,'(A,F9.1)') 'RR_1=',RR_1
      END IF
      IF (WW.NE.rvind.AND.WW.LT.200) THEN ! 508-511 and w1w1 (in 333 9 group) ignored here
         IF (NINT(WW).LT.100) THEN
            WRITE(*,'(A,I11)') 'WW=',NINT(WW)
         ELSE
            WRITE(*,'(A,I9)') 'WAWA=',NINT(WW-100)
         END IF
      END IF
      IF (W1.NE.rvind) THEN
         IF (NINT(W1).LT.10) THEN
            WRITE(*,'(A,I11)') 'W1=',NINT(W1)
         ELSE
            WRITE(*,'(A,I10)') 'WA1=',NINT(W1-10)
         END IF
      END IF
      IF (W2.NE.rvind) THEN
         IF (NINT(W2).LT.10) THEN
            WRITE(*,'(A,I11)') 'W2=',NINT(W2)
         ELSE
            WRITE(*,'(A,I10)') 'WA2=',NINT(W2-10)
         END IF
      END IF
      IF (EE.NE.rvind) THEN
         WRITE(*,'(A,I11)') 'EE=',NINT(EE)
      END IF
C     Use the BUFR some special values for SA, *not* the Kvalobs special values
C     So 'trace' is showed as SA=-1, 'Snow cover not continuous' as SA=-2
C     Note that conversion from synop to BUFR normally will set SA=0 if E < 10
      IF (SA.NE.rvind) THEN
         WRITE(*,'(A,I11)') 'SA=',NINT(SA*100)
      END IF
      IF (SS_24.NE.rvind) THEN
         WRITE(*,'(A,I8)') 'SS_24=',NINT(SS_24*100)
      END IF
      IF (VV.NE.rvind) THEN
         WRITE(*,'(A,I11)') 'VV=',NINT(VV)
      END IF
      IF (NN.NE.rvind) THEN
         WRITE(*,'(A,I11)') 'NN=',NNtoWMO_N(NINT(NN))
      END IF
      IF (HL.NE.rvind) THEN
         WRITE(*,'(A,I11)') 'HL=',NINT(HL)
      END IF
      IF (NH.NE.rvind) THEN
         WRITE(*,'(A,I11)') 'NH=',NINT(NH)
      END IF
C     Convert 020012 Cloud type in BUFR into one digit CL (0513), CM
C     (0515) and CH (0509) in TAC
      IF (CL.NE.rvind) THEN
         IF (NINT(CL).GE.30.AND.NINT(CL).LT.40) THEN
            WRITE(*,'(A,I11)') 'CL=',NINT(CL) - 30
         END IF
      END IF
      IF (CM.NE.rvind) THEN
         IF (NINT(CM).GE.20.AND.NINT(CM).LT.30) THEN
            WRITE(*,'(A,I11)') 'CM=',NINT(CM) - 20
         END IF
      END IF
      IF (CH.NE.rvind) THEN
         IF (NINT(CH).GE.10.AND.NINT(CH).LT.20) THEN
            WRITE(*,'(A,I11)') 'CH=',NINT(CH) - 10
         END IF
      END IF
C     Cloud layer data
      DO i=1,num_cloud_layers
         IF (NS(i).NE.rvind) THEN
            WRITE(*,'(A,I1,A,I10)') 'NS',i,'=',NINT(NS(i))
         END IF
         IF (CC(i).NE.rvind) THEN
            WRITE(*,'(A,I1,A,I10)') 'CC',i,'=',NINT(CC(i))
         END IF
         IF (HS(i).NE.rvind) THEN
            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
         WRITE(*,'(A,F11.1)') 'HW=',HW
      END IF
      IF (DW1.NE.rvind) THEN
         WRITE(*,'(A,I10)') 'DW1=',NINT(DW1)
      END IF
      IF (PW1.NE.rvind) THEN
         WRITE(*,'(A,I10)') 'PW1=',NINT(PW1)
      END IF
      IF (HW1.NE.rvind) THEN
         WRITE(*,'(A,F10.1)') 'HW1=',HW1
      END IF
      IF (DW2.NE.rvind) THEN
         WRITE(*,'(A,I10)') 'DW2=',NINT(DW2)
      END IF
      IF (PW2.NE.rvind) THEN
         WRITE(*,'(A,I10)') 'PW2=',NINT(PW2)
      END IF
      IF (HW2.NE.rvind) THEN
         WRITE(*,'(A,F10.1)') 'HW2=',HW2
      END IF
      IF (XIS.NE.rvind) THEN
         WRITE(*,'(A,F10.1)') 'XIS=',XIS
      END IF
      IF (ES.NE.rvind) THEN
         WRITE(*,'(A,F11.1)') 'ES=',ES/10
      END IF
      IF (RS.NE.rvind) THEN
         WRITE(*,'(A,I11)') 'RS=',NINT(RS)
      END IF
      IF (CI.NE.rvind) THEN
         WRITE(*,'(A,I11)') 'CI=',NINT(CI)
      END IF
      IF (SI.NE.rvind) THEN
         WRITE(*,'(A,I11)') 'SI=',NINT(SI)
      END IF
      IF (BI.NE.rvind) THEN
         WRITE(*,'(A,I11)') 'BI=',NINT(BI)
      END IF
      IF (DI.NE.rvind) THEN
         WRITE(*,'(A,I11)') 'DI=',degree2dir(DI)
      END IF
      IF (ZI.NE.rvind) THEN
         WRITE(*,'(A,I11)') 'ZI=',NINT(ZI)
      END IF
      IF (OT_1.NE.rvind) THEN
         WRITE(*,'(A,I9)') 'OT_1=',NINT(OT_1)
      END IF
      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
         WRITE(*,'(A,I9)') 'EV_1=',NINT(EV_1)
      END IF
      IF (EV_24.NE.rvind) THEN
         WRITE(*,'(A,I8)') 'EV_24=',NINT(EV_24)
      END IF
 
      END SUBROUTINE print_surface_values
 
C     -----------------------------------------------------------------
 
      SUBROUTINE print_oceanographic_values(ksub,kxelem,ktdexl,ktdexp,
     +     values,cvals,rectangle,verbose)
C     Identify pressure, temperature etc and print parameter=value to screen
      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)
 
C     Parameters defined in Kvalobs
      REAL*8 HW,PW,TW,
C     Other parameters
     +     year,month,day,hour,minute,
     +     buoy_id,latitude,longitude,sal,
     +     wmo_region_number,wmo_region_subarea
      INTEGER idx,cidx,desc,ind,maxlevel,n,numlevels
      PARAMETER(maxlevel=2000) ! Largest number seen was 499 in a bathy
      CHARACTER*9 call_sign,missing9
      CHARACTER one_bits
      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
C     (--rectangle option set)
      REAL*8 x1,y1,x2,y2
      COMMON /COM_RECTANGLE/  x1,y1,x2,y2
 
C     Functions
      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
      HW = rvind
      PW = rvind
      TW = rvind
      year = rvind
      month = rvind
      day = rvind
      hour = rvind
      minute = rvind
      latitude = rvind
      longitude = rvind
      buoy_id = rvind
      sal = rvind
      wmo_region_number = rvind
      wmo_region_subarea = rvind
      DO n=1,maxlevel
         z(n) = rvind
         w(n) = rvind
         T(n) = rvind
         s(n) = rvind
         d(n) = rvind
         v(n) = rvind
      END DO
 
C     Loop through all expanded descriptors
      n = 0 ! Numbering the subsurface levels
      DO idx=1,ktdexl
         desc = ktdexp(idx)
         value = values(idx + (ksub-1)*kxelem)
 
         IF (ABS(value - rvind)/rvind.LE.0.001) THEN
C     Missing value, nothing to do
            CYCLE
         END IF
 
C     Continue the loop for non missing values
         IF (desc.EQ.4001) THEN ! Year
            IF (year.EQ.rvind) THEN
               year = value
           END IF
         ELSE IF (desc.EQ.4002) THEN ! Month
            IF (month.EQ.rvind) THEN
               month = value
            END IF
         ELSE IF (desc.EQ.4003) THEN ! Day
            IF (day.EQ.rvind) THEN
               day = value
            END IF
         ELSE IF (desc.EQ.4004) THEN ! Hour
            IF (hour.EQ.rvind) THEN
               hour = value
            END IF
         ELSE IF (desc.EQ.4005) THEN ! Minute
            IF (minute.EQ.rvind) THEN
               minute = value
            END IF
         ELSE IF (desc.EQ.5001.OR.  ! Latitude (high accuracy)
     +           desc.EQ.5002) THEN ! Latitude (coarse accuracy)
            IF (latitude.EQ.rvind) THEN
               latitude = value
            END IF
         ELSE IF (desc.EQ.6001.OR.  ! Longitude (high accuracy)
     +           desc.EQ.6002) THEN ! Longitude (coarse accuracy)
            IF (longitude.EQ.rvind) THEN
               longitude = value
            END IF
         ELSE IF (desc.EQ.31001) THEN
C     Delayed descriptor replication factor; this should be number of
C     subsurface levels
            IF (value.EQ.rvind) THEN
               WRITE(*,*) 'WARNING: delayed descriptor replication'
     +              // ' factor 31001 undefined!!!'
               RETURN
            END IF
         ELSE IF (desc.EQ.7062) THEN ! Depth below sea/water surface
            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.7065) THEN ! Water pressure
            n = n + 1           ! new level
            IF (n.GT.maxlevel) THEN
               n = maxlevel
               WRITE(*,*) 'Too many levels! Skipping rest of message'
               GOTO 120
            END IF
            w(n) = value
         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 (n.EQ.0 .AND. TW.EQ.rvind) THEN
               TW = value
            END IF
         ELSE IF (desc.EQ.22064.OR.  ! Salinity [part per thousand] (17 bits)
     +           desc.EQ.22062) THEN ! Salinity [part per thousand] (14 bits)
            IF (n.GT.0 .AND. n.LE.maxlevel) THEN
               s(n) = value
            ELSE IF (n.EQ.0 .AND. sal.EQ.rvind) THEN
               sal = value
            END IF
         ELSE IF (desc.EQ.22004) THEN ! Direction of current
               IF (n.GT.0 .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
               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 .AND. n.LE.maxlevel) THEN
               s(n) = value
            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 (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 (HW.EQ.rvind) THEN
               HW = value
            END IF
         ELSE IF (desc.EQ.1005) THEN ! Buoy/platform identifier
            IF (buoy_id.EQ.rvind) THEN
               buoy_id = value
            END IF
         ELSE IF (desc.EQ.1003) THEN ! WMO region number/geographical area
            IF (wmo_region_number.EQ.rvind) THEN
               wmo_region_number = value
            END IF
         ELSE IF (desc.EQ.1020) THEN ! WMO region sub-area
            IF (wmo_region_subarea.EQ.rvind) THEN
               wmo_region_subarea = value
            END IF
         ELSE IF (desc.EQ.1011) THEN ! Ship or mobile land station identifier
            IF (value.NE.rvind) THEN
               cidx = int(value/1000)
               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
      END DO
 120  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 (buoy_id.NE.rvind.AND.wmo_region_number.NE.rvind
     +     .AND.wmo_region_subarea.NE.rvind) THEN
         WRITE(*,*)
         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(*,*)
         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
         IF (ind.EQ.-1) ind = 9
C     Remove trailing NULL characters, which some centres erronously
C     insert
         DO WHILE (IACHAR(call_sign(ind:ind)).EQ.0)
            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
         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)
 900     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 (longitude.NE.rvind) THEN
         WRITE(*,'(A,F10.5)') 'lon=',longitude
      END IF
      IF (latitude.NE.rvind) THEN
         WRITE(*,'(A,F10.5)') 'lat=',latitude
      END IF
      IF (PW.NE.rvind) THEN
         WRITE(*,'(A,I11)') 'PW=',NINT(PW)
      END IF
      IF (HW.NE.rvind) THEN
         WRITE(*,'(A,F11.1)') 'HW=',HW
      END IF
      IF (TW.NE.rvind) THEN
         WRITE(*,'(A,F11.1)') 'TW=',TW-273.15
      END IF
      IF (sal.NE.rvind) THEN
         WRITE(*,'(A,F11.2)') 'ss=',sal
      END IF
 
      DO n=1,numlevels
C     Do not print level if empty (might happen with last levels in
C     compressed messages, or simply due to bad practice. Setting
C     delayed description replication factor = 1 when ought to be 0 is
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
     +        .AND. w(n).EQ.rvind) THEN
            CYCLE
         END IF
         WRITE(*,'(A,I12)'),'n=',n
         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
            WRITE(*,'(A,F11.2)') 'tt=',t(n)-273.15
         END IF
         IF (s(n).NE.rvind) THEN
            WRITE(*,'(A,F11.2)') 'ss=',s(n)
         END IF
         IF (d(n).NE.rvind) THEN
            WRITE(*,'(A,I11)') 'dc=',NINT(d(n))
         END IF
         IF (v(n).NE.rvind) THEN
            WRITE(*,'(A,F11.1)') 'vc=',v(n)
         END IF
      END DO
 
      END SUBROUTINE print_oceanographic_values
 
C     -----------------------------------------------------------------
 
      SUBROUTINE print_amdar_values(ksub,kxelem,ktdexl,ktdexp,values,
     +     cvals,rectangle,verbose)
C     Identify pressure, temperature etc and print parameter=value to screen
      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)
 
      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,mixing_ratio
      CHARACTER*3 cphase,missing3
      INTEGER idx,cidx,ind
      CHARACTER one_bits
      REAL*8 value
      INTEGER desc
 
C     Variables used for geographical filtering av observations
C     (--rectangle option set)
      REAL*8 x1,y1,x2,y2
      COMMON /COM_RECTANGLE/  x1,y1,x2,y2
 
C     Functions
      INTEGER lenstr
      CHARACTER*8 ctrim ! length must be >= longest variable ctrim is used for
      CHARACTER*3 phase_8004,phase_8009
 
      one_bits = CHAR(255)
      WRITE(missing3,'(3A)') 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
 
C     Initialize all parameters to missing values
      aircraft = missing8
      flight_number = missing8
      cphase = missing3
      osn = rvind
      DD = rvind
      FF = rvind
      TT = rvind
      PP = rvind
      year = rvind
      month = rvind
      day = rvind
      hour = rvind
      minute = rvind
      second = rvind
      latitude = rvind
      longitude = rvind
      flight_level = rvind
      phase = rvind
      mixing_ratio = rvind
 
C     Loop through all expanded descriptors
      DO idx=1,ktdexl
         desc = ktdexp(idx)
         value = values(idx + (ksub-1)*kxelem)
 
         IF (ABS(value - rvind)/rvind.LE.0.001) THEN
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
            cidx = int(value/1000)
            aircraft = cvals(cidx) ! CCITTIA5 data
            aircraft = ctrim(aircraft,8,missing8)
         ELSE IF (desc.EQ.1006) THEN ! Aircraft flight number
            cidx = int(value/1000)
            flight_number = cvals(cidx) ! CCITTIA5 data
            flight_number = ctrim(flight_number,8,missing8)
         ELSE IF (desc.EQ.1023) THEN ! Observation sequence number
            IF (osn.EQ.rvind) osn = value
         ELSE IF (desc.EQ.4001) THEN ! Year
            IF (year.EQ.rvind) THEN
               year = value
            END IF
         ELSE IF (desc.EQ.4002) THEN ! Month
            IF (month.EQ.rvind) THEN
               month = value
            END IF
         ELSE IF (desc.EQ.4003) THEN ! Day
            IF (day.EQ.rvind) THEN
               day = value
            END IF
         ELSE IF (desc.EQ.4004) THEN ! Hour
            IF (hour.EQ.rvind) THEN
               hour = value
            END IF
         ELSE IF (desc.EQ.4005) THEN ! Minute
            IF (minute.EQ.rvind) THEN
               minute = value
            END IF
         ELSE IF (desc.EQ.4006) THEN ! Second
            IF (second.EQ.rvind) THEN
               second = value
            END IF
         ELSE IF (desc.EQ.5001.OR.  ! Latitude (high accuracy)
     +           desc.EQ.5002) THEN ! Latitude (coarse accuracy)
            IF (latitude.EQ.rvind) THEN
               latitude = value
            END IF
         ELSE IF (desc.EQ.6001.OR.  ! Longitude (high accuracy)
     +           desc.EQ.6002) THEN ! Longitude (coarse accuracy)
            IF (longitude.EQ.rvind) THEN
               longitude = value
            END IF
         ELSE IF (desc.EQ.7010.OR.   ! Flight level
     +           desc.EQ.7002.OR.    ! Height or altitude
     +           desc.EQ.10070) THEN ! Indicated aircraft altitude
            IF (flight_level.EQ.rvind) THEN
               flight_level = value
            END IF
         ELSE IF (desc.EQ.8009.OR.  ! Detailed phase of aircraft flight
     +           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 (phase.EQ.rvind) THEN
               phase = value
               cphase = phase_8004(NINT(phase),missing3)
            END IF
         ELSE IF (desc.EQ.12104.OR. ! Dry bulb temperature at 2m (data width 16 bits)
     +          desc.EQ.12004.OR.   ! Dry bulb temperature at 2m (12 bits)
     +          desc.EQ.12101.OR.   ! Temperature/dry bulb temperature (16 bits)
     +          desc.EQ.12001) THEN ! Temperature/dry bulb temperature (12 bits)
            IF (TT.EQ.rvind) THEN
               TT = value
            END IF
         ELSE IF (desc.EQ.11011.OR.  ! Wind direction at 10 m
     +           desc.EQ.11001) THEN ! Wind direction
            IF (DD.EQ.rvind) THEN
               DD = value
            END IF
         ELSE IF (desc.EQ.11012.OR.  ! Wind speed at 10 m
     +           desc.EQ.11002) THEN ! Wind speed
            IF (FF.EQ.rvind) THEN
               FF = value
            END IF
         ELSE IF (desc.EQ.7004) THEN ! Pressure
            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
      END DO
 
      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 (aircraft.NE.missing8) THEN
         WRITE(*,*)
         WRITE(*,'(A,A)') 'aircraft=',aircraft(1:lenstr(aircraft,1))
      ELSE IF (flight_number.NE.missing8) THEN
         WRITE(*,*)
         WRITE(*,'(A,A)') 'aircraft=',
     +        flight_number(1:lenstr(flight_number,1))
      ELSE
         IF (verbose .GT. 1) THEN
            WRITE(*,*)
            WRITE(*,*) 'Aircraft (001008/001006) is missing!!!'
         END IF
         RETURN
      END IF
 
 901  FORMAT('obstime=',I4,'-',I2.2,'-',I2.2,' ',I2.2,':',I2.2,':',I2.2)
 902  FORMAT('obstime=',I4,'-',I2.2,'-',I2.2,' ',I2.2,':',I2.2,':00')
 
      IF (year.NE.rvind.AND.month.NE.rvind.AND.day.NE.rvind
     +     .AND.hour.NE.rvind.AND.minute.NE.rvind) THEN
         IF (second.NE.rvind) THEN
            WRITE(*,901),NINT(year),NINT(month),NINT(day),
     +           NINT(hour),NINT(minute),NINT(second)
         ELSE
            WRITE(*,902),NINT(year),NINT(month),NINT(day),
     +           NINT(hour),NINT(minute)
         END IF
      ELSE IF (verbose .GT. 1) THEN
         WRITE(*,*) 'obstime is missing!!!'
         RETURN
      END IF
      IF (longitude.NE.rvind) THEN
         WRITE(*,'(A,F15.5)') 'lon=',longitude
      END IF
      IF (latitude.NE.rvind) THEN
         WRITE(*,'(A,F15.5)') 'lat=',latitude
      END IF
      IF (osn.NE.rvind) THEN
         WRITE(*,'(A,I12)') 'seqnum=',NINT(osn)
      END IF
      IF (flight_level.NE.rvind) THEN
         WRITE(*,'(A,I6)') 'flight_level=',NINT(flight_level)
      END IF
      IF (phase.NE.rvind) THEN
         WRITE(*,'(A,A3)') 'phase_of_flight=',cphase
      END IF
      IF (DD.NE.rvind) THEN
         WRITE(*,'(A,I16)') 'DD=',NINT(DD)
      END IF
      IF (FF.NE.rvind) THEN
         WRITE(*,'(A,F16.1)') 'FF=',FF
      END IF
      IF (TT.NE.rvind) THEN
         WRITE(*,'(A,F16.1)') 'TT=',TT-273.15
      END IF
      IF (PP.NE.rvind) THEN
         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
 
      END SUBROUTINE print_amdar_values
 
C     -----------------------------------------------------------------
 
C     Remove trailing blanks and nulls from string, returning missing if
C     nothing left
      CHARACTER*(*) FUNCTION ctrim(string,dim,missing)
      IMPLICIT NONE
      CHARACTER*(*) string,missing
      INTEGER dim,ind,ascii
 
      ind = dim
      ascii = IACHAR(string(ind:ind))
      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.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
         ctrim = missing
      END IF
 
      RETURN
      END FUNCTION ctrim
 
 
C     -----------------------------------------------------------------
 
      CHARACTER*3 FUNCTION phase_8004(phase,missing3)
      IMPLICIT NONE
      INTEGER phase
      CHARACTER*3 missing3,indicator(6)
 
      DATA indicator/'res','UNS','LVR','LVW','ASC','DES'/
C                      1     2     3     4     5     6
 
      IF (phase.LE.1 .OR. phase.GT.6) THEN
         phase_8004 = missing3
      ELSE
         phase_8004 = indicator(phase)
      END IF
 
      RETURN
      END FUNCTION phase_8004
 
C     -----------------------------------------------------------------
 
      CHARACTER*3 FUNCTION phase_8009(phase,missing3)
      IMPLICIT NONE
      INTEGER phase
      CHARACTER*3 missing3,indicator(15)
 
C                      0     1     2     3     4     5     6
      DATA indicator/'LVR','LVW','UNS','LVR','LVW','ASC','DES',
     +     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 + 1)
      END IF
 
      RETURN
      END FUNCTION phase_8009
 
C     -----------------------------------------------------------------
 
C     Convert value of NN (in percent) into WMO code (table 2700)
      INTEGER FUNCTION NNtoWMO_N(NN)
      IMPLICIT NONE
      INTEGER NN,N
      IF (NN.EQ.0) THEN
         N = 0
      ELSE IF (NN.LE.15) THEN  ! 1/10 or less
         N = 1
      ELSE IF (NN.LE.30) THEN  ! 2/10 - 3/10
         N = 2
      ELSE IF (NN.LE.45) THEN  ! 4/10
         N = 3
      ELSE IF (NN.LE.55) THEN  ! 5/10
         N = 4
      ELSE IF (NN.LE.65) THEN  ! 6/10
         N = 5
      ELSE IF (NN.LE.80) THEN  ! 7/10 - 8/10
         N = 6
      ELSE IF (NN.LE.99) THEN  ! 9/10 or more, but not 10/10
         N = 7
      ELSE IF (NN.GT.100) THEN ! Sky obscured by fog or other meteorological phenomena
         N = 9                 ! NN will be 113 in WMO BUFR, 109 in surf files
      ELSE
         N = 8
      END IF
 
      NNtoWMO_N = N
      END FUNCTION NNtoWMO_N
 
C     -----------------------------------------------------------------
 
      INTEGER FUNCTION degree2dir(rd)
      IMPLICIT NONE
      REAL*8 rd
      INTEGER id
      IF (NINT(rd).EQ.0) THEN
         id = 0
      ELSE IF (rd.LT.0 .OR. rd.GT.360) THEN
         id = 9
      ELSE IF (rd.LT.22.5 .OR.rd.GE.337.5) THEN
         id = 8                 ! N
      ELSE IF (rd.LT.67.5) THEN
         id = 1                 ! NE
      ELSE IF (rd.LT.112.5) THEN
         id = 2                 ! E
      ELSE IF (rd.LT.157.5) THEN
         id = 3                 ! SE
      ELSE IF (rd.LT.202.5) THEN
         id = 4                 ! S
      ELSE IF (rd.LT.247.5) THEN
         id = 5                 ! SW
      ELSE IF (rd.LT.292.5) THEN
         id = 6                 ! W
      ELSE IF (rd.LT.337.5) THEN
         id = 7                 ! NW
      END IF
 
      degree2dir = id
      END FUNCTION degree2dir
 
C     -----------------------------------------------------------------
 
      SUBROUTINE vss_8001(vss,significance)
C     Returns vertical sounding significance as right justified string:
C     s means surface level
C     S       standard level
C     t       tropopause level
C     M       maximum wind level
C     T       significant level, temperature and/or relative humidity
C     W       significant level, wind
      implicit none
      integer vss,jj
      character*8 significance
 
      significance = '        '
      jj = 8
      IF (btest(vss,7-6)) THEN
         significance(jj:jj) = 'W' ! significant level, wind
         jj = jj-1
      END IF
      IF (btest(vss,7-5)) THEN
         significance(jj:jj) = 'T' ! significant level, temperature and/or relative humidity
         jj = jj-1
      END IF
      IF (btest(vss,7-4)) THEN
         significance(jj:jj) = 'M' ! maximum wind level
         jj = jj-1
      END IF
      IF (btest(vss,7-3)) THEN
         significance(jj:jj) = 't' ! tropopause level
         jj = jj-1
      END IF
      IF (btest(vss,7-2)) THEN
         significance(jj:jj) = 'S' ! standard level
         jj = jj-1
      END IF
      IF (btest(vss,7-1)) THEN
         significance(jj:jj) = 's' ! surface level
         jj = jj-1
      END IF
 
      END SUBROUTINE vss_8001
 
C     -----------------------------------------------------------------
 
      SUBROUTINE vss_8042(vss,significance)
C     Returns vertical sounding significance as right justified string:
C     s means surface level
C     S       standard level
C     t       tropopause level
C     M       maximum wind level
C     T       significant level temperature
C     H       significant level humidity
C     W       significant level wind
C     r       level determined by regional decision
C     These should be the most interesting levels in 008042 (but we might
C     consider including more later)
      implicit none
      integer vss,jj
      character*8 significance
 
      significance = '        '
      jj = 8
      IF (btest(vss,18-15)) THEN
         significance(jj:jj) = 'r' ! level determined by regional decision
         jj = jj-1
      END IF
      IF (btest(vss,18-7)) THEN
         significance(jj:jj) = 'W' ! significant level wind
         jj = jj-1
      END IF
      IF (btest(vss,18-6)) THEN
         significance(jj:jj) = 'H' ! significant level humidity
         jj = jj-1
      END IF
      IF (btest(vss,18-5)) THEN
         significance(jj:jj) = 'T' ! significant level temperature
         jj = jj-1
      END IF
      IF (btest(vss,18-4)) THEN
         significance(jj:jj) = 'M' ! maximum wind level
         jj = jj-1
      END IF
      IF (btest(vss,18-3)) THEN
         significance(jj:jj) = 't' ! tropopause level
         jj = jj-1
      END IF
      IF (btest(vss,18-2)) THEN
         significance(jj:jj) = 'S' ! standard level
         jj = jj-1
      END IF
      IF (btest(vss,18-1)) THEN
         significance(jj:jj) = 's' ! surface level
         jj = jj-1
      END IF
 
      END SUBROUTINE vss_8042
 
C     -----------------------------------------------------------------
 
      integer function lenstr(text,minlen)
c
c  NAME:
c     lenstr
c
c  PURPOSE:
c     Return actual length of text string, excl. trailing blanks.
c
c  SYNOPSIS:
c     integer function lenstr(text,minlen)
c     character*(*) text
c     integer       minlen
c
c  INPUT:
c     text   - text string
c     minlen - the minimum length returned (see NOTES below)
c
c  OUTPUT:
c     lenstr - the string length
c
c  NOTES:
c     Use minlen=1 if you use lenstr to print a text string without
c     trailing blanks. E.g. filenames declared with a large maximum
c     length which hardly ever is longer than one line, but always
c     is printed with trailing blanks on two ore more lines.
c     Example:  write(6,*) filename(1:lenstr(filename,1))
c     (filename(1:0) may abort the program if filename=' ')
c
c-----------------------------------------------------------------------
c  DNMI/FoU  15.03.1995  Anstein Foss
c-----------------------------------------------------------------------
c
      character*(*) text
      integer       k,l,lt,minlen
c
      lt=len(text)
      l=0
      do k=1,lt
        if(text(k:k).ne.' ') l=k
      end do
c
      if(l.lt.minlen) l=minlen
c
      lenstr=l
c
      return
      end
comfilter.f
C (C) Copyright 2010, met.no
C
C This program is free software; you can redistribute it and/or modify
C it under the terms of the GNU General Public License as published by
C the Free Software Foundation; either version 2 of the License, or
C (at your option) any later version.
C
C This program is distributed in the hope that it will be useful, but
C WITHOUT ANY WARRANTY; without even the implied warranty of
C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
C General Public License for more details.
C
C You should have received a copy of the GNU General Public License
C along with this program; if not, write to the Free Software
C Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
C 02110-1301, USA.
 
C     Variables used for filtering av observations (--filter option set
C     in bufrread or bufrdump)
 
      INTEGER dim1_fid, dim2_fid, dim_fiv
      PARAMETER (dim1_fid=10, dim2_fid=5, dim_fiv=1000)
 
      INTEGER fid(1:dim1_fid,1:dim2_fid)
                                ! filter descriptors arrays. fid(i,j) is descriptor
                                ! j in filter criterium i
      CHARACTER*10 fidformat(1:dim1_fid,1:dim2_fid)
                                ! filter descriptor format (I2.2, A9 etc)
                                ! for fid(i,j)
      INTEGER nd_fid(1:dim1_fid)  ! number of descriptors for each filter criterium
                                ! (max value of j for fid(i,j))
      INTEGER nvl_fid(1:dim_fiv) ! nvl_fid(i) is number av lines with filter
                                ! values described by filter descriptor line i
      INTEGER nfidlines         ! number of filter descriptor lines
                                ! (max value of i for fid(i,j))
 
      INTEGER fivI(1:dim_fiv,1:dim2_fid) ! the filter values of integer type 
      CHARACTER*32 fivC(1:dim_fiv,1:dim2_fid) ! the filter values of character type
 
      INTEGER nfivlines         ! number of lines with filter values
 
      COMMON /COM_FILTER/ fid,nd_fid,nvl_fid,nfidlines,nfivlines,fivI
 
      COMMON /COM_FILTERC/ fidformat,fivC
bufr.pm/bufrdump.txt · Last modified: 2018-10-24 09:19:30 by pals