;----------------------------------------------------------------------
;-- Script:          unrotate_grid.ncl
;--
;-- Description:     REMO - unrotate rotated grid to original grid
;--                  Using functions from the library file 
;--                    functions_for_rotated_grids.ncl
;--
;-- 11.11.14 Karin Meier-Fleischer, DKRZ
;----------------------------------------------------------------------
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"

;------------------------------------------------------------
;-- set global constants
;------------------------------------------------------------
pi       =  4.0*atan(1.)
deg2rad  =  pi/180.
rad2deg  =  45./atan(1.)
fillval  = -99999.9

;------------------------------------------------------------
;-- Function:      create_lon2d(y,x)
;-- Description:   create 2d lon array from 1D rotlon, rotlat
;------------------------------------------------------------
undef("create_lon2d")
function create_lon2d(rotlat[*]:numeric, rotlon[*]:numeric)
local x, i
begin
  x = new((/dimsizes(rotlat),dimsizes(rotlon)/),typeof(rotlat))
;  print("Function create_lon2d: dimsizes of x: "+dimsizes(x))
  do i=0,dimsizes(rotlat)-1
    x(i,:) = rotlon
  end do
  return(x)
end

;------------------------------------------------------------
;-- Function:      create_lat2d(y,x)
;-- Description:   create 2d lat array from 1D rotlon,rotlat
;------------------------------------------------------------
undef("create_lat2d")
function create_lat2d(rotlat[*]:numeric, rotlon[*]:numeric)
local y, i
begin
  y = new((/dimsizes(rotlat),dimsizes(rotlon)/),typeof(rotlat))
;  print("Function create_lat2d: dimsizes of y: "+dimsizes(y))
  do i=0,dimsizes(rotlon)-1
    y(:,i) = rotlat
  end do
  return(y)
end

;------------------------------------------------------------
;-- Function:      (rotlat,rotlon,pollat,pollon)
;-- Description:   transform rotated longitude to longitude
;------------------------------------------------------------
undef("unrot_lon")
function unrot_lon( rotlat:numeric, rotlon:numeric, pollat[1]:numeric, pollon[1]:numeric )
local rotlat, rotlon, pollon, pollat, lon, s1, c1, s2, c2, rlo, rla, i, tmp1, tmp2
begin
  lon = fillval
  lon@_FillValue = fillval
  
  if (any(dimsizes(rotlat) .ne. dimsizes(rotlon)) .and. \
      (dimsizes(dimsizes(rotlat)).ne.1 .or. dimsizes(dimsizes(rotlon)).ne.1)) then
      print("Function unrot_lon: unrot_lon:  rotlat and rotlon dimensions do not match")
    return(lon)
  end if
  
  if (dimsizes(dimsizes(rotlat)).eq.1 .and. dimsizes(dimsizes(rotlon)).eq.1) then
    rla = create_lat2d(rotlat,rotlon)    ;-- create 2D latitude array
    rlo = create_lon2d(rotlat,rotlon)    ;-- create 2D longitude array
  else
    rla = rotlat
    rlo = rotlon
  end if
  
  rla = rla*deg2rad                      ;-- convert from degree to radians
  rlo = rlo*deg2rad                      ;-- convert from degree to radians
    
  lon := (/rlo/)                         ;-- reassign lon
  lon@_FillValue=fillval
  
  s1   = sin(pollat*deg2rad)
  c1   = cos(pollat*deg2rad)
  s2   = sin(pollon*deg2rad)
  c2   = cos(pollon*deg2rad)
  
  tmp1 = s2*(-s1*cos(rlo)*cos(rla)+c1*sin(rla))-c2*sin(rlo)*cos(rla)
  tmp2 = c2*(-s1*cos(rlo)*cos(rla)+c1*sin(rla))+s2*sin(rlo)*cos(rla)
  
  lon  = atan(tmp1/tmp2)*rad2deg
;  lon  = where(lon.lt.0, lon+360,lon)
;  lon  = where(lon.gt.180, lon-180,lon)
  
  lon@units = "degrees_east"
  print("Function unrot_lon: min/max     "+sprintf("%8.4f", min(lon(0,:)))+"  "+sprintf("%8.4f", max(lon(0,:))))
  
  delete([/rlo,rlo,c1,s1,c2,s2,tmp1,tmp2/])
  
  return(lon)
end

;------------------------------------------------------------
;-- Function:      unrot_lat(rotlat,rotlon,pollat,pollon)
;-- Description:   transform rotated latitude to latitude
;------------------------------------------------------------
undef("unrot_lat")
function unrot_lat( rotlat:numeric, rotlon:numeric, pollat[1]:numeric, pollon[1]:numeric )
local rotlat, rotlon, pollon, pollat, lat, s1, c1, rlo, rla, i
begin
  lat = fillval
  lat@_FillValue = fillval
  
  if (any(dimsizes(rotlat) .ne. dimsizes(rotlon)) .and. \
      (dimsizes(dimsizes(rotlat)).ne.1 .or. dimsizes(dimsizes(rotlon)).ne.1)) then
      print("Function unrot_lat:  rotlat and rotlon dimensions do not match")
    return(lat)
  end if
  
  if (dimsizes(dimsizes(rotlat)).eq.1 .and. dimsizes(dimsizes(rotlon)).eq.1) then
    rla = create_lat2d(rotlat,rotlon)    ;-- create 2D latitude array
    rlo = create_lon2d(rotlat,rotlon)    ;-- create 2D longitude array
  else
    rla = rotlat
    rlo = rotlon
  end if
  
  rla = rla*deg2rad                      ;-- convert from degree to radians
  rlo = rlo*deg2rad                      ;-- convert from degree to radians
  
  lat := (/rla/)                         ;-- reassign lat
  lat@_FillValue=fillval

  s1  = sin(pollat*deg2rad)
  c1  = cos(pollat*deg2rad)
  
  lat = s1*sin(rla)+c1*cos(rla)*cos(rlo)
  lat = asin(lat)*rad2deg
  
  lat@units = "degrees_north"
  print("Function unrot_lat: min/max     "+sprintf("%8.4f", min(lat(:,0)))+"  "+sprintf("%8.4f", max(lat(:,0))))
  
  delete([/rlo,rla,c1,s1/])
  
  return(lat)
end

;-- MAIN
begin
   f      =  addfile("tas_mon_1961-1990_rectilinear_grid_2D.nc","r")
   var    =  f->tas(0,0,:,:)                            ;-- read var to plot
   rlat   =  f->rlat                                    ;-- get 1D latitude array
   rlon   =  f->rlon                                    ;-- get 1D longitude array
   pollon =  f->rotated_pole@grid_north_pole_longitude  ;-- rotated pole longitude
   pollat =  f->rotated_pole@grid_north_pole_latitude   ;-- rotated pole latitude

   delete(var@units)
   delete(var@long_name)
   delete(var@standard_name)

;-- convert the 1D latitude and longitude arrays to 2D arrays
   lon2d  =  unrot_lon(rlat, rlon, pollat, pollon)
   lat2d  =  unrot_lat(rlat, rlon, pollat, pollon)

   nlat  = dimsizes(lat2d(:,0))
   nlon  = dimsizes(lon2d(0,:))


;-- open a workstation
   wks = gsn_open_wks("x11","plot_grids")

;-- set common plot resources for both plots
   res                       =  True                    ;-- set resources
   res@gsnAddCyclic          =  False                   ;-- data are not global, don't add lon cyclic point
   res@gsnMaximize           =  True                    ;-- maximize plot output
   res@gsnDraw               =  False                   ;-- will panel plots later
   res@gsnFrame              =  False
   
   res@cnFillOn              =  True                    ;-- turn on contour fill
   res@cnLinesOn             =  False                   ;-- turn off contour lines
   res@cnLevelSpacingF       =  1.125                   ;-- NCL chose 2.5

   res@pmTickMarkDisplayMode = "Always"                 ;-- draw tickmarks

   res@mpDataBaseVersion     = "MediumRes"              ;-- use finer database

   res_r  = res
   res_ur = res

;-- create the plot using the calculated 2D unrotated lat/lon arrays along with the native projection
   res_ur@mpCenterLonF          =  180. - pollon           ;-- set map center longitude
   res_ur@mpCenterLatF          =   90. - pollat           ;-- set map center latitude
   res_ur@tfDoNDCOverlay        = True
   res_ur@mpLimitMode           = "Corners"
   res_ur@mpLeftCornerLatF      = lat2d(0,0)
   res_ur@mpLeftCornerLonF      = lon2d(0,0)
   res_ur@mpRightCornerLatF     = lat2d(nlat-1,nlon-1)      
   res_ur@mpRightCornerLonF     = lon2d(nlat-1,nlon-1)      
   res_ur@tiMainString          = "calculated 2D lat/lon arrays"

   plot_ur = gsn_csm_contour_map(wks,(/var/),res_ur)
   
;-- create the plot using the original rotated lat/lon coordinate arrays
   var&rlat@units = "degrees_north"
   var&rlon@units = "degrees_east"

   res_r@mpMinLatF             =  min(var&rlat)
   res_r@mpMaxLatF             =  max(var&rlat)
   res_r@mpMinLonF             =  min(var&rlon)
   res_r@mpMaxLonF             =  max(var&rlon)
   res_r@tiMainString          = "rlat/rlon coordinate arrays"

   plot_r = gsn_csm_contour_map(wks,var,res_r)

;-- panel both plots
   gsn_panel(wks,(/plot_r,plot_ur/),(/1,2/),False)
end