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


SUBROUTINE read_station(infile,mxstalo,latsta,lonsta,                   & 2
           nstatyp,nstapro,nsta,sname,state,sitena,nelev)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  This subroutine will read external staion information.
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
!  AUTHOR:
!    Min Zou (6/1/97)
!
!  Modification history:
!
!-----------------------------------------------------------------------
!
!  Variable Declarations
!
!-----------------------------------------------------------------------
!
  CHARACTER (LEN=*) :: infile
  INTEGER :: nsta,nstapro(mxstalo),nstatyp(mxstalo)
  REAL :: latsta(mxstalo), lonsta(mxstalo)
  CHARACTER (LEN=5) :: sname(mxstalo)
  CHARACTER (LEN=2) :: state(mxstalo)
  CHARACTER (LEN=20) :: sitena(mxstalo)
  INTEGER :: nelev(mxstalo)
  CHARACTER (LEN=132) :: line

  OPEN(1,IOSTAT=ios,FILE=infile,STATUS='old',                           &
        FORM='formatted')
  IF(ios /= 0) THEN     ! error during read
    istatus = -1
    WRITE(6,650) infile
    650       FORMAT(' +++ ERROR opening: ',a70,' +++')
    WRITE(6,651) ios
    651       FORMAT('     IOS code = ',i5)
    RETURN
  END IF

  nsta = 0

! Read only lines that begin with A-Z, a-z, or 0-9 -- treat the rest as
! comments

  DO i=1,mxstalo
    READ(1,'(A)',END=999) line
    IF ( ( line(:1) >= 'A' .AND. line(:1) <= 'Z' ) .OR.                 &
           ( line(:1) >= 'a' .AND. line(:1) <= 'z' ) .OR.               &
           ( line(:1) >= '0' .AND. line(:1) <= '9' ) ) THEN
      nsta = nsta + 1
      READ(line,101,ERR=999)sname(nsta),state(nsta),sitena(nsta),       &
          latsta(nsta),lonsta(nsta),nelev(nsta),nstatyp(nsta),          &
          nstapro(nsta)
    END IF
    101   FORMAT(a5,2X,a2,1X,a20,4X,f8.3,1X,f8.3,1X,i5,2X,i2,i1)
    102   FORMAT(a5,2X,a2,1X,a20,4X,f8.3,1X,f8.3,1X,i5,2X,i2,1X,i1)
  END DO
  999   CONTINUE
  CLOSE(1)

  RETURN
END SUBROUTINE read_station
!
!##################################################################
!##################################################################
!######                                                      ######
!######                SUBROUTINE INTERP_P                   ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE interp_p (pbar,zpc,ibgn,iend,nx,jbgn,jend,ny,                & 1
           kbgn,kend,nz)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  This subroutine will interpolate pressure for draw pressure bar.
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
!  AUTHOR:
!    Min Zou (6/1/97)
!
!  Modification history:
!
!-----------------------------------------------------------------------
!
!  Variable Declarations
!
!-----------------------------------------------------------------------
!

  REAL :: pbar(nx,ny,nz), zpc(nx,ny,nz)
  INTEGER :: presaxis_no
  REAL :: pres_val(20), pres_z(20)
  COMMON /pressbar_par/presaxis_no,pres_val,pres_z
  REAL :: tmp
  REAL :: pres_val1(20)
  REAL :: pz(100), pb(100)
! estimate a value for pz(kbgn-1) and pb(kbgn-1)
  kbgn1 = kbgn-1
  pb(kbgn1)=1010.
  pz(kbgn1)=0.5
! get max pressure
  DO k=kbgn,kend
    pb(k)=0.
    pz(k)=0.
    ini = 0
    inj = 0
    DO j=jbgn, jend
      DO i=ibgn, iend
        IF(0.01*pbar(i,j,k) > pb(k)) THEN
          ini = i
          inj = j
        END IF
      END DO
    END DO
    IF(ini /= 0 .AND. inj /= 0) THEN
      pb(k) = 0.01*pbar(ini,inj,k)
      pz(k) = zpc(ini, inj,k)
    ELSE                     !! missing values
      WRITE(6,'(a,i2,a)')                                               &
          'Warning: Missing pressure value on level ',k,                &
          ' Using previous level instead of.'
      pb(k) = pb(k-1)
      pz(k) = pz(k-1)
    END IF
  END DO

  k=0
  DO j=1,presaxis_no
    IF (pres_val(j) > pb(kbgn) ) THEN
      k  = 1
!       print*,'pres_val(j)',pres_val(j)
!       print*, '>= pb(kbgn)',pb(kbgn)
      pres_val1(k) = pb(kbgn)
      pres_z(k) = pz(kbgn)
    END IF
  END DO

  DO j=1,presaxis_no
    DO i=kbgn,kend-1
      tmp = pres_val(j)
      IF(tmp <= pb(i).AND.tmp > pb(i+1))THEN !find pressure interpolate
        k = k+1
        a1 = ALOG(pb(i))-ALOG(tmp)
        a2 = ALOG(pb(i))-ALOG(pb(i+1))
        a3 = pz(i+1)-pz(i)
        pres_z(k) = pz(i)+a3*a1/a2
        pres_val1(k) = pres_val(j)
      END IF
    END DO
  END DO
  k=k+1
  pres_val1(k) = pb(kend)
  pres_z(k) = pz(kend)


  presaxis_no = k
  DO i=1,presaxis_no
    pres_val(i) = pres_val1(i)
  END DO
  RETURN
END SUBROUTINE interp_p
!
!##################################################################
!##################################################################
!######                                                      ######
!######                SUBROUTINE GET_MULOVRLAY              ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE get_mulovrlay(var,LEN,num,ovrname,sovrlay) 78
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  This subroutine will find out the overlay multiple plots.
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
!  AUTHOR:
!    Min Zou (6/1/97)
!
!  Modification history:
!
!-----------------------------------------------------------------------
!
!  Variable Declarations
!
!-----------------------------------------------------------------------
!
  INTEGER :: num, sovrlay
  CHARACTER (LEN=*) :: var
  CHARACTER (LEN=*) :: ovrname(10)

  INTEGER :: i,LEN

  sovrlay = 0

  DO i = 1,num
    IF(var(1:LEN) == ovrname(i)(1:LEN) ) sovrlay=1
  END DO

  RETURN
END SUBROUTINE get_mulovrlay

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


SUBROUTINE arps_ct (cct,nx,ny,nz,p,t,td,ppccl,wrk1,wrk2,wrk3) 1,2

!-----------------------------------------------------------------------
!  Purpose:
!  Calculate the convective temperature (celsius)
!
!  AUTHOR:  Min Zou
!    07/10/1997
!
!  MODIFICATIONS:
!
!
!-----------------------------------------------------------------------
!
!  INPUT:
!    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 vertical
!    p        pressure (Pascals)
!    t        temperature(degrees Kelvin)
!    td       dew-point temperature (degrees Kelvin)
!
!  OUTPUT:
!
!    cct      convective temperature (Celsius)
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE

  INTEGER :: nx,ny,nz
  REAL :: p(nx,ny,nz)    ! pressure (Pa)
  REAL :: t(nx,ny,nz)    ! temperature(K)
  REAL :: td(nx,ny,nz)   ! dew-point temperature (K)

  REAL :: cct(nx,ny)     ! convective temperature (C)

!-----------------------------------------------------------------------
!
!  Misc. temporary variables
!
!-----------------------------------------------------------------------

  REAL :: ppccl(nx,ny)   ! pressure (millibars) at the convective
                         !condensation level

  REAL :: wrk1(nz), wrk2(nz), wrk3(nz)

!-----------------------------------------------------------------------
!
!  Misc. local variables
!
!-----------------------------------------------------------------------
!
  INTEGER :: i,j
  REAL :: pm             ! pressure(millibars) at upper boundary of the
                         ! layer for computing the mean mixing ratio.
  REAL :: mrbar          ! mean mixing ratio (g/kg) in the layer bounded
                         ! by pressures at the p bottom and the pm at the top
!
!-----------------------------------------------------------------------
!
!  Function f_pccl and f_ct and inline directive for Cray PVP
!
!-----------------------------------------------------------------------
!
  REAL :: f_pccl, f_ct

!fpp$ expand (f_pccl)
!fpp$ expand (f_ct)
!dir$ inline always f_pccl,f_ct
!*$*  inline routine (f_pccl,f_ct)

!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!

  DO j = 1,ny
    DO i = 1,nx
      pm = p(i,j,2) - 5000.

      ppccl(i,j) = f_pccl(pm,p(i,j,1),t(i,j,1),td(i,j,1),mrbar,nz)

      cct(i,j) = f_ct(mrbar, ppccl(i,j), p(i,j,2))

    END DO
  END DO

  RETURN
END SUBROUTINE arps_ct


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


SUBROUTINE temper ( nx,ny,nz,theta, ppert, pbar, t ) 6

!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!    Using a version of Poisson's formula, calculate temperature.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Joe Bradley
!    12/05/91
!
!  MODIFICATIONS:
!    Modified by Ming Xue so that arrays are only defined at
!             one time level.
!    6/09/92  Added full documentation and phycst include file for
!             rddcp=Rd/Cp  (K. Brewster)
!
!-----------------------------------------------------------------------
!
!  INPUT:
!    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 vertical
!
!    theta    Potential temperature (degrees Kelvin)
!    ppert    Perturbation pressure (Pascals)
!    pbar     Base state pressure (Pascals)
!
!  OUTPUT:
!
!    t        Temperature (degrees Kelvin)
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
!
  INTEGER :: nx,ny,nz
!
  REAL :: theta(nx,ny,nz)      ! potential temperature (degrees Kelvin)
  REAL :: ppert(nx,ny,nz)      ! perturbation pressure (Pascals)
  REAL :: pbar (nx,ny,nz)      ! base state pressure (Pascals)
!
  REAL :: t    (nx,ny,nz)      ! temperature (degrees Kelvin)
!
!-----------------------------------------------------------------------
!
!  Include file
!
!-----------------------------------------------------------------------
!
  INCLUDE 'phycst.inc'
!
!-----------------------------------------------------------------------
!
!  Misc. local variables
!
!-----------------------------------------------------------------------
!
  INTEGER :: i,j,k
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!
!-----------------------------------------------------------------------
!
!  Calculate the temperature using Poisson's formula.
!
!-----------------------------------------------------------------------
!
  DO k=1,nz-1
    DO j=1,ny-1
      DO i=1,nx-1

        t(i,j,k) = theta(i,j,k) *                                       &
             (((ppert(i,j,k) + pbar(i,j,k)) / p0) ** rddcp)

      END DO
    END DO
  END DO

  RETURN
END SUBROUTINE temper

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


SUBROUTINE intepo(ns,xpos,ypos,vpos,m,n,x,y,var) 2
!
!  Does bilinear interpolation of variables (nvar variables)
!  to position (xpos,ypos) given a grid of variables "var", their
!  x and y positions of the grid (x,y) and the index of the
!  grid point (i,j) which is to the "lower-left" of xpos and ypos.
!  That means (xpos,ypos) is between (i,j) and (i+1,j+1).
!
!  A saftey feature in case xpos and ypos are outside the grid
!  is built-in.
!
  IMPLICIT NONE
!
!  Arguments, input
!
  INTEGER :: m,n
  REAL :: var(m,n)
  REAL :: x(m,n),y(m,n)
  INTEGER :: ns
  REAL :: xpos(ns),ypos(ns), vpos(ns)
  INTEGER :: i,j
  REAL :: a1,a2
!
!  Arguments output
!
!  Misc Internal Variables
!
  REAL :: dxpos,dypos,dx,dy
  INTEGER :: k
!

  DO k=1,ns
    DO j=1,n-1
      DO i=1,m-1
        IF( (xpos(k) >= x(i,j) .AND. xpos(k) < x(i+1,j+1)) .AND.        &
              (ypos(k) >= y(i,j) .AND. ypos(k) < y(i+1,j+1)) )  THEN

          IF(var(i,j) /= -9999. .AND. var(i+1,j) /= -9999. .AND.        &
                var(i,j+1) /= -9999. .AND. var(i+1,j+1) /= -9999.) THEN
            dx=x(i+1,j)-x(i,j)
            dy=y(i,j+1)-y(i,j)
            dxpos=(xpos(k)-x(i,j))/dx
            dypos=(ypos(k)-y(i,j))/dy

            a1=(1.-dxpos)*var(i,j) + dxpos*var(i+1,j)
            a2=(1.-dxpos)*var(i,j+1) + dxpos*var(i+1,j+1)

            vpos(k) = (1.-dypos)*a1 + dypos*a2
          ELSE
            vpos(k) = -9999.
          END IF
!
        END IF
      END DO
    END DO
  END DO
  RETURN
END SUBROUTINE intepo

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


SUBROUTINE unigrid(nx,nz,f,z,fdata,zdata,fprof,zprof) 6,1
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  E
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!    nx       first dimension of f
!    nz       second dimension of f
!    f        2-dimension array of variable
!    z        z coordinate of grid points in physcal space (m)
!    fdata    1-D array defined at levels in zdat to be
!               interpolated to levels defined by zpro.
!    zdata    The height of the input data given by fdat.
!    zprof    The grid level height to which data are interpolated
!
!  OUTPUT:
!
!    fprof    The number of interpolated data levels
!    f        2-dimension array of variable
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
!
  INTEGER :: nx,nz
  INTEGER :: i,k
  REAL :: f(nx,nz)
  REAL :: z(nx,nz)
  REAL :: fdata(nz),zdata(nz),fprof(nz),zprof(nz)
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  DO i=1,nx

    DO k=1,nz
      fdata(k)=f(i,k)
      zdata(k)=z(i,k)
    END DO

    CALL inte1d(fdata,zdata,nz,fprof,zprof,nz)

    DO k=1,nz
      f(i,k)=fprof(k)
      IF(zprof(k) < zdata(1) .OR. zprof(k) > zdata(nz) ) f(i,k)=-9999.0
    END DO

  END DO

  RETURN
END SUBROUTINE unigrid
!
!
!##################################################################
!##################################################################
!######                                                      ######
!######                SUBROUTINE FILZERO                    ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE filzero( a, n) 5

!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!    Fill vector a with zeros.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: St. Paul, Dead Sea Scrolls
!
!  MODIFICATIONS:
!    6/09/92  Added full documentation (K. Brewster)
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
!
  INTEGER :: n
  REAL :: a(n)
!
!-----------------------------------------------------------------------
!
!  Misc. local variables
!
!-----------------------------------------------------------------------
!
  INTEGER :: i
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  DO i=1,n
    a(i)=0.0
  END DO

  RETURN
END SUBROUTINE filzero
!
!##################################################################
!##################################################################
!######                                                      ######
!######                SUBROUTINE IFILZERO                   ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE ifilzero( ia, n )

!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!    Fill vector a with zeros.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: St. Paul, Dead Sea Scrolls
!
!  MODIFICATIONS:
!    6/09/92  Added full documentation (K. Brewster)
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
!
  INTEGER :: n
  INTEGER :: ia(n)
!
!-----------------------------------------------------------------------
!
!  Misc. local variables
!
!-----------------------------------------------------------------------
!
  INTEGER :: i
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  DO i=1,n
    ia(i)=0
  END DO

  RETURN
END SUBROUTINE ifilzero

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


SUBROUTINE slength ( string, length ) 27

!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Return the length of the non-blank part of a string.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Ming Xue
!
!  MODIFICATION HISTORY:
!
!  6/09/92 (K. Brewster)
!  Added full documentation and streamlined logic
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!    string   character string to be sized
!
!  INPUT/OUTPUT:
!
!    length   on input, full size of character string
!             on output, true length of string as measured by the
!                        location of the last non-blank character
!
!
!-----------------------------------------------------------------------
!

!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
  INTEGER :: length
  CHARACTER (LEN=*) :: string
!
!-----------------------------------------------------------------------
!
!  Misc. local variables
!
!-----------------------------------------------------------------------
!
  INTEGER :: i
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  DO i = length,1,-1
    IF(string(i:i) /= ' ') EXIT
  END DO
  101   CONTINUE

  length = i

  RETURN
END SUBROUTINE slength
!
!
!
!##################################################################
!##################################################################
!######                                                      ######
!######                 SUBROUTINES cal_t                    ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!  Calculate T value.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Min Zou
!  3/2/98
!
!  MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!


SUBROUTINE cal_t ( tem,tz,nx,ny,nz,                                     & 1
           tob,label,length,units )

  IMPLICIT NONE
  INCLUDE 'arpsplt.inc'

  INTEGER :: ovrobs,obsset,obscol,obs_marktyp
  REAL :: obs_marksz
  COMMON /obspar/ ovrobs,obsset,obscol,obs_marktyp, obs_marksz

  INTEGER :: nobs
  COMMON /sfc_obs1/ nobs
  REAL :: latob(mxsfcob),lonob(mxsfcob)
  REAL :: obs1(mxsfcob),obs2(mxsfcob)
  COMMON /sfc_obs2/ latob,lonob,obs1,obs2

  INTEGER :: nx,ny,nz, length
!   real tem(*), tz(*), tob(*)
  REAL :: tem(nx,ny,nz), tz(nx,ny,nz), tob(*)
  CHARACTER (LEN=*) :: units
  CHARACTER (LEN=*) :: label

  INTEGER :: i,j,k, iob
  INTEGER :: ibgn,iend, jbgn, jend, kbgn,kend, isize,jsize,ksize
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  ibgn = 1
  iend = nx-1
  jbgn = 1
  jend = ny-1
  kbgn = 1
  kend = nz-1
  isize = (iend-ibgn)+1
  jsize = (jend-jbgn)+1
  ksize = (kend-kbgn)+1

  IF (units(1:1) == 'F' .OR. units(1:1) == 'f') THEN
    DO k= kbgn, kend
      DO j= jbgn, jend
        DO i= ibgn, iend
!       ij = i-ibgn+1 + (j-jbgn)*isize + (k-kbgn)*jsize
!       tem(ii) = 32.0 + 1.8*(tz(ii) - 273.15)
          tem(i,j,k) = 32.0 + 1.8*(tz(i,j,k) - 273.15)
        END DO
      END DO
    END DO
    label = 'T (F)'
    length = 5

    IF(ovrobs == 1 .AND. nobs > 0) THEN
      DO iob=1,nobs
        obs1(iob)=tob(iob)
      END DO
      obsset=1
    END IF

  ELSE                 ! default units is C or c
    DO k= kbgn, kend
      DO j= jbgn, jend
        DO i= ibgn, iend
!       ij = i-ibgn+1 + (j-jbgn)*isize + (k-kbgn)*jsize
!       tem(ii) = tz(ii) - 273.15
          tem(i,j,k) = tz(i,j,k) - 273.15
        END DO
      END DO
    END DO
    label = 'T (C)'
    length = 5

    IF(ovrobs == 1 .AND. nobs > 0) THEN
      DO iob=1,nobs
        obs1(iob)=(tob(iob)-32.)*5./9.
      END DO
      obsset=1
    END IF
  END IF

  RETURN
END SUBROUTINE cal_t


!
!
!##################################################################
!##################################################################
!######                                                      ######
!######                 SUBROUTINES cal_vh                   ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!  Calculate vh value.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Min Zou
!  3/2/98
!
!  MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!


SUBROUTINE cal_vh(tem9,u,v,nx,ny,nz,vhunits,label,length,tem8) 1,1

  IMPLICIT NONE
  INTEGER :: nx,ny,nz
  REAL :: tem9(nx,ny,nz),u(nx,ny,nz), v(nx,ny,nz)
  REAL :: tem8(nx,ny,nz)
  INTEGER :: vhunits, length
  CHARACTER (LEN=*) :: label

  INTEGER :: i,j,k, onvf
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  DO k=1,nz
    DO j=1,ny
      DO i=1,nx
        tem8(i,j,k) = SQRT(      END DO
    END DO
  END DO

  onvf = 0
  CALL avgy(tem8 , onvf,                                                &
       nx,ny,nz, 1,nx-1, 1,ny-1, 1,nz-1, tem9)
  IF(vhunits == 1) THEN
    label = 'Horiz. wind (m/s)'
    length =17
  ELSE IF( vhunits == 2) THEN
    label = 'Horiz. wind (kts)'
    length = 17
    DO k=1,nz
      DO j=1,ny
        DO i=1,nx
          tem9(i,j,k) = tem9(i,j,k)*1.943844
        END DO
      END DO
    END DO
  ELSE IF (vhunits == 3) THEN
    label = 'Horiz. wind (MPH)'
    length = 17
    DO k=1,nz
      DO j=1,ny
        DO i=1,nx
          tem9(i,j,k) = tem9(i,j,k)*2.236936
        END DO
      END DO
    END DO
  END IF
  RETURN
END SUBROUTINE cal_vh


!
!
!##################################################################
!##################################################################
!######                                                      ######
!######                 SUBROUTINES cal_qw                   ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!  Calculate qw value.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Min Zou
!  3/2/98
!
!  MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!


SUBROUTINE cal_qw(tem9,qc,qr,qi,qs,qh, nx,ny,nz) 1

  IMPLICIT NONE
  INTEGER :: nx,ny,nz
  REAL :: tem9(nx,ny,nz)
  REAL :: qc(nx,ny,nz), qr(nx,ny,nz), qi(nx,ny,nz), qs(nx,ny,nz)
  REAL :: qh(nx,ny,nz)

  INTEGER :: i,j,k
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  DO k=1,nz
    DO j=1,ny
      DO i=1,nx
        tem9(i,j,k)=qc(i,j,k)+qr(i,j,k)+qi(i,j,k)+qs(i,j,k)+qh(i,j,k)
      END DO
    END DO
  END DO

  RETURN
END SUBROUTINE cal_qw

!
!
!##################################################################
!##################################################################
!######                                                      ######
!######                 SUBROUTINES cal_rh                   ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!  Calculate rh value.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Min Zou
!  3/2/98
!
!  MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!


SUBROUTINE cal_rh(tem9,pt, pprt ,pbar,qv,tem1,tem2,nx,ny,nz) 1,2

  IMPLICIT NONE
  INTEGER :: nx,ny,nz
  REAL :: tem9(nx,ny,nz), tem1(nx,ny,nz), tem2(nx,ny,nz)
  REAL :: pt(nx,ny,nz), pprt(nx,ny,nz), pbar(nx,ny,nz), qv(nx,ny,nz)

  INTEGER :: i,j,k
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!

  CALL temper (nx,ny,nz,pt, pprt ,pbar,tem1)
  CALL getqvs(nx,ny,nz, 1,nx-1,1,ny-1,1,nz-1, pbar,tem1,tem2)

  DO k=1,nz-1
    DO j=1,ny-1
      DO i=1,nx-1
        tem9(i,j,k) = qv(i,j,k)/tem2(i,j,k)
      END DO
    END DO
  END DO

  RETURN
END SUBROUTINE cal_rh
!
!##################################################################
!##################################################################
!######                                                      ######
!######                 SUBROUTINES cal_td                   ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!  Calculate td value.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Min Zou
!  3/2/98
!
!  MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!


SUBROUTINE cal_td(tem9,td,nx,ny,nz,tdunits,label, length) 1

  IMPLICIT NONE
  INTEGER :: nx,ny,nz
  REAL :: tem9(nx,ny,nz), td(nx,ny,nz)
  INTEGER :: length
  CHARACTER (LEN=*) :: label
  CHARACTER (LEN=*) :: tdunits

  INTEGER :: i,j,k
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!

  IF (tdunits(1:1) == 'F' .OR. tdunits(1:1) == 'f') THEN
    DO k=1,nz-1
      DO j=1,ny-1
        DO i=1,nx-1
          tem9(i,j,k) = 32.0 + 1.8*(td(i,j,k) - 273.15)
        END DO
      END DO
    END DO
    label = 'Td (F)'
    length = 6
  ELSE
    DO k=1,nz-1
      DO j=1,ny-1
        DO i=1,nx-1
          tem9(i,j,k) = td(i,j,k) - 273.15
        END DO
      END DO
    END DO
    label = 'Td (C)'
    length = 6
  END IF

  RETURN
END SUBROUTINE cal_td

!
!
!##################################################################
!##################################################################
!######                                                      ######
!######                 SUBROUTINES cal_tdobs                ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!  Calculate td observation value.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Min Zou
!  3/2/98
!
!  MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!


SUBROUTINE cal_tdobs(tdob, tdunits) 1

  IMPLICIT NONE
  INCLUDE 'arpsplt.inc'

  INTEGER :: ovrobs,obsset,obscol,obs_marktyp
  REAL :: obs_marksz
  COMMON /obspar/ ovrobs,obsset,obscol,obs_marktyp, obs_marksz

  INTEGER :: nobs
  COMMON /sfc_obs1/ nobs
  REAL :: latob(mxsfcob),lonob(mxsfcob)
  REAL :: obs1(mxsfcob),obs2(mxsfcob)
  COMMON /sfc_obs2/ latob,lonob,obs1,obs2

  REAL :: tdob(mxsfcob)

  INTEGER :: iob
  CHARACTER (LEN=*) :: tdunits
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  IF (tdunits(1:1) == 'F' .OR. tdunits(1:1) == 'f') THEN
    IF(ovrobs == 1 .AND. nobs > 0) THEN
      DO iob=1,nobs
        obs1(iob)=tdob(iob)
      END DO
      obsset=1
    END IF
  ELSE
    IF(ovrobs == 1 .AND. nobs > 0) THEN
      DO iob=1,nobs
        obs1(iob)=(tdob(iob)-32.)*5./9.
      END DO
      obsset=1
    END IF
  END IF

  RETURN
END SUBROUTINE cal_tdobs

!
!
!##################################################################
!##################################################################
!######                                                      ######
!######                 SUBROUTINES cal_rfc                  ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!  Calculate rfc value.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Min Zou
!  3/2/98
!
!  MODIFICATION HISTORY:
!
!  Ming Xue (10/16/2001)
!  Now passing in precalculated reflectivity field instead of calculating
!  it inside.
!
!-----------------------------------------------------------------------
!


SUBROUTINE cal_rfc(nx, ny, nz, ref, refc) 1

  IMPLICIT NONE
  INTEGER :: nx,ny,nz
  REAL, INTENT(IN ) :: ref (nx,ny,nz) ! Reflectivity
  REAL, INTENT(OUT) :: refc(nx,ny,nz) ! Composite reflectivity

  INTEGER :: i,j,k
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!

  DO j=1,ny
    DO i=1,nx
      refc(i,j,1)= ref(i,j,1)
      DO k=2,nz-1
        refc(i,j,1) = MAX(refc(i,j,1),ref(i,j,k))
      END DO
    END DO
  END DO

  DO j=1,ny
    DO i=1,nx
      DO k=2,nz-1
        refc(i,j,k) = refc(i,j,1)
      END DO
    END DO
  END DO


  RETURN
END SUBROUTINE cal_rfc

!
!
!##################################################################
!##################################################################
!######                                                      ######
!######                 SUBROUTINES cal_vorp                 ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!  Calculate Vort*10^5 (1/s) value.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Min Zou
!  3/2/98
!
!  MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!


SUBROUTINE cal_vorp(tem9,u,v,x,y,nx,ny,nz) 1

  IMPLICIT NONE
  INTEGER :: nx,ny,nz
  REAL :: tem9(nx,ny,nz)
  REAL :: u(nx,ny,nz), v(nx,ny,nz)
  REAL :: x(nx), y(ny)

  INTEGER :: i,j,k
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!

  DO k=2,nz-2
    DO  j=2,ny-2
      DO  i=2,nx-2
        tem9(i,j,k)= 1.0E5*(                                            &
               (v(i+1,j,k)-v(i-1,j,k)+v(i+1,j+1,k)-v(i-1,j+1,k))/       &
               (4*(x(i+1)-x(i)))-                                       &
               (u(i,j+1,k)-u(i,j-1,k)+u(i+1,j+1,k)-u(i+1,j-1,k))/       &
               (4*(y(j+1)-y(j))) )
      END DO
    END DO
  END DO

  DO j=2,ny-2
    DO i=2,nx-2
      tem9(i,j,   1)=tem9(i,j,   2)
      tem9(i,j,nz-1)=tem9(i,j,nz-2)
    END DO
  END DO

  DO k=1,nz-1
    DO j=2,ny-2
      tem9(   1,j,k)=tem9(   2,j,k)
      tem9(nx-1,j,k)=tem9(nx-2,j,k)
    END DO
  END DO

  DO  k=1,nz-1
    DO i=1,nx-1
      tem9(i,   1,k)=tem9(i,   2,k)
      tem9(i,ny-1,k)=tem9(i,ny-2,k)
    END DO
  END DO

  RETURN
END SUBROUTINE cal_vorp

!
!##################################################################
!##################################################################
!######                                                      ######
!######                  SUBROUTINE cal_div                  ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!  Calculate 1000.*Divergence (1/s)
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Min Zou
!  3/2/98
!
!  MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!


SUBROUTINE cal_div(tem9,u,v,x,y,nx,ny,nz) 1

  IMPLICIT NONE
  INTEGER :: nx,ny,nz
  REAL :: tem9(nx,ny,nz)
  REAL :: u(nx,ny,nz), v(nx,ny,nz)
  REAL :: x(nx), y(ny)

  INTEGER :: i,j,k
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  DO k=2,nz-1
    DO j=2,ny-1
      DO i=2,nx-1
        tem9(i,j,k)=1000.*( (u(i+1,j,k)-u(i,j,k))/(x(i+1)-x(i))         &
                          + (v(i,j+1,k)-v(i,j,k))/(y(j+1)-y(j)) )
      END DO
    END DO
  END DO

  DO j=2,ny-2
    DO i=2,nx-2
      tem9(i,j,   1)=tem9(i,j,   2)
      tem9(i,j,nz-1)=tem9(i,j,nz-2)
    END DO
  END DO

  RETURN
END SUBROUTINE cal_div
!
!##################################################################
!##################################################################
!######                                                      ######
!######                 SUBROUTINES cal_divq                 ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!  Calculate Moist Conv.*1000. (g/kg/s)
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Min Zou
!  3/2/98
!
!  MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!


SUBROUTINE cal_divq(tem9,u,v,qv,x,y,nx,ny,nz) 1

  IMPLICIT NONE
  INTEGER :: nx,ny,nz
  REAL :: tem9(nx,ny,nz)
  REAL :: u(nx,ny,nz), v(nx,ny,nz), qv(nx,ny,nz)
  REAL :: x(nx), y(ny)

  INTEGER :: i,j,k
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  DO k=2,nz-1
    DO j=2,ny-1
      DO i=2,nx-1
        tem9(i,j,k)= -1. * 1000.0 * 1000.0 * 0.5 * (                    &
              ( u(i+1,j,k)*(qv(i,j,k)+qv(i+1,j,k))                      &
              -u(i,j,k)*(qv(i-1,j,k)+qv(i,j,k)) ) /(x(i+1)-x(i))        &
              + ( v(i,j+1,k)*(qv(i,j,k)+qv(i,j+1,k))                    &
              -v(i,j,k)*(qv(i,j-1,k)+qv(i,j,k)) ) /(y(j+1)-y(j)) )
      END DO
    END DO
  END DO

  DO j=2,ny-2
    DO i=2,nx-2
      tem9(i,j,   1)=tem9(i,j,   2)
      tem9(i,j,nz-1)=tem9(i,j,nz-2)
    END DO
  END DO

  RETURN
END SUBROUTINE cal_divq

!
!##################################################################
!##################################################################
!######                                                      ######
!######                 SUBROUTINES cal_vtp                  ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!  Calculate perturbation wind vectors.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Min Zou
!  3/2/98
!
!  MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!


SUBROUTINE cal_vtp(tem7,tem8,tem9,uprt,vprt,wprt,nx,ny,nz,              & 1,3
           vtpunits,label,length)

  IMPLICIT NONE
  INTEGER :: nx,ny,nz
  REAL :: tem7(nx,ny,nz), tem8(nx,ny,nz), tem9(nx,ny,nz)
  REAL :: uprt(nx,ny,nz),vprt(nx,ny,nz),wprt(nx,ny,nz)
  INTEGER :: vtpunits,length
  CHARACTER (LEN=*) :: label

  INTEGER :: i,j,k, onvf
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  onvf = 0
  CALL avgx(uprt , onvf,                                                &
           nx,ny,nz, 1,nx-1, 1,ny-1, 1,nz-1, tem7)
  CALL avgy(vprt , onvf,                                                &
           nx,ny,nz, 1,nx-1, 1,ny-1, 1,nz-1, tem8)
  CALL avgz(wprt , onvf,                                                &
           nx,ny,nz, 1,nx-1, 1,ny-1, 1,nz-1, tem9)
  IF(vtpunits == 1) THEN
    label = '(m/s)'
    length =5
  ELSE IF(vtpunits == 2) THEN
    label = '(kts)'
    length = 5
    DO k=1,nz
      DO j=1,ny
        DO i=1,nx
          tem7(i,j,k) = tem7(i,j,k)*1.943844
          tem8(i,j,k) = tem8(i,j,k)*1.943844
          tem9(i,j,k) = tem9(i,j,k)*1.943844
        END DO
      END DO
    END DO
  ELSE IF(vtpunits == 3) THEN
    label = '(MPH)'
    length = 5
    DO k=1,nz
      DO j=1,ny
        DO i=1,nx
          tem7(i,j,k) = tem7(i,j,k)*2.236936
          tem8(i,j,k) = tem8(i,j,k)*2.236936
          tem9(i,j,k) = tem9(i,j,k)*2.236936
        END DO
      END DO
    END DO
  END IF
  RETURN
END SUBROUTINE cal_vtp
!
!##################################################################
!##################################################################
!######                                                      ######
!######                 SUBROUTINES cal_vtrstrm              ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!  Calculate wind streamline
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Min Zou
!  3/2/98
!
!  MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!


SUBROUTINE cal_vtrstrm(tem7,tem8,tem9,u,v,w,nx,ny,nz,aspratio) 1,3

  IMPLICIT NONE
  INTEGER :: nx,ny,nz
  REAL :: tem7(nx,ny,nz), tem8(nx,ny,nz),tem9(nx,ny,nz)
  REAL :: u(nx,ny,nz),v(nx,ny,nz),w(nx,ny,nz)
  REAL :: aspratio

  INTEGER :: i,j,k,onvf
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  onvf = 0
  CALL avgx(u , onvf,                                                   &
           nx,ny,nz, 1,nx-1, 1,ny-1, 1,nz-1, tem7)
  CALL avgy(v , onvf,                                                   &
           nx,ny,nz, 1,nx-1, 1,ny-1, 1,nz-1, tem8)
  CALL avgz(w , onvf,                                                   &
           nx,ny,nz, 1,nx-1, 1,ny-1, 1,nz-1, tem9)

  DO k=1,nz
    DO j=1,ny
      DO i=1,nx
        tem9(i,j,k)=aspratio*tem9(i,j,k)
      END DO
    END DO
  END DO
  RETURN
END SUBROUTINE cal_vtrstrm
!
!##################################################################
!##################################################################
!######                                                      ######
!######                 SUBROUTINES cal_vtpstrm              ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!  Calculate the perturbation of wind streamlins.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Min Zou
!  3/2/98
!
!  MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!


SUBROUTINE cal_vtpstrm(tem7,tem8,tem9,uprt,vprt,wprt,nx,ny,nz,          & 1,3
           aspratio)

  IMPLICIT NONE
  INTEGER :: nx,ny,nz
  REAL :: tem7(nx,ny,nz), tem8(nx,ny,nz),tem9(nx,ny,nz)
  REAL :: uprt(nx,ny,nz),vprt(nx,ny,nz),wprt(nx,ny,nz)
  REAL :: aspratio

  INTEGER :: i,j,k,onvf
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!

  onvf = 0
  CALL avgx(uprt , onvf,                                                &
           nx,ny,nz, 1,nx-1, 1,ny-1, 1,nz-1, tem7)
  CALL avgy(vprt , onvf,                                                &
           nx,ny,nz, 1,nx-1, 1,ny-1, 1,nz-1, tem8)
  CALL avgz(wprt , onvf,                                                &
           nx,ny,nz, 1,nx-1, 1,ny-1, 1,nz-1, tem9)


  DO k=1,nz
    DO j=1,ny
      DO i=1,nx
        tem9(i,j,k)=aspratio*tem9(i,j,k)
      END DO
    END DO
  END DO

  RETURN
END SUBROUTINE cal_vtpstrm

!
!
!##################################################################
!##################################################################
!######                                                      ######
!######                 SUBROUTINES cal_ vs                  ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!  Calculate Vertical wind shear*1000(1/s)
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Min Zou
!  3/2/98
!
!  MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!


SUBROUTINE cal_vs(tem9,u,v,zp,tem7,tem8,nx,ny,nz) 1,2

  IMPLICIT NONE
  INTEGER :: nx,ny,nz
  REAL :: tem9(nx,ny,nz)
  REAL :: zp(nx,ny,nz),u(nx,ny,nz), v(nx,ny,nz)
  REAL :: tem7(nx,ny,nz), tem8(nx,ny,nz)

  INTEGER :: i,j,k, onvf
  REAL :: tmp1, tmp2, tmp3
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!

  onvf = 0
  CALL avgx(u , onvf,                                                   &
           nx,ny,nz, 1,nx-1, 1,ny-1, 1,nz-1, tem7)
  onvf = 0
  CALL avgy(v , onvf,                                                   &
           nx,ny,nz, 1,nx-1, 1,ny-1, 1,nz-1, tem8)
  DO k=2,nz-2
    DO j=1,ny
      DO i=1,nx
        tmp1 = ( zp(i,j,k+2) + zp(i,j,k+1) -                            &
                       zp(i,j,k) - zp(i,j,k-1) )*0.5
        tmp2 = tem7(i,j,k+1) - tem7(i,j,k-1)
        tmp3 = tem8(i,j,k+1) - tem8(i,j,k-1)
        tem9(i,j,k) = 1000.*SQRT((tmp2/tmp1)**2+(tmp3/tmp1)**2)

      END DO
    END DO
  END DO
  RETURN
END SUBROUTINE cal_vs

!
!
!##################################################################
!##################################################################
!######                                                      ######
!######                 SUBROUTINES cal_gric                 ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!  Calculate Richardson Number
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Min Zou
!  3/2/98
!
!  MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!


SUBROUTINE cal_gric(tem9,u,v,zp,pt,tem7,tem8,nx,ny,nz) 1,2

  IMPLICIT NONE
  INTEGER :: nx,ny,nz
  REAL :: tem9(nx,ny,nz)
  REAL :: zp(nx,ny,nz),pt(nx,ny,nz),u(nx,ny,nz), v(nx,ny,nz)
  REAL :: tem7(nx,ny,nz), tem8(nx,ny,nz)

  INTEGER :: i,j,k, onvf
  REAL :: tmp1, tmp2, tmp3, tmp4
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  onvf = 0
  CALL avgx(u , onvf,                                                   &
           nx,ny,nz, 1,nx-1, 1,ny-1, 1,nz-1, tem7)
  onvf = 0
  CALL avgy(v , onvf,                                                   &
           nx,ny,nz, 1,nx-1, 1,ny-1, 1,nz-1, tem8)
  DO k=2,nz-2
    DO j=1,ny
      DO i=1,nx
        tmp1 = ( zp(i,j,k+2) + zp(i,j,k+1) -                            &
                       zp(i,j,k) - zp(i,j,k-1) )*0.5
        tmp2 = tem7(i,j,k+1) - tem7(i,j,k-1)
        tmp3 = tem8(i,j,k+1) - tem8(i,j,k-1)
        tmp4 =  (tmp2/tmp1)**2 + (tmp3/tmp1)**2
        tmp2 = pt(i,j,k+1)-pt(i,j,k-1)
        tmp3 = 9.8*tmp2/pt(i,j,k)/tmp1/(tmp4+1.e-20)
        tem9(i,j,k) = SIGN( MIN(ABS(tmp3),10.0), tmp3 )
      END DO
    END DO
  END DO
  RETURN
END SUBROUTINE cal_gric

!
!
!##################################################################
!##################################################################
!######                                                      ######
!######                 SUBROUTINES cal_avor                 ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!  Calculate absolute Vort*10^5 (1/s)
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Min Zou
!  3/2/98
!
!  MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!


SUBROUTINE cal_avor(tem9,u,v,x,y,nx,ny,nz,mode,flagsin,omega,           & 1,1
           sinlat,tem1,tem2,tem3)

  IMPLICIT NONE
  INTEGER :: nx,ny,nz
  REAL :: tem9(nx,ny,nz),sinlat(nx,ny)
  REAL :: u(nx,ny,nz), v(nx,ny,nz)
  REAL :: x(nx), y(ny)
  REAL :: tem1(nx,ny,nz), tem2(nx,ny,nz), tem3(nx,ny,nz)
  REAL :: omega
  INTEGER :: mode,flagsin

  INTEGER :: i,j,k
  REAL :: tmp1
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  DO k=2,nz-2
    DO j=2,ny-2
      DO i=2,nx-2
        tem9(i,j,k)= 1.0E5*(                                            &
               (v(i+1,j,k)-v(i-1,j,k)+v(i+1,j+1,k)-v(i-1,j+1,k))/       &
               (4*(x(i+1)-x(i)))-                                       &
               (u(i,j+1,k)-u(i,j-1,k)+u(i+1,j+1,k)-u(i+1,j-1,k))/       &
               (4*(y(j+1)-y(j))) )
      END DO
    END DO
  END DO

  DO j=2,ny-2
    DO i=2,nx-2
      tem9(i,j,   1)=tem9(i,j,   2)
      tem9(i,j,nz-1)=tem9(i,j,nz-2)
    END DO
  END DO

  DO k=1,nz-1
    DO j=2,ny-2
      tem9(   1,j,k)=tem9(   2,j,k)
      tem9(nx-1,j,k)=tem9(nx-2,j,k)
    END DO
  END DO

  DO k=1,nz-1
    DO i=1,nx-1
      tem9(i,   1,k)=tem9(i,   2,k)
      tem9(i,ny-1,k)=tem9(i,ny-2,k)
    END DO
  END DO

  IF(mode == 1.OR.mode == 4.OR.mode == 6.OR.mode == 7) THEN
    IF( flagsin == 0) THEN
      CALL gtsinlat(nx,ny,x,y,sinlat,tem1,tem2, tem3)
      tmp1 = 2.0* omega
      DO j=1,ny
        DO i=1,nx
          sinlat(i,j) = tmp1 * sinlat(i,j)*1.0E5
        END DO
      END DO
      flagsin=1
    END IF
    DO k=1,nz
      DO j=1,ny
        DO i=1,nx
          tem9(i,j,k) = tem9(i,j,k) + sinlat(i,j)
        END DO
      END DO
    END DO
  END IF
  RETURN
END SUBROUTINE cal_avor

!
!##################################################################
!##################################################################
!######                                                      ######
!######                 SUBROUTINES cal_qt                   ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!  Calculate Total water & vapor (g/kg)
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Min Zou
!  3/2/98
!
!  MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!


SUBROUTINE cal_qt(tem9,qc,qr,qi,qs,qh,qv,nx,ny,nz) 1

  IMPLICIT NONE
  INTEGER :: nx,ny,nz
  REAL :: tem9(nx,ny,nz)
  REAL :: qc(nx,ny,nz), qr(nx,ny,nz), qi(nx,ny,nz)
  REAL :: qs(nx,ny,nz), qh(nx,ny,nz), qv(nx,ny,nz)

  INTEGER :: i,j,k
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  DO k=1,nz
    DO j=1,ny
      DO i=1,nx
        tem9(i,j,k)=qc(i,j,k)+qr(i,j,k)+qi(i,j,k)+qs(i,j,k)+            &
                      qh(i,j,k) + qv(i,j,k)
      END DO
    END DO
  END DO
  RETURN
END SUBROUTINE cal_qt


!
!
!##################################################################
!##################################################################
!######                                                      ######
!######                 SUBROUTINES cal_rhi                  ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!  Calculate rhi value.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Min Zou
!  3/2/98
!
!  MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!


SUBROUTINE cal_rhi(tem9,pt, pprt ,pbar,qv,tem1,tem2,nx,ny,nz) 1,2

  IMPLICIT NONE
  INTEGER :: nx,ny,nz
  REAL :: tem9(nx,ny,nz), tem1(nx,ny,nz), tem2(nx,ny,nz)
  REAL :: pt(nx,ny,nz), pprt(nx,ny,nz), pbar(nx,ny,nz), qv(nx,ny,nz)

  INTEGER :: i,j,k
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!

  CALL temper (nx,ny,nz,pt, pprt ,pbar,tem1)
  CALL getqvs(nx,ny,nz, 1,nx-1,1,ny-1,1,nz-1, pbar,tem1,tem2)

  DO k=1,nz-1
    DO j=1,ny-1
      DO i=1,nx-1
        tem9(i,j,k) = qv(i,j,k)/tem2(i,j,k)
      END DO
    END DO
  END DO

  RETURN
END SUBROUTINE cal_rhi

!
!##################################################################
!##################################################################
!######                                                      ######
!######                 SUBROUTINES cal_xuv                  ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!  Calculate wind vector (xuv)
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Min Zou
!  3/2/98
!
!  MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!


SUBROUTINE cal_xuv(tem7,tem8,tem9,u,v,w,nx,ny,nz,xuvunits,              & 1,3
           label,length)

  IMPLICIT NONE
  INTEGER :: nx,ny,nz
  REAL :: tem7(nx,ny,nz), tem8(nx,ny,nz),tem9(nx,ny,nz)
  REAL :: u(nx,ny,nz),v(nx,ny,nz),w(nx,ny,nz)
  INTEGER :: xuvunits,length
  CHARACTER (LEN=*) :: label

  INTEGER :: i,j,k,onvf
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  onvf = 0
  CALL avgx(u , onvf,                                                   &
           nx,ny,nz, 1,nx-1, 1,ny-1, 1,nz-1, tem7)
  CALL avgy(v , onvf,                                                   &
           nx,ny,nz, 1,nx-1, 1,ny-1, 1,nz-1, tem8)
  CALL avgz(w , onvf,                                                   &
           nx,ny,nz, 1,nx-1, 1,ny-1, 1,nz-1, tem9)

  IF(xuvunits == 1) THEN
    label = '(m/s)'
    length =5
  ELSE IF(xuvunits == 2) THEN
    label = '(kts)'
    length = 5
    DO k=1,nz
      DO j=1,ny
        DO i=1,nx
          tem7(i,j,k) = tem7(i,j,k)*1.943844
          tem8(i,j,k) = tem8(i,j,k)*1.943844
          tem9(i,j,k) = tem9(i,j,k)*1.943844
        END DO
      END DO
    END DO
  ELSE IF(xuvunits == 3) THEN
    label = '(MPH)'
    length = 5
    DO k=1,nz
      DO j=1,ny
        DO i=1,nx
          tem7(i,j,k) = tem7(i,j,k)*2.236936
          tem8(i,j,k) = tem8(i,j,k)*2.236936
          tem9(i,j,k) = tem9(i,j,k)*2.236936
        END DO
      END DO
    END DO
  END IF
  RETURN
END SUBROUTINE cal_xuv

!
!##################################################################
!##################################################################
!######                                                      ######
!######                 SUBROUTINES cal_vtr                  ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!  Calculate wind vector
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Min Zou
!  3/2/98
!
!  MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!


SUBROUTINE cal_vtr(tem7,tem8,tem9,u,v,w,nx,ny,nz,vtrunits,label,        & 2,3
           length)

  IMPLICIT NONE
  INTEGER :: nx,ny,nz
  REAL :: tem7(nx,ny,nz), tem8(nx,ny,nz),tem9(nx,ny,nz)
  REAL :: u(nx,ny,nz),v(nx,ny,nz),w(nx,ny,nz)
  INTEGER :: vtrunits,length
  CHARACTER (LEN=*) :: label

  INTEGER :: i,j,k,onvf
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  onvf = 0
  CALL avgx(u , onvf,                                                   &
           nx,ny,nz, 1,nx-1, 1,ny-1, 1,nz-1, tem7)
  CALL avgy(v , onvf,                                                   &
           nx,ny,nz, 1,nx-1, 1,ny-1, 1,nz-1, tem8)
  CALL avgz(w , onvf,                                                   &
           nx,ny,nz, 1,nx-1, 1,ny-1, 1,nz-1, tem9)

  IF(vtrunits == 1) THEN
    label = '(m/s)'
    length =5
  ELSE IF(vtrunits == 2) THEN
    label = '(kts)'
    length = 5
    DO k=1,nz
      DO j=1,ny
        DO i=1,nx
          tem7(i,j,k) = tem7(i,j,k)*1.943844
          tem8(i,j,k) = tem8(i,j,k)*1.943844
          tem9(i,j,k) = tem9(i,j,k)*1.943844
        END DO
      END DO
    END DO
  ELSE IF(vtrunits == 3) THEN
    label = '(MPH)'
    length = 5
    DO k=1,nz
      DO j=1,ny
        DO i=1,nx
          tem7(i,j,k) = tem7(i,j,k)*2.236936
          tem8(i,j,k) = tem8(i,j,k)*2.236936
          tem9(i,j,k) = tem9(i,j,k)*2.236936
        END DO
      END DO
    END DO
  END IF
  RETURN
END SUBROUTINE cal_vtr
!
!##################################################################
!##################################################################
!######                                                      ######
!######                 SUBROUTINES cal_vtrobs               ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!  Calculate wind observation value.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Min Zou
!  3/2/98
!
!  MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!


SUBROUTINE cal_vtrobs(dd,ff,drot, vtrunits) 1,1

  IMPLICIT NONE
  INCLUDE 'arpsplt.inc'

  INTEGER :: ovrobs,obsset,obscol,obs_marktyp
  REAL :: obs_marksz
  COMMON /obspar/ ovrobs,obsset,obscol,obs_marktyp, obs_marksz

  INTEGER :: nobs
  COMMON /sfc_obs1/ nobs
  REAL :: latob(mxsfcob),lonob(mxsfcob)
  REAL :: obs1(mxsfcob),obs2(mxsfcob)
  COMMON /sfc_obs2/ latob,lonob,obs1,obs2

  REAL :: dd(mxsfcob),ff(mxsfcob)
  REAL :: drot

  INTEGER :: iob, vtrunits
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  DO iob=1,nobs
    IF(dd(iob) >= 0. .AND. dd(iob) < 360. .AND.                         &
          ff(iob) >= 0. .AND. ff(iob) < 98.) THEN
      CALL ddrotuv(1,lonob(iob),dd(iob),ff(iob),                        &
                    drot,obs1(iob),obs2(iob))
      obs1(iob)=0.51444*obs1(iob)
      obs2(iob)=0.51444*obs2(iob)
    ELSE
      obs1(iob)=-999.
      obs2(iob)=-999.
    END IF
  END DO
  obsset=1

  IF(vtrunits == 2) THEN    !! kts
    DO iob=1,nobs
      IF(obs1(iob) /= -999.) obs1(iob)= obs1(iob)*1.943844
      IF(obs2(iob) /= -999.) obs2(iob)= obs2(iob)*1.943844
    END DO
  ELSE IF(vtrunits == 3) THEN   !! MPH
    DO iob=1,nobs
      IF(obs1(iob) /= -999.) obs1(iob) = obs1(iob)*2.236936
      IF(obs1(iob) /= -999.) obs2(iob) = obs2(iob)*2.236936
    END DO
  END IF
  RETURN
END SUBROUTINE cal_vtrobs

!
!##################################################################
!##################################################################
!######                                                      ######
!######                 SUBROUTINES cal_viqc                 ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!  Calculate viqc
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Min Zou
!  3/2/98
!
!  MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!


SUBROUTINE cal_viqc(tem7, qc, rhobar, zp, nx,ny,nz) 1

  IMPLICIT NONE
  INTEGER :: nx,ny,nz
  REAL :: tem7(nx,ny,nz)
  REAL :: qc(nx,ny,nz), rhobar(nx,ny,nz), zp(nx,ny,nz)

  INTEGER :: i,j,k
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  DO j=1,ny
    DO i=1,nx
      tem7(i,j,1)=0.
    END DO
  END DO

  DO k=2,nz-2
    DO j=1,ny
      DO i=1,nx
        tem7(i,j,1) = tem7(i,j,1)+qc(i,j,k)*rhobar(i,j,k)               &
                        *(zp(i,j,k+1)-zp(i,j,k))
      END DO
    END DO
  END DO
  RETURN
END SUBROUTINE cal_viqc


!
!##################################################################
!##################################################################
!######                                                      ######
!######                 SUBROUTINES cal_viqr                 ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!  Calculate viqr
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Min Zou
!  3/2/98
!
!  MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!


SUBROUTINE cal_viqr(tem7, qr, rhobar, zp, nx,ny,nz) 1

  IMPLICIT NONE
  INTEGER :: nx,ny,nz
  REAL :: tem7(nx,ny,nz)
  REAL :: qr(nx,ny,nz), rhobar(nx,ny,nz), zp(nx,ny,nz)

  INTEGER :: i,j,k
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  DO j=1,ny
    DO i=1,nx
      tem7(i,j,1)=0.
    END DO
  END DO

  DO k=2,nz-2
    DO j=1,ny
      DO i=1,nx
        tem7(i,j,1) = tem7(i,j,1)+qr(i,j,k)*rhobar(i,j,k)               &
                        *(zp(i,j,k+1)-zp(i,j,k))
      END DO
    END DO
  END DO
  RETURN
END SUBROUTINE cal_viqr


!
!##################################################################
!##################################################################
!######                                                      ######
!######                 SUBROUTINES cal_viqi                 ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!  Calculate viqi
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Min Zou
!  3/2/98
!
!  MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!


SUBROUTINE cal_viqi(tem7, qi, rhobar, zp, nx,ny,nz) 1

  IMPLICIT NONE
  INTEGER :: nx,ny,nz
  REAL :: tem7(nx,ny,nz)
  REAL :: qi(nx,ny,nz), rhobar(nx,ny,nz), zp(nx,ny,nz)

  INTEGER :: i,j,k
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  DO j=1,ny
    DO i=1,nx
      tem7(i,j,1)=0.
    END DO
  END DO

  DO k=2,nz-2
    DO j=1,ny
      DO i=1,nx
        tem7(i,j,1) = tem7(i,j,1)+qi(i,j,k)*rhobar(i,j,k)               &
                        *(zp(i,j,k+1)-zp(i,j,k))
      END DO
    END DO
  END DO
  RETURN
END SUBROUTINE cal_viqi

!
!##################################################################
!##################################################################
!######                                                      ######
!######                 SUBROUTINES cal_viqs                 ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!  Calculate viqs
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Min Zou
!  3/2/98
!
!  MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!


SUBROUTINE cal_viqs(tem7, qs, rhobar, zp, nx,ny,nz) 1

  IMPLICIT NONE
  INTEGER :: nx,ny,nz
  REAL :: tem7(nx,ny,nz)
  REAL :: qs(nx,ny,nz), rhobar(nx,ny,nz), zp(nx,ny,nz)

  INTEGER :: i,j,k
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  DO j=1,ny
    DO i=1,nx
      tem7(i,j,1)=0.
    END DO
  END DO

  DO k=2,nz-2
    DO j=1,ny
      DO i=1,nx
        tem7(i,j,1) = tem7(i,j,1)+qs(i,j,k)*rhobar(i,j,k)               &
                        *(zp(i,j,k+1)-zp(i,j,k))
      END DO
    END DO
  END DO
  RETURN
END SUBROUTINE cal_viqs

!
!##################################################################
!##################################################################
!######                                                      ######
!######                 SUBROUTINES cal_viqh                 ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!  Calculate viqh
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Min Zou
!  3/2/98
!
!  MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!


SUBROUTINE cal_viqh(tem7, qh, rhobar, zp, nx,ny,nz) 1

  IMPLICIT NONE
  INTEGER :: nx,ny,nz
  REAL :: tem7(nx,ny,nz)
  REAL :: qh(nx,ny,nz), rhobar(nx,ny,nz), zp(nx,ny,nz)

  INTEGER :: i,j,k
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  DO j=1,ny
    DO i=1,nx
      tem7(i,j,1)=0.
    END DO
  END DO

  DO k=2,nz-2
    DO j=1,ny
      DO i=1,nx
        tem7(i,j,1) = tem7(i,j,1)+qh(i,j,k)*rhobar(i,j,k)               &
                        *(zp(i,j,k+1)-zp(i,j,k))
      END DO
    END DO
  END DO
  RETURN
END SUBROUTINE cal_viqh

!##################################################################
!##################################################################
!######                                                      ######
!######                 SUBROUTINES cal_vil                  ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!  Calculate Vert. Integ Liquid (kg/m2)
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Min Zou
!  3/2/98
!
!  MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!


SUBROUTINE cal_vil(tem7,qc,qr,rhobar,zp, nx,ny,nz,tem6) 1

  IMPLICIT NONE
  INTEGER :: nx,ny,nz
  REAL :: tem7(nx,ny,nz)
  REAL :: qc(nx,ny,nz),qr(nx,ny,nz), rhobar(nx,ny,nz), zp(nx,ny,nz)
  REAL :: tem6(nx,ny,nz)

  INTEGER :: i,j,k
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  DO j=1,ny
    DO i=1,nx
      tem7(i,j,1)=0.
    END DO
  END DO

  DO k=2,nz-2
    DO j=1,ny
      DO i=1,nx
        tem6(i,j,k) = qc(i,j,k) + qr(i,j,k)
        tem7(i,j,1) = tem7(i,j,1)+tem6(i,j,k)*rhobar(i,j,k)             &
                  *(zp(i,j,k+1)-zp(i,j,k))
      END DO
    END DO
  END DO

  RETURN
END SUBROUTINE cal_vil


!##################################################################
!##################################################################
!######                                                      ######
!######                 SUBROUTINES cal_vii                  ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!  Calculate Vert. Integrated ice (kg/m2)
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Min Zou
!  3/2/98
!
!  MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!


SUBROUTINE cal_vii(tem7,qi,qs,qh,rhobar,zp, nx,ny,nz,tem6) 1

  IMPLICIT NONE
  INTEGER :: nx,ny,nz
  REAL :: tem7(nx,ny,nz)
  REAL :: qi(nx,ny,nz), qs(nx,ny,nz), qh(nx,ny,nz)
  REAL :: rhobar(nx,ny,nz), zp(nx,ny,nz)
  REAL :: tem6(nx,ny,nz)

  INTEGER :: i,j,k
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!

  DO j=1,ny
    DO i=1,nx
      tem7(i,j,1)=0.
    END DO
  END DO

  DO k=2,nz-2
    DO j=1,ny
      DO i=1,nx
        tem6(i,j,k) = qi(i,j,k) + qs(i,j,k) + qh(i,j,k)
        tem7(i,j,1) = tem7(i,j,1)+tem6(i,j,k)*rhobar(i,j,k)             &
                  *(zp(i,j,k+1)-zp(i,j,k))
      END DO
    END DO
  END DO
  RETURN
END SUBROUTINE cal_vii

!##################################################################
!##################################################################
!######                                                      ######
!######                 SUBROUTINES cal_vic                  ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!  Calculate vic
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Min Zou
!  3/2/98
!
!  MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!


SUBROUTINE cal_vic(tem7,qc,qr,qi,qs,qh,rhobar,zp,nx,ny,nz,tem6) 1

  IMPLICIT NONE
  INTEGER :: nx,ny,nz
  REAL :: tem7(nx,ny,nz)
  REAL :: qc(nx,ny,nz), qr(nx,ny,nz), qi(nx,ny,nz), qs(nx,ny,nz)
  REAL :: qh(nx,ny,nz), rhobar(nx,ny,nz), zp(nx,ny,nz)
  REAL :: tem6(nx,ny,nz)

  INTEGER :: i,j,k
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  DO j=1,ny
    DO i=1,nx
      tem7(i,j,1)=0.
    END DO
  END DO

  DO k=2,nz-2
    DO j=1,ny
      DO i=1,nx
        tem6(i,j,k) = qc(i,j,k) + qr(i,j,k) + qi(i,j,k) +               &
                      qs(i,j,k) + qh(i,j,k)
        tem7(i,j,1) = tem7(i,j,1)+tem6(i,j,k)*rhobar(i,j,k)             &
                  *(zp(i,j,k+1)-zp(i,j,k))
      END DO
    END DO
  END DO

  RETURN
END SUBROUTINE cal_vic

!##################################################################
!##################################################################
!######                                                      ######
!######                 SUBROUTINES cal_vit                  ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!  Calculate vit ( vertically intergrated total water(kg/m**2))
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Min Zou
!  5/10/98
!
!  MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!


SUBROUTINE cal_vit(tem7,qv,qc,qr,qi,qs,qh,rhobar,zp,nx,ny,nz,tem6) 1

  IMPLICIT NONE
  INTEGER :: nx,ny,nz
  REAL :: tem7(nx,ny,nz)
  REAL :: qc(nx,ny,nz), qr(nx,ny,nz), qi(nx,ny,nz), qs(nx,ny,nz)
  REAL :: qh(nx,ny,nz), qv(nx,ny,nz), rhobar(nx,ny,nz), zp(nx,ny,nz)
  REAL :: tem6(nx,ny,nz)

  INTEGER :: i,j,k
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  DO j=1,ny
    DO i=1,nx
      tem7(i,j,1)=0.
    END DO
  END DO

  DO k=2,nz-2
    DO j=1,ny
      DO i=1,nx
        tem6(i,j,k) = qc(i,j,k) + qr(i,j,k) + qi(i,j,k) +               &
                      qs(i,j,k) + qh(i,j,k) + qv(i,j,k)
        tem7(i,j,1) = tem7(i,j,1)+tem6(i,j,k)*rhobar(i,j,k)             &
                  *(zp(i,j,k+1)-zp(i,j,k))
      END DO
    END DO
  END DO

  RETURN
END SUBROUTINE cal_vit
!
!##################################################################
!##################################################################
!######                                                      ######
!######                 SUBROUTINES cal_pw                   ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE cal_pw(tem7,qv,rhobar,zp,nx,ny,nz,tem6) 1
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!  Calculate pw ( precipitable water vapor(cm)
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Min Zou
!  5/10/98
!
!  MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
  INTEGER :: nx,ny,nz
  REAL :: tem7(nx,ny,nz)
  REAL :: qv(nx,ny,nz), rhobar(nx,ny,nz), zp(nx,ny,nz)
  REAL :: tem6(nx,ny,nz)

  INTEGER :: i,j,k
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  DO j=1,ny
    DO i=1,nx
      tem7(i,j,1)=0.
    END DO
  END DO

  DO k=2,nz-2
    DO j=1,ny
      DO i=1,nx
        tem6(i,j,k) = qv(i,j,k)
        tem7(i,j,1) = tem7(i,j,1)                                       &
                    + tem6(i,j,k)*rhobar(i,j,k)*(zp(i,j,k+1)-zp(i,j,k))
      END DO
    END DO
  END DO

! change kg/m**2 to cm

  DO j=1,ny
    DO i=1,nx
      tem7(i,j,1) = 0.1*tem7(i,j,1)
    END DO
  END DO

  RETURN
END SUBROUTINE cal_pw
!
!
!##################################################################
!##################################################################
!######                                                      ######
!######                 SUBROUTINES cal_tpr                  ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!  Calculate tpr ( total precipatation rate )
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Min Zou
!  5/15/98
!
!  MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!


SUBROUTINE cal_tpr(tem7,prcrate,nx,ny,nz,tprunits,label,length) 1

  IMPLICIT NONE
  INTEGER :: nx,ny,nz
  REAL :: tem7(nx,ny,nz)
  REAL :: prcrate(nx,ny)
  INTEGER :: tprunits, length
  CHARACTER (LEN=*) :: label

  INTEGER :: i,j
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  DO j=1,ny
    DO i=1,nx
      tem7(i,j,1)=0.
    END DO
  END DO

  IF(tprunits == 1) THEN
    label = 'Total precip. rate(mm/h)'
    length =24
    DO j=1,ny
      DO i=1,nx
           !RLC 1998/05/18
!tem7(i,j,1) = prcrate(i,j)/100./3600. !!kg/(m**2*s) -> mm/h
        tem7(i,j,1) = prcrate(i,j)*3600. !!kg/(m**2*s) -> mm/h
      END DO
    END DO
  ELSE IF(tprunits == 2) THEN
                          !! kg/(m**2*s) -> in/h
    label = 'Total precip. rate(in/h)'
    length =24
    DO j=1,ny
      DO i=1,nx
        tem7(i,j,1) = prcrate(i,j)*0.039370079*3600.
      END DO
    END DO
  END IF

  RETURN
END SUBROUTINE cal_tpr
!
!##################################################################
!##################################################################
!######                                                      ######
!######                 SUBROUTINES cal_gpr                  ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!  Calculate gpr ( Grid-scale precip. rate)
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Min Zou
!  5/15/98
!
!  MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!


SUBROUTINE cal_gpr(tem7,prcrate,prcrate1,nx,ny,nz,                      & 1
           gprunits,label,length)

  IMPLICIT NONE
  INTEGER :: nx,ny,nz
  REAL :: tem7(nx,ny,nz)
  REAL :: prcrate(nx,ny), prcrate1(nx,ny)
  INTEGER :: gprunits, length
  CHARACTER (LEN=*) :: label

  INTEGER :: i,j
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  DO j=1,ny
    DO i=1,nx
      tem7(i,j,1)=0.
    END DO
  END DO

  IF(gprunits == 1) THEN
    label = 'Grid-scale precip. rate(mm/h)'
    length =29
    DO j=1,ny
      DO i=1,nx
        tem7(i,j,1) = (prcrate(i,j)+prcrate1(i,j))*3600. !kg/(m**2*s) -> mm/h
      END DO
    END DO
  ELSE IF(gprunits == 2) THEN
                          !! kg/(m**2*s) -> in/h
    label = 'Grid-scale precip. rate(in/h)'
    length =29
    DO j=1,ny
      DO i=1,nx
        tem7(i,j,1) = (prcrate(i,j)+prcrate1(i,j)) *0.039370079*3600.
      END DO
    END DO
  END IF

  RETURN
END SUBROUTINE cal_gpr
!
!##################################################################
!##################################################################
!######                                                      ######
!######                 SUBROUTINES cal_cpr                  ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!  Calculate cpr ( Convective precip. rate )
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Min Zou
!  5/15/98
!
!  MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!


SUBROUTINE cal_cpr(tem7,prcrate,nx,ny,nz,cprunits,label,length) 1

  IMPLICIT NONE
  INTEGER :: nx,ny,nz
  REAL :: tem7(nx,ny,nz)
  REAL :: prcrate(nx,ny)
  INTEGER :: cprunits, length
  CHARACTER (LEN=*) :: label

  INTEGER :: i,j
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  DO j=1,ny
    DO i=1,nx
      tem7(i,j,1)=0.
    END DO
  END DO

  IF(cprunits == 1) THEN
    label = 'Convective precip. rate(mm/h)'
    length =29
    DO j=1,ny
      DO i=1,nx
        tem7(i,j,1) = prcrate(i,j)*3600. !!kg/(m**2*s) -> mm/h
      END DO
    END DO
  ELSE IF(cprunits == 2) THEN
                          !! kg/(m**2*s) -> in/h
    label = 'Convective precip. rate(in/h)'
    length =29
    DO j=1,ny
      DO i=1,nx
        tem7(i,j,1) = prcrate(i,j)*0.039370079*3600.
      END DO
    END DO
  END IF

  RETURN
END SUBROUTINE cal_cpr
!
!##################################################################
!##################################################################
!######                                                      ######
!######                 SUBROUTINES cal_strm                 ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!  Calculate strom motion
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Min Zou
!  3/2/98
!
!  MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!


SUBROUTINE cal_strm(tem7, tem8,ustrm,vstrm,strmunits,nx,ny,nz,          & 1
           label,length)

  IMPLICIT NONE
  INTEGER :: nx,ny,nz
  REAL :: tem7(nx,ny,nz), tem8(nx,ny,nz)
  REAL :: ustrm(nx,ny) ,vstrm(nx,ny)
  INTEGER :: strmunits, length
  CHARACTER (LEN=*) :: label

  INTEGER :: i,j
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  DO j=1,ny
    DO i=1,nx-1
      tem7(i,j,1)=(ustrm(i,j)+ustrm(i+1,j))*0.5
    END DO
  END DO
  DO j=1,ny-1
    DO i=1,nx
      tem8(i,j,1)=(vstrm(i,j)+vstrm(i,j+1))*0.5
    END DO
  END DO

  IF(strmunits == 1) THEN
    label = 'strom motion (m/s)'
    length =18
  ELSE IF(strmunits == 2) THEN
    label = 'storm motion (kts)'
    length =18
    DO j=1,ny
      DO i=1,nx
        tem7(i,j,1) = tem7(i,j,1)*1.943844
        tem8(i,j,1) = tem8(i,j,1)*1.943844
      END DO
    END DO
  ELSE IF(strmunits == 3) THEN
    label = 'storm motion (MPH)'
    length = 18
    DO j=1,ny
      DO i=1,nx
        tem7(i,j,1) = tem7(i,j,1)*2.236936
        tem8(i,j,1) = tem8(i,j,1)*2.236936
      END DO
    END DO
  END IF

  RETURN
END SUBROUTINE cal_strm

!##################################################################
!##################################################################
!######                                                      ######
!######                 SUBROUTINES cal_dist                 ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Min Zou
!  3/2/98
!
!  MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!


SUBROUTINE cal_dist(haxisu,dx,dy,x01,y01,x02,y02,slicopt,                  & 3
           tmpx,tmpy,distc)


  IMPLICIT NONE
  INTEGER :: haxisu, slicopt
  REAL :: dx,dy
  REAL :: x101, y101, x102,y102
  COMMON /slicev1/x101, y101, x102,y102

  REAL :: x01,y01                  ! the first  point of interpolation
  REAL :: x02,y02                  ! the second point of interpolation
  REAL :: tmpx, tmpy
  CHARACTER (LEN=*) :: distc
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!

  IF ( slicopt == 2 .OR. slicopt == 3 .OR. slicopt == 5 .OR.   & 
       slicopt == 10 .OR. slicopt == 11) THEN
    IF (haxisu == 0) THEN
      tmpx = dx          !!!km
      tmpy = dy         !!!km
      x101 = x01
      y101 = y01
      x102 = x02
      y102 = y02
      WRITE(distc,'('' KM'')')
    ELSE IF(haxisu == 1 ) THEN
      tmpx = dx*0.62137 !!! mile
      tmpy = dy*0.62137 !!! mile
      x101 = x01*0.62137
      y101 = y01*0.62137
      x102 = x02*0.62137
      y102 = y02*0.62137
      WRITE(distc,'('' MILE'')')
    ELSE IF(haxisu == 2) THEN
      tmpx = dx*0.53997 !!!naut mile
      tmpy = dy*0.53997 !!!naut mile
      x101 = x01*0.53997
      y101 = y01*0.53997
      x102 = x02*0.53997
      y102 = y02*0.53997
      WRITE(distc,'('' NAUT MILE'')')
    ELSE IF(haxisu == 3) THEN
      tmpx = dx*3.28084 !!!kft
      tmpy = dy*3.28084 !!!kft
      x101 = x01*3.28084
      y101 = y01*3.28084
      x102 = x02*3.28084
      y102 = y02*3.28084
      WRITE(distc,'('' KFT'')')
    END IF
  END IF
!  IF(mode.eq.5) THEN
!    length=len(distc)
!    CALL strmin(distc,length)
!    write(distc,'(a)') '('//distc(2:length)
!  ENDIF
  RETURN
END SUBROUTINE cal_dist

!##################################################################
!##################################################################
!######                                                      ######
!######                 SUBROUTINES setcords                 ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE setcords(xl,xr,yb,yt,dx,dy, slicopt,                            & 1
           x1,x2,y1,y2,xlabel,ylabel,xstep,ystep,xmstep,ymstep)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!  Set coordinate related variables according to desired length unit
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Min Zou
!  3/2/98
!
!  MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
!
!  INPUT:
!  xl,xr,yb,yt,dx,dy, slicopt
!  OUTPUT:
!  x1,x2,y1,y2,xlabel,ylabel,xstep,ystep,xmstep,ymstep
!
!-----------------------------------------------------------------------

  IMPLICIT NONE
  REAL :: x1,x2,y1,y2,xl,xr,yb,yt, dx,dy
  REAL :: xstep,ystep, xmstep, ymstep
  CHARACTER (LEN=*) :: ylabel
  CHARACTER (LEN=*) :: xlabel
  INTEGER :: slicopt

  INTEGER :: xfont   ! the font of character
  INTEGER :: haxisu, vaxisu
  INTEGER :: lbaxis
  INTEGER :: tickopt
  INTEGER :: xfmat
  REAL :: hmintick,vmajtick,vmintick,hmajtick
  COMMON /var_par/ xfont,haxisu,vaxisu,lbaxis,tickopt,hmintick,         &
          vmajtick, vmintick,hmajtick,xfmat
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!

  IF(slicopt == 1 .OR. slicopt == 4 .OR. slicopt == 6 .OR.  & 
     slicopt == 7 .OR. slicopt ==9) THEN

    IF(haxisu == 0) THEN
      x1=xl
      x2=xr
      y1=yb
      y2=yt
      WRITE(xlabel,'(a15)')'(km)'
      WRITE(ylabel,'(a15)')'(km)'
    ELSE IF(haxisu == 1 ) THEN
      x1=xl*0.62137
      x2=xr*0.62137
      y1=yb*0.62137
      y2=yt*0.62137
      WRITE(xlabel,'(a15)')'(mile)'
      WRITE(ylabel,'(a15)')'(mile)'
    ELSE IF(haxisu == 2) THEN
      x1=xl*0.53997
      x2=xr*0.53997
      y1=yb*0.53997
      y2=yt*0.53997
      WRITE(xlabel,'(a15)')'(naut mile)'
      WRITE(ylabel,'(a15)')'(naut mile)'
    ELSE IF(haxisu == 3) THEN
      x1=xl*3.28084
      x2=xr*3.28084
      y1=yb*3.28084
      y2=yt*3.28084
      WRITE(xlabel,'(a15)')'(kft)'
      WRITE(ylabel,'(a15)')'(kft)'
    END IF
    IF(tickopt == 0) THEN
      xstep = dx
      ystep = dy
      xmstep = 0.
      ymstep = 0.
    ELSE IF(tickopt == 1) THEN
      xstep = hmintick
      ystep = hmintick
      xmstep = hmajtick
      ymstep = hmajtick
    END IF

  END IF

  IF(slicopt == 2 .OR. slicopt == 3 .OR. slicopt == 5) THEN
    IF(haxisu == 0) THEN
      x1=xl
      x2=xr
      WRITE(xlabel,'(a15)')'(km)'
    ELSE IF(haxisu == 1 ) THEN
      x1=xl*0.62137
      x2=xr*0.62137
      WRITE(xlabel,'(a15)')'(mile)'
    ELSE IF(haxisu == 2) THEN
      x1=xl*0.53997
      x2=xr*0.53997
      WRITE(xlabel,'(a15)')'(naut mile)'
    ELSE IF(haxisu == 3) THEN
      x1=xl*3.28084
      x2=xr*3.28084
      WRITE(xlabel,'(a15)')'(kft)'
    END IF
    IF(vaxisu == 0) THEN
      y1=yb
      y2=yt
      WRITE(ylabel,'(a15)')'(km)'
    ELSE IF(vaxisu == 1 ) THEN
      y1=yb*0.62137
      y2=yt*0.62137
      WRITE(ylabel,'(a15)')'(mile)'
    ELSE IF(vaxisu == 2) THEN
      y1=yb*0.53997
      y2=yt*0.53997
      WRITE(ylabel,'(a15)')'(naut mile)'
    ELSE IF(vaxisu == 3) THEN
      y1=yb*3.28084
      y2=yt*3.28084
      WRITE(ylabel,'(a15)')'(kft)'
    ELSE IF(vaxisu == 4) THEN
      y1=yb*3.28084
      y2=yt*3.28084
      WRITE(ylabel,'(a15)')'(presure)'
    END IF
    IF(tickopt == 0) THEN
      xstep = dx
      ystep = dy
      xmstep = 0.
      ymstep = 0.
    ELSE IF(tickopt == 1) THEN
      xstep = hmintick
      ystep = vmintick
      xmstep = hmajtick
      ymstep = vmajtick
    END IF
  END IF

  IF(slicopt == 10 .OR. slicopt == 11) THEN
    IF(haxisu == 0) THEN
      x1=xl
      x2=xr
      WRITE(xlabel,'(a15)')'(km)'
    ELSE IF(haxisu == 1 ) THEN
      x1=xl*0.62137
      x2=xr*0.62137
      WRITE(xlabel,'(a15)')'(mile)'
    ELSE IF(haxisu == 2) THEN
      x1=xl*0.53997
      x2=xr*0.53997
      WRITE(xlabel,'(a15)')'(naut mile)'
    ELSE IF(haxisu == 3) THEN
      x1=xl*3.28084
      x2=xr*3.28084
      WRITE(xlabel,'(a15)')'(kft)'
    END IF 
    IF(vaxisu == 0) THEN
      y1=yb
      y2=yt
      WRITE(ylabel,'(a15)')'(cm)'
    END IF 
    IF(tickopt == 0) THEN
      xstep = dx
      ystep = dy
      xmstep = 0.
      ymstep = 0.
    ELSE IF(tickopt == 1) THEN
      xstep = hmintick
      ystep = vmintick
      xmstep = hmajtick
      ymstep = vmajtick
    END IF
  END IF

  RETURN
END SUBROUTINE setcords



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


SUBROUTINE get_contour(ncont, tcont) 1,1
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!    Set-up ccontour values for when plot option set to 11 .
!    right now only work for several variables:
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Min Zou
!
!  MODIFICATION HISTORY:
!  2/18/97
!
!-----------------------------------------------------------------------
!
!  OUTPUT:
!
!  tcon     the total number of contours.
!  ncon     the contour array.
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  INCLUDE 'arpsplt.inc'

  CHARACTER (LEN=12) :: varname
  COMMON /varplt1/ varname

  INTEGER :: setcontopt, setcontnum
  CHARACTER (LEN=12) :: setcontvar(maxuneva)
  REAL :: setconts(maxunevm,maxuneva)
  COMMON /setcont_var/setcontvar
  COMMON /setcon_par/setcontopt,setcontnum,setconts

  INTEGER :: i, LEN, var, ncont
  REAL :: tcont(maxunevm)

  ncont = 0
  var = 0


  LEN=12
  CALL strlnth( varname, LEN)
  IF(setcontopt == 1) THEN
    DO i = 1,setcontnum
      IF(varname(1:LEN) == setcontvar(i)(1:LEN)) THEN
        var = i
        GO TO 10
      END IF
    END DO
  END IF

  IF(var == 0) RETURN

  10    CONTINUE

  DO i=1,maxunevm
    IF( setconts(i,var) == -9999.0) THEN
      GO TO 100
    END IF
  END DO
  100   ncont = i-1


  DO i=1,ncont
    tcont(i) = setconts(i,var)
  END DO

  IF(var /= 0) PRINT*, setcontvar(var)(1:LEN), ncont, (tcont(i),i=1,ncont)

  RETURN
END SUBROUTINE get_contour