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


SUBROUTINE ctr3d(b,x,y,z, x1,x2,dx,y1,y2,dy,z1,z2,dz,                   & 41,43
           nx,ibgn,iend, ny,jbgn,jend, nz,kbgn,kend,                    &
           label,time,slicopt, kslice, jslice, islice,                  &
           n,xp,yp,axy2d,av2d,zp, runname, factor,tem1,tem2,tem3,       &
           tem4,bb,tem5,hterain,pltopt)

!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!    Set-up 2-d slices of a 3-d data array to contour with
!      subroutine ctr2d.
!
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Ming Xue
!
!  MODIFICATION HISTORY:
!    6/08/92  Added full documentation (K. Brewster)
!
!  12/25/1992 M. Xue and H. Jin
!    Added capability to plot arbitary cross sections.
!
!  8/28/1994 M. Zou
!    Added color shader to contour plot,add full documentation
!
!  3/25/96 (K. Brewster)
!    Added variables isize,jsize,ksize
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!    b        3-dimension array of variable
!    x        x-coord of scalar point (km)
!    y        y-coord of scalar point (km)
!    z        z-coord of scalar point in computation space (km)
!    label    character string describing the contents of a plot
!    time     model runing time step
!    slicopt  slice orientation indicator
!             = 1, x-y slice of at k=kslice is plotted.
!             = 2, x-z slice of at j=jslice is plotted.
!             = 3, y-z slice of at i=islice is plotted.
!             = 4, horizontal slice at z index islice is plotted.
!             = 5, xy-z cross section of wind islice is plotted.
!             = 6, data field on constant p-level is plotted.
!             = 0, all of the three slices above are plotted.
!    axy2d    2d x-y array
!    av2d     2D array for the vertical slice
!    xp       x-coordinate of grid points on arbitary vertical
!               cross-section
!    yp       y-coordinate of grid points on arbitary vertical
!               cross-section
!    zp       z-coordinate of grid points on arbitary vertical
!               cross-section
!    runname  character string describing the model run
!    factor   scaling factor
!    hterain  2-D terrain data for contour
!    trnplt   flag to plot terrain (0/1)
!  WORK ARRAY:
!    tem1     Temporary work array.
!    tem2     Temporary work array.
!    tem3     Temporary work array.
!    tem4     Temporary work array.
!    tem5     Temporary work array.
!
!  (These arrays are defined and used locally (i.e. inside this
!   subroutine), they may also be passed into routines called by
!   this one. Exiting the call to this subroutine, these temporary
!   work arrays may be used for other purposes therefore their
!   contents overwritten. Please examine the usage of work arrays
!   before you alter the code.)
!
!   pp01      The pressure (mb) value at the specific p-level
!   ercpl     reciprocal of exponent
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
!
  INTEGER :: nx,ny,nz
  INTEGER :: n

  REAL :: b(nx,ny,nz)
  REAL :: x(nx,ny,nz)
  REAL :: y(nx,ny,nz)
  REAL :: z(nx,ny,nz)
!
  REAL :: axy2d(nx,ny)
  REAL :: av2d(n,nz),zp(n,nz)
  REAL :: xp(n),yp(n)

  REAL :: x1,x2,dx,y1,y2,dy,z1,z2,dz
  INTEGER :: ibgn,iend,jbgn,jend,kbgn,kend,length
!
  CHARACTER (LEN=6) :: timhms
  CHARACTER (LEN=*) :: label
!
  REAL :: time
!
  INTEGER :: slicopt
  INTEGER :: kslice,jslice,islice
!
  CHARACTER (LEN=*) :: runname
!
  REAL :: factor
!
  INTEGER :: trnplt            ! plot terrain option (0/1/2/3)
  INTEGER :: pltopt            ! plot variable option (0/1/2/3)
  REAL :: hterain(nx,ny)       ! The height of the terrain.
!
  REAL :: tem1(*)
  REAL :: tem2(*)
  REAL :: tem3(*)
  REAL :: tem4(*)
  REAL :: bb(nx,ny,nz)
  REAL :: tem5(*)
!
!-----------------------------------------------------------------------
!
!  Some constants
!
!-----------------------------------------------------------------------
!
  REAL :: pp01, ercpl
  PARAMETER (ercpl=0.3678794)              ! exp(-1.0)
!
!-----------------------------------------------------------------------
!
!  Common blocks for plotting control parameters
!
!-----------------------------------------------------------------------
!
  REAL :: x01,y01                  ! the first  point of interpolation
  REAL :: x02,y02                  ! the second point of interpolation
  REAL :: zlevel                   ! the given height of the slice
  REAL :: sinaf,cosaf,dist,sqrtdxy
  COMMON /slicev/x01,y01,x02,y02,sinaf,cosaf,dist,sqrtdxy
  COMMON /sliceh/zlevel

  INTEGER :: ovrtrn               ! overlay terrain option (0/1)
  REAL :: trninc,trnmin, trnmax   ! terrain interval minimum, maximum
  REAL :: ztmin,ztmax
  COMMON /trnpar/ trnplt,ovrtrn,trninc,trnmin, trnmax,ztmin,ztmax

  INTEGER :: smooth
  COMMON /smoothopt/smooth
!
  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
  CHARACTER (LEN=4) :: stem2
  CHARACTER (LEN=1) :: stem1
  REAL :: x_tmp
  COMMON /tmphc2/ x_tmp

  REAL :: tmpx, tmpy
  CHARACTER (LEN=20) :: distc
  REAL :: x101, y101, x102,y102
  COMMON /slicev1/x101, y101, x102,y102
!
!-----------------------------------------------------------------------
!
!  Misc. local variables
!
!-----------------------------------------------------------------------
!
  INTEGER :: i,j,k,ij,ik,jk,isize,jsize,ksize, llabel
  CHARACTER (LEN=120) :: label_copy
  CHARACTER (LEN=120) :: title

  INTEGER :: wrtflag
  CHARACTER (LEN=80) :: levlab
  CHARACTER (LEN=50) :: timelab
  CHARACTER (LEN=25) :: timestring
  COMMON /timelev/wrtflag, timelab, levlab, timestring

!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  isize = (iend-ibgn)+1
  jsize = (jend-jbgn)+1
  ksize = (kend-kbgn)+1

  label_copy = label
  llabel = 120
  CALL xstrlnth(label_copy, llabel)

  CALL xpscmnt('Start plotting '//label_copy(1:llabel))
!
!-----------------------------------------------------------------------
!
!  Set up terrain, if needed.
!
!-----------------------------------------------------------------------
!
  IF(trnplt == 1 .OR.trnplt == 2 .OR. ovrtrn == 1)  THEN
!   IF(ovrtrn.eq.1)  THEN
    DO j=jbgn,jend
      DO i=ibgn,iend
        ij = i-ibgn+1 + (j-jbgn)*isize
        tem4(ij)=hterain(i,j)
      END DO
    END DO
  END IF

  CALL cvttim( time, timhms )

  IF( timhms(1:1) == '0' ) timhms(1:1)=' '

  WRITE(timelab,'(''T='',F8.1,A)') time,                                &
      ' s ('//timhms(1:2)//':'//timhms(3:4)//':'//timhms(5:6)//')'
  CALL get_time_string ( time, timestring)

  IF ( slicopt == 2 .OR. slicopt == 3 .OR. slicopt == 5 .OR. & 
       slicopt == 10 .OR. slicopt == 11) THEN
    CALL cal_dist(haxisu,dx,dy,x01,y01,x02,y02,                         &
                  slicopt,tmpx,tmpy,distc)
  END IF


  IF(slicopt == 1 .OR. slicopt == 0 ) THEN

    k = kslice
    DO j=jbgn,jend
      DO i=ibgn,iend
        ij = i-ibgn+1 + (j-jbgn)*isize
        tem1(ij) = -9999.0
        IF(        tem2(ij)=x(i,j,k)
        tem3(ij)=y(i,j,k)
      END DO
    END DO


    IF (k /= 2) THEN
      WRITE(levlab,'(''GRID LEVEL='',I3)')k
      WRITE(title,'(a)') label
    ELSE
      WRITE(levlab,'(''FIRST LEVEL ABOVE GROUND (SURFACE)'')')
      WRITE(title,'(a)') label
    END IF


    length = 120
    CALL strlnth( title, length)
    CALL strmin ( title, length)

    DO i=1,smooth
      CALL smooth9pmv(tem1,isize,jsize,1,isize,1,jsize,tem5)
    END DO

    CALL ctr2d(tem1,tem2,tem3, x1,x2,dx, y1,y2,dy,                      &
        isize,jsize,title(1:length),runname,                            &
        tem4,slicopt,pltopt)

!
!-----------------------------------------------------------------------
!
!  slicopt=2   Plot x-z cross-section
!
!-----------------------------------------------------------------------
!
  ELSE IF (slicopt == 2 .OR. slicopt == 0 ) THEN

    x_tmp = y(1,jslice,1)

    j = jslice
    DO k=kbgn,kend
      DO i=ibgn,iend
        ik = i-ibgn+1 + (k-kbgn)*isize
        tem1(ik) = -9999.0
        IF(        tem2(ik)=x(i,j,k)
        tem3(ik)=z(i,j,k)
      END DO
    END DO


    dist = (j-1.5)*tmpy
    length=20
    CALL strmin ( distc, length)
    WRITE(levlab,'(''X-Z PLANE AT Y='',F8.1,A)')dist,distc(1:length)

    WRITE(title,'( a )') label

    length = 120
    CALL strlnth( title, length)
    CALL strmin ( title, length)

    DO i=1,smooth
      CALL smooth9pmv(tem1,isize,ksize,1,isize,1,ksize,tem5)
    END DO

    CALL ctr2d(tem1,tem2,tem3, x1,x2,dx, z1,z2,dz,                      &
        isize,ksize,title(1:length),runname,                            &
        tem4,slicopt,pltopt)

!
!-----------------------------------------------------------------------
!
!  slicopt=3   Plot y-z cross-section
!
!-----------------------------------------------------------------------
!
  ELSE IF ( slicopt == 3 .OR. slicopt == 0) THEN

    x_tmp = x(islice,1,1)

    i = islice
    DO k=kbgn,kend
      DO j=jbgn,jend
        jk = j-jbgn+1 + (k-kbgn)*jsize
        tem1(jk) = -9999.0
        IF(        tem2(jk)=y(i,j,k)
        tem3(jk)=z(i,j,k)
      END DO
    END DO

    dist = (i-1.5)*tmpx
    length=20
    CALL strmin ( distc, length)
    WRITE(levlab,'(''Y-Z PLANE AT X='',F8.1,A)')dist,distc(1:length)

    WRITE(title,'( a )' ) label

    length = 120
    CALL strlnth( title, length)
    CALL strmin ( title, length)

    DO i=1,smooth
      CALL smooth9pmv(tem1,jsize,ksize,1,jsize,1,ksize,tem5)
    END DO

    CALL ctr2d(tem1,tem2,tem3, y1,y2,dy, z1,z2,dz,                      &
        jsize,ksize,title(1:length),runname,                            &
        tem4,slicopt,pltopt)
!
!-----------------------------------------------------------------------
!
!  slicopt=4   Plot horizontal slice at given height
!  slicopt=6   Plot constant pressure slice at given pressure(mb)
!  slicopt=7   Plot isentropic surfaces
!
!-----------------------------------------------------------------------
!
  ELSE IF( slicopt == 4.OR.slicopt == 6.OR.slicopt == 7 ) THEN

    DO k=kbgn,kend
      DO j=jbgn,jend
        DO i=ibgn,iend
          bb(i,j,k) = -9999.0
          IF(        END DO
      END DO
    END DO

    CALL hintrp1(nx,ny,nz,kbgn,kend,bb,z,zlevel,axy2d)

    DO j=jbgn,jend
      DO i=ibgn,iend
        ij = i-ibgn+1 + (j-jbgn)*isize
        tem1(ij)=axy2d(i,j)
        tem2(ij)=x(i,j,2)
        tem3(ij)=y(i,j,2)
      END DO
    END DO

    IF(slicopt == 4) THEN
      WRITE(levlab,'(''Z='',F7.3,'' KM MSL'')')                         &
            zlevel
    ELSE IF(slicopt == 6) THEN
      pp01 = 0.01*ercpl**zlevel
      WRITE(levlab,'(''P='',F7.2,A)') pp01, ' MB'
    ELSE
      WRITE(levlab,'(''THETA='',F5.1,A)') zlevel, ' (K)'
    END IF

    WRITE(title,'(a)') label

    length = 120
    CALL strlnth( title, length)
    CALL strmin ( title, length)

    DO i=1,smooth
      CALL smooth9pmv(tem1,isize,jsize,1,isize,1,jsize,tem5)
    END DO

    CALL ctr2d(tem1,tem2,tem3, x1,x2,dx, y1,y2,dy,                      &
        isize,jsize,title(1:length),runname,                            &
        tem4,slicopt,pltopt)

!
!-----------------------------------------------------------------------
!
!  slicopt=5   Plot vectical slice through two given points
!
!-----------------------------------------------------------------------
!
  ELSE IF( slicopt == 5 ) THEN

    CALL sectvrt(nx,ny,nz,b,x,y,z,av2d,zp,n,xp,yp)

    DO k=kbgn,kend
      DO i=ibgn,iend
        ik = i-ibgn+1 + (k-kbgn)*isize
        tem1(ik) = -9999.0
        IF(av2d(i,k) /= -9999.0) tem1(ik)=av2d(i,k)*factor
        tem2(ik)=x1+(i-ibgn)* sqrtdxy
        tem3(ik)=zp(i,k)
      END DO
    END DO


    IF(xfmat == -1) THEN
      length=20
      CALL strmin ( distc, length)
      WRITE(levlab,                                                     &
          '(''VERTICAL PLANE FROM '',4(A,F8.1),A,A)')                   &
          '(',x101,',',y101,') through (',x102,',',y102,') ',           &
          distc(1:length)
      WRITE(title,'(a)') label
      WRITE(title,'(a)') label
    ELSE IF(xfmat == 0) THEN
      length=20
      CALL strmin ( distc, length)
      WRITE(levlab,                                                     &
          '(''VERTICAL PLANE FROM '',4(A,I5),A,A)')                     &
          '(',INT(x101),',',INT(y101),') through (',INT(x102),','       &
          ,INT(y102),') ', distc(1:length)
      WRITE(title,'(a)') label
    ELSE
      WRITE(title,'(a)') label
      WRITE(stem1,'(i1)')xfmat
      WRITE(stem2,43)stem1
      43       FORMAT('f8.',a1)
      WRITE(title,'(a)') label
    END IF


    length = 120
    CALL strlnth( title, length)
    CALL strmin ( title, length)

    DO i=1,smooth
      CALL smooth9pmv(tem1,isize,ksize,1,isize,1,ksize,tem5)
    END DO

    CALL ctr2d(tem1,tem2,tem3, x1,x2,sqrtdxy, z1,z2,dz,                 &
        isize,ksize,title(1:length),runname,                            &
        tem4,slicopt,pltopt)
!
!-----------------------------------------------------------------------
!
!  slicopt=9   Plot x-z cross-section of the soil model 
! 
!  06/03/2002 Zuwen He 
!
!  slicopt (9) is the same as slicopt (1), except that 
!  the labels. 
!
!-----------------------------------------------------------------------
!
  ELSE IF(slicopt == 9) THEN  

    k = kslice
    DO j=jbgn,jend
      DO i=ibgn,iend
        ij = i-ibgn+1 + (j-jbgn)*isize
        tem1(ij) = -9999.0
        IF(        tem2(ij)=x(i,j,k)
        tem3(ij)=y(i,j,k)
      END DO
    END DO

    WRITE(levlab,'(''GRID LEVEL (SOIL) ='',I3)')k
    WRITE(title,'(a)') label

    length = 120
    CALL strlnth( title, length)
    CALL strmin ( title, length)

    DO i=1,smooth
      CALL smooth9pmv(tem1,isize,jsize,1,isize,1,jsize,tem5)
    END DO

    CALL ctr2d(tem1,tem2,tem3, x1,x2,dx, y1,y2,dy,                      &
        isize,jsize,title(1:length),runname,                            &
        tem4,slicopt,pltopt)
!
!-----------------------------------------------------------------------
!
! Zuwen He, 06/06/2002 
!
!  slicopt=10  Plot x-z cross-section of the soil model.
!
!-----------------------------------------------------------------------
!
  ELSE IF (slicopt == 10) THEN

    x_tmp = y(1,jslice,1)

    j = jslice
    DO k=kbgn,kend
      DO i=ibgn,iend
        ik = i-ibgn+1 + (k-kbgn)*isize
        tem1(ik) = -9999.0
        IF(        tem2(ik)=x(i,j,k)
        tem3(ik)=z(i,j,k)
      END DO
    END DO

    dist = (j-1.5)*tmpy
    length=20
    CALL strmin ( distc, length)
    WRITE(levlab,'(''X-Z PLANE (SOIL) AT Y='',F8.1,A)')dist,distc(1:length)

    WRITE(title,'( a )') label

    length = 120
    CALL strlnth( title, length)
    CALL strmin ( title, length)

    DO i=1,smooth
      CALL smooth9pmv(tem1,isize,ksize,1,isize,1,ksize,tem5)
    END DO

    CALL ctr2d(tem1,tem2,tem3, x1,x2,dx, z1,z2,dz,                      &
        isize,ksize,title(1:length),runname,                            &
        tem4,slicopt,pltopt)
!
!-----------------------------------------------------------------------
!
!  slicopt=11   Plot y-z cross-section of the soil model
!
!-----------------------------------------------------------------------
!
  ELSE IF ( slicopt == 11) THEN

    x_tmp = x(islice,1,1)

    i = islice
    DO k=kbgn,kend
      DO j=jbgn,jend
        jk = j-jbgn+1 + (k-kbgn)*jsize
        tem1(jk) = -9999.0
        IF(        tem2(jk)=y(i,j,k)
        tem3(jk)=z(i,j,k)
      END DO
    END DO

    dist = (i-1.5)*tmpx
    length=20
    CALL strmin ( distc, length)
    WRITE(levlab,'(''Y-Z PLANE (SOIL) AT X='',F8.1,A)')dist,distc(1:length)

    write (*,*) "levlab", levlab

    WRITE(title,'( a )' ) label

    length = 120
    CALL strlnth( title, length)
    CALL strmin ( title, length)

    DO i=1,smooth
      CALL smooth9pmv(tem1,jsize,ksize,1,jsize,1,ksize,tem5)
    END DO

    CALL ctr2d(tem1,tem2,tem3, y1,y2,dy, z1,z2,dz,                      &
        jsize,ksize,title(1:length),runname,                            &
        tem4,slicopt,pltopt)

  END IF

  CALL xpscmnt('End plotting '//label_copy(1:llabel))

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


SUBROUTINE ctr2d(a,x,y,xl,xr,dx,yb,yt,dy,m,n,title,runname,             & 9,26
           hterain,slicopt,pltopt)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!    Generate contour plots of 2-d field A given its coordinates
!      using ZXPLOT package..
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Ming Xue
!
!  MODIFICATION HISTORY:
!
!  6/08/92 (K. Brewster)
!  Added full documentation.
!
!  8/28/94 (M. Zou)
!  Added color routing , overlay terrain.
!
!  1/24/96 (J. Zong and M. Xue)
!  Fixed a problem related to finding the minimum and maximum of the
!  2D array, a, when there exist missing data. Initial min. and max.
!  should be set to values other than the missing value, -9999.0.
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!    a        2-dimensional slice of data to contour
!
!    x        x coordinate of grid points in plot space (over on page)
!    y        y coordinate of grid points in plot space (up on page)
!
!
!    xl       Left bound of the physical domain
!    xr       Right bound of the physical domain
!    dx    Spacing between x-axis tick marks
!    yb       Bottom bound of the physical domain.
!    yt       Top bound of the physical domain.
!    dy    Spacing between y-axis tick marks
!
!    m        first dimension of a
!    n        second dimension of a
!
!    title    character string describing the contents of a
!    runname  character string describing the model run
!
!    hterain  2-D terrain data to contour
!    slicopt  slice orientation indicator
!             = 1, x-y slice of at k=kslice is plotted.
!             = 2, x-z slice of at j=jslice is plotted.
!             = 3, y-z slice of at i=islice is plotted.
!             = 4, horizontal slice at z index islice is plotted.
!             = 5, xy-z cross section of wind islice is plotted.
!             = 6, data field on constant p-level is plotted.
!             = 0, all of the three slices above are plotted.
!    plot      variable plot option (0/1/2/3)
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
  INCLUDE 'arpsplt.inc'
!
!
  INTEGER :: m,n
!
  REAL :: a(m,n)
  REAL :: x(m,n)
  REAL :: y(m,n)
  REAL :: xl,xr,dx,yb,yt,dy
  REAL :: hterain(m,n)

!
  CHARACTER (LEN=*) :: runname
  CHARACTER (LEN=*) :: title
!
!-----------------------------------------------------------------------
!
!  Plotting control common blocks
!
!-----------------------------------------------------------------------
!
  INTEGER :: layover
  COMMON /laypar/ layover

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

!
  INTEGER :: ovrstaopt
  INTEGER :: ovrstam,staset,ovrstan,ovrstav,stacol,markprio,wrtstax
  INTEGER :: nsta_typ,sta_typ(10),sta_marktyp(10),sta_markcol(10)
  REAL :: sta_marksz(10),wrtstad
  COMMON /sta_par/ ovrstaopt,ovrstam,staset,ovrstan,ovrstav,stacol,     &
         markprio,nsta_typ,sta_typ,sta_marktyp,                         &
         sta_markcol,sta_marksz,wrtstax,wrtstad
!
  REAL :: ctinc,ctmin,ctmax,vtunt  ! contour interval and vector unit
  COMMON /incunt/ ctinc,ctmin,ctmax,vtunt
!
  INTEGER :: icolor,icolor1,lbcolor,trcolor                ! required color
  COMMON /recolor/icolor,icolor1,lbcolor,trcolor
!
  INTEGER :: flag
  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
!
  REAL :: yxratio
  COMMON /yratio/ yxratio       ! the scaling factor the y/x ratio.
!
  INTEGER :: ntitle,titcol, nxpic, nypic, wpltime
  REAL :: titsiz
  CHARACTER (LEN=132) :: ptitle(3), footer_l, footer_c, footer_r

  COMMON /titpar1/ptitle, footer_l, footer_c, footer_r
  COMMON /titpar2/ntitle,titcol,wpltime, nxpic, nypic
  COMMON /titpar3/titsiz
!
!-----------------------------------------------------------------------
!
!  Misc. local variables
!
!-----------------------------------------------------------------------
!
  CHARACTER (LEN=150) :: ch1
  CHARACTER (LEN=150) :: ch

  INTEGER, ALLOCATABLE :: iwrk(:)
  REAL   , ALLOCATABLE :: xwk(:),ywk(:)
  INTEGER :: istatus

  INTEGER :: i,j
  REAL :: cl(500)       ! contour levels
  REAL :: pl,pr,pb,pt   ! plot space left, right, bottom, top coordinate
  REAL :: px,py         ! plot space left-right length and up-down height
  REAL :: pxc,pyc       ! plot space left-right center and
!                                   up-down    center
  REAL :: xs,ys         ! real space left-right length and up-down height
  REAL :: zinc          ! contour interval
  REAL :: zmin,zmax     ! max and min of data array
  INTEGER :: ncl,slicopt,mode1

  REAL :: zlevel
  COMMON/sliceh/zlevel
  INTEGER :: pltopt     ! variavle plot option (0/1/2/3)

  INTEGER :: timeovr
  COMMON /timover/ timeovr

  REAL :: lblmag, ctrlbsiz, axlbsiz
  COMMON /labmag/ lblmag, ctrlbsiz, axlbsiz
  REAL :: xfinc
!
  INTEGER :: col_table,pcolbar
  COMMON /coltable/col_table,pcolbar
!
  INTEGER :: LEN,len1
!
  CHARACTER (LEN=12) :: varname
  COMMON /varplt1/ varname
!
  CHARACTER (LEN=150) :: f_ch
!
  INTEGER :: setcontopt, setcontnum
  REAL :: setconts(maxunevm,maxuneva)
  COMMON /setcon_par/setcontopt,setcontnum,setconts
  INTEGER :: ncont
  REAL :: tcont(maxunevm)

  INTEGER :: wrtflag
  CHARACTER (LEN=25) :: timestring
  CHARACTER (LEN=80) :: levlab
  CHARACTER (LEN=50) :: timelab
  COMMON /timelev/wrtflag,timelab, levlab, timestring

  CHARACTER (LEN=80) :: prestr
  INTEGER :: preflag
  COMMON /preinfo/ prestr,preflag

  REAL :: x101, y101, x102,y102
  COMMON /slicev1/x101, y101, x102,y102

  REAL :: ytmp     !! local temporary variable
  !wdt update
  REAL :: f_cputime,cpu1, cpu2
  DOUBLE PRECISION :: f_walltime,second1,second2
  REAL :: hatch_angle
!
  INTEGER :: missval_colind, missfill_opt    ! miss value color index
  COMMON /multi_value/ missfill_opt,missval_colind
  INTEGER :: missfill
  DATA missfill/0/
  INTEGER :: mxset

  INTEGER :: xnwpic_called
  COMMON /callnwpic/xnwpic_called

  INTEGER :: iclfrq
  INTEGER :: ctrlbfrq
  COMMON /clb_frq/ ctrlbfrq

!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!
!-----------------------------------------------------------------------
!
!  Check for adequate room in work array
!
!-----------------------------------------------------------------------
!
  second1= f_walltime()
  cpu1 = f_cputime()

  ncont = 0
  ncl = 1

  ALLOCATE(iwrk(m*n),STAT=istatus) 
  ALLOCATE(xwk(m*n),ywk(m*n),STAT=istatus) 

  WRITE(6,'(//a,a)') 'Plotting ',title

  IF( layover == 0  .OR. xnwpic_called == 0) THEN
    CALL xnwpic
    xnwpic_called = 1
    timeovr=0              ! set overlayer terrain agian
    wrtflag = 0             !
    preflag = 0
    prestr = levlab
    len1=80
    CALL strmin(prestr,len1)
  ELSE
    timeovr=1
    wrtflag = wrtflag + 1
  END IF
!
!-----------------------------------------------------------------------
!
!  Get plotting space variables
!
!-----------------------------------------------------------------------
!

  CALL xqpspc( pl, pr, pb, pt)
  px = pr - pl
  py = pt - pb
  pxc = (pr+pl)/2
  pyc = (pb+pt)/2

  xs = xr-xl
  ys = yt-yb
!
!-----------------------------------------------------------------------
!
!  Let the longest lenth determine size scaling of plot
!
!-----------------------------------------------------------------------
!
  IF( py/px >= (ys*yxratio)/xs ) THEN
    py = (ys*yxratio)/xs*px
    CALL xpspac(pl, pr, pyc-py/2, pyc+py/2 )
  ELSE
    px = xs/(ys*yxratio)*py
    CALL xpspac(pxc-px/2, pxc+px/2, pb, pt)
  END IF
!
!-----------------------------------------------------------------------
!
!  Set the real distance to plot distance scaling
!
!-----------------------------------------------------------------------
!
  CALL xmap( xl, xr, yb, yt )
!
  CALL xlbsiz( ctrlbsiz*(yt-yb)*lblmag )
!
!-----------------------------------------------------------------------
!
!  Find max and min of data array
!
!-----------------------------------------------------------------------
!
  mxset = 0
  missfill = 0
  DO j=1,n
    DO i=1,m
      IF(ABS(        missfill = 1
        CYCLE
      END IF
      IF( mxset == 0) THEN
        zmax= a(i,j)
        zmin= a(i,j)
        mxset = 1
      ELSE
        zmax= MAX (zmax,a(i,j))
        zmin= MIN (zmin,a(i,j))
      END IF
    END DO
  END DO

  IF (missfill == 1 .AND. missfill_opt == 1)                            &
      CALL fillmissval ( m,n,xl, xr, yb,yt )
  CALL xcolor(lbcolor)
!
!-----------------------------------------------------------------------
!
!  Find proper contour interval and then contour field
!     using ZXPLOT routine xconta
!
!-----------------------------------------------------------------------

  cl(1)=0.0
  IF( zmax-zmin > 1.0E-20 ) THEN
!
!-----------------------------------------------------------------------
!
!    Check to see if user defined contour levels is available for the
!    current variable.
!
!-----------------------------------------------------------------------


    CALL get_contour (ncont, tcont)
    iclfrq = ctrlbfrq

    IF(setcontopt > 0 .AND.ncont > 0) THEN
      ch1(1:11)=' contours: '
      LEN=11
      DO i =1,ncont
        CALL xrch1(tcont(i),f_ch,len1)
        WRITE(ch1,'(a,a,'' '')')ch1(1:LEN), f_ch(1:len1)
        LEN=LEN+len1+1
      END DO
      DO i=1,ncont
        cl(i)=tcont(i)
      END DO
      ncl = ncont
      mode1 = 4
      iclfrq = 1
      GO TO 150
    END IF

    IF( ctinc == 0.0) THEN
      cl(2)=cl(1)+ xfinc(zmax-zmin)/2
      IF(cl(2)-cl(1) == 0.0) cl(2)=cl(1)+1.0
      mode1=1
      CALL xnctrs( 8,20)
    ELSE IF ( ctinc == -9999.) THEN
      CALL xnctrs( 8,20)
      CALL set_interval(a, m,n,zmin,zmax,ctmin,ctmax,cl)
      ctinc = cl(2)-cl(1)
      zinc = ctinc
      mode1=1
    ELSE
      cl(2)=cl(1)+ctinc
      CALL xnctrs(1,300)
      mode1=1
    END IF

    150     CONTINUE
    zinc = cl(2)-cl(1)


!-----------------------------------------------------------------------
!
!  Plot contour or color filled contour fields
!
!-----------------------------------------------------------------------

    CALL xwindw(xl, xr, yb, yt)

    CALL xctrlim(ctmin,ctmax)
    CALL xclfrq(iclfrq)


    IF(pltopt == 1) THEN
      CALL xctrclr(icolor, icolor)
      CALL xconta(a,x,y,iwrk,m,m,n,cl,ncl,mode1)
    ELSE IF( pltopt == 2) THEN
      CALL xctrclr(icolor, icolor1)
      CALL xcolfil(a,x,y,iwrk,xwk,ywk,m,m,n,cl,ncl,mode1)
      CALL xchsiz(0.025*(yt-yb))
      CALL xcpalet(pcolbar)
    ELSE IF(pltopt == 4) THEN
      CALL xctrclr(icolor, icolor1)
      CALL xconta(a,x,y,iwrk,m,m,n,cl,ncl,mode1)
    ELSE IF(pltopt == 5) THEN
      CALL xctrclr(icolor, icolor1)
      CALL xcolfil(a,x,y,iwrk,xwk,ywk,m,m,n,cl,ncl,mode1)
      CALL xchsiz(0.025*(yt-yb))
      CALL xcpalet(pcolbar)
      CALL xctrclr(lbcolor, lbcolor)
      CALL xconta(a,x,y,iwrk,m,m,n,cl,ncl,mode1)
    ELSE IF(pltopt == 6) THEN
      CALL xctrclr(icolor, icolor)

      CALL xdhtch(0.003)

      CALL xctrclr(icolor, icolor)
      ncl = 2
      mode1 = 4
      cl(1) = ctmin
      cl(2) = ctmax
      CALL xclfrq(1)
      CALL xhilit(0)
      CALL xconta(a,x,y,iwrk,m,m,n,cl,ncl,mode1)
      CALL xhilit(1)

      hatch_angle = 45.0
      CALL xdhtch(0.004)
      CALL xhatcha(a,x,y,xwk,ywk,m,m,n,ctmin,1.0E10,hatch_angle)

      CALL xdhtch(0.002)
      CALL xhatcha(a,x,y,xwk,ywk,m,m,n,ctmax,1.0E10,hatch_angle)
    END IF

    CALL xclfrq(2)

  ELSE
    cl(2)=1.0
    ncl=2
  END IF

  IF(ctinc == 0.0) THEN
    zinc = cl(2) - cl(1)
  ELSE
    zinc = ctinc
  END IF
!
!-----------------------------------------------------------------------
!
!  Plot map, boxes and polygons.
!
!-----------------------------------------------------------------------
!
  CALL pltextra(slicopt, pltopt)

!-----------------------------------------------------------------------
!
!  Terrain outline in vertical slices.
!
!-----------------------------------------------------------------------

  IF(slicopt == 2 .OR. slicopt == 3 .OR. slicopt == 5 .OR.  & 
     slicopt == 10 .OR. slicopt == 11) THEN
    CALL xcolor(trcolor)
    CALL xthick(3)
    CALL xpenup( x(1,1), y(1,1)-0.5*(y(1,2)-y(1,1)) )
    DO i=2,m
      CALL xpendn( x(i,1), y(i,1)-0.5*(y(i,2)-y(i,1)) )
    END DO
    CALL xthick(1)
  END IF
!
!-----------------------------------------------------------------------
!
!  Overlay terrain contour if required in x-y level
!  or Plot terrain outline in slice zlevel
!
!-----------------------------------------------------------------------
!
  IF ( timeovr == 0 ) CALL plttrn(hterain,x,y,m,n,slicopt)
!
!-----------------------------------------------------------------------
!
!  Plot observations
!
!-----------------------------------------------------------------------
!
  IF(obsset == 1 .AND. ovrobs == 1 .AND.  & 
     (slicopt == 1 .OR. slicopt == 4 .OR. slicopt == 6 .OR. & 
      slicopt == 7 .OR. slicopt == 8 .OR. slicopt ==  9) ) THEN
    CALL xchsiz(0.025*ys * lblmag)
    CALL pltobs(1)
    obsset=0
  END IF
!
!-----------------------------------------------------------------------
!
!  Plot station labels
!
!-----------------------------------------------------------------------
!
  IF( ovrstaopt == 1 .AND. wrtstax == 1                                 &
        .AND. (timeovr == 0 .OR. (timeovr == 1 .AND. pltopt == 2))      &
        .AND. (slicopt == 2 .OR. slicopt == 3 .OR. slicopt == 5 .OR.   & 
               slicopt == 10 .OR. slicopt == 11) ) THEN
    CALL xchsiz(0.025*ys * lblmag)
    flag=1
    CALL pltsta(a,a,x,y,m,n,flag,slicopt)
  END IF

  IF(ovrstaopt == 1 .AND. staset == 1 .AND.                             &
        (ovrstam == 1 .OR. ovrstan == 1 .OR. ovrstav == 1)              &
        .AND.(slicopt == 1 .OR. slicopt == 4 .OR. slicopt == 6 .OR.     & 
              slicopt == 7 .OR. slicopt == 8 .OR. slicopt == 9)         &
        .AND.(timeovr == 0.OR.(timeovr == 1.AND.pltopt == 2))) THEN
    CALL xchsiz(0.025*ys * lblmag)
    flag=0
    CALL pltsta(a,a,x,y,m,n,flag,slicopt)
!    staset=0
  END IF

  CALL xwdwof
!
!-----------------------------------------------------------------------
!
!  Plot axes with tick marks
!
!-----------------------------------------------------------------------
!
  CALL pltaxes(slicopt,dx,dy)

!
  IF(ntitle > 0 .AND. nxpic == 1 .AND. nypic == 1                       &
                   .AND. timeovr == 0 ) THEN
    CALL xcolor(titcol)
    CALL xchsiz(0.025*ys * titsiz)
    DO i=1,ntitle
      LEN=132
      CALL strlnth(ptitle(i),LEN)
      CALL xchori(0.)
      CALL xcharc( xl+xs/2,yt+(0.25-(i-1)*0.06)*ys,ptitle(i)(1:LEN))
    END DO
    CALL xcolor(lbcolor)
  END IF

  CALL xchsiz( 0.030*ys * lblmag )

! plot time and level label
  IF ( layover < 1 ) THEN
    IF(levlab /= ' ') THEN
      len1=80
      CALL strmin(levlab,len1)
      CALL xcharc((xl+xr)*0.5,yt+0.015*ys, levlab(1:len1))
      preflag = 1
    END IF
    len1=50
    CALL strmin(timelab,len1)
    CALL xcharc((xl+xr)*0.5,yt+0.06*ys,                                 &
                 timestring(1:25)//'  '//timelab(1:len1))
  END IF

  IF(preflag == 0 .AND. levlab /= ' ') THEN
    len1=80
    CALL strmin(levlab,len1)
    CALL xcharc((xl+xr)*0.5,yt+0.015*ys, levlab(1:len1))
    preflag = 1
  END IF

  LEN = 120
  CALL strmin(title, LEN)
  IF( title(LEN-1:LEN-1) == ')' ) THEN
    title(LEN-1:LEN-1)=','
  ELSE
    title(LEN:LEN+1)=' ('
    LEN = LEN+1
  END IF

  IF(pltopt == 2) THEN
    WRITE(f_ch, '(a,''SHADED)'')')title(1:LEN)
  ELSE IF( pltopt == 5 ) THEN
    WRITE(f_ch, '(a,''SHADED/CONTOUR)'')')title(1:LEN)
  ELSE
    WRITE(f_ch, '(a,''CONTOUR)'')')title(1:LEN)
  END IF

! if first levlab is not equal second levlab then attatch levlab on f_ch
  LEN=120
  CALL strmin(f_ch, LEN)
  len1=80
  CALL strmin(levlab,len1)
!
!mx
!  IF (preflag.eq.1 .and. prestr(1:len1).ne.levlab(1:len1)
!    :   .and. prestr(1:1).ne.' '
!    :   .and.layover.ne.0 .and. levlab(1:1).ne.' ') THEN
!    write(f_ch,'(a,a)') f_ch(1:len),levlab(1:len1)
!  ENDIF

!   IF(pltopt.eq.1 .and. layover.ge.1 ) CALL xcolor(icolor)
  IF(pltopt == 1) CALL xcolor(icolor)
! plot variable name
  IF (preflag == 1 .AND. prestr /= levlab .AND. prestr /= ' '           &
        .AND.layover /= 0 .AND. levlab /= ' ') THEN
    CALL xchsiz( 0.018*ys * lblmag )
  ELSE
    CALL xchsiz( 0.028*ys * lblmag )
  END IF
! save for next time use
  IF(prestr(1:1) == ' ' .AND.layover /= 0 ) prestr=levlab

  IF(lbaxis == 1 ) THEN
    IF( wrtstax == 0) THEN
      ytmp = 0.08
    ELSE
      ytmp = 0.14
    END IF
  ELSE
    ytmp = 0.12
  END IF

  LEN=150
  CALL strmin(f_ch,LEN)

  CALL xchsiz( 0.025*ys * lblmag )
!mx
  CALL xcolor(lbcolor)

  CALL xcharl(xl-0.20*(xr-xl), yb-(ytmp+wrtflag*0.030)*ys ,             &
       f_ch(1:LEN))

  IF (pltopt == 1.OR.pltopt == 3.OR.pltopt == 4.OR.pltopt == 5)THEN
    IF(ABS(zmin-zmax) <= 1.e-15 .OR. ncont > 0)  THEN
      WRITE(ch,'(''MIN='',G9.3E2,'' MAX='',G9.3E2)')zmin,zmax
    ELSE
      WRITE(ch,'(''MIN='',G10.4E2,'' MAX='',G10.4E2,                    &
      &   '' inc='',g10.4E2)')zmin,zmax,zinc
    END IF
  ELSE IF( pltopt == 2 ) THEN
    WRITE(ch,'(''MIN='',G9.3E2,'' MAX='',G9.3E2)')zmin,zmax
  END IF

  LEN=150
  CALL strmin(ch,LEN)
  CALL xcharr(xr+0.20*(xr-xl), yb-(ytmp+wrtflag*0.030)*ys ,             &
       ch(1:LEN))
  IF (ncont > 1 .AND. (pltopt == 1 .OR. pltopt == 4) ) THEN
    wrtflag = wrtflag+1
    len1=150
    CALL strmin(ch1,len1)
    CALL xcharr(xr+0.20*(xr-xl),yb-(ytmp+wrtflag*0.030)*ys,             &
         ch1(1:len1))
  END IF

!-----------------------------------------------------------------------
!
!  Plot additional text below the figure
!
!-----------------------------------------------------------------------

  CALL label2d(runname)

  CALL xpspac(pl, pr, pb, pt)  ! set frame back

  cpu2 = f_cputime()
  second2 = f_walltime()
!  write(6,*)'!!!!  total cpu time for one CTR2D  :',
!    :    cpu2-cpu1,' PLOT:',varname

  DEALLOCATE(iwrk)
  DEALLOCATE(xwk,ywk)

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


SUBROUTINE ctrinc( ctinc0, ctmin0, ctmax0 ) 3

!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!    Set the contour interval for field to plotted by CTR2D.
!
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Ming Xue
!
!  MODIFICATION HISTORY:
!    6/08/92  Added full documentation (K. Brewster)
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!    ctinc0    Contour interval
!              If CTINC0 = 0.0, the interval is internally determined.
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
  REAL :: ctinc0,ctmin0,ctmax0
!
!-----------------------------------------------------------------------
!
!  Plotting control common blocks
!
!-----------------------------------------------------------------------
!
  REAL :: ctinc,ctmin,ctmax,vtunt   ! contour interval and vector unit
  COMMON /incunt/ ctinc,ctmin,ctmax,vtunt
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!

  ctinc = ctinc0
  ctmin = ctmin0
  ctmax = ctmax0

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


SUBROUTINE label2d(runname) 2,3
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!  Plot certain text labels for VTR2D and CTR2d.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Ming Xue
!
!  MODIFICATION HISTORY:
!  Taked from former CTR2D.
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!    xl       Left bound of the physical domain
!    xr       Right bound of the physical domain
!    yb       Bottom bound of the physical domain.
!    yt       Top bound of the physical domain.
!
!    runname  character string describing the model run
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE

  CHARACTER (LEN=*) :: runname
  INTEGER :: layover
  COMMON /laypar/ layover

  INTEGER :: icolor,icolor1,lbcolor,trcolor                ! required color
  COMMON /recolor/icolor,icolor1,lbcolor,trcolor
!
  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
!
  INTEGER :: ntitle,titcol, nxpic, nypic, wpltime
  REAL :: titsiz
  CHARACTER (LEN=132) :: ptitle(3), footer_l, footer_c, footer_r

  COMMON /titpar1/ptitle, footer_l, footer_c, footer_r
  COMMON /titpar2/ntitle,titcol,wpltime, nxpic, nypic
  COMMON /titpar3/titsiz

  REAL :: xl,xr,yb,yt
!
!-----------------------------------------------------------------------
!
!  Misc. local variables
!
!-----------------------------------------------------------------------
!
  INTEGER :: nopic
  REAL :: xlimit, ylimit, rotang
  INTEGER :: nhpic, nvpic,ifont

  INTEGER :: ovrtrn ,trnplt       ! overlay terrain option (0/1)
  REAL :: trninc,trnmin, trnmax   ! terrain interval minimum, maximum
  REAL :: ztmin,ztmax
  COMMON /trnpar/ trnplt,ovrtrn,trninc,trnmin, trnmax,ztmin,ztmax
  INTEGER :: timeovr
  COMMON /timover/ timeovr

  REAL :: lblmag, ctrlbsiz, axlbsiz
  COMMON /labmag/ lblmag, ctrlbsiz, axlbsiz
!
  INTEGER :: col_table,pcolbar
  COMMON /coltable/col_table,pcolbar
!
  CHARACTER (LEN=24) :: tzstring
  CHARACTER (LEN=24) :: tz
  CHARACTER (LEN=132) :: datetimestr

  INTEGER :: lnblnk, len1, len2, len3
  CHARACTER (LEN=120) :: string_l, string_c, string_r

  CHARACTER (LEN=8) :: tzone
  CHARACTER (LEN=10) :: cur_time
  CHARACTER (LEN=8) :: cur_date
  INTEGER :: t_values(8)

  REAL :: ytmp, hch
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!

  CALL xqmap(xl,xr,yb,yt)

  CALL xcolor(lbcolor)

  CALL xqnpic(nopic)
  CALL xqspac(nhpic, nvpic, rotang, xlimit, ylimit)

  CALL xchsiz( 0.021*(yt-yb) * lblmag )

  IF(timeovr == 0) THEN
    IF(nopic == nhpic*(nvpic-1)+1 ) THEN

      IF ( wpltime == 1) THEN

        CALL date_and_time(cur_date,cur_time,tzone,t_values)

        IF(t_values(4) == 0) THEN
          tzstring = ' UTC'
        ELSE
          tzstring = ' Local Time'
        END IF

        WRITE (datetimestr,999) 'Plotted ',                             &
            t_values(1),t_values(2),t_values(3),                        &
            t_values(5),t_values(6),tzstring
        999        FORMAT (a, i4.4,'/',i2.2,'/',i2.2,' ',i2.2,':',i2.2,a)
      END IF

      IF ( footer_l == ' ') THEN
        string_l = 'ARPSPLT/ZXPLOT '
      ELSE
        string_l = footer_l
      END IF

      IF( footer_c == '  ') THEN
        string_c = runname
      ELSE
        string_c = footer_c
      END IF

      IF(wpltime == 1 ) THEN
        string_r = datetimestr(:lnblnk(datetimestr))
      ELSE
        string_r = footer_r
      END IF

      CALL xqcfnt(ifont)
      CALL xcfont(xfont)

      ytmp = 0.29

      CALL xqchsz(hch)

      IF ( layover < 1) THEN

        len1=120
        CALL strmin(string_l, len1)
        len2=120
        CALL strmin(string_c, len2)
        len3=120
        CALL strmin(string_r, len3)

        CALL xcharc(xl+0.5*(xr-xl),                                     &
             yb-(ytmp+layover*0.03)*(yt-yb),                            &
             string_l(1:len1)//'  '//string_c(1:len2)//'  '//           &
             string_r(1:len3))

      END IF
      CALL xcfont(ifont)
    END IF
    timeovr=1
  END IF

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


SUBROUTINE vtr3d(u,v,w, x,y,z, xw,xe,dx, ys,yn,dy, zb,zt,dz,            & 3,38
           nx,ibgn,iend,ist, ny,jbgn,jend,jst, nz,kbgn,kend,kst,        &
           kslice, jslice, islice, label,time, runname, factor,         &
           slicopt,n,xp,yp,zp,u1,v1,u2,v2,w2,                           &
           tem1,tem2,tem3,tem4,                                         &
           tem5,tem6,hterain)

!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!    Plot vector fields in 2-d slices
!
!  AUTHOR: Ming Xue
!
!  MODIFICATION HISTORY:
!    6/08/92  Added full documentation (K. Brewster)
!
!  3/25/96 (K. Brewster)
!    Added variables isize,jsize,ksize
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!    u        3-dimensional array of u wind components (m/s)
!    v        3-dimensional array of v wind components (m/s)
!    w        3-dimensional array of w wind components (m/s)
!
!    x        x coordinate of grid points in physical/comp. space (m)
!    y        y coordinate of grid points in physical/comp. space (m)
!    z        z coordinate of grid points in physical space (m)
!
!    xw       value of x for first i grid point to plot
!    xe       value of x for last i grid point to plot
!    ys       value of y for first j grid point to plot
!    yn       value of y for last j grid point to plot
!    zb       value of z for first k grid point to plot
!    zt       value of z for last k grid point to plot
!
!    nx       first dimension of b
!    ibgn     index of first i grid point to plot
!    iend     index of last  i grid point to plot
!
!    ny       second dimension of b
!    jbgn     index of first j grid point to plot
!    jend     index of last  j grid point to plot
!
!    nz       third dimension of b
!    kbgn     index of first k grid point to plot
!    kend     index of last  k grid point to plot
!
!    ist      step size in x direction
!    jst      step size in y direction
!    kst      step size in z direction
!
!    time     time of data in seconds
!
!    kslice   k index of plane for slicopt=1 x-y slice
!    jslice   j index of plane for slicopt=2 x-z slice
!    islice   i index of plane for slicopt=1 y-z slice
!
!    runname  character string decribing run
!
!    factor   scaling factor for winds
!             V*factor wind vectors are plotted
!
!    slicopt  slice orientation indicator
!             = 1, x-y slice of at k=kslice is plotted.
!             = 2, x-z slice of at j=jslice is plotted.
!             = 3, y-z slice of at i=islice is plotted.
!             = 4, horizontal slice at z index islice is plotted.
!             = 5, xy-z cross section of wind islice is plotted.
!             = 6, data field on constant p-level is plotted.
!             = 0, all of the three slices above are plotted.
!
!  WORK ARRAYS:
!
!    tem1     Temporary work array.
!    tem2     Temporary work array.
!    tem3     Temporary work array.
!    tem4     Temporary work array.
!    tem5     Temporary work array.
!
!  (These arrays are defined and used locally (i.e. inside this
!   subroutine), they may also be passed into routines called by
!   this one. Exiting the call to this subroutine, these temporary
!   work arrays may be used for other purposes therefore their
!   contents overwritten. Please examine the usage of work arrays
!   before you alter the code.)
!
!   pp01      The pressure (mb) value at the specific p-level
!   ercpl     reciprocal of exponent
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
!
  INTEGER :: nx,ny,nz
  INTEGER :: n
!
  REAL :: u(nx,ny,nz)
  REAL :: v(nx,ny,nz)
  REAL :: w(nx,ny,nz)

  REAL :: x(nx,ny,nz)
  REAL :: y(nx,ny,nz)
  REAL :: z(nx,ny,nz)
!
  REAL :: u1(nx,ny),v1(nx,ny)
  REAL :: u2(n,nz),v2(n,nz),w2(n,nz),zp(n,nz)
  REAL :: xp(n),yp(n)

  INTEGER :: kslice,jslice,islice
  CHARACTER (LEN=*) :: runname
  CHARACTER (LEN=*) :: label

  REAL :: xw,xe,dx,ys,yn,dy,zb,zt,dz
  INTEGER :: ibgn,iend,ist, jbgn,jend,jst, kbgn,kend,kst

  REAL :: time,factor
  INTEGER :: slicopt

  INTEGER :: iunits, itype
  COMMON /windvtr/iunits, itype

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

  REAL :: xw1,xe1,ys1,yn1
  COMMON /xuvpar/xw1,xe1,ys1,yn1

!
!-----------------------------------------------------------------------
!
!  Some constants
!
!-----------------------------------------------------------------------
!
  REAL :: pp01, ercpl
  PARAMETER (ercpl=0.3678794)              ! exp(-1.0)
!
!-----------------------------------------------------------------------
!
!  Work arrays: tem1,tem2,tem3,tem4,tem5 of size at least
!          max( nx*ny, nx*nz, ny*nz).
!
!-----------------------------------------------------------------------
!
  REAL :: tem1(*),tem2(*),tem3(*),tem4(*),tem5(*)
  REAL :: tem6(*)
!
!-----------------------------------------------------------------------
!
!  Common blocks for plotting control parameters
!
!-----------------------------------------------------------------------
!
  REAL :: x01,y01                  ! the first  point of interpolation
  REAL :: x02,y02                  ! the second point of interpolation
  REAL :: zlevel                   ! the given height of the slice
  REAL :: sinaf,cosaf,dist,sqrtdxy
  COMMON /slicev/x01,y01,x02,y02,sinaf,cosaf,dist,sqrtdxy
  COMMON /sliceh/zlevel

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

!
!-----------------------------------------------------------------------
!
!  Misc. local Variables
!
!-----------------------------------------------------------------------
!
  INTEGER :: i,j,k,ij,ik,jk,istep,jstep,length,isize,jsize,ksize
  REAL :: uunit
  CHARACTER (LEN=6) :: timhms
  CHARACTER (LEN=120) :: title
!
  INTEGER :: icolor,icolor1,lbcolor,trcolor        ! required color
  COMMON /recolor/icolor,icolor1,lbcolor,trcolor
!
  INTEGER :: trnplt                ! flag to plot terain (1 or 0)
  REAL :: hterain(nx,ny)           ! The height of the terrain.

  INTEGER :: ovrtrn         ! overlay terrain option (0/1)

  REAL :: trninc,trnmin, trnmax    ! terrain interval minimum, maximum
  REAL :: ztmin,ztmax
  COMMON /trnpar/ trnplt,ovrtrn,trninc,trnmin, trnmax,ztmin,ztmax
!
  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
  CHARACTER (LEN=6) :: stem2
  CHARACTER (LEN=1) :: stem1

  INTEGER :: smooth
  COMMON /smoothopt/smooth

  INTEGER :: id

  REAL :: x_tmp
  COMMON /tmphc2/ x_tmp

  INTEGER :: wrtflag
  CHARACTER (LEN=80) :: levlab
  CHARACTER (LEN=50) :: timelab
  CHARACTER (LEN=25) :: timestring
  COMMON /timelev/wrtflag, timelab, levlab, timestring

  REAL :: tmpx, tmpy
  CHARACTER (LEN=20) :: distc
  REAL :: x101, y101, x102,y102
  COMMON /slicev1/x101, y101, x102,y102
  INTEGER :: llabel
  CHARACTER (LEN=120) :: label_copy
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  isize=(iend-ibgn)+1
  jsize=(jend-jbgn)+1
  ksize=(kend-kbgn)+1

  label_copy = label
  llabel = 120
  CALL xstrlnth(label_copy, llabel)
  CALL xpscmnt('Start plotting '//label_copy(1:llabel))
!
!-----------------------------------------------------------------------
!
!  slicopt=1   Plot u-v field
!
!-----------------------------------------------------------------------
!
  CALL cvttim ( time, timhms)

  IF( timhms(1:1) == '0' ) timhms(1:1)=' '

  WRITE(timelab,'(''T='',F8.1,A)') time,                                &
      ' s ('//timhms(1:2)//':'//timhms(3:4)//':'//timhms(5:6)//')'
  CALL get_time_string ( time, timestring)

  IF ( slicopt == 2 .OR. slicopt == 3  .OR. slicopt == 5) THEN
    CALL cal_dist(haxisu,dx,dy,x01,y01,x02,y02,slicopt,                 &
                  tmpx,tmpy,distc)
  END IF

!
!-----------------------------------------------------------------------
!
!  Set up terrain, if needed.
!
!-----------------------------------------------------------------------
!
  IF(trnplt == 1 .OR.trnplt == 2 .OR. ovrtrn == 1)  THEN
    DO j=jbgn,jend
      DO i=ibgn,iend
        ij = i-ibgn+1 + (j-jbgn)*isize
        tem5(ij)=hterain(i,j)
      END DO
    END DO
  END IF


  IF( slicopt == 1 .OR. slicopt == 0 ) THEN

    k = kslice
    DO j=jbgn,jend
      DO i=ibgn,iend
        ij = i-ibgn+1 + (j-jbgn)*isize
        tem1(ij) = -9999.0
        tem2(ij) = -9999.0
        IF(        IF(        tem3(ij)=x(i,j,k)
        tem4(ij)=y(i,j,k)
      END DO
    END DO

    IF (k /= 2) THEN
      WRITE(levlab,'(''GRID LEVEL='',I3)')k
      WRITE(title,'(''U-V '',A)')label
    ELSE
      WRITE(levlab,'(''FIRST LEVEL ABOVE GROUND (SURFACE)'')')
      WRITE(title,'(''U-V '',A)') label
    END IF

    length = 120
    CALL strlnth( title, length )
    CALL strmin ( title, length)

    uunit = 10.0
    CALL xvmode(1)
    istep = ist
    jstep = jst

    DO i=1,smooth
      CALL smooth9pmv(tem1,isize,jsize,1,isize,1,jsize,tem6)
      CALL smooth9pmv(tem2,isize,jsize,1,isize,1,jsize,tem6)
    END DO

    CALL vtr2d(tem1,tem2,tem3,tem4, uunit, xw,xe,dx,ys,yn,dy,           &
        isize,istep,jsize,jstep,                                        &
        title(1:length),runname, 1,                                     &
        tem5,slicopt)
!
!-----------------------------------------------------------------------
!
!  slicopt=2   Plot u-w field
!
!-----------------------------------------------------------------------
!
  ELSE IF( slicopt == 2 .OR. slicopt == 0 ) THEN

    x_tmp = y(1,jslice,1)

    j = jslice

    dist = (j-1.5)*tmpy
    length=20
    CALL strmin ( distc, length)
    WRITE(levlab,'(''X-Z PLANE AT Y='',F8.1,A)')dist,distc(1:length)

    IF(varname(1:6) == 'xuvplt') THEN
      xw1=xw
      xe1=xe
      ys1=ys
      yn1=yn
      id=4
      DO k=kbgn,kend
        DO i=ibgn,iend
          ik = i-ibgn+1 + (k-kbgn)*isize
          tem1(ik) = -9999.0
          tem2(ik) = -9999.0
          IF(          IF(          tem3(ik)=x(i,j,k)
          tem4(ik)=z(i,j,k)
        END DO
      END DO
      WRITE(title,'(''U-V '',A)')label
    ELSE
      id=2
      DO k=kbgn,kend
        DO i=ibgn,iend
          ik = i-ibgn+1 + (k-kbgn)*isize
          tem1(ik) = -9999.0
          tem2(ik) = -9999.0
          IF(          IF(          tem3(ik)=x(i,j,k)
          tem4(ik)=z(i,j,k)
        END DO
      END DO
      WRITE(title,'(''U-W '',A)')label
    END IF

    length = 120
    CALL strlnth( title, length )
    CALL strmin ( title, length)

    DO i=1,smooth
      CALL smooth9pmv(tem1,isize,ksize,1,isize,1,ksize,tem6)
      CALL smooth9pmv(tem2,isize,ksize,1,isize,1,ksize,tem6)
    END DO

    uunit = 10.0
    CALL xvmode(1)
    istep = ist
    jstep = kst
    CALL vtr2d(tem1,tem2,tem3,tem4,uunit, xw,xe,dx,zb,zt,dz,            &
         isize,istep,ksize,jstep,                                       &
         title(1:length),runname, id,                                   &
         tem5,slicopt)


!
!-----------------------------------------------------------------------
!
!  slicopt=3   Plot v-w field
!
!-----------------------------------------------------------------------
!
  ELSE IF( slicopt == 3 .OR. slicopt == 0 ) THEN

    x_tmp = y(1,jslice,1)

    i = islice

    dist = (i-1.5)*tmpx
    length=20
    CALL strmin ( distc, length)
    WRITE(levlab,'(''Y-Z PLANE AT X='',F8.1,A)')dist,distc(1:length)

    IF(varname(1:6) == 'xuvplt') THEN
      xw1=xw
      xe1=xe
      ys1=ys
      yn1=yn
      id=4
      DO k=kbgn,kend
        DO j=jbgn,jend
          jk = j-jbgn+1 + (k-kbgn)*jsize
          tem1(jk) = -9999.0
          tem2(jk) = -9999.0
          IF(          IF(          tem3(jk)=y(i,j,k)
          tem4(jk)=z(i,j,k)
        END DO
      END DO
      WRITE(title,'(''U-V '',A)')label
    ELSE
      id=3
      DO k=kbgn,kend
        DO j=jbgn,jend
          jk = j-jbgn+1 + (k-kbgn)*jsize
          tem1(jk) = -9999.0
          tem2(jk) = -9999.0
          IF(          IF(          tem3(jk)=y(i,j,k)
          tem4(jk)=z(i,j,k)
        END DO
      END DO
      WRITE(title,'(''V-W '',A)')label
    END IF

    length = 120
    CALL strlnth( title, length )
    CALL strmin ( title, length)

    DO i=1,smooth
      CALL smooth9pmv(tem1,jsize,ksize,1,jsize,1,ksize,tem6)
      CALL smooth9pmv(tem2,jsize,ksize,1,jsize,1,ksize,tem6)
    END DO

    uunit = 10.0
    CALL xvmode(1)
    istep = jst
    jstep = kst
    CALL vtr2d(tem1,tem2,tem3,tem4,uunit, ys,yn,dy,zb,zt,dz,            &
         jsize,istep,ksize,jstep,                                       &
         title(1:length),runname, id,                                   &
         tem5,slicopt)

!
!-----------------------------------------------------------------------
!
!  slicopt=4   Plot u-v field on constant z levels
!  slicopt=6   Plot u-v field on constant pressure levels
!  slicopt=7   Plot u-v field on constant PT levels
!
!-----------------------------------------------------------------------
!
  ELSE IF( slicopt == 4.OR.slicopt == 6.OR.slicopt == 7 ) THEN

!    CALL hintrp(nx,ny,nz,u,z,zlevel,u1)
!    CALL hintrp(nx,ny,nz,v,z,zlevel,v1)

    CALL hintrp1(nx,ny,nz,kbgn,kend,u,z,zlevel,u1)
    CALL hintrp1(nx,ny,nz,kbgn,kend,v,z,zlevel,v1)

    DO j=jbgn,jend
      DO i=ibgn,iend
        ij = i-ibgn+1 + (j-jbgn)*isize
        tem1(ij) = -9999.0
        tem2(ij) = -9999.0
        IF(u1(i,j) /= -9999.0) tem1(ij)=u1(i,j)*factor
        IF(v1(i,j) /= -9999.0) tem2(ij)=v1(i,j)*factor
        tem3(ij)=x(i,j,2)
        tem4(ij)=y(i,j,2)
      END DO
    END DO

    IF( slicopt == 4) THEN
      WRITE(levlab,'(''Z='',F7.3,'' KM MSL'')')                         &
            zlevel
    ELSE IF( slicopt == 6) THEN
      pp01 = 0.01*ercpl**zlevel
      WRITE(levlab,'(''P='',F7.2,A)') pp01, ' MB'
    ELSE
      WRITE(levlab,'(''THETA='',F5.1,A)') zlevel, ' (K)'
    END IF

    WRITE(title,'(''U-V '',A)') label

    length = 120
    CALL strlnth( title, length )
    CALL strmin ( title, length)

    uunit = 10.0
    CALL xvmode(1)
    istep = ist
    jstep = jst

    DO i=1,smooth
      CALL smooth9pmv(tem1,isize,jsize,1,isize,1,jsize,tem6)
      CALL smooth9pmv(tem2,isize,jsize,1,isize,1,jsize,tem6)
    END DO

    CALL vtr2d(tem1,tem2,tem3,tem4, uunit, xw,xe,dx,ys,yn,dy,           &
        isize,istep,jsize,jstep,                                        &
        title(1:length),runname, 1,                                     &
        tem5,slicopt)
!
!-----------------------------------------------------------------------
!
!  slicopt=5   Plot u-v field
!
!-----------------------------------------------------------------------
!
  ELSE IF( slicopt == 5 ) THEN

    CALL sectvrt(nx,ny,nz,u,x,y,z,u2,zp,n,xp,yp)
    CALL sectvrt(nx,ny,nz,v,x,y,z,v2,zp,n,xp,yp)
    CALL sectvrt(nx,ny,nz,w,x,y,z,w2,zp,n,xp,yp)


    IF(varname(1:6) == 'xuvplt') THEN
      xw1=xw
      xe1=xe
      ys1=ys
      yn1=yn
      id=4
      DO k=kbgn,kend
        DO i=ibgn,iend
          ik = i-ibgn+1 + (k-kbgn)*isize
          tem1(ik) = -9999.0
          tem2(ik) = -9999.0
          IF(u2(i,k) /= -9999.0) tem1(ik)= u2(i,k)*factor
          IF(v2(i,k) /= -9999.0) tem2(ik)= v2(i,k)*factor
          tem3(ik)=xw+(i-ibgn)* sqrtdxy
          tem4(ik)=zp(i,k)
        END DO
      END DO
    ELSE
      id=2
      DO k=kbgn,kend
        DO i=ibgn,iend
          ik = i-ibgn+1 + (k-kbgn)*isize
          tem1(ik) = -9999.0
          tem2(ik) = -9999.0
          IF(u2(i,k) /= -9999.0 .AND. v2(i,k) /= -9999.0)               &
               tem1(ik)=(u2(i,k)*cosaf+v2(i,k)*sinaf)*factor
          IF(w2(i,k) /= -9999.0) tem2(ik)=w2(i,k)*factor
          tem3(ik)=xw+(i-ibgn)* sqrtdxy
          tem4(ik)=zp(i,k)
        END DO
      END DO
    END IF

    IF(xfmat == -1) THEN
      length=20
      CALL strmin(distc,length)
      IF(varname(1:6) == 'xuvplt') THEN
        length=20
        CALL strmin(distc,length)
        WRITE(title,'(''U-V '',A)') label
        WRITE(levlab,'(''XY-Z PLOT FROM '',4(A,F5.1),A,A)')             &
            '(',x101,',',y101,') through (',x102,',',y102,') ',         &
            distc(1:length)
      ELSE
        WRITE(title,'(''V-W '',A)') label
        WRITE(levlab,                                                   &
            '(''VERTICAL PLANE FROM '',4(A,F8.1),A,A)')                 &
            '(',x101,',',y101,') through (',x102,',',y102,') ',         &
            distc(1:length)
      END IF
    ELSE IF(xfmat == 0) THEN
      length=20
      CALL strmin(distc,length)
      IF(varname(1:6) == 'xuvplt') THEN
        WRITE(title,'(''U-V '',A)') label
        WRITE(levlab,'(''XY-Z PLOT FROM '',4(A,F5.1),A,A)')             &
            '(',x101,',',y101,') through (',x102,',',y102,') ',         &
            distc(1:length)
      ELSE
        WRITE(title,'(''V-W '',A)') label
        WRITE(levlab,                                                   &
            '(''VERTICAL PLANE FROM '',4(A,I5),A,A)')                   &
            '(',INT(x101),',',INT(y101),') through (',INT(x102),',',    &
            INT(y102),') ', distc(1:length)
      END IF
    ELSE
      WRITE(stem1,'(i1)')xfmat
      WRITE(stem2,43)stem1
      43       FORMAT('f8.',a1)
    END IF

    length = 120
    CALL strlnth( title, length )
    CALL strmin ( title, length)

    DO i=1,smooth
      CALL smooth9pmv(tem1,isize,ksize,1,isize,1,ksize,tem6)
      CALL smooth9pmv(tem2,isize,ksize,1,isize,1,ksize,tem6)
    END DO

    uunit = 10.0
    CALL xvmode(1)
    istep = ist
    jstep = kst
    CALL vtr2d(tem1,tem2,tem3,tem4,uunit, xw,xe,sqrtdxy,zb,zt,dz,       &
         isize,istep,ksize,jstep,                                       &
         title(1:length),runname, id,                                   &
         tem5,slicopt)

  END IF

  CALL xpscmnt('End plotting '//label_copy(1:llabel))
  RETURN
END SUBROUTINE vtr3d
!
!##################################################################
!##################################################################
!######                                                      ######
!######                SUBROUTINE VTR2D                      ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE vtr2d(u,v,x,y,uunit1, xl,xr,dx,yb,yt,dy,                     & 6,20
           m,istep,n,jstep,char1,char2, vpltmod,                        &
           hterain,slicopt)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!    Plot 2-d wind (u1,u2) vector field defined on grid points (x,y)
!    using ZXPLOT package..
!
!-----------------------------------------------------------------------
!
!  AUTHOR:
!
!  MODIFICATION HISTORY:
!
!  1/24/96 (J. Zong and M. Xue)
!  Fixed a problem related to finding the minima and maxima of u & v
!  when there exist missing data. The initial min. and max. should be
!  set to values other than the missing value, -9999.0.
!
!-----------------------------------------------------------------------
!
!
!  INPUT:
!    u        m by n 2-dimensional array of u (left-to-right)
!               wind components (m/s)
!    v        m by n 2-dimensional array of v (down-to-up)
!               wind components (m/s)
!
!    x        x coordinate of grid points in physical/comp. space (m)
!    y        y coordinate of grid points in physical/comp. space (m)
!
!    uunit1
!
!    xl,xr    The left and right bound of the physical domain.
!    dx       Spacing between the x-axis tick marks
!    yb,yt    Bottom and top bound of the physical domain.
!    dy       Spacing between the y-axis tick marks
!
!    m        First dimension of vector component array
!    istep    Step increment for plotting in x direction
!
!    n        Second dimension of vector component array
!    jstep    Step increment for plotting in y direction
!
!    char1    First character string to plot (title)
!    char2    Second character string to plot (runname)
!
!    vpltmod  vpltmod = 1 for u-v vector   (u=u, v=v in model space)
!             vpltmod = 2 for u-w vector   (u=u, v=w in model space)
!             vpltmod = 3 for v-w vector   (u=v, v=w in model space)
!    hterain  the height of terrain
!    slicopt  slice orientation indicator
!             = 1, x-y slice of at k=kslice is plotted.
!             = 2, x-z slice of at j=jslice is plotted.
!             = 3, y-z slice of at i=islice is plotted.
!             = 4, horizontal slice at z index islice is plotted.
!             = 5, xy-z cross section of wind islice is plotted.
!             = 6, data field on constant p-level is plotted.
!             = 0, all of the three slices above are plotted.
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
  INCLUDE 'arpsplt.inc'
!
  INTEGER :: m,n
!
  REAL :: u(m,n)
  REAL :: v(m,n)
  REAL :: x(m,n)
  REAL :: y(m,n)

  REAL :: uunit1
  REAL :: xl,xr,dx,yb,yt,dy
  INTEGER :: istep,jstep
!
  CHARACTER (LEN=*) :: char2
  CHARACTER (LEN=*) :: char1
!
  INTEGER :: slicopt,vpltmod
!
!-----------------------------------------------------------------------
!
!  Plotting control common blocks
!
!-----------------------------------------------------------------------
!
  INTEGER :: layover
  COMMON /laypar/ layover

  REAL :: ctinc,ctmin,ctmax,vtunt    !contour interval and vector unit
  COMMON /incunt/ ctinc,ctmin,ctmax,vtunt

  REAL :: xleng,vunit
  COMMON /vecscl/ xleng,vunit
!
  REAL :: yxratio                  !the scaling factor the y/x ratio.
  COMMON /yratio/ yxratio
!
  INTEGER :: iunits, itype
  COMMON /windvtr/iunits, itype
!
  INTEGER :: ovrstaopt
  INTEGER :: ovrstam,staset,ovrstan,ovrstav,stacol,markprio,wrtstax
  INTEGER :: nsta_typ,sta_typ(10),sta_marktyp(10),sta_markcol(10)
  REAL :: sta_marksz(10),wrtstad
  COMMON /sta_par/ ovrstaopt,ovrstam,staset,ovrstan,ovrstav,stacol,     &
         markprio,nsta_typ,sta_typ,sta_marktyp,                         &
         sta_markcol,sta_marksz,wrtstax,wrtstad
!
  INTEGER :: icolor,icolor1,lbcolor,trcolor                ! required color
  COMMON /recolor/icolor,icolor1,lbcolor,trcolor

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

  REAL :: lblmag, ctrlbsiz, axlbsiz
  COMMON /labmag/ lblmag, ctrlbsiz, axlbsiz
!
  INTEGER :: flag
  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
!
  REAL :: ubarb(200,200), vbarb(200,200)
  COMMON /windtmp/ubarb, vbarb
!
!
!-----------------------------------------------------------------------
!
!  Misc. local variables
!
!-----------------------------------------------------------------------
!
  INTEGER :: i,j,key
  REAL :: pl,pr,pb,pt   ! plot space left, right, bottom, top coordinate
  REAL :: px,py         ! plot space left-right length and up-down height
  REAL :: xs,ys         ! real space left-right length and up-down height
  REAL :: pxc,pyc       ! plot space left-right center and
!                                   up-down    center
  REAL :: x0,y0
  REAL :: umax,umin     ! max and min of u component
  REAL :: vmax,vmin     ! max and min of v component
  REAL :: uunit
  REAL :: am
!
  REAL :: zlevel
  COMMON/sliceh/zlevel

  INTEGER, ALLOCATABLE :: iwrk(:)
  REAL   , ALLOCATABLE :: xwk(:),ywk(:)
  INTEGER :: istatus
!
  REAL :: hterain(m,n)             ! The height of the terrain.

  INTEGER :: timeovr
  COMMON /timover/ timeovr

!
  INTEGER :: LEN, len1
  REAL :: xleng0,istand
  INTEGER :: iunits0
  CHARACTER (LEN=15) :: ichar2

  CHARACTER (LEN=150) :: f_char1
  CHARACTER (LEN=150) :: ch

  INTEGER :: ntitle,titcol, nxpic, nypic, wpltime
  REAL :: titsiz
  CHARACTER (LEN=132) :: ptitle(3), footer_l, footer_c, footer_r

  COMMON /titpar1/ptitle, footer_l, footer_c, footer_r
  COMMON /titpar2/ntitle,titcol,wpltime, nxpic, nypic
  COMMON /titpar3/titsiz

  INTEGER :: col_table,pcolbar
  COMMON /coltable/col_table,pcolbar
!
  CHARACTER (LEN=12) :: varname
  COMMON /varplt1/ varname

  REAL :: xw1,xe1,ys1,yn1
  COMMON /xuvpar/xw1,xe1,ys1,yn1

  INTEGER :: wrtflag
  CHARACTER (LEN=80) :: levlab
  CHARACTER (LEN=50) :: timelab
  CHARACTER (LEN=25) :: timestring
  COMMON /timelev/wrtflag, timelab, levlab, timestring

  CHARACTER (LEN=80) :: prestr
  INTEGER :: preflag
  COMMON /preinfo/ prestr,preflag

  REAL :: x101, y101, x102,y102
  COMMON /slicev1/x101, y101, x102,y102

  REAL :: ytmp   !!local temporary variable
  !wdt update
  REAL :: f_cputime,cpu1,cpu2
  DOUBLE PRECISION :: f_walltime,second1,second2

  INTEGER :: xnwpic_called
  COMMON /callnwpic/xnwpic_called

!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  ALLOCATE(iwrk(m*n),STAT=istatus) 
  ALLOCATE(xwk(m*n),ywk(m*n),STAT=istatus) 

  second1= f_walltime()
  cpu1 = f_cputime()

  WRITE(6,'(/1x,a,a)') ' Plotting ',char1

  IF( layover == 0 .OR. xnwpic_called == 0) THEN
    CALL xnwpic
    xnwpic_called=1
    timeovr=0
    wrtflag = 0
    preflag = 0
    prestr = levlab
    len1=80
    CALL strmin(prestr,len1)
  ELSE
    timeovr=1
    wrtflag = wrtflag + 1
  END IF
!
!-----------------------------------------------------------------------
!
!  Get plotting space variables
!
!-----------------------------------------------------------------------
!
  CALL xqpspc( pl, pr, pb, pt)
  px = pr - pl
  py = pt - pb

  xs = xr-xl
  ys = yt-yb

  pxc = (pr+pl)/2
  pyc = (pb+pt)/2
!
!-----------------------------------------------------------------------
!
!  Let the longest lenth determine size scaling of plot
!
!-----------------------------------------------------------------------
!
  IF( py/px >= ys*yxratio/xs ) THEN
    py = ys*yxratio/xs*px
    CALL xpspac(pl, pr, pyc-py/2, pyc+py/2 )
  ELSE
    px = xs/(ys*yxratio)*py
    CALL xpspac(pxc-px/2, pxc+px/2, pb, pt)
  END IF
!
!-----------------------------------------------------------------------
!
!  Set the real distance to plot distance scaling
!
!-----------------------------------------------------------------------
!
  CALL xmap( xl, xr, yb,yt)
!
!-----------------------------------------------------------------------
!
!  Plot maps, boxes, and polygons
!
!-----------------------------------------------------------------------
!
  CALL xcolor(lbcolor)

  CALL pltextra(slicopt, 1 )
!
!-----------------------------------------------------------------------
!
!  Find max and min of data array
!
!-----------------------------------------------------------------------
!
  DO j=1,n
    DO i=1,m
      IF(      umin = u(i,j)
      vmin = v(i,j)
      GO TO 110
    END DO
  END DO
  110   CONTINUE

  umax=umin
  vmax=vmin

  DO j=1,n
    DO i=1,m
      IF(      IF(      IF(      IF(    END DO
  END DO
!
!-----------------------------------------------------------------------
!
!  Fill various labels
!
!-----------------------------------------------------------------------
!
  CALL xchsiz( 0.030*(yt-yb)  * lblmag )
  CALL xcolor(lbcolor)

  IF ( layover < 1 ) THEN
    len1=50
    CALL strmin(timelab,len1)
    CALL xcharl(xl,yt+0.07*ys, timestring(1:25))
    CALL xcharr(xr+0.05*(xr-xl),yt+0.07*ys, timelab(1:len1))
    IF(levlab /= ' ') THEN
      len1=80
      CALL strmin(levlab,len1)
      CALL xcharc(xl+xs*0.5,yt+0.015*ys, levlab(1:len1))
      preflag = 1
    END IF
!     len1=80
!     CALL strmin(levlab,len1)
!     CALL xcharc(xl+xs*0.5,yt+0.015*ys, levlab(1:len1))
  END IF
  IF(preflag == 0 .AND. levlab /= ' ') THEN
    len1=80
    CALL strmin(levlab,len1)
    CALL xcharc(xl+xs*0.5,yt+0.015*ys, levlab(1:len1))
    preflag = 1
  END IF

  IF( vpltmod == 1 .OR. vpltmod == 4 ) THEN
    WRITE(ch,'(''Umin='',F7.2,'' Umax='',F7.2,                          &
    &        '' Vmin='',f7.2,'' Vmax='',f7.2)') umin,umax,vmin,vmax
  ELSE IF( vpltmod == 2 ) THEN
    WRITE(ch,'(''Umin='',F7.2,'' Umax='',F7.2,                          &
    &        '' Wmin='',f7.2,'' Wmax='',f7.2)') umin,umax,vmin,vmax
  ELSE
    WRITE(ch,'(''Vmin='',F7.2,'' Vmax='',F7.2,                          &
    &        '' Wmin='',f7.2,'' Wmax='',f7.2)') umin,umax,vmin,vmax
  END IF

  LEN=150
  CALL strmin(char1,LEN)
  IF( char1(LEN-1:LEN-1) == ')' ) char1(LEN-1:LEN-1)=','
  IF(itype == 1) THEN
    WRITE(f_char1, '(a,''VECTOR)'')')char1(1:LEN)
  ELSE IF(itype == 2) THEN
    WRITE(f_char1, '(a, ''BARB)'')')char1(1:LEN)
  END IF

! if first levlab is not equal second levlab then attatch levlab on f_ch
!
!mx

  LEN=150
  CALL strmin(f_char1,LEN)
  len1=80
  CALL strmin(levlab,len1)
!  IF (preflag.eq.1 .and. prestr(1:len1).ne.levlab(1:len1)
!    :   .and. prestr(1:1).ne.' '
!    :   .and.layover.ne.0 .and. levlab(1:1).ne.' ') THEN
!    write(f_char1,'(a,a)') f_char1(1:len),levlab(1:len1)
!  ENDIF

  WRITE(6,'(1x,a51)') ch(1:51)
  CALL xcolor(icolor)

  IF(lbaxis == 1) THEN
    IF(wrtstax == 0) THEN
      ytmp = 0.08
    ELSE
      ytmp =0.14
    END IF
  ELSE
    ytmp = 0.12
  END IF
  LEN=150
  CALL strmin(f_char1,LEN)

  CALL xchsiz(0.025*ys * lblmag )
  CALL xcharl(xl-0.20*(xr-xl), yb-(yt-yb)*(ytmp+wrtflag*0.030),         &
           f_char1(1:LEN))

  len1=150
  CALL strmin(ch,len1)
  CALL xcharr(xr+0.20*(xr-xl), yb-(yt-yb)*(ytmp+wrtflag*0.030),         &
           ch(1:len1))
!
!-----------------------------------------------------------------------
!
!  Set vector unit and plot vectors.
!
!-----------------------------------------------------------------------
!
  uunit=uunit1
  IF( vtunt /= 0.0 ) THEN
    uunit=vtunt
    CALL xvmode(2)
  END IF

! Set parameter for barb

  xleng0 = (pr-pl)/(m-1) * istep * 0.65
  IF(iunits == 1 .AND. itype == 2) THEN
    iunits0=1
    istand = 5.
    WRITE(ichar2,'(a15)')'5 m/s'
  ELSE IF(iunits == 2 .AND. itype == 2) THEN
    iunits0=2
    istand = 10.
    WRITE(ichar2,'(a15)')'10 knots'
  ELSE IF (iunits == 3 .AND. itype == 2) THEN
    iunits0=2
    istand = 10.
    WRITE(ichar2,'(a15)')'10 MPH'
  END IF

  IF(layover >= 1) CALL xcolor(icolor)

  CALL xmap(xl,xr, yb,yt)
  CALL xvectu(u,v,m,m,istep,n,jstep,xleng,uunit)

  CALL xcolor(icolor)
  CALL xwindw(xl, xr, yb, yt)

  IF(itype == 1) THEN
    CALL xvectr(u,v,x,y,m,m,istep,n,jstep,xleng,uunit)
  ELSE IF(itype == 2) THEN
    CALL xbarbs(u,v,x,y,m,m,istep,n,jstep,iunits0,xleng*0.65,2)
  END IF

  CALL xwdwof
!
!-----------------------------------------------------------------------
!
!  Plot axes with tick marks
!
!-----------------------------------------------------------------------
!
  CALL pltaxes(slicopt,dx,dy)

  vunit=uunit
  x0=xl-(xr-xl)*0.08
  y0=yb-(yt-yb)*0.07
  key=0
  am=0.5
  IF( ((m-1)/istep) > 30 ) am=1.0
  IF(itype == 1) THEN
    IF(varname(1:6) == 'xuvplt') CALL xmap(xl, xr,yb,yt)
    CALL xvectk(x0,y0,xleng*am,uunit*am, key)
    CALL xmap(xl, xr, yb, yt)
  END IF
!
!-----------------------------------------------------------------------
!
!  Plot terrain profile in vertical slices
!
!-----------------------------------------------------------------------
!
  IF(slicopt == 2 .OR. slicopt == 3 .OR.slicopt == 5) THEN
    CALL xcolor(trcolor)
    CALL xthick(2)
    CALL xpenup( x(1,1), y(1,1)-0.5*(y(1,2)-y(1,1)) )
    DO i=2,m
      CALL xpendn( x(i,1), y(i,1)-0.5*(y(i,2)-y(i,1)) )
    END DO
    CALL xthick(1)
  END IF
!
!-----------------------------------------------------------------------
!
!  Overlay terrain contour if required in x-y level
!      or Plot terrain outline in this slice zlevel .
!
!-----------------------------------------------------------------------
!
  IF(timeovr == 0) CALL plttrn(hterain,x,y,m,n,slicopt)

  CALL xcolor(lbcolor)

  CALL xwindw(xl, xr, yb, yt)

!
!-----------------------------------------------------------------------
!
!  Plot observations
!
!-----------------------------------------------------------------------
!
  IF(ovrobs == 1 .AND. obsset == 1.AND. (slicopt == 1.OR.slicopt == 4.OR. &
        slicopt == 6.OR.slicopt == 7.OR.slicopt == 8)) THEN
    CALL pltobs(3)
    obsset=0
  END IF
!
!-----------------------------------------------------------------------
!
!  Plot station labels
!
!-----------------------------------------------------------------------
!
  CALL xcolor(lbcolor)
  IF(ovrstaopt == 1 .AND. staset == 1 .AND.                             &
        (ovrstam == 1.OR.ovrstan == 1.OR.ovrstav == 1).AND.             &
        (slicopt == 1.OR.slicopt == 4.OR.slicopt == 6                   &
        .OR.slicopt == 7.OR.slicopt == 8)                               &
        .AND.timeovr == 0  ) THEN
    CALL xchsiz(0.025*ys * lblmag)
    CALL pltsta(u,v,x,y,m,n,0,slicopt)
!    staset=0
  END IF
  IF( ovrstaopt == 1 .AND. wrtstax == 1 .AND. timeovr == 0              &
        .AND.(slicopt == 2.OR.slicopt == 3.OR. slicopt == 5) ) THEN
    CALL xchsiz(0.025*ys * lblmag)
    flag=1
    CALL pltsta(u,v,x,y,m,n,flag,slicopt)
  END IF

  CALL xwdwof

!-----------------------------------------------------------------------
!
!  Plot additional text below the figure
!
!-----------------------------------------------------------------------

  CALL label2d(char2)

  cpu2 = f_cputime()
  second2 = f_walltime()
!  write(6,*)'!!!!  total cpu time for one VTR2D  :',
!    :    cpu2-cpu1,' PLOT:',varname

  DEALLOCATE(iwrk)
  DEALLOCATE(xwk,ywk)

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


SUBROUTINE vtrunt( vtunt0 ) 6

!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!    Set the wind vector unit for wind field to be plotted by VTR2D.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Ming Xue
!
!  MODIFICATION HISTORY:
!    6/08/92  Added full documentation (K. Brewster)
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!    vtunt0    Unit vector
!              If VTUNT0 = 0.0, the unit is internally determined.
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
  REAL :: vtunt0
!
!-----------------------------------------------------------------------
!
!  Plotting control common blocks
!
!-----------------------------------------------------------------------
!
  REAL :: ctinc,ctmin,ctmax,vtunt   ! contour interval and vector unit
  COMMON /incunt/ ctinc,ctmin,ctmax,vtunt
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  vtunt = vtunt0

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


SUBROUTINE strm3d(u,v,w, x,y,z, xw,xe,dx, ys,yn,dy, zb,zt,dz,           & 2,43
           nx,ibgn,iend,ist, ny,jbgn,jend,jst, nz,kbgn,kend,kst,        &
           kslice, jslice, islice, time, runname,factor,slicopt,        &
           n,xp,yp,zp,u1,v1,u2,v2,w2,                                   &
           tem1,tem2,tem3,tem4,tem5,                                    &
           tem6,hterain)

!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!    Plot a streamline field in 2-d slices
!
!  AUTHOR: Ming Xue
!    1/16/1992
!
!  MODIFICATION HISTORY:
!
!  3/25/96 (K. Brewster)
!    Added variables isize,jsize,ksize
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!    u        3-dimensional array of u wind components (m/s)
!    v        3-dimensional array of v wind components (m/s)
!    w        3-dimensional array of w wind components (m/s)
!
!    x        x coordinate of grid points in physical/comp. space (m)
!    y        y coordinate of grid points in physical/comp. space (m)
!    z        z coordinate of grid points in physcal space (m)
!
!    xw       value of x for first i grid point to plot
!    xe       value of x for last i grid point to plot
!    ys       value of y for first j grid point to plot
!    yn       value of y for last j grid point to plot
!    zb       value of z for first k grid point to plot
!    zt       value of z for last k grid point to plot
!
!    nx       first dimension of b
!    ibgn     index of first i grid point to plot
!    iend     index of last  i grid point to plot
!
!    ny       second dimension of b
!    jbgn     index of first j grid point to plot
!    jend     index of last  j grid point to plot
!
!    nz       third dimension of b
!    kbgn     index of first k grid point to plot
!    kend     index of last  k grid point to plot
!
!    ist      step size in x direction
!    jst      step size in y direction
!    kst      step size in z direction
!
!    time     time of data in seconds
!
!    kslice   k index of plane for slicopt=1 x-y slice
!    jslice   j index of plane for slicopt=2 x-z slice
!    islice   i index of plane for slicopt=1 y-z slice
!
!    runname  character string decribing run
!
!    factor   scaling factor for winds
!             V*factor wind vectors are plotted
!
!    slicopt  slice orientation indicator
!             = 1, x-y slice of at k=kslice is plotted.
!             = 2, x-z slice of at j=jslice is plotted.
!             = 3, y-z slice of at i=islice is plotted.
!             = 4, horizontal slice at z index islice is plotted.
!             = 5, xy-z cross section of wind islice is plotted.
!             = 6, data field on constant p-level is plotted.
!             = 0, all of the three slices above are plotted.
!
!  WORK ARRAYS:
!
!    tem1     Temporary work array.
!    tem2     Temporary work array.
!    tem3     Temporary work array.
!    tem4     Temporary work array.
!    tem5     Temporary work array.
!    tem6     Temporary work array.
!
!  (These arrays are defined and used locally (i.e. inside this
!   subroutine), they may also be passed into routines called by
!   this one. Exiting the call to this subroutine, these temporary
!   work arrays may be used for other purposes therefore their
!   contents overwritten. Please examine the usage of work arrays
!   before you alter the code.)
!
!   pp01      The pressure (mb) value at the specific p-level
!   ercpl     reciprocal of exponent
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
!
  INTEGER :: nx,ny,nz
  INTEGER :: n
!
  REAL :: u(nx,ny,nz)
  REAL :: v(nx,ny,nz)
  REAL :: w(nx,ny,nz)

  REAL :: x(nx,ny,nz)
  REAL :: y(nx,ny,nz)
  REAL :: z(nx,ny,nz)
!
  REAL :: u1(nx,ny),v1(nx,ny)
  REAL :: u2(n,nz),v2(n,nz),w2(n,nz),zp(n,nz)
  REAL :: xp(n),yp(n)

  INTEGER :: kslice,jslice,islice
  CHARACTER (LEN=*) :: runname

  REAL :: xw,xe,dx,ys,yn,dy,zb,zt,dz
  INTEGER :: ibgn,iend,ist, jbgn,jend,jst, kbgn,kend,kst

  REAL :: time,factor
  INTEGER :: slicopt

  REAL :: x_tmp
  COMMON /tmphc2/ x_tmp
!
!-----------------------------------------------------------------------
!
!  Some constants
!
!-----------------------------------------------------------------------
!
  REAL :: pp01, ercpl
  PARAMETER (ercpl=0.3678794)              ! exp(-1.0)
!
!-----------------------------------------------------------------------
!
!  Work arrays: tem1,tem2,tem3 of size at least
!               max( nx*ny, nx*nz, ny*nz).
!
!-----------------------------------------------------------------------
!
  REAL :: tem1(*),tem2(*),tem3(*),tem4(*),tem5(*),tem6(*)

  INTEGER :: nzmax
  PARAMETER (nzmax = 300)
  REAL :: fdata(nzmax),zdata(nzmax),fprof(nzmax),zprof(nzmax)
!
!-----------------------------------------------------------------------
!
!  Common blocks for plotting control parameters
!
!-----------------------------------------------------------------------
!
  REAL :: x01,y01                  ! the first  point of interpolation
  REAL :: x02,y02                  ! the second point of interpolation
  REAL :: zlevel                   ! the given height of the slice
  REAL :: sinaf,cosaf,dist,sqrtdxy
  COMMON /slicev/x01,y01,x02,y02,sinaf,cosaf,dist,sqrtdxy
  COMMON /sliceh/zlevel
!
!-----------------------------------------------------------------------
!
!  Misc. local Variables
!
!-----------------------------------------------------------------------
!
  INTEGER :: i,j,k,ij,ik,jk,length,isize,jsize,ksize
  CHARACTER (LEN=6) :: timhms
  CHARACTER (LEN=120) :: title
!
  INTEGER :: icolor,icolor1,lbcolor,trcolor                ! required color
  COMMON /recolor/icolor,icolor1,lbcolor,trcolor
!
  INTEGER :: trnplt                ! flag to plot terain (1 or 0)
  REAL :: hterain(nx,ny)           ! The height of the terrain.

  INTEGER :: ovrtrn                ! overlay terrain option (0/1)
  REAL :: trninc,trnmin, trnmax    ! terrain interval minimum, maximum
  REAL :: ztmin,ztmax
  COMMON /trnpar/ trnplt,ovrtrn,trninc,trnmin, trnmax,ztmin,ztmax
!
  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
  CHARACTER (LEN=4) :: stem2
  CHARACTER (LEN=1) :: stem1

  INTEGER :: smooth
  COMMON /smoothopt/smooth

  INTEGER :: wrtflag
  CHARACTER (LEN=80) :: levlab
  CHARACTER (LEN=50) :: timelab
  CHARACTER (LEN=25) :: timestring
  COMMON /timelev/wrtflag, timelab, levlab, timestring

  REAL :: tmpx, tmpy
  CHARACTER (LEN=20) :: distc
  REAL :: x101, y101, x102,y102
  COMMON /slicev1/x101, y101, x102,y102

!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  isize=(iend-ibgn)+1
  jsize=(jend-jbgn)+1
  ksize=(kend-kbgn)+1
!
!-----------------------------------------------------------------------
!
!  setup  time label
!
!-----------------------------------------------------------------------
!
  CALL cvttim ( time, timhms)
  IF( timhms(1:1) == '0' ) timhms(1:1)=' '
  WRITE(timelab,'(''T='',F8.1,A)') time,                                &
      ' s ('//timhms(1:2)//':'//timhms(3:4)//':'//timhms(5:6)//')'
  CALL get_time_string ( time, timestring)
!
!-----------------------------------------------------------------------
!
!  Set up terrain, if needed.
!
!-----------------------------------------------------------------------
!
  IF(trnplt == 1 .OR.trnplt == 2 .OR. ovrtrn == 1)  THEN
    DO j=jbgn,jend
      DO i=ibgn,iend
        ij = i-ibgn+1 + (j-jbgn)*isize
        tem6(ij)=hterain(i,j)
      END DO
    END DO
  END IF

  IF ( slicopt == 2 .OR. slicopt == 3  .OR. slicopt == 5) THEN
    CALL cal_dist(haxisu,dx,dy,x01,y01,x02,y02,slicopt,                 &
                  tmpx,tmpy,distc)
  END IF

!
!-----------------------------------------------------------------------
!
!  slicopt=1   Plot u-v field
!
!-----------------------------------------------------------------------
!
  IF( slicopt == 1 .OR. slicopt == 0 ) THEN

    k = kslice
    DO j=jbgn,jend
      DO i=ibgn,iend
        ij = i-ibgn+1 + (j-jbgn)*isize
        tem1(ij) = -9999.0
        tem2(ij) = -9999.0
        IF(        IF(        tem3(ij)=x(i,j,k)
        tem5(ij)=y(i,j,k)
      END DO
    END DO

    IF (k /= 2) THEN
      WRITE(title,'(''U-V STREAMLINE'')')
      WRITE(levlab,'(''X-Y CROSS SECTION THROUGH K='',I3)')k
    ELSE
      WRITE(title,'(''U-V STREAMLINE'')')
      WRITE(levlab,'(''X-Y CROSS SECTION THROUGH K=2 (SURFACE)'')')
    END IF

    length = 120
    CALL strlnth( title, length )
    CALL strmin ( title, length)

    DO i=1,smooth
      CALL smooth9pmv(tem1,isize,jsize,1,isize,1,jsize,tem4)
      CALL smooth9pmv(tem2,isize,jsize,1,isize,1,jsize,tem4)
    END DO

    CALL strm2d(tem1,tem2, xw,xe,ys,yn, dx, dy,                         &
         isize,jsize,                                                   &
         title(1:length),runname, tem3,tem5,                            &
         tem6,slicopt)

!
!-----------------------------------------------------------------------
!
!  slicopt=2   Plot u-w streamline
!
!-----------------------------------------------------------------------
!
  ELSE IF( slicopt == 2 .OR. slicopt == 0 ) THEN

    x_tmp = y(1,jslice,1)

    j = jslice
    DO k=kbgn,kend
      DO i=ibgn,iend
        ik = i-ibgn+1 + (k-kbgn)*isize
        tem1(ik) = -9999.0
        tem2(ik) = -9999.0
        IF(        IF(        tem4(ik)=z(i,j,k)
        tem3(ik)=x(i,j,k)
      END DO
    END DO

    IF( nzmax < ksize) THEN
      WRITE(6,'(1x,a)')                                                 &
          'nzmax given in STRM3D too small. Job stopped.'
      STOP
    END IF

    DO k=1,ksize
      zprof(k)= zb+(k-1)*(zt-zb)/(kend-kbgn)
    END DO

    CALL unigrid(isize,ksize,tem1,tem4,                                 &
         fdata,zdata,fprof,zprof)
    CALL unigrid(isize,ksize,tem2,tem4,                                 &
         fdata,zdata,fprof,zprof)

    WRITE(title,'(''U-W STREAMLINE '')')
    dist = (j-1)*tmpy
    length=20
    CALL strmin ( distc, length)
    WRITE(levlab,'(''X-Z CROSS SECTION THROUGH J='',I3,                 &
    &   '' (y = '',f8.1,a)')j,dist,distc(1:length)

    length = 120
    CALL strlnth( title, length )
    CALL strmin ( title, length)


    DO i=1,smooth
      CALL smooth9pmv(tem1,isize,ksize,1,isize,1,ksize,tem4)
      CALL smooth9pmv(tem2,isize,ksize,1,isize,1,ksize,tem4)
    END DO

    CALL strm2d(tem1,tem2, xw,xe,zb,zt, dx, dz,                         &
         isize,ksize,                                                   &
         title(1:length),runname, tem3 ,tem4,                           &
         tem6,slicopt)

!
!-----------------------------------------------------------------------
!
!  slicopt=3   Plot v-w field
!
!-----------------------------------------------------------------------
!
  ELSE IF( slicopt == 3 .OR. slicopt == 0 ) THEN

    x_tmp = y(1,jslice,1)

    i = islice
    DO k=kbgn,kend
      DO j=jbgn,jend
        jk = j-jbgn+1 + (k-kbgn)*jsize
        tem1(jk) = -9999.0
        tem2(jk) = -9999.0
        IF(        IF(        tem4(jk)=z(i,j,k)
        tem5(jk)=y(i,j,k)
      END DO
    END DO

    IF( nzmax < ksize) THEN
      WRITE(6,'(1x,a)')                                                 &
          'nzmax given in STRM3D too small. Job stopped.'
      STOP
    END IF

    DO k=1,ksize
      zprof(k)= zb+(k-1)*(zt-zb)/(kend-kbgn)
    END DO

    CALL unigrid(jsize,ksize,tem1,tem4,                                 &
         fdata,zdata,fprof,zprof)
    CALL unigrid(jsize,ksize,tem2,tem4,                                 &
         fdata,zdata,fprof,zprof)

    dist = (i-1)*tmpx
    length=20
    CALL strmin ( distc, length)
    WRITE(levlab,'(''Y-Z CROSS SECTION THROUGH I='',I3,                 &
    &   '' ( x='',f8.1,a )')i,dist, distc(1:length)
    WRITE(title,'(''V-W STREAMLINE'')')

    length = 120
    CALL strlnth( title, length )
    CALL strmin ( title, length)

    DO i=1,smooth
      CALL smooth9pmv(tem1,jsize,ksize,1,jsize,1,ksize,tem4)
      CALL smooth9pmv(tem2,jsize,ksize,1,jsize,1,ksize,tem4)
    END DO

    CALL strm2d(tem1,tem2, ys,yn,zb,zt, dy, dz,                         &
         jsize,ksize,                                                   &
         title(1:length),runname, tem5 ,tem4,                           &
         tem6,slicopt)

!
!-----------------------------------------------------------------------
!
!  slicopt=4   Plot u-v streamlines on constant z levels
!  slicopt=6   Plot u-v streamlines on constant pressure levels
!  slicopt=7   Plot u-v streamlines on constant PT levels
!
!-----------------------------------------------------------------------
!
  ELSE IF( slicopt == 4.OR.slicopt == 6.OR.slicopt == 7 ) THEN

!    CALL hintrp(nx,ny,nz,u,z,zlevel,u1)
!    CALL hintrp(nx,ny,nz,v,z,zlevel,v1)

    CALL hintrp1(nx,ny,nz,kbgn,kend,u,z,zlevel,u1)
    CALL hintrp1(nx,ny,nz,kbgn,kend,v,z,zlevel,v1)

    DO j=jbgn,jend
      DO i=ibgn,iend
        ij = i-ibgn+1 + (j-jbgn)*isize
        tem1(ij) = -9999.0
        tem2(ij) = -9999.0
        IF(u1(i,j) /= -9999.0) tem1(ij)=u1(i,j)*factor
        IF(v1(i,j) /= -9999.0) tem2(ij)=v1(i,j)*factor
        tem3(ij)=x(i,j,2)
        tem5(ij)=y(i,j,2)
      END DO
    END DO

    IF( slicopt == 4) THEN
      WRITE(levlab,'(''Z='',F7.3,A,'' MSL'')')zlevel,' KM'
    ELSE IF( slicopt == 6) THEN
      pp01=0.01*ercpl**zlevel
      WRITE(levlab,'(''P='',F7.2,A)') pp01, ' MB'
    ELSE
      WRITE(levlab,'(''THETA='',F5.1,A)') zlevel, ' (K)'
    END IF

    WRITE(title,'(''U-V STREAMLINE'')')

    length = 120
    CALL strlnth( title, length )
    CALL strmin ( title, length)

    DO i=1,smooth
      CALL smooth9pmv(tem1,isize,jsize,1,isize,1,jsize,tem4)
      CALL smooth9pmv(tem2,isize,jsize,1,isize,1,jsize,tem4)
    END DO

    CALL strm2d(tem1,tem2, xw,xe,ys,yn, dx, dy,                         &
         isize,jsize,                                                   &
         title(1:length),runname, tem3 ,tem5,                           &
         tem6,slicopt)
!
!-----------------------------------------------------------------------
!
!  slicopt=5   Plot V-w field
!
!-----------------------------------------------------------------------
!
  ELSE IF( slicopt == 5 ) THEN

    CALL sectvrt(nx,ny,nz,u,x,y,z,u2,zp,n,xp,yp)
    CALL sectvrt(nx,ny,nz,v,x,y,z,v2,zp,n,xp,yp)
    CALL sectvrt(nx,ny,nz,w,x,y,z,w2,zp,n,xp,yp)

    DO k=kbgn,kend
      DO i=ibgn,iend
        ik = i-ibgn+1 + (k-kbgn)*isize
        tem1(ik) = -9999.0
        tem2(ik) = -9999.0
        IF(u2(i,k) /= -9999.0 .AND. v2(i,k) /= -9999.0)                 &
              tem1(ik)=(u2(i,k)*cosaf+v2(i,k)*sinaf)*factor
        IF(w2(i,k) /= -9999.0) tem2(ik)=w2(i,k)*factor
        tem3(ik)=xw+(i-ibgn)*sqrtdxy
        tem4(ik)=zp(i,k)
      END DO
    END DO

    IF( nzmax < ksize) THEN
      WRITE(6,'(1x,a)')                                                 &
          'nzmax given in STRM3D too small. Job stopped.'
      STOP
    END IF

    DO k=1,ksize
      zprof(k)= zb+(k-1)*(zt-zb)/(kend-kbgn)
    END DO

    CALL unigrid(isize,ksize,tem1,tem4,                                 &
         fdata,zdata,fprof,zprof)
    CALL unigrid(isize,ksize,tem2,tem4,                                 &
         fdata,zdata,fprof,zprof)

    IF(xfmat == -1) THEN
      length=20
      CALL strmin(distc,length)
      WRITE(title,'(''V-W STREAMLINE'')')
      WRITE(levlab,                                                     &
          '(''VERT CROSS SECTION THROUGH '',4(A,F8.1),A,A)')            &
          '(',x101,',',y101,') (',x102,',',y102,')',distc(1:length)
    ELSE IF(xfmat == 0 ) THEN
      length=20
      CALL strmin(distc,length)
      WRITE(title,'(''V-W STREAMLINE'')')
      WRITE(levlab,                                                     &
          '(''VERT CROSS SECTION THROUGH '',4(A,I5),A,A)')              &
          '(',INT(x101),',',INT(y101),') (',INT(x102),',',INT(y102),    &
          ')',distc(1:length)
    ELSE
      WRITE(stem1,'(i1)')xfmat
      WRITE(stem2,43)stem1
      43       FORMAT('f8.',a1)
!       write(title,'(''V-w streamline at t='',f8.1,a,
!  :       '' through '',4(a,stem2),a)')
!  :       time,' s ('//timhms(1:2)//':'//timhms(3:4)//':'//
!  :       timhms(5:6)//')','(',x01,',',y01,') (',x02,',',y02,')'
    END IF

    length = 120
    CALL strlnth( title, length )
    CALL strmin ( title, length)

    DO i=1,smooth
      CALL smooth9pmv(tem1,isize,ksize,1,isize,1,ksize,tem4)
      CALL smooth9pmv(tem2,isize,ksize,1,isize,1,ksize,tem4)
    END DO

    CALL strm2d(tem1,tem2, xw,xe,zb,zt, sqrtdxy, dz,                    &
         isize,ksize,                                                   &
         title(1:length),runname, tem3,tem4,                            &
         tem6,slicopt)

  END IF

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


SUBROUTINE strm2d(u,v,xl,xr,yb,yt,dx,dy,m,n,char1,char2, x,y,           & 5,6
           hterain,slicopt)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!    Plot streamlines of a 2-d wind (u1,u2) field using ncargraphic
!    subroutine strmln
!
!  INPUT:
!    u        m by n 2-dimensional array of u (left-to-right)
!               wind components (m/s)
!    v        m by n 2-dimensional array of v (down-to-up)
!               wind components (m/s)
!
!    xl,xr    The left and right bound of the physical domain.
!    yb,yt    Bottom and top bound of the physical domain.
!    dx,dy    Grid interval in x and y direction (km)
!
!    m        First dimension of vector component array
!    n        Second dimension of vector component array
!
!    char1    First character string to plot (title)
!    char2    Second character string to plot (runname)
!
!    x        x coordinate of grid points in plot space (over on page)
!    y        y coordinate of grid points in plot space (up on page)
!
!    hterain  the height of terrain
!    slicopt  slice orientation indicator
!       slicopt = 1, x-y slice of u,v at z index kslice is plotted.
!       slicopt = 2, x-z slice of u,w at y index jslice is plotted.
!       slicopt = 3, y-z slice of v,w at x index islice is plotted.
!       slicopt = 4, x-y slice of u,v at z index islice is plotted.
!       slicopt = 5, xy-z cross section of wind islice is plotted.
!       slicopt = 6, data field on constant p-level is plotted.
!       slicopt = 0, all of the three slices above are plotted.
!
!  WORK ARRAY
!    tem1      A work array of size at least m*n*2
!    tem2      A work array of size at least m*n*2
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
  INCLUDE 'arpsplt.inc'
!
  INTEGER :: m,n
  INTEGER :: i
!
  REAL :: u(m,n)
  REAL :: v(m,n)

  REAL :: x(m,n)
  REAL :: y(m,n)
  REAL :: xl,xr,yb,yt,dx,dy
!
  CHARACTER (LEN=*) :: char2
  CHARACTER (LEN=*) :: char1
  INTEGER :: ierror
!
!-----------------------------------------------------------------------
!
!  Plotting control common blocks
!
!-----------------------------------------------------------------------
!
  INTEGER :: layover
  COMMON /laypar/ layover

  REAL :: ctinc,ctmin,ctmax,vtunt  !contour interval and vector unit
  COMMON /incunt/ ctinc,ctmin,ctmax,vtunt
!
  INTEGER :: icolor,icolor1,lbcolor,trcolor                ! required color
  COMMON /recolor/icolor,icolor1,lbcolor,trcolor
!
  INTEGER :: flag
  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
!
  INTEGER :: ovrstaopt
  INTEGER :: ovrstam,staset,ovrstan,ovrstav,stacol,markprio,wrtstax
  INTEGER :: nsta_typ,sta_typ(10),sta_marktyp(10),sta_markcol(10)
  REAL :: sta_marksz(10),wrtstad
  COMMON /sta_par/ ovrstaopt,ovrstam,staset,ovrstan,ovrstav,stacol,     &
         markprio,nsta_typ,sta_typ,sta_marktyp,                         &
         sta_markcol,sta_marksz,wrtstax,wrtstad
!
  REAL :: yxratio                  !the scaling factor the y/x ratio.
  COMMON /yratio/ yxratio

  INTEGER :: col_table,pcolbar
  COMMON /coltable/col_table,pcolbar
!
!-----------------------------------------------------------------------
!
!  Misc. local variables
!
!-----------------------------------------------------------------------
!
  INTEGER :: nopic,nhpic,nvpic,ifont
  REAL :: pl,pr,pb,pt   ! plot space left, right, bottom, top coordinate
  REAL :: px,py         ! plot space left-right length and up-down height
  REAL :: xs,ys         ! real space left-right length and up-down height
  REAL :: pxc,pyc       ! plot space left-right center and
!                                   up-down    center
  REAL :: xlimit,ylimit
  REAL :: rotang
  REAL :: xp1,xp2,yp1,yp2
  REAL :: xd1,xd2,yd1,yd2,xpos1,xpos2,ypos1,ypos2
!
  REAL :: zlevel
  COMMON/sliceh/zlevel
!
  REAL :: hterain(m,n)       ! The height of the terrain.
!
  INTEGER :: slicopt

  INTEGER :: timeovr
  COMMON /timover/ timeovr

  REAL :: lblmag, ctrlbsiz, axlbsiz
  COMMON /labmag/ lblmag, ctrlbsiz, axlbsiz
!
  INTEGER :: len1

  INTEGER :: wrtflag
  CHARACTER (LEN=80) :: levlab
  CHARACTER (LEN=50) :: timelab
  CHARACTER (LEN=25) :: timestring
  COMMON /timelev/wrtflag, timelab, levlab, timestring

  REAL :: x101, y101, x102,y102
  COMMON /slicev1/x101, y101, x102,y102

  INTEGER :: xnwpic_called
  COMMON /callnwpic/xnwpic_called

  REAL :: ytmp   !! local temporary variable
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  WRITE(6,'(/1x,a,a)') ' Plotting ',char1

  IF( layover == 0 .OR. xnwpic_called == 0) THEN
    CALL xnwpic
    xnwpic_called=1
    timeovr=0
    wrtflag = 0
  ELSE
    timeovr=1
    wrtflag = wrtflag + 1
  END IF
!
!-----------------------------------------------------------------------
!
!  Get plotting space variables
!
!-----------------------------------------------------------------------
!
  CALL xqpspc( pl, pr, pb, pt)
  px = pr - pl
  py = pt - pb

  xs = xr-xl
  ys = yt-yb

  pxc = (pr+pl)/2
  pyc = (pb+pt)/2
!
!-----------------------------------------------------------------------
!
!  Let the longest lenth determine size scaling of plot
!
!-----------------------------------------------------------------------
!
  IF( py/px >= ys*yxratio/xs ) THEN
    py = ys*yxratio/xs*px
    CALL xpspac(pl, pr, pyc-py/2, pyc+py/2 )
  ELSE
    px = xs/(ys*yxratio)*py
    CALL xpspac(pxc-px/2, pxc+px/2, pb, pt)
  END IF
!
!-----------------------------------------------------------------------
!
!  Set the real distance to plot distance scaling
!
!-----------------------------------------------------------------------
!
  CALL xmap( xl, xr, yb,yt)
!
!-----------------------------------------------------------------------
!
!  Plot map, boxes and polygons.
!
!-----------------------------------------------------------------------
!
  CALL xcolor(lbcolor)

  CALL pltextra(slicopt, 1 )

  xpos1 = xl
  xpos2 = xr
  ypos1 = yb
  ypos2 = yt

  CALL xtrans(xpos1,ypos1)
  CALL xtrans(xpos2,ypos2)
  CALL xzx2ncar(xpos1,ypos1)
  CALL xzx2ncar(xpos2,ypos2)
!
  IF(slicopt == 2 .OR. slicopt == 3 .OR.slicopt == 5) THEN
    CALL xcolor(trcolor)
    CALL xthick(3)
    CALL xpenup( x(1,1), y(1,1)-0.5*(y(1,2)-y(1,1)) )
    DO i=2,m
      CALL xpendn( x(i,1), y(i,1)-0.5*(y(i,2)-y(i,1)) )
    END DO
    CALL xthick(1)
  END IF
!
!-----------------------------------------------------------------------
!
!  Overlay terrain contour if required in x-y level
!      or Plot terrain outline in this slice zlevel .
!
!-----------------------------------------------------------------------
!
  IF( timeovr == 0)CALL plttrn(hterain,x,y,m,n,slicopt)

  CALL xcolor(lbcolor)
!
!-----------------------------------------------------------------------
!
!  Plot station labels
!
!-----------------------------------------------------------------------
!
  IF(ovrstaopt == 1 .AND. staset == 1 .AND.                             &
        (ovrstam == 1.OR.ovrstan == 1.OR.ovrstav == 1).AND.             &
        (slicopt == 1.OR.slicopt == 4.OR.slicopt == 6                   &
        .OR.slicopt == 7.OR.slicopt == 8)                               &
        .AND.timeovr == 0 ) THEN
    CALL xchsiz(0.025*ys * lblmag)
    CALL pltsta(u,v,x,y,m,n,0,slicopt)
!    staset=0
  END IF
!
!-----------------------------------------------------------------------
!
!  Plot observations
!
!-----------------------------------------------------------------------
!
  IF( ovrstaopt == 1 .AND. wrtstax == 1 .AND. timeovr == 0              &
        .AND.(slicopt == 2.OR.slicopt == 3.OR. slicopt == 5) ) THEN
    CALL xchsiz(0.025*ys * lblmag)
    flag=1
    CALL pltsta(u,v,x,y,m,n,flag,slicopt)
  END IF
!
!-----------------------------------------------------------------------
!
!  Plot streamlines
!
!-----------------------------------------------------------------------
!
  CALL xcolor(lbcolor)
!
  CALL xqset(xp1,xp2,yp1,yp2, xd1,xd2,yd1,yd2)
  CALL set(xpos1,xpos2,ypos1,ypos2, 1.0, FLOAT(m), 1.0, FLOAT(n),1)

  CALL xcolor(icolor)

  CALL strmln(u,v,y ,m,m,n, 1, ierror)
  CALL set(xp1,xp2,yp1,yp2, xd1,xd2,yd1,yd2, 1)
!
!-----------------------------------------------------------------------
!
!  Plot axes with tick marks
!
!-----------------------------------------------------------------------
!
  CALL pltaxes(slicopt,dx,dy)

!-----------------------------------------------------------------------
!
!  Plot labels
!
!-----------------------------------------------------------------------

  CALL xcolor(lbcolor)

  CALL xqnpic(nopic)
  CALL xqspac(nhpic, nvpic, rotang, xlimit, ylimit)

! write time and level
  CALL xchsiz( 0.030*ys  * lblmag )
  IF ( layover < 1 ) THEN
    len1=50
    CALL strmin(timelab,len1)
    CALL xcharl(xl,yt+0.07*ys, timestring(1:25))
    CALL xcharr(xr+0.05*(xr-xl),yt+0.07*ys, timelab(1:len1))
    len1=80
    CALL strmin(levlab,len1)
    CALL xcharc(xl+xs*0.5,yt+0.015*ys, levlab(1:len1))
  END IF

! write variable label
  CALL xcolor(icolor)
  IF(lbaxis == 1) THEN
    IF(wrtstax == 0) THEN
      ytmp = 0.08
    ELSE
      ytmp = 0.14
    END IF
  ELSE
    ytmp =0.12
  END IF

  CALL xchsiz( 0.025*ys  * lblmag )
  CALL xcharl(xl-0.20*(xr-xl), yb-(yt-yb)*(ytmp+layover*0.030),         &
       char1)

  CALL xcolor(lbcolor)
  IF (timeovr == 0) THEN
    IF(nopic == nhpic*(nvpic-1)+1 ) THEN
      ytmp =0.25
      IF(layover < 1) CALL xcharl(xl,yb-(ytmp+layover*0.03)*(yt-yb), char2 )

      CALL xqcfnt(ifont)
      CALL xcfont(xfont)
      ytmp = 0.20
      IF(layover < 1) CALL xcharl(xl,yb-(0.20+layover*0.03)*(yt-yb),    &
          'CAPS/ARPS ' )
!    :     'Project Hub-CAPS Experimental ' )         !Hub-CAPS+

      CALL xcfont(ifont)
    END IF
    timeovr=1
  END IF

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


SUBROUTINE ctrsfc(a,x,y,x1,x2,dx,y1,y2,dy,                              & 41,6
           nx,ibgn,iend, ny,jbgn,jend,                                  &
           label,time, runname, factor,tem1,tem2,tem3,                  &
           tem4,tem5,hterain,slicopt,pltopt)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!    To plot a contour map for a 2-d surface array.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Ming Xue
!    4/20/1994
!
!  MODIFICATION HISTORY:
!
!    9/27/95 (Yuhe Liu)
!    Fixed a bug in call of smth. Added the temporary array tem5 to
!    the argument list.
!
!    3/25/96 (Keith Brewster)
!    Added variables isize,jsize and replaced smth with smooth9p
!
!-----------------------------------------------------------------------
!
!  INPUT :
!
!    a        2-d surface array.
!
!    x        x coordinate of grid points in physical/comp. space (m)
!    y        y coordinate of grid points in physical/comp. space (m)
!
!    x1       value of x for first i grid point to plot
!    x2       value of x for last i grid point to plot
!    dx
!    y1       value of y for first j grid point to plot
!    y2       value of y for last j grid point to plot
!    dy
!
!    nx       first dimension of a
!    ibgn     index of first i grid point to plot
!    iend     index of last  i grid point to plot
!
!    ny       second dimension of a
!    jbgn     index of first j grid point to plot
!    jend     index of last  j grid point to plot
!
!    label    character string describing the contents of a
!
!    time     time of data in seconds
!
!    runname  character string decribing run
!
!    factor   scaling factor for data
!             contours are labelled a*factor
!    slicopt  slice orientation indicator
!       slicopt = 1, x-y slice of u,v at z index kslice is plotted.
!       slicopt = 2, x-z slice of u,w at y index jslice is plotted.
!       slicopt = 3, y-z slice of v,w at x index islice is plotted.
!       slicopt = 4, x-y slice of u,v at z index islice is plotted.
!       slicopt = 5, xy-z cross section of wind islice is plotted.
!       slicopt = 6, data field on constant p-level is plotted.
!       slicopt = 0, all of the three slices above are plotted.
!    plot     variable plot option (0/1/2/3)
!
!  WORK ARRAYS:
!
!    tem1     Temporary work array.
!    tem2     Temporary work array.
!    tem3     Temporary work array.
!    tem4     Temporary work array.
!    tem5     Temporary work array.
!
!
!    hterain  The height of the terrain.
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
!
  INTEGER :: nx,ny

  REAL :: a(nx,ny)
  REAL :: x(nx,ny)
  REAL :: y(nx,ny)

  REAL :: x1,x2,dx,y1,y2,dy
  INTEGER :: ibgn,iend,jbgn,jend,length

  CHARACTER (LEN=6) :: timhms
  CHARACTER (LEN=*) :: label
  CHARACTER (LEN=*) :: runname

  REAL :: time
  REAL :: factor

  REAL :: tem1(*)
  REAL :: tem2(*)
  REAL :: tem3(*)
  REAL :: tem4(*)
  REAL :: tem5(*)

  REAL :: hterain(nx,ny)

  INTEGER :: slicopt
  INTEGER :: pltopt       ! variable plot option (0/1/2/3)

  INTEGER :: ovrtrn,trnplt               ! overlay terrain option (0/1)
  REAL :: trninc,trnmin, trnmax   ! terrain interval minimum, maximum
  REAL :: ztmin,ztmax
  COMMON /trnpar/ trnplt,ovrtrn,trninc,trnmin, trnmax,ztmin,ztmax

  INTEGER :: smooth
  COMMON /smoothopt/smooth

  INTEGER :: wrtflag
  CHARACTER (LEN=120) :: label_copy
  CHARACTER (LEN=80) :: levlab
  CHARACTER (LEN=50) :: timelab
  CHARACTER (LEN=25) :: timestring
  COMMON /timelev/wrtflag, timelab, levlab, timestring

!
!-----------------------------------------------------------------------
!
!  Misc. local variables
!
!-----------------------------------------------------------------------
!
  INTEGER :: i,j,ij,isize,jsize,llabel
  CHARACTER (LEN=120) :: title
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  isize=(iend-ibgn)+1
  jsize=(jend-jbgn)+1

  label_copy = label
  llabel = 120
  CALL xstrlnth(label_copy, llabel)
  CALL xpscmnt('Start plotting '//label_copy(1:llabel))
!
!-----------------------------------------------------------------------
!
!  Set up terrain, if needed.
!
!-----------------------------------------------------------------------
!
  IF(ovrtrn == 1)  THEN
    DO j=jbgn,jend
      DO i=ibgn,iend
        ij = i-ibgn+1 + (j-jbgn)*isize
        tem4(ij)=hterain(i,j)
      END DO
    END DO
  END IF
!
  CALL cvttim( time, timhms )

  IF( timhms(1:1) == '0' ) timhms(1:1)=' '
  WRITE(timelab,'(''T='',F8.1,A)') time,                                &
      ' s ('//timhms(1:2)//':'//timhms(3:4)//':'//timhms(5:6)//')'
  CALL get_time_string ( time, timestring)

  DO j=jbgn,jend
    DO i=ibgn,iend
      ij = i-ibgn+1 + (j-jbgn)*isize
      tem1(ij) = -9999.0
      IF(      tem2(ij)=x(i,j)
      tem3(ij)=y(i,j)
    END DO
  END DO

  levlab=' '
  WRITE(title,'(a)') label

  length = 120
  CALL strlnth( title, length)
  CALL strmin ( title, length)

  DO i=1,smooth
    CALL smooth9pmv(tem1,isize,jsize,1,isize,1,jsize,tem5)
  END DO

  CALL ctr2d(tem1,tem2,tem3, x1,x2,dx, y1,y2,dy,                        &
       isize,jsize,title(1:length),runname,                             &
       tem4,slicopt,pltopt)

  CALL xpscmnt('End plotting '//label_copy(1:llabel))

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


SUBROUTINE overlay (layovr) 14

!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!    Set the layover counter parameter in the laypar common block
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Ming Xue
!
!  MODIFICATION HISTORY:
!    6/08/92  Added full documentation (K. Brewster)
!
!    8/08/93 (MX)
!    Automatically set the overlay parameter when input is not zero.
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!    layovr   The 'overlay' parameter.
!             If layover .ne. 0, the following 2-d contour plot will be
!             superimposed on the previous plot.
!             layover =1, 2, ... indicating this is the
!                layover'th (1st or 2nd ...) plot to be overlayed
!                on the previous one.
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  INTEGER :: layovr
!
!-----------------------------------------------------------------------
!
!  Plotting control common blocks
!
!-----------------------------------------------------------------------
!
  INTEGER :: layover, first_frame
  COMMON /laypar/ layover
  DATA first_frame /1/
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  IF( first_frame == 1 .OR. layovr == 0 ) THEN
    layover = 0
  ELSE
    layover = layover + 1
  END IF

  first_frame = 0

  RETURN
END SUBROUTINE overlay

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


SUBROUTINE styxrt( yxrt )
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!    Set the scaling factor of the y/x ratio of the plot.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Ming Xue
!
!  MODIFICATION HISTORY:
!    6/08/92  Added full documentation (K. Brewster)
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!    yxrt      Ratio of height to length of plot space
!              Default is set in the main program to 1.0
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
  REAL :: yxrt
!
!-----------------------------------------------------------------------
!
!  Plotting control common blocks
!
!-----------------------------------------------------------------------
!
  REAL :: yxratio
  COMMON /yratio/ yxratio       ! the scaling factor the y/x ratio.
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  yxratio = yxrt

  RETURN
END SUBROUTINE styxrt
!
!##################################################################
!##################################################################
!######                                                      ######
!######                 FUNCTION XFINC                       ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


  REAL FUNCTION xfinc(x)

!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!    Automatically divide domain (0,x) to a number of subdomain
!    with interval xfinc which is >=4 and =<16 for fold=1.0
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Ming Xue
!    sometime
!
!  MODIFICATIONS:
!    6/09/92  Added full documentation (K. Brewster)
!
!-----------------------------------------------------------------------
!
!  INPUT:
!     x       not sure
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
!
  REAL :: x
!
!-----------------------------------------------------------------------
!
!  Misc. local variables
!
!-----------------------------------------------------------------------
!
  INTEGER :: ipower
  REAL :: d,fold
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  ipower= INT( ALOG10(x) )
  d= INT(x/(10.0**ipower))
  fold=1.0
  xfinc=0.1*x
  IF( d >= 0.0 .AND. d < 3.0 ) THEN
    xfinc=2.0*10.0**(ipower-1)
  ELSE IF( d >= 3.0 .AND. d < 7.0 ) THEN
    xfinc=5.0*10.0**(ipower-1)*fold
  ELSE IF( d >= 7.0 .AND. d < 10. ) THEN
    xfinc=1.0*10.0** ipower*fold
  END IF
  IF(xfinc == 0.0) xfinc=x*0.1
  RETURN
  END FUNCTION xfinc
!
!##################################################################
!##################################################################
!######                                                      ######
!######                SUBROUTINE CLIPWD                     ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE clipwd(x1,y1,x2,y2,idispl) 1,2
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Return the portion of a line that is within a given window
!  (xw1,xw2,yw1,yw2)
!
!  If the given line is completely outside the window,
!  idispl=0, otherwise, idispl=1.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Ming Xue
!  3/6/93
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!    x1       value of x for first i grid point to plot
!    x2       value of x for last i grid point to plot
!    y1       value of y for first j grid point to plot
!    y2       value of y for last j grid point to plot
!
!    idispl   line orientation indicator
!       idispl = 0, the given line is completely outside the window
!       idispl = 1, the given line is partly inside the window
!
!-----------------------------------------------------------------------
!

!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
!
  REAL :: x1,x2,y1,y2
  INTEGER :: idispl
!
!-----------------------------------------------------------------------
!
!  Common blocks for plotting control parameters
!
!-----------------------------------------------------------------------
!
  REAL :: xw1,xw2,yw1,yw2
  COMMON /pltwdw/ xw1,xw2,yw1,yw2
  INTEGER :: ic1(4),ic2(4)
!
!-----------------------------------------------------------------------
!
!  Misc. local Variables
!
!-----------------------------------------------------------------------
!
  INTEGER :: i,knt,isw
  REAL :: x0,y0
  REAL :: isum1,isum2,ic01,ic02,ic03,ic04
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!

  knt = 0
  5     knt = knt+1
  CALL encodwd(x1,y1,ic1)
  CALL encodwd(x2,y2,ic2)

  isum1=ic1(1)+ic1(2)+ic1(3)+ic1(4)
  isum2=ic2(1)+ic2(2)+ic2(3)+ic2(4)
  idispl=1
  IF(isum1+isum2 == 0) GO TO 999

  idispl=0
  DO i=1,4
    20    IF(ic1(i)+ic2(i) == 2) GO TO 999
  END DO
!
!-----------------------------------------------------------------------
!
!  Make sure (x1,y1) is outside the window
!
!-----------------------------------------------------------------------
!
  isw=0
  15    IF(isum1 == 0) THEN
    ic01=ic1(1)
    ic02=ic1(2)
    ic03=ic1(3)
    ic04=ic1(4)
    DO i=1,4
      ic1(i)=ic2(i)
    END DO
    ic2(1)=ic01
    ic2(2)=ic02
    ic2(3)=ic03
    ic2(4)=ic04
    x0=x1
    y0=y1
    x1=x2
    y1=y2
    x2=x0
    y2=y0
    isw=1
  END IF

  IF(ic1(1) == 1) THEN
    y1=y1+(xw1-x1)*(y2-y1)/(x2-x1)
    x1=xw1
  ELSE IF(ic1(2) == 1) THEN
    y1=y1+(xw2-x1)*(y2-y1)/(x2-x1)
    x1=xw2
  ELSE IF(ic1(3) == 1) THEN
    x1=x1+(yw1-y1)*(x2-x1)/(y2-y1)
    y1=yw1
  ELSE IF(ic1(4) == 1) THEN
    x1=x1+(yw2-y1)*(x2-x1)/(y2-y1)
    y1=yw2
  END IF

  IF(isw == 1) THEN
    x0=x1
    y0=y1
    x1=x2
    y1=y2
    x2=x0
    y2=y0
  END IF

  idispl=1

  IF(knt > 10) THEN
    WRITE(6,*)'Dead loop encountered in CLIPWD, job stopped.'
    STOP 991
  END IF

  GO TO 5

  999   RETURN
END SUBROUTINE clipwd
!
!
!##################################################################
!##################################################################
!######                                                      ######
!######                SUBROUTINE ENCODWD                    ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE encodwd(x,y,ic) 2
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Encode a line section for window clipping purpose.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Ming Xue
!  3/6/93
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!    x       value of x
!    y       value of y
!    ic
!
!-----------------------------------------------------------------------
!

!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
!
  REAL :: x,y
  INTEGER :: ic(4)
!
!-----------------------------------------------------------------------
!
!  Common blocks for plotting control parameters
!
!-----------------------------------------------------------------------
!
  REAL :: xw1,xw2,yw1,yw2
  COMMON /pltwdw/ xw1,xw2,yw1,yw2
!
!-----------------------------------------------------------------------
!
!  Misc. local Variables
!
!-----------------------------------------------------------------------
!
  INTEGER :: i
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  DO i=1,4
    ic(i)=0
  END DO

  IF(x < xw1) ic(1)=1
  IF(x > xw2) ic(2)=1
  IF(y < yw1) ic(3)=1
  IF(y > yw2) ic(4)=1

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


SUBROUTINE ctrcol (icol,icol0) 9

!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!    Set the color  for field to plotted by CTR2D.
!
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Min Zou
!
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!    icol    begin color index
!    icol0   end color index
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
  INTEGER :: icol,icol0
!
!-----------------------------------------------------------------------
!
!  Plotting control common blocks
!
!-----------------------------------------------------------------------
!
  INTEGER :: icolor,icolor1,lbcolor,trcolor
  COMMON /recolor/icolor,icolor1,lbcolor,trcolor

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

  icolor=icol
  icolor1=icol0

  RETURN
END SUBROUTINE ctrcol


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


SUBROUTINE ctrvtr (units0,type0) 6

!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!    Set the units and type for plot wind
!
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Min Zou
!
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!    units0   the units of wind
!    type0    the type of wind
!    wcolor0  the color index
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
  INTEGER :: units0,type0
!
!-----------------------------------------------------------------------
!
!  Plotting control common blocks
!
!-----------------------------------------------------------------------
!
  INTEGER :: iunits, itype
  COMMON /windvtr/iunits, itype
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!

  iunits = units0
  itype = type0

  RETURN
END SUBROUTINE ctrvtr


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


SUBROUTINE varplt( var ) 13

!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!    Set the variable plot name for xconta.
!
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Min Zou
!
!  MODIFICATION HISTORY:
!    3/28/96
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!    var   variable plot name
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
  CHARACTER (LEN=*) :: var
  CHARACTER (LEN=12) :: varname

  COMMON /varplt1/ varname

  varname=var

  RETURN
END SUBROUTINE varplt

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


SUBROUTINE vtrsfc(u,v, x,y, xw,xe,dx, ys,yn,dy,                         & 3,7
           nx,ibgn,iend,ist, ny,jbgn,jend,jst,                          &
           label,time, runname, factor, slicopt,                        &
           tem1,tem2,tem3,tem4,                                         &
           tem5,tem6,hterain)

!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!    Plot vector fields in 2-d array
!
!  AUTHOR: Min Zou
!  4/28/97
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!    u        2-dimensional array of u wind components (m/s)
!    v        2-dimensional array of v wind components (m/s)
!
!    x        x coordinate of grid points in physical/comp. space (m)
!    y        y coordinate of grid points in physical/comp. space (m)
!    z        z coordinate of grid points in physical space (m)
!
!    xw       value of x for first i grid point to plot
!    xe       value of x for last i grid point to plot
!    ys       value of y for first j grid point to plot
!    yn       value of y for last j grid point to plot
!
!    nx       first dimension of b
!    ibgn     index of first i grid point to plot
!    iend     index of last  i grid point to plot
!
!    ny       second dimension of b
!    jbgn     index of first j grid point to plot
!    jend     index of last  j grid point to plot
!
!
!    time     time of data in seconds
!
!    runname  character string decribing run
!
!    factor   scaling factor for winds
!             V*factor wind vectors are plotted
!    slicopt  slice orientation indicator
!       slicopt = 1, x-y slice of u,v at z index kslice is plotted.
!       slicopt = 2, x-z slice of u,w at y index jslice is plotted.
!       slicopt = 3, y-z slice of v,w at x index islice is plotted.
!       slicopt = 4, x-y slice of u,v at z index islice is plotted.
!       slicopt = 5, xy-z cross section of wind islice is plotted.
!       slicopt = 6, data field on constant p-level  is plotted.
!       slicopt = 0, all of the three slices above are plotted.
!
!  WORK ARRAYS:
!
!    tem1     Temporary work array.
!    tem2     Temporary work array.
!    tem3     Temporary work array.
!    tem4     Temporary work array.
!    tem5     Temporary work array.
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
!
  INTEGER :: nx,ny
!
  REAL :: u(nx,ny)
  REAL :: v(nx,ny)

  REAL :: x(nx,ny)
  REAL :: y(nx,ny)
!

  CHARACTER (LEN=*) :: runname
  CHARACTER (LEN=*) :: label

  REAL :: xw,xe,dx,ys,yn,dy
  INTEGER :: ibgn,iend,ist, jbgn,jend,jst

  REAL :: time,factor
  INTEGER :: slicopt

  INTEGER :: iunits, itype
  COMMON /windvtr/iunits, itype

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

  REAL :: xw1,xe1,ys1,yn1
  COMMON /xuvpar/xw1,xe1,ys1,yn1
!
!-----------------------------------------------------------------------
!
!  Work arrays: tem1,tem2,tem3,tem4,tem5 of size at least
!          max( nx*ny, nx*nz, ny*nz).
!
!-----------------------------------------------------------------------
!
  REAL :: tem1(*),tem2(*),tem3(*),tem4(*),tem5(*)
  REAL :: tem6(*)
!
!-----------------------------------------------------------------------
!
!  Common blocks for plotting control parameters
!
!-----------------------------------------------------------------------
!
  REAL :: x01,y01                  ! the first  point of interpolation
  REAL :: x02,y02                  ! the second point of interpolation
  REAL :: zlevel                   ! the given height of the slice
  REAL :: sinaf,cosaf,dist,sqrtdxy
  COMMON /slicev/x01,y01,x02,y02,sinaf,cosaf,dist,sqrtdxy
  COMMON /sliceh/zlevel

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

!
!-----------------------------------------------------------------------
!
!  Misc. local Variables
!
!-----------------------------------------------------------------------
!
  INTEGER :: i,j,ij,istep,jstep,length,isize,jsize
  REAL :: uunit
  CHARACTER (LEN=6) :: timhms
  CHARACTER (LEN=120) :: title
!
  INTEGER :: icolor,icolor1,lbcolor,trcolor        ! required color
  COMMON /recolor/icolor,icolor1,lbcolor,trcolor
!
  INTEGER :: trnplt                ! flag to plot terain (1 or 0)
  REAL :: hterain(nx,ny)           ! The height of the terrain.

  INTEGER :: ovrtrn         ! overlay terrain option (0/1)

  REAL :: trninc,trnmin, trnmax    ! terrain interval minimum, maximum
  REAL :: ztmin,ztmax
  COMMON /trnpar/ trnplt,ovrtrn,trninc,trnmin, trnmax,ztmin,ztmax
!

  INTEGER :: smooth
  COMMON /smoothopt/smooth

  INTEGER :: wrtflag, llabel
  CHARACTER (LEN=80) :: levlab
  CHARACTER (LEN=50) :: timelab
  CHARACTER (LEN=25) :: timestring
  COMMON /timelev/wrtflag,timelab, levlab, timestring
  CHARACTER (LEN=120) :: label_copy
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!
  isize=(iend-ibgn)+1
  jsize=(jend-jbgn)+1
!
!
!-----------------------------------------------------------------------
!
!  Set up terrain, if needed.
!
!-----------------------------------------------------------------------
!
  label_copy = label
  llabel = 120
  CALL xstrlnth(label_copy, llabel)

  CALL xpscmnt('Start plotting '//label_copy(1:llabel))

  IF(trnplt == 1 .OR.trnplt == 2 .OR. ovrtrn == 1)  THEN
    DO j=jbgn,jend
      DO i=ibgn,iend
        ij = i-ibgn+1 + (j-jbgn)*isize
        tem5(ij)=hterain(i,j)
      END DO
    END DO
  END IF

  CALL cvttim ( time, timhms)

  IF( timhms(1:1) == '0' ) timhms(1:1)=' '

  WRITE(timelab,'(''T='',F8.1,A)') time,                                &
      ' s ('//timhms(1:2)//':'//timhms(3:4)//':'//timhms(5:6)//')'
  CALL get_time_string ( time, timestring)

!   length=50
!   CALL strmin(timelab,length)
!   write(timelab,'(a,'' '',a)') timestring(1:21), timelab(1:length)
!   print*,'in vtrsfc', timelab


  DO j=jbgn,jend
    DO i=ibgn,iend
      ij = i-ibgn+1 + (j-jbgn)*isize
      tem1(ij) = -9999.0
      tem2(ij) = -9999.0
      IF(      IF(      tem3(ij)=x(i,j)
      tem4(ij)=y(i,j)
    END DO
  END DO

  levlab = 'First level above ground (surface)'
  WRITE(title,'(''U-V '',A)') label

  length = 120
  CALL strlnth( title, length )
  CALL strmin ( title, length)

  uunit = 10.0
  CALL xvmode(1)
  istep = ist
  jstep = jst

  DO i=1,smooth
    CALL smooth9pmv(tem1,isize,jsize,1,isize,1,jsize,tem6)
    CALL smooth9pmv(tem2,isize,jsize,1,isize,1,jsize,tem6)
  END DO

  CALL vtr2d(tem1,tem2,tem3,tem4, uunit, xw,xe,dx,ys,yn,dy,             &
        isize,istep,jsize,jstep,                                        &
        title(1:length),runname, 1,                                     &
        tem5,slicopt)

  CALL xpscmnt('End plotting '//label_copy(1:llabel))

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


SUBROUTINE set_interval(z, m,n,zmin1,zmax1,ctmin,ctmax,cl) 1
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!    Limited contour interval when uinc = -9999.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Min Zou
!
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!    z        2-D array
!    m,n      dimension of 2-d array
!    zmin1    the minimum value of 2-D array
!    zmax1    The maximum value of 2-D array
!    ctmin    the input minimum value
!    ctmax    the input maximum value
!    cl       the intervals
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!

  REAL :: z(m,n),cl(*)
  REAL :: zmin, zmax
  REAL :: zinc
  COMMON /xclm19/ nmin, nmax
  COMMON /xcrf17/clref,lcptn,labtyp,iclf,lhilit,ihlf,kct0
  COMMON /zchole/ nhole,specia,nvtrbadv
  COMMON /xoutch/ nch

  IF(ctmin == 0.0 .AND. ctmax == 0.0) THEN
    cl(2)=cl(1)+ xfinc(zmax1-zmin1)/2
  ELSE
    cl(2)=cl(1)+ xfinc(ctmax-ctmin)/2
  END IF
  IF(cl(2)-cl(1) == 0.0) cl(2)=cl(1)+1.0
  zinc = cl(2)-cl(1)

  ncmin=nmin
  ncmax=nmax
  diff=ctmax-ctmin
  IF(diff-ABS( zinc)*1.0E-6  > 0.0) THEN
    GO TO     4
  END IF
  999  WRITE(nch,'(a,a)')                                               &
       ' Bad first guess of contour increment or field is constant',    &
       ', number of contours is one.'
  ncnt=1
  cl(1)= ctmin
  899  RETURN
  4    kcount=0
  1    CONTINUE
  eps=0.1*zinc
  kcount=kcount+1
  IF( kcount > 20) GO TO 998
  kzinc=(ctmin-clref)/zinc
  zmin=kzinc*zinc+clref
  kzinc=(ctmax-clref)/zinc
  zmax=kzinc*zinc+clref
  IF(ctmin-clref > 0.0) zmin=zmin+zinc
  IF(ctmax-clref < 0.0) zmax=zmax-zinc
!
  clv=zmin-zinc
  ncnt=0
  6    clv=clv+zinc
  IF(clv-zmax-eps > 0.0) THEN
    GO TO     8
  END IF
  3    ncnt=ncnt+1
  IF(ncnt > ncmax) THEN
    zinc=zinc*2
    WRITE(nch,1000) ncmax, zinc
    1000 FORMAT(' Number of contours > ',i3,' ,Zinc is doubled. Zinc='  &
            ,e10.3)
    GO TO 1
  END IF
  IF( ABS( clv-clref ) < eps ) clv=clref
  cl(ncnt)=clv
  GO TO 6
  8    CONTINUE

  IF( ncnt < ncmin) THEN
    zinc=zinc/2
    WRITE(nch,2000) ncmin,zinc
    2000 FORMAT(' Number of contours < ',i3,' ,Zinc is halved. Zinc='   &
           ,e10.3)
    GO TO 1
  END IF
  WRITE(nch,'('' * NUMBER OF CONTOURS= '',I5,''  MIN='',E12.4,          &
  &   '' MAX='', e12.4,'' inc='',e12.5 )')                              &
      ncnt,ctmin,ctmax,zinc

  IF( zmin1 >= ctmin .AND. zmax1 <= ctmax) THEN
    zinc = cl(2) - cl(1)
    WRITE(nch,'(''SET MINIMUM CONTOUR INTERVAL IS'',E12.4,              &
    &   '' ctmin='',e12.4,'' ctmax='',e12.4 )')zinc,ctmin,ctmax
    CALL xctref(zinc)
    CALL xnctrs( 1,300)
  ELSE
    WRITE(nch,'(''NO NEED SET MINIMUM CONTOUR INTERVAL'')' )
    WRITE(nch,'(''CNTOUR INTERVAL IS SET AUTOMATICALLY'')' )
    cl(2)=cl(1)+ xfinc(zmax1-zmin1)/2
    IF(cl(2)-cl(1) == 0.0) cl(2)=cl(1)+1.0
  END IF
  RETURN
  998  WRITE(nch,*)' Contour levels can not be selected by XCNTLV.'
  WRITE(nch,*)                                                          &
      ' Plz alter input contour interval or limits of contour number'

END SUBROUTINE set_interval

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


SUBROUTINE xrch1( r,ch,lch) 1,3

! Return real number R as a character string in automatically set format
  REAL :: r
  CHARACTER (LEN=20) :: str
  CHARACTER (LEN=20) :: ch

  CALL get_format(r,str)
  IF(ABS(r-0.0) < 1.e-20) THEN
    WRITE(ch,'(F3.1)') r
  ELSE
    WRITE(ch,str) r
  END IF
  lch=20
  CALL strlnth( ch, lch)
  CALL strmin ( ch, lch)
  RETURN
END SUBROUTINE xrch1

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


SUBROUTINE get_format(r,ch) 1
  INTEGER :: npoz
  CHARACTER (LEN=1) :: FORM,ndrob
  CHARACTER (LEN=20) :: ch
  WRITE(ch,10)r
  10   FORMAT(g11.4)
  DO i=20,1,-1
    IF(ch(i:i) == '0'.OR.ch(i:i) == ' ') THEN
      ch(i:i)=' '
    ELSE
      GO TO 1
    END IF
  END DO
  1    CONTINUE
  npoz=0
  ndot=0
  nmant=0
  ndrob=' '
  FORM='F'
  DO i = 1,20
    IF(ch(i:i) /= ' ' ) npoz=npoz+1
    IF(ch(i:i) == 'E') FORM='E'
    IF(ndrob == '.'.AND.ch(i:i) /= ' ') ndot=ndot+1
    IF(ch(i:i) == '.') ndrob='.'
    IF(FORM /= 'E') nmant=npoz
  END DO
  npoz=npoz
  IF(FORM == 'F') THEN
    IF(ndot /= 0) THEN
      WRITE(ch,20) '(',FORM,npoz,'.',ndot,')'
    ELSE
      WRITE(ch,20) '(',FORM,npoz,'.',ndot,')'
    END IF
  ELSE IF(FORM == 'E') THEN
    ch = '(1PE20.2)'
  ELSE
    WRITE(ch,20) '(',FORM,npoz,'.',nmant,')'
  END IF
  20   FORMAT(a1,a1,i1,a1,i1,a1)
  RETURN
END SUBROUTINE get_format

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


SUBROUTINE drawmap(nunit) 1
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  This subroutine will plot the map
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Ming Xue
!
!  MODIFICATION HISTORY:
!    6/2/97 Min Zou
!    Read multiple mapfile only once. Using differnt line style to
!    plot mapdata.
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!    nunit     the channel of the mapfile data
!    mapfile   character of map file name
!
!-----------------------------------------------------------------------
!
!  Variable Declarations
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE

  INCLUDE 'arpsplt.inc'

  INTEGER :: nunit,i
  INTEGER :: lmapfile

  CHARACTER (LEN=132) :: mapfile(maxmap)
  INTEGER :: mapgrid,mapgridcol, kolor
  REAL :: latgrid,longrid
  INTEGER :: nmapfile,mapcol(maxmap),mapline_style(maxmap)
  COMMON /mappar1/nmapfile,mapcol,mapline_style,mapfile
  COMMON /mappar2/mapgrid,mapgridcol,latgrid,longrid

  REAL :: x1,x2,y1,y2

  CALL xpscmnt('Start of map plotting ')

  CALL xqmap (x1,x2,y1,y2)
  CALL xwindw(x1,x2,y1,y2)
  CALL xqcolor(kolor)

  DO i=1,nmapfile
    CALL xcolor(mapcol(i))
    IF(mapline_style(i) == 1) THEN
      CALL xthick(1)
      CALL xbrokn(6,3,6,3)
    ELSE IF(mapline_style(i) == 2) THEN
      CALL xthick(1)
    ELSE IF(mapline_style(i) == 3) THEN
      CALL xthick(3)
      CALL xfull
    END IF

    lmapfile=132
    CALL xstrlnth(mapfile(i), lmapfile)

!    print*,'plotting ',mapfile(i)(1:lmapfile)

    CALL xdrawmap(nunit,mapfile(i)(1:lmapfile),latgrid,longrid)

  END DO

  CALL xcolor(kolor)
  CALL xfull
  CALL xwdwof

  CALL xpscmnt('End of map plotting ')

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


SUBROUTINE pltobs(obopt) 2
!
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!  Plots observations on an arpsplt contour map.
!
!-----------------------------------------------------------------------
!
!  INPUT:
!  obopt    Plotting option
!           1  Plot data in obs1 as characters
!           2  Plot data in obs1 and obs2 as characters
!           3  Plot wind arrows with obs1 as u and obs2 as v.
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!

  IMPLICIT NONE

  INCLUDE 'arpsplt.inc'
!
!  Arguments
!
  INTEGER :: obopt
!
  INTEGER :: nobs
  REAL :: latob(mxsfcob)
  REAL :: lonob(mxsfcob)
  REAL :: obs1(mxsfcob)
  REAL :: obs2(mxsfcob)
!
  COMMON /sfc_obs1/ nobs
  COMMON /sfc_obs2/ latob,lonob,obs1,obs2
!
!  Plotting parameters
!
  CHARACTER (LEN=1) :: cross
  PARAMETER(cross='+')
!
  INTEGER :: ovrobs,obsset,obscol,obs_marktyp
  REAL :: obs_marksz
  COMMON /obspar/ ovrobs,obsset,obscol,obs_marktyp, obs_marksz
!
!  Plotting common blocks
!
  INTEGER :: icolor,icolor1,lbcolor,trcolor       ! required color
  COMMON /recolor/icolor,icolor1,lbcolor,trcolor

  REAL :: ctinc,ctmin,ctmax,vtunt  ! contour interval and vector unit
  COMMON /incunt/ ctinc,ctmin,ctmax,vtunt

  REAL :: xleng,vunit
  COMMON /vecscl/ xleng,vunit

  INTEGER :: iunits, itype
  COMMON /windvtr/iunits, itype
!
!  Misc local variables
!
  INTEGER :: iob
  REAL :: orgmag,obmag,yoff,yoff2
  REAL :: x1,x2,y1,y2
  REAL :: xob,yob
  CHARACTER (LEN=4) :: chplot
  INTEGER :: imkrfil
!
!  Set-up plotting space and zxplot variables
!
  CALL xcolor(lbcolor)
  CALL xqmap (x1,x2,y1,y2)
  CALL xwindw(x1,x2,y1,y2)
  CALL xchori(0.0)
  CALL xqchsz(orgmag)

  obmag=0.8*orgmag
  yoff=0.5*orgmag
  yoff2=2.*yoff
  CALL xchsiz(obmag)

  IF(obopt == 1) THEN
    CALL xcolor(obscol)
    CALL xmrksz(obs_marksz)
    DO iob=1,nobs
      IF(obs1(iob) > -98. .AND. obs1(iob) < 500.) THEN
        CALL xlltoxy(1,1,latob(iob),lonob(iob),xob,yob)
        WRITE(chplot,810) nint(obs1(iob))
        810       FORMAT(i4)
        CALL xcharc((0.001*xob),(0.001*yob+yoff),chplot)
!        call XCHARC((0.001*xob),(0.001*yob),cross)
        CALL xmarker((0.001*xob),(0.001*yob),obs_marktyp)
        IF(obs_marktyp > 5) THEN
          CALL xqmkrfil(imkrfil)
          CALL xmkrfil(1)
          CALL xmarker((0.001*xob),(0.001*yob),MOD(obs_marktyp,5))
          CALL xmkrfil(imkrfil)
        END IF
      END IF
    END DO
  ELSE IF(obopt == 2) THEN
    CALL xcolor(obscol)
    CALL xmrksz(obs_marksz)
    DO iob=1,nobs
      CALL xlltoxy(1,1,latob(iob),lonob(iob),xob,yob)
      IF(obs1(iob) > -98. .AND. obs1(iob) < 500.) THEN
        WRITE(chplot,810) nint(obs1(iob))
        CALL xcharc((0.001*xob),(0.001*yob+yoff),chplot)
!        call XCHARC((0.001*xob),(0.001*yob),cross)
        IF(obs_marktyp > 5) THEN
          CALL xqmkrfil(imkrfil)
          CALL xmkrfil(1)
          CALL xmarker((0.001*xob),(0.001*yob),MOD(obs_marktyp,5))
          CALL xmkrfil(imkrfil)
        END IF
      END IF
      IF(obs2(iob) > -98. .AND. obs2(iob) < 500.) THEN
        WRITE(chplot,810) nint(obs2(iob))
        CALL xcharc((0.001*xob),(0.001*yob+yoff),chplot)
!         call XCHARC((0.001*xob),(0.001*yob),cross)
        CALL xmarker((0.001*xob),(0.001*yob),obs_marktyp)
        IF(obs_marktyp > 5) THEN
          CALL xqmkrfil(imkrfil)
          CALL xmkrfil(1)
          CALL xmarker((0.001*xob),(0.001*yob),MOD(obs_marktyp,5))
          CALL xmkrfil(imkrfil)
        END IF
      END IF
    END DO
  ELSE IF(obopt == 3) THEN
    CALL xcolor(obscol)
    CALL xmrksz(obs_marksz)
    DO iob=1,nobs
      IF(obs1(iob) > -98. .AND. obs1(iob) < 500. .AND.                  &
            obs2(iob) > -98. .AND. obs2(iob) < 500.) THEN
        CALL xlltoxy(1,1,latob(iob),lonob(iob),xob,yob)
        xob=0.001*xob
        yob=0.001*yob
        IF( xob > x1 .AND. xob < x2 .AND. yob > y1 .AND. yob < y2 ) THEN
          CALL xarrow(obs1(iob),obs2(iob),xob,yob,xleng,vunit)
          CALL xmarker(xob,yob,obs_marktyp)
          IF(obs_marktyp > 5) THEN
            CALL xqmkrfil(imkrfil)
            CALL xmkrfil(1)
            CALL xmarker(xob,yob,MOD(obs_marktyp,5))
            CALL xmkrfil(imkrfil)
          END IF
        END IF
      END IF
    END DO
  END IF

  CALL xcolor(lbcolor)
  CALL xchsiz(orgmag)
  CALL xfull
  CALL xwdwof

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


SUBROUTINE pltsta(a,b,x,y,m,n,flag,slicopt) 6,5
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  This subroutine will plot some station information.
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
!  AUTHOR:
!    Min Zou (6/1/97)
!
!  Modification history:
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!    a, b    2-dimension array of Variable
!    m, n    Array dimensions
!    x, y    x-coord and y-coord of the staions
!    flag    a flag for different plot
!    slicopt  slice orientation indicator
!
!-----------------------------------------------------------------------
!
!  Variable Declarations
!
!-----------------------------------------------------------------------
!

  IMPLICIT NONE
  INCLUDE 'arpsplt.inc'

  INTEGER :: m,n
  REAL :: a(m,n)
  REAL :: b(m,n)
  REAL :: x(m,n)
  REAL :: y(m,n)
  INTEGER :: nsta,nstapro(mxstalo),nstatyp(mxstalo)
  REAL :: latsta(mxstalo), lonsta(mxstalo)
  CHARACTER (LEN=5) :: s_name(mxstalo)
  INTEGER :: ovrstaopt
  INTEGER :: ovrstam,staset,ovrstan,ovrstav,stacol,markprio,wrtstax
  INTEGER :: nsta_typ,sta_typ(10),sta_marktyp(10), sta_markcol(10)
  REAL :: sta_marksz(10)
  REAL :: wrtstad
  COMMON /sta_par/ ovrstaopt,ovrstam,staset,ovrstan,ovrstav,stacol,     &
         markprio, nsta_typ,sta_typ,sta_marktyp,                        &
         sta_markcol,sta_marksz,wrtstax,wrtstad
  COMMON /sta_loc/latsta,lonsta,nstatyp,nstapro,nsta
  COMMON /sta_loc1/s_name
  REAL :: xob(mxstalo), yob(mxstalo),aob(mxstalo),bob(mxstalo)
  COMMON /xob_yob/xob, yob
  INTEGER :: LEN,i,j
!
  REAL :: x01,x02,y01,y02
  REAL :: sinaf,cosaf,dist,sqrtdxy
  COMMON /slicev/x01,y01,x02,y02,sinaf,cosaf,dist,sqrtdxy
!
  INTEGER :: icolor,icolor1,lbcolor,trcolor       ! required color
  COMMON /recolor/icolor,icolor1,lbcolor,trcolor
!
  INTEGER :: layover
  COMMON /laypar/ layover
!
  CHARACTER (LEN=12) :: varname
  COMMON /varplt1/ varname
!
  REAL :: xori1,xori2,yori1,yori2,zori1,zori2
  COMMON /tmphc1/xori1,xori2,yori1,yori2,zori1,zori2
!
  REAL :: xleng,vunit
  COMMON /vecscl/ xleng,vunit

  INTEGER :: iunits, itype
  COMMON /windvtr/iunits, itype

  REAL :: x_tmp
  COMMON /tmphc2/ x_tmp
!
  REAL :: x1,x2,y1,y2
  REAL :: orgmag,obmag,yoff,yoff2,xoff, xoff2
  CHARACTER (LEN=30) :: ctmp
!
  INTEGER :: flag,slicopt,fg
  REAL :: xdist, ydist,xd0,yd0,xa,xb
  SAVE fg

  REAL :: xleng0, spd, dir, istand
  INTEGER :: iunits0, imkrfil
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! calculate xob and yob on the xy plane
  CALL xwindw(xori1,xori2,yori1,yori2)
  IF(fg == 0) THEN
    CALL xlltoxy(nsta,1,latsta,lonsta,xob,yob)
    DO i= 1,nsta
      xob(i) = xob(i)*0.001
      yob(i) = yob(i)*0.001
    END DO
    fg=1
  END IF
!
  CALL xqmap (x1,x2,y1,y2)
  CALL xwindw(x1,x2,y1,y2)
  CALL xchori(0.0)
  CALL xqchsz(orgmag)

  obmag=0.8*orgmag
  yoff=0.5*orgmag
  yoff2=3.*yoff
  xoff = 0.001*(x2-x1)
  xoff2 = 4.*xoff
  CALL xchsiz(obmag)


  IF(ovrstav == 1 ) THEN  ! interpolation
    CALL intepo (nsta,xob,yob,aob,m,n,x,y,a)
    IF(varname(1:6) == 'vtrplt' .OR. varname(1:6) == 'vtpplt')          &
        CALL intepo (nsta,xob,yob,bob,m,n,x,y,b)
  END IF

  IF (flag == 1) THEN
    CALL xcolor(stacol)
    IF(slicopt == 5) THEN
      xa = (y02-y01)/(x02-x01)
      xb = y01 - xa*x01
    END IF
    DO i = 1,nsta
      IF(nstapro(i) <= markprio) THEN
        IF(wrtstax == 1 )THEN
          CALL xwindw(x1,x2-0.005*(x2-x1),                              &
                    y1-0.3*(y2-y1),y1)
          LEN=5
          CALL strlnth(s_name(i),LEN)
          CALL xchori(90.)
          IF(slicopt == 2 .OR. slicopt == 10) THEN
            CALL xwindw(xori1,xori2-0.005*(xori2-xori1),                &
                    y1-0.3*(y2-y1),y1)
            IF ( (xob(i) <= xori2.AND.xob(i) >= xori1)                  &
                  .AND. (yob(i) <= yori2.AND.yob(i) >= yori1) ) THEN
              IF( ABS(yob(i)-x_tmp) <= wrtstad ) THEN
!               CALL XCHARR((xob(i)),y1-1.75*yoff2,
                CALL xcharr((xob(i)),y1-1.50*yoff2,                     &
                    s_name(i)(1:LEN))
              END IF
            END IF
          ELSE IF(slicopt == 3 .OR. slicopt == 11) THEN
            CALL xwindw(yori1,yori2-0.005*(yori2-yori1),                &
                    y1-0.3*(y2-y1),y1)
            IF ( (xob(i) <= xori2.AND.xob(i) >= xori1)                  &
                  .AND. (yob(i) <= yori2.AND.yob(i) >= yori1) ) THEN
              IF( ABS(xob(i)-x_tmp) <= wrtstad ) THEN
                CALL xcharr((yob(i)),y1-1.75*yoff2,                     &
                    s_name(i)(1:LEN))
              END IF
            END IF
          ELSE IF( slicopt == 5) THEN
            CALL xwindw(x1,x2-0.005*(x2-x1),                            &
                    y1-0.3*(y2-y1),y1)
            IF ( (xob(i) <= xori2.AND.xob(i) >= xori1)                  &
                  .AND. (yob(i) <= yori2.AND.yob(i) >= yori1) ) THEN
              xd0 = 1./(xa*xa+1.0)*((yob(i)-xb)*xa+xob(i))
              yd0 = xa*xd0+xb
              xdist = SQRT((x01-xd0)*(x01-xd0) + (y01-yd0)*(y01-yd0))
              ydist = SQRT((xob(i)-xd0)*(xob(i)-xd0)+                   &
                         (yob(i)-yd0)*(yob(i)-yd0))
              xdist = xdist+x1
              IF(ABS(ydist) <= wrtstad ) THEN
                CALL xcharr(xdist,y1-1.75*yoff2,                        &
                    s_name(i)(1:LEN))
              END IF
            END IF
          END IF
          CALL xchori(0.)
        END IF
      END IF
    END DO

  ELSE IF(flag == 0) THEN
    CALL xwindw(x1,x2,y1,y2)
    DO i = 1,nsta
      IF( (xob(i) >= x1.AND.xob(i) <= x2) .AND.                         &
            (yob(i) >= y1.AND.yob(i) <= y2) ) THEN
        IF(nstapro(i) <= markprio) THEN
          IF(ovrstan == 1) THEN
            LEN=5
            CALL strlnth(s_name(i),LEN)
            CALL xcharc((xob(i)),(yob(i)-yoff2),                        &
                  s_name(i)(1:LEN))
          END IF
          IF(ovrstam == 1) THEN
            DO j=1,nsta_typ
              CALL xmrksz(sta_marksz(j))
              CALL xcolor(sta_markcol(j))
              IF(nstatyp(i) == sta_typ(j)) THEN
                CALL xmarker((xob(i)),(yob(i)),                         &
                     sta_marktyp(j))

                IF(sta_marktyp(j) > 5) THEN
                  CALL xqmkrfil(imkrfil)
                  CALL xmkrfil(1)
                  CALL xmarker((xob(i)),(yob(i)),                       &
                               MOD(sta_marktyp(j) ,5))
                  CALL xmkrfil(imkrfil)
                END IF

                IF(ovrstan == 1) THEN
                  LEN=5
                  CALL strlnth(s_name(i),LEN)
                  CALL xcharc((xob(i)),(yob(i)-yoff2),                  &
                      s_name(i)(1:LEN))
                END IF
              END IF
            END DO
          END IF
          IF(ovrstav == 1) THEN
            CALL xcolor(stacol)
            IF(varname(1:6) == 'vtrplt' .OR. varname(1:6) == 'vtpplt') THEN
!          IF(i.eq.1) THEN
              xleng0=xleng*0.0004
              IF(iunits == 2 ) THEN
                iunits0=2
                istand = 10.
                WRITE(ctmp,'(a30)')'10 knots'
              ELSE IF (iunits == 3) THEN
                iunits0=2
                istand = 10.
                WRITE(ctmp,'(a30)')'10 MPH'
              ELSE IF(iunits == 1) THEN
                iunits0=1
                istand = 5.
                WRITE(ctmp,'(a30)')'5 m/s'
              END IF
!          ENDIF
              IF(aob(i) /= -9999. .AND. bob(i) /= -9999.) THEN
                spd = SQRT(aob(i)*aob(i)+bob(i)*bob(i))
                dir = ATAN2(-1.*aob(i),-1.*bob(i))*180./3.1415926
                IF(dir <= 0.) dir = 360.+dir
!            CALL barb((xob(i)),
!    :                   (yob(i)),dir,spd,iunits0-1, xleng0)
                CALL xbarb(aob(i),bob(i),xob(i),yob(i),                 &
                     iunits0,xleng*0.65,2)

              END IF
            ELSE
              IF(aob(i) /= -9999.) THEN
                CALL xrch(aob(i),ctmp,LEN)
                IF(layover == 0) CALL xcharr((xob(i)-xoff2),            &
                    (yob(i)+yoff),ctmp(1:LEN))
                IF(layover == 1) CALL xcharl((xob(i)+xoff2),            &
                    (yob(i)+yoff) ,ctmp(1:LEN))
                IF(layover == 2) CALL xcharc((xob(i)+xoff2),            &
                    (yob(i)-yoff),ctmp(1:LEN))
                IF(layover == 3) CALL xcharl((xob(i))+xoff2,            &
                    (yob(i)-yoff) ,ctmp(1:LEN))
              END IF
            END IF
          END IF

        END IF
      END IF
    END DO
  END IF

  CALL xcolor(lbcolor)
  CALL xchsiz(orgmag)
  CALL xfull
  CALL xwdwof

  RETURN
END SUBROUTINE pltsta
!
!
!##################################################################
!##################################################################
!######                                                      ######
!######                SUBROUTINE RUNLAB                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE runlab(runname) 38
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!    Plot a run label at the lower left cornor of the picture frame.
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!    runname      character string of run label
!
!-----------------------------------------------------------------------
!
!  Variables Declarations
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
  CHARACTER (LEN=*) :: runname
  REAL :: xl, xr, yb, yt, rotang, xlimit, ylimit
  INTEGER :: nopic, nxpic, nypic
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  CALL xqmap(xl,xr,yb,yt)
  CALL xqnpic(nopic)
  CALL xqspac(nxpic, nypic, rotang, xlimit, ylimit)

  IF( rotang == 0.0 ) THEN

    IF(nopic == nxpic*nypic -(nxpic-1)) THEN
      CALL xcharl( xl, yb-0.15*(yt-yb), runname )
    END IF

  ELSE

    IF(nopic == nypic*nxpic -(nypic-1)) THEN
      CALL xcharl( xl, yb-0.15*(yt-yb), runname )
    END IF

  END IF

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


SUBROUTINE vprofil(nx,ny,nz,nzprofc,var,xc,yc,zpc,plwr,pupr,            & 27,2
           xpnt,ypnt,npoints,zlwr,zupr,xcaptn,ycaptn,npicprof,          &
           profil,height)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  This subroutine will plot the vertical profiles of a given
!  variable through points (xpnt(i),ypnt(i),i=1,npoints).
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
!  AUTHOR:
!    Adwait Sathye (2/28/94)
!
!  Modification history:
!
!    4/18/94, (Ming Xue)
!    Major overhaul. Many temporary arrays removed. New frame option
!    added.
!
!    9/18/1995 (Ming Xue)
!    Fixed a problem in the code that determines kbgn and kend.
!
!    10/8/1996 (Y. Richardson)
!    Corrected a bug in the interpolation weights.
!
!-----------------------------------------------------------------------
!

!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!    nx, ny, nz    Array dimensions
!    nzprofc       the maximum vertical index in height (zpc/zpsoilc)
!                  and variables to be profiled when calling vprofil
!                  subroutine. 06/10/2002, Zuwen He
!                  In the atmosphere model, the vertical index is
!                  typically nz-1, while in the soil model, it's nzsoil.
!    var           Variable data array
!    xc,yc,zpc     The coordinate of input data var.
!    plwr,pupr     Lower and upper bounds for the horiz. axis of profile
!    xpnt, ypnt    Arrays containing the X and Y locations of the
!                  mulitple profiles to be plotted
!    npoints       Number of profile points to be plotted
!    zlwr, zupr    Bounds in the vertical direction
!    xcaptn        Caption for the X axis
!    ycaptn        Caption for the Y axis
!
!  Work arrays:
!
!    profil,height Temporary arrays
!
!-----------------------------------------------------------------------
!

!
!-----------------------------------------------------------------------
!
!  Variable Declarations
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
!
!-----------------------------------------------------------------------
!
!  Variables passed in
!
!-----------------------------------------------------------------------
!
  INTEGER :: nx, ny, nz, nzprofc

  REAL :: var(nx,ny,nz)
  REAL :: xc(nx,ny,nz), yc(nx,ny,nz), zpc(nx,ny,nz)

  REAL :: plwr,pupr
  REAL :: zlwr, zupr


  INTEGER :: npoints
  REAL :: xpnt(npoints), ypnt(npoints)

  CHARACTER (LEN=*) :: xcaptn
  CHARACTER (LEN=*) :: ycaptn

  INTEGER :: npicprof

  REAL :: profil(nz,npoints)
  REAL :: height(nz,npoints)

  LOGICAL :: multiprof
!
!-----------------------------------------------------------------------
!
!  Temporary local variables
!
!-----------------------------------------------------------------------
!
  REAL :: lower, upper, zmin, zmax
  REAL :: x1, x2, y1, y2
  REAL :: a1, a2, a3, a4
  INTEGER :: i, j, k, ix, jy,kbgn,kend,ip,lchar
  REAL :: dx,dy,temp,hmaxk,hmink
  CHARACTER (LEN=80) :: ch

  REAL :: lblmag, ctrlbsiz, axlbsiz
  COMMON /labmag/ lblmag, ctrlbsiz, axlbsiz
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  multiprof = .false.
  IF( npicprof == 0 .AND. npoints > 1 ) multiprof = .true.
!
!-----------------------------------------------------------------------
!
!  If lower boundary is bigger, then swap boundaries.
!
!-----------------------------------------------------------------------
!
  IF (plwr > pupr) THEN
    lower = pupr
    upper = plwr
  ELSE
    lower = plwr
    upper = pupr
  END IF
!
!-----------------------------------------------------------------------
!
!  Find corresponding coordinates for the boundaries in the Z
!  dimension. IF both have been set to 0, use the boundary values.
!  Else, loop through the zpc array and find the location of the
!  point.
!
!-----------------------------------------------------------------------
!
  dx = xc(2,1,1)-xc(1,1,1)
  dy = yc(1,2,1)-yc(1,1,1)

  zmin = zpc(1,1,1)
  zmax = zpc(1,1,nzprofc)
  DO j=1,ny-1
    DO i=1,nx-1
      zmin=MIN( zmin, zpc(i,j,1))
      zmax=MAX( zmax, zpc(i,j,nzprofc))
      zmin=MIN( zmin, zpc(i,j,nzprofc))  ! because of soil model, Zuwen
      zmax=MAX( zmax, zpc(i,j,1))    ! because of soil model, Zuwen 
    END DO
  END DO

  IF( zlwr /= zupr ) THEN
    zmin = MAX(zmin, zlwr)
    zmax = MIN(zmax, zupr)
  END IF

  IF( zmax < zmin) WRITE(6,'(a, f10.3, a, f10.3)')                      &
      'Warning: zmax is less then zmin. Check input data zprofbgn',     &
      zlwr , 'and zprofbgn',zupr

  DO ip = 1, npoints

    ix = INT ( (xpnt(ip) - xc(1,1,1))/dx ) + 1
    jy = INT ( (ypnt(ip) - yc(1,1,1))/dy ) + 1
    ix = MIN(MAX(1,ix), nx-2)
    jy = MIN(MAX(1,jy), ny-2)

    WRITE(6,'(1x,2a,2(a,f10.3),2(a,i4) )')                              &
        'Plotting ',xcaptn,' profile through (',                        &
        xpnt(ip),',',ypnt(ip),') km, at i=',ix,' j=',jy
!
!-----------------------------------------------------------------------
!
!  Interpolate the data value and its height to the specified point.
!
!-----------------------------------------------------------------------
!
    a1 = ABS ((xc(ix  ,jy  ,1)-xpnt(ip))*(yc(ix  ,jy  ,1)-ypnt(ip)))
    a2 = ABS ((xc(ix+1,jy  ,1)-xpnt(ip))*(yc(ix+1,jy  ,1)-ypnt(ip)))
    a3 = ABS ((xc(ix+1,jy+1,1)-xpnt(ip))*(yc(ix+1,jy+1,1)-ypnt(ip)))
    a4 = ABS ((xc(ix  ,jy+1,1)-xpnt(ip))*(yc(ix  ,jy+1,1)-ypnt(ip)))

    DO k = 1,nzprofc

      profil(k,ip)= (a3*var(ix  ,jy,k)+a2*var(ix  ,jy+1,k)+             &
                     a4*var(ix+1,jy,k)+a1*var(ix+1,jy+1,k))             &
                    /(a1 + a2 + a3 + a4)

      height(k,ip)= (a3*zpc(ix  ,jy,k)+a2*zpc(ix  ,jy+1,k)+             &
                     a4*zpc(ix+1,jy,k)+a1*zpc(ix+1,jy+1,k))             &
                    /(a1 + a2 + a3 + a4)

    END DO

  END DO

  kbgn = nzprofc
  DO k=nzprofc,1,-1

    hmaxk = height(k,1)
    DO ip=1,npoints
      hmaxk = MAX(hmaxk,height(k,ip))
    END DO

    IF( hmaxk >= zmin) kbgn = k

  END DO

  kend = 1
  DO k=1,nzprofc

    hmink = height(k,1)
    DO ip=1,npoints
      hmink = MIN(hmink,height(k,ip))
    END DO

    IF( hmink <= zmax) kend=k

  END DO

!
!-----------------------------------------------------------------------
!
!  If input bounds for the profile are zero, use the min. and max.
!  in the profile as the lower and upper bounds for the horizontal
!  axis.
!
!-----------------------------------------------------------------------
!
  IF( plwr == 0.0 .AND. pupr == 0.0 ) THEN

    lower = profil(kbgn,1)
    upper = profil(kend,1)
    DO ip=1,npoints
      DO k = kbgn,kend
        lower = MIN(lower, profil(k,ip))
        upper = MAX(upper, profil(k,ip))
      END DO
    END DO

  ELSE

    lower = plwr
    upper = pupr

  END IF

!
!-----------------------------------------------------------------------
!
!    If the lower and upper bounds are equal, set the horizontal
!    axis scale to 1.0.
!
!-----------------------------------------------------------------------
!

  IF ((lower == 0.0 .AND. upper == 0.0).OR.upper == lower) upper = lower+1.0

!
!-----------------------------------------------------------------------
!
!  Start to plot the profile...
!
!-----------------------------------------------------------------------
!
  DO ip=1,npoints

    IF( (.NOT.multiprof) .OR. (multiprof.AND.ip == 1) ) THEN

      CALL xnwpic
      CALL xaxtik(1, 1)
      CALL xaxant(-1, -1)
      CALL xmap (lower, upper, zmin, zmax)

      CALL xaxnsz ( axlbsiz*(zmax-zmin)*lblmag )

      CALL xqmap(x1,x2,y1,y2)
      CALL xchsiz(0.03*(y2-y1)*lblmag)
      CALL xchori(0.0)

      IF( .NOT.multiprof ) THEN
        lchar = LEN( xcaptn)
        ch = xcaptn
        WRITE(ch(lchar+1:lchar+33), '(a,f13.3,a,f13.3,a)')              &
            ' at (',xpnt(ip),',',ypnt(ip),')'
        lchar = lchar+33
        CALL strmin(ch(1:lchar), lchar)
        CALL xcharc((x1+x2)*0.5, y1-(y2-y1)*0.10, ch(1:lchar))
      ELSE
        CALL xcharc((x1+x2)*0.5, y1-(y2-y1)*0.10, xcaptn )
      END IF
!
!-----------------------------------------------------------------------
!
!  Check if the points lie on one side of the axis. If the points are
!  all positive, draw the y-axis on the left border, if all points are
!  negative, draw the y-axis on the right border. If points lie on both
!  sides, draw the y-axis through x=0.0.
!
!-----------------------------------------------------------------------
!
      temp = lower * upper

      IF (temp > 0.0) THEN

        IF (lower > 0.0) THEN
          CALL xaxes(lower,0.0,zmin,0.0)
          CALL xchori(90.0)
          CALL xcharc(x1-0.12*(x2-x1), (y1+y2)*0.5, ycaptn)
        ELSE
          CALL xaxant(-1, 1)
          CALL xaxtik(1, -1)
          CALL xaxes(upper,0.0,zmin,0.0)
          CALL xchori(90.0)
          CALL xcharc(x1-0.05*(x2-x1), (y1+y2)*0.5, ycaptn)
        END IF

      ELSE
        CALL xaxes(0.0,0.0,zmin,0.0)
        CALL xchori(90.0)
        CALL xcharc(x1-0.10*(x2-x1), (y1+y2)*0.5, ycaptn)
      END IF

    END IF

    CALL xchori(0.0)
    CALL xbordr
    CALL xfull
!
!-----------------------------------------------------------------------
!
!  The first plot is labeled `A'. The subsequent plots will be `B'...
!
!-----------------------------------------------------------------------
!
    IF( multiprof ) THEN
      CALL xlbon
      CALL xlabel(CHAR(64+ip))

      ch(1:1) =  CHAR(64+ip)
      ch(2:2) =  ' '
      lchar = 2
      WRITE(ch(lchar+1:lchar+33), '(a,f13.2,a,f13.2,a)')                &
          ' at (',xpnt(ip),',',ypnt(ip),')'
      lchar = lchar+33
      CALL strmin(ch(1:lchar), lchar)

      CALL xqmap(x1,x2,y1,y2)
      CALL xchsiz(0.025*(y2-y1)*lblmag)
      CALL xchori(0.0)
      CALL xcharl(x1+(x2-x1)*0.03, y2-(y2-y1)*(0.03+0.035*ip),          &
          ch(1:lchar))

      CALL xlbsiz( ctrlbsiz*(y2-y1)*lblmag )
    ELSE
      CALL xlboff
    END IF

    CALL xwindw(lower, upper, zmin, zmax)

    CALL xqmap(x1,x2,y1,y2)

    CALL xcurve(profil(kbgn,ip),height(kbgn,ip),kend-kbgn+1,0)
    CALL xwdwof

  END DO

  RETURN
END SUBROUTINE vprofil


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


SUBROUTINE spltpara(inc,MIN,MAX,ovr,hlf,zro,col1,col2,pltvar),4

!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!  Set some parameters for one plot.
!
!-----------------------------------------------------------------------
!
!  AUTHOR:
!    Min Zou (3/2/98)
!
!  Modification history:
!
!-----------------------------------------------------------------------
!
!  INPUT:
!    inc      interval of the contour
!    min      the minimum value for the contour
!    max      the maximum valur foj the contour
!    ovr      overlay option
!    hlf      the contour highlight frequency
!    zro      define the attributes of zero contours
!    col1     the start color index for contour
!    col2     the end color index for contour
!    pltvar   the plot name
!    len      the length of pltvar
!
!-----------------------------------------------------------------------
!
  INTEGER :: ovr,hlf,zro,col1,col2
  REAL :: inc, MIN, MAX
  CHARACTER (LEN=12) :: pltvar

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

  CALL ctrinc( inc, MIN, MAX )
  CALL overlay(ovr)
  CALL xhlfrq (hlf)
  CALL xczero (zro)
  CALL ctrcol( col1,col2)
  CALL varplt( pltvar )

  RETURN
END SUBROUTINE spltpara


SUBROUTINE fillmissval (m, n, xl, xr, yb,yt ) 1

  REAL :: x1,x2,y1,y2
  REAL :: xra(4), yra(4)
  INTEGER :: missval_colind, missfill_opt    ! miss value color index
  COMMON /multi_value/ missfill_opt,missval_colind

  x1 = xl + (xr-xl)/REAL(m)*0.5
  x2 = xr - (xr-xl)/REAL(m)*0.5
  y1 = yb + (yt-yb)/REAL(n)*0.5
  y2 = yt - (yt-yb)/REAL(n)*0.5

  xra(1) = x1
  xra(2) = x2
  xra(3) = x2
  xra(4) = x1
  yra(1) = y1
  yra(2) = y1
  yra(3) = y2
  yra(4) = y2

  CALL xcolor(missval_colind)
  CALL xfilarea(xra, yra, 4)

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


SUBROUTINE hintrp(nx,ny,nz,a3din,z3d,zlevel, a2dout) 7
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Interpolate a 3-D array to horizontal level z=zlevel.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Ming Xue
!  Based on original SECTHRZ.
!  12/10/98.
!
!  MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
!  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
!
!    a3din    3-d input array
!    z3d      z-coordinate of data in a3din
!    zlevel   Level to which data is interpolated.
!
!  OUTPUT:
!    a2dout   2-d output array interpolated to zlevel
!
!-----------------------------------------------------------------------
!
!  Parameters of output
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
  INTEGER :: nx,ny,nz

  REAL :: a3din(nx,ny,nz) ! 3-d input array
  REAL :: z3d  (nx,ny,nz) ! z-coordinate of data in a3din
  REAL :: zlevel          ! Level to which data is interpolated.

  REAL :: a2dout(nx,ny)   ! 2-d output array interpolated to zlevel

  INTEGER :: i,j,k
!
!-----------------------------------------------------------------------
!
!  Find index for interpolation
!
!-----------------------------------------------------------------------
!
  DO i=1,nx-1
    DO j=1,ny-1
      IF(zlevel <= z3d(i,j,1)) GO TO 11
      IF(zlevel >= z3d(i,j,nz-1)) GO TO 12
      DO k=2,nz-2
        IF(zlevel >= z3d(i,j,k).AND.zlevel < z3d(i,j,k+1)) GO TO 15
      END DO

      11    k=1
      GO TO 15
      12    k=nz-1
      GO TO 15

      15    a2dout(i,j)=a3din(i,j,k)+(a3din(i,j,k+1)-a3din(i,j,k))*     &
                        (zlevel-z3d(i,j,k))/(z3d(i,j,k+1)-z3d(i,j,k))

!-----------------------------------------------------------------------
!
!  If the data point is below the ground level, set the
!  data value to the missing value.
!
!-----------------------------------------------------------------------

      IF( zlevel < z3d(i,j,2)   ) a2dout(i,j) = -9999.0
      IF( zlevel > z3d(i,j,nz-1)) a2dout(i,j) = -9999.0

    END DO
  END DO

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


SUBROUTINE hintrp1(nx,ny,nz, kbgn,kend,a3din,z3d,zlevel, a2dout) 8
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Interpolate a 3-D array to horizontal level z=zlevel.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Ming Xue
!  Based on original SECTHRZ.
!  12/10/98.
!
!  MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
!  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
!    kbgn
!    kend
!
!    a3din    3-d input array
!    z3d      z-coordinate of data in a3din
!    zlevel   Level to which data is interpolated.
!
!  OUTPUT:
!    a2dout   2-d output array interpolated to zlevel
!
!-----------------------------------------------------------------------
!
!  Parameters of output
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
  INTEGER :: nx,ny,nz
  INTEGER :: kbgn, kend

  REAL :: a3din(nx,ny,nz) ! 3-d input array
  REAL :: z3d  (nx,ny,nz) ! z-coordinate of data in a3din
  REAL :: zlevel          ! Level to which data is interpolated.

  REAL :: a2dout(nx,ny)   ! 2-d output array interpolated to zlevel

  INTEGER :: i,j,k
!
!-----------------------------------------------------------------------
!
!  Find index for interpolation
!
!-----------------------------------------------------------------------
!
  DO i=1,nx-1
    DO j=1,ny-1
      IF(zlevel <= z3d(i,j,kbgn)) GO TO 11
      IF(zlevel >= z3d(i,j,kend)) GO TO 12
      DO k=kbgn,kend-1
        IF(zlevel >= z3d(i,j,k).AND.zlevel < z3d(i,j,k+1)) GO TO 15
      END DO

      11    k=kbgn
      GO TO 15
      12    k=kend-1
      GO TO 15

      15    a2dout(i,j)=a3din(i,j,k)+(a3din(i,j,k+1)-a3din(i,j,k))*     &
                        (zlevel-z3d(i,j,k))/(z3d(i,j,k+1)-z3d(i,j,k))

!-----------------------------------------------------------------------
!
!  If the data point is below the ground level, set the
!  data value to the missing value.
!
!-----------------------------------------------------------------------

      IF( zlevel < z3d(i,j,kbgn) ) a2dout(i,j) = -9999.0
      IF( zlevel > z3d(i,j,kend) ) a2dout(i,j) = -9999.0

    END DO
  END DO

  RETURN
END SUBROUTINE hintrp1



SUBROUTINE indxbnds(xc,yc,zpc,zpsoilc,nx,ny,nz,nzsoil,                  & 1
           xbgn,xend,ybgn,yend,zbgn,zend,zsoilbgn,zsoilend,             &
           ibgn,iend,jbgn,jend,kbgn,kend,ksoilbgn,ksoilend)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!  Return index bounds of the domain to be plotted
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
  INTEGER :: nx,ny,nz
  INTEGER :: nzsoil
  REAL :: xc   (nx,ny,nz)   ! x-coor of sacalar point (km)
  REAL :: yc   (nx,ny,nz)   ! y-coor of sacalar point (km)
  REAL :: zpc  (nx,ny,nz)   ! z-coor of sacalar point in physical
                            ! space (km)
  REAL :: zpsoilc(nx,ny,nzsoil)   ! z-coor of sacalar point in physical
                            ! space (m) for soil model
  REAL :: xbgn,xend,ybgn,yend,zbgn,zend,zsoilbgn,zsoilend
  INTEGER :: ibgn,iend,jbgn,jend,kbgn,kend,ksoilbgn,ksoilend
  INTEGER :: i,j,k

  IF(xbgn /= xend) THEN
    DO i = 1,nx-1
      IF(xc(i,2,2) >= xbgn) EXIT
    END DO
    ibgn=MAX(i,1)

    DO i = nx-1,1,-1
      IF(xc(i,2,2) <= xend ) EXIT
    END DO
    iend=MIN(i,nx-1)

  ELSE

    ibgn = 2
    iend = nx- 2
    xbgn = xc(ibgn,2,2)
    xend = xc(iend,2,2)

  END IF

  IF(ybgn /= yend) THEN
    DO j = 1,ny-1
      IF(yc(2,j,2) >= ybgn) EXIT
    END DO
    jbgn=MAX(j,1)

    DO j = ny-1,1,-1
      IF(yc(2,j,2) <= yend) EXIT
    END DO
    jend=MIN(j,ny-1)

  ELSE

    jbgn = 2
    jend = ny-2
    ybgn = yc(2,jbgn,2)
    yend = yc(2,jend,2)

  END IF

  IF(zbgn /= zend) THEN

    kend = 2
    DO k = 2,nz-1

      DO j = jbgn,jend 
        DO i = ibgn,iend 
          IF(zpc(i,j,k) < zend) THEN
            kend=k
            GO TO 225
          END IF
        END DO
      END DO

      GO TO 235
225   CONTINUE

    END DO

235 kend = MIN(kend, nz-1)

    kbgn= nz-1
    DO k = nz-1,2,-1

      DO j = jbgn,jend 
        DO i = ibgn,iend 
          IF(zpc(i,j,k) > zbgn) THEN
            kbgn=k
            GO TO 250
          END IF
        END DO
      END DO

      GO TO 245
250   CONTINUE

    END DO

245 kbgn = MAX(kbgn,2)

  ELSE

    kbgn = 2
    kend = nz-2

  END IF

  IF(zsoilbgn /= zsoilend) THEN
!
! 05/31/2002 Zuwen He
!
! Note: k is 1 at the surface in the soil model, 
!       and k increase when zpsoilc decrease.
!       zpsoilc=zpsoil(k)-zpsoil(1) < 0. 
!
    ksoilend = 1
    DO k = 1,nzsoil

      DO j = jbgn,jend 
        DO i = ibgn,iend 
          IF(zpsoilc(i,j,k) > zsoilend) THEN
            ksoilend=k
            GO TO 325
          END IF
        END DO
      END DO

      GO TO 335
325   CONTINUE

    END DO

335 ksoilend = MIN(ksoilend+1, nzsoil)

    ksoilbgn= nzsoil
    DO k = nzsoil,1,-1

      DO j = jbgn,jend 
        DO i = ibgn,iend 
          IF(zpsoilc(i,j,k) < zsoilbgn) THEN
            ksoilbgn=k
            GO TO 350
          END IF
        END DO
      END DO

      GO TO 345
350   CONTINUE

    END DO

345 ksoilbgn = MAX(ksoilbgn-1,1)

  ELSE

    ksoilbgn = 1
    ksoilend = nzsoil

  END IF

  WRITE(6,'(/1x,a,i3,a,i3)') 'ibgn =',ibgn,', iend =',iend
  IF(iend < ibgn) THEN
    WRITE(6,'(1x,a,/1x,a)')                                             &
        'iend was found smaller than ibgn. Check the input',            &
        'domain bounds in x direction. Program stopped.'
    STOP
  END IF

  WRITE(6,'(1x,a,i3,a,i3)') 'jbgn =',jbgn,', jend =',jend
  IF(jend < jbgn) THEN
    WRITE(6,'(1x,a,/1x,a)')                                             &
        'jend was found smaller than jbgn. Check the input',            &
        'domain bounds in y direction. Program stopped.'
    STOP
  END IF

  WRITE(6,'(1x,a,i3,a,i3)') 'kbgn =',kbgn,', kend =',kend
  IF(kend < kbgn) THEN
    WRITE(6,'(1x,a,/1x,a)')                                             &
        'kend was found smaller than kbgn. Check the input',            &
        'domain bounds in z direction. Program stopped.'
    STOP
  END IF

  WRITE(6,'(1x,a,i3,a,i3)') 'ksoilbgn =',ksoilbgn,', ksoilend =',ksoilend
  IF(ksoilend < ksoilbgn) THEN
    WRITE(6,'(1x,a,/1x,a)')                                             &
        'ksoilend was found smaller than ksoilbgn. Check the input',            &
        'domain bounds in zpsoil direction. Program stopped.'
    STOP
  END IF

  RETURN
END SUBROUTINE indxbnds


SUBROUTINE ctrsetup(zinc,zminc,zmaxc,                                   & 82,4
           zovr,zhlf,zzro,zcol1,zcol2,zlabel)

  IMPLICIT NONE
  REAL :: zinc,zminc,zmaxc
  INTEGER :: zovr,zhlf,zzro,zcol1,zcol2
  CHARACTER (LEN=*) :: zlabel

  CALL ctrinc (zinc,zminc,zmaxc )
  CALL overlay(zovr)
  CALL xhlfrq (zhlf)
  CALL xczero (zzro)
  CALL ctrcol (zcol1,zcol2 )
  CALL varplt (zlabel)

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


SUBROUTINE plttrn(hterain,x,y,m,n,slicopt) 2,1
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Generate terrain contours
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Ming Xue
!
!  MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!    hterain  2-D terrain data to contour
!    x        x coordinate of grid points in plot space (over on page)
!    y        y coordinate of grid points in plot space (up on page)
!    m        first dimension
!    n        second dimension
!
!    slicopt  slice orientation indicator
!       slicopt = 1, x-y slice of u,v at z index kslice is plotted.
!       slicopt = 2, x-z slice of u,w at y index jslice is plotted.
!       slicopt = 3, y-z slice of v,w at x index islice is plotted.
!       slicopt = 4, x-y slice of u,v at z index islice is plotted.
!       slicopt = 5, xy-z cross section of wind islice is plotted.
!       slicopt = 6, data field on constant p-level is plotted.
!       slicopt = 0, all of the three slices above are plotted.
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
  INTEGER :: m,n
  REAL :: hterain(m,n)
  REAL :: x(m,n)
  REAL :: y(m,n)
  INTEGER :: slicopt
!
!-----------------------------------------------------------------------
!
!  Plotting control common blocks
!
!-----------------------------------------------------------------------
!
  REAL :: ctinc,ctmin,ctmax,vtunt  ! contour interval and vector unit
  COMMON /incunt/ ctinc,ctmin,ctmax,vtunt
!
  INTEGER :: icolor,icolor1,lbcolor,trcolor                ! required color
  COMMON /recolor/icolor,icolor1,lbcolor,trcolor
!
  REAL :: ztmin,ztmax
  INTEGER :: ovrtrn ,trnplt       ! overlay terrain option (0/1)
  REAL :: trninc,trnmin, trnmax   ! terrain interval minimum, maximum
  COMMON /trnpar/ trnplt,ovrtrn,trninc,trnmin, trnmax,ztmin,ztmax
  REAL :: zlevel
  COMMON /sliceh/zlevel
  INTEGER :: col_table,pcolbar
  COMMON /coltable/col_table,pcolbar
!
!-----------------------------------------------------------------------
!
!  Misc. local variables
!
!-----------------------------------------------------------------------

  REAL :: cl(500)
  INTEGER :: ncl

  REAL :: z02,xl,xr,yt,yb,xfinc
  INTEGER :: mode1

  INTEGER, ALLOCATABLE :: iwrk(:)
  REAL   , ALLOCATABLE :: xwk(:),ywk(:)
  INTEGER :: istatus
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  ALLOCATE(iwrk(m*n),STAT=istatus) 
  ALLOCATE(xwk(m*n),ywk(m*n),STAT=istatus) 

  IF( ovrtrn == 0 .OR. ztmax-ztmin < 1.0E-20  ) RETURN
!
!-----------------------------------------------------------------------
!
!  Overlay terrain contour if required in x-y level
!  or Plot terrain outline in slice zlevel
!
!-----------------------------------------------------------------------
!
  CALL xqmap(xl,xr,yb,yt)

  cl(1)=0.0
  IF(slicopt == 1 .OR. slicopt == 8 .OR. slicopt == 9)  THEN

    CALL ctrinc( trninc,trnmin, trnmax )
    IF( trninc == 0.0) THEN
      cl(2)=cl(1)+ xfinc(ztmax-ztmin)/2
      IF(cl(2)-cl(1) == 0.0) cl(2)=cl(1)+1.0
      mode1=1
      CALL xnctrs(6,18)
    ELSE
      cl(2)=cl(1)+trninc
      CALL xnctrs(1,300)
      IF(ztmin == 0.0 .AND. ztmax == 0.0) THEN
        mode1=1
      ELSE
        mode1=3
      END IF
    END IF

    CALL xctrlim(ctmin,ctmax)
    IF (trnplt == 1) THEN
      CALL xthick(2)
      CALL xctrclr(trcolor, trcolor)
      IF(mode1 == 3) THEN
        ncl=(ztmax-ztmin)/trninc+1
        cl(1)=ztmin
        cl(2)=cl(1)+trninc
      END IF
      CALL xconta(hterain,x,y,iwrk,m,m,n,cl,ncl,mode1)
    ELSE IF (trnplt == 2) THEN
      CALL xctrclr(icolor, icolor1)
      CALL xcolfil(hterain,x,y,iwrk,xwk,ywk,m,m,n,cl,ncl,mode1)
      CALL xchsiz(0.025*(yt-yb))
      CALL xcpalet(pcolbar)
    ELSE IF (trnplt == 4) THEN
      CALL xctrclr(icolor, icolor1)
      CALL xconta(hterain,x,y,iwrk,m,m,n,cl,ncl,mode1)
    END IF
  ELSE IF(slicopt == 4.OR.slicopt == 6.OR.slicopt == 7) THEN
    CALL xcolor(trcolor)
    z02=zlevel*1000.
    CALL xthick(2)
    CALL xcontr(hterain,x,y,iwrk,m,m,n,z02)
    CALL xthick(1)
  END IF

  DEALLOCATE(iwrk)
  DEALLOCATE(xwk,ywk)

  RETURN
END SUBROUTINE plttrn

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


SUBROUTINE pltaxes(slicopt,dx,dy) 3,3
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!-----------------------------------------------------------------------
!
!  AUTHOR:  M. Xue
!
!  MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
!
!  INPUT:
!    slicopt slice orientation indicator
!         = 1, x-y slice of at k=kslice is plotted.
!         = 2, x-z slice of at j=jslice is plotted.
!         = 3, y-z slice of at i=islice is plotted.
!         = 4, horizontal slice at z index islice is plotted.
!         = 5, xy-z cross section of wind islice is plotted.
!         = 6, data field on constant p-level is plotted.
!         = 0, all of the three slices above are plotted.
!    dx   Spacing between the x-axis tick marks
!    dy   Spacing between the y-axis tick marks
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
  REAL :: dx,dy
  INTEGER :: slicopt
!
!-----------------------------------------------------------------------
!
!  Plotting control common blocks
!
!-----------------------------------------------------------------------
!
  INTEGER :: layover
  COMMON /laypar/ layover

  INTEGER :: presaxis_no
  REAL :: pres_val(20), pres_z(20)
  COMMON /pressbar_par/presaxis_no,pres_val,pres_z

  REAL :: lblmag, ctrlbsiz, axlbsiz
  COMMON /labmag/ lblmag, ctrlbsiz, axlbsiz

  INTEGER :: timeovr
  COMMON /timover/ timeovr

  REAL :: x101, y101, x102,y102
  COMMON /slicev1/x101, y101, x102,y102

  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

  INTEGER :: icolor,icolor1,lbcolor,trcolor       ! required color
  COMMON /recolor/icolor,icolor1,lbcolor,trcolor
!
!-----------------------------------------------------------------------
!
!  Misc. local variables
!
!-----------------------------------------------------------------------
!
  REAL :: axnmag
  REAL :: xl,xr,yb,yt,pl,pr,pb,pt

  REAL :: xtem1, xtem2  !local temporary variable
  REAL :: x1,x2, y1,y2, xstep, ystep, xmstep, ymstep
  CHARACTER (LEN=15) :: ylabel
  CHARACTER (LEN=15) :: xlabel
  INTEGER :: LEN
!
!-----------------------------------------------------------------------
!
!  Set-up variables for tick marks and draw axes
!
!-----------------------------------------------------------------------
!
  CALL xqmap(xl,xr, yb,yt)

  CALL setcords(xl,xr,yb,yt,dx,dy, slicopt,                             &
       x1,x2,y1,y2,xlabel,ylabel,xstep,ystep,xmstep,ymstep)

  CALL xqpspc( pl, pr, pb, pt)
  axnmag = axlbsiz*MIN(pt-pb, pr-pl)*lblmag

  CALL xaxnmg( axnmag )

  IF(slicopt == 5) THEN
    IF( ABS(y101-y102) <= 1.0E-3 ) THEN
      xtem1 = x101
      xtem2 = x102
    ELSE IF(ABS(x101-x102) <= 1.0E-3 ) THEN
      xtem1 = y101
      xtem2 = y102
    ELSE
      xtem1 = SQRT(x101*x101 + y101*y101)
      xtem2 = xtem1 + SQRT( (y102-y101)*(y102-y101) +                   &
                            (x102-x101)*(x102-x101) )
    END IF
  ELSE
    xtem1 = x1
    xtem2 = x2
  END IF

  CALL xmap(xtem1,xtem2,y1,y2)
  IF( layover == 0) THEN
    CALL xaxsor(0.0, 0.0)
!    call xthick(2)
    CALL xaxsca1( xtem1,xtem2,xstep,xmstep, y1,y2,ystep,ymstep )
    CALL xthick(1)
  END IF

!
!  Plot pressure axis
!
  CALL xqpspc(pl, pr, pb, pt)
  IF(presaxis_no > 0 .AND.timeovr == 0 .AND.                            &
     (slicopt == 2 .OR. slicopt == 3 .OR. slicopt == 5 .OR.  & 
      slicopt == 10 .OR. slicopt ==11) ) THEN
    x1 = pl - (pr-pl)*0.25
    x2 = pl
    y1 = pb
    y2 = pt
    CALL xpspac(x1,x2,y1,y2)
    y1 = yb
    y2 = yt
    CALL xmap(x1,x2,y1,y2)
    CALL xaxfmt('(I4)')
    CALL xyaxis(x1+0.40*(x2-x1),pres_z,pres_val,presaxis_no)
    CALL xchori(90.)
    CALL xcharc(x1,yb+(yt-yb)*0.5 ,'Pressure(mb)')
    CALL xchori(0.)
  END IF
!
!  Restore the original plotting scape
!
  CALL xpspac( pl, pr, pb, pt)
  CALL xmap(xl,xr, yb,yt)

  IF(layover > 1) THEN
    CALL xchsiz( 0.018*(yt-yb) * lblmag )
  ELSE
    CALL xchsiz( 0.020*(yt-yb) * lblmag )
  END IF

  IF(lbaxis == 1 .AND. timeovr == 0) THEN
    CALL xcolor(lbcolor)
    LEN=15
    CALL strmin(xlabel,LEN)
    CALL xcharc( xl+(xr-xl)*0.5,yb-0.08*(yt-yb),xlabel(1:LEN))
    LEN=15
    CALL strmin(ylabel,LEN)
    CALL xchori(90.)
    CALL xcharc(xl-0.10*(xr-xl),yb+(yt-yb)*0.5,ylabel(1:LEN))
    CALL xchori(0.)
  END IF

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


SUBROUTINE pltextra(slicopt, pltopt) 3,2
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Plot extra things such as map, boxes, polygons and stations
!  in a 2D plot
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Ming Xue
!
!  MODIFICATION HISTORY:
!
!  6/08/92 (K. Brewster)
!  Added full documentation.
!
!  8/28/94 (M. Zou)
!  Added color routing , overlay terrain.
!
!  1/24/96 (J. Zong and M. Xue)
!  Fixed a problem related to finding the minimum and maximum of the
!  2D array, a, when there exist missing data. Initial min. and max.
!  should be set to values other than the missing value, -9999.0.
!
!-----------------------------------------------------------------------
!
!  INPUT:
!    slicopt     slice orientation indicator
!             = 1, x-y slice of at k=kslice is plotted.
!             = 2, x-z slice of at j=jslice is plotted.
!             = 3, y-z slice of at i=islice is plotted.
!             = 4, horizontal slice at z index islice is plotted.
!             = 5, xy-z cross section of wind islice is plotted.
!             = 6, data field on constant p-level is plotted.
!             = 0, all of the three slices above are plotted.
!    plot
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
  INTEGER :: slicopt, pltopt

  INCLUDE 'arpsplt.inc'
!
!-----------------------------------------------------------------------
!
!  Plotting control common blocks
!
!-----------------------------------------------------------------------
!
  INTEGER :: ovrstaopt
  INTEGER :: ovrstam,staset,ovrstan,ovrstav,stacol,markprio,wrtstax
  INTEGER :: nsta_typ,sta_typ(10),sta_marktyp(10),sta_markcol(10)
  REAL :: sta_marksz(10),wrtstad
  COMMON /sta_par/ ovrstaopt,ovrstam,staset,ovrstan,ovrstav,stacol,     &
         markprio,nsta_typ,sta_typ,sta_marktyp,                         &
         sta_markcol,sta_marksz,wrtstax,wrtstad

  INTEGER :: icolor,icolor1,lbcolor,trcolor                ! required color
  COMMON /recolor/icolor,icolor1,lbcolor,trcolor
!
!-----------------------------------------------------------------------
!
!  Misc. local variables
!
!-----------------------------------------------------------------------
!
  INTEGER :: timeovr
  COMMON /timover/ timeovr

  INTEGER :: ovrmap
  COMMON /mappar / ovrmap

  INTEGER :: number_of_boxes,boxcol
  REAL :: bx1(10),bx2(10),by1(10),by2(10)
  COMMON /boxesopt/number_of_boxes,boxcol,bx1,bx2,by1,by2

  INTEGER :: num_of_verts
  INTEGER :: number_of_polys,polycol
  REAL :: vertx(max_verts,max_polys),verty(max_verts,max_polys)
  COMMON /polysopt/number_of_polys,polycol,vertx,verty

  REAL :: lblmag, ctrlbsiz, axlbsiz
  COMMON /labmag/ lblmag, ctrlbsiz, axlbsiz

  INTEGER :: nunit
  INTEGER :: i,j
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!
!-----------------------------------------------------------------------
!
!  Plot map boundaries.
!
!-----------------------------------------------------------------------
!
  CALL xcolor(lbcolor)

  IF((slicopt == 1 .OR. slicopt == 4 .OR. slicopt == 6 .OR.             & 
      slicopt == 7 .OR. slicopt == 8 .OR. slicopt == 9)                 &
        .AND.ovrmap == 1                                                &
        .AND.(timeovr == 0 .OR. (timeovr == 1 .AND. pltopt == 2) ))THEN
    CALL getunit(nunit)
    CALL drawmap(nunit)


    CLOSE (UNIT=nunit)
    CALL retunit(nunit)
    CALL xthick(1)
  END IF
!
!-----------------------------------------------------------------------
!
!  Draw boxes
!
!-----------------------------------------------------------------------
!
  IF(number_of_boxes /= 0 .AND.  & 
     (slicopt == 1 .OR. slicopt == 4 .OR. slicopt == 6 .OR.  & 
      slicopt == 7 .OR. slicopt == 8 .OR. slicopt == 9)                &
        .AND. timeovr == 0 ) THEN
    CALL xthick(1)
    CALL xcolor(boxcol)
    CALL xbrokn(6,3,6,3)
    DO i=1,number_of_boxes
      CALL xbox(bx1(i),bx2(i),by1(i),by2(i))
    END DO
    CALL xthick(1)
    CALL xfull
  END IF
!
!-----------------------------------------------------------------------
!
!  Draw polylines
!
!-----------------------------------------------------------------------
!
  IF(number_of_polys /= 0 .AND.  & 
     (slicopt == 1 .OR. slicopt == 4 .OR. slicopt == 6 .OR.  & 
      slicopt == 7 .OR. slicopt == 8 .OR. slicopt == 9)                &
        .AND. timeovr == 0 ) THEN
    CALL xthick(2)
    CALL xcolor(polycol)
!    CALL xbrokn(6,3,6,3)
    DO j=1,number_of_polys
      num_of_verts=0
      DO i=1,max_verts
        IF(vertx(i,j) /= -9999. .AND. verty(i,j) /= -9999.)             &
                  num_of_verts = num_of_verts +1
      END DO
      IF(num_of_verts /= 0 ) CALL xcurve(vertx(1,j),verty(1,j),num_of_verts, 0)
    END DO
    CALL xthick(1)
    CALL xfull
  END IF

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


SUBROUTINE smooth9pmv( arr, nx,ny,ibgn,iend,jbgn,jend, tem1 ) 31
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!                                        1 2 1
!  Smooth a 2-D array by the filter of { 2 4 2 }
!                                        1 2 1
!
!-----------------------------------------------------------------------
!
!  AUTHOR:       Yuhe Liu
!
!  5/3/94
!
!  Modification History
!  8/20/1995 (M. Xue)
!  Fixed errors in the index bound of loops 100 and 200.
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!  nx       Number of grid points in the x-direction
!  ny       Number of grid points in the y-direction
!  ibgn     First index in x-direction in the soomthing region.
!  iend     Last  index in x-direction in the soomthing region.
!  jbgn     First index in j-direction in the soomthing region.
!  jend     Last  index in j-direction in the soomthing region.
!
!  arr    2-D array
!
!  OUTPUT:
!
!  arr    2-D array
!
!  TEMPORARY:
!
!  tem1     Temporary 2-D array
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
!
  INTEGER :: nx         ! Number of grid points in the x-direction
  INTEGER :: ny         ! Number of grid points in the y-direction
  INTEGER :: ibgn
  INTEGER :: iend
  INTEGER :: jbgn
  INTEGER :: jend
!
  REAL :: arr (nx,ny)   ! 2-D array
!
  REAL :: tem1(nx,ny)   ! Temporary array
!
!-----------------------------------------------------------------------
!
!  Misc. local variables:
!
!-----------------------------------------------------------------------
!
  INTEGER :: i,j,ip,im,jp,jm
  REAL :: wtf,mv
!
!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------
!

!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  wtf    = 1.0/16.0

  mv = -9999.0 ! missing value flag

  DO i=1,nx
    DO j=1,ny
      IF( ABS(arr(i,j)-mv) <= 0.1) arr(i,j)=mv
    END DO
  END DO

  DO j = jbgn,jend
    DO i = ibgn,iend
      ip=MIN(nx,i+1)
      im=MAX( 1,i-1)
      jp=MIN(ny,j+1)
      jm=MAX( 1,j-1)

      tem1(i,j) = wtf                                                   &
          * (    arr(im,jm) + 2.*arr(i,jm) +    arr(ip,jm)              &
          + 2.*arr(im,j ) + 4.*arr(i,j ) + 2.*arr(ip,j )                &
          +    arr(im,jp) + 2.*arr(i,jp) +    arr(ip,jp))

      IF(arr(im,jm) == mv.OR.arr(i,jm) == mv.OR.arr(ip,jm) == mv.OR.    &
            arr(im,j ) == mv.OR.arr(i,j ) == mv.OR.arr(ip,j ) == mv.OR. &
            arr(im,jp) == mv.OR.arr(i,jp) == mv.OR.arr(ip,jp) == mv)THEN
        tem1(i,j)=mv
      END IF

    END DO
  END DO

  DO j = jbgn,jend
    DO i = ibgn,iend
      arr(i,j) = tem1(i,j)
    END DO
  END DO

  RETURN
END SUBROUTINE smooth9pmv


SUBROUTINE buoycy_plt(nx,ny,nz,ptprt,pprt,qv,qc,qr,qi,qs,qh,            & 1
           ptbar,pbar,rhobar,qvbar, wbuoy, tem1)
!
!-----------------------------------------------------------------------
!
!     PURPOSE:
!
!     Calculate the total buoyancy including liquid and solid water
!     loading.
!
!-----------------------------------------------------------------------
!
!     AUTHOR: Ming Xue
!     10/10/91.
!
!     MODIFICATION HISTORY:
!
!     5/05/92 (M. Xue)
!     Added full documentation.
!
!     3/10/93 (M. Xue)
!     The buoyancy term is reformulated. The previous formula was
!     in error. The water loading was calculated wrong, resulting in
!     a value of the water loading that is typically an order of
!     magnitude too small.
!
!     3/25/94 (G. Bassett & M. Xue)
!     The buoyancy terms are reformulated for better numerical accuracy.
!     Instead of storing numbers which had the form (1+eps)*(1+eps1)
!     (eps << 1 and eps1 <<1), terms were expanded out, and most of the
!     high order terms neglected, except for the second order terms
!     in ptprt, pprt and qvbar.
!
!     9/10/94 (D. Weber & Y. Lu)
!     Cleaned up documentation.
!
!     6/21/95 (Alan Shapiro)
!     Fixed bug involving missing qvpert term in buoyancy formulation.
!
!     10/15/97 (Donghai Wang)
!     Added a new option for including the second order terms.
!
!     11/05/97 (D. Weber)
!     Changed lower loop bounds in DO LOOP 400 for computing the
!     buoyancy term from k=3,nz-2 to k=2,nz-1.  Level k=2 data will be
!     used in the hydrostatic pprt lower boundary condition (removed
!     DO LOOP 410 used to set wbuoy = 0.0 for k= 2 and nz-1).
!
!-----------------------------------------------------------------------
!
!     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 direction.
!
!       ptprt    Perturbation potential temperature at a time level (K)
!       pprt     Perturbation pressure at a given time level (Pascal)
!       qv       Water vapor specific humidity at a given time level
!                (kg/kg)
!       qc       Cloud water mixing ratio at a given time level (kg/kg)
!       qr       Rainwater mixing ratio at a given time level (kg/kg)
!       qi       Cloud ice mixing ratio at a given time level (kg/kg)
!       qs       Snow mixing ratio at a given time level (kg/kg)
!       qh       Hail mixing ratio at a given time level (kg/kg)
!
!       ptbar    Base state potential temperature (K)
!       pbar     Base state pressure (Pascal)
!       rhobar   Base state density rhobar
!       qvbar    Base state water vapor specific humidity (kg/kg)
!
!     OUTPUT:
!
!       wbuoy    The total buoyancy force (kg/(m*s)**2)
!
!     WORK ARRAYS:
!
!       tem1     Temporary work array.
!
!-----------------------------------------------------------------------
!
!     Variable Declarations
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE

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

  REAL :: ptprt (nx,ny,nz)     ! Perturbation potential temperature
                               ! at a given time level (K)
  REAL :: pprt  (nx,ny,nz)     ! Perturbation pressure at a given time
                               ! level (Pascal)
  REAL :: qv    (nx,ny,nz)     ! Water vapor specific humidity (kg/kg)
  REAL :: qc    (nx,ny,nz)     ! Cloud water mixing ratio (kg/kg)
  REAL :: qr    (nx,ny,nz)     ! Rain water mixing ratio (kg/kg)
  REAL :: qi    (nx,ny,nz)     ! Cloud ice mixing ratio (kg/kg)
  REAL :: qs    (nx,ny,nz)     ! Snow mixing ratio (kg/kg)
  REAL :: qh    (nx,ny,nz)     ! Hail mixing ratio (kg/kg)

  REAL :: ptbar (nx,ny,nz)     ! Base state potential temperature (K)
  REAL :: pbar  (nx,ny,nz)     ! Base state pressure (Pascal).
  REAL :: rhobar(nx,ny,nz)     ! Base state density rhobar
  REAL :: qvbar (nx,ny,nz)     ! Base state water vapor specific
                               ! humidity(kg/kg)

  REAL :: wbuoy(nx,ny,nz)      ! Total buoyancy in w-eq. (kg/(m*s)**2)

  REAL :: tem1  (nx,ny,nz)     ! Temporary work array.

!
!-----------------------------------------------------------------------
!
!     Misc. local variables:
!
!-----------------------------------------------------------------------
!
  INTEGER :: i,j,k
  REAL :: g5
  REAL :: pttem,tema
!
!-----------------------------------------------------------------------
!
!     Include files:
!
!-----------------------------------------------------------------------
!
  INCLUDE 'globcst.inc'     ! Global model control parameters
  INCLUDE 'phycst.inc'      ! Physical constants
!
!-----------------------------------------------------------------------
!
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!     Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!
!-----------------------------------------------------------------------
!
!     The total buoyancy
!
!       wbuoy = rhobar*g ( ptprt/ptbar-pprt/(rhobar*csndsq)+
!       qvprt/(rddrv+qvbar)-(qvprt+qc+qr+qs+qi+qh)/(1+qvbar)
!       -(ptprt*ptprt)/(ptbar*ptbar)                        !2nd-order
!       +0.5*(ptprt*pprt)/(cpdcv*ptbar*pbar))               !2nd-order
!
!     and rddrv=rd/rv, cp, cv, rd and rv are defined in phycst.inc.
!
!     Here, the contribution from pprt (i.e., term pprt/(rhobar*csndsq))
!     is evaluated inside the small time steps, therefore wbuoy
!     does not include this part.
!
!     The contribution from ptprt is calculated inside the small time
!     steps if the potential temperature equation is solved inside
!     small time steps, i.e., if ptsmlstp=1.
!
!-----------------------------------------------------------------------
!
  tema = 1.0/cpdcv
  DO k=1,nz-1
    DO j=1,ny-1
      DO i=1,nx-1
        pttem = ptprt(i,j,k)/ptbar(i,j,k)
        tem1(i,j,k) = pttem*                                            &
            (1.0-pttem+0.5*pprt(i,j,k)*(tema/pbar(i,j,k)))
      END DO
    END DO
  END DO
!
!-----------------------------------------------------------------------
!
!     Add on the contributions to the buoyancy from the water vapor
!     content and the liquid and ice water loading.
!
!-----------------------------------------------------------------------
!

  DO k=1,nz-1
    DO j=1,ny-1
      DO i=1,nx-1
        tem1(i,j,k) = tem1(i,j,k)                                       &
            + (qv(i,j,k) - qvbar(i,j,k))/(rddrv + qvbar(i,j,k))         &
            - (qv(i,j,k) - qvbar(i,j,k) + qc(i,j,k) + qr(i,j,k) +       &
            qs(i,j,k) + qi(i,j,k) + qh(i,j,k))/(1 + qvbar(i,j,k))
      END DO
    END DO
  END DO
!
!-----------------------------------------------------------------------
!
!     Then the total buoyancy:
!
!       wbuoy = tem1 * rhobar * g
!
!     averged to the w-point on the staggered grid.
!
!-----------------------------------------------------------------------
!
  DO k=1,nz-1
    DO j=1,ny-1
      DO i=1,nx-1
        wbuoy(i,j,k)= tem1(i,j, k )*rhobar(i,j, k )*g
      END DO
    END DO
  END DO

  RETURN
END SUBROUTINE buoycy_plt