Differences

This shows you the differences between two versions of the page.

--- bufr.pm:bufrdump [2018-04-25 08:39:51]
pals
+++ bufr.pm:bufrdump [2022-05-31 09:29:31] (current)
@@ Line 1: / Line 1: @@
-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
+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>
@@ Line 7: / Line 7: @@
 <code fortran bufrdump.F>
-</code>
+C (C) Copyright 2010-2016 MET Norway
-C (C) Copyright 2010-2018 MET Norway
 C
 C This program is free software; you can redistribute it and/or modify
@@ Line 162: / Line 161: @@
       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
@@ Line 209: / Line 209: @@
       kerr = 0
+      metar = .FALSE.
 C     Do a partial expansion without quality control
       kreq(1) = 1               ! All original elements without quality control
@@ Line 227: / Line 228: @@
 C
-C     Using parameter kelem in call to BUFREX might be too big for
+C     Using parameter kelem in call to BUFREX might be too big for some
-C     multisubset messages. Have copied the method used in decode_bufr.F
+C     multisubset messages (or too small for other messages). Have
-C     in libbufr, first calling BUS012 in order to get number of subsets
+C     copied the method used in decode_bufr.F in libbufr, first calling
-C     ksup(6)
+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
@@ Line 239: / Line 240: @@
       END IF
       kxelem = kvals/ksup(6)
-      IF (kxelem .GT. kelem) kxelem = kelem
+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,
@@ Line 260: / Line 263: @@
       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
@@ Line 295: / Line 304: @@
             IF (obstype.EQ.'surface') THEN
                CALL print_surface_values(ksub,kxelem,ktdexl,ktdexp,
-     +              values,cvals,rectangle,verbose)
+     +              values,cvals,rectangle,metar,verbose)
             ELSE IF (obstype(1:8).EQ.'sounding') THEN
                CALL print_sounding_values(ksub,kxelem,ktdexl,ktdexp,
@@ Line 307: / Line 316: @@
             ELSE IF (ksec1(6).LE.1) THEN
                CALL print_surface_values(ksub,kxelem,ktdexl,ktdexp,
-     +              values,cvals,rectangle,verbose)
+     +              values,cvals,rectangle,metar,verbose)
             ELSE IF (ksec1(6).EQ.2) THEN
                CALL print_sounding_values(ksub,kxelem,ktdexl,ktdexp,
@@ Line 335: / Line 344: @@
             IF (obstype.EQ.'surface') THEN
                CALL print_surface_values(ksub,kxelem,ktdexl,ktdexp,
-     +              values,cvals,rectangle,verbose)
+     +              values,cvals,rectangle,metar,verbose)
             ELSE IF (obstype(1:8).EQ.'sounding') THEN
                CALL print_sounding_values(ksub,kxelem,ktdexl,ktdexp,
@@ Line 347: / Line 356: @@
             ELSE IF (ksec1(6).LE.1) THEN
                CALL print_surface_values(ksub,kxelem,ktdexl,ktdexp,
-     +              values,cvals,rectangle,verbose)
+     +              values,cvals,rectangle,metar,verbose)
             ELSE IF (ksec1(6).EQ.2) THEN
                CALL print_sounding_values(ksub,kxelem,ktdexl,ktdexp,
@@ Line 833: / Line 842: @@
       REAL*8 value
       INTEGER idx,cidx,desc,n,maxlevel,numlevels,i,ind,ind2,ind3
-      PARAMETER(maxlevel=10000)
+      PARAMETER(maxlevel=100000)
       REAL*8 P(maxlevel),D(maxlevel),F(maxlevel),
@@ Line 1344: / Line 1353: @@
       SUBROUTINE print_surface_values(ksub,kxelem,ktdexl,ktdexp,values,
-     +     cvals,rectangle,verbose)
+     +     cvals,rectangle,metar,verbose)
 C     Identify pressure, temperature etc and print parameter=value to screen
       IMPLICIT NONE
@@ Line 1355: / Line 1364: @@
       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
@@ Line 1382: / Line 1392: @@
       INTEGER idx,cidx
       INTEGER cloud_type_count  ! Will be increased by one for each 020012
-                                ! (cloud type) encountered (0 initially)
+                                ! (cloud type) encountered (0 initially).
-      INTEGER num_cloud_layers  ! Number of individual cloud layers,
+                                ! Not used for metar
-                                ! set to value of 031001 (delayed
+      INTEGER num_cloud_layers  ! Number of individual cloud layers, set
-                                ! descriptor) if this is met immediately
+                                ! to value of 031001 (delayed descriptor
-                                ! after a 020012 descriptor (-1 initially)
+                                ! 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
@@ Line 1840: / Line 1853: @@
 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
@@ Line 1867: / Line 1881: @@
             END IF
          ELSE IF (desc.EQ.20011 .AND. .NOT.bad_cloud_data) THEN ! Cloud amount
-            IF (cloud_type_count.EQ.0) THEN ! First occurrence
+            IF (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
@@ Line 1886: / Line 1907: @@
             END IF
          ELSE IF (desc.EQ.20012 .AND. .NOT.bad_cloud_data) THEN ! Cloud type
-            cloud_type_count = cloud_type_count + 1
+            IF (metar) THEN
-            IF (cloud_type_count.GT.3) THEN
+               CC(num_cloud_layers) = value
-               cloud_layer = cloud_type_count - 3
+            ELSE
-               IF (num_cloud_layers .GT.-1) THEN
+               cloud_type_count = cloud_type_count + 1
-                  IF (value < 10.0 ! Accept one digit values only
+               IF (cloud_type_count.GT.3) THEN
-     +                 .AND. cloud_layer.LE.num_cloud_layers) 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_layer.LT.5) THEN ! rdb-files always have 0 or 4 cloud layers
+               ELSE
-                  CC(cloud_layer) = value
+                  IF (cloud_type_count.EQ.1) THEN
-               END IF
+                     IF (CL.EQ.rvind) THEN
-            ELSE
+                        CL = value
-               IF (cloud_type_count.EQ.1) THEN
+                     END IF
-                  IF (CL.EQ.rvind) THEN
+                  ELSE IF (cloud_type_count.EQ.2) THEN
-                     CL = value
+                     IF (CM.EQ.rvind) THEN
-                  END IF
+                        CM = value
-               ELSE IF (cloud_type_count.EQ.2) THEN
+                     END IF
-                  IF (CM.EQ.rvind) THEN
+                  ELSE IF (cloud_type_count.EQ.3) THEN
-                     CM = value
+                     IF (CH.EQ.rvind) THEN
-                  END IF
+                        CH = value
-               ELSE IF (cloud_type_count.EQ.3) THEN
+                     END IF
-                  IF (CH.EQ.rvind) THEN
-                     CH = value
                   END IF
                END IF
             END IF
          ELSE IF (desc.EQ.20013 .AND. .NOT.bad_cloud_data) THEN ! Height of base of cloud
-            IF (cloud_type_count.EQ.0) THEN ! First occurrence
+            IF (metar) THEN
+               HS(num_cloud_layers) = value
+            ELSE IF (cloud_type_count.EQ.0) THEN ! First occurrence
                IF (HL.EQ.rvind) THEN
                   HL = value