Cleaned up cape fortran routine. Fixed issue with generators not being copied correctly due to changes in python (the python works as expected now). Fixed a uniqueness problem with coordinate caching which was causing problems in jupyter notebook when files were changed. Fixed an issue with the cache test script failing due to unitialized thread local data in child threads. Fixes #34. Fixes #14.

fortran/eqthecalc.f90

@@ -32,7 +32,6 @@ SUBROUTINE DEQTHECALC(qvp, tmk, prs, eth, miy, mjx, mkzh)
     REAL(KIND=8) :: tlcl
     INTEGER :: i, j, k
 
-    ! Note:  removing temporaries did not improve performance for this algorithm.
     !$OMP PARALLEL DO COLLAPSE(3) PRIVATE(i, j, k, q, t, p, e, tlcl)
     DO k = 1,mkzh
         DO j = 1,mjx

fortran/rip_cape.f90

@@ -227,6 +227,7 @@ SUBROUTINE DPFCALC(prs, sfp, pf, mix, mjy, mkzh, ter_follow)
 
     INTEGER :: i,j,k
     
+    !$OMP PARALLEL DO COLLAPSE(3)
     DO j = 1,mjy
         DO i = 1,mix
             DO k = 1,mkzh
@@ -245,6 +246,8 @@ SUBROUTINE DPFCALC(prs, sfp, pf, mix, mjy, mkzh, ter_follow)
         END DO
     END DO
 
+    !$OMP END PARALLEL DO
+
     RETURN
 
 END SUBROUTINE DPFCALC
@@ -351,10 +354,7 @@ SUBROUTINE DCAPECALC3D(prs,tmk,qvp,ght,ter,sfp,cape,cin,&
     !  calculated the pressure at full sigma levels (a set of pressure
     !  levels that bound the layers represented by the vertical grid points)
 
-    !CALL cpu_time(t1)
-    !CALL OMP_SET_NUM_THREADS(16)
-
-    !$OMP PARALLEL DO
+    !$OMP PARALLEL DO COLLAPSE(3)
     DO j = 1,mjy
         DO i = 1,mix
             DO k = 1,mkzh
@@ -403,10 +403,6 @@ SUBROUTINE DCAPECALC3D(prs,tmk,qvp,ght,ter,sfp,cape,cin,&
                 
                 zlcl = ght_new(kpar,i,j) + (tmk_new(kpar,i,j) - tlcl)/(G/CP * (1.D0 + CPMD*qvp_new(kpar,i,j)))
 
-           !   DO k = kpar,1,-1
-           !       tmklift_new(k) = TONPSADIABAT(ethpari, prs_new(k,i,j), psadithte, psadiprs,&
-           !                                  psaditmk, GAMMA, errstat, errmsg)
-           !   END DO
                 ! Calculate buoyancy and relative height of lifted parcel at
                 ! all levels, and store in bottom up arrays.  add a level at the lcl,
                 ! and at all points where buoyancy is zero.
@@ -428,7 +424,8 @@ SUBROUTINE DCAPECALC3D(prs,tmk,qvp,ght,ter,sfp,cape,cin,&
 
                     ! Model level is below lcl
                     IF (ght_new(k,i,j) .LT. zlcl) THEN
-                      tmklift = tmk_new(kpar,i,j) - G/(CP * (1.D0 + CPMD*qvp_new(kpar,i,j))) * (ght_new(k,i,j) - ght_new(kpar,i,j))
+                        tmklift = tmk_new(kpar,i,j) - G/(CP * (1.D0 + CPMD*qvp_new(kpar,i,j)))*&
+                                  (ght_new(k,i,j) - ght_new(kpar,i,j))
                         tvenv = tmk_new(k,i,j)*(EPS + qvp_new(k,i,j))/(EPS*(1.D0 + qvp_new(k,i,j)))
                         tvlift = tmklift*(EPS + qvp_new(kpar,i,j))/(EPS*(1.D0 + qvp_new(kpar,i,j)))
                         ghtlift = ght_new(k,i,j)
@@ -436,8 +433,10 @@ SUBROUTINE DCAPECALC3D(prs,tmk,qvp,ght,ter,sfp,cape,cin,&
                         ! This model level and previous model level straddle the lcl,
                         ! so first create a new level in the bottom-up array, at the lcl.
                         facden = 1.0/(ght_new(k,i,j) - ght_new(k+1,i,j))
-                      tmkenv = tmk_new(k+1,i,j)*((ght_new(k,i,j)-zlcl)*facden) + tmk_new(k,i,j)*((zlcl-ght_new(k+1,i,j))*facden)
-                      qvpenv = qvp_new(k+1,i,j)*((ght_new(k,i,j)-zlcl)*facden) + qvp_new(k,i,j)*((zlcl-ght_new(k+1,i,j))*facden)
+                        tmkenv = tmk_new(k+1,i,j)*((ght_new(k,i,j)-zlcl)*facden) + tmk_new(k,i,j)*&
+                                 ((zlcl-ght_new(k+1,i,j))*facden)
+                        qvpenv = qvp_new(k+1,i,j)*((ght_new(k,i,j)-zlcl)*facden) + qvp_new(k,i,j)*&
+                                 ((zlcl-ght_new(k+1,i,j))*facden)
                         tvenv = tmkenv* (EPS + qvpenv) / (EPS * (1.D0 + qvpenv))
                         tvlift = tlcl* (EPS + qvp_new(kpar,i,j)) / (EPS *(1.D0 + qvp_new(kpar,i,j)))
                         ghtlift = zlcl
@@ -570,8 +569,7 @@ SUBROUTINE DCAPECALC3D(prs,tmk,qvp,ght,ter,sfp,cape,cin,&
         END DO
     END DO
     !$OMP END PARALLEL DO
-    !CALL cpu_time(t2)
-    !print *,'Time taken in seconds ',(t2-t1)
+
     RETURN
 END SUBROUTINE DCAPECALC3D
 
@@ -590,8 +588,6 @@ END SUBROUTINE DCAPECALC3D
 !   the cape and cin arrays.  Also, LCL and LFC heights
 !   are put in the k=mkzh-1 and k=mkzh-2 slabs of the cin array.
 !
-
-
 ! Important!  The z-indexes must be arranged so that mkzh (max z-index) is the
 ! surface pressure.  So, pressure must be ordered in ascending order before
 ! calling this routine.  Other variables must be ordered the same (p,tk,q,z).
@@ -711,10 +707,6 @@ SUBROUTINE DCAPECALC2D(prs,tmk,qvp,ght,ter,sfp,cape,cin,&
         RETURN
     END IF
 
-    !CALL OMP_SET_NUM_THREADS(16)
-    !nthreads = omp_get_num_threads()
-
-
     !$OMP PARALLEL DO COLLAPSE(2) PRIVATE(tlcl, ethpari, &
     !$OMP zlcl, kk, ilcl, klcl, tmklift, tvenv, tvlift, ghtlift, &
     !$OMP facden, tmkenv, qvpenv, eslift, qvplift, buoy, benamin, &
@@ -784,7 +776,6 @@ SUBROUTINE DCAPECALC2D(prs,tmk,qvp,ght,ter,sfp,cape,cin,&
             tmkpari = (totthe/totprs)*&
                       (prs_new(kpar1,i,j)/1000.D0)**(GAMMA*(1.D0+GAMMAMD*qvp_new(kpar1,i,j)))
 
-!CALL CPU_TIME(t3)
             DO kpar = kpar1, kpar2
 
                 ! Calculate temperature and moisture properties of parcel
@@ -825,7 +816,8 @@ SUBROUTINE DCAPECALC2D(prs,tmk,qvp,ght,ter,sfp,cape,cin,&
 
                     ! Model level is below lcl
                     IF (ght_new(k,i,j) .LT. zlcl) THEN
-                      tmklift = tmk_new(kpar,i,j) - G/(CP * (1.D0 + CPMD*qvp_new(kpar,i,j))) * (ght_new(k,i,j) - ght_new(kpar,i,j))
+                        tmklift = tmk_new(kpar,i,j) - G/(CP * (1.D0 + CPMD*qvp_new(kpar,i,j)))*&
+                                  (ght_new(k,i,j) - ght_new(kpar,i,j))
                         tvenv = tmk_new(k,i,j)*(EPS + qvp_new(k,i,j))/(EPS*(1.D0 + qvp_new(k,i,j)))
                         tvlift = tmklift*(EPS + qvp_new(kpar,i,j))/(EPS*(1.D0 + qvp_new(kpar,i,j)))
                         ghtlift = ght_new(k,i,j)
@@ -833,8 +825,10 @@ SUBROUTINE DCAPECALC2D(prs,tmk,qvp,ght,ter,sfp,cape,cin,&
                         ! This model level and previous model level straddle the lcl,
                         ! so first create a new level in the bottom-up array, at the lcl.
                         facden = 1/(ght_new(k,i,j) - ght_new(k+1,i,j))
-                      tmkenv = tmk_new(k+1,i,j)*((ght_new(k,i,j)-zlcl)*facden) + tmk_new(k,i,j)*((zlcl-ght_new(k+1,i,j))*facden)
-                      qvpenv = qvp_new(k+1,i,j)*((ght_new(k,i,j)-zlcl)*facden) + qvp_new(k,i,j)*((zlcl-ght_new(k+1,i,j))*facden)
+                        tmkenv = tmk_new(k+1,i,j)*((ght_new(k,i,j)-zlcl)*facden) + tmk_new(k,i,j)*&
+                                 ((zlcl-ght_new(k+1,i,j))*facden)
+                        qvpenv = qvp_new(k+1,i,j)*((ght_new(k,i,j)-zlcl)*facden) + qvp_new(k,i,j)*&
+                                 ((zlcl-ght_new(k+1,i,j))*facden)
                         tvenv = tmkenv* (EPS + qvpenv) / (EPS * (1.D0 + qvpenv))
                         tvlift = tlcl* (EPS + qvp_new(kpar,i,j)) / (EPS *(1.D0 + qvp_new(kpar,i,j)))
                         ghtlift = zlcl
@@ -866,7 +860,6 @@ SUBROUTINE DCAPECALC2D(prs,tmk,qvp,ght,ter,sfp,cape,cin,&
                         ilcl = 2
                         CYCLE
                     END IF
-
                 END DO
 
                 kmax = kk
@@ -980,5 +973,6 @@ SUBROUTINE DCAPECALC2D(prs,tmk,qvp,ght,ter,sfp,cape,cin,&
         END DO
     END DO
     !$OMP END PARALLEL DO
+
     RETURN
 END SUBROUTINE DCAPECALC2D

src/wrf/config.py

@@ -2,13 +2,19 @@ from __future__ import (absolute_import, division, print_function,
                         unicode_literals)
 
 from threading import local
+import wrapt
 
 _local_config = local()
+
+def _init_local():
+    global _local_config
+    
     _local_config.xarray_enabled = True
     _local_config.cartopy_enabled = True
     _local_config.basemap_enabled = True
     _local_config.pyngl_enabled = True
     _local_config.cache_size = 20
+    _local_config.initialized = True
     
     try:
         from xarray import DataArray
@@ -31,6 +37,29 @@ except ImportError:
         _local_config.pyngl_enabled = False
         
 
+# Initialize the main thread's configuration
+_init_local()
+
+
+def init_local():
+    """A decorator that initializes thread local data if necessary."""
+    @wrapt.decorator
+    def func_wrapper(wrapped, instance, args, kwargs):
+        global _local_config
+        try:
+            init = _local_config.init
+        except AttributeError:
+            _init_local()
+        else:
+            if not init:
+                _init_local() 
+                
+        return wrapped(*args, **kwargs)
+        
+    return func_wrapper
+
+
+@init_local()
 def xarray_enabled():
     """Return True if xarray is installed and enabled.
     
@@ -43,18 +72,21 @@ def xarray_enabled():
     return _local_config.xarray_enabled
 
 
+@init_local()
 def disable_xarray():
     """Disable xarray."""
     global _local_config
     _local_config.xarray_enabled = False
     
 
+@init_local()
 def enable_xarray():
     """Enable xarray."""
     global _local_config
     _local_config.xarray_enabled = True
     
 
+@init_local()
 def cartopy_enabled():
     """Return True if cartopy is installed and enabled.
     
@@ -67,18 +99,21 @@ def cartopy_enabled():
     return _local_config.cartopy_enabled
 
 
+@init_local()
 def enable_cartopy():
     """Enable cartopy."""
     global _local_config
     _local_config.cartopy_enabled = True
     
 
+@init_local()
 def disable_cartopy():
     """Disable cartopy."""
     global _local_config
     _local_config.cartopy_enabled = True
     
 
+@init_local()
 def basemap_enabled():
     """Return True if basemap is installed and enabled.
     
@@ -91,17 +126,21 @@ def basemap_enabled():
     return _local_config.basemap_enabled
 
 
+@init_local()
 def enable_basemap():
     """Enable basemap."""
     global _local_config
     _local_config.basemap_enabled = True
     
 
+@init_local()
 def disable_basemap():
     """Disable basemap."""
     global _local_config
     _local_config.basemap_enabled = True
 
+
+@init_local()
 def pyngl_enabled():
     """Return True if pyngl is installed and enabled.
     
@@ -114,18 +153,21 @@ def pyngl_enabled():
     return _local_config.pyngl_enabled
 
 
+@init_local()
 def enable_pyngl():
     """Enable pyngl."""
     global _local_config
     _local_config.pyngl_enabled = True
     
 
+@init_local()  
 def disable_pyngl():
     """Disable pyngl."""
     global _local_config
     _local_config.pyngl_enabled = True
     
 
+@init_local() 
 def set_cache_size(size):
     """Set the maximum number of items that the threadlocal cache can retain.
     
@@ -144,6 +186,7 @@ def set_cache_size(size):
     _local_config.cache_size = size
     
 
+@init_local()  
 def get_cache_size():
     """Return the maximum number of items that the threadlocal cache can retain.
     

src/wrf/util.py

@@ -209,10 +209,29 @@ def _generator_copy(gen):
     module = getmodule(gen.gi_frame)
     
     if module is not None:
+        try:
+            try:
+                argd = {key:argvals.locals[key] for key in argvals.args}
+                res = module.get(funcname)(**argd)
+            except AttributeError:
+                res = getattr(module, funcname)(**argd)
+        except:
+            # This is the old way it used to work, but it looks like this was
+            # fixed by Python.  
+            try:
                 res = module.get(funcname)(**argvals.locals)
+            except AttributeError:
+                res = getattr(module, funcname)(**argvals.locals)    
     else:
         # Created in jupyter or the python interpreter
         import __main__
+        
+        try:
+            argd = {key:argvals.locals[key] for key in argvals.args}
+            res = getattr(__main__, funcname)(**argd)
+        except:
+            # This was the old way it used to work, but appears to have 
+            # been fixed by Python.
             res = getattr(__main__, funcname)(**argvals.locals)
         
     return res
@@ -2583,26 +2602,6 @@ def get_proj_params(wrfin):#, timeidx=0, varname=None):
                                                 "DX", "DY"))
     
     return proj_params
-#     multitime = is_multi_time_req(timeidx)
-#     if not multitime:
-#         time_idx_or_slice = timeidx
-#     else:
-#         time_idx_or_slice = slice(None)
-#     
-#     if varname is not None:
-#         if not is_coordvar(varname):
-#             coord_names = getattr(wrfin.variables[varname], 
-#                                   "coordinates").split()
-#             lon_coord = coord_names[0]
-#             lat_coord = coord_names[1]
-#         else:
-#             lat_coord, lon_coord = get_coord_pairs(varname)
-#     else:
-#         lat_coord, lon_coord = latlon_coordvars(wrfin.variables)
-#     
-#     return (wrfin.variables[lat_coord][time_idx_or_slice,:],
-#             wrfin.variables[lon_coord][time_idx_or_slice,:],
-#             proj_params)
         
 
 def from_args(func, argnames, *args, **kwargs):
@@ -2936,18 +2935,31 @@ def psafilepath():
     return os.path.join(os.path.dirname(__file__), "data", "psadilookup.dat")
 
 
-def get_id(obj):
-    """Return the object id.
+def get_filepath(obj):
     
-    The object id is used as a caching key for various routines.  If the
+    try:
+        path = obj.filepath()
+    except AttributeError:
+        try:
+            path = obj.file.path
+        except:
+            raise ValueError("file contains no path information")
+    
+    return path
+
+def get_id(obj, prefix=''):
+    """Return the cache id.
+    
+    The cache id is used as a caching key for various routines. If the
     object type is a mapping, then the result will also be a 
-    mapping of each key to the object id for the value.  Otherwise, only the 
-    object id is returned.
+    mapping of each key to the object id for the value.
     
     Args:
     
         obj (:obj:`object`): Any object type.
         
+        prefix (:obj:`str`): A string to help with recursive calls.
+        
     Returns:
     
         :obj:`int` or :obj:`dict`: If the *obj* parameter is not a mapping, 
@@ -2955,12 +2967,18 @@ def get_id(obj):
         key to the object id for the value is returned.
     
     """
+    if not is_multi_file(obj):
+        return hash(prefix + get_filepath(obj))
+    
+    # For sequences, the hashing string will be the list ID and the 
+    # path for the first file in the sequence
     if not is_mapping(obj):
-        return id(obj)
+        _next = next(iter(obj))
+        return get_id(_next, prefix + str(id(obj)))
     
     # For each key in the mapping, recursively call get_id until
     # until a non-mapping is found
-    return {key : get_id(val) for key,val in viewitems(obj)}
+    return {key : get_id(val, prefix) for key,val in viewitems(obj)}
 
 
 def geo_bounds(var=None, wrfin=None, varname=None, timeidx=0, method="cat",

test/cachetest.py

@@ -2,7 +2,10 @@ from __future__ import (absolute_import, division, print_function,
                         unicode_literals)
 
 from threading import Thread
+try:
     from Queue import Queue
+except ImportError:
+    from queue import Queue
 from collections import OrderedDict
 
 import unittest as ut

test/generator_test.py

@@ -0,0 +1,17 @@
+from __future__ import (absolute_import, division, print_function, unicode_literals)
+
+from wrf import getvar
+from netCDF4 import Dataset as nc
+#ncfile = nc("/Users/ladwig/Documents/wrf_files/wrfout_d01_2016-02-25_18_00_00")
+ncfile = nc("/Users/ladwig/Documents/wrf_files/wrfout_d01_2016-10-07_00_00_00")
+
+def gen_seq():
+    wrfseq = [ncfile, ncfile, ncfile]
+    for wrf in wrfseq:
+        yield wrf
+        
+p_gen = getvar(gen_seq(), "P", method="join")
+
+print(p_gen)
+del p_gen
+