Hi all,
Just wanted to follow up that this question was very helpfully answered
by a number of people, and I really appreciate all the information that
was provided.
Just so a problem solution gets archived, I'll attach one of the scripts
that was generously sent my way in case anyone else has this problem
down the road.
Thanks,
Kelly Mahoney
ncl-talk-request@ucar.edu wrote:
> Send ncl-talk mailing list submissions to
> ncl-talk@ucar.edu
>
> To subscribe or unsubscribe via the World Wide Web, visit
> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
> or, via email, send a message with subject or body 'help' to
> ncl-talk-request@ucar.edu
>
> You can reach the person managing the list at
> ncl-talk-owner@ucar.edu
>
> When replying, please edit your Subject line so it is more specific
> than "Re: Contents of ncl-talk digest..."
>
>
> Today's Topics:
>
> 1. Re: About Wrf_cape_3d.shtml function (Dennis Shea)
> 2. Re: Hourly precipitation plots for hourly WRF output (Mary Haley)
> 3. Re: do loop + if statement problems with CALIPSO data (Mary Haley)
> 4. Re: fatal:XyPlotChanges: error setting overlay object view (Kim)
>
>
> ----------------------------------------------------------------------
>
> Message: 1
> Date: Tue, 23 Feb 2010 19:46:18 -0700
> From: Dennis Shea <shea@ucar.edu>
> Subject: Re: About Wrf_cape_3d.shtml function
> To: z3303149@unsw.edu.au
> Cc: ncl-talk@ucar.edu
> Message-ID: <4B8492FA.3060603@ucar.edu>
> Content-Type: text/plain; charset=ISO-8859-1; format=flowed
>
> The documentation says the fortran codes used are from RIP.
> Regardless of the data source, if you input in the correct
> variables *with the correct units*, the function should
> calculate CAPE.
>
> The CAPE code used by the skewT script is specific
> to the skewT script. It can not be used independently.
>
> --
> Kerry Emanuel [MIT] makes his CAPE fortran code available at
>
>
> http://wind.mit.edu/~emanuel/ftpacc.html
>
> Perhaps, you can call that from NCL
>
> http://www.ncl.ucar.edu/Document/Manuals/Getting_Started/beyondbasics.shtml#IncorporateYourOwn
>
> Good luck
>
>
>
> z3303149@unsw.edu.au wrote:
>
>> Hello Dennis
>>
>> i used these functions before but i got wrong result and when i asked the wrf
>> help about the problem, they told me this functions was developed for wrf model
>> and there is a large change that this will not work for other model output. i
>> work in MM5 model.
>>
>> Thanks
>>
>> Ali
>>
>>
>> Hi group
>>
>> i'm trying to calculate the CAPE from skew T digram , all my variables 4D
>> (time, level,lat,lon) i attachment the figure, it work just in (eps) and i
>> can't open it as (x11) or (pdf),but i got the figure without the red line of the
>> parcel lapse rate like in example no. 2 in skew examples because i want to
>> calculate the CAPE, if anyone can help me please.
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>
>> load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
>> load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
>> load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/skewt_func.ncl"
>>
>>
>>
>> begin
>>
>>
>>
>>
>>
>>
>>
>> colmin = 40
>> colmax = 47
>> rowmin = 33
>> rowmax = 40
>>
>>
>>
>> class = 6
>>
>>
>>
>> ;1999
>> diri = "/srv/science/maths1/z3236814/data02/swap/co2/fc/"
>>
>>
>>
>> f = addfile(diri+"evb_minclass_200noev_2095.nc","r") ; open data netcdf file
>>
>> topo = f->terrain
>> w = f->w(class-1,:,:,:,:)
>> u = f->u(class-1,:,:,:,:)
>> v = f->v(class-1,:,:,:,:)
>> t = f->t(class-1,:,:,:,:)
>> q = f->q(class-1,:,:,:,:)
>> pp = f->pp(class-1,:,:,:,:)
>> sigma = f->sigma_level
>> ps = f->pstarcrs
>> time = f->time
>> lat = f->latitcrs
>> lon = f->longicrs
>> ptop = f->ptop
>>
>>
>>
>>
>> dsize = dimsizes(q)
>> pres = new(dsize,float)
>>
>>
>> ;to calculate the pressure for each grid point
>> pres = pp + ptop +conform(pp,ps,(/2,3/))*conform(pp,sigma,1)
>>
>>
>> rel = relhum(t,q,pres)
>>
>>
>>
>> ;to calculates the dewpoint temperature
>> dewpoint= dewtemp_trh(t,rel)
>>
>>
>> dewpoint = dewpoint -273 ;returned depoint temperature in (C)
>>
>>
>> ; to define the dewpoint
>> dewpoint!0 = "time"
>> dewpoint!1 = "sigma_level"
>> dewpoint!2 = "i_cross"
>> dewpoint!3 = "j_cross"
>>
>> pres = pres*0.01 ; to converte it to (hPa)
>>
>> ;calculate virtual temperature
>> vir_temp = new(dsize,float)
>> vir_temp = t*(1.+q*0.61)
>>
>>
>>
>> ; to define the pres
>> pres!0 = "time"
>> pres!1 = "sigma_level"
>> pres!2 = "i_cross"
>> pres!3 = "j_cross"
>>
>> ; to define the virtual temperature
>> vir_temp!0 = "time"
>> vir_temp!1 = "sigma_level"
>> vir_temp!2 = "i_cross"
>> vir_temp!3 = "j_cross"
>>
>>
>>
>>
>>
>> t = t-273 ; to converte it to (C)
>>
>>
>>
>>
>> gh = hydro(pres(time|:,i_cross|:,j_cross|:,sigma_level|::-1), \
>> vir_temp(time|:,i_cross|:,j_cross|:,sigma_level|::-1), \
>> conform(vir_temp(:,0,:,:),topo,(/1,2/)))
>>
>>
>>
>>
>> ; to define the geopotential height
>> gh!0 = "time"
>> gh!1 = "i_cross"
>> gh!2 = "j_cross"
>> gh!3 = "sigma_level"
>>
>>
>>
>>
>>
>>
>>
>> ;to calculates the meteorological wind direction
>> r2d = 45.0/atan(1.0) ;to conversion from radians to degrees because the return
>> value in radians
>> wdir = atan2(u,v)* r2d + 180
>>
>>
>>
>> wdir!0 = "time"
>> wdir!1 = "sigma_level"
>> wdir!2 = "i_cross"
>> wdir!3 = "j_cross"
>>
>>
>> wks = gsn_open_wks ("eps", "skewclass62095")
>>
>>
>> pres1=
>> dim_avg(dim_avg(pres(time|:,sigma_level|:,i_cross|rowmin:rowmax,j_cross|colmin:colmax)))
>> z =
>> dim_avg(dim_avg(gh(time|:,sigma_level|:,i_cross|rowmin:rowmax,j_cross|colmin:colmax)))
>> wind =
>> dim_avg(dim_avg(v(time|:,sigma_level|:,i_cross|rowmin:rowmax,j_cross|colmin:colmax)))
>> windd=
>> dim_avg(dim_avg(wdir(time|:,sigma_level|:,i_cross|rowmin:rowmax,j_cross|colmin:colmax)))
>> temp =
>> dim_avg(dim_avg(t(time|:,sigma_level|:,i_cross|rowmin:rowmax,j_cross|colmin:colmax)))
>> dew1 =
>> dim_avg(dim_avg(dewpoint(time|:,sigma_level|:,i_cross|rowmin:rowmax,j_cross|colmin:colmax)))
>>
>>
>>
>>
>>
>>
>> skewtOpts = True
>> skewtOpts@DrawFahrenheit = False
>> skewtOpts@DrawWind = True
>> skewtOpts@DrawHeightScale = True
>>
>>
>> dataOpts = True ; options describing data and ploting
>> dataOpts@linePatternDewPt = 1
>> dataOpts@colDewPt = "red"
>>
>> ;add legend
>> legend = create "Legend" legendClass wks
>> "vpXF" : 0.15
>> "vpYF" : 0.88
>> "vpWidthF" : 0.24 ; width
>> "vpHeightF" : 0.09 ; height
>> "lgPerimOn" : True ; perimeter
>> "lgItemCount" : 2 ; how many
>> "lgLabelStrings" : (/"Temperature","Dew point Temperature"/)
>> "lgLineThicknessF" : 3.0
>> "lgBoxMinorExtentF" : 0.2
>> end create
>>
>>
>>
>>
>> skewt_bkgd = skewT_BackGround (wks, skewtOpts)
>>
>> skewt_data=skewT_PlotData(wks,skewt_bkgd,pres1(15,:),temp(15,:),dew1(15,:),z(15,:),wind(15,:),windd(15,:),dataOpts)
>>
>>
>>
>>
>>
>> draw (skewt_data)
>> draw (skewt_bkgd)
>>
>> draw(legend)
>> frame(wks)
>>
>>
>>
>>
>> end
>>
>> ----- End forwarded message -----
>>
>>
>>
>>
>> ------------------------------------------------------------------------
>>
>> _______________________________________________
>> ncl-talk mailing list
>> List instructions, subscriber options, unsubscribe:
>> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
>>
>
>
> ------------------------------
>
> Message: 2
> Date: Tue, 23 Feb 2010 22:05:28 -0700
> From: Mary Haley <haley@ucar.edu>
> Subject: Re: Hourly precipitation plots for hourly WRF
> output
> To: ncl-talk forum <ncl-talk@ucar.edu>
> Message-ID: <19AF8976-2071-4B4D-B747-B1A61DF14F94@ucar.edu>
> Content-Type: text/plain; charset=US-ASCII; format=flowed; delsp=yes
>
> [This message was forwarded to wrfhelp@ucar.edu.]
>
> On Feb 23, 2010, at 2:51 PM, Kelly Mahoney wrote:
>
>
>> Hi,
>>
>> I'm a new NCL user and am having some trouble figuring out the best
>> way
>> to generate plots of hourly precipitation from my hourly wrfout_*
>> files.
>>
>> There is a sample script online (
>> http://www.mmm.ucar.edu/wrf/users/graphics/NCL/Examples/Precip/wrf_Precip.ncl)
>> that almost gets the job done, but I'm having trouble because I have a
>> list of individual files with hourly output at one time, instead of
>> one
>> file with several times in it. [As an aside, I also tried using NCO to
>> concatenate all of the files into one huge file, and running the
>> example
>> script on it, but I then got errors referencing issues with the
>> WRFUserARW.ncl file; an example is pasted below.]
>> I have been trying to use addfiles to generate my list but then have
>> lots of subsequent issues with functions like wrf_user_list_times,
>> wrf_map, wrf_user_getvar, etc, that require a single file instead of a
>> variable of type "list" as is generated with addfiles.
>>
>> Ideally, I would love to have a script that loops over each time/
>> file in
>> my model run and generates a separate plot of hourly precipitation for
>> that time (I can do this for total precipitation easily, as that only
>> needs the time in the file itself, whereas hourly precip needs to be
>> calculated by subtracting off the total precip from the time before).
>>
>> If anyone has a similar script that they wouldn't mind sharing, or any
>> advice on how this would be most easily done, I would be very
>> appreciative!
>>
>> Thanks so much,
>> Kelly Mahoney
>>
>>
>> ---------------------------------------------------------
>> errors when trying to run sample script with one large, NCO-
>> concatenated
>> file:
>>
>>> Coordinates:
>>> fatal:Number of subscripts do not match number of dimensions of
>>> variable,(2) Subscripts used, (3) Subscripts expected
>>> fatal:Execute: Error occurred at or near line 1373 in file
>>> $NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl
>>>
>>>
>> _______________________________________________
>> ncl-talk mailing list
>> List instructions, subscriber options, unsubscribe:
>> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
>>
>
>
>
> ------------------------------
>
> Message: 3
> Date: Tue, 23 Feb 2010 22:21:41 -0700
> From: Mary Haley <haley@ucar.edu>
> Subject: Re: do loop + if statement problems with CALIPSO
> data
> To: Mike <mtosca@uci.edu>
> Cc: ncl-talk@ucar.edu
> Message-ID: <ADA5E76C-38B9-475C-B174-C660BE2C47C6@ucar.edu>
> Content-Type: text/plain; charset="us-ascii"
>
> Mike,
>
> I think you need to put in some debug writes to determine where the
> problem is.
>
> For one, do you know if the code is ever making it inside the "if
> (.not.ismissing(latN_lgth(k)) .and. latN_lgth(k) .eq. 134) then"
> statement? If not, then layer_top will always be set to -9999.
>
> First, though, I think we can clean up this code so there aren't as
> many nested "do" loops. Do loops run slow in NCL (as well as other
> interpreted languages).
>
> For example, instead of:
>
>
>> do l = 0,3
>> do m=0,134
>> if (ftr_chk(k,m,l) .eq. 6) then
>> layer_top_(k,m,l) = layer_top(k,m,l)
>> else
>> layer_top_(k,m,l) = -9999
>> end if
>> end do
>> end do
>>
>>
>
> you can use the "where" statement:
>
> layer_top_ = where(ftr_chk(k,:,:).eq.6,layer_top(k,:,:),-9999)
>
> To make sure that ftr_chk(k,:,:) actually contains some values that
> equal 6, you can print this:
>
> print("# of values where ftr_chk is 6 = " + num(ftr_chk(k,:,:).eq.
> 6))
>
> to get the count of how many values are equal to 6.
>
> Please put some debug statements like this in your code to see if you
> can narrow down the problem.
>
> --Mary
>
> On Feb 23, 2010, at 3:02 PM, Mike wrote:
>
>
>> I am using 5km CALIPSO aerosol layer data and the aerosol
>> classification and layer top variables. What I want to do is assign
>> a value in the array called "layer_top_" from the array "layer_top"
>> which I read from the file. However, I only want this value to be
>> assigned to "layer_top_" if the aerosol is classified as "smoke",
>> that is when ftr_chk (the feature classification, unpacked integer
>> value) = 6. The code runs fine, but when I go back and check to see
>> if everywhere the aerosol was classified as smoke, the corresponding
>> "layer top" value has been put into the array "layer_top_" I notice
>> that it didn't work correctly. Namely, ftr_chk(k,m,l) = 6, but
>> layer_top_(k,m,l) = -9999 even though layer_top(k,m,l) = *value* (as
>> reported in the original dataset). I've attached the code below (the
>> latN_top and latN_btm are just there so I can just remove values for
>> the latitude range I'm interested in):
>>
>> do k = 0,37
>> if (.not.ismissing(latN_lgth(k)) .and. latN_lgth(k) .eq. 134) then
>> ftr_flg(k,:,:) = ftr_flg_tst(k,latN_top(k):latN_btm(k),:)
>> ftr_chk(k,:,:) = dim_gbits(ftr_flg(k,:,:),4,3,13,8)
>> layer_top(k,:,:) = layer_top_tst(k, latN_top(k):latN_btm(k),:)
>> lon1(k,:) = lonN(k, latN_top(k):latN_btm(k))
>> do l = 0,3
>> do m=0,134
>> if (ftr_chk(k,m,l) .eq. 6) then
>> layer_top_(k,m,l) = layer_top(k,m,l)
>> else
>> layer_top_(k,m,l) = -9999
>> end if
>> end do
>> end do
>> else
>> ftr_flg(k,:,:) = -999
>> layer_top(k,:,:) = -9999
>> end if
>> end do
>>
>> Thanks for your help!
>>
>> -Mike
>>
>> Michael Tosca
>> Ph.D. Candidate
>> Dept. Earth System Science
>> University of California, Irvine
>> Irvine, CA 92697
>> http://dust.ess.uci.edu/mtosca/web
>>
>> _______________________________________________
>> ncl-talk mailing list
>> List instructions, subscriber options, unsubscribe:
>> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
>>
>
> -------------- next part --------------
> An HTML attachment was scrubbed...
> URL: http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20100223/1edacf67/attachment.html
>
> ------------------------------
>
> Message: 4
> Date: Wed, 24 Feb 2010 04:06:53 -0800 (PST)
> From: Kim <r4rid@yahoo.com>
> Subject: Re: fatal:XyPlotChanges: error setting overlay
> object view
> To: Mary Haley <haley@ucar.edu>
> Cc: ncl-talk@ucar.edu
> Message-ID: <719474.56190.qm@web37106.mail.mud.yahoo.com>
> Content-Type: text/plain; charset="utf-8"
>
> Dear Mary, thanks you so much spenting more time on my problem, but unfortunately it does not working in the light of your suggested methods.
> thanks again
> kim
>
> --- On Tue, 2/23/10, Mary Haley <haley@ucar.edu> wrote:
>
>
> From: Mary Haley <haley@ucar.edu>
> Subject: Re: fatal:XyPlotChanges: error setting overlay object view
> To: "Kim" <r4rid@yahoo.com>
> Cc: ncl-talk@ucar.edu
> Date: Tuesday, February 23, 2010, 7:23 AM
>
>
>
>
> I'm not sure if this will work, but instead of setting trXMinF and trXMaxF explicitly,
> try to index your time array for whichever values you want to actually plot.?
>
>
>
> For example, if you wanted to start at time "19900101" and end at "20001231", then
> you would do something like this:
>
>
>
> ??yrs = floattoint(utc(:,0)) ??
> ??ibeg = where(yrs.eq.1990.and.mons.eq.1.and.days.eq.1)
> ??iend = where(yrs.eq.2000.and.mons.eq.12.and.days.eq.31)
>
>
> Then:?
>
>
> ??plotA = gsn_csm_xy(wks, time(ibeg:iend), zRegion(ibeg:iend), res)
> ??. . .
> ??plotB = gsn_csm_xy(wks, time(ibeg_iend), z1Region(ibeg:iend), res)
>
>
>
>
>
> The other possibility is to use ut_calendar with an option of -2 to convert your time to YYYYMMDD, and then plot against these new values instead of your original time values.
>
>
> For example:
>
>
> ??new_time = ut_calendar(time, -2)
> ??. . .
> ??tm_major_values = new_time(day_1)
>
>
> Change these to the appropriate YYYYMMDD format:
>
>
>
>
>
>
>
>
>
> ???res@trXMinF = 17552088????????; set minimum X-axis value
> ???res@trXMaxF = 17555040???????; set maximum X-axis value
>
> and then:
>
>
>
> ??plotA = gsn_csm_xy(wks, new_time, zRegion, res)
> ??. . .
> ??plotB = gsn_csm_xy(wks, new_time, z1Region, res)
>
>
> Without knowing what your original time values look like, I'm not sure what's the best method.
>
>
>
>
> --Mary
>
>
>
> On Feb 22, 2010, at 8:02 PM, Kim wrote:
>
>
>
>
>
>
> Dear Mary,
> Please have a look on?my code is given below. If I don't set trXMinF and trXMaxF it is working but the figure is not correct the lines start?away from the y-axis. Please your help in this regard will be appreciable.
> thanks again,
> kim
> *******************
> --------------------------------------
> --------------------------------------
> ?? f = addfile("shtfl.sfc.gauss.2003.nc","r")
> ?? f1? = addfile("shtfl.sfc.gauss.2002.nc","r")
> ?? z = short2flt(f->shtfl(120:243,:,:))????
> ?? z1= short2flt(f1->shtfl(120:243,:,:))
> ?? printVarSummary(z)
> ?? printVarSummary(z1)
> ?? time = z&time
> ?? time = z1&time
> ?? utc? = ut_calendar(time,0)??? ; yr, mo, dy, hr, mn, sc
> ?? mons = floattoint(utc(:, 1))? ; Values from 1 to 12
> ?? days = floattoint(utc(:, 2))? ; Values from 1 to 31
> ?? day_1? = ind (days.eq.1)
> ?? mon_1? = mons(day_1)?????
> ???smonths? = (/"","J", "F", "M", "A", "MAY", "JUN",\
> ?? "JUL", "AUG", "SEP", "OCT", "N", "D"/)
> ?? tm_major_values = time(day_1)
> ?? tm_major_labels = smonths(mon_1)
> ?
> ?? zRegion = wgt_areaave_Wrap(z(:,{15:40},{50:100}),1.0,1.0,0)
> ?? z1Region = wgt_areaave_Wrap(z1(:,{15:40},{50:100}),1.0,1.0,0)
> ? ; printVarSummary(z1Region)
> ? ; printMinMax(z1Region, True)
> ? ;**********************************
> ?? wks?? = gsn_open_wks ("x11","sh")
> ?? res?????????????????? = True
> ?? res@gsnDraw? = False????????????????????????????
> ?? res@gsnFrame = False
> ?? res@tmLabelAutoStride = True
> ?? res@xyLineColors????? =(/"red"/)
> ?? res@xyLineThicknesses = 2.0
> ?? res@tiMainString????? = "?
> ?? res@tmXBMode????????? = "Explicit"?? ; Define own tick mark labels.
> ?? res@tmXBValues =?? tm_major_values
> ?? res@tmXBLabels =?? tm_major_labels
> ?? res@tmXTOn? = False?? ; Turn off tickmarks and labels
> ?? res@tmYROn? = False
> ? ?res@trXMinF = 17552088?????? ?; set minimum X-axis value
> ?? res@trXMaxF = 17555040? ?????; set maximum X-axis value
> ?? res@tiYAxisString??????? = "SHFL(W/m^2)"
> ? ?res@vpHeightF= 0.4
> ?? res@vpWidthF = 0.8
> ?? plotA? = gsn_csm_xy(wks, time,zRegion,res)
> ?;**************************************
> ?? res@xyLineColors????? = (/"blue"/)
> ??res@trXMinF = 17543328???????? ?; set minimum X-axis value
> ?? res@trXMaxF = 17546280????? ; set maximum X-axis value
> ?? plotB? = gsn_csm_xy(wks, time,z1Region,res)
> ?? overlay(plotA, plotB)
> ? ?draw(plotA)
> ?? frame(wks)
> ?? end
>
>
> --- On Mon, 2/22/10, Mary Haley <haley@ucar.edu> wrote:
>
>
> From: Mary Haley <haley@ucar.edu>
> Subject: Re: fatal:XyPlotChanges: error setting overlay object view
> To: "Kim" <r4rid@yahoo.com>
> Cc: ncl-talk@ucar.edu
> Date: Monday, February 22, 2010, 5:55 AM
>
>
> Dear Kim,
>
>
> How are you doing the overlay, and when are you setting trXMinF and trXMaxF? It would help if I could see a script.
>
>
> --Mary
>
>
>
> On Feb 21, 2010, at 9:57 PM, Kim wrote:
>
>
>
>
>
>
> Dear NCL users
> I?m trying to xy-plotA & plotB. Individually and also when I overlay (plotA, plotB) without trXMinF and ?trXMax resources it is working properly. When I set trXMinF and trXMax during their overlaying it does not work. If I not set these recourses then the lines are not starting at the initial point. Is there any resource or function available in NCL.
> Thanks
> kim
> ?**************
> fatal:["MultiText.c":940]:[errno=12]
> fatal:TickMarkSetValues: Could not set MultiText item values for bottom tick mark labels, cannot continue
> fatal:TickMarkSetValues: A fatal error was detected while computing tick marks labels,cannot continue
> fatal:PlotManagerSetValues: Error setting TickMark values
> fatal:XyPlotChanges: error setting overlay object view
> fatal:["MultiText.c":940]:[errno=12]
> fatal:TickMarkSetValues: Could not set MultiText item values for bottom tick mark labels, cannot continue
> fatal:TickMarkSetValues: A fatal error was detected while computing tick marks labels,cannot continue
> fatal:XyPlotChanges: error setting overlay object view
> fatal:["MultiText.c":940]:[errno=12]
> fatal:TickMarkSetValues: Could not set MultiText item values for bottom tick mark labels, cannot continue
> fatal:TickMarkSetValues: A fatal error was detected while computing tick marks labels,cannot continue
> fatal:XyPlotChanges: error setting overlay object view
> fatal:["MultiText.c":940]:[errno=12]
> fatal:TickMarkSetValues: Could not set MultiText item values for bottom tick mark labels, cannot continue
> fatal:TickMarkSetValues: A fatal error was detected while computing tick marks labels,cannot continue
> fatal:PlotManagerSetValues: Error setting TickMark values
> fatal:XyPlotChanges: error setting overlay object view
> _______________________________________________
> ncl-talk mailing list
> List instructions, subscriber options, unsubscribe:
> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
>
>
>
>
>
>
>
> -------------- next part --------------
> An HTML attachment was scrubbed...
> URL: http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20100224/62f6698c/attachment.html
>
> ------------------------------
>
> _______________________________________________
> ncl-talk mailing list
> ncl-talk@ucar.edu
> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
>
>
> End of ncl-talk Digest, Vol 75, Issue 62
> ****************************************
>
; Example script to produce plots for a WRF real-data run,
; with the ARW coordinate dynamics option.
; In this example we first get the entire field over time, which will
; make it easier to calculate tendencies
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"
begin
;
; The WRF ARW input file.
; This needs to have a ".nc" appended, so just do it.
diri = getenv("DIRI")
domain = getenv("DOMAIN")
Files = systemfunc("ls "+diri+"/wrfout_d0"+domain+"*00:00")
print(Files)
; We generate plots, but what kind do we prefer?
; type = "x11"
; type = "pdf"
type = "ps"
; type = "ncgm"
; Set some basic resources
res = True
res@MainTitle = "REAL-TIME WRF"
pltres = True
mpres = True
mpres@mpGeophysicalLineColor = "Black"
mpres@mpNationalLineColor = "Black"
mpres@mpUSStateLineColor = "Black"
mpres@mpGridLineColor = "Black"
mpres@mpLimbLineColor = "Black"
mpres@mpPerimLineColor = "Black"
mpres@mpFillOn = True
mpres@mpFillDrawOrder = "PreDraw"
mpres@mpLandFillColor = "DarkOliveGreen3"
mpres@mpOceanFillColor = -1
mpres@mpInlandWaterFillColor = -1
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; What times and how many time steps are in the data set?
ntimes = dimsizes(Files) ; number of times in the file
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
times = new(ntimes,string)
a = addfile(Files(0)+".nc","r")
times(0) = wrf_user_list_times(a) ; get times in the file
mpres@mpDataBaseVersion = "Ncarg4_1" ; higher res data base
if (a@GRID_ID .ge. 3) then
mpres@mpDataBaseVersion = "HighRes"
end if
do it = 1,ntimes-1,1
a = addfile(Files(it)+".nc","r")
a_last = addfile(Files(it-1)+".nc","r")
ff = floattointeger(a->XTIME/60.)
PlotName = "pp_d0"+domain+"_"+sprinti("%0.2i",ff)
print("PlotName: "+PlotName)
wks = gsn_open_wks(type,PlotName)
colors = (/"white","black","white","yellow","orange","DarkOrange",\
"OrangeRed","Red1","Red2","Red4","DarkOrchid1","purple",\
"MediumPurple3","Blue2","Blue4","DarkOliveGreen3"/)
gsn_define_colormap(wks, colors)
times(it) = wrf_user_list_times(a) ; get times in the file
print("Working on time: " + times(0) )
res@TimeLabel = times(it) ; Set Valid time to use on plots
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; First get the variables we will need
slp = wrf_user_getvar(a,"slp",0) ; slp
wrf_smooth_2d( slp, 3 ) ; smooth slp
; Get non-convective, convective
; Calculate total precipitation
rain_exp = a->RAINNC(0,:,:)
rain_con = a->RAINC(0,:,:)
rain_exp = rain_exp - a_last->RAINNC(0,:,:)
rain_con = rain_con - a_last->RAINC(0,:,:)
rain_tot = rain_exp + rain_con
rain_tot@description = "Total Precipitation"
; Plotting options for Sea Level Pressure
opts_psl = res
opts_psl@ContourParameters = (/ 900., 1100., 2. /)
opts_psl@cnLineColor = "Blue"
opts_psl@cnInfoLabelOn = False
opts_psl@cnLineLabelFontHeightF = 0.01
opts_psl@cnLineLabelPerimOn = False
opts_psl@gsnContourLineThicknessesScale = 1.5
contour_psl = wrf_contour(a,wks,slp,opts_psl)
delete(opts_psl)
; Plotting options for Precipitation
opts_r = res
opts_r@UnitLabel = "mm"
opts_r@cnLevelSelectionMode = "ExplicitLevels"
opts_r@cnLevels = (/ .1, .5, 1., 2., 4., 8. /)
opts_r@cnFillColors = (/"transparent","Yellow","orange",\
"DarkOrange",\
"OrangeRed","Red1","Violet"/)
opts_r@cnInfoLabelOn = False
opts_r@cnConstFLabelOn = False
opts_r@cnFillOn = True
; Total Precipitation (color fill)
;; contour_tot = wrf_contour(a,wks, rain_tot, opts_r)
; Precipitation Tendencies
opts_r@SubFieldTitle = "from " + times(it-1) + " to " + times(it)
contour_tend = wrf_contour(a,wks, rain_tot,opts_r) ; total (color)
delete(opts_r)
; MAKE PLOTS
plot = wrf_map_overlays(a,wks,(/contour_tend,contour_psl/),pltres,mpres)
delete(a)
delete(a_last)
end do ; END OF TIME LOOP
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
end
_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Thu Feb 25 13:41:41 2010
This archive was generated by hypermail 2.1.8 : Mon Mar 01 2010 - 08:49:37 MST