;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; This script regrids NCEP/NCAR Reanalysis Data from a 
; GAUSSIAN GRID (T62, 94 by 192) 
; 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("lhflx_ncep_to_foam.nc","c")

; Open NCEP/NCAR Data File 
; This file contains a long-term mean of monthly latent heat flux on a gaussian grid

data = addfile("/grove4/sara/ncep/lhflx.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 "g2gsh" function to regrid our GAUSSIAN-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 g2gsh function expects data that has lat in ASCENDING order (i.e., -90 to 90)
; In the NCEP/NCAR files, lat is in DESCENDING order (i.e., 90 to -90)
; Therefore, we must re-order the lats so our data will be properly oriented

; Create a variable to hold the re-ordered data

lhflx = new((/12,94,192/),float)

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

k=93
do i = 0,93
	lhflx(:,i,:) = data->lhtfl(:,k,:)
	k = k-1
end do

; Now we regrid the data using the g2gsh function since our data is on a gaussian grid
; function g2gsh(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 lat 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
lh = g2gsh(lhflx(:,:,:),(/jlat,ilon/),0)


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

; Assign metadata to variables
; Since the g2gsh function does NOT create coordinate variables for our output variable "lh",
; 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)

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

;write out data to file

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

end