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