Merge pull request #105 from csyhuang/release1.2.0

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)