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
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
C [--help]
C [--filter ]
C [--obstype]
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 [options]\n\n'
+ // 'Will print section 4 in BUFR messages in '
+ // ' 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 \n'
+ // '\t\t\tDecode observations meeting criteria in '
+ // ' only\n'
+ // '\t--obstype \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. 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
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