readwritesb.F90

Go to the documentation of this file.

 
 
recursive subroutine readsb(lunit,iret)
 
  use bufrlib
 
  use modv_vars, only: im8b
 
  use moda_msgcwd
  use moda_unptyp
  use moda_bitbuf
  use moda_bitmaps
  use moda_stcode
  use moda_borts
 
  implicit none
 
  integer, intent(in) :: lunit
  integer, intent(out) :: iret
  integer my_lunit, lun, il, im, ier, nbyt
 
  ! Check for I8 integers
 
  if(im8b) then
    im8b=.false.
    call x84(lunit,my_lunit,1)
    call readsb(my_lunit,iret)
    call x48(iret,iret,1)
    im8b=.true.
    return
  endif
 
  ! If we're catching bort errors, set a target return location if one doesn't already exist.
 
  if (bort_target_is_unset) then
    bort_target_is_unset = .false.
    caught_str_len = 0
    call catch_bort_readsb_c(lunit,iret)
    bort_target_is_unset = .true.
    return
  endif
 
  iret = -1
 
  ! Check the file status
 
  call status(lunit,lun,il,im)
  if(il==0) call bort('BUFRLIB: READSB - INPUT BUFR FILE IS CLOSED, IT MUST BE OPEN FOR INPUT')
  if(il>0) call bort('BUFRLIB: READSB - INPUT BUFR FILE IS OPEN FOR OUTPUT, IT MUST BE OPEN FOR INPUT')
  if(im==0) return
 
  ! See if there is another subset in the message
 
  if(nsub(lun)==msub(lun)) return
  nsub(lun) = nsub(lun) + 1
 
  ! Read the next subset and reset the pointers
 
  nbtm = 0
  lstnod = 0
  lstnodct = 0
  iscodes(lun) = 0
  linbtm = .false.
 
  if(msgunp(lun)==0) then
    ibit = mbyt(lun)*8
    call upb(nbyt,16,mbay(1,lun),ibit)
    call rdtree(lun,ier)
    if(ier/=0) return
    mbyt(lun) = mbyt(lun) + nbyt
  elseif(msgunp(lun)==1) then
    ! message with "standard" Section 3
    ibit = mbyt(lun)
    call rdtree(lun,ier)
    if(ier/=0) return
    mbyt(lun) = ibit
  else
    ! compressed message
    call rdcmps(lun)
    if (iscodes(lun) /= 0) return
  endif
 
  iret = 0
 
  return
end subroutine readsb
 
recursive integer function ireadsb(lunit) result(iret)
 
  use modv_vars, only: im8b
 
  implicit none
 
  integer, intent(in) :: lunit
  integer my_lunit
 
  ! Check for I8 integers.
 
  if(im8b) then
    im8b=.false.
    call x84(lunit,my_lunit,1)
    iret=ireadsb(my_lunit)
    im8b=.true.
    return
  endif
 
  call readsb(lunit,iret)
 
  return
end function ireadsb
 
recursive subroutine readns(lunit,subset,jdate,iret)
 
  use bufrlib
 
  use modv_vars, only: im8b, lendat
 
  use moda_msgcwd
  use moda_tables
  use moda_borts
 
  implicit none
 
  integer, intent(in) :: lunit
  integer, intent(out) :: jdate, iret
  integer my_lunit, lun, il, im
 
  character*8, intent(out) :: subset
  character*9 csubset
 
  ! Check for I8 integers
 
  if(im8b) then
    im8b=.false.
    call x84(lunit,my_lunit,1)
    call readns(my_lunit,subset,jdate,iret)
    call x48(jdate,jdate,1)
    call x48(iret,iret,1)
    im8b=.true.
    return
  endif
 
  ! If we're catching bort errors, set a target return location if one doesn't already exist.
 
  if (bort_target_is_unset) then
    bort_target_is_unset = .false.
    caught_str_len = 0
    call catch_bort_readns_c(lunit,csubset,jdate,len(csubset),iret)
    subset(1:8) = csubset(1:8)
    bort_target_is_unset = .true.
    return
  endif
 
  ! Refresh the subset and jdate parameters
 
  call status(lunit,lun,il,im)
  if(il==0) call bort('BUFRLIB: READNS - INPUT BUFR FILE IS CLOSED, IT MUST BE OPEN FOR INPUT')
  if(il>0) call bort('BUFRLIB: READNS - INPUT BUFR FILE IS OPEN FOR OUTPUT, IT MUST BE OPEN FOR INPUT')
  if(inode(lun)==0) then
    subset = '        '
  else
    subset = tag(inode(lun))(1:8)
  endif
  jdate = idate(lun)
  if (lendat/=10) jdate = mod(jdate,10**8)
 
  ! Read the next subset in the BUFR file
 
  do while (.true.)
    call readsb(lunit,iret)
    if (iret==0) exit
    call readmg(lunit,subset,jdate,iret)
    if (iret/=0) exit
  enddo
 
  return
end subroutine readns
 
recursive integer function ireadns(lunit,subset,idate) result(iret)
 
  use modv_vars, only: im8b
 
  implicit none
 
  integer, intent(in) :: lunit
  integer, intent(out) :: idate
  integer my_lunit
 
  character*8, intent(out) :: subset
 
  ! Check for I8 integers.
 
  if(im8b) then
     im8b=.false.
     call x84(lunit,my_lunit,1)
     iret=ireadns(my_lunit,subset,idate)
     call x48(idate,idate,1)
     im8b=.true.
     return
  endif
 
  call readns(lunit,subset,idate,iret)
 
  return
end function ireadns
 
recursive subroutine writsb(lunit)
 
  use modv_vars, only: im8b
 
  use moda_msgcmp
 
  implicit none
 
  integer, intent(in) :: lunit
  integer my_lunit, lun, il, im
 
  ! Check for I8 integers
 
  if(im8b) then
    im8b=.false.
 
    call x84 ( lunit, my_lunit, 1 )
    call writsb ( my_lunit )
 
    im8b=.true.
    return
  endif
 
  ! Check the file status
 
  call status(lunit,lun,il,im)
  if(il==0) call bort('BUFRLIB: WRITSB - OUTPUT BUFR FILE IS CLOSED, IT MUST BE OPEN FOR OUTPUT')
  if(il<0) call bort('BUFRLIB: WRITSB - OUTPUT BUFR FILE IS OPEN FOR INPUT, IT MUST BE OPEN FOR OUTPUT')
  if(im==0) call bort('BUFRLIB: WRITSB - A MESSAGE MUST BE OPEN IN OUTPUT BUFR FILE, NONE ARE')
 
  ! Pack up the subset and put it into the message
 
  call wrtree(lun)
  if( ccmf=='Y' ) then
    call wrcmps(lunit)
  else
    call msgupd(lunit,lun)
  endif
 
  return
end subroutine writsb
 
recursive subroutine writsa(lunxx,lmsgt,msgt,msgl)
 
  use modv_vars, only: im8b
 
  use moda_bufrmg
  use moda_msgcmp
 
  implicit none
 
  integer, intent(in) :: lunxx, lmsgt
  integer, intent(out) :: msgt(*), msgl
  integer my_lunxx, my_lmsgt, lunit, lun, il, im, n
 
  ! Check for I8 integers
 
  if(im8b) then
    im8b=.false.
 
    call x84 ( lunxx, my_lunxx, 1 )
    call x84 ( lmsgt, my_lmsgt, 1 )
    call writsa ( my_lunxx, my_lmsgt*2, msgt, msgl )
    msgl = msgl/2
    call x48 ( msgl, msgl, 1 )
 
    im8b=.true.
    return
  endif
 
  lunit = abs(lunxx)
 
  ! Check the file status
 
  call status(lunit,lun,il,im)
  if(il==0) call bort('BUFRLIB: WRITSA - OUTPUT BUFR FILE IS CLOSED, IT MUST BE OPEN FOR OUTPUT')
  if(il<0) call bort('BUFRLIB: WRITSA - OUTPUT BUFR FILE IS OPEN FOR INPUT, IT MUST BE OPEN FOR OUTPUT')
  if(im==0) call bort('BUFRLIB: WRITSA - A MESSAGE MUST BE OPEN IN OUTPUT BUFR FILE, NONE ARE')
 
  ! If lunxx < 0, force memory msg to be written (w/o any current subset)
 
  if(lunxx<0) call closmg(lunit)
 
  ! Is there a completed BUFR message to be returned?
 
  if(msglen(lun)>0) then
    if(msglen(lun)>lmsgt) call bort('BUFRLIB: WRITSA - OVERFLOW OF OUTPUT BUFR MESSAGE ARRAY; TRY A LARGER '// &
      'DIMENSION FOR THIS ARRAY')
    msgl = msglen(lun)
    do n=1,msgl
      msgt(n) = msgtxt(n,lun)
    enddo
    msglen(lun) = 0
  else
    msgl = 0
  endif
 
  if(lunxx<0) return
 
  ! Pack up the subset and put it into the message
 
  call wrtree(lun)
  if( ccmf=='Y' ) then
    call wrcmps(lunit)
  else
    call msgupd(lunit,lun)
  endif
 
  ! If the just-completed call to wrcmps() or msgupd() for this subset caused a message to be flushed to abs(lunxx), then
  ! attempt to retrieve and return that message now.  Otherwise, we run the risk that the next call to openmb() or openmg()
  ! might cause another message to be flushed, and thus overwrite the current message within array msgtxt before we
  ! had the chance to retrieve it during the next call to writsa().
 
  ! Also note that, in rare instances (e.g. if the byte count of the most recent subset is > 65530), we could end up with
  ! two BUFR messages available to be returned from this one call to writsa().  If sufficient space is available in the
  ! msgt array, then go ahead and return both messages now.
 
  if( (msglen(lun)>0) .and. (msgl+msglen(lun)<=lmsgt) ) then
    do n = 1,msglen(lun)
      msgt(msgl+n) = msgtxt(n,lun)
    enddo
    msgl = msgl+msglen(lun)
    msglen(lun) = 0
  endif
 
  return
end subroutine writsa
 
recursive subroutine rdmgsb(lunit,imsg,isub)
 
  use modv_vars, only: im8b
 
  use moda_msgcwd
  use moda_bitbuf
 
  implicit none
 
  integer, intent(in) :: lunit, imsg, isub
  integer my_lunit, my_imsg, my_isub, lun, il, im, i, jdate, iret
 
  character*128 bort_str
  character*8 subset
 
  ! Check for I8 integers
 
  if(im8b) then
    im8b=.false.
 
    call x84(lunit,my_lunit,1)
    call x84(imsg,my_imsg,1)
    call x84(isub,my_isub,1)
    call rdmgsb(my_lunit,my_imsg,my_isub)
 
    im8b=.true.
    return
  endif
 
  ! Open the file and skip to message #imsg
 
  call openbf(lunit,'IN',lunit)
  call status(lunit,lun,il,im)
 
  ! Note that we need to use subroutine readmg() to actually read in all of the messages (including the
  ! first (imsg-1) messages!), just in case there are any embedded dictionary messages in the file.
 
  do i=1,imsg
    call readmg(lunit,subset,jdate,iret)
    if(iret<0) then
      write(bort_str,'("BUFRLIB: RDMGSB - HIT END OF FILE BEFORE READING REQUESTED MESSAGE NO.",I5," IN '//&
        'BUFR FILE CONNECTED TO UNIT",I4)') imsg,lunit
      call bort(bort_str)
    endif
  enddo
 
  ! Position at subset #isub
 
  do i=1,isub
    call readsb(lunit,iret)
    if(iret<0) then
      write(bort_str,'("BUFRLIB: RDMGSB - ALL SUBSETS READ BEFORE READING REQ. SUBSET NO.",I3," IN '// &
        'REQ. MSG NO.",I5," IN BUFR FILE CONNECTED TO UNIT",I4)') isub,imsg,lunit
      call bort(bort_str)
    endif
  enddo
 
  return
end subroutine rdmgsb
 
subroutine msgupd(lunit,lun)
 
  use modv_vars, only: iprt, nby0, nby1, nby2, nby3
 
  use moda_msgcwd
  use moda_bitbuf
  use moda_h4wlc
 
  implicit none
 
  integer, intent(in) :: lunit, lun
  integer ibyt, lbyt, lbit, nbyt, ii, iupb
 
  logical msgfull
 
  character*128 errstr
 
  ! Pad the subset buffer
 
  call pad(ibay,ibit,ibyt,8)
 
  ! Check whether the new subset should be written into the currently open message
 
  if(msgfull(mbyt(lun),ibyt,maxbyt) .or. ((ibyt>65530).and.(nsub(lun)>0))) then
    ! No it should not, either because it doesn't fit
    !      OR
    ! It has byte count > 65530 (sufficiently close to the upper limit for the 16 bit byte counter placed at the beginning
    ! of each subset), and the current message has at least one subset in it
    !
    ! In either of these cases, we need to write out the current message and then create a new one to hold the current subset
    call msgwrt(lunit,mbay(1,lun),mbyt(lun))
    call msgini(lun)
  endif
 
  if(msgfull(mbyt(lun),ibyt,maxbyt)) then
    ! This is an overlarge subset that won't fit in any message given the current value of maxbyt, so discard the subset
    ! and exit gracefully.
    if(iprt>=0) then
      call errwrt('+++++++++++++++++++++WARNING+++++++++++++++++++++++')
      write ( unit=errstr, fmt='(A,A,I7,A)') 'BUFRLIB: MSGUPD - SUBSET LONGER THAN ANY POSSIBLE MESSAGE ', &
        '{MAXIMUM MESSAGE LENGTH = ', maxbyt, '}'
      call errwrt(errstr)
      call errwrt('>>>>>>>OVERLARGE SUBSET DISCARDED FROM FILE<<<<<<<<')
      call errwrt('+++++++++++++++++++++WARNING+++++++++++++++++++++++')
      call errwrt(' ')
    endif
    call usrtpl(lun,1,1)
    return
  endif
 
  ! Set a byte count and transfer the subset buffer into the message
 
  lbit = 0
  call pkb(ibyt,16,ibay,lbit)
 
  ! Note that we want to append the data for this subset to the end of Section 4, but the value in mbyt(lun) already includes
  ! the length of Section 5 (i.e. 4 bytes).  Therefore, we need to begin writing at the point 3 bytes prior to the byte
  ! currently pointed to by mbyt(lun).
 
  call mvb(ibay,1,mbay(1,lun),mbyt(lun)-3,ibyt)
 
  ! Update the subset and byte counters
 
  mbyt(lun) = mbyt(lun) + ibyt
  nsub(lun) = nsub(lun) + 1
 
  lbit = (nby0+nby1+nby2+4)*8
  call pkb(nsub(lun),16,mbay(1,lun),lbit)
 
  lbyt = nby0+nby1+nby2+nby3
  nbyt = iupb(mbay(1,lun),lbyt+1,24)
  lbit = lbyt*8
  call pkb(nbyt+ibyt,24,mbay(1,lun),lbit)
 
  ! If any long character strings are being held internally for storage into this subset, store them now
 
  if(nh4wlc>0) then
    do ii = 1, nh4wlc
      call writlc(luh4wlc(ii),chh4wlc(ii),sth4wlc(ii))
    enddo
    nh4wlc = 0
  endif
 
  ! If the subset byte count is > 65530, then give it its own one-subset message (cannot have any other subsets in this
  ! message because their beginning would be beyond the upper limit of 65535 in the 16-bit byte counter, meaning they
  ! could not be located!)
 
  if(ibyt>65530) then
    if(iprt>=1) then
      call errwrt('+++++++++++++++++++++WARNING+++++++++++++++++++++++')
      write ( unit=errstr, fmt='(A,I7,A,A)') 'BUFRLIB: MSGUPD - SUBSET HAS BYTE COUNT = ',ibyt,' > UPPER LIMIT OF 65535'
      call errwrt(errstr)
      call errwrt('>>>>>>>WILL BE WRITTEN INTO ITS OWN MESSAGE<<<<<<<<')
      call errwrt('+++++++++++++++++++++WARNING+++++++++++++++++++++++')
      call errwrt(' ')
    endif
    call msgwrt(lunit,mbay(1,lun),mbyt(lun))
    call msgini(lun)
  endif
 
  ! Reset the user arrays
 
  call usrtpl(lun,1,1)
 
  return
end subroutine msgupd
 
subroutine pad(ibay,ibit,ibyt,ipadb)
 
  implicit none
 
  integer, intent(inout) :: ibay(*), ibit
  integer, intent(in) :: ipadb
  integer, intent(out) :: ibyt
  integer ipad
 
  character*128 bort_str
 
  ! Pad the subset to an ipadb bit boundary
 
  ipad = ipadb - mod(ibit+8,ipadb)
  ! First pack the # of bits being padded (this is a delayed replication factor)
  call pkb(ipad,8,ibay,ibit)
  ! Now pad with zeroes to the byte boundary
  call pkb(0,ipad,ibay,ibit)
  ibyt = ibit/8
 
  if(mod(ibit,8)/=0) then
    write(bort_str,'("BUFRLIB: PAD - THE NUMBER OF BITS IN A PACKED'// &
      ' SUBSET AFTER PADDING (",I8,") IS NOT A MULTIPLE OF 8")') ibit
    call bort(bort_str)
  endif
 
  return
end subroutine pad
 
recursive integer function lcmgdf(lunit,subset) result(iret)
 
  use modv_vars, only: im8b
 
  use moda_tables
 
  implicit none
 
  integer, intent(in) :: lunit
  integer my_lunit, lun, il, im, mtyp, msbt, inod, nte, i
 
  character*8, intent(in) :: subset
 
  ! Check for I8 integers.
 
  if(im8b) then
    im8b=.false.
 
    call x84(lunit,my_lunit,1)
    iret=lcmgdf(my_lunit,subset)
 
    im8b=.true.
    return
  endif
 
  iret = 0
 
  ! Get lun from lunit.
 
  call status(lunit,lun,il,im)
  if (il==0) call bort('BUFRLIB: LCMGDF - INPUT BUFR FILE IS CLOSED, IT MUST BE OPEN')
 
  ! Confirm that subset is defined for this logical unit.
 
  call nemtba(lun,subset,mtyp,msbt,inod)
 
  ! Check if there's a long character string in the definition.
 
  nte = isc(inod)-inod
 
  do i = 1, nte
    if ( (typ(inod+i)=='CHR') .and. (ibt(inod+i)>64) ) then
      iret = 1
      return
    endif
  enddo
 
  iret = 0
 
  return
end function lcmgdf
 
recursive subroutine ufbpos(lunit,irec,isub,subset,jdate)
 
  use bufrlib
 
  use modv_vars, only: im8b
 
  use moda_msgcwd
  use moda_bitbuf
 
  implicit none
 
  integer, intent(in) :: lunit, irec, isub
  integer, intent(out) :: jdate
  integer my_lunit, my_irec, my_isub, lun, il, im, jrec, jsub, iret
 
  character*128 bort_str
  character*8, intent(out) :: subset
 
  ! Check for I8 integers
 
  if(im8b) then
    im8b=.false.
    call x84(lunit,my_lunit,1)
    call x84(irec,my_irec,1)
    call x84(isub,my_isub,1)
    call ufbpos(my_lunit,my_irec,my_isub,subset,jdate)
    call x48(jdate,jdate,1)
    im8b=.true.
    return
  endif
 
  !  Make sure a file is open for input
 
  call status(lunit,lun,il,im)
  if(il==0) call bort('BUFRLIB: UFBPOS - INPUT BUFR FILE IS CLOSED, IT MUST BE OPEN FOR INPUT')
  if(il>0) call bort('BUFRLIB: UFBPOS - INPUT BUFR FILE IS OPEN FOR OUTPUT, IT MUST BE OPEN FOR INPUT')
 
  if(irec<=0) then
    write(bort_str,'("BUFRLIB: UFBPOS - REQUESTED MESSAGE NUMBER TO READ IN (",I5,") IS NOT VALID")') irec
    call bort(bort_str)
  endif
  if(isub<=0) then
    write(bort_str,'("BUFRLIB: UFBPOS - REQUESTED SUBSET NUMBER TO READ IN (",I5,") IS NOT VALID")') isub
    call bort(bort_str)
  endif
 
  !  See where pointers are currently located
 
  call ufbcnt(lunit,jrec,jsub)
 
  !  Rewind file if requested pointers are behind current pointers
 
  if(irec<jrec .or. (irec==jrec.and.isub<jsub)) then
    call cewind_c(lun)
    nmsg(lun) = 0
    nsub(lun) = 0
    call ufbcnt(lunit,jrec,jsub)
  endif
 
  ! Read subset #isub from message #irec from file
 
  do while (irec>jrec)
    call readmg(lunit,subset,jdate,iret)
    if(iret<0) then
      write(bort_str,'("BUFRLIB: UFBPOS - REQUESTED MESSAGE NUMBER '// &
        'TO READ IN (",I5,") EXCEEDS THE NUMBER OF MESSAGES IN THE FILE (",I5,")")') irec, jrec
      call bort(bort_str)
    endif
    call ufbcnt(lunit,jrec,jsub)
  enddo
 
  do while (isub>jsub)
    call readsb(lunit,iret)
    if(iret/=0) then
      write(bort_str,'("BUFRLIB: UFBPOS - REQ. SUBSET NUMBER TO READ'// &
        ' IN (",I5,") EXCEEDS THE NUMBER OF SUBSETS (",I5,") IN THE REQ. MESSAGE (",I5,")")') isub, jsub, irec
      call bort(bort_str)
    endif
    call ufbcnt(lunit,jrec,jsub)
  enddo
 
  return
end subroutine ufbpos
 
subroutine rdtree(lun,iret)
 
  use modv_vars, only: bmiss
 
  use moda_usrint
  use moda_usrbit
  use moda_ival
  use moda_bitbuf
  use moda_tables
 
  implicit none
 
  integer, intent(in) :: lun
  integer, intent(out) :: iret
  integer ier, n, node, kbit, nbt, icbfms, igetrfel
 
  character*8 cval
 
  real*8 rval, ups
 
  equivalence(cval,rval)
 
  iret = 0
 
  ! Cycle through a subset setting up the template
 
  mbit(1) = ibit
  nbit(1) = 0
  call rcstpl(lun,ier)
  if(ier/=0) then
    iret = -1
    return
  endif
 
  ! Loop through each element of the subset, unpacking each value and then converting it to the proper type
 
  do n=1,nval(lun)
    call upb8(ival(n),nbit(n),mbit(n),mbay(1,lun))
    node = inv(n,lun)
    if(itp(node)==1) then
      ! The unpacked value is a delayed descriptor replication factor.
      val(n,lun) = ival(n)
    elseif(itp(node)==2) then
      ! The unpacked value is a real.
      nrfelm(n,lun) = igetrfel(n,lun)
      if (ival(n)<2_8**ibt(node)-1) then
        val(n,lun) = ups(ival(n),node)
      else
        val(n,lun) = bmiss
      endif
    elseif(itp(node)==3) then
      ! The value is a character string, so unpack it using an equivalenced real*8 value.  Note that a maximum of 8 characters
      ! will be unpacked here, so a separate subsequent call to subroutine readlc() will be needed to fully unpack any string
      ! longer than 8 characters.
      cval = ' '
      kbit = mbit(n)
      nbt = min(8,nbit(n)/8)
      call upc(cval,nbt,mbay(1,lun),kbit,.true.)
      if (nbit(n)<=64 .and. icbfms(cval,nbt)/=0) then
        val(n,lun) = bmiss
      else
        val(n,lun) = rval
      endif
    endif
  enddo
 
  ibit = nbit(nval(lun))+mbit(nval(lun))
 
  return
end subroutine rdtree
 
subroutine wrtree(lun)
 
  use moda_usrint
  use moda_ival
  use moda_ufbcpl
  use moda_bitbuf
  use moda_tables
 
  implicit none
 
  integer, intent(in) :: lun
  integer*8 ipks
  integer n, node, nbit, ncr, numchr, jj, ibfms, igetrfel, imrkopr
 
  character*120 lstr
  character*8 cval
 
  real*8 rval
 
  equivalence(cval,rval)
 
  ! Convert user numbers into scaled integers
 
  do n=1,nval(lun)
    node = inv(n,lun)
    nrfelm(n,lun) = igetrfel(n,lun)
    if(itp(node)==1) then
      ival(n) = nint(val(n,lun))
    elseif(typ(node)=='NUM') then
      if( (ibfms(val(n,lun))==1) .or. (val(n,lun)/=val(n,lun)) ) then
        ! The user number is either "missing" or NaN.
        ival(n) = -1
      else
        ival(n) = ipks(val(n,lun),node)
      endif
      call strbtm(n,lun,int(ival(n)))
    endif
  enddo
 
  ! Pack the user array into the subset buffer
 
  ibit = 16
 
  do n=1,nval(lun)
    node = inv(n,lun)
    if(itp(node)<3) then
      ! The value to be packed is numeric.
      if ( imrkopr(tag(node)) == 1 ) then
        nbit = ibt(inv(nrfelm(n,lun),lun))
      else
        nbit = ibt(node)
      endif
      call pkb8(ival(n),nbit,ibay,ibit)
    else
      ! The value to be packed is a character string.
      ncr=ibt(node)/8
      if ( ncr>8 .and. luncpy(lun)/=0 ) then
        ! The string is longer than 8 characters and there was a preceeding call to ufbcpy() involving this output unit,
        ! so read the long string with readlc() and then write it into the output buffer using pkc().
        call readlc(luncpy(lun),lstr,tag(node))
        call pkc(lstr,ncr,ibay,ibit)
      else
        rval = val(n,lun)
        if(ibfms(rval)/=0) then
          ! The value is "missing", so set all bits to 1 before packing the field as a character string.
          numchr = min(ncr,len(lstr))
          do jj = 1, numchr
            call ipkm(lstr(jj:jj),1,255)
          enddo
          call pkc(lstr,numchr,ibay,ibit)
        else
          ! The value is not "missing", so pack the equivalenced character string.  Note that a maximum of 8 characters
          ! will be packed here, so a separate subsequent call to subroutine writlc() will be needed to fully encode any
          ! string longer than 8 characters.
          call pkc(cval,ncr,ibay,ibit)
        endif
      endif
    endif
  enddo
 
  ! Reset ufbcpy() file pointer
 
  luncpy(lun)=0
 
  return
end subroutine wrtree
 
subroutine rcstpl(lun,iret)
 
  use modv_vars, only: maxjl, maxss, maxrcr, iprt
 
  use moda_usrint
  use moda_usrbit
  use moda_msgcwd
  use moda_bitbuf
  use moda_tables
  use moda_usrtmp
 
  implicit none
 
  character*128 bort_str
 
  integer, intent(in) :: lun
  integer, intent(out) :: iret
  integer nbmp(2,maxrcr), newn(2,maxrcr), knx(maxrcr), nodi, node, mbmp, nr, i, j, n, nn, n1, n2, new, ivob, igetrfel
 
  iret = 0
 
  ! Set the initial values for the template
 
  inv(1,lun) = inode(lun)
  val(1,lun) = 0
  nbmp(1,1) = 1
  nbmp(2,1) = 1
  nodi = inode(lun)
  node = inode(lun)
  mbmp = 1
  nval(lun) = 1
  nr = 0
  knx(1:maxrcr) = 0
 
  outer: do while (.true.)
 
    ! Set up the parameters for a level of recursion
 
    nr = nr+1
    if(nr>maxrcr) then
      write(bort_str,'("BUFRLIB: RCSTPL - THE NUMBER OF RECURSION LEVELS EXCEEDS THE LIMIT (",I3,")")') maxrcr
      call bort(bort_str)
    endif
    nbmp(1,nr) = 1
    nbmp(2,nr) = mbmp
 
    n1 = iseq(node,1)
    n2 = iseq(node,2)
    if(n1==0) then
      write(bort_str,'("BUFRLIB: RCSTPL - UNSET EXPANSION SEGMENT ",A)') tag(nodi)
      call bort(bort_str)
    endif
    if(n2-n1+1>maxjl) then
      if(iprt>=0) then
        call errwrt('++++++++++++++BUFR ARCHIVE LIBRARY+++++++++++++++++')
        call errwrt('BUFRLIB: RCSTPL - MAXJL OVERFLOW; SUBSET SKIPPED')
        call errwrt('++++++++++++++BUFR ARCHIVE LIBRARY+++++++++++++++++')
      endif
      iret = -1
      return
    endif
    newn(1,nr) = 1
    newn(2,nr) = n2-n1+1
 
    do n=1,newn(2,nr)
      nn = jseq(n+n1-1)
      iutmp(n,nr) = nn
      vutmp(n,nr) = vali(nn)
    enddo
 
    do while (.true.)
 
      ! Store nodes at some recursion level
 
      do i=nbmp(1,nr),nbmp(2,nr)
        if(knx(nr)==0) knx(nr) = nval(lun)
        if(i>nbmp(1,nr)) newn(1,nr) = 1
        do j=newn(1,nr),newn(2,nr)
          if(nval(lun)+1>maxss) then
            if(iprt>=0) then
              call errwrt('++++++++++++++BUFR ARCHIVE LIBRARY+++++++++++++++++')
              call errwrt('BUFRLIB: RCSTPL - MAXSS OVERFLOW; SUBSET SKIPPED')
              call errwrt('++++++++++++++BUFR ARCHIVE LIBRARY+++++++++++++++++')
            endif
            iret = -1
            return
          endif
          nval(lun) = nval(lun)+1
          node = iutmp(j,nr)
          ! inv is positional index in internal jump/link table for packed subset element nval(lun) in mbay
          inv(nval(lun),lun) = node
          ! mbit is the bit in mbay pointing to where the packed subset element nval(lun) begins
          mbit(nval(lun)) = mbit(nval(lun)-1)+nbit(nval(lun)-1)
          ! nbit is the number of bits in mbay occupied by packed subset element nval(lun)
          nrfelm(nval(lun),lun) = igetrfel(nval(lun),lun)
          nbit(nval(lun)) = ibt(node)
          if(nbit(nval(lun))==1) then
            ! Check whether this is a bitmap entry
            call upbb(ivob,nbit(nval(lun)),mbit(nval(lun)),mbay(1,lun))
            call strbtm(nval(lun),lun,ivob)
          endif
          ! Actual unpacked subset values are initialized here
          val(nval(lun),lun) = vutmp(j,nr)
          if(itp(node)==1) then
            call upbb(mbmp,nbit(nval(lun)),mbit(nval(lun)),mbay(1,lun))
            newn(1,nr) = j+1
            nbmp(1,nr) = i
            cycle outer
          endif
        enddo
        new = nval(lun)-knx(nr)
        val(knx(nr)+1,lun) = val(knx(nr)+1,lun) + new
        knx(nr) = 0
      enddo
 
      ! Check if we need to continue one recursion level back
 
      if(nr-1 == 0) exit outer
      nr = nr-1
    enddo
 
  enddo outer
 
  return
end subroutine rcstpl
 
subroutine usrtpl(lun,invn,nbmp)
 
  use modv_vars, only: maxjl, maxss, iprt
 
  use moda_usrint
  use moda_msgcwd
  use moda_tables
  use moda_ivttmp
  use moda_stcode
 
  implicit none
 
  integer, intent(in) :: lun, invn, nbmp
  integer i, j, ival, jval, n, n1, n2, nodi, node, newn, invr, knvn
 
  character*128 bort_str, errstr
 
  logical drp, drs, drb, drx
 
  if(iprt>=2) then
    call errwrt('++++++++++++++BUFR ARCHIVE LIBRARY+++++++++++++++++')
    write ( unit=errstr, fmt='(A,I3,A,I7,A,I5,A,A10)' ) &
      'BUFRLIB: USRTPL - LUN:INVN:NBMP:TAG(INODE(LUN)) = ', lun, ':', invn, ':', nbmp, ':', tag(inode(lun))
    call errwrt(errstr)
    call errwrt('++++++++++++++BUFR ARCHIVE LIBRARY+++++++++++++++++')
    call errwrt(' ')
  endif
 
  if(nbmp<=0) then
    if(iprt>=1) then
      call errwrt('+++++++++++++++++++++WARNING+++++++++++++++++++++++')
      call errwrt(.LE.'BUFRLIB: USRTPL - NBMP  0 - IMMEDIATE RETURN')
      call errwrt('+++++++++++++++++++++WARNING+++++++++++++++++++++++')
      call errwrt(' ')
    endif
    return
  endif
 
  drp = .false.
  drs = .false.
  drx = .false.
 
  ! Set up a node expansion
 
  if(invn==1) then
    ! The node is a Table A mnemonic
    nodi = inode(lun)
    inv(1,lun) = nodi
    nval(lun) = 1
    if(nbmp/=1) then
      write(bort_str,'("BUFRLIB: USRTPL - THIRD ARGUMENT (INPUT) = ",'// &
        'I4,", MUST BE 1 WHEN SECOND ARGUMENT (INPUT) IS 1 (SUBSET NODE) (",A,")")') nbmp, tag(nodi)
      call bort(bort_str)
    endif
  elseif(invn>0 .and. invn<=nval(lun)) then
    ! The node is (hopefully) a delayed replication factor
    nodi = inv(invn,lun)
    drp = typ(nodi) == 'DRP'
    drs = typ(nodi) == 'DRS'
    drb = typ(nodi) == 'DRB'
    drx = drp .or. drs .or. drb
    ival = nint(val(invn,lun))
    jval = 2**ibt(nodi)-1
    val(invn,lun) = ival+nbmp
    if(drb.and.nbmp/=1) then
      write(bort_str,'("BUFRLIB: USRTPL - THIRD ARGUMENT (INPUT) = ",'// &
        'I4,", MUST BE 1 WHEN NODE IS DRB (1-BIT DELAYED REPL. FACTOR) (",A,")")') nbmp, tag(nodi)
      call bort(bort_str)
    endif
    if(.not.drx) then
      write(bort_str,'("BUFRLIB: USRTPL - NODE IS OF TYPE ",A," - IT '// &
        'MUST BE EITHER A SUBSET OR DELAYED REPL. FACTOR (",A,")")') typ(nodi), tag(nodi)
      call bort(bort_str)
    endif
    if(ival<0) then
      write(bort_str,'("BUFRLIB: USRTPL - REPLICATION FACTOR IS NEGATIVE (=",I5,") (",A,")")') ival, tag(nodi)
      call bort(bort_str)
    endif
    if(ival+nbmp>jval) then
      write(bort_str,'("BUFRLIB: USRTPL - REPLICATION FACTOR OVERFLOW (EXCEEDS MAXIMUM OF",I6," (",A,")")') jval, tag(nodi)
      call errwrt(bort_str)
      iscodes(lun) = 1
      return
    endif
  else
    write(bort_str,'("BUFRLIB: USRTPL - INVENTORY INDEX {FIRST '// &
      'ARGUMENT (INPUT)} OUT OF BOUNDS (=",I5,", RANGE IS 1 TO",I6,") ")') invn, nval(lun)
    call bort(bort_str)
  endif
 
  ! Recall a pre-fab node expansion segment
 
  newn = 0
  n1 = iseq(nodi,1)
  n2 = iseq(nodi,2)
 
  if(n1==0) then
    write(bort_str,'("BUFRLIB: USRTPL - UNSET EXPANSION SEGMENT (",A,")")') tag(nodi)
    call bort(bort_str)
  endif
  if(n2-n1+1>maxjl) then
    write(bort_str,'("BUFRLIB: USRTPL - TEMPLATE ARRAY OVERFLOW, EXCEEDS THE LIMIT (",I6,") (",A,")")') maxjl, tag(nodi)
    call bort(bort_str)
  endif
 
  do n=n1,n2
    newn = newn+1
    itmp(newn) = jseq(n)
    vtmp(newn) = vali(jseq(n))
  enddo
 
  ! Move old nodes and store new ones
 
  if(nval(lun)+newn*nbmp>maxss) then
    write(bort_str,'("BUFRLIB: USRTPL - INVENTORY OVERFLOW (",I6,"), EXCEEDS THE LIMIT (",I6,") (",A,")")') &
      nval(lun)+newn*nbmp, maxss, tag(nodi)
    call bort(bort_str)
  endif
 
  do j=nval(lun),invn+1,-1
    inv(j+newn*nbmp,lun) = inv(j,lun)
    val(j+newn*nbmp,lun) = val(j,lun)
  enddo
 
  if(drp.or.drs) vtmp(1) = newn
  knvn = invn
 
  do i=1,nbmp
    do j=1,newn
      knvn = knvn+1
      inv(knvn,lun) = itmp(j)
      val(knvn,lun) = vtmp(j)
    enddo
  enddo
 
  ! Reset pointers and counters
 
  nval(lun) = nval(lun) + newn*nbmp
 
  if(iprt>=2) then
    call errwrt('++++++++++++++BUFR ARCHIVE LIBRARY+++++++++++++++++')
    write ( unit=errstr, fmt='(A,A,A10,2(A,I5),A,I7)' ) 'BUFRLIB: USRTPL - TAG(INV(INVN,LUN)):NEWN:NBMP:', &
      'NVAL(LUN) = ', tag(inv(invn,lun)), ':', newn, ':', nbmp, ':', nval(lun)
    call errwrt(errstr)
    do i=1,newn
      write ( unit=errstr, fmt='(2(A,I5),A,A10)' ) 'For I = ', i, ', ITMP(I) = ', itmp(i), ', TAG(ITMP(I)) = ', tag(itmp(i))
      call errwrt(errstr)
    enddo
    call errwrt('++++++++++++++BUFR ARCHIVE LIBRARY+++++++++++++++++')
    call errwrt(' ')
  endif
 
  if(drx) then
    node = nodi
    invr = invn
    outer: do while (.true.)
      node = jmpb(node)
      if(node<=0) exit
      if(itp(node)==0) then
        do invr=invr-1,1,-1
          if(inv(invr,lun)==node) then
            val(invr,lun) = val(invr,lun)+newn*nbmp
            cycle outer
          endif
        enddo
        write(bort_str,'("BUFRLIB: USRTPL - BAD BACKUP STRATEGY (",A,")")') tag(nodi)
        call bort(bort_str)
      else
        cycle
      endif
    enddo outer
  endif
 
  return
end subroutine usrtpl
 
recursive subroutine invmrg(lubfi,lubfj)
 
  use modv_vars, only: im8b
 
  use moda_usrint
  use moda_tables
  use moda_mrgcom
 
  implicit none
 
  integer, intent(in) :: lubfi, lubfj
  integer my_lubfi, my_lubfj, luni, il, im, lunj, jl, jm, is, js, node, nodj, ityp, iwrds, jwrds, &
    n, ioff, nwords, ibfms
 
  character*128 bort_str
 
  logical herei, herej, missi, missj, samei
 
  ! Check for I8 integers
 
  if(im8b) then
    im8b=.false.
    call x84(lubfi,my_lubfi,1)
    call x84(lubfj,my_lubfj,1)
    call invmrg(my_lubfi,my_lubfj)
    im8b=.true.
    return
  endif
 
  is = 1
  js = 1
 
  ! Get the unit pointers
 
  call status(lubfi,luni,il,im)
  call status(lubfj,lunj,jl,jm)
 
  ! Step through the buffers comparing the inventory and merging data
 
  do while(is<=nval(luni))
    ! Confirm we're at the same node in each buffer
    node = inv(is,luni)
    nodj = inv(js,lunj)
    if(node/=nodj) then
      write(bort_str,'("BUFRLIB: INVMRG - NODE FROM INPUT BUFR FILE '// &
        '(",I7,") DOES NOT EQUAL NODE FROM OUTPUT BUFR FILE (",I7,"), TABULAR MISMATCH")') node, nodj
      call bort(bort_str)
    endif
 
    ityp = itp(node)
    if(ityp==1) then
      ! Do an entire sequence replacement
      if(typ(node)=='DRB') then
        ioff = 0
      else
        ioff = 1
      endif
      iwrds = nwords(is,luni)+ioff
      jwrds = nwords(js,lunj)+ioff
      if(iwrds>ioff .and. jwrds==ioff) then
        do n=nval(lunj),js+1,-1
          inv(n+iwrds-jwrds,lunj) = inv(n,lunj)
          val(n+iwrds-jwrds,lunj) = val(n,lunj)
        enddo
        do n=0,iwrds
          inv(js+n,lunj) = inv(is+n,luni)
          val(js+n,lunj) = val(is+n,luni)
        enddo
        nval(lunj) = nval(lunj)+iwrds-jwrds
        jwrds = iwrds
        nrpl = nrpl+1
      endif
      is = is+iwrds
      js = js+jwrds
    elseif((ityp==2).or.(ityp==3)) then
      ! Fill missing values
      herei = ibfms(val(is,luni))==0
      herej = ibfms(val(js,lunj))==0
      missi = .not.(herei)
      missj = .not.(herej)
      samei = val(is,luni)==val(js,lunj)
      if(herei.and.missj) then
        val(js,lunj) = val(is,luni)
        nmrg = nmrg+1
      elseif(herei.and.herej.and..not.samei) then
        namb = namb+1
      endif
    endif
 
    ! Bump the counters and go check the next pair
    is = is + 1
    js = js + 1
  enddo
 
  ntot = ntot+1
 
  return
end subroutine invmrg
 
integer function nwords(n,lun) result(iret)
 
  use moda_usrint
 
  implicit none
 
  integer, intent(in) :: n, lun
  integer k
 
  iret = 0
 
  do k=1,nint(val(n,lun))
    iret = iret + nint(val(iret+n+1,lun))
  enddo
 
  return
end function nwords