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

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

     2 ; add another model to plot

     3 ; add panel plot to 22.lines.ncl

     4 ; add zone plot to 21.lines.ncl

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

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

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

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

     9 ;************************************************

    10 begin

    11 ;************************************************

    12 ; read in data: observed

    13 ;************************************************

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

    15  fili  = "co2_globalView_98.nc"

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

    17  val   = g->CO2_SEAS

    18  lon   = g->LON

    19  lat   = g->LAT

    20  sta   = chartostring(g->STATION)

    21  delete (g)

    23 ;print (sta(0))

    25  ncase = dimsizes(lat)

    26 ;print (ncase)

    28 ;**************************************************************

    29 ; get only the lowest level at each station

    30 ;**************************************************************

    31  lat_tmp = lat

    32  lat_tmp@_FillValue = 1.e+36

    34  do n = 0,ncase-1

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

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

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

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

    39        end if

    40        delete (indexes)

    41     end if

    42  end do

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

    45 ;print (dimsizes(indexes))

    46 ;print (indexes)

    48  lat_ob = lat(indexes)

    49  lon_ob = lon(indexes)

    50  val_ob = val(indexes,:)

    51 ;printVarSummary (val_ob)

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

    54 ;************************************************

    55 ; read in model data

    56 ;************************************************

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

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

    59   fili3 = "b30.061m_401_425_MONS_climo_atm.nc"

    60 ;--------------------------------------------

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

    62   x1   = g->CO2

    63   xi   = g->lon

    64   yi   = g->lat

    65   delete (g)

    67   xdim  = dimsizes(x1)

    68   nlev  = xdim(1)

    69   y1     = x1(:,0,:,:)

    70 ; printVarSummary (y1)

    72 ; get the co2 at the lowest level

    73   y1     = x1(:,nlev-1,:,:)

    74   delete (x1)

    75 ;---------------------------------------------

    76   g     = addfile(diri2+fili3,"r")

    77   x2    = g->CO2

    78   delete (g)

    79   y2     = x2(:,0,:,:)

    80   y2     = x2(:,nlev-1,:,:)

    81   delete (x2)

    82 ;---------------------------------------------

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

    84 ; Model_units [=] kgCO2 / kgDryAir

    85 ; 28.966 = molecular weight of dry air

    86 ; 44.       = molecular weight of CO2

    87 ; u mol = 1e-6 mol

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

    90 ;---------------------------------------------

    91   y1 = y1 * factor

    92   y1@_FillValue = 1.e36

    93   y1@units      = "u mol/mol"

    94 ; y1 = where(y0 .lt. 287.,1y@_FillValue,y1)

    95 ; printVarSummary (y1)

    96 ; print (min(y1)+"/"+max(y1))

    97 ;---------------------------------------------

    98   y2 = y2 * factor

    99   y2@_FillValue = 1.e36

   100   y2@units      = "u mol/mol"

   101 ;---------------------------------------------

   102 ; interpolate into observed station

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

   104 ;       ob    is -180-180, 90S-90N

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

   107   print (yi(0))

   108   yi(0) = -90.

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

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

   112 ;----------------------------------------------------------------

   113   yo = linint2_points_Wrap(xi,yi,y1,True,lon_ob,lat_ob,0)

   115   val_model1 = yo(pts|:,time|:)

   116   val_model1_0 = val_model1

   117 ; printVarSummary (val_model1)

   118 ; print (min(val_model1)+"/"+max(val_model1))

   120 ; remove annual mean

   121   val_model1 = val_model1 - conform(val_model1,dim_avg(val_model1),0)

   122 ; print (min(val_model1)+"/"+max(val_model1))

   123   delete (yo)

   124 ;-----------------------------------------------------------------

   125   yo = linint2_points_Wrap(xi,yi,y2,True,lon_ob,lat_ob,0)

   127   val_model2 = yo(pts|:,time|:)

   128   val_model2_0 = val_model2

   129 ; printVarSummary (val_model2)

   130 ; print (min(val_model2)+"/"+max(val_model2))

   132 ; remove annual mean

   133   val_model2 = val_model2 - conform(val_model2,dim_avg(val_model2),0)

   134 ; print (min(val_model2)+"/"+max(val_model2))

   135   delete (yo)

   136 ;-----------------------------------------------------------------

   137 ; for x-axis in xyplot

   138   mon = ispan(1,12,1)

   139   mon@long_name = "month"

   140 ;-----------------------------------------------------------------

   141   plot_type = "ps"

   142   plot_type_new = "png"

   144   res                   = True                      ; plot mods desired

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

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

   147 ;------------------------------------------------------------------

   148 ; Add a boxed legend using the more simple method

   150   res@pmLegendDisplayMode    = "Always"

   151 ; res@pmLegendWidthF         = 0.1

   152   res@pmLegendWidthF         = 0.08

   153   res@pmLegendHeightF        = 0.06

   154 ; res@pmLegendOrthogonalPosF = -1.17

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

   156   res@pmLegendOrthogonalPosF = -0.30  ;(downward)

   158 ; res@pmLegendParallelPosF   =  0.18

   159   res@pmLegendParallelPosF   =  0.23  ;(rightward)

   161 ; res@lgPerimOn             = False

   162   res@lgLabelFontHeightF     = 0.015

   163   res@xyExplicitLegendLabels = (/"b30.061n","b30.061m","observed"/)

   164 ;-------------------------------------------------------------------

   166   nzone = 4

   167   do z = 0,nzone-1

   169   if (z .eq. 0) then

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

   171      score_max = 5.0

   172      zone = "60N-90N"

   173 ;    index of stations in this zone

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

   175 ;    print (ind_z)

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

   177 ;    print (val_ob(ind_z,:))

   178 ;    print (val_model(ind_z,:))

   179   end if

   181   if (z .eq. 1) then

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

   183      score_max = 5.0

   184      zone = "30N-60N"

   185 ;    index of stations in this zone

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

   187 ;    print (ind_z)

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

   189 ;    print (val_ob(ind_z,:))

   190 ;    print (val_model(ind_z,:))

   191   end if

   193   if (z .eq. 2) then

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

   195      score_max = 5.0

   196      zone = "EQ-30N"

   197 ;    index of stations in this zone

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

   199 ;    print (ind_z)

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

   201 ;    print (val_ob(ind_z,:))

   202 ;    print (val_model(ind_z,:))

   203   end if

   205   if (z .eq. 3) then

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

   207      score_max = 5.0

   208      zone = "90S-EQ"

   209 ;    index of stations in this zone

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

   211 ;    print (ind_z)

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

   213 ;    print (val_ob(ind_z,:))

   214 ;    print (val_model(ind_z,:))

   215   end if

   217   npts = dimsizes(ind_z)

   218   print (npts)

   220   npts_str = ""

   221   npts_str = npts

   222 ; print (npts_str)

   224   plot_data   = new((/3,12,npts/),float)

   225   plot_data_0 = new((/2,12,npts/),float)

   227   plot_data!0 = "case"

   228   plot_data!1 = "month"

   229   plot_data!2 = "pts"

   230   plot_data@long_name   = "CO2 Seasonal"

   232   plot_data_0!0 = "month"

   233   plot_data_0!1 = "pts"

   234   plot_data_0@long_name = "CO2"

   236   do n=0,npts-1

   238     plot_data(0,:,n) = (/val_model1(ind_z(n),:)/)

   239     plot_data(1,:,n) = (/val_model2(ind_z(n),:)/)

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

   242     plot_data_0(0,:,n) = (/val_model1_0(ind_z(n),:)/)

   243     plot_data_0(1,:,n) = (/val_model2_0(ind_z(n),:)/)

   245 ;*****************************************************

   246 ; create plot with 3 lines for each observed station

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

   248     plot_name = sta(ind_z(n))

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

   250 ;   print (title)

   251 ;   print (plot_name)

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

   254 ;---------------------------------------------------------------------------

   255 ;   for panel plot

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

   258     res@gsnFrame     = False                   ; Do not draw plot

   259     res@gsnDraw      = False                   ; Do not advance frame

   261     pres                            = True     ; panel plot mods desired

   262     pres@gsnPanelYWhiteSpacePercent = 5        ; increase white space around

   263                                                ; indiv. plots in panel

   264     pres@gsnMaximize                = True     ; fill the page

   265 ;---------------------------------------------------------------------------

   266     res@tiMainString = title                           ; add title

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

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

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

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

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

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

   277     clear (wks)

   278  end do

   279 ;*****************************************************************************

   280 ; plot zone average

   281 ;*****************************************************************************

   282     plot_name = "All_"+npts_str

   283     title = plot_name + " in "+ zone

   284 ;   print (title)

   285 ;   print (plot_name)

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

   288 ;---------------------------------------------------------------------------

   289 ;   for panel plot

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

   292     res@gsnFrame     = False                   ; Do not draw plot

   293     res@gsnDraw      = False                   ; Do not advance frame

   295     pres                            = True     ; panel plot mods desired

   296     pres@gsnPanelYWhiteSpacePercent = 5        ; increase white space around

   297                                                ; indiv. plots in panel

   298     pres@gsnMaximize                = True     ; fill the page

   299 ;---------------------------------------------------------------------------

   300     res@tiMainString = title                                     ; add title

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

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

   306     xx = dim_avg_Wrap(plot_data_0)

   307     plot(1) = gsn_csm_xy (wks,mon,xx(0,:),res) ; create plot 2

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

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

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

   314     clear (wks)

   315     delete (ind_z)

   316     delete (plot_data)

   317     delete (plot_data_0)

   318  end do

   319 end

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