UPP (develop)
Loading...
Searching...
No Matches
CALWXT_RAMER.f
1!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2!
3! DoPhase is a subroutine written and provided by Jim Ramer at NOAA/FSL
4!
5! Ramer, J, 1993: An empirical technique for diagnosing precipitation
6! type from model output. Preprints, 5th Conf. on Aviation
7! Weather Systems, Vienna, VA, Amer. Meteor. Soc., 227-230.
8!
9! CODE ADAPTED FOR WRF POST 24 AUGUST 2005 G MANIKIN
10!
11! PROGRAM HISTORY LOG:
12! 10-30-19 Bo CUI - Remove "GOTO" statement
13! 21-10-31 JESSE MENG - 2D DECOMPOSITION
14!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
15!
16 SUBROUTINE calwxt_ramer_post(T,Q,PMID,PINT,LMH,PREC,PTYP)
17
18! SUBROUTINE dophase(pq, ! input pressure sounding mb
19! + t, ! input temperature sounding K
20! + pmid, ! input pressure
21! + pint, ! input interface pressure
22! + q, ! input spec humidityfraction
23! + lmh, ! input number of levels in sounding
24! + prec, ! input amount of precipitation
25! + ptyp) ! output(2) phase 2=Rain, 3=Frzg, 4=Solid,
26! 6=IP JC 9/16/99
27 use params_mod, only: pq0, a2, a3, a4
28 use ctlblk_mod, only: me, im, jsta_2l, jend_2u, lm, lp1, jsta, jend, pthresh,&
29 ista, iend, ista_2l, iend_2u
30!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
31 implicit none
32!
33 LOGICAL,parameter :: trace = .false.
34! real,PARAMETER :: RCP=0.2857141,LECP=1572.5
35 real,PARAMETER :: twice=266.55,rhprcp=0.80,deltag=1.02,prcpmin=0.3, &
36 & emelt=0.045,rlim=0.04,slim=0.85
37 real,PARAMETER :: twmelt=273.15,tz=273.15,efac=1.0 ! specify in params now
38!
39 INTEGER*4 i, k1, lll, k2, toodry, iflag, nq
40!
41 REAL xxx ,mye, icefrac,flg,flag
42 real,DIMENSION(ista_2l:iend_2u,jsta_2l:jend_2u,LM), intent(in) :: T,Q,PMID
43 real,DIMENSION(ista_2l:iend_2u,jsta_2l:jend_2u,LP1),intent(in) :: PINT
44 real,DIMENSION(ista_2l:iend_2u,jsta_2l:jend_2u), intent(in) :: LMH,PREC
45 integer,DIMENSION(ista:iend,jsta:jend), intent(inout) :: PTYP
46!
47 real,DIMENSION(ista_2l:iend_2u,jsta_2l:jend_2u,LM) :: P,TQ,PQ,RHQ
48 real,DIMENSION(ista:iend,jsta:jend,LM) :: TWQ
49! REAL, ALLOCATABLE :: TWET(:,:,:)
50!
51 integer J,L,LEV,LNQ,LMHK,ii
52 real RHMAX,TWMAX,PTOP,dpdrh,twtop,rhtop,wgt1,wgt2, &
53 rhavg,dtavg,dpk,ptw,rate,pbot,qc, b
54 real,external :: xmytw_post,esat,tdofesat
55!
56 DATA iflag / -9/
57!
58! Initialize.
59 IF (trace) WRITE (20,*) '******* NEW STATION ******'
60 IF (trace) WRITE (20,*) 'Twmelt,Twice,rhprcp,Emelt'
61 IF (trace) WRITE (20,*) twmelt, twice, rhprcp, emelt
62 flag=0.;flg=0.
63 icefrac = flag
64!
65 DO j=jsta,jend
66 DO i=ista,iend
67 ptyp(i,j) = 0
68 nq=lmh(i,j)
69 DO l = 1,nq
70 lev = nq-l+1
71 p(i,j,l) = pmid(i,j,l)
72 qc = pq0/p(i,j,l) * exp(a2*(t(i,j,l)-a3)/(t(i,j,l)-a4))
73 tq(i,j,lev) = t(i,j,l)
74 pq(i,j,lev) = p(i,j,l)/100.
75 rhq(i,j,lev) = max(0.0, q(i,j,l)/qc)
76 enddo
77 enddo
78 enddo
79
80! BIG LOOP
81 DO 800 j=jsta,jend
82 DO 800 i=ista,iend
83!
84! SKIP THIS POINT IF NO PRECIP THIS TIME STEP
85!
86 IF (prec(i,j)<=pthresh) cycle
87 lmhk=nint(lmh(i,j))
88
89!
90!
91!CC RATE RESTRICTION REMOVED BY JOHN CORTINAS 3/16/99
92!
93! Construct wet-bulb sounding, locate generating level.
94 twmax = -999.0
95 rhmax = 0.0
96 k1 = 0 ! top of precip generating layer
97 k2 = 0 ! layer of maximum rh
98!
99 IF (trace) WRITE (20,*) 'rhq(1)', rhq(i,j,1),'me=',me
100 IF (rhq(i,j,1)<rhprcp) THEN
101 toodry = 1
102 ELSE
103 toodry = 0
104 END IF
105!
106! toodry=((Rhq(I,J,1)<rhprcp).and.1)
107 pbot = pq(i,j,1)
108 nq=lmh(i,j)
109 DO 10 l = 1, nq
110! xxx = tdofesat(esat(tq(I,J,L),flag,flg)*rhq(I,J,L),flag,flg)
111 xxx = max(0.0,min(pq(i,j,l),esat(tq(i,j,l),flag,flg))*rhq(i,j,l))
112 xxx = tdofesat(xxx,flag,flg)
113 twq(i,j,l) = xmytw_post(tq(i,j,l),xxx,pq(i,j,l))
114
115!` IF(I == 324 .and. J == 390) THEN
116! print *, 'tw ramer ', L, Twq(I,J,L),'me=',me
117! ENDIF
118 IF (trace) WRITE (*,*) 'Twq(I,J,L),L ', twq(i,j,l), l,'me=',me
119 twmax = max(twq(i,j,l),twmax)
120 IF (trace) WRITE (*,*) 'Tw,Rh,P ', twq(i,j,l) - 273.15, &
121 rhq(i,j,l), pq(i,j,l),'me=',me
122 IF (pq(i,j,l)>=400.0) THEN
123 IF (rhq(i,j,l)>rhmax) THEN
124 rhmax = rhq(i,j,l)
125 k2 = l
126 IF (trace) WRITE (*,*) 'rhmax,k2,L', rhmax, k2, l
127 END IF
128!
129 IF (l/=1) THEN
130 IF (trace) WRITE (*,*) 'ME: toodry,L', toodry, l
131 IF (rhq(i,j,l)>=rhprcp.or.toodry==0) THEN
132 IF (toodry/=0) THEN
133 dpdrh = alog(pq(i,j,l)/pq(i,j,l-1)) / &
134 (rhq(i,j,l)-rhq(i,j,l-1))
135 pbot = exp(alog(pq(i,j,l))+(rhprcp-rhq(i,j,l))*dpdrh)
136!
137!lin dpdrh=(Pq(I,J,L)-Pq(I,J,L-1))/
138!lin (Rhq(I,J,L)-Rhq(I,J,L-1))
139!lin pbot=Pq(I,J,L)+(rhprcp-Rhq(I,J,L))*dpdrh
140 ptw = pq(i,j,l)
141 toodry = 0
142 IF (trace) WRITE (*,*) &
143 'dpdrh,pbot,rhprcp-rhq(I,J,L),L,ptw,toodry', &
144 dpdrh, pbot, rhprcp - rhq(i,j,l), &
145 l,ptw,toodry
146 ELSE IF (rhq(i,j,l)>=rhprcp) THEN
147 ptw = pq(i,j,l)
148 IF (trace) WRITE (*,*) 'HERE1: ptw,toodry',ptw, toodry
149 ELSE
150 toodry = 1
151 dpdrh = alog(pq(i,j,l)/pq(i,j,l-1)) / &
152 (rhq(i,j,l)-rhq(i,j,l-1))
153 ptw = exp(alog(pq(i,j,l))+(rhprcp-rhq(i,j,l))*dpdrh)
154 IF (trace) WRITE (*,*)'HERE2:dpdrh,pbot,L,ptw,toodry', &
155 dpdrh, pbot, l, ptw, toodry
156!lin dpdrh=(Pq(i)-Pq(i-1))/(Rhq(i)-Rhq(i-1))
157!lin ptw=Pq(i)+(rhprcp-Rhq(i))*dpdrh
158!
159 END IF
160!
161 IF (trace) WRITE (*,*) 'HERE3:pbot,ptw,deltag', &
162 & pbot, ptw, deltag
163 IF (pbot/ptw>=deltag) THEN
164!lin If (pbot-ptw<deltag) Goto 2003
165 k1 = l
166 ptop = ptw
167 END IF
168 END IF
169 END IF
170 END IF
171!
172 10 CONTINUE
173
174!
175! Gross checks for liquid and solid precip which dont require generating level.
176!
177 IF (twq(i,j,1)>=273.15+2.0) THEN
178 ptyp(i,j) = 8 ! liquid
179 IF (trace) print *, 'liquid',i,j,'me=',me
180 icefrac = 0.0
181 cycle
182 END IF
183!
184 IF (twmax<=twice) THEN
185 icefrac = 1.0
186 ptyp(i,j) = 1 ! solid
187 cycle
188 END IF
189!
190! Check to see if we had no success with locating a generating level.
191!
192 IF (trace) WRITE (*,*) 'HERE6: k1,ptyp', k1, ptyp(i,j),'me=',me
193 IF (k1==0) THEN
194 rate = flag
195 cycle
196 END IF
197!
198 IF (ptop==pq(i,j,k1)) THEN
199 twtop = twq(i,j,k1)
200 rhtop = rhq(i,j,k1)
201 k2 = k1
202 k1 = k1 - 1
203 ELSE
204 k2 = k1
205 k1 = k1 - 1
206 wgt1 = alog(ptop/pq(i,j,k2)) / alog(pq(i,j,k1)/pq(i,j,k2))
207!lin wgt1=(ptop-Pq(I,J,k2))/(Pq(I,J,k1)-Pq(I,J,k2))
208 wgt2 = 1.0 - wgt1
209 twtop = twq(i,j,k1) * wgt1 + twq(i,j,k2) * wgt2
210 rhtop = rhq(i,j,k1) * wgt1 + rhq(i,j,k2) * wgt2
211 END IF
212!
213
214! Calculate temp and wet-bulb ranges below precip generating level.
215 DO 20 l = 1, k1
216 twmax = amax1(twq(i,j,l),twmax)
217 20 CONTINUE
218!
219! Gross check for solid precip, initialize ice fraction.
220 IF (i==1.and.j==1) WRITE (*,*) 'twmax=',twmax,twice,'twtop=',twtop
221 IF (twtop<=twice) THEN
222 icefrac = 1.0
223 IF (twmax<=twmelt) THEN ! gross check for solid precip.
224 IF (trace) print *, 'solid'
225 ptyp(i,j) = 1 ! solid precip
226 cycle
227 END IF
228 lll = 0
229 ELSE
230 icefrac = 0.0
231 lll = 1
232 END IF
233!
234! Loop downward through sounding from highest precip generating level.
235 30 CONTINUE
236!
237 IF (trace) print *, ptop, twtop - 273.15, icefrac,'me=',me
238 IF (trace) WRITE (*,*) 'P,Tw,frac,twq(I,J,k1)', ptop, &
239 & twtop - 273.15, icefrac, twq(i,j,k1),'me=',me
240 IF (icefrac>=1.0) THEN ! starting as all ice
241 IF (trace) WRITE (*,*) 'ICEFRAC=1', icefrac
242! print *, 'twq twmwelt twtop ', twq(I,J,k1), twmelt, twtop
243 IF (twq(i,j,k1)<twmelt) GO TO 40 ! cannot commence melting
244 IF (twq(i,j,k1)==twtop) GO TO 40 ! both equal twmelt, nothing h
245 wgt1 = (twmelt-twq(i,j,k1)) / (twtop-twq(i,j,k1))
246 rhavg = rhq(i,j,k1) + wgt1 * (rhtop-rhq(i,j,k1)) / 2
247 dtavg = (twmelt-twq(i,j,k1)) / 2
248 dpk = wgt1 * alog(pq(i,j,k1)/ptop) !lin dpk=wgt1*(Pq(k1)-Ptop)
249! mye=emelt*(1.0-(1.0-Rhavg)*efac)
250 mye = emelt * rhavg ** efac
251 icefrac = icefrac + dpk * dtavg / mye
252 IF (trace) WRITE (*,*) &
253 & 'HERE8: wgt1,rhavg,dtavg,dpk,mye,icefrac', wgt1, rhavg, &
254 & dtavg, dpk, mye, icefrac,'me=',me
255 ELSE IF (icefrac<=0.0) THEN ! starting as all liquid
256 IF (trace) WRITE (*,*) 'HERE9: twtop,twq(I,J,k1),k1,lll' &
257 & , twtop, twq(i,j,k1), k1, lll
258 lll = 1
259! If (Twq(I,J,k1)<=Twice) icefrac=1.0 ! autoconvert
260! Goto 1020
261 IF (twq(i,j,k1)>twice) GO TO 40 ! cannot commence freezing
262 IF (twq(i,j,k1)==twtop) THEN
263 wgt1 = 0.5
264 ELSE
265 wgt1 = (twice-twq(i,j,k1)) / (twtop-twq(i,j,k1))
266 END IF
267 rhavg = rhq(i,j,k1) + wgt1 * (rhtop-rhq(i,j,k1)) / 2
268 dtavg = twmelt - (twq(i,j,k1)+twice) / 2
269 dpk = wgt1 * alog(pq(i,j,k1)/ptop) !lin dpk=wgt1*(Pq(k1)-Ptop)
270! mye=emelt*(1.0-(1.0-Rhavg)*efac)
271 mye = emelt * rhavg ** efac
272 icefrac = icefrac + dpk * dtavg / mye
273 IF (trace) WRITE (*,*) 'HERE10: wgt1,rhtop,rhq(I,J,k1),dtavg', &
274 wgt1, rhtop, rhq(i,j,k1), dtavg,'me=',me
275 ELSE IF ((twq(i,j,k1)<=twmelt).and.(twq(i,j,k1)<twmelt)) THEN ! mix
276 rhavg = (rhq(i,j,k1)+rhtop) / 2
277 dtavg = twmelt - (twq(i,j,k1)+twtop) / 2
278 dpk = alog(pq(i,j,k1)/ptop) !lin dpk=Pq(I,J,k1)-Ptop
279! mye=emelt*(1.0-(1.0-Rhavg)*efac)
280 mye = emelt * rhavg ** efac
281 icefrac = icefrac + dpk * dtavg / mye
282
283 IF (trace) WRITE (*,*) 'HERE11: twq(i,j,K1),twtop', &
284 twq(i,j,k1),twtop,'me=',me
285 ELSE ! mix where Tw curve crosses twmelt in layer
286 IF (twq(i,j,k1)==twtop) GO TO 40 ! both equal twmelt, nothing h
287 wgt1 = (twmelt-twq(i,j,k1)) / (twtop-twq(i,j,k1))
288 wgt2 = 1.0 - wgt1
289 rhavg = rhtop + wgt2 * (rhq(i,j,k1)-rhtop) / 2
290 dtavg = (twmelt-twtop) / 2
291 dpk = wgt2 * alog(pq(i,j,k1)/ptop) !lin dpk=wgt2*(Pq(k1)-Ptop)
292! mye=emelt*(1.0-(1.0-Rhavg)*efac)
293 mye = emelt * rhavg ** efac
294 icefrac = icefrac + dpk * dtavg / mye
295 icefrac = amin1(1.0,amax1(icefrac,0.0))
296 IF (trace) WRITE (*,*) 'HERE12: twq(I,J,k1),twtop,icefrac,wgt1,wg'// &
297 't2,rhavg,rhtop,rhq(I,J,k1),dtavg,k1', &
298 twq(i,j,k1), twtop, &
299 icefrac,wgt1,wgt2, rhavg, rhtop, rhq(i,j,k1), dtavg, k1 ,'me=',me
300 IF (icefrac<=0.0) THEN
301! If (Twq(I,J,k1)<=Twice) icefrac=1.0 ! autoconvert
302! Goto 1020
303 IF (twq(i,j,k1)>twice) GO TO 40 ! cannot commence freezin
304 wgt1 = (twice-twq(i,j,k1)) / (twtop-twq(i,j,k1))
305 dtavg = twmelt - (twq(i,j,k1)+twice) / 2
306 IF (trace) WRITE (*,*) 'IN IF','me=',me
307 ELSE
308 dtavg = (twmelt-twq(i,j,k1)) / 2
309 IF (trace) WRITE (*,*) 'IN ELSE','me=',me
310 END IF
311 IF (trace) WRITE (*,*) 'NEW ICE FRAC CALC','me=',me
312 rhavg = rhq(i,j,k1) + wgt1 * (rhtop-rhq(i,j,k1)) / 2
313 dpk = wgt1 * alog(pq(i,j,k1)/ptop) !lin dpk=wgt1*(Pq(k1)-Ptop)
314! mye=emelt*(1.0-(1.0-Rhavg)*efac)
315 mye = emelt * rhavg ** efac
316 icefrac = icefrac + dpk * dtavg / mye
317 IF (trace) WRITE (*,*) 'HERE13: icefrac,k1,dtavg,rhavg', &
318 icefrac, k1, dtavg, rhavg,'me=',me
319 END IF
320!
321 icefrac = amin1(1.0,amax1(icefrac,0.0))
322 IF (i==1.and.j==1) WRITE (*,*) 'NEW ICEFRAC:', icefrac, icefrac,'me=',me
323!
324! Get next level down if there is one, loop back.
325 40 IF (k1>1) THEN
326 IF (trace) WRITE (*,*) 'LOOPING BACK','me=',me
327 twtop = twq(i,j,k1)
328 ptop = pq(i,j,k1)
329 rhtop = rhq(i,j,k1)
330 k1 = k1 - 1
331 GO TO 30
332 END IF
333!
334!
335! Determine precip type based on snow fraction and surface wet-bulb.
336!
337 IF (trace) WRITE (*,*) 'P,Tw,frac,lll', pq(i,j,k1), &
338 twq(i,j,k2) - 273.15, icefrac, lll,'me=',me
339!
340 IF (icefrac>=slim) THEN
341 IF (lll/=0) THEN
342 ptyp(i,j) = 2 ! Ice Pellets JC 9/16/99
343 IF (trace) WRITE (*,*) 'frozen',i,j,'me=',me
344 ELSE
345 ptyp(i,j) = 1 ! Snow
346 IF (trace) WRITE (*,*) 'snow',i,j,'me=',me
347 END IF
348 ELSE IF (icefrac<=rlim) THEN
349 IF (twq(i,j,1)<tz) THEN
350 ptyp(i,j) = 4 ! Freezing Precip
351 IF (trace) WRITE (*,*) 'freezing',i,j,'me=',me
352 ELSE
353 ptyp(i,j) = 8 ! Rain
354 IF (trace) WRITE (*,*) 'liquid',i,j,'me=',me
355 END IF
356 ELSE
357 IF (trace) WRITE (*,*) 'Mix',i,j
358 IF (twq(i,j,1)<tz) THEN
359 IF (trace) WRITE (*,*) 'freezing',i,j
360!GSM not sure what to do when 'mix' is predicted; In previous
361!GSM versions of this code for which I had to have an answer,
362!GSM I chose sleet. Here, though, since we have 4 other
363!GSM algorithms to provide an answer, I will not declare a
364!GSM type from the Ramer in this situation and allow the
365!GSM other algorithms to make the call.
366
367 ptyp(i,j) = 0 ! don't know
368! ptyp = 5 ! Mix
369 ELSE
370! ptyp = 5 ! Mix
371 ptyp(i,j) = 0 ! don't know
372 END IF
373 END IF
374 IF (trace) WRITE (*,*) "Returned ptyp is:ptyp,lll ", ptyp, lll,'me=',me
375 IF (trace) WRITE (*,*) "Returned icefrac is: ", icefrac,'me=',me
376 800 CONTINUE
377
378 RETURN
379!
380 END
381!
382!--------------------------------------------------------------------------
383 REAL*4 FUNCTION esat(t,flag,flg)
384!
385!* Calculates saturation vapor pressure in millibars as a function of
386!* either Kelvin of Celcius temperature.
387!
388 IMPLICIT NONE
389!
390 real*4 t, k
391!
392 real*4 flag, flg
393!
394! Account for both Celsius and Kelvin.
395 k = t
396 IF (k<100.) k = k + 273.15
397!
398! Flag ridiculous values.
399 IF (k<0.0.or.k>373.15) THEN
400 esat = flag
401 RETURN
402 END IF
403!
404! Avoid floating underflow.
405 IF (k<173.15) THEN
406 esat = 3.777647e-05
407 RETURN
408 END IF
409!
410! Calculation for normal range of values.
411 esat = exp(26.660820-0.0091379024*k-6106.3960/k)
412!
413 RETURN
414 END
415!--------------------------------------------------------------------------
416 REAL*4 FUNCTION tdofesat(es,flag,flg)
417!
418!* As a function of saturation vapor pressure in millibars, returns
419!* dewpoint in degrees K.
420!
421 IMPLICIT NONE
422!
423 real*4 es, lim1, lim2, b
424!
425 DATA lim1, lim2 /3.777647e-05, 980.5386/
426!
427 real*4 flag, flg
428! COMMON /flagflg/ flag, flg
429!
430! Flag ridiculous values.
431 IF (es<0.0.or.es>lim2) THEN
432 tdofesat = flag
433 RETURN
434 END IF
435!
436! Avoid floating underflow.
437 IF (es<lim1) THEN
438 tdofesat = 173.15
439 RETURN
440 END IF
441!
442! Calculations for normal range of values.
443 b = 26.66082 - alog(es)
444 tdofesat = (b-sqrt(b*b-223.1986)) / 0.0182758048
445!
446 RETURN
447 END
448!
449!--------------------------------------------------------------------------
450! REAL*4 FUNCTION mytw(t,td,p)
451 FUNCTION xmytw_post(t,td,p)
452!
453 IMPLICIT NONE
454!
455 INTEGER*4 cflag, l,i,j
456 real*4 f, c0, c1, c2, k, kd, kw, ew, t, td, p, ed, fp, s, &
457 & de, xmytw_post
458 DATA f, c0, c1, c2 /0.0006355, 26.66082, 0.0091379024, 6106.3960/
459!
460!
461 xmytw_post= (t+td) / 2
462 IF (td>=t) RETURN
463!
464 IF (t<100.0) THEN
465 k = t + 273.15
466 kd = td + 273.15
467 IF (kd>=k) RETURN
468 cflag = 1
469 ELSE
470 k = t
471 kd = td
472 cflag = 0
473 END IF
474 if (kd == 0.0) write(*,*)' kd=',kd,' t=',t,' p=',p,' td=',td
475!
476 ed = c0 - c1 * kd - c2 / kd
477 IF (ed<-14.0.or.ed>7.0) RETURN
478 ed = exp(ed)
479 ew = c0 - c1 * k - c2 / k
480 IF (ew<-14.0.or.ew>7.0) RETURN
481 ew = exp(ew)
482 fp = p * f
483 s = (ew-ed) / (k-kd)
484 kw = (k*fp+kd*s) / (fp+s)
485!
486 DO 10 l = 1, 5
487 ew = c0 - c1 * kw - c2 / kw
488 IF (ew<-14.0.or.ew>7.0) RETURN
489 ew = exp(ew)
490 de = fp * (k-kw) + ed - ew
491 IF (abs(de/ew)<1e-5) exit
492 s = ew * (c1-c2/(kw*kw)) - fp
493 kw = kw - de / s
494 10 CONTINUE
495!
496! print *, 'kw ', kw
497 IF (cflag/=0) THEN
498 xmytw_post= kw - 273.15
499 ELSE
500 xmytw_post = kw
501 END IF
502!
503 RETURN
504 END