!
!##################################################################
!##################################################################
!######                                                      ######
!######                SUBROUTINE ARPS_BE                    ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE arps_be(nx,ny,nz,                                            & 3,10
           pres,hgt,tk,wmr,                                             &
           lcl,lfc,el,twdf,li,cape,mcape,cin,tcap,                      &
           p1d,ht1d,t1d,tv1d,td1d,wmr1d,                                &
           partem,buoy,wload,mbuoy,pbesnd,mbesnd,tem1)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Calculate the lifting condensation level (lcl), level of free
!  convection (lfc), equilibrium level (el), max wet-bulb potential
!  temperature difference (twdf), lifted index (LI), Convective
!  Available Potential Energy (CAPE), Moist Convective Potential
!  Energy (MCAPE, includes water loading), convective inhibition
!  (CIN) and lid strength (tcap)  over the ARPS grid.
!
!-----------------------------------------------------------------------
!
!
!  AUTHOR: Keith Brewster
!  April, 1995   Based on OLAPS, hence LAPS, version of same.
!
!  MODIFICATION HISTORY:
!  3/11/1996  (Keith Brewster)
!  Cleaned-up, removed OLAPS artifacts.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
  INTEGER :: nx,ny,nz
!
!-----------------------------------------------------------------------
!
!  3-D input variables
!
!-----------------------------------------------------------------------
!
  REAL :: pres(nx,ny,nz)
  REAL :: hgt(nx,ny,nz)
  REAL :: tk(nx,ny,nz)
  REAL :: wmr(nx,ny,nz)
!
!-----------------------------------------------------------------------
!
!  Output variables (2d)
!
!-----------------------------------------------------------------------
!
  REAL :: lcl(nx,ny)
  REAL :: lfc(nx,ny)
  REAL :: el(nx,ny)
  REAL :: twdf(nx,ny)
  REAL :: li(nx,ny)
  REAL :: cape(nx,ny)
  REAL :: mcape(nx,ny)
  REAL :: cin(nx,ny)
  REAL :: tcap(nx,ny)
!
!-----------------------------------------------------------------------
!
!  Scratch space for calculations
!
!-----------------------------------------------------------------------
!
  REAL :: p1d(nz),ht1d(nz)
  REAL :: t1d(nz),tv1d(nz),td1d(nz),wmr1d(nz)
  REAL :: partem(nz),buoy(nz),wload(nz)
  REAL :: mbuoy(nz),pbesnd(nz),mbesnd(nz)
  REAL :: tem1(nx,ny)
!
!-----------------------------------------------------------------------
!
!  Functions
!
!-----------------------------------------------------------------------
!
  REAL :: wmr2td,tctotv
!
!-----------------------------------------------------------------------
!
!  Misc internal variables
!
!-----------------------------------------------------------------------
!
  INTEGER :: i,j,k,n,imid,jmid,nlevel
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  imid=(nx/2) + 1
  jmid=(ny/2) + 1
!
!-----------------------------------------------------------------------
!
!  Loop over all columns
!
!-----------------------------------------------------------------------
!
  DO j=1,ny-1
    DO i=1,nx-1
!
      DO k=2,nz-1
        n = k-1
        p1d(n)  = 0.01 * pres(i,j,k)
        ht1d(n) = hgt(i,j,k)*1000.
        wmr1d(n)= 1000. * wmr(i,j,k)
        t1d(n)  = tk(i,j,k) - 273.15
        td1d(n) = wmr2td(p1d(n),wmr1d(n))
        tv1d(n) = tctotv(t1d(n),wmr1d(n))
      END DO
!
      nlevel = nz-2
!
      CALL sindex(nz,nlevel,p1d,ht1d,t1d,tv1d,td1d,wmr1d,               &
                  partem,buoy,wload,mbuoy,pbesnd,mbesnd,                &
                  lcl(i,j),lfc(i,j),el(i,j),twdf(i,j),                  &
                  li(i,j),cape(i,j),mcape(i,j),                         &
                  cin(i,j),tcap(i,j))


      IF(i == imid .AND. j == jmid) THEN
        WRITE(6,*)                                                      &
            ' n    p         t      tpar      buoy      wmr       be'
        DO n = 1,nlevel
          WRITE(6,801) n,p1d(n),t1d(n),partem(n),buoy(n),               &
                      wmr1d(n),pbesnd(n)
          801       FORMAT(1X,i2,f8.1,f10.2,f10.2,f10.4,f10.4,f10.1)
        END DO
      END IF

    END DO
  END DO
!
!-----------------------------------------------------------------------
!
!  Smooth the output arrays
!
!-----------------------------------------------------------------------
!
  CALL smooth9p(  lcl,nx,ny,1,nx-1,1,ny-1,tem1)
  CALL smooth9p(  lfc,nx,ny,1,nx-1,1,ny-1,tem1)
  CALL smooth9p(   el,nx,ny,1,nx-1,1,ny-1,tem1)
  CALL smooth9p( twdf,nx,ny,1,nx-1,1,ny-1,tem1)
  CALL smooth9p(   li,nx,ny,1,nx-1,1,ny-1,tem1)
  CALL smooth9p( cape,nx,ny,1,nx-1,1,ny-1,tem1)
  CALL smooth9p(mcape,nx,ny,1,nx-1,1,ny-1,tem1)
  CALL smooth9p(  cin,nx,ny,1,nx-1,1,ny-1,tem1)
  CALL smooth9p( tcap,nx,ny,1,nx-1,1,ny-1,tem1)
!
  RETURN
END SUBROUTINE arps_be
!
!
!##################################################################
!##################################################################
!######                                                      ######
!######                SUBROUTINE SINDEX                     ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE sindex(maxlev,nlevel,p,ht,t,tv,td,w,                         & 1,5
           partem,buoy,wload,mbuoy,pbesnd,mbesnd,                       &
           lcl_pbe,lfc,el,twdf,li,cape,mcape,cin,tcap)
!
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Calculate Convective Available Potential Energy (CAPE)
!  in each column of ARPS grid.
!
!-----------------------------------------------------------------------
!
!
!  AUTHOR: Keith Brewster
!  April, 1995   Based on OLAPS, hence LAPS, version of same.
!
!  MODIFICATION HISTORY:
!  3/11/1996  Cleaned-up, removed OLAPS artifacts.
!  5/13/1996  Added cap strength.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
!
  INTEGER :: maxlev,nlevel
!
  REAL :: p(maxlev),ht(maxlev),t(maxlev),tv(maxlev),                    &
       td(maxlev),w(maxlev)
  REAL :: partem(maxlev),buoy(maxlev),wload(maxlev),mbuoy(maxlev)
  REAL :: pbesnd(maxlev),mbesnd(maxlev)
!
!  Returned from sindx
!
  REAL :: lfc,el,twdf,li,cape,mcape,cin,mcin,tcap
!
!  Potbe variables
!
  REAL :: plcl_pbe,tlcl_pbe,lcl_pbe,thepcl
  REAL :: velneg,mvelneg
!
!  Functions
!
  REAL :: oe
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!
  thepcl=oe(!
!  print *, ' theta-e of parcel: ',thepcl
!
  CALL ptlcl(!
!  Find height of LCL
!
  CALL intrpr(maxlev,nlevel,p,ht,plcl_pbe,lcl_pbe)
!
!  Calculate the CAPE and such
!
  CALL potbe(maxlev,nlevel,p(1),t(1),w(1),                              &
             thepcl,plcl_pbe,tlcl_pbe,lcl_pbe,                          &
             p,ht,t,tv,td,w,                                            &
             partem,buoy,wload,mbuoy,pbesnd,mbesnd,                     &
             cin,velneg,cape,mcin,mvelneg,mcape,lfc,el,tcap)
!
!  Calculate Lifted Index
!
  CALL calcli(maxlev,nlevel,thepcl,p,t,li)
!
!  Calculate max and min wet bulb potential temperature
!
  CALL thwxn(maxlev,nlevel,p,ht,t,td,ht(1),twdf)
!
  RETURN
END SUBROUTINE sindex
!
!##################################################################
!##################################################################
!######                                                      ######
!######                SUBROUTINE POTBE                      ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!
!


SUBROUTINE potbe(maxlev,nlevel,pmean,tmean,wmean,                       & 1,3
           blthte,plcl,tlcl,lcl,                                        &
           p,ht,t,tv,td,w,                                              &
           partem,buoy,wload,mbuoy,pbesnd,mbesnd,                       &
           pbeneg,velneg,pos_max,                                       &
           mbeneg,mvelneg,mpos_max,lfc,el,tcap)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Calculate Convective Available Potential Energy (CAPE)
!  in each column of ARPS grid.
!
!-----------------------------------------------------------------------
!
!
!  AUTHOR: Keith Brewster
!  February, 1994  Based on OLAPS, hence LAPS, version of same.
!                  from FSL, by Steve Albers 1991
!
!  MODIFICATION HISTORY:
!  3/11/1996  Cleaned-up, removed OLAPS artifacts.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
!
  INTEGER :: maxlev,nlevel
  REAL :: pmean,tmean,wmean,blthte,plcl,tlcl,lcl
  REAL :: p(maxlev),ht(maxlev)
  REAL :: t(maxlev),tv(maxlev),td(maxlev),w(maxlev)
  REAL :: partem(maxlev),buoy(maxlev),wload(maxlev),mbuoy(maxlev)
  REAL :: pbesnd(maxlev),mbesnd(maxlev)
  REAL :: pbeneg,velneg,pos_max
  REAL :: mbeneg,mvelneg,mpos_max,lfc,el,tcap
!
!  Parameters
!
  REAL :: g,gamma
  PARAMETER (g=9.80665,                                                 &
             gamma = .009760)   ! Dry Adiabatic Lapse Rate Deg/m
!
!  Functions
!
  REAL :: tsa_fast,tctotv
!
!  Misc internal variables
!
  INTEGER :: n,nel
  REAL :: deltah,delta_ht_dry,delta_ht_wet
  REAL :: sntlcl,buoy_lcl,wsat,partv
  REAL :: nbe_min,pbe_wet,pbe_dry,pos_area
  REAL :: wlow,htzero,adjeng
!
!-----------------------------------------------------------------------
!
!  Function f_mrsat and inline directive for Cray PVP
!
!-----------------------------------------------------------------------
!
  REAL :: f_mrsat

!fpp$ expand (f_mrsat)
!dir$ inline always f_mrsat
!*$*  inline routine (f_mrsat)

!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!
!  Reset output variables.
!
!  These should be the same as what is assigned
!  no positive area found.  They are limited in
!  range to allow for contouring when positive
!  areas exist in some columns of a domain and
!  not in others.
!
  pbeneg=-400.
  velneg=20.
  pos_max=0.
  mbeneg=-400.
  mvelneg=20.
  mpos_max=0.
  lfc=10000.
  el=0.
  tcap=0.
!
!  Inxtialize parcel path arrays
!
  partem(1) = t(1)
  buoy(1) = 0.
  wload(1) = 0.
  mbuoy(1) = 0.
  pbesnd(1) = 0.
  mbesnd(1) = 0.

!  WRITE(6,810)pmean,tmean,wmean,plcl,tlcl,lcl
! 810 format(' pmean,tmean,wmean,plcl,tlcl,lcl',2F10.2,F10.5,2F10.2
!    +   ,F5.1)

  DO n=2,nlevel
    deltah = ht(n) - ht(n-1)
    IF(plcl < p(n-1))THEN ! lower level is below LCL
      IF(plcl < p(n))THEN ! upper level is below LCL
!        WRITE(6,*)' DRY CASE'
        partem(n)=partem(n-1)-gamma*deltah
        partv=tctotv(partem(n),w(1))
        buoy(n)=(partv-tv(n))/tv(n)
        pbesnd(n)=pbesnd(n-1)+g*0.5*(buoy(n)+buoy(n-1))*deltah
        wload(n)=0.
        mbuoy(n)=buoy(n)
        mbesnd(n)=pbesnd(n)
        IF((p(1)-p(n)) < 300.) tcap=AMAX1(tcap,(tv(n)-partv))

      ELSE ! Upper level is above LCL
!
!  BRACKETING CASE AROUND lcl - DRY ADIABATIC PART
!
!        WRITE(6,*)' DRY ADIABATIC PART'
        delta_ht_dry = lcl - ht(n-1)
!        WRITE(6,307)tlcl
!307        format(' PARCEL TEMP AT lcl= ',F10.3)
        CALL intrpr(maxlev,nlevel,p,tv,plcl,sntlcl)
        partv=tctotv(tlcl,w(1))
        buoy_lcl=(partv-sntlcl)/sntlcl
        pbe_dry=g*0.5*(buoy_lcl+buoy(n-1))*delta_ht_dry
        IF((p(1)-plcl) < 300.) tcap=AMAX1(tcap,(sntlcl-partv))
!        WRITE(6,777)N,P(N),tlcl,sntlcl,buoy_lcl
!#          ,buoy(N-1),delta_ht_dry,HT(N),pbesnd(N-1)+pbe_dry
!
!        MOIST ADIABATIC PART
!
!        WRITE(6,*)' MOIST ADIABATIC PART'
        delta_ht_wet=deltah-delta_ht_dry

        partem(n) = tsa_fast(blthte,p(n))
        wsat=1000.*f_mrsat( p(n)*100., partem(n)+273.15 )
        partv=tctotv(partem(n),wsat)
        buoy(n)=(partv-tv(n))/tv(n)
        pbe_wet = g*0.5*(buoy(n)+buoy_lcl)*delta_ht_wet
        pbesnd(n)=pbesnd(n-1) + pbe_dry + pbe_wet
!
        wload(n)=0.001*(w(1)-wsat)
        mbuoy(n)=buoy(n) - wload(n)
        pbe_wet = g*0.5*(mbuoy(n)+buoy_lcl)*delta_ht_wet
        mbesnd(n)=mbesnd(n-1) + pbe_dry + pbe_wet
        IF((p(1)-plcl) < 300.) tcap=AMAX1(tcap,(tv(n)-partv))

      END IF ! Upper level below LCL (Dry or bracket)
    ELSE ! Lower Level is above LCL
!      WRITE(6,*)' GETTING PARCEL TEMPERATURE FOR MOIST CASE'
      partem(n) = tsa_fast(blthte,p(n))

      wsat=1000.*f_mrsat( p(n)*100., partem(n)+273.15 )
      partv=tctotv(partem(n),wsat)
      buoy(n)=(partv-tv(n))/tv(n)
      pbesnd(n)=pbesnd(n-1)+g*0.5*(buoy(n)+buoy(n-1))*deltah
!
      wload(n)=0.001*(w(1)-wsat)
      mbuoy(n)=buoy(n) - wload(n)
      mbesnd(n)=mbesnd(n-1)+g*0.5*(mbuoy(n)+mbuoy(n-1))*deltah
      IF((p(1)-p(n)) < 300.) tcap=AMAX1(tcap,(tv(n)-partv))

    END IF

!    WRITE(6,777)N,P(N),partem(N),T(N),(buoy(n)*1000.),pbesnd(n)
!777    format(' PBE: P,partem,t,b,pbe=',I3,F6.1,4F8.2)
  END DO
!
!  DETERMINE ENERGY EXTREMA
!  Find heights with nuetral buoyancy
!
  pos_area=0.
  nbe_min=0.
  DO n=2,nlevel
!    WRITE(6,940)N
!940    format(
!    :' LOOKING FOR NEUTRAL BUOYANCY - ENERGY EXTREMUM, LEVEL',I3)

    IF((buoy(n)*buoy(n-1)) < 0.)THEN
      wlow=buoy(n)/(buoy(n)-buoy(n-1))
      htzero=ht(n)*(1.-wlow) + wlow*ht(n-1)
      deltah=htzero-ht(n-1)
      adjeng=pbesnd(n-1)+g*0.5*buoy(n-1)*deltah
!
      IF (p(n) >= 500.)  THEN
        nbe_min=AMIN1(adjeng,nbe_min)
      END IF
!
      pos_area=adjeng-nbe_min
      pos_max=AMAX1(pos_area,pos_max)
    END IF
  END DO

!  WRITE(6,464)ICP,ICT,N1,NLEVEL
!464  format(' ICP,ICT,N1,NLEVEL',4I5)
!
!  Case when equlibrium level is above top of domain
!
  pos_area=pbesnd(nlevel)-nbe_min
  pos_max=AMAX1(pos_area,pos_max)
!
!  At least one region of positive area in sounding
!  Make sure there is at least 1 J/kg to avoid some
!  round-off errors esp near LCL.
!
  IF(pos_max > 1.0)THEN
    pbeneg=AMAX1(nbe_min,-400.)
    velneg=SQRT(2.0*ABS(pbeneg))
    velneg=AMIN1(velneg,20.)
  ELSE ! Case when no positive area exists anywhere in sounding
    pos_max=0.0
    pbeneg =-400.
    velneg = 20.
  END IF
!  WRITE(6,485)pos_max,PBENEG,VELNEG
!485  format(' pos_max',F10.1,' PBENEG',F10.1,' VELNEG',F10.1)

!
!  DETERMINE ENERGY EXTREMA FOR MOIST BUOYANCY
!  Find heights with nuetral buoyancy
!
  pos_area=0.
  nbe_min=0.
  DO n=2,nlevel
!    WRITE(6,940)N

    IF((mbuoy(n)*mbuoy(n-1)) < 0.)THEN
      wlow=mbuoy(n)/(mbuoy(n)-mbuoy(n-1))
      htzero=ht(n)*(1.-wlow) + wlow*ht(n-1)
      deltah=htzero-ht(n-1)
      adjeng=mbesnd(n-1)+g*0.5*mbuoy(n-1)*deltah
!
      IF (p(n) >= 500.)  THEN
        nbe_min=AMIN1(adjeng,nbe_min)
      END IF
!
      pos_area=adjeng-nbe_min
      mpos_max=AMAX1(pos_area,mpos_max)
    END IF
  END DO

!  WRITE(6,464)ICP,ICT,N1,NLEVEL
!
!  Case when equlibrium level is above top of domain
!
  pos_area=mbesnd(nlevel)-nbe_min
  mpos_max=AMAX1(pos_area,mpos_max)
!
!  At least one region of positive area in sounding
!  Make sure there is at least 1 J/kg to
!  spurious pos energy due to round off.
!
  IF(mpos_max > 1.0)THEN
    mbeneg=AMAX1(nbe_min,-400.)
    mvelneg=SQRT(2.0*ABS(pbeneg))
    mvelneg=AMIN1(mvelneg,20.)
  ELSE ! Case when no positive area exists anywhere in sounding
    mpos_max=0.0
    mbeneg =-400.
    mvelneg = 20.
  END IF
!  WRITE(6,486)mpos_max,PBENEG,VELNEG
!486  format(' Mpos_max',F10.1,' MBENEG',F10.1,' mVELNEG',F10.1)
!
!    Case when equlibrium level is above top of domain
!
  mpos_max = MAX(mpos_max,(mbesnd(nlevel) - nbe_min))
!
!  Find EL and LFC
!  Unxts are set to km ASL
!
  IF(pos_max > 1.0) THEN
    IF(buoy(nlevel) > 0.) THEN
      nel=nlevel
      el=0.001*ht(nlevel)
    ELSE
      DO  n=nlevel-1,2,-1
        IF(buoy(n) > 0.) EXIT
      END DO
!      1201     CONTINUE
      nel=n
      wlow=buoy(n+1)/(buoy(n+1)-buoy(n))
      el=0.001 * (ht(n+1)*(1.-wlow) + ht(n)*wlow)
    END IF
!
    DO n=nel,1,-1
      IF(buoy(n) < 0.) EXIT
    END DO
!    1301   CONTINUE
    IF(n > 0) THEN
      wlow=buoy(n+1)/(buoy(n+1)-buoy(n))
      lfc=ht(n+1)*(1.-wlow) + ht(n)*wlow
    ELSE
      lfc=ht(1)
    END IF
  ELSE
    el=0.
    lfc=10000.
  END IF
  lfc=AMIN1(lfc,10000.)
  RETURN
END SUBROUTINE potbe
!
!
!##################################################################
!##################################################################
!######                                                      ######
!######                 FUNCTION TCTOTV                      ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!
!


  FUNCTION tctotv(tt,ww)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Virtual Temperature
!
!  Given T in Celcius and mixing ratio in g/kg
!  find the virtual temperature.
!
!-----------------------------------------------------------------------
!
!
!  AUTHOR: Keith Brewster
!  April, 1995   Based on OLAPS, hence LAPS, version of same.
!
!  MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
!
  REAL :: tctotv,tt,ww
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  tctotv=(tt+273.15)*(1.+0.0006*ww)
  RETURN
  END FUNCTION tctotv
!
!##################################################################
!##################################################################
!######                                                      ######
!######                 SUBROUTINE THWXN                     ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!
!


SUBROUTINE thwxn(maxlev,nlevel,p,ht,t,td,elev,twdf) 1
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!
!-----------------------------------------------------------------------
!
!
!  AUTHOR: Keith Brewster
!  April, 1995   Based on OLAPS, hence LAPS, version of same.
!
!  MODIFICATION HISTORY:
!  3/11/1996  Cleaned-up, removed OLAPS artifacts.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
!
!  Input variables
!
  INTEGER :: maxlev,nlevel
  REAL :: p(maxlev),ht(maxlev),t(maxlev),td(maxlev)
  REAL :: elev
!
!  Output variables
!
  REAL :: twdf
!
!  Functions
!
  REAL :: oe,tsa_fast
!
!  Misc internal variables
!
  INTEGER :: n
  REAL :: h3km,thaec,thw,twx,twn
!
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  twx=-999.
  twn=999.
  h3km=elev+3000.
  DO n=1,nlevel
    IF(ht(n) >= elev) THEN
      IF(ht(n) > h3km) EXIT
      thaec=oe(      thw=tsa_fast(thaec,1000.)
      twx=AMAX1(twx,thw)
      twn=AMIN1(twn,thw)
    END IF
  END DO
!  101 CONTINUE
!
!  Find difference between max and min
!
  twdf=twx-twn
  RETURN
END SUBROUTINE thwxn
!
!##################################################################
!##################################################################
!######                                                      ######
!######                SUBROUTINE CALCLI                     ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE calcli(maxlev,nlevel,thepcl,p,t,li) 1
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Calculate Convective Available Potential Energy (CAPE)
!  in each column of ARPS grid.
!
!-----------------------------------------------------------------------
!
!
!  AUTHOR: Keith Brewster
!  April, 1995   Based on OLAPS, hence LAPS, version of same.
!
!  MODIFICATION HISTORY:
!  3/11/1996  Cleaned-up, removed OLAPS artifacts.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
!
!  Input variables
!
  INTEGER :: maxlev,nlevel
  REAL :: thepcl
  REAL :: p(maxlev),t(maxlev)
!
!  Output variable
!
  REAL :: li
!
!  Functions
!
  REAL :: tsa_fast
!
!  Misc internal variables
!
  INTEGER :: n
  REAL :: dp,wlow,t500,par500
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!
  DO n=2,nlevel-1
    IF( p(n) <= 500.) EXIT
  END DO
!  101 CONTINUE
  dp=ALOG(  wlow=ALOG(500./p(n))/dp
  t500=t(n)*(1.-wlow) + t(n-1)*wlow
  par500=tsa_fast(thepcl,500.)
!
  li=t500-par500
!
  RETURN
END SUBROUTINE calcli
!
!##################################################################
!##################################################################
!######                                                      ######
!######                 FUNCTION WMR2TD                      ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!
!


  FUNCTION wmr2td(pres,wmr)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Calculate Convective Available Potential Energy (CAPE)
!  in each column of ARPS grid.
!
!-----------------------------------------------------------------------
!
!
!  AUTHOR: Keith Brewster
!  April, 1995   Based on GEMPAK routine of same name.
!
!  MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
  REAL :: wmr2td
  REAL :: pres
  REAL :: wmr
  REAL :: wkgkg,e,evap
!
  wkgkg = 0.001 * wmr
  wkgkg = AMAX1(wkgkg,0.00005)
  e= (pres*wkgkg) / (0.62197 + wkgkg)
  evap = e /(1.001 + (( pres - 100.) /900.) * 0.0034)
  wmr2td = ALOG(evap/6.112) * 243.5 /( 17.67 - ALOG (evap/6.112))

  RETURN
  END FUNCTION wmr2td
!
!##################################################################
!##################################################################
!######                                                      ######
!######                SUBROUTINE INTRPR                     ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!
!


SUBROUTINE intrpr(nz,nlev,p,var,plvl,varatp) 2
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Interpolate variable "var" linearly in log-pressure (p).
!
!-----------------------------------------------------------------------
!
!
!  AUTHOR: Keith Brewster
!  April, 1995   Based on OLAPS, hence LAPS, version of same.
!
!  MODIFICATION HISTORY:
!  3/11/1996  Cleaned-up, removed OLAPS artifacts.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
  INTEGER :: nz,nlev
  REAL :: p(nz),var(nz)
  REAL :: plvl
  REAL :: varatp
!
  INTEGER :: k
  REAL :: w1
!
  DO k=2,nlev-1
    IF(  END DO
!  101 CONTINUE
!
  w1=ALOG(  varatp = w1*var(k-1) + (1.-w1)*var(k)
!
  RETURN
END SUBROUTINE intrpr

!
!##################################################################
!##################################################################
!######                                                      ######
!######                SUBROUTINE CALCSHR                    ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE calcshr(nx,ny,nz,x,y,zp,sigma,                               & 2,11
           p_pa,t3d,u3d,v3d,cape,                                       &
           shr37,ustrm,vstrm,srlfl,srmfl,helicity,brn,brnu,blcon,       &
           tem1,tem2,tem3,elev,hgt3d)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Calculate various wind shear parameters useful for gauging
!  the potential for severe storms.
!
!-----------------------------------------------------------------------
!
!
!  AUTHOR: Keith Brewster
!  February, 1994   Based on OLAPS, hence LAPS, version of same.
!
!  MODIFICATION HISTORY:
!  3/11/1996  Keith Brewster
!  Added storm-relative flows, general clean-up to
!  meet ARPS coding standards.
!
!  3/22/1996  (Keith Brewster)
!  Fixed some bugs, added smoothing at the end.
!
!  06/20/2000 (Eric Kemp and Keith Brewster)
!  Changed BRN Shear to be the denominator of BRN, instead of wind
!  speed, now has units of speed squared.
!
!-----------------------------------------------------------------------
!
!  Calculates some of the shear related variables from the
!  ARPS 3D wind field.
!
!  shr37         Magnitude of wind shear between 3 and 7 km AGL
!  ustrm,vstrm   Estimated storm motion (modified from Bob Johns)
!  srlfl         Low-level storm-relative wind
!  srmfl         Mid-level storm-relative wind
!  helicity      Helicity, storm relative
!  brn           Bulk Richardson Number (Weisman and Klemp)
!  brnu          Shear parameter of BRN, "U"
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
!
!-----------------------------------------------------------------------
!
!  Input variables
!
!-----------------------------------------------------------------------
!
  INTEGER :: nx,ny,nz
  REAL :: x(nx)
  REAL :: y(ny)
  REAL :: zp(nx,ny,nz)
  REAL :: sigma(nx,ny)
  REAL :: p_pa(nx,ny,nz)
  REAL :: t3d(nx,ny,nz)
  REAL :: u3d(nx,ny,nz)
  REAL :: v3d(nx,ny,nz)
  REAL :: cape(nx,ny)
!
!-----------------------------------------------------------------------
!
!  Output variables
!
!-----------------------------------------------------------------------
!
  REAL :: shr37(nx,ny)      ! 7km - 3km wind shear
  REAL :: ustrm(nx,ny)
  REAL :: vstrm(nx,ny)  ! Estimated storm motion (Bob Johns)
  REAL :: srlfl(nx,ny)
  REAL :: srmfl(nx,ny)
  REAL :: helicity(nx,ny)   ! Helicity, storm relative
  REAL :: brn(nx,ny)        ! Bulk Richardson Number (Weisman and Klemp)
  REAL :: brnu(nx,ny)       ! Shear parameter of BRN, "U"
  REAL :: blcon(nx,ny)
!
!-----------------------------------------------------------------------
!
!  Temporary variables
!
!-----------------------------------------------------------------------
!
  REAL :: tem1(nx,ny)
  REAL :: tem2(nx,ny)
  REAL :: tem3(nx,ny)
  REAL :: elev(nx,ny)
  REAL :: hgt3d(nx,ny,nz)
!
!-----------------------------------------------------------------------
!
!  Misc internal variables
!
!-----------------------------------------------------------------------
!
  INTEGER :: imid,jmid
  INTEGER :: i,j,k,ksfc,k2km,k3km
  REAL :: h3km,u3km,v3km,h7km,u7km,v7km,u2,v2
  REAL :: p2km,t2km,h2km,u2km,v2km,p9km,t9km,h9km,u9km,v9km
  REAL :: p500m,t500m,h500m,u500m,v500m,p6km,t6km,h6km,u6km,v6km
  REAL :: sumu,sumv,sump,sumh,wlow,whigh,dx,dy,dz,dp,dx2,dy2
  REAL :: rhohi,rholo,rhoinv
  REAL :: dirmean,spmean,ushr,vshr,ddir,perc
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  imid=nx/2
  jmid=nx/2
!
  DO j=1,ny-1
    DO i=1,nx-1
      elev(i,j)=zp(i,j,2)
    END DO
  END DO

  DO k=1,nz-1
    DO j=1,ny-1
      DO i=1,nx-1
        hgt3d(i,j,k)=0.5*(zp(i,j,k)+zp(i,j,k+1))
      END DO
    END DO
  END DO
!
!-----------------------------------------------------------------------
!
!  Main loop over all horizontal scalar points
!
!-----------------------------------------------------------------------
!
  ksfc=2
  DO j=1,ny-1
    DO i=1,nx-1
!
!-----------------------------------------------------------------------
!
!  Find mass weighted mean wind in first 500m
!
!-----------------------------------------------------------------------
!
      sumu=0.
      sumv=0.
      sump=0.
      h500m=elev(i,j)+500.
      DO k=ksfc+1,nz-1
        IF( hgt3d(i,j,k) < h500m ) THEN
          dp=p_pa(i,j,k-1)-p_pa(i,j,k)
          rhohi=p_pa(i,j,k)/t3d(i,j,k)
          rholo=p_pa(i,j,k-1)/t3d(i,j,k-1)
          rhoinv=1./(rhohi+rholo)
          sumu=sumu+dp*rhoinv*(rholo*u3d(i,j,k-1)+rhohi*u3d(i,j,k))
          sumv=sumv+dp*rhoinv*(rholo*v3d(i,j,k-1)+rhohi*v3d(i,j,k))
          sump=sump+dp
        ELSE
          dz=hgt3d(i,j,k)-hgt3d(i,j,k-1)
          wlow=(hgt3d(i,j,k)-h500m)/dz
          u500m=u3d(i,j,k)*(1.-wlow) + u3d(i,j,k-1)*wlow
          v500m=v3d(i,j,k)*(1.-wlow) + v3d(i,j,k-1)*wlow
          p500m=p_pa(i,j,k)*(1.-wlow) + p_pa(i,j,k-1)*wlow
          t500m=t3d(i,j,k)*(1.-wlow) + t3d(i,j,k-1)*wlow
          dp=p_pa(i,j,k-1)-p500m
          rhohi=p500m/t500m
          rholo=p_pa(i,j,k-1)/t3d(i,j,k-1)
          rhoinv=1./(rhohi+rholo)
!        print *, ' p500,t500 = ',p500m,t500m
!        print *, ' u500,v500 = ',u500m,v500m
          sumu=sumu+dp*rhoinv*(rholo*u3d(i,j,k-1)+rhohi*u500m)
          sumv=sumv+dp*rhoinv*(rholo*v3d(i,j,k-1)+rhohi*v500m)
          sump=sump+dp
          EXIT
        END IF
      END DO
!      121    CONTINUE
      u500m=sumu/sump
      v500m=sumv/sump
!
!-----------------------------------------------------------------------
!
!  Find mass weighted mean wind sfc-2km AGL
!
!-----------------------------------------------------------------------
!
      sumu=0.
      sumv=0.
      sump=0.
      h2km=elev(i,j)+2000.
      DO k=ksfc+1,nz-1
        IF( hgt3d(i,j,k) < h2km ) THEN
          dp=p_pa(i,j,k-1)-p_pa(i,j,k)
          rhohi=p_pa(i,j,k)/t3d(i,j,k)
          rholo=p_pa(i,j,k-1)/t3d(i,j,k-1)
          rhoinv=1./(rhohi+rholo)
          sumu=sumu+dp*rhoinv*(rholo*u3d(i,j,k-1)+rhohi*u3d(i,j,k))
          sumv=sumv+dp*rhoinv*(rholo*v3d(i,j,k-1)+rhohi*v3d(i,j,k))
          sump=sump+dp
        ELSE
          dz=hgt3d(i,j,k)-hgt3d(i,j,k-1)
          wlow=(hgt3d(i,j,k)-h2km)/dz
          u2km=u3d(i,j,k)*(1.-wlow) + u3d(i,j,k-1)*wlow
          v2km=v3d(i,j,k)*(1.-wlow) + v3d(i,j,k-1)*wlow
          p2km=p_pa(i,j,k)*(1.-wlow) + p_pa(i,j,k-1)*wlow
          t2km=t3d(i,j,k)*(1.-wlow) + t3d(i,j,k-1)*wlow
          dp=p_pa(i,j,k-1)-p2km
          rhohi=p2km/t2km
          rholo=p_pa(i,j,k-1)/t3d(i,j,k-1)
          rhoinv=1./(rhohi+rholo)
          sumu=sumu+dp*rhoinv*(rholo*u3d(i,j,k-1)+rhohi*u2km)
          sumv=sumv+dp*rhoinv*(rholo*v3d(i,j,k-1)+rhohi*v2km)
          sump=sump+dp
          EXIT
        END IF
      END DO
!      141    CONTINUE
      u2km=sumu/sump
      v2km=sumv/sump
      tem2(i,j)=u2km
      tem3(i,j)=v2km
!
!-----------------------------------------------------------------------
!
!  Find mass weighted mean wind sfc-6km AGL
!
!-----------------------------------------------------------------------
!
      sumu=0.
      sumv=0.
      sump=0.
      h6km=elev(i,j)+6000.
      DO k=ksfc+1,nz-1
        IF( hgt3d(i,j,k) < h6km ) THEN
          dp=p_pa(i,j,k-1)-p_pa(i,j,k)
          rhohi=p_pa(i,j,k)/t3d(i,j,k)
          rholo=p_pa(i,j,k-1)/t3d(i,j,k-1)
          rhoinv=1./(rhohi+rholo)
          sumu=sumu+dp*rhoinv*(rholo*u3d(i,j,k-1)+rhohi*u3d(i,j,k))
          sumv=sumv+dp*rhoinv*(rholo*v3d(i,j,k-1)+rhohi*v3d(i,j,k))
          sump=sump+dp
        ELSE
          dz=hgt3d(i,j,k)-hgt3d(i,j,k-1)
          wlow=(hgt3d(i,j,k)-h6km)/dz
          u6km=u3d(i,j,k)*(1.-wlow) + u3d(i,j,k-1)*wlow
          v6km=v3d(i,j,k)*(1.-wlow) + v3d(i,j,k-1)*wlow
          p6km=p_pa(i,j,k)*(1.-wlow) + p_pa(i,j,k-1)*wlow
          t6km=t3d(i,j,k)*(1.-wlow) + t3d(i,j,k-1)*wlow
          dp=p_pa(i,j,k-1)-p6km
          rhohi=p6km/t6km
          rholo=p_pa(i,j,k-1)/t3d(i,j,k-1)
          rhoinv=1./(rhohi+rholo)
          sumu=sumu+dp*rhoinv*(rholo*u3d(i,j,k-1)+rhohi*u6km)
          sumv=sumv+dp*rhoinv*(rholo*v3d(i,j,k-1)+rhohi*v6km)
          sump=sump+dp
          EXIT
        END IF
      END DO
      u6km=sumu/sump
      v6km=sumv/sump
!
!-----------------------------------------------------------------------
!
!  Find mass weighted mean wind 2km-9km AGL
!
!-----------------------------------------------------------------------
!
      sumu=0.
      sumv=0.
      sump=0.
      h9km=elev(i,j)+9000.
      DO k=ksfc+1,nz-2
        IF( hgt3d(i,j,k) > h2km ) EXIT
      END DO
!      181   CONTINUE
      k2km=k
      dz=hgt3d(i,j,k)-hgt3d(i,j,k-1)
      wlow=(hgt3d(i,j,k)-h2km)/dz
      u2=u3d(i,j,k)*(1.-wlow) + u3d(i,j,k-1)*wlow
      v2=v3d(i,j,k)*(1.-wlow) + v3d(i,j,k-1)*wlow
      p2km=p_pa(i,j,k)*(1.-wlow) + p_pa(i,j,k-1)*wlow
      t2km=t3d(i,j,k)*(1.-wlow) + t3d(i,j,k-1)*wlow
      dp=p2km-p_pa(i,j,k)
      rholo=p2km/t2km
      rhohi=p_pa(i,j,k)/t3d(i,j,k)
      rhoinv=1./(rhohi+rholo)
      sumu=sumu+dp*rhoinv*(rholo*u2+rhohi*u3d(i,j,k))
      sumv=sumv+dp*rhoinv*(rholo*v2+rhohi*v3d(i,j,k))
      sump=sump+dp
      DO k=k2km+1,nz-1
        IF( hgt3d(i,j,k) < h9km ) THEN
          dp=p_pa(i,j,k-1)-p_pa(i,j,k)
          rhohi=p_pa(i,j,k)/t3d(i,j,k)
          rholo=p_pa(i,j,k-1)/t3d(i,j,k-1)
          rhoinv=1./(rhohi+rholo)
          sumu=sumu+dp*rhoinv*(rholo*u3d(i,j,k-1)+rhohi*u3d(i,j,k))
          sumv=sumv+dp*rhoinv*(rholo*v3d(i,j,k-1)+rhohi*v3d(i,j,k))
          sump=sump+dp
        ELSE
          dz=hgt3d(i,j,k)-hgt3d(i,j,k-1)
          wlow=(hgt3d(i,j,k)-h9km)/dz
          u9km=u3d(i,j,k)*(1.-wlow) + u3d(i,j,k-1)*wlow
          v9km=v3d(i,j,k)*(1.-wlow) + v3d(i,j,k-1)*wlow
          p9km=p_pa(i,j,k)*(1.-wlow) + p_pa(i,j,k-1)*wlow
          t9km=t3d(i,j,k)*(1.-wlow) + t3d(i,j,k-1)*wlow
          dp=p_pa(i,j,k-1)-p9km
          rhohi=p9km/t9km
          rholo=p_pa(i,j,k-1)/t3d(i,j,k-1)
          rhoinv=1./(rhohi+rholo)
          sumu=sumu+dp*rhoinv*(rholo*u3d(i,j,k-1)+rhohi*u9km)
          sumv=sumv+dp*rhoinv*(rholo*v3d(i,j,k-1)+rhohi*v9km)
          sump=sump+dp
          EXIT
        END IF
      END DO
!      191    CONTINUE
      u9km=sumu/sump
      v9km=sumv/sump
!
!-----------------------------------------------------------------------
!
!  Storm motion estimation
!  From Davies and Johns, 1993
!  "Some wind and instability parameters associated With
!  strong and violent tornadoes."
!  AGU Monograph 79, The Tornado...(Page 575)
!
!  Becuase of the discontinuity produced by that method
!  at the 15.5 m/s cutoff, their rules have been modified
!  to provide a gradual transition, and accomodate all the
!  data they mention in the article.
!
!-----------------------------------------------------------------------
!
      CALL uv2ddff(u6km,v6km,dirmean,spmean)
      IF(spmean >= 20.0) THEN
        dirmean=dirmean+18.
        IF(dirmean > 360.) dirmean=dirmean-360.
        spmean=spmean*0.89
      ELSE IF (spmean > 8.0) THEN
        whigh=(spmean - 8.0)/12.
        wlow =1.-whigh
        ddir=wlow*32.0 + whigh*18.0
        perc=wlow*0.75 + whigh*0.89
        dirmean=dirmean+ddir
        IF(dirmean > 360.) dirmean=dirmean-360.
        spmean=spmean*perc
      ELSE
        dirmean=dirmean+32.
        IF(dirmean > 360.) dirmean=dirmean-360.
        spmean=spmean*0.75
      END IF
      CALL ddff2uv(dirmean,spmean,ustrm(i,j),vstrm(i,j))
!
!-----------------------------------------------------------------------
!
!  Storm-relative low-level flow
!
!-----------------------------------------------------------------------
!
      srlfl(i,j)=SQRT((ustrm(i,j)-u2km)*(ustrm(i,j)-u2km) +             &
                      (vstrm(i,j)-v2km)*(vstrm(i,j)-v2km))
!
!-----------------------------------------------------------------------
!
!  Storm relative mid-level flow
!
!-----------------------------------------------------------------------
!
      srmfl(i,j)=SQRT((ustrm(i,j)-u9km)*(ustrm(i,j)-u9km) +             &
                      (vstrm(i,j)-v9km)*(vstrm(i,j)-v9km))
!
!-----------------------------------------------------------------------
!
!  Shear parameter for Bulk Richardson number
!
!-----------------------------------------------------------------------
!
!    print *, ' density-weight mean 0-500 m ',u500m,v500m
!    print *, ' density-weight mean 0-6  km ',u6km,v6km
!
      brnu(i,j)=0.5*( (u6km-u500m)*(u6km-u500m) +                       &
                      (v6km-v500m)*(v6km-v500m) )
!
!-----------------------------------------------------------------------
!
!  Bulk Richardson number
!  A limit of 200 is imposed, since this could
!  go to inifinity.
!
!-----------------------------------------------------------------------
!
      IF(brnu(i,j) > 0.) THEN
        brn(i,j)=cape(i,j)/brnu(i,j)
        brn(i,j)=AMIN1(brn(i,j),200.)
      ELSE
        brn(i,j)=200.
      END IF
!
    END DO
  END DO
!
!-----------------------------------------------------------------------
!
!  Calculate Helicity and 3km to 7km shear
!  since both involve the 3km wind.
!
!  For more efficient computation the Helicity is
!  computed for zero storm motion and the storm
!  motion is accounted for by adding a term at the end.
!  This is mathematically equivalent to accounting
!  for the storm motion at each level.
!
!-----------------------------------------------------------------------
!
  DO j=1,ny-1
    DO i=1,nx-1
      h3km=elev(i,j)+3000.
      h7km=elev(i,j)+7000.
!
!-----------------------------------------------------------------------
!
!  Find level just above 3 km AGL
!  Note, it is assumed here that there is at least
!  one level between the sfc and 3 km.
!
!-----------------------------------------------------------------------
!
      sumh=0.
      DO k=ksfc+1,nz-2
        IF(hgt3d(i,j,k) > h3km) EXIT
        sumh=sumh +                                                     &
            ( u3d(i,j,k)*v3d(i,j,k-1) ) -                               &
            ( v3d(i,j,k)*u3d(i,j,k-1) )
      END DO
!      240    CONTINUE
      k3km=k
      dz=hgt3d(i,j,k)-hgt3d(i,j,k-1)
      wlow=(hgt3d(i,j,k)-h3km)/dz
      u3km=u3d(i,j,k)*(1.-wlow) + u3d(i,j,k-1)*wlow
      v3km=v3d(i,j,k)*(1.-wlow) + v3d(i,j,k-1)*wlow
      sumh=sumh +                                                       &
            ( u3km*v3d(i,j,k-1) ) -                                     &
            ( v3km*u3d(i,j,k-1) )
      ushr=u3km-u3d(i,j,2)
      vshr=v3km-v3d(i,j,2)
      helicity(i,j)=sumh + vshr*ustrm(i,j) - ushr*vstrm(i,j)
!
!-----------------------------------------------------------------------
!
!  Now Find 7km wind for 3-to-7km shear
!
!-----------------------------------------------------------------------
!
      DO k=k3km,nz-1
        IF(hgt3d(i,j,k) > h7km) EXIT
      END DO
!      260   CONTINUE
      dz=hgt3d(i,j,k)-hgt3d(i,j,k-1)
      wlow=(hgt3d(i,j,k)-h7km)/dz
      u7km=u3d(i,j,k)*(1.-wlow) + u3d(i,j,k-1)*wlow
      v7km=v3d(i,j,k)*(1.-wlow) + v3d(i,j,k-1)*wlow
      shr37(i,j)=(SQRT( (u7km-u3km)*(u7km-u3km) +                       &
                        (v7km-v3km)*(v7km-v3km) ))/4000.
    END DO
  END DO

!-----------------------------------------------------------------------
!
!  Calculate the low-level convergence
!
!-----------------------------------------------------------------------
!
  dx=(x(2)-x(1))
  dy=(y(2)-y(1))
  dx2=2.*dx
  dy2=2.*dy
  DO j=2,ny-2
    DO i=2,nx-2
      blcon(i,j)=-1000.*(                                               &
                 (tem2(i+1,j)-tem2(i-1,j))/(sigma(i,j)*dx2)+            &
                 (tem3(i,j+1)-tem3(i,j-1))/(sigma(i,j)*dy2) )
    END DO
  END DO

  DO j=2,ny-2
    blcon(1,j)=-1000.*(                                                 &
                (tem2(2,j)-tem2(1,j))/(sigma(1,j)*dx)+                  &
                (tem3(1,j+1)-tem3(1,j-1))/(sigma(1,j)*dy2) )
    blcon(nx-1,j)=-1000.*(                                              &
           (tem2(nx-1,j)-tem2(nx-2,j))/(sigma(nx-1,j)*dx)+              &
           (tem3(nx-1,j+1)-tem3(nx-1,j-1))/(sigma(nx-1,j)*dy2) )
  END DO
  DO i=2,nx-2
    blcon(i,1)=-1000.*(                                                 &
                (tem2(i+1,1)-tem2(i-1,1))/(sigma(i,1)*dx2)+             &
                (tem3(i,2)-tem3(i,1))/(sigma(i,1)*dy) )
    blcon(i,ny-1)=-1000.*(                                              &
            (tem2(i+1,ny-1)-tem2(i-1,ny-1))/(sigma(i,ny-1)*dx2)+        &
            (tem3(i,ny-1)-tem3(i,ny-2))/(sigma(i,ny-1)*dy) )
  END DO
  blcon(1,1)=blcon(2,2)
  blcon(nx-1,ny-1)=blcon(nx-2,ny-2)
!
!-----------------------------------------------------------------------
!
!  Smooth the output arrays
!
!-----------------------------------------------------------------------
!
  CALL smooth9p(   shr37,nx,ny,1,nx-1,1,ny-1,tem1)
  CALL smooth9p(   ustrm,nx,ny,1,nx-1,1,ny-1,tem1)
  CALL smooth9p(   vstrm,nx,ny,1,nx-1,1,ny-1,tem1)
  CALL smooth9p(   srlfl,nx,ny,1,nx-1,1,ny-1,tem1)
  CALL smooth9p(   srmfl,nx,ny,1,nx-1,1,ny-1,tem1)
  CALL smooth9p(helicity,nx,ny,1,nx-1,1,ny-1,tem1)
  CALL smooth9p(     brn,nx,ny,1,nx-1,1,ny-1,tem1)
  CALL smooth9p(    brnu,nx,ny,1,nx-1,1,ny-1,tem1)
  CALL smooth9p(   blcon,nx,ny,1,nx-1,1,ny-1,tem1)
!
  RETURN
END SUBROUTINE calcshr
!
!##################################################################
!##################################################################
!######                                                      ######
!######                FUNCTION TSA_FAST                     ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


  FUNCTION tsa_fast(os,p),1
!
!   THIS FUNCTION RETURNS THE TEMPERATURE TSA (CELSIUS) ON A SATURATION
!   ADIABAT AT PRESSURE P (MILLIBARS). OS IS THE EQUIVALENT POTENTIAL
!   TEMPERATURE OF THE PARCEL (CELSIUS). SIGN(A,B) REPLACES THE
!   ALGEBRAIC SIGN OF A WITH THAT OF B.
!
!    BAKER,SCHLATTER 17-MAY-1982     Original version
!    Modification for better convergence, Keith Brewster, Feb 1994.
!
!   B IS AN EMPIRICAL CONSTANT APPROXIMATELY EQUAL TO THE LATENT HEAT
!   OF VAPORIZATION FOR WATER DIVIDED BY THE SPECIFIC HEAT AT CONSTANT
!   PRESSURE FOR DRY AIR.
!
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Calculate Convective Available Potential Energy (CAPE)
!  in each column of ARPS grid.
!
!-----------------------------------------------------------------------
!
!
!  AUTHOR: Keith Brewster
!  April, 1995   Based on OLAPS, hence LAPS, version of same.
!
!  MODIFICATION HISTORY:
!  3/11/1996  Cleaned-up, removed OLAPS artifacts.
!
!-----------------------------------------------------------------------
!
  REAL :: b
!  PARAMETER (B=2.6518986)
  PARAMETER (b=2651.8986)
  a= os+273.15
!
!   Above 200 mb figure all the moisture is wrung-out, so
!   the temperature is that which has potential temp of theta-e.
!   Otherwise iterate to find combo of moisture and temp corresponding
!   to thetae.
!
  IF( p < 200.) THEN
    tq=a*((p/1000.)**.286)
  ELSE
!   D IS AN INITIAL VALUE USED IN THE ITERATION BELOW.
    d= 120.
!   TQ IS THE FIRST GUESS FOR TSA.
    tq= 253.15
    x = 0.
!
!   ITERATE TO OBTAIN SUFFICIENT ACCURACY....SEE TABLE 1, P.8
!   OF STIPANUK (1973) FOR EQUATION USED IN ITERATION.
    DO i= 1,25
      d= 0.5*d

      x_last = x

      x= a*EXP(-b*f_mrsat(p*100.,tq)/tq)-tq*((1000./p)**.286)

      IF (ABS(x) < 1E-3) GO TO 2 
!
      IF (x_last * x < 0.) THEN
        slope = (x-x_last) / (tq - tq_last)
        delta = - x / slope
        ad = AMIN1(ABS(delta),d)
        tq_last = tq
        tq = tq + SIGN(ad,delta)
      ELSE
        tq_last = tq
        tq= tq+SIGN(d,x)
      END IF

    END DO
  END IF
2 tsa_fast = tq-273.15
  RETURN
  END FUNCTION tsa_fast
!
!##################################################################
!##################################################################
!######                                                      ######
!######                SUBROUTINE UV2DDFF                    ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE uv2ddff(u,v,dd,ff) 1
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Calculate direction and speed from u and v.
!
!-----------------------------------------------------------------------
!
!
!  AUTHOR: Keith Brewster
!  3/11/1996
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
  REAL :: u,v,dd,ff
  REAL :: dlon
  REAL :: r2deg
  PARAMETER (r2deg=180./3.141592654)
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!
  ff = SQRT(u*u + v*v)

  IF(v > 0.) THEN
    dlon=r2deg*ATAN(u/v)
  ELSE IF(v < 0.) THEN
    dlon=180. + r2deg*ATAN(u/v)
  ELSE IF(u >= 0.) THEN
    dlon=90.
  ELSE
    dlon=-90.
  END IF

  dd= dlon + 180.
  dd= dd-360.*(nint(dd)/360)
  RETURN
END SUBROUTINE uv2ddff
!
!
!##################################################################
!##################################################################
!######                                                      ######
!######                SUBROUTINE DDFF2UV                    ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE ddff2uv(dd,ff,u,v) 1
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Calculate u and v wind components from direction and speed.
!
!-----------------------------------------------------------------------
!
!
!  AUTHOR: Keith Brewster
!  3/11/1996
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
  REAL :: u,v,dd,ff
  REAL :: arg
  REAL :: d2rad
  PARAMETER (d2rad=3.141592654/180.)
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  arg = (dd * d2rad)
  u = -ff * SIN(arg)
  v = -ff * COS(arg)
  RETURN
END SUBROUTINE ddff2uv