CMAX: [Version 1.0 (6/22/93)]
CMAX: ----------------------------------------------------
CMAX: Transformation of SEARCH from search.f
CMAX: ----------------------------------------------------
      SUBROUTINE search(lu,ld,datalist,nreps,data_mm,nobs)
      
      INTEGER MaxObservPoints
      PARAMETER (MaxObservPoints = 150000)
C------------------------------------------------------------------C
C                                                                  !
C   Program Name: Search              Date:  May 16, 1992          !
C   Programmer  : J. Pfaendtner                                    !
C   Modifications: M. Seablom, D. Lamich                           !
C                                                                  !
C   NASA / GMSB  Optimum Interpolation Version 1.3                 !
C   ----------------------------------------------                 !
C                                                                  !
C  Performs data selection for all volumes in a vertical column.   !
C  Selection algorithm is designed to search for mostly rawinsondes!
C  first, then other data types at mandatory levels, then left     !
C  over data at significant levels.                                !
C                                                                  !
C                                                                  !
C------------------------------------------------------------------C
c
c     Maximal observation points
c
c
c
      PARAMETER (maxnrep = 30000, mxreps = 150, ntrig = 4, percentA = .2
     &, percentB = .2, percentS = .5)
      REAL alat (maxnrep), alon (maxnrep), data (MaxObservPoints * 10), 
     &datalist (mxreps, 50,NumMinVols), dheap (maxnrep), dist (maxnrep), 
     &pres (maxnrep), trigtab (450 * ntrig + 1)
      INTEGER iflag (maxnrep, 0 : 4), iheap (maxnrep), indexd (maxnrep),
     & itype (maxnrep), ixkxs (maxnrep), kheap (maxnrep), knt (0 : 4), m
     &arkoff (maxnrep), nreps (50), num (0 : 4)
      LOGICAL first, onestage, qprint
c     save first
      CHARACTER * 60 message
      COMMON / sparam / lsolve, ns, nr, nv, nh (13104), nq, kq (3), volc
     &entre_x (13104), ngv, volcentre_y (13104), glatup (13104), glonup 
     &(13104)
      COMMON / DATAPTRS / IKX, ILAT, ILON, IPRS, ITYP, IDEL, IFCE
      COMMON / scan / maxobs, xscan, onestage, MaxData
      COMMON / GRIDXX / IM, JNP, MLEV
      COMMON / OFFDATA / KOFFLV (200)
      COMMON / PMH / PMAND (18)
      COMMON / TRIGCOM / NTRIGS, TRIGCON, trigtab
!hpf$ distribute alat(block)
!hpf$ distribute alon(block)
!hpf$ distribute data(block)
!hpf$ distribute datalist(block,block)
!hpf$ distribute dheap(block)
!hpf$ distribute dist(block)
!hpf$ distribute glatup(block)
!hpf$ distribute glonup(block)
!hpf$ distribute iflag(block,block)
!hpf$ distribute iheap(block)
!hpf$ distribute indexd(block)
!hpf$ distribute itype(block)
!hpf$ distribute ixkxs(block)
!hpf$ distribute kheap(block)
!hpf$ distribute knt(block)
!hpf$ distribute kofflv(block)
!hpf$ distribute kq(block)
!hpf$ distribute markoff(block)
!hpf$ distribute nh(block)
!hpf$ distribute nreps(block)
!hpf$ distribute num(block)
!hpf$ distribute pmand(block)
!hpf$ distribute pres(block)
!hpf$ distribute trigtab(block)
!hpf$ distribute volcentre_x(block)
!hpf$ distribute volcentre_y(block)
      DATA first / .true. /
      INTEGER Ti1, Ti2
c      save rscan, rade, maxrep

      CALL getconh('RSCAN   ',rscan)
      CALL getconh('RADE    ',rade)
      CALL getconi('MAXREP  ',maxrep)
C  Gather all data within RSCAN kilometres of this data volume.

      CALL dgsort(data_mm,first,NumMinVols,nobs,nrep,maxnrep,
     &     rscan,rade,indexd,dist,alat,alon,itype,pres,ixkxs)

      where(nrep.eq.0) nreps ( lu : ld : nv, : ) = 0
      iwarn=any (mask=nrep .GT. maxnrep) 
      if(iwarn) THEN
        message = 'SEARCH: NREP exceeds MAXNREP.'
        CALL flash(5,message,1)
      ENDIF

C  Set MARKOFF = 1 for off-level and other secondary data types,
C  otherwise zero.

      forall(igp=1:NumMinVols,j=1:maxnrep,j.le.nrep(igp)) 
     &             markoff(igp,j)=kofflv(ixkxs(igp,j))

C  LOOP OVER THE GRID VOLUME LEVELS
      DO 7012 lev = lu , ld , nv
        nrepok(lev,:NumMinVols) = 0
        nreps (lev,:NumMinVols) = 0
        maxdat = maxrep * 3
        iflag(:,:,:) = 0

        IF (lev .LE. mlev) THEN
C  START UPPER AIR DATA SELECTION
          numcat = 4
          pabove = 0.0
          IF (lev .GT. 1) THEN
            pabove(:NumMinVols) = pmand(lev - 1)
          ENDIF
          pbelow(:NumMinVols) = 2000.0
          IF ((lev + nv) .LE. mlev) THEN
            pbelow(:NumMinVols) = pmand(lev + nv)
          ENDIF
C  Initialize flags for selecting upper air data.
C
C    0    for rawinsonde data
C    1    for data at pres = pabove with markoff = 0
C    2    for data at pres = pbelow with markoff = 0
C    3    for data between pres = pabove and  pres = pbelow   with
C         markoff = 0
C    4    for data between pres = pabove and  pres = pbelow   with
C         markoff = 1

C  Count and flag the available data.
          where(ixkxs.eq.7)
            where((pres.ge.pabove).and.(pres.le.pbelow))
                 iflag(:,0,:) = 1-markoff
            else
              where(pres.eq.pabove)
                 iflag(:,1,:) = 1-markoff
              else
                where(pres.eq.pbelow)
                   iflag(:,2,:) = 1-markoff
                else
                   where(pres.gt.pabove.and.pres.lt.pbelow)
                       iflag(:,3,:) = 1-markoff
                       iflag(:,4,:) = markoff
                   endwhere
                endwhere
              endwhere
            endwhere
          endwhere

          knt(0,:)= sum(iflag(:,0)*itype(:))
          knt(1,:)= sum(iflag(:,1)*itype(:))
          knt(2,:)= sum(iflag(:,2)*itype(:))
          knt(3,:)= sum(iflag(:,3)*itype(:))
          knt(4,:)= sum(iflag(:,4)*itype(:))

C Initially request 60% of the data to be rawinsondes.
          num(0,:) = min0(knt(0,:),maxdat)
          num(0,:) = num(0,:) * .6
          maxpa = (maxdat - num(0,:)) * percenta
          maxpb = (maxdat - num(0,:)) * percentb
          maxps = (maxdat - num(0,:)) * percents
          maxsn = (maxdat - num(0,:)) - maxpa - maxpb - maxps
          num(1,:) = min0(knt(1,:),maxpa)
          num(2,:) = min0(knt(2,:),maxpb)
          num(3,:) = min0(knt(3,:),maxps)
          num(4,:) = min0(knt(4,:),maxsn)
          ngot = num(1,:) + num(2,:) + num(3,:) + num(4,:) + num(0,:)

C  ADD SOME MORE DATA WITH MARKOFF=0 BETWEEN PABOVE AND PBELOW
          WHERE (ngot .LT. maxdat) 
            need = maxdat - ngot
            ngot = ngot - num(3,:)
            num(3,:) = min0(knt(3,:),need + num(3,:))
            ngot = ngot + num(3,:)
          ENDWHERE

C  ADD SOME MORE DATE WITH MARKOFF=0 AT PABOVE
          WHERE (ngot .LT. maxdat)
            need = maxdat - ngot
            ngot = ngot - num(1,:)
            num(1,:) = min0(knt(1,:),need + num(1,:))
            ngot = ngot + num(1,:)
          ENDWHERE

C  ADD SOME MORE DATA WITH MARKOFF=0 AT BBELOW
          WHERE (ngot .LT. maxdat)
            need = maxdat - ngot
            ngot = ngot - num(2,:)
            num(2,:) = min0(knt(2,:),need + num(2,:))
            ngot = ngot + num(2,:)
          ENDWHERE

C  ADD SOME MORE DATA WITH MARKOFF=1 BETWEEN PABOVE AND PBELOW
          WHERE (ngot .LT. maxdat) THEN
            need = maxdat - ngot
            num(4,:) = min0(knt(4,:),need + num(4,:))
          ENDWHERE

        ELSE

C************************************************
C  SEA LEVEL DATA SELECTION
C************************************************

          numcat = 1
          num(1) = 0
          num(1,:)= sum(iflag(:,1)*itype(:))
          num(1,:) = min0(num(1,:),maxdat)

        ENDIF

C  LOOP TO CHOOSE THE DATA FROM EACH CATEGORY
C  CHECK TO SEE IF THERE ARE ANY USABLE OBSERVATIONS IN CATEGORY ICAT
      prevcat=0
      do 280 icat = 0, numcat
          imask=iflag(:,icat,:)
          call cmf_scan_add(imask,imask,imask,1,
     &                      CMF_EXCLUSIVE,CMF_NONE,imask)
          nflag = maxval(imask,dim=1)
          maxnflag=maxval(nflag)
          forall(igp=1:NumMinVols,j=1:maxnrep,imask(igp,j).ne.0)
     &       kheap(igp,imask(igp,j))=j
          forall(igp=1:NumMinVols,j=1:maxnrep,imask(igp,j).ne.0)
     &       dheap(igp,imask(igp,j))=dist(igp,j)
          forall(igp=1:NumMinVols,j=1:maxnrep,imask(igp,j).ne.0)
     &       theap(igp,imask(igp,j))=itype(igp,j)
          forall(igp=1:NumMinVols,j=1:maxnrep,imask(igp,j).ne.0)
     &       xheap(igp,imask(igp,j))=indexd(igp,j)

          rmask=.false.
          forall(igp=1:NumMinVols,j=1:maxnflag,j.le.nflag(igp))
     &       rmask(igp,j)=.true.
          call cmf_order (ranking,dheap,2,rmask)

          forall(igp=1:NumMinVols,j=1:maxnflag,rmask(igp,j).ne.0)
     &       heapd(igp,ranking(igp,j))=dheap(igp,j)
          forall(igp=1:NumMinVols,j=1:maxnflag,rmask(igp,j).ne.0)
     &       heapk(igp,ranking(igp,j))=kheap(igp,j)
          forall(igp=1:NumMinVols,j=1:maxnflag,rmask(igp,j).ne.0)
     &       typei(igp,ranking(igp,j))=theap(igp,j)
          forall(igp=1:NumMinVols,j=1:maxnflag,rmask(igp,j).ne.0)
     &       heapx(igp,ranking(igp,j))=xheap(igp,j)

          call cmf_scan_add(tempcat(icat,:,:),typei,dummy,1,
     &                      CMF_UPWARD,CMF_EXCLUSIVE,CMF_NONE,rmask)
          forall(igp=1:NumMinVols,j=1:maxnflag, 
     &            tempcat(icat,igp,j).le.num(icat,igp)) 
            catmask(igp,j)=1
          call cmf_scan_add(catmask,catmask,dummy,1,
     &                      CMF_UPWARD,CMF_EXCLUSIVE,CMF_NONE,rmask)
          forall(igp=1:NumMinVols,j=1:maxnflag, 
     &            tempcat(icat,igp,j).le.num(icat,igp)) 
     &      datalist(j+prevcat(igp),lev,igp)=heapx(igp,j)
          prevcat=prevcat+maxval(catmask,dim=2)
  280     CONTINUE
        nreps(lev,:) = prevcat
 7012   CONTINUE
      RETURN
      END
