forrest@0: ;---------------------------------------------------------------
forrest@0: ; convert from 1x1 to T31/T42 
forrest@0: ;
forrest@0: ; use model T31/T42 lat and lon
forrest@0: ; output model area, too.
forrest@0: ;
forrest@0: ;-------------------------------------------------------------
forrest@0: load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
forrest@0: load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
forrest@0: load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
forrest@0: undef("copyatt")
forrest@0: load "/fis/cgd/cseg/people/jeff/ncl/remap_areaave.ncl"
forrest@0: begin
forrest@0: ;************************************************
forrest@0: ; output data
forrest@0: ;************************************************
forrest@0:   diro  = "/fis/cgd/cseg/people/jeff/fire_data/ob/GFEDv2_C/"
forrest@0:   filo  = "Fire_C_1997-2006_monthly_T42.nc"
forrest@0:   c = addfile(diro+filo,"c")
forrest@0:   filedimdef(c,"year",-1,True)
forrest@0: 
forrest@0: ;******************************************************
forrest@0: ; input grid: 1deg x 1deg, N->S, 180W-180E
forrest@0: ;******************************************************
forrest@0:   diri  = "/fis/cgd/cseg/people/jeff/fire_data/ob/GFEDv2_C/"
forrest@0:   fili  = "Fire_C_1997-2006_monthly_1x1.nc"
forrest@0:   g     = addfile (diri+fili,"r")
forrest@0:   bi    = g->FIRE_C   
forrest@0:   xi    = g->lon 
forrest@0:   yi    = g->lat
forrest@0: 
forrest@0:   dsizes = dimsizes(bi)
forrest@0:   nyear  = dsizes(0)
forrest@0:   nmonth = dsizes(1)
forrest@0:   nlatx  = dsizes(2)
forrest@0:   nlonx  = dsizes(3) 
forrest@0: 
forrest@0: ;===================================================
forrest@0: ; change from 180W-180E, 90N-90S to  0-360E, 90S-90N
forrest@0: ;===================================================
forrest@0:  
forrest@0:  yi    = (/ yi(::-1) /)
forrest@0:  bi    = (/ bi(:,:,::-1,:) /)
forrest@0: 
forrest@0:  b2    = bi
forrest@0:  x2    = xi   
forrest@0:  
forrest@0:  do i= 0,nlonx-1
forrest@0:     if (i .lt. 180) then
forrest@0:        p = i + 180
forrest@0:        xi(p) = x2(i) + 360.      
forrest@0:     else
forrest@0:        p = i - 180
forrest@0:        xi(p) = x2(i)
forrest@0:     end if
forrest@0:     bi(:,:,:,p)= b2(:,:,:,i) 
forrest@0:  end do
forrest@0: 
forrest@0:  bi&lat =  yi
forrest@0:  bi&lon =  xi
forrest@0: 
forrest@0:  print (xi)
forrest@0:  print (yi)
forrest@0: 
forrest@0: ;==================================================
forrest@0: ; edges 
forrest@0: ;==================================================
forrest@0:  
forrest@0:   xi_edge = fspan(0.0, 360.0, nlonx+1)
forrest@0:   xi_edge@units = "degrees_east"
forrest@0:   yi_edge = fspan(-90.0, 90.0, nlatx+1)
forrest@0:   yi_edge@units = "degrees_north"
forrest@0: 
forrest@0: ;==================================================
forrest@0: ; area
forrest@0: ;==================================================
forrest@0: 
forrest@0:   deg2rad = 4.0*atan(1.0)/180.0
forrest@0:   earth_rad = 6.37122e6
forrest@0:   dx = earth_rad * deg2rad * (xi_edge(1)-xi_edge(0))
forrest@0:   dsin = sin(deg2rad*yi_edge(1:180)) - sin(deg2rad*yi_edge(0:179))
forrest@0:   area_in = new((/ nlatx, nlonx /), float)
forrest@0:   area_in = conform(area_in, earth_rad*dx*dsin,0)
forrest@0:   area_in@units = "m^2"
forrest@0: 
forrest@0: ;**************************************************************
forrest@0: ; output grid: (0-360), (90S-90N)
forrest@0: ;**************************************************************
forrest@0: ; interpolate into model T31/T42, using model lat and lon
forrest@0: 
forrest@0:   diri  = "/fis/cgd/cseg/people/jeff/surface_data/"
forrest@0:   fili  = "lnd_T42.nc"   
forrest@0:   b = addfile(diri+fili,"r")
forrest@0:  
forrest@0:   yo   = b->lat
forrest@0:   xo   = b->lon
forrest@0:   nlat = dimsizes(yo)
forrest@0:   nlon = dimsizes(xo)       
forrest@0: 
forrest@0:   xo_edge = fspan(0.0, 360.0, nlon+1) - 0.5*360.0/nlon
forrest@0:   xo_edge@units = "degrees_east"
forrest@0: 
forrest@0:   yo_edge = new(nlat+1, float)
forrest@0:   gau_info = gaus(nlat/2)
forrest@0:   yo_edge(0) = -90.0
forrest@0:   yo_edge(1:nlat-1) = doubletofloat(0.5*(gau_info(0:nlat-2,0)+gau_info(1:nlat-1,0)))
forrest@0:   yo_edge(nlat) = 90.0
forrest@0:   yo_edge@units = "degrees_north" 
forrest@0: 
forrest@0: ; yo = 0.5*(yo_edge(0:nlat-1) + yo_edge(1:nlat))
forrest@0: ; yo@units = yo_edge@units
forrest@0: ; xo = 0.5*(xo_edge(0:nlon-1) + xo_edge(1:nlon))
forrest@0: ; xo@units = xo_edge@units
forrest@0:   
forrest@0: ;===================================================
forrest@0: ; output
forrest@0: ;===================================================
forrest@0: 
forrest@0:   w  = new((/nyear,nmonth,nlat,nlon/),float)
forrest@0: 
forrest@0:   do m = 0,nyear-1
forrest@0:   do n = 0,nmonth-1
forrest@0:      w(m,n,:,:) = remap_areaave(xi_edge, yi_edge, bi(m,n,:,:), xo_edge, yo_edge)
forrest@0:   end do
forrest@0:   end do
forrest@0: 
forrest@0:   w!0   = "year"
forrest@0:   w!1   = "month"               
forrest@0:   w!2   = "lat" 
forrest@0:   w&lat = yo            
forrest@0:   w!3   = "lon"
forrest@0:   w&lon = xo  
forrest@0:   w@units      = bi@units
forrest@0:   w@long_name  = bi@long_name
forrest@0:   w@_FillValue = bi@_FillValue   
forrest@0:  
forrest@0:   c->FIRE_C = w
forrest@0:   c->date   = g->date
forrest@0:   c->area   = b->area
forrest@0:  
forrest@0: end