subroutine sdpar(sdir,sdiro,rlat,th1,th2,f1,f2,spar)
c     This routine calculates the vertical shear 
c     direction parameter for SHIPS and LGEM.
c
c     Input: 
c       sdir - The direction of the vertical shear in deg clockwise from N,
c              For example sdir=45 deg is shear from the NE. 
c       sdiro - The optimal shear direction at latitude th1
c       rlat  - The storm latitude (deg N)
c       th1   - The lower latitude for the scaling of spar
c       th2   - The upper latitude for the scaling of spar
c       f1    - The scale factor for spar at lat=th1
c       f2    - The scale factor for spar at lat=th2
c
c     Output:
c       spar - The shear direction parameter
c
c     Calcuate the difference between the shear direction and the optimal direction
c     and shift it so that it lies between -90 and +90 deg. 
      adif = abs(sdir-sdiro)
      if (adif .gt. 180.0) adif = abs(360.0-adif)
      adif = adif - 90.0
c
c     Scale the difference by latitude using a cubic polynomial
      if (rlat .ge. 0.0) then
         if (rlat .le. th1) then
            fscale = f1
         elseif (rlat .ge. th2) then
            fscale = f2
         else
            z = (rlat-th1)/(th2-th1)
            fscale = f1 + 3.0*(f2-f1)*z*z - 2.0*(f2-f1)*z*z*z
         endif
      else
         if (rlat .ge. th1) then
            fscale = f1
         elseif (rlat .le. th2) then
            fscale = f2
         else
            z = (rlat-th1)/(th2-th1)
            fscale = f1 + 3.0*(f2-f1)*z*z - 2.0*(f2-f1)*z*z*z
         endif
      endif
c
      spar = fscale*adif
c
      return
      end