nilyin
/
wrf-python
forked from 3rdparty/wrf-python
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity
Browse Source

Added changes for faster cape and twb routines done by @supreethms1809

lon0
Bill Ladwig 8 years ago
parent
0be864a194
commit
74ded8f0f9
4 changed files with 620 additions and 197 deletions
Whitespace
Show all changes Ignore whitespace when comparing lines Ignore changes in amount of whitespace Ignore changes in whitespace at EOL
Unified View
Diff Options
Show Stats Download Patch File Download Diff File
  1. 704
      fortran/rip_cape.f90
  2. 64
      fortran/wrf_rip_phys_routines.f90
  3. 31
      fortran/wrffortran.pyf
  4. 12
      src/wrf/extension.py

704
fortran/rip_cape.f90
Unescape View File

@ -37,7 +37,8 @@ END FUNCTION TVIRTUAL
REAL(KIND=8) FUNCTION TONPSADIABAT(thte, prs, psadithte, psadiprs, psaditmk, gamma,& REAL(KIND=8) FUNCTION TONPSADIABAT(thte, prs, psadithte, psadiprs, psaditmk, gamma,&
errstat, errmsg) errstat, errmsg)
USE wrf_constants, ONLY : ALGERR USE wrf_constants, ONLY : ALGERR
!$OMP DECLARE SIMD (TONPSADIABAT)
!!uniform(thte,prs,psadithte,psadiprs,psaditmk)
!f2py threadsafe !f2py threadsafe
!f2py intent(in,out) :: cape, cin !f2py intent(in,out) :: cape, cin
@ -58,7 +59,8 @@ REAL(KIND=8) FUNCTION TONPSADIABAT(thte, prs, psadithte, psadiprs, psaditmk, gam
REAL(KIND=8) :: fracip REAL(KIND=8) :: fracip
REAL(KIND=8) :: fracip2 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 ! This function gives the temperature (in K) on a moist adiabat
! (specified by thte in K) given pressure in hPa. It uses a ! (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. ! Otherwise, look for the given thte/prs point in the lookup table.
jt = -1 jt = -1
DO jtch = 1, 150-1 l1 = 1
IF (thte .GE. psadithte(jtch) .AND. thte .LT. psadithte(jtch+1)) THEN h1 = 149
jt = jtch rang1 = h1 - l1
EXIT mid1 = (h1 + l1) / 2
!GO TO 213 DO WHILE(rang1 .GT. 1)
END IF if(thte .GE. psadithte(mid1)) then
l1 = mid1
else
h1 = mid1
end if
rang1 = h1 - l1
mid1 = (h1 + l1) / 2
END DO END DO
jt = l1
ip = -1
DO ipch = 1, 150-1 ! DO jtch = 1, 150-1
IF (prs .LE. psadiprs(ipch) .AND. prs .GT. psadiprs(ipch+1)) THEN ! IF (thte .GE. psadithte(jtch) .AND. thte .LT. psadithte(jtch+1)) THEN
ip = ipch ! jt = jtch
EXIT ! EXIT
!GO TO 215 ! !GO TO 213
END IF ! 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 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 IF (jt .EQ. -1 .OR. ip .EQ. -1) THEN
! Set the error and return ! Set the error and return
@ -187,28 +220,26 @@ END SUBROUTINE DLOOKUP_TABLE
! level is above. ! level is above.
SUBROUTINE DPFCALC(prs, sfp, pf, miy, mjx, mkzh, ter_follow) 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), 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, INTENT(IN) :: ter_follow,miy,mjx,mkzh
INTEGER :: i,j,k INTEGER :: i,j,k
! do j=1,mjx-1 Artifact of MM5
DO j = 1,mjx DO j = 1,mjx
! do i=1,miy-1 staggered grid
DO i = 1,miy DO i = 1,miy
DO k = 1,mkzh DO k = 1,mkzh
IF (k .EQ. mkzh) THEN IF (k .EQ. mkzh) THEN
! terrain-following data ! terrain-following data
IF (ter_follow .EQ. 1) THEN IF (ter_follow .EQ. 1) THEN
pf(i,j,k) = sfp(i,j) pf(k,i,j) = sfp(i,j)
! pressure-level data ! pressure-level data
ELSE 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 END IF
ELSE 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 IF
END DO END DO
END DO END DO
@ -224,17 +255,338 @@ END SUBROUTINE DPFCALC
! !
! !DESCRIPTION: ! !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 ! or J/kg) for every grid point in the entire 3D domain (treating
! each grid point as a parcel). If i3dflag=0, then it ! each grid point as a parcel).
! calculates CAPE and CIN only for the parcel with max theta-e in !
! 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", ! the column, (i.e. something akin to Colman's MCAPE). By "parcel",
! we mean a 500-m deep parcel, with actual temperature and moisture ! we mean a 500-m deep parcel, with actual temperature and moisture
! averaged over that depth. ! 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 ! 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. ! 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 ! 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). ! 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 ! Also also, Pressure must be hPa
! NCLFORTSTART ! NCLFORTSTART
SUBROUTINE DCAPECALC3D(prs,tmk,qvp,ght,ter,sfp,cape,cin,& SUBROUTINE DCAPECALC2D(prs,tmk,qvp,ght,ter,sfp,cape,cin,&
cmsg,miy,mjx,mkzh,i3dflag,ter_follow,& cmsg,miy,mjx,mkzh,ter_follow,&
psafile, errstat, errmsg) psafile, errstat, errmsg)
USE wrf_constants, ONLY : ALGERR, CELKEL, G, EZERO, ESLCON1, ESLCON2, & USE wrf_constants, ONLY : ALGERR, CELKEL, G, EZERO, ESLCON1, ESLCON2, &
EPS, RD, CP, GAMMA, CPMD, RGASMD, GAMMAMD, TLCLC1, & EPS, RD, CP, GAMMA, CPMD, RGASMD, GAMMAMD, TLCLC1, &
TLCLC2, TLCLC3, TLCLC4, THTECON1, THTECON2, THTECON3 TLCLC2, TLCLC3, TLCLC4, THTECON1, THTECON2, THTECON3
!USE omp_lib
IMPLICIT NONE IMPLICIT NONE
!f2py threadsafe !f2py threadsafe
!f2py intent(in,out) :: cape, cin !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) :: prs
REAL(KIND=8), DIMENSION(miy,mjx,mkzh), INTENT(IN) :: tmk 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) :: qvp
@ -275,26 +629,35 @@ SUBROUTINE DCAPECALC3D(prs,tmk,qvp,ght,ter,sfp,cape,cin,&
! NCLFORTEND ! NCLFORTEND
! local variables ! local variables
INTEGER :: i, j, k, ilcl, kel, kk, klcl, klev, klfc, kmax, kpar, kpar1, kpar2 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) :: pavg, tvirtual, p1, p2, pp1, pp2, th, totthe, totqvp, totprs
REAL(KIND=8) :: cpm, deltap, ethpari, gammam, ghtpari, qvppari, prspari, tmkpari REAL(KIND=8) :: cpm, deltap, ethpari, gammam, qvppari, tmkpari
REAL(KIND=8) :: facden, fac1, fac2, qvplift, tmklift, tvenv, tvlift, ghtlift REAL(KIND=8) :: facden, qvplift, tmklift, tvenv, tvlift, ghtlift, fac1, fac2
REAL(KIND=8) :: eslift, tmkenv, qvpenv, tonpsadiabat REAL(KIND=8) :: eslift, tmkenv, qvpenv, tonpsadiabat
REAL(KIND=8) :: benamin, dz, pup, pdn REAL(KIND=8) :: benamin, dz, pup, pdn
REAL(KIND=8), DIMENSION(150) :: buoy, zrel, benaccum 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) :: psadithte, psadiprs
REAL(KIND=8), DIMENSION(150,150) :: psaditmk REAL(KIND=8), DIMENSION(150,150) :: psaditmk
LOGICAL :: elfound 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 ! To remove compiler warnings
errstat = 0
tmkpari = 0 tmkpari = 0
qvppari = 0 qvppari = 0
klev = 0 klev = 0
klcl = 0 klcl = 0
kel = 0 kel = 0
deltap = 0
! the comments were taken from a mark stoelinga email, 23 apr 2007, ! 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) ! 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 ! calculated the pressure at full sigma levels (a set of pressure
! levels that bound the layers represented by the vertical grid points) ! levels that bound the layers represented by the vertical grid points)
CALL DPFCALC(prs_new, sfp, prsf, miy, mjx, mkzh, ter_follow)
CALL DPFCALC(prs, sfp, prsf, miy, mjx, mkzh, ter_follow)
! before looping, set lookup table for getting temperature on ! before looping, set lookup table for getting temperature on
! a pseudoadiabat. ! a pseudoadiabat.
@ -334,148 +710,144 @@ SUBROUTINE DCAPECALC3D(prs,tmk,qvp,ght,ter,sfp,cape,cin,&
RETURN RETURN
END IF 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 j = 1,mjx
! do i=1,miy-1
DO i = 1,miy DO i = 1,miy
cape(i,j,1) = 0.d0 cape(i,j,1) = 0.d0
cin(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 = klev
kpar1 = 2 kpar2 = klev
kpar2 = mkzh
ELSE
! find parcel with max theta-e in lowest 3 km agl. ! Establish average properties of that parcel
ethmax = -1.d0 ! (over depth of approximately davg meters)
DO k = mkzh,1,-1
IF (ght(i,j,k)-ter(i,j) .LT. 3000.d0) THEN !davg = 500.d0
q = MAX(qvp(i,j,k), 1.d-15) pavg = 500.d0 * prs_new(kpar1,i,j)*&
t = tmk(i,j,k) G/(RD*tvirtual(tmk_new(kpar1,i,j), qvp_new(kpar1,i,j)))
p = prs(i,j,k) p2 = MIN(prs_new(kpar1,i,j)+.5d0*pavg, prsf(mkzh,i,j))
e = q*p/(EPS + q) p1 = p2 - pavg
tlcl = TLCLC1 / (LOG(t**TLCLC2/e)-TLCLC3) + TLCLC4 totthe = 0.D0
eth = t * (1000.d0/p)**(GAMMA*(1.d0 + GAMMAMD*q))*& totqvp = 0.D0
EXP((THTECON1/tlcl - THTECON2)*q*(1.d0 + THTECON3*q)) totprs = 0.D0
IF (eth .GT. ethmax) THEN DO k = mkzh,2,-1
klev = k IF (prsf(k,i,j) .LE. p1) EXIT !GOTO 35
ethmax = eth IF (prsf(k-1,i,j) .GE. p2) CYCLE !GOTO 34
END IF p = prs_new(k,i,j)
END IF pup = prsf(k,i,j)
END DO pdn = prsf(k-1,i,j)
kpar1 = klev !q = MAX(qvp_new(k,i,j),1.D-15)
kpar2 = klev 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)
! Establish average properties of that parcel pp2 = MIN(p2,pup)
! (over depth of approximately davg meters) IF (pp2 .GT. pp1) THEN
! deltap = pp2 - pp1
! davg=.1 totqvp = totqvp + MAX(qvp_new(k,i,j),1.D-15)*(pp2 - pp1)
davg = 500.d0 totthe = totthe + th*(pp2 - pp1)
pavg = davg*prs(i,j,kpar1)*& totprs = totprs + (pp2 - pp1)
G/(RD*tvirtual(tmk(i,j,kpar1), qvp(i,j,kpar1))) END IF
p2 = MIN(prs(i,j,kpar1)+.5d0*pavg, prsf(i,j,mkzh)) END DO
p1 = p2 - pavg qvppari = totqvp/totprs
totthe = 0.D0 tmkpari = (totthe/totprs)*&
totqvp = 0.D0 (prs_new(kpar1,i,j)/1000.D0)**(GAMMA*(1.D0+GAMMAMD*qvp_new(kpar1,i,j)))
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
!CALL CPU_TIME(t3)
DO kpar = kpar1, kpar2 DO kpar = kpar1, kpar2
! Calculate temperature and moisture properties of parcel ! 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 !prspari = prs_new(kpar,i,j)
qvppari = qvp(i,j,kpar) !ghtpari = ght_new(kpar,i,j)
tmkpari = tmk(i,j,kpar)
END IF
prspari = prs(i,j,kpar)
ghtpari = ght(i,j,kpar)
gammam = GAMMA * (1.D0 + GAMMAMD*qvppari) gammam = GAMMA * (1.D0 + GAMMAMD*qvppari)
cpm = CP * (1.D0 + CPMD*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 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)) 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 ! 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. ! and at all points where buoyancy is zero.
! !
!
! For arrays that go bottom to top ! For arrays that go bottom to top
kk = 0 kk = 0
ilcl = 0 ilcl = 0
IF (ghtpari .GE. zlcl) THEN IF (ght_new(kpar,i,j) .GE. zlcl) THEN
! Initial parcel already saturated or supersaturated. ! Initial parcel already saturated or supersaturated.
ilcl = 2 ilcl = 2
klcl = 1 klcl = 1
END IF END IF
k = kpar k = kpar
DO WHILE (k .GE. 1)!k = kpar, 1, -1 DO k = kpar,1,-1
!DO k = kpar, 1, -1 ! For arrays that go bottom to top
! For arrays that go bottom to top
! 33 kk = kk + 1
kk = kk + 1 kk = kk + 1
! Model level is below lcl ! Model level is below lcl
IF (ght(i,j,k) .LT. zlcl) THEN IF (ght_new(k,i,j) .LT. zlcl) THEN
qvplift = qvppari 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))
tmklift = tmkpari - G/cpm*(ght(i,j,k) - ghtpari) tvenv = tmk_new(k,i,j)*(EPS + qvp_new(k,i,j))/(EPS*(1.D0 + qvp_new(k,i,j)))
tvenv = tvirtual(tmk(i,j,k), qvp(i,j,k)) tvlift = tmklift*(EPS + qvp_new(kpar,i,j))/(EPS*(1.D0 + qvp_new(kpar,i,j)))
tvlift = tvirtual(tmklift, qvplift) ghtlift = ght_new(k,i,j)
ghtlift = ght(i,j,k)
ELSE IF (ght(i,j,k) .GE. zlcl .AND. ilcl .EQ. 0) THEN ELSE IF (ght(i,j,k) .GE. zlcl .AND. ilcl .EQ. 0) THEN
! This model level and previous model level straddle the lcl, ! This model level and previous model level straddle the lcl,
! so first create a new level in the bottom-up array, at the lcl. ! so first create a new level in the bottom-up array, at the lcl.
tmklift = tlcl facden = 1/(ght_new(k,i,j) - ght_new(k+1,i,j))
qvplift = qvppari 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)
facden = ght(i,j,k) - ght(i,j,k+1) 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)
fac1 = (zlcl-ght(i,j,k+1))/facden tvenv = tmkenv* (EPS + qvpenv) / (EPS * (1.D0 + qvpenv))
fac2 = (ght(i,j,k)-zlcl)/facden tvlift = tlcl* (EPS + qvp_new(kpar,i,j)) / (EPS *(1.D0 + qvp_new(kpar,i,j)))
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)
ghtlift = zlcl ghtlift = zlcl
ilcl = 1 ilcl = 1
ELSE ELSE
tmklift = TONPSADIABAT(ethpari, prs(i,j,k), psadithte, psadiprs,& tmklift = TONPSADIABAT(ethpari, prs_new(k,i,j), psadithte, psadiprs,&
psaditmk, GAMMA, errstat, errmsg) psaditmk, GAMMA, errstat, errmsg)
eslift = EZERO*EXP(ESLCON1*(tmklift - CELKEL)/(tmklift - ESLCON2)) eslift = EZERO*EXP(ESLCON1*(tmklift - CELKEL)/(tmklift - ESLCON2))
qvplift = EPS*eslift/(prs(i,j,k) - eslift) qvplift = EPS*eslift/(prs_new(k,i,j) - eslift)
tvenv = tvirtual(tmk(i,j,k), qvp(i,j,k)) tvenv = tmk_new(k,i,j) * (EPS + qvp_new(k,i,j)) / (EPS * (1.D0 + qvp_new(k,i,j)))
tvlift = tvirtual(tmklift, qvplift) tvlift = tmklift*(EPS + qvplift) / (EPS * (1.D0 + qvplift))
ghtlift = ght(i,j,k) ghtlift = ght_new(k,i,j)
END IF END IF
! Buoyancy ! Buoyancy
buoy(kk) = G*(tvlift - tvenv)/tvenv 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 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 ! Parcel ascent curve crosses sounding curve, so create a new level
! in the bottom-up array at the crossing. ! 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)/& zrel(kk-1) = zrel(kk-2) + buoy(kk-2)/&
(buoy(kk-2) - buoy(kk))*(zrel(kk) - zrel(kk-2)) (buoy(kk-2) - buoy(kk))*(zrel(kk) - zrel(kk-2))
END IF END IF
IF (ilcl .EQ. 1) THEN IF (ilcl .EQ. 1) THEN
klcl = kk klcl = kk
ilcl = 2 ilcl = 2
!GOTO 33
CYCLE CYCLE
END IF END IF
k = k - 1
END DO END DO
kmax = kk kmax = kk
IF (kmax .GT. 150) THEN ! IF (kmax .GT. 150) THEN
errstat = ALGERR ! errstat = ALGERR
WRITE(errmsg, *) 'capecalc3d: kmax got too big. kmax=',kmax ! WRITE(errmsg, *) 'capecalc3d: kmax got too big. kmax=',kmax
RETURN ! RETURN
END IF ! END IF
! If no lcl was found, set klcl to kmax. it is probably not really ! 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 ! at kmax, but this will make the rest of the routine behave
! properly. ! properly.
IF (ilcl .EQ. 0) klcl=kmax IF (ilcl .EQ. 0) klcl=kmax
@ -520,7 +890,6 @@ SUBROUTINE DCAPECALC3D(prs,tmk,qvp,ght,ter,sfp,cape,cin,&
benamin = benaccum(k) benamin = benaccum(k)
END IF END IF
END DO END DO
! Determine equilibrium level (el), which we define as the highest ! Determine equilibrium level (el), which we define as the highest
! level of non-negative buoyancy above the lcl. note, this may be ! level of non-negative buoyancy above the lcl. note, this may be
! the top level if the parcel is still buoyant there. ! 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 kel = k
elfound = .TRUE. elfound = .TRUE.
EXIT EXIT
!GOTO 50
END IF END IF
END DO 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 ! 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 ! that the zero contour in either the cin or cape fields will
! circumscribe regions of non-zero cape. ! 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 ! value of this function. at that time we decided to change -0.1
! to a more appropriate missing value, which is passed into this ! to a more appropriate missing value, which is passed into this
! routine as cmsg. ! routine as cmsg.
! cape(i,j,kpar) = -0.1D0
! cin(i,j,kpar) = -0.1D0
IF (.NOT. elfound) THEN IF (.NOT. elfound) THEN
cape(i,j,kpar) = cmsg cape(i,j,kpar) = cmsg
cin(i,j,kpar) = cmsg cin(i,j,kpar) = cmsg
@ -558,17 +923,20 @@ SUBROUTINE DCAPECALC3D(prs,tmk,qvp,ght,ter,sfp,cape,cin,&
CYCLE CYCLE
END IF END IF
! GOTO 102
! 50 CONTINUE ! If there is an equilibrium level, then cape is positive.
! we'll
! If there is an equilibrium level, then cape is positive. we'll ! define the level of free convection (lfc) as the point below
! define the level of free convection (lfc) as the point below the ! the
! el, but at or above the lcl, where accumulated buoyant energy is a ! el, but at or above the lcl, where accumulated buoyant energy
! minimum. the net positive area (accumulated buoyant energy) from ! 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 ! the lfc up to the el will be defined as the cape, and the net
! negative area (negative of accumulated buoyant energy) from the ! negative area (negative of accumulated buoyant energy) from
! parcel starting point to the lfc will be defined as the convective ! the
! parcel starting point to the lfc will be defined as the
! convective
! inhibition (cin). ! inhibition (cin).
! First get the lfc according to the above definition. ! 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 ! to a more appropriate missing value, which is passed into this
! routine as cmsg. ! 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 IF (cape(i,j,kpar) .LT. 100.D0) cin(i,j,kpar) = cmsg
! 102 CONTINUE
END DO END DO
IF (i3dflag .EQ. 0) THEN
cape(i,j,mkzh) = cape(i,j,kpar1) cape(i,j,mkzh) = cape(i,j,kpar1)
cin(i,j,mkzh) = cin(i,j,kpar1) cin(i,j,mkzh) = cin(i,j,kpar1)
! meters agl ! 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 ! 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
END DO END DO
!$OMP END PARALLEL DO
RETURN RETURN
END SUBROUTINE DCAPECALC3D END SUBROUTINE DCAPECALC2D

64
fortran/wrf_rip_phys_routines.f90
Unescape View File

@ -50,13 +50,15 @@ SUBROUTINE WETBULBCALC(prs, tmk, qvp, twb, nx, ny, nz, psafile, errstat, errmsg)
!NCLEND !NCLEND
INTEGER :: i, j, k INTEGER :: i, j, k
INTEGER :: jtch, jt, ipch, ip INTEGER :: jt, ip
REAL(KIND=8) :: q, t, p, e, tlcl, eth REAL(KIND=8) :: q, t, p, e, tlcl, eth
REAL(KIND=8) :: fracip, fracip2, fracjt, fracjt2 REAL(KIND=8) :: fracip, fracip2, fracjt, fracjt2
REAL(KIND=8), DIMENSION(150) :: PSADITHTE, PSADIPRS REAL(KIND=8), DIMENSION(150) :: PSADITHTE, PSADIPRS
REAL(KIND=8), DIMENSION(150,150) :: PSADITMK REAL(KIND=8), DIMENSION(150,150) :: PSADITMK
REAL(KIND=8) :: tonpsadiabat 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 ! Before looping, set lookup table for getting temperature on
! a pseudoadiabat. ! a pseudoadiabat.
@ -92,21 +94,53 @@ SUBROUTINE WETBULBCALC(prs, tmk, qvp, twb, nx, ny, nz, psafile, errstat, errmsg)
tonpsadiabat = eth*(p/1000.)**GAMMA tonpsadiabat = eth*(p/1000.)**GAMMA
ELSE ELSE
! Otherwise, look for the given thte/prs point in the lookup table. ! Otherwise, look for the given thte/prs point in the lookup table.
jt=-1 jt = -1
DO jtch=1,150-1 l1 = 1
IF (eth .GE. PSADITHTE(jtch) .AND. eth .LT. PSADITHTE(jtch+1)) THEN h1 = 149
jt = jtch rang1 = h1 - l1
EXIT mid1 = (h1 + l1) / 2
ENDIF DO WHILE(rang1 .GT. 1)
END DO if(eth .GE. psadithte(mid1)) then
l1 = mid1
ip=-1 else
DO ipch=1,150-1 h1 = mid1
IF (p .LE. PSADIPRS(ipch) .AND. p .GT. PSADIPRS(ipch+1)) THEN end if
ip = ipch rang1 = h1 - l1
EXIT mid1 = (h1 + l1) / 2
ENDIF
END DO 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 IF (jt .EQ. -1 .OR. ip .EQ. -1) THEN
errstat = ALGERR errstat = ALGERR

31
fortran/wrffortran.pyf
Unescape View File

@ -64,16 +64,38 @@ python module _wrffortran ! in
character*(*) intent(inout) :: errmsg character*(*) intent(inout) :: errmsg
end subroutine dlookup_table end subroutine dlookup_table
subroutine dpfcalc(prs,sfp,pf,miy,mjx,mkzh,ter_follow) ! in :_wrffortran:rip_cape.f90 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) :: 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),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,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,1)==mjx),depend(prs) :: mjx=shape(prs,1)
integer, optional,intent(in),check(shape(prs,2)==mkzh),depend(prs) :: mkzh=shape(prs,2) integer, optional,intent(in),check(shape(prs,2)==mkzh),depend(prs) :: mkzh=shape(prs,2)
integer intent(in) :: ter_follow integer intent(in) :: ter_follow
end subroutine dpfcalc character*(*) intent(in) :: psafile
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 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 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 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) :: 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) :: 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,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,1)==mjx),depend(prs) :: mjx=shape(prs,1)
integer, optional,intent(in),check(shape(prs,2)==mkzh),depend(prs) :: mkzh=shape(prs,2) integer, optional,intent(in),check(shape(prs,2)==mkzh),depend(prs) :: mkzh=shape(prs,2)
integer intent(in) :: i3dflag
integer intent(in) :: ter_follow integer intent(in) :: ter_follow
character*(*) intent(in) :: psafile character*(*) intent(in) :: psafile
integer intent(inout) :: errstat integer intent(inout) :: errstat
character*(*) intent(inout) :: errmsg 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 subroutine dcloudfrac(pres,rh,lowc,midc,highc,nz,ns,ew) ! in :_wrffortran:wrf_cloud_fracf.f90
threadsafe threadsafe
real(kind=8) dimension(ew,ns,nz),intent(in) :: pres real(kind=8) dimension(ew,ns,nz),intent(in) :: pres

12
src/wrf/extension.py
Unescape View File

@ -7,8 +7,8 @@ from .constants import Constants
from ._wrffortran import (dcomputetk, dinterp3dz, dinterp2dxy, dinterp1d, from ._wrffortran import (dcomputetk, dinterp3dz, dinterp2dxy, dinterp1d,
dcomputeseaprs, dfilter2d, dcomputerh, dcomputeuvmet, dcomputeseaprs, dfilter2d, dcomputerh, dcomputeuvmet,
dcomputetd, dcapecalc3d, dcloudfrac, wrfcttcalc, dcomputetd, dcapecalc2d, dcapecalc3d, dcloudfrac,
calcdbz, dcalrelhl, dcalcuh, dcomputepv, wrfcttcalc, calcdbz, dcalrelhl, dcalcuh, dcomputepv,
dcomputeabsvort, dlltoij, dijtoll, deqthecalc, dcomputeabsvort, dlltoij, dijtoll, deqthecalc,
omgcalc, virtual_temp, wetbulbcalc, dcomputepw, omgcalc, virtual_temp, wetbulbcalc, dcomputepw,
wrf_monotonic, wrf_vintrp, dcomputewspd, 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) errstat = np.array(0)
errmsg = np.zeros(Constants.ERRLEN, "c") 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 # note that p_hpa, tk, qv, and ht have the vertical flipped
result = dcapecalc3d(p_hpa, result = cape_routine(p_hpa,
tk, tk,
qv, qv,
ht, ht,
@ -607,7 +612,6 @@ def _cape(p_hpa, tk, qv, ht, ter, sfp, missing, i3dflag, ter_follow,
capeview, capeview,
cinview, cinview,
missing, missing,
i3dflag,
ter_follow, ter_follow,
psafile, psafile,
errstat, errstat,

Write Preview
Loading…
Cancel
Save