UPP  001
 All Data Structures Files Functions Pages
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 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
14 !
15  SUBROUTINE calwxt_ramer_post(T,Q,PMID,PINT,LMH,PREC,PTYP)
16 
17 ! SUBROUTINE dophase(pq, ! input pressure sounding mb
18 ! + t, ! input temperature sounding K
19 ! + pmid, ! input pressure
20 ! + pint, ! input interface pressure
21 ! + q, ! input spec humidityfraction
22 ! + lmh, ! input number of levels in sounding
23 ! + prec, ! input amount of precipitation
24 ! + ptyp) ! output(2) phase 2=Rain, 3=Frzg, 4=Solid,
25 ! 6=IP JC 9/16/99
26  use params_mod, only: pq0, a2, a3, a4
27  use ctlblk_mod, only: me, im, jsta_2l, jend_2u, lm, lp1, jsta, jend, pthresh
28 !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
29  implicit none
30 !
31  LOGICAL,parameter :: trace = .false.
32 ! real,PARAMETER :: RCP=0.2857141,LECP=1572.5
33  real,PARAMETER :: twice=266.55,rhprcp=0.80,deltag=1.02,prcpmin=0.3, &
34  & emelt=0.045,rlim=0.04,slim=0.85
35  real,PARAMETER :: twmelt=273.15,tz=273.15,efac=1.0 ! specify in params now
36 !
37  INTEGER*4 i, k1, lll, k2, toodry, iflag, nq
38 !
39  REAL xxx ,mye, icefrac,flg,flag
40  real,DIMENSION(IM,jsta_2l:jend_2u,LM), intent(in) :: t,q,pmid
41  real,DIMENSION(IM,jsta_2l:jend_2u,LP1),intent(in) :: pint
42  real,DIMENSION(IM,jsta_2l:jend_2u), intent(in) :: lmh,prec
43  integer,DIMENSION(IM,jsta:jend), intent(inout) :: ptyp
44 !
45  real,DIMENSION(IM,jsta_2l:jend_2u,LM) :: p,tq,pq,rhq
46  real,DIMENSION(IM,jsta:jend,LM) :: twq
47 ! REAL, ALLOCATABLE :: TWET(:,:,:)
48 !
49  integer j,l,lev,lnq,lmhk,ii
50  real rhmax,twmax,ptop,dpdrh,twtop,rhtop,wgt1,wgt2, &
51  rhavg,dtavg,dpk,ptw,rate,pbot,qc, b
52  real,external :: xmytw_post,esat,tdofesat
53 !
54  DATA iflag / -9/
55 !
56 ! Initialize.
57  IF (trace) WRITE (20,*) '******* NEW STATION ******'
58  IF (trace) WRITE (20,*) 'Twmelt,Twice,rhprcp,Emelt'
59  IF (trace) WRITE (20,*) twmelt, twice, rhprcp, emelt
60  flag=0.;flg=0.
61  icefrac = flag
62 !
63  DO j=jsta,jend
64  DO i=1,im
65  ptyp(i,j) = 0
66  nq=lmh(i,j)
67  DO l = 1,nq
68  lev = nq-l+1
69  p(i,j,l) = pmid(i,j,l)
70  qc = pq0/p(i,j,l) * exp(a2*(t(i,j,l)-a3)/(t(i,j,l)-a4))
71  tq(i,j,lev) = t(i,j,l)
72  pq(i,j,lev) = p(i,j,l)/100.
73  rhq(i,j,lev) = max(0.0, q(i,j,l)/qc)
74  enddo
75  enddo
76  enddo
77 
78 ! BIG LOOP
79  DO 800 j=jsta,jend
80  DO 800 i=1,im
81 !
82 ! SKIP THIS POINT IF NO PRECIP THIS TIME STEP
83 !
84  IF (prec(i,j)<=pthresh) cycle
85  lmhk=nint(lmh(i,j))
86 
87 !
88 !
89 !CC RATE RESTRICTION REMOVED BY JOHN CORTINAS 3/16/99
90 !
91 ! Construct wet-bulb sounding, locate generating level.
92  twmax = -999.0
93  rhmax = 0.0
94  k1 = 0 ! top of precip generating layer
95  k2 = 0 ! layer of maximum rh
96 !
97  IF (trace) WRITE (20,*) 'rhq(1)', rhq(i,j,1),'me=',me
98  IF (rhq(i,j,1)<rhprcp) THEN
99  toodry = 1
100  ELSE
101  toodry = 0
102  END IF
103 !
104 ! toodry=((Rhq(I,J,1)<rhprcp).and.1)
105  pbot = pq(i,j,1)
106  nq=lmh(i,j)
107  DO 10 l = 1, nq
108 ! xxx = tdofesat(esat(tq(I,J,L),flag,flg)*rhq(I,J,L),flag,flg)
109  xxx = max(0.0,min(pq(i,j,l),esat(tq(i,j,l),flag,flg))*rhq(i,j,l))
110  xxx = tdofesat(xxx,flag,flg)
111  twq(i,j,l) = xmytw_post(tq(i,j,l),xxx,pq(i,j,l))
112 
113 !` IF(I == 324 .and. J == 390) THEN
114 ! print *, 'tw ramer ', L, Twq(I,J,L),'me=',me
115 ! ENDIF
116  IF (trace) WRITE (*,*) 'Twq(I,J,L),L ', twq(i,j,l), l,'me=',me
117  twmax = max(twq(i,j,l),twmax)
118  IF (trace) WRITE (*,*) 'Tw,Rh,P ', twq(i,j,l) - 273.15, &
119  rhq(i,j,l), pq(i,j,l),'me=',me
120  IF (pq(i,j,l)>=400.0) THEN
121  IF (rhq(i,j,l)>rhmax) THEN
122  rhmax = rhq(i,j,l)
123  k2 = l
124  IF (trace) WRITE (*,*) 'rhmax,k2,L', rhmax, k2, l
125  END IF
126 !
127  IF (l/=1) THEN
128  IF (trace) WRITE (*,*) 'ME: toodry,L', toodry, l
129  IF (rhq(i,j,l)>=rhprcp.or.toodry==0) THEN
130  IF (toodry/=0) THEN
131  dpdrh = alog(pq(i,j,l)/pq(i,j,l-1)) / &
132  (rhq(i,j,l)-rhq(i,j,l-1))
133  pbot = exp(alog(pq(i,j,l))+(rhprcp-rhq(i,j,l))*dpdrh)
134 !
135 !lin dpdrh=(Pq(I,J,L)-Pq(I,J,L-1))/
136 !lin (Rhq(I,J,L)-Rhq(I,J,L-1))
137 !lin pbot=Pq(I,J,L)+(rhprcp-Rhq(I,J,L))*dpdrh
138  ptw = pq(i,j,l)
139  toodry = 0
140  IF (trace) WRITE (*,*) &
141  'dpdrh,pbot,rhprcp-rhq(I,J,L),L,ptw,toodry', &
142  dpdrh, pbot, rhprcp - rhq(i,j,l), &
143  l,ptw,toodry
144  ELSE IF (rhq(i,j,l)>=rhprcp) THEN
145  ptw = pq(i,j,l)
146  IF (trace) WRITE (*,*) 'HERE1: ptw,toodry',ptw, toodry
147  ELSE
148  toodry = 1
149  dpdrh = alog(pq(i,j,l)/pq(i,j,l-1)) / &
150  (rhq(i,j,l)-rhq(i,j,l-1))
151  ptw = exp(alog(pq(i,j,l))+(rhprcp-rhq(i,j,l))*dpdrh)
152  IF (trace) WRITE (*,*)'HERE2:dpdrh,pbot,L,ptw,toodry', &
153  dpdrh, pbot, l, ptw, toodry
154 !lin dpdrh=(Pq(i)-Pq(i-1))/(Rhq(i)-Rhq(i-1))
155 !lin ptw=Pq(i)+(rhprcp-Rhq(i))*dpdrh
156 !
157  END IF
158 !
159  IF (trace) WRITE (*,*) 'HERE3:pbot,ptw,deltag', &
160  & pbot, ptw, deltag
161  IF (pbot/ptw>=deltag) THEN
162 !lin If (pbot-ptw<deltag) Goto 2003
163  k1 = l
164  ptop = ptw
165  END IF
166  END IF
167  END IF
168  END IF
169 !
170  10 CONTINUE
171 
172 !
173 ! Gross checks for liquid and solid precip which dont require generating level.
174 !
175  IF (twq(i,j,1)>=273.15+2.0) THEN
176  ptyp(i,j) = 8 ! liquid
177  IF (trace) print *, 'liquid',i,j,'me=',me
178  icefrac = 0.0
179  cycle
180  END IF
181 !
182  IF (twmax<=twice) THEN
183  icefrac = 1.0
184  ptyp(i,j) = 1 ! solid
185  cycle
186  END IF
187 !
188 ! Check to see if we had no success with locating a generating level.
189 !
190  IF (trace) WRITE (*,*) 'HERE6: k1,ptyp', k1, ptyp(i,j),'me=',me
191  IF (k1==0) THEN
192  rate = flag
193  cycle
194  END IF
195 !
196  IF (ptop==pq(i,j,k1)) THEN
197  twtop = twq(i,j,k1)
198  rhtop = rhq(i,j,k1)
199  k2 = k1
200  k1 = k1 - 1
201  ELSE
202  k2 = k1
203  k1 = k1 - 1
204  wgt1 = alog(ptop/pq(i,j,k2)) / alog(pq(i,j,k1)/pq(i,j,k2))
205 !lin wgt1=(ptop-Pq(I,J,k2))/(Pq(I,J,k1)-Pq(I,J,k2))
206  wgt2 = 1.0 - wgt1
207  twtop = twq(i,j,k1) * wgt1 + twq(i,j,k2) * wgt2
208  rhtop = rhq(i,j,k1) * wgt1 + rhq(i,j,k2) * wgt2
209  END IF
210 !
211 
212 ! Calculate temp and wet-bulb ranges below precip generating level.
213  DO 20 l = 1, k1
214  twmax = amax1(twq(i,j,l),twmax)
215  20 CONTINUE
216 !
217 ! Gross check for solid precip, initialize ice fraction.
218  IF (i==1.and.j==1) WRITE (*,*) 'twmax=',twmax,twice,'twtop=',twtop
219  IF (twtop<=twice) THEN
220  icefrac = 1.0
221  IF (twmax<=twmelt) THEN ! gross check for solid precip.
222  IF (trace) print *, 'solid'
223  ptyp(i,j) = 1 ! solid precip
224  cycle
225  END IF
226  lll = 0
227  ELSE
228  icefrac = 0.0
229  lll = 1
230  END IF
231 !
232 ! Loop downward through sounding from highest precip generating level.
233  30 CONTINUE
234 !
235  IF (trace) print *, ptop, twtop - 273.15, icefrac,'me=',me
236  IF (trace) WRITE (*,*) 'P,Tw,frac,twq(I,J,k1)', ptop, &
237  & twtop - 273.15, icefrac, twq(i,j,k1),'me=',me
238  IF (icefrac>=1.0) THEN ! starting as all ice
239  IF (trace) WRITE (*,*) 'ICEFRAC=1', icefrac
240 ! print *, 'twq twmwelt twtop ', twq(I,J,k1), twmelt, twtop
241  IF (twq(i,j,k1)<twmelt) go to 40 ! cannot commence melting
242  IF (twq(i,j,k1)==twtop) go to 40 ! both equal twmelt, nothing h
243  wgt1 = (twmelt-twq(i,j,k1)) / (twtop-twq(i,j,k1))
244  rhavg = rhq(i,j,k1) + wgt1 * (rhtop-rhq(i,j,k1)) / 2
245  dtavg = (twmelt-twq(i,j,k1)) / 2
246  dpk = wgt1 * alog(pq(i,j,k1)/ptop) !lin dpk=wgt1*(Pq(k1)-Ptop)
247 ! mye=emelt*(1.0-(1.0-Rhavg)*efac)
248  mye = emelt * rhavg ** efac
249  icefrac = icefrac + dpk * dtavg / mye
250  IF (trace) WRITE (*,*) &
251  & 'HERE8: wgt1,rhavg,dtavg,dpk,mye,icefrac', wgt1, rhavg, &
252  & dtavg, dpk, mye, icefrac,'me=',me
253  ELSE IF (icefrac<=0.0) THEN ! starting as all liquid
254  IF (trace) WRITE (*,*) 'HERE9: twtop,twq(I,J,k1),k1,lll' &
255  & , twtop, twq(i,j,k1), k1, lll
256  lll = 1
257 ! If (Twq(I,J,k1)<=Twice) icefrac=1.0 ! autoconvert
258 ! Goto 1020
259  IF (twq(i,j,k1)>twice) go to 40 ! cannot commence freezing
260  IF (twq(i,j,k1)==twtop) THEN
261  wgt1 = 0.5
262  ELSE
263  wgt1 = (twice-twq(i,j,k1)) / (twtop-twq(i,j,k1))
264  END IF
265  rhavg = rhq(i,j,k1) + wgt1 * (rhtop-rhq(i,j,k1)) / 2
266  dtavg = twmelt - (twq(i,j,k1)+twice) / 2
267  dpk = wgt1 * alog(pq(i,j,k1)/ptop) !lin dpk=wgt1*(Pq(k1)-Ptop)
268 ! mye=emelt*(1.0-(1.0-Rhavg)*efac)
269  mye = emelt * rhavg ** efac
270  icefrac = icefrac + dpk * dtavg / mye
271  IF (trace) WRITE (*,*) 'HERE10: wgt1,rhtop,rhq(I,J,k1),dtavg', &
272  wgt1, rhtop, rhq(i,j,k1), dtavg,'me=',me
273  ELSE IF ((twq(i,j,k1)<=twmelt).and.(twq(i,j,k1)<twmelt)) THEN ! mix
274  rhavg = (rhq(i,j,k1)+rhtop) / 2
275  dtavg = twmelt - (twq(i,j,k1)+twtop) / 2
276  dpk = alog(pq(i,j,k1)/ptop) !lin dpk=Pq(I,J,k1)-Ptop
277 ! mye=emelt*(1.0-(1.0-Rhavg)*efac)
278  mye = emelt * rhavg ** efac
279  icefrac = icefrac + dpk * dtavg / mye
280 
281  IF (trace) WRITE (*,*) 'HERE11: twq(i,j,K1),twtop', &
282  twq(i,j,k1),twtop,'me=',me
283  ELSE ! mix where Tw curve crosses twmelt in layer
284  IF (twq(i,j,k1)==twtop) go to 40 ! both equal twmelt, nothing h
285  wgt1 = (twmelt-twq(i,j,k1)) / (twtop-twq(i,j,k1))
286  wgt2 = 1.0 - wgt1
287  rhavg = rhtop + wgt2 * (rhq(i,j,k1)-rhtop) / 2
288  dtavg = (twmelt-twtop) / 2
289  dpk = wgt2 * alog(pq(i,j,k1)/ptop) !lin dpk=wgt2*(Pq(k1)-Ptop)
290 ! mye=emelt*(1.0-(1.0-Rhavg)*efac)
291  mye = emelt * rhavg ** efac
292  icefrac = icefrac + dpk * dtavg / mye
293  icefrac = amin1(1.0,amax1(icefrac,0.0))
294  IF (trace) WRITE (*,*) 'HERE12: twq(I,J,k1),twtop,icefrac,wgt1,wg'// &
295  't2,rhavg,rhtop,rhq(I,J,k1),dtavg,k1', &
296  twq(i,j,k1), twtop, &
297  icefrac,wgt1,wgt2, rhavg, rhtop, rhq(i,j,k1), dtavg, k1 ,'me=',me
298  IF (icefrac<=0.0) THEN
299 ! If (Twq(I,J,k1)<=Twice) icefrac=1.0 ! autoconvert
300 ! Goto 1020
301  IF (twq(i,j,k1)>twice) go to 40 ! cannot commence freezin
302  wgt1 = (twice-twq(i,j,k1)) / (twtop-twq(i,j,k1))
303  dtavg = twmelt - (twq(i,j,k1)+twice) / 2
304  IF (trace) WRITE (*,*) 'IN IF','me=',me
305  ELSE
306  dtavg = (twmelt-twq(i,j,k1)) / 2
307  IF (trace) WRITE (*,*) 'IN ELSE','me=',me
308  END IF
309  IF (trace) WRITE (*,*) 'NEW ICE FRAC CALC','me=',me
310  rhavg = rhq(i,j,k1) + wgt1 * (rhtop-rhq(i,j,k1)) / 2
311  dpk = wgt1 * alog(pq(i,j,k1)/ptop) !lin dpk=wgt1*(Pq(k1)-Ptop)
312 ! mye=emelt*(1.0-(1.0-Rhavg)*efac)
313  mye = emelt * rhavg ** efac
314  icefrac = icefrac + dpk * dtavg / mye
315  IF (trace) WRITE (*,*) 'HERE13: icefrac,k1,dtavg,rhavg', &
316  icefrac, k1, dtavg, rhavg,'me=',me
317  END IF
318 !
319  icefrac = amin1(1.0,amax1(icefrac,0.0))
320  IF (i==1.and.j==1) WRITE (*,*) 'NEW ICEFRAC:', icefrac, icefrac,'me=',me
321 !
322 ! Get next level down if there is one, loop back.
323  40 IF (k1>1) THEN
324  IF (trace) WRITE (*,*) 'LOOPING BACK','me=',me
325  twtop = twq(i,j,k1)
326  ptop = pq(i,j,k1)
327  rhtop = rhq(i,j,k1)
328  k1 = k1 - 1
329  go to 30
330  END IF
331 !
332 !
333 ! Determine precip type based on snow fraction and surface wet-bulb.
334 !
335  IF (trace) WRITE (*,*) 'P,Tw,frac,lll', pq(i,j,k1), &
336  twq(i,j,k2) - 273.15, icefrac, lll,'me=',me
337 !
338  IF (icefrac>=slim) THEN
339  IF (lll/=0) THEN
340  ptyp(i,j) = 2 ! Ice Pellets JC 9/16/99
341  IF (trace) WRITE (*,*) 'frozen',i,j,'me=',me
342  ELSE
343  ptyp(i,j) = 1 ! Snow
344  IF (trace) WRITE (*,*) 'snow',i,j,'me=',me
345  END IF
346  ELSE IF (icefrac<=rlim) THEN
347  IF (twq(i,j,1)<tz) THEN
348  ptyp(i,j) = 4 ! Freezing Precip
349  IF (trace) WRITE (*,*) 'freezing',i,j,'me=',me
350  ELSE
351  ptyp(i,j) = 8 ! Rain
352  IF (trace) WRITE (*,*) 'liquid',i,j,'me=',me
353  END IF
354  ELSE
355  IF (trace) WRITE (*,*) 'Mix',i,j
356  IF (twq(i,j,1)<tz) THEN
357  IF (trace) WRITE (*,*) 'freezing',i,j
358 !GSM not sure what to do when 'mix' is predicted; In previous
359 !GSM versions of this code for which I had to have an answer,
360 !GSM I chose sleet. Here, though, since we have 4 other
361 !GSM algorithms to provide an answer, I will not declare a
362 !GSM type from the Ramer in this situation and allow the
363 !GSM other algorithms to make the call.
364 
365  ptyp(i,j) = 0 ! don't know
366 ! ptyp = 5 ! Mix
367  ELSE
368 ! ptyp = 5 ! Mix
369  ptyp(i,j) = 0 ! don't know
370  END IF
371  END IF
372  IF (trace) WRITE (*,*) "Returned ptyp is:ptyp,lll ", ptyp, lll,'me=',me
373  IF (trace) WRITE (*,*) "Returned icefrac is: ", icefrac,'me=',me
374  800 CONTINUE
375  DO 900 j=jsta,jend
376  DO 900 i=1,im
377  900 CONTINUE
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(0,*)' 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
params_mod
Definition: params.F:1