fire/10.1x1_to_1.9.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
; convert from 1x1 to T31/T42 
forrest@0
     3
;
forrest@0
     4
; use model T31/T42 lat and lon
forrest@0
     5
; output model area, too.
forrest@0
     6
;
forrest@0
     7
;-------------------------------------------------------------
forrest@0
     8
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
forrest@0
     9
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
forrest@0
    10
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
forrest@0
    11
undef("copyatt")
forrest@0
    12
load "/fis/cgd/cseg/people/jeff/ncl/remap_areaave.ncl"
forrest@0
    13
begin
forrest@0
    14
;************************************************
forrest@0
    15
; output data
forrest@0
    16
;************************************************
forrest@0
    17
  diro  = "/fis/cgd/cseg/people/jeff/fire_data/ob/GFEDv2_C/"
forrest@0
    18
  filo  = "Fire_C_1997-2006_monthly_1.9.nc"
forrest@0
    19
  c = addfile(diro+filo,"c")
forrest@0
    20
  filedimdef(c,"year",-1,True)
forrest@0
    21
forrest@0
    22
;******************************************************
forrest@0
    23
; input grid: 1deg x 1deg, N->S, 180W-180E
forrest@0
    24
;******************************************************
forrest@0
    25
  diri  = "/fis/cgd/cseg/people/jeff/fire_data/ob/GFEDv2_C/"
forrest@0
    26
  fili  = "Fire_C_1997-2006_monthly_1x1.nc"
forrest@0
    27
  g     = addfile (diri+fili,"r")
forrest@0
    28
  bi    = g->FIRE_C   
forrest@0
    29
  xi    = g->lon 
forrest@0
    30
  yi    = g->lat
forrest@0
    31
forrest@0
    32
  dsizes = dimsizes(bi)
forrest@0
    33
  nyear  = dsizes(0)
forrest@0
    34
  nmonth = dsizes(1)
forrest@0
    35
  nlatx  = dsizes(2)
forrest@0
    36
  nlonx  = dsizes(3) 
forrest@0
    37
forrest@0
    38
;===================================================
forrest@0
    39
; change from 180W-180E, 90N-90S to  0-360E, 90S-90N
forrest@0
    40
;===================================================
forrest@0
    41
 
forrest@0
    42
 yi    = (/ yi(::-1) /)
forrest@0
    43
 bi    = (/ bi(:,:,::-1,:) /)
forrest@0
    44
forrest@0
    45
 b2    = bi
forrest@0
    46
 x2    = xi   
forrest@0
    47
 
forrest@0
    48
 do i= 0,nlonx-1
forrest@0
    49
    if (i .lt. 180) then
forrest@0
    50
       p = i + 180
forrest@0
    51
       xi(p) = x2(i) + 360.      
forrest@0
    52
    else
forrest@0
    53
       p = i - 180
forrest@0
    54
       xi(p) = x2(i)
forrest@0
    55
    end if
forrest@0
    56
    bi(:,:,:,p)= b2(:,:,:,i) 
forrest@0
    57
 end do
forrest@0
    58
forrest@0
    59
 bi&lat =  yi
forrest@0
    60
 bi&lon =  xi
forrest@0
    61
forrest@0
    62
 print (xi)
forrest@0
    63
 print (yi)
forrest@0
    64
forrest@0
    65
;==================================================
forrest@0
    66
; edges 
forrest@0
    67
;==================================================
forrest@0
    68
 
forrest@0
    69
  xi_edge = fspan(0.0, 360.0, nlonx+1)
forrest@0
    70
  xi_edge@units = "degrees_east"
forrest@0
    71
  yi_edge = fspan(-90.0, 90.0, nlatx+1)
forrest@0
    72
  yi_edge@units = "degrees_north"
forrest@0
    73
forrest@0
    74
;==================================================
forrest@0
    75
; area
forrest@0
    76
;==================================================
forrest@0
    77
forrest@0
    78
  deg2rad = 4.0*atan(1.0)/180.0
forrest@0
    79
  earth_rad = 6.37122e6
forrest@0
    80
  dx = earth_rad * deg2rad * (xi_edge(1)-xi_edge(0))
forrest@0
    81
  dsin = sin(deg2rad*yi_edge(1:180)) - sin(deg2rad*yi_edge(0:179))
forrest@0
    82
  area_in = new((/ nlatx, nlonx /), float)
forrest@0
    83
  area_in = conform(area_in, earth_rad*dx*dsin,0)
forrest@0
    84
  area_in@units = "m^2"
forrest@0
    85
forrest@0
    86
;**************************************************************
forrest@0
    87
; output grid: (0-360), (90S-90N)
forrest@0
    88
;**************************************************************
forrest@0
    89
; interpolate into model T31/T42, using model lat and lon
forrest@0
    90
forrest@0
    91
  diri  = "/fis/cgd/cseg/people/jeff/surface_data/"
forrest@0
    92
  fili  = "lnd_1.9.nc"   
forrest@0
    93
  b = addfile(diri+fili,"r")
forrest@0
    94
 
forrest@0
    95
  yo   = b->lat
forrest@0
    96
  xo   = b->lon
forrest@0
    97
  nlat = dimsizes(yo)
forrest@0
    98
  nlon = dimsizes(xo)       
forrest@0
    99
forrest@0
   100
  xo_edge = fspan(0.0, 360.0, nlon+1) - 0.5*360.0/nlon
forrest@0
   101
  xo_edge@units = "degrees_east"
forrest@0
   102
forrest@0
   103
  yo_edge = new(nlat+1, float)
forrest@0
   104
  gau_info = gaus(nlat/2)
forrest@0
   105
  yo_edge(0) = -90.0
forrest@0
   106
  yo_edge(1:nlat-1) = doubletofloat(0.5*(gau_info(0:nlat-2,0)+gau_info(1:nlat-1,0)))
forrest@0
   107
  yo_edge(nlat) = 90.0
forrest@0
   108
  yo_edge@units = "degrees_north" 
forrest@0
   109
forrest@0
   110
; yo = 0.5*(yo_edge(0:nlat-1) + yo_edge(1:nlat))
forrest@0
   111
; yo@units = yo_edge@units
forrest@0
   112
; xo = 0.5*(xo_edge(0:nlon-1) + xo_edge(1:nlon))
forrest@0
   113
; xo@units = xo_edge@units
forrest@0
   114
  
forrest@0
   115
;===================================================
forrest@0
   116
; output
forrest@0
   117
;===================================================
forrest@0
   118
forrest@0
   119
  w  = new((/nyear,nmonth,nlat,nlon/),float)
forrest@0
   120
forrest@0
   121
  do m = 0,nyear-1
forrest@0
   122
  do n = 0,nmonth-1
forrest@0
   123
     w(m,n,:,:) = remap_areaave(xi_edge, yi_edge, bi(m,n,:,:), xo_edge, yo_edge)
forrest@0
   124
  end do
forrest@0
   125
  end do
forrest@0
   126
forrest@0
   127
  w!0   = "year"
forrest@0
   128
  w!1   = "month"               
forrest@0
   129
  w!2   = "lat" 
forrest@0
   130
  w&lat = yo            
forrest@0
   131
  w!3   = "lon"
forrest@0
   132
  w&lon = xo  
forrest@0
   133
  w@units      = bi@units
forrest@0
   134
  w@long_name  = bi@long_name
forrest@0
   135
  w@_FillValue = bi@_FillValue   
forrest@0
   136
 
forrest@0
   137
  c->FIRE_C = w
forrest@0
   138
  c->date   = g->date
forrest@0
   139
  c->area   = b->area
forrest@0
   140
 
forrest@0
   141
end