      PROGRAM ARPSPLT
c
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                PROGRAM ARPSPLT                       ######
c     ######                                                      ######
c     ######                Copyright (c) 1992                    ######
c     ######  Center for the Analysis and Prediction of Storms    ######
c     ######    University of Oklahoma.  All rights reserved.     ######
c     ######                                                      ######
c     ##################################################################
c     ##################################################################
c

c
c#######################################################################
c
c     This is a graphic analysis/plotting program for model ARPS 3.0 
c
c     It is based on Ming Xue's plotting package ZXPLOT.
c
c     ZXPLOT is an independent graphic package. 
c
c     ZXPLOT can be interfaced with NCAR graphics to produce metafile 
c     output consistent with NCAR graphics viewing facilities
c     or linked to the Postscript driver (pure fortran program)
c     to produce Postscript output. A few other drivers are also
c     available. 
c
c     It shares with the model the include file dims.inc 
c     for definition of dimensions and domain size. 
c
c     The input data format is flagged by parameter hinfmt .
c
c     hinfmt  = 1 for machine dependent unformatted binary dump,
c     hinfmt  = 2 for formatted ascii dump,
c     hinfmt  = 3 for NCSA HDF file format.
c     hinfmt  = 4 for machine dependent packed unformatted binary dump
c
c     Parameters grdin, basin, mstin and icein indicate
c     if the data set includes grid, base state variables, 
c     water variables qv, qc, qr  and icein variables qi, qs and qh.
c     The values of these parameters read in from the data sets 
c     themselves.
c
c#######################################################################
c
c     AUTHOR: Ming Xue, CAPS/OU. 
c       2/19/1992.     
c
c     MODIFICATION HISTORY:
C     4/1/92 version 2.1.
c       wind vector plotting calls added.
c
c     4/7/92
c       3-D wireframe plotting of w, qc+qr fields added 
c       (borrow NCARgraphic routines).
c
c     4/13/92, version 2.2
c       Modifed to handle multi-dataset.
c
c     4/14/92 Version 2.3.
c       plt3d and plt2d names changed to ctr3d and ctr2d.
c       Physical space and mapping definitions changed in ctr2d and vtr2d. 
c       Plotting of a portion of domain now possible. 
c       Function for overlaying of fields added.
c
c       Control over the contour interval via subroutine ctrinc added.
c       Control over the wind vector unit via subroutine vtrunt added.
c       The aspect ratio of the 2-d plots can now be modified.
c
c     5/05/92 (M. Xue) and 06/04/92 (K. Brewster)
c       Added full documentation.
c
c     6/20/92 K. Brewster
c       Cosmetic changes to the output.
c
c       Added runname variable to ascread.
c
c     6/29/92 M. Xue
c       Removed external statement for getpte (caused problem on CRAY).
c       Removed close(unit=...) statement before open(unit...), several
c       places (caused trouble on CRAY).
c
c       Vertical aspect ratio increased to 2.0
c
c       Added a parameter that controls the plotting layout, this change
c       affects the input file.
c
c     7/01/92 M. Xue
c       Plot comments on the wireframe plots.
c
c     8/02/92 Keith Brewster
c       Changed reading sequence to allow for a separate file containing
c       the grid and base state information.
c
c       Changed structure so that the files are not re-read for each loop
c       over mode (DO 130).
c
c       Removed start and end time parameters.  It is assumed that each file
c       contains data from a single time.  Times desired are controlled by
c       naming the files to plot (the end of the input file).
c
c       These changes affect the input file.
c
c     8/09/92 M. Xue
c       uprt , vprt , wprt , ptprt and qvprt  are read in from the data
c       instead of u, v, w, pt and qv.
c
c     8/16/92 M. Xue
c       Grouped data read into subroutine dtaread.
c
c     9/04/92 M. Xue
c       Reading of parameters grdin, basin, mstin and icein from parameter
c       input file was removed.  
c       They are now read in directly from the data set.  
c
c     9/20/92 M. Xue
c       parameter kmplot added to the input parameter reading list and 
c       subroutine call to plot km added.
c
c#######################################################################
c
c     DATA ARRAYS READ IN:
c
c       x        x coordinate of grid points in physical/comp. space (m)
c       y        y coordinate of grid points in physical/comp. space (m)
c       zc       z coordinate of grid points in computational space (km)
c       z        z coordinate of grid points in computational space (m)
c
c       uprt     x component of perturbation velocity (m/s)
c       vprt     y component of perturbation velocity (m/s)
c       wprt     vertical component of perturbation velocity in Cartesian
c                coordinates (m/s).
c
c       ptprt    perturbation potential temperature (K)
c       pprt     perturbation pressure (Pascal)
c
c       qvprt    perturbation water vapor mixing ratio (kg/kg)
c       qc       Cloud water mixing ratio (kg/kg)
c       qr       Rainwater mixing ratio (kg/kg)
c       qi       Cloud ice mixing ratio (kg/kg)
c       qs       Snow mixing ratio (kg/kg)
c       qh       Hail mixing ratio (kg/kg)
c       km       Turbulent mixing coefficient
c
c       ubar     Base state x velocity component (m/s)
c       vbar     Base state y velocity component (m/s)
c       wbar     Base state z velocity component (m/s)
c       ptbar    Base state potential temperature (K)
c       pbar     Base state pressure (Pascal)
c       rhobar   Base state air density (kg/m**3)
c       qvbar    Base state water vapor mixing ratio (kg/kg)
c
c     CALCULATED DATA ARRAYS:
c
c       u        x component of velocity (m/s)
c       v        y component of velocity (m/s)
c       w        z component of velocity (m/s)
c       pt       potential temperature (K)
c       qv       water vapor mixing ratio (kg/kg)
c
c     WORK ARRAYS:
c
c       tem1     Temporary work array.
c       tem2     Temporary work array.
c       tem3     Temporary work array.
c       tem4     Temporary work array.
c
c     (These arrays are defined and used locally (i.e. inside this
c      subroutine), they may also be passed into routines called by
c      this one. Exiting the call to this subroutine, these temporary
c      work arrays may be used for other purposes, and therefore their
c      contents may be overwritten. Please examine the usage of work 
c      arrays before you make any change to the code.)
c
c#######################################################################
c
c     Variable Declarations.
c
c#######################################################################
c
      implicit none 
c
c#######################################################################
c
c     Include files: 
c
c#######################################################################
c
c     This include file defines the dimension of the arrays
c     i.e., nx,ny and nz
c
      include 'dims.inc'
      include 'indtflg.inc'
      include 'globcst.inc'
c
c#######################################################################
c
c     Arrays to be read in:
c
c#######################################################################
c
      integer hinfmt ! input data format

      real x    (nx,ny,nz)
      real y    (nx,ny,nz)
      real zc   (nx,ny,nz)
      real z    (nx,ny,nz)

      real u    (nx,ny,nz)
      real v    (nx,ny,nz)
      real w    (nx,ny,nz) 
      real pt   (nx,ny,nz)
      real pprt (nx,ny,nz)
      real qv   (nx,ny,nz)
      real qc   (nx,ny,nz)
      real qr   (nx,ny,nz)
      real qi   (nx,ny,nz)
      real qs   (nx,ny,nz)
      real qh   (nx,ny,nz)
      real km   (nx,ny,nz)

      real ubar  (nx,ny,nz)
      real vbar  (nx,ny,nz)
      real wbar  (nx,ny,nz)
      real ptbar (nx,ny,nz)
      real rhobar(nx,ny,nz)
      real pbar  (nx,ny,nz)
      real qvbar (nx,ny,nz)
c
c#######################################################################
c
c    Perturbation winds and potential temperature
c
c#######################################################################
c
      real uprt  (nx,ny,nz)
      real vprt  (nx,ny,nz)
      real wprt  (nx,ny,nz)
      real ptprt (nx,ny,nz)
      real qvprt (nx,ny,nz)
c
c#######################################################################
c
c     Temporary work arrays
c
c#######################################################################
c
      real tem1(nx,ny,nz)
      real tem2(nx,ny,nz)
      real tem3(nx,ny,nz)
      real tem4(nx,ny,nz)

c
c#######################################################################
c
c     Plotting control parameters entered by user
c
c#######################################################################
c
      integer vtrplt, vtpplt
      integer uplot, vplot, wplot, ptplot, pplot
      integer qvplot,qcplot,qrplot,qiplot,qsplot,qhplot,kmplot,
     :        rhplot,pteplt
c
      integer upplot,vpplot,wpplot,ptpplt,ppplot,qvpplt
      integer w3dplt, q3dplt
      real wisosf,qisosf
c
c#######################################################################
c
c     Common blocks for plotting control parameters
c
c#######################################################################
c
      integer layover               ! set by OVERLAY, used in ctr2d,vtr2d
      common /laypar/ layover
c
      real ctinc,vtunt
      common /incunt/ ctinc, vtunt  ! contour interval and vector unit 
c
      real yxratio
      common /yratio/ yxratio       ! the scaling factor the y/x ratio. 
c
C#######################################################################
C
C     Misc. local variables
C
C#######################################################################
c
      character grdbasfn*80, filename*80, label*25
      integer iplot(3)
      integer ireply,iret,ireturn
      integer mode,lengbf,lenfil,level
      integer i,j,k,ibgn,iend,jbgn,jend,kbgn,kend,ist,jst,kst
      integer nxpic,nypic,islice,jslice,kslice,layout
      real time,tbgnplt,tendplt,angl,aspect
      real xr,yr,zr,x1,x2,y1,y2,z1,z2,xlimit,ylimit
      real factor
c
C#######################################################################
C
C     ZXPLOT Routines Called
C
C#######################################################################
c
      external xdevic
      external xdspac
      external xspace
      external xpmagn
      external xlbint
      external xczero
      external xhlfrq
      external wirfrm
      external xgrend
c
c#######################################################################
c
c     Other Routines Called
c
c#######################################################################
c
      integer dsfirst
      external dsfirst
      external filzero
c     external getpte
      external temper
      external satmr
      external slength
      external ctr3d
      external vtr3d
      external vtrunt
      external ctrinc
      external overlay
c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c
      layover = 0
      ctinc = 0.0
      vtunt = 0.0
      yxratio = 1.0
c
c#######################################################################
c
c     Obtain from user the control parameters for
c     water, ice, and base state variable reading.
c
c#######################################################################
c
      write(6,'(/ 16(/5x,a)//)')

     :'###############################################################',
     :'###############################################################',
     :'####                                                       ####',
     :'####                 Welcome to ARPSPLT                    ####',
     :'####                                                       ####',
     :'####       A graphic analysis program for model ARPS 3.0   ####',
     :'####                                                       ####',
     :'####                    Version 2.3                        ####',
     :'####                                                       ####',
     :'####            The graphic plotting is based              ####',
     :'####              on graphic package ZXPLOT.               ####',
     :'####                                                       ####',
     :'####               By Ming Xue CAPS/OU                     ####',
     :'####                                                       ####',
     :'###############################################################',
     :'###############################################################'

c
      write(6,'(/a/)')' Please tell the program about the data sets.'
      write(6,'(a)')' Please input the data format flag value 1/2/3/4 ' 
      read(5,*) hinfmt  

c
c#######################################################################
c
c     Inquire whether default plotting parameters are desired
c
c#######################################################################
c
      write(6,'(/a,/a,/a,/a/)')
     :' Do you wish to let the program determine the rest of the',
     :' plotting parameters? If so, the program will plot 2-D' ,
     :' x-y, x-z and y-z slices of most fields through the center',
     :' of model domain. Type in 1 or 0 for yes or no.'
      read(5,*) ireply
c
c#######################################################################
c
c     Inquire about which plotting parameters are desired
c
c#######################################################################
c
      IF( ireply.ne.1) THEN
 
        write(6,'(/a,/a/)')
     : ' Then tell the program what variables to plot,',
     : ' by answering 1 or 0 to the following questions.'

        write(6,'(/a/)') ' First, the total quantities.'

        write(6,'(a)') ' Do you want to plot wind vector field?' 
        read(5,*) vtrplt

        write(6,'(a)') ' Do you want to plot total u-velocity?'
        read(5,*) uplot

        write(6,'(a)') ' Do you want to plot total v-velocity?'
        read(5,*) vplot

        write(6,'(a)') ' Do you want to plot total w-velocity?'
        read(5,*) wplot

        write(6,'(a)') 
     :  ' Do you want to plot total potential temperature?'
        read(5,*) ptplot

        write(6,'(a)') ' Do you want to plot total pressure?'
        read(5,*) pplot

        write(6,'(a)') 
     :  ' Do you want to plot total water vapor mixing ratio?'
        read(5,*) qvplot

        write(6,'(a)') 
     :  ' Do you want to plot total cloud water mixing ratio?'
        read(5,*) qcplot

        write(6,'(a)') 
     :  ' Do you want to plot total rain water mixing ratio?'
        read(5,*) qrplot

        write(6,'(a)') ' Do you want to plot total ice mixing ratio?'
        read(5,*) qiplot

        write(6,'(a)') ' Do you want to plot total snow mixing ratio?'
        read(5,*) qsplot

        write(6,'(a)') ' Do you want to plot total hail mixing ratio?'
        read(5,*) qhplot

        write(6,'(a)') ' Do you want to plot total hail mixing ratio?'
        read(5,*) kmplot

        write(6,'(a)') ' Do you want to plot relative humidity?'
        read(5,*) rhplot

        write(6,'(a)') 
     :' Do you want to plot equivalent potential temperature?'
        read(5,*) pteplt 

        write(6,'(/a/)') ' Then the perturbations from the base state:'

        write(6,'(a)') 
     :      ' Do you want to plot perturbation wind vector field?' 
        read(5,*) vtpplt

        write(6,'(a)') ' Do you want to plot perturbation u velocity?'
        read(5,*) upplot

        write(6,'(a)') ' Do you want to plot perturbation v velocity?'
        read(5,*) vpplot

        write(6,'(a)') ' Do you want to plot perturbation w velocity?'
        read(5,*) wpplot

        write(6,'(a)') 
     :  ' Do you want to plot perturbation potential temperature?'
        read(5,*) ptpplt

        write(6,'(a)') ' Do you want to plot perturbation pressure?'
        read(5,*) ppplot

        write(6,'(a)') ' Do you want to plot perturbation water vapor?'
        read(5,*) qvpplt

        write(6,'(/a/)') ' Now tell the program which slices you want.'

        write(6,'(a)') ' Do you want to plot x-y slices? (1 or 0)'
        read(5,*) iplot(1)

        write(6,'(a)') ' Do you want to plot x-z slices? (1 or 0)'
        read(5,*) iplot(2)

        write(6,'(a)') ' Do you want to plot y-z slices? (1 or 0)'
        read(5,*) iplot(3)

        write(6,'(/a)') 
     :' Now input the i, j or k index for the slice to be selected.'
        write(6,'(a/)') 
     :' if i<0, i=(nx-1)/1+1 is assumed. Same is for j and k.'

        write(6,'(a,a,i3,a)') 
     :' Input the index k for the x-y slice to be plotted.',
     :' (1=< k =<',nz-1,')'
        read(5,*) kslice

        write(6,'(a,a,i3,a)') 
     :' Input the index j for the x-z slice to be plotted.',
     :' (1=< j =<',ny-1,')'
        read(5,*) jslice

        write(6,'(a,a,i3,a)') 
     :' Input the index i for the y-z slice to be plotted.',
     :' (1=< i =<',nx-1,')'
        read(5,*) islice

        IF( kslice.lt.0) kslice = (nz-2)/2+1
        IF( jslice.lt.0) jslice = (ny-2)/2+1
        IF( islice.lt.0) islice = (nx-2)/2+1

        write(6,'(/a,/a/)') 
     :' Input the control parameters for 3-D wireframe',
     :' plots of iso-surfaces.'

        write(6,'(a)') 
     :' Do you want to plot wireframe isosurface of w? (1/0)'
        read(5,*) w3dplt

        write(6,'(a)') ' Input the value of this isosurface.'
        read(5,*) wisosf

        write(6,'(a)') 
     :' Do you want to plot wireframe isosurface of qc+qr? (1/0)'
        read(5,*) q3dplt

        write(6,'(a)') ' Input the value of this isosurface.'
        read(5,*) qisosf
c
c#######################################################################
c
c     Page set-up parameters
c
c#######################################################################
c
        write(6,'(/a,/a,/a/,a/)') 
     :' Now tell the program how to arrange your plots (pictures)',
     :' in each plotting frame. These pictures can be arranged in',
     :' multiple rows and multiple columns in a single frame,',
     :' and the portrait or landscape ploting mode can be selected.'

        write(6,'(a)')
     :' Input the plotting mode (1 for portrait, 2 for landscape):'
        read(5,*) layout

        write(6,'(a)')' Input the number of columns:'
        read(5,*) nxpic

        write(6,'(a)')' Input the number of rows:'     
        read(5,*) nypic
c
c#######################################################################
c
c     Default plotting parameters set up
c
c#######################################################################
c
      ELSE                            ! automatic plot selection

        w3dplt = 1
        q3dplt = 1

        vtrplt = 1
        uplot  = 1
        vplot  = 1
        wplot  = 1
        ptplot = 0
        pplot  = 0
        qvplot = 0
        qcplot = 0
        qrplot = 0
        qiplot = 0
        qsplot = 0
        qhplot = 0
        kmplot = 0
        rhplot = 0
        pteplt = 0

        vtpplt = 0
        upplot = 0
        vpplot = 0
        wpplot = 0
        ptpplt = 0
        ppplot = 1
        qvpplt = 0
c
c#######################################################################
c
c     set the plotting parameter for each mode
c
c#######################################################################
c
        iplot(1)=1       ! control for plotting x-y slice
        iplot(2)=0       ! control for plotting x-z slice
        iplot(3)=0       ! control for plotting y-z slice

        kslice = (nz-2)/2+1
        jslice = (ny-2)/2+1
        islice = (nx-2)/2+1

        nxpic = 2  ! number of pictures in horizontal in each frame
        nypic = 2  ! number of pictures in vertical in each frame

      END IF

c
c#######################################################################
c
c     Initialize plotting package and set plotting space parameters
c
c#######################################################################
c
      CALL xdevic 

      IF( nxpic.eq.1 .and. nypic.eq.1) call xdspac(0.80)

      IF( layout .eq. 1 ) THEN
        angl = 00.0
      ELSE
        angl = 90.0
      ENDIF

      CALL xspace(nxpic, nypic, angl , xlimit,ylimit)
c
c#######################################################################
c
c     Set the margins between pictures in nondimensional units.
c
c#######################################################################
c
      CALL xpmagn(0.05,0.05)

c
c#######################################################################
c
c     Zero-out data arrays
c
c#######################################################################
c
      CALL filzero(u     ,nxyz)
      CALL filzero(v     ,nxyz) 
      CALL filzero(w     ,nxyz) 
      CALL filzero(pt     ,nxyz) 
      CALL filzero(qv    ,nxyz) 
      CALL filzero(qc    ,nxyz) 
      CALL filzero(qr    ,nxyz) 
      CALL filzero(qi    ,nxyz) 
      CALL filzero(qs    ,nxyz) 
      CALL filzero(qh    ,nxyz) 
      CALL filzero(km    ,nxyz) 

      CALL filzero(uprt  ,nxyz)
      CALL filzero(vprt  ,nxyz) 
      CALL filzero(wprt  ,nxyz) 
      CALL filzero(ptprt  ,nxyz) 
      CALL filzero(pprt  ,nxyz) 
      CALL filzero(qvprt ,nxyz) 

      CALL filzero(ubar  ,nxyz) 
      CALL filzero(vbar  ,nxyz) 
      CALL filzero(wbar  ,nxyz) 
      CALL filzero(ptbar ,nxyz) 
      CALL filzero(pbar  ,nxyz) 
      CALL filzero(qvbar ,nxyz) 
c
c#######################################################################
c
c     Get the name of the input data set.
c
c#######################################################################
c

      write(6,'(/a/a/a)')
     :' Please give the name of the file',
     :' containing the grid and base state information.',
     :' This is not used if these fields are stored in each data file.'
      read(5,*) grdbasfn

      lengbf = 80
      CALL slength( filename, lengbf)
      write(6,'(a,a/)')' The grid/base name is ', grdbasfn(1:lengbf) 

10    CONTINUE

      write(6,'(/a)')' Please give the name of the input data file.'
      read(5,*) filename

      lenfil = 80
      CALL slength( filename, lenfil)
      write(6,'(a,a/)')' The data set name is ', filename(1:lenfil) 
c
c#######################################################################
c
c     Read all input data arrays 
c
c#######################################################################
c
      CALL dtaread(nx,ny,nz, hinfmt ,grdbasfn(1:lengbf),lengbf,
     :            filename(1:lenfil),lenfil,time,
     :            x,y,zc,z, uprt ,vprt ,wprt ,ptprt, pprt ,
     :            qvprt, qc, qr, qi, qs, qh, km,
     :            ubar, vbar, wbar, ptbar, pbar, rhobar, qvbar,
     :            ireturn, tem1)

      curtim = time
c
c#######################################################################
c 
c     ireturn = 0 for a successful read
c
c#######################################################################
c
      IF( ireturn.eq.0 ) THEN   ! successful read

c
c#######################################################################
c
c     Loop through three plotting modes
c     producing x-y, x-z and y-z slices of 2-d plotting each time.
c     If iplot(mode)=0, plotting is skipped. 
c
c#######################################################################
c
        DO 130 mode = 1, 3

          IF( iplot(mode).ne.0 ) THEN

            IF( mode.eq. 1) yxratio = 1.0
            IF( mode.eq. 2) yxratio = 2.0
            IF( mode.eq. 3) yxratio = 2.0

            ibgn = 2
            jbgn = 2
            kbgn = 2

            iend = nx-2
            jend = ny-2
            kend = nz-2
 
            IF ( nx.gt.90 ) THEN
               ist = 3
            ELSE IF( nx.gt.45) THEN
               ist = 2
            ELSE
               ist = 1
            END IF

            IF ( ny.gt.90 ) THEN
               jst = 3
            ELSE IF( ny.gt.45) THEN
               jst = 2
            ELSE
               jst = 1
            END IF

            IF ( nz.gt.90 ) THEN
               kst = 3
            ELSE IF( nz.gt.45) THEN
               kst = 2
            ELSE
               kst = 1
            END IF

            dx = x(3,2,2)-x(2,2,2)
            dy = y(2,3,2)-y(2,2,2)
            dz = z(2,2,3)-z(2,2,2)

            xr = (iend-ibgn+1)*dx 
            yr = (jend-jbgn+1)*dy
            zr = (kend-kbgn+1)*dz 

            x1 = (ibgn-2)*dx 
            x2 = x1 + xr
            y1 = (jbgn-2)*dy
            y2 = y1 + yr
            z1 = (kbgn-2)*dz
            z2 = z1 + zr
c
c#######################################################################
c
c      Set an appropriate yxratio
c
c#######################################################################
c
            IF(mode .eq. 1) THEN
              aspect=xr/yr
            ELSE IF (mode .eq. 2) THEN
              aspect=xr/zr
            ELSE
              aspect=yr/zr
            ENDIF
c
            IF( aspect .lt. 1.5 ) THEN
              yxratio=1.
            ELSE IF( aspect .lt. 3.) THEN
              yxratio=2.
            ELSE IF( aspect .lt. 7.) THEN
              yxratio=5.
            ELSE IF( aspect .lt. 12.) THEN
              yxratio=10.
            ELSE
              yxratio=20.
            ENDIF
c
c
c#######################################################################
c
c      Plot vector field
c
c#######################################################################
c
            DO 310 k=1,nz
            DO 310 j=1,ny
            DO 310 i=1,nx
              u(i,j,k)=uprt(i,j,k)+ubar(i,j,k)
              v(i,j,k)=vprt(i,j,k)+vbar(i,j,k)
              w(i,j,k)=wprt(i,j,k)+wbar(i,j,k)
              qv(i,j,k)=qvprt(i,j,k)+qvbar(i,j,k)
310         CONTINUE

            IF(vtrplt.eq.1) THEN
              CALL vtrunt ( 10.0 )
              CALL vtr3d( u,v,w,x,y,z,
     :          x1,x2,y1,y2,z1,z2, nx,ibgn,iend,ist, 
     :          ny,jbgn,jend,jst, nz,kbgn,kend,kst,
     :          kslice,jslice,islice, time, runname,1.0, mode, 
     :          tem1,tem2,tem3,tem4)
              CALL vtrunt (0.0)
            ENDIF
c
c#######################################################################
c
c      Plot wind components
c
c#######################################################################
c

            CALL xlbint(5)
            CALL xczero(3)

            IF(uplot.eq.1) THEN 
              label = 'u     (m/s)   '
              CALL ctrinc( 2.0 ) 
              CALL ctr3d( u, x,y,z,x1,x2,y1,y2,z1,z2, 
     :          nx,ibgn,iend, ny,jbgn,jend, nz,kbgn,kend,
     :          label(1:14),time,mode, kslice,jslice,islice,
     :          runname,1.0,tem1,tem2,tem3)
              CALL ctrinc( 0.0 ) 
            ENDIF

            IF(vplot.eq.1) THEN 
              label = 'v     (m/s)   '
              CALL ctrinc( 2.0 ) 
              CALL ctr3d( v, x,y,z,x1,x2,y1,y2,z1,z2,
     :          nx,ibgn,iend, ny,jbgn,jend, nz,kbgn,kend,
     :          label(1:14),time,mode, kslice,jslice,islice,
     :          runname,1.0,tem1,tem2,tem3)
              CALL ctrinc( 0.0 ) 
            ENDIF

            IF(wplot.eq.1) THEN 
              label = 'w     (m/s)   '
              CALL overlay (0)
              CALL xczero(0)
              CALL xhlfrq(4)

              CALL ctr3d( w,x,y,z,x1,x2,y1,y2,z1,z2,
     :           nx,ibgn,iend, ny,jbgn,jend, nz,kbgn,kend,
     :           label(1:14),time,mode, kslice,jslice,islice,
     :           runname,1.0,tem1,tem2,tem3)

              CALL overlay (0)
              CALL xczero(3)
              CALL xhlfrq(4)

            ENDIF
c
c#######################################################################
c
c      Plot scalars
c
c#######################################################################
c
            IF(ptplot.eq.1) THEN 
              DO 315 k=1,nz
              DO 315 j=1,ny
              DO 315 i=1,nx
                pt(i,j,k)=ptprt(i,j,k)+ptbar(i,j,k)
315           CONTINUE

              label = 'ptemp  (K)    '
              factor=1.
              CALL ctr3d( pt ,x,y,z,x1,x2,y1,y2,z1,z2,
     :          nx,ibgn,iend, ny,jbgn,jend, nz,kbgn,kend,
     :          label(1:14),time,mode, kslice,jslice,islice,
     :          runname,factor,tem1,tem2,tem3)
            ENDIF

            IF(pplot.eq.1) THEN 
              DO 60 k=1,nz-1
              DO 60 j=1,ny-1
              DO 60 i=1,nx-1
                tem4(i,j,k)=pbar(i,j,k)+pprt (i,j,k)
60            CONTINUE
              label = 'p (PASCAL)   '
              factor=1.
              CALL ctr3d( tem4 ,x,y,z,x1,x2,y1,y2,z1,z2, 
     :          nx,ibgn,iend, ny,jbgn,jend, nz,kbgn,kend,
     :          label(1:14),time,mode, kslice,jslice,islice,
     :          runname,factor,tem1,tem2,tem3)
            ENDIF

            IF(qvplot.eq.1) THEN 
              label = 'qv     (g/kg) '
              factor=1000.
              CALL ctr3d( qv,x,y,z,x1,x2,y1,y2,z1,z2,
     :           nx,ibgn,iend, ny,jbgn,jend, nz,kbgn,kend,
     :          label(1:14),time,mode, kslice,jslice,islice,
     :          runname,factor,tem1,tem2,tem3)
            ENDIF

            IF(qcplot.eq.1) THEN 
              label = 'qc     (g/kg) '
              factor=1000.
              CALL ctr3d( qc,x,y,z,x1,x2,y1,y2,z1,z2, 
     :          nx,ibgn,iend, ny,jbgn,jend, nz,kbgn,kend,
     :          label(1:14),time,mode, kslice,jslice,islice,
     :          runname,factor,tem1,tem2,tem3)
            ENDIF

            IF(qrplot.eq.1) THEN 
              label = 'qr     (g/kg) '
              factor=1000.
              CALL ctr3d( qr,x,y,z,x1,x2,y1,y2,z1,z2, 
     :          nx,ibgn,iend, ny,jbgn,jend, nz,kbgn,kend,
     :          label(1:14),time,mode, kslice,jslice,islice,
     :          runname,factor,tem1,tem2,tem3)
            ENDIF

            IF(qiplot.eq.1) THEN 
              label = 'qi     (g/kg) '
              factor=1000.
              CALL ctr3d( qi,x,y,z,x1,x2,y1,y2,z1,z2, 
     :          nx,ibgn,iend, ny,jbgn,jend, nz,kbgn,kend,
     :          label(1:14),time,mode, kslice,jslice,islice,
     :          runname,factor,tem1,tem2,tem3)
            ENDIF

            IF(qsplot.eq.1) THEN 
              label = 'qs     (g/kg) '
              factor=1000.
              CALL ctr3d( qs,x,y,z,x1,x2,y1,y2,z1,z2, 
     :          nx,ibgn,iend, ny,jbgn,jend, nz,kbgn,kend,
     :          label(1:14),time,mode, kslice,jslice,islice,
     :          runname,factor,tem1,tem2,tem3)
            ENDIF

            IF(qhplot.eq.1) THEN 
              label = 'qh     (g/kg) '
              factor=1000.
              CALL ctr3d( qh,x,y,z,x1,x2,y1,y2,z1,z2,
     :          nx,ibgn,iend, ny,jbgn,jend, nz,kbgn,kend,
     :          label(1:14),time,mode, kslice,jslice,islice,
     :          runname,factor,tem1,tem2,tem3)
            ENDIF

            IF(kmplot.eq.1) THEN 
              label = 'km (m**2/s)   '
              factor=1.
              CALL ctr3d( km,x,y,z,x1,x2,y1,y2,z1,z2,
     :          nx,ibgn,iend, ny,jbgn,jend, nz,kbgn,kend,
     :          label(1:14),time,mode, kslice,jslice,islice,
     :          runname,factor,tem1,tem2,tem3)
            ENDIF

            IF(rhplot.eq.1) THEN 
c
c#######################################################################
c
c     Calculate relative humidity
c        Here  tem1 = temperature, tem2 = saturation qv
c
c#######################################################################
c
              CALL temper (nx,ny,nz,pt, pprt ,pbar,tem1)
              CALL satmr (nx,ny,nz,pbar, tem1, tem2) 
              DO 102 k=1,nz-1
              DO 102 j=1,ny-1
              DO 102 i=1,nx-1
                tem4(i,j,k) = (qvbar(i,j,k)+qvprt(i,j,k))/tem2(i,j,k)
102           CONTINUE

              label = 'RH'    
              factor=10.
              CALL ctr3d( tem4,x,y,z,x1,x2,y1,y2,z1,z2, 
     :          nx,ibgn,iend, ny,jbgn,jend, nz,kbgn,kend,
     :          label(1:14),time,mode, kslice,jslice,islice,
     :          runname,factor,tem1,tem2,tem3)
            ENDIF

            IF(pteplt.eq.1) THEN 
              print*,'plotting of eq. pot. temp. not implemented yet.'
c
c#######################################################################
c
c     Calculate equivalent potential temperature. 
c
c#######################################################################
c
              Level =1
c             CALL getpte() 

              label = 'thetaE (K)    '    
              factor=10.
c
c not implemented yet.
c
c             CALL ctr3d( tem4,x,y,z,x1,x2,y1,y2,z1,z2, 
c    :          nx,ibgn,iend, ny,jbgn,jend, nz,kbgn,kend,
c    :          label(1:14),time,mode, kslice,jslice,islice,
c    :          runname,factor,tem1,tem2,tem3)
            ENDIF
c
c#######################################################################
c
c     Calculate and plot perturbation wind components
c
c#######################################################################
c
            IF(vtpplt.eq.1) THEN
              CALL vtr3d( uprt , vprt , wprt, x,y,z,
     :          x1,x2,y1,y2,z1,z2, nx,ibgn,iend,ist, 
     :          ny,jbgn,jend,jst, nz,kbgn,kend,kst,
     :          kslice,jslice,islice, time, runname,1.0, mode, 
     :          tem1,tem2,tem3,tem4)
            ENDIF

            IF(upplot.eq.1) THEN 
              label = 'uprt  (m/s)   '
              factor=1.
              CALL ctr3d( uprt ,X,y,z,x1,x2,y1,y2,z1,z2,
     :          nx,ibgn,iend, ny,jbgn,jend, nz,kbgn,kend,
     :          label(1:14),time,mode, kslice,jslice,islice,
     :          runname,factor,tem1,tem2,tem3)
            ENDIF

            IF(vpplot.eq.1) THEN 
              label = 'vprt  (m/s)   '
              factor=1.
              CALL ctr3d( vprt ,x,y,z,x1,x2,y1,y2,z1,z2,
     :           nx,ibgn,iend, ny,jbgn,jend, nz,kbgn,kend,
     :           label(1:14),time,mode, kslice,jslice,islice,
     :           runname,factor,tem1,tem2,tem3)
            ENDIF

            IF(wpplot.eq.1) THEN 
              label = 'wprt  (m/s)   '
              factor=1.
              CALL ctr3d( wprt ,X,y,z,x1,x2,y1,y2,z1,z2,
     :          nx,ibgn,iend, ny,jbgn,jend, nz,kbgn,kend,
     :          label(1:14),time,mode, kslice,jslice,islice,
     :          runname,factor,tem1,tem2,tem3)
            ENDIF
c
c#######################################################################
c
c     Plot perturbation scalars, calculating and storing
c     perturbations in tem4, where necessary
c
c#######################################################################
c
            IF(ptpplt.eq.1) THEN 
              label = 'ptprt  (K)    '
              factor=1.
              CALL ctr3d(ptprt,x,y,z,x1,x2,y1,y2,z1,z2,
     :           nx,ibgn,iend, ny,jbgn,jend, nz,kbgn,kend,
     :           label(1:14),time,mode, kslice,jslice,islice,
     :           runname,factor,tem1,tem2,tem3)
             ENDIF

            IF(ppplot.eq.1) THEN 
              label = 'pprt  (PASCAL)'
              factor=1.
              CALL ctr3d(pprt ,x,y,z,x1,x2,y1,y2,z1,z2,
     :           nx,ibgn,iend, ny,jbgn,jend, nz,kbgn,kend,
     :           label(1:14),time,mode, kslice,jslice,islice,
     :           runname,factor,tem1,tem2,tem3)
            ENDIF

            IF(qvpplt.eq.1) THEN 
              label = 'qvprt  (g/kg) '
              factor=1000.
              CALL ctr3d(qvprt,x,y,z,x1,x2,y1,y2,z1,z2,
     :           nx,ibgn,iend, ny,jbgn,jend, nz,kbgn,kend,
     :           label(1:14),time,mode, kslice,jslice,islice,
     :           runname,factor,tem1,tem2,tem3)
            ENDIF 
c
c#######################################################################
c
c     3-D wireframe plots of isosurfaces.
c
c#######################################################################
c
            IF( w3dplt.eq.1) THEN 
c
c     Plot w=-wisosf isosurface.
c 
              CALL wirfrm(w,nx,ibgn,iend,ny,jbgn,jend,nz,kbgn,kend,
     :                    wisosf, tem2)

              CALL xqpspc(x1,x2,y1,y2)
              CALL xchmag(0.02*(y2-y1))
              CALL xqmap (x1,x2,y1,y2)
              CALL xcharl(x1,y1-0.02*(y2-y1),runname)

              write(label,'(a,f4.1,a,f6.1,a)') 
     :        'w=',wisosf,' m/s isosurface at t=',time/60.0,' min'

              CALL xcharl(x1,y1-0.04*(y2-y1),label )
c
c     Plot w=-wisosf isosurface.
c 
              CALL wirfrm(w,nx,ibgn,iend,ny,jbgn,jend,nz,kbgn,kend,
     :                    -wisosf, tem2)

              CALL xqpspc(x1,x2,y1,y2)
              CALL xchmag(0.02*(y2-y1))
              CALL xqmap (x1,x2,y1,y2)
              CALL xcharl(x1,y1-0.02*(y2-y1),runname)

              write(label,'(a,f4.1,a,f6.1,a)') 
     :        'w=',-wisosf,' m/s isosurface at t=',time/60.0,' min'

              CALL xcharl(x1,y1-0.04*(y2-y1),label )
  

            ENDIF

            IF( q3dplt.eq.1) THEN 
              DO 450 k=1,nz
              DO 450 j=1,ny
              DO 450 i=1,nx
                tem1(i,j,k)=(qr(i,j,k)+qc(i,j,k))*1000.0
450           CONTINUE

              CALL wirfrm(tem1, 
     :               nx,ibgn,iend,ny,jbgn,jend,nz,kbgn,kend,qisosf,tem2)

              CALL xqpspc(x1,x2,y1,y2)
              CALL xchmag(0.02*(y2-y1))
              CALL xqmap (x1,x2,y1,y2)
              CALL xcharl(x1,y1-0.02*(y2-y1),runname)

              write(label,'(a,f4.1,a,f6.1,a)') 
     :        'qr+qc=',qisosf,' kg/kg isosurface at t=',time/60.0,' min'

              CALL xcharl(x1,y1-0.04*(y2-y1),label )
  

            ENDIF
c
            IF( nxpic.eq.1 .and. nypic.eq.1) call xdspac(0.80)

            CALL xnwfrm
          ENDIF
130   CONTINUE

        write(6,'(/a/)')
     :  ' Do you want to process another data set? (0 or 1)'
        read(5,*) ireply 

        IF( ireply .eq. 1) THEN
          write(6,'(/1x,a/)') 'Reply was 1, another data will be read.'
          GO TO 10   
        ELSE
          write(6,'(/1x,a/)') 'Reply was 0, job is stopped here.'
        ENDIF

      ELSE
          write(6,'(a,/a,i3,a)') 
     :        ' Bad return status from data file reading',
     :        ' ireturn = ', ireturn,' program stopped.'
      ENDIF
      
      CALL xgrend 

      STOP
      END
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE CTR3D                      ######
c     ######                                                      ######
c     ######                Copyright (c) 1992                    ######
c     ######  Center for the Analysis and Prediction of Storms    ######
c     ######    University of Oklahoma.  All rights reserved.     ######
c     ######                                                      ######
c     ##################################################################
c     ##################################################################
c

      SUBROUTINE ctr3d(b,x,y,z, x1,x2,y1,y2,z1,z2,
     :  nx,ibgn,iend, ny,jbgn,jend, nz,kbgn,kend, label,time,mode, 
     :  kslice, jslice, islice, runname, factor, tem1,tem2,tem3)

c
c#######################################################################
c
c     PURPOSE:
c
c       Set-up 2-d slices of a 3-d data array to contour with
c         subroutine ctr2d.
c 
c
c#######################################################################
c
c     AUTHOR: Ming Xue
c
c     MODIFICATION HISTORY:
c       6/08/92  Added full documentation (K. Brewster)
C
c#######################################################################
c
c     INPUT : 
c
c       b        3-dimensional array of data to contour
c
c       x        x coordinate of grid points in physical/comp. space (m)
c       y        y coordinate of grid points in physical/comp. space (m)
c       z        z coordinate of grid points in computational space (m)
c
c       x1       value of x for first i grid point to plot
c       x2       value of x for last i grid point to plot
c       y1       value of y for first j grid point to plot
c       y2       value of y for last j grid point to plot
c       z1       value of z for first k grid point to plot
c       z2       value of z for last k grid point to plot
c
c       nx       first dimension of b
c       ibgn     index of first i grid point to plot
c       iend     index of last  i grid point to plot
c
c       ny       second dimension of b
c       jbgn     index of first j grid point to plot
c       jend     index of last  j grid point to plot
c
c       nz       third dimension of b
c       kbgn     index of first k grid point to plot
c       kend     index of last  k grid point to plot
c
c       label    character string describing the contents of b
c
c       time     time of data in seconds
c
c       mode     slice orientation indicator
c          mode = 1, x-y slice at z index kslice is plotted. 
c          mode = 2, x-z slice at y index jslice is plotted. 
c          mode = 3, y-z slice at x index islice is plotted. 
c          mode = 0, all of the three slices above are plotted. 
c 
c       kslice   k index of plane for mode=1 x-y slice
c       jslice   j index of plane for mode=2 x-z slice
c       islice   i index of plane for mode=1 y-z slice
c
c       runname  character string decribing run
c
c       factor   scaling factor for data  
c                contours are labelled b*factor
c
c     WORK ARRAYS:
c
c       tem1     Temporary work array.
c       tem2     Temporary work array.
c       tem3     Temporary work array.
c
c     (These arrays are defined and used locally (i.e. inside this
c      subroutine), they may also be passed into routines called by
c      this one. Exiting the call to this subroutine, these temporary
c      work arrays may be used for other purposes therefore their
c      contents overwritten. Please examine the usage of work arrays
c      before you alter the code.)
c
c#######################################################################
c
c     Variable Declarations.
c
c#######################################################################
c
      implicit none 
c
      integer nx,ny,nz
c
      real b(nx,ny,nz)
      real x(nx,ny,nz)
      real y(nx,ny,nz)
      real z(nx,ny,nz)
c
      real x1,x2,y1,y2,z1,z2
      integer ibgn,iend,jbgn,jend,kbgn,kend
c
      character label*(*)
c
      real time
c
      integer mode
      integer kslice,jslice,islice
c
      character runname*(*)
c
      real factor
c
      real tem1(*)
      real tem2(*)
      real tem3(*)
c
c#######################################################################
c
c     Misc. local variables
c
c#######################################################################
c
      integer i,j,k,ij,ik,jk
      character title*80
c
C#######################################################################
C
C     Routines Called
C
C#######################################################################
c
      external ctr2d
c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c
c
c#######################################################################
c
c     mode=1   Plot x-y cross-section
c
c#######################################################################
c
      IF(mode.eq.1 .or. mode.eq.0 ) THEN

        k = kslice
        DO 10 j=jbgn,jend
        DO 10 i=ibgn,iend
          ij = i-ibgn+1 + (j-jbgn)*(iend-ibgn+1)
          tem1(ij)=b(i,j,k)*factor
          tem2(ij)=x(i,j,k)
          tem3(ij)=y(i,j,k)
10      CONTINUE
        write(title,'(''x-y slice of '',a,'' at t='',f8.1
     :    ,'' (s) at k='',i2)') label,time,k
        CALL ctr2d(tem1,tem2,tem3, x1,x2, y1,y2, 
     :  iend-ibgn+1,jend-jbgn+1,title(1:54),runname)

c
c#######################################################################
c
c     mode=2   Plot x-z cross-section
c
c#######################################################################
c
      ELSE IF (mode.eq.2 .or. mode.eq.0 ) THEN

        j = jslice
        DO 20 k=kbgn,kend
        DO 20 i=ibgn,iend
          ik = i-ibgn+1 + (k-kbgn)*(iend-ibgn+1)
          tem1(ik)=b(i,j,k)*factor
          tem2(ik)=x(i,j,k)
          tem3(ik)=z(i,j,k)
20      CONTINUE

        write(title,'(''x-z slice of '',a,'' at t='',f8.1
     :    ,'' (s) at j='',i2)') label,time,j
        CALL ctr2d(tem1,tem2,tem3, x1,x2, z1,z2, 
     :   iend-ibgn+1,kend-kbgn+1,title(1:54),runname)

c
c#######################################################################
c
c     mode=3   Plot y-z cross-section
c
c#######################################################################
c
      ELSE IF ( mode.eq.3 .or. mode.eq.0) THEN 

        i = islice
        DO 30 k=kbgn,kend
        DO 30 j=jbgn,jend
          jk = j-jbgn+1 + (k-kbgn)*(jend-jbgn+1)
          tem1(jk)=b(i,j,k)*factor
          tem2(jk)=y(i,j,k)
          tem3(jk)=z(i,j,k)
30      CONTINUE

        write(title,'(''y-z slice of '',a,'' at t='',f8.1
     :    ,'' (s) at i='',i2)') label,time,i
        CALL ctr2d(tem1,tem2,tem3, y1,y2, z1,z2, 
     :  jend-jbgn+1,kend-kbgn+1,title(1:54),runname)

      ENDIF

      RETURN
      END
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE CTR2D                      ######
c     ######                                                      ######
c     ######                Copyright (c) 1992                    ######
c     ######  Center for the Analysis and Prediction of Storms    ######
c     ######    University of Oklahoma.  All rights reserved.     ######
c     ######                                                      ######
c     ##################################################################
c     ##################################################################
c
      SUBROUTINE ctr2d(a,x,y,xl,xr,yb,yt,m,n,title,runname)
c
c#######################################################################
c
c     PURPOSE:
c
c       Generate contour plots of 2-d field A given its coordinates 
c         using ZXPLOT package.. 
c
c#######################################################################
c
c     AUTHOR: Ming Xue
c
c     MODIFICATION HISTORY:
c       6/08/92  Added full documentation (K. Brewster)
C
c#######################################################################
c
c     INPUT: 
c
c       a        2-dimensional slice of data to contour
c
c       x        x coordinate of grid points in plot space (over on page)
c       y        y coordinate of grid points in plot space (up on page)
c
c
c       xl       Left bound of the physical domain
c       xr       Right bound of the physical domain
c       yb       Bottom bound of the physical domain.  
c       yt       Top bound of the physical domain.  
c
c       m        first dimension of a
c       n        second dimension of a
c
c       title    character string describing the contents of a
c       runname  character string describing the model run
c
c#######################################################################
c
c     Variable Declarations.
c
c#######################################################################
c
      implicit none 
c
      integer m,n
c
      real a(m,n)
      real x(m,n)
      real y(m,n)
      real xl,xr,yb,yt
c
      character title*(*), runname*(*)
c
c#######################################################################
c
c     Plotting control common blocks
c
c#######################################################################
c
      integer layover
      common /laypar/ layover
c
      real ctinc,vtunt
      common /incunt/ ctinc, vtunt  ! contour interval and vector unit 
c
      real yxratio
      common /yratio/ yxratio       ! the scaling factor the y/x ratio. 
c
c#######################################################################
c
c     Misc. local variables
c
c#######################################################################
c
      character ch*50
      integer iwrk(20000)
      integer i,j
      real cl(150)       ! 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
c                                      up-down    center
      real xs,ys         ! real space left-right length and up-down height
      real axnmag        ! size of x tick marks
      real xstep,ystep   ! real space distance between tick marks, x and y axis
      real zinc          ! contour interval
      real zmin,zmax     ! max and min of data array
      integer ncl,mode,nopic
      real xlimit, ylimit, rotang
      integer nhpic, nvpic
c
C#######################################################################
C
C     ZXPLOT Library routines called
C
C#######################################################################
c
      real xfinc
      external xnwpic
      external xqpspc
      external xpspac
      external xmap
      external xaxnmg
      external xaxsca
      external xnctrs
      external xfinc
      external xconta
      external xchsiz
      external xcharc
      external xcharl
c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c
c
C#######################################################################
C
C     Check for adequate room in work array
C
C#######################################################################
c
      IF( (m*n).gt.20000 ) THEN
        write(6,'(a)') 
     :     '  Work array IW in plotting routine ctr2d too small'
        write(6,'(a,i6,a,i6)')
     :     '  Space needed:',m*n,', allocated:',20000
        write(6,'(a)') 
     :     '  No plotting performed in this routine.'
        RETURN
      ENDIF

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

      IF( layover.eq.0 ) CALL xnwpic                      
c
C#######################################################################
C
C     Get plotting space variables
C
C#######################################################################
c

      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
      
c
C#######################################################################
C
C     Let the longest lenth determine size scaling of plot
C
C#######################################################################
c
      IF( py/px .ge. (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) 
      ENDIF
c
C#######################################################################
C
C     Set the real distance to plot distance scaling
C
C#######################################################################
c
      CALL xmap( xl, xr, yb, yt )
c
C#######################################################################
C
C     Set-up variables for tick marks and draw axes
C
C#######################################################################
c
      axnmag = 0.03*min(py, px) 

      xstep = (x(m,1)-x(1,1))/(m-1)
      ystep = (y(1,n)-y(1,1))/(n-1)

      CALL xaxnmg( axnmag )

      IF( layover.eq.0) call xaxsca( xl,xr,xstep, yb,yt,ystep )
c
C#######################################################################
C
C     Find max and min of data array
C
C#######################################################################
c
      zmin=a(1,1)
      zmax=zmin

      DO 10 j=1,n
      DO 10 i=1,m
       zmax=max(zmax,a(i,j))
10     zmin=min(zmin,a(i,j))

c
C#######################################################################
C
C     Find proper contour interval and then contour field
c        using ZXPLOT routine xconta
C
C#######################################################################
c
      cl(1)=0.0
      IF( zmax-zmin.gt.1.0e-8 ) THEN
        IF( ctinc.eq.0.0) THEN 
          cl(2)=cl(1)+ xfinc(zmax-zmin)/2
          IF(cl(2)-cl(1).eq.0.0) cl(2)=cl(1)+1.0
          mode=1
          CALL xnctrs(10,30)
        ELSE
          cl(2)=cl(1)+ctinc 
          mode=1
          CALL xnctrs(1,300)
        ENDIF
        CALL xconta(a,x,y,iwrk,m,m,n,cl,ncl,mode)
      ELSE
        cl(2)=1.0
        ncl=2
      ENDIF
      zinc = cl(min(2,ncl))-cl(1)
c
C#######################################################################
C
C     Fill label strings and write on plot
C
C#######################################################################
c
      write(ch,'(''Min='',g9.3e2,'' Max='',g9.3e2,
     :    '' Inc='',g9.3e2)')zmin,zmax,zinc

      CALL xchsiz( 0.025*ys)
      CALL xcharc( xl+xs/2, yb-(0.08+layover*0.03)*ys ,ch(1:41) )
      CALL xcharc( xl+xs/2, yt+(0.03+layover*0.03)*ys, title )

      CALL xqnpic(nopic)
      CALL xqspac(nhpic, nvpic, rotang, xlimit, ylimit)
      IF(nopic.eq.nhpic*(nvpic-1)+1 ) THEN
        CALL xcharl(xl,yb-(0.15+layover*0.03)*(yt-yb), runname)
      ENDIF

      RETURN
      END
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE CTRINC                     ######
c     ######                                                      ######
c     ######                Copyright (c) 1992                    ######
c     ######  Center for the Analysis and Prediction of Storms    ######
c     ######    University of Oklahoma.  All rights reserved.     ######
c     ######                                                      ######
c     ##################################################################
c     ##################################################################
c

      SUBROUTINE ctrinc( ctinc0 )  

c
c#######################################################################
c
c     PURPOSE:
c
c       Set the contour interval for field to plotted by CTR2D.
c
c
c#######################################################################
c
c     AUTHOR: Ming Xue
c
c     MODIFICATION HISTORY:
c       6/08/92  Added full documentation (K. Brewster)
C
c#######################################################################
c
c     INPUT: 
c
c       ctinc0    Contour interval 
c                 If CTINC0 = 0.0, the interval is internally determined.
c
c#######################################################################
c
c     Variable Declarations.
c
c#######################################################################
c
      implicit none
      real ctinc0
c
c#######################################################################
c
c     Plotting control common blocks
c
c#######################################################################
c
      real ctinc,vtunt
      common /incunt/ ctinc, vtunt  ! contour interval and vector unit 
c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c

      ctinc = ctinc0

      RETURN
      END
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE VTR3D                      ######
c     ######                                                      ######
c     ######                Copyright (c) 1992                    ######
c     ######  Center for the Analysis and Prediction of Storms    ######
c     ######    University of Oklahoma.  All rights reserved.     ######
c     ######                                                      ######
c     ##################################################################
c     ##################################################################
c

      SUBROUTINE vtr3d(u,v,w, x,y,z, xw,xe, ys,yn, zb,zt, 
     :           nx,ibgn,iend,ist, ny,jbgn,jend,jst, nz,kbgn,kend,kst,
     :           kslice, jslice, islice, time, runname, factor, mode,
     :           tem1,tem2,tem3,tem4)

c
c#######################################################################
c
c     PURPOSE:
c
c       Plot vector fields in 2-d slices 
c        
c     AUTHOR: Ming Xue
c
c     MODIFICATION HISTORY:
c       6/08/92  Added full documentation (K. Brewster)
C
c#######################################################################
c
c     INPUT: 
c
c       u        3-dimensional array of u wind components (m/s)
c       v        3-dimensional array of v wind components (m/s)
c       w        3-dimensional array of w wind components (m/s)
c
c       x        x coordinate of grid points in physical/comp. space (m)
c       y        y coordinate of grid points in physical/comp. space (m)
c       z        z coordinate of grid points in computational space (m)
c
c       xw       value of x for first i grid point to plot
c       xe       value of x for last i grid point to plot
c       ys       value of y for first j grid point to plot
c       yn       value of y for last j grid point to plot
c       zb       value of z for first k grid point to plot
c       zt       value of z for last k grid point to plot
c
c       nx       first dimension of b
c       ibgn     index of first i grid point to plot
c       iend     index of last  i grid point to plot
c
c       ny       second dimension of b
c       jbgn     index of first j grid point to plot
c       jend     index of last  j grid point to plot
c
c       nz       third dimension of b
c       kbgn     index of first k grid point to plot
c       kend     index of last  k grid point to plot
c    
c       ist      step size in x direction
c       jst      step size in y direction
c       kst      step size in z direction
c
c       time     time of data in seconds
c
c       kslice   k index of plane for mode=1 x-y slice
c       jslice   j index of plane for mode=2 x-z slice
c       islice   i index of plane for mode=1 y-z slice
c
c       runname  character string decribing run
c
c       factor   scaling factor for winds 
c                V*factor wind vectors are plotted
c
c       mode     slice orientation indicator
c          mode = 1, x-y slice of u,v at z index kslice is plotted. 
c          mode = 2, x-z slice of u,w at y index jslice is plotted. 
c          mode = 3, y-z slice of v,w at x index islice is plotted. 
c          mode = 0, all of the three slices above are plotted. 
c 
c     WORK ARRAYS:
c
c       tem1     Temporary work array.
c       tem2     Temporary work array.
c       tem3     Temporary work array.
c       tem4     Temporary work array.
c
c     (These arrays are defined and used locally (i.e. inside this
c      subroutine), they may also be passed into routines called by
c      this one. Exiting the call to this subroutine, these temporary
c      work arrays may be used for other purposes therefore their
c      contents overwritten. Please examine the usage of work arrays
c      before you alter the code.)
c
c#######################################################################
c
c     Variable Declarations.
c
c#######################################################################
c
      implicit none
c
      integer nx,ny,nz
c
      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)
c
      integer kslice,jslice,islice
      character runname*(*)

      real xw,xe,ys,yn,zb,zt
      integer ibgn,iend,ist, jbgn,jend,jst, kbgn,kend,kst

      real time,factor
      integer mode
c
c#######################################################################
c
c     Work arrays: tem1,tem2,tem3,tem4 of size at least
c             max( nx*ny, nx*nz, ny*nz).
c
c#######################################################################
c
      real tem1(*),tem2(*),tem3(*),tem4(*) 
c
c#######################################################################
c
c     Misc. local Variables
c
c#######################################################################
c
      integer i,j,k,ij,ik,jk,istep,jstep
      real uunit
      character title*80
c
C#######################################################################
C
C     ZXPLOT Library routines called
C
C#######################################################################
c
      external xvmode
c
C#######################################################################
C
C     Other routines called
C
C#######################################################################
c
      external vtr2d
c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c
C#######################################################################
C
C     mode=1   Plot u-v field
C
C#######################################################################
c
      IF( mode.eq.1 .or. mode.eq.0 ) THEN

        k = kslice
        DO 10 j=jbgn,jend
        DO 10 i=ibgn,iend
          ij = i-ibgn+1 + (j-jbgn)*(iend-ibgn+1)
          tem1(ij)=u(i,j,k)*factor
          tem2(ij)=v(i,j,k)*factor
          tem3(ij)=x(i,j,k)
          tem4(ij)=y(i,j,k)
10      CONTINUE

        write(title,
     :    '(''u-v field at t='',f8.1,'' (s) at k='',i2)') time,k

        uunit = 10.0
        CALL xvmode(1)
        istep = ist
        jstep = jst
        CALL vtr2d(tem1,tem2,tem3,tem4, uunit, xw,xe,ys,yn,
     :      iend-ibgn+1,istep,jend-jbgn+1,jstep,
     :      title(1:35),runname, 1)

c
C#######################################################################
C
C     mode=2   Plot u-w field
C
C#######################################################################
c
      ELSE IF( mode.eq.2 .or. mode.eq.0 ) THEN

        j = jslice
        DO 20 k=kbgn,kend
        DO 20 i=ibgn,iend
          ik = i-ibgn+1 + (k-kbgn)*(iend-ibgn+1)
          tem1(ik)=u(i,j,k)*factor
          tem2(ik)=w(i,j,k)*factor
          tem3(ik)=x(i,j,k)
          tem4(ik)=z(i,j,k)
20      CONTINUE

        write(title,
     :    '(''u-w field at t='',f8.1,'' (s) at j='',i2)') time,j

        uunit = 10.0 
        CALL xvmode(1)
        istep = ist
        jstep = kst
        CALL vtr2d(tem1,tem2,tem3,tem4,uunit, xw,xe,zb,zt,
     :       iend-ibgn+1,istep,kend-kbgn+1,jstep,
     :       title(1:35),runname, 2)

c
C#######################################################################
C
C     mode=3   Plot v-w field
C
C#######################################################################
c
      ELSE IF( mode.eq.3 .or. mode.eq.0 ) THEN 

        i = islice
        DO 30 k=kbgn,kend
        DO 30 j=jbgn,jend
          jk = j-jbgn+1 + (k-kbgn)*(jend-jbgn+1)
          tem1(jk)=v(i,j,k)*factor
          tem2(jk)=w(i,j,k)*factor
          tem3(jk)=y(i,j,k)
          tem4(jk)=z(i,j,k)
30      CONTINUE

        write(title,'(''v-w field at t='',f8.1,'' (s) at i='',i2)') 
     :       time,i

        uunit = 10.0 
        CALL xvmode(1)
        istep = jst
        jstep = kst
        CALL vtr2d(tem1,tem2,tem3,tem4,uunit, ys,yn,zb,zt,
     :       jend-jbgn+1,istep,kend-kbgn+1,jstep,
     :       title(1:35),runname, 3)
      ENDIF
      RETURN
      END
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE VTR2D                      ######
c     ######                                                      ######
c     ######                Copyright (c) 1992                    ######
c     ######  Center for the Analysis and Prediction of Storms    ######
c     ######    University of Oklahoma.  All rights reserved.     ######
c     ######                                                      ######
c     ##################################################################
c     ##################################################################
c
      SUBROUTINE vtr2d(u,v,x,y,uunit1,
     :           xl,xr,yb,yt,m,istep,n,jstep,char1,char2, mode)
c
c#######################################################################
c
c     PURPOSE:
c
c       Plot 2-d wind (u1,u2) vector field defined on grid points (x,y) 
c       using ZXPLOT package.. 
c
c     INPUT:
c       u        m by n 2-dimensional array of u (left-to-right)
c                  wind components (m/s)
c       v        m by n 2-dimensionsl array of v (down-to-up)
c                  wind components (m/s)
c
c       x        x coordinate of grid points in physical/comp. space (m)
c       y        y coordinate of grid points in physical/comp. space (m)
c
c       uunit1
c
c       xl,xr    The left and right bound of the physical domain.  
c       yb,yt    Bottom and top bound of the physical domain.  
c
c       m        First dimension of vector component array
c       istep    Step increament for plotting in x direction
c
c       n        Second dimension of vector component array
c       jstep    Step increament for plotting in y direction
c
c       char1    First character string to plot (title)
c       char2    Second character string to plot (runname)
c
c       mode     mode = 1 for u-v vector   (u=u, v=v in model space)
c                mode = 2 for u-w vector   (u=u, v=w in model space)
c                mode = 3 for v-w vector   (u=v, v=w in model space)
c
c#######################################################################
c
c     Variable Declarations.
c
c#######################################################################
c
      implicit none
c
      integer m,n
c
      real u(m,n)
      real v(m,n)
      real x(m,n)
      real y(m,n)

      real uunit1
      real xl,xr,yb,yt
      integer istep,jstep
c
      character char1*(*), char2*(*)
c
      integer mode
c
c#######################################################################
c
c     Plotting control common blocks
c
c#######################################################################
c
      integer layover
      common /laypar/ layover

      real ctinc,vtunt              !contour interval and vector unit 
      common /incunt/ ctinc, vtunt
c
      real yxratio                  !the scaling factor the y/x ratio. 
      common /yratio/ yxratio
c
c#######################################################################
c
c     Misc. local variables
c
c#######################################################################
c
      integer i,j,key
      integer nopic,nhpic,nvpic
      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
c                                      up-down    center
      real x0,xlimit,xleng
      real y0,ylimit
      real xstep,ystep   ! real space distance between tick marks, x and y axis
      real axnmag        ! size of x tick marks
      real umax,umin     ! max and min of u component
      real vmax,vmin     ! max and min of v component
      real uunit
      real am
      real rotang
      character*60 ch
c
c#######################################################################
c
c     ZXPLOT Library routines called
C
C#######################################################################
c
      external xnwpic
      external xqpspc
      external xpspac
      external xmap
      external xaxnmg
      external xaxsca
      external xchsiz
      external xcharc
      external xvmode
      external xvectu
      external xvectr
      external xqnpic
      external xqspac
      external xcharl
c
C#######################################################################
C
C     Other routines called
C
C#######################################################################
c
      external xvectk1
c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c
      write(6,'(/1x,a,a)') ' Plotting ',char1

      IF( layover.eq.0) call xnwpic                      
c
C#######################################################################
C
C     Get plotting space variables
C
C#######################################################################
c
      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
c
C#######################################################################
C
C     Let the longest lenth determine size scaling of plot
C
C#######################################################################
c
      IF( py/px .ge. 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) 
      ENDIF
c
C#######################################################################
C
C     Set the real distance to plot distance scaling
C
C#######################################################################
c
      CALL xmap( xl, xr, yb,yt) 
c
C#######################################################################
C
C     Set-up variables for tick marks and draw axes
C
C#######################################################################
c
      axnmag = 0.03*min(py, px) 
      CALL xaxnmg( axnmag )

      xstep = (x(m,1)-x(1,1))/(m-1)
      ystep = (y(1,n)-y(1,1))/(n-1)

      IF( layover.eq.0) call xaxsca( xl,xr,xstep, yb,yt,ystep )
c
C#######################################################################
C
C     Find max and min of data array
C
C#######################################################################
c
      umax=u(1,1)
      umin=u(1,1)
      vmax=v(1,1)
      vmin=v(1,1)

      DO 200 j=1,n
      DO 200 i=1,m
        umax = max (umax, u(i,j))
        umin = min (umin, u(i,j))
        vmax = max (vmax, v(i,j))
        vmin = min (vmin, v(i,j))
200   CONTINUE
c
C#######################################################################
C
C     Fill label strings and write on plot
C
C#######################################################################
c
      IF( mode.eq.1 ) THEN 
        write(ch,'(''Umin='',f7.2,'' Umax='',f7.2,
     :           '' Vmin='',f7.2,'' Vmax='',f7.2)') umin,umax,vmin,vmax 
      ELSE IF( mode.eq.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 
      ENDIF 

      write(6,'(1x,a51)') ch(1:51)

      CALL xchsiz( 0.025*(yt-yb) )
      CALL xcharc( (xl+xr)/2, yb-(yt-yb)*(0.08+layover*0.03),ch(1:51) )
      CALL xcharc( (xl+xr)/2, yt+(yt-yb)*(0.03+layover*0.03),char1 )
c
c#######################################################################
c
c     Set vector unit and plot vectors.
c
c#######################################################################
c
      uunit=uunit1
      IF( vtunt.ne.0.0 ) THEN 
        uunit=vtunt
        CALL xvmode(2)
      ENDIF

      CALL xvectu(u,v,m,m,istep,n,jstep,xleng,uunit)

      CALL xvectr(u,v,x,y,m,m,istep,n,jstep,xleng,uunit)

      x0=xl-(xr-xl)*0.08
      y0=yb-(yt-yb)*0.1
      key=0
      am=0.5
      IF( ((m-1)/istep).gt.30 ) am=1.0
      CALL xvectk1(x0,y0,xleng*am,uunit*am, key)

      CALL xqnpic(nopic)
      CALL xqspac(nhpic, nvpic, rotang, xlimit, ylimit)
      IF(nopic.eq.nhpic*(nvpic-1)+1 ) THEN 
        CALL xcharl(xl,yb-(0.15+layover*0.03)*(yt-yb), char2 )
      ENDIF
     
      RETURN
      END
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE VTRUNT                     ######
c     ######                                                      ######
c     ######                Copyright (c) 1992                    ######
c     ######  Center for the Analysis and Prediction of Storms    ######
c     ######    University of Oklahoma.  All rights reserved.     ######
c     ######                                                      ######
c     ##################################################################
c     ##################################################################
c

      SUBROUTINE vtrunt( vtunt0 )  

c
c#######################################################################
c
c     PURPOSE:
c
c       Set the wind vector unit for wind field to plotted by VTR2D.
c
c#######################################################################
c
c     AUTHOR: Ming Xue
c
c     MODIFICATION HISTORY:
c       6/08/92  Added full documentation (K. Brewster)
C
c#######################################################################
c
c     INPUT: 
c
c       vtunt0    Contour interval 
c                 If VTUNT0 = 0.0, the unit is internally determined.
c
c#######################################################################
c
c     Variable Declarations.
c
c#######################################################################
c
      implicit none
      real vtunt0
c
c#######################################################################
c
c     Plotting control common blocks
c
c#######################################################################
c
      real ctinc,vtunt
      common /incunt/ ctinc, vtunt  ! contour interval and vector unit 
c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c
      vtunt = vtunt0

      RETURN
      END
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE OVERLAY                    ######
c     ######                                                      ######
c     ######                Copyright (c) 1992                    ######
c     ######  Center for the Analysis and Prediction of Storms    ######
c     ######    University of Oklahoma.  All rights reserved.     ######
c     ######                                                      ######
c     ##################################################################
c     ##################################################################
c

      SUBROUTINE overlay (layovr)

c
c#######################################################################
c
c     PURPOSE:
c
c       Set the layover counter parameter in the laypar common block
c
c#######################################################################
c
c     AUTHOR: Ming Xue
c
c     MODIFICATION HISTORY:
c       6/08/92  Added full documentation (K. Brewster)
C
c#######################################################################
c
c     INPUT: 
c
c       layovr   The 'overlay' parameter. 
c                If layover .ne. 0, the following 2-d contour plot will be  
c                superimposed on the previous plot. 
c                layover =1, 2, ... indicating this is the 
c                   layover'th (1st or 2nd ...) plot to be overlayed 
c                   on the previous one.
c
c#######################################################################
c
c     Variable Declarations.
c
c#######################################################################
c
      integer layovr
c
c#######################################################################
c
c     Plotting control common blocks
c
c#######################################################################
c
      integer layover
      common /laypar/ layover
c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c
      layover = layovr
      RETURN
      END
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE STYXRT                     ######
c     ######                                                      ######
c     ######                Copyright (c) 1992                    ######
c     ######  Center for the Analysis and Prediction of Storms    ######
c     ######    University of Oklahoma.  All rights reserved.     ######
c     ######                                                      ######
c     ##################################################################
c     ##################################################################
c

      SUBROUTINE styxrt( yxrt )  

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

      RETURN
      END
