Re: LCL and CAPE in skewt_func.ncl

From: James Boyle (boyle5 AT llnl.gov)
Date: Tue Oct 12 2004 - 11:20:01 MDT

  • Next message: Takeshi Enomoto: "Re: combining data"

    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