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_0.25deg.ncl.0
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 0.25degree

     3 ; no change in lat and lon, to be in consistent with 

     4 ; land class data

     5 ; ***********************************************

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

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

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

     9 ;************************************************

    10 begin

    11 

    12   year = 2003

    13 

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

    15 ; output data

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

    17   diro  = "/fis/cgd/cseg/people/jeff/clamp_data/lai/"

    18   filo  = "LAI_2000-2005_mean_0.25.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/lai/"

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

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

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

    29   bi    = g->LAI   

    30   xi    = g->lon 

    31   yi    = g->lat

    32 

    33 ;************************************************

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

    35 ; Observed NPP*scale_factor

    36 ;************************************************

    37  

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

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

    40  printVarSummary(bi)

    41 

    42  b2    = bi

    43  x2    = xi   

    44  

    45  nx = dimsizes(xi)

    46  do i= 0,nx-1

    47     if (i .lt. 3600) then

    48        p = i + 3600

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

    50     else

    51        p = i - 3600

    52        xi(p) = x2(i)

    53     end if

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

    55  end do

    56 

    57  bi&lat =  yi

    58  bi&lon =  xi

    59 

    60 ;print (xi)

    61 ;print (yi)

    62 ;exit

    63 

    64 ;************************************************

    65 ; read in model data

    66 ;************************************************

    67 ;fili2  = "i01.03cn_1545-1569_ANN_climo.nc"

    68  fili2  = ""

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

    70 

    71  lon    = f->lon     

    72  lat    = f->lat

    73  nlon   = dimsizes(lon)

    74  nlat   = dimsizes(lat)      

    75 

    76 ; print (xi)

    77 ; print (yi)

    78 ; print (xo)

    79 ; print (yo)

    80 

    81 ;(A)

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

    83 

    84 ;(B)

    85 

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

    87 ; clat  = lat

    88 ; clat  = cos(lat*rad)

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

    90 ; delete(clat@units)

    91 ; printVarSummary(clat)

    92 

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

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

    95  do m = 0,11

    96     do j=0,nlat-1

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

    98           if (j.eq.0) then

    99              LATS = -90.          

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

   101           end if

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

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

   104              LATN = 90.                  

   105           end if

   106        else

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

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

   109        end if

   110  

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

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

   113  

   114       do i=0,nlon-1

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

   116           if (i.eq.0) then

   117              LONL = 0.          

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

   119           end if

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

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

   122              LONR = 360.                 

   123           end if

   124        else

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

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

   127        end if

   128 

   129 ;print (LATS)

   130 ;print (LATN)

   131 ;print (LONL)

   132 ;print (LONR)

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

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

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

   136       end do

   137  

   138 ;     delete(CLAT)

   139 ;     delete(TEMP) 

   140     end do

   141   end do

   142 

   143   bo!0   = "time"

   144   bo!1   = "lat"

   145   bo!2   = "lon"

   146   bo&time= bi&time

   147 ; bo&time= 1

   148   bo&lat = lat

   149   bo&lon = lon

   150 ; bo@units      = bi@units

   151 ; bo@long_name  = bi@long_name

   152   bo@units      = "none"

   153   bo@long_name  = "Leaf Area Index"

   154   bo@_FillValue = bi@_FillValue

   155 

   156   c->LAI  = bo

   157 end
clamp