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",129,"string")
      delete(parameter@_FillValue)
   else
      if (isfilepresent("~/palm/current_version/trunk/SRIPTS/NCL/.ncl_preferences")) then
         parameter = asciiread("~/palm/current_version/trunk/SRIPTS/NCL/.ncl_preferences",129,"string")
         delete(parameter@_FillValue)
      else
         print(" ")
         print("'.ncl_preferences' does not exist in '~/palm/current_version/trunk/SRIPTS/NCL/'")
         print(" ")
         exit
      end if
   end if

   if ( .not. isvar("file_1") ) then     		  
      if (parameter(7) .EQ. "File in") then
         print(" ")
         print("Please provide 1st input file 'file_1=' either in prompt or parameter_list")
         print(" ")
         exit
      else
         file_in = parameter(7)
      end if   
   else
      file_in = file_1   
   end if
   if (.not. isfilepresent(file_in)) then
      print(" ")
      print("Your 1st input file: '"+file_in+"' does not exist")
      print(" ")
      exit
   end if 

   if ( .not. isvar("format_out") ) then   		
      format_out = "x11"
      if (parameter(9) .NE. "x11") then
         format_out = parameter(9)  
         if (format_out .NE. "x11" .AND. format_out .NE. "pdf" .AND. format_out .NE. "eps" .AND. format_out .NE. "ps" .AND. format_out .NE. "epsi" .AND. format_out .NE. "ncgm")then
            print(" ")
            print("Your 'format_out = "+format_out+"' is invalid and set to'x11'")
            print(" ")
         end if  
      end if 
   else
      if (format_out .NE. "x11" .AND. format_out .NE. "pdf" .AND. format_out .NE. "eps" .AND. format_out .NE. "ps" .AND. format_out .NE. "epsi" .AND. format_out .NE. "ncgm")then
         print(" ")
         print("Your 'format_out = "+format_out+"' is invalid and set to'x11'")
         print(" ")
      end if   
   end if

   if ( .not. isvar("file_out") ) then    	 
      file_out = "test"
      if (parameter(11) .NE. "test_ts") then
         file_out = parameter(11)  
     end if     
   end if
   if ( .not. isvar("no_columns") ) then    	
      no_columns = 1
      if (parameter(17) .NE. "1") then
         no_columns = stringtointeger(parameter(17))  
      end if
   end if
   if ( .not. isvar("no_lines") ) then     		
      no_lines = 2
      if (parameter(19) .NE. "2") then
         no_lines = stringtointeger(parameter(19))  
      end if
   end if 

   if (.not. isvar("logy"))then
      logy = 0
      if (stringtointeger(parameter(77)) .NE. 0) then
         logy = stringtointeger(parameter(77))
         if (stringtointeger(parameter(77)) .NE. 1) then
            print(" ")
            print("Your 'logy'= "+logy+" is invalid and set to 0")
            print(" ")
            logy = 0
         end if   
      end if
   else
      if (logy .NE. 0 .AND. logy .NE. 1)then
         print(" ")
         print("Your 'logy'= "+logy+" is invalid and set to 0")
         print(" ")
         logy = 0
      end if    
   end if

   if ( .not. isvar("sort") ) then   			
      sort = "layer"
      if (parameter(51) .NE. "layer") then
         sort = parameter(51)  
         if (sort .NE. "time") then
            print(" ")
            print("Your 'sort'= "+sort+" is invalid and set to 'layer'")
            print(" ")
            sort = "layer"  
         end if
      end if
   else
      if (sort .NE. "time" .OR. sort .NE. "layer")then
         print(" ")
         print("Your 'sort'= "+sort+" is invalid and set to 'layer'")
         print(" ")
         sort = "layer"    
      end if
   end if 

   if ( .not. isvar("black") ) then  		 	
      black = 0
      if (parameter(31) .NE. "0") then
         black = stringtointeger(parameter(31)) 
         if (stringtointeger(parameter(31)) .NE. 1) then
            print(" ")
            print("Your 'black'= "+black+" is invalid and set to 0")
            print(" ")
            black = 0
         end if
      end if
   else
      if (black .NE. 0 .AND. black .NE. 1)then
         print(" ")
         print("Your 'black'= "+black+" is invalid and set to 0")
         print(" ")
         black = 0
      end if
   end if
   if ( .not. isvar("dash") ) then   		 	
      dash = 0
      if (parameter(29) .NE. "0") then
         dash = stringtointeger(parameter(29))
         if (stringtointeger(parameter(29)) .NE. 1) then
            print(" ")
            print("Your 'dash'= "+dash+" is invalid and set to 0")
            print(" ")
            dash = 0
         end if  
      end if
   else
      if (dash .NE. 0 .AND. dash .NE. 1)then
         print(" ")
         print("Your 'dash'= "+dash+" is invalid and set to 0")
         print(" ")
         dash = 0
      end if
   end if

   if ( .not. isvar("norm") ) then     		
      norm = 1.0
      if (parameter(79) .NE. "1") then
         norm = stringtofloat(parameter(79)) 
         if (stringtofloat(parameter(79)) .EQ. 0) then 
            print(" ")
            print("You cannot normalise with 0, 'norm' is set to 1")
            print(" ")
            norm = 1.0
         end if
      end if
   else
      if (norm .EQ. 0) then 
         print(" ")
         print("You cannot normalise with 0, 'norm' is set to 1")
         print(" ")
         norm = 1.0
      end if
   end if 

   f=addfile(file_in,"r")
   
   vNam = getfilevarnames(f) 
   print(" ")
   print("Variable in input file: " + vNam)
   print(" ")
   dim = dimsizes(vNam)
   vDim = getfiledimsizes(f) 
 
   t_all = f->time
   nt    = dimsizes(t_all)
   delta_t=t_all(nt-1)/nt
   
   do i=0,dim-1
      if (vNam(i) .EQ. "zu_sp")then
         zu=f->zu_sp
         z=zu
         dimz=dimsizes(zu)
      else 
         if (vNam(i) .EQ. "zw_sp")then
            zw=f->zw_sp
            z=zw
            dimz=dimsizes(zw)
         end if
      end if
   end do
   
   ; ****************************************************	
   ; start of time step and different types of mistakes that could be done
   ; ****************************************************
   
   if ( .not. isvar("start_time_step") ) then   	
      start_time_step=t_all(0)/3600
      if (parameter(13) .NE. "t(0)") then
         if (stringtodouble(parameter(13)) .GT. t_all(nt-1)/3600)then
            print(" ")
            print("'start_time_step' = "+ parameter(13) +"h is greater than last time step = " + t_all(nt-1)+"s = "+t_all(nt-1)/3600+"h")
            print(" ")
            print("Please select another 'start_time_step'")
            print(" ")
            exit
         end if
         if (stringtofloat(parameter(13)) .LT. t_all(0)/3600)then
            print(" ")
            print("'start_time_step' = "+ parameter(13) +"h is lower than first time step = " + t_all(0)+"s = "+t_all(0)/3600+"h")           
            exit
         end if
         start_time_step=stringtodouble(parameter(13)) 
      end if
   else 
      if (start_time_step .GT. t_all(nt-1)/3600)then
         print(" ")
         print("'start_time_step' = "+ start_time_step +"h is greater than last time step = " + t_all(nt-1)+"s = "+t_all(nt-1)/3600+"h")
         print(" ")
         print("Please select another 'start_time_step'")
         print(" ")
         exit
      end if
      if (start_time_step .LT. t_all(0)/3600)then
         print(" ")
         print("'start_time_step' = "+ start_time_step +"h is lower than first time step = " + t_all(0)+"s = "+t_all(0)/3600+"h")        
         exit
      end if
   end if
   start_time_step = start_time_step*3600

   do i=0,nt-1     
      if (start_time_step .GE. t_all(i)-delta_t/2 .AND. start_time_step .LT. t_all(i)+delta_t/2)then
         st=i
         break
      end if
   end do
   
   ; ****************************************************
   ; end of time step and different types of mistakes that could be done
   ; ****************************************************

   if ( .not. isvar("end_time_step") ) then     	
      end_time_step = t_all(nt-1)/3600
      if (parameter(15) .NE. "t(end)") then
         if (stringtodouble(parameter(15)) .GT. t_all(nt-1)/3600)then
            print(" ")
            print("'end_time_step' = "+ parameter(15) +"h is greater than last time step = " + t_all(nt-1)+"s = "+t_all(nt-1)/3600+"h")
            print(" ")
            print("Please select another 'end_time_step'")  
            print(" ")
            exit
         end if
         if (stringtodouble(parameter(15)) .LT. start_time_step/3600)then
            print(" ")
            print("'end_time_step' = "+ parameter(15) +"h is lower than 'start_time_step' = "+start_time_step/3600+"h")
            print(" ")
            print("Please select another 'start_time_step' or 'end_time_step'")
            print(" ")
            exit
         end if
         end_time_step = stringtodouble(parameter(15))
      end if    
   else
      if (end_time_step .GT. t_all(nt-1)/3600)then
         print(" ")
         print("'end_time_step' = "+ end_time_step +"h is greater than last time step = " + t_all(nt-1)+"s = "+t_all(nt-1)/3600+"h")
         print(" ")
         print("Please select another 'end_time_step'")  
         print(" ")
         exit
      end if
      if (end_time_step .LT. start_time_step/3600)then
         print(" ")
         print("'end_time_step' = "+ end_time_step +"h is lower than 'start_time_step' = "+start_time_step+"h")
         print(" ")
         print("Please select another 'start_time_step' or 'end_time_step'")
         print(" ")
         exit
      end if
   end if
   end_time_step = end_time_step*3600

   do i=0,nt-1     
      if (end_time_step .GE. t_all(i)-delta_t/2 .AND. end_time_step .LT. t_all(i)+delta_t/2)then
         et=i
         break
      end if
   end do 

   delete(start_time_step)
   start_time_step=round(st,3)
   delete(end_time_step)
   end_time_step=round(et,3)

   print(" ")
   print("Output of time steps from "+t_all(start_time_step)/3600+" h = "+t_all(start_time_step)+" s => index = "+start_time_step)
   print("                     till "+t_all(end_time_step)/3600+" h = "+t_all(end_time_step)+" s => index = "+end_time_step)
   print(" ")

   dimt = end_time_step-start_time_step+1
 
   ; ***************************************************
   ; 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@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@pmLegendParallelPosF    = 1.15
   res@pmLegendOrthogonalPosF  = -1.0
   res@pmLegendWidthF          = 0.12
   if (sort .EQ. "time")then
      res@pmLegendHeightF         = 0.04*(end_time_step-start_time_step+1)
   else
      res@pmLegendHeightF         = 0.04*dimz
   end if
   res@lgLabelFontHeightF     = .02 
   res@lgTitleFontHeightF     = .02
   res@txFontHeightF      = 0.02
   res@tiXAxisFontHeightF = 0.02
   res@tiYAxisFontHeightF = 0.02
   
   if (logy .EQ. 1) then
      res@trXLog = True 
      res@trYLog = True 
   end if

   legend_label=new(dimt,double)
   legend_label_zu=new(dimz,double)
   legend_label_zw=new(dimz,double)
   legend_label_z=new(dimz,double)
   do p=start_time_step,end_time_step
      if (t_all(p)/3600 .LT. 1) then
         legend_label(p-start_time_step)=decimalPlaces(t_all(p)/3600,2,True)
      else
         legend_label(p-start_time_step)=decimalPlaces(t_all(p)/3600,0,True)
      end if
   end do  
   if (sort .EQ. "time")
      plot  = new(dim*dimz,graphic)
      plot_ = new(dim*dimz,graphic)
      np=dimt
      res@lgTitleString = "Time [h]"
   else
      plot  = new(dim*dimt,graphic)
      plot_ = new(dim*dimt,graphic)
      np=dimz
      res@lgTitleString = "Height [m]"  
      do p=0,dimz-1
         legend_label_zu(p)=round(zu(p),3)
         legend_label_zw(p)=round(zw(p),3)
      end do
   end if

   if ( black .eq. 0 ) then  
      res@xyLineColors = ispan(2,237,235/np)
   end if
   if ( dash .eq. 0 ) then
      res@xyMonoDashPattern       = True 
   end if

   wks=gsn_open_wks(format_out,file_out)
   gsn_define_colormap(wks,"rainbow+white") 

   n=0
   do varn =dim-1,0,1 

      check = True

      if ( isStrSubset( vNam(varn), "time") .OR. isStrSubset( vNam(varn), "zu_sp") .OR. isStrSubset( vNam(varn), "zw_sp") .OR. isStrSubset( vNam(varn), "k_x") .OR. isStrSubset( vNam(varn), "k_y")) then
            check = False
      end if

      if (.not. isvar("var")) then
         if (parameter(21) .NE. "variables") then
            var=parameter(21)
            check = isStrSubset( var,","+vNam(varn)+"," )
         end if
      else          
         check = isStrSubset( var,","+vNam(varn)+"," )
      end if

      if(check) then

         temp = f->$vNam(varn)$(start_time_step:end_time_step,0:dimz-1,:)

         temp=temp/norm  ;SMOOTHING: #temp=smth9(temp/norm, 0.50, -0.25, False)#
   
         a=getvardims(temp)
         b=dimsizes(a)            
         do i=0,b-1            
            if (isStrSubset( a(i),"zu_sp" ))then
               legend_label_z=legend_label_zu
            else
               if (isStrSubset( a(i),"zw_sp" ))then
                  legend_label_z=legend_label_zw     
               end if
            end if
         end do  

         if (isStrSubset(vNam(varn),"x"))then
            x_axis = f->k_x
            res@tiXAxisString = "k_x"
         else
            x_axis = f->k_y 
            res@tiXAxisString = "k_y"
         end if
        
         if (sort .EQ. "time")
            do p=dimz-1,0,1  
               do q=0,dimt-1
                  do r=0,dimsizes(x_axis)-1 
                     if (temp(q,p,r) .EQ. 0)then
                        st=p+start_time_step
                        print(" ")
                        print("'"+vNam(varn)+"("+st+","+q+","+r+")' is equal 0; Logarithmic scale for y-axis cannot be used")
                        print(" ")
                        res@trYLog = False
                     else
                        res@trYLog = True
                     end if
                  end do
               end do
               res@trYMinF = min(temp(:,p,:))
               res@trYMaxF = max(temp(:,p,:))
               res@trXMinF = min(x_axis)
               res@trXMaxF = max(x_axis)
               res@gsnLeftString      = vNam(varn)
               res@gsnRightString     = "Height: "+z(p)+"m"
               if (norm .NE. 1)then
                  res@tiYAxisString      = vNam(varn)+" / "+norm
               else
                  res@tiYAxisString      = vNam(varn)
               end if
               res@xyExplicitLegendLabels  = legend_label
               plot(n)  = gsn_csm_xy(wks,x_axis,temp(:,p,:),res)
               n=n+1
            end do
         else
            if (sort .EQ. "layer")
               do p=dimt-1,0,1 
                  do q=0,dimz-1
                     do r=0,dimsizes(x_axis)-1 
                        if (temp(p,q,r) .EQ. 0)then
                           st=p+start_time_step
                           print(" ")
                           print("'"+vNam(varn)+"("+st+","+q+","+r+")' is equal 0; Logarithmic scale for y-axis cannot be used")
                           print(" ")
                           res@trYLog = False
                        else
                           res@trYLog = True
                        end if
                     end do
                  end do
                  res@trYMinF = min(temp(p,:,:))
                  res@trYMaxF = max(temp(p,:,:))
                  res@trXMinF = min(x_axis)
                  res@trXMaxF = max(x_axis)
                  res@gsnLeftString      = vNam(varn)
                  res@gsnRightString     = "Time: "+legend_label(p)+"h"
                  if (norm .NE. 1)then
                     res@tiYAxisString      = vNam(varn)+" / "+norm
                  else
                     res@tiYAxisString      = vNam(varn)
                  end if
                  res@xyExplicitLegendLabels  = legend_label_z
                  plot(n)  = gsn_csm_xy(wks,x_axis,temp(p,:,:),res)
                  n=n+1
               end do
            else
               print(" ")
               print("Please choose 'sort' either equal 'time' or 'height'")
               print(" ")
               exit
            end if
         end if
         delete(temp)
      end if
   end do

   if (n .EQ. 0) then
      print(" ")
      print("The variables 'var=°"+var+"°' do not exist on your input file")
      print(" ")
      exit
   end if
  
   ; ***************************************************
   ; merge plots onto one page 
   ; ***************************************************

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

   do m=0,n-1
      plot_(m)=plot(n-1-m)
   end do

   if (format_out .EQ. "eps" .OR. format_out .EQ. "epsi") then
      gsn_panel(wks,plot_,(/n,1/),resP)
   else    
      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
   end if

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