78 itf,ktf,its,ite, kts,kte &
81 ,fv,r_d & ! ratio of vapor to dry air gas constants minus one
82 ,dicycle & ! diurnal cycle flag
83 ,ichoice & ! choice of closure, use "0" for ensemble average
84 ,ipr & ! this flag can be used for debugging prints
85 ,ccn & ! not well tested yet
87 ,dtime & ! dt over which forcing is applied
88 ,imid & ! flag to turn on mid level convection
89 ,kpbl & ! level of boundary layer height
90 ,dhdt & ! boundary layer forcing (one closure for shallow)
154 ,do_capsuppress,cap_suppress_j &
162 nranflag,itf,ktf,its,ite, kts,kte,ipr,imid
163 integer,
intent (in ) :: &
165 real(kind=kind_phys),
dimension (its:,:) &
166 ,
intent (in ) :: rand_clos
167 real(kind=kind_phys),
dimension (its:) &
168 ,
intent (in ) :: rand_mom,rand_vmas
170 real(kind=kind_phys),
intent(in),
dimension (its:),
optional :: ca_deep(:)
171 integer,
intent(in) :: do_capsuppress
172 real(kind=kind_phys),
intent(in),
dimension(:) :: cap_suppress_j
177 real(kind=kind_phys),
dimension (its:ite,1:maxens3) :: xf_ens,pr_ens
183 real(kind=kind_phys),
dimension (its:,kts:) &
184 ,
intent (inout ) :: &
185 cnvwt,outu,outv,outt,outq,outqc,cupclw
186 real(kind=kind_phys),
dimension (its:) &
189 real(kind=kind_phys),
dimension (its:,kts:) &
192 real(kind=kind_phys),
dimension (its:,kts:) &
193 ,
intent (in ),
optional :: &
194 qmicro, sigmain, forceqv_spechum
195 real(kind=kind_phys),
dimension (its:) &
196 ,
intent (inout ) :: &
199 real(kind=kind_phys),
dimension (its:) &
201 hfx,qfx,xmbm_in,xmbs_in
203 integer,
dimension (its:) &
204 ,
intent (inout ) :: &
207 integer,
dimension (its:) &
216 real(kind=kind_phys),
dimension (its:,kts:) &
218 dhdt,rho,t,po,us,vs,tn,delp
220 real(kind=kind_phys),
dimension (its:,kts:) &
221 ,
intent (inout ) :: &
224 real(kind=kind_phys),
dimension (its:,kts:) &
228 real(kind=kind_phys),
dimension (its:,kts:) &
229 ,
intent (out),
optional :: &
231 real(kind=kind_phys),
dimension (its:) &
235 real(kind=kind_phys),
dimension (its:) &
236 ,
intent (inout ) :: &
241 real(kind=kind_phys) &
243 dtime,ccnclean,fv,r_d,betascu,betamcu,betadcu
249 real(kind=kind_phys),
dimension (its:ite,1) :: &
251 real(kind=kind_phys),
dimension (its:ite,1) :: &
253 real(kind=kind_phys),
dimension (its:ite,kts:kte,1) :: &
254 dellat_ens,dellaqc_ens,dellaq_ens,pwo_ens
322 real(kind=kind_phys),
dimension (its:ite,kts:kte) :: &
323 entr_rate_2d,mentrd_rate_2d,he,hes,qes,z, heo,heso,qeso,zo, &
324 xhe,xhes,xqes,xz,xt,xq,qes_cup,q_cup,he_cup,hes_cup,z_cup, &
325 p_cup,gamma_cup,t_cup, qeso_cup,qo_cup,heo_cup,heso_cup, &
326 zo_cup,po_cup,gammao_cup,tn_cup, &
327 xqes_cup,xq_cup,xhe_cup,xhes_cup,xz_cup, &
328 xt_cup, dby,hc,zu,clw_all, &
329 dbyo,qco,qrcdo,pwdo,pwo,hcdo,qcdo,dbydo,hco,qrco, &
338 cd,cdd,dellah,dellaq,dellat,dellaqc, &
339 u_cup,v_cup,uc,vc,ucd,vcd,dellu,dellv, &
341 wu2,omega_u,zeta,zdqca,dbyo1,del
361 real(kind=kind_phys),
dimension (its:ite) :: &
362 edt,edto,edtm,aa1,aa0,xaa0,hkb, &
365 pwevo,bu,bud,cap_max,wc,omegac,sigmab, &
366 cap_max_increment,closure_n,psum,psumh,sig,sigd
367 real(kind=kind_phys),
dimension (its:ite) :: &
368 axx,edtmax,edtmin,entr_rate
369 integer,
dimension (its:ite) :: &
370 kzdown,kdet,k22,jmin,kstabi,kstabm,k22x,xland1, &
371 ktopdby,kbconx,ierr2,ierr3,kbmax
381 integer,
dimension (its:),
intent(inout) :: ierr
382 integer,
dimension (its:),
intent(in),
optional :: csum
383 logical,
intent(in) :: do_ca, progsigma
384 logical,
intent(in) :: flag_init, flag_restart
387 iloop,nens3,ki,kk,i,k
388 real(kind=kind_phys) :: &
389 dz,dzo,mbdt,radius, &
390 zcutdown,depth_min,zkbmax,z_detr,zktop, &
391 dh,cap_maxs,trash,trash2,frh,sig_thresh
392 real(kind=kind_phys),
dimension (its:ite) :: pefc
393 real(kind=kind_phys) entdo,dp,subin,detdo,entup, &
394 detup,subdown,entdoj,entupk,detupk,totmas
395 real(kind=kind_phys) :: &
396 sigmind,sigminm,sigmins
397 parameter(sigmind=0.005,sigmins=0.03,sigminm=0.01)
399 real(kind=kind_phys),
dimension (its:ite) :: lambau,flux_tun,zws,ztexec,zqexec
402 integer :: jprnt,jmini,start_k22
403 logical :: keep_going,flg(its:ite),cnvflg(its:ite)
404 logical :: flag_shallow,flag_mid
408 character*50 :: ierrc(its:ite)
409 character*4 :: cumulus
410 real(kind=kind_phys),
dimension (its:ite,kts:kte) :: &
411 up_massentr,up_massdetr,c1d &
412 ,up_massentro,up_massdetro,dd_massentro,dd_massdetro
413 real(kind=kind_phys),
dimension (its:ite,kts:kte) :: &
414 up_massentru,up_massdetru,dd_massentru,dd_massdetru
417 real(kind=kind_phys) c1_max,buo_flux,pgcon,pgc,blqe
419 real(kind=kind_phys) :: xff_mid(its:ite,2)
421 integer :: iversion=1
422 real(kind=kind_phys) :: denom,h_entr,umean,t_star,dq
423 integer,
intent(in) :: dicycle
424 real(kind=kind_phys),
dimension (its:ite) :: aa1_bl,hkbo_bl,tau_bl,tau_ecmwf,wmean
425 real(kind=kind_phys),
dimension (its:ite,kts:kte) :: tn_bl, qo_bl, qeso_bl, heo_bl, heso_bl &
426 ,qeso_cup_bl,qo_cup_bl, heo_cup_bl,heso_cup_bl &
427 ,gammao_cup_bl,tn_cup_bl,hco_bl,dbyo_bl
428 real(kind=kind_phys),
dimension(its:ite) :: xf_dicycle,xf_progsigma
433 real(kind=kind_phys),
intent(inout),
dimension(its:,:) :: forcing
435 integer :: turn,pmin_lev(its:ite),start_level(its:ite),ktopkeep(its:ite)
436 real(kind=kind_phys),
dimension (its:ite,kts:kte) :: dtempdz
437 integer,
dimension (its:ite,kts:kte) :: k_inv_layers
438 real(kind=kind_phys),
dimension (its:ite) :: c0
442 real(kind=kind_phys) zuh2(40)
443 real(kind=kind_phys),
dimension (its:ite) :: rntot,delqev,delq2,qevap,rn,qcond
445 real(kind=kind_phys) :: rain,t1,q1,elocp,evef,el2orc,evfact,evfactl,g_rain,e_dn,c_up
446 real(kind=kind_phys) :: pgeoh,dts,fp,fpi,pmin,x_add,beta,beta_u
447 real(kind=kind_phys) :: cbeg,cmid,cend,const_a,const_b,const_c
450 real(kind=kind_phys),
dimension (its:ite,kts:kte) :: p_liq_ice,melting_layer,melting
457 melting_layer(:,:)=0.
463 if(imid.eq.1)cumulus=
'mid'
465 if(imid.eq.1)pmin=75.
471 el2orc=xlv*xlv/(r_v*cp)
484 if(imid.eq.1)lambau(:)=2.0
486 if(nranflag == 1)
then
487 lambau(:)=1.5+rand_mom(:)
499 xland1(i)=int(xland(i)+.0001)
500 if(xland(i).gt.1.5 .or. xland(i).lt.0.5)
then
503 if(xland1(i).eq.1)c0(i)=0.002
518 buo_flux= (hfx(i)/cp+0.608*t(i,1)*qfx(i)/xlv)/rho(i,1)
521 zws(i) = max(0.,flux_tun(i)*0.41*buo_flux*zo(i,2)*g/t(i,1))
522 if(zws(i) > tiny(pgeoh))
then
524 zws(i) = 1.2*zws(i)**.3333
526 ztexec(i) = max(flux_tun(i)*hfx(i)/(rho(i,1)*zws(i)*cp),0.0)
528 zqexec(i) = max(flux_tun(i)*qfx(i)/xlv/(rho(i,1)*zws(i)),0.)
532 zws(i) = max(0.,.001-flux_tun(i)*0.41*buo_flux*zo(i,kpbl(i))*g/t(i,kpbl(i)))
533 zws(i) = 1.2*zws(i)**.3333
534 zws(i) = zws(i)*rho(i,kpbl(i))
544 cap_max_increment(i)=20.
548 if (xland1(i)==0)
then
549 cap_max_increment(i)=20.
551 if(ztexec(i).gt.0.)cap_max(i)=cap_max(i)+25.
552 if(ztexec(i).lt.0.)cap_max(i)=cap_max(i)-25.
559 if(use_excess == 0 )
then
565 if(do_capsuppress == 1)
then
569 if (abs(cap_suppress_j(i) - 1.0 ) < 0.1 )
then
570 cap_max(i)=cap_maxs+75.
571 elseif (abs(cap_suppress_j(i) - 0.0 ) < 0.1 )
then
589 entr_rate(i)=7.e-5 - min(20.,float(csum(i))) * 3.e-6
590 if(xland1(i) == 0)entr_rate(i)=7.e-5
591 if(dx(i)<dx_thresh) entr_rate(i)=2.e-4
592 if(imid.eq.1)entr_rate(i)=3.e-4
593 radius=.2/entr_rate(i)
594 frh=min(1.,3.14*radius*radius/dx(i)/dx(i))
595 if(frh > frh_thresh)
then
597 radius=sqrt(frh*dx(i)*dx(i)/3.14)
598 entr_rate(i)=.2/radius
602 if(forcing(i,7).eq.0.)sig(i)=1.
603 frh_out(i) = frh*sig(i)
606 sig_thresh = (1.-frh_thresh)**2
624 cd(i,k)=.1*entr_rate(i)
626 if(imid.eq.1)cd(i,k)=.5*entr_rate(i)
648 if(dx(its)<dx_thresh)depth_min=5000.
649 if(imid.eq.1)depth_min=2500.
670 if(imid.eq.1)zkbmax=2000.
695 call cup_env(z,qes,he,hes,t,q,po,z1, &
696 psur,ierr,tcrit,-1, &
699 call cup_env(zo,qeso,heo,heso,tn,qo,po,z1, &
700 psur,ierr,tcrit,-1, &
707 call cup_env_clev(t,qes,q,he,hes,z,po,qes_cup,q_cup,he_cup, &
708 hes_cup,z_cup,p_cup,gamma_cup,t_cup,psur, &
712 call cup_env_clev(tn,qeso,qo,heo,heso,zo,po,qeso_cup,qo_cup, &
713 heo_cup,heso_cup,zo_cup,po_cup,gammao_cup,tn_cup,psur, &
720 itf,ktf,its,ite,kts,kte,cumulus)
725 if(kpbl(i).gt.5 .and. imid.eq.1)cap_max(i)=po_cup(i,kpbl(i))
726 u_cup(i,kts)=us(i,kts)
727 v_cup(i,kts)=vs(i,kts)
729 u_cup(i,k)=.5*(us(i,k-1)+us(i,k))
730 v_cup(i,k)=.5*(vs(i,k-1)+vs(i,k))
737 if(zo_cup(i,k).gt.zkbmax+z1(i))
then
747 if(zo_cup(i,k).gt.z_detr+z1(i))
then
767 k22(i)=maxloc(heo_cup(i,start_k22:kbmax(i)+2),1)+start_k22-1
768 if(k22(i).ge.kbmax(i))
then
771 ierrc(i)=
"could not find k22"
788 x_add = xlv*zqexec(i)+cp*ztexec(i)
789 call get_cloud_bc(kte,he_cup(i,1:kte),hkb(i),k22(i),x_add)
790 call get_cloud_bc(kte,heo_cup(i,1:kte),hkbo(i),k22(i),x_add)
798 call cup_kbcon(ierrc,cap_max_increment,iloop,k22,kbcon,heo_cup,heso_cup, &
799 hkbo,ierr,kbmax,po_cup,cap_max, &
803 z_cup,entr_rate,heo,imid)
807 call cup_minimi(heso_cup,kbcon,kstabm,kstabi,ierr, &
813 frh = min(qo_cup(i,kbcon(i))/qeso_cup(i,kbcon(i)),1.)
814 if(frh >= rh_thresh .and. sig(i) <= sig_thresh )
then
824 if(po(i,kbcon(i))-po(i,k) > pmin+x_add)
then
832 start_level(i)=k22(i)
833 x_add = xlv*zqexec(i)+cp*ztexec(i)
834 call get_cloud_bc(kte,he_cup(i,1:kte),hkb(i),k22(i),x_add)
844 kbcon,kstabi,dtempdz,itf,ktf,its,ite, kts,kte)
848 if(kstabi(i).lt.kbcon(i))
then
853 entr_rate_2d(i,k)=entr_rate(i)
856 kbcon(i)=max(2,kbcon(i))
858 frh = min(qo_cup(i,k)/qeso_cup(i,k),1.)
859 entr_rate_2d(i,k)=entr_rate(i) *(1.3-frh)
862 if(k_inv_layers(i,2).gt.0 .and. &
863 (po_cup(i,k22(i))-po_cup(i,k_inv_layers(i,2))).lt.500.)
then
865 ktop(i)=min(kstabi(i),k_inv_layers(i,2))
870 if((po_cup(i,k22(i))-po_cup(i,k)).gt.500.)
then
890 call rates_up_pdf(rand_vmas,ipr,
'mid',ktop,ierr,po_cup,entr_rate_2d,hkbo,heo,heso_cup,zo_cup, &
891 xland1,kstabi,k22,kbcon,its,ite,itf,kts,kte,ktf,zuo,kpbl,ktopdby,csum,pmin_lev)
893 call rates_up_pdf(rand_vmas,ipr,
'deep',ktop,ierr,po_cup,entr_rate_2d,hkbo,heo,heso_cup,zo_cup, &
894 xland1,kstabi,k22,kbcon,its,ite,itf,kts,kte,ktf,zuo,kbcon,ktopdby,csum,pmin_lev)
930 ,ierr,ktop,zo_cup,zuo,cd,entr_rate_2d &
931 ,up_massentro, up_massdetro ,up_massentr, up_massdetr &
932 ,3,kbcon,k22,up_massentru,up_massdetru,lambau)
935 ,ierr,ktop,zo_cup,zuo,cd,entr_rate_2d &
936 ,up_massentro, up_massdetro ,up_massentr, up_massdetr &
937 ,1,kbcon,k22,up_massentru,up_massdetru,lambau)
958 do k=1,start_level(i)
962 do k=1,start_level(i)-1
964 hco(i,k)=heo_cup(i,k)
981 if(ierr(i) /= 0) cycle
984 do k=start_level(i) +1,ktop(i)
986 denom=zuo(i,k-1)-.5*up_massdetro(i,k-1)+up_massentro(i,k-1)
987 if(denom.lt.1.e-8)
then
991 hco(i,k)=(hco(i,k-1)*zuo(i,k-1)-.5*up_massdetro(i,k-1)*hco(i,k-1)+ &
992 up_massentro(i,k-1)*heo(i,k-1)) / &
993 (zuo(i,k-1)-.5*up_massdetro(i,k-1)+up_massentro(i,k-1))
994 dbyo(i,k)=hco(i,k)-heso_cup(i,k)
998 do k=ktop(i)-1,kbcon(i),-1
999 if(dbyo(i,k).gt.0.)
then
1010 if(ierr(i).eq.0)
then
1011 zktop=(zo_cup(i,ktop(i))-z1(i))*.6
1012 if(imid.eq.1)zktop=(zo_cup(i,ktop(i))-z1(i))*.4
1013 zktop=min(zktop+z1(i),zcutdown+z1(i))
1016 if(zo_cup(i,k).gt.zktop)
then
1018 kzdown(i)=min(kzdown(i),kstabi(i)-1)
1029 call cup_minimi(heso_cup,k22,kzdown,jmin,ierr, &
1034 if(ierr(i).eq.0)
then
1045 do while ( keep_going )
1046 keep_going = .false.
1047 if ( jmini - 1 .lt. kdet(i) ) kdet(i) = jmini-1
1048 if ( jmini .ge. ktop(i)-1 ) jmini = ktop(i) - 2
1050 hcdo(i,ki)=heso_cup(i,ki)
1051 dz=zo_cup(i,ki+1)-zo_cup(i,ki)
1055 hcdo(i,k)=heso_cup(i,jmini)
1056 dz=zo_cup(i,k+1)-zo_cup(i,k)
1057 dh=dh+dz*(hcdo(i,k)-heso_cup(i,k))
1060 if ( jmini .gt. 5 )
then
1065 ierrc(i) =
"could not find jmini9"
1073 if ( jmini .le. 5 )
then
1076 ierrc(i) =
"could not find jmini4"
1082 if(ierr(i) /= 0) cycle
1085 hco(i,k)=heso_cup(i,k)
1095 p_cup,kbcon,ktop,dbyo,clw_all,xland1, &
1096 qo,gammao_cup,zuo,qeso_cup,k22,qo_cup,c0, &
1097 zqexec,ccn,ccnclean,rho,c1d,tn_cup,autoconv,up_massentr,up_massdetr,psum,psumh, &
1102 p_cup,kbcon,ktop,dbyo,clw_all,xland1, &
1103 qo,gammao_cup,zuo,qeso_cup,k22,qo_cup,c0, &
1104 zqexec,ccn,ccnclean,rho,c1d,tn_cup,autoconv,up_massentr,up_massdetr,psum,psumh, &
1115 if(ierr(i) /= 0) cycle
1118 do k=start_level(i) +1,ktop(i)
1120 denom=zuo(i,k-1)-.5*up_massdetro(i,k-1)+up_massentro(i,k-1)
1121 if(denom.lt.1.e-8)
then
1126 hc(i,k)=(hc(i,k-1)*zu(i,k-1)-.5*up_massdetr(i,k-1)*hc(i,k-1)+ &
1127 up_massentr(i,k-1)*he(i,k-1)) / &
1128 (zu(i,k-1)-.5*up_massdetr(i,k-1)+up_massentr(i,k-1))
1129 uc(i,k)=(uc(i,k-1)*zu(i,k-1)-.5*up_massdetru(i,k-1)*uc(i,k-1)+ &
1130 up_massentru(i,k-1)*us(i,k-1) &
1131 -pgcon*.5*(zu(i,k)+zu(i,k-1))*(u_cup(i,k)-u_cup(i,k-1))) / &
1132 (zu(i,k-1)-.5*up_massdetru(i,k-1)+up_massentru(i,k-1))
1133 vc(i,k)=(vc(i,k-1)*zu(i,k-1)-.5*up_massdetru(i,k-1)*vc(i,k-1)+ &
1134 up_massentru(i,k-1)*vs(i,k-1) &
1135 -pgcon*.5*(zu(i,k)+zu(i,k-1))*(v_cup(i,k)-v_cup(i,k-1))) / &
1136 (zu(i,k-1)-.5*up_massdetru(i,k-1)+up_massentru(i,k-1))
1137 dby(i,k)=hc(i,k)-hes_cup(i,k)
1138 hco(i,k)=(hco(i,k-1)*zuo(i,k-1)-.5*up_massdetro(i,k-1)*hco(i,k-1)+ &
1139 up_massentro(i,k-1)*heo(i,k-1)) / &
1140 (zuo(i,k-1)-.5*up_massdetro(i,k-1)+up_massentro(i,k-1))
1145 hc(i,k)= hc(i,k)+(1.-p_liq_ice(i,k))*qrco(i,k)*xlf
1146 hco(i,k)= hco(i,k)+(1.-p_liq_ice(i,k))*qrco(i,k)*xlf
1148 dby(i,k)=hc(i,k)-hes_cup(i,k)
1150 dbyo(i,k)=hco(i,k)-heso_cup(i,k)
1151 dz=zo_cup(i,k+1)-zo_cup(i,k)
1152 dbyt(i,k)=dbyt(i,k-1)+dbyo(i,k)*dz
1156 kk=maxloc(dbyt(i,:),1)
1157 ki=maxloc(zuo(i,:),1)
1160 do k=ktop(i)-1,kbcon(i),-1
1161 if(dbyo(i,k).gt.0.)
then
1172 if(ierr(i) /= 0) cycle
1174 hc(i,k)=hes_cup(i,k)
1177 hco(i,k)=heso_cup(i,k)
1183 entr_rate_2d(i,k)=0.
1186 up_massentro(i,k)=0.
1187 up_massdetro(i,k)=0.
1193 if(ktop(i).lt.kbcon(i)+2)
then
1196 ierrc(i)=
'ktop too small deep'
1208 if(ierr(i).eq.0)
then
1209 if ( jmin(i) - 1 .lt. kdet(i) ) kdet(i) = jmin(i)-1
1210 if(-zo_cup(i,kbcon(i))+zo_cup(i,ktop(i)).lt.depth_min)
then
1213 ierrc(i)=
"cloud depth very shallow"
1229 dd_massentro(i,k)=0.
1230 dd_massdetro(i,k)=0.
1231 dd_massentru(i,k)=0.
1232 dd_massdetru(i,k)=0.
1233 hcdo(i,k)=heso_cup(i,k)
1237 mentrd_rate_2d(i,k)=entr_rate(i)
1245 beta=max(.025,.055-float(csum(i))*.0015)
1246 if(imid.eq.1)beta=.025
1248 cdd(i,1:jmin(i))=.1*entr_rate(i)
1250 dd_massdetro(i,:)=0.
1251 dd_massentro(i,:)=0.
1252 call get_zu_zd_pdf_fim(0,po_cup(i,:),rand_vmas(i),0.0_kind_phys,ipr,xland1(i),zuh2,4, &
1253 ierr(i),kdet(i),jmin(i)+1,zdo(i,:),kts,kte,ktf,beta,kpbl(i),csum(i),pmin_lev(i))
1254 if(zdo(i,jmin(i)) .lt.1.e-8)
then
1257 cdd(i,jmin(i):ktf)=0.
1258 zdo(i,jmin(i)+1:ktf)=0.
1259 if(zdo(i,jmin(i)) .lt.1.e-8)
then
1265 itemp = maxloc(zdo(i,:),1)
1266 do ki=jmin(i) , itemp,-1
1268 dzo=zo_cup(i,ki+1)-zo_cup(i,ki)
1269 dd_massdetro(i,ki)=cdd(i,ki)*dzo*zdo(i,ki+1)
1270 dd_massentro(i,ki)=zdo(i,ki)-zdo(i,ki+1)+dd_massdetro(i,ki)
1271 if(dd_massentro(i,ki).lt.0.)
then
1272 dd_massentro(i,ki)=0.
1273 dd_massdetro(i,ki)=zdo(i,ki+1)-zdo(i,ki)
1274 if(zdo(i,ki+1).gt.0.)cdd(i,ki)=dd_massdetro(i,ki)/(dzo*zdo(i,ki+1))
1276 if(zdo(i,ki+1).gt.0.)mentrd_rate_2d(i,ki)=dd_massentro(i,ki)/(dzo*zdo(i,ki+1))
1278 mentrd_rate_2d(i,1)=0.
1281 dzo=zo_cup(i,ki+1)-zo_cup(i,ki)
1282 dd_massentro(i,ki)=mentrd_rate_2d(i,ki)*dzo*zdo(i,ki+1)
1283 dd_massdetro(i,ki) = zdo(i,ki+1)+dd_massentro(i,ki)-zdo(i,ki)
1284 if(dd_massdetro(i,ki).lt.0.)
then
1285 dd_massdetro(i,ki)=0.
1286 dd_massentro(i,ki)=zdo(i,ki)-zdo(i,ki+1)
1287 if(zdo(i,ki+1).gt.0.)mentrd_rate_2d(i,ki)=dd_massentro(i,ki)/(dzo*zdo(i,ki+1))
1289 if(zdo(i,ki+1).gt.0.)cdd(i,ki)= dd_massdetro(i,ki)/(dzo*zdo(i,ki+1))
1294 dd_massentru(i,k-1)=dd_massentro(i,k-1)+lambau(i)*dd_massdetro(i,k-1)
1295 dd_massdetru(i,k-1)=dd_massdetro(i,k-1)+lambau(i)*dd_massdetro(i,k-1)
1297 dbydo(i,jmin(i))=hcdo(i,jmin(i))-heso_cup(i,jmin(i))
1298 bud(i)=dbydo(i,jmin(i))*(zo_cup(i,jmin(i)+1)-zo_cup(i,jmin(i)))
1299 ucd(i,jmin(i)+1)=.5*(uc(i,jmin(i)+1)+u_cup(i,jmin(i)+1))
1302 dzo=zo_cup(i,ki+1)-zo_cup(i,ki)
1303 h_entr=.5*(heo(i,ki)+.5*(hco(i,ki)+hco(i,ki+1)))
1304 ucd(i,ki)=(ucd(i,ki+1)*zdo(i,ki+1) &
1305 -.5*dd_massdetru(i,ki)*ucd(i,ki+1)+ &
1306 dd_massentru(i,ki)*us(i,ki) &
1307 -pgcon*zdo(i,ki+1)*(us(i,ki+1)-us(i,ki))) / &
1308 (zdo(i,ki+1)-.5*dd_massdetru(i,ki)+dd_massentru(i,ki))
1309 vcd(i,ki)=(vcd(i,ki+1)*zdo(i,ki+1) &
1310 -.5*dd_massdetru(i,ki)*vcd(i,ki+1)+ &
1311 dd_massentru(i,ki)*vs(i,ki) &
1312 -pgcon*zdo(i,ki+1)*(vs(i,ki+1)-vs(i,ki))) / &
1313 (zdo(i,ki+1)-.5*dd_massdetru(i,ki)+dd_massentru(i,ki))
1314 hcdo(i,ki)=(hcdo(i,ki+1)*zdo(i,ki+1) &
1315 -.5*dd_massdetro(i,ki)*hcdo(i,ki+1)+ &
1316 dd_massentro(i,ki)*h_entr) / &
1317 (zdo(i,ki+1)-.5*dd_massdetro(i,ki)+dd_massentro(i,ki))
1318 dbydo(i,ki)=hcdo(i,ki)-heso_cup(i,ki)
1319 bud(i)=bud(i)+dbydo(i,ki)*dzo
1325 ierrc(i)=
'downdraft is not negatively buoyant '
1335 pwdo,qo_cup,zo_cup,dd_massentro,dd_massdetro,jmin,ierr,gammao_cup, &
1336 pwevo,bu,qrcdo,po_cup,qo,heo,1, &
1344 dp=100.*(po_cup(i,1)-po_cup(i,2))
1345 cupclw(i,k)=qrco(i,k)
1346 cnvwt(i,k)=zuo(i,k)*cupclw(i,k)*g/dp
1355 if(k > kbcon(i) .and. k < ktop(i))
then
1356 dbyo1(i,k)=hco(i,k)-heso_cup(i,k)
1365 call cup_up_aa0(aa0,z,zu,dby,gamma_cup,t_cup, &
1369 call cup_up_aa0(aa1,zo,zuo,dbyo,gammao_cup,tn_cup, &
1377 if(aa1(i).eq.0.)
then
1380 ierrc(i)=
"cloud work function zero"
1389 k22,kbcon,ktop,zo,entr_rate_2d,cd,fv,r_d,el2orc,qeso,tn,qo,po,dbyo, &
1390 clw_all,qrco,delp,zu,wu2,omega_u,zeta,wc,omegac,zdqca)
1411 if(ierr(i).eq.0)
then
1414 if(imid.eq.1)wmean(i) = 3.0
1416 tau_ecmwf(i)=( zo_cup(i,ktop(i))- zo_cup(i,kbcon(i)) ) / wmean(i)
1417 tau_ecmwf(i)=max(tau_ecmwf(i),720.)
1418 tau_ecmwf(i)= tau_ecmwf(i) * (1.0061 + 1.23e-2 * (dx(i)/1000.))
1425 if(dicycle == 1)
then
1429 if(ierr(i).eq.0)
then
1430 if(xland1(i) == 0 )
then
1432 umean= 2.0+sqrt(0.5*(us(i,1)**2+vs(i,1)**2+us(i,kbcon(i))**2+vs(i,kbcon(i))**2))
1433 tau_bl(i) = (zo_cup(i,kbcon(i))- z1(i)) /umean
1436 tau_bl(i) =( zo_cup(i,ktopdby(i))- zo_cup(i,kbcon(i)) ) / wmean(i)
1445 tn_bl(i,:)=0.;qo_bl(i,:)=0.
1446 if ( ierr(i) == 0 )
then
1448 tn_bl(i,1:kbcon(i)) = tn(i,1:kbcon(i))
1449 qo_bl(i,1:kbcon(i)) = qo(i,1:kbcon(i))
1451 tn_bl(i,kbcon(i)+1:ktf) = t(i,kbcon(i)+1:ktf)
1452 qo_bl(i,kbcon(i)+1:ktf) = q(i,kbcon(i)+1:ktf)
1457 call cup_env(zo,qeso_bl,heo_bl,heso_bl,tn_bl,qo_bl,po,z1, &
1458 psur,ierr,tcrit,-1, &
1459 itf,ktf, its,ite, kts,kte)
1461 call cup_env_clev(tn_bl,qeso_bl,qo_bl,heo_bl,heso_bl,zo,po,qeso_cup_bl,qo_cup_bl, &
1462 heo_cup_bl,heso_cup_bl,zo_cup,po_cup,gammao_cup_bl,tn_cup_bl,psur,&
1464 itf,ktf,its,ite, kts,kte)
1466 if(iversion == 1)
then
1473 zo_cup,zuo,dbyo_bl,gammao_cup_bl,tn_cup_bl, &
1475 itf,ktf,its,ite, kts,kte)
1479 if(ierr(i).eq.0)
then
1487 aa1_bl(i) = ( aa1_bl(i)/t_star)* tau_bl(i)
1499 if(ierr(i).eq.0)
then
1500 hkbo_bl(i)=heo_cup_bl(i,k22(i))
1510 if(ierr(i).eq.0)
then
1512 hco_bl(i,k)=hkbo_bl(i)
1515 hco_bl(i,k)=hkbo_bl(i)
1516 dbyo_bl(i,k)=hkbo_bl(i) - heso_cup_bl(i,k)
1522 if(ierr(i).eq.0)
then
1523 do k=kbcon(i)+1,ktop(i)
1524 hco_bl(i,k)=(hco_bl(i,k-1)*zuo(i,k-1)-.5*up_massdetro(i,k-1)*hco_bl(i,k-1)+ &
1525 up_massentro(i,k-1)*heo_bl(i,k-1)) / &
1526 (zuo(i,k-1)-.5*up_massdetro(i,k-1)+up_massentro(i,k-1))
1527 dbyo_bl(i,k)=hco_bl(i,k)-heso_cup_bl(i,k)
1530 hco_bl(i,k)=heso_cup_bl(i,k)
1537 call cup_up_aa0(aa1_bl,zo,zuo,dbyo_bl,gammao_cup_bl,tn_cup_bl, &
1539 itf,ktf,its,ite, kts,kte)
1543 if(ierr(i).eq.0)
then
1545 aa1_bl(i) = aa1_bl(i) - aa0(i)
1551 aa1_bl(i) = aa1_bl(i)* tau_bl(i)/ dtime
1569 call cup_dd_edt(ierr,us,vs,zo,ktop,kbcon,edt,po,pwavo, &
1570 pwo,ccn,ccnclean,pwevo,edtmax,edtmin,edtc,psum,psumh, &
1571 rho,aeroevap,pefc,xland1,itf,ktf, &
1580 ,pwo,edto,pwdo,melting &
1581 ,itf,ktf,its,ite,kts,kte,cumulus)
1585 dellat_ens(i,k,1)=0.
1586 dellaq_ens(i,k,1)=0.
1587 dellaqc_ens(i,k,1)=0.
1649 dp=100.*(po_cup(i,1)-po_cup(i,2))
1650 dellu(i,1)=pgcd*(edto(i)*zdo(i,2)*ucd(i,2) &
1651 -edto(i)*zdo(i,2)*u_cup(i,2))*g/dp &
1652 -zuo(i,2)*(uc(i,2)-u_cup(i,2)) *g/dp
1653 dellv(i,1)=pgcd*(edto(i)*zdo(i,2)*vcd(i,2) &
1654 -edto(i)*zdo(i,2)*v_cup(i,2))*g/dp &
1655 -zuo(i,2)*(vc(i,2)-v_cup(i,2)) *g/dp
1661 if(k == k22(i)-1) entupk=zuo(i,k+1)
1665 detdo=edto(i)*dd_massdetro(i,k)
1666 entdo=edto(i)*dd_massentro(i,k)
1668 entup=up_massentro(i,k)
1669 detup=up_massdetro(i,k)
1671 subin=-zdo(i,k+1)*edto(i)
1672 subdown=-zdo(i,k)*edto(i)
1674 if(k.eq.ktop(i))
then
1675 detupk=zuo(i,ktop(i))
1683 totmas=subin-subdown+detup-entup-entdo+ &
1684 detdo-entupk-entdoj+detupk+zuo(i,k+1)-zuo(i,k)
1685 if(abs(totmas).gt.1.e-6)
then
1687 write(0,123)
'totmas=',k22(i),kbcon(i),k,entup,detup,edto(i),zdo(i,k+1),dd_massdetro(i,k),dd_massentro(i,k)
1688123
format(a7,1x,3i3,2e12.4,1(1x,f5.2),3e12.4)
1691 dp=100.*(po_cup(i,k)-po_cup(i,k+1))
1693 if(k.ge.ktop(i))pgc=0.
1695 dellu(i,k) =-(zuo(i,k+1)*(uc(i,k+1)-u_cup(i,k+1) ) - &
1696 zuo(i,k )*(uc(i,k )-u_cup(i,k ) ) )*g/dp &
1697 +(zdo(i,k+1)*(ucd(i,k+1)-u_cup(i,k+1) ) - &
1698 zdo(i,k )*(ucd(i,k )-u_cup(i,k ) ) )*g/dp*edto(i)*pgcd
1699 dellv(i,k) =-(zuo(i,k+1)*(vc(i,k+1)-v_cup(i,k+1) ) - &
1700 zuo(i,k )*(vc(i,k )-v_cup(i,k ) ) )*g/dp &
1701 +(zdo(i,k+1)*(vcd(i,k+1)-v_cup(i,k+1) ) - &
1702 zdo(i,k )*(vcd(i,k )-v_cup(i,k ) ) )*g/dp*edto(i)*pgcd
1711 if(ierr(i).eq.0)
then
1713 dp=100.*(po_cup(i,1)-po_cup(i,2))
1715 dellah(i,1)=(edto(i)*zdo(i,2)*hcdo(i,2) &
1716 -edto(i)*zdo(i,2)*heo_cup(i,2))*g/dp &
1717 -zuo(i,2)*(hco(i,2)-heo_cup(i,2))*g/dp
1719 dellaq(i,1)=(edto(i)*zdo(i,2)*qcdo(i,2) &
1720 -edto(i)*zdo(i,2)*qo_cup(i,2))*g/dp &
1721 -zuo(i,2)*(qco(i,2)-qo_cup(i,2))*g/dp
1723 g_rain= 0.5*(pwo(i,1)+pwo(i,2))*g/dp
1724 e_dn = -0.5*(pwdo(i,1)+pwdo(i,2))*g/dp*edto(i)
1725 dellaq(i,1) = dellaq(i,1)+ e_dn-g_rain
1735 dp=100.*(po_cup(i,k)-po_cup(i,k+1))
1738 dellah(i,k) =-(zuo(i,k+1)*(hco(i,k+1)-heo_cup(i,k+1) ) - &
1739 zuo(i,k )*(hco(i,k )-heo_cup(i,k ) ) )*g/dp &
1740 +(zdo(i,k+1)*(hcdo(i,k+1)-heo_cup(i,k+1) ) - &
1741 zdo(i,k )*(hcdo(i,k )-heo_cup(i,k ) ) )*g/dp*edto(i)
1745 dellah(i,k) = dellah(i,k) + xlf*((1.-p_liq_ice(i,k))*0.5*(qrco(i,k+1)+qrco(i,k)) &
1746 - melting(i,k))*g/dp
1755 detup=up_massdetro(i,k)
1756 dz=zo_cup(i,k)-zo_cup(i,k-1)
1757 if(k.lt.ktop(i))
then
1758 dellaqc(i,k) = zuo(i,k)*c1d(i,k)*qrco(i,k)*dz/dp*g
1760 dellaqc(i,k)= detup*0.5*(qrco(i,k+1)+qrco(i,k)) *g/dp
1763 g_rain= 0.5*(pwo(i,k)+pwo(i,k+1))*g/dp
1764 e_dn = -0.5*(pwdo(i,k)+pwdo(i,k+1))*g/dp*edto(i)
1768 c_up = dellaqc(i,k)+(zuo(i,k+1)* qrco(i,k+1) - &
1769 zuo(i,k )* qrco(i,k ) )*g/dp + g_rain
1774 dellaq(i,k) =-(zuo(i,k+1)*(qco(i,k+1)-qo_cup(i,k+1) ) - &
1775 zuo(i,k )*(qco(i,k )-qo_cup(i,k ) ) )*g/dp &
1776 +(zdo(i,k+1)*(qcdo(i,k+1)-qo_cup(i,k+1) ) - &
1777 zdo(i,k )*(qcdo(i,k )-qo_cup(i,k ) ) )*g/dp*edto(i) &
1788444
format(1x,i2,1x,7e12.4)
1799 if(ierr(i).eq.0)
then
1801 xhe(i,k)=dellah(i,k)*mbdt+heo(i,k)
1803 xq(i,k)=max(1.e-16,dellaq(i,k)*mbdt+qo(i,k))
1804 dellat(i,k)=(1./cp)*(dellah(i,k)-xlv*dellaq(i,k))
1806 xt(i,k)= dellat(i,k)*mbdt+tn(i,k)
1807 xt(i,k)=max(190.,xt(i,k))
1812 xt(i,k)=tn(i,k)+0.25*(dellat(i,k-1) + 2.*dellat(i,k) + dellat(i,k+1)) * mbdt
1813 xt(i,k)=max(190.,xt(i,k))
1814 xq(i,k)=max(1.e-16, qo(i,k)+0.25*(dellaq(i,k-1) + 2.*dellaq(i,k) + dellaq(i,k+1)) * mbdt)
1815 xhe(i,k)=heo(i,k)+0.25*(dellah(i,k-1) + 2.*dellah(i,k) + dellah(i,k+1)) * mbdt
1820 if(ierr(i).eq.0)
then
1821 xhe(i,ktf)=heo(i,ktf)
1830 call cup_env(xz,xqes,xhe,xhes,xt,xq,po,z1, &
1831 psur,ierr,tcrit,-1, &
1837 call cup_env_clev(xt,xqes,xq,xhe,xhes,xz,po,xqes_cup,xq_cup, &
1838 xhe_cup,xhes_cup,xz_cup,po_cup,gamma_cup,xt_cup,psur, &
1858 if(ierr(i).eq.0)
then
1859 x_add = xlv*zqexec(i)+cp*ztexec(i)
1860 call get_cloud_bc(kte,xhe_cup(i,1:kte),xhkb(i),k22(i),x_add)
1861 do k=1,start_level(i)-1
1862 xhc(i,k)=xhe_cup(i,k)
1873 if(ierr(i).eq.0)
then
1875 do k=start_level(i) +1,ktop(i)
1876 xhc(i,k)=(xhc(i,k-1)*xzu(i,k-1)-.5*up_massdetro(i,k-1)*xhc(i,k-1) + &
1877 up_massentro(i,k-1)*xhe(i,k-1)) / &
1878 (xzu(i,k-1)-.5*up_massdetro(i,k-1)+up_massentro(i,k-1))
1885 xhc(i,k)= xhc(i,k)+ xlf*(1.-p_liq_ice(i,k))*qrco(i,k)
1888 xdby(i,k)=xhc(i,k)-xhes_cup(i,k)
1892 xhc(i,k)=xhes_cup(i,k)
1901 call cup_up_aa0(xaa0,xz,xzu,xdby,gamma_cup,xt_cup, &
1907 if(ierr(i).eq.0)
then
1908 xaa0_ens(i,1)=xaa0(i)
1915 pr_ens(i,nens3)=pr_ens(i,nens3) &
1916 +pwo(i,k)+edto(i)*pwdo(i,k)
1918 else if(nens3.eq.8)
then
1919 pr_ens(i,nens3)=pr_ens(i,nens3)+ &
1920 pwo(i,k)+edto(i)*pwdo(i,k)
1922 else if(nens3.eq.9)
then
1923 pr_ens(i,nens3)=pr_ens(i,nens3) &
1924 + pwo(i,k)+edto(i)*pwdo(i,k)
1926 pr_ens(i,nens3)=pr_ens(i,nens3)+ &
1927 pwo(i,k) +edto(i)*pwdo(i,k)
1931 if(pr_ens(i,7).lt.1.e-6)
then
1934 ierrc(i)=
"total normalized condensate too small"
1941 if(pr_ens(i,nens3).lt.1.e-5)
then
1968 call cup_kbcon(ierrc,cap_max_increment,iloop,k22x,kbconx,heo_cup, &
1969 heso_cup,hkbo,ierr2,kbmax,po_cup,cap_max, &
1973 z_cup,entr_rate,heo,imid)
1975 call cup_kbcon(ierrc,cap_max_increment,iloop,k22x,kbconx,heo_cup, &
1976 heso_cup,hkbo,ierr3,kbmax,po_cup,cap_max, &
1980 z_cup,entr_rate,heo,imid)
1990 dq=(qo_cup(i,k+1)-qo_cup(i,k))
1992 mconv(i)=mconv(i)+omeg(i,k)*dq/g
2001 flag_shallow = .false.
2007 del(i,k) = delp(i,k)*0.001
2018 call progsigma_calc(itf,ktf,flag_init,flag_restart,flag_shallow, &
2019 flag_mid,del,tmf,qmicro,dbyo1,zdqca,omega_u,zeta,xlv,dtime, &
2020 forceqv_spechum,kbcon,ktop,cnvflg,betascu,betamcu,betadcu, &
2021 sigmind,sigminm,sigmins,sigmain,sigmaout,sigmab)
2026 ierr,ierr2,ierr3,xf_ens,axx,forcing,
progsigma, &
2027 maxens3,mconv,rand_clos, &
2028 po_cup,ktop,omeg,zdo,zdm,k22,zuo,pr_ens,edto,edtm,kbcon, &
2029 ichoice,omegac,sigmab, &
2032 dicycle,tau_ecmwf,aa1_bl,xf_dicycle,xf_progsigma)
2038 if(ierr(i).eq.0)
then
2039 dellat_ens(i,k,1)=dellat(i,k)
2040 dellaq_ens(i,k,1)=dellaq(i,k)
2041 dellaqc_ens(i,k,1)=dellaqc(i,k)
2042 pwo_ens(i,k,1)=pwo(i,k) +edto(i)*pwdo(i,k)
2044 dellat_ens(i,k,1)=0.
2045 dellaq_ens(i,k,1)=0.
2046 dellaqc_ens(i,k,1)=0.
2057 if(imid.eq.1 .and. ichoice .le.2)
then
2063 if(ierr(i).eq.0)
then
2066 if(k22(i).lt.kpbl(i)+1)
then
2068 blqe=blqe+100.*dhdt(i,k)*(po_cup(i,k)-po_cup(i,k+1))/g
2070 trash=max((hco(i,kbcon(i))-heo_cup(i,kbcon(i))),1.e1)
2071 xff_mid(i,1)=max(0.,blqe/trash)
2072 xff_mid(i,1)=min(0.1,xff_mid(i,1))
2074 xff_mid(i,2)=min(0.1,.03*zws(i))
2075 forcing(i,1)=xff_mid(i,1)
2076 forcing(i,2)=xff_mid(i,2)
2084 outq,outqc,zuo,pre,pwo_ens,xmb,ktop,
progsigma, &
2085 edto,pwdo,
'deep',ierr2,ierr3, &
2086 po_cup,pr_ens,maxens3, &
2087 sig,closure_n,xland1,xmbm_in,xmbs_in, &
2088 ichoice,imid,ipr,itf,ktf, &
2089 its,ite, kts,kte,dx,sigmab, &
2090 dicycle,xf_dicycle,xf_progsigma)
2097 if(ierr(i).eq.0 .and.pre(i).gt.0.)
then
2099 pre(i)=max(pre(i),0.)
2101 outu(i,1)=dellu(i,1)*xmb(i)
2102 outv(i,1)=dellv(i,1)*xmb(i)
2104 outu(i,k)=.25*(dellu(i,k-1)+2.*dellu(i,k)+dellu(i,k+1))*xmb(i)
2105 outv(i,k)=.25*(dellv(i,k-1)+2.*dellv(i,k)+dellv(i,k+1))*xmb(i)
2107 elseif(ierr(i).ne.0 .or. pre(i).eq.0.)
then
2121 if(irainevap.eq.1)
then
2130 if(ierr(i).eq.0)
then
2132 do k = ktop(i), 1, -1
2133 rain = pwo(i,k) + edto(i) * pwdo(i,k)
2135 rntot(i) = rntot(i) + rain * xmb(i)* .001 * dtime
2142 if(ierr(i).eq.0)
then
2143 evef = edt(i) * evfact * sig(i)**2
2144 if(xland(i).gt.0.5 .and. xland(i).lt.1.5) evef = edt(i) * evfactl * sig(i)**2
2146 do k = ktop(i), 1, -1
2147 rain = pwo(i,k) + edto(i) * pwdo(i,k)
2148 rn(i) = rn(i) + rain * xmb(i) * .001 * dtime
2149 if(k.gt.jmin(i))
then
2151 q1=qo(i,k)+(outq(i,k))*dtime
2152 t1=tn(i,k)+(outt(i,k))*dtime
2153 qcond(i) = evef * (q1 - qeso(i,k)) &
2154 & / (1. + el2orc * qeso(i,k) / t1**2)
2155 dp = -100.*(p_cup(i,k+1)-p_cup(i,k))
2156 if(rn(i).gt.0. .and. qcond(i).lt.0.)
then
2157 qevap(i) = -qcond(i) * (1.-exp(-.32*sqrt(dtime*rn(i))))
2158 qevap(i) = min(qevap(i), rn(i)*1000.*g/dp)
2159 delq2(i) = delqev(i) + .001 * qevap(i) * dp / g
2161 if(rn(i).gt.0..and.qcond(i).lt.0..and. &
2162 & delq2(i).gt.rntot(i))
then
2163 qevap(i) = 1000.* g * (rntot(i) - delqev(i)) / dp
2166 if(rn(i).gt.0..and.qevap(i).gt.0.)
then
2167 outq(i,k) = outq(i,k) + qevap(i)/dtime
2168 outt(i,k) = outt(i,k) - elocp * qevap(i)/dtime
2169 rn(i) = max(0.,rn(i) - .001 * qevap(i) * dp / g)
2170 pre(i) = pre(i) - qevap(i) * dp /g/dtime
2171 pre(i)=max(pre(i),0.)
2172 delqev(i) = delqev(i) + .001*dp*qevap(i)/g
2182 rainevap(i)=delqev(i)
2190 if(ierr(i).eq.0)
then
2191 if(aeroevap.gt.1)
then
2193 ccnloss(i)=ccn(i)*pefc(i)*xmb(i)
2194 ccn(i) = ccn(i) - ccnloss(i)*scav_factor
2205 if(ierr(i).eq.0)
then
2209 dp=(po_cup(i,k)-po_cup(i,k+1))*100.
2211 dts= dts -(outu(i,k)*us(i,k)+outv(i,k)*vs(i,k))*dp/g
2213 fpi = fpi +sqrt(outu(i,k)*outu(i,k) + outv(i,k)*outv(i,k))*dp
2217 fp= sqrt((outu(i,k)*outu(i,k)+outv(i,k)*outv(i,k)))/fpi
2218 outt(i,k)=outt(i,k)+fp*dts*g/cp