1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/fire/10.1x1_to_T42.ncl Mon Jan 26 22:08:20 2009 -0500
1.3 @@ -0,0 +1,141 @@
1.4 +;---------------------------------------------------------------
1.5 +; convert from 1x1 to T31/T42
1.6 +;
1.7 +; use model T31/T42 lat and lon
1.8 +; output model area, too.
1.9 +;
1.10 +;-------------------------------------------------------------
1.11 +load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
1.12 +load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
1.13 +load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
1.14 +undef("copyatt")
1.15 +load "/fis/cgd/cseg/people/jeff/ncl/remap_areaave.ncl"
1.16 +begin
1.17 +;************************************************
1.18 +; output data
1.19 +;************************************************
1.20 + diro = "/fis/cgd/cseg/people/jeff/fire_data/ob/GFEDv2_C/"
1.21 + filo = "Fire_C_1997-2006_monthly_T42.nc"
1.22 + c = addfile(diro+filo,"c")
1.23 + filedimdef(c,"year",-1,True)
1.24 +
1.25 +;******************************************************
1.26 +; input grid: 1deg x 1deg, N->S, 180W-180E
1.27 +;******************************************************
1.28 + diri = "/fis/cgd/cseg/people/jeff/fire_data/ob/GFEDv2_C/"
1.29 + fili = "Fire_C_1997-2006_monthly_1x1.nc"
1.30 + g = addfile (diri+fili,"r")
1.31 + bi = g->FIRE_C
1.32 + xi = g->lon
1.33 + yi = g->lat
1.34 +
1.35 + dsizes = dimsizes(bi)
1.36 + nyear = dsizes(0)
1.37 + nmonth = dsizes(1)
1.38 + nlatx = dsizes(2)
1.39 + nlonx = dsizes(3)
1.40 +
1.41 +;===================================================
1.42 +; change from 180W-180E, 90N-90S to 0-360E, 90S-90N
1.43 +;===================================================
1.44 +
1.45 + yi = (/ yi(::-1) /)
1.46 + bi = (/ bi(:,:,::-1,:) /)
1.47 +
1.48 + b2 = bi
1.49 + x2 = xi
1.50 +
1.51 + do i= 0,nlonx-1
1.52 + if (i .lt. 180) then
1.53 + p = i + 180
1.54 + xi(p) = x2(i) + 360.
1.55 + else
1.56 + p = i - 180
1.57 + xi(p) = x2(i)
1.58 + end if
1.59 + bi(:,:,:,p)= b2(:,:,:,i)
1.60 + end do
1.61 +
1.62 + bi&lat = yi
1.63 + bi&lon = xi
1.64 +
1.65 + print (xi)
1.66 + print (yi)
1.67 +
1.68 +;==================================================
1.69 +; edges
1.70 +;==================================================
1.71 +
1.72 + xi_edge = fspan(0.0, 360.0, nlonx+1)
1.73 + xi_edge@units = "degrees_east"
1.74 + yi_edge = fspan(-90.0, 90.0, nlatx+1)
1.75 + yi_edge@units = "degrees_north"
1.76 +
1.77 +;==================================================
1.78 +; area
1.79 +;==================================================
1.80 +
1.81 + deg2rad = 4.0*atan(1.0)/180.0
1.82 + earth_rad = 6.37122e6
1.83 + dx = earth_rad * deg2rad * (xi_edge(1)-xi_edge(0))
1.84 + dsin = sin(deg2rad*yi_edge(1:180)) - sin(deg2rad*yi_edge(0:179))
1.85 + area_in = new((/ nlatx, nlonx /), float)
1.86 + area_in = conform(area_in, earth_rad*dx*dsin,0)
1.87 + area_in@units = "m^2"
1.88 +
1.89 +;**************************************************************
1.90 +; output grid: (0-360), (90S-90N)
1.91 +;**************************************************************
1.92 +; interpolate into model T31/T42, using model lat and lon
1.93 +
1.94 + diri = "/fis/cgd/cseg/people/jeff/surface_data/"
1.95 + fili = "lnd_T42.nc"
1.96 + b = addfile(diri+fili,"r")
1.97 +
1.98 + yo = b->lat
1.99 + xo = b->lon
1.100 + nlat = dimsizes(yo)
1.101 + nlon = dimsizes(xo)
1.102 +
1.103 + xo_edge = fspan(0.0, 360.0, nlon+1) - 0.5*360.0/nlon
1.104 + xo_edge@units = "degrees_east"
1.105 +
1.106 + yo_edge = new(nlat+1, float)
1.107 + gau_info = gaus(nlat/2)
1.108 + yo_edge(0) = -90.0
1.109 + yo_edge(1:nlat-1) = doubletofloat(0.5*(gau_info(0:nlat-2,0)+gau_info(1:nlat-1,0)))
1.110 + yo_edge(nlat) = 90.0
1.111 + yo_edge@units = "degrees_north"
1.112 +
1.113 +; yo = 0.5*(yo_edge(0:nlat-1) + yo_edge(1:nlat))
1.114 +; yo@units = yo_edge@units
1.115 +; xo = 0.5*(xo_edge(0:nlon-1) + xo_edge(1:nlon))
1.116 +; xo@units = xo_edge@units
1.117 +
1.118 +;===================================================
1.119 +; output
1.120 +;===================================================
1.121 +
1.122 + w = new((/nyear,nmonth,nlat,nlon/),float)
1.123 +
1.124 + do m = 0,nyear-1
1.125 + do n = 0,nmonth-1
1.126 + w(m,n,:,:) = remap_areaave(xi_edge, yi_edge, bi(m,n,:,:), xo_edge, yo_edge)
1.127 + end do
1.128 + end do
1.129 +
1.130 + w!0 = "year"
1.131 + w!1 = "month"
1.132 + w!2 = "lat"
1.133 + w&lat = yo
1.134 + w!3 = "lon"
1.135 + w&lon = xo
1.136 + w@units = bi@units
1.137 + w@long_name = bi@long_name
1.138 + w@_FillValue = bi@_FillValue
1.139 +
1.140 + c->FIRE_C = w
1.141 + c->date = g->date
1.142 + c->area = b->area
1.143 +
1.144 +end