co2/12.test.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
forrest@0
     1
; ***********************************************
forrest@0
     2
; xy_4.ncl
forrest@0
     3
; ***********************************************
forrest@0
     4
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
forrest@0
     5
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
forrest@0
     6
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
forrest@0
     7
;************************************************
forrest@0
     8
begin
forrest@0
     9
;************************************************
forrest@0
    10
; read in data: observed
forrest@0
    11
;************************************************
forrest@0
    12
 diri  = "/fis/cgd/cseg/people/jeff/clamp_data/co2/"
forrest@0
    13
 fili  = "co2_globalView_222.nc"
forrest@0
    14
 g     = addfile (diri+fili,"r")
forrest@0
    15
 val   = g->CO2_SEAS  
forrest@0
    16
 lon   = g->LON 
forrest@0
    17
 lat   = g->LAT 
forrest@0
    18
 delete (g)
forrest@0
    19
 
forrest@0
    20
 ncase = dimsizes(lat)
forrest@0
    21
;print (ncase)
forrest@0
    22
forrest@0
    23
;**************************************************************
forrest@0
    24
; get only the lowest level at each station 
forrest@0
    25
;**************************************************************
forrest@0
    26
 lat_tmp = lat
forrest@0
    27
 lat_tmp@_FillValue = 1.e+36
forrest@0
    28
 
forrest@0
    29
 do n = 0,ncase-1
forrest@0
    30
    if (.not. ismissing(lat_tmp(n))) then 
forrest@0
    31
       indexes = ind(lat(n) .eq. lat .and. lon(n) .eq. lon)
forrest@0
    32
       if (dimsizes(indexes) .gt. 1) then
forrest@0
    33
          lat_tmp(indexes(1:)) = lat_tmp@_FillValue
forrest@0
    34
       end if
forrest@0
    35
       delete (indexes)
forrest@0
    36
    end if
forrest@0
    37
 end do
forrest@0
    38
forrest@0
    39
 indexes = ind(.not. ismissing(lat_tmp))
forrest@0
    40
;print (indexes)
forrest@0
    41
 
forrest@0
    42
 lat_ob = lat(indexes)
forrest@0
    43
 lon_ob = lon(indexes)
forrest@0
    44
 val_ob = val(indexes,:)
forrest@0
    45
;printVarSummary (val_ob)
forrest@0
    46
;print (lat_ob +"/"+lon_ob)
forrest@0
    47
forrest@0
    48
;************************************************
forrest@0
    49
; read in model data
forrest@0
    50
;************************************************
forrest@0
    51
  diri2 = "/fis/cgd/cseg/people/jeff/clamp_data/model/"
forrest@0
    52
  fili2 = "b30.061m_401_425_MONS_climo.nc"
forrest@0
    53
forrest@0
    54
  g     = addfile(diri2+fili2,"r")
forrest@0
    55
  x     = g->CO2
forrest@0
    56
  xi    = g->lon
forrest@0
    57
  yi    = g->lat
forrest@0
    58
  xdim  = dimsizes(x)
forrest@0
    59
  nlev  = xdim(1)
forrest@0
    60
  y     = x(:,0,:,:)
forrest@0
    61
; printVarSummary (y)
forrest@0
    62
  
forrest@0
    63
; get the co2 at the lowest level
forrest@0
    64
  y     = x(:,nlev-1,:,:)
forrest@0
    65
forrest@0
    66
; change to unit of observed (u mol/mol)
forrest@0
    67
; Model_units [=] kgCO2 / kgDryAir
forrest@0
    68
; 28.966 = molecular weight of dry air
forrest@0
    69
; 44.       = molecular weight of CO2
forrest@0
    70
; u mol = 1e-6 mol
forrest@0
    71
forrest@0
    72
  factor = (28.966/44.) * 1e6
forrest@0
    73
  y      = y * factor
forrest@0
    74
forrest@0
    75
  y@_FillValue = 1.e36
forrest@0
    76
  y@units      = "u mol/mol"
forrest@0
    77
; y = where(y0 .lt. 287.,y@_FillValue,y)
forrest@0
    78
; printVarSummary (y)
forrest@0
    79
; print (min(y)+"/"+max(y))
forrest@0
    80
forrest@0
    81
; interpolate into observed station
forrest@0
    82
forrest@0
    83
  yo = linint2_points_Wrap(xi,yi,y,True,lon_ob,lat_ob,0)
forrest@0
    84
forrest@0
    85
  val_model = val_ob
forrest@0
    86
  val_model = yo(pts|:,time|:)
forrest@0
    87
  printVarSummary (val_model)
forrest@0
    88
  print (min(val_model)+"/"+max(val_model))
forrest@0
    89
forrest@0
    90
; remove annual mean
forrest@0
    91
  val_model = val_model - conform(val_model,dim_avg(val_model),0)
forrest@0
    92
  print (min(val_model)+"/"+max(val_model))
forrest@0
    93
exit
forrest@0
    94
;**************************************************************
forrest@0
    95
; get stations of 60N-90N
forrest@0
    96
;**************************************************************
forrest@0
    97
forrest@0
    98
 ind_1 = ind(lat_new .gt. 60.)
forrest@0
    99
 print (ind_1)
forrest@0
   100
 print (lat_new(ind_1)+"/"+lon_new(ind_1))
forrest@0
   101
 print (val_new(ind_1(0),:))
forrest@0
   102
exit
forrest@0
   103
forrest@0
   104
end