1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/fire/21.table+tseries.ncl Mon Jan 26 22:08:20 2009 -0500
1.3 @@ -0,0 +1,258 @@
1.4 +;********************************************************
1.5 +;using model biome vlass
1.6 +;
1.7 +; required command line input parameters:
1.8 +; ncl 'model_name="10cn" model_grid="T42" dirm="/.../ film="..."' 01.npp.ncl
1.9 +;
1.10 +; histogram normalized by rain and compute correleration
1.11 +;**************************************************************
1.12 +load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
1.13 +load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
1.14 +load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
1.15 +;**************************************************************
1.16 +procedure set_line(lines:string,nline:integer,newlines:string)
1.17 +begin
1.18 +; add line to ascci/html file
1.19 +
1.20 + nnewlines = dimsizes(newlines)
1.21 + if(nline+nnewlines-1.ge.dimsizes(lines))
1.22 + print("set_line: bad index, not setting anything.")
1.23 + return
1.24 + end if
1.25 + lines(nline:nline+nnewlines-1) = newlines
1.26 +; print ("lines = " + lines(nline:nline+nnewlines-1))
1.27 + nline = nline + nnewlines
1.28 + return
1.29 +end
1.30 +;**************************************************************
1.31 +; Main code.
1.32 +begin
1.33 +
1.34 + plot_type = "ps"
1.35 + plot_type_new = "png"
1.36 +
1.37 +;---------------------------------------------------
1.38 +; model name and grid
1.39 +
1.40 + model_grid = "T42"
1.41 +
1.42 + model_name = "cn"
1.43 + model_name1 = "i01.06cn"
1.44 + model_name2 = "i01.10cn"
1.45 +
1.46 +;---------------------------------------------------
1.47 +; get biome data: model
1.48 +
1.49 + biome_name_mod = "Model PFT Class"
1.50 +
1.51 + dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/"
1.52 + film = "class_pft_"+model_grid+".nc"
1.53 + fm = addfile(dirm+film,"r")
1.54 +
1.55 + classmod = fm->CLASS_PFT
1.56 +
1.57 + delete (fm)
1.58 +
1.59 +; model data has 17 land-type classes
1.60 +
1.61 + nclass_mod = 17
1.62 +
1.63 +;--------------------------------------------------
1.64 +; get model data: landfrac and area
1.65 +
1.66 + dirm = "/fis/cgd/cseg/people/jeff/surface_data/"
1.67 + film = "lnd_T42.nc"
1.68 + fm = addfile (dirm+film,"r")
1.69 +
1.70 + landmask = fm->landmask
1.71 + landfrac = fm->landfrac
1.72 + area = fm->area
1.73 +
1.74 + delete (fm)
1.75 +
1.76 +; change area from km**2 to m**2
1.77 + area = area * 1.e6
1.78 +
1.79 +;----------------------------------------------------
1.80 +; read data: time series, model
1.81 +
1.82 + dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/"
1.83 + film = model_name2 + "_Fire_C_1979-2004_monthly.nc"
1.84 + fm = addfile (dirm+film,"r")
1.85 +
1.86 + data_mod = fm->COL_FIRE_CLOSS(18:25,:,:,:)
1.87 +
1.88 + delete (fm)
1.89 +
1.90 +; Units for these variables are:
1.91 +; g C/m^2/s
1.92 +
1.93 +; change unit to g C/m^2/month
1.94 +
1.95 + nsec_per_month = 60*60*24*30
1.96 +
1.97 + data_mod = data_mod * nsec_per_month
1.98 +
1.99 + data_mod@unit = "gC/m2/month"
1.100 +;----------------------------------------------------
1.101 +; read data: time series, observed
1.102 +
1.103 + dirm = "/fis/cgd/cseg/people/jeff/fire_data/ob/GFEDv2_C/"
1.104 + film = "Fire_C_1997-2006_monthly_"+ model_grid+".nc"
1.105 + fm = addfile (dirm+film,"r")
1.106 +
1.107 + data_ob = fm->FIRE_C(0:7,:,:,:)
1.108 +
1.109 + delete (fm)
1.110 +
1.111 + ob_name = "GFEDv2"
1.112 +
1.113 +; Units for these variables are:
1.114 +; g C/m^2/month
1.115 +
1.116 +;---------------------------------------------------
1.117 +; take into account landfrac
1.118 +
1.119 + data_mod = data_mod * conform(data_mod, landfrac, (/2,3/))
1.120 + data_ob = data_ob * conform(data_ob, landfrac, (/2,3/))
1.121 +
1.122 +;---------------------------------------------------
1.123 +; get time-mean
1.124 +
1.125 + x = dim_avg_Wrap(data_mod(lat|:,lon|:,month|:,year|:))
1.126 + data_mod_m = dim_avg_Wrap( x(lat|:,lon|:,month|:))
1.127 + delete (x)
1.128 +
1.129 + x = dim_avg_Wrap( data_ob(lat|:,lon|:,month|:,year|:))
1.130 + data_ob_m = dim_avg_Wrap( x(lat|:,lon|:,month|:))
1.131 + delete (x)
1.132 +
1.133 +; printVarSummary(y)
1.134 +
1.135 +;----------------------------------------------------
1.136 +; compute correlation coef
1.137 +
1.138 + landmask_1d = ndtooned(landmask)
1.139 + data_mod_1d = ndtooned(data_mod_m)
1.140 + data_ob_1d = ndtooned(data_ob_m)
1.141 +
1.142 + good = ind(landmask_1d .gt. 0.)
1.143 +; print (dimsizes(good))
1.144 +
1.145 + cc = esccr(data_mod_1d(good),data_ob_1d(good),0)
1.146 +; print (cc)
1.147 +
1.148 + delete (landmask_1d)
1.149 + delete (data_mod_1d)
1.150 + delete (data_ob_id)
1.151 +
1.152 +;----------------------------------------------------
1.153 +; compute M_global
1.154 +
1.155 + score_max = 1.
1.156 +
1.157 + Mscore = cc * cc * score_max
1.158 +
1.159 + M_global = sprintf("%.2f", Mscore)
1.160 +
1.161 +;----------------------------------------------------
1.162 +; global res
1.163 +
1.164 + resg = True ; Use plot options
1.165 + resg@cnFillOn = True ; Turn on color fill
1.166 + resg@gsnSpreadColors = True ; use full colormap
1.167 + resg@cnLinesOn = False ; Turn off contourn lines
1.168 + resg@mpFillOn = False ; Turn off map fill
1.169 + resg@cnLevelSelectionMode = "ManualLevels" ; Manual contour invtervals
1.170 +
1.171 +;----------------------------------------------------
1.172 +; global contour: model vs ob
1.173 +
1.174 + plot_name = "global_model_vs_ob"
1.175 +
1.176 + wks = gsn_open_wks (plot_type,plot_name)
1.177 + gsn_define_colormap(wks,"gui_default")
1.178 +
1.179 + plot=new(3,graphic) ; create graphic array
1.180 +
1.181 + resg@gsnFrame = False ; Do not draw plot
1.182 + resg@gsnDraw = False ; Do not advance frame
1.183 +
1.184 +;----------------------
1.185 +; plot correlation coef
1.186 +
1.187 + gRes = True
1.188 + gRes@txFontHeightF = 0.02
1.189 + gRes@txAngleF = 90
1.190 +
1.191 + correlation_text = "(correlation coef = "+sprintf("%.2f", cc)+")"
1.192 +
1.193 + gsn_text_ndc(wks,correlation_text,0.20,0.50,gRes)
1.194 +
1.195 +;-----------------------
1.196 +; plot ob
1.197 +
1.198 + data_ob_m = where(landmask .gt. 0., data_ob_m, data_ob_m@_FillValue)
1.199 +
1.200 + title = ob_name
1.201 + resg@tiMainString = title
1.202 +
1.203 + resg@cnMinLevelValF = 1.
1.204 + resg@cnMaxLevelValF = 10.
1.205 + resg@cnLevelSpacingF = 1.
1.206 +
1.207 + plot(0) = gsn_csm_contour_map_ce(wks,data_ob_m,resg)
1.208 +
1.209 +;-----------------------
1.210 +; plot model
1.211 +
1.212 + data_mod_m = where(landmask .gt. 0., data_mod_m, data_mod_m@_FillValue)
1.213 +
1.214 + title = "Model "+ model_name
1.215 + resg@tiMainString = title
1.216 +
1.217 + resg@cnMinLevelValF = 1.
1.218 + resg@cnMaxLevelValF = 10.
1.219 + resg@cnLevelSpacingF = 1.
1.220 +
1.221 + plot(1) = gsn_csm_contour_map_ce(wks,data_mod_m,resg)
1.222 +
1.223 +;-----------------------
1.224 +; plot model-ob
1.225 +
1.226 + resg@cnMinLevelValF = -8.
1.227 + resg@cnMaxLevelValF = 2.
1.228 + resg@cnLevelSpacingF = 1.
1.229 +
1.230 + zz = data_ob_m
1.231 + zz = data_mod_m - data_ob_m
1.232 + title = "Model_"+model_name+" - Observed"
1.233 + resg@tiMainString = title
1.234 +
1.235 + plot(2) = gsn_csm_contour_map_ce(wks,zz,resg)
1.236 +
1.237 +; plot panel
1.238 +
1.239 + pres = True ; panel plot mods desired
1.240 + pres@gsnMaximize = True ; fill the page
1.241 +
1.242 + gsn_panel(wks,plot,(/3,1/),pres) ; create panel plot
1.243 +
1.244 +exit
1.245 +
1.246 +; system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
1.247 +; "rm "+plot_name+"."+plot_type)
1.248 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
1.249 +
1.250 + clear (wks)
1.251 + delete (plot)
1.252 +
1.253 + delete (data_ob_m)
1.254 + delete (data_mod_m)
1.255 + delete (zz)
1.256 +
1.257 + resg@gsnFrame = True ; Do advance frame
1.258 + resg@gsnDraw = True ; Do draw plot
1.259 +
1.260 +end
1.261 +