biomass/30.mask_to_T31-T42_1.9.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.
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
; read model grid
forrest@0
    12
;************************************************
forrest@0
    13
forrest@0
    14
 model_grid = "T31"
forrest@0
    15
;model_grid = "T42"
forrest@0
    16
;model_grid = "1.9"
forrest@0
    17
forrest@0
    18
 diri   = "/fis/cgd/cseg/people/jeff/clamp_data/surface_model/"
forrest@0
    19
 fili   = "lnd_"+ model_grid +".nc"
forrest@0
    20
 f      = addfile (diri+fili,"r")
forrest@0
    21
forrest@0
    22
 lon    = f->lon     
forrest@0
    23
 lat    = f->lat
forrest@0
    24
 nlon   = dimsizes(lon)
forrest@0
    25
 nlat   = dimsizes(lat)      
forrest@0
    26
forrest@0
    27
;************************************************
forrest@0
    28
; output data
forrest@0
    29
;************************************************
forrest@0
    30
  diro  = "/fis/cgd/cseg/people/jeff/clamp_data/observed/biomass/"
forrest@0
    31
  filo  = "amazon_mask_"+model_grid + ".nc"
forrest@0
    32
  c = addfile(diro+filo,"c")
forrest@0
    33
forrest@0
    34
;************************************************
forrest@0
    35
; read 1x1 amazon_mask data
forrest@0
    36
;************************************************
forrest@0
    37
  diri  = "/fis/cgd/cseg/people/jeff/clamp_data/observed/biomass/"
forrest@0
    38
  fili  = "amazon_mask_1x1.nc"
forrest@0
    39
  g = addfile(diri+fili,"r")
forrest@0
    40
forrest@0
    41
  bi    = g->mask_amazon   
forrest@0
    42
  xi    = g->lon 
forrest@0
    43
  yi    = g->lat
forrest@0
    44
forrest@0
    45
;************************************************
forrest@0
    46
; change from -180-180, 90N-90S into 0-360E, 90S-90N
forrest@0
    47
; Observed NPP*scale_factor
forrest@0
    48
;************************************************
forrest@0
    49
 
forrest@0
    50
 yi    = (/ yi(::-1) /)
forrest@0
    51
 bi    =  (/ bi(::-1,:) /)
forrest@0
    52
 b2    = bi
forrest@0
    53
 x2    = xi   
forrest@0
    54
 
forrest@0
    55
 nx = dimsizes(xi)
forrest@0
    56
 do i= 0,nx-1
forrest@0
    57
    if (i .lt. 180) then
forrest@0
    58
       p = i + 180
forrest@0
    59
       xi(p) = x2(i) + 360.      
forrest@0
    60
    else
forrest@0
    61
       p = i - 180
forrest@0
    62
       xi(p) = x2(i)
forrest@0
    63
    end if
forrest@0
    64
    bi(:,p)= b2(:,i) 
forrest@0
    65
 end do
forrest@0
    66
forrest@0
    67
 bi&lat =  yi
forrest@0
    68
 bi&lon =  xi
forrest@0
    69
forrest@0
    70
;print (xi)
forrest@0
    71
;print (yi)
forrest@0
    72
forrest@0
    73
;(A)
forrest@0
    74
;bo = linint2_Wrap(xi,yi,bi,True,xo,yo,0)
forrest@0
    75
forrest@0
    76
;(B)
forrest@0
    77
forrest@0
    78
  rad   = 4.*atan(1.)/180.
forrest@0
    79
  clat  = lat
forrest@0
    80
  clat  = cos(lat*rad)
forrest@0
    81
forrest@0
    82
  clat@long_name = "cos(latitude)"
forrest@0
    83
  delete(clat@units)
forrest@0
    84
; printVarSummary(clat)
forrest@0
    85
forrest@0
    86
 bo = new ((/nlat,nlon/),"float")
forrest@0
    87
forrest@0
    88
    do j=0,nlat-1
forrest@0
    89
       if (j.eq.0 .or. j.eq.nlat-1) then
forrest@0
    90
          if (j.eq.0) then
forrest@0
    91
             LATS = -90.          
forrest@0
    92
             LATN = lat(j)+0.5*(lat(j+1)-lat(j))
forrest@0
    93
          end if
forrest@0
    94
          if (j.eq.nlat-1) then
forrest@0
    95
             LATS = lat(j)-0.5*(lat(j)-lat(j-1))
forrest@0
    96
             LATN = 90.                  
forrest@0
    97
          end if
forrest@0
    98
       else
forrest@0
    99
          LATS = lat(j)-0.5*(lat(j)-lat(j-1))
forrest@0
   100
          LATN = lat(j)+0.5*(lat(j+1)-lat(j))
forrest@0
   101
       end if
forrest@0
   102
 
forrest@0
   103
;      CLAT = clat({LATS:LATN})      ; do once for *slight* efficiency
forrest@0
   104
;      TEMP = bi({LATS:LATN},:)      ; 2D [lat,lon]
forrest@0
   105
 
forrest@0
   106
      do i=0,nlon-1
forrest@0
   107
       if (i.eq.0 .or. i.eq.nlon-1) then
forrest@0
   108
          if (i.eq.0) then
forrest@0
   109
             LONL = 0.          
forrest@0
   110
             LONR = lon(i)+0.5*(lon(i+1)-lon(i))
forrest@0
   111
          end if
forrest@0
   112
          if (i.eq.nlon-1) then
forrest@0
   113
             LONL = lon(i)-0.5*(lon(i)-lon(i-1))
forrest@0
   114
             LONR = 360.                 
forrest@0
   115
          end if
forrest@0
   116
       else
forrest@0
   117
          LONL = lon(i)-0.5*(lon(i)-lon(i-1))
forrest@0
   118
          LONR = lon(i)+0.5*(lon(i+1)-lon(i))
forrest@0
   119
       end if
forrest@0
   120
forrest@0
   121
;print (LATS)
forrest@0
   122
;print (LATN)
forrest@0
   123
;print (LONL)
forrest@0
   124
;print (LONR)
forrest@0
   125
forrest@0
   126
         bo(j,i) = avg(bi({LATS:LATN},{LONL:LONR}))  
forrest@0
   127
;        bo(j,i) = wgt_areaave(TEMP(:,{LONL:LONR}), CLAT, 1.0, 0)
forrest@0
   128
forrest@0
   129
      end do
forrest@0
   130
 
forrest@0
   131
;     delete(CLAT)
forrest@0
   132
;     delete(TEMP) 
forrest@0
   133
    end do
forrest@0
   134
forrest@0
   135
; bo = where(bo.ge.0.5, 1.0,0.)
forrest@0
   136
forrest@0
   137
  bo!0   = "lat"
forrest@0
   138
  bo!1   = "lon"
forrest@0
   139
  bo&lat =  lat
forrest@0
   140
  bo&lon =  lon
forrest@0
   141
; bo@units      = bi@units
forrest@0
   142
  bo@long_name  = "amazon mask"
forrest@0
   143
; bo@_FillValue = bi@_FillValue
forrest@0
   144
  bo@_FillValue = 1.e+36
forrest@0
   145
forrest@0
   146
  c->mask_amazon  = bo
forrest@0
   147
end