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.land_class_to_T42.ncl
author Forrest Hoffman <forrest@climatemodeling.org>
Thu, 26 Mar 2009 14:02:21 -0400
changeset 1 4be95183fbcd
permissions -rw-r--r--
Modifications to scoring and graphics production for the final version of code for the C-LAMP paper in GCB. 
     1 ; ***********************************************

     2 ; interpolate into model grids (T42)

     3 ; output: lat: S->N    lon: 0->360

     4 ; input : lat: N->S    lon: -180->180 

     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   nclass = 20

    13 

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

    15 ; output data

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

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

    18 ; filo  = "land_class_T42.nc"

    19   filo  = "land_class_T42_new.nc"

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

    21 

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

    23 ; read in observed data

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

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

    26   fili  = "land_class_0.25deg.nc"

    27   fili  = "land_class_0.25deg_new.nc"

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

    29   bi    = g->LAND_CLASS   

    30   xi    = g->lon 

    31   yi    = g->lat

    32 

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

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

    35 ;************************************************

    36  

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

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

    39  printVarSummary(bi)

    40 

    41  b2    = bi

    42  x2    = xi   

    43  

    44  nx  = dimsizes(xi)

    45  nx2 = nx/2

    46 

    47  do i= 0,nx-1

    48     if (i .lt. nx2) then

    49        p = i + nx2

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

    51     else

    52        p = i - nx2

    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  diri2  = "/fis/cgd/cseg/people/jeff/clamp_data/model/"

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

    70  f     = addfile (diri2+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  bo = new((/nlat,nlon/),integer)

    83  count = new((/nclass/),integer)

    84 

    85   do j=0,nlat-1

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

    87         if (j.eq.0) then

    88            LATS = -90.          

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

    90         end if

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

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

    93            LATN = 90.                  

    94         end if

    95      else

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

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

    98      end if

    99  

   100      do i=0,nlon-1

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

   102            if (i.eq.0) then

   103               LONL = 0.          

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

   105            end if

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

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

   108               LONR = 360.                 

   109            end if

   110         else

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

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

   113         end if

   114 

   115 ;print (LATS)

   116 ;print (LATN)

   117 ;print (LONL)

   118 ;print (LONR)

   119 

   120         count = 0

   121 

   122         do n = 0,nclass-1

   123            count(n) = num (bi({LATS:LATN},{LONL:LONR}).eq.n)

   124         end do

   125  

   126         bo(j,i) = maxind(count)  

   127      end do     

   128   end do

   129 

   130   bo!0   = "lat"

   131   bo!1   = "lon"

   132   bo&lat = lat

   133   bo&lon = lon

   134   bo@units      = bi@units

   135   bo@long_name  = bi@long_name

   136   bo@_FillValue = bi@_FillValue

   137 

   138   print (min(bo) + "/" + max(bo))

   139 

   140   c->LAND_CLASS  = bo

   141 end
clamp