diff --git a/fortran/eqthecalc.f90 b/fortran/eqthecalc.f90 index d805be9..33cad88 100644 --- a/fortran/eqthecalc.f90 +++ b/fortran/eqthecalc.f90 @@ -32,8 +32,7 @@ SUBROUTINE DEQTHECALC(qvp, tmk, prs, eth, miy, mjx, mkzh) REAL(KIND=8) :: tlcl INTEGER :: i, j, k - ! Note: removing temporaries did not improve performance for this algorithm. - !$OMP PARALLEL DO COLLAPSE(3) PRIVATE(i,j,k,q,t,p,e,tlcl) + !$OMP PARALLEL DO COLLAPSE(3) PRIVATE(i, j, k, q, t, p, e, tlcl) DO k = 1,mkzh DO j = 1,mjx DO i = 1,miy diff --git a/fortran/rip_cape.f90 b/fortran/rip_cape.f90 index d474702..d6a5fec 100644 --- a/fortran/rip_cape.f90 +++ b/fortran/rip_cape.f90 @@ -227,24 +227,27 @@ SUBROUTINE DPFCALC(prs, sfp, pf, mix, mjy, mkzh, ter_follow) INTEGER :: i,j,k + !$OMP PARALLEL DO COLLAPSE(3) DO j = 1,mjy - DO i = 1,mix - DO k = 1,mkzh - IF (k .EQ. mkzh) THEN - ! terrain-following data - IF (ter_follow .EQ. 1) THEN - pf(k,i,j) = sfp(i,j) - ! pressure-level data - ELSE - pf(k,i,j) = .5D0 * (3.D0*prs(k,i,j) - prs(k-1,i,j)) - END IF - ELSE - pf(k,i,j) = .5D0 * (prs(k+1,i,j) + prs(k,i,j)) - END IF - END DO - END DO + DO i = 1,mix + DO k = 1,mkzh + IF (k .EQ. mkzh) THEN + ! terrain-following data + IF (ter_follow .EQ. 1) THEN + pf(k,i,j) = sfp(i,j) + ! pressure-level data + ELSE + pf(k,i,j) = .5D0 * (3.D0*prs(k,i,j) - prs(k-1,i,j)) + END IF + ELSE + pf(k,i,j) = .5D0 * (prs(k+1,i,j) + prs(k,i,j)) + END IF + END DO + END DO END DO + !$OMP END PARALLEL DO + RETURN END SUBROUTINE DPFCALC @@ -351,19 +354,16 @@ SUBROUTINE DCAPECALC3D(prs,tmk,qvp,ght,ter,sfp,cape,cin,& ! calculated the pressure at full sigma levels (a set of pressure ! levels that bound the layers represented by the vertical grid points) - !CALL cpu_time(t1) - !CALL OMP_SET_NUM_THREADS(16) - - !$OMP PARALLEL DO + !$OMP PARALLEL DO COLLAPSE(3) DO j = 1,mjy - DO i = 1,mix - DO k = 1,mkzh - prs_new(k,i,j) = prs(i,j,k) - tmk_new(k,i,j) = tmk(i,j,k) - qvp_new(k,i,j) = qvp(i,j,k) - ght_new(k,i,j) = ght(i,j,k) - END DO - END DO + DO i = 1,mix + DO k = 1,mkzh + prs_new(k,i,j) = prs(i,j,k) + tmk_new(k,i,j) = tmk(i,j,k) + qvp_new(k,i,j) = qvp(i,j,k) + ght_new(k,i,j) = ght(i,j,k) + END DO + END DO END DO !$OMP END PARALLEL DO @@ -383,195 +383,193 @@ SUBROUTINE DCAPECALC3D(prs,tmk,qvp,ght,ter,sfp,cape,cin,& !$OMP facden, tmkenv, qvpenv, eslift, qvplift, buoy, benamin, & !$OMP benaccum, zrel, kmax, dz, elfound, & !$OMP kel, klfc, & - !$OMP i,j,k,kpar) + !$OMP i, j, k, kpar) DO j = 1,mjy - DO i = 1,mix - cape(i,j,1) = 0.D0 - cin(i,j,1) = 0.D0 + DO i = 1,mix + cape(i,j,1) = 0.D0 + cin(i,j,1) = 0.D0 -!!$OMP SIMD - DO kpar = 2, mkzh + !!$OMP SIMD + DO kpar = 2, mkzh - ! Calculate temperature and moisture properties of parcel - ! (note, qvppari and tmkpari already calculated above for 2d case.) + ! Calculate temperature and moisture properties of parcel + ! (note, qvppari and tmkpari already calculated above for 2d case.) - tlcl = TLCLC1/(LOG(tmk_new(kpar,i,j)**TLCLC2/(MAX(1.D-20,qvp_new(kpar,i,j)*prs_new(kpar,i,j)/ & - (EPS + qvp_new(kpar,i,j))))) - TLCLC3) + TLCLC4 + tlcl = TLCLC1/(LOG(tmk_new(kpar,i,j)**TLCLC2/(MAX(1.D-20,qvp_new(kpar,i,j)*prs_new(kpar,i,j)/ & + (EPS + qvp_new(kpar,i,j))))) - TLCLC3) + TLCLC4 - ethpari = tmk_new(kpar,i,j)*(1000.D0/prs_new(kpar,i,j))**(GAMMA*(1.D0 + GAMMAMD*qvp_new(kpar,i,j)))* & - EXP((THTECON1/tlcl - THTECON2)*qvp_new(kpar,i,j)*(1.D0 + THTECON3*qvp_new(kpar,i,j))) + ethpari = tmk_new(kpar,i,j)*(1000.D0/prs_new(kpar,i,j))**(GAMMA*(1.D0 + GAMMAMD*qvp_new(kpar,i,j)))* & + EXP((THTECON1/tlcl - THTECON2)*qvp_new(kpar,i,j)*(1.D0 + THTECON3*qvp_new(kpar,i,j))) - zlcl = ght_new(kpar,i,j) + (tmk_new(kpar,i,j) - tlcl)/(G/CP * (1.D0 + CPMD*qvp_new(kpar,i,j))) - - ! DO k = kpar,1,-1 - ! tmklift_new(k) = TONPSADIABAT(ethpari, prs_new(k,i,j), psadithte, psadiprs,& - ! psaditmk, GAMMA, errstat, errmsg) - ! END DO - ! Calculate buoyancy and relative height of lifted parcel at - ! all levels, and store in bottom up arrays. add a level at the lcl, - ! and at all points where buoyancy is zero. - ! - ! For arrays that go bottom to top - kk = 0 - ilcl = 0 - - IF (ght_new(kpar,i,j) .GE. zlcl) THEN - ! Initial parcel already saturated or supersaturated. - ilcl = 2 - klcl = 1 - END IF - - !!$OMP SIMD lastprivate(qvplift,tmklift,ghtlift,tvlift,tmkenv,qvpenv,tvenv,eslift,facden) - DO k = kpar,1,-1 - ! For arrays that go bottom to top - kk = kk + 1 - - ! Model level is below lcl - IF (ght_new(k,i,j) .LT. zlcl) THEN - tmklift = tmk_new(kpar,i,j) - G/(CP * (1.D0 + CPMD*qvp_new(kpar,i,j))) * (ght_new(k,i,j) - ght_new(kpar,i,j)) - tvenv = tmk_new(k,i,j)*(EPS + qvp_new(k,i,j))/(EPS*(1.D0 + qvp_new(k,i,j))) - tvlift = tmklift*(EPS + qvp_new(kpar,i,j))/(EPS*(1.D0 + qvp_new(kpar,i,j))) - ghtlift = ght_new(k,i,j) - ELSE IF (ght(i,j,k) .GE. zlcl .AND. ilcl .EQ. 0) THEN - ! This model level and previous model level straddle the lcl, - ! so first create a new level in the bottom-up array, at the lcl. - facden = 1.0/(ght_new(k,i,j) - ght_new(k+1,i,j)) - tmkenv = tmk_new(k+1,i,j)*((ght_new(k,i,j)-zlcl)*facden) + tmk_new(k,i,j)*((zlcl-ght_new(k+1,i,j))*facden) - qvpenv = qvp_new(k+1,i,j)*((ght_new(k,i,j)-zlcl)*facden) + qvp_new(k,i,j)*((zlcl-ght_new(k+1,i,j))*facden) - tvenv = tmkenv* (EPS + qvpenv) / (EPS * (1.D0 + qvpenv)) - tvlift = tlcl* (EPS + qvp_new(kpar,i,j)) / (EPS *(1.D0 + qvp_new(kpar,i,j))) - ghtlift = zlcl - ilcl = 1 - ELSE - tmklift = TONPSADIABAT(ethpari, prs_new(k,i,j), psadithte, psadiprs,& + zlcl = ght_new(kpar,i,j) + (tmk_new(kpar,i,j) - tlcl)/(G/CP * (1.D0 + CPMD*qvp_new(kpar,i,j))) + + ! Calculate buoyancy and relative height of lifted parcel at + ! all levels, and store in bottom up arrays. add a level at the lcl, + ! and at all points where buoyancy is zero. + ! + ! For arrays that go bottom to top + kk = 0 + ilcl = 0 + + IF (ght_new(kpar,i,j) .GE. zlcl) THEN + ! Initial parcel already saturated or supersaturated. + ilcl = 2 + klcl = 1 + END IF + + !!$OMP SIMD lastprivate(qvplift,tmklift,ghtlift,tvlift,tmkenv,qvpenv,tvenv,eslift,facden) + DO k = kpar,1,-1 + ! For arrays that go bottom to top + kk = kk + 1 + + ! Model level is below lcl + IF (ght_new(k,i,j) .LT. zlcl) THEN + tmklift = tmk_new(kpar,i,j) - G/(CP * (1.D0 + CPMD*qvp_new(kpar,i,j)))*& + (ght_new(k,i,j) - ght_new(kpar,i,j)) + tvenv = tmk_new(k,i,j)*(EPS + qvp_new(k,i,j))/(EPS*(1.D0 + qvp_new(k,i,j))) + tvlift = tmklift*(EPS + qvp_new(kpar,i,j))/(EPS*(1.D0 + qvp_new(kpar,i,j))) + ghtlift = ght_new(k,i,j) + ELSE IF (ght(i,j,k) .GE. zlcl .AND. ilcl .EQ. 0) THEN + ! This model level and previous model level straddle the lcl, + ! so first create a new level in the bottom-up array, at the lcl. + facden = 1.0/(ght_new(k,i,j) - ght_new(k+1,i,j)) + tmkenv = tmk_new(k+1,i,j)*((ght_new(k,i,j)-zlcl)*facden) + tmk_new(k,i,j)*& + ((zlcl-ght_new(k+1,i,j))*facden) + qvpenv = qvp_new(k+1,i,j)*((ght_new(k,i,j)-zlcl)*facden) + qvp_new(k,i,j)*& + ((zlcl-ght_new(k+1,i,j))*facden) + tvenv = tmkenv* (EPS + qvpenv) / (EPS * (1.D0 + qvpenv)) + tvlift = tlcl* (EPS + qvp_new(kpar,i,j)) / (EPS *(1.D0 + qvp_new(kpar,i,j))) + ghtlift = zlcl + ilcl = 1 + ELSE + tmklift = TONPSADIABAT(ethpari, prs_new(k,i,j), psadithte, psadiprs,& psaditmk, GAMMA, errstat, errmsg) - eslift = EZERO*EXP(ESLCON1*(tmklift - CELKEL)/(tmklift - ESLCON2)) - qvplift = EPS*eslift/(prs_new(k,i,j) - eslift) - tvenv = tmk_new(k,i,j) * (EPS + qvp_new(k,i,j)) / (EPS * (1.D0 + qvp_new(k,i,j))) - tvlift = tmklift*(EPS + qvplift) / (EPS * (1.D0 + qvplift)) - ghtlift = ght_new(k,i,j) - END IF - ! Buoyancy - buoy(kk) = G*(tvlift - tvenv)/tvenv - zrel(kk) = ghtlift - ght_new(kpar,i,j) - IF ((kk .GT. 1) .AND. (buoy(kk)*buoy(kk-1) .LT. 0.0D0)) THEN - ! Parcel ascent curve crosses sounding curve, so create a new level - ! in the bottom-up array at the crossing. - kk = kk + 1 - buoy(kk) = buoy(kk-1) - zrel(kk) = zrel(kk-1) - buoy(kk-1) = 0.D0 - zrel(kk-1) = zrel(kk-2) + buoy(kk-2)/& - (buoy(kk-2) - buoy(kk))*(zrel(kk) - zrel(kk-2)) - END IF - IF (ilcl .EQ. 1) THEN - klcl = kk - ilcl = 2 - CYCLE - END IF - - END DO + eslift = EZERO*EXP(ESLCON1*(tmklift - CELKEL)/(tmklift - ESLCON2)) + qvplift = EPS*eslift/(prs_new(k,i,j) - eslift) + tvenv = tmk_new(k,i,j) * (EPS + qvp_new(k,i,j)) / (EPS * (1.D0 + qvp_new(k,i,j))) + tvlift = tmklift*(EPS + qvplift) / (EPS * (1.D0 + qvplift)) + ghtlift = ght_new(k,i,j) + END IF + ! Buoyancy + buoy(kk) = G*(tvlift - tvenv)/tvenv + zrel(kk) = ghtlift - ght_new(kpar,i,j) + IF ((kk .GT. 1) .AND. (buoy(kk)*buoy(kk-1) .LT. 0.0D0)) THEN + ! Parcel ascent curve crosses sounding curve, so create a new level + ! in the bottom-up array at the crossing. + kk = kk + 1 + buoy(kk) = buoy(kk-1) + zrel(kk) = zrel(kk-1) + buoy(kk-1) = 0.D0 + zrel(kk-1) = zrel(kk-2) + buoy(kk-2)/& + (buoy(kk-2) - buoy(kk))*(zrel(kk) - zrel(kk-2)) + END IF + IF (ilcl .EQ. 1) THEN + klcl = kk + ilcl = 2 + CYCLE + END IF + + END DO - kmax = kk - ! IF (kmax .GT. 150) THEN - ! print *,'kmax got too big' - ! errstat = ALGERR - ! WRITE(errmsg, *) 'capecalc3d: kmax got too big. kmax=',kmax - ! RETURN - ! END IF - - ! If no lcl was found, set klcl to kmax. it is probably not really - ! at kmax, but this will make the rest of the routine behave - ! properly. - IF (ilcl .EQ. 0) klcl=kmax - - ! Get the accumulated buoyant energy from the parcel's starting - ! point, at all levels up to the top level. - benaccum(1) = 0.0D0 - benamin = 9d9 - DO k = 2,kmax - dz = zrel(k) - zrel(k-1) - benaccum(k) = benaccum(k-1) + .5D0*dz*(buoy(k-1) + buoy(k)) - IF (benaccum(k) .LT. benamin) THEN - benamin = benaccum(k) - END IF - END DO - ! Determine equilibrium level (el), which we define as the highest - ! level of non-negative buoyancy above the lcl. note, this may be - ! the top level if the parcel is still buoyant there. - - elfound = .FALSE. - DO k = kmax,klcl,-1 - IF (buoy(k) .GE. 0.D0) THEN - ! k of equilibrium level - kel = k - elfound = .TRUE. - EXIT - END IF - END DO - - ! If we got through that loop, then there is no non-negative - ! buoyancy above the lcl in the sounding. in these situations, - ! both cape and cin will be set to -0.1 j/kg. (see below about - ! missing values in v6.1.0). also, where cape is - ! non-zero, cape and cin will be set to a minimum of +0.1 j/kg, so - ! that the zero contour in either the cin or cape fields will - ! circumscribe regions of non-zero cape. - - ! In v6.1.0 of ncl, we added a _fillvalue attribute to the return - ! value of this function. at that time we decided to change -0.1 - ! to a more appropriate missing value, which is passed into this - ! routine as cmsg. - - IF (.NOT. elfound) THEN - !print *,'el not found' - cape(i,j,kpar) = cmsg - cin(i,j,kpar) = cmsg - klfc = kmax - CYCLE - END IF - - ! If there is an equilibrium level, then cape is positive. we'll - ! define the level of free convection (lfc) as the point below the - ! el, but at or above the lcl, where accumulated buoyant energy is a - ! minimum. the net positive area (accumulated buoyant energy) from - ! the lfc up to the el will be defined as the cape, and the net - ! negative area (negative of accumulated buoyant energy) from the - ! parcel starting point to the lfc will be defined as the convective - ! inhibition (cin). - - ! First get the lfc according to the above definition. - benamin = 9D9 - klfc = kmax - DO k = klcl,kel - IF (benaccum(k) .LT. benamin) THEN - benamin = benaccum(k) - klfc = k - END IF - END DO - - ! Now we can assign values to cape and cin - - cape(i,j,kpar) = MAX(benaccum(kel)-benamin, 0.1D0) - cin(i,j,kpar) = MAX(-benamin, 0.1D0) - - ! cin is uninteresting when cape is small (< 100 j/kg), so set - ! cin to -0.1 (see note about missing values in v6.1.0) in - ! that case. - - ! In v6.1.0 of ncl, we added a _fillvalue attribute to the return - ! value of this function. at that time we decided to change -0.1 - ! to a more appropriate missing value, which is passed into this - ! routine as cmsg. - - IF (cape(i,j,kpar) .LT. 100.D0) cin(i,j,kpar) = cmsg - - END DO - END DO + kmax = kk + ! IF (kmax .GT. 150) THEN + ! print *,'kmax got too big' + ! errstat = ALGERR + ! WRITE(errmsg, *) 'capecalc3d: kmax got too big. kmax=',kmax + ! RETURN + ! END IF + + ! If no lcl was found, set klcl to kmax. it is probably not really + ! at kmax, but this will make the rest of the routine behave + ! properly. + IF (ilcl .EQ. 0) klcl=kmax + + ! Get the accumulated buoyant energy from the parcel's starting + ! point, at all levels up to the top level. + benaccum(1) = 0.0D0 + benamin = 9d9 + DO k = 2,kmax + dz = zrel(k) - zrel(k-1) + benaccum(k) = benaccum(k-1) + .5D0*dz*(buoy(k-1) + buoy(k)) + IF (benaccum(k) .LT. benamin) THEN + benamin = benaccum(k) + END IF + END DO + ! Determine equilibrium level (el), which we define as the highest + ! level of non-negative buoyancy above the lcl. note, this may be + ! the top level if the parcel is still buoyant there. + + elfound = .FALSE. + DO k = kmax,klcl,-1 + IF (buoy(k) .GE. 0.D0) THEN + ! k of equilibrium level + kel = k + elfound = .TRUE. + EXIT + END IF + END DO + + ! If we got through that loop, then there is no non-negative + ! buoyancy above the lcl in the sounding. in these situations, + ! both cape and cin will be set to -0.1 j/kg. (see below about + ! missing values in v6.1.0). also, where cape is + ! non-zero, cape and cin will be set to a minimum of +0.1 j/kg, so + ! that the zero contour in either the cin or cape fields will + ! circumscribe regions of non-zero cape. + + ! In v6.1.0 of ncl, we added a _fillvalue attribute to the return + ! value of this function. at that time we decided to change -0.1 + ! to a more appropriate missing value, which is passed into this + ! routine as cmsg. + + IF (.NOT. elfound) THEN + !print *,'el not found' + cape(i,j,kpar) = cmsg + cin(i,j,kpar) = cmsg + klfc = kmax + CYCLE + END IF + + ! If there is an equilibrium level, then cape is positive. we'll + ! define the level of free convection (lfc) as the point below the + ! el, but at or above the lcl, where accumulated buoyant energy is a + ! minimum. the net positive area (accumulated buoyant energy) from + ! the lfc up to the el will be defined as the cape, and the net + ! negative area (negative of accumulated buoyant energy) from the + ! parcel starting point to the lfc will be defined as the convective + ! inhibition (cin). + + ! First get the lfc according to the above definition. + benamin = 9D9 + klfc = kmax + DO k = klcl,kel + IF (benaccum(k) .LT. benamin) THEN + benamin = benaccum(k) + klfc = k + END IF + END DO + + ! Now we can assign values to cape and cin + + cape(i,j,kpar) = MAX(benaccum(kel)-benamin, 0.1D0) + cin(i,j,kpar) = MAX(-benamin, 0.1D0) + + ! cin is uninteresting when cape is small (< 100 j/kg), so set + ! cin to -0.1 (see note about missing values in v6.1.0) in + ! that case. + + ! In v6.1.0 of ncl, we added a _fillvalue attribute to the return + ! value of this function. at that time we decided to change -0.1 + ! to a more appropriate missing value, which is passed into this + ! routine as cmsg. + + IF (cape(i,j,kpar) .LT. 100.D0) cin(i,j,kpar) = cmsg + + END DO + END DO END DO -!$OMP END PARALLEL DO - !CALL cpu_time(t2) - !print *,'Time taken in seconds ',(t2-t1) + !$OMP END PARALLEL DO + RETURN END SUBROUTINE DCAPECALC3D @@ -590,8 +588,6 @@ END SUBROUTINE DCAPECALC3D ! the cape and cin arrays. Also, LCL and LFC heights ! are put in the k=mkzh-1 and k=mkzh-2 slabs of the cin array. ! - - ! Important! The z-indexes must be arranged so that mkzh (max z-index) is the ! surface pressure. So, pressure must be ordered in ascending order before ! calling this routine. Other variables must be ordered the same (p,tk,q,z). @@ -686,14 +682,14 @@ SUBROUTINE DCAPECALC2D(prs,tmk,qvp,ght,ter,sfp,cape,cin,& !$OMP PARALLEL DO COLLAPSE(3) DO j = 1,mjy - DO i = 1,mix - DO k = 1,mkzh - prs_new(k,i,j) = prs(i,j,k) - tmk_new(k,i,j) = tmk(i,j,k) - qvp_new(k,i,j) = qvp(i,j,k) - ght_new(k,i,j) = ght(i,j,k) - END DO - END DO + DO i = 1,mix + DO k = 1,mkzh + prs_new(k,i,j) = prs(i,j,k) + tmk_new(k,i,j) = tmk(i,j,k) + qvp_new(k,i,j) = qvp(i,j,k) + ght_new(k,i,j) = ght(i,j,k) + END DO + END DO END DO !$OMP END PARALLEL DO @@ -711,274 +707,272 @@ SUBROUTINE DCAPECALC2D(prs,tmk,qvp,ght,ter,sfp,cape,cin,& RETURN END IF - !CALL OMP_SET_NUM_THREADS(16) - !nthreads = omp_get_num_threads() - - !$OMP PARALLEL DO COLLAPSE(2) PRIVATE(tlcl, ethpari, & !$OMP zlcl, kk, ilcl, klcl, tmklift, tvenv, tvlift, ghtlift, & !$OMP facden, tmkenv, qvpenv, eslift, qvplift, buoy, benamin, & !$OMP benaccum, zrel, kmax, dz, elfound, & !$OMP kel, klfc, pavg, p2, p1, totthe, totqvp, totprs, & - !$OMP i,j,k,kpar, kpar1, kpar2, qvppari, tmkpari,p, pup, pdn, q, th, & + !$OMP i, j, k, kpar, kpar1, kpar2, qvppari, tmkpari,p, pup, pdn, q, th, & !$OMP pp1, pp2, ethmax, eth_temp, klev) DO j = 1,mjy - DO i = 1,mix - cape(i,j,1) = 0.D0 - cin(i,j,1) = 0.D0 - ! find parcel with max theta-e in lowest 3 km agl. - ethmax = -1.D0 - eth_temp = -1.D0 - DO k = 1, mkzh - IF (ght_new(k,i,j)-ter(i,j) .LT. 3000.D0) THEN - tlcl = TLCLC1 / (LOG(tmk_new(k,i,j)**TLCLC2/& - (MAX(qvp_new(k,i,j), 1.d-15)*prs_new(k,i,j)/(EPS+MAX(qvp_new(k,i,j), 1.d-15))))-TLCLC3)+& - TLCLC4 - eth_temp(k) = tmk_new(k,i,j) * (1000.D0/prs_new(k,i,j))**& - (GAMMA*(1.D0 + GAMMAMD*(MAX(qvp_new(k,i,j), 1.d-15))))*& - EXP((THTECON1/tlcl - THTECON2)*(MAX(qvp_new(k,i,j), 1.d-15))*& - (1.D0 + THTECON3*(MAX(qvp_new(k,i,j), 1.d-15)))) - END IF - END DO - klev = mkzh - DO k = 1,mkzh - IF (eth_temp(k) .GT. ethmax) THEN - klev = k - ethmax = eth_temp(k) - END IF - END DO - - kpar1 = klev - kpar2 = klev - - - ! Establish average properties of that parcel - ! (over depth of approximately davg meters) - - !davg = 500.D0 - pavg = 500.D0 * prs_new(kpar1,i,j)*& - G/(RD*tvirtual(tmk_new(kpar1,i,j), qvp_new(kpar1,i,j))) - p2 = MIN(prs_new(kpar1,i,j)+.5d0*pavg, prsf(mkzh,i,j)) - p1 = p2 - pavg - totthe = 0.D0 - totqvp = 0.D0 - totprs = 0.D0 - DO k = mkzh,2,-1 - IF (prsf(k,i,j) .LE. p1) EXIT !GOTO 35 - IF (prsf(k-1,i,j) .GE. p2) CYCLE !GOTO 34 - p = prs_new(k,i,j) - pup = prsf(k,i,j) - pdn = prsf(k-1,i,j) - !q = MAX(qvp_new(k,i,j),1.D-15) - th = tmk_new(k,i,j)*(1000.D0/prs_new(k,i,j))**(GAMMA*(1.D0 + GAMMAMD*MAX(qvp_new(k,i,j),1.D-15))) - pp1 = MAX(p1,pdn) - pp2 = MIN(p2,pup) - IF (pp2 .GT. pp1) THEN - ! deltap = pp2 - pp1 - totqvp = totqvp + MAX(qvp_new(k,i,j),1.D-15)*(pp2 - pp1) - totthe = totthe + th*(pp2 - pp1) - totprs = totprs + (pp2 - pp1) - END IF - END DO - qvppari = totqvp/totprs - tmkpari = (totthe/totprs)*& - (prs_new(kpar1,i,j)/1000.D0)**(GAMMA*(1.D0+GAMMAMD*qvp_new(kpar1,i,j))) - -!CALL CPU_TIME(t3) - DO kpar = kpar1, kpar2 - - ! Calculate temperature and moisture properties of parcel - ! (note, qvppari and tmkpari already calculated above for 2d - ! case.) - - !prspari = prs_new(kpar,i,j) - !ghtpari = ght_new(kpar,i,j) - gammam = GAMMA * (1.D0 + GAMMAMD*qvppari) - cpm = CP * (1.D0 + CPMD*qvppari) - - e = MAX(1.D-20,qvppari*prs_new(kpar,i,j)/(EPS + qvppari)) - tlcl = TLCLC1/(LOG(tmkpari**TLCLC2/e) - TLCLC3) + TLCLC4 - ethpari = tmkpari*(1000.D0/prs_new(kpar,i,j))**(GAMMA*(1.D0 + GAMMAMD*qvppari))*& - EXP((THTECON1/tlcl - THTECON2)*qvppari*(1.D0 + THTECON3*qvppari)) - zlcl = ght_new(kpar,i,j) + (tmkpari - tlcl)/(G/cpm) - - ! Calculate buoyancy and relative height of lifted parcel at - ! all levels, and store in bottom up arrays. add a level at the - ! lcl, - ! and at all points where buoyancy is zero. - ! - ! - ! For arrays that go bottom to top - kk = 0 - ilcl = 0 - - IF (ght_new(kpar,i,j) .GE. zlcl) THEN - ! Initial parcel already saturated or supersaturated. - ilcl = 2 - klcl = 1 - END IF - - k = kpar - DO k = kpar,1,-1 - ! For arrays that go bottom to top - kk = kk + 1 - - ! Model level is below lcl - IF (ght_new(k,i,j) .LT. zlcl) THEN - tmklift = tmk_new(kpar,i,j) - G/(CP * (1.D0 + CPMD*qvp_new(kpar,i,j))) * (ght_new(k,i,j) - ght_new(kpar,i,j)) - tvenv = tmk_new(k,i,j)*(EPS + qvp_new(k,i,j))/(EPS*(1.D0 + qvp_new(k,i,j))) - tvlift = tmklift*(EPS + qvp_new(kpar,i,j))/(EPS*(1.D0 + qvp_new(kpar,i,j))) - ghtlift = ght_new(k,i,j) - ELSE IF (ght(i,j,k) .GE. zlcl .AND. ilcl .EQ. 0) THEN - ! This model level and previous model level straddle the lcl, - ! so first create a new level in the bottom-up array, at the lcl. - facden = 1/(ght_new(k,i,j) - ght_new(k+1,i,j)) - tmkenv = tmk_new(k+1,i,j)*((ght_new(k,i,j)-zlcl)*facden) + tmk_new(k,i,j)*((zlcl-ght_new(k+1,i,j))*facden) - qvpenv = qvp_new(k+1,i,j)*((ght_new(k,i,j)-zlcl)*facden) + qvp_new(k,i,j)*((zlcl-ght_new(k+1,i,j))*facden) - tvenv = tmkenv* (EPS + qvpenv) / (EPS * (1.D0 + qvpenv)) - tvlift = tlcl* (EPS + qvp_new(kpar,i,j)) / (EPS *(1.D0 + qvp_new(kpar,i,j))) - ghtlift = zlcl - ilcl = 1 - ELSE - tmklift = TONPSADIABAT(ethpari, prs_new(k,i,j), psadithte, psadiprs,& - psaditmk, GAMMA, errstat, errmsg) - eslift = EZERO*EXP(ESLCON1*(tmklift - CELKEL)/(tmklift - ESLCON2)) - qvplift = EPS*eslift/(prs_new(k,i,j) - eslift) - tvenv = tmk_new(k,i,j) * (EPS + qvp_new(k,i,j)) / (EPS * (1.D0 + qvp_new(k,i,j))) - tvlift = tmklift*(EPS + qvplift) / (EPS * (1.D0 + qvplift)) - ghtlift = ght_new(k,i,j) - END IF - ! Buoyancy - buoy(kk) = G*(tvlift - tvenv)/tvenv - zrel(kk) = ghtlift - ght_new(kpar,i,j) - IF ((kk .GT. 1) .AND. (buoy(kk)*buoy(kk-1) .LT. 0.0D0)) THEN - ! Parcel ascent curve crosses sounding curve, so create a new level - ! in the bottom-up array at the crossing. - kk = kk + 1 - buoy(kk) = buoy(kk-1) - zrel(kk) = zrel(kk-1) - buoy(kk-1) = 0.D0 - zrel(kk-1) = zrel(kk-2) + buoy(kk-2)/& - (buoy(kk-2) - buoy(kk))*(zrel(kk) - zrel(kk-2)) - END IF - IF (ilcl .EQ. 1) THEN - klcl = kk - ilcl = 2 - CYCLE - END IF - - END DO - - kmax = kk - ! IF (kmax .GT. 150) THEN - ! errstat = ALGERR - ! WRITE(errmsg, *) 'capecalc3d: kmax got too big. kmax=',kmax - ! RETURN - ! END IF - - ! If no lcl was found, set klcl to kmax. it is probably not - ! really - ! at kmax, but this will make the rest of the routine behave - ! properly. - IF (ilcl .EQ. 0) klcl=kmax - - ! Get the accumulated buoyant energy from the parcel's starting - ! point, at all levels up to the top level. - benaccum(1) = 0.0D0 - benamin = 9d9 - DO k = 2,kmax - dz = zrel(k) - zrel(k-1) - benaccum(k) = benaccum(k-1) + .5D0*dz*(buoy(k-1) + buoy(k)) - IF (benaccum(k) .LT. benamin) THEN - benamin = benaccum(k) - END IF - END DO - ! Determine equilibrium level (el), which we define as the highest - ! level of non-negative buoyancy above the lcl. note, this may be - ! the top level if the parcel is still buoyant there. - - elfound = .FALSE. - DO k = kmax,klcl,-1 - IF (buoy(k) .GE. 0.D0) THEN - ! k of equilibrium level - kel = k - elfound = .TRUE. - EXIT - END IF - END DO - - ! If we got through that loop, then there is no non-negative - ! buoyancy above the lcl in the sounding. in these situations, - ! both cape and cin will be set to -0.1 j/kg. (see below about - ! missing values in v6.1.0). also, where cape is - ! non-zero, cape and cin will be set to a minimum of +0.1 j/kg, so - ! that the zero contour in either the cin or cape fields will - ! circumscribe regions of non-zero cape. - ! In v6.1.0 of ncl, we added a _fillvalue attribute to the return - ! value of this function. at that time we decided to change -0.1 - ! to a more appropriate missing value, which is passed into this - ! routine as cmsg. - - IF (.NOT. elfound) THEN - cape(i,j,kpar) = cmsg - cin(i,j,kpar) = cmsg - klfc = kmax - CYCLE - END IF - - - ! If there is an equilibrium level, then cape is positive. - ! we'll - ! define the level of free convection (lfc) as the point below - ! the - ! el, but at or above the lcl, where accumulated buoyant energy - ! is a - ! minimum. the net positive area (accumulated buoyant energy) - ! from - ! the lfc up to the el will be defined as the cape, and the net - ! negative area (negative of accumulated buoyant energy) from - ! the - ! parcel starting point to the lfc will be defined as the - ! convective - ! inhibition (cin). - - ! First get the lfc according to the above definition. - benamin = 9D9 - klfc = kmax - DO k = klcl,kel - IF (benaccum(k) .LT. benamin) THEN - benamin = benaccum(k) - klfc = k - END IF - END DO - - ! Now we can assign values to cape and cin - - cape(i,j,kpar) = MAX(benaccum(kel)-benamin, 0.1D0) - cin(i,j,kpar) = MAX(-benamin, 0.1D0) - - ! cin is uninteresting when cape is small (< 100 j/kg), so set - ! cin to -0.1 (see note about missing values in v6.1.0) in - ! that case. - - ! In v6.1.0 of ncl, we added a _fillvalue attribute to the return - ! value of this function. at that time we decided to change -0.1 - ! to a more appropriate missing value, which is passed into this - ! routine as cmsg. - - IF (cape(i,j,kpar) .LT. 100.D0) cin(i,j,kpar) = cmsg - - END DO - - cape(i,j,mkzh) = cape(i,j,kpar1) - cin(i,j,mkzh) = cin(i,j,kpar1) - ! meters agl - cin(i,j,mkzh-1) = zrel(klcl) + ght_new(kpar,i,j) - ter(i,j) - ! meters agl - cin(i,j,mkzh-2) = zrel(klfc) + ght_new(kpar,i,j) - ter(i,j) - - END DO + DO i = 1,mix + cape(i,j,1) = 0.D0 + cin(i,j,1) = 0.D0 + ! find parcel with max theta-e in lowest 3 km agl. + ethmax = -1.D0 + eth_temp = -1.D0 + DO k = 1, mkzh + IF (ght_new(k,i,j)-ter(i,j) .LT. 3000.D0) THEN + tlcl = TLCLC1 / (LOG(tmk_new(k,i,j)**TLCLC2/& + (MAX(qvp_new(k,i,j), 1.d-15)*prs_new(k,i,j)/(EPS+MAX(qvp_new(k,i,j), 1.d-15))))-TLCLC3)+& + TLCLC4 + eth_temp(k) = tmk_new(k,i,j) * (1000.D0/prs_new(k,i,j))**& + (GAMMA*(1.D0 + GAMMAMD*(MAX(qvp_new(k,i,j), 1.d-15))))*& + EXP((THTECON1/tlcl - THTECON2)*(MAX(qvp_new(k,i,j), 1.d-15))*& + (1.D0 + THTECON3*(MAX(qvp_new(k,i,j), 1.d-15)))) + END IF + END DO + klev = mkzh + DO k = 1,mkzh + IF (eth_temp(k) .GT. ethmax) THEN + klev = k + ethmax = eth_temp(k) + END IF + END DO + + kpar1 = klev + kpar2 = klev + + + ! Establish average properties of that parcel + ! (over depth of approximately davg meters) + + !davg = 500.D0 + pavg = 500.D0 * prs_new(kpar1,i,j)*& + G/(RD*tvirtual(tmk_new(kpar1,i,j), qvp_new(kpar1,i,j))) + p2 = MIN(prs_new(kpar1,i,j)+.5d0*pavg, prsf(mkzh,i,j)) + p1 = p2 - pavg + totthe = 0.D0 + totqvp = 0.D0 + totprs = 0.D0 + DO k = mkzh,2,-1 + IF (prsf(k,i,j) .LE. p1) EXIT !GOTO 35 + IF (prsf(k-1,i,j) .GE. p2) CYCLE !GOTO 34 + p = prs_new(k,i,j) + pup = prsf(k,i,j) + pdn = prsf(k-1,i,j) + !q = MAX(qvp_new(k,i,j),1.D-15) + th = tmk_new(k,i,j)*(1000.D0/prs_new(k,i,j))**(GAMMA*(1.D0 + GAMMAMD*MAX(qvp_new(k,i,j),1.D-15))) + pp1 = MAX(p1,pdn) + pp2 = MIN(p2,pup) + IF (pp2 .GT. pp1) THEN + ! deltap = pp2 - pp1 + totqvp = totqvp + MAX(qvp_new(k,i,j),1.D-15)*(pp2 - pp1) + totthe = totthe + th*(pp2 - pp1) + totprs = totprs + (pp2 - pp1) + END IF + END DO + qvppari = totqvp/totprs + tmkpari = (totthe/totprs)*& + (prs_new(kpar1,i,j)/1000.D0)**(GAMMA*(1.D0+GAMMAMD*qvp_new(kpar1,i,j))) + + DO kpar = kpar1, kpar2 + + ! Calculate temperature and moisture properties of parcel + ! (note, qvppari and tmkpari already calculated above for 2d + ! case.) + + !prspari = prs_new(kpar,i,j) + !ghtpari = ght_new(kpar,i,j) + gammam = GAMMA * (1.D0 + GAMMAMD*qvppari) + cpm = CP * (1.D0 + CPMD*qvppari) + + e = MAX(1.D-20,qvppari*prs_new(kpar,i,j)/(EPS + qvppari)) + tlcl = TLCLC1/(LOG(tmkpari**TLCLC2/e) - TLCLC3) + TLCLC4 + ethpari = tmkpari*(1000.D0/prs_new(kpar,i,j))**(GAMMA*(1.D0 + GAMMAMD*qvppari))*& + EXP((THTECON1/tlcl - THTECON2)*qvppari*(1.D0 + THTECON3*qvppari)) + zlcl = ght_new(kpar,i,j) + (tmkpari - tlcl)/(G/cpm) + + ! Calculate buoyancy and relative height of lifted parcel at + ! all levels, and store in bottom up arrays. add a level at the + ! lcl, + ! and at all points where buoyancy is zero. + ! + ! + ! For arrays that go bottom to top + kk = 0 + ilcl = 0 + + IF (ght_new(kpar,i,j) .GE. zlcl) THEN + ! Initial parcel already saturated or supersaturated. + ilcl = 2 + klcl = 1 + END IF + + k = kpar + DO k = kpar,1,-1 + ! For arrays that go bottom to top + kk = kk + 1 + + ! Model level is below lcl + IF (ght_new(k,i,j) .LT. zlcl) THEN + tmklift = tmk_new(kpar,i,j) - G/(CP * (1.D0 + CPMD*qvp_new(kpar,i,j)))*& + (ght_new(k,i,j) - ght_new(kpar,i,j)) + tvenv = tmk_new(k,i,j)*(EPS + qvp_new(k,i,j))/(EPS*(1.D0 + qvp_new(k,i,j))) + tvlift = tmklift*(EPS + qvp_new(kpar,i,j))/(EPS*(1.D0 + qvp_new(kpar,i,j))) + ghtlift = ght_new(k,i,j) + ELSE IF (ght(i,j,k) .GE. zlcl .AND. ilcl .EQ. 0) THEN + ! This model level and previous model level straddle the lcl, + ! so first create a new level in the bottom-up array, at the lcl. + facden = 1/(ght_new(k,i,j) - ght_new(k+1,i,j)) + tmkenv = tmk_new(k+1,i,j)*((ght_new(k,i,j)-zlcl)*facden) + tmk_new(k,i,j)*& + ((zlcl-ght_new(k+1,i,j))*facden) + qvpenv = qvp_new(k+1,i,j)*((ght_new(k,i,j)-zlcl)*facden) + qvp_new(k,i,j)*& + ((zlcl-ght_new(k+1,i,j))*facden) + tvenv = tmkenv* (EPS + qvpenv) / (EPS * (1.D0 + qvpenv)) + tvlift = tlcl* (EPS + qvp_new(kpar,i,j)) / (EPS *(1.D0 + qvp_new(kpar,i,j))) + ghtlift = zlcl + ilcl = 1 + ELSE + tmklift = TONPSADIABAT(ethpari, prs_new(k,i,j), psadithte, psadiprs,& + psaditmk, GAMMA, errstat, errmsg) + eslift = EZERO*EXP(ESLCON1*(tmklift - CELKEL)/(tmklift - ESLCON2)) + qvplift = EPS*eslift/(prs_new(k,i,j) - eslift) + tvenv = tmk_new(k,i,j) * (EPS + qvp_new(k,i,j)) / (EPS * (1.D0 + qvp_new(k,i,j))) + tvlift = tmklift*(EPS + qvplift) / (EPS * (1.D0 + qvplift)) + ghtlift = ght_new(k,i,j) + END IF + ! Buoyancy + buoy(kk) = G*(tvlift - tvenv)/tvenv + zrel(kk) = ghtlift - ght_new(kpar,i,j) + IF ((kk .GT. 1) .AND. (buoy(kk)*buoy(kk-1) .LT. 0.0D0)) THEN + ! Parcel ascent curve crosses sounding curve, so create a new level + ! in the bottom-up array at the crossing. + kk = kk + 1 + buoy(kk) = buoy(kk-1) + zrel(kk) = zrel(kk-1) + buoy(kk-1) = 0.D0 + zrel(kk-1) = zrel(kk-2) + buoy(kk-2)/& + (buoy(kk-2) - buoy(kk))*(zrel(kk) - zrel(kk-2)) + END IF + IF (ilcl .EQ. 1) THEN + klcl = kk + ilcl = 2 + CYCLE + END IF + END DO + + kmax = kk + ! IF (kmax .GT. 150) THEN + ! errstat = ALGERR + ! WRITE(errmsg, *) 'capecalc3d: kmax got too big. kmax=',kmax + ! RETURN + ! END IF + + ! If no lcl was found, set klcl to kmax. it is probably not + ! really + ! at kmax, but this will make the rest of the routine behave + ! properly. + IF (ilcl .EQ. 0) klcl=kmax + + ! Get the accumulated buoyant energy from the parcel's starting + ! point, at all levels up to the top level. + benaccum(1) = 0.0D0 + benamin = 9d9 + DO k = 2,kmax + dz = zrel(k) - zrel(k-1) + benaccum(k) = benaccum(k-1) + .5D0*dz*(buoy(k-1) + buoy(k)) + IF (benaccum(k) .LT. benamin) THEN + benamin = benaccum(k) + END IF + END DO + ! Determine equilibrium level (el), which we define as the highest + ! level of non-negative buoyancy above the lcl. note, this may be + ! the top level if the parcel is still buoyant there. + + elfound = .FALSE. + DO k = kmax,klcl,-1 + IF (buoy(k) .GE. 0.D0) THEN + ! k of equilibrium level + kel = k + elfound = .TRUE. + EXIT + END IF + END DO + + ! If we got through that loop, then there is no non-negative + ! buoyancy above the lcl in the sounding. in these situations, + ! both cape and cin will be set to -0.1 j/kg. (see below about + ! missing values in v6.1.0). also, where cape is + ! non-zero, cape and cin will be set to a minimum of +0.1 j/kg, so + ! that the zero contour in either the cin or cape fields will + ! circumscribe regions of non-zero cape. + ! In v6.1.0 of ncl, we added a _fillvalue attribute to the return + ! value of this function. at that time we decided to change -0.1 + ! to a more appropriate missing value, which is passed into this + ! routine as cmsg. + + IF (.NOT. elfound) THEN + cape(i,j,kpar) = cmsg + cin(i,j,kpar) = cmsg + klfc = kmax + CYCLE + END IF + + + ! If there is an equilibrium level, then cape is positive. + ! we'll + ! define the level of free convection (lfc) as the point below + ! the + ! el, but at or above the lcl, where accumulated buoyant energy + ! is a + ! minimum. the net positive area (accumulated buoyant energy) + ! from + ! the lfc up to the el will be defined as the cape, and the net + ! negative area (negative of accumulated buoyant energy) from + ! the + ! parcel starting point to the lfc will be defined as the + ! convective + ! inhibition (cin). + + ! First get the lfc according to the above definition. + benamin = 9D9 + klfc = kmax + DO k = klcl,kel + IF (benaccum(k) .LT. benamin) THEN + benamin = benaccum(k) + klfc = k + END IF + END DO + + ! Now we can assign values to cape and cin + + cape(i,j,kpar) = MAX(benaccum(kel)-benamin, 0.1D0) + cin(i,j,kpar) = MAX(-benamin, 0.1D0) + + ! cin is uninteresting when cape is small (< 100 j/kg), so set + ! cin to -0.1 (see note about missing values in v6.1.0) in + ! that case. + + ! In v6.1.0 of ncl, we added a _fillvalue attribute to the return + ! value of this function. at that time we decided to change -0.1 + ! to a more appropriate missing value, which is passed into this + ! routine as cmsg. + + IF (cape(i,j,kpar) .LT. 100.D0) cin(i,j,kpar) = cmsg + + END DO + + cape(i,j,mkzh) = cape(i,j,kpar1) + cin(i,j,mkzh) = cin(i,j,kpar1) + ! meters agl + cin(i,j,mkzh-1) = zrel(klcl) + ght_new(kpar,i,j) - ter(i,j) + ! meters agl + cin(i,j,mkzh-2) = zrel(klfc) + ght_new(kpar,i,j) - ter(i,j) + + END DO END DO -!$OMP END PARALLEL DO + !$OMP END PARALLEL DO + RETURN END SUBROUTINE DCAPECALC2D diff --git a/src/wrf/config.py b/src/wrf/config.py index c7f28b2..556f632 100644 --- a/src/wrf/config.py +++ b/src/wrf/config.py @@ -2,35 +2,64 @@ from __future__ import (absolute_import, division, print_function, unicode_literals) from threading import local +import wrapt _local_config = local() -_local_config.xarray_enabled = True -_local_config.cartopy_enabled = True -_local_config.basemap_enabled = True -_local_config.pyngl_enabled = True -_local_config.cache_size = 20 - -try: - from xarray import DataArray -except ImportError: - _local_config.xarray_enabled = False - -try: - from cartopy import crs -except ImportError: - _local_config.cartopy_enabled = False + +def _init_local(): + global _local_config -try: - from mpl_toolkits.basemap import Basemap -except ImportError: - _local_config.basemap_enabled = False + _local_config.xarray_enabled = True + _local_config.cartopy_enabled = True + _local_config.basemap_enabled = True + _local_config.pyngl_enabled = True + _local_config.cache_size = 20 + _local_config.initialized = True -try: - from Ngl import Resources -except ImportError: - _local_config.pyngl_enabled = False + try: + from xarray import DataArray + except ImportError: + _local_config.xarray_enabled = False + + try: + from cartopy import crs + except ImportError: + _local_config.cartopy_enabled = False + + try: + from mpl_toolkits.basemap import Basemap + except ImportError: + _local_config.basemap_enabled = False + + try: + from Ngl import Resources + except ImportError: + _local_config.pyngl_enabled = False + +# Initialize the main thread's configuration +_init_local() + +def init_local(): + """A decorator that initializes thread local data if necessary.""" + @wrapt.decorator + def func_wrapper(wrapped, instance, args, kwargs): + global _local_config + try: + init = _local_config.init + except AttributeError: + _init_local() + else: + if not init: + _init_local() + + return wrapped(*args, **kwargs) + + return func_wrapper + + +@init_local() def xarray_enabled(): """Return True if xarray is installed and enabled. @@ -43,18 +72,21 @@ def xarray_enabled(): return _local_config.xarray_enabled +@init_local() def disable_xarray(): """Disable xarray.""" global _local_config _local_config.xarray_enabled = False - + +@init_local() def enable_xarray(): """Enable xarray.""" global _local_config _local_config.xarray_enabled = True - + +@init_local() def cartopy_enabled(): """Return True if cartopy is installed and enabled. @@ -67,18 +99,21 @@ def cartopy_enabled(): return _local_config.cartopy_enabled +@init_local() def enable_cartopy(): """Enable cartopy.""" global _local_config _local_config.cartopy_enabled = True - + +@init_local() def disable_cartopy(): """Disable cartopy.""" global _local_config _local_config.cartopy_enabled = True - + +@init_local() def basemap_enabled(): """Return True if basemap is installed and enabled. @@ -91,17 +126,21 @@ def basemap_enabled(): return _local_config.basemap_enabled +@init_local() def enable_basemap(): """Enable basemap.""" global _local_config _local_config.basemap_enabled = True - + +@init_local() def disable_basemap(): """Disable basemap.""" global _local_config _local_config.basemap_enabled = True - + + +@init_local() def pyngl_enabled(): """Return True if pyngl is installed and enabled. @@ -114,18 +153,21 @@ def pyngl_enabled(): return _local_config.pyngl_enabled +@init_local() def enable_pyngl(): """Enable pyngl.""" global _local_config _local_config.pyngl_enabled = True - + +@init_local() def disable_pyngl(): """Disable pyngl.""" global _local_config _local_config.pyngl_enabled = True - + +@init_local() def set_cache_size(size): """Set the maximum number of items that the threadlocal cache can retain. @@ -143,7 +185,8 @@ def set_cache_size(size): global _local_config _local_config.cache_size = size - + +@init_local() def get_cache_size(): """Return the maximum number of items that the threadlocal cache can retain. diff --git a/src/wrf/util.py b/src/wrf/util.py index 6c520e1..fff88ec 100644 --- a/src/wrf/util.py +++ b/src/wrf/util.py @@ -209,11 +209,30 @@ def _generator_copy(gen): module = getmodule(gen.gi_frame) if module is not None: - res = module.get(funcname)(**argvals.locals) + try: + try: + argd = {key:argvals.locals[key] for key in argvals.args} + res = module.get(funcname)(**argd) + except AttributeError: + res = getattr(module, funcname)(**argd) + except: + # This is the old way it used to work, but it looks like this was + # fixed by Python. + try: + res = module.get(funcname)(**argvals.locals) + except AttributeError: + res = getattr(module, funcname)(**argvals.locals) else: # Created in jupyter or the python interpreter import __main__ - res = getattr(__main__, funcname)(**argvals.locals) + + try: + argd = {key:argvals.locals[key] for key in argvals.args} + res = getattr(__main__, funcname)(**argd) + except: + # This was the old way it used to work, but appears to have + # been fixed by Python. + res = getattr(__main__, funcname)(**argvals.locals) return res @@ -2583,26 +2602,6 @@ def get_proj_params(wrfin):#, timeidx=0, varname=None): "DX", "DY")) return proj_params -# multitime = is_multi_time_req(timeidx) -# if not multitime: -# time_idx_or_slice = timeidx -# else: -# time_idx_or_slice = slice(None) -# -# if varname is not None: -# if not is_coordvar(varname): -# coord_names = getattr(wrfin.variables[varname], -# "coordinates").split() -# lon_coord = coord_names[0] -# lat_coord = coord_names[1] -# else: -# lat_coord, lon_coord = get_coord_pairs(varname) -# else: -# lat_coord, lon_coord = latlon_coordvars(wrfin.variables) -# -# return (wrfin.variables[lat_coord][time_idx_or_slice,:], -# wrfin.variables[lon_coord][time_idx_or_slice,:], -# proj_params) def from_args(func, argnames, *args, **kwargs): @@ -2936,18 +2935,31 @@ def psafilepath(): return os.path.join(os.path.dirname(__file__), "data", "psadilookup.dat") -def get_id(obj): - """Return the object id. +def get_filepath(obj): - The object id is used as a caching key for various routines. If the + try: + path = obj.filepath() + except AttributeError: + try: + path = obj.file.path + except: + raise ValueError("file contains no path information") + + return path + +def get_id(obj, prefix=''): + """Return the cache id. + + The cache id is used as a caching key for various routines. If the object type is a mapping, then the result will also be a - mapping of each key to the object id for the value. Otherwise, only the - object id is returned. + mapping of each key to the object id for the value. Args: obj (:obj:`object`): Any object type. + prefix (:obj:`str`): A string to help with recursive calls. + Returns: :obj:`int` or :obj:`dict`: If the *obj* parameter is not a mapping, @@ -2955,12 +2967,18 @@ def get_id(obj): key to the object id for the value is returned. """ + if not is_multi_file(obj): + return hash(prefix + get_filepath(obj)) + + # For sequences, the hashing string will be the list ID and the + # path for the first file in the sequence if not is_mapping(obj): - return id(obj) + _next = next(iter(obj)) + return get_id(_next, prefix + str(id(obj))) # For each key in the mapping, recursively call get_id until # until a non-mapping is found - return {key : get_id(val) for key,val in viewitems(obj)} + return {key : get_id(val, prefix) for key,val in viewitems(obj)} def geo_bounds(var=None, wrfin=None, varname=None, timeidx=0, method="cat", diff --git a/test/cachetest.py b/test/cachetest.py index cc7efee..fad5e88 100644 --- a/test/cachetest.py +++ b/test/cachetest.py @@ -2,7 +2,10 @@ from __future__ import (absolute_import, division, print_function, unicode_literals) from threading import Thread -from Queue import Queue +try: + from Queue import Queue +except ImportError: + from queue import Queue from collections import OrderedDict import unittest as ut @@ -62,4 +65,4 @@ class CacheTest(ut.TestCase): if __name__ == "__main__": - ut.main() \ No newline at end of file + ut.main() diff --git a/test/generator_test.py b/test/generator_test.py new file mode 100644 index 0000000..6a45d89 --- /dev/null +++ b/test/generator_test.py @@ -0,0 +1,17 @@ +from __future__ import (absolute_import, division, print_function, unicode_literals) + +from wrf import getvar +from netCDF4 import Dataset as nc +#ncfile = nc("/Users/ladwig/Documents/wrf_files/wrfout_d01_2016-02-25_18_00_00") +ncfile = nc("/Users/ladwig/Documents/wrf_files/wrfout_d01_2016-10-07_00_00_00") + +def gen_seq(): + wrfseq = [ncfile, ncfile, ncfile] + for wrf in wrfseq: + yield wrf + +p_gen = getvar(gen_seq(), "P", method="join") + +print(p_gen) +del p_gen +