;======================================================== load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl" load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl" ;======================================================== begin ; YYYY = 2010 ; MM = 12 ; CC = 1200 ; p and q used for the 9 point smoother in smth9_Wrap call p = 0.50 q = 0.25 ; heavy smoothing ;q = -0.25 ; light smoothing if (.not. isvar("YYYY")) then print (" ") print (" ") print ("Enter the year on the command line.") print ("ncl -n YYYY=2015 MM=12 CC=1200 read_narr_ncar-help.ncl") print (" ") print (" ") exit end if ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;This uses just two files for input while testing, the 30 and 31st of the month. ;Use with months having 31 days. ; myFiles = "Archive/"+YYYY+""+MM+"/narr-a_221_"+YYYY+""+MM+"3?_"+CC+"_000.sub.nc4" myFiles = "Archive/"+YYYY+""+MM+"/narr-a_221_"+YYYY+""+MM+"??_"+CC+"_000.sub.nc4" ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; print ("Input is from ") print (myFiles) fileList = systemfunc("/bin/ls "+myFiles) fin = addfiles (fileList,"r") print("File list has been added") ; print(fin) ; printVarSummary (fin) print("Number of files to open: "+ ListCount(fin)) ; ListSetType (fin, "cat") ListSetType (fin, "join") varnames = getfilevarnames(fin[0]) print (varnames) if(.not.any(ismissing(varnames))) then do i=0,dimsizes(varnames)-1 ; printFileVarSummary (fin[0],varnames(i)) ; print ("") end do end if ; Get lat/lon/time info first lat = fin[0]->gridlat_221 print ("Done reading lat") lon = fin[0]->gridlon_221 print ("Done reading lon") hgt3 = smth9_Wrap(fin[:]->HGT_221_ISBL,p,q,False) tmp3 = smth9_Wrap(fin[:]->TMP_221_ISBL,p,q,False) ; printVarSummary(hgt3) ; if(any(ismissing(hgt3(:,0,:,:)))) then ; numMissing = num(ismissing(hgt3(:,0,:,:))) ; print("There are "+numMissing+" missing data") ; ; katie1D = ndtooned(hgt3(:,0,:,:)) ; print(ind(ismissing(katie1D))) ; ; ; end if ; ; printVarSummary (katie) ;exit hgt_4k =hgt3(:,2,:,:) hgt_8k =hgt3(:,1,:,:) hgt_12k=hgt3(:,0,:,:) hgt_4k_avg =dim_avg_n(hgt_4k,0) hgt_8k_avg =dim_avg_n(hgt_8k,0) hgt_12k_avg=dim_avg_n(hgt_12k,0) ; create metadata for temperature tmp_4k =tmp3(:,2,:,:) tmp_8k =tmp3(:,1,:,:) tmp_12k=tmp3(:,0,:,:) ; printVarSummary (hgt_4k) ; printVarSummary (hgt_8k) ; printVarSummary (hgt_12k) pmsl = smth9_Wrap(fin[:]->PRMSL_221_MSL,p,q,False) pmsl_avg = dim_avg_n(pmsl, 0) pres = smth9_Wrap(smth9_Wrap(fin[:]->PRES_221_SFC,p,q,False),p,q,False) pres_avg = dim_avg_n(pres, 0) ; ; Read suface height from static data staticIn = addfile("AWIP32-fixed.grb","r") stnhgt2d = smth9_Wrap(smth9_Wrap(staticIn->HGT_221_SFC,p,q,False),p,q,False) stnhgt = conform_dims(dimsizes(pres), stnhgt2d, (/1,2/)) printVarSummary (stnhgt) printVarSummary (pres) ;============================================================================== ; Calculations ;============================================================================== ; Constants m2f = 3.28084 ; meters to feet conversion Pa2mb = 0.01 ; pascals to millibars Rd = 8314.32/28.9644 ; gas constant for dry air from US Std Atmos values g = 9.80665 ; gravity To=288.15 ; initial temperature L = -0.0065 ; [negative of] the lapse rate (Celsius/Kelvin per meter) c = -Rd*L/g ; exponent constant po = 101325.0 ; initial pressure in Pa p_12k = 65000. ; pressure at 12,000 feet p_8k = 75000. ; pressure at 8,000 feet p_4k = 87500. ; pressure at 4,000 feet ; Calculate altimeter setting from surface pressure altstg = smth9_Wrap(smth9_Wrap(pres *((1-((po/pres)^c) * (L*stnhgt/To))^(1/c)),p,q,False),p,q,False) printVarSummary (altstg) ; Determining Pressure Altitude PA_12k = (To/L)*((p_12k/po)^c - 1) PA_8k = (To/L)*((p_8k/po)^c - 1) PA_4k = (To/L)*((p_4k/po)^c - 1) ; Calculating corresponding temp. at altitude in Kelvin Tstd_12k = To + L*PA_12k Tstd_8k = To + L*PA_8k Tstd_4k = To + L*PA_4k ; calculate difference between actual temperature and standard temperature dT_12k = tmp_12k - Tstd_12k dT_8k = tmp_8k - Tstd_8k dT_4k = tmp_4k - Tstd_4k ; calculating indicated altitude - old way ; hia_4k = (To/L)*((p_4k/po)^c - (pmsl/po)^c) ; hia_8k = (To/L)*((p_8k/po)^c - (pmsl/po)^c) ; hia_12k =(To/L)*((p_12k/po)^c - (pmsl/po)^c) ; Use altstg (altimeter setting not PMSL hia_4k = (To/L)*((p_4k/po)^c - (altstg/po)^c) hia_8k = (To/L)*((p_8k/po)^c - (altstg/po)^c) hia_12k =(To/L)*((p_12k/po)^c - (altstg/po)^c) ; calculating Rule of Thumb at different pressure levels ROT_12k = (hia_12k - stnhgt)* 0.004*dT_12k * m2f ROT_8k = (hia_8k - stnhgt)* 0.004*dT_8k * m2f ROT_4k = (hia_4k - stnhgt)* 0.004*dT_4k * m2f ; Mask for terrain ROT_4k = mask(ROT_4k, (hia_4k - stnhgt).lt. 0.0,False) ; calculating corrected D-value at different pressure levels Dc_4k = m2f*(hgt_4k - hia_4k) ; Mask areas where the terrain (stnhgt) is above the indicated altitude Dc_4k = mask(Dc_4k,(hia_4k - stnhgt).lt.0.0,False) ; Still need to mask for these Dc_8k = m2f*(hgt_8k - hia_8k) Dc_12k = m2f*(hgt_12k - hia_12k) printVarSummary(Dc_4k) if(any(ismissing(Dc_4k))) then numMissing = num(ismissing(Dc_4k)) print("There are "+numMissing+" missing data") end if ;exit ; averaging the corrected D-value Dc_4k_avg = dim_avg_n(Dc_4k,0) Dc_8k_avg = dim_avg_n(Dc_8k,0) Dc_12k_avg = dim_avg_n(Dc_12k,0) ; averageing the Rule of Thumb ROT_4k_avg = dim_avg_n(ROT_4k,0) ROT_8k_avg = dim_avg_n(ROT_8k,0) ROT_12k_avg = dim_avg_n(ROT_12k,0) ;printVarSummary(ROT_12k) ;printVarSummary(ROT_12k_avg) ; calculating difference between Dc and ROT delta_4k = Dc_4k - ROT_4k delta_8k = Dc_8k - ROT_8k delta_12k = Dc_12k - ROT_12k ; averaging delta delta_4k_avg = dim_avg_n_Wrap(delta_4k,0) delta_8k_avg = dim_avg_n_Wrap(delta_8k,0) delta_12k_avg = dim_avg_n_Wrap(delta_12k,0) ; calculating rmsd rmsd_4k_avg = dim_rmsd_n_Wrap(Dc_4k,ROT_4k,0) rmsd_8k_avg = dim_rmsd_n_Wrap(Dc_8k,ROT_8k,0) rmsd_12k_avg = dim_rmsd_n_Wrap(Dc_12k,ROT_12k,0) ; printVarSummary(rmsd_12k_avg) ;=============================================================================== ; Metadata ;=============================================================================== ; Creating Meta Data for Dc_4k_avg Dc_4k_avg!0 = "gridx_221" Dc_4k_avg!1 = "gridy_221" Dc_4k_avg@initial_time = pmsl@initial_time Dc_4k_avg@forecast_time_units = "hours" Dc_4k_avg@forecast_time = 0 Dc_4k_avg@model = "North American Regional Reanalysis (NARR)" Dc_4k_avg@grid_number = 221 Dc_4k_avg@level_indicator = 100 Dc_4k_avg@coordinates = "gridlat_221, gridlon_221" Dc_4k_avg@_FillValue = 1e+20 Dc_4k_avg@units = "gpm" Dc_4k_avg@long_name = "Corrected D-Value 4k" Dc_4k_avg@center = "US National Weather Service - NCEP (WMC)" Dc_4k_avg@sub_center = "The North American Regional Reanalysis (NARR) Project" ; printVarSummary(Dc_4k_avg) ; Creating Meta Data for Dc_8k_avg Dc_8k_avg!0 = "gridx_221" Dc_8k_avg!1 = "gridy_221" Dc_8k_avg@initial_time = pmsl@initial_time Dc_8k_avg@forecast_time_units = "hours" Dc_8k_avg@forecast_time = 0 Dc_8k_avg@model = "North American Regional Reanalysis (NARR)" Dc_8k_avg@grid_number = 221 Dc_8k_avg@level_indicator = 100 Dc_8k_avg@coordinates = "gridlat_221 gridlon_221" Dc_8k_avg@_FillValue = 1e+20 Dc_8k_avg@units = "gpf" Dc_8k_avg@long_name = "Corrected D-Value 8k" Dc_8k_avg@center = "US National Weather Service - NCEP (WMC)" Dc_8k_avg@sub_center = "The North American Regional Reanalysis (NARR) Project" ; printVarSummary(Dc_8k_avg) ; Creating Meta Data for Dc_12k_avg Dc_12k_avg!0 = "gridx_221" Dc_12k_avg!1 = "gridy_221" Dc_12k_avg@initial_time = pmsl@initial_time Dc_12k_avg@forecast_time_units = "hours" Dc_12k_avg@forecast_time = 0 Dc_12k_avg@model = "North American Regional Reanalysis (NARR)" Dc_12k_avg@grid_number = 221 Dc_12k_avg@level_indicator = 100 Dc_12k_avg@coordinates ="gridlat_221 gridlon_221" Dc_12k_avg@_FillValue = 1e+20 Dc_12k_avg@units = "gpf" Dc_12k_avg@long_name = "Corrected D-Value 12k" Dc_12k_avg@center = "US National Weather Service - NCEP (WMC)" Dc_12k_avg@sub_center = "The North American Regional Reanalysis (NARR) Project" ; printVarSummary (Dc_12k_avg) ; Creating Metadata for ROT_4k ROT_4k_avg!0 = "gridx_221" ROT_4k_avg!1 = "gridy_221" ROT_4k_avg@initial_time = pmsl@initial_time ROT_4k_avg@forecast_time_units = "hours" ROT_4k_avg@forecast_time = 0 ROT_4k_avg@model = "North American Regional Reanalysis (NARR)" ROT_4k_avg@grid_number = 221 ROT_4k_avg@level_indicator = 100 ROT_4k_avg@coordinates = "gridlat_221, gridlon_221" ROT_4k_avg@_FillValue = 1e+20 ROT_4k_avg@units = "gpf" ROT_4k_avg@long_name = "Rule of Thumb 4k Height Correction" ROT_4k_avg@center = "US National Weather Service - NCEP (WMC)" ROT_4k_avg@sub_center = "The North American Regional Reanalysis (NARR) Project" ; printVarSummary(ROT_4k_avg) ; Creating Metadata for ROT_8k ROT_8k_avg!0 = "gridx_221" ROT_8k_avg!1 = "gridy_221" ROT_8k_avg@initial_time = pmsl@initial_time ROT_8k_avg@forecast_time_units = "hours" ROT_8k_avg@forecast_time = 0 ROT_8k_avg@model = "North American Regional Reanalysis (NARR)" ROT_8k_avg@grid_number = 221 ROT_8k_avg@level_indicator = 100 ROT_8k_avg@coordinates = "gridlat_221 gridlon_221" ROT_8k_avg@_FillValue = 1e+20 ROT_8k_avg@units = "gpf" ROT_8k_avg@long_name = "Rule of Thumb 8k Height Correction" ROT_8k_avg@center = "US National Weather Service - NCEP (WMC)" ROT_8k_avg@sub_center = "The North American Regional Reanalysis (NARR) Project" ; printVarSummary(ROT_8k_avg) ; Creating Metadata for ROT_12k ROT_12k_avg!0 = "gridx_221" ROT_12k_avg!1 = "gridy_221" ROT_12k_avg@initial_time = pmsl@initial_time ROT_12k_avg@forecast_time_units = "hours" ROT_12k_avg@forecast_time = 0 ROT_12k_avg@model = "North American Regional Reanalysis (NARR)" ROT_12k_avg@grid_number = 221 ROT_12k_avg@level_indicator = 100 ROT_12k_avg@coordinates ="gridlat_221 gridlon_221" ROT_12k_avg@_FillValue = 1e+20 ROT_12k_avg@units = "gpf" ROT_12k_avg@long_name = "Rule of Thumb 12k Height Correction" ROT_12k_avg@center = "US National Weather Service - NCEP (WMC)" ROT_12k_avg@sub_center = "The North American Regional Reanalysis (NARR) Project" ; printVarSummary (ROT_12k_avg) ; Creating Meta Data for delta_4k_avg delta_4k_avg!0 = "gridx_221" delta_4k_avg!1 = "gridy_221" delta_4k_avg@initial_time = pmsl@initial_time delta_4k_avg@forecast_time_units = "hours" delta_4k_avg@forecast_time = 0 delta_4k_avg@model = "North American Regional Reanalysis (NARR)" delta_4k_avg@grid_number = 221 delta_4k_avg@level_indicator = 100 delta_4k_avg@coordinates = "gridlat_221, gridlon_221" delta_4k_avg@_FillValue = 1e+20 delta_4k_avg@units = "gpf" delta_4k_avg@long_name = "D-value minus ROT 4k" delta_4k_avg@center = "US National Weather Service - NCEP (WMC)" delta_4k_avg@sub_center = "The North American Regional Reanalysis (NARR) Project" ; printVarSummary(delta_4k_avg) ; Creating Meta Data for delta_8k_avg delta_8k_avg!0 = "gridx_221" delta_8k_avg!1 = "gridy_221" delta_8k_avg@initial_time = pmsl@initial_time delta_8k_avg@forecast_time_units = "hours" delta_8k_avg@forecast_time = 0 delta_8k_avg@model = "North American Regional Reanalysis (NARR)" delta_8k_avg@grid_number = 221 delta_8k_avg@level_indicator = 100 delta_8k_avg@coordinates = "gridlat_221 gridlon_221" delta_8k_avg@_FillValue = 1e+20 delta_8k_avg@units = "gpf" delta_8k_avg@long_name = "D-value minus ROT 8k" delta_8k_avg@center = "US National Weather Service - NCEP (WMC)" delta_8k_avg@sub_center = "The North American Regional Reanalysis (NARR) Project" ; printVarSummary(delta_8k_avg) ; Creating Meta Data for delta_12k_avg delta_12k_avg!0 = "gridx_221" delta_12k_avg!1 = "gridy_221" delta_12k_avg@initial_time = pmsl@initial_time delta_12k_avg@forecast_time_units = "hours" delta_12k_avg@forecast_time = 0 delta_12k_avg@model = "North American Regional Reanalysis (NARR)" delta_12k_avg@grid_number = 221 delta_12k_avg@level_indicator = 100 delta_12k_avg@coordinates ="gridlat_221 gridlon_221" delta_12k_avg@_FillValue = 1e+20 delta_12k_avg@units = "gpf" delta_12k_avg@long_name = "D-value minus ROT 12k" delta_12k_avg@center = "US National Weather Service - NCEP (WMC)" delta_12k_avg@sub_center = "The North American Regional Reanalysis (NARR) Project" ; printVarSummary (delta_12k_avg) ; Creating Meta Data for rmsd_4k_avg rmsd_4k_avg!0 = "gridx_221" rmsd_4k_avg!1 = "gridy_221" rmsd_4k_avg@initial_time = pmsl@initial_time rmsd_4k_avg@forecast_time_units = "hours" rmsd_4k_avg@forecast_time = 0 rmsd_4k_avg@model = "North American Regional Reanalysis (NARR)" rmsd_4k_avg@grid_number = 221 rmsd_4k_avg@level_indicator = 100 rmsd_4k_avg@coordinates = "gridlat_221, gridlon_221" rmsd_4k_avg@_FillValue = 1e+20 rmsd_4k_avg@units = "gpf" rmsd_4k_avg@long_name = "RMSD of D-Value and ROT 4k" rmsd_4k_avg@center = "US National Weather Service - NCEP (WMC)" rmsd_4k_avg@sub_center = "The North American Regional Reanalysis (NARR) Project" ; printVarSummary(rmsd_4k_avg) ; Creating Meta Data for rmsd_8k_avg rmsd_8k_avg!0 = "gridx_221" rmsd_8k_avg!1 = "gridy_221" rmsd_8k_avg@initial_time = pmsl@initial_time rmsd_8k_avg@forecast_time_units = "hours" rmsd_8k_avg@forecast_time = 0 rmsd_8k_avg@model = "North American Regional Reanalysis (NARR)" rmsd_8k_avg@grid_number = 221 rmsd_8k_avg@level_indicator = 100 rmsd_8k_avg@coordinates = "gridlat_221 gridlon_221" rmsd_8k_avg@_FillValue = 1e+20 rmsd_8k_avg@units = "gpf" rmsd_8k_avg@long_name = "RMSD of D-Value and ROT 8k" rmsd_8k_avg@center = "US National Weather Service - NCEP (WMC)" rmsd_8k_avg@sub_center = "The North American Regional Reanalysis (NARR) Project" ; printVarSummary(rmsd_8k_avg) ; Creating Meta Data for rmsd_12k_avg rmsd_12k_avg!0 = "gridx_221" rmsd_12k_avg!1 = "gridy_221" rmsd_12k_avg@initial_time = pmsl@initial_time rmsd_12k_avg@forecast_time_units = "hours" rmsd_12k_avg@forecast_time = 0 rmsd_12k_avg@model = "North American Regional Reanalysis (NARR)" rmsd_12k_avg@grid_number = 221 rmsd_12k_avg@level_indicator = 100 rmsd_12k_avg@coordinates ="gridlat_221 gridlon_221" rmsd_12k_avg@_FillValue = 1e+20 rmsd_12k_avg@units = "gpf" rmsd_12k_avg@long_name = "RMSD of D-Value and ROT 12k" rmsd_12k_avg@center = "US National Weather Service - NCEP (WMC)" rmsd_12k_avg@sub_center = "The North American Regional Reanalysis (NARR) Project" ; printVarSummary (rmsd_12k_avg) ; Creating Meta Data for altstg ; altstg!0 = "time" altstg!1 = "gridx_221" altstg!2 = "gridy_221" ; altstg!0 = "gridx_221" ; altstg!1 = "gridy_221" altstg@initial_time = pmsl@initial_time altstg@forecast_time_units = "hours" altstg@forecast_time = 0 altstg@model = "North American Regional Reanalysis (NARR)" altstg@grid_number = 221 altstg@level_indicator = 100 altstg@coordinates ="gridlat_221 gridlon_221" altstg@_FillValue = 1e+20 altstg@units = "Pa" altstg@long_name = "Altimeter setting - calculated" altstg@center = "US National Weather Service - NCEP (WMC)" altstg@sub_center = "The North American Regional Reanalysis (NARR) Project" printVarSummary (altstg) ; Creating Meta Data for pres ; pres!0 = "time" pres!1 = "gridx_221" pres!2 = "gridy_221" ; pres!0 = "gridx_221" ; pres!1 = "gridy_221" pres@initial_time = pmsl@initial_time pres@forecast_time_units = "hours" pres@forecast_time = 0 pres@model = "North American Regional Reanalysis (NARR)" pres@grid_number = 221 pres@level_indicator = 100 pres@coordinates ="gridlat_221 gridlon_221" pres@_FillValue = 1e+20 pres@units = "Pa" pres@long_name = "Station Pres" pres@center = "US National Weather Service - NCEP (WMC)" pres@sub_center = "The North American Regional Reanalysis (NARR) Project" printVarSummary (pres) ; Creating Meta Data for stnhgt2d ; stnhgt2d!0 = "time" ; stnhgt2d!1 = "gridx_221" ; stnhgt2d!2 = "gridy_221" stnhgt2d!0 = "gridx_221" stnhgt2d!1 = "gridy_221" stnhgt2d@initial_time = pmsl@initial_time stnhgt2d@forecast_time_units = "hours" stnhgt2d@forecast_time = 0 stnhgt2d@model = "North American Regional Reanalysis (NARR)" stnhgt2d@grid_number = 221 stnhgt2d@level_indicator = 100 stnhgt2d@coordinates ="gridlat_221 gridlon_221" stnhgt2d@_FillValue = 1e+20 stnhgt2d@units = "gpm" stnhgt2d@long_name = "Surface Elevation" stnhgt2d@center = "US National Weather Service - NCEP (WMC)" stnhgt2d@sub_center = "The North American Regional Reanalysis (NARR) Project" printVarSummary (stnhgt2d) ;;;;;;;;;;;;;;;;;;;;;;; ; this block needs to be debugged statsFile = "Stats_"+MM+"/D-Value_Monthly_Stats."+YYYY+"."+MM+".averages."+CC+".nc" system("/bin/rm -f "+statsFile) ; remove any pre-existing file myNCDF = addfile(statsFile,"c") ; Create output file ;=================================================================== ; create global attributes of the file ;=================================================================== fAtt = True ; assign file attributes fAtt@title = "Monthly statistics for period beginning at "+pmsl@initial_time fAtt@valid_date = pmsl@initial_time fAtt@source_file = statsFile fAtt@creation_date = systemfunc ("date") fileattdef( myNCDF, fAtt ) ; copy file attributes ;=================================================================== ; make time an UNLIMITED dimension; recommended for most applications ;=================================================================== filedimdef(myNCDF,"time",-1,True) ;=================================================================== ; output variables directly; NCL will call appropriate functions ; to write the meta data associated with each variable ;=================================================================== myNCDF->gridlat_221 = lat myNCDF->gridlon_221 = lon ; myNCDF->time = time0 myNCDF->Dc_4k_avg = Dc_4k_avg myNCDF->Dc_8k_avg = Dc_8k_avg myNCDF->Dc_12k_avg = Dc_12k_avg myNCDF->ROT_4k_avg = ROT_4k_avg myNCDF->ROT_8k_avg = ROT_8k_avg myNCDF->ROT_12k_avg = ROT_12k_avg myNCDF->delta_4k_avg = delta_4k_avg myNCDF->delta_8k_avg = delta_8k_avg myNCDF->delta_12k_avg = delta_12k_avg myNCDF->rmsd_4k_avg = rmsd_4k_avg myNCDF->rmsd_8k_avg = rmsd_8k_avg myNCDF->rmsd_12k_avg = rmsd_12k_avg ; myNCDF->altstg = altstg(0,:,:) ; myNCDF->pres = pres(0,:,:) ; myNCDF->stnhgt2d = stnhgt2d ;exit ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;;; ;;; create plot ;;;; ; wks_type = "png" ; wks_type@wkWidth = 1200 ; wks_type@wkHeight = 1200 ; wks_type@wkOrientation = "landscape" ; outfile = "sfc_hgt" ; wks = gsn_open_wks(wks_type,outfile) ; open wk station ; gsn_define_colormap(wks,"BkBlAqGrYeOrReViWh200") ; load color table lat2d = fin[0]->gridlat_221 lon2d = fin[0]->gridlon_221 dimlc = dimsizes(lat2d) nlat = dimlc(0) mlon = dimlc(1) ; United States Boundary (rough estimate) latS = 20 latN = 60 lonW = -120 lonE = -75 ji = region_ind(lat2d, lon2d, latS, latN, lonW, lonE) jStrt = ji(0) ; lat start jLast = ji(1) ; lat last iStrt = ji(2) ; lon start iLast = ji(3) ; lon last ; Create subset of full data for plotting Dc_4k_small_avg = Dc_4k_avg(jStrt:jLast, iStrt:iLast) Dc_8k_small_avg = Dc_8k_avg(jStrt:jLast, iStrt:iLast) Dc_12k_small_avg = Dc_12k_avg(jStrt:jLast, iStrt:iLast) ROT_4k_small_avg = ROT_4k_avg(jStrt:jLast, iStrt:iLast) ROT_8k_small_avg = ROT_8k_avg(jStrt:jLast, iStrt:iLast) ROT_12k_small_avg = ROT_12k_avg(jStrt:jLast, iStrt:iLast) ;---------------------------------------------------- ; create plots ;---------------------------------------------------- wks = gsn_open_wks ("png", "narr-paper") ; send graphics to PNG file cmap = read_colormap_file("BlueWhiteOrangeRed") ; read color data res = True ; plot mods desired for original grid res@cnFillOn = True ; color fill res@cnFillPalette = cmap(2:,:) ; set color map res@cnLinesOn = True ; add contour lines ;res@cnMinLevelValF = -400. ; sets minimum value for contour lines ; res@cnMaxLevelValF = 200. ; sets maximum value for contour lines ;res@cnLevelSpacingF = 50. ; sets contour interval res@cnLevelSelectionMode = "ExplicitLevels" res@cnLevels = (/ -400, -350, -300, -250, -200, -150, -100, -50, 0, 50, 100, 150, 200, 250, 300, 350, 400/) res@mpGridAndLimbOn = False res@mpDataBaseVersion = "MediumRes" ; include National/US Boundaries res@mpOutlineBoundarySets = "GeophysicalAndUSStates" ; draws area outlines of US res@pmTickMarkDisplayMode = "Always" ; turn on tickmarks res@tmXTOn = False res@gsnAddCyclic = False ; regional data res@gsnMaximize = True res@mpLimitMode = "Corners" ; choose range of map res@mpLeftCornerLatF = lat2d(jStrt,iStrt) res@mpLeftCornerLonF = lon2d(jStrt,iStrt) res@mpRightCornerLatF = lat2d(jLast,iLast) res@mpRightCornerLonF = lon2d(jLast,iLast) res@tfDoNDCOverlay = True ;************* Is this the mapping issue? ; res@NhlTOverlayMode = "NDCDataExtent" res@mpProjection = "LambertConformal" res@mpLambertMeridianF = (lonW-lonE)/2+lonE res@gsnLeftString = "Original grid" plot = new(6,graphic) ; upper left res@tiMainString = "Corrected D-Value 4,000ft" plot(0) = gsn_csm_contour_map(wks,Dc_4k_small_avg,res) ; upper right res@tiMainString = "Corrected D-Value 8,000ft" plot(2) = gsn_csm_contour_map(wks,Dc_8k_small_avg,res) ; middle left res@tiMainString = "Corrected D-Value 12,000ft" plot(4) = gsn_csm_contour_map(wks,Dc_12k_small_avg,res) ; middle right res@tiMainString = "Rule of Thumb 4,000ft" plot(1) = gsn_csm_contour_map(wks,ROT_4k_small_avg,res) ; lower left res@tiMainString = "Rule of Thumb 8,000ft" plot(3) = gsn_csm_contour_map(wks,ROT_8k_small_avg,res) ; lower right res@tiMainString = "Rule of Thumb 12,000ft" plot(5) = gsn_csm_contour_map(wks,ROT_12k_small_avg,res) ; draw the panel plot gsn_panel(wks,plot,(/3,2/),False) end exit ;========================================================================== ; Plot Troubleshooting ;========================================================================== ;res@mpLimitMode = "Corners" ; choose range of map ;res@mpLeftCornerLatF = lat2d(60,110) ;res@mpLeftCornerLonF = lon2d(60,110) ;res@mpRightCornerLatF = lat2d(nlat-120,mlon-75) ; res@mpRightCornerLonF = lon2d(nlat-120,mlon-75) ; res@mpLeftCornerLatF = lat2d(0,0) ; res@mpLeftCornerLonF = lon2d(0,0) ; res@mpRightCornerLatF = lat2d(nlat-1,mlon-1) ; res@mpRightCornerLonF = lon2d(nlat-1,mlon-1) ; res@tfDoNDCOverlay = True ; res@mpProjection = "LambertConformal" ; res@mpLambertParallel1F = lat2d@mpLambertParallel1F ; res@mpLambertParallel2F = lat2d@mpLambertParallel2F ; res@mpLambertMeridianF = lat2d@mpLambertMeridianF ;res@mpLambertMeridianF = (lonW-lonE)/2+lonE ; res@gsnCenterString = "U@"+x&lv_ISBL0(10) + "hPa" ; draw center subtitle ;res@gsnLeftString = "Original grid" ; draw left subtitle ; plot = gsn_csm_contour_map(wks,Dc_4k_small_avg,res) ; Draw original grid on map ; plot = gsn_csm_contour_map(wks,Dc_8k_small_avg,res) ; Draw original grid on map ; plot = gsn_csm_contour_map(wks,Dc_12k_small_avg,res) ; Draw original grid on map ; plot = gsn_csm_contour_map(wks,hgt_8k_small*m2f,res) ; Draw original grid on map ; plot = gsn_csm_contour_map(wks,ROT_4k_small_avg,res) ; Draw original grid on map ; plot = gsn_csm_contour_map(wks,ROT_8k_small_avg,res) ; Draw original grid on map ;plot = gsn_csm_contour_map(wks,ROT_12k_small_avg,res) ; Draw original grid on map ; end ; plot = gsn_csm_contour_map(wks,hgt,res1) ; create plot ; end