PageRenderTime 93ms CodeModel.GetById 26ms RepoModel.GetById 0ms app.codeStats 1ms

/wrfv2_fire/phys/module_cu_gd.F

http://github.com/jbeezley/wrf-fire
FORTRAN Legacy | 4555 lines | 3062 code | 317 blank | 1176 comment | 192 complexity | f99ad882e1f460753dbe05b590e319d4 MD5 | raw file
Possible License(s): AGPL-1.0
  1. !WRF:MODEL_LAYER:PHYSICS
  2. !
  3. MODULE module_cu_gd
  4. CONTAINS
  5. !-------------------------------------------------------------
  6. SUBROUTINE GRELLDRV( &
  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 &
  16. ,GDC,GDC2 &
  17. ,ensdim,maxiens,maxens,maxens2,maxens3 &
  18. ,ids,ide, jds,jde, kds,kde &
  19. ,ims,ime, jms,jme, kms,kme &
  20. ,its,ite, jts,jte, kts,kte &
  21. ,periodic_x,periodic_y &
  22. ,RQVCUTEN,RQCCUTEN,RQICUTEN &
  23. ,RQVFTEN,RQVBLTEN &
  24. ,RTHFTEN,RTHCUTEN,RTHRATEN,RTHBLTEN &
  25. ,F_QV ,F_QC ,F_QR ,F_QI ,F_QS &
  26. ,CFU1,CFD1,DFU1,EFU1,DFD1,EFD1,f_flux )
  27. !-------------------------------------------------------------
  28. IMPLICIT NONE
  29. !-------------------------------------------------------------
  30. INTEGER, INTENT(IN ) :: &
  31. ids,ide, jds,jde, kds,kde, &
  32. ims,ime, jms,jme, kms,kme, &
  33. its,ite, jts,jte, kts,kte
  34. LOGICAL periodic_x,periodic_y
  35. integer, intent (in ) :: &
  36. ensdim,maxiens,maxens,maxens2,maxens3
  37. INTEGER, INTENT(IN ) :: STEPCU, ITIMESTEP
  38. LOGICAL, INTENT(IN ) :: warm_rain
  39. REAL, INTENT(IN ) :: XLV, R_v
  40. REAL, INTENT(IN ) :: CP,G
  41. REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , &
  42. INTENT(IN ) :: &
  43. U, &
  44. V, &
  45. W, &
  46. pi, &
  47. t, &
  48. q, &
  49. p, &
  50. dz8w, &
  51. p8w, &
  52. rho
  53. REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , &
  54. OPTIONAL , &
  55. INTENT(INOUT ) :: &
  56. GDC,GDC2
  57. REAL, DIMENSION( ims:ime , jms:jme ),INTENT(IN) :: GSW,HT,XLAND
  58. !
  59. REAL, INTENT(IN ) :: DT, DX
  60. !
  61. REAL, DIMENSION( ims:ime , jms:jme ), &
  62. INTENT(INOUT) :: RAINCV, PRATEC, MASS_FLUX, &
  63. APR_GR,APR_W,APR_MC,APR_ST,APR_AS, &
  64. APR_CAPMA,APR_CAPME,APR_CAPMI,htop,hbot
  65. !+lxz
  66. ! REAL, DIMENSION( ims:ime , jms:jme ) :: & !, INTENT(INOUT) :: &
  67. ! HTOP, &! highest model layer penetrated by cumulus since last reset in radiation_driver
  68. ! HBOT ! lowest model layer penetrated by cumulus since last reset in radiation_driver
  69. ! ! HBOT>HTOP follow physics leveling convention
  70. LOGICAL, DIMENSION( ims:ime , jms:jme ), &
  71. INTENT(INOUT) :: CU_ACT_FLAG
  72. !
  73. ! Optionals
  74. !
  75. REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
  76. OPTIONAL, &
  77. INTENT(INOUT) :: RTHFTEN, &
  78. RQVFTEN
  79. REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
  80. OPTIONAL, &
  81. INTENT(IN ) :: &
  82. RTHRATEN, &
  83. RTHBLTEN, &
  84. RQVBLTEN
  85. REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
  86. OPTIONAL, &
  87. INTENT(INOUT) :: &
  88. RTHCUTEN, &
  89. RQVCUTEN, &
  90. RQCCUTEN, &
  91. RQICUTEN
  92. REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
  93. OPTIONAL, &
  94. INTENT(INOUT) :: &
  95. CFU1, &
  96. CFD1, &
  97. DFU1, &
  98. EFU1, &
  99. DFD1, &
  100. EFD1
  101. !
  102. ! Flags relating to the optional tendency arrays declared above
  103. ! Models that carry the optional tendencies will provdide the
  104. ! optional arguments at compile time; these flags all the model
  105. ! to determine at run-time whether a particular tracer is in
  106. ! use or not.
  107. !
  108. LOGICAL, OPTIONAL :: &
  109. F_QV &
  110. ,F_QC &
  111. ,F_QR &
  112. ,F_QI &
  113. ,F_QS
  114. LOGICAL, intent(in), OPTIONAL :: f_flux
  115. ! LOCAL VARS
  116. real, dimension ( ims:ime , jms:jme , 1:ensdim) :: &
  117. massfln,xf_ens,pr_ens
  118. real, dimension (its:ite,kts:kte+1) :: &
  119. OUTT,OUTQ,OUTQC,phh,cupclw, &
  120. outCFU1,outCFD1,outDFU1,outEFU1,outDFD1,outEFD1
  121. logical :: l_flux
  122. real, dimension (its:ite) :: &
  123. pret, ter11, aa0, fp
  124. !+lxz
  125. integer, dimension (its:ite) :: &
  126. kbcon, ktop
  127. !.lxz
  128. integer, dimension (its:ite,jts:jte) :: &
  129. iact_old_gr
  130. integer :: ichoice,iens,ibeg,iend,jbeg,jend
  131. !
  132. ! basic environmental input includes moisture convergence (mconv)
  133. ! omega (omeg), windspeed (us,vs), and a flag (aaeq) to turn off
  134. ! convection for this call only and at that particular gridpoint
  135. !
  136. real, dimension (its:ite,kts:kte+1) :: &
  137. T2d,TN,q2d,qo,PO,P2d,US,VS,omeg
  138. real, dimension (its:ite) :: &
  139. Z1,PSUR,AAEQ,direction,mconv,cuten,umean,vmean,pmean
  140. INTEGER :: i,j,k,ICLDCK,ipr,jpr
  141. REAL :: tcrit,dp,dq
  142. INTEGER :: itf,jtf,ktf
  143. REAL :: rkbcon,rktop !-lxz
  144. l_flux=.FALSE.
  145. if (present(f_flux)) l_flux=f_flux
  146. if (l_flux) then
  147. l_flux = l_flux .and. present(cfu1) .and. present(cfd1) &
  148. .and. present(dfu1) .and. present(efu1) &
  149. .and. present(dfd1) .and. present(efd1)
  150. endif
  151. ichoice=0
  152. iens=1
  153. ipr=0
  154. jpr=0
  155. IF ( periodic_x ) THEN
  156. ibeg=max(its,ids)
  157. iend=min(ite,ide-1)
  158. ELSE
  159. ibeg=max(its,ids+4)
  160. iend=min(ite,ide-5)
  161. END IF
  162. IF ( periodic_y ) THEN
  163. jbeg=max(jts,jds)
  164. jend=min(jte,jde-1)
  165. ELSE
  166. jbeg=max(jts,jds+4)
  167. jend=min(jte,jde-5)
  168. END IF
  169. tcrit=258.
  170. itf=MIN(ite,ide-1)
  171. ktf=MIN(kte,kde-1)
  172. jtf=MIN(jte,jde-1)
  173. !
  174. DO 100 J = jts,jtf
  175. DO I= its,itf
  176. cuten(i)=0.
  177. iact_old_gr(i,j)=0
  178. mass_flux(i,j)=0.
  179. pratec(i,j) = 0.
  180. raincv(i,j)=0.
  181. CU_ACT_FLAG(i,j) = .true.
  182. ENDDO
  183. DO k=1,ensdim
  184. DO I= its,itf
  185. massfln(i,j,k)=0.
  186. ENDDO
  187. ENDDO
  188. #if ( EM_CORE == 1 )
  189. DO k= kts,ktf
  190. DO I= its,itf
  191. RTHFTEN(i,k,j)=(RTHFTEN(i,k,j)+RTHRATEN(i,k,j)+RTHBLTEN(i,k,j))*pi(i,k,j)
  192. RQVFTEN(i,k,j)=RQVFTEN(i,k,j)+RQVBLTEN(i,k,j)
  193. ENDDO
  194. ENDDO
  195. #endif
  196. ! put hydrostatic pressure on half levels
  197. DO K=kts,ktf
  198. DO I=ITS,ITF
  199. phh(i,k) = p(i,k,j)
  200. ENDDO
  201. ENDDO
  202. DO I=ITS,ITF
  203. PSUR(I)=p8w(I,1,J)*.01
  204. TER11(I)=HT(i,j)
  205. mconv(i)=0.
  206. aaeq(i)=0.
  207. direction(i)=0.
  208. pret(i)=0.
  209. umean(i)=0.
  210. vmean(i)=0.
  211. pmean(i)=0.
  212. ENDDO
  213. DO K=kts,ktf
  214. DO I=ITS,ITF
  215. omeg(i,k)=0.
  216. ! cupclw(i,k)=0.
  217. po(i,k)=phh(i,k)*.01
  218. P2d(I,K)=PO(i,k)
  219. US(I,K) =u(i,k,j)
  220. VS(I,K) =v(i,k,j)
  221. T2d(I,K)=t(i,k,j)
  222. q2d(I,K)=q(i,k,j)
  223. omeg(I,K)= -g*rho(i,k,j)*w(i,k,j)
  224. TN(I,K)=t2d(i,k)+RTHFTEN(i,k,j)*dt
  225. IF(TN(I,K).LT.200.)TN(I,K)=T2d(I,K)
  226. QO(I,K)=q2d(i,k)+RQVFTEN(i,k,j)*dt
  227. IF(Q2d(I,K).LT.1.E-08)Q2d(I,K)=1.E-08
  228. IF(QO(I,K).LT.1.E-08)QO(I,K)=1.E-08
  229. OUTT(I,K)=0.
  230. OUTQ(I,K)=0.
  231. OUTQC(I,K)=0.
  232. ! RTHFTEN(i,k,j)=0.
  233. ! RQVFTEN(i,k,j)=0.
  234. ENDDO
  235. ENDDO
  236. do k= kts+1,ktf-1
  237. DO I = its,itf
  238. if((p2d(i,1)-p2d(i,k)).gt.150.and.p2d(i,k).gt.300)then
  239. dp=-.5*(p2d(i,k+1)-p2d(i,k-1))
  240. umean(i)=umean(i)+us(i,k)*dp
  241. vmean(i)=vmean(i)+vs(i,k)*dp
  242. pmean(i)=pmean(i)+dp
  243. endif
  244. enddo
  245. enddo
  246. DO I = its,itf
  247. if(pmean(i).gt.0)then
  248. umean(i)=umean(i)/pmean(i)
  249. vmean(i)=vmean(i)/pmean(i)
  250. direction(i)=(atan2(umean(i),vmean(i))+3.1415926)*57.29578
  251. if(direction(i).gt.360.)direction(i)=direction(i)-360.
  252. endif
  253. ENDDO
  254. DO K=kts,ktf-1
  255. DO I = its,itf
  256. dq=(q2d(i,k+1)-q2d(i,k))
  257. mconv(i)=mconv(i)+omeg(i,k)*dq/g
  258. ENDDO
  259. ENDDO
  260. DO I = its,itf
  261. if(mconv(i).lt.0.)mconv(i)=0.
  262. ENDDO
  263. !
  264. !---- CALL CUMULUS PARAMETERIZATION
  265. !
  266. CALL CUP_enss(outqc,j,AAEQ,T2d,Q2d,TER11,TN,QO,PO,PRET, &
  267. P2d,OUTT,OUTQ,DT,PSUR,US,VS,tcrit,iens, &
  268. mconv,massfln,iact_old_gr,omeg,direction,MASS_FLUX, &
  269. maxiens,maxens,maxens2,maxens3,ensdim, &
  270. APR_GR,APR_W,APR_MC,APR_ST,APR_AS, &
  271. APR_CAPMA,APR_CAPME,APR_CAPMI,kbcon,ktop, &
  272. xf_ens,pr_ens,XLAND,gsw,cupclw, &
  273. xlv,r_v,cp,g,ichoice,ipr,jpr, &
  274. outCFU1,outCFD1,outDFU1,outEFU1,outDFD1,outEFD1,l_flux,&
  275. ids,ide, jds,jde, kds,kde, &
  276. ims,ime, jms,jme, kms,kme, &
  277. its,ite, jts,jte, kts,kte )
  278. CALL neg_check(dt,q2d,outq,outt,outqc,pret,its,ite,kts,kte,itf,ktf)
  279. if(j.ge.jbeg.and.j.le.jend)then
  280. DO I=its,itf
  281. ! cuten(i)=0.
  282. if(i.ge.ibeg.and.i.le.iend)then
  283. if(pret(i).gt.0.)then
  284. pratec(i,j)=pret(i)
  285. raincv(i,j)=pret(i)*dt
  286. cuten(i)=1.
  287. rkbcon = kte+kts - kbcon(i)
  288. rktop = kte+kts - ktop(i)
  289. if (ktop(i) > HTOP(i,j)) HTOP(i,j) = ktop(i)+.001
  290. if (kbcon(i) < HBOT(i,j)) HBOT(i,j) = kbcon(i)+.001
  291. endif
  292. else
  293. pret(i)=0.
  294. endif
  295. ENDDO
  296. DO K=kts,ktf
  297. DO I=its,itf
  298. RTHCUTEN(I,K,J)=outt(i,k)*cuten(i)/pi(i,k,j)
  299. RQVCUTEN(I,K,J)=outq(i,k)*cuten(i)
  300. ENDDO
  301. ENDDO
  302. IF(PRESENT(RQCCUTEN)) THEN
  303. IF ( F_QC ) THEN
  304. DO K=kts,ktf
  305. DO I=its,itf
  306. RQCCUTEN(I,K,J)=outqc(I,K)*cuten(i)
  307. IF ( PRESENT( GDC ) ) GDC(I,K,J)=CUPCLW(I,K)*cuten(i)
  308. IF ( PRESENT( GDC2 ) ) GDC2(I,K,J)=0.
  309. ENDDO
  310. ENDDO
  311. ENDIF
  312. ENDIF
  313. !...... QSTEN STORES GRAUPEL TENDENCY IF IT EXISTS, OTHERISE SNOW (V2)
  314. IF(PRESENT(RQICUTEN).AND.PRESENT(RQCCUTEN))THEN
  315. IF (F_QI) THEN
  316. DO K=kts,ktf
  317. DO I=its,itf
  318. if(t2d(i,k).lt.258.)then
  319. RQICUTEN(I,K,J)=outqc(I,K)*cuten(i)
  320. RQCCUTEN(I,K,J)=0.
  321. IF ( PRESENT( GDC2 ) ) GDC2(I,K,J)=CUPCLW(I,K)*cuten(i)
  322. else
  323. RQICUTEN(I,K,J)=0.
  324. RQCCUTEN(I,K,J)=outqc(I,K)*cuten(i)
  325. IF ( PRESENT( GDC ) ) GDC(I,K,J)=CUPCLW(I,K)*cuten(i)
  326. endif
  327. ENDDO
  328. ENDDO
  329. ENDIF
  330. ENDIF
  331. if (l_flux) then
  332. DO K=kts,ktf
  333. DO I=its,itf
  334. cfu1(i,k,j)=outcfu1(i,k)*cuten(i)
  335. cfd1(i,k,j)=outcfd1(i,k)*cuten(i)
  336. dfu1(i,k,j)=outdfu1(i,k)*cuten(i)
  337. efu1(i,k,j)=outefu1(i,k)*cuten(i)
  338. dfd1(i,k,j)=outdfd1(i,k)*cuten(i)
  339. efd1(i,k,j)=outefd1(i,k)*cuten(i)
  340. enddo
  341. enddo
  342. endif
  343. endif !jbeg,jend
  344. 100 continue
  345. END SUBROUTINE GRELLDRV
  346. SUBROUTINE CUP_enss(OUTQC,J,AAEQ,T,Q,Z1, &
  347. TN,QO,PO,PRE,P,OUTT,OUTQ,DTIME,PSUR,US,VS, &
  348. TCRIT,iens,mconv,massfln,iact, &
  349. omeg,direction,massflx,maxiens, &
  350. maxens,maxens2,maxens3,ensdim, &
  351. APR_GR,APR_W,APR_MC,APR_ST,APR_AS, &
  352. APR_CAPMA,APR_CAPME,APR_CAPMI,kbcon,ktop, & !-lxz
  353. xf_ens,pr_ens,xland,gsw,cupclw, &
  354. xl,rv,cp,g,ichoice,ipr,jpr, &
  355. outCFU1,outCFD1,outDFU1,outEFU1,outDFD1,outEFD1,l_flux, &
  356. ids,ide, jds,jde, kds,kde, &
  357. ims,ime, jms,jme, kms,kme, &
  358. its,ite, jts,jte, kts,kte )
  359. IMPLICIT NONE
  360. integer &
  361. ,intent (in ) :: &
  362. ids,ide, jds,jde, kds,kde, &
  363. ims,ime, jms,jme, kms,kme, &
  364. its,ite, jts,jte, kts,kte,ipr,jpr
  365. integer, intent (in ) :: &
  366. j,ensdim,maxiens,maxens,maxens2,maxens3,ichoice,iens
  367. !
  368. !
  369. !
  370. real, dimension (ims:ime,jms:jme,1:ensdim) &
  371. ,intent (inout) :: &
  372. massfln,xf_ens,pr_ens
  373. real, dimension (ims:ime,jms:jme) &
  374. ,intent (inout ) :: &
  375. APR_GR,APR_W,APR_MC,APR_ST,APR_AS,APR_CAPMA, &
  376. APR_CAPME,APR_CAPMI,massflx
  377. real, dimension (ims:ime,jms:jme) &
  378. ,intent (in ) :: &
  379. xland,gsw
  380. integer, dimension (its:ite,jts:jte) &
  381. ,intent (in ) :: &
  382. iact
  383. ! outtem = output temp tendency (per s)
  384. ! outq = output q tendency (per s)
  385. ! outqc = output qc tendency (per s)
  386. ! pre = output precip
  387. real, dimension (its:ite,kts:kte) &
  388. ,intent (out ) :: &
  389. OUTT,OUTQ,OUTQC,CUPCLW, &
  390. outCFU1,outCFD1,outDFU1,outEFU1,outDFD1,outEFD1
  391. logical, intent(in) :: l_flux
  392. real, dimension (its:ite) &
  393. ,intent (out ) :: &
  394. pre
  395. !+lxz
  396. integer, dimension (its:ite) &
  397. ,intent (out ) :: &
  398. kbcon,ktop
  399. !.lxz
  400. !
  401. ! basic environmental input includes moisture convergence (mconv)
  402. ! omega (omeg), windspeed (us,vs), and a flag (aaeq) to turn off
  403. ! convection for this call only and at that particular gridpoint
  404. !
  405. real, dimension (its:ite,kts:kte) &
  406. ,intent (in ) :: &
  407. T,TN,PO,P,US,VS,omeg
  408. real, dimension (its:ite,kts:kte) &
  409. ,intent (inout) :: &
  410. Q,QO
  411. real, dimension (its:ite) &
  412. ,intent (in ) :: &
  413. Z1,PSUR,AAEQ,direction,mconv
  414. real &
  415. ,intent (in ) :: &
  416. dtime,tcrit,xl,cp,rv,g
  417. !
  418. ! local ensemble dependent variables in this routine
  419. !
  420. real, dimension (its:ite,1:maxens) :: &
  421. xaa0_ens
  422. real, dimension (1:maxens) :: &
  423. mbdt_ens
  424. real, dimension (1:maxens2) :: &
  425. edt_ens
  426. real, dimension (its:ite,1:maxens2) :: &
  427. edtc
  428. real, dimension (its:ite,kts:kte,1:maxens2) :: &
  429. dellat_ens,dellaqc_ens,dellaq_ens,pwo_ens
  430. real, dimension (its:ite,kts:kte,1:maxens2) :: &
  431. CFU1_ens,CFD1_ens,DFU1_ens,EFU1_ens,DFD1_ens,EFD1_ens
  432. !
  433. !
  434. !
  435. !***************** the following are your basic environmental
  436. ! variables. They carry a "_cup" if they are
  437. ! on model cloud levels (staggered). They carry
  438. ! an "o"-ending (z becomes zo), if they are the forced
  439. ! variables. They are preceded by x (z becomes xz)
  440. ! to indicate modification by some typ of cloud
  441. !
  442. ! z = heights of model levels
  443. ! q = environmental mixing ratio
  444. ! qes = environmental saturation mixing ratio
  445. ! t = environmental temp
  446. ! p = environmental pressure
  447. ! he = environmental moist static energy
  448. ! hes = environmental saturation moist static energy
  449. ! z_cup = heights of model cloud levels
  450. ! q_cup = environmental q on model cloud levels
  451. ! qes_cup = saturation q on model cloud levels
  452. ! t_cup = temperature (Kelvin) on model cloud levels
  453. ! p_cup = environmental pressure
  454. ! he_cup = moist static energy on model cloud levels
  455. ! hes_cup = saturation moist static energy on model cloud levels
  456. ! gamma_cup = gamma on model cloud levels
  457. !
  458. !
  459. ! hcd = moist static energy in downdraft
  460. ! zd normalized downdraft mass flux
  461. ! dby = buoancy term
  462. ! entr = entrainment rate
  463. ! zd = downdraft normalized mass flux
  464. ! entr= entrainment rate
  465. ! hcd = h in model cloud
  466. ! bu = buoancy term
  467. ! zd = normalized downdraft mass flux
  468. ! gamma_cup = gamma on model cloud levels
  469. ! mentr_rate = entrainment rate
  470. ! qcd = cloud q (including liquid water) after entrainment
  471. ! qrch = saturation q in cloud
  472. ! pwd = evaporate at that level
  473. ! pwev = total normalized integrated evaoprate (I2)
  474. ! entr= entrainment rate
  475. ! z1 = terrain elevation
  476. ! entr = downdraft entrainment rate
  477. ! jmin = downdraft originating level
  478. ! kdet = level above ground where downdraft start detraining
  479. ! psur = surface pressure
  480. ! z1 = terrain elevation
  481. ! pr_ens = precipitation ensemble
  482. ! xf_ens = mass flux ensembles
  483. ! massfln = downdraft mass flux ensembles used in next timestep
  484. ! omeg = omega from large scale model
  485. ! mconv = moisture convergence from large scale model
  486. ! zd = downdraft normalized mass flux
  487. ! zu = updraft normalized mass flux
  488. ! dir = "storm motion"
  489. ! mbdt = arbitrary numerical parameter
  490. ! dtime = dt over which forcing is applied
  491. ! iact_gr_old = flag to tell where convection was active
  492. ! kbcon = LFC of parcel from k22
  493. ! k22 = updraft originating level
  494. ! icoic = flag if only want one closure (usually set to zero!)
  495. ! dby = buoancy term
  496. ! ktop = cloud top (output)
  497. ! xmb = total base mass flux
  498. ! hc = cloud moist static energy
  499. ! hkb = moist static energy at originating level
  500. ! mentr_rate = entrainment rate
  501. real, dimension (its:ite,kts:kte) :: &
  502. he,hes,qes,z, &
  503. heo,heso,qeso,zo, &
  504. xhe,xhes,xqes,xz,xt,xq, &
  505. qes_cup,q_cup,he_cup,hes_cup,z_cup,p_cup,gamma_cup,t_cup, &
  506. qeso_cup,qo_cup,heo_cup,heso_cup,zo_cup,po_cup,gammao_cup, &
  507. tn_cup, &
  508. xqes_cup,xq_cup,xhe_cup,xhes_cup,xz_cup,xp_cup,xgamma_cup, &
  509. xt_cup, &
  510. dby,qc,qrcd,pwd,pw,hcd,qcd,dbyd,hc,qrc,zu,zd,clw_all, &
  511. dbyo,qco,qrcdo,pwdo,pwo,hcdo,qcdo,dbydo,hco,qrco,zuo,zdo, &
  512. xdby,xqc,xqrcd,xpwd,xpw,xhcd,xqcd,xhc,xqrc,xzu,xzd, &
  513. ! cd = detrainment function for updraft
  514. ! cdd = detrainment function for downdraft
  515. ! dellat = change of temperature per unit mass flux of cloud ensemble
  516. ! dellaq = change of q per unit mass flux of cloud ensemble
  517. ! dellaqc = change of qc per unit mass flux of cloud ensemble
  518. cd,cdd,scr1,DELLAH,DELLAQ,DELLAT,DELLAQC, &
  519. CFU1,CFD1,DFU1,EFU1,DFD1,EFD1
  520. ! aa0 cloud work function for downdraft
  521. ! edt = epsilon
  522. ! aa0 = cloud work function without forcing effects
  523. ! aa1 = cloud work function with forcing effects
  524. ! xaa0 = cloud work function with cloud effects (ensemble dependent)
  525. ! edt = epsilon
  526. real, dimension (its:ite) :: &
  527. edt,edto,edtx,AA1,AA0,XAA0,HKB,HKBO,aad,XHKB,QKB,QKBO, &
  528. XMB,XPWAV,XPWEV,PWAV,PWEV,PWAVO,PWEVO,BU,BUO,cap_max,xland1, &
  529. cap_max_increment,closure_n
  530. integer, dimension (its:ite) :: &
  531. kzdown,KDET,K22,KB,JMIN,kstabi,kstabm,K22x, & !-lxz
  532. KBCONx,KBx,KTOPx,ierr,ierr2,ierr3,KBMAX
  533. integer :: &
  534. nall,iedt,nens,nens3,ki,I,K,KK,iresult
  535. real :: &
  536. day,dz,mbdt,entr_rate,radius,entrd_rate,mentr_rate,mentrd_rate, &
  537. zcutdown,edtmax,edtmin,depth_min,zkbmax,z_detr,zktop, &
  538. massfld,dh,cap_maxs
  539. integer :: itf,jtf,ktf
  540. integer :: jmini
  541. logical :: keep_going
  542. itf=MIN(ite,ide-1)
  543. ktf=MIN(kte,kde-1)
  544. jtf=MIN(jte,jde-1)
  545. !sms$distribute end
  546. day=86400.
  547. do i=its,itf
  548. closure_n(i)=16.
  549. xland1(i)=1.
  550. if(xland(i,j).gt.1.5)xland1(i)=0.
  551. cap_max_increment(i)=25.
  552. enddo
  553. !
  554. !--- specify entrainmentrate and detrainmentrate
  555. !
  556. if(iens.le.4)then
  557. radius=14000.-float(iens)*2000.
  558. else
  559. radius=12000.
  560. endif
  561. !
  562. !--- gross entrainment rate (these may be changed later on in the
  563. !--- program, depending what your detrainment is!!)
  564. !
  565. entr_rate=.2/radius
  566. !
  567. !--- entrainment of mass
  568. !
  569. mentrd_rate=0.
  570. mentr_rate=entr_rate
  571. !
  572. !--- initial detrainmentrates
  573. !
  574. do k=kts,ktf
  575. do i=its,itf
  576. cupclw(i,k)=0.
  577. cd(i,k)=0.1*entr_rate
  578. cdd(i,k)=0.
  579. enddo
  580. enddo
  581. !
  582. !--- max/min allowed value for epsilon (ratio downdraft base mass flux/updraft
  583. ! base mass flux
  584. !
  585. edtmax=.8
  586. edtmin=.2
  587. !
  588. !--- minimum depth (m), clouds must have
  589. !
  590. depth_min=500.
  591. !
  592. !--- maximum depth (mb) of capping
  593. !--- inversion (larger cap = no convection)
  594. !
  595. cap_maxs=75.
  596. !sms$to_local(grid_dh: <1, mix :size>, <2, mjx :size>) begin
  597. DO 7 i=its,itf
  598. kbmax(i)=1
  599. aa0(i)=0.
  600. aa1(i)=0.
  601. aad(i)=0.
  602. edt(i)=0.
  603. kstabm(i)=ktf-1
  604. IERR(i)=0
  605. IERR2(i)=0
  606. IERR3(i)=0
  607. if(aaeq(i).lt.-1.)then
  608. ierr(i)=20
  609. endif
  610. 7 CONTINUE
  611. !
  612. !--- first check for upstream convection
  613. !
  614. do i=its,itf
  615. cap_max(i)=cap_maxs
  616. ! if(tkmax(i,j).lt.2.)cap_max(i)=25.
  617. if(gsw(i,j).lt.1.)cap_max(i)=25.
  618. iresult=0
  619. ! massfld=0.
  620. ! call cup_direction2(i,j,direction,iact, &
  621. ! cu_mfx,iresult,0,massfld, &
  622. ! ids,ide, jds,jde, kds,kde, &
  623. ! ims,ime, jms,jme, kms,kme, &
  624. ! its,ite, jts,jte, kts,kte)
  625. ! cap_max(i)=cap_maxs
  626. if(iresult.eq.1)then
  627. cap_max(i)=cap_maxs+20.
  628. endif
  629. ! endif
  630. enddo
  631. !
  632. !--- max height(m) above ground where updraft air can originate
  633. !
  634. zkbmax=4000.
  635. !
  636. !--- height(m) above which no downdrafts are allowed to originate
  637. !
  638. zcutdown=3000.
  639. !
  640. !--- depth(m) over which downdraft detrains all its mass
  641. !
  642. z_detr=1250.
  643. !
  644. do nens=1,maxens
  645. mbdt_ens(nens)=(float(nens)-3.)*dtime*1.e-3+dtime*5.E-03
  646. enddo
  647. do nens=1,maxens2
  648. edt_ens(nens)=.95-float(nens)*.01
  649. enddo
  650. ! if(j.eq.jpr)then
  651. ! print *,'radius ensemble ',iens,radius
  652. ! print *,mbdt_ens
  653. ! print *,edt_ens
  654. ! endif
  655. !
  656. !--- environmental conditions, FIRST HEIGHTS
  657. !
  658. do i=its,itf
  659. if(ierr(i).ne.20)then
  660. do k=1,maxens*maxens2*maxens3
  661. xf_ens(i,j,(iens-1)*maxens*maxens2*maxens3+k)=0.
  662. pr_ens(i,j,(iens-1)*maxens*maxens2*maxens3+k)=0.
  663. enddo
  664. endif
  665. enddo
  666. !
  667. !--- calculate moist static energy, heights, qes
  668. !
  669. call cup_env(z,qes,he,hes,t,q,p,z1, &
  670. psur,ierr,tcrit,0,xl,cp, &
  671. ids,ide, jds,jde, kds,kde, &
  672. ims,ime, jms,jme, kms,kme, &
  673. its,ite, jts,jte, kts,kte)
  674. call cup_env(zo,qeso,heo,heso,tn,qo,po,z1, &
  675. psur,ierr,tcrit,0,xl,cp, &
  676. ids,ide, jds,jde, kds,kde, &
  677. ims,ime, jms,jme, kms,kme, &
  678. its,ite, jts,jte, kts,kte)
  679. !
  680. !--- environmental values on cloud levels
  681. !
  682. call cup_env_clev(t,qes,q,he,hes,z,p,qes_cup,q_cup,he_cup, &
  683. hes_cup,z_cup,p_cup,gamma_cup,t_cup,psur, &
  684. ierr,z1,xl,rv,cp, &
  685. ids,ide, jds,jde, kds,kde, &
  686. ims,ime, jms,jme, kms,kme, &
  687. its,ite, jts,jte, kts,kte)
  688. call cup_env_clev(tn,qeso,qo,heo,heso,zo,po,qeso_cup,qo_cup, &
  689. heo_cup,heso_cup,zo_cup,po_cup,gammao_cup,tn_cup,psur, &
  690. ierr,z1,xl,rv,cp, &
  691. ids,ide, jds,jde, kds,kde, &
  692. ims,ime, jms,jme, kms,kme, &
  693. its,ite, jts,jte, kts,kte)
  694. do i=its,itf
  695. if(ierr(i).eq.0)then
  696. !
  697. do k=kts,ktf-2
  698. if(zo_cup(i,k).gt.zkbmax+z1(i))then
  699. kbmax(i)=k
  700. go to 25
  701. endif
  702. enddo
  703. 25 continue
  704. !
  705. !
  706. !--- level where detrainment for downdraft starts
  707. !
  708. do k=kts,ktf
  709. if(zo_cup(i,k).gt.z_detr+z1(i))then
  710. kdet(i)=k
  711. go to 26
  712. endif
  713. enddo
  714. 26 continue
  715. !
  716. endif
  717. enddo
  718. !
  719. !
  720. !
  721. !------- DETERMINE LEVEL WITH HIGHEST MOIST STATIC ENERGY CONTENT - K22
  722. !
  723. CALL cup_MAXIMI(HEO_CUP,3,KBMAX,K22,ierr, &
  724. ids,ide, jds,jde, kds,kde, &
  725. ims,ime, jms,jme, kms,kme, &
  726. its,ite, jts,jte, kts,kte)
  727. DO 36 i=its,itf
  728. IF(ierr(I).eq.0.)THEN
  729. IF(K22(I).GE.KBMAX(i))ierr(i)=2
  730. endif
  731. 36 CONTINUE
  732. !
  733. !--- DETERMINE THE LEVEL OF CONVECTIVE CLOUD BASE - KBCON
  734. !
  735. call cup_kbcon(cap_max_increment,1,k22,kbcon,heo_cup,heso_cup, &
  736. ierr,kbmax,po_cup,cap_max, &
  737. ids,ide, jds,jde, kds,kde, &
  738. ims,ime, jms,jme, kms,kme, &
  739. its,ite, jts,jte, kts,kte)
  740. ! call cup_kbcon_cin(1,k22,kbcon,heo_cup,heso_cup,z,tn_cup, &
  741. ! qeso_cup,ierr,kbmax,po_cup,cap_max,xl,cp,&
  742. ! ids,ide, jds,jde, kds,kde, &
  743. ! ims,ime, jms,jme, kms,kme, &
  744. ! its,ite, jts,jte, kts,kte)
  745. !
  746. !--- increase detrainment in stable layers
  747. !
  748. CALL cup_minimi(HEso_cup,Kbcon,kstabm,kstabi,ierr, &
  749. ids,ide, jds,jde, kds,kde, &
  750. ims,ime, jms,jme, kms,kme, &
  751. its,ite, jts,jte, kts,kte)
  752. do i=its,itf
  753. IF(ierr(I).eq.0.)THEN
  754. if(kstabm(i)-1.gt.kstabi(i))then
  755. do k=kstabi(i),kstabm(i)-1
  756. cd(i,k)=cd(i,k-1)+1.5*entr_rate
  757. if(cd(i,k).gt.10.0*entr_rate)cd(i,k)=10.0*entr_rate
  758. enddo
  759. ENDIF
  760. ENDIF
  761. ENDDO
  762. !
  763. !--- calculate incloud moist static energy
  764. !
  765. call cup_up_he(k22,hkb,z_cup,cd,mentr_rate,he_cup,hc, &
  766. kbcon,ierr,dby,he,hes_cup, &
  767. ids,ide, jds,jde, kds,kde, &
  768. ims,ime, jms,jme, kms,kme, &
  769. its,ite, jts,jte, kts,kte)
  770. call cup_up_he(k22,hkbo,zo_cup,cd,mentr_rate,heo_cup,hco, &
  771. kbcon,ierr,dbyo,heo,heso_cup, &
  772. ids,ide, jds,jde, kds,kde, &
  773. ims,ime, jms,jme, kms,kme, &
  774. its,ite, jts,jte, kts,kte)
  775. !--- DETERMINE CLOUD TOP - KTOP
  776. !
  777. call cup_ktop(1,dbyo,kbcon,ktop,ierr, &
  778. ids,ide, jds,jde, kds,kde, &
  779. ims,ime, jms,jme, kms,kme, &
  780. its,ite, jts,jte, kts,kte)
  781. DO 37 i=its,itf
  782. kzdown(i)=0
  783. if(ierr(i).eq.0)then
  784. zktop=(zo_cup(i,ktop(i))-z1(i))*.6
  785. zktop=min(zktop+z1(i),zcutdown+z1(i))
  786. do k=kts,kte
  787. if(zo_cup(i,k).gt.zktop)then
  788. kzdown(i)=k
  789. go to 37
  790. endif
  791. enddo
  792. endif
  793. 37 CONTINUE
  794. !
  795. !--- DOWNDRAFT ORIGINATING LEVEL - JMIN
  796. !
  797. call cup_minimi(HEso_cup,K22,kzdown,JMIN,ierr, &
  798. ids,ide, jds,jde, kds,kde, &
  799. ims,ime, jms,jme, kms,kme, &
  800. its,ite, jts,jte, kts,kte)
  801. DO 100 i=its,ite
  802. IF(ierr(I).eq.0.)THEN
  803. !
  804. !--- check whether it would have buoyancy, if there where
  805. !--- no entrainment/detrainment
  806. !
  807. !jm begin 20061212: the following code replaces code that
  808. ! was too complex and causing problem for optimization.
  809. ! Done in consultation with G. Grell.
  810. jmini = jmin(i)
  811. keep_going = .TRUE.
  812. DO WHILE ( keep_going )
  813. keep_going = .FALSE.
  814. if ( jmini - 1 .lt. kdet(i) ) kdet(i) = jmini-1
  815. if ( jmini .ge. ktop(i)-1 ) jmini = ktop(i) - 2
  816. ki = jmini
  817. hcdo(i,ki)=heso_cup(i,ki)
  818. DZ=Zo_cup(i,Ki+1)-Zo_cup(i,Ki)
  819. dh=0.
  820. DO k=ki-1,1,-1
  821. hcdo(i,k)=heso_cup(i,jmini)
  822. DZ=Zo_cup(i,K+1)-Zo_cup(i,K)
  823. dh=dh+dz*(HCDo(i,K)-heso_cup(i,k))
  824. IF(dh.gt.0.)THEN
  825. jmini=jmini-1
  826. IF ( jmini .gt. 3 ) THEN
  827. keep_going = .TRUE.
  828. ELSE
  829. ierr(i) = 9
  830. EXIT
  831. ENDIF
  832. ENDIF
  833. ENDDO
  834. ENDDO
  835. jmin(i) = jmini
  836. IF ( jmini .le. 3 ) THEN
  837. ierr(i)=4
  838. ENDIF
  839. !jm end 20061212
  840. ENDIF
  841. 100 CONTINUE
  842. !
  843. ! - Must have at least depth_min m between cloud convective base
  844. ! and cloud top.
  845. !
  846. do i=its,itf
  847. IF(ierr(I).eq.0.)THEN
  848. IF(-zo_cup(I,KBCON(I))+zo_cup(I,KTOP(I)).LT.depth_min)then
  849. ierr(i)=6
  850. endif
  851. endif
  852. enddo
  853. !
  854. !c--- normalized updraft mass flux profile
  855. !
  856. call cup_up_nms(zu,z_cup,mentr_rate,cd,kbcon,ktop,ierr,k22, &
  857. ids,ide, jds,jde, kds,kde, &
  858. ims,ime, jms,jme, kms,kme, &
  859. its,ite, jts,jte, kts,kte)
  860. call cup_up_nms(zuo,zo_cup,mentr_rate,cd,kbcon,ktop,ierr,k22, &
  861. ids,ide, jds,jde, kds,kde, &
  862. ims,ime, jms,jme, kms,kme, &
  863. its,ite, jts,jte, kts,kte)
  864. !
  865. !c--- normalized downdraft mass flux profile,also work on bottom detrainment
  866. !--- in this routine
  867. !
  868. call cup_dd_nms(zd,z_cup,cdd,mentrd_rate,jmin,ierr, &
  869. 0,kdet,z1, &
  870. ids,ide, jds,jde, kds,kde, &
  871. ims,ime, jms,jme, kms,kme, &
  872. its,ite, jts,jte, kts,kte)
  873. call cup_dd_nms(zdo,zo_cup,cdd,mentrd_rate,jmin,ierr, &
  874. 1,kdet,z1, &
  875. ids,ide, jds,jde, kds,kde, &
  876. ims,ime, jms,jme, kms,kme, &
  877. its,ite, jts,jte, kts,kte)
  878. !
  879. !--- downdraft moist static energy
  880. !
  881. call cup_dd_he(hes_cup,zd,hcd,z_cup,cdd,mentrd_rate, &
  882. jmin,ierr,he,dbyd,he_cup, &
  883. ids,ide, jds,jde, kds,kde, &
  884. ims,ime, jms,jme, kms,kme, &
  885. its,ite, jts,jte, kts,kte)
  886. call cup_dd_he(heso_cup,zdo,hcdo,zo_cup,cdd,mentrd_rate, &
  887. jmin,ierr,heo,dbydo,he_cup,&
  888. ids,ide, jds,jde, kds,kde, &
  889. ims,ime, jms,jme, kms,kme, &
  890. its,ite, jts,jte, kts,kte)
  891. !
  892. !--- calculate moisture properties of downdraft
  893. !
  894. call cup_dd_moisture(zd,hcd,hes_cup,qcd,qes_cup, &
  895. pwd,q_cup,z_cup,cdd,mentrd_rate,jmin,ierr,gamma_cup, &
  896. pwev,bu,qrcd,q,he,t_cup,2,xl, &
  897. ids,ide, jds,jde, kds,kde, &
  898. ims,ime, jms,jme, kms,kme, &
  899. its,ite, jts,jte, kts,kte)
  900. call cup_dd_moisture(zdo,hcdo,heso_cup,qcdo,qeso_cup, &
  901. pwdo,qo_cup,zo_cup,cdd,mentrd_rate,jmin,ierr,gammao_cup, &
  902. pwevo,bu,qrcdo,qo,heo,tn_cup,1,xl, &
  903. ids,ide, jds,jde, kds,kde, &
  904. ims,ime, jms,jme, kms,kme, &
  905. its,ite, jts,jte, kts,kte)
  906. !
  907. !--- calculate moisture properties of updraft
  908. !
  909. call cup_up_moisture(ierr,z_cup,qc,qrc,pw,pwav, &
  910. kbcon,ktop,cd,dby,mentr_rate,clw_all, &
  911. q,GAMMA_cup,zu,qes_cup,k22,q_cup,xl, &
  912. ids,ide, jds,jde, kds,kde, &
  913. ims,ime, jms,jme, kms,kme, &
  914. its,ite, jts,jte, kts,kte)
  915. do k=kts,ktf
  916. do i=its,itf
  917. cupclw(i,k)=qrc(i,k)
  918. enddo
  919. enddo
  920. call cup_up_moisture(ierr,zo_cup,qco,qrco,pwo,pwavo, &
  921. kbcon,ktop,cd,dbyo,mentr_rate,clw_all, &
  922. qo,GAMMAo_cup,zuo,qeso_cup,k22,qo_cup,xl,&
  923. ids,ide, jds,jde, kds,kde, &
  924. ims,ime, jms,jme, kms,kme, &
  925. its,ite, jts,jte, kts,kte)
  926. !
  927. !--- calculate workfunctions for updrafts
  928. !
  929. call cup_up_aa0(aa0,z,zu,dby,GAMMA_CUP,t_cup, &
  930. kbcon,ktop,ierr, &
  931. ids,ide, jds,jde, kds,kde, &
  932. ims,ime, jms,jme, kms,kme, &
  933. its,ite, jts,jte, kts,kte)
  934. call cup_up_aa0(aa1,zo,zuo,dbyo,GAMMAo_CUP,tn_cup, &
  935. kbcon,ktop,ierr, &
  936. ids,ide, jds,jde, kds,kde, &
  937. ims,ime, jms,jme, kms,kme, &
  938. its,ite, jts,jte, kts,kte)
  939. do i=its,itf
  940. if(ierr(i).eq.0)then
  941. if(aa1(i).eq.0.)then
  942. ierr(i)=17
  943. endif
  944. endif
  945. enddo
  946. !
  947. !--- DETERMINE DOWNDRAFT STRENGTH IN TERMS OF WINDSHEAR
  948. !
  949. call cup_dd_edt(ierr,us,vs,zo,ktop,kbcon,edt,po,pwavo, &
  950. pwevo,edtmax,edtmin,maxens2,edtc, &
  951. ids,ide, jds,jde, kds,kde, &
  952. ims,ime, jms,jme, kms,kme, &
  953. its,ite, jts,jte, kts,kte)
  954. do 250 iedt=1,maxens2
  955. do i=its,itf
  956. if(ierr(i).eq.0)then
  957. edt(i)=edtc(i,iedt)
  958. edto(i)=edtc(i,iedt)
  959. edtx(i)=edtc(i,iedt)
  960. endif
  961. enddo
  962. do k=kts,ktf
  963. do i=its,itf
  964. dellat_ens(i,k,iedt)=0.
  965. dellaq_ens(i,k,iedt)=0.
  966. dellaqc_ens(i,k,iedt)=0.
  967. pwo_ens(i,k,iedt)=0.
  968. enddo
  969. enddo
  970. if (l_flux) then
  971. do k=kts,ktf
  972. do i=its,itf
  973. cfu1_ens(i,k,iedt)=0.
  974. cfd1_ens(i,k,iedt)=0.
  975. dfu1_ens(i,k,iedt)=0.
  976. efu1_ens(i,k,iedt)=0.
  977. dfd1_ens(i,k,iedt)=0.
  978. efd1_ens(i,k,iedt)=0.
  979. enddo
  980. enddo
  981. endif
  982. !
  983. do i=its,itf
  984. aad(i)=0.
  985. enddo
  986. ! do i=its,itf
  987. ! if(ierr(i).eq.0)then
  988. ! eddt(i,j)=edt(i)
  989. ! EDTX(I)=EDT(I)
  990. ! BU(I)=0.
  991. ! BUO(I)=0.
  992. ! endif
  993. ! enddo
  994. !
  995. !---downdraft workfunctions
  996. !
  997. ! call cup_dd_aa0(edt,ierr,aa0,jmin,gamma_cup,t_cup, &
  998. ! hcd,hes_cup,z,zd, &
  999. ! ids,ide, jds,jde, kds,kde, &
  1000. ! ims,ime, jms,jme, kms,kme, &
  1001. ! its,ite, jts,jte, kts,kte)
  1002. ! call cup_dd_aa0(edto,ierr,aad,jmin,gammao_cup,tn_cup, &
  1003. ! hcdo,heso_cup,zo,zdo, &
  1004. ! ids,ide, jds,jde, kds,kde, &
  1005. ! ims,ime, jms,jme, kms,kme, &
  1006. ! its,ite, jts,jte, kts,kte)
  1007. !
  1008. !--- change per unit mass that a model cloud would modify the environment
  1009. !
  1010. !--- 1. in bottom layer
  1011. !
  1012. call cup_dellabot(ipr,jpr,heo_cup,ierr,zo_cup,po,hcdo,edto, &
  1013. zuo,zdo,cdd,heo,dellah,j,mentrd_rate,zo,g, &
  1014. CFU1,CFD1,DFU1,EFU1,DFD1,EFD1,l_flux, &
  1015. ids,ide, jds,jde, kds,kde, &
  1016. ims,ime, jms,jme, kms,kme, &
  1017. its,ite, jts,jte, kts,kte)
  1018. call cup_dellabot(ipr,jpr,qo_cup,ierr,zo_cup,po,qrcdo,edto, &
  1019. zuo,zdo,cdd,qo,dellaq,j,mentrd_rate,zo,g,&
  1020. CFU1,CFD1,DFU1,EFU1,DFD1,EFD1,.FALSE., &
  1021. ids,ide, jds,jde, kds,kde, &
  1022. ims,ime, jms,jme, kms,kme, &
  1023. its,ite, jts,jte, kts,kte)
  1024. !
  1025. !--- 2. everywhere else
  1026. !
  1027. call cup_dellas(ierr,zo_cup,po_cup,hcdo,edto,zdo,cdd, &
  1028. heo,dellah,j,mentrd_rate,zuo,g, &
  1029. cd,hco,ktop,k22,kbcon,mentr_rate,jmin,heo_cup,kdet, &
  1030. k22,ipr,jpr,'deep', &
  1031. CFU1,CFD1,DFU1,EFU1,DFD1,EFD1,l_flux, &
  1032. ids,ide, jds,jde, kds,kde, &
  1033. ims,ime, jms,jme, kms,kme, &
  1034. its,ite, jts,jte, kts,kte)
  1035. !
  1036. !-- take out cloud liquid water for detrainment
  1037. !
  1038. !?? do k=kts,ktf
  1039. do k=kts,ktf-1
  1040. do i=its,itf
  1041. scr1(i,k)=0.
  1042. dellaqc(i,k)=0.
  1043. if(ierr(i).eq.0)then
  1044. ! print *,'in vupnewg, after della ',ierr(i),aa0(i),i,j
  1045. scr1(i,k)=qco(i,k)-qrco(i,k)
  1046. if(k.eq.ktop(i)-0)dellaqc(i,k)= &
  1047. .01*zuo(i,ktop(i))*qrco(i,ktop(i))* &
  1048. 9.81/(po_cup(i,k)-po_cup(i,k+1))
  1049. if(k.lt.ktop(i).and.k.gt.kbcon(i))then
  1050. dz=zo_cup(i,k+1)-zo_cup(i,k)
  1051. dellaqc(i,k)=.01*9.81*cd(i,k)*dz*zuo(i,k) &
  1052. *.5*(qrco(i,k)+qrco(i,k+1))/ &
  1053. (po_cup(i,k)-po_cup(i,k+1))
  1054. endif
  1055. endif
  1056. enddo
  1057. enddo
  1058. call cup_dellas(ierr,zo_cup,po_cup,qrcdo,edto,zdo,cdd, &
  1059. qo,dellaq,j,mentrd_rate,zuo,g, &
  1060. cd,scr1,ktop,k22,kbcon,mentr_rate,jmin,qo_cup,kdet, &
  1061. k22,ipr,jpr,'deep', &
  1062. CFU1,CFD1,DFU1,EFU1,DFD1,EFD1,.FALSE., &
  1063. ids,ide, jds,jde, kds,kde, &
  1064. ims,ime, jms,jme, kms,kme, &
  1065. its,ite, jts,jte, kts,kte )
  1066. !
  1067. !--- using dellas, calculate changed environmental profiles
  1068. !
  1069. ! do 200 nens=1,maxens
  1070. mbdt=mbdt_ens(2)
  1071. do i=its,itf
  1072. xaa0_ens(i,1)=0.
  1073. xaa0_ens(i,2)=0.
  1074. xaa0_ens(i,3)=0.
  1075. enddo
  1076. ! mbdt=mbdt_ens(nens)
  1077. ! do i=its,itf
  1078. ! xaa0_ens(i,nens)=0.
  1079. ! enddo
  1080. do k=kts,ktf
  1081. do i=its,itf
  1082. dellat(i,k)=0.
  1083. if(ierr(i).eq.0)then
  1084. XHE(I,K)=DELLAH(I,K)*MBDT+HEO(I,K)
  1085. XQ(I,K)=DELLAQ(I,K)*MBDT+QO(I,K)
  1086. DELLAT(I,K)=(1./cp)*(DELLAH(I,K)-xl*DELLAQ(I,K))
  1087. XT(I,K)= DELLAT(I,K)*MBDT+TN(I,K)
  1088. IF(XQ(I,K).LE.0.)XQ(I,K)=1.E-08
  1089. ! if(i.eq.ipr.and.j.eq.jpr)then
  1090. ! print *,k,DELLAH(I,K),DELLAQ(I,K),DELLAT(I,K)
  1091. ! endif
  1092. ENDIF
  1093. enddo
  1094. enddo
  1095. do i=its,itf
  1096. if(ierr(i).eq.0)then
  1097. XHE(I,ktf)=HEO(I,ktf)
  1098. XQ(I,ktf)=QO(I,ktf)
  1099. XT(I,ktf)=TN(I,ktf)
  1100. IF(XQ(I,ktf).LE.0.)XQ(I,ktf)=1.E-08
  1101. endif
  1102. enddo
  1103. !
  1104. !--- calculate moist static energy, heights, qes
  1105. !
  1106. call cup_env(xz,xqes,xhe,xhes,xt,xq,po,z1, &
  1107. psur,ierr,tcrit,2,xl,cp, &
  1108. ids,ide, jds,jde, kds,kde, &
  1109. ims,ime, jms,jme, kms,kme, &
  1110. its,ite, jts,jte, kts,kte)
  1111. !
  1112. !--- environmental values on cloud levels
  1113. !
  1114. call cup_env_clev(xt,xqes,xq,xhe,xhes,xz,po,xqes_cup,xq_cup, &
  1115. xhe_cup,xhes_cup,xz_cup,po_cup,gamma_cup,xt_cup,psur, &
  1116. ierr,z1,xl,rv,cp, &
  1117. ids,ide, jds,jde, kds,kde, &
  1118. ims,ime, jms,jme, kms,kme, &
  1119. its,ite, jts,jte, kts,kte)
  1120. !
  1121. !
  1122. !**************************** static control
  1123. !
  1124. !--- moist static energy inside cloud
  1125. !
  1126. do i=its,itf
  1127. if(ierr(i).eq.0)then
  1128. xhkb(i)=xhe(i,k22(i))
  1129. endif
  1130. enddo
  1131. call cup_up_he(k22,xhkb,xz_cup,cd,mentr_rate,xhe_cup,xhc, &
  1132. kbcon,ierr,xdby,xhe,xhes_cup, &
  1133. ids,ide, jds,jde, kds,kde, &
  1134. ims,ime, jms,jme, kms,kme, &
  1135. its,ite, jts,jte, kts,kte)
  1136. !
  1137. !c--- normalized mass flux profile
  1138. !
  1139. call cup_up_nms(xzu,xz_cup,mentr_rate,cd,kbcon,ktop,ierr,k22, &
  1140. ids,ide, jds,jde, kds,kde, &
  1141. ims,ime, jms,jme, kms,kme, &
  1142. its,ite, jts,jte, kts,kte)
  1143. !
  1144. !--- moisture downdraft
  1145. !
  1146. call cup_dd_nms(xzd,xz_cup,cdd,mentrd_rate,jmin,ierr, &
  1147. 1,kdet,z1, &
  1148. ids,ide, jds,jde, kds,kde, &
  1149. ims,ime, jms,jme, kms,kme, &
  1150. its,ite, jts,jte, kts,kte)
  1151. call cup_dd_he(xhes_cup,xzd,xhcd,xz_cup,cdd,mentrd_rate, &
  1152. jmin,ierr,xhe,dbyd,xhe_cup,&
  1153. ids,ide, jds,jde, kds,kde, &
  1154. ims,ime, jms,jme, kms,kme, &
  1155. its,ite, jts,jte, kts,kte)
  1156. call cup_dd_moisture(xzd,xhcd,xhes_cup,xqcd,xqes_cup, &
  1157. xpwd,xq_cup,xz_cup,cdd,mentrd_rate,jmin,ierr,gamma_cup, &
  1158. xpwev,bu,xqrcd,xq,xhe,xt_cup,3,xl, &
  1159. ids,ide, jds,jde, kds,kde, &
  1160. ims,ime, jms,jme, kms,kme, &
  1161. its,ite, jts,jte, kts,kte)
  1162. !
  1163. !------- MOISTURE updraft
  1164. !
  1165. call cup_up_moisture(ierr,xz_cup,xqc,xqrc,xpw,xpwav, &
  1166. kbcon,ktop,cd,xdby,mentr_rate,clw_all, &
  1167. xq,GAMMA_cup,xzu,xqes_cup,k22,xq_cup,xl, &
  1168. ids,ide, jds,jde, kds,kde, &
  1169. ims,ime, jms,jme, kms,kme, &
  1170. its,ite, jts,jte, kts,kte)
  1171. !
  1172. !--- workfunctions for updraft
  1173. !
  1174. call cup_up_aa0(xaa0,xz,xzu,xdby,GAMMA_CUP,xt_cup, &
  1175. kbcon,ktop,ierr, &
  1176. ids,ide, jds,jde, kds,kde, &
  1177. ims,ime, jms,jme, kms,kme, &
  1178. its,ite, jts,jte, kts,kte)
  1179. !
  1180. !--- workfunctions for downdraft
  1181. !
  1182. !
  1183. ! call cup_dd_aa0(edtx,ierr,xaa0,jmin,gamma_cup,xt_cup, &
  1184. ! xhcd,xhes_cup,xz,xzd, &
  1185. ! ids,ide, jds,jde, kds,kde, &
  1186. ! ims,ime, jms,jme, kms,kme, &
  1187. ! its,ite, jts,jte, kts,kte)
  1188. do 200 nens=1,maxens
  1189. do i=its,itf
  1190. if(ierr(i).eq.0)then
  1191. xaa0_ens(i,nens)=xaa0(i)
  1192. nall=(iens-1)*maxens3*maxens*maxens2 &
  1193. +(iedt-1)*maxens*maxens3 &
  1194. +(nens-1)*maxens3
  1195. do k=kts,ktf
  1196. if(k.le.ktop(i))then
  1197. do nens3=1,maxens3
  1198. if(nens3.eq.7)then
  1199. !--- b=0
  1200. pr_ens(i,j,nall+nens3)=pr_ens(i,j,nall+nens3)+ &
  1201. pwo(i,k)
  1202. ! +edto(i)*pwdo(i,k)
  1203. !--- b=beta
  1204. else if(nens3.eq.8)then
  1205. pr_ens(i,j,nall+nens3)=pr_ens(i,j,nall+nens3)+ &
  1206. pwo(i,k)
  1207. !--- b=beta/2
  1208. else if(nens3.eq.9)then
  1209. pr_ens(i,j,nall+nens3)=pr_ens(i,j,nall+nens3)+ &
  1210. pwo(i,k)
  1211. ! +.5*edto(i)*pwdo(i,k)
  1212. else
  1213. pr_ens(i,j,nall+nens3)=pr_ens(i,j,nall+nens3)+ &
  1214. pwo(i,k)+edto(i)*pwdo(i,k)
  1215. endif
  1216. enddo
  1217. endif
  1218. enddo
  1219. if(pr_ens(i,j,nall+7).lt.1.e-6)then
  1220. ierr(i)=18
  1221. do nens3=1,maxens3
  1222. pr_ens(i,j,nall+nens3)=0.
  1223. enddo
  1224. endif
  1225. do nens3=1,maxens3
  1226. if(pr_ens(i,j,nall+nens3).lt.1.e-4)then
  1227. pr_ens(i,j,nall+nens3)=0.
  1228. endif
  1229. enddo
  1230. endif
  1231. enddo
  1232. 200 continue
  1233. !
  1234. !--- LARGE SCALE FORCING
  1235. !
  1236. !
  1237. !------- CHECK wether aa0 should have been zero
  1238. !
  1239. !
  1240. CALL cup_MAXIMI(HEO_CUP,3,KBMAX,K22x,ierr, &
  1241. ids,ide, jds,jde, kds,kde, &
  1242. ims,ime, jms,jme, kms,kme, &
  1243. its,ite, jts,jte, kts,kte)
  1244. do i=its,itf
  1245. ierr2(i)=ierr(i)
  1246. ierr3(i)=ierr(i)
  1247. enddo
  1248. call cup_kbcon(cap_max_increment,2,k22x,kbconx,heo_cup, &
  1249. heso_cup,ierr2,kbmax,po_cup,cap_max, &
  1250. ids,ide, jds,jde, kds,kde, &
  1251. ims,ime, jms,jme, kms,kme, &
  1252. its,ite, jts,jte, kts,kte)
  1253. call cup_kbcon(cap_max_increment,3,k22x,kbconx,heo_cup, &
  1254. heso_cup,ierr3,kbmax,po_cup,cap_max, &
  1255. ids,ide, jds,jde, kds,kde, &
  1256. ims,ime, jms,jme, kms,kme, &
  1257. its,ite, jts,jte, kts,kte)
  1258. !
  1259. !--- DETERMINE THE LEVEL OF CONVECTIVE CLOUD BASE - KBCON
  1260. !
  1261. call cup_forcing_ens(closure_n,xland1,aa0,aa1,xaa0_ens,mbdt_ens,dtime, &
  1262. ierr,ierr2,ierr3,xf_ens,j,'deeps', &
  1263. maxens,iens,iedt,maxens2,maxens3,mconv, &
  1264. po_cup,ktop,omeg,zdo,k22,zuo,pr_ens,edto,kbcon, &
  1265. massflx,iact,direction,ensdim,massfln,ichoice, &
  1266. ids,ide, jds,jde, kds,kde, &
  1267. ims,ime, jms,jme, kms,kme, &
  1268. its,ite, jts,jte, kts,kte)
  1269. !
  1270. do k=kts,ktf
  1271. do i=its,itf
  1272. if(ierr(i).eq.0)then
  1273. dellat_ens(i,k,iedt)=dellat(i,k)
  1274. dellaq_ens(i,k,iedt)=dellaq(i,k)
  1275. dellaqc_ens(i,k,iedt)=dellaqc(i,k)
  1276. pwo_ens(i,k,iedt)=pwo(i,k)+edt(i)*pwdo(i,k)
  1277. else
  1278. dellat_ens(i,k,iedt)=0.
  1279. dellaq_ens(i,k,iedt)=0.
  1280. dellaqc_ens(i,k,iedt)=0.
  1281. pwo_ens(i,k,iedt)=0.
  1282. endif
  1283. ! if(i.eq.ipr.and.j.eq.jpr)then
  1284. ! print *,iens,iedt,dellat(i,k),dellat_ens(i,k,iedt), &
  1285. ! dellaq(i,k), dellaqc(i,k)
  1286. ! endif
  1287. enddo
  1288. enddo
  1289. if (l_flux) then
  1290. do k=kts,ktf
  1291. do i=its,itf
  1292. if(ierr(i).eq.0)then
  1293. cfu1_ens(i,k,iedt)=cfu1(i,k)
  1294. cfd1_ens(i,k,iedt)=cfd1(i,k)
  1295. dfu1_ens(i,k,iedt)=dfu1(i,k)
  1296. efu1_ens(i,k,iedt)=efu1(i,k)
  1297. dfd1_ens(i,k,iedt)=dfd1(i,k)
  1298. efd1_ens(i,k,iedt)=efd1(i,k)
  1299. else
  1300. cfu1_ens(i,k,iedt)=0.
  1301. cfd1_ens(i,k,iedt)=0.
  1302. dfu1_ens(i,k,iedt)=0.
  1303. efu1_ens(i,k,iedt)=0.
  1304. dfd1_ens(i,k,iedt)=0.
  1305. efd1_ens(i,k,iedt)=0.
  1306. end if
  1307. end do
  1308. end do
  1309. end if
  1310. 250 continue
  1311. !
  1312. !--- FEEDBACK
  1313. !
  1314. call cup_output_ens(xf_ens,ierr,dellat_ens,dellaq_ens, &
  1315. dellaqc_ens,outt,outq,outqc,pre,pwo_ens,xmb,ktop, &
  1316. j,'deep',maxens2,maxens,iens,ierr2,ierr3, &
  1317. pr_ens,maxens3,ensdim,massfln, &
  1318. APR_GR,APR_W,APR_MC,APR_ST,APR_AS, &
  1319. APR_CAPMA,APR_CAPME,APR_CAPMI,closure_n,xland1, &
  1320. outCFU1,outCFD1,outDFU1,outEFU1,outDFD1,outEFD1, &
  1321. CFU1_ens,CFD1_ens,DFU1_ens,EFU1_ens,DFD1_ens,EFD1_ens, &
  1322. l_flux, &
  1323. ids,ide, jds,jde, kds,kde, &
  1324. ims,ime, jms,jme, kms,kme, &
  1325. its,ite, jts,jte, kts,kte)
  1326. do i=its,itf
  1327. PRE(I)=MAX(PRE(I),0.)
  1328. enddo
  1329. !
  1330. !---------------------------done------------------------------
  1331. !
  1332. END SUBROUTINE CUP_enss
  1333. SUBROUTINE cup_dd_aa0(edt,ierr,aa0,jmin,gamma_cup,t_cup, &
  1334. hcd,hes_cup,z,zd, &
  1335. ids,ide, jds,jde, kds,kde, &
  1336. ims,ime, jms,jme, kms,kme, &
  1337. its,ite, jts,jte, kts,kte )
  1338. IMPLICIT NONE
  1339. !
  1340. ! on input
  1341. !
  1342. ! only local wrf dimensions are need as of now in this routine
  1343. integer &
  1344. ,intent (in ) :: &
  1345. ids,ide, jds,jde, kds,kde, &
  1346. ims,ime, jms,jme, kms,kme, &
  1347. its,ite, jts,jte, kts,kte
  1348. ! aa0 cloud work function for downdraft
  1349. ! gamma_cup = gamma on model cloud levels
  1350. ! t_cup = temperature (Kelvin) on model cloud levels
  1351. ! hes_cup = saturation moist static energy on model cloud levels
  1352. ! hcd = moist static energy in downdraft
  1353. ! edt = epsilon
  1354. ! zd normalized downdraft mass flux
  1355. ! z = heights of model levels
  1356. ! ierr error value, maybe modified in this routine
  1357. !
  1358. real, dimension (its:ite,kts:kte) &
  1359. ,intent (in ) :: &
  1360. z,zd,gamma_cup,t_cup,hes_cup,hcd
  1361. real, dimension (its:ite) &
  1362. ,intent (in ) :: &
  1363. edt
  1364. integer, dimension (its:ite) &
  1365. ,intent (in ) :: &
  1366. jmin
  1367. !
  1368. ! input and output
  1369. !
  1370. integer, dimension (its:ite) &
  1371. ,intent (inout) :: &
  1372. ierr
  1373. real, dimension (its:ite) &
  1374. ,intent (out ) :: &
  1375. aa0
  1376. !
  1377. ! local variables in this routine
  1378. !
  1379. integer :: &
  1380. i,k,kk
  1381. real :: &
  1382. dz
  1383. !
  1384. integer :: itf, ktf
  1385. !
  1386. itf=MIN(ite,ide-1)
  1387. ktf=MIN(kte,kde-1)
  1388. !
  1389. !?? DO k=kts,kte-1
  1390. DO k=kts,ktf-1
  1391. do i=its,itf
  1392. IF(ierr(I).eq.0.and.k.lt.jmin(i))then
  1393. KK=JMIN(I)-K
  1394. !
  1395. !--- ORIGINAL
  1396. !
  1397. DZ=(Z(I,KK)-Z(I,KK+1))
  1398. AA0(I)=AA0(I)+zd(i,kk)*EDT(I)*DZ*(9.81/(1004.*T_cup(I,KK))) &
  1399. *((hcd(i,kk)-hes_cup(i,kk))/(1.+GAMMA_cup(i,kk)))
  1400. endif
  1401. enddo
  1402. enddo
  1403. END SUBROUTINE CUP_dd_aa0
  1404. SUBROUTINE cup_dd_edt(ierr,us,vs,z,ktop,kbcon,edt,p,pwav, &
  1405. pwev,edtmax,edtmin,maxens2,edtc, &
  1406. ids,ide, jds,jde, kds,kde, &
  1407. ims,ime, jms,jme, kms,kme, &
  1408. its,ite, jts,jte, kts,kte )
  1409. IMPLICIT NONE
  1410. integer &
  1411. ,intent (in ) :: &
  1412. ids,ide, jds,jde, kds,kde, &
  1413. ims,ime, jms,jme, kms,kme, &
  1414. its,ite, jts,jte, kts,kte
  1415. integer, intent (in ) :: &
  1416. maxens2
  1417. !
  1418. ! ierr error value, maybe modified in this routine
  1419. !
  1420. real, dimension (its:ite,kts:kte) &
  1421. ,intent (in ) :: &
  1422. us,vs,z,p
  1423. real, dimension (its:ite,1:maxens2) &
  1424. ,intent (out ) :: &
  1425. edtc
  1426. real, dimension (its:ite) &
  1427. ,intent (out ) :: &
  1428. edt
  1429. real, dimension (its:ite) &
  1430. ,intent (in ) :: &
  1431. pwav,pwev
  1432. real &
  1433. ,intent (in ) :: &
  1434. edtmax,edtmin
  1435. integer, dimension (its:ite) &
  1436. ,intent (in ) :: &
  1437. ktop,kbcon
  1438. integer, dimension (its:ite) &
  1439. ,intent (inout) :: &
  1440. ierr
  1441. !
  1442. ! local variables in this routine
  1443. !
  1444. integer i,k,kk
  1445. real einc,pef,pefb,prezk,zkbc
  1446. real, dimension (its:ite) :: &
  1447. vshear,sdp,vws
  1448. integer :: itf, ktf
  1449. itf=MIN(ite,ide-1)
  1450. ktf=MIN(kte,kde-1)
  1451. !
  1452. !--- DETERMINE DOWNDRAFT STRENGTH IN TERMS OF WINDSHEAR
  1453. !
  1454. ! */ calculate an average wind shear over the depth of the cloud
  1455. !
  1456. do i=its,itf
  1457. edt(i)=0.
  1458. vws(i)=0.
  1459. sdp(i)=0.
  1460. vshear(i)=0.
  1461. enddo
  1462. do kk = kts,ktf-1
  1463. do 62 i=its,itf
  1464. IF(ierr(i).ne.0)GO TO 62
  1465. if (kk .le. min0(ktop(i),ktf) .and. kk .ge. kbcon(i)) then
  1466. vws(i) = vws(i)+ &
  1467. (abs((us(i,kk+1)-us(i,kk))/(z(i,kk+1)-z(i,kk))) &
  1468. + abs((vs(i,kk+1)-vs(i,kk))/(z(i,kk+1)-z(i,kk)))) * &
  1469. (p(i,kk) - p(i,kk+1))
  1470. sdp(i) = sdp(i) + p(i,kk) - p(i,kk+1)
  1471. endif
  1472. if (kk .eq. ktf-1)vshear(i) = 1.e3 * vws(i) / sdp(i)
  1473. 62 continue
  1474. end do
  1475. do i=its,itf
  1476. IF(ierr(i).eq.0)then
  1477. pef=(1.591-.639*VSHEAR(I)+.0953*(VSHEAR(I)**2) &
  1478. -.00496*(VSHEAR(I)**3))
  1479. if(pef.gt.edtmax)pef=edtmax
  1480. if(pef.lt.edtmin)pef=edtmin
  1481. !
  1482. !--- cloud base precip efficiency
  1483. !
  1484. zkbc=z(i,kbcon(i))*3.281e-3
  1485. prezk=.02
  1486. if(zkbc.gt.3.)then
  1487. prezk=.96729352+zkbc*(-.70034167+zkbc*(.162179896+zkbc &
  1488. *(- 1.2569798E-2+zkbc*(4.2772E-4-zkbc*5.44E-6))))
  1489. endif
  1490. if(zkbc.gt.25)then
  1491. prezk=2.4
  1492. endif
  1493. pefb=1./(1.+prezk)
  1494. if(pefb.gt.edtmax)pefb=edtmax
  1495. if(pefb.lt.edtmin)pefb=edtmin
  1496. EDT(I)=1.-.5*(pefb+pef)
  1497. !--- edt here is 1-precipeff!
  1498. ! einc=(1.-edt(i))/float(maxens2)
  1499. ! einc=edt(i)/float(maxens2+1)
  1500. !--- 20 percent
  1501. einc=.2*edt(i)
  1502. do k=1,maxens2
  1503. edtc(i,k)=edt(i)+float(k-2)*einc
  1504. enddo
  1505. endif
  1506. enddo
  1507. do i=its,itf
  1508. IF(ierr(i).eq.0)then
  1509. do k=1,maxens2
  1510. EDTC(I,K)=-EDTC(I,K)*PWAV(I)/PWEV(I)
  1511. IF(EDTC(I,K).GT.edtmax)EDTC(I,K)=edtmax
  1512. IF(EDTC(I,K).LT.edtmin)EDTC(I,K)=edtmin
  1513. enddo
  1514. endif
  1515. enddo
  1516. END SUBROUTINE cup_dd_edt
  1517. SUBROUTINE cup_dd_he(hes_cup,zd,hcd,z_cup,cdd,entr, &
  1518. jmin,ierr,he,dby,he_cup, &
  1519. ids,ide, jds,jde, kds,kde, &
  1520. ims,ime, jms,jme, kms,kme, &
  1521. its,ite, jts,jte, kts,kte )
  1522. IMPLICIT NONE
  1523. !
  1524. ! on input
  1525. !
  1526. ! only local wrf dimensions are need as of now in this routine
  1527. integer &
  1528. ,intent (in ) :: &
  1529. ids,ide, jds,jde, kds,kde, &
  1530. ims,ime, jms,jme, kms,kme, &
  1531. its,ite, jts,jte, kts,kte
  1532. ! hcd = downdraft moist static energy
  1533. ! he = moist static energy on model levels
  1534. ! he_cup = moist static energy on model cloud levels
  1535. ! hes_cup = saturation moist static energy on model cloud levels
  1536. ! dby = buoancy term
  1537. ! cdd= detrainment function
  1538. ! z_cup = heights of model cloud levels
  1539. ! entr = entrainment rate
  1540. ! zd = downdraft normalized mass flux
  1541. !
  1542. real, dimension (its:ite,kts:kte) &
  1543. ,intent (in ) :: &
  1544. he,he_cup,hes_cup,z_cup,cdd,zd
  1545. ! entr= entrainment rate
  1546. real &
  1547. ,intent (in ) :: &
  1548. entr
  1549. integer, dimension (its:ite) &
  1550. ,intent (in ) :: &
  1551. jmin
  1552. !
  1553. ! input and output
  1554. !
  1555. ! ierr error value, maybe modified in this routine
  1556. integer, dimension (its:ite) &
  1557. ,intent (inout) :: &
  1558. ierr
  1559. real, dimension (its:ite,kts:kte) &
  1560. ,intent (out ) :: &
  1561. hcd,dby
  1562. !
  1563. ! local variables in this routine
  1564. !
  1565. integer :: &
  1566. i,k,ki
  1567. real :: &
  1568. dz
  1569. integer :: itf, ktf
  1570. itf=MIN(ite,ide-1)
  1571. ktf=MIN(kte,kde-1)
  1572. do k=kts+1,ktf
  1573. do i=its,itf
  1574. dby(i,k)=0.
  1575. IF(ierr(I).eq.0)then
  1576. hcd(i,k)=hes_cup(i,k)
  1577. endif
  1578. enddo
  1579. enddo
  1580. !
  1581. do 100 i=its,itf
  1582. IF(ierr(I).eq.0)then
  1583. k=jmin(i)
  1584. hcd(i,k)=hes_cup(i,k)
  1585. dby(i,k)=hcd(i,jmin(i))-hes_cup(i,k)
  1586. !
  1587. do ki=jmin(i)-1,1,-1
  1588. DZ=Z_cup(i,Ki+1)-Z_cup(i,Ki)
  1589. HCD(i,Ki)=(HCD(i,Ki+1)*(1.-.5*CDD(i,Ki)*DZ) &
  1590. +entr*DZ*HE(i,Ki) &
  1591. )/(1.+entr*DZ-.5*CDD(i,Ki)*DZ)
  1592. dby(i,ki)=HCD(i,Ki)-hes_cup(i,ki)
  1593. enddo
  1594. !
  1595. endif
  1596. !--- end loop over i
  1597. 100 continue
  1598. END SUBROUTINE cup_dd_he
  1599. SUBROUTINE cup_dd_moisture(zd,hcd,hes_cup,qcd,qes_cup, &
  1600. pwd,q_cup,z_cup,cdd,entr,jmin,ierr, &
  1601. gamma_cup,pwev,bu,qrcd, &
  1602. q,he,t_cup,iloop,xl, &
  1603. ids,ide, jds,jde, kds,kde, &
  1604. ims,ime, jms,jme, kms,kme, &
  1605. its,ite, jts,jte, kts,kte )
  1606. IMPLICIT NONE
  1607. integer &
  1608. ,intent (in ) :: &
  1609. ids,ide, jds,jde, kds,kde, &
  1610. ims,ime, jms,jme, kms,kme, &
  1611. its,ite, jts,jte, kts,kte
  1612. ! cdd= detrainment function
  1613. ! q = environmental q on model levels
  1614. ! q_cup = environmental q on model cloud levels
  1615. ! qes_cup = saturation q on model cloud levels
  1616. ! hes_cup = saturation h on model cloud levels
  1617. ! hcd = h in model cloud
  1618. ! bu = buoancy term
  1619. ! zd = normalized downdraft mass flux
  1620. ! gamma_cup = gamma on model cloud levels
  1621. ! mentr_rate = entrainment rate
  1622. ! qcd = cloud q (including liquid water) after entrainment
  1623. ! qrch = saturation q in cloud
  1624. ! pwd = evaporate at that level
  1625. ! pwev = total normalized integrated evaoprate (I2)
  1626. ! entr= entrainment rate
  1627. !
  1628. real, dimension (its:ite,kts:kte) &
  1629. ,intent (in ) :: &
  1630. zd,t_cup,hes_cup,hcd,qes_cup,q_cup,z_cup,cdd,gamma_cup,q,he
  1631. real &
  1632. ,intent (in ) :: &
  1633. entr,xl
  1634. integer &
  1635. ,intent (in ) :: &
  1636. iloop
  1637. integer, dimension (its:ite) &
  1638. ,intent (in ) :: &
  1639. jmin
  1640. integer, dimension (its:ite) &
  1641. ,intent (inout) :: &
  1642. ierr
  1643. real, dimension (its:ite,kts:kte) &
  1644. ,intent (out ) :: &
  1645. qcd,qrcd,pwd
  1646. real, dimension (its:ite) &
  1647. ,intent (out ) :: &
  1648. pwev,bu
  1649. !
  1650. ! local variables in this routine
  1651. !
  1652. integer :: &
  1653. i,k,ki
  1654. real :: &
  1655. dh,dz,dqeva
  1656. integer :: itf, ktf
  1657. itf=MIN(ite,ide-1)
  1658. ktf=MIN(kte,kde-1)
  1659. do i=its,itf
  1660. bu(i)=0.
  1661. pwev(i)=0.
  1662. enddo
  1663. do k=kts,ktf
  1664. do i=its,itf
  1665. qcd(i,k)=0.
  1666. qrcd(i,k)=0.
  1667. pwd(i,k)=0.
  1668. enddo
  1669. enddo
  1670. !
  1671. !
  1672. !
  1673. do 100 i=its,itf
  1674. IF(ierr(I).eq.0)then
  1675. k=jmin(i)
  1676. DZ=Z_cup(i,K+1)-Z_cup(i,K)
  1677. qcd(i,k)=q_cup(i,k)
  1678. ! qcd(i,k)=.5*(qes_cup(i,k)+q_cup(i,k))
  1679. qrcd(i,k)=qes_cup(i,k)
  1680. pwd(i,jmin(i))=min(0.,qcd(i,k)-qrcd(i,k))
  1681. pwev(i)=pwev(i)+pwd(i,jmin(i))
  1682. qcd(i,k)=qes_cup(i,k)
  1683. !
  1684. DH=HCD(I,k)-HES_cup(I,K)
  1685. bu(i)=dz*dh
  1686. do ki=jmin(i)-1,1,-1
  1687. DZ=Z_cup(i,Ki+1)-Z_cup(i,Ki)
  1688. QCD(i,Ki)=(qCD(i,Ki+1)*(1.-.5*CDD(i,Ki)*DZ) &
  1689. +entr*DZ*q(i,Ki) &
  1690. )/(1.+entr*DZ-.5*CDD(i,Ki)*DZ)
  1691. !
  1692. !--- to be negatively buoyant, hcd should be smaller than hes!
  1693. !
  1694. DH=HCD(I,ki)-HES_cup(I,Ki)
  1695. bu(i)=bu(i)+dz*dh
  1696. QRCD(I,Ki)=qes_cup(i,ki)+(1./XL)*(GAMMA_cup(i,ki) &
  1697. /(1.+GAMMA_cup(i,ki)))*DH
  1698. dqeva=qcd(i,ki)-qrcd(i,ki)
  1699. if(dqeva.gt.0.)dqeva=0.
  1700. pwd(i,ki)=zd(i,ki)*dqeva
  1701. qcd(i,ki)=qrcd(i,ki)
  1702. pwev(i)=pwev(i)+pwd(i,ki)
  1703. ! if(iloop.eq.1.and.i.eq.102.and.j.eq.62)then
  1704. ! print *,'in cup_dd_moi ', hcd(i,ki),HES_cup(I,Ki),dh,dqeva
  1705. ! endif
  1706. enddo
  1707. !
  1708. !--- end loop over i
  1709. if(pwev(I).eq.0.and.iloop.eq.1)then
  1710. ! print *,'problem with buoy in cup_dd_moisture',i
  1711. ierr(i)=7
  1712. endif
  1713. if(BU(I).GE.0.and.iloop.eq.1)then
  1714. ! print *,'problem with buoy in cup_dd_moisture',i
  1715. ierr(i)=7
  1716. endif
  1717. endif
  1718. 100 continue
  1719. END SUBROUTINE cup_dd_moisture
  1720. SUBROUTINE cup_dd_nms(zd,z_cup,cdd,entr,jmin,ierr, &
  1721. itest,kdet,z1, &
  1722. ids,ide, jds,jde, kds,kde, &
  1723. ims,ime, jms,jme, kms,kme, &
  1724. its,ite, jts,jte, kts,kte )
  1725. IMPLICIT NONE
  1726. !
  1727. ! on input
  1728. !
  1729. ! only local wrf dimensions are need as of now in this routine
  1730. integer &
  1731. ,intent (in ) :: &
  1732. ids,ide, jds,jde, kds,kde, &
  1733. ims,ime, jms,jme, kms,kme, &
  1734. its,ite, jts,jte, kts,kte
  1735. ! z_cup = height of cloud model level
  1736. ! z1 = terrain elevation
  1737. ! entr = downdraft entrainment rate
  1738. ! jmin = downdraft originating level
  1739. ! kdet = level above ground where downdraft start detraining
  1740. ! itest = flag to whether to calculate cdd
  1741. real, dimension (its:ite,kts:kte) &
  1742. ,intent (in ) :: &
  1743. z_cup
  1744. real, dimension (its:ite) &
  1745. ,intent (in ) :: &
  1746. z1
  1747. real &
  1748. ,intent (in ) :: &
  1749. entr
  1750. integer, dimension (its:ite) &
  1751. ,intent (in ) :: &
  1752. jmin,kdet
  1753. integer &
  1754. ,intent (in ) :: &
  1755. itest
  1756. !
  1757. ! input and output
  1758. !
  1759. ! ierr error value, maybe modified in this routine
  1760. integer, dimension (its:ite) &
  1761. ,intent (inout) :: &
  1762. ierr
  1763. ! zd is the normalized downdraft mass flux
  1764. ! cdd is the downdraft detrainmen function
  1765. real, dimension (its:ite,kts:kte) &
  1766. ,intent (out ) :: &
  1767. zd
  1768. real, dimension (its:ite,kts:kte) &
  1769. ,intent (inout) :: &
  1770. cdd
  1771. !
  1772. ! local variables in this routine
  1773. !
  1774. integer :: &
  1775. i,k,ki
  1776. real :: &
  1777. a,perc,dz
  1778. integer :: itf, ktf
  1779. itf=MIN(ite,ide-1)
  1780. ktf=MIN(kte,kde-1)
  1781. !
  1782. !--- perc is the percentage of mass left when hitting the ground
  1783. !
  1784. perc=.03
  1785. do k=kts,ktf
  1786. do i=its,itf
  1787. zd(i,k)=0.
  1788. if(itest.eq.0)cdd(i,k)=0.
  1789. enddo
  1790. enddo
  1791. a=1.-perc
  1792. !
  1793. !
  1794. !
  1795. do 100 i=its,itf
  1796. IF(ierr(I).eq.0)then
  1797. zd(i,jmin(i))=1.
  1798. !
  1799. !--- integrate downward, specify detrainment(cdd)!
  1800. !
  1801. do ki=jmin(i)-1,1,-1
  1802. DZ=Z_cup(i,Ki+1)-Z_cup(i,Ki)
  1803. if(ki.le.kdet(i).and.itest.eq.0)then
  1804. cdd(i,ki)=entr+(1.- (a*(z_cup(i,ki)-z1(i)) &
  1805. +perc*(z_cup(i,kdet(i))-z1(i)) ) &
  1806. /(a*(z_cup(i,ki+1)-z1(i)) &
  1807. +perc*(z_cup(i,kdet(i))-z1(i))))/dz
  1808. endif
  1809. zd(i,ki)=zd(i,ki+1)*(1.+(entr-cdd(i,ki))*dz)
  1810. enddo
  1811. !
  1812. endif
  1813. !--- end loop over i
  1814. 100 continue
  1815. END SUBROUTINE cup_dd_nms
  1816. SUBROUTINE cup_dellabot(ipr,jpr,he_cup,ierr,z_cup,p_cup, &
  1817. hcd,edt,zu,zd,cdd,he,della,j,mentrd_rate,z,g, &
  1818. CFU1,CFD1,DFU1,EFU1,DFD1,EFD1,l_flux, &
  1819. ids,ide, jds,jde, kds,kde, &
  1820. ims,ime, jms,jme, kms,kme, &
  1821. its,ite, jts,jte, kts,kte )
  1822. IMPLICIT NONE
  1823. integer &
  1824. ,intent (in ) :: &
  1825. ids,ide, jds,jde, kds,kde, &
  1826. ims,ime, jms,jme, kms,kme, &
  1827. its,ite, jts,jte, kts,kte
  1828. integer, intent (in ) :: &
  1829. j,ipr,jpr
  1830. !
  1831. ! ierr error value, maybe modified in this routine
  1832. !
  1833. real, dimension (its:ite,kts:kte) &
  1834. ,intent (out ) :: &
  1835. della
  1836. real, dimension (its:ite,kts:kte) &
  1837. ,intent (in ) :: &
  1838. z_cup,p_cup,hcd,zu,zd,cdd,he,z,he_cup
  1839. real, dimension (its:ite) &
  1840. ,intent (in ) :: &
  1841. edt
  1842. real &
  1843. ,intent (in ) :: &
  1844. g,mentrd_rate
  1845. integer, dimension (its:ite) &
  1846. ,intent (inout) :: &
  1847. ierr
  1848. real, dimension (its:ite,kts:kte) &
  1849. ,intent (inout ) :: &
  1850. CFU1,CFD1,DFU1,EFU1,DFD1,EFD1
  1851. logical, intent(in) :: l_flux
  1852. !
  1853. ! local variables in this routine
  1854. !
  1855. integer i
  1856. real detdo,detdo1,detdo2,entdo,dp,dz,subin, &
  1857. totmas
  1858. !
  1859. integer :: itf, ktf
  1860. itf=MIN(ite,ide-1)
  1861. ktf=MIN(kte,kde-1)
  1862. !
  1863. !
  1864. ! if(j.eq.jpr)print *,'in cup dellabot '
  1865. do 100 i=its,itf
  1866. if (l_flux) then
  1867. cfu1(i,1)=0.
  1868. cfd1(i,1)=0.
  1869. cfu1(i,2)=0.
  1870. cfd1(i,2)=0.
  1871. dfu1(i,1)=0.
  1872. efu1(i,1)=0.
  1873. dfd1(i,1)=0.
  1874. efd1(i,1)=0.
  1875. endif
  1876. della(i,1)=0.
  1877. if(ierr(i).ne.0)go to 100
  1878. dz=z_cup(i,2)-z_cup(i,1)
  1879. DP=100.*(p_cup(i,1)-P_cup(i,2))
  1880. detdo1=edt(i)*zd(i,2)*CDD(i,1)*DZ
  1881. detdo2=edt(i)*zd(i,1)
  1882. entdo=edt(i)*zd(i,2)*mentrd_rate*dz
  1883. subin=-EDT(I)*zd(i,2)
  1884. detdo=detdo1+detdo2-entdo+subin
  1885. DELLA(I,1)=(detdo1*.5*(HCD(i,1)+HCD(i,2)) &
  1886. +detdo2*hcd(i,1) &
  1887. +subin*he_cup(i,2) &
  1888. -entdo*he(i,1))*g/dp
  1889. if (l_flux) then
  1890. cfd1(i,2) = -edt(i)*zd(i,2) !only contribution to subin, subdown=0
  1891. dfd1(i,1) = detdo1+detdo2
  1892. efd1(i,1) = -entdo
  1893. endif
  1894. 100 CONTINUE
  1895. END SUBROUTINE cup_dellabot
  1896. SUBROUTINE cup_dellas(ierr,z_cup,p_cup,hcd,edt,zd,cdd, &
  1897. he,della,j,mentrd_rate,zu,g, &
  1898. cd,hc,ktop,k22,kbcon,mentr_rate,jmin,he_cup,kdet,kpbl, &
  1899. ipr,jpr,name, &
  1900. CFU1,CFD1,DFU1,EFU1,DFD1,EFD1,l_flux, &
  1901. ids,ide, jds,jde, kds,kde, &
  1902. ims,ime, jms,jme, kms,kme, &
  1903. its,ite, jts,jte, kts,kte )
  1904. IMPLICIT NONE
  1905. integer &
  1906. ,intent (in ) :: &
  1907. ids,ide, jds,jde, kds,kde, &
  1908. ims,ime, jms,jme, kms,kme, &
  1909. its,ite, jts,jte, kts,kte
  1910. integer, intent (in ) :: &
  1911. j,ipr,jpr
  1912. !
  1913. ! ierr error value, maybe modified in this routine
  1914. !
  1915. real, dimension (its:ite,kts:kte) &
  1916. ,intent (out ) :: &
  1917. della
  1918. real, dimension (its:ite,kts:kte) &
  1919. ,intent (in ) :: &
  1920. z_cup,p_cup,hcd,zd,cdd,he,hc,cd,zu,he_cup
  1921. real, dimension (its:ite) &
  1922. ,intent (in ) :: &
  1923. edt
  1924. real &
  1925. ,intent (in ) :: &
  1926. g,mentrd_rate,mentr_rate
  1927. integer, dimension (its:ite) &
  1928. ,intent (in ) :: &
  1929. kbcon,ktop,k22,jmin,kdet,kpbl
  1930. integer, dimension (its:ite) &
  1931. ,intent (inout) :: &
  1932. ierr
  1933. character *(*), intent (in) :: &
  1934. name
  1935. real, dimension (its:ite,kts:kte) &
  1936. ,intent (inout ) :: &
  1937. CFU1,CFD1,DFU1,EFU1,DFD1,EFD1
  1938. logical, intent(in) :: l_flux
  1939. !
  1940. ! local variables in this routine
  1941. !
  1942. integer i,k
  1943. real detdo1,detdo2,entdo,dp,dz,subin,detdo,entup, &
  1944. detup,subdown,entdoj,entupk,detupk,totmas
  1945. !
  1946. integer :: itf, ktf
  1947. itf=MIN(ite,ide-1)
  1948. ktf=MIN(kte,kde-1)
  1949. !
  1950. !
  1951. i=ipr
  1952. ! if(j.eq.jpr)then
  1953. ! print *,'in dellas kpbl(i),k22(i),kbcon(i),ktop(i),jmin(i)'
  1954. ! print *,kpbl(i),k22(i),kbcon(i),ktop(i),jmin(i)
  1955. ! endif
  1956. DO K=kts+1,ktf
  1957. do i=its,itf
  1958. della(i,k)=0.
  1959. enddo
  1960. enddo
  1961. if (l_flux) then
  1962. DO K=kts+1,ktf-1
  1963. do i=its,itf
  1964. cfu1(i,k+1)=0.
  1965. cfd1(i,k+1)=0.
  1966. enddo
  1967. enddo
  1968. DO K=kts+1,ktf
  1969. do i=its,itf
  1970. dfu1(i,k)=0.
  1971. efu1(i,k)=0.
  1972. dfd1(i,k)=0.
  1973. efd1(i,k)=0.
  1974. enddo
  1975. enddo
  1976. endif
  1977. !
  1978. DO 100 k=kts+1,ktf-1
  1979. DO 100 i=its,ite
  1980. IF(ierr(i).ne.0)GO TO 100
  1981. IF(K.Gt.KTOP(I))GO TO 100
  1982. !
  1983. !--- SPECIFY DETRAINMENT OF DOWNDRAFT, HAS TO BE CONSISTENT
  1984. !--- WITH ZD CALCULATIONS IN SOUNDD.
  1985. !
  1986. DZ=Z_cup(I,K+1)-Z_cup(I,K)
  1987. detdo=edt(i)*CDD(i,K)*DZ*ZD(i,k+1)
  1988. entdo=edt(i)*mentrd_rate*dz*zd(i,k+1)
  1989. subin=zu(i,k+1)-zd(i,k+1)*edt(i)
  1990. entup=0.
  1991. detup=0.
  1992. if(k.ge.kbcon(i).and.k.lt.ktop(i))then
  1993. entup=mentr_rate*dz*zu(i,k)
  1994. detup=CD(i,K+1)*DZ*ZU(i,k)
  1995. endif
  1996. subdown=(zu(i,k)-zd(i,k)*edt(i))
  1997. entdoj=0.
  1998. entupk=0.
  1999. detupk=0.
  2000. !
  2001. if(k.eq.jmin(i))then
  2002. entdoj=edt(i)*zd(i,k)
  2003. endif
  2004. if(k.eq.k22(i)-1)then
  2005. entupk=zu(i,kpbl(i))
  2006. endif
  2007. if(k.gt.kdet(i))then
  2008. detdo=0.
  2009. endif
  2010. if(k.eq.ktop(i)-0)then
  2011. detupk=zu(i,ktop(i))
  2012. subin=0.
  2013. endif
  2014. if(k.lt.kbcon(i))then
  2015. detup=0.
  2016. endif
  2017. if (l_flux) then
  2018. ! z_cup(k+1): zu(k+1), -zd(k+1) ==> subin ==> cf[du]1 (k+1) (full-eta level k+1)
  2019. !
  2020. ! z(k) : detup, detdo, entup, entdo ==> [de]f[du]1 (k) (half-eta level k )
  2021. !
  2022. ! z_cup(k) : zu(k), -zd(k) ==> subdown ==> cf[du]1 (k) (full-eta level k )
  2023. ! Store downdraft/updraft mass fluxes at full eta level k (z_cup(k)) in cf[ud]1(k):
  2024. cfu1(i,k+1) = zu(i,k+1)
  2025. cfd1(i,k+1) = -edt(i)*zd(i,k+1)
  2026. ! Store detrainment/entrainment mass fluxes at half eta level k (z(k)) in [de]f[du]1(k):
  2027. dfu1(i,k) = detup+detupk
  2028. efu1(i,k) = -(entup+entupk)
  2029. dfd1(i,k) = detdo
  2030. efd1(i,k) = -(entdo+entdoj)
  2031. endif
  2032. !C
  2033. !C--- CHANGED DUE TO SUBSIDENCE AND ENTRAINMENT
  2034. !C
  2035. totmas=subin-subdown+detup-entup-entdo+ &
  2036. detdo-entupk-entdoj+detupk
  2037. ! if(j.eq.jpr.and.i.eq.ipr)print *,'k,totmas,sui,sud = ',k,
  2038. ! 1 totmas,subin,subdown
  2039. ! if(j.eq.jpr.and.i.eq.ipr)print *,'updr stuff = ',detup,
  2040. ! 1 entup,entupk,detupk
  2041. ! if(j.eq.jpr.and.i.eq.ipr)print *,'dddr stuff = ',entdo,
  2042. ! 1 detdo,entdoj
  2043. if(abs(totmas).gt.1.e-6)then
  2044. ! print *,'*********************',i,j,k,totmas,name
  2045. ! print *,kpbl(i),k22(i),kbcon(i),ktop(i)
  2046. !c print *,'updr stuff = ',subin,
  2047. !c 1 subdown,detup,entup,entupk,detupk
  2048. !c print *,'dddr stuff = ',entdo,
  2049. !c 1 detdo,entdoj
  2050. ! call wrf_error_fatal ( 'totmas .gt.1.e-6' )
  2051. endif
  2052. dp=100.*(p_cup(i,k-1)-p_cup(i,k))
  2053. della(i,k)=(subin*he_cup(i,k+1) &
  2054. -subdown*he_cup(i,k) &
  2055. +detup*.5*(HC(i,K+1)+HC(i,K)) &
  2056. +detdo*.5*(HCD(i,K+1)+HCD(i,K)) &
  2057. -entup*he(i,k) &
  2058. -entdo*he(i,k) &
  2059. -entupk*he_cup(i,k22(i)) &
  2060. -entdoj*he_cup(i,jmin(i)) &
  2061. +detupk*hc(i,ktop(i)) &
  2062. )*g/dp
  2063. ! if(i.eq.ipr.and.j.eq.jpr)then
  2064. ! print *,k,della(i,k),subin*he_cup(i,k+1),subdown*he_cup(i,k),
  2065. ! 1 detdo*.5*(HCD(i,K+1)+HCD(i,K))
  2066. ! print *,k,detup*.5*(HC(i,K+1)+HC(i,K)),detupk*hc(i,ktop(i)),
  2067. ! 1 entup*he(i,k),entdo*he(i,k)
  2068. ! print *,k,he_cup(i,k+1),he_cup(i,k),entupk*he_cup(i,k)
  2069. ! endif
  2070. 100 CONTINUE
  2071. END SUBROUTINE cup_dellas
  2072. SUBROUTINE cup_direction2(i,j,dir,id,massflx, &
  2073. iresult,imass,massfld, &
  2074. ids,ide, jds,jde, kds,kde, &
  2075. ims,ime, jms,jme, kms,kme, &
  2076. its,ite, jts,jte, kts,kte )
  2077. IMPLICIT NONE
  2078. integer &
  2079. ,intent (in ) :: &
  2080. ids,ide, jds,jde, kds,kde, &
  2081. ims,ime, jms,jme, kms,kme, &
  2082. its,ite, jts,jte, kts,kte
  2083. integer, intent (in ) :: &
  2084. i,j,imass
  2085. integer, intent (out ) :: &
  2086. iresult
  2087. !
  2088. ! ierr error value, maybe modified in this routine
  2089. !
  2090. integer, dimension (ims:ime,jms:jme) &
  2091. ,intent (in ) :: &
  2092. id
  2093. real, dimension (ims:ime,jms:jme) &
  2094. ,intent (in ) :: &
  2095. massflx
  2096. real, dimension (its:ite) &
  2097. ,intent (inout) :: &
  2098. dir
  2099. real &
  2100. ,intent (out ) :: &
  2101. massfld
  2102. !
  2103. ! local variables in this routine
  2104. !
  2105. integer k,ia,ja,ib,jb
  2106. real diff
  2107. !
  2108. !
  2109. !
  2110. if(imass.eq.1)then
  2111. massfld=massflx(i,j)
  2112. endif
  2113. iresult=0
  2114. ! return
  2115. diff=22.5
  2116. if(dir(i).lt.22.5)dir(i)=360.+dir(i)
  2117. if(id(i,j).eq.1)iresult=1
  2118. ! ja=max(2,j-1)
  2119. ! ia=max(2,i-1)
  2120. ! jb=min(mjx-1,j+1)
  2121. ! ib=min(mix-1,i+1)
  2122. ja=j-1
  2123. ia=i-1
  2124. jb=j+1
  2125. ib=i+1
  2126. if(dir(i).gt.90.-diff.and.dir(i).le.90.+diff)then
  2127. !--- steering flow from the east
  2128. if(id(ib,j).eq.1)then
  2129. iresult=1
  2130. if(imass.eq.1)then
  2131. massfld=max(massflx(ib,j),massflx(i,j))
  2132. endif
  2133. return
  2134. endif
  2135. else if(dir(i).gt.135.-diff.and.dir(i).le.135.+diff)then
  2136. !--- steering flow from the south-east
  2137. if(id(ib,ja).eq.1)then
  2138. iresult=1
  2139. if(imass.eq.1)then
  2140. massfld=max(massflx(ib,ja),massflx(i,j))
  2141. endif
  2142. return
  2143. endif
  2144. !--- steering flow from the south
  2145. else if(dir(i).gt.180.-diff.and.dir(i).le.180.+diff)then
  2146. if(id(i,ja).eq.1)then
  2147. iresult=1
  2148. if(imass.eq.1)then
  2149. massfld=max(massflx(i,ja),massflx(i,j))
  2150. endif
  2151. return
  2152. endif
  2153. !--- steering flow from the south west
  2154. else if(dir(i).gt.225.-diff.and.dir(i).le.225.+diff)then
  2155. if(id(ia,ja).eq.1)then
  2156. iresult=1
  2157. if(imass.eq.1)then
  2158. massfld=max(massflx(ia,ja),massflx(i,j))
  2159. endif
  2160. return
  2161. endif
  2162. !--- steering flow from the west
  2163. else if(dir(i).gt.270.-diff.and.dir(i).le.270.+diff)then
  2164. if(id(ia,j).eq.1)then
  2165. iresult=1
  2166. if(imass.eq.1)then
  2167. massfld=max(massflx(ia,j),massflx(i,j))
  2168. endif
  2169. return
  2170. endif
  2171. !--- steering flow from the north-west
  2172. else if(dir(i).gt.305.-diff.and.dir(i).le.305.+diff)then
  2173. if(id(ia,jb).eq.1)then
  2174. iresult=1
  2175. if(imass.eq.1)then
  2176. massfld=max(massflx(ia,jb),massflx(i,j))
  2177. endif
  2178. return
  2179. endif
  2180. !--- steering flow from the north
  2181. else if(dir(i).gt.360.-diff.and.dir(i).le.360.+diff)then
  2182. if(id(i,jb).eq.1)then
  2183. iresult=1
  2184. if(imass.eq.1)then
  2185. massfld=max(massflx(i,jb),massflx(i,j))
  2186. endif
  2187. return
  2188. endif
  2189. !--- steering flow from the north-east
  2190. else if(dir(i).gt.45.-diff.and.dir(i).le.45.+diff)then
  2191. if(id(ib,jb).eq.1)then
  2192. iresult=1
  2193. if(imass.eq.1)then
  2194. massfld=max(massflx(ib,jb),massflx(i,j))
  2195. endif
  2196. return
  2197. endif
  2198. endif
  2199. END SUBROUTINE cup_direction2
  2200. SUBROUTINE cup_env(z,qes,he,hes,t,q,p,z1, &
  2201. psur,ierr,tcrit,itest,xl,cp, &
  2202. ids,ide, jds,jde, kds,kde, &
  2203. ims,ime, jms,jme, kms,kme, &
  2204. its,ite, jts,jte, kts,kte )
  2205. IMPLICIT NONE
  2206. integer &
  2207. ,intent (in ) :: &
  2208. ids,ide, jds,jde, kds,kde, &
  2209. ims,ime, jms,jme, kms,kme, &
  2210. its,ite, jts,jte, kts,kte
  2211. !
  2212. ! ierr error value, maybe modified in this routine
  2213. ! q = environmental mixing ratio
  2214. ! qes = environmental saturation mixing ratio
  2215. ! t = environmental temp
  2216. ! tv = environmental virtual temp
  2217. ! p = environmental pressure
  2218. ! z = environmental heights
  2219. ! he = environmental moist static energy
  2220. ! hes = environmental saturation moist static energy
  2221. ! psur = surface pressure
  2222. ! z1 = terrain elevation
  2223. !
  2224. !
  2225. real, dimension (its:ite,kts:kte) &
  2226. ,intent (in ) :: &
  2227. p,t
  2228. real, dimension (its:ite,kts:kte) &
  2229. ,intent (out ) :: &
  2230. he,hes,qes
  2231. real, dimension (its:ite,kts:kte) &
  2232. ,intent (inout) :: &
  2233. z,q
  2234. real, dimension (its:ite) &
  2235. ,intent (in ) :: &
  2236. psur,z1
  2237. real &
  2238. ,intent (in ) :: &
  2239. xl,cp
  2240. integer, dimension (its:ite) &
  2241. ,intent (inout) :: &
  2242. ierr
  2243. integer &
  2244. ,intent (in ) :: &
  2245. itest
  2246. !
  2247. ! local variables in this routine
  2248. !
  2249. integer :: &
  2250. i,k,iph
  2251. real, dimension (1:2) :: AE,BE,HT
  2252. real, dimension (its:ite,kts:kte) :: tv
  2253. real :: tcrit,e,tvbar
  2254. integer :: itf, ktf
  2255. itf=MIN(ite,ide-1)
  2256. ktf=MIN(kte,kde-1)
  2257. HT(1)=XL/CP
  2258. HT(2)=2.834E6/CP
  2259. BE(1)=.622*HT(1)/.286
  2260. AE(1)=BE(1)/273.+ALOG(610.71)
  2261. BE(2)=.622*HT(2)/.286
  2262. AE(2)=BE(2)/273.+ALOG(610.71)
  2263. ! print *, 'TCRIT = ', tcrit,its,ite
  2264. DO k=kts,ktf
  2265. do i=its,itf
  2266. if(ierr(i).eq.0)then
  2267. !Csgb - IPH is for phase, dependent on TCRIT (water or ice)
  2268. IPH=1
  2269. IF(T(I,K).LE.TCRIT)IPH=2
  2270. ! print *, 'AE(IPH),BE(IPH) = ',AE(IPH),BE(IPH),AE(IPH)-BE(IPH),T(i,k),i,k
  2271. E=EXP(AE(IPH)-BE(IPH)/T(I,K))
  2272. ! print *, 'P, E = ', P(I,K), E
  2273. QES(I,K)=.622*E/(100.*P(I,K)-E)
  2274. IF(QES(I,K).LE.1.E-08)QES(I,K)=1.E-08
  2275. IF(Q(I,K).GT.QES(I,K))Q(I,K)=QES(I,K)
  2276. TV(I,K)=T(I,K)+.608*Q(I,K)*T(I,K)
  2277. endif
  2278. enddo
  2279. enddo
  2280. !
  2281. !--- z's are calculated with changed h's and q's and t's
  2282. !--- if itest=2
  2283. !
  2284. if(itest.ne.2)then
  2285. do i=its,itf
  2286. if(ierr(i).eq.0)then
  2287. Z(I,1)=max(0.,Z1(I))-(ALOG(P(I,1))- &
  2288. ALOG(PSUR(I)))*287.*TV(I,1)/9.81
  2289. endif
  2290. enddo
  2291. ! --- calculate heights
  2292. DO K=kts+1,ktf
  2293. do i=its,itf
  2294. if(ierr(i).eq.0)then
  2295. TVBAR=.5*TV(I,K)+.5*TV(I,K-1)
  2296. Z(I,K)=Z(I,K-1)-(ALOG(P(I,K))- &
  2297. ALOG(P(I,K-1)))*287.*TVBAR/9.81
  2298. endif
  2299. enddo
  2300. enddo
  2301. else
  2302. do k=kts,ktf
  2303. do i=its,itf
  2304. if(ierr(i).eq.0)then
  2305. z(i,k)=(he(i,k)-1004.*t(i,k)-2.5e6*q(i,k))/9.81
  2306. z(i,k)=max(1.e-3,z(i,k))
  2307. endif
  2308. enddo
  2309. enddo
  2310. endif
  2311. !
  2312. !--- calculate moist static energy - HE
  2313. ! saturated moist static energy - HES
  2314. !
  2315. DO k=kts,ktf
  2316. do i=its,itf
  2317. if(ierr(i).eq.0)then
  2318. if(itest.eq.0)HE(I,K)=9.81*Z(I,K)+1004.*T(I,K)+2.5E06*Q(I,K)
  2319. HES(I,K)=9.81*Z(I,K)+1004.*T(I,K)+2.5E06*QES(I,K)
  2320. IF(HE(I,K).GE.HES(I,K))HE(I,K)=HES(I,K)
  2321. ! if(i.eq.2)then
  2322. ! print *,k,z(i,k),t(i,k),p(i,k),he(i,k),hes(i,k)
  2323. ! endif
  2324. endif
  2325. enddo
  2326. enddo
  2327. END SUBROUTINE cup_env
  2328. SUBROUTINE cup_env_clev(t,qes,q,he,hes,z,p,qes_cup,q_cup, &
  2329. he_cup,hes_cup,z_cup,p_cup,gamma_cup,t_cup,psur, &
  2330. ierr,z1,xl,rv,cp, &
  2331. ids,ide, jds,jde, kds,kde, &
  2332. ims,ime, jms,jme, kms,kme, &
  2333. its,ite, jts,jte, kts,kte )
  2334. IMPLICIT NONE
  2335. integer &
  2336. ,intent (in ) :: &
  2337. ids,ide, jds,jde, kds,kde, &
  2338. ims,ime, jms,jme, kms,kme, &
  2339. its,ite, jts,jte, kts,kte
  2340. !
  2341. ! ierr error value, maybe modified in this routine
  2342. ! q = environmental mixing ratio
  2343. ! q_cup = environmental mixing ratio on cloud levels
  2344. ! qes = environmental saturation mixing ratio
  2345. ! qes_cup = environmental saturation mixing ratio on cloud levels
  2346. ! t = environmental temp
  2347. ! t_cup = environmental temp on cloud levels
  2348. ! p = environmental pressure
  2349. ! p_cup = environmental pressure on cloud levels
  2350. ! z = environmental heights
  2351. ! z_cup = environmental heights on cloud levels
  2352. ! he = environmental moist static energy
  2353. ! he_cup = environmental moist static energy on cloud levels
  2354. ! hes = environmental saturation moist static energy
  2355. ! hes_cup = environmental saturation moist static energy on cloud levels
  2356. ! gamma_cup = gamma on cloud levels
  2357. ! psur = surface pressure
  2358. ! z1 = terrain elevation
  2359. !
  2360. !
  2361. real, dimension (its:ite,kts:kte) &
  2362. ,intent (in ) :: &
  2363. qes,q,he,hes,z,p,t
  2364. real, dimension (its:ite,kts:kte) &
  2365. ,intent (out ) :: &
  2366. qes_cup,q_cup,he_cup,hes_cup,z_cup,p_cup,gamma_cup,t_cup
  2367. real, dimension (its:ite) &
  2368. ,intent (in ) :: &
  2369. psur,z1
  2370. real &
  2371. ,intent (in ) :: &
  2372. xl,rv,cp
  2373. integer, dimension (its:ite) &
  2374. ,intent (inout) :: &
  2375. ierr
  2376. !
  2377. ! local variables in this routine
  2378. !
  2379. integer :: &
  2380. i,k
  2381. integer :: itf, ktf
  2382. itf=MIN(ite,ide-1)
  2383. ktf=MIN(kte,kde-1)
  2384. do k=kts+1,ktf
  2385. do i=its,itf
  2386. if(ierr(i).eq.0)then
  2387. qes_cup(i,k)=.5*(qes(i,k-1)+qes(i,k))
  2388. q_cup(i,k)=.5*(q(i,k-1)+q(i,k))
  2389. hes_cup(i,k)=.5*(hes(i,k-1)+hes(i,k))
  2390. he_cup(i,k)=.5*(he(i,k-1)+he(i,k))
  2391. if(he_cup(i,k).gt.hes_cup(i,k))he_cup(i,k)=hes_cup(i,k)
  2392. z_cup(i,k)=.5*(z(i,k-1)+z(i,k))
  2393. p_cup(i,k)=.5*(p(i,k-1)+p(i,k))
  2394. t_cup(i,k)=.5*(t(i,k-1)+t(i,k))
  2395. gamma_cup(i,k)=(xl/cp)*(xl/(rv*t_cup(i,k) &
  2396. *t_cup(i,k)))*qes_cup(i,k)
  2397. endif
  2398. enddo
  2399. enddo
  2400. do i=its,itf
  2401. if(ierr(i).eq.0)then
  2402. qes_cup(i,1)=qes(i,1)
  2403. q_cup(i,1)=q(i,1)
  2404. hes_cup(i,1)=hes(i,1)
  2405. he_cup(i,1)=he(i,1)
  2406. z_cup(i,1)=.5*(z(i,1)+z1(i))
  2407. p_cup(i,1)=.5*(p(i,1)+psur(i))
  2408. t_cup(i,1)=t(i,1)
  2409. gamma_cup(i,1)=xl/cp*(xl/(rv*t_cup(i,1) &
  2410. *t_cup(i,1)))*qes_cup(i,1)
  2411. endif
  2412. enddo
  2413. END SUBROUTINE cup_env_clev
  2414. SUBROUTINE cup_forcing_ens(closure_n,xland,aa0,aa1,xaa0,mbdt,dtime,ierr,ierr2,ierr3,&
  2415. xf_ens,j,name,maxens,iens,iedt,maxens2,maxens3,mconv, &
  2416. p_cup,ktop,omeg,zd,k22,zu,pr_ens,edt,kbcon,massflx, &
  2417. iact_old_gr,dir,ensdim,massfln,icoic, &
  2418. ids,ide, jds,jde, kds,kde, &
  2419. ims,ime, jms,jme, kms,kme, &
  2420. its,ite, jts,jte, kts,kte )
  2421. IMPLICIT NONE
  2422. integer &
  2423. ,intent (in ) :: &
  2424. ids,ide, jds,jde, kds,kde, &
  2425. ims,ime, jms,jme, kms,kme, &
  2426. its,ite, jts,jte, kts,kte
  2427. integer, intent (in ) :: &
  2428. j,ensdim,maxens,iens,iedt,maxens2,maxens3
  2429. !
  2430. ! ierr error value, maybe modified in this routine
  2431. ! pr_ens = precipitation ensemble
  2432. ! xf_ens = mass flux ensembles
  2433. ! massfln = downdraft mass flux ensembles used in next timestep
  2434. ! omeg = omega from large scale model
  2435. ! mconv = moisture convergence from large scale model
  2436. ! zd = downdraft normalized mass flux
  2437. ! zu = updraft normalized mass flux
  2438. ! aa0 = cloud work function without forcing effects
  2439. ! aa1 = cloud work function with forcing effects
  2440. ! xaa0 = cloud work function with cloud effects (ensemble dependent)
  2441. ! edt = epsilon
  2442. ! dir = "storm motion"
  2443. ! mbdt = arbitrary numerical parameter
  2444. ! dtime = dt over which forcing is applied
  2445. ! iact_gr_old = flag to tell where convection was active
  2446. ! kbcon = LFC of parcel from k22
  2447. ! k22 = updraft originating level
  2448. ! icoic = flag if only want one closure (usually set to zero!)
  2449. ! name = deep or shallow convection flag
  2450. !
  2451. real, dimension (ims:ime,jms:jme,1:ensdim) &
  2452. ,intent (inout) :: &
  2453. pr_ens
  2454. real, dimension (ims:ime,jms:jme,1:ensdim) &
  2455. ,intent (out ) :: &
  2456. xf_ens,massfln
  2457. real, dimension (ims:ime,jms:jme) &
  2458. ,intent (in ) :: &
  2459. massflx
  2460. real, dimension (its:ite,kts:kte) &
  2461. ,intent (in ) :: &
  2462. omeg,zd,zu,p_cup
  2463. real, dimension (its:ite,1:maxens) &
  2464. ,intent (in ) :: &
  2465. xaa0
  2466. real, dimension (its:ite) &
  2467. ,intent (in ) :: &
  2468. aa1,edt,dir,mconv,xland
  2469. real, dimension (its:ite) &
  2470. ,intent (inout) :: &
  2471. aa0,closure_n
  2472. real, dimension (1:maxens) &
  2473. ,intent (in ) :: &
  2474. mbdt
  2475. real &
  2476. ,intent (in ) :: &
  2477. dtime
  2478. integer, dimension (its:ite,jts:jte) &
  2479. ,intent (in ) :: &
  2480. iact_old_gr
  2481. integer, dimension (its:ite) &
  2482. ,intent (in ) :: &
  2483. k22,kbcon,ktop
  2484. integer, dimension (its:ite) &
  2485. ,intent (inout) :: &
  2486. ierr,ierr2,ierr3
  2487. integer &
  2488. ,intent (in ) :: &
  2489. icoic
  2490. character *(*), intent (in) :: &
  2491. name
  2492. !
  2493. ! local variables in this routine
  2494. !
  2495. real, dimension (1:maxens3) :: &
  2496. xff_ens3
  2497. real, dimension (1:maxens) :: &
  2498. xk
  2499. integer :: &
  2500. i,k,nall,n,ne,nens,nens3,iresult,iresultd,iresulte,mkxcrt,kclim
  2501. parameter (mkxcrt=15)
  2502. real :: &
  2503. a1,massfld,xff0,xomg,aclim1,aclim2,aclim3,aclim4
  2504. real, dimension(1:mkxcrt) :: &
  2505. pcrit,acrit,acritt
  2506. integer :: itf,nall2
  2507. itf=MIN(ite,ide-1)
  2508. DATA PCRIT/850.,800.,750.,700.,650.,600.,550.,500.,450.,400., &
  2509. 350.,300.,250.,200.,150./
  2510. DATA ACRIT/.0633,.0445,.0553,.0664,.075,.1082,.1521,.2216, &
  2511. .3151,.3677,.41,.5255,.7663,1.1686,1.6851/
  2512. ! GDAS DERIVED ACRIT
  2513. DATA ACRITT/.203,.515,.521,.566,.625,.665,.659,.688, &
  2514. .743,.813,.886,.947,1.138,1.377,1.896/
  2515. !
  2516. nens=0
  2517. !--- LARGE SCALE FORCING
  2518. !
  2519. DO 100 i=its,itf
  2520. ! if(i.eq.ipr.and.j.eq.jpr)print *,'ierr = ',ierr(i)
  2521. if(name.eq.'deeps'.and.ierr(i).gt.995)then
  2522. ! print *,i,j,ierr(i),aa0(i)
  2523. aa0(i)=0.
  2524. ierr(i)=0
  2525. endif
  2526. IF(ierr(i).eq.0)then
  2527. ! kclim=0
  2528. do k=mkxcrt,1,-1
  2529. if(p_cup(i,ktop(i)).lt.pcrit(k))then
  2530. kclim=k
  2531. go to 9
  2532. endif
  2533. enddo
  2534. if(p_cup(i,ktop(i)).ge.pcrit(1))kclim=1
  2535. 9 continue
  2536. kclim=max(kclim,1)
  2537. k=max(kclim-1,1)
  2538. aclim1=acrit(kclim)*1.e3
  2539. aclim2=acrit(k)*1.e3
  2540. aclim3=acritt(kclim)*1.e3
  2541. aclim4=acritt(k)*1.e3
  2542. ! print *,'p_cup(ktop(i)),kclim,pcrit(kclim)'
  2543. ! print *,p_cup(i,ktop(i)),kclim,pcrit(kclim)
  2544. ! print *,'aclim1,aclim2,aclim3,aclim4'
  2545. ! print *,aclim1,aclim2,aclim3,aclim4
  2546. ! print *,dtime,name,ierr(i),aa1(i),aa0(i)
  2547. ! print *,dtime,name,ierr(i),aa1(i),aa0(i)
  2548. !
  2549. !--- treatment different for this closure
  2550. !
  2551. if(name.eq.'deeps')then
  2552. !
  2553. xff0= (AA1(I)-AA0(I))/DTIME
  2554. xff_ens3(1)=(AA1(I)-AA0(I))/dtime
  2555. xff_ens3(2)=.9*xff_ens3(1)
  2556. xff_ens3(3)=1.1*xff_ens3(1)
  2557. !
  2558. !--- more original Arakawa-Schubert (climatologic value of aa0)
  2559. !
  2560. !
  2561. !--- omeg is in bar/s, mconv done with omeg in Pa/s
  2562. ! more like Brown (1979), or Frank-Cohen (199?)
  2563. !
  2564. xff_ens3(4)=-omeg(i,k22(i))/9.81
  2565. xff_ens3(5)=-omeg(i,kbcon(i))/9.81
  2566. xff_ens3(6)=-omeg(i,1)/9.81
  2567. do k=2,kbcon(i)-1
  2568. xomg=-omeg(i,k)/9.81
  2569. if(xomg.gt.xff_ens3(6))xff_ens3(6)=xomg
  2570. enddo
  2571. !
  2572. !--- more like Krishnamurti et al.
  2573. !
  2574. xff_ens3(7)=mconv(i)
  2575. xff_ens3(8)=mconv(i)
  2576. xff_ens3(9)=mconv(i)
  2577. !
  2578. !--- more like Fritsch Chappel or Kain Fritsch (plus triggers)
  2579. !
  2580. xff_ens3(10)=AA1(I)/(60.*20.)
  2581. xff_ens3(11)=AA1(I)/(60.*30.)
  2582. xff_ens3(12)=AA1(I)/(60.*40.)
  2583. !
  2584. !--- more original Arakawa-Schubert (climatologic value of aa0)
  2585. !
  2586. xff_ens3(13)=max(0.,(AA1(I)-aclim1)/dtime)
  2587. xff_ens3(14)=max(0.,(AA1(I)-aclim2)/dtime)
  2588. xff_ens3(15)=max(0.,(AA1(I)-aclim3)/dtime)
  2589. xff_ens3(16)=max(0.,(AA1(I)-aclim4)/dtime)
  2590. ! if(ierr2(i).gt.0.or.ierr3(i).gt.0)then
  2591. ! xff_ens3(10)=0.
  2592. ! xff_ens3(11)=0.
  2593. ! xff_ens3(12)=0.
  2594. ! xff_ens3(13)=0.
  2595. ! xff_ens3(14)=0.
  2596. ! xff_ens3(15)=0.
  2597. ! xff_ens3(16)=0.
  2598. ! endif
  2599. do nens=1,maxens
  2600. XK(nens)=(XAA0(I,nens)-AA1(I))/MBDT(2)
  2601. if(xk(nens).le.0.and.xk(nens).gt.-1.e-6) &
  2602. xk(nens)=-1.e-6
  2603. if(xk(nens).gt.0.and.xk(nens).lt.1.e-6) &
  2604. xk(nens)=1.e-6
  2605. enddo
  2606. !
  2607. !--- add up all ensembles
  2608. !
  2609. do 350 ne=1,maxens
  2610. !
  2611. !--- for every xk, we have maxens3 xffs
  2612. !--- iens is from outermost ensemble (most expensive!
  2613. !
  2614. !--- iedt (maxens2 belongs to it)
  2615. !--- is from second, next outermost, not so expensive
  2616. !
  2617. !--- so, for every outermost loop, we have maxens*maxens2*3
  2618. !--- ensembles!!! nall would be 0, if everything is on first
  2619. !--- loop index, then ne would start counting, then iedt, then iens....
  2620. !
  2621. iresult=0
  2622. iresultd=0
  2623. iresulte=0
  2624. nall=(iens-1)*maxens3*maxens*maxens2 &
  2625. +(iedt-1)*maxens*maxens3 &
  2626. +(ne-1)*maxens3
  2627. !
  2628. ! over water, enfor!e small cap for some of the closures
  2629. !
  2630. if(xland(i).lt.0.1)then
  2631. if(ierr2(i).gt.0.or.ierr3(i).gt.0)then
  2632. ! - ierr2 - 75 mb cap thickness, ierr3 - 125 cap thickness
  2633. ! - for larger cap, set Grell closure to zero
  2634. xff_ens3(1) =0.
  2635. massfln(i,j,nall+1)=0.
  2636. xff_ens3(2) =0.
  2637. massfln(i,j,nall+2)=0.
  2638. xff_ens3(3) =0.
  2639. massfln(i,j,nall+3)=0.
  2640. closure_n(i)=closure_n(i)-1.
  2641. xff_ens3(7) =0.
  2642. massfln(i,j,nall+7)=0.
  2643. xff_ens3(8) =0.
  2644. massfln(i,j,nall+8)=0.
  2645. xff_ens3(9) =0.
  2646. ! massfln(i,j,nall+9)=0.
  2647. closure_n(i)=closure_n(i)-1.
  2648. endif
  2649. !
  2650. ! also take out some closures in general
  2651. !
  2652. xff_ens3(4) =0.
  2653. massfln(i,j,nall+4)=0.
  2654. xff_ens3(5) =0.
  2655. massfln(i,j,nall+5)=0.
  2656. xff_ens3(6) =0.
  2657. massfln(i,j,nall+6)=0.
  2658. closure_n(i)=closure_n(i)-3.
  2659. xff_ens3(10)=0.
  2660. massfln(i,j,nall+10)=0.
  2661. xff_ens3(11)=0.
  2662. massfln(i,j,nall+11)=0.
  2663. xff_ens3(12)=0.
  2664. massfln(i,j,nall+12)=0.
  2665. if(ne.eq.1)closure_n(i)=closure_n(i)-3
  2666. xff_ens3(13)=0.
  2667. massfln(i,j,nall+13)=0.
  2668. xff_ens3(14)=0.
  2669. massfln(i,j,nall+14)=0.
  2670. xff_ens3(15)=0.
  2671. massfln(i,j,nall+15)=0.
  2672. massfln(i,j,nall+16)=0.
  2673. if(ne.eq.1)closure_n(i)=closure_n(i)-4
  2674. endif
  2675. !
  2676. ! end water treatment
  2677. !
  2678. !--- check for upwind convection
  2679. ! iresult=0
  2680. massfld=0.
  2681. ! call cup_direction2(i,j,dir,iact_old_gr, &
  2682. ! massflx,iresult,1, &
  2683. ! massfld, &
  2684. ! ids,ide, jds,jde, kds,kde, &
  2685. ! ims,ime, jms,jme, kms,kme, &
  2686. ! its,ite, jts,jte, kts,kte )
  2687. ! if(i.eq.ipr.and.j.eq.jpr.and.iedt.eq.1.and.ne.eq.1)then
  2688. ! if(iedt.eq.1.and.ne.eq.1)then
  2689. ! print *,massfld,ne,iedt,iens
  2690. ! print *,xk(ne),xff_ens3(1),xff_ens3(2),xff_ens3(3)
  2691. ! endif
  2692. ! print *,i,j,massfld,aa0(i),aa1(i)
  2693. IF(XK(ne).lt.0.and.xff0.gt.0.)iresultd=1
  2694. iresulte=max(iresult,iresultd)
  2695. iresulte=1
  2696. if(iresulte.eq.1)then
  2697. !
  2698. !--- special treatment for stability closures
  2699. !
  2700. if(xff0.gt.0.)then
  2701. xf_ens(i,j,nall+1)=max(0.,-xff_ens3(1)/xk(ne)) &
  2702. +massfld
  2703. xf_ens(i,j,nall+2)=max(0.,-xff_ens3(2)/xk(ne)) &
  2704. +massfld
  2705. xf_ens(i,j,nall+3)=max(0.,-xff_ens3(3)/xk(ne)) &
  2706. +massfld
  2707. xf_ens(i,j,nall+13)=max(0.,-xff_ens3(13)/xk(ne)) &
  2708. +massfld
  2709. xf_ens(i,j,nall+14)=max(0.,-xff_ens3(14)/xk(ne)) &
  2710. +massfld
  2711. xf_ens(i,j,nall+15)=max(0.,-xff_ens3(15)/xk(ne)) &
  2712. +massfld
  2713. xf_ens(i,j,nall+16)=max(0.,-xff_ens3(16)/xk(ne)) &
  2714. +massfld
  2715. else
  2716. xf_ens(i,j,nall+1)=massfld
  2717. xf_ens(i,j,nall+2)=massfld
  2718. xf_ens(i,j,nall+3)=massfld
  2719. xf_ens(i,j,nall+13)=massfld
  2720. xf_ens(i,j,nall+14)=massfld
  2721. xf_ens(i,j,nall+15)=massfld
  2722. xf_ens(i,j,nall+16)=massfld
  2723. endif
  2724. !
  2725. !--- if iresult.eq.1, following independent of xff0
  2726. !
  2727. xf_ens(i,j,nall+4)=max(0.,xff_ens3(4) &
  2728. +massfld)
  2729. xf_ens(i,j,nall+5)=max(0.,xff_ens3(5) &
  2730. +massfld)
  2731. xf_ens(i,j,nall+6)=max(0.,xff_ens3(6) &
  2732. +massfld)
  2733. a1=max(1.e-3,pr_ens(i,j,nall+7))
  2734. xf_ens(i,j,nall+7)=max(0.,xff_ens3(7) &
  2735. /a1)
  2736. a1=max(1.e-3,pr_ens(i,j,nall+8))
  2737. xf_ens(i,j,nall+8)=max(0.,xff_ens3(8) &
  2738. /a1)
  2739. a1=max(1.e-3,pr_ens(i,j,nall+9))
  2740. xf_ens(i,j,nall+9)=max(0.,xff_ens3(9) &
  2741. /a1)
  2742. if(XK(ne).lt.0.)then
  2743. xf_ens(i,j,nall+10)=max(0., &
  2744. -xff_ens3(10)/xk(ne)) &
  2745. +massfld
  2746. xf_ens(i,j,nall+11)=max(0., &
  2747. -xff_ens3(11)/xk(ne)) &
  2748. +massfld
  2749. xf_ens(i,j,nall+12)=max(0., &
  2750. -xff_ens3(12)/xk(ne)) &
  2751. +massfld
  2752. else
  2753. xf_ens(i,j,nall+10)=massfld
  2754. xf_ens(i,j,nall+11)=massfld
  2755. xf_ens(i,j,nall+12)=massfld
  2756. endif
  2757. if(icoic.ge.1)then
  2758. closure_n(i)=0.
  2759. xf_ens(i,j,nall+1)=xf_ens(i,j,nall+icoic)
  2760. xf_ens(i,j,nall+2)=xf_ens(i,j,nall+icoic)
  2761. xf_ens(i,j,nall+3)=xf_ens(i,j,nall+icoic)
  2762. xf_ens(i,j,nall+4)=xf_ens(i,j,nall+icoic)
  2763. xf_ens(i,j,nall+5)=xf_ens(i,j,nall+icoic)
  2764. xf_ens(i,j,nall+6)=xf_ens(i,j,nall+icoic)
  2765. xf_ens(i,j,nall+7)=xf_ens(i,j,nall+icoic)
  2766. xf_ens(i,j,nall+8)=xf_ens(i,j,nall+icoic)
  2767. xf_ens(i,j,nall+9)=xf_ens(i,j,nall+icoic)
  2768. xf_ens(i,j,nall+10)=xf_ens(i,j,nall+icoic)
  2769. xf_ens(i,j,nall+11)=xf_ens(i,j,nall+icoic)
  2770. xf_ens(i,j,nall+12)=xf_ens(i,j,nall+icoic)
  2771. xf_ens(i,j,nall+13)=xf_ens(i,j,nall+icoic)
  2772. xf_ens(i,j,nall+14)=xf_ens(i,j,nall+icoic)
  2773. xf_ens(i,j,nall+15)=xf_ens(i,j,nall+icoic)
  2774. xf_ens(i,j,nall+16)=xf_ens(i,j,nall+icoic)
  2775. endif
  2776. !
  2777. ! replace 13-16 for now with other stab closures
  2778. ! (13 gave problems for mass model)
  2779. !
  2780. ! xf_ens(i,j,nall+13)=xf_ens(i,j,nall+1)
  2781. if(icoic.eq.0)xf_ens(i,j,nall+14)=xf_ens(i,j,nall+13)
  2782. ! xf_ens(i,j,nall+15)=xf_ens(i,j,nall+11)
  2783. ! xf_ens(i,j,nall+16)=xf_ens(i,j,nall+11)
  2784. ! xf_ens(i,j,nall+7)=xf_ens(i,j,nall+4)
  2785. ! xf_ens(i,j,nall+8)=xf_ens(i,j,nall+5)
  2786. ! xf_ens(i,j,nall+9)=xf_ens(i,j,nall+6)
  2787. !
  2788. !--- store new for next time step
  2789. !
  2790. do nens3=1,maxens3
  2791. massfln(i,j,nall+nens3)=edt(i) &
  2792. *xf_ens(i,j,nall+nens3)
  2793. massfln(i,j,nall+nens3)=max(0., &
  2794. massfln(i,j,nall+nens3))
  2795. enddo
  2796. !
  2797. !
  2798. !--- do some more on the caps!!! ne=1 for 175, ne=2 for 100,....
  2799. !
  2800. ! do not care for caps here for closure groups 1 and 5,
  2801. ! they are fine, do not turn them off here
  2802. !
  2803. !
  2804. if(ne.eq.2.and.ierr2(i).gt.0)then
  2805. xf_ens(i,j,nall+1) =0.
  2806. xf_ens(i,j,nall+2) =0.
  2807. xf_ens(i,j,nall+3) =0.
  2808. xf_ens(i,j,nall+4) =0.
  2809. xf_ens(i,j,nall+5) =0.
  2810. xf_ens(i,j,nall+6) =0.
  2811. xf_ens(i,j,nall+7) =0.
  2812. xf_ens(i,j,nall+8) =0.
  2813. xf_ens(i,j,nall+9) =0.
  2814. xf_ens(i,j,nall+10)=0.
  2815. xf_ens(i,j,nall+11)=0.
  2816. xf_ens(i,j,nall+12)=0.
  2817. xf_ens(i,j,nall+13)=0.
  2818. xf_ens(i,j,nall+14)=0.
  2819. xf_ens(i,j,nall+15)=0.
  2820. xf_ens(i,j,nall+16)=0.
  2821. massfln(i,j,nall+1)=0.
  2822. massfln(i,j,nall+2)=0.
  2823. massfln(i,j,nall+3)=0.
  2824. massfln(i,j,nall+4)=0.
  2825. massfln(i,j,nall+5)=0.
  2826. massfln(i,j,nall+6)=0.
  2827. massfln(i,j,nall+7)=0.
  2828. massfln(i,j,nall+8)=0.
  2829. massfln(i,j,nall+9)=0.
  2830. massfln(i,j,nall+10)=0.
  2831. massfln(i,j,nall+11)=0.
  2832. massfln(i,j,nall+12)=0.
  2833. massfln(i,j,nall+13)=0.
  2834. massfln(i,j,nall+14)=0.
  2835. massfln(i,j,nall+15)=0.
  2836. massfln(i,j,nall+16)=0.
  2837. endif
  2838. if(ne.eq.3.and.ierr3(i).gt.0)then
  2839. xf_ens(i,j,nall+1) =0.
  2840. xf_ens(i,j,nall+2) =0.
  2841. xf_ens(i,j,nall+3) =0.
  2842. xf_ens(i,j,nall+4) =0.
  2843. xf_ens(i,j,nall+5) =0.
  2844. xf_ens(i,j,nall+6) =0.
  2845. xf_ens(i,j,nall+7) =0.
  2846. xf_ens(i,j,nall+8) =0.
  2847. xf_ens(i,j,nall+9) =0.
  2848. xf_ens(i,j,nall+10)=0.
  2849. xf_ens(i,j,nall+11)=0.
  2850. xf_ens(i,j,nall+12)=0.
  2851. xf_ens(i,j,nall+13)=0.
  2852. xf_ens(i,j,nall+14)=0.
  2853. xf_ens(i,j,nall+15)=0.
  2854. xf_ens(i,j,nall+16)=0.
  2855. massfln(i,j,nall+1)=0.
  2856. massfln(i,j,nall+2)=0.
  2857. massfln(i,j,nall+3)=0.
  2858. massfln(i,j,nall+4)=0.
  2859. massfln(i,j,nall+5)=0.
  2860. massfln(i,j,nall+6)=0.
  2861. massfln(i,j,nall+7)=0.
  2862. massfln(i,j,nall+8)=0.
  2863. massfln(i,j,nall+9)=0.
  2864. massfln(i,j,nall+10)=0.
  2865. massfln(i,j,nall+11)=0.
  2866. massfln(i,j,nall+12)=0.
  2867. massfln(i,j,nall+13)=0.
  2868. massfln(i,j,nall+14)=0.
  2869. massfln(i,j,nall+15)=0.
  2870. massfln(i,j,nall+16)=0.
  2871. endif
  2872. endif
  2873. 350 continue
  2874. ! ne=1, cap=175
  2875. !
  2876. nall=(iens-1)*maxens3*maxens*maxens2 &
  2877. +(iedt-1)*maxens*maxens3
  2878. ! ne=2, cap=100
  2879. !
  2880. nall2=(iens-1)*maxens3*maxens*maxens2 &
  2881. +(iedt-1)*maxens*maxens3 &
  2882. +(2-1)*maxens3
  2883. xf_ens(i,j,nall+4) = xf_ens(i,j,nall2+4)
  2884. xf_ens(i,j,nall+5) =xf_ens(i,j,nall2+5)
  2885. xf_ens(i,j,nall+6) =xf_ens(i,j,nall2+6)
  2886. xf_ens(i,j,nall+7) =xf_ens(i,j,nall2+7)
  2887. xf_ens(i,j,nall+8) =xf_ens(i,j,nall2+8)
  2888. xf_ens(i,j,nall+9) =xf_ens(i,j,nall2+9)
  2889. xf_ens(i,j,nall+10)=xf_ens(i,j,nall2+10)
  2890. xf_ens(i,j,nall+11)=xf_ens(i,j,nall2+11)
  2891. xf_ens(i,j,nall+12)=xf_ens(i,j,nall2+12)
  2892. go to 100
  2893. endif
  2894. elseif(ierr(i).ne.20.and.ierr(i).ne.0)then
  2895. do n=1,ensdim
  2896. xf_ens(i,j,n)=0.
  2897. massfln(i,j,n)=0.
  2898. enddo
  2899. endif
  2900. 100 continue
  2901. END SUBROUTINE cup_forcing_ens
  2902. SUBROUTINE cup_kbcon(cap_inc,iloop,k22,kbcon,he_cup,hes_cup, &
  2903. ierr,kbmax,p_cup,cap_max, &
  2904. ids,ide, jds,jde, kds,kde, &
  2905. ims,ime, jms,jme, kms,kme, &
  2906. its,ite, jts,jte, kts,kte )
  2907. IMPLICIT NONE
  2908. !
  2909. ! only local wrf dimensions are need as of now in this routine
  2910. integer &
  2911. ,intent (in ) :: &
  2912. ids,ide, jds,jde, kds,kde, &
  2913. ims,ime, jms,jme, kms,kme, &
  2914. its,ite, jts,jte, kts,kte
  2915. !
  2916. !
  2917. !
  2918. ! ierr error value, maybe modified in this routine
  2919. !
  2920. real, dimension (its:ite,kts:kte) &
  2921. ,intent (in ) :: &
  2922. he_cup,hes_cup,p_cup
  2923. real, dimension (its:ite) &
  2924. ,intent (in ) :: &
  2925. cap_max,cap_inc
  2926. integer, dimension (its:ite) &
  2927. ,intent (in ) :: &
  2928. kbmax
  2929. integer, dimension (its:ite) &
  2930. ,intent (inout) :: &
  2931. kbcon,k22,ierr
  2932. integer &
  2933. ,intent (in ) :: &
  2934. iloop
  2935. !
  2936. ! local variables in this routine
  2937. !
  2938. integer :: &
  2939. i
  2940. real :: &
  2941. pbcdif,plus
  2942. integer :: itf
  2943. itf=MIN(ite,ide-1)
  2944. !
  2945. !--- DETERMINE THE LEVEL OF CONVECTIVE CLOUD BASE - KBCON
  2946. !
  2947. DO 27 i=its,itf
  2948. kbcon(i)=1
  2949. IF(ierr(I).ne.0)GO TO 27
  2950. KBCON(I)=K22(I)
  2951. GO TO 32
  2952. 31 CONTINUE
  2953. KBCON(I)=KBCON(I)+1
  2954. IF(KBCON(I).GT.KBMAX(i)+2)THEN
  2955. if(iloop.lt.4)ierr(i)=3
  2956. ! if(iloop.lt.4)ierr(i)=997
  2957. GO TO 27
  2958. ENDIF
  2959. 32 CONTINUE
  2960. IF(HE_cup(I,K22(I)).LT.HES_cup(I,KBCON(I)))GO TO 31
  2961. ! cloud base pressure and max moist static energy pressure
  2962. ! i.e., the depth (in mb) of the layer of negative buoyancy
  2963. if(KBCON(I)-K22(I).eq.1)go to 27
  2964. PBCDIF=-P_cup(I,KBCON(I))+P_cup(I,K22(I))
  2965. plus=max(25.,cap_max(i)-float(iloop-1)*cap_inc(i))
  2966. if(iloop.eq.4)plus=cap_max(i)
  2967. IF(PBCDIF.GT.plus)THEN
  2968. K22(I)=K22(I)+1
  2969. KBCON(I)=K22(I)
  2970. GO TO 32
  2971. ENDIF
  2972. 27 CONTINUE
  2973. END SUBROUTINE cup_kbcon
  2974. SUBROUTINE cup_kbcon_cin(iloop,k22,kbcon,he_cup,hes_cup, &
  2975. z,tmean,qes,ierr,kbmax,p_cup,cap_max,xl,cp, &
  2976. ids,ide, jds,jde, kds,kde, &
  2977. ims,ime, jms,jme, kms,kme, &
  2978. its,ite, jts,jte, kts,kte )
  2979. IMPLICIT NONE
  2980. !
  2981. ! only local wrf dimensions are need as of now in this routine
  2982. integer &
  2983. ,intent (in ) :: &
  2984. ids,ide, jds,jde, kds,kde, &
  2985. ims,ime, jms,jme, kms,kme, &
  2986. its,ite, jts,jte, kts,kte
  2987. !
  2988. !
  2989. !
  2990. ! ierr error value, maybe modified in this routine
  2991. !
  2992. real, dimension (its:ite,kts:kte) &
  2993. ,intent (in ) :: &
  2994. he_cup,hes_cup,p_cup,z,tmean,qes
  2995. real, dimension (its:ite) &
  2996. ,intent (in ) :: &
  2997. cap_max
  2998. real &
  2999. ,intent (in ) :: &
  3000. xl,cp
  3001. integer, dimension (its:ite) &
  3002. ,intent (in ) :: &
  3003. kbmax
  3004. integer, dimension (its:ite) &
  3005. ,intent (inout) :: &
  3006. kbcon,k22,ierr
  3007. integer &
  3008. ,intent (in ) :: &
  3009. iloop
  3010. !
  3011. ! local variables in this routine
  3012. !
  3013. integer :: &
  3014. i,k
  3015. real :: &
  3016. cin,cin_max,dh,tprim,gamma
  3017. !
  3018. integer :: itf
  3019. itf=MIN(ite,ide-1)
  3020. !
  3021. !
  3022. !
  3023. !--- DETERMINE THE LEVEL OF CONVECTIVE CLOUD BASE - KBCON
  3024. !
  3025. DO 27 i=its,itf
  3026. cin_max=-cap_max(i)
  3027. kbcon(i)=1
  3028. cin = 0.
  3029. IF(ierr(I).ne.0)GO TO 27
  3030. KBCON(I)=K22(I)
  3031. GO TO 32
  3032. 31 CONTINUE
  3033. KBCON(I)=KBCON(I)+1
  3034. IF(KBCON(I).GT.KBMAX(i)+2)THEN
  3035. if(iloop.eq.1)ierr(i)=3
  3036. ! if(iloop.eq.2)ierr(i)=997
  3037. GO TO 27
  3038. ENDIF
  3039. 32 CONTINUE
  3040. dh = HE_cup(I,K22(I)) - HES_cup(I,KBCON(I))
  3041. if (dh.lt. 0.) then
  3042. GAMMA=(xl/cp)*(xl/(461.525*(Tmean(I,K22(i))**2)))*QES(I,K22(i))
  3043. tprim = dh/(cp*(1.+gamma))
  3044. cin = cin + 9.8066 * tprim &
  3045. *(z(i,k22(i))-z(i,k22(i)-1)) / tmean(i,k22(i))
  3046. go to 31
  3047. end if
  3048. ! If negative energy in negatively buoyant layer
  3049. ! exceeds convective inhibition (CIN) threshold,
  3050. ! then set K22 level one level up and see if that
  3051. ! will work.
  3052. IF(cin.lT.cin_max)THEN
  3053. ! print *,i,cin,cin_max,k22(i),kbcon(i)
  3054. K22(I)=K22(I)+1
  3055. KBCON(I)=K22(I)
  3056. GO TO 32
  3057. ENDIF
  3058. 27 CONTINUE
  3059. END SUBROUTINE cup_kbcon_cin
  3060. SUBROUTINE cup_ktop(ilo,dby,kbcon,ktop,ierr, &
  3061. ids,ide, jds,jde, kds,kde, &
  3062. ims,ime, jms,jme, kms,kme, &
  3063. its,ite, jts,jte, kts,kte )
  3064. IMPLICIT NONE
  3065. !
  3066. ! on input
  3067. !
  3068. ! only local wrf dimensions are need as of now in this routine
  3069. integer &
  3070. ,intent (in ) :: &
  3071. ids,ide, jds,jde, kds,kde, &
  3072. ims,ime, jms,jme, kms,kme, &
  3073. its,ite, jts,jte, kts,kte
  3074. ! dby = buoancy term
  3075. ! ktop = cloud top (output)
  3076. ! ilo = flag
  3077. ! ierr error value, maybe modified in this routine
  3078. !
  3079. real, dimension (its:ite,kts:kte) &
  3080. ,intent (inout) :: &
  3081. dby
  3082. integer, dimension (its:ite) &
  3083. ,intent (in ) :: &
  3084. kbcon
  3085. integer &
  3086. ,intent (in ) :: &
  3087. ilo
  3088. integer, dimension (its:ite) &
  3089. ,intent (out ) :: &
  3090. ktop
  3091. integer, dimension (its:ite) &
  3092. ,intent (inout) :: &
  3093. ierr
  3094. !
  3095. ! local variables in this routine
  3096. !
  3097. integer :: &
  3098. i,k
  3099. !
  3100. integer :: itf, ktf
  3101. itf=MIN(ite,ide-1)
  3102. ktf=MIN(kte,kde-1)
  3103. !
  3104. !
  3105. DO 42 i=its,itf
  3106. ktop(i)=1
  3107. IF(ierr(I).EQ.0)then
  3108. DO 40 K=KBCON(I)+1,ktf-1
  3109. IF(DBY(I,K).LE.0.)THEN
  3110. KTOP(I)=K-1
  3111. GO TO 41
  3112. ENDIF
  3113. 40 CONTINUE
  3114. if(ilo.eq.1)ierr(i)=5
  3115. ! if(ilo.eq.2)ierr(i)=998
  3116. GO TO 42
  3117. 41 CONTINUE
  3118. do k=ktop(i)+1,ktf
  3119. dby(i,k)=0.
  3120. enddo
  3121. endif
  3122. 42 CONTINUE
  3123. END SUBROUTINE cup_ktop
  3124. SUBROUTINE cup_MAXIMI(ARRAY,KS,KE,MAXX,ierr, &
  3125. ids,ide, jds,jde, kds,kde, &
  3126. ims,ime, jms,jme, kms,kme, &
  3127. its,ite, jts,jte, kts,kte )
  3128. IMPLICIT NONE
  3129. !
  3130. ! on input
  3131. !
  3132. ! only local wrf dimensions are need as of now in this routine
  3133. integer &
  3134. ,intent (in ) :: &
  3135. ids,ide, jds,jde, kds,kde, &
  3136. ims,ime, jms,jme, kms,kme, &
  3137. its,ite, jts,jte, kts,kte
  3138. ! array input array
  3139. ! x output array with return values
  3140. ! kt output array of levels
  3141. ! ks,kend check-range
  3142. real, dimension (its:ite,kts:kte) &
  3143. ,intent (in ) :: &
  3144. array
  3145. integer, dimension (its:ite) &
  3146. ,intent (in ) :: &
  3147. ierr,ke
  3148. integer &
  3149. ,intent (in ) :: &
  3150. ks
  3151. integer, dimension (its:ite) &
  3152. ,intent (out ) :: &
  3153. maxx
  3154. real, dimension (its:ite) :: &
  3155. x
  3156. real :: &
  3157. xar
  3158. integer :: &
  3159. i,k
  3160. integer :: itf
  3161. itf=MIN(ite,ide-1)
  3162. DO 200 i=its,itf
  3163. MAXX(I)=KS
  3164. if(ierr(i).eq.0)then
  3165. X(I)=ARRAY(I,KS)
  3166. !
  3167. DO 100 K=KS,KE(i)
  3168. XAR=ARRAY(I,K)
  3169. IF(XAR.GE.X(I)) THEN
  3170. X(I)=XAR
  3171. MAXX(I)=K
  3172. ENDIF
  3173. 100 CONTINUE
  3174. endif
  3175. 200 CONTINUE
  3176. END SUBROUTINE cup_MAXIMI
  3177. SUBROUTINE cup_minimi(ARRAY,KS,KEND,KT,ierr, &
  3178. ids,ide, jds,jde, kds,kde, &
  3179. ims,ime, jms,jme, kms,kme, &
  3180. its,ite, jts,jte, kts,kte )
  3181. IMPLICIT NONE
  3182. !
  3183. ! on input
  3184. !
  3185. ! only local wrf dimensions are need as of now in this routine
  3186. integer &
  3187. ,intent (in ) :: &
  3188. ids,ide, jds,jde, kds,kde, &
  3189. ims,ime, jms,jme, kms,kme, &
  3190. its,ite, jts,jte, kts,kte
  3191. ! array input array
  3192. ! x output array with return values
  3193. ! kt output array of levels
  3194. ! ks,kend check-range
  3195. real, dimension (its:ite,kts:kte) &
  3196. ,intent (in ) :: &
  3197. array
  3198. integer, dimension (its:ite) &
  3199. ,intent (in ) :: &
  3200. ierr,ks,kend
  3201. integer, dimension (its:ite) &
  3202. ,intent (out ) :: &
  3203. kt
  3204. real, dimension (its:ite) :: &
  3205. x
  3206. integer :: &
  3207. i,k,kstop
  3208. integer :: itf
  3209. itf=MIN(ite,ide-1)
  3210. DO 200 i=its,itf
  3211. KT(I)=KS(I)
  3212. if(ierr(i).eq.0)then
  3213. X(I)=ARRAY(I,KS(I))
  3214. KSTOP=MAX(KS(I)+1,KEND(I))
  3215. !
  3216. DO 100 K=KS(I)+1,KSTOP
  3217. IF(ARRAY(I,K).LT.X(I)) THEN
  3218. X(I)=ARRAY(I,K)
  3219. KT(I)=K
  3220. ENDIF
  3221. 100 CONTINUE
  3222. endif
  3223. 200 CONTINUE
  3224. END SUBROUTINE cup_MINIMI
  3225. SUBROUTINE cup_output_ens(xf_ens,ierr,dellat,dellaq,dellaqc, &
  3226. outtem,outq,outqc,pre,pw,xmb,ktop, &
  3227. j,name,nx,nx2,iens,ierr2,ierr3,pr_ens, &
  3228. maxens3,ensdim,massfln, &
  3229. APR_GR,APR_W,APR_MC,APR_ST,APR_AS, &
  3230. APR_CAPMA,APR_CAPME,APR_CAPMI,closure_n,xland1, &
  3231. outCFU1,outCFD1,outDFU1,outEFU1,outDFD1,outEFD1, &
  3232. CFU1_ens,CFD1_ens,DFU1_ens,EFU1_ens,DFD1_ens,EFD1_ens, &
  3233. l_flux, &
  3234. ids,ide, jds,jde, kds,kde, &
  3235. ims,ime, jms,jme, kms,kme, &
  3236. its,ite, jts,jte, kts,kte)
  3237. IMPLICIT NONE
  3238. !
  3239. ! on input
  3240. !
  3241. ! only local wrf dimensions are need as of now in this routine
  3242. integer &
  3243. ,intent (in ) :: &
  3244. ids,ide, jds,jde, kds,kde, &
  3245. ims,ime, jms,jme, kms,kme, &
  3246. its,ite, jts,jte, kts,kte
  3247. integer, intent (in ) :: &
  3248. j,ensdim,nx,nx2,iens,maxens3
  3249. ! xf_ens = ensemble mass fluxes
  3250. ! pr_ens = precipitation ensembles
  3251. ! dellat = change of temperature per unit mass flux of cloud ensemble
  3252. ! dellaq = change of q per unit mass flux of cloud ensemble
  3253. ! dellaqc = change of qc per unit mass flux of cloud ensemble
  3254. ! outtem = output temp tendency (per s)
  3255. ! outq = output q tendency (per s)
  3256. ! outqc = output qc tendency (per s)
  3257. ! pre = output precip
  3258. ! xmb = total base mass flux
  3259. ! xfac1 = correction factor
  3260. ! pw = pw -epsilon*pd (ensemble dependent)
  3261. ! ierr error value, maybe modified in this routine
  3262. !
  3263. real, dimension (ims:ime,jms:jme,1:ensdim) &
  3264. ,intent (inout) :: &
  3265. xf_ens,pr_ens,massfln
  3266. real, dimension (ims:ime,jms:jme) &
  3267. ,intent (inout) :: &
  3268. APR_GR,APR_W,APR_MC,APR_ST,APR_AS,APR_CAPMA, &
  3269. APR_CAPME,APR_CAPMI
  3270. real, dimension (its:ite,kts:kte) &
  3271. ,intent (out ) :: &
  3272. outtem,outq,outqc
  3273. real, dimension (its:ite) &
  3274. ,intent (out ) :: &
  3275. pre,xmb
  3276. real, dimension (its:ite) &
  3277. ,intent (inout ) :: &
  3278. closure_n,xland1
  3279. real, dimension (its:ite,kts:kte,1:nx) &
  3280. ,intent (in ) :: &
  3281. dellat,dellaqc,dellaq,pw
  3282. integer, dimension (its:ite) &
  3283. ,intent (in ) :: &
  3284. ktop
  3285. integer, dimension (its:ite) &
  3286. ,intent (inout) :: &
  3287. ierr,ierr2,ierr3
  3288. real, dimension (its:ite,kts:kte,1:ensdim) &
  3289. ,intent (in ) :: &
  3290. CFU1_ens,CFD1_ens,DFU1_ens,EFU1_ens,DFD1_ens,EFD1_ens
  3291. real, dimension (its:ite,kts:kte) &
  3292. ,intent (out) :: &
  3293. outCFU1,outCFD1,outDFU1,outEFU1,outDFD1,outEFD1
  3294. logical, intent(in) :: l_flux
  3295. !
  3296. ! local variables in this routine
  3297. !
  3298. integer :: &
  3299. i,k,n,ncount
  3300. real :: &
  3301. outtes,ddtes,dtt,dtq,dtqc,dtpw,tuning,prerate,clos_wei
  3302. real, dimension (its:ite) :: &
  3303. xfac1
  3304. real, dimension (its:ite):: &
  3305. xmb_ske,xmb_ave,xmb_std,xmb_cur,xmbweight
  3306. real, dimension (its:ite):: &
  3307. pr_ske,pr_ave,pr_std,pr_cur
  3308. real, dimension (its:ite,jts:jte):: &
  3309. pr_gr,pr_w,pr_mc,pr_st,pr_as,pr_capma, &
  3310. pr_capme,pr_capmi
  3311. integer :: iedt, kk
  3312. !
  3313. character *(*), intent (in) :: &
  3314. name
  3315. !
  3316. integer :: itf, ktf
  3317. itf=MIN(ite,ide-1)
  3318. ktf=MIN(kte,kde-1)
  3319. tuning=0.
  3320. !
  3321. !
  3322. DO k=kts,ktf
  3323. do i=its,itf
  3324. outtem(i,k)=0.
  3325. outq(i,k)=0.
  3326. outqc(i,k)=0.
  3327. enddo
  3328. enddo
  3329. do i=its,itf
  3330. pre(i)=0.
  3331. xmb(i)=0.
  3332. xfac1(i)=1.
  3333. xmbweight(i)=1.
  3334. enddo
  3335. do i=its,itf
  3336. IF(ierr(i).eq.0)then
  3337. do n=(iens-1)*nx*nx2*maxens3+1,iens*nx*nx2*maxens3
  3338. if(pr_ens(i,j,n).le.0.)then
  3339. xf_ens(i,j,n)=0.
  3340. endif
  3341. enddo
  3342. endif
  3343. enddo
  3344. !
  3345. !--- calculate ensemble average mass fluxes
  3346. !
  3347. call massflx_stats(xf_ens,ensdim,nx2,nx,maxens3, &
  3348. xmb_ave,xmb_std,xmb_cur,xmb_ske,j,ierr,1, &
  3349. APR_GR,APR_W,APR_MC,APR_ST,APR_AS, &
  3350. APR_CAPMA,APR_CAPME,APR_CAPMI, &
  3351. pr_gr,pr_w,pr_mc,pr_st,pr_as, &
  3352. pr_capma,pr_capme,pr_capmi, &
  3353. ids,ide, jds,jde, kds,kde, &
  3354. ims,ime, jms,jme, kms,kme, &
  3355. its,ite, jts,jte, kts,kte )
  3356. call massflx_stats(pr_ens,ensdim,nx2,nx,maxens3, &
  3357. pr_ave,pr_std,pr_cur,pr_ske,j,ierr,2, &
  3358. APR_GR,APR_W,APR_MC,APR_ST,APR_AS, &
  3359. APR_CAPMA,APR_CAPME,APR_CAPMI, &
  3360. pr_gr,pr_w,pr_mc,pr_st,pr_as, &
  3361. pr_capma,pr_capme,pr_capmi, &
  3362. ids,ide, jds,jde, kds,kde, &
  3363. ims,ime, jms,jme, kms,kme, &
  3364. its,ite, jts,jte, kts,kte )
  3365. !
  3366. !-- now do feedback
  3367. !
  3368. ddtes=200.
  3369. ! if(name.eq.'shal')ddtes=200.
  3370. do i=its,itf
  3371. if(ierr(i).eq.0)then
  3372. if(xmb_ave(i).le.0.)then
  3373. ierr(i)=13
  3374. xmb_ave(i)=0.
  3375. endif
  3376. ! xmb(i)=max(0.,xmb_ave(i))
  3377. xmb(i)=max(.1*xmb_ave(i),xmb_ave(i)-tuning*xmb_std(i))
  3378. ! --- Now use proper count of how many closures were actually
  3379. ! used in cup_forcing_ens (including screening of some
  3380. ! closures over water) to properly normalize xmb
  3381. clos_wei=16./max(1.,closure_n(i))
  3382. if (xland1(i).lt.0.5)xmb(i)=xmb(i)*clos_wei
  3383. if(xmb(i).eq.0.)then
  3384. ierr(i)=19
  3385. endif
  3386. if(xmb(i).gt.100.)then
  3387. ierr(i)=19
  3388. endif
  3389. xfac1(i)=xmb(i)
  3390. endif
  3391. xfac1(i)=xmb_ave(i)
  3392. ENDDO
  3393. DO k=kts,ktf
  3394. do i=its,itf
  3395. dtt=0.
  3396. dtq=0.
  3397. dtqc=0.
  3398. dtpw=0.
  3399. IF(ierr(i).eq.0.and.k.le.ktop(i))then
  3400. do n=1,nx
  3401. dtt=dtt+dellat(i,k,n)
  3402. dtq=dtq+dellaq(i,k,n)
  3403. dtqc=dtqc+dellaqc(i,k,n)
  3404. dtpw=dtpw+pw(i,k,n)
  3405. enddo
  3406. outtes=dtt*XMB(I)*86400./float(nx)
  3407. IF((OUTTES.GT.2.*ddtes.and.k.gt.2))THEN
  3408. XMB(I)= 2.*ddtes/outtes * xmb(i)
  3409. outtes=1.*ddtes
  3410. endif
  3411. if (outtes .lt. -ddtes) then
  3412. XMB(I)= -ddtes/outtes * xmb(i)
  3413. outtes=-ddtes
  3414. endif
  3415. if (outtes .gt. .5*ddtes.and.k.le.2) then
  3416. XMB(I)= ddtes/outtes * xmb(i)
  3417. outtes=.5*ddtes
  3418. endif
  3419. OUTTEM(I,K)=XMB(I)*dtt/float(nx)
  3420. OUTQ(I,K)=XMB(I)*dtq/float(nx)
  3421. OUTQC(I,K)=XMB(I)*dtqc/float(nx)
  3422. PRE(I)=PRE(I)+XMB(I)*dtpw/float(nx)
  3423. endif
  3424. enddo
  3425. enddo
  3426. do i=its,itf
  3427. if(ierr(i).eq.0)then
  3428. prerate=pre(i)*3600.
  3429. if(prerate.lt.0.1)then
  3430. if(ierr2(i).gt.0.or.ierr3(i).gt.0)then
  3431. pre(i)=0.
  3432. ierr(i)=221
  3433. do k=kts,ktf
  3434. outtem(i,k)=0.
  3435. outq(i,k)=0.
  3436. outqc(i,k)=0.
  3437. enddo
  3438. do k=(iens-1)*nx*nx2*maxens3+1,iens*nx*nx2*maxens3
  3439. massfln(i,j,k)=0.
  3440. xf_ens(i,j,k)=0.
  3441. enddo
  3442. endif
  3443. endif
  3444. endif
  3445. ENDDO
  3446. do i=its,itf
  3447. if(ierr(i).eq.0)then
  3448. xfac1(i)=xmb(i)/xfac1(i)
  3449. do k=(iens-1)*nx*nx2*maxens3+1,iens*nx*nx2*maxens3
  3450. massfln(i,j,k)=massfln(i,j,k)*xfac1(i)
  3451. xf_ens(i,j,k)=xf_ens(i,j,k)*xfac1(i)
  3452. enddo
  3453. endif
  3454. ENDDO
  3455. if (l_flux) then
  3456. if (iens .eq. 1) then ! Only do deep convection mass fluxes
  3457. do k=kts,ktf
  3458. do i=its,itf
  3459. outcfu1(i,k)=0.
  3460. outcfd1(i,k)=0.
  3461. outdfu1(i,k)=0.
  3462. outefu1(i,k)=0.
  3463. outdfd1(i,k)=0.
  3464. outefd1(i,k)=0.
  3465. if (ierr(i) .eq. 0) then
  3466. do iedt=1,nx
  3467. do kk=1,nx2*maxens3
  3468. n=(iens-1)*nx*nx2*maxens3 + &
  3469. (iedt-1)*nx2*maxens3 + kk
  3470. outcfu1(i,k)=outcfu1(i,k)+cfu1_ens(i,k,iedt)*xf_ens(i,j,n)
  3471. outcfd1(i,k)=outcfd1(i,k)+cfd1_ens(i,k,iedt)*xf_ens(i,j,n)
  3472. outdfu1(i,k)=outdfu1(i,k)+dfu1_ens(i,k,iedt)*xf_ens(i,j,n)
  3473. outefu1(i,k)=outefu1(i,k)+efu1_ens(i,k,iedt)*xf_ens(i,j,n)
  3474. outdfd1(i,k)=outdfd1(i,k)+dfd1_ens(i,k,iedt)*xf_ens(i,j,n)
  3475. outefd1(i,k)=outefd1(i,k)+efd1_ens(i,k,iedt)*xf_ens(i,j,n)
  3476. end do
  3477. end do
  3478. outcfu1(i,k)=outcfu1(i,k)/ensdim
  3479. outcfd1(i,k)=outcfd1(i,k)/ensdim
  3480. outdfu1(i,k)=outdfu1(i,k)/ensdim
  3481. outefu1(i,k)=outefu1(i,k)/ensdim
  3482. outdfd1(i,k)=outdfd1(i,k)/ensdim
  3483. outefd1(i,k)=outefd1(i,k)/ensdim
  3484. end if !ierr
  3485. end do !i
  3486. end do !k
  3487. end if !iens .eq. 1
  3488. end if !l_flux
  3489. END SUBROUTINE cup_output_ens
  3490. SUBROUTINE cup_up_aa0(aa0,z,zu,dby,GAMMA_CUP,t_cup, &
  3491. kbcon,ktop,ierr, &
  3492. ids,ide, jds,jde, kds,kde, &
  3493. ims,ime, jms,jme, kms,kme, &
  3494. its,ite, jts,jte, kts,kte )
  3495. IMPLICIT NONE
  3496. !
  3497. ! on input
  3498. !
  3499. ! only local wrf dimensions are need as of now in this routine
  3500. integer &
  3501. ,intent (in ) :: &
  3502. ids,ide, jds,jde, kds,kde, &
  3503. ims,ime, jms,jme, kms,kme, &
  3504. its,ite, jts,jte, kts,kte
  3505. ! aa0 cloud work function
  3506. ! gamma_cup = gamma on model cloud levels
  3507. ! t_cup = temperature (Kelvin) on model cloud levels
  3508. ! dby = buoancy term
  3509. ! zu= normalized updraft mass flux
  3510. ! z = heights of model levels
  3511. ! ierr error value, maybe modified in this routine
  3512. !
  3513. real, dimension (its:ite,kts:kte) &
  3514. ,intent (in ) :: &
  3515. z,zu,gamma_cup,t_cup,dby
  3516. integer, dimension (its:ite) &
  3517. ,intent (in ) :: &
  3518. kbcon,ktop
  3519. !
  3520. ! input and output
  3521. !
  3522. integer, dimension (its:ite) &
  3523. ,intent (inout) :: &
  3524. ierr
  3525. real, dimension (its:ite) &
  3526. ,intent (out ) :: &
  3527. aa0
  3528. !
  3529. ! local variables in this routine
  3530. !
  3531. integer :: &
  3532. i,k
  3533. real :: &
  3534. dz,da
  3535. !
  3536. integer :: itf, ktf
  3537. itf = MIN(ite,ide-1)
  3538. ktf = MIN(kte,kde-1)
  3539. do i=its,itf
  3540. aa0(i)=0.
  3541. enddo
  3542. DO 100 k=kts+1,ktf
  3543. DO 100 i=its,itf
  3544. IF(ierr(i).ne.0)GO TO 100
  3545. IF(K.LE.KBCON(I))GO TO 100
  3546. IF(K.Gt.KTOP(I))GO TO 100
  3547. DZ=Z(I,K)-Z(I,K-1)
  3548. da=zu(i,k)*DZ*(9.81/(1004.*( &
  3549. (T_cup(I,K)))))*DBY(I,K-1)/ &
  3550. (1.+GAMMA_CUP(I,K))
  3551. IF(K.eq.KTOP(I).and.da.le.0.)go to 100
  3552. AA0(I)=AA0(I)+da
  3553. if(aa0(i).lt.0.)aa0(i)=0.
  3554. 100 continue
  3555. END SUBROUTINE cup_up_aa0
  3556. SUBROUTINE cup_up_he(k22,hkb,z_cup,cd,entr,he_cup,hc, &
  3557. kbcon,ierr,dby,he,hes_cup, &
  3558. ids,ide, jds,jde, kds,kde, &
  3559. ims,ime, jms,jme, kms,kme, &
  3560. its,ite, jts,jte, kts,kte )
  3561. IMPLICIT NONE
  3562. !
  3563. ! on input
  3564. !
  3565. ! only local wrf dimensions are need as of now in this routine
  3566. integer &
  3567. ,intent (in ) :: &
  3568. ids,ide, jds,jde, kds,kde, &
  3569. ims,ime, jms,jme, kms,kme, &
  3570. its,ite, jts,jte, kts,kte
  3571. ! hc = cloud moist static energy
  3572. ! hkb = moist static energy at originating level
  3573. ! he = moist static energy on model levels
  3574. ! he_cup = moist static energy on model cloud levels
  3575. ! hes_cup = saturation moist static energy on model cloud levels
  3576. ! dby = buoancy term
  3577. ! cd= detrainment function
  3578. ! z_cup = heights of model cloud levels
  3579. ! entr = entrainment rate
  3580. !
  3581. real, dimension (its:ite,kts:kte) &
  3582. ,intent (in ) :: &
  3583. he,he_cup,hes_cup,z_cup,cd
  3584. ! entr= entrainment rate
  3585. real &
  3586. ,intent (in ) :: &
  3587. entr
  3588. integer, dimension (its:ite) &
  3589. ,intent (in ) :: &
  3590. kbcon,k22
  3591. !
  3592. ! input and output
  3593. !
  3594. ! ierr error value, maybe modified in this routine
  3595. integer, dimension (its:ite) &
  3596. ,intent (inout) :: &
  3597. ierr
  3598. real, dimension (its:ite,kts:kte) &
  3599. ,intent (out ) :: &
  3600. hc,dby
  3601. real, dimension (its:ite) &
  3602. ,intent (out ) :: &
  3603. hkb
  3604. !
  3605. ! local variables in this routine
  3606. !
  3607. integer :: &
  3608. i,k
  3609. real :: &
  3610. dz
  3611. !
  3612. integer :: itf, ktf
  3613. itf = MIN(ite,ide-1)
  3614. ktf = MIN(kte,kde-1)
  3615. !
  3616. !--- moist static energy inside cloud
  3617. !
  3618. do i=its,itf
  3619. if(ierr(i).eq.0.)then
  3620. hkb(i)=he_cup(i,k22(i))
  3621. do k=1,k22(i)
  3622. hc(i,k)=he_cup(i,k)
  3623. DBY(I,K)=0.
  3624. enddo
  3625. do k=k22(i),kbcon(i)-1
  3626. hc(i,k)=hkb(i)
  3627. DBY(I,K)=0.
  3628. enddo
  3629. k=kbcon(i)
  3630. hc(i,k)=hkb(i)
  3631. DBY(I,Kbcon(i))=Hkb(I)-HES_cup(I,K)
  3632. endif
  3633. enddo
  3634. do k=kts+1,ktf
  3635. do i=its,itf
  3636. if(k.gt.kbcon(i).and.ierr(i).eq.0.)then
  3637. DZ=Z_cup(i,K)-Z_cup(i,K-1)
  3638. HC(i,K)=(HC(i,K-1)*(1.-.5*CD(i,K)*DZ)+entr* &
  3639. DZ*HE(i,K-1))/(1.+entr*DZ-.5*cd(i,k)*dz)
  3640. DBY(I,K)=HC(I,K)-HES_cup(I,K)
  3641. endif
  3642. enddo
  3643. enddo
  3644. END SUBROUTINE cup_up_he
  3645. SUBROUTINE cup_up_moisture(ierr,z_cup,qc,qrc,pw,pwav, &
  3646. kbcon,ktop,cd,dby,mentr_rate,clw_all, &
  3647. q,GAMMA_cup,zu,qes_cup,k22,qe_cup,xl, &
  3648. ids,ide, jds,jde, kds,kde, &
  3649. ims,ime, jms,jme, kms,kme, &
  3650. its,ite, jts,jte, kts,kte )
  3651. IMPLICIT NONE
  3652. !
  3653. ! on input
  3654. !
  3655. ! only local wrf dimensions are need as of now in this routine
  3656. integer &
  3657. ,intent (in ) :: &
  3658. ids,ide, jds,jde, kds,kde, &
  3659. ims,ime, jms,jme, kms,kme, &
  3660. its,ite, jts,jte, kts,kte
  3661. ! cd= detrainment function
  3662. ! q = environmental q on model levels
  3663. ! qe_cup = environmental q on model cloud levels
  3664. ! qes_cup = saturation q on model cloud levels
  3665. ! dby = buoancy term
  3666. ! cd= detrainment function
  3667. ! zu = normalized updraft mass flux
  3668. ! gamma_cup = gamma on model cloud levels
  3669. ! mentr_rate = entrainment rate
  3670. !
  3671. real, dimension (its:ite,kts:kte) &
  3672. ,intent (in ) :: &
  3673. q,zu,gamma_cup,qe_cup,dby,qes_cup,z_cup,cd
  3674. ! entr= entrainment rate
  3675. real &
  3676. ,intent (in ) :: &
  3677. mentr_rate,xl
  3678. integer, dimension (its:ite) &
  3679. ,intent (in ) :: &
  3680. kbcon,ktop,k22
  3681. !
  3682. ! input and output
  3683. !
  3684. ! ierr error value, maybe modified in this routine
  3685. integer, dimension (its:ite) &
  3686. ,intent (inout) :: &
  3687. ierr
  3688. ! qc = cloud q (including liquid water) after entrainment
  3689. ! qrch = saturation q in cloud
  3690. ! qrc = liquid water content in cloud after rainout
  3691. ! pw = condensate that will fall out at that level
  3692. ! pwav = totan normalized integrated condensate (I1)
  3693. ! c0 = conversion rate (cloud to rain)
  3694. real, dimension (its:ite,kts:kte) &
  3695. ,intent (out ) :: &
  3696. qc,qrc,pw,clw_all
  3697. real, dimension (its:ite) &
  3698. ,intent (out ) :: &
  3699. pwav
  3700. !
  3701. ! local variables in this routine
  3702. !
  3703. integer :: &
  3704. iall,i,k
  3705. real :: &
  3706. dh,qrch,c0,dz,radius
  3707. !
  3708. integer :: itf, ktf
  3709. itf = MIN(ite,ide-1)
  3710. ktf = MIN(kte,kde-1)
  3711. iall=0
  3712. c0=.002
  3713. !
  3714. !--- no precip for small clouds
  3715. !
  3716. if(mentr_rate.gt.0.)then
  3717. radius=.2/mentr_rate
  3718. if(radius.lt.900.)c0=0.
  3719. ! if(radius.lt.900.)iall=0
  3720. endif
  3721. do i=its,itf
  3722. pwav(i)=0.
  3723. enddo
  3724. do k=kts,ktf
  3725. do i=its,itf
  3726. pw(i,k)=0.
  3727. if(ierr(i).eq.0)qc(i,k)=qes_cup(i,k)
  3728. clw_all(i,k)=0.
  3729. qrc(i,k)=0.
  3730. enddo
  3731. enddo
  3732. do i=its,itf
  3733. if(ierr(i).eq.0.)then
  3734. do k=k22(i),kbcon(i)-1
  3735. qc(i,k)=qe_cup(i,k22(i))
  3736. enddo
  3737. endif
  3738. enddo
  3739. DO 100 k=kts+1,ktf
  3740. DO 100 i=its,itf
  3741. IF(ierr(i).ne.0)GO TO 100
  3742. IF(K.Lt.KBCON(I))GO TO 100
  3743. IF(K.Gt.KTOP(I))GO TO 100
  3744. DZ=Z_cup(i,K)-Z_cup(i,K-1)
  3745. !
  3746. !------ 1. steady state plume equation, for what could
  3747. !------ be in cloud without condensation
  3748. !
  3749. !
  3750. QC(i,K)=(QC(i,K-1)*(1.-.5*CD(i,K)*DZ)+mentr_rate* &
  3751. DZ*Q(i,K-1))/(1.+mentr_rate*DZ-.5*cd(i,k)*dz)
  3752. !
  3753. !--- saturation in cloud, this is what is allowed to be in it
  3754. !
  3755. QRCH=QES_cup(I,K)+(1./XL)*(GAMMA_cup(i,k) &
  3756. /(1.+GAMMA_cup(i,k)))*DBY(I,K)
  3757. !
  3758. !------- LIQUID WATER CONTENT IN cloud after rainout
  3759. !
  3760. clw_all(i,k)=QC(I,K)-QRCH
  3761. QRC(I,K)=(QC(I,K)-QRCH)/(1.+C0*DZ*zu(i,k))
  3762. if(qrc(i,k).lt.0.)then
  3763. qrc(i,k)=0.
  3764. endif
  3765. !
  3766. !------- 3.Condensation
  3767. !
  3768. PW(i,k)=c0*dz*QRC(I,K)*zu(i,k)
  3769. if(iall.eq.1)then
  3770. qrc(i,k)=0.
  3771. pw(i,k)=(QC(I,K)-QRCH)*zu(i,k)
  3772. if(pw(i,k).lt.0.)pw(i,k)=0.
  3773. endif
  3774. !
  3775. !----- set next level
  3776. !
  3777. QC(I,K)=QRC(I,K)+qrch
  3778. !
  3779. !--- integrated normalized ondensate
  3780. !
  3781. PWAV(I)=PWAV(I)+PW(I,K)
  3782. 100 CONTINUE
  3783. END SUBROUTINE cup_up_moisture
  3784. SUBROUTINE cup_up_nms(zu,z_cup,entr,cd,kbcon,ktop,ierr,k22, &
  3785. ids,ide, jds,jde, kds,kde, &
  3786. ims,ime, jms,jme, kms,kme, &
  3787. its,ite, jts,jte, kts,kte )
  3788. IMPLICIT NONE
  3789. !
  3790. ! on input
  3791. !
  3792. ! only local wrf dimensions are need as of now in this routine
  3793. integer &
  3794. ,intent (in ) :: &
  3795. ids,ide, jds,jde, kds,kde, &
  3796. ims,ime, jms,jme, kms,kme, &
  3797. its,ite, jts,jte, kts,kte
  3798. ! cd= detrainment function
  3799. real, dimension (its:ite,kts:kte) &
  3800. ,intent (in ) :: &
  3801. z_cup,cd
  3802. ! entr= entrainment rate
  3803. real &
  3804. ,intent (in ) :: &
  3805. entr
  3806. integer, dimension (its:ite) &
  3807. ,intent (in ) :: &
  3808. kbcon,ktop,k22
  3809. !
  3810. ! input and output
  3811. !
  3812. ! ierr error value, maybe modified in this routine
  3813. integer, dimension (its:ite) &
  3814. ,intent (inout) :: &
  3815. ierr
  3816. ! zu is the normalized mass flux
  3817. real, dimension (its:ite,kts:kte) &
  3818. ,intent (out ) :: &
  3819. zu
  3820. !
  3821. ! local variables in this routine
  3822. !
  3823. integer :: &
  3824. i,k
  3825. real :: &
  3826. dz
  3827. integer :: itf, ktf
  3828. itf = MIN(ite,ide-1)
  3829. ktf = MIN(kte,kde-1)
  3830. !
  3831. ! initialize for this go around
  3832. !
  3833. do k=kts,ktf
  3834. do i=its,itf
  3835. zu(i,k)=0.
  3836. enddo
  3837. enddo
  3838. !
  3839. ! do normalized mass budget
  3840. !
  3841. do i=its,itf
  3842. IF(ierr(I).eq.0)then
  3843. do k=k22(i),kbcon(i)
  3844. zu(i,k)=1.
  3845. enddo
  3846. DO K=KBcon(i)+1,KTOP(i)
  3847. DZ=Z_cup(i,K)-Z_cup(i,K-1)
  3848. ZU(i,K)=ZU(i,K-1)*(1.+(entr-cd(i,k))*DZ)
  3849. enddo
  3850. endif
  3851. enddo
  3852. END SUBROUTINE cup_up_nms
  3853. !====================================================================
  3854. SUBROUTINE gdinit(RTHCUTEN,RQVCUTEN,RQCCUTEN,RQICUTEN, &
  3855. MASS_FLUX,cp,restart, &
  3856. P_QC,P_QI,P_FIRST_SCALAR, &
  3857. RTHFTEN, RQVFTEN, &
  3858. APR_GR,APR_W,APR_MC,APR_ST,APR_AS, &
  3859. APR_CAPMA,APR_CAPME,APR_CAPMI, &
  3860. allowed_to_read, &
  3861. ids, ide, jds, jde, kds, kde, &
  3862. ims, ime, jms, jme, kms, kme, &
  3863. its, ite, jts, jte, kts, kte )
  3864. !--------------------------------------------------------------------
  3865. IMPLICIT NONE
  3866. !--------------------------------------------------------------------
  3867. LOGICAL , INTENT(IN) :: restart,allowed_to_read
  3868. INTEGER , INTENT(IN) :: ids, ide, jds, jde, kds, kde, &
  3869. ims, ime, jms, jme, kms, kme, &
  3870. its, ite, jts, jte, kts, kte
  3871. INTEGER , INTENT(IN) :: P_FIRST_SCALAR, P_QI, P_QC
  3872. REAL, INTENT(IN) :: cp
  3873. REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(OUT) :: &
  3874. RTHCUTEN, &
  3875. RQVCUTEN, &
  3876. RQCCUTEN, &
  3877. RQICUTEN
  3878. REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(OUT) :: &
  3879. RTHFTEN, &
  3880. RQVFTEN
  3881. REAL, DIMENSION( ims:ime , jms:jme ) , INTENT(OUT) :: &
  3882. APR_GR,APR_W,APR_MC,APR_ST,APR_AS, &
  3883. APR_CAPMA,APR_CAPME,APR_CAPMI, &
  3884. MASS_FLUX
  3885. INTEGER :: i, j, k, itf, jtf, ktf
  3886. jtf=min0(jte,jde-1)
  3887. ktf=min0(kte,kde-1)
  3888. itf=min0(ite,ide-1)
  3889. IF(.not.restart)THEN
  3890. DO j=jts,jtf
  3891. DO k=kts,ktf
  3892. DO i=its,itf
  3893. RTHCUTEN(i,k,j)=0.
  3894. RQVCUTEN(i,k,j)=0.
  3895. ENDDO
  3896. ENDDO
  3897. ENDDO
  3898. DO j=jts,jtf
  3899. DO k=kts,ktf
  3900. DO i=its,itf
  3901. RTHFTEN(i,k,j)=0.
  3902. RQVFTEN(i,k,j)=0.
  3903. ENDDO
  3904. ENDDO
  3905. ENDDO
  3906. IF (P_QC .ge. P_FIRST_SCALAR) THEN
  3907. DO j=jts,jtf
  3908. DO k=kts,ktf
  3909. DO i=its,itf
  3910. RQCCUTEN(i,k,j)=0.
  3911. ENDDO
  3912. ENDDO
  3913. ENDDO
  3914. ENDIF
  3915. IF (P_QI .ge. P_FIRST_SCALAR) THEN
  3916. DO j=jts,jtf
  3917. DO k=kts,ktf
  3918. DO i=its,itf
  3919. RQICUTEN(i,k,j)=0.
  3920. ENDDO
  3921. ENDDO
  3922. ENDDO
  3923. ENDIF
  3924. DO j=jts,jtf
  3925. DO i=its,itf
  3926. mass_flux(i,j)=0.
  3927. ENDDO
  3928. ENDDO
  3929. ENDIF
  3930. DO j=jts,jtf
  3931. DO i=its,itf
  3932. APR_GR(i,j)=0.
  3933. APR_ST(i,j)=0.
  3934. APR_W(i,j)=0.
  3935. APR_MC(i,j)=0.
  3936. APR_AS(i,j)=0.
  3937. APR_CAPMA(i,j)=0.
  3938. APR_CAPME(i,j)=0.
  3939. APR_CAPMI(i,j)=0.
  3940. ENDDO
  3941. ENDDO
  3942. END SUBROUTINE gdinit
  3943. SUBROUTINE massflx_stats(xf_ens,ensdim,maxens,maxens2,maxens3, &
  3944. xt_ave,xt_std,xt_cur,xt_ske,j,ierr,itest, &
  3945. APR_GR,APR_W,APR_MC,APR_ST,APR_AS, &
  3946. APR_CAPMA,APR_CAPME,APR_CAPMI, &
  3947. pr_gr,pr_w,pr_mc,pr_st,pr_as, &
  3948. pr_capma,pr_capme,pr_capmi, &
  3949. ids,ide, jds,jde, kds,kde, &
  3950. ims,ime, jms,jme, kms,kme, &
  3951. its,ite, jts,jte, kts,kte)
  3952. IMPLICIT NONE
  3953. integer, intent (in ) :: &
  3954. j,ensdim,maxens3,maxens,maxens2,itest
  3955. INTEGER, INTENT(IN ) :: &
  3956. ids,ide, jds,jde, kds,kde, &
  3957. ims,ime, jms,jme, kms,kme, &
  3958. its,ite, jts,jte, kts,kte
  3959. real, dimension (its:ite) &
  3960. , intent(inout) :: &
  3961. xt_ave,xt_cur,xt_std,xt_ske
  3962. integer, dimension (its:ite), intent (in) :: &
  3963. ierr
  3964. real, dimension (ims:ime,jms:jme,1:ensdim) &
  3965. , intent(in ) :: &
  3966. xf_ens
  3967. real, dimension (ims:ime,jms:jme) &
  3968. , intent(inout) :: &
  3969. APR_GR,APR_W,APR_MC,APR_ST,APR_AS, &
  3970. APR_CAPMA,APR_CAPME,APR_CAPMI
  3971. real, dimension (its:ite,jts:jte) &
  3972. , intent(inout) :: &
  3973. pr_gr,pr_w,pr_mc,pr_st,pr_as, &
  3974. pr_capma,pr_capme,pr_capmi
  3975. !
  3976. ! local stuff
  3977. !
  3978. real, dimension (its:ite , 1:maxens3 ) :: &
  3979. x_ave,x_cur,x_std,x_ske
  3980. real, dimension (its:ite , 1:maxens ) :: &
  3981. x_ave_cap
  3982. integer, dimension (1:maxens3) :: nc1
  3983. integer :: i,k
  3984. integer :: num,kk,num2,iedt
  3985. real :: a3,a4
  3986. num=ensdim/maxens3
  3987. num2=ensdim/maxens
  3988. if(itest.eq.1)then
  3989. do i=its,ite
  3990. pr_gr(i,j) = 0.
  3991. pr_w(i,j) = 0.
  3992. pr_mc(i,j) = 0.
  3993. pr_st(i,j) = 0.
  3994. pr_as(i,j) = 0.
  3995. pr_capma(i,j) = 0.
  3996. pr_capme(i,j) = 0.
  3997. pr_capmi(i,j) = 0.
  3998. enddo
  3999. endif
  4000. do k=1,maxens
  4001. do i=its,ite
  4002. x_ave_cap(i,k)=0.
  4003. enddo
  4004. enddo
  4005. do k=1,maxens3
  4006. do i=its,ite
  4007. x_ave(i,k)=0.
  4008. x_std(i,k)=0.
  4009. x_ske(i,k)=0.
  4010. x_cur(i,k)=0.
  4011. enddo
  4012. enddo
  4013. do i=its,ite
  4014. xt_ave(i)=0.
  4015. xt_std(i)=0.
  4016. xt_ske(i)=0.
  4017. xt_cur(i)=0.
  4018. enddo
  4019. do kk=1,num
  4020. do k=1,maxens3
  4021. do i=its,ite
  4022. if(ierr(i).eq.0)then
  4023. x_ave(i,k)=x_ave(i,k)+xf_ens(i,j,maxens3*(kk-1)+k)
  4024. endif
  4025. enddo
  4026. enddo
  4027. enddo
  4028. do iedt=1,maxens2
  4029. do k=1,maxens
  4030. do kk=1,maxens3
  4031. do i=its,ite
  4032. if(ierr(i).eq.0)then
  4033. x_ave_cap(i,k)=x_ave_cap(i,k) &
  4034. +xf_ens(i,j,maxens3*(k-1)+(iedt-1)*maxens*maxens3+kk)
  4035. endif
  4036. enddo
  4037. enddo
  4038. enddo
  4039. enddo
  4040. do k=1,maxens
  4041. do i=its,ite
  4042. if(ierr(i).eq.0)then
  4043. x_ave_cap(i,k)=x_ave_cap(i,k)/float(num2)
  4044. endif
  4045. enddo
  4046. enddo
  4047. do k=1,maxens3
  4048. do i=its,ite
  4049. if(ierr(i).eq.0)then
  4050. x_ave(i,k)=x_ave(i,k)/float(num)
  4051. endif
  4052. enddo
  4053. enddo
  4054. do k=1,maxens3
  4055. do i=its,ite
  4056. if(ierr(i).eq.0)then
  4057. xt_ave(i)=xt_ave(i)+x_ave(i,k)
  4058. endif
  4059. enddo
  4060. enddo
  4061. do i=its,ite
  4062. if(ierr(i).eq.0)then
  4063. xt_ave(i)=xt_ave(i)/float(maxens3)
  4064. endif
  4065. enddo
  4066. !
  4067. !--- now do std, skewness,curtosis
  4068. !
  4069. do kk=1,num
  4070. do k=1,maxens3
  4071. do i=its,ite
  4072. if(ierr(i).eq.0.and.x_ave(i,k).gt.0.)then
  4073. ! print *,i,j,k,kk,x_std(i,k),xf_ens(i,j,maxens3*(kk-1)+k),x_ave(i,k)
  4074. x_std(i,k)=x_std(i,k)+(xf_ens(i,j,maxens3*(kk-1)+k)-x_ave(i,k))**2
  4075. x_ske(i,k)=x_ske(i,k)+(xf_ens(i,j,maxens3*(kk-1)+k)-x_ave(i,k))**3
  4076. x_cur(i,k)=x_cur(i,k)+(xf_ens(i,j,maxens3*(kk-1)+k)-x_ave(i,k))**4
  4077. endif
  4078. enddo
  4079. enddo
  4080. enddo
  4081. do k=1,maxens3
  4082. do i=its,ite
  4083. if(ierr(i).eq.0.and.xt_ave(i).gt.0.)then
  4084. xt_std(i)=xt_std(i)+(x_ave(i,k)-xt_ave(i))**2
  4085. xt_ske(i)=xt_ske(i)+(x_ave(i,k)-xt_ave(i))**3
  4086. xt_cur(i)=xt_cur(i)+(x_ave(i,k)-xt_ave(i))**4
  4087. endif
  4088. enddo
  4089. enddo
  4090. do k=1,maxens3
  4091. do i=its,ite
  4092. if(ierr(i).eq.0.and.x_std(i,k).gt.0.)then
  4093. x_std(i,k)=x_std(i,k)/float(num)
  4094. a3=max(1.e-6,x_std(i,k))
  4095. x_std(i,k)=sqrt(a3)
  4096. a3=max(1.e-6,x_std(i,k)**3)
  4097. a4=max(1.e-6,x_std(i,k)**4)
  4098. x_ske(i,k)=x_ske(i,k)/float(num)/a3
  4099. x_cur(i,k)=x_cur(i,k)/float(num)/a4
  4100. endif
  4101. ! print*,' '
  4102. ! print*,'Some statistics at gridpoint i,j, ierr',i,j,ierr(i)
  4103. ! print*,'statistics for closure number ',k
  4104. ! print*,'Average= ',x_ave(i,k),' Std= ',x_std(i,k)
  4105. ! print*,'Skewness= ',x_ske(i,k),' Curtosis= ',x_cur(i,k)
  4106. ! print*,' '
  4107. enddo
  4108. enddo
  4109. do i=its,ite
  4110. if(ierr(i).eq.0.and.xt_std(i).gt.0.)then
  4111. xt_std(i)=xt_std(i)/float(maxens3)
  4112. a3=max(1.e-6,xt_std(i))
  4113. xt_std(i)=sqrt(a3)
  4114. a3=max(1.e-6,xt_std(i)**3)
  4115. a4=max(1.e-6,xt_std(i)**4)
  4116. xt_ske(i)=xt_ske(i)/float(maxens3)/a3
  4117. xt_cur(i)=xt_cur(i)/float(maxens3)/a4
  4118. ! print*,' '
  4119. ! print*,'Total ensemble independent statistics at i =',i
  4120. ! print*,'Average= ',xt_ave(i),' Std= ',xt_std(i)
  4121. ! print*,'Skewness= ',xt_ske(i),' Curtosis= ',xt_cur(i)
  4122. ! print*,' '
  4123. !
  4124. ! first go around: store massflx for different closures/caps
  4125. !
  4126. if(itest.eq.1)then
  4127. pr_gr(i,j) = .333*(x_ave(i,1)+x_ave(i,2)+x_ave(i,3))
  4128. pr_w(i,j) = .333*(x_ave(i,4)+x_ave(i,5)+x_ave(i,6))
  4129. pr_mc(i,j) = .333*(x_ave(i,7)+x_ave(i,8)+x_ave(i,9))
  4130. pr_st(i,j) = .333*(x_ave(i,10)+x_ave(i,11)+x_ave(i,12))
  4131. pr_as(i,j) = .25*(x_ave(i,13)+x_ave(i,14)+x_ave(i,15) &
  4132. + x_ave(i,16))
  4133. pr_capma(i,j) = x_ave_cap(i,1)
  4134. pr_capme(i,j) = x_ave_cap(i,2)
  4135. pr_capmi(i,j) = x_ave_cap(i,3)
  4136. !
  4137. ! second go around: store preciprates (mm/hour) for different closures/caps
  4138. !
  4139. else if (itest.eq.2)then
  4140. APR_GR(i,j)=.333*(x_ave(i,1)+x_ave(i,2)+x_ave(i,3))* &
  4141. 3600.*pr_gr(i,j) +APR_GR(i,j)
  4142. APR_W(i,j)=.333*(x_ave(i,4)+x_ave(i,5)+x_ave(i,6))* &
  4143. 3600.*pr_w(i,j) +APR_W(i,j)
  4144. APR_MC(i,j)=.333*(x_ave(i,7)+x_ave(i,8)+x_ave(i,9))* &
  4145. 3600.*pr_mc(i,j) +APR_MC(i,j)
  4146. APR_ST(i,j)=.333*(x_ave(i,10)+x_ave(i,11)+x_ave(i,12))* &
  4147. 3600.*pr_st(i,j) +APR_ST(i,j)
  4148. APR_AS(i,j)=.25*(x_ave(i,13)+x_ave(i,14)+x_ave(i,15) &
  4149. + x_ave(i,16))* &
  4150. 3600.*pr_as(i,j) +APR_AS(i,j)
  4151. APR_CAPMA(i,j) = x_ave_cap(i,1)* &
  4152. 3600.*pr_capma(i,j) +APR_CAPMA(i,j)
  4153. APR_CAPME(i,j) = x_ave_cap(i,2)* &
  4154. 3600.*pr_capme(i,j) +APR_CAPME(i,j)
  4155. APR_CAPMI(i,j) = x_ave_cap(i,3)* &
  4156. 3600.*pr_capmi(i,j) +APR_CAPMI(i,j)
  4157. endif
  4158. endif
  4159. enddo
  4160. END SUBROUTINE massflx_stats
  4161. SUBROUTINE neg_check(dt,q,outq,outt,outqc,pret,its,ite,kts,kte,itf,ktf)
  4162. INTEGER, INTENT(IN ) :: its,ite,kts,kte,itf,ktf
  4163. real, dimension (its:ite,kts:kte ) , &
  4164. intent(inout ) :: &
  4165. q,outq,outt,outqc
  4166. real, dimension (its:ite ) , &
  4167. intent(inout ) :: &
  4168. pret
  4169. real &
  4170. ,intent (in ) :: &
  4171. dt
  4172. real :: thresh,qmem,qmemf,qmem2,qtest,qmem1
  4173. !
  4174. ! first do check on vertical heating rate
  4175. !
  4176. thresh=200.01
  4177. do i=its,itf
  4178. qmemf=1.
  4179. qmem=0.
  4180. do k=kts,ktf
  4181. qmem=outt(i,k)*86400.
  4182. if(qmem.gt.2.*thresh)then
  4183. qmem2=2.*thresh/qmem
  4184. qmemf=min(qmemf,qmem2)
  4185. !
  4186. !
  4187. ! print *,'1',' adjusted massflux by factor ',i,k,qmem,qmem2,qmemf
  4188. endif
  4189. if(qmem.lt.-thresh)then
  4190. qmem2=-thresh/qmem
  4191. qmemf=min(qmemf,qmem2)
  4192. !
  4193. !
  4194. ! print *,'2',' adjusted massflux by factor ',i,k,qmem,qmem2,qmemf
  4195. endif
  4196. enddo
  4197. do k=kts,ktf
  4198. outq(i,k)=outq(i,k)*qmemf
  4199. outt(i,k)=outt(i,k)*qmemf
  4200. outqc(i,k)=outqc(i,k)*qmemf
  4201. enddo
  4202. pret(i)=pret(i)*qmemf
  4203. enddo
  4204. !
  4205. ! check whether routine produces negative q's. This can happen, since
  4206. ! tendencies are calculated based on forced q's. This should have no
  4207. ! influence on conservation properties, it scales linear through all
  4208. ! tendencies
  4209. !
  4210. thresh=1.e-10
  4211. do i=its,itf
  4212. qmemf=1.
  4213. do k=kts,ktf
  4214. qmem=outq(i,k)
  4215. if(abs(qmem).gt.0.)then
  4216. qtest=q(i,k)+outq(i,k)*dt
  4217. if(qtest.lt.thresh)then
  4218. !
  4219. ! qmem2 would be the maximum allowable tendency
  4220. !
  4221. qmem1=outq(i,k)
  4222. qmem2=(thresh-q(i,k))/dt
  4223. qmemf=min(qmemf,qmem2/qmem1)
  4224. ! qmemf=max(0.,qmemf)
  4225. ! print *,'4 adjusted tendencies ',i,k,qmem,qmem2,qmemf
  4226. endif
  4227. endif
  4228. enddo
  4229. do k=kts,ktf
  4230. outq(i,k)=outq(i,k)*qmemf
  4231. outt(i,k)=outt(i,k)*qmemf
  4232. outqc(i,k)=outqc(i,k)*qmemf
  4233. enddo
  4234. pret(i)=pret(i)*qmemf
  4235. enddo
  4236. END SUBROUTINE neg_check
  4237. !-------------------------------------------------------
  4238. END MODULE module_cu_gd