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/33.contour_diff_phase.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 ; 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 = "Month of 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) = maxind(s) + 1

    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 = "Month of 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) = maxind(s) + 1

    60   end do

    61   end do

    62 

    63   print (min(x)+"/"+max(x))

    64   d = x              

    65   d = x - y

    66   d = where((d .gt. 6.),12.-d,d)

    67   d = where((d .lt. -6.),-12.-d,d)

    68   print (min(d)+"/"+max(d))

    69 

    70   delete (z)

    71   delete (s)

    72 

    73   delta = 0.000001

    74   y = where(ismissing(y).and.(ismissing(x).or.(x.lt.delta)),0.,y)

    75 ;************************************************

    76 ; create default plot

    77 ;************************************************

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

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

    80 

    81   res                     = True             ; Use plot options

    82   res@cnFillOn            = True             ; Turn on color fill

    83   res@gsnSpreadColors     = True             ; use full colormap

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

    85 ; res@lbLabelAutoStride   = True

    86   res@cnLinesOn           = False            ; Turn off contourn lines

    87   res@mpFillOn            = False            ; Turn off map fill

    88 

    89   res@gsnSpreadColors      = True            ; use full colormap

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

    91   res@cnMinLevelValF       = 1.              ; Min level

    92   res@cnMaxLevelValF       = 12.             ; Max level

    93   res@cnLevelSpacingF      = 1.              ; interval

    94 

    95   pres                     = True            ; panel plot mods desired

    96   pres@gsnMaximize         = True            ; fill the page

    97 

    98   plot=new(3,graphic)                        ; create graphic array

    99 

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

   101   plot(0) = gsn_csm_contour_map_ce(wks,y,res)   

   102 

   103   res@tiMainString        = "Model i01.03cn"

   104 ; res@tiMainString        = "Model i01.04casa"

   105   plot(1) = gsn_csm_contour_map_ce(wks,x,res)   

   106  

   107   res@tiMainString        = "(Model i01.03cn) - (Observed)"

   108 ; res@tiMainString        = "(Model i01.04casa) - (Observed)"    

   109   res@cnMinLevelValF       = -6.              ; Min level

   110   res@cnMaxLevelValF       =  6.             ; Max level

   111   res@cnLevelSpacingF      =  1.              ; interval

   112   plot(2) = gsn_csm_contour_map_ce(wks,d,res)   ; for observed

   113 

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

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

   116 end
clamp