245 & prsi,p_phy,t_phy, &
247 & LOWLYR,SR,TRAIN_PHY, &
248 & F_ICE_PHY,F_RAIN_PHY,F_RIMEF_PHY, &
258 INTEGER,
INTENT(IN) :: D_SS,IMS,IME,LM,DX1
259 REAL,
INTENT(IN) :: DT,RHgrd
260 INTEGER,
INTENT(IN) :: THREADS
261 REAL,
INTENT(IN),
DIMENSION(ims:ime, lm+1):: &
263 REAL,
INTENT(IN),
DIMENSION(ims:ime, lm):: &
265 REAL,
INTENT(INOUT),
DIMENSION(ims:ime, lm):: &
267 REAL,
INTENT(INOUT),
DIMENSION(ims:ime, lm ):: &
269 REAL,
INTENT(INOUT),
DIMENSION(ims:ime, lm) :: &
270 & F_ICE_PHY,F_RAIN_PHY,F_RIMEF_PHY
271 REAL,
INTENT(OUT),
DIMENSION(ims:ime, lm) :: &
273 REAL,
INTENT(INOUT),
DIMENSION(ims:ime) :: &
275 REAL,
INTENT(OUT),
DIMENSION(ims:ime):: SR
276 REAL,
INTENT(OUT),
DIMENSION( ims:ime, lm ) :: &
279 INTEGER,
DIMENSION( ims:ime ),
INTENT(INOUT) :: LOWLYR
289 REAL,
DIMENSION(ims:ime):: APREC,PREC,ACPREC
292 REAL :: wc, RDIS, BETA6
296 INTEGER :: LSFC,I_index,J_index,L
297 INTEGER,
DIMENSION(ims:ime) :: LMH
298 REAL :: TC,QI,QRdum,QW,Fice,Frain,DUM,ASNOW,ARAIN
299 REAL,
DIMENSION(lm) :: P_col,Q_col,T_col,WC_col, &
300 RimeF_col,QI_col,QR_col,QW_col, THICK_col,DPCOL,pcond1d, &
301 pidep1d,piacw1d,piacwi1d,piacwr1d,piacr1d,picnd1d,pievp1d, &
302 pimlt1d,praut1d,pracw1d,prevp1d,pisub1d,pevap1d,DBZ_col, &
303 NR_col,NS_col,vsnow1d,vrain11d,vrain21d,vci1d,NSmICE1d, &
304 INDEXS1d,INDEXR1d,RFlag1d,RHC_col
320 ciacw=dt*0.25*pi*0.5*(1.e5)**c1
337 beta6=( (1.+3.*rdis*rdis)*(1.+4.*rdis*rdis)*(1.+5.*rdis*rdis)/ &
338 & ((1.+rdis*rdis)*(1.+2.*rdis*rdis) ) )
340 braut=dt*1.1e10*beta6/ncw
367 dpcol(k)=prsi(i,k)-prsi(i,k+1)
374 IF (qt(i,l) .LE. epsq) qt(i,l)=epsq
383 thick_col(l)=dpcol(l)*rgrav
386 q_col(l)=max(epsq, q(i,l))
387 IF (qt(i,l) .LE. epsq1)
THEN
389 IF (tc .LT. t_ice)
THEN
401IF (wc_col(l)>qtwarn .AND. p_col(l)<pwarn .AND. tc==tc)
THEN
402 WRITE(0,*)
'WARN4: >1 g/kg condensate in stratosphere; I,L,TC,P,QT=', &
403 i,l,tc,.01*p_col(l),1000.*wc_col(l)
404 qtwarn=max(wc_col(l),10.*qtwarn)
405 pwarn=min(p_col(l),0.5*pwarn)
408IF (warn5 .AND. tc/=tc)
THEN
409 WRITE(0,*)
'WARN5: NaN temperature; I,L,P=',i,l,.01*p_col(l)
414 IF (t_ice<=-100.) f_ice_phy(i,l)=0.
424 frain=f_rain_phy(i,l)
426 IF (fice .GE. 1.)
THEN
428 ELSE IF (fice .LE. 0.)
THEN
435 IF (qw.GT.0. .AND. frain.GT.0.)
THEN
436 IF (frain .GE. 1.)
THEN
444 IF (qi .LE. 0.) f_rimef_phy(i,l)=1.
445 rimef_col(l)=f_rimef_phy(i,l)
454 IF(dx1 .GE. 10000 .AND. p_col(l)<p_rhgrd_out)
THEN
469 & i_index, j_index, lsfc, &
470 & p_col, qi_col, qr_col, q_col, qw_col, rimef_col, t_col, &
471 & thick_col, wc_col,lm,pcond1d,pidep1d, &
472 & piacw1d,piacwi1d,piacwr1d,piacr1d,picnd1d,pievp1d,pimlt1d, &
473 & praut1d,pracw1d,prevp1d,pisub1d,pevap1d, dbz_col,nr_col,ns_col, &
474 & vsnow1d,vrain11d,vrain21d,vci1d,nsmice1d,indexs1d,indexr1d, &
481 train_phy(i,l)=(t_col(l)-t_phy(i,l))/dt
492 IF (qi_col(l) .LE. epsq)
THEN
494 IF (t_col(l) .LT. t_icek) f_ice_phy(i,l)=1.
497 f_ice_phy(i,l)=max( 0., min(1., qi_col(l)/wc_col(l)) )
498 f_rimef_phy(i,l)=max(1., rimef_col(l))
500 IF (qr_col(l) .LE. epsq)
THEN
503 dum=qr_col(l)/(qr_col(l)+qw_col(l))
506 refl_10cm(i,l)=dbz_col(l)
514 aprec(i)=(arain+asnow)*rrhol
515 prec(i)=prec(i)+aprec(i)
516 acprec(i)=acprec(i)+aprec(i)
517 IF(aprec(i) .LT. 1.e-8)
THEN
520 sr(i)=rrhol*asnow/aprec(i)
539 IF(f_ice_phy(i,k)>=1.)
THEN
541 ELSEIF(f_ice_phy(i,k)<=0.)
THEN
544 qs(i,k)=f_ice_phy(i,k)*wc
548 IF(qc(i,k)>0..AND.f_rain_phy(i,k)>0.)
THEN
549 IF(f_rain_phy(i,k).GE.1.)
THEN
553 qr(i,k)=f_rain_phy(i,k)*qc(i,k)
554 qc(i,k)=qc(i,k)-qr(i,k)
563 rainnc(i)=aprec(i)*1000.+rainnc(i)
564 rainncv(i)=aprec(i)*1000.
663 & I_index, J_index, LSFC, &
664 & P_col, QI_col, QR_col, Q_col, QW_col, RimeF_col, T_col, &
665 & THICK_col, WC_col ,LM,pcond1d,pidep1d, &
666 & piacw1d,piacwi1d,piacwr1d,piacr1d,picnd1d,pievp1d,pimlt1d, &
667 & praut1d,pracw1d,prevp1d,pisub1d,pevap1d, DBZ_col,NR_col,NS_col, &
668 & vsnow1d,vrain11d,vrain21d,vci1d,NSmICE1d,INDEXS1d,INDEXR1d, & !jul28
750 INTEGER,
INTENT(IN) :: LM,I_index, J_index, LSFC,DX1
751 REAL,
INTENT(IN) :: DTPH
752 REAL,
INTENT(INOUT) :: ARAIN, ASNOW
753 REAL,
DIMENSION(LM),
INTENT(INOUT) :: P_col, QI_col,QR_col &
754 & ,Q_col ,QW_col, RimeF_col, T_col, THICK_col,WC_col,pcond1d &
755 & ,pidep1d,piacw1d,piacwi1d,piacwr1d,piacr1d,picnd1d,pievp1d &
756 & ,pimlt1d,praut1d,pracw1d,prevp1d,pisub1d,pevap1d,DBZ_col,NR_col &
757 & ,NS_col,vsnow1d,vrain11d,vrain21d,vci1d,NSmICE1d,INDEXS1d &
758 & ,INDEXR1d,RFlag1d,RHC_col
802 REAL,
PARAMETER :: TOLER=5.e-7, c2=1./6., rho0=1.194, &
803 xratio=.025, zmin=0.01, dbzmin=-20.
815 REAL,
PARAMETER :: BLEND=1.
819 LOGICAL,
PARAMETER :: PRINT_diag=.false.
825 REAL :: EMAIRI, N0r, NLICE, NSmICE, NInuclei, Nrain, Nsnow, Nmix
827 LOGICAL :: CLEAR, ICE_logical, DBG_logical, RAIN_logical, &
829 INTEGER :: INDEX_MY,INDEXR,INDEXR1,INDEXR2,INDEXS,IPASS,ITDX,IXRF,&
830 & IXS,LBEF,L,INDEXRmin,INDEXS0C,IDR
833 REAL :: ABI,ABW,AIEVP,ARAINnew,ASNOWnew,BLDTRH,BUDGET, &
834 & CREVP,DELI,DELR,DELT,DELV,DELW,DENOMF, &
835 & DENOMI,DENOMW,DENOMWI,DIDEP, &
836 & DIEVP,DIFFUS,DLI,DTRHO,DUM,DUM1,DUM2,DUM3, &
837 & DWV0,DWVI,DWVR,DYNVIS,ESI,ESW,FIR,FLIMASS, &
838 & FWR,FWS,GAMMAR,GAMMAS, &
839 & PCOND,PIACR,PIACW,PIACWI,PIACWR,PICND,PIDEP,PIDEP_max, &
840 & PIEVP,PILOSS,PIMLT,PINIT,PP,PRACW,PRAUT,PREVP,PRLOSS, &
841 & QI,QInew,QLICE,QR,QRnew,QSI,QSIgrd,QSInew,QSW,QSW0, &
842 & QSWgrd,QSWnew,QT,QTICE,QTnew,QTRAIN,Q,QW,QWnew,Rcw, &
843 & RFACTOR,RFmx,RFrime,RHO,RIMEF,RIMEF1,RQR,RR,RRHO,SFACTOR, &
844 & TC,TCC,TFACTOR,THERM_COND,THICK,TK,TK2,TNEW, &
845 & TOT_ICE,TOT_ICEnew,TOT_RAIN,TOT_RAINnew, &
846 & VEL_INC,VENTR,VENTIL,VENTIS,VRAIN1,VRAIN2,VRIMEF,VSNOW, &
847 & VSNOW1,WC,WCnew,WSgrd,WS,WSnew,WV,WVnew, &
849 & NSImax,QRdum,QSmICE,QLgIce,RQLICE,VCI,TIMLT, &
850 & RQSnew,RQRnew,Zrain,Zsnow,Ztot,RHOX0C,RFnew,PSDEP,DELS
851 REAL,
SAVE :: Revised_LICE=1.e-3
872big_loop:
DO l=lm,1,-1
904 IF (arain .GT. climit)
THEN
917 IF (asnow .GT. climit)
THEN
942 esw=min(1000.*fpvs0(tk),0.99*pp)
947 esi=min(1000.*fpvs(tk),0.99*pp)
960 IF (wv.GT.wsgrd .OR. wc.GT.epsq) clear=.false.
982 rho=pp/(rd*tk*(1.+eps1*q))
1023 dum=xmimax*exp(xmiexp*dum1)
1024 indexs=min(mdimax, max(mdimin, int(dum) ) )
1059 tfactor=sqrt(tk*tk*tk)/(tk+120.)
1060 dynvis=1.496e-6*tfactor
1061 therm_cond=2.116e-3*tfactor
1062 diffus=8.794e-5*tk**1.81/pp
1067 gammas=min(1.5, (1.e5/pp)**c1)
1071 gammar=(rho0/rho)**.4
1079 IF (tc.LT.0. .OR. qi.GT. epsq .OR. asnow.GT.climit)
THEN
1086 IF (t_ice <= -100.)
THEN
1094 rain_logical=.false.
1095 IF (arain.GT.climit .OR. qr.GT.epsq) rain_logical=.true.
1097ice_test:
IF (ice_logical)
THEN
1137 nsimax=max(nsi_max, 0.1*rho*qi/massi(mdimin) )
1148 ninuclei=min(0.01*exp(-0.6*tc), nsimax)
1149 ninuclei=min(5.*exp(-0.304*tc), ninuclei)
1151 dum=rrho*massi(mdimin)
1152 nsmice=min(ninuclei, qi/dum)
1156 init_ice:
IF (qi<=epsq .AND. asnow<=climit)
THEN
1166 xlimass=rimef1*massi(indexs)
1181 tot_ice=thick*qi+blend*asnow
1182 piloss=-tot_ice/thick
1183 qlgice=max(0., qi-qsmice)
1184 vci=gammas*vsnowi(mdimin)
1189 rimef1=(rimef_col(l)*thick*qlgice &
1190 & +rimef_col(lbef)*blend*asnow)/tot_ice
1194 dum3=rh_ngc*(rho*qlgice)**c1
1195 dum2=rh_ngt*(rho*qlgice)**c1
1196 IF (rimef1>=10.)
THEN
1198 ELSE IF (rimef1>=5.)
THEN
1200 dum=dum3*dum+dum2*(1.-dum)
1203 dum=0.33333*(rimef1-2.)
1204 dum=dum2*dum+dum1*(1.-dum)
1206 indexs=min(mdimax, max(mdimin, int(dum) ) )
1209 emairi=thick+bldtrh*vsnow1
1210 qlice=(thick*qlgice+blend*asnow)/emairi
1220 two_pass:
DO ipass=1,2
1226 dum=1.e-6*real(indexs)
1227 rfmx=rfmx1*dum*dum*dum/massi(indexs)
1228 rimef1=min(rimef1, rfmx)
1230 vel_rime:
IF (rimef1<=1.)
THEN
1235 rimef1=min(rimef1, rfmax)
1236 ixs=max(2, min(indexs/100, 9))
1237 xrf=10.492*alog(rimef1)
1238 ixrf=max(0, min(int(xrf), nrime))
1239 IF (ixrf .GE. nrime)
THEN
1240 vrimef=vel_rf(ixs,nrime)
1242 vrimef=vel_rf(ixs,ixrf)+(xrf-float(ixrf))* &
1243 & (vel_rf(ixs,ixrf+1)-vel_rf(ixs,ixrf))
1245 vrimef=max(1., vrimef)
1249 vel_inc=gammas*vrimef
1250 vsnow=vel_inc*vsnowi(indexs)
1251 IF (rimef1>=5. .AND. indexs==mdimax .AND. rqlice>rqhail)
THEN
1253 dum=gammas*avhail*rqlice**bvhail
1256 vel_inc=vsnow/vsnowi(indexs)
1259 xlimass=rimef1*massi(indexs)
1260 nlice=rqlice/xlimass
1266 (nlice>=nlimin .AND. nlice<=nsi_max))
EXIT two_pass
1273 nlice=max(nlimin, min(nsi_max, nlice) )
1276 xli=rqlice/(nlice*rimef1)
1277new_size:
IF (xli<=massi(mdimin) )
THEN
1279 ELSE IF (xli<=massi(450) )
THEN new_size
1280 dli=9.5885e5*xli**.42066
1281 indexs=min(mdimax, max(mdimin, int(dli) ) )
1282 ELSE IF (xli<massi(mdimax) )
THEN new_size
1283 dli=3.9751e6*xli**.49870
1284 indexs=min(mdimax, max(mdimin, int(dli) ) )
1306 rhox0c=22.5*max(1.,min(rimef1,40.))
1308 IF(.NOT.rain_logical)
THEN
1311 IF(.NOT.ice_logical) timlt=0.
1314 IF(timlt>epsq .AND. rhox0c<=225.)
THEN
1320 IF(strat .AND. indexrmin<=mdrmin)
THEN
1321 indexrmin=indexs0c*(0.001*rhox0c)**c1
1322 indexrmin=max(mdrmin, min(indexrmin, indexrstrmax) )
1327 IF(strat .OR. timlt>epsq)
THEN
1342 IF (qw.GT.epsq .AND. tc.GE.t_ice)
THEN
1345 IF (dum2>qaut0)
THEN
1350 praut=min(qw, dum*(1.-exp(-dum1*dum1)) )
1352 IF (qlice .GT. epsq)
THEN
1357 fws=min(1., ciacw*vel_inc*nlice*accri(indexs) )
1370ice_only:
IF (tc.LT.t_ice .AND. (wv.GT.qswgrd .OR. qw.GT.epsq))
THEN
1381 dum=min(1000.*fpvs(dum1),0.99*pp)
1382 dum=rhgrd*eps*dum/(pp-dum)
1385IF (warn1 .AND. dum1<xmin)
THEN
1386 WRITE(0,*)
'WARN1: Water saturation T<180K; I,J,L,TC,P=', &
1387 i_index,j_index,l,dum1-t0c,.01*pp
1393 & pidep=
deposit(pp,dum1,dum2,rhgrd,i_index,j_index,l)
1397 ELSE IF (tc<0.)
THEN ice_only
1402 denomi=1.+xls3*qsi*tk2
1405 pidep_max=max(piloss, dwvi/denomi)
1407 IF (qtice .GT. 0.)
THEN
1416 abi=1./(rho*xls2*qsi*tk2/therm_cond+1./diffus)
1424 dum=(rho/(diffus*diffus*dynvis))**c2
1425 sfactor=sqrt(gammas)*dum
1426 ventis=(venti1(mdimin)+sfactor*venti2(mdimin))*nsmice
1427 sfactor=sqrt(vel_inc)*dum
1428 ventil=(venti1(indexs)+sfactor*venti2(indexs))*nlice
1429 didep=abi*(ventil+ventis)*dtph
1450 IF (wv>qswgrd .AND. tc<t_ice_init .AND. nsmice<ninuclei)
THEN
1451 index_my=max(my_t1, min( int(.5-tc), my_t2 ) )
1453 dum=max(0., ninuclei-nsmice)
1454 pinit=max(0., dum*my_growth_nmm(index_my)*rrho)
1460 pidep=min(dwvi*didep+pinit, pidep_max)
1461 ELSE IF (dwvi<0.)
THEN
1462 pidep=max(dwvi*didep, pidep_max)
1470 IF (tc>=t_ice .AND. (qw>epsq .OR. wv>qswgrd))
THEN
1473 dum2=tk+xls1*pidep+xlf1*piacwi
1475 pcond=
condense(pp,dum1,dum2,dum3,rhgrd,i_index,j_index,l)
1482 tcc=tc+xlv1*pcond+xls1*pidep+xlf1*piacwi
1492 IF (tc.GT.0. .AND. tcc.GT.0. .AND. ice_logical)
THEN
1499 sfactor=sqrt(vel_inc)*(rho/(diffus*diffus*dynvis))**c2
1500 ventil=nlice*(venti1(indexs)+sfactor*venti2(indexs))
1501 aievp=ventil*diffus*dtph
1502 IF (aievp .LT. xratio)
THEN
1505 dievp=1.-exp(-aievp)
1508 IF (wv.LT.qsw .AND. dum.LE.epsq)
THEN
1513 dum=min(esw0, 0.99*pp)
1514 qsw0=max(epsq, eps*dum/(pp-dum) )
1515 dwv0=min(wv,qsw)-qsw0
1517 pievp=max( min(0., dum), piloss)
1520 pimlt=therm_cond*tcc*ventil*rrho*dtph/xlf
1524 dum1=max( 0., (tcc+xlv1*pievp)/xlf1 )
1525 pimlt=min(pimlt, dum1)
1530 IF (dum .LT. piloss)
THEN
1555do_rain:
IF (rain_logical)
THEN
1556 tot_rain=thick*qr+blend*arain
1557 prloss=-tot_rain/thick
1559 qtrain=tot_rain/(thick+bldtrh*vrain1)
1561 IF(strat .AND. rqr>1.e-3) strat=.false.
1562 IF(drzl .AND. pimlt>epsq) drzl=.false.
1564dsd1:
IF (rqr<=rqr_drmin)
THEN
1567 n0r=rqr/massr(indexr)
1568 ELSE IF (drzl)
THEN dsd1
1570 dum=max(1.0, 0.5e-3/rqr)
1571 n0r=min(1.e9, n0r0*dum*sqrt(dum) )
1572 dum1=rqr/(pi*rhol*n0r)
1573 dum=1.e6*sqrt(sqrt(dum1))
1574 indexr=max(xmrmin, min(dum, xmrmax) )
1575 n0r=rqr/massr(indexr)
1577 ELSE IF (rqr>=rqr_drmax)
THEN dsd1
1580 n0r=rqr/massr(indexr)
1584 dum=cn0r0*sqrt(sqrt(rqr))
1585 indexr=max(xmrmin, min(dum, xmrmax) )
1586 IF (strat .AND. indexr<indexrmin)
THEN
1589 n0r=rqr/massr(indexr)
1593 vrain2=gammar*vrain(indexr)
1596 IF (tc .LT. t_ice)
THEN
1601 IF (dwvr.LT.0. .AND. dum.LE.epsq)
THEN
1612 rfactor=sqrt(gammar)*(rho/(diffus*diffus*dynvis))**c2
1613 abw=1./(rho*xlv2*qsw*tk2/therm_cond+1./diffus)
1618 ventr=n0r*(ventr1(indexr)+rfactor*ventr2(indexr))
1619 crevp=abw*ventr*dtph
1620 prevp=max(dwvr*crevp, prloss)
1621 ELSE IF (qw .GT. epsq)
THEN
1622 fwr=cracw*gammar*n0r*accrr(indexr)
1623 pracw=min(1.,fwr)*qw
1626 IF (ice_logical .AND. tc<0. .AND. tcc<0.)
THEN
1631 dum=.001*float(indexr)
1632 dum=(exp(abfr*tc)-1.)*dum*dum*dum*dum*dum*dum*dum
1633 piacr=min(cbfr*n0r*rrho*dum, qtrain)
1637 IF (qlice .GT. epsq)
THEN
1641 dum=gammar*vrain(indexr)
1648 dum2=sqrt(dum1*dum1+.04*dum*vsnow)
1649 dum1=5.e-12*indexr*indexr+2.e-12*indexr*indexs &
1650 & +.5e-12*indexs*indexs
1651 fir=min(1., ciacr*nlice*dum1*dum2)
1655 piacr=min(piacr+fir*qtrain, qtrain)
1658 If (dum .LT. prloss)
THEN
1680 dum1=piacw+praut+pracw-min(0.,pcond)
1681 IF (dum1 .GT. qw)
THEN
1687 IF (pcond .LT. 0.) pcond=dum*pcond
1693 delw=pcond-piacw-praut-pracw
1695 IF (qwnew .LE. epsq) qwnew=0.
1696 IF (qw.GT.0. .AND. qwnew.NE.0.)
THEN
1698 IF (dum .LT. toler) qwnew=0.
1703 delt= xlv1*(pcond+pievp+picnd+prevp) &
1704 & +xls1*pidep+xlf1*(piacwi+piacr-pimlt)
1707 delv=-pcond-pidep-pievp-picnd-prevp
1734 IF (ice_logical)
THEN
1735 deli=pidep+pievp+piacwi+piacr-pimlt
1736 tot_icenew=tot_ice+thick*deli
1737 IF (tot_ice.GT.0. .AND. tot_icenew.NE.0.)
THEN
1738 dum=tot_icenew/tot_ice
1739 IF (dum .LT. toler) tot_icenew=0.
1741 IF (tot_icenew .LE. climit)
THEN
1751 IF (pinit<=epsq)
THEN
1752 psdep=max(0., pidep)
1757 dels=piacwi+piacr+psdep
1758 IF (dels<=epsq .OR. tnew>=t0c)
THEN
1809 dum=1.e-6*real(indexs)
1810 dum1=dum*dum*dum/massi(indexs)
1812 IF (piacwi>0. .AND. rcw>0.)
THEN
1813 dum=-rcw*vsnow/min(-0.5,tc)
1814 dum=min(12.14, max(0.275, dum) )
1816 dum=min(900., max(170., dum) )
1823 rfnew=(psdep+rfrime*piacwi+rfmx*piacr)/dels
1824 dum2=tot_ice+thick*dels
1825 dum3=rimef1*tot_ice+rfnew*thick*dels
1826 IF (dum2<=climit)
THEN
1829 rimef=min(rfmx, max(1., dum3/dum2) )
1834 dum1=bldtrh*(flimass*vsnow+(1.-flimass)*vci)
1835 qinew=tot_icenew/(thick+dum1)
1836 IF (qinew .LE. epsq) qinew=0.
1837 IF (qi.GT.0. .AND. qinew.NE.0.)
THEN
1839 IF (dum .LT. toler) qinew=0.
1842 IF (asnow.GT.0. .AND. asnownew.NE.0.)
THEN
1844 IF (dum .LT. toler) asnownew=0.
1846 rqsnew=flimass*rho*qinew
1849 nsnow=rqsnew/(rimef*massi(indexs))
1850 IF (rqsnew>0.0025 .AND. rimef>5.)
THEN
1852 IF (rqsnew<=0.005)
THEN
1853 dum=min(1., 400.*rqsnew-1.)
1854 nsnow=1.e3*dum+nsnow*(1.-dum)
1860 dum=180.*(rqsnew-0.005)
1861 nsnow=1.e3*(1.-min(dum,0.8))
1864 zsnow=cdry*rqsnew*rqsnew/nsnow
1882 timlt=timlt+pimlt*thick
1883 delr=praut+pracw+piacwr-piacr+pimlt+prevp+picnd
1884 tot_rainnew=tot_rain+thick*delr
1885 IF (tot_rain.GT.0. .AND. tot_rainnew.NE.0.)
THEN
1886 dum=tot_rainnew/tot_rain
1887 IF (dum .LT. toler) tot_rainnew=0.
1889update_rn:
IF (tot_rainnew .LE. climit)
THEN
1896 qrnew=tot_rainnew/(thick+bldtrh*vrain2)
1900 IF(strat .AND. rqrnew>1.e-3) strat=.false.
1901 IF(drzl .AND. timlt>epsq) drzl=.false.
1906rain_pass:
DO ipass=1,3
1907dsd2:
IF (rqrnew<=rqr_drmin)
THEN
1910 n0r=rqrnew/massr(indexr)
1911 ELSE IF (drzl)
THEN dsd2
1913 dum=max(1.0, 0.5e-3/rqrnew)
1914 n0r=min(1.e9, n0r0*dum*sqrt(dum) )
1915 dum1=rqrnew/(pi*rhol*n0r)
1916 dum=1.e6*sqrt(sqrt(dum1))
1917 indexr=max(xmrmin, min(dum, xmrmax) )
1918 n0r=rqrnew/massr(indexr)
1919 ELSE IF (rqrnew>=rqr_drmax)
THEN
1922 n0r=rqrnew/massr(indexr)
1926 dum=cn0r0*sqrt(sqrt(rqrnew))
1927 indexr=max(xmrmin, min(dum, xmrmax) )
1928 IF (strat .AND. indexr<indexrmin)
THEN
1931 n0r=rqrnew/massr(indexr)
1934 vrain2=gammar*vrain(indexr)
1935 qrnew=tot_rainnew/(thick+bldtrh*vrain2)
1937 idr=abs(indexr-indexr2)
1939 IF(idr<=5)
EXIT rain_pass
1943 IF (qrnew .LE. epsq) qrnew=0.
1944 IF (qr.GT.0. .AND. qrnew.NE.0.)
THEN
1946 IF (dum .LT. toler) qrnew=0.
1948 arainnew=bldtrh*vrain2*qrnew
1949 IF (arain.GT.0. .AND. arainnew.NE.0.)
THEN
1951 IF (dum .LT. toler) arainnew=0.
1955 IF (rqrnew>epsq)
THEN
1957 nrain=n0r*1.e-6*real(indexr2)
1958 zrain=crain*rqrnew*rqrnew/nrain
1963 wcnew=qwnew+qrnew+qinew
1971 qt=thick*(wv+wc)+arain+asnow
1972 qtnew=thick*(wvnew+wcnew)+arainnew+asnownew
1979 IF (dum .GT. toler)
THEN
1980 dum=dum/min(qt, qtnew)
1981 IF (dum .GT. toler) dbg_logical=.true.
1984 IF (print_diag)
THEN
1985 esw=min(1000.*fpvs0(tnew),0.99*pp)
1986 qswnew=(rhgrd+0.001)*eps*esw/(pp-esw)
1987 IF( (qwnew>epsq .OR. qrnew>epsq .OR. wvnew>qswnew) &
1988 & .AND. tc<t_ice) dbg_logical=.true.
1991 IF ((wvnew.LT.epsq .OR. dbg_logical) .AND. print_diag)
THEN
1993 WRITE(0,
"(/2(a,i4),2(a,i2))")
'{} i=',i_index,
' j=',j_index,&
1994 &
' L=',l,
' LSFC=',lsfc
1996 esw=min(1000.*fpvs0(tnew),0.99*pp)
1997 qswnew=eps*esw/(pp-esw)
1998 IF (tc.LT.0. .OR. tnew .LT. 0.)
THEN
1999 esi=min(1000.*fpvs(tnew),0.99*pp)
2000 qsinew=eps*esi/(pp-esi)
2006 WRITE(0,
"(4(a12,g11.4,1x))") &
2007 &
'{} TCold=',tc,
'TCnew=',tnew-t0c,
'P=',.01*pp,
'RHO=',rho, &
2008 &
'{} THICK=',thick,
'RHold=',wv/ws,
'RHnew=',wvnew/wsnew, &
2010 &
'{} RHWold=',wv/qsw,
'RHWnew=',wvnew/qswnew,
'RHIold=',wv/qsi, &
2011 &
'RHInew=',wvnew/qsinew, &
2012 &
'{} QSWold=',qsw,
'QSWnew=',qswnew,
'QSIold=',qsi,
'QSInew=',qsinew, &
2013 &
'{} WSold=',ws,
'WSnew=',wsnew,
'WVold=',wv,
'WVnew=',wvnew, &
2014 &
'{} WCold=',wc,
'WCnew=',wcnew,
'QWold=',qw,
'QWnew=',qwnew, &
2015 &
'{} QIold=',qi,
'QInew=',qinew,
'QRold=',qr,
'QRnew=',qrnew, &
2016 &
'{} ARAINold=',arain,
'ARAINnew=',arainnew,
'ASNOWold=',asnow, &
2017 &
'ASNOWnew=',asnownew, &
2018 &
'{} TOT_RAIN=',tot_rain,
'TOT_RAINnew=',tot_rainnew, &
2019 &
'TOT_ICE=',tot_ice,
'TOT_ICEnew=',tot_icenew, &
2020 &
'{} BUDGET=',budget,
'QTold=',qt,
'QTnew=',qtnew
2022 WRITE(0,
"(4(a12,g11.4,1x))") &
2023 &
'{} DELT=',delt,
'DELV=',delv,
'DELW=',delw,
'DELI=',deli, &
2024 &
'{} DELR=',delr,
'PCOND=',pcond,
'PIDEP=',pidep,
'PIEVP=',pievp, &
2025 &
'{} PICND=',picnd,
'PREVP=',prevp,
'PRAUT=',praut,
'PRACW=',pracw, &
2026 &
'{} PIACW=',piacw,
'PIACWI=',piacwi,
'PIACWR=',piacwr,
'PIMLT=', &
2030 WRITE(0,
"(4(a15,L2))") &
2031 &
'{} ICE_logical=',ice_logical,
'RAIN_logical=',rain_logical, &
2032 &
'STRAT=',strat,
'DRZL=',drzl
2034 IF (ice_logical)
WRITE(0,
"(4(a12,g11.4,1x))") &
2035 &
'{} RimeF1=',rimef1,
'GAMMAS=',gammas,
'VrimeF=',vrimef, &
2037 &
'{} QSmICE=',qsmice,
'XLIMASS=',xlimass,
'RQLICE=',rqlice, &
2039 &
'{} NLICE=',nlice,
'NSmICE=',nsmice,
'PILOSS=',piloss, &
2040 &
'EMAIRI=',emairi, &
2041 &
'{} RimeF=',rimef,
'VCI=',vci,
'INDEXS=',float(indexs), &
2042 &
'FLIMASS=',flimass, &
2043 &
'{} Nsnow=',nsnow,
'Zsnow=',zsnow,
'TIMLT=',timlt,
'RHOX0C=',rhox0c, &
2044 &
'{} INDEXS0C=',float(indexs0c)
2046 IF (rain_logical)
WRITE(0,
"(4(a12,g11.4,1x))") &
2047 &
'{} INDEXR1=',float(indexr1),
'INDEXR=',float(indexr), &
2048 &
'GAMMAR=',gammar,
'N0r=',n0r, &
2049 &
'{} VRAIN1=',vrain1,
'VRAIN2=',vrain2,
'QTRAIN=',qtrain,
'RQR=',rqr, &
2050 &
'{} PRLOSS=',prloss,
'VOLR1=',thick+bldtrh*vrain1, &
2051 &
'VOLR2=',thick+bldtrh*vrain2,
'INDEXR2=',indexr2, &
2052 &
'{} Nrain=',nrain,
'Zrain=',zrain
2054 IF (pracw .GT. 0.)
WRITE(0,
"(a12,g11.4,1x)")
'{} FWR=',fwr
2056 IF (piacr .GT. 0.)
WRITE(0,
"(a12,g11.4,1x)")
'{} FIR=',fir
2058 dum=pimlt+picnd-prevp-pievp
2059 IF (dum.GT.0. .or. dwvi.NE.0.) &
2060 &
WRITE(0,
"(4(a12,g11.4,1x))") &
2061 &
'{} TFACTOR=',tfactor,
'DYNVIS=',dynvis, &
2062 &
'THERM_CON=',therm_cond,
'DIFFUS=',diffus
2064 IF (prevp .LT. 0.)
WRITE(0,
"(4(a12,g11.4,1x))") &
2065 &
'{} RFACTOR=',rfactor,
'ABW=',abw,
'VENTR=',ventr,
'CREVP=',crevp, &
2066 &
'{} DWVr=',dwvr,
'DENOMW=',denomw
2068 IF (pidep.NE.0. .AND. dwvi.NE.0.) &
2069 &
WRITE(0,
"(4(a12,g11.4,1x))") &
2070 &
'{} DWVi=',dwvi,
'DENOMI=',denomi,
'PIDEP_max=',pidep_max, &
2071 &
'SFACTOR=',sfactor, &
2072 &
'{} ABI=',abi,
'VENTIL=',ventil,
'VENTIL1=',venti1(indexs), &
2073 &
'VENTIL2=',sfactor*venti2(indexs), &
2074 &
'{} VENTIS=',ventis,
'DIDEP=',didep
2076 IF (pidep.GT.0. .AND. pcond.NE.0.) &
2077 &
WRITE(0,
"(4(a12,g11.4,1x))") &
2078 &
'{} DENOMW=',denomw,
'DENOMWI=',denomwi,
'DENOMF=',denomf, &
2079 &
'DUM2=',pcond-piacw
2081 IF (fws .GT. 0.)
WRITE(0,
"(4(a12,g11.4,1x))") &
2084 dum=pimlt+picnd-pievp
2085 IF (dum.GT. 0.)
WRITE(0,
"(4(a12,g11.4,1x))") &
2086 &
'{} SFACTOR=',sfactor,
'VENTIL=',ventil,
'VENTIL1=',venti1(indexs), &
2087 &
'VENTIL2=',sfactor*venti2(indexs), &
2088 &
'{} AIEVP=',aievp,
'DIEVP=',dievp,
'QSW0=',qsw0,
'DWV0=',dwv0
2097 q_col(l)=max(epsq, wvnew/(1.+wvnew))
2098 wc_col(l)=max(epsq, wcnew)
2099 qi_col(l)=max(epsq, qinew)
2100 qr_col(l)=max(epsq, qrnew)
2101 qw_col(l)=max(epsq, qwnew)
2108 ztot=max(ztot, zrain+zsnow)
2110 IF (ztot>zmin) dbz_col(l)=10.*alog10(ztot)
2120 elseif(pcond.lt.0)
then
2125 elseif(pidep.lt.0)
then
2138 if (qinew>epsq)
then
2141 if (flimass<1.)
then
2146 indexs1d(l)=float(indexs)
2148 if (qrnew>epsq)
then
2151 indexr1d(l)=float(indexr2)
2186 REAL (kind=kind_phys),
PARAMETER :: &
2187 & rhlimit=.001, rhlimit1=-rhlimit
2188 REAL (kind=kind_phys) :: cond, ssat, wcdum
2190 REAL,
INTENT(IN) :: qw,pp,wv,tk,rhgrd
2191 REAL wvdum,tdum,dwv,ws,esw
2193integer,
INTENT(IN) :: i,j,l
2195real :: dwvinp,ssatinp
2203 esw=min(1000.*fpvs0(tdum),0.99*pp)
2204 ws=rhgrd*eps*esw/(pp-esw)
2212 DO niter=1,max_iterations
2213 cond=dwv/(1.+xlv3*ws/(tdum*tdum))
2214 cond=max(cond, -wcdum)
2219 esw=min(1000.*fpvs0(tdum),0.99*pp)
2220 ws=rhgrd*eps*esw/(pp-esw)
2223 IF (ssat>=rhlimit1 .AND. ssat<=rhlimit)
EXIT
2224 IF (ssat<rhlimit1 .AND. wcdum<=epsq)
EXIT
2228IF (niter>max_iterations)
THEN
2231 write(0,*)
'WARN3: Too many iterations in function CONDENSE; ', &
2232 ' I,J,L,TC,SSAT,QW,DWV=',i,j,l,tk-t0c,ssatinp,1000.*qw,dwvinp
2252 use machine,
only: high_pres => kind_dbl_prec
2256 REAL (kind=high_pres),
PARAMETER :: rhlimit=.001, &
2258 REAL (kind=high_pres) :: dep, ssat
2260 real,
INTENT(IN) :: pp,rhgrd
2261 real,
INTENT(INOUT) :: wvdum,tdum
2264integer,
INTENT(IN) :: i,j,l
2266real :: tinp,dwvinp,ssatinp
2271 esi=min(1000.*fpvs(tdum),0.99*pp)
2272 ws=rhgrd*eps*esi/(pp-esi)
2281 DO niter=1,max_iterations
2282 dep=dwv/(1.+xls3*ws/(tdum*tdum))
2286 esi=min(1000.*fpvs(tdum),0.99*pp)
2287 ws=rhgrd*eps*esi/(pp-esi)
2290 IF (ssat>=rhlimit1 .AND. ssat<=rhlimit)
EXIT
2294IF (niter>max_iterations)
THEN
2297 write(0,*)
'WARN2: Too many iterations in function DEPOSIT; ', &
2298 ' I,J,L,TC,SSAT,DWV=',i,j,l,tinp-t0c,ssatinp,dwvinp
2314 INTEGER FUNCTION get_indexr(RR)
2316 REAL,
INTENT(IN) :: rr
2317 IF (rr .LE. rr_drmin)
THEN
2323 ELSE IF (rr .LE. rr_dr1)
THEN
2332 get_indexr=int( 1.123e3*rr**.1947 + .5 )
2333 get_indexr=max( mdrmin, min(get_indexr, mdr1) )
2334 ELSE IF (rr .LE. rr_dr2)
THEN
2343 get_indexr=int( 1.225e3*rr**.2017 + .5 )
2344 get_indexr=max( mdr1, min(get_indexr, mdr2) )
2345 ELSE IF (rr .LE. rr_dr3)
THEN
2354 get_indexr=int( 1.3006e3*rr**.2083 + .5 )
2355 get_indexr=max( mdr2, min(get_indexr, mdr3) )
2356 ELSE IF (rr .LE. rr_dr4)
THEN
2365 get_indexr=int( 1.354e3*rr**.2143 + .5 )
2366 get_indexr=max( mdr3, min(get_indexr, mdr4) )
2367 ELSE IF (rr .LE. rr_dr5)
THEN
2376 get_indexr=int( 1.404e3*rr**.2213 + .5 )
2377 get_indexr=max( mdr4, min(get_indexr, mdr5) )
2378 ELSE IF (rr .LE. rr_drmax)
THEN
2387 get_indexr=int( 1.4457e3*rr**.2303 + .5 )
2388 get_indexr=max( mdr5, min(get_indexr, mdrmax) )
2397 END FUNCTION get_indexr