load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"

start_day = 27
end_day   = 28

do dy = start_day,end_day
  if (dy .lt. 10) then
    dy2 = "0"+dy  
  else
    dy2 = dy
  end if

  if (dy .eq. 27) then
    start_hr = 0
    end_hr   = 23
  end if

  if (dy .eq. 28) then
    start_hr = 0
    end_hr   = 6
  end if

  do hr = start_hr,end_hr
    if (hr .lt. 10) then 
      hr2 = "0"+hr
    else
      hr2 = hr
    end if
    
    ;Change year/month here:
    print("Processing 04/"+dy2+" "+hr2+":00")

    ;Rename the below file directory and name to what your file directory and names are.
    ;You can run this on both d01 and d02
    ;Note that "w" is used at end of addfile call so that we can write new dimensions and variables to the wrfout NetCDF
    a = addfile("/Volumes/ready-01-2015/Home/Rutgers/Project/Sea_Breeze/ruwrf/Case03/wrf/wrfout_d01_2013-04-"+dy2+"_"+hr2+":00:00_r2.nc","w")

    ur = a->Ur(0,:,:,:) ;Ur(time,z,y,x)
    vr = a->Vr(0,:,:,:) ;Vr(time,z,y,x)
    wr = a->Wr(0,:,:,:) ;Wr(time,z,y,x)
    terrainhgt = a->HGT(0,:,:) ;pull terrain height from wrfout NetCDF

    time = -1

    z = wrf_user_getvar(a,"z",time) ;z on mass points
    z2 = z(0,:,:,:) ;z(time,z,y,x)
    nheight = conform(z2,terrainhgt,(/1,2/)) ; z2 is a 3d array and terrainhgt2 is a 2d array
    zagl = z2 - nheight  ;subtract terrain height from model z values to get above ground level heights
     
    u3d = new((/20,299,299/),float) ;initialize 3d u, v, and w matrices
    v3d = new((/20,299,299/),float) ;v3d(z,y,x)
    w3d = new((/20,299,299/),float)
    
    u3d!0 = "Vertical3d" ;name first dimension of u3d, v3d, and w3d as "Vertical3d"
    v3d!0 = "Vertical3d"
    w3d!0 = "Vertical3d"

    i = 0
    do k = 100,2000,100     ;loop through all heights, from 100m to 2000m at 100m intervals
       u3d(i,:,:) = wrf_user_intrp3d(ur,zagl,"h",k,0.,False) ;vertically interpolate u, v, and w winds at each height
       v3d(i,:,:) = wrf_user_intrp3d(vr,zagl,"h",k,0.,False)
       w3d(i,:,:) = wrf_user_intrp3d(wr,zagl,"h",k,0.,False)
       i = i + 1
    end do
   
    dzdx = dimsizes(u3d) ;pull dimsizes of u3d so that we can define a new dimension called "Vertical3d"
    nvert = dzdx(0) ;first dimsize of u3d is number of vertical layers

    dim_names = (/ "Vertical3d" /) ;write new dimension called "Vertical3d" to wrfout NetCDF
    dim_sizes = (/ nvert  /)
    dimUnlim = (/ False  /)
    filedimdef(a, dim_names, dim_sizes, dimUnlim )
 
    filevardef(a, "U3d", typeof(u3d), getvardims(u3d)) ;write new variable called U3d to wrfout NetCDF
    a->U3d = (/u3d/)
    U3d = a->U3d
    U3d@description = "u wind rotated to earth coord. (m/s) at 100,200,300...2000m ABOVE GROUND LEVEL (AGL)"
    U3d@units = "m/s"
    a->U3d = U3d
    
    filevardef(a, "V3d", typeof(v3d), getvardims(v3d)) ;write new variable called V3d to wrfout NetCDF
    a->V3d = (/v3d/)
    V3d = a->V3d
    V3d@description = "v wind rotated to earth coord. (m/s) at 100,200,300...2000m ABOVE GROUND LEVEL (AGL)"
    V3d@units = "m/s"
    a->V3d = V3d
    
    filevardef(a, "W3d", typeof(w3d), getvardims(w3d)) ;write new variable called W3d to wrfout NetCDF
    a->W3d = (/w3d/)
    W3d = a->W3d
    W3d@description = "w wind rotated to earth coord. (m/s) at 100,200,300...2000m ABOVE GROUND LEVEL (AGL)"
    W3d@units = "m/s"
    a->W3d = W3d
 
  end do
end do