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



begin

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

; open file to write interpolated variable(s) to

out = addfile("foam_test.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 FOAM data file you will be comparing with the NCEP/NCAR data
; (need to use coordinate variables from this file to put into NCEP file)

foam = addfile("/grove3/iesuw/foam/ha.F1.0k.mone486605.nc","r")

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

; Enter dimensions of foam variables
; FOAM is on a Gaussian R15 grid with dims 40 lat by 48 lon

jlat = 40 
ilon = 48

; We will use the "f2gsh" function to regrid our FIXED-GRID NCEP data to the FOAM 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)


; Since the data is R15, we set twave=0
airt = f2gsh(airlat(:,:,:),(/jlat,ilon/),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 FOAM 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 FOAM data)

airt!0 = "time"
airt&time = foam->time
airt!1 = "lat"
airt&lat = foam->lat 
airt!2 = "lon"
airt&lon = foam->lon 

;write out data to file

out->airt = airt 
out->time = foam->time
out->lat = foam->lat 
out->lon = foam->lon 

end