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

bufrdump: bufrdump.F comfilter.f
$(FC) $(FCFLAGS) -o bufrdump $< -L $(LDIR) -lbufr

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 [--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.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 LOGICAL stop_on_error ! TRUE if program should stop if a libbufr call returns an error 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, + stop_on_error) 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,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,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 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 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 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 IF (stop_on_error) STOP 1 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 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 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 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 (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 IF (stop_on_error) STOP 1 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,stop_on_error) 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=80) argument ! Buffer for next argument CHARACTER(LEN=850) Usage CHARACTER(LEN=900) 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--stop_on_error\tIf a libbufr call returns an ' + // 'error, bufrdump will return\n\t\t\timmediately with ' + // 'exit status 1\n' + // '\t--rectangle x1 y1 x2 y2\n' + // '\t\t\tDecode observations inside rectangle with lower ' + // 'left coordinates\n\t\t\t(x1,y1) (i.e.lon,lat) ' + // 'and upper right coordinates (x2,y2) only\n' + // 'As an alternative to --rectangle you can use one or ' + // 'more of\n' + // '\t--lon1 x1' + // '\tDecode observations with longitude >= x1 only\n' + // '\t--lat1 y1' + // '\tDecode observations with latitude >= y1 only\n' + // '\t--lon2 x2' + // '\tDecode observations with longitude <= x2 only\n' + // '\t--lat2 y2' + // '\tDecode observations with latitude <= y2 only\n' + // '\t\t\tx1,y1,x2,y2 should be degrees or decimal degrees' Help = '\nYou should probably set\n' + // '\texport BUFR_TABLES=/usr/local/lib/bufrtables/\n' + // '\texport PRINT_TABLE_NAMES=false\n' + // '\nUsing --filter will decode only those observations ' + // 'that meet the criteria in\n<filter file>. An example ' + // 'of a filter file is\n\n' + // 'D: 001001 I2.2\n' + // '01\n' + // 'D: 001001 I2.2 001002 I3.3\n' + // '03 895\n' + // '06 252\n' + // 'D: 001011 A9\n' + // 'LMEL \n\n' + // 'which decodes all observations with block number 01, ' + // 'two other specific wmo\nstations and one ship. ' + // 'So far only filtering on exact matches on integer and\n' + // 'character valued BUFR descriptors has been implemented.' + // '\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. 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. '--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 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 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 (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 call_sign,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 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 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) C The missing value returned from libbufr might not be exactly equal to rvind IF (ABS(value - rvind)/rvind.LE.0.001) value = rvind 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) call_sign = 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 (call_sign.NE.missing9) THEN WRITE(*,'(A,A)') 'call_sign=', + call_sign(1:index(call_sign,' ')-1) ELSE IF (verbose .GT. 1) THEN 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 (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 call_sign,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,DG_X,DI,DW1,DW2, + E,EM,ES,EV_1,EV_24,FF,FG,FG_010,FG_1,FG_X,FX,FX_1,FX_X,HL, + 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,TAN,TAX,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 call_sign = 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 DG_X = 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 FG_X = rvind FX = rvind FX_1 = rvind FX_X = 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 TAN = rvind TAX = 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) C The missing value returned from libbufr might not be exactly equal to rvind IF (ABS(value - rvind)/rvind.LE.0.001) value = rvind 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 C DG is treated same way as FG IF (value.NE.rvind + .AND. 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 (value.NE.rvind + .AND. 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 (value.NE.rvind + .AND. minute_p.NE.rvind .AND. hour.NE.rvind) THEN C FX is "Vindhastighet, maks. 10 minutt glidende middel siden C forrige hovedobservasjon, m/s". Assume this is reported for C termins 0,3,6,9... only mm = NINT(minute_p) hh = NINT(hour) 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) 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) 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) 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) 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 (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) THEN C Don't check for SA.EQ.rvind, because SA might earlier have been set to 0 if C E < 10, which probably means that 20062 has been wrongly encoded, not 13013 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 call_sign = 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 (call_sign.NE.missing9) THEN WRITE(*,'(A,A)') 'call_sign=', + call_sign(1:index(call_sign,' ')-1) ELSE IF (buoy_id.NE.rvind.AND.wmo_region_number.NE.rvind + .AND.wmo_region_subarea.NE.rvind) THEN WRITE(*,'(A,I5)') 'buoy_id=',NINT(wmo_region_number)*10000 + + 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_id=',NINT(buoy_id) ELSE IF (verbose .GT. 1) THEN WRITE(*,*) 'Both wmonr, call_sign 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 (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 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=',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,'=',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.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 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