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