forrest@0: ; *********************************************** forrest@0: ; using gsn_table for all forrest@0: ; combine 19.metric_plot.ncl and 24.lines.ncl forrest@0: ; *********************************************** forrest@0: load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl" forrest@0: load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl" forrest@0: load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl" forrest@0: ;************************************************ forrest@0: procedure set_line(lines:string,nline:integer,newlines:string) forrest@0: begin forrest@0: ; add line to ascci/html file forrest@0: forrest@0: nnewlines = dimsizes(newlines) forrest@0: if(nline+nnewlines-1.ge.dimsizes(lines)) forrest@0: print("set_line: bad index, not setting anything.") forrest@0: return forrest@0: end if forrest@0: lines(nline:nline+nnewlines-1) = newlines forrest@0: ; print ("lines = " + lines(nline:nline+nnewlines-1)) forrest@0: nline = nline + nnewlines forrest@0: return forrest@0: end forrest@0: ;**************************************************************************** forrest@0: forrest@0: begin forrest@0: forrest@0: plot_type = "ps" forrest@0: plot_type_new = "png" forrest@0: forrest@0: ;************************************************ forrest@0: ; read model data forrest@0: ;************************************************ forrest@0: forrest@0: ; from command line inputs forrest@0: forrest@0: ;-------------------------------------------------- forrest@0: ; edit table.html of current model for movel1_vs_model2 forrest@0: forrest@0: if (isvar("compare")) then forrest@0: html_name2 = compare+"/table.html" forrest@0: html_new2 = html_name2 +".new" forrest@0: forrest@0: ; the following only needs to execute once forrest@0: ; system("sed 1,/model_nameA/s//"+model_name+"/ "+html_name2+" > "+html_new2+";"+ \ forrest@0: ; "mv -f "+html_new2+" "+html_name2+";"+ \ forrest@0: ; "sed 1,/model_nameB/s//"+model_name+"/ "+html_name2+" > "+html_new2+";"+ \ forrest@0: ; "mv -f "+html_new2+" "+html_name2+";"+ \ forrest@0: ; "sed s#"+modeln+"#"+model_name+"# "+html_name2+" > "+html_new2+";"+ \ forrest@0: ; "mv -f "+html_new2+" "+html_name2) forrest@0: end if forrest@0: forrest@0: ;------------------------------------- forrest@0: ; edit table.html for current model forrest@0: forrest@0: html_name = model_name+"/table.html" forrest@0: html_new = html_name +".new" forrest@0: forrest@0: ; the following only needs to execute once forrest@0: ; system("sed s#model_name#"+model_name+"# "+html_name+" > "+html_new+";"+ \ forrest@0: ; "mv -f "+html_new+" "+html_name) forrest@0: ;------------------------------------------------ forrest@0: fm = addfile(dirm+film3,"r") forrest@0: forrest@0: x = fm->CO2 forrest@0: xi = fm->lon forrest@0: yi = fm->lat forrest@0: forrest@0: xdim = dimsizes(x) forrest@0: nlev = xdim(1) forrest@0: y = x(:,0,:,:) forrest@0: forrest@0: ; get co2 at the lowest level forrest@0: y = x(:,nlev-1,:,:) forrest@0: forrest@0: ; change to unit of observed (u mol/mol) forrest@0: ; Model_units [=] kgCO2 / kgDryAir forrest@0: ; 28.966 = molecular weight of dry air forrest@0: ; 44. = molecular weight of CO2 forrest@0: ; u mol = 1e-6 mol forrest@0: forrest@0: factor = (28.966/44.) * 1e6 forrest@0: y = y * factor forrest@0: forrest@0: y@_FillValue = 1.e36 forrest@0: y@units = "u mol/mol" forrest@0: ;************************************************ forrest@0: ; read data: observed forrest@0: ;************************************************ forrest@0: diri = "/fis/cgd/cseg/people/jeff/clamp_data/co2/ob/" forrest@0: fili = "co2_globalView_98.nc" forrest@0: g = addfile (diri+fili,"r") forrest@0: val = g->CO2_SEAS forrest@0: lon = g->LON forrest@0: lat = g->LAT forrest@0: sta = chartostring(g->STATION) forrest@0: delete (g) forrest@0: forrest@0: ncase = dimsizes(lat) forrest@0: ;************************************************************** forrest@0: ; get only the lowest level at each station forrest@0: ;************************************************************** forrest@0: lat_tmp = lat forrest@0: lat_tmp@_FillValue = 1.e+36 forrest@0: forrest@0: do n = 0,ncase-1 forrest@0: if (.not. ismissing(lat_tmp(n))) then forrest@0: indexes = ind(lat(n) .eq. lat .and. lon(n) .eq. lon) forrest@0: if (dimsizes(indexes) .gt. 1) then forrest@0: lat_tmp(indexes(1:)) = lat_tmp@_FillValue forrest@0: end if forrest@0: delete (indexes) forrest@0: end if forrest@0: end do forrest@0: forrest@0: indexes = ind(.not. ismissing(lat_tmp)) forrest@0: forrest@0: lat_ob = lat(indexes) forrest@0: lon_ob = lon(indexes) forrest@0: val_ob = val(indexes,:) forrest@0: ;************************************************************ forrest@0: ; interpolate model data into observed station forrest@0: ; note: model is 0-360E, 90S-90N forrest@0: ;************************************************************ forrest@0: ; to be able to handle observation at (-89.98,-24.80) forrest@0: yi(0) = -90. forrest@0: forrest@0: i = ind(lon_ob .lt. 0.) forrest@0: lon_ob(i) = lon_ob(i) + 360. forrest@0: forrest@0: yo = linint2_points_Wrap(xi,yi,y,True,lon_ob,lat_ob,0) forrest@0: forrest@0: val_model = yo(pts|:,time|:) forrest@0: val_model_0 = val_model forrest@0: ;************************************************************ forrest@0: ; remove annual mean forrest@0: ;************************************************************ forrest@0: val_model = val_model - conform(val_model,dim_avg(val_model),0) forrest@0: forrest@0: ;******************************************************************* forrest@0: ; res for station line plot forrest@0: ;******************************************************************* forrest@0: ; for x-axis in xyplot forrest@0: mon = ispan(1,12,1) forrest@0: mon@long_name = "month" forrest@0: forrest@0: res = True ; plot mods desired forrest@0: res@xyLineThicknesses = (/2.0,2.0,2.0/) ; make 2nd lines thicker forrest@0: res@xyLineColors = (/"red","black"/) ; change line color forrest@0: forrest@0: ; Add a boxed legend using the more simple method forrest@0: res@pmLegendDisplayMode = "Always" forrest@0: ; res@pmLegendWidthF = 0.1 forrest@0: res@pmLegendWidthF = 0.08 forrest@0: res@pmLegendHeightF = 0.06 forrest@0: ; res@pmLegendOrthogonalPosF = -1.17 forrest@0: ; res@pmLegendOrthogonalPosF = -1.00 ;(downward) forrest@0: res@pmLegendOrthogonalPosF = -0.30 ;(downward) forrest@0: forrest@0: ; res@pmLegendParallelPosF = 0.18 forrest@0: res@pmLegendParallelPosF = 0.23 ;(rightward) forrest@0: forrest@0: ; res@lgPerimOn = False forrest@0: res@lgLabelFontHeightF = 0.015 forrest@0: res@xyExplicitLegendLabels = (/model_name,"observed"/) forrest@0: ;************************************************************ forrest@0: ; number of latitude zone forrest@0: ;************************************************************ forrest@0: nzone = 4 forrest@0: forrest@0: ; saving data for zone forrest@0: ; number of rows for zone table (with data) forrest@0: nrow_zone = nzone forrest@0: forrest@0: ; number of columns for zone table forrest@0: ncol_zone = 6 forrest@0: forrest@0: text4 = new((/nrow_zone,ncol_zone/),string) forrest@0: forrest@0: do z = 0,nzone-1 forrest@0: forrest@0: if (z .eq. 0) then forrest@0: zone = "60N-90N" forrest@0: score_max = 5.0 forrest@0: ind_z = ind(lat_ob .ge. 60.) forrest@0: end if forrest@0: forrest@0: if (z .eq. 1) then forrest@0: zone = "30N-60N" forrest@0: score_max = 5.0 forrest@0: ind_z = ind(lat_ob .ge. 30. .and. lat_ob .lt. 60.) forrest@0: end if forrest@0: forrest@0: if (z .eq. 2) then forrest@0: zone = "EQ-30N" forrest@0: score_max = 5.0 forrest@0: ind_z = ind(lat_ob .ge. 0. .and. lat_ob .lt. 30.) forrest@0: end if forrest@0: forrest@0: if (z .eq. 3) then forrest@0: zone = "90S-EQ" forrest@0: score_max = 5.0 forrest@0: ind_z = ind(lat_ob .lt. 0. ) forrest@0: end if forrest@0: forrest@0: npts = dimsizes(ind_z) forrest@0: forrest@0: ;------------------------------------------------------ forrest@0: ; for metric table computation forrest@0: amp_ob = new((/npts/),float) forrest@0: amp_model = new((/npts/),float) forrest@0: forrest@0: amp_ratio_sta = new((/npts/),float) forrest@0: ccr_sta = new((/npts/),float) forrest@0: M_sta = new((/npts/),float) forrest@0: score_sta = new((/npts/),float) forrest@0: ;----------------------------------------------------- forrest@0: ; for station line plot forrest@0: npts_str = "" forrest@0: npts_str = npts forrest@0: forrest@0: plot_data = new((/2,12,npts/),float) forrest@0: plot_data_0 = new((/12,npts/),float) forrest@0: forrest@0: plot_data!0 = "case" forrest@0: plot_data!1 = "month" forrest@0: plot_data!2 = "pts" forrest@0: plot_data@long_name = "CO2 Seasonal" forrest@0: forrest@0: plot_data_0!0 = "month" forrest@0: plot_data_0!1 = "pts" forrest@0: plot_data_0@long_name = "CO2" forrest@0: ;-------------------------------------------------------------------------- forrest@0: do n=0,npts-1 forrest@0: forrest@0: amp_ob(n) = max(val_ob(ind_z(n),:)) - min(val_ob(ind_z(n),:)) forrest@0: amp_model(n) = max(val_model(ind_z(n),:)) - min(val_model(ind_z(n),:)) forrest@0: forrest@0: amp_ratio_sta(n) = amp_model(n)/amp_ob(n) forrest@0: ccr_sta(n) = esccr(val_ob(ind_z(n),:),val_model(ind_z(n),:),0) forrest@0: M_sta(n) = 1.-abs(amp_ratio_sta(n)-1.) forrest@0: score_sta(n) = (ccr_sta(n)*ccr_sta(n) + M_sta(n))*0.5 * score_max forrest@0: forrest@0: ;---------------------------------------------------------------------- forrest@0: ; for station line plot forrest@0: forrest@0: plot_data(0,:,n) = (/val_model(ind_z(n),:)/) forrest@0: plot_data(1,:,n) = (/val_ob(ind_z(n),:)/) forrest@0: forrest@0: plot_data_0(:,n) = (/val_model_0(ind_z(n),:)/) forrest@0: forrest@0: plot_name = sta(ind_z(n)) forrest@0: title = plot_name+"("+lat(ind_z(n))+","+lon(ind_z(n))+")" forrest@0: forrest@0: wks = gsn_open_wks (plot_type,plot_name) ; open workstation forrest@0: ;------------------------------------------ forrest@0: ; for panel plot forrest@0: forrest@0: plot=new(2,graphic) ; create graphic array forrest@0: res@gsnFrame = False ; Do not draw plot forrest@0: res@gsnDraw = False ; Do not advance frame forrest@0: forrest@0: pres = True ; panel plot mods desired forrest@0: pres@gsnPanelYWhiteSpacePercent = 5 ; increase white space around forrest@0: ; indiv. plots in panel forrest@0: pres@gsnMaximize = True ; fill the page forrest@0: ;------------------------------------------ forrest@0: res@tiMainString = title ; add title forrest@0: forrest@0: plot(0)=gsn_csm_xy(wks,mon,plot_data(:,:,n),res) ; create plot 1 forrest@0: forrest@0: plot(1)=gsn_csm_xy(wks,mon,plot_data_0(:,n),res) ; create plot 2 forrest@0: forrest@0: gsn_panel(wks,plot,(/2,1/),pres) ; create panel plot forrest@0: forrest@0: system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new) forrest@0: system("rm "+plot_name+"."+plot_type) forrest@0: forrest@0: clear (wks) forrest@0: ;--------------------------------------------------------------------------- forrest@0: end do forrest@0: ;------------------------------------------------------------------------- forrest@0: ; for line plot in a zone forrest@0: forrest@0: plot_name = "All_"+npts_str forrest@0: title = plot_name + " in "+ zone forrest@0: forrest@0: wks = gsn_open_wks (plot_type,plot_name) ; open workstation forrest@0: ;----------------------------------------- forrest@0: ; for panel plot forrest@0: forrest@0: plot=new(2,graphic) ; create graphic array forrest@0: res@gsnFrame = False ; Do not draw plot forrest@0: res@gsnDraw = False ; Do not advance frame forrest@0: forrest@0: pres = True ; panel plot mods desired forrest@0: pres@gsnPanelYWhiteSpacePercent = 5 ; increase white space around forrest@0: ; indiv. plots in panel forrest@0: pres@gsnMaximize = True ; fill the page forrest@0: ;----------------------------------------- forrest@0: res@tiMainString = title ; add title forrest@0: forrest@0: plot(0) = gsn_csm_xy (wks,mon,dim_avg_Wrap(plot_data),res) ; create plot 1 forrest@0: forrest@0: plot(1) = gsn_csm_xy (wks,mon,dim_avg_Wrap(plot_data_0),res) ; create plot 2 forrest@0: forrest@0: gsn_panel(wks,plot,(/2,1/),pres) ; create panel plot forrest@0: forrest@0: system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new) forrest@0: system("rm "+plot_name+"."+plot_type) forrest@0: forrest@0: clear (wks) forrest@0: ; delete (ind_z) forrest@0: delete (plot_data) forrest@0: delete (plot_data_0) forrest@0: ;--------------------------------------------------------------------------- forrest@0: ; values saved for zone table forrest@0: forrest@0: amp_ratio_zone = avg(amp_ratio_sta) forrest@0: ccr_zone = avg(ccr_sta) forrest@0: M_zone = 1.- (sum(abs(amp_model-amp_ob)/(amp_model+amp_ob))/npts) forrest@0: score_zone = (ccr_zone*ccr_zone + M_zone)*0.5 * score_max forrest@0: forrest@0: text4(z,0) = zone forrest@0: text4(z,1) = sprintf("%.0f", npts) forrest@0: text4(z,2) = sprintf("%.2f", amp_ratio_zone) forrest@0: text4(z,3) = sprintf("%.2f", ccr_zone) forrest@0: text4(z,4) = sprintf("%.2f", M_zone) forrest@0: text4(z,5) = sprintf("%.2f", score_zone) forrest@0: forrest@0: ;******************************************************************* forrest@0: ; html table -- station forrest@0: ;******************************************************************* forrest@0: output_html = "score+line_"+zone+".html" forrest@0: forrest@0: header = (/"" \ forrest@0: ,"" \ forrest@0: ,"CLAMP metrics" \ forrest@0: ,"" \ forrest@0: ,"

CO2 in Zone "+zone+": Model "+model_name+"

" \ forrest@0: /) forrest@0: footer = "" forrest@0: forrest@0: table_header = (/ \ forrest@0: "" \ forrest@0: ,"" \ forrest@0: ," " \ forrest@0: ," " \ forrest@0: ," " \ forrest@0: ," " \ forrest@0: ," " \ forrest@0: ," " \ forrest@0: ," " \ forrest@0: ,"" \ forrest@0: /) forrest@0: table_footer = "
Site NameLatitudeLongitudeAmplitude RatioCorreleration CoefM ScoreCombined Score
" forrest@0: row_header = "" forrest@0: row_footer = "" forrest@0: forrest@0: lines = new(50000,string) forrest@0: nline = 0 forrest@0: forrest@0: set_line(lines,nline,header) forrest@0: set_line(lines,nline,table_header) forrest@0: ;----------------------------------------------- forrest@0: ; row of table forrest@0: forrest@0: do n = 0,npts-1 forrest@0: set_line(lines,nline,row_header) forrest@0: forrest@0: txt0 = sta(ind_z(n)) forrest@0: txt1 = sprintf("%5.2f", (/lat(ind_z(n))/)) forrest@0: txt2 = sprintf("%5.2f", (/lon(ind_z(n))/)) forrest@0: txt3 = sprintf("%5.2f", (/amp_ratio_sta(n)/)) forrest@0: txt4 = sprintf("%5.2f", (/ccr_sta(n)/)) forrest@0: txt5 = sprintf("%5.2f", (/M_sta(n)/)) forrest@0: txt6 = sprintf("%5.2f", (/score_sta(n)/)) forrest@0: forrest@0: set_line(lines,nline,""+txt0+"") forrest@0: set_line(lines,nline,""+txt1+"") forrest@0: set_line(lines,nline,""+txt2+"") forrest@0: set_line(lines,nline,""+txt3+"") forrest@0: set_line(lines,nline,""+txt4+"") forrest@0: set_line(lines,nline,""+txt5+"") forrest@0: set_line(lines,nline,""+txt6+"") forrest@0: forrest@0: set_line(lines,nline,row_footer) forrest@0: end do forrest@0: forrest@0: ; last row, summary forrest@0: set_line(lines,nline,row_header) forrest@0: forrest@0: txt0 = "All_"+sprintf("%.0f", (/npts/)) forrest@0: txt1 = "-" forrest@0: txt2 = "-" forrest@0: txt3 = sprintf("%5.2f", (/amp_ratio_zone/)) forrest@0: txt4 = sprintf("%5.2f", (/ccr_zone/)) forrest@0: txt5 = sprintf("%5.2f", (/M_zone/)) forrest@0: txt6 = sprintf("%5.2f", (/score_zone/)) forrest@0: forrest@0: set_line(lines,nline,""+txt0+"") forrest@0: set_line(lines,nline,""+txt1+"") forrest@0: set_line(lines,nline,""+txt2+"") forrest@0: set_line(lines,nline,""+txt3+"") forrest@0: set_line(lines,nline,""+txt4+"") forrest@0: set_line(lines,nline,""+txt5+"") forrest@0: set_line(lines,nline,""+txt6+"") forrest@0: forrest@0: set_line(lines,nline,row_footer) forrest@0: ;----------------------------------------------- forrest@0: set_line(lines,nline,table_footer) forrest@0: set_line(lines,nline,footer) forrest@0: forrest@0: ; Now write to an HTML file. forrest@0: idx = ind(.not.ismissing(lines)) forrest@0: if(.not.any(ismissing(idx))) then forrest@0: asciiwrite(output_html,lines(idx)) forrest@0: else forrest@0: print ("error?") forrest@0: end if forrest@0: delete (idx) forrest@0: ;----------------------------------------------------------------- forrest@0: forrest@0: delete (ind_z) forrest@0: delete (amp_model) forrest@0: delete (amp_ob) forrest@0: delete (amp_ratio_sta) forrest@0: delete (ccr_sta) forrest@0: delete (M_sta) forrest@0: delete (score_sta) forrest@0: clear (wks) forrest@0: end do forrest@0: forrest@0: ;******************************************************************* forrest@0: ; html table -- zone forrest@0: ;******************************************************************* forrest@0: output_html = "score+line_vs_ob.html" forrest@0: forrest@0: header = (/"" \ forrest@0: ,"" \ forrest@0: ,"CLAMP metrics" \ forrest@0: ,"" \ forrest@0: ,"

CO2 in Latitude Zone: Model "+model_name+"

" \ forrest@0: /) forrest@0: footer = "" forrest@0: forrest@0: delete (table_header) forrest@0: table_header = (/ \ forrest@0: "" \ forrest@0: ,"" \ forrest@0: ," " \ forrest@0: ," " \ forrest@0: ," " \ forrest@0: ," " \ forrest@0: ," " \ forrest@0: ," " \ forrest@0: ,"" \ forrest@0: /) forrest@0: table_footer = "
ZoneNumber of SiteAmplitide RatioCorreleration CoefM ScoreCombined Score
" forrest@0: row_header = "" forrest@0: row_footer = "" forrest@0: forrest@0: lines = new(50000,string) forrest@0: nline = 0 forrest@0: forrest@0: set_line(lines,nline,header) forrest@0: set_line(lines,nline,table_header) forrest@0: ;----------------------------------------------- forrest@0: ;row of table forrest@0: forrest@0: do n = 0,nrow_zone-1 forrest@0: set_line(lines,nline,row_header) forrest@0: forrest@0: set_line(lines,nline,""+text4(n,0)+"") forrest@0: set_line(lines,nline,""+text4(n,1)+"") forrest@0: set_line(lines,nline,""+text4(n,2)+"") forrest@0: set_line(lines,nline,""+text4(n,3)+"") forrest@0: set_line(lines,nline,""+text4(n,4)+"") forrest@0: set_line(lines,nline,""+text4(n,5)+"") forrest@0: forrest@0: set_line(lines,nline,row_footer) forrest@0: end do forrest@0: forrest@0: ; for the last row forrest@0: forrest@0: txt0 = "All" forrest@0: txt1 = sum(stringtofloat(text4(0:3,1))) forrest@0: txt2 = "-" forrest@0: txt3 = "-" forrest@0: txt4 = "-" forrest@0: txt5 = sum(stringtofloat(text4(0:3,5))) forrest@0: forrest@0: set_line(lines,nline,row_header) forrest@0: forrest@0: set_line(lines,nline,""+txt0+"") forrest@0: set_line(lines,nline,""+txt1+"") forrest@0: set_line(lines,nline,""+txt2+"") forrest@0: set_line(lines,nline,""+txt3+"") forrest@0: set_line(lines,nline,""+txt4+"") forrest@0: set_line(lines,nline,""+txt5+"") forrest@0: forrest@0: set_line(lines,nline,row_footer) forrest@0: ;----------------------------------------------- forrest@0: set_line(lines,nline,table_footer) forrest@0: set_line(lines,nline,footer) forrest@0: forrest@0: ; Now write to an HTML file. forrest@0: idx = ind(.not.ismissing(lines)) forrest@0: if(.not.any(ismissing(idx))) then forrest@0: asciiwrite(output_html,lines(idx)) forrest@0: else forrest@0: print ("error?") forrest@0: end if forrest@0: delete (idx) forrest@0: ;-------------------------------------------------------------------------- forrest@0: forrest@0: M_co2 = txt5 forrest@0: forrest@0: if (isvar("compare")) then forrest@0: system("sed 1,/M_co2/s//"+M_co2+"/ "+html_name2+" > "+html_new2+";"+ \ forrest@0: "mv -f "+html_new2+" "+html_name2) forrest@0: end if forrest@0: forrest@0: system("sed s#M_co2#"+M_co2+"# "+html_name+" > "+html_new+";"+ \ forrest@0: "mv -f "+html_new+" "+html_name) forrest@0: ;*************************************************************************** forrest@0: ; output plots forrest@0: ;*************************************************************************** forrest@0: output_dir = model_name+"/co2" forrest@0: forrest@0: system("mv *.png *.html " + output_dir) forrest@0: ;*************************************************************************** forrest@0: end