c *******************************************************************
c **
c **    these routines are from tom schlatter (profs)
c **                   14 may 1991
c **
c **    i have changed the names of several functions to prevent
c **    problems with my codes:
c *******************************************************************
c **
c **    i have added routines at the end
c **
c *******************************************************************
 
 
C NCLFORTSTART
	real function alclThermo(t,td,p)
c       NCL: aLcl = alclThermo(t:float, td:float, p:float)
c NCLEND
 
c	include 'lib_dev:[gudoc]edfvaxbox.for/list'
c	g.s. stipanuk     1973      	  original version.
c	reference stipanuk paper entitled:
c            "algorithms for generating a skew-t, log p
c	     diagram and computing selected meteorological
c	     quantities."
c	     atmospheric sciences laboratory
c	     u.s. army electronics command
c	     white sands missile range, new mexico 88002
c	     33 pages
c	baker, schlatter  17-may-1982	
 
c   this function returns the pressure alcl (mb) of the lifting conden-
c   sation level (lcl) for a parcel initially at temperature t (celsius)
c   dew point td (celsius) and pressure p (millibars). alcl is computed
c   by an iterative procedure described by eqs. 8-12 in stipanuk (1973),
c   pp.13-14.
c   determine the mixing ratio line through td and p.
 
	aw = wx(td,p)
 
c   determine the dry adiabat through t and p.
 
	ao = ox(t,p)         ! was ao = o(t,p)
 
c   iterate to locate pressure pi at the intersection of the two
c   curves. pi has been set to p for the initial guess.
 
 3	continue
	   pi = p
	   do 4 i= 1,10
	      x= .02*(tmr(aw,pi)-tda(ao,pi))
	      if (abs(x).lt.0.01) go to 5
 4	      pi= pi*(2.**(x))
 5	   alcl= pi
	return
 
	entry alclm(t,td,p)
c   for entry alclm only, t is the mean potential temperature (celsius)
c   and td is the mean mixing ratio (g/kg) of the layer containing the
c   parcel.
 
	aw = td
	ao = t
	go to 3
	end
 
c ---------------------------------------------------------------------
 
C NCLFORTSTART
	real function ctThermo(wbar,pc,ps)
c       NCL: ct = ctThermo(wbar:float, pc:float, ps:float)
c NCLEND
 
c	include 'lib_dev:[gudoc]edfvaxbox.for/list'
c	g.s. stipanuk     1973      	  original version.
c	reference stipanuk paper entitled:
c            "algorithms for generating a skew-t, log p
c	     diagram and computing selected meteorological
c	     quantities."
c	     atmospheric sciences laboratory
c	     u.s. army electronics command
c	     white sands missile range, new mexico 88002
c	     33 pages
c	baker, schlatter  17-may-1982	
 
c   this function returns the convective temperature ct (celsius)
c   given the mean mixing ratio wbar (g/kg) in the surface layer,
c   the pressure pc (mb) at the convective condensation level (ccl)
c   and the surface pressure ps (mb).
c   compute the temperature (celsius) at the ccl.
 
	tc= tmr(wbar,pc)
 
c   compute the potential temperature (celsius), i.e., the dry
c   adiabat ao through the ccl.
 
	ao= ox(tc,pc)          ! was ao = o(tc,pc)
 
c   compute the surface temperature on the same dry adiabat ao.
 
	ct= tda(ao,ps)
	return
	end
 
c ---------------------------------------------------------------------
 
C NCLFORTSTART
	real function dewpt1Thermo(ew)
C       NCL: dpT = dewpt1Thermo (ew:float)
c NCLEND
 
c   this function yields the dew point dewpt (celsius), given the
c   water vapor pressure ew (millibars).
c   the empirical formula appears in bolton, david, 1980:
c   "the computation of equivalent potential temperature,"
c   monthly weather review, vol. 108, no. 7 (july), p. 1047, eq.(11).
c   the quoted accuracy is 0.03c or less for -35 < dewpt < 35c.
 
c	include 'lib_dev:[gudoc]edfvaxbox.for/list'
c	baker, schlatter  17-may-1982	  original version.
 
	enl = alog(ew)
	dewpt = (243.5*enl-440.8)/(19.48-enl)
	return
	end
 
c ---------------------------------------------------------------------
 
C NCLFORTSTART
	real function dpt2Thermo(ew)
C       NCL: dpT = dewTew2Thermo (ew:float)
c NCLEND
 
c   this function returns the dew point dpt (celsius), given the
c   water vapor pressure ew (millibars).
 
c	include 'lib_dev:[gudoc]edfvaxbox.for/list'
c	baker, schlatter  17-may-1982	  original version.
 
	data es0/6.1078/
 
c   es0 = saturation vapor pressure (mb) over water at 0c
c   return a flag value if the vapor pressure is out of range.
 
	if (ew.gt..06.and.ew.lt.1013.) go to 5
	dpt = 9999.
	return
    5	continue
 
c   approximate dew point by means of teten's formula.
c   the formula appears as eq.(8) in bolton, david, 1980:
c   "the computation of equivalent potential temperature,"
c   monthly weather review, vol 108, no. 7 (july), p.1047.
c   the formula is ew(t) = es0*10**(7.5*t/(t+237.3))
c            or    ew(t) = es0*exp(17.269388*t/(t+237.3))
c   the inverse formula is used below.
 
	x = alog(ew/es0)
	dnm = 17.269388-x
	t = 237.3*x/dnm
	fac = 1./(ew*dnm)
 
c   loop for iterative improvement of the estimate of dew point
 
   10	continue
 
c   get the precise vapor pressure corresponding to t.
 
	edp = esw(t)
 
c   estimate the change in temperature corresponding to (ew-edp)
c   assume that the derivative of temperature with respect to
c   vapor pressure (dtdew) is given by the derivative of the
c   inverse teten formula.
 
	dtdew = (t+237.3)*fac
	dt = dtdew*(ew-edp)
	t = t+dt
	if (abs(dt).gt.1.e-04) go to 10
	dpt = t
	return
	end
 
c ---------------------------------------------------------------------
 
C NCLFORTSTART
	real function dwpt3Thermo(t,rh)
C       NCL: dpT = dewpt3Thermo (t:float, rh:float)
c NCLEND
 
c	include 'lib_dev:[gudoc]edfvaxbox.for/list'
c	baker, schlatter  17-may-1982	  original version.
 
c   this function returns the dew point (celsius) given the temperature
c   (celsius) and relative humidity (%). the formula is used in the
c   processing of u.s. rawinsonde data and is referenced in parry, h.
c   dean, 1969: "the semiautomatic computation of rawinsondes,"
c   technical memorandum wbtm edl 10, u.s. department of commerce,
c   environmental science services administration, weather bureau,
c   office of systems development, equipment development laboratory,
c   silver spring, md (october), page 9 and page ii-4, line 460.
 
	x = 1.-0.01*rh
 
c   compute dew point depression.
 
	dpd =(14.55+0.114*t)*x+((2.5+0.007*t)*x)**3+(15.9+0.117*t)*x**14
	dwpt = t-dpd
	return
	end
 
c ---------------------------------------------------------------------
 
C NCLFORTSTART
	real function eptThermo(t,td,p)
C       NCL: epT = eptThermo (t:float, td:float, p:float)
c NCLEND
 
c	include 'lib_dev:[gudoc]edfvaxbox.for/list'
c	baker, schlatter  17-may-1982	  original version.
 
c   this function returns the equivalent potential temperature ept
c   (celsius) for a parcel of air initially at temperature t (celsius),
c   dew point td (celsius) and pressure p (millibars). the formula used
c   is eq.(43) in bolton, david, 1980: "the computation of equivalent
c   potential temperature," monthly weather review, vol. 108, no. 7
c   (july), pp. 1046-1053. the maximum error in ept in 0.3c.  in most
c   cases the error is less than 0.1c.
c
c   compute the mixing ratio (grams of water vapor per kilogram of
c   dry air).
 
	w = wmr(p,td)
 
c   compute the temperature (celsius) at the lifting condensation level.
 
	tlcl = tcon(t,td)
	tk = t+273.15
	tl = tlcl+273.15
	pt = tk*(1000./p)**(0.2854*(1.-0.00028*w))
	eptk = pt*exp((3.376/tl-0.00254)*w*(1.+0.00081*w))
	ept= eptk-273.15
	return
	end
 
c ---------------------------------------------------------------------
 
C NCLFORTSTART
	real function esatThermo(t)
C       NCL: esat = esatThermo (t:float)
C NCLEND
 
c	include 'lib_dev:[gudoc]edfvaxbox.for/list'
c	g.s. stipanuk     1973      	  original version.
c	reference stipanuk paper entitled:
c            "algorithms for generating a skew-t, log p
c	     diagram and computing selected meteorological
c	     quantities."
c	     atmospheric sciences laboratory
c	     u.s. army electronics command
c	     white sands missile range, new mexico 88002
c	     33 pages
c	baker, schlatter  17-may-1982	
 
c   this function returns the saturation vapor pressure over
c   water (mb) given the temperature (celsius).
c   the algorithm is due to nordquist, w.s.,1973: "numerical approxima-
c   tions of selected meteorlolgical parameters for cloud physics prob-
c   lems," ecom-5475, atmospheric sciences laboratory, u.s. army
c   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
 
c ---------------------------------------------------------------------
 
C NCLFORTSTART
	real function esggThermo(t)
C       NCL: esgg = esggThermo (t:float)
c NCLEND
 
c	include 'lib_dev:[gudoc]edfvaxbox.for/list'
c	baker, schlatter  17-may-1982	  original version.
 
c   this function returns the saturation vapor pressure over liquid
c   water esgg (millibars) given the temperature t (celsius). the
c   formula used, due to goff and gratch, appears on p. 350 of the
c   smithsonian meteorological tables, sixth revised edition, 1963,
c   by roland list.
 
	data cta,ews,ts/273.15,1013.246,373.15/
 
c   cta = difference between kelvin and celsius temperatures
c   ews = saturation vapor pressure (mb) over liquid water at 100c
c   ts = boiling point of water (k)
 
	data c1,      c2,      c3,      c4,       c5,       c6
     1	/ 7.90298, 5.02808, 1.3816e-7, 11.344, 8.1328e-3, 3.49149 /
	tk = t+cta
 
c   goff-gratch formula
 
	rhs = -c1*(ts/tk-1.)+c2*alog10(ts/tk)-c3*(10.**(c4*(1.-tk/ts))
     1	      -1.)+c5*(10.**(-c6*(ts/tk-1.))-1.)+alog10(ews)
	esw = 10.**rhs
	if (esw.lt.0.) esw = 0.
	esgg = esw
	return
	end
 
c ---------------------------------------------------------------------
 
C NCLFORTSTART
	real function esiceThermo(t)
C       NCL: esice = esiceThermo (t:float)
c NCLEND
 
c	include 'lib_dev:[gudoc]edfvaxbox.for/list'
c	baker, schlatter  17-may-1982	  original version.
 
c   this function returns the saturation vapor pressure with respect to
c   ice esice (millibars) given the temperature t (celsius).
c   the formula used is based upon the integration of the clausius-
c   clapeyron equation by goff and gratch.  the formula appears on p.350
c   of the smithsonian meteorological tables, sixth revised edition,
c   1963.
 
	data cta,eis/273.15,6.1071/
 
c   cta = difference between kelvin and celsius temperature
c   eis = saturation vapor pressure (mb) over a water-ice mixture at 0c
 
	data c1,c2,c3/9.09718,3.56654,0.876793/
 
c   c1,c2,c3 = empirical coefficients in the goff-gratch formula
 
	if (t.le.0.) go to 5
	esice = 99999.
	write(6,3)esice
    3	format(' saturation vapor pressure for ice cannot be computed',
     1	       /' for temperature > 0c. esice =',f7.0)
	return
    5	continue
 
c   freezing point of water (k)
 
	tf = cta
	tk = t+cta
 
c   goff-gratch formula
 
	rhs = -c1*(tf/tk-1.)-c2*alog10(tf/tk)+c3*(1.-tk/tf)+alog10(eis)
	esi = 10.**rhs
	if (esi.lt.0.) esi = 0.
	esice = esi
	return
	end
 
c ---------------------------------------------------------------------
 
C NCLFORTSTART
	real function esiloThermo(t)
C       NCL: esilo = esiloThermo (t:float)
c NCLEND
 
c	include 'lib_dev:[gudoc]edfvaxbox.for/list'
c	baker, schlatter  17-may-1982	  original version.
 
c   this function returns the saturation vapor pressure over ice
c   esilo (millibars) given the temperature t (celsius). the formula
c   is due to lowe, paul r., 1977: an approximating polynomial for
c   the computation of saturation vapor pressure, journal of applied
c   meteorology, vol. 16, no. 1 (january), pp. 100-103.
c   the polynomial coefficients are a0 through a6.
 
	data a0,a1,a2,a3,a4,a5,a6
     1	/6.109177956,     5.034698970e-01, 1.886013408e-02,
     2	 4.176223716e-04, 5.824720280e-06, 4.838803174e-08,
     3	 1.838826904e-10/
	if (t.le.0.) go to 5
	esilo = 9999.
	write(6,3)esilo
    3	format(' saturation vapor pressure over ice is undefined for',
     1	/' temperature > 0c. esilo =',f6.0)
	return
    5	continue
	esilo = a0+t*(a1+t*(a2+t*(a3+t*(a4+t*(a5+a6*t)))))
	return
	end
 
c ---------------------------------------------------------------------
 
C NCLFORTSTART
	real function esloThermo(t)
C       NCL: eslo = esloThermo (t:float)
c NCLEND
 
c	include 'lib_dev:[gudoc]edfvaxbox.for/list'
c	baker, schlatter  17-may-1982	  original version.
 
c   this function returns the saturation vapor pressure over liquid
c   water eslo (millibars) given the temperature t (celsius). the
c   formula is due to lowe, paul r.,1977: an approximating polynomial
c   for the computation of saturation vapor pressure, journal of applied
c   meteorology, vol 16, no. 1 (january), pp. 100-103.
c   the polynomial coefficients are a0 through a6.
 
	data a0,a1,a2,a3,a4,a5,a6
     1  /6.107799961,     4.436518521e-01, 1.428945805e-02,
     2   2.650648471e-04, 3.031240396e-06, 2.034080948e-08,
     3	 6.136820929e-11/
	es = a0+t*(a1+t*(a2+t*(a3+t*(a4+t*(a5+a6*t)))))
	if (es.lt.0.) es = 0.
	eslo = es
	return
	end
 
c ---------------------------------------------------------------------
 
C NCLFORTSTART
	real function esrwThermo(t)
C       NCL: esrw = esrwThermo (t:float)
c NCLEND
 
c	include 'lib_dev:[gudoc]edfvaxbox.for/list'
c	baker, schlatter  17-may-1982	  original version.
 
c   this function returns the saturation vapor pressure over liquid
c   water esrw (millibars) given the temperature t (celsius). the
c   formula used is due to richards, j.m., 1971: simple expression
c   for the saturation vapour pressure of water in the range -50 to
c   140c, british journal of applied physics, vol. 4, pp.l15-l18.
c   the formula was quoted more recently by wigley, t.m.l.,1974:
c   comments on 'a simple but accurate formula for the saturation
c   vapor pressure over liquid water,' journal of applied meteorology,
c   vol. 13, no. 5 (august) p.606.
 
	data cta,ts,ews/273.15,373.15,1013.25/
 
c   cta = difference between kelvin and celsius temperature
c   ts = temperature of the boiling point of water (k)
c   ews = saturation vapor pressure over liquid water at 100c
 
	data c1,     c2,     c3,     c4
     1	/ 13.3185,-1.9760,-0.6445,-0.1299 /
	tk = t+cta
	x = 1.-ts/tk
	px = x*(c1+x*(c2+x*(c3+c4*x)))
	vp = ews*exp(px)
	if (vp.lt.0) vp = 0.
	esrw = vp
	return
	end
 
c ---------------------------------------------------------------------
 
C NCLFORTSTART
	real function eswThermo(t)
C       NCL: esw = eswThermo (t:float)
c NCLEND
 
c   this function returns the saturation vapor pressure esw (millibars)
c   over liquid water given the temperature t (celsius). the polynomial
c   approximation below is due to herman wobus, a mathematician who
c   worked at the navy weather research facility, norfolk, virginia,
c   but who is now retired. the coefficients of the polynomial were
c   chosen to fit the values in table 94 on pp. 351-353 of the smith-
c   sonian meteorological tables by roland list (6th edition). the
c   approximation is valid for -50 < t < 100c.
 
c	include 'lib_dev:[gudoc]edfvaxbox.for/list'
c	baker, schlatter  17-may-1982	  original version.
 
c
c   es0 = saturation vapor ressure over liquid water at 0c
 
	data es0/6.1078/
	pol = 0.99999683       + t*(-0.90826951e-02 +
     1	   t*(0.78736169e-04   + t*(-0.61117958e-06 +
     2     t*(0.43884187e-08   + t*(-0.29883885e-10 +
     3     t*(0.21874425e-12   + t*(-0.17892321e-14 +
     4     t*(0.11112018e-16   + t*(-0.30994571e-19)))))))))
	esw = es0/pol**8
	return
	end
 
c ---------------------------------------------------------------------
 
C NCLFORTSTART
	real function esThermo(t)
C       NCL: es = esThermo (t:float)
c NCLEND
 
c   this function returns the saturation vapor pressure es (mb) over
c   liquid water given the temperature t (celsius). the formula appears
c   in bolton, david, 1980: "the computation of equivalent potential
c   temperature," monthly weather review, vol. 108, no. 7 (july),
c   p. 1047, eq.(10). the quoted accuracy is 0.3% or better for
c   -35 < t < 35c.
 
c	include 'lib_dev:[gudoc]edfvaxbox.for/list'
c	baker, schlatter  17-may-1982	  original version.
 
c   es0 = saturation vapor pressure over liquid water at 0c
 
	data es0/6.1121/
	es = es0*exp(17.67*t/(t+243.5))
	return
	end
 
c ---------------------------------------------------------------------
 
C NCLFORTSTART
	real function heatlThermo(key,t)
C       NCL: heatl = heatlThermo (ket:integer, t:float)
c NCLEND
 
c	include 'lib_dev:[gudoc]edfvaxbox.for/list'
c	baker, schlatter  17-may-1982	  original version.
 
c   this function returns the latent heat of
c		evaporation/condensation         for key=1
c		melting/freezing		 for key=2
c		sublimation/deposition		 for key=3
c   for water. the latent heat heatl (joules per kilogram) is a
c   function of temperature t (celsius). the formulas are polynomial
c   approximations to the values in table 92, p. 343 of the smithsonian
c   meteorological tables, sixth revised edition, 1963 by roland list.
c   the approximations were developed by eric smith at colorado state
c   university.
c   polynomial coefficients
 
	data a0,a1,a2/ 3337118.5,-3642.8583, 2.1263947/
	data b0,b1,b2/-1161004.0, 9002.2648,-12.931292/
	data c0,c1,c2/ 2632536.8, 1726.9659,-3.6248111/
	hltnt = 0.
	tk = t+273.15
	if (key.eq.1) hltnt=a0+a1*tk+a2*tk*tk
	if (key.eq.2) hltnt=b0+b1*tk+b2*tk*tk
	if (key.eq.3) hltnt=c0+c1*tk+c2*tk*tk
	heatl = hltnt
	return
	end
 
c ---------------------------------------------------------------------
 
C NCLFORTSTART
	real function humThermo(t,td)
C       NCL: hum = humThermo (t:float, td:float)
c NCLEND
 
c	include 'lib_dev:[gudoc]edfvaxbox.for/list'
c	g.s. stipanuk     1973      	  original version.
c	reference stipanuk paper entitled:
c            "algorithms for generating a skew-t, log p
c	     diagram and computing selected meteorological
c	     quantities."
c	     atmospheric sciences laboratory
c	     u.s. army electronics command
c	     white sands missile range, new mexico 88002
c	     33 pages
c	baker, schlatter  17-may-1982	
 
c   this function returns relative humidity (%) given the
c   temperature t and dew point td (celsius).  as calculated here,
c   relative humidity is the ratio of the actual vapor pressure to
c   the saturation vapor pressure.
 
	hum= 100.*(esat(td)/esat(t))
	return
	end
 
c ---------------------------------------------------------------------
 
C NCLFORTSTART
	real function oeThermo(t,td,p)
C       NCL: oe = oeThermo (t:float, td:float, p:float)
c NCLEND
 
c	include 'lib_dev:[gudoc]edfvaxbox.for/list'
c	g.s. stipanuk     1973      	  original version.
c	reference stipanuk paper entitled:
c            "algorithms for generating a skew-t, log p
c	     diagram and computing selected meteorological
c	     quantities."
c	     atmospheric sciences laboratory
c	     u.s. army electronics command
c	     white sands missile range, new mexico 88002
c	     33 pages
c	baker, schlatter  17-may-1982	
 
c   this function returns equivalent potential temperature oe (celsius)
c   of a parcel of air given its temperature t (celsius), dew point
c   td (celsius) and pressure p (millibars).
c   find the wet bulb temperature of the parcel.
 
	atw = tw(t,td,p)
 
c   find the equivalent potential temperature.
 
	oe = os(atw,p)
	return
	end
 
c ---------------------------------------------------------------------
 
C NCLFORTSTART
	real function osThermo(t,p)
C       NCL: os = osThermo (t:float, p:float)
c NCLEND
 
c	include 'lib_dev:[gudoc]edfvaxbox.for/list'
c	g.s. stipanuk     1973      	  original version.
c	reference stipanuk paper entitled:
c            "algorithms for generating a skew-t, log p
c	     diagram and computing selected meteorological
c	     quantities."
c	     atmospheric sciences laboratory
c	     u.s. army electronics command
c	     white sands missile range, new mexico 88002
c	     33 pages
c	baker, schlatter  17-may-1982	
 
c   this function returns the equivalent potential temperature os
c   (celsius) for a parcel of air saturated at temperature t (celsius)
c   and pressure p (millibars).
	data b/2.6518986/
c   b is an empirical constant approximately equal to the latent heat
c   of vaporization for water divided by the specific heat at constant
c   pressure for dry air.
 
	tk = t+273.15
	osk= tk*((1000./p)**.286)*(exp(b*wx(t,p)/tk))
	os= osk-273.15
	return
	end
 
c ---------------------------------------------------------------------
 
C NCLFORTSTART
	real function owThermo(t,td,p)
C       NCL: ow = owThermo (t:float, td:float, p:float)
c NCLEND
 
c	include 'lib_dev:[gudoc]edfvaxbox.for/list'
c	g.s. stipanuk     1973      	  original version.
c	reference stipanuk paper entitled:
c            "algorithms for generating a skew-t, log p
c	     diagram and computing selected meteorological
c	     quantities."
c	     atmospheric sciences laboratory
c	     u.s. army electronics command
c	     white sands missile range, new mexico 88002
c	     33 pages
c	baker, schlatter  17-may-1982	
 
c   this function returns the wet-bulb potential temperature ow
c   (celsius) given the temperature t (celsius), dew point td
c   (celsius), and pressure p (millibars).  the calculation for ow is
c   very similar to that for wet bulb temperature. see p.13 stipanuk (1973).
c   find the wet bulb temperature of the parcel
 
	atw = tw(t,td,p)
 
c   find the equivalent potential temperature of the parcel.
 
	aos= os(atw,p)
 
c   find the wet-bulb potential temperature.
 
	ow= tsa(aos,1000.)
	return
	end
 
c ---------------------------------------------------------------------
 
C NCLFORTSTART
	real function oxThermo(t,p)
C       NCL: ox = oxThermo (t:float, p:float)
c       function o(t,p)           ! name changed
c NCLEND
 
c	include 'lib_dev:[gudoc]edfvaxbox.for/list'
c	g.s. stipanuk     1973      	  original version.
c	reference stipanuk paper entitled:
c            "algorithms for generating a skew-t, log p
c	     diagram and computing selected meteorological
c	     quantities."
c	     atmospheric sciences laboratory
c	     u.s. army electronics command
c	     white sands missile range, new mexico 88002
c	     33 pages
c	baker, schlatter  17-may-1982	
 
c   this function returns potential temperature (celsius) given
c   temperature t (celsius) and pressure p (mb) by solving the poisson
c   equation.
 
	tk= t+273.15
	ok= tk*((1000./p)**.286)
	ox= ok-273.15
	return
	end
 
c ---------------------------------------------------------------------
 
C NCLFORTSTART
	real function pcclThermo(pm,p,t,td,wbar,n)
	dimension t(n),p(n),td(n)
C       NCL: pccl = pcclThermo (pm:float, p[*]:float, t[*]:float, td[*]:float\
C                              ,wbar:float)
C            Note: n= dimsizes(p)
c NCLEND
 
c	include 'lib_dev:[gudoc]edfvaxbox.for/list'
c	g.s. stipanuk     1973      	  original version.
c	reference stipanuk paper entitled:
c            "algorithms for generating a skew-t, log p
c	     diagram and computing selected meteorological
c	     quantities."
c	     atmospheric sciences laboratory
c	     u.s. army electronics command
c	     white sands missile range, new mexico 88002
c	     33 pages
c	baker, schlatter  17-may-1982	
 
c   this function returns the pressure at the convective condensation
c   level given the appropriate sounding data.
c   on input:
c	p = pressure (millibars). note that p(i).gt.p(i+1).
c	t = temperature (celsius)
c	td = dew point (celsius)
c	n = number of levels in the sounding and the dimension of
c	    p, t and td
c	pm = pressure (millibars) at upper boundary of the layer for
c	     computing the mean mixing ratio. p(1) is the lower
c	     boundary.
c   on output:
c	pccl = pressure (millibars) at the convective condensation level
c	wbar = mean mixing ratio (g/kg) in the layer bounded by
c	       pressures p(1) at the bottom and pm at the top
c   the algorithm is decribed on p.17 of stipanuk, g.s.,1973:
c   "algorithms for generating a skew-t log p diagram and computing
c   selected meteorological quantities," atmospheric sciences labora-
c   tory, u.s. army electronics command, white sands missile range, new
c   mexico 88002.
 
	if (pm.ne.p(1)) go to 5
	wbar= wx(td(1),p(1))
	pc= pm
	if (abs(t(1)-td(1)).lt.0.05) go to 45
	go to 25
    5	continue
	wbar= 0.
	k= 0
   10	continue
	k = k+1
	if (p(k).gt.pm) go to 10
	k= k-1
	j= k-1
	if(j.lt.1) go to 20
 
c   compute the average mixing ratio....alog = natural log
 
	do 15 i= 1,j
	   l= i+1
   15      wbar= (wx(td(i),p(i))+wx(td(l),p(l)))*alog(p(i)/p(l))
     *          +wbar
   20	continue
	l= k+1
 
c   estimate the dew point at pressure pm.
 
	tq= td(k)+(td(l)-td(k))*(alog(pm/p(k)))/(alog(p(l)/p(k)))
	wbar= wbar+(wx(td(k),p(k))+wx(tq,pm))*alog(p(k)/pm)
	wbar= wbar/(2.*alog(p(1)/pm))
 
c   find level at which the mixing ratio line wbar crosses the
c   environmental temperature profile.
 
   25	continue
	do 30 j= 1,n
	   i= n-j+1
	   if (p(i).lt.300.) go to 30
 
c   tmr = temperature (celsius) at pressure p (mb) along a mixing
c   	  ratio line given by wbar (g/kg)
 
	   x= tmr(wbar,p(i))-t(i)
	   if (x.le.0.) go to 35
   30	continue
	pccl= 0.0
	return
 
c  set up bisection routine
 
   35	l = i
	i= i+1
	del= p(l)-p(i)
	pc= p(i)+.5*del
	a= (t(i)-t(l))/alog(p(l)/p(i))
	do 40 j= 1,10
	   del= del/2.
	   x= tmr(wbar,pc)-t(l)-a*(alog(p(l)/pc))
 
c   the sign function replaces the sign of the first argument
c   with that of the second.
 
   40	pc= pc+sign(del,x)
   45	pccl = pc
	return
	end
 
c ---------------------------------------------------------------------
 
C NCLFORTSTART
	real function pconThermo(p,t,tc)
C       NCL: pcon = pconThermo (p:float, t:float, tc:float)
c NCLEND
 
c   this function returns the pressure pcon (mb) at the lifted condensa-
c   tion level, given the initial pressure p (mb) and temperature t
c   (celsius) of the parcel and the temperature tc (celsius) at the
c   lcl. the algorithm is exact.  it makes use of the formula for the
c   potential temperatures corresponding to t at p and tc at pcon.
c   these two potential temperatures are equal.
 
c	include 'lib_dev:[gudoc]edfvaxbox.for/list'
c	baker, schlatter  17-may-1982	  original version.
 
c
	data akapi/3.5037/
 
c   akapi = (specific heat at constant pressure for dry air) /
c	    (gas constant for dry air)
 
c   convert t and tc to kelvin temperatures.
 
	tk = t+273.15
	tck = tc+273.15
	pcon = p*(tck/tk)**akapi
	return
	end
 
c ---------------------------------------------------------------------
 
C NCLFORTSTART
	real function powtThermo(t,p,td)
C       NCL: powt = powtThermo (t:float, p:float, td:float)
c NCLEND
 
c   this function yields wet-bulb potential temperature powt
c   (celsius), given the following input:
c	   t = temperature (celsius)
c	   p = pressure (millibars)
c	   td = dew point (celsius)
 
c	include 'lib_dev:[gudoc]edfvaxbox.for/list'
c	baker, schlatter  17-may-1982	  original version.
 
	data cta,akap/273.15,0.28541/
 
c   cta = difference between kelvin and celsius temperatures
c   akap = (gas constant for dry air) / (specific heat at
c	   constant pressure for dry air)
c   compute the potential temperature (celsius)
 
	pt = (t+cta)*(1000./p)**akap-cta
 
c   compute the lifting condensation level (lcl).
 
	tc = tcon(t,td)
 
c   for the origin of the following approximation, see the documen-
c   tation for the wobus function.
 
	powt = pt-wobf(pt)+wobf(tc)
	return
	end
 
c ---------------------------------------------------------------------
 
C NCLFORTSTART
	real function precpwThermo(td,p,n)
	dimension td(n),p(n)
C       NCL: precpw = precpwThermo (td[*]:float, p[*]:float)
C            Note: n= dimsizes(p)
c NCLEND
 
c	include 'lib_dev:[gudoc]edfvaxbox.for/list'
c	baker, schlatter  17-may-1982	  original version.
 
c   this function computes total precipitable water precpw (cm) in a
c   vertical column of air based upon sounding data at n levels:
c	   td = dew point (celsius)
c	   p = pressure (millibars)
c   calculations are done in cgs units.
 
 
c   g = acceleration due to the earth's gravity (cm/s**2)
 
	data g/980.616/
 
c   initialize value of precipitable water
 
	pw = 0.
	nl = n-1
 
c   calculate the mixing ratio at the lowest level.
 
	wbot = wmr(p(1),td(1))
	do 5 i=1,nl
	wtop = wmr(p(i+1),td(i+1))
 
c   calculate the layer-mean mixing ratio (g/kg).
 
	w = 0.5*(wtop+wbot)
 
c   make the mixing ratio dimensionless.
 
	wl = .001*w
 
c   calculate the specific humidity.
 
	ql = wl/(wl+1.)
 
c   the factor of 1000. below converts from millibars to dynes/cm**2.
 
	dp = 1000.*(p(i)-p(i+1))
	pw = pw+(ql/g)*dp
	wbot = wtop
    5	continue
	precpw = pw
	return
	end
 
c ---------------------------------------------------------------------
 
C NCLFORTSTART
	subroutine ptlclThermo(p,t,td,pc,tc)
C       NCL: ptlcl = ptlclThermo (p:float, t:float, td:float, pc:float, tc:float)
c NCLEND
 
c   this subroutine estimates the pressure pc (mb) and the temperature
c   tc (celsius) at the lifted condensation level (lcl), given the
c   initial pressure p (mb), temperature t (celsius) and dew point
c   (celsius) of the parcel.  the approximation is that lines of
c   constant potential temperature and constant mixing ratio are
c   straight on the skew t/log p chart.
c   teten's formula for saturation vapor pressure as a function of
c   pressure was used in the derivation of the formula below.  for
c   additional details, see math notes by t. schlatter dated 8 sep 81.
c   t. schlatter, noaa/erl/profs program office, boulder, colorado,
c   wrote this subroutine.
 
c	include 'lib_dev:[gudoc]edfvaxbox.for/list'
c	baker, schlatter  17-may-1982	  original version.
 
c   akap = (gas constant for dry air) / (specific heat at constant
c	   pressure for dry air)
c   cta = difference between kelvin and celsius temperatures
 
	data akap,cta/0.28541,273.15/
 
	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
 
c ---------------------------------------------------------------------
 
C NCLFORTSTART
	real function satlftThermo(thw,p)
C       NCL: satlft = satlftThermo (thw:float, p:float)
c NCLEND
 
c	include 'lib_dev:[gudoc]edfvaxbox.for/list'
c	baker, schlatter  17-may-1982	  original version.
 
c   input:  thw = wet-bulb potential temperature (celsius).
c		  thw defines a moist adiabat.
c	    p = pressure (millibars)
c   output: satlft = temperature (celsius) where the moist adiabat
c		  crosses p
 
	data cta,akap/273.15,0.28541/
 
c   cta = difference between kelvin and celsius temperatures
c   akap = (gas constant for dry air) / (specific heat at constant
c	    pressure for dry air)
 
c        the algorithm below can best be understood by referring to a
c   skew-t/log p chart.  it was devised by herman wobus, a mathemati-
c   cian formerly at the navy weather research facility but now retired.
c   the value returned by satlft can be checked by referring to table
c   78, pp.319-322, smithsonian meteorological tables, by roland list
c   (6th revised edition).
c
 
	if (p.ne.1000.) go to 5
	satlft = thw
	return
    5	continue
 
c   compute tone, the temperature where the dry adiabat with value thw
c   (celsius) crosses p.
 
	pwrp = (p/1000.)**akap
	tone = (thw+cta)*pwrp-cta
 
c   consider the moist adiabat ew1 through tone at p.  using the defini-
c   tion of the wobus function (see documentation on wobf), it can be
c   shown that eone = ew1-thw.
 
	eone = wobf(tone)-wobf(thw)
	rate = 1.
	go to 15
 
c   in the loop below, the estimate of satlft is iteratively improved.
 
   10	continue
 
c   rate is the ratio of a change in t to the corresponding change in
c   e.  its initial value was set to 1 above.
 
	rate = (ttwo-tone)/(etwo-eone)
	tone = ttwo
	eone = etwo
   15	continue
 
c   ttwo is an improved estimate of satlft.
 
	ttwo = tone-eone*rate
 
c   pt is the potential temperature (celsius) corresponding to ttwo at p
 
	pt = (ttwo+cta)/pwrp-cta
 
c   consider the moist adiabat ew2 through ttwo at p. using the defini-
c   tion of the wobus function, it can be shown that etwo = ew2-thw.
 
	etwo = pt+wobf(ttwo)-wobf(pt)-thw
 
c   dlt is the correction to be subtracted from ttwo.
 
	dlt = etwo*rate
	if (abs(dlt).gt.0.1) go to 10
	satlft = ttwo-dlt
	return
	end
 
c ---------------------------------------------------------------------
 
C NCLFORTSTART
	real function sshThermo(p,t)
C       NCL: sh = sshThermo (p:float, t:float)
c NCLEND
 
c	include 'lib_dev:[gudoc]edfvaxbox.for/list'
c	baker, schlatter  17-may-1982	  original version.
 
c   this function returns saturation specific humidity ssh (grams of
c   water vapor per kilogram of moist air) given the pressure p
c   (millibars) and the temperature t (celsius). the equation is given
c   in standard meteorological texts. if t is dew point (celsius), then
c   ssh returns the actual specific humidity.
c   compute the dimensionless mixing ratio.
 
	w = .001*wmr(p,t)
 
c   compute the dimensionless saturation specific humidity.
 
	q = w/(1.+w)
	ssh = 1000.*q
	return
	end
 
c ---------------------------------------------------------------------
 
C NCLFORTSTART
	real function tconThermo(t,d)
C       NCL: tcon = tconThermo (t:float, d:float)
c NCLEND
 
c   this function returns the temperature tcon (celsius) at the lifting
c   condensation level, given the temperature t (celsius) and the
c   dew point d (celsius).
 
c	include 'lib_dev:[gudoc]edfvaxbox.for/list'
c	baker, schlatter  17-may-1982	  original version.
 
c
c   compute the dew point depression s.
	s = t-d
c   the approximation below, a third order polynomial in s and t,
c   is due to herman wobus. the source of data for fitting the
c   polynomial is unknown.
 
	dlt = s*(1.2185+1.278e-03*t+
     1        s*(-2.19e-03+1.173e-05*s-5.2e-06*t))
	tcon = t-dlt
	return
	end
 
c ---------------------------------------------------------------------
 
C NCLFORTSTART
	real function tdaThermo(o,p)
C       NCL: tda = tdaThermo (o:float, p:float)
c NCLEND
 
c	include 'lib_dev:[gudoc]edfvaxbox.for/list'
c	g.s. stipanuk     1973      	  original version.
c	reference stipanuk paper entitled:
c            "algorithms for generating a skew-t, log p
c	     diagram and computing selected meteorological
c	     quantities."
c	     atmospheric sciences laboratory
c	     u.s. army electronics command
c	     white sands missile range, new mexico 88002
c	     33 pages
c	baker, schlatter  17-may-1982	
 
c   this function returns the temperature tda (celsius) on a dry adiabat
c   at pressure p (millibars). the dry adiabat is given by
c   potential temperature o (celsius). the computation is based on
c   poisson's equation.
 
	ok= o+273.15
	tdak= ok*((p*.001)**.286)
	tda= tdak-273.15
	return
	end
 
c ---------------------------------------------------------------------
 
C NCLFORTSTART
	real function teThermo(t,td,p)
C       NCL: te = teThermo (t:float, td:float, p:float)
c NCLEND
 
c	include 'lib_dev:[gudoc]edfvaxbox.for/list'
c	g.s. stipanuk     1973      	  original version.
c	reference stipanuk paper entitled:
c            "algorithms for generating a skew-t, log p
c	     diagram and computing selected meteorological
c	     quantities."
c	     atmospheric sciences laboratory
c	     u.s. army electronics command
c	     white sands missile range, new mexico 88002
c	     33 pages
c	baker, schlatter  17-may-1982	
 
c   this function returns the equivalent temperature te (celsius) of a
c   parcel of air given its temperature t (celsius), dew point (celsius)
c   and pressure p (millibars).
c   calculate equivalent potential temperature.
 
	aoe = oe(t,td,p)
 
c   use poissons's equation to calculate equivalent temperature.
 
	te= tda(aoe,p)
	return
	end
 
c ---------------------------------------------------------------------
 
C NCLFORTSTART
	real function thmThermo(t,p)
C       NCL: thm = thmThermo (t:float, p:float)
c NCLEND
 
c	include 'lib_dev:[gudoc]edfvaxbox.for/list'
c	baker, schlatter  17-may-1982	  original version.
 
c   this function returns the wet-bulb potential temperature thm
c   (celsius) corresponding to a parcel of air saturated at
c   temperature t (celsius) and pressure p (millibars).
 
	f(x) =   1.8199427e+01+x*( 2.1640800e-01+x*( 3.0716310e-04+x*
     1	       (-3.8953660e-06+x*( 1.9618200e-08+x*( 5.2935570e-11+x*
     2	       ( 7.3995950e-14+x*(-4.1983500e-17)))))))
	thm = t
	if (p.eq.1000.) return
 
c   compute the potential temperature (celsius).
 
	thd = (t+273.15)*(1000./p)**.286-273.15
	thm = thd+6.071*(exp(t/f(t))-exp(thd/f(thd)))
	return
	end
 
c ---------------------------------------------------------------------
 
C NCLFORTSTART
	real function tlcl1Thermo(t,td)
C       NCL: tlcl1 = tlcl1Thermo (t:float, td:float)
c NCLEND
 
c	include 'lib_dev:[gudoc]edfvaxbox.for/list'
c	baker, schlatter  17-may-1982	  original version.
 
c   this function returns the temperature tlcl1 (celsius) of the lifting
c   condensation level (lcl) given the initial temperature t (celsius)
c   and dew point td (celsius) of a parcel of air.
c   eric smith at colorado state university has used the formula
c   below, but its origin is unknown.
 
	data cta/273.15/
 
c   cta = difference between kelvin and celsius temperature
 
	tk = t+cta
 
c   compute the parcel vapor pressure (mb).
	es = eslo(td)
	tlcl = 2840./(3.5*alog(tk)-alog(es)-4.805)+55.
	tlcl1 = tlcl-cta
	return
	end
 
c ---------------------------------------------------------------------
 
C NCLFORTSTART
	real function tlclThermo(t,td)
C       NCL: tlcl = tlclThermo (t:float, td:float)
c NCLEND
c   this function yields the temperature tlcl (celsius) of the lifting
c   condensation level, given the temperature t (celsius) and the
c   dew point td (celsius).  the formula used is in bolton, david,
c   1980: "the computation of equivalent potential temperature,"
c   monthly weather review, vol. 108, no. 7 (july), p. 1048, eq.(15).
 
c	include 'lib_dev:[gudoc]edfvaxbox.for/list'
c	baker, schlatter  17-may-1982	  original version.
 
c   convert from celsius to kelvin degrees.
 
	tk = t+273.15
	tdk = td+273.15
	a = 1./(tdk-56.)
	b = alog(tk/tdk)/800.
	tc = 1./(a+b)+56.
	tlcl = tc-273.15
	return
	end
 
c ---------------------------------------------------------------------
 
C NCLFORTSTART
	real function tmlapsThermo(thetae,p)
C       NCL: tmlaps = tmlapsThermo (thetae:float, p:float)
c NCLEND
 
c	include 'lib_dev:[gudoc]edfvaxbox.for/list'
c	baker, schlatter  17-may-1982	  original version.
 
c   this function returns the temperature tmlaps (celsius) at pressure
c   p (millibars) along the moist adiabat corresponding to an equivalent
c   potential temperature thetae (celsius).
c   the algorithm was written by eric smith at colorado state
c   university.
 
	data crit/0.1/
c   cta = difference between kelvin and celsius temperatures.
c   crit = convergence criterion (degrees kelvin)
 
	eq0 = thetae
 
c   initial guess for solution
 
	tlev = 25.
 
c   compute the saturation equivalent potential temperature correspon-
c   ding to temperature tlev and pressure p.
 
	eq1 = ept(tlev,tlev,p)
	dif = abs(eq1-eq0)
	if (dif.lt.crit) go to 3
	if (eq1.gt.eq0) go to 1
 
c   dt is the initial stepping increment.
 
	dt = 10.
	i = -1
	go to 2
    1	dt = -10.
	i = 1
    2	tlev = tlev+dt
	eq1 = ept(tlev,tlev,p)
	dif = abs(eq1-eq0)
	if (dif.lt.crit) go to 3
	j = -1
	if (eq1.gt.eq0) j=1
	if (i.eq.j) go to 2
 
c   the solution has been passed. reverse the direction of search
c   and decrease the stepping increment.
 
	tlev = tlev-dt
	dt = dt/10.
	go to 2
    3	tmlaps = tlev
	return
	end
 
c ---------------------------------------------------------------------
 
C NCLFORTSTART
	real function tmrThermo(w,p)
C       NCL: tmr = tmrThermo (w:float, p:float)
c NCLEND
 
c	include 'lib_dev:[gudoc]edfvaxbox.for/list'
c	g.s. stipanuk     1973      	  original version.
c	reference stipanuk paper entitled:
c            "algorithms for generating a skew-t, log p
c	     diagram and computing selected meteorological
c	     quantities."
c	     atmospheric sciences laboratory
c	     u.s. army electronics command
c	     white sands missile range, new mexico 88002
c	     33 pages
c	baker, schlatter  17-may-1982	
 
c   this function returns the temperature (celsius) on a mixing
c   ratio line w (g/kg) at pressure p (mb). the formula is given in
c   table 1 on page 7 of stipanuk (1973).
c
c   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
 
c ---------------------------------------------------------------------
 
C NCLFORTSTART
	real function tsaThermo(os,p)
C       NCL: tsa = tsaThermo (os:float, p:float)
c NCLEND
 
c	include 'lib_dev:[gudoc]edfvaxbox.for/list'
c	g.s. stipanuk     1973      	  original version.
c	reference stipanuk paper entitled:
c            "algorithms for generating a skew-t, log p
c	     diagram and computing selected meteorological
c	     quantities."
c	     atmospheric sciences laboratory
c	     u.s. army electronics command
c	     white sands missile range, new mexico 88002
c	     33 pages
c	baker, schlatter  17-may-1982	
 
c   this function returns the temperature tsa (celsius) on a saturation
c   adiabat at pressure p (millibars). os is the equivalent potential
c   temperature of the parcel (celsius). sign(a,b) replaces the
c   algebraic sign of a with that of b.
c   b is an empirical constant approximately equal to 0.001 of the latent
c   heat of vaporization for water divided by the specific heat at constant
c   pressure for dry air.
 
	data b/2.6518986/
	a= os+273.15
 
c   tq is the first guess for tsa.
 
	tq= 253.15
 
c   d is an initial value used in the iteration below.
 
	d= 120.
 
c   iterate to obtain sufficient accuracy....see table 1, p.8
c   of stipanuk (1973) for equation used in iteration.
 
	do 1 i= 1,12
	   tqk= tq-273.15
	   d= d/2.
	   x= a*exp(-b*wx(tqk,p)/tq)-tq*((1000./p)**.286)
	   if (abs(x).lt.1e-7) go to 2
	   tq= tq+sign(d,x)
 1	continue
 2	tsa= tq-273.15
	return
	end
 
c ---------------------------------------------------------------------
 
C NCLFORTSTART
	real function tvirtThermo(t,td,p)
c c c   real function tv(t,td,p)         ! name changed
C       NCL: tvirt = tvirtThermo (t:float, td:float, p:float)
c NCLEND
 
c	include 'lib_dev:[gudoc]edfvaxbox.for/list'
c	baker, schlatter  17-may-1982	  original version.
 
c   this function returns the virtual temperature tvirt (celsius) of
c   a parcel of air at temperature t (celsius), dew point td
c   (celsius), and pressure p (millibars). the equation appears
c   in most standard meteorological texts.
 
	data cta,eps/273.15,0.62197/
 
c   cta = difference between kelvin and celsius temperatures.
c   eps = ratio of the mean molecular weight of water (18.016 g/mole)
c	  to that of dry air (28.966 g/mole)
 
	tk = t+cta
 
c   calculate the dimensionless mixing ratio.
 
	w = .001*wmr(p,td)
	tvirt = tk*(1.+w/eps)/(1.+w)-cta
        return
        end
 
c ---------------------------------------------------------------------
 
C NCLFORTSTART
	real function twThermo(t,td,p)
C       NCL: tw = twThermo (t:float, td:float, p:float)
c NCLEND
 
c	include 'lib_dev:[gudoc]edfvaxbox.for/list'
c	g.s. stipanuk     1973      	  original version.
c	reference stipanuk paper entitled:
c            "algorithms for generating a skew-t, log p
c	     diagram and computing selected meteorological
c	     quantities."
c	     atmospheric sciences laboratory
c	     u.s. army electronics command
c	     white sands missile range, new mexico 88002
c	     33 pages
c	baker, schlatter  17-may-1982	
 
c   this function returns the wet-bulb temperature tw (celsius)
c   given the temperature t (celsius), dew point td (celsius)
c   and pressure p (mb).  see p.13 in stipanuk (1973), referenced
c   above, for a description of the technique.
c
c
c   determine the mixing ratio line thru td and p.
 
	aw = wx(td,p)
c
c   determine the dry adiabat thru t and p.
 
	ao = ox(t,p)               ! was ao = o(t,p)
	pi = p
 
c   iterate to locate pressure pi at the intersection of the two
c   curves .  pi has been set to p for the initial guess.
 
	do 4 i= 1,10
	   x= .02*(tmr(aw,pi)-tda(ao,pi))
	   if (abs(x).lt.0.01) go to 5
 4	   pi= pi*(2.**(x))
 
c   find the temperature on the dry adiabat ao at pressure pi.
 
 5	ti= tda(ao,pi)
 
c   the intersection has been located...now, find a saturation
c   adiabat thru this point. function os returns the equivalent
c   potential temperature (c) of a parcel saturated at temperature
c   ti and pressure pi.
 
	aos= os(ti,pi)
 
c   function tsa returns the wet-bulb temperature (c) of a parcel at
c   pressure p whose equivalent potential temperature is aos.
 
	tw = tsa(aos,p)
	return
	end
 
c ---------------------------------------------------------------------
 
C NCLFORTSTART
	real function wmrThermo(p,t)
C       NCL: wmr = wmrThermo (p:float, t:float)
c NCLEND
 
c   this function approximates the mixing ratio wmr (grams of water
c   vapor per kilogram of dry air) given the pressure p (mb) and the
c   temperature t (celsius). the formula used is given on p. 302 of the
c   smithsonian meteorological tables by roland list (6th edition).
 
c	baker, schlatter  17-may-1982	  original version.
 
c   eps = ratio of the mean molecular weight of water (18.016 g/mole)
c         to that of dry air (28.966 g/mole)
 
	data eps/0.62197/
 
c   the next two lines contain a formula by herman wobus for the
c   correction factor wfw for the departure of the mixture of air
c   and water vapor from the ideal gas law. the formula fits values
c   in table 89, p. 340 of the smithsonian meteorological tables,
c   but only for temperatures and pressures normally encountered in
c   in the atmosphere.
 
	x = 0.02*(t-12.5+7500./p)
	wfw = 1.+4.5e-06*p+1.4e-03*x*x
	fwesw = wfw*esw(t)
	r = eps*fwesw/(p-fwesw)
 
c   convert r from a dimensionless ratio to grams/kilogram.
 
	wmr = 1000.*r
	return
	end
 
c ---------------------------------------------------------------------
 
C NCLFORTSTART
	real function wobfThermo(t)
C       NCL: wobf = wobfThermo (t:float)
c NCLEND
 
c   this function calculates the difference of the wet-bulb potential
c   temperatures for saturated and dry air given the temperature.
 
c	include 'lib_dev:[gudoc]edfvaxbox.for/list'
c	baker, schlatter  17-may-1982	  original version.
 
c        let wbpts = wet-bulb potential temperature for saturated
c   air at temperature t (celsius). let wbptd = wet-bulb potential
c   temperature for completely dry air at the same temperature t.
c   the wobus function wobf (in degrees celsius) is defined by
c                      wobf(t) = wbpts-wbptd.
c   although wbpts and wbptd are functions of both pressure and
c   temperature, their difference is a function of temperature only.
 
c        to understand why, consider a parcel of dry air at tempera-
c   ture t and pressure p. the thermodynamic state of the parcel is
c   represented by a point on a pseudoadiabatic chart. the wet-bulb
c   potential temperature curve (moist adiabat) passing through this
c   point is wbpts. now t is the equivalent temperature for another
c   parcel saturated at some lower temperature tw, but at the same
c   pressure p.  to find tw, ascend along the dry adiabat through
c   (t,p). at a great height, the dry adiabat and some moist
c   adiabat will nearly coincide. descend along this moist adiabat
c   back to p. the parcel temperature is now tw. the wet-bulb
c   potential temperature curve (moist adiabat) through (tw,p) is wbptd.
c   the difference (wbpts-wbptd) is proportional to the heat imparted
c   to a parcel saturated at temperature tw if all its water vapor
c   were condensed. since the amount of water vapor a parcel can
c   hold depends upon temperature alone, (wbptd-wbpts) must depend
c   on temperature alone.
 
c        the wobus function is useful for evaluating several thermo-
c   dynamic quantities.  by definition:
c		    wobf(t) = wbpts-wbptd.               (1)
c   if t is at 1000 mb, then t is a potential temperature pt and
c   wbpts = pt. thus
c		    wobf(pt) = pt-wbptd.                 (2)
c   if t is at the condensation level, then t is the condensation
c   temperature tc and wbpts is the wet-bulb potential temperature
c   wbpt. thus
c		    wobf(tc) = wbpt-wbptd.               (3)
c   if wbptd is eliminated from (2) and (3), there results
c		    wbpt = pt-wobf(pt)+wobf(tc).
c   if wbptd is eliminated from (1) and (2), there results
c		    wbpts = pt-wobf(pt)+wobf(t).
 
c        if t is an equivalent potential temperature ept (implying
c   that the air at 1000 mb is completely dry), then wbpts = ept
c   and wbptd = wbpt. thus
c		    wobf(ept) = ept-wbpt.
c   this form is the basis for a polynomial approximation to wobf.
c   in table 78 on pp.319-322 of the smithsonian meteorological
c   tables by roland list (6th revised edition), one finds wet-bulb
c   potential temperatures and the corresponding equivalent potential
c   temperatures listed together. herman wobus, a mathematician for-
c   merly at the navy weather research facility, norfolk, virginia,
c   and now retired, computed the coefficients for the polynomial
c   approximation from numbers in this table.
c
c                                    notes by t.w. schlatter
c                                    noaa/erl/profs program office
c                                    august 1981
 
	x = t-20.
	if (x.gt.0.) go to 10
	pol = 1.                 +x*(-8.8416605e-03
     1	     +x*( 1.4714143e-04  +x*(-9.6719890e-07
     2       +x*(-3.2607217e-08  +x*(-3.8598073e-10)))))
	wobf = 15.130/pol**4
	return
   10	continue
	pol = 1.                 +x*( 3.6182989e-03
     1       +x*(-1.3603273e-05  +x*( 4.9618922e-07
     2       +x*(-6.1059365e-09  +x*( 3.9401551e-11
     3       +x*(-1.2588129e-13  +x*( 1.6688280e-16)))))))
	wobf = 29.930/pol**4+0.96*x-14.8
	return
	end
 
c ---------------------------------------------------------------------
 
C NCLFORTSTART
	real function wxThermo(t,p)
C       NCL: wx = wxThermo (t:float, p:float)
c       real function w(t,p)          ! name changed
c NCLEND
 
c	include 'lib_dev:[gudoc]edfvaxbox.for/list'
c	g.s. stipanuk     1973      	  original version.
c	reference stipanuk paper entitled:
c            "algorithms for generating a skew-t, log p
c	     diagram and computing selected meteorological
c	     quantities."
c	     atmospheric sciences laboratory
c	     u.s. army electronics command
c	     white sands missile range, new mexico 88002
c	     33 pages
c	baker, schlatter  17-may-1982	
 
c  this function returns the mixing ratio (grams of water vapor per
c  kilogram of dry air) given the dew point (celsius) and pressure
c  (millibars). if the temperture  is input instead of the
c  dew point, then saturation mixing ratio (same units) is returned.
c  the formula is found in most meteorological texts.
 
	x = esat(t)
	wx= 622.*x/(p-x)
	return
	end
 
c ---------------------------------------------------------------------
 
C NCLFORTSTART
	real function zthkThermo(pt,p,t,td,n)
c       real function z(pt,p,t,td,n)     ! name changed
	dimension t(n),p(n),td(n)
C       NCL: zthk = zthkThermo (pt:float, p[*]:float, t[*]:float \
C                              ,td[*]:float)
C            Note: n= dimsizes(p)
c NCLEND
 
c	include 'lib_dev:[gudoc]edfvaxbox.for/list'
c	g.s. stipanuk     1973      	  original version.
c	reference stipanuk paper entitled:
c            "algorithms for generating a skew-t, log p
c	     diagram and computing selected meteorological
c	     quantities."
c	     atmospheric sciences laboratory
c	     u.s. army electronics command
c	     white sands missile range, new mexico 88002
c	     33 pages
c	baker, schlatter  17-may-1982	
 
c   this function returns the thickness of a layer bounded by pressure
c   p(1) at the bottom and pressure pt at the top.
c   on input:
c	p = pressure (mb).  note that p(i).gt.p(i+1).
c	t = temperature (celsius)
c	td = dew point (celsius)
c	n = number of levels in the sounding and the dimension of
c	    p, t and td
c   on output:
c	z = geometric thickness of the layer (m)
c   the algorithm involves numerical integration of the hydrostatic
c   equation from p(1) to pt. it is described on p.15 of stipanuk
c   (1973).
 
	dimension tk(100)
 
c	c1 = .001*(1./eps-1.) where eps = .62197 is the ratio of the
c			      molecular weight of water to that of
c			      dry air. the factor 1000. converts the
c			      mixing ratio w from g/kg to a dimension-
c			      less ratio.
c	c2 = r/(2.*g) where r is the gas constant for dry air
c		      (287 kg/joule/deg k) and g is the acceleration
c		      due to the earth's gravity (9.8 m/s**2). the
c		      factor of 2 is used in averaging two virtual
c		      temperatures.
 
	data c1/.0006078/,c2/14.64285/
	do 5 i= 1,n
	   tk(i)= t(i)+273.15
    5	continue
        z   = 0.0              ! ? djs ?
	zthk= 0.0
	if (pt.lt.p(n)) go to 20
	i= 0
   10	i= i+1
	j= i+1
	if (pt.ge.p(j)) go to 15
	a1= tk(j)*(1.+c1*wx(td(j),p(j)))
	a2= tk(i)*(1.+c1*wx(td(i),p(i)))
	z= z+c2*(a1+a2)*(alog(p(i)/p(j)))
        zthk = z
	go to 10
   15	continue
	a1= tk(j)*(1.+c1*wx(td(j),p(j)))
	a2= tk(i)*(1.+c1*wx(td(i),p(i)))
	z= z+c2*(a1+a2)*(alog(p(i)/pt))
        zthk = z
	return
 20	z= -1.0
        zthk = z
 
	return
	end
c ------------- Shea Routine---------------------
C NCLFORTSTART
      subroutine qcfmThermo (nlev,plev,t,rh,q,iswit,xmsg)
      integer nlev, iswit
      real plev(nlev), t(nlev), rh(nlev), q(nlev), xmsg
C     NCL: qcfm = qcfmThermo (plev[*]:float, t[*]:float, rh[*]:float \
C                            ,q[*]:float, iswit:integer)
C                             xmsg = T@_FillValue
C          note: nlev=dimsizes(plev)
c NCLEND
 
c****         this routine computes q (mixing ratio or specific humidity)
c****         from rh and temp.
 
c****     definition of mixing ratio
c****     if,
c****         es  - is the saturation mixing ratio
c****         ep  - is the ratio of the molecular weights of water
c****               vapor to dry air
c****         p   - is the atmospheric pressure
c****         rh  - is the relative humidity (given as a percent)
 
c****         q =   rh * ( (ep*es)/(p-es) )
 
c****     on input-
c****         nlev   - number of levels
c****         plev   - pressure (mb)
c****         t      - temperature (c)
c****         rh     - rel hum in percent
c****         iswit  - if iswit=0 q will contain mix ratio
c****                  otherwise  q will contaib spc humidity
c****         xmsg   - msg code
 
c****     on output-
c****         q      - mixing ratio or specific hum
 
      data t0,ep,es0,a,b /273.15,.622,6.11,17.269,35.86/
 
      onemep=1.0-ep
 
 
      do 10 k=1,nlev
      if (plev(k).ne.xmsg .and. t(k).ne.xmsg .and. rh(k).ne.xmsg) then
          p  = plev(k)
          tk = t(k)+273.15
          est=es0*exp((a*( tk -t0))/( tk -b))
          qst=(ep*est)/(p-onemep*est)
          q(k)=qst*rh(k)*.01
          if (iswit.ne.0) q(k) = q(k)/(1.+q(k))
c c c     q(k) = q(k)*1000.                    !   kg/kg
      else
          q(k) = xmsg
      endif
   10 continue
 
      return
      end
c ------------- Shea Routine---------------------
C NCLFORTSTART
      real function dewtempThermo (rh,tc)    ! dew point temperature [c]
C     NCL: dewtemp = dewtempThermo(rh:float, tc:float)
c NCLEND
 
      implicit none
      real     rh,tc                   ! input
      real     gc,gcx,tk,tdk,tnot,lhv  ! local
 
c calculate the dew pt temperature [c] given:
c     rh - relative humidty [%]
c     tc - temperature [c]
 
                                       ! gc is gas constant for water vapor
      parameter (gc =461.5)            ! [j/{kg-k}]
      parameter (gcx=gc/(1000.*4.186)) ! [cal/{g-k}]
      parameter (tnot=273.15)
 
 
      tk      = tc+tnot                       ! c ==> k
                                              ! lhv=latent heat vap
      lhv     = ( 597.3-0.57*(tk-273.) )/gcx  ! dutton top p273
      tdk     = tk*lhv/(lhv-tk*alog(rh*0.01)) ! dutton p274 (6) [k]
      dewtemp = tdk-tnot                      ! k ==> c
 
      return
      end
c --- ---------- Shea Routine (test)-------------------
C NCLFORTSTART
      real function rhCCMThermo (t,w,p)  ! relative humidity [%]
C     NCL: rh_ccm = rhCCMThermo(t:float, w:float, p:float)
C NCLEND
c c c implicit none                ! implicit not allowed for ncl

c NCL function: rh = relhum_ccm (t,w,p)
                                   ! input
      real  t                      ! temperature  [K]
      real  w                      ! mixing ratio [kg/kg]
      real  p                      ! pressure     [Pa]  

c This is built from CCM Processor routines:
c     subroutine crlhmp1(w,t,mlon,mlon,klvl,p,rh)
c     subroutine estbl0b(ES,T)
C
C** NAME=ESBL0B
C
C     TABLE OF ES FROM -100 TO +102 C IN ONE-DEGREE INCREMENTS.
C
      real  tabl1(82),tabl2(122)
      real  table(204)
      equivalence (tabl1(1),table(1)),(tabl2(1),table(83))
      data tabl1            /.01403,.01719,.02101,.02561,.03117,.03784,
     1.04584,.05542,.06685,.08049,.09672,.1160,.1388,.1658,.1977,.2353,
     2.2796,.3316,.3925,.4638,.5472,.6444,.7577,.8894,1.042,1.220,1.425,
     31.662,1.936,2.252,2.615,3.032,3.511,4.060,4.688,5.406,6.225,7.159,
     48.223,9.432,10.80,12.36,14.13,16.12,18.38,20.92,23.80,27.03,30.67,
     534.76,39.35,44.49,50.26,56.71,63.93,71.98,80.97,90.98,102.1,114.5,
     6128.3,143.6,160.6,179.4,200.2,223.3,248.8,276.9,307.9,342.1,379.8,
     7421.3,466.9,517.0,572.0,632.3,698.5,770.9,850.2,937.0,1032.0,
     8 1146.6/
      data tabl2                     /1272.0,1408.1,1556.7,1716.9,1890.
     13,2077.6,2279.6,2496.7,2729.8,2980.0,3247.8,3534.1,3839.8,4164.8,
     24510.5,4876.9,5265.1,5675.2,6107.8,6566.2,7054.7,7575.3,8129.4,87
     319.2,9346.5,10013.,10722.,11474.,12272.,13119.,14017.,14969.,1597
     47.,17044.,18173.,19367.,20630.,21964.,23373.,24861.,26430.,28086.
     5,29831.,31671.,33608.,35649.,37796.,40055.,42430.,44927.,47551.,
     650307.,53200.,56236.,59422.,62762.,66264.,69934.,73777.,77802.,82
     7015.,86423.,91034.,95855.,100890.,106160.,111660.,117400.,123400.
     8,129650.,136170.,142980.,150070.,157460.,165160.,173180.,181530.,
     9190220.,199260.,
     X                          208670.,218450.,228610.,239180.,250160.,
     1261560.,273400.,285700.,298450.,311690.,325420.,339650.,354410.,
     2369710.,385560.,401980.,418980.,436590.,454810.,473670.,493170.,
     3513350.,534220.,555800.,578090.,601130.,624940.,649530.,674920.,
     4701130.,728190.,756110.,784920.,814630.,845280.,876880.,909450.,
     5943020.,977610.,1013250.,1049940.,1087740.,1087740./

      real    tp, t2, es, pp 
      integer it

      tp = t
      if(tp .gt. 375.16) tp = 375.16
      if(tp .lt. 173.16) tp = 173.16
      it = tp-173.16
      t2 = 173.16+it
      es =(t2+1.-tp)*table(it+1)+(tp-t2)*table(it+2)
      es = es*1.0e-01

      relhum = (w*(p-0.378*es)/ (0.622*es))*100.

      return
      end
c ---------------From skewt.f ---------------------------
c ---------------This was a subroutine [pvalue]
c ---------------I made this a function with a new name [pvalf] 
C NCLFORTSTART
      real function pvalfThermo (pres,p,np,param)
C Input arguments.
      real p(np),param(np)
      real pres
      integer np
C     NCL: pvalf = pvalfThermo(pres:float, p[*]:float, param[*]:float)
C                  np = dimsizes (p)
c c c implicit none
C ###########################################################################
C Statement of purpose.
C ---------------------
C Compute the value of some parameter in a sounding at a given pressure level
C using log pressure interpolation.
C
C History.
C --------                    
C D. Baker    25 Mar 87     Original version.
C
C Description of input and output.
C --------------------------------
C On input:
C ---------                
C pres        Real          Pressure level where parameter value desired.
C p           Real Array    Sounding pressures (mb).
C np          Integer       Number of pressure levels input.
C param       Real Array    Values of parameter at each pressure level.
C
C On output:
C ----------               
C value       Real          Value of parameter at desired pressure level.
C pvalf replace value
C ###########################################################################
C NCLEND
C Output argument.
      real value
C Local variables.
      integer i
      real p1,p2,p3
C Local constants.
      real flag
      parameter (flag=99999.)
C External functions.
      real interp
C Make sure that the computation is possible.
      value=flag
      if (pres.lt.p(np) .or. pres.gt.p(1)) go to 9999
C Check and see if the lowest pressure is equal to the desired pressure.
      if (abs(p(1)-pres).lt.0.1) then
       value=param(1)
       go to 9999
      end if
C Determine value of parameter at level "pres" in the sounding.  
      do 100 i=2,np
       if (p(i).le.pres) then
        p1=alog(p(i-1))
        p2=alog(pres)
        p3=alog(p(i))
        value=interp(param(i-1),param(i),p1,p2,p3)
        go to 9999
       end if
100   continue
C Return to calling program. Shea added "pvalf=value"
9999  pvalf=value             
      return
      end
c -------------------- From skewt.f; I made this a function-------
C NCLFORTSTART
      REAL FUNCTION SHOWALThermo(P,T,TD,NLVLS)
C     NCL: showal = showalThermo(p[*]:float, t[*]:float, td[*]:float)
c                   nlvls = dimsizes(p)
C C C SUBROUTINE SHOWAL(P,T,TD,NLVLS,SHWLTR)
C C C IMPLICIT NONE
      INTEGER NLVLS
      REAL P(NLVLS),T(NLVLS),TD(NLVLS)
C NCLEND
C
C Statement of purpose.
C ---------------------
C This routine computes the Showalter stability index for a sounding.
C
C History.
C --------                    
C Don Baker      12 May 85    Original version.
C
C Description of input and output.
C --------------------------------
C On input:
C ---------                
C P           Real Array    Sounding pressures (mb).
C T           Real Array    Sounding temperatures (C).
C TD          Real Array    Sounding dew points (C).
C NLVLS       Integer       Number of sounding levels passed.
C
C On output:
C ----------  
C SHWLTR      Real          Showalter stability index (C).
C
C
C Input arguments.
C
C
C Output arguments.
C
      REAL SHWLTR
C
C Internal variables.
C
      REAL P1,P2,P3,T850,TD850,T500,TP500
      REAL TMOIST,PRES
      INTEGER I
C
C External functions.
C
      REAL INTERP,TSA,EPT      
C
C Subroutine constants.
C
      PARAMETER FLAG=99999.
C
C Initialize showalter index to flag.  Exit if surface pressure too low.
C
      SHWLTR=FLAG
      IF (P(1).LE.820.) GO TO 999
C
C Determine 850 mb temperature and dew point from sounding.  If the surface
C is above 900 mb, use 800 mb temp and dew point.
C
      IF (P(1).LT.900.) THEN
         PRES=800.
      ELSE
         PRES=850.
      ENDIF
      DO 1 I=1,NLVLS
         IF (P(I).LE.PRES) THEN
            P1=ALOG(P(I-1))
            P2=ALOG(PRES)
            P3=ALOG(P(I))
            T850=INTERP(T(I-1),T(I),P1,P2,P3)
            TD850=INTERP(TD(I-1),TD(I),P1,P2,P3)
            GO TO 2
         ENDIF
 1    CONTINUE
 2    CONTINUE
C
C Compute the moist adiabat (given by an equivalent potential temperature)
C through the 850 mb based LCL.
C
      TMOIST=EPT(T850,TD850,PRES)
C
C Compute the parcel temperature along this moist adiabat at 500 mb.
C
      TP500=TSA(TMOIST,500.)
C
C Determine 500 mb temperature from the sounding.
C
      PRES=500.
      DO 3 I=1,NLVLS
         IF (P(I).LE.PRES) THEN
            P1=ALOG(P(I-1))
            P2=ALOG(PRES)
            P3=ALOG(P(I))
            T500=INTERP(T(I-1),T(I),P1,P2,P3)
            GO TO 4
         ENDIF
 3    CONTINUE
 4    CONTINUE
C
C Compute the showalter index.
C
      SHWLTR=T500-TP500
      showal=SHWLTR    ! shea
C
C Exit.
C
 999  CONTINUE
      RETURN
      END
c ----------------------------------------------------------
C NCLFORTSTART
      real function interpThermo(y1,y3,x1,x2,x3)
c c c implicit none
C     NCL: interp = interpThermo(y1:float, y3:float, x1:float, x2:float, x3:float)
C NCLEND
C ###########################################################################
C Statement of purpose.
C ---------------------
C Given Y1 AT X1, Y3 AT X3, and X2, calculate Y2 (INTERP) AT X2 via linear
C interpolation.
C
C History.
C --------                    
C D. Baker    01 Jul 84     Original version.
C
C Description of input and output.
C --------------------------------
C On input:
C ---------                
C y1          Real          Value of variable being interpolated at point x1.
C y3          Real          Value of variable being interpolated at point x3.
C x1          Real          Point corresponding to y1.
C x2          Real          Point to which interpolation is desired.
C x3          Real          Point corresponding to y3.
C
C On output:
C ----------               
C interp (y2) Real          Value of interpolated y variable at point x2.
C ###########################################################################
C Input arguments
      real y1,y3,x1,x2,x3
C Perform simple linear interpolation.
      interp=y1+((y3-y1)*((x2-x1)/(x3-x1)))   !interp=y2
C Return to calling routine.
      return
      end
c -------------------------------------------------------------------
C NCLFORTSTART
        real function capeThermo(penv,tenv,nlvl,lclmb,iprnt,tparcel,tmsg
     *                        ,jlcl,jlfc,jcross)  ! <== ncl coords

	integer nlvl,iprnt
	real penv(nlvl)         ! mb [hPa]
        real tenv(nlvl)         ! C
        real lclmb              ! lifting condensation lvl [mb]
        real tparcel(nlvl)      ! returned parcel temp [C]
        real tmsg               ! missing code
c c c   real cape_ncl           ! J
C       NCL: cape = capeThermo(penv[*]:float, tenv[*]:float, lclmb:float\
C                             ,iprnt:integer)
C                   nlvl = dimsizes(penv)
C                tparcel = returned as vector attribute of cape
C                   jlcl = returned as attribute of cape
C                          (integer) subscript of lifting condensation level
C                   jlfc = returned as attribute of cape
C                          (integer) subscript of level of free convection
C                   jlfc = returned as attribute of cape
C                          (integer) subscript where tparcel crosses env sounding
C NCLEND
	parameter(NLMAX=500)	! maximum number of sounding levels
	real tpar(NLMAX)        ! must match "tpar" in cape_djs


        ier   = 0
        iplot = 0
        tnot  = 273.15

        do nl=1,nlvl
           tparcel(nl) = tmsg   ! initialize
        end do

        do nl=1,nlvl
           if (tenv(nl).eq.tmsg) then
               print *,'fortran: cape_ncl: no missing values allowed'
               cape_ncl = tmsg  ! return missing code
               return
           end if
        end do

        do nl=1,nlvl
            tenv(nl) = tenv(nl)+tnot    ! cape_djs requires K
        end do

        do n=1,nlmax
           tpar(n) = tmsg               ! initialize
        end do

	call cape_djs(penv,tenv,nlvl,lclmb,iplot,iprnt,capeval,tpar 
     *                    ,jlcl,jlfc,jcross,ier)
 
        do nl=1,nlvl
           tenv(nl)    = tenv(nl)-tnot   ! return to C
           tparcel(nl) = tpar(nl)-tnot   ! return C
c          write (*,"(' cape_ncl: ',i5,4f10.2)") nl, penv(nl),tenv(nl)
c    *                          ,tparcel(nl)
c    *                          ,(tparcel(nl)-tenv(nl))
        end do
                                         ! return in NCL subscripts
        jlcl   = jlcl-1
        jlfc   = jlfc-1
        jcross = jcross-1

        cape_ncl = capeval
c       print *,"cape_ncl=",cape_ncl

        return 
        end
c ------cape_djs not called directly by NCL--------------------------
cNAME
c	cape - calculates CAPE, given temperature/pressure
c		 data from a vertical sounding
c
cSYNOPSIS
c	call cape_djs(penv,tenv,nlvl,lclmb,iplot,iprnt,capeval,ier)         old
c	call cape_djs(penv,tenv,nlvl,lclmb,iplot,iprnt,capeval,tparcel,ier) change for NCL
c
cDESCRIPTION
c	cape will compute the convective available potential energy
c	given vertical soundings of temperature and pressure.
c	optional plots of environmental temperature, parcel
c	temperature, and cape as a function of altitude are available
c
c	definition of CAPE: the integral wrt -ln(p/p0) of R*(Tparcel-Tenv)	
c		from the Level of Free Convection 'up' to the
c		'crossing level'.
c
cPARAMETERS
c	Parameter	Type	Description
c       
c	penv		(i)	environmental pressure sounding array in mb;
c				the pressure array is assumed to be
c				loaded in a monotonically decreasing sequence.
c
c	tenv		(i)	temperature sounding array in Kelvin at
c				the pressure levels
c
c	nlvl		(i)	number of sounding levels
c
c	lclmb		(i)	approximate pressure at lifting condensation
c				level (mb). the pressure ACTUALLY USED will
c				be the first value on the sounding grid greater
c				than or equal to lclmb.
c
c	iplot		(i)	plotting flag: plots (=1); no plots (=other).
c				the parcel and environmental temperature
c				profiles are plotted throughout the sounding
c				depth; the CAPE profile is plotted only
c				between the LFC and the crossing level.
c
c	iprnt		(i)	printing flag: print (>0); no print (=0). DJS
c
c	capeval		(o)	CAPE at the 'crossing level',
c				defined to be the height where the
c				parcel temp has decreased enough to
c				again equal the environmental temp.
c				if the Level of Free Convection DOES NOT exist
c                               (so that CAPE will be zero or negative)
c				capeval will be assigned ZERO.
c				if the Level of Free Convection DOES exist
c				but the 'crossing level' DOES NOT exist,
c				capeval will be assigned the value of 
c				CAPE at the "top" of the
c				sounding (i.e. smallest pressure)
c
c	tparcel		(o)	temperature of parcel
c
c	ier		(o)	error status:
c				zero on successful return;
c				minus one on error;
c				minus two if 'crossing level' does not
c				exist;
cTECHNICAL NOTES
c	The parcel used to compute CAPE is assumed to initially lie at the
c	level of the largest pressure value in the sounding.
c
cGENERAL NOTES
c	No input parameters are modified by this routine.
c
c	No error checking for missing data is done in this routine.
c	It is the user's responsibility to ensure that no missing data 
c	is present.
c
c       Opening and closing NCAR graphics (via opgks and clsgks) must be
c       performed in the calling program, if iplot = 1
cAUTHOR
c	GLF 9/95

	subroutine cape_djs(penv,tenv,nlvl,lclmb
     *                     ,iplot,iprnt,capeval,tparcel
     *                     ,jlcl,jlfc,jcross,ier)

	implicit none
	integer nlvl,iplot,iprnt,ier,iflag 
	integer NLMAX,ii,jlcl,jlcl1,ndim,jstrt,jcross,jlfc
	real penv(nlvl),tenv(nlvl),lclmb,capeval
	real CENTKEL,G,CP,RD,plcl,tlcl,fret,pmb,mdlnp
	real gammadry,ppar0,tpar0,parcelmb,pelfc

        real tparcel(*)         ! special for NCL return [dimensioned NLMAX]

	parameter(NLMAX=500)	! maximum number of sounding levels
	parameter(CENTKEL=273.15)
	parameter(G=9.81)
	parameter(CP=1004.)
	parameter(RD=287.)	
	
	real zlev,tpar,pe,tempenv,zzlev,ppe
	dimension zlev(NLMAX),tpar(NLMAX),pe(NLMAX),tempenv(NLMAX)
	dimension zzlev(NLMAX),ppe(NLMAX)
	real tmlaps

c	arrays and equivalences for NCAR graphics   [ NOT USED]
        real work,xi,the
        integer idsh
        dimension work(NLMAX,2),idsh(2),xi(5),the(5)
        equivalence (work(1,1),tempenv(1)),(work(1,2),tpar(1))
c       character*16 agdshn
c       data idsh/ O"177777",O"177777"/

	if (nlvl.gt.NLMAX) then
		print *,"cape.f: too many sounding levels"
		ier = -1
		return
	endif

c	determine approximate location of lcl

	jlcl=1
	do ii=2,nlvl
		if (lclmb.le.penv(ii)) then
			jlcl=ii
		endif	
	enddo
        if (iprnt.gt.0) then
            print *,"cape.f: LCL pressure: original ",lclmb," mb"
            print *,"cape.f: adjusted to ",penv(jlcl)," mb" 
            print *,"cape.f: and lies at level ",jlcl
            do ii=1,nlvl
               print *,' data=',penv(ii),tenv(ii)
            end do
        endif
	plcl=penv(jlcl)

c	to compute cape, the parcel is assumed to lie initially at the level
c	of largest pressure

	jstrt=1
	parcelmb=penv(jstrt)
        if (iprnt.gt.0) then
            print *,"cape.f: initial parcel pressure: "
     *             ,         penv(jstrt)," mb"
            print *,"cape.f: parcel temp is ", tenv(jstrt), " K"
            print *,"cape.f: initial P and T lie at level ",jstrt
        endif
	ppar0=penv(jstrt)
	tpar0=tenv(jstrt)

c	assign log-pressure levels for input data,
c	referenced to parcel's initial pressure.

	do ii=1,nlvl
		zlev(ii)=-alog(penv(ii)/ppar0)
	enddo

c	compute parcel temperature at dry lapse rate until the LCL is reached

	gammadry = G/CP
	do ii=1,jlcl
		tpar(ii)=tpar0*exp((RD*gammadry/G)*(-zlev(ii)))
	enddo
	jlcl1=jlcl+1

c	above LCL compute temperature of saturated parcel

	tlcl=tpar(jlcl)-CENTKEL		! convert to C
	the(1)=62.
	xi(1)=1.
	ndim=1

c	linmin determines theta_e

     
	call linmin(the,xi,ndim,fret,plcl,tlcl)
	do ii=jlcl1,nlvl
		pmb=ppar0*exp(-zlev(ii))

c		tmlaps computes temperature given theta_e and pressure

		tpar(ii)=CENTKEL+tmlaps(the(1),pmb)
	enddo
        if (iprnt.gt.0) then
            print * ," "
            write (*,"('cape.f:',16x,'P',6x,'zlev'
     *                          , 6x,'tpar',6x,'tenv')")
            do ii=1,nlvl
               write (*,"(i4,7x,4f10.2)") 
     *                    ii, penv(ii), zlev(ii),tpar(ii),tenv(ii)
            enddo
            print * ," "
        endif

c	determine approximate location of 1) the upper limit of 
c	CAPE integration (the so-called 'crossing level')
c	and 2) the Level of Free Convection

	if (tpar(nlvl).gt.tenv(nlvl)) then
		jcross=nlvl		! no crossing level exists
		ier = -2		! so assign one
		iflag=0
                do ii=nlvl-1,1,-1	! now locate LFC
                        if((tpar(ii).eq.tenv(ii)).and.(iflag.eq.0))then
				jlfc=ii
				iflag=1
			endif
                        if((tpar(ii).lt.tenv(ii)).and.(iflag.eq.0))then
				jlfc=ii+1
				iflag=1
			endif
		enddo
		if (iflag.eq.0) then	! if LFC couldn't be found
			jlfc=1		! assign LFC to surface
		endif	
	else
		iflag=0	
		do ii=nlvl,1,-1		! locate crossing level 
			if((tpar(ii).ge.tenv(ii)).and.(iflag.eq.0))then
				jcross=ii
				iflag=1
			endif
		enddo
		if ((iflag.eq.0).or.(jcross.eq.1)) then	! no crossing level 
			capeval=0.0			! hence no LFC
			ier=0
			return
		endif
		iflag=0
		do ii=jcross-1,1,-1	! crossing level exists so find LFC
                        if((tpar(ii).eq.tenv(ii)).and.(iflag.eq.0))then
				jlfc=ii
				iflag=1
			endif
			if((tpar(ii).lt.tenv(ii)).and.(iflag.eq.0))then
                                jlfc=ii+1
                                iflag=1
                        endif
                enddo
		if (iflag.eq.0) then	! if LFC couldn't be found
			jlfc=1		! assign LFC to surface
		endif	
	endif
        if (iprnt.gt.0) then
            print *,"LFC pressure= ",penv(jlfc)
     &	           ," Crossing pressure= ",penv(jcross)
        endif

c	determine CAPE profile

	pe(1)=0.0
	do ii=2,jcross
 		mdlnp=zlev(ii)-zlev(ii-1)
		pe(ii)=pe(ii-1)+(RD*mdlnp*((tpar(ii)+tpar(ii-1))-
     &		(tenv(ii)+tenv(ii-1)))/2.)
 	enddo
	pelfc=pe(jlfc)
	do ii=1,jcross
		pe(ii)=pe(ii)-pelfc
	enddo

c	save CAPE computed by integrating from LFC to crossing level 

	capeval = pe(jcross)


        do ii=jlfc,jcross
	   zzlev(ii+1-jlfc)=zlev(ii)
	   ppe(ii+1-jlfc)=pe(ii)
	enddo
	do ii=1,nlvl
	   tparcel(ii)=tpar(ii)
        enddo

        if (iprnt.gt.0) then               ! djs added
            write (*,"(/,' jlfc, jcross=',2i5)") jlfc, jcross
            write (*,"(/,'Env Temp & Parcel Temp',/)")   
            do ii=1,jcross
               write (*,"(i4,f7.1,4f9.2,f6.2)") ii 
     *                           ,penv(ii),pe(ii),zlev(ii)
     *                           ,(tempenv(ii)-273.15)
     *                           ,(tpar(ii)-273.15)
     *                           ,(tpar(ii)-tempenv(ii))
c              write (*,"(i4,f7.1,4f9.2,f6.2)") ii 
c    *                           ,penv(ii),pe(ii),zlev(ii)
c    *                           ,(work(ii,1)-273.15)
c    *                           ,(work(ii,2)-273.15)
c    *                           ,(work(ii,2)-work(ii,1))
            enddo
            write (*,"(  '----------------------',/)") 
        endif

	return
	end
c ----------Not called by NCL directly--------------------------
      subroutine linmin(p,xi,n,fret,plcl,tlcl)
c
c  LINMIN: given the n-dim point P and n-dim direction XI, move and
c      reset P to whee the function 'dfd' takes on a minimum along the 
c      direction XI from P, and replaces XI by the actual vector 
c      displacement that P was moved. Also returns as FRET the value
c      of the function at P.
c=================================================================
      save
      parameter(NMAX=50,TOL=1.e-5)
      common /f1com/ pcom(NMAX),xicom(NMAX),plclx,tlclx,ncom
      dimension p(n),xi(n)
      external func
c
c  set common block
c---------------------
      plclx=plcl
      tlclx=tlcl
      ncom=n
      do 100 j=1,n
      pcom(j)=p(j)
      xicom(j)=xi(j)
  100 continue
c
c  find minimum point
c---------------------
      ax=0.0
      xx=1.0
      bx=2.0
      call mnbrak(ax,xx,bx,fa,fx,fb,func)
      fret=brent(ax,xx,bx,func,TOL,xmin)
c
c  construct vector results
c---------------------------
      do 200 j=1,n
      xi(j)=xmin*xi(j)
      p(j)=p(j)+xi(j)
  200 continue
c
c  end routine
c---------------
      return
      end
c ----------Not called by NCL directly--------------------------
      subroutine mnbrak(ax,bx,cx,fa,fb,fc,func)
c
c  MNBRAK: given the subroutine 'func' which returns a function 
c     value, and given distinct initail points AX and BX, this routine
c     searches in the downhill direction (defined by the function as
c     evaluated at the intial points) and returns new points ZX, BX, CX
c     which bracket a minimum of the funtion. Also returned are the
c     function values at the three points, FA, FB, FC.
c===================================================================
      save
ccc   parameter(GOLD=1.618034d0,GLIMIT=100.d0,TINY=1.d-20)
cc    parameter(GOLD=1.6180340,GLIMIT=100.0,TINY=1.e-7)
      parameter(GOLD=1.6180340,GLIMIT=100.0,TINY=1.e-5)
      external func
c
c  get intial function values switch roles of A and B if needed.
c------------------------------------------------------------------
      fa=func(ax)
      fb=func(bx)
      if(fb.gt.fa) then
      dum=ax
      ax=bx
      bx=dum
      dum=fb
      fb=fa
      fa=dum
      endif
c
c  first guess
c-------------
      cx=bx+gold*(bx-ax)
      fc=func(cx)
c
c  loop point for bracketing
c=============================
  100 continue
      if(fb.ge.fc) then
       r=(bx-ax)*(fb-fc)
       q=(bx-cx)*(fb-fa)
       u=bx-((bx-cx)*q-(bx-ax)*r)/(2.0*sign(max(abs(q-r),TINY),q-r))
       ulim=bx+GLIMIT*(cx-bx)
       if((bx-u)*(u-cx).gt.0.0) then
        fu=func(u)
        if(fu.lt.fc) then
         ax=bx
         fa=fb
         bx=u
         fb=fu
         goto 100
        elseif(fu.gt.fb) then
         cx=u
         fc=fu
         goto 100
        endif
        u=cx+GOLD*(cx-bx)
        fu=func(u)
       elseif((cx-u)*(u-ulim).gt.0.0) then
        fu=func(u)
        if(fu.lt.fc) then
         bx=cx
         cx=u
         u=cx+GOLD*(cx-bx)
         fb=fc
         fc=fu
         fu=func(u)
        endif
       elseif((u-ulim)*(ulim-cx).ge.0.0) then
        u=ulim
        fu=func(u)
       else
        u=cx+GOLD*(cx-bx)
        fu=func(u)
       endif
       ax=bx
       bx=cx
       cx=u
       fa=fb
       fb=fc
       fc=fu
       goto 100
      endif
c
c  end routine
c--------------
      return
      end
c ----------Not called by NCL directly--------------------------
      function brent(ax,bx,cx,func,tol,xmin)
c
c  BRENT: given a subroutine 'func' which returns a function value, 
c    and given a bracketing triplet of abscissas AX, BX, CX (such that 
c    BX is between AX and CX, and F(BX) is less than both F(AX) and 
c    F(CX) ), this routine isolates the minimum to a frational precision 
c    of about TOL using brent's method.
c    the abscissa of the minimum is returned as XMIN, and the minimum
c    function value is returned as BRENT.
c======================================================================
      save
ccc   parameter(ITMAX=100,CGOLD=.3819660,ZEPS=1.e-10)
      parameter(ITMAX=100,CGOLD=.3819660,ZEPS=1.e-3)
      external func
c
c  set initial values
c-------------------------
      a=min(ax,cx)
      b=max(ax,cx)
cc    print *,' ax=',ax,' bx=',bx,' cx=',cx,' a=',a,' b=',b
      v=bx
      w=v
      x=v
      e=0.0
      fx=func(x)
      fv=fx
      fw=fx
c
c  loop over the maximum number of iterations
c================================================
      do 300 iter=1,ITMAX
      xm=0.50*(a+b)
      tol1=tol*abs(x)+ZEPS
      tol2=2.0*tol1
      if(abs(x-xm).le.(tol2-.5*(b-a))) goto 400
cc    print *,'|x-xm|=',abs(x-xm),' t2-(b-a)/2=',(tol2-.5e0*(b-a))
      if(abs(e).gt.tol1) then
       r=(x-w)*(fx-fv)
       q=(x-v)*(fx-fw)
       p=(x-v)*q-(x-w)*r
       q=2.*(q-r)
       if(q.gt.0.) p=-p
       q=abs(q)
       etemp=e
       e=d

       if((abs(p).gt.abs(0.5*q*etemp)).or.(p.le.q*(a-x))
     $                                .or.(p.ge.q*(b-x))) goto 100
       d=p/q
       u=x+d
       if(((u-a).lt.tol2).or.((b-u).lt.tol2)) d=sign(tol1,xm-x)
       goto 200
      endif
  100 continue
      if(x.ge.xm) then
      e=a-x
      else
      e=b-x
      endif
      d=CGOLD*e
  200 continue
      if(abs(d).ge.tol1) then
       u=x+d
      else
       u=x+sign(tol1,d)
      endif
      fu=func(u)
      if(fu.le.fx) then
       if(u.ge.x) then
       a=x
       else
       b=x
       endif
       v=w
       fv=fw
       w=x
       fw=fx
       x=u
       fx=fu
      else
       if(u.lt.x) then
       a=u
       else
       b=u
       endif
       if((fu.le.fw).or.(w.eq.x)) then
        v=w
        fv=fw
        w=u
        fw=fu
       elseif((fu.le.fv).or.(v.eq.x).or.(v.eq.w)) then
        v=u
        fv=fu
       endif
      endif
  300 continue
      print *,' BRENT ERROR STOP: Maximum Number of Iterations'
      print *,'      reached, ITMAX=',ITMAX
      stop
c
c  end routine
c-----------------
  400 continue
      xmin=x
      brent=fx
      return
      end
c ----------Not called by NCL directly--------------------------
      function func(x)
c
      save
      parameter(NMAX=50)
      common /f1com/ p(NMAX),xi(NMAX),plcl,tlcl,ncom
      dimension xt(NMAX)
c
      do 100 j=1,ncom
      xt(j)=p(j)+(x*xi(j))
  100 continue
c
c  end routine
c--------------
      func=abs(tlcl-tmlaps(xt,plcl))
      return
      end
