Merge pull request #128 from csyhuang/pr_branch_v2.1.0a

Pre-release v2.1.0a
csyhuang · Aug 30, 2024 · 640960e · 640960e
2 parents 57d5f45 + 5bfee5b
commit 640960e
Show file tree

Hide file tree

Showing 8 changed files with 114 additions and 32 deletions.
diff --git a/docs/source/conf.py b/docs/source/conf.py
@@ -76,9 +76,9 @@
 # built documents.
 #
 # The short X.Y version.
-version = u'2.0.0'
+version = u'2.1.0a'
 # The full version, including alpha/beta/rc tags.
-release = u'2.0.0'
+release = u'2.1.0a'
 
 # The language for content autogenerated by Sphinx. Refer to documentation
 # for a list of supported languages.
@@ -170,7 +170,7 @@
 #  dir menu entry, description, category)
 texinfo_documents = [
     (master_doc, 'falwa', u'falwa Documentation',
-     author, 'falwa', 'Python package to compute finite-amplitude wave activity.',
+     author, 'falwa', 'Python package to compute finite-amplitude local wave activity diagnostics.',
      'Miscellaneous'),
 ]
 

diff --git a/pyproject.toml b/pyproject.toml
@@ -8,10 +8,10 @@ authors = [
     { name="Clare S. Y. Huang", email="csyhuang@uchicago.edu" },
     { name="Christopher Polster", email="cpolster@uni-mainz.de" },
 ]
-description = "python package to compute finite-amplitude local wave activity (Huang and Nakamura 2016, JAS)"
+description = "Python package to compute finite-amplitude local wave activity diagnostics (Huang and Nakamura 2016, JAS)"
 readme = "readme.md"
 license = { file="LICENSE.txt" }
-version = "2.0.0"
+version = "2.1.0a"
 requires-python = ">=3.10"
 classifiers = [
     "License :: OSI Approved :: MIT License",

diff --git a/readme.md b/readme.md
@@ -1,4 +1,4 @@
-## Python Library: falwa (v2.0.0)
+## Python Library: falwa (v2.1.0a)
 
 [![Python package](https://github.com/csyhuang/hn2016_falwa/actions/workflows/python-build-test.yml/badge.svg)](https://github.com/csyhuang/hn2016_falwa/actions/workflows/python-build-test.yml)[![codecov.io](https://codecov.io/gh/csyhuang/hn2016_falwa/branch/master/graph/badge.svg)](https://codecov.io/gh/csyhuang/hn2016_falwa)[![Documentation Status](https://readthedocs.org/projects/hn2016-falwa/badge/?version=latest)](http://hn2016-falwa.readthedocs.io/en/latest/?badge=latest)[![DOI](https://zenodo.org/badge/63908662.svg)](https://zenodo.org/badge/latestdoi/63908662)
 

diff --git a/src/falwa/__init__.py b/src/falwa/__init__.py
@@ -4,7 +4,7 @@
 Author: Clare Huang, Christopher Polster
 """
 
-__version__ = "2.0.0"
+__version__ = "2.1.0a"
 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

diff --git a/src/falwa/data_storage.py b/src/falwa/data_storage.py
@@ -228,6 +228,7 @@ def __init__(self, pydim: Union[int, Tuple], fdim: Union[int, Tuple],
         self.convergence_zonal_advective_flux = np.zeros(self.pydim)
         self.divergence_eddy_momentum_flux = np.zeros(self.pydim)
         self.meridional_heat_flux = np.zeros(self.pydim)
+        self.ncforce_baro = np.zeros(self.pydim)
 
 
 class BarotropicFluxTermsStorage(DerivedQuantityStorage):
@@ -247,6 +248,7 @@ def __init__(self, pydim: Union[int, Tuple], fdim: Union[int, Tuple],
         self.ep2baro = np.zeros(self.fdim)
         self.ep3baro = np.zeros(self.fdim)
         self.ep4 = np.zeros(self.fdim)
+        self.ncforce_baro = np.zeros(self.fdim)
         self.u_baro = np.zeros(self.fdim)
         self.lwa_baro = np.zeros(self.fdim)  # This is barotropic LWA (astarbaro)
 
@@ -268,6 +270,9 @@ def __init__(self, pydim: Union[int, Tuple], fdim: Union[int, Tuple],
     ep4_nhem = NHemProperty("ep4", (1, 0))
     ep4_shem = SHemProperty("ep4", (1, 0))
 
+    ncforce_nhem = NHemProperty("ncforce_baro", (1, 0))
+    ncforce_shem = SHemProperty("ncforce_baro", (1, 0))
+
     u_baro_nhem = NHemProperty("u_baro", (1, 0))
     u_baro_shem = SHemProperty("u_baro", (1, 0))
 

diff --git a/src/falwa/f90_modules/compute_flux_dirinv.f90 b/src/falwa/f90_modules/compute_flux_dirinv.f90
@@ -1,18 +1,22 @@
-SUBROUTINE compute_flux_dirinv_nshem(pv,uu,vv,pt,tn0,qref,uref,tref,&
-        imax, JMAX, kmax, nd, jb, jd, is_nhem, &
-        a, om, dz, h, rr, cp, prefac,&
-        astarbaro,ubaro,urefbaro,ua1baro,ua2baro,ep1baro,ep2baro,ep3baro,ep4,astar1,astar2)
-
-  REAL, INTENT(IN) :: pv(imax,jmax,kmax),uu(imax,jmax,kmax),vv(imax,jmax,kmax),pt(imax,jmax,kmax),&
+SUBROUTINE compute_flux_dirinv_nshem( &
+    pv, uu, vv, pt, ncforce, tn0, qref, uref, tref, &
+    imax, JMAX, kmax, nd, jb, jd, is_nhem, &
+    a, om, dz, h, rr, cp, prefac,&
+    astarbaro,ubaro,urefbaro,ua1baro,ua2baro,ep1baro,ep2baro,ep3baro,ep4,astar1,astar2,ncforcebaro)
+
+  REAL, INTENT(IN) :: pv(imax,jmax,kmax),uu(imax,jmax,kmax),vv(imax,jmax,kmax),pt(imax,jmax,kmax), &
+          ncforce(imax, jmax, kmax), &
           tn0(kmax),qref(nd,kmax),uref(jd,kmax),tref(jd,kmax)
   INTEGER, INTENT(IN) :: imax, JMAX, kmax, nd, jb, jd
   LOGICAL, INTENT(IN) :: is_nhem
   REAL, INTENT(IN) :: a, om, dz, h, rr, cp, prefac
   REAL, INTENT(OUT) :: astarbaro(imax,nd),ubaro(imax,nd),urefbaro(nd),ua1baro(imax,nd),ua2baro(imax,nd),&
-          ep1baro(imax,nd),ep2baro(imax,nd),ep3baro(imax,nd),ep4(imax,nd),astar1(imax,nd,kmax),astar2(imax,nd,kmax)
+          ep1baro(imax,nd),ep2baro(imax,nd),ep3baro(imax,nd),ep4(imax,nd),astar1(imax,nd,kmax),astar2(imax,nd,kmax),&
+          ncforcebaro(imax,nd)
 
   REAL :: tg(kmax),ua1(imax,nd),ua2(imax,nd),ep1(imax,nd),ep2(imax,nd),ep3(imax,nd)
   REAL :: qe(imax,nd),ue(imax,nd)
+  REAL :: ncforce2d(imax,nd),ncforce3d(imax,nd,kmax)
   REAL :: z(kmax)
   REAL :: aa, ab
   INTEGER :: jstart, jend
@@ -39,9 +43,14 @@ SUBROUTINE compute_flux_dirinv_nshem(pv,uu,vv,pt,tn0,qref,uref,tref,&
 ! **** hemispheric-mean potential temperature ****
   tg(:) = tn0(:)
 
+! **** Computing ncforce (testing) ***
+  ncforce2d(:,:) = 0.0
+  ncforce3d(:,:,:) = 0.0
+
 ! **** wave activity and nonlinear zonal flux F2 ****
 
   astarbaro(:,:) = 0.
+  ncforcebaro(:,:) = 0.
   ubaro(:,:) = 0.
   urefbaro(:) = 0.
   ua1baro(:,:) = 0.
@@ -67,6 +76,7 @@ SUBROUTINE compute_flux_dirinv_nshem(pv,uu,vv,pt,tn0,qref,uref,tref,&
       do j = jstart, jend
         astar1(i,j,k) = 0.       ! LWA*cos(phi)
         astar2(i,j,k) = 0.       ! LWA*cos(phi)
+        ncforce3d(i,j,k) = 0.
         ua2(i,j) = 0.          !F2
         if (is_nhem) then        !latitude
           phi0 = dp*float(j-1)
@@ -79,10 +89,12 @@ SUBROUTINE compute_flux_dirinv_nshem(pv,uu,vv,pt,tn0,qref,uref,tref,&
           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
+            ncforce2d(i,jj) = ncforce(i,jj+nd-1,k)
             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
+            ncforce2d(i,jj) = ncforce(i,jj,k)
             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
@@ -92,6 +104,7 @@ SUBROUTINE compute_flux_dirinv_nshem(pv,uu,vv,pt,tn0,qref,uref,tref,&
             else  ! Southern Hemisphere
               astar1(i,j,k)=astar1(i,j,k)-qe(i,jj)*aa  !cyclonic
             endif
+            ncforce3d(i,j,k)=ncforce3d(i,j,k)-ncforce2d(i,jj)*aa
             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
@@ -100,6 +113,7 @@ SUBROUTINE compute_flux_dirinv_nshem(pv,uu,vv,pt,tn0,qref,uref,tref,&
             else  ! Southern Hemisphere
               astar2(i,j,k)=astar2(i,j,k)+qe(i,jj)*aa  !anticyclonic
             endif
+            ncforce3d(i,j,k)=ncforce3d(i,j,k)+ncforce2d(i,jj)*aa
             ua2(i,j) = ua2(i,j)+qe(i,jj)*ue(i,jj)*ab
           endif
         enddo
@@ -172,6 +186,7 @@ SUBROUTINE compute_flux_dirinv_nshem(pv,uu,vv,pt,tn0,qref,uref,tref,&
 
     astarbaro(:,:) = astarbaro(:,:)+(astar1(:,:,k)    &
     + astar2(:,:,k))*exp(-zk/h)*dc
+    ncforcebaro(:,:) = ncforcebaro(:,:) + ncforce3d(:,:,k)*exp(-zk/h)*dc
     ua1baro(:,:) = ua1baro(:,:)+ua1(:,:)*exp(-zk/h)*dc
     ua2baro(:,:) = ua2baro(:,:)+ua2(:,:)*exp(-zk/h)*dc
     ep1baro(:,:) = ep1baro(:,:)+ep1(:,:)*exp(-zk/h)*dc

diff --git a/src/falwa/f90_modules/compute_lwa_and_barotropic_fluxes.f90 b/src/falwa/f90_modules/compute_lwa_and_barotropic_fluxes.f90
@@ -1,22 +1,25 @@
 SUBROUTINE compute_lwa_and_barotropic_fluxes(nlon, nlat, kmax, jd, &
-                          pv,uu,vv,pt,qref,uref,tref,&
+                         pv, uu, vv, pt, ncforce, qref, uref, tref, &
                          a, om, dz, h, r, cp, prefactor, &
-                         astar,astarbaro,ua1baro,ubaro,ua2baro,ep1baro,ep2baro,ep3baro,ep4)
+                         astar, astarbaro, ua1baro, ubaro, ua2baro, ep1baro, ep2baro, ep3baro, ep4, ncforcebaro)
 
     implicit none
 
     integer, intent(in) :: nlon, nlat, kmax, jd
     real, intent(in) :: pv(nlon,nlat,kmax),uu(nlon,nlat,kmax),vv(nlon,nlat,kmax), &
-                        pt(nlon,nlat,kmax),qref(jd,kmax),uref(jd,kmax),tref(jd,kmax)
+                        pt(nlon,nlat,kmax),ncforce(nlon,nlat,kmax), &
+                        qref(jd,kmax),uref(jd,kmax),tref(jd,kmax)
+    ! Comment from Clare on 2023/10/10:
+    ! ncforce is an optional argument in QGField interface. If not required, simply put in zero array.
     real, intent(in) ::  a, om, dz, h, r, cp, prefactor
-    real, intent(out) :: astar(nlon,jd,kmax),astarbaro(nlon,jd),ua1baro(nlon,jd),ubaro(nlon,jd),&
-                         ua2baro(nlon,jd),ep1baro(nlon,jd),ep2baro(nlon,jd),ep3baro(nlon,jd),ep4(nlon,jd)
+    real, intent(out) :: astar(nlon,jd,kmax),astarbaro(nlon,jd),ua1baro(nlon,jd),ubaro(nlon,jd), &
+                         ua2baro(nlon,jd),ep1baro(nlon,jd),ep2baro(nlon,jd),ep3baro(nlon,jd),ep4(nlon,jd), &
+                         ncforcebaro(nlon,jd)
     ! === dummy variables ===
     integer :: i, j, jj, k
     real :: qe(nlon,nlat),ue(nlon,nlat)
-    real :: ua1(nlon,jd,kmax)
-    real :: ua2(nlon,jd,kmax),ep1(nlon,jd,kmax)
-    real :: ep2(nlon,jd,kmax),ep3(nlon,jd,kmax)
+    real :: ua1(nlon,jd,kmax), ua2(nlon,jd,kmax), ep1(nlon,jd,kmax), ep2(nlon,jd,kmax), ep3(nlon,jd,kmax)
+    real :: ncforce3d(nlon,jd,kmax)
     real :: tg(kmax)
     real :: aa, ab, cor, dc, ep11, ep41, ep42, ep43, wt, zk, zz, dp, rkappa
     real :: cosphi(0:jd+1), sinphi(0:jd+1)
@@ -45,12 +48,17 @@ SUBROUTINE compute_lwa_and_barotropic_fluxes(nlon, nlat, kmax, jd, &
     enddo
     tg(:) = tg(:)/wt       ! averaging
 
+    ! **** Computing ncforce (testing) ***
+    ncforcebaro(:,:) = 0.0
+    ncforce3d(:,:,:) = 0.0
+
     ! **** wave activity and nonlinear zonal flux F2 ****
 
     do k = 2,kmax-1
         do i = 1,nlon
             do j = 1,jd-1            ! 13.5N and higher latitude
                 astar(i,j,k) = 0.       ! LWA*cos(phi)
+                ncforce3d(i,j,k) = 0.
                 ua2(i,j,k) = 0.         ! F2
                 cor = 2.*om*sinphi(j)          !Coriolis parameter
                 ab = a*dp !constant length element
@@ -61,6 +69,7 @@ SUBROUTINE compute_lwa_and_barotropic_fluxes(nlon, nlat, kmax, jd, &
                     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
+                        ncforce3d(i,j,k)=ncforce3d(i,j,k)+ncforce(i,jj+jd-1,k)*aa
                         ua2(i,j,k) = ua2(i,j,k)+qe(i,jj)*ue(i,jj)*ab
                     endif  
                 enddo
@@ -71,6 +80,7 @@ SUBROUTINE compute_lwa_and_barotropic_fluxes(nlon, nlat, kmax, jd, &
                     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
+                        ncforce3d(i,j,k)=ncforce3d(i,j,k)-ncforce(i,jj+jd-1,k)*aa
                         ua2(i,j,k) = ua2(i,j,k)-qe(i,jj)*ue(i,jj)*ab
                     endif  
                 enddo
@@ -109,6 +119,7 @@ SUBROUTINE compute_lwa_and_barotropic_fluxes(nlon, nlat, kmax, jd, &
     ! ******** Column average: (25) of SI-HN17 ********
 
     astarbaro(:,:) = 0.
+    ncforcebaro(:,:) = 0.
     ubaro(:,:) = 0.
     ua1baro(:,:) = 0.
     ua2baro(:,:) = 0.
@@ -120,6 +131,7 @@ SUBROUTINE compute_lwa_and_barotropic_fluxes(nlon, nlat, kmax, jd, &
     do k = 2,kmax-1
         zk = dz*float(k-1)
         astarbaro(:,:) = astarbaro(:,:)+astar(:,:,k)*exp(-zk/h)*dc
+        ncforcebaro(:,:) = ncforcebaro(:,:)+ncforce3d(:,:,k)*exp(-zk/h)*dc
         ua1baro(:,:) = ua1baro(:,:)+ua1(:,:,k)*exp(-zk/h)*dc
         ua2baro(:,:) = ua2baro(:,:)+ua2(:,:,k)*exp(-zk/h)*dc
         ep1baro(:,:) = ep1baro(:,:)+ep1(:,:,k)*exp(-zk/h)*dc