;***************************************************************
; godas_3.ncl
;***************************************************************
;
; Concepts illustrated:
;   - Reading multiple GODAS netCDF files
;   - Calculating a monthly climatology using user specified years
;   - Calculating global anomalies at each grid point   
;   - Calculating a regional (area averaged) anomaly time series
;   - Drawing a time series plot
;
;***************************************************************
;
; These files are loaded by default in NCL V6.2.0 and newer
; 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 the variable which spans 30 netCDF files [1980-2009]
; Use 'short2flt' to unpack the variable which is type 'short'
;***************************************************************
  diri = "./"
  fils = systemfunc ("cd "+diri+" ; ls pottmp.*.nc") ; file paths
  fils = systemfunc ("cd "+diri+" ; ls pottmp.*.nc") ; file paths
  f    = addfiles (fils, "r")   
  print(fils)

  x    = short2flt( f[:]->pottmp(:,{5},:,:) ) ; read/unpack from all yearly files
  printVarSummary (x)             ; [time| 354]x[lat| 418]x[lon| 360]

;***************************************************************
; Create a 1995-2004 monthly climatology
; Use cd_calendar to return 'human understandable' units      
;***************************************************************
  yrStrt = 1995
  yrLast = 2004

  time   = x&time                   ; time associated with "x"
  yyyymm = cd_calendar(time, -1)    ; units 'yyyymm' 
  yyyy   = yyyymm/100               ; truncate to get yyyy only
  ii     = ind(yyyy.ge.yrStrt .and. yyyy.le.yrLast)

  xClm   = clmMonTLL( x(ii,:,:) )
  xClm@time_span = yrStrt+"-"+yrLast
  printVarSummary(xClm)
  printMinMax(xClm, True)
   
;***************************************************************
; At each grid point calculate anomalies from 'xClm'
;***************************************************************
  jj     = ind(yyyy.le.2008)
  xAnom  = calcMonAnomTLL ( x(jj,:,:), xClm)

;***************************************************************
; "Trick" overwrite the "days since ..." with yyyymm
;***************************************************************
  xAnom&time = yyyymm(jj)

  printVarSummary(xAnom)          ; note time units
  printMinMax(xAnom, True)

;***************************************************************
; Specify a region and compute an areal average anomaly
; Use simple cosine weighting
;***************************************************************
  latS   = 10
  latN   = 20
  lonL   = 280
  lonR   = 360

  clat   = f[0]->lat
  clat   = cos(0.01745329*clat)
  
  xAnomRegion = wgt_areaave_Wrap(xAnom(:,{latS:latN},{lonL:lonR}) \
                                ,clat({latS:latN}), 1.0, 0)
  printVarSummary(xAnomRegion)
  printMinMax(xAnomRegion, True)

;********************************
; create map plots
;********************************    
  nmo = 0                                 ; for demo, only plot Jan
  wks   = gsn_open_wks ("png", "godas"  ) ; send graphics to PNG file

  res   = True                                     ; plot mods desired
  res@gsnMaximize         = True
;;res@gsnPaperOrientation = "portrait"             ; force portrait

  res@cnFillOn            = True                   ; turn on color
  res@cnFillPalette       = "amwg"                 ; set color map
  res@cnLinesOn           = False
  res@cnLineLabelsOn      = False
  res@mpLandFillColor     = "grey"                 ; color of land
  res@mpFillDrawOrder     = "PostDraw"

  res@tiMainString         = "January Clm: "+yrStrt+"-"+yrLast
  plot = gsn_csm_contour_map(wks,xClm(nmo,:,:), res)  

;********************************
; create time series plot 
;********************************    
  yrfrac         = yyyymm_to_yyyyfrac(yyyymm(jj), 0.0)  ; plot axis

  resxy          = True                   ; plot mods desired
  resxy@gsnMaximize= True                 ; force ps/eps/pdf as large as possible
  resxy@vpHeightF= 0.4                    ; change aspect ratio of plot
  resxy@vpWidthF = 0.75                 
  
 ;resxy@trYMinF  = -1.0                   ; min value on y-axis
 ;resxy@trYMaxF  =  1.0                   ; max value on y-axis
  resxy@trXMinF  = 1980
  resxy@trXMaxF  = 2010                   ; Force to 2012          

  resxy@gsnLeftString   = "pottmp: Region 10N-20N  , 80W-0" 
  resxy@gsnRightString  =  x@units
  resxy@gsnYRefLine     = 0.0             ; create a reference line   

  resxy@gsnAboveYRefLineColor = "red"              ; above ref line fill red
  resxy@gsnBelowYRefLineColor = "blue"             ; below ref line fill blue

  plot = gsn_csm_xy (wks,yrfrac,xAnomRegion,resxy) ; create plot

end