PR for release 1.3 (#123)

* test squash merge * leave ncforce for next release * add compute_lwa_only * add docstring for nonconvergent_uref * Attempt to fix sphinx build * retry with sphinx build on readthedoc * add jinja<3.1 into .readthedocs.yaml * specify jinja2 version in doc/requirements.txt * convert back to python3.9 for sphinx build * update sphinx version * Update requirements.txt * fix documentation format * try specifying explicit package version * change pip3 back to pip * change sequence of installation to test again * update build status badge * add test coverage for QGField.compute_lwa_only * test for qgfield.nonconvergent_uref * add docstring * remove test branch from pr * Fix misformed link in docs * Remove duplicate _DataArrayCollectors --------- Co-authored-by: Christopher Polster <[email protected]>
csyhuang · Jul 15, 2024 · bb95b87 · bb95b87
1 parent e7a57dd
commit bb95b87
Show file tree

Hide file tree

Showing 41 changed files with 9,753 additions and 1,383 deletions.
diff --git a/.github/workflows/python-build-test.yml b/.github/workflows/python-build-test.yml
@@ -5,7 +5,7 @@ name: Python package
 
 on:
   push:
-    branches: [ "master", "draft_branch_for_v1.3c" ]
+    branches: [ "master" ]
   pull_request:
     branches: [ "master" ]
 
@@ -31,10 +31,10 @@ jobs:
         python -m pip install -U pytest
         python -m pip install --upgrade pip setuptools wheel
         sudo apt-get install gfortran
-        python -m pip install numpy==1.22.3
-        python -m pip install scipy
-        python -m pip install xarray
-        python setup.py develop        
+        python -m pip install scipy==1.11.2
+        python -m pip install xarray==2024.1.1
+        python -m pip install numpy==1.24.4
+        python setup.py develop
     - name: Test with pytest
       run: |
         coverage run -m pytest -k "not output_results"

diff --git a/developers/run_through_notebooks.sh b/developers/run_through_notebooks.sh
@@ -0,0 +1,7 @@
+jupyter nbconvert --execute --to notebook --inplace ../notebooks/nh2018_science/demo_script_for_nh2018.ipynb
+jupyter nbconvert --execute --to notebook --inplace ../notebooks/nh2018_science/demo_script_for_nh2018_with_xarray.ipynb
+jupyter nbconvert --execute --to notebook --inplace ../notebooks/nh2018_science/nhn22_reference_states.ipynb
+jupyter nbconvert --execute --to notebook --inplace ../notebooks/sandy/sandy_lwa_plots.ipynb
+jupyter nbconvert --execute --to notebook --inplace ../notebooks/simple/BarotropicField_example.ipynb
+jupyter nbconvert --execute --to notebook --inplace ../notebooks/simple/Example_barotropic.ipynb
+jupyter nbconvert --execute --to notebook --inplace ../notebooks/simple/Example_qgpv.ipynb
diff --git a/docs/source/Plot Utilities.rst b/docs/source/Plot Utilities.rst
@@ -1,10 +1,10 @@
-
-
 Plot Utilities
 ===============
 
 .. warning::
-    This module uses `matplotlib`, which is not included in the dependencies of this package.
+    This module uses `matplotlib` and `cartopy`, which are not included in the dependencies of this package.
 
 .. automodule:: plot_utils
    :members:
+
+
diff --git a/docs/source/Preprocessing.rst b/docs/source/Preprocessing.rst
@@ -0,0 +1,11 @@
+
+
+Preprocessing
+==============
+
+Since `QGField` requires all grid points of the input fields to have finite values, to process model output with grid points below topography masked as `NaN`, users have to fill in reasonable values for those grid points.
+
+A pragmatic implementation was published in `Lubis et al (2018) <https://journals.ametsoc.org/view/journals/clim/31/10/jcli-d-17-0382.1.xml>`_, in which the grids are filled level-by-level by solving Poisson equation, such that no extremum values are filled in. (In contrast, extrapolation over vertical direction would introduce unrealistic extreme values.) Such gridfill procedure is implemented by the python package `gridfill <https://github.com/ajdawson/gridfill>`_ published by Andrew Dawson. After computing the finite-amplitude wave diagnostics, regions with gridfilled values are masked out.
+
+.. automodule:: preprocessing
+   :members:
diff --git a/docs/source/Statistics Utilities.rst b/docs/source/Statistics Utilities.rst
@@ -0,0 +1,7 @@
+
+
+Statistics Utilities
+====================
+
+.. automodule:: stat_utils
+   :members:
diff --git a/docs/source/conf.py b/docs/source/conf.py
@@ -28,6 +28,7 @@
 autodoc_mock_imports = [
     "matplotlib",  # for plot_utils.py
     "falwa",
+    "cartopy",
     "interpolate_fields",
     "compute_reference_states",
     "compute_lwa_and_barotropic_fluxes",
@@ -78,9 +79,9 @@
 # built documents.
 #
 # The short X.Y version.
-version = u'1.2.1'
+version = u'1.3.0'
 # The full version, including alpha/beta/rc tags.
-release = u'1.2.1'
+release = u'1.3.0'
 
 # The language for content autogenerated by Sphinx. Refer to documentation
 # for a list of supported languages.

diff --git a/docs/source/index.rst b/docs/source/index.rst
@@ -73,16 +73,18 @@ Modules
 .. toctree::
    :maxdepth: 2
 
+   Example Notebooks
    Object Oriented Interface
+   Data Storage
    Xarray Interface
+   Preprocessing
    Barotropic Field
-   Basis Functions
-   Plot Utilities
    Utility Functions
+   Plot Utilities
+   netCDF Utilities
+   Statistics Utilities
+   Basis Functions
    Wrapper Functions
-   Data Storage
-   Example Notebooks
-
 
 Indices and tables
 ==================

diff --git a/docs/source/netCDF Utilities.rst b/docs/source/netCDF Utilities.rst
@@ -0,0 +1,7 @@
+
+
+netCDF Utilities
+=================
+
+.. automodule:: netcdf_utils
+   :members:
diff --git a/environment.yml b/environment.yml
@@ -11,6 +11,10 @@ dependencies:
 - pytest=7.4.0
 - xarray=2023.2.0
 - netCDF4=1.5.8
+- gridfill
 - jupyter
 - matplotlib
 - twine
+- cartopy
+# Other useful packages to visualize plots in notebooks
+# - cartopy
diff --git a/falwa/__init__.py b/falwa/__init__.py
@@ -4,13 +4,14 @@
 Author: Clare Huang, Christopher Polster
 """
 
-__version__ = "1.2.1"
-from .interpolate_fields import interpolate_fields
-from .interpolate_fields_direct_inv import interpolate_fields_direct_inv
+__version__ = "1.3.0"
+from .compute_qgpv import compute_qgpv
+from .compute_qgpv_direct_inv import compute_qgpv_direct_inv
 from .compute_qref_and_fawa_first import compute_qref_and_fawa_first
 from .matrix_b4_inversion import matrix_b4_inversion
 from .matrix_after_inversion import matrix_after_inversion
 from .upward_sweep import upward_sweep
 from .compute_reference_states import compute_reference_states
 from .compute_flux_dirinv import compute_flux_dirinv_nshem
 from .compute_lwa_and_barotropic_fluxes import compute_lwa_and_barotropic_fluxes
+from .compute_lwa_only_nhn22 import compute_lwa_only_nhn22
diff --git a/falwa/data_storage.py b/falwa/data_storage.py
@@ -1,7 +1,7 @@
 """
 ------------------------------------------
 File name: data_storage.py
-Author: Clare Huang
+Author: Clare Huang, Christopher Polster
 """
 from typing import Tuple, Optional, Union, List, NamedTuple
 import numpy as np

diff --git a/falwa/f90_modules/compute_lwa_only_nhn22.f90 b/falwa/f90_modules/compute_lwa_only_nhn22.f90
@@ -0,0 +1,107 @@
+SUBROUTINE compute_lwa_only_nhn22(pv,uu,qref, &
+        imax, JMAX, kmax, nd, jb, is_nhem, &
+        a, om, dz, h, rr, cp, prefac, &
+        astarbaro, ubaro, astar1, astar2)
+  ! This was from compute_flux_dirinv.f90
+
+  REAL, INTENT(IN) :: pv(imax,jmax,kmax),uu(imax,jmax,kmax),qref(nd,kmax)
+  INTEGER, INTENT(IN) :: imax, JMAX, kmax, nd, jb
+  LOGICAL, INTENT(IN) :: is_nhem
+  REAL, INTENT(IN) :: a, om, dz, h, rr, cp, prefac
+  REAL, INTENT(OUT) :: astarbaro(imax, nd), ubaro(imax, nd), astar1(imax,nd,kmax), astar2(imax,nd,kmax)
+
+  REAL :: qe(imax,nd)
+  REAL :: z(kmax)
+  REAL :: aa, ab
+  INTEGER :: jstart, jend
+
+  pi = acos(-1.)
+  dp = pi/float(jmax-1)
+  rkappa = rr/cp
+
+  ! *** Default values for boundary ***
+  !jb = 5
+  !jd = 86 ! nd - lower bounding latitude
+
+  do k = 1,kmax
+   z(k) = dz*float(k-1)
+  enddo
+
+
+! **** wave activity and nonlinear zonal flux F2 ****
+
+  astarbaro(:,:) = 0.
+  ubaro(:,:) = 0.
+  dc = dz/prefac
+
+  ! Bounds of y-indices in N/SHem
+  if (is_nhem) then  ! 5N and higher latitude
+    jstart = jb+1  ! 6
+    jend = nd-1
+  else
+    jstart = 2
+    jend = nd-jb   ! nd - 5
+  endif
+
+  do k = 2,kmax-1
+    zk = dz*float(k-1)
+    do i = 1,imax
+      do j = jstart, jend
+        astar1(i,j,k) = 0.       ! LWA*cos(phi)
+        astar2(i,j,k) = 0.       ! LWA*cos(phi)
+        if (is_nhem) then        !latitude
+          phi0 = dp*float(j-1)
+        else
+          phi0 = dp*float(j-1)-0.5*pi
+        endif
+        cor = 2.*om*sin(phi0)          !Coriolis parameter
+        ab = a*dp  !constant length element
+        do jj = 1,nd
+          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
+          else  ! Southern Hemisphere
+            phi1 = dp*float(jj-1)-0.5*pi
+            qe(i,jj) = pv(i,jj,k)-qref(j,k)     !qe; Q = qref
+          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  ! Northern Hemisphere
+              astar2(i,j,k)=astar2(i,j,k)-qe(i,jj)*aa  !anticyclonic
+            else  ! Southern Hemisphere
+              astar1(i,j,k)=astar1(i,j,k)-qe(i,jj)*aa  !cyclonic
+            endif
+          endif
+          if((qe(i,jj).gt.0.).and.(jj.lt.j)) then
+            if (is_nhem) then  ! Northern Hemisphere
+              astar1(i,j,k)=astar1(i,j,k)+qe(i,jj)*aa  !cyclonic
+            else  ! Southern Hemisphere
+              astar2(i,j,k)=astar2(i,j,k)+qe(i,jj)*aa  !anticyclonic
+            endif
+          endif
+        enddo
+
+      enddo
+      phip = dp*jb
+      phi0 = dp*(jb-1)
+      cosp = cos(phip)          ! cosine for one grid north
+      cos0 = cos(phi0)          ! cosine for latitude grid
+    enddo
+
+    ! ******** Column average: (25) of SI-HN17 ********
+
+    astarbaro(:,:) = astarbaro(:,:)+(astar1(:,:,k)    &
+    + astar2(:,:,k))*exp(-zk/h)*dc
+
+    if (is_nhem) then
+      do j = jstart,jend
+        ubaro(:,j) = ubaro(:,j)+uu(:,nd-1+j,k)*exp(-zk/h)*dc
+      enddo
+    else
+      do j = jstart,jend
+        ubaro(:,j) = ubaro(:,j)+uu(:,j,k)*exp(-zk/h)*dc
+      enddo
+    endif
+  enddo
+
+END SUBROUTINE
diff --git a/falwa/f90_modules/interpolate_fields.f90 → falwa/f90_modules/compute_qgpv.f90 b/falwa/f90_modules/interpolate_fields.f90 → falwa/f90_modules/compute_qgpv.f90
@@ -1,19 +1,18 @@
-SUBROUTINE interpolate_fields(nlon, nlat, nlev, kmax, uu, vv, temp, plev, height, &
-                              aa, omega, dz, hh, rr, cp, pv, ut, vt, avort, theta, stat)
+SUBROUTINE compute_qgpv(nlon, nlat, kmax, ut, vt, theta, height, t0, stat, &
+                        aa, omega, dz, hh, rr, cp, &
+                        pv, avort)
 
 
-    INTEGER, INTENT(IN) :: nlon, nlat, nlev, kmax
-    REAL, INTENT(IN) :: uu(nlon,nlat,nlev), vv(nlon,nlat,nlev), temp(nlon,nlat,nlev), &
-                        plev(nlev), height(kmax)
+    INTEGER, INTENT(IN) :: nlon, nlat, kmax
+    REAL, INTENT(IN) :: ut(nlon,nlat,kmax), vt(nlon,nlat,kmax), theta(nlon,nlat,kmax), &
+                        height(kmax), t0(kmax), stat(kmax)
     REAL, INTENT(in) :: aa, omega, hh, rr, cp
-    REAL, INTENT(out) :: pv(nlon,nlat,kmax), ut(nlon,nlat,kmax), vt(nlon,nlat,kmax), avort(nlon,nlat,kmax)
-    REAL, INTENT(out) :: theta(nlon,nlat,kmax), stat(kmax)
+    REAL, INTENT(out) :: pv(nlon,nlat,kmax), avort(nlon,nlat,kmax)
 
 
-    REAL ::  tt(nlon,nlat,nlev),tz(nlat,kmax),tzd(nlat,kmax)
+    REAL ::  tz(nlat,kmax),tzd(nlat,kmax)
     REAL ::  uz(nlat,kmax),uzd(nlat,kmax)
     REAL ::  vz(nlat,kmax),vzd(nlat,kmax)
-    REAL ::  t0(kmax),zlev(nlev)
     REAL ::  st(nlon,nlat),zmst(nlat)
     REAL ::  zmav(nlat,kmax)
     REAL ::  zmpv(nlat,kmax)
@@ -24,42 +23,6 @@ SUBROUTINE interpolate_fields(nlon, nlat, nlev, kmax, uu, vv, temp, plev, height
     pi = acos(-1.)
     dphi = pi/float(nlat-1)
 
-
-    ! ====  vertical interpolation ====
-    do k = 1,nlev
-        zlev(k) = -hh*alog(plev(k)/1000.)
-    enddo
-
-    do i = 1,nlon
-        do j = 1,nlat
-
-            do k = 1,nlev
-                tt(i,j,k) = temp(i,j,k)*((1000./plev(k))**rkappa)
-            enddo
-
-
-
-            st(i,j) = tt(i,j,1)      ! surface pot. temp
-
-            do kk = 1,kmax   ! vertical interpolation 
-                ttt = height(kk)
-                do k = 1,nlev-1
-                    tt2 = zlev(k+1)
-                    tt1 = zlev(k)
-                    if((ttt.ge.tt1).and.(ttt.lt.tt2)) then
-                        dd1 = (ttt-tt1)/(tt2-tt1)
-                        dd2 = 1.-dd1
-                        theta(i,j,kk) = tt(i,j,k)*dd2 + tt(i,j,k+1)*dd1
-                        ut(i,j,kk) = uu(i,j,k)*dd2 + uu(i,j,k+1)*dd1
-                        vt(i,j,kk) = vv(i,j,k)*dd2 + vv(i,j,k+1)*dd1
-                        goto 348
-                    endif
-                enddo
-                348 continue
-            enddo
-        enddo
-    enddo
-
     !  **** compute zonal mean ****
     tz = 0.
     uz = 0.
@@ -83,30 +46,7 @@ SUBROUTINE interpolate_fields(nlon, nlat, nlev, kmax, uu, vv, temp, plev, height
         enddo
     enddo
 
-
-    ! reference theta
-    do kk = 1,kmax
-        t0(kk) = 0.
-        csm = 0.
-        do j = 1,nlat
-            phi0 = -90.+float(j-1)*180./float(nlat-1)
-            phi0 = phi0*pi/180.
-            t0(kk) = t0(kk) + tzd(j,kk)*cos(phi0)
-            csm = csm + cos(phi0)
-        enddo
-        t0(kk) = t0(kk)/csm
-    enddo
-
-    ! static stability
-    do kk = 2,kmax-1
-        stat(kk) = (t0(kk+1)-t0(kk-1))/(height(kk+1)-height(kk-1))
-    enddo
-    stat(kmax) = (t0(kmax)-t0(kmax-1))/(height(kmax)-height(kmax-1))
-    stat(1) = (t0(2)-t0(1))/(height(2)-height(1))
-
-
     ! surface temp
-
     do j = 1,nlat
         zmst(j) = 0.
         do i = 1,nlon
@@ -194,12 +134,4 @@ SUBROUTINE interpolate_fields(nlon, nlat, nlev, kmax, uu, vv, temp, plev, height
         enddo
     enddo
 
-
-!     write(6,*) 'U ',maxval(uz),minval(uz)
-!     write(6,*) 'V ',maxval(vz),minval(vz)
-!     write(6,*) 'T ',maxval(tz),minval(tz)    
-!     write(6,*) 'AVORT',maxval(avort),minval(avort)    
-!     write(6,*) 'PV ',maxval(zmpv),minval(zmpv)
-!     write(6,*) 'zmpv', sum(zmpv, dim=1)
-
 END SUBROUTINE