285 &, ccwf, area, dxmin, dxinv &
286 &, psauras, prauras, wminras, dlqf, flipv &
287 &, me, rannum, nrcm, mp_phys, mp_phys_mg &
289 &, tin, qin, uin, vin, ccin, fscav &
290 &, prsi, prsl, prsik, prslk, phil, phii &
291 &, KPBL, CDRAG, RAINC, kbot, ktop, kcnv &
292 &, DDVEL, ud_mf, dd_mf, dt_mf &
293 &, QLCN, QICN, w_upi, cf_upi, CNV_MFD &
294 &, CNV_DQLDT,CLCN,CNV_FICE,CNV_NDROP,CNV_NICE &
317 logical,
intent(in) :: flipv
319 integer,
intent(in) :: im, k, itc, ntc, ntr, me, nrcm, ntk, kdt &
320 &, mp_phys, mp_phys_mg
321 integer,
dimension(:),
intent(out) :: kbot, ktop
322 integer,
dimension(:),
intent(inout) :: kcnv
323 integer,
dimension(:),
intent(in) :: kpbl
325 real(kind=kind_phys),
intent(in) :: dxmin, dxinv, ccwf(:) &
326 &, psauras(:), prauras(:) &
327 &, wminras(:), dlqf(:)
329 real(kind=kind_phys),
dimension(:,:),
intent(in) :: prsi, prsik, phii
331 real(kind=kind_phys),
dimension(:,:),
intent(inout) :: tin, qin, uin, vin
332 real(kind=kind_phys),
dimension(:,:),
intent(in) :: prsl, prslk, phil &
334 real(kind=kind_phys),
dimension(:,:),
intent(out),
optional :: ud_mf
335 real(kind=kind_phys),
dimension(:,:),
intent(out) :: dd_mf, dt_mf
336 real(kind=kind_phys),
dimension(:,:),
intent(inout),
optional :: qlcn, qicn, w_upi &
339 &, cnv_fice, cnv_ndrop &
341 real(kind=kind_phys),
dimension(:) ,
intent(in) :: area, cdrag
342 real(kind=kind_phys),
dimension(:) ,
intent(out) :: rainc
343 real(kind=kind_phys),
dimension(:) ,
intent(out),
optional :: ddvel
344 real(kind=kind_phys),
dimension(:,:),
intent(in) :: rannum
345 real(kind=kind_phys),
intent(inout) :: ccin(:,:,:)
346 real(kind=kind_phys),
intent(in) :: dt, dtf
350 real(kind=kind_phys),
intent(in) :: fscav(:)
353 character(len=*),
intent(out) :: errmsg
354 integer,
intent(out) :: errflg
358 real(kind=kind_phys) :: trcmin(ntr+2)
359 real(kind=kind_phys),
dimension(k) :: toi, qoi, tcu, qcu &
360 &, pcu, clw, cli, qii, qli&
361 &, phi_l, prsm,psjm &
364 &, qoi_l, qli_l, qii_l
365 real(kind=kind_phys),
dimension(k+1) :: prs, psj, phi_h, flx, flxd
368 integer,
dimension(100) :: ic
369 real(kind=kind_phys),
parameter :: clwmin=1.0e-10_kp
371 real(kind=kind_phys),
allocatable :: alfint(:,:), uvi(:,:) &
372 &, trcfac(:,:), rcu(:,:)
373 real(kind=kind_phys) dtvd(2,4)
375 real(kind=kind_phys) cfac, tem, sgc, ccwfac, tem1, tem2, rain &
376 &, wfnc,tla,pl,qiid,qlid, c0, c0i, dlq_fac, sumq&
380 Integer kcr, kfx, ncmx, nc, ktem, i, ii, lm1, l &
381 &, ntrc, ll, km1, kp1, ipt, lv, KBL, n &
382 &, KRMIN, KRMAX, KFMAX, kblmx, irnd,ib &
383 &, kblmn, ksfc, ncrnd
384 real(kind=kind_phys) sgcs(k)
391 if (itc > 0 .and. ntc > 0)
then
393 if (n <= ntr + 2)
then
394 fscav_(itc:n) = fscav
396 errmsg =
'Error in rascnv_run: test ntr >= itc + ntc - 3 FAILED'
402 if (ntk-2 > 0) trcmin(ntk-2) = 1.0e-4_kp
422 if (.not.
allocated(trcfac))
allocate (trcfac(k,ntrc))
423 if (.not.
allocated(uvi))
allocate (uvi(k,ntrc))
424 if (.not.
allocated(rcu))
allocate (rcu(k,ntrc))
434 if(mp_phys == mp_phys_mg)
then
443 cnv_dqldt(i,l) = zero
446 cnv_ndrop(i,l) = zero
452 if (.not.
allocated(alfint))
allocate(alfint(k,ntrc+4))
467 tem1 = max(zero, min(one, (log(area(ipt)) - dxmin) * dxinv))
469 ccwfac = ccwf(1)*tem1 + ccwf(2)*tem2
470 dlq_fac = dlqf(1)*tem1 + dlqf(2)*tem2
472 c0i = (psauras(1)*tem1 + psauras(2)*tem2) * tem
473 c0 = (prauras(1)*tem1 + prauras(2)*tem2) * tem
474 if (ccwfac == zero) ccwfac = half
483 tem = one / prsi(ipt,ksfc)
493 if (flipv) ll = kp1 -l
494 sgc = prsl(ipt,ll) * tem
496 IF (sgc <= 0.050_kp) krmin = l
499 IF (sgc <= 0.760_kp) krmax = l
501 IF (sgc <= 0.970_kp) kfmax = l
503 IF (sgc <= 0.600_kp) kblmx = l
505 IF (sgc <= 0.980_kp) kblmn = l
510 if (fix_ncld_hr)
then
512 ncrnd = min(nrcmax, (krmax-krmin+1)) * (dtf/1800) + 0.10001_kp
520 ncrnd = min(nrcmax, (krmax-krmin+1))
522 ncrnd = min(nrcm,max(ncrnd, 1))
531 ic(nc) = ktem + 1 - nc
535 ic(nc) = ktem + 1 - nc
543 ii = mod(i-1,nrcm) + 1
544 irnd = (rannum(ipt,ii)-0.0005_kp)*(kcr-krmin+1)
545 ic(kfx+i) = irnd + krmin
576 prsm(l) = prsl(ipt,ll) * facmb
577 psjm(l) = prslk(ipt,ll)
578 phi_l(l) = phil(ipt,ll)
579 rhc_l(l) = rhc(ipt,ll)
582 uvi(l,ntr+1) = uin(ipt,ll)
583 uvi(l,ntr+2) = vin(ipt,ll)
588 uvi(l,n) = ccin(ipt,ll,n+2)
589 if (abs(uvi(l,n)) < 1.0e-20_kp) uvi(l,n) = zero
595 prs(ll) = prsi(ipt,l) * facmb
596 psj(ll) = prsik(ipt,l)
597 phi_h(ll) = phii(ipt,l)
600 if (ccin(ipt,1,2) <= -998.0_kp)
then
603 tem = ccin(ipt,ll,1) &
604 & * max(zero, min(one, (tcr-toi(l))*tcrf))
605 ccin(ipt,ll,2) = ccin(ipt,ll,1) - tem
612 qii(l) = ccin(ipt,ll,1)
613 qli(l) = ccin(ipt,ll,2)
616 kbl = max(min(k, kp1-kpbl(ipt)), k/2)
625 prsm(l) = prsl(ipt, l) * facmb
626 psjm(l) = prslk(ipt,l)
627 phi_l(l) = phil(ipt,l)
628 rhc_l(l) = rhc(ipt,l)
631 uvi(l,ntr+1) = uin(ipt,l)
632 uvi(l,ntr+2) = vin(ipt,l)
637 uvi(l,n) = ccin(ipt,l,n+2)
638 if (abs(uvi(l,n)) < 1.0e-20_kp) uvi(l,n) = zero
643 prs(l) = prsi(ipt,l) * facmb
644 psj(l) = prsik(ipt,l)
645 phi_h(l) = phii(ipt,l)
648 if (ccin(ipt,1,2) <= -998.0_kp)
then
650 tem = ccin(ipt,l,1) &
651 & * max(zero, min(one, (tcr-toi(l))*tcrf))
652 ccin(ipt,l,2) = ccin(ipt,l,1) - tem
658 qii(l) = ccin(ipt,l,1)
659 qli(l) = ccin(ipt,l,2)
680 dtvd(1,1) = cp*(toi(l)-toi(lm1)) + phi_l(l)-phi_l(lm1) &
681 & + alhl*(qoi(l)-qoi(lm1))
682 dtvd(1,2) = qoi(l) - qoi(lm1)
683 dtvd(1,3) = qli(l) - qli(lm1)
684 dtvd(1,4) = qii(l) - qii(lm1)
691 dtvd(2,1) = cp*(toi(l)-toi(lm1)) + phi_l(l)-phi_l(lm1) &
692 & + alhl*(qoi(l)-qoi(lm1))
696 if (abs(dtvd(2,1)) > 1.0e-10_kp)
then
697 tem1 = dtvd(1,1) / dtvd(2,1)
699 alfint(l,1) = one - half*(tem1 + tem2)/(one + tem2)
704 dtvd(1,1) = dtvd(2,1)
706 dtvd(2,2) = qoi(l) - qoi(lm1)
710 if (abs(dtvd(2,2)) > 1.0e-10_kp)
then
711 tem1 = dtvd(1,2) / dtvd(2,2)
713 alfint(l,2) = one - half*(tem1 + tem2)/(one + tem2)
715 dtvd(1,2) = dtvd(2,2)
717 dtvd(2,3) = qli(l) - qli(lm1)
721 if (abs(dtvd(2,3)) > 1.0e-10_kp)
then
722 tem1 = dtvd(1,3) / dtvd(2,3)
724 alfint(l,3) = one - half*(tem1 + tem2)/(one + tem2)
726 dtvd(1,3) = dtvd(2,3)
728 dtvd(2,4) = qii(l) - qii(lm1)
732 if (abs(dtvd(2,4)) > 1.0e-10_kp)
then
733 tem1 = dtvd(1,4) / dtvd(2,4)
735 alfint(l,4) = one - half*(tem1 + tem2)/(one + tem2)
737 dtvd(1,4) = dtvd(2,4)
743 dtvd(1,1) = uvi(l,n) - uvi(l-1,n)
745 dtvd(2,1) = uvi(l,n) - uvi(l-1,n)
749 if (abs(dtvd(2,1)) > 1.0e-10_kp)
then
750 tem1 = dtvd(1,1) / dtvd(2,1)
752 alfint(l,n+4) = one - half*(tem1 + tem2)/(one + tem2)
754 dtvd(1,1) = dtvd(2,1)
769 tem = one - max(pgfbot, min(pgftop, pgftop+pgfgrad*prsm(l)))
770 trcfac(l,ntr+1) = tem
771 trcfac(l,ntr+2) = tem
780 kbl = min(kbl, kblmn)
786 if (ib > kbl-1) cycle
864 CALL cloud(k, kp1, ib, ntrc, kblmx, kblmn &
865 &, frac, max_neg_bouy, vsmooth, aw_scal &
866 &, revap, wrkfun, calkbl, crtfun &
867 &, dt, kdt, tla, dpd &
868 &, alfint, rhfacl, rhfacs, area(ipt) &
869 &, ccwfac, cdrag(ipt), trcfac &
870 &, alfind, rhc_l, phi_l, phi_h, prs, prsm,sgcs &
871 &, toi, qoi, uvi, qli, qii, kbl, ddvel(ipt) &
872 &, tcu, qcu, rcu, pcu, flx, flxd, rain, wfnc, fscav_ &
873 &, trcmin, ntk-2, c0, wminras(1), c0i, wminras(2) &
881 ud_mf(ipt,ll) = ud_mf(ipt,ll) + flx(l+1)
882 dd_mf(ipt,ll) = dd_mf(ipt,ll) + flxd(l+1)
885 dt_mf(ipt,ll) = dt_mf(ipt,ll) + flx(ib)
887 if (mp_phys == mp_phys_mg)
then
889 cnv_mfd(ipt,ll) = cnv_mfd(ipt,ll) + flx(ib)/dt
893 cnv_dqldt(ipt,ll) = cnv_dqldt(ipt,ll) + flx(ib)* &
894 & max(0.,(qli(ib)+qii(ib)-qiid-qlid))/dt
896 if(flx(ib)<0)
write(*,*)
"AAA666", flx(ib),qli(ib),qii(ib) &
903 ud_mf(ipt,l) = ud_mf(ipt,l) + flx(l+1)
904 dd_mf(ipt,l) = dd_mf(ipt,l) + flxd(l+1)
906 dt_mf(ipt,ib) = dt_mf(ipt,ib) + flx(ib)
908 if (mp_phys == mp_phys_mg)
then
909 cnv_mfd(ipt,ib) = cnv_mfd(ipt,ib) + flx(ib)/dt
912 cnv_dqldt(ipt,ib) = cnv_dqldt(ipt,ib) + flx(ib)* &
913 & max(zero,(qli(ib)+qii(ib)-qiid-qlid))/dt
915 if(flx(ib)<0)
write(*,*)
"AAA666", flx(ib),qli(ib),qii(ib) &
927 if (.not. advcld)
then
929 clw(l) = clw(l) + qli(l)
930 cli(l) = cli(l) + qii(l)
938 rainc(ipt) = rain * 0.001_kp
939 if (rainc(ipt) < rain_min) rainc(ipt) = zero
953 if (sgcs(l) < 0.93_kp .and. abs(tcu(l)) > one_m10)
then
959 if (clw(l)+cli(l) > zero .OR. &
960 & qli(l)+qii(l) > clwmin) ktop(ipt) = l
963 if (clw(l)+cli(l) > zero .OR. &
964 & qli(l)+qii(l) > clwmin) kbot(ipt) = l
972 uin(ipt,ll) = uvi(l,ntr+1)
973 vin(ipt,ll) = uvi(l,ntr+2)
976 if (mp_phys == mp_phys_mg)
then
978 qlcn(ipt,ll) = max(qli(l)-ccin(ipt,ll,2), zero)
979 qicn(ipt,ll) = max(qii(l)-ccin(ipt,ll,1), zero)
980 cnv_fice(ipt,ll) = qicn(ipt,ll) &
981 & / max(1.0e-10_kp,qlcn(ipt,ll)+qicn(ipt,ll))
983 qlcn(ipt,ll) = qli(l)
984 qicn(ipt,ll) = qii(l)
985 cnv_fice(ipt,ll) = qii(l)/max(1.0e-10_kp,qii(l)+qli(l))
987 cf_upi(ipt,ll) = max(zero,min(0.02_kp*log(one+ &
988 & 500.0_kp*ud_mf(ipt,ll)/dt), cfmax))
990 clcn(ipt,ll) = cf_upi(ipt,ll)
991 w_upi(ipt,ll) = ud_mf(ipt,ll)*toi(l)*rgas / &
992 & (dt*max(cf_upi(ipt,ll),1.0e-12_kp)*prsl(ipt,ll))
997 ccin(ipt,ll,n+2) = uvi(l,n)
1004 ccin(ipt,ll,1) = qii(l)
1005 ccin(ipt,ll,2) = qli(l)
1010 ccin(ipt,ll,1) = ccin(ipt,ll,1) + cli(l)
1011 ccin(ipt,ll,2) = ccin(ipt,ll,2) + clw(l)
1015 ktop(ipt) = kp1 - ktop(ipt)
1016 kbot(ipt) = kp1 - kbot(ipt)
1023 uin(ipt,l) = uvi(l,ntr+1)
1024 vin(ipt,l) = uvi(l,ntr+2)
1027 if (mp_phys == mp_phys_mg)
then
1029 qlcn(ipt,l) = max(qli(l)-ccin(ipt,l,2), zero)
1030 qicn(ipt,l) = max(qii(l)-ccin(ipt,l,1), zero)
1031 cnv_fice(ipt,l) = qicn(ipt,l) &
1032 & / max(1.0e-10_kp,qlcn(ipt,l)+qicn(ipt,l))
1034 qlcn(ipt,l) = qli(l)
1035 qicn(ipt,l) = qii(l)
1036 cnv_fice(ipt,l) = qii(l)/max(1.0e-10_kp,qii(l)+qli(l))
1041 cf_upi(ipt,l) = max(zero,min(0.02_kp*log(one+ &
1042 & 500.0_kp*ud_mf(ipt,l)/dt), cfmax))
1044 clcn(ipt,l) = cf_upi(ipt,l)
1045 w_upi(ipt,l) = ud_mf(ipt,l)*toi(l)*rgas / &
1046 & (dt*max(cf_upi(ipt,l),1.0e-12_kp)*prsl(ipt,l))
1051 ccin(ipt,l,n+2) = uvi(l,n)
1057 ccin(ipt,l,1) = qii(l)
1058 ccin(ipt,l,2) = qli(l)
1062 ccin(ipt,l,1) = ccin(ipt,l,1) + cli(l)
1063 ccin(ipt,l,2) = ccin(ipt,l,2) + clw(l)
1070 ddvel(ipt) = ddvel(ipt) * ddfac * grav / (prs(kp1)-prs(k))
1074 deallocate (alfint, uvi, trcfac, rcu)
1081 & K, KP1, KD, NTRC, KBLMX, kblmn &
1082 &, FRACBL, MAX_NEG_BOUY, vsmooth, aw_scal &
1083 &, REVAP, WRKFUN, CALKBL, CRTFUN &
1084 &, DT, KDT, TLA, DPD &
1085 &, ALFINT, RHFACL, RHFACS, area, ccwf, cd, trcfac &
1086 &, alfind, rhc_ls, phil, phih, prs, prsm, sgcs &
1087 &, TOI, QOI, ROI, QLI, QII, KPBL, DSFC &
1088 &, TCU, QCU, RCU, PCU, FLX, FLXD, CUP, WFNC,fscav_ &
1089 &, trcmin, ntk, c0, qw0, c0i, qi0, dlq_fac)
1151 real (kind=kind_phys),
parameter :: rhmax=1.0_kp &
1152 &, quad_lam=1.0_kp &
1155 &, hcritd=4000.0_kp &
1157 &, hcrits=2500.0_kp &
1158 &, pcrit_lcl=250.0_kp &
1162 &, qudfac=quad_lam*half &
1168 &, dpnegcr = 150.0_kp
1172 real(kind=kind_phys),
parameter :: errmin=0.0001_kp &
1173 &, errmi2=0.1_kp*errmin &
1175 &, rainmin=1.0e-8_kp &
1176 &, oneopt9=one/0.09_kp &
1177 &, oneopt4=one/0.04_kp
1178 real(kind=kind_phys),
parameter :: almax=1.0e-2_kp &
1179 &, almin1=0.0_kp, almin2=0.0_kp
1180 real(kind=kind_phys),
parameter :: bldmax=300.0_kp, bldmin=25.0_kp
1185 LOGICAL REVAP, WRKFUN, CALKBL, CRTFUN, CALCUP
1186 logical vsmooth, aw_scal
1187 INTEGER K, KP1, KD, NTRC, kblmx, kblmn, ntk
1190 real(kind=kind_phys),
dimension(K) :: toi, qoi, prsm, qli, qii&
1191 &, phil, sgcs, rhc_ls &
1193 real(kind=kind_phys),
dimension(KP1) :: prs, phih
1194 real(kind=kind_phys),
dimension(K,NTRC) :: roi, trcfac
1195 real(kind=kind_phys),
dimension(ntrc) :: trcmin
1196 real(kind=kind_phys) :: cd, dsfc
1197 INTEGER :: KPBL, KBL, KB1, kdt
1199 real(kind=kind_phys) alfint(k,ntrc+4)
1200 real(kind=kind_phys) fracbl, max_neg_bouy, dpd &
1201 &, rhfacl, rhfacs, area, ccwf &
1202 &, c0, qw0, c0i, qi0, dlq_fac
1206 real(kind=kind_phys),
dimension(K) :: tcu, qcu, tcd, qcd, pcu
1207 real(kind=kind_phys),
dimension(KP1) :: flx, flxd
1208 real(kind=kind_phys),
dimension(K,NTRC) :: rcu
1209 real(kind=kind_phys) :: cup
1213 real(kind=kind_phys),
dimension(KD:K) :: hol, qol, hst, qst &
1214 &, tol, gmh, akt, akc, bkc, ltl, rnn &
1215 &, fco, pri, qil, qll, zet, xi, rns &
1216 &, q0u, q0d, vtf, cil, cll, etai, dlq &
1217 &, wrk1, wrk2, dhdp, qrb, qrt, evp &
1218 &, ghd, gsd, etz, cldfr, sigf, rho
1220 real(kind=kind_phys),
dimension(KD:KP1) :: gaf, gms, gam, dlb &
1221 &, dlt, eta, prl, buy, etd, hod, qod, wvl
1222 real(kind=kind_phys),
dimension(KD:K-1) :: etzi
1224 real(kind=kind_phys) fscav_(ntrc)
1226 LOGICAL ep_wfn, cnvflg, LOWEST, DDFT, UPDRET
1228 real(kind=kind_phys) alm, det, hcc, clp &
1229 &, hsu, hsd, qtl, qtv &
1230 &, akm, wfn, hos, qos &
1231 &, amb, tx1, tx2, tx3 &
1232 &, tx4, tx5, qis, qls &
1233 &, hbl, qbl, rbl(ntrc), wcbase &
1236 &, tx7, tx8, tx9, rhc &
1237 &, hstkd, qstkd, ltlkd, q0ukd, q0dkd, dlbkd &
1238 &, qtp, qw00, qi00, qrbkd &
1239 &, hstold, rel_fac, prism &
1240 &, tl, pl, ql, qs, dqs, st1, sgn, tau, &
1241 & qtvp, hb, qb, tb, qqq, &
1242 & hccp, ds, dh, ambmax, x00, epp, qtlp, &
1243 & dpi, dphib, dphit, del_eta, detp, &
1244 & tem, tem1, tem2, tem3, tem4, &
1245 & st2, st3, st4, st5, &
1246 & errh, errw, erre, tem5, &
1247 & tem6, hbd, qbd, st1s, shal_fac, hmax, hmin, &
1248 & dhdpmn, avt, avq, avr, avh &
1249 &, train, dof, cldfrd, tla, gmf &
1250 &, fac, rsum1, rsum2, rsum3, dpneg, hcrit &
1251 &, actevap,arearat,deltaq,mass,massinv,potevap &
1252 &, teq,qsteq,dqdt,qeq &
1253 &, clfrac, dt, clvfr, delzkm, fnoscav, delp
1257 INTEGER I, L, N, KD1, II, iwk, idh, lcon &
1258 &, IT, KM1, KTEM, KK, KK1, LM1, LL, LP1, kbls, kmxh &
1259 &, kblh, kblm, kblpmn, kmax, kmaxm1, kmaxp1, klcl, kmin, kmxb
1300 tol(l) = pt25*wrk1(l-1) + half*wrk1(l) + pt25*wrk1(l+1)
1301 qol(l) = pt25*wrk2(l-1) + half*wrk2(l) + pt25*wrk2(l+1)
1306 dpi = one / (prl(l+1) - prl(l))
1307 pri(l) = gravfac * dpi
1312 rho(l) = cmb2pa * pl / (rgas*tl*(one+nu*qol(l)))
1314 akt(l) = (prl(l+1) - pl) * dpi
1316 CALL qsatcn(tl, pl, qs, dqs)
1319 gam(l) = dqs * elocp
1321 gaf(l) = oneoalhl * gam(l) / st1
1323 ql = max(min(qs*rhmax,qol(l)), one_m10)
1327 ltl(l) = tem * st1 / (one+nu*(qst(l)+tl*dqs))
1328 vtf(l) = one + nu * ql
1329 eta(l) = one / (ltl(l) * vtf(l))
1331 hol(l) = tem + ql * alhl
1332 hst(l) = tem + qs * alhl
1348 dphib = phil(l) - phih(l+1)
1349 dphit = phih(l) - phil(l)
1351 dlb(l) = dphib * eta(l)
1352 dlt(l) = dphit * eta(l)
1357 eta(l) = eta(l+1) + dphib
1359 hol(l) = hol(l) + eta(l)
1361 hst(l) = hst(l) + eta(l)
1363 eta(l) = eta(l) + dphit
1370 dphib = phil(l) - phih(l+1)
1372 dlb(l) = dphib * eta(l)
1377 eta(l) = eta(l+1) + dphib
1379 hol(l) = hol(l) + eta(l)
1380 hst(l) = hst(l) + eta(l)
1385 if (sgcs(kd) < 0.5_kp)
then
1387 elseif (sgcs(kd) > 0.65_kp)
then
1390 hcrit = (hcrits*(sgcs(kd)-0.5_kp) + hcritd*(0.65_kp-sgcs(kd)))&
1394 ktem = max(kd+1, kblmx)
1398 if (hmin > hol(l))
then
1403 if (kmin == k)
return
1407 if (hmax < hol(l))
then
1413 if (kmax < kmin)
then
1417 elseif (kmax < k)
then
1420 if (abs(hol(kmax)-hol(l)) > hcrit)
then
1433 dhdp(l) = (hol(l)-hol(l+1)) / (prl(l+2)-prl(l))
1434 if (dhdp(l) < dhdpmn)
then
1437 elseif (dhdp(l) > zero .and. l <= kmin)
then
1442 if (kblpmn < kmax)
then
1444 if (hmax-hol(l) < half*hcrit)
then
1452 if (kmax > kd1)
then
1454 if (hmax > hst(l))
then
1465 if (kmax < kmxb)
then
1466 kmax = max(kd1, min(kmxb,k))
1476 tem = min(50.0_kp,max(10.0_kp,(prl(kmaxp1)-prl(kd))*0.10_kp))
1477 if (prl(kmaxp1) - prl(ii) > tem .and. ii > kbl) kbl = ii
1479 if (kbl .ne. ii)
then
1480 if (prl(kmaxp1)-prl(kbl) > bldmax) kbl = max(kbl,ii)
1483 if (hol(ii)-hol(ii-1) > half*hcrit) kbl = ii
1486 if (prl(kbl) - prl(klcl) > pcrit_lcl)
return
1489 kbl = min(kblmn, max(kbl,kblmx))
1515 kbl = min(kmax,max(kbl,kd+2))
1518 if (prl(kmaxp1)-prl(kbl) > bldmax .or. kb1 <= kd )
then
1523 pris = one / (prl(kp1)-prl(kbl))
1524 prism = one / (prl(kmaxp1)-prl(kbl))
1533 if (prl(kbl)-prl(kd) < 350.0_kp .and. kmax == k) shal_fac = shalfac
1536 eta(l) = (prl(kmaxp1)-prl(l)) * prism
1538 zet(l) = (eta(l) - tx1) * onebg
1539 xi(l) = zet(l) * zet(l) * (qudfac*shal_fac)
1540 eta(l) = zet(l) - zet(l+1)
1541 gms(l) = xi(l) - xi(l+1)
1550 hbl = hol(kmax) * eta(kmax)
1551 qbl = qol(kmax) * eta(kmax)
1552 qlb = cll(kmax) * eta(kmax)
1553 qib = cil(kmax) * eta(kmax)
1554 tx1 = qst(kmax) * eta(kmax)
1557 tem = eta(l) - eta(l+1)
1558 hbl = hbl + hol(l) * tem
1559 qbl = qbl + qol(l) * tem
1560 qlb = qlb + cll(l) * tem
1561 qib = qib + cil(l) * tem
1562 tx1 = tx1 + qst(l) * tem
1574 IF (hol(l) < tx2)
THEN
1583 tem1 = hbl - hol(kd)
1584 tem = hbl - hst(kd1) - ltl(kd1) * nu *(qol(kd1)-qst(kd1))
1589 if (hbl >= hst(l))
then
1595 if (lcon == kd .or. kbl <= kd .or. prl(kbl)-prsm(lcon) > 150.0_kp) &
1598 tx1 = rhfacs - qbl / tx1
1600 cnvflg = (tem > zero .OR. (lowest .AND. tem1 >= zero)) &
1603 IF (.NOT. cnvflg)
RETURN
1605 rhc = max(zero, min(one, exp(-20.0_kp*tx1) ))
1610 rbl(n) = roi(kmax,n) * eta(kmax)
1614 rbl(n) = rbl(n) + roi(l,n)*(eta(l)-eta(l+1))
1620 if (rbl(ntk) > zero)
then
1621 wcbase = min(two, max(wcbase, sqrt(twoo3*rbl(ntk))))
1631 tx3 = qst(kbl) - gaf(kbl) * hst(kbl)
1633 qil(l) = max(zero, min(one, (tcr-tcl-tol(l))*tcrf))
1638 tem = qst(l) - gaf(l) * hst(l)
1639 tem1 = (tx3 + tem) * half
1640 st2 = (gaf(l)+gaf(lp1)) * half
1642 fco(lp1) = tem1 + st2 * hbl
1644 rnn(lp1) = zet(lp1) * tem1 + st2 * tx4
1645 gmh(lp1) = xi(lp1) * tem1 + st2 * tx5
1648 tx4 = tx4 + eta(l) * hol(l)
1649 tx5 = tx5 + gms(l) * hol(l)
1651 qil(l) = max(zero, min(one, (tcr-tcl-tol(l))*tcrf))
1652 qll(lp1) = (half*alhf) * st2 * (qil(l)+qil(lp1)) + one
1660 tem = qst(l) - gaf(l) * hst(l)
1661 tem1 = (tx3 + tem) * half
1662 st2 = (gaf(l)+gaf(lp1)) * half
1664 fco(lp1) = tem1 + st2 * hbl
1665 rnn(lp1) = zet(lp1) * tem1 + st2 * tx4
1666 gmh(lp1) = xi(lp1) * tem1 + st2 * tx5
1668 fco(l) = tem + gaf(l) * hbl
1669 rnn(l) = tem * zet(l) + (tx4 + eta(l)*hol(l)) * gaf(l)
1670 gmh(l) = tem * xi(l) + (tx5 + gms(l)*hol(l)) * gaf(l)
1678 qil(kd) = max(zero, min(one, (tcr-tcl-tol(kd))*tcrf))
1679 qll(kd1) = (half*alhf) * st2 * (qil(kd) + qil(kd1)) + one
1680 qll(kd ) = alhf * gaf(kd) * qil(kd) + one
1684 if (c0ifac > 1.0e-6_kp) st2 = st2 * exp(c0ifac*min(tol(kd)-t0c,zero))
1685 tem = c0 * (one-st1)
1686 tem2 = st2*qi0 + tem*qw0
1690 tx2 = akt(l) * eta(l)
1693 fco(l) = fco(lp1) - fco(l) + tx1
1694 rnn(l) = rnn(lp1) - rnn(l) &
1695 & + eta(l)*(qol(l)+cll(l)+cil(l)) + tx1*zet(l)
1696 gmh(l) = gmh(lp1) - gmh(l) &
1697 & + gms(l)*(qol(l)+cll(l)+cil(l)) + tx1*xi(l)
1699 tem1 = (one-akt(l)) * eta(l)
1701 akt(l) = qll(l) + (st2 + tem) * tx2
1703 akc(l) = one / akt(l)
1705 st1 = half * (qil(l)+qil(lp1))
1707 if (c0ifac > 1.0e-6_kp) st2 = st2 * exp(c0ifac*min(tol(lp1)-t0c,zero))
1708 tem = c0 * (one-st1)
1709 tem2 = st2*qi0 + tem*qw0
1711 bkc(l) = qll(lp1) - (st2 + tem) * tem1
1715 fco(l) = fco(l) + tx1
1716 rnn(l) = rnn(l) + tx1*zet(lp1)
1717 gmh(l) = gmh(l) + tx1*xi(lp1)
1727 tx3 = bkc(kb1) * (qib + qlb)
1731 tem = bkc(l-1) * akc(l)
1732 tx3 = (tx3 + fco(l)) * tem
1733 tx4 = (tx4 + rnn(l)) * tem
1734 tx5 = (tx5 + gmh(l)) * tem
1737 hsd = hst(kd1) + ltl(kd1) * nu *(qol(kd1)-qst(kd1))
1742 tx3 = (tx3 + fco(kd)) * akc(kd)
1743 tx4 = (tx4 + rnn(kd)) * akc(kd)
1744 tx5 = (tx5 + gmh(kd)) * akc(kd)
1745 alm = alhf*qil(kd) - ltl(kd) * vtf(kd)
1747 hsu = hst(kd) + ltl(kd) * nu * (qol(kd)-qst(kd))
1756 tx1 = tx1 + tau * eta(l)
1757 tx2 = tx2 + tau * gms(l)
1762 hsu = hsu - alm * tx3
1771 cnvflg = hbl > hsu .and. abs(tx1) > 1.0e-4_kp
1775 st1 = half*(hsu + hsd)
1786 if (tx2 == zero)
then
1788 if (alm > almax) alm = -100.0_kp
1791 epp = tx1 * tx1 - (x00+x00)*st2
1792 if (epp > zero)
then
1795 tem1 = (-tx1-tem)*x00
1796 tem2 = (-tx1+tem)*x00
1797 if (tem1 > almax) tem1 = -100.0_kp
1798 if (tem2 > almax) tem2 = -100.0_kp
1799 alm = max(tem1,tem2)
1808 ELSEIF (hbl <= hsu .AND. hbl > st1)
THEN
1814 IF (almin1 > zero)
THEN
1815 IF (alm >= almin1) cnvflg = .false.
1818 IF ( (alm > zero) .OR. &
1819 & (.NOT. lowest .AND. alm == zero) ) cnvflg = .false.
1825 IF (ii > 0 .or. (qw00 == zero .and. qi00 == zero))
RETURN
1832 IF(clp > zero .AND. clp < one)
THEN
1833 st1 = half*(one+clp)
1844 hst(kd) = hst(kd)*st1 + hst(kd1)*st2
1845 hos = hol(kd)*st1 + hol(kd1)*st2
1846 qst(kd) = qst(kd)*st1 + qst(kd1)*st2
1847 qos = qol(kd)*st1 + qol(kd1)*st2
1848 qls = cll(kd)*st1 + cll(kd1)*st2
1849 qis = cil(kd)*st1 + cil(kd1)*st2
1850 ltl(kd) = ltl(kd)*st1 + ltl(kd1)*st2
1852 dlb(kd) = dlb(kd)*clp
1853 qrb(kd) = qrb(kd)*clp
1854 eta(kd) = eta(kd)*clp
1855 gms(kd) = gms(kd)*clp
1856 q0u(kd) = q0u(kd)*clp
1857 q0d(kd) = q0d(kd)*clp
1866 tem = prl(kd1) - (prl(kd1)-prl(kd)) * clp * half
1867 tx1 = prl(kbl) - tem
1868 tx2 = min(900.0_kp, max(tx1,100.0_kp))
1869 tem1 = log(tx2*0.01_kp) * oneolog10
1871 if ( kdt == 1 )
then
1873 rel_fac = (dt * facdt) / (tem1*6.0_kp + tem2*adjts_s)
1875 rel_fac = (dt * facdt) / (tem1*adjts_d + tem2*adjts_s)
1879 rel_fac = max(zero, min(half,rel_fac))
1882 iwk = tem*0.02_kp - 0.999999999_kp
1883 iwk = max(1, min(iwk, 16))
1884 acr = tx1 * (ac(iwk) + tem * ad(iwk)) * ccwf
1895 eta(l) = eta(l+1) + alm * (eta(l) + alm * gms(l))
1896 etai(l) = one / eta(l)
1907 qtl = qst(kb1) - gaf(kb1)*hst(kb1)
1919 del_eta = eta(l) - eta(lp1)
1920 hccp = hcc + del_eta*hol(l)
1922 qtlp = qst(lm1) - gaf(lm1)*hst(lm1)
1923 qtvp = half * ((qtlp+qtl)*eta(l) &
1924 & + (gaf(l)+gaf(lm1))*hccp)
1925 st1 = eta(l)*q0u(l) + eta(lp1)*q0d(l)
1926 detp = (bkc(l)*det - (qtvp-qtv) &
1927 & + del_eta*(qol(l)+cll(l)+cil(l)) + st1) * akc(l)
1929 tem1 = akt(l) - qll(l)
1930 tem2 = qll(lp1) - bkc(l)
1931 rns(l) = tem1*detp + tem2*det - st1
1933 qtp = half * (qil(l)+qil(lm1))
1934 tem2 = min(qtp*(detp-eta(l)*qw00), &
1935 & (one-qtp)*(detp-eta(l)*qi00))
1939 IF (rns(l) < zero .or. st1 < zero) ep_wfn = .true.
1940 IF (detp <= zero) cnvflg = .true.
1942 st1 = hst(l) - ltl(l)*nu*(qst(l)-qol(l))
1945 tem2 = hccp + detp * qtp * alhf
1947 st2 = ltl(l) * vtf(l)
1948 tem5 = cll(l) + cil(l)
1949 tem3 = (tx1 - eta(lp1)*st1 - st2*(det-tem5*eta(lp1))) * dlb(l)
1950 tem4 = (tem2 - eta(l )*st1 - st2*(detp-tem5*eta(l))) * dlt(l)
1955 akm = akm - min(st1,zero)
1957 if (st1 < zero .and. wfn < zero)
then
1958 dpneg = dpneg + prl(lp1) - prl(l)
1961 buy(l) = half * (tem3/(eta(lp1)*qrb(l)) + tem4/(eta(l)*qrt(l)))
1971 del_eta = eta(kd) - eta(kd1)
1972 hccp = hcc + del_eta*hos
1974 qtlp = qst(kd) - gaf(kd)*hst(kd)
1975 qtvp = qtlp*eta(kd) + gaf(kd)*hccp
1976 st1 = eta(kd)*q0u(kd) + eta(kd1)*q0d(kd)
1977 detp = (bkc(kd)*det - (qtvp-qtv) &
1978 & + del_eta*(qos+qls+qis) + st1) * akc(kd)
1980 tem1 = akt(kd) - qll(kd)
1981 tem2 = qll(kd1) - bkc(kd)
1982 rns(kd) = tem1*detp + tem2*det - st1
1984 IF (rns(kd) < zero) ep_wfn = .true.
1985 IF (detp <= zero) cnvflg = .true.
1990 IF ((qw00 == zero .and. qi00 == zero))
RETURN
1993 if (clp > zero .and. clp < one)
then
2004 fco(l) = fco(l) - q0u(l) - q0d(l)
2005 rnn(l) = rnn(l) - q0u(l)*zet(l) - q0d(l)*zet(lp1)
2006 gmh(l) = gmh(l) - q0u(l)*xi(l) - q0d(l)*zet(lp1)
2007 eta(l) = zet(l) - zet(lp1)
2008 gms(l) = xi(l) - xi(lp1)
2025 st1 = hst(kd) - ltl(kd)*nu*(qst(kd)-qos)
2026 st2 = ltl(kd) * vtf(kd)
2027 tem5 = (qls + qis) * eta(kd1)
2028 st1 = half * (tx1-eta(kd1)*st1-st2*(det-tem5))*dlb(kd)
2031 akm = akm - min(st1,zero)
2034 buy(kd) = st1 / (eta(kd1)*qrb(kd))
2044 IF (wfn >= zero) wfnc = wfn
2046 ELSEIF (.NOT. crtfun)
THEN
2054 tem = max(0.05_kp, min(cd*200.0_kp, max_neg_bouy))
2055 IF (.not. cnvflg .and. wfn > acr .and. &
2056 & dpneg < dpnegcr .and. akm <= tem) calcup = .true.
2060 IF (.NOT. calcup)
RETURN
2074 tem = one / (one + dlq_fac)
2076 rnn(l) = rns(l) * tem
2077 dlq(l) = rns(l) * tem * dlq_fac
2087 IF (dpd > zero)
THEN
2089 IF (clp > zero)
THEN
2091 train = train + rnn(l)
2095 pl = (prl(kd1) + prl(kd))*half
2096 IF (train > 1.0e-4_kp .AND. pl <= dpd*prl(kp1)) ddft = .true.
2102 &, tla, alfind, wcbase &
2103 &, tol, qol, hol, prl, qst, hst, gam, gaf &
2105 &, qrb, qrt, buy, kbl, idh, eta, rnn, etai &
2106 &, alm, wfn, train, ddft &
2107 &, etd, hod, qod, evp, dof, cldfr, etz &
2108 &, gms, gsd, ghd, wvl)
2115 IF (.NOT. ddft)
THEN
2148 st1 = (hcc+alhf*tem-eta(kd)*hst(kd)) / (one+gam(kd))
2149 ds = eta(kd1) * (hos- hol(kd)) - alhl*(qos - qol(kd))
2150 dh = eta(kd1) * (hos- hol(kd))
2153 gms(kd) = (ds + st1 - tem1*det*alhl-tem*alhf) * pri(kd)
2154 gmh(kd) = pri(kd) * (hcc-eta(kd)*hos + dh)
2158 qll(kd) = (tem2*(det-tem) + eta(kd1)*(qls-cll(kd)) &
2159 & + (one-qil(kd))*dlq(kd) - eta(kd)*qls ) * pri(kd)
2161 qil(kd) = (tem2*tem + eta(kd1)*(qis-cil(kd)) &
2162 & + qil(kd)*dlq(kd) - eta(kd)*qis ) * pri(kd)
2169 st1 = one - alfint(l,1)
2170 st2 = one - alfint(l,2)
2171 st3 = one - alfint(l,3)
2172 st4 = one - alfint(l,4)
2173 st5 = one - alfind(l)
2174 hb = alfint(l,1)*hol(lm1) + st1*hol(l)
2175 qb = alfint(l,2)*qol(lm1) + st2*qol(l)
2177 tem = alfint(l,4)*cil(lm1) + st4*cil(l)
2178 tem2 = alfint(l,3)*cll(lm1) + st3*cll(l)
2180 tem1 = eta(l) * (tem - cil(l))
2181 tem3 = eta(l) * (tem2 - cll(l))
2183 hbd = alfind(l)*hol(lm1) + st5*hol(l)
2184 qbd = alfind(l)*qol(lm1) + st5*qol(l)
2186 tem5 = etd(l) * (hod(l) - hbd)
2187 tem6 = etd(l) * (qod(l) - qbd)
2189 dh = eta(l) * (hb - hol(l)) + tem5
2190 ds = dh - alhl * (eta(l) * (qb - qol(l)) + tem6)
2192 gmh(l) = dh * pri(l)
2193 gms(l) = ds * pri(l)
2195 ghd(l) = tem5 * pri(l)
2196 gsd(l) = (tem5 - alhl * tem6) * pri(l)
2198 qll(l) = (tem3 + (one-qil(l))*dlq(l)) * pri(l)
2199 qil(l) = (tem1 + qil(l)*dlq(l)) * pri(l)
2201 tem1 = eta(l) * (cil(lm1) - tem)
2202 tem3 = eta(l) * (cll(lm1) - tem2)
2204 dh = eta(l) * (hol(lm1) - hb) - tem5
2205 ds = dh - alhl * eta(l) * (qol(lm1) - qb) &
2206 & + alhl * (tem6 - evp(lm1))
2208 gmh(lm1) = gmh(lm1) + dh * pri(lm1)
2209 gms(lm1) = gms(lm1) + ds * pri(lm1)
2211 ghd(lm1) = ghd(lm1) - tem5 * pri(lm1)
2212 gsd(lm1) = gsd(lm1) - (tem5-alhl*(tem6-evp(lm1))) * pri(lm1)
2214 qil(lm1) = qil(lm1) + tem1 * pri(lm1)
2215 qll(lm1) = qll(lm1) + tem3 * pri(lm1)
2217 avh = avh + gmh(lm1)*(prs(l)-prs(lm1))
2223 tem5 = etd(kp1) * (hod(kp1) - hbd)
2224 tem6 = etd(kp1) * (qod(kp1) - qbd)
2226 ds = dh + alhl * tem6
2228 tem2 = (ds - alhl * evp(k)) * pri(k)
2229 gmh(k) = gmh(k) + tem1
2230 gms(k) = gms(k) + tem2
2231 ghd(k) = ghd(k) + tem1
2232 gsd(k) = gsd(k) + tem2
2235 avh = avh + gmh(k)*(prs(kp1)-prs(k))
2237 tem4 = - gravfac * pris
2242 gmh(l) = gmh(l) + tx1
2243 gms(l) = gms(l) + tx2
2244 ghd(l) = ghd(l) + tx1
2245 gsd(l) = gsd(l) + tx2
2247 avh = avh + tx1*(prs(l+1)-prs(l))
2261 hol(l) = hol(l) + tem1*testmb
2262 qol(l) = qol(l) + (tem1-tem2) * testmboalhl
2263 hst(l) = hst(l) + tem2*(one+gam(l))*testmb
2264 qst(l) = qst(l) + tem2*gam(l) * testmboalhl
2265 cll(l) = cll(l) + qll(l) * testmb
2266 cil(l) = cil(l) + qil(l) * testmb
2269 if (alm > zero)
then
2270 hos = hos + gmh(kd) * testmb
2271 qos = qos + (gmh(kd)-gms(kd)) * testmboalhl
2272 qls = qls + qll(kd) * testmb
2273 qis = qis + qil(kd) * testmb
2276 hos = hos + (st1s*gmh(kd)+st2*gmh(kd1)) * testmb
2277 qos = qos + (st1s * (gmh(kd)-gms(kd)) &
2278 & + st2 * (gmh(kd1)-gms(kd1))) * testmboalhl
2279 hst(kd) = hst(kd) + (st1s*gms(kd)*(one+gam(kd)) &
2280 & + st2*gms(kd1)*(one+gam(kd1))) * testmb
2281 qst(kd) = qst(kd) + (st1s*gms(kd)*gam(kd) &
2282 & + st2*gms(kd1)*gam(kd1)) * testmboalhl
2284 qls = qls + (st1s*qll(kd)+st2*qll(kd1)) * testmb
2285 qis = qis + (st1s*qil(kd)+st2*qil(kd1)) * testmb
2289 tem = prl(kmaxp1) - prl(kmax)
2290 hbl = hol(kmax) * tem
2291 qbl = qol(kmax) * tem
2292 qlb = cll(kmax) * tem
2293 qib = cil(kmax) * tem
2295 tem = prl(l+1) - prl(l)
2296 hbl = hbl + hol(l) * tem
2297 qbl = qbl + qol(l) * tem
2298 qlb = qlb + cll(l) * tem
2299 qib = qib + cil(l) * tem
2317 qtl = qst(kb1) - gaf(kb1)*hst(kb1)
2321 tx4 = (alhf*half)*max(zero,min(one,(tcr-tcl-tol(kb1))*tcrf))
2330 del_eta = eta(l) - eta(lp1)
2331 hccp = hcc + del_eta*hol(l)
2333 qtlp = qst(lm1) - gaf(lm1)*hst(lm1)
2334 qtvp = half * ((qtlp+qtl)*eta(l) + (gaf(l)+gaf(lm1))*hccp)
2336 detp = (bkc(l)*tx1 - (qtvp-qtv) &
2337 & + del_eta*(qol(l)+cll(l)+cil(l)) &
2338 & + eta(l)*q0u(l) + eta(lp1)*q0d(l)) * akc(l)
2339 IF (detp <= zero) cnvflg = .true.
2341 st1 = hst(l) - ltl(l)*nu*(qst(l)-qol(l))
2343 tem2 = (alhf*half)*max(zero,min(one,(tcr-tcl-tol(lm1))*tcrf))
2344 tem1 = hccp + detp * (tem2+tx4)
2346 st2 = ltl(l) * vtf(l)
2347 tem5 = cll(l) + cil(l)
2349 & ( (tx2 -eta(lp1)*st1-st2*(tx1-tem5*eta(lp1))) * dlb(l) &
2350 & + (tem1 -eta(l )*st1-st2*(detp-tem5*eta(l))) * dlt(l) )
2367 st1 = hst(kd) - ltl(kd)*nu*(qst(kd)-qos)
2368 st2 = ltl(kd) * vtf(kd)
2369 tem5 = (qls + qis) * eta(kd1)
2370 akm = akm + half * (tx2-eta(kd1)*st1-st2*(tx1-tem5)) * dlb(kd)
2372 akm = (akm - wfn) * testmbi
2379 amb = - (wfn-acr) / akm
2383 amb = amb * clp * rel_fac
2389 ambmax = (prl(kmaxp1)-prl(kbl))*(fracbl*gravcon)
2390 amb = max(min(amb, ambmax),zero)
2398 if (amb > zero)
then
2406 tx1 = (0.2_kp / max(alm, 1.0e-5_kp))
2407 tx2 = one - min(one, pi * tx1 * tx1 / area)
2432 tx1 = max(1.0e-6_kp, abs(gms(kd) * onebcp * sigf(kd)))
2433 amb = min(tx1*amb, tfrac_max*toi(kd)) / tx1
2438 avr = dof * sigf(kbl)
2440 dsfc = dsfc + amb * etd(k) * (one/dt) * sigf(kbl)
2443 pcu(l) = pcu(l) + amb*rnn(l)*sigf(l)
2444 avr = avr + rnn(l) * sigf(l)
2446 pcu(k) = pcu(k) + amb * dof * sigf(kbl)
2451 tx2 = amb * oneoalhl
2453 delp = prs(l+1) - prs(l)
2455 st1 = gms(l) * tx1 * sigf(l)
2456 toi(l) = toi(l) + st1
2457 tcu(l) = tcu(l) + st1
2458 tcd(l) = tcd(l) + gsd(l) * tx1 * sigf(l)
2460 st1 = st1 - elocp * (qil(l) + qll(l)) * tx3
2462 avt = avt + st1 * delp
2464 flx(l) = flx(l) + eta(l) * tx3
2465 flxd(l) = flxd(l) + etd(l) * tx3
2467 qii(l) = qii(l) + qil(l) * tx3
2470 qli(l) = qli(l) + qll(l) * tx3 + tem
2472 st1 = (gmh(l)-gms(l)) * tx2 * sigf(l)
2474 qoi(l) = qoi(l) + st1
2475 qcu(l) = qcu(l) + st1
2476 qcd(l) = qcd(l) + (ghd(l)-gsd(l)) * tx2 * sigf(l)
2478 avq = avq + (st1 + (qll(l)+qil(l))*tx3) * delp
2481 avr = avr + (qll(l) + qil(l)) * delp * sigf(l) * gravcon
2484 if (qii(l) < zero)
then
2485 tem = qii(l) * elfocp
2486 qoi(l) = qoi(l) + qii(l)
2487 qcu(l) = qcu(l) + qii(l)
2488 toi(l) = toi(l) - tem
2489 tcu(l) = tcu(l) - tem
2492 if (qli(l) < zero)
then
2493 tem = qli(l) * elocp
2494 qoi(l) = qoi(l) + qli(l)
2495 qcu(l) = qcu(l) + qli(l)
2496 toi(l) = toi(l) - tem
2497 tcu(l) = tcu(l) - tem
2529 IF (l < idh .or. (.not. ddft))
THEN
2530 tem = tem + amb * rnn(l) * sigf(l)
2533 tem = tem + amb * dof * sigf(kbl)
2534 tem = tem * (3600.0_kp/dt)
2535 tem1 = sqrt(max(one, min(100.0_kp,(6.25e10_kp/max(area,one)))))
2537 clfrac = max(zero, min(half, rknob*
clf(tem)*tem1))
2542 IF (l >= idh .AND. ddft)
THEN
2544 tx2 = tx2 + tem * rnn(l)
2545 cldfrd = min(tem*cldfr(l), clfrac)
2547 tx1 = tx1 + amb * rnn(l) * sigf(l)
2549 tx4 = zfac * phil(l)
2550 tx4 = (one - tx4 * (one - half*tx4)) * afc
2552 IF (tx1 > zero .OR. tx2 > zero)
THEN
2555 pl = half * (prl(l+1)+prl(l))
2557 st1 = max(zero, min(one, (tcr-teq)*tcrf))
2558 st2 = st1*elfocp + (one-st1)*elocp
2560 CALL qsatcn ( teq,pl,qsteq,dqdt)
2562 deltaq = half * (qsteq*rhc_ls(l)-qeq) / (one+st2*dqdt)
2565 teq = teq - deltaq*st2
2567 tem1 = max(zero, min(one, (tcr-teq)*tcrf))
2568 tem2 = tem1*elfocp + (one-tem1)*elocp
2570 CALL qsatcn ( teq,pl,qsteq,dqdt)
2572 deltaq = (qsteq*rhc_ls(l)-qeq) / (one+tem2*dqdt)
2575 teq = teq - deltaq*tem2
2577 IF (qeq > qoi(l))
THEN
2578 potevap = (qeq-qoi(l))*(prl(l+1)-prl(l))*gravcon
2582 & tem4 = potevap * (one - exp( tx4*tx1**0.57777778_kp ))
2583 actevap = min(tx1, tem4*clfrac)
2585 if (tx1 < rainmin*dt) actevap = min(tx1, potevap)
2589 & tem4 = potevap * (one - exp( tx4*tx2**0.57777778_kp ))
2590 tem4 = min(min(tx2, tem4*cldfrd), potevap-actevap)
2591 if (tx2 < rainmin*dt) tem4 = min(tx2, potevap-actevap)
2595 st1 = (actevap+tem4) * pri(l)
2596 qoi(l) = qoi(l) + st1
2597 qcu(l) = qcu(l) + st1
2601 toi(l) = toi(l) - st1
2602 tcu(l) = tcu(l) - st1
2607 cup = cup + tx1 + tx2 + dof * amb * sigf(kbl)
2610 tx1 = tx1 + amb * rnn(l) * sigf(l)
2612 cup = cup + tx1 + dof * amb * sigf(kbl)
2619 if (etz(l) /= zero) etzi(l) = one / etz(l)
2632 st1 = one - alfind(l)
2633 hb = alfind(l) * hol(lm1) + st1 * hol(l)
2634 IF (etz(lm1) /= zero)
THEN
2636 IF (etd(l) > etd(lm1))
THEN
2637 hod(l) = (etd(lm1)*(hod(lm1)-hol(lm1)) &
2638 & + etd(l) *(hol(lm1)-hb) + etz(lm1)*hb) * tem
2640 hod(l) = (etd(lm1)*(hod(lm1)-hb) + etz(lm1)*hb) * tem
2648 hcc = hcc + (eta(l)-eta(l+1))*hol(l)
2656 if (fscav_(n) > zero)
then
2657 delzkm = ( phil(kd) - phih(kd1) ) *(onebg*0.001_kp)
2658 fnoscav = exp(- fscav_(n) * delzkm)
2663 gmh(kd) = pri(kd) * (hcc-eta(kd)*hol(kd)) * trcfac(kd,n) &
2666 if (fscav_(n) > zero)
then
2667 delzkm = ( phil(kd) - phih(l+1) ) *(onebg*0.001_kp)
2668 fnoscav = exp(- fscav_(n) * delzkm)
2671 st1 = one - alfint(l,n+4)
2672 st2 = one - alfind(l)
2673 hb = alfint(l,n+4) * hol(lm1) + st1 * hol(l)
2674 hbd = alfind(l) * hol(lm1) + st2 * hol(l)
2675 tem5 = etd(l) * (hod(l) - hbd)
2676 dh = eta(l) * (hb - hol(l)) * fnoscav + tem5
2677 gmh(l ) = dh * pri(l) * trcfac(l,n)
2678 dh = eta(l) * (hol(lm1) - hb) * fnoscav - tem5
2679 gmh(lm1) = gmh(lm1) + dh * pri(lm1) * trcfac(l,n)
2684 st1 = gmh(l)*amb*sigf(l) + st2
2686 st2 = st3 - trcmin(n)
2687 if (st2 < zero)
then
2688 roi(l,n) = trcmin(n)
2689 rcu(l,n) = rcu(l,n) + st1
2691 & st2 = st2 * (prl(l+1)-prl(l))*pri(l+1) * (cmb2pa/grav)
2694 rcu(l,n) = rcu(l,n) + st1
2709 &, TLA, ALFIND, wcbase &
2710 &, TOL, QOL, HOL, PRL, QST, HST, GAM, GAF &
2711! &, TOL, QOL, HOL, PRL, QST, HST, GAM, GAF, HBL, QBL&
2712 &, QRB, QRT, BUY, KBL, IDH, ETA, RNN, ETAI &
2713 &, ALM, WFN, TRAIN, DDFT &
2714 &, ETD, HOD, QOD, EVP, DOF, CLDFRD, WCB &
2715 &, GMS, GSD, GHD, wvlu)
2742 INTEGER K, KP1, KD, KBL
2743 real(kind=kind_phys) alfind(k), wcbase
2745 real(kind=kind_phys),
dimension(kd:k) :: hol, qol, hst, qst &
2746 &, tol, qrb, qrt, rnn &
2748 real(kind=kind_phys),
dimension(kd:kp1) :: gaf, buy, gam, eta &
2756 real(kind=kind_phys),
dimension(KD:K) :: rnf, wcb, evp, stlt &
2757 &, ghd, gsd, cldfrd &
2758 &, gqw, qrpi, qrps, bud
2760 real(kind=kind_phys),
dimension(KD:KP1) :: qrp, wvl, wvlu, etd &
2761 &, hod, qod, ror, gms
2763 real(kind=kind_phys) tl, pl, ql, qs, dqs, st1 &
2764 &, qqq, del_eta, hb, qb, tb &
2765 &, tem, tem1, tem2, tem3, tem4, st2 &
2766 &, errmin, errmi2, errh, errw, erre, tem5 &
2767 &, tem6, hbd, qbd, tx1, tx2, tx3 &
2768 &, tx4, tx5, tx6, tx7, tx8, tx9 &
2770 &, train, gmf, onpg, ctla, vtrm &
2771 &, rpart, qrmin, aa1, bb1, cc1, dd1 &
2773 &, wc2min, wcmin, f2, f3, f5 &
2774 &, gmf1, gmf5, qraf, qrbf, del_tla &
2775 &, tla, stla, ctl2, ctl3 &
2778 &, rnt, rnb, errq, rntp &
2779 &, edz, ddz, ce, qhs, fac, facg &
2780 &, rsum1, rsum2, rsum3, cee, dof, dofw
2783 INTEGER I, L, N, IX, KD1, II, kb1, IP1, JJ, ntla &
2784 &, IT, KM1, KTEM, KK, KK1, LM1, LL, LP1 &
2785 &, IDW, IDH, IDN(K), idnm, itr
2787 parameter(errmin=0.0001_kp, errmi2=0.1_kp*errmin)
2792 parameter(onpg=one+half, gmf=one/onpg, rpart=zero)
2797 parameter(aa1=1.0_kp, bb1=1.0_kp, cc1=1.0_kp, dd1=1.0_kp, &
2798 & f3=cc1, f5=1.0_kp)
2799 parameter(qrmin=1.0e-6_kp, wc2min=0.01_kp, gmf1=gmf/aa1, gmf5=gmf/f5)
2801 parameter(wcmin=sqrt(wc2min))
2804 integer,
parameter :: itrmu=25, itrmd=25 &
2805 &, itrmin=15, itrmnd=12, numtla=2
2808 real(kind=kind_phys),
parameter :: ddunc1=0.25_kp &
2809 &, ddunc2=one-ddunc1 &
2812 &, vtpexp=-0.3636_kp &
2813 &, vtp=36.34_kp*sqrt(1.2_kp)*(0.001_kp)**0.1364_kp
2816 real(kind=kind_phys) elm(k), aa(kd:k,kd:kp1), qw(kd:k,kd:k) &
2817 &, vt(2), vrw(2), trw(2), qa(3), wa(3)
2819 LOGICAL SKPUP, cnvflg, DDFT, UPDRET, DDLGK
2843 tx1 = (prl(kd) + prl(kd1)) * half
2844 ror(kd) = cmpor*tx1 / (tol(kd)*(one+nu*qol(kd)))
2846 gms(kd) = vtp *
vtpf(ror(kd))
2850 tem = tol(k) * (one + nu * qol(k))
2851 ror(kp1) = half * cmpor * (prl(kp1)+prl(k)) / tem
2852 gms(kp1) = vtp *
vtpf(ror(kp1))
2857 tem = half * (tol(l)+tol(l-1)) &
2858 & * (one + (half*nu) * (qol(l)+qol(l-1)))
2859 ror(l) = cmpor * prl(l) / tem
2861 gms(l) = vtp *
vtpf(ror(l))
2863 if (buy(l) <= zero .and. kk == kbl)
then
2869 buy(l) = 0.9_kp * buy(l-1)
2877 buy(l) = 0.25_kp * (qrpi(l-1)+qrpi(l)+qrpi(l)+qrpi(l+1))
2882 tx1 = 1000.0_kp + tx1 - prl(kp1)
2884 CALL angrad(tx1, alm, al2, tla)
2888 f2 = (bb1+bb1)*onebg/(pi*0.2_kp)
2894 del_tla = tla * 0.3_kp
2914 qw(kd,kd) = -qrb(kd) * gmf1
2915 ghd(kd) = eta(kd) * eta(kd)
2916 gqw(kd) = qw(kd,kd) * ghd(kd)
2917 gsd(kd) = etai(kd) * etai(kd)
2919 gqw(kk) = - qrb(kk-1) * (gmf1+gmf1)
2921 wcb(kk) = wcbase * wcbase
2929 ghd(l) = eta(l) * eta(l)
2930 gsd(l) = etai(l) * etai(l)
2931 gqw(l) = - ghd(l) * (qrb(l-1)+qrt(l)) * tem
2932 qw(l,l) = - qrt(l) * tem
2934 st1 = half * (eta(l) + eta(l+1))
2935 tx1 = tx1 + buy(l) * tem * (qrb(l)+qrt(l)) * st1 * st1
2936 wcb(l) = tx1 * gsd(l)
2939 tem1 = (qrb(kd) + qrt(kd1) + qrt(kd1)) * gmf1
2940 gqw(kd1) = - ghd(kd1) * tem1
2941 qw(kd1,kd1) = - qrt(kd1) * tem
2942 st1 = half * (eta(kd) + eta(kd1))
2943 wcb(kd) = (tx1 + buy(kd)*tem*qrb(kd)*st1*st1) * gsd(kd)
2947 qw(n,l) = gqw(l) * gsd(n)
2955 if (errq < 0.1_kp .or. tla > 45.0_kp) cycle
2958 stla = sin(tla*deg2rad)
2959 ctl2 = one - stla * stla
2961 stla = f2 * stla * al2
2963 ctl3 = 0.1364_kp * ctl2
2974 stlt(kbl) = one / wcbase
2985 IF (.NOT. skpup)
THEN
2991 qrpi(kbl) = one / qrp(kbl)
2993 tx1 = tx1 + qrp(l+1)*gqw(l+1)
2994 st1 = wcb(l) + qw(l,l)*qrp(l) + tx1*gsd(l)
2996 if (st1 > zero)
then
2997 wvl(l) = max(ddunc1*sqrt(st1) + ddunc2*wvl(l), wcmin)
3006 wvl(l) = max(wvl(l),wcmin)
3007 qrp(l) = (wvl(l)*wvl(l) - wcb(l) - tx1*gsd(l))/qw(l,l)
3011 qrp(l) = max(qrp(l), qrmin)
3013 stlt(l) = one / wvl(l)
3014 qrpi(l) = one / qrp(l)
3020 vt(1) = gms(kd) *
qrpf(qrp(kd))
3021 trw(1) = eta(kd) * qrp(kd) * stlt(kd)
3022 tx6 = trw(1) * vt(1)
3023 vrw(1) = f3*wvl(kd) - ctl2*vt(1)
3024 bud(kd) = stla * tx6 * qrb(kd) * half
3026 dof = 1.1364_kp * bud(kd) * qrpi(kd)
3027 dofw = -bud(kd) * stlt(kd)
3029 rnt = trw(1) * vrw(1)
3036 IF (rnt >= zero)
THEN
3037 tx3 = (rnt-ctl3*tx6) * qrpi(kd)
3038 tx5 = ctl2 * tx6 * stlt(kd)
3051 vt(2) = gms(l) *
qrpf(qrp(l))
3052 trw(2) = eta(l) * qrp(l) * stlt(l)
3053 vrw(2) = f3*wvl(l) - ctl2*vt(2)
3054 qqq = stla * trw(2) * vt(2)
3056 bud(l) = qqq * (st1 + qrt(l))
3060 dof = 1.1364_kp * bud(l) * qrpi(l)
3061 dofw = -bud(l) * stlt(l)
3063 rnf(ll) = rnf(ll) + qqq * st1
3064 rnf(l) = qqq * qrt(l)
3066 tem3 = vrw(1) + vrw(2)
3067 tem4 = trw(1) + trw(2)
3069 tx6 = pt25 * tem3 * tem4
3075 tem1 = pt25*(trw(1)*tem3 - tem4*vt(1))*qrpi(ll)
3076 st1 = pt25*(trw(1)*(ctl2*vt(1)-vrw(2)) &
3077 & * stlt(ll) + f3*trw(2))
3079 tem2 = pt25*(trw(2)*tem3 - tem4*vt(2))*qrpi(l)
3080 st2 = pt25*(trw(2)*(ctl2*vt(2)-vrw(1)) &
3081 & * stlt(l) + f3*trw(1))
3086 qa(2) = qa(2) + tx3 - tem1
3090 wa(2) = wa(2) + tx5 - st1
3102 IF (wvl(ktem) == wcmin) wa(1) = zero
3103 IF (wvl(ll) == wcmin) wa(2) = zero
3104 IF (wvl(l) == wcmin) wa(3) = zero
3106 aa(ll,n) = (wa(1)*qw(ktem,n) * stlt(ktem) &
3107 & + wa(2)*qw(ll,n) * stlt(ll) &
3108 & + wa(3)*qw(l,n) * stlt(l) ) * half
3110 aa(ll,ktem) = aa(ll,ktem) + qa(1)
3111 aa(ll,ll) = aa(ll,ll) + qa(2)
3112 aa(ll,l) = aa(ll,l) + qa(3)
3113 bud(ll) = (tx8 + rnn(ll)) * half &
3114 & - rnb + tx6 - bud(ll)
3115 aa(ll,kbl+1) = bud(ll)
3123 vt(2) = gms(l) *
qrpf(qrp(l))
3124 trw(2) = eta(l) * qrp(l) * stlt(l)
3125 vrw(2) = f3*wvl(l) - ctl2*vt(2)
3126 st1 = stla * trw(2) * vt(2) * qrb(ll)
3131 dof = 1.1364_kp * bud(l) * qrpi(l)
3132 dofw = -bud(l) * stlt(l)
3134 rnf(ll) = rnf(ll) + st1
3136 tem3 = vrw(1) + vrw(2)
3137 tem4 = trw(1) + trw(2)
3139 tx6 = pt25 * tem3 * tem4
3144 tem1 = pt25*(trw(1)*tem3 - tem4*vt(1))*qrpi(ll)
3145 st1 = pt25*(trw(1)*(ctl2*vt(1)-vrw(2)) &
3146 & * stlt(ll) + f3*trw(2))
3148 tem2 = pt25*(trw(2)*tem3 - tem4*vt(2))*qrpi(l)
3149 st2 = pt25*(trw(2)*(ctl2*vt(2)-vrw(1)) &
3150 & * stlt(l) + f3*trw(1))
3155 qa(2) = qa(2) + tx3 - tem1
3159 wa(2) = wa(2) + tx5 - st1
3169 IF (wvl(idw) == wcmin) wa(1) = zero
3170 IF (wvl(ll) == wcmin) wa(2) = zero
3171 IF (wvl(l) == wcmin) wa(3) = zero
3175 aa(ll,n) = (wa(1)*qw(kk,n) * stlt(kk) &
3176 & + wa(2)*qw(ll,n) * stlt(ll) &
3177 & + wa(3)*qw(l,n) * stlt(l) ) * half
3181 aa(ll,idw) = aa(ll,idw) + qa(1)
3182 aa(ll,ll) = aa(ll,ll) + qa(2)
3183 aa(ll,l) = aa(ll,l) + qa(3)
3184 bud(ll) = (tx8+rnn(ll)) * half - rnb + tx6 - bud(ll)
3186 aa(ll,l+1) = bud(ll)
3188 rnb = trw(2) * vrw(2)
3192 IF (rnb < zero)
THEN
3194 tem = vt(2) * trw(2)
3195 qa(2) = (rnb - ctl3*tem) * qrpi(kk)
3196 wa(2) = ctl2 * tem * stlt(kk)
3204 qa(2) = dof + tx3 - qa(2)
3208 wa(2) = dofw + tx5 - wa(2)
3212 IF (wvl(kk-1) == wcmin) wa(1) = zero
3213 IF (wvl(kk) == wcmin) wa(2) = zero
3217 aa(kk,n) = (wa(1)*qw(kk-1,n) * stlt(kk-1) &
3218 & + wa(2)*qw(kk,n) * stlt(kk)) * half
3224 aa(ll,lm1) = aa(ll,lm1) + qa(1)
3225 aa(ll,ll) = aa(ll,ll) + qa(2)
3226 bud(ll) = half*rnn(lm1) - tx6 + rnb - bud(ll)
3227 aa(ll,ll+1) = bud(ll)
3233 cnvflg = abs(aa(lm1,lm1)) < abs(aa(l,lm1))
3241 tx1 = aa(l,lm1) / aa(lm1,lm1)
3243 aa(l,n) = aa(l,n) - tx1 * aa(lm1,n)
3251 aa(kk,kk1) = aa(kk,kk1) / aa(kk,kk)
3252 tx2 = abs(aa(kk,kk1)) * qrpi(kk)
3259 tx1 = tx1 + aa(l,n) * aa(n,kk)
3261 aa(l,kk) = (aa(l,kk) - tx1) / aa(l,l)
3262 tx2 = max(tx2, abs(aa(l,kk))*qrpi(l))
3266 if (tx2 > one .and. abs(errq-tx2) > 0.1_kp)
then
3276 qrp(l) = max(qrp(l)+aa(l,kbl+1)*tem, qrmin)
3279 IF (itr < itrmin)
THEN
3281 IF (tem >= errmi2 .AND. tx2 >= errmin)
THEN
3290 IF (tem < zero .AND. errq > 0.5_kp)
THEN
3296 ELSEIF (tx2 < errmin)
THEN
3299 elseif (tem < zero .and. errq < 0.1_kp)
then
3318 IF (errq < 0.1_kp)
THEN
3338 tx1 = train / (tx1+rnt+rnb)
3339 IF (abs(tx1-one) < 0.2_kp)
THEN
3340 rnt = max(rnt*tx1,zero)
3343 rnf(l) = rnf(l) * tx1
3354 IF (.NOT. ddft)
then
3363 wvlu(kd:kp1) = wvl(kd:kp1)
3384 stlt(l) = eta(l) * stlt(l) * tem / ror(l)
3416 IF (rnt > zero)
THEN
3417 if (tx1 > zero)
THEN
3418 qrp(kd) = (rpart*rnt / (ror(kd)*tx1*gms(kd))) &
3419 & ** (one/1.1364_kp)
3421 tx1 = rpart*rnt / (ror(kd)*gms(kd)*qrp(kd)**1.1364_kp)
3423 rntp = (one - rpart) * rnt
3424 buy(kd) = - ror(kd) * tx1 * qrp(kd)
3436 ddlgk = idn(idnm) == idnmax
3437 if (.not. ddlgk) cycle
3439 st1 = one - alfind(l)
3440 wa(1) = alfind(l)*hol(l-1) + st1*hol(l)
3441 wa(2) = alfind(l)*qol(l-1) + st1*qol(l)
3442 wa(3) = alfind(l)*tol(l-1) + st1*tol(l)
3443 qa(2) = alfind(l)*hst(l-1) + st1*hst(l)
3444 qa(3) = alfind(l)*qst(l-1) + st1*qst(l)
3454 IF (l == kd1) fac = one
3456 facg = fac * half * gmf5
3463 wvl(l) = max(wvl(l-1),one_m1)
3465 qrp(l) = max(qrp(l-1),qrp(l))
3468 vt(1) = gms(l-1) *
qrpf(qrp(l-1))
3469 rnt = ror(l-1) * (wvl(l-1)+vt(1))*qrp(l-1)
3472 tem = max(alm,one_m6) * max(eta(l), one)
3474 trw(1) = picon*tem*(qrb(l-1)+qrt(l-1))
3475 trw(2) = one / trw(1)
3477 vrw(1) = half * (gam(l-1) + gam(l))
3478 vrw(2) = one / (vrw(1) + vrw(1))
3480 tx4 = (qrt(l-1)+qrb(l-1))*(onebg*fac*500.0_kp*eknob)
3482 dofw = one / (wa(3) * (one + nu*wa(2)))
3492 tx3 = stlt(l-1) * qrt(l-1) * (half*fac)
3493 tx8 = stlt(l) * qrb(l-1) * (half*fac)
3501 tem = wvl(l-1) + vt(1)
3502 IF (tem > zero)
THEN
3503 tem1 = one / (tem*ror(l-1))
3504 tx3 = vt(1) * tem1 * ror(l-1) * tx3
3511 IF (rnt > zero)
THEN
3512 tem = max(qrp(l-1),qrp(l))
3513 wvl(l) = tx1 * tem * qrb(l-1)*(facg*5.0_kp)
3515 wvl(l) = max(one_m2, wvl(l))
3516 trw(1) = trw(1) * half
3517 trw(2) = trw(2) + trw(2)
3519 IF (ddlgk) evp(l-1) = evp(l-2)
3536 rntp = rntp + rnt * tx1
3549 cnvflg = errq > errmin
3553 vt(1) = gms(l) *
qrpf(qrp(l))
3554 tem = wvl(l) + vt(1)
3556 IF (tem > zero)
THEN
3557 st1 = ror(l) * tem * qrp(l) + rnt
3558 IF (st1 /= zero) st1 = two * evp(l-1) / st1
3559 tem1 = one / (tem*ror(l))
3560 tem2 = vt(1) * tem1 * ror(l) * tx8
3568 tem = ror(l)*wvl(l) - ror(l-1)*wvl(l-1)
3569 if (tem > zero)
then
3570 tx5 = (tx1 - st1 + tem2 + tx3)/(one+tem*tem1)
3572 tx5 = tx1 - tem*tx6 - st1 + tem2 + tx3
3575 tx5 = half * (tx5 + st2)
3586 etd(l) = ror(l) * tx5 * max(wvl(l),zero)
3588 if (etd(l) > zero) etd(l) = half * (etd(l) + tem1)
3591 del_eta = etd(l) - etd(l-1)
3601 erre = etd(l) - tem1
3603 tem = max(abs(del_eta), trw(1))
3604 tem2 = del_eta / tem
3605 tem1 = sqrt(max((tem+del_eta)*(tem-del_eta),zero))
3608 edz = (half + asin(tem2)*piinv)*del_eta + tem1*piinv
3611 wcb(l-1) = etd(l) + ddz
3614 IF (del_eta > zero)
THEN
3615 qqq = one / (etd(l) + ddz)
3616 hod(l) = (etd(l-1)*hod(l-1) + del_eta*hol(l-1) &
3617 & + ddz*wa(1)) * qqq
3618 qod(l) = (etd(l-1)*qod(l-1) + del_eta*qol(l-1) &
3619 & + ddz*wa(2)) * qqq
3620 ELSEif((etd(l-1) + edz) > zero)
then
3621 qqq = one / (etd(l-1) + edz)
3622 hod(l) = (etd(l-1)*hod(l-1) + edz*wa(1)) * qqq
3623 qod(l) = (etd(l-1)*qod(l-1) + edz*wa(2)) * qqq
3625 errh = hod(l) - tem1
3626 errq = abs(errh/hod(l)) + abs(erre/max(etd(l),one_m5))
3635 qhs = qa(3) + half * (gaf(l-1)+gaf(l)) * (hod(l)-qa(2))
3639 st2 = prl(l) * (qhs + tem1 * (qhs-qod(l)))
3640 tem2 = ror(l) * qrp(l)
3641 CALL qrabf(tem2,qraf,qrbf)
3642 tem6 = tx5 * (1.6_kp + 124.9_kp * qraf) * qrbf * tx4
3644 ce = tem6 * st2 / ((5.4e5_kp*st2 + 2.55e6_kp)*(etd(l)+ddz))
3646 tem2 = - ((one+tem1)*(qhs+ce) + tem1*qod(l))
3647 tem3 = (one + tem1) * qhs * (qod(l)+ce)
3648 tem = max(tem2*tem2 - four*tem1*tem3,zero)
3649 qod(l) = max(tem4, (- tem2 - sqrt(tem)) * vrw(2))
3654 st2 = prl(l) * (qhs + tem1 * (qhs-qod(l)))
3655 ce = tem6 * st2 / ((5.4e5_kp*st2 + 2.55e6_kp)*(etd(l)+ddz))
3660 tem2 = - ((one+tem1)*(qhs+ce) + tem1*tem4)
3661 tem3 = (one + tem1) * qhs * (tem4+ce)
3662 tem = max(tem2*tem2 - four*tem1*tem3,zero)
3663 qod(l) = max(tem4, (- tem2 - sqrt(tem)) * vrw(2))
3666 evp(l-1) = (qod(l)-tem4) * (etd(l)+ddz)
3668 qa(1) = tx1*rnt + rnf(l-1) - evp(l-1)
3670 if (qa(1) > zero)
then
3671 IF (etd(l) > zero)
THEN
3672 tem = qa(1) / (etd(l)+ror(l)*tx5*vt(1))
3673 qrp(l) = max(tem,zero)
3674 ELSEIF (tx5 > zero)
THEN
3675 qrp(l) = (max(zero,qa(1)/(ror(l)*tx5*gms(l)))) &
3676 & ** (one/1.1364_kp)
3681 qrp(l) = half * qrp(l)
3684 tem1 = wa(3) + (hod(l)-wa(1)-alhl*(qod(l)-wa(2))) &
3687 tem1 = tem1 * (one + nu*qod(l))
3688 ror(l) = cmpor * prl(l) / tem1
3692 buy(l) = (tem1 - one - qrp(l)) * ror(l) * tx5
3696 wvl(l) = vt(2) * (etd(l-1)*wvl(l-1) - facg &
3697 & * (buy(l-1)*qrt(l-1)+buy(l)*qrb(l-1)))
3700 if (wvl(l) < zero)
then
3707 wvl(l) = half*(wvl(l)+tem1)
3713 errw = wvl(l) - tem1
3715 errq = errq + abs(errw/max(wvl(l),one_m5))
3719 IF (itr >= min(itrmnd,itrmd/2))
THEN
3721 IF (etd(l-1) == zero .AND. errq > 0.2_kp)
THEN
3733 tx5 = (stlt(kb1) * qrt(kb1) &
3734 & + stlt(kbl) * qrb(kb1)) * (0.5_kp*fac)
3738 tem = max(tx1*rnt+rnf(l-1),zero)
3739 qa(1) = tem - evp(l-1)
3742 qrp(l) = (qa(1) / (ror(l)*tx5*gms(l))) &
3743 & ** (one/1.1364_kp)
3745 buy(l) = - ror(l) * tx5 * qrp(l)
3749 del_eta = etd(l) - etd(l-1)
3750 IF(del_eta < zero .AND. errq > 0.1_kp)
THEN
3757 del_eta = - etd(l-1)
3768 qhs = qa(3) + half * (gaf(l-1)+gaf(l)) &
3774 st2 = prl(l) * (qhs + tem1 * (qhs-qod(l)))
3775 tem2 = ror(l) * qrp(l-1)
3776 CALL qrabf(tem2,qraf,qrbf)
3777 tem6 = tx5 * (1.6_kp + 124.9_kp * qraf) * qrbf * tx4
3779 ce = tem6*st2/((5.4e5_kp*st2 + 2.55e6_kp)*(etd(l)+ddz))
3782 tem2 = - ((one+tem1)*(qhs+ce) + tem1*qod(l))
3783 tem3 = (one + tem1) * qhs * (qod(l)+ce)
3784 tem = max(tem2*tem2 -four*tem1*tem3,zero)
3785 qod(l) = max(tem4, (- tem2 - sqrt(tem)) * vrw(2))
3789 st2 = prl(l) * (qhs + tem1 * (qhs-qod(l)))
3790 ce = tem6*st2/((5.4e5_kp*st2 + 2.55e6_kp)*(etd(l)+ddz))
3795 tem2 = - ((one+tem1)*(qhs+ce) + tem1*tem4)
3796 tem3 = (one + tem1) * qhs * (tem4+ce)
3797 tem = max(tem2*tem2 -four*tem1*tem3,zero)
3798 qod(l) = max(tem4, (- tem2 - sqrt(tem)) * vrw(2))
3802 evp(l-1) = (qod(l)-tem4) * (etd(l)+ddz)
3807 qa(1) = tx1*rnt + rnf(l-1)
3808 evp(l-1) = min(evp(l-1), qa(1))
3809 qa(1) = qa(1) - evp(l-1)
3843 IF (etd(l-1) == zero .AND. errq > 0.1_kp .and. l <= kbl)
THEN
3853 qa(1) = tx1*rnt + rnf(l-1)
3854 evp(l-1) = min(evp(l-1), qa(1))
3855 qa(1) = qa(1) - evp(l-1)
3865 qrp(l) = (qa(1) / (ror(l)*tx5*gms(l))) &
3866 & ** (one/1.1364_kp)
3870 st1 = one - alfind(l)
3873 buy(l) = - ror(l) * tx5 * qrp(l)
3878 ll = min(idn(idnm), kp1)
3879 IF (errq < one .AND. l <= ll)
THEN
3880 IF (etd(l-1) > zero .AND. etd(l) == zero)
THEN
3883 if (l < kbl .or. tx5 > zero) idnm = idnm + 1
3886 if (etd(l) == zero .and. l > kbl)
then
3888 if (tx5 > zero) idnm = idnm + 1
3896 IF (errq > 0.1_kp .AND. idn(idnm) == idnmax) ddft = .false.
3899 IF (.NOT. ddft)
RETURN
3921 dof = max(dof, zero)
3924 tx2 = rntp + rnb + dof
3927 IF (ii >= kd1+1)
THEN
3928 rnn(kd) = rnn(kd) + rnf(kd)
3936 IF (l > kd1 .AND. l < ii)
THEN
3938 tx2 = tx2 + rnn(l-1)
3939 ELSEIF (l >= ii .AND. l < idn(idnm))
THEN
3943 ELSEIF (l >= idn(idnm))
THEN
3948 rnn(l) = rnf(l) + rns(l)
3957 rnn(l) = rnn(l) + rnb
3965 IF (tx1 > zero)
THEN
3966 tx1 = (train - tx2) / tx1
3972 evp(l) = evp(l) * tx1