Merge pull request #44 from csyhuang/release_for_nhn_grl2022

Release v0.6.0

docs/source/conf.py

@@ -59,17 +59,17 @@
 
 # General information about the project.
 project = u'hn2016_falwa'
-copyright = u'2020, Clare S. Y. Huang'
+copyright = u'2022, Clare S. Y. Huang'
 author = u'Clare S. Y. Huang'
 
 # The version info for the project you're documenting, acts as replacement for
 # |version| and |release|, also used in various other places throughout the
 # built documents.
 #
 # The short X.Y version.
-version = u'0.5.0'
+version = u'0.6.0'
 # The full version, including alpha/beta/rc tags.
-release = u'0.5.0'
+release = u'0.6.0'
 
 # The language for content autogenerated by Sphinx. Refer to documentation
 # for a list of supported languages.

docs/source/index.rst

@@ -95,7 +95,7 @@ Issues Reporting
 
 Please make inquiries about / report issues / with the package and suggest feature extensions on the `Issues page <https://github.com/csyhuang/hn2016_falwa/issues>`_. 
 
-If you need help analyzing output from particular model/analysis with our techniques, feel free to email me *csyhuang@uchicago.edu* with sample datasets and I can configure the code for you.
+If you need help analyzing output from particular model/analysis with our techniques, feel free to email me *csyhuang@protonmail.com* with sample datasets and I can configure the code for you.
 
 
 Indices and tables

examples/nh2018_science/demo_script_for_nh2018.ipynb

@@ -601,7 +601,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.7.6"
+   "version": "3.7.10"
   }
  },
  "nbformat": 4,

hn2016_falwa/compute_flux_dirinv.f90

@@ -0,0 +1,132 @@
+SUBROUTINE compute_flux_dirinv(pv,uu,vv,pt,tn0,ts0,statn,stats,qref,uref,tref,fawa,ubar,tbar,&
+        imax, JMAX, kmax, nd, nnd, jb, jd, &
+        a, om, dz, h, rr, cp, prefac,&
+        astarbaro,ubaro,urefbaro,ua1baro,ua2baro,ep1baro,ep2baro,ep3baro,ep4,astar1,astar2)
+
+  INTEGER, INTENT(IN) :: imax, JMAX, kmax, nd, nnd, jb, jd
+  REAL, INTENT(IN) :: pv(imax,jmax,kmax),uu(imax,jmax,kmax),vv(imax,jmax,kmax),pt(imax,jmax,kmax),&
+          tn0(kmax),ts0(kmax),statn(kmax),stats(kmax),qref(nd,kmax),uref(jd,kmax),tref(jd,kmax),&
+          fawa(nd,kmax),ubar(nd,kmax),tbar(nd,kmax)
+  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)
+
+  REAL :: tb(kmax),tg(kmax)
+  REAL :: ua1(imax,nd),ua2(imax,nd),ep1(imax,nd)
+  REAL :: ep2(imax,nd),ep3(imax,nd)
+  REAL :: qe(imax,nd),ue(imax,nd)
+  REAL :: z(kmax)
+  REAL :: u(nd,kmax)
+
+  !a = 6378000.
+  pi = acos(-1.)
+  !om = 7.29e-5
+  dp = pi/float(jmax-1)
+  !dz = 500.
+  !h = 7000.
+  !r = 287.
+  rkappa = rr/cp
+  !prefac = 6745.348
+
+  ! *** Default values for boundary ***
+  !jb = 5
+  !jd = 86 ! nd - lower bounding latitude
+
+  do k = 1,kmax
+   z(k) = dz*float(k-1)
+  enddo
+
+
+! **** hemispheric-mean potential temperature ****
+  tg(:) = tn0(:)
+
+! **** wave activity and nonlinear zonal flux F2 ****
+
+  astarbaro(:,:) = 0.
+  ubaro(:,:) = 0.
+  urefbaro(:) = 0.
+  ua1baro(:,:) = 0.
+  ua2baro(:,:) = 0.
+  ep1baro(:,:) = 0.
+  ep2baro(:,:) = 0.
+  ep3baro(:,:) = 0.
+  ep4(:,:) = 0.
+  dc = dz/prefac
+
+  do k = 2,kmax-1
+    zk = dz*float(k-1)
+    do i = 1,imax
+      do j = jb+1,nd-1            ! 5N and higher latitude
+        astar1(i,j,k) = 0.       ! LWA*cos(phi)
+        astar2(i,j,k) = 0.       ! LWA*cos(phi)
+        ua2(i,j) = 0.          !F2
+        phi0 = dp*float(j-1)           !latitude
+        cor = 2.*om*sin(phi0)          !Coriolis parameter
+        do jj = 1,nd
+          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
+          aa = a*dp*cos(phi1)                      !length element
+          if((qe(i,jj).le.0.).and.(jj.ge.j)) then  !LWA*cos and F2
+            astar2(i,j,k)=astar2(i,j,k)-qe(i,jj)*aa  !anticyclonic
+            ua2(i,j) = ua2(i,j)-qe(i,jj)*ue(i,jj)*aa
+          endif
+          if((qe(i,jj).gt.0.).and.(jj.lt.j)) then
+            astar1(i,j,k)=astar1(i,j,k)+qe(i,jj)*aa  !cyclonic
+            ua2(i,j) = ua2(i,j)+qe(i,jj)*ue(i,jj)*aa
+          endif
+        enddo
+
+        !  ******** Other fluxes ******
+
+        ua1(i,j) = uref(j-jb,k)*(astar1(i,j,k) +  &
+        astar2(i,j,k))            !F1
+        ep1(i,j) = -0.5*(uu(i,j+nd-1,k)-uref(j-jb,k))**2  !F3a
+        ep1(i,j) = ep1(i,j)+0.5*vv(i,j+nd-1,k)**2    !F3a+b
+        ep11 = 0.5*(pt(i,j+nd-1,k)-tref(j-jb,k))**2         !F3c
+        zz = dz*float(k-1)
+        ep11 = ep11*(rr/h)*exp(-rkappa*zz/h)
+        ep11 = ep11*2.*dz/(tg(k+1)-tg(k-1))
+        ep1(i,j) = ep1(i,j)-ep11                   !F3
+        phip = dp*float(j)
+        cosp = cos(phip)          ! cosine for one grid north
+        phi0 = dp*float(j-1)
+        cos0 = cos(phi0)          ! cosine for latitude grid
+        sin0 = sin(phi0)          ! sine for latitude grid
+        phim = dp*float(j-2)
+        cosm = cos(phim)          ! cosine for one grid south
+        ep1(i,j) = ep1(i,j)*cos0 ! correct for cosine factor
+
+
+        ! meridional eddy momentum flux one grid north and south
+        ep2(i,j)=(uu(i,j+91,k)-uref(j-jb,k))*vv(i,j+91,k)*cosp*cosp
+        ep3(i,j)=(uu(i,j+89,k)-uref(j-jb,k))*vv(i,j+89,k)*cosm*cosm
+
+        ! low-level meridional eddy heat flux
+        if(k.eq.2) then     ! (26) of SI-HN17
+          ep41 = 2.*om*sin0*cos0*dz/prefac       ! prefactor
+          ep42 = exp(-dz/h)*vv(i,j+nd-1,2)*(pt(i,j+nd-1,2)-tref(j-jb,2))
+          ep42 = ep42/(tg(3)-tg(1))
+          ep43 = vv(i,j+nd-1,1)*(pt(i,j+nd-1,1)-tref(j-jb,1))
+          ep43 = 0.5*ep43/(tg(2)-tg(1))
+          ep4(i,j) = ep41*(ep42+ep43)   ! low-level heat flux
+        endif
+      enddo
+    enddo
+
+    ! ******** Column average: (25) of SI-HN17 ********
+
+    astarbaro(:,:) = astarbaro(:,:)+(astar1(:,:,k)    &
+    + astar2(:,:,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
+    ep2baro(:,:) = ep2baro(:,:)+ep2(:,:)*exp(-zk/h)*dc
+    ep3baro(:,:) = ep3baro(:,:)+ep3(:,:)*exp(-zk/h)*dc
+    do j = jb+1,nd  ! ### yet to be multiplied by cosine
+      ubaro(:,j) = ubaro(:,j)+uu(:,j+nd-1,k)*exp(-zk/h)*dc
+      urefbaro(j) = urefbaro(j)+uref(j-jb,k)*exp(-zk/h)*dc
+    enddo
+  enddo
+
+END SUBROUTINE

hn2016_falwa/compute_qref_and_fawa_first.f90

@@ -0,0 +1,181 @@
+SUBROUTINE compute_qref_and_fawa_first(pv, uu, vort, pt, tn0, ts0, statn, stats, imax, JMAX, kmax, nd, nnd, jb, jd, &
+        a, omega, dz, h, rr, cp, &
+        qref,u,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),&
+          stats(kmax),statn(kmax),ts0(kmax),tn0(kmax)
+  REAL, INTENT(OUT) :: qref(nd,kmax),u(jd,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)
+
+  !   **** take QG analysis and compute Q_ref and invert for U_ref & Theta_ref for NH (Direct solver) ***
+
+  !integer,parameter :: imax = 360, JMAX = 181, KMAX = 97
+  !integer,parameter :: nd = 91,nnd=181
+  !integer,parameter :: jb = 5   ! lower bounding latitude
+  !integer,parameter :: jd = 86  ! nd - lower bounding latitude
+
+  REAL :: pv2(imax,jmax)
+  REAL :: vort2(imax,jmax)
+  REAL :: qn(nnd),an(nnd),aan(nnd),tb(kmax),tg(kmax)
+  REAL :: cn(nnd),ccn(nnd),cref(nd,kmax)
+  REAL :: alat(nd),phi(nd),z(kmax),cbar(nd,kmax)
+  REAL :: tj(jd-2),rj(jd-2)
+  REAL :: qjj(jd-2,jd-2),cjj(jd-2,jd-2)
+  REAL :: xjj(jd-2,jd-2),yj(jd-2)
+  REAL :: djj(jd-2,jd-2),sjj(jd-2,jd-2)
+  REAL :: pjk(jd-2,kmax),pj(jd-2)
+  REAL :: qbar(nd,kmax)
+
+  pi = acos(-1.)
+  dp = pi/float(jmax-1)
+  rkappa = rr/cp
+
+  do nn = 1,nd
+    phi(nn) = dp*float(nn-1)
+    alat(nn) = 2.*pi*a*a*(1.-sin(phi(nn)))
+  enddo
+
+  do k = 1,kmax
+    z(k) = dz*float(k-1)
+  enddo
+
+
+  ! **** Zonal-mean field ****
+  do j = 91,jmax
+    qbar(j-90,:) = 0.
+    tbar(j-90,:) = 0.
+    ubar(j-90,:) = 0.
+    do i = 1,imax
+      qbar(j-90,:) = qbar(j-90,:)+pv(i,j,:)/float(imax)
+      tbar(j-90,:) = tbar(j-90,:)+pt(i,j,:)/float(imax)
+      ubar(j-90,:) = ubar(j-90,:)+uu(i,j,:)/float(imax)
+    enddo
+  enddo
+
+  ! **** hemispheric-mean potential temperature ****
+  tb(:) = tn0(:)
+
+  do k = 2,kmax-1
+    pv2(:,:) = pv(:,:,k)
+    vort2(:,:) = vort(:,:,k)
+
+    !  **** compute qref via area analysis ****
+    qmax = maxval(pv2)
+    qmin = minval(pv2)
+    dq = (qmax-qmin)/float(nnd-1)
+    qn(:) = 0.
+    an(:) = 0.
+    cn(:) = 0.
+    do nn = 1,nnd
+      qn(nn) = qmax - dq*float(nn-1)
+    enddo
+    do j = 1,jmax
+      phi0 = -0.5*pi+dp*float(j-1)
+      do i = 1,imax
+        ind = 1+int((qmax-pv2(i,j))/dq)
+        da = a*a*dp*dp*cos(phi0)
+        an(ind) = an(ind) + da
+        cn(ind) = cn(ind) + da*pv2(i,j)
+      enddo
+    enddo
+    aan(1) = 0.
+    ccn(1) = 0.
+    do nn = 2,nnd
+      aan(nn) = aan(nn-1)+an(nn)
+      ccn(nn) = ccn(nn-1)+cn(nn)
+    enddo
+    do j = 1,nd-1
+      do nn = 1,nnd-1
+        if(aan(nn).le.alat(j).and.aan(nn+1).gt.alat(j)) then
+          dd = (alat(j)-aan(nn))/(aan(nn+1)-aan(nn))
+          qref(j,k) = qn(nn)*(1.-dd)+qn(nn+1)*dd
+          cref(j,k) = ccn(nn)*(1.-dd)+ccn(nn+1)*dd
+        endif
+      enddo
+    enddo
+
+    qref(nd,k) = qmax
+
+    cbar(nd,k) = 0.
+    do j=nd-1,1,-1
+      phi0 = dp*(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)
+    enddo
+
+    ! **** compute Kelvin's circulation based on absolute vorticity (for b.c.) ****
+
+
+    qmax = maxval(vort2)
+    qmin = minval(vort2)
+    dq = (qmax-qmin)/float(nnd-1)
+    qn(:) = 0.
+    an(:) = 0.
+    cn(:) = 0.
+    do nn = 1,nnd
+      qn(nn) = qmax - dq*float(nn-1)
+    enddo
+    do j = 1,jmax
+      phi0 = -0.5*pi+dp*float(j-1)
+      do i = 1,imax
+      ind = 1+int((qmax-vort2(i,j))/dq)
+      da = a*a*dp*dp*cos(phi0)
+      an(ind) = an(ind) + da
+      cn(ind) = cn(ind) + da*vort2(i,j)
+      enddo
+    enddo
+    aan(1) = 0.
+    ccn(1) = 0.
+    do nn = 2,nnd
+      aan(nn) = aan(nn-1)+an(nn)
+      ccn(nn) = ccn(nn-1)+cn(nn)
+    enddo
+    do j = 1,nd-1
+      do nn = 1,nnd-1
+        if(aan(nn).le.alat(j).and.aan(nn+1).gt.alat(j)) then
+          dd = (alat(j)-aan(nn))/(aan(nn+1)-aan(nn))
+          ckref(j,k) = ccn(nn)*(1.-dd)+ccn(nn+1)*dd
+        endif
+      enddo
+    enddo
+
+  enddo
+
+  ! ***** normalize QGPV by sine (latitude)  ****
+
+  do j = 2,nd
+    phi0 = dp*float(j-1)
+    cor = sin(phi0)
+    qref(j,:) = qref(j,:)/cor
+  enddo
+
+  do k = 2,kmax-1
+    qref(1,k) = 2.*qref(2,k)-qref(3,k)
+  enddo
+
+  ! ***** FAWA *****
+  fawa(:,:) = (cref(:,:)-cbar(:,:))/(2.*pi*a)
+
+  !  ***** Direct solver to invert Q_ref *****
+
+  !  *** downward sweep ***
+
+  ! **** top boundary condition (Eqs. 24-25) *****
+  tjk(:,:) = 0.
+  sjk(:,:,:) = 0.
+  do jj = jb+2,90
+    j = jj-jb
+    phi0 = float(jj-1)*dp
+    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)
+    sjk(j-1,j-1,kmax-1) = 1.
+  enddo
+END

hn2016_falwa/interpolate_fields.f90

@@ -1,12 +1,12 @@
 SUBROUTINE interpolate_fields(nlon, nlat, nlev, kmax, uu, vv, temp, plev, height, &
-                              aa, omega, dz, hh, rr, cp, pv, ut, vt, theta, stat)
+                              aa, omega, dz, hh, rr, cp, pv, ut, vt, avort, theta, stat)
 
 
     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)
     REAL, INTENT(in) :: aa, omega, dz,hh, rr, cp
-    REAL, INTENT(out) :: pv(nlon,nlat,kmax), ut(nlon,nlat,kmax), vt(nlon,nlat,kmax)
+    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)
 
 
@@ -18,7 +18,7 @@ SUBROUTINE interpolate_fields(nlon, nlat, nlev, kmax, uu, vv, temp, plev, height
     REAL ::  xlon(nlon),ylat(nlat)
     REAL ::  t0(kmax),zlev(nlev)
     REAL ::  st(nlon,nlat),zmst(nlat)
-    REAL ::  avort(nlon,nlat,kmax),zmav(nlat,kmax)
+    REAL ::  zmav(nlat,kmax)
     REAL ::  zmpv(nlat,kmax)
     REAL :: rkappa, pi, dphi