CCPP SciDoc v7.0.0  v7.0.0
Common Community Physics Package Developed at DTC
 
Loading...
Searching...
No Matches
zhaocarr_gscond.f
1
5
8
9 implicit none
10 public :: zhaocarr_gscond_init, zhaocarr_gscond_run
11 private
12 logical :: is_initialized = .false.
13 contains
14
15
16! \brief Brief description of the subroutine
17!
20 subroutine zhaocarr_gscond_init (imp_physics, &
21 & imp_physics_zhao_carr, &
22 & imp_physics_zhao_carr_pdf, &
23 & errmsg, errflg)
24 implicit none
25
26 ! Interface variables
27 integer, intent(in ) :: imp_physics
28 integer, intent(in ) :: imp_physics_zhao_carr, &
29 & imp_physics_zhao_carr_pdf
30
31 ! CCPP error handling
32 character(len=*), intent( out) :: errmsg
33 integer, intent( out) :: errflg
34
35 ! Initialize the CCPP error handling variables
36 errmsg = ''
37 errflg = 0
38
39 if (is_initialized) return
40
41 ! Consistency checks
42 if (imp_physics/=imp_physics_zhao_carr .and. &
43 & imp_physics/=imp_physics_zhao_carr_pdf) then
44 write(errmsg,'(*(a))') "Logic error: namelist choice of &
45 & microphysics is different from Zhao-Carr MP"
46 errflg = 1
47 return
48 end if
49
50 is_initialized = .true.
51 end subroutine zhaocarr_gscond_init
52
53
59#if 0
60
63#endif
64
73 subroutine zhaocarr_gscond_run (im,km,dt,dtf,prsl,ps,q,clw1 &
74 &, clw2, cwm, t, tp, qp, psp &
75 &, psat,hvap,grav,hfus,ttp,rd,cp,eps,epsm1,rv &
76 &, tp1, qp1, psp1, u, lprnt, ipr, errmsg, errflg)
77
78!
79! ******************************************************************
80! * *
81! * subroutine for grid-scale condensation & evaporation *
82! * for the mrf model at ncep. *
83! * *
84! ******************************************************************
85! * *
86! * created by: q. zhao jan. 1995 *
87! * modified by: h.-l. pan sep. 1998 *
88! * modified by: s. moorthi aug. 1998, 1999, 2000 *
89! * *
90! * references: *
91! * *
92! ******************************************************************
93!
94 use machine , only : kind_phys
95 use funcphys , only : fpvs
96! use namelist_def, only: nsdfi,fhdfi
97 implicit none
98!
99! Interface variables
100 integer, intent(in) :: im, km, ipr
101 real(kind=kind_phys), intent(in) :: dt, dtf
102 real(kind=kind_phys), intent(in) :: prsl(:,:), ps(:)
103 real(kind=kind_phys), intent(inout) :: q(:,:), t(:,:)
104 real(kind=kind_phys), intent(in) :: clw1(:,:), clw2(:,:)
105 real(kind=kind_phys), intent(out) :: cwm(:,:)
106 real(kind=kind_phys), intent(inout), optional :: &
107 &, tp(:,:), qp(:,:), psp(:) &
108 &, tp1(:,:), qp1(:,:), psp1(:)
109 real(kind=kind_phys), intent(in) :: u(:,:)
110 logical, intent(in) :: lprnt
111 real(kind=kind_phys), intent(in) :: psat, hvap, grav, hfus &
112 &, ttp, rd, cp, eps, epsm1, rv
113!
114 character(len=*), intent(out) :: errmsg
115 integer, intent(out) :: errflg
116!
117! Local variables
118 real (kind=kind_phys) h1, d00, elwv, eliv
119 &, epsq, r, cpr, rcp
120!
121 parameter(h1=1.e0, d00=0.e0, epsq=2.e-12)
122!
123 real(kind=kind_phys), parameter :: cons_0=0.0, cons_m15=-15.0
124!
125 real (kind=kind_phys) qi(im), qint(im), ccrik, e0
126 &, cond, rdt, us, cclimit, climit
127 &, tmt0, tmt15, qik, cwmik
128 &, ai, qw, u00ik, tik, pres, pp0, fi
129 &, at, aq, ap, fiw, elv, qc, rqik
130 &, rqikk, tx1, tx2, tx3, es, qs
131 &, tsq, delq, condi, cone0, us00, ccrik1
132 &, aa, ab, ac, ad, ae, af, ag
133 &, el2orc, albycp
134! real (kind=kind_phys) vprs(im)
135 integer iw(im,km), i, k, iwik
136!
137 ! Initialize CCPP error handling variables
138 errmsg = ''
139 errflg = 0
140!
141!-----------------GFS interstitial in driver ----------------------------
142 do i = 1,im
143 do k= 1,km
144 cwm(i,k) = clw1(i,k)+clw2(i,k)
145 enddo
146 enddo
147!-----------------prepare constants for later uses-----------------
148!
149 elwv = hvap
150 eliv = hvap+hfus
151 r = rd
152 cpr = cp*r
153 rcp = h1/cp
154 el2orc = hvap*hvap / (rv*cp)
155 albycp = hvap / cp
156!
157 rdt = h1/dt
158 us = h1
159 cclimit = 1.0e-3
160 climit = 1.0e-20
161!
162 do i = 1, im
163 iw(i,km) = d00
164 enddo
165!
166! check for first time step
167!
168! if (tp(1,1) < 1.) then
169! do k = 1, km
170! do i = 1, im
171! tp(i,k) = t(i,k)
172! qp(i,k) = max(q(i,k),epsq)
173! tp1(i,k) = t(i,k)
174! qp1(i,k) = max(q(i,k),epsq)
175! enddo
176! enddo
177! do i = 1, im
178! psp(i) = ps(i)
179! psp1(i) = ps(i)
180! enddo
181! endif
182!
183!*************************************************************
186!*************************************************************
187!
188! do k = km-1,2,-1
189 do k = km,1,-1
190! vprs(:) = 0.001 * fpvs(t(:,k)) ! fpvs in pa
191!-----------------------------------------------------------------------
192!------------------qw, qi and qint--------------------------------------
193 do i = 1, im
194 tmt0 = t(i,k)-273.16
195 tmt15 = min(tmt0,cons_m15)
196 qik = max(q(i,k),epsq)
197 cwmik = max(cwm(i,k),climit)
198!
199! ai = 0.008855
200! bi = 1.0
201! if (tmt0 .lt. -20.0) then
202! ai = 0.007225
203! bi = 0.9674
204! end if
205!
206! the global qsat computation is done in pa
207 pres = prsl(i,k)
208!
209! qw = vprs(i)
210 qw = min(pres, fpvs(t(i,k)))
211!
212 qw = eps * qw / (pres + epsm1 * qw)
213 qw = max(qw,epsq)
214! qi(i) = qw *(bi+ai*min(tmt0,cons_0))
215! qint(i) = qw *(1.-0.00032*tmt15*(tmt15+15.))
216 qi(i) = qw
217 qint(i) = qw
218! if (tmt0 .le. -40.) qint(i) = qi(i)
219
236
237!-------------------ice-water id number iw------------------------------
238 if(tmt0.lt.-15.0) then
239 u00ik = u(i,k)
240 fi = qik - u00ik*qi(i)
241 if(fi > d00.or.cwmik > climit) then
242 iw(i,k) = 1
243 else
244 iw(i,k) = 0
245 end if
246 end if
247!
248 if(tmt0.ge.0.0) then
249 iw(i,k) = 0
250 end if
251!
252 if (tmt0 < 0.0 .and. tmt0 >= -15.0) then
253 iw(i,k) = 0
254 if (k < km) then
255 if (iw(i,k+1) == 1 .and. cwmik > climit) iw(i,k) = 1
256 endif
257 end if
258 enddo
260!--------------condensation and evaporation of cloud--------------------
261 do i = 1, im
274!------------------------at, aq and dp/dt-------------------------------
275 qik = max(q(i,k),epsq)
276 cwmik = max(cwm(i,k),climit)
277 iwik = iw(i,k)
278 u00ik = u(i,k)
279 tik = t(i,k)
280 pres = prsl(i,k)
281 pp0 = (pres / ps(i)) * psp(i)
282 at = (tik-tp(i,k)) * rdt
283 aq = (qik-qp(i,k)) * rdt
284 ap = (pres-pp0) * rdt
287!----------------the satuation specific humidity------------------------
288 fiw = float(iwik)
289 elv = (h1-fiw)*elwv + fiw*eliv
290 qc = (h1-fiw)*qint(i) + fiw*qi(i)
291! if (lprnt) print *,' qc=',qc,' qint=',qint(i),' qi=',qi(i)
292!----------------the relative humidity----------------------------------
293 if(qc.le.1.0e-10) then
294 rqik=d00
295 else
296 rqik = qik/qc
297 endif
298
382
383!----------------cloud cover ratio ccrik--------------------------------
384 if (rqik .lt. u00ik) then
385 ccrik = d00
386 elseif(rqik.ge.us) then
387 ccrik = us
388 else
389 rqikk = min(us,rqik)
390 ccrik = h1-sqrt((us-rqikk)/(us-u00ik))
391 endif
392!-----------correct ccr if it is too small in large cwm regions--------
393! if(ccrik.ge.0.01.and.ccrik.le.0.2.and
394! & .cwmik.ge.0.2e-3) then
395! ccrik=min(1.0,cwmik*1.0e3)
396! end if
397!----------------------------------------------------------------------
398! if no cloud exists then evaporate any existing cloud condensate
399!----------------evaporation of cloud water-----------------------------
400 e0 = d00
401 if (ccrik <= cclimit.and. cwmik > climit) then
402!
403! first iteration - increment halved
404!
405 tx1 = tik
406 tx3 = qik
407!
408 es = min(pres, fpvs(tx1))
409 qs = u00ik * eps * es / (pres + epsm1*es)
410 tsq = tx1 * tx1
411 delq = 0.5 * (qs - tx3) * tsq / (tsq + el2orc * qs)
412!
413 tx2 = delq
414 tx1 = tx1 - delq * albycp
415 tx3 = tx3 + delq
416!
417! second iteration
418!
419 es = min(pres, fpvs(tx1))
420 qs = u00ik * eps * es / (pres + epsm1*es)
421 tsq = tx1 * tx1
422 delq = (qs - tx3) * tsq / (tsq + el2orc * qs)
423!
424 tx2 = tx2 + delq
425 tx1 = tx1 - delq * albycp
426 tx3 = tx3 + delq
427!
428! third iteration
429!
430 es = min(pres, fpvs(tx1))
431 qs = u00ik * eps * es / (pres + epsm1*es)
432 tsq = tx1 * tx1
433 delq = (qs - tx3) * tsq / (tsq + el2orc * qs)
434 tx2 = tx2 + delq
435!
436 e0 = max(tx2*rdt, cons_0)
437! if (lprnt .and. i .eq. ipr .and. k .eq. 34)
438! & print *,' tx2=',tx2,' qc=',qc,' u00ik=',u00ik,' rqik=',rqik
439! &,' cwmik=',cwmik,' e0',e0
440
441! e0 = max(qc*(u00ik-rqik)*rdt, cons_0)
442 e0 = min(cwmik*rdt, e0)
443 e0 = max(cons_0,e0)
444 end if
445! if cloud cover > 0.2 condense water vapor in to cloud condensate
446!-----------the eqs. for cond. has been reorganized to reduce cpu------
447 cond = d00
448! if (ccrik .gt. 0.20 .and. qc .gt. epsq) then
449 if (ccrik .gt. cclimit .and. qc .gt. epsq) then
450 us00 = us - u00ik
451 ccrik1 = 1.0 - ccrik
452 aa = eps*elv*pres*qik
453 ab = ccrik*ccrik1*qc*us00
454 ac = ab + 0.5*cwmik
455 ad = ab * ccrik1
456 ae = cpr*tik*tik
457 af = ae * pres
458 ag = aa * elv
459 ai = cp * aa
460 cond = (ac-ad)*(af*aq-ai*at+ae*qik*ap)/(ac*(af+ag))
461!-----------check & correct if over condensation occurs-----------------
462 condi = (qik -u00ik *qc*1.0)*rdt
463 cond = min(cond, condi)
464!----------check & correct if supersatuation is too high----------------
465! qtemp=qik-max(0.,(cond-e0))*dt
466! if(qc.le.1.0e-10) then
467! rqtmp=0.0
468! else
469! rqtmp=qtemp/qc
470! end if
471! if(rqtmp.ge.1.10) then
472! cond=(qik-1.10*qc)*rdt
473! end if
474!-----------------------------------------------------------------------
475 cond = max(cond, d00)
476!-------------------update of t, q and cwm------------------------------
477 end if
478 cone0 = (cond-e0) * dt
479 cwm(i,k) = cwm(i,k) + cone0
480! if (lprnt .and. i .eq. ipr) print *,' t=',t(i,k),' cone0',cone0
481! &,' cond=',cond,' e0=',e0,' elv=',elv,' rcp=',rcp,' k=',k
482! &,' cwm=',cwm(i,k)
483 t(i,k) = t(i,k) + elv*rcp*cone0
484 q(i,k) = q(i,k) - cone0
485 enddo ! end of i-loop!
486 enddo ! end of k-loop!
487!
488!*********************************************************************
490!*********************************************************************
491!
493
494 if (dt > dtf+0.001) then ! three time level
495 do k = 1, km
496 do i = 1, im
497 tp(i,k) = tp1(i,k)
498 qp(i,k) = qp1(i,k)
499!
500 tp1(i,k) = t(i,k)
501 qp1(i,k) = max(q(i,k),epsq)
502 enddo
503 enddo
504 do i = 1, im
505 psp(i) = psp1(i)
506 psp1(i) = ps(i)
507 enddo
508 else ! two time level scheme - tp1, qp1, psp1 not used
509 do k = 1, km
510! write(0,*)' in gscond k=',k,' im=',im,' km=',km
511 do i = 1, im
512! write(0,*)' in gscond i=',i
513 tp(i,k) = t(i,k)
514 qp(i,k) = max(q(i,k),epsq)
515! qp(i,k) = q(i,k)
516 tp1(i,k) = tp(i,k)
517 qp1(i,k) = qp(i,k)
518 enddo
519 enddo
520 do i = 1, im
521 psp(i) = ps(i)
522 psp1(i) = ps(i)
523 enddo
524 endif
525!-----------------------------------------------------------------------
526 return
527 end subroutine zhaocarr_gscond_run
530 end module zhaocarr_gscond
subroutine, public zhaocarr_gscond_run(im, km, dt, dtf, prsl, ps, q, clw1, clw2, cwm, t, tp, qp, psp, psat, hvap, grav, hfus, ttp, rd, cp, eps, epsm1, rv, tp1, qp1, psp1, u, lprnt, ipr, 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_gscond scheme.