diff --git a/fortran/rip_cape.f90 b/fortran/rip_cape.f90 index 08cf9a0..945b98e 100644 --- a/fortran/rip_cape.f90 +++ b/fortran/rip_cape.f90 @@ -37,7 +37,8 @@ END FUNCTION TVIRTUAL REAL(KIND=8) FUNCTION TONPSADIABAT(thte, prs, psadithte, psadiprs, psaditmk, gamma,& errstat, errmsg) USE wrf_constants, ONLY : ALGERR - +!$OMP DECLARE SIMD (TONPSADIABAT) +!!uniform(thte,prs,psadithte,psadiprs,psaditmk) !f2py threadsafe !f2py intent(in,out) :: cape, cin @@ -57,8 +58,9 @@ REAL(KIND=8) FUNCTION TONPSADIABAT(thte, prs, psadithte, psadiprs, psaditmk, gam REAL(KIND=8) :: fracjt2 REAL(KIND=8) :: fracip REAL(KIND=8) :: fracip2 - - INTEGER :: ip, ipch, jt, jtch + + INTEGER :: l1, h1, mid1, rang1, l2, h2, mid2, rang2 + INTEGER :: ip, jt ! This function gives the temperature (in K) on a moist adiabat ! (specified by thte in K) given pressure in hPa. It uses a @@ -78,22 +80,53 @@ REAL(KIND=8) FUNCTION TONPSADIABAT(thte, prs, psadithte, psadiprs, psaditmk, gam ! Otherwise, look for the given thte/prs point in the lookup table. jt = -1 - DO jtch = 1, 150-1 - IF (thte .GE. psadithte(jtch) .AND. thte .LT. psadithte(jtch+1)) THEN - jt = jtch - EXIT - !GO TO 213 - END IF + l1 = 1 + h1 = 149 + rang1 = h1 - l1 + mid1 = (h1 + l1) / 2 + DO WHILE(rang1 .GT. 1) + if(thte .GE. psadithte(mid1)) then + l1 = mid1 + else + h1 = mid1 + end if + rang1 = h1 - l1 + mid1 = (h1 + l1) / 2 END DO - - ip = -1 - DO ipch = 1, 150-1 - IF (prs .LE. psadiprs(ipch) .AND. prs .GT. psadiprs(ipch+1)) THEN - ip = ipch - EXIT - !GO TO 215 - END IF + jt = l1 + + ! DO jtch = 1, 150-1 + ! IF (thte .GE. psadithte(jtch) .AND. thte .LT. psadithte(jtch+1)) THEN + ! jt = jtch + ! EXIT + ! !GO TO 213 + ! END IF + ! END DO + + ip = -1 + l2 = 1 + h2 = 149 + rang2 = h2 - l2 + mid2 = (h2 + l2) / 2 + DO WHILE(rang2 .GT. 1) + if(prs .LE. psadiprs(mid2)) then + l2 = mid2 + else + h2 = mid2 + end if + rang2 = h2 - l2 + mid2 = (h2 + l2) / 2 END DO + ip = l2 + + ! ip = -1 + ! DO ipch = 1, 150-1 + ! IF (prs .LE. psadiprs(ipch) .AND. prs .GT. psadiprs(ipch+1)) THEN + ! ip = ipch + ! EXIT + ! !GO TO 215 + ! END IF + ! END DO IF (jt .EQ. -1 .OR. ip .EQ. -1) THEN ! Set the error and return @@ -187,28 +220,26 @@ END SUBROUTINE DLOOKUP_TABLE ! level is above. SUBROUTINE DPFCALC(prs, sfp, pf, miy, mjx, mkzh, ter_follow) - REAL(KIND=8), DIMENSION(miy,mjx,mkzh), INTENT(IN) :: prs + REAL(KIND=8), DIMENSION(mkzh,miy,mjx), INTENT(IN) :: prs REAL(KIND=8), DIMENSION(miy,mjx), INTENT(IN) :: sfp - REAL(KIND=8), DIMENSION(miy,mjx,mkzh), INTENT(OUT) :: pf + REAL(KIND=8), DIMENSION(mkzh,miy,mjx), INTENT(OUT) :: pf INTEGER, INTENT(IN) :: ter_follow,miy,mjx,mkzh INTEGER :: i,j,k - - ! do j=1,mjx-1 Artifact of MM5 + DO j = 1,mjx - ! do i=1,miy-1 staggered grid DO i = 1,miy DO k = 1,mkzh IF (k .EQ. mkzh) THEN ! terrain-following data IF (ter_follow .EQ. 1) THEN - pf(i,j,k) = sfp(i,j) + pf(k,i,j) = sfp(i,j) ! pressure-level data ELSE - pf(i,j,k) = .5D0 * (3.D0*prs(i,j,k) - prs(i,j,k-1)) + pf(k,i,j) = .5D0 * (3.D0*prs(k,i,j) - prs(k-1,i,j)) END IF ELSE - pf(i,j,k) = .5D0 * (prs(i,j,k+1) + prs(i,j,k)) + pf(k,i,j) = .5D0 * (prs(k+1,i,j) + prs(k,i,j)) END IF END DO END DO @@ -224,17 +255,338 @@ END SUBROUTINE DPFCALC ! ! !DESCRIPTION: ! -! If i3dflag=1, this routine calculates CAPE and CIN (in m**2/s**2, +! This routine calculates CAPE and CIN (in m**2/s**2, ! or J/kg) for every grid point in the entire 3D domain (treating -! each grid point as a parcel). If i3dflag=0, then it -! calculates CAPE and CIN only for the parcel with max theta-e in +! each grid point as a parcel). +! + + +! 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). + +! Also, be advised that missing data values are not checked during the computation. +! Also also, Pressure must be hPa + +! NCLFORTSTART +SUBROUTINE DCAPECALC3D(prs,tmk,qvp,ght,ter,sfp,cape,cin,& + cmsg,miy,mjx,mkzh,ter_follow,& + psafile, errstat, errmsg) + USE wrf_constants, ONLY : ALGERR, CELKEL, G, EZERO, ESLCON1, ESLCON2, & + EPS, RD, CP, GAMMA, CPMD, RGASMD, GAMMAMD, TLCLC1, & + TLCLC2, TLCLC3, TLCLC4, THTECON1, THTECON2, THTECON3 + + !USE omp_lib + IMPLICIT NONE + + !f2py threadsafe + !f2py intent(in,out) :: cape, cin + + INTEGER, INTENT(IN) :: miy, mjx, mkzh, ter_follow + REAL(KIND=8), DIMENSION(miy,mjx,mkzh), INTENT(IN) :: prs + REAL(KIND=8), DIMENSION(miy,mjx,mkzh), INTENT(IN) :: tmk + REAL(KIND=8), DIMENSION(miy,mjx,mkzh), INTENT(IN) :: qvp + REAL(KIND=8), DIMENSION(miy,mjx,mkzh), INTENT(IN) :: ght + REAL(KIND=8), DIMENSION(miy,mjx), INTENT(IN) :: ter + REAL(KIND=8), DIMENSION(miy,mjx), INTENT(IN) ::sfp + REAL(KIND=8), DIMENSION(miy,mjx,mkzh), INTENT(OUT) :: cape + REAL(KIND=8), DIMENSION(miy,mjx,mkzh), INTENT(OUT) :: cin + REAL(KIND=8), INTENT(IN) :: cmsg + CHARACTER(LEN=*), INTENT(IN) :: psafile + INTEGER, INTENT(INOUT) :: errstat + CHARACTER(LEN=*), INTENT(INOUT) :: errmsg + +! NCLFORTEND + + ! local variables + INTEGER :: i, j, k, ilcl, kel, kk, klcl, klev, klfc, kmax, kpar + REAL(KIND=8) :: tlcl, zlcl + REAL(KIND=8) :: ethpari, qvppari, tmkpari + REAL(KIND=8) :: facden, qvplift, tmklift, tvenv, tvlift, ghtlift + REAL(KIND=8) :: eslift, tmkenv, qvpenv, tonpsadiabat + REAL(KIND=8) :: benamin, dz + REAL(KIND=8), DIMENSION(150) :: buoy, zrel, benaccum + REAL(KIND=8), DIMENSION(mkzh,miy,mjx) :: prsf + REAL(KIND=8), DIMENSION(150) :: psadithte, psadiprs + REAL(KIND=8), DIMENSION(150,150) :: psaditmk + LOGICAL :: elfound + REAL :: t1,t2 + + REAL(KIND=8), DIMENSION(mkzh,miy,mjx) :: prs_new + REAL(KIND=8), DIMENSION(mkzh,miy,mjx) :: tmk_new + REAL(KIND=8), DIMENSION(mkzh,miy,mjx) :: qvp_new + REAL(KIND=8), DIMENSION(mkzh,miy,mjx) :: ght_new + + ! To remove compiler warnings + tmkpari = 0 + qvppari = 0 + klev = 0 + klcl = 0 + kel = 0 + + + ! the comments were taken from a mark stoelinga email, 23 apr 2007, + ! in response to a user getting the "outside of lookup table bounds" + ! error message. + + ! tmkpari - initial temperature of parcel, k + ! values of 300 okay. (not sure how much from this you can stray.) + + ! prspari - initial pressure of parcel, hpa + ! values of 980 okay. (not sure how much from this you can stray.) + + ! thtecon1, thtecon2, thtecon3 + ! these are all constants, the first in k and the other two have + ! no units. values of 3376, 2.54, and 0.81 were stated as being + ! okay. + + ! tlcl - the temperature at the parcel's lifted condensation level, k + ! should be a reasonable atmospheric temperature around 250-300 k + ! (398 is "way too high") + + ! qvppari - the initial water vapor mixing ratio of the parcel, + ! kg/kg (should range from 0.000 to 0.025) + ! + + ! 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 + DO i = 1,mjx + DO j = 1,miy + DO k = 1,mkzh + prs_new(k,j,i) = prs(j,i,k) + tmk_new(k,j,i) = tmk(j,i,k) + qvp_new(k,j,i) = qvp(j,i,k) + ght_new(k,j,i) = ght(j,i,k) + END DO + END DO + END DO +!$OMP END PARALLEL DO + + CALL DPFCALC(prs_new, sfp, prsf, miy, mjx, mkzh, ter_follow) + + ! before looping, set lookup table for getting temperature on + ! a pseudoadiabat. + + CALL DLOOKUP_TABLE(psadithte, psadiprs, psaditmk, psafile, errstat, errmsg) + + IF (errstat .NE. 0) THEN + RETURN + END IF + +!$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, & +!$OMP i,j,k,kpar) + DO j = 1,mjx + DO i = 1,miy + cape(i,j,1) = 0.d0 + cin(i,j,1) = 0.d0 + +!$OMP SIMD + DO kpar = 2, mkzh + + ! 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 + + 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,& + 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 + ! 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) + RETURN +END SUBROUTINE DCAPECALC3D + +!====================================================================== +! +! !IROUTINE: capecalc2d -- Calculate CAPE and CIN +! +! !DESCRIPTION: +! +! Calculates CAPE and CIN only for the parcel with max theta-e in ! the column, (i.e. something akin to Colman's MCAPE). By "parcel", ! we mean a 500-m deep parcel, with actual temperature and moisture ! averaged over that depth. ! -! In the case of i3dflag=0, ! CAPE and CIN are 2D fields that are placed in the k=mkzh slabs of -! the cape and cin arrays. Also, if i3dflag=0, LCL and LFC heights +! 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. ! @@ -243,23 +595,25 @@ END SUBROUTINE DPFCALC ! 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). -! Also, be advised that missing data values are not checked during the computation. +! Also, be advised that missing data values are not checked during the +! computation. ! Also also, Pressure must be hPa ! NCLFORTSTART -SUBROUTINE DCAPECALC3D(prs,tmk,qvp,ght,ter,sfp,cape,cin,& - cmsg,miy,mjx,mkzh,i3dflag,ter_follow,& +SUBROUTINE DCAPECALC2D(prs,tmk,qvp,ght,ter,sfp,cape,cin,& + cmsg,miy,mjx,mkzh,ter_follow,& psafile, errstat, errmsg) USE wrf_constants, ONLY : ALGERR, CELKEL, G, EZERO, ESLCON1, ESLCON2, & EPS, RD, CP, GAMMA, CPMD, RGASMD, GAMMAMD, TLCLC1, & TLCLC2, TLCLC3, TLCLC4, THTECON1, THTECON2, THTECON3 + !USE omp_lib IMPLICIT NONE !f2py threadsafe !f2py intent(in,out) :: cape, cin - INTEGER, INTENT(IN) :: miy, mjx, mkzh, i3dflag, ter_follow + INTEGER, INTENT(IN) :: miy, mjx, mkzh, ter_follow REAL(KIND=8), DIMENSION(miy,mjx,mkzh), INTENT(IN) :: prs REAL(KIND=8), DIMENSION(miy,mjx,mkzh), INTENT(IN) :: tmk REAL(KIND=8), DIMENSION(miy,mjx,mkzh), INTENT(IN) :: qvp @@ -275,26 +629,35 @@ SUBROUTINE DCAPECALC3D(prs,tmk,qvp,ght,ter,sfp,cape,cin,& ! NCLFORTEND + ! local variables INTEGER :: i, j, k, ilcl, kel, kk, klcl, klev, klfc, kmax, kpar, kpar1, kpar2 - REAL(KIND=8) :: davg, ethmax, q, t, p, e, eth, tlcl, zlcl + REAL(KIND=8) :: ethmax, q, p, e, tlcl, zlcl REAL(KIND=8) :: pavg, tvirtual, p1, p2, pp1, pp2, th, totthe, totqvp, totprs - REAL(KIND=8) :: cpm, deltap, ethpari, gammam, ghtpari, qvppari, prspari, tmkpari - REAL(KIND=8) :: facden, fac1, fac2, qvplift, tmklift, tvenv, tvlift, ghtlift + REAL(KIND=8) :: cpm, deltap, ethpari, gammam, qvppari, tmkpari + REAL(KIND=8) :: facden, qvplift, tmklift, tvenv, tvlift, ghtlift, fac1, fac2 REAL(KIND=8) :: eslift, tmkenv, qvpenv, tonpsadiabat REAL(KIND=8) :: benamin, dz, pup, pdn REAL(KIND=8), DIMENSION(150) :: buoy, zrel, benaccum - REAL(KIND=8), DIMENSION(miy,mjx,mkzh) :: prsf + REAL(KIND=8), DIMENSION(mkzh,miy,mjx) :: prsf REAL(KIND=8), DIMENSION(150) :: psadithte, psadiprs REAL(KIND=8), DIMENSION(150,150) :: psaditmk LOGICAL :: elfound + INTEGER :: nthreads + REAL(KIND=8), DIMENSION(mkzh) :: eth_temp + REAL(KIND=8), DIMENSION(mkzh,miy,mjx) :: prs_new + REAL(KIND=8), DIMENSION(mkzh,miy,mjx) :: tmk_new + REAL(KIND=8), DIMENSION(mkzh,miy,mjx) :: qvp_new + REAL(KIND=8), DIMENSION(mkzh,miy,mjx) :: ght_new ! To remove compiler warnings + errstat = 0 tmkpari = 0 qvppari = 0 klev = 0 klcl = 0 kel = 0 + deltap = 0 ! the comments were taken from a mark stoelinga email, 23 apr 2007, @@ -320,10 +683,23 @@ SUBROUTINE DCAPECALC3D(prs,tmk,qvp,ght,ter,sfp,cape,cin,& ! kg/kg (should range from 0.000 to 0.025) ! +!$OMP PARALLEL DO + DO i = 1,mjx + DO j = 1,miy + DO k = 1,mkzh + prs_new(k,j,i) = prs(j,i,k) + tmk_new(k,j,i) = tmk(j,i,k) + qvp_new(k,j,i) = qvp(j,i,k) + ght_new(k,j,i) = ght(j,i,k) + END DO + END DO + END DO +!$OMP END PARALLEL DO + + ! calculated the pressure at full sigma levels (a set of pressure ! levels that bound the layers represented by the vertical grid points) - - CALL DPFCALC(prs, sfp, prsf, miy, mjx, mkzh, ter_follow) + CALL DPFCALC(prs_new, sfp, prsf, miy, mjx, mkzh, ter_follow) ! before looping, set lookup table for getting temperature on ! a pseudoadiabat. @@ -334,148 +710,144 @@ SUBROUTINE DCAPECALC3D(prs,tmk,qvp,ght,ter,sfp,cape,cin,& RETURN END IF - ! do j=1,mjx-1 + !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, qvppari, tmkpari,p, pup, pdn, q, th, & +!$OMP pp1, pp2) DO j = 1,mjx - ! do i=1,miy-1 DO i = 1,miy 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 - IF (i3dflag .EQ. 1) THEN - kpar1 = 2 - kpar2 = mkzh - ELSE - ! find parcel with max theta-e in lowest 3 km agl. - ethmax = -1.d0 - DO k = mkzh,1,-1 - IF (ght(i,j,k)-ter(i,j) .LT. 3000.d0) THEN - q = MAX(qvp(i,j,k), 1.d-15) - t = tmk(i,j,k) - p = prs(i,j,k) - e = q*p/(EPS + q) - tlcl = TLCLC1 / (LOG(t**TLCLC2/e)-TLCLC3) + TLCLC4 - eth = t * (1000.d0/p)**(GAMMA*(1.d0 + GAMMAMD*q))*& - EXP((THTECON1/tlcl - THTECON2)*q*(1.d0 + THTECON3*q)) - IF (eth .GT. ethmax) THEN - klev = k - ethmax = eth - END IF - END IF - END DO - kpar1 = klev - kpar2 = klev - - ! Establish average properties of that parcel - ! (over depth of approximately davg meters) - - ! davg=.1 - davg = 500.d0 - pavg = davg*prs(i,j,kpar1)*& - G/(RD*tvirtual(tmk(i,j,kpar1), qvp(i,j,kpar1))) - p2 = MIN(prs(i,j,kpar1)+.5d0*pavg, prsf(i,j,mkzh)) - p1 = p2 - pavg - totthe = 0.D0 - totqvp = 0.D0 - totprs = 0.D0 - DO k = mkzh,2,-1 - IF (prsf(i,j,k) .LE. p1) EXIT !GOTO 35 - IF (prsf(i,j,k-1) .GE. p2) CYCLE !GOTO 34 - p = prs(i,j,k) - pup = prsf(i,j,k) - pdn = prsf(i,j,k-1) - q = MAX(qvp(i,j,k),1.D-15) - th = tmk(i,j,k)*(1000.D0/p)**(GAMMA*(1.D0 + GAMMAMD*q)) - pp1 = MAX(p1,pdn) - pp2 = MIN(p2,pup) - IF (pp2 .GT. pp1) THEN - deltap = pp2 - pp1 - totqvp = totqvp + q*deltap - totthe = totthe + th*deltap - totprs = totprs + deltap - END IF -! 34 CONTINUE - END DO -! 35 CONTINUE - qvppari = totqvp/totprs - tmkpari = (totthe/totprs)*& - (prs(i,j,kpar1)/1000.D0)**(GAMMA*(1.D0+GAMMAMD*qvp(i,j,kpar1))) - END IF - + 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.) + ! (note, qvppari and tmkpari already calculated above for 2d + ! case.) - IF (i3dflag .EQ. 1) THEN - qvppari = qvp(i,j,kpar) - tmkpari = tmk(i,j,kpar) - END IF - prspari = prs(i,j,kpar) - ghtpari = ght(i,j,kpar) + !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*prspari/(EPS + 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/prspari)**(GAMMA*(1.D0 + GAMMAMD*qvppari))*& + ethpari = tmkpari*(1000.D0/prs_new(kpar,i,j))**(GAMMA*(1.D0 + GAMMAMD*qvppari))*& EXP((THTECON1/tlcl - THTECON2)*qvppari*(1.D0 + THTECON3*qvppari)) - zlcl = ghtpari + (tmkpari - tlcl)/(G/cpm) + 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, + ! 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 (ghtpari .GE. zlcl) THEN + IF (ght_new(kpar,i,j) .GE. zlcl) THEN ! Initial parcel already saturated or supersaturated. ilcl = 2 klcl = 1 END IF k = kpar - DO WHILE (k .GE. 1)!k = kpar, 1, -1 - !DO k = kpar, 1, -1 - ! For arrays that go bottom to top -! 33 kk = kk + 1 + DO k = kpar,1,-1 + ! For arrays that go bottom to top kk = kk + 1 ! Model level is below lcl - IF (ght(i,j,k) .LT. zlcl) THEN - qvplift = qvppari - tmklift = tmkpari - G/cpm*(ght(i,j,k) - ghtpari) - tvenv = tvirtual(tmk(i,j,k), qvp(i,j,k)) - tvlift = tvirtual(tmklift, qvplift) - ghtlift = ght(i,j,k) + 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. - tmklift = tlcl - qvplift = qvppari - facden = ght(i,j,k) - ght(i,j,k+1) - fac1 = (zlcl-ght(i,j,k+1))/facden - fac2 = (ght(i,j,k)-zlcl)/facden - tmkenv = tmk(i,j,k+1)*fac2 + tmk(i,j,k)*fac1 - qvpenv = qvp(i,j,k+1)*fac2 + qvp(i,j,k)*fac1 - tvenv = tvirtual(tmkenv, qvpenv) - tvlift = tvirtual(tmklift, qvplift) + 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(i,j,k), psadithte, psadiprs,& + 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(i,j,k) - eslift) - tvenv = tvirtual(tmk(i,j,k), qvp(i,j,k)) - tvlift = tvirtual(tmklift, qvplift) - ghtlift = ght(i,j,k) + 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 - ghtpari - + 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. @@ -486,25 +858,23 @@ SUBROUTINE DCAPECALC3D(prs,tmk,qvp,ght,ter,sfp,cape,cin,& 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 - !GOTO 33 CYCLE END IF - k = k - 1 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 + ! 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 @@ -520,7 +890,6 @@ SUBROUTINE DCAPECALC3D(prs,tmk,qvp,ght,ter,sfp,cape,cin,& 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. @@ -532,7 +901,6 @@ SUBROUTINE DCAPECALC3D(prs,tmk,qvp,ght,ter,sfp,cape,cin,& kel = k elfound = .TRUE. EXIT - !GOTO 50 END IF END DO @@ -543,14 +911,11 @@ SUBROUTINE DCAPECALC3D(prs,tmk,qvp,ght,ter,sfp,cape,cin,& ! 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 + ! 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. - ! cape(i,j,kpar) = -0.1D0 - ! cin(i,j,kpar) = -0.1D0 IF (.NOT. elfound) THEN cape(i,j,kpar) = cmsg cin(i,j,kpar) = cmsg @@ -558,17 +923,20 @@ SUBROUTINE DCAPECALC3D(prs,tmk,qvp,ght,ter,sfp,cape,cin,& CYCLE END IF -! GOTO 102 -! 50 CONTINUE - - ! 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 + ! 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 + ! 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. @@ -595,23 +963,19 @@ SUBROUTINE DCAPECALC3D(prs,tmk,qvp,ght,ter,sfp,cape,cin,& ! 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) = -0.1D0 IF (cape(i,j,kpar) .LT. 100.D0) cin(i,j,kpar) = cmsg -! 102 CONTINUE END DO - IF (i3dflag .EQ. 0) THEN 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) + ghtpari - ter(i,j) + cin(i,j,mkzh-1) = zrel(klcl) + ght_new(kpar,i,j) - ter(i,j) ! meters agl - cin(i,j,mkzh-2) = zrel(klfc) + ghtpari - ter(i,j) + cin(i,j,mkzh-2) = zrel(klfc) + ght_new(kpar,i,j) - ter(i,j) - ENDIF END DO END DO - +!$OMP END PARALLEL DO RETURN -END SUBROUTINE DCAPECALC3D +END SUBROUTINE DCAPECALC2D diff --git a/fortran/wrf_rip_phys_routines.f90 b/fortran/wrf_rip_phys_routines.f90 index 350fb75..0b307f2 100644 --- a/fortran/wrf_rip_phys_routines.f90 +++ b/fortran/wrf_rip_phys_routines.f90 @@ -50,13 +50,15 @@ SUBROUTINE WETBULBCALC(prs, tmk, qvp, twb, nx, ny, nz, psafile, errstat, errmsg) !NCLEND INTEGER :: i, j, k - INTEGER :: jtch, jt, ipch, ip + INTEGER :: jt, ip REAL(KIND=8) :: q, t, p, e, tlcl, eth REAL(KIND=8) :: fracip, fracip2, fracjt, fracjt2 REAL(KIND=8), DIMENSION(150) :: PSADITHTE, PSADIPRS REAL(KIND=8), DIMENSION(150,150) :: PSADITMK REAL(KIND=8) :: tonpsadiabat + INTEGER :: l1, h1, mid1, rang1, l2, h2, mid2, rang2 + !INTEGER :: ip, ipch, jt, jtch ! Before looping, set lookup table for getting temperature on ! a pseudoadiabat. @@ -92,21 +94,53 @@ SUBROUTINE WETBULBCALC(prs, tmk, qvp, twb, nx, ny, nz, psafile, errstat, errmsg) tonpsadiabat = eth*(p/1000.)**GAMMA ELSE ! Otherwise, look for the given thte/prs point in the lookup table. - jt=-1 - DO jtch=1,150-1 - IF (eth .GE. PSADITHTE(jtch) .AND. eth .LT. PSADITHTE(jtch+1)) THEN - jt = jtch - EXIT - ENDIF - END DO - - ip=-1 - DO ipch=1,150-1 - IF (p .LE. PSADIPRS(ipch) .AND. p .GT. PSADIPRS(ipch+1)) THEN - ip = ipch - EXIT - ENDIF + jt = -1 + l1 = 1 + h1 = 149 + rang1 = h1 - l1 + mid1 = (h1 + l1) / 2 + DO WHILE(rang1 .GT. 1) + if(eth .GE. psadithte(mid1)) then + l1 = mid1 + else + h1 = mid1 + end if + rang1 = h1 - l1 + mid1 = (h1 + l1) / 2 END DO + jt = l1 + +! jt=-1 +! DO jtch=1,150-1 +! IF (eth .GE. PSADITHTE(jtch) .AND. eth .LT. PSADITHTE(jtch+1)) THEN +! jt = jtch +! EXIT +! ENDIF +! END DO + + ip = -1 + l2 = 1 + h2 = 149 + rang2 = h2 - l2 + mid2 = (h2 + l2) / 2 + DO WHILE(rang2 .GT. 1) + if(p .LE. psadiprs(mid2)) then + l2 = mid2 + else + h2 = mid2 + end if + rang2 = h2 - l2 + mid2 = (h2 + l2) / 2 + END DO + ip = l2 + +! ip=-1 +! DO ipch=1,150-1 +! IF (p .LE. PSADIPRS(ipch) .AND. p .GT. PSADIPRS(ipch+1)) THEN +! ip = ipch +! EXIT +! ENDIF +! END DO IF (jt .EQ. -1 .OR. ip .EQ. -1) THEN errstat = ALGERR diff --git a/fortran/wrffortran.pyf b/fortran/wrffortran.pyf index bf99a1e..eed6320 100644 --- a/fortran/wrffortran.pyf +++ b/fortran/wrffortran.pyf @@ -64,16 +64,38 @@ python module _wrffortran ! in character*(*) intent(inout) :: errmsg end subroutine dlookup_table subroutine dpfcalc(prs,sfp,pf,miy,mjx,mkzh,ter_follow) ! in :_wrffortran:rip_cape.f90 + real(kind=8) dimension(mkzh,miy,mjx),intent(in) :: prs + real(kind=8) dimension(miy,mjx),intent(in),depend(miy,mjx) :: sfp + real(kind=8) dimension(mkzh,miy,mjx),intent(out),depend(mkzh,miy,mjx) :: pf + integer, optional,intent(in),check(shape(prs,1)==miy),depend(prs) :: miy=shape(prs,1) + integer, optional,intent(in),check(shape(prs,2)==mjx),depend(prs) :: mjx=shape(prs,2) + integer, optional,intent(in),check(shape(prs,0)==mkzh),depend(prs) :: mkzh=shape(prs,0) + integer intent(in) :: ter_follow + end subroutine dpfcalc + subroutine dcapecalc3d(prs,tmk,qvp,ght,ter,sfp,cape,cin,cmsg,miy,mjx,mkzh,ter_follow,psafile,errstat,errmsg) ! in :_wrffortran:rip_cape.f90 + threadsafe + use omp_lib + use wrf_constants, only: tlclc2,gamma,tlclc1,rgasmd,tlclc4,g,tlclc3,thtecon3,eps,rd,cpmd,celkel,gammamd,eslcon2,eslcon1,cp,thtecon1,algerr,ezero,thtecon2 real(kind=8) dimension(miy,mjx,mkzh),intent(in) :: prs + real(kind=8) dimension(miy,mjx,mkzh),intent(in),depend(miy,mjx,mkzh) :: tmk + real(kind=8) dimension(miy,mjx,mkzh),intent(in),depend(miy,mjx,mkzh) :: qvp + real(kind=8) dimension(miy,mjx,mkzh),intent(in),depend(miy,mjx,mkzh) :: ght + real(kind=8) dimension(miy,mjx),intent(in),depend(miy,mjx) :: ter real(kind=8) dimension(miy,mjx),intent(in),depend(miy,mjx) :: sfp - real(kind=8) dimension(miy,mjx,mkzh),intent(out),depend(miy,mjx,mkzh) :: pf + real(kind=8) dimension(miy,mjx,mkzh),intent(out,in),depend(miy,mjx,mkzh) :: cape + real(kind=8) dimension(miy,mjx,mkzh),intent(out,in),depend(miy,mjx,mkzh) :: cin + real(kind=8) intent(in) :: cmsg integer, optional,intent(in),check(shape(prs,0)==miy),depend(prs) :: miy=shape(prs,0) integer, optional,intent(in),check(shape(prs,1)==mjx),depend(prs) :: mjx=shape(prs,1) integer, optional,intent(in),check(shape(prs,2)==mkzh),depend(prs) :: mkzh=shape(prs,2) integer intent(in) :: ter_follow - end subroutine dpfcalc - subroutine dcapecalc3d(prs,tmk,qvp,ght,ter,sfp,cape,cin,cmsg,miy,mjx,mkzh,i3dflag,ter_follow,psafile,errstat,errmsg) ! in :_wrffortran:rip_cape.f90 + character*(*) intent(in) :: psafile + integer intent(inout) :: errstat + character*(*) intent(inout) :: errmsg + end subroutine dcapecalc3d + subroutine dcapecalc2d(prs,tmk,qvp,ght,ter,sfp,cape,cin,cmsg,miy,mjx,mkzh,ter_follow,psafile,errstat,errmsg) ! in :_wrffortran:rip_cape.f90 threadsafe + use omp_lib use wrf_constants, only: tlclc2,gamma,tlclc1,rgasmd,tlclc4,g,tlclc3,thtecon3,eps,rd,cpmd,celkel,gammamd,eslcon2,eslcon1,cp,thtecon1,algerr,ezero,thtecon2 real(kind=8) dimension(miy,mjx,mkzh),intent(in) :: prs real(kind=8) dimension(miy,mjx,mkzh),intent(in),depend(miy,mjx,mkzh) :: tmk @@ -87,12 +109,11 @@ python module _wrffortran ! in integer, optional,intent(in),check(shape(prs,0)==miy),depend(prs) :: miy=shape(prs,0) integer, optional,intent(in),check(shape(prs,1)==mjx),depend(prs) :: mjx=shape(prs,1) integer, optional,intent(in),check(shape(prs,2)==mkzh),depend(prs) :: mkzh=shape(prs,2) - integer intent(in) :: i3dflag integer intent(in) :: ter_follow character*(*) intent(in) :: psafile integer intent(inout) :: errstat character*(*) intent(inout) :: errmsg - end subroutine dcapecalc3d + end subroutine dcapecalc2d subroutine dcloudfrac(pres,rh,lowc,midc,highc,nz,ns,ew) ! in :_wrffortran:wrf_cloud_fracf.f90 threadsafe real(kind=8) dimension(ew,ns,nz),intent(in) :: pres diff --git a/src/wrf/extension.py b/src/wrf/extension.py index 5f3fcea..90bb8cc 100755 --- a/src/wrf/extension.py +++ b/src/wrf/extension.py @@ -7,8 +7,8 @@ from .constants import Constants from ._wrffortran import (dcomputetk, dinterp3dz, dinterp2dxy, dinterp1d, dcomputeseaprs, dfilter2d, dcomputerh, dcomputeuvmet, - dcomputetd, dcapecalc3d, dcloudfrac, wrfcttcalc, - calcdbz, dcalrelhl, dcalcuh, dcomputepv, + dcomputetd, dcapecalc2d, dcapecalc3d, dcloudfrac, + wrfcttcalc, calcdbz, dcalrelhl, dcalcuh, dcomputepv, dcomputeabsvort, dlltoij, dijtoll, deqthecalc, omgcalc, virtual_temp, wetbulbcalc, dcomputepw, wrf_monotonic, wrf_vintrp, dcomputewspd, @@ -597,8 +597,13 @@ def _cape(p_hpa, tk, qv, ht, ter, sfp, missing, i3dflag, ter_follow, errstat = np.array(0) errmsg = np.zeros(Constants.ERRLEN, "c") + if i3dflag: + cape_routine = dcapecalc3d + else: + cape_routine = dcapecalc2d + # note that p_hpa, tk, qv, and ht have the vertical flipped - result = dcapecalc3d(p_hpa, + result = cape_routine(p_hpa, tk, qv, ht, @@ -607,7 +612,6 @@ def _cape(p_hpa, tk, qv, ht, ter, sfp, missing, i3dflag, ter_follow, capeview, cinview, missing, - i3dflag, ter_follow, psafile, errstat,