clamp: all/02.lai.ncl@4be95183fbcd

     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 MOD15A2 (2000-2005)"

    98   ob_name = "MODIS MOD15A2 (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 -trim -density 100 "+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 -trim -density 100 "+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      score_max = 4.0

   262   end if

   264 ; (B) Max

   266   if (m .eq. 1) then

   268 ;    observed

   269      data_ob = laiob(0,:,:)

   270      data_ob@long_name = "Leaf Area Index Max"

   272      do j = 0,nlat-1

   273      do i = 0,nlon-1

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

   275      end do

   276      end do

   278 ;    model

   279      data_mod = laimod(0,:,:)

   280      data_mod@long_name = "Leaf Area Index Max"

   282      do j = 0,nlat-1

   283      do i = 0,nlon-1

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

   285      end do

   286      end do

   288      score_max = 5.0

   290   end if

   292 ; (C) phase

   294   if (m .eq. 2) then

   296 ;    observed

   297      data_ob = laiob(0,:,:)

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

   300      do j = 0,nlat-1

   301      do i = 0,nlon-1

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

   303      end do

   304      end do

   306 ;    model

   307      data_mod = laimod(0,:,:)

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

   310      do j = 0,nlat-1

   311      do i = 0,nlon-1

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

   313      end do

   314      end do

   316      score_max = 6.0

   318   end if

   320 ;==============================

   321 ; put data into bins

   322 ;==============================

   324 ; Loop through each range, using base

   326   do i=0,nx-1

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

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

   330      else

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

   332      end if

   334 ;    loop through each dataset

   336      do n = 0,data_n-1

   338         if (n .eq. 0) then

   339            data = ndtooned(data_ob)

   340         end if

   342         if (n .eq. 1) then

   343            data = ndtooned(data_mod)

   344         end if

   346 ;       Calculate average

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

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

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

   351         else

   352            yvalues(n,i) = yvalues@_FillValue

   353            count(n,i)   = 0

   354         end if

   356 ;#############################################################

   357 ; using model biome class:

   358 ;

   359 ;     set the following 4 classes to _FillValue:

   360 ;     (3)Needleleaf Deciduous Boreal Tree,

   361 ;     (8)Broadleaf Deciduous Boreal Tree,

   362 ;     (9)Broadleaf Evergreen Shrub,

   363 ;     (16)Wheat

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

   366          yvalues(n,i) = yvalues@_FillValue

   367          count(n,i)   = 0

   368       end if

   369 ;#############################################################

   371       delete(data)

   372     end do                 ; n-loop

   374     delete(idx)

   375   end do                   ; i-loop

   377 ;=====================================

   378 ; compute correlation coef and M score

   379 ;=====================================

   381   ; FMH: Max scores are depend on which component is being analyzed, so this

   382   ; is now set in each if statement above.

   383   ;score_max = 5.0

   385   u       = yvalues(0,:)

   386   v       = yvalues(1,:)

   387   u_count = count(0,:)

   388   v_count = count(1,:)

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

   392   uu = u(good)

   393   vv = v(good)

   394   uu_count = u_count(good)

   395   vv_count = v_count(good)

   397 ; compute correlation coef

   398   cc = esccr(uu,vv,0)

   400   if (n .eq. 2) then

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

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

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

   404   else

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

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

   407   end if

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

   411 ; compute M_total

   413   M_total = M_total + Mscore

   415 ;================================

   416 ; output M_score to score sheet

   417 ;===============================

   419   if (isvar("compare")) then

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

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

   422   end if

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

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

   427 ;==============================

   428 ; output M_score to html table

   429 ;==============================

   431   n = m*3

   433   do i=0,nrow-2

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

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

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

   437   end do

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

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

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

   442   delete (u)

   443   delete (v)

   444   delete (uu)

   445   delete (vv)

   446   delete (u_count)

   447   delete (v_count)

   448   delete (uu_count)

   449   delete (vv_count)

   450   delete (good)

   452 ;========================================

   453 ; global res changes for each component

   454 ;========================================

   455   delta = 0.00001

   457   if (m .eq. 0) then

   458      resg@cnMinLevelValF       = 0.

   459      resg@cnMaxLevelValF       = 10.

   460      resg@cnLevelSpacingF      = 1.

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

   463   end if

   465   if (m .eq. 1) then

   466      resg@cnMinLevelValF       = 0.

   467      resg@cnMaxLevelValF       = 10.

   468      resg@cnLevelSpacingF      = 1.

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

   471   end if

   473   if (m .eq. 2) then

   474      resg@cnMinLevelValF       = 1.

   475      resg@cnMaxLevelValF       = 12.

   476      resg@cnLevelSpacingF      = 1.

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

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

   481   end if

   483 ;=========================

   484 ; global contour : ob

   485 ;=========================

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

   488   title     = ob_name

   489   resg@tiMainString  = title

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

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

   494   plot = gsn_csm_contour_map_ce(wks,data_ob,resg)

   496   delete (wks)

   497   delete (plot)

   499   system("convert -trim -density 100 "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \

   500          "rm "+plot_name+"."+plot_type)

   501   ;system("convert -trim -density 100 "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)

   504 ;============================

   505 ; global contour : model

   506 ;============================

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

   509   title     = "Model " + model_name + " (2000-2004)"

   510   resg@tiMainString  = title

   512   wks = gsn_open_wks (plot_type,plot_name)

   513   gsn_define_colormap(wks,"gui_default")

   515   plot = gsn_csm_contour_map_ce(wks,data_mod,resg)

   517   delete (wks)

   518   delete (plot)

   520   system("convert -trim -density 100 "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \

   521          "rm "+plot_name+"."+plot_type)

   522   ;system("convert -trim -density 100 "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)

   524 ;================================

   525 ; global contour: model vs ob

   526 ;================================

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

   530   wks = gsn_open_wks (plot_type,plot_name)

   531   gsn_define_colormap(wks,"gui_default")

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

   535   resg@gsnFrame             = False          ; Do not draw plot

   536   resg@gsnDraw              = False          ; Do not advance frame

   538 ; plot correlation coef

   540   gRes               = True

   541   gRes@txFontHeightF = 0.02

   542   gRes@txAngleF      = 90

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

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

   548 ; plot ob

   550   title     = ob_name

   551   resg@tiMainString  = title

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

   555 ; plot model

   557   title     = "Model "+ model_name + " (2000-2004)"

   558   resg@tiMainString  = title

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

   562 ; plot model-ob

   564   if (m .eq. 0) then

   565      resg@cnMinLevelValF  = -2.

   566      resg@cnMaxLevelValF  =  2.

   567      resg@cnLevelSpacingF =  0.4

   568   end if

   570   if (m .eq. 1) then

   571      resg@cnMinLevelValF  = -6.

   572      resg@cnMaxLevelValF  =  6.

   573      resg@cnLevelSpacingF =  1.

   574   end if

   576   if (m .eq. 2) then

   577      resg@cnMinLevelValF  = -6.

   578      resg@cnMaxLevelValF  =  6.

   579      resg@cnLevelSpacingF =  1.

   580   end if

   582   zz = data_mod

   583   zz = data_mod - data_ob

   584   title = "Model_"+model_name+" - MODIS MOD15A2"

   585   resg@tiMainString    = title

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

   589 ; plot panel

   591   pres                            = True        ; panel plot mods desired

   592   pres@gsnMaximize                = True        ; fill the page

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

   596   delete (wks)

   597   delete (plot)

   599   system("convert -trim -density 150 "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \

   600          "rm "+plot_name+"."+plot_type)

   601   ;system("convert -trim -density 150 "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)

   603   delete (data_ob)

   604   delete (data_mod)

   606   resg@gsnFrame             = True          ; Do advance frame

   607   resg@gsnDraw              = True          ; Do draw plot

   609  end do    ; m-loop

   611 ;**************************************************

   612 ; html table

   613 ;**************************************************

   614   output_html = "table_model_vs_ob.html"

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

   618   header = (/"<HTML>" \

   619             ,"<HEAD>" \

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

   621             ,"</HEAD>" \

   622             ,header_text \

   623             /)

   624   footer = "</HTML>"

   626   table_header = (/ \

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

   628        ,"<tr>" \

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

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

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

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

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

   634        ,"</tr>" \

   635        ,"<tr>" \

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

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

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

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

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

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

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

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

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

   645        ,"</tr>" \

   646        /)

   647   table_footer = "</table>"

   648   row_header = "<tr>"

   649   row_footer = "</tr>"

   651   lines = new(50000,string)

   652   nline = 0

   654   set_line(lines,nline,header)

   655   set_line(lines,nline,table_header)

   656 ;-----------------------------------------------

   657 ;row of table

   659   do n = 0,nrow-1

   660      set_line(lines,nline,row_header)

   662      txt1  = row_head(n)

   663      txt2  = text(n,0)

   664      txt3  = text(n,1)

   665      txt4  = text(n,2)

   666      txt5  = text(n,3)

   667      txt6  = text(n,4)

   668      txt7  = text(n,5)

   669      txt8  = text(n,6)

   670      txt9  = text(n,7)

   671      txt10 = text(n,8)

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

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

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

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

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

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

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

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

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

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

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

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

   686      set_line(lines,nline,row_footer)

   687   end do

   688 ;-----------------------------------------------

   689   set_line(lines,nline,table_footer)

   690   set_line(lines,nline,footer)

   692 ; Now write to an HTML file

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

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

   696     asciiwrite(output_html,lines(idx))

   697   else

   698    print ("error?")

   699   end if

   701   delete (yvalues)

   702   delete (count)

   703   delete (idx)

   705 ;************************************************************************

   706 ; go through all ntime

   707 ;************************************************************************

   709 ; for 2 data: model and observed

   710   data_n = 2

   712 ; using model biome class

   714   base = ndtooned(classmod)

   716 ; output

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

   720 ;==============================

   721 ; put data into bins

   722 ;==============================

   724 ; Loop through each range, using base

   726   do i=0,nx-1

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

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

   730      else

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

   732      end if

   734 ;    loop through each dataset

   736      do n = 0,data_n-1

   738         do m = 0,ntime-1

   740         if (n .eq. 0) then

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

   742         end if

   744         if (n .eq. 1) then

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

   746         end if

   748 ;       Calculate average

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

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

   752         else

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

   754         end if

   756 ;#############################################################

   757 ; using model biome class:

   758 ;

   759 ;     set the following 4 classes to _FillValue:

   760 ;     (3)Needleleaf Deciduous Boreal Tree,

   761 ;     (8)Broadleaf Deciduous Boreal Tree,

   762 ;     (9)Broadleaf Evergreen Shrub,

   763 ;     (16)Wheat

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

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

   767       end if

   768 ;#############################################################

   770       end do               ; m-loop

   772       delete(data)

   773     end do                 ; n-loop

   775     delete(idx)

   776   end do                   ; i-loop

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

   780   n_biome = dimsizes(good)

   782 ;----------------------------------------------------------------

   783 ; data for tseries plot

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

   787   yvalues_g@units = ""

   789   yvalues_g = yvalues(:,good,:)

   791   delete (good)

   793 ;*******************************************************************

   794 ; res for line plot

   795 ;*******************************************************************

   796 ; for x-axis in xyplot

   797   mon = ispan(1,12,1)

   798   mon@long_name = "month"

   800   res                   = True                      ; plot mods desired

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

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

   804 ; res@tiMainFontHeightF = 0.025                     ; size of title

   806 ; Add a boxed legend using the more simple method

   807   res@pmLegendDisplayMode    = "Always"

   808 ; res@pmLegendWidthF         = 0.1

   809   res@pmLegendWidthF         = 0.08

   810   res@pmLegendHeightF        = 0.06

   811 ; res@pmLegendOrthogonalPosF = -1.17

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

   813   res@pmLegendOrthogonalPosF = -0.30  ;(downward)

   815 ; res@pmLegendParallelPosF   =  0.18

   816   res@pmLegendParallelPosF   =  0.23  ;(rightward)

   818 ; res@lgPerimOn             = False

   819   res@lgLabelFontHeightF     = 0.015

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

   822 ;************************************************************

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

   825   plot_data@long_name = "Leaf Area Index"

   827 ;----------------------------------------------

   828 ; time series plot : per biome

   830   do m = 0, n_biome-1

   832      plot_name = "annual_biome_"+ m

   834      wks = gsn_open_wks (plot_type,plot_name)

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

   837      res@tiMainString = title

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

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

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

   844      delete (wks)

   845      delete (plot)

   847      system("convert -trim -density 100 "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \

   848             "rm "+plot_name+"."+plot_type)

   849   end do

   851 ;--------------------------------------------

   852 ; time series plot: all biome

   854      plot_name = "annual_biome_"+ n_biome

   856      wks = gsn_open_wks (plot_type,plot_name)

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

   859      res@tiMainString = title

   861      do k = 0,ntime-1

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

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

   864      end do

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

   868      delete (wks)

   869      delete (plot)

   871      system("convert -trim -density 100 "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \

   872             "rm "+plot_name+"."+plot_type)

   874   delete (plot_data)

   876 ;***************************************************************************

   877 ; write total score to file

   878 ;***************************************************************************

   880   asciiwrite("M_save.lai", M_total)

   882 ;***************************************************************************

   883 ; output plots

   884 ;***************************************************************************

   885   output_dir = model_name+"/lai"

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

   888 ;***************************************************************************

   889 end

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