lai/03.land_class_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
forrest@0
     1
; ***********************************************
forrest@0
     2
; interpolate into model grids (T42)
forrest@0
     3
; output: lat: S->N    lon: 0->360
forrest@0
     4
; input : lat: N->S    lon: -180->180 
forrest@0
     5
; ***********************************************
forrest@0
     6
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
forrest@0
     7
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
forrest@0
     8
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
forrest@0
     9
;************************************************
forrest@0
    10
begin
forrest@0
    11
forrest@0
    12
  nclass = 20
forrest@0
    13
forrest@0
    14
;************************************************
forrest@0
    15
; output data
forrest@0
    16
;************************************************
forrest@0
    17
  diro  = "/fis/cgd/cseg/people/jeff/clamp_data/lai/"
forrest@0
    18
; filo  = "land_class_T42.nc"
forrest@0
    19
  filo  = "land_class_T42_new.nc"
forrest@0
    20
  c = addfile(diro+filo,"c")
forrest@0
    21
forrest@0
    22
;************************************************
forrest@0
    23
; read in observed data
forrest@0
    24
;************************************************
forrest@0
    25
  diri  = "/fis/cgd/cseg/people/jeff/clamp_data/lai/"
forrest@0
    26
  fili  = "land_class_0.25deg.nc"
forrest@0
    27
  fili  = "land_class_0.25deg_new.nc"
forrest@0
    28
  g     = addfile (diri+fili,"r")
forrest@0
    29
  bi    = g->LAND_CLASS   
forrest@0
    30
  xi    = g->lon 
forrest@0
    31
  yi    = g->lat
forrest@0
    32
forrest@0
    33
;************************************************
forrest@0
    34
; change into 0-360E, 90S-90N
forrest@0
    35
;************************************************
forrest@0
    36
 
forrest@0
    37
 yi    = (/ yi(::-1) /)
forrest@0
    38
 bi    = (/ bi(::-1,:) /)
forrest@0
    39
 printVarSummary(bi)
forrest@0
    40
forrest@0
    41
 b2    = bi
forrest@0
    42
 x2    = xi   
forrest@0
    43
 
forrest@0
    44
 nx  = dimsizes(xi)
forrest@0
    45
 nx2 = nx/2
forrest@0
    46
forrest@0
    47
 do i= 0,nx-1
forrest@0
    48
    if (i .lt. nx2) then
forrest@0
    49
       p = i + nx2
forrest@0
    50
       xi(p) = x2(i) + 360.      
forrest@0
    51
    else
forrest@0
    52
       p = i - nx2
forrest@0
    53
       xi(p) = x2(i)
forrest@0
    54
    end if
forrest@0
    55
    bi(:,p)= b2(:,i) 
forrest@0
    56
 end do
forrest@0
    57
forrest@0
    58
 bi&lat =  yi
forrest@0
    59
 bi&lon =  xi
forrest@0
    60
forrest@0
    61
 print (xi)
forrest@0
    62
 print (yi)
forrest@0
    63
;exit
forrest@0
    64
forrest@0
    65
;************************************************
forrest@0
    66
; read in model data
forrest@0
    67
;************************************************
forrest@0
    68
 diri2  = "/fis/cgd/cseg/people/jeff/clamp_data/model/"
forrest@0
    69
 fili2  = "i01.03cn_1545-1569_ANN_climo.nc"
forrest@0
    70
 f     = addfile (diri2+fili2,"r")
forrest@0
    71
forrest@0
    72
 lon    = f->lon     
forrest@0
    73
 lat    = f->lat
forrest@0
    74
 nlon   = dimsizes(lon)
forrest@0
    75
 nlat   = dimsizes(lat)      
forrest@0
    76
forrest@0
    77
; print (xi)
forrest@0
    78
; print (yi)
forrest@0
    79
; print (xo)
forrest@0
    80
; print (yo)
forrest@0
    81
forrest@0
    82
 bo = new((/nlat,nlon/),integer)
forrest@0
    83
 count = new((/nclass/),integer)
forrest@0
    84
forrest@0
    85
  do j=0,nlat-1
forrest@0
    86
     if (j.eq.0 .or. j.eq.nlat-1) then
forrest@0
    87
        if (j.eq.0) then
forrest@0
    88
           LATS = -90.          
forrest@0
    89
           LATN = lat(j)+0.5*(lat(j+1)-lat(j))
forrest@0
    90
        end if
forrest@0
    91
        if (j.eq.nlat-1) then
forrest@0
    92
           LATS = lat(j)-0.5*(lat(j)-lat(j-1))
forrest@0
    93
           LATN = 90.                  
forrest@0
    94
        end if
forrest@0
    95
     else
forrest@0
    96
        LATS = lat(j)-0.5*(lat(j)-lat(j-1))
forrest@0
    97
        LATN = lat(j)+0.5*(lat(j+1)-lat(j))
forrest@0
    98
     end if
forrest@0
    99
 
forrest@0
   100
     do i=0,nlon-1
forrest@0
   101
        if (i.eq.0 .or. i.eq.nlon-1) then
forrest@0
   102
           if (i.eq.0) then
forrest@0
   103
              LONL = 0.          
forrest@0
   104
              LONR = lon(i)+0.5*(lon(i+1)-lon(i))
forrest@0
   105
           end if
forrest@0
   106
           if (i.eq.nlon-1) then
forrest@0
   107
              LONL = lon(i)-0.5*(lon(i)-lon(i-1))
forrest@0
   108
              LONR = 360.                 
forrest@0
   109
           end if
forrest@0
   110
        else
forrest@0
   111
           LONL = lon(i)-0.5*(lon(i)-lon(i-1))
forrest@0
   112
           LONR = lon(i)+0.5*(lon(i+1)-lon(i))
forrest@0
   113
        end if
forrest@0
   114
forrest@0
   115
;print (LATS)
forrest@0
   116
;print (LATN)
forrest@0
   117
;print (LONL)
forrest@0
   118
;print (LONR)
forrest@0
   119
forrest@0
   120
        count = 0
forrest@0
   121
forrest@0
   122
        do n = 0,nclass-1
forrest@0
   123
           count(n) = num (bi({LATS:LATN},{LONL:LONR}).eq.n)
forrest@0
   124
        end do
forrest@0
   125
 
forrest@0
   126
        bo(j,i) = maxind(count)  
forrest@0
   127
     end do     
forrest@0
   128
  end do
forrest@0
   129
forrest@0
   130
  bo!0   = "lat"
forrest@0
   131
  bo!1   = "lon"
forrest@0
   132
  bo&lat = lat
forrest@0
   133
  bo&lon = lon
forrest@0
   134
  bo@units      = bi@units
forrest@0
   135
  bo@long_name  = bi@long_name
forrest@0
   136
  bo@_FillValue = bi@_FillValue
forrest@0
   137
forrest@0
   138
  print (min(bo) + "/" + max(bo))
forrest@0
   139
forrest@0
   140
  c->LAND_CLASS  = bo
forrest@0
   141
end