!   ***************************************************************
!   * Storm Relative Helicity (SRH) is a measure of the           *
!   * streamwise vorticity within the inflow environment of a     *
!   * convective storm. It is calculated by multiplying the       *
!   * storm-relative inflow velocity vector (Vh-C) by the         *
!   * streamwise vorticity (Zh) and integrating this quantity     *
!   * over the inflow depth (lowest 1-3 km layers above ground    *
!   * level). It describes the extent to which corkscrew-like     *
!   * motion occurs (similar to the spiraling motion of an        *
!   * American football). SRH corresponds to the transfer of      *
!   * vorticity from the environment to an air parcel in          *
!   * convective motion and is used to predict the potential      *
!   * for tornadic development (cyclonic updraft rotation) in     *
!   * right-moving supercells.                                    *
!   *                                                             *
!   * There is no clear threshold value for SRH when forecasting  *
!   * supercells, since the formation of supercells appears to be *
!   * related more strongly to the deeper layer vertical shear.   *
!   * Larger values of 0-3-km SRH (greater than 250 m**2/s**2)    *
!   * and 0-1-km SRH (greater than 100 m**2/s**2), suggest an     *
!   * increased threat of tornadoes with supercells. For SRH,     *
!   * larger values are generally better, but there are no clear  *
!   * "boundaries" between non-tornadic and significant tornadic  *
!   * supercells.                                                 *
!   *                                                             *
!   * SRH < 100 (lowest 1 km): cutoff value                       *
!   * SRH = 150-299: supercells possible with weak tornadoes      *
!   * SRH = 300-499: very favorable to supercell development and  *
!   *                strong tornadoes                             *
!   * SRH > 450    : violent tornadoes                            *
!   ***************************************************************

! NCLFORTSTART
SUBROUTINE DCALRELHL(u, v, ght, ter, top, sreh, miy, mjx, mkzh)
    USE wrf_constants, ONLY : PI, RAD_PER_DEG, DEG_PER_RAD

    IMPLICIT NONE

    !f2py threadsafe
    !f2py intent(in,out) :: sreh

    INTEGER, INTENT(IN) :: miy, mjx, mkzh
    REAL(KIND=8), DIMENSION(miy,mjx,mkzh), INTENT(IN) :: u, v, ght
    REAL(KIND=8), INTENT(IN) :: top
    REAL(KIND=8), DIMENSION(miy,mjx), INTENT(IN) :: ter
    REAL(KIND=8), DIMENSION(miy,mjx), INTENT(OUT) :: sreh

! NCLEND

    ! This helicity code was provided by Dr. Craig Mattocks, and
    ! verified by Cindy Bruyere to produce results equivalent to
    ! those generated by RIP4. (The code came from RIP4?)

    REAL(KIND=8) :: dh, sdh, su, sv, ua, va, asp, adr, bsp, bdr
    REAL(KIND=8) :: cu, cv, x, sum
    INTEGER :: i, j, k, k10, k3, ktop
    !REAL(KIND=8), PARAMETER :: DTR=PI/180.d0, DPR=180.d0/PI

    !$OMP PARALLEL DO COLLAPSE(2) PRIVATE(i,j,k,k10,k3,ktop, cu, cv, x, &
    !$OMP sum, dh, sdh, su, sv, ua, va, asp, adr, bsp, bdr)
    DO j=1, mjx
        DO i=1, miy
            sdh = 0.D0
            su = 0.D0
            sv = 0.D0
            k3 = 0
            k10 = 0
            ktop = 0
            DO k = mkzh, 2, -1
                IF (((ght(i,j,k) - ter(i,j)) .GT. 10000.D0) .AND. (k10 .EQ. 0)) THEN
                    k10 = k
                    EXIT
                ENDIF
                IF (((ght(i,j,k) - ter(i,j)) .GT. top) .AND. (ktop .EQ. 0)) THEN
                    ktop = k
                ENDIF
                IF (((ght(i,j,k) - ter(i,j)) .GT. 3000.D0) .AND. (k3 .EQ. 0)) THEN
                    k3 = k
                ENDIF
            END DO

            IF (k10 .EQ. 0) THEN
                k10 = 2
            ENDIF

            DO k = k3, k10, -1
                dh = ght(i,j,k-1) - ght(i,j,k)
                sdh = sdh + dh
                su = su + 0.5D0*dh*(u(i,j,k-1) + u(i,j,k))
                sv = sv + 0.5D0*dh*(v(i,j,k-1) + v(i,j,k))
            END DO

            ua = su/sdh
            va = sv/sdh
            asp = SQRT(ua*ua + va*va)
            IF (ua .EQ. 0.D0 .AND. va .EQ. 0.D0) THEN
                adr = 0.D0
            ELSE
                adr = DEG_PER_RAD * (PI + ATAN2(ua,va))
            ENDIF

            bsp = 0.75D0*asp
            bdr = adr + 30.D0
            IF (bdr .GT. 360.D0) THEN
                bdr = bdr - 360.D0
            ENDIF

            cu = -bsp*SIN(bdr * RAD_PER_DEG)
            cv = -bsp*COS(bdr * RAD_PER_DEG)
            sum = 0.D0
            DO k = mkzh-1, ktop, -1
                x = ((u(i,j,k) - cu)*(v(i,j,k) - v(i,j,k+1))) - &
                    ((v(i,j,k) - cv)*(u(i,j,k) - u(i,j,k+1)))
                sum = sum + x
            END DO
            sreh(i,j) = -sum
        END DO
    END DO
    !$OMP END PARALLEL DO

    RETURN

END SUBROUTINE DCALRELHL