UPP  V11.0.0
 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 ! 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(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