This is an old revision of the document!
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
bufrdump: bufrdump.F comfilter.f $(FC) -o bufrdump $< -L $(LDIR) -lbufr
bufrdump.F:
C (C) Copyright 2010, met.no
C
C This program is free software; you can redistribute it and/or modify
C it under the terms of the GNU General Public License as published by
C the Free Software Foundation; either version 2 of the License, or
C (at your option) any later version.
C
C This program is distributed in the hope that it will be useful, but
C WITHOUT ANY WARRANTY; without even the implied warranty of
C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
C General Public License for more details.
C
C You should have received a copy of the GNU General Public License
C along with this program; if not, write to the Free Software
C Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
C 02110-1301, USA.
C Extract BUFR messages from bufr file(s) using the Fortran program
C bufrdump and print the data as 'parameter=value' lines. See
C usage_verbose for explanation of the options allowed.
C
C Usage: bufrdump <bufr file>
C [--help]
C [--filter <filter file>]
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/emos/bufrtables/
C
C Author: P.Sannes met.no 2010
PROGRAM BUFRDUMP
IMPLICIT NONE
CHARACTER(LEN=80) 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
CHARACTER(LEN=80) filter_file ! File containing the filter criteria
CHARACTER(LEN=80) 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)
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,
+ 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,nlibbufr_errors)
IMPLICIT NONE
INTEGER ibuff(*)
INTEGER ilen ! size in WORDS (not bytes) of the BUFR message read
LOGICAL filter
LOGICAL rectangle
INTEGER nlibbufr_errors ! Number of errors encountered in libbufr calls
INTEGER verbose
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
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
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
C multisubset messages. Have copied the method used in decode_bufr.F
C in libbufr, first calling BUS012 in order to get number of subsets
C ksup(6)
CALL BUS012(ilen,ibuff,ksup,ksec0,ksec1,ksec2,kerr)
IF (kerr.NE.0) THEN
WRITE(*,*) 'ERROR IN BUS012: KERR= ',kerr
nlibbufr_errors = nlibbufr_errors + 1
RETURN
END IF
kxelem = kvals/ksup(6)
IF (kxelem .GT. kelem) 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
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 only. 0 =
C Surface data - land, 1 = Surface data - sea, 2 = Vertical sounding
C (other than satellite)
IF (ksec1(6).GT.2) RETURN
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
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
WRITE(*,*) 'ERROR IN BUSEL2: KERR= ',kerr
nlibbufr_errors = nlibbufr_errors + 1
RETURN
END IF
IF (ksec1(6).LE.1) THEN
CALL print_surface_values(ksub,kxelem,ktdexl,ktdexp,
+ values,cvals,rectangle,verbose)
ELSE IF (ksec1(6).EQ.2) THEN
CALL print_sounding_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
WRITE(*,*) 'ERROR IN BUSEL2: KERR= ',kerr
nlibbufr_errors = nlibbufr_errors + 1
RETURN
END IF
IF (ksec1(6).LE.1) THEN
CALL print_surface_values(ksub,kxelem,ktdexl,ktdexp,
+ values,cvals,rectangle,verbose)
ELSE IF (ksec1(6).EQ.2) THEN
CALL print_sounding_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)
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
CHARACTER(LEN=*) filter_file ! File containing the filter criteria
CHARACTER(LEN=80) argument ! Buffer for next argument
CHARACTER(LEN=850) Usage
CHARACTER(LEN=800) Help
INTEGER iargc,iarg
C Variables used for geographical filtering av observations
REAL*8 x1,y1,x2,y2
COMMON /COM_RECTANGLE/ x1,y1,x2,y2
x1 = -180.00001
y1 = -90.00001
x2 = 180.00001
y2 = 90.00001
Usage = '\nUsage: bufrdump <bufr file> [options]\n\n'
+ // 'Will print section 4 in BUFR messages in <bufr file>'
+ // ' as "parameter=value" lines.'
+ // '\nOptions (may be abbreviated to one character)'
+ // ' are:\n'
+ // '\t--help\t\tPrint this Usage plus some more help\n'
+ // '\t--filter <filter file>\n'
+ // '\t\t\tDecode observations meeting criteria in '
+ // '<filter file> only\n'
+ // '\t--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/emos/bufrtables/\n'
+ // '\texport PRINT_TABLE_NAMES=false\n'
+ // '\nUsing --filter will decode only those observations '
+ // 'that meet the criteria in\n<filter file>. An example '
+ // 'of a filter file is\n\n'
+ // 'D: 001001 I2.2\n'
+ // '01\n'
+ // 'D: 001001 I2.2 001002 I3.3\n'
+ // '03 895\n'
+ // '06 252\n'
+ // 'D: 001011 A9\n'
+ // 'LMEL \n\n'
+ // 'which decodes all observations with block number 01, '
+ // 'two other specific wmo\nstations and one ship. '
+ // 'So far only filtering on excact matches on integer and\n'
+ // 'character valued BUFR descriptors has been implemented.'
+ // '\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.\n'
C Default values
filter = .FALSE.
rectangle = .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. '--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
END SUBROUTINE get_arguments
C -----------------------------------------------------------------
SUBROUTINE err_msg(msg)
IMPLICIT NONE
CHARACTER*80 msg
WRITE(*,*) msg
STOP
END SUBROUTINE err_msg
C -----------------------------------------------------------------
SUBROUTINE err_msg1(msg,kerr)
IMPLICIT NONE
CHARACTER*80 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=80) 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
C Allow comment line
IF (inpline(1:1).EQ.'#') CYCLE filter_loop
i1 = 0
i2 = 0
CALL advance_word(inpline,inplen,i1,i2)
IF (i1 .EQ. 0) EXIT filter_loop
IF (inpline(1:2).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 (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
ifiv = 0
C loop through all subsets:
DO nsub = 1,num_subsets
C loop through all different conditions:
DO i1 = 1,nfidlines
C loop through all filter values 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 values 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
+ GO TO 300 ! 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
300 CONTINUE
END DO ! all filter values 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 one example (high resolution data) where
C 7002 HEIGHT OR ALTITUDE was used as vertical coordinate instead of
C 7004 PRESSURE (10004 was used for pressure). Not able to handle
C that in present code.
SUBROUTINE print_sounding_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 DDDD,missing9
CHARACTER one_bits
REAL*8 value
INTEGER idx,cidx,desc,n,maxlevel,numlevels
PARAMETER(maxlevel=1000)
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
C Variables used for geographical filtering av observations
REAL*8 x1,y1,x2,y2
COMMON /COM_RECTANGLE/ x1,y1,x2,y2
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
DDDD = 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
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) = ' '
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)
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) THEN ! Ship or mobile land station identifier
IF (value.NE.rvind) THEN
cidx = int(value/1000)
DDDD = cvals(cidx) ! CCITTIA5 data
END IF
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
ELSE IF (desc.EQ.5001) THEN ! Latitude (high accuracy)
latitude = value
ELSE IF (desc.EQ.5002) THEN ! Latitude (coarse accuracy)
latitude = value
ELSE IF (desc.EQ.6001) THEN ! Longitude (high accuracy)
longitude = value
ELSE IF (desc.EQ.6002) THEN ! Longitude (coarse accuracy)
longitude = value
ELSE IF (desc.EQ.7001) THEN ! Height of station
height = value
ELSE IF (desc.EQ.7004) THEN ! Pressure
n = n + 1 ! new level
IF (n.GT.maxlevel) THEN
WRITE(*,*) 'Too many levels! Skipping rest of message'
GOTO 110
END IF
P(n) = value
ELSE 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.4086) THEN ! Long time period or displacement [second]
C In WMO templates 004086 comes before 7004 pressure
tp(n+1) = value
ELSE IF (desc.EQ.8001 .AND. value.NE.rvind) THEN ! Vertical sounding significance
C In HIRLAM templates 008001 comes AFTER 7004 pressure
CALL vss_8001(NINT(value),vss(n))
ELSE IF (desc.EQ.8042 .AND. value.NE.rvind) THEN ! Extended vertical sounding significance
C In WMO templates 008042 comes BEFORE 7004 pressure
CALL vss_8042(NINT(value),vss(n+1))
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 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
WRITE(*,*)
IF (II.NE.rvind .AND. iii.NE.rvind) THEN
WRITE(*,'(A,I5.5)') 'wmonr=',NINT(II)*1000 + NINT(iii)
ELSE IF (DDDD.NE.missing9) THEN
WRITE(*,'(A,A)') 'DDDD=',DDDD
ELSE
IF (verbose .GT. 1) THEN
WRITE(*,*) 'Both wmonr and call sign (DDDD)',
+ ' 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 (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
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
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 -----------------------------------------------------------------
SUBROUTINE print_surface_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 icao_id,missing8,spc8
CHARACTER*9 DDDD,missing9,spc9
CHARACTER*20 name,missing20,spc20
C Parameters defined in Kvalobs
REAL*8 AA,BI,CH,CI,CL,CM,DD,DG,DG_010,DG_1,DI,DW1,DW2,
+ E,EM,ES,EV_1,EV_24,FF,FG,FG_010,FG_1,FX,HL,
+ HW,HW1,HW2,NH,NN,OT_1,OT_24,PH,PO,PP,PR,PW,PW1,PW2,
+ RR_1,RR_3,RR_6,RR_12,RR_24,RS,SA,SG,SI,
+ TA,TAN_12,TAX_12,TD,TGN_12,TW,UU,VV,W1,W2,WW,XIS,ZI,
C Other parameters
+ year,month,day,hour,minute,
+ a3,buoy_id,ds,height,hhh,hour_p,II,iii,ix,latitude,longitude,
+ minute_p,TbTbTb,vert_sign_first,vs,wmo_region_number,
+ wmo_region_subarea
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)
INTEGER num_cloud_layers ! Number of individual cloud layers,
! set to value of 031001 (delayed
! descriptor) if this is met immediately
! after a 020012 descriptor (-1 initially)
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
CHARACTER one_bits
REAL*8 value
INTEGER desc,i,mm,hh
INTEGER degree2dir,HStoWMO_HSHS,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
one_bits = CHAR(255)
WRITE(missing8,'(8A)') one_bits,one_bits,one_bits,one_bits,
+ one_bits,one_bits,one_bits,one_bits
missing9 = missing8 // one_bits
missing20 = missing9 // missing9 // one_bits // one_bits
spc8 = ' '
spc9 = ' '
spc20 = ' '
cloud_type_count = 0
bad_cloud_data = .FALSE.
num_cloud_layers = -1
surface_data = .TRUE.
C Initialize all parameters to missing values
hour_p = rvind
minute_p = rvind
icao_id = missing8
DDDD = missing9
name = missing20
AA = rvind
BI = rvind
CH = rvind
CI = rvind
CL = rvind
CM = rvind
DD = rvind
DG = rvind
DG_010 = rvind
DG_1 = rvind
DI = rvind
ds = rvind
DW1 = rvind
DW2 = rvind
ES = rvind
E = rvind
EM = rvind
EV_1 = rvind
EV_24 = rvind
FF = rvind
FG = rvind
FG_010 = rvind
FG_1 = rvind
FX = rvind
HL = rvind
HW = rvind
HW1 = rvind
HW2 = 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
RR_1 = rvind
RR_3 = rvind
RR_6 = rvind
RR_12 = rvind
RR_24 = rvind
RS = rvind
SA = rvind
SG = rvind
SI = rvind
TA = rvind
TAN_12 = rvind
TAX_12 = rvind
TD = rvind
TGN_12 = rvind
TW = rvind
UU = rvind
VV = rvind
vs = 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
vert_sign_first = rvind
II = rvind
iii = rvind
ix = rvind
a3 = rvind
hhh = rvind
ds = rvind
vs = rvind
TbTbTb = rvind
buoy_id = rvind
wmo_region_number = rvind
wmo_region_subarea = 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 (desc.EQ.4024) THEN ! Time period or displacement [hour]
hour_p = value
ELSE IF (desc.EQ.4025) THEN ! Time period or displacement [minute]
minute_p = value
ELSE IF (desc.EQ.1001) THEN ! WMO block number
IF (value.NE.rvind .AND. II.EQ.rvind) THEN
II = value
END IF
ELSE IF (desc.EQ.1002) THEN ! WMO station number
IF (value.NE.rvind .AND. iii.EQ.rvind) THEN
iii = value
END IF
ELSE IF (desc.EQ.1015) THEN ! Station or site name
IF (value.NE.rvind) THEN
cidx = int(value/1000)
IF (cvals(cidx).NE.spc20) THEN
name = cvals(cidx) ! CCITTIA5 data
END IF
END IF
ELSE IF (desc.EQ.2001) THEN ! Type of station
IF (value.NE.rvind .AND. ix.EQ.rvind) THEN
ix = value
END IF
ELSE IF (desc.EQ.4001) THEN ! Year
IF (value.NE.rvind .AND. year.EQ.rvind) THEN
year = value
END IF
ELSE IF (desc.EQ.4002) THEN ! Month
IF (value.NE.rvind .AND. month.EQ.rvind) THEN
month = value
END IF
ELSE IF (desc.EQ.4003) THEN ! Day
IF (value.NE.rvind .AND. day.EQ.rvind) THEN
day = value
END IF
ELSE IF (desc.EQ.4004) THEN ! Hour
IF (value.NE.rvind .AND. hour.EQ.rvind) THEN
hour = value
END IF
ELSE IF (desc.EQ.4005) THEN ! Minute
IF (value.NE.rvind .AND. minute.EQ.rvind) THEN
minute = value
END IF
ELSE IF (desc.EQ.5001) THEN ! Latitude (high accuracy)
IF (value.NE.rvind .AND. latitude.EQ.rvind) THEN
latitude = value
END IF
ELSE IF (desc.EQ.5002) THEN ! Latitude (coarse accuracy)
IF (value.NE.rvind .AND. latitude.EQ.rvind) THEN
latitude = value
END IF
ELSE IF (desc.EQ.6001) THEN ! Longitude (high accuracy)
IF (value.NE.rvind .AND. longitude.EQ.rvind) THEN
longitude = value
END IF
ELSE IF (desc.EQ.6002) THEN ! Longitude (coarse accuracy)
IF (value.NE.rvind .AND. longitude.EQ.rvind) THEN
longitude = value
END IF
ELSE IF (desc.EQ.7001) THEN ! Height of station
IF (value.NE.rvind .AND. height.EQ.rvind) THEN
height = value
END IF
ELSE IF (desc.EQ.7030) THEN ! Height of station ground above mean sea level
IF (value.NE.rvind .AND. height.EQ.rvind) THEN
height = value
END IF
ELSE IF (desc.EQ.10004) THEN ! Pressure
IF (value.NE.rvind .AND. PO.EQ.rvind) THEN
PO = value
END IF
ELSE IF (desc.EQ.10051) THEN ! Pressure reduced to mean sea level
IF (value.NE.rvind .AND. PR.EQ.rvind) THEN
PR = value
END IF
ELSE IF (desc.EQ.10061) THEN ! 3-hour pressure change
IF (value.NE.rvind .AND. PP.EQ.rvind) THEN
PP = value
END IF
ELSE IF (desc.EQ.10063) THEN ! Characteristic of pressure tendency
IF (value.NE.rvind .AND. AA.EQ.rvind) THEN
AA = value
END IF
ELSE IF (desc.EQ.11011) THEN ! Wind direction at 10 m
IF (value.NE.rvind .AND. DD.EQ.rvind) THEN
DD = value
END IF
ELSE IF (desc.EQ.11001) THEN ! Wind direction
IF (value.NE.rvind .AND. DD.EQ.rvind) THEN
DD = value
END IF
ELSE IF (desc.EQ.11012) THEN ! Wind speed at 10 m
IF (value.NE.rvind .AND. FF.EQ.rvind) THEN
FF = value
END IF
ELSE IF (desc.EQ.11002) THEN ! Wind speed
IF (value.NE.rvind .AND. FF.EQ.rvind) THEN
FF = value
END IF
ELSE IF (desc.EQ.11043) THEN ! Maximum wind gust direction
IF (value.NE.rvind .AND. minute_p.NE.rvind) THEN
mm = NINT(minute_p)
IF ((NINT(hour).EQ.0 .OR. NINT(hour).EQ.6
+ .OR. NINT(hour).EQ.12 .OR. NINT(hour).EQ.18)
+ .AND. mm.EQ.-360)
+ THEN
DG = value
ELSE IF (mm.EQ.-10) THEN
DG_010 = value
ELSE IF (mm.EQ.-60) THEN
DG_1 = value
END IF
END IF
ELSE IF (desc.EQ.11041) THEN ! Maximum wind gust speed
IF (value.NE.rvind .AND. minute_p.NE.rvind) THEN
mm = NINT(minute_p)
IF ((NINT(hour).EQ.0 .OR. NINT(hour).EQ.6
+ .OR. NINT(hour).EQ.12 .OR. NINT(hour).EQ.18)
+ .AND. mm.EQ.-360)
+ THEN
FG = value
ELSE IF (mm.EQ.-10) THEN
FG_010 = value
ELSE IF (mm.EQ.-60) THEN
FG_1 = value
END IF
END IF
ELSE IF (desc.EQ.11042) THEN ! Maximum wind speed (10-min mean wind)
IF (value.NE.rvind .AND. FX.EQ.rvind
+ .AND. (NINT(hour).EQ.0 .OR. NINT(hour).EQ.6
+ .OR. NINT(hour).EQ.12 .OR. NINT(hour).EQ.18)
+ .AND. NINT(minute_p).EQ.-360)
+ FX = value
ELSE IF (desc.EQ.12104) THEN ! Dry bulb temperature at 2m (data width 16 bits)
IF (value.NE.rvind .AND. TA.EQ.rvind) THEN
TA = value
END IF
ELSE IF (desc.EQ.12004) THEN ! Dry bulb temperature at 2m (12 bits)
IF (value.NE.rvind .AND. TA.EQ.rvind) THEN
TA = value
END IF
ELSE IF (desc.EQ.12101) THEN ! Temperature/dry bulb temperature (16 bits)
IF (value.NE.rvind .AND. TA.EQ.rvind) THEN
TA = value
END IF
ELSE IF (desc.EQ.12001) THEN ! Temperature/dry bulb temperature (12 bits)
IF (value.NE.rvind .AND. TA.EQ.rvind) THEN
TA = value
END IF
ELSE IF (desc.EQ.12106) THEN ! Dew-point temperature at 2m (16 bits)
IF (value.NE.rvind .AND. TD.EQ.rvind) THEN
TD = value
END IF
ELSE IF (desc.EQ.12006) THEN ! Dew-point temperature at 2m (12 bits)
IF (value.NE.rvind .AND. TD.EQ.rvind) THEN
TD = value
END IF
ELSE IF (desc.EQ.12103) THEN ! Dew-point temperature (16 bits)
IF (value.NE.rvind .AND. TD.EQ.rvind) THEN
TD = value
END IF
ELSE IF (desc.EQ.12003) THEN ! Dew-point temperature (12 bits)
IF (value.NE.rvind .AND. TD.EQ.rvind) THEN
TD = value
END IF
ELSE IF (desc.EQ.12113) THEN ! Ground minimum temperature at 2m (data width 16 bits)
IF (value.NE.rvind .AND. TGN_12.EQ.rvind) THEN
TGN_12 = value
END IF
ELSE IF (desc.EQ.12013) THEN ! Ground minimum temperature at 2m (12 bits)
IF (value.NE.rvind .AND. TGN_12.EQ.rvind) THEN
TGN_12 = value
END IF
ELSE IF (desc.EQ.12114) THEN ! Maximum temperature at 2m, past 12 hours (16 bits)
IF (value.NE.rvind .AND. TAX_12.EQ.rvind) THEN
TAX_12 = value
END IF
ELSE IF (desc.EQ.12014) THEN ! Maximum temperature at 2m, past 12 hours (12 bits)
IF (value.NE.rvind .AND. TAX_12.EQ.rvind) THEN
TAX_12 = value
END IF
ELSE IF (desc.EQ.12111) THEN ! Maximum temperature at height and over period specified
IF (value.NE.rvind .AND. TAX_12.EQ.rvind .AND. idx.GT.2
+ .AND. ktdexp(idx-1).EQ.4024
+ .AND. NINT(values(idx-1 + (ksub-1)*kxelem)).EQ.0
+ .AND. ktdexp(idx-2).EQ.4024
+ .AND. NINT(values(idx-2 + (ksub-1)*kxelem)).EQ.-12)THEN
TAX_12 = value
END IF
C For TAX_12: do we also need to consider 12021 'Maximum temperature at 2m'?
ELSE IF (desc.EQ.12115) THEN ! Minimum temperature at 2m, past 12 hours (16 bits)
IF (value.NE.rvind .AND. TAN_12.EQ.rvind) THEN
TAN_12 = value
END IF
ELSE IF (desc.EQ.12015) THEN ! Minimum temperature at 2m, past 12 hours (12 bits)
IF (value.NE.rvind .AND. TAN_12.EQ.rvind) THEN
TAN_12 = value
END IF
ELSE IF (desc.EQ.12112) THEN ! Minimum temperature at height and over period specified
IF (value.NE.rvind .AND. TAN_12.EQ.rvind .AND. idx.GT.2
+ .AND. ktdexp(idx-1).EQ.4024
+ .AND. values(idx-1 + (ksub-1)*kxelem).EQ.0
+ .AND. ktdexp(idx-2).EQ.4024
+ .AND. values(idx-2 + (ksub-1)*kxelem).EQ.-12) THEN
TAN_12 = value
END IF
C For TAN_12: do we also need to consider 12022 'Minimum temperature at 2m'?
ELSE IF (desc.EQ.13003) THEN ! Relative humidity
IF (value.NE.rvind .AND. UU.EQ.rvind) THEN
UU = value
END IF
ELSE IF (desc.EQ.20001) THEN ! Horizontal visibility
IF (value.NE.rvind .AND. VV.EQ.rvind) THEN
VV = value
END IF
ELSE IF (desc.EQ.20003) THEN ! Present weather
IF (value.NE.rvind .AND. WW.EQ.rvind) THEN
WW = value
END IF
ELSE IF (desc.EQ.20004) THEN ! Past weather (1)
IF (value.NE.rvind .AND. W1.EQ.rvind) THEN
W1 = value
END IF
ELSE IF (desc.EQ.20005) THEN ! Past weather (2)
IF (value.NE.rvind .AND. W2.EQ.rvind) THEN
W2 = value
END IF
ELSE IF (desc.EQ.20010) THEN ! Cloud cover (total)
IF (value.NE.rvind .AND. 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
IF (ktdexp(idx - 1).EQ.20012) THEN
IF (value.EQ.rvind) THEN
WRITE(*,*) 'WARNING: delayed descriptor replication'
+ // ' factor after 020012 undefined !!!'
bad_cloud_data = .TRUE.
ELSE
num_cloud_layers = NINT(value)
END IF
END IF
ELSE IF (desc.EQ.8002 .AND. .NOT.bad_cloud_data) THEN ! Vertical significance (surface observations)
IF (cloud_type_count.EQ.0) THEN ! First occurrence
IF (value.NE.rvind .AND. vert_sign_first.EQ.rvind) THEN
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
IF (value.NE.rvind) THEN
vert_sign(cloud_layer) = value
END IF
END IF
ELSE ! rdb-files always have 0 or 4 cloud layers
cloud_layer = cloud_type_count - 2
IF (cloud_layer.LT.5) THEN
IF (value.NE.rvind) THEN
vert_sign(cloud_layer) = value
END IF
END IF
END IF
ELSE IF (desc.EQ.20011 .AND. .NOT.bad_cloud_data) THEN ! Cloud amount
IF (cloud_type_count.EQ.0) THEN ! First occurrence
IF (value.NE.rvind .AND. 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
IF (value.NE.rvind) THEN
NS(cloud_layer) = value
END IF
END IF
ELSE ! rdb-files always have 0 or 4 cloud layers
cloud_layer = cloud_type_count - 2
IF (cloud_layer.LT.5) THEN
IF (value.NE.rvind) THEN
NS(cloud_layer) = value
END IF
END IF
END IF
ELSE IF (desc.EQ.20012 .AND. .NOT.bad_cloud_data) THEN ! Cloud type
cloud_type_count = cloud_type_count + 1
IF (cloud_type_count.GT.3) THEN
cloud_layer = cloud_type_count - 3
IF (num_cloud_layers .GT.-1) THEN
IF (cloud_layer.LE.num_cloud_layers) THEN
IF (value.NE.rvind) THEN
CC(cloud_layer) = value
END IF
END IF
ELSE IF (cloud_layer.LT.5) THEN ! rdb-files always have 0 or 4 cloud layers
IF (value.NE.rvind) THEN
CC(cloud_layer) = value
END IF
END IF
ELSE
IF (cloud_type_count.EQ.1) THEN
IF (value.NE.rvind .AND. CL.EQ.rvind) THEN
CL = value
END IF
ELSE IF (cloud_type_count.EQ.2) THEN
IF (value.NE.rvind .AND. CM.EQ.rvind) THEN
CM = value
END IF
ELSE IF (cloud_type_count.EQ.3) THEN
IF (value.NE.rvind .AND. CH.EQ.rvind) THEN
CH = value
END IF
END IF
END IF
ELSE IF (desc.EQ.20013 .AND. .NOT.bad_cloud_data) THEN ! Height of base of cloud
IF (cloud_type_count.EQ.0) THEN ! First occurrence
IF (value.NE.rvind .AND. 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
IF (value.NE.rvind) THEN
HS(cloud_layer) = value
END IF
END IF
ELSE ! rdb-files always have 0 or 4 cloud layers
IF (value.NE.rvind) THEN
cloud_layer = cloud_type_count - 3
IF (cloud_layer.LT.5) THEN
IF (value.NE.rvind) THEN
HS(cloud_layer) = value
END IF
END IF
END IF
END IF
ELSE IF (desc.EQ.13023) THEN ! Total precipitation past 24 hours
IF (value.NE.rvind .AND. RR_24.EQ.rvind) THEN
RR_24 = value
END IF
ELSE IF (desc.EQ.13022) THEN ! Total precipitation past 12 hours
IF (value.NE.rvind .AND. RR_12.EQ.rvind) THEN
RR_12 = value
END IF
ELSE IF (desc.EQ.13021) THEN ! Total precipitation past 6 hours
IF (value.NE.rvind .AND. RR_6.EQ.rvind) THEN
RR_6 = value
END IF
ELSE IF (desc.EQ.13020) THEN ! Total precipitation past 3 hours
IF (value.NE.rvind .AND. RR_3.EQ.rvind) THEN
RR_3 = value
END IF
ELSE IF (desc.EQ.13019) THEN ! Total precipitation past 1 hour
IF (value.NE.rvind .AND. RR_1.EQ.rvind) THEN
RR_1 = value
END IF
ELSE IF (desc.EQ.13011) THEN ! Total precipitation/total water equivalent
IF (value.NE.rvind .AND. hour_p.NE.rvind) THEN
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
IF (value.NE.rvind .AND. SA.EQ.rvind) THEN
SA = value
END IF
ELSE IF (desc.EQ.20062) THEN ! State of the ground (with or without snow)
IF (value.NE.rvind) THEN
IF (NINT(value).LE.10) THEN ! E in TAC (3Ejjj)
E = value
IF (SA.EQ.rvind) THEN
SA = 0
END IF
ELSE IF (NINT(value).LE.20) THEN ! E' in TAC (4Esss)
EM = value
END IF
END IF
ELSE IF (desc.EQ.14031) THEN ! Total sunshine
IF (value.NE.rvind .AND. 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 (value.NE.rvind .AND. 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
C Special for high altitude stations
ELSE IF (desc.EQ.7004) THEN ! Pressure (location class)
IF (value.NE.rvind .AND. a3.EQ.rvind) THEN
a3 = value
END IF
ELSE IF (desc.EQ.10009) THEN ! Geopotential height
IF (value.NE.rvind .AND. hhh.EQ.rvind) THEN
hhh = value
END IF
ELSE IF (desc.EQ.10008) THEN ! Geopotential (20 bits)
IF (value.NE.rvind .AND. hhh.EQ.rvind) THEN
hhh = value * 9.8
END IF
ELSE IF (desc.EQ.10003) THEN ! Geopotential (17 bits)
IF (value.NE.rvind .AND. hhh.EQ.rvind) THEN
hhh = value * 9.8
END IF
C Special for ship or marine stations
ELSE IF (desc.EQ.1011) THEN ! Ship or mobile land station identifier
IF (value.NE.rvind) THEN
cidx = int(value/1000)
IF (cvals(cidx).NE.spc9) THEN
DDDD = cvals(cidx) ! CCITTIA5 data
END IF
END IF
ELSE IF (desc.EQ.1012) THEN ! Direction of motion of moving observing platform
IF (value.NE.rvind .AND. ds.EQ.rvind) THEN
ds = value
END IF
ELSE IF (desc.EQ.1013) THEN ! Speed of motion of moving observing platform
IF (value.NE.rvind .AND. vs.EQ.rvind) THEN
vs = value
END IF
ELSE IF (desc.EQ.7062) THEN ! Depth below sea/water surface
IF (value.NE.rvind) THEN
C Some buoy reports starts with depth 1.5 m, others starts with 0 m
C then 1.5 m and always have same sea/water temperature for these 2
C levels, so it seems like 0 m should be considered equivalent with 1.5 m
IF (value.LT.1.6) THEN
surface_data = .TRUE.
ELSE
surface_data = .FALSE.
END IF
END IF
ELSE IF (desc.EQ.22043) THEN ! Sea/water temperature (15 bits)
IF (value.NE.rvind .AND. TW.EQ.rvind .AND. surface_data)THEN
TW = value
END IF
ELSE IF (desc.EQ.22042) THEN ! Sea/water temperature (12 bits)
IF (value.NE.rvind .AND. TW.EQ.rvind
+ .AND. surface_data) THEN
TW = value
END IF
ELSE IF (desc.EQ.12102) THEN ! Wet-bulb temperature (16 bits)
IF (value.NE.rvind .AND. TbTbTb.EQ.rvind) THEN
TbTbTb = value
END IF
ELSE IF (desc.EQ.12005) THEN ! Wet-bulb temperature (12 bits)
IF (value.NE.rvind .AND. TbTbTb.EQ.rvind) THEN
TbTbTb = value
END IF
ELSE IF (desc.EQ.22011) THEN ! Period of waves (BUFR doesn't distinguish between PW and PWA)
IF (value.NE.rvind .AND. PW.EQ.rvind) THEN
PW = value
END IF
ELSE IF (desc.EQ.22021) THEN ! Heigth of waves (BUFR doesn't distinguish between HW and HWA)
IF (value.NE.rvind .AND. HW.EQ.rvind) THEN
HW = value
END IF
ELSE IF (desc.EQ.22003) THEN ! Direction of swell waves
IF (value.NE.rvind) THEN
IF (DW1.EQ.rvind) THEN
DW1 = value
ELSE
DW2 = value
END IF
END IF
ELSE IF (desc.EQ.22013) THEN ! Period of swell waves
IF (value.NE.rvind) THEN
IF (PW1.EQ.rvind) THEN
PW1 = value
ELSE
PW2 = value
END IF
END IF
ELSE IF (desc.EQ.22023) THEN ! Height of swell waves
IF (value.NE.rvind) THEN
IF (HW1.EQ.rvind) THEN
HW1 = value
ELSE
HW2 = value
END IF
END IF
ELSE IF (desc.EQ.22061) THEN ! State of the sea
IF (value.NE.rvind .AND. SG.EQ.rvind) THEN
SG = value
END IF
ELSE IF (desc.EQ.20033) THEN ! Cause of ice accretion
IF (value.NE.rvind .AND. XIS.EQ.rvind) THEN
XIS = value
END IF
ELSE IF (desc.EQ.20031) THEN ! Ice deposit (thickness)
IF (value.NE.rvind .AND. ES.EQ.rvind) THEN
ES = value
END IF
ELSE IF (desc.EQ.20032) THEN ! Rate of ice accretion
IF (value.NE.rvind .AND. RS.EQ.rvind) THEN
RS = value
END IF
ELSE IF (desc.EQ.20034) THEN ! Sea ice concentration
IF (value.NE.rvind .AND. CI.EQ.rvind) THEN
CI = value
END IF
ELSE IF (desc.EQ.20037) THEN ! Ice development
IF (value.NE.rvind .AND. SI.EQ.rvind) THEN
SI = value
END IF
ELSE IF (desc.EQ.20035) THEN ! Amount and type of ice
IF (value.NE.rvind .AND. BI.EQ.rvind) THEN
BI = value
END IF
ELSE IF (desc.EQ.20038) THEN ! Bearing of ice edge
IF (value.NE.rvind .AND. DI.EQ.rvind) THEN
DI = value
END IF
ELSE IF (desc.EQ.20036) THEN ! Ice situation
IF (value.NE.rvind .AND. ZI.EQ.rvind) THEN
ZI = value
END IF
ELSE IF (desc.EQ.1005) THEN ! Buoy/platform identifier
IF (value.NE.rvind .AND. buoy_id.EQ.rvind) THEN
buoy_id = value
END IF
ELSE IF (desc.EQ.1003) THEN ! WMO region number/geographical area
IF (value.NE.rvind .AND. wmo_region_number.EQ.rvind) THEN
wmo_region_number = value
END IF
ELSE IF (desc.EQ.1020) THEN ! WMO region sub-area
IF (value.NE.rvind .AND. wmo_region_subarea.EQ.rvind) THEN
wmo_region_subarea = value
END IF
C Special for metar
ELSE IF (desc.EQ.1063) THEN ! ICAO location indicator
IF (value.NE.rvind) THEN
cidx = int(value/1000)
IF (cvals(cidx).NE.spc8) THEN
icao_id = cvals(cidx) ! CCITTIA5 data
END IF
END IF
ELSE IF (desc.EQ.10052) THEN ! Altimeter setting (QNH)
IF (value.NE.rvind .AND. 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
WRITE(*,*)
IF (II.NE.rvind .AND. iii.NE.rvind) THEN
WRITE(*,'(A,I5.5)') 'wmonr=',NINT(II)*1000 + NINT(iii)
ELSE IF (DDDD.NE.missing9) THEN
WRITE(*,'(A,A)') 'DDDD=',DDDD
ELSE IF (buoy_id.NE.rvind.AND.wmo_region_number.NE.rvind
+ .AND.wmo_region_subarea.NE.rvind) THEN
WRITE(*,'(A,I5)') 'buoy=',NINT(wmo_region_number)*10000
+ + NINT(wmo_region_subarea)*1000 + NINT(buoy_id)
ELSE IF (buoy_id.NE.rvind.AND.buoy_id.GT.1000) THEN
C Old drau files (wrongly) includes wmo_region_number and
C wmo_region_subarea in buoy_id
WRITE(*,'(A,I5)') 'buoy=',NINT(buoy_id)
ELSE
IF (verbose .GT. 1) THEN
WRITE(*,*) 'Both wmonr, call sign (DDDD) and buoy_id',
+ ' 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 (icao_id.NE.missing8) THEN
WRITE(*,'(A,A)') 'icao_id=',icao_id
END IF
IF (name.NE.missing20) THEN
WRITE(*,'(A,A)') 'name=',name
END IF
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
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 (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 (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 (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 (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 (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 Precipitation = -0.1 means "trace" (less than 0.05 kg/m2), is
C represented as 0 in Kvalobs
IF (RR_24.NE.rvind) THEN
IF (NINT(RR_24*10).EQ.-1) THEN
WRITE(*,'(A,F8.1)') 'RR_24=',0.0
ELSE
WRITE(*,'(A,F8.1)') 'RR_24=',RR_24
END IF
END IF
IF (RR_12.NE.rvind) THEN
IF (NINT(RR_12*10).EQ.-1) THEN
WRITE(*,'(A,F8.1)') 'RR_12=',0.0
ELSE
WRITE(*,'(A,F8.1)') 'RR_12=',RR_12
END IF
END IF
IF (RR_6.NE.rvind) THEN
IF (NINT(RR_6*10).EQ.-1) THEN
WRITE(*,'(A,F9.1)') 'RR_6=',0.0
ELSE
WRITE(*,'(A,F9.1)') 'RR_6=',RR_6
END IF
END IF
IF (RR_3.NE.rvind) THEN
IF (NINT(RR_3*10).EQ.-1) THEN
WRITE(*,'(A,F9.1)') 'RR_3=',0.0
ELSE
WRITE(*,'(A,F9.1)') 'RR_3=',RR_3
END IF
END IF
IF (RR_1.NE.rvind) THEN
IF (NINT(RR_1*10).EQ.-1) THEN
WRITE(*,'(A,F9.1)') 'RR_1=',0.0
ELSE
WRITE(*,'(A,F9.1)') 'RR_1=',RR_1
END IF
END IF
IF (WW.NE.rvind.AND.WW.LT.200) THEN ! 508-511 and w1w1 (in 333 9 group) ignored here
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 (E.NE.rvind) THEN
WRITE(*,'(A,I12)') 'E=',NINT(E)
END IF
IF (EM.NE.rvind) THEN
WRITE(*,'(A,I11)') 'EM=',NINT(EM)-10
END IF
C SA has some special values in BUFR as well as in Kvalobs
C Note that conversion from synop to BUFR normally will set SA=0 if E < 10
c$$$ IF (E.NE.rvind .AND. NINT(E).LE.10
c$$$ + .AND. (SA.EQ.rvind .OR. NINT(SA).EQ.0) ) THEN
c$$$ WRITE(*,'(A,I11)') 'SA=',-1 ! No snow
IF (SA.NE.rvind) THEN
IF (NINT(SA*100).EQ.-1) THEN ! Trace: less than 0.5 cm snow
WRITE(*,'(A,I11)') 'SA=',0
ELSE IF (NINT(SA*100).EQ.-2) THEN ! Snow cover not continuos
WRITE(*,'(A,I11)') 'SA=',-1
ELSE IF (NINT(SA*100).EQ.0) THEN ! 0 snow coded as -1 in Kvalobs. Stupid but true
WRITE(*,'(A,I11)') 'SA=',-1
ELSE
WRITE(*,'(A,I11)') 'SA=',NINT(SA*100)
END IF
END IF
IF (VV.NE.rvind) THEN
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=',NNtoWMO_N(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(CL).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,'=',NNtoWMO_N(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,'=',HStoWMO_HSHS(HS(i))
END IF
END DO
IF (SG.NE.rvind) THEN
WRITE(*,'(A,I11)') 'SG=',NINT(SG)
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 (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 (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 -----------------------------------------------------------------
C Convert value of HS (in meter) into WMO code for hshs (table 1677)
INTEGER FUNCTION HStoWMO_HSHS(HS)
IMPLICIT NONE
REAL*8 HS
INTEGER HSHS
IF (NINT(HS).LT.1800) THEN
HSHS=NINT(HS/30)-1
ELSE IF (NINT(HS).LT.10500) THEN
HSHS=NINT((HS-1800)/300)+55
ELSE IF (NINT(HS).LE.21000) THEN
HSHS=NINT((HS-10500)/1500)+80
ELSE
HSHS=89
END IF
HStoWMO_HSHS = HSHS
RETURN
END FUNCTION HStoWMO_HSHS
C -----------------------------------------------------------------
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.10) 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 minlen
c
lt=len(text)
l=0
do k=1,lt
if(text(k:k).ne.' ') l=k
end do
c
if(l.lt.minlen) l=minlen
c
lenstr=l
c
return
end
comfilter.f:
C (C) Copyright 2010, met.no
C
C This program is free software; you can redistribute it and/or modify
C it under the terms of the GNU General Public License as published by
C the Free Software Foundation; either version 2 of the License, or
C (at your option) any later version.
C
C This program is distributed in the hope that it will be useful, but
C WITHOUT ANY WARRANTY; without even the implied warranty of
C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
C General Public License for more details.
C
C You should have received a copy of the GNU General Public License
C along with this program; if not, write to the Free Software
C Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
C 02110-1301, USA.
C Variables used for filtering av observations (--filter option set
C in readbufr 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