clamp: co2/23.lines.ncl@4be95183fbcd

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

     2 ; add panel plot to 22.lines.ncl

     3 ; add zone plot to 21.lines.ncl

     4 ; ***********************************************

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

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

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

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

     9 begin

    10 ;************************************************

    11 ; read in data: observed

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

    13  diri  = "/fis/cgd/cseg/people/jeff/clamp_data/co2/"

    14  fili  = "co2_globalView_98.nc"

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

    16  val   = g->CO2_SEAS

    17  lon   = g->LON

    18  lat   = g->LAT

    19  sta   = chartostring(g->STATION)

    20  delete (g)

    22 ;print (sta(0))

    24  ncase = dimsizes(lat)

    25 ;print (ncase)

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

    28 ; get only the lowest level at each station

    29 ;**************************************************************

    30  lat_tmp = lat

    31  lat_tmp@_FillValue = 1.e+36

    33  do n = 0,ncase-1

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

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

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

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

    38        end if

    39        delete (indexes)

    40     end if

    41  end do

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

    44 ;print (dimsizes(indexes))

    45 ;print (indexes)

    47  lat_ob = lat(indexes)

    48  lon_ob = lon(indexes)

    49  val_ob = val(indexes,:)

    50 ;printVarSummary (val_ob)

    51 ;print (lat_ob +"/"+lon_ob)

    53 ;************************************************

    54 ; read in model data

    55 ;************************************************

    56   diri2 = "/fis/cgd/cseg/people/jeff/clamp_data/model/"

    57 ; fili2 = "b30.061m_401_425_MONS_climo.nc"

    58   fili2 = "b30.061n_1995-2004_MONS_climo.nc"

    60   g     = addfile(diri2+fili2,"r")

    61   x     = g->CO2

    62   xi    = g->lon

    63   yi    = g->lat

    64   xdim  = dimsizes(x)

    65   nlev  = xdim(1)

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

    67 ; printVarSummary (y)

    69 ; get the co2 at the lowest level

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

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

    73 ; Model_units [=] kgCO2 / kgDryAir

    74 ; 28.966 = molecular weight of dry air

    75 ; 44.       = molecular weight of CO2

    76 ; u mol = 1e-6 mol

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

    79   y      = y * factor

    81   y@_FillValue = 1.e36

    82   y@units      = "u mol/mol"

    83 ; y = where(y0 .lt. 287.,y@_FillValue,y)

    84 ; printVarSummary (y)

    85 ; print (min(y)+"/"+max(y))

    87 ; interpolate into observed station

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

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

    91   print (yi(0))

    92   yi(0) = -90.

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

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

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

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

   100   val_model_0 = val_model

   101 ; printVarSummary (val_model)

   102 ; print (min(val_model)+"/"+max(val_model))

   104 ; remove annual mean

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

   106 ; print (min(val_model)+"/"+max(val_model))

   108   mon = ispan(1,12,1)

   109   mon@long_name = "month"

   111   plot_type = "ps"

   112   plot_type_new = "png"

   114   res                   = True                      ; plot mods desired

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

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

   118 ;------------------------------------------------------------------

   119 ; Add a boxed legend using the more simple method, which won't have

   120 ; vertical lines going through the markers.

   122   res@pmLegendDisplayMode    = "Always"

   123 ; res@pmLegendWidthF         = 0.1

   124   res@pmLegendWidthF         = 0.08

   125   res@pmLegendHeightF        = 0.05

   126 ; res@pmLegendOrthogonalPosF = -1.17

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

   128   res@pmLegendOrthogonalPosF = -0.30  ;(downward)

   130 ; res@pmLegendParallelPosF   =  0.18

   131   res@pmLegendParallelPosF   =  0.23  ;(rightward)

   133 ; res@lgPerimOn             = False

   134   res@lgLabelFontHeightF     = 0.015

   135 ; res@xyExplicitLegendLabels = (/"model_b30.061m","observed"/)

   136   res@xyExplicitLegendLabels = (/"model_b30.061n","observed"/)

   137 ;-------------------------------------------------------------------

   139   nzone = 4

   140   do z = 0,nzone-1

   142   if (z .eq. 0) then

   143 ;    maximum score for the zone, 60N-90N

   144      score_max = 5.0

   145      zone = "60N-90N"

   146 ;    index of stations in this zone

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

   148 ;    print (ind_z)

   149 ;    print (lat_ob(ind_z)+"/"+lon_ob(ind_z))

   150 ;    print (val_ob(ind_z,:))

   151 ;    print (val_model(ind_z,:))

   152   end if

   154   if (z .eq. 1) then

   155 ;    maximum score for the zone, 30N-60N

   156      score_max = 5.0

   157      zone = "30N-60N"

   158 ;    index of stations in this zone

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

   160 ;    print (ind_z)

   161 ;    print (lat_ob(ind_z)+"/"+lon_ob(ind_z))

   162 ;    print (val_ob(ind_z,:))

   163 ;    print (val_model(ind_z,:))

   164   end if

   166   if (z .eq. 2) then

   167 ;    maximum score for the zone, EQ-30N

   168      score_max = 5.0

   169      zone = "EQ-30N"

   170 ;    index of stations in this zone

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

   172 ;    print (ind_z)

   173 ;    print (lat_ob(ind_z)+"/"+lon_ob(ind_z))

   174 ;    print (val_ob(ind_z,:))

   175 ;    print (val_model(ind_z,:))

   176   end if

   178   if (z .eq. 3) then

   179 ;    maximum score for the zone, 90S-EQ

   180      score_max = 5.0

   181      zone = "90S-EQ"

   182 ;    index of stations in this zone

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

   184 ;    print (ind_z)

   185 ;    print (lat_ob(ind_z)+"/"+lon_ob(ind_z))

   186 ;    print (val_ob(ind_z,:))

   187 ;    print (val_model(ind_z,:))

   188   end if

   190   npts = dimsizes(ind_z)

   191   print (npts)

   193   npts_str = ""

   194   npts_str = npts

   195 ; print (npts_str)

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

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

   200   plot_data!0 = "case"

   201   plot_data!1 = "month"

   202   plot_data!2 = "pts"

   203   plot_data@long_name   = "CO2 Seasonal"

   205   plot_data_0!0 = "month"

   206   plot_data_0!1 = "pts"

   207   plot_data_0@long_name = "CO2"

   209   do n=0,npts-1

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

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

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

   215     title1 = sta(ind_z(n))+"("+lat(ind_z(n))+","+lon(ind_z(n))+")"

   216     plot_name = sta(ind_z(n))

   218 ;   print (title1)

   219 ;   print (plot_name)

   220 ;***********************************************

   221 ; create panel plot

   222 ;***********************************************

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

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

   225     res@gsnFrame             = False           ; Do not draw plot

   226     res@gsnDraw              = False           ; Do not advance frame

   228     res@tiMainString      = title1                   ; add title

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

   231 ;   title2 = "Model b30.061m"

   232     title2 = "Model b30.061n"

   233     res@tiMainString      = title2                   ; add title

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

   236     pres                            = True           ; panel plot mods desired

   237     pres@gsnPanelYWhiteSpacePercent = 5              ; increase white space around

   238                                                      ; indiv. plots in panel

   239     pres@gsnMaximize                = True           ; fill the page

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

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

   244     system("rm "+plot_name+"."+plot_type)

   245 ;   system("rm "+plot_name+"-1."+plot_type_new)

   246 ;   system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)

   247     clear (wks)

   248  end do

   250 ; plot zone average

   252     plot_name = "All_"+npts_str

   253     title = plot_name + " in "+ zone

   255     print (title)

   256     print (plot_name)

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

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

   260     res@gsnFrame             = False                ; Do not draw plot

   261     res@gsnDraw              = False                ; Do not advance frame

   263     res@tiMainString      = title                   ; add title

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

   266 ;   title2 = "Model b30.061m"

   267     title2 = "Model b30.061n"

   268     res@tiMainString      = title2

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

   271     pres                            = True           ; panel plot mods desired

   272     pres@gsnPanelYWhiteSpacePercent = 5              ; increase white space around

   273                                                      ; indiv. plots in panel

   274     pres@gsnMaximize                = True           ; fill the page

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

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

   279     system("rm "+plot_name+"."+plot_type)

   280 ;   system("rm "+plot_name+"-1."+plot_type_new)

   281 ;   system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)

   283     clear (wks)

   284     delete (ind_z)

   285     delete (plot_data)

   286     delete (plot_data_0)

   287  end do

   288 end

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