biomass/03.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
; interpolate into model grids (T31)
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
forrest@0
    11
;************************************************
forrest@0
    12
; read in observed data
forrest@0
    13
;************************************************
forrest@0
    14
  diri  = "/fis/cgd/cseg/people/jeff/clamp_data/biomass/ob/"
forrest@0
    15
  fili  = "amazon_biomass.nc"
forrest@0
    16
  g     = addfile (diri+fili,"r")
forrest@0
    17
  bi    = g->BIOMASS   
forrest@0
    18
  xi    = g->lon 
forrest@0
    19
  yi    = g->lat
forrest@0
    20
forrest@0
    21
forrest@0
    22
;************************************************
forrest@0
    23
; read in model data
forrest@0
    24
;************************************************
forrest@0
    25
 diri2  = "/fis/cgd/cseg/people/jeff/clamp_data/model/"
forrest@0
    26
 fili2  = "b30.061n_1995-2004_ANN_climo_lnd.nc"
forrest@0
    27
 f      = addfile (diri2+fili2,"r")
forrest@0
    28
forrest@0
    29
 lon    = f->lon     
forrest@0
    30
 lat    = f->lat
forrest@0
    31
 nlon   = dimsizes(lon)
forrest@0
    32
 nlat   = dimsizes(lat)      
forrest@0
    33
forrest@0
    34
; print (xi)
forrest@0
    35
; print (yi)
forrest@0
    36
forrest@0
    37
 ind_lonL = min(ind(lon .ge. xi(0)))
forrest@0
    38
 ind_lonR = max(ind(lon .le. xi(dimsizes(xi)-1)))
forrest@0
    39
forrest@0
    40
print (ind_lonL)
forrest@0
    41
print (xi(0))
forrest@0
    42
print (lon(ind_lonL))
forrest@0
    43
forrest@0
    44
print (ind_lonR)
forrest@0
    45
print (xi(dimsizes(xi)-1))
forrest@0
    46
print (lon(ind_lonR))
forrest@0
    47
forrest@0
    48
 ind_latS = min(ind(lat .ge. yi(0)))
forrest@0
    49
 ind_latN = max(ind(lat .le. yi(dimsizes(yi)-1)))
forrest@0
    50
forrest@0
    51
print (ind_latS)
forrest@0
    52
print (yi(0))
forrest@0
    53
print (lat(ind_latS))
forrest@0
    54
forrest@0
    55
print (ind_latN)
forrest@0
    56
print (yi(dimsizes(yi)-1))
forrest@0
    57
print (lat(ind_latN))
forrest@0
    58
exit
forrest@0
    59
 bo = new((/nlat,nlon/),float)
forrest@0
    60
    do j=0,nlat-1
forrest@0
    61
       if (j.eq.0 .or. j.eq.nlat-1) then
forrest@0
    62
          if (j.eq.0) then
forrest@0
    63
             LATS = -90.          
forrest@0
    64
             LATN = lat(j)+0.5*(lat(j+1)-lat(j))
forrest@0
    65
          end if
forrest@0
    66
          if (j.eq.nlat-1) then
forrest@0
    67
             LATS = lat(j)-0.5*(lat(j)-lat(j-1))
forrest@0
    68
             LATN = 90.                  
forrest@0
    69
          end if
forrest@0
    70
       else
forrest@0
    71
          LATS = lat(j)-0.5*(lat(j)-lat(j-1))
forrest@0
    72
          LATN = lat(j)+0.5*(lat(j+1)-lat(j))
forrest@0
    73
       end if
forrest@0
    74
 
forrest@0
    75
;      CLAT = clat({LATS:LATN})      ; do once for *slight* efficiency
forrest@0
    76
;      TEMP = bi({LATS:LATN},:)      ; 2D [lat,lon]
forrest@0
    77
 
forrest@0
    78
      do i=0,nlon-1
forrest@0
    79
       if (i.eq.0 .or. i.eq.nlon-1) then
forrest@0
    80
          if (i.eq.0) then
forrest@0
    81
             LONL = 0.          
forrest@0
    82
             LONR = lon(i)+0.5*(lon(i+1)-lon(i))
forrest@0
    83
          end if
forrest@0
    84
          if (i.eq.nlon-1) then
forrest@0
    85
             LONL = lon(i)-0.5*(lon(i)-lon(i-1))
forrest@0
    86
             LONR = 360.                 
forrest@0
    87
          end if
forrest@0
    88
       else
forrest@0
    89
          LONL = lon(i)-0.5*(lon(i)-lon(i-1))
forrest@0
    90
          LONR = lon(i)+0.5*(lon(i+1)-lon(i))
forrest@0
    91
       end if
forrest@0
    92
forrest@0
    93
;print (LATS)
forrest@0
    94
;print (LATN)
forrest@0
    95
;print (LONL)
forrest@0
    96
;print (LONR)
forrest@0
    97
         bo(j,i) = avg(bi({LATS:LATN},{LONL:LONR}))  
forrest@0
    98
;        bo(j,i) = wgt_areaave(TEMP(:,{LONL:LONR}), CLAT, 1.0, 0)
forrest@0
    99
      end do
forrest@0
   100
 
forrest@0
   101
;     delete(CLAT)
forrest@0
   102
;     delete(TEMP) 
forrest@0
   103
    end do
forrest@0
   104
forrest@0
   105
  bo!0   = "lat"
forrest@0
   106
  bo!1   = "lon"
forrest@0
   107
  bo&lat =  lat
forrest@0
   108
  bo&lon =  lon
forrest@0
   109
  bo@units      = bi@units
forrest@0
   110
  bo@long_name  = bi@long_name
forrest@0
   111
; bo@_FillValue = bi@_FillValue
forrest@0
   112
  bo@_FillValue = 1.e+36
forrest@0
   113
forrest@0
   114
  c->NPP  = bo
forrest@0
   115
end