Mercurial Mercurial > repos > public_software > clamp / file revision
summary | shortlog | changelog | graph | tags | bookmarks | branches | files | changeset | file | latest | revisions | annotate | diff | comparison | raw | help
lai/03.lai_ob_to_T42.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 
     1 ; ***********************************************

     2 ; interpolate into model grids (T31)

     3 ; ***********************************************

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

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

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

     7 ;************************************************

     8 begin

     9 

    10   year = 2003

    11 

    12 ;************************************************

    13 ; output data

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

    15   diro  = "/fis/cgd/cseg/people/jeff/clamp_data/"

    16 ; filo  = "LAI_" + year + "_monthly_T42.nc"

    17 ; filo  = "LAI_2000-2005_mean_T42.nc"

    18   filo  = "LAI_2000-2005_ensemble_T42.nc"

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

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

    21 

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

    23 ; read in observed data

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

    25   diri  = "/fis/cgd/cseg/people/jeff/clamp_data/"

    26 ; fili  = "LAI_" + year + "_monthly.nc"

    27 ; fili  = "LAI_2000-2005_mean_0.05deg.nc"

    28   fili  = "LAI_2000-2005_ensemble_0.05deg.nc"

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

    30   bi    = g->LAI   

    31   xi    = g->lon 

    32   yi    = g->lat

    33 

    34 ;************************************************

    35 ; change into 0-360E, 90S-90N

    36 ; Observed NPP*scale_factor

    37 ;************************************************

    38  

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

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

    41  printVarSummary(bi)

    42 

    43  b2    = bi

    44  x2    = xi   

    45  

    46  nx = dimsizes(xi)

    47  do i= 0,nx-1

    48     if (i .lt. 3600) then

    49        p = i + 3600

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

    51     else

    52        p = i - 3600

    53        xi(p) = x2(i)

    54     end if

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

    56  end do

    57 

    58  bi&lat =  yi

    59  bi&lon =  xi

    60 

    61 ;print (xi)

    62 ;print (yi)

    63 ;exit

    64 

    65 ;************************************************

    66 ; read in model data

    67 ;************************************************

    68 ;fili2  = "i01.04casa_1605-1629_ANN_climo.nc"

    69  fili2  = "i01.03cn_1545-1569_ANN_climo.nc"

    70  f     = addfile (diri+fili2,"r")

    71 

    72  lon    = f->lon     

    73  lat    = f->lat

    74  nlon   = dimsizes(lon)

    75  nlat   = dimsizes(lat)      

    76 

    77 ; print (xi)

    78 ; print (yi)

    79 ; print (xo)

    80 ; print (yo)

    81 

    82 ;(A)

    83 ;bo = linint2_Wrap(xi,yi,bi,True,xo,yo,0)

    84 

    85 ;(B)

    86 

    87 ; rad   = 4.*atan(1.)/180.

    88 ; clat  = lat

    89 ; clat  = cos(lat*rad)

    90 ; clat@long_name = "cos(latitude)"

    91 ; delete(clat@units)

    92 ; printVarSummary(clat)

    93 

    94  bo = new((/12,nlat,nlon/),float)

    95 ;bo = new((/1,nlat,nlon/),float)

    96  do m = 0,11

    97     do j=0,nlat-1

    98        if (j.eq.0 .or. j.eq.nlat-1) then

    99           if (j.eq.0) then

   100              LATS = -90.          

   101              LATN = lat(j)+0.5*(lat(j+1)-lat(j))

   102           end if

   103           if (j.eq.nlat-1) then

   104              LATS = lat(j)-0.5*(lat(j)-lat(j-1))

   105              LATN = 90.                  

   106           end if

   107        else

   108           LATS = lat(j)-0.5*(lat(j)-lat(j-1))

   109           LATN = lat(j)+0.5*(lat(j+1)-lat(j))

   110        end if

   111  

   112 ;      CLAT = clat({LATS:LATN})      ; do once for *slight* efficiency

   113 ;      TEMP = bi(:,{LATS:LATN},:)      ; 2D [lat,lon]

   114  

   115       do i=0,nlon-1

   116        if (i.eq.0 .or. i.eq.nlon-1) then

   117           if (i.eq.0) then

   118              LONL = 0.          

   119              LONR = lon(i)+0.5*(lon(i+1)-lon(i))

   120           end if

   121           if (i.eq.nlon-1) then

   122              LONL = lon(i)-0.5*(lon(i)-lon(i-1))

   123              LONR = 360.                 

   124           end if

   125        else

   126           LONL = lon(i)-0.5*(lon(i)-lon(i-1))

   127           LONR = lon(i)+0.5*(lon(i+1)-lon(i))

   128        end if

   129 

   130 ;print (LATS)

   131 ;print (LATN)

   132 ;print (LONL)

   133 ;print (LONR)

   134 ;        bo(:,j,i) = avg(bi(:,{LATS:LATN},{LONL:LONR}))

   135          bo(m,j,i) = avg(bi(m,{LATS:LATN},{LONL:LONR}))  

   136 ;        bo(:,j,i) = wgt_areaave(TEMP(:,{LONL:LONR}), CLAT, 1.0, 0)

   137       end do

   138  

   139 ;     delete(CLAT)

   140 ;     delete(TEMP) 

   141     end do

   142   end do

   143 

   144   bo!0   = "time"

   145   bo!1   = "lat"

   146   bo!2   = "lon"

   147   bo&time= bi&time

   148 ; bo&time= 1

   149   bo&lat = lat

   150   bo&lon = lon

   151 ; bo@units      = bi@units

   152 ; bo@long_name  = bi@long_name

   153   bo@units      = "none"

   154   bo@long_name  = "Leaf Area Index"

   155   bo@_FillValue = bi@_FillValue

   156 

   157   c->LAI  = bo

   158 end
clamp