Re: LCL and CAPE in skewt_func.ncl

From: Christian Pagé (page.christian AT XXXXXX)
Date: Tue Oct 12 2004 - 10:30:38 MDT


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 AT 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 AT 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 AT ucar.edu http://mailman.ucar.edu/mailman/listinfo/ncl-talk



This archive was generated by hypermail 2b29 : Tue Oct 12 2004 - 17:06:26 MDT