c -----------------------------------------------------------------------
c ------- DJS created the subroutine interface and most of documentation
c -------
c -------    *** NOT COMPLETE or TESTED ****************

C NCLFORTSTART
      subroutine cape (p,tt,q,nlvl,xmsg,iprflg       ! input
     *                ,cape_pos,cape_neg,dcape,ier)  ! output

c Given standard data from a raob sounding calculate
c .   various properties of samples of air raised or lowered
c .   to different levels.

c The original code was provided by Natalie Mahowald.
c .   She got it from Kerry Emmanual in 1993.
c .   Loop 250 has some problems for soundings extending
c .        into the stratosphere.

c No provision is made for missing data. Check BEFORE CALLING
c .   to make sure no missing data is present.

c Input Nomenclature:
c .   p        - pressure (mb or hPa)
c .   tt       - temperature (C)
c .   q        - mixing ratio (g/kg)
c .   xmsg     - code for undefined quantities (eg, 1.e+36)
c .   nlvl     - number of sounding levels
c .   iprflg   - print flag: =0 for no print; otherwise print results
c Output Nomenclature:
c .   cape_pos - integrated positive CAPE (joules/kg)
c .   cape_neg - 
c .   dcape    - down draft CAPE
c .   ier      - error code (not really used)

      implicit none
      integer nlvl, iprflg, ier                 ! input
      real    p(nlvl), tt(nlvl), q(nlvl), xmsg  ! input
      real    cape_pos, cape_neg, dcape(nlvl)   ! output
C NCLEND

c local arrays: these contain various quantities 

      integer     nlvmax
      parameter  (nlvmax=200)      
      real r     (nlvmax)          ! dimensionless mix ratio
     *   , t     (nlvmax)          ! temperature (K)
     *   , ev    (nlvmax)          ! vapor pressure
     *   , es    (nlvmax)          ! saturation vapor pressure
     *   , tvd   (nlvmax)
     *   , caper (nlvmax)          ! cape reversible
     *   , capep (nlvmax)          ! cape pseudo adiabat
     *   , par   (nlvmax)          ! positive area reversible
     *   , pap   (nlvmax)          ! positive area pseudo-adiabat
     *   , nar   (nlvmax)          ! negative area reversible
     *   , nap   (nlvmax)          ! negative area pseudo-adiabat
ccccc*   , dcape (nlvmax)          ! downdraft cape (pass in)

c note the double dimensions

      real lw    (nlvmax,nlvmax)
     *   , tlr   (nlvmax,nlvmax)
     *   , tlp   (nlvmax,nlvmax)
     *   , tlvr  (nlvmax,nlvmax)
     *   , tlvp  (nlvmax,nlvmax)
     *   , tvrdif(nlvmax,nlvmax)
     *   , tvpdif(nlvmax,nlvmax)

c constants & temporaries

      real sp, spd, spg, sum1, sum2, tc, tg, tg0, tgd0, tcdifm, tvm   
     *   , ah, ahg, alv, alv0, alv1, chi, cl, cpd, cpv, cpvmcl
     *   , cpw, eg, em, enew, eps, rd, pm, rs, rv, s, sg, sl
     *   , slope, slp, tvdifm, plcl, rg, rg0, rgd0, rh 

      integer i, j, k, n, nl, ml, icb, imax, inbp, inbr

      ier = 0
      cape_pos = xmsg

c initilize all local arrays to xmsg

      call cape_set (r,t,ev,es,tvd,caper,capep,dcape
     *              ,par,pap,nar,nap,lw,tlr,tlp
     *              ,tlvr,tlvp,tvrdif,tvpdif,nlvmax,xmsg)

      if (nlvl.gt.nlvmax) then
          ier = 1
          write (*,'(/,'' SUB CAPE: nlvl, nlvmax=''
     *                ,2i5,'': increase nlvmax'',/)') nlvl, nlvmax 
          return 
      endif

c thermodynamic  and other constants 

      cl     = 4190.0     ! specific heat of liquid water ??? [=4187 J/(kg-deg)]
      rv     = 461.5      ! gas constant for water vapor 
      rd     = 287.04     ! gas const dry air (J/{kg-K})
      cpd    = 1005.7     ! specific heat for dry air at 288K
                          ! =1004.6 for dry air at 273K
      cpv    = 1870.0
      eps    = rd/rv
      alv0   = 2.501e6
      cpvmcl = 2320.0

c change units + calculate some basic quantities 

      do nl=1,nlvl
         r (nl) = q(nl)*0.001
         ev(nl) = r(nl)*p(nl)/(eps+r(nl))
         es(nl) = 6.112*exp(17.67*tt(nl)/(243.5+tt(nl)))
         t (nl) = tt(nl)+273.15
      enddo

      if (iprflg.ne.0) then
          write (*,'(/,25(''>''),1x,'' SUB CAPE: SOUNDING INFO: '',/)') 
          do nl=1,nlvl
             write (*,'(i5,f7.1,2f7.1,'' <==> '',f7.1,f11.6,2f10.4)') 
     *              nl,p(nl),tt(nl),q(nl),t(nl),r(nl),ev(nl),es(nl) 
          enddo
      endif

      n = nlvl 
c
c   ***  begin outer loop, which cycles through parcel origin levels i ***
c  
	do 500 i=1,n
c
c   ***  define various conserved parcel quantities: reversible   ***
c   ***        entropy, s, pseudo-adiabatic entropy, sp,          *** 
c   ***                   and enthalpy, ah                        ***
c
	rs =eps*es(i)/(p(i)-es(i))
	alv=alv0-cpvmcl*(t(i)-273.15)
	em =max(ev(i),1.0e-6) 
	s  =(cpd+r(i)*cl)*log(t(i))-rd*log(p(i)-ev(i))+
     1       alv*r(i)/t(i)-r(i)*rv*log(em/es(i))
	sp =cpd*log(t(i))-rd*log(p(i)-ev(i))+
     1       alv*r(i)/t(i)-r(i)*rv*log(em/es(i))
	ah =(cpd+r(i)*cl)*t(i)+alv*r(i)
c      
c   ***  find the temperature and mixing ratio of the parcel at   ***
c   ***    level i saturated by a wet bulb process                ***
c
	slope=cpd+alv*alv*rs/(rv*t(i)*t(i))
	tg   =t(i)
	rg   =rs  
	do 100 j=1,20 
	 alv1=alv0-cpvmcl*(tg-273.15)
	 ahg =(cpd+cl*rg)*tg+alv1*rg
	 tg  =tg+(ah-ahg)/slope
	 tc  =tg-273.15
	 enew=6.112*exp(17.67*tc/(243.5+tc))
	 rg  =eps*enew/(p(i)-enew)
  100	continue
c   
c   ***  calculate conserved variable at top of downdraft   ***
c
	eg =rg*p(i)/(eps+rg)
	spd=cpd*log(tg)-rd*log(p(i)-eg)+alv1*rg/tg
	tvd(i)=tg*(1.+rg/eps)/(1.+rg)-t(i)*(1.+r(i)/eps)/
     1        (1.+r(i))
	if(p(i).lt.100.0)tvd(i)=0.0
	rgd0=rg
	tgd0=tg
c
c   ***   find lifted condensation pressure     ***
c
	rh  =r(i)/rs
	rh  =min(rh,1.0)
	chi =t(i)/(1669.0-122.0*rh-t(i))
	plcl=1.0
	if(rh.gt.0.0)then
	 plcl=p(i)*(rh**chi)
	end if
c
c   ***  begin updraft loop   ***
c
	sum1=0.0
	rg0=r(i)
	tg0=t(i)
	do 200 j=i,n
c
c   ***  calculate estimates of the rates of change of the entropies  ***
c   ***           with temperature at constant pressure               ***
c  
	 rs =eps*es(j)/(p(j)-es(j))
	 alv=alv0-cpvmcl*(t(j)-273.15)
	 sl =(cpd+r(i)*cl+alv*alv*rs/(rv*t(j)*t(j)))/t(j)
	 slp=(cpd+rs*cl+alv*alv*rs/(rv*t(j)*t(j)))/t(j)
c   
c   ***  calculate lifted parcel temperature below its lcl   ***
c
	 if(p(j).ge.plcl)then
	  tlr(i,j)=t(i)*(p(j)/p(i))**(rd/cpd)
	  tlp(i,j)=tlr(i,j) 
	  lw(i,j)=0.0
	  tlvr(i,j)=tlr(i,j)*(1.+r(i)/eps)/(1.+r(i))
	  tlvp(i,j)=tlvr(i,j)
	  tvrdif(i,j)=tlvr(i,j)-t(j)*(1.+r(j)/eps)/(1.+r(j))
	  tvpdif(i,j)=tvrdif(i,j)
	 else
c
c   ***  iteratively calculate lifted parcel temperature and mixing   ***
c   ***    ratios for both reversible and pseudo-adiabatic ascent     ***
c
	 tg=t(j)
	 rg=rs
	 do 150 k=1,20
	  em=rg*p(j)/(eps+rg)
	  alv=alv0-cpvmcl*(tg-273.15)
	  sg=(cpd+r(i)*cl)*log(tg)-rd*log(p(j)-em)+alv*rg/tg
	  tg=tg+(s-sg)/sl  
	  tc=tg-273.15
	  enew=6.112*exp(17.67*tc/(243.5+tc))
	  rg=eps*enew/(p(j)-enew)           
  150	 continue
	 tlr(i,j)=tg
	 tlvr(i,j)=tg*(1.+rg/eps)/(1.+r(i))
	 lw(i,j)=r(i)-rg
	 lw(i,j)=max(0.0,lw(i,j))
	 tvrdif(i,j)=tlvr(i,j)-t(j)*(1.+r(j)/eps)/(1.+r(j))
c
c   ***   now do pseudo-adiabatic ascent   ***
c
	 tg=t(j)
	 rg=rs
	 do 180 k=1,20 
	  cpw=0.0
	  if(j.gt.1)then
	   cpw=sum1+cl*0.5*(rg0+rg)*(log(tg)-log(tg0))
	  end if
	  em=rg*p(j)/(eps+rg)
	  alv=alv0-cpvmcl*(tg-273.15)
	  spg=cpd*log(tg)-rd*log(p(j)-em)+cpw+
     1      alv*rg/tg
	  tg=tg+(sp-spg)/slp  
	  tc=tg-273.15
	  enew=6.112*exp(17.67*tc/(243.5+tc))
	  rg=eps*enew/(p(j)-enew)           
  180	 continue
	 tlp(i,j)=tg
	 tlvp(i,j)=tg*(1.+rg/eps)/(1.+rg)
	 tvpdif(i,j)=tlvp(i,j)-t(j)*(1.+r(j)/eps)/(1.+r(j))
	 rg0=rg
	 tg0=tg
	 sum1=cpw
 	end if
  200	continue
	if(i.eq.1)goto 500
c
c   ***  begin downdraft loop   ***
c
	sum2=0.0
	do 300 j=i-1,1,-1
c
c   ***  calculate estimate of the rate of change of entropy          ***
c   ***           with temperature at constant pressure               ***
c  
	 rs=eps*es(j)/(p(j)-es(j))
	 alv=alv0-cpvmcl*(t(j)-273.15)
	 slp=(cpd+rs*cl+alv*alv*rs/(rv*t(j)*t(j)))/t(j)
	 tg=t(j)
	 rg=rs
c
c   ***  do iteration to find downdraft temperature   ***
c
c >>>>>>>>>>>>>>>> Dan Davidoff [MIT] says the following
c >>>>>>>>>>>>>>>> loop which causes problems when the sounding
c >>>>>>>>>>>>>>>> extends to the stratosphere may be skipped.
c >>>>>>>>>>>>>>>> It does not affect CAPE computations.
C C C    GO TO 500 ! D Davidoff says this will not hurt cape computation

	 do 250 k=1,20
          IF (TG.LE.0. .OR. TGD0.LE.0.) THEN      ! djs added
c c c         WRITE (*,'('' SUB CAPE: LOOP 250: K,I,J,P(J),TG,TGD0=''
c c c*                  ,3i5,f8.1,2F10.3)') K,I,J,P(J),TG,TGD0 
              GO TO 300
          ENDIF
	  cpw=sum2+cl*0.5*(rgd0+rg)*(log(tg)-log(tgd0))
	  em =rg*p(j)/(eps+rg)
	  alv=alv0-cpvmcl*(tg-273.15)
          IF ((P(J)-EM).LE.0.) THEN            ! djs added
c             WRITE (*,'('' SUB CAPE: LOOP 250: K,I,J,P(J),EM=''
c    *                  ,3i5,f8.1,2F10.3)') K,I,J,P(J),EM
              GO TO 300
          ENDIF
	  spg =cpd*log(tg)-rd*log(p(j)-em)+cpw+alv*rg/tg
	  tg  =tg+(spd-spg)/slp  
	  tc  =tg-273.15
	  enew=6.112*exp(17.67*tc/(243.5+tc))
	  rg  =eps*enew/(p(j)-enew)           
  250	 continue
	 sum2=cpw
	 tgd0=tg
	 rgd0=rg
	 tlp(i,j)=tg
	 tlvp(i,j)=tg*(1.+rg/eps)/(1.+rg)
	 tvpdif(i,j)=tlvp(i,j)-t(j)*(1.+r(j)/eps)/(1.+r(j))
	 if(p(i).lt.100.0)tvpdif(i,j)=0.0
	 tvpdif(i,j)=min(tvpdif(i,j),0.0)
	 tlr(i,j)=t(j)
	 tlvr(i,j)=t(j)
	 tvrdif(i,j)=0.0
	 lw(i,j)=0.0
  300	continue
  500   continue
c
c  ***  begin loop to find cape, pa, and na from reversible and ***
c  ***            pseudo-adiabatic ascent, and dcape            ***
c
	do 800 i=1,n
	 caper(i)=0.0
	 capep(i)=0.0
	 dcape(i)=0.0
	 pap(i)=0.0
	 par(i)=0.0
	 nap(i)=0.0
	 nar(i)=0.0
c
c   ***   find lifted condensation pressure     ***
c
	rs=eps*es(i)/(p(i)-es(i))
	rh=r(i)/rs
	rh=min(rh,1.0)
	chi=t(i)/(1669.0-122.0*rh-t(i))
	plcl=1.0
	if(rh.gt.0.0)then
	 plcl=p(i)*(rh**chi)
	end if
c
c   ***  find lifted condensation level and maximum level   ***
c   ***               of positive buoyancy                  ***
c
	 icb=n
	 inbr=1
	 inbp=1
	 do 550 j=n,i,-1
	  if(p(j).lt.plcl)icb=min(icb,j)
	  if(tvrdif(i,j).gt.0.0)inbr=max(inbr,j)
	  if(tvpdif(i,j).gt.0.0)inbp=max(inbp,j)
  550	 continue
	  imax=max(inbr,i)
	   do 555 j=imax,n
	    tvrdif(i,j)=0.0
  555	   continue
	  imax=max(inbp,i)
	   do 565 j=imax,n
	    tvpdif(i,j)=0.0
  565	   continue
c
c   ***  do updraft loops        ***
c
	 if(inbr.gt.i)then
	  do 600 j=i+1,inbr
	   tvm=0.5*(tvrdif(i,j)+tvrdif(i,j-1))
	   pm=0.5*(p(j)+p(j-1))
	   if(tvm.le.0.0)then
	    nar(i)=nar(i)-rd*tvm*(p(j-1)-p(j))/pm
	   else
	    par(i)=par(i)+rd*tvm*(p(j-1)-p(j))/pm
	   end if
  600	  continue
	  caper(i)=par(i)-nar(i)
	 end if
	 if(inbp.gt.i)then
	  do 650 j=i+1,inbp
	   tvm=0.5*(tvpdif(i,j)+tvpdif(i,j-1))
	   pm=0.5*(p(j)+p(j-1))
	   if(tvm.le.0.0)then
	    nap(i)=nap(i)-rd*tvm*(p(j-1)-p(j))/pm
	   else
	    pap(i)=pap(i)+rd*tvm*(p(j-1)-p(j))/pm
	   end if
  650	  continue
	  capep(i)=pap(i)-nap(i)
	 end if
c
c  ***       find dcape     ***
c
	 if(i.eq.1)goto 800
	 do 700 j=i-1,1,-1
	  tvdifm=tvpdif(i,j+1)
	  if(i.eq.(j+1))tvdifm=tvd(i)
	  tvm=0.5*(tvpdif(i,j)+tvdifm)
	  pm=0.5*(p(j)+p(j+1))
	  if(tvm.lt.0.0)then
	   dcape(i)=dcape(i)-rd*tvm*(p(j)-p(j+1))/pm
	  end if
  700	 continue
  800	continue

c
c  ***  write values of pa, na, cape and dcape   ***
c
      if (iprflg.ne.0) then
          write(*,'(/,1x,''SUB CAPE: all areas in units of j/kg'')')
          write(*,'(1x,''origin'',4x,''rev.'',4x,''p.a.''
     *             ,4x,''rev.'',4x,''p.a.'',4x,''rev.'',4x,''p.a.'')')
          write(*,'(1x,''p (mb)'',4x,'' pa '',4x,'' pa ''
     *             ,4x,'' na '',4x,'' na '',4x,''cape''
     *             ,4x,''cape'',4x,''dcape'')')
          write(*,'(1x,''------'',4x,''----'',4x,''----'',4x,''----''
     *             ,4x,''----'',4x,''----'',4x,''----'',4x,''-----'')')
          do nl=1,nlvl
             write(*,'(1x,f6.1,7(f8.1))')
     *                 p(nl),par(nl),pap(nl),nar(nl),nap(nl)
     *              ,  caper(nl),capep(nl),dcape(nl)
          enddo
      endif
     

      return
      end
c -----------------------------------------------------------------------
      subroutine cape_set (r,t,ev,es,tvd,caper,capep,dcape
     *                    ,par,pap,nar,nap,lw,tlr,tlp
     *                    ,tlvr,tlvp,tvrdif,tvpdif,nlvmax,value)
      implicit none

c set all array elements to a particular value.

      integer     nlvmax
      real r     (nlvmax)
     *   , t     (nlvmax)
     *   , ev    (nlvmax)
     *   , es    (nlvmax)
     *   , tvd   (nlvmax)
     *   , caper (nlvmax)
     *   , capep (nlvmax)
     *   , dcape (nlvmax)
     *   , par   (nlvmax)
     *   , pap   (nlvmax)
     *   , nar   (nlvmax)
     *   , nap   (nlvmax)
      real lw    (nlvmax,nlvmax)
     *   , tlr   (nlvmax,nlvmax)
     *   , tlp   (nlvmax,nlvmax)
     *   , tlvr  (nlvmax,nlvmax)
     *   , tlvp  (nlvmax,nlvmax)
     *   , tvrdif(nlvmax,nlvmax)
     *   , tvpdif(nlvmax,nlvmax)
      real value                      ! value to which array elements
                                      ! are to be set/initialized
      integer nl, ml

      do nl=1,nlvmax
         r     (nl) = value
         t     (nl) = value
         ev    (nl) = value
         es    (nl) = value
         tvd   (nl) = value
         caper (nl) = value
         capep (nl) = value
         dcape (nl) = value
         par   (nl) = value
         pap   (nl) = value
         nar   (nl) = value
         nap   (nl) = value
      enddo

      do nl=1,nlvmax
        do ml=1,nlvmax
           lw    (ml,nl) = value
           tlr   (ml,nl) = value
           tlp   (ml,nl) = value
           tlvr  (ml,nl) = value
           tlvp  (ml,nl) = value
           tvrdif(ml,nl) = value
           tvpdif(ml,nl) = value
        enddo
      enddo

      return
      end