Re: Vertical cross sections in Mirage

Forgot to mention the error I get when I run it within Mirage.

warning:ContourPlotInitialize: scalar field is constant; ContourPlot not possible:[errno=1102]

Thank you!

On Dec 17, 2012, at 2:16 PM, Cecille Villanueva-Birriel wrote:

> Hello,
>
> I am trying to make vertical cross sections using a WRF output file located in analysis cluster Mirage. Everything works fine if I run this script on my computer. But for some reason it doesn't work in Mirage. I would for the code to work in Mirage since I will need to analyze big WRF data files. I won't be able to do it in my computer due to the low memory allocations. Is there something I can do? Any help will be greatly appreciated. Here is the ncl code.
>
> Thanks,
> Cecille
>
>
> ; Example script to produce plots for a WRF real-data run,
> ; with the ARW coordinate dynamics option.
> ; Plot data on a cross section
> ; This script will plot data from a a given point A to point B
> ; Vertical coordinate is height
>
> ;load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
> load "./gsn_code.ncl"
> ;load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"
> load "./WRFUserARW.ncl"
>
> begin
>
> scheme = "MOR"
> city = "Jasper"
> dir = "ICCP_project/"
> period = "fut"
> grid = "1km"
>
> path = "/glade/scratch/cvillanu/Morrison/"
> ;
> ; The WRF ARW input file.
> ; This needs to have a ".nc" appended, so just do it.
> a = addfile(path+"wrfout_d01_0001-01-01_00:00:00_MOR_1km.nc","r")
>
> ; We generate plots, but what kind do we prefer?
> ; type = "x11"
> type = "pdf"
> ; type = "ps"
> ; type = "ncgm"
> wks = gsn_open_wks(type,"W_vcross_1km")
>
> ; Set some basic resources
> res = True
> res@MainTitle = "REAL-TIME WRF"
>
> pltres = True
>
> ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
>
> FirstTime = True
> times = wrf_user_list_times(a) ; get times in the file
> ntimes = dimsizes(times) ; number of times in the file
>
> mdims = getfilevardimsizes(a,"P") ; get some dimension sizes for the file
> nd = dimsizes(mdims)
>
> ;---------------------------------------------------------------
>
> ; do it = 0,ntimes-1,2 ; TIME LOOP
>
> do it = 0,ntimes-1
> ; do it = 5, 13
>
> print("Working on time: " + times(it) )
> res@TimeLabel = times(it) ; Set Valid time to use on plots
>
> tc = wrf_user_getvar(a,"tc",it) ; T in C
> rh = wrf_user_getvar(a,"rh",it) ; relative humidity
> z = wrf_user_getvar(a, "z",it) ; grid point height
> w = wrf_user_getvar(a,"wa",it) ; vertical velocity
> tk = wrf_user_getvar(a,"tk",it)
>
> p = wrf_user_getvar(a,"PB",it)
> rgas = 287.04 ; J Kg-1 K-1
>
> qv = wrf_user_getvar(a, "QVAPOR",it) ;
>
> if(isfilevar(a,"QCLOUD"))
> qc = wrf_user_getvar(a, "QCLOUD",it)
> end if
> if(isfilevar(a,"QRAIN"))
> qr = wrf_user_getvar(a, "QRAIN",it)
> end if
> if(isfilevar(a,"QSNOW"))
> qs = wrf_user_getvar(a, "QSNOW",it)
> end if
> if(isfilevar(a,"QICE"))
> qi = wrf_user_getvar(a, "QICE",it)
> end if
> if(isfilevar(a,"QGRAUP"))
> qg = wrf_user_getvar(a, "QGRAUP",it)
> end if
>
> rho = p/(rgas*tk*(1+0.61*qv))
>
> PWC = qr*rho
> LWC = (qr+qc)*rho
>
> if ( FirstTime ) then ; get height info for labels
> zmin = 0.
> zmax = max(z)/1000.
> nz = floattoint(zmax/2 + 1)
> FirstTime = False
> end if
>
> ;---------------------------------------------------------------
>
> do ip = 1, 24 ; we are doing 3 plots, specifying different start and end points
>
> opts = True ; setting start and end times
> plane = new(4,float)
>
> if(ip .eq. 1) then
> plane = (/ 0,16, 100,16 /) ; start x;y & end x;y point
> end if
> if(ip .eq. 2) then
> plane = (/ 0,17, 100,17 /) ; start x;y & end x;y point
> end if
> if(ip .eq. 3) then
> plane = (/ 0,18, 100,18 /) ; start x;y & end x;y point
> end if
> if(ip .eq. 4) then
> plane = (/ 0,19, 100,19 /) ; start x;y & end x;y point
> end if
> if(ip .eq. 5) then
> plane = (/ 0,20, 100,20 /) ; start x;y & end x;y point
> end if
> if(ip .eq. 6) then
> plane = (/ 0,21, 100,21 /) ; start x;y & end x;y point
> end if
> if(ip .eq. 7) then
> plane = (/ 0,22, 100,22 /) ; start x;y & end x;y point
> end if
> if(ip .eq. 8) then
> plane = (/ 0,23, 100,23 /) ; start x;y & end x;y point
> end if
> if(ip .eq. 9) then
> plane = (/ 0,24, 100,24 /) ; start x;y & end x;y point
> end if
> if(ip .eq. 10) then
> plane = (/ 0,25, 100,25 /) ; start x;y & end x;y point
> end if
> if(ip .eq. 11) then
> plane = (/ 0,26, 100,26 /) ; start x;y & end x;y point
> end if
> if(ip .eq. 12) then
> plane = (/ 0,27, 100,27 /) ; start x;y & end x;y point
> end if
> if(ip .eq. 13) then
> plane = (/ 0,28, 100,28 /) ; start x;y & end x;y point
> end if
> if(ip .eq. 14) then
> plane = (/ 0,29, 100,29 /) ; start x;y & end x;y point
> end if
> if(ip .eq. 15) then
> plane = (/ 0,30, 100,30 /) ; start x;y & end x;y point
> end if
> if(ip .eq. 16) then
> plane = (/ 0,31, 100,31 /) ; start x;y & end x;y point
> end if
> if(ip .eq. 17) then
> plane = (/ 0,32, 100,32 /) ; start x;y & end x;y point
> end if
> if(ip .eq. 18) then
> plane = (/ 0,33, 100,33 /) ; start x;y & end x;y point
> end if
> if(ip .eq. 19) then
> plane = (/ 0,34, 100,34 /) ; start x;y & end x;y point
> end if
> if(ip .eq. 20) then
> plane = (/ 0,35, 100,35 /) ; start x;y & end x;y point
> end if
> if(ip .eq. 21) then
> plane = (/ 0,36, 100,36 /) ; start x;y & end x;y point
> end if
> if(ip .eq. 22) then
> plane = (/ 0,37, 100,37 /) ; start x;y & end x;y point
> end if
> if(ip .eq. 23) then
> plane = (/ 0,38, 100,38 /) ; start x;y & end x;y point
> end if
> if(ip .eq. 24) then
> plane = (/ 0,39, 100,39 /) ; start x;y & end x;y point
> end if
>
>
> gsn_define_colormap(wks, "BlWhRe") ; Define a new color map
> ; Vertical Velocity
> w_plane = wrf_user_intrp3d(w,z,"v",plane,0.,opts)
> ; w_plane = 100.*w_plane
> ;res@TimeLabel = it
> ;res@vpWidthF = plot_width
> ;res@vpHeightF = plot_height
> opts_w = res
> ;opts_w@FieldTitle = w@description
> opts_w@UnitLabel = "m/s"
> ; opts_w@PlotLevelID = 0.001*height + " km"
> opts_w@cnFillOn = True
> opts_w@cnLevelSelectionMode = "ManualLevels" ; set manual contour levels
> opts_w@cnMinLevelValF = -38 ; set min contour level
> opts_w@cnMaxLevelValF = 38. ; set max contour level
> opts_w@cnLevelSpacingF = 2 ; set contour spacing
>
>
>
> ; Basic Options for all cloud plots
> opts_clouds = res
> opts_clouds@cnFillOn = True
> opts_clouds@UnitLabel = "g/kg"
>
> if (isvar("qv"))
> qv_plane = wrf_user_intrp3d(qv,z,"v",plane,0.,opts)
> qv_plane = qv_plane*1000.
> opts_qv = opts_clouds
> end if
>
> if (isvar("qc"))
> qc_plane = wrf_user_intrp3d(qc,z,"v",plane,0.,opts)
> qc_plane = qc_plane*1000.
> opts_qc = opts_clouds
> end if
>
> if (isvar("qr"))
> qr_plane = wrf_user_intrp3d( qr,z,"v",plane,0.,opts)
> qr_plane = qr_plane*1000.
> opts_qr = opts_clouds
> end if
>
> if (isvar("PWC"))
> ; opts_clouds@UnitLabel = "Precipitation Mass (g/m3)"
> PWC_plane = wrf_user_intrp3d( PWC,z,"v",plane,0.,opts)
> PWC_plane = PWC_plane*1000.
> opts_PWC = opts_clouds
> end if
>
> if (isvar("LWC"))
> ; opts_clouds@UnitLabel = "Liquid Water Content (g/m3)"
> LWC_plane = wrf_user_intrp3d( LWC,z,"v",plane,0.,opts)
> LWC_plane = LWC_plane*1000.
> opts_LWC = opts_clouds
> end if
>
>
> if (isvar("qs"))
> qs_plane = wrf_user_intrp3d( qs,z,"v",plane,0.,opts)
> qs_plane = qs_plane*1000.
> opts_qs = opts_clouds
> end if
>
> if (isvar("qi"))
> qi_plane = wrf_user_intrp3d( qi,z,"v",plane,0.,opts)
> qi_plane = qi_plane*1000.
> opts_qi = opts_clouds
> end if
>
> if (isvar("qg"))
> qg_plane = wrf_user_intrp3d( qg,z,"v",plane,0.,opts)
> qg_plane = qg_plane*1000.
> opts_qg = opts_clouds
> end if
>
>
>
> tc_plane = wrf_user_intrp3d(tc,z,"v",plane,0.,opts)
>
> dim = dimsizes(qv_plane) ; Find the data span - for use in labels
> zspan = dim(0)
>
>
> ; Options for XY Plots
> opts_xy = res
> opts_xy@tiYAxisString = "Height (km)"
> opts_xy@AspectRatio = 0.75
> opts_xy@cnMissingValPerimOn = True
> opts_xy@cnMissingValFillColor = 0
> opts_xy@cnMissingValFillPattern = 11
> opts_xy@tmYLMode = "Explicit"
> opts_xy@tmYLValues = fspan(0,zspan,nz) ; Create tick marks
> opts_xy@tmYLLabels = sprintf("%.1f",fspan(zmin,zmax,nz)) ; Create labels
> opts_xy@tiXAxisFontHeightF = 0.020
> opts_xy@tiYAxisFontHeightF = 0.020
> opts_xy@tmXBMajorLengthF = 0.02
> opts_xy@tmYLMajorLengthF = 0.02
> opts_xy@tmYLLabelFontHeightF = 0.015
> opts_xy@PlotOrientation = tc_plane@Orientation
>
>
> ; Plotting options for RH
> opts_clouds = opts_xy
> opts_clouds@pmLabelBarOrthogonalPosF = -0.07
> ;opts_clouds@ContourParameters = (/ 10., 90., 10. /)
> ;opts_clouds@ContourParameters = (/ 1., 20., 1. /)
> opts_w = opts_xy
> opts_PWC = opts_xy
> opts_LWC = opts_xy
>
> opts_clouds@cnFillOn = True
> opts_clouds@cnFillColors = (/"White","White","White", \
> "White","Chartreuse","Green", \
> "Green3","Green4", \
> "ForestGreen","PaleGreen4"/)
>
> ; Plotting options for Temperature
> opts_tc = opts_xy
> opts_tc@cnInfoLabelOrthogonalPosF = 0.00
> opts_tc@ContourParameters = (/ 5. /)
>
>
> ; Get the contour info for the rh and temp
> contour_tc = wrf_contour(a,wks,tc_plane,opts_tc)
> contour_w = wrf_contour(a,wks, w_plane,opts_w)
> contour_qv = wrf_contour(a,wks,qv_plane,opts_clouds)
> contour_qr = wrf_contour(a,wks, qr_plane,opts_clouds)
> contour_qc = wrf_contour(a,wks, qc_plane,opts_clouds)
> contour_qs = wrf_contour(a,wks, qs_plane,opts_clouds)
> contour_qi = wrf_contour(a,wks, qi_plane,opts_clouds)
> contour_qg = wrf_contour(a,wks, qg_plane,opts_clouds)
> contour_PWC = wrf_contour(a,wks, PWC_plane,opts_PWC)
> contour_LWC = wrf_contour(a,wks, LWC_plane,opts_LWC)
>
> ; MAKE PLOTS
> plot = wrf_overlays(a,wks,(/contour_w,contour_tc/),pltres)
>
> ; Delete options and fields, so we don't have carry over
> delete(opts_tc)
> delete(opts_clouds)
> delete(tc_plane)
> delete(qv_plane)
> delete(qc_plane)
> delete(qr_plane)
> delete(qs_plane)
> delete(qi_plane)
> delete(qg_plane)
> delete(w_plane)
> delete(PWC_plane)
> delete(LWC_plane)
>
> end do ; make next cross section
>
> ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
>
> end do ; END OF TIME LOOP
>
> end
>> ********************************************************
>> Cecille M. Villanueva Birriel
>> Ph.D. Candidate
>> Cloud Microphysics Research Group
>> Purdue University
>> Department of Earth, Atmospheric, & Planetary Sciences
>> 550 Stadium Mall Drive, West Lafayette, IN 47907-2051
>> email: cvillanu@purdue.edu
>> ********************************************************
>

_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Mon Dec 17 12:22:29 2012