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/24.check.ncl
author Forrest Hoffman <forrest@climatemodeling.org>
Mon, 26 Jan 2009 22:08:20 -0500
changeset 0 0c6405ab2ff4
permissions -rw-r--r--
Initial commit of C-LAMP Diagnostics from Jeff Lee 
     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)

    22  

    23 ;print (sta(0))

    24 

    25  ncase = dimsizes(lat)

    26 ;print (ncase)

    27 

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

    29 ; get only the lowest level at each station 

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

    31  lat_tmp = lat

    32  lat_tmp@_FillValue = 1.e+36

    33  

    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

    43 

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

    45 ;print (dimsizes(indexes))

    46 ;print (indexes)

    47  

    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)

    53 

    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)

    66 

    67   xdim  = dimsizes(x1)

    68   nlev  = xdim(1)

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

    70   y1     = x1

    71 ; printVarSummary (y1)

    72   

    73 ; get the co2 at the lowest level

    74 ; y1     = x1(:,nlev-1,:,:)

    75   delete (x1)

    76 ;---------------------------------------------

    77 ; g     = addfile(diri2+fili3,"r")

    78 ; x2    = g->CO2

    79 ; delete (g)

    80 ; y2     = x2(:,0,:,:)

    81 ; y2     = x2(:,nlev-1,:,:)

    82 ; delete (x2)

    83 ;---------------------------------------------

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

    85 ; Model_units [=] kgCO2 / kgDryAir

    86 ; 28.966 = molecular weight of dry air

    87 ; 44.       = molecular weight of CO2

    88 ; u mol = 1e-6 mol

    89 

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

    91 ;---------------------------------------------

    92   y1 = y1 * factor

    93   y1@_FillValue = 1.e36

    94   y1@units      = "u mol/mol"

    95 ; y1 = where(y0 .lt. 287.,y1@_FillValue,y1)

    96 ; printVarSummary (y1)

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

    98 ;---------------------------------------------

    99 ; y2 = y2 * factor

   100 ; y2@_FillValue = 1.e36

   101 ; y2@units      = "u mol/mol"

   102 ;---------------------------------------------

   103 ; interpolate into observed station

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

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

   106 

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

   108 ; print (yi(0))

   109   yi(0) = -90.

   110 

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

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

   113 ;----------------------------------------------------------------

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

   115   printVarSummary (yo)

   116 ; yo:[time | 12] x [lev | 26] x [pts | 98]

   117 

   118   val_model1 = yo(pts|:,lev|:,time|:)

   119   delete (yo)

   120   val_model1_0 = val_model1

   121 ; printVarSummary (val_model1)

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

   123 

   124 ; remove annual mean

   125   val_model1 = val_model1 - conform(val_model1,dim_avg(val_model1),(/0,1/))

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

   127 

   128 ;-----------------------------------------------------------------

   129 ;    index of station Barrow, Alaska (71.32,-156.60)

   130      ind_z = ind(lat_ob .eq. 71.32)

   131      print (ind_z)

   132      print (lat_ob(ind_z)+"/"+lon_ob(ind_z))

   133      print ("observation at Barrow, Alaska (71.32,-156.60)")

   134      print (val_ob(ind_z,:))

   135      print ("model top atm  at Barrow, Alaska (71.32,-156.60)")

   136      print (val_model1_0(ind_z,0,:))

   137      print ("model surface at Barrow, Alaska (71.32,-156.60)")

   138      print (val_model1_0(ind_z,nlev-1,:))

   139      print ("model top atm  at Barrow, Alaska (71.32,-156.60)")

   140      print (val_model1(ind_z,0,:))

   141      print ("model surface at Barrow, Alaska (71.32,-156.60)")

   142      print (val_model1(ind_z,nlev-1,:))

   143 end
clamp