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