!
!
!##################################################################
!##################################################################
!######                                                      ######
!######                SUBROUTINE CELTRK                     ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE celtrk(nx,ny,nz,ibeg,iend,jbeg,jend,kbeg,kend,               & 1,5
           w, qr, rhobar, x,y,zp, ntrcell, xcw, ycw,                    &
           ireturn,                                                     &
           ihead,itail,jseg,kcell,icmin,icmax,                          &
           xs,ys,tem1,tem2)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Driver for cell tracking subroutines.
!  Part of the ARPS cell tracking subsystem.
!
!  AUTHOR: Keith Brewster
!   22 April 1993
!
!-----------------------------------------------------------------------
!
!  INPUT VARIABLES
!
!    x        x coordinate of grid points in physical/comp. space (m)
!    y        y coordinate of grid points in physical/comp. space (m)
!    zp       z coordinate of grid points in physical space (m)
!
!    w        vertical component of velocity in Cartesian
!             coordinates (m/s).
!
!    qr       Rain water mixing ratio (kg/kg)
!
!    rhobar   Mean density (kg/m^3)
!
!  OUTPUT:
!
!  Writes tracking information to a file, and returns the following
!  variables:
!
!  ntrcell    Number of cells found at this time
!  xcw        Weighted central x position of all cells at this time
!  ycw        Weighted central y position of all cells at this time
!
!  ireturn    return status
!             =0  center of "mass" successfully determined
!             =-1 no cells found so center of "mass" only estimated based
!                 on previous location
!
!  WORK ARRAYS:
!
!    tem1     Temporary work array.
!    tem2     Temporary work array.
!
!  (These arrays are defined and used locally (i.e. inside this
!   subroutine), they may also be passed into routines called by
!   this one. Exiting the call to this subroutine, these temporary
!   work arrays may be used for other purposes, and therefore their
!   contents may be overwritten. Please examine the usage of work
!   arrays before you make any change to the code.)
!
!-----------------------------------------------------------------------
!
!  Variable Declarations:
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
!
!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------
!
  INCLUDE 'globcst.inc'
  INCLUDE 'grid.inc'          ! Grid parameters
!
!-----------------------------------------------------------------------
!
!  Arguments
!
!-----------------------------------------------------------------------
!
  INTEGER :: nx,ny,nz
  INTEGER :: ibeg,iend,jbeg,jend,kbeg,kend
!
  REAL :: x(nx)
  REAL :: y(ny)
  REAL :: zp(nx,ny,nz)

  REAL :: w     (nx,ny,nz)
  REAL :: qr    (nx,ny,nz)
  REAL :: rhobar(nx,ny,nz)
!
  INTEGER :: ntrcell
  REAL :: xcw,ycw
  INTEGER :: ireturn
!
!-----------------------------------------------------------------------
!
!  Work Arrays
!
!-----------------------------------------------------------------------
!
  REAL :: tem1(nx,ny,nz)
  REAL :: tem2(nx,ny,nz)
  REAL :: xs(nx),ys(ny)
!
!-----------------------------------------------------------------------
!
!  Array size parameters.
!
!  maxcell: maximum number of cells
!           not all may be present at one time
!
!  maxtime: number of cell locations over time to save
!
!-----------------------------------------------------------------------
!
  INTEGER :: maxcell,maxtime
  PARAMETER (maxcell=26, maxtime=15)
!
!-----------------------------------------------------------------------
!
!  Important adjustable parameters.
!
!  maxhgt:  maximum height to look for cells.
!  thresh:  minimum value of vertical velocity used to define a cell
!     eps:  small margin of error to allow to pass as meeting threshold
!           if thresh is not exactly met
!  maxblw:  maximum number of points meeting thresh-eps threshold but
!           not thresh threshold in a segment
!  minlen:  minimum length of a segment (meters)
!    minarea:  minimum area of a horizontal cell (m*m)
!  minvol:  minimum volume of a 3-d cell (m*m*m)
!   minqr:  minimum rainwater maximum within a 3-d cell (kg/kg)
!  disthr:  maximum distance from an extrapolated cell location
!           to match with a current cell location (m)
!
!
!-----------------------------------------------------------------------
!
  INTEGER :: maxblw
  REAL :: thresh,eps,maxhgt,minlen,minqr,minarea,minvol,disthr
  PARAMETER (maxhgt=12000.,    & ! m
         thresh=5.0,       & ! m/s
            eps=0.5,       & ! m/s
         maxblw=2,         & ! grid points
         minlen=2000.,     & ! m
        minarea=9.0E06,    & ! m*m
         minvol=50.0E09,   & ! m*m*m
          minqr=0.001,     & ! kg/kg maximum in cell
         disthr=4000.)     ! m
!
!
!-----------------------------------------------------------------------
!
!  Variables to determine cell location at a given time.
!
!-----------------------------------------------------------------------
!

  INTEGER :: ihead(nx*ny),itail(nx*ny),jseg(nx*ny),                     &
          kcell(nx*ny)

  REAL :: xcent(maxcell,2)
  REAL :: ycent(maxcell,2)
  REAL :: zcent(maxcell,2)
  REAL :: qrmax(maxcell,2)
  REAL :: wgt(maxcell,2)
  REAL :: area(maxcell,2)
  INTEGER :: icmin(maxcell,ny,2)
  INTEGER :: icmax(maxcell,ny,2)
  INTEGER :: ncell
!
!-----------------------------------------------------------------------
!
!  Location-to-motion variables
!
!-----------------------------------------------------------------------
!
  CHARACTER (LEN=1) :: cname(maxcell)
  INTEGER :: ktime,ncalls
  INTEGER :: nctim(maxtime)
  REAL :: alltime(maxtime)
  REAL :: xpos(maxtime,maxcell)
  REAL :: ypos(maxtime,maxcell)
  REAL :: zpos(maxtime,maxcell)
  REAL :: xnot(maxcell)
  REAL :: ynot(maxcell)
  REAL :: allvol(maxtime,maxcell)
  REAL :: allwgt(maxtime,maxcell)
  REAL :: umotion(maxtime,maxcell)
  REAL :: vmotion(maxtime,maxcell)
  REAL :: wmotion(maxtime,maxcell)
  REAL :: growth(maxtime,maxcell)
!
!-----------------------------------------------------------------------
!
!  Save cname and all variables that are a function of time.
!
!-----------------------------------------------------------------------
!
  SAVE cname
  SAVE ktime,ncalls
  SAVE nctim
  SAVE alltime,xpos,ypos,zpos,xnot,ynot
  SAVE allvol,allwgt,umotion,vmotion,wmotion,growth
!
  DATA cname /'A','B','C','D','E','F','G','H','I','J',                  &
              'K','L','M','N','O','P','Q','R','S','T',                  &
              'U','V','W','X','Y','Z'/
  DATA ktime,ncalls /0,0/
!
!
!-----------------------------------------------------------------------
!
!  Misc. local variables
!
!-----------------------------------------------------------------------
!
  REAL :: prtmis
  PARAMETER (prtmis= 999.)
  CHARACTER (LEN=1) :: comma
  PARAMETER (comma=',')
  REAL :: wgtsum,xsum,ysum
  INTEGER :: maxseg
  INTEGER :: istr,ifin,jstr,jfin,kstr,kfin
  INTEGER :: i,j,kc,mtime,mpast
!
  CHARACTER (LEN=80  ) :: trkfn  ! Name of max./min. file
  INTEGER :: ltrkfn
  INTEGER :: istat, nchtrk
  SAVE nchtrk
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  maxseg = nx*ny

  IF( ((nx*nz)/2) < maxcell) THEN
    WRITE(6,'(1x,a,/1x,a,I4,a,/1x,a)')                                  &
    'The size of temporary arrays passed to CELTRK not big enough& !',  &
    '(nx*nz)/2 must be greater than ',maxcell,                          &
    ' for CELTRK to work properly.',                                    &
    ' No cell tracking was done.'
    RETURN
  END IF

  IF(ncalls == 0) THEN

    ncalls=1

    trkfn = runname(1:lfnkey)//'.track'
    ltrkfn = lfnkey + 6

    CALL getunit(nchtrk)

    WRITE(6,'(1x,a,a,a/,1x,a)')                                         &
        'Check to see if file ',trkfn(1:ltrkfn),' already exists.',     &
        'If so, append a version number to the filename.'

    CALL fnversn( trkfn, ltrkfn )

    OPEN(nchtrk,FORM='formatted',STATUS='new',                          &
                FILE=trkfn(1:ltrkfn),IOSTAT=istat)

    IF( istat /= 0) THEN

      WRITE(6,'(/a,i2,/a/)')                                            &
          ' Error occured when opening file '//trkfn(1:ltrkfn)//        &
          ' using FORTRAN unit ',nchtrk,' Program stopped in trkdriv.'
      CALL arpsstop ("arpstop called from celtrk",1)

    END IF

    WRITE(nchtrk,'(a)') ''''//runname//''''

    WRITE(nchtrk,'(/a,a)')                                              &
        '   time  name  xloc   yloc   zloc',                            &
        '    u      v      w     vol       wgt  '
  END IF
!
!-----------------------------------------------------------------------
!
!    Check the limits that were passed in.
!    For bookkeeping reasons (avoid array out of bounds and some
!    storage tricks), as well as physical sense, the range of i,j,k
!    should not involve the artificial boundary points.
!
!    Use istr,ifin,jstr,jfin,kstr,kfin in place of ibeg,iend, etc.
!    for the call to loccell.
!
!-----------------------------------------------------------------------
!
  istr=MAX(ibeg,2)
  ifin=MIN(iend,nx-1)
  jstr=MAX(jbeg,2)
  jfin=MIN(jend,ny-1)
  kstr=MAX(kbeg,2)
  kfin=MIN(kend,nz-1)
!
!-----------------------------------------------------------------------
!
!  Locate cells in grid.
!
!-----------------------------------------------------------------------
!
  DO i=1,nx-1
    xs(i)=(x(i)+x(i+1))*0.5
  END DO
  xs(nx)=2.*xs(nx-1)-xs(nx-2)
  DO j=1,ny-1
    ys(j)=(y(j)+y(j+1))*0.5
  END DO
  ys(ny)=2.*ys(ny-1)-ys(ny-2)

  CALL loccell(nx,ny,nz, maxseg,maxcell,ncell,                          &
              istr,ifin,jstr,jfin,kstr,kfin,                            &
              w,qr, rhobar, xs,ys,zp,                                   &
              ihead,itail,kcell,jseg,                                   &
              xcent,ycent,zcent,icmin,icmax,                            &
              qrmax,area,wgt,maxhgt,                                    &
              thresh,eps,maxblw,minlen,minqr,minarea,minvol,            &
              ireturn, tem1,tem2)
!
!-----------------------------------------------------------------------
!
!  Add movement of grid to the cell locations.
!  Compute a weighted mean location of all cells at this time.
!
!-----------------------------------------------------------------------
!
  xsum=0.
  ysum=0.
  wgtsum=0.
  DO kc=1,ncell
    xcent(kc,1)=xcent(kc,1)+xgrdorg
    ycent(kc,1)=ycent(kc,1)+ygrdorg
    wgtsum=wgtsum+wgt(kc,1)
    xsum=xsum+xcent(kc,1)*wgt(kc,1)
    ysum=ysum+ycent(kc,1)*wgt(kc,1)
  END DO
!
  ntrcell=ncell
!
  IF(wgtsum > 0.) THEN
    ireturn=0
    xcw=xsum/wgtsum
    ycw=ysum/wgtsum
  ELSE
    ireturn=-1
    IF( ktime > 1) THEN
      xcw=xcw+umove*tceltrk
      ycw=ycw+vmove*tceltrk
    ELSE
      xcw=(xs(1)+xs(nx))*0.5
      ycw=(ys(1)+ys(ny))*0.5
    END IF
  END IF
!
!-----------------------------------------------------------------------
!
!  Spit out some diagnostics to the list file.
!
!-----------------------------------------------------------------------
!
  WRITE(6,'(//a,i4,/a,a)')                                              &
      ' Cell positions for whole vol  Total cells: ',ncell,             &
      ' kc   xcabs   ycabs    xcrel    ycrel  ',                        &
      '  zcent     vol       wgt     qrmax'
  DO kc=1,ncell
    WRITE(6,'(1x,i4,2(F9.2,F9.2),F9.2,F9.0,F12.0)')                     &
           kc,(xcent(kc,1)*0.001),(ycent(kc,1)*0.001)                   &
             ,((xcent(kc,1)-xgrdorg)*0.001)                             &
             ,((ycent(kc,1)-ygrdorg)*0.001)                             &
             ,(zcent(kc,1)*0.001)                                       &
             ,( area(kc,1)*1.0E-09),wgt(kc,1)                           &
             ,(qrmax(kc,1)*1000.)
  END DO
!
!-----------------------------------------------------------------------
!
!  Match current cell volumes with previous.
!
!-----------------------------------------------------------------------
!
!  print *, ' Calling matctim'
  CALL matctim(maxcell,maxtime,ktime,mtime,mpast,                       &
              disthr,curtim,umove,vmove,                                &
              ncell,xcent(1,1),ycent(1,1),zcent(1,1),                   &
              area(1,1),wgt(1,1),                                       &
              cname,nctim,alltime,                                      &
              xpos,ypos,zpos,xnot,ynot,                                 &
              allvol,allwgt,                                            &
              umotion,vmotion,wmotion)
  WRITE(6,'(//a,i4,/a)')                                                &
      ' Cells combined in time arrays: ',nctim(mtime),                  &
      '   jc   xpos    ypos   zpos   volume     wgt'
  DO kc=1,nctim(mtime)
    IF(xpos(mtime,kc) > -1.e30) THEN
      WRITE(6,'(3x,a1,F9.2,F9.2,F9.2,F9.0,F12.0)')                      &
           cname(kc),(xpos(mtime,kc)*0.001),                            &
              (ypos(mtime,kc)*0.001),                                   &
              (zpos(mtime,kc)*0.001),                                   &
            (allvol(mtime,kc)*1.0E-09),                                 &
             allwgt(mtime,kc)
    END IF
  END DO
!
!-----------------------------------------------------------------------
!
!  Calculate cell motion vector and growth rate.
!
!-----------------------------------------------------------------------
!
  CALL getvec(maxcell,maxtime,ktime,mtime,                              &
              cname,nctim,alltime,                                      &
              xpos,ypos,zpos,xnot,ynot,                                 &
              allvol,allwgt,                                            &
              umotion,vmotion,wmotion,growth)
!
!-----------------------------------------------------------------------
!
!  Spit out motion vector results.
!
!-----------------------------------------------------------------------
!
  IF(nctim(mtime) > 0) THEN
    DO kc=1,nctim(mtime)
      IF(xpos(mtime,kc) > -1.0E30) THEN
        IF(umotion(mtime,kc) > -1.0E30) THEN
          WRITE(nchtrk,'(1x,F7.0,2x,a1,1x,f7.2,                         &
  &           f7.2,f7.2,f7.2,f7.2,f7.2,f7.0,f10.0)')                    &
              alltime(mtime),cname(kc),                                 &
              (xpos(mtime,kc)*0.001),                                   &
              (ypos(mtime,kc)*0.001),                                   &
              (zpos(mtime,kc)*0.001),                                   &
              umotion(mtime,kc),                                        &
              vmotion(mtime,kc),                                        &
              wmotion(mtime,kc),                                        &
              (allvol(mtime,kc)*1.0E-09),                               &
              allwgt(mtime,kc)
        ELSE
          WRITE(nchtrk,'(1x,F7.0,2x,a1,1x,                              &
  &           f7.2,f7.2,f7.2,f7.2,f7.2,f7.2,f7.0,f10.0)')               &
              alltime(mtime),cname(kc),                                 &
              (xpos(mtime,kc)*0.001),                                   &
              (ypos(mtime,kc)*0.001),                                   &
              (zpos(mtime,kc)*0.001),                                   &
                  prtmis,                                               &
                  prtmis,                                               &
                  prtmis,                                               &
              (allvol(mtime,kc)*1.0E-09),                               &
              allwgt(mtime,kc)
        END IF
      END IF
    END DO
  END IF
  WRITE(6,'(/a,i8,f9.0,f9.0//)') ' End of tracking for time:',          &
             mtime,alltime(mtime),curtim
  RETURN
END SUBROUTINE celtrk
!
!##################################################################
!##################################################################
!######                                                      ######
!######                SUBROUTINE LOCCELL                    ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE loccell(nx,ny,nz, maxseg,maxcell,ncell,                      & 1,11
           ibeg,iend,jbeg,jend,kbeg,kend,                               &
           tvar,qr,rhobar,x,y,zp,                                       &
           ihead,itail,kcell,jseg,                                      &
           xcent,ycent,zcent,icmin,icmax,                               &
           qrmax,area,wgt,maxhgt,                                       &
           thresh,eps,maxblw,minlen,minqr,minarea,minvol,               &
           ireturn, tem1,tem2)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Builds cells from tracking variable (tvar, i.e. w) field and
!  establishes a central postion for each which is used for
!  storm tracking.
!
!  A module of the ARPS cell tracking subsystem.
!
!  AUTHOR: Keith Brewster
!   21 Mar 1993
!
!-----------------------------------------------------------------------
!
!  DATA ARRAYS:
!
!    x        x coordinate of grid points in physical/comp. space (m)
!    y        y coordinate of grid points in physical/comp. space (m)
!    zp       z coordinate of grid points in physical space (m)
!
!    w        z component of velocity (m/s)
!
!    qr       Rainwater mixing ratio (kg/kg)
!
!    rhobar   Mean density (kg/m^3)
!
!  WORK ARRAYS:
!
!    tem1     Temporary work array.
!    tem2     Temporary work array.
!
!  (These arrays are defined and used locally (i.e. inside this
!   subroutine), they may also be passed into routines called by
!   this one. Exiting the call to this subroutine, these temporary
!   work arrays may be used for other purposes, and therefore their
!   contents may be overwritten. Please examine the usage of work
!   arrays before you make any change to the code.)
!
!-----------------------------------------------------------------------
!
!  Variable Declarations:
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
  INTEGER :: nx,ny,nz
  INTEGER :: maxseg
  INTEGER :: maxcell
  INTEGER :: ncell
  INTEGER :: ibeg,iend,jbeg,jend,kbeg,kend
  REAL :: x(nx)
  REAL :: y(ny)
  REAL :: zp(nx,ny,nz)
  REAL :: tvar(nx,ny,nz)
  REAL :: qr(nx,ny,nz)
  REAL :: rhobar(nx,ny,nz)
  REAL :: thresh
  INTEGER :: ihead(maxseg)
  INTEGER :: itail(maxseg)
  INTEGER :: jseg(maxseg)
  INTEGER :: kcell(maxseg)
  REAL :: xcent(maxcell,2)
  REAL :: ycent(maxcell,2)
  REAL :: zcent(maxcell,2)
  REAL :: wgt(maxcell,2)
  REAL :: qrmax(maxcell,2)
  REAL :: area(maxcell,2)
  INTEGER :: icmin(maxcell,ny,2)
  INTEGER :: icmax(maxcell,ny,2)
  INTEGER :: maxblw
  REAL :: eps,maxhgt,minlen,minqr,minarea,minvol
  REAL :: tem1(nx,ny,nz)
  REAL :: tem2(nx,ny,nz)
  INTEGER :: ireturn
!
!-----------------------------------------------------------------------
!
!  Misc. local variables
!
!-----------------------------------------------------------------------
!
  INTEGER :: nseg,onvf
  INTEGER :: i,j,k,kc,ncellk
  REAL :: tmin,tmax
  INTEGER :: ilmax,jlmax,klmax,ilmin,jlmin,klmin
  INTEGER :: icengd,jcengd,iseg,ihd,itl,kntblw
  REAL :: xlen,thmin,zsfc
  LOGICAL :: onseg
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!
!-----------------------------------------------------------------------
!
!  Average qr in the vertical to get qr at w points.
!  This is stored in tem1.
!
!-----------------------------------------------------------------------
!
  onvf=1
  CALL avgz(qr, onvf, nx,ny,nz,                                         &
             1,nx-1,1,ny-1,2,nz-1, tem1)

!
!-----------------------------------------------------------------------
!
!  Average rhobar in the vertical to get rhobar at w points.
!  This is stored in tem2.
!
!-----------------------------------------------------------------------
!
  onvf=1
  CALL avgz(rhobar, onvf, nx,ny,nz,                                     &
             1,nx-1,1,ny-1,2,nz-1, tem2)
!
!-----------------------------------------------------------------------
!
!  Find the w max and report it.
!  For now this information is only used for diagnostics,
!  but in the future the cell matching could begin at the level
!  of max w and work up and down from there.
!  For now, it begins at the top and works down.
!
!-----------------------------------------------------------------------
!
  CALL a3dmax(tvar, 1,nx,ibeg,iend, 1,ny,jbeg,jend,                     &
                    1,nz,kbeg,kend, tmax,tmin,                          &
                    ilmax,jlmax,klmax,ilmin,jlmin,klmin)
  PRINT *, ' w max = ',tmax
!  print *, ' imax,jmax,kmax = ',ilmax,jlmax,klmax
!  print *, ' thresh= ',thresh
!  print *, ' minlen= ',minlen
!
!-----------------------------------------------------------------------
!
!  Search for segments, strings of points meeting the threshold.
!  The variable onseg is a flag indicating whether we are on a
!  segment looking for its end (onseg=true) or looking for the
!  beginning of the next segment (onseg=false).
!
!  When the end of a segment is reached matchsg is called which
!  matches the current segment to any previously found segments
!  on this level.
!
!  Variables ihead and itail describe the beginning and end of
!  the each segment.
!
!  At the end of each level, matchsg is called in reverse
!  order to see if any further matching of cells can be made
!  (since the matching process can be directional dependent).
!  In the following example segments A and B would not be matched
!  to each other going from 1-to-3 but on the reverse matching
!  (3-to-1) B would be matched to C which had originally been
!  matched to A.
!
!  3      DDDD     EEEE
!  2        CCCCCCCCCCCCC
!  1         AAA      BBBB
!
!-----------------------------------------------------------------------
!
  icengd=((nx-3)/2) + 1
  jcengd=((ny-3)/2) + 1
  zsfc=0.5*(zp(icengd,jcengd,1)+zp(icengd,jcengd,2))
  thmin=thresh-eps
  ncell=0
  DO k=kend,kbeg,-1
!    print *, ' k = ',k
    IF((zp(icengd,jcengd,k)-zsfc) <= maxhgt) THEN
      nseg=0
      ncellk=0
      IF(tmax > thresh) THEN
        DO j=jbeg,jend
          onseg=.false.
          DO i=ibeg,iend
            IF(onseg) THEN
              IF(tvar(i,j,k) < thresh) THEN
                kntblw=kntblw+1
                IF(tvar(i,j,k) >= thmin) THEN
                  IF(kntblw > maxblw) THEN
                    itl=(i-kntblw)
                    xlen=ABS(                    IF(xlen >= minlen) THEN
                      IF(nseg < maxseg) THEN
                        nseg=nseg+1
                        ihead(nseg)=ihd
                        itail(nseg)=itl
                        jseg(nseg)=j
                        kcell(nseg)=-99
                        CALL matchsg(maxseg,maxcell,                    &
                            nseg, 1,(nseg-1),1,ihd,itl,j,               &
                            ihead,itail,jseg,kcell,ncellk)
                      END IF
                    END IF
                    onseg=.false.
                  END IF
                ELSE  ! LT.thresh and LT.thmin
                  itl=(i-kntblw)
                  xlen=ABS(                  IF(xlen >= minlen) THEN
                    IF(nseg < maxseg) THEN
                      nseg=nseg+1
                      ihead(nseg)=ihd
                      itail(nseg)=itl
                      jseg(nseg)=j
                      kcell(nseg)=-99
                      CALL matchsg(maxseg,maxcell,                      &
                          nseg, 1,(nseg-1),1,ihd,itl,j,                 &
                          ihead,itail,jseg,kcell,ncellk)
                    END IF
                  END IF
                  onseg=.false.
                END IF
!
!-----------------------------------------------------------------------
!
!  Onseg and tvar is greater than thresh
!  Reset kntblw to zero,
!  then resume search for the end of the segment.
!
!-----------------------------------------------------------------------
!
              ELSE
                kntblw=0
              END IF
!
!-----------------------------------------------------------------------
!
!  Not "onseg"
!  See if the threshold is exceeded, if so, initialize
!  segment pointers, ihd and kntblw, and set onseg to true.
!
!-----------------------------------------------------------------------
!
            ELSE
              IF(tvar(i,j,k) >= thresh) THEN
                ihd=i
                kntblw=0
                onseg=.true.
              END IF
            END IF
          END DO
!
!-----------------------------------------------------------------------
!
!  If you've reached the end of the i indices and your are
!  still "onseg", then close this segment.
!
!-----------------------------------------------------------------------
!
          IF(onseg) THEN
            itl=iend-kntblw
            xlen=ABS(            IF(xlen >= minlen) THEN
              IF(nseg < maxseg) THEN
                nseg=nseg+1
                ihead(nseg)=ihd
                itail(nseg)=itl
                jseg(nseg)=j
                kcell(nseg)=-99
                CALL matchsg(maxseg,maxcell,                            &
                    nseg, 1,(nseg-1),1,ihd,itl,j,                       &
                    ihead,itail,jseg,kcell,ncellk)
              END IF
            END IF
            onseg=.false.
          END IF
        END DO
      END IF     ! tmax .gt. thresh
!
!-----------------------------------------------------------------------
!
!  Do backwards matching of segments
!
!-----------------------------------------------------------------------
!
!    print *, '  before backward matching ncellk = ',ncellk
      DO iseg=(nseg-1),1,-1
        CALL matchsg(maxseg,maxcell,iseg, nseg,(iseg+1),-1,             &
                     ihead(iseg),itail(iseg),jseg(iseg),                &
                     ihead,itail,jseg,kcell,ncellk)
      END DO
!    print *, '  after backward matching ncell = ',ncellk
!
!-----------------------------------------------------------------------
!
!  Calculate cell areas and centroids.
!
!-----------------------------------------------------------------------
!
      IF(ncellk > 0) THEN
        CALL ctrwgt(nx,ny,maxseg,maxcell,nseg,ncellk,                   &
                  tvar(1,1,k),tem1(1,1,k),tem2(1,1,k),                  &
                  x,y,zp(1,1,k),                                        &
                  ihead,itail,jseg,kcell,                               &
                  xcent(1,2),ycent(1,2),zcent(1,2),                     &
                  qrmax(1,2),area(1,2),wgt(1,2))
!
!-----------------------------------------------------------------------
!
!  Create array describing min and max of i for each j column
!  in each cell.
!
!  This is required for efficient vertical matching of horizontal
!  cells to create 3d cells.
!
!-----------------------------------------------------------------------
!
        DO j=1,ny
          DO kc=1,maxcell
            icmin(kc,j,2)=nx+1
            icmax(kc,j,2)=-99
          END DO
        END DO
        DO iseg=1,nseg
          kc=kcell(iseg)
          icmin(kc,jseg(iseg),2)=                                       &
                         MIN(icmin(kc,jseg(iseg),2),ihead(iseg))
          icmax(kc,jseg(iseg),2)=                                       &
                         MAX(icmax(kc,jseg(iseg),2),itail(iseg))
        END DO
!
!-----------------------------------------------------------------------
!
!  Clean up cell list, eliminate those with small area.
!  Or those not meeting the minimum rainwater (qr) requirement.
!
!-----------------------------------------------------------------------
!
!      print *, ' inside loccell minqr,minarea= ',minqr,minarea
!      print *, ' Before sieve2d, ncell= ',ncellk
        CALL sieve2d(maxcell,ny,ncellk,                                 &
                    xcent(1,2),ycent(1,2),zcent(1,2),                   &
                    qrmax(1,2),area(1,2),wgt(1,2),                      &
                    icmin(1,1,2),icmax(1,1,2),                          &
                    minarea)
!      print *, ' After sieving, ncell= ',ncellk
!
!  tell us about icmin and icmax
!      do 256 kc=1,ncellk
!       print *, ' k, kc = ',k,kc
!      do 255 j=1,ny
!        print *,' j= ',j,' imin= ',icmin(kc,j,2),
!    :                       ' imax= ',icmax(kc,j,2)
!255       continue
!256       continue
      END IF
!
!    IF( ncellk.gt.0 ) THEN
!      write(6,'(//a,i4,a,i4,/a)')
!    : '  Cells on level: ',k,'  Total cells: ',ncellk,
!    : '   kc       xcent       ycent      zcent      area        wgt'
!      DO 300 kc=1,ncellk
!        write(6,'(1x,i4,F9.2,F9.2,F9.2,F9.0,F12.0)')
!    :      kc,xcent(kc,2),ycent(kc,2),zcent(kc,2),area(kc,2),wgt(kc,2)
!00       CONTINUE
!    END IF
!
!-----------------------------------------------------------------------
!
!  Combine cells  with cells found at previous levels.
!  The combined cell info is put in the level one index.
!
!-----------------------------------------------------------------------
!
      CALL linkcell(maxcell,nx,ny,nz,k,ncell,ncellk,                    &
                  x,y,zp,xcent,ycent,zcent,icmin,icmax,                 &
                  qrmax,area,wgt)
    END IF    ! below maxhgt?
  END DO
!
!-----------------------------------------------------------------------
!
!  Compute the central x,y,z
!  for each cell volume from the info collected on plane 1.
!
!-----------------------------------------------------------------------
!
  DO kc=1,ncell
    IF(wgt(kc,1) /= 0.0) THEN
      xcent(kc,1)=xcent(kc,1)/wgt(kc,1)
      ycent(kc,1)=ycent(kc,1)/wgt(kc,1)
      zcent(kc,1)=zcent(kc,1)/wgt(kc,1)
    END IF
  END DO
!
!-----------------------------------------------------------------------
!
!  Eliminate cells which are too small.
!
!-----------------------------------------------------------------------
  PRINT *, ' before sieve3d, ncell= ',ncell
  CALL sieve3d(maxcell,ncell,                                           &
               xcent(1,1),ycent(1,1),zcent(1,1),                        &
               qrmax(1,1),area(1,1),wgt(1,1),                           &
               minqr,minvol)
  PRINT *, ' after sieve3d, ncell= ',ncell
!
  RETURN
END SUBROUTINE loccell
!
!
!##################################################################
!##################################################################
!######                                                      ######
!######                SUBROUTINE MATCHSG                    ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE matchsg(maxseg,maxcell,                                      & 4
           kseg, mbeg,mend,mdir,ihd,itl,j,                              &
           ihead,itail,jseg,kcell,ncell)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Matches line segments to form horizontal cells.
!  A module of the ARPS cell tracking subsystem.
!
!  AUTHOR: Keith Brewster
!   21 Mar 1993
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
!  Variable Declarations:
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
  INTEGER :: maxseg,maxcell,kseg
  INTEGER :: ihd,itl,j
  INTEGER :: mbeg,mend,mdir
  INTEGER :: ihead(maxseg),itail(maxseg),jseg(maxseg),                  &
          kcell(maxseg)
  INTEGER :: ncell
!
!-----------------------------------------------------------------------
!
!  Misc. local variables
!
!-----------------------------------------------------------------------
!
  INTEGER :: mseg
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!
!-----------------------------------------------------------------------
!
!  Form cells from line segments
!
!-----------------------------------------------------------------------
!
  DO mseg=mbeg,mend,mdir
    IF(ABS(jseg(mseg)-j) == 1) THEN
      IF(ihd >= ihead(mseg) .AND. ihd <= itail(mseg)) GO TO 120
      IF(itl >= ihead(mseg) .AND. itl <= itail(mseg)) GO TO 120
      IF(itl >= itail(mseg) .AND. ihd <= ihead(mseg)) GO TO 120
    END IF
  END DO
!
!-----------------------------------------------------------------------
!
!  No cell found so this may be a new one.
!
!-----------------------------------------------------------------------
!
  IF(kcell(kseg) < 0) THEN
    IF(ncell < maxcell) THEN
      ncell=ncell+1
      kcell(kseg)=ncell
!      print *, ' kseg j,ihd,itl: ',kseg,j,ihd,itl
!      print *, ' new cell: ',ncell,kcell(kseg)
    ELSE
      PRINT *, ' WARNING: ran out of space for cells.'
    END IF
  END IF
  RETURN
!
!-----------------------------------------------------------------------
!
!  Cell found, it belongs with kcell(mseg)
!
!-----------------------------------------------------------------------
!
  120   CONTINUE
  IF(kcell(mseg) > 0) THEN
    kcell(kseg)=kcell(mseg)
  ELSE
    IF(ncell < maxcell) THEN
      ncell=ncell+1
      kcell(kseg)=ncell
      kcell(mseg)=ncell
    END IF
  END IF
!  print *, ' kseg j,ihd,itl: ',kseg,j,ihd,itl
!  print *, ' matched with cell: ',kcell(mseg),j,
!    :             ihead(kcell(mseg)),itail(kcell(mseg))
  RETURN
END SUBROUTINE matchsg
!
!
!##################################################################
!##################################################################
!######                                                      ######
!######                SUBROUTINE CTRWGT                     ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE ctrwgt(nx,ny,maxseg,maxcell,nseg,ncell,                      & 1
           tvar,qratw,rhoatw,x,y,zp,                                    &
           ihead,itail,jseg,kcell,                                      &
           xcent,ycent,zcent,qrmax,area,wgt)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Compute a central location for all cells formed by matchsg.
!  A module of the ARPS cell tracking subsystem.
!
!  AUTHOR: Keith Brewster
!   21 Mar 1993
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
!  Variable Declarations:
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
  INTEGER :: nx,ny
  INTEGER :: maxseg
  INTEGER :: maxcell
  REAL :: tvar(nx,ny)
  REAL :: qratw(nx,ny)
  REAL :: rhoatw(nx,ny)
  REAL :: x(nx)
  REAL :: y(ny)
  REAL :: zp(nx,ny)
  INTEGER :: ihead(maxseg)
  INTEGER :: itail(maxseg)
  INTEGER :: jseg(maxseg)
  INTEGER :: kcell(maxseg)
  REAL :: xcent(maxcell)
  REAL :: ycent(maxcell)
  REAL :: zcent(maxcell)
  REAL :: qrmax(maxcell)
  REAL :: area(maxcell)
  REAL :: wgt(maxcell)
  INTEGER :: nseg
  INTEGER :: ncell
!
!-----------------------------------------------------------------------
!
!  Misc. local variables
!
!-----------------------------------------------------------------------
!
  REAL :: ddx,ddy,ddxy,wayt
  INTEGER :: i,j,m,kc
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  ddx=x(2)-x(1)
  ddy=y(2)-y(1)
  ddxy=ABS(ddx*ddy)
!
!-----------------------------------------------------------------------
!
!  Initializations.
!
!-----------------------------------------------------------------------
!
  DO kc=1,ncell
    wgt(kc)=0.
    qrmax(kc)=-999.
    area(kc)=0.
    xcent(kc)=0.
    ycent(kc)=0.
    zcent(kc)=0.
  END DO
!
!-----------------------------------------------------------------------
!
!  Go through all segments.  Build sums for each cell.
!
!-----------------------------------------------------------------------
!
  DO m=1,nseg
    IF(kcell(m) > 0) THEN
      j=jseg(m)
      kc=kcell(m)
      DO i=ihead(m),itail(m)
        wayt=rhoatw(i,j)*tvar(i,j)*tvar(i,j)
!         wayt=tvar(i,j)*tvar(i,j)
        wgt(kc)=wgt(kc)+wayt
        xcent(kc)=xcent(kc)+wayt*x(i)
        ycent(kc)=ycent(kc)+wayt*y(j)
        zcent(kc)=zcent(kc)+wayt*zp(i,j)
        qrmax(kc)=AMAX1(qrmax(kc),qratw(i,j))
        area(kc)=area(kc)+ddxy
      END DO
    ELSE
      PRINT *, ' Warning: segment ',m,' not assigned to a cell'
    END IF
  END DO
!
!-----------------------------------------------------------------------
!
!  Calculated cell centroid, a weight average of the locations
!  of all the grid points in the cell.
!
!-----------------------------------------------------------------------
!
  DO kc=1,ncell
    IF(wgt(kc) > 0.) THEN
      xcent(kc)=xcent(kc)/wgt(kc)
      ycent(kc)=ycent(kc)/wgt(kc)
      zcent(kc)=zcent(kc)/wgt(kc)
    END IF
  END DO
!
  RETURN
END SUBROUTINE ctrwgt
!
!##################################################################
!##################################################################
!######                                                      ######
!######               SUBROUTINE SIEVE2D                     ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE sieve2d(maxcell,ny,ncell,                                    & 1
           xcent,ycent,zcent,qrmax,area,wgt,                            &
           icmin,icmax,                                                 &
           minarea)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Discard cells less than a given areal size.
!  A module of the ARPS cell tracking subsystem.
!
!  AUTHOR: Keith Brewster
!   21 Mar 1993
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
!  Variable Declarations:
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
  INTEGER :: maxcell,ny,ncell
  REAL :: xcent(maxcell)
  REAL :: ycent(maxcell)
  REAL :: zcent(maxcell)
  REAL :: qrmax(maxcell)
  REAL :: area(maxcell)
  REAL :: wgt(maxcell)
  INTEGER :: icmin(maxcell,ny)
  INTEGER :: icmax(maxcell,ny)
  REAL :: minarea
!
!-----------------------------------------------------------------------
!
!  Misc. local variables
!
!-----------------------------------------------------------------------
!
  INTEGER :: ic,j,kc,nn
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  nn=ncell
  DO ic=1,ncell
    20      CONTINUE
    IF(ic > nn) EXIT
!    print *, ' ic area, qrmax= ',ic,area(ic),qrmax(ic)
    IF(area(ic) < minarea) THEN
      nn=nn-1
      DO kc=ic,nn
        wgt(kc)=wgt(kc+1)
        area(kc)=area(kc+1)
        qrmax(kc)=qrmax(kc+1)
        xcent(kc)=xcent(kc+1)
        ycent(kc)=ycent(kc+1)
        zcent(kc)=zcent(kc+1)
      END DO
      DO j=1,ny
        DO kc=ic,nn
          icmin(kc,j)=icmin(kc+1,j)
          icmax(kc,j)=icmax(kc+1,j)
        END DO
      END DO
      GO TO 20
    END IF
  END DO
  ncell=nn
  RETURN
END SUBROUTINE sieve2d
!
!
!##################################################################
!##################################################################
!######                                                      ######
!######               SUBROUTINE SIEVE3D                     ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE sieve3d(maxcell,ncell,                                       & 1
           xcent,ycent,zcent,qrmax,area,wgt,                            &
           minqr,minvol)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Discard 3-dimensional cells less than a given volume size.
!  A module of the ARPS cell tracking subsystem.
!
!  AUTHOR: Keith Brewster
!   21 Mar 1993
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
!  Variable Declarations:
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
  INTEGER :: maxcell
  INTEGER :: ncell
  REAL :: xcent(maxcell)
  REAL :: ycent(maxcell)
  REAL :: zcent(maxcell)
  REAL :: qrmax(maxcell)
  REAL :: area(maxcell)
  REAL :: wgt(maxcell)
  REAL :: minqr,minvol
!
!-----------------------------------------------------------------------
!
!  Misc. local variables
!
!-----------------------------------------------------------------------
!
  INTEGER :: ic,kc,nn
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  nn=ncell
  DO ic=1,ncell
    20      CONTINUE
    IF(ic > nn) EXIT
    IF(area(ic) < minvol .OR. qrmax(ic) < minqr) THEN
      nn=nn-1
      DO kc=ic,nn
        wgt(kc)=wgt(kc+1)
        area(kc)=area(kc+1)
        qrmax(kc)=qrmax(kc+1)
        xcent(kc)=xcent(kc+1)
        ycent(kc)=ycent(kc+1)
        zcent(kc)=zcent(kc+1)
      END DO
      GO TO 20
    END IF
  END DO
  ncell=nn
  RETURN
END SUBROUTINE sieve3d
!
!##################################################################
!##################################################################
!######                                                      ######
!######                SUBROUTINE LINKCELL                   ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE linkcell(maxcell,nx,ny,nz,k,ncell,ncellk,                    & 1
           x,y,zp,xcent,ycent,zcent,icmin,icmax,                        &
           qrmax,area,wgt)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Link cells on adjacent horizontal surfaces to form 3-d cell.
!  Combine location info to get 3D location info.
!  A module of the ARPS cell tracking subsystem.
!
!  AUTHOR: Keith Brewster
!   21 Mar 1993
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
!  Variable Declarations:
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
  INTEGER :: maxcell,nx,ny,nz
  INTEGER :: k,ncell,ncellk
  REAL :: x(nx)
  REAL :: y(ny)
  REAL :: zp(nx,ny,nz)
  REAL :: xcent(maxcell,2)
  REAL :: ycent(maxcell,2)
  REAL :: zcent(maxcell,2)
  REAL :: wgt(maxcell,2)
  REAL :: qrmax(maxcell,2)
  REAL :: area(maxcell,2)
  INTEGER :: icmin(maxcell,ny,2)
  INTEGER :: icmax(maxcell,ny,2)
!
!-----------------------------------------------------------------------
!
!  Misc. local variables
!
!-----------------------------------------------------------------------
!
  INTEGER :: ic,jc,j,iloc,jloc
  REAL :: dx,dy
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  NOTE:
!  3-d info is in 3rd index=1,
!  2-d info for a single level is in 3rd index=2.
!
!  Find a cell that has an overlap with this cell.
!  For this nested loop, 2 is the plane being checked
!                     1 is the plane being searched for possible matches
!                    ic is the cell being checked
!                    jc is the cell being compared on plane 1 (3-d info)
!
!-----------------------------------------------------------------------
!
  dx=x(2)-x(1)
  dy=y(2)-y(1)
  DO ic=1,ncellk
    DO jc=1,ncell
!
!-----------------------------------------------------------------------
!
!  Check for overlap between these two cells
!
!-----------------------------------------------------------------------
!
      DO j=1,ny
        IF(icmax(ic,j,2) > 0 .AND. icmax(jc,j,1) > 0) THEN
          IF((icmax(ic,j,2) >= icmin(jc,j,1)) .AND.                     &
              (icmin(ic,j,2) <= icmax(jc,j,1)))                         &
              GO TO 250
        END IF
      END DO
    END DO
!
!-----------------------------------------------------------------------
!
!  No match found, must be a new cell
!
!-----------------------------------------------------------------------
!
    ncell=ncell+1
    jc=ncell
!
    wgt(jc,1)=wgt(ic,2)
    xcent(jc,1)=wgt(ic,2)*xcent(ic,2)
    ycent(jc,1)=wgt(ic,2)*ycent(ic,2)
    zcent(jc,1)=wgt(ic,2)*zcent(ic,2)
    qrmax(jc,1)=qrmax(ic,2)
    iloc=nint((xcent(ic,2)/dx) + 1.5)
    jloc=nint((ycent(ic,2)/dy) + 1.5)
    area(jc,1)=area(ic,2)*                                              &
               0.5*(zp(iloc,jloc,k+1)-zp(iloc,jloc,k-1))
!
!-----------------------------------------------------------------------
!
!    Set icmin, icmax
!
!    1 index is added to the max and 1 index is subtracted from
!    the min so that the comparisons in do loop 180 accept cells
!    that touch rather than just overlap.  Adding now saves time
!    of adding inside the IF statement.
!
!-----------------------------------------------------------------------
!
    DO j=1,ny
      icmin(jc,j,1)=icmin(ic,j,2)-1
      icmax(jc,j,1)=icmax(ic,j,2)+1
    END DO
!
    CYCLE
!
!-----------------------------------------------------------------------
!
!    Match found.
!    Update weighted sums
!
!-----------------------------------------------------------------------
!
    250     CONTINUE
    wgt(jc,1)=wgt(jc,1)+wgt(ic,2)
    xcent(jc,1)=xcent(jc,1)+wgt(ic,2)*xcent(ic,2)
    ycent(jc,1)=ycent(jc,1)+wgt(ic,2)*ycent(ic,2)
    zcent(jc,1)=zcent(jc,1)+wgt(ic,2)*zcent(ic,2)
    qrmax(jc,1)=AMAX1(qrmax(jc,1),qrmax(ic,2))
    iloc=nint((xcent(ic,2)/dx) + 1.5)
    jloc=nint((ycent(ic,2)/dy) + 1.5)
    area(jc,1)=area(jc,1)+area(ic,2)*                                   &
               0.5*(zp(iloc,jloc,k+1)-zp(iloc,jloc,k-1))
!
!-----------------------------------------------------------------------
!
!    Match found.
!    Update icmin,icmax
!
!-----------------------------------------------------------------------
!
    DO j=1,ny
      icmin(jc,j,1)=MIN(icmin(jc,j,1),(icmin(ic,j,2)-1))
      icmax(jc,j,1)=MAX(icmax(jc,j,1),(icmax(ic,j,2)+1))
    END DO
  END DO
!
  RETURN
END SUBROUTINE linkcell
!
!
!##################################################################
!##################################################################
!######                                                      ######
!######               SUBROUTINE MATCTIM                     ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE matctim(maxcell,maxtime,ktime,mtime,mpast,                   & 1
           disthr,time,umove,vmove,                                     &
           ncell,xcent,ycent,zcent,                                     &
           area,wgt,                                                    &
           cname,nctim,alltime,                                         &
           xpos,ypos,zpos,xnot,ynot,                                    &
           allvol,allwgt,                                               &
           umotion,vmotion,wmotion)
!
  IMPLICIT NONE
  INTEGER :: maxcell,maxtime,ktime,mtime
  REAL :: time,disthr,umove,vmove
  INTEGER :: ncell
  REAL :: xcent(maxcell),ycent(maxcell),zcent(maxcell),                 &
       wgt(maxcell),area(maxcell)
  CHARACTER (LEN=2) :: cname(maxcell)
  INTEGER :: nctim(maxtime)
  REAL :: alltime(maxtime)
  REAL :: xpos(maxtime,maxcell),ypos(maxtime,maxcell),                  &
       zpos(maxtime,maxcell)
  REAL :: xnot(maxcell),ynot(maxcell)
  REAL :: allvol(maxtime,maxcell),allwgt(maxtime,maxcell)
  REAL :: umotion(maxtime,maxcell),vmotion(maxtime,maxcell),            &
       wmotion(maxtime,maxcell)
!
!-----------------------------------------------------------------------
!
!  Misc. local variables
!
!-----------------------------------------------------------------------
!
  INTEGER :: icell,jcell,inear,mpast
  REAL :: xxtrap,yxtrap
  REAL :: d2thr,deltat,ddx,ddy,dist2,distmin,d2rat,drmin
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!-----------------------------------------------------------------------
!
!  Index control:   ktime is count of times processed here
!                   mtime is current index in arrays saving locations
!                   mpast is past index in arrays saving locations
!
!-----------------------------------------------------------------------
!
  ktime=ktime+1
  mtime=MOD(ktime,maxtime)
  IF(mtime == 0) mtime=maxtime
  mpast=MOD((ktime-1),maxtime)
  IF(mpast == 0) mpast=maxtime
  PRINT *, ' inside matctim '
  PRINT *, ' ncell,ktime,mtime = ',ncell,ktime,mtime
  PRINT *, ' mpast,nctim(mpast)= ',mpast,nctim(mpast)
!
  d2thr=disthr*disthr
  alltime(mtime)=time
!
  DO icell=1,maxcell
    nctim(mtime)=0
    xpos(mtime,icell)=-1.0E30
    ypos(mtime,icell)=-1.0E30
    zpos(mtime,icell)=-1.0E30
    allvol(mtime,icell)=-1.0E30
    allwgt(mtime,icell)=-1.0E30
    umotion(mtime,icell)=-1.0E30
    vmotion(mtime,icell)=-1.0E30
    wmotion(mtime,icell)=-1.0E30
  END DO
!
  IF(ktime == 1) THEN
    nctim(mtime)=ncell
    DO icell=1,ncell
      xpos(mtime,icell)=xcent(icell)
      ypos(mtime,icell)=ycent(icell)
      zpos(mtime,icell)=zcent(icell)
      allvol(mtime,icell)=area(icell)
      allwgt(mtime,icell)=wgt(icell)
    END DO
  ELSE      ! not first time, need to match up with past cells
!
!-----------------------------------------------------------------------
!
!    Try to match location of "current" cells with extrapolated
!    positions of past cells.  It is assumed throughout that we are
!    dealing with absolute (not grid relative locations).  This
!    is necessary for the linear fit to work right under the condition
!    that the grid motion (umove,vmove) are changing with time.
!
!
!    For extrapolation use, (in the following order of priority,
!     according to the available data for each cell):
!
!   1) Least-squares line from previous u,v fitting
!      (umotion,vmotion and xnot,ynot)
!
!   2) Last umotion and vmotion and last location
!
!   3) Guess motion (zero relative to grid, e.g.)
!      and last location
!
!-----------------------------------------------------------------------
!
    nctim(mtime)=nctim(mpast)
    deltat=time-alltime(mpast)
    DO jcell=1,nctim(mpast)
      IF(umotion(mpast,jcell) /= -1.0E30 .AND.                          &
            vmotion(mpast,jcell) /= -1.0E30 ) THEN
        PRINT *, ' Past abs umotion,vmotion (m/s):',                    &
                  umotion(mpast,jcell),                                 &
                  vmotion(mpast,jcell)
        IF(xnot(jcell) /= -1.0E30 .AND. ynot(jcell) /= -1.0E30 ) THEN
          xxtrap=xnot(jcell)+time*umotion(mpast,jcell)
          yxtrap=ynot(jcell)+time*vmotion(mpast,jcell)
        ELSE
          xxtrap=xpos(mpast,jcell)+deltat*umotion(mpast,jcell)
          yxtrap=ypos(mpast,jcell)+deltat*vmotion(mpast,jcell)
        END IF
      ELSE
        xxtrap=xpos(mpast,jcell)+deltat*umove
        yxtrap=ypos(mpast,jcell)+deltat*vmove
      END IF
      PRINT *, ' past x,xxtrap = ',xpos(mpast,jcell),xxtrap
      PRINT *, ' past y,yxtrap = ',ypos(mpast,jcell),yxtrap
!
!-----------------------------------------------------------------------
!
!    Find the location among the current cells (icell) that is the
!    closest to this past cell (jcell).  When a cell has been used
!    already, its wgt is set to -1.
!
!-----------------------------------------------------------------------
!
      inear=0
      drmin=1.0E30
      DO icell=1,ncell
        IF(wgt(icell) > 0.) THEN
          ddx=xcent(icell)-xxtrap
          ddy=ycent(icell)-yxtrap
          dist2=ddx*ddx + ddy*ddy
          IF(dist2 < d2thr) THEN
            d2rat=dist2/area(icell)
!            d2rat=dist2/wgt(icell)
            IF(d2rat < drmin) THEN
              inear=icell
              drmin=d2rat
            END IF
          END IF
        END IF
      END DO
      IF(inear > 0) THEN
        PRINT *, ' time match jcell,inear= ',jcell,inear
        distmin=SQRT(drmin*area(inear))
        PRINT *, ' min distance (m) = ',distmin
        xpos(mtime,jcell)=xcent(inear)
        ypos(mtime,jcell)=ycent(inear)
        zpos(mtime,jcell)=zcent(inear)
        allvol(mtime,jcell)=area(inear)
        allwgt(mtime,jcell)=wgt(inear)
        wgt(inear)=-1.
      END IF
    END DO
!
!-----------------------------------------------------------------------
!
!    Now make new locations for any cells not matched in the
!    above process.
!
!-----------------------------------------------------------------------
!
    DO icell=1,ncell
      IF(wgt(icell) > 0.) THEN
        PRINT *, ' New cell in time/motion matching: ', icell
        IF(nctim(mtime) < maxcell) THEN
          nctim(mtime)=nctim(mtime)+1
          PRINT *, ' New count is ',nctim(mtime)
          jcell=nctim(mtime)
          xpos(mtime,jcell)=xcent(icell)
          ypos(mtime,jcell)=ycent(icell)
          zpos(mtime,jcell)=zcent(icell)
          allvol(mtime,jcell)=area(icell)
          allwgt(mtime,jcell)=wgt(icell)
        ELSE
          PRINT *, ' WARNING ran out space for cells'
          PRINT *, '    in subroutine matctim. '
        END IF
      END IF
    END DO
  END IF

!
  PRINT *, ' Leaving matctim, nctim = ',nctim(mtime)
  RETURN
END SUBROUTINE matctim
!
!
!##################################################################
!##################################################################
!######                                                      ######
!######                 SUBROUTINE GETVEC                    ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE getvec(maxcell,maxtime,ktime,mtime,                          & 1,4
           cname,nctim,alltime,                                         &
           xpos,ypos,zpos,xnot,ynot,                                    &
           allvol,allwgt,                                               &
           umotion,vmotion,wmotion,growth)
!
  IMPLICIT NONE
  INTEGER :: maxcell,maxtime,ktime,mtime
  CHARACTER (LEN=2) :: cname(maxcell)
  INTEGER :: nctim(maxtime)
  REAL :: alltime(maxtime)
  REAL :: xpos(maxtime,maxcell),ypos(maxtime,maxcell),                  &
       zpos(maxtime,maxcell)
  REAL :: xnot(maxcell),ynot(maxcell)
  REAL :: allvol(maxtime,maxcell),allwgt(maxtime,maxcell)
  REAL :: umotion(maxtime,maxcell),vmotion(maxtime,maxcell),            &
       wmotion(maxtime,maxcell),growth(maxtime,maxcell)
!
!-----------------------------------------------------------------------
!
!  Parameters and variables for least squares fitting.
!
!-----------------------------------------------------------------------
!
  REAL :: twparm
  PARAMETER (twparm=900.*900.)    ! secs*secs
  REAL :: tvalid,ufit,vfit
  INTEGER :: ireturn
!
!-----------------------------------------------------------------------
!
!  Misc. local variables
!
!-----------------------------------------------------------------------
!
  INTEGER :: icell,mpast,itime
  REAL :: deltat,znot,dummy
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  itime=MIN(ktime,maxtime)
  mpast=MOD((ktime-1),maxtime)
  IF(mpast == 0) mpast=maxtime
  PRINT *, '  Inside getvec: ktime= ',ktime
  PRINT *, '     mtime ',mtime,'  nctim(mtime) = ',nctim(mtime)
  PRINT *, '     mpast ',mpast,'  nctim(mpast) = ',nctim(mpast)
!
!-----------------------------------------------------------------------
!
!  Initialize all motions to a wacky missing value.
!
!-----------------------------------------------------------------------
!
  DO icell=1,maxcell
    xnot(icell)=-1.0E30
    ynot(icell)=-1.0E30
    umotion(mtime,icell)=-1.0E30
    vmotion(mtime,icell)=-1.0E30
    wmotion(mtime,icell)=-1.0E30
    growth(mtime,icell)=-1.0E30
  END DO
!
!-----------------------------------------------------------------------
!
!   Use the least squares fit of a line to "itime" number of
!   points to get the velocity.
!
!   In the event a solution could not be found, ireturn.ne.0,
!   use the current position less the past position to get the
!   velocity.
!
!-----------------------------------------------------------------------
!
  IF(ktime > 1) THEN
    deltat=alltime(mtime)-alltime(mpast)
    tvalid=0.5*(alltime(mtime)+alltime(mpast))
    DO icell=1,nctim(mtime)
      IF(xpos(mtime,icell) /= -1.0E30) THEN
        CALL leastsq(xpos(1,icell),alltime,                             &
            tvalid,twparm,itime,xnot(icell),ufit,ireturn)
        IF(ireturn == 0) THEN
          umotion(mtime,icell)=ufit
        ELSE IF(xpos(mpast,icell) /= -1.0E30) THEN
          umotion(mtime,icell)=                                         &
              (xpos(mtime,icell)-xpos(mpast,icell))/deltat
        END IF
        CALL leastsq(ypos(1,icell),alltime,                             &
            tvalid,twparm,itime,ynot(icell),vfit,ireturn)
        IF(ireturn == 0) THEN
          vmotion(mtime,icell)=vfit
        ELSE IF(ypos(mpast,icell) /= -1.0E30) THEN
          vmotion(mtime,icell)=                                         &
              (ypos(mtime,icell)-ypos(mpast,icell))/deltat
        END IF
        CALL leastsq(zpos(1,icell),alltime,                             &
                 tvalid,twparm,itime,                                   &
                 znot,wmotion(mtime,icell),ireturn)
        IF(ireturn /= 0 .AND. zpos(mpast,icell) /= -1.0E30) THEN
          wmotion(mtime,icell)=                                         &
              (zpos(mtime,icell)-zpos(mpast,icell))/deltat
        END IF
        CALL leastsq(allvol(1,icell),alltime,                           &
                 tvalid,twparm,itime,                                   &
                 dummy,growth(mtime,icell),ireturn)
        IF(ireturn /= 0 .AND. allvol(mpast,icell) /= -1.0E30) THEN
          growth(mtime,icell)=                                          &
              (allvol(mtime,icell)-allvol(mpast,icell))/deltat
        END IF
      END IF
    END DO
  END IF
  RETURN
END SUBROUTINE getvec
!
!##################################################################
!##################################################################
!######                                                      ######
!######                SUBROUTINE LEASTSQ                    ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE leastsq(xo,TO,tvalid,twparm,npts,xnot,uo,ireturn) 4
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Modified least-squares estimate for the following equations:
!
!                  X = xnot + Uo*t
!
!  Modification is that what is minimized is not the X deviation
!  squared, rather X deviation squared times a weight specified
!  in the time weighting function tweight.  The idea is to
!  make the line fit the observations near the time specified as
!  the "valid time", tvalid.
!
!-----------------------------------------------------------------------
!
!  AUTHOR:  Steven Lazarus
!           4/12/93
!
!  MODIFICATIONS:
!    4/13/93  Keith Brewster
!             Added weighting based on time.
!             Changed to do one variables at a time.
!    4/19/93  Keith Brewster
!             Added ireturn variable.
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!    npts     Number of points in the Sample
!    xo       Observed storm position (predictand) for X = xnot + Uo*t
!    to       Observed time for xo(to),yo(to) (t=to)
!    tvalid   "Valid" time of estimation, used for time weighting
!             Closer fit to data near vaild time.
!    twparm   Time weighting parameter used as denominator in exponential
!             time weighting function  (secs*secs)
!
!   OUTPUT:
!
!    xnot   Origin of regression line    X =Xnot + Uo*t, (f(xo,to))
!    uo     Slope of the regression line X =Xnot + Uo*t, (f(xo,to))
!    ireturn  Return status indicator.
!
!-----------------------------------------------------------------------
!
!   Variable Declarations
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE

  INTEGER :: npts       ! Number of points in the sample

  REAL :: xo  (npts)    ! Observed independent variable
  REAL :: TO  (npts)    ! Observed time for xo,yo
  REAL :: tvalid        ! valid time
  REAL :: twparm        ! time weighting parameter (secs*secs)
  REAL :: xnot          ! Least-Squares fit for x=xnot+uo*t
  REAL :: uo            ! Least-Squares fit for x=xnot+uo*t
  INTEGER :: ireturn    ! Return status indicator
                        ! =0  succesful completion
                        ! =1  not enough data
                        ! =2  divide by zero avoided, nothing computed

!
!-----------------------------------------------------------------------
!
!   Misc. Local Variables
!
!-----------------------------------------------------------------------
!
  REAL :: sumx,sumt,sumw,sumt2
  REAL :: sumxt
  REAL :: dt,twgt
  REAL :: denomx,denomu

  INTEGER :: i,knt
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!   Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  knt    = 0
  sumw   = 0.0
  sumx   = 0.0
  sumt   = 0.0
  sumxt  = 0.0
  sumt2  = 0.0

  DO i=1,npts
    IF(xo(i) /= -1.0E30) THEN
      knt=knt+1
      dt=tvalid-TO(i)
      twgt=EXP(-dt*dt/twparm)
      sumw  = sumw  + twgt
      sumx  = sumx  + xo(i)*twgt
      sumt  = sumt  + TO(i)*twgt
      sumt2 = sumt2 + TO(i)*TO(i)*twgt
      sumxt = sumxt + xo(i)*TO(i)*twgt
    END IF
  END DO

  IF( knt > 2) THEN
    denomx = sumt*sumt - sumt2*sumw
    denomu = sumw*sumt2 - sumt*sumt
    IF( denomx /= 0.  .AND.  denomu /= 0.) THEN
      ireturn=0
      xnot  = (sumt*sumxt - sumt2*sumx) / denomx
      uo    = (sumw*sumxt - sumx*sumt)  / denomu
    ELSE
      ireturn=2
      xnot  = -1.e30
      uo    = -1.e30
    END IF
  ELSE
    ireturn=1
    xnot  = -1.e30
    uo    = -1.e30
  END IF

  RETURN
END SUBROUTINE leastsq