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