load "/opt/local/ncl-6.3.0/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "/opt/local/ncl-6.3.0/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "/opt/local/ncl-6.3.0/lib/ncarg/nclscripts/csm/contributed.ncl"
load "/opt/local/ncl-6.3.0/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"

; --------------  BEGINING OF NCL SCRIPT ----------------


begin
;********************************************************
; read in netCDF file and make a loop for all time steps
;********************************************************
  in     = addfile("wrfout_d01_0503-01-01_01:00:00.nc","w")


;********************************************************
; Process all the time steps
;********************************************************

times = wrf_user_list_times(in)     ; get times in the file 
ntimes = dimsizes(times)            ; number of times in the file 

hfx100 = new(ntimes,float)          ; creation of a hfx vector at each time step
;print(ntimes)

do it = 0,10 ;ntimes-1                  ;Loop for the time: it= starting time

time = it

;*************************************************************************************
;  - extract HFX and height coordinates-  *
;*************************************************************************************

; sensible heat flux and Height
    hfx  = wrf_user_getvar(in,"HFX",time)        ; u averaged to mass points
    height  = wrf_user_getvar(in, "z",time) ; height is our vertical coordinate
    ter = wrf_user_getvar(in, "ter",time) ; model terrain height (HGT_M, HGT_U, HGT_V)
 
; Conform data to Terrain Height   
   nheight = conform(height,ter,(/1,2/)) ; assuming height is a 3d array and ter is a 2d array
   height = height - nheight



;*******************************************************************************
;     - Interpolate sensible heat flux at 100m height -          
;*******************************************************************************

	; Interpolate HFX to 100 Meters
    hfx_plane  = wrf_interp_1d(hfx, height, 100)
    hfx100(it) = hfx_plane(x,y)
     
  end do

;************************************************************
;  - Print the variables at the screen -                    * 
;************************************************************

npts=ntimes
fName = "HFX100.txt"
data  = new( npts, "string")
       
     
    print("  Time      HFX_100m		")
     
    do it = 0,10 ;ntimes-1

        print (sprintf("%5.0f",it)    +" " \
          +sprintf("%21.2f", hfx100(it)) +"  " )
                
 end do                     ; end of time loop

	asciiwrite (fName , data)

end