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

From: Dennis Shea <shea_at_nyahnyahspammersnyahnyah>
Date: Thu Mar 27 2014 - 07:37:19 MDT

[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
>
_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Thu Mar 27 07:37:28 2014

This archive was generated by hypermail 2.1.8 : Mon Mar 31 2014 - 11:47:09 MDT