
! --------------------------------------------------------------------------
!  PROGRAMA gera_sondagens
!
!           Autor: ERNANI DE LIMA NASCIMENTO
!           Data: OUTUBRO de 2006
!
!  Este programa lê arquivo binário do grads gerado pelo script SOUND3.EXEC
!  e gera um arquivo de saída ASCII contendo as 3264 sondagens do ETA para o caso
!  do tornado de Indaiatuba. A formatação desta saída satisfaz a leitura para
!  o programa indices_severos5.f90
! ---------------------------------------------------------------------------
  program gera_sondagens
! Variáveis de trabalho
  integer i,j,k,irec,geop,umrl,station,istart,nlev
  dimension station(3270),istart(3270),nlev(3270)
! Variáveis de entrada
  real :: variab,lev
  dimension  variab(20,9,3270),lev(20)
!
! Início do programa executável...
!
! Lendo o arquivo de entrada contendo as saída binária do GrADS:
! O arquivo binário gerado pelo SOUND3.EXEC para o caso Indaiatuba tem 9 variáveis
! à supercície mais 9 variáveis em 19 níveis de pressão em matrizes de 64 X 51 pontos.
!
  open (2, file='sonds_2400_2418.le.bin',FORM='UNFORMATTED',status=   &
        'UNKNOWN',ACCESS='DIRECT',recl=3264*9*20*4)
 irec=1
 read (2,rec=irec) (((variab(i,j,k),k=1,3264),j=1,9),i=1,20)
! i= 20 níveis (19 níveis de pressão mais o nível de superfície)
! j= 9 variáveis
! k= 3264 pontos da matriz (64 x 51)
!
! Agora checamos o número de níveis de cada perfil. Isto vai depender de quantos
! níveis estão "abaixo" do solo em cada ponto de grade. A sondagem vai sempre
! começar do nível de superfície, não do nível de 1000hPa.
!
 do k=1,3264
  if (variab(1,2,k) > 1000.0) then
   nlev(k)=20
   istart(k)=2
  else
   if (variab(1,2,k) > 925.0) then
    nlev(k)=19
    istart(k)=3
   else
    if (variab(1,2,k) > 900.0) then
     nlev(k)=18
     istart(k)=4
    else
     if (variab(1,2,k) > 850.0) then
      nlev(k)=17
      istart(k)=5
     else
      if (variab(1,2,k) > 800.0) then
       nlev(k)=16
       istart(k)=6
      else
       if (variab(1,2,k) > 750.0) then
        nlev(k)=15
        istart(k)=7
       else
        if (variab(1,2,k) > 700.0) then
         nlev(k)=14
         istart(k)=8
        else
         if (variab(1,2,k) > 650.0) then
          nlev(k)=13
          istart(k)=9
         else
          if (variab(1,2,k) > 600.0) then
           nlev(k)=12
           istart(k)=10
          else
           if (variab(1,2,k) > 550.0) then
            nlev(k)=11
            istart(k)=11
           else
            nlev(k)=0
           endif
          endif
         endif
        endif
       endif
      endif
     endif
    endif
   endif
  endif
 enddo
 lev(2)=1000.0
 lev(3)=925.0
 lev(4)=900.0
 do i=5,20
  lev(i)=lev(i-1)-50.0
 enddo
 geop=0
 umrl=0
! Definindo po arquivo de saída contendo as sondagens para serem
  open (3, file='eta.sound_2400_2418.txt', FORM='FORMATTED',status='UNKNOWN')
! A LINHA ACIMA TEM QUE SER EDITADA TODA VEZ QUE MUDARMOS O ARQUIVO DE ENTRADA.
!
  station(1)=75000
  write(3,*)' 3264'
  do k=1,3264   ! Loop pelo número de pontos de grade.
   write(3,*) station(k),' ',nlev(k),' 18 24 05 2005'
! A LINHA ACIMA TEM QUE SER EDITADA TODA VEZ QUE MUDARMOS O ARQUIVO DE SONDAGEM.
   if (nlev(k) == 0) then
    print*,' Ponto sem sondagem: ', k
   else
    geop=nint(variab(1,1,k))
    umrl=nint(variab(1,5,k))
! Escreve o primeiro nível do ponto de grade k (sempre o nível de superfície).
    write(3,301) variab(1,2,k),geop,variab(1,3,k),variab(1,4,k),  &
           umrl,variab(1,6,k),variab(1,7,k),variab(1,8,k),variab(1,9,k)
! Loop pelo número de níveis do ponto de grade k.
    do i=istart(k),20
     geop=nint(variab(i,1,k))
     umrl=nint(variab(i,5,k))
     write(3,301) lev(i),geop,variab(i,3,k),variab(i,4,k),  &
           umrl,variab(i,6,k),variab(i,7,k),variab(i,8,k),variab(i,9,k)
    enddo
    station(k+1)=station(k)+1
   endif
  enddo
!300  format (1X,I5,3I4,I5)
301  format (1X,F7.1,I7,2F7.1,I7,F7.2,2F8.2,F9.3)
305  format (1X,I3)
  stop
! END OF PROGRAM GERA SONDAGENS
  end



! --------------------------------------------------------------------------
!  PROGRAMA índices_severos5
!
!           Autor: ERNANI DE LIMA NASCIMENTO
!           Data: Outubro de 2006
!
!  Este programa calcula índices de tempo severo a partir de perfis
!  verticais atmosféricos gerados por modelos numéricos.
!  NOTA 1: ESTA VERSÃO CALCULA ÍNDICES SEVEROS TERMODINÂMICOS E CINEMÁTICOS.
!  NOTA 2: ESTA VERSÃO É IDÊNTICA À VERSÃO índices_severos4, EXCETO QUE LÊ
!          UM ARQUIVO DE ENTRADA COM CONTEÚDO DIFERENTE E GERA UM ARQUIVO
!          BINÁRIO PARA O GRADS, ALÉM DO ARQUIVO ASCII PARA O ARPS.
!  NOTA 3: ESTA VERSÃO GERA SAÍDA DE SONDAGENS EM FORMATO PARA EXECUÇÃO NO
!          ARPS. POR EXEMPLO, LÊ O ARQUIVO eta.sound_0212_0318.txt E GERA
!          O ARQUIVO arps.eta.sound_0212_0318.txt. DEPOIS ESTE ARQUIVO
!          DEVE SER EDITADO PARA INCLUIR O CABEÇALHO DESEJADO E ENTÃO
!          SER SALVO COM EXTENSÃO .snd PARA SER UTILIZADO COM O ARPS
!          (LEIA O ARQUIVO LEIA-ME PARA MAIORES DETALHES).
!
!
!  Exemplo de arquivo de entrada necessário para este programa (este formato
!  é aquele gerado pelo programa gera_sond_2.f90):
!
! 78122          20 18 24 05 2005
!  1000.0    125   25.4   19.7     70  14.65   -1.24   -7.64  298.581
!   925.0    804   20.3   16.4     79  12.83   -1.26   -7.91  300.023
!   900.0   1040   18.3   15.0     81  11.97   -1.31   -8.07  300.422
!   850.0   1528   15.4   11.6     78  10.18   -2.30  -10.65  302.267
!   800.0   2041   13.0    7.5     69   8.17   -2.44  -10.54  305.057
!   750.0   2581    9.8    3.8     67   6.72   -1.68   -7.88  307.192
!   700.0   3151    6.3    0.2     65   5.54   -0.90   -5.52  309.413
!   650.0   3757    3.1   -4.2     59   4.30   -0.38   -4.89  312.487
!
!   ...    ....    ...   ...     ...   ...    ...     ...   ....   ....  ....
!  etc.....
!
! O que significa cada número no exemplo acima:
! 78122 => identificador da "estação" (irrelevante neste programa)
! 20 => número de níveis do perfil
! 18 => horário UTC
! 24 => dia
! 05 => mês
! 2005 => ano
! Colunas: pressão (hPa), altitude (m), temperatura (C), temperatura do ponto
! de orvalho (C), umidade relativa (%), razão de mistura (g/kg), componente zonal
! do vento (m/s), componente meridional do vento (m/s), temperatura potencial (K).
!
!
! Índices calculados: CAPE de superfície, CAPE com correção de water loading,
! indice de levantamento, nível de condensação por levantamento, nível de convecção
! espontânea, inibição convectiva, componentes u e v do movimento esperado de tempestades,
! escoamento relativo à tempestade em níveis baixos e médios, helicidade relativa à
! tempestade, número de Richardson volumétrico, denominador do número de Richardson
! volumétrico, índice de energia-helicidade.
!
! Com exceção da rotina de cálculo do índice de energia-helicidade (interiramente nova),
! todas as rotinas foram modificadas a partir do ARPS5.0_IHOP6 e ADAS.
!
!
! ---------------------------------------------------------------------------
  program indices_severos
! Variáveis miscelâneas
  integer i,j,k
! Variáveis de entrada
  integer :: nrad,estacao,hora,dia,mes,ano
  integer :: altm,urporc,dddeg,vvkt,nlev,strm_opt
  dimension :: estacao(3270),hora(3270),dia(3270),mes(3270),ano(3270),  &
               nlev(3270)
  dimension  altm(102,3270),urporc(102,3270)
!  dddeg(100,3270),vvkt(100,3270) vvmps(100,3270)
  real :: presshpa,tc,tdc,wvgkg,thetae,thetav,presspa,tk,wvkgkg, &
          u,v,arg,u2dd,v2dd,theta,specif
  dimension  presshpa(102,3270),tc(102,3270),tdc(102,3270),wvgkg(102,3270), &
             presspa(102,3270),tk(102,3270),wvkgkg(102,3270),u(102,3270),   &
             v(102,3270),theta(102,3270),specif(102,3270)
! Variáveis para os índices de tempo severo
  real :: cape,mcape,li,ncl,nce,cin,ustrm,vstrm,llsrm,mlsrm,heli,nrv,  &
          dnrv,ieh,sup,tcap,thetab,bl,blcon,shr37,ffstrm,ddstrm,dnrv2km
  dimension  cape(3270),mcape(3270),li(3270),ncl(3270),nce(3270),cin(3270),el(3270), &
             ustrm(3270),vstrm(3270),llsrm(3270),mlsrm(3270),heli(3270),nrv(3270), &
             dnrv(3270),ieh(3270),tcap(3270),thetab(3270),bl(3270),shr37(3270),    &
             blcon(3270),sup(3270),ffstrm(3270),ddstrm(3270),dnrv2km(3270)
! Obs: tcap = lid strength; thetab = twdf on arpsplt.f (max wet bulb potential
! temperature difference)
!
! Variáveis para cálculos em coluna
  real :: p1d,ht1d,t1d,tv1d,td1d,wmr1d,partem,buoy,wload,mbuoy,pbesnd, &
          mbesnd
  dimension  p1d(102),ht1d(102),t1d(102),tv1d(102),td1d(102),          &
             wmr1d(102),partem(102),buoy(102),wload(102),mbuoy(102),   &
             pbesnd(102),mbesnd(102)
! Variável para nome do arquivo de entrada
  character(60) inname,outname1,outname2
!
! Início do programa executável...
!
! Lendo o arquivo de entrada contendo as sondagens:
  write(*,*) 'Este programa calcula índices de tempo severo.'
  write(*,*) 'Ernani L. Nascimento (OUT/2006) (modificado do ARPS). '
  write(*,*) 'Entre com o nome do arquivo de entrada (até 60 caracteres):'
  read(*,*) inname
  write(*,*) 'Utilizar deslocamento observado da tempestade? (1=sim 5=não)'
  read(*,*) strm_opt
  if ((strm_opt /= 1).AND.(strm_opt /= 5)) then
   write(*,*) 'Opção não disponível. Programa interrompido.'
   stop
  endif
  outname1='grads.'//inname(1:42)
  outname2='arps.'//inname(1:42)
  open (10,file=inname,status='old')
  read(10,*) nrad
  do k=1,nrad
   read (10,*) estacao(k),nlev(k),hora(k),dia(k),mes(k),ano(k)
   if (nlev(k) /= 0) then
    do i=1,nlev(k)
     read (10,*) presshpa(i,k),altm(i,k),tc(i,k),tdc(i,k),urporc(i,k),   &
                wvgkg(i,k),u(i,k),v(i,k),theta(i,k)
     presspa(i,k)=presshpa(i,k)*100.
     tk(i,k)=tc(i,k)+273.15
     specif(i,k)=(wvgkg(i,k)/1000.)/(1+(wvgkg(i,k)/1000.))
!      print *,'specif(1,k)= ',specif(1,k), k
     arg=0.0
    end do
!     write(*,*) k,' ',presspa(1,k)

! As linhas abaixo geram 2 níveis adicionais na sondagem. O ETA possui 19 níveis, mas
! para as rodadas do ARPS tornou-se conveniente acrescentar dois níveis a mais para
! garantir um topo de domínio mais alto. Estas linhas podem ser convenientemente comentadas.
    altm(nlev(k)+1,k)=altm(nlev(k),k)+(altm(nlev(k),k)-altm(nlev(k)-1,k))
    altm(nlev(k)+2,k)=altm(nlev(k)+1,k)+(altm(nlev(k),k)-altm(nlev(k)-1,k))
    specif(nlev(k)+1,k)=specif(nlev(k),k)
    specif(nlev(k)+2,k)=specif(nlev(k),k)
    u(nlev(k)+1,k)=u(nlev(k),k)
    u(nlev(k)+2,k)=u(nlev(k),k)
    v(nlev(k)+1,k)=v(nlev(k),k)
    v(nlev(k)+2,k)=v(nlev(k),k)
   endif
!   print *,'presshpa(19,k)= ',presshpa(19,k),' ', k
  end do
  close(10)
!
!     if (j==1) then
!      write(*,*) 'dddeg(1,k)= ',dddeg(j,k),' para k= ',k
!      write(*,*) 'vvmps(1,k)= ',vvmps(j,k),' para k= ',k
!     endif
!     call ddff2uv(dddeg,vvmps,uu,vv,j)
!
! Transformando graus e magnitude do vento em componentes zonal e meridional
!     arg = (dddeg(j,k) * (3.141592654/180.))
!     uu = -vvmps(j,k) * sin(arg)
!     vv = -vvmps(j,k) * cos(arg)
!     u(j,k)=uu
!     v(j,k)=vv
!     if (j==1) then
!      write(*,*) 'u(1,k)= ',u(j,k),' para k= ',k
!      write(*,*) 'v(1,k)= ',v(j,k),' para k= ',k
!     endif

!
!   Chama rotina de cálculo de índices termodinâmicos
!
!   print *,'hgt(m) antes: ', altm(1,1)
!   do k=1,nrad
!     print *,'presshpa(19,k)= ',presshpa(19,k),' ', k
!   enddo
   write(*,*) k,' PASSEI AQUI 1.'
   call arps_be(nrad,presshpa,altm,tc,tdc,wvgkg,ncl,nce,el,thetab,li,cape,   &
               mcape,cin,tcap,p1d,ht1d,t1d,tv1d,td1d,wmr1d,partem,          &
               buoy,wload,mbuoy,pbesnd,mbesnd,nlev,estacao,hora,dia,        &
               mes,ano)
!
!  Cálculo de índices cinemáticos
!
!
   write(*,*) k,' PASSEI AQUI 2'
   call calcshr(nrad,altm,bl,presspa,tk,u,v,cape,shr37,ustrm,vstrm,llsrm,  &
               mlsrm,heli,nrv,dnrv,dnrv2km,blcon,u1,v1,nlev,estacao,hora,dia,  &
               mes,ano,strm_opt)
   call calcehi_sup(nrad,nlev,presshpa,cape,heli,dnrv,ieh,sup)

!  Arquivo de saída: binário para o grads
  open(12,file=outname1,FORM='UNFORMATTED',status='unknown',ACCESS='DIRECT',   &
       recl=3264*17*4)
!   write(12,*) ' EST UTC DIA MÊS ANO SHR37 ustrm vstrm ddstrm ffstrm helic BRN BRNSHR B2km IEH SUP'
  irec=1
  do k=1,nrad
   u2dd=ustrm(k)
   v2dd=vstrm(k)
   call uv2ddff(u2dd,v2dd,dd,ff)
   ddstrm(k)=dd
   ffstrm(k)=ff
   write(12,rec=irec) ncl(k),nce(k),el(k),thetab(k),li(k),cape(k),mcape(k),cin(k),shr37(k),       &
                ustrm(k),vstrm(k),heli(k),nrv(k),dnrv(k),dnrv2km(k),ieh(k),sup(k)
!    if ((ieh(k)<-1.0).and.(sup(k)<-1.0)) then
!      write(*,*) 'EHI AND SUP LESS THAN -1 (k EHI SUP): ',k,ieh(k),sup(k)
!    endif
  enddo
  close(12)

!  Arquivo de saída: sondagem com formato para o ARPS
!   open(13,file=outname2,status='unknown')
!   write(13,304) nrad
!   do k=1,nrad
!    write(13,300) estacao(k),hora(k),dia(k),mes(k),ano(k)
!    write (13,*) ' ZSND   THSND   QVSND    USND    VSND '
!  As duas linhas abaixo acrescentam dois níveis para theta.
!    theta(nlev(k)+1,k)=theta(nlev(k),k)+(theta(nlev(k),k)-theta(nlev(k)-1,k))
!    theta(nlev(k)+2,k)=theta(nlev(k)+1,k)+(theta(nlev(k),k)-theta(nlev(k)-1,k))
!  O loop abaixo acrescenta os dois níveis adicionais para o ARPS.
!    do i=1,nlev(k)+2
!     write(13,303) float(altm(nlev(k)+3-i,k)),theta(nlev(k)+3-i,k),specif(nlev(k)+3-i,k),u(nlev(k)+3-i,k), &
!                  v(nlev(k)+3-i,k)
!    enddo
!   enddo
!   close(13)
300  format (1X,I5,3I4,I5)
301  format (1X,I5,3I4,I5,F8.1,1X,2F8.1,2F7.1,3F8.1)
302  format (1X,I5,3I4,I5,F8.5,2F7.3,F8.1,F7.1,F9.2,F6.1,2F7.1,2F6.1)
303  format (1X,F9.3,F12.5,F12.8,2F12.5)
304  format (1X,I4)
  stop
! END OF PROGRAM ÍNDICES SEVEROS
  end
!
!
! BELOW: FUNCTIONS AND SUBROUTINES USED BY PROGRAM ÍNDICES SEVEROS
!
!##################################################################
!##################################################################
!######                                                      ######
!######                SUBROUTINE ARPS_BE                    ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######    and modified by Ernani L. Nascimento at Simepar   ######
!######                                                      ######
!##################################################################
!##################################################################
!
SUBROUTINE arps_be(nrad,pres,hgt,tc,tdc,wmr,lcl,lfc,el,twdf,li,cape,mcape,   &
          cin,tcap,p1d,ht1d,t1d,tv1d,td1d,wmr1d,partem,buoy,wload,      &
          mbuoy,pbesnd,mbesnd,nlev,estacao,hora,dia,mes,ano)
!
!-----------------------------------------------------------------------
!
!  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.
!
!  FEB 13-14 2004 (ERNANI L. NASCIMENTO)
!  Modified from ARPS to compute indices from rawinsonde data.
!
!-----------------------------------------------------------------------
!
 IMPLICIT NONE
 INTEGER :: nrad
!
!-----------------------------------------------------------------------
!
!  "1-D" input variables (Ernani L. Nascimento)
!
!-----------------------------------------------------------------------
!
 INTEGER :: nlev(nrad),estacao(nrad),hora(nrad),dia(nrad),mes(nrad),   &
            ano(nrad)
!
!-----------------------------------------------------------------------
!
!  "2-D" input variables (Ernani L. Nascimento)
!
!-----------------------------------------------------------------------
!
 REAL :: pres(102,nrad)
! hgt is integer! Different from the original code: Ernani L. Nascimento
 INTEGER :: hgt(102,nrad)
 REAL :: tc(102,nrad)
 REAL :: wmr(102,nrad)
 REAL :: tdc(102,nrad)
!
!-----------------------------------------------------------------------
!
!  Output variables ("1-D") (Ernani L. Nascimento)
!
!-----------------------------------------------------------------------
!
 REAL :: lcl(nrad)
 REAL :: lfc(nrad)
 REAL :: el(nrad)
 REAL :: twdf(nrad)
 REAL :: li(nrad)
 REAL :: cape(nrad)
 REAL :: mcape(nrad)
 REAL :: cin(nrad)
 REAL :: tcap(nrad)
!
!-----------------------------------------------------------------------
!
!  Scratch space for calculations
!
!-----------------------------------------------------------------------
!
 REAL :: p1d(102),ht1d(102),tv1d(102),td1d(102)
 REAL :: partem(102),buoy(102),wload(102)
 REAL :: mbuoy(102),pbesnd(102),mbesnd(102)
 REAL :: wmr1d(102),t1d(102)
!
!-----------------------------------------------------------------------
!
!  Functions
!
!-----------------------------------------------------------------------
!
 REAL :: wmr2td,tctotv
!
!-----------------------------------------------------------------------
!
!  Misc internal variables
!
!-----------------------------------------------------------------------
!
 INTEGER :: i,j,k,n,nlevel,bla
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!-----------------------------------------------------------------------
!  Loop over all soundings (Ernani L. Nascimento)
!-----------------------------------------------------------------------
!  print *, 'hgt(m) depois: ',hgt(1,1)
 DO k=1,nrad
  if (nlev(k)==0) then
   lcl(k)=-999.9
   lfc(k)=-999.9
   el(k)=-999.9
   twdf(k)=-999.9
   li(k)=-999.9
   cape(k)=-999.9
   mcape(k)=-999.9
   cin(k)=-999.9
   tcap(k)=-999.9
  else
   j=nlev(k)
   if (pres(j,k)>300.0) then
    write(*,*) 'Sondagem abaixo não chegou aos 300hPa (estacao,hora,dia,mes,ano):'
    write(*,*) estacao(k),hora(k),dia(k),mes(k),ano(k)
    write(*,*) pres(j,k),nlev(k)
    write(*,*) 'Buscando próxima sondagem.'
    cycle
   endif
!-----------------------------------------------------------------------
! Loop over all levels of the sounding
!-----------------------------------------------------------------------
   DO i=1,nlev(k)
     p1d(i) = pres(i,k)
     ht1d(i) = float(hgt(i,k))
     wmr1d(i) = wmr(i,k)
     t1d(i) = tc(i,k)
!      td1d(i) = wmr2td(p1d(i),wmr1d(i))
     td1d(i) = tdc(i,k)
     tv1d(i) = tctotv(t1d(i),wmr1d(i))
   END DO
!
!    print *, 'ht1d(1): ', ht1d(1)
   nlevel = nlev(k)
!
   CALL sindex(nlevel,p1d,ht1d,t1d,tv1d,td1d,wmr1d,partem,buoy,wload,    &
               mbuoy,pbesnd,mbesnd,lcl(k),lfc(k),el(k),twdf(k),li(k),    &
               cape(k),mcape(k),cin(k),tcap(k))
  endif
 END DO
RETURN
END SUBROUTINE arps_be
!
!
!##################################################################
!##################################################################
!######                                                      ######
!######                SUBROUTINE SINDEX                     ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######    and modified by Ernank L. Nascimento at Simepar   ######
!######                                                      ######
!##################################################################
!##################################################################
!

SUBROUTINE sindex(nlevel,p,ht,t,tv,td,w,                                &
          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.
!  Feb 13 2004 ERNANI L. NASCIMENTO: removed variable maxlev
!
!-----------------------------------------------------------------------
!
 IMPLICIT NONE
!
 INTEGER :: nlevel
!
 REAL :: p(nlevel+2),ht(nlevel+2),t(nlevel+2),tv(nlevel+2),td(nlevel+2),w(nlevel+2)
 REAL :: partem(nlevel+2),buoy(nlevel+2),wload(nlevel+2),mbuoy(nlevel+2)
 REAL :: pbesnd(nlevel+2),mbesnd(nlevel+2)
!
!  Returned from sindex
!
 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(t(1),td(1),p(1))
!
!  print *, ' theta-e of parcel: ',thepcl
!
 CALL ptlcl(p(1),t(1),td(1),plcl_pbe,tlcl_pbe)
!  print *, ' press and temp at LCL: ',plcl_pbe,tlcl_pbe
!
!  Find height of LCL
!

 CALL intrpr(nlevel,p,ht,plcl_pbe,lcl_pbe)
!  print *, ' NCL: ', lcl_pbe
!
!  Calculate the CAPE and such
!
 CALL potbe(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(nlevel,thepcl,p,t,li)
!
!  Calculate max and min wet bulb potential temperature
!
 CALL thwxn(nlevel,p,ht,t,td,ht(1),twdf)
!
 RETURN
END SUBROUTINE sindex
!
!
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
SUBROUTINE ptlcl(p,t,td,pc,tc)
!
!   this subroutine estimates the pressure pc (mb) and the temperature
!   tc (celsius) at the lifted condensation level (lcl), given the
!   initial pressure p (mb), temperature t (celsius) and dew point
!   (celsius) of the parcel.  the approximation is that lines of
!   constant potential temperature and constant mixing ratio are
!   straight on the skew t/log p chart.
!
!    baker,schlatter   17-may-1982   original version
!
!   teten's formula for saturation vapor pressure as a function of
!   pressure was used in the derivation of the formula below.  for
!   additional details, see math notes by t. schlatter dated 8 sep 81.
!   t. schlatter, noaa/erl/profs program office, boulder, colorado,
!   wrote this subroutine.
!
!   akap = (gas constant for dry air) / (specific heat at constant
!       pressure for dry air)
!   cta = difference between kelvin and celsius temperatures
!
 DATA akap,cta/0.28541,273.16/
 c1 = 4098.026/(td+237.3)**2
 c2 = 1./(akap*(t+cta))
 pc = p*EXP(c1*c2*(t-td)/(c2-c1))
 tc = t+c1*(t-td)/(c2-c1)
 RETURN
END SUBROUTINE ptlcl
!
!
!
!##################################################################
!##################################################################
!######                                                      ######
!######                SUBROUTINE INTRPR                     ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######   and modified by Ernani L. Nascimento at Simepar    ######
!######                                                      ######
!##################################################################
!##################################################################
!
!
SUBROUTINE intrpr(nlev,p,var,plvl,varatp)
!
!-----------------------------------------------------------------------
!
!  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.
!
!  Feb 13 2004  Ernani L. Nascimento
!  Removed variable maxlev
!
!-----------------------------------------------------------------------
!
 IMPLICIT NONE
 INTEGER :: nlev
 REAL :: p(nlev+2),var(nlev+2)
 REAL :: plvl
 REAL :: varatp
!
 INTEGER :: k
 REAL :: w1
!
 DO k=2,nlev
   IF(p(k) < plvl) EXIT
 END DO
!
 w1=ALOG(p(k)/plvl)/ALOG(p(k)/p(k-1))
 varatp = w1*var(k-1) + (1.-w1)*var(k)
!
 RETURN
END SUBROUTINE intrpr
!
!
!##################################################################
!##################################################################
!######                                                      ######
!######                SUBROUTINE POTBE                      ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######   and modified by Ernani L. Nascimento at Simepar    ######
!######                                                      ######
!##################################################################
!##################################################################
!
!
SUBROUTINE potbe(nlevel,pmean,tmean,wmean,                              &
          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.
!
!  Feb/14/2004  Ernani L. Nascimento: removed variable maxlev
!-----------------------------------------------------------------------
!
 IMPLICIT NONE
!
 INTEGER :: nlevel
 REAL :: pmean,tmean,wmean,blthte,plcl,tlcl,lcl
 REAL :: p(nlevel+2),ht(nlevel+2)
 REAL :: t(nlevel+2),tv(nlevel+2),td(nlevel+2),w(nlevel+2)
 REAL :: partem(nlevel+2),buoy(nlevel+2),wload(nlevel+2),mbuoy(nlevel+2)
 REAL :: pbesnd(nlevel+2),mbesnd(nlevel+2)
 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
!
!  Note from Ernani L. Nascimento: this directive is ignored whenever
!  this code is compiled in a machine different from Cray. Thus, no need to
!  delete the piece of code !fpp$, !dir$, !*$* below.
!
!-----------------------------------------------------------------------
!
 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.
!
!  Initialize 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(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
!
!
!
!##################################################################
!##################################################################
!######                                                      ######
!######                SUBROUTINE CALCLI                     ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######  and modified by Ernani L. Nascimento at Simepar     ######
!######                                                      ######
!##################################################################
!##################################################################
!
SUBROUTINE calcli(nlevel,thepcl,p,t,li)
!
!-----------------------------------------------------------------------
!
!  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.
!
!  Fev/14/2004  Ernani L. Nascimento: removed variable maxlev
!
!-----------------------------------------------------------------------
!
 IMPLICIT NONE
!
!  Input variables
!
 INTEGER :: nlevel
 REAL :: thepcl
 REAL :: p(nlevel+2),t(nlevel+2)
!
!  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(p(n-1)/p(n))
 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
!
!
!
!##################################################################
!##################################################################
!######                                                      ######
!######                 SUBROUTINE THWXN                     ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######   and modified by Ernani L. Nascimento at Simepar    ######
!######                                                      ######
!##################################################################
!##################################################################
!
!
SUBROUTINE thwxn(nlevel,p,ht,t,td,elev,twdf)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!
!-----------------------------------------------------------------------
!
!
!  AUTHOR: Keith Brewster
!  April, 1995   Based on OLAPS, hence LAPS, version of same.
!
!  MODIFICATION HISTORY:
!  3/11/1996  Cleaned-up, removed OLAPS artifacts.
!
!  Feb/14/2004 Ernani L. Nascimento: removed variable maxlev
!
!-----------------------------------------------------------------------
!
 IMPLICIT NONE
!
!  Input variables
!
 INTEGER :: nlevel
 REAL :: p(nlevel+2),ht(nlevel+2),t(nlevel+2),td(nlevel+2)
 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(t(n),td(n),p(n))
     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
!
!
!##################################################################
!##################################################################
!######                                                      ######
!######                 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.  => WRONG!!
!
!  CONVERTS WATER VAPOR MIXING RATIO INTO DEW POINT TEMPERATURE
!  (Ernani L. Nascimento)
!
!-----------------------------------------------------------------------
!
!
!  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(wmr,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
!
!
!##################################################################
!##################################################################
!######                                                      ######
!######                 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
!
!
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
 FUNCTION oe(t,td,p)
!
!    g.s. stipanuk     1973          original version.
!    reference stipanuk paper entitled:
!         "algorithms for generating a skew-t, log p
!         diagram and computing selected meteorological
!         quantities."
!         atmospheric sciences laboratory
!         u.s. army electronics command
!         white sands missile range, new mexico 88002
!         33 pages
!    baker, schlatter  17-may-1982

!   this function returns equivalent potential temperature oe (celsius)
!   of a parcel of air given its temperature t (celsius), dew point
!   td (celsius) and pressure p (millibars).
!   find the wet bulb temperature of the parcel.

 atw = tw(t,td,p)

!   find the equivalent potential temperature.

 oe = os(atw,p)
 RETURN
 END FUNCTION oe
!
!
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
 FUNCTION os(t,p)
!
!    g.s. stipanuk     1973          original version.
!    reference stipanuk paper entitled:
!         "algorithms for generating a skew-t, log p
!         diagram and computing selected meteorological
!         quantities."
!         atmospheric sciences laboratory
!         u.s. army electronics command
!         white sands missile range, new mexico 88002
!         33 pages
!    baker, schlatter  17-may-1982

!   this function returns the equivalent potential temperature os
!   (celsius) for a parcel of air saturated at temperature t (celsius)
!   and pressure p (millibars).
 DATA b/2.6518986/
!   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.

 tk = t+273.15
 osk= tk*((1000./p)**.286)*(EXP(b*w(t,p)/tk))
 os= osk-273.15
 RETURN
 END FUNCTION os
!
!
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
 FUNCTION w(t,p)
!
!    g.s. stipanuk     1973              original version.
!    reference stipanuk paper entitled:
!         "algorithms for generating a skew-t, log p
!         diagram and computing selected meteorological
!         quantities."
!         atmospheric sciences laboratory
!         u.s. army electronics command
!         white sands missile range, new mexico 88002
!         33 pages
!    baker, schlatter  17-may-1982
!
!  this function returns the mixing ratio (grams of water vapor per
!  kilogram of dry air) given the dew point (celsius) and pressure
!  (millibars). if the temperture  is input instead of the
!  dew point, then saturation mixing ratio (same units) is returned.
!  the formula is found in most meteorological texts.

 x= esat(t)
 w= 622.*x/(p-x)
 RETURN
 END FUNCTION w
!
!
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
 FUNCTION esat(t)
!
!    g.s. stipanuk     1973           original version.
!    reference stipanuk paper entitled:
!         "algorithms for generating a skew-t, log p
!         diagram and computing selected meteorological
!         quantities."
!         atmospheric sciences laboratory
!         u.s. army electronics command
!         white sands missile range, new mexico 88002
!         33 pages
!    baker, schlatter  17-may-1982
!
!   this function returns the saturation vapor pressure over
!   water (mb) given the temperature (celsius).
!   the algorithm is due to nordquist, w.s.,1973: "numerical approxima-
!   tions of selected meteorlolgical parameters for cloud physics prob-
!   lems," ecom-5475, atmospheric sciences laboratory, u.s. army
!   electronics command, white sands missile range, new mexico 88002.

 tk = t+273.15
 p1 = 11.344-0.0303998*tk
 p2 = 3.49149-1302.8844/tk
 c1 = 23.832241-5.02808*ALOG10(tk)
 esat = 10.**(c1-1.3816E-7*10.**p1+8.1328E-3*10.**p2-2949.076/tk)
 RETURN
 END FUNCTION esat
!
!
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
 FUNCTION tw(t,td,p)

!    g.s. stipanuk     1973           original version.
!    reference stipanuk paper entitled:
!         "algorithms for generating a skew-t, log p
!         diagram and computing selected meteorological
!         quantities."
!         atmospheric sciences laboratory
!         u.s. army electronics command
!         white sands missile range, new mexico 88002
!         33 pages
!    baker, schlatter  17-may-1982

!   this function returns the wet-bulb temperature tw (celsius)
!   given the temperature t (celsius), dew point td (celsius)
!   and pressure p (mb).  see p.13 in stipanuk (1973), referenced
!   above, for a description of the technique.
!
!
!   determine the mixing ratio line thru td and p.

 aw = w(td,p)
!
!   determine the dry adiabat thru t and p.

 ao = o(t,p)
 pi = p

!   iterate to locate pressure pi at the intersection of the two
!   curves .  pi has been set to p for the initial guess.

 DO i= 1,10
   x= .02*(tmr(aw,pi)-tda(ao,pi))
   IF (ABS(x) < 0.01) EXIT
   pi= pi*(2.**(x))
 END DO

!   find the temperature on the dry adiabat ao at pressure pi.

 ti= tda(ao,pi)

!   the intersection has been located...now, find a saturation
!   adiabat thru this point. function os returns the equivalent
!   potential temperature (c) of a parcel saturated at temperature
!   ti and pressure pi.

 aos= os(ti,pi)

!   function tsa returns the wet-bulb temperature (c) of a parcel at
!   pressure p whose equivalent potential temperature is aos.

 tw = tsa(aos,p)
 RETURN
 END FUNCTION tw
!
!
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
 FUNCTION o(t,p)
!
!    g.s. stipanuk     1973          original version.
!    reference stipanuk paper entitled:
!         "algorithms for generating a skew-t, log p
!         diagram and computing selected meteorological
!         quantities."
!         atmospheric sciences laboratory
!         u.s. army electronics command
!         white sands missile range, new mexico 88002
!         33 pages
!    baker, schlatter  17-may-1982

!   this function returns potential temperature (celsius) given
!   temperature t (celsius) and pressure p (mb) by solving the poisson
!   equation.

 tk= t+273.15
 ok= tk*((1000./p)**.286)
 o= ok-273.15
 RETURN
 END FUNCTION o
!
!
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
 FUNCTION tmr(w,p)
!
!    g.s. stipanuk     1973           original version.
!    reference stipanuk paper entitled:
!         "algorithms for generating a skew-t, log p
!         diagram and computing selected meteorological
!         quantities."
!         atmospheric sciences laboratory
!         u.s. army electronics command
!         white sands missile range, new mexico 88002
!         33 pages
!    baker, schlatter  17-may-1982

!   this function returns the temperature (celsius) on a mixing
!   ratio line w (g/kg) at pressure p (mb). the formula is given in
!   table 1 on page 7 of stipanuk (1973).
!
!   initialize constants

 DATA c1/.0498646455/,c2/2.4082965/,c3/7.07475/
 DATA c4/38.9114/,c5/.0915/,c6/1.2035/

 x= ALOG10(w*p/(622.+w))
 tmrk= 10.**(c1*x+c2)-c3+c4*((10.**(c5*x)-c6)**2.)
 tmr= tmrk-273.15
 RETURN
 END FUNCTION tmr
!
!
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
 FUNCTION tda(o,p)
!
!    g.s. stipanuk     1973           original version.
!    reference stipanuk paper entitled:
!         "algorithms for generating a skew-t, log p
!         diagram and computing selected meteorological
!         quantities."
!         atmospheric sciences laboratory
!         u.s. army electronics command
!         white sands missile range, new mexico 88002
!         33 pages
!    baker, schlatter  17-may-1982

!   this function returns the temperature tda (celsius) on a dry adiabat
!   at pressure p (millibars). the dry adiabat is given by
!   potential temperature o (celsius). the computation is based on
!   poisson's equation.

 ok= o+273.15
 tdak= ok*((p*.001)**.286)
 tda= tdak-273.15
 RETURN
 END FUNCTION tda
!
!
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
 FUNCTION tsa(os,p)
!
!    g.s. stipanuk     1973           original version.
!    reference stipanuk paper entitled:
!         "algorithms for generating a skew-t, log p
!         diagram and computing selected meteorological
!         quantities."
!         atmospheric sciences laboratory
!         u.s. army electronics command
!         white sands missile range, new mexico 88002
!         33 pages
!    baker, schlatter  17-may-1982

!   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.
!   b is an empirical constant approximately equal to 0.001 of the latent
!   heat of vaporization for water divided by the specific heat at constant
!   pressure for dry air.

 DATA b/2.6518986/
 a= os+273.15

!   tq is the first guess for tsa.

 tq= 253.15

!   d is an initial value used in the iteration below.

 d= 120.

!   iterate to obtain sufficient accuracy....see table 1, p.8
!   of stipanuk (1973) for equation used in iteration.

 DO i= 1,12
   tqk= tq-273.15
   d= d/2.
   x= a*EXP(-b*w(tqk,p)/tq)-tq*((1000./p)**.286)
   IF (ABS(x) < 1E-7) GOTO 2
   tq= tq+SIGN(d,x)
 END DO
2 tsa= tq-273.15
 RETURN
 END FUNCTION tsa
!
!
!
!##################################################################
!##################################################################
!######                                                      ######
!######                FUNCTION TSA_FAST                     ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!
FUNCTION tsa_fast(os,p)
!
!   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
!
!
!##################################################################
!##################################################################
!######                                                      ######
!######                   FUNCTION F_MRSAT                   ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!
FUNCTION f_mrsat( p, t )
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Calculate the saturation water vapor mixing ratio using enhanced
!  Teten's formula.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Yuhe Liu
!  01/08/1998
!
!  MODIFICATION HISTORY:
!
!  16 FEB 2004
!  Ernani L. Nascimento: modified to include constant rddrv without
!  INCLUDE command.
!
!-----------------------------------------------------------------------
!
!  INPUT :
!
!    p        Pressure (Pascal)
!    t        Temperature (K)
!
!  OUTPUT:
!
!    f_mrsat  Saturation water vapor mixing ratio (kg/kg).
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
 IMPLICIT NONE

 REAL :: p         ! Pressure (Pascal)
 REAL :: t         ! Temperature (K)
 REAL :: f_mrsat   ! Saturation water vapor mixing ratio (kg/kg)

 REAL :: rd        ! Gas constant for dry air  (m**2/(s**2*K))
 PARAMETER( rd     = 287.0 )

 REAL :: rv        ! Gas constant for water vapor  (m**2/(s**2*K)).
 PARAMETER( rv     = 461.0 )

 REAL :: rddrv
 PARAMETER( rddrv  = rd/rv )
!
!-----------------------------------------------------------------------
!
!  Misc. local variables:
!
!-----------------------------------------------------------------------
!
 REAL :: fes
!
!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------
!
! Removed line below. ERNANI L. NASCIMENTO
!  INCLUDE 'phycst.inc'
!
!-----------------------------------------------------------------------
!
!  Function f_es and inline directive for Cray PVP
!
!-----------------------------------------------------------------------
!
 REAL :: f_es
!fpp$ expand (f_es)
!dir$ inline always f_es
!*$*  inline routine (f_es)
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
 fes = f_es( p,t )
 f_mrsat = rddrv * fes / (p-fes)

 RETURN
END FUNCTION f_mrsat
!
!
!
!##################################################################
!##################################################################
!######                                                      ######
!######                   FUNCTION F_ES                      ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!
FUNCTION f_es( p, t )
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Calculate the saturation specific humidity using enhanced Teten's
!  formula.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Yuhe Liu
!  01/08/1998
!
!  MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
!  INPUT :
!
!    p        Pressure (Pascal)
!    t        Temperature (K)
!
!  OUTPUT:
!
!    f_es     Saturation water vapor pressure (Pa)
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
 IMPLICIT NONE

 REAL :: p         ! Pressure (Pascal)
 REAL :: t         ! Temperature (K)
 REAL :: f_es      ! Saturation water vapor pressure (Pa)
!
!-----------------------------------------------------------------------
!
!  Function f_es and inline directive for Cray PVP
!
!-----------------------------------------------------------------------
!
 REAL :: f_esl, f_esi

!fpp$ expand (f_esl)
!fpp$ expand (f_esi)
!dir$ inline always f_esl, f_esi
!*$*  inline routine (f_esl, f_esi)

!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
 IF ( t >= 273.15 ) THEN      ! for water
   f_es = f_esl( p,t )
 ELSE                            ! for ice
   f_es = f_esi( p,t )
 END IF

 RETURN
END FUNCTION f_es
!
!
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
FUNCTION f_esl( p, t )
!-----------------------------------------------------------------------
!  Calculate the saturation water vapor over liquid water using
!  enhanced Teten's formula.
!
!  Feb 16 2004:
!  Modified by Ernani L. Nascimento to include constants satfwa,satfwb,
!  satewa,satewb,satewc without INCLUDE command.
!-----------------------------------------------------------------------
 IMPLICIT NONE

 REAL :: p         ! Pressure (Pascal)
 REAL :: t         ! Temperature (K)
 REAL :: f_esl     ! Saturation water vapor pressure over liquid water

 REAL :: f

 REAL :: satfwa, satfwb
 PARAMETER ( satfwa = 1.0007 )
 PARAMETER ( satfwb = 3.46E-8 )  ! for p in Pa

 REAL :: satewa, satewb, satewc
 PARAMETER ( satewa = 611.21 )   ! es in Pa
 PARAMETER ( satewb = 17.502 )
 PARAMETER ( satewc = 32.18 )

!  Removed line below. ERNANI L. NASCIMENTO
!  INCLUDE 'globcst.inc'
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
 f = satfwa + satfwb * p
 f_esl = f * satewa * EXP( satewb*(t-273.15)/(t-satewc) )

 RETURN
END FUNCTION f_esl
!
!
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
FUNCTION f_esi( p, t )
!-----------------------------------------------------------------------
!  Calculate the saturation water vapor over ice using enhanced
!  Teten's formula.
!
!  Feb 16 2004:
!  Modified by Ernani L. Nascimento to include constants satfia,satfib,
!  sateia,sateib,sateic without INCLUDE command.
!-----------------------------------------------------------------------
 IMPLICIT NONE

 REAL :: p         ! Pressure (Pascal)
 REAL :: t         ! Temperature (K)
 REAL :: f_esi     ! Saturation water vapor pressure over ice (Pa)

 REAL :: f

 REAL :: satfia, satfib
 PARAMETER ( satfia = 1.0003 )
 PARAMETER ( satfib = 4.18E-8 )  ! for p in Pa

 REAL :: sateia, sateib, sateic
 PARAMETER ( sateia = 611.15 )   ! es in Pa
 PARAMETER ( sateib = 22.452 )
 PARAMETER ( sateic = 0.6 )

!  Removed line below. ERNANI L. NASCIMENTO
!  INCLUDE 'globcst.inc'
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
 f = satfia + satfib * p
 f_esi = f * sateia * EXP( sateib*(t-273.15)/(t-sateic) )

 RETURN
END FUNCTION f_esi

!
!
!     ##################################################################
!     ##################################################################
!     ######                                                      ######
!     ######                SUBROUTINE DDFF2UV                    ######
!     ######                                                      ######
!     ######                     Developed by                     ######
!     ######     Center for Analysis and Prediction of Storms     ######
!     ######                University of Oklahoma                ######
!     ######                                                      ######
!     ##################################################################
!     ##################################################################
!
     SUBROUTINE ddff2uv(dd,ff,u,v)
!
!#######################################################################
!
!     PURPOSE:
!
!     Calculate u and v wind components from direction and speed.
!
!#######################################################################
!
!
!     AUTHOR: Keith Brewster
!     3/11/1996
!
!     09/feb/2004 (ERNANI L. NASCIMENTO)
!       THIS SUBROUTINE IS EXACTLY LIKE SUBROUTINE  ddff2uv IN
!     /src/adas/thermo3d.f  I´VE JUST CHANGED THE NAME.
!
!#######################################################################
!
  implicit none
  integer j,k
  real :: u,v,dd,ff
!   dimension  u(100,35),v(100,35),dd(100,35),ff(100,35)
  real :: arg
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!     Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!   arg = (dd(j,k) * (3.141592654/180.))
  arg = (dd * (3.141592654/180.))
!   if (m==1) then
!    write(*,*) 'Dentro do ddff2uv: ',dd,ff,arg
!   endif
!   u(j,k) = -ff(j,k) * sin(arg)
  u = -ff * sin(arg)
!   v(j,k) = -ff(j,k) * cos(arg)
  v = -ff * cos(arg)
  RETURN
END SUBROUTINE ddff2uv
!
!
!##################################################################
!##################################################################
!######                                                      ######
!######                SUBROUTINE UV2DDFF                    ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!

SUBROUTINE uv2ddff(u,v,dd,ff)
!
!-----------------------------------------------------------------------
!
!  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 CALCSHR                    ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!

SUBROUTINE calcshr(nrad,zp,sigma,                                       &
          p_pa,t3d,u3d,v3d,cape,                                       &
          shr37,ustrm,vstrm,srlfl,srmfl,helicity,brn,brnu,brnu2km,blcon, &
          tem2,tem3,nlev,estacao,hora,dia,mes,ano,strm_opt)
!
!-----------------------------------------------------------------------
!
!  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.
!
!  04/19/2004 (Ernani L. Nascimento)
!  Modified the code to fit in índices_severos3.f90, and adapted
!  the computation of expected storm-motion for the Southern Hemisphere.
!
!-----------------------------------------------------------------------
!
!  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
 INTEGER :: nrad,strm_opt
!
!-----------------------------------------------------------------------
!
!  "1-D" input variables (Ernani L. Nascimento)
!
!-----------------------------------------------------------------------
!
 INTEGER :: nlev(nrad),estacao(nrad),hora(nrad),dia(nrad),mes(nrad),   &
            ano(nrad)
 REAL :: sigma(nrad)
 REAL :: cape(nrad)
!
!-----------------------------------------------------------------------
!
!  "2-D" input variables (Ernani L. Nascimento)
!
!-----------------------------------------------------------------------
!
 INTEGER :: zp(102,nrad)
 REAL :: p_pa(102,nrad)   ! Pressure in Pascals
 REAL :: t3d(102,nrad)    ! Temperature in Kelvin
 REAL :: u3d(102,nrad)
 REAL :: v3d(102,nrad)

!
!-----------------------------------------------------------------------
!
!  Output variables
!
!-----------------------------------------------------------------------
!
 REAL :: shr37(nrad)      ! 7km - 3km wind shear
 REAL :: ustrm(nrad)
 REAL :: vstrm(nrad)  ! Estimated storm motion (Bob Johns)
 REAL :: srlfl(nrad)
 REAL :: srmfl(nrad)
 REAL :: helicity(nrad)   ! Helicity, storm relative
 REAL :: brn(nrad)        ! Bulk Richardson Number (Weisman and Klemp)
 REAL :: brnu(nrad)       ! Shear parameter of BRN, "U"
 REAL :: brnu2km(nrad)    ! 2km shear parameter of BRN
 REAL :: blcon(nrad)
!
!-----------------------------------------------------------------------
!
!  Temporary variables
!
!-----------------------------------------------------------------------
!

 REAL :: tem2(nrad)
 REAL :: tem3(nrad)
!
!-----------------------------------------------------------------------
!
!  Misc internal variables
!
!-----------------------------------------------------------------------
!
!  INTEGER :: imid,jmid
 REAL :: elev(nrad),hgt3d(102,nrad)
 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,arg,new_ustrm,new_vstrm
 REAL :: dirmean,spmean,ushr,vshr,ddir,perc,obs_ddstrm,obs_ffstrm
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!-----------------------------------------------------------------------
!  Loop over all soundings (Ernani L. Nascimento)
!-----------------------------------------------------------------------
!
 write(*,*) 'Entrada da subrotina calcshr'
 DO k=1,nrad
  if (nlev(k)==0) then
   shr37(k)=-999.9
   ustrm(k)=-999.9
   vstrm(k)=-999.9
   srlfl(k)=-999.9
   srmfl(k)=-999.9
   helicity(k)=-999.9
   brn(k)=-999.9
   brnu(k)=-999.9
   brnu2km(k)=-999.9
   blcon(k)=-999.9
  else
   print*, estacao(k),hora(k),dia(k),mes(k),ano(k)
   j=nlev(k)
   if ((p_pa(j,k)*0.01)>300.0) then
    cycle
   endif

   arg=0.0
   elev(k)=real(zp(1,k))
!    write(*,*) 'elev(k)= ',elev(k),' para k= ',k
   DO j=1,nlev(k)
     hgt3d(j,k)=real(zp(j,k))
   ENDDO
   write(*,*) 'PASSEI AQUI 3, com k= ',k
!
!-----------------------------------------------------------------------
!
!  Find mass weighted mean wind in first 500m
!
!-----------------------------------------------------------------------
!
     sumu=0.
     sumv=0.
     sump=0.
     h500m=elev(k)+500.
     DO j=2,nlev(k)
       IF( hgt3d(j,k) < h500m ) THEN
         dp=p_pa(j-1,k)-p_pa(j,k)
         rhohi=p_pa(j,k)/t3d(j,k)
         rholo=p_pa(j-1,k)/t3d(j-1,k)
         rhoinv=1./(rhohi+rholo)
         sumu=sumu+dp*rhoinv*(rholo*u3d(j-1,k)+rhohi*u3d(j,k))
         sumv=sumv+dp*rhoinv*(rholo*v3d(j-1,k)+rhohi*v3d(j,k))
         sump=sump+dp
       ELSE
         dz=hgt3d(j,k)-hgt3d(j-1,k)
         wlow=(hgt3d(j,k)-h500m)/dz
         u500m=u3d(j,k)*(1.-wlow) + u3d(j-1,k)*wlow
         v500m=v3d(j,k)*(1.-wlow) + v3d(j-1,k)*wlow
         p500m=p_pa(j,k)*(1.-wlow) + p_pa(j-1,k)*wlow
         t500m=t3d(j,k)*(1.-wlow) + t3d(j-1,k)*wlow
         dp=p_pa(j-1,k)-p500m
         rhohi=p500m/t500m
         rholo=p_pa(j-1,k)/t3d(j-1,k)
         rhoinv=1./(rhohi+rholo)
         sumu=sumu+dp*rhoinv*(rholo*u3d(j-1,k)+rhohi*u500m)
         sumv=sumv+dp*rhoinv*(rholo*v3d(j-1,k)+rhohi*v500m)
         sump=sump+dp
!          print *, ' sumu,sumv,sump = ',sumu,sumv,sump,'para k= ',k
         EXIT
       END IF
     END DO
!      121    CONTINUE
     u500m=sumu/sump
     v500m=sumv/sump
     write(*,*) 'PASSEI AQUI 4, com k= ',k
!      write(*,*) '(u500m,v500m)= ',u500m,v500m,' para k= ',k
!
!-----------------------------------------------------------------------
!
!  Find mass weighted mean wind sfc-2km AGL
!
!-----------------------------------------------------------------------
!
     sumu=0.
     sumv=0.
     sump=0.
     h2km=elev(k)+2000.
     DO j=2,nlev(k)
       IF( hgt3d(j,k) < h2km ) THEN
         dp=p_pa(j-1,k)-p_pa(j,k)
         rhohi=p_pa(j,k)/t3d(j,k)
         rholo=p_pa(j-1,k)/t3d(j-1,k)
         rhoinv=1./(rhohi+rholo)
         sumu=sumu+dp*rhoinv*(rholo*u3d(j-1,k)+rhohi*u3d(j,k))
         sumv=sumv+dp*rhoinv*(rholo*v3d(j-1,k)+rhohi*v3d(j,k))
         sump=sump+dp
       ELSE
         dz=hgt3d(j,k)-hgt3d(j-1,k)
         wlow=(hgt3d(j,k)-h2km)/dz
         u2km=u3d(j,k)*(1.-wlow) + u3d(j-1,k)*wlow
         v2km=v3d(j,k)*(1.-wlow) + v3d(j-1,k)*wlow
         p2km=p_pa(j,k)*(1.-wlow) + p_pa(j-1,k)*wlow
         t2km=t3d(j,k)*(1.-wlow) + t3d(j-1,k)*wlow
         dp=p_pa(j-1,k)-p2km
         rhohi=p2km/t2km
         rholo=p_pa(j-1,k)/t3d(j-1,k)
         rhoinv=1./(rhohi+rholo)
         sumu=sumu+dp*rhoinv*(rholo*u3d(j-1,k)+rhohi*u2km)
         sumv=sumv+dp*rhoinv*(rholo*v3d(j-1,k)+rhohi*v2km)
         sump=sump+dp
         EXIT
       END IF
     END DO
!      141    CONTINUE
     write(*,*) 'PASSEI AQUI 4.5, com k= ',k
     u2km=sumu/sump
     v2km=sumv/sump
     write(*,*) 'u2km,v2km: ',u2km,v2km
     tem2(k)=u2km
     tem3(k)=v2km
     write(*,*) 'tem2(k),tem3(k): ',tem2(k),tem3(k)
     write(*,*) 'PASSEI AQUI 5, com k= ',k
!
!-----------------------------------------------------------------------
!
!  Find mass weighted mean wind sfc-6km AGL
!
!-----------------------------------------------------------------------
!
     sumu=0.
     sumv=0.
     sump=0.
     h6km=elev(k)+6000.
     DO j=2,nlev(k)
       IF( hgt3d(j,k) < h6km ) THEN
         dp=p_pa(j-1,k)-p_pa(j,k)
         rhohi=p_pa(j,k)/t3d(j,k)
         rholo=p_pa(j-1,k)/t3d(j-1,k)
         rhoinv=1./(rhohi+rholo)
         sumu=sumu+dp*rhoinv*(rholo*u3d(j-1,k)+rhohi*u3d(j,k))
         sumv=sumv+dp*rhoinv*(rholo*v3d(j-1,k)+rhohi*v3d(j,k))
         sump=sump+dp
       ELSE
         dz=hgt3d(j,k)-hgt3d(j-1,k)
         wlow=(hgt3d(j,k)-h6km)/dz
         u6km=u3d(j,k)*(1.-wlow) + u3d(j-1,k)*wlow
         v6km=v3d(j,k)*(1.-wlow) + v3d(j-1,k)*wlow
         p6km=p_pa(j,k)*(1.-wlow) + p_pa(j-1,k)*wlow
         t6km=t3d(j,k)*(1.-wlow) + t3d(j-1,k)*wlow
         dp=p_pa(j-1,k)-p6km
         rhohi=p6km/t6km
         rholo=p_pa(j-1,k)/t3d(j-1,k)
         rhoinv=1./(rhohi+rholo)
         sumu=sumu+dp*rhoinv*(rholo*u3d(j-1,k)+rhohi*u6km)
         sumv=sumv+dp*rhoinv*(rholo*v3d(j-1,k)+rhohi*v6km)
         sump=sump+dp
         EXIT
       END IF
     END DO
     u6km=sumu/sump
     v6km=sumv/sump
     write(*,*) 'PASSEI AQUI 6, com k= ',k,' e (u6km,v6km)= ',u6km,v6km
!
!-----------------------------------------------------------------------
!
!  Find mass weighted mean wind 2km-9km AGL
!
!-----------------------------------------------------------------------
!
     sumu=0.
     sumv=0.
     sump=0.
     h9km=elev(k)+9000.
     DO j=2,nlev(k)
       IF( hgt3d(j,k) > h2km ) EXIT
     END DO
!      181   CONTINUE
     k2km=j
     dz=hgt3d(j,k)-hgt3d(j-1,k)
     wlow=(hgt3d(j,k)-h2km)/dz
     u2=u3d(j,k)*(1.-wlow) + u3d(j-1,k)*wlow
     v2=v3d(j,k)*(1.-wlow) + v3d(j-1,k)*wlow
     p2km=p_pa(j,k)*(1.-wlow) + p_pa(j-1,k)*wlow
     t2km=t3d(j,k)*(1.-wlow) + t3d(j-1,k)*wlow
     dp=p2km-p_pa(j,k)
     rholo=p2km/t2km
     rhohi=p_pa(j,k)/t3d(j,k)
     rhoinv=1./(rhohi+rholo)
     sumu=sumu+dp*rhoinv*(rholo*u2+rhohi*u3d(j,k))
     sumv=sumv+dp*rhoinv*(rholo*v2+rhohi*v3d(j,k))
     sump=sump+dp
     DO j=k2km+1,nlev(k)
       IF( hgt3d(j,k) < h9km ) THEN
         dp=p_pa(j-1,k)-p_pa(j,k)
         rhohi=p_pa(j,k)/t3d(j,k)
         rholo=p_pa(j-1,k)/t3d(j-1,k)
         rhoinv=1./(rhohi+rholo)
         sumu=sumu+dp*rhoinv*(rholo*u3d(j-1,k)+rhohi*u3d(j,k))
         sumv=sumv+dp*rhoinv*(rholo*v3d(j-1,k)+rhohi*v3d(j,k))
         sump=sump+dp
       ELSE
         dz=hgt3d(j,k)-hgt3d(j-1,k)
         wlow=(hgt3d(j,k)-h9km)/dz
         u9km=u3d(j,k)*(1.-wlow) + u3d(j-1,k)*wlow
         v9km=v3d(j,k)*(1.-wlow) + v3d(j-1,k)*wlow
         p9km=p_pa(j,k)*(1.-wlow) + p_pa(j-1,k)*wlow
         t9km=t3d(j,k)*(1.-wlow) + t3d(j-1,k)*wlow
         dp=p_pa(j-1,k)-p9km
         rhohi=p9km/t9km
         rholo=p_pa(j-1,k)/t3d(j-1,k)
         rhoinv=1./(rhohi+rholo)
         sumu=sumu+dp*rhoinv*(rholo*u3d(j-1,k)+rhohi*u9km)
         sumv=sumv+dp*rhoinv*(rholo*v3d(j-1,k)+rhohi*v9km)
         sump=sump+dp
         EXIT
       END IF
     END DO
!      191    CONTINUE
     u9km=sumu/sump
     v9km=sumv/sump
     write(*,*) 'PASSEI AQUI 7, com k= ',k
!
!-----------------------------------------------------------------------
!
!  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.
!
!  (04/19/2004) Modified by Ernani L. Nascimento. Adapting for the
!  Southern Hemisphere.
!
!-----------------------------------------------------------------------
!
     CALL uv2ddff(u6km,v6km,dirmean,spmean)
!      write(*,*) 'PASSAGEM 5.1, com k= ',k
!      write(*,*) 'PASSAGEM 5.1, com k= ',k,' e spmean= ',spmean
!      write(*,*) 'PASSAGEM 5.1, com k= ',k,' e dirmean= ',dirmean
     IF(spmean >= 20.0) THEN
!        write(*,*) 'Hei, passei aqui!'
       dirmean=dirmean-18.
       IF(dirmean <= 0.) dirmean=360.+dirmean
       spmean=spmean*0.89
     ELSE IF (spmean > 8.0) THEN
!        write(*,*) 'Não! Passei foi aqui!'
       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 <= 0.) dirmean=360.+dirmean
       spmean=spmean*perc
     ELSE
!        write(*,*) 'Na verdade, passei foi aqui!'
       dirmean=dirmean-32.
       IF(dirmean <= 0.) dirmean=360.+dirmean
       spmean=spmean*0.75
     END IF
!      write(*,*) 'PASSAGEM 5.2, com k= ',k
     arg = (dirmean * (3.141592654/180.))
     ustrm(k) = -spmean * sin(arg)
     vstrm(k) = -spmean * cos(arg)
! Utilizando movimento estimado da célula da esquerda em 9 de outubro de 2003
     IF (strm_opt == 1) then
      IF ((estacao(k)==83827).and.(hora(k)==00).and.(dia(k)==09).and. &
          (mes(k)==10).and.(ano(k)==2003)) THEN
       obs_ddstrm=270.0
       obs_ffstrm=10.7
       CALL ddff2uv(obs_ddstrm,obs_ffstrm,new_ustrm,new_vstrm)
       ustrm(k)=new_ustrm
       vstrm(k)=new_vstrm
      ENDIF
! Utilizando movimento estimado da célula da célula de 24 de maio de 2005
!       IF ((estacao(k)==99999).and.(hora(k)==18).and.(dia(k)==24).and. &
!           (mes(k)==05).and.(ano(k)==2005)) THEN
!        obs_ddstrm=297.0
!        obs_ffstrm=17.8
!        CALL ddff2uv(obs_ddstrm,obs_ffstrm,new_ustrm,new_vstrm)
!        ustrm(k)=new_ustrm
!        vstrm(k)=new_vstrm
!       ENDIF
     ENDIF
!      CALL ddff2uv(dirmean,spmean,ustrm(k),vstrm(k),5)
!      write(*,*) 'PASSAGEM 6, com k= ',k
!
!-----------------------------------------------------------------------
!
!  Storm-relative low-level flow
!
!-----------------------------------------------------------------------
!
     srlfl(k)=SQRT((ustrm(k)-u2km)*(ustrm(k)-u2km) +             &
                     (vstrm(k)-v2km)*(vstrm(k)-v2km))
!
!-----------------------------------------------------------------------
!
!  Storm relative mid-level flow
!
!-----------------------------------------------------------------------
!
     srmfl(k)=SQRT((ustrm(k)-u9km)*(ustrm(k)-u9km) +             &
                     (vstrm(k)-v9km)*(vstrm(k)-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(k)=0.5*( (u6km-u500m)*(u6km-u500m) +                       &
                     (v6km-v500m)*(v6km-v500m) )
! Added the 2km BRNSHR (Ernani L. Nascimento, 03/feb/2005
     brnu2km(k)=0.5*( (u2km-u500m)*(u2km-u500m) +                    &
                     (v2km-v500m)*(v2km-v500m) )
!     write(*,*) 'PASSAGEM 7, com k= ',k
!
!-----------------------------------------------------------------------
!
!  Bulk Richardson number
!  A limit of 200 is imposed, since this could
!  go to inifinity.
!
!-----------------------------------------------------------------------
!
     IF(brnu(k) > 0.) THEN
       brn(k)=cape(k)/brnu(k)
       brn(k)=AMIN1(brn(k),200.)
     ELSE
       brn(k)=200.
     END IF
    write(*,*) 'PASSEI AQUI 8, com k= ',k
!
!
!-----------------------------------------------------------------------
!
!  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.
!
!-----------------------------------------------------------------------
!
     h3km=elev(k)+3000.
     h7km=elev(k)+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 j=2,nlev(k)
       IF(hgt3d(j,k) > h3km) EXIT
       sumh=sumh +                                                     &
           ( u3d(j,k)*v3d(j-1,k) ) -                                   &
           ( v3d(j,k)*u3d(j-1,k) )
     END DO
!      240    CONTINUE
     k3km=j
     dz=hgt3d(j,k)-hgt3d(j-1,k)
     wlow=(hgt3d(j,k)-h3km)/dz
     u3km=u3d(j,k)*(1.-wlow) + u3d(j-1,k)*wlow
     v3km=v3d(j,k)*(1.-wlow) + v3d(j-1,k)*wlow
     sumh=sumh +                                                       &
           ( u3km*v3d(j-1,k) ) -                                       &
           ( v3km*u3d(j-1,k) )
     ushr=u3km-u3d(1,k)
     vshr=v3km-v3d(1,k)
     helicity(k)=sumh + vshr*ustrm(k) - ushr*vstrm(k)
!      write(*,*) 'PASSAGEM 9, com k= ',k
!
!-----------------------------------------------------------------------
!
!  Now Find 7km wind for 3-to-7km shear
!
!-----------------------------------------------------------------------
!
     DO j=k3km,nlev(k)
       IF(hgt3d(j,k) > h7km) EXIT
     END DO
!      260   CONTINUE
     dz=hgt3d(j,k)-hgt3d(j-1,k)
     wlow=(hgt3d(j,k)-h7km)/dz
     u7km=u3d(j,k)*(1.-wlow) + u3d(j-1,k)*wlow
     v7km=v3d(j,k)*(1.-wlow) + v3d(j-1,k)*wlow
     shr37(k)=(SQRT( (u7km-u3km)*(u7km-u3km) +                       &
                       (v7km-v3km)*(v7km-v3km) ))/4000.
  endif
 END DO
!  write(*,*) 'Saída da subrotina calcshr.'

!-----------------------------------------------------------------------
!
!  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
!
!
!
!     ##################################################################
!     ##################################################################
!     ######                                                      ######
!     ######                SUBROUTINE CALCEHI_SUP                ######
!     ######                                                      ######
!     ######                     Developed by                     ######
!     ######     Center for Analysis and Prediction of Storms     ######
!     ######                University of Oklahoma                ######
!     ######        and by Ernani L. Nascimento at SIMEPAR        ######
!     ######                                                      ######
!     ##################################################################
!     ##################################################################
!
SUBROUTINE calcehi_sup(nrad,nlev,presshpa,cape,heli,dnrv,ehi,sup)
!
!#######################################################################
!
!     PURPOSE:
!
!     Calculate the energy helicity-index (EHI) (DAVIES 1993; RASMUSSEN AND
!     BLANCHARD 1998), and the supercell composite index (THOMPSON et al
!     2003).
!
!#######################################################################
!
!     AUTHOR: Ernani L. Nascimento, adapting code from arpspltderive.f
!     09 February, 2004
!
!     MODIFICATION HISTORY:
!     THIS IS A NEW SUBROUTINE, NOT AVAILABLE IN THE ORIGINAL ARPS CODE
!
!     (04/19/2004) (Ernani L. Nascimento)
!     Code was modified to fit in índices_severos3.f90
!
!     (08/02/2004) (Ernani L. Nascimento)
!     Code was modified to include the computation of the supercell
!     composite index.
!
!#######################################################################
!
!     Calculates EHI and SUP from cape, helicity and bulk Richardson
!     shear.
!
!
!     nrad          Number of soundings for which EHI will be computed
!     cape          Convective available potential energy (J/kg)
!     heli          Helicity, storm relative (m²/s²)
!     dnrv          Bulk Richardson number shear (m²/s²)
!     ehi           Energy-helicity index (non-dimensional)
!     sup           Supercell composite index
!
!#######################################################################
!
   implicit none
!
!#######################################################################
!
!     Input variables
!
!#######################################################################
!
   integer :: nrad,k,j
   integer :: nlev(nrad)
   real :: cape(nrad), heli(nrad), dnrv(nrad)
   real :: presshpa(102,nrad)
!
!#######################################################################
!
!     Output variables
!
!#######################################################################
!
   real :: ehi(nrad)      ! Energy-helicity index
   real :: sup(nrad)      ! Supercell composite index

!#######################################################################
!     Work variable (constant)
!#######################################################################

   real :: denomi

!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!     Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
   denomi=160000.0
   DO  k=1,nrad
    if (nlev(k)==0) then
     ehi(k)=-999.9
     sup(k)=-999.9
    else
     j=nlev(k)
     if (presshpa(j,k)>300.0) then
      cycle
     endif
     ehi(k) = ( cape(k)*heli(k) )/denomi
     sup(k) = ( (cape(k)/1000.)*(heli(k)/150.)*(dnrv(k)/40.) )
    endif
   ENDDO

 RETURN
END SUBROUTINE calcehi_sup
