; Example NCL code to read CoRTADv3 via OPeNDAP and
; convert to netCDF.
; Provided to NODC by Ruben van Hooidonk, NOAA AOML
; February 2011
; Please contact Ruben with any questions about the code:
; Ruben.van.Hooidonk@noaa.gov

; ***********************************************
; get_cortad_tsa_dhw.ncl
;
; Example NCL script for reading CoRTAD version 3
; data using NODC's Hyrax OPeNDAP server and 
; saving the results to a local netCDF file.
;
; This script was provided as-is by Ruben 
; van.Hooidonk (Ruben.van.Hooidonk@noaa.gov) on
; 15 February 2011.
;
; ***********************************************

; ***********************************************
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
;************************************************


begin
;************************************************
; specify output directory
diro = "~/cortad"

;************************************************
; loop over desired rows and columns
do row =4,5
do col_index =9,10
col = sprinti("%0.2i", col_index) 		; make sure col is of the shape 00, 01, ... 15  

;************************************************
; specify the filename and desired variables, here Latitude, Longitude, TSA_DHW and time
;
; These are the variables (and dimensions) available in the cortad v3 files:
;     AllBad[512][512];
;     ClimSST[512][512][53];
;     FilledSST[512][512][1462];
;     Harmonics[512][512][5];
;     Land[512][512];
;     Latitude[1][512];
;     Longitude[1][512];
;     MedfillSST[512][512][1462];
;     NumberGood[512][512];
;     SSTA[512][512][1462];
;     SSTA_DHW[512][512][1462];
;     SSTA_DHW_Stats[512][512][3];
;     SSTA_Frequency[512][512][1462];
;     SSTA_Frequency_Stats[512][512][3];
;     SSTA_Stats[512][512][5];
;     SST_Stats[512][512][4];
;     TSA[512][512][1462];
;     TSA_DHW[512][512][1462];
;     TSA_DHW_Stats[512][512][3];
;     TSA_Frequency[512][512][1462];
;     TSA_Frequency_Stats[512][512][3];
;     TSA_Stats[512][512][3];
;     Time[1][1462];
;     WeeklySST[512][512][1462];



filename = "http://data.nodc.noaa.gov/opendap/cortad/Version3/cortadv3_row0"+row+"_col"+col+".h5?Latitude[0:1:0][0:1:511],Longitude[0:1:0][0:1:511],TSA_DHW[0:1:511][0:1:511][0:1:1461],Time[0:1:0][0:1:1461]"


a =addfile(filename,"r") 			; read in the file name
tsa_dhw = a->TSA_DHW 				; read in actual TSA_DHW data
lat = a->Latitude				; read in Latitudes
lon = a->Longitude				; read in Longitudes
time = a->Time					; read in Time

;************************************************
;  add lat, lon, and time to the tsa_dhw variable as coordinate variables / dimensions
tsa_dhw!0="lat"
tsa_dhw!1="lon"
tsa_dhw!2="time"
tsa_dhw&lat=lat(0,:)
tsa_dhw&lon=lon(0,:)
tsa_dhw&time=time(0,:)
c = tsa_dhw({time|:},{lat|:},{lon|:})  		; just a re-ordering of the dimensions so that we can later make time an unlimited dimension

;************************************************
; create a netcdf file
ncdf = addfile(diro+"/cortad_v3_tsa_dhw_r"+row+"_c"+col+".nc","c")
	filedimdef (ncdf, "time", -1, True) 	; create unlimited time dimension
ncdf->dhw=c 					; write out the re-ordered variable

;************************************************
; clean up 
delete(a)
delete(tsa_dhw)
delete(c)
delete(ncdf)
delete(lat)
delete(lon)
delete(time)

end do  				; close loop
end do					; close loop
end