From 71fe678e916ab632ad8fb2b99ddc19a34bf7aa9d Mon Sep 17 00:00:00 2001 From: Bill Ladwig Date: Fri, 9 Mar 2018 11:42:19 -0700 Subject: [PATCH] Moved work arrays out of fortran. Fixes #47 --- fortran/rip_cape.f90 | 28 ++++++++++++++++--------- fortran/wrf_fctt.f90 | 6 ++++-- fortran/wrf_vinterp.f90 | 6 ++++-- src/wrf/extension.py | 45 +++++++++++++++++++++++++++++------------ 4 files changed, 58 insertions(+), 27 deletions(-) diff --git a/fortran/rip_cape.f90 b/fortran/rip_cape.f90 index 5e9cdac..1fc3b23 100644 --- a/fortran/rip_cape.f90 +++ b/fortran/rip_cape.f90 @@ -272,6 +272,7 @@ END SUBROUTINE DPFCALC ! NCLFORTSTART SUBROUTINE DCAPECALC3D(prs,tmk,qvp,ght,ter,sfp,cape,cin,& + prsf, prs_new, tmk_new, qvp_new, ght_new,& cmsg,mix,mjy,mkzh,ter_follow,& psafile, errstat, errmsg) USE wrf_constants, ONLY : CELKEL, G, EZERO, ESLCON1, ESLCON2, & @@ -293,6 +294,14 @@ SUBROUTINE DCAPECALC3D(prs,tmk,qvp,ght,ter,sfp,cape,cin,& REAL(KIND=8), DIMENSION(mix,mjy), INTENT(IN) ::sfp REAL(KIND=8), DIMENSION(mix,mjy,mkzh), INTENT(OUT) :: cape REAL(KIND=8), DIMENSION(mix,mjy,mkzh), INTENT(OUT) :: cin + + REAL(KIND=8), DIMENSION(mkzh,mix,mjy), INTENT(INOUT) :: prsf + REAL(KIND=8), DIMENSION(mkzh,mix,mjy), INTENT(INOUT) :: prs_new + REAL(KIND=8), DIMENSION(mkzh,mix,mjy), INTENT(INOUT) :: tmk_new + REAL(KIND=8), DIMENSION(mkzh,mix,mjy), INTENT(INOUT) :: qvp_new + REAL(KIND=8), DIMENSION(mkzh,mix,mjy), INTENT(INOUT) :: ght_new + + REAL(KIND=8), INTENT(IN) :: cmsg CHARACTER(LEN=*), INTENT(IN) :: psafile INTEGER, INTENT(INOUT) :: errstat @@ -308,15 +317,10 @@ SUBROUTINE DCAPECALC3D(prs,tmk,qvp,ght,ter,sfp,cape,cin,& 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,mix,mjy) :: prsf REAL(KIND=8), DIMENSION(150) :: psadithte, psadiprs REAL(KIND=8), DIMENSION(150,150) :: psaditmk LOGICAL :: elfound - REAL(KIND=8), DIMENSION(mkzh,mix,mjy) :: prs_new - REAL(KIND=8), DIMENSION(mkzh,mix,mjy) :: tmk_new - REAL(KIND=8), DIMENSION(mkzh,mix,mjy) :: qvp_new - REAL(KIND=8), DIMENSION(mkzh,mix,mjy) :: ght_new ! To remove compiler warnings tmkpari = 0 @@ -597,6 +601,7 @@ END SUBROUTINE DCAPECALC3D ! NCLFORTSTART SUBROUTINE DCAPECALC2D(prs,tmk,qvp,ght,ter,sfp,cape,cin,& + prsf, prs_new, tmk_new, qvp_new, ght_new,& cmsg,mix,mjy,mkzh,ter_follow,& psafile, errstat, errmsg) USE wrf_constants, ONLY : CELKEL, G, EZERO, ESLCON1, ESLCON2, & @@ -618,6 +623,13 @@ SUBROUTINE DCAPECALC2D(prs,tmk,qvp,ght,ter,sfp,cape,cin,& REAL(KIND=8), DIMENSION(mix,mjy), INTENT(IN) ::sfp REAL(KIND=8), DIMENSION(mix,mjy,mkzh), INTENT(OUT) :: cape REAL(KIND=8), DIMENSION(mix,mjy,mkzh), INTENT(OUT) :: cin + + REAL(KIND=8), DIMENSION(mkzh,mix,mjy), INTENT(INOUT) :: prsf + REAL(KIND=8), DIMENSION(mkzh,mix,mjy), INTENT(INOUT) :: prs_new + REAL(KIND=8), DIMENSION(mkzh,mix,mjy), INTENT(INOUT) :: tmk_new + REAL(KIND=8), DIMENSION(mkzh,mix,mjy), INTENT(INOUT) :: qvp_new + REAL(KIND=8), DIMENSION(mkzh,mix,mjy), INTENT(INOUT) :: ght_new + REAL(KIND=8), INTENT(IN) :: cmsg CHARACTER(LEN=*), INTENT(IN) :: psafile INTEGER, INTENT(INOUT) :: errstat @@ -635,15 +647,11 @@ SUBROUTINE DCAPECALC2D(prs,tmk,qvp,ght,ter,sfp,cape,cin,& 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(mkzh,mix,mjy) :: prsf REAL(KIND=8), DIMENSION(150) :: psadithte, psadiprs REAL(KIND=8), DIMENSION(150,150) :: psaditmk LOGICAL :: elfound REAL(KIND=8), DIMENSION(mkzh) :: eth_temp - REAL(KIND=8), DIMENSION(mkzh,mix,mjy) :: prs_new - REAL(KIND=8), DIMENSION(mkzh,mix,mjy) :: tmk_new - REAL(KIND=8), DIMENSION(mkzh,mix,mjy) :: qvp_new - REAL(KIND=8), DIMENSION(mkzh,mix,mjy) :: ght_new + ! To remove compiler warnings errstat = 0 diff --git a/fortran/wrf_fctt.f90 b/fortran/wrf_fctt.f90 index 2ef7ecb..0be18ae 100644 --- a/fortran/wrf_fctt.f90 +++ b/fortran/wrf_fctt.f90 @@ -1,5 +1,5 @@ !NCLFORTSTART -SUBROUTINE wrfcttcalc(prs, tk, qci, qcw, qvp, ght, ter, ctt, haveqci, nz, ns, ew) +SUBROUTINE wrfcttcalc(prs, tk, qci, qcw, qvp, ght, ter, ctt, pf, haveqci, nz, ns, ew) USE wrf_constants, ONLY : EPS, USSALR, RD, G, ABSCOEFI, ABSCOEF, CELKEL IMPLICIT NONE @@ -12,6 +12,8 @@ SUBROUTINE wrfcttcalc(prs, tk, qci, qcw, qvp, ght, ter, ctt, haveqci, nz, ns, ew REAL(KIND=8), DIMENSION(ew,ns), INTENT(IN) :: ter REAL(KIND=8), DIMENSION(ew,ns), INTENT(OUT) :: ctt + REAL(KIND=8), DIMENSION(ew,ns,nz), INTENT(INOUT) :: pf + !NCLEND ! REAL(KIND=8) :: znfac(nz) @@ -20,7 +22,7 @@ SUBROUTINE wrfcttcalc(prs, tk, qci, qcw, qvp, ght, ter, ctt, haveqci, nz, ns, ew INTEGER i,j,k,ripk REAL(KIND=8) :: opdepthu, opdepthd, dp, arg1, fac, prsctt, ratmix REAL(KIND=8) :: arg2, agl_hgt, vt - REAL(KIND=8), DIMENSION(ew,ns,nz) :: pf + REAL(KIND=8) :: p1, p2 !$OMP PARALLEL diff --git a/fortran/wrf_vinterp.f90 b/fortran/wrf_vinterp.f90 index c34efd7..d81270b 100644 --- a/fortran/wrf_vinterp.f90 +++ b/fortran/wrf_vinterp.f90 @@ -127,7 +127,7 @@ END FUNCTION wrf_intrp_value !NCLFORTSTART SUBROUTINE wrf_vintrp(datain, dataout, pres, tk, qvp, ght, terrain,& sfp, smsfp, vcarray, interp_levels, numlevels,& - icase, ew, ns, nz, extrap, vcor, logp, rmsg,& + icase, ew, ns, nz, extrap, vcor, logp, tempout, rmsg,& errstat, errmsg) USE wrf_constants, ONLY : ALGERR, SCLHT, EXPON, EXPONI, GAMMA, GAMMAMD, TLCLC1, & TLCLC2, TLCLC3, TLCLC4, THTECON1, THTECON2, THTECON3, & @@ -146,6 +146,9 @@ SUBROUTINE wrf_vintrp(datain, dataout, pres, tk, qvp, ght, terrain,& REAL(KIND=8), DIMENSION(ew,ns,numlevels), INTENT(OUT) :: dataout REAL(KIND=8), DIMENSION(ew,ns,nz), INTENT(IN) :: vcarray REAL(KIND=8), DIMENSION(numlevels), INTENT(IN) :: interp_levels + + REAL(KIND=8), DIMENSION(ew,ns), INTENT(INOUT) :: tempout + REAL(KIND=8), INTENT(IN) :: rmsg INTEGER, INTENT(INOUT) :: errstat CHARACTER(LEN=*), INTENT(INOUT) :: errmsg @@ -156,7 +159,6 @@ SUBROUTINE wrf_vintrp(datain, dataout, pres, tk, qvp, ght, terrain,& INTEGER :: i, j, k, kupper !itriv INTEGER :: ifound, isign !miy,mjx INTEGER :: log_errcnt, interp_errcnt, interp_errstat - REAL(KIND=8), DIMENSION(ew,ns) :: tempout REAL(KIND=8) :: rlevel, vlev, diff REAL(KIND=8) :: tmpvlev REAL(KIND=8) :: vcp1, vcp0, valp0, valp1 diff --git a/src/wrf/extension.py b/src/wrf/extension.py index 167d27c..5bcccfa 100755 --- a/src/wrf/extension.py +++ b/src/wrf/extension.py @@ -618,20 +618,33 @@ def _cape(p_hpa, tk, qv, ht, ter, sfp, missing, i3dflag, ter_follow, else: cape_routine = dcapecalc2d + # Work arrays + k_left_shape = (p_hpa.shape[2], p_hpa.shape[0], p_hpa.shape[1]) + prsf = np.empty(k_left_shape, np.float64, order="F") + prs_new = np.empty(k_left_shape, np.float64, order="F") + tmk_new = np.empty(k_left_shape, np.float64, order="F") + qvp_new = np.empty(k_left_shape, np.float64, order="F") + ght_new = np.empty(k_left_shape, np.float64, order="F") + # note that p_hpa, tk, qv, and ht have the vertical flipped result = cape_routine(p_hpa, - tk, - qv, - ht, - ter, - sfp, - capeview, - cinview, - missing, - ter_follow, - psafile, - errstat, - errmsg) + tk, + qv, + ht, + ter, + sfp, + capeview, + cinview, + prsf, + prs_new, + tmk_new, + qvp_new, + ght_new, + missing, + ter_follow, + psafile, + errstat, + errmsg) if int(errstat) != 0: raise DiagnosticError("".join(npbytes_to_str(errmsg)).strip()) @@ -766,6 +779,8 @@ def _ctt(p_hpa, tk, qice, qcld, qv, ght, ter, haveqci, outview=None): if outview is None: outview = np.empty_like(ter) + pf = np.empty(p_hpa.shape[0:3], np.float64, order="F") + result = wrfcttcalc(p_hpa, tk, qice, @@ -774,6 +789,7 @@ def _ctt(p_hpa, tk, qice, qcld, qv, ght, ter, haveqci, outview=None): ght, ter, outview, + pf, haveqci) return result @@ -858,7 +874,9 @@ def _vintrp(field, pres, tk, qvp, ght, terrain, sfp, smsfp, if outview is None: outdims = field.shape[0:2] + interp_levels.shape outview = np.empty(outdims, field.dtype, order="F") - + + tempout = np.zeros(field.shape[0:2], np.float64, order="F") + errstat = np.array(0) errmsg = np.zeros(Constants.ERRLEN, "c") @@ -877,6 +895,7 @@ def _vintrp(field, pres, tk, qvp, ght, terrain, sfp, smsfp, extrap, vcor, logp, + tempout, missing, errstat, errmsg)