lai/03.lai_ob_to_0.25deg.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 (T31)
forrest@0
     3
; ***********************************************
forrest@0
     4
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
forrest@0
     5
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
forrest@0
     6
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
forrest@0
     7
;************************************************
forrest@0
     8
begin
forrest@0
     9
forrest@0
    10
  year = 2000
forrest@0
    11
forrest@0
    12
;************************************************
forrest@0
    13
; output data
forrest@0
    14
;************************************************
forrest@0
    15
  diro  = "/fis/cgd/cseg/people/jeff/clamp_data/lai/"
forrest@0
    16
; filo  = "LAI_" + year + "_monthly_0.25deg.nc"
forrest@0
    17
  filo  = "LAI_2000-2005_mean_0.25deg.nc"
forrest@0
    18
; filo  = "LAI_2000-2005_ensemble_0.25deg.nc"
forrest@0
    19
  c = addfile(diro+filo,"c")
forrest@0
    20
  filedimdef(c,"time",-1,True)
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  = "LAI_" + year + "_monthly.nc"
forrest@0
    27
  fili  = "LAI_2000-2005_mean_0.05deg.nc"
forrest@0
    28
; fili  = "LAI_2000-2005_ensemble_0.05deg.nc"
forrest@0
    29
  g     = addfile (diri+fili,"r")
forrest@0
    30
  bi    = g->LAI   
forrest@0
    31
  xi    = g->lon 
forrest@0
    32
  yi    = g->lat
forrest@0
    33
forrest@0
    34
;************************************************
forrest@0
    35
; change into 0-360E, 90S-90N
forrest@0
    36
; Observed NPP*scale_factor
forrest@0
    37
;************************************************
forrest@0
    38
 
forrest@0
    39
 yi    = (/ yi(::-1) /)
forrest@0
    40
 bi    = (/ bi(:,::-1,:) /)
forrest@0
    41
 printVarSummary(bi)
forrest@0
    42
forrest@0
    43
 b2    = bi
forrest@0
    44
 x2    = xi   
forrest@0
    45
 
forrest@0
    46
 nx = dimsizes(xi)
forrest@0
    47
 do i= 0,nx-1
forrest@0
    48
    if (i .lt. 3600) then
forrest@0
    49
       p = i + 3600
forrest@0
    50
       xi(p) = x2(i) + 360.      
forrest@0
    51
    else
forrest@0
    52
       p = i - 3600
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
; create 0.25deg lat and lon
forrest@0
    67
;************************************************
forrest@0
    68
  nlon   = 360*4
forrest@0
    69
  nlat   = 180*4
forrest@0
    70
forrest@0
    71
  lat  = latGlobeFo(nlat, "lat", "latitude", "degrees_north") ; S->N
forrest@0
    72
; lat  = (/ lat(::-1) /)                                      ; N->S  
forrest@0
    73
  lat&lat = lat
forrest@0
    74
  lon  = lonGlobeFo(nlon, "lon", "longitude", "degrees_east") ; 0->360 
forrest@0
    75
; lon  = (/ lon - 180. /) ; subtract 180 from all values      ; 180W-180E
forrest@0
    76
  lon&lon = lon           ; update coordinates
forrest@0
    77
forrest@0
    78
; print (lon)
forrest@0
    79
; print (lat)
forrest@0
    80
forrest@0
    81
; rad   = 4.*atan(1.)/180.
forrest@0
    82
; clat  = lat
forrest@0
    83
; clat  = cos(lat*rad)
forrest@0
    84
; clat@long_name = "cos(latitude)"
forrest@0
    85
; delete(clat@units)
forrest@0
    86
; printVarSummary(clat)
forrest@0
    87
forrest@0
    88
;bo = new((/12,nlat,nlon/),float)
forrest@0
    89
 bo = new((/1,nlat,nlon/),float)
forrest@0
    90
forrest@0
    91
;do m = 0,11
forrest@0
    92
 do m = 0,0
forrest@0
    93
    do j=0,nlat-1
forrest@0
    94
       if (j.eq.0 .or. j.eq.nlat-1) then
forrest@0
    95
          if (j.eq.0) then
forrest@0
    96
             LATS = -90.          
forrest@0
    97
             LATN = lat(j)+0.5*(lat(j+1)-lat(j))
forrest@0
    98
          end if
forrest@0
    99
          if (j.eq.nlat-1) then
forrest@0
   100
             LATS = lat(j)-0.5*(lat(j)-lat(j-1))
forrest@0
   101
             LATN = 90.                  
forrest@0
   102
          end if
forrest@0
   103
       else
forrest@0
   104
          LATS = lat(j)-0.5*(lat(j)-lat(j-1))
forrest@0
   105
          LATN = lat(j)+0.5*(lat(j+1)-lat(j))
forrest@0
   106
       end if
forrest@0
   107
 
forrest@0
   108
;      CLAT = clat({LATS:LATN})      ; do once for *slight* efficiency
forrest@0
   109
;      TEMP = bi(:,{LATS:LATN},:)      ; 2D [lat,lon]
forrest@0
   110
 
forrest@0
   111
      do i=0,nlon-1
forrest@0
   112
       if (i.eq.0 .or. i.eq.nlon-1) then
forrest@0
   113
          if (i.eq.0) then
forrest@0
   114
             LONL = 0.          
forrest@0
   115
             LONR = lon(i)+0.5*(lon(i+1)-lon(i))
forrest@0
   116
          end if
forrest@0
   117
          if (i.eq.nlon-1) then
forrest@0
   118
             LONL = lon(i)-0.5*(lon(i)-lon(i-1))
forrest@0
   119
             LONR = 360.                 
forrest@0
   120
          end if
forrest@0
   121
       else
forrest@0
   122
          LONL = lon(i)-0.5*(lon(i)-lon(i-1))
forrest@0
   123
          LONR = lon(i)+0.5*(lon(i+1)-lon(i))
forrest@0
   124
       end if
forrest@0
   125
forrest@0
   126
;print (LATS)
forrest@0
   127
;print (LATN)
forrest@0
   128
;print (LONL)
forrest@0
   129
;print (LONR)
forrest@0
   130
forrest@0
   131
         bo(m,j,i) = avg(bi(m,{LATS:LATN},{LONL:LONR}))  
forrest@0
   132
;        bo(m,j,i) = wgt_areaave(TEMP(m,{LONL:LONR}), CLAT, 1.0, 0)
forrest@0
   133
      end do
forrest@0
   134
 
forrest@0
   135
;     delete(CLAT)
forrest@0
   136
;     delete(TEMP) 
forrest@0
   137
    end do
forrest@0
   138
  end do
forrest@0
   139
forrest@0
   140
  bo!0   = "time"
forrest@0
   141
  bo!1   = "lat"
forrest@0
   142
  bo!2   = "lon"
forrest@0
   143
  bo&time= bi&time
forrest@0
   144
; bo&time= 1
forrest@0
   145
  bo&lat = lat
forrest@0
   146
  bo&lon = lon
forrest@0
   147
; bo@units      = bi@units
forrest@0
   148
; bo@long_name  = bi@long_name
forrest@0
   149
  bo@units      = "none"
forrest@0
   150
  bo@long_name  = "Leaf Area Index"
forrest@0
   151
  bo@_FillValue = bi@_FillValue
forrest@0
   152
forrest@0
   153
  c->LAI  = bo
forrest@0
   154
end