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_grow.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   day_of_data = (/31,28,31,30,31,30,31,31,30,31,30,31/)

    11 

    12 ;************************************************

    13 ; read in observed data

    14 ;************************************************

    15   diri1  = "/fis/cgd/cseg/people/jeff/clamp_data/lai/"

    16   fili1  = "LAI_2000-2005_ensemble_T42.nc"

    17   f     = addfile(diri1+fili1,"r")

    18 

    19   z     = f->LAI

    20   y     = z(0,:,:)

    21   y@long_name = "Days of Growing Season"

    22  

    23   dsizes_z = dimsizes(z)

    24   ntime    = dsizes_z(0)

    25   nlat     = dsizes_z(1)

    26   nlon     = dsizes_z(2)

    27   

    28   do j = 0,nlat-1

    29   do i = 0,nlon-1

    30      nday = 0.

    31      do k = 0,ntime-1

    32         if (.not. ismissing(z(k,j,i)) .and. z(k,j,i) .gt. 1.0) then

    33            nday = nday + day_of_data(k)

    34         end if

    35      end do

    36      if (nday .ne. 0.) then

    37         y(j,i) = nday

    38      end if

    39   end do

    40   end do

    41 

    42 ; printVarSummary(y)

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

    44 

    45   delete (z)

    46 ;************************************************

    47 ; read in data: model       

    48 ;************************************************

    49  diri2  = "/fis/cgd/cseg/people/jeff/clamp_data/model/"

    50 ;fili3  = "i01.03cn_1545-1569_MONS_climo.nc"

    51  fili3  = "i01.04casa_1605-1629_MONS_climo.nc"

    52  data_file_model = addfile(diri2+fili3,"r")

    53 

    54   z     = data_file_model->TLAI

    55   x     = z(0,:,:)

    56   x@long_name = "Days of Growing Season"  

    57  

    58   dsizes_z = dimsizes(z)

    59   ntime    = dsizes_z(0)

    60   nlat     = dsizes_z(1)

    61   nlon     = dsizes_z(2)

    62   

    63   do j = 0,nlat-1

    64   do i = 0,nlon-1

    65      nday = 0.

    66      do k = 0,ntime-1

    67         if (.not. ismissing(z(k,j,i)) .and. z(k,j,i) .gt. 1.0) then

    68            nday = nday + day_of_data(k)

    69         end if

    70      end do

    71     if (nday .ne. 0.) then

    72        x(j,i) = nday

    73     end if

    74   end do

    75   end do

    76 

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

    78 

    79   d = x              

    80   d = x - y

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

    82 

    83   delete (z)

    84 

    85   delta = 0.000001

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

    87 ;************************************************

    88 ; create default plot

    89 ;************************************************

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

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

    92 

    93   res                     = True             ; Use plot options

    94   res@cnFillOn            = True             ; Turn on color fill

    95   res@gsnSpreadColors     = True             ; use full colormap

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

    97 ; res@lbLabelAutoStride   = True

    98   res@cnLinesOn           = False            ; Turn off contourn lines

    99   res@mpFillOn            = False            ; Turn off map fill

   100 

   101   res@gsnSpreadColors      = True            ; use full colormap

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

   103   res@cnMinLevelValF       = 0.              ; Min level

   104   res@cnMaxLevelValF       = 365.             ; Max level

   105   res@cnLevelSpacingF      = 30.              ; interval

   106 

   107   pres                     = True            ; panel plot mods desired

   108   pres@gsnMaximize         = True            ; fill the page

   109 

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

   111 

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

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

   114 

   115 ; res@tiMainString        = "Model i01.03cn"

   116   res@tiMainString        = "Model i01.04casa"

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

   118  

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

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

   121   res@cnMinLevelValF       = -120.              ; Min level

   122   res@cnMaxLevelValF       =  120.             ; Max level

   123   res@cnLevelSpacingF      =  20.              ; interval

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

   125 

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

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

   128 end
clamp