Updated NCL script used for computation testing

test/ncl_get_var.ncl

@@ -22,7 +22,7 @@
                 "geopt", "helicity", "lat", "lon", "omg", "p", "pressure", \
                 "pvo", "pw", "rh2", "rh", "slp", "ter", "td2", "td", "tc", \
                 "theta", "tk", "tv", "twb", "updraft_helicity", "ua", "va", \
-                "wa", "uvmet10", "uvmet", "z", "cfrac"/]
+                "wa", "uvmet10", "uvmet", "z", "cfrac", "height_agl" /]
     
     unique_dimname_list = NewList("fifo")
     unique_dimsize_list = NewList("fifo")
@@ -40,91 +40,11 @@
     ; variables with unique names.  This is the reason for all of the 
     ; name attribute stuff below.     
     do i = 0, ListCount(wrf_vars) - 1
-       
        print("working on: " + wrf_vars[i])
        v := wrf_user_getvar(input_file, wrf_vars[i], time)
-       
+       fout->$wrf_vars[i]$ = v
-       ;if (wrf_vars[i] .eq. "avo") then
-       ;    print(v)
-       ;end if
-       
-       ; pw is written in pure NCL and does not contain dimension names
-       ; so manually creating the dimension names here
-       if (wrf_vars[i] .eq. "pw") then
-           dim_names := (/"south_north", "west_east"/)
-           dim_sizes := dimsizes(v)
-       else
-           dim_names := getvardims(v)
-           dim_sizes := dimsizes(v)
-       end if 
-       
-       vardata@$wrf_vars[i]$ := v
-       vardimnamedata@$wrf_vars[i]$ := dim_names
-       ListAppend(vardata_list,vardata@$wrf_vars[i]$) 
-       ListAppend(full_vardimname_list, vardimnamedata@$wrf_vars[i]$)
-       ;print(vardata_list)
-       
-       dimname=True
-       dimsize=True
-       
-       ; Determine the unique dimensions names, which will be used when 
-       ; creating the output NetCDF file
-       do j=0, dimsizes(dim_sizes)-1
-           ;print(dim_names)
-           ;print(dim_names(j))
-           
-           name_id = sprintf("dimname_%i",i*j)
-           size_id = sprintf("dimsize_%i",i*j)
-           
-           dimname@$name_id$ = dim_names(j) 
-           dimsize@$size_id$ = dim_sizes(j)
-           
-           has_name = False
-           do k=0, ListCount(unique_dimname_list)-1
-               if ((/unique_dimname_list[k]/) .eq. (/dimname@$name_id$/)) then
-                   has_name = True
-               end if
-           end do
-
-           if (.not. has_name) then 
-               ;print("inserting: " + dimname@$name_id$)
-               ListAppend(unique_dimname_list, dimname@$name_id$)
-               ListAppend(unique_dimsize_list, dimsize@$size_id$)
-           end if
-           
-       end do
     end do 
     
-    setfileoption(fout,"DefineMode",True)
-    
-    ; Set global attributes
-    f_att               = True            ; assign file attributes
-    f_att@title         = "NCL generated netCDF file"
-    f_att@Conventions   = "None"
-    fileattdef(fout, f_att)      ; copy file attributes
-    
-    ; Set up the NetCDF dimensions
-    d_names = new(ListCount(unique_dimname_list), string)
-    d_sizes = new(ListCount(unique_dimname_list), integer)
-    d_unlim  = new(ListCount(unique_dimname_list), logical)
-    
-    ; Note: Need to do this copy since NCL can't coerce the list data to 
-    ; array data
-    do i=0, ListCount(unique_dimname_list) - 1
-        d_names(i) = unique_dimname_list[i]
-        d_sizes(i) = unique_dimsize_list[i]
-        d_unlim(i) = False
-    end do
-    
-    filedimdef(fout, d_names, d_sizes, d_unlim)
-    
-    ; Save the variables to the NetCDF file
-    do i=0, ListCount(vardata_list)-1       
-       d := vardata_list[i]
-       filevardef(fout, wrf_vars[i], typeof(d), full_vardimname_list[i])
-       filevarattdef(fout,wrf_vars[i], d)
-       fout->$wrf_vars[i]$ = (/d/)
-    end do
     
     ; Do the interpolation routines manually
     ; 3D vertical cross section
@@ -136,14 +56,71 @@
     
     ht_cross = wrf_user_intrp3d(z,p,"v",pivot,90.0,False)
     p_cross = wrf_user_intrp3d(p,z,"v",pivot,90.0,False)
-    
+    p_cross!0 = "Vertical_p"
-    ht_cross_dims = dimsizes(ht_cross)
+    p_cross!1 = "Horizontal_p"
-    p_cross_dims = dimsizes(p_cross)
+    
+    fout->ht_cross = ht_cross
+    fout->p_cross = p_cross
+    
+    ; For the new cross section routine
+    xopt = True
+    xopt@use_pivot = True
+    xopt@angle = 90.0
+    xopt@file_handle = input_file
+    xopt@timeidx = 0
+    xopt@linecoords = True
+    
+    ht_vertcross1 = wrf_user_vertcross(z, p, pivot, xopt)
+    
+    fout->ht_vertcross1 = ht_vertcross1
+    
+    ; For the new cross section routine
+    xopt := True
+    xopt@use_pivot = True
+    xopt@angle = 90.0
+    xopt@levels = (/1000., 850., 700., 500., 250./)
+    xopt@file_handle = input_file
+    xopt@timeidx = 0
+    xopt@linecoords = True
+    
+    ht_vertcross2 = wrf_user_vertcross(z, p, pivot, xopt)
+    ht_vertcross2!0 = "vertical2"
+    ht_vertcross2!1 = "cross_line_idx2"
+    
+    
+    fout->ht_vertcross2 = ht_vertcross2
+    
+    ; Note for the default file:
+    ; >>> np.amin(lats)
+    ; 32.77895
+    ; >>> np.amax(lats)
+    ; 40.209736
+    ; >>> np.amin(lons)
+    ; -104.71738
+    ; >>> np.amax(lons)
+    ; -95.153076
+    
+    start_end = (/ -102.0, 35.0, -98.0, 37.0 /)
+    ; For the new cross section routine
+    xopt := True
+    xopt@use_pivot = False
+    xopt@latlon = True
+    xopt@file_handle = input_file
+    xopt@timeidx = 0
+    xopt@linecoords = True
+    
+    ht_vertcross3 = wrf_user_vertcross(z, p, start_end, xopt)
+    ht_vertcross3!0 = "vertical3"
+    ht_vertcross3!1 = "cross_line_idx3"
+    
+    fout->ht_vertcross3 = ht_vertcross3
     
     ; 3D horizontal interpolation
     plev = 500.   ; 500 MB
     z_500 = wrf_user_intrp3d(z,p,"h",plev,0.,False)
-    z_500_dims = dimsizes(z_500)
+    ;z_500_dims = dimsizes(z_500)
+    
+    fout->z_500 = z_500
     
     
     ; 2D interpolation along line
@@ -152,30 +129,7 @@
     pivot = (/ dimst2(1)/2, dimst2(0)/2 /)
     
     t2_line = wrf_user_intrp2d(t2, pivot, 90.0, False)
-    t2_line_dims = dimsizes(t2_line) 
+    
-    
-    filedimdef(fout, (/"ht_cross_vert", "ht_cross_horiz", "p_cross_vert", "p_cross_horiz"/), \
-                     (/ht_cross_dims(0), ht_cross_dims(1), p_cross_dims(0), p_cross_dims(1)/), \ 
-                     (/False,False,False,False/))
-    filedimdef(fout, (/"z500_vert", "z500_horiz"/), \
-                     (/z_500_dims(0), z_500_dims(1) /), \
-                     (/False, False/))                 
-    filedimdef(fout, (/"t2_line_horiz"/), \
-                     (/t2_line_dims(0) /), \
-                     (/False/))
-    
-    filevardef(fout, "ht_cross", typeof(ht_cross), (/"ht_cross_vert", "ht_cross_horiz"/))
-    filevardef(fout, "p_cross", typeof(p_cross), (/"p_cross_vert", "p_cross_horiz"/))
-    filevardef(fout, "z_500", typeof(z_500), (/"z500_vert", "z500_horiz"/))
-    filevardef(fout, "t2_line", typeof(t2_line), (/"t2_line_horiz"/))
-    
-    filevarattdef(fout,"ht_cross", ht_cross)
-    filevarattdef(fout,"p_cross", p_cross)
-    filevarattdef(fout, "z_500", z_500)
-    filevarattdef(fout, "t2_line", t2_line)
-    fout->ht_cross = (/ht_cross/)
-    fout->p_cross = (/p_cross/)
-    fout->z_500 = (/z_500/)
     fout->t2_line = (/t2_line/)
     
     
@@ -191,32 +145,17 @@
     opts@logP        = True 
     
     fld1_intrp = wrf_user_vert_interp(input_file,fld1,vert_coord,interp_levels,opts)
+    fld1_intrp!1 = "interp_levels1"
     
-    fld1_dims = dimsizes(fld1_intrp)
+    fout->fld_tk_theta = fld1_intrp
-    
-    filedimdef(fout, (/"fld1_time", "fld1_levs", "fld1_sn", "fld1_we"/), \
-                     (/fld1_dims(0), fld1_dims(1), fld1_dims(2), fld1_dims(3)/), \ 
-                     (/False,False,False,False/))
-                     
-    filevardef(fout, "fld_tk_theta", typeof(fld1_intrp), (/"fld1_time", "fld1_levs", "fld1_sn", "fld1_we"/))
-    filevarattdef(fout,"fld_tk_theta", fld1_intrp)
-    fout->fld_tk_theta = (/fld1_intrp/)
     
     ; interp t to theta-e
     fld2 = fld1
     vert_coord := "theta-e"
     fld2_intrp = wrf_user_vert_interp(input_file,fld2,vert_coord,interp_levels,opts)
+    fld2_intrp!1 = "interp_levels2"
     
-    
+    fout->fld_tk_theta_e = fld2_intrp
-    fld2_dims = dimsizes(fld2_intrp)
-    
-    filedimdef(fout, (/"fld2_time", "fld2_levs", "fld2_sn", "fld2_we"/), \
-                     (/fld2_dims(0), fld2_dims(1), fld2_dims(2), fld2_dims(3)/), \ 
-                     (/False,False,False,False/))
-                     
-    filevardef(fout, "fld_tk_theta_e", typeof(fld2_intrp), (/"fld2_time", "fld2_levs", "fld2_sn", "fld2_we"/))
-    filevarattdef(fout,"fld_tk_theta_e", fld2_intrp)
-    fout->fld_tk_theta_e = (/fld2_intrp/)
     
     
     ; interp t to pressure
@@ -224,34 +163,19 @@
     vert_coord := "pressure"
     interp_levels    := (/850,500/)
     fld3_intrp = wrf_user_vert_interp(input_file,fld3,vert_coord,interp_levels,opts)
+    fld3_intrp!1 = "interp_levels3"
     
+    fout->fld_tk_pres = fld3_intrp
     
-    fld3_dims = dimsizes(fld3_intrp)
-    
-    filedimdef(fout, (/"fld3_time", "fld3_levs", "fld3_sn", "fld3_we"/), \
-                     (/fld3_dims(0), fld3_dims(1), fld3_dims(2), fld3_dims(3)/), \ 
-                     (/False,False,False,False/))
-                     
-    filevardef(fout, "fld_tk_pres", typeof(fld3_intrp), (/"fld3_time", "fld3_levs", "fld3_sn", "fld3_we"/))
-    filevarattdef(fout,"fld_tk_pres", fld3_intrp)
-    fout->fld_tk_pres = (/fld3_intrp/)
     
     ; interp t to ght_msl
     fld4 = fld1
     vert_coord := "ght_msl"
     interp_levels    := (/1,2/)
     fld4_intrp = wrf_user_vert_interp(input_file,fld4,vert_coord,interp_levels,opts)
+    fld4_intrp!1 = "interp_levels4"
     
-    
+    fout->fld_tk_ght_msl = fld4_intrp
-    fld4_dims = dimsizes(fld4_intrp)
-    
-    filedimdef(fout, (/"fld4_time", "fld4_levs", "fld4_sn", "fld4_we"/), \
-                     (/fld4_dims(0), fld4_dims(1), fld4_dims(2), fld4_dims(3)/), \ 
-                     (/False,False,False,False/))
-                     
-    filevardef(fout, "fld_tk_ght_msl", typeof(fld4_intrp), (/"fld4_time", "fld4_levs", "fld4_sn", "fld4_we"/))
-    filevarattdef(fout,"fld_tk_ght_msl", fld4_intrp)
-    fout->fld_tk_ght_msl = (/fld4_intrp/)
     
     
     ; interp t to ght_agl
@@ -259,17 +183,9 @@
     vert_coord := "ght_agl"
     interp_levels    := (/1,2/)
     fld5_intrp = wrf_user_vert_interp(input_file,fld1,vert_coord,interp_levels,opts)
+    fld5_intrp!1 = "interp_levels5"
     
-    
+    fout->fld_tk_ght_agl = fld5_intrp
-    fld5_dims = dimsizes(fld5_intrp)
-    
-    filedimdef(fout, (/"fld5_time", "fld5_levs", "fld5_sn", "fld5_we"/), \
-                     (/fld5_dims(0), fld5_dims(1), fld5_dims(2), fld5_dims(3)/), \ 
-                     (/False,False,False,False/))
-                     
-    filevardef(fout, "fld_tk_ght_agl", typeof(fld5_intrp), (/"fld5_time", "fld5_levs", "fld5_sn", "fld5_we"/))
-    filevarattdef(fout,"fld_tk_ght_agl", fld5_intrp)
-    fout->fld_tk_ght_agl = (/fld5_intrp/)
     
     ; interp ht to pres
     fld6 = wrf_user_getvar(input_file, "height", -1)
@@ -277,17 +193,9 @@
     opts@field_type  = "ght"
     interp_levels    := (/500,50/)
     fld6_intrp = wrf_user_vert_interp(input_file,fld6,vert_coord,interp_levels,opts)
+    fld6_intrp!1 = "interp_levels6"
     
-    
+    fout->fld_ht_pres = fld6_intrp
-    fld6_dims = dimsizes(fld6_intrp)
-    
-    filedimdef(fout, (/"fld6_time", "fld6_levs", "fld6_sn", "fld6_we"/), \
-                     (/fld6_dims(0), fld6_dims(1), fld6_dims(2), fld6_dims(3)/), \ 
-                     (/False,False,False,False/))
-                     
-    filevardef(fout, "fld_ht_pres", typeof(fld6_intrp), (/"fld6_time", "fld6_levs", "fld6_sn", "fld6_we"/))
-    filevarattdef(fout,"fld_ht_pres", fld6_intrp)
-    fout->fld_ht_pres = (/fld6_intrp/)
     
     
     ; interp pres to theta
@@ -296,35 +204,21 @@
     opts@field_type  = "pressure"
     interp_levels    := (/200,300,500,1000/)
     fld7_intrp = wrf_user_vert_interp(input_file,fld7,vert_coord,interp_levels,opts)
+    fld7_intrp!1 = "interp_levels7"
     
-    
-    fld7_dims = dimsizes(fld7_intrp)
-    
-    filedimdef(fout, (/"fld7_time", "fld7_levs", "fld7_sn", "fld7_we"/), \
-                     (/fld7_dims(0), fld7_dims(1), fld7_dims(2), fld7_dims(3)/), \ 
-                     (/False,False,False,False/))
-                     
-    filevardef(fout, "fld_pres_theta", typeof(fld7_intrp), (/"fld7_time", "fld7_levs", "fld7_sn", "fld7_we"/))
-    filevarattdef(fout,"fld_pres_theta", fld7_intrp)
     fout->fld_pres_theta = (/fld7_intrp/)
     
+    
     ; interp theta-e to pressure
     fld8 = wrf_user_getvar(input_file, "eth", -1)
     vert_coord := "pressure"
     opts@field_type  = "T"
     interp_levels    := (/850,500,5/)
     fld8_intrp = wrf_user_vert_interp(input_file,fld8,vert_coord,interp_levels,opts)
+    fld8_intrp!1 = "interp_levels8"
     
+    fout->fld_thetae_pres = fld8_intrp
     
-    fld8_dims = dimsizes(fld8_intrp)
-    
-    filedimdef(fout, (/"fld8_time", "fld8_levs", "fld8_sn", "fld8_we"/), \
-                     (/fld8_dims(0), fld8_dims(1), fld8_dims(2), fld8_dims(3)/), \ 
-                     (/False,False,False,False/))
-                     
-    filevardef(fout, "fld_thetae_pres", typeof(fld8_intrp), (/"fld8_time", "fld8_levs", "fld8_sn", "fld8_we"/))
-    filevarattdef(fout,"fld_thetae_pres", fld8_intrp)
-    fout->fld_thetae_pres = (/fld8_intrp/)
     
     ; lat/lon to x/y and x/y to lat/lon routines
     lats = (/-55, -60, -65 /)
@@ -334,24 +228,8 @@
     ij = wrf_user_ll_to_ij(input_file, lons, lats, True) 
     ll = wrf_user_ij_to_ll(input_file, i_s, j_s, True)
     
-    ij_dims = dimsizes(ij)
+    fout->ij = ij
-    ll_dims = dimsizes(ll)
+    fout->ll = ll
-    
-    filedimdef(fout, (/"i_j", "ij_idx"/), \
-                     (/ij_dims(0), ij_dims(1)/), \ 
-                     (/False,False/))
-                     
-    filedimdef(fout, (/"lat_lon", "ll_idx"/), \
-                     (/ll_dims(0), ll_dims(1)/), \ 
-                     (/False,False/))
-                     
-    filevardef(fout, "ij", typeof(ij), (/"i_j", "ij_idx"/))
-    filevardef(fout, "ll", typeof(ll), (/"lat_lon", "ll_idx"/))
-    filevarattdef(fout,"ij", ij)
-    filevarattdef(fout,"ll", ll)
-    
-    fout->ij = (/ij/)
-    fout->ll = (/ll/)
     
     delete(fout)