1.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     1.2 +++ b/lai/33.contour_diff_phase.ncl	Mon Jan 26 22:08:20 2009 -0500
     1.3 @@ -0,0 +1,116 @@
     1.4 +;*************************************************
     1.5 +; 3 plots, model - ob
     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 +load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"  
    1.10 +;************************************************
    1.11 +begin
    1.12 +;************************************************
    1.13 +; read in observed data
    1.14 +;************************************************
    1.15 +  diri1  = "/fis/cgd/cseg/people/jeff/clamp_data/lai/"
    1.16 +  fili1  = "LAI_2000-2005_ensemble_T42.nc"
    1.17 +  f     = addfile(diri1+fili1,"r")
    1.18 +
    1.19 +  z     = f->LAI
    1.20 +  y     = z(0,:,:)
    1.21 +  y@long_name = "Month of Leaf Area Index Max"
    1.22 +  s     = z(:,0,0)
    1.23 + 
    1.24 +  dsizes_z = dimsizes(z)
    1.25 +  ntime    = dsizes_z(0)
    1.26 +  nlat     = dsizes_z(1)
    1.27 +  nlon     = dsizes_z(2)
    1.28 +  
    1.29 +  do j = 0,nlat-1
    1.30 +  do i = 0,nlon-1
    1.31 +     s      = z(:,j,i) 
    1.32 +     y(j,i) = maxind(s) + 1
    1.33 +  end do
    1.34 +  end do
    1.35 +
    1.36 +; printVarSummary(y)
    1.37 +  print (min(y)+"/"+max(y))
    1.38 +
    1.39 +  delete (s)
    1.40 +  delete (z)
    1.41 +;************************************************
    1.42 +; read in data: model       
    1.43 +;************************************************
    1.44 + diri2  = "/fis/cgd/cseg/people/jeff/clamp_data/model/"
    1.45 + fili3  = "i01.03cn_1545-1569_MONS_climo.nc"
    1.46 +;fili3  = "i01.04casa_1605-1629_MONS_climo.nc"
    1.47 + data_file_model = addfile(diri2+fili3,"r")
    1.48 +
    1.49 +  z     = data_file_model->TLAI
    1.50 +  x     = z(0,:,:)
    1.51 +  x@long_name = "Month of Leaf Area Index Max"
    1.52 +  s     = z(:,0,0)  
    1.53 + 
    1.54 +  dsizes_z = dimsizes(z)
    1.55 +  ntime    = dsizes_z(0)
    1.56 +  nlat     = dsizes_z(1)
    1.57 +  nlon     = dsizes_z(2)
    1.58 +  
    1.59 +  do j = 0,nlat-1
    1.60 +  do i = 0,nlon-1
    1.61 +     s      = z(:,j,i) 
    1.62 +     x(j,i) = maxind(s) + 1
    1.63 +  end do
    1.64 +  end do
    1.65 +
    1.66 +  print (min(x)+"/"+max(x))
    1.67 +  d = x              
    1.68 +  d = x - y
    1.69 +  d = where((d .gt. 6.),12.-d,d)
    1.70 +  d = where((d .lt. -6.),-12.-d,d)
    1.71 +  print (min(d)+"/"+max(d))
    1.72 +
    1.73 +  delete (z)
    1.74 +  delete (s)
    1.75 +
    1.76 +  delta = 0.000001
    1.77 +  y = where(ismissing(y).and.(ismissing(x).or.(x.lt.delta)),0.,y)
    1.78 +;************************************************
    1.79 +; create default plot
    1.80 +;************************************************
    1.81 +  wks = gsn_open_wks("ps","xy")              ; open a ps file
    1.82 +  gsn_define_colormap(wks,"gui_default")     ; choose colormap
    1.83 +
    1.84 +  res                     = True             ; Use plot options
    1.85 +  res@cnFillOn            = True             ; Turn on color fill
    1.86 +  res@gsnSpreadColors     = True             ; use full colormap
    1.87 +; res@cnFillMode          = "RasterFill"     ; Turn on raster color
    1.88 +; res@lbLabelAutoStride   = True
    1.89 +  res@cnLinesOn           = False            ; Turn off contourn lines
    1.90 +  res@mpFillOn            = False            ; Turn off map fill
    1.91 +
    1.92 +  res@gsnSpreadColors      = True            ; use full colormap
    1.93 +  res@cnLevelSelectionMode = "ManualLevels"  ; Manual contour invtervals
    1.94 +  res@cnMinLevelValF       = 1.              ; Min level
    1.95 +  res@cnMaxLevelValF       = 12.             ; Max level
    1.96 +  res@cnLevelSpacingF      = 1.              ; interval
    1.97 +
    1.98 +  pres                     = True            ; panel plot mods desired
    1.99 +  pres@gsnMaximize         = True            ; fill the page
   1.100 +
   1.101 +  plot=new(3,graphic)                        ; create graphic array
   1.102 +
   1.103 +  res@tiMainString        = "MODIS MOD 15A2 2000-2005"
   1.104 +  plot(0) = gsn_csm_contour_map_ce(wks,y,res)   
   1.105 +
   1.106 +  res@tiMainString        = "Model i01.03cn"
   1.107 +; res@tiMainString        = "Model i01.04casa"
   1.108 +  plot(1) = gsn_csm_contour_map_ce(wks,x,res)   
   1.109 + 
   1.110 +  res@tiMainString        = "(Model i01.03cn) - (Observed)"
   1.111 +; res@tiMainString        = "(Model i01.04casa) - (Observed)"    
   1.112 +  res@cnMinLevelValF       = -6.              ; Min level
   1.113 +  res@cnMaxLevelValF       =  6.             ; Max level
   1.114 +  res@cnLevelSpacingF      =  1.              ; interval
   1.115 +  plot(2) = gsn_csm_contour_map_ce(wks,d,res)   ; for observed
   1.116 +
   1.117 +  gsn_panel(wks,plot,(/3,1/),pres)           ; create panel plot
   1.118 +  system("convert xy.ps xy.png")
   1.119 +end
   1.120 \ No newline at end of file