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
;; This script goes through the 3 output files of the SCM model
;; and creats 3 .ps files containing all the 1D variables (lev,time)
;; it finds (and are not constant).
;; Daniel Klocke 2012 Nov - created
;;--------------------------------------------------------------------------------
;; define some values:

plot_top = 150.		; top of plot in hPa
plot_bot = 1000.	; bottom of plot in hPa

plotp  = True 		; plot output feom progvar the names are hardcoded for now
plotd  = True		; plot output from diagvar
plotd2 = True		; plot output from diagvar2

;;--------------------------------------------------------------------------------
;; Some talking
print("Top pressure shown on plot will be "+plot_top+" hPa")
print("Bot pressure shown on plot will be "+plot_bot+" hPa")
print(" ")
print("Files where the script will look for 1D SCM output variables are:")
print("		progvar.nc  = " +plotp)
print("		diagvar.nc  = " +plotd)
print("		diagvar2.nc = " +plotd2)

;; INPUT files
  f  = addfile("progvar.nc","r")
  simID = str_split(f@simulationID," ")
  modID = str_split(f@modelID," ")

  if (plotd) then
    f1  = addfile("diagvar.nc","r")
  end if
  if (plotd2) then
    f2  = addfile("diagvar2.nc","r")
  end if
;; some necessary variables read from progvar  
  time   = f->time/86400
  nlev   = f->nlev
  nlevp1 = max(nlev)+1
  pres_f = f->pressure_f(0,:)/100.
  pres_h = f->pressure_h(0,:)/100.
;; creating an index for which level to plot between
;; (adds one level on either side for contouring and to be safe with the half/full levs)
  idxt = 0
  idxb = 0
  do i=0,dimsizes(nlev)-1
    if (plot_top .gt. pres_f(i)) then
      idxt = idxt + 1
    end if
    if (plot_bot .gt. pres_f(i)) then
       idxb = idxb + 1  
    end if
  end do
  idxb = idxb + 1
  idxt = idxt - 1
;; read variable names from files
  if (plotp) then
    var_names = getfilevarnames(f)
  end if

  if (plotd) then
    var_names1 = getfilevarnames(f1)
  end if

  if (plotd2) then
    var_names2 = getfilevarnames(f2)
  end if

;; Set plotting resources
  res                  = True         			; plot mods desired
  res@tmXTBorderOn     = False        			; top border off
  res@tmXBBorderOn     = False        			; bottom border off
  res@tmYLBorderOn     = False        			; left border off
  res@tmYRBorderOn     = False        			; right border off
  res@tmXTOn           = False        			; top tickmarks off
  res@tmYROn           = False        			; right tickmarks off
  res@tmXBMinorOn      = False        			; minor X tickmarks off
  res@tmYLMinorOn      = False        			; minor Y tickmarks off
  res@tmXBMajorThicknessF = 4.0
  res@tmYLMajorThicknessF = 4.0
  res@xyLineThicknessF = 3.0
  res@tiXAxisString    = "Time in days"			; set X-Axxis title

;;--- Plotting progvar variables
if (plotp) then
  print(" ")
  print("*-------- Plotting from PROGVAR ---------*")
  print("*----------------------------------------*")
  wks  = gsn_open_wks ("ps", ("progvar_1d_" +modID+"_" +simID))
  if (.not.all(ismissing(var_names))) then
    do i = 0,dimsizes(var_names)-1

      var = f->$var_names(i)$ 
      dsize = (dimsizes(var))      

      if(dimsizes(dimsizes(var)) .eq. 1 .and. (max(var) .ne. min(var)) .and. var_names(i) .ne. "time" .and. var_names(i) .ne. "nlev" .and. var_names(i) .ne. "nlevs" .and. var_names(i) .ne. "nlevp1" .and. var_names(i) .ne. "t_snow") then 
        res@gsnLeftString    =  var_names(i)
        res@tiYAxisString    = ""	
        plot = gsn_csm_xy(wks,time, var, res)
      end if

      delete(var)
      delete(dsize)

    end do
  end if
end if ; plotp

;;--- Plotting diagvar variables

if (plotd) then
  print(" ")
  print("*-------- Plotting from DIAGVAR ---------*")
  print("*----------------------------------------*")
  if (.not.all(ismissing(var_names1))) then
    wks  = gsn_open_wks ("ps", "diagvar_1d_" +modID+"_" +simID)
    do i = 0,dimsizes(var_names1)-1
      var = f1->$var_names1(i)$
      dsize = (dimsizes(var))
      if(dimsizes(dimsizes(var)) .eq. 1 .and. dsize(0) .eq. dimsizes(time) .and. max(var) .ne. min(var) .and. var_names1(i) .ne. "time") then 
        if (var_names1(i) .eq. "conv_type") then
          resxy = True
          resxy = res
          resxy@trYMinF          = 0 
          resxy@trYMaxF          = 3
          resxy@tmYLMode         = "Explicit"
          resxy@tmYLValues       = (/0,1,2,3/)
          resxy@tmYLLabels       = (/"None","Deep","Shallow","Midlevel"/)
          resxy@gsnLeftString    = ("Convection Type")
          resxy@xyMarkLineMode   = "Markers"
          resxy@xyMarker         = 16   
          resxy@xyMarkerColor    = "Red"
          resxy@tiYAxisString    = ""	
          plot = gsn_csm_xy(wks,time, var, resxy)
        else
          res@gsnLeftString    =  var_names1(i)
          res@tiYAxisString    = ""	
          plot = gsn_csm_xy(wks,time, var, res)
        end if
      end if 

      delete(var)
      delete(dsize)

    end do
  end if
end if ; plotd

;;--- Plotting diagvar variables

if (plotd2) then
  print(" ")
  print("*-------- Plotting from DIAGVAR2 --------*")
  print("*----------------------------------------*")
  if (.not.all(ismissing(var_names2))) then
    wks  = gsn_open_wks ("ps", "diagvar2_1d_" +modID+"_" +simID )
    do i = 0,dimsizes(var_names2)-1

      var = f2->$var_names2(i)$
      dsize = (dimsizes(var))

      if(dimsizes(dimsizes(var)) .eq. 1 .and. dsize(0) .eq. dimsizes(time) .and. max(var) .ne. min(var) .and. var_names2(i) .ne. "time") then 
        res@gsnLeftString    =  var_names2(i)
        res@tiYAxisString    = ""	
        plot = gsn_csm_xy(wks,time, var, res)
      end if 

      delete(var)
      delete(dsize)

    end do
  end if
end if ; plotd2

print(" ")
print("*----------------------------------------*")
print("Done plotting 1D variables found in files")

end