GFS Operational Physics Documentation  Revision: 81451
precpd.f
Go to the documentation of this file.
1 
4 
50 
78  subroutine precpd (im,ix,km,dt,del,prsl,q,cwm,t,rn,sr
79  &, rainp,u00k,psautco,prautco,evpco,wminco
80  &, lprnt,jpr)
81 !
82 !
83 ! ******************************************************************
84 ! * *
85 ! * subroutine for precipitation processes *
86 ! * from suspended cloud water/ice *
87 ! * *
88 ! ******************************************************************
89 ! * *
90 ! * originally created by q. zhao jan. 1995 *
91 ! * ------- *
92 ! * modified and rewritten by shrinivas moorthi oct. 1998 *
93 ! * ----------------- *
94 ! * and hua-lu pan *
95 ! * ---------- *
96 ! * *
97 ! * references: *
98 ! * *
99 ! * zhao and carr (1997), monthly weather review (august) *
100 ! * sundqvist et al., (1989) monthly weather review. (august) *
101 ! * chuang 2013, modify sr to define frozen precipitation fraction*
102 ! ******************************************************************
103 !
104 ! in this code vertical indexing runs from surface to top of the
105 ! model
106 !
107 ! argument list:
108 ! --------------
109 ! im : inner dimension over which calculation is made
110 ! ix : maximum inner dimension
111 ! km : number of vertical levels
112 ! dt : time step in seconds
113 ! del(km) : pressure layer thickness (bottom to top)
114 ! prsl(km) : pressure values for model layers (bottom to top)
115 ! q(ix,km) : specific humidity (updated in the code)
116 ! cwm(ix,km) : condensate mixing ratio (updated in the code)
117 ! t(ix,km) : temperature (updated in the code)
118 ! rn(im) : precipitation over one time-step dt (m/dt)
119 !old sr(im) : index (=-1 snow, =0 rain/snow, =1 rain)
120 !new sr(im) : "snow ratio", ratio of snow to total precipitation
121 ! cll(ix,km) : cloud cover
122 !hchuang rn(im) unit in m per time step
123 ! precipitation rate conversion 1 mm/s = 1 kg/m2/s
124 !
125  use machine , only : kind_phys
126  use funcphys , only : fpvs
127  use physcons, grav => con_g, hvap => con_hvap, hfus => con_hfus
128  &, ttp => con_ttp, cp => con_cp
129  &, eps => con_eps, epsm1 => con_epsm1
130  implicit none
131 ! include 'constant.h'
132 !
133  real (kind=kind_phys) g, h1, h1000
134  &, d00
135  &, elwv, eliv, row
136  &, epsq, eliw
137  &, rcp, rrow
138  parameter(g=grav, h1=1.e0, h1000=1000.0
139  &, d00=0.e0
140  &, elwv=hvap, eliv=hvap+hfus, row=1.e3
141  &, epsq=2.e-12
142  &, eliw=eliv-elwv, rcp=h1/cp, rrow=h1/row)
143 !
144  real(kind=kind_phys), parameter :: cons_0=0.0, cons_p01=0.01
145  &, cons_20=20.0
146  &, cons_m30=-30.0, cons_50=50.0
147 !
148  integer im, ix, km, jpr
149  real (kind=kind_phys) q(ix,km), t(ix,km), cwm(ix,km)
150  &, del(ix,km), prsl(ix,km)
151  &, rn(im), sr(im)
152  &, dt
153  &, rainp(im,km), rnp(im),
154  & psautco(im), prautco(im), evpco, wminco(2)
155 !
156 !
157  real (kind=kind_phys) err(im), ers(im), precrl(im)
158  &, precsl(im), precrl1(im), precsl1(im)
159  &, rq(im), condt(im)
160  &, conde(im), rconde(im), tmt0(im)
161  &, wmin(im,km), wmink(im), pres(im)
162  &, wmini(im,km), ccr(im)
163  &, tt(im), qq(im), ww(im)
164  &, u00k(im,km)
165  &, zaodt
166  real (kind=kind_phys) cclim(km)
167 !
168  integer iw(im,km), ipr(im), iwl(im), iwl1(im)
169 !
170  logical comput(im)
171  logical lprnt
172 !
173  real (kind=kind_phys) ke, rdt, us, climit, cws, csm1
174  &, crs1, crs2, cr, aa2, dtcp, c00, cmr
175  &, tem, c1, c2, wwn
176 ! &, tem, c1, c2, u00b, u00t, wwn
177  &, precrk, precsk, pres1, qk, qw, qi
178  &, qint, fiw, wws, cwmk, expf
179  &, psaut, psaci, amaxcm, tem1, tem2
180  &, tmt0k, psm1, psm2, ppr
181  &, rprs, erk, pps, sid, rid, amaxps
182  &, praut, fi, qc, amaxrq, rqkll
183  integer i, k, ihpr, n
184 !
185 !-----------------------preliminaries ---------------------------------
186 !
187 ! do k=1,km
188 ! do i=1,im
189 ! cll(i,k) = 0.0
190 ! enddo
191 ! enddo
192 !
193  rdt = h1 / dt
194 ! ke = 2.0e-5 ! commented on 09/10/99 -- opr value
195 ! ke = 2.0e-6
196 ! ke = 1.0e-5
197 !!! ke = 5.0e-5
198 !! ke = 7.0e-5
199  ke = evpco
200 ! ke = 7.0e-5
201  us = h1
202  climit = 1.0e-20
203  cws = 0.025
204 !
205  zaodt = 800.0 * rdt
206 !
207  csm1 = 5.0000e-8 * zaodt
208  crs1 = 5.00000e-6 * zaodt
209  crs2 = 6.66600e-10 * zaodt
210  cr = 5.0e-4 * zaodt
211  aa2 = 1.25e-3 * zaodt
212 !
213  ke = ke * sqrt(rdt)
214 ! ke = ke * sqrt(zaodt)
215 !
216  dtcp = dt * rcp
217 !
218 ! c00 = 1.5e-1 * dt
219 ! c00 = 10.0e-1 * dt
220 ! c00 = 3.0e-1 * dt !05/09/2000
221 ! c00 = 1.0e-4 * dt !05/09/2000
222 ! c00 = prautco * dt !05/09/2000
223  cmr = 1.0 / 3.0e-4
224 ! cmr = 1.0 / 5.0e-4
225 ! c1 = 100.0
226  c1 = 300.0
227  c2 = 0.5
228 !
229 !
230 !--------calculate c0 and cmr using lc at previous step-----------------
231 !
232  do k=1,km
233  do i=1,im
234  tem = (prsl(i,k)*0.00001)
235 ! tem = sqrt(tem)
236  iw(i,k) = 0.0
237 ! wmin(i,k) = 1.0e-5 * tem
238 ! wmini(i,k) = 1.0e-5 * tem ! testing for ras
239 !
240 
241  wmin(i,k) = wminco(1) * tem
242  wmini(i,k) = wminco(2) * tem
243 
244 
245  rainp(i,k) = 0.0
246 
247  enddo
248  enddo
249  do i=1,im
250 ! c0(i) = 1.5e-1
251 ! cmr(i) = 3.0e-4
252 !
253  iwl1(i) = 0
254  precrl1(i) = d00
255  precsl1(i) = d00
256  comput(i) = .false.
257  rn(i) = d00
258  sr(i) = d00
259  ccr(i) = d00
260 !
261  rnp(i) = d00
262  enddo
274 
275 !------------select columns where rain can be produced--------------
276  do k=1, km-1
277  do i=1,im
278  tem = min(wmin(i,k), wmini(i,k))
279  if (cwm(i,k) > tem) comput(i) = .true.
280  enddo
281  enddo
282  ihpr = 0
283  do i=1,im
284  if (comput(i)) then
285  ihpr = ihpr + 1
286  ipr(ihpr) = i
287  endif
288  enddo
289 !***********************************************************************
290 !-----------------begining of precipitation calculation-----------------
291 !***********************************************************************
292 ! do k=km-1,2,-1
293  do k=km,1,-1
294  do n=1,ihpr
295  precrl(n) = precrl1(n)
296  precsl(n) = precsl1(n)
297  err(n) = d00
298  ers(n) = d00
299  iwl(n) = 0
300 !
301  i = ipr(n)
302  tt(n) = t(i,k)
303  qq(n) = q(i,k)
304  ww(n) = cwm(i,k)
305  wmink(n) = wmin(i,k)
306  pres(n) = prsl(i,k)
307 !
308  precrk = max(cons_0, precrl1(n))
309  precsk = max(cons_0, precsl1(n))
310  wwn = max(ww(n), climit)
311 ! if (wwn .gt. wmink(n) .or. (precrk+precsk) .gt. d00) then
312  if (wwn > climit .or. (precrk+precsk) > d00) then
313  comput(n) = .true.
314  else
315  comput(n) = .false.
316  endif
317  enddo
318 !
319 ! es(1:ihpr) = fpvs(tt(1:ihpr))
320  do n=1,ihpr
321  if (comput(n)) then
322  i = ipr(n)
323  conde(n) = (dt/g) * del(i,k)
324  condt(n) = conde(n) * rdt
325  rconde(n) = h1 / conde(n)
326  qk = max(epsq, qq(n))
327  tmt0(n) = tt(n) - 273.16
328  wwn = max(ww(n), climit)
329 !
330 ! pl = pres(n) * 0.01
331 ! call qsatd(tt(n), pl, qc)
332 ! rq(n) = max(qq(n), epsq) / max(qc, 1.0e-10)
333 ! rq(n) = max(1.0e-10, rq(n)) ! -- relative humidity---
334 !
335 ! the global qsat computation is done in pa
336  pres1 = pres(n)
337 ! qw = es(n)
338  qw = min(pres1, fpvs(tt(n)))
339  qw = eps * qw / (pres1 + epsm1 * qw)
340  qw = max(qw,epsq)
341 !
342 ! tmt15 = min(tmt0(n), cons_m15)
343 ! ai = 0.008855
344 ! bi = 1.0
345 ! if (tmt0(n) .lt. -20.0) then
346 ! ai = 0.007225
347 ! bi = 0.9674
348 ! endif
349 ! qi = qw * (bi + ai*min(tmt0(n),cons_0))
350 ! qint = qw * (1.-0.00032*tmt15*(tmt15+15.))
351 !
352  qi = qw
353  qint = qw
354 ! if (tmt0(n).le.-40.) qint = qi
355 !
356 !-------------------ice-water id number iw------------------------------
359  if(tmt0(n) < -15.) then
360  fi = qk - u00k(i,k)*qi
361  if(fi > d00 .or. wwn > climit) then
362  iwl(n) = 1
363  else
364  iwl(n) = 0
365  endif
366 ! endif
367  elseif (tmt0(n) >= 0.) then
368  iwl(n) = 0
369 !
370 ! if(tmt0(n).lt.0.0.and.tmt0(n).ge.-15.0) then
371  else
372  iwl(n) = 0
373  if(iwl1(n) == 1 .and. wwn > climit) iwl(n) = 1
374  endif
375 !
376 ! if(tmt0(n).ge.0.) then
377 ! iwl(n) = 0
378 ! endif
379 !----------------the satuation specific humidity------------------------
380  fiw = float(iwl(n))
381  qc = (h1-fiw)*qint + fiw*qi
382 !----------------the relative humidity----------------------------------
383  if(qc <= 1.0e-10) then
384  rq(n) = d00
385  else
386  rq(n) = qk / qc
387  endif
388 !----------------cloud cover ratio ccr----------------------------------
390  if(rq(n) < u00k(i,k)) then
391  ccr(n) = d00
392  elseif(rq(n) >= us) then
393  ccr(n) = us
394  else
395  rqkll = min(us,rq(n))
396  ccr(n) = h1-sqrt((us-rqkll)/(us-u00k(i,k)))
397  endif
398 !
399  endif
400  enddo
401 !-------------------ice-water id number iwl------------------------------
402 ! do n=1,ihpr
403 ! if (comput(n) .and. (ww(n) .gt. climit)) then
404 ! if (tmt0(n) .lt. -15.0
405 ! * .or. (tmt0(n) .lt. 0.0 .and. iwl1(n) .eq. 1))
406 ! * iwl(n) = 1
407 ! cll(ipr(n),k) = 1.0 ! cloud cover!
408 ! cll(ipr(n),k) = min(1.0, ww(n)*cclim(k)) ! cloud cover!
409 ! endif
410 ! enddo
411 !
443 !--- precipitation production -- auto conversion and accretion
444 !
445  do n=1,ihpr
446  if (comput(n) .and. ccr(n) > 0.0) then
447  wws = ww(n)
448  cwmk = max(cons_0, wws)
449  i = ipr(n)
450 ! amaxcm = max(cons_0, cwmk - wmink(n))
451  if (iwl(n) == 1) then ! ice phase
452  amaxcm = max(cons_0, cwmk - wmini(i,k))
453  expf = dt * exp(0.025*tmt0(n))
454  psaut = min(cwmk, psautco(i)*expf*amaxcm)
455  ww(n) = ww(n) - psaut
456  cwmk = max(cons_0, ww(n))
457 ! cwmk = max(cons_0, ww(n)-wmini(i,k))
458  psaci = min(cwmk, aa2*expf*precsl1(n)*cwmk)
459 
460  ww(n) = ww(n) - psaci
461  precsl(n) = precsl(n) + (wws - ww(n)) * condt(n)
462  else ! liquid water
463 !
472 ! for using sundqvist precip formulation of rain
473 !
474  amaxcm = max(cons_0, cwmk - wmink(n))
475 !! amaxcm = cwmk
476  tem1 = precsl1(n) + precrl1(n)
477  tem2 = min(max(cons_0, 268.0-tt(n)), cons_20)
478  tem = (1.0+c1*sqrt(tem1*rdt)) * (1+c2*sqrt(tem2))
479 !
480  tem2 = amaxcm * cmr * tem / max(ccr(n),cons_p01)
481  tem2 = min(cons_50, tem2*tem2)
482 ! praut = c00 * tem * amaxcm * (1.0-exp(-tem2))
483  praut = (prautco(i)*dt) * tem * amaxcm
484  & * (1.0-exp(-tem2))
485  praut = min(praut, cwmk)
486  ww(n) = ww(n) - praut
487 !
497 ! below is for zhao's precip formulation (water)
498 !
499 ! amaxcm = max(cons_0, cwmk - wmink(n))
500 ! praut = min(cwmk, c00*amaxcm*amaxcm)
501 ! ww(n) = ww(n) - praut
502 !
503 ! cwmk = max(cons_0, ww(n))
504 ! tem1 = precsl1(n) + precrl1(n)
505 ! pracw = min(cwmk, cr*dt*tem1*cwmk)
506 ! ww(n) = ww(n) - pracw
507 !
508  precrl(n) = precrl(n) + (wws - ww(n)) * condt(n)
509 !
510 !hchuang code change [+1l] : add record to record information in vertical
511 ! turn rnp in unit of ww (cwm and q, kg/kg ???)
512  rnp(n) = rnp(n) + (wws - ww(n))
513  endif
514  endif
515  enddo
538 !
539 !-----evaporation of precipitation-------------------------
540 !**** err & ers positive--->evaporation-- negtive--->condensation
541 !
542  do n=1,ihpr
543  if (comput(n)) then
544  i = ipr(n)
545  qk = max(epsq, qq(n))
546  tmt0k = max(cons_m30, tmt0(n))
547  precrk = max(cons_0, precrl(n))
548  precsk = max(cons_0, precsl(n))
549  amaxrq = max(cons_0, u00k(i,k)-rq(n)) * conde(n)
550 !----------------------------------------------------------------------
551 ! increase the evaporation for strong/light prec
552 !----------------------------------------------------------------------
553  ppr = ke * amaxrq * sqrt(precrk)
554 ! ppr = ke * amaxrq * sqrt(precrk*rdt)
555  if (tmt0(n) .ge. 0.) then
556  pps = 0.
557  else
558  pps = (crs1+crs2*tmt0k) * amaxrq * precsk / u00k(i,k)
559  end if
560 !---------------correct if over-evapo./cond. occurs--------------------
561  erk=precrk+precsk
562  if(rq(n).ge.1.0e-10) erk = amaxrq * qk * rdt / rq(n)
563  if (ppr+pps .gt. abs(erk)) then
564  rprs = erk / (precrk+precsk)
565  ppr = precrk * rprs
566  pps = precsk * rprs
567  endif
568  ppr = min(ppr, precrk)
569  pps = min(pps, precsk)
570  err(n) = ppr * rconde(n)
571  ers(n) = pps * rconde(n)
572  precrl(n) = precrl(n) - ppr
573 !hchuang code change [+1l] : add record to record information in vertical
574 ! use err for kg/kg/dt not the ppr (mm/dt=kg/m2/dt)
575 !
576  rnp(n) = rnp(n) - err(n)
577 !
578  precsl(n) = precsl(n) - pps
579  endif
580  enddo
612 !--------------------melting of the snow--------------------------------
613  do n=1,ihpr
614  if (comput(n)) then
615  if (tmt0(n) .gt. 0.) then
616  amaxps = max(cons_0, precsl(n))
617  psm1 = csm1 * tmt0(n) * tmt0(n) * amaxps
618  psm2 = cws * cr * max(cons_0, ww(n)) * amaxps
619  ppr = (psm1 + psm2) * conde(n)
620  if (ppr .gt. amaxps) then
621  ppr = amaxps
622  psm1 = amaxps * rconde(n)
623  endif
624  precrl(n) = precrl(n) + ppr
625 !
626 !hchuang code change [+1l] : add record to record information in vertical
627 ! turn ppr (mm/dt=kg/m2/dt) to kg/kg/dt -> ppr/air density (kg/m3)
628  rnp(n) = rnp(n) + ppr * rconde(n)
629 !
630  precsl(n) = precsl(n) - ppr
631  else
632  psm1 = d00
633  endif
634 !
635 !---------------update t and q------------------------------------------
643 
644  tt(n) = tt(n) - dtcp * (elwv*err(n)+eliv*ers(n)+eliw*psm1)
645  qq(n) = qq(n) + dt * (err(n)+ers(n))
646  endif
647  enddo
648 !
649  do n=1,ihpr
650  iwl1(n) = iwl(n)
651  precrl1(n) = max(cons_0, precrl(n))
652  precsl1(n) = max(cons_0, precsl(n))
653  i = ipr(n)
654  t(i,k) = tt(n)
655  q(i,k) = qq(n)
656  cwm(i,k) = ww(n)
657  iw(i,k) = iwl(n)
658 !hchuang code change [+1l] : add record to record information in vertical
659 ! rnp = precrl1*rconde(n) unit in kg/kg/dt
660 !
661  rainp(i,k) = rnp(n)
662  enddo
663 !
664 ! move water from vapor to liquid should the liquid amount be negative
665 !
666  do i = 1, im
667  if (cwm(i,k) < 0.) then
668  tem = q(i,k) + cwm(i,k)
669  if (tem >= 0.0) then
670  q(i,k) = tem
671  t(i,k) = t(i,k) - elwv * rcp * cwm(i,k)
672  cwm(i,k) = 0.
673  elseif (q(i,k) > 0.0) then
674  cwm(i,k) = tem
675  t(i,k) = t(i,k) + elwv * rcp * q(i,k)
676  q(i,k) = 0.0
677  endif
678  endif
679  enddo
680 !
681  enddo ! k loop ends here!
682 !**********************************************************************
683 !-----------------------end of precipitation processes-----------------
684 !**********************************************************************
685 !
694  do n=1,ihpr
695  i = ipr(n)
696  rn(i) = (precrl1(n) + precsl1(n)) * rrow ! precip at surface
697 !
698 !----sr=1 if sfc prec is rain ; ----sr=-1 if sfc prec is snow
699 !----sr=0 for both of them or no sfc prec
700 !
701 ! rid = 0.
702 ! sid = 0.
703 ! if (precrl1(n) .ge. 1.e-13) rid = 1.
704 ! if (precsl1(n) .ge. 1.e-13) sid = -1.
705 ! sr(i) = rid + sid ! sr=1 --> rain, sr=-1 -->snow, sr=0 -->both
706 ! chuang, june 2013: change sr to define fraction of frozen precipitation instead
707 ! because wpc uses it in their winter experiment
708 
709  rid = precrl1(n) + precsl1(n)
710  if (rid < 1.e-13) then
711  sr(i) = 0.
712  else
713  sr(i) = precsl1(n)/rid
714  endif
715  enddo
716 !
717  return
718  end
719 !! @}
real(kind=kind_phys), parameter con_g
gravity ( )
Definition: physcons.f:59
subroutine precpd(im, ix, km, dt, del, prsl, q, cwm, t, rn, sr , rainp, u00k, psautco, prautco, evpco, wminco , lprnt, jpr)
Definition: precpd.f:81
real(kind=kind_phys), parameter con_hfus
lat heat H2O fusion ( )
Definition: physcons.f:92
real(kind=kind_phys), parameter con_cp
spec heat air at p ( )
Definition: physcons.f:80
real(kind=kind_phys), parameter con_hvap
lat heat H2O cond ( )
Definition: physcons.f:90
real(kind=kind_phys), parameter con_ttp
temp at H2O 3pt (K)
Definition: physcons.f:98