;;; Program to compute composite climatologies for Q-vector, height tendency, and ;;; topographic forcing of Bay Area storms. ;;; ;;; Philip Martin-Edwards ;;; METR 206 Final Project load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl" load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl" load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl" begin ;;; Read files and variables wks = gsn_open_wks("ps", "PMartinFinalProj") ; open workstation gsn_define_colormap(wks,"gui_default") ; choose colormap fils = systemfunc ("ls storm*.nc") ; list of file paths from this directory f = addfiles (fils, "r") ft = addfile ("/home/pmartin/metr206/finalproj/terrain.nc","r") ListSetType (f, "cat") ; *concatenates elements from leftmost dimension (time) from T = short2flt(f[:]->t) ; all files, all other dimensions remain the same phi =(short2flt(f[:]->z)) / 9.8 ; variable(time,level,lat,lon) U = short2flt(f[:]->u) ; *short2flt transforms short type to floating decimal and V = short2flt(f[:]->v) ; applies offsets, etc. rvo = short2flt(f[:]->vo) H = (short2flt(ft->z)) / 9.8 ;convert to meters lat = f[1]->latitude ; time invariant lon = f[1]->longitude levs = f[1]->level times = f[:]->time yyyymmddhh = cd_calendar(times, -3) ; print(times) ; printVarSummary (H) wspd925 = (sqrt(U(:,3,:,:)^2 + V(:,3,:,:)^2)) * 2.23694 ;windspeed mph wspd925@units = "mph" copy_VarCoords(U(:,3,:,:),wspd925) copy_VarCoords(T,phi) ; printVarSummary(times) ;;; Constants nlev = dimsizes(levs) nlat = dimsizes(lat) nlon = dimsizes(lon) dlat = 0.25 * 0.0174533 ;in rads dlon = 0.25 * 0.0174533 a = 6.378e6 f0 = 0.001 p0 = 100000. ; in Pa R = 287.058 P = 70000. omega = 7.292e-5 pvo = 2.*omega*sin(lat) ;planetary vorticity te0inds = (/ind(yyyymmddhh.eq.1982033112),ind(yyyymmddhh.eq.1982122218),ind(yyyymmddhh.eq.1986021418), \ ind(yyyymmddhh.eq.1986021800),ind(yyyymmddhh.eq.1995121212), ind(yyyymmddhh.eq.2000021412), \ ind(yyyymmddhh.eq.2002121606),ind(yyyymmddhh.eq.2005123112), ind(yyyymmddhh.eq.2008010418), \ ind(yyyymmddhh.eq.2009101318)/) tm36inds = te0inds - 6 dy = 2. * a * dlat ;dx will require a loop ;;; Plot 925 heights and wind speed for T=0 plot = new(10,graphic) ; create graphic array res = True res@mpDataBaseVersion = "Ncarg4_1" ; choose more recent database ; res@mpDataSetName = "Earth..4" ; high resolution res@mpOutlineOn = True res@mpOutlineBoundarySets = "National" ; National borders res@mpNationalLineColor = "Black" ; National borders color res@mpGeophysicalLineColor = "Black" res@cnFillOn = False ; turn on color fill res@gsnDraw = False ; do not draw picture res@gsnFrame = False ; do not advance frame res@lbOrientation = "Vertical" ; vertical label bar res@gsnAddCyclic = False res@mpGeophysicalLineThicknessF = 2.0 res@mpMinLatF = min(lat) ; range to zoom in on res@mpMaxLatF = max(lat) res@mpMinLonF = min(lon) res@mpMaxLonF = max(lon) wres = True wres@cnFillOn = True ; wres@cnFillDrawOrder = "PreDraw" wres@cnLinesOn = False wres@gsnAddCyclic = False wres@gsnDraw = False ; do not draw picture wres@gsnFrame = False ; do not advance frame wres@cnFillOpacityF = 0.4 vres = True vres@gsnDraw = False ; don't draw vres@gsnFrame = False ; don't advance frame vres@vcGlyphStyle = "CurlyVector" ; curly vectors vres@vcRefMagnitudeF = 20 ; define vector ref mag vres@vcRefLengthF = 0.045 ; define length of vec ref vres@gsnRightString = " " ; turn off right string vres@gsnLeftString = " " ; turn off left string vres@tiXAxisString = " " ; turn off axis label do i=0,dimsizes(te0inds)-1 plotA = gsn_csm_contour_map(wks,phi(te0inds(i),3,:,:),res) plotB = gsn_csm_contour(wks,wspd925(te0inds(i),:,:),wres) overlay(plotA,plotB) ; result will be plotA plot(i) = plotA delete(plotA) ;otherwise all 10 will overlay on each other delete(plotB) end do ;;; Compute composites of geopotential height, vorticity, temperature, and wind fields @ ;;; 500 hPa, 700 hPa, and 925 hPa across all 10 storms at time intervals T-36, T-24, T-12, ;;; T=0, T+12. ;;; while loop n=0 do while(n.le.4) ;do n=0,4 cT = dim_avg_n_Wrap(T(tm36inds,:,:,:),0) ;computes average over time dimension cphi = dim_avg_n_Wrap(phi(tm36inds,:,:,:),0) ;using only indices of given storm cU = dim_avg_n_Wrap(U(tm36inds,:,:,:),0) ;intervals, starting at T-36hr cV = dim_avg_n_Wrap(V(tm36inds,:,:,:),0) ;result(lev,lat,lon) crvo = dim_avg_n_Wrap(rvo(tm36inds,:,:,:),0) cavo = new((/nlev,nlat,nlon/),"float") do i=0, nlat-1 cavo(:,i:i,:) = crvo(:,i:i,:) + pvo(i) end do copy_VarCoords(crvo,cavo) ;printVarSummary(cavo) ;print(cavo(:,100,100)) ctheta = cT*((p0 / P)^(0.286)) ;;; Compute Q-vector components + divergence ; center_finite_diff_n(u,x,True=cyclic False=one-sided bounds,unused,difference this u dim) dudx = new ((/nlev,nlat,nlon/), typeof(cU)) dvdx = new((/nlev,nlat,nlon/), typeof(cV)) dthetadx = new((/nlev,nlat,nlon/), typeof(ctheta)) davodx = new((/nlev,nlat,nlon/), typeof(cavo)) do nl=0,nlat-1 ;dx changes depending on lat dx = a * cos(0.0174533 * lat(nl)) * dlon dudx(:,nl:nl,:) = (center_finite_diff_n(cU(:,nl:nl,:),dx,False,0,2)) ;NCL automatically truncates single-element dimensions dvdx(:,nl:nl,:) = center_finite_diff_n(cV(:,nl:nl,:),dx,False,0,2) ;when performing calculations and loops, dthetadx(:,nl:nl,:) = center_finite_diff_n(ctheta(:,nl:nl,:),dx,False,0,2) ;"nl:nl" preserves the dimension davodx(:,nl:nl,:) = center_finite_diff_n(cavo(:,nl:nl,:),dx,False,0,2) end do dudy = (center_finite_diff_n(cU,lat,False,0,1)) / dy dvdy = (center_finite_diff_n(cV,lat,False,0,1)) / dy dthetady = (center_finite_diff_n(ctheta,lat,False,0,1)) / dy davody = (center_finite_diff_n(cavo,lat,False,0,1)) / dy q1 = (-1.)*(R/P)*((dudx(1,:,:) * dthetadx(1,:,:)) + (dvdx(1,:,:) * dthetady(1,:,:))) ;(nlat,nlon) q2 = (-1.)*(R/P)*((dudy(1,:,:) * dthetadx(1,:,:)) + (dvdy(1,:,:) * dthetady(1,:,:))) dq1dx = new((/nlat,nlon/), typeof(q1)) dHdx = new((/nlat,nlon/), typeof(H)) do nl=0,nlat-1 ;dx changes depending on lat dx = a * cos(0.0174533 * lat(nl)) * dlon dq1dx(nl:nl,:) = center_finite_diff_n(q1(nl:nl,:),dx,False,0,1) dHdx(nl:nl,:) = center_finite_diff_n(H(0,nl:nl,:),dx,False,0,1) end do dq2dy = (center_finite_diff_n(q2,lat,False,0,0)) / dy divQ = dq1dx + dq2dy copy_VarCoords(U(1,1,:,:),divQ) ;printVarSummary(divQ) ;;; Compute height tendency voadv850 = (-1.)*f0*((cU(2,:,:) * davodx(2,:,:)) + (cV(2,:,:) * davody(2,:,:))) ;printVarSummary(voadv850) Tadv700 = (-1.)*(cU(1,:,:) * dthetadx(1,:,:)) - (cV(1,:,:) * dthetady(1,:,:)) Tadv925 = (-1.)*(cU(3,:,:) * dthetadx(3,:,:)) - (cV(3,:,:) * dthetady(3,:,:)) dTadv850 = (Tadv700 - Tadv925) / 22500. chi850 = voadv850 - dTadv850 copy_VarCoords(divQ,chi850) ;printVarSummary(chi850) ;;; Compute terrain forcing dHdy = (center_finite_diff_n(H(0,:,:),lat,False,0,0)) / dy Wtopo = (-1.)*((cU(3,:,:) * dHdx) + (cV(3,:,:) * dHdy)) copy_VarCoords(divQ,Wtopo) ;printVarSummary(Wtopo) ;;; plot stuff c500plt1 = new(5,graphic) c500plt2 = new(5,graphic) c700plt1 = new(5,graphic) c700plt2 = new(5,graphic) c925plt1 = new(5,graphic) c925plt2 = new(5,graphic) cdivQplt = new(5,graphic) cchiplt = new(5,graphic) cWtopoplt= new(5,graphic) cres=True cres@gsnDraw = False cres@gsnFrame = False cres@gsnAddCyclic = False cres@mpMinLatF = min(lat) ; range to zoom in on cres@mpMaxLatF = max(lat) cres@mpMinLonF = min(lon) cres@mpMaxLonF = max(lon) plotC = gsn_csm_contour_map(wks,cphi(0,:,:),cres) plotD = gsn_csm_contour(wks,cavo(0,:,:),wres) overlay(plotC, plotD) c500plt1(n) = plotC delete(plotC) delete(plotD) plotE = gsn_csm_contour_map(wks,cT(0,:,:),res) plotF = gsn_csm_vector(wks,cU(0,::4,::4),cV(0,::4,::4),vres) overlay(plotE, plotF) c500plt2(n) = plotE delete(plotE) delete(plotF) plotG = gsn_csm_contour_map(wks,cphi(1,:,:),res) plotH = gsn_csm_contour(wks,cavo(1,:,:),wres) overlay(plotG, plotH) c700plt1(n) = plotG delete(plotG) delete(plotH) plotI = gsn_csm_contour_map(wks,cT(1,:,:),res) plotJ = gsn_csm_vector(wks,cU(1,::4,::4),cV(1,::4,::4),vres) overlay(plotI, plotJ) c700plt2(n) = plotI delete(plotI) delete(plotJ) plotK = gsn_csm_contour_map(wks,cphi(3,:,:),res) plotL = gsn_csm_contour(wks,cavo(3,:,:),wres) overlay(plotK, plotL) c925plt1(n) = plotK delete(plotK) delete(plotL) plotM = gsn_csm_contour_map(wks,cT(3,:,:),res) plotN = gsn_csm_vector(wks,cU(3,::4,::4),cV(3,::4,::4),vres) overlay(plotM, plotN) c925plt2(n) = plotM delete(plotM) delete(plotN) plotO = gsn_csm_contour_map(wks,divQ(:,:),res) cdivQplt(n) = plotO delete(plotO) plotP = gsn_csm_contour_map(wks,chi850(:,:),res) cchiplt(n) = plotP delete(plotP) plotQ = gsn_csm_contour_map(wks,Wtopo(:,:),res) cWtopoplt(n) = plotQ delete(plotQ) delete([/cT,cphi,cU,cV,crvo,cavo,ctheta,dudx,dvdx,dthetadx,davodx,dudy,dvdy,dthetady,davody,q1,q2,\ dq1dx,dHdx,dq2dy,divQ,voadv850,Tadv700,Tadv925,dTadv850,chi850,dHdy,Wtopo/]) print(n) tm36inds = tm36inds + 2 ;advance to next time interval (+12 hrs) n=n+1 end do printVarSummary(c500plt1) ; draw panel with white space added resP = True ; resP@gsnPanelYWhiteSpacePercent = 5 ; resP@gsnPanelXWhiteSpacePercent = 5 ; resP@gsnMaximize = True ; resP@txString = "ERA-Interim Reanalysis for T=0, 925 hPa" gsn_panel(wks,plot,(/5,2/),resP) gsn_panel(wks,c500plt1,(/3,2/),False) gsn_panel(wks,c500plt2,(/3,2/),resP) gsn_panel(wks,c700plt1,(/3,2/),resP) gsn_panel(wks,c700plt2,(/3,2/),resP) gsn_panel(wks,c925plt1,(/3,2/),resP) gsn_panel(wks,c925plt2,(/3,2/),resP) gsn_panel(wks,cdivQplt,(/3,2/),resP) gsn_panel(wks,cchiplt,(/3,2/),resP) gsn_panel(wks,cWtopoplt,(/3,2/),resP) end