clamp: all/03.co2.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 load dirscript + "taylor_diagram.ncl"

     6 ;************************************************

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

     8 begin

     9 ; add line to ascci/html file

    11   nnewlines = dimsizes(newlines)

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

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

    14     return

    15   end if

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

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

    18   nline = nline + nnewlines

    19   return

    20 end

    21 ;****************************************************************************

    23 begin

    25   plot_type = "ps"

    26   plot_type_new = "png"

    28 ;-----------------------------------------------------

    29 ; edit table.html of current model for movel1_vs_model2

    31   if (isvar("compare")) then

    32      html_name2 = compare+"/table.html"

    33      html_new2  = html_name2 +".new"

    34   end if

    36 ;------------------------------------------------------

    37 ; edit table.html for current model

    39   html_name = model_name+"/table.html"

    40   html_new  = html_name +".new"

    42 ;------------------------------------------------------

    43 ; read model data

    45   fm    = addfile(dirm+film3,"r")

    47   x     = fm->CO2

    48   xi    = fm->lon

    49   yi    = fm->lat

    51   delete (fm)

    53   xdim  = dimsizes(x)

    54   nlev  = xdim(1)

    55   y     = x(:,0,:,:)

    57 ; get co2 at the lowest level

    58   y     = x(:,nlev-1,:,:)

    60 ; change to unit of observed (u mol/mol)

    61 ; Model_units [=] kgCO2 / kgDryAir

    62 ; 28.966 = molecular weight of dry air

    63 ; 44.       = molecular weight of CO2

    64 ; u mol = 1e-6 mol

    66   factor = (28.966/44.) * 1e6

    67   y      = y * factor

    69   y@_FillValue = 1.e36

    70   y@units      = "u mol/mol"

    72 ;************************************************

    73 ; read data: observed

    74 ;************************************************

    75   diri  = diro + "co2/"

    76   fili  = "co2_globalView_98.nc"

    77   g     = addfile (diri+fili,"r")

    79   val   = g->CO2_SEAS

    80   lon   = g->LON

    81   lat   = g->LAT

    82   sta   = chartostring(g->STATION)

    84   delete (g)

    86   ncase = dimsizes(lat)

    87 ;**************************************************************

    88 ; get only the lowest level at each station

    89 ;**************************************************************

    90   lat_tmp = lat

    91   lat_tmp@_FillValue = 1.e+36

    93   do n = 0,ncase-1

    94      if (.not. ismissing(lat_tmp(n))) then

    95         indexes = ind(lat(n) .eq. lat .and. lon(n) .eq. lon)

    96         if (dimsizes(indexes) .gt. 1) then

    97            lat_tmp(indexes(1:)) = lat_tmp@_FillValue

    98         end if

    99         delete (indexes)

   100      end if

   101   end do

   103   indexes = ind(.not. ismissing(lat_tmp))

   105   lat_ob = lat(indexes)

   106   lon_ob = lon(indexes)

   107   val_ob = val(indexes,:)

   108 ;************************************************************

   109 ; interpolate model data into observed station

   110 ; note: model is 0-360E, 90S-90N

   111 ;************************************************************

   112 ; to be able to handle observation at (-89.98,-24.80)

   113   yi(0) = -90.

   115   i = ind(lon_ob .lt. 0.)

   116   lon_ob(i) = lon_ob(i) + 360.

   118   yo = linint2_points_Wrap(xi,yi,y,True,lon_ob,lat_ob,0)

   120   val_model = yo(pts|:,time|:)

   121   val_model_0 = val_model

   122 ;************************************************************

   123 ; remove annual mean

   124 ;************************************************************

   125   val_model = val_model - conform(val_model,dim_avg(val_model),0)

   127 ;*******************************************************************

   128 ; res for station line plot

   129 ;*******************************************************************

   130 ; for x-axis in xyplot

   131   mon = ispan(1,12,1)

   132   mon@long_name = "month"

   134   res                   = True                      ; plot mods desired

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

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

   138 ; Add a boxed legend using the more simple method

   139   res@pmLegendDisplayMode    = "Always"

   140 ; res@pmLegendWidthF         = 0.1

   141   res@pmLegendWidthF         = 0.08

   142   res@pmLegendHeightF        = 0.06

   143 ; res@pmLegendOrthogonalPosF = -1.17

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

   145   res@pmLegendOrthogonalPosF = -0.30  ;(downward)

   147 ; res@pmLegendParallelPosF   =  0.18

   148   res@pmLegendParallelPosF   =  0.23  ;(rightward)

   150 ; res@lgPerimOn             = False

   151   res@lgLabelFontHeightF     = 0.015

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

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

   154 ; number of latitude zone

   155 ;************************************************************

   156   nzone = 4

   158 ; saving data for zone

   159 ; number of rows for zone table (with data)

   160   nrow_zone = nzone

   162 ; number of columns for zone table

   163   ncol_zone = 7

   165   text = new((/nrow_zone,ncol_zone/),string)

   167 do z = 0,nzone-1

   169   if (z .eq. 0) then

   170      zone = "60N-90N"

   171      score_max = 5.0

   172      ind_z = ind(lat_ob .ge. 60.)

   173   end if

   175   if (z .eq. 1) then

   176      zone = "30N-60N"

   177      score_max = 5.0

   178      ind_z = ind(lat_ob .ge. 30. .and. lat_ob .lt. 60.)

   179   end if

   181   if (z .eq. 2) then

   182      zone = "EQ-30N"

   183      score_max = 5.0

   184      ind_z = ind(lat_ob .ge. 0. .and. lat_ob .lt. 30.)

   185   end if

   187   if (z .eq. 3) then

   188      zone = "90S-EQ"

   189      score_max = 5.0

   190      ind_z = ind(lat_ob .lt. 0. )

   191   end if

   193   npts = dimsizes(ind_z)

   195 ;------------------------------------------------------

   196 ; for metric table computation

   197   amp_ob        = new((/npts/),float)

   198   amp_model     = new((/npts/),float)

   200   amp_ratio_sta = new((/npts/),float)

   201   ccr_sta       = new((/npts/),float)

   202   M_sta         = new((/npts/),float)

   203   score_sta     = new((/npts/),float)

   205   var_sta       = new((/npts/),float)

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

   207 ; for station line plot

   208   npts_str = ""

   209   npts_str = npts

   211   plot_data   = new((/2,12,npts/),float)

   212   plot_data_0 = new((/12,npts/),float)

   214   plot_data!0 = "case"

   215   plot_data!1 = "month"

   216   plot_data!2 = "pts"

   217   plot_data@long_name   = "CO2 Annual Cycle (ppm)"

   219   plot_data_0!0 = "month"

   220   plot_data_0!1 = "pts"

   221   plot_data_0@long_name = "CO2 (ppm)"

   222 ;--------------------------------------------------------------------------

   223   do n=0,npts-1

   225      amp_ob(n)    = max(val_ob(ind_z(n),:)) - min(val_ob(ind_z(n),:))

   226      amp_model(n) = max(val_model(ind_z(n),:)) - min(val_model(ind_z(n),:))

   228      amp_ratio_sta(n) = amp_model(n)/amp_ob(n)

   229      ccr_sta(n) = esccr(val_ob(ind_z(n),:),val_model(ind_z(n),:),0)

   230      M_sta(n) = 1.-abs(amp_ratio_sta(n)-1.)

   231      score_sta(n) = (ccr_sta(n)*ccr_sta(n) + M_sta(n))*0.5 * score_max

   233      var_model = stddev(val_model(ind_z(n),:))

   234      var_ob    = stddev(val_ob(ind_z(n),:))

   235      var_sta(n)= var_model/var_ob

   236 ;----------------------------------------------------------------------

   237 ; for station line plot

   239      plot_data(0,:,n) = (/val_model(ind_z(n),:)/)

   240      plot_data(1,:,n) = (/val_ob(ind_z(n),:)/)

   242      plot_data_0(:,n) = (/val_model_0(ind_z(n),:)/)

   244      plot_name = sta(ind_z(n))

   245      title = plot_name+" ("+lat(ind_z(n))+","+lon(ind_z(n))+")"

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

   248 ;------------------------------------------

   249 ;    for panel plot

   251      plot=new(2,graphic)                        ; create graphic array

   252      res@gsnFrame     = False                   ; Do not draw plot

   253      res@gsnDraw      = False                   ; Do not advance frame

   255      pres                            = True     ; panel plot mods desired

   256      pres@gsnPanelYWhiteSpacePercent = 5        ; increase white space around

   257                                                ; indiv. plots in panel

   258      pres@gsnMaximize                = True     ; fill the page

   259 ;------------------------------------------

   260      res@tiMainString = title                           ; add title

   262      plot(0)=gsn_csm_xy(wks,mon,plot_data(:,:,n),res)   ; create plot 1

   264      plot(1)=gsn_csm_xy(wks,mon,plot_data_0(:,n),res) ; create plot 2

   266      gsn_panel(wks,plot,(/2,1/),pres)                 ; create panel plot

   268      delete (wks)

   269      delete (plot)

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

   272             "rm "+plot_name+"."+plot_type)

   274 ;---------------------------------------------------------------------------

   275   end do

   277 ;-------------------------------------------------------------------------

   278 ; for Taylor plot in a zone

   280 ; Cases [Model]

   281   case      = (/ model_name /)

   282   nCase     = dimsizes(case )                 ; # of Cases [Cases]

   284 ; station compared

   285   var       = sta(ind_z)

   286   nVar      = dimsizes(var)                   ; # of stations

   288 ; arrays to be passed to taylor plot

   289   ratio      = new ((/nCase, nVar/),typeof(var_sta) )

   290   cc         = new ((/nCase, nVar/),typeof(ccr_sta) )

   292   ratio(0,:) = var_sta

   293   cc(0,:)    = ccr_sta

   295 ;---------------------------------

   296 ; create plot

   298   rest   = True                           ; default taylor diagram

   300   rest@Markers      = (/16/)               ; make all solid fill

   301   rest@Colors       = (/"red" /)

   302 ; rest@varLabels    = var

   303   rest@caseLabels   = case

   305   rOpts  = True

   306   rOpts@caseLabels = True

   307 ; rOpts@caseLabelsFontHeightF = 0.05

   308   rOpts@caseLabelsFontHeightF = 0.15

   310   plot_name = "taylor_diagram_"+ zone

   311   title = "CO2 Annual Cycle:  "+ zone

   312   rest@tiMainString = title

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

   315   plot = taylor_diagram(wks,ratio,cc,rest)

   317   delete (wks)

   318   delete (plot)

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

   321          "rm "+plot_name+"."+plot_type)

   323   delete (ratio)

   324   delete (cc)

   325   delete (var_sta)

   326   delete (var)

   328 ;-------------------------------------------------------------------------

   329 ; for line plot in a zone

   331   plot_name = "All_"+npts_str

   332   title = plot_name + " in "+ zone

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

   335 ;-----------------------------------------

   336 ; for panel plot

   338   plot=new(2,graphic)                        ; create graphic array

   339   res@gsnFrame     = False                   ; Do not draw plot

   340   res@gsnDraw      = False                   ; Do not advance frame

   342   pres                            = True     ; panel plot mods desired

   343   pres@gsnPanelYWhiteSpacePercent = 5        ; increase white space around

   344                                                ; indiv. plots in panel

   345   pres@gsnMaximize                = True     ; fill the page

   346 ;-----------------------------------------

   347   res@tiMainString = title                                     ; add title

   349   plot(0) = gsn_csm_xy (wks,mon,dim_avg_Wrap(plot_data),res)   ; create plot 1

   351   plot(1) = gsn_csm_xy (wks,mon,dim_avg_Wrap(plot_data_0),res) ; create plot 2

   353   gsn_panel(wks,plot,(/2,1/),pres)                 ; create panel plot

   355   delete (wks)

   356   delete (plot)

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

   359          "rm "+plot_name+"."+plot_type)

   361   delete (plot_data)

   362   delete (plot_data_0)

   363 ;---------------------------------------------------------------------------

   364 ; values saved for zone table

   366   amp_ratio_zone = avg(amp_ratio_sta)

   367   ccr_zone       = avg(ccr_sta)

   368   M_zone         = 1.- (sum(abs(amp_model-amp_ob)/(amp_model+amp_ob))/npts)

   369   score_zone     = (ccr_zone*ccr_zone + M_zone)*0.5 * score_max

   371   text(z,0) = zone

   372   text(z,1) = sprintf("%.0f", npts)

   373   text(z,2) = sprintf("%.2f", amp_ratio_zone)

   374   text(z,3) = sprintf("%.2f", ccr_zone)

   375   text(z,4) = sprintf("%.2f", M_zone)

   376   text(z,5) = sprintf("%.2f", score_zone)

   377   text(z,6) = zone

   379 ;*******************************************************************

   380 ; html table -- station

   381 ;*******************************************************************

   382   output_html = "score+line_"+zone+".html"

   384   header = (/"<HTML>" \

   385             ,"<HEAD>" \

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

   387             ,"</HEAD>" \

   388             ,"<H1>Latitude Zone "+zone+": Model "+model_name+"</H1>" \

   389             /)

   390   footer = "</HTML>"

   392   table_header = (/ \

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

   394        ,"<tr>" \

   395        ,"   <th bgcolor=DDDDDD >Site Name</th>" \

   396        ,"   <th bgcolor=DDDDDD >Latitude</th>" \

   397        ,"   <th bgcolor=DDDDDD >Longitude</th>" \

   398        ,"   <th bgcolor=DDDDDD >model vs obs.<br>amplitude ratio</th>" \

   399        ,"   <th bgcolor=DDDDDD >Correlation Coef.</th>" \

   400        ,"   <th bgcolor=DDDDDD >M Score</th>" \

   401        ,"   <th bgcolor=DDDDDD >Combined Score</th>" \

   402        ,"</tr>" \

   403        /)

   404   table_footer = "</table>"

   405   row_header = "<tr>"

   406   row_footer = "</tr>"

   408   lines = new(50000,string)

   409   nline = 0

   411   set_line(lines,nline,header)

   412   set_line(lines,nline,table_header)

   413 ;-----------------------------------------------

   414 ; row of table

   416   do n = 0,npts-1

   417      set_line(lines,nline,row_header)

   419      txt0 = sta(ind_z(n))

   420      txt1 = sprintf("%5.2f", (/lat(ind_z(n))/))

   421      txt2 = sprintf("%5.2f", (/lon(ind_z(n))/))

   422      txt3 = sprintf("%5.2f", (/amp_ratio_sta(n)/))

   423      txt4 = sprintf("%5.2f", (/ccr_sta(n)/))

   424      txt5 = sprintf("%5.2f", (/M_sta(n)/))

   425      txt6 = sprintf("%5.2f", (/score_sta(n)/))

   427      set_line(lines,nline,"<th><a href="+txt0+".png>"+txt0+"</a></th>")

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

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

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

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

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

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

   435      set_line(lines,nline,row_footer)

   436   end do

   438 ; last row, summary

   439   set_line(lines,nline,row_header)

   441   txt0 = "All_"+sprintf("%.0f", (/npts/))

   442   txt1 = "-"

   443   txt2 = "-"

   444   txt3 = sprintf("%5.2f", (/amp_ratio_zone/))

   445   txt4 = sprintf("%5.2f", (/ccr_zone/))

   446   txt5 = sprintf("%5.2f", (/M_zone/))

   447   txt6 = sprintf("%5.2f", (/score_zone/))

   449   set_line(lines,nline,"<th><a href="+txt0+".png>"+txt0+"</a></th>")

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

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

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

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

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

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

   457   set_line(lines,nline,row_footer)

   458 ;-----------------------------------------------

   459   set_line(lines,nline,table_footer)

   460   set_line(lines,nline,footer)

   462 ; Now write to an HTML file.

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

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

   465     asciiwrite(output_html,lines(idx))

   466   else

   467    print ("error?")

   468   end if

   469   delete (idx)

   470 ;-----------------------------------------------------------------

   472   delete (ind_z)

   473   delete (amp_model)

   474   delete (amp_ob)

   475   delete (amp_ratio_sta)

   476   delete (ccr_sta)

   477   delete (M_sta)

   478   delete (score_sta)

   479 end do

   481 ;*******************************************************************

   482 ; html table -- zone

   483 ;*******************************************************************

   484   output_html = "score+line_vs_ob.html"

   486   header = (/"<HTML>" \

   487             ,"<HEAD>" \

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

   489             ,"</HEAD>" \

   490             ,"<H1>CO2 Seasonal Cycle Comparisons by Latitude Zone: Model "+model_name+"</H1>" \

   491             /)

   492   footer = "</HTML>"

   494   delete (table_header)

   495   table_header = (/ \

   496         "<table border=1 cellspacing=0 cellpadding=3 width=80%>" \

   497        ,"<tr>" \

   498        ,"   <th bgcolor=DDDDDD >Zone</th>" \

   499        ,"   <th bgcolor=DDDDDD >Number of Site</th>" \

   500        ,"   <th bgcolor=DDDDDD >model vs obs.<br>amplitide ratio</th>" \

   501        ,"   <th bgcolor=DDDDDD >Correlation Coef.</th>" \

   502        ,"   <th bgcolor=DDDDDD >M Score</th>" \

   503        ,"   <th bgcolor=DDDDDD >Combined Score</th>" \

   504        ,"   <th bgcolor=DDDDDD >Taylor diagram</th>" \

   505        ,"</tr>" \

   506        /)

   507   table_footer = "</table>"

   508   row_header = "<tr>"

   509   row_footer = "</tr>"

   511   lines = new(50000,string)

   512   nline = 0

   514   set_line(lines,nline,header)

   515   set_line(lines,nline,table_header)

   516 ;-----------------------------------------------

   517 ;row of table

   519   do n = 0,nrow_zone-1

   520      set_line(lines,nline,row_header)

   522      set_line(lines,nline,"<th><a href=score+line_"+text(n,0)+".html>"+text(n,0)+"</th>")

   523      set_line(lines,nline,"<th>"+text(n,1)+"</th>")

   524      set_line(lines,nline,"<th>"+text(n,2)+"</th>")

   525      set_line(lines,nline,"<th>"+text(n,3)+"</th>")

   526      set_line(lines,nline,"<th>"+text(n,4)+"</th>")

   527      set_line(lines,nline,"<th>"+text(n,5)+"</th>")

   528      set_line(lines,nline,"<th><a href=taylor_diagram_"+text(n,6)+".png>Taylor_diagram</th>")

   530      set_line(lines,nline,row_footer)

   531   end do

   533 ; for the last row

   535      txt0 = "All"

   536      txt1 = sum(stringtofloat(text(0:3,1)))

   537      txt2 = "-"

   538      txt3 = "-"

   539      txt4 = "-"

   540      txt5 = sum(stringtofloat(text(0:3,5)))

   541      txt6 = "-"

   543      set_line(lines,nline,row_header)

   545      set_line(lines,nline,"<th>"+txt0+"</th>")

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

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

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

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

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

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

   553      set_line(lines,nline,row_footer)

   554 ;-----------------------------------------------

   555   set_line(lines,nline,table_footer)

   556   set_line(lines,nline,footer)

   558 ; Now write to an HTML file.

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

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

   561     asciiwrite(output_html,lines(idx))

   562   else

   563    print ("error?")

   564   end if

   565   delete (idx)

   566 ;--------------------------------------------------------------------------

   568   M_co2 = txt5

   570   if (isvar("compare")) then

   572      system("sed -e '1,/M_co2/s/M_co2/"+M_co2+"/' "+html_name2+" > "+html_new2+";"+ \

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

   574   end if

   576   system("sed s#M_co2#"+M_co2+"# "+html_name+" > "+html_new+";"+ \

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

   579 ;***************************************************************************

   580 ; add total score and write to file

   581 ;***************************************************************************

   582   M_total = M_co2

   584   asciiwrite("M_save.co2", M_total)

   586 ;***************************************************************************

   587 ; output plots

   588 ;***************************************************************************

   589   output_dir = model_name+"/co2"

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

   592 ;***************************************************************************

   593 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--