c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE FILZERO                    ######
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 filzero( a, n)

c
c#######################################################################
c
c     PURPOSE:
c       Fill vector a with zeros.
c
c#######################################################################
c
c     AUTHOR: St. Paul, Dead Sea Scrolls
c
c     MODIFICATIONS:
c       6/09/92  Added full documentation (K. Brewster)
c
c#######################################################################
c
      implicit none
c
      integer n
      real a(n)
c
c#######################################################################
c
c     Misc. local variables
c
c#######################################################################
c
      integer i
c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c
      DO 100 i=1,n
        a(i)=1.0E29
100   CONTINUE

      RETURN
      END
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE TEMPER                     ######
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 TEMPER ( nx,ny,nz,theta, ppert, pbar, T )

c
c#######################################################################
c
c     PURPOSE:
c
c       Using a version of Poisson's formula, calculate temperature.
c
c#######################################################################
c
c     AUTHOR: Joe Bradley
c       12/05/91
c
c     MODIFICATIONS:
c       Modified by Ming Xue so that arrays are only defined at
c                one time level.
c       6/09/92  Added full documentation and phycst include file for
c                rddcp=Rd/Cp  (K. Brewster)
c
c#######################################################################
c
c     INPUT:
c       nx       Number of grid points in the x-direction (east/west)
c       ny       Number of grid points in the y-direction (north/south)
c       nz       Number of grid points in the vertical
c
c	theta    Potential temperature (degrees Kelvin)
c	ppert    Perturbation pressure (Pascals)
c	pbar     Base state pressure (Pascals)
c
c     OUTPUT:
c
c	t        Temperature (degrees Kelvin)
c
c#######################################################################
c
c     Variable Declarations.
c
c#######################################################################
c
      implicit none
c
      integer nx,ny,nz
c
      real theta(nx,ny,nz)	! potential temperature (degrees Kelvin)
      real ppert(nx,ny,nz)	! perturbation pressure (Pascals)
      real pbar (nx,ny,nz)	! base state pressure (Pascals)
c
      real T    (nx,ny,nz)	! temperature (degrees Kelvin)
c
c#######################################################################
c
c     Include file
c
c#######################################################################
c
      include 'phycst.inc'
c
c#######################################################################
c
c     Misc. local variables
c
c#######################################################################
c
      integer i,j,k
c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c
c
c#######################################################################
c
c     Calculate the temperature using Poisson's formula.
c
c#######################################################################
c
      DO 10 k=1,nz-1
      DO 10 j=1,ny-1
      DO 10 i=1,nx-1

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

10    CONTINUE

      RETURN
      END
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE SATMR                      ######
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 SATMR (nx,ny,nz,pbar, T, qvS)

c
c#######################################################################
c
c     PURPOSE:
c
c     Calculate the saturation mixing ratio (saturation specific humidity)
c     qvs using Teten's formula.
c
c
c#######################################################################
c
c     AUTHOR: Joe Bradley
c       12/05/91
c
c     MODIFICATIONS:
c       Modified by Ming Xue so that arrays are only defined at
c                one time level.
c       6/09/92  Added full documentation and phycst include file for
c                rddcp=Rd/Cp  (K. Brewster)
c
c#######################################################################
c
c     INPUT:
c
c       nx       Number of grid points in the x-direction (east/west)
c       ny       Number of grid points in the y-direction (north/south)
c       nz       Number of grid points in the vertical
c
c	pbar	 Base state pressure (Pascals)
c	T	 Temperature (degrees Kelvin)
c
c     OUTPUT:
c
c	qvs	saturation mixing ratio (saturation specific humidity)
c		(kilograms/kilogram)
c
c#######################################################################
c
c     Variable Declarations.
c
c#######################################################################
c
      implicit none
c
      integer nx,ny,nz
c
      real pbar (nx,ny,nz)	! base state pressure (Pascals)
      real T    (nx,ny,nz)	! temperature (degrees Kelvin)
      real qvs  (nx,ny,nz)	! saturation specific humidity
				! (saturation mixing ratio) kg/kg
c
c#######################################################################
c
c     Misc. local variables
c
c#######################################################################
c
      integer i,j,k
c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

c
c#######################################################################
c
c     Teten's formula:
c 
c       qvs = (380./pbar) # exp(17.27 # (T-273.15)/(T-35.86))
c 
c                         T is temperature in degrees Kelvin.
c
c#######################################################################
c
      DO 30 k=1,nz-1
      DO 30 j=1,ny-1
      DO 30 i=1,nx-1

	qvs(i,j,k) =
     :		  (380. / pbar(I,J,K)) *
     :		  exp( 17.27 * (t(I,J,K) - 273.15) /
     :	               (t(I,J,K) - 35.86) )

30    CONTINUE

      RETURN
      END
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE SLENGTH                    ######
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 slength ( string, length )

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

c
c#######################################################################
c
c     Variable Declarations.
c
c#######################################################################
c
      implicit none
      integer length
      character string*(*)
c
c#######################################################################
c
c     Misc. local variables
c
c#######################################################################
c
      integer i
c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c
      DO 100 i = length,1,-1
        IF(string(i:i) .ne. ' ') GO TO 101
100   CONTINUE
101   CONTINUE

      length = i

      RETURN
      END
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                 FUNCTION XFINC                       ######
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

      REAL FUNCTION xfinc(x)

c
c#######################################################################
c
c     PURPOSE:
c
c       Automatically divide domain (0,x) to a number of subdomain
c       with interval xfinc which is >=4 and =<16 for fold=1.0
c
c#######################################################################
c
c     AUTHOR: Ming Xue
c       sometime
c
c     MODIFICATIONS:
c       6/09/92  Added full documentation (K. Brewster)
c
c#######################################################################
c
c     INPUT:
c        x       not sure
c
c#######################################################################
c
c     Variable Declarations.
c
c#######################################################################
c
      implicit none
c
      real x
c
c#######################################################################
c
c     Misc. local variables
c
c#######################################################################
c
      integer ipower
      real d,fold
c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c
      ipower= int( alog10(x) )
      d= int(x/(10.0**ipower))
      fold=1.0
      xfinc=0.1*x
      IF( d.ge.0.0 .and. d.lt.3.0 ) THEN
        xfinc=2.0*10.0**(ipower-1)
      ELSE IF( d.ge.3.0 .and. d.lt.7.0 ) THEN
        xfinc=5.0*10.0**(ipower-1)*fold
      ELSE IF( d.ge.7.0 .and. d.lt.10. ) THEN
        xfinc=1.0*10.0** ipower*fold
      ENDIF
      IF(xfinc .eq. 0.0) xfinc=x*0.1
      RETURN
      END  
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE XVECTK1                    ######
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 xvectk1(x0,y0, xleng, uunit, key)

c
c#######################################################################
c
c     PURPOSE:
c
c     Plot unit vectors  starting at (X0,Y0)
C     KEY=-1, 0, 1, 2, for none,in both X and Y-direction,X only, Y only
c
c#######################################################################
c
c     AUTHOR: Ming Xue
c
c     MODIFICATION HISTORY:
c
c     6/09/92 (K. Brewster)
c     Added full documentation.
c 
c#######################################################################
c
c     INPUT:
c
c       x0
c       y0
c       xleng
c       uunit
c       key
c    
c
c#######################################################################
c

c
c#######################################################################
c
c     Variable Declarations.
c
c#######################################################################
c
      implicit none
c
      real x0,y0
      real xleng,uunit
      integer key
c
c#######################################################################
c
c     Misc. local Variables
c
c#######################################################################
c
      real pi,angle1,angle2
      parameter(pi=3.14159,
     :      angle1=( 15./180+1)*pi,
     :      angle2=(-15./180+1)*pi)
      real sina1,cosa1,sina2,cosa2
      parameter(sina1=-.17365,
     :          cosa1=-.98481,
     :          sina2=-sina1,
     :          cosa2=cosa1)
      real xrg,xl,xr,xscale,xf
      real yrg,yb,yt,yscale,yf
      real dx,dx0,dx1,dx2
      real dy,dy0,dy1,dy2
      real px0,px1,px2
      real py0,py1,py2
      real dph,dpv
      real costa,sinta
      real vunit
      real ctrsiz
      real arrow
      real asym
      real xang
      real yang
      integer lch
      character ch*20
c
C#######################################################################
C
C     ZXPLOT Library routines called
C
C#######################################################################
c
      external xqrang
      external xqmap
      external xchmag
      external xtrans
      external xtpnup
      external xtpndn
      external xchlj
      external xqobag
      external xcharl
      external xchori
c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c
      CALL xqrang(xrg,yrg)
      CALL xqmap(xl,xr,yb,yt)
      xscale=xr-xl
      yscale=yt-yb
      yf=0.4+( min(xrg,yrg)    -0.4)*0.5
      ctrsiz=yf*0.025
      CALL xchmag(ctrsiz)

      vunit=uunit
      dx=xleng
      dy=xleng
      px0=x0
      py0=y0
      px1=x0+dx
      py1=y0
      px2=x0
      py2=y0+dy
      CALL xtrans(px0,py0)
      CALL xtrans(px1,py1)
      CALL xtrans(px2,py2)
      dph=sqrt( (px1-px0)**2+(py1-py0)**2  )
      dpv=sqrt( (px2-px0)**2+(py2-py0)**2  )
 5    IF( dpv.gt.(1.5*dph) ) THEN
        dpv=dpv*0.5
        dy=dy*0.5
        vunit=vunit*0.5
        GO TO 5
      ENDIF
 6    IF( dpv.lt.(0.75*dph) ) THEN
        dpv=dpv*2
        dy=dy*2
        vunit=vunit*2
        GO TO 6
      ENDIF
      px2=x0
      py2=y0+dy
      CALL xtrans(px2,py2)

      IF( key.eq.0.or.key.eq.1) THEN
        arrow=0.30*dph
        costa=(px1-px0)/dph
        sinta=(py1-py0)/dph
        dx1=arrow*(costa*cosa1-sinta*sina1)
        dy1=arrow*(sinta*cosa1+costa*sina1)
        dx2=arrow*(costa*cosa2-sinta*sina2)
        dy2=arrow*(sinta*cosa2+costa*sina2)
        CALL xtpnup(px0,py0)
        CALL xtpndn(px1    , py1)
        CALL xtpndn(px1+dx1, py1+dy1)
        CALL xtpnup(px1    , py1    )
        CALL xtpndn(px1+dx2, py1+dy2)
        write(ch,'(f6.1)') uunit
        lch=6
        CALL xchlj (ch,lch)
        CALL xcharl(x0+dx   +0.01*xscale,y0,ch(1:lch) )
      ENDIF
c
      IF( key.eq.0.or.key.eq.2 ) THEN
        arrow=0.20*dph
        costa=(px2-px0)/dpv
        sinta=(py2-py0)/dpv
        dx1=arrow*(costa*cosa1-sinta*sina1)
        dy1=arrow*(sinta*cosa1+costa*sina1)
        dx2=arrow*(costa*cosa2-sinta*sina2)
        dy2=arrow*(sinta*cosa2+costa*sina2)
        CALL xtpnup(px0,py0)
        CALL xtpndn(px2    , py2)
        CALL xtpndn(px2+dx1, py2+dy1)
        CALL xtpnup(px2    , py2    )
        CALL xtpndn(px2+dx2, py2+dy2)
        CALL xqobag( xang, yang )
        CALL xqchor( asym )
        CALL xchori(90.0+ yang- xang)
        write(ch,'(f6.1)') vunit
        lch=6
        CALL xchlj (ch,lch)
        CALL xcharl(x0-.01*xscale ,y0 ,ch(1:lch) )
        CALL xchori( asym )
      ENDIF
      RETURN
      END
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE KNTARY                     ######
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 kntary(fmtver,
     :            grdout,basout,varout,mstout,iceout,trbout,
     :            knt,kntscal)
c
c#######################################################################
c
c     PURPOSE:
c
c     Count the number of grids that should be in a grid file
c     depending on the switches that have been set and the version.
c     For now, we have one version so fmtwer is not consulted.
c
c#######################################################################
c
c     AUTHOR: Keith Brewster
c     7/14/92
c
c     MODIFICATION HISTORY:
c
c#######################################################################
c
c     INPUT:
c
c     fmtver  Format/Version number a character string describing the
c             format and version number of the grid writing routine
c
c     grdout =0 or 1. If grdout=0, grid variables are not printed. 
c     basout =0 or 1. If basout=0, base state variables are not dumped.
c     varout =0 or 1. If varout=0, model perturbation variables are not printed. 
c     mstout =0 or 1. If mstout=0, water variables are not printed. 
c     iceout =0 or 1. If iceout=0, printing of qi, qs and qh is skipped.
c     trbout =0 or 1. If trbout=0, printing of turbulence param km is skipped.
c
c     OUTPUT:
c
c     knt     Total number of 3-D grids in gridded data file
c     kntscal Total number of scalar 3-D grids in gridded data file
c
c#######################################################################
c
c     Variable Declarations.
c
c#######################################################################
c
      implicit none
c
      character*40 fmtver
      integer grdout,basout,varout,mstout,iceout,trbout
      integer knt,kntscal
c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c
      knt=0
      kntscal=0
      IF( grdout.eq.1 ) knt=knt+4
      IF( basout.eq.1 ) THEN
        knt=knt+7
        kntscal=kntscal+4
      END IF
      IF( varout.eq.1 ) THEN
        knt=knt+5
        kntscal=kntscal+2
      END IF
      IF( mstout.eq.1 ) THEN
        knt=knt+3
        kntscal=kntscal+3
        IF( iceout.eq.1 ) THEN
          knt=knt+3
          kntscal=kntscal+4
        END IF
      END IF
      IF( trbout.eq.1 ) THEN
        knt=knt+1
        kntscal=kntscal+1
      END IF

      RETURN
      END
c
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE A3DMAX0                    ######
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 A3DMAX0(a,m1,m2,i1,i2,n1,n2,j1,j2,l1,l2,k1,k2,
     :           amax,amin)
c
c#######################################################################
c
c     PURPOSE:
c
C     Find the maximum and minimum of array A in a specified subdomain.
c
c#######################################################################
c
c     AUTHOR: Ming Xue
c     10/01/91.   
c
c     MODIFICATION HISTORY:
c
c     5/06/92 (M. Xue)
c     Added full documentation.
C
c#######################################################################
c
c     INPUT :
c
c       a        A 3-D array whose max. and min. values is sought here.
c
c       m1,m2    The dimension bounds of the first  index of array a.
c       i1,i2    The bounds of 1st index defining a subdomain of array
c                a inside which the max. and min. of the array is sought
c
c       n1,n2    The dimension bounds of the second index of array a.
c       j1,j2    The bounds of 2nd index defining a subdomain of array
c                a inside which the max. and min. of the array is sought
c
c       l1,l2    The dimension bounds of the third  index of array a.
c       k1,k2    The bounds of 3rd index defining a subdomain of array
c                a inside which the max. and min. of the array is sought
c
c
c     OUTPUT:
c
c       amax     The maximum value of array a in a specified subdomain.
c       amin     The minimum value of array a in a specified subdomain.
c
c#######################################################################
c

c
c#######################################################################
c
c     Variable Declarations.
c
c#######################################################################
c
      implicit none 

      integer m1,n1,l1,m2,n2,l2
      real    a(m1:m2,n1:n2,l1:l2)

      real amin, amax
c
c#######################################################################
c
c     Misc. local variables: 
c
c#######################################################################
c
      integer i1,i2,j1,j2,k1,k2,i,j,k
c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c


      amax=a(i1,j1,k1)

      DO 10 k=k1,k2
      DO 10 j=j1,j2
      DO 10 i=i1,i2
        amax = max(amax, a(i,j,k))
 10   CONTINUE

      amin=a(i1,j1,k1)

      DO 20 k=k1,k2
      DO 20 j=j1,j2
      DO 20 i=i1,i2
        amin = min(amin, a(i,j,k))
 20   CONTINUE

      RETURN
      END
