Merge pull request #93 from csyhuang/pre-release-1.0.0a0

Pre release 1.0.0a0

docs/Makefile

@@ -4,7 +4,7 @@
 # You can set these variables from the command line.
 SPHINXOPTS    =
 SPHINXBUILD   = sphinx-build
-SPHINXPROJ    = hn2016_falwa
+SPHINXPROJ    = falwa
 SOURCEDIR     = source
 BUILDDIR      = build
 

docs/source/conf.py

@@ -1,6 +1,6 @@
 # -*- coding: utf-8 -*-
 #
-# hn2016_falwa documentation build configuration file, created by
+# falwa documentation build configuration file, created by
 # sphinx-quickstart on Mon Aug 21 13:58:59 2017.
 #
 # This file is execfile()d with the current directory set to its
@@ -20,14 +20,14 @@
 import sys
 # Need to fix the issue below
 sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)) + "/../../")
-sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)) + "/../../hn2016_falwa/")
+sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)) + "/../../falwa/")
 
 # import sphinx_rtd_theme
 
 
 autodoc_mock_imports = [
     "matplotlib",  # for plot_utils.py
-    "hn2016_falwa",
+    "falwa",
     "interpolate_fields",
     "compute_reference_states",
     "compute_lwa_and_barotropic_fluxes",
@@ -69,7 +69,7 @@
 master_doc = 'index'
 
 # General information about the project.
-project = u'hn2016_falwa'
+project = u'falwa'
 copyright = u'2022, Clare S. Y. Huang'
 author = u'Clare S. Y. Huang'
 
@@ -78,9 +78,9 @@
 # built documents.
 #
 # The short X.Y version.
-version = u'0.7.2'
+version = u'1.0.0a0'
 # The full version, including alpha/beta/rc tags.
-release = u'0.7.2'
+release = u'1.0.0a0'
 
 # The language for content autogenerated by Sphinx. Refer to documentation
 # for a list of supported languages.
@@ -150,7 +150,7 @@
 # (source start file, target name, title,
 #  author, documentclass [howto, manual, or own class]).
 latex_documents = [
-    (master_doc, 'hn2016_falwa.tex', u'hn2016\\_falwa Documentation',
+    (master_doc, 'falwa.tex', u'hn2016falwa Documentation',
      u'Shao Ying (Clare) Huang', 'manual'),
 ]
 
@@ -160,7 +160,7 @@
 # One entry per manual page. List of tuples
 # (source start file, name, description, authors, manual section).
 man_pages = [
-    (master_doc, 'hn2016_falwa', u'hn2016_falwa Documentation',
+    (master_doc, 'falwa', u'falwa Documentation',
      [author], 1)
 ]
 
@@ -171,8 +171,8 @@
 # (source start file, target name, title, author,
 #  dir menu entry, description, category)
 texinfo_documents = [
-    (master_doc, 'hn2016_falwa', u'hn2016_falwa Documentation',
-     author, 'hn2016_falwa', 'Python package to compute finite-amplitude wave activity.',
+    (master_doc, 'falwa', u'falwa Documentation',
+     author, 'falwa', 'Python package to compute finite-amplitude wave activity.',
      'Miscellaneous'),
 ]
 

docs/source/index.rst

@@ -1,14 +1,14 @@
-.. hn2016_falwa documentation master file, created by
+.. falwa documentation master file, created by
    sphinx-quickstart on Mon Aug 21 13:58:59 2017.
    You can adapt this file completely to your liking, but it should at least
    contain the root `toctree` directive.
    
-hn2016_falwa: Finite-amplitude local wave activity
+falwa: Finite-amplitude local wave activity
 ==================================================
 
 .. image:: https://github.com/csyhuang/csyhuang.github.io/raw/master/assets/img/hn2016_falwa_diagram.png
 
-`hn2016_falwa <https://github.com/csyhuang/hn2016_falwa>`_ is a python library that provides tools to measure and study life cycle of large-scale
+`falwa <https://github.com/csyhuang/hn2016_falwa>`_ is a python library that provides tools to measure and study life cycle of large-scale
 extreme weather events. It implements the finite-amplitude local wave activity and flux diagnostic introduced in:
 
 - `Huang and Nakamura (2016, JAS) <https://journals.ametsoc.org/doi/abs/10.1175/JAS-D-15-0194.1>`_
@@ -26,15 +26,15 @@ Since the package is still being actively developed, please use the *develop* mo
 To install the package for the first time, clone the GitHub repo and install via `develop` mode:::
 
  git clone https://github.com/csyhuang/hn2016_falwa.git
- cd hn2016_falwa
+ cd falwa
  python setup.py develop
 
 
 To incorporate updates, pull the new version of the code from GitHub. Remove any existing f2py modules and recompile.::
 
- # Assume you are already in the hn2016_falwa/ repo
+ # Assume you are already in the falwa/ repo
  git pull
- rm hn2016_falwa/*.so
+ rm falwa/*.so
  python setup.py develop
  pytest tests/ # to check if the package can be run properly
 

hn2016_falwa/__init__.py → falwa/__init__.py

@@ -4,7 +4,7 @@
 Author: Clare Huang, Christopher Polster
 """
 
-__version__ = "0.7.2"
+__version__ = "1.0.0a0"
 from .interpolate_fields import interpolate_fields
 from .interpolate_fields_direct_inv import interpolate_fields_direct_inv
 from .compute_qref_and_fawa_first import compute_qref_and_fawa_first

hn2016_falwa/barotropic_field.py → falwa/barotropic_field.py

@@ -7,7 +7,7 @@
 
 import numpy as np
 
-from hn2016_falwa import basis
+from falwa import basis
 
 
 class BarotropicField(object):

hn2016_falwa/f90_modules/compute_flux_dirinv.f90 → falwa/f90_modules/compute_flux_dirinv.f90

@@ -74,30 +74,30 @@ SUBROUTINE compute_flux_dirinv_nshem(pv,uu,vv,pt,tn0,qref,uref,tref,&
           phi0 = dp*float(j-1)-0.5*pi
         endif
         cor = 2.*om*sin(phi0)          !Coriolis parameter
-        ab = a*dp*cos(phi0)
+        ab = a*dp  !constant length element
         do jj = 1,nd
-          if (is_nhem) then
+          if (is_nhem) then  ! Northern Hemisphere
             phi1 = dp*float(jj-1)
             qe(i,jj) = pv(i,jj+nd-1,k)-qref(j,k)    !qe; Q = qref
-            ue(i,jj) = uu(i,jj+nd-1,k)-uref(j-jb,k) !ue; shift uref 5N
-          else
+            ue(i,jj) = uu(i,jj+nd-1,k)*cos(phi0)-uref(j-jb,k)*cos(phi1) !ue; shift uref 5N
+          else  ! Southern Hemisphere
             phi1 = dp*float(jj-1)-0.5*pi
             qe(i,jj) = pv(i,jj,k)-qref(j,k)     !qe; Q = qref
-            ue(i,jj) = uu(i,jj,k)-uref(j-jb,k)  !ue;
+            ue(i,jj) = uu(i,jj,k)*cos(phi0)-uref(j,k)*cos(phi1)  !ue;
           endif
           aa = a*dp*cos(phi1)   !cosine factor in the meridional integral
           if((qe(i,jj).le.0.).and.(jj.ge.j)) then  !LWA*cos and F2
-            if (is_nhem) then
+            if (is_nhem) then  ! Northern Hemisphere
               astar2(i,j,k)=astar2(i,j,k)-qe(i,jj)*aa  !anticyclonic
-            else
+            else  ! Southern Hemisphere
               astar1(i,j,k)=astar1(i,j,k)-qe(i,jj)*aa  !cyclonic
             endif
             ua2(i,j) = ua2(i,j)-qe(i,jj)*ue(i,jj)*ab
           endif
           if((qe(i,jj).gt.0.).and.(jj.lt.j)) then
-            if (is_nhem) then
+            if (is_nhem) then  ! Northern Hemisphere
               astar1(i,j,k)=astar1(i,j,k)+qe(i,jj)*aa  !cyclonic
-            else
+            else  ! Southern Hemisphere
               astar2(i,j,k)=astar2(i,j,k)+qe(i,jj)*aa  !anticyclonic
             endif
             ua2(i,j) = ua2(i,j)+qe(i,jj)*ue(i,jj)*ab

hn2016_falwa/f90_modules/compute_lwa_and_barotropic_fluxes.f90 → falwa/f90_modules/compute_lwa_and_barotropic_fluxes.f90

@@ -53,11 +53,11 @@ SUBROUTINE compute_lwa_and_barotropic_fluxes(nlon, nlat, kmax, jd, &
                 astar(i,j,k) = 0.       ! LWA*cos(phi)
                 ua2(i,j,k) = 0.         ! F2
                 cor = 2.*om*sinphi(j)          !Coriolis parameter
-                ab = a*dp*cosphi(j)
+                ab = a*dp !constant length element
                 ! South of the current latitude
                 do jj = 1,j-1
                     qe(i,jj) = pv(i,jj+jd-1,k)-qref(j,k)*cor   !qe; Q = qref*cor
-                    ue(i,jj) = uu(i,jj+jd-1,k)-uref(j,k)       !ue
+                    ue(i,jj) = uu(i,jj+jd-1,k)*cosphi(j)-uref(j,k)*cosphi(jj)       !ue
                     aa = a*dp*cosphi(jj)                       !length element
                     if(qe(i,jj).gt.0.) then                    !LWA*cos(phi) and F2
                         astar(i,j,k)=astar(i,j,k)+qe(i,jj)*aa
@@ -67,7 +67,7 @@ SUBROUTINE compute_lwa_and_barotropic_fluxes(nlon, nlat, kmax, jd, &
                 ! North of the current latitude
                 do jj = j,jd
                     qe(i,jj) = pv(i,jj+jd-1,k)-qref(j,k)*cor   !qe; Q = qref*cor
-                    ue(i,jj) = uu(i,jj+jd-1,k)-uref(j,k)       !ue
+                    ue(i,jj) = uu(i,jj+jd-1,k)*cosphi(j)-uref(j,k)*cosphi(jj)       !ue
                     aa = a*dp*cosphi(jj)                       !length element
                     if(qe(i,jj).le.0.) then                    !LWA*cos(phi) and F2
                         astar(i,j,k)=astar(i,j,k)-qe(i,jj)*aa

hn2016_falwa/oopinterface.py → falwa/oopinterface.py

@@ -11,13 +11,13 @@
 from scipy.interpolate import interp1d
 from scipy.linalg.lapack import dgetrf, dgetri
 
-from hn2016_falwa import utilities
-from hn2016_falwa.constant import P_GROUND, SCALE_HEIGHT, CP, DRY_GAS_CONSTANT, EARTH_RADIUS, EARTH_OMEGA
-from hn2016_falwa.data_storage import InterpolatedFieldsStorage, DomainAverageStorage, ReferenceStatesStorage, \
+from falwa import utilities
+from falwa.constant import P_GROUND, SCALE_HEIGHT, CP, DRY_GAS_CONSTANT, EARTH_RADIUS, EARTH_OMEGA
+from falwa.data_storage import InterpolatedFieldsStorage, DomainAverageStorage, ReferenceStatesStorage, \
     LWAStorage, BarotropicFluxTermsStorage, OutputBarotropicFluxTermsStorage
 
 # *** Import f2py modules ***
-from hn2016_falwa import interpolate_fields, interpolate_fields_direct_inv, compute_qref_and_fawa_first,\
+from falwa import interpolate_fields, interpolate_fields_direct_inv, compute_qref_and_fawa_first,\
     matrix_b4_inversion, matrix_after_inversion, upward_sweep, compute_flux_dirinv_nshem, compute_reference_states,\
     compute_lwa_and_barotropic_fluxes
 from collections import namedtuple
@@ -139,12 +139,12 @@ def __init__(self, xlon, ylat, plev, u_field, v_field, t_field, kmax=49, maxit=1
 
         # === Do Interpolation on latitude grid if needed ===
         if self.need_latitude_interpolation:
-            interp_u = interp1d(self.ylat_no_equator, u_field, axis=1, fill_value="extrapolate")
-            interp_v = interp1d(self.ylat_no_equator, v_field, axis=1, fill_value="extrapolate")
-            interp_t = interp1d(self.ylat_no_equator, t_field, axis=1, fill_value="extrapolate")
-            self.u_field = interp_u(self.ylat)
-            self.v_field = interp_v(self.ylat)
-            self.t_field = interp_t(self.ylat)
+            interp_u = interp1d(self._input_ylat, u_field, axis=1, fill_value="extrapolate")
+            interp_v = interp1d(self._input_ylat, v_field, axis=1, fill_value="extrapolate")
+            interp_t = interp1d(self._input_ylat, t_field, axis=1, fill_value="extrapolate")
+            self.u_field = interp_u(self._ylat)
+            self.v_field = interp_v(self._ylat)
+            self.t_field = interp_t(self._ylat)
         else:
             self.u_field = u_field
             self.v_field = v_field
@@ -283,25 +283,29 @@ def _check_and_flip_ylat(self, ylat):
         # Check if ylat is in ascending order and include the equator
         if np.diff(ylat)[0] < 0:
             raise TypeError("ylat must be in ascending order")
+        # Save ylat input by user first
+        self._input_ylat = ylat
         if (ylat.size % 2 == 0) & (sum(ylat == 0.0) == 0):
             # Even grid
             self.need_latitude_interpolation = True
-            self.ylat_no_equator = ylat
-            self.ylat = np.linspace(-90., 90., ylat.size+1, endpoint=True)
+            self._ylat = np.linspace(-90., 90., ylat.size+1, endpoint=True)
             self.equator_idx = \
-                np.argwhere(self.ylat == 0)[0][0] + 1
+                np.argwhere(self._ylat == 0)[0][0] + 1
             # Fortran indexing starts from 1
         elif sum(ylat == 0) == 1:
             # Odd grid
             self.need_latitude_interpolation = False
-            self.ylat_no_equator = None
-            self.ylat = ylat
+            self._ylat = ylat
             self.equator_idx = np.argwhere(ylat == 0)[0][0] + 1 # Fortran indexing starts from 1
         else:
             raise TypeError(
                 "There are more than 1 grid point with latitude 0."
             )
-        self.clat = np.abs(np.cos(np.deg2rad(self.ylat)))
+        self.clat = np.abs(np.cos(np.deg2rad(self._ylat)))
+
+    @property
+    def ylat(self):
+        return self._input_ylat
 
     @staticmethod
     def _check_dimension_of_fields(field, field_name, expected_dim):
@@ -319,7 +323,7 @@ def _interp_back(self, field, interp_from, interp_to, which_axis=1):
         Private function to interpolate the results from odd grid to even grid.
         If the initial input to the QGField object is an odd grid, error will be raised.
         """
-        if self.ylat_no_equator is None:
+        if self._input_ylat is None:
             raise TypeError("No need for such interpolation.")
         else:
             return interp1d(
@@ -348,11 +352,11 @@ def _return_interp_variables(self, variable, interp_axis):
         if self.need_latitude_interpolation:
             if self.northern_hemisphere_results_only:
                 return self._interp_back(
-                    variable, self.ylat[-(self.nlat//2+1):],
-                    self.ylat_no_equator[-(self.nlat//2):],
+                    variable, self._ylat[-(self.nlat//2+1):],
+                    self._input_ylat[-(self.nlat // 2):],
                     which_axis=interp_axis)
             else:
-                return self._interp_back(variable, self.ylat, self.ylat_no_equator, which_axis=interp_axis)
+                return self._interp_back(variable, self._ylat, self._input_ylat, which_axis=interp_axis)
         else:
             return variable
 
@@ -646,8 +650,8 @@ def ylat_ref_states(self) -> np.array:
         Latitude dimension of reference state
         """
         if self.northern_hemisphere_results_only:
-            return self.ylat[-(self.nlat//2+1):]
-        return self.ylat
+            return self._ylat[-(self.nlat//2+1):]
+        return self._ylat
 
     @property
     def northern_hemisphere_results_only(self) -> bool:
@@ -829,10 +833,7 @@ def get_latitude_dim(self):
         """
         Return the latitude dimension of the input data.
         """
-        if self.need_latitude_interpolation:
-            return self.ylat_no_equator.size
-        else:
-            return self.nlat
+        return self._input_ylat.size
 
 
 class QGFieldNH18(QGFieldBase):
@@ -1226,7 +1227,7 @@ def _compute_intermediate_flux_terms(self):
 
         # Turn qref back to correct unit
 
-        ylat_input = self.ylat[-self.equator_idx:] if self.northern_hemisphere_results_only else self.ylat
+        ylat_input = self._ylat[-self.equator_idx:] if self.northern_hemisphere_results_only else self._ylat
         qref_correct_unit = self._reference_states_storage.qref_correct_unit(
             ylat=ylat_input, omega=self.omega, python_indexing=False)
 

hn2016_falwa/wrapper.py → falwa/wrapper.py

@@ -4,8 +4,8 @@
 Author: Clare Huang
 """
 import numpy as np
-from hn2016_falwa.constant import EARTH_RADIUS
-from hn2016_falwa import basis
+from falwa.constant import EARTH_RADIUS
+from falwa import basis
 
 
 def barotropic_eqlat_lwa(ylat, vort, area, dmu, n_points, planet_radius=EARTH_RADIUS):