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


SUBROUTINE vrhole(nx,ny,nz,x,y,z,zp,                                    & 1,3
           u,v,w,                                                       &
           vrom,ubar,vbar,ptol,                                         &
           tem4,tem5,tem6,tem7,tem8,tem9)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:This program is designed to fill-in the holes for missing
!          single-Doppler data on a Cartesian grid.
!
!-----------------------------------------------------------------------
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Steve Weygandt
!
!
!  MODIFICATION HISTORY:
!
!  1/96 (Steve Lazarus)
!  Created this sbroutine from a block of code generated by SW.
!
!  23/02/96 (Limin Zhao)
!  Adapted the code into the data assimilation package.
!
!  06/03/96 (Limin Zhao)
!  Added the include file 'assim.inc'
!
!  08/03/96 (Limin Zhao)
!  Added checks for processing real data.
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
!  INPUT ARRAYS:
!
!
!    vrom     Observed radial velocity on model grid scalar point.
!
!    ubar     Base-state u wind component
!    vbar     Base-state v wind component
!
!    nx       Number of grid points in the x-direction (east/west)
!    ny       Number of grid points in the y-direction (north/south)
!    nz       Number of grid points in the z-direction (vertical)
!
!    dx       Grid spacing in the x-direction (east/west)
!    dy       Grid spacing in the y-direction (north/south)
!    dz       Grid spacing in the z-direction (vertical)
!
!  OUTPUT ARRAYS:
!
!    tem9     Output hole-filled vr (observed)
!
!
!  WORK ARRAYS:
!
!    tem4
!    tem5
!    tem6
!    tem7
!    tem8
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
!  Variable Declarations:
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE             ! Force explicit declarations

  INTEGER :: nx, ny, nz        ! Number of grid points in 3 directions

  REAL :: u(nx,ny,nz)          !Note: u,v,w are working arrays here.
  REAL :: v(nx,ny,nz)
  REAL :: w(nx,ny,nz)

  REAL :: vrom  (nx,ny,nz)     ! Observed radial velocity on model grid

  REAL :: ubar  (nx,ny,nz)     ! U component mean wind
  REAL :: vbar  (nx,ny,nz)     ! V component mean wind

  REAL :: tem4  (nx,ny,nz)     ! Work array
  REAL :: tem5  (nx,ny,nz)     ! Work array
  REAL :: tem6  (nx,ny,nz)     ! Work array
  REAL :: tem7  (nx,ny,nz)     ! Work array
  REAL :: tem8  (nx,ny,nz)     ! Work array
  REAL :: tem9  (nx,ny,nz)     ! Hole-filled vr

  REAL :: x     (nx)           ! Work array
  REAL :: y     (ny)           ! Work array
  REAL :: z     (nz)           ! Work array
  REAL :: zp    (nx,ny,nz)     ! Work array

  REAL :: assimtim (100)

  REAL :: ptol

!
!-----------------------------------------------------------------------
!
!  Miscellaneous local variables:
!
!-----------------------------------------------------------------------
!
  INTEGER :: i,j,k             ! Do-loop indices

  INTEGER :: count

  REAL :: rad,xs,ys,zs
  REAL :: xmove,ymove

!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------
!
  INCLUDE 'globcst.inc'     ! Model control constants
  INCLUDE 'assim.inc'       ! Assim/Retr control parameters
  INCLUDE 'grid.inc'
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  WRITE(6,*) 'code in vrhole; spval=',spval
  WRITE(6,*) 'adas mean wind ubar(5,5,5),vbar(5,5,5): ',                &
                             ubar(5,5,5),vbar(5,5,5)
  WRITE(6,*) 'xshift,yshift,zshift: ', xshift,yshift,zshift

  umove=0.0
  vmove=0.0
  xmove= xshift - umove*(curtim-assimtim(1))
  ymove= yshift - vmove*(curtim-assimtim(1))

  count = 0

  DO k=1,nz
    DO j=1,ny
      DO i=1,nx
        tem5(i,j,k) = 0.0
        tem4(i,j,k) = 0.0
        tem6(i,j,k) = 0.0
        tem7(i,j,k) = 0.0
        tem8(i,j,k) = 0.0
        tem9(i,j,k) = 0.0
      END DO
    END DO
  END DO

  DO k=1,nz-1
    DO j=1,ny-1
      DO i=1,nx-1
        tem6(i,j,k) = 0.0
        tem7(i,j,k) = 0.0
        tem8(i,j,k) = 0.0
        xs    = 0.5*(x(i)+x(i+1)) - xmove
        ys    = 0.5*(y(j)+y(j+1)) - ymove
        zs    = 0.5*(z(k)+z(k+1)) - zshift
        rad   = SQRT(xs**2+ys**2+zs**2)
        IF(vrom(i,j,k) /= spval) THEN
          tem5(i,j,k) = vrom(i,j,k)*(xs/rad)-u(i,j,k)
          tem4(i,j,k) = 0.0
        ELSE IF(i == 1.OR.i == nx-1.OR.j == 1.OR.j == ny-1) THEN
          tem5(i,j,k) = 0.0
          tem4(i,j,k) = 0.0
        ELSE                          ! Fill vru outside rain regions
          tem5(i,j,k) = 0.0
          tem4(i,j,k) = spval
          count = count + 1
        END IF
      END DO
    END DO
  END DO

  PRINT *,'On call to POIS3D, there are ',count,' filled vrU values'


  CALL pois3d(nx,ny,nz,dx,dy,dz,ptol,2.0,tem8,tem4,tem5,tem6,tem7)


  DO k=1,nz-1
    DO j=1,ny-1
      DO i=1,nx-1
        xs    = 0.5*(x(i)+x(i+1)) - xmove
        ys    = 0.5*(y(j)+y(j+1)) - ymove
        zs    = 0.5*(z(k)+z(k+1)) - zshift
        rad   = SQRT(xs**2+ys**2+zs**2)
        tem9(i,j,k) = tem5(i,j,k)*xs/rad + tem9(i,j,k)
      END DO
    END DO
  END DO

  count = 0
  DO k=1,nz-1
    DO j=1,ny-1
      DO i=1,nx-1
        tem6(i,j,k) = 0.0
        tem7(i,j,k) = 0.0
        tem8(i,j,k) = 0.0
        xs    = 0.5*(x(i)+x(i+1)) - xmove
        ys    = 0.5*(y(j)+y(j+1)) - ymove
        zs    = 0.5*(z(k)+z(k+1)) - zshift
        rad   = SQRT(xs**2+ys**2+zs**2)
        IF(vrom(i,j,k) /= spval) THEN
          tem5(i,j,k) = vrom(i,j,k)*(ys/rad)-v(i,j,k)
          tem4(i,j,k) = 0.0
        ELSE IF(i == 1.OR.i == nx-1.OR.j == 1.OR.j == ny-1) THEN
          tem5(i,j,k) = 0.0
          tem4(i,j,k) = 0.0
        ELSE                         ! Fill vrV outside rain regions
          tem5(i,j,k) = 0.0
          tem4(i,j,k) = spval
          count = count + 1
        END IF
      END DO
    END DO
  END DO

  PRINT *,'On call to POIS3D, there are ',count,' filled vrV values'

  CALL pois3d(nx,ny,nz,dx,dy,dz,ptol,2.0,tem8,tem4,tem5,tem6,tem7)

  DO k=1,nz-1
    DO j=1,ny-1
      DO i=1,nx-1
        xs    = 0.5*(x(i)+x(i+1)) - xmove
        ys    = 0.5*(y(j)+y(j+1)) - ymove
        zs    = 0.5*(z(k)+z(k+1)) - zshift
        rad   = SQRT(xs**2+ys**2+zs**2)
        tem9(i,j,k) = tem5(i,j,k)*ys/rad + tem9(i,j,k)
      END DO
    END DO
  END DO


  count = 0
  DO k=1,nz-1
    DO j=1,ny-1
      DO i=1,nx-1
        tem6(i,j,k) = 0.0
        tem7(i,j,k) = 0.0
        tem8(i,j,k) = 0.0
        xs    = 0.5*(x(i)+x(i+1)) - xmove
        ys    = 0.5*(y(j)+y(j+1)) - ymove
        zs    = 0.5*(z(k)+z(k+1)) - zshift
        rad   = SQRT(xs**2+ys**2+zs**2)
        IF(vrom(i,j,k) /= spval) THEN
          tem5(i,j,k) = vrom(i,j,k)*(zs/rad)
          tem4(i,j,k) = 0.0
        ELSE IF(k == 2.OR.k == nz-1) THEN
          tem5(i,j,k) = 0.0
          tem4(i,j,k) = 0.0
        ELSE                          ! Fill w outside rain regions
          tem5(i,j,k) = 0.0
          tem4(i,j,k) = spval
          count = count + 1
        END IF
      END DO
    END DO
  END DO

  PRINT *,'On call to POIS3D, there are ',count,' filled vrW values'

  CALL pois3d(nx,ny,nz,dx,dy,dz,ptol,2.0,tem8,tem4,tem5,tem6,tem7)

  DO k=1,nz-1
    DO j=1,ny-1
      DO i=1,nx-1
        rad   = SQRT(xs**2+ys**2+zs**2)
        xs    = 0.5*(x(i)+x(i+1)) - xmove
        ys    = 0.5*(y(j)+y(j+1)) - ymove
        zs    = 0.5*(z(k)+z(k+1)) - zshift
        rad   = SQRT(xs**2+ys**2+zs**2)
        tem9(i,j,k) = tem5(i,j,k)*zs/rad + tem9(i,j,k)
      END DO
    END DO
  END DO

  DO k=1,nz-1
    DO j=1,ny-1
      DO i=1,nx-1
        xs    = 0.5*(x(i)+x(i+1)) - xmove
        ys    = 0.5*(y(j)+y(j+1)) - ymove
        zs    = 0.5*(z(k)+z(k+1)) - zshift
        rad   = SQRT(xs**2+ys**2+zs**2)
        tem9(i,j,k) = tem9(i,j,k) + u(i,j,k)*(xs/rad)                   &
                                 + v(i,j,k)*(ys/rad)
        vrom(i,j,k) = tem9(i,j,k)
      END DO
    END DO
  END DO


  RETURN

END SUBROUTINE vrhole