C NCLFORTSTART
      subroutine xy2pdf77(nxyz,x,y,z, xmsg,ymsg,zmsg, nby,mbx,kbz, pdf
     *                   ,mbxp1,nbyp1,kbzp1, binxbnd,binybnd,binzbnd
     *                   ,ipcnt,ier)
      implicit none
      integer  nxyz, nby,mbx,kbz, mbxp1,nbyp1,kbzp1, ipcnt,ier
      real x(nxyz), y(nxyz), z(nxyz), xmsg, ymsg, zmsg
     *    ,binxbnd(mbxp1), binybnd(nbyp1), binzbnd(kbzp1)
      real pdf(mbx,nby,kbz)
C NCLEND
      integer k, m, n, kb, mb, nb, kxyz, units, msgflg 

      ier = 0
      do kb=1,kbz
        do nb=1,nby
          do mb=1,mbx
             pdf(mb,nb,kb) = 0.0
          end do
        end do
      end do


      kxyz = 0
      do n=1,nxyz 
         if (x(n).eq.xmsg .or. y(n).eq.ymsg .or. 
     *       z(n).eq.zmsg) kxyz = kxyz+1
      end do

c Valid data check

      if (kxyz.eq.0) then
          ier = 1
          return
      end if

c Very large arrays can result in slow execution times. 
c .   To minimize unnecessary do loop checks, see if missing
c .   values are present.

      msgflg = 1
      if (kxyz.eq.nxyz) msgflg = 0

c Binning

      if (msgflg.gt.0) then
          do kb=1,kbz
            do nb=1,nby
              do mb=1,mbx
                 do n=1,nxyz 
                    if (x(n).ne.xmsg .and. y(n).ne.ymsg .and. 
     *                  z(n).ne.zmsg) then
                        if (x(n).ge.binxbnd(mb)  .and. 
     *                      x(n).lt.binxbnd(mb+1).and. 
     *                      y(n).ge.binybnd(nb)  .and.
     *                      y(n).lt.binybnd(nb+1).and. 
     *                      z(n).ge.binzbnd(kb)  .and.
     *                      z(n).lt.binzbnd(kb+1)) then 
                            pdf(mb,nb,kb) = pdf(mb,nb,kb) + 1.0
                        end if
                    end if
                 end do   
              end do     
            end do      
          end do       
      else
          do kb=1,kbz
            do nb=1,nby
              do mb=1,mbx
                do n=1,nxyz 
                   if (x(n).ge.binxbnd(mb)  .and. 
     *                 x(n).lt.binxbnd(mb+1).and. 
     *                 y(n).ge.binybnd(nb)  .and.
     *                 y(n).lt.binybnd(nb+1).and.
     *                 z(n).ge.binzbnd(kb)  .and.
     *                 z(n).lt.binzbnd(kb+1)) then 
                       pdf(mb,nb,kb) = pdf(mb,nb,kb) + 1.0
                   end if
                end do 
              end do   
            end do    
          end do      
      end if

c change return units to % if flag set

      units = 1
      if (ipcnt.eq.1) then
          units = 100
      end if

      do kb=1,kbz
        do nb=1,nby
          do mb=1,mbx
             pdf(mb,nb,kb) = units*(pdf(mb,nb,kb)/kxyz)    
          end do
        end do
      end do

      return
      end