;************************************************
; zonal_2.ncl
;
; Concepts illustrated:
;   - Drawing a zonal means plot
;   - Using dim_avg_Wrap to calculate a zonal average
;************************************************
;load "/usr/lib/ncarg/nclscripts/csm/gsn_code.ncl" 
;load "/usr/lib/ncarg/nclscripts/csm/gsn_csm.ncl" 
;load "/usr/lib/ncarg/nclscripts/csm/contributed.ncl" 
;load "/usr/lib/ncarg/nclscripts/contrib/time_axis_labels.ncl"
;************************************************
begin

;----------------------------------------------------------------------
; Converts YYYYMM time to "days since" time.
;----------------------------------------------------------------------
function yyyymm_to_time(time)
local ntime, stime, year, month, day, hour, minute, second, units
begin
  ntime  = dimsizes(time)
  stime  = tostring(time)
  year   = toint(str_get_cols(stime,0,3))
  month  = toint(str_get_cols(stime,4,5))
  day    = new(ntime,integer)
  hour   = new(ntime,integer)
  minute = new(ntime,integer)
  second = new(ntime,double)
 
;----------------------------------------------------------------------
; Note: you may need to change the day, hour, minute values depending
; on your definition of "YYYYMM".  Is this the first day and the 
; first hour of the month? Is it the middle of the month?
;----------------------------------------------------------------------
  day    = 15
  hour   = 0
  minute = 0
  second = 0.d
  units  = "days since " + min(year) + "-01-15 00:00"

  newtime = cd_inv_calendar(year,month,day,hour,minute,second,units,0)

  return(newtime)
end

time = (/200001,200002,200003,200004,200005,200006,200007,200008,\
           200009,200010,200011,200012/)

  time!0     = "time"
  time@units = "YYYYMM"


  
  

;************************************************
; variable and file handling
;************************************************
  files = systemfunc("ls MAM_GEOS_8y_rest.cam.h0.2001-*.nc")
;  files = systemfunc("more liste") 
  f = addfiles(files,"r")

  ListSetType (f, "cat")

begin
;When a variable is read from a series of files in "cat" mode, then the leftmost dimension of each variable will be COncatenated to form a new leftmost dimension whose size is the sum of the sizes of the leftmost dimension in each individual file.
;
  dsizes = getfiledimsizes(f[0])
  

   var = f[:]->NOX(:,:,{15:30},{65:105})	 ; select variable to ave

   u        = f[:]->U(:,:,{15:30},{65:105})
   v        = f[:]->V(:,:,{15:30},{65:105})

  date = f[:]->date
  ntim   = dimsizes(time)
  
  ;;************************************************
; read needed variables from file
;************************************************
   T = f[:]->T                                    ; select variable to ave
    p0   = 1000.
   hyam = f[:]->hyam                             ; get a coefficiants
   hybm = f[:]->hybm                              ; get b coefficiants
   PS   = f[:]->PS                                ; get pressure
   ;P0  = f[:]->lev
  ;p    = pres_hybrid_ccm (PS, p0, hyam, hybm)  ; p(ntim,klvl,nlat,mlon)
  ;p@long_name = "pressure at mid levels"
  ;p@units     =  PS@units
  ; print(P0)                               ; for illustration convert to hPa (mb)
  ;
;print(T)
;;************************************************
;************************************************
; define other arguments required by vinth2p
;************************************************
; type of interpolation: 1 = linear, 2 = log, 3 = loglog
   interp = 2 

; is extrapolation desired if data is outside the range of PS
   extrap = False

; create an array of desired pressure levels:
   lev = (/400.,300.,200.0,150.0,100.0,80.,60.,50. /)
   lev!0         = "lev"                  ; variable/dim name 
   lev&lev     =  lev                   ; create coordinate variable
   lev@long_name = "pressure"               ; attach some attributes
   lev@units     = "hPa"
   lev@positive  = "down"        
   time!0	= "month"
      
;************************************************
; calculate var on pressure levels
;************************************************
   TonP = vinth2p(var(:,:,:,:),hyam,hybm,lev,PS(:,{15:30},{65:105}),interp,p0,1,extrap)

   u_2 = vinth2p(u(:,:,:,:),hyam,hybm,lev,PS(:,{15:30},{65:105}),interp,p0,1,extrap)
   v_2 = vinth2p(v(:,:,:,:),hyam,hybm,lev,PS(:,{15:30},{65:105}),interp,p0,1,extrap)
;************************************************
printVarSummary(TonP)

TonP_avg =  dim_avg_n_Wrap(TonP, (/2,3/))  ; average time  
TonP_avg_reorder =  TonP_avg(lev|:,time|:) ; reordena las dimensiones para poder graficar en niveles de presion
u_avg =  dim_avg_n_Wrap(u_2, (/2,3/))  ; average time
v_avg =  dim_avg_n_Wrap(v_2, (/2,3/))  ; average time

u_avg_reorder =  u_avg(lev|:,time|:)
v_avg_reorder =  v_avg(lev|:,time|:)  ; average time
nlev = dimsizes(lev)
ntime = dimsizes(time)

printVarSummary(TonP_avg)
printVarSummary(u)
printVarSummary(v)
printVarSummary(TonP_avg_reorder)


;************************************************
; plot the graphic
;************************************************
 wks = gsn_open_wks("x11","NOX_vertical distribution_2001")
  wks1 = gsn_open_wks("x11","NOX_vertical_2001_wind")

  res                 = True                      ; plot mods desired
  res@gsnDraw        = True
  res@cnFillOn        = True                     ; turn on color fill
  ;res@cnFillPalette   = "wh-bl-gr-ye-re"          ; set color map
  res1 = True

  res1@gsnDraw  = True                          ; don't draw
  res1@gsnFrame = True                          ; don't advance frame
  res1@cnFillOn             = True               ; turn on color fill
  res1@cnFillPalette        = "gui_default"      ; set color map   
  res@cnLinesOn      = True              ; default is True
  ;res@cnLineLabelsOn = True             ; default is True, pone los valores sobre las lineas de el grafico
  ;res@cnLevelSelectionMode = "ManualLevels"       ; set manual contour levels
  res@tiYAxisString   = "Pressure (hPa)"
  res@cnFillPalette        = "gui_default"      ; set color map
  ;res@trYMaxF           =  400.               ; set max
  ;res@trYMinF            =  50.               ; set min

  ;res@cnLevelSelectionMode = "ExplicitLevels"
  ;res@cnLevels     = (/20,40,80, 120,160,200,240,280,320/)
  ;res@tmYUseLeft = False
  ;res@tmYRmode = "Manual"
  ;res@tmYRTickStartF = 1
  ;res@tmYRTickEndF = 6
  ;res@tmYRTickSpacingF = 1
  
  
;---Convert the time values and reassign as x's time coordinate array
 ; newtime = yyyymm_to_time(time)
 ; zavepp&time  := newtime
 ; TonP_avg_reorder&time = newtime
  
  ;u_avg_reorder&time = newtime
  ;v_avg_reorder&time = newtime
 ; printVarSummary(zavepp)
 
;;set vector;;
        res_vc                            = res
        res_vc@vcGlyphStyle               = "LineArrow"
        res_vc@vcLineArrowThicknessF      = 5
        res_vc@vcMinDistanceF             = 0.02
        
        res_vc@vcRefLengthF               = 0.03
  
;---Set resources necessary to nicely format X axis
 ; restick           = True
 ; restick@ttmFormat = "%c %Y"              ; Jan 2010, Feb 2010, etc.
 ; time_axis_labels(newtime,res,restick)
 ; time_axis_labels(newtime,res_vc,restick)
;***************************************
;   adjust label bar
;***************************************
;
; Retrieve contour levels.
;
  res@tiMainString = "NOX(ppt) avg (lat=15N:30N, lon=65E:105E)"
  res_vc@tiMainString = "NOX(ppt) avg (lat=15N:30N, lon=65E:105E)"
  plot = gsn_csm_pres_hgt (wks,TonP_avg_reorder,res)
  vector  = gsn_csm_pres_hgt_vector(wks1,TonP_avg_reorder,u_avg_reorder,v_avg_reorder,res_vc) 
 ;plot = gsn_csm_pres_hgt (wks,TonP_avg_reorder,res)
;overlay(plot,vector)
;***************************************
;   adjust label bar
;***************************************
;
; Retrieve contour levels.
  
getvalues plot@contour
    "cnLevels" : levels
  end getvalues

  ;res@lbLabelStrings       = sprintf("%4.2g",levels)   ; Format the labels
	
;plot = gsn_csm_pres_hgt (wks,zavepp,res)
;;drawing "map" will draw everything: map, contours, vectors, and text;;
;        draw(plot)
;        frame(wks1)

end
end


