Mercurial Mercurial > repos > public_software > clamp / file revision
summary | shortlog | changelog | graph | tags | bookmarks | branches | files | changeset | file | latest | revisions | annotate | diff | comparison | raw | help
co2/22.lines.ncl
author Forrest Hoffman <forrest@climatemodeling.org>
Thu, 26 Mar 2009 14:02:21 -0400
changeset 1 4be95183fbcd
permissions -rw-r--r--
Modifications to scoring and graphics production for the final version of code for the C-LAMP paper in GCB. 
     1 ; ***********************************************

     2 ; add zone plot to 21.lines.ncl

     3 ; ***********************************************

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

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

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

     7 ;************************************************

     8 begin

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

    10 ; read in data: observed

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

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

    13  fili  = "co2_globalView_98.nc"

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

    15  val   = g->CO2_SEAS  

    16  lon   = g->LON 

    17  lat   = g->LAT

    18  sta   = chartostring(g->STATION) 

    19  delete (g)

    20  

    21 ;print (sta(0))

    22 

    23  ncase = dimsizes(lat)

    24 ;print (ncase)

    25 

    26 ;**************************************************************

    27 ; get only the lowest level at each station 

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

    29  lat_tmp = lat

    30  lat_tmp@_FillValue = 1.e+36

    31  

    32  do n = 0,ncase-1

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

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

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

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

    37        end if

    38        delete (indexes)

    39     end if

    40  end do

    41 

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

    43 ;print (dimsizes(indexes))

    44 ;print (indexes)

    45  

    46  lat_ob = lat(indexes)

    47  lon_ob = lon(indexes)

    48  val_ob = val(indexes,:)

    49 ;printVarSummary (val_ob)

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

    51 

    52 ;************************************************

    53 ; read in model data

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

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

    56   fili2 = "b30.061m_401_425_MONS_climo.nc"

    57 

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

    59   x     = g->CO2

    60   xi    = g->lon

    61   yi    = g->lat

    62   xdim  = dimsizes(x)

    63   nlev  = xdim(1)

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

    65 ; printVarSummary (y)

    66   

    67 ; get the co2 at the lowest level

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

    69 

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

    71 ; Model_units [=] kgCO2 / kgDryAir

    72 ; 28.966 = molecular weight of dry air

    73 ; 44.       = molecular weight of CO2

    74 ; u mol = 1e-6 mol

    75 

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

    77   y      = y * factor

    78 

    79   y@_FillValue = 1.e36

    80   y@units      = "u mol/mol"

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

    82 ; printVarSummary (y)

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

    84 

    85 ; interpolate into observed station

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

    87 

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

    89   print (yi(0))

    90   yi(0) = -90.

    91 

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

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

    94 

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

    96 

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

    98 ; printVarSummary (val_model)

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

   100 

   101 ; remove annual mean

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

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

   104 

   105   plot_sta      = new((/2,12/),float)

   106   plot_sta@long_name = "Seasonal CO2"

   107 

   108   mon = ispan(1,12,1)

   109   mon@long_name = "month"

   110 

   111   plot_type = "ps"

   112   plot_type_new = "png"

   113 

   114   res                   = True                      ; plot mods desired

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

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

   117 

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

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

   120 ; vertical lines going through the markers.

   121 

   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)

   129 

   130 ; res@pmLegendParallelPosF   =  0.18

   131   res@pmLegendParallelPosF   =  0.23  ;(rightward)

   132 

   133 ; res@lgPerimOn             = False

   134   res@lgLabelFontHeightF     = 0.015

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

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

   137 

   138   nzone = 4

   139   do z = 0,nzone-1

   140 

   141   if (z .eq. 0) then 

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

   143      score_max = 5.0

   144      zone = "60N-90N"

   145 ;    index of stations in this zone

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

   147 ;    print (ind_z)

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

   149 ;    print (val_ob(ind_z,:))

   150 ;    print (val_model(ind_z,:))

   151   end if

   152 

   153   if (z .eq. 1) then 

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

   155      score_max = 5.0

   156      zone = "30N-60N"

   157 ;    index of stations in this zone

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

   159 ;    print (ind_z)

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

   161 ;    print (val_ob(ind_z,:))

   162 ;    print (val_model(ind_z,:))

   163   end if

   164 

   165   if (z .eq. 2) then 

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

   167      score_max = 5.0

   168      zone = "EQ-30N"

   169 ;    index of stations in this zone

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

   171 ;    print (ind_z)

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

   173 ;    print (val_ob(ind_z,:))

   174 ;    print (val_model(ind_z,:))

   175   end if

   176 

   177   if (z .eq. 3) then 

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

   179      score_max = 5.0

   180      zone = "90S-EQ"     

   181 ;    index of stations in this zone

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

   183 ;    print (ind_z)

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

   185 ;    print (val_ob(ind_z,:))

   186 ;    print (val_model(ind_z,:))

   187   end if

   188 

   189  npts = dimsizes(ind_z)

   190  print (npts)

   191 

   192  npts_str = ""

   193  npts_str = npts

   194 ;print (npts_str)

   195 

   196  plot_zon      = new((/2,12,npts/),float)

   197 

   198  do n=0,npts-1

   199     plot_zon(0,:,n) = (/val_ob(ind_z(n),:)/)

   200     plot_zon(1,:,n) = (/val_model(ind_z(n),:)/)

   201 

   202 ;   plot_sta(0,:) = (/val_ob(ind_z(n),:)/)

   203 ;   plot_sta(1,:) = (/val_model(ind_z(n),:)/)

   204     

   205 ;   title = sta(ind_z(n))+"("+lat(ind_z(n))+","+lon(ind_z(n))+")"

   206 ;   plot_name = sta(ind_z(n))

   207     

   208 ;   print (title)

   209 ;   print (plot_name)

   210 

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

   212 ;   res@tiMainString      = title                     ; add title

   213  

   214 ;   plot  = gsn_csm_xy (wks,mon,plot_sta,res)         ; create plot

   215 ;   frame(wks)

   216 

   217 ;   system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)

   218 ;   system("rm "+plot_name+"."+plot_type)

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

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

   221 ;   clear (wks)

   222  end do

   223     

   224     plot_sta = dim_avg(plot_zon)

   225   

   226     plot_name = "All_"+npts_str

   227     title = plot_name + " in "+ zone

   228 

   229     print (title)

   230     print (plot_name)

   231 

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

   233     res@tiMainString      = title                     ; add title

   234  

   235     plot  = gsn_csm_xy (wks,mon,plot_sta,res)         ; create plot

   236     frame(wks)

   237 

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

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

   240     system("rm "+plot_name+"-1."+plot_type_new)

   241     system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)

   242     clear (wks)

   243     delete (ind_z)

   244     delete (plot_zon)

   245 

   246  end do

   247 end
clamp