source: palm/trunk/SCRIPTS/NCL/timeseries.ncl @ 155

load "$NCARG_ROOT/lib/ncarg/nclex/gsun/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"

begin
 
   ; ***************************************************
   ; read parameter_list
   ; ***************************************************
   
   if (isfilepresent("~/.ncl_preferences")) then
      parameter = asciiread("~/.ncl_preferences",63,"string")
      delete(parameter@_FillValue)
   else
      print(" ")
      print("Please copy '.ncl_preferences' into your $home dircetory")
      print(" ")
      exit
   end if

   ; ***************************************************
   ; set up default parameter values and strings if not assigned in prompt or parameter list
   ; ***************************************************
   
   if ( .not. isvar("file_in") ) then                   ; path+name of input file      
      if (parameter(7) .EQ. "input file") then
         print(" ")
         print("Please provide input file 'file_in = ' either in prompt or parameter_list")
         print(" ")
         exit
      else
         file_in = parameter(7)
      end if      
   end if
   if ( .not. isvar("format_out") ) then                ; format of output file 
      format_out = "x11"
      if (parameter(9) .NE. "x11") then
         format_out = parameter(9)  
      end if 
   end if
   if ( .not. isvar("file_out") ) then                  ; path+name of output file
      file_out = "test"
      if (parameter(11) .NE. "test") then
         file_out = parameter(11)  
      end if     
   end if
   if ( .not. isvar("no_columns") ) then                ; number of plots in one row
      no_columns = 1
      if (parameter(17) .NE. "1") then
         no_columns = stringtointeger(parameter(17))  
      end if
   end if
   if ( .not. isvar("no_lines") ) then                  ; number of plot-lines on one sheet
      no_lines = 2
      if (parameter(19) .NE. "2") then
         no_lines = stringtointeger(parameter(19))  
      end if
   end if  
   if ( .not. isvar("var") ) then                       ; variable name
      check = True
   end if 
  
   ; ***************************************************
   ; open input file
   ; ***************************************************

   f  = addfile(file_in , "r" )

   vNam  = getfilevarnames(f) 
   print(" ")
   print("Variable on netCDF file: " + vNam)
   print(" ")
   dim   = dimsizes(vNam)
   t_all = f->time
   nt    = dimsizes(t_all)

   ; ****************************************************       
   ; start of time step and different types of mistakes that could be done
   ; ****************************************************

   if ( .not. isvar("start_time_step") ) then           
      start_time_step = 1
      if (parameter(13) .LE. "1") then
         if (parameter(13) .EQ. "0")
            print(" ")
            print("Begin with time step 1")
            print(" ")
            exit
         end if
         if (stringtointeger(parameter(13)) .GE. nt)
            print(" ")
            print("'start_time_step' = "+ parameter(13) +" is greater than available time steps = " + (nt-1))
            print(" ")
            exit
         end if
         start_time_step = stringtointeger(parameter(13))  
      end if
   else 
      if (start_time_step .LE. 0)
         print(" ")
         print("Begin with time step 1")
         print(" ")
         exit
      end if
      if (start_time_step .GE. nt)
         print(" ")
         print("'start_time_step' = "+ start_time_step +" is greater than available time steps = " + (nt-1))
         print(" ")
         exit
      end if
   end if

   ; ****************************************************
   ; end of time step and different types of mistakes that could be done
   ; ****************************************************

   if ( .not. isvar("end_time_step") ) then             
      end_time_step = nt-1
      if (parameter(15) .NE. "nt-1") then
         if (parameter(15) .LE. "0")
            print(" ")
            print("'end_time_step' = "+parameter(15)+ " is too small; 'end_time_step' should be at least 1 ")
            print(" ")
            exit
         end if
         if (stringtointeger(parameter(15)) .GE. nt)
            print(" ")
            print("'end_time_step' = "+ parameter(15) +" is greater than available time steps = " + (nt-1))  
            print(" ")
            exit
         end if
         if (stringtointeger(parameter(15)) .LT. stringtointeger(parameter(13)) )
            print(" ")
            print("'end_time_step' = "+ parameter(15) +" is lower than 'start_time_step' = "+parameter(13))
            print(" ")
            exit
         end if
         end_time_step = stringtointeger(parameter(15))  
      end if    
   else
      if (end_time_step .LE. 0)
         print(" ")
         print("'end_time_step' = "+end_time_step+ " is too small; 'end_time_step' should be at least 1 ")
         print(" ")
         exit
      end if
      if (end_time_step .GE. nt)
         print(" ")
         print("'end_time_step' = "+ end_time_step +" is greater than available time steps = "+(nt-1))
         print(" ")
         exit
      end if
      if (end_time_step .LT. start_time_step)
            print(" ")
            print("'end_time_step' = "+end_time_step +" is lower than 'start_time_step' = "+start_time_step)
            print(" ")
            exit
      end if
   end if  

   t = f->time(start_time_step:end_time_step)

  
   ; ***************************************************
   ; set up recourses
   ; ***************************************************

   res                         = True
   res@gsnDraw                 = False
   res@gsnFrame                = False
   res@gsnPaperOrientation     = "portrait"
   res@gsnPaperWidth           = 8.27
   res@gsnPaperHeight          = 11.69
   res@gsnPaperMargin          = 0.79
   res@tmXBMode                = True
   res@tmYLMode                = True
   res@txFont                  = "helvetica"
   res@tiMainFont              = "helvetica"
   res@tiXAxisFont             = "helvetica"
   res@tiYAxisFont             = "helvetica"
   res@tmXBLabelFont           = "helvetica"
   res@tmYLLabelFont           = "helvetica"  
   res@xyLineColors            = (/237/)
   resP                        = True
   resP@txFont                 = "helvetica"
   resP@txString               = f@title+" time series "
   resP@txFuncCode             = "~" 
   resP@txFontHeightF          = 0.015

   res@vpWidthF=4


   ; ***************************************************
   ; read data and create plots
   ; ***************************************************
   
   wks_ps = gsn_open_wks(format_out,file_out)           
   gsn_define_colormap(wks_ps,"rainbow+white") 
   plot_ps=new(dim,graphic)                              
   n=0 
 
   do varn = 0, dim-1  
      if( isvar("var") ) then                   
         check = isStrSubset( var,vNam(varn)+"," )
      end if
      if (parameter(21) .NE. "variables") then
         var = parameter(21) 
         check = isStrSubset( var,vNam(varn)+"," ) 
      end if 
    
      if( isStrSubset (vNam(varn), "time") )
         check = False
      end if
  
      if(check) then
 
         n=n+1
         data = f ->$vNam(varn)$          
         print("plot of " + vNam(varn))

         res@gsnRightString       = vNam(varn)
         res@tiXAxisString        = " time [s] "
         res@tiXAxisFontHeightF   = 0.07
         res@txFontHeightF        = 0.07
         res@tiYAxisFontHeightF   = 0.07
      
         plot_ps(n) = gsn_csm_xy(wks_ps,t,data(start_time_step:end_time_step),res)
   
      end if
 
   end do

   ; ***************************************************
   ; merge plots onto one page 
   ; ***************************************************

   do np = 1,n,no_lines*no_columns
    
      if ( np + no_lines*no_columns .gt. n) then   
         gsn_panel(wks_ps, plot_ps(np:n),(/no_lines,no_columns/),resP)
      else
         gsn_panel(wks_ps, plot_ps(np:np+no_lines*no_columns-1),(/no_lines,no_columns/),resP)
      end if
     
   end do

   print(" ")
   print("Output to: " + file_out +"."+ format_out)
   print(" ")
 
end