;---------------------------------------------------------------
;-- DKRZ NCL Light Hands-On                   EGU 2014   SPM1.12
;--
;-- Example script:       ex_regrid_1.ncl    
;-- 
;-- Settings:             regrid data on a rotated grid to a 
;--                       rectilinear grid,
;--                       filled contour plot
;-- 
;-- 2014-03-22   Karin Meier-Fleischer (meier-fleischer@dkrz.de)
;---------------------------------------------------------------
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"
;load "$NCARG_ROOT/lib/ncarg/nclscripts/esmf/ESMF_regridding.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/contrib/calendar_decode2.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/crop.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"

begin
;-- open file
 a = addfile("E7_tasmax2000.nc","r")             ;read the nc file
 b = addfile("E7_tasmin2000.nc","r")
 c = addfile("$NCARG_ROOT/lib/ncarg/data/cdf/landsea.nc","r")   ; read netCDF for landmask

  lsdata = c->LSMASK                                    ; 1x1 degree of baseline land/sea mask
  tasmax = a->tasmax                                   ; read daily tasmax values from its nc file
  tasmin = b->tasmin                                   ; read daily tasmin values from its nc file
 
  time = a->time                                       ; read time values
  time_bnds = a->time_bnds
  lat = a->lat
  lon = a->lon 
  lon_bnds = a->lon_bnds
  lat_bnds = a->lat_bnds
  rlon = a->rlon
  rlat = a->rlat
  rotated = a->rotated_latitude_longitude

  lsm = landsea_mask(lsdata,lat,lon)

  tmax = tasmax + 273.15                               ; convert K to Celsius degree
  tmin = tasmin + 273.15                               ; convert K to Celsius degree

  time2 = calendar_decode2(time,-5)                    ; convert the time codes (days since 1949,12,1) to yyyy,mm,dd in integer
  doy = day_of_year(time2(:,0),time2(:,1),time2(:,2))  ; calculate day of year

  radext = radext_fao56(doy, lat, 0)             ; Compute extraterrestrial radiation for daily periods
  radext = where(ismissing(radext), 0, radext)   ; convert the missing value due to limited validity at high latitudes to 0
  evtH_0 = refevt_hargreaves_fao56( tmin, tmax, radext, (/0,0,0/) )         ; Hargreaves ETo equation in mm/day
    out1 = addfile("evtH_31.nc","c")
    out1->evtH = evtH_0
    out1->evtH!0 = "time"
    out1->evtH!1 = "rlat"
    out1->evtH!2 = "rlon"
    out1->lat = lat
    out1->lat_bnds = lat_bnds
    out1->lon = lon
    out1->lon_bnds = lon_bnds
    out1->rlon = rlon
    out1->rlat = rlat
    out1->rotated_latitude_longitude = rotated
    out1->time = time
    out1->time_bnds = time_bnds
    
    
   evtH_m = dim_avg_n_Wrap(evtH_0,0)

  lat2d  = lat                 ;-- 2D latitude array
  lon2d  = lon                 ;-- 2D longitude array
  minlat = min(lat2d)             ;-- retrieve minimum latitude value
  minlon = min(lon2d)             ;-- retrieve maximum latitude value
  maxlat = max(lat2d)             ;-- retrieve minimum longitude value
  maxlon = max(lon2d)             ;-- retrieve maximum longitude value
  

  var = evtH_m(:,:)
  var@lat2d = lat2d               ;-- set 2D latitudes
  var@lon2d = lon2d               ;-- set 2D longitudes
  var = mask(var,lsm.eq.0,False)  ;-- set the values in land 


;-- begin the graphics

  wks  = gsn_open_wks("pdf", "PET")             ;-- open workstation

  plot = new(1,graphic)                                 ;-- assign plot array for the two plots

   cmap = read_colormap_file("BlGrYeOrReVi200")

;-- global resources 
  res                        =  True
  res@gsnDraw                =  False                   ;-- don't draw the plot yet
  res@gsnFrame               =  False                   ;-- don't advance the frame yet
  res@gsnAddCyclic           =  False                   ;-- don't add cyclic point
  res@gsnMaximize            =  True
  res@cnFillOn               =  True                    ;-- contour fill
  res@cnFillPalette          =  cmap                    ;-- choose colormap
  res@cnLinesOn              =  False                   ;-- don't draw contour lines
  res@cnLevelSelectionMode   = "ManualLevels"           ;-- use same values for both plots
  res@cnMinLevelValF         =  toint(min(var))         ;-- minimum contour level
  res@cnMaxLevelValF         =  toint(max(var))         ;-- maximum contour level
  res@cnLevelSpacingF        =  1                       ;-- contour level spacing

  
  res@mpMinLatF              =  minlat              ;-- minimum lat
  res@mpMaxLatF              =  maxlat              ;-- maximum lat
  res@mpMinLonF              =  minlon              ;-- minimum lon
  res@mpMaxLonF              =  maxlon              ;-- maximum lon
  
  res@sfXArray               =  lon2d                   ;-- x-axis location of grid cell center
  res@sfYArray               =  lat2d                   ;-- y-axis location of grid cell center
;---aridity_label
  res@lbFillColors           = cmap
  res@lbLabelFontHeightF     = 0.007
  res@lbLabelAlignment       = "BoxCenters"
;  res@lbBoxCount             = 5

;-- create the plots
  plot = gsn_csm_contour_map(wks,var,res)            ;-- create plot on original curvilinear grid
  
  draw(plot)
  frame(wks)  
 
  
end