clamp: fire/21.table+tseries.ncl@4be95183fbcd

     1 ;********************************************************

     2 ;using model biome vlass

3 ;

     4 ; required command line input parameters:

     5 ;  ncl 'model_name="10cn" model_grid="T42" dirm="/.../ film="..."' 01.npp.ncl

6 ;

     7 ; histogram normalized by rain and compute correleration

     8 ;**************************************************************

     9 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"

    10 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"

    11 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"

    12 ;**************************************************************

    13 procedure set_line(lines:string,nline:integer,newlines:string)

    14 begin

    15 ; add line to ascci/html file

    17   nnewlines = dimsizes(newlines)

    18   if(nline+nnewlines-1.ge.dimsizes(lines))

    19     print("set_line: bad index, not setting anything.")

    20     return

    21   end if

    22   lines(nline:nline+nnewlines-1) = newlines

    23 ;  print ("lines = " + lines(nline:nline+nnewlines-1))

    24   nline = nline + nnewlines

    25   return

    26 end

    27 ;**************************************************************

    28 ; Main code.

    29 begin

    31   plot_type     = "ps"

    32   plot_type_new = "png"

    34 ;---------------------------------------------------

    35 ; model name and grid

    37   model_grid = "T42"

    39   model_name  = "cn"

    40   model_name1 = "i01.06cn"

    41   model_name2 = "i01.10cn"

    43 ;---------------------------------------------------

    44 ; get biome data: model

    46   biome_name_mod = "Model PFT Class"

    48   dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/"

    49   film = "class_pft_"+model_grid+".nc"

    50   fm   = addfile(dirm+film,"r")

    52   classmod = fm->CLASS_PFT

    54   delete (fm)

    56 ; model data has 17 land-type classes

    58   nclass_mod = 17

    60 ;--------------------------------------------------

    61 ; get model data: landfrac and area

    63   dirm = "/fis/cgd/cseg/people/jeff/surface_data/"

    64   film = "lnd_T42.nc"

    65   fm   = addfile (dirm+film,"r")

    67   landmask = fm->landmask

    68   landfrac = fm->landfrac

    69   area     = fm->area

    71   delete (fm)

    73 ; change area from km**2 to m**2

    74   area = area * 1.e6

    76 ;----------------------------------------------------

    77 ; read data: time series, model

    79  dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/"

    80  film = model_name2 + "_Fire_C_1979-2004_monthly.nc"

    81  fm   = addfile (dirm+film,"r")

    83  data_mod = fm->COL_FIRE_CLOSS(18:25,:,:,:)

    85  delete (fm)

    87 ; Units for these variables are:

    88 ; g C/m^2/s

    90 ; change unit to g C/m^2/month

    92   nsec_per_month = 60*60*24*30

    94   data_mod = data_mod * nsec_per_month

    96   data_mod@unit = "gC/m2/month"

    97 ;----------------------------------------------------

    98 ; read data: time series, observed

   100  dirm = "/fis/cgd/cseg/people/jeff/fire_data/ob/GFEDv2_C/"

   101  film = "Fire_C_1997-2006_monthly_"+ model_grid+".nc"

   102  fm   = addfile (dirm+film,"r")

   104  data_ob = fm->FIRE_C(0:7,:,:,:)

   106  delete (fm)

   108  ob_name = "GFEDv2"

   110 ; Units for these variables are:

   111 ; g C/m^2/month

   113 ;---------------------------------------------------

   114 ; take into account landfrac

   116  data_mod = data_mod * conform(data_mod, landfrac, (/2,3/))

   117  data_ob  = data_ob  * conform(data_ob,  landfrac, (/2,3/))

   119 ;---------------------------------------------------

   120 ; get time-mean

   122   x          = dim_avg_Wrap(data_mod(lat|:,lon|:,month|:,year|:))

   123   data_mod_m = dim_avg_Wrap(       x(lat|:,lon|:,month|:))

   124   delete (x)

   126   x          = dim_avg_Wrap( data_ob(lat|:,lon|:,month|:,year|:))

   127   data_ob_m  = dim_avg_Wrap(       x(lat|:,lon|:,month|:))

   128   delete (x)

   130 ; printVarSummary(y)

   132 ;----------------------------------------------------

   133 ; compute correlation coef

   135   landmask_1d = ndtooned(landmask)

   136   data_mod_1d = ndtooned(data_mod_m)

   137   data_ob_1d  = ndtooned(data_ob_m)

   139   good = ind(landmask_1d .gt. 0.)

   140 ; print (dimsizes(good))

   142   cc = esccr(data_mod_1d(good),data_ob_1d(good),0)

   143 ; print (cc)

   145   delete (landmask_1d)

   146   delete (data_mod_1d)

   147   delete (data_ob_id)

   149 ;----------------------------------------------------

   150 ; compute M_global

   152   score_max = 1.

   154   Mscore = cc * cc * score_max

   156   M_global = sprintf("%.2f", Mscore)

   158 ;----------------------------------------------------

   159 ; global res

   161   resg                      = True             ; Use plot options

   162   resg@cnFillOn             = True             ; Turn on color fill

   163   resg@gsnSpreadColors      = True             ; use full colormap

   164   resg@cnLinesOn            = False            ; Turn off contourn lines

   165   resg@mpFillOn             = False            ; Turn off map fill

   166   resg@cnLevelSelectionMode = "ManualLevels"   ; Manual contour invtervals

   168 ;----------------------------------------------------

   169 ; global contour: model vs ob

   171   plot_name = "global_model_vs_ob"

   173   wks = gsn_open_wks (plot_type,plot_name)

   174   gsn_define_colormap(wks,"gui_default")

   176   plot=new(3,graphic)                        ; create graphic array

   178   resg@gsnFrame             = False          ; Do not draw plot

   179   resg@gsnDraw              = False          ; Do not advance frame

   181 ;----------------------

   182 ; plot correlation coef

   184   gRes               = True

   185   gRes@txFontHeightF = 0.02

   186   gRes@txAngleF      = 90

   188   correlation_text = "(correlation coef = "+sprintf("%.2f", cc)+")"

   190   gsn_text_ndc(wks,correlation_text,0.20,0.50,gRes)

   192 ;-----------------------

   193 ; plot ob

   195   data_ob_m = where(landmask .gt. 0., data_ob_m, data_ob_m@_FillValue)

   197   title     = ob_name

   198   resg@tiMainString  = title

   200   resg@cnMinLevelValF       = 1.

   201   resg@cnMaxLevelValF       = 10.

   202   resg@cnLevelSpacingF      = 1.

   204   plot(0) = gsn_csm_contour_map_ce(wks,data_ob_m,resg)

   206 ;-----------------------

   207 ; plot model

   209   data_mod_m = where(landmask .gt. 0., data_mod_m, data_mod_m@_FillValue)

   211   title     = "Model "+ model_name

   212   resg@tiMainString  = title

   214   resg@cnMinLevelValF       = 1.

   215   resg@cnMaxLevelValF       = 10.

   216   resg@cnLevelSpacingF      = 1.

   218   plot(1) = gsn_csm_contour_map_ce(wks,data_mod_m,resg)

   220 ;-----------------------

   221 ; plot model-ob

   223      resg@cnMinLevelValF  = -8.

   224      resg@cnMaxLevelValF  =  2.

   225      resg@cnLevelSpacingF =  1.

   227   zz = data_ob_m

   228   zz = data_mod_m - data_ob_m

   229   title = "Model_"+model_name+" - Observed"

   230   resg@tiMainString    = title

   232   plot(2) = gsn_csm_contour_map_ce(wks,zz,resg)

   234 ; plot panel

   236   pres                            = True        ; panel plot mods desired

   237   pres@gsnMaximize                = True        ; fill the page

   239   gsn_panel(wks,plot,(/3,1/),pres)              ; create panel plot

   241 exit

   243 ; system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \

   244 ;        "rm "+plot_name+"."+plot_type)

   245   system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)

   247   clear (wks)

   248   delete (plot)

   250   delete (data_ob_m)

   251   delete (data_mod_m)

   252   delete (zz)

   254   resg@gsnFrame             = True          ; Do advance frame

   255   resg@gsnDraw              = True          ; Do draw plot

   257 end

author	Forrest Hoffman <forrest@climatemodeling.org>
	Thu, 26 Mar 2009 14:02:21 -0400
changeset 1	4be95183fbcd
permissions	-rw-r--r--