clamp: all/02.lai.ncl@0c6405ab2ff4

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

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

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

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

     5 ;**************************************************************

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

     7 begin

     8 ; add line to ascci/html file

    10   nnewlines = dimsizes(newlines)

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

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

    13     return

    14   end if

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

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

    17   nline = nline + nnewlines

    18   return

    19 end

    20 ;**************************************************************

    21 ; Main code.

    22 begin

    24   plot_type     = "ps"

    25   plot_type_new = "png"

    27 ;-----------------------------------------------------

    28 ; edit table.html of current model for movel1_vs_model2

    30   if (isvar("compare")) then

    31      html_name2 = compare+"/table.html"

    32      html_new2  = html_name2 +".new"

    33   end if

    35 ;------------------------------------------------------

    36 ; edit table.html for current model

    38   html_name = model_name+"/table.html"

    39   html_new  = html_name +".new"

    41 ;------------------------------------------------------

    42 ; read data: model

    44   fm = addfile(dirm+film10,"r")

    45   laimod = fm->TLAI

    47   delete (fm)

    49   dsizes = dimsizes(laimod)

    50   ntime  = dsizes(0)

    51   nlat   = dsizes(1)

    52   nlon   = dsizes(2)

    54 ;-----------------------------------

    55 ; get landfrac data

    57  film_l   = "lnd_"+ model_grid + ".nc"

    58  fm_l     = addfile (dirs+film_l,"r")

    59  landfrac = fm_l->landfrac

    61  delete (fm_l)

    62 ;----------------------------------

    63 ; read biome data: model

    65   biome_name_mod = "Model PFT Class"

    67   film_c   = "class_pft_"+model_grid+".nc"

    68   fm_c     = addfile(dirs+film_c,"r")

    69   classmod = fm_c->CLASS_PFT

    71   delete (fm_c)

    73 ; model data has 17 land-type classes

    74   nclass_mod = 17

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

    77 ; read data: ob

    79 ;----------------------------------

    80 ; read biome data: observed

    82   biome_name_ob = "MODIS LandCover"

    84   dir_c = diro + "lai/"

    85   filo_c = "land_class_"+model_grid+".nc"

    86   fo = addfile(dir_c+filo_c,"r")

    87   classob = tofloat(fo->LAND_CLASS)

    89   delete (fo)

    91 ; input observed data has 20 land-type classes

    92   nclass_ob = 20

    94 ;---------------------------------

    95 ; read lai data: observed

    97   ;ob_name = "MODIS MOD 15A2 2000-2005"

    98   ob_name = "MODIS MOD 15A2 2000-2004"

   100   dir_l = diro + "lai/"

   101   ;filo = "LAI_2000-2005_MONS_"+model_grid+".nc"

   102   filo = "LAI_2000-2004_MONS_"+model_grid+".nc"

   103   fo = addfile(dir_l+filo,"r")

   104   laiob = fo->LAI

   106   delete (fo)

   108 ;-------------------------------------------------

   109 ; take into account landfrac

   111   laimod = laimod * conform(laimod,landfrac,(/1,2/))

   112   laiob  = laiob  * conform(laiob ,landfrac,(/1,2/))

   114   delete (landfrac)

   116 ;************************************************

   117 ; global res

   118 ;************************************************

   119   resg                      = True             ; Use plot options

   120   resg@cnFillOn             = True             ; Turn on color fill

   121   resg@gsnSpreadColors      = True             ; use full colormap

   122   resg@cnLinesOn            = False            ; Turn off contourn lines

   123   resg@mpFillOn             = False            ; Turn off map fill

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

   126 ;************************************************

   127 ; plot global biome class: (1) observed

   128 ;************************************************

   130   resg@cnMinLevelValF       = 1.              ; Min level

   131   resg@cnMaxLevelValF       = 19.             ; Max level

   132   resg@cnLevelSpacingF      = 1.              ; interval

   134   classob@_FillValue = 1.e+36

   135   classob = where(classob.eq.0,classob@_FillValue,classob)

   137   plot_name = "global_class_ob"

   138   title     = biome_name_ob

   139   resg@tiMainString  = title

   141   wks = gsn_open_wks (plot_type,plot_name)   ; open workstation

   142   gsn_define_colormap(wks,"gui_default")     ; choose colormap

   144   plot = gsn_csm_contour_map_ce(wks,classob,resg)

   146   delete (wks)

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

   149          "rm "+plot_name+"."+plot_type)

   151 ;************************************************

   152 ; plot global biome class: (2) model

   153 ;************************************************

   155   resg@cnMinLevelValF       = 0.              ; Min level

   156   resg@cnMaxLevelValF       = 16.             ; Max level

   157   resg@cnLevelSpacingF      = 1.              ; interval

   159   plot_name = "global_class_model"

   160   title     = biome_name_mod

   161   resg@tiMainString  = title

   163   wks = gsn_open_wks (plot_type,plot_name)   ; open workstation

   164   gsn_define_colormap(wks,"gui_default")     ; choose colormap

   166   plot = gsn_csm_contour_map_ce(wks,classmod,resg)

   168   delete (wks)

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

   171          "rm "+plot_name+"."+plot_type)

   173 ;*******************************************************************

   174 ; for html table : all 3 components (Mean, Max, Phase)

   175 ;*******************************************************************

   177 ; column (not including header column)

   179   component = (/"Mean","Max","Phase"/)

   181   col_head  = (/"observed",model_name,"M_score" \

   182                ,"observed",model_name,"M_score" \

   183                ,"observed",model_name,"M_score" \

   184                /)

   186   n_comp = dimsizes(component)

   187   ncol   = dimsizes(col_head )

   189 ; row (not including header row)

   191 ;----------------------------------------------------

   192 ; using model biome class:

   193   row_head  = (/"Not Vegetated" \

   194                ,"Needleleaf Evergreen Temperate Tree" \

   195                ,"Needleleaf Evergreen Boreal Tree" \

   196 ;              ,"Needleleaf Deciduous Boreal Tree" \

   197                ,"Broadleaf Evergreen Tropical Tree" \

   198                ,"Broadleaf Evergreen Temperate Tree" \

   199                ,"Broadleaf Deciduous Tropical Tree" \

   200                ,"Broadleaf Deciduous Temperate Tree" \

   201 ;              ,"Broadleaf Deciduous Boreal Tree" \

   202 ;              ,"Broadleaf Evergreen Shrub" \

   203                ,"Broadleaf Deciduous Temperate Shrub" \

   204                ,"Broadleaf Deciduous Boreal Shrub" \

   205                ,"C3 Arctic Grass" \

   206                ,"C3 Non-Arctic Grass" \

   207                ,"C4 Grass" \

   208                ,"Corn" \

   209 ;              ,"Wheat" \

   210                ,"All Biomes" \

   211                /)

   213   nrow = dimsizes(row_head)

   215 ; arrays to be passed to table.

   216   text = new ((/nrow, ncol/),string )

   218 ; total M_score

   219   M_total = 0.

   221 ;********************************************************************

   222 ; use land-type class to bin the data in equally spaced ranges

   223 ;********************************************************************

   225 ; using model biome class

   226   nclass = nclass_mod

   228   range  = fspan(0,nclass,nclass+1)

   230 ; Use this range information to grab all the values in a

   231 ; particular range, and then take an average.

   233   nx = dimsizes(range) - 1

   235 ; for 2 data: model and observed

   236   data_n = 2

   238 ; using model biome class

   240   base = ndtooned(classmod)

   242 ; output

   244   yvalues = new((/data_n,nx/),float)

   245   count   = new((/data_n,nx/),float)

   247 ;************************************************************************

   248 ; go through all components

   249 ;************************************************************************

   251  do m = 0,n_comp-1

   253 ;===================

   254 ; get data:

   255 ;===================

   256 ; (A) Mean

   258   if (m .eq. 0) then

   259      data_ob  = dim_avg_Wrap(laiob (lat|:,lon|:,time|:))

   260      data_mod = dim_avg_Wrap(laimod(lat|:,lon|:,time|:))

   261   end if

   263 ; (B) Max

   265   if (m .eq. 1) then

   267 ;    observed

   268      data_ob = laiob(0,:,:)

   269      data_ob@long_name = "Leaf Area Index Max"

   271      do j = 0,nlat-1

   272      do i = 0,nlon-1

   273         data_ob(j,i) = max(laiob(:,j,i))

   274      end do

   275      end do

   277 ;    model

   278      data_mod = laimod(0,:,:)

   279      data_mod@long_name = "Leaf Area Index Max"

   281      do j = 0,nlat-1

   282      do i = 0,nlon-1

   283         data_mod(j,i) = max(laimod(:,j,i))

   284      end do

   285      end do

   287   end if

   289 ; (C) phase

   291   if (m .eq. 2) then

   293 ;    observed

   294      data_ob = laiob(0,:,:)

   295      data_ob@long_name = "Leaf Area Index Max Month"

   297      do j = 0,nlat-1

   298      do i = 0,nlon-1

   299         data_ob(j,i) = maxind(laiob(:,j,i)) + 1

   300      end do

   301      end do

   303 ;    model

   304      data_mod = laimod(0,:,:)

   305      data_mod@long_name = "Leaf Area Index Max Month"

   307      do j = 0,nlat-1

   308      do i = 0,nlon-1

   309         data_mod(j,i) = maxind(laimod(:,j,i)) + 1

   310      end do

   311      end do

   313   end if

   315 ;==============================

   316 ; put data into bins

   317 ;==============================

   319 ; Loop through each range, using base

   321   do i=0,nx-1

   323      if (i.ne.(nx-1)) then

   324         idx = ind((base.ge.range(i)).and.(base.lt.range(i+1)))

   325      else

   326         idx = ind(base.ge.range(i))

   327      end if

   329 ;    loop through each dataset

   331      do n = 0,data_n-1

   333         if (n .eq. 0) then

   334            data = ndtooned(data_ob)

   335         end if

   337         if (n .eq. 1) then

   338            data = ndtooned(data_mod)

   339         end if

   341 ;       Calculate average

   343         if (.not.any(ismissing(idx))) then

   344            yvalues(n,i) = avg(data(idx))

   345            count(n,i)   = dimsizes(idx)

   346         else

   347            yvalues(n,i) = yvalues@_FillValue

   348            count(n,i)   = 0

   349         end if

   351 ;#############################################################

   352 ; using model biome class:

   353 ;

   354 ;     set the following 4 classes to _FillValue:

   355 ;     (3)Needleleaf Deciduous Boreal Tree,

   356 ;     (8)Broadleaf Deciduous Boreal Tree,

   357 ;     (9)Broadleaf Evergreen Shrub,

   358 ;     (16)Wheat

   360       if (i.eq.3 .or. i.eq.8 .or. i.eq.9 .or. i.eq.16) then

   361          yvalues(n,i) = yvalues@_FillValue

   362          count(n,i)   = 0

   363       end if

   364 ;#############################################################

   366       delete(data)

   367     end do                 ; n-loop

   369     delete(idx)

   370   end do                   ; i-loop

   372 ;=====================================

   373 ; compute correlation coef and M score

   374 ;=====================================

   376   score_max = 5.0

   378   u       = yvalues(0,:)

   379   v       = yvalues(1,:)

   380   u_count = count(0,:)

   381   v_count = count(1,:)

   383   good = ind(.not.ismissing(u) .and. .not.ismissing(v))

   385   uu = u(good)

   386   vv = v(good)

   387   uu_count = u_count(good)

   388   vv_count = v_count(good)

   390 ; compute correlation coef

   391   cc = esccr(uu,vv,0)

   393   if (n .eq. 2) then

   394      bias   = avg(abs(vv-uu))

   395      bias   = where((bias.gt. 6.),12.-bias,bias)

   396      Mscore = ((6. - bias)/6.)*score_max

   397   else

   398      bias  = sum(abs(vv-uu)/abs(vv+uu))

   399      Mscore = (1.- (bias/dimsizes(uu)))*score_max

   400   end if

   402   M_score = sprintf("%.2f", Mscore)

   404 ; compute M_total

   406   M_total = M_total + Mscore

   408 ;================================

   409 ; output M_score to score sheet

   410 ;===============================

   412   if (isvar("compare")) then

   413      system("sed -e '1,/M_lai_"+component(m)+"/s/M_lai_"+component(m)+"/"+M_score+"/' "+html_name2+" > "+html_new2+";"+ \

   414             "mv -f "+html_new2+" "+html_name2)

   415   end if

   417   system("sed s#M_lai_"+component(m)+"#"+M_score+"# "+html_name+" > "+html_new+";"+ \

   418          "mv -f "+html_new+" "+html_name)

   420 ;==============================

   421 ; output M_score to html table

   422 ;==============================

   424   n = m*3

   426   do i=0,nrow-2

   427      text(i,n)   = sprintf("%.1f",uu(i))

   428      text(i,n+1) = sprintf("%.1f",vv(i))

   429      text(i,n+2) = "-"

   430   end do

   431   text(nrow-1,n)   = sprintf("%.1f",sum(uu*uu_count)/sum(uu_count))

   432   text(nrow-1,n+1) = sprintf("%.1f",sum(vv*vv_count)/sum(vv_count))

   433   text(nrow-1,n+2) = M_score

   435   delete (u)

   436   delete (v)

   437   delete (uu)

   438   delete (vv)

   439   delete (u_count)

   440   delete (v_count)

   441   delete (uu_count)

   442   delete (vv_count)

   443   delete (good)

   445 ;========================================

   446 ; global res changes for each component

   447 ;========================================

   448   delta = 0.00001

   450   if (m .eq. 0) then

   451      resg@cnMinLevelValF       = 0.

   452      resg@cnMaxLevelValF       = 10.

   453      resg@cnLevelSpacingF      = 1.

   455      data_ob = where(ismissing(data_ob).and.(ismissing(data_mod).or.(data_mod.lt.delta)),0.,data_ob)

   456   end if

   458   if (m .eq. 1) then

   459      resg@cnMinLevelValF       = 0.

   460      resg@cnMaxLevelValF       = 10.

   461      resg@cnLevelSpacingF      = 1.

   463      data_ob = where(ismissing(data_ob).and.(ismissing(data_mod).or.(data_mod.lt.delta)),0.,data_ob)

   464   end if

   466   if (m .eq. 2) then

   467      resg@cnMinLevelValF       = 1.

   468      resg@cnMaxLevelValF       = 12.

   469      resg@cnLevelSpacingF      = 1.

   471      data_mod = where(ismissing(data_ob).or.(data_mod.lt.delta),data_mod@_FillValue,data_mod)

   472      data_ob = where(ismissing(data_ob).and.(ismissing(data_mod).or.(data_mod.lt.delta)),0.,data_ob)

   474   end if

   476 ;=========================

   477 ; global contour : ob

   478 ;=========================

   480   plot_name = "global_"+component(m)+"_ob"

   481   title     = ob_name

   482   resg@tiMainString  = title

   484   wks = gsn_open_wks (plot_type,plot_name)   ; open workstation

   485   gsn_define_colormap(wks,"gui_default")     ; choose colormap

   487   plot = gsn_csm_contour_map_ce(wks,data_ob,resg)

   489   delete (wks)

   490   delete (plot)

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

   493          "rm "+plot_name+"."+plot_type)

   495 ;============================

   496 ; global contour : model

   497 ;============================

   499   plot_name = "global_"+component(m)+"_model"

   500   title     = "Model " + model_name

   501   resg@tiMainString  = title

   503   wks = gsn_open_wks (plot_type,plot_name)

   504   gsn_define_colormap(wks,"gui_default")

   506   plot = gsn_csm_contour_map_ce(wks,data_mod,resg)

   508   delete (wks)

   509   delete (plot)

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

   512          "rm "+plot_name+"."+plot_type)

   514 ;================================

   515 ; global contour: model vs ob

   516 ;================================

   518   plot_name = "global_"+component(m)+"_model_vs_ob"

   520   wks = gsn_open_wks (plot_type,plot_name)

   521   gsn_define_colormap(wks,"gui_default")

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

   525   resg@gsnFrame             = False          ; Do not draw plot

   526   resg@gsnDraw              = False          ; Do not advance frame

   528 ; plot correlation coef

   530   gRes               = True

   531   gRes@txFontHeightF = 0.02

   532   gRes@txAngleF      = 90

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

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

   538 ; plot ob

   540   title     = ob_name

   541   resg@tiMainString  = title

   543   plot(0) = gsn_csm_contour_map_ce(wks,data_ob,resg)

   545 ; plot model

   547   title     = "Model "+ model_name

   548   resg@tiMainString  = title

   550   plot(1) = gsn_csm_contour_map_ce(wks,data_mod,resg)

   552 ; plot model-ob

   554   if (m .eq. 0) then

   555      resg@cnMinLevelValF  = -2.

   556      resg@cnMaxLevelValF  =  2.

   557      resg@cnLevelSpacingF =  0.4

   558   end if

   560   if (m .eq. 1) then

   561      resg@cnMinLevelValF  = -6.

   562      resg@cnMaxLevelValF  =  6.

   563      resg@cnLevelSpacingF =  1.

   564   end if

   566   if (m .eq. 2) then

   567      resg@cnMinLevelValF  = -6.

   568      resg@cnMaxLevelValF  =  6.

   569      resg@cnLevelSpacingF =  1.

   570   end if

   572   zz = data_mod

   573   zz = data_mod - data_ob

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

   575   resg@tiMainString    = title

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

   579 ; plot panel

   581   pres                            = True        ; panel plot mods desired

   582   pres@gsnMaximize                = True        ; fill the page

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

   586   delete (wks)

   587   delete (plot)

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

   590          "rm "+plot_name+"."+plot_type)

   592   delete (data_ob)

   593   delete (data_mod)

   595   resg@gsnFrame             = True          ; Do advance frame

   596   resg@gsnDraw              = True          ; Do draw plot

   598  end do    ; m-loop

   600 ;**************************************************

   601 ; html table

   602 ;**************************************************

   603   output_html = "table_model_vs_ob.html"

   605   header_text = "<H1>LAI: Model "+model_name+" vs. MODIS observations</H1>"

   607   header = (/"<HTML>" \

   608             ,"<HEAD>" \

   609             ,"<TITLE>CLAMP metrics</TITLE>" \

   610             ,"</HEAD>" \

   611             ,header_text \

   612             /)

   613   footer = "</HTML>"

   615   table_header = (/ \

   616         "<table border=1 cellspacing=0 cellpadding=3 width=100%>" \

   617        ,"<tr>" \

   618        ,"   <th bgcolor=DDDDDD rowspan=2>Biome Class</th>" \

   619        ,"   <th bgcolor=DDDDDD colspan=3>"+component(0)+" (m2/m2)</th>" \

   620        ,"   <th bgcolor=DDDDDD colspan=3>"+component(1)+" (m2/m2)</th>" \

   621        ,"   <th bgcolor=DDDDDD colspan=3>"+component(2)+" (months)</th>" \

   622        ,"   <th bgcolor=DDDDDD rowspan=2>Annual Plot</th>" \

   623        ,"</tr>" \

   624        ,"<tr>" \

   625        ,"   <th bgcolor=DDDDDD >observed</th>" \

   626        ,"   <th bgcolor=DDDDDD >"+model_name+"</th>" \

   627        ,"   <th bgcolor=DDDDDD >M_score</th>" \

   628        ,"   <th bgcolor=DDDDDD >observed</th>" \

   629        ,"   <th bgcolor=DDDDDD >"+model_name+"</th>" \

   630        ,"   <th bgcolor=DDDDDD >M_score</th>" \

   631        ,"   <th bgcolor=DDDDDD >observed</th>" \

   632        ,"   <th bgcolor=DDDDDD >"+model_name+"</th>" \

   633        ,"   <th bgcolor=DDDDDD >M_score</th>" \

   634        ,"</tr>" \

   635        /)

   636   table_footer = "</table>"

   637   row_header = "<tr>"

   638   row_footer = "</tr>"

   640   lines = new(50000,string)

   641   nline = 0

   643   set_line(lines,nline,header)

   644   set_line(lines,nline,table_header)

   645 ;-----------------------------------------------

   646 ;row of table

   648   do n = 0,nrow-1

   649      set_line(lines,nline,row_header)

   651      txt1  = row_head(n)

   652      txt2  = text(n,0)

   653      txt3  = text(n,1)

   654      txt4  = text(n,2)

   655      txt5  = text(n,3)

   656      txt6  = text(n,4)

   657      txt7  = text(n,5)

   658      txt8  = text(n,6)

   659      txt9  = text(n,7)

   660      txt10 = text(n,8)

   661      txt11 = "<a href=./annual_biome_"+n+".png>model_vs_ob</a>"

   663      set_line(lines,nline,"<th>"+txt1+"</th>")

   664      set_line(lines,nline,"<th>"+txt2+"</th>")

   665      set_line(lines,nline,"<th>"+txt3+"</th>")

   666      set_line(lines,nline,"<th>"+txt4+"</th>")

   667      set_line(lines,nline,"<th>"+txt5+"</th>")

   668      set_line(lines,nline,"<th>"+txt6+"</th>")

   669      set_line(lines,nline,"<th>"+txt7+"</th>")

   670      set_line(lines,nline,"<th>"+txt8+"</th>")

   671      set_line(lines,nline,"<th>"+txt9+"</th>")

   672      set_line(lines,nline,"<th>"+txt10+"</th>")

   673      set_line(lines,nline,"<th>"+txt11+"</th>")

   675      set_line(lines,nline,row_footer)

   676   end do

   677 ;-----------------------------------------------

   678   set_line(lines,nline,table_footer)

   679   set_line(lines,nline,footer)

   681 ; Now write to an HTML file

   683   idx = ind(.not.ismissing(lines))

   684   if(.not.any(ismissing(idx))) then

   685     asciiwrite(output_html,lines(idx))

   686   else

   687    print ("error?")

   688   end if

   690   delete (yvalues)

   691   delete (count)

   692   delete (idx)

   694 ;************************************************************************

   695 ; go through all ntime

   696 ;************************************************************************

   698 ; for 2 data: model and observed

   699   data_n = 2

   701 ; using model biome class

   703   base = ndtooned(classmod)

   705 ; output

   707   yvalues = new((/data_n,nx,ntime/),float)

   709 ;==============================

   710 ; put data into bins

   711 ;==============================

   713 ; Loop through each range, using base

   715   do i=0,nx-1

   717      if (i.ne.(nx-1)) then

   718         idx = ind((base.ge.range(i)).and.(base.lt.range(i+1)))

   719      else

   720         idx = ind(base.ge.range(i))

   721      end if

   723 ;    loop through each dataset

   725      do n = 0,data_n-1

   727         do m = 0,ntime-1

   729         if (n .eq. 0) then

   730            data = ndtooned(laiob (m,:,:))

   731         end if

   733         if (n .eq. 1) then

   734            data = ndtooned(laimod(m,:,:))

   735         end if

   737 ;       Calculate average

   739         if (.not.any(ismissing(idx))) then

   740            yvalues(n,i,m) = avg(data(idx))

   741         else

   742            yvalues(n,i,m) = yvalues@_FillValue

   743         end if

   745 ;#############################################################

   746 ; using model biome class:

   747 ;

   748 ;     set the following 4 classes to _FillValue:

   749 ;     (3)Needleleaf Deciduous Boreal Tree,

   750 ;     (8)Broadleaf Deciduous Boreal Tree,

   751 ;     (9)Broadleaf Evergreen Shrub,

   752 ;     (16)Wheat

   754       if (i.eq.3 .or. i.eq.8 .or. i.eq.9 .or. i.eq.16) then

   755          yvalues(n,i,m) = yvalues@_FillValue

   756       end if

   757 ;#############################################################

   759       end do               ; m-loop

   761       delete(data)

   762     end do                 ; n-loop

   764     delete(idx)

   765   end do                   ; i-loop

   767   good = ind(.not.ismissing(yvalues(0,:,0)))

   769   n_biome = dimsizes(good)

   771 ;----------------------------------------------------------------

   772 ; data for tseries plot

   774   yvalues_g = new((/data_n,n_biome,ntime/),float)

   776   yvalues_g@units = ""

   778   yvalues_g = yvalues(:,good,:)

   780   delete (good)

   782 ;*******************************************************************

   783 ; res for line plot

   784 ;*******************************************************************

   785 ; for x-axis in xyplot

   786   mon = ispan(1,12,1)

   787   mon@long_name = "month"

   789   res                   = True                      ; plot mods desired

   790   res@xyLineThicknesses = (/2.0,2.0,2.0/)           ; make 2nd lines thicker

   791   res@xyLineColors      = (/"blue","red"/)         ; change line color

   793 ; res@tiMainFontHeightF = 0.025                     ; size of title

   795 ; Add a boxed legend using the more simple method

   796   res@pmLegendDisplayMode    = "Always"

   797 ; res@pmLegendWidthF         = 0.1

   798   res@pmLegendWidthF         = 0.08

   799   res@pmLegendHeightF        = 0.06

   800 ; res@pmLegendOrthogonalPosF = -1.17

   801 ; res@pmLegendOrthogonalPosF = -1.00  ;(downward)

   802   res@pmLegendOrthogonalPosF = -0.30  ;(downward)

   804 ; res@pmLegendParallelPosF   =  0.18

   805   res@pmLegendParallelPosF   =  0.23  ;(rightward)

   807 ; res@lgPerimOn             = False

   808   res@lgLabelFontHeightF     = 0.015

   809   res@xyExplicitLegendLabels = (/"observed",model_name/)

   811 ;************************************************************

   813   plot_data = new((/2,ntime/),float)

   814   plot_data@long_name = "Leaf Area Index"

   816 ;----------------------------------------------

   817 ; time series plot : per biome

   819   do m = 0, n_biome-1

   821      plot_name = "annual_biome_"+ m

   823      wks = gsn_open_wks (plot_type,plot_name)

   825      title = "LAI : "+ row_head(m)

   826      res@tiMainString = title

   828      plot_data(0,:) = yvalues_g(0,m,:)

   829      plot_data(1,:) = yvalues_g(1,m,:)

   831      plot = gsn_csm_xy(wks,mon,plot_data,res)

   833      delete (wks)

   834      delete (plot)

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

   837             "rm "+plot_name+"."+plot_type)

   838   end do

   840 ;--------------------------------------------

   841 ; time series plot: all biome

   843      plot_name = "annual_biome_"+ n_biome

   845      wks = gsn_open_wks (plot_type,plot_name)

   847      title = "LAI : "+ row_head(n_biome)

   848      res@tiMainString = title

   850      do k = 0,ntime-1

   851         plot_data(0,k) = avg(yvalues_g(0,:,k))

   852         plot_data(1,k) = avg(yvalues_g(1,:,k))

   853      end do

   855      plot = gsn_csm_xy(wks,mon,plot_data,res)

   857      delete (wks)

   858      delete (plot)

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

   861             "rm "+plot_name+"."+plot_type)

   863   delete (plot_data)

   865 ;***************************************************************************

   866 ; write total score to file

   867 ;***************************************************************************

   869   asciiwrite("M_save.lai", M_total)

   871 ;***************************************************************************

   872 ; output plots

   873 ;***************************************************************************

   874   output_dir = model_name+"/lai"

   876   system("mv *.png *.html " + output_dir)

   877 ;***************************************************************************

   878 end

author	Forrest Hoffman <forrest@climatemodeling.org>
	Mon, 26 Jan 2009 22:08:20 -0500
changeset 0	0c6405ab2ff4
child 1	4be95183fbcd
permissions	-rw-r--r--