FUNCTION alcl(t,td,p)
! g.s. stipanuk 1973 original version.
! reference stipanuk paper entitled:
! "algorithms for generating a skew-t, log p
! diagram and computing selected meteorological
! quantities."
! atmospheric sciences laboratory
! u.s. army electronics command
! white sands missile range, new mexico 88002
! 33 pages
! baker, schlatter 17-may-1982
! this function returns the pressure alcl (mb) of the lifting conden-
! sation level (lcl) for a parcel initially at temperature t (celsius)
! dew point td (celsius) and pressure p (millibars). alcl is computed
! by an iterative procedure described by eqs. 8-12 in stipanuk (1973),
! pp.13-14.
! determine the mixing ratio line through td and p.
aw = w(td,p)
! determine the dry adiabat through t and p.
ao = o(t,p)
! iterate to locate pressure pi at the intersection of the two
! curves. pi has been set to p for the initial guess.
3 CONTINUE
pi = p
DO i= 1,10
x= .02*(tmr(aw,pi)-tda(ao,pi))
IF (ABS(x) < 0.01) EXIT
pi= pi*(2.**(x))
END DO
alcl= pi
RETURN
ENTRY alclm(t,td,p)
! for entry alclm only, t is the mean potential temperature (celsius)
! and td is the mean mixing ratio (g/kg) of the layer containing the
! parcel.
aw = td
ao = t
GO TO 3
END FUNCTION alcl
FUNCTION ct(wbar,pc,ps)
!
! g.s. stipanuk 1973 original version.
! reference stipanuk paper entitled:
! "algorithms for generating a skew-t, log p
! diagram and computing selected meteorological
! quantities."
! atmospheric sciences laboratory
! u.s. army electronics command
! white sands missile range, new mexico 88002
! 33 pages
! baker, schlatter 17-may-1982
! this function returns the convective temperature ct (celsius)
! given the mean mixing ratio wbar (g/kg) in the surface layer,
! the pressure pc (mb) at the convective condensation level (ccl)
! and the surface pressure ps (mb).
! compute the temperature (celsius) at the ccl.
tc= tmr(wbar,pc)
! compute the potential temperature (celsius), i.e., the dry
! adiabat ao through the ccl.
ao= o(tc,pc)
! compute the surface temperature on the same dry adiabat ao.
ct= tda(ao,ps)
RETURN
END FUNCTION ct
FUNCTION dewpt(ew)
!
! this function yields the dew point dewpt (celsius), given the
! water vapor pressure ew (millibars).
! the empirical formula appears in bolton, david, 1980:
! "the computation of equivalent potential temperature,"
! monthly weather review, vol. 108, no. 7 (july), p. 1047, eq.(11).
! the quoted accuracy is 0.03c or less for -35 < dewpt < 35c.
!
! baker, schlatter 17-may-1982 original version.
enl = ALOG(ew)
dewpt = (243.5*enl-440.8)/(19.48-enl)
RETURN
END FUNCTION dewpt
FUNCTION dpt(ew)
!
! this function returns the dew point dpt (celsius), given the
! water vapor pressure ew (millibars).
!
! baker, schlatter 17-may-1982 original version.
DATA es0/6.1078/
! es0 = saturation vapor pressure (mb) over water at 0c
! return a flag value if the vapor pressure is out of range.
IF (ew > .06.AND.ew < 1013.) GO TO 5
dpt = 9999.
RETURN
5 CONTINUE
! approximate dew point by means of teten's formula.
! the formula appears as eq.(8) in bolton, david, 1980:
! "the computation of equivalent potential temperature,"
! monthly weather review, vol 108, no. 7 (july), p.1047.
! the formula is ew(t) = es0*10**(7.5*t/(t+237.3))
! or ew(t) = es0*exp(17.269388*t/(t+237.3))
! the inverse formula is used below.
x = ALOG(ew/es0)
dnm = 17.269388-x
t = 237.3*x/dnm
fac = 1./(ew*dnm)
! loop for iterative improvement of the estimate of dew point
10 CONTINUE
! get the precise vapor pressure corresponding to t.
edp = esw(t)
! estimate the change in temperature corresponding to (ew-edp)
! assume that the derivative of temperature with respect to
! vapor pressure (dtdew) is given by the derivative of the
! inverse teten formula.
dtdew = (t+237.3)*fac
dt = dtdew*(ew-edp)
t = t+dt
IF (ABS(dt) > 1.e-04) GO TO 10
dpt = t
RETURN
END FUNCTION dpt
FUNCTION dwpt(t,rh)
!
! baker, schlatter 17-may-1982 original version.
!
! this function returns the dew point (celsius) given the temperature
! (celsius) and relative humidity (%). the formula is used in the
! processing of u.s. rawinsonde data and is referenced in parry, h.
! dean, 1969: "the semiautomatic computation of rawinsondes,"
! technical memorandum wbtm edl 10, u.s. department of commerce,
! environmental science services administration, weather bureau,
! office of systems development, equipment development laboratory,
! silver spring, md (october), page 9 and page ii-4, line 460.
x = 1.-0.01*rh
! 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 FUNCTION dwpt
FUNCTION ept(t,td,p)
!
! baker, schlatter 17-may-1982 original version.
!
! this function returns the equivalent potential temperature ept
! (celsius) for a parcel of air initially at temperature t (celsius),
! dew point td (celsius) and pressure p (millibars). the formula used
! is eq.(43) in bolton, david, 1980: "the computation of equivalent
! potential temperature," monthly weather review, vol. 108, no. 7
! (july), pp. 1046-1053. the maximum error in ept in 0.3c. in most
! cases the error is less than 0.1c.
!
! compute the mixing ratio (grams of water vapor per kilogram of
! dry air).
w = wmr(p,td)
! 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 FUNCTION ept
FUNCTION es(t)
!
! this function returns the saturation vapor pressure es (mb) over
! liquid water given the temperature t (celsius). the formula appears
! in bolton, david, 1980: "the computation of equivalent potential
! temperature," monthly weather review, vol. 108, no. 7 (july),
! p. 1047, eq.(10). the quoted accuracy is 0.3% or better for
! -35 < t < 35c.
!
! baker, schlatter 17-may-1982 original version.
!
! es0 = saturation vapor pressure over liquid water at 0c
DATA es0/6.1121/
es = es0*EXP(17.67*t/(t+243.5))
RETURN
END FUNCTION es
FUNCTION esat(t)
!
! g.s. stipanuk 1973 original version.
! reference stipanuk paper entitled:
! "algorithms for generating a skew-t, log p
! diagram and computing selected meteorological
! quantities."
! atmospheric sciences laboratory
! u.s. army electronics command
! white sands missile range, new mexico 88002
! 33 pages
! baker, schlatter 17-may-1982
!
! this function returns the saturation vapor pressure over
! water (mb) given the temperature (celsius).
! the algorithm is due to nordquist, w.s.,1973: "numerical approxima-
! tions of selected meteorlolgical parameters for cloud physics prob-
! lems," ecom-5475, atmospheric sciences laboratory, u.s. army
! electronics command, white sands missile range, new mexico 88002.
tk = t+273.15
p1 = 11.344-0.0303998*tk
p2 = 3.49149-1302.8844/tk
c1 = 23.832241-5.02808*ALOG10(tk)
esat = 10.**(c1-1.3816E-7*10.**p1+8.1328E-3*10.**p2-2949.076/tk)
RETURN
END FUNCTION esat
FUNCTION esgg(t)
!
! baker, schlatter 17-may-1982 original version.
!
! this function returns the saturation vapor pressure over liquid
! water esgg (millibars) given the temperature t (celsius). the
! formula used, due to goff and gratch, appears on p. 350 of the
! smithsonian meteorological tables, sixth revised edition, 1963,
! by roland list.
DATA cta,ews,ts/273.15,1013.246,373.15/
! cta = difference between kelvin and celsius temperatures
! ews = saturation vapor pressure (mb) over liquid water at 100c
! ts = boiling point of water (k)
DATA c1, c2, c3, c4, c5, c6 &
/ 7.90298, 5.02808, 1.3816E-7, 11.344, 8.1328E-3, 3.49149 /
tk = t+cta
! goff-gratch formula
rhs = -c1*(ts/tk-1.)+c2*ALOG10(ts/tk)-c3*(10.**(c4*(1.-tk/ts)) &
-1.)+c5*(10.**(-c6*(ts/tk-1.))-1.)+ALOG10(ews)
esw = 10.**rhs
IF (esw < 0.) esw = 0.
esgg = esw
RETURN
END FUNCTION esgg
FUNCTION esice(t)
!
! baker, schlatter 17-may-1982 original version.
!
! this function returns the saturation vapor pressure with respect to
! ice esice (millibars) given the temperature t (celsius).
! the formula used is based upon the integration of the clausius-
! clapeyron equation by goff and gratch. the formula appears on p.350
! of the smithsonian meteorological tables, sixth revised edition,
! 1963.
DATA cta,eis/273.15,6.1071/
! cta = difference between kelvin and celsius temperature
! eis = saturation vapor pressure (mb) over a water-ice mixture at 0c
DATA c1,c2,c3/9.09718,3.56654,0.876793/
! c1,c2,c3 = empirical coefficients in the goff-gratch formula
IF (t <= 0.) GO TO 5
esice = 99999.
WRITE(6,3)esice
3 FORMAT(' saturation vapor pressure for ice cannot be computed', &
/' for temperature > 0c. esice =',f7.0)
RETURN
5 CONTINUE
! freezing point of water (k)
tf = cta
tk = t+cta
! goff-gratch formula
rhs = -c1*(tf/tk-1.)-c2*ALOG10(tf/tk)+c3*(1.-tk/tf)+ALOG10(eis)
esi = 10.**rhs
IF (esi < 0.) esi = 0.
esice = esi
RETURN
END FUNCTION esice
FUNCTION esilo(t)
!
! baker, schlatter 17-may-1982 original version.
!
! this function returns the saturation vapor pressure over ice
! esilo (millibars) given the temperature t (celsius). the formula
! is due to lowe, paul r., 1977: an approximating polynomial for
! the computation of saturation vapor pressure, journal of applied
! meteorology, vol. 16, no. 1 (january), pp. 100-103.
! the polynomial coefficients are a0 through a6.
DATA a0,a1,a2,a3,a4,a5,a6 &
/6.109177956, 5.034698970E-01, 1.886013408E-02, &
4.176223716E-04, 5.824720280E-06, 4.838803174E-08, &
1.838826904E-10/
IF (t <= 0.) GO TO 5
esilo = 9999.
WRITE(6,3)esilo
3 FORMAT(' saturation vapor pressure over ice is undefined for', &
/' temperature > 0c. esilo =',f6.0)
RETURN
5 CONTINUE
esilo = a0+t*(a1+t*(a2+t*(a3+t*(a4+t*(a5+a6*t)))))
RETURN
END FUNCTION esilo
FUNCTION eslo(t)
!
! baker, schlatter 17-may-1982 original version.
!
! this function returns the saturation vapor pressure over liquid
! water eslo (millibars) given the temperature t (celsius). the
! formula is due to lowe, paul r.,1977: an approximating polynomial
! for the computation of saturation vapor pressure, journal of applied
! meteorology, vol 16, no. 1 (january), pp. 100-103.
! the polynomial coefficients are a0 through a6.
DATA a0,a1,a2,a3,a4,a5,a6 &
/6.107799961, 4.436518521E-01, 1.428945805E-02, &
2.650648471E-04, 3.031240396E-06, 2.034080948E-08, &
6.136820929E-11/
es = a0+t*(a1+t*(a2+t*(a3+t*(a4+t*(a5+a6*t)))))
IF (es < 0.) es = 0.
eslo = es
RETURN
END FUNCTION eslo
FUNCTION esrw(t)
!
! baker, schlatter 17-may-1982 original version.
!
! this function returns the saturation vapor pressure over liquid
! water esrw (millibars) given the temperature t (celsius). the
! formula used is due to richards, j.m., 1971: simple expression
! for the saturation vapour pressure of water in the range -50 to
! 140c, british journal of applied physics, vol. 4, pp.l15-l18.
! the formula was quoted more recently by wigley, t.m.l.,1974:
! comments on 'a simple but accurate formula for the saturation
! vapor pressure over liquid water,' journal of applied meteorology,
! vol. 13, no. 5 (august) p.606.
DATA cta,ts,ews/273.15,373.15,1013.25/
! cta = difference between kelvin and celsius temperature
! ts = temperature of the boiling point of water (k)
! ews = saturation vapor pressure over liquid water at 100c
DATA c1, c2, c3, c4 &
/ 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 < 0) vp = 0.
esrw = vp
RETURN
END FUNCTION esrw
FUNCTION esw(t)
!
! this function returns the saturation vapor pressure esw (millibars)
! over liquid water given the temperature t (celsius). the polynomial
! approximation below is due to herman wobus, a mathematician who
! worked at the navy weather research facility, norfolk, virginia,
! but who is now retired. the coefficients of the polynomial were
! chosen to fit the values in table 94 on pp. 351-353 of the smith-
! sonian meteorological tables by roland list (6th edition). the
! approximation is valid for -50 < t < 100c.
!
! baker, schlatter 17-may-1982 original version.
!
! es0 = saturation vapor ressure over liquid water at 0c
DATA es0/6.1078/
pol = 0.99999683 + t*(-0.90826951E-02 + &
t*(0.78736169E-04 + t*(-0.61117958E-06 + &
t*(0.43884187E-08 + t*(-0.29883885E-10 + &
t*(0.21874425E-12 + t*(-0.17892321E-14 + &
t*(0.11112018E-16 + t*(-0.30994571E-19)))))))))
esw = es0/pol**8
RETURN
END FUNCTION esw
FUNCTION heatl(key,t)
!
! baker, schlatter 17-may-1982 original version.
!
! this function returns the latent heat of
! evaporation/condensation for key=1
! melting/freezing for key=2
! sublimation/deposition for key=3
! for water. the latent heat heatl (joules per kilogram) is a
! function of temperature t (celsius). the formulas are polynomial
! approximations to the values in table 92, p. 343 of the smithsonian
! meteorological tables, sixth revised edition, 1963 by roland list.
! the approximations were developed by eric smith at colorado state
! university.
! 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 == 1) hltnt=a0+a1*tk+a2*tk*tk
IF (key == 2) hltnt=b0+b1*tk+b2*tk*tk
IF (key == 3) hltnt=c0+c1*tk+c2*tk*tk
heatl = hltnt
RETURN
END FUNCTION heatl
FUNCTION hum(t,td)
!
! g.s. stipanuk 1973 original version.
! reference stipanuk paper entitled:
! "algorithms for generating a skew-t, log p
! diagram and computing selected meteorological
! quantities."
! atmospheric sciences laboratory
! u.s. army electronics command
! white sands missile range, new mexico 88002
! 33 pages
! baker, schlatter 17-may-1982
!
! this function returns relative humidity (%) given the
! temperature t and dew point td (celsius). as calculated here,
! relative humidity is the ratio of the actual vapor pressure to
! the saturation vapor pressure.
hum= 100.*(esat(td)/esat(t))
RETURN
END FUNCTION hum
FUNCTION o(t,p)
!
! g.s. stipanuk 1973 original version.
! reference stipanuk paper entitled:
! "algorithms for generating a skew-t, log p
! diagram and computing selected meteorological
! quantities."
! atmospheric sciences laboratory
! u.s. army electronics command
! white sands missile range, new mexico 88002
! 33 pages
! baker, schlatter 17-may-1982
! this function returns potential temperature (celsius) given
! temperature t (celsius) and pressure p (mb) by solving the poisson
! equation.
tk= t+273.15
ok= tk*((1000./p)**.286)
o= ok-273.15
RETURN
END FUNCTION o
FUNCTION oe(t,td,p)
!
! g.s. stipanuk 1973 original version.
! reference stipanuk paper entitled:
! "algorithms for generating a skew-t, log p
! diagram and computing selected meteorological
! quantities."
! atmospheric sciences laboratory
! u.s. army electronics command
! white sands missile range, new mexico 88002
! 33 pages
! baker, schlatter 17-may-1982
! this function returns equivalent potential temperature oe (celsius)
! of a parcel of air given its temperature t (celsius), dew point
! td (celsius) and pressure p (millibars).
! find the wet bulb temperature of the parcel.
atw = tw(t,td,p)
! find the equivalent potential temperature.
oe = os(atw,p)
RETURN
END FUNCTION oe
FUNCTION os(t,p)
!
! g.s. stipanuk 1973 original version.
! reference stipanuk paper entitled:
! "algorithms for generating a skew-t, log p
! diagram and computing selected meteorological
! quantities."
! atmospheric sciences laboratory
! u.s. army electronics command
! white sands missile range, new mexico 88002
! 33 pages
! baker, schlatter 17-may-1982
! this function returns the equivalent potential temperature os
! (celsius) for a parcel of air saturated at temperature t (celsius)
! and pressure p (millibars).
DATA b/2.6518986/
! b is an empirical constant approximately equal to the latent heat
! of vaporization for water divided by the specific heat at constant
! pressure for dry air.
tk = t+273.15
osk= tk*((1000./p)**.286)*(EXP(b*w(t,p)/tk))
os= osk-273.15
RETURN
END FUNCTION os
FUNCTION ow(t,td,p)
!
! g.s. stipanuk 1973 original version.
! reference stipanuk paper entitled:
! "algorithms for generating a skew-t, log p
! diagram and computing selected meteorological
! quantities."
! atmospheric sciences laboratory
! u.s. army electronics command
! white sands missile range, new mexico 88002
! 33 pages
! baker, schlatter 17-may-1982
! this function returns the wet-bulb potential temperature ow
! (celsius) given the temperature t (celsius), dew point td
! (celsius), and pressure p (millibars). the calculation for ow is
! very similar to that for wet bulb temperature. see p.13 stipanuk (1973).
! find the wet bulb temperature of the parcel
atw = tw(t,td,p)
! find the equivalent potential temperature of the parcel.
aos= os(atw,p)
! find the wet-bulb potential temperature.
ow= tsa(aos,1000.)
RETURN
END FUNCTION ow
FUNCTION pccl(pm,p,t,td,wbar,n)
!
! g.s. stipanuk 1973 original version.
! reference stipanuk paper entitled:
! "algorithms for generating a skew-t, log p
! diagram and computing selected meteorological
! quantities."
! atmospheric sciences laboratory
! u.s. army electronics command
! white sands missile range, new mexico 88002
! 33 pages
! baker, schlatter 17-may-1982
!
! this function returns the pressure at the convective condensation
! level given the appropriate sounding data.
! on input:
! p = pressure (millibars). note that p(i).gt.p(i+1).
! t = temperature (celsius)
! td = dew point (celsius)
! n = number of levels in the sounding and the dimension of
! p, t and td
! pm = pressure (millibars) at upper boundary of the layer for
! computing the mean mixing ratio. p(1) is the lower
! boundary.
! on output:
! pccl = pressure (millibars) at the convective condensation level
! wbar = mean mixing ratio (g/kg) in the layer bounded by
! pressures p(1) at the bottom and pm at the top
! the algorithm is decribed on p.17 of stipanuk, g.s.,1973:
! "algorithms for generating a skew-t log p diagram and computing
! selected meteorological quantities," atmospheric sciences labora-
! tory, u.s. army electronics command, white sands missile range, new
! mexico 88002.
DIMENSION t(1),p(1),td(1)
IF (pm /= p(1)) GO TO 5
wbar= w(td(1),p(1))
pc= pm
IF (ABS( GO TO 25
5 CONTINUE
wbar= 0.
k= 0
10 CONTINUE
k = k+1
IF (p(k) > pm) GO TO 10
k= k-1
j= k-1
IF(j < 1) GO TO 20
! compute the average mixing ratio....alog = natural log
DO i= 1,j
l= i+1
wbar= (w(td(i),p(i))+w(td(l),p(l)))*ALOG( +wbar
END DO
20 CONTINUE
l= k+1
! estimate the dew point at pressure pm.
tq= td(k)+(td(l)-td(k))*(ALOG(pm/p(k)))/(ALOG( wbar= wbar+(w(td(k),p(k))+w(tq,pm))*ALOG( wbar= wbar/(2.*ALOG(
! find level at which the mixing ratio line wbar crosses the
! environmental temperature profile.
25 CONTINUE
DO j= 1,n
i= n-j+1
IF (p(i) < 300.) CYCLE
! tmr = temperature (celsius) at pressure p (mb) along a mixing
! ratio line given by wbar (g/kg)
x= tmr(wbar,p(i))-t(i)
IF (x <= 0.) GO TO 35
END DO
pccl= 0.0
RETURN
! 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( DO j= 1,10
del= del/2.
x= tmr(wbar,pc)-t(l)-a*(ALOG(
! the sign function replaces the sign of the first argument
! with that of the second.
pc= pc+SIGN(del,x)
END DO
45 pccl = pc
RETURN
END FUNCTION pccl
FUNCTION pcon(p,t,tc)
!
! this function returns the pressure pcon (mb) at the lifted condensa-
! tion level, given the initial pressure p (mb) and temperature t
! (celsius) of the parcel and the temperature tc (celsius) at the
! lcl. the algorithm is exact. it makes use of the formula for the
! potential temperatures corresponding to t at p and tc at pcon.
! these two potential temperatures are equal.
! baker, schlatter 17-may-1982 original version.
!
DATA akapi/3.5037/
! akapi = (specific heat at constant pressure for dry air) /
! (gas constant for dry air)
! convert t and tc to kelvin temperatures.
tk = t+273.15
tck = tc+273.15
pcon = p*(tck/tk)**akapi
RETURN
END FUNCTION pcon
FUNCTION powt(t,p,td)
!
! this function yields wet-bulb potential temperature powt
! (celsius), given the following input:
! t = temperature (celsius)
! p = pressure (millibars)
! td = dew point (celsius)
!
! baker, schlatter 17-may-1982 original version.
DATA cta,akap/273.15,0.28541/
! cta = difference between kelvin and celsius temperatures
! akap = (gas constant for dry air) / (specific heat at
! constant pressure for dry air)
! compute the potential temperature (celsius)
pt = (t+cta)*(1000./p)**akap-cta
! compute the lifting condensation level (lcl).
tc = tcon(t,td)
! for the origin of the following approximation, see the documen-
! tation for the wobus function.
powt = pt-wobf(pt)+wobf(tc)
RETURN
END FUNCTION powt
FUNCTION precpw(td,p,n)
!
! baker, schlatter 17-may-1982 original version.
!
! this function computes total precipitable water precpw (cm) in a
! vertical column of air based upon sounding data at n levels:
! td = dew point (celsius)
! p = pressure (millibars)
! calculations are done in cgs units.
DIMENSION td(n),p(n)
! g = acceleration due to the earth's gravity (cm/s**2)
DATA g/980.616/
! initialize value of precipitable water
pw = 0.
nl = n-1
! calculate the mixing ratio at the lowest level.
wbot = wmr( DO i=1,nl
wtop = wmr(
! calculate the layer-mean mixing ratio (g/kg).
w = 0.5*(wtop+wbot)
! make the mixing ratio dimensionless.
wl = .001*w
! calculate the specific humidity.
ql = wl/(wl+1.)
! 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
END DO
precpw = pw
RETURN
END FUNCTION precpw
SUBROUTINE ptlcl(p,t,td,pc,tc) 2
!
! this subroutine estimates the pressure pc (mb) and the temperature
! tc (celsius) at the lifted condensation level (lcl), given the
! initial pressure p (mb), temperature t (celsius) and dew point
! (celsius) of the parcel. the approximation is that lines of
! constant potential temperature and constant mixing ratio are
! straight on the skew t/log p chart.
!
! baker,schlatter 17-may-1982 original version
!
! teten's formula for saturation vapor pressure as a function of
! pressure was used in the derivation of the formula below. for
! additional details, see math notes by t. schlatter dated 8 sep 81.
! t. schlatter, noaa/erl/profs program office, boulder, colorado,
! wrote this subroutine.
!
! akap = (gas constant for dry air) / (specific heat at constant
! pressure for dry air)
! cta = difference between kelvin and celsius temperatures
!
DATA akap,cta/0.28541,273.16/
c1 = 4098.026/(td+237.3)**2
c2 = 1./(akap*(t+cta))
pc = p*EXP(c1*c2*(t-td)/(c2-c1))
tc = t+c1*(t-td)/(c2-c1)
RETURN
END SUBROUTINE ptlcl
FUNCTION satlft(thw,p)
!
! baker, schlatter 17-may-1982 original version.
!
! input: thw = wet-bulb potential temperature (celsius).
! thw defines a moist adiabat.
! p = pressure (millibars)
! output: satlft = temperature (celsius) where the moist adiabat
! crosses p
DATA cta,akap/273.15,0.28541/
! cta = difference between kelvin and celsius temperatures
! akap = (gas constant for dry air) / (specific heat at constant
! pressure for dry air)
! the algorithm below can best be understood by referring to a
! skew-t/log p chart. it was devised by herman wobus, a mathemati-
! cian formerly at the navy weather research facility but now retired.
! the value returned by satlft can be checked by referring to table
! 78, pp.319-322, smithsonian meteorological tables, by roland list
! (6th revised edition).
!
IF (p /= 1000.) GO TO 5
satlft = thw
RETURN
5 CONTINUE
! compute tone, the temperature where the dry adiabat with value thw
! (celsius) crosses p.
pwrp = (p/1000.)**akap
tone = (thw+cta)*pwrp-cta
! consider the moist adiabat ew1 through tone at p. using the defini-
! tion of the wobus function (see documentation on wobf), it can be
! shown that eone = ew1-thw.
eone = wobf(tone)-wobf(thw)
rate = 1.
GO TO 15
! in the loop below, the estimate of satlft is iteratively improved.
10 CONTINUE
! rate is the ratio of a change in t to the corresponding change in
! e. its initial value was set to 1 above.
rate = (ttwo-tone)/(etwo-eone)
tone = ttwo
eone = etwo
15 CONTINUE
! ttwo is an improved estimate of satlft.
ttwo = tone-eone*rate
! pt is the potential temperature (celsius) corresponding to ttwo at p
pt = (ttwo+cta)/pwrp-cta
! consider the moist adiabat ew2 through ttwo at p. using the defini-
! tion of the wobus function, it can be shown that etwo = ew2-thw.
etwo = pt+wobf(ttwo)-wobf(pt)-thw
! dlt is the correction to be subtracted from ttwo.
dlt = etwo*rate
IF (ABS(dlt) > 0.1) GO TO 10
satlft = ttwo-dlt
RETURN
END FUNCTION satlft
FUNCTION ssh(p,t)
!
! baker, schlatter 17-may-1982 original version.
!
! this function returns saturation specific humidity ssh (grams of
! water vapor per kilogram of moist air) given the pressure p
! (millibars) and the temperature t (celsius). the equation is given
! in standard meteorological texts. if t is dew point (celsius), then
! ssh returns the actual specific humidity.
! compute the dimensionless mixing ratio.
w = .001*wmr(p,t)
! compute the dimensionless saturation specific humidity.
q = w/(1.+w)
ssh = 1000.*q
RETURN
END FUNCTION ssh
FUNCTION tcon(t,d)
!
! this function returns the temperature tcon (celsius) at the lifting
! condensation level, given the temperature t (celsius) and the
! dew point d (celsius).
!
! baker, schlatter 17-may-1982 original version.
!
! compute the dew point depression s.
s = t-d
! the approximation below, a third order polynomial in s and t,
! is due to herman wobus. the source of data for fitting the
! polynomial is unknown.
dlt = s*(1.2185+1.278E-03*t+ &
s*(-2.19E-03+1.173E-05*s-5.2E-06*t))
tcon = t-dlt
RETURN
END FUNCTION tcon
FUNCTION tda(o,p)
!
! g.s. stipanuk 1973 original version.
! reference stipanuk paper entitled:
! "algorithms for generating a skew-t, log p
! diagram and computing selected meteorological
! quantities."
! atmospheric sciences laboratory
! u.s. army electronics command
! white sands missile range, new mexico 88002
! 33 pages
! baker, schlatter 17-may-1982
! this function returns the temperature tda (celsius) on a dry adiabat
! at pressure p (millibars). the dry adiabat is given by
! potential temperature o (celsius). the computation is based on
! poisson's equation.
ok= o+273.15
tdak= ok*((p*.001)**.286)
tda= tdak-273.15
RETURN
END FUNCTION tda
FUNCTION te(t,td,p)
!
! g.s. stipanuk 1973 original version.
! reference stipanuk paper entitled:
! "algorithms for generating a skew-t, log p
! diagram and computing selected meteorological
! quantities."
! atmospheric sciences laboratory
! u.s. army electronics command
! white sands missile range, new mexico 88002
! 33 pages
! baker, schlatter 17-may-1982
! this function returns the equivalent temperature te (celsius) of a
! parcel of air given its temperature t (celsius), dew point (celsius)
! and pressure p (millibars).
! calculate equivalent potential temperature.
aoe = oe(t,td,p)
! use poissons's equation to calculate equivalent temperature.
te= tda(aoe,p)
RETURN
END FUNCTION te
FUNCTION thm(t,p)
!
! baker, schlatter 17-may-1982 original version.
!
! this function returns the wet-bulb potential temperature thm
! (celsius) corresponding to a parcel of air saturated at
! temperature t (celsius) and pressure p (millibars).
f(x) = 1.8199427E+01+x*( 2.1640800E-01+x*( 3.0716310E-04+x* &
(-3.8953660E-06+x*( 1.9618200E-08+x*( 5.2935570E-11+x* &
( 7.3995950E-14+x*(-4.1983500E-17)))))))
thm = t
IF (p == 1000.) RETURN
! 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 FUNCTION thm
FUNCTION tlcl(t,td)
!
! this function yields the temperature tlcl (celsius) of the lifting
! condensation level, given the temperature t (celsius) and the
! dew point td (celsius).
!
! baker,schlatter 17-may-1982 original version
!
! the formula used is n bolton, david,
! 1980: "the computation of equivalent potential temperature,"
! monthly weather review, vol. 108, no. 7 (july), p. 1048, eq.(15).
!
! convert from celsius to kelvin degrees.
tk = t+273.16
tdk = td+273.16
a = 1./(tdk-56.)
b = ALOG(tk/tdk)/800.
tc = 1./(a+b)+56.
tlcl = tc-273.16
RETURN
END FUNCTION tlcl
FUNCTION tlcl1(t,td)
!
! baker, schlatter 17-may-1982 original version.
!
! this function returns the temperature tlcl1 (celsius) of the lifting
! condensation level (lcl) given the initial temperature t (celsius)
! and dew point td (celsius) of a parcel of air.
! eric smith at colorado state university has used the formula
! below, but its origin is unknown.
DATA cta/273.15/
! cta = difference between kelvin and celsius temperature
tk = t+cta
! 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 FUNCTION tlcl1
FUNCTION tmlaps(thetae,p)
!
! baker, schlatter 17-may-1982 original version.
!
! this function returns the temperature tmlaps (celsius) at pressure
! p (millibars) along the moist adiabat corresponding to an equivalent
! potential temperature thetae (celsius).
! the algorithm was written by eric smith at colorado state
! university.
DATA crit/0.1/
! cta = difference between kelvin and celsius temperatures.
! crit = convergence criterion (degrees kelvin)
eq0 = thetae
! initial guess for solution
tlev = 25.
! compute the saturation equivalent potential temperature correspon-
! ding to temperature tlev and pressure p.
eq1 = ept(tlev,tlev,p)
dif = ABS(eq1-eq0)
IF (dif < crit) GO TO 3
IF (eq1 > eq0) GO TO 1
! 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 < crit) GO TO 3
j = -1
IF (eq1 > eq0) j=1
IF (i == j) GO TO 2
! the solution has been passed. reverse the direction of search
! and decrease the stepping increment.
tlev = tlev-dt
dt = dt/10.
GO TO 2
3 tmlaps = tlev
RETURN
END FUNCTION tmlaps
FUNCTION tmr(w,p)
!
! g.s. stipanuk 1973 original version.
! reference stipanuk paper entitled:
! "algorithms for generating a skew-t, log p
! diagram and computing selected meteorological
! quantities."
! atmospheric sciences laboratory
! u.s. army electronics command
! white sands missile range, new mexico 88002
! 33 pages
! baker, schlatter 17-may-1982
! this function returns the temperature (celsius) on a mixing
! ratio line w (g/kg) at pressure p (mb). the formula is given in
! table 1 on page 7 of stipanuk (1973).
!
! initialize constants
DATA c1/.0498646455/,c2/2.4082965/,c3/7.07475/
DATA c4/38.9114/,c5/.0915/,c6/1.2035/
x= ALOG10(w*p/(622.+w))
tmrk= 10.**(c1*x+c2)-c3+c4*((10.**(c5*x)-c6)**2.)
tmr= tmrk-273.15
RETURN
END FUNCTION tmr
FUNCTION tsa(os,p)
!
! g.s. stipanuk 1973 original version.
! reference stipanuk paper entitled:
! "algorithms for generating a skew-t, log p
! diagram and computing selected meteorological
! quantities."
! atmospheric sciences laboratory
! u.s. army electronics command
! white sands missile range, new mexico 88002
! 33 pages
! baker, schlatter 17-may-1982
! this function returns the temperature tsa (celsius) on a saturation
! adiabat at pressure p (millibars). os is the equivalent potential
! temperature of the parcel (celsius). sign(a,b) replaces the
! algebraic sign of a with that of b.
! b is an empirical constant approximately equal to 0.001 of the latent
! heat of vaporization for water divided by the specific heat at constant
! pressure for dry air.
DATA b/2.6518986/
a= os+273.15
! tq is the first guess for tsa.
tq= 253.15
! d is an initial value used in the iteration below.
d= 120.
! iterate to obtain sufficient accuracy....see table 1, p.8
! of stipanuk (1973) for equation used in iteration.
DO i= 1,12
tqk= tq-273.15
d= d/2.
x= a*EXP(-b*w(tqk,p)/tq)-tq*((1000./p)**.286)
IF (ABS(x) < 1E-7) GOTO 2
tq= tq+SIGN(d,x)
END DO
2 tsa= tq-273.15
RETURN
END FUNCTION tsa
FUNCTION tv(t,td,p)
!
! baker, schlatter 17-may-1982 original version.
!
! this function returns the virtual temperature tv (celsius) of
! a parcel of air at temperature t (celsius), dew point td
! (celsius), and pressure p (millibars). the equation appears
! in most standard meteorological texts.
DATA cta,eps/273.15,0.62197/
! cta = difference between kelvin and celsius temperatures.
! eps = ratio of the mean molecular weight of water (18.016 g/mole)
! to that of dry air (28.966 g/mole)
tk = t+cta
! calculate the dimensionless mixing ratio.
w = .001*wmr(p,td)
tv = tk*(1.+w/eps)/(1.+w)-cta
RETURN
END FUNCTION tv
FUNCTION tw(t,td,p)
! g.s. stipanuk 1973 original version.
! reference stipanuk paper entitled:
! "algorithms for generating a skew-t, log p
! diagram and computing selected meteorological
! quantities."
! atmospheric sciences laboratory
! u.s. army electronics command
! white sands missile range, new mexico 88002
! 33 pages
! baker, schlatter 17-may-1982
! this function returns the wet-bulb temperature tw (celsius)
! given the temperature t (celsius), dew point td (celsius)
! and pressure p (mb). see p.13 in stipanuk (1973), referenced
! above, for a description of the technique.
!
!
! determine the mixing ratio line thru td and p.
aw = w(td,p)
!
! determine the dry adiabat thru t and p.
ao = o(t,p)
pi = p
! iterate to locate pressure pi at the intersection of the two
! curves . pi has been set to p for the initial guess.
DO i= 1,10
x= .02*(tmr(aw,pi)-tda(ao,pi))
IF (ABS(x) < 0.01) EXIT
pi= pi*(2.**(x))
END DO
! find the temperature on the dry adiabat ao at pressure pi.
ti= tda(ao,pi)
! the intersection has been located...now, find a saturation
! adiabat thru this point. function os returns the equivalent
! potential temperature (c) of a parcel saturated at temperature
! ti and pressure pi.
aos= os(ti,pi)
! function tsa returns the wet-bulb temperature (c) of a parcel at
! pressure p whose equivalent potential temperature is aos.
tw = tsa(aos,p)
RETURN
END FUNCTION tw
FUNCTION w(t,p)
!
! g.s. stipanuk 1973 original version.
! reference stipanuk paper entitled:
! "algorithms for generating a skew-t, log p
! diagram and computing selected meteorological
! quantities."
! atmospheric sciences laboratory
! u.s. army electronics command
! white sands missile range, new mexico 88002
! 33 pages
! baker, schlatter 17-may-1982
!
! this function returns the mixing ratio (grams of water vapor per
! kilogram of dry air) given the dew point (celsius) and pressure
! (millibars). if the temperture is input instead of the
! dew point, then saturation mixing ratio (same units) is returned.
! the formula is found in most meteorological texts.
x= esat(t)
w= 622.*x/(p-x)
RETURN
END FUNCTION w
FUNCTION wmr(p,t)
!
! this function approximates the mixing ratio wmr (grams of water
! vapor per kilogram of dry air) given the pressure p (mb) and the
! temperature t (celsius). the formula used is given on p. 302 of the
! smithsonian meteorological tables by roland list (6th edition).
!
! baker, schlatter 17-may-1982 original version.
!
! eps = ratio of the mean molecular weight of water (18.016 g/mole)
! to that of dry air (28.966 g/mole)
DATA eps/0.62197/
! the next two lines contain a formula by herman wobus for the
! correction factor wfw for the departure of the mixture of air
! and water vapor from the ideal gas law. the formula fits values
! in table 89, p. 340 of the smithsonian meteorological tables,
! but only for temperatures and pressures normally encountered in
! 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)
! convert r from a dimensionless ratio to grams/kilogram.
wmr = 1000.*r
RETURN
END FUNCTION wmr
FUNCTION wobf(t)
! this function calculates the difference of the wet-bulb potential
! temperatures for saturated and dry air given the temperature.
!
! baker, schlatter 17-may-1982 original version.
!
! let wbpts = wet-bulb potential temperature for saturated
! air at temperature t (celsius). let wbptd = wet-bulb potential
! temperature for completely dry air at the same temperature t.
! the wobus function wobf (in degrees celsius) is defined by
! wobf(t) = wbpts-wbptd.
! although wbpts and wbptd are functions of both pressure and
! temperature, their difference is a function of temperature only.
! to understand why, consider a parcel of dry air at tempera-
! ture t and pressure p. the thermodynamic state of the parcel is
! represented by a point on a pseudoadiabatic chart. the wet-bulb
! potential temperature curve (moist adiabat) passing through this
! point is wbpts. now t is the equivalent temperature for another
! parcel saturated at some lower temperature tw, but at the same
! pressure p. to find tw, ascend along the dry adiabat through
! (t,p). at a great height, the dry adiabat and some moist
! adiabat will nearly coincide. descend along this moist adiabat
! back to p. the parcel temperature is now tw. the wet-bulb
! potential temperature curve (moist adiabat) through (tw,p) is wbptd.
! the difference (wbpts-wbptd) is proportional to the heat imparted
! to a parcel saturated at temperature tw if all its water vapor
! were condensed. since the amount of water vapor a parcel can
! hold depends upon temperature alone, (wbptd-wbpts) must depend
! on temperature alone.
! the wobus function is useful for evaluating several thermo-
! dynamic quantities. by definition:
! wobf(t) = wbpts-wbptd. (1)
! if t is at 1000 mb, then t is a potential temperature pt and
! wbpts = pt. thus
! wobf(pt) = pt-wbptd. (2)
! if t is at the condensation level, then t is the condensation
! temperature tc and wbpts is the wet-bulb potential temperature
! wbpt. thus
! wobf(tc) = wbpt-wbptd. (3)
! if wbptd is eliminated from (2) and (3), there results
! wbpt = pt-wobf(pt)+wobf(tc).
! if wbptd is eliminated from (1) and (2), there results
! wbpts = pt-wobf(pt)+wobf(t).
! if t is an equivalent potential temperature ept (implying
! that the air at 1000 mb is completely dry), then wbpts = ept
! and wbptd = wbpt. thus
! wobf(ept) = ept-wbpt.
! this form is the basis for a polynomial approximation to wobf.
! in table 78 on pp.319-322 of the smithsonian meteorological
! tables by roland list (6th revised edition), one finds wet-bulb
! potential temperatures and the corresponding equivalent potential
! temperatures listed together. herman wobus, a mathematician for-
! merly at the navy weather research facility, norfolk, virginia,
! and now retired, computed the coefficients for the polynomial
! approximation from numbers in this table.
!
! notes by t.w. schlatter
! noaa/erl/profs program office
! august 1981
x = t-20.
IF (x > 0.) GO TO 10
pol = 1. +x*(-8.8416605E-03 &
+x*( 1.4714143E-04 +x*(-9.6719890E-07 &
+x*(-3.2607217E-08 +x*(-3.8598073E-10)))))
wobf = 15.130/pol**4
RETURN
10 CONTINUE
pol = 1. +x*( 3.6182989E-03 &
+x*(-1.3603273E-05 +x*( 4.9618922E-07 &
+x*(-6.1059365E-09 +x*( 3.9401551E-11 &
+x*(-1.2588129E-13 +x*( 1.6688280E-16)))))))
wobf = 29.930/pol**4+0.96*x-14.8
RETURN
END FUNCTION wobf
FUNCTION z(pt,p,t,td,n)
!
! g.s. stipanuk 1973 original version.
! reference stipanuk paper entitled:
! "algorithms for generating a skew-t, log p
! diagram and computing selected meteorological
! quantities."
! atmospheric sciences laboratory
! u.s. army electronics command
! white sands missile range, new mexico 88002
! 33 pages
! baker, schlatter 17-may-1982
!
! this function returns the thickness of a layer bounded by pressure
! p(1) at the bottom and pressure pt at the top.
! on input:
! p = pressure (mb). note that p(i).gt.p(i+1).
! t = temperature (celsius)
! td = dew point (celsius)
! n = number of levels in the sounding and the dimension of
! p, t and td
! on output:
! z = geometric thickness of the layer (m)
! the algorithm involves numerical integration of the hydrostatic
! equation from p(1) to pt. it is described on p.15 of stipanuk
! (1973).
DIMENSION t(1),p(1),td(1),tk(100)
! c1 = .001*(1./eps-1.) where eps = .62197 is the ratio of the
! molecular weight of water to that of
! dry air. the factor 1000. converts the
! mixing ratio w from g/kg to a dimension-
! less ratio.
! c2 = r/(2.*g) where r is the gas constant for dry air
! (287 kg/joule/deg k) and g is the acceleration
! due to the earth's gravity (9.8 m/s**2). the
! factor of 2 is used in averaging two virtual
! temperatures.
DATA c1/.0006078/,c2/14.64285/
DO i= 1,n
tk(i)= t(i)+273.15
END DO
z= 0.0
IF (pt < p(n)) GO TO 20
i= 0
10 i= i+1
j= i+1
IF (pt >= p(j)) GO TO 15
a1= tk(j)*(1.+c1*w(td(j),p(j)))
a2= tk(i)*(1.+c1*w(td(i),p(i)))
z= z+c2*(a1+a2)*(ALOG( GO TO 10
15 CONTINUE
a1= tk(j)*(1.+c1*w(td(j),p(j)))
a2= tk(i)*(1.+c1*w(td(i),p(i)))
z= z+c2*(a1+a2)*(ALOG( RETURN
20 z= -1.0
RETURN
END FUNCTION z