;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; This script regrids NCEP/NCAR Reanalysis Data from a
; FIXED GRID (2.5 degrees by 2.5 degrees, 73 by 144)
; to the CCM3 T42 GAUSSIAN GRID (64 by 128) and outputs it to
; a netCDF file that can be read with GrADS
;
; surface data example
;
; questions? Email Sara Rauscher at saraamy@hotmail.com
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; OPEN DATA FILES
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

begin

; open file to write interpolated variable(s) to
out = addfile("ccm3.nc","c")

; Open NCEP/NCAR Data File
; This file contains a long-term mean of monthly air temperature on a fixed grid

data = addfile("/grove4/selin/ncep/air.mon.ltm.nc", "r")

; Open CCM3  data file you will be comparing with the NCEP/NCAR data
; (need to use coordinate variables from this file to put into NCEP file)

ccm3 = addfile("/grove4/iesuw/ccm3_jja/CCM3-0KPMIPEX-MONE1317.modnoshift.nc","r")

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; START the REGRIDDING PROCESS
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

; Enter dimensions of foam variables
; CCM3 is on a Gaussian T42 grid with dims 64 lat by 128 lon

; set number of lats and lons for output file
jlatitude = 64
ilongitude = 128

; We will use the "f2gsh" function to regrid our FIXED-GRID NCEP data to the CCM3 Gaussian Grid
; Information on this function is available at
; http://ngwww.ucar.edu/ngdoc/ng/ref/ncl/functions/g2g.html

; The f2gsh function expects data that has latitude in ASCENDING order (i.e., -90 to 90)
; In the NCEP/NCAR files, latitude is in DESCENDING order (i.e., 90 to -90)
; Therefore, we must re-order the latitudes so our data will be properly oriented

; Create a variable to hold the re-ordered data

airlat = new((/12,73,144/),float)

; Using a do loop, re-order the data (from 0...n-1 to n-1...0)
; There are 73 latitudes in the NCEP file; since NCL indexes from 0...n-1,
; our do loop goes from 0 to 72

k=72
do i = 0,72
        airlat(:,i,:) = data->air(:,k,:)
        k = k-1
end do

; Now we regrid the data using the f2gsh function since our data is on a fixed grid
; function f2gsh(grid : numeric, outdims[2] : integer, twave[1] : integer)
; Arguments:
; grid: input array of 2 or more dimensions (last two dimensions must be nlata x nlona, values must be in ascending latitude order)
; outdims : indicates dimensions of rightmost two dimensions of output grid (outdims[0] = nlatb, outdims[1] = nlonb)
; twave : scalar integer indicating the optional wave truncation
;         twave = 0 => exact interpolation
;         twave > 0 => exact interpolation and triangular truncation at twave
;         twave < 0 => exact interpolation, triangular truncation at twave and spectral coefficient tapering (the
;         effect is to smooth the data)


airt = f2gsh(airlat(:,:,:),(/jlatitude,ilongitude/), 0)

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;Write out data to file
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

; Assign metadata to variables
; Since the f2gsh function does NOT create coordinate variables for our output variable "airt",
; we must assign the coordinate variables from our foam file and copy them into our output file

; Assign coordinate variables to our new regridded variable
; Note that we put the variable "time" from our CCM3 file into our new output file
; so that grads can difference the two files (otherwise, if we kept the date
; from the original NCEP file, we would have to use xdfopen to compare this new
; output file and our CCM3 data)

airt!0 = "time"
airt&time = ccm3->time
airt!1 = "latitude"
airt&latitude = ccm3->latitude 
airt!2 = "longitude"
airt&longitude = ccm3->longitude 

;write out data to file

out->airt =airt 
out->time = ccm3->time
out->latitude = ccm3->latitude 
out->longitude = ccm3->longitude 

end