Re: Vertical Integration of WRF PM10 with Known Thickness Lengths; Problems with Integration and Coordination of Dimensions

Hello,

These were all great tips! Thank you much!

Mike

On Thu, Mar 27, 2014 at 5:37 AM, Dennis Shea <shea@ucar.edu> wrote:

> [1] Eliminate the entire 'do it=' loop
> [2] You do not need: thick = new((/27/),float) ; define a new array
> NCL can determine the size of 'thick automatically'
> [3] Read documentation for conform_dims (conform, also), dim_sum_n
> and getfilevardimsizes
> [4] It is recommended to use the dim_*_n over the older dim_* functions.
> The latter are kept for backward compatibility only.
>
> ;--------
> a = addfile(...,"r")
> dimpm = getfilevardimsizes(a,"PM10") ; (0,1,2,3)
> thick = (/ ... /) ; (klev) => (27)
> THICK = conform_dims(dimpm, thick, 1) ; dim #1 is klev
> printVarSummary(THICK) ; (ntim,klev,nlat,mlon)
>
> pm = a->PM10 ; (ntim,klev,nlat,mlon)
> printVarSummary(pm)
>
> vertint = dim_sum_n(pm*THICK,1) ; (ntim,nlat,mlon)
> printVarSummary(vertint)
> delete(THICK) ; not necessary
>
> copy_VarCoords(pm(:,0,:,:), vertint) ; copy appropriate coord info
> vertint@long_name = "..." ; assign attributes
> vertint@units = "..."
> printVarSummary(vertint)
> ;--------
>
>
> On 3/26/14, 4:45 PM, James Madden wrote:
>
>> Hello,
>>
>> I am trying to perform a vertical integration of PM10 in layers across the
>> entire domain. I already have the thicknesses of each layer.
>>
>> However, I am unable to perform the integration, and coordinate the arrays
>> into something manageable for plotting.
>>
>> I have attached the code below. In the code, I have notated where I have
>> stumbled into problems, and I give explanations at those spots.
>>
>> Am I missing something much simpler? I feel like I am making this problem
>> way too complex.
>>
>> I have a dimension issue, but I have handled these problems before, and if
>> I introduce more loops to handle the latitude and longitude grid cells, my
>> plots will take way too long to form.
>>
>> The problem is later in the code, after all of the resources.
>>
>> Thank you much,
>>
>> Mike
>>
>> 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/wrf/WRF_contributed.ncl"
>> load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"
>> load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"
>>
>> begin
>>
>> wks = gsn_open_wks("x11","vert_int_plot")
>> gsn_define_colormap(wks,"BlAqGrYeOrReVi200")
>>
>> monthname = "June"
>> month = "06"
>> day = "24"
>> year = "2009"
>>
>> ; Read WRF File
>> a =
>> addfile("/import/c/w/jmmadden/WRFruns/wrfout_d01_2009-"+
>> month+"-"+day+"_00:00:00"+".nc","r")
>>
>> ; Set Resources
>> ;----------------------------------------------------------------------
>> ; Begin plotting
>> ;----------------------------------------------------------------------
>> res = True ; plot mods desired
>> res@gsnMaximize = True ; uncomment to maximize size
>> res@gsnSpreadColors = True ; use full range of colormap
>> res@cnFillOn = True ; color plot desired
>> res@cnLinesOn = False ; turn off contour lines
>> res@cnLineLabelsOn = False ; turn off contour labels
>> ; res@cnLevelSelectionMode = "ExplicitLevels"
>> ; Levels = (/3.,4.,5.,6.5,8.,10.,12.,15.,18.,21.,25.,28.,34.,39.,45/)
>> ; res@cnLevels = Levels
>> res@cnFillMode = "RasterFill" ; activate raster mode
>> res@cnRasterSmoothingOn = True
>> res@lbLabelAutoStride = True
>> res@NoHeaderFooter = True
>> ;----------------------------------------------------------------------
>> ; Use WRF_contributed procedure to set map resources
>> ;----------------------------------------------------------------------
>> mpres = True
>> pltres = True
>> pltres@NoTitles = True
>> mpres@mpGeophysicalLineColor = 1
>> mpres@mpGridLineColor = 1
>> mpres@mpCenterLatF = a@CEN_LAT
>> mpres@mpCenterLonF = a@CEN_LON
>> ;----------------------------------------------------------------------
>> ; set True for native projection (faster)
>> ;----------------------------------------------------------------------
>>
>>
>> ; HOURLY LOOP OF PM 10 PLOTS
>>
>> do it = 0,23
>> itp1=it+1
>> print("I am working on hour "+itp1)
>>
>> ; INSERT EXAMPLE THICKNESS LENGTHS
>> thick = new((/27/),float) ; define a new array
>> thick = (/ 7.53, 26.39, 45.39, 72.2, 107.417, 120.2, 133.4, 159.0,
>> 177.8,
>> 212.1, 213.2, 214.6, 263.764, 289.3, 326.5, 367.2, 402.1, 457.47, 809.45,
>> 916.6, 695.72, 750.3, 1049.4, 1282.5, 1582.8, 2874.39, 2875.0/)
>>
>> pm = a->PM10(it,:,:,:)
>> pm!0 = "lev"
>> pm!1 = "lat"
>> pm!2 = "lon"
>>
>> ;;;;;;;;;;;PROBLEM IS HERE;;;;;;;;;;;;;;;;
>> ; integrated = (thick) * (pm(:,0,0)) ;********
>> integrated = (thick) * (pm) ; *problem is here* ;*********
>>
>> vertint = dim_sum(integrated(lat|:,lon|:,lev|:)) ;*********
>>
>> ; I understand here that the dimensions do not match, and I have solved
>> problems like this before. However, the latitude and longitude must be
>> included, as I am vertically integrating the PM10 with each grid-cell.
>>
>> ; However, I cannot create loops of "lat = 0,159" and "lon = 0,199". The
>> plots will take way too long to generate. I believe I need to multiply
>> the
>> first dimension of "pm" by the thicknesses, but I am not sure how to do
>> so.
>>
>> ; Is there something simple that I am missing? I believe there is a much
>> simpler option. All I need to do is vertically integrate each layer, and
>> I
>> already have the thicknesses to do so.
>>
>>
>> vertint@description="Vertically Integrated PM~B~10~N~ Concentration
>> "+itp1+"00 UTC "+monthname+" "+day+", "+year
>> vertint@units="~F33~m~F21~g ~N~m~S~-2~N~"
>>
>> ref_contour = wrf_contour(a,wks,vertint(5:154,5:194),res)
>> plot1 = wrf_map_overlays(a,wks,(/ref_contour/),pltres,mpres)
>>
>> end do
>> end
>>
>>
>> Thank you much once again.
>>
>> Mike
>>
>>
>>
>> _______________________________________________
>> ncl-talk mailing list
>> List instructions, subscriber options, unsubscribe:
>> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
>>
>>

-- 
*Mike Madden*
*Graduate Research Assistant *
*Department of Atmospheric Sciences *
*University of Alaska Fairbanks*
*IARC 338N*
*Office:  (907) 474-7618*
*Cell: (417)
439-2830---------------------------------------------------------------------"Buy
the ticket, take the ride."  Hunter S. Thompson*

_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Thu Mar 27 11:56:47 2014