CCPP SciDoc v7.0.0  v7.0.0
Common Community Physics Package Developed at DTC
 
Loading...
Searching...
No Matches
zhaocarr_precpd.f
1
4
7
8 implicit none
9 public :: zhaocarr_precpd_init, zhaocarr_precpd_run
10 private
11 logical :: is_initialized = .false.
12 contains
13
14 subroutine zhaocarr_precpd_init (imp_physics, &
15 & imp_physics_zhao_carr, &
16 & imp_physics_zhao_carr_pdf, &
17 & errmsg, errflg)
18 implicit none
19
20 ! Interface variables
21 integer, intent(in ) :: imp_physics
22 integer, intent(in ) :: imp_physics_zhao_carr, &
23 & imp_physics_zhao_carr_pdf
24 ! CCPP error handling
25 character(len=*), intent( out) :: errmsg
26 integer, intent( out) :: errflg
27
28 ! Initialize the CCPP error handling variables
29 errmsg = ''
30 errflg = 0
31
32 if (is_initialized) return
33
34 ! Consistency checks
35 if (imp_physics/=imp_physics_zhao_carr .and. &
36 & imp_physics/=imp_physics_zhao_carr_pdf) then
37 write(errmsg,'(*(a))') "Logic error: namelist choice of &
38 & microphysics is different from Zhao-Carr MP"
39 errflg = 1
40 return
41 end if
42
43 is_initialized = .true.
44 end subroutine zhaocarr_precpd_init
45
78 subroutine zhaocarr_precpd_run (im,km,dt,del,prsl,q,cwm,t,rn &
79 &, grav, hvap, hfus, ttp, cp, eps, epsm1 &
80 &, sr,rainp,u00k,psautco,prautco,evpco,wminco &
81 &, wk1,lprnt,jpr,errmsg,errflg)
82
83!
84! ******************************************************************
85! * *
86! * subroutine for precipitation processes *
87! * from suspended cloud water/ice *
88! * *
89! ******************************************************************
90! * *
91! * originally created by q. zhao jan. 1995 *
92! * ------- *
93! * modified and rewritten by shrinivas moorthi oct. 1998 *
94! * ----------------- *
95! * and hua-lu pan *
96! * ---------- *
97! * *
98! * references: *
99! * *
100! * zhao and carr (1997), monthly weather review (august) *
101! * sundqvist et al., (1989) monthly weather review. (august) *
102! * chuang 2013, modify sr to define frozen precipitation fraction*
103! ******************************************************************
104!
105! in this code vertical indexing runs from surface to top of the
106! model
107!
108! argument list:
109! --------------
110! im : inner dimension over which calculation is made
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(im,km) : specific humidity (updated in the code)
116! cwm(im,km) : condensate mixing ratio (updated in the code)
117! t(im,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(im,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
128 implicit none
129! include 'constant.h'
130!
131! Interface variables
132 integer, intent(in) :: im, km, jpr
133 real (kind=kind_phys), intent(in) :: grav, hvap, hfus, ttp, cp, &
134 & eps, epsm1
135 real (kind=kind_phys), intent(in) :: dt
136 real (kind=kind_phys), intent(in) :: del(:,:), prsl(:,:)
137 real (kind=kind_phys), intent(inout) :: q(:,:), t(:,:), &
138 & cwm(:,:)
139 real (kind=kind_phys), intent(out) :: rn(:), sr(:), rainp(:,:)
140 real (kind=kind_phys), intent(in) :: u00k(:,:)
141 real (kind=kind_phys), intent(in) :: psautco(:), prautco(:), &
142 & evpco, wminco(:), wk1(:)
143 logical, intent(in) :: lprnt
144 character(len=*), intent(out) :: errmsg
145 integer, intent(out) :: errflg
146!
147! Local variables
148 real (kind=kind_phys) g, h1, h1000
149 &, d00
150 &, elwv, eliv, row
151 &, epsq, eliw
152 &, rcp, rrow
153 parameter( h1=1.e0, h1000=1000.0 &
154 &, d00=0.e0, row=1.e3 &
155 &, epsq=2.e-12)
156!
157
158 real(kind=kind_phys), parameter :: cons_0=0.0, cons_p01=0.01 &
159 &, cons_20=20.0 &
160 &, cons_m30=-30.0, cons_50=50.0
161!
162 real (kind=kind_phys) rnp(im), psautco_l(im), prautco_l(im) &
163 &, wk2(im)
164!
165 real (kind=kind_phys) err(im), ers(im), precrl(im) &
166 &, precsl(im), precrl1(im), precsl1(im) &
167 &, rq(im), condt(im) &
168 &, conde(im), rconde(im), tmt0(im) &
169 &, wmin(im,km), wmink(im), pres(im) &
170 &, wmini(im,km), ccr(im) &
171 &, tt(im), qq(im), ww(im) &
172 &, zaodt
173 real (kind=kind_phys) cclim(km)
174!
175 integer iw(im,km), ipr(im), iwl(im), iwl1(im)
176!
177 logical comput(im)
178!
179 real (kind=kind_phys) ke, rdt, us, climit, cws, csm1
180 &, crs1, crs2, cr, aa2, dtcp, c00, cmr
181 &, tem, c1, c2, wwn
182! &, tem, c1, c2, u00b, u00t, wwn
183 &, precrk, precsk, pres1, qk, qw, qi
184 &, qint, fiw, wws, cwmk, expf
185 &, psaut, psaci, amaxcm, tem1, tem2
186 &, tmt0k, psm1, psm2, ppr
187 &, rprs, erk, pps, sid, rid, amaxps
188 &, praut, fi, qc, amaxrq, rqkll
189 integer i, k, ihpr, n
190!
191 ! Initialize CCPP error handling variables
192 errmsg = ''
193 errflg = 0
194!-------------- GFS psautco/prautco interstitial ----------------
195 do i=1, im
196 wk2(i) = 1.0-wk1(i)
197 psautco_l(i) = psautco(1)*wk1(i) + psautco(2)*wk2(i)
198 prautco_l(i) = prautco(1)*wk1(i) + prautco(2)*wk2(i)
199 enddo
200!-----------------------preliminaries ---------------------------------
201!
202! do k=1,km
203! do i=1,im
204! cll(i,k) = 0.0
205! enddo
206! enddo
207!
208 g=grav
209 elwv=hvap
210 eliv=hvap+hfus
211 eliw=eliv-elwv
212 rcp=h1/cp
213 rrow=h1/row
214 rdt = h1 / dt
215! ke = 2.0e-5 ! commented on 09/10/99 -- opr value
216! ke = 2.0e-6
217! ke = 1.0e-5
218!!! ke = 5.0e-5
219!! ke = 7.0e-5
220 ke = evpco
221! ke = 7.0e-5
222 us = h1
223 climit = 1.0e-20
224 cws = 0.025
225!
226 zaodt = 800.0 * rdt
227!
228 csm1 = 5.0000e-8 * zaodt
229 crs1 = 5.00000e-6 * zaodt
230 crs2 = 6.66600e-10 * zaodt
231 cr = 5.0e-4 * zaodt
232 aa2 = 1.25e-3 * zaodt
233!
234 ke = ke * sqrt(rdt)
235! ke = ke * sqrt(zaodt)
236!
237 dtcp = dt * rcp
238!
239! c00 = 1.5e-1 * dt
240! c00 = 10.0e-1 * dt
241! c00 = 3.0e-1 * dt !05/09/2000
242! c00 = 1.0e-4 * dt !05/09/2000
243! c00 = prautco * dt !05/09/2000
244 cmr = 1.0 / 3.0e-4
245! cmr = 1.0 / 5.0e-4
246! c1 = 100.0
247 c1 = 300.0
248 c2 = 0.5
249!
250!
251!--------calculate c0 and cmr using lc at previous step-----------------
252!
253 do k=1,km
254 do i=1,im
255 tem = (prsl(i,k)*0.00001)
256! tem = sqrt(tem)
257 iw(i,k) = 0.0
258! wmin(i,k) = 1.0e-5 * tem
259! wmini(i,k) = 1.0e-5 * tem ! testing for ras
260!
261
262 wmin(i,k) = wminco(1) * tem
263 wmini(i,k) = wminco(2) * tem
264
265
266 rainp(i,k) = 0.0
267
268 enddo
269 enddo
270 do i=1,im
271! c0(i) = 1.5e-1
272! cmr(i) = 3.0e-4
273!
274 iwl1(i) = 0
275 precrl1(i) = d00
276 precsl1(i) = d00
277 comput(i) = .false.
278 rn(i) = d00
279 sr(i) = d00
280 ccr(i) = d00
281!
282 rnp(i) = d00
283 enddo
295
296!------------select columns where rain can be produced--------------
297 do k=1, km-1
298 do i=1,im
299 tem = min(wmin(i,k), wmini(i,k))
300 if (cwm(i,k) > tem) comput(i) = .true.
301 enddo
302 enddo
303 ihpr = 0
304 do i=1,im
305 if (comput(i)) then
306 ihpr = ihpr + 1
307 ipr(ihpr) = i
308 endif
309 enddo
310!***********************************************************************
311!-----------------begining of precipitation calculation-----------------
312!***********************************************************************
313! do k=km-1,2,-1
314 do k=km,1,-1
315 do n=1,ihpr
316 precrl(n) = precrl1(n)
317 precsl(n) = precsl1(n)
318 err(n) = d00
319 ers(n) = d00
320 iwl(n) = 0
321!
322 i = ipr(n)
323 tt(n) = t(i,k)
324 qq(n) = q(i,k)
325 ww(n) = cwm(i,k)
326 wmink(n) = wmin(i,k)
327 pres(n) = prsl(i,k)
328!
329 precrk = max(cons_0, precrl1(n))
330 precsk = max(cons_0, precsl1(n))
331 wwn = max(ww(n), climit)
332! if (wwn .gt. wmink(n) .or. (precrk+precsk) .gt. d00) then
333 if (wwn > climit .or. (precrk+precsk) > d00) then
334 comput(n) = .true.
335 else
336 comput(n) = .false.
337 endif
338 enddo
339!
340! es(1:ihpr) = fpvs(tt(1:ihpr))
341 do n=1,ihpr
342 if (comput(n)) then
343 i = ipr(n)
344 conde(n) = (dt/g) * del(i,k)
345 condt(n) = conde(n) * rdt
346 rconde(n) = h1 / conde(n)
347 qk = max(epsq, qq(n))
348 tmt0(n) = tt(n) - 273.16
349 wwn = max(ww(n), climit)
350!
351! pl = pres(n) * 0.01
352! call qsatd(tt(n), pl, qc)
353! rq(n) = max(qq(n), epsq) / max(qc, 1.0e-10)
354! rq(n) = max(1.0e-10, rq(n)) ! -- relative humidity---
355!
356! the global qsat computation is done in pa
357 pres1 = pres(n)
358! qw = es(n)
359 qw = min(pres1, fpvs(tt(n)))
360 qw = eps * qw / (pres1 + epsm1 * qw)
361 qw = max(qw,epsq)
362!
363! tmt15 = min(tmt0(n), cons_m15)
364! ai = 0.008855
365! bi = 1.0
366! if (tmt0(n) .lt. -20.0) then
367! ai = 0.007225
368! bi = 0.9674
369! endif
370! qi = qw * (bi + ai*min(tmt0(n),cons_0))
371! qint = qw * (1.-0.00032*tmt15*(tmt15+15.))
372!
373 qi = qw
374 qint = qw
375! if (tmt0(n).le.-40.) qint = qi
376!
377!-------------------ice-water id number iw------------------------------
380 if(tmt0(n) < -15.) then
381 fi = qk - u00k(i,k)*qi
382 if(fi > d00 .or. wwn > climit) then
383 iwl(n) = 1
384 else
385 iwl(n) = 0
386 endif
387! endif
388 elseif (tmt0(n) >= 0.) then
389 iwl(n) = 0
390!
391! if(tmt0(n).lt.0.0.and.tmt0(n).ge.-15.0) then
392 else
393 iwl(n) = 0
394 if(iwl1(n) == 1 .and. wwn > climit) iwl(n) = 1
395 endif
396!
397! if(tmt0(n).ge.0.) then
398! iwl(n) = 0
399! endif
400!----------------the satuation specific humidity------------------------
401 fiw = float(iwl(n))
402 qc = (h1-fiw)*qint + fiw*qi
403!----------------the relative humidity----------------------------------
404 if(qc <= 1.0e-10) then
405 rq(n) = d00
406 else
407 rq(n) = qk / qc
408 endif
409!----------------cloud cover ratio ccr----------------------------------
411 if(rq(n) < u00k(i,k)) then
412 ccr(n) = d00
413 elseif(rq(n) >= us) then
414 ccr(n) = us
415 else
416 rqkll = min(us,rq(n))
417 ccr(n) = h1-sqrt((us-rqkll)/(us-u00k(i,k)))
418 endif
419!
420 endif
421 enddo
422!-------------------ice-water id number iwl------------------------------
423! do n=1,ihpr
424! if (comput(n) .and. (ww(n) .gt. climit)) then
425! if (tmt0(n) .lt. -15.0
426! * .or. (tmt0(n) .lt. 0.0 .and. iwl1(n) .eq. 1))
427! * iwl(n) = 1
428! cll(ipr(n),k) = 1.0 ! cloud cover!
429! cll(ipr(n),k) = min(1.0, ww(n)*cclim(k)) ! cloud cover!
430! endif
431! enddo
432!
464!--- precipitation production -- auto conversion and accretion
465!
466 do n=1,ihpr
467 if (comput(n) .and. ccr(n) > 0.0) then
468 wws = ww(n)
469 cwmk = max(cons_0, wws)
470 i = ipr(n)
471! amaxcm = max(cons_0, cwmk - wmink(n))
472 if (iwl(n) == 1) then ! ice phase
473 amaxcm = max(cons_0, cwmk - wmini(i,k))
474 expf = dt * exp(0.025*tmt0(n))
475 psaut = min(cwmk, psautco_l(i)*expf*amaxcm)
476 ww(n) = ww(n) - psaut
477 cwmk = max(cons_0, ww(n))
478! cwmk = max(cons_0, ww(n)-wmini(i,k))
479 psaci = min(cwmk, aa2*expf*precsl1(n)*cwmk)
480
481 ww(n) = ww(n) - psaci
482 precsl(n) = precsl(n) + (wws - ww(n)) * condt(n)
483 else ! liquid water
484!
493! for using sundqvist precip formulation of rain
494!
495 amaxcm = max(cons_0, cwmk - wmink(n))
496!! amaxcm = cwmk
497 tem1 = precsl1(n) + precrl1(n)
498 tem2 = min(max(cons_0, 268.0-tt(n)), cons_20)
499 tem = (1.0+c1*sqrt(tem1*rdt)) * (1+c2*sqrt(tem2))
500!
501 tem2 = amaxcm * cmr * tem / max(ccr(n),cons_p01)
502 tem2 = min(cons_50, tem2*tem2)
503! praut = c00 * tem * amaxcm * (1.0-exp(-tem2))
504 praut = (prautco_l(i)*dt) * tem * amaxcm
505 & * (1.0-exp(-tem2))
506 praut = min(praut, cwmk)
507 ww(n) = ww(n) - praut
508!
509! - Calculate the accretion of cloud water by rain \f$P_{racw}\f$,
510! can be expressed using the cloud mixing ratio \f$cwm\f$ and rainfall
511! rate \f$P_{r}\f$:
512!\f[
513! P_{racw}=C_{r}cwmP_{r}
514!\f]
515! where \f$C_{r}=5.0\times10^{-4}m^{2}kg^{-1}s^{-1}\f$ is the
516! collection coeffiecient. Note that this process is not included in
517! current operational physcics.
518! below is for zhao's precip formulation (water)
519!
520! amaxcm = max(cons_0, cwmk - wmink(n))
521! praut = min(cwmk, c00*amaxcm*amaxcm)
522! ww(n) = ww(n) - praut
523!
524! cwmk = max(cons_0, ww(n))
525! tem1 = precsl1(n) + precrl1(n)
526! pracw = min(cwmk, cr*dt*tem1*cwmk)
527! ww(n) = ww(n) - pracw
528!
529 precrl(n) = precrl(n) + (wws - ww(n)) * condt(n)
530!
531!hchuang code change [+1l] : add record to record information in vertical
532! turn rnp in unit of ww (cwm and q, kg/kg ???)
533 rnp(n) = rnp(n) + (wws - ww(n))
534 endif
535 endif
536 enddo
558!
559!-----evaporation of precipitation-------------------------
560!**** err & ers positive--->evaporation-- negtive--->condensation
561!
562 do n=1,ihpr
563 if (comput(n)) then
564 i = ipr(n)
565 qk = max(epsq, qq(n))
566 tmt0k = max(cons_m30, tmt0(n))
567 precrk = max(cons_0, precrl(n))
568 precsk = max(cons_0, precsl(n))
569 amaxrq = max(cons_0, u00k(i,k)-rq(n)) * conde(n)
570!----------------------------------------------------------------------
571! increase the evaporation for strong/light prec
572!----------------------------------------------------------------------
573 ppr = ke * amaxrq * sqrt(precrk)
574! ppr = ke * amaxrq * sqrt(precrk*rdt)
575 if (tmt0(n) .ge. 0.) then
576 pps = 0.
577 else
578 pps = (crs1+crs2*tmt0k) * amaxrq * precsk / u00k(i,k)
579 end if
580!---------------correct if over-evapo./cond. occurs--------------------
581 erk=precrk+precsk
582 if(rq(n).ge.1.0e-10) erk = amaxrq * qk * rdt / rq(n)
583 if (ppr+pps .gt. abs(erk)) then
584 rprs = erk / (precrk+precsk)
585 ppr = precrk * rprs
586 pps = precsk * rprs
587 endif
588 ppr = min(ppr, precrk)
589 pps = min(pps, precsk)
590 err(n) = ppr * rconde(n)
591 ers(n) = pps * rconde(n)
592 precrl(n) = precrl(n) - ppr
593!hchuang code change [+1l] : add record to record information in vertical
594! use err for kg/kg/dt not the ppr (mm/dt=kg/m2/dt)
595!
596 rnp(n) = rnp(n) - err(n)
597!
598 precsl(n) = precsl(n) - pps
599 endif
600 enddo
632!--------------------melting of the snow--------------------------------
633 do n=1,ihpr
634 if (comput(n)) then
635 if (tmt0(n) .gt. 0.) then
636 amaxps = max(cons_0, precsl(n))
637 psm1 = csm1 * tmt0(n) * tmt0(n) * amaxps
638 psm2 = cws * cr * max(cons_0, ww(n)) * amaxps
639 ppr = (psm1 + psm2) * conde(n)
640 if (ppr .gt. amaxps) then
641 ppr = amaxps
642 psm1 = amaxps * rconde(n)
643 endif
644 precrl(n) = precrl(n) + ppr
645!
646!hchuang code change [+1l] : add record to record information in vertical
647! turn ppr (mm/dt=kg/m2/dt) to kg/kg/dt -> ppr/air density (kg/m3)
648 rnp(n) = rnp(n) + ppr * rconde(n)
649!
650 precsl(n) = precsl(n) - ppr
651 else
652 psm1 = d00
653 endif
654!
655!---------------update t and q------------------------------------------
663
664 tt(n) = tt(n) - dtcp * (elwv*err(n)+eliv*ers(n)+eliw*psm1)
665 qq(n) = qq(n) + dt * (err(n)+ers(n))
666 endif
667 enddo
668!
669 do n=1,ihpr
670 iwl1(n) = iwl(n)
671 precrl1(n) = max(cons_0, precrl(n))
672 precsl1(n) = max(cons_0, precsl(n))
673 i = ipr(n)
674 t(i,k) = tt(n)
675 q(i,k) = qq(n)
676 cwm(i,k) = ww(n)
677 iw(i,k) = iwl(n)
678!hchuang code change [+1l] : add record to record information in vertical
679! rnp = precrl1*rconde(n) unit in kg/kg/dt
680!
681 rainp(i,k) = rnp(n)
682 enddo
683!
684! move water from vapor to liquid should the liquid amount be negative
685!
686 do i = 1, im
687 if (cwm(i,k) < 0.) then
688 tem = q(i,k) + cwm(i,k)
689 if (tem >= 0.0) then
690 q(i,k) = tem
691 t(i,k) = t(i,k) - elwv * rcp * cwm(i,k)
692 cwm(i,k) = 0.
693 elseif (q(i,k) > 0.0) then
694 cwm(i,k) = tem
695 t(i,k) = t(i,k) + elwv * rcp * q(i,k)
696 q(i,k) = 0.0
697 endif
698 endif
699 enddo
700!
701 enddo ! k loop ends here!
702!**********************************************************************
703!-----------------------end of precipitation processes-----------------
704!**********************************************************************
705!
714 do n=1,ihpr
715 i = ipr(n)
716 rn(i) = (precrl1(n) + precsl1(n)) * rrow ! precip at surface
717!
718!----sr=1 if sfc prec is rain ; ----sr=-1 if sfc prec is snow
719!----sr=0 for both of them or no sfc prec
720!
721! rid = 0.
722! sid = 0.
723! if (precrl1(n) .ge. 1.e-13) rid = 1.
724! if (precsl1(n) .ge. 1.e-13) sid = -1.
725! sr(i) = rid + sid ! sr=1 --> rain, sr=-1 -->snow, sr=0 -->both
726! chuang, june 2013: change sr to define fraction of frozen precipitation instead
727! because wpc uses it in their winter experiment
728
729 rid = precrl1(n) + precsl1(n)
730 if (rid < 1.e-13) then
731 sr(i) = 0.
732 else
733 sr(i) = precsl1(n)/rid
734 endif
735 enddo
736!
737 return
738 end subroutine zhaocarr_precpd_run
740
741 end module zhaocarr_precpd
subroutine, public zhaocarr_precpd_run(im, km, dt, del, prsl, q, cwm, t, rn, grav, hvap, hfus, ttp, cp, eps, epsm1, sr, rainp, u00k, psautco, prautco, evpco, wminco, wk1, lprnt, jpr, errmsg, errflg)
This module provides an Application Program Interface (API) for computing basic thermodynamic physics...
Definition funcphys.f90:26
This module contains the CCPP-compliant zhao_carr_precpd scheme.