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

Code cleanup

main
Bill Ladwig 9 years ago
parent
5a8b146d38
commit
66881a1ad0
1 changed files with 125 additions and 134 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. 169
      wrf_open/var/src/python/wrf/var/wrfcape.f90

169
wrf_open/var/src/python/wrf/var/wrfcape.f90
Unescape View File

@ -24,7 +24,7 @@ REAL(KIND=8) FUNCTION tvirtual(temp,ratmix)
REAL(KIND=8),INTENT(IN) :: temp,ratmix REAL(KIND=8),INTENT(IN) :: temp,ratmix
REAL(KIND=8),PARAMETER :: EPS = .622D0 REAL(KIND=8),PARAMETER :: EPS = .622D0
tvirtual = temp*(EPS+ratmix)/(EPS*(1.D0+ratmix)) tvirtual = temp*(EPS + ratmix)/(EPS*(1.D0 + ratmix))
RETURN RETURN
END FUNCTION tvirtual END FUNCTION tvirtual
@ -65,7 +65,7 @@ REAL(KIND=8) FUNCTION tonpsadiabat(thte,prs,PSADITHTE,PSADIPRS,PSADITMK,GAMMA,&
jt = -1 jt = -1
DO jtch = 1,150 - 1 DO jtch = 1,150 - 1
IF (thte.GE.PSADITHTE(jtch) .AND. thte.LT.PSADITHTE(jtch+1)) THEN IF (thte .GE. PSADITHTE(jtch) .AND. thte .LT. PSADITHTE(jtch+1)) THEN
jt = jtch jt = jtch
EXIT EXIT
!GO TO 213 !GO TO 213
@ -94,8 +94,8 @@ REAL(KIND=8) FUNCTION tonpsadiabat(thte,prs,PSADITHTE,PSADIPRS,PSADITMK,GAMMA,&
fracjt2 = 1.D0 - fracjt fracjt2 = 1.D0 - fracjt
fracip = (PSADIPRS(ip)-prs) / (PSADIPRS(ip)-PSADIPRS(ip+1)) fracip = (PSADIPRS(ip)-prs) / (PSADIPRS(ip)-PSADIPRS(ip+1))
fracip2 = 1.D0 - fracip fracip2 = 1.D0 - fracip
IF (PSADITMK(ip,jt).GT.1D9 .OR. PSADITMK(ip+1,jt).GT.1D9 .OR. & IF (PSADITMK(ip,jt) .GT. 1D9 .OR. PSADITMK(ip+1,jt) .GT. 1D9 .OR. &
PSADITMK(ip,jt+1).GT.1D9 .OR. PSADITMK(ip+1,jt+1).GT.1D9) THEN PSADITMK(ip,jt+1) .GT. 1D9 .OR. PSADITMK(ip+1,jt+1) .GT. 1D9) THEN
CALL throw_exception('capecalc3d: ','Tried to access missing temperature in lookup table.',& CALL throw_exception('capecalc3d: ','Tried to access missing temperature in lookup table.',&
'Prs and Thte probably unreasonable. prs,thte=',prs,thte) 'Prs and Thte probably unreasonable. prs,thte=',prs,thte)
!STOP !STOP
@ -137,16 +137,16 @@ SUBROUTINE dpfcalc(prs,sfp,pf,miy,mjx,mkzh,ter_follow)
! do i=1,miy-1 staggered grid ! 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(i,j,k) = 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(i,j,k) = .5D0 * (3.D0*prs(i,j,k) - prs(i,j,k-1))
END IF END IF
ELSE ELSE
pf(i,j,k) = .5D0* (prs(i,j,k+1)+prs(i,j,k)) pf(i,j,k) = .5D0 * (prs(i,j,k+1) + prs(i,j,k))
END IF END IF
END DO END DO
END DO END DO
@ -289,24 +289,22 @@ SUBROUTINE f_computecape(prs,tmk,qvp,ght,ter,sfp,cape,cin,&
cape(i,j,1) = 0.d0 cape(i,j,1) = 0.d0
cin(i,j,1) = 0.d0 cin(i,j,1) = 0.d0
IF (i3dflag.eq.1) THEN IF (i3dflag .EQ. 1) THEN
kpar1 = 2 kpar1 = 2
kpar2 = mkzh kpar2 = mkzh
ELSE ELSE
! find parcel with max theta-e in lowest 3 km agl. ! find parcel with max theta-e in lowest 3 km agl.
ethmax = -1.d0 ethmax = -1.d0
DO k = mkzh,1,-1 DO k = mkzh,1,-1
IF (ght(i,j,k)-ter(i,j).lt.3000.d0) then IF (ght(i,j,k)-ter(i,j) .LT. 3000.d0) THEN
q = MAX(qvp(i,j,k),1.d-15) q = MAX(qvp(i,j,k), 1.d-15)
t = tmk(i,j,k) t = tmk(i,j,k)
p = prs(i,j,k) p = prs(i,j,k)
e = q*p/ (EPS+q) e = q*p/(EPS + q)
tlcl = TLCLC1/ (log(t**TLCLC2/e)-TLCLC3) + TLCLC4 tlcl = TLCLC1/ (LOG(t**TLCLC2/e)-TLCLC3) + TLCLC4
eth = t* (1000.d0/p)**(GAMMA* (1.d0+GAMMAMD*q))*& eth = t * (1000.d0/p)**(GAMMA*(1.d0 + GAMMAMD*q))*&
EXP((THTECON1/tlcl-THTECON2)*q*(1.d0+THTECON3*q)) EXP((THTECON1/tlcl - THTECON2)*q*(1.d0 + THTECON3*q))
IF (eth.gt.ethmax) then IF (eth .GT. ethmax) THEN
klev = k klev = k
ethmax = eth ethmax = eth
END IF END IF
@ -315,29 +313,29 @@ SUBROUTINE f_computecape(prs,tmk,qvp,ght,ter,sfp,cape,cin,&
kpar1 = klev kpar1 = klev
kpar2 = klev kpar2 = klev
! establish average properties of that parcel ! Establish average properties of that parcel
! (over depth of approximately davg meters) ! (over depth of approximately davg meters)
! davg=.1 ! davg=.1
davg = 500.d0 davg = 500.d0
pavg = davg*prs(i,j,kpar1)*& pavg = davg*prs(i,j,kpar1)*&
GRAV/(RGAS*tvirtual(tmk(i,j,kpar1),qvp(i,j,kpar1))) GRAV/(RGAS*tvirtual(tmk(i,j,kpar1), qvp(i,j,kpar1)))
p2 = MIN(prs(i,j,kpar1)+.5d0*pavg,prsf(i,j,mkzh)) p2 = MIN(prs(i,j,kpar1)+.5d0*pavg, prsf(i,j,mkzh))
p1 = p2 - pavg p1 = p2 - pavg
totthe = 0.d0 totthe = 0.d0
totqvp = 0.d0 totqvp = 0.d0
totprs = 0.d0 totprs = 0.d0
DO k = mkzh,2,-1 DO k = mkzh,2,-1
IF (prsf(i,j,k).le.p1) GOTO 35 IF (prsf(i,j,k) .LE. p1) GOTO 35
IF (prsf(i,j,k-1).ge.p2) GOTO 34 IF (prsf(i,j,k-1) .GE. p2) GOTO 34
p = prs(i,j,k) p = prs(i,j,k)
pup = prsf(i,j,k) pup = prsf(i,j,k)
pdn = prsf(i,j,k-1) pdn = prsf(i,j,k-1)
q = MAX(qvp(i,j,k),1.d-15) q = MAX(qvp(i,j,k),1.d-15)
th = tmk(i,j,k)* (1000.d0/p)**(GAMMA* (1.d0+GAMMAMD*q)) th = tmk(i,j,k)*(1000.d0/p)**(GAMMA*(1.d0 + GAMMAMD*q))
pp1 = MAX(p1,pdn) pp1 = MAX(p1,pdn)
pp2 = MIN(p2,pup) pp2 = MIN(p2,pup)
IF (pp2.gt.pp1) then IF (pp2 .GT. pp1) THEN
deltap = pp2 - pp1 deltap = pp2 - pp1
totqvp = totqvp + q*deltap totqvp = totqvp + q*deltap
totthe = totthe + th*deltap totthe = totthe + th*deltap
@ -353,53 +351,52 @@ SUBROUTINE f_computecape(prs,tmk,qvp,ght,ter,sfp,cape,cin,&
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 IF (i3dflag .EQ. 1) THEN
qvppari = qvp(i,j,kpar) qvppari = qvp(i,j,kpar)
tmkpari = tmk(i,j,kpar) tmkpari = tmk(i,j,kpar)
END IF END IF
prspari = prs(i,j,kpar) prspari = prs(i,j,kpar)
ghtpari = ght(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*prspari/(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/prspari)**(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)/ (GRAV/cpm) zlcl = ghtpari + (tmkpari - tlcl)/(GRAV/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
! initial parcel already saturated or supersaturated.
IF (ghtpari .GE. zlcl) THEN
! Initial parcel already saturated or supersaturated.
ilcl = 2 ilcl = 2
klcl = 1 klcl = 1
END IF END IF
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 33 kk = kk + 1
! model level is below lcl
! Model level is below lcl
IF (ght(i,j,k).lt.zlcl) THEN IF (ght(i,j,k).lt.zlcl) THEN
qvplift = qvppari qvplift = qvppari
tmklift = tmkpari - GRAV/cpm*(ght(i,j,k)-ghtpari) tmklift = tmkpari - GRAV/cpm*(ght(i,j,k) - ghtpari)
tvenv = tvirtual(tmk(i,j,k),qvp(i,j,k)) tvenv = tvirtual(tmk(i,j,k), qvp(i,j,k))
tvlift = tvirtual(tmklift,qvplift) tvlift = tvirtual(tmklift, qvplift)
ghtlift = ght(i,j,k) 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 tmklift = tlcl
qvplift = qvppari qvplift = qvppari
facden = ght(i,j,k) - ght(i,j,k+1) facden = ght(i,j,k) - ght(i,j,k+1)
@ -407,83 +404,78 @@ SUBROUTINE f_computecape(prs,tmk,qvp,ght,ter,sfp,cape,cin,&
fac2 = (ght(i,j,k)-zlcl)/facden fac2 = (ght(i,j,k)-zlcl)/facden
tmkenv = tmk(i,j,k+1)*fac2 + tmk(i,j,k)*fac1 tmkenv = tmk(i,j,k+1)*fac2 + tmk(i,j,k)*fac1
qvpenv = qvp(i,j,k+1)*fac2 + qvp(i,j,k)*fac1 qvpenv = qvp(i,j,k+1)*fac2 + qvp(i,j,k)*fac1
tvenv = tvirtual(tmkenv,qvpenv) tvenv = tvirtual(tmkenv, qvpenv)
tvlift = tvirtual(tmklift,qvplift) tvlift = tvirtual(tmklift, qvplift)
ghtlift = zlcl ghtlift = zlcl
ilcl = 1 ilcl = 1
ELSE ELSE
tmklift = tonpsadiabat(ethpari,prs(i,j,k),PSADITHTE,PSADIPRS,PSADITMK,GAMMA,throw_exception) tmklift = tonpsadiabat(ethpari, prs(i,j,k), PSADITHTE, PSADIPRS,&
PSADITMK, GAMMA, throw_exception)
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(i,j,k) - eslift)
tvenv = tvirtual(tmk(i,j,k),qvp(i,j,k)) tvenv = tvirtual(tmk(i,j,k), qvp(i,j,k))
tvlift = tvirtual(tmklift,qvplift) tvlift = tvirtual(tmklift, qvplift)
ghtlift = ght(i,j,k) ghtlift = ght(i,j,k)
END IF END IF
! buoyancy ! Buoyancy
buoy(kk) = GRAV* (tvlift-tvenv)/tvenv buoy(kk) = GRAV*(tvlift - tvenv)/tvenv
zrel(kk) = ghtlift - ghtpari zrel(kk) = ghtlift - ghtpari
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.
kk = kk + 1 kk = kk + 1
buoy(kk) = buoy(kk-1) buoy(kk) = buoy(kk-1)
zrel(kk) = zrel(kk-1) zrel(kk) = zrel(kk-1)
buoy(kk-1) = 0.d0 buoy(kk-1) = 0.d0
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 GOTO 33
END IF END IF
END DO END DO
kmax = kk kmax = kk
IF (kmax.gt.150) THEN IF (kmax .GT. 150) THEN
! Need an exception here ! Need an exception here
CALL throw_exception('capecalc3d: kmax got too big. kmax=',kmax) CALL throw_exception('capecalc3d: kmax got too big. kmax=',kmax)
!STOP !STOP
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 ! Get the accumulated buoyant energy from the parcel's starting
! get the accumulated buoyant energy from the parcel's starting
! point, at all levels up to the top level. ! point, at all levels up to the top level.
benaccum(1) = 0.0d0 benaccum(1) = 0.0d0
benamin = 9d9 benamin = 9d9
DO k = 2,kmax DO k = 2,kmax
dz = zrel(k) - zrel(k-1) dz = zrel(k) - zrel(k-1)
benaccum(k) = benaccum(k-1) +.5d0*dz* (buoy(k-1)+buoy(k)) benaccum(k) = benaccum(k-1) + .5d0*dz*(buoy(k-1) + buoy(k))
IF (benaccum(k).lt.benamin) then IF (benaccum(k) .LT. benamin) THEN
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.
DO k = kmax,klcl,-1 DO k = kmax,klcl,-1
IF (buoy(k).ge.0.d0) THEN IF (buoy(k) .GE. 0.d0) THEN
! k of equilibrium level ! k of equilibrium level
kel = k kel = k
GOTO 50 GOTO 50
END IF END IF
END DO END DO
! if we got through that loop, then there is no non-negative ! If we got through that loop, then there is no non-negative
! buoyancy above the lcl in the sounding. in these situations, ! 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 ! 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 ! missing values in v6.1.0). also, where cape is
@ -491,7 +483,7 @@ SUBROUTINE f_computecape(prs,tmk,qvp,ght,ter,sfp,cape,cin,&
! 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.
@ -506,7 +498,7 @@ SUBROUTINE f_computecape(prs,tmk,qvp,ght,ter,sfp,cape,cin,&
50 CONTINUE 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 the ! 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 ! el, but at or above the lcl, where accumulated buoyant energy is a
! minimum. the net positive area (accumulated buoyant energy) from ! minimum. the net positive area (accumulated buoyant energy) from
@ -515,38 +507,37 @@ SUBROUTINE f_computecape(prs,tmk,qvp,ght,ter,sfp,cape,cin,&
! parcel starting point to the lfc will be defined as the convective ! 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.
benamin = 9d9 benamin = 9d9
klfc = kmax klfc = kmax
DO k = klcl,kel DO k = klcl,kel
IF (benaccum(k).lt.benamin) THEN IF (benaccum(k) .LT. benamin) THEN
benamin = benaccum(k) benamin = benaccum(k)
klfc = k klfc = k
END IF END IF
END DO END DO
! now we can assign values to cape and cin ! Now we can assign values to cape and cin
cape(i,j,kpar) = MAX(benaccum(kel)-benamin,0.1d0) cape(i,j,kpar) = MAX(benaccum(kel)-benamin, 0.1d0)
cin(i,j,kpar) = MAX(-benamin,0.1d0) cin(i,j,kpar) = MAX(-benamin, 0.1d0)
! cin is uninteresting when cape is small (< 100 j/kg), so set ! 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 ! cin to -0.1 (see note about missing values in v6.1.0) in
! that case. ! that case.
! 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.
! 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) = -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 102 CONTINUE
END DO END DO
IF (i3dflag.eq.0) THEN 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

Write Preview
Loading…
Cancel
Save