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 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("sort") ) then   			; sort of plots
      sort = "time"
      if (parameter(37) .NE. "time") then
         sort = parameter(37)  
      end if
   end if 
   if ( .not. isvar("mode") ) then   			; pattern of contour plots
      mode = "Fill"
      if (parameter(39) .NE. "Fill") then
         mode = parameter(39)  
      end if
   end if
   if ( .not. isvar("fill_mode") ) then 		; pattern of filling 
      fill_mode = "AreaFill"
      if (parameter(41) .NE. "AreaFill") then
         mode = parameter(41)  
      end if 
   end if
   if ( .not. isvar("shape") ) then   			; keeping of aspect ratio
      shape = 1
      if (parameter(43) .NE. "1") then
         shape = stringtointeger(parameter(43)) 
         if (stringtointeger(parameter(43)) .NE. 0) then
            print(" ")
            print("Please set 'shape' to 0 or 1")
            print(" ")
            exit
         end if 
      end if  
   end if  
   if ( .not. isvar("var") ) then     			; set output of all variables
      check = True
   end if

   if ( .not. isvar("xyc") ) then   			; turn xy cross-section on or off
      xyc = 0
      if (parameter(45) .NE. "0") then
         xyc = stringtointeger(parameter(45)) 
         if (stringtointeger(parameter(45)) .NE. 1) then
            print(" ")
            print("Please set 'xyc' to 0 or 1")
            print(" ")
            exit
         end if 
      end if
   end if 
   if ( .not. isvar("xzc") ) then   			; turn xz cross-section on or off 
      xzc = 0
      if (parameter(47) .NE. "0") then
         xzc = stringtointeger(parameter(47)) 
         if (stringtointeger(parameter(47)) .NE. 1) then
            print(" ")
            print("Please set 'xzc' to 0 or 1")
            print(" ")
            exit
         end if 
      end if
   end if
   if ( .not. isvar("yzc") ) then   			; turn yz cross-section on or off 
      yzc = 0
      if (parameter(49) .NE. "0") then
         yzc = stringtointeger(parameter(49)) 
         if (stringtointeger(parameter(49)) .NE. 1) then
            print(" ")
            print("Please set 'yzc' to 0 or 1")
            print(" ")
            exit
         end if 
      end if
   end if
   if ( .not. isvar("xs") ) then     			; start of x-coordinate range 
      xs = 0
      if (parameter(51) .NE. "0") then
         xs = stringtointeger(parameter(51))
      end if
   end if 
   if ( .not. isvar("ys") ) then     			; start of y-coordinate range
      ys = 0
      if (parameter(55) .NE. "0") then
         ys = stringtointeger(parameter(55))
      end if
   end if
   if ( .not. isvar("zs") ) then     			; start of z-coordinate range
      zs = 0
      if (parameter(59) .NE. "0") then
         zs = stringtointeger(parameter(59))
      end if
   end if  

   if (xyc .EQ. 0 .AND. xzc .EQ. 0 .AND. yzc .EQ. 0) then
      print(" ")
      print("Please select one crossection (xyc=1 or xzc=1 or yzc=1)")
      print(" ")
      exit
   end if
   if (xyc .EQ. 1 ) then 
      if (xzc .EQ. 1 .OR. yzc .EQ. 1) then
         print(" ")
         print("Please select just one crossection (xyc=1 or xzc=1 or yzc=1)")
         print(" ") 
         exit
      end if
   end if
   if (xzc .EQ. 1) then
      if (xyc .EQ. 1 .OR. yzc .EQ. 1) then
         print(" ")
         print("Please select just one crossection (xyc=1 or xzc=1 or yzc=1)")
         print(" ")
         exit
      end if
   end if
   if (yzc .EQ. 1) then
      if (xyc .EQ. 1 .OR. xzc .EQ. 1) then
         print(" ")
         print("Please select just one crossection (xyc=1 or xzc=1 or yzc=1)")
         print(" ")
         exit
      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)
  
   ; ***************************************************
   ; set up recourses
   ; ***************************************************

   cs_res                          = True
   cs_res@gsnYAxisIrregular2Linear = True  

   if( shape .eq. 1 ) then          			; keep the shape
      cs_res@gsnShape = True
   end if
  
   cs_res@gsnDraw                 = False
   cs_res@gsnFrame                = False
   cs_res@gsnMaximize             = True 
   cs_res@gsnPaperOrientation     = "portrait"
   cs_res@gsnPaperWidth           = 8.27
   cs_res@gsnPaperHeight          = 11.69
   cs_res@gsnPaperMargin          = 0.79
   cs_res@tmXBLabelFontHeightF   = .02
   cs_res@tmYLLabelFontHeightF   = .02
   cs_res@tiXAxisFontHeightF     = .02
   cs_res@tiYAxisFontHeightF     = .02
   cs_res@tmXBMode                ="Automatic"
   cs_res@tmYLMode                ="Automatic"
   cs_res@lgTitleFontHeightF     = .02 
   cs_res@lgLabelFontHeightF     = .02
   cs_res@txFontHeightF          = .02
   cs_resP = True
   cs_resP@txString               = f@title
   cs_resP@txFuncCode             = "~"         	; necessary for title
   cs_resP@txFontHeightF          = .02
  
   if ( mode .eq. "Fill" ) then
      cs_res@cnFillOn             = True
      cs_res@gsnSpreadColors      = True
      cs_res@cnFillMode           = fill_mode
      cs_res@lbOrientation        = "Vertical"  	; vertical label bar
      cs_res@cnLinesOn            = False	
      cs_res@cnLineLabelsOn       = False
   end if

   if ( mode .eq. "Both" ) then
      cs_res@cnFillOn             = True
      cs_res@gsnSpreadColors      = True
      cs_res@cnFillMode           = fill_mode
      cs_res@lbOrientation        = "Vertical"  	; vertical label bar
      cs_res@cnLinesOn            = True
      cs_res@cnLineLabelsOn       = True
   end if

   ; ***************************************************
   ; open workstation(s) 
   ; ***************************************************

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

   plot=new((/no_columns*no_lines/),graphic)

   page = 0

   nc1 = no_columns
   nl1 = no_lines

   print("plots sorted by '"+sort+"'")
   print(" ")

   do varn=0,dim-1    

      data_all = f->$vNam(varn)$

      if ( vNam(varn) .eq. "zw" .or. vNam(varn) .eq. "yv" .or. vNam(varn) .eq. "xu" ) then
         check = False
      end if
      if ( vNam(varn) .eq. "zwwi" .or. vNam(varn) .eq. "zusi") then
         check = False
      end if
      
      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(check) then

         print("plot of "+vNam(varn))
         print(" ")

         ; *********************************************
         ; set up of start_time_step and end_time_step
         ; *********************************************

         t = f->time
         nt = dimsizes(t)
         
         ; *********************************************	
         ; 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(15))
                  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 ranges of x-, y- and z-coordinates
         ; ****************************************************         
                    
         xdim = dimsizes(data_all(0,0,0,:)) - 1    
         ydim = dimsizes(data_all(0,0,:,0)) - 1        
         zdim = dimsizes(data_all(0,:,0,0)) - 1   
    
         if ( .not. isvar("xe")) then     ; output x-coordinate range end (in index)
            xe = xdim
            if (parameter(53) .NE. "xdim") then
               if (stringtointeger(parameter(53)) .GT. xdim) then
                  print(" ")
                  print("range end for x-coordinate = "+parameter(53)+" is higher than available dimensions = "+xdim)
                  print(" ") 
                  exit
               end if
               if (stringtointeger(parameter(53)) .LT. 0 .OR. stringtointeger(parameter(53)) .LT. xs) then
                  print(" ")
                  print("range end for x-coordinate = "+parameter(53)+" is too small")
                  print(" ")
                  exit
               end if
               xe = stringtointeger(parameter(53)) 
            end if 
         else
           if (xe .GT. xdim) then
              print(" ")
              print("range end for x-coordinate = "+xe+" is higher than available dimensions = "+xdim)
              print(" ")
              exit
           end if
           if (xe .LT. 0 .OR. xe .LT. xs) then
              print(" ")
              print("range end for x-coordinate = "+xe+" is too small")
              print(" ")
              exit
           end if                  
         end if
         if ( .not. isvar("ye")) then     ; output y-coordinate range end (in index)
            ye = ydim
            if (parameter(57) .NE. "ydim") then
               if (stringtointeger(parameter(57)) .GT. ydim) then
                  print(" ")
                  print("range end for y-coordinate = "+parameter(57)+" is higher than available dimensions = "+ydim)
                  print(" ")
                  exit
               end if
               if (stringtointeger(parameter(57)) .LT. 0 .OR. stringtointeger(parameter(57)) .LT. ys) then
                  print(" ")
                  print("range end for y-coordinate = "+parameter(57)+" is too small")
                  print(" ")
                  exit
               end if
               ye = stringtointeger(parameter(57)) 
            end if 
         else
            if (ye .GT. ydim) then
               print(" ")
               print("range end for y-coordinate = "+ye+" is higher than available dimensions = "+ydim)
               print(" ")
               exit
            end if
            if (ye .LT. 0 .OR. ye .LT. ys) then
               print(" ")
               print("range end for y-coordinate = "+ye+" is too small")
               print(" ")
               exit
            end if
         end if
         if ( .not. isvar("ze")) then     ; output z-coordinate range end (in index)
            ze = zdim
            if (parameter(61) .NE. "zdim") then
               if (stringtointeger(parameter(61)) .GT. zdim) then
                  print(" ")
                  print("range end for z-coordinate = "+parameter(61)+" is higher than available dimensions = "+zdim)
                  print(" ")
                  exit
               end if
               if (stringtointeger(parameter(61)) .LT. 0 .OR. stringtointeger(parameter(61)) .LT. zs) then
                  print(" ")
                  print("range end for z-coordinate = "+parameter(61)+" is too small")
                  print(" ")
                  exit
               end if
               ze = stringtointeger(parameter(61)) 
            end if 
         else
            if (ze .GT. zdim) then
               print(" ")
               print("range end for z-coordinate = "+ze+" is higher than available dimensions = "+zdim)
               print(" ") 
               exit
            end if
            if (ze .LT. 0 .OR. ze .LT. zs) then
               print(" ")
               print("range end for z-coordinate = "+ze+" is too small")
               print(" ")
               exit
            end if
         end if
  
         ; Defining inner and outer loop depending on variable
         ; "sort" (values: time or layers) 
    
         if ( xyc .eq. 1 ) then
            lays = zs
            laye = ze
         end if
         if ( xzc .eq. 1 ) then
            lays = ys
            laye = ye
         end if
         if ( yzc .eq. 1) then
            lays = xs
            laye = xe
         end if

         if ( sort .eq. "time" ) then
            lis = start_time_step
            lie = end_time_step
            los = lays 
            loe = laye
         else
            lis = lays
            lie = laye
            los = start_time_step 
            loe = end_time_step
         end if 


         ; loops over time and layer (order depends on "sort")

         do lo = los, loe           				; lo: loop outer
            do li = lis, lie, no_lines*no_columns		; li: loop inner

               print("page "+page)
               page = page + 1

               ; Defining indices for reading data needed to plot one page
	       ; The indices depend on the parameter "sort"

               if ( sort .eq. "time" ) then
                  ta = li
                  tb = min( (/li + no_lines*no_columns - 1, end_time_step/) )
                  if ( xyc .eq. 1 ) then
                     xa = xs
                     xb = xe
                     ya = ys
                     yb = ye
                     za = lo
                     zb = lo
                  end if
                  if ( xzc .eq. 1 ) then
                     xa = xs
                     xb = xe
                     ya = lo
                     yb = lo
                     za = zs
                     zb = ze
                  end if
                  if ( yzc .eq. 1 ) then
                     xa = lo
                     xb = lo
                     ya = ys
                     yb = ye
                     za = zs
                     zb = ze
                  end if
               end if

               if ( sort .eq. "layer" ) then
                  ta = lo
                  tb = lo
                  if ( xyc .eq. 1 ) then
                     xa = xs 
                     xb = xe
                     ya = ys
                     yb = ye
                     za = li
                     zb = min( (/li + no_lines*no_columns - 1, ze/) )  
                  end if
                  if ( xzc .eq. 1 ) then
                     xa = xs 
                     xb = xe
                     ya = li
                     yb = min( (/li + no_lines*no_columns - 1, ye/) ) 
                     za = zs
                     zb = ze
                  end if
                  if ( yzc .eq. 1 ) then
                     xa = li 
                     xb = min( (/li + no_lines*no_columns - 1, xe/) ) 
                     ya = ys
                     yb = ye
                     za = zs
                     zb = ze
                  end if
               end if
                            
               data = data_all(ta:tb,za:zb,ya:yb,xa:xb)
       
               data!0 = "t"
               data!1 = "z"
               data!2 = "y"
               data!3 = "x"

               ; set up labels and plot

               ; xy cross section

               if ( xyc .eq. 1 ) then
                  cs_res@tiXAxisString   = "x [m]"
                  cs_res@tiYAxisString   = "y [m]"

                  if ( sort .eq. "time" ) then
                     if ( data&z(0) .eq. -1 ) then
                        level = "-average"
                     else 
                        level = "=" + data&z(0) + "m"
                     end if
                     do n = 0, tb-ta
                        cs_res@gsnCenterString = "t=" + data&t(n) + "s  z"+ level
                        plot(n) = gsn_csm_contour(wks_ps,data(n,0,:,:),cs_res)  
                     end do
                  end if

                  if ( sort .eq. "layer" ) then
                     do n = 0, zb-za
                        if ( data&z(n) .eq. -1 ) then
                           level = "-average"
                        else 
                           level = "=" + data&z(n) + "m"
                        end if
                        cs_res@gsnCenterString = "t=" + data&t(0) + "s  z"+ level
                        plot(n) = gsn_csm_contour(wks_ps,data(0,n,:,:),cs_res)
                     end do
                  end if
               end if

               ; xz cross section

               if ( xzc .eq. 1 ) then
                  cs_res@tiXAxisString   = "x [m]"
                  cs_res@tiYAxisString   = "z [m]"
  
                  if ( sort .eq. "time" ) then
                     if ( data&y(0) .eq. -1 ) then
                        level = "-average"
                     else 
                        level = "=" + data&y(0) + "m"
                     end if
                     do n = 0, tb-ta
                        cs_res@gsnCenterString = "t=" + data&t(n) + "s  z"+ level
                        plot(n) = gsn_csm_contour(wks_ps,data(n,:,0,:),cs_res)  
                     end do
                  end if

                  if ( sort .eq. "layer" ) then
                     do n = 0, yb-ya
                        if ( data&y(n) .eq. -1 ) then
                           level = "-average"
                        else 
                           level = "=" + data&y(n) + "m"
                        end if
                        cs_res@gsnCenterString = "t=" + data&t(0) + "s  z"+ level
                        plot(n) = gsn_csm_contour(wks_ps,data(0,:,n,:),cs_res)
                     end do
                  end if
               end if

               ; yz cross section

               if ( yzc .eq. 1 ) then
                  cs_res@tiXAxisString   = "y [m]"
                  cs_res@tiYAxisString   = "z [m]"
  
                  if ( sort .eq. "time" ) then
                     if ( data&x(0) .eq. -1 ) then
                        level = "-average"
                     else 
                        level = "=" + data&x(0) + "m"
                     end if
                     do n = 0, tb-ta
                        cs_res@gsnCenterString = "t=" + data&t(n) + "s  x"+ level
                        plot(n) = gsn_csm_contour(wks_ps,data(n,:,:,0),cs_res)  
                     end do
                  end if

                  if ( sort .eq. "layer" ) then
                     do n = 0, xb-xa
                        if ( data&x(n) .eq. -1 ) then
                           level = "-average"
                        else 
                           level = "=" + data&x(n) + "m"
                        end if
                        cs_res@gsnCenterString = "t=" + data&t(0) + "s  x"+ level
                        plot(n) = gsn_csm_contour(wks_ps,data(0,:,:,n),cs_res)
                     end do
                  end if
               end if

               delete(data)

               ; ***************************************************
               ; merge plots onto one page 
               ; ***************************************************
         
               gsn_panel(wks_ps, plot(0:n-1),(/no_lines,no_columns/),cs_resP)

               ; To hold no_lines and no_columns constant, it must be defined again:

               no_lines   = nl1
               no_columns = nc1

            end do 	
         end do  
      end if  
  
      delete(data_all)

   end do   

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