      SUBROUTINE PUV_PREP_COV(mv_begin,mv_end,NumMinVols,
     &             nobs,LEVLINDX,nrep,data,
     &             FCERR,glatup,glonup,nh,volcentre_x,volcentre_y,
     &             gplat,gplon,plat,plon,pstd,gstd,pdel,pkx,wlat,
     &             wlon,wstd,udel,vdel,wkx,glat,glon,np,nw)

c...................................................................
! PUV_PREP_COV gathers data indexed by LEVLINDX into arrays indexed
! by mini-volume label, which will be shipped to local PUVCOV routines.
c..................................................................c
!                                                                  !
! NASA Goddard Space Flight Center       Data Assimilation Office  !
!                                                                  !
! Goddard Earth Observing System                                   !
! Data Analysis System - Version 1.3.2                             !
! September, 1993                                                  !
!                                                                  !
c------------------------------------------------------------------c
C**********************************************************************
C  Data parallel implementation (CM Fortran) : M. Makivic (NPAC)      *
C  Date:       4/94                                                   *
C**********************************************************************

      include 'params.h'
      PARAMETER (maxdat = 150)

      real,dimension(NumMinVols,mxreps):: LEVLINDX,typemm
!hpf$ distribute levlindx(block,block)
!hpf$ distribute typemm(block,block)
      integer,dimension(NumMinVols)::nrep
!hpf$ distribute nrep(block)
      real,dimension(im,jnp):: fcerr     
!hpf$ distribute fcerr(block,block)
      real,dimension(nobs,NRECORDS):: data
!hpf$ distribute data(block,block)
 
      integer mv_begin, mv_end
      integer,dimension(NumMinVols,700):: gplat,gplon
!hpf$ distribute gplat(block,block)
!hpf$ distribute gplon(block,block)

      real,dimension(NumMinVols,700):: plat,plon,pstd,
     &     gstd,pdel,pkx,wlat,wlon,wstd,wkx,
     &     glat,glon,plstmm,wlstmm,udel,vdel

!hpf$ distribute     udel(block,block)
!hpf$ distribute     vdel(block,block)
!hpf$ distribute     plat(block,block)
!hpf$ distribute     plon(block,block)
!hpf$ distribute     pstd(block,block)
!hpf$ distribute     gstd(block,block)
!hpf$ distribute     pdel(block,block)
!hpf$ distribute     pkx(block,block)
!hpf$ distribute     wlat(block,block)
!hpf$ distribute     wlon(block,block)
!hpf$ distribute     wstd(block,block)
!hpf$ distribute     wkx(block,block)
!hpf$ distribute     plstmm(block,block)
!hpf$ distribute     wlstmm(block,block)
!hpf$ distribute     glat(block,block)
!hpf$ distribute     glon(block,block)

      integer,dimension(NumMinVols)::np,nw
!hpf$ distribute np(block)
!hpf$ distribute nw(block)

      COMMON /SFCDRAG/ DRAG(imr,jmr+1),IWATER(13104)
!hpf$ distribute drag(block,block)
!hpf$ distribute iwater(block)


      integer,dimension(NumMinVols,150)::flag,icount,valid
!hpf$ distribute flag(block,block)
!hpf$ distribute icount(block,block)
!hpf$ distribute valid(block,block)
      integer,dimension(NumMinVols,150)::mmcount1
!hpf$ distribute mmcount1(block,block)

      integer,dimension(NumMinVols):: index
!hpf$ distribute index(block)
      real,dimension(NumMinVols):: levind
!hpf$ distribute levind(block)

      integer,dimension(NumMinVols)::nhmask,kxx,iigpx,iigpy
!hpf$ distribute nhmask(block)
!hpf$ distribute kxx(block)
!hpf$ distribute iigpx(block)
!hpf$ distribute iigpy(block)
      parameter (zero=0.0,one=1.0,two=2.0,onept1=1.10)
      parameter (ninety=90.0,one80=180.0,tenth=0.10,pt05=0.05)
      parameter (tenm4=0.0001)


      real,dimension(13104):: volcentre_x,volcentre_y
      integer,dimension(13104):: nh
      real,dimension(13104):: glatup,glonup
!hpf$ distribute volcentre_x(block)
!hpf$ distribute volcentre_y(block)
!hpf$ distribute glatup(block)
!hpf$ distribute glonup(block)
!hpf$ distribute nh(block) 
 
      character*60 message
      integer,dimension(NumMinVols):: nhgp,nh_count
!hpf$ distribute nhgp(block)
!hpf$ distribute nh_count(block)
      real,dimension(NumMinVols):: glatx,glonx
!hpf$ distribute glatx(block)
!hpf$ distribute glonx(block)
      
      LOGICAL,dimension(NumMinVols):: WATER
!hpf$ distribute water(block)
      COMMON /GRIDXX/   IM, JNP, MLEV

      logical first
      data first/.true./

      message(1:60) = ' '
 

      if(mv_end.ne.(mv_begin+NumMinVols-1)) then
             print *, 'Mini Volume indices messed up!'
             stop
      end if
      CALL GETCONI('MAXREP  ',MAXREP)
      if(MAXREP*3.GT.MAXDAT) then
          message = ' PUV_PREP_COV: Value of MAXREP exceeds MAXDAT.'
          call flash( 5,message,1 )
      end if
 
      call getconh('DLAT    ',dlat)
      call getconh('DLON    ',dlon)

      nhgp=nh(mv_begin:mv_end)
      max_nhgp=maxval(nhgp)
      max_nhgp2=maxval(nhgp/2)
 
!  Locate grid points within a volume.

      nh_count(:) = 0
      glatx=0
      glonx=0
      glon=0
      glat=0
      gstd=0
      gplat=0
      gplon=0

      do j = 1,2
       do i = 1,max_nhgp2

c         print *, 'within the lat-lon loop'
         where(i.le.nhgp/2)
           nhmask=1
         else where
           nhmask=0
         end where

         forall(igp=1:NumMinVols,nhmask(igp).eq.1)
     &     nh_count(igp) = nh_count(igp) + 1

         glatx(:) = glatup(mv_begin:mv_end) + ( j-1 )*dlat
         where ( glatx .gt. 90. ) glatx = 90.
         forall(igp=1:NumMinVols,nhmask(igp).eq.1)
     &     glat( igp, nh_count(igp) ) = glatx(igp)

         glonx(:) = glonup(mv_begin:mv_end ) + ( i-1 )*dlon
         forall(igp=1:NumMinVols,nhmask(igp).eq.1)
     &     glon( igp,nh_count(igp) ) = glonx(igp)

c         print *, 'before iigpx'
         iigpx=(GLATX+NINETY)/DLAT + ONEPT1
         iigpy=(glonx+one80)/dlon + onept1
c         print *, 'after iigpy'

         forall(igp=1:NumMinVols,nhmask(igp).eq.1)
     &     gplat(igp,nh_count(igp)) = iigpx(igp)
         forall(igp=1:NumMinVols,nhmask(igp).eq.1)
     &     gplon(igp,nh_count(igp))=  iigpy(igp)

         forall(igp=1:NumMinVols,nhmask(igp).eq.1)
     &     gstd(igp,nh_count(igp))= fcerr(iigpx(igp),iigpy(igp))
     

       end do   !        <-- End looping over latitudes.
      end do      !        <-- End looping over longitudes.

                                                                                
 
C  MAKE LISTS OF PRESSURE AND WIND DATA
C  ------------------------------------
 
 
c      print *, 'ready to gather'
c      print *, NumMinVols,nobs,mxreps
c      print *, 'DATA(:,ind_type)'
c      print *, DATA(:,ind_type)
c      print *, 'LEVLINDEX'
c      print *, LEVLINDEX
c      print *, 'TYPEMM'
c      print *, TYPEMM
c      print *, 'NREP'
c      print *, NREP
      CALL MMGATHER(NumMinVols,nobs,mxreps,
     &              DATA(:,ind_type),LEVLINDX,TYPEMM,NREP)
      max_nrep=maxval(nrep)
 
      valid=0
      forall(igp=1:NumMinVols,j=1:mxreps,j.le.nrep(igp))
     &   valid(igp,j)=1

      where(valid.eq.1.and.TYPEMM .NE. TWO) 
             flag=1
      end where
c      print *, 'setting flag'
      where(valid.eq.1.and.(.not.(TYPEMM .NE. TWO)))
             flag=0
      end where
      where(valid.ne.1)
          flag=0
      end where

      icount=0
      icount(:,1)=flag(:,1)
      do ipt=2,max_nrep
        icount(:,ipt)=icount(:,ipt-1)+flag(:,ipt)
      end do
      do ipt=1,max_nrep
        icount(:,ipt)=icount(:,ipt)*flag(:,ipt)
      end do

c      print *, 'done with icount'
      plstmm=0.
c      do ipt=1,max_nrep
c       index(:NumMinVols)=icount(:NumMinVols,ipt)
c       levind(:NumMinVols)=int(levlindx(:NumMinVols,ipt))
c       forall(igp=1:NumMinVols,index(igp).ne.0)
c     &      plstmm(igp,index(igp))=levind(igp)
c      end do

c      print *, 'Damn loop!!!!!!'
      do ipt=1,max_nrep
       do igp=1,NumMinVols
        if(icount(igp,ipt).ne.0) then
c            print *, igp,ipt,icount(igp,ipt),levlindx(igp,ipt)
            plstmm(igp,icount(igp,ipt))=levlindx(igp,ipt)
        end if
       enddo
      end do

c      print *, 'done with plstmm'
      np(:)=maxval(icount(:,:),dim=2)
c      print *, 'done with np'

      where(valid.eq.1.and.TYPEMM .NE. ONE)
             flag=1
      else where
             flag=0
      end where
c      print *, 'finished pressure list'

      icount=0
      wlstmm=0.
      
c      forall(ipt=1:max_nrep)
c     &   icount(:,ipt)=sum(flag(:,1:ipt),dim=2)*flag(:,ipt)
c      index=0
c      do ipt=1,max_nrep
c        do igp=1,NumMinVols
c           index(igp)=icount(igp,ipt)
c           if(index(igp).ne.0) then
c              wlst(igp,index(igp))=levlindx(igp,ipt)
c           endif
c        enddo
c      end do

      icount=0
      icount(:,1)=flag(:,1)
      do ipt=2,max_nrep
        icount(:,ipt)=icount(:,ipt-1)+flag(:,ipt)
      end do
      do ipt=1,max_nrep
        icount(:,ipt)=icount(:,ipt)*flag(:,ipt)
      end do

c      forall(igp=1:NumMinVols,ipt=1:max_nrep)
c     & mmcount1(igp,ipt)=sum(flag(igp,1:ipt),dim=2)*flag(igp,ipt)
c      forall(igp=1:NumMinVols,ipt=1:max_nrep,mmcount1(igp,ipt).ne.0)
c     &      wlstmm(igp,mmcount1(igp,ipt))=levlindx(igp,ipt)

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c      Compiler bug here!!!!
c      It cannot compile correctly this piece of code
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c      forall(ipt=1:max_nrep)
c     &   icount(:,ipt)=sum(flag(:,1:ipt),dim=2)*flag(:,ipt)
c      forall(igp=1:NumMinVols,ipt=1:max_nrep,icount(igp,ipt).ne.0)
c     &      wlst(igp,icount(igp,ipt))=levlindx(igp,ipt)
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c      print *, 'Damn loop!!!!!!'
      do ipt=1,max_nrep
       do igp=1,NumMinVols
        if(icount(igp,ipt).ne.0) then
            wlstmm(igp,icount(igp,ipt))=levlindx(igp,ipt)
        end if
       enddo
      end do

 
      nw(:)=maxval(icount(:,:),dim=2)

c      print *, 'finished wind list'
C.......................................................
C
C     IF THE MINI-VOLUME IS NOT A WATER POINT
C     THEN ONLY PERFORM PRESSURE ANALYSIS
C
C.......................................................
      kxx (:)  = (gplat(:,1) - 1)*IM + gplon(:,1)
      WATER(:) = IWATER(KXX).EQ.1
      WHERE(.NOT.WATER(:)) NW(:) = 0
 
 
c      print *, 'finished water'
C  GATHER STUFF FOR PRESSURE-PRESSURE COVARIANCE
C  ---------------------------------------------
c   Note the ordering of 1st data gather then model copy
c   The reason is that target arrays are initialized to 1e+16
c   in mmgather for debugging purposes
cccccccccccccccccccccccccccccccccccccccccccccccccccccc

      CALL MMGATHER(NumMinVols,nobs,700,
     &              DATA(:,ind_lat),PLSTMM,PLAT,NP)
      CALL MMGATHER(NumMinVols,nobs,700,
     &              DATA(:,ind_lon),PLSTMM,PLON,NP)
      CALL MMGATHER(NumMinVols,nobs,700,
     &              DATA(:,ind_fcer),PLSTMM,PSTD,NP)
      CALL MMGATHER(NumMinVols,nobs,700,
     &              DATA(:,ind_slp),PLSTMM,PDEL,NP)
      CALL MMGATHER(NumMinVols,nobs,700,
     &              DATA(:,ind_kx ),PLSTMM,PKX ,NP)

c      print *, 'finished final pressure gathers'
       forall(igp=1:NumMinVols)
     &   PLAT(igp,NP(igp)+1:NP(igp)+max_nhgp) = GLAT(igp,1:max_nhgp)
       forall(igp=1:NumMinVols)
     &   PLON(igp,NP(igp)+1:NP(igp)+max_nhgp) = GLON(igp,1:max_nhgp)
       forall(igp=1:NumMinVols)
     &   PSTD(igp,NP(igp)+1:NP(igp)+max_nhgp) = GSTD(igp,1:max_nhgp)
C  GATHER STUFF FOR COVARIANCES WITH WIND TERMS
C  --------------------------------------------
 
      CALL MMGATHER(NumMinVols,nobs,700,
     &              DATA(:,ind_lat),WLSTMM,WLAT,NW)
      CALL MMGATHER(NumMinVols,nobs,700,
     &              DATA(:,ind_lon),WLSTMM,WLON,NW)
      CALL MMGATHER(NumMinVols,nobs,700,
     &              DATA(:,ind_fcer),WLSTMM,WSTD,NW)
      CALL MMGATHER(NumMinVols,nobs,700,
     &              DATA(:,ind_uwnd),WLSTMM,UDEL,NW)
      CALL MMGATHER(NumMinVols,nobs,700,
     &              DATA(:,ind_vwnd),WLSTMM,VDEL,NW)
      CALL MMGATHER(NumMinVols,nobs,700,
     &              DATA(:,ind_kx ),WLSTMM,WKX ,NW)
 
c      print *, 'finished final wind gathers'
       forall(igp=1:NumMinVols)
     &   WLAT(igp,NW(igp)+1:NW(igp)+max_nhgp) = GLAT(igp,1:max_nhgp)
       forall(igp=1:NumMinVols)
     &   WLON(igp,NW(igp)+1:NW(igp)+max_nhgp) = GLON(igp,1:max_nhgp)
       forall(igp=1:NumMinVols)
     &   WSTD(igp,NW(igp)+1:NW(igp)+max_nhgp) = GSTD(igp,1:max_nhgp)
 
      print *, 'from puv_prep_cov'
      print *, np
      print *, nw
      RETURN
      END
