PageRenderTime 73ms CodeModel.GetById 32ms RepoModel.GetById 0ms app.codeStats 1ms

/wrfv2_fire/phys/module_cu_g3.F

http://github.com/jbeezley/wrf-fire
FORTRAN Legacy | 5313 lines | 3667 code | 294 blank | 1352 comment | 275 complexity | 39bd7e8616606edac6d7d990e6a7ef28 MD5 | raw file
Possible License(s): AGPL-1.0

Large files files are truncated, but you can click here to view the full file

  1. !WRF:MODEL_LAYER:PHYSICS
  2. !
  3. MODULE module_cu_g3
  4. CONTAINS
  5. !-------------------------------------------------------------
  6. SUBROUTINE G3DRV( &
  7. DT,itimestep,DX &
  8. ,rho,RAINCV,PRATEC &
  9. ,U,V,t,W,q,p,pi &
  10. ,dz8w,p8w,XLV,CP,G,r_v &
  11. ,STEPCU,htop,hbot &
  12. ,CU_ACT_FLAG,warm_rain &
  13. ,APR_GR,APR_W,APR_MC,APR_ST,APR_AS &
  14. ,APR_CAPMA,APR_CAPME,APR_CAPMI &
  15. ,MASS_FLUX,XF_ENS,PR_ENS,HT,XLAND,gsw,edt_out &
  16. ,GDC,GDC2 ,kpbl,k22_shallow,kbcon_shallow &
  17. ,ktop_shallow,xmb_shallow &
  18. ,cugd_tten,cugd_qvten ,cugd_qcten &
  19. ,cugd_ttens,cugd_qvtens,cugd_avedx,imomentum &
  20. ,ensdim,maxiens,maxens,maxens2,maxens3,ichoice &
  21. ,ishallow_g3,ids,ide, jds,jde, kds,kde &
  22. ,ims,ime, jms,jme, kms,kme &
  23. ,ips,ipe, jps,jpe, kps,kpe &
  24. ,its,ite, jts,jte, kts,kte &
  25. ,periodic_x,periodic_y &
  26. ,RQVCUTEN,RQCCUTEN,RQICUTEN &
  27. ,RQVFTEN,RTHFTEN,RTHCUTEN &
  28. ,rqvblten,rthblten &
  29. ,F_QV ,F_QC ,F_QR ,F_QI ,F_QS &
  30. #if ( WRF_DFI_RADAR == 1 )
  31. ! Optional CAP suppress option
  32. ,do_capsuppress,cap_suppress_loc &
  33. #endif
  34. )
  35. !-------------------------------------------------------------
  36. IMPLICIT NONE
  37. !-------------------------------------------------------------
  38. INTEGER, INTENT(IN ) :: &
  39. ids,ide, jds,jde, kds,kde, &
  40. ims,ime, jms,jme, kms,kme, &
  41. ips,ipe, jps,jpe, kps,kpe, &
  42. its,ite, jts,jte, kts,kte
  43. LOGICAL periodic_x,periodic_y
  44. integer, parameter :: ens4_spread = 3 ! max(3,cugd_avedx)
  45. integer, parameter :: ens4=ens4_spread*ens4_spread
  46. integer, intent (in ) :: &
  47. ensdim,maxiens,maxens,maxens2,maxens3,ichoice
  48. INTEGER, INTENT(IN ) :: STEPCU, ITIMESTEP,cugd_avedx, &
  49. ishallow_g3,imomentum
  50. LOGICAL, INTENT(IN ) :: warm_rain
  51. REAL, INTENT(IN ) :: XLV, R_v
  52. REAL, INTENT(IN ) :: CP,G
  53. REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , &
  54. INTENT(IN ) :: &
  55. U, &
  56. V, &
  57. W, &
  58. pi, &
  59. t, &
  60. q, &
  61. p, &
  62. dz8w, &
  63. p8w, &
  64. rho
  65. REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , &
  66. OPTIONAL , &
  67. INTENT(INOUT ) :: &
  68. GDC,GDC2
  69. REAL, DIMENSION( ims:ime , jms:jme ),INTENT(IN) :: GSW,HT,XLAND
  70. INTEGER, DIMENSION( ims:ime , jms:jme ),INTENT(IN) :: KPBL
  71. INTEGER, DIMENSION( ims:ime , jms:jme ),INTENT(INOUT) :: k22_shallow, &
  72. kbcon_shallow,ktop_shallow
  73. !
  74. REAL, INTENT(IN ) :: DT, DX
  75. !
  76. REAL, DIMENSION( ims:ime , jms:jme ), &
  77. INTENT(INOUT) :: pratec,RAINCV, MASS_FLUX, &
  78. APR_GR,APR_W,APR_MC,APR_ST,APR_AS, &
  79. edt_out,APR_CAPMA,APR_CAPME,APR_CAPMI, &
  80. htop,hbot,xmb_shallow
  81. !+lxz
  82. ! REAL, DIMENSION( ims:ime , jms:jme ) :: & !, INTENT(INOUT) :: &
  83. ! HTOP, &! highest model layer penetrated by cumulus since last reset in radiation_driver
  84. ! HBOT ! lowest model layer penetrated by cumulus since last reset in radiation_driver
  85. ! ! HBOT>HTOP follow physics leveling convention
  86. LOGICAL, DIMENSION( ims:ime , jms:jme ), &
  87. INTENT(INOUT) :: CU_ACT_FLAG
  88. !
  89. ! Optionals
  90. !
  91. REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
  92. OPTIONAL, &
  93. INTENT(INOUT) :: RTHFTEN, &
  94. cugd_tten,cugd_qvten,cugd_qcten, &
  95. cugd_ttens,cugd_qvtens, &
  96. RQVFTEN
  97. REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
  98. OPTIONAL, &
  99. INTENT(INOUT) :: &
  100. RTHCUTEN, &
  101. RQVCUTEN, &
  102. RQVBLTEN, &
  103. RTHBLTEN, &
  104. RQCCUTEN, &
  105. RQICUTEN
  106. !
  107. ! Flags relating to the optional tendency arrays declared above
  108. ! Models that carry the optional tendencies will provdide the
  109. ! optional arguments at compile time; these flags all the model
  110. ! to determine at run-time whether a particular tracer is in
  111. ! use or not.
  112. !
  113. LOGICAL, OPTIONAL :: &
  114. F_QV &
  115. ,F_QC &
  116. ,F_QR &
  117. ,F_QI &
  118. ,F_QS
  119. #if ( WRF_DFI_RADAR == 1 )
  120. !
  121. ! option of cap suppress:
  122. ! do_capsuppress = 1 do
  123. ! do_capsuppress = other don't
  124. !
  125. !
  126. INTEGER, INTENT(IN ) ,OPTIONAL :: do_capsuppress
  127. REAL, DIMENSION( ims:ime, jms:jme ),INTENT(IN ),OPTIONAL :: cap_suppress_loc
  128. REAL, DIMENSION( its:ite ) :: cap_suppress_j
  129. #endif
  130. ! LOCAL VARS
  131. real, dimension(ims:ime,jms:jme,1:ensdim),intent(inout) :: &
  132. xf_ens,pr_ens
  133. real, dimension ( its:ite , jts:jte , 1:ensdim) :: &
  134. massflni,xfi_ens,pri_ens
  135. REAL, DIMENSION( its:ite , jts:jte ) :: MASSI_FLX, &
  136. APRi_GR,APRi_W,APRi_MC,APRi_ST,APRi_AS, &
  137. edti_out,APRi_CAPMA,APRi_CAPME,APRi_CAPMI,gswi
  138. real, dimension (its:ite,kts:kte) :: &
  139. SUBT,SUBQ,OUTT,OUTQ,OUTQC,phh,subm,cupclw,dhdt, &
  140. outts,outqs
  141. real, dimension (its:ite) :: &
  142. pret, ter11, aa0, fp,xlandi
  143. !+lxz
  144. integer, dimension (its:ite) :: &
  145. kbcon, ktop,kpbli,k22s,kbcons,ktops
  146. !.lxz
  147. integer, dimension (its:ite,jts:jte) :: &
  148. iact_old_gr
  149. integer :: iens,ibeg,iend,jbeg,jend,n,nn,ens4n
  150. integer :: ibegh,iendh,jbegh,jendh
  151. integer :: ibegc,iendc,jbegc,jendc
  152. !
  153. ! basic environmental input includes moisture convergence (mconv)
  154. ! omega (omeg), windspeed (us,vs), and a flag (aaeq) to turn off
  155. ! convection for this call only and at that particular gridpoint
  156. !
  157. real, dimension (its:ite,kts:kte) :: &
  158. T2d,q2d,PO,P2d,US,VS,tn,qo,tshall,qshall
  159. real, dimension (ips-2:ipe+2,kps:kpe,jps-2:jpe+2) :: &
  160. ave_f_t,ave_f_q
  161. real, dimension (its:ite,kts:kte,1:ens4) :: &
  162. omeg,tx,qx
  163. real, dimension (its:ite) :: &
  164. Z1,PSUR,AAEQ,direction,cuten,umean,vmean,pmean,xmbs
  165. real, dimension (its:ite,1:ens4) :: &
  166. mconv
  167. INTEGER :: i,j,k,ICLDCK,ipr,jpr
  168. REAL :: tcrit,tscl_KF,dp,dq,sub_spread,subcenter
  169. INTEGER :: itf,jtf,ktf,iss,jss,nbegin,nend
  170. INTEGER :: high_resolution
  171. REAL :: rkbcon,rktop !-lxz
  172. ! ruc variable
  173. real, dimension (its:ite) :: tkm
  174. ! A. Betts for shallow convection: suggestion for the KF timescale < DELTAX / 25 m/s
  175. tscl_kf=dx/25.
  176. !
  177. ! write(0,*)'ishallow = ',ishallow_g3
  178. high_resolution=0
  179. if(cugd_avedx.gt.1) high_resolution=1
  180. subcenter=0.
  181. ! subcenter=1./float(cugd_avedx)
  182. sub_spread=max(1.,float(cugd_avedx*cugd_avedx-1))
  183. sub_spread=(1.-subcenter)/sub_spread
  184. iens=1
  185. ipr=43
  186. jpr=1
  187. ipr=0
  188. jpr=0
  189. ! if(itimestep.eq.8)then
  190. ! ipr=37
  191. ! jpr=16
  192. ! endif
  193. IF ( periodic_x ) THEN
  194. ibeg=max(its,ids)
  195. iend=min(ite,ide-1)
  196. ibegc=max(its,ids)
  197. iendc=min(ite,ide-1)
  198. ELSE
  199. ibeg=max(its,ids)
  200. iend=min(ite,ide-1)
  201. ibegc=max(its,ids+4)
  202. iendc=min(ite,ide-5)
  203. END IF
  204. IF ( periodic_y ) THEN
  205. jbeg=max(jts,jds)
  206. jend=min(jte,jde-1)
  207. jbegc=max(jts,jds)
  208. jendc=min(jte,jde-1)
  209. ELSE
  210. jbeg=max(jts,jds)
  211. jend=min(jte,jde-1)
  212. jbegc=max(jts,jds+4)
  213. jendc=min(jte,jde-5)
  214. END IF
  215. do j=jts,jte
  216. do i=its,ite
  217. k22_shallow(i,j)=0
  218. kbcon_shallow(i,j)=0
  219. ktop_shallow(i,j)=0
  220. xmb_shallow(i,j)=0
  221. enddo
  222. enddo
  223. tcrit=258.
  224. ave_f_t=0.
  225. ave_f_q=0.
  226. itf=MIN(ite,ide-1)
  227. ktf=MIN(kte,kde-1)
  228. jtf=MIN(jte,jde-1)
  229. !
  230. #if ( EM_CORE == 1 )
  231. if(high_resolution.eq.1)then
  232. !
  233. ! calculate these on the halo...the incominh tendencies have been exchanged on a 24pt halo
  234. ! only neede for high resolution run
  235. !
  236. ibegh=its
  237. jbegh=jts
  238. iendh=ite
  239. jendh=jte
  240. if(its.eq.ips)ibegh=max(its-1,ids)
  241. if(jts.eq.jps)jbegh=max(jts-1,jds)
  242. if(jte.eq.jpe)jendh=min(jte+1,jde-1)
  243. if(ite.eq.ipe)iendh=min(ite+1,ide-1)
  244. DO J = jbegh,jendh
  245. DO k= kts,ktf
  246. DO I= ibegh,iendh
  247. ave_f_t(i,k,j)=(rthften(i-1,k,j-1)+rthften(i-1,k,j) + rthften(i-1,k,j+1)+ &
  248. rthften(i,k,j-1) +rthften(i,k,j) +rthften(i,k,j+1)+ &
  249. rthften(i+1,k,j-1) +rthften(i+1,k,j) +rthften(i+1,k,j+1))/9.
  250. ave_f_q(i,k,j)=(rqvften(i-1,k,j-1)+rqvften(i-1,k,j) + rqvften(i-1,k,j+1)+ &
  251. rqvften(i,k,j-1) +rqvften(i,k,j) +rqvften(i,k,j+1)+ &
  252. rqvften(i+1,k,j-1) +rqvften(i+1,k,j) +rqvften(i+1,k,j+1))/9.
  253. ! ave_f_t(i,k,j)=rthften(i,k,j)
  254. ! ave_f_q(i,k,j)=rqvften(i,k,j)
  255. ENDDO
  256. ENDDO
  257. ENDDO
  258. endif
  259. #endif
  260. DO 100 J = jts,jtf
  261. DO n= 1,ensdim
  262. DO I= its,itf
  263. xfi_ens(i,j,n)=0.
  264. pri_ens(i,j,n)=0.
  265. ! xfi_ens(i,j,n)=xf_ens(i,j,n)
  266. ! pri_ens(i,j,n)=pr_ens(i,j,n)
  267. ENDDO
  268. ENDDO
  269. DO I= its,itf
  270. kbcon(i)=0
  271. ktop(i)=0
  272. tkm(i)=0.
  273. HBOT(I,J) =REAL(KTE)
  274. HTOP(I,J) =REAL(KTS)
  275. iact_old_gr(i,j)=0
  276. mass_flux(i,j)=0.
  277. massi_flx(i,j)=0.
  278. raincv(i,j)=0.
  279. pratec (i,j)=0.
  280. edt_out(i,j)=0.
  281. edti_out(i,j)=0.
  282. gswi(i,j)=gsw(i,j)
  283. xlandi(i)=xland(i,j)
  284. APRi_GR(i,j)=apr_gr(i,j)
  285. APRi_w(i,j)=apr_w(i,j)
  286. APRi_mc(i,j)=apr_mc(i,j)
  287. APRi_st(i,j)=apr_st(i,j)
  288. APRi_as(i,j)=apr_as(i,j)
  289. APRi_capma(i,j)=apr_capma(i,j)
  290. APRi_capme(i,j)=apr_capme(i,j)
  291. APRi_capmi(i,j)=apr_capmi(i,j)
  292. CU_ACT_FLAG(i,j) = .true.
  293. ENDDO
  294. do k=kts,kte
  295. DO I= its,itf
  296. cugd_tten(i,k,j)=0.
  297. cugd_ttens(i,k,j)=0.
  298. cugd_qvten(i,k,j)=0.
  299. cugd_qvtens(i,k,j)=0.
  300. cugd_qcten(i,k,j)=0.
  301. ENDDO
  302. ENDDO
  303. DO n=1,ens4
  304. DO I= its,itf
  305. mconv(i,n)=0.
  306. ENDDO
  307. do k=kts,kte
  308. DO I= its,itf
  309. omeg(i,k,n)=0.
  310. tx(i,k,n)=0.
  311. qx(i,k,n)=0.
  312. ENDDO
  313. ENDDO
  314. ENDDO
  315. DO k=1,ensdim
  316. DO I= its,itf
  317. massflni(i,j,k)=0.
  318. ENDDO
  319. ENDDO
  320. ! put hydrostatic pressure on half levels
  321. DO K=kts,ktf
  322. DO I=ITS,ITF
  323. phh(i,k) = p(i,k,j)
  324. ENDDO
  325. ENDDO
  326. DO I=ITS,ITF
  327. PSUR(I)=p8w(I,1,J)*.01
  328. ! PSUR(I)=p(I,1,J)*.01
  329. TER11(I)=HT(i,j)
  330. aaeq(i)=0.
  331. direction(i)=0.
  332. pret(i)=0.
  333. umean(i)=0.
  334. vmean(i)=0.
  335. pmean(i)=0.
  336. kpbli(i)=kpbl(i,j)
  337. ENDDO
  338. ! if(j.eq.jpr)write(0,*)'psur(ipr),ter11(ipr),kpbli(ipr)'
  339. ! if(j.eq.jpr)write(0,*)psur(ipr),ter11(ipr),kpbli(ipr),r_v
  340. DO K=kts,ktf
  341. DO I=ITS,ITF
  342. po(i,k)=phh(i,k)*.01
  343. subm(i,k)=0.
  344. P2d(I,K)=PO(i,k)
  345. US(I,K) =u(i,k,j)
  346. VS(I,K) =v(i,k,j)
  347. T2d(I,K)=t(i,k,j)
  348. q2d(I,K)=q(i,k,j)
  349. IF(Q2d(I,K).LT.1.E-08)Q2d(I,K)=1.E-08
  350. SUBT(I,K)=0.
  351. SUBQ(I,K)=0.
  352. OUTT(I,K)=0.
  353. OUTQ(I,K)=0.
  354. OUTQC(I,K)=0.
  355. OUTTS(I,K)=0.
  356. OUTQS(I,K)=0.
  357. TN(I,K)=t2d(i,k)+RTHFTEN(i,k,j)*dt
  358. QO(I,K)=q2d(i,k)+RQVFTEN(i,k,j)*dt
  359. TSHALL(I,K)=t2d(i,k)+RTHBLTEN(i,k,j)*pi(i,k,j)*dt
  360. DHDT(I,K)=cp*RTHBLTEN(i,k,j)*pi(i,k,j)+ XLV*RQVBLTEN(i,k,j)
  361. QSHALL(I,K)=q2d(i,k)+RQVBLTEN(i,k,j)*dt
  362. if(high_resolution.eq.1)then
  363. TN(I,K)=t2d(i,k)+ave_f_t(i,k,j)*dt
  364. QO(I,K)=q2d(i,k)+ave_f_q(i,k,j)*dt
  365. endif
  366. IF(TN(I,K).LT.200.)TN(I,K)=T2d(I,K)
  367. IF(QO(I,K).LT.1.E-08)QO(I,K)=1.E-08
  368. ! if(i.eq.ipr.and.j.eq.jpr)then
  369. ! write(0,123)k,p2d(i,k),t2d(i,k),tn(i,k),q2d(i,k),QO(i,k),RTHBLTEN(i,k,j),RQVBLTEN(i,k,j)
  370. ! endif
  371. ENDDO
  372. ENDDO
  373. 123 format(1x,i2,f8.0,1x,2(1x,f8.3),4(1x,e12.4))
  374. ens4n=0
  375. nbegin=0
  376. nend=0
  377. if(ens4_spread.gt.1)then
  378. nbegin=-ens4_spread/2
  379. nend=ens4_spread/2
  380. endif
  381. do nn=nbegin,nend,1
  382. jss=max(j+nn,jds+0)
  383. jss=min(jss,jde-1)
  384. do n=nbegin,nend,1
  385. ens4n=ens4n+1
  386. DO K=kts,ktf
  387. DO I=ITS,ITF
  388. iss=max(i+n,ids+0)
  389. iss=min(iss,ide-1)
  390. omeg(I,K,ens4n)= -g*rho(i,k,j)*w(iss,k,jss)
  391. ! omeg(I,K,ens4n)= -g*rho(i,k,j)*w(i,k,j)
  392. Tx(I,K,ens4n)=t2d(i,k)+RTHFTEN(iss,k,jss)*dt
  393. ! Tx(I,K,ens4n)=t2d(i,k)+RTHFTEN(i,k,j)*dt
  394. if(high_resolution.eq.1)Tx(I,K,ens4n)=t2d(i,k)+ave_f_t(iss,k,jss)*dt
  395. IF(Tx(I,K,ens4n).LT.200.)Tx(I,K,ens4n)=T2d(I,K)
  396. Qx(I,K,ens4n)=q2d(i,k)+RQVFTEN(iss,k,jss)*dt
  397. Qx(I,K,ens4n)=q2d(i,k)+RQVFTEN(i,k,j)*dt
  398. if(high_resolution.eq.1)qx(I,K,ens4n)=q2d(i,k)+ave_f_q(iss,k,jss)*dt
  399. IF(Qx(I,K,ens4n).LT.1.E-08)Qx(I,K,ens4n)=1.E-08
  400. enddo
  401. enddo
  402. enddo !n
  403. enddo !nn
  404. do k= kts+1,ktf-1
  405. DO I = its,itf
  406. if((p2d(i,1)-p2d(i,k)).gt.150.and.p2d(i,k).gt.300)then
  407. dp=-.5*(p2d(i,k+1)-p2d(i,k-1))
  408. umean(i)=umean(i)+us(i,k)*dp
  409. vmean(i)=vmean(i)+vs(i,k)*dp
  410. pmean(i)=pmean(i)+dp
  411. endif
  412. enddo
  413. enddo
  414. DO I = its,itf
  415. umean(i)=umean(i)/pmean(i)
  416. vmean(i)=vmean(i)/pmean(i)
  417. direction(i)=(atan2(umean(i),vmean(i))+3.1415926)*57.29578
  418. if(direction(i).gt.360.)direction(i)=direction(i)-360.
  419. ENDDO
  420. do n=1,ens4
  421. DO K=kts,ktf-1
  422. DO I = its,itf
  423. dq=(q2d(i,k+1)-q2d(i,k))
  424. mconv(i,n)=mconv(i,n)+omeg(i,k,n)*dq/g
  425. enddo
  426. ENDDO
  427. ENDDO
  428. do n=1,ens4
  429. DO I = its,itf
  430. if(mconv(i,n).lt.0.)mconv(i,n)=0.
  431. ENDDO
  432. ENDDO
  433. !
  434. !---- CALL CUMULUS PARAMETERIZATION
  435. !
  436. #if ( WRF_DFI_RADAR == 1 )
  437. if(do_capsuppress == 1 ) then
  438. DO I= its,itf
  439. cap_suppress_j(i)=cap_suppress_loc(i,j)
  440. ENDDO
  441. endif
  442. #endif
  443. CALL CUP_enss_3d(outqc,j,AAEQ,T2d,Q2d,TER11,subm,TN,QO,PO,PRET,&
  444. P2d,OUTT,OUTQ,DT,itimestep,tkm,PSUR,US,VS,tcrit,iens,tx,qx, &
  445. tshall,qshall,kpbli,DHDT,outts,outqs,tscl_kf, &
  446. k22s,kbcons,ktops,xmbs, &
  447. mconv,massflni,iact_old_gr,omeg,direction,MASSi_FLX, &
  448. maxiens,maxens,maxens2,maxens3,ensdim, &
  449. APRi_GR,APRi_W,APRi_MC,APRi_ST,APRi_AS, &
  450. APRi_CAPMA,APRi_CAPME,APRi_CAPMI,kbcon,ktop,cupclw, &
  451. xfi_ens,pri_ens,XLANDi,gswi,edti_out,subt,subq, &
  452. ! ruc lv_p,rv_p,cpd_p,g0_p,ichoice,ipr,jpr, &
  453. xlv,r_v,cp,g,ichoice,ipr,jpr,ens4,high_resolution, &
  454. ishallow_g3,itf,jtf,ktf, &
  455. its,ite, jts,jte, kts,kte &
  456. #if ( WRF_DFI_RADAR == 1 )
  457. ,do_capsuppress,cap_suppress_j &
  458. #endif
  459. )
  460. if(j.lt.jbegc.or.j.gt.jendc)go to 100
  461. DO I=ibegc,iendc
  462. xmb_shallow(i,j)=xmbs(i)
  463. k22_shallow(i,j)=k22s(i)
  464. kbcon_shallow(i,j)=kbcons(i)
  465. ktop_shallow(i,j)=ktops(i)
  466. cuten(i)=0.
  467. if(pret(i).gt.0.)then
  468. cuten(i)=1.
  469. ! raincv(i,j)=pret(i)*dt
  470. endif
  471. ENDDO
  472. ! if(j.eq.jpr)write(0,*)'precip,ktop,kbcon = ',pret(ipr),ktop(ipr),kbcon(ipr)
  473. DO I=ibegc,iendc
  474. DO K=kts,ktf
  475. cugd_ttens(I,K,J)=subt(i,k)*cuten(i)*sub_spread
  476. cugd_qvtens(I,K,J)=subq(i,k)*cuten(i)*sub_spread
  477. ! cugd_tten(I,K,J)=outt(i,k)*cuten(i)
  478. ! cugd_qvten(I,K,J)=outq(i,k)*cuten(i)
  479. cugd_tten(I,K,J)=outts(i,k)+outt(i,k)*cuten(i)
  480. cugd_qvten(I,K,J)=outqs(i,k)+outq(i,k)*cuten(i)
  481. cugd_qcten(I,K,J)=outqc(i,k)*cuten(i)
  482. ! if(i.eq.ipr.and.j.eq.jpr)then
  483. ! write(0,*)subt(i,k)+outt(i,k),subq(i,k)+outq(i,k),outts(i,k),outqs(i,k)
  484. ! endif
  485. ENDDO
  486. ENDDO
  487. DO I=ibegc,iendc
  488. if(pret(i).gt.0.)then
  489. raincv(i,j)=pret(i)*dt
  490. pratec(i,j)=pret(i)
  491. rkbcon = kte+kts - kbcon(i)
  492. rktop = kte+kts - ktop(i)
  493. if (ktop(i) > HTOP(i,j)) HTOP(i,j) = ktop(i)+.001
  494. if (kbcon(i) < HBOT(i,j)) HBOT(i,j) = kbcon(i)+.001
  495. endif
  496. ENDDO
  497. DO n= 1,ensdim
  498. DO I= ibegc,iendc
  499. xf_ens(i,j,n)=xfi_ens(i,j,n)
  500. pr_ens(i,j,n)=pri_ens(i,j,n)
  501. ENDDO
  502. ENDDO
  503. DO I= ibegc,iendc
  504. APR_GR(i,j)=apri_gr(i,j)
  505. APR_w(i,j)=apri_w(i,j)
  506. APR_mc(i,j)=apri_mc(i,j)
  507. APR_st(i,j)=apri_st(i,j)
  508. APR_as(i,j)=apri_as(i,j)
  509. APR_capma(i,j)=apri_capma(i,j)
  510. APR_capme(i,j)=apri_capme(i,j)
  511. APR_capmi(i,j)=apri_capmi(i,j)
  512. mass_flux(i,j)=massi_flx(i,j)
  513. edt_out(i,j)=edti_out(i,j)
  514. ENDDO
  515. IF(PRESENT(RQCCUTEN)) THEN
  516. IF ( F_QC ) THEN
  517. DO K=kts,ktf
  518. DO I=ibegc,iendc
  519. RQCCUTEN(I,K,J)=outqc(I,K)*cuten(i)
  520. IF ( PRESENT( GDC ) ) GDC(I,K,J)=CUPCLW(I,K)*cuten(i)
  521. IF ( PRESENT( GDC2 ) ) GDC2(I,K,J)=0.
  522. ENDDO
  523. ENDDO
  524. ENDIF
  525. ENDIF
  526. !...... QSTEN STORES GRAUPEL TENDENCY IF IT EXISTS, OTHERISE SNOW (V2)
  527. IF(PRESENT(RQICUTEN).AND.PRESENT(RQCCUTEN))THEN
  528. IF (F_QI) THEN
  529. DO K=kts,ktf
  530. DO I=ibegc,iendc
  531. if(t2d(i,k).lt.258.)then
  532. RQICUTEN(I,K,J)=outqc(I,K)*cuten(i)
  533. cugd_qcten(i,k,j)=0.
  534. RQCCUTEN(I,K,J)=0.
  535. IF ( PRESENT( GDC2 ) ) GDC2(I,K,J)=CUPCLW(I,K)*cuten(i)
  536. else
  537. RQICUTEN(I,K,J)=0.
  538. RQCCUTEN(I,K,J)=outqc(I,K)*cuten(i)
  539. IF ( PRESENT( GDC ) ) GDC(I,K,J)=CUPCLW(I,K)*cuten(i)
  540. endif
  541. ENDDO
  542. ENDDO
  543. ENDIF
  544. ENDIF
  545. 100 continue
  546. END SUBROUTINE G3DRV
  547. SUBROUTINE CUP_enss_3d(OUTQC,J,AAEQ,T,Q,Z1,sub_mas, &
  548. TN,QO,PO,PRE,P,OUTT,OUTQ,DTIME,ktau,tkmax,PSUR,US,VS, &
  549. TCRIT,iens,tx,qx, &
  550. tshall,qshall,kpbl,dhdt,outts,outqs,tscl_kf, &
  551. k23,kbcon3,ktop3,xmb3, &
  552. mconv,massfln,iact, &
  553. omeg,direction,massflx,maxiens, &
  554. maxens,maxens2,maxens3,ensdim, &
  555. APR_GR,APR_W,APR_MC,APR_ST,APR_AS, &
  556. APR_CAPMA,APR_CAPME,APR_CAPMI,kbcon,ktop,cupclw, & !-lxz
  557. xf_ens,pr_ens,xland,gsw,edt_out,subt,subq, &
  558. xl,rv,cp,g,ichoice,ipr,jpr,ens4,high_resolution, &
  559. ishallow_g3,itf,jtf,ktf, &
  560. its,ite, jts,jte, kts,kte &
  561. #if ( WRF_DFI_RADAR == 1 )
  562. ! Optional CAP suppress option
  563. ,do_capsuppress,cap_suppress_j &
  564. #endif
  565. )
  566. IMPLICIT NONE
  567. integer &
  568. ,intent (in ) :: &
  569. itf,jtf,ktf,ktau, &
  570. its,ite, jts,jte, kts,kte,ipr,jpr,ens4,high_resolution
  571. integer, intent (in ) :: &
  572. j,ensdim,maxiens,ishallow_g3,maxens,maxens2,maxens3,ichoice,iens
  573. !
  574. !
  575. !
  576. real, dimension (its:ite,jts:jte,1:ensdim) &
  577. ,intent (inout) :: &
  578. massfln,xf_ens,pr_ens
  579. real, dimension (its:ite,jts:jte) &
  580. ,intent (inout ) :: &
  581. APR_GR,APR_W,APR_MC,APR_ST,APR_AS,APR_CAPMA, &
  582. APR_CAPME,APR_CAPMI,massflx,edt_out
  583. real, dimension (its:ite,jts:jte) &
  584. ,intent (in ) :: &
  585. gsw
  586. integer, dimension (its:ite,jts:jte) &
  587. ,intent (in ) :: &
  588. iact
  589. ! outtem = output temp tendency (per s)
  590. ! outq = output q tendency (per s)
  591. ! outqc = output qc tendency (per s)
  592. ! pre = output precip
  593. real, dimension (its:ite,kts:kte) &
  594. ,intent (inout ) :: &
  595. DHDT,OUTT,OUTQ,OUTQC,subt,subq,sub_mas,cupclw,outts,outqs
  596. real, dimension (its:ite) &
  597. ,intent (out ) :: &
  598. pre,xmb3
  599. integer, dimension (its:ite) &
  600. ,intent (out ) :: &
  601. kbcon,ktop,k23,kbcon3,ktop3
  602. integer, dimension (its:ite) &
  603. ,intent (in ) :: &
  604. kpbl
  605. !
  606. ! basic environmental input includes moisture convergence (mconv)
  607. ! omega (omeg), windspeed (us,vs), and a flag (aaeq) to turn off
  608. ! convection for this call only and at that particular gridpoint
  609. !
  610. real, dimension (its:ite,kts:kte) &
  611. ,intent (in ) :: &
  612. T,PO,P,US,VS,tn,tshall,qshall
  613. real, dimension (its:ite,kts:kte,1:ens4) &
  614. ,intent (inout ) :: &
  615. omeg,tx,qx
  616. real, dimension (its:ite,kts:kte) &
  617. ,intent (inout) :: &
  618. Q,QO
  619. real, dimension (its:ite) &
  620. ,intent (in ) :: &
  621. Z1,PSUR,AAEQ,direction,tkmax,xland
  622. real, dimension (its:ite,1:ens4) &
  623. ,intent (in ) :: &
  624. mconv
  625. real &
  626. ,intent (in ) :: &
  627. dtime,tcrit,xl,cp,rv,g,tscl_kf
  628. #if ( WRF_DFI_RADAR == 1 )
  629. !
  630. ! option of cap suppress:
  631. ! do_capsuppress = 1 do
  632. ! do_capsuppress = other don't
  633. !
  634. !
  635. INTEGER, INTENT(IN ) ,OPTIONAL :: do_capsuppress
  636. REAL, DIMENSION( its:ite ),INTENT(IN ) ,OPTIONAL :: cap_suppress_j
  637. #endif
  638. !
  639. ! local ensemble dependent variables in this routine
  640. !
  641. real, dimension (its:ite,1:maxens) :: &
  642. xaa0_ens
  643. real, dimension (1:maxens) :: &
  644. mbdt_ens
  645. real, dimension (1:maxens2) :: &
  646. edt_ens
  647. real, dimension (its:ite,1:maxens2) :: &
  648. edtc
  649. real, dimension (its:ite,kts:kte,1:maxens2) :: &
  650. dellat_ens,dellaqc_ens,dellaq_ens,pwo_ens,subt_ens,subq_ens
  651. !
  652. !
  653. !
  654. !***************** the following are your basic environmental
  655. ! variables. They carry a "_cup" if they are
  656. ! on model cloud levels (staggered). They carry
  657. ! an "o"-ending (z becomes zo), if they are the forced
  658. ! variables. They are preceded by x (z becomes xz)
  659. ! to indicate modification by some typ of cloud
  660. !
  661. ! z = heights of model levels
  662. ! q = environmental mixing ratio
  663. ! qes = environmental saturation mixing ratio
  664. ! t = environmental temp
  665. ! p = environmental pressure
  666. ! he = environmental moist static energy
  667. ! hes = environmental saturation moist static energy
  668. ! z_cup = heights of model cloud levels
  669. ! q_cup = environmental q on model cloud levels
  670. ! qes_cup = saturation q on model cloud levels
  671. ! t_cup = temperature (Kelvin) on model cloud levels
  672. ! p_cup = environmental pressure
  673. ! he_cup = moist static energy on model cloud levels
  674. ! hes_cup = saturation moist static energy on model cloud levels
  675. ! gamma_cup = gamma on model cloud levels
  676. !
  677. !
  678. ! hcd = moist static energy in downdraft
  679. ! zd normalized downdraft mass flux
  680. ! dby = buoancy term
  681. ! entr = entrainment rate
  682. ! zd = downdraft normalized mass flux
  683. ! entr= entrainment rate
  684. ! hcd = h in model cloud
  685. ! bu = buoancy term
  686. ! zd = normalized downdraft mass flux
  687. ! gamma_cup = gamma on model cloud levels
  688. ! mentr_rate = entrainment rate
  689. ! qcd = cloud q (including liquid water) after entrainment
  690. ! qrch = saturation q in cloud
  691. ! pwd = evaporate at that level
  692. ! pwev = total normalized integrated evaoprate (I2)
  693. ! entr= entrainment rate
  694. ! z1 = terrain elevation
  695. ! entr = downdraft entrainment rate
  696. ! jmin = downdraft originating level
  697. ! kdet = level above ground where downdraft start detraining
  698. ! psur = surface pressure
  699. ! z1 = terrain elevation
  700. ! pr_ens = precipitation ensemble
  701. ! xf_ens = mass flux ensembles
  702. ! massfln = downdraft mass flux ensembles used in next timestep
  703. ! omeg = omega from large scale model
  704. ! mconv = moisture convergence from large scale model
  705. ! zd = downdraft normalized mass flux
  706. ! zu = updraft normalized mass flux
  707. ! dir = "storm motion"
  708. ! mbdt = arbitrary numerical parameter
  709. ! dtime = dt over which forcing is applied
  710. ! iact_gr_old = flag to tell where convection was active
  711. ! kbcon = LFC of parcel from k22
  712. ! k22 = updraft originating level
  713. ! icoic = flag if only want one closure (usually set to zero!)
  714. ! dby = buoancy term
  715. ! ktop = cloud top (output)
  716. ! xmb = total base mass flux
  717. ! hc = cloud moist static energy
  718. ! hkb = moist static energy at originating level
  719. ! mentr_rate = entrainment rate
  720. real, dimension (its:ite,kts:kte) :: &
  721. he3,hes3,qes3,z3,zdo3,zu3_0,hc3_0,dby3_0,gamma3_0_cup, &
  722. qes3_cup,q3_cup,he3_cup,hes3_cup,z3_cup,gamma3_cup,t3_cup, &
  723. xhe3,xhes3,xqes3,xz3,xt3,xq3, &
  724. xqes3_cup,xq3_cup,xhe3_cup,xhes3_cup,xz3_cup,xgamma3_cup, &
  725. xt3_cup, &
  726. xdby3,xqc3,xhc3,xqrc3,xzu3, &
  727. dby3,qc3,pw3,hc3,qrc3,zu3,cd3,DELLAH3,DELLAQ3, &
  728. dsubt3,dsubq3,DELLAT3,DELLAQC3
  729. real, dimension (its:ite,kts:kte) :: &
  730. he,hes,qes,z, &
  731. heo,heso,qeso,zo, &
  732. xhe,xhes,xqes,xz,xt,xq, &
  733. qes_cup,q_cup,he_cup,hes_cup,z_cup,p_cup,gamma_cup,t_cup, &
  734. qeso_cup,qo_cup,heo_cup,heso_cup,zo_cup,po_cup,gammao_cup, &
  735. tn_cup, &
  736. xqes_cup,xq_cup,xhe_cup,xhes_cup,xz_cup,xp_cup,xgamma_cup, &
  737. xt_cup, &
  738. dby,qc,qrcd,pwd,pw,hcd,qcd,dbyd,hc,qrc,zu,zd,clw_all, &
  739. dbyo,qco,qrcdo,pwdo,pwo,hcdo,qcdo,dbydo,hco,qrco,zuo,zdo, &
  740. xdby,xqc,xqrcd,xpwd,xpw,xhcd,xqcd,xhc,xqrc,xzu,xzd, &
  741. ! cd = detrainment function for updraft
  742. ! cdd = detrainment function for downdraft
  743. ! dellat = change of temperature per unit mass flux of cloud ensemble
  744. ! dellaq = change of q per unit mass flux of cloud ensemble
  745. ! dellaqc = change of qc per unit mass flux of cloud ensemble
  746. cd,cdd,scr1,DELLAH,DELLAQ,DELLAT,DELLAQC,dsubt,dsubq
  747. ! aa0 cloud work function for downdraft
  748. ! edt = epsilon
  749. ! aa0 = cloud work function without forcing effects
  750. ! aa1 = cloud work function with forcing effects
  751. ! xaa0 = cloud work function with cloud effects (ensemble dependent)
  752. ! edt = epsilon
  753. real, dimension (its:ite) :: &
  754. aa3_0,aa3,hkb3,qkb3,pwav3,bu3,xaa3,xhkb3, &
  755. hkb3_0,edt,edto,edtx,AA1,AA0,XAA0,HKB, &
  756. HKBO,aad,XHKB,QKB,QKBO,edt3, &
  757. XMB,XPWAV,XPWEV,PWAV,PWEV,PWAVO, &
  758. PWEVO,BU,BUO,cap_max,xland1, &
  759. cap_max_increment,closure_n,cap_max3
  760. real, dimension (its:ite,1:ens4) :: &
  761. axx
  762. integer, dimension (its:ite) :: &
  763. kzdown,KDET,K22,KB,JMIN,kstabi,kstabm,K22x,jmin3,kdet3, & !-lxz
  764. KBCONx,KBx,KTOPx,ierr,ierr2,ierr3,KBMAX,ierr5,ierr5_0
  765. integer :: &
  766. nall,iedt,nens,nens3,ki,I,K,KK,iresult
  767. real :: &
  768. day,dz,mbdt,mbdt_s,entr_rate,radius,entrd_rate,mentr_rate,mentrd_rate, &
  769. zcutdown,edtmax,edtmin,depth_min,zkbmax,z_detr,zktop, &
  770. massfld,dh,cap_maxs,trash,entr_rate3,mentr_rate3
  771. integer :: jmini
  772. logical :: keep_going
  773. real xff_shal(9),blqe,xkshal
  774. day=86400.
  775. do i=its,itf
  776. xmb3(i)=0.
  777. closure_n(i)=16.
  778. xland1(i)=1.
  779. if(xland(i).gt.1.5)xland1(i)=0.
  780. ! cap_max_increment(i)=50.
  781. cap_max_increment(i)=25.
  782. enddo
  783. !
  784. !--- specify entrainmentrate and detrainmentrate
  785. !
  786. if(iens.le.4)then
  787. radius=14000.-float(iens)*2000.
  788. else
  789. radius=12000.
  790. endif
  791. !
  792. !--- gross entrainment rate (these may be changed later on in the
  793. !--- program, depending what your detrainment is!!)
  794. !
  795. entr_rate =.2/radius
  796. entr_rate3=.2/200.
  797. !
  798. !--- entrainment of mass
  799. !
  800. mentrd_rate=0.
  801. mentr_rate=entr_rate
  802. mentr_rate3=entr_rate3
  803. !
  804. !--- initial detrainmentrates
  805. !
  806. do k=kts,ktf
  807. do i=its,itf
  808. cupclw(i,k)=0.
  809. cd(i,k)=0.01*entr_rate
  810. cd3(i,k)=entr_rate3
  811. cdd(i,k)=0.
  812. zdo3(i,k)=0.
  813. hcdo(i,k)=0.
  814. qrcdo(i,k)=0.
  815. dellaqc(i,k)=0.
  816. enddo
  817. enddo
  818. !
  819. !--- max/min allowed value for epsilon (ratio downdraft base mass flux/updraft
  820. ! base mass flux
  821. !
  822. edtmax=1.
  823. edtmin=.2
  824. !
  825. !--- minimum depth (m), clouds must have
  826. !
  827. depth_min=500.
  828. !
  829. !--- maximum depth (mb) of capping
  830. !--- inversion (larger cap = no convection)
  831. !
  832. ! cap_maxs=125.
  833. cap_maxs=75.
  834. DO i=its,itf
  835. kbmax(i)=1
  836. jmin3(i)=0
  837. kdet3(i)=0
  838. aa0(i)=0.
  839. aa3_0(i)=0.
  840. aa1(i)=0.
  841. aa3(i)=0.
  842. aad(i)=0.
  843. edt(i)=0.
  844. edt3(i)=0.
  845. kstabm(i)=ktf-1
  846. IERR(i)=0
  847. IERR2(i)=0
  848. IERR3(i)=0
  849. IERR5(i)=0
  850. IERR5_0(i)=0
  851. enddo
  852. !
  853. !--- first check for upstream convection
  854. !
  855. #if ( WRF_DFI_RADAR == 1 )
  856. if(do_capsuppress == 1) then
  857. do i=its,itf
  858. cap_max(i)=cap_maxs
  859. cap_max3(i)=25.
  860. if(gsw(i,j).lt.1.or.high_resolution.eq.1)cap_max(i)=25.
  861. if (abs(cap_suppress_j(i) - 1.0 ) < 0.1 ) then
  862. cap_max(i)=cap_maxs+75.
  863. elseif (abs(cap_suppress_j(i) - 0.0 ) < 0.1 ) then
  864. cap_max(i)=10.0
  865. endif
  866. iresult=0
  867. enddo
  868. else
  869. do i=its,itf
  870. cap_max(i)=cap_maxs
  871. cap_max3(i)=25.
  872. if(gsw(i,j).lt.1.or.high_resolution.eq.1)cap_max(i)=25.
  873. iresult=0
  874. enddo
  875. endif
  876. do i=its,itf
  877. edt_out(i,j)=cap_max(i)
  878. enddo
  879. #else
  880. do i=its,itf
  881. cap_max(i)=cap_maxs
  882. cap_max3(i)=25.
  883. if(gsw(i,j).lt.1.or.high_resolution.eq.1)cap_max(i)=25.
  884. iresult=0
  885. enddo
  886. #endif
  887. !
  888. !--- max height(m) above ground where updraft air can originate
  889. !
  890. zkbmax=4000.
  891. !
  892. !--- height(m) above which no downdrafts are allowed to originate
  893. !
  894. zcutdown=3000.
  895. !
  896. !--- depth(m) over which downdraft detrains all its mass
  897. !
  898. z_detr=1250.
  899. !
  900. do nens=1,maxens
  901. mbdt_ens(nens)=(float(nens)-3.)*dtime*1.e-3+dtime*5.E-03
  902. enddo
  903. do nens=1,maxens2
  904. edt_ens(nens)=.95-float(nens)*.01
  905. enddo
  906. !
  907. !--- environmental conditions, FIRST HEIGHTS
  908. !
  909. do i=its,itf
  910. if(ierr(i).ne.20)then
  911. do k=1,maxens*maxens2*maxens3
  912. xf_ens(i,j,(iens-1)*maxens*maxens2*maxens3+k)=0.
  913. pr_ens(i,j,(iens-1)*maxens*maxens2*maxens3+k)=0.
  914. enddo
  915. endif
  916. enddo
  917. !
  918. !--- calculate moist static energy, heights, qes
  919. !
  920. call cup_env(z,qes,he,hes,t,q,p,z1, &
  921. psur,ierr,tcrit,0,xl,cp, &
  922. itf,jtf,ktf, &
  923. its,ite, jts,jte, kts,kte)
  924. call cup_env(zo,qeso,heo,heso,tn,qo,po,z1, &
  925. psur,ierr,tcrit,0,xl,cp, &
  926. itf,jtf,ktf, &
  927. its,ite, jts,jte, kts,kte)
  928. !
  929. !--- environmental values on cloud levels
  930. !
  931. call cup_env_clev(t,qes,q,he,hes,z,p,qes_cup,q_cup,he_cup, &
  932. hes_cup,z_cup,p_cup,gamma_cup,t_cup,psur, &
  933. ierr,z1,xl,rv,cp, &
  934. itf,jtf,ktf, &
  935. its,ite, jts,jte, kts,kte)
  936. call cup_env_clev(tn,qeso,qo,heo,heso,zo,po,qeso_cup,qo_cup, &
  937. heo_cup,heso_cup,zo_cup,po_cup,gammao_cup,tn_cup,psur, &
  938. ierr,z1,xl,rv,cp, &
  939. itf,jtf,ktf, &
  940. its,ite, jts,jte, kts,kte)
  941. do i=its,itf
  942. if(aaeq(i).lt.-0.1)then
  943. ierr(i)=20
  944. endif
  945. ! if(ierr(i).eq.0)then
  946. !
  947. do k=kts,ktf
  948. if(zo_cup(i,k).gt.zkbmax+z1(i))then
  949. kbmax(i)=k
  950. go to 25
  951. endif
  952. enddo
  953. 25 continue
  954. !
  955. !--- level where detrainment for downdraft starts
  956. !
  957. do k=kts,ktf
  958. if(zo_cup(i,k).gt.z_detr+z1(i))then
  959. kdet(i)=k
  960. go to 26
  961. endif
  962. enddo
  963. 26 continue
  964. !
  965. ! endif
  966. enddo
  967. !
  968. !
  969. !
  970. !------- DETERMINE LEVEL WITH HIGHEST MOIST STATIC ENERGY CONTENT - K22
  971. !
  972. CALL cup_MAXIMI(HEO_CUP,3,KBMAX,K22,ierr, &
  973. itf,jtf,ktf, &
  974. its,ite, jts,jte, kts,kte)
  975. DO 36 i=its,itf
  976. IF(ierr(I).eq.0.)THEN
  977. IF(K22(I).GE.KBMAX(i))ierr(i)=2
  978. endif
  979. 36 CONTINUE
  980. !
  981. !--- DETERMINE THE LEVEL OF CONVECTIVE CLOUD BASE - KBCON
  982. !
  983. call cup_kbcon(cap_max_increment,1,k22,kbcon,heo_cup,heso_cup, &
  984. ierr,kbmax,po_cup,cap_max, &
  985. itf,jtf,ktf, &
  986. its,ite, jts,jte, kts,kte)
  987. !
  988. !--- increase detrainment in stable layers
  989. !
  990. CALL cup_minimi(HEso_cup,Kbcon,kstabm,kstabi,ierr, &
  991. itf,jtf,ktf, &
  992. its,ite, jts,jte, kts,kte)
  993. do i=its,itf
  994. IF(ierr(I).eq.0.)THEN
  995. if(kstabm(i)-1.gt.kstabi(i))then
  996. do k=kstabi(i),kstabm(i)-1
  997. cd(i,k)=cd(i,k-1)+.15*entr_rate
  998. if(cd(i,k).gt.1.0*entr_rate)cd(i,k)=1.0*entr_rate
  999. enddo
  1000. ENDIF
  1001. ENDIF
  1002. ENDDO
  1003. !
  1004. !--- calculate incloud moist static energy
  1005. !
  1006. call cup_up_he(k22,hkb,z_cup,cd,mentr_rate,he_cup,hc, &
  1007. kbcon,ierr,dby,he,hes_cup,'deep', &
  1008. itf,jtf,ktf, &
  1009. its,ite, jts,jte, kts,kte)
  1010. call cup_up_he(k22,hkbo,zo_cup,cd,mentr_rate,heo_cup,hco, &
  1011. kbcon,ierr,dbyo,heo,heso_cup,'deep', &
  1012. itf,jtf,ktf, &
  1013. its,ite, jts,jte, kts,kte)
  1014. !--- DETERMINE CLOUD TOP - KTOP
  1015. !
  1016. call cup_ktop(1,dbyo,kbcon,ktop,ierr, &
  1017. itf,jtf,ktf, &
  1018. its,ite, jts,jte, kts,kte)
  1019. DO 37 i=its,itf
  1020. kzdown(i)=0
  1021. if(ierr(i).eq.0)then
  1022. zktop=(zo_cup(i,ktop(i))-z1(i))*.6
  1023. zktop=min(zktop+z1(i),zcutdown+z1(i))
  1024. do k=kts,kte
  1025. if(zo_cup(i,k).gt.zktop)then
  1026. kzdown(i)=k
  1027. go to 37
  1028. endif
  1029. enddo
  1030. endif
  1031. 37 CONTINUE
  1032. !
  1033. !--- DOWNDRAFT ORIGINATING LEVEL - JMIN
  1034. !
  1035. call cup_minimi(HEso_cup,K22,kzdown,JMIN,ierr, &
  1036. itf,jtf,ktf, &
  1037. its,ite, jts,jte, kts,kte)
  1038. DO 100 i=its,ite
  1039. IF(ierr(I).eq.0.)THEN
  1040. !
  1041. !--- check whether it would have buoyancy, if there where
  1042. !--- no entrainment/detrainment
  1043. !
  1044. jmini = jmin(i)
  1045. keep_going = .TRUE.
  1046. do while ( keep_going )
  1047. keep_going = .FALSE.
  1048. if ( jmini - 1 .lt. kdet(i) ) kdet(i) = jmini-1
  1049. if ( jmini .ge. ktop(i)-1 ) jmini = ktop(i) - 2
  1050. ki = jmini
  1051. hcdo(i,ki)=heso_cup(i,ki)
  1052. DZ=Zo_cup(i,Ki+1)-Zo_cup(i,Ki)
  1053. dh=0.
  1054. do k=ki-1,1,-1
  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))
  1058. if(dh.gt.0.)then
  1059. jmini=jmini-1
  1060. if ( jmini .gt. 3 ) then
  1061. keep_going = .TRUE.
  1062. else
  1063. ierr(i) = 9
  1064. exit
  1065. endif
  1066. endif
  1067. enddo
  1068. enddo
  1069. jmin(i) = jmini
  1070. if ( jmini .le. 3 ) then
  1071. ierr(i)=4
  1072. endif
  1073. ENDIF
  1074. 100 continue
  1075. !
  1076. ! - Must have at least depth_min m between cloud convective base
  1077. ! and cloud top.
  1078. !
  1079. do i=its,itf
  1080. IF(ierr(I).eq.0.)THEN
  1081. IF(-zo_cup(I,KBCON(I))+zo_cup(I,KTOP(I)).LT.depth_min)then
  1082. ierr(i)=6
  1083. endif
  1084. endif
  1085. enddo
  1086. !
  1087. !c--- normalized updraft mass flux profile
  1088. !
  1089. call cup_up_nms(zu,z_cup,mentr_rate,cd,kbcon,ktop,ierr,k22, &
  1090. itf,jtf,ktf, &
  1091. its,ite, jts,jte, kts,kte)
  1092. call cup_up_nms(zuo,zo_cup,mentr_rate,cd,kbcon,ktop,ierr,k22, &
  1093. itf,jtf,ktf, &
  1094. its,ite, jts,jte, kts,kte)
  1095. !
  1096. !c--- normalized downdraft mass flux profile,also work on bottom detrainment
  1097. !--- in this routine
  1098. !
  1099. call cup_dd_nms(zd,z_cup,cdd,mentrd_rate,jmin,ierr, &
  1100. 0,kdet,z1, &
  1101. itf,jtf,ktf, &
  1102. its,ite, jts,jte, kts,kte)
  1103. call cup_dd_nms(zdo,zo_cup,cdd,mentrd_rate,jmin,ierr, &
  1104. 1,kdet,z1, &
  1105. itf,jtf,ktf, &
  1106. its,ite, jts,jte, kts,kte)
  1107. !
  1108. !--- downdraft moist static energy
  1109. !
  1110. call cup_dd_he(hes_cup,zd,hcd,z_cup,cdd,mentrd_rate, &
  1111. jmin,ierr,he,dbyd,he_cup, &
  1112. itf,jtf,ktf, &
  1113. its,ite, jts,jte, kts,kte)
  1114. call cup_dd_he(heso_cup,zdo,hcdo,zo_cup,cdd,mentrd_rate, &
  1115. jmin,ierr,heo,dbydo,he_cup,&
  1116. itf,jtf,ktf, &
  1117. its,ite, jts,jte, kts,kte)
  1118. !
  1119. !--- calculate moisture properties of downdraft
  1120. !
  1121. call cup_dd_moisture_3d(zd,hcd,hes_cup,qcd,qes_cup, &
  1122. pwd,q_cup,z_cup,cdd,mentrd_rate,jmin,ierr,gamma_cup, &
  1123. pwev,bu,qrcd,q,he,t_cup,2,xl,high_resolution, &
  1124. itf,jtf,ktf, &
  1125. its,ite, jts,jte, kts,kte)
  1126. call cup_dd_moisture_3d(zdo,hcdo,heso_cup,qcdo,qeso_cup, &
  1127. pwdo,qo_cup,zo_cup,cdd,mentrd_rate,jmin,ierr,gammao_cup, &
  1128. pwevo,bu,qrcdo,qo,heo,tn_cup,1,xl,high_resolution, &
  1129. itf,jtf,ktf, &
  1130. its,ite, jts,jte, kts,kte)
  1131. !
  1132. !--- calculate moisture properties of updraft
  1133. !
  1134. call cup_up_moisture('deep',ierr,z_cup,qc,qrc,pw,pwav, &
  1135. kbcon,ktop,cd,dby,mentr_rate,clw_all, &
  1136. q,GAMMA_cup,zu,qes_cup,k22,q_cup,xl, &
  1137. itf,jtf,ktf, &
  1138. its,ite, jts,jte, kts,kte)
  1139. do k=kts,ktf
  1140. do i=its,itf
  1141. cupclw(i,k)=qrc(i,k)
  1142. enddo
  1143. enddo
  1144. call cup_up_moisture('deep',ierr,zo_cup,qco,qrco,pwo,pwavo, &
  1145. kbcon,ktop,cd,dbyo,mentr_rate,clw_all, &
  1146. qo,GAMMAo_cup,zuo,qeso_cup,k22,qo_cup,xl,&
  1147. itf,jtf,ktf, &
  1148. its,ite, jts,jte, kts,kte)
  1149. !
  1150. !--- calculate workfunctions for updrafts
  1151. !
  1152. call cup_up_aa0(aa0,z,zu,dby,GAMMA_CUP,t_cup, &
  1153. kbcon,ktop,ierr, &
  1154. itf,jtf,ktf, &
  1155. its,ite, jts,jte, kts,kte)
  1156. call cup_up_aa0(aa1,zo,zuo,dbyo,GAMMAo_CUP,tn_cup, &
  1157. kbcon,ktop,ierr, &
  1158. itf,jtf,ktf, &
  1159. its,ite, jts,jte, kts,kte)
  1160. do i=its,itf
  1161. if(ierr(i).eq.0)then
  1162. if(aa1(i).eq.0.)then
  1163. ierr(i)=17
  1164. endif
  1165. endif
  1166. enddo
  1167. !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  1168. ! NEXT section for shallow convection
  1169. !
  1170. if(ishallow_g3.eq.1)then
  1171. ! write(0,*)'now do shallow for j = ',j
  1172. call cup_env(z3,qes3,he3,hes3,tshall,qshall,po,z1, &
  1173. psur,ierr5,tcrit,0,xl,cp, &
  1174. itf,jtf,ktf, &
  1175. its,ite, jts,jte, kts,kte)
  1176. call cup_env_clev(tshall,qes3,qshall,he3,hes3,z3,po,qes3_cup,q3_cup, &
  1177. he3_cup,hes3_cup,z3_cup,po_cup,gamma3_cup,t3_cup,psur, &
  1178. ierr5,z1,xl,rv,cp, &
  1179. itf,jtf,ktf, &
  1180. its,ite, jts,jte, kts,kte)
  1181. CALL cup_MAXIMI(HE3_CUP,1,kbmax,K23,ierr5, &
  1182. itf,jtf,ktf, &
  1183. its,ite, jts,jte, kts,kte)
  1184. DO i=its,itf
  1185. if(kpbl(i).gt.5)cap_max3(i)=po_cup(i,kpbl(i))
  1186. IF(ierr5(I).eq.0.)THEN
  1187. IF(K23(I).Gt.Kbmax(i))ierr5(i)=2
  1188. if(kpbl(i).gt.5)k23(i)=kpbl(i)
  1189. endif
  1190. ierr5_0(i)=ierr5(i)
  1191. ENDDO
  1192. call cup_kbcon(cap_max_increment,5,k23,kbcon3,he3_cup,hes3_cup, &
  1193. ierr5,kbmax,po_cup,cap_max3, &
  1194. ! ierr5,kpbl,po_cup,cap_max3, &
  1195. itf,jtf,ktf, &
  1196. its,ite, jts,jte, kts,kte)
  1197. call cup_up_he(k23,hkb3,z3_cup,cd3,mentr_rate3,he3_cup,hc3, &
  1198. kbcon3,ierr5,dby3,he3,hes3_cup,'shallow', &
  1199. itf,jtf,ktf, &
  1200. its,ite, jts,jte, kts,kte)
  1201. call cup_up_he(k23,hkb3_0,z_cup,cd3,mentr_rate3,he_cup,hc3_0, &
  1202. kbcon3,ierr5,dby3_0,he,hes_cup,'shallow', &
  1203. itf,jtf,ktf, &
  1204. its,ite, jts,jte, kts,kte)
  1205. call cup_ktop(1,dby3,kbcon3,ktop3,ierr5, &
  1206. itf,jtf,ktf, &
  1207. its,ite, jts,jte, kts,kte)
  1208. call cup_up_nms(zu3,z3_cup,mentr_rate3,cd3,kbcon3,ktop3, &
  1209. ierr5,k23, &
  1210. itf,jtf,ktf, &
  1211. its,ite, jts,jte, kts,kte)
  1212. call cup_up_nms(zu3_0,z_cup,mentr_rate3,cd3,kbcon3,ktop3, &
  1213. ierr5,k23, &
  1214. itf,jtf,ktf, &
  1215. its,ite, jts,jte, kts,kte)
  1216. !
  1217. ! first calculate aa3_0_cup
  1218. !
  1219. call cup_up_aa0(aa3_0,z,zu3_0,dby3_0,GAMMA3_0_CUP,t_cup, &
  1220. kbcon3,ktop3,ierr5, &
  1221. itf,jtf,ktf, &
  1222. its,ite, jts,jte, kts,kte)
  1223. !
  1224. ! now what is necessary for aa3 and feedbacks
  1225. !
  1226. call cup_up_moisture('shallow',ierr5,z3_cup,qc3,qrc3,pw3,pwav3, &
  1227. kbcon3,ktop3,cd3,dby3,mentr_rate3,clw_all, &
  1228. qshall,GAMMA3_cup,zu3,qes3_cup,k23,q3_cup,xl,&
  1229. itf,jtf,ktf, &
  1230. its,ite, jts,jte, kts,kte)
  1231. call cup_up_aa0(aa3,z3,zu3,dby3,GAMMA3_CUP,t3_cup, &
  1232. kbcon3,ktop3,ierr5, &
  1233. itf,jtf,ktf, &
  1234. its,ite, jts,jte, kts,kte)
  1235. ! do i=its,itf
  1236. ! if(ierr5(i).eq.0)then
  1237. ! if(aa3(i).eq.0.)then
  1238. ! ierr5(i)=17
  1239. ! endif
  1240. ! endif
  1241. ! enddo
  1242. ! call cup_dellabot('shallow',ipr,jpr,q3_cup,ierr5,z3_cup,po,qrcdo,edto, &
  1243. ! zdo,cdd,q3,dellaq3,dsubq,j,mentrd_rate,z3,g,&
  1244. ! itf,jtf,ktf, &
  1245. ! its,ite, jts,jte, kts,kte)
  1246. call cup_dellas_3d(ierr5,z3_cup,po_cup,hcdo,edt3,zdo3,cdd, &
  1247. he3,dellah3,dsubt3,j,mentrd_rate,zu3,g, &
  1248. cd3,hc3,ktop3,k23,kbcon3,mentr_rate3,jmin,he3_cup,kdet, &
  1249. k23,ipr,jpr,'shallow',0, &
  1250. itf,jtf,ktf, &
  1251. its,ite, jts,jte, kts,kte)
  1252. call cup_dellas_3d(ierr5,z3_cup,po_cup,qrcdo,edt3,zdo3,cdd, &
  1253. qshall,dellaq3,dsubq3,j,mentrd_rate,zu3,g, &
  1254. cd3,qc3,ktop3,k23,kbcon3,mentr_rate3,jmin,q3_cup,kdet, &
  1255. k23,ipr,jpr,'shallow',0, &
  1256. itf,jtf,ktf, &
  1257. its,ite, jts,jte, kts,kte )
  1258. mbdt_s=1.e-1*mbdt_ens(1)
  1259. do k=kts,ktf
  1260. do i=its,itf
  1261. dellat3(i,k)=0.
  1262. if(ierr5(i).eq.0)then
  1263. trash=dsubt3(i,k)
  1264. XHE3(I,K)=(dsubt3(i,k)+DELLAH3(I,K))*MBDT_S+HE3(I,K)
  1265. XQ3(I,K)=(dsubq3(i,k)+DELLAQ3(I,K))*MBDT_S+QSHALL(I,K)
  1266. DELLAT3(I,K)=(1./cp)*(DELLAH3(I,K)-xl*DELLAQ3(I,K))
  1267. dSUBT3(I,K)=(1./cp)*(dsubt3(i,k)-xl*dsubq3(i,k))
  1268. XT3(I,K)= (DELLAT3(I,K)+dsubt3(i,k))*MBDT_S+TSHALL(I,K)
  1269. IF(XQ3(I,K).LE.0.)XQ3(I,K)=1.E-08
  1270. ! if(i.eq.ipr.and.j.eq.jpr)then
  1271. ! write(0,*)k,trash,DELLAQ3(I,K),dsubq3(I,K),dsubt3(i,k)
  1272. ! endif
  1273. ENDIF
  1274. enddo
  1275. enddo
  1276. do i=its,itf
  1277. if(ierr5(i).eq.0)then
  1278. XHE3(I,ktf)=HE3(I,ktf)
  1279. XQ3(I,ktf)=QSHALL(I,ktf)
  1280. XT3(I,ktf)=TSHALL(I,ktf)
  1281. IF(XQ3(I,ktf).LE.0.)XQ3(I,ktf)=1.E-08
  1282. endif
  1283. enddo
  1284. !
  1285. !--- calculate moist static energy, heights, qes
  1286. !
  1287. call cup_env(xz3,xqes3,xhe3,xhes3,xt3,xq3,po,z1, &
  1288. psur,ierr5,tcrit,2,xl,cp, &
  1289. itf,jtf,ktf, &
  1290. its,ite, jts,jte, kts,kte)
  1291. !
  1292. !--- environmental values on cloud levels
  1293. !
  1294. call cup_env_clev(xt3,xqes3,xq3,xhe3,xhes3,xz3,po,xqes3_cup,xq3_cup, &
  1295. xhe3_cup,xhes3_cup,xz3_cup,po_cup,gamma3_cup,xt3_cup,psur, &
  1296. ierr5,z1,xl,rv,cp, &
  1297. itf,jtf,ktf, &
  1298. its,ite, jts,jte, kts,kte)
  1299. !
  1300. !
  1301. !**************************** static control
  1302. !
  1303. !--- moist static energy inside cloud
  1304. !
  1305. do i=its,itf
  1306. if(ierr5(i).eq.0)then
  1307. xhkb3(i)=xhe3(i,k23(i))
  1308. endif
  1309. enddo
  1310. call cup_up_he(k23,xhkb3,xz3_cup,cd3,mentr_rate3,xhe3_cup,xhc3, &
  1311. kbcon3,ierr5,xdby3,xhe3,xhes3_cup,'shallow', &
  1312. itf,jtf,ktf, &
  1313. its,ite, jts,jte, kts,kte)
  1314. !
  1315. !c--- normalized mass flux profile and CWF
  1316. !
  1317. call cup_up_nms(xzu3,xz3_cup,mentr_rate3,cd3,kbcon3,ktop3,ierr5,k23, &
  1318. itf,jtf,ktf, &
  1319. its,ite, jts,jte, kts,kte)
  1320. call cup_up_aa0(xaa3,xz3,xzu3,xdby3,GAMMA3_CUP,xt3_cup, &
  1321. kbcon3,ktop3,ierr5, &
  1322. itf,jtf,ktf, &
  1323. its,ite, jts,jte, kts,kte)
  1324. !
  1325. ! now for shallow forcing
  1326. !
  1327. do i=its,itf
  1328. xmb3(i)=0.
  1329. xff_shal(1:9)=0.
  1330. if(ierr5(i).eq.0)then
  1331. xkshal=(xaa3(i)-aa3(i))/mbdt_s
  1332. if(xkshal.ge.0.)xkshal=+1.e6
  1333. if(xkshal.gt.-1.e-4 .and. xkshal.lt.0.)xkshal=-1.e-4
  1334. xff_shal(1)=max(0.,-(aa3(i)-aa3_0(i))/(xkshal*dtime))
  1335. xff_shal(2)=max(0.,-(aa3(i)-aa3_

Large files files are truncated, but you can click here to view the full file