clamp: fire/10.1x1_to_T31.ncl@4be95183fbcd

     1 ;---------------------------------------------------------------

     2 ; convert from 1x1 to T31/T42

3 ;

     4 ; use model T31/T42 lat and lon

     5 ; output model area, too.

6 ;

     7 ;-------------------------------------------------------------

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

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

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

    11 undef("copyatt")

    12 load "/fis/cgd/cseg/people/jeff/ncl/remap_areaave.ncl"

    13 begin

    14 ;************************************************

    15 ; output data

    16 ;************************************************

    17   diro  = "/fis/cgd/cseg/people/jeff/fire_data/ob/GFEDv2_C/"

    18   filo  = "Fire_C_1997-2006_monthly_T31.nc"

    19   c = addfile(diro+filo,"c")

    20   filedimdef(c,"year",-1,True)

    22 ;******************************************************

    23 ; input grid: 1deg x 1deg, N->S, 180W-180E

    24 ;******************************************************

    25   diri  = "/fis/cgd/cseg/people/jeff/fire_data/ob/GFEDv2_C/"

    26   fili  = "Fire_C_1997-2006_monthly_1x1.nc"

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

    28   bi    = g->FIRE_C

    29   xi    = g->lon

    30   yi    = g->lat

    32   dsizes = dimsizes(bi)

    33   nyear  = dsizes(0)

    34   nmonth = dsizes(1)

    35   nlatx  = dsizes(2)

    36   nlonx  = dsizes(3)

    38 ;===================================================

    39 ; change from 180W-180E, 90N-90S to  0-360E, 90S-90N

    40 ;===================================================

    42  yi    = (/ yi(::-1) /)

    43  bi    = (/ bi(:,:,::-1,:) /)

    45  b2    = bi

    46  x2    = xi

    48  do i= 0,nlonx-1

    49     if (i .lt. 180) then

    50        p = i + 180

    51        xi(p) = x2(i) + 360.

    52     else

    53        p = i - 180

    54        xi(p) = x2(i)

    55     end if

    56     bi(:,:,:,p)= b2(:,:,:,i)

    57  end do

    59  bi&lat =  yi

    60  bi&lon =  xi

    62  print (xi)

    63  print (yi)

    65 ;==================================================

    66 ; edges

    67 ;==================================================

    69   xi_edge = fspan(0.0, 360.0, nlonx+1)

    70   xi_edge@units = "degrees_east"

    71   yi_edge = fspan(-90.0, 90.0, nlatx+1)

    72   yi_edge@units = "degrees_north"

    74 ;==================================================

    75 ; area

    76 ;==================================================

    78   deg2rad = 4.0*atan(1.0)/180.0

    79   earth_rad = 6.37122e6

    80   dx = earth_rad * deg2rad * (xi_edge(1)-xi_edge(0))

    81   dsin = sin(deg2rad*yi_edge(1:180)) - sin(deg2rad*yi_edge(0:179))

    82   area_in = new((/ nlatx, nlonx /), float)

    83   area_in = conform(area_in, earth_rad*dx*dsin,0)

    84   area_in@units = "m^2"

    86 ;**************************************************************

    87 ; output grid: (0-360), (90S-90N)

    88 ;**************************************************************

    89 ; interpolate into model T31/T42, using model data

    90   diri  = "/fis/cgd/cseg/people/jeff/fire_data/model/"

    91   fili  = "b30.061n_1980-1990_ANN_climo.nc"

    92   b = addfile(diri+fili,"r")

    94   yo   = b->lat

    95   xo   = b->lon

    96   nlat = dimsizes(yo)

    97   nlon = dimsizes(xo)

    99   xo_edge = fspan(0.0, 360.0, nlon+1) - 0.5*360.0/nlon

   100   xo_edge@units = "degrees_east"

   102   yo_edge = new(nlat+1, float)

   103   gau_info = gaus(nlat/2)

   104   yo_edge(0) = -90.0

   105   yo_edge(1:nlat-1) = doubletofloat(0.5*(gau_info(0:nlat-2,0)+gau_info(1:nlat-1,0)))

   106   yo_edge(nlat) = 90.0

   107   yo_edge@units = "degrees_north"

   109 ; yo = 0.5*(yo_edge(0:nlat-1) + yo_edge(1:nlat))

   110 ; yo@units = yo_edge@units

   111 ; xo = 0.5*(xo_edge(0:nlon-1) + xo_edge(1:nlon))

   112 ; xo@units = xo_edge@units

   114 ;===================================================

   115 ; output

   116 ;===================================================

   118   w  = new((/nyear,nmonth,nlat,nlon/),float)

   120   do m = 0,nyear-1

   121   do n = 0,nmonth-1

   122      w(m,n,:,:) = remap_areaave(xi_edge, yi_edge, bi(m,n,:,:), xo_edge, yo_edge)

   123   end do

   124   end do

   126   w!0   = "year"

   127   w!1   = "month"

   128   w!2   = "lat"

   129   w&lat = yo

   130   w!3   = "lon"

   131   w&lon = xo

   132   w@units      = bi@units

   133   w@long_name  = bi@long_name

   134   w@_FillValue = bi@_FillValue

   136   c->FIRE_C = w

   137   c->date   = g->date

   138   c->area   = b->area

   140 end

author	Forrest Hoffman <forrest@climatemodeling.org>
	Thu, 26 Mar 2009 14:02:21 -0400
changeset 1	4be95183fbcd
permissions	-rw-r--r--