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/csm/contributed.ncl" load "$NCARG_ROOT/lib/ncarg/nclscripts/contrib/ut_string.ncl" begin yr = 1982 imageDirectory = "/keyserlab_rit/kbiernat/ms_paper_1/Jan_1982/DT_theta_tracks_400km_circle_Jan1982_reducedCB_paperRes_redBlue" ;*********************************************************************************************************************** ; Time Configurations ;*********************************************************************************************************************** syyyy = yr ; Start Date smm = 1 sdd = 1 shh = 00 eyyyy = yr ; End Date emm = 1 edd = 13 ehh = 18 timeUnits = "hours since 1800-01-01 00:00:00" ; Time is represented as "hours since 1800-01-01 00:00:00" startTime = ut_inv_calendar(syyyy,smm,sdd,shh,00,00,timeUnits,0) ; Convert start time to the above mentioned units endTime = ut_inv_calendar(eyyyy,emm,edd,ehh,00,00,timeUnits,0) ; Convert end time to the above mentioned units ntimes = toint(((endTime - startTime)/6) + 1) ; Find the number of times from the start time to the end time. ; E.g., if startTime = 700024 and endTime = 700048, and the data is ; available every 6 hours, we'd expect ntimes = 5. ; We see that (700048 - 700024)/6 + 1 = 4+1 = 5 timestepfactor = 4 ;in 6hr intervals times = new( (/ntimes/), "float") ; Create a new 1-dimensional float array of size nTimes times@units = timeUnits ; Give the new times array the same units as timeUnits (in this case: ; "hours since 1800-01-01 00:00:00") do n=0,ntimes-1 ; This do loop fills the times array with the times from the startTime to the endTime times(n) = tofloat(startTime + (6*n)) ; in 6 hour increments, in units of "hours since 1800-01-01 00:00:00" end do ;************************************************ ; create pointer to file and read in data ;************************************************ fileName = "/keyserlab_rit/kbiernat/erai_ecmwf/"+yr+"/u_pv."+yr+".nc" inFile_u_pv = addfile( fileName, "r" ) u_pv = inFile_u_pv->u_pv({times},:,:) fileName = "/keyserlab_rit/kbiernat/erai_ecmwf/"+yr+"/v_pv."+yr+".nc" inFile_v_pv = addfile( fileName, "r" ) v_pv = inFile_v_pv->v_pv({times},:,:) fileName = "/keyserlab_rit/kbiernat/erai_ecmwf/"+yr+"/pt_pv."+yr+".nc" inFile_theta_pv = addfile( fileName, "r" ) theta_pv = inFile_theta_pv->pt_pv({times},:,:) ws = (u_pv*u_pv + v_pv*v_pv)^(0.5) ; Compute Wind Speed (m/s) copy_VarCoords(u_pv,ws) ; Copies all named dimensions and coordinate variables from variable "u" to variable "wind_speed" copy_VarMeta(u_pv,ws) ;print(ws) ws = smth9(ws, 0.50, 0.25, False) ws = smth9(ws, 0.50, 0.25, False) ;*********************************************************************************************************** ; Specify Attributes for Map Settings ;*********************************************************************************************************** res = True ; Turn on resources for map res@gsnDraw = False ; Do not draw plot until specified otherwise res@gsnFrame = False ; Do not advance frame until specified otherwise ;res@gsnPolar = "NH" ; specify the hemisphere ;res@mpMinLatF = 30 ; minimum lat to plot ;res@pmTickMarkDisplayMode = "Always" res@mpProjection = "Stereographic" res@mpCenterLatF = 90 res@mpCenterLonF = -95 res@mpGridPolarLonSpacingF = 0 res@mpLimitMode = "Corners" res@mpLeftCornerLatF = 17 res@mpRightCornerLatF = 55 res@mpLeftCornerLonF = -118 res@mpRightCornerLonF =-32 res@mpGridAndLimbOn = True res@mpGridLatSpacingF = 20 res@mpGridLonSpacingF = 20 ;res@mpCenterLonF = -40 res@mpGridLineColor = (/0.25,0.25,0.25/) res@mpGridLineThicknessF = 4.0 res@mpGridLineDashPattern = 10. res@mpGridAndLimbDrawOrder = "Draw" res@mpFillOn = False ; Turn off map fill ;res@mpOutlineDrawOrder = "Draw" res@mpFillDrawOrder = "PreDraw" res@mpOutlineOn = True ; Turn the map outline on res@mpOutlineBoundarySets = "National" ; draw national boundaries res@mpOutlineSpecifiers = (/"Canada", "United States"/) + " : States" ; draw US States, Canadian provinces res@mpOutlineDrawOrder = "Draw" res@mpGeophysicalLineThicknessF = 8.0 ; Change Thickness of the Geophysical Boundary Lines res@mpNationalLineThicknessF = 8.0 ; Change Thickness of National Boundary Lines res@mpUSStateLineThicknessF = 8.0 ; Change Thickness of US State and Canadian Province Boundaries res@mpGeophysicalLineColor = (/0.22,0.22,0.22/) res@mpNationalLineColor = res@mpGeophysicalLineColor res@mpUSStateLineColor = res@mpGeophysicalLineColor res@mpDataSetName = "Earth..4" ; Basically, this dataset draws more accurate coastlines/boundaries ; as well as state/provinces outlines. res@mpDataBaseVersion = 1 ; medium resolution res@mpDataResolution = 2 ; medium data resolution ;res@gsnMaximize = True ; make as large as possible res@gsnPaperOrientation = "landscape" ; Change orientation of display to landscape ;*********************************************************************************************************** ; Attributes for shaded theta on dynamic tropopause ;*********************************************************************************************************** thtres = True thtres@gsnDraw = False ; Do not draw plot until specified otherwise thtres@gsnFrame = False ; Do not advance frame until specified otherwise thtres@gsnLeftString = "" thtres@gsnRightString = "" thtres@gsnAddCyclic = True thtres@cnLevelSelectionMode = "ManualLevels" thtres@cnMinLevelValF = 252. ; set the minimum contour level thtres@cnMaxLevelValF = 384. ; set the maximum contour level thtres@cnLevelSpacingF = 6. ; set the contour interval thtres@cnFillOn = True ; turn on contour shading thtres@cnLinesOn = False ; don't draw the contour lines thtres@cnInfoLabelOn = False ; don't draw the info label thtres@cnLineLabelsOn = False ; don't draw the contour line labels thtres@tiMainFontHeightF = 0.012 ; set the main title font height thtres@cnFillDrawOrder = "PreDraw" thtres@lbLabelBarOn= False thtres@lbAutoManage = False thtres@lbBoxLineColor = 1 thtres@lbLabelAutoStride = True ;thtres@lbTopMarginF = 1.3 ;thtres@lbBottomMarginF = -1.3 thtres@lbLabelOffsetF = 0.075 thtres@lbOrientation = "Horizontal" thtres@pmLabelBarSide = "Right" thtres@pmLabelBarOrthogonalPosF = 0.005 ; Moves label bar further down away from plot ;thtres@lbTitleString = "" thtres@lbLabelFont = "helvetica" ;thtres@lbTitleFont = "helvetica" thtres@lbLabelFontHeightF = 0.009 ;thtres@lbTitleFontHeightF = 0.0095 ;thtres@pmLabelBarWidthF = 0.70 thtres@pmLabelBarHeightF = 0.04 ;thtres@lbTitlePosition = "Right" ;thtres@lbTitleDirection = "Across" ;thtres@lbTitleExtentF = 0.005 ;thtres@gsnSpreadColors = True thtres@cnFillColors = (/(/0.90,0.90,0.90/),\ (/0.75,0.75,0.75/),\ (/0.60,0.60,0.60/),\ (/0.44,0.45,0.45/),\ ;start gray (/0.58,0.11,0.38/),\ ;end pink (/0.74,0.35,0.58/),\ (/0.92,0.59,0.78/),\ (/0.38,0.20,0.54/),\ ; end purple (/0.54,0.38,0.69/),\ (/0.72,0.59,0.84/),\ ; start purple (/0.13,0.41,0.68/),\ ; end blue (/0.34,0.60,0.76/),\ (/0.61,0.81,0.87/),\ ; end blue (/0.98,0.89,0.76/),\ (/0.95,0.75,0.61/),\ (/0.94,0.62,0.48/),\ (/0.92,0.45,0.36/),\ (/0.93,0.20,0.19/),\ (/0.80,0.04,0.03/),\ (/0.86,0.20,0.37/),\ (/0.88,0.35,0.52/),\ (/0.89,0.50,0.67/),\ (/0.89,0.65,0.82/),\ (/0.90,0.78,0.94/)/) ;*********************************************************************************************************** ; Specify Attributes for DT wind speed ;*********************************************************************************************************** wsres = True wsres@gsnDraw = False ; Do not draw plot until specified otherwise wsres@gsnFrame = False ; Do not advance frame until specified otherwise wsres@gsnLeftString = "" ; Make blank the string above the plot to the left (usually a string representing variable name) wsres@gsnRightString = "" ; Make blank the string above the plot to the right (usually a string representing variable units) wsres@cnLinesOn = True ; Dispalay contour lines wsres@cnLevelSelectionMode = "ExplicitLevels" ; use explicit levelsp wsres@cnLevels = ispan(50,100,10) ; set the contour levels wsres@cnLineColor = 1 ; wsres@cnFillOn = False ; Turn off contour fill. wsres@cnLineLabelsOn = True ; Turn on line labels. wsres@cnInfoLabelOn = False ; Turn off info label. wsres@cnLineThicknessF = 7.5 wsres@cnLabelMasking = True ; Contour lines are masked where contour labels appear wsres@cnLineLabelDensityF = 1.0 wsres@cnLineLabelAngleF = 0.0 ;wsres@cnLineLabelPlacementMode = 2 wsres@cnLineLabelBackgroundColor = -1 ; set background box of label to transparent wsres@cnMonoLineLabelFontColor = True ; Allow font colors to be displayed on all labels wsres@cnLineLabelFont = "Helvetica" wsres@cnLineLabelInterval = 1 wsres@cnLineLabelFontColor = wsres@cnLineColor ; set color of labels on contours (can only use color table) wsres@cnLineLabelFontThicknessF = 1.5 ; Line Label Font thickness (multiplier of unit thickness of 1.0) wsres@cnLineLabelFontHeightF = 0.012 ; Change Font Size of Contour Labels wsres@cnLineLabelPerimSpaceF = 0.00 ; Essentially, limit the margins surrounding each label bo wsres@cnLineLabelPlacementMode = "Computed" wsres@cnLineDrawOrder = "Draw" wsres@cnLabelDrawOrder = "Draw" ;*********************************************************************************************************** ; Attributes for wind barbs on the dynamic tropopause ;*********************************************************************************************************** wbres = True ; Resources for PWAT desired wbres@gsnDraw = False ; Do not draw plot until specified otherwise wbres@gsnFrame = False ; Do not advance frame until specified otherwise wbres@gsnAddCyclic = True wbres@gsnLeftString = "" wbres@gsnRightString = "" wbres@vcGlyphStyle = "WindBarb" wbres@vcWindBarbColor = (/0.25,0.25,0.25/) wbres@vcWindBarbCalmCircleSizeF = 0 wbres@vcWindBarbLineThicknessF = 6.0 wbres@vcRefLengthF = 0.024 wbres@vcRefMagnitudeF = 10.0 wbres@vcRefAnnoOn = False wbres@vcRefAnnoString2On = False wbres@vcMinDistanceF = 0.03 wbres@vcRefAnnoOrthogonalPosF = -1.08 wbres@vcWindBarbTickLengthF = 0.37 wbres@vcWindBarbTickSpacingF = 0.15 wbres@vcWindBarbTickAngleF = 55.0 wbres@vcWindBarbScaleFactorF = 2.0 ;************************************************ ; create plots ;************************************************ ;; TPV Tracks ypts1 = (/72.0,69.0,67.0,66.0,65.5,66.5,67.0,67.5,68.5,69.0,70.5,71.0,73.5,73.0,72.0,71.5,\ 71.5,71.5,72.0,72.0,73.0,73.5,76.5,78.0,78.5,81.0,81.5,80.0,77.5,76.0,74.0,72.0,\ 70.5,69.5,69.0,69.5,70.5,71.5,70.5,69.0,68.5,67.0,69.0,68.0,68.0,67.0,65.5,64.0,\ 63.5,63.5,63.5,63.5,62.5,62.0,62.0,62.0,62.0,60.0,59.5,56.5,55.5,56.5,57.5,59.0,\ 60.5,60.0,60.0,61.0,62.5,64.0,66.0,67.0,67.0,66.5,66.0,66.5,68.0,69.0,70.5,70.5,\ 68.5,68.0,68.5,69.0,70.5,71.0,71.0,70.0,68.5,67.0,65.5,64.5,64.5,65.0,64.0,64.5,\ 64.0,63.5,61.5,58.5,56.5,53.0,49.5,46.0,42.0,39.0,39.0,42.5,48.0,50.0,51.0,59.0,\ 63.5,64.5,64.0,65.5/) xpts1 = (/288.0,290.0,290.0,294.0,301.5,306.0,308.0,313.5,308.0,324.0,328.0,332.0,313.0,\ 310.5,306.0,308.0,310.0,314.0,318.0,320.0,324.0,326.0,323.5,321.0,318.0,308.0,\ 284.0,270.0,268.0,267.5,261.0,263.0,266.5,268.5,273.5,277.5,277.5,274.0,270.0,\ 268.5,267.0,268.5,263.0,261.5,263.5,260.0,258.0,261.0,263.5,263.5,262.0,260.5,\ 259.5,259.0,259.0,258.5,257.5,259.0,257.5,262.5,270.0,273.5,276.0,277.5,270.5,\ 270.0,270.0,272.0,273.0,273.0,268.5,265.0,260.0,259.5,260.5,262.5,263.0,263.0,\ 258.5,253.0,250.0,251.5,253.5,253.5,251.5,248.0,244.0,242.0,243.5,245.0,247.5,\ 252.0,256.5,259.5,260.5,261.0,261.0,259.5,258.5,260.0,260.0,260.0,262.0,265.0,\ 270.0,280.0,283.0,292.5,295.5,297.0,300.0,317.5,319.5,319.5,322.5,325.5/) ;; Cold Pool Track ypts2 = (/82.0,82.0,82.0,81.5,79.5,79.0,79.0,77.0,75.5,73.5,71.5,69.5,69.0,69.5,70.5,70.5,\ 70.0,69.0,69.5,69.0,68.5,67.5,66.5,65.5,64.5,64.5,63.5,63.5,62.5,62.0,62.0,62.0,\ 62.0,62.5,63.0,58.0,57.5,57.0,57.0,57.0,58.0,59.0,60.0,60.5,61.0,62.5,63.5,64.5,\ 66.0,66.5,66.5,66.5,66.5,67.0,67.0,68.0,69.0,69.5,70.0,70.0,70.0,69.5,70.0,70.5,\ 70.5,70.0,69.0,68.0,66.0,65.0,63.5,63.0,61.5,61.5,61.5,63.0,62.5,60.5,57.5,54.0,\ 50.0,46.5,43.5,45.0,45.0,46.5,49.5,50.0,48.5,47.0,46.0,45.5,45.5,45.0,45.0,48.5,49.5/) xpts2 = (/309.5,305.5,304.0,299.5,289.5,284.5,273.5,267.5,267.0,263.5,262.5,266.5,271.5,275.0,\ 274.0,272.5,270.5,269.0,267.5,265.5,264.0,262.5,262.0,261.5,262.0,261.5,262.5,260.0,\ 259.5,260.5,261.5,261.5,261.5,261.0,259.5,262.5,264.5,266.0,268.5,271.0,272.5,273.0,\ 272.0,271.5,271.0,271.5,271.0,269.5,266.0,264.0,261.0,260.0,259.5,259.5,259.0,258.0,\ 256.5,255.0,253.5,252.0,250.5,250.5,249.5,249.0,246.5,244.0,243.5,244.5,246.0,249.5,\ 252.0,256.0,258.5,262.0,264.0,263.0,261.0,259.5,260.5,263.0,265.5,268.5,271.5,282.0,\ 285.5,289.5,291.0,288.0,290.0,290.0,290.5,294.5,297.5,300.0,303.5,309.5,312.0/) do i=0, ntimes-1 ; Begin a do loop that creates a plot for all times in the file, one by one. ; i is an arbitrary counter variable wks_type = "png" wks_type@wkWidth = 2500 wks_type@wkHeight = 2500 j = i+0 imageName = "DT_theta_tracks_circle_reducedCB_400km_paperRes_redBlue_" + ut_string(times(j),"%Y%N%D_%H%M") wks = gsn_open_wks(wks_type, imageName) ; open a ps file with a unique name according to the time/date cmap = (/(/1.00,1.00,1.00/),\ (/0.00,0.00,0.00/)/) gsn_define_colormap(wks, cmap) ; Define Color Table res@gsnStringFont = "Helvetica" res@gsnStringFont = "Helvetica" res@gsnLeftString = "" ; cd_string uses cd_calendar to initially convert the time format ; such that nicely formatted strings van be customized ; %c is an informat character meaning small month abbreviation (e.g., Mar) ; An example output string is "1200 UTC 14 Mar 1993" res@gsnRightString = "" res@gsnLeftStringFontHeightF = 0.011 res@gsnRightStringFontHeightF = 0.009 ;res@gsnCenterString = date_str(i) ; Change title to reflect the time/date in file plot_map = gsn_csm_map(wks,res) plot_theta = gsn_csm_contour(wks,theta_pv(i,:,:),thtres) ;plot_ws = gsn_csm_contour(wks,ws(i,:,:),wsres) wbres@vcWindBarbTickAngleF = abs(wbres@vcWindBarbTickAngleF) plot_nh_wb = gsn_csm_vector(wks,u_pv(i,{0:90},:),v_pv(i,{0:90},:),wbres) overlay(plot_map,plot_theta) ;overlay(plot_map,plot_ws) overlay(plot_map,plot_nh_wb) earth_rad = 6371.0 rad_to_deg = 57.29578 k=i+67 ;if (k.gt.0) then if(k.le.115) rad_km_to_deg = (400./earth_rad)*rad_to_deg circ_lat = new (100,float) circ_lon = new (100,float) nggcog(ypts1(k),xpts1(k),rad_km_to_deg,circ_lat,circ_lon) rescl = True ; polyline mods desired ;rescl@gsLineColor = (/0.05,0.38,0.02/) ; color of lines ;rescl@gsLineColor = (/0.65,0.00,0.00/) ; color of lines ;rescl@gsLineColor = (/0.00,0.00,0.00/) ; color of lines rescl@gsLineThicknessF = 30.0 ; thickness of lines cpolyres = True ; polygon mods desired ;cpolyres@gsFillColor = (/0.58,0.88,0.54/) cpolyres@gsFillColor = (/1.00,0.40,0.37/) cpolyres@gsFillColor = (/0.5,0.5,0.5/) cpolyres@gsFillOpacityF = 0.4 ;dum_circ_1 = gsn_add_polygon(wks,plot_map,circ_lon,circ_lat,cpolyres) ; add polygon dum_circ_2 = gsn_add_polyline(wks,plot_map,circ_lon,circ_lat,rescl) ; add polygon end if print(i) ;k=i+0 k=i+67 ;if (k.gt.0) then ;if (k.gt.67) then n_gc_pts_1 = k ypts1_gc = new((/n_gc_pts_1,2/),"float") do j = 0, k-1 ypts1_gc(j,0) = ypts1(j) ypts1_gc(j,1) = ypts1(j+1) end do xpts1_gc = new((/n_gc_pts_1,2/),"float") do j = 0, k-1 xpts1_gc(j,0) = xpts1(j) xpts1_gc(j,1) = xpts1(j+1) end do npts = 5 ; number of points to interpolate along great circle path between two points ;---Define array to hold box. trackLat1 = new(n_gc_pts_1*npts,float) trackLon1 = new(n_gc_pts_1*npts,float) do j=0, n_gc_pts_1-1 jbeg = j*npts jend = jbeg+npts-1 dist = gc_latlon(ypts1_gc(j,0),xpts1_gc(j,0), \ ypts1_gc(j,1),xpts1_gc(j,1),npts,2) trackLat1(jbeg:jend) = dist@gclat trackLon1(jbeg:jend) = dist@gclon end do resp = True ; polyline mods desired resp@gsLineColor = (/0.00,0.00,0.00/) ; color of lines resp@gsLineThicknessF = 45.0 ; thickness of lines resp@gsLineLabelString= "" ; adds a line label string dumTrackBlack1 = gsn_add_polyline(wks, plot_map, trackLon1, trackLat1, resp) resp@gsLineColor = (/1.00,0.23,0.22/) resp@gsLineThicknessF = 23.0 dumTrackColor1 = gsn_add_polyline(wks, plot_map, trackLon1, trackLat1, resp) delete(ypts1_gc) delete(xpts1_gc) delete(trackLat1) delete(trackLon1) ;end if if (i.le.48) then dum_mk_1 = new(1,graphic) dum_mk_hollow_1 = new(1,graphic) mk_res = True mk_res@gsMarkerIndex = 16 ; Solid Circle mk_res@gsMarkerSizeF = 16 mk_res@gsMarkerColor = (/1.00,0.00,0.00/) dum_mk_1(0)=gsn_add_polymarker(wks,plot_map,xpts1(k),ypts1(k),mk_res) mk_res@gsMarkerIndex = 4 ; Hollow dots mk_res@gsMarkerThicknessF = 14.0 mk_res@gsMarkerColor = (/0.00,0.00,0.00/) dum_mk_hollow_1(0)=gsn_add_polymarker(wks,plot_map,xpts1(k),ypts1(k),mk_res) end if ;end if ;if (i.gt.22) then ; k=i-23 ; if (k.gt.0) then k=i+45 ;if (k.gt.45) then n_gc_pts_2 = k ypts2_gc = new((/n_gc_pts_2,2/),"float") do j = 0, k-1 ypts2_gc(j,0) = ypts2(j) ypts2_gc(j,1) = ypts2(j+1) end do xpts2_gc = new((/n_gc_pts_2,2/),"float") do j = 0, k-1 xpts2_gc(j,0) = xpts2(j) xpts2_gc(j,1) = xpts2(j+1) end do npts = 5 ; number of points to interpolate along great circle path between two points ;---Define array to hold box. trackLat2 = new(n_gc_pts_2*npts,float) trackLon2 = new(n_gc_pts_2*npts,float) do j=0, n_gc_pts_2-1 jbeg = j*npts jend = jbeg+npts-1 dist = gc_latlon(ypts2_gc(j,0),xpts2_gc(j,0), \ ypts2_gc(j,1),xpts2_gc(j,1),npts,2) trackLat2(jbeg:jend) = dist@gclat trackLon2(jbeg:jend) = dist@gclon end do resp = True ; polyline mods desired resp@gsLineColor = (/0.00,0.00,0.00/) ; color of lines resp@gsLineThicknessF = 45.0 ; thickness of lines resp@gsLineLabelString= "" ; adds a line label string dumTrackBlack2 = gsn_add_polyline(wks, plot_map, trackLon2, trackLat2, resp) resp@gsLineColor = (/0.02,0.54,1.00/) resp@gsLineThicknessF = 23.0 dumTrackColor2 = gsn_add_polyline(wks, plot_map, trackLon2, trackLat2, resp) delete(ypts2_gc) delete(xpts2_gc) delete(trackLat2) delete(trackLon2) ;end if dum_mk_2 = new(1,graphic) dum_mk_hollow_2 = new(1,graphic) mk_res = True mk_res@gsMarkerIndex = 16 ; Solid Circle mk_res@gsMarkerSizeF = 16 mk_res@gsMarkerColor = (/0.00,0.00,1.00/) dum_mk_2(0)=gsn_add_polymarker(wks,plot_map,xpts2(k),ypts2(k),mk_res) mk_res@gsMarkerIndex = 4 ; Hollow dots mk_res@gsMarkerThicknessF = 14.0 mk_res@gsMarkerColor = (/0.00,0.00,0.00/) dum_mk_hollow_2(0)=gsn_add_polymarker(wks,plot_map,xpts2(k),ypts2(k),mk_res) ;end if res@gsnDraw = False ; Do not draw plot until specified otherwise res@gsnFrame = False ; Do not advance frame until specified otherwise maximize_output(wks,res) draw(plot_map) frame(wks) system("convert -trim +repage "+imageName+".png "+ \ imageDirectory+"/"+imageName+".png") system("rm *.png") end do end