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

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