67 us,vs,zo,t,q,z1,tn,qo,po,psur,dhdt,kpbl,rho, & ! input variables, must be supplied
68 hfx,qfx,xland,ichoice,tcrit,dtime, &
69 zuo,xmb_out,kbcon,ktop,k22,ierr,ierrc, &
70 flag_init, flag_restart,fv,r_d,delp,tmf,qmicro, &
71 forceqv_spechum,betascu,betamcu,betadcu,sigmain,&
72 sigmaout,progsigma,dx, &
73 outt,outq,outqc,outu,outv,cnvwt,pre,cupclw, & ! output tendencies
74 itf,ktf,its,ite, kts,kte,ipr,tropics)
88 logical,
intent(in) :: flag_init, flag_restart, progsigma
89 logical :: make_calc_for_xk = .true.
90 integer,
intent (in ) :: &
99 real(kind=kind_phys),
dimension (its:,kts:) &
100 ,
intent (inout ) :: &
101 cnvwt,outt,outq,outqc,cupclw,zuo,outu,outv
103 real(kind=kind_phys),
dimension (its:,kts:) &
106 real(kind=kind_phys),
dimension (its:,kts:) &
107 ,
intent (in ),
optional :: &
108 qmicro, sigmain, forceqv_spechum
110 real(kind=kind_phys),
dimension (its:) &
113 integer,
dimension (its:) &
114 ,
intent (inout ) :: &
116 integer,
dimension (its:) &
119 integer,
dimension (its:) &
127 real(kind=kind_phys),
dimension (its:,kts:) &
129 t,po,tn,dhdt,rho,us,vs,delp
130 real(kind=kind_phys),
dimension (its:,kts:) &
133 real(kind=kind_phys),
dimension (its:) &
135 xland,z1,psur,hfx,qfx,dx
137 real(kind=kind_phys) &
139 dtime,tcrit,fv,r_d,betascu,betamcu,betadcu
141 real(kind=kind_phys),
dimension (its:,kts:) &
142 ,
intent (out),
optional :: &
187 real(kind=kind_phys),
dimension (its:ite,kts:kte) :: &
188 entr_rate_2d,he,hes,qes,z, &
190 xhe,xhes,xqes,xz,xt,xq, &
191 qes_cup,q_cup,he_cup,hes_cup,z_cup,p_cup,gamma_cup,t_cup, &
192 qeso_cup,qo_cup,heo_cup,heso_cup,zo_cup,po_cup,gammao_cup, &
194 xqes_cup,xq_cup,xhe_cup,xhes_cup,xz_cup, &
196 dbyo,qco,pwo,hco,qrco, &
204 cd,dellah,dellaq,dellat,dellaqc,uc,vc,dellu,dellv,u_cup,v_cup, &
205 wu2,omega_u,zeta,zdqca,del,clw_all
225 real(kind=kind_phys),
dimension (its:ite) :: &
226 zws,ztexec,zqexec,pre,aa1,aa0,xaa0,hkb, &
227 flux_tun,hkbo,xhkb, &
228 rand_vmas,xmbmax,xmb, &
230 cap_max_increment,lambau,wc,omegac,sigmab, &
232 integer,
dimension (its:ite) :: &
233 kstabi,xland1,kbmax,ktopx
242 logical :: flag_shallow,flag_mid
243 logical,
dimension(its:ite) :: cnvflg
246 real(kind=kind_phys) :: &
248 cap_maxs,trash,trash2,frh,el2orc,gravinv
250 real(kind=kind_phys) buo_flux,pgeoh,dp,entup,detup,totmas
251 real(kind=kind_phys) :: &
252 sigmind,sigminm,sigmins
253 parameter(sigmind=0.005,sigmins=0.03,sigminm=0.01)
255 real(kind=kind_phys) xff_shal(3),blqe,xkshal
256 character*50 :: ierrc(its:)
257 real(kind=kind_phys),
dimension (its:ite,kts:kte) :: &
258 up_massentr,up_massdetr,up_massentro,up_massdetro,up_massentru,up_massdetru
260 real(kind=kind_phys) :: c_up,x_add,qaver,dts,fp,fpi
261 real(kind=kind_phys),
dimension (its:ite,kts:kte) :: c1d,dtempdz
262 integer,
dimension (its:ite,kts:kte) :: k_inv_layers
263 integer,
dimension (its:ite) :: start_level, pmin_lev
266 real(kind=kind_phys),
parameter :: zero = 0
277 el2orc=xlv*xlv/(r_v*cp)
281 xland1(i)=int(xland(i)+.001)
283 if(xland(i).gt.1.5 .or. xland(i).lt.0.5)
then
291 cap_max_increment(i)=25.
318 cd(i,k)=.75*entr_rate(i)
348 buo_flux= (hfx(i)/cp+0.608*t(i,1)*qfx(i)/xlv)/rho(i,1)
351 zws(i) = max(0.,flux_tun(i)*0.41*buo_flux*zo(i,2)*g/t(i,1))
352 if(zws(i) > tiny(pgeoh))
then
354 zws(i) = 1.2*zws(i)**.3333
356 ztexec(i) = max(flux_tun(i)*hfx(i)/(rho(i,1)*zws(i)*cp),0.0)
358 zqexec(i) = max(flux_tun(i)*qfx(i)/xlv/(rho(i,1)*zws(i)),0.)
362 zws(i) = max(0.,flux_tun(i)*0.41*buo_flux*zo(i,kpbl(i))*g/t(i,kpbl(i)))
363 zws(i) = 1.2*zws(i)**.3333
364 zws(i) = zws(i)*rho(i,kpbl(i))
375 call cup_env(z,qes,he,hes,t,q,po,z1, &
376 psur,ierr,tcrit,-1, &
379 call cup_env(zo,qeso,heo,heso,tn,qo,po,z1, &
380 psur,ierr,tcrit,-1, &
387 call cup_env_clev(t,qes,q,he,hes,z,po,qes_cup,q_cup,he_cup, &
388 hes_cup,z_cup,p_cup,gamma_cup,t_cup,psur, &
392 call cup_env_clev(tn,qeso,qo,heo,heso,zo,po,qeso_cup,qo_cup, &
393 heo_cup,heso_cup,zo_cup,po_cup,gammao_cup,tn_cup,psur, &
401 u_cup(i,kts)=us(i,kts)
402 v_cup(i,kts)=vs(i,kts)
404 u_cup(i,k)=.5*(us(i,k-1)+us(i,k))
405 v_cup(i,k)=.5*(vs(i,k-1)+vs(i,k))
415 if(zo_cup(i,k).gt.zkbmax+z1(i))
then
422 kbmax(i)=min(kbmax(i),ktf/2)
434 if(kpbl(i).gt.3)cap_max(i)=po_cup(i,kpbl(i))
436 k22(i)=maxloc(heo_cup(i,2:kbmax(i)),1)
438 if(k22(i).gt.kbmax(i))
then
441 ierrc(i)=
"could not find k22"
457 x_add = xlv*zqexec(i)+cp*ztexec(i)
458 call get_cloud_bc(kte,he_cup(i,1:kte),hkb(i),k22(i),x_add)
459 call get_cloud_bc(kte,heo_cup(i,1:kte),hkbo(i),k22(i),x_add)
476 call cup_kbcon(ierrc,cap_max_increment,5,k22,kbcon,heo_cup,heso_cup, &
477 hkbo,ierr,kbmax,po_cup,cap_max, &
481 z_cup,entr_rate,heo,0)
484 call cup_minimi(heso_cup,kbcon,kbmax,kstabi,ierr, &
489 kbcon,kstabi,dtempdz,itf,ktf,its,ite, kts,kte)
494 entr_rate_2d(i,:)=entr_rate(i)
496 start_level(i)=k22(i)
497 x_add = xlv*zqexec(i)+cp*ztexec(i)
498 call get_cloud_bc(kte,he_cup(i,1:kte),hkb(i),k22(i),x_add)
499 if(kbcon(i).gt.ktf-4)
then
503 frh = 2.*min(qo_cup(i,k)/qeso_cup(i,k),1.)
504 entr_rate_2d(i,k)=entr_rate(i)
505 cd(i,k)=.75*entr_rate_2d(i,k)
512 if(kpbl(i).lt.5)kstart=kbcon(i)
515 if(k_inv_layers(i,1).gt.0 .and. &
516 (po_cup(i,kstart)-po_cup(i,k_inv_layers(i,1))).lt.200.)
then
517 ktop(i)=k_inv_layers(i,1)
520 if((po_cup(i,kstart)-po_cup(i,k)).gt.200.)
then
530 call rates_up_pdf(rand_vmas,ipr,
'shallow',ktop,ierr,po_cup,entr_rate_2d,hkbo,heo,heso_cup,zo_cup, &
531 xland1,kstabi,k22,kbcon,its,ite,itf,kts,kte,ktf,zuo,kpbl,ktopx,kbcon,pmin_lev)
551 do k=maxloc(zuo(i,:),1),ktop(i)
552 if(zuo(i,k).lt.1.e-6)
then
576 ,ierr,ktop,zo_cup,zuo,cd,entr_rate_2d &
577 ,up_massentro, up_massdetro ,up_massentr, up_massdetr &
578 ,2,kbcon,k22,up_massentru,up_massdetru,lambau)
591 if(ierr(i) /= 0) cycle
592 do k=1,start_level(i)
596 do k=1,start_level(i)-1
598 hco(i,k)=heo_cup(i,k)
611 if(ierr(i) /= 0) cycle
613 do k=start_level(i)+1,ktop(i)
614 hc(i,k)=(hc(i,k-1)*zu(i,k-1)-.5*up_massdetr(i,k-1)*hc(i,k-1)+ &
615 up_massentr(i,k-1)*he(i,k-1)) / &
616 (zu(i,k-1)-.5*up_massdetr(i,k-1)+up_massentr(i,k-1))
617 uc(i,k)=(uc(i,k-1)*zu(i,k-1)-.5*up_massdetr(i,k-1)*uc(i,k-1)+ &
618 up_massentr(i,k-1)*us(i,k-1)) / &
619 (zu(i,k-1)-.5*up_massdetr(i,k-1)+up_massentr(i,k-1))
620 vc(i,k)=(vc(i,k-1)*zu(i,k-1)-.5*up_massdetr(i,k-1)*vc(i,k-1)+ &
621 up_massentr(i,k-1)*vs(i,k-1))/ &
622 (zu(i,k-1)-.5*up_massdetr(i,k-1)+up_massentr(i,k-1))
623 dby(i,k)=max(0.,hc(i,k)-hes_cup(i,k))
624 hco(i,k)=(hco(i,k-1)*zuo(i,k-1)-.5*up_massdetro(i,k-1)*hco(i,k-1)+ &
625 up_massentro(i,k-1)*heo(i,k-1)) / &
626 (zuo(i,k-1)-.5*up_massdetro(i,k-1)+up_massentro(i,k-1))
627 dbyo(i,k)=hco(i,k)-heso_cup(i,k)
628 dz=zo_cup(i,k+1)-zo_cup(i,k)
629 if(k.ge.kbcon(i))dbyt(i,k)=dbyt(i,k-1)+dbyo(i,k)*dz
631 ki=maxloc(dbyt(i,:),1)
632 if(ktop(i).gt.ki+1)
then
634 zuo(i,ktop(i)+1:ktf)=0.
635 zu(i,ktop(i)+1:ktf)=0.
636 cd(i,ktop(i)+1:ktf)=0.
637 up_massdetro(i,ktop(i))=zuo(i,ktop(i))
639 up_massentro(i,ktop(i):ktf)=0.
640 up_massdetro(i,ktop(i)+1:ktf)=0.
641 entr_rate_2d(i,ktop(i)+1:ktf)=0.
646 if(ktop(i).lt.kbcon(i)+1)
then
649 ierrc(i)=
'ktop is less than kbcon+1'
653 if(ktop(i).gt.ktf-2)
then
656 ierrc(i)=
"ktop is larger than ktf-2"
661 call get_cloud_bc(kte,qo_cup(i,1:kte),qaver,k22(i),zero)
662 qaver = qaver + zqexec(i)
663 do k=1,start_level(i)-1
664 qco(i,k)= qo_cup(i,k)
670 do k=start_level(i)+1,ktop(i)
671 trash=qeso_cup(i,k)+(1./xlv)*(gammao_cup(i,k) &
672 /(1.+gammao_cup(i,k)))*dbyo(i,k)
675 qco(i,k)= (trash2* ( zuo(i,k-1)-0.5*up_massdetr(i,k-1)) + &
676 up_massentr(i,k-1)*qo(i,k-1)) / &
677 (zuo(i,k-1)-.5*up_massdetr(i,k-1)+up_massentr(i,k-1))
679 if(qco(i,k)>=trash )
then
680 dz=z_cup(i,k)-z_cup(i,k-1)
683 clw_all(i,k)=max(0._kind_phys,qco(i,k)-trash)
684 qrco(i,k)= (qco(i,k)-trash)/(1.+(c0_shal+c1d(i,k))*dz)
685 if(qrco(i,k).lt.0.)
then
689 pwo(i,k)=c0_shal*dz*qrco(i,k)*zuo(i,k)
691 qco(i,k)= trash+qrco(i,k)
696 cupclw(i,k)=qrco(i,k)
701 do k=k22(i)+1,ktop(i)
702 dp=100.*(po_cup(i,k)-po_cup(i,k+1))
703 cnvwt(i,k)=zuo(i,k)*cupclw(i,k)*g/dp
705 trash2=trash2+entr_rate_2d(i,k)
707 qco(i,k)=qco(i,k)-qrco(i,k)
710 do k=k22(i)+1,max(kbcon(i),k22(i)+1)
712 trash=trash+entr_rate_2d(i,k)
717 hco(i,k)=heso_cup(i,k)
718 qco(i,k)=qeso_cup(i,k)
733 if(make_calc_for_xk)
then
734 call cup_up_aa0(aa0,z,zu,dby,gamma_cup,t_cup, &
736 itf,ktf, its,ite, kts,kte)
737 call cup_up_aa0(aa1,zo,zuo,dbyo,gammao_cup,tn_cup, &
739 itf,ktf, its,ite, kts,kte)
746 ierrc(i)=
"cloud work function zero"
756 k22,kbcon,ktop,zo,entr_rate_2d,cd,fv,r_d,el2orc,qeso,tn,qo,po,dbyo, &
757 clw_all,qrco,delp,zu,wu2,omega_u,zeta,wc,omegac,zdqca)
823 dp=100.*(po_cup(i,1)-po_cup(i,2))
824 dellu(i,1)= -zuo(i,2)*(uc(i,2)-u_cup(i,2)) *g/dp
825 dellv(i,1)= -zuo(i,2)*(vc(i,2)-v_cup(i,2)) *g/dp
826 dellah(i,1)=-zuo(i,2)*(hco(i,2)-heo_cup(i,2))*g/dp
828 dellaq(i,1)=-zuo(i,2)*(qco(i,2)-qo_cup(i,2))*g/dp
832 entup=up_massentro(i,k)
833 detup=up_massdetro(i,k)
834 totmas=detup-entup+zuo(i,k+1)-zuo(i,k)
836 if(abs(totmas).gt.1.e-6)
then
837 write(0,*)
'*********************',i,k,totmas
838 write(0,*)k22(i),kbcon(i),ktop(i)
841 dp=100.*(po_cup(i,k)-po_cup(i,k+1))
842 dellah(i,k) =-(zuo(i,k+1)*(hco(i,k+1)-heo_cup(i,k+1) )- &
843 zuo(i,k )*(hco(i,k )-heo_cup(i,k ) ))*g/dp
846 dz=zo_cup(i,k+1)-zo_cup(i,k)
847 if(k.lt.ktop(i) .and. c1d(i,k) > 0)
then
848 dellaqc(i,k)= zuo(i,k)*c1d(i,k)*qrco(i,k)*dz/dp*g
850 dellaqc(i,k)=detup*0.5*(qrco(i,k+1)+qrco(i,k)) *g/dp
856 c_up = dellaqc(i,k)+(zuo(i,k+1)* qrco(i,k+1) - &
857 zuo(i,k )* qrco(i,k ) )*g/dp
861 dellaq(i,k) =-(zuo(i,k+1)*(qco(i,k+1)-qo_cup(i,k+1) ) - &
862 zuo(i,k )*(qco(i,k )-qo_cup(i,k ) ) )*g/dp &
863 - c_up - 0.5*(pwo(i,k)+pwo(i,k+1))*g/dp
864 dellu(i,k) =-(zuo(i,k+1)*(uc(i,k+1)-u_cup(i,k+1) ) - &
865 zuo(i,k )*(uc(i,k )-u_cup(i,k ) ) )*g/dp
866 dellv(i,k) =-(zuo(i,k+1)*(vc(i,k+1)-v_cup(i,k+1) ) - &
867 zuo(i,k )*(vc(i,k )-v_cup(i,k ) ) )*g/dp
883 xhe(i,k)=dellah(i,k)*mbdt+heo(i,k)
884 xq(i,k)=max(1.e-16,(dellaq(i,k)+dellaqc(i,k))*mbdt+qo(i,k))
885 dellat(i,k)=(1./cp)*(dellah(i,k)-xlv*(dellaq(i,k)))
886 xt(i,k)= (-dellaqc(i,k)*xlv/cp+dellat(i,k))*mbdt+tn(i,k)
887 xt(i,k)= max(190.,xt(i,k))
894 xhe(i,ktf)=heo(i,ktf)
902 if(make_calc_for_xk)
then
906 call cup_env(xz,xqes,xhe,xhes,xt,xq,po,z1, &
907 psur,ierr,tcrit,-1, &
913 call cup_env_clev(xt,xqes,xq,xhe,xhes,xz,po,xqes_cup,xq_cup, &
914 xhe_cup,xhes_cup,xz_cup,po_cup,gamma_cup,xt_cup,psur, &
933 x_add = xlv*zqexec(i)+cp*ztexec(i)
934 call get_cloud_bc(kte,xhe_cup(i,1:kte),xhkb(i),k22(i),x_add)
935 do k=1,start_level(i)-1
936 xhc(i,k)=xhe_cup(i,k)
948 xzu(i,1:ktf)=zuo(i,1:ktf)
950 do k=start_level(i)+1,ktop(i)
951 xhc(i,k)=(xhc(i,k-1)*xzu(i,k-1)-.5*up_massdetro(i,k-1)*xhc(i,k-1)+ &
952 up_massentro(i,k-1)*xhe(i,k-1)) / &
953 (xzu(i,k-1)-.5*up_massdetro(i,k-1)+up_massentro(i,k-1))
954 xdby(i,k)=xhc(i,k)-xhes_cup(i,k)
958 xhc(i,k)=xhes_cup(i,k)
970 flag_shallow = .true.
974 del(i,k) = delp(i,k)*0.001
985 call progsigma_calc(itf,ktf,flag_init,flag_restart,flag_shallow, &
986 flag_mid,del,tmf,qmicro,dbyo,zdqca,omega_u,zeta,xlv,dtime, &
987 forceqv_spechum,kbcon,ktop,cnvflg,betascu,betamcu,betadcu, &
988 sigmind,sigminm,sigmins,sigmain,sigmaout,sigmab)
994 call cup_up_aa0(xaa0,xz,xzu,xdby,gamma_cup,xt_cup, &
1012 xmb(i) = sigmab(i)*((-1.0*omegac(i))*gravinv)
1013 if (dx(i) < 10.e3)
then
1014 scaldfunc(i)=(1.-sigmab(i))*(1.-sigmab(i))
1015 scaldfunc(i) = max(min(scaldfunc(i), 1.0), 0.)
1019 xmb(i)=scaldfunc(i)*xmb(i)
1025 if(ierr(i).eq.0)
then
1030 xkshal=(xaa0(i)-aa1(i))/mbdt
1031 if(xkshal.le.0.and.xkshal.gt.-.01*mbdt) &
1033 if(xkshal.gt.0.and.xkshal.lt.1.e-2) &
1036 xff_shal(1)=max(0.,-(aa1(i)-aa0(i))/(xkshal*dtime))
1039 xff_shal(2)=.03*zws(i)
1044 blqe=blqe+100.*dhdt(i,k)*(po_cup(i,k)-po_cup(i,k+1))/g
1046 trash=max((hc(i,kbcon(i))-he_cup(i,kbcon(i))),1.e1)
1047 xff_shal(3)=max(0.,blqe/trash)
1048 xff_shal(3)=min(xmbmax(i),xff_shal(3))
1050 xmb(i)=(xff_shal(1)+xff_shal(2)+xff_shal(3))/3.
1051 xmb(i)=min(xmbmax(i),xmb(i))
1052 if(ichoice > 0)xmb(i)=min(xmbmax(i),xff_shal(ichoice))
1053 if(xmb(i) <= 0.)
then
1063 if(ierr(i).ne.0)
then
1073 else if(ierr(i).eq.0)
then
1081 outt(i,k)= dellat(i,k)*xmb(i)
1082 outq(i,k)= dellaq(i,k)*xmb(i)
1083 outqc(i,k)= dellaqc(i,k)*xmb(i)
1085 pre(i) = pre(i)+pwo(i,k)*xmb(i)
1087 outt(i,1)= dellat(i,1)*xmb(i)
1088 outq(i,1)= dellaq(i,1)*xmb(i)
1089 outu(i,1)=dellu(i,1)*xmb(i)
1090 outv(i,1)=dellv(i,1)*xmb(i)
1092 outu(i,k)=.25*(dellu(i,k-1)+2.*dellu(i,k)+dellu(i,k+1))*xmb(i)
1093 outv(i,k)=.25*(dellv(i,k-1)+2.*dellv(i,k)+dellv(i,k+1))*xmb(i)
1103 if(ierr(i).eq.0)
then
1107 dp=(po_cup(i,k)-po_cup(i,k+1))*100.
1109 dts= dts -(outu(i,k)*us(i,k)+outv(i,k)*vs(i,k))*dp/g
1111 fpi = fpi +sqrt(outu(i,k)*outu(i,k) + outv(i,k)*outv(i,k))*dp
1115 fp= sqrt((outu(i,k)*outu(i,k)+outv(i,k)*outv(i,k)))/fpi
1116 outt(i,k)=outt(i,k)+fp*dts*g/cp