forrest@0: ;*************************************************
forrest@0: ; 3 plots, model - ob
forrest@0: ;************************************************
forrest@0: load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"  
forrest@0: load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
forrest@0: load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"  
forrest@0: ;************************************************
forrest@0: begin
forrest@0: ;************************************************
forrest@0: ; read in observed data
forrest@0: ;************************************************
forrest@0:   diri1  = "/fis/cgd/cseg/people/jeff/clamp_data/lai/"
forrest@0:   fili1  = "LAI_2000-2005_ensemble_T42.nc"
forrest@0:   f     = addfile(diri1+fili1,"r")
forrest@0: 
forrest@0:   z     = f->LAI
forrest@0:   y     = z(0,:,:)
forrest@0:   y@long_name = "Leaf Area Index Mean"
forrest@0:   s     = z(:,0,0)
forrest@0:  
forrest@0:   dsizes_z = dimsizes(z)
forrest@0:   ntime    = dsizes_z(0)
forrest@0:   nlat     = dsizes_z(1)
forrest@0:   nlon     = dsizes_z(2)
forrest@0:   
forrest@0:   do j = 0,nlat-1
forrest@0:   do i = 0,nlon-1
forrest@0:      s      = z(:,j,i) 
forrest@0:      y(j,i) = avg(s)
forrest@0:   end do
forrest@0:   end do
forrest@0: 
forrest@0: ; printVarSummary(y)
forrest@0:   print (min(y)+"/"+max(y))
forrest@0: 
forrest@0:   delete (s)
forrest@0:   delete (z)
forrest@0: ;************************************************
forrest@0: ; read in data: model       
forrest@0: ;************************************************
forrest@0:  diri2  = "/fis/cgd/cseg/people/jeff/clamp_data/model/"
forrest@0:  fili3  = "i01.03cn_1545-1569_MONS_climo.nc"
forrest@0: ;fili3  = "i01.04casa_1605-1629_MONS_climo.nc"
forrest@0:  data_file_model = addfile(diri2+fili3,"r")
forrest@0: 
forrest@0:   z     = data_file_model->TLAI
forrest@0:   x     = z(0,:,:)
forrest@0:   x@long_name = "Leaf Area Index Mean"
forrest@0:   s     = z(:,0,0)  
forrest@0:  
forrest@0:   dsizes_z = dimsizes(z)
forrest@0:   ntime    = dsizes_z(0)
forrest@0:   nlat     = dsizes_z(1)
forrest@0:   nlon     = dsizes_z(2)
forrest@0:   
forrest@0:   do j = 0,nlat-1
forrest@0:   do i = 0,nlon-1
forrest@0:      s      = z(:,j,i) 
forrest@0:      x(j,i) = avg(s)
forrest@0:   end do
forrest@0:   end do
forrest@0: 
forrest@0:   print (min(x)+"/"+max(x))
forrest@0:   d = x              
forrest@0:   d = x - y
forrest@0:   print (min(d)+"/"+max(d))
forrest@0: 
forrest@0:   delete (z)
forrest@0:   delete (s)
forrest@0: 
forrest@0:   delta = 0.000001
forrest@0:   y = where(ismissing(y).and.(ismissing(x).or.(x.lt.delta)),0.,y)
forrest@0: ;************************************************
forrest@0: ; create default plot
forrest@0: ;************************************************
forrest@0:   wks = gsn_open_wks("ps","xy")              ; open a ps file
forrest@0:   gsn_define_colormap(wks,"gui_default")     ; choose colormap
forrest@0: 
forrest@0:   res                     = True             ; Use plot options
forrest@0:   res@cnFillOn            = True             ; Turn on color fill
forrest@0:   res@gsnSpreadColors     = True             ; use full colormap
forrest@0: ; res@cnFillMode          = "RasterFill"     ; Turn on raster color
forrest@0: ; res@lbLabelAutoStride   = True
forrest@0:   res@cnLinesOn           = False            ; Turn off contourn lines
forrest@0:   res@mpFillOn            = False            ; Turn off map fill
forrest@0: 
forrest@0:   res@gsnSpreadColors      = True            ; use full colormap
forrest@0:   res@cnLevelSelectionMode = "ManualLevels"  ; Manual contour invtervals
forrest@0:   res@cnMinLevelValF       = 0.              ; Min level
forrest@0:   res@cnMaxLevelValF       = 6.             ; Max level
forrest@0:   res@cnLevelSpacingF      = 0.5             ; interval
forrest@0: 
forrest@0:   pres                     = True            ; panel plot mods desired
forrest@0:   pres@gsnMaximize         = True            ; fill the page
forrest@0: 
forrest@0:   plot=new(3,graphic)                        ; create graphic array
forrest@0: 
forrest@0:   res@tiMainString        = "MODIS MOD 15A2 2000-2005"
forrest@0:   plot(0) = gsn_csm_contour_map_ce(wks,y,res)   
forrest@0: 
forrest@0:   res@tiMainString        = "Model i01.03cn"
forrest@0: ; res@tiMainString        = "Model i01.04casa"
forrest@0:   plot(1) = gsn_csm_contour_map_ce(wks,x,res)   
forrest@0:  
forrest@0:   res@tiMainString        = "(Model i01.03cn) - (Observed)"
forrest@0: ; res@tiMainString        = "(Model i01.04casa) - (Observed)"    
forrest@0:   res@cnMinLevelValF       = -2.              ; Min level
forrest@0:   res@cnMaxLevelValF       =  2.             ; Max level
forrest@0:   res@cnLevelSpacingF      =  0.4              ; interval
forrest@0:   plot(2) = gsn_csm_contour_map_ce(wks,d,res)   ; for observed
forrest@0: 
forrest@0:   gsn_panel(wks,plot,(/3,1/),pres)           ; create panel plot
forrest@0:   system("convert xy.ps xy.png")
forrest@0: end