Cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
C
C DoPhase is a subroutine written and provided by Jim Ramer at NOAA/FSL
C
C Ramer, J, 1993: An empirical technique for diagnosing precipitation
C type from model output. Preprints, 5th Conf. on Aviation
C Weather Systems, Vienna, VA, Amer. Meteor. Soc., 227-230.
C
Cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
C
SUBROUTINE CALWXT_RAMER(ttq,qq,ppq,pint,nq,lm,ppt,ptyp)
c SUBROUTINE dophase(pq, ! input pressure sounding mb
c + tq, ! input temperature sounding K
c + pq, | input pressure
c + qq, ! input spec humidityfraction
c + nq, ! input number of levels in sounding
c + twq, ! output wet-bulb sounding K
c + icefrac, ! output ice fraction
c + ptyp) ! output(2) phase 2=Rain, 3=Frzg, 4=Solid,
C 6=IP JC 9/16/99
PARAMETER (A2=17.2693882,A3=273.16,A4=35.86,PQ0=379.90516)
PARAMETER (G=9.80665,CP=1004.686,RCP=0.2857141,LECP=1572.5)
PARAMETER (twice=266.55,rhprcp=0.80,deltag=1.02,prcpmin=0.3,
* emelt=0.045,rlim=0.04,slim=0.85)
PARAMETER (twmelt=273.15,tz=273.15,efac=1.0,PTHRES=0.02)
C
INTEGER*4 i, k1, lll, k2, toodry, iflag, nq
C
INTEGER ptyp
C
REAL rcp, flg, flag, xxx, pq(lm), tq(lm), twq(lm), rhq(lm), mye,
* qq(lm), icefrac, tqtmp(lm), pqtmp(lm), rhqtmp(lm),
* pint(lm+1), ttq(lm), ppq(lm), ppint(lm+1),
* pinttmp(lm+1)
C
COMMON /flagflg/ flag, flg
DATA iflag / -9/
C
C Initialize.
icefrac = flag
ptyp = 0
IF (PPT.LE.PTHRES) RETURN
C
C GSM compute RH, convert pressure to mb, and reverse order
DO 88 i = 1, nq
LEV=NQ-I+1
QC=PQ0/PPQ(I) * EXP(A2*(TTQ(I)-A3)/(TTQ(I)-A4))
RHQTMP(LEV)=QQ(I)/QC
PQTMP(LEV)=PPQ(I)/100.
TQTMP(LEV)=TTQ(I)
PINTTMP(LEV)=PINT(I)
88 CONTINUE
do 92 i=1,nq
TQ(I)=TQTMP(I)
PQ(I)=PQTMP(I)
RHQ(I)=RHQTMP(I)
PPINT(I)=PINTTMP(I)
92 continue
C See if there was too little precip reported.
C
CCC RATE RESTRICTION REMOVED BY JOHN CORTINAS 3/16/99
C
C Construct wet-bulb sounding, locate generating level.
twmax = -999.0
rhmax = 0.0
k1 = 0 ! top of precip generating layer
k2 = 0 ! layer of maximum rh
C
IF (rhq(1).lt.rhprcp) THEN
toodry = 1
ELSE
toodry = 0
END IF
C
C toodry=((Rhq(1).lt.rhprcp).and.1)
lmhk=nq
pbot = ppint(lmhk+1)
DO 10 i = 1, nq
xxx = tdofesat(esat(tq(i))*rhq(i))
twq(i) = xmytw(tq(i),xxx,pq(i))
twmax = amax1(twq(i),twmax)
IF (pq(i).ge.400.0) THEN
IF (rhq(i).gt.rhmax) THEN
rhmax = rhq(i)
k2 = i
END IF
C
IF (i.ne.1) THEN
IF (rhq(i).ge.rhprcp.or.toodry.eq.0) THEN
IF (toodry.ne.0) THEN
dpdrh = alog(pq(i)/pq(i-1)) / (rhq(i)-
+ rhq(i-1))
pbot = exp(alog(pq(i))+(rhprcp-rhq(i))*dpdrh)
C
ptw = pq(i)
toodry = 0
ELSE IF (rhq(i).ge.rhprcp) THEN
ptw = pq(i)
ELSE
toodry = 1
dpdrh = alog(pq(i)/pq(i-1)) / (rhq(i)-
+ rhq(i-1))
ptw = exp(alog(pq(i))+(rhprcp-rhq(i))*dpdrh)
C
END IF
C
IF (pbot/ptw.ge.deltag) THEN
k1 = i
ptop = ptw
END IF
END IF
END IF
END IF
C
10 CONTINUE
C Gross checks for liquid and solid precip which dont require generating level.
C
IF (twq(1).ge.273.15+2.0) THEN
ptyp = 8 ! liquid
icefrac = 0.0
RETURN
END IF
C
IF (twmax.le.twice) THEN
icefrac = 1.0
ptyp = 1 ! solid
RETURN
END IF
C
C Check to see if we had no success with locating a generating level.
C
IF (k1.eq.0) THEN
RETURN
END IF
C
IF (ptop.eq.pq(k1)) THEN
twtop = twq(k1)
rhtop = rhq(k1)
k2 = k1
k1 = k1 - 1
ELSE
k2 = k1
k1 = k1 - 1
wgt1 = alog(ptop/pq(k2)) / alog(pq(k1)/pq(k2))
Clin wgt1=(ptop-Pq(k2))/(Pq(k1)-Pq(k2))
wgt2 = 1.0 - wgt1
twtop = twq(k1) * wgt1 + twq(k2) * wgt2
rhtop = rhq(k1) * wgt1 + rhq(k2) * wgt2
END IF
C
C Calculate temp and wet-bulb ranges below precip generating level.
DO 20 i = 1, k1
twmax = amax1(twq(i),twmax)
20 CONTINUE
C
C Gross check for solid precip, initialize ice fraction.
IF (twtop.le.twice) THEN
icefrac = 1.0
IF (twmax.le.twmelt) THEN ! gross check for solid precip.
ptyp = 1 ! solid precip
RETURN
END IF
lll = 0
ELSE
icefrac = 0.0
lll = 1
END IF
C
C Loop downward through sounding from highest precip generating level.
30 CONTINUE
C
IF (icefrac.ge.1.0) THEN ! starting as all ice
IF (twq(k1).lt.twmelt) GO TO 40 ! cannot commence melting
IF (twq(k1).eq.twtop) GO TO 40 ! both equal twmelt, nothing h
wgt1 = (twmelt-twq(k1)) / (twtop-twq(k1))
rhavg = rhq(k1) + wgt1 * (rhtop-rhq(k1)) / 2
dtavg = (twmelt-twq(k1)) / 2
dpk = wgt1 * alog(pq(k1)/ptop) !lin dpk=wgt1*(Pq(k1)-Ptop)
C mye=emelt*(1.0-(1.0-Rhavg)*efac)
mye = emelt * rhavg ** efac
icefrac = icefrac + dpk * dtavg / mye
ELSE IF (icefrac.le.0.0) THEN ! starting as all liquid
lll = 1
C If (Twq(k1).le.Twice) icefrac=1.0 ! autoconvert
C Goto 1020
IF (twq(k1).gt.twice) GO TO 40 ! cannot commence freezing
IF (twq(k1).eq.twtop) THEN
wgt1 = 0.5
ELSE
wgt1 = (twice-twq(k1)) / (twtop-twq(k1))
END IF
rhavg = rhq(k1) + wgt1 * (rhtop-rhq(k1)) / 2
dtavg = twmelt - (twq(k1)+twice) / 2
dpk = wgt1 * alog(pq(k1)/ptop) !lin dpk=wgt1*(Pq(k1)-Ptop)
C mye=emelt*(1.0-(1.0-Rhavg)*efac)
mye = emelt * rhavg ** efac
icefrac = icefrac + dpk * dtavg / mye
ELSE IF ((twq(k1).le.twmelt).and.(twq(k1).lt.twmelt)) THEN ! mix
rhavg = (rhq(k1)+rhtop) / 2
dtavg = twmelt - (twq(k1)+twtop) / 2
dpk = alog(pq(k1)/ptop) !lin dpk=Pq(k1)-Ptop
C mye=emelt*(1.0-(1.0-Rhavg)*efac)
mye = emelt * rhavg ** efac
icefrac = icefrac + dpk * dtavg / mye
ELSE ! mix where Tw curve crosses twmelt in layer
IF (twq(k1).eq.twtop) GO TO 40 ! both equal twmelt, nothing h
wgt1 = (twmelt-twq(k1)) / (twtop-twq(k1))
wgt2 = 1.0 - wgt1
rhavg = rhtop + wgt2 * (rhq(k1)-rhtop) / 2
dtavg = (twmelt-twtop) / 2
dpk = wgt2 * alog(pq(k1)/ptop) !lin dpk=wgt2*(Pq(k1)-Ptop)
C mye=emelt*(1.0-(1.0-Rhavg)*efac)
mye = emelt * rhavg ** efac
icefrac = icefrac + dpk * dtavg / mye
icefrac = amin1(1.0,amax1(icefrac,0.0))
IF (icefrac.le.0.0) THEN
C If (Twq(k1).le.Twice) icefrac=1.0 ! autoconvert
C Goto 1020
IF (twq(k1).gt.twice) GO TO 40 ! cannot commence freezin
wgt1 = (twice-twq(k1)) / (twtop-twq(k1))
dtavg = twmelt - (twq(k1)+twice) / 2
ELSE
dtavg = (twmelt-twq(k1)) / 2
END IF
rhavg = rhq(k1) + wgt1 * (rhtop-rhq(k1)) / 2
dpk = wgt1 * alog(pq(k1)/ptop) !lin dpk=wgt1*(Pq(k1)-Ptop)
C mye=emelt*(1.0-(1.0-Rhavg)*efac)
mye = emelt * rhavg ** efac
icefrac = icefrac + dpk * dtavg / mye
END IF
C
icefrac = amin1(1.0,amax1(icefrac,0.0))
C
C Get next level down if there is one, loop back.
40 IF (k1.gt.1) THEN
twtop = twq(k1)
ptop = pq(k1)
rhtop = rhq(k1)
k1 = k1 - 1
GO TO 30
END IF
C
C
C Determine precip type based on snow fraction and surface wet-bulb.
C
C
IF (icefrac.ge.slim) THEN
IF (lll.ne.0) THEN
ptyp = 2 ! Ice Pellets JC 9/16/99
ELSE
ptyp = 1 ! Snow
END IF
ELSE IF (icefrac.le.rlim) THEN
IF (twq(1).lt.tz) THEN
ptyp = 4 ! Freezing Precip
ELSE
ptyp = 8 ! Rain
END IF
ELSE
IF (twq(1).lt.tz) THEN
cGSM not sure what to do when 'mix' is predicted; I chose sleet as
cGSM a shaky best option. The use of dominant type should
CGSM help with these situation. Future users may want to
CGSM make the type 0 in this situation and leave the decision
CGSM up to the other algorithms.
ptyp = 2 ! Ice Pellets
c ptyp = 5 ! Mix
ELSE
c ptyp = 5 ! Mix
ptyp = 2 ! Ice Pellets
END IF
END IF
RETURN
C
END
C
REAL*4 FUNCTION esat(t)
C
C* Calculates saturation vapor pressure in millibars as a function of
C* either Kelvin of Celcius temperature.
C
IMPLICIT NONE
C
REAL*4 t, k
C
REAL*4 flag, flg
COMMON /flagflg/ flag, flg
C
C Account for both Celsius and Kelvin.
k = t
IF (k.lt.100.) k = k + 273.15
C
C Flag ridiculous values.
IF (k.lt.0.0.or.k.gt.373.15) THEN
esat = flag
RETURN
END IF
C
C Avoid floating underflow.
IF (k.lt.173.15) THEN
esat = 3.777647E-05
RETURN
END IF
C
C Calculation for normal range of values.
esat = exp(26.660820-0.0091379024*k-6106.3960/k)
C
RETURN
END
C
REAL*4 FUNCTION tdofesat(es)
C
C* As a function of saturation vapor pressure in millibars, returns
C* dewpoint in degrees K.
C
IMPLICIT NONE
C
REAL*4 es, lim1, lim2, b
C
DATA lim1, lim2 /3.777647E-05, 980.5386/
C
REAL*4 flag, flg
COMMON /flagflg/ flag, flg
C
C Flag ridiculous values.
IF (es.lt.0.0.or.es.gt.lim2) THEN
tdofesat = flag
RETURN
END IF
C
C Avoid floating underflow.
IF (es.lt.lim1) THEN
tdofesat = 173.15
RETURN
END IF
C
C Calculations for normal range of values.
b = 26.66082 - alog(es)
tdofesat = (b-sqrt(b*b-223.1986)) / 0.0182758048
C
RETURN
END
C
c REAL*4 FUNCTION mytw(t,td,p)
FUNCTION xmytw(t,td,p)
C
IMPLICIT NONE
C
INTEGER*4 cflag, l
REAL*4 f, c0, c1, c2, k, kd, kw, ew, t, td, p, ed, fp, s,
* de, xmytw
DATA f, c0, c1, c2 /0.0006355, 26.66082, 0.0091379024, 6106.3960/
C
C
xmytw = (t+td) / 2
IF (td.ge.t) RETURN
C
IF (t.lt.100.0) THEN
k = t + 273.15
kd = td + 273.15
IF (kd.ge.k) RETURN
cflag = 1
ELSE
k = t
kd = td
cflag = 0
END IF
C
ed = c0 - c1 * kd - c2 / kd
IF (ed.lt.-14.0.or.ed.gt.7.0) RETURN
ed = exp(ed)
ew = c0 - c1 * k - c2 / k
IF (ew.lt.-14.0.or.ew.gt.7.0) RETURN
ew = exp(ew)
fp = p * f
s = (ew-ed) / (k-kd)
kw = (k*fp+kd*s) / (fp+s)
C
DO 10 l = 1, 5
ew = c0 - c1 * kw - c2 / kw
IF (ew.lt.-14.0.or.ew.gt.7.0) RETURN
ew = exp(ew)
de = fp * (k-kw) + ed - ew
IF (abs(de/ew).lt.1E-5) GO TO 20
s = ew * (c1-c2/(kw*kw)) - fp
kw = kw - de / s
10 CONTINUE
20 CONTINUE
C
IF (cflag.ne.0) THEN
xmytw = kw - 273.15
ELSE
xmytw = kw
END IF
C
RETURN
END