regrid a rectilinear grid to a POP ocean grid.

regrid_standard_to_POP.ncl

;----------------------------------------------------------------------
; ESMF_regrid_28.ncl
;
; Author: Riley X. Brady
; Date: February 5th, 2019
; *This script was adapted from the following example:
; http://www.ncl.ucar.edu/Applications/Scripts/ESMF_regrid_28.ncl
;
; - takes rectilinear global grid (e.g. 180x360) and regrids to POP 1 degree.
; - plots a single time slice of both to make sure the comparison looks good.
;
; Opposite approach (POP to world 1 deg):
; http://www.ncl.ucar.edu/Applications/Scripts/ESMF_regrid_24.ncl
;
; Unstructured to POP:
; http://www.ncl.ucar.edu/Applications/Scripts/ESMF_regrid_19.ncl
;----------------------------------------------------------------------
load "$NCARG_ROOT/lib/ncarg/nclscripts/esmf/ESMF_regridding.ncl"

begin
;;---Config options
remap_method = "bilinear" ; "bilinear", "patch", "conserve"
var_name     = "sDIC"
f_in         = "../JMA.oceanCO2.199001-201712_-180to180.nc"
f_out        = "sDIC.nc"


;---Specify remap method
;    remap_method = "bilinear"     ; "bilinear" , "patch", "conserve"

;---Specify name of rectilinear (RL) grid
    RL_res       = "1degGlobal"

;---Specify name of destination POP grid
    DstGridName  = "gx1v6"
    POP_res      = DstGridName

;---Specify name of weight file to be generated; name of destination grid; dst directory
    WgtFileDir   = "../weight_files/"
    WgtFileName  = "wgt.RL_"+RL_res+"_to_"+DstGridName+"."+remap_method+".nc"

;---Data file containing source grid lat and lon
;    f_in = "../JMA.oceanCO2.199001-201712_-180to180.nc"
    sfile    = addfile(f_in, "r")

;---Data file containing destination POP grid information
    POP_dir  = "../grids/"
    POP_file = "gx1v6_ocn.nc"
    dfile    = addfile(POP_dir+POP_file, "r")

;---Get variable to regrid from source file (rectilinear)
;    var_name = "ALK"                                  ; type short must convert to float or double
    var_in   = sfile->$var_name$(:,:,:) ; just need grid
    src_lat  = sfile->lat
    src_lon  = sfile->lon

;---Sample plot options
    pltDir   = "../plots/"
    pltType  = "png"                  ; send graphics to PNG file
    pltName  = var_name + "_RL_" + RL_res + "_to_" + DstGridName

;---Set up regridding options
    Opt                = True
    Opt@InterpMethod   = remap_method

    Opt@SrcGridLat     = src_lat
    Opt@SrcGridLon     = src_lon

    Opt@DstGridLat     = dfile->TLAT                   ; scalar grid locations
    Opt@DstGridLon     = dfile->TLONG
    Opt@DstGridType    = "curvilinear"

    Opt@RemovePETLog   = True
    Opt@RemoveSrcFile  = True ; only preserve weight file
    Opt@RemoveDstFile  = True ; only preserve weight file

    POP_var            = dfile->TEMP(0,0,:,:)          ; read POP variable to create mask
                                                       ; first time and 'near' sfc level
    Opt@DstMask2D      = where(ismissing(POP_var),0,1)

    Opt@ForceOverwrite = True
    Opt@PrintTimings   = True

    wgtfile            = "../weight_files/" + WgtFileName
    Opt@WgtFileName    = wgtfile
;---Call the regridding function; Generate weight file
    var_regrid = ESMF_regrid(var_in,Opt)
    printVarSummary(var_regrid)

    dims_regrid = dimsizes(var_regrid)
    rank_regrid = dimsizes(dims_regrid)
    var_regrid!(rank_regrid-1) = "ncol"

    delete(var_regrid@lon1d)
    delete(var_regrid@lat1d)
    var_regrid!1 = "nlat"
    var_regrid!2 = "nlon"
    printVarSummary(var_regrid)

;----Write out regridded file
; https://www.ncl.ucar.edu/Applications/method_2.shtml
    ; output file name
;    f_out = "ALK.nc"
    system("/bin/rm -f " + f_out)    ; remove if exists
    fout = addfile(f_out, "c")

    ; set up netCDF file
    setfileoption(fout, "DefineMode", True)

    ; global attributes
    fAtt = True
    fAtt@creation_date = systemfunc("date")
    fAtt@weight_file = WgtFileName
    fAtt@src_file = f_in
    fAtt@author = "Riley X. Brady"
    fAtt@contact = "riley.brady@colorado.edu"
    fileattdef(fout, fAtt)

    ; set up dimensions
    dimNames = (/"time", "nlat","nlon"/)
    dimSizes = (/-1, 384, 320/)
    dimUnlim = (/True, False, False/)
    filedimdef(fout,dimNames,dimSizes,dimUnlim)

    ; create variables
    TLAT = dfile->TLAT
    TLONG = dfile->TLONG
    TAREA = dfile->TAREA
    filevardef(fout,var_name,typeof(var_regrid),getvardims(var_regrid))
    filevardef(fout,"TLAT",typeof(TLAT),getvardims(TLAT))
    filevardef(fout,"TLONG",typeof(TLONG),getvardims(TLONG))
    filevardef(fout,"TAREA",typeof(TAREA),getvardims(TAREA))

    ; copy attributes from variables
    filevarattdef(fout,var_name,var_in)
    remapAtt=0
    remapAtt@remap_method = remap_method
    filevarattdef(fout,var_name,remapAtt)
    filevarattdef(fout,"TLAT",TLAT)
    filevarattdef(fout,"TLONG",TLONG)
    filevarattdef(fout,"TAREA",TAREA)
    setfileoption(fout,"DefineMode",False)

    ; add to netCDF file
    fout->$var_name$ = (/var_regrid/)
    fout->TLAT = (/TLAT/)
    fout->TLONG = (/TLONG/)
    fout->TAREA = (/TAREA/)

;----------------------------------------------------------------------
; Plotting section
;
; This section creates filled contour plots of both the original
; data and the regridded data, and panels them.
;----------------------------------------------------------------------
   ; just need one time slice for plotting
    orig_grid = var_in(0, :, :)

    dims_in = dimsizes(orig_grid)
    rank_in = dimsizes(dims_in)

    wks = gsn_open_wks(pltType,pltDir+pltName)

;---Resources to share between both plots
    res                      = True              ; Plot modes desired.

    res@gsnDraw              = False
    res@gsnFrame             = False

    res@gsnMaximize          = True              ; Maximize plot

    res@cnFillOn             = True              ; color plot desired
    res@cnLinesOn            = False             ; turn off contour lines
    res@cnLineLabelsOn       = False             ; turn off contour labels
    res@cnFillMode           = "RasterFill"      ; turn raster on

    res@lbLabelBarOn         = False       ; Will turn on in panel later

    res@mpFillOn             = False

   ;res@gsnLeftString        = var_name    ; long_name is too long!
    res@gsnLeftString        = orig_grid@long_name

;---Resources for plotting regridded data
    res@gsnAddCyclic  = True               ; default is True

    dims = tostring(dimsizes(orig_grid))
    rank = dimsizes(dims)
    res@tiMainString = "Original rectilinear grid: " + str_join(dims(rank-2:)," x ")

    plot_orig = gsn_csm_contour_map(wks,orig_grid,res)

;---Resources for contouring the above grid; apply to next grid

    res@cnLevelSelectionMode = "ManualLevels"
    getvalues plot_orig@contour
      "cnMinLevelValF"       : res@cnMinLevelValF
      "cnMaxLevelValF"       : res@cnMaxLevelValF
      "cnMaxLevelValF"       : res@cnMaxLevelValF
      "cnLevelSpacingF"      : res@cnLevelSpacingF
    end getvalues

;---Resources for plotting destination (POP) data
    ; just one slice for plotting
    regrid = var_regrid(0,:,:)
    regrid@lat2d = dfile->TLAT
    regrid@lon2d = dfile->TLONG
    res@gsnAddCyclic = True                 ; default is True

    res@tiMainString  = "Rectilinear to POP grid using "+Opt@InterpMethod + " method"

    plot_regrid = gsn_csm_contour_map(wks,regrid,res)

;---Draw both plots in a panel
    pres                   = True
    pres@gsnMaximize       = True
    pres@gsnPanelLabelBar  = True

    gsn_panel(wks,(/plot_orig,plot_regrid/),(/2,1/),pres)
end