Updated tests for new horizontal interpolation

test/ncl_get_var.ncl

@@ -6,6 +6,7 @@
     
     if (.not. isvar("in_file")) then
         in_file = "/Users/ladwig/Documents/wrf_files/wrfout_d01_2010-06-13_21:00:00.nc"
+        ;in_file = "/Users/ladwig/Documents/wrf_files/wrf_vortex_multi/wrfout_d02_2005-08-28_00:00:00.nc"
     end if 
     
     if (.not. isvar("out_file")) then
@@ -118,10 +119,51 @@
     fout->ht_vertcross3 = ht_vertcross3
     
     ; 3D horizontal interpolation
+    
+    ; First, check backwards compat
     plev = 500.   ; 500 MB
+    hlev = 5000; ; 5000 m
+    
     z_500 = wrf_user_intrp3d(z,p,"h",plev,0.,False)
+    p_5000 = wrf_user_intrp3d(p,z,"h",hlev,0.,False)
     
     fout->z_500 = z_500
+    fout->p_5000 = p_5000
+    
+    plev := (/1000., 850., 500., 250./) 
+    hlev := (/500., 2500., 5000., 10000. /)
+    z_multi = wrf_user_intrp3d(z,p,"h",plev,0.,False)
+    p_multi = wrf_user_intrp3d(p,z,"h",hlev,0.,False)
+    
+    fout->z_multi = z_multi
+    fout->p_multi = p_multi
+    
+    ; Now check the new routine
+    
+    plev := 500.   ; 500 MB
+    hlev := 5000   ; 5000 m
+    
+    z2_500 = wrf_user_interplevel(z,p,plev,False)
+    p2_5000 = wrf_user_interplevel(p,z,hlev,False)
+    
+    fout->z2_500 = z2_500
+    fout->p2_5000 = p2_5000
+    
+    
+    plev := (/1000., 850., 500., 250./) 
+    hlev := (/500., 2500., 5000., 10000. /)
+    z2_multi = wrf_user_interplevel(z,p,plev,False)
+    p2_multi = wrf_user_interplevel(p,z,hlev,False)
+    
+    fout->z2_multi = z2_multi
+    fout->p2_multi = p2_multi
+    
+    pblh = wrf_user_getvar(input_file, "PBLH", 0)
+    opts := False
+    opts@inc2dlevs = True
+    p_lev2d = wrf_user_interplevel(p, z, pblh, opts)
+    
+    fout->p_lev2d = p_lev2d
     
     
     ; 2D interpolation along line

test/utests.py

@@ -300,6 +300,17 @@ def make_interp_test(varname, wrf_in, referent, multi=False,
         
         if (varname == "interplevel"):
             ref_ht_500 = _get_refvals(referent, "z_500", repeat, multi)
+            ref_p_5000 = _get_refvals(referent, "p_5000", repeat, multi)
+            ref_ht_multi = _get_refvals(referent, "z_multi", repeat, multi)
+            ref_p_multi = _get_refvals(referent, "p_multi", repeat, multi)
+            
+            ref_ht2_500 = _get_refvals(referent, "z2_500", repeat, multi)
+            ref_p2_5000 = _get_refvals(referent, "p2_5000", repeat, multi)
+            ref_ht2_multi = _get_refvals(referent, "z2_multi", repeat, multi)
+            ref_p2_multi = _get_refvals(referent, "p2_multi", repeat, multi)
+            
+            ref_p_lev2d = _get_refvals(referent, "p_lev2d", repeat, multi)
+            
             hts = getvar(in_wrfnc, "z", timeidx=timeidx)
             p = getvar(in_wrfnc, "pressure", timeidx=timeidx)
             wspd_wdir = getvar(in_wrfnc, "wspd_wdir", timeidx=timeidx)
@@ -308,40 +319,67 @@ def make_interp_test(varname, wrf_in, referent, multi=False,
             hts_500 = interplevel(to_np(hts), to_np(p), 500)
             hts_500 = interplevel(hts, p, 500)
             
-            print(hts_500)
+            # Note: the '*2*' versions in the reference are testing 
+            # against the new version of interplevel in NCL
+            nt.assert_allclose(to_np(hts_500), ref_ht_500)
+            nt.assert_allclose(to_np(hts_500), ref_ht2_500)
             
-            #nt.assert_allclose(to_np(hts_500), ref_ht_500)
+            # Make sure the numpy versions work first
+            p_5000 = interplevel(to_np(p), to_np(hts), 5000)
+            p_5000 = interplevel(p, hts, 5000)
             
-            hts_multi= interplevel(to_np(hts), to_np(p), [1000,500,250])
-            hts_multi = interplevel(hts, p, [1000,500,250])
             
-            print(hts_multi)
+            nt.assert_allclose(to_np(p_5000), ref_p_5000)
+            nt.assert_allclose(to_np(p_5000), ref_p2_5000)
             
-            pblh = getvar(in_wrfnc, "PBLH", timeidx=timeidx)
+            hts_multi= interplevel(to_np(hts), to_np(p), 
-            hts_pblh = interplevel(p, hts, pblh)
+                                   [1000., 850., 500., 250.])
+            hts_multi = interplevel(hts, p, [1000., 850., 500., 250.])
+            
+            nt.assert_allclose(to_np(hts_multi), ref_ht_multi)
+            nt.assert_allclose(to_np(hts_multi), ref_ht2_multi)
             
-            print(hts_pblh)
+            p_multi= interplevel(to_np(p), to_np(hts), 
+                                   [500., 2500., 5000., 10000. ])
+            p_multi = interplevel(p, hts, [500., 2500., 5000., 10000. ])
             
-            #nt.assert_allclose(to_np(hts_500), ref_ht_500)
+            nt.assert_allclose(to_np(p_multi), ref_p_multi)
+            nt.assert_allclose(to_np(p_multi), ref_p2_multi)
             
+            pblh = getvar(in_wrfnc, "PBLH", timeidx=timeidx)
+            p_lev2d = interplevel(to_np(p), to_np(hts), to_np(pblh))
+            p_lev2d = interplevel(p, hts, pblh)
+            nt.assert_allclose(to_np(p_lev2d), ref_p_lev2d)
+            
+            # Just make sure these run below
             wspd_wdir_500 = interplevel(to_np(wspd_wdir), to_np(p), 500)
             wspd_wdir_500 = interplevel(wspd_wdir, p, 500)
-            print(wspd_wdir_500)
+            #print(wspd_wdir_500)
             
             wspd_wdir_multi= interplevel(to_np(wspd_wdir), 
                                          to_np(p), [1000,500,250])
             wdpd_wdir_multi = interplevel(wspd_wdir, p, [1000,500,250])
             
-            print(wdpd_wdir_multi)
             
             wspd_wdir_pblh = interplevel(to_np(wspd_wdir), to_np(hts), pblh)
             wspd_wdir_pblh = interplevel(wspd_wdir, hts, pblh)
-            print(wspd_wdir_pblh)
             
-            wspd_wdir_pblh = interplevel(to_np(wspd_wdir), 
+            if multi:
+                wspd_wdir_pblh_2 = interplevel(to_np(wspd_wdir), 
                                          to_np(hts), pblh[0,:])
-            wspd_wdir_pblh = interplevel(wspd_wdir, hts, pblh[0,:])
+                wspd_wdir_pblh_2 = interplevel(wspd_wdir, hts, pblh[0,:])
-            print(wspd_wdir_pblh)
+                
+                # Since PBLH doesn't change in this case, it should match
+                # the 0 time from previous computation. Note that this 
+                # only works when the data has 1 time step that is repeated.
+                # If you use a different case with multiple times, 
+                # this will probably fail.
+                nt.assert_allclose(to_np(wspd_wdir_pblh_2[:,0,:]), 
+                                   to_np(wspd_wdir_pblh[:,0,:]))
+                
+                nt.assert_allclose(to_np(wspd_wdir_pblh_2[:,-1,:]), 
+                                   to_np(wspd_wdir_pblh[:,0,:]))
+                
             
         elif (varname == "vertcross"):
             ref_ht_cross = _get_refvals(referent, "ht_cross", repeat, multi)