Bugfix and hybrid coordinate example from Pragallva + additional contribution + readme files

CONTRIBUTING.md

@@ -0,0 +1,15 @@
+# How to contribute to this repo?
+
+1. Submit an issue to notify us with your plan in advance. 
+2. Fork this repo to your own GitHub account. The detailed logistics to fork the repo, track the original repository as a remote of the fork, etc., can be found in this comprehensive guide by [Jake Jarvis](https://jarv.is/): [How To Fork a GitHub Repository & Submit a Pull Request](https://jarv.is/notes/how-to-pull-request-fork-github/)
+3. Set this repo as the `upstream` remote and pull the latest commits to your own repo. 
+4. Create a new branch from `master` and make changes there.
+5. Before submitting a pull request, make sure you:
+   1. no longer have new commits to push,
+   2. run through all unit tests `pytest tests/` and make sure they all pass (if your changes cause some tests to fail, check and fix them), and
+   3. have merged all commits from `master` of this repo onto your branch.
+6. Please include the maintainers of this repo as reviewers of your pull request:
+   - Clare Huang `@csyhuang`
+   - Christopher Polster `@chpolste`
+
+ This page is subject to changes. Suggestions on improvement are welcome.

LICENSE.txt

@@ -1,4 +1,4 @@
-Copyright (c) 2015-2016 Clare S. Y. Huang
+Copyright (c) 2015 Clare S. Y. Huang
 
 Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:
 

docs/source/conf.py

@@ -78,9 +78,9 @@
 # built documents.
 #
 # The short X.Y version.
-version = u'1.1.0'
+version = u'1.2.0'
 # The full version, including alpha/beta/rc tags.
-release = u'1.1.0'
+release = u'1.2.0'
 
 # The language for content autogenerated by Sphinx. Refer to documentation
 # for a list of supported languages.

falwa/__init__.py

@@ -4,7 +4,7 @@
 Author: Clare Huang, Christopher Polster
 """
 
-__version__ = "1.1.0"
+__version__ = "1.2.0"
 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

falwa/f90_modules/compute_qref_and_fawa_first.f90

@@ -1,13 +1,13 @@
 SUBROUTINE compute_qref_and_fawa_first(pv, uu, vort, pt, tn0, imax, JMAX, kmax, nd, nnd, jb, jd, &
-        a, omega, dz, h, rr, cp, &
+        a, omega, dz, h, dphi, dlambda, rr, cp, &
         qref, ubar, tbar, fawa, ckref, tjk, sjk)
 
 
   !USE mkl95_LAPACK, ONLY: GETRF,GETRI
 
-  INTEGER, INTENT(IN) :: imax, JMAX, kmax, nd, nnd, jb, jd
-  REAL, INTENT(in) :: a, omega, dz, h, rr, cp
   REAL, INTENT(IN) :: pv(imax,jmax,kmax),uu(imax,jmax,kmax),vort(imax,jmax,kmax),pt(imax,jmax,kmax),tn0(kmax)
+  INTEGER, INTENT(IN) :: imax, JMAX, kmax, nd, nnd, jb, jd
+  REAL, INTENT(in) :: a, omega, dz, h, dphi, dlambda, rr, cp
   REAL, INTENT(OUT) :: qref(nd,kmax),ubar(nd,kmax),tbar(nd,kmax),fawa(nd,kmax),ckref(nd,kmax),&
           tjk(jd-2,kmax-1),sjk(jd-2,jd-2,kmax-1)
 
@@ -26,11 +26,12 @@ SUBROUTINE compute_qref_and_fawa_first(pv, uu, vort, pt, tn0, imax, JMAX, kmax,
   REAL :: qbar(nd,kmax)
 
   pi = acos(-1.)
-  dp = pi/float(jmax-1)
+  !dphi = pi/float(jmax-1)  !!!  This is dlat
+  !dlambda = 2*pi/float(imax)  !!!! pragallva- correction added on Dec 13, 2023. This is dlon
   rkappa = rr/cp
 
   do nn = 1,nd
-    phi(nn) = dp*float(nn-1)
+    phi(nn) = dphi*float(nn-1)
     alat(nn) = 2.*pi*a*a*(1.-sin(phi(nn)))
   enddo
 
@@ -69,10 +70,10 @@ SUBROUTINE compute_qref_and_fawa_first(pv, uu, vort, pt, tn0, imax, JMAX, kmax,
       qn(nn) = qmax - dq*float(nn-1)
     enddo
     do j = 1,jmax
-      phi0 = -0.5*pi+dp*float(j-1)
+      phi0 = -0.5*pi+dphi*float(j-1)
       do i = 1,imax
         ind = 1+int((qmax-pv2(i,j))/dq)
-        da = a*a*dp*dp*cos(phi0)
+        da = a*a*dphi*dlambda*cos(phi0)
         an(ind) = an(ind) + da
         cn(ind) = cn(ind) + da*pv2(i,j)
       enddo
@@ -97,9 +98,9 @@ SUBROUTINE compute_qref_and_fawa_first(pv, uu, vort, pt, tn0, imax, JMAX, kmax,
 
     cbar(nd,k) = 0.
     do j=nd-1,1,-1
-      phi0 = dp*(float(j)-0.5)
+      phi0 = dphi*(float(j)-0.5)
       cbar(j,k) = cbar(j+1,k)+0.5*(qbar(j+1,k)+qbar(j,k)) &
-      *a*dp*2.*pi*a*cos(phi0)
+      *a*dphi*2.*pi*a*cos(phi0)
     enddo
 
     ! **** compute Kelvin's circulation based on absolute vorticity (for b.c.) ****
@@ -115,10 +116,10 @@ SUBROUTINE compute_qref_and_fawa_first(pv, uu, vort, pt, tn0, imax, JMAX, kmax,
       qn(nn) = qmax - dq*float(nn-1)
     enddo
     do j = 1,jmax
-      phi0 = -0.5*pi+dp*float(j-1)
+      phi0 = -0.5*pi+dphi*float(j-1)
       do i = 1,imax
       ind = 1+int((qmax-vort2(i,j))/dq)
-      da = a*a*dp*dp*cos(phi0)
+      da = a*a*dphi*dlambda*cos(phi0)
       an(ind) = an(ind) + da
       cn(ind) = cn(ind) + da*vort2(i,j)
       enddo
@@ -143,7 +144,7 @@ SUBROUTINE compute_qref_and_fawa_first(pv, uu, vort, pt, tn0, imax, JMAX, kmax,
   ! ***** normalize QGPV by sine (latitude)  ****
 
   do j = 2,nd
-    phi0 = dp*float(j-1)
+    phi0 = dphi*float(j-1)
     cor = sin(phi0)
     qref(j,:) = qref(j,:)/cor
   enddo
@@ -164,12 +165,12 @@ SUBROUTINE compute_qref_and_fawa_first(pv, uu, vort, pt, tn0, imax, JMAX, kmax,
   sjk(:,:,:) = 0.
   do jj = jb+2,(nd-1)
     j = jj-jb
-    phi0 = float(jj-1)*dp
+    phi0 = float(jj-1)*dphi
     cos0 = cos(phi0)
     sin0 = sin(phi0)
     tjk(j-1,kmax-1) = -dz*rr*cos0*exp(-z(kmax-1)*rkappa/h)
     tjk(j-1,kmax-1) = tjk(j-1,kmax-1)*(tbar(j+1,kmax)-tbar(j-1,kmax))
-    tjk(j-1,kmax-1) = tjk(j-1,kmax-1)/(4.*omega*sin0*dp*h*a)
+    tjk(j-1,kmax-1) = tjk(j-1,kmax-1)/(4.*omega*sin0*dphi*h*a)
     sjk(j-1,j-1,kmax-1) = 1.
   enddo
 END

falwa/f90_modules/compute_reference_states.f90

@@ -1,10 +1,10 @@
 SUBROUTINE compute_reference_states(pv,uu,pt,stat,nlon,nlat,kmax,jd,npart,maxits,&
-               a,om,dz,eps,h,r,cp,rjac,qref,u,tref,num_of_iter)
+               a,om,dz,eps,h,dphi,dlambda,r,cp,rjac,qref,u,tref,num_of_iter)
 
 
     integer, intent(in) :: nlon,nlat,kmax,jd,npart,maxits
     real, intent(in) :: pv(nlon,nlat,kmax),uu(nlon,nlat,kmax),pt(nlon,nlat,kmax),stat(kmax)
-    real, intent(in) :: a, om, dz,eps,h, r, cp, rjac
+    real, intent(in) :: a, om, dz,eps, h, dphi, dlambda, r, cp, rjac
     real, intent(out) :: qref(jd,kmax),u(jd,kmax),tref(jd,kmax)
     integer, intent(out) :: num_of_iter
 
@@ -20,10 +20,11 @@ SUBROUTINE compute_reference_states(pv,uu,pt,stat,nlon,nlat,kmax,jd,npart,maxits
 
     real :: cref(jd,kmax)
     real :: qbar(jd,kmax),ubar(jd,kmax),tbar(jd,kmax)
-    real :: pi, dp, zero, half, qtr, one, rkappa
+    real :: pi, zero, half, qtr, one, rkappa
 
     pi = acos(-1.)
-    dp = pi/float(nlat-1)
+    !dphi = pi/float(nlat-1)
+    !dlambda = 2*pi/float(nlon)
     zero = 0.
     half = 0.5
     qtr = 0.25
@@ -32,7 +33,7 @@ SUBROUTINE compute_reference_states(pv,uu,pt,stat,nlon,nlat,kmax,jd,npart,maxits
 
 
     do nn = 1,jd
-        phi(nn) = dp*float(nn-1)
+        phi(nn) = dphi*float(nn-1)
         alat(nn) = 2.*pi*a*a*(1.-sin(phi(nn)))
     enddo
 
@@ -57,7 +58,7 @@ SUBROUTINE compute_reference_states(pv,uu,pt,stat,nlon,nlat,kmax,jd,npart,maxits
     tb(:) = 0.
     wt = 0.
     do j = jd,nlat
-        phi0 = dp*float(j-1)-0.5*pi
+        phi0 = dphi*float(j-1)-0.5*pi
         !tb(:) = tb(:)+cos(phi0)*tbar(j,:)
         tb(:) = tb(:)+cos(phi0)*tbar(j-(jd-1),:)
         wt = wt + cos(phi0)
@@ -80,10 +81,10 @@ SUBROUTINE compute_reference_states(pv,uu,pt,stat,nlon,nlat,kmax,jd,npart,maxits
             qn(nn) = qmax - dq*float(nn-1)
         enddo            
         do j = 1,nlat
-            phi0 = -0.5*pi+dp*float(j-1)
+            phi0 = -0.5*pi+dphi*float(j-1)
             do i = 1,nlon
                 ind = 1+int((qmax-pv2(i,j))/dq)
-                da = a*a*dp*dp*cos(phi0)
+                da = a*a*dphi*dlambda*cos(phi0)
                 an(ind) = an(ind) + da
                 cn(ind) = cn(ind) + da*pv2(i,j)
             enddo
@@ -108,9 +109,9 @@ SUBROUTINE compute_reference_states(pv,uu,pt,stat,nlon,nlat,kmax,jd,npart,maxits
 
         cbar(jd,k) = 0.
         do j=(jd-1),1,-1
-            phi0 = dp*(float(j)-0.5)
+            phi0 = dphi*(float(j)-0.5)
             cbar(j,k) = cbar(j+1,k)+0.5*(qbar(j+1,k)+qbar(j,k)) &
-            *a*dp*2.*pi*a*cos(phi0)
+            *a*dphi*2.*pi*a*cos(phi0)
         enddo 
 
     enddo
@@ -119,7 +120,7 @@ SUBROUTINE compute_reference_states(pv,uu,pt,stat,nlon,nlat,kmax,jd,npart,maxits
     ! ***** normalize QGPV by the Coriolis parameter ****
 
     do j = 2,jd
-        phi0 = dp*float(j-1)
+        phi0 = dphi*float(j-1)
         cor = 2.*om*sin(phi0)
         qref(j,:) = qref(j,:)/cor
     enddo
@@ -132,9 +133,9 @@ SUBROUTINE compute_reference_states(pv,uu,pt,stat,nlon,nlat,kmax,jd,npart,maxits
     ! **** SOR (elliptic solver a la Numerical Recipes) ****
 
     do j = 2,(jd-1)
-        phi0 = float(j-1)*dp 
-        phip = (float(j)-0.5)*dp
-        phim = (float(j)-1.5)*dp
+        phi0 = float(j-1)*dphi
+        phip = (float(j)-0.5)*dphi
+        phim = (float(j)-1.5)*dphi
         cos0 = cos(phi0)
         cosp = cos(phip)
         cosm = cos(phim)
@@ -146,7 +147,7 @@ SUBROUTINE compute_reference_states(pv,uu,pt,stat,nlon,nlat,kmax,jd,npart,maxits
             zm = 0.5*(z(k-1)+z(k))
             statp = 0.5*(stat(k+1)+stat(k))
             statm = 0.5*(stat(k-1)+stat(k))
-            fact = 4.*om*om*h*a*a*sin0*dp*dp/(dz*dz*r*cos0)
+            fact = 4.*om*om*h*a*a*sin0*dphi*dphi/(dz*dz*r*cos0)
             amp = exp(-zp/h)*exp(rkappa*zp/h)/statp
             amp = amp*fact*exp(z(k)/h)
             amm = exp(-zm/h)*exp(rkappa*zm/h)/statm
@@ -156,7 +157,7 @@ SUBROUTINE compute_reference_states(pv,uu,pt,stat,nlon,nlat,kmax,jd,npart,maxits
             cjk(j,k) = amp
             djk(j,k) = amm
             ejk(j,k) = -ajk(j,k)-bjk(j,k)-cjk(j,k)-djk(j,k)
-            fjk(j,k) = -om*a*dp*(qref(j+1,k)-qref(j-1,k))
+            fjk(j,k) = -om*a*dphi*(qref(j+1,k)-qref(j-1,k))
         enddo
     enddo
 
@@ -184,9 +185,9 @@ SUBROUTINE compute_reference_states(pv,uu,pt,stat,nlon,nlat,kmax,jd,npart,maxits
                 endif
             enddo
             u(j,1) = 0.
-            phi0 = dp*float(j-1)
+            phi0 = dphi*float(j-1)
             uz = dz*r*cos(phi0)*exp(-z(KMAX-1)*rkappa/h)
-            uz = uz*(tbar(j+1,KMAX-1)-tbar(j-1,KMAX-1))/(2.*om*sin(phi0)*dp*h*a)
+            uz = uz*(tbar(j+1,KMAX-1)-tbar(j-1,KMAX-1))/(2.*om*sin(phi0)*dphi*h*a)
             u(j,KMAX) = u(j,KMAX-2)-uz
         enddo
 
@@ -221,7 +222,7 @@ SUBROUTINE compute_reference_states(pv,uu,pt,stat,nlon,nlat,kmax,jd,npart,maxits
     num_of_iter = nnn
 
     do j = 2,(jd-1)
-        phi0 = dp*float(j-1)
+        phi0 = dphi*float(j-1)
         u(j,:) = u(j,:)/cos(phi0)
     enddo
     u(1,:) = ubar(1,:)
@@ -235,17 +236,17 @@ SUBROUTINE compute_reference_states(pv,uu,pt,stat,nlon,nlat,kmax,jd,npart,maxits
         tref(1,k) = t00
         tref(2,k) = t00
         do j = 2,(jd-1)
-            phi0 = dp*float(j-1)
+            phi0 = dphi*float(j-1)
             cor = 2.*om*sin(phi0)  
             uz = (u(j,k+1)-u(j,k-1))/(2.*dz)
             ty = -cor*uz*a*h*exp(rkappa*zz/h)
             ty = ty/r
-            tref(j+1,k) = tref(j-1,k)+2.*ty*dp
+            tref(j+1,k) = tref(j-1,k)+2.*ty*dphi
         enddo
         tg(k) = 0.
         wt = 0.
         do j = 1,jd
-            phi0 = dp*float(j-1)
+            phi0 = dphi*float(j-1)
             tg(k) = tg(k)+cos(phi0)*tref(j,k)
             wt = wt + cos(phi0)
         enddo

falwa/oopinterface.py

@@ -154,8 +154,8 @@ def __init__(self, xlon, ylat, plev, u_field, v_field, t_field, kmax=49, maxit=1
         self._nlat_analysis = self._ylat.size  # This is the number of latitude grid point used in analysis
 
         # === Coordinate-related ===
-        self.dphi = np.deg2rad(180./(self._nlat_analysis-1))
-        self.dlambda = np.deg2rad(self.xlon[1] - self.xlon[0])
+        self.dphi = np.deg2rad(180./(self._nlat_analysis-1))     # F90 code: dphi = pi/float(nlat-1)
+        self.dlambda = np.deg2rad(360./self.nlon)                # F90 code: dlambda = 2*pi/float(nlon)
         self.npart = npart if npart is not None else self._nlat_analysis
         self.kmax = kmax
         self.height = np.array([i * dz for i in range(kmax)])
@@ -947,21 +947,23 @@ def _compute_reference_state_wrapper(self, qgpv, u, theta):
             PTref(numpy.ndarray): Reference state of potential temperature of dimension [nlat, kmax]
         """
         return compute_reference_states(
-            qgpv,
-            u,
-            theta,
-            self._domain_average_storage.static_stability,
-            self.equator_idx,
-            self.npart,
-            self.maxit,
-            self.planet_radius,
-            self.omega,
-            self.dz,
-            self.tol,
-            self.scale_height,
-            self.dry_gas_constant,
-            self.cp,
-            self.rjac,
+            pv=qgpv,
+            uu=u,
+            pt=theta,
+            stat=self._domain_average_storage.static_stability,
+            jd=self.equator_idx,
+            npart=self.npart,
+            maxits=self.maxit,
+            a=self.planet_radius,
+            om=self.omega,
+            dz=self.dz,
+            eps=self.tol,
+            h=self.scale_height,
+            dphi=self.dphi,
+            dlambda=self.dlambda,
+            r=self.dry_gas_constant,
+            cp=self.cp,
+            rjac=self.rjac,
         )
 
     def _compute_intermediate_flux_terms(self):
@@ -1147,6 +1149,8 @@ def _compute_reference_states_nhn22_hemispheric_wrapper(self, qgpv, u, avort, th
             omega=self.omega,
             dz=self.dz,
             h=self.scale_height,
+            dphi=self.dphi,
+            dlambda=self.dlambda,
             rr=self.dry_gas_constant,
             cp=self.cp)