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