      SUBROUTINE HELIC(PW,UW,VW,HW,NLVLS,AVSPD,AVDIR,ELEV,
     1                 HELICITY,STSPD,STDIR)
C
C Based on code developed by Erik Rasmussen (1990-1991)
C Modified by David Blanchard (1992).  Computes from 0-4 km AGL.
C
C This function computes the integrated "streamwise component" of horizontal
C vorticity flux.  This parameter is commonly called helicity and is
C discussed in Davies-Jones, JAS 41, p. 2991- (1984).  Several research
C efforts have shown that this integral, from 0-3km (a "typical" inflow
C layer) is highly correlated with the development of rotation in storm
C updrafts, and subsequent tornadoes.  For additional details and later
C refs, refer to Davies-Jones, Burgess, and Foster p 588 of the Kannanasskis
C Park AMS Conf on Sev Loc Storms preprints (1990).
C
C Since actual sounding data are used, rather than data interpolated to
C constant height intervals, improving both speed and accuracy, the equation
C for helicity is:
C
C  HELICITY = INTGRL((u-c_u) * dv/dz - (v-c_v) * du/dz)Dz
C
C Note the geometric interpretation of helicity: minus twice the area swept
C out by the hodograph with vertices at the storm motion.
C
C INPUT:
C PW pressure levels
C UW u-winds
C VW v-winds
C HW heights of winds
C NLVLS number of wind levels
C AVSPD mean wind speed in cloud layer
C AVDIR mean wind direction in cloud layer
C ELEV station elevation in meters
C
C OUTPUT:
C HELICITY storm relative helicity
C STSPD storm motion speed
C STDIR storm motion direction
C
      REAL PW(1),UW(1), VW(1), HW(1), UW1(200), VW1(200)
      REAL INTERP
C
C Initialize values
C
      HELICITY = 0.
      U4 = 0.
      V4 = 0.
C
C Compute storm speed.  Use .75 the value of the 3-10 km mean wind
C and 30 deg to the right of the mean direction (i.e., 30R75).
C If the wind speed or direction are missing, exit.
C
      IF(AVDIR .EQ. 99999. .OR. AVSPD .EQ. 99999.) THEN
        HELICITY  = 99999.
        STDIR =   999.
        STSPD =    99.
        RETURN
      ELSE
        STSPD = .75*AVSPD
        STDIR = AVDIR+30.
        IF(STDIR .GT. 360) STDIR=STDIR-360.
        CALL UVCOMP(STDIR,STSPD,ST_U,ST_V,1)
      END IF
C
C Get the winds at 4 km above the ground
C
      ELEVP4=ELEV+4000.
      DO 200 NN=1,NLVLS
        IF(HW(NN) .GE. ELEVP4) THEN
          NUMLEVS=NN-1
          PT=PW(NN)
          PB=PW(NN-1)
          PLOG1=ALOG(PW(NN))
          PLOG3=ALOG(PW(NN-1))
C
C Iterate to find the level exactly.  Keep cutting level in half
C until the level is found.
C
          DO 100 ICOUNT=1,100
            P0=0.5*(PT+PB)
            PLOG2=ALOG(P0)
            HGHT=INTERP(HW(NN),HW(NN-1),PLOG1,PLOG2,PLOG3)
            IF(ABS(HGHT-ELEVP4) .LE. 10.) THEN
              U4=INTERP(UW(NN),UW(NN-1),PLOG1,PLOG2,PLOG3)
              V4=INTERP(VW(NN),VW(NN-1),PLOG1,PLOG2,PLOG3)
              GO TO 225
            ENDIF
          IF(HGHT .GT. ELEVP4) PT=P0
          IF(HGHT .LT. ELEVP4) PB=P0
100       CONTINUE
        END IF
200   CONTINUE
C
C  If we get to this section of code, we did not find the winds at
C  4 km AGL and must assume that the sounding did not go high enough.
C  Exit and set missing value flags.
C
        IF(U4 .EQ. 0. .AND. V4 .EQ. 0.) THEN
          HELICITY  = 99999.
          STDIR =   999.
          STSPD =    99.
          RETURN
        END IF
C
C We now have winds from the surface to 4 km AGL.  Load the array
C of u- and v-components into a local array.
C
225   DO 250 I=1,NUMLEVS-1
        UW1(I) = UW(I)
        VW1(I) = VW(I)
250   CONTINUE
      UW1(NUMLEVS) = U4
      VW1(NUMLEVS) = V4
C
C Compute the helicity
C
      DO  300 I = 2, NUMLEVS
        U_SHR = UW1(I) - UW1(I-1)
        V_SHR = VW1(I) - VW1(I-1)
        HELICITY = HELICITY + (UW1(I)-ST_U)*V_SHR - (VW1(I)-ST_V)*U_SHR
 300  CONTINUE
 400  CONTINUE
C
      HELICITY=-1.*HELICITY
      RETURN
      END
      SUBROUTINE UVCOMP(DIR,SPD,U,V,NLVLS)
C**++      IMPLICIT NONE
C
C Statement of purpose.
C ---------------------
C This subroutine computes rectangular wind components given wind direction
C and speed.
C
C History.
C --------                    
C Don Baker      01 Jul 84    Original version.
C Don Baker      15 May 85    Modified for CWP.
C
C Description of input and output.
C --------------------------------
C On input:
C ---------                
C DIR         Real Array    Wind directions (deg clockwise from north).
C SPD         Real Array    Wind speeds (m/s).
C NLVLS       Integer       Number of winds passed.
C
C On output:
C ----------               
C U           Real Array    U (east-west) wind component (m/s).
C V           Real Array    V (north-south) wind component (m/s).
C
C User notes:
C -----------
C 1) This routine can be used to return components for just one wind
C    report if NLVLS is passed as 1.
C
C
C Input arguments.
C
      REAL DIR(1),SPD(1)
      INTEGER NLVLS
C
C Output arguments.
C
      REAL U(1),V(1)
C
C Internal variables.
C
      REAL ANGLE
      INTEGER I
C
C Subroutine constants.
C
      PARAMETER (FLAG=99999.)
      PARAMETER (RPD=0.0174533)
C
C Calculate components.  Assign flag if direction or speed is bad.
C
      DO 100 I=1,NLVLS
         IF (DIR(I).LT.0. .OR. DIR(I).GT.360. .OR.
     +       SPD(I).LT.0. .OR. SPD(I).GT.250.) THEN
            U(I)=FLAG
            V(I)=FLAG
         ELSE
            ANGLE=RPD*DIR(I)
            U(I)=-SPD(I)*SIN(ANGLE)
            V(I)=-SPD(I)*COS(ANGLE)
         ENDIF
 100  CONTINUE
C
C EXIT.
C
      RETURN
      END
      real function interp(y1,y3,x1,x2,x3)
C**++      implicit none
C##############################################################################
C Statement of purpose.
C ---------------------
C Given Y1 AT X1, Y3 AT X3, and X2, calculate Y2 (INTERP) AT X2 via linear
C interpolation.
C
C History.
C --------                    
C D. Baker    01 Jul 84     Original version.
C
C Description of input and output.
C --------------------------------
C On input:
C ---------                
C y1          Real          Value of variable being interpolated at point x1.
C y3          Real          Value of variable being interpolated at point x3.
C x1          Real          Point corresponding to y1.
C x2          Real          Point to which interpolation is desired.
C x3          Real          Point corresponding to y3.
C
C On output:
C ----------               
C interp (y2) Real          Value of interpolated y variable at point x2.
C##############################################################################
C Input arguments
      real y1,y3,x1,x2,x3
C Perform simple linear interpolation.
      interp=y1+((y3-y1)*((x2-x1)/(x3-x1)))   !interp=y2
C Return to calling routine.
      return
      end
