lai/33.contour_diff_max.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
; 3 plots, model - ob
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
; read in observed data
forrest@0
    11
;************************************************
forrest@0
    12
  diri1  = "/fis/cgd/cseg/people/jeff/clamp_data/lai/"
forrest@0
    13
  fili1  = "LAI_2000-2005_ensemble_T42.nc"
forrest@0
    14
  f     = addfile(diri1+fili1,"r")
forrest@0
    15
forrest@0
    16
  z     = f->LAI
forrest@0
    17
  y     = z(0,:,:)
forrest@0
    18
  y@long_name = "Leaf Area Index Max"
forrest@0
    19
  s     = z(:,0,0)
forrest@0
    20
 
forrest@0
    21
  dsizes_z = dimsizes(z)
forrest@0
    22
  ntime    = dsizes_z(0)
forrest@0
    23
  nlat     = dsizes_z(1)
forrest@0
    24
  nlon     = dsizes_z(2)
forrest@0
    25
  
forrest@0
    26
  do j = 0,nlat-1
forrest@0
    27
  do i = 0,nlon-1
forrest@0
    28
     s      = z(:,j,i) 
forrest@0
    29
     y(j,i) = max(s)
forrest@0
    30
  end do
forrest@0
    31
  end do
forrest@0
    32
forrest@0
    33
; printVarSummary(y)
forrest@0
    34
  print (min(y)+"/"+max(y))
forrest@0
    35
forrest@0
    36
  delete (s)
forrest@0
    37
  delete (z)
forrest@0
    38
;************************************************
forrest@0
    39
; read in data: model       
forrest@0
    40
;************************************************
forrest@0
    41
 diri2  = "/fis/cgd/cseg/people/jeff/clamp_data/model/"
forrest@0
    42
;fili3  = "i01.03cn_1545-1569_MONS_climo.nc"
forrest@0
    43
 fili3  = "i01.04casa_1605-1629_MONS_climo.nc"
forrest@0
    44
 data_file_model = addfile(diri2+fili3,"r")
forrest@0
    45
forrest@0
    46
  z     = data_file_model->TLAI
forrest@0
    47
  x     = z(0,:,:)
forrest@0
    48
  x@long_name = "Leaf Area Index Max"
forrest@0
    49
  s     = z(:,0,0)  
forrest@0
    50
 
forrest@0
    51
  dsizes_z = dimsizes(z)
forrest@0
    52
  ntime    = dsizes_z(0)
forrest@0
    53
  nlat     = dsizes_z(1)
forrest@0
    54
  nlon     = dsizes_z(2)
forrest@0
    55
  
forrest@0
    56
  do j = 0,nlat-1
forrest@0
    57
  do i = 0,nlon-1
forrest@0
    58
     s      = z(:,j,i) 
forrest@0
    59
     x(j,i) = max(s)
forrest@0
    60
  end do
forrest@0
    61
  end do
forrest@0
    62
forrest@0
    63
  print (min(x)+"/"+max(x))
forrest@0
    64
  d = x              
forrest@0
    65
  d = x - y
forrest@0
    66
  print (min(d)+"/"+max(d))
forrest@0
    67
forrest@0
    68
  delete (z)
forrest@0
    69
  delete (s)
forrest@0
    70
forrest@0
    71
  delta = 0.000001
forrest@0
    72
  y = where(ismissing(y).and.(ismissing(x).or.(x.lt.delta)),0.,y)
forrest@0
    73
;************************************************
forrest@0
    74
; create default plot
forrest@0
    75
;************************************************
forrest@0
    76
  wks = gsn_open_wks("ps","xy")              ; open a ps file
forrest@0
    77
  gsn_define_colormap(wks,"gui_default")     ; choose colormap
forrest@0
    78
forrest@0
    79
  res                     = True             ; Use plot options
forrest@0
    80
  res@cnFillOn            = True             ; Turn on color fill
forrest@0
    81
  res@gsnSpreadColors     = True             ; use full colormap
forrest@0
    82
; res@cnFillMode          = "RasterFill"     ; Turn on raster color
forrest@0
    83
; res@lbLabelAutoStride   = True
forrest@0
    84
  res@cnLinesOn           = False            ; Turn off contourn lines
forrest@0
    85
  res@mpFillOn            = False            ; Turn off map fill
forrest@0
    86
forrest@0
    87
  res@gsnSpreadColors      = True            ; use full colormap
forrest@0
    88
  res@cnLevelSelectionMode = "ManualLevels"  ; Manual contour invtervals
forrest@0
    89
  res@cnMinLevelValF       = 0.              ; Min level
forrest@0
    90
  res@cnMaxLevelValF       = 10.             ; Max level
forrest@0
    91
  res@cnLevelSpacingF      = 1.              ; interval
forrest@0
    92
forrest@0
    93
  pres                     = True            ; panel plot mods desired
forrest@0
    94
  pres@gsnMaximize         = True            ; fill the page
forrest@0
    95
forrest@0
    96
  plot=new(3,graphic)                        ; create graphic array
forrest@0
    97
forrest@0
    98
  res@tiMainString        = "MODIS MOD 15A2 2000-2005"
forrest@0
    99
  plot(0) = gsn_csm_contour_map_ce(wks,y,res)   
forrest@0
   100
forrest@0
   101
; res@tiMainString        = "Model i01.03cn"
forrest@0
   102
  res@tiMainString        = "Model i01.04casa"
forrest@0
   103
  plot(1) = gsn_csm_contour_map_ce(wks,x,res)   
forrest@0
   104
 
forrest@0
   105
; res@tiMainString        = "(Model i01.03cn) - (Observed)"
forrest@0
   106
  res@tiMainString        = "(Model i01.04casa) - (Observed)"    
forrest@0
   107
  res@cnMinLevelValF       = -6.              ; Min level
forrest@0
   108
  res@cnMaxLevelValF       =  6.             ; Max level
forrest@0
   109
  res@cnLevelSpacingF      =  1.              ; interval
forrest@0
   110
  plot(2) = gsn_csm_contour_map_ce(wks,d,res)   ; for observed
forrest@0
   111
forrest@0
   112
  gsn_panel(wks,plot,(/3,1/),pres)           ; create panel plot
forrest@0
   113
  system("convert xy.ps xy.png")
forrest@0
   114
end