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/31.contour_ob_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 
     1 ;*************************************************

     2 ; ce_1.ncl

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

    13   fili  = "LAI_2000-2005_ensemble_T42.nc"

    14 ; fili  = "LAI_2000_monthly_T42.nc"

    15   f     = addfile(diri+fili,"r")

    16 

    17   z     = f->LAI

    18   y     = z(0,:,:)

    19   y@long_name = "Leaf Area Index Max"

    20   p     = z(0,:,:)

    21   p@long_name = "Leaf Area Index Max Month"

    22   s     = z(:,0,0)

    23 ; printVarSummary(s)

    24  

    25   dsizes_z = dimsizes(z)

    26   ntime    = dsizes_z(0)

    27   nlat     = dsizes_z(1)

    28   nlon     = dsizes_z(2)

    29   

    30   do j = 0,nlat-1

    31   do i = 0,nlon-1

    32      s      = z(:,j,i) 

    33      y(j,i) = max(s)

    34      p(j,i) = maxind(s)+ 1

    35 

    36 ;    iMax   = ind(s.eq.y(j,i)) + 1      ; may be multiple max

    37 ;    if (.not.ismissing(iMax(0))) then 

    38 ;        p(j,i) =  iMax(0)          ; no meta data

    39 ;    end if

    40 ;    delete(iMax)

    41 

    42   end do

    43   end do

    44 

    45 ; printVarSummary(p)

    46   print (min(y)+"/"+max(y))

    47   print (min(p)+"/"+max(p))

    48 

    49 ;************************************************

    50 ; read in model data

    51 ;************************************************

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

    53   g     = addfile(diri+fili2,"r")

    54   x     = g->TLAI

    55 

    56   delta = 0.00000000001

    57   x0    = x(0,:,:)

    58   y = where(ismissing(y).and.(ismissing(x0).or.(x0.lt.delta)),0.,y)

    59 ;************************************************

    60 ; create default plot

    61 ;************************************************

    62   

    63 ; setvalues NhlGetWorkspaceObjectId()

    64 ;   "wsMaximumSize" : 199999999

    65 ; end setvalues

    66 

    67   wks = gsn_open_wks("ps","xy")         ; open a ps file

    68 ; gsn_define_colormap(wks,"wgne15")          ; choose colormap

    69   gsn_define_colormap(wks,"gui_default")     ; choose colormap

    70 

    71   res                     = True             ; Use plot options

    72   res@cnFillOn            = True             ; Turn on color fill

    73   res@gsnSpreadColors      = True            ; use full colormap

    74 ; res@cnFillMode          = "RasterFill"     ; Turn on raster color

    75 ; res@lbLabelAutoStride   = True

    76   res@cnLinesOn           = False            ; Turn off contourn lines

    77   res@mpFillOn            = False            ; Turn off map fill

    78   res@tiMainString        = "MODIS MOD 15A2 2000-2005"

    79 

    80   res@gsnSpreadColors     = True             ; use full colormap

    81   res@cnLevelSelectionMode = "ManualLevels"  ; Manual contour invtervals

    82   res@cnMinLevelValF       = 0.              ; Min level

    83   res@cnMaxLevelValF       = 10.             ; Max level

    84   res@cnLevelSpacingF      = 1.              ; interval

    85 

    86 ; pres                            = True     ; panel plot mods desired

    87 ; pres@gsnMaximize                = True     ; fill the page

    88 

    89   plot = gsn_csm_contour_map_ce(wks,y,res)   ; for observed

    90 

    91   res@cnMinLevelValF       = 1.              ; Min level

    92   res@cnMaxLevelValF       = 12.             ; Max level

    93   res@cnLevelSpacingF      = 1.              ; interval

    94   plot = gsn_csm_contour_map_ce(wks,p,res)   ; for observed

    95 

    96 ; gsn_panel(wks,plot,(/1,1/),pres)           ; create panel plot

    97   system("convert xy.ps xy.png")

    98 end
clamp