diff --git a/docs/source/conf.py b/docs/source/conf.py index 6b79b68..40ec938 100644 --- a/docs/source/conf.py +++ b/docs/source/conf.py @@ -59,7 +59,7 @@ # 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 @@ -67,9 +67,9 @@ # 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. diff --git a/docs/source/index.rst b/docs/source/index.rst index 63e1fad..b3f24d1 100644 --- a/docs/source/index.rst +++ b/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 `_. -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 diff --git a/examples/nh2018_science/demo_script_for_nh2018.ipynb b/examples/nh2018_science/demo_script_for_nh2018.ipynb index 0566516..6e82a6f 100644 --- a/examples/nh2018_science/demo_script_for_nh2018.ipynb +++ b/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, diff --git a/hn2016_falwa/compute_flux_dirinv.f90 b/hn2016_falwa/compute_flux_dirinv.f90 new file mode 100644 index 0000000..6a3813d --- /dev/null +++ b/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 diff --git a/hn2016_falwa/compute_qref_and_fawa_first.f90 b/hn2016_falwa/compute_qref_and_fawa_first.f90 new file mode 100644 index 0000000..2fe9e16 --- /dev/null +++ b/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 diff --git a/hn2016_falwa/interpolate_fields.f90 b/hn2016_falwa/interpolate_fields.f90 index 9e6a625..b42f3c2 100644 --- a/hn2016_falwa/interpolate_fields.f90 +++ b/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 diff --git a/hn2016_falwa/interpolate_fields_dirinv.f90 b/hn2016_falwa/interpolate_fields_dirinv.f90 new file mode 100644 index 0000000..116dd60 --- /dev/null +++ b/hn2016_falwa/interpolate_fields_dirinv.f90 @@ -0,0 +1,226 @@ +SUBROUTINE interpolate_fields_direct_inv(nlon, nlat, nlev, kmax, jd, uu, vv, tt, plev, & + aa, omega, dz, hh, rr, cp, & + pv, uq, vq, avort, tq, statn, stats, tn0, ts0) + + + INTEGER, INTENT(IN) :: nlon, nlat, nlev, kmax + REAL, INTENT(IN) :: uu(nlon,nlat,nlev), vv(nlon,nlat,nlev), tt(nlon,nlat,nlev), & + plev(nlev) + REAL, INTENT(in) :: aa, omega, dz,hh, rr, cp + REAL, INTENT(out) :: pv(nlon,nlat,kmax), uq(nlon,nlat,kmax), vq(nlon,nlat,kmax), avort(nlon,nlat,kmax) + REAL, INTENT(out) :: tq(nlon,nlat,kmax), statn(kmax), stats(kmax), tn0(kmax), ts0(kmax) + + real :: tzd(nlat,kmax) + real :: xlon(nlon),ylat(nlat) + real :: height(kmax) + real :: zlev(nlev) + real :: st(nlon,nlat),zmst(nlat) + real :: tt0(nlon,nlat,kmax) + real :: zmav(nlat,kmax) + real :: zmpv(nlat,kmax) + integer :: dsadata,mm(12),inverse(12,nlev) + integer :: k0(kmax),kp(kmax) + real :: dd2(kmax),dd1(kmax),pks(kmax) + character*5 :: yy + character*4 :: y0(44),y00 + character*3 :: mn(12) + character*8 :: yr + + rkappa = rr/cp + pi = acos(-1.) + dphi = pi/float(nlat-1) + +! ====== Assign pseudoheight ===== + + do k = 1,kmax + height(k) = float(k-1)*dz + pks(k) = exp(rkappa*height(k)/hh) + enddo + + do k = 1,nlev + zlev(k) = -hh*alog(plev(k)/1000.) + enddo + + do kk = 2,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 + k0(kk) = k + kp(kk) = k+1 + dd1(kk) = (ttt-tt1)/(tt2-tt1) + dd2(kk) = 1.-dd1(kk) + endif + enddo + enddo + +! ==== vertical interpolation ==== + + do i = 1,nlon + do j = 1,nlat + + st(i,j) = tt(i,j,1) ! surface pot. temp + + do kk = 2,kmax ! vertical interpolation + uq(i,j,kk) = uu(i,j,k0(kk))*dd2(kk) + uu(i,j,kp(kk))*dd1(kk) + vq(i,j,kk) = vv(i,j,k0(kk))*dd2(kk) + vv(i,j,kp(kk))*dd1(kk) + ! wq(i,j,kk) = ww(i,j,k0(kk))*dd2(kk) + ww(i,j,kp(kk))*dd1(kk) + tq(i,j,kk) = tt(i,j,k0(kk))*dd2(kk) + tt(i,j,kp(kk))*dd1(kk) + tq(i,j,kk) = tq(i,j,kk)*pks(kk) ! potential temperature + ! zq(i,j,kk) = zz(i,j,k0(kk))*dd2(kk) + zz(i,j,kp(kk))*dd1(kk) + enddo + + tq(i,j,1) = tt(i,j,1) + uq(i,j,1) = uu(i,j,1) + vq(i,j,1) = vv(i,j,1) + ! wq(i,j,1) = ww(i,j,1) + ! zq(i,j,1) = zz(i,j,1) + enddo + enddo + +! **** compute zonal mean **** + + tzd = 0. + + do j = 1,nlat + do k = 1,kmax + do i = 1,nlon + tzd(j,k) = tzd(j,k) + tq(i,j,k)/float(nlon) + enddo + enddo + enddo + + + ! reference theta + do kk = 1,kmax + ts0(kk) = 0. + tn0(kk) = 0. + csm = 0. + cnm = 0. + do j = 1,jd + phi0 = -90.+float(j-1) + phi0 = phi0*pi/float(nlat-1) + ts0(kk) = ts0(kk) + tzd(j,kk)*cos(phi0) + csm = csm + cos(phi0) + enddo + ts0(kk) = ts0(kk)/csm + do j = jd,nlat + phi0 = -90.+float(j-1) + phi0 = phi0*pi/float(nlat-1) + tn0(kk) = tn0(kk) + tzd(j,kk)*cos(phi0) + cnm = cnm + cos(phi0) + enddo + tn0(kk) = tn0(kk)/cnm + enddo + + ! static stability + do kk = 2,kmax-1 + stats(kk) = (ts0(kk+1)-ts0(kk-1))/(height(kk+1)-height(kk-1)) + statn(kk) = (tn0(kk+1)-tn0(kk-1))/(height(kk+1)-height(kk-1)) + enddo + stats(kmax) = 2.*stats(kmax-1)-stats(kmax-2) + statn(kmax) = 2.*statn(kmax-1)-statn(kmax-2) + stats(1) = 2.*stats(2)-stats(3) + statn(1) = 2.*statn(2)-statn(3) + + ! surface temp + + do j = 1,nlat + zmst(j) = 0. + do i = 1,nlon + zmst(j) = zmst(j) + st(i,j)/float(nlon) + enddo + enddo + +! interior abs. vort + + do kk = 1,kmax + do j = 2,nlat-1 + phi0 = -90.+float(j-1) + phi0 = phi0*pi/float(nlat-1) + phim = -90.+float(j-2) + phim = phim*pi/float(nlat-1) + phip = -90.+float(j) + phip = phip*pi/float(nlat-1) + + do i = 2,359 + av1 = 2.*omega*sin(phi0) + av2 = (vq(i+1,j,kk)-vq(i-1,j,kk))/(2.*aa*cos(phi0)*dphi) + av3 = -(uq(i,j+1,kk)*cos(phip)-uq(i,j-1,kk)*cos(phim))/(2.*aa*cos(phi0)*dphi) + avort(i,j,kk) = av1+av2+av3 + enddo + + av1 = 2.*omega*sin(phi0) + av2 = (vq(2,j,kk)-vq(nlon,j,kk))/(2.*aa*cos(phi0)*dphi) + av3 = -(uq(1,j+1,kk)*cos(phip)-uq(1,j-1,kk)*cos(phim))/(2.*aa*cos(phi0)*dphi) + avort(1,j,kk) = av1+av2+av3 + av4 = 2.*omega*sin(phi0) + av5 = (vq(1,j,kk)-vq(359,j,kk))/(2.*aa*cos(phi0)*dphi) + av6 = & + -(uq(nlon,j+1,kk)*cos(phip)-uq(nlon,j-1,kk)*cos(phim))/(2.*aa*cos(phi0)*dphi) + avort(nlon,j,kk) = av4+av5+av6 + enddo + + avs = 0. + avn = 0. + do i = 1,nlon + avs = avs + avort(i,2,kk)/float(nlon) + avn = avn + avort(i,nlat-1,kk)/float(nlon) + enddo + avort(:,1,kk) = avs + avort(:,nlat,kk) = avn + enddo + + ! zonal mean vort + + do kk = 1,kmax + do j = 1,nlat + zmav(j,kk) = 0. + do i = 1,nlon + zmav(j,kk) = zmav(j,kk)+avort(i,j,kk)/float(nlon) + enddo + enddo + enddo + + ! interior pv + + do kk = 2,kmax-1 + do j = 1,nlat + phi0 = -90.+float(j-1) + phi0 = phi0*pi/float(nlat-1) + f = 2.*omega*sin(phi0) + if (j .le. jd) then + statp = stats(kk+1) + statm = stats(kk-1) + t00p = ts0(kk+1) + t00m = ts0(kk-1) + else + statp = statn(kk+1) + statm = statn(kk-1) + t00p = tn0(kk+1) + t00m = tn0(kk-1) + endif + + do i = 1,nlon + thetap = tq(i,j,kk+1) + thetam = tq(i,j,kk-1) + altp = exp(-height(kk+1)/hh)*(thetap-t00p)/statp + altm = exp(-height(kk-1)/hh)*(thetam-t00m)/statm + strc = (altp-altm)*f/(height(kk+1)-height(kk-1)) + pv(i,j,kk) = avort(i,j,kk) + exp(height(kk)/hh)*strc + enddo + enddo + enddo + +! zonal mean pv + + do kk = 1,kmax + do j = 1,nlat + zmpv(j,kk) = 0. + do i = 1,nlon + zmpv(j,kk) = zmpv(j,kk)+pv(i,j,kk)/float(nlon) + enddo + enddo + enddo +END SUBROUTINE diff --git a/hn2016_falwa/matrix_after_inversion.f90 b/hn2016_falwa/matrix_after_inversion.f90 new file mode 100644 index 0000000..f126b5b --- /dev/null +++ b/hn2016_falwa/matrix_after_inversion.f90 @@ -0,0 +1,39 @@ +SUBROUTINE matrix_after_inversion(k,kmax,jb,jd,qjj,djj,cjj,tj,rj,sjk,tjk) + + INTEGER, INTENT(in) :: k, kmax, jb, jd + REAL, INTENT(in) :: qjj(jd-2,jd-2),djj(jd-2,jd-2),cjj(jd-2,jd-2),rj(jd-2) + REAL, INTENT(INOUT) :: sjk(jd-2,jd-2,kmax-1),tjk(jd-2,kmax-1),tj(jd-2) ! Note that tj is not used in subsequent modules + + integer :: i, j + real :: xjj(jd-2,jd-2),yj(jd-2) + + + do i = 1,jd-2 + do j = 1,jd-2 + xjj(i,j) = 0. + do kk = 1,jd-2 + xjj(i,j) = xjj(i,j)+qjj(i,kk)*djj(kk,j) + enddo + sjk(i,j,k-1) = -xjj(i,j) + enddo + enddo + + ! **** Evaluate rk - Ck Tk **** + do i = 1,jd-2 + yj(i) = 0. + do kk = 1,jd-2 + yj(i) = yj(i)+cjj(i,kk)*tj(kk) + enddo + yj(i) = rj(i)-yj(i) + enddo + + ! ***** Evaluate Eq. 23 ******* + do i = 1,jd-2 + tj(i) = 0. + do kk = 1,jd-2 + tj(i) = tj(i)+qjj(i,kk)*yj(kk) + enddo + tjk(i,k-1) = tj(i) + enddo + +END SUBROUTINE matrix_after_inversion diff --git a/hn2016_falwa/matrix_b4_inversion.f90 b/hn2016_falwa/matrix_b4_inversion.f90 new file mode 100644 index 0000000..cb32def --- /dev/null +++ b/hn2016_falwa/matrix_b4_inversion.f90 @@ -0,0 +1,94 @@ +SUBROUTINE matrix_b4_inversion(k,jmax,kmax,nd,jb,jd,z,statn,qref,ckref,& + a, om, dz, h, rr, cp, & + qjj,djj,cjj,rj,tj,u,sjk,tjk) + + integer, INTENT(in) :: k, jmax, kmax, nd, jb, jd + REAL, INTENT(in) :: z(kmax),statn(kmax),qref(nd,kmax),ckref(nd,kmax) + REAL, INTENT(in) :: a, om, dz, h, rr, cp + REAL, INTENT(OUT) :: qjj(jd-2,jd-2),djj(jd-2,jd-2),cjj(jd-2,jd-2),rj(jd-2),tj(jd-2) + REAL, INTENT(INOUT) :: u(jd,kmax),sjk(jd-2,jd-2,kmax-1),tjk(jd-2,kmax-1) + + REAL :: stats(kmax),ts0(kmax),tn0(kmax) + REAL :: qn(jmax),an(jmax),aan(jmax),tb(kmax),tg(kmax) + REAL :: cn(jmax),ccn(jmax),tref(jd,kmax) + REAL :: alat(nd),phi(nd),cbar(nd,kmax) + REAL :: xjj(jd-2,jd-2) + REAL :: sjj(jd-2,jd-2) + + rkappa = rr/cp + pi = acos(-1.) + dp = pi/float(jmax-1) + + zp = 0.5*(z(k+1)+z(k)) + zm = 0.5*(z(k-1)+z(k)) + statp = 0.5*(statn(k+1)+statn(k)) + statm = 0.5*(statn(k-1)+statn(k)) + cjj(:,:) = 0. + djj(:,:) = 0. + qjj(:,:) = 0. + sjj(:,:) = sjk(:,:,k) + tj(:) = tjk(:,k) + do jj = jb+2,90 + j = jj - jb + phi0 = float(jj-1)*dp + phip = (float(jj)-0.5)*dp + phim = (float(jj)-1.5)*dp + cos0 = cos(phi0) + cosp = cos(phip) + cosm = cos(phim) + sin0 = sin(phi0) + sinp = sin(phip) + sinm = sin(phim) + + fact = 4.*om*om*h*a*a*sin0*dp*dp/(dz*dz*rr*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 + amm = amm*fact*exp(z(k)/h) + + ! ***** Specify A, B, C, D, E, F (Eqs. 4-9) ***** + ajk = 1./(sinp*cosp) + bjk = 1./(sinm*cosm) + cjk = amp + djk = amm + ejk = ajk+bjk+cjk+djk + fjk = -0.5*a*dp*(qref(jj+1,k)-qref(jj-1,k)) + + ! ***** Specify rk (Eq. 15) **** + + ! **** North-south boundary conditions **** + u(jd,k) = 0. + phi0 = dp*float(jb) + ! u(1,k) = ubar(jb+1,k)*cos(phi0) + u(1,k) = ckref(jb+1,k)/(2.*pi*a)-om*a*cos(phi0) + + rj(j-1) = fjk + if(j.eq.2) rj(j-1) = fjk - bjk*u(1,k) + if(j.eq.jd-1) rj(j-1) = fjk - ajk*u(jd,k) + + ! ***** Specify Ck & Dk (Eqs. 18-19) ***** + cjj(j-1,j-1) = cjk + djj(j-1,j-1) = djk + + ! **** Specify Qk (Eq. 17) ***** + qjj(j-1,j-1) = -ejk + if(j-1.ge.1.and.j-1.lt.jd-2) then + qjj(j-1,j) = ajk + endif + if(j-1.gt.1.and.j-1.le.jd-2) then + qjj(j-1,j-2) = bjk + endif + enddo + + ! **** Compute Qk + Ck Sk ******* + do i = 1,jd-2 + do j = 1,jd-2 + xjj(i,j) = 0. + do kk = 1,jd-2 + xjj(i,j) = xjj(i,j)+cjj(i,kk)*sjj(kk,j) + enddo + qjj(i,j) = qjj(i,j)+xjj(i,j) + enddo + enddo + +END SUBROUTINE matrix_b4_inversion diff --git a/hn2016_falwa/oopinterface.py b/hn2016_falwa/oopinterface.py index 0686740..97c0cf8 100644 --- a/hn2016_falwa/oopinterface.py +++ b/hn2016_falwa/oopinterface.py @@ -2,10 +2,18 @@ import warnings import numpy as np from scipy.interpolate import interp1d +from scipy.linalg.lapack import dgetrf, dgetri from hn2016_falwa import utilities from hn2016_falwa.constant import * 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 +from matrix_b4_inversion import matrix_b4_inversion +from matrix_after_inversion import matrix_after_inversion +from upward_sweep import upward_sweep +from compute_flux_dirinv import compute_flux_dirinv + from compute_reference_states import compute_reference_states from compute_lwa_and_barotropic_fluxes import compute_lwa_and_barotropic_fluxes from collections import namedtuple @@ -63,7 +71,10 @@ class QGField(object): planet_radius : float, optional Radius of the planet in meters. Default = 6.378e+6 (Earth's radius). - + eq_boundary_index: int, optional + The improved inversion algorithm of reference states allow modification of equatorward boundary + to be the absolute vorticity. This parameter specify the location of grid point (from equator) + which will be used as boundary. Default = 5. Examples -------- @@ -73,7 +84,7 @@ class QGField(object): def __init__(self, xlon, ylat, plev, u_field, v_field, t_field, kmax=49, maxit=100000, dz=1000., npart=None, tol=1.e-5, rjac=0.95, scale_height=SCALE_HEIGHT, cp=CP, dry_gas_constant=DRY_GAS_CONSTANT, - omega=EARTH_OMEGA, planet_radius=EARTH_RADIUS, prefactor=None): + omega=EARTH_OMEGA, planet_radius=EARTH_RADIUS, prefactor=None, eq_boundary_index=5): """ Create a QGField object. @@ -153,6 +164,7 @@ def __init__(self, xlon, ylat, plev, u_field, v_field, t_field, kmax=49, maxit=1 self.dz = dz self.tol = tol self.rjac = rjac + self.eq_boundary_index = eq_boundary_index # === Constants === self.scale_height = scale_height @@ -174,8 +186,17 @@ def __init__(self, xlon, ylat, plev, u_field, v_field, t_field, kmax=49, maxit=1 self._qgpv_temp = None self._interpolated_u_temp = None self._interpolated_v_temp = None + self._interpolated_avort_temp = None self._interpolated_theta_temp = None self._static_stability = None + self._static_stability_n = None + self._static_stability_s = None + self._tn0 = None + self._ts0 = None + + # === Temporary solution for direct inv routine + self._f_qref, self._f_u, self._f_tref, self._f_ubar, self._f_tbar, self._f_fawa, self._f_ckref = \ + None, None, None, None, None, None, None # Computation from computer_reference_states self._qref_stemp = None @@ -204,6 +225,14 @@ def __init__(self, xlon, ylat, plev, u_field, v_field, t_field, kmax=49, maxit=1 self._lwa = None self._divergence_eddy_momentum_flux = None + # Temporary solution for GRL computation + self._ua1baro_nhem = None + self._ua2baro_nhem = None + self._ep1baro_nhem = None + self._ep2baro_nhem = None + self._ep3baro_nhem = None + self._ep4_nhem = None + def _compute_prefactor(self): """ Private function. Compute prefactor for normalization by evaluating @@ -211,8 +240,11 @@ def _compute_prefactor(self): using rectangular rule consistent with the integral evaluation in compute_lwa_and_barotropic_fluxes.f90. TODO: evaluate numerical integration scheme used in the fortran module. """ - self.prefactor = sum([math.exp(-k * self.dz / self.scale_height) * self.dz for k in range(1, self.kmax-1)]) - return self.prefactor + self._prefactor = sum([math.exp(-k * self.dz / self.scale_height) * self.dz for k in range(1, self.kmax-1)]) + + @property + def prefactor(self): + return self._prefactor def _check_valid_plev(self, plev, scale_height, kmax, dz): """ @@ -417,6 +449,7 @@ def interpolate_fields(self): self._qgpv_temp, \ self._interpolated_u_temp, \ self._interpolated_v_temp, \ + self._interpolated_avort_temp, \ self._interpolated_theta_temp, \ self._static_stability = \ interpolate_fields( @@ -440,7 +473,7 @@ def interpolate_fields(self): self._interpolated_theta_temp, 0, 2 ) - # Construct a named tuple + # Construct a named tuple # TODO: add absolute vorticity here? But only after testing Interpolated_fields = namedtuple('Interpolated_fields', ['QGPV', 'U', 'V', 'Theta', 'Static_stability']) interpolated_fields = Interpolated_fields( self.qgpv, @@ -626,6 +659,14 @@ def compute_lwa_and_barotropic_fluxes(self, northern_hemisphere_results_only=Fal self._uref_ntemp, self._ptref_ntemp) + # === Access barotropic components of ua1, ua2, ep1, ep2, ep3, ep4: for the use of nhn GLR paper only === + self._ua1baro_nhem = ua1baro_nhem + self._ua2baro_nhem = ua2baro_nhem + self._ep1baro_nhem = ep1baro_nhem + self._ep2baro_nhem = ep2baro_nhem + self._ep3baro_nhem = ep3baro_nhem + self._ep4_nhem = ep4_nhem + # === Compute barotropic flux terms (SHem) === if not northern_hemisphere_results_only: lwa_shem, astarbaro_shem, ua1baro_shem, ubaro_shem, ua2baro_shem,\ @@ -727,6 +768,167 @@ def compute_lwa_and_barotropic_fluxes(self, northern_hemisphere_results_only=Fal self.divergence_eddy_momentum_flux, self.meridional_heat_flux, self.lwa_baro, self.u_baro, self.lwa) return lwa_and_fluxes + # *** Added in Release 0.6.0 *** + # The following internal functions are used to compute results in NHN (2022, GRL): + # - _interpolate_field_dirinv + # - _compute_qref_fawa_and_bc + # - _compute_lwa_flux_dirinv + # They will be refactored in the upcoming releases. + def _interpolate_field_dirinv(self): + """ + Added for NHN 2022 GRL + :return: + """ + self._qgpv_temp, \ + self._interpolated_u_temp, \ + self._interpolated_v_temp, \ + self._interpolated_avort_temp, \ + self._interpolated_theta_temp, \ + self._static_stability_n, \ + self._static_stability_s,\ + self._tn0, self._ts0 = \ + interpolate_fields_direct_inv( + self.kmax, + self.nlat // 2 + self.nlat % 2, + np.swapaxes(self.u_field, 0, 2), + np.swapaxes(self.v_field, 0, 2), + np.swapaxes(self.t_field, 0, 2), + self.plev, + self.planet_radius, + self.omega, + self.dz, + self.scale_height, + self.dry_gas_constant, + self.cp) + + self._check_nan("self._qgpv_temp", self._qgpv_temp) + self._check_nan("self._interpolated_u_temp", self._interpolated_u_temp) + self._check_nan("self._interpolated_v_temp", self._interpolated_v_temp) + self._check_nan("self._interpolated_avort_temp", self._interpolated_avort_temp) + self._check_nan("self._interpolated_theta_temp", self._interpolated_theta_temp) + self._check_nan("self._static_stability_n", self._static_stability_n) + self._check_nan("self._static_stability_s", self._static_stability_s) + self._check_nan("self._tn0", self._tn0) + self._check_nan("self._ts0", self._ts0) + + return self._qgpv_temp, self._interpolated_u_temp, self._interpolated_v_temp, self._interpolated_avort_temp, \ + self._interpolated_theta_temp, self._static_stability_n, self._static_stability_s, self._tn0, self._ts0 + + @staticmethod + def _check_nan(name, var): + nan_num = np.count_nonzero(np.isnan(var)) + if nan_num > 0: + print(f"num of nan in {name}: {np.count_nonzero(np.isnan(var))}.") + + def _compute_qref_fawa_and_bc(self): + """ + Added for NHN 2022 GRL + :return: + """ + # ans = compute_qref_and_fawa_first( + # pv, uu, vort, pt, tn0, ts0, statn, stats, nd, nnd, jb, jd, aa, omega, dz, h, rr, cp) + ans = compute_qref_and_fawa_first( + pv=self._qgpv_temp, + uu=self._interpolated_u_temp, + vort=self._interpolated_avort_temp, + pt=self._interpolated_theta_temp, + tn0=self._tn0, + ts0=self._ts0, + statn=self._static_stability_n, + stats=self._static_stability_s, + nd=self.nlat//2 + self.nlat % 2, # 91 + nnd=self.nlat, # 181 + jb=self.eq_boundary_index, # 5 + jd=self.nlat//2 + self.nlat % 2 - self.eq_boundary_index, # 86 TODO fix its formula + a=self.planet_radius, + omega=self.omega, + dz=self.dz, + h=self.scale_height, + rr=self.dry_gas_constant, + cp=self.cp) + + qref_over_cor, u, ubar, tbar, fawa, ckref, tjk, sjk = ans # unpack tuple + + self._check_nan("qref_over_cor", qref_over_cor) + self._check_nan("u", u) + self._check_nan("ubar", ubar) + self._check_nan("tbar", tbar) + self._check_nan("fawa", fawa) + self._check_nan("ckref", ckref) + self._check_nan("tjk", tjk) + self._check_nan("sjk", sjk) + + for k in range(self.kmax-1, 1, -1): # Fortran indices + ans = matrix_b4_inversion( + k=k, + jmax=self.nlat, + jb=self.eq_boundary_index, # 5 + jd=self.nlat // 2 + self.nlat % 2 - self.eq_boundary_index, # 86 + z=np.arange(0, self.kmax*self.dz, self.dz), + statn=self._static_stability_n, + qref=qref_over_cor, + ckref=ckref, + a=self.planet_radius, + om=self.omega, + dz=self.dz, + h=self.scale_height, + rr=self.dry_gas_constant, + cp=self.cp, + u=u, + sjk=sjk, + tjk=tjk) + qjj, djj, cjj, rj, tj = ans + + # TODO: The inversion algorithm is the bottleneck of the computation + # SciPy is very slow compared to MKL in Fortran... + lu, piv, info = dgetrf(qjj) + qjj, info = dgetri(lu, piv) + + _ = matrix_after_inversion( + k=k, + jb=self.eq_boundary_index, + qjj=qjj, + djj=djj, + cjj=cjj, + tj=tj, + rj=rj, + sjk=sjk, + tjk=tjk) + + tref, qref = upward_sweep( + jmax=self.nlat, + nnd=self.nlat, + jb=self.eq_boundary_index, + sjk=sjk, + tjk=tjk, + ckref=ckref, + tb=self._tn0, + qref_over_cor=qref_over_cor, + u=u, + a=self.planet_radius, + om=self.omega, + dz=self.dz, + h=self.scale_height, + rr=self.dry_gas_constant, + cp=self.cp) + + return qref, u, tref, fawa, ubar, tbar # uref = u + + def _compute_lwa_flux_dirinv(self, qref, uref, tref, fawa, ubar, tbar): + """ + Added for NHN 2022 GRL + :return: + """ + ans = compute_flux_dirinv(pv=self._qgpv_temp, uu=self._interpolated_u_temp, vv=self._interpolated_v_temp, + pt=self._interpolated_theta_temp, tn0=self._tn0, ts0=self._ts0, + statn=self._static_stability_n, stats=self._static_stability_s, + qref=qref, uref=uref, tref=tref, fawa=fawa, ubar=ubar, tbar=tbar, + nnd=self.nlat, jb=self.eq_boundary_index, a=self.planet_radius, om=self.omega, + dz=self.dz, h=self.scale_height, rr=self.dry_gas_constant, cp=self.cp, + prefac=self.prefactor) + # astarbaro, ubaro, urefbaro, ua1baro, ua2baro, ep1baro, ep2baro, ep3baro, ep4, astar1, astar2 = ans + return ans + @property def qgpv(self): """ diff --git a/hn2016_falwa/upward_sweep.f90 b/hn2016_falwa/upward_sweep.f90 new file mode 100644 index 0000000..c53ff78 --- /dev/null +++ b/hn2016_falwa/upward_sweep.f90 @@ -0,0 +1,95 @@ +SUBROUTINE upward_sweep(jmax, kmax, nd, nnd, jb, jd, sjk, tjk, ckref, tb, qref_over_cor, u, tref, qref, a, om, dz, h, rr, cp) + + INTEGER, INTENT(IN) :: jmax, kmax, nd, nnd, jb, jd + REAL, INTENT(IN) :: sjk(jd-2,jd-2,kmax-1),tjk(jd-2,kmax-1),ckref(nd,kmax),tb(kmax),qref_over_cor(nd,kmax) + REAL, INTENT(IN) :: a, om, dz, h, rr, cp + REAL, INTENT(INOUT) :: u(jd,kmax) + REAL, INTENT(OUT) :: qref(nd,kmax), tref(jd,kmax) + real :: tg(kmax) + real :: pjk(jd-2,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 :: sjj(jd-2,jd-2) + real :: pj(jd-2) + + rkappa = rr/cp + pi = acos(-1.) + dp = pi/float(jmax-1) + + + pjk(:,1) = 0. + do k = 1,kmax-1 + pj(:) = pjk(:,k) + sjj(:,:) = sjk(:,:,k) + tj(:) = tjk(:,k) + + do i = 1,jd-2 + yj(i) = 0. + do kk = 1,jd-2 + yj(i) = yj(i)+sjj(i,kk)*pj(kk) + enddo + pjk(i,k+1) = yj(i)+tj(i) + enddo + ! call gemv(sjj,pj,yj) + ! pjk(:,k+1) = yj(:) + tj(:) + enddo + + ! **** Recover u ***** + do k = 1,kmax + do j = 2,jd-1 + u(j,k) = pjk(j-1,k) + enddo + enddo + + ! *** Corner boundary conditions *** + u(1,1) = 0. + u(jd,1) = 0. + ! u(1,kmax) = ubar(1+jb,kmax)*cos(dp*float(jb)) + u(1,kmax) = ckref(1+jb,kmax)/(2.*pi*a)-om*a*cos(dp*float(jb)) + u(jd,kmax) = 0. + + ! *** Divide by cos phi to revover Uref **** + do jj = jb+1,nd-1 + j = jj-jb + phi0 = dp*float(jj-1) + u(j,:) = u(j,:)/cos(phi0) + enddo + u(jd,:) = 2.*u(jd-1,:)-u(jd-2,:) + + ! ******** compute tref ******* + qref(:, :) = qref_over_cor(:, :) ! modify for f2py wrapping purpose + do k = 2,kmax-1 + t00 = 0. + zz = dz*float(k-1) + tref(1,k) = t00 + tref(2,k) = t00 + do j = 2,jd-1 + phi0 = dp*float(j-1+jb) + 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/rr + tref(j+1,k) = tref(j-1,k)+2.*ty*dp + enddo + do j = 1,nd + phi0 = dp*float(j-1) + qref(j,k) = qref_over_cor(j,k)*sin(phi0) + enddo + + tg(k) = 0. + wt = 0. + do jj = jb+1,nd + j = jj-jb + phi0 = dp*float(jj-1) + tg(k) = tg(k)+cos(phi0)*tref(j,k) + wt = wt + cos(phi0) + enddo + tg(k) = tg(k)/wt + tres = tb(k)-tg(k) + tref(:,k) = tref(:,k)+tres + enddo + tref(:,1) = tref(:,2)-tb(2)+tb(1) + tref(:,kmax) = tref(:,kmax-1)-tb(kmax-1)+tb(kmax) + +END SUBROUTINE upward_sweep \ No newline at end of file diff --git a/readme.md b/readme.md index e62c072..cda1055 100644 --- a/readme.md +++ b/readme.md @@ -1,4 +1,4 @@ -## Python Library: hn2016_falwa (v0.5.0) +## Python Library: hn2016_falwa (v0.6.0) [![Build Status](https://github.com/csyhuang/hn2016_falwa/actions/workflows/workflow.yml/badge.svg)](https://github.com/csyhuang/hn2016_falwa/actions/workflows/workflow.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) diff --git a/scripts/nhn_grl2022/graph_plot_module.py b/scripts/nhn_grl2022/graph_plot_module.py new file mode 100644 index 0000000..eb8d643 --- /dev/null +++ b/scripts/nhn_grl2022/graph_plot_module.py @@ -0,0 +1,1087 @@ +""" +This module contains plot functions to reproduce the graphs in NHN GRL2021 +""" +import numpy as np +from cartopy import crs as ccrs +from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter +from matplotlib import pyplot as plt +from netCDF4._netCDF4 import Dataset + + +# *** Shared variables *** +date_stamp = [f'00 UTC {i} June 2021' for i in range(20, 31)] + + +def plot_figure1a(z_filename, u_filename, v_filename): + ncin1 = Dataset(z_filename, 'r', format='NETCDF4') + ncin2 = Dataset(u_filename, 'r', format='NETCDF4') + ncin3 = Dataset(v_filename, 'r', format='NETCDF4') + + zmean = ncin1.variables['z'] + zmean = (np.array(zmean)) + umean = ncin2.variables['u'] + umean = (np.array(umean)) + vmean = ncin3.variables['v'] + vmean = (np.array(vmean)) + + print(zmean.shape) + + zz = np.zeros((44,181,360)) + z = np.zeros((181,360)) + uu = np.zeros((44,181,360)) + u = np.zeros((181,360)) + vv = np.zeros((44,181,360)) + v = np.zeros((181,360)) + m = 76 + while(m < 120): + zz[m-76,:,:] = zmean[m,16,:,:] + uu[m-76,:,:] = umean[m,16,:,:] + vv[m-76,:,:] = vmean[m,16,:,:] + m = m+1 + b = [f'06{i}.png' for i in range(20, 31)] + for n in range(0, 11): + nn = n*4 + for j in range(0, 181): + z[j,:]=zz[nn,180-j,:]/9.81 + u[j,:]=uu[nn,180-j,:] + v[j,:]=vv[nn,180-j,:] + j = j+1 + cl1 = np.arange(9600,11300,100) + cl2 = np.arange(0,95,5) + plt.rcParams.update({'font.size': 15}) + x = np.arange(0,360) + y = np.arange(0,181)-90. + plt.rcParams.update({'font.size':15, 'text.usetex': False}) + fig = plt.figure(figsize=(8, 4)) + ax5 = fig.add_subplot(1,1,1,projection=ccrs.PlateCarree(180)) + plt.xlim(140, 280) + plt.ylim(10, 80) + # plot.title(''+day[n]) + if(n > 9): + plt.xlabel('Longitude', fontsize = 22) + plt.ylabel('Latitude', fontsize = 22) + ax5.set_extent([-220, -80, 10, 80], ccrs.PlateCarree()) + ax5.coastlines(color='white',alpha = 0.7) + ax5.set_aspect('auto', adjustable=None) + if(n > 9): + ax5.set_xticks([140,160,180,200,220,240,260,280], crs=ccrs.PlateCarree()) + ax5.set_yticks([10, 20, 30, 40, 50, 60, 70, 80], crs=ccrs.PlateCarree()) + lon_formatter = LongitudeFormatter(zero_direction_label=True) + lat_formatter = LatitudeFormatter() + ax5.xaxis.set_major_formatter(lon_formatter) + ax5.yaxis.set_major_formatter(lat_formatter) + ott = ax5.contourf(x,y,np.sqrt(u*u+v*v),levels=cl2,transform=ccrs.PlateCarree(),cmap='rainbow') + fig.colorbar(ott,ax=ax5,label='wind speed (m/s)') + ott = ax5.contour(x,y,z,levels=cl1,colors='black',transform=ccrs.PlateCarree(),linewidths=1) + ax5.clabel(ott, ott.levels,fmt='%5i') + plt.savefig(b[n], bbox_inches='tight', dpi =600) + plt.close() + + +def plot_figure1b(t_filename): + ncin1 = Dataset(t_filename, 'r', format='NETCDF4') + + tmean = ncin1.variables['t'] + tmean = (np.array(tmean)) + + print(tmean.shape) + + tt = np.zeros((44,181,360)) + t = np.zeros((181,360)) + + r = 287. + cp = 1004. + kappa = r/cp + + for m in range(76, 120): + tt[m-76,:,:] = tmean[m,20,:,:]*np.power((1000./450.),kappa) + + b = [f'T_06{i}.png' for i in range(20, 31)] + + for n in range(0, 11): + nn = n*4 + for j in range(0, 181): + t[j,:]=tt[nn,180-j,:] + cl1 = np.arange(290,342,2) # 450 hPa + cl2 = np.arange(0,95,5) + x = np.arange(0,360) + y = np.arange(0,181)-90. + plt.rcParams.update({'font.size': 15}) + fig = plt.figure(figsize=(8, 4)) + ax5 = fig.add_subplot(1,1,1,projection=ccrs.PlateCarree(180)) + plt.xlim(140, 280) + plt.ylim(10, 80) + if(n > 9): + plt.xlabel('Longitude', fontsize = 22) + ax5.set_extent([-220, -80, 10, 80], ccrs.PlateCarree()) + ax5.coastlines(color='black',alpha = 0.7) + ax5.set_aspect('auto', adjustable=None) + if(n > 9): + ax5.set_xticks([140,160,180,200,220,240,260,280], crs=ccrs.PlateCarree()) + lon_formatter = LongitudeFormatter(zero_direction_label=True) + lat_formatter = LatitudeFormatter() + ax5.xaxis.set_major_formatter(lon_formatter) + ax5.yaxis.set_major_formatter(lat_formatter) + ott = ax5.contourf(x,y,t,levels=cl1,transform=ccrs.PlateCarree(),cmap='rainbow') + fig.colorbar(ott,ax=ax5,label='Kelvin') + plt.savefig(b[n], bbox_inches='tight', dpi =600) + plt.close() + + +def plot_figure1c(t2m_filename): + ncin1 = Dataset(t2m_filename, 'r', format='NETCDF4') + tm = ncin1.variables['t2m'] + tm = (np.array(tm)) + + print(tm.shape) + + tt = np.zeros((44,181,360)) + t = np.zeros((181,360)) + + r = 287. + cp = 1004. + kappa = r/cp + + for m in range(76, 120): + tt[m-76,:,:] = tm[m,:,:] + + b = [f'2T_06{i}.png' for i in range(20, 31)] + n = 0 + for n in range(0, 11): + nn = n*4 + + for j in range(0, 181): + t[j,:]=tt[nn,180-j,:] + + cl1 = np.arange(250,325,5) + cl2 = np.arange(0,95,5) + plt.rcParams.update({'font.size': 15}) + x = np.arange(0,360) + y = np.arange(0,181)-90. + fig = plt.figure(figsize=(8, 4)) + ax5 = fig.add_subplot(1,1,1,projection=ccrs.PlateCarree(180)) + plt.xlim(140, 280) + plt.ylim(10, 80) + if(n > 9): + plt.xlabel('Longitude', fontsize=22) + ax5.set_extent([-220, -80, 10, 80], ccrs.PlateCarree()) + ax5.coastlines(color='black',alpha = 0.7) + ax5.set_aspect('auto', adjustable=None) + if(n > 9): + ax5.set_xticks([140,160,180,200,220,240,260,280], crs=ccrs.PlateCarree()) + lon_formatter = LongitudeFormatter(zero_direction_label=True) + lat_formatter = LatitudeFormatter() + ax5.xaxis.set_major_formatter(lon_formatter) + ax5.yaxis.set_major_formatter(lat_formatter) + ott = ax5.contourf(x,y,t,levels=cl1,transform=ccrs.PlateCarree(),cmap='rainbow') + fig.colorbar(ott,ax=ax5,label='Kelvin') + plt.savefig(b[n], bbox_inches='tight', dpi =600) + plt.close() + + +def plot_figure1d_2a(t_filename): + ncin1 = Dataset(t_filename, 'r', format='NETCDF4') + + tmean = ncin1.variables['t'] + tmean = (np.array(tmean)) + + ttheta = np.zeros((44, 37, 360)) + theta = np.zeros((37, 360)) + z = np.zeros(37) + p = np.array([1., 2., 3., 5., 7., 10., 20., 30., 50., + 70., 100., 125., 150., 175., 200., 225., 250., 300., + 350., 400., 450., 500., 550., 600., 650., 700., 750., + 775., 800., 825., 850., 875., 900., 925., 950., 975., 1000.]) + r = 287. + cp = 1004. + kappa = r / cp + for m in range(76, 120): + for k in range(37): + z[36 - k] = -8000. * np.log(p[k] / 1000.) # pseudoheight + ttheta[m - 76, 36 - k, :] = tmean[m, k, 41, :] + ttheta[m - 76, 36 - k, :] = ttheta[m - 76, 36 - k, :] * np.power((1000. / p[k]), kappa) # potential temp + + b = [f'THETA_06{i}.png' for i in range(20, 31)] + for n in range(11): + nn = n * 4 + theta[:, :] = ttheta[nn, :, :] + cl1 = np.arange(250, 365, 5) + x = np.arange(0, 360) + plt.rcParams.update({'font.size': 15, 'text.usetex': False}) + fig = plt.figure(figsize=(8, 4)) + ax5 = fig.add_subplot(1, 1, 1, projection=ccrs.PlateCarree(180)) + plt.xlim(140, 280) + plt.ylim(0, 10) + # plot.title('49$^\circ$N '+day[n]) + if (n > 9): + plt.xlabel('Longitude', fontsize=22) + plt.ylabel('pseudoheight (km)', fontsize=22) + ax5.set_extent([-220, -80, 0, 10], ccrs.PlateCarree()) + ax5.set_aspect('auto', adjustable=None) + if (n > 9): + ax5.set_xticks([140, 160, 180, 200, 220, 240, 260, 280], crs=ccrs.PlateCarree()) + ax5.set_yticks([0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10], crs=ccrs.PlateCarree()) + lon_formatter = LongitudeFormatter(zero_direction_label=True) + lat_formatter = LatitudeFormatter() + ax5.xaxis.set_major_formatter(lon_formatter) + ott = ax5.contourf(x, z / 1000., theta, levels=cl1, transform=ccrs.PlateCarree(), cmap='rainbow') + fig.colorbar(ott, ax=ax5, label='Kelvin') + ott = ax5.contour(x, z / 1000., theta, levels=cl1, transform=ccrs.PlateCarree(), colors='black', linewidths=0.5) + ott = ax5.contour(x, z / 1000., theta, levels=np.arange(320, 325, 5), transform=ccrs.PlateCarree(), + colors='black', linewidths=1) + ax5.clabel(ott, ott.levels, fmt='%5i') + plt.savefig(b[n], bbox_inches='tight', dpi=600) + plt.close() + + plt.rcParams.update({'font.size': 16}) + fig = plt.figure(figsize=(6, 4)) + ax5.set_yticks([0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10]) + plt.title('49$^\circ$N 119$^\circ$W 00 UTC') + plt.xlabel('Kelvin') + plt.ylabel('pseudoheight (km)') + plt.xlim(290, 360) + plt.ylim(0, 10) + + lcolor = np.array(['blue', 'red', 'green', 'black']) + lstyle = np.array(['dotted', 'dashed', 'dashdot', 'solid']) + lalpha = np.array([1, 1, 1, 1]) + + for n in range(2, 6): + thetaz = np.zeros(37) + nn = n * 8 + thetaz[:] = ttheta[nn, :, 241] + + for i in range(37): + if (z[i] < 1000.): + thetaz[i] = np.nan + + fig = plt.plot(thetaz, z / 1000., color=lcolor[n - 2], linestyle=lstyle[n - 2], alpha=lalpha[n - 2]) + + plt.savefig('t_profile.png', bbox_inches='tight', dpi=600) + plt.close() + + +def plot_figure3_and_S1(lwa_flux_filename): + ncin1 = Dataset(lwa_flux_filename, 'r', format='NETCDF4') + + lwa = ncin1.variables['lwa'] + lwa = (np.array(lwa)) + + z = np.zeros((91,360)) + z[:,:] = lwa[100,:,:]-lwa[76,:,:] # m = 100 is 00 UTC 26 June 2021, m = 76 is 00 UTC 20 June 2021 + + zs = np.zeros((91,360)) # smoothed z # + + #### smoothing in longitude #### + n = 5 # smoothing width # + j = 0 + while(j < 91): + zx = np.zeros(360) + zx[:] = z[j,:] + nn = -n + while(nn < n+1): + zy = np.roll(zx,nn) + zs[j,:] = zs[j,:] + zy[:]/(2*n+1) + nn = nn+1 + j = j+1 + + + cl2 = np.arange(-80,90,10) + x = np.arange(0,360) + y = np.arange(0,91) + plt.rcParams.update({'font.size': 16}) + fig = plt.figure(figsize=(10, 5)) + ax5 = fig.add_subplot(1,1,1,projection=ccrs.PlateCarree(180)) + plt.xlim(140, 280) + plt.ylim(10, 80) + plt.title('Column LWA Change 00 UTC 20 - 26 June 2021') + #plot.xlabel('Longitude') + plt.ylabel('Latitude', fontsize=22) + ax5.set_extent([-220, -80, 10, 80], ccrs.PlateCarree()) + ax5.coastlines(alpha = 0.7) + ax5.set_aspect('auto', adjustable=None) + ax5.set_xticks([140,160,180,200,220,240,260,280], crs=ccrs.PlateCarree()) + ax5.set_yticks([10, 20, 30, 40, 50, 60, 70, 80], crs=ccrs.PlateCarree()) + lon_formatter = LongitudeFormatter(zero_direction_label=True) + lat_formatter = LatitudeFormatter() + ax5.xaxis.set_major_formatter(lon_formatter) + ax5.yaxis.set_major_formatter(lat_formatter) + ott = ax5.contourf(x,y,zs,levels=cl2,transform=ccrs.PlateCarree(),cmap='rainbow') + fig.colorbar(ott,ax=ax5,label='LWA (m/s)') + plt.savefig('dLWA_0.png', bbox_inches='tight', dpi =600) + + ncin1 = Dataset(lwa_flux_filename, 'r', format='NETCDF4') + ua1 = ncin1.variables['ua1'] + ua1 = (np.array(ua1)) + ncin2 = Dataset(lwa_flux_filename, 'r', format='NETCDF4') + ua2 = ncin2.variables['ua2'] + ua2 = (np.array(ua2)) + ncin3 = Dataset(lwa_flux_filename, 'r', format='NETCDF4') + ep1 = ncin3.variables['ep1'] + ep1 = (np.array(ep1)) + ncin4 = Dataset(lwa_flux_filename, 'r', format='NETCDF4') + ep2 = ncin4.variables['ep2'] + ep2 = (np.array(ep2)) + ncin5 = Dataset(lwa_flux_filename, 'r', format='NETCDF4') + ep3 = ncin5.variables['ep3'] + ep3 = (np.array(ep3)) + ncin6 = Dataset(lwa_flux_filename, 'r', format='NETCDF4') + ep4 = ncin6.variables['ep4'] + ep4 = (np.array(ep4)) + + f1 = np.zeros((91,360)) + f2 = np.zeros((91,360)) + f11 = np.zeros((91,360)) + f22 = np.zeros((91,360)) + + z1 = np.zeros((91,360)) + z2 = np.zeros((91,360)) + z3 = np.zeros((91,360)) + dt = 3600.*6. + a = 6378000. + dl = 2.*np.pi/360. + dp = 2.*np.pi/360. + m = 76 # m = 76 is 20 June 2021 00 UTC + while(m < 100): # m = 100 is 26 June 2021 00 UTC + #m = 52 # m = 52 is 14 June 2021 00 UTC + #while(m < 76): # m = 76 is 20 June 2021 00 UTC + + z3[:,:] = z3[:,:]+0.5*dt*ep4[m,:,:] + z3[:,:] = z3[:,:]+0.5*dt*ep4[m+1,:,:] + f1[:,:] = f1[:,:]+(0.5/24.)*(ua1[m,:,:]+ua2[m,:,:]+ep1[m,:,:]) + f1[:,:] = f1[:,:]+(0.5/24.)*(ua1[m+1,:,:]+ua2[m+1,:,:]+ep1[m+1,:,:]) + f11[:,:] = f11[:,:]+(0.5/24.)*(ua1[m-24,:,:]+ua2[m-24,:,:]+ep1[m-24,:,:]) + f11[:,:] = f11[:,:]+(0.5/24.)*(ua1[m-23,:,:]+ua2[m-23,:,:]+ep1[m-23,:,:]) + j = 0 + while(j < 90): + phi = dp*j + const = 0.5*dt/(2.*a*np.cos(phi)*dl) + z2[j,:]=z2[j,:]+const*(ep2[m,j,:]-ep3[m,j,:]) + z2[j,:]=z2[j,:]+const*(ep2[m+1,j,:]-ep3[m+1,j,:]) + f2[j,:] = f2[j,:]+(0.25/24.)*(ep2[m,j,:]+ep3[m,j,:])/np.cos(phi) + f2[j,:] = f2[j,:]+(0.25/24.)*(ep2[m+1,j,:]+ep3[m+1,j,:])/np.cos(phi) + f22[j,:] = f22[j,:]+(0.25/24.)*(ep2[m-24,j,:]+ep3[m-24,j,:])/np.cos(phi) + f22[j,:] = f22[j,:]+(0.25/24.)*(ep2[m-23,j,:]+ep3[m-23,j,:])/np.cos(phi) + i = 1 + while(i < 359): + z1[j,i] = z1[j,i]-const*(ua1[m,j,i+1]+ua2[m,j,i+1]+ep1[m,j,i+1]-ua1[m,j,i-1]-ua2[m,j,i-1]-ep1[m,j,i-1]) + z1[j,i] = z1[j,i]-const*(ua1[m+1,j,i+1]+ua2[m+1,j,i+1]+ep1[m+1,j,i+1]-ua1[m+1,j,i-1]-ua2[m+1,j,i-1]-ep1[m+1,j,i-1]) + i = i+1 + z1[j,0] = z1[j,0]-const*(ua1[m,j,1]+ua2[m,j,1]+ep1[m,j,1]-ua1[m,j,359]-ua2[m,j,359]-ep1[m,j,359]) + z1[j,0] = z1[j,0]-const*(ua1[m+1,j,1]+ua2[m+1,j,1]+ep1[m+1,j,1]-ua1[m+1,j,359]-ua2[m+1,j,359]-ep1[m+1,j,359]) + z1[j,359] = z1[j,359]-const*(ua1[m,j,0]+ua2[m,j,0]+ep1[m,j,0]-ua1[m,j,358]-ua2[m,j,358]-ep1[m,j,358]) + z1[j,359] = z1[j,359]-const*(ua1[m+1,j,0]+ua2[m+1,j,0]+ep1[m+1,j,0]-ua1[m+1,j,358]-ua2[m+1,j,358]-ep1[m+1,j,358]) + j = j+1 + m = m+1 + + z1s = np.zeros((91,360)) # smoothed z1 # + z2s = np.zeros((91,360)) # smoothed z2 # + z3s = np.zeros((91,360)) # smoothed z3 # + + #### smoothing in longitude #### + j = 0 + while(j < 91): + z1x = np.zeros(360) + z1x[:] = z1[j,:] + z2x = np.zeros(360) + z2x[:] = z2[j,:] + z3x = np.zeros(360) + z3x[:] = z3[j,:] + nn = -n + while(nn < n+1): + z1y = np.roll(z1x,nn) + z1s[j,:] = z1s[j,:] + z1y[:]/(2*n+1) + z2y = np.roll(z2x,nn) + z2s[j,:] = z2s[j,:] + z2y[:]/(2*n+1) + z3y = np.roll(z3x,nn) + z3s[j,:] = z3s[j,:] + z3y[:]/(2*n+1) + nn = nn+1 + j = j+1 + + ##### Wind vectors ###### + + x1 = np.arange(0,24)*15.+5. + y1 = np.arange(0,30)*3. + xx,yy = np.meshgrid(x1,y1) + uu = np.zeros((30,24)) + vv = np.zeros((30,24)) + + j = 0 + while(j < 30): + i = 0 + while(i < 24): + uu[j,i] = f1[j*3,i*15+5]-f11[j*3,i*15+5] + vv[j,i] = f2[j*3,i*15+5]-f22[j*3,i*15+5] + i = i+1 + j = j+1 + + + cl1 = np.arange(-200,220,20) + x = np.arange(0,360) + y = np.arange(0,91) + plt.rcParams.update({'font.size': 16}) + fig = plt.figure(figsize=(10, 5)) + ax6 = fig.add_subplot(1,1,1,projection=ccrs.PlateCarree(180)) + plt.xlim(0, 360) + plt.ylim(10, 80) + plt.title('Integrated column -div (Fx + Fy) 20 - 26 June 2021') + plt.ylabel('Latitude', fontsize=22) + ax6.set_extent([-220, -80, 10, 80], ccrs.PlateCarree()) + ax6.coastlines(alpha = 0.7) + ax6.set_aspect('auto', adjustable=None) + ax6.set_xticks([140,160,180,200,220,240,260,280], crs=ccrs.PlateCarree()) + ax6.set_yticks([10, 20, 30, 40, 50, 60, 70, 80], crs=ccrs.PlateCarree()) + lon_formatter = LongitudeFormatter(zero_direction_label=True) + lat_formatter = LatitudeFormatter() + ax6.xaxis.set_major_formatter(lon_formatter) + ax6.yaxis.set_major_formatter(lat_formatter) + ott = ax6.contourf(x,y,z1s+z2s,levels=cl1,transform=ccrs.PlateCarree(),cmap='rainbow') + fig.colorbar(ott,ax=ax6,label='(m/s)') + ax6.quiver(xx[2:-2,:],yy[2:-2,:],uu[2:-2, :],vv[2:-2, :],transform=ccrs.PlateCarree()) + plt.savefig('divFx+Fy_0.png', bbox_inches='tight', dpi =600) + + cl1 = np.arange(-200,220,20) + x = np.arange(0,360) + y = np.arange(0,91) + plt.rcParams.update({'font.size': 16}) + fig = plt.figure(figsize=(10, 5)) + ax6 = fig.add_subplot(1,1,1,projection=ccrs.PlateCarree(180)) + plt.xlim(0, 360) + plt.ylim(10, 80) + plt.title('Integrated column -div Fy 20 - 26 June 2021') + plt.xlabel('Longitude', fontsize=22) + plt.ylabel('Latitude', fontsize=22) + ax6.set_extent([-220, -80, 10, 80], ccrs.PlateCarree()) + ax6.coastlines(alpha = 0.7) + ax6.set_aspect('auto', adjustable=None) + ax6.set_xticks([140,160,180,200,220,240,260,280], crs=ccrs.PlateCarree()) + ax6.set_yticks([10, 20, 30, 40, 50, 60, 70, 80], crs=ccrs.PlateCarree()) + lon_formatter = LongitudeFormatter(zero_direction_label=True) + lat_formatter = LatitudeFormatter() + ax6.xaxis.set_major_formatter(lon_formatter) + ax6.yaxis.set_major_formatter(lat_formatter) + ott = ax6.contourf(x,y,z2s,levels=cl1,transform=ccrs.PlateCarree(),cmap='rainbow') + fig.colorbar(ott,ax=ax6,label='(m/s)') + ax6.quiver(xx[2:-2,:],yy[2:-2,:],uu[2:-2, :],vv[2:-2, :],transform=ccrs.PlateCarree()) + plt.savefig('divFy_0.png', bbox_inches='tight', dpi =600) + + cl1 = np.arange(-200,220,20) + x = np.arange(0,360) + y = np.arange(0,91) + plt.rcParams.update({'font.size': 16}) + fig = plt.figure(figsize=(10, 5)) + ax6 = fig.add_subplot(1,1,1,projection=ccrs.PlateCarree(180)) + plt.xlim(0, 360) + plt.ylim(10, 80) + plt.title('Integrated low-level source 20 - 26 June 2021') + ax6.set_extent([-220, -80, 10, 80], ccrs.PlateCarree()) + ax6.coastlines(alpha = 0.7) + ax6.set_aspect('auto', adjustable=None) + ax6.set_xticks([140,160,180,200,220,240,260,280], crs=ccrs.PlateCarree()) + ax6.set_yticks([10, 20, 30, 40, 50, 60, 70, 80], crs=ccrs.PlateCarree()) + lon_formatter = LongitudeFormatter(zero_direction_label=True) + lat_formatter = LatitudeFormatter() + ax6.xaxis.set_major_formatter(lon_formatter) + ax6.yaxis.set_major_formatter(lat_formatter) + ott = ax6.contourf(x,y,z3s,levels=cl1,transform=ccrs.PlateCarree(),cmap='rainbow') + fig.colorbar(ott,ax=ax6,label='(m/s)') + plt.savefig('EP4_0.png', bbox_inches='tight', dpi =600) + + cl1 = np.arange(-200,220,20) + x = np.arange(0,360) + y = np.arange(0,91) + fig = plt.figure(figsize=(10, 5)) + ax6 = fig.add_subplot(1,1,1,projection=ccrs.PlateCarree(180)) + plt.xlim(0, 360) + plt.ylim(10, 80) + plt.title('Integrated residual 20 - 26 June 2021') + ax6.set_extent([-220, -80, 10, 80], ccrs.PlateCarree()) + ax6.coastlines(alpha = 0.7) + ax6.set_aspect('auto', adjustable=None) + ax6.set_xticks([140,160,180,200,220,240,260,280], crs=ccrs.PlateCarree()) + ax6.set_yticks([10, 20, 30, 40, 50, 60, 70, 80], crs=ccrs.PlateCarree()) + lon_formatter = LongitudeFormatter(zero_direction_label=True) + lat_formatter = LatitudeFormatter() + ax6.xaxis.set_major_formatter(lon_formatter) + ax6.yaxis.set_major_formatter(lat_formatter) + ott = ax6.contourf(x,y,zs-z1s-z2s-z3s,levels=cl1,transform=ccrs.PlateCarree(),cmap='rainbow') + fig.colorbar(ott,ax=ax6,label='(m/s)') + ax6.quiver(xx[2:-2,:],yy[2:-2,:],uu[2:-2, :],vv[2:-2, :],transform=ccrs.PlateCarree()) + plt.savefig('Residual_0.png', bbox_inches='tight', dpi =600) + plt.close("all") + + +def plot_figure3e(mtnlwrf_filename, mtnlwrfcs_filename): + """ + :param mtnlwrf_filename: netCDF fileof net OLR data + :param mtnlwrfcs_filename: netCDF fileof clear sky OLR data + :return: + """ + ncin1 = Dataset(mtnlwrf_filename, 'r', format='NETCDF4') + ncin2 = Dataset(mtnlwrfcs_filename, 'r', format='NETCDF4') + olr = ncin1.variables['mtnlwrf'] + olr = (np.array(olr)) + olrcs = ncin2.variables['mtnlwrfcs'] + olrcs = (np.array(olrcs)) + + tt = np.zeros((44,181,360)) + t = np.zeros((181,360)) + m = 77 + for m in range(77,120): + tt[m-77,:,:] = olr[m,:,:] + + b = [f'OLR_06{i}.png' for i in range(20, 31)] + + for n in range(0, 11): + nn = n*4 + + for j in range(0, 181): + t[j,:]=tt[nn,180-j,:] + + cl1 = np.arange(80,390,10) + x = np.arange(0,360) + y = np.arange(0,181)-90. + plt.rcParams.update({'font.size':16, 'text.usetex': False}) + fig = plt.figure(figsize=(10, 5)) + ax5 = fig.add_subplot(1,1,1,projection=ccrs.PlateCarree(180)) + plt.xlim(140, 280) + plt.ylim(10, 80) + plt.title('OLR ' + date_stamp[n]) + plt.xlabel('Longitude', fontsize=22) + plt.ylabel('Latitude', fontsize=22) + ax5.set_extent([-220, -80, 10, 80], ccrs.PlateCarree()) + ax5.coastlines(color='black',alpha = 0.7) + ax5.set_aspect('auto', adjustable=None) + ax5.set_xticks([140,160,180,200,220,240,260,280], crs=ccrs.PlateCarree()) + ax5.set_yticks([10, 20, 30, 40, 50, 60, 70, 80], crs=ccrs.PlateCarree()) + lon_formatter = LongitudeFormatter(zero_direction_label=True) + lat_formatter = LatitudeFormatter() + ax5.xaxis.set_major_formatter(lon_formatter) + ax5.yaxis.set_major_formatter(lat_formatter) + ott = ax5.contourf(x,y,-t,levels=cl1,transform=ccrs.PlateCarree(),cmap='rainbow') + fig.colorbar(ott,ax=ax5,label='W/m$^2$') + plt.savefig(b[n], bbox_inches='tight', dpi =600) + plt.close() + + +def plot_figure3f(tcw_filename, tcwv_filename, sp_filename): + """ + :param tcw_filename: filename of netCDF file with total column water (kg/m^2) + :param tcwv_filename: filename of netCDF file with total column water vapor (kg/m^2) + :param sp_filename: filename of netCDF file with sea level pressure (hPa) + :return: + """ + ncin1 = Dataset(tcw_filename, 'r', format='NETCDF4') + ncin2 = Dataset(tcwv_filename, 'r', format='NETCDF4') + ncin3 = Dataset(sp_filename, 'r', format='NETCDF4') + cw = ncin1.variables['tcw'] + cw = (np.array(cw)) + cwv = ncin2.variables['tcwv'] + cwv = (np.array(cwv)) + sp = ncin3.variables['sp'] + sp = (np.array(sp)) + + tt = np.zeros((44,181,360)) + pp = np.zeros((44,181,360)) + t = np.zeros((181,360)) + p = np.zeros((181,360)) + + for m in range(77, 120): + tt[m-77,:,:] = cw[m,:,:]-cwv[m,:,:] + pp[m-77,:,:] = sp[m,:,:]/100. + + b = [f'CW_06{i}.png' for i in range(20, 31)] + + for n in range(0, 11): + nn = n*4 + + for j in range(0, 181): + t[j,:]=tt[nn,180-j,:] + p[j,:]=pp[nn,180-j,:] + + cl1 = np.arange(-0.1,3.6,0.1) + c12 = np.arange(980,1032,4) + x = np.arange(0,360) + y = np.arange(0,181)-90. + plt.rcParams.update({'font.size':16, 'text.usetex': False}) + fig = plt.figure(figsize=(10, 5)) + ax5 = fig.add_subplot(1,1,1,projection=ccrs.PlateCarree(180)) + plt.xlim(140, 280) + plt.ylim(10, 80) + plt.title('Column water ' + date_stamp[n]) + plt.xlabel('Longitude', fontsize=22) + #plot.ylabel('Latitude') + ax5.set_extent([-220, -80, 10, 80], ccrs.PlateCarree()) + ax5.coastlines(color='white',alpha = 0.7) + ax5.set_aspect('auto', adjustable=None) + ax5.set_xticks([140,160,180,200,220,240,260,280], crs=ccrs.PlateCarree()) + ax5.set_yticks([10, 20, 30, 40, 50, 60, 70, 80], crs=ccrs.PlateCarree()) + lon_formatter = LongitudeFormatter(zero_direction_label=True) + lat_formatter = LatitudeFormatter() + ax5.xaxis.set_major_formatter(lon_formatter) + ax5.yaxis.set_major_formatter(lat_formatter) + ott = ax5.contourf(x,y,t,levels=cl1,transform=ccrs.PlateCarree(),cmap='rainbow') + fig.colorbar(ott,ax=ax5,label='kg/m$^2$') + ott = ax5.contour(x,y,p,levels=c12,colors='white',alpha=0.5,transform=ccrs.PlateCarree(),linewidths=1) + ax5.clabel(ott, ott.levels,fmt='%5i') + plt.savefig(b[n], bbox_inches='tight', dpi =600) + plt.close() + + +def plot_figure4(lwa_flux_filename): + dt = 3600. * 6. + a = 6378000. + dl = 2. * np.pi / 360. + dp = 2. * np.pi / 360. + + jm = 49 + jp = 49 + + ncin1 = Dataset(lwa_flux_filename, 'r', format='NETCDF4') + + lwa = ncin1.variables['lwa'] + lwa = (np.array(lwa)) + + z = np.zeros((124,360)) # wave activity tendency + w = np.zeros((124,360)) # wave activity + + for m in range(1, 120): + + cos0 = 0 + for j in range(jm, jp+1): + + phi = dp*j + z[m,:] = z[m,:]+np.cos(phi)*(lwa[m+1,j,:]-lwa[m-1,j,:])/(2.*dt) # wave activity tendency at 49 N + w[m,:] = w[m,:]+np.cos(phi)*lwa[m,j,:] + cos0 = cos0 + np.cos(phi) + + z[m,:] = z[m,:]/cos0 + w[m,:] = w[m,:]/cos0 + + ncin1 = Dataset(lwa_flux_filename, 'r', format='NETCDF4') + ua1 = ncin1.variables['ua1'] + ua2 = ncin1.variables['ua2'] + ep1 = ncin1.variables['ep1'] + ep2 = ncin1.variables['ep2'] + ep3 = ncin1.variables['ep3'] + ep4 = (np.array(ncin1.variables['ep4'])) + + f1 = np.zeros((124,360)) + + z1 = np.zeros((124,360)) + z2 = np.zeros((124,360)) + z3 = np.zeros((124,360)) + z4 = np.zeros((124,360)) + + + for i in range(1, 359): + + for j in range(jm, jp+1): + phi = dp*j + const = 1./(2.*a*dl) + z1[:,i] = z1[:,i]-const*(ua1[:,j,i+1]+ua2[:,j,i+1]+ep1[:,j,i+1]-ua1[:,j,i-1]-ua2[:,j,i-1]-ep1[:,j,i-1]) + + z1[:,i] = z1[:,i]/cos0 + + + for i in range(0, 360): + + for j in range(jm, jp+1): + phi = dp*j + const = 1./(2.*a*dl) + f1[:,i]=f1[:,i]+(ua1[:,j,i]+ua2[:,j,i]+ep1[:,j,i])*np.cos(phi) # zonal wave activity flux + z2[:,i]=z2[:,i]+const*(ep2[:,j,i]-ep3[:,j,i]) # metisional convergence of wave activity flux + z3[:,i]=z3[:,i]+np.cos(phi)*ep4[:,j,i] # bottom source + + f1[:,i] = f1[:,i]/cos0 + z2[:,i] = z2[:,i]/cos0 + z3[:,i] = z3[:,i]/cos0 + + + zs = np.zeros((124,360)) # smoothed z # + ws = np.zeros((124,360)) # smoothed w # + f1s = np.zeros((124,360)) # smoothed f1 # + z1s = np.zeros((124,360)) # smoothed z1 # + z2s = np.zeros((124,360)) # smoothed z2 # + z3s = np.zeros((124,360)) # smoothed z3 # + + #### smoothing in longitude #### + n = 4 # smoothing width = 2n+1 + + for j in range(0, 124): + zx = np.zeros(360) + zx[:] = z[j,:] + wx = np.zeros(360) + wx[:] = w[j,:] + f1x = np.zeros(360) + f1x[:] = f1[j,:] + z1x = np.zeros(360) + z1x[:] = z1[j,:] + z2x = np.zeros(360) + z2x[:] = z2[j,:] + z3x = np.zeros(360) + z3x[:] = z3[j,:] + + for nn in range(-n, n+1): + + wy = np.roll(wx,nn) + ws[j,:] = ws[j,:] + wy[:]/(2*n+1) + f1y = np.roll(f1x,nn) + f1s[j,:] = f1s[j,:] + f1y[:]/(2*n+1) + zy = np.roll(zx,nn) + zs[j,:] = zs[j,:] + zy[:]/(2*n+1) + z1y = np.roll(z1x,nn) + z1s[j,:] = z1s[j,:] + z1y[:]/(2*n+1) + z2y = np.roll(z2x,nn) + z2s[j,:] = z2s[j,:] + z2y[:]/(2*n+1) + z3y = np.roll(z3x,nn) + z3s[j,:] = z3s[j,:] + z3y[:]/(2*n+1) + + z4[:,:]=zs[:,:]-z1s[:,:]-z2s[:,:]-z3s[:,:] # residual + + ############################################## + + cl2 = np.arange(0,135,5) + x = np.arange(0,360) + y = np.arange(0,124)*0.25 + plt.rcParams.update({'font.size': 28}) + fig,ax5 = plt.subplots(1, figsize=(8, 8)) + plt.xlim(140, 280) + plt.ylim(20, 29.75) + plt.title('Column LWA 49$^\circ$N') + plt.xlabel('Longitude') + plt.ylabel('Day') + ott = ax5.contourf(x,y,w,levels=cl2,cmap='rainbow') + fig.colorbar(ott,ax=ax5,label='(ms$^{-1}$)') + plt.savefig('HovmollerLWA.png', bbox_inches='tight', dpi =600) + + ############################################## + + cl2 = np.arange(-100,1500,100) + x = np.arange(0,360) + y = np.arange(0,124)*0.25 + plt.rcParams.update({'font.size': 28}) + fig,ax5 = plt.subplots(1, figsize=(8, 8)) + plt.xlim(140, 280) + plt.ylim(20, 29.75) + plt.title(' 49$^\circ$N') + plt.xlabel('Longitude') + plt.ylabel('Day') + ott = ax5.contourf(x,y,f1,levels=cl2,cmap='rainbow') + fig.colorbar(ott,ax=ax5,label='(m$^2$s$^{-2}$)') + plt.savefig('HovmollerFx.png', bbox_inches='tight', dpi =600) + + + ############################################## + + cl2 = np.arange(-10.,9,1) + x = np.arange(0,360) + y = np.arange(0,124)*0.25 + plt.rcParams.update({'font.size': 28}) + fig,ax5 = plt.subplots(1, figsize=(8, 8)) + plt.xlim(140, 280) + plt.ylim(20, 29.5) + plt.title('Term (IV) 49$^\circ$N') + plt.xlabel('Longitude') + plt.ylabel('Day') + ott = ax5.contourf(x,y,z4*10000,levels=cl2,cmap='rainbow') + fig.colorbar(ott,ax=ax5,label='(10$^{-4}$ ms$^{-2}$)') + plt.savefig('HovmollerRes.png', bbox_inches='tight', dpi =600) + #plot.show() + + + ############################################## + + cl2 = np.arange(-10.,9,1) + x = np.arange(0,360) + y = np.arange(0,124)*0.25 + plt.rcParams.update({'font.size': 28}) + fig,ax5 = plt.subplots(1, figsize=(8, 8)) + plt.xlim(140, 280) + plt.ylim(20, 29.75) + plt.title('Terms (I)+(II) 49$^\circ$N') + plt.xlabel('Longitude') + plt.ylabel('Day') + ott = ax5.contourf(x,y,(z1s+z2s+z3s)*10000,levels=cl2,cmap='rainbow') + fig.colorbar(ott,ax=ax5,label='(10$^{-4}$ ms$^{-2}$)') + plt.savefig('HovmollerFxy.png', bbox_inches='tight', dpi =600) + + +def plot_figure5(lwa_flux_filename): + + ncin1 = Dataset(lwa_flux_filename, 'r', format='NETCDF4') + + lwa = ncin1.variables['lwa'] + lwa = (np.array(lwa)) + + print(lwa.shape) + + z = np.zeros((120,360)) + for m in range(0, 120): + z[m,:] = lwa[m,49,:] # 49N LWA for June 1-30 + + zs = np.zeros((120,360)) # smoothed z # + + #### smoothing in longitude #### + n = 5 # smoothing width # + for m in range(0, 120): + zx = np.zeros(360) + zx[:] = z[m,:] + + for nn in range(-n, n+1): + zy = np.roll(zx,nn) + zs[m,:] = zs[m,:] + zy[:]/(2*n+1) + + zx[:] = zs[100,:]-zs[76,:] + zy[:] = 0 + + ncin1 = Dataset(lwa_flux_filename, 'r', format='NETCDF4') + ua1 = ncin1.variables['ua1'] + ua1 = (np.array(ua1)) + ncin2 = Dataset(lwa_flux_filename, 'r', format='NETCDF4') + ua2 = ncin2.variables['ua2'] + ua2 = (np.array(ua2)) + ncin3 = Dataset(lwa_flux_filename, 'r', format='NETCDF4') + ep1 = ncin3.variables['ep1'] + ep1 = (np.array(ep1)) + ncin4 = Dataset(lwa_flux_filename, 'r', format='NETCDF4') + ep2 = ncin4.variables['ep2'] + ep2 = (np.array(ep2)) + ncin5 = Dataset(lwa_flux_filename, 'r', format='NETCDF4') + ep3 = ncin5.variables['ep3'] + ep3 = (np.array(ep3)) + ncin6 = Dataset(lwa_flux_filename, 'r', format='NETCDF4') + ep4 = ncin6.variables['ep4'] + ep4 = (np.array(ep4)) + + f1 = np.zeros((120,360)) + f2 = np.zeros((120,360)) + f11 = np.zeros((120,360)) + f22 = np.zeros((120,360)) + + z1 = np.zeros((120,360)) + z2 = np.zeros((120,360)) + z3 = np.zeros((120,360)) + dt = 3600.*6. + a = 6378000. + dl = 2.*np.pi/360. + dp = 2.*np.pi/360. + # m = 76 is 20 June 2021 00 UTC + for m in range(0, 120): # m = 100 is 26 June 2021 00 UTC + z3[m,:] = ep4[m,49,:] + f1[m,:] = ua1[m,49,:]+ua2[m,49,:]+ep1[m,49,:] + j = 49 + phi = dp*j + const = 1./(2.*a*np.cos(phi)*dl) + z2[m,:]=const*(ep2[m,49,:]-ep3[m,49,:]) + f2[m,:] = 0.5*(ep2[m,49,:]+ep3[m,49,:])/np.cos(phi) + + for i in range(1, 359): + z1[m,i] = -const*(f1[m,i+1]-f1[m,i-1]) + + z1[m,0] = -const*(f1[m,1]-f1[m,359]) + z1[m,359] = -const*(f1[m,0]-f1[m,358]) + + + z1s = np.zeros((120,360)) # smoothed z1 # + z2s = np.zeros((120,360)) # smoothed z2 # + z3s = np.zeros((120,360)) # smoothed z3 # + f1s = np.zeros((120,360)) # smoothed f1 # + + #### smoothing in longitude #### + + for m in range(0, 120): + z1x = np.zeros(360) + z1x[:] = z1[m,:] + z2x = np.zeros(360) + z2x[:] = z2[m,:] + z3x = np.zeros(360) + z3x[:] = z3[m,:] + f1x = np.zeros(360) + f1x[:] = f1[m,:] + + nn = -n + for nn in range(-n, n+1): + z1y = np.roll(z1x,nn) + z1s[m,:] = z1s[m,:] + z1y[:]/(2*n+1) + z2y = np.roll(z2x,nn) + z2s[m,:] = z2s[m,:] + z2y[:]/(2*n+1) + z3y = np.roll(z3x,nn) + z3s[m,:] = z3s[m,:] + z3y[:]/(2*n+1) + f1y = np.roll(f1x,nn) + f1s[m,:] = f1s[m,:] + f1y[:]/(2*n+1) + + + + ##### Time integration test ###### + + lwa = np.zeros(360) + res = np.zeros((120,360)) + dlwa = np.zeros((120,360)) + gama = np.zeros((120,360)) + gama1 = np.zeros((120,360)) + gama2 = np.zeros((120,360)) + cgx = np.zeros((120,360)) + + m = 0 + for m in range(0, 119): + dlwa[m,:] = zs[m+1,:]-zs[m,:] + cgx[m,:] = (f1s[m,:]+f1s[m+1,:])/(zs[m,:]+zs[m+1,:]) + res[m,:] = (dlwa[m,:] - 0.5*dt*(z1s[m,:]+z1s[m+1,:]+z2s[m,:]+z2s[m+1,:]+z3s[m,:]+z3s[m+1,:]))/dt + gama[m,:] = 2.*res[m,:]/(zs[m,:]+zs[m+1,:]) + + + for m in range(76, 100): + lwa[:] = lwa[:] + 0.5*dt*(z1s[m,:]+z1s[m+1,:]+z2s[m,:]+z2s[m+1,:]+z3s[m,:]+z3s[m+1,:]) + dt*res[m,:] + m = m+1 + + j = 49 + phi = dp*j + const = 1./(2.*a*np.cos(phi)*dl) + lwa0 = np.zeros((120,360)) + lwa1 = np.zeros((120,360)) + + lwa0[:,:] = zs[:,:] + lwa1[:,:] = zs[:,:] + + + for m in range(76, 100): + + for i in range(1, 359): + d0 = -const*0.5*((lwa0[m,i+1]+lwa0[m+1,i+1])*cgx[m,i+1]-(lwa0[m,i-1]+lwa0[m+1,i-1])*cgx[m,i-1]) + d2 = 0.5*(z2s[m+1,i]+z2s[m,i]) + d3 = 0.5*(z3s[m+1,i]+z3s[m,i]) + d4 = gama[m,i]*0.5*(zs[m,i]+zs[m+1,i]) + lwa1[m+1,i] = lwa1[m,i]+dt*(d0+d2+d3+d4) + + d0 = -const*0.5*((lwa0[m,0]+lwa0[m+1,0])*cgx[m,0]-(lwa0[m,358]+lwa0[m+1,358])*cgx[m,358]) + d2 = 0.5*(z2s[m+1,359]+z2s[m,359]) + d3 = 0.5*(z3s[m+1,359]+z3s[m,359]) + d4 = gama[m,359]*0.5*(zs[m,359]+zs[m+1,359]) + lwa1[m+1,359] = lwa1[m,359]+dt*(d0+d2+d3+d4) + + d0 = -const*0.5*((lwa0[m,1]+lwa0[m+1,1])*cgx[m,1]-(lwa0[m,359]+lwa0[m+1,359])*cgx[m,359]) + d2 = 0.5*(z2s[m+1,0]+z2s[m,0]) + d3 = 0.5*(z3s[m+1,0]+z3s[m,0]) + d4 = gama[m,0]*0.5*(zs[m,0]+zs[m+1,0]) + lwa1[m+1,0] = lwa1[m,0]+dt*(d0+d2+d3+d4) + + + #### Modify gamma #### + + gama1[:,:] = gama[:,:] + gama2[:,:] = gama[:,:] + + for m in range(76, 100): + if((m >= 84) and (m <= 92)): + i = 200 + for i in range(200, 221): + if(gama[m,i] > 0): + gama1[m,i] = gama[m,i]*0.7 + gama2[m,i] = gama[m,i]*0. + + + gama1[:,:] = gama1[:,:]-gama[:,:] + gama2[:,:] = gama2[:,:]-gama[:,:] + + ##### Interpolate gama, gama1, cgx onto finer tiem mesh dt = 30 min ##### + + gamap = np.zeros((8640,360)) + gama1p = np.zeros((8640,360)) + gama2p = np.zeros((8640,360)) + cgxp = np.zeros((8640,360)) + forcep = np.zeros((8640,360)) + lwap = np.zeros((8640,360)) + + for m in range(76, 101): + for i in range(0, 360): + mm = (m-76)*360+i-180 + if(mm >= 0): + if(mm < 8640): + gamap[mm,:] = gama[m,:]*(i/360.)+gama[m-1,:]*(1.-i/360.) + gama1p[mm,:] = gama1[m,:]*(i/360.)+gama1[m-1,:]*(1.-i/360.) + gama2p[mm,:] = gama2[m,:]*(i/360.)+gama2[m-1,:]*(1.-i/360.) + cgxp[mm,:] = cgx[m,:]*(i/360.)+cgx[m-1,:]*(1.-i/360.) + + for m in range(76, 100): + for i in range(0, 360): + mm = (m-76)*360+i + forcep[mm,:] = (z2s[m+1,:]+z3s[m+1,:])*(i/360.)+(z2s[m,:]+z3s[m,:])*(1.-i/360.) + lwap[mm,:] = zs[m+1,:]*(i/360.)+zs[m,:]*(1.-i/360.) + + + ##### Time integration #### + + dt = 60. + j = 49 + phi = dp*j + const = 1./(2.*a*np.cos(phi)*dl) + al = 0. + diff = 200000. + dlwap = np.zeros((8641,360)) + dlwap2 = np.zeros((8641,360)) + + + for m in range(1, 8640): + for i in range(1, 359): + d1 = -const*(dlwap[m,i+1]*(cgxp[m,i+1]-al*(dlwap[m,i+1]+lwap[m,i+1]))-dlwap[m,i-1]*(cgxp[m,i-1]-al*(dlwap[m,i-1]+lwap[m,i-1]))) + d2 = gama1p[m,i]*dlwap[m,i] + d3 = gamap[m,i]*dlwap[m,i] + d4 = gama1p[m,i]*lwap[m,i] + d5 = diff*(dlwap[m-1,i+1]+dlwap[m-1,i-1]-2.*dlwap[m-1,i])/(a*a*dl*dl*np.cos(phi)*np.cos(phi)) + dlwap[m+1,i] = dlwap[m-1,i]+2.*dt*(d1+d2+d3+d4+d5) + + d12 = -const*(dlwap2[m,i+1]*(cgxp[m,i+1]-al*(dlwap2[m,i+1]+lwap[m,i+1]))-dlwap2[m,i-1]*(cgxp[m,i-1]-al*(dlwap2[m,i-1]+lwap[m,i-1]))) + d22 = gama2p[m,i]*dlwap2[m,i] + d32 = gamap[m,i]*dlwap2[m,i] + d42 = gama2p[m,i]*lwap[m,i] + d52 = diff*(dlwap2[m-1,i+1]+dlwap2[m-1,i-1]-2.*dlwap2[m-1,i])/(a*a*dl*dl*np.cos(phi)*np.cos(phi)) + dlwap2[m+1,i] = dlwap2[m-1,i]+2.*dt*(d12+d22+d32+d42+d52) + + d1 = -const*(dlwap[m,0]*(cgxp[m,0]-al*(dlwap[m,0]+lwap[m,0]))-dlwap[m,358]*(cgxp[m,358]-al*(dlwap[m,358]+lwap[m,358]))) + d2 = gama1p[m,359]*dlwap[m,359] + d3 = gamap[m,359]*dlwap[m,359] + d4 = gama1p[m,359]*lwap[m,359] + d5 = diff*(dlwap[m-1,0]+dlwap[m-1,358]-2.*dlwap[m-1,359])/(a*a*dl*dl*np.cos(phi)*np.cos(phi)) + dlwap[m+1,359] = dlwap[m-1,359]+2.*dt*(d1+d2+d3+d4+d5) + + d12 = -const*(dlwap2[m,0]*(cgxp[m,0]-al*(dlwap2[m,0]+lwap[m,0]))-dlwap2[m,358]*(cgxp[m,358]-al*(dlwap2[m,358]+lwap[m,358]))) + d22 = gama2p[m,359]*dlwap2[m,359] + d32 = gamap[m,359]*dlwap2[m,359] + d42 = gama2p[m,359]*lwap[m,359] + d52 = diff*(dlwap2[m-1,0]+dlwap2[m-1,358]-2.*dlwap2[m-1,359])/(a*a*dl*dl*np.cos(phi)*np.cos(phi)) + dlwap2[m+1,359] = dlwap[m-1,359]+2.*dt*(d12+d22+d32+d42+d52) + + d1 = -const*(dlwap[m,1]*(cgxp[m,1]-al*(dlwap[m,1]+lwap[m,1]))-dlwap[m,359]*(cgxp[m,359]-al*(dlwap[m,359]+lwap[m,359]))) + d2 = gama1p[m,0]*dlwap[m,0] + d3 = gamap[m,0]*dlwap[m,0] + d4 = gama1p[m,0]*lwap[m,0] + d5 = diff*(dlwap[m-1,1]+dlwap[m-1,359]-2.*dlwap[m-1,0])/(a*a*dl*dl*np.cos(phi)*np.cos(phi)) + dlwap[m+1,0] = dlwap[m-1,0]+2.*dt*(d1+d2+d3+d4+d5) + + d12 = -const*(dlwap2[m,1]*(cgxp[m,1]-al*(dlwap2[m,1]+lwap[m,1]))-dlwap2[m,359]*(cgxp[m,359]-al*(dlwap2[m,359]+lwap[m,359]))) + d22 = gama2p[m,0]*dlwap2[m,0] + d32 = gamap[m,0]*dlwap2[m,0] + d42 = gama2p[m,0]*lwap[m,0] + d52 = diff*(dlwap2[m-1,1]+dlwap2[m-1,359]-2.*dlwap2[m-1,0])/(a*a*dl*dl*np.cos(phi)*np.cos(phi)) + dlwap2[m+1,0] = dlwap2[m-1,0]+2.*dt*(d12+d22+d32+d42+d52) + + gm = np.zeros(360) + gm2 = np.zeros(360) + gm[:] = lwa[:]+dlwap[8640,:] + gm2[:] = lwa[:]+dlwap2[8640,:] + x = np.arange(0,360) + plt.rcParams.update({'font.size':14}) + fig = plt.figure(figsize=(8, 4)) + plt.xlim(140, 280) + plt.ylim(-80, 80) + plt.title('Column LWA Change 49$^\circ$N June 20 - 26 00 UTC') + plt.xlabel('Longitude') + plt.ylabel('$\Delta$LWA (m/s)') + ott = plt.plot(x, zx) + ott = plt.plot(x, zy, color='black', alpha = 0.3) + ott = plt.plot(x, gm, color='red', alpha = 0.5) + ott = plt.plot(x, gm2, 'r--', alpha = 0.5) + plt.savefig('dLWAp.png', bbox_inches='tight', dpi =600) + + diff --git a/scripts/nhn_grl2022/readme.md b/scripts/nhn_grl2022/readme.md new file mode 100644 index 0000000..3b9cf83 --- /dev/null +++ b/scripts/nhn_grl2022/readme.md @@ -0,0 +1,35 @@ +## Demo script for the analyses done in Neal et al. (submitted to GRL) + +This repo contains the sample script to reproduce the plots in Neal et al. "The 2021 Pacific Northwest heat wave and +associated blocking: Meteorology and the role of an upstream cyclone as a diabatic source of wave activity". + +To produce all the plots, required ERA5 data (6-hourly, 1° x 1° spatial resolution at all pressure levels +(37 available levels), from June 1-30, 2021): + +- geopotential height (z) +- zonal wind (u) +- meridional wind (v) +- temperature (t) +- 2-m temperature (t2m) +- net OLR (mtnlwrf) +- clear-sky OLR (mtnlwrfcs) +- total column water (tcw) +- total column water vapor (tcwv) +- sea level pressure (sp) + +The run script `sample_run_script.py` called functions in `QGField` object (see `hn2016_falwa/oopinterface.py`) +to compute LWA and fluxes, while `scripts/graph_plot_module.py` contains the graph plotting functions to reproduce +all the figures in the paper. + +Note that there is a difference in computing reference state in this most recent manuscript. In the past, and also +the analysis in NH18 Science, we used equator as the latitudinal boundary. In this current version of the script +(release 0.6.0), we use absolute vorticity at 5°N as boundary condition such that the solution is no longer sensitive +to fields at the equator, which improves the quality of the analysis. + +We are in a hurry submitting the manuscript, so this version of code (v0.6.0) has not been properly refactored yet. +If you have any questions, please [submit an issue](https://github.com/csyhuang/hn2016_falwa/issues) or email me. +I'll get back ASAP. + +Thanks! + +Clare (csyhuang@protonmail.com) \ No newline at end of file diff --git a/scripts/nhn_grl2022/sample_run_script.py b/scripts/nhn_grl2022/sample_run_script.py new file mode 100644 index 0000000..45feb74 --- /dev/null +++ b/scripts/nhn_grl2022/sample_run_script.py @@ -0,0 +1,158 @@ +import os +import sys +import numpy as np +from math import pi +from netCDF4 import Dataset +from hn2016_falwa.oopinterface import QGField + +sys.path.insert(0, os.getcwd()) +from graph_plot_module import plot_figure1a, plot_figure1b, plot_figure1c, plot_figure1d_2a, plot_figure3_and_S1, \ + plot_figure3e, plot_figure3f, plot_figure4, plot_figure5 + +data_dir = "grl2021_data/" +to_generate_data = True + +# --- Load the zonal wind and QGPV at 240hPa --- # +u_file = Dataset(data_dir + '2021_06_u.nc', mode='r') +v_file = Dataset(data_dir + '2021_06_v.nc', mode='r') +t_file = Dataset(data_dir + '2021_06_t.nc', mode='r') + +time_array = u_file.variables['time'][:] +time_units = u_file.variables['time'].units +time_calendar = u_file.variables['time'].calendar +ntimes = time_array.shape[0] + +print('Dimension of time: {}'.format(time_array.size)) + +# --- Longitude, latitude and pressure grid --- +xlon = u_file.variables['longitude'][:] + +# latitude has to be in ascending order +ylat = u_file.variables['latitude'][:] +if np.diff(ylat)[0]<0: + print('Flip ylat.') + ylat = ylat[::-1] + +# pressure level has to be in descending order (ascending height) +plev = u_file.variables['level'][:] +if np.diff(plev)[0]>0: + print('Flip plev.') + plev = plev[::-1] + +nlon = xlon.size +nlat = ylat.size +nlev = plev.size + +# --- Coordinates --- +clat = np.cos(np.deg2rad(ylat)) # cosine latitude +p0 = 1000. # surface pressure [hPa] +kmax = 97 # number of grid points for vertical extrapolation (dimension of height) +dz = 500. # differential height element +height = np.arange(0,kmax)*dz # pseudoheight [m] +dphi = np.diff(ylat)[0]*pi/180. # differential latitudinal element +dlambda = np.diff(xlon)[0]*pi/180. # differential latitudinal element +hh = 7000. # scale height +cp = 1004. # heat capacity of dry air +rr = 287. # gas constant +omega = 7.29e-5 # rotation rate of the earth +aa = 6.378e+6 # earth radius +prefactor = np.array([np.exp(-z/hh) for z in height[1:]]).sum() # integrated sum of density from the level + #just above the ground (z=1km) to aloft +npart = nlat # number of partitions to construct the equivalent latitude grids +maxits = 100000 # maximum number of iteration in the SOR solver to solve for reference state +tol = 1.e-5 # tolerance that define convergence of solution +rjac = 0.95 # spectral radius of the Jacobi iteration in the SOR solver. +nd = nlat//2+1 # (one plus) index of latitude grid point with value 0 deg + # This is to be input to fortran code. The index convention is different. + + +# --- Outputing files --- +output_fname = '2021-06-01_to_2021-06-30_output.nc' +print(f"output_fname: {output_fname}") + +# --- Generate analysis results --- +if to_generate_data: + output_file = Dataset(output_fname, 'w') + output_file.createDimension('latitude',nlat//2+1) + output_file.createDimension('longitude',nlon) + output_file.createDimension('time',ntimes+4) + lats = output_file.createVariable('latitude',np.dtype('float32').char,('latitude',)) # Define the coordinate variables + lons = output_file.createVariable('longitude',np.dtype('float32').char,('longitude',)) + times = output_file.createVariable('time',np.dtype('int').char,('time',)) + lats.units = 'degrees_north' + lons.units = 'degrees_east' + times.units = time_units + times.calendar = time_calendar + lats[:] = ylat[90:] + lons[:] = xlon + # times[:] = time_array + lwa_baro = output_file.createVariable('lwa', np.dtype('float32').char, ('time', 'latitude', 'longitude')) + lwa_baro.units = 'm/s' + ua1 = output_file.createVariable('ua1', np.dtype('float32').char, ('time', 'latitude', 'longitude')) + ua1.units = 'm/s' + ua2 = output_file.createVariable('ua2', np.dtype('float32').char, ('time', 'latitude', 'longitude')) + ua2.units = 'm/s' + ep1 = output_file.createVariable('ep1', np.dtype('float32').char, ('time', 'latitude', 'longitude')) + ep1.units = 'm/s' + ep2 = output_file.createVariable('ep2', np.dtype('float32').char, ('time', 'latitude', 'longitude')) + ep2.units = 'm/s' + ep3 = output_file.createVariable('ep3', np.dtype('float32').char, ('time', 'latitude', 'longitude')) + ep3.units = 'm/s' + ep4 = output_file.createVariable('ep4', np.dtype('float32').char, ('time', 'latitude', 'longitude')) + ep4.units = 'm/s' + + + # --- Compute LWA + fluxes and save the data into netCDF file --- + for tstep in range(ntimes): # or ntimes + + uu = u_file.variables['u'][tstep, ::-1, ::-1, :].data + vv = v_file.variables['v'][tstep, ::-1, ::-1, :].data + tt = t_file.variables['t'][tstep, ::-1, ::-1, :].data + + qgfield_object = QGField(xlon, ylat, plev, uu, vv, tt, kmax=kmax, dz=dz, eq_boundary_index=5) + + qgpv_temp, interpolated_u_temp, interpolated_v_temp, interpolated_avort_temp, interpolated_theta_temp, \ + static_stability_n, static_stability_s, tn0, ts0 = qgfield_object._interpolate_field_dirinv() + + qref, uref, tref, fawa, ubar, tbar = qgfield_object._compute_qref_fawa_and_bc() + + astarbaro, ubaro, urefbaro, ua1baro, ua2baro, ep1baro, ep2baro, ep3baro, ep4baro, astar1, astar2 = \ + qgfield_object._compute_lwa_flux_dirinv(qref, uref, tref, fawa, ubar, tbar) + + lwa_baro[tstep, :, :] = np.swapaxes(astarbaro, 0, 1) + ua1[tstep, :, :] = np.swapaxes(ua1baro, 0, 1) + ua2[tstep, :, :] = np.swapaxes(ua2baro, 0, 1) + ep1[tstep, :, :] = np.swapaxes(ep1baro, 0, 1) + ep2[tstep, :, :] = np.swapaxes(ep2baro, 0, 1) + ep3[tstep, :, :] = np.swapaxes(ep3baro, 0, 1) + ep4[tstep, :, :] = np.swapaxes(ep4baro, 0, 1) + + print(f'tstep = {tstep}/{ntimes}.') + output_file.close() + +# --- Graph plotting for GRL2021 --- + +# Location of data files +data_dir = "grl2021_data/" +z_filename = data_dir + "2021_06_z.nc" # geopotential height +u_filename = data_dir + "2021_06_u.nc" # u +v_filename = data_dir + "2021_06_v.nc" # v +t_filename = data_dir + "2021_06_t.nc" # temperature +t2m_filename = data_dir + "2021_06_2t.nc" # t2m +mtnlwrf_filename = data_dir + "2021_06_mtnlwrf.nc" # net OLR +mtnlwrfcs_filename = data_dir + "2021_06_mtnlwrfcs.nc" # OLR clear sky +tcw_filename = data_dir + "2021_06_tcw.nc" # total column water (kg/m^2) +tcwv_filename = data_dir + "2021_06_tcwv.nc" # total column water vapor (kg/m^2) +sp_filename = data_dir + "2021_06_sp.nc" # sea level pressure (hPa) +lwa_flux_filename = output_fname + +# Execute graph plotting functions +plot_figure1a(z_filename, u_filename, v_filename) +plot_figure1b(t_filename) +plot_figure1c(t2m_filename) +plot_figure1d_2a(t_filename) +plot_figure3_and_S1(lwa_flux_filename) +plot_figure3e(mtnlwrf_filename, mtnlwrfcs_filename) +plot_figure3f(tcw_filename, tcwv_filename, sp_filename) +plot_figure4(lwa_flux_filename) +plot_figure5(lwa_flux_filename) diff --git a/setup.py b/setup.py index 7ad1898..6ff51d3 100644 --- a/setup.py +++ b/setup.py @@ -8,7 +8,10 @@ Huang and Nakamura (2016, JAS): http://dx.doi.org/10.1175/JAS-D-15-0194.1 Huang and Nakamura (2017, GRL): http://onlinelibrary.wiley.com/doi/10.1002/2017GL073760/full Nakamura and Huang (2018, Science): https://doi.org/10.1126/science.aat0721 - The current version of the library handles calculation of FALWA in a spherical barotropic model and QGPV fields on isobaric surfaces. + Neal et al (submitted to GRL.) + + The current version of the library handles calculation of FALWA in a spherical barotropic model and QGPV fields on + isobaric surfaces. The functions in this library can compute the tracer equivalent-latitude relationship proposed in Nakamura (1996) (Also, see Allen and Nakamura (2003)) and the (zonal mean) @@ -31,19 +34,43 @@ sources=['hn2016_falwa/compute_lwa_and_barotropic_fluxes.f90'], f2py_options=['--quiet']) +# *** Extensions 4-9 are used by the direct inversion algorithm *** +ext4 = Extension(name='interpolate_fields_direct_inv', + sources=['hn2016_falwa/interpolate_fields_dirinv.f90'], + f2py_options=['--quiet']) + +ext5 = Extension(name='compute_qref_and_fawa_first', + sources=['hn2016_falwa/compute_qref_and_fawa_first.f90'], + f2py_options=['--quiet']) + +ext6 = Extension(name='matrix_b4_inversion', + sources=['hn2016_falwa/matrix_b4_inversion.f90'], + f2py_options=['--quiet']) + +ext7 = Extension(name='matrix_after_inversion', + sources=['hn2016_falwa/matrix_after_inversion.f90'], + f2py_options=['--quiet']) + +ext8 = Extension(name='upward_sweep', + sources=['hn2016_falwa/upward_sweep.f90'], + f2py_options=['--quiet']) + +ext9 = Extension(name='compute_flux_dirinv', + sources=['hn2016_falwa/compute_flux_dirinv.f90'], + f2py_options=['--quiet']) setup( name='hn2016_falwa', - version='0.5.0', + version='0.6.0', description='python package to compute finite-amplitude local wave activity (Huang and Nakamura 2016, JAS)', long_description=LONG_DESCRIPTION, url='https://github.com/csyhuang/hn2016_falwa', author='Clare S. Y. Huang', - author_email='csyhuang@uchicago.edu', + author_email='csyhuang@protonmail.com', license='MIT', packages=find_packages(), setup_requires=['pytest-runner'], tests_require=['pytest'], - ext_modules=[ext1, ext2, ext3], + ext_modules=[ext1, ext2, ext3, ext4, ext5, ext6, ext7, ext8, ext9], zip_safe=False )