load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.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_ts") 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("combine") ) then  		 	; color of lines
      combine = 0
      if (parameter(23) .NE. "0") then
         combine = stringtointeger(parameter(23)) 
         if (stringtointeger(parameter(23)) .NE. 1) then
            print(" ")
            print("Please set 'combine' to 0 or 1")
            print(" ")
            exit
         end if
      end if
   end if
   if (combine .EQ. 1) then
      if( .not. isvar("c_var") ) then
         c_var=" "
         if (parameter(27) .NE. "c_variables") then
            c_var=parameter(27)
         end if
      end if
   end if
   if ( .not. isvar("black") ) then  		 	; color of lines
      black = 0
      if (parameter(31) .NE. "0") then
         black = stringtointeger(parameter(31)) 
         if (stringtointeger(parameter(31)) .NE. 1) then
            print(" ")
            print("Please set 'black' to 0 or 1")
            print(" ")
            exit
         end if
      end if
   end if
   if ( .not. isvar("dash") ) then   		 	; pattern of lines
      dash = 0
      if (parameter(29) .NE. "0") then
         dash = stringtointeger(parameter(29))
         if (stringtointeger(parameter(29)) .NE. 1) then
            print(" ")
            print("Please set 'dash' to 0 or 1")
            print(" ")
            exit
         end if  
      end if
   end if

   ; ***************************************************
   ; open input file
   ; ***************************************************

   f=addfile( file_in,"r")
   
   vNam  = getfilevarnames(f) 
   print(" ")
   print("Variable on netCDF file: " + vNam)
   print(" ")
   dim   = dimsizes(vNam)
   z_pr  = f->zpt
   dimz  = dimsizes(z_pr)
   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) .NE. "1") then
         if (parameter(13) .LE. "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

   ; ****************************************************
   ; set up legend and colors 
   ; ****************************************************
   
   legend_label=new(nt-1,double)
   do p=start_time_step,end_time_step
      legend_label(p-start_time_step)=t_all(p)
   end do
   
   np = end_time_step-start_time_step  
   if ( black .eq. 0 ) then
      res=True      
      res@xyLineColors = ispan(2,237,235/np)
   end if

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

   res@gsnDraw                 = False
   res@gsnFrame                = False
   res@gsnPaperOrientation     = "portrait"
   ;  res@gsnmaximize             = True
   res@gsnPaperWidth           = 8.27
   res@gsnPaperHeight          = 11.69
   res@gsnPaperMargin          = 0.79
   res@txFont                  = "helvetica"
   res@tiMainFont              = "helvetica"
   res@tiXAxisFont             = "helvetica"
   res@tiYAxisFont             = "helvetica"
   res@tmXBLabelFont           = "helvetica"
   res@tmYLLabelFont           = "helvetica"
   res@lgLabelFont             = "helvetica"
   res@tmLabelAutoStride       = True
   res@pmLegendDisplayMode     = "Always"
   res@pmLegendSide            = "Top"
   res@xyExplicitLegendLabels  = legend_label
   res@pmLegendParallelPosF    = 1.4
   res@pmLegendOrthogonalPosF  = -1.085
   res@pmLegendWidthF          = 0.12
   res@pmLegendHeightF         = 0.3
   res@lgLabelFontHeightF     = .02 
   ;  res@XBLabelConstantSpacingF = 1.0

   if ( dash .eq. 0 ) then
      res@xyMonoDashPattern       = True  
   end if 

   resP                        = True
   resP@txFont                 = "helvetica"
   resP@txString               = f@title
   resP@txFuncCode             = "~"
   resP@txFontHeightF          = 0.014

   ; ***************************************************
   ; set up colors and recourses for combined plot
   ; ***************************************************
   
   if (combine .EQ. 1) then

      if (.not. isvar("number_comb")) then
         m=0
         m=stringtointeger(parameter(25))
         if(m .EQ. 0) then
            print(" ")
            print("Please indicate the number of variables you would like to combine ('number_comb')")
            print(" ")
            exit
         end if
      end if
    
      ores=True
      
      colors=new(m*(nt-1),integer)
      do j=0,m-1
         colors(j+j*(nt-2):(j+1)*(nt-1)-1) = ispan(2,237,235/np)
      end do   
      ores@xyLineColors = colors

      dash_oplot=new(m*(nt-1),integer)
      do j=0,m-1
         dash_oplot(j+j*(nt-2):(j+1)*(nt-1)-1)=j
      end do
      ores@xyDashPatterns = dash_oplot

      ores                         = True
      ores@gsnDraw                 = False
      ores@gsnFrame                = False
      ores@gsnPaperOrientation     = "portrait"
      ;  ores@gsnmaximize             = True
      ores@gsnPaperWidth           = 8.27
      ores@gsnPaperHeight          = 11.69
      ores@gsnPaperMargin          = 0.79
      ores@txFont                  = "helvetica"
      ores@tiMainFont              = "helvetica"
      ores@tiXAxisFont             = "helvetica"
      ores@tiYAxisFont             = "helvetica"
      ores@tmXBLabelFont           = "helvetica"
      ores@tmYLLabelFont           = "helvetica"
      ores@lgLabelFont             = "helvetica"
      ores@tmLabelAutoStride       = True
      ores@pmLegendDisplayMode     = "Always"
      ores@pmLegendSide            = "Top"
      ores@pmLegendParallelPosF    = 1.4
      ores@pmLegendOrthogonalPosF  = -1.085
      ores@pmLegendWidthF          = 0.15
      ores@pmLegendHeightF         = 0.60
      ores@lgLabelFontHeightF     = .02 
      ;  ores@XBLabelConstantSpacingF = 1.0
 
   end if

   ; ***************************************************
   ; set up graphics for plot
   ; ***************************************************

   plot  = new(dim,graphic)

   wks=gsn_open_wks(format_out,file_out)
   gsn_define_colormap(wks,"rainbow+white") 
 
   ; ***************************************************
   ; indicate plot number
   ; ***************************************************

   if (combine .EQ. 1) then
      n = 1 
   else
      n = 0
   end if 
  
   ; ***************************************************
   ; set up minimum and maximum height
   ; ***************************************************

   if (.not. isvar("min_z"))
      min_z=0
      if (stringtointeger(parameter(33)) .NE. 0) then
         if (stringtointeger(parameter(33)) .GE. max(z_pr) ) then
            print(" ")
            print("Minimum of height ('min_z'="+stringtointeger(parameter(33))+") is greater than available heights (="+max(z_pr)+")")
            print(" ")
            exit
         end if
         if (stringtointeger(parameter(33)) .LT. 0 ) then
            print(" ")
            print("Begin minimum of height 'min_z' with 0")
            print(" ")
            exit
         end if 
         min_z=stringtointeger(parameter(33))
      end if
   else
      if (min_z .GE. max(z_pr) ) then
         print(" ")
         print("Minimum of height ('min_z'="+min_z+") is greater than available heights (="+max(z_pr)+")")
         print(" ")
         exit
      end if
      if (min_z .LT. 0 ) then
         print(" ")
         print("Begin minimum of height 'min_z' with 0")
         print(" ")
         exit
      end if
   end if

   if (.not. isvar("max_z"))
      max_z=max(z_pr)
      if ((parameter(35)) .NE. "max(z_pr)") then
         if (stringtofloat(parameter(35)) .GE. max(z_pr) ) then
            print(" ")
            print("Maximum of height ('max_z'="+parameter(35)+") is greater than available heights (="+max(z_pr)+")")
            print(" ")
            exit
         end if
         if (stringtointeger(parameter(35)) .LE. 0 ) then
            print(" ")
            print("Maximum of height 'max_z' should be at least 1")
            print(" ")
            exit
         end if 
         max_z=stringtointeger(parameter(35))
      end if
   else
      if (max_z .GE. max(z_pr) ) then
         print(" ")
         print("Maximum of height ('max_z'="+max_z+") is greater than available heights (="+max(z_pr)+")")
         print(" ")
         exit
      end if
      if (max_z .LE. 0 ) then
         print(" ")
         print("Maximum of height 'max_z' should be at least 1")
         print(" ")
         exit
      end if
   end if

   ; ***************************************************
   ; read data and create plots
   ; ***************************************************
      
   do ti = start_time_step, end_time_step
      if( t_all(ti) .lt. 10^36) then
         start_time_step = ti
         break
      end if
   end do  
   
   if (combine .EQ. 1) then
      data_o=new((/m,end_time_step,dimz/),float)
      data_o!0 = "e"
      legend_label_oplot=new(m*(nt-1),string)
      mini=new(m,float)
      maxi=new(m,float)
      e=-1
   end if  
    
   do varn = 0, dim-1   
     
      if ( isStrSubset( vNam(varn), "NORM") .or. isStrSubset( vNam(varn), "time") )
         check = False
      else
         if (.not. isvar("var")) then
            check = True
            if (parameter(21) .NE. "variables") then
               var=parameter(21)
               check = isStrSubset( var,vNam(varn)+"," )
            end if
         else          
            check = isStrSubset( var,vNam(varn)+"," )
         end if
      end if
     
      if (combine .EQ. 1) then
         if (c_var .EQ. " ") then
            print(" ")
            print("Please indicate the variables you would like to combine ('c_var')")
            print(" ")
            exit
         end if               
         com=isStrSubset(c_var,vNam(varn)+",")      
         if (com) then
            e=e+1
            check = False
            temp = f->$vNam(varn)$            
            data_o!1 = temp!0
            data_o!2 = temp!1           
            data_o(e,:,:) = temp(start_time_step:end_time_step,:)
            mini(e) = min(data_o(e,:,:)) 
            maxi(e) = max(data_o(e,:,:))
            
            do j=e,e
               do p=start_time_step,end_time_step
                  legend_label_oplot(e*end_time_step+p-start_time_step)=" "+flt2string(doubletofloat(t_all(p)))+" "+vNam(varn)
               end do
            end do
           
            delete(temp)
            varn = varn + 1
    
         end if        
      end if  
 
      if(check) then

         temp = f->$vNam(varn)$      
         data = temp(start_time_step:end_time_step,:)      
         print("         plot of " + vNam(varn))
         data!1 = "z"
  
         res@gsnLeftString      = "PROFILE plot of"
         res@gsnRightString     = vNam(varn)
         res@txFontHeightF      = 0.02
         res@tiXAxisFontHeightF = 0.02
         res@tiYAxisFontHeightF = 0.02
         res@tiXAxisString      = " "    
         res@tiYAxisString      = "Height [z]"
         res@trXMinF            = min(data(:,:))
         res@trXMaxF            = max(data(:,:))
         res@trYMinF            = min_z
         res@trYMaxF            = max_z    
        
         plot(n) = gsn_csm_xy(wks,data,z_pr,res)  
      
         delete(temp)
         delete(data)
         n = n + 1 
         varn = varn + 1
       
      end if      

   end do

   if (combine .EQ. 1) then
      data_all=new((/m*end_time_step,dimz/),float)
      do j=0,e
         data_all(start_time_step-1+j*end_time_step:(j+1)*end_time_step-1,:) = data_o(j,:,:)
      end do  

      ores@gsnLeftString      = "combined PROFILE plot of" 
      ores@gsnRightString     = c_var
      ores@txFontHeightF      = 0.02
      ores@tiXAxisFontHeightF = 0.02
      ores@tiYAxisFontHeightF = 0.02
      ores@tiXAxisString      = " "    
      ores@tiYAxisString      = "Height [z]"
      ores@trXMinF            = min(mini(:))
      ores@trXMaxF            = max(maxi(:))
      ores@trYMinF            = min_z
      ores@trYMaxF            = max_z   
      ores@xyExplicitLegendLabels  = legend_label_oplot

      print(" ")
      print("combined plot of " + c_var) 
 
      plot(0) = gsn_csm_xy(wks,data_all,z_pr,ores) 
      
   end if

   ; ***************************************************
   ; merge plots onto one page 
   ; ***************************************************
  
   do i = 0, n-1, no_lines*no_columns 
      if( (i+no_lines*no_columns) .gt. (n-1)) then
        gsn_panel(wks,plot(i:n-1),(/no_lines,no_columns/),resP)
      else
        gsn_panel(wks,plot(i:i+no_lines*no_columns-1),(/no_lines,no_columns/),resP)
      end if
   end do

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

end