|
1 ;************************************************* |
|
2 ; 3 plots, model - ob |
|
3 ;************************************************ |
|
4 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl" |
|
5 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl" |
|
6 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl" |
|
7 ;************************************************ |
|
8 begin |
|
9 |
|
10 day_of_data = (/31,28,31,30,31,30,31,31,30,31,30,31/) |
|
11 |
|
12 ;************************************************ |
|
13 ; read in observed data |
|
14 ;************************************************ |
|
15 diri1 = "/fis/cgd/cseg/people/jeff/clamp_data/lai/" |
|
16 fili1 = "LAI_2000-2005_ensemble_T42.nc" |
|
17 f = addfile(diri1+fili1,"r") |
|
18 |
|
19 z = f->LAI |
|
20 y = z(0,:,:) |
|
21 y@long_name = "Days of Growing Season" |
|
22 |
|
23 dsizes_z = dimsizes(z) |
|
24 ntime = dsizes_z(0) |
|
25 nlat = dsizes_z(1) |
|
26 nlon = dsizes_z(2) |
|
27 |
|
28 do j = 0,nlat-1 |
|
29 do i = 0,nlon-1 |
|
30 nday = 0. |
|
31 do k = 0,ntime-1 |
|
32 if (.not. ismissing(z(k,j,i)) .and. z(k,j,i) .gt. 1.0) then |
|
33 nday = nday + day_of_data(k) |
|
34 end if |
|
35 end do |
|
36 if (nday .ne. 0.) then |
|
37 y(j,i) = nday |
|
38 end if |
|
39 end do |
|
40 end do |
|
41 |
|
42 ; printVarSummary(y) |
|
43 print (min(y)+"/"+max(y)) |
|
44 |
|
45 delete (z) |
|
46 ;************************************************ |
|
47 ; read in data: model |
|
48 ;************************************************ |
|
49 diri2 = "/fis/cgd/cseg/people/jeff/clamp_data/model/" |
|
50 ;fili3 = "i01.03cn_1545-1569_MONS_climo.nc" |
|
51 fili3 = "i01.04casa_1605-1629_MONS_climo.nc" |
|
52 data_file_model = addfile(diri2+fili3,"r") |
|
53 |
|
54 z = data_file_model->TLAI |
|
55 x = z(0,:,:) |
|
56 x@long_name = "Days of Growing Season" |
|
57 |
|
58 dsizes_z = dimsizes(z) |
|
59 ntime = dsizes_z(0) |
|
60 nlat = dsizes_z(1) |
|
61 nlon = dsizes_z(2) |
|
62 |
|
63 do j = 0,nlat-1 |
|
64 do i = 0,nlon-1 |
|
65 nday = 0. |
|
66 do k = 0,ntime-1 |
|
67 if (.not. ismissing(z(k,j,i)) .and. z(k,j,i) .gt. 1.0) then |
|
68 nday = nday + day_of_data(k) |
|
69 end if |
|
70 end do |
|
71 if (nday .ne. 0.) then |
|
72 x(j,i) = nday |
|
73 end if |
|
74 end do |
|
75 end do |
|
76 |
|
77 print (min(x)+"/"+max(x)) |
|
78 |
|
79 d = x |
|
80 d = x - y |
|
81 print (min(d)+"/"+max(d)) |
|
82 |
|
83 delete (z) |
|
84 |
|
85 delta = 0.000001 |
|
86 y = where(ismissing(y).and.(ismissing(x).or.(x.lt.delta)),0.,y) |
|
87 ;************************************************ |
|
88 ; create default plot |
|
89 ;************************************************ |
|
90 wks = gsn_open_wks("ps","xy") ; open a ps file |
|
91 gsn_define_colormap(wks,"gui_default") ; choose colormap |
|
92 |
|
93 res = True ; Use plot options |
|
94 res@cnFillOn = True ; Turn on color fill |
|
95 res@gsnSpreadColors = True ; use full colormap |
|
96 ; res@cnFillMode = "RasterFill" ; Turn on raster color |
|
97 ; res@lbLabelAutoStride = True |
|
98 res@cnLinesOn = False ; Turn off contourn lines |
|
99 res@mpFillOn = False ; Turn off map fill |
|
100 |
|
101 res@gsnSpreadColors = True ; use full colormap |
|
102 res@cnLevelSelectionMode = "ManualLevels" ; Manual contour invtervals |
|
103 res@cnMinLevelValF = 0. ; Min level |
|
104 res@cnMaxLevelValF = 365. ; Max level |
|
105 res@cnLevelSpacingF = 30. ; interval |
|
106 |
|
107 pres = True ; panel plot mods desired |
|
108 pres@gsnMaximize = True ; fill the page |
|
109 |
|
110 plot=new(3,graphic) ; create graphic array |
|
111 |
|
112 res@tiMainString = "MODIS MOD 15A2 2000-2005" |
|
113 plot(0) = gsn_csm_contour_map_ce(wks,y,res) |
|
114 |
|
115 ; res@tiMainString = "Model i01.03cn" |
|
116 res@tiMainString = "Model i01.04casa" |
|
117 plot(1) = gsn_csm_contour_map_ce(wks,x,res) |
|
118 |
|
119 ; res@tiMainString = "(Model i01.03cn) - (Observed)" |
|
120 res@tiMainString = "(Model i01.04casa) - (Observed)" |
|
121 res@cnMinLevelValF = -120. ; Min level |
|
122 res@cnMaxLevelValF = 120. ; Max level |
|
123 res@cnLevelSpacingF = 20. ; interval |
|
124 plot(2) = gsn_csm_contour_map_ce(wks,d,res) ; for observed |
|
125 |
|
126 gsn_panel(wks,plot,(/3,1/),pres) ; create panel plot |
|
127 system("convert xy.ps xy.png") |
|
128 end |