;*************************************************
; profile_master.ncl
; zp/height vs. x plot. It includes
; the local time axis, the special hao labeling, etc.
;************************************************
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/shea_util.ncl" 
load "$HAO_NCL_LIB/hao_contrib.ncl"
load "$HAO_NCL_LIB/halfzp2full.ncl"
load "$HAO_NCL_LIB/labeling.ncl"
load "$HAO_NCL_LIB/axis.ncl"
load "$HAO_NCL_LIB/derive.ncl"
;************************************************
begin
;************************************************
; user parameters

; fields must be 4D
;************************************************
 ncfile   = "/chewy/d/murphys/files/astrid.nc"

 fields   = (/"ED1M3D","TN"/) 
 mtime    = (/40,22,0/)      ; select time
 
 
 type     = "slt"           ; "lon","slt","zmean","gmean"
 lon2plt  = 90              ; lon to plot
 slt      = 12              ; slt to plot
 
 lat2plt  = -22
 haxis    = True
 
; interp y-axis to height
 height = False
 top    = 100
 bott   = 40
 dz     = 30
 
; limit axis
;  ymin     = -13
;  ymax     = -5
;  xmin     = -13
;  xmax     = -5
;  xlog     = True
 
; convert units
;  density_units = "CM3"  ; "MMR","CM3","CM3-MR","GM/CM3"
 
 plt_type = "ps"            ; ps,eps,epsi,x11,or ncgm
 plt_name = "profile_master"  ; name of x11 window or output plot
;************************************************
; no user changes after this point
;************************************************
  in = addfile(ncfile,"r")
; read in these variables for the hao labeling  
  ut = in->ut
  
; get time index from model time  
  t2plt = model_times2idx(mtime,in)           

 ; convert slt to longitude if necessary  
  if(type.eq."slt")then
    arg     = True
    arg@slt = slt
    if(isvar("lon2plt"))then
      delete(lon2plt)
    end if
    lon2plt = slt2lon(ut(t2plt),slt)
  end if 
  
; plot parameters

  wks = gsn_open_wks(plt_type,plt_name)      ; open a ps file
  print("creating "+plt_name+"."+plt_type)     ; so can cut & paste fm window

  res                          = True        ; plot mods desired
  res@xyComputeXMin            = True
  res@gsnDraw                  = False       ; don't draw yet
  res@gsnFrame                 = False       ; don't advance frame yet
  res@tiXAxisOn                = False
  res@tmLabelAutoStride        = True

  if(isvar("colormap"))then
    gsn_define_colormap(wks,colormap)             ; choose colormap
  end if  

  if(isvar("xlog"))then
    res@trXLog = True
  end if
  if(isvar("ymin"))then
    res@trYMinF = ymin
  end if
  if(isvar("ymax"))then
    res@trYMaxF = ymax
  end if
  
  if(isvar("xmin"))then
    res@trXMinF = xmin
  end if
  if(isvar("xmax"))then
    res@trXMaxF = xmax
  end if
  
  
; main loop  
  do i = 0,dimsizes(fields)-1
    if(isfilevar(in,fields(i)))then
      if(isfilevardim(in,fields(i),"mlat"))then  ; magnetic data
        x      = mag2latlon(in->$fields(i)$)
        y      = in->mlev
        height = False
        haxis  = False
      else
        x   = in->$fields(i)$
        y      = in->lev
        height = False
        if(.not.isatt(in,"tgcmproc"))then
          x   = halfzp(fields(i),x)                ; interp to full zp 
        end if        
      end if

; denisty conversion if necessary
      if(isvar("density_units"))then
        x@name = fields(i)
        x = den_convert_units(x,in,density_units)
      end if

; check to see if 4D
      if(dimsizes(dimsizes(x)).ne.4)then
        print("Variable "+fields(i)+" does not have 4D, skipping")
        break
      end if

; assign _FillValue
      if(isatt(in,"missing_value"))then
        x@_FillValue=in@missing_value
      end if

; prepare top text. this must be before height interpolation  
      arg      = True
      arg@time = t2plt
      arg@mag  = isMag(x)
      arg@file = in
      if(type.eq."lon".or.type.eq."slt")then
        arg@lon   = lon2plt
        arg@lat   = lat2plt
      end if
      if(isatt(x,"differences"))then
        arg@diff = True
        haxis    = False
        height   = False
      end if

; set height and haxis to false if already interpolated  
      if(x&lev@units.eq."km")then
        if(height)then
          res@trYMinF = bott
          res@trYMaxF = top
        end if
        haxis       = False
        height      = False
      end if

; interpolate to height if necessary
     
      if(height)then 
        if(bott .gt. top )then
          print("bott must be less than top when specifying heights")
          height = False
        else
          arg      = True
          arg@bott = bott
          arg@top  = top
          arg@dz   = dz
          haxis    = False
          z        = in->Z
          tmp      = zplev_to_height(z,x,arg)
          delete(x)
          x        = tmp
          delete(tmp)

          delete(y)            ; remove org lev array
          y = x&lev            ; reasign new interpolated level array
        end if
      end if

; calc zonal mean if necessary
      if(type.eq."zmean")then
        zx = x(:,:,:,0)                        ; trick to retain coordinates
        zx = dim_avg(x)                        ; calc zonal average
        delete(x)
        x = zx
        delete(zx)
        arg@zmean = True
        arg@lat   = lat2plt
      end if

; calc global mean if necessary
      if(type.eq."gmean")then
        tmp = x(:,:,:,0)               ; trick to retain coordinates
        gx  = tmp(:,:,0)
        tmp = dim_avg(x)          ; calc global average
        gx = dim_avg(tmp)
        delete(x)
        x = gx
        delete(gx)
        delete(tmp)
        arg@gmean = True
      end if

; add attribute to x so text_bott knows it is a profile
      x@profile = True

; get title depending on availability of info
      res   = set_title(x,fields(i),res)         
;*********************************
; create plot
;*********************************
; global mean
      if(dimsizes(dimsizes(x)).eq.2)then
        if(isvar("haxis").and. haxis.eq.True)then
          z    = in->Z
          tmp  = z(:,:,:,0)
          zx   = tmp(:,:,0)
          tmp  = dim_avg(z)
          zx   = dim_avg(tmp)
          res  = h_axis(zx(t2plt,:),res)
        end if
        res = textTop_prof(wks,arg,res)
        plot = gsn_csm_xy(wks,x(t2plt,:),y,res)
        textBott(wks,plot,x(t2plt,:),t2plt,ncfile,in,res)
        draw(plot)
        frame(wks)
      end if
; zonal mean
      if(dimsizes(dimsizes(x)).eq.3)then  
        if(isvar("haxis").and. haxis.eq.True)then
          z    = in->Z
          dimz = z(:,:,:,0)
          dimz = dim_avg(z)
          res  = h_axis(dimz(t2plt,:,{lat2plt}),res)
        end if
        res = textTop_prof(wks,arg,res)
        plot = gsn_csm_xy(wks,x(t2plt,:,{lat2plt}),y, res)
        textBott(wks,plot,x(t2plt,:,{lat2plt}),t2plt,ncfile,in,res)
        draw(plot)
        frame(wks)
      end if
; full variable
      if(dimsizes(dimsizes(x)).eq.4)then 
        if(isvar("haxis").and.haxis.eq.True)then
          z    = in->Z
          res  = h_axis(z(t2plt,:,{lat2plt},{lon2plt}),res)
        end if
        res  = textTop_prof(wks,arg,res)
        plot = gsn_csm_xy(wks,x(t2plt,:,{lat2plt},{lon2plt}),y, res)
        textBott(wks,plot,x(t2plt,:,:,{lon2plt}),t2plt,ncfile,in,res)
        draw(plot)
        frame(wks)
      end if
      delete(x)
      delete(y)
      delete(arg)
    else
      print("Variable "+fields(i)+" does not exist on file, skipping")
    end if
  end do
end