1.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     1.2 +++ b/lai/32.contour_diff.ncl	Mon Jan 26 22:08:20 2009 -0500
     1.3 @@ -0,0 +1,87 @@
     1.4 +;*************************************************
     1.5 +; ce_1.ncl
     1.6 +;************************************************
     1.7 +load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"  
     1.8 +load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"  
     1.9 +;************************************************
    1.10 +begin
    1.11 +;************************************************
    1.12 +; read in observed data
    1.13 +;************************************************
    1.14 + f     = addfile ("Npp_T42_mean.nc","r")
    1.15 + y     = f->NPP
    1.16 +;************************************************
    1.17 +; read in model data
    1.18 +;************************************************
    1.19 + g     = addfile ("i01.03cn_1545-1569_ANN_climo.nc","r")
    1.20 +;g     = addfile ("i01.04casa_1605-1629_ANN_climo.nc","r")
    1.21 + x     = g->NPP
    1.22 +
    1.23 + nsec_per_year = 86400.*365.
    1.24 + x = x * nsec_per_year
    1.25 + 
    1.26 + x@units = "gC/m^2/year" 
    1.27 +;************************************************
    1.28 +;fill ob missing grid the same as model
    1.29 +;************************************************
    1.30 + delta = 0.00001
    1.31 + x0    = x(0,:,:)
    1.32 + y = where(ismissing(y).and.(ismissing(x0).or.(x0.lt.delta)),0.,y)
    1.33 +
    1.34 + y@units = "gC/m^2/year"
    1.35 +;************************************************
    1.36 +;model - observed
    1.37 +;************************************************ 
    1.38 + z = x
    1.39 + z = x(0,:,:) - y(:,:)
    1.40 +;************************************************
    1.41 +; create 3 plots
    1.42 +;************************************************
    1.43 +  wks = gsn_open_wks("ps","xy")             ; open a ps file
    1.44 +  gsn_define_colormap(wks,"gui_default")     ; choose colormap
    1.45 +
    1.46 +  res                     = True             ; Use plot options
    1.47 +  res@cnFillOn            = True             ; Turn on color fill
    1.48 +  res@gsnSpreadColors      = True            ; use full colormap
    1.49 +; res@cnFillMode          = "RasterFill"     ; Turn on raster color
    1.50 +; res@lbLabelAutoStride   = True
    1.51 +  res@cnLinesOn           = False            ; Turn off contourn lines
    1.52 +  res@mpFillOn            = False            ; Turn off map fill
    1.53 +  res@tiMainString        = "Observed MODIS MOD 17"
    1.54 +
    1.55 +  res@gsnSpreadColors     = True             ; use full colormap
    1.56 +  res@cnLevelSelectionMode = "ManualLevels"  ; Manual contour invtervals
    1.57 +  res@cnMinLevelValF       = 0.              ; Min level
    1.58 +  res@cnMaxLevelValF       = 2200.           ; Max level
    1.59 +  res@cnLevelSpacingF      = 200.            ; interval
    1.60 +
    1.61 +  res@gsnFrame             = False           ; Do not draw plot 
    1.62 +  res@gsnDraw              = False           ; Do not advance frame
    1.63 +
    1.64 +  plot=new(3,graphic)                        ; create graphic array
    1.65 +
    1.66 +  plot(0) = gsn_csm_contour_map_ce(wks,y,res)        ; for observed
    1.67 +
    1.68 +  res@tiMainString        = "Model i01.03cn"
    1.69 +; res@tiMainString        = "Model i01.04casa"
    1.70 +  plot(1) = gsn_csm_contour_map_ce(wks,x(0,:,:),res) ; for model
    1.71 +
    1.72 +  res@cnMinLevelValF       = -500           ; Min level
    1.73 +  res@cnMaxLevelValF       =  500.          ; Max level
    1.74 +  res@cnLevelSpacingF      = 50.            ; interval
    1.75 +  res@tiMainString        = "(Model i01.03cn) - (observed)"
    1.76 +; res@tiMainString        = "(Model i01.04casa) - (observed)"
    1.77 +  plot(2) = gsn_csm_contour_map_ce(wks,z(0,:,:),res) ; for model - ob
    1.78 +
    1.79 +;***********************************************
    1.80 +; create panel plot
    1.81 +;***********************************************
    1.82 +  pres                            = True        ; panel plot mods desired
    1.83 +  pres@gsnPanelYWhiteSpacePercent = 5           ; increase white space around
    1.84 +                                                ; indiv. plots in panel
    1.85 +  pres@gsnMaximize                = True        ; fill the page
    1.86 +
    1.87 +  gsn_panel(wks,plot,(/3,1/),pres)              ; create panel plot
    1.88 +
    1.89 +  system("convert xy.ps xy.png")
    1.90 +end
    1.91 \ No newline at end of file