tocgrib.F90

Go to the documentation of this file.

 
 
PROGRAM tocgrib
  parameter(mxsiz3=1000000)
 
  INTEGER dum
  INTEGER jgds(200)
  INTEGER mpds(200)
  INTEGER kgds(200)
  INTEGER kpds(200)
  INTEGER mapnum
  INTEGER nbits
  INTEGER nbul
  INTEGER nparm
  INTEGER punum
  INTEGER,dimension(28) :: hexpds
 
  CHARACTER * 6 bulhed
  CHARACTER * 100 cparm
  CHARACTER * 20 desc
  CHARACTER * 3 eoml
  CHARACTER * 1 grib(mxsiz3)
  CHARACTER * 200 fileb,filei,fileo
  CHARACTER * 6 envvar
  CHARACTER * 4 kwbx
  CHARACTER * 1 pds(28)
  CHARACTER * 1 pdsl(28)
  CHARACTER * 1 csep(80)
  CHARACTER * 132 title
  integer,parameter :: lenhead=21
  CHARACTER * 1 wmohdr(lenhead)
 
  LOGICAL iw3pds
 
  hexpds = 0
  lugb = 11
  lugi = 31
  lugo = 51
 
  ! Get parm field with up to 100 characters. Parm field should
  ! contain the originating center part of the WMO Header.
  cparm = '    '
  kwbx  = 'KWBC'
  CALL w3as00(nparm, cparm, ier)
  IF (ier .EQ. 0) THEN
     IF (nparm .EQ. 0 .OR. cparm(1:4) .EQ. '    ') THEN
        print *,'THERE IS A PARM FIELD BUT IT IS EMPTY'
        print *,'OR BLANK, I WILL USE THE DEFAULT KWBC'
     ELSE
        kwbx(1:4) = cparm(1:4)
     END IF
  ELSE IF (ier .EQ. 2 .OR. ier .EQ. 3) THEN
     print *,'W3AS00 ERROR = ', ier
     print *,'THERE IS NO PARM FIELD, I USED DEFAULT KWBC'
  ELSE
     print *,'W3AS00 ERROR = ', ier
  END IF
  print *,'NPARM = ',nparm
  print *,'CPARM = ',cparm(1:4)
  print *,'KWBX  = ',kwbx(1:4)
 
  ! Read GRIB data and index file names from the FORTnn
  ! environment variables, and open the files.
  envvar='FORT  '
  write(envvar(5:6),fmt='(I2)') lugb
  call getenv(envvar,fileb)
  write(envvar(5:6),fmt='(I2)') lugi
  call getenv(envvar,filei)
 
  call baopenr(lugb,fileb,iret1)
  if ( iret1  .ne. 0 ) then
     write(6,fmt='(" Error opening GRIB file: ", A200)') fileb
     write(6,fmt='(" baopenr error = ",I5)') iret1
     stop 10
  endif
 
  ! Open GRIB index file. If doesn't open, use just the data
  ! file.
  call baopenr(lugi,filei,iret2)
  if ( iret2  .ne. 0 ) then
     lugi=0
  endif
 
  ! Read output GRIB bulletin file name from FORTnn
  ! environment variable, and open file.
  write(envvar(5:6),fmt='(I2)') lugo
  call getenv(envvar,fileo)
  call baopenw(lugo,fileo,iret1)
  if ( iret1  .ne. 0 ) then
     write(6,fmt='(" Error opening GRIB file: ",A200)') fileo
     write(6,fmt='(" baopenw error = ",I5)') iret1
     stop 20
  endif
 
  iret = 0
  iopt=2
  insize=19
  nbul = 0
 
  ! loop through input control records.
  nrec = 0
  foreachelement: do
 
     READ (*,66,iostat=ios) (hexpds(j),j=1,12), &
          (hexpds(j),j=17,20), punum, desc
66   FORMAT(3(2x,4z2.2),3x,4z2.2,6x,i3,1x,a20)
     if ( ios .ne. 0 ) then
        write(6,fmt='(" Error reading PDS from input file. iostat = ", i5)') ios
        exit
     endif
     pds=char(hexpds)
 
     ! exit loop, if no more bulletins in input cards
     IF ( mova2i(pds(1)) .EQ. 255) exit
 
     nrec = nrec + 1
     WRITE(6,fmt='(/,''***********************************'', ''********************************************'')')
     print *,'Start new record no. = ',nrec
     WRITE (6,fmt='('' INPUT PDS, PUNUM'', '' & DESC...DESIRED GRIB MAPS LISTED ON FOLLOWING '', ' // &
          ' ''LINES...'',/,4X,3(2X,4z2.2),3X,4z2.2,6X,I3,1X, A20)') (hexpds(j),j=1,12), &
          (hexpds(j),j=17,20), punum, desc
 
     ! Read WNO Header associated with this element
     READ (*,iostat=ios,fmt='(4X,I3,2X,I2,2X,A6,1X,I3,24X,A3)') &
          mapnum,nbits, bulhed, dum, eoml
     WRITE (6,fmt='(4X,I3,2X,I2,2X,A6,1X,I3,24X,A3)') &
          mapnum,nbits, bulhed, dum, eoml
     if ( ios .ne. 0 ) then
        write(6,fmt='(" Error reading PDS from input file. iostat =", i6)') ios
     endif
 
     ! Set up 25 word PDS array of GRIB field to read
     jrew = 0
     jgds = -1
     mpds = -1
     mpds(3) = mova2i(pds(7))
     mpds(5) = mova2i(pds(9))
     mpds(6) = mova2i(pds(10))
     mpds(7) = mova2i(pds(11)) * 256 + mova2i(pds(12))
     mpds(14) = mova2i(pds(19))
     mpds(15) = mova2i(pds(20))
 
     ! Read and return packed GRIB field
     CALL getgbp(lugb,lugi,mxsiz3,jrew,mpds,jgds, &
          itot,krew,kpds,kgds,grib,iret)
     IF (iret.NE.0) THEN
        IF (iret.LT.96) print *,'GETGB-W3FI63: ERROR = ',iret
        IF (iret.EQ.96) print *,'GETGB: ERROR READING INDEX FILE'
        IF (iret.EQ.97) print *,'GETGB: ERROR READING GRIB FILE'
        IF (iret.EQ.98) THEN
           print *,'GETGB ERROR: NUM. OF DATA POINTS GREATER THAN JF'
        END IF
        IF (iret.EQ.99) print *,'GETGB ERROR: REQUEST NOT FOUND'
        IF (iret.GT.99) print *,'GETGB ERROR = ',iret
        cycle
     END IF
     print *,'RECORD NO. OF GRIB RECORD IN INPUT FILE = ',krew
 
     ! COMPARE RECORD (GRIB) TO CONTROL CARD (PDS), THEY SHOULD MATCH
     pdsl(1:28)=grib(9:36)
     key = 2
     IF (.NOT.iw3pds(pdsl,pds,key)) THEN
        print 2900, nrec,(mova2i(pdsl(j)),j=1,28), &
             (mova2i(pds(j)),j=1,28)
2900    FORMAT ( 1x,i4,' (PDS) IN RECORD DOES NOT MATCH (PDS) IN ', 'CONTROL CARD ',/,7(1x,4z2.2), /,7(1x,4z2.2))
        cycle
     END IF
 
     ! Print PDS
     print 2, (mova2i(pdsl(j)),j=1,28)
2    FORMAT (' PDS = ',7(4z2.2,1x))
 
     ! Construct and print Description of GRIB field
     CALL w3fp11 (grib,pdsl,title,ier)
     IF (ier.NE.0) print *,'W3FP11 ERROR = ',ier
     print *,title(1:86)
 
     print *,' Size of GRIB Field = ',itot
 
     ! MAKE Flag Field Separator block
     call mkfldsep(csep,iopt,insize,itot+lenhead,lenout)
 
     ! MAKE WMO HEADER
     ! Get system date and time
     CALL makwmo (bulhed,kpds(10),kpds(11),kwbx,wmohdr)
 
     ! write out Separator block, Abbreviated WMO Heading,
     ! and GRIB field to output file.
     call wryte(lugo,lenout,csep)
     call wryte(lugo,lenhead,wmohdr)
     call wryte(lugo,itot,grib)
     nbul=nbul+1
  enddo foreachelement
 
  ! CLOSING SECTION
  IF (nbul .EQ. 0 ) THEN
     WRITE (6,fmt='('' SOMETHING WRONG WITH DATA CARDS...'', ''NOTHING WAS PROCESSED'')')
     stop 19
  ELSE
     CALL baclose (lugb,iret)
     CALL baclose (lugi,iret)
     CALL baclose (lugo,iret)
     WRITE (6,fmt='(//,'' ******** RECAP OF THIS EXECUTION '', ''********'',/,5X,''READ  '',I6,'' INDIVIDUAL IDS'', ' // &
          ' /,5X,''WROTE '',I6,'' BULLETINS OUT FOR TRANSMISSION'', //)') nrec, nbul
  ENDIF
 
  ! TEST TO SEE IF ANY BULLETINS MISSING
  mbul = nrec - nbul
  IF (mbul.NE.0) THEN
     print *,'BULLETINS MISSING = ',mbul
     stop 30
  END IF
 
  stop
END PROGRAM tocgrib
 
subroutine mkfldsep(csep,iopt,lenin,lenbull,lenout)
  !
  character*(*),intent(out) :: csep
  integer,intent(in) :: iopt,lenin,lenbull
  integer,intent(out) :: lenout
  !
  character(len=4),parameter :: cstar='****',clb='####'
  !
  if (iopt.eq.1) then
     if ( lenin .le. 18 .and. lenbull .le. 999999 ) then
        ! Create OPTION 1 separator block
        csep(1:4)=clb
        csep(5:7)='018'
        write(csep(8:13),fmt='(I6.6)') lenbull
        csep(14:17)=clb
        csep(18:18)=char(10)
        lenout=18
     else                         ! Create OPTION 1a separator block
        nnn=lenin
        if ( nnn.lt.23 ) nnn=23
        csep(1:4)=clb
        write(csep(5:7),fmt='(I3.3)') nnn
        write(csep(8:18),fmt='(I11.11)') lenbull
        csep(19:nnn-5)='0'
        csep(nnn-4:nnn-1)=clb
        csep(nnn:nnn)=char(10)
        lenout=nnn
     endif
  elseif (iopt.eq.2) then         !  Create OPTION 2 separator block
     csep(1:4)=cstar
     write(csep(5:14),fmt='(I10.10)') lenbull
     csep(15:18)=cstar
     csep(19:19)=char(10)
     lenout=19
  else
     print *,"mkfldsep: Option ",iopt," not recognized."
     csep(1:lenin)=' '
  endif
  !
  return
end subroutine mkfldsep