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
-- *Mike Madden*
_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Wed Mar 26 16:46:04 2014
This archive was generated by hypermail 2.1.8 : Mon Mar 31 2014 - 11:47:09 MDT