Kerry Emanuel provides FORTRAN 77 code to accompanying his book
'Atmospheric Convection'.
The code computes CAPE in several different ways - and the theory is
described in the book.
I have found this useful and enlightening.
The code can be obtained - free - off the page:
http://wind.mit.edu/~emanuel/ftpacc.html
The fortran would be easier to wrap for ncl than translating idl.
I have wrapped in for python with no problems.
--Jim
On Oct 12, 2004, at 9:30 AM, Christian Pagé wrote:
> Hi Dennis,
>
> There is an IDL subroutine available here :
> http://www.iac.ethz.ch/staff/dominik/idltools/atmos_phys/cape_sound.pro
> based on fortran code, but greatly enhanced per the comments.
>
> Christian
>
> On Tue, 12 Oct 2004 10:10:13 -0600 (MDT), Dennis Shea
> <shea@cgd.ucar.edu> wrote:
>> Hi Christian
>>
>> off-line
>>
>> The "cape_thermo" is one of several included but not supported
>> functions. That is why there is no description in the NCL
>> functions and procedures page.
>>
>> I created skewt as a way of learning NCL.
>> It was for my own use. Then others saw the plots
>> and it became a defacto standard in the NCL distribution.
>> It has not changed since 1997.
>>
>> I included several test codes. One was cape thermo.
>> I am attaching it. It was written by a person at the U of Colo.
>> It was/is not supported because I could not verify that the
>> returned "cape" value was correct. Perhaps, the way the LCL
>> is used/computed is a reason why the cape valu never
>> verified with other code.
>>
>> The LCL as calculated by the PROFS function
>> http://ngwww.ucar.edu/ngdoc/ng/ref/ncl/functions/lclvl.html
>> is the only supported way to calculate the LCL.
>>
>> Do u have a verifiable "cape" code [C or fortran]?
>>
>> THX
>> D
>>
>>
>>>
>>> Hi,
>>>
>>> I have a question. What is the use of the subroutine
>>> cape = cape_thermo(p,tc,plcl,iprnt)
>>> in skewt_func.ncl ?
>>> It recomputes the plcl, but in a weird way in my sense.
>>> In my case, it recomputes the LCL to be at the first level, which is
>>> not saturated.
>>> The saturation occurs at 951.9 mb, not at 994 mb. Why does this
>>> subroutine recomputes the LCL pressure??
>>>
>>> Many thanks
>>>
>>> cape.f: LCL pressure: original 951.8922729492188 mb
>>> cape.f: adjusted to 994.0000000000000 mb
>>> cape.f: and lies at level 1
>>> data= 994.0000000000000 287.7500003814697
>>> data= 925.0000000000000 282.7500003814697
>>> data= 907.0000000000000 281.3799995422363
>>> data= 866.0000000000000 278.1500000000000
>>> data= 850.0000000000000 277.5500000953674
>>> data= 811.0000000000000 275.7499999046325
>>> data= 753.0000000000000 273.1500000000000
>>> data= 725.0000000000000 271.0500000953674
>>> data= 700.0000000000000 269.2499999046325
>>> data= 673.0000000000000 267.2499999046325
>>> data= 646.0000000000000 265.0499996185303
>>> data= 596.0000000000000 260.6500000000000
>>> data= 582.0000000000000 260.8499998092651
>>> data= 500.0000000000000 253.6500000000000
>>> data= 492.0000000000000 252.8500007629394
>>> data= 458.0000000000000 248.4499992370605
>>> data= 446.0000000000000 246.7600006103515
>>> data= 400.0000000000000 239.8500007629394
>>> data= 359.0000000000000 233.2499984741211
>>> data= 304.0000000000000 224.2499984741211
>>> data= 300.0000000000000 223.6500000000000
>>> data= 282.0000000000000 220.9999984741211
>>> data= 250.0000000000000 215.8500007629394
>>> data= 245.0000000000000 215.0500015258789
>>> data= 223.0000000000000 216.4400009155273
>>> data= 200.0000000000000 218.0500015258789
>>> data= 167.0000000000000 220.8500007629394
>>> data= 163.0000000000000 220.4000000000000
>>> data= 150.0000000000000 218.8500007629394
>>> data= 134.0000000000000 216.8500007629394
>>> data= 107.0000000000000 218.0500015258789
>>> data= 100.0000000000000 217.8500007629394
>>> cape.f: initial parcel pressure: 994.0000000000000 mb
>>> cape.f: parcel temp is 287.7500003814697 K
>>> cape.f: initial P and T lie at level 1
>>>
>>> cape.f: P zlev tpar tenv
>>> 1 994.00 0.00 287.75 287.75
>>> 2 925.00 0.07 284.82 282.75
>>> 3 907.00 0.09 283.95 281.38
>>> 4 866.00 0.14 282.05 278.15
>>> 5 850.00 0.16 281.25 277.55
>>> 6 811.00 0.20 279.25 275.75
>>> 7 753.00 0.28 275.95 273.15
>>> 8 725.00 0.32 274.25 271.05
>>> 9 700.00 0.35 272.65 269.25
>>> 10 673.00 0.39 270.75 267.25
>>> 11 646.00 0.43 268.85 265.05
>>> 12 596.00 0.51 264.85 260.65
>>> 13 582.00 0.54 263.65 260.85
>>> 14 500.00 0.69 255.65 253.65
>>> 15 492.00 0.70 254.75 252.85
>>> 16 458.00 0.77 250.65 248.45
>>> 17 446.00 0.80 249.15 246.76
>>> 18 400.00 0.91 242.75 239.85
>>> 19 359.00 1.02 236.15 233.25
>>> 20 304.00 1.18 225.85 224.25
>>> 21 300.00 1.20 225.05 223.65
>>> 22 282.00 1.26 221.25 221.00
>>> 23 250.00 1.38 213.95 215.85
>>> 24 245.00 1.40 212.75 215.05
>>> 25 223.00 1.49 207.15 216.44
>>> 26 200.00 1.60 200.85 218.05
>>> 27 167.00 1.78 190.75 220.85
>>> 28 163.00 1.81 189.45 220.40
>>> 29 150.00 1.89 185.05 218.85
>>> 30 134.00 2.00 179.15 216.85
>>> 31 107.00 2.23 168.05 218.05
>>> 32 100.00 2.30 164.85 217.85
>>>
>>> LFC pressure= 994.0000000000000 Crossing pressure=
>>> 282.0000000000000
>>>
>>>
>>> jlfc, jcross= 1 22
>>>
>>> Env Temp & Parcel Temp
>>>
>>> 1 994.0 0.00 0.00 -273.15 14.60287.75
>>> 2 925.0 21.37 0.07 -273.15 11.67284.82
>>> 3 907.0 34.46 0.09 -273.15 10.80283.95
>>> 4 866.0 77.40 0.14 -273.15 8.90282.05
>>> 5 850.0 97.74 0.16 -273.15 8.10281.25
>>> 6 811.0 146.27 0.20 -273.15 6.10279.25
>>> 7 753.0 213.35 0.28 -273.15 2.80275.95
>>> 8 725.0 245.98 0.32 -273.15 1.10274.25
>>> 9 700.0 279.21 0.35 -273.15 -0.50272.65
>>> 10 673.0 318.16 0.39 -273.15 -2.40270.75
>>> 11 646.0 361.05 0.43 -273.15 -4.30268.85
>>> 12 596.0 453.53 0.51 -273.15 -8.30264.85
>>> 13 582.0 477.41 0.54 -273.15 -9.50263.65
>>> 14 500.0 582.01 0.69 -273.15 -17.50255.65
>>> 15 492.0 591.04 0.70 -273.15 -18.40254.75
>>> 16 458.0 633.17 0.77 -273.15 -22.50250.65
>>> 17 446.0 650.66 0.80 -273.15 -24.00249.15
>>> 18 400.0 733.29 0.91 -273.15 -30.40242.75
>>> 19 359.0 823.30 1.02 -273.15 -37.00236.15
>>> 20 304.0 930.69 1.18 -273.15 -47.30225.85
>>> 21 300.0 936.39 1.20 -273.15 -48.10225.05
>>> 22 282.0 951.04 1.26 -273.15 -51.90221.25
>>>
>>>
>>> --
>>> Christian Pagé
>>> http://meteocentre.com/ http://meteoalerte.com/
>>> Etudiant au Doctorat en Sciences de l'environnement UQAM
>>>
>>> _______________________________________________
>>> ncl-talk mailing list
>>> ncl-talk@ucar.edu
>>> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
>>
>>
>> c -------------------------------------------------------------------
>> C NCLFORTSTART
>> real function
>> cape_ncl(penv,tenv,nlvl,lclmb,iprnt,tparcel,tmsg)
>>
>> 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 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 -------------------------------------------------------------------
>> 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))
>> 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
>>
>> c CAPE
>>
>> c plots
>> c if (iplot.eq.1) then
>> c call agseti('INVERT.',1)
>> c call agsetf('GRAPH/LEFT.',0.1)
>> c call agsetf('GRAPH/RIGHT.',0.9)
>> c call agsetf('GRAPH/BOTTOM.',0.0)
>> c call agsetf('GRAPH/TOP.',1.0)
>> c call agseti('RIGHT/FUNCTION.',1)
>> c call agsetf('RIGHT/NUMERIC/TYPE.',1.e36)
>> c call ezxy(zzlev,ppe,(jcross-jlfc)+1,'CAPE (joules)')
>>
>> c environmental temp and parcel temp
>>
>> c call agseti('TOP/FUNCTION.',2)
>> c call agsetf('TOP/NUMERIC/TYPE.',1.e36)
>> c call agseti('DASH/SELECTOR.',3)
>> c call agseti('DASH/LENGTH.',8)
>> c call agseti(agdshn(1),idsh(1))
>> c call agseti(agdshn(2),idsh(2))
>> c call ezmxy(zlev,work,NLMAX,2,jcross
>> c * ,'Env Temp & Parcel Temp')
>> c endif
>>
>> 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 --------------------------------------------------------------
>> 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 --------------------------------------------
>> 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
>> 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 -----------------------------
>> 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
>>
>>
>>
>
>
> --
> Christian Pagé
> http://meteocentre.com/ http://meteoalerte.com/
> Etudiant au Doctorat en Sciences de l'environnement UQAM
>
> _______________________________________________
> ncl-talk mailing list
> ncl-talk@ucar.edu
> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
>
_______________________________________________
ncl-talk mailing list
ncl-talk@ucar.edu
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
This archive was generated by hypermail 2b29 : Tue Oct 12 2004 - 17:06:34 MDT