;load built-in functions 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/csm/shea_util.ncl" begin a1 = addfile("/ptmp/notaro/africa/ha.foam1-5lpj.ctl-pi-test.sfcU.6351034.nc","r") ; datafile 1 lat = a1->lat ; foam lat lon = a1->lon ; foam lon time = a1->time ; time b = addfile("/ptmp/notaro/africa/hl.foam1-5lpj.cpt.fpar.moni6351034.correct.nc","r") ; datafile 2 lat1 = b->lat ; lpj lat lon1 = b->lon ; lpj lon b1 = addfile("/ptmp/notaro/run/foam.oro.nc","r") ; read oro for foam orofoam = b1->ORO(3,:,:) orofoam!0="lat" orofoam!1="lon" orofoam&lat=lat orofoam&lon=lon b2 = addfile("/ptmp/notaro/run/hl.oro.nc","r") ; read oro for lpj orolpj = b2->ORO(3,:,:) orolpj!0="lat" orolpj!1="lon" orolpj&lat=lat1 orolpj&lon=lon1 ;;;;;;;;;;;;;;;;; speed = new((/4800,40,48/),float) ; 4800 months speed = (a1->U(:,0,:,:)) ; speed = U wind speed!0 = "time" speed!1 = "lat" speed!2 = "lon" speed&lat = lat speed&lon = lon speed@time = time speed@_FillValue = 1e+35 speed=rmMonAnnCycTLL(speed) ; remove annual cycle from monthly data anomspeed = new((/4800,40,48/),float) anomspeed = mask (speed, conform(speed, orofoam, (/1,2/)), 1) ; apply land mask copy_VarCoords(speed,anomspeed) veg = new((/4800,128,128/),float) veg = (b->fpar) ; veg = fpar veg!0 = "time" veg!1 = "lat" veg!2 = "lon" veg&lat = lat1 veg&lon = lon1 veg@time = time veg@_FillValue = 1e+35 veg=rmMonAnnCycTLL(veg) ; remove annual cycle from monthly data anomfor = new((/4800,128,128/),float) anomfor = mask (veg, conform(veg, orolpj, (/1,2/)), 1) ; apply land mask copy_VarCoords(veg,anomfor) ; Area av for 8-22N, 18W to 40E ts1=new((/4800/),float) ts2=new((/4800/),float) do it=0,4799 ts1(it)=(7.*avg(anomfor(it,69:79,121:127))+15.*avg(anomfor(it,69:79,0:14)))/22. ; ts1=fpar ts2(it)=(3.*avg(anomspeed(it,21:25,45:47))+6.*avg(anomspeed(it,21:25,0:5)))/9. ; ts2=u wind end do times=ispan(-10,10,1) correl=new((/21/),float) ; compute lead/lag correk between u and fpar do it=-10,0 correl(it+10)=escorc(ts1(0-it:4799),ts2(0:4799+it)) ; u leads end do do it=1,10 correl(it+10)=escorc(ts1(0:4799-it),ts2(0+it:4799)) ; fpar leads end do ;;;;;;;;;;;;;;; ;Auto-correlation calculations auto1=new((/21/),float) ; fpar auto2=new((/21/),float) ; speed do it=-10,0 auto1(it+10)=escorc(ts1(0-it:4799),ts1(0:4799+it)) auto2(it+10)=escorc(ts2(0-it:4799),ts2(0:4799+it)) end do do it=1,10 auto1(it+10)=escorc(ts1(0:4799-it),ts1(0+it:4799)) auto2(it+10)=escorc(ts2(0:4799-it),ts2(0+it:4799)) end do data=new((/3,21/),float) ; want to plot correl and two auto-correl, so assign to variable data data(0,:)=correl data(1,:)=auto1 data(2,:)=auto2 ;;;;;;;;;;;;;;; wks = gsn_open_wks("ps","autocorr_correl") ; open a ps file res = True ; plot mods desired res@gsnDraw = False ; don't draw res@gsnFrame = False ; don't advance frame res@tiMainString = "Lead/Lag Correls-N.Africa FPAR vs. Sfc U-Wind Mn Anom (8-22N,18W-40E)" res@tiMainFontHeightF = 0.011 res@trXMinF = -10.0 ; set min x value res@trXMaxF = 10.0 ; set max x value res@trYMinF = min(data) ; set min x value res@trYMaxF = max(data) ; set max x value res@tmXBMode = "Explicit" ; Define tick mark labels. res@tmXBValues = (/-10.,-8.,-6.,-4.,-2.,0.,2.,4.,6.,8.,10./) res@tmXBLabels = (/"-10","-8","-6","-4","-2","0","2","4","6","8","10"/) res@tmXBLabelFontHeightF = 0.01 ; resize tick labels res@tmYLLabelFontHeightF = 0.01 ; resize tick labels res@xyMarkLineModes = (/"MarkLines","MarkLines","MarkLines"/) ; line style res@xyLineThicknesses = (/6.,3.,3./) ; line thickness res@xyDashPatterns = (/0.,0.,1./) res@xyLineColors = (/"red","blue","green"/) res@tiXAxisString = "WIND Lead VEG Lead" res@tiYAxisString = "Correl Coeff" ; y-axis label txres = True ; Set up variable for TextItem resources. txres@txFontHeightF = 0.0055 ; Set the font height. res@tfPolyDrawOrder = "Predraw" ; put line on top res@pmLegendDisplayMode = "Always" ; turn on legend res@pmLegendZone = 0 res@pmLegendSide = "Top" ; Change location of res@lgJustification = "BottomRight" res@pmLegendOrthogonalPosF = 0.2 ; more neg = down res@pmLegendParallelPosF = -0.17 res@pmLegendWidthF = 0.18 ; Change width and res@pmLegendHeightF = 0.1 ; height of legend. res@lgLabelFontHeightF = .025 ; change font height res@lgPerimOn = True ; no box around res@xyExplicitLegendLabels = (/" FPAR,U"," FPAR"," U"/) plot = gsn_xy (wks,times,data,res) ; plots correl and auto-correls txres = True ; Set up variable for TextItem resources. txres@txFontHeightF = 0.015 gsres = True ; now draw significance lines xp = new( (/2/), float ) yp = new( (/2/), float ) xp = (/-10.,10./) yp = (/.04,.04/) ; if n=4800 months, need r=.04 for .90 signif gsres@gsLineColor = "brown" gsres@gsLineThicknessF = 2.2 dummy = gsn_add_polyline (wks,plot,xp,yp,gsres) xxp = new( (/2/), float ) yyp = new( (/2/), float ) xxp = (/-10.,10./) yyp = (/-.04,-.04/) ; if n=4800 months, need r=.04 for .90 signif dummy2 = gsn_add_polyline (wks,plot,xxp,yyp,gsres) xxxp = new((/2/),float) yyyp = new((/2/),float) xxxp=(/0.,0./) yyyp=(/-1.,1./) dummy3 = gsn_add_polyline (wks,plot,xxxp,yyyp,gsres) pres = True ; mod panel plot pres@gsnDraw=True pres@gsnFrame=True pres@gsnMaximize = True ; blow up plot pres@gsnPaperWidth=8.5 pres@gsnPaperHeight=11.0 pres@gsnPaperOrientation="portrait" pres@gsnPaperMargin=0.07 pres@gsnPanelLeft = 0.04 pres@gsnPanelRight = 0.96 pres@gsnPanelBottom = 0.04 pres@gsnPanelTop = 0.96 pres@gsnPanelYWhiteSpacePercent = 3 pres@gsnPanelXWhiteSpacePercent = 3 pres@gsnPanelRowSpec=True gsn_panel(wks,plot,(/1/),pres) ; create panel plot end