PageRenderTime 95ms CodeModel.GetById 16ms RepoModel.GetById 1ms app.codeStats 1ms

/wrfv2_fire/phys/module_cu_sas.F

http://github.com/jbeezley/wrf-fire
FORTRAN Legacy | 5568 lines | 4146 code | 107 blank | 1315 comment | 325 complexity | a86d7ae44bfb23c8b96a4392f35319e5 MD5 | raw file
Possible License(s): AGPL-1.0
  1. !!
  2. MODULE module_cu_sas
  3. CONTAINS
  4. !-----------------------------------------------------------------
  5. SUBROUTINE CU_SAS(DT,ITIMESTEP,STEPCU, &
  6. RTHCUTEN,RQVCUTEN,RQCCUTEN,RQICUTEN, &
  7. RUCUTEN,RVCUTEN, &
  8. RAINCV,PRATEC,HTOP,HBOT, &
  9. U3D,V3D,W,T3D,QV3D,QC3D,QI3D,PI3D,RHO3D, &
  10. DZ8W,PCPS,P8W,XLAND,CU_ACT_FLAG, &
  11. P_QC, &
  12. MOMMIX, & ! gopal's doing
  13. PGCON,sas_mass_flux, &
  14. shalconv,shal_pgcon, &
  15. HPBL2D,EVAP2D,HEAT2D, & !Kwon for shallow convection
  16. P_QI,P_FIRST_SCALAR, &
  17. CUDT, CURR_SECS, ADAPT_STEP_FLAG, &
  18. CUDTACTTIME, &
  19. ids,ide, jds,jde, kds,kde, &
  20. ims,ime, jms,jme, kms,kme, &
  21. its,ite, jts,jte, kts,kte )
  22. !-------------------------------------------------------------------
  23. USE MODULE_GFS_MACHINE , ONLY : kind_phys
  24. USE MODULE_GFS_FUNCPHYS , ONLY : gfuncphys
  25. USE MODULE_GFS_PHYSCONS, grav => con_g, CP => con_CP, HVAP => con_HVAP &
  26. &, RV => con_RV, FV => con_fvirt, T0C => con_T0C &
  27. &, CVAP => con_CVAP, CLIQ => con_CLIQ &
  28. &, EPS => con_eps, EPSM1 => con_epsm1 &
  29. &, ROVCP => con_rocp, RD => con_rd
  30. !-------------------------------------------------------------------
  31. IMPLICIT NONE
  32. !-------------------------------------------------------------------
  33. !-- U3D 3D u-velocity interpolated to theta points (m/s)
  34. !-- V3D 3D v-velocity interpolated to theta points (m/s)
  35. !-- TH3D 3D potential temperature (K)
  36. !-- T3D temperature (K)
  37. !-- QV3D 3D water vapor mixing ratio (Kg/Kg)
  38. !-- QC3D 3D cloud mixing ratio (Kg/Kg)
  39. !-- QI3D 3D ice mixing ratio (Kg/Kg)
  40. !-- P8w 3D pressure at full levels (Pa)
  41. !-- Pcps 3D pressure (Pa)
  42. !-- PI3D 3D exner function (dimensionless)
  43. !-- rr3D 3D dry air density (kg/m^3)
  44. !-- RUBLTEN U tendency due to
  45. ! PBL parameterization (m/s^2)
  46. !-- RVBLTEN V tendency due to
  47. ! PBL parameterization (m/s^2)
  48. !-- RTHBLTEN Theta tendency due to
  49. ! PBL parameterization (K/s)
  50. !-- RQVBLTEN Qv tendency due to
  51. ! PBL parameterization (kg/kg/s)
  52. !-- RQCBLTEN Qc tendency due to
  53. ! PBL parameterization (kg/kg/s)
  54. !-- RQIBLTEN Qi tendency due to
  55. ! PBL parameterization (kg/kg/s)
  56. !
  57. !-- MOMMIX MOMENTUM MIXING COEFFICIENT (can be set in the namelist)
  58. !-- RUCUTEN U tendency due to Cumulus Momentum Mixing (gopal's doing for SAS)
  59. !-- RVCUTEN V tendency due to Cumulus Momentum Mixing (gopal's doing for SAS)
  60. !
  61. !-- CP heat capacity at constant pressure for dry air (J/kg/K)
  62. !-- GRAV acceleration due to gravity (m/s^2)
  63. !-- ROVCP R/CP
  64. !-- RD gas constant for dry air (J/kg/K)
  65. !-- ROVG R/G
  66. !-- P_QI species index for cloud ice
  67. !-- dz8w dz between full levels (m)
  68. !-- z height above sea level (m)
  69. !-- PSFC pressure at the surface (Pa)
  70. !-- UST u* in similarity theory (m/s)
  71. !-- PBL PBL height (m)
  72. !-- PSIM similarity stability function for momentum
  73. !-- PSIH similarity stability function for heat
  74. !-- HFX upward heat flux at the surface (W/m^2)
  75. !-- QFX upward moisture flux at the surface (kg/m^2/s)
  76. !-- TSK surface temperature (K)
  77. !-- GZ1OZ0 log(z/z0) where z0 is roughness length
  78. !-- WSPD wind speed at lowest model level (m/s)
  79. !-- BR bulk Richardson number in surface layer
  80. !-- DT time step (s)
  81. !-- rvovrd R_v divided by R_d (dimensionless)
  82. !-- EP1 constant for virtual temperature (R_v/R_d - 1) (dimensionless)
  83. !-- KARMAN Von Karman constant
  84. !-- ids start index for i in domain
  85. !-- ide end index for i in domain
  86. !-- jds start index for j in domain
  87. !-- jde end index for j in domain
  88. !-- kds start index for k in domain
  89. !-- kde end index for k in domain
  90. !-- ims start index for i in memory
  91. !-- ime end index for i in memory
  92. !-- jms start index for j in memory
  93. !-- jme end index for j in memory
  94. !-- kms start index for k in memory
  95. !-- kme end index for k in memory
  96. !-- its start index for i in tile
  97. !-- ite end index for i in tile
  98. !-- jts start index for j in tile
  99. !-- jte end index for j in tile
  100. !-- kts start index for k in tile
  101. !-- kte end index for k in tile
  102. !-------------------------------------------------------------------
  103. INTEGER :: ICLDCK
  104. INTEGER, INTENT(IN) :: ids,ide, jds,jde, kds,kde, &
  105. ims,ime, jms,jme, kms,kme, &
  106. its,ite, jts,jte, kts,kte, &
  107. ITIMESTEP, & !NSTD
  108. P_FIRST_SCALAR, &
  109. P_QC, &
  110. P_QI, &
  111. STEPCU
  112. REAL, INTENT(IN) :: &
  113. DT
  114. REAL, OPTIONAL, INTENT(IN) :: PGCON,sas_mass_flux,shal_pgcon
  115. INTEGER, OPTIONAL, INTENT(IN) :: shalconv
  116. REAL(kind=kind_phys) :: PGCON_USE,SHAL_PGCON_USE,massf
  117. INTEGER :: shalconv_use
  118. REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: &
  119. RQCCUTEN, &
  120. RQICUTEN, &
  121. RQVCUTEN, &
  122. RTHCUTEN
  123. REAL, DIMENSION(ims:ime, jms:jme, kms:kme), INTENT(INOUT) :: &
  124. RUCUTEN, &
  125. RVCUTEN
  126. REAL, OPTIONAL, INTENT(IN) :: MOMMIX
  127. REAL, DIMENSION( ims:ime , jms:jme ), OPTIONAL, &
  128. INTENT(IN) :: HPBL2D,EVAP2D,HEAT2D !Kwon for sha
  129. REAL, DIMENSION(ims:ime, jms:jme), INTENT(IN) :: &
  130. XLAND
  131. REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: &
  132. RAINCV, PRATEC
  133. REAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT) :: &
  134. HBOT, &
  135. HTOP
  136. LOGICAL, DIMENSION(IMS:IME,JMS:JME), INTENT(INOUT) :: &
  137. CU_ACT_FLAG
  138. REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(IN) :: &
  139. DZ8W, &
  140. P8w, &
  141. Pcps, &
  142. PI3D, &
  143. QC3D, &
  144. QI3D, &
  145. QV3D, &
  146. RHO3D, &
  147. T3D, &
  148. U3D, &
  149. V3D, &
  150. W
  151. ! Adaptive time-step variables
  152. REAL, INTENT(IN ) :: CUDT
  153. REAL, INTENT(IN ) :: CURR_SECS
  154. LOGICAL,INTENT(IN ) , OPTIONAL :: ADAPT_STEP_FLAG
  155. REAL, INTENT (INOUT) :: CUDTACTTIME
  156. !--------------------------- LOCAL VARS ------------------------------
  157. REAL, DIMENSION(ims:ime, jms:jme) :: &
  158. PSFC
  159. REAL (kind=kind_phys) :: &
  160. DELT, &
  161. DPSHC, &
  162. RDELT, &
  163. RSEED
  164. REAL (kind=kind_phys), DIMENSION(its:ite) :: &
  165. CLDWRK, &
  166. PS, &
  167. RCS, &
  168. RN, &
  169. SLIMSK, &
  170. HPBL,EVAP,HEAT !Kwon for shallow convection
  171. REAL (kind=kind_phys), DIMENSION(its:ite, kts:kte+1) :: &
  172. PRSI
  173. REAL (kind=kind_phys), DIMENSION(its:ite, kts:kte) :: &
  174. DEL, &
  175. DOT, &
  176. PHIL, &
  177. PRSL, &
  178. PRSLK, &
  179. Q1, &
  180. T1, &
  181. U1, &
  182. V1, &
  183. ZI, &
  184. ZL
  185. REAL (kind=kind_phys), DIMENSION(its:ite, kts:kte, 2) :: &
  186. QL
  187. INTEGER, DIMENSION(its:ite) :: &
  188. KBOT, &
  189. KTOP, &
  190. KCNV
  191. INTEGER :: &
  192. I, &
  193. IGPVS, &
  194. IM, &
  195. J, &
  196. JCAP, &
  197. K, &
  198. KM, &
  199. KP, &
  200. KX, &
  201. NCLOUD
  202. LOGICAL :: run_param , doing_adapt_dt , decided
  203. DATA IGPVS/0/
  204. !-----------------------------------------------------------------------
  205. !
  206. !*** CHECK TO SEE IF THIS IS A CONVECTION TIMESTEP
  207. !
  208. ! Initialization for adaptive time step.
  209. doing_adapt_dt = .FALSE.
  210. IF ( PRESENT(adapt_step_flag) ) THEN
  211. IF ( adapt_step_flag ) THEN
  212. doing_adapt_dt = .TRUE.
  213. IF ( cudtacttime .EQ. 0. ) THEN
  214. cudtacttime = curr_secs + cudt*60.
  215. END IF
  216. END IF
  217. END IF
  218. ! Do we run through this scheme or not?
  219. ! Test 1: If this is the initial model time, then yes.
  220. ! ITIMESTEP=1
  221. ! Test 2: If the user asked for the cumulus to be run every time step, then yes.
  222. ! CUDT=0 or STEPCU=1
  223. ! Test 3: If not adaptive dt, and this is on the requested cumulus frequency, then yes.
  224. ! MOD(ITIMESTEP,STEPCU)=0
  225. ! Test 4: If using adaptive dt and the current time is past the last requested activate cumulus time, then yes.
  226. ! CURR_SECS >= CUDTACTTIME
  227. ! If we do run through the scheme, we set the flag run_param to TRUE and we set the decided flag
  228. ! to TRUE. The decided flag says that one of these tests was able to say "yes", run the scheme.
  229. ! We only proceed to other tests if the previous tests all have left decided as FALSE.
  230. ! If we set run_param to TRUE and this is adaptive time stepping, we set the time to the next
  231. ! cumulus run.
  232. decided = .FALSE.
  233. run_param = .FALSE.
  234. IF ( ( .NOT. decided ) .AND. &
  235. ( itimestep .EQ. 1 ) ) THEN
  236. run_param = .TRUE.
  237. decided = .TRUE.
  238. END IF
  239. IF ( ( .NOT. decided ) .AND. &
  240. ( ( cudt .EQ. 0. ) .OR. ( stepcu .EQ. 1 ) ) ) THEN
  241. run_param = .TRUE.
  242. decided = .TRUE.
  243. END IF
  244. IF ( ( .NOT. decided ) .AND. &
  245. ( .NOT. doing_adapt_dt ) .AND. &
  246. ( MOD(itimestep,stepcu) .EQ. 0 ) ) THEN
  247. run_param = .TRUE.
  248. decided = .TRUE.
  249. END IF
  250. IF ( ( .NOT. decided ) .AND. &
  251. ( doing_adapt_dt ) .AND. &
  252. ( curr_secs .GE. cudtacttime ) ) THEN
  253. run_param = .TRUE.
  254. decided = .TRUE.
  255. cudtacttime = curr_secs + cudt*60
  256. END IF
  257. !-----------------------------------------------------------------------
  258. if(present(shalconv)) then
  259. shalconv_use=shalconv
  260. else
  261. #if (NMM_CORE==1)
  262. shalconv_use=0
  263. #else
  264. #if (EM_CORE==1)
  265. shalconv_use=1
  266. #else
  267. shalconv_use=0
  268. #endif
  269. #endif
  270. endif
  271. if(present(pgcon)) then
  272. pgcon_use = pgcon
  273. else
  274. ! pgcon_use = 0.7 ! Gregory et al. (1997, QJRMS)
  275. pgcon_use = 0.55 ! Zhang & Wu (2003,JAS), used in GFS (25km res spectral)
  276. ! pgcon_use = 0.2 ! HWRF, for model tuning purposes
  277. ! pgcon_use = 0.3 ! GFDL, or so I am told
  278. ! For those attempting to tune pgcon:
  279. ! The value of 0.55 comes from an observational study of
  280. ! synoptic-scale deep convection and 0.7 came from an
  281. ! incorrect fit to the same data. That value is likely
  282. ! correct for deep convection at gridscales near that of GFS,
  283. ! but is questionable in shallow convection, or for scales
  284. ! much finer than synoptic scales.
  285. ! Then again, the assumptions of SAS break down when the
  286. ! gridscale is near the convection scale anyway. In a large
  287. ! storm such as a hurricane, there is often no environment to
  288. ! detrain into since adjancent gridsquares are also undergoing
  289. ! active convection. Each gridsquare will no longer have many
  290. ! updrafts and downdrafts. At sub-convective timescales, you
  291. ! will find unstable columns for many (say, 5 second length)
  292. ! timesteps in a real atmosphere during a convection cell's
  293. ! lifetime, so forcing it to be neutrally stable is unphysical.
  294. ! Hence, in scales near the convection scale (cells have
  295. ! ~0.5-4km diameter in hurricanes), this parameter is more of a
  296. ! tuning parameter to get a scheme that is inappropriate for
  297. ! that resolution to do a reasonable job.
  298. ! Your mileage might vary.
  299. ! - Sam Trahan
  300. endif
  301. if(present(sas_mass_flux)) then
  302. massf=sas_mass_flux
  303. ! Use this to reduce the fluxes added by SAS to prevent
  304. ! computational instability as a result of large fluxes.
  305. else
  306. massf=9e9 ! large number to disable check
  307. endif
  308. if(present(shal_pgcon)) then
  309. if(shal_pgcon>=0) then
  310. shal_pgcon_use = shal_pgcon
  311. else
  312. ! shal_pgcon<0 means use deep pgcon
  313. shal_pgcon_use = pgcon_use
  314. endif
  315. else
  316. ! Default: Same as deep convection pgcon
  317. shal_pgcon_use = pgcon_use
  318. ! Read the warning above though. It may be advisable for
  319. ! these to be different.
  320. endif
  321. if_run_param: IF(run_param) THEN
  322. DO J=JTS,JTE
  323. DO I=ITS,ITE
  324. CU_ACT_FLAG(I,J)=.TRUE.
  325. ENDDO
  326. ENDDO
  327. IM=ITE-ITS+1
  328. KX=KTE-KTS+1
  329. JCAP=126
  330. DPSHC=30_kind_phys
  331. DELT=DT*STEPCU
  332. RDELT=1./DELT
  333. NCLOUD=1
  334. DO J=jms,jme
  335. DO I=ims,ime
  336. PSFC(i,j)=P8w(i,kms,j)
  337. ENDDO
  338. ENDDO
  339. if(igpvs.eq.0) CALL GFUNCPHYS
  340. igpvs=1
  341. !------------- J LOOP (OUTER) --------------------------------------------------
  342. big_outer_j_loop: DO J=jts,jte
  343. ! --------------- compute zi and zl -----------------------------------------
  344. DO i=its,ite
  345. ZI(I,KTS)=0.0
  346. ENDDO
  347. DO k=kts+1,kte
  348. KM=K-1
  349. DO i=its,ite
  350. ZI(I,K)=ZI(I,KM)+dz8w(i,km,j)
  351. ENDDO
  352. ENDDO
  353. DO k=kts+1,kte
  354. KM=K-1
  355. DO i=its,ite
  356. ZL(I,KM)=(ZI(I,K)+ZI(I,KM))*0.5
  357. ENDDO
  358. ENDDO
  359. DO i=its,ite
  360. ZL(I,KTE)=2.*ZI(I,KTE)-ZL(I,KTE-1)
  361. ENDDO
  362. ! --------------- end compute zi and zl -------------------------------------
  363. DO i=its,ite
  364. PS(i)=PSFC(i,j)*.001
  365. RCS(i)=1.
  366. SLIMSK(i)=ABS(XLAND(i,j)-2.)
  367. ENDDO
  368. #if (NMM_CORE == 1)
  369. if(shalconv_use==1) then
  370. DO i=its,ite
  371. HPBL(I) = HPBL2D(I,J) !kwon for shallow convection
  372. EVAP(I) = EVAP2D(I,J) !kwon for shallow convection
  373. HEAT(I) = HEAT2D(I,J) !kwon for shallow convection
  374. ENDDO
  375. endif
  376. #endif
  377. DO i=its,ite
  378. PRSI(i,kts)=PS(i)
  379. ENDDO
  380. DO k=kts,kte
  381. kp=k+1
  382. DO i=its,ite
  383. PRSL(I,K)=Pcps(i,k,j)*.001
  384. PHIL(I,K)=ZL(I,K)*GRAV
  385. DOT(i,k)=-5.0E-4*GRAV*rho3d(i,k,j)*(w(i,k,j)+w(i,kp,j))
  386. ENDDO
  387. ENDDO
  388. DO k=kts,kte
  389. DO i=its,ite
  390. DEL(i,k)=PRSL(i,k)*GRAV/RD*dz8w(i,k,j)/T3D(i,k,j)
  391. U1(i,k)=U3D(i,k,j)
  392. V1(i,k)=V3D(i,k,j)
  393. Q1(i,k)=QV3D(i,k,j)/(1.+QV3D(i,k,j))
  394. T1(i,k)=T3D(i,k,j)
  395. QL(i,k,1)=QI3D(i,k,j)/(1.+QI3D(i,k,j))
  396. QL(i,k,2)=QC3D(i,k,j)/(1.+QC3D(i,k,j))
  397. PRSLK(I,K)=(PRSL(i,k)*.01)**ROVCP
  398. ENDDO
  399. ENDDO
  400. DO k=kts+1,kte+1
  401. km=k-1
  402. DO i=its,ite
  403. PRSI(i,k)=PRSI(i,km)-del(i,km)
  404. ENDDO
  405. ENDDO
  406. CALL SASCNVN(IM,IM,KX,JCAP,DELT,DEL,PRSL,PS,PHIL, &
  407. QL,Q1,T1,U1,V1,RCS,CLDWRK,RN,KBOT, &
  408. KTOP,KCNV,SLIMSK,DOT,NCLOUD,PGCON_USE,massf)
  409. if_shallow_conv: if(shalconv_use==1) then
  410. #if (NMM_CORE == 1)
  411. ! NMM calls the new shallow convection developed by J Han
  412. ! (Added to WRF by Y.Kwon)
  413. call shalcnv(im,im,kx,jcap,delt,del,prsl,ps,phil,ql, &
  414. & q1,t1,u1,v1,rcs,rn,kbot,ktop,kcnv,slimsk, &
  415. & dot,ncloud,hpbl,heat,evap,shal_pgcon_use)
  416. #else
  417. #if (EM_CORE == 1)
  418. ! NOTE: ARW should be able to call the new shalcnv here, but
  419. ! they need to add the three new variables, so I'm leaving the
  420. ! old shallow convection call here - Sam Trahan
  421. CALL OLD_ARW_SHALCV(IM,IM,KX,DELT,DEL,PRSI,PRSL,PRSLK,KCNV,Q1,T1,DPSHC)
  422. #else
  423. ! Shallow convection is untested for other cores.
  424. #endif
  425. #endif
  426. endif if_shallow_conv
  427. DO I=ITS,ITE
  428. RAINCV(I,J)=RN(I)*1000./STEPCU
  429. PRATEC(I,J)=RN(I)*1000./(STEPCU * DT)
  430. HBOT(I,J)=KBOT(I)
  431. HTOP(I,J)=KTOP(I)
  432. ENDDO
  433. DO K=KTS,KTE
  434. DO I=ITS,ITE
  435. RTHCUTEN(I,K,J)=(T1(I,K)-T3D(I,K,J))/PI3D(I,K,J)*RDELT
  436. RQVCUTEN(I,K,J)=(Q1(I,K)/(1.-q1(i,k))-QV3D(I,K,J))*RDELT
  437. ENDDO
  438. ENDDO
  439. !===============================================================================
  440. ! ADD MOMENTUM MIXING TERM AS TENDENCIES. This is gopal's doing for SAS
  441. ! MOMMIX is the reduction factor set to 0.7 by default. Because NMM has
  442. ! divergence damping term, a reducion factor for cumulum mixing may be
  443. ! required otherwise storms were too weak.
  444. !===============================================================================
  445. !
  446. #if (NMM_CORE == 1)
  447. DO K=KTS,KTE
  448. DO I=ITS,ITE
  449. ! RUCUTEN(I,J,K)=MOMMIX*(U1(I,K)-U3D(I,K,J))*RDELT
  450. ! RVCUTEN(I,J,K)=MOMMIX*(V1(I,K)-V3D(I,K,J))*RDELT
  451. RUCUTEN(I,J,K)=(U1(I,K)-U3D(I,K,J))*RDELT
  452. RVCUTEN(I,J,K)=(V1(I,K)-V3D(I,K,J))*RDELT
  453. ENDDO
  454. ENDDO
  455. #endif
  456. IF(P_QC .ge. P_FIRST_SCALAR)THEN
  457. DO K=KTS,KTE
  458. DO I=ITS,ITE
  459. RQCCUTEN(I,K,J)=(QL(I,K,2)/(1.-ql(i,k,2))-QC3D(I,K,J))*RDELT
  460. ENDDO
  461. ENDDO
  462. ENDIF
  463. IF(P_QI .ge. P_FIRST_SCALAR)THEN
  464. DO K=KTS,KTE
  465. DO I=ITS,ITE
  466. RQICUTEN(I,K,J)=(QL(I,K,1)/(1.-ql(i,k,1))-QI3D(I,K,J))*RDELT
  467. ENDDO
  468. ENDDO
  469. ENDIF
  470. ENDDO big_outer_j_loop ! Outer most J loop
  471. ENDIF if_run_param
  472. END SUBROUTINE CU_SAS
  473. !====================================================================
  474. SUBROUTINE sasinit(RTHCUTEN,RQVCUTEN,RQCCUTEN,RQICUTEN, &
  475. RUCUTEN,RVCUTEN, &
  476. RESTART,P_QC,P_QI,P_FIRST_SCALAR, &
  477. allowed_to_read, &
  478. ids, ide, jds, jde, kds, kde, &
  479. ims, ime, jms, jme, kms, kme, &
  480. its, ite, jts, jte, kts, kte )
  481. !--------------------------------------------------------------------
  482. IMPLICIT NONE
  483. !--------------------------------------------------------------------
  484. LOGICAL , INTENT(IN) :: allowed_to_read,restart
  485. INTEGER , INTENT(IN) :: ids, ide, jds, jde, kds, kde, &
  486. ims, ime, jms, jme, kms, kme, &
  487. its, ite, jts, jte, kts, kte
  488. INTEGER , INTENT(IN) :: P_FIRST_SCALAR, P_QI, P_QC
  489. REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(OUT) :: &
  490. RTHCUTEN, &
  491. RQVCUTEN, &
  492. RQCCUTEN, &
  493. RQICUTEN
  494. REAL, DIMENSION( ims:ime , jms:jme , kms:kme ) , INTENT(OUT) :: &
  495. RUCUTEN, & ! gopal's doing for SAS
  496. RVCUTEN
  497. INTEGER :: i, j, k, itf, jtf, ktf
  498. jtf=min0(jte,jde-1)
  499. ktf=min0(kte,kde-1)
  500. itf=min0(ite,ide-1)
  501. #ifdef HWRF
  502. !zhang's doing
  503. IF(.not.restart .or. .not.allowed_to_read)THEN
  504. !end of zhang's doing
  505. #else
  506. IF(.not.restart)THEN
  507. #endif
  508. DO j=jts,jtf
  509. DO k=kts,ktf
  510. DO i=its,itf
  511. RTHCUTEN(i,k,j)=0.
  512. RQVCUTEN(i,k,j)=0.
  513. RUCUTEN(i,j,k)=0.
  514. RVCUTEN(i,j,k)=0.
  515. ENDDO
  516. ENDDO
  517. ENDDO
  518. IF (P_QC .ge. P_FIRST_SCALAR) THEN
  519. DO j=jts,jtf
  520. DO k=kts,ktf
  521. DO i=its,itf
  522. RQCCUTEN(i,k,j)=0.
  523. ENDDO
  524. ENDDO
  525. ENDDO
  526. ENDIF
  527. IF (P_QI .ge. P_FIRST_SCALAR) THEN
  528. DO j=jts,jtf
  529. DO k=kts,ktf
  530. DO i=its,itf
  531. RQICUTEN(i,k,j)=0.
  532. ENDDO
  533. ENDDO
  534. ENDDO
  535. ENDIF
  536. ENDIF
  537. END SUBROUTINE sasinit
  538. ! ------------------------------------------------------------------------
  539. SUBROUTINE SASCNV(IM,IX,KM,JCAP,DELT,DEL,PRSL,PS,PHIL,QL, &
  540. ! SUBROUTINE SASCNV(IM,IX,KM,JCAP,DLT,DEL,PRSL,PHIL,QL, &
  541. & Q1,T1,U1,V1,RCS,CLDWRK,RN,KBOT,KTOP,KUO,SLIMSK, &
  542. & DOT,XKT2,ncloud)
  543. ! for cloud water version
  544. ! parameter(ncloud=0)
  545. ! SUBROUTINE SASCNV(KM,JCAP,DELT,DEL,SL,SLK,PS,QL,
  546. ! & Q1,T1,U1,V1,RCS,CLDWRK,RN,KBOT,KTOP,KUO,SLIMSK,
  547. ! & DOT,xkt2,ncloud)
  548. !
  549. USE MODULE_GFS_MACHINE , ONLY : kind_phys
  550. USE MODULE_GFS_FUNCPHYS ,ONLY : fpvs
  551. USE MODULE_GFS_PHYSCONS, grav => con_g, CP => con_CP, HVAP => con_HVAP &
  552. &, RV => con_RV, FV => con_fvirt, T0C => con_T0C &
  553. &, CVAP => con_CVAP, CLIQ => con_CLIQ &
  554. &, EPS => con_eps, EPSM1 => con_epsm1
  555. implicit none
  556. !
  557. ! include 'constant.h'
  558. !
  559. integer IM, IX, KM, JCAP, ncloud, &
  560. & KBOT(IM), KTOP(IM), KUO(IM), J
  561. real(kind=kind_phys) DELT
  562. real(kind=kind_phys) PS(IM), DEL(IX,KM), PRSL(IX,KM), &
  563. ! real(kind=kind_phys) DEL(IX,KM), PRSL(IX,KM),
  564. & QL(IX,KM,2), Q1(IX,KM), T1(IX,KM), &
  565. & U1(IX,KM), V1(IX,KM), RCS(IM), &
  566. & CLDWRK(IM), RN(IM), SLIMSK(IM), &
  567. & DOT(IX,KM), XKT2(IM), PHIL(IX,KM)
  568. !
  569. integer I, INDX, jmn, k, knumb, latd, lond, km1
  570. !
  571. real(kind=kind_phys) adw, alpha, alphal, alphas, &
  572. & aup, beta, betal, betas, &
  573. & c0, cpoel, dellat, delta, &
  574. & desdt, deta, detad, dg, &
  575. & dh, dhh, dlnsig, dp, &
  576. & dq, dqsdp, dqsdt, dt, &
  577. & dt2, dtmax, dtmin, dv1, &
  578. & dv1q, dv2, dv2q, dv1u, &
  579. & dv1v, dv2u, dv2v, dv3u, &
  580. & dv3v, dv3, dv3q, dvq1, &
  581. & dz, dz1, e1, edtmax, &
  582. & edtmaxl, edtmaxs, el2orc, elocp, &
  583. & es, etah, &
  584. & evef, evfact, evfactl, fact1, &
  585. & fact2, factor, fjcap, fkm, &
  586. & fuv, g, gamma, onemf, &
  587. & onemfu, pdetrn, pdpdwn, pprime, &
  588. & qc, qlk, qrch, qs, &
  589. & rain, rfact, shear, tem1, &
  590. & tem2, terr, val, val1, &
  591. & val2, w1, w1l, w1s, &
  592. & w2, w2l, w2s, w3, &
  593. & w3l, w3s, w4, w4l, &
  594. & w4s, xdby, xpw, xpwd, &
  595. & xqc, xqrch, xlambu, mbdt, &
  596. & tem
  597. !
  598. !
  599. integer JMIN(IM), KB(IM), KBCON(IM), KBDTR(IM), &
  600. & KT2(IM), KTCON(IM), LMIN(IM), &
  601. & kbm(IM), kbmax(IM), kmax(IM)
  602. !
  603. real(kind=kind_phys) AA1(IM), ACRT(IM), ACRTFCT(IM), &
  604. & DELHBAR(IM), DELQ(IM), DELQ2(IM), &
  605. & DELQBAR(IM), DELQEV(IM), DELTBAR(IM), &
  606. & DELTV(IM), DTCONV(IM), EDT(IM), &
  607. & EDTO(IM), EDTX(IM), FLD(IM), &
  608. & HCDO(IM), HKBO(IM), HMAX(IM), &
  609. & HMIN(IM), HSBAR(IM), UCDO(IM), &
  610. & UKBO(IM), VCDO(IM), VKBO(IM), &
  611. & PBCDIF(IM), PDOT(IM), PO(IM,KM), &
  612. & PWAVO(IM), PWEVO(IM), &
  613. ! & PSFC(IM), PWAVO(IM), PWEVO(IM), &
  614. & QCDO(IM), QCOND(IM), QEVAP(IM), &
  615. & QKBO(IM), RNTOT(IM), VSHEAR(IM), &
  616. & XAA0(IM), XHCD(IM), XHKB(IM), &
  617. & XK(IM), XLAMB(IM), XLAMD(IM), &
  618. & XMB(IM), XMBMAX(IM), XPWAV(IM), &
  619. & XPWEV(IM), XQCD(IM), XQKB(IM)
  620. !
  621. ! PHYSICAL PARAMETERS
  622. PARAMETER(G=grav)
  623. PARAMETER(CPOEL=CP/HVAP,ELOCP=HVAP/CP, &
  624. & EL2ORC=HVAP*HVAP/(RV*CP))
  625. PARAMETER(TERR=0.,C0=.002,DELTA=fv)
  626. PARAMETER(FACT1=(CVAP-CLIQ)/RV,FACT2=HVAP/RV-FACT1*T0C)
  627. ! LOCAL VARIABLES AND ARRAYS
  628. real(kind=kind_phys) PFLD(IM,KM), TO(IM,KM), QO(IM,KM), &
  629. & UO(IM,KM), VO(IM,KM), QESO(IM,KM)
  630. ! cloud water
  631. real(kind=kind_phys) QLKO_KTCON(IM), DELLAL(IM), TVO(IM,KM), &
  632. & DBYO(IM,KM), ZO(IM,KM), SUMZ(IM,KM), &
  633. & SUMH(IM,KM), HEO(IM,KM), HESO(IM,KM), &
  634. & QRCD(IM,KM), DELLAH(IM,KM), DELLAQ(IM,KM),&
  635. & DELLAU(IM,KM), DELLAV(IM,KM), HCKO(IM,KM), &
  636. & UCKO(IM,KM), VCKO(IM,KM), QCKO(IM,KM), &
  637. & ETA(IM,KM), ETAU(IM,KM), ETAD(IM,KM), &
  638. & QRCDO(IM,KM), PWO(IM,KM), PWDO(IM,KM), &
  639. & RHBAR(IM), TX1(IM)
  640. !
  641. LOGICAL TOTFLG, CNVFLG(IM), DWNFLG(IM), DWNFLG2(IM), FLG(IM)
  642. !
  643. real(kind=kind_phys) PCRIT(15), ACRITT(15), ACRIT(15)
  644. ! SAVE PCRIT, ACRITT
  645. DATA PCRIT/850.,800.,750.,700.,650.,600.,550.,500.,450.,400., &
  646. & 350.,300.,250.,200.,150./
  647. DATA ACRITT/.0633,.0445,.0553,.0664,.075,.1082,.1521,.2216, &
  648. & .3151,.3677,.41,.5255,.7663,1.1686,1.6851/
  649. ! GDAS DERIVED ACRIT
  650. ! DATA ACRITT/.203,.515,.521,.566,.625,.665,.659,.688, &
  651. ! & .743,.813,.886,.947,1.138,1.377,1.896/
  652. !
  653. real(kind=kind_phys) TF, TCR, TCRF, RZERO, RONE
  654. parameter (TF=233.16, TCR=263.16, TCRF=1.0/(TCR-TF))
  655. parameter (RZERO=0.0,RONE=1.0)
  656. !-----------------------------------------------------------------------
  657. !
  658. km1 = km - 1
  659. ! INITIALIZE ARRAYS
  660. !
  661. DO I=1,IM
  662. RN(I)=0.
  663. KBOT(I)=KM+1
  664. KTOP(I)=0
  665. KUO(I)=0
  666. CNVFLG(I) = .TRUE.
  667. DTCONV(I) = 3600.
  668. CLDWRK(I) = 0.
  669. PDOT(I) = 0.
  670. KT2(I) = 0
  671. QLKO_KTCON(I) = 0.
  672. DELLAL(I) = 0.
  673. ENDDO
  674. !!
  675. DO K = 1, 15
  676. ACRIT(K) = ACRITT(K) * (975. - PCRIT(K))
  677. ENDDO
  678. DT2 = DELT
  679. !cmr dtmin = max(dt2,1200.)
  680. val = 1200.
  681. dtmin = max(dt2, val )
  682. !cmr dtmax = max(dt2,3600.)
  683. val = 3600.
  684. dtmax = max(dt2, val )
  685. ! MODEL TUNABLE PARAMETERS ARE ALL HERE
  686. MBDT = 10.
  687. EDTMAXl = .3
  688. EDTMAXs = .3
  689. ALPHAl = .5
  690. ALPHAs = .5
  691. BETAl = .15
  692. betas = .15
  693. BETAl = .05
  694. betas = .05
  695. ! change for hurricane model
  696. BETAl = .5
  697. betas = .5
  698. ! EVEF = 0.07
  699. evfact = 0.3
  700. evfactl = 0.3
  701. ! change for hurricane model
  702. evfact = 0.6
  703. evfactl = .6
  704. #if ( EM_CORE == 1 )
  705. ! HAWAII TEST - ZCX
  706. ALPHAl = .5
  707. ALPHAs = .75
  708. BETAl = .05
  709. betas = .05
  710. evfact = 0.5
  711. evfactl = 0.5
  712. #endif
  713. PDPDWN = 0.
  714. PDETRN = 200.
  715. xlambu = 1.e-4
  716. fjcap = (float(jcap) / 126.) ** 2
  717. !cmr fjcap = max(fjcap,1.)
  718. val = 1.
  719. fjcap = max(fjcap,val)
  720. fkm = (float(km) / 28.) ** 2
  721. !cmr fkm = max(fkm,1.)
  722. fkm = max(fkm,val)
  723. W1l = -8.E-3
  724. W2l = -4.E-2
  725. W3l = -5.E-3
  726. W4l = -5.E-4
  727. W1s = -2.E-4
  728. W2s = -2.E-3
  729. W3s = -1.E-3
  730. W4s = -2.E-5
  731. !CCCC IF(IM.EQ.384) THEN
  732. LATD = 92
  733. lond = 189
  734. !CCCC ELSEIF(IM.EQ.768) THEN
  735. !CCCC LATD = 80
  736. !CCCC ELSE
  737. !CCCC LATD = 0
  738. !CCCC ENDIF
  739. !
  740. ! DEFINE TOP LAYER FOR SEARCH OF THE DOWNDRAFT ORIGINATING LAYER
  741. ! AND THE MAXIMUM THETAE FOR UPDRAFT
  742. !
  743. DO I=1,IM
  744. KBMAX(I) = KM
  745. KBM(I) = KM
  746. KMAX(I) = KM
  747. TX1(I) = 1.0 / PS(I)
  748. ENDDO
  749. !
  750. DO K = 1, KM
  751. DO I=1,IM
  752. IF (prSL(I,K)*tx1(I) .GT. 0.45) KBMAX(I) = K + 1
  753. IF (prSL(I,K)*tx1(I) .GT. 0.70) KBM(I) = K + 1
  754. IF (prSL(I,K)*tx1(I) .GT. 0.04) KMAX(I) = MIN(KM,K + 1)
  755. ENDDO
  756. ENDDO
  757. DO I=1,IM
  758. KBMAX(I) = MIN(KBMAX(I),KMAX(I))
  759. KBM(I) = MIN(KBM(I),KMAX(I))
  760. ENDDO
  761. !
  762. ! CONVERT SURFACE PRESSURE TO MB FROM CB
  763. !
  764. !!
  765. DO K = 1, KM
  766. DO I=1,IM
  767. if (K .le. kmax(i)) then
  768. PFLD(I,k) = PRSL(I,K) * 10.0
  769. PWO(I,k) = 0.
  770. PWDO(I,k) = 0.
  771. TO(I,k) = T1(I,k)
  772. QO(I,k) = Q1(I,k)
  773. UO(I,k) = U1(I,k)
  774. VO(I,k) = V1(I,k)
  775. DBYO(I,k) = 0.
  776. SUMZ(I,k) = 0.
  777. SUMH(I,k) = 0.
  778. endif
  779. ENDDO
  780. ENDDO
  781. !
  782. ! COLUMN VARIABLES
  783. ! P IS PRESSURE OF THE LAYER (MB)
  784. ! T IS TEMPERATURE AT T-DT (K)..TN
  785. ! Q IS MIXING RATIO AT T-DT (KG/KG)..QN
  786. ! TO IS TEMPERATURE AT T+DT (K)... THIS IS AFTER ADVECTION AND TURBULAN
  787. ! QO IS MIXING RATIO AT T+DT (KG/KG)..Q1
  788. !
  789. DO K = 1, KM
  790. DO I=1,IM
  791. if (k .le. kmax(i)) then
  792. !jfe QESO(I,k) = 10. * FPVS(T1(I,k))
  793. !
  794. QESO(I,k) = 0.01 * fpvs(T1(I,K)) ! fpvs is in Pa
  795. !
  796. QESO(I,k) = EPS * QESO(I,k) / (PFLD(I,k) + EPSM1*QESO(I,k))
  797. !cmr QESO(I,k) = MAX(QESO(I,k),1.E-8)
  798. val1 = 1.E-8
  799. QESO(I,k) = MAX(QESO(I,k), val1)
  800. !cmr QO(I,k) = max(QO(I,k),1.e-10)
  801. val2 = 1.e-10
  802. QO(I,k) = max(QO(I,k), val2 )
  803. ! QO(I,k) = MIN(QO(I,k),QESO(I,k))
  804. TVO(I,k) = TO(I,k) + DELTA * TO(I,k) * QO(I,k)
  805. endif
  806. ENDDO
  807. ENDDO
  808. !
  809. ! HYDROSTATIC HEIGHT ASSUME ZERO TERR
  810. !
  811. DO K = 1, KM
  812. DO I=1,IM
  813. ZO(I,k) = PHIL(I,k) / G
  814. ENDDO
  815. ENDDO
  816. ! COMPUTE MOIST STATIC ENERGY
  817. DO K = 1, KM
  818. DO I=1,IM
  819. if (K .le. kmax(i)) then
  820. ! tem = G * ZO(I,k) + CP * TO(I,k)
  821. tem = PHIL(I,k) + CP * TO(I,k)
  822. HEO(I,k) = tem + HVAP * QO(I,k)
  823. HESO(I,k) = tem + HVAP * QESO(I,k)
  824. ! HEO(I,k) = MIN(HEO(I,k),HESO(I,k))
  825. endif
  826. ENDDO
  827. ENDDO
  828. !
  829. ! DETERMINE LEVEL WITH LARGEST MOIST STATIC ENERGY
  830. ! THIS IS THE LEVEL WHERE UPDRAFT STARTS
  831. !
  832. DO I=1,IM
  833. HMAX(I) = HEO(I,1)
  834. KB(I) = 1
  835. ENDDO
  836. !!
  837. DO K = 2, KM
  838. DO I=1,IM
  839. if (k .le. kbm(i)) then
  840. IF(HEO(I,k).GT.HMAX(I).AND.CNVFLG(I)) THEN
  841. KB(I) = K
  842. HMAX(I) = HEO(I,k)
  843. ENDIF
  844. endif
  845. ENDDO
  846. ENDDO
  847. ! DO K = 1, KMAX - 1
  848. ! TOL(k) = .5 * (TO(I,k) + TO(I,k+1))
  849. ! QOL(k) = .5 * (QO(I,k) + QO(I,k+1))
  850. ! QESOL(I,k) = .5 * (QESO(I,k) + QESO(I,k+1))
  851. ! HEOL(I,k) = .5 * (HEO(I,k) + HEO(I,k+1))
  852. ! HESOL(I,k) = .5 * (HESO(I,k) + HESO(I,k+1))
  853. ! ENDDO
  854. DO K = 1, KM1
  855. DO I=1,IM
  856. if (k .le. kmax(i)-1) then
  857. DZ = .5 * (ZO(I,k+1) - ZO(I,k))
  858. DP = .5 * (PFLD(I,k+1) - PFLD(I,k))
  859. !jfe ES = 10. * FPVS(TO(I,k+1))
  860. !
  861. ES = 0.01 * fpvs(TO(I,K+1)) ! fpvs is in Pa
  862. !
  863. PPRIME = PFLD(I,k+1) + EPSM1 * ES
  864. QS = EPS * ES / PPRIME
  865. DQSDP = - QS / PPRIME
  866. DESDT = ES * (FACT1 / TO(I,k+1) + FACT2 / (TO(I,k+1)**2))
  867. DQSDT = QS * PFLD(I,k+1) * DESDT / (ES * PPRIME)
  868. GAMMA = EL2ORC * QESO(I,k+1) / (TO(I,k+1)**2)
  869. DT = (G * DZ + HVAP * DQSDP * DP) / (CP * (1. + GAMMA))
  870. DQ = DQSDT * DT + DQSDP * DP
  871. TO(I,k) = TO(I,k+1) + DT
  872. QO(I,k) = QO(I,k+1) + DQ
  873. PO(I,k) = .5 * (PFLD(I,k) + PFLD(I,k+1))
  874. endif
  875. ENDDO
  876. ENDDO
  877. !
  878. DO K = 1, KM1
  879. DO I=1,IM
  880. if (k .le. kmax(I)-1) then
  881. !jfe QESO(I,k) = 10. * FPVS(TO(I,k))
  882. !
  883. QESO(I,k) = 0.01 * fpvs(TO(I,K)) ! fpvs is in Pa
  884. !
  885. QESO(I,k) = EPS * QESO(I,k) / (PO(I,k) + EPSM1*QESO(I,k))
  886. !cmr QESO(I,k) = MAX(QESO(I,k),1.E-8)
  887. val1 = 1.E-8
  888. QESO(I,k) = MAX(QESO(I,k), val1)
  889. !cmr QO(I,k) = max(QO(I,k),1.e-10)
  890. val2 = 1.e-10
  891. QO(I,k) = max(QO(I,k), val2 )
  892. ! QO(I,k) = MIN(QO(I,k),QESO(I,k))
  893. HEO(I,k) = .5 * G * (ZO(I,k) + ZO(I,k+1)) + &
  894. & CP * TO(I,k) + HVAP * QO(I,k)
  895. HESO(I,k) = .5 * G * (ZO(I,k) + ZO(I,k+1)) + &
  896. & CP * TO(I,k) + HVAP * QESO(I,k)
  897. UO(I,k) = .5 * (UO(I,k) + UO(I,k+1))
  898. VO(I,k) = .5 * (VO(I,k) + VO(I,k+1))
  899. endif
  900. ENDDO
  901. ENDDO
  902. ! k = kmax
  903. ! HEO(I,k) = HEO(I,k)
  904. ! hesol(k) = HESO(I,k)
  905. ! IF(LAT.EQ.LATD.AND.lon.eq.lond.and.CNVFLG(I)) THEN
  906. ! PRINT *, ' HEO ='
  907. ! PRINT 6001, (HEO(I,K),K=1,KMAX)
  908. ! PRINT *, ' HESO ='
  909. ! PRINT 6001, (HESO(I,K),K=1,KMAX)
  910. ! PRINT *, ' TO ='
  911. ! PRINT 6002, (TO(I,K)-273.16,K=1,KMAX)
  912. ! PRINT *, ' QO ='
  913. ! PRINT 6003, (QO(I,K),K=1,KMAX)
  914. ! PRINT *, ' QSO ='
  915. ! PRINT 6003, (QESO(I,K),K=1,KMAX)
  916. ! ENDIF
  917. !
  918. ! LOOK FOR CONVECTIVE CLOUD BASE AS THE LEVEL OF FREE CONVECTION
  919. !
  920. DO I=1,IM
  921. IF(CNVFLG(I)) THEN
  922. INDX = KB(I)
  923. HKBO(I) = HEO(I,INDX)
  924. QKBO(I) = QO(I,INDX)
  925. UKBO(I) = UO(I,INDX)
  926. VKBO(I) = VO(I,INDX)
  927. ENDIF
  928. FLG(I) = CNVFLG(I)
  929. KBCON(I) = KMAX(I)
  930. ENDDO
  931. !!
  932. DO K = 1, KM
  933. DO I=1,IM
  934. if (k .le. kbmax(i)) then
  935. IF(FLG(I).AND.K.GT.KB(I)) THEN
  936. HSBAR(I) = HESO(I,k)
  937. IF(HKBO(I).GT.HSBAR(I)) THEN
  938. FLG(I) = .FALSE.
  939. KBCON(I) = K
  940. ENDIF
  941. ENDIF
  942. endif
  943. ENDDO
  944. ENDDO
  945. DO I=1,IM
  946. IF(CNVFLG(I)) THEN
  947. PBCDIF(I) = -PFLD(I,KBCON(I)) + PFLD(I,KB(I))
  948. PDOT(I) = 10.* DOT(I,KBCON(I))
  949. IF(PBCDIF(I).GT.150.) CNVFLG(I) = .FALSE.
  950. IF(KBCON(I).EQ.KMAX(I)) CNVFLG(I) = .FALSE.
  951. ENDIF
  952. ENDDO
  953. !!
  954. TOTFLG = .TRUE.
  955. DO I=1,IM
  956. TOTFLG = TOTFLG .AND. (.NOT. CNVFLG(I))
  957. ENDDO
  958. IF(TOTFLG) RETURN
  959. ! FOUND LFC, CAN DEFINE REST OF VARIABLES
  960. 6001 FORMAT(2X,-2P10F12.2)
  961. 6002 FORMAT(2X,10F12.2)
  962. 6003 FORMAT(2X,3P10F12.2)
  963. !
  964. ! DETERMINE ENTRAINMENT RATE BETWEEN KB AND KBCON
  965. !
  966. DO I = 1, IM
  967. alpha = alphas
  968. if(SLIMSK(I).eq.1.) alpha = alphal
  969. IF(CNVFLG(I)) THEN
  970. IF(KB(I).EQ.1) THEN
  971. DZ = .5 * (ZO(I,KBCON(I)) + ZO(I,KBCON(I)-1)) - ZO(I,1)
  972. ELSE
  973. DZ = .5 * (ZO(I,KBCON(I)) + ZO(I,KBCON(I)-1)) &
  974. & - .5 * (ZO(I,KB(I)) + ZO(I,KB(I)-1))
  975. ENDIF
  976. IF(KBCON(I).NE.KB(I)) THEN
  977. !cmr XLAMB(I) = -ALOG(ALPHA) / DZ
  978. XLAMB(I) = - LOG(ALPHA) / DZ
  979. ELSE
  980. XLAMB(I) = 0.
  981. ENDIF
  982. ENDIF
  983. ENDDO
  984. ! DETERMINE UPDRAFT MASS FLUX
  985. DO K = 1, KM
  986. DO I = 1, IM
  987. if (k .le. kmax(i) .and. CNVFLG(I)) then
  988. ETA(I,k) = 1.
  989. ETAU(I,k) = 1.
  990. ENDIF
  991. ENDDO
  992. ENDDO
  993. DO K = KM1, 2, -1
  994. DO I = 1, IM
  995. if (k .le. kbmax(i)) then
  996. IF(CNVFLG(I).AND.K.LT.KBCON(I).AND.K.GE.KB(I)) THEN
  997. DZ = .5 * (ZO(I,k+1) - ZO(I,k-1))
  998. ETA(I,k) = ETA(I,k+1) * EXP(-XLAMB(I) * DZ)
  999. ETAU(I,k) = ETA(I,k)
  1000. ENDIF
  1001. endif
  1002. ENDDO
  1003. ENDDO
  1004. DO I = 1, IM
  1005. IF(CNVFLG(I).AND.KB(I).EQ.1.AND.KBCON(I).GT.1) THEN
  1006. DZ = .5 * (ZO(I,2) - ZO(I,1))
  1007. ETA(I,1) = ETA(I,2) * EXP(-XLAMB(I) * DZ)
  1008. ETAU(I,1) = ETA(I,1)
  1009. ENDIF
  1010. ENDDO
  1011. !
  1012. ! WORK UP UPDRAFT CLOUD PROPERTIES
  1013. !
  1014. DO I = 1, IM
  1015. IF(CNVFLG(I)) THEN
  1016. INDX = KB(I)
  1017. HCKO(I,INDX) = HKBO(I)
  1018. QCKO(I,INDX) = QKBO(I)
  1019. UCKO(I,INDX) = UKBO(I)
  1020. VCKO(I,INDX) = VKBO(I)
  1021. PWAVO(I) = 0.
  1022. ENDIF
  1023. ENDDO
  1024. !
  1025. ! CLOUD PROPERTY BELOW CLOUD BASE IS MODIFIED BY THE ENTRAINMENT PROCES
  1026. !
  1027. DO K = 2, KM1
  1028. DO I = 1, IM
  1029. if (k .le. kmax(i)-1) then
  1030. IF(CNVFLG(I).AND.K.GT.KB(I).AND.K.LE.KBCON(I)) THEN
  1031. FACTOR = ETA(I,k-1) / ETA(I,k)
  1032. ONEMF = 1. - FACTOR
  1033. HCKO(I,k) = FACTOR * HCKO(I,k-1) + ONEMF * &
  1034. & .5 * (HEO(I,k) + HEO(I,k+1))
  1035. UCKO(I,k) = FACTOR * UCKO(I,k-1) + ONEMF * &
  1036. & .5 * (UO(I,k) + UO(I,k+1))
  1037. VCKO(I,k) = FACTOR * VCKO(I,k-1) + ONEMF * &
  1038. & .5 * (VO(I,k) + VO(I,k+1))
  1039. DBYO(I,k) = HCKO(I,k) - HESO(I,k)
  1040. ENDIF
  1041. IF(CNVFLG(I).AND.K.GT.KBCON(I)) THEN
  1042. HCKO(I,k) = HCKO(I,k-1)
  1043. UCKO(I,k) = UCKO(I,k-1)
  1044. VCKO(I,k) = VCKO(I,k-1)
  1045. DBYO(I,k) = HCKO(I,k) - HESO(I,k)
  1046. ENDIF
  1047. endif
  1048. ENDDO
  1049. ENDDO
  1050. ! DETERMINE CLOUD TOP
  1051. DO I = 1, IM
  1052. FLG(I) = CNVFLG(I)
  1053. KTCON(I) = 1
  1054. ENDDO
  1055. ! DO K = 2, KMAX
  1056. ! KK = KMAX - K + 1
  1057. ! IF(DBYO(I,kK).GE.0..AND.FLG(I).AND.KK.GT.KBCON(I)) THEN
  1058. ! KTCON(I) = KK + 1
  1059. ! FLG(I) = .FALSE.
  1060. ! ENDIF
  1061. ! ENDDO
  1062. DO K = 2, KM
  1063. DO I = 1, IM
  1064. if (k .le. kmax(i)) then
  1065. IF(DBYO(I,k).LT.0..AND.FLG(I).AND.K.GT.KBCON(I)) THEN
  1066. KTCON(I) = K
  1067. FLG(I) = .FALSE.
  1068. ENDIF
  1069. endif
  1070. ENDDO
  1071. ENDDO
  1072. DO I = 1, IM
  1073. IF(CNVFLG(I).AND.(PFLD(I,KBCON(I)) - PFLD(I,KTCON(I))).LT.150.) &
  1074. & CNVFLG(I) = .FALSE.
  1075. ENDDO
  1076. TOTFLG = .TRUE.
  1077. DO I = 1, IM
  1078. TOTFLG = TOTFLG .AND. (.NOT. CNVFLG(I))
  1079. ENDDO
  1080. IF(TOTFLG) RETURN
  1081. !
  1082. ! SEARCH FOR DOWNDRAFT ORIGINATING LEVEL ABOVE THETA-E MINIMUM
  1083. !
  1084. DO I = 1, IM
  1085. HMIN(I) = HEO(I,KBCON(I))
  1086. LMIN(I) = KBMAX(I)
  1087. JMIN(I) = KBMAX(I)
  1088. ENDDO
  1089. DO I = 1, IM
  1090. DO K = KBCON(I), KBMAX(I)
  1091. IF(HEO(I,k).LT.HMIN(I).AND.CNVFLG(I)) THEN
  1092. LMIN(I) = K + 1
  1093. HMIN(I) = HEO(I,k)
  1094. ENDIF
  1095. ENDDO
  1096. ENDDO
  1097. !
  1098. ! Make sure that JMIN(I) is within the cloud
  1099. !
  1100. DO I = 1, IM
  1101. IF(CNVFLG(I)) THEN
  1102. JMIN(I) = MIN(LMIN(I),KTCON(I)-1)
  1103. XMBMAX(I) = .1
  1104. JMIN(I) = MAX(JMIN(I),KBCON(I)+1)
  1105. ENDIF
  1106. ENDDO
  1107. !
  1108. ! ENTRAINING CLOUD
  1109. !
  1110. do k = 2, km1
  1111. DO I = 1, IM
  1112. if (k .le. kmax(i)-1) then
  1113. if(CNVFLG(I).and.k.gt.JMIN(I).and.k.le.KTCON(I)) THEN
  1114. SUMZ(I,k) = SUMZ(I,k-1) + .5 * (ZO(I,k+1) - ZO(I,k-1))
  1115. SUMH(I,k) = SUMH(I,k-1) + .5 * (ZO(I,k+1) - ZO(I,k-1)) &
  1116. & * HEO(I,k)
  1117. ENDIF
  1118. endif
  1119. enddo
  1120. enddo
  1121. !!
  1122. DO I = 1, IM
  1123. IF(CNVFLG(I)) THEN
  1124. ! call random_number(XKT2)
  1125. ! call srand(fhour)
  1126. ! XKT2(I) = rand()
  1127. KT2(I) = nint(XKT2(I)*float(KTCON(I)-JMIN(I))-.5)+JMIN(I)+1
  1128. ! KT2(I) = nint(sqrt(XKT2(I))*float(KTCON(I)-JMIN(I))-.5) + JMIN(I) + 1
  1129. ! KT2(I) = nint(ranf() *float(KTCON(I)-JMIN(I))-.5) + JMIN(I) + 1
  1130. tem1 = (HCKO(I,JMIN(I)) - HESO(I,KT2(I)))
  1131. tem2 = (SUMZ(I,KT2(I)) * HESO(I,KT2(I)) - SUMH(I,KT2(I)))
  1132. if (abs(tem2) .gt. 0.000001) THEN
  1133. XLAMB(I) = tem1 / tem2
  1134. else
  1135. CNVFLG(I) = .false.
  1136. ENDIF
  1137. ! XLAMB(I) = (HCKO(I,JMIN(I)) - HESO(I,KT2(I)))
  1138. ! & / (SUMZ(I,KT2(I)) * HESO(I,KT2(I)) - SUMH(I,KT2(I)))
  1139. XLAMB(I) = max(XLAMB(I),RZERO)
  1140. XLAMB(I) = min(XLAMB(I),2.3/SUMZ(I,KT2(I)))
  1141. ENDIF
  1142. ENDDO
  1143. !!
  1144. DO I = 1, IM
  1145. DWNFLG(I) = CNVFLG(I)
  1146. DWNFLG2(I) = CNVFLG(I)
  1147. IF(CNVFLG(I)) THEN
  1148. if(KT2(I).ge.KTCON(I)) DWNFLG(I) = .false.
  1149. if(XLAMB(I).le.1.e-30.or.HCKO(I,JMIN(I))-HESO(I,KT2(I)).le.1.e-30)&
  1150. & DWNFLG(I) = .false.
  1151. do k = JMIN(I), KT2(I)
  1152. if(DWNFLG(I).and.HEO(I,k).gt.HESO(I,KT2(I))) DWNFLG(I)=.false.
  1153. enddo
  1154. ! IF(CNVFLG(I).AND.(PFLD(KBCON(I))-PFLD(KTCON(I))).GT.PDETRN)
  1155. ! & DWNFLG(I)=.FALSE.
  1156. IF(CNVFLG(I).AND.(PFLD(I,KBCON(I))-PFLD(I,KTCON(I))).LT.PDPDWN) &
  1157. & DWNFLG2(I)=.FALSE.
  1158. ENDIF
  1159. ENDDO
  1160. !!
  1161. DO K = 2, KM1
  1162. DO I = 1, IM
  1163. if (k .le. kmax(i)-1) then
  1164. IF(DWNFLG(I).AND.K.GT.JMIN(I).AND.K.LE.KT2(I)) THEN
  1165. DZ = .5 * (ZO(I,k+1) - ZO(I,k-1))
  1166. ! ETA(I,k) = ETA(I,k-1) * EXP( XLAMB(I) * DZ)
  1167. ! to simplify matter, we will take the linear approach here
  1168. !
  1169. ETA(I,k) = ETA(I,k-1) * (1. + XLAMB(I) * dz)
  1170. ETAU(I,k) = ETAU(I,k-1) * (1. + (XLAMB(I)+xlambu) * dz)
  1171. ENDIF
  1172. endif
  1173. ENDDO
  1174. ENDDO
  1175. !!
  1176. DO K = 2, KM1
  1177. DO I = 1, IM
  1178. if (k .le. kmax(i)-1) then
  1179. ! IF(.NOT.DWNFLG(I).AND.K.GT.JMIN(I).AND.K.LE.KT2(I)) THEN
  1180. IF(.NOT.DWNFLG(I).AND.K.GT.JMIN(I).AND.K.LE.KTCON(I)) THEN
  1181. DZ = .5 * (ZO(I,k+1) - ZO(I,k-1))
  1182. ETAU(I,k) = ETAU(I,k-1) * (1. + xlambu * dz)
  1183. ENDIF
  1184. endif
  1185. ENDDO
  1186. ENDDO
  1187. ! IF(LAT.EQ.LATD.AND.lon.eq.lond.and.CNVFLG(I)) THEN
  1188. ! PRINT *, ' LMIN(I), KT2(I)=', LMIN(I), KT2(I)
  1189. ! PRINT *, ' KBOT, KTOP, JMIN(I) =', KBCON(I), KTCON(I), JMIN(I)
  1190. ! ENDIF
  1191. ! IF(LAT.EQ.LATD.AND.lon.eq.lond) THEN
  1192. ! print *, ' xlamb =', xlamb
  1193. ! print *, ' eta =', (eta(k),k=1,KT2(I))
  1194. ! print *, ' ETAU =', (ETAU(I,k),k=1,KT2(I))
  1195. ! print *, ' HCKO =', (HCKO(I,k),k=1,KT2(I))
  1196. ! print *, ' SUMZ =', (SUMZ(I,k),k=1,KT2(I))
  1197. ! print *, ' SUMH =', (SUMH(I,k),k=1,KT2(I))
  1198. ! ENDIF
  1199. DO I = 1, IM
  1200. if(DWNFLG(I)) THEN
  1201. KTCON(I) = KT2(I)
  1202. ENDIF
  1203. ENDDO
  1204. !
  1205. ! CLOUD PROPERTY ABOVE CLOUD Base IS MODIFIED BY THE DETRAINMENT PROCESS
  1206. !
  1207. DO K = 2, KM1
  1208. DO I = 1, IM
  1209. if (k .le. kmax(i)-1) then
  1210. !jfe
  1211. IF(CNVFLG(I).AND.K.GT.KBCON(I).AND.K.LE.KTCON(I)) THEN
  1212. !jfe IF(K.GT.KBCON(I).AND.K.LE.KTCON(I)) THEN
  1213. FACTOR = ETA(I,k-1) / ETA(I,k)
  1214. ONEMF = 1. - FACTOR
  1215. fuv = ETAU(I,k-1) / ETAU(I,k)
  1216. onemfu = 1. - fuv
  1217. HCKO(I,k) = FACTOR * HCKO(I,k-1) + ONEMF * &
  1218. & .5 * (HEO(I,k) + HEO(I,k+1))
  1219. UCKO(I,k) = fuv * UCKO(I,k-1) + ONEMFu * &
  1220. & .5 * (UO(I,k) + UO(I,k+1))
  1221. VCKO(I,k) = fuv * VCKO(I,k-1) + ONEMFu * &
  1222. & .5 * (VO(I,k) + VO(I,k+1))
  1223. DBYO(I,k) = HCKO(I,k) - HESO(I,k)
  1224. ENDIF
  1225. endif
  1226. ENDDO
  1227. ENDDO
  1228. ! IF(LAT.EQ.LATD.AND.lon.eq.lond.and.CNVFLG(I)) THEN
  1229. ! PRINT *, ' UCKO=', (UCKO(I,k),k=KBCON(I)+1,KTCON(I))
  1230. ! PRINT *, ' uenv=', (.5*(UO(I,k)+UO(I,k-1)),k=KBCON(I)+1,KTCON(I))
  1231. ! ENDIF
  1232. DO I = 1, IM
  1233. if(CNVFLG(I).and.DWNFLG2(I).and.JMIN(I).le.KBCON(I)) &
  1234. & THEN
  1235. CNVFLG(I) = .false.
  1236. DWNFLG(I) = .false.
  1237. DWNFLG2(I) = .false.
  1238. ENDIF
  1239. ENDDO
  1240. !!
  1241. TOTFLG = .TRUE.
  1242. DO I = 1, IM
  1243. TOTFLG = TOTFLG .AND. (.NOT. CNVFLG(I))
  1244. ENDDO
  1245. IF(TOTFLG) RETURN
  1246. !!
  1247. !
  1248. ! COMPUTE CLOUD MOISTURE PROPERTY AND PRECIPITATION
  1249. !
  1250. DO I = 1, IM
  1251. AA1(I) = 0.
  1252. RHBAR(I) = 0.
  1253. ENDDO
  1254. DO K = 1, KM
  1255. DO I = 1, IM
  1256. if (k .le. kmax(i)) then
  1257. IF(CNVFLG(I).AND.K.GT.KB(I).AND.K.LT.KTCON(I)) THEN
  1258. DZ = .5 * (ZO(I,k+1) - ZO(I,k-1))
  1259. DZ1 = (ZO(I,k) - ZO(I,k-1))
  1260. GAMMA = EL2ORC * QESO(I,k) / (TO(I,k)**2)
  1261. QRCH = QESO(I,k) &
  1262. & + GAMMA * DBYO(I,k) / (HVAP * (1. + GAMMA))
  1263. FACTOR = ETA(I,k-1) / ETA(I,k)
  1264. ONEMF = 1. - FACTOR
  1265. QCKO(I,k) = FACTOR * QCKO(I,k-1) + ONEMF * &
  1266. & .5 * (QO(I,k) + QO(I,k+1))
  1267. DQ = ETA(I,k) * QCKO(I,k) - ETA(I,k) * QRCH
  1268. RHBAR(I) = RHBAR(I) + QO(I,k) / QESO(I,k)
  1269. !
  1270. ! BELOW LFC CHECK IF THERE IS EXCESS MOISTURE TO RELEASE LATENT HEAT
  1271. !
  1272. IF(DQ.GT.0.) THEN
  1273. ETAH = .5 * (ETA(I,k) + ETA(I,k-1))
  1274. QLK = DQ / (ETA(I,k) + ETAH * C0 * DZ)
  1275. AA1(I) = AA1(I) - DZ1 * G * QLK
  1276. QC = QLK + QRCH
  1277. PWO(I,k) = ETAH * C0 * DZ * QLK
  1278. QCKO(I,k) = QC
  1279. PWAVO(I) = PWAVO(I) + PWO(I,k)
  1280. ENDIF
  1281. ENDIF
  1282. endif
  1283. ENDDO
  1284. ENDDO
  1285. DO I = 1, IM
  1286. RHBAR(I) = RHBAR(I) / float(KTCON(I) - KB(I) - 1)
  1287. ENDDO
  1288. !
  1289. ! this section is ready for cloud water
  1290. !
  1291. if(ncloud.gt.0) THEN
  1292. !
  1293. ! compute liquid and vapor separation at cloud top
  1294. !
  1295. DO I = 1, IM
  1296. k = KTCON(I)
  1297. IF(CNVFLG(I)) THEN
  1298. GAMMA = EL2ORC * QESO(I,K) / (TO(I,K)**2)
  1299. QRCH = QESO(I,K) &
  1300. & + GAMMA * DBYO(I,K) / (HVAP * (1. + GAMMA))
  1301. DQ = QCKO(I,K-1) - QRCH
  1302. !
  1303. ! CHECK IF THERE IS EXCESS MOISTURE TO RELEASE LATENT HEAT
  1304. !
  1305. IF(DQ.GT.0.) THEN
  1306. QLKO_KTCON(I) = dq
  1307. QCKO(I,K-1) = QRCH
  1308. ENDIF
  1309. ENDIF
  1310. ENDDO
  1311. ENDIF
  1312. !
  1313. ! CALCULATE CLOUD WORK FUNCTION AT T+DT
  1314. !
  1315. DO K = 1, KM
  1316. DO I = 1, IM
  1317. if (k .le. kmax(i)) then
  1318. IF(CNVFLG(I).AND.K.GT.KBCON(I).AND.K.LE.KTCON(I)) THEN
  1319. DZ1 = ZO(I,k) - ZO(I,k-1)
  1320. GAMMA = EL2ORC * QESO(I,k-1) / (TO(I,k-1)**2)
  1321. RFACT = 1. + DELTA * CP * GAMMA &
  1322. & * TO(I,k-1) / HVAP
  1323. AA1(I) = AA1(I) + &
  1324. & DZ1 * (G / (CP * TO(I,k-1))) &
  1325. & * DBYO(I,k-1) / (1. + GAMMA) &
  1326. & * RFACT
  1327. val = 0.
  1328. AA1(I)=AA1(I)+ &
  1329. & DZ1 * G * DELTA * &
  1330. !cmr & MAX( 0.,(QESO(I,k-1) - QO(I,k-1))) &
  1331. & MAX(val,(QESO(I,k-1) - QO(I,k-1)))
  1332. ENDIF
  1333. endif
  1334. ENDDO
  1335. ENDDO
  1336. DO I = 1, IM
  1337. IF(CNVFLG(I).AND.AA1(I).LE.0.) DWNFLG(I) = .FALSE.
  1338. IF(CNVFLG(I).AND.AA1(I).LE.0.) DWNFLG2(I) = .FALSE.
  1339. IF(CNVFLG(I).AND.AA1(I).LE.0.) CNVFLG(I) = .FALSE.
  1340. ENDDO
  1341. !!
  1342. TOTFLG = .TRUE.
  1343. DO I = 1, IM
  1344. TOTFLG = TOTFLG .AND. (.NOT. CNVFLG(I))
  1345. ENDDO
  1346. IF(TOTFLG) RETURN
  1347. !!
  1348. !cccc IF(LAT.EQ.LATD.AND.lon.eq.lond.and.CNVFLG(I)) THEN
  1349. !cccc PRINT *, ' AA1(I) BEFORE DWNDRFT =', AA1(I)
  1350. !cccc ENDIF
  1351. !
  1352. !------- DOWNDRAFT CALCULATIONS
  1353. !
  1354. !
  1355. !--- DETERMINE DOWNDRAFT STRENGTH IN TERMS OF WINDSHEAR
  1356. !
  1357. DO I = 1, IM
  1358. IF(CNVFLG(I)) THEN
  1359. VSHEAR(I) = 0.
  1360. ENDIF
  1361. ENDDO
  1362. DO K = 1, KM
  1363. DO I = 1, IM
  1364. if (k .le. kmax(i)) then
  1365. IF(K.GE.KB(I).AND.K.LE.KTCON(I).AND.CNVFLG(I)) THEN
  1366. shear=rcs(I) * sqrt((UO(I,k+1)-UO(I,k)) ** 2 &
  1367. & + (VO(I,k+1)-VO(I,k)) ** 2)
  1368. VSHEAR(I) = VSHEAR(I) + SHEAR
  1369. ENDIF
  1370. endif
  1371. ENDDO
  1372. ENDDO
  1373. DO I = 1, IM
  1374. EDT(I) = 0.
  1375. IF(CNVFLG(I)) THEN
  1376. KNUMB = KTCON(I) - KB(I) + 1
  1377. KNUMB = MAX(KNUMB,1)
  1378. VSHEAR(I) = 1.E3 * VSHEAR(I) / (ZO(I,KTCON(I))-ZO(I,KB(I)))
  1379. E1=1.591-.639*VSHEAR(I) &
  1380. & +.0953*(VSHEAR(I)**2)-.00496*(VSHEAR(I)**3)
  1381. EDT(I)=1.-E1
  1382. !cmr EDT(I) = MIN(EDT(I),.9)
  1383. val = .9
  1384. EDT(I) = MIN(EDT(I),val)
  1385. !cmr EDT(I) = MAX(EDT(I),.0)
  1386. val = .0
  1387. EDT(I) = MAX(EDT(I),val)
  1388. EDTO(I)=EDT(I)
  1389. EDTX(I)=EDT(I)
  1390. ENDIF
  1391. ENDDO
  1392. ! DETERMINE DETRAINMENT RATE BETWEEN 1 AND KBDTR
  1393. DO I = 1, IM
  1394. KBDTR(I) = KBCON(I)
  1395. beta = betas
  1396. if(SLIMSK(I).eq.1.) beta = betal
  1397. IF(CNVFLG(I)) THEN
  1398. KBDTR(I) = KBCON(I)
  1399. KBDTR(I) = MAX(KBDTR(I),1)
  1400. XLAMD(I) = 0.
  1401. IF(KBDTR(I).GT.1) THEN
  1402. DZ = .5 * ZO(I,KBDTR(I)) + .5 * ZO(I,KBDTR(I)-1) &
  1403. & - ZO(I,1)
  1404. XLAMD(I) = LOG(BETA) / DZ
  1405. ENDIF
  1406. ENDIF
  1407. ENDDO
  1408. ! DETERMINE DOWNDRAFT MASS FLUX
  1409. DO K = 1, KM
  1410. DO I = 1, IM
  1411. IF(k .le. kmax(i)) then
  1412. IF(CNVFLG(I)) THEN
  1413. ETAD(I,k) = 1.
  1414. ENDIF
  1415. QRCDO(I,k) = 0.
  1416. endif
  1417. ENDDO
  1418. ENDDO
  1419. DO K = KM1, 2, -1
  1420. DO I = 1, IM
  1421. if (k .le. kbmax(i)) then
  1422. IF(CNVFLG(I).AND.K.LT.KBDTR(I)) THEN
  1423. DZ = .5 * (ZO(I,k+1) - ZO(I,k-1))
  1424. ETAD(I,k) = ETAD(I,k+1) * EXP(XLAMD(I) * DZ)
  1425. ENDIF
  1426. endif
  1427. ENDDO
  1428. ENDDO
  1429. K = 1
  1430. DO I = 1, IM
  1431. IF(CNVFLG(I).AND.KBDTR(I).GT.1) THEN
  1432. DZ = .5 * (ZO(I,2) - ZO(I,1))
  1433. ETAD(I,k) = ETAD(I,k+1) * EXP(XLAMD(I) * DZ)
  1434. ENDIF
  1435. ENDDO
  1436. !
  1437. !--- DOWNDRAFT MOISTURE PROPERTIES
  1438. !
  1439. DO I = 1, IM
  1440. PWEVO(I) = 0.
  1441. FLG(I) = CNVFLG(I)
  1442. ENDDO
  1443. DO I = 1, IM
  1444. IF(CNVFLG(I)) THEN
  1445. JMN = JMIN(I)
  1446. HCDO(I) = HEO(I,JMN)
  1447. QCDO(I) = QO(I,JMN)
  1448. QRCDO(I,JMN) = QESO(I,JMN)
  1449. UCDO(I) = UO(I,JMN)
  1450. VCDO(I) = VO(I,JMN)
  1451. ENDIF
  1452. ENDDO
  1453. DO K = KM1, 1, -1
  1454. DO I = 1, IM
  1455. if (k .le. kmax(i)-1) then
  1456. IF(CNVFLG(I).AND.K.LT.JMIN(I)) THEN
  1457. DQ = QESO(I,k)
  1458. DT = TO(I,k)
  1459. GAMMA = EL2ORC * DQ / DT**2
  1460. DH = HCDO(I) - HESO(I,k)
  1461. QRCDO(I,k) = DQ+(1./HVAP)*(GAMMA/(1.+GAMMA))*DH
  1462. DETAD = ETAD(I,k+1) - ETAD(I,k)
  1463. PWDO(I,k) = ETAD(I,k+1) * QCDO(I) - &
  1464. & ETAD(I,k) * QRCDO(I,k)
  1465. PWDO(I,k) = PWDO(I,k) - DETAD * &
  1466. & .5 * (QRCDO(I,k) + QRCDO(I,k+1))
  1467. QCDO(I) = QRCDO(I,k)
  1468. PWEVO(I) = PWEVO(I) + PWDO(I,k)
  1469. ENDIF
  1470. endif
  1471. ENDDO
  1472. ENDDO
  1473. ! IF(LAT.EQ.LATD.AND.lon.eq.lond.and.DWNFLG(I)) THEN
  1474. ! PRINT *, ' PWAVO(I), PWEVO(I) =', PWAVO(I), PWEVO(I)
  1475. ! ENDIF
  1476. !
  1477. !--- FINAL DOWNDRAFT STRENGTH DEPENDENT ON PRECIP
  1478. !--- EFFICIENCY (EDT), NORMALIZED CONDENSATE (PWAV), AND
  1479. !--- EVAPORATE (PWEV)
  1480. !
  1481. DO I = 1, IM
  1482. edtmax = edtmaxl
  1483. if(SLIMSK(I).eq.0.) edtmax = edtmaxs
  1484. IF(DWNFLG2(I)) THEN
  1485. IF(PWEVO(I).LT.0.) THEN
  1486. EDTO(I) = -EDTO(I) * PWAVO(I) / PWEVO(I)
  1487. EDTO(I) = MIN(EDTO(I),EDTMAX)
  1488. ELSE
  1489. EDTO(I) = 0.
  1490. ENDIF
  1491. ELSE
  1492. EDTO(I) = 0.
  1493. ENDIF
  1494. ENDDO
  1495. !
  1496. !
  1497. !--- DOWNDRAFT CLOUDWORK FUNCTIONS
  1498. !
  1499. !
  1500. DO K = KM1, 1, -1
  1501. DO I = 1, IM
  1502. if (k .le. kmax(i)-1) then
  1503. IF(DWNFLG2(I).AND.K.LT.JMIN(I)) THEN
  1504. GAMMA = EL2ORC * QESO(I,k+1) / TO(I,k+1)**2
  1505. DHH=HCDO(I)
  1506. DT=TO(I,k+1)
  1507. DG=GAMMA
  1508. DH=HESO(I,k+1)
  1509. DZ=-1.*(ZO(I,k+1)-ZO(I,k))
  1510. AA1(I)=AA1(I)+EDTO(I)*DZ*(G/(CP*DT))*((DHH-DH)/(1.+DG)) &
  1511. & *(1.+DELTA*CP*DG*DT/HVAP)
  1512. val=0.
  1513. AA1(I)=AA1(I)+EDTO(I)* &
  1514. !cmr & DZ*G*DELTA*MAX( 0.,(QESO(I,k+1)-QO(I,k+1))) &
  1515. & DZ*G*DELTA*MAX(val,(QESO(I,k+1)-QO(I,k+1)))
  1516. ENDIF
  1517. endif
  1518. ENDDO
  1519. ENDDO
  1520. !cccc IF(LAT.EQ.LATD.AND.lon.eq.lond.and.DWNFLG2(I)) THEN
  1521. !cccc PRINT *, ' AA1(I) AFTER DWNDRFT =', AA1(I)
  1522. !cccc ENDIF
  1523. DO I = 1, IM
  1524. IF(AA1(I).LE.0.) CNVFLG(I) = .FALSE.
  1525. IF(AA1(I).LE.0.) DWNFLG(I) = .FALSE.
  1526. IF(AA1(I).LE.0.) DWNFLG2(I) = .FALSE.
  1527. ENDDO
  1528. !!
  1529. TOTFLG = .TRUE.
  1530. DO I = 1, IM
  1531. TOTFLG = TOTFLG .AND. (.NOT. CNVFLG(I))
  1532. ENDDO
  1533. IF(TOTFLG) RETURN
  1534. !!
  1535. !
  1536. !
  1537. !--- WHAT WOULD THE CHANGE BE, THAT A CLOUD WITH UNIT MASS
  1538. !--- WILL DO TO THE ENVIRONMENT?
  1539. !
  1540. DO K = 1, KM
  1541. DO I = 1, IM
  1542. IF(k .le. kmax(i) .and. CNVFLG(I)) THEN
  1543. DELLAH(I,k) = 0.
  1544. DELLAQ(I,k) = 0.
  1545. DELLAU(I,k) = 0.
  1546. DELLAV(I,k) = 0.
  1547. ENDIF
  1548. ENDDO
  1549. ENDDO
  1550. DO I = 1, IM
  1551. IF(CNVFLG(I)) THEN
  1552. DP = 1000. * DEL(I,1)
  1553. DELLAH(I,1) = EDTO(I) * ETAD(I,1) * (HCDO(I) &
  1554. & - HEO(I,1)) * G / DP
  1555. DELLAQ(I,1) = EDTO(I) * ETAD(I,1) * (QCDO(I) &
  1556. & - QO(I,1)) * G / DP
  1557. DELLAU(I,1) = EDTO(I) * ETAD(I,1) * (UCDO(I) &
  1558. & - UO(I,1)) * G / DP
  1559. DELLAV(I,1) = EDTO(I) * ETAD(I,1) * (VCDO(I) &
  1560. & - VO(I,1)) * G / DP
  1561. ENDIF
  1562. ENDDO
  1563. !
  1564. !--- CHANGED DUE TO SUBSIDENCE AND ENTRAINMENT
  1565. !
  1566. DO K = 2, KM1
  1567. DO I = 1, IM
  1568. if (k .le. kmax(i)-1) then
  1569. IF(CNVFLG(I).AND.K.LT.KTCON(I)) THEN
  1570. AUP = 1.
  1571. IF(K.LE.KB(I)) AUP = 0.
  1572. ADW = 1.
  1573. IF(K.GT.JMIN(I)) ADW = 0.
  1574. DV1= HEO(I,k)
  1575. DV2 = .5 * (HEO(I,k) + HEO(I,k+1))
  1576. DV3= HEO(I,k-1)
  1577. DV1Q= QO(I,k)
  1578. DV2Q = .5 * (QO(I,k) + QO(I,k+1))
  1579. DV3Q= QO(I,k-1)
  1580. DV1U= UO(I,k)
  1581. DV2U = .5 * (UO(I,k) + UO(I,k+1))
  1582. DV3U= UO(I,k-1)
  1583. DV1V= VO(I,k)
  1584. DV2V = .5 * (VO(I,k) + VO(I,k+1))
  1585. DV3V= VO(I,k-1)
  1586. DP = 1000. * DEL(I,K)
  1587. DZ = .5 * (ZO(I,k+1) - ZO(I,k-1))
  1588. DETA = ETA(I,k) - ETA(I,k-1)
  1589. DETAD = ETAD(I,k) - ETAD(I,k-1)
  1590. DELLAH(I,k) = DELLAH(I,k) + &
  1591. & ((AUP * ETA(I,k) - ADW * EDTO(I) * ETAD(I,k)) * DV1 &
  1592. & - (AUP * ETA(I,k-1) - ADW * EDTO(I) * ETAD(I,k-1))* DV3 &
  1593. & - AUP * DETA * DV2 &
  1594. & + ADW * EDTO(I) * DETAD * HCDO(I)) * G / DP
  1595. DELLAQ(I,k) = DELLAQ(I,k) + &
  1596. & ((AUP * ETA(I,k) - ADW * EDTO(I) * ETAD(I,k)) * DV1Q &
  1597. & - (AUP * ETA(I,k-1) - ADW * EDTO(I) * ETAD(I,k-1))* DV3Q &
  1598. & - AUP * DETA * DV2Q &
  1599. & +ADW*EDTO(I)*DETAD*.5*(QRCDO(I,k)+QRCDO(I,k-1))) * G / DP
  1600. DELLAU(I,k) = DELLAU(I,k) + &
  1601. & ((AUP * ETA(I,k) - ADW * EDTO(I) * ETAD(I,k)) * DV1U &
  1602. & - (AUP * ETA(I,k-1) - ADW * EDTO(I) * ETAD(I,k-1))* DV3U &
  1603. & - AUP * DETA * DV2U &
  1604. & + ADW * EDTO(I) * DETAD * UCDO(I) &
  1605. & ) * G / DP
  1606. DELLAV(I,k) = DELLAV(I,k) + &
  1607. & ((AUP * ETA(I,k) - ADW * EDTO(I) * ETAD(I,k)) * DV1V &
  1608. & - (AUP * ETA(I,k-1) - ADW * EDTO(I) * ETAD(I,k-1))* DV3V &
  1609. & - AUP * DETA * DV2V &
  1610. & + ADW * EDTO(I) * DETAD * VCDO(I) &
  1611. & ) * G / DP
  1612. ENDIF
  1613. endif
  1614. ENDDO
  1615. ENDDO
  1616. !
  1617. !------- CLOUD TOP
  1618. !
  1619. DO I = 1, IM
  1620. IF(CNVFLG(I)) THEN
  1621. INDX = KTCON(I)
  1622. DP = 1000. * DEL(I,INDX)
  1623. DV1 = HEO(I,INDX-1)
  1624. DELLAH(I,INDX) = ETA(I,INDX-1) * &
  1625. & (HCKO(I,INDX-1) - DV1) * G / DP
  1626. DVQ1 = QO(I,INDX-1)
  1627. DELLAQ(I,INDX) = ETA(I,INDX-1) * &
  1628. & (QCKO(I,INDX-1) - DVQ1) * G / DP
  1629. DV1U = UO(I,INDX-1)
  1630. DELLAU(I,INDX) = ETA(I,INDX-1) * &
  1631. & (UCKO(I,INDX-1) - DV1U) * G / DP
  1632. DV1V = VO(I,INDX-1)
  1633. DELLAV(I,INDX) = ETA(I,INDX-1) * &
  1634. & (VCKO(I,INDX-1) - DV1V) * G / DP
  1635. !
  1636. ! cloud water
  1637. !
  1638. DELLAL(I) = ETA(I,INDX-1) * QLKO_KTCON(I) * g / dp
  1639. ENDIF
  1640. ENDDO
  1641. !
  1642. !------- FINAL CHANGED VARIABLE PER UNIT MASS FLUX
  1643. !
  1644. DO K = 1, KM
  1645. DO I = 1, IM
  1646. if (k .le. kmax(i)) then
  1647. IF(CNVFLG(I).and.k.gt.KTCON(I)) THEN
  1648. QO(I,k) = Q1(I,k)
  1649. TO(I,k) = T1(I,k)
  1650. UO(I,k) = U1(I,k)
  1651. VO(I,k) = V1(I,k)
  1652. ENDIF
  1653. IF(CNVFLG(I).AND.K.LE.KTCON(I)) THEN
  1654. QO(I,k) = DELLAQ(I,k) * MBDT + Q1(I,k)
  1655. DELLAT = (DELLAH(I,k) - HVAP * DELLAQ(I,k)) / CP
  1656. TO(I,k) = DELLAT * MBDT + T1(I,k)
  1657. !cmr QO(I,k) = max(QO(I,k),1.e-10)
  1658. val = 1.e-10
  1659. QO(I,k) = max(QO(I,k), val )
  1660. ENDIF
  1661. endif
  1662. ENDDO
  1663. ENDDO
  1664. !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  1665. !
  1666. !--- THE ABOVE CHANGED ENVIRONMENT IS NOW USED TO CALULATE THE
  1667. !--- EFFECT THE ARBITRARY CLOUD (WITH UNIT MASS FLUX)
  1668. !--- WOULD HAVE ON THE STABILITY,
  1669. !--- WHICH THEN IS USED TO CALCULATE THE REAL MASS FLUX,
  1670. !--- NECESSARY TO KEEP THIS CHANGE IN BALANCE WITH THE LARGE-SCALE
  1671. !--- DESTABILIZATION.
  1672. !
  1673. !--- ENVIRONMENTAL CONDITIONS AGAIN, FIRST HEIGHTS
  1674. !
  1675. DO K = 1, KM
  1676. DO I = 1, IM
  1677. IF(k .le. kmax(i) .and. CNVFLG(I)) THEN
  1678. !jfe QESO(I,k) = 10. * FPVS(TO(I,k))
  1679. !
  1680. QESO(I,k) = 0.01 * fpvs(TO(I,K)) ! fpvs is in Pa
  1681. !
  1682. QESO(I,k) = EPS * QESO(I,k) / (PFLD(I,k)+EPSM1*QESO(I,k))
  1683. !cmr QESO(I,k) = MAX(QESO(I,k),1.E-8)
  1684. val = 1.E-8
  1685. QESO(I,k) = MAX(QESO(I,k), val )
  1686. TVO(I,k) = TO(I,k) + DELTA * TO(I,k) * QO(I,k)
  1687. ENDIF
  1688. ENDDO
  1689. ENDDO
  1690. DO I = 1, IM
  1691. IF(CNVFLG(I)) THEN
  1692. XAA0(I) = 0.
  1693. XPWAV(I) = 0.
  1694. ENDIF
  1695. ENDDO
  1696. !
  1697. ! HYDROSTATIC HEIGHT ASSUME ZERO TERR
  1698. !
  1699. ! DO I = 1, IM
  1700. ! IF(CNVFLG(I)) THEN
  1701. ! DLNSIG = LOG(PRSL(I,1)/PS(I))
  1702. ! ZO(I,1) = TERR - DLNSIG * RD / G * TVO(I,1)
  1703. ! ENDIF
  1704. ! ENDDO
  1705. ! DO K = 2, KM
  1706. ! DO I = 1, IM
  1707. ! IF(k .le. kmax(i) .and. CNVFLG(I)) THEN
  1708. ! DLNSIG = LOG(PRSL(I,K) / PRSL(I,K-1))
  1709. ! ZO(I,k) = ZO(I,k-1) - DLNSIG * RD / G
  1710. ! & * .5 * (TVO(I,k) + TVO(I,k-1))
  1711. ! ENDIF
  1712. ! ENDDO
  1713. ! ENDDO
  1714. !
  1715. !--- MOIST STATIC ENERGY
  1716. !
  1717. DO K = 1, KM1
  1718. DO I = 1, IM
  1719. IF(k .le. kmax(i)-1 .and. CNVFLG(I)) THEN
  1720. DZ = .5 * (ZO(I,k+1) - ZO(I,k))
  1721. DP = .5 * (PFLD(I,k+1) - PFLD(I,k))
  1722. !jfe ES = 10. * FPVS(TO(I,k+1))
  1723. !
  1724. ES = 0.01 * fpvs(TO(I,K+1)) ! fpvs is in Pa
  1725. !
  1726. PPRIME = PFLD(I,k+1) + EPSM1 * ES
  1727. QS = EPS * ES / PPRIME
  1728. DQSDP = - QS / PPRIME
  1729. DESDT = ES * (FACT1 / TO(I,k+1) + FACT2 / (TO(I,k+1)**2))
  1730. DQSDT = QS * PFLD(I,k+1) * DESDT / (ES * PPRIME)
  1731. GAMMA = EL2ORC * QESO(I,k+1) / (TO(I,k+1)**2)
  1732. DT = (G * DZ + HVAP * DQSDP * DP) / (CP * (1. + GAMMA))
  1733. DQ = DQSDT * DT + DQSDP * DP
  1734. TO(I,k) = TO(I,k+1) + DT
  1735. QO(I,k) = QO(I,k+1) + DQ
  1736. PO(I,k) = .5 * (PFLD(I,k) + PFLD(I,k+1))
  1737. ENDIF
  1738. ENDDO
  1739. ENDDO
  1740. DO K = 1, KM1
  1741. DO I = 1, IM
  1742. IF(k .le. kmax(i)-1 .and. CNVFLG(I)) THEN
  1743. !jfe QESO(I,k) = 10. * FPVS(TO(I,k))
  1744. !
  1745. QESO(I,k) = 0.01 * fpvs(TO(I,K)) ! fpvs is in Pa
  1746. !
  1747. QESO(I,k) = EPS * QESO(I,k) / (PO(I,k) + EPSM1 * QESO(I,k))
  1748. !cmr QESO(I,k) = MAX(QESO(I,k),1.E-8)
  1749. val1 = 1.E-8
  1750. QESO(I,k) = MAX(QESO(I,k), val1)
  1751. !cmr QO(I,k) = max(QO(I,k),1.e-10)
  1752. val2 = 1.e-10
  1753. QO(I,k) = max(QO(I,k), val2 )
  1754. ! QO(I,k) = MIN(QO(I,k),QESO(I,k))
  1755. HEO(I,k) = .5 * G * (ZO(I,k) + ZO(I,k+1)) + &
  1756. & CP * TO(I,k) + HVAP * QO(I,k)
  1757. HESO(I,k) = .5 * G * (ZO(I,k) + ZO(I,k+1)) + &
  1758. & CP * TO(I,k) + HVAP * QESO(I,k)
  1759. ENDIF
  1760. ENDDO
  1761. ENDDO
  1762. DO I = 1, IM
  1763. k = kmax(i)
  1764. IF(CNVFLG(I)) THEN
  1765. HEO(I,k) = G * ZO(I,k) + CP * TO(I,k) + HVAP * QO(I,k)
  1766. HESO(I,k) = G * ZO(I,k) + CP * TO(I,k) + HVAP * QESO(I,k)
  1767. ! HEO(I,k) = MIN(HEO(I,k),HESO(I,k))
  1768. ENDIF
  1769. ENDDO
  1770. DO I = 1, IM
  1771. IF(CNVFLG(I)) THEN
  1772. INDX = KB(I)
  1773. XHKB(I) = HEO(I,INDX)
  1774. XQKB(I) = QO(I,INDX)
  1775. HCKO(I,INDX) = XHKB(I)
  1776. QCKO(I,INDX) = XQKB(I)
  1777. ENDIF
  1778. ENDDO
  1779. !
  1780. !
  1781. !**************************** STATIC CONTROL
  1782. !
  1783. !
  1784. !------- MOISTURE AND CLOUD WORK FUNCTIONS
  1785. !
  1786. DO K = 2, KM1
  1787. DO I = 1, IM
  1788. if (k .le. kmax(i)-1) then
  1789. ! IF(CNVFLG(I).AND.K.GT.KB(I).AND.K.LE.KBCON(I)) THEN
  1790. IF(CNVFLG(I).AND.K.GT.KB(I).AND.K.LE.KTCON(I)) THEN
  1791. FACTOR = ETA(I,k-1) / ETA(I,k)
  1792. ONEMF = 1. - FACTOR
  1793. HCKO(I,k) = FACTOR * HCKO(I,k-1) + ONEMF * &
  1794. & .5 * (HEO(I,k) + HEO(I,k+1))
  1795. ENDIF
  1796. ! IF(CNVFLG(I).AND.K.GT.KBCON(I)) THEN
  1797. ! HEO(I,k) = HEO(I,k-1)
  1798. ! ENDIF
  1799. endif
  1800. ENDDO
  1801. ENDDO
  1802. DO K = 2, KM1
  1803. DO I = 1, IM
  1804. if (k .le. kmax(i)-1) then
  1805. IF(CNVFLG(I).AND.K.GT.KB(I).AND.K.LT.KTCON(I)) THEN
  1806. DZ = .5 * (ZO(I,k+1) - ZO(I,k-1))
  1807. GAMMA = EL2ORC * QESO(I,k) / (TO(I,k)**2)
  1808. XDBY = HCKO(I,k) - HESO(I,k)
  1809. !cmr XDBY = MAX(XDBY,0.)
  1810. val = 0.
  1811. XDBY = MAX(XDBY,val)
  1812. XQRCH = QESO(I,k) &
  1813. & + GAMMA * XDBY / (HVAP * (1. + GAMMA))
  1814. FACTOR = ETA(I,k-1) / ETA(I,k)
  1815. ONEMF = 1. - FACTOR
  1816. QCKO(I,k) = FACTOR * QCKO(I,k-1) + ONEMF * &
  1817. & .5 * (QO(I,k) + QO(I,k+1))
  1818. DQ = ETA(I,k) * QCKO(I,k) - ETA(I,k) * XQRCH
  1819. IF(DQ.GT.0.) THEN
  1820. ETAH = .5 * (ETA(I,k) + ETA(I,k-1))
  1821. QLK = DQ / (ETA(I,k) + ETAH * C0 * DZ)
  1822. XAA0(I) = XAA0(I) - (ZO(I,k) - ZO(I,k-1)) * G * QLK
  1823. XQC = QLK + XQRCH
  1824. XPW = ETAH * C0 * DZ * QLK
  1825. QCKO(I,k) = XQC
  1826. XPWAV(I) = XPWAV(I) + XPW
  1827. ENDIF
  1828. ENDIF
  1829. ! IF(CNVFLG(I).AND.K.GT.KBCON(I).AND.K.LT.KTCON(I)) THEN
  1830. IF(CNVFLG(I).AND.K.GT.KBCON(I).AND.K.LE.KTCON(I)) THEN
  1831. DZ1 = ZO(I,k) - ZO(I,k-1)
  1832. GAMMA = EL2ORC * QESO(I,k-1) / (TO(I,k-1)**2)
  1833. RFACT = 1. + DELTA * CP * GAMMA &
  1834. & * TO(I,k-1) / HVAP
  1835. XDBY = HCKO(I,k-1) - HESO(I,k-1)
  1836. XAA0(I) = XAA0(I) &
  1837. & + DZ1 * (G / (CP * TO(I,k-1))) &
  1838. & * XDBY / (1. + GAMMA) &
  1839. & * RFACT
  1840. val=0.
  1841. XAA0(I)=XAA0(I)+ &
  1842. & DZ1 * G * DELTA * &
  1843. !cmr & MAX( 0.,(QESO(I,k-1) - QO(I,k-1))) &
  1844. & MAX(val,(QESO(I,k-1) - QO(I,k-1)))
  1845. ENDIF
  1846. endif
  1847. ENDDO
  1848. ENDDO
  1849. !cccc IF(LAT.EQ.LATD.AND.lon.eq.lond.and.CNVFLG(I)) THEN
  1850. !cccc PRINT *, ' XAA BEFORE DWNDRFT =', XAA0(I)
  1851. !cccc ENDIF
  1852. !
  1853. !------- DOWNDRAFT CALCULATIONS
  1854. !
  1855. !
  1856. !--- DOWNDRAFT MOISTURE PROPERTIES
  1857. !
  1858. DO I = 1, IM
  1859. XPWEV(I) = 0.
  1860. ENDDO
  1861. DO I = 1, IM
  1862. IF(DWNFLG2(I)) THEN
  1863. JMN = JMIN(I)
  1864. XHCD(I) = HEO(I,JMN)
  1865. XQCD(I) = QO(I,JMN)
  1866. QRCD(I,JMN) = QESO(I,JMN)
  1867. ENDIF
  1868. ENDDO
  1869. DO K = KM1, 1, -1
  1870. DO I = 1, IM
  1871. if (k .le. kmax(i)-1) then
  1872. IF(DWNFLG2(I).AND.K.LT.JMIN(I)) THEN
  1873. DQ = QESO(I,k)
  1874. DT = TO(I,k)
  1875. GAMMA = EL2ORC * DQ / DT**2
  1876. DH = XHCD(I) - HESO(I,k)
  1877. QRCD(I,k)=DQ+(1./HVAP)*(GAMMA/(1.+GAMMA))*DH
  1878. DETAD = ETAD(I,k+1) - ETAD(I,k)
  1879. XPWD = ETAD(I,k+1) * QRCD(I,k+1) - &
  1880. & ETAD(I,k) * QRCD(I,k)
  1881. XPWD = XPWD - DETAD * &
  1882. & .5 * (QRCD(I,k) + QRCD(I,k+1))
  1883. XPWEV(I) = XPWEV(I) + XPWD
  1884. ENDIF
  1885. endif
  1886. ENDDO
  1887. ENDDO
  1888. !
  1889. DO I = 1, IM
  1890. edtmax = edtmaxl
  1891. if(SLIMSK(I).eq.0.) edtmax = edtmaxs
  1892. IF(DWNFLG2(I)) THEN
  1893. IF(XPWEV(I).GE.0.) THEN
  1894. EDTX(I) = 0.
  1895. ELSE
  1896. EDTX(I) = -EDTX(I) * XPWAV(I) / XPWEV(I)
  1897. EDTX(I) = MIN(EDTX(I),EDTMAX)
  1898. ENDIF
  1899. ELSE
  1900. EDTX(I) = 0.
  1901. ENDIF
  1902. ENDDO
  1903. !
  1904. !
  1905. !
  1906. !--- DOWNDRAFT CLOUDWORK FUNCTIONS
  1907. !
  1908. !
  1909. DO K = KM1, 1, -1
  1910. DO I = 1, IM
  1911. if (k .le. kmax(i)-1) then
  1912. IF(DWNFLG2(I).AND.K.LT.JMIN(I)) THEN
  1913. GAMMA = EL2ORC * QESO(I,k+1) / TO(I,k+1)**2
  1914. DHH=XHCD(I)
  1915. DT= TO(I,k+1)
  1916. DG= GAMMA
  1917. DH= HESO(I,k+1)
  1918. DZ=-1.*(ZO(I,k+1)-ZO(I,k))
  1919. XAA0(I)=XAA0(I)+EDTX(I)*DZ*(G/(CP*DT))*((DHH-DH)/(1.+DG)) &
  1920. & *(1.+DELTA*CP*DG*DT/HVAP)
  1921. val=0.
  1922. XAA0(I)=XAA0(I)+EDTX(I)* &
  1923. !cmr & DZ*G*DELTA*MAX( 0.,(QESO(I,k+1)-QO(I,k+1))) &
  1924. & DZ*G*DELTA*MAX(val,(QESO(I,k+1)-QO(I,k+1)))
  1925. ENDIF
  1926. endif
  1927. ENDDO
  1928. ENDDO
  1929. !cccc IF(LAT.EQ.LATD.AND.lon.eq.lond.and.DWNFLG2(I)) THEN
  1930. !cccc PRINT *, ' XAA AFTER DWNDRFT =', XAA0(I)
  1931. !cccc ENDIF
  1932. !
  1933. ! CALCULATE CRITICAL CLOUD WORK FUNCTION
  1934. !
  1935. DO I = 1, IM
  1936. ACRT(I) = 0.
  1937. IF(CNVFLG(I)) THEN
  1938. ! IF(CNVFLG(I).AND.SLIMSK(I).NE.1.) THEN
  1939. IF(PFLD(I,KTCON(I)).LT.PCRIT(15))THEN
  1940. ACRT(I)=ACRIT(15)*(975.-PFLD(I,KTCON(I))) &
  1941. & /(975.-PCRIT(15))
  1942. ELSE IF(PFLD(I,KTCON(I)).GT.PCRIT(1))THEN
  1943. ACRT(I)=ACRIT(1)
  1944. ELSE
  1945. !cmr K = IFIX((850. - PFLD(I,KTCON(I)))/50.) + 2
  1946. K = int((850. - PFLD(I,KTCON(I)))/50.) + 2
  1947. K = MIN(K,15)
  1948. K = MAX(K,2)
  1949. ACRT(I)=ACRIT(K)+(ACRIT(K-1)-ACRIT(K))* &
  1950. & (PFLD(I,KTCON(I))-PCRIT(K))/(PCRIT(K-1)-PCRIT(K))
  1951. ENDIF
  1952. ! ELSE
  1953. ! ACRT(I) = .5 * (PFLD(I,KBCON(I)) - PFLD(I,KTCON(I)))
  1954. ENDIF
  1955. ENDDO
  1956. DO I = 1, IM
  1957. ACRTFCT(I) = 1.
  1958. IF(CNVFLG(I)) THEN
  1959. if(SLIMSK(I).eq.1.) THEN
  1960. w1 = w1l
  1961. w2 = w2l
  1962. w3 = w3l
  1963. w4 = w4l
  1964. else
  1965. w1 = w1s
  1966. w2 = w2s
  1967. w3 = w3s
  1968. w4 = w4s
  1969. ENDIF
  1970. !C IF(CNVFLG(I).AND.SLIMSK(I).EQ.1.) THEN
  1971. ! ACRTFCT(I) = PDOT(I) / W3
  1972. !
  1973. ! modify critical cloud workfunction by cloud base vertical velocity
  1974. !
  1975. IF(PDOT(I).LE.W4) THEN
  1976. ACRTFCT(I) = (PDOT(I) - W4) / (W3 - W4)
  1977. ELSEIF(PDOT(I).GE.-W4) THEN
  1978. ACRTFCT(I) = - (PDOT(I) + W4) / (W4 - W3)
  1979. ELSE
  1980. ACRTFCT(I) = 0.
  1981. ENDIF
  1982. !cmr ACRTFCT(I) = MAX(ACRTFCT(I),-1.)
  1983. val1 = -1.
  1984. ACRTFCT(I) = MAX(ACRTFCT(I),val1)
  1985. !cmr ACRTFCT(I) = MIN(ACRTFCT(I),1.)
  1986. val2 = 1.
  1987. ACRTFCT(I) = MIN(ACRTFCT(I),val2)
  1988. ACRTFCT(I) = 1. - ACRTFCT(I)
  1989. !
  1990. ! modify ACRTFCT(I) by colume mean rh if RHBAR(I) is greater than 80 percent
  1991. !
  1992. ! if(RHBAR(I).ge..8) THEN
  1993. ! ACRTFCT(I) = ACRTFCT(I) * (.9 - min(RHBAR(I),.9)) * 10.
  1994. ! ENDIF
  1995. !
  1996. ! modify adjustment time scale by cloud base vertical velocity
  1997. !
  1998. DTCONV(I) = DT2 + max((1800. - DT2),RZERO) * &
  1999. & (PDOT(I) - W2) / (W1 - W2)
  2000. ! DTCONV(I) = MAX(DTCONV(I), DT2)
  2001. ! DTCONV(I) = 1800. * (PDOT(I) - w2) / (w1 - w2)
  2002. DTCONV(I) = max(DTCONV(I),dtmin)
  2003. DTCONV(I) = min(DTCONV(I),dtmax)
  2004. ENDIF
  2005. ENDDO
  2006. !
  2007. !--- LARGE SCALE FORCING
  2008. !
  2009. DO I= 1, IM
  2010. FLG(I) = CNVFLG(I)
  2011. IF(CNVFLG(I)) THEN
  2012. ! F = AA1(I) / DTCONV(I)
  2013. FLD(I) = (AA1(I) - ACRT(I) * ACRTFCT(I)) / DTCONV(I)
  2014. IF(FLD(I).LE.0.) FLG(I) = .FALSE.
  2015. ENDIF
  2016. CNVFLG(I) = FLG(I)
  2017. IF(CNVFLG(I)) THEN
  2018. ! XAA0(I) = MAX(XAA0(I),0.)
  2019. XK(I) = (XAA0(I) - AA1(I)) / MBDT
  2020. IF(XK(I).GE.0.) FLG(I) = .FALSE.
  2021. ENDIF
  2022. !
  2023. !--- KERNEL, CLOUD BASE MASS FLUX
  2024. !
  2025. CNVFLG(I) = FLG(I)
  2026. IF(CNVFLG(I)) THEN
  2027. XMB(I) = -FLD(I) / XK(I)
  2028. XMB(I) = MIN(XMB(I),XMBMAX(I))
  2029. ENDIF
  2030. ENDDO
  2031. ! IF(LAT.EQ.LATD.AND.lon.eq.lond.and.CNVFLG(I)) THEN
  2032. ! print *, ' RHBAR(I), ACRTFCT(I) =', RHBAR(I), ACRTFCT(I)
  2033. ! PRINT *, ' A1, XA =', AA1(I), XAA0(I)
  2034. ! PRINT *, ' XMB(I), ACRT =', XMB(I), ACRT
  2035. ! ENDIF
  2036. TOTFLG = .TRUE.
  2037. DO I = 1, IM
  2038. TOTFLG = TOTFLG .AND. (.NOT. CNVFLG(I))
  2039. ENDDO
  2040. IF(TOTFLG) RETURN
  2041. !
  2042. ! restore t0 and QO to t1 and q1 in case convection stops
  2043. !
  2044. do k = 1, km
  2045. DO I = 1, IM
  2046. if (k .le. kmax(i)) then
  2047. TO(I,k) = T1(I,k)
  2048. QO(I,k) = Q1(I,k)
  2049. !jfe QESO(I,k) = 10. * FPVS(T1(I,k))
  2050. !
  2051. QESO(I,k) = 0.01 * fpvs(T1(I,K)) ! fpvs is in Pa
  2052. !
  2053. QESO(I,k) = EPS * QESO(I,k) / (PFLD(I,k) + EPSM1*QESO(I,k))
  2054. !cmr QESO(I,k) = MAX(QESO(I,k),1.E-8)
  2055. val = 1.E-8
  2056. QESO(I,k) = MAX(QESO(I,k), val )
  2057. endif
  2058. enddo
  2059. enddo
  2060. !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  2061. !
  2062. !--- FEEDBACK: SIMPLY THE CHANGES FROM THE CLOUD WITH UNIT MASS FLUX
  2063. !--- MULTIPLIED BY THE MASS FLUX NECESSARY TO KEEP THE
  2064. !--- EQUILIBRIUM WITH THE LARGER-SCALE.
  2065. !
  2066. DO I = 1, IM
  2067. DELHBAR(I) = 0.
  2068. DELQBAR(I) = 0.
  2069. DELTBAR(I) = 0.
  2070. QCOND(I) = 0.
  2071. ENDDO
  2072. DO K = 1, KM
  2073. DO I = 1, IM
  2074. if (k .le. kmax(i)) then
  2075. IF(CNVFLG(I).AND.K.LE.KTCON(I)) THEN
  2076. AUP = 1.
  2077. IF(K.Le.KB(I)) AUP = 0.
  2078. ADW = 1.
  2079. IF(K.GT.JMIN(I)) ADW = 0.
  2080. DELLAT = (DELLAH(I,k) - HVAP * DELLAQ(I,k)) / CP
  2081. T1(I,k) = T1(I,k) + DELLAT * XMB(I) * DT2
  2082. Q1(I,k) = Q1(I,k) + DELLAQ(I,k) * XMB(I) * DT2
  2083. U1(I,k) = U1(I,k) + DELLAU(I,k) * XMB(I) * DT2
  2084. V1(I,k) = V1(I,k) + DELLAV(I,k) * XMB(I) * DT2
  2085. DP = 1000. * DEL(I,K)
  2086. DELHBAR(I) = DELHBAR(I) + DELLAH(I,k)*XMB(I)*DP/G
  2087. DELQBAR(I) = DELQBAR(I) + DELLAQ(I,k)*XMB(I)*DP/G
  2088. DELTBAR(I) = DELTBAR(I) + DELLAT*XMB(I)*DP/G
  2089. ENDIF
  2090. endif
  2091. ENDDO
  2092. ENDDO
  2093. DO K = 1, KM
  2094. DO I = 1, IM
  2095. if (k .le. kmax(i)) then
  2096. IF(CNVFLG(I).AND.K.LE.KTCON(I)) THEN
  2097. !jfe QESO(I,k) = 10. * FPVS(T1(I,k))
  2098. !
  2099. QESO(I,k) = 0.01 * fpvs(T1(I,K)) ! fpvs is in Pa
  2100. !
  2101. QESO(I,k) = EPS * QESO(I,k)/(PFLD(I,k) + EPSM1*QESO(I,k))
  2102. !cmr QESO(I,k) = MAX(QESO(I,k),1.E-8)
  2103. val = 1.E-8
  2104. QESO(I,k) = MAX(QESO(I,k), val )
  2105. !
  2106. ! cloud water
  2107. !
  2108. if(ncloud.gt.0.and.CNVFLG(I).and.k.eq.KTCON(I)) THEN
  2109. tem = DELLAL(I) * XMB(I) * dt2
  2110. tem1 = MAX(RZERO, MIN(RONE, (TCR-t1(I,K))*TCRF))
  2111. if (QL(I,k,2) .gt. -999.0) then
  2112. QL(I,k,1) = QL(I,k,1) + tem * tem1 ! Ice
  2113. QL(I,k,2) = QL(I,k,2) + tem *(1.0-tem1) ! Water
  2114. else
  2115. tem2 = QL(I,k,1) + tem
  2116. QL(I,k,1) = tem2 * tem1 ! Ice
  2117. QL(I,k,2) = tem2 - QL(I,k,1) ! Water
  2118. endif
  2119. ! QL(I,k) = QL(I,k) + DELLAL(I) * XMB(I) * dt2
  2120. dp = 1000. * del(i,k)
  2121. DELLAL(I) = DELLAL(I) * XMB(I) * dp / g
  2122. ENDIF
  2123. ENDIF
  2124. endif
  2125. ENDDO
  2126. ENDDO
  2127. ! IF(LAT.EQ.LATD.AND.lon.eq.lond.and.CNVFLG(I) ) THEN
  2128. ! PRINT *, ' DELHBAR, DELQBAR, DELTBAR ='
  2129. ! PRINT *, DELHBAR, HVAP*DELQBAR, CP*DELTBAR
  2130. ! PRINT *, ' DELLBAR ='
  2131. ! PRINT 6003, HVAP*DELLbar
  2132. ! PRINT *, ' DELLAQ ='
  2133. ! PRINT 6003, (HVAP*DELLAQ(I,k)*XMB(I),K=1,KMAX)
  2134. ! PRINT *, ' DELLAT ='
  2135. ! PRINT 6003, (DELLAH(i,k)*XMB(I)-HVAP*DELLAQ(I,k)*XMB(I), &
  2136. ! & K=1,KMAX)
  2137. ! ENDIF
  2138. DO I = 1, IM
  2139. RNTOT(I) = 0.
  2140. DELQEV(I) = 0.
  2141. DELQ2(I) = 0.
  2142. FLG(I) = CNVFLG(I)
  2143. ENDDO
  2144. DO K = KM, 1, -1
  2145. DO I = 1, IM
  2146. if (k .le. kmax(i)) then
  2147. IF(CNVFLG(I).AND.K.LE.KTCON(I)) THEN
  2148. AUP = 1.
  2149. IF(K.Le.KB(I)) AUP = 0.
  2150. ADW = 1.
  2151. IF(K.GT.JMIN(I)) ADW = 0.
  2152. rain = AUP * PWO(I,k) + ADW * EDTO(I) * PWDO(I,k)
  2153. RNTOT(I) = RNTOT(I) + rain * XMB(I) * .001 * dt2
  2154. ENDIF
  2155. endif
  2156. ENDDO
  2157. ENDDO
  2158. DO K = KM, 1, -1
  2159. DO I = 1, IM
  2160. if (k .le. kmax(i)) then
  2161. DELTV(I) = 0.
  2162. DELQ(I) = 0.
  2163. QEVAP(I) = 0.
  2164. IF(CNVFLG(I).AND.K.LE.KTCON(I)) THEN
  2165. AUP = 1.
  2166. IF(K.Le.KB(I)) AUP = 0.
  2167. ADW = 1.
  2168. IF(K.GT.JMIN(I)) ADW = 0.
  2169. rain = AUP * PWO(I,k) + ADW * EDTO(I) * PWDO(I,k)
  2170. RN(I) = RN(I) + rain * XMB(I) * .001 * dt2
  2171. ENDIF
  2172. IF(FLG(I).AND.K.LE.KTCON(I)) THEN
  2173. evef = EDT(I) * evfact
  2174. if(SLIMSK(I).eq.1.) evef=EDT(I) * evfactl
  2175. ! if(SLIMSK(I).eq.1.) evef=.07
  2176. ! if(SLIMSK(I).ne.1.) evef = 0.
  2177. QCOND(I) = EVEF * (Q1(I,k) - QESO(I,k)) &
  2178. & / (1. + EL2ORC * QESO(I,k) / T1(I,k)**2)
  2179. DP = 1000. * DEL(I,K)
  2180. IF(RN(I).GT.0..AND.QCOND(I).LT.0.) THEN
  2181. QEVAP(I) = -QCOND(I) * (1.-EXP(-.32*SQRT(DT2*RN(I))))
  2182. QEVAP(I) = MIN(QEVAP(I), RN(I)*1000.*G/DP)
  2183. DELQ2(I) = DELQEV(I) + .001 * QEVAP(I) * dp / g
  2184. ENDIF
  2185. if(RN(I).gt.0..and.QCOND(I).LT.0..and. &
  2186. & DELQ2(I).gt.RNTOT(I)) THEN
  2187. QEVAP(I) = 1000.* g * (RNTOT(I) - DELQEV(I)) / dp
  2188. FLG(I) = .false.
  2189. ENDIF
  2190. IF(RN(I).GT.0..AND.QEVAP(I).gt.0.) THEN
  2191. Q1(I,k) = Q1(I,k) + QEVAP(I)
  2192. T1(I,k) = T1(I,k) - ELOCP * QEVAP(I)
  2193. RN(I) = RN(I) - .001 * QEVAP(I) * DP / G
  2194. DELTV(I) = - ELOCP*QEVAP(I)/DT2
  2195. DELQ(I) = + QEVAP(I)/DT2
  2196. DELQEV(I) = DELQEV(I) + .001*dp*QEVAP(I)/g
  2197. ENDIF
  2198. DELLAQ(I,k) = DELLAQ(I,k) + DELQ(I) / XMB(I)
  2199. DELQBAR(I) = DELQBAR(I) + DELQ(I)*DP/G
  2200. DELTBAR(I) = DELTBAR(I) + DELTV(I)*DP/G
  2201. ENDIF
  2202. endif
  2203. ENDDO
  2204. ENDDO
  2205. ! IF(LAT.EQ.LATD.AND.lon.eq.lond.and.CNVFLG(I) ) THEN
  2206. ! PRINT *, ' DELLAH ='
  2207. ! PRINT 6003, (DELLAH(k)*XMB(I),K=1,KMAX)
  2208. ! PRINT *, ' DELLAQ ='
  2209. ! PRINT 6003, (HVAP*DELLAQ(I,k)*XMB(I),K=1,KMAX)
  2210. ! PRINT *, ' DELHBAR, DELQBAR, DELTBAR ='
  2211. ! PRINT *, DELHBAR, HVAP*DELQBAR, CP*DELTBAR
  2212. ! PRINT *, ' PRECIP =', HVAP*RN(I)*1000./DT2
  2213. !CCCC PRINT *, ' DELLBAR ='
  2214. !CCCC PRINT *, HVAP*DELLbar
  2215. ! ENDIF
  2216. !
  2217. ! PRECIPITATION RATE CONVERTED TO ACTUAL PRECIP
  2218. ! IN UNIT OF M INSTEAD OF KG
  2219. !
  2220. DO I = 1, IM
  2221. IF(CNVFLG(I)) THEN
  2222. !
  2223. ! IN THE EVENT OF UPPER LEVEL RAIN EVAPORATION AND LOWER LEVEL DOWNDRAF
  2224. ! MOISTENING, RN CAN BECOME NEGATIVE, IN THIS CASE, WE BACK OUT OF TH
  2225. ! HEATING AND THE MOISTENING
  2226. !
  2227. if(RN(I).lt.0..and..not.FLG(I)) RN(I) = 0.
  2228. IF(RN(I).LE.0.) THEN
  2229. RN(I) = 0.
  2230. ELSE
  2231. KTOP(I) = KTCON(I)
  2232. KBOT(I) = KBCON(I)
  2233. KUO(I) = 1
  2234. CLDWRK(I) = AA1(I)
  2235. ENDIF
  2236. ENDIF
  2237. ENDDO
  2238. DO K = 1, KM
  2239. DO I = 1, IM
  2240. if (k .le. kmax(i)) then
  2241. IF(CNVFLG(I).AND.RN(I).LE.0.) THEN
  2242. T1(I,k) = TO(I,k)
  2243. Q1(I,k) = QO(I,k)
  2244. ENDIF
  2245. endif
  2246. ENDDO
  2247. ENDDO
  2248. !!
  2249. RETURN
  2250. END SUBROUTINE SASCNV
  2251. ! ------------------------------------------------------------------------
  2252. SUBROUTINE OLD_ARW_SHALCV(IM,IX,KM,DT,DEL,PRSI,PRSL,PRSLK,KUO,Q,T,DPSHC)
  2253. !
  2254. USE MODULE_GFS_MACHINE , ONLY : kind_phys
  2255. USE MODULE_GFS_PHYSCONS, grav => con_g, CP => con_CP, HVAP => con_HVAP &
  2256. &, RD => con_RD
  2257. implicit none
  2258. !
  2259. ! include 'constant.h'
  2260. !
  2261. integer IM, IX, KM, KUO(IM)
  2262. real(kind=kind_phys) DEL(IX,KM), PRSI(IX,KM+1), PRSL(IX,KM), &
  2263. & PRSLK(IX,KM), &
  2264. & Q(IX,KM), T(IX,KM), DT, DPSHC
  2265. !
  2266. ! Locals
  2267. !
  2268. real(kind=kind_phys) ck, cpdt, dmse, dsdz1, dsdz2, &
  2269. & dsig, dtodsl, dtodsu, eldq, g, &
  2270. & gocp, rtdls
  2271. !
  2272. integer k,k1,k2,kliftl,kliftu,kt,N2,I,iku,ik1,ik,ii
  2273. integer INDEX2(IM), KLCL(IM), KBOT(IM), KTOP(IM),kk &
  2274. &, KTOPM(IM)
  2275. !!
  2276. ! PHYSICAL PARAMETERS
  2277. PARAMETER(G=GRAV, GOCP=G/CP)
  2278. ! BOUNDS OF PARCEL ORIGIN
  2279. PARAMETER(KLIFTL=2,KLIFTU=2)
  2280. LOGICAL LSHC(IM)
  2281. real(kind=kind_phys) Q2(IM*KM), T2(IM*KM), &
  2282. & PRSL2(IM*KM), PRSLK2(IM*KM), &
  2283. & AL(IM*(KM-1)), AD(IM*KM), AU(IM*(KM-1))
  2284. !-----------------------------------------------------------------------
  2285. ! COMPRESS FIELDS TO POINTS WITH NO DEEP CONVECTION
  2286. ! AND MOIST STATIC INSTABILITY.
  2287. DO I=1,IM
  2288. LSHC(I)=.FALSE.
  2289. ENDDO
  2290. DO K=1,KM-1
  2291. DO I=1,IM
  2292. IF(KUO(I).EQ.0) THEN
  2293. ELDQ = HVAP*(Q(I,K)-Q(I,K+1))
  2294. CPDT = CP*(T(I,K)-T(I,K+1))
  2295. RTDLS = (PRSL(I,K)-PRSL(I,K+1)) / &
  2296. & PRSI(I,K+1)*RD*0.5*(T(I,K)+T(I,K+1))
  2297. DMSE = ELDQ+CPDT-RTDLS
  2298. LSHC(I) = LSHC(I).OR.DMSE.GT.0.
  2299. ENDIF
  2300. ENDDO
  2301. ENDDO
  2302. N2 = 0
  2303. DO I=1,IM
  2304. IF(LSHC(I)) THEN
  2305. N2 = N2 + 1
  2306. INDEX2(N2) = I
  2307. ENDIF
  2308. ENDDO
  2309. IF(N2.EQ.0) RETURN
  2310. DO K=1,KM
  2311. KK = (K-1)*N2
  2312. DO I=1,N2
  2313. IK = KK + I
  2314. ii = index2(i)
  2315. Q2(IK) = Q(II,K)
  2316. T2(IK) = T(II,K)
  2317. PRSL2(IK) = PRSL(II,K)
  2318. PRSLK2(IK) = PRSLK(II,K)
  2319. ENDDO
  2320. ENDDO
  2321. do i=1,N2
  2322. ktopm(i) = KM
  2323. enddo
  2324. do k=2,KM
  2325. do i=1,N2
  2326. ii = index2(i)
  2327. if (prsi(ii,1)-prsi(ii,k) .le. dpshc) ktopm(i) = k
  2328. enddo
  2329. enddo
  2330. !-----------------------------------------------------------------------
  2331. ! COMPUTE MOIST ADIABAT AND DETERMINE LIMITS OF SHALLOW CONVECTION.
  2332. ! CHECK FOR MOIST STATIC INSTABILITY AGAIN WITHIN CLOUD.
  2333. CALL MSTADBT3(N2,KM-1,KLIFTL,KLIFTU,PRSL2,PRSLK2,T2,Q2, &
  2334. & KLCL,KBOT,KTOP,AL,AU)
  2335. DO I=1,N2
  2336. KBOT(I) = min(KLCL(I)-1, ktopm(i)-1)
  2337. KTOP(I) = min(KTOP(I)+1, ktopm(i))
  2338. LSHC(I) = .FALSE.
  2339. ENDDO
  2340. DO K=1,KM-1
  2341. KK = (K-1)*N2
  2342. DO I=1,N2
  2343. IF(K.GE.KBOT(I).AND.K.LT.KTOP(I)) THEN
  2344. IK = KK + I
  2345. IKU = IK + N2
  2346. ELDQ = HVAP * (Q2(IK)-Q2(IKU))
  2347. CPDT = CP * (T2(IK)-T2(IKU))
  2348. RTDLS = (PRSL2(IK)-PRSL2(IKU)) / &
  2349. & PRSI(index2(i),K+1)*RD*0.5*(T2(IK)+T2(IKU))
  2350. DMSE = ELDQ + CPDT - RTDLS
  2351. LSHC(I) = LSHC(I).OR.DMSE.GT.0.
  2352. AU(IK) = G/RTDLS
  2353. ENDIF
  2354. ENDDO
  2355. ENDDO
  2356. K1=KM+1
  2357. K2=0
  2358. DO I=1,N2
  2359. IF(.NOT.LSHC(I)) THEN
  2360. KBOT(I) = KM+1
  2361. KTOP(I) = 0
  2362. ENDIF
  2363. K1 = MIN(K1,KBOT(I))
  2364. K2 = MAX(K2,KTOP(I))
  2365. ENDDO
  2366. KT = K2-K1+1
  2367. IF(KT.LT.2) RETURN
  2368. !-----------------------------------------------------------------------
  2369. ! SET EDDY VISCOSITY COEFFICIENT CKU AT SIGMA INTERFACES.
  2370. ! COMPUTE DIAGONALS AND RHS FOR TRIDIAGONAL MATRIX SOLVER.
  2371. ! EXPAND FINAL FIELDS.
  2372. KK = (K1-1) * N2
  2373. DO I=1,N2
  2374. IK = KK + I
  2375. AD(IK) = 1.
  2376. ENDDO
  2377. !
  2378. ! DTODSU=DT/DEL(K1)
  2379. DO K=K1,K2-1
  2380. ! DTODSL=DTODSU
  2381. ! DTODSU= DT/DEL(K+1)
  2382. ! DSIG=SL(K)-SL(K+1)
  2383. KK = (K-1) * N2
  2384. DO I=1,N2
  2385. ii = index2(i)
  2386. DTODSL = DT/DEL(II,K)
  2387. DTODSU = DT/DEL(II,K+1)
  2388. DSIG = PRSL(II,K) - PRSL(II,K+1)
  2389. IK = KK + I
  2390. IKU = IK + N2
  2391. IF(K.EQ.KBOT(I)) THEN
  2392. CK=1.5
  2393. ELSEIF(K.EQ.KTOP(I)-1) THEN
  2394. CK=1.
  2395. ELSEIF(K.EQ.KTOP(I)-2) THEN
  2396. CK=3.
  2397. ELSEIF(K.GT.KBOT(I).AND.K.LT.KTOP(I)-2) THEN
  2398. CK=5.
  2399. ELSE
  2400. CK=0.
  2401. ENDIF
  2402. DSDZ1 = CK*DSIG*AU(IK)*GOCP
  2403. DSDZ2 = CK*DSIG*AU(IK)*AU(IK)
  2404. AU(IK) = -DTODSL*DSDZ2
  2405. AL(IK) = -DTODSU*DSDZ2
  2406. AD(IK) = AD(IK)-AU(IK)
  2407. AD(IKU) = 1.-AL(IK)
  2408. T2(IK) = T2(IK)+DTODSL*DSDZ1
  2409. T2(IKU) = T2(IKU)-DTODSU*DSDZ1
  2410. ENDDO
  2411. ENDDO
  2412. IK1=(K1-1)*N2+1
  2413. CALL TRIDI2T3(N2,KT,AL(IK1),AD(IK1),AU(IK1),Q2(IK1),T2(IK1), &
  2414. & AU(IK1),Q2(IK1),T2(IK1))
  2415. DO K=K1,K2
  2416. KK = (K-1)*N2
  2417. DO I=1,N2
  2418. IK = KK + I
  2419. Q(INDEX2(I),K) = Q2(IK)
  2420. T(INDEX2(I),K) = T2(IK)
  2421. ENDDO
  2422. ENDDO
  2423. !-----------------------------------------------------------------------
  2424. RETURN
  2425. END SUBROUTINE OLD_ARW_SHALCV
  2426. !-----------------------------------------------------------------------
  2427. SUBROUTINE TRIDI2T3(L,N,CL,CM,CU,R1,R2,AU,A1,A2)
  2428. !yt INCLUDE DBTRIDI2;
  2429. !!
  2430. USE MODULE_GFS_MACHINE , ONLY : kind_phys
  2431. implicit none
  2432. integer k,n,l,i
  2433. real(kind=kind_phys) fk
  2434. !!
  2435. real(kind=kind_phys) &
  2436. & CL(L,2:N),CM(L,N),CU(L,N-1),R1(L,N),R2(L,N), &
  2437. & AU(L,N-1),A1(L,N),A2(L,N)
  2438. !-----------------------------------------------------------------------
  2439. DO I=1,L
  2440. FK=1./CM(I,1)
  2441. AU(I,1)=FK*CU(I,1)
  2442. A1(I,1)=FK*R1(I,1)
  2443. A2(I,1)=FK*R2(I,1)
  2444. ENDDO
  2445. DO K=2,N-1
  2446. DO I=1,L
  2447. FK=1./(CM(I,K)-CL(I,K)*AU(I,K-1))
  2448. AU(I,K)=FK*CU(I,K)
  2449. A1(I,K)=FK*(R1(I,K)-CL(I,K)*A1(I,K-1))
  2450. A2(I,K)=FK*(R2(I,K)-CL(I,K)*A2(I,K-1))
  2451. ENDDO
  2452. ENDDO
  2453. DO I=1,L
  2454. FK=1./(CM(I,N)-CL(I,N)*AU(I,N-1))
  2455. A1(I,N)=FK*(R1(I,N)-CL(I,N)*A1(I,N-1))
  2456. A2(I,N)=FK*(R2(I,N)-CL(I,N)*A2(I,N-1))
  2457. ENDDO
  2458. DO K=N-1,1,-1
  2459. DO I=1,L
  2460. A1(I,K)=A1(I,K)-AU(I,K)*A1(I,K+1)
  2461. A2(I,K)=A2(I,K)-AU(I,K)*A2(I,K+1)
  2462. ENDDO
  2463. ENDDO
  2464. !-----------------------------------------------------------------------
  2465. RETURN
  2466. END SUBROUTINE TRIDI2T3
  2467. !-----------------------------------------------------------------------
  2468. SUBROUTINE MSTADBT3(IM,KM,K1,K2,PRSL,PRSLK,TENV,QENV, &
  2469. & KLCL,KBOT,KTOP,TCLD,QCLD)
  2470. !yt INCLUDE DBMSTADB;
  2471. !!
  2472. USE MODULE_GFS_MACHINE, ONLY : kind_phys
  2473. USE MODULE_GFS_FUNCPHYS, ONLY : FTDP, FTHE, FTLCL, STMA
  2474. USE MODULE_GFS_PHYSCONS, EPS => con_eps, EPSM1 => con_epsm1, FV => con_FVirt
  2475. implicit none
  2476. !!
  2477. ! include 'constant.h'
  2478. !!
  2479. integer k,k1,k2,km,i,im
  2480. real(kind=kind_phys) pv,qma,slklcl,tdpd,thelcl,tlcl
  2481. real(kind=kind_phys) tma,tvcld,tvenv
  2482. !!
  2483. real(kind=kind_phys) PRSL(IM,KM), PRSLK(IM,KM), TENV(IM,KM), &
  2484. & QENV(IM,KM), TCLD(IM,KM), QCLD(IM,KM)
  2485. INTEGER KLCL(IM), KBOT(IM), KTOP(IM)
  2486. ! LOCAL ARRAYS
  2487. real(kind=kind_phys) SLKMA(IM), THEMA(IM)
  2488. !-----------------------------------------------------------------------
  2489. ! DETERMINE WARMEST POTENTIAL WET-BULB TEMPERATURE BETWEEN K1 AND K2.
  2490. ! COMPUTE ITS LIFTING CONDENSATION LEVEL.
  2491. !
  2492. DO I=1,IM
  2493. SLKMA(I) = 0.
  2494. THEMA(I) = 0.
  2495. ENDDO
  2496. DO K=K1,K2
  2497. DO I=1,IM
  2498. PV = 1000.0 * PRSL(I,K)*QENV(I,K)/(EPS-EPSM1*QENV(I,K))
  2499. TDPD = TENV(I,K)-FTDP(PV)
  2500. IF(TDPD.GT.0.) THEN
  2501. TLCL = FTLCL(TENV(I,K),TDPD)
  2502. SLKLCL = PRSLK(I,K)*TLCL/TENV(I,K)
  2503. ELSE
  2504. TLCL = TENV(I,K)
  2505. SLKLCL = PRSLK(I,K)
  2506. ENDIF
  2507. THELCL=FTHE(TLCL,SLKLCL)
  2508. IF(THELCL.GT.THEMA(I)) THEN
  2509. SLKMA(I) = SLKLCL
  2510. THEMA(I) = THELCL
  2511. ENDIF
  2512. ENDDO
  2513. ENDDO
  2514. !-----------------------------------------------------------------------
  2515. ! SET CLOUD TEMPERATURES AND HUMIDITIES WHEREVER THE PARCEL LIFTED UP
  2516. ! THE MOIST ADIABAT IS BUOYANT WITH RESPECT TO THE ENVIRONMENT.
  2517. DO I=1,IM
  2518. KLCL(I)=KM+1
  2519. KBOT(I)=KM+1
  2520. KTOP(I)=0
  2521. ENDDO
  2522. DO K=1,KM
  2523. DO I=1,IM
  2524. TCLD(I,K)=0.
  2525. QCLD(I,K)=0.
  2526. ENDDO
  2527. ENDDO
  2528. DO K=K1,KM
  2529. DO I=1,IM
  2530. IF(PRSLK(I,K).LE.SLKMA(I)) THEN
  2531. KLCL(I)=MIN(KLCL(I),K)
  2532. CALL STMA(THEMA(I),PRSLK(I,K),TMA,QMA)
  2533. ! TMA=FTMA(THEMA(I),PRSLK(I,K),QMA)
  2534. TVCLD=TMA*(1.+FV*QMA)
  2535. TVENV=TENV(I,K)*(1.+FV*QENV(I,K))
  2536. IF(TVCLD.GT.TVENV) THEN
  2537. KBOT(I)=MIN(KBOT(I),K)
  2538. KTOP(I)=MAX(KTOP(I),K)
  2539. TCLD(I,K)=TMA-TENV(I,K)
  2540. QCLD(I,K)=QMA-QENV(I,K)
  2541. ENDIF
  2542. ENDIF
  2543. ENDDO
  2544. ENDDO
  2545. !-----------------------------------------------------------------------
  2546. RETURN
  2547. END SUBROUTINE MSTADBT3
  2548. subroutine sascnvn(im,ix,km,jcap,delt,del,prsl,ps,phil,ql, &
  2549. & q1,t1,u1,v1,rcs,cldwrk,rn,kbot,ktop,kcnv,slimsk, &
  2550. & dot,ncloud,pgcon,sas_mass_flux)
  2551. ! & dot,ncloud,ud_mf,dd_mf,dt_mf)
  2552. ! & dot,ncloud,ud_mf,dd_mf,dt_mf,me)
  2553. !
  2554. ! use machine , only : kind_phys
  2555. ! use funcphys , only : fpvs
  2556. ! use physcons, grav => con_g, cp => con_cp, hvap => con_hvap &
  2557. USE MODULE_GFS_MACHINE, ONLY : kind_phys
  2558. USE MODULE_GFS_FUNCPHYS, ONLY : fpvs
  2559. USE MODULE_GFS_PHYSCONS, grav => con_g, cp => con_cp &
  2560. &, hvap => con_hvap &
  2561. &, rv => con_rv, fv => con_fvirt, t0c => con_t0c &
  2562. &, cvap => con_cvap, cliq => con_cliq &
  2563. &, eps => con_eps, epsm1 => con_epsm1
  2564. implicit none
  2565. !
  2566. integer im, ix, km, jcap, ncloud, &
  2567. & kbot(im), ktop(im), kcnv(im)
  2568. ! &, me
  2569. real(kind=kind_phys) delt,sas_mass_flux
  2570. real(kind=kind_phys) ps(im), del(ix,km), prsl(ix,km), &
  2571. & ql(ix,km,2),q1(ix,km), t1(ix,km), &
  2572. & u1(ix,km), v1(ix,km), rcs(im), &
  2573. & cldwrk(im), rn(im), slimsk(im), &
  2574. & dot(ix,km), phil(ix,km)
  2575. ! hchuang code change mass flux output
  2576. ! &, ud_mf(im,km),dd_mf(im,km),dt_mf(im,km)
  2577. !
  2578. integer i, j, indx, jmn, k, kk, latd, lond, km1
  2579. !
  2580. real(kind=kind_phys) clam, cxlamu, xlamde, xlamdd
  2581. !
  2582. real(kind=kind_phys) adw, aup, aafac, &
  2583. & beta, betal, betas, &
  2584. & c0, cpoel, dellat, delta, &
  2585. & desdt, deta, detad, dg, &
  2586. & dh, dhh, dlnsig, dp, &
  2587. & dq, dqsdp, dqsdt, dt, &
  2588. & dt2, dtmax, dtmin, dv1h, &
  2589. & dv1q, dv2h, dv2q, dv1u, &
  2590. & dv1v, dv2u, dv2v, dv3q, &
  2591. & dv3h, dv3u, dv3v, &
  2592. & dz, dz1, e1, edtmax, &
  2593. & edtmaxl, edtmaxs, el2orc, elocp, &
  2594. & es, etah, cthk, dthk, &
  2595. & evef, evfact, evfactl, fact1, &
  2596. & fact2, factor, fjcap, fkm, &
  2597. & g, gamma, pprime, &
  2598. & qlk, qrch, qs, c1, &
  2599. & rain, rfact, shear, tem1, &
  2600. & tem2, terr, val, val1, &
  2601. & val2, w1, w1l, w1s, &
  2602. & w2, w2l, w2s, w3, &
  2603. & w3l, w3s, w4, w4l, &
  2604. & w4s, xdby, xpw, xpwd, &
  2605. & xqrch, mbdt, tem, &
  2606. & ptem, ptem1
  2607. !
  2608. real(kind=kind_phys), intent(in) :: pgcon
  2609. integer kb(im), kbcon(im), kbcon1(im), &
  2610. & ktcon(im), ktcon1(im), &
  2611. & jmin(im), lmin(im), kbmax(im), &
  2612. & kbm(im), kmax(im)
  2613. !
  2614. real(kind=kind_phys) aa1(im), acrt(im), acrtfct(im), &
  2615. & delhbar(im), delq(im), delq2(im), &
  2616. & delqbar(im), delqev(im), deltbar(im), &
  2617. & deltv(im), dtconv(im), edt(im), &
  2618. & edto(im), edtx(im), fld(im), &
  2619. & hcdo(im,km), hmax(im), hmin(im), &
  2620. & ucdo(im,km), vcdo(im,km),aa2(im), &
  2621. & pbcdif(im), pdot(im), po(im,km), &
  2622. & pwavo(im), pwevo(im), xlamud(im), &
  2623. & qcdo(im,km), qcond(im), qevap(im), &
  2624. & rntot(im), vshear(im), xaa0(im), &
  2625. & xk(im), xlamd(im), &
  2626. & xmb(im), xmbmax(im), xpwav(im), &
  2627. & xpwev(im), delubar(im),delvbar(im)
  2628. !cj
  2629. real(kind=kind_phys) cincr, cincrmax, cincrmin
  2630. real(kind=kind_phys) xmbmx1
  2631. !cj
  2632. !c physical parameters
  2633. parameter(g=grav)
  2634. parameter(cpoel=cp/hvap,elocp=hvap/cp, &
  2635. & el2orc=hvap*hvap/(rv*cp))
  2636. parameter(terr=0.,c0=.002,c1=.002,delta=fv)
  2637. parameter(fact1=(cvap-cliq)/rv,fact2=hvap/rv-fact1*t0c)
  2638. parameter(cthk=150.,cincrmax=180.,cincrmin=120.,dthk=25.)
  2639. !c local variables and arrays
  2640. real(kind=kind_phys) pfld(im,km),to(im,km), qo(im,km), &
  2641. & uo(im,km), vo(im,km), qeso(im,km)
  2642. !c cloud water
  2643. real(kind=kind_phys)qlko_ktcon(im),dellal(im,km),tvo(im,km), &
  2644. & dbyo(im,km), zo(im,km), xlamue(im,km), &
  2645. & fent1(im,km),fent2(im,km), frh(im,km), &
  2646. & heo(im,km), heso(im,km), &
  2647. & qrcd(im,km), dellah(im,km), dellaq(im,km), &
  2648. & dellau(im,km),dellav(im,km), hcko(im,km), &
  2649. & ucko(im,km), vcko(im,km), qcko(im,km), &
  2650. & eta(im,km), etad(im,km), zi(im,km), &
  2651. & qrcdo(im,km),pwo(im,km), pwdo(im,km), &
  2652. & tx1(im), sumx(im)
  2653. ! &, rhbar(im)
  2654. !
  2655. logical totflg, cnvflg(im), flg(im)
  2656. !
  2657. real(kind=kind_phys) pcrit(15), acritt(15), acrit(15)
  2658. ! save pcrit, acritt
  2659. data pcrit/850.,800.,750.,700.,650.,600.,550.,500.,450.,400.,&
  2660. & 350.,300.,250.,200.,150./
  2661. data acritt/.0633,.0445,.0553,.0664,.075,.1082,.1521,.2216, &
  2662. & .3151,.3677,.41,.5255,.7663,1.1686,1.6851/
  2663. !c gdas derived acrit
  2664. !c data acritt/.203,.515,.521,.566,.625,.665,.659,.688,
  2665. !c & .743,.813,.886,.947,1.138,1.377,1.896/
  2666. real(kind=kind_phys) tf, tcr, tcrf
  2667. parameter (tf=233.16, tcr=263.16, tcrf=1.0/(tcr-tf))
  2668. !
  2669. !c-----------------------------------------------------------------------
  2670. !
  2671. km1 = km - 1
  2672. !c
  2673. !c initialize arrays
  2674. !c
  2675. do i=1,im
  2676. kcnv(i)=0
  2677. cnvflg(i) = .true.
  2678. rn(i)=0.
  2679. kbot(i)=km+1
  2680. ktop(i)=0
  2681. kbcon(i)=km
  2682. ktcon(i)=1
  2683. dtconv(i) = 3600.
  2684. cldwrk(i) = 0.
  2685. pdot(i) = 0.
  2686. pbcdif(i)= 0.
  2687. lmin(i) = 1
  2688. jmin(i) = 1
  2689. qlko_ktcon(i) = 0.
  2690. edt(i) = 0.
  2691. edto(i) = 0.
  2692. edtx(i) = 0.
  2693. acrt(i) = 0.
  2694. acrtfct(i) = 1.
  2695. aa1(i) = 0.
  2696. aa2(i) = 0.
  2697. xaa0(i) = 0.
  2698. pwavo(i)= 0.
  2699. pwevo(i)= 0.
  2700. xpwav(i)= 0.
  2701. xpwev(i)= 0.
  2702. vshear(i) = 0.
  2703. enddo
  2704. ! hchuang code change
  2705. ! do k = 1, km
  2706. ! do i = 1, im
  2707. ! ud_mf(i,k) = 0.
  2708. ! dd_mf(i,k) = 0.
  2709. ! dt_mf(i,k) = 0.
  2710. ! enddo
  2711. ! enddo
  2712. !c
  2713. do k = 1, 15
  2714. acrit(k) = acritt(k) * (975. - pcrit(k))
  2715. enddo
  2716. dt2 = delt
  2717. val = 1200.
  2718. dtmin = max(dt2, val )
  2719. val = 3600.
  2720. dtmax = max(dt2, val )
  2721. !c model tunable parameters are all here
  2722. mbdt = 10.
  2723. edtmaxl = .3
  2724. edtmaxs = .3
  2725. clam = .1
  2726. aafac = .1
  2727. ! betal = .15
  2728. ! betas = .15
  2729. betal = .05
  2730. betas = .05
  2731. !c evef = 0.07
  2732. evfact = 0.3
  2733. evfactl = 0.3
  2734. #if ( EM_CORE == 1 )
  2735. ! HAWAII TEST - ZCX
  2736. BETAl = .05
  2737. betas = .05
  2738. evfact = 0.5
  2739. evfactl = 0.5
  2740. #endif
  2741. !
  2742. cxlamu = 1.0e-4
  2743. xlamde = 1.0e-4
  2744. xlamdd = 1.0e-4
  2745. !
  2746. fjcap = (float(jcap) / 126.) ** 2
  2747. val = 1.
  2748. fjcap = max(fjcap,val)
  2749. fkm = (float(km) / 28.) ** 2
  2750. fkm = max(fkm,val)
  2751. w1l = -8.e-3
  2752. w2l = -4.e-2
  2753. w3l = -5.e-3
  2754. w4l = -5.e-4
  2755. w1s = -2.e-4
  2756. w2s = -2.e-3
  2757. w3s = -1.e-3
  2758. w4s = -2.e-5
  2759. !c
  2760. !c define top layer for search of the downdraft originating layer
  2761. !c and the maximum thetae for updraft
  2762. !c
  2763. do i=1,im
  2764. kbmax(i) = km
  2765. kbm(i) = km
  2766. kmax(i) = km
  2767. tx1(i) = 1.0 / ps(i)
  2768. enddo
  2769. !
  2770. do k = 1, km
  2771. do i=1,im
  2772. IF (prSL(I,K)*tx1(I) .GT. 0.04) KMAX(I) = MIN(KM,K + 1)
  2773. !2011bugfix if (prsl(i,k)*tx1(i) .gt. 0.04) kmax(i) = k + 1
  2774. if (prsl(i,k)*tx1(i) .gt. 0.45) kbmax(i) = k + 1
  2775. if (prsl(i,k)*tx1(i) .gt. 0.70) kbm(i) = k + 1
  2776. enddo
  2777. enddo
  2778. do i=1,im
  2779. kbmax(i) = min(kbmax(i),kmax(i))
  2780. kbm(i) = min(kbm(i),kmax(i))
  2781. enddo
  2782. !c
  2783. !c hydrostatic height assume zero terr and initially assume
  2784. !c updraft entrainment rate as an inverse function of height
  2785. !c
  2786. do k = 1, km
  2787. do i=1,im
  2788. zo(i,k) = phil(i,k) / g
  2789. enddo
  2790. enddo
  2791. do k = 1, km1
  2792. do i=1,im
  2793. zi(i,k) = 0.5*(zo(i,k)+zo(i,k+1))
  2794. xlamue(i,k) = clam / zi(i,k)
  2795. enddo
  2796. enddo
  2797. !c
  2798. !c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  2799. !c convert surface pressure to mb from cb
  2800. !c
  2801. do k = 1, km
  2802. do i = 1, im
  2803. if (k .le. kmax(i)) then
  2804. pfld(i,k) = prsl(i,k) * 10.0
  2805. eta(i,k) = 1.
  2806. fent1(i,k)= 1.
  2807. fent2(i,k)= 1.
  2808. frh(i,k) = 0.
  2809. hcko(i,k) = 0.
  2810. qcko(i,k) = 0.
  2811. ucko(i,k) = 0.
  2812. vcko(i,k) = 0.
  2813. etad(i,k) = 1.
  2814. hcdo(i,k) = 0.
  2815. qcdo(i,k) = 0.
  2816. ucdo(i,k) = 0.
  2817. vcdo(i,k) = 0.
  2818. qrcd(i,k) = 0.
  2819. qrcdo(i,k)= 0.
  2820. dbyo(i,k) = 0.
  2821. pwo(i,k) = 0.
  2822. pwdo(i,k) = 0.
  2823. dellal(i,k) = 0.
  2824. to(i,k) = t1(i,k)
  2825. qo(i,k) = q1(i,k)
  2826. uo(i,k) = u1(i,k) * rcs(i)
  2827. vo(i,k) = v1(i,k) * rcs(i)
  2828. endif
  2829. enddo
  2830. enddo
  2831. !c
  2832. !c column variables
  2833. !c p is pressure of the layer (mb)
  2834. !c t is temperature at t-dt (k)..tn
  2835. !c q is mixing ratio at t-dt (kg/kg)..qn
  2836. !c to is temperature at t+dt (k)... this is after advection and turbulan
  2837. !c qo is mixing ratio at t+dt (kg/kg)..q1
  2838. !c
  2839. do k = 1, km
  2840. do i=1,im
  2841. if (k .le. kmax(i)) then
  2842. qeso(i,k) = 0.01 * fpvs(to(i,k)) ! fpvs is in pa
  2843. qeso(i,k) = eps * qeso(i,k) / (pfld(i,k) + epsm1*qeso(i,k))
  2844. val1 = 1.e-8
  2845. qeso(i,k) = max(qeso(i,k), val1)
  2846. val2 = 1.e-10
  2847. qo(i,k) = max(qo(i,k), val2 )
  2848. ! qo(i,k) = min(qo(i,k),qeso(i,k))
  2849. ! tvo(i,k) = to(i,k) + delta * to(i,k) * qo(i,k)
  2850. endif
  2851. enddo
  2852. enddo
  2853. !c
  2854. !c compute moist static energy
  2855. !c
  2856. do k = 1, km
  2857. do i=1,im
  2858. if (k .le. kmax(i)) then
  2859. ! tem = g * zo(i,k) + cp * to(i,k)
  2860. tem = phil(i,k) + cp * to(i,k)
  2861. heo(i,k) = tem + hvap * qo(i,k)
  2862. heso(i,k) = tem + hvap * qeso(i,k)
  2863. !c heo(i,k) = min(heo(i,k),heso(i,k))
  2864. endif
  2865. enddo
  2866. enddo
  2867. !c
  2868. !c determine level with largest moist static energy
  2869. !c this is the level where updraft starts
  2870. !c
  2871. do i=1,im
  2872. hmax(i) = heo(i,1)
  2873. kb(i) = 1
  2874. enddo
  2875. do k = 2, km
  2876. do i=1,im
  2877. if (k .le. kbm(i)) then
  2878. if(heo(i,k).gt.hmax(i)) then
  2879. kb(i) = k
  2880. hmax(i) = heo(i,k)
  2881. endif
  2882. endif
  2883. enddo
  2884. enddo
  2885. !c
  2886. do k = 1, km1
  2887. do i=1,im
  2888. if (k .le. kmax(i)-1) then
  2889. dz = .5 * (zo(i,k+1) - zo(i,k))
  2890. dp = .5 * (pfld(i,k+1) - pfld(i,k))
  2891. es = 0.01 * fpvs(to(i,k+1)) ! fpvs is in pa
  2892. pprime = pfld(i,k+1) + epsm1 * es
  2893. qs = eps * es / pprime
  2894. dqsdp = - qs / pprime
  2895. desdt = es * (fact1 / to(i,k+1) + fact2 / (to(i,k+1)**2))
  2896. dqsdt = qs * pfld(i,k+1) * desdt / (es * pprime)
  2897. gamma = el2orc * qeso(i,k+1) / (to(i,k+1)**2)
  2898. dt = (g * dz + hvap * dqsdp * dp) / (cp * (1. + gamma))
  2899. dq = dqsdt * dt + dqsdp * dp
  2900. to(i,k) = to(i,k+1) + dt
  2901. qo(i,k) = qo(i,k+1) + dq
  2902. po(i,k) = .5 * (pfld(i,k) + pfld(i,k+1))
  2903. endif
  2904. enddo
  2905. enddo
  2906. !
  2907. do k = 1, km1
  2908. do i=1,im
  2909. if (k .le. kmax(i)-1) then
  2910. qeso(i,k) = 0.01 * fpvs(to(i,k)) ! fpvs is in pa
  2911. qeso(i,k) = eps * qeso(i,k) / (po(i,k) + epsm1*qeso(i,k))
  2912. val1 = 1.e-8
  2913. qeso(i,k) = max(qeso(i,k), val1)
  2914. val2 = 1.e-10
  2915. qo(i,k) = max(qo(i,k), val2 )
  2916. ! qo(i,k) = min(qo(i,k),qeso(i,k))
  2917. val1 = 1.0
  2918. frh(i,k) = 1. - min(qo(i,k)/qeso(i,k), val1)
  2919. heo(i,k) = .5 * g * (zo(i,k) + zo(i,k+1)) + &
  2920. & cp * to(i,k) + hvap * qo(i,k)
  2921. heso(i,k) = .5 * g * (zo(i,k) + zo(i,k+1)) + &
  2922. & cp * to(i,k) + hvap * qeso(i,k)
  2923. uo(i,k) = .5 * (uo(i,k) + uo(i,k+1))
  2924. vo(i,k) = .5 * (vo(i,k) + vo(i,k+1))
  2925. endif
  2926. enddo
  2927. enddo
  2928. !c
  2929. !c look for the level of free convection as cloud base
  2930. !c
  2931. do i=1,im
  2932. flg(i) = .true.
  2933. kbcon(i) = kmax(i)
  2934. enddo
  2935. do k = 1, km1
  2936. do i=1,im
  2937. if (flg(i).and.k.le.kbmax(i)) then
  2938. if(k.gt.kb(i).and.heo(i,kb(i)).gt.heso(i,k)) then
  2939. kbcon(i) = k
  2940. flg(i) = .false.
  2941. endif
  2942. endif
  2943. enddo
  2944. enddo
  2945. !c
  2946. do i=1,im
  2947. if(kbcon(i).eq.kmax(i)) cnvflg(i) = .false.
  2948. enddo
  2949. !!
  2950. totflg = .true.
  2951. do i=1,im
  2952. totflg = totflg .and. (.not. cnvflg(i))
  2953. enddo
  2954. if(totflg) return
  2955. !!
  2956. !c
  2957. !c determine critical convective inhibition
  2958. !c as a function of vertical velocity at cloud base.
  2959. !c
  2960. do i=1,im
  2961. if(cnvflg(i)) then
  2962. pdot(i) = 10.* dot(i,kbcon(i))
  2963. endif
  2964. enddo
  2965. do i=1,im
  2966. if(cnvflg(i)) then
  2967. if(slimsk(i).eq.1.) then
  2968. w1 = w1l
  2969. w2 = w2l
  2970. w3 = w3l
  2971. w4 = w4l
  2972. else
  2973. w1 = w1s
  2974. w2 = w2s
  2975. w3 = w3s
  2976. w4 = w4s
  2977. endif
  2978. if(pdot(i).le.w4) then
  2979. tem = (pdot(i) - w4) / (w3 - w4)
  2980. elseif(pdot(i).ge.-w4) then
  2981. tem = - (pdot(i) + w4) / (w4 - w3)
  2982. else
  2983. tem = 0.
  2984. endif
  2985. val1 = -1.
  2986. tem = max(tem,val1)
  2987. val2 = 1.
  2988. tem = min(tem,val2)
  2989. tem = 1. - tem
  2990. tem1= .5*(cincrmax-cincrmin)
  2991. cincr = cincrmax - tem * tem1
  2992. pbcdif(i) = pfld(i,kb(i)) - pfld(i,kbcon(i))
  2993. if(pbcdif(i).gt.cincr) then
  2994. cnvflg(i) = .false.
  2995. endif
  2996. endif
  2997. enddo
  2998. !!
  2999. totflg = .true.
  3000. do i=1,im
  3001. totflg = totflg .and. (.not. cnvflg(i))
  3002. enddo
  3003. if(totflg) return
  3004. !!
  3005. !c
  3006. !c assume that updraft entrainment rate above cloud base is
  3007. !c same as that at cloud base
  3008. !c
  3009. do k = 2, km1
  3010. do i=1,im
  3011. if(cnvflg(i).and. &
  3012. & (k.gt.kbcon(i).and.k.lt.kmax(i))) then
  3013. xlamue(i,k) = xlamue(i,kbcon(i))
  3014. endif
  3015. enddo
  3016. enddo
  3017. !c
  3018. !c assume the detrainment rate for the updrafts to be same as
  3019. !c the entrainment rate at cloud base
  3020. !c
  3021. do i = 1, im
  3022. if(cnvflg(i)) then
  3023. xlamud(i) = xlamue(i,kbcon(i))
  3024. endif
  3025. enddo
  3026. !c
  3027. !c functions rapidly decreasing with height, mimicking a cloud ensemble
  3028. !c (Bechtold et al., 2008)
  3029. !c
  3030. do k = 2, km1
  3031. do i=1,im
  3032. if(cnvflg(i).and. &
  3033. & (k.gt.kbcon(i).and.k.lt.kmax(i))) then
  3034. tem = qeso(i,k)/qeso(i,kbcon(i))
  3035. fent1(i,k) = tem**2
  3036. fent2(i,k) = tem**3
  3037. endif
  3038. enddo
  3039. enddo
  3040. !c
  3041. !c final entrainment rate as the sum of turbulent part and organized entrainment
  3042. !c depending on the environmental relative humidity
  3043. !c (Bechtold et al., 2008)
  3044. !c
  3045. do k = 2, km1
  3046. do i=1,im
  3047. if(cnvflg(i).and. &
  3048. & (k.ge.kbcon(i).and.k.lt.kmax(i))) then
  3049. tem = cxlamu * frh(i,k) * fent2(i,k)
  3050. xlamue(i,k) = xlamue(i,k)*fent1(i,k) + tem
  3051. endif
  3052. enddo
  3053. enddo
  3054. !c
  3055. !c determine updraft mass flux for the subcloud layers
  3056. !c
  3057. do k = km1, 1, -1
  3058. do i = 1, im
  3059. if (cnvflg(i)) then
  3060. if(k.lt.kbcon(i).and.k.ge.kb(i)) then
  3061. dz = zi(i,k+1) - zi(i,k)
  3062. ptem = 0.5*(xlamue(i,k)+xlamue(i,k+1))-xlamud(i)
  3063. eta(i,k) = eta(i,k+1) / (1. + ptem * dz)
  3064. endif
  3065. endif
  3066. enddo
  3067. enddo
  3068. !c
  3069. !c compute mass flux above cloud base
  3070. !c
  3071. do k = 2, km1
  3072. do i = 1, im
  3073. if(cnvflg(i))then
  3074. if(k.gt.kbcon(i).and.k.lt.kmax(i)) then
  3075. dz = zi(i,k) - zi(i,k-1)
  3076. ptem = 0.5*(xlamue(i,k)+xlamue(i,k-1))-xlamud(i)
  3077. eta(i,k) = eta(i,k-1) * (1 + ptem * dz)
  3078. endif
  3079. endif
  3080. enddo
  3081. enddo
  3082. !c
  3083. !c compute updraft cloud properties
  3084. !c
  3085. do i = 1, im
  3086. if(cnvflg(i)) then
  3087. indx = kb(i)
  3088. hcko(i,indx) = heo(i,indx)
  3089. ucko(i,indx) = uo(i,indx)
  3090. vcko(i,indx) = vo(i,indx)
  3091. pwavo(i) = 0.
  3092. endif
  3093. enddo
  3094. !c
  3095. !c cloud property is modified by the entrainment process
  3096. !c
  3097. do k = 2, km1
  3098. do i = 1, im
  3099. if (cnvflg(i)) then
  3100. if(k.gt.kb(i).and.k.lt.kmax(i)) then
  3101. dz = zi(i,k) - zi(i,k-1)
  3102. tem = 0.5 * (xlamue(i,k)+xlamue(i,k-1)) * dz
  3103. tem1 = 0.5 * xlamud(i) * dz
  3104. factor = 1. + tem - tem1
  3105. ptem = 0.5 * tem + pgcon
  3106. ptem1= 0.5 * tem - pgcon
  3107. hcko(i,k) = ((1.-tem1)*hcko(i,k-1)+tem*0.5* &
  3108. & (heo(i,k)+heo(i,k-1)))/factor
  3109. ucko(i,k) = ((1.-tem1)*ucko(i,k-1)+ptem*uo(i,k) &
  3110. & +ptem1*uo(i,k-1))/factor
  3111. vcko(i,k) = ((1.-tem1)*vcko(i,k-1)+ptem*vo(i,k) &
  3112. & +ptem1*vo(i,k-1))/factor
  3113. dbyo(i,k) = hcko(i,k) - heso(i,k)
  3114. endif
  3115. endif
  3116. enddo
  3117. enddo
  3118. !c
  3119. !c taking account into convection inhibition due to existence of
  3120. !c dry layers below cloud base
  3121. !c
  3122. do i=1,im
  3123. flg(i) = cnvflg(i)
  3124. kbcon1(i) = kmax(i)
  3125. enddo
  3126. do k = 2, km1
  3127. do i=1,im
  3128. if (flg(i).and.k.lt.kmax(i)) then
  3129. if(k.ge.kbcon(i).and.dbyo(i,k).gt.0.) then
  3130. kbcon1(i) = k
  3131. flg(i) = .false.
  3132. endif
  3133. endif
  3134. enddo
  3135. enddo
  3136. do i=1,im
  3137. if(cnvflg(i)) then
  3138. if(kbcon1(i).eq.kmax(i)) cnvflg(i) = .false.
  3139. endif
  3140. enddo
  3141. do i=1,im
  3142. if(cnvflg(i)) then
  3143. tem = pfld(i,kbcon(i)) - pfld(i,kbcon1(i))
  3144. if(tem.gt.dthk) then
  3145. cnvflg(i) = .false.
  3146. endif
  3147. endif
  3148. enddo
  3149. !!
  3150. totflg = .true.
  3151. do i = 1, im
  3152. totflg = totflg .and. (.not. cnvflg(i))
  3153. enddo
  3154. if(totflg) return
  3155. !!
  3156. !c
  3157. !c determine first guess cloud top as the level of zero buoyancy
  3158. !c
  3159. do i = 1, im
  3160. flg(i) = cnvflg(i)
  3161. ktcon(i) = 1
  3162. enddo
  3163. do k = 2, km1
  3164. do i = 1, im
  3165. if (flg(i).and.k .lt. kmax(i)) then
  3166. if(k.gt.kbcon1(i).and.dbyo(i,k).lt.0.) then
  3167. ktcon(i) = k
  3168. flg(i) = .false.
  3169. endif
  3170. endif
  3171. enddo
  3172. enddo
  3173. !c
  3174. do i = 1, im
  3175. if(cnvflg(i)) then
  3176. tem = pfld(i,kbcon(i))-pfld(i,ktcon(i))
  3177. if(tem.lt.cthk) cnvflg(i) = .false.
  3178. endif
  3179. enddo
  3180. !!
  3181. totflg = .true.
  3182. do i = 1, im
  3183. totflg = totflg .and. (.not. cnvflg(i))
  3184. enddo
  3185. if(totflg) return
  3186. !!
  3187. !c
  3188. !c search for downdraft originating level above theta-e minimum
  3189. !c
  3190. do i = 1, im
  3191. if(cnvflg(i)) then
  3192. hmin(i) = heo(i,kbcon1(i))
  3193. lmin(i) = kbmax(i)
  3194. jmin(i) = kbmax(i)
  3195. endif
  3196. enddo
  3197. do k = 2, km1
  3198. do i = 1, im
  3199. if (cnvflg(i) .and. k .le. kbmax(i)) then
  3200. if(k.gt.kbcon1(i).and.heo(i,k).lt.hmin(i)) then
  3201. lmin(i) = k + 1
  3202. hmin(i) = heo(i,k)
  3203. endif
  3204. endif
  3205. enddo
  3206. enddo
  3207. !c
  3208. !c make sure that jmin(i) is within the cloud
  3209. !c
  3210. do i = 1, im
  3211. if(cnvflg(i)) then
  3212. jmin(i) = min(lmin(i),ktcon(i)-1)
  3213. jmin(i) = max(jmin(i),kbcon1(i)+1)
  3214. if(jmin(i).ge.ktcon(i)) cnvflg(i) = .false.
  3215. endif
  3216. enddo
  3217. !c
  3218. !c specify upper limit of mass flux at cloud base
  3219. !c
  3220. do i = 1, im
  3221. if(cnvflg(i)) then
  3222. ! xmbmax(i) = .1
  3223. !
  3224. k = kbcon(i)
  3225. dp = 1000. * del(i,k)
  3226. xmbmax(i) = dp / (g * dt2)
  3227. xmbmax(i) = min(sas_mass_flux,xmbmax(i))
  3228. !
  3229. ! tem = dp / (g * dt2)
  3230. ! xmbmax(i) = min(tem, xmbmax(i))
  3231. endif
  3232. enddo
  3233. !c
  3234. !c compute cloud moisture property and precipitation
  3235. !c
  3236. do i = 1, im
  3237. if (cnvflg(i)) then
  3238. aa1(i) = 0.
  3239. qcko(i,kb(i)) = qo(i,kb(i))
  3240. ! rhbar(i) = 0.
  3241. endif
  3242. enddo
  3243. do k = 2, km1
  3244. do i = 1, im
  3245. if (cnvflg(i)) then
  3246. if(k.gt.kb(i).and.k.lt.ktcon(i)) then
  3247. dz = zi(i,k) - zi(i,k-1)
  3248. gamma = el2orc * qeso(i,k) / (to(i,k)**2)
  3249. qrch = qeso(i,k) &
  3250. & + gamma * dbyo(i,k) / (hvap * (1. + gamma))
  3251. !cj
  3252. tem = 0.5 * (xlamue(i,k)+xlamue(i,k-1)) * dz
  3253. tem1 = 0.5 * xlamud(i) * dz
  3254. factor = 1. + tem - tem1
  3255. qcko(i,k) = ((1.-tem1)*qcko(i,k-1)+tem*0.5* &
  3256. & (qo(i,k)+qo(i,k-1)))/factor
  3257. !cj
  3258. dq = eta(i,k) * (qcko(i,k) - qrch)
  3259. !c
  3260. ! rhbar(i) = rhbar(i) + qo(i,k) / qeso(i,k)
  3261. !c
  3262. !c check if there is excess moisture to release latent heat
  3263. !c
  3264. if(k.ge.kbcon(i).and.dq.gt.0.) then
  3265. etah = .5 * (eta(i,k) + eta(i,k-1))
  3266. if(ncloud.gt.0..and.k.gt.jmin(i)) then
  3267. dp = 1000. * del(i,k)
  3268. qlk = dq / (eta(i,k) + etah * (c0 + c1) * dz)
  3269. dellal(i,k) = etah * c1 * dz * qlk * g / dp
  3270. else
  3271. qlk = dq / (eta(i,k) + etah * c0 * dz)
  3272. endif
  3273. aa1(i) = aa1(i) - dz * g * qlk
  3274. qcko(i,k) = qlk + qrch
  3275. pwo(i,k) = etah * c0 * dz * qlk
  3276. pwavo(i) = pwavo(i) + pwo(i,k)
  3277. endif
  3278. endif
  3279. endif
  3280. enddo
  3281. enddo
  3282. !c
  3283. ! do i = 1, im
  3284. ! if(cnvflg(i)) then
  3285. ! indx = ktcon(i) - kb(i) - 1
  3286. ! rhbar(i) = rhbar(i) / float(indx)
  3287. ! endif
  3288. ! enddo
  3289. !c
  3290. !c calculate cloud work function
  3291. !c
  3292. do k = 2, km1
  3293. do i = 1, im
  3294. if (cnvflg(i)) then
  3295. if(k.ge.kbcon(i).and.k.lt.ktcon(i)) then
  3296. dz1 = zo(i,k+1) - zo(i,k)
  3297. gamma = el2orc * qeso(i,k) / (to(i,k)**2)
  3298. rfact = 1. + delta * cp * gamma &
  3299. & * to(i,k) / hvap
  3300. aa1(i) = aa1(i) + &
  3301. & dz1 * (g / (cp * to(i,k))) &
  3302. & * dbyo(i,k) / (1. + gamma) &
  3303. & * rfact
  3304. val = 0.
  3305. aa1(i)=aa1(i)+ &
  3306. & dz1 * g * delta * &
  3307. & max(val,(qeso(i,k) - qo(i,k)))
  3308. endif
  3309. endif
  3310. enddo
  3311. enddo
  3312. do i = 1, im
  3313. if(cnvflg(i).and.aa1(i).le.0.) cnvflg(i) = .false.
  3314. enddo
  3315. !!
  3316. totflg = .true.
  3317. do i=1,im
  3318. totflg = totflg .and. (.not. cnvflg(i))
  3319. enddo
  3320. if(totflg) return
  3321. !!
  3322. !c
  3323. !c estimate the onvective overshooting as the level
  3324. !c where the [aafac * cloud work function] becomes zero,
  3325. !c which is the final cloud top
  3326. !c
  3327. do i = 1, im
  3328. if (cnvflg(i)) then
  3329. aa2(i) = aafac * aa1(i)
  3330. endif
  3331. enddo
  3332. !c
  3333. do i = 1, im
  3334. flg(i) = cnvflg(i)
  3335. ktcon1(i) = kmax(i) - 1
  3336. enddo
  3337. do k = 2, km1
  3338. do i = 1, im
  3339. if (flg(i)) then
  3340. if(k.ge.ktcon(i).and.k.lt.kmax(i)) then
  3341. dz1 = zo(i,k+1) - zo(i,k)
  3342. gamma = el2orc * qeso(i,k) / (to(i,k)**2)
  3343. rfact = 1. + delta * cp * gamma &
  3344. & * to(i,k) / hvap
  3345. aa2(i) = aa2(i) + &
  3346. & dz1 * (g / (cp * to(i,k))) &
  3347. & * dbyo(i,k) / (1. + gamma) &
  3348. & * rfact
  3349. if(aa2(i).lt.0.) then
  3350. ktcon1(i) = k
  3351. flg(i) = .false.
  3352. endif
  3353. endif
  3354. endif
  3355. enddo
  3356. enddo
  3357. !c
  3358. !c compute cloud moisture property, detraining cloud water
  3359. !c and precipitation in overshooting layers
  3360. !c
  3361. do k = 2, km1
  3362. do i = 1, im
  3363. if (cnvflg(i)) then
  3364. if(k.ge.ktcon(i).and.k.lt.ktcon1(i)) then
  3365. dz = zi(i,k) - zi(i,k-1)
  3366. gamma = el2orc * qeso(i,k) / (to(i,k)**2)
  3367. qrch = qeso(i,k) &
  3368. & + gamma * dbyo(i,k) / (hvap * (1. + gamma))
  3369. !cj
  3370. tem = 0.5 * (xlamue(i,k)+xlamue(i,k-1)) * dz
  3371. tem1 = 0.5 * xlamud(i) * dz
  3372. factor = 1. + tem - tem1
  3373. qcko(i,k) = ((1.-tem1)*qcko(i,k-1)+tem*0.5* &
  3374. & (qo(i,k)+qo(i,k-1)))/factor
  3375. !cj
  3376. dq = eta(i,k) * (qcko(i,k) - qrch)
  3377. !c
  3378. !c check if there is excess moisture to release latent heat
  3379. !c
  3380. if(dq.gt.0.) then
  3381. etah = .5 * (eta(i,k) + eta(i,k-1))
  3382. if(ncloud.gt.0.) then
  3383. dp = 1000. * del(i,k)
  3384. qlk = dq / (eta(i,k) + etah * (c0 + c1) * dz)
  3385. dellal(i,k) = etah * c1 * dz * qlk * g / dp
  3386. else
  3387. qlk = dq / (eta(i,k) + etah * c0 * dz)
  3388. endif
  3389. qcko(i,k) = qlk + qrch
  3390. pwo(i,k) = etah * c0 * dz * qlk
  3391. pwavo(i) = pwavo(i) + pwo(i,k)
  3392. endif
  3393. endif
  3394. endif
  3395. enddo
  3396. enddo
  3397. !c
  3398. !c exchange ktcon with ktcon1
  3399. !c
  3400. do i = 1, im
  3401. if(cnvflg(i)) then
  3402. kk = ktcon(i)
  3403. ktcon(i) = ktcon1(i)
  3404. ktcon1(i) = kk
  3405. endif
  3406. enddo
  3407. !c
  3408. !c this section is ready for cloud water
  3409. !c
  3410. if(ncloud.gt.0) then
  3411. !c
  3412. !c compute liquid and vapor separation at cloud top
  3413. !c
  3414. do i = 1, im
  3415. if(cnvflg(i)) then
  3416. k = ktcon(i) - 1
  3417. gamma = el2orc * qeso(i,k) / (to(i,k)**2)
  3418. qrch = qeso(i,k) &
  3419. & + gamma * dbyo(i,k) / (hvap * (1. + gamma))
  3420. dq = qcko(i,k) - qrch
  3421. !c
  3422. !c check if there is excess moisture to release latent heat
  3423. !c
  3424. if(dq.gt.0.) then
  3425. qlko_ktcon(i) = dq
  3426. qcko(i,k) = qrch
  3427. endif
  3428. endif
  3429. enddo
  3430. endif
  3431. !c
  3432. !ccccc if(lat.eq.latd.and.lon.eq.lond.and.cnvflg(i)) then
  3433. !ccccc print *, ' aa1(i) before dwndrft =', aa1(i)
  3434. !ccccc endif
  3435. !c
  3436. !c------- downdraft calculations
  3437. !c
  3438. !c--- compute precipitation efficiency in terms of windshear
  3439. !c
  3440. do i = 1, im
  3441. if(cnvflg(i)) then
  3442. vshear(i) = 0.
  3443. endif
  3444. enddo
  3445. do k = 2, km
  3446. do i = 1, im
  3447. if (cnvflg(i)) then
  3448. if(k.gt.kb(i).and.k.le.ktcon(i)) then
  3449. shear= sqrt((uo(i,k)-uo(i,k-1)) ** 2 &
  3450. & + (vo(i,k)-vo(i,k-1)) ** 2)
  3451. vshear(i) = vshear(i) + shear
  3452. endif
  3453. endif
  3454. enddo
  3455. enddo
  3456. do i = 1, im
  3457. if(cnvflg(i)) then
  3458. vshear(i) = 1.e3 * vshear(i) / (zi(i,ktcon(i))-zi(i,kb(i)))
  3459. e1=1.591-.639*vshear(i) &
  3460. & +.0953*(vshear(i)**2)-.00496*(vshear(i)**3)
  3461. edt(i)=1.-e1
  3462. val = .9
  3463. edt(i) = min(edt(i),val)
  3464. val = .0
  3465. edt(i) = max(edt(i),val)
  3466. edto(i)=edt(i)
  3467. edtx(i)=edt(i)
  3468. endif
  3469. enddo
  3470. !c
  3471. !c determine detrainment rate between 1 and kbcon
  3472. !c
  3473. do i = 1, im
  3474. if(cnvflg(i)) then
  3475. sumx(i) = 0.
  3476. endif
  3477. enddo
  3478. do k = 1, km1
  3479. do i = 1, im
  3480. if(cnvflg(i).and.k.ge.1.and.k.lt.kbcon(i)) then
  3481. dz = zi(i,k+1) - zi(i,k)
  3482. sumx(i) = sumx(i) + dz
  3483. endif
  3484. enddo
  3485. enddo
  3486. do i = 1, im
  3487. beta = betas
  3488. if(slimsk(i).eq.1.) beta = betal
  3489. if(cnvflg(i)) then
  3490. dz = (sumx(i)+zi(i,1))/float(kbcon(i))
  3491. tem = 1./float(kbcon(i))
  3492. xlamd(i) = (1.-beta**tem)/dz
  3493. endif
  3494. enddo
  3495. !c
  3496. !c determine downdraft mass flux
  3497. !c
  3498. do k = km1, 1, -1
  3499. do i = 1, im
  3500. if (cnvflg(i) .and. k .le. kmax(i)-1) then
  3501. if(k.lt.jmin(i).and.k.ge.kbcon(i)) then
  3502. dz = zi(i,k+1) - zi(i,k)
  3503. ptem = xlamdd - xlamde
  3504. etad(i,k) = etad(i,k+1) * (1. - ptem * dz)
  3505. else if(k.lt.kbcon(i)) then
  3506. dz = zi(i,k+1) - zi(i,k)
  3507. ptem = xlamd(i) + xlamdd - xlamde
  3508. etad(i,k) = etad(i,k+1) * (1. - ptem * dz)
  3509. endif
  3510. endif
  3511. enddo
  3512. enddo
  3513. !c
  3514. !c--- downdraft moisture properties
  3515. !c
  3516. do i = 1, im
  3517. if(cnvflg(i)) then
  3518. jmn = jmin(i)
  3519. hcdo(i,jmn) = heo(i,jmn)
  3520. qcdo(i,jmn) = qo(i,jmn)
  3521. qrcdo(i,jmn)= qeso(i,jmn)
  3522. ucdo(i,jmn) = uo(i,jmn)
  3523. vcdo(i,jmn) = vo(i,jmn)
  3524. pwevo(i) = 0.
  3525. endif
  3526. enddo
  3527. !cj
  3528. do k = km1, 1, -1
  3529. do i = 1, im
  3530. if (cnvflg(i) .and. k.lt.jmin(i)) then
  3531. dz = zi(i,k+1) - zi(i,k)
  3532. if(k.ge.kbcon(i)) then
  3533. tem = xlamde * dz
  3534. tem1 = 0.5 * xlamdd * dz
  3535. else
  3536. tem = xlamde * dz
  3537. tem1 = 0.5 * (xlamd(i)+xlamdd) * dz
  3538. endif
  3539. factor = 1. + tem - tem1
  3540. ptem = 0.5 * tem - pgcon
  3541. ptem1= 0.5 * tem + pgcon
  3542. hcdo(i,k) = ((1.-tem1)*hcdo(i,k+1)+tem*0.5* &
  3543. & (heo(i,k)+heo(i,k+1)))/factor
  3544. ucdo(i,k) = ((1.-tem1)*ucdo(i,k+1)+ptem*uo(i,k+1) &
  3545. & +ptem1*uo(i,k))/factor
  3546. vcdo(i,k) = ((1.-tem1)*vcdo(i,k+1)+ptem*vo(i,k+1) &
  3547. & +ptem1*vo(i,k))/factor
  3548. dbyo(i,k) = hcdo(i,k) - heso(i,k)
  3549. endif
  3550. enddo
  3551. enddo
  3552. !c
  3553. do k = km1, 1, -1
  3554. do i = 1, im
  3555. if (cnvflg(i).and.k.lt.jmin(i)) then
  3556. gamma = el2orc * qeso(i,k) / (to(i,k)**2)
  3557. qrcdo(i,k) = qeso(i,k)+ &
  3558. & (1./hvap)*(gamma/(1.+gamma))*dbyo(i,k)
  3559. ! detad = etad(i,k+1) - etad(i,k)
  3560. !cj
  3561. dz = zi(i,k+1) - zi(i,k)
  3562. if(k.ge.kbcon(i)) then
  3563. tem = xlamde * dz
  3564. tem1 = 0.5 * xlamdd * dz
  3565. else
  3566. tem = xlamde * dz
  3567. tem1 = 0.5 * (xlamd(i)+xlamdd) * dz
  3568. endif
  3569. factor = 1. + tem - tem1
  3570. qcdo(i,k) = ((1.-tem1)*qcdo(i,k+1)+tem*0.5* &
  3571. & (qo(i,k)+qo(i,k+1)))/factor
  3572. !cj
  3573. ! pwdo(i,k) = etad(i,k+1) * qcdo(i,k+1) -
  3574. ! & etad(i,k) * qrcdo(i,k)
  3575. ! pwdo(i,k) = pwdo(i,k) - detad *
  3576. ! & .5 * (qrcdo(i,k) + qrcdo(i,k+1))
  3577. !cj
  3578. pwdo(i,k) = etad(i,k+1) * (qcdo(i,k) - qrcdo(i,k))
  3579. qcdo(i,k) = qrcdo(i,k)
  3580. pwevo(i) = pwevo(i) + pwdo(i,k)
  3581. endif
  3582. enddo
  3583. enddo
  3584. !c
  3585. !c--- final downdraft strength dependent on precip
  3586. !c--- efficiency (edt), normalized condensate (pwav), and
  3587. !c--- evaporate (pwev)
  3588. !c
  3589. do i = 1, im
  3590. edtmax = edtmaxl
  3591. if(slimsk(i).eq.0.) edtmax = edtmaxs
  3592. if(cnvflg(i)) then
  3593. if(pwevo(i).lt.0.) then
  3594. edto(i) = -edto(i) * pwavo(i) / pwevo(i)
  3595. edto(i) = min(edto(i),edtmax)
  3596. else
  3597. edto(i) = 0.
  3598. endif
  3599. endif
  3600. enddo
  3601. !c
  3602. !c--- downdraft cloudwork functions
  3603. !c
  3604. do k = km1, 1, -1
  3605. do i = 1, im
  3606. if (cnvflg(i) .and. k .lt. jmin(i)) then
  3607. gamma = el2orc * qeso(i,k) / to(i,k)**2
  3608. dhh=hcdo(i,k)
  3609. dt=to(i,k)
  3610. dg=gamma
  3611. dh=heso(i,k)
  3612. dz=-1.*(zo(i,k+1)-zo(i,k))
  3613. aa1(i)=aa1(i)+edto(i)*dz*(g/(cp*dt))*((dhh-dh)/(1.+dg)) &
  3614. & *(1.+delta*cp*dg*dt/hvap)
  3615. val=0.
  3616. aa1(i)=aa1(i)+edto(i)* &
  3617. & dz*g*delta*max(val,(qeso(i,k)-qo(i,k)))
  3618. endif
  3619. enddo
  3620. enddo
  3621. do i = 1, im
  3622. if(cnvflg(i).and.aa1(i).le.0.) then
  3623. cnvflg(i) = .false.
  3624. endif
  3625. enddo
  3626. !!
  3627. totflg = .true.
  3628. do i=1,im
  3629. totflg = totflg .and. (.not. cnvflg(i))
  3630. enddo
  3631. if(totflg) return
  3632. !!
  3633. !c
  3634. !c--- what would the change be, that a cloud with unit mass
  3635. !c--- will do to the environment?
  3636. !c
  3637. do k = 1, km
  3638. do i = 1, im
  3639. if(cnvflg(i) .and. k .le. kmax(i)) then
  3640. dellah(i,k) = 0.
  3641. dellaq(i,k) = 0.
  3642. dellau(i,k) = 0.
  3643. dellav(i,k) = 0.
  3644. endif
  3645. enddo
  3646. enddo
  3647. do i = 1, im
  3648. if(cnvflg(i)) then
  3649. dp = 1000. * del(i,1)
  3650. dellah(i,1) = edto(i) * etad(i,1) * (hcdo(i,1) &
  3651. & - heo(i,1)) * g / dp
  3652. dellaq(i,1) = edto(i) * etad(i,1) * (qcdo(i,1) &
  3653. & - qo(i,1)) * g / dp
  3654. dellau(i,1) = edto(i) * etad(i,1) * (ucdo(i,1) &
  3655. & - uo(i,1)) * g / dp
  3656. dellav(i,1) = edto(i) * etad(i,1) * (vcdo(i,1) &
  3657. & - vo(i,1)) * g / dp
  3658. endif
  3659. enddo
  3660. !c
  3661. !c--- changed due to subsidence and entrainment
  3662. !c
  3663. do k = 2, km1
  3664. do i = 1, im
  3665. if (cnvflg(i).and.k.lt.ktcon(i)) then
  3666. aup = 1.
  3667. if(k.le.kb(i)) aup = 0.
  3668. adw = 1.
  3669. if(k.gt.jmin(i)) adw = 0.
  3670. dp = 1000. * del(i,k)
  3671. dz = zi(i,k) - zi(i,k-1)
  3672. !c
  3673. dv1h = heo(i,k)
  3674. dv2h = .5 * (heo(i,k) + heo(i,k-1))
  3675. dv3h = heo(i,k-1)
  3676. dv1q = qo(i,k)
  3677. dv2q = .5 * (qo(i,k) + qo(i,k-1))
  3678. dv3q = qo(i,k-1)
  3679. dv1u = uo(i,k)
  3680. dv2u = .5 * (uo(i,k) + uo(i,k-1))
  3681. dv3u = uo(i,k-1)
  3682. dv1v = vo(i,k)
  3683. dv2v = .5 * (vo(i,k) + vo(i,k-1))
  3684. dv3v = vo(i,k-1)
  3685. !c
  3686. tem = 0.5 * (xlamue(i,k)+xlamue(i,k-1))
  3687. tem1 = xlamud(i)
  3688. !c
  3689. if(k.le.kbcon(i)) then
  3690. ptem = xlamde
  3691. ptem1 = xlamd(i)+xlamdd
  3692. else
  3693. ptem = xlamde
  3694. ptem1 = xlamdd
  3695. endif
  3696. !cj
  3697. dellah(i,k) = dellah(i,k) + &
  3698. & ((aup*eta(i,k)-adw*edto(i)*etad(i,k))*dv1h &
  3699. & - (aup*eta(i,k-1)-adw*edto(i)*etad(i,k-1))*dv3h &
  3700. & - (aup*tem*eta(i,k-1)+adw*edto(i)*ptem*etad(i,k))*dv2h*dz &
  3701. & + aup*tem1*eta(i,k-1)*.5*(hcko(i,k)+hcko(i,k-1))*dz &
  3702. & + adw*edto(i)*ptem1*etad(i,k)*.5*(hcdo(i,k)+hcdo(i,k-1)) &
  3703. & *dz) *g/dp
  3704. !cj
  3705. dellaq(i,k) = dellaq(i,k) + &
  3706. & ((aup*eta(i,k)-adw*edto(i)*etad(i,k))*dv1q &
  3707. & - (aup*eta(i,k-1)-adw*edto(i)*etad(i,k-1))*dv3q &
  3708. & - (aup*tem*eta(i,k-1)+adw*edto(i)*ptem*etad(i,k))*dv2q*dz &
  3709. & + aup*tem1*eta(i,k-1)*.5*(qcko(i,k)+qcko(i,k-1))*dz &
  3710. & + adw*edto(i)*ptem1*etad(i,k)*.5*(qrcdo(i,k)+qrcdo(i,k-1)) &
  3711. & *dz) *g/dp
  3712. !23456789012345678901234567890123456789012345678901234567890123456789012
  3713. !cj
  3714. dellau(i,k) = dellau(i,k) + &
  3715. & ((aup*eta(i,k)-adw*edto(i)*etad(i,k))*dv1u &
  3716. & - (aup*eta(i,k-1)-adw*edto(i)*etad(i,k-1))*dv3u &
  3717. & - (aup*tem*eta(i,k-1)+adw*edto(i)*ptem*etad(i,k))*dv2u*dz &
  3718. & + aup*tem1*eta(i,k-1)*.5*(ucko(i,k)+ucko(i,k-1))*dz &
  3719. & + adw*edto(i)*ptem1*etad(i,k)*.5*(ucdo(i,k)+ucdo(i,k-1))*dz &
  3720. & - pgcon*(aup*eta(i,k-1)-adw*edto(i)*etad(i,k))*(dv1u-dv3u) &
  3721. & ) *g/dp
  3722. !cj
  3723. dellav(i,k) = dellav(i,k) + &
  3724. & ((aup*eta(i,k)-adw*edto(i)*etad(i,k))*dv1v &
  3725. & - (aup*eta(i,k-1)-adw*edto(i)*etad(i,k-1))*dv3v &
  3726. & - (aup*tem*eta(i,k-1)+adw*edto(i)*ptem*etad(i,k))*dv2v*dz &
  3727. & + aup*tem1*eta(i,k-1)*.5*(vcko(i,k)+vcko(i,k-1))*dz &
  3728. & + adw*edto(i)*ptem1*etad(i,k)*.5*(vcdo(i,k)+vcdo(i,k-1))*dz &
  3729. & - pgcon*(aup*eta(i,k-1)-adw*edto(i)*etad(i,k))*(dv1v-dv3v) &
  3730. & ) *g/dp
  3731. !cj
  3732. endif
  3733. enddo
  3734. enddo
  3735. !c
  3736. !c------- cloud top
  3737. !c
  3738. do i = 1, im
  3739. if(cnvflg(i)) then
  3740. indx = ktcon(i)
  3741. dp = 1000. * del(i,indx)
  3742. dv1h = heo(i,indx-1)
  3743. dellah(i,indx) = eta(i,indx-1) * &
  3744. & (hcko(i,indx-1) - dv1h) * g / dp
  3745. dv1q = qo(i,indx-1)
  3746. dellaq(i,indx) = eta(i,indx-1) * &
  3747. & (qcko(i,indx-1) - dv1q) * g / dp
  3748. dv1u = uo(i,indx-1)
  3749. dellau(i,indx) = eta(i,indx-1) * &
  3750. & (ucko(i,indx-1) - dv1u) * g / dp
  3751. dv1v = vo(i,indx-1)
  3752. dellav(i,indx) = eta(i,indx-1) * &
  3753. & (vcko(i,indx-1) - dv1v) * g / dp
  3754. !c
  3755. !c cloud water
  3756. !c
  3757. dellal(i,indx) = eta(i,indx-1) * &
  3758. & qlko_ktcon(i) * g / dp
  3759. endif
  3760. enddo
  3761. !c
  3762. !c------- final changed variable per unit mass flux
  3763. !c
  3764. do k = 1, km
  3765. do i = 1, im
  3766. if (cnvflg(i).and.k .le. kmax(i)) then
  3767. if(k.gt.ktcon(i)) then
  3768. qo(i,k) = q1(i,k)
  3769. to(i,k) = t1(i,k)
  3770. endif
  3771. if(k.le.ktcon(i)) then
  3772. qo(i,k) = dellaq(i,k) * mbdt + q1(i,k)
  3773. dellat = (dellah(i,k) - hvap * dellaq(i,k)) / cp
  3774. to(i,k) = dellat * mbdt + t1(i,k)
  3775. val = 1.e-10
  3776. qo(i,k) = max(qo(i,k), val )
  3777. endif
  3778. endif
  3779. enddo
  3780. enddo
  3781. !c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  3782. !c
  3783. !c--- the above changed environment is now used to calulate the
  3784. !c--- effect the arbitrary cloud (with unit mass flux)
  3785. !c--- would have on the stability,
  3786. !c--- which then is used to calculate the real mass flux,
  3787. !c--- necessary to keep this change in balance with the large-scale
  3788. !c--- destabilization.
  3789. !c
  3790. !c--- environmental conditions again, first heights
  3791. !c
  3792. do k = 1, km
  3793. do i = 1, im
  3794. if(cnvflg(i) .and. k .le. kmax(i)) then
  3795. qeso(i,k) = 0.01 * fpvs(to(i,k)) ! fpvs is in pa
  3796. qeso(i,k) = eps * qeso(i,k) / (pfld(i,k)+epsm1*qeso(i,k))
  3797. val = 1.e-8
  3798. qeso(i,k) = max(qeso(i,k), val )
  3799. ! tvo(i,k) = to(i,k) + delta * to(i,k) * qo(i,k)
  3800. endif
  3801. enddo
  3802. enddo
  3803. !c
  3804. !c--- moist static energy
  3805. !c
  3806. do k = 1, km1
  3807. do i = 1, im
  3808. if(cnvflg(i) .and. k .le. kmax(i)-1) then
  3809. dz = .5 * (zo(i,k+1) - zo(i,k))
  3810. dp = .5 * (pfld(i,k+1) - pfld(i,k))
  3811. es = 0.01 * fpvs(to(i,k+1)) ! fpvs is in pa
  3812. pprime = pfld(i,k+1) + epsm1 * es
  3813. qs = eps * es / pprime
  3814. dqsdp = - qs / pprime
  3815. desdt = es * (fact1 / to(i,k+1) + fact2 / (to(i,k+1)**2))
  3816. dqsdt = qs * pfld(i,k+1) * desdt / (es * pprime)
  3817. gamma = el2orc * qeso(i,k+1) / (to(i,k+1)**2)
  3818. dt = (g * dz + hvap * dqsdp * dp) / (cp * (1. + gamma))
  3819. dq = dqsdt * dt + dqsdp * dp
  3820. to(i,k) = to(i,k+1) + dt
  3821. qo(i,k) = qo(i,k+1) + dq
  3822. po(i,k) = .5 * (pfld(i,k) + pfld(i,k+1))
  3823. endif
  3824. enddo
  3825. enddo
  3826. do k = 1, km1
  3827. do i = 1, im
  3828. if(cnvflg(i) .and. k .le. kmax(i)-1) then
  3829. qeso(i,k) = 0.01 * fpvs(to(i,k)) ! fpvs is in pa
  3830. qeso(i,k) = eps * qeso(i,k) / (po(i,k) + epsm1 * qeso(i,k))
  3831. val1 = 1.e-8
  3832. qeso(i,k) = max(qeso(i,k), val1)
  3833. val2 = 1.e-10
  3834. qo(i,k) = max(qo(i,k), val2 )
  3835. ! qo(i,k) = min(qo(i,k),qeso(i,k))
  3836. heo(i,k) = .5 * g * (zo(i,k) + zo(i,k+1)) + &
  3837. & cp * to(i,k) + hvap * qo(i,k)
  3838. heso(i,k) = .5 * g * (zo(i,k) + zo(i,k+1)) + &
  3839. & cp * to(i,k) + hvap * qeso(i,k)
  3840. endif
  3841. enddo
  3842. enddo
  3843. do i = 1, im
  3844. if(cnvflg(i)) then
  3845. k = kmax(i)
  3846. heo(i,k) = g * zo(i,k) + cp * to(i,k) + hvap * qo(i,k)
  3847. heso(i,k) = g * zo(i,k) + cp * to(i,k) + hvap * qeso(i,k)
  3848. !c heo(i,k) = min(heo(i,k),heso(i,k))
  3849. endif
  3850. enddo
  3851. !c
  3852. !c**************************** static control
  3853. !c
  3854. !c------- moisture and cloud work functions
  3855. !c
  3856. do i = 1, im
  3857. if(cnvflg(i)) then
  3858. xaa0(i) = 0.
  3859. xpwav(i) = 0.
  3860. endif
  3861. enddo
  3862. !c
  3863. do i = 1, im
  3864. if(cnvflg(i)) then
  3865. indx = kb(i)
  3866. hcko(i,indx) = heo(i,indx)
  3867. qcko(i,indx) = qo(i,indx)
  3868. endif
  3869. enddo
  3870. do k = 2, km1
  3871. do i = 1, im
  3872. if (cnvflg(i)) then
  3873. if(k.gt.kb(i).and.k.le.ktcon(i)) then
  3874. dz = zi(i,k) - zi(i,k-1)
  3875. tem = 0.5 * (xlamue(i,k)+xlamue(i,k-1)) * dz
  3876. tem1 = 0.5 * xlamud(i) * dz
  3877. factor = 1. + tem - tem1
  3878. hcko(i,k) = ((1.-tem1)*hcko(i,k-1)+tem*0.5* &
  3879. & (heo(i,k)+heo(i,k-1)))/factor
  3880. endif
  3881. endif
  3882. enddo
  3883. enddo
  3884. do k = 2, km1
  3885. do i = 1, im
  3886. if (cnvflg(i)) then
  3887. if(k.gt.kb(i).and.k.lt.ktcon(i)) then
  3888. dz = zi(i,k) - zi(i,k-1)
  3889. gamma = el2orc * qeso(i,k) / (to(i,k)**2)
  3890. xdby = hcko(i,k) - heso(i,k)
  3891. xqrch = qeso(i,k) &
  3892. & + gamma * xdby / (hvap * (1. + gamma))
  3893. !cj
  3894. tem = 0.5 * (xlamue(i,k)+xlamue(i,k-1)) * dz
  3895. tem1 = 0.5 * xlamud(i) * dz
  3896. factor = 1. + tem - tem1
  3897. qcko(i,k) = ((1.-tem1)*qcko(i,k-1)+tem*0.5* &
  3898. & (qo(i,k)+qo(i,k-1)))/factor
  3899. !cj
  3900. dq = eta(i,k) * (qcko(i,k) - xqrch)
  3901. !c
  3902. if(k.ge.kbcon(i).and.dq.gt.0.) then
  3903. etah = .5 * (eta(i,k) + eta(i,k-1))
  3904. if(ncloud.gt.0..and.k.gt.jmin(i)) then
  3905. qlk = dq / (eta(i,k) + etah * (c0 + c1) * dz)
  3906. else
  3907. qlk = dq / (eta(i,k) + etah * c0 * dz)
  3908. endif
  3909. if(k.lt.ktcon1(i)) then
  3910. xaa0(i) = xaa0(i) - dz * g * qlk
  3911. endif
  3912. qcko(i,k) = qlk + xqrch
  3913. xpw = etah * c0 * dz * qlk
  3914. xpwav(i) = xpwav(i) + xpw
  3915. endif
  3916. endif
  3917. if(k.ge.kbcon(i).and.k.lt.ktcon1(i)) then
  3918. dz1 = zo(i,k+1) - zo(i,k)
  3919. gamma = el2orc * qeso(i,k) / (to(i,k)**2)
  3920. rfact = 1. + delta * cp * gamma &
  3921. & * to(i,k) / hvap
  3922. xaa0(i) = xaa0(i) &
  3923. & + dz1 * (g / (cp * to(i,k))) &
  3924. & * xdby / (1. + gamma) &
  3925. & * rfact
  3926. val=0.
  3927. xaa0(i)=xaa0(i)+ &
  3928. & dz1 * g * delta * &
  3929. & max(val,(qeso(i,k) - qo(i,k)))
  3930. endif
  3931. endif
  3932. enddo
  3933. enddo
  3934. !c
  3935. !c------- downdraft calculations
  3936. !c
  3937. !c--- downdraft moisture properties
  3938. !c
  3939. do i = 1, im
  3940. if(cnvflg(i)) then
  3941. jmn = jmin(i)
  3942. hcdo(i,jmn) = heo(i,jmn)
  3943. qcdo(i,jmn) = qo(i,jmn)
  3944. qrcd(i,jmn) = qeso(i,jmn)
  3945. xpwev(i) = 0.
  3946. endif
  3947. enddo
  3948. !cj
  3949. do k = km1, 1, -1
  3950. do i = 1, im
  3951. if (cnvflg(i) .and. k.lt.jmin(i)) then
  3952. dz = zi(i,k+1) - zi(i,k)
  3953. if(k.ge.kbcon(i)) then
  3954. tem = xlamde * dz
  3955. tem1 = 0.5 * xlamdd * dz
  3956. else
  3957. tem = xlamde * dz
  3958. tem1 = 0.5 * (xlamd(i)+xlamdd) * dz
  3959. endif
  3960. factor = 1. + tem - tem1
  3961. hcdo(i,k) = ((1.-tem1)*hcdo(i,k+1)+tem*0.5* &
  3962. & (heo(i,k)+heo(i,k+1)))/factor
  3963. endif
  3964. enddo
  3965. enddo
  3966. !cj
  3967. do k = km1, 1, -1
  3968. do i = 1, im
  3969. if (cnvflg(i) .and. k .lt. jmin(i)) then
  3970. dq = qeso(i,k)
  3971. dt = to(i,k)
  3972. gamma = el2orc * dq / dt**2
  3973. dh = hcdo(i,k) - heso(i,k)
  3974. qrcd(i,k)=dq+(1./hvap)*(gamma/(1.+gamma))*dh
  3975. ! detad = etad(i,k+1) - etad(i,k)
  3976. !cj
  3977. dz = zi(i,k+1) - zi(i,k)
  3978. if(k.ge.kbcon(i)) then
  3979. tem = xlamde * dz
  3980. tem1 = 0.5 * xlamdd * dz
  3981. else
  3982. tem = xlamde * dz
  3983. tem1 = 0.5 * (xlamd(i)+xlamdd) * dz
  3984. endif
  3985. factor = 1. + tem - tem1
  3986. qcdo(i,k) = ((1.-tem1)*qcdo(i,k+1)+tem*0.5* &
  3987. & (qo(i,k)+qo(i,k+1)))/factor
  3988. !cj
  3989. ! xpwd = etad(i,k+1) * qcdo(i,k+1) -
  3990. ! & etad(i,k) * qrcd(i,k)
  3991. ! xpwd = xpwd - detad *
  3992. ! & .5 * (qrcd(i,k) + qrcd(i,k+1))
  3993. !cj
  3994. xpwd = etad(i,k+1) * (qcdo(i,k) - qrcd(i,k))
  3995. qcdo(i,k)= qrcd(i,k)
  3996. xpwev(i) = xpwev(i) + xpwd
  3997. endif
  3998. enddo
  3999. enddo
  4000. !c
  4001. do i = 1, im
  4002. edtmax = edtmaxl
  4003. if(slimsk(i).eq.0.) edtmax = edtmaxs
  4004. if(cnvflg(i)) then
  4005. if(xpwev(i).ge.0.) then
  4006. edtx(i) = 0.
  4007. else
  4008. edtx(i) = -edtx(i) * xpwav(i) / xpwev(i)
  4009. edtx(i) = min(edtx(i),edtmax)
  4010. endif
  4011. endif
  4012. enddo
  4013. !c
  4014. !c
  4015. !c--- downdraft cloudwork functions
  4016. !c
  4017. !c
  4018. do k = km1, 1, -1
  4019. do i = 1, im
  4020. if (cnvflg(i) .and. k.lt.jmin(i)) then
  4021. gamma = el2orc * qeso(i,k) / to(i,k)**2
  4022. dhh=hcdo(i,k)
  4023. dt= to(i,k)
  4024. dg= gamma
  4025. dh= heso(i,k)
  4026. dz=-1.*(zo(i,k+1)-zo(i,k))
  4027. xaa0(i)=xaa0(i)+edtx(i)*dz*(g/(cp*dt))*((dhh-dh)/(1.+dg)) &
  4028. & *(1.+delta*cp*dg*dt/hvap)
  4029. val=0.
  4030. xaa0(i)=xaa0(i)+edtx(i)* &
  4031. & dz*g*delta*max(val,(qeso(i,k)-qo(i,k)))
  4032. endif
  4033. enddo
  4034. enddo
  4035. !c
  4036. !c calculate critical cloud work function
  4037. !c
  4038. do i = 1, im
  4039. if(cnvflg(i)) then
  4040. if(pfld(i,ktcon(i)).lt.pcrit(15))then
  4041. acrt(i)=acrit(15)*(975.-pfld(i,ktcon(i))) &
  4042. & /(975.-pcrit(15))
  4043. else if(pfld(i,ktcon(i)).gt.pcrit(1))then
  4044. acrt(i)=acrit(1)
  4045. else
  4046. k = int((850. - pfld(i,ktcon(i)))/50.) + 2
  4047. k = min(k,15)
  4048. k = max(k,2)
  4049. acrt(i)=acrit(k)+(acrit(k-1)-acrit(k))* &
  4050. & (pfld(i,ktcon(i))-pcrit(k))/(pcrit(k-1)-pcrit(k))
  4051. endif
  4052. endif
  4053. enddo
  4054. do i = 1, im
  4055. if(cnvflg(i)) then
  4056. if(slimsk(i).eq.1.) then
  4057. w1 = w1l
  4058. w2 = w2l
  4059. w3 = w3l
  4060. w4 = w4l
  4061. else
  4062. w1 = w1s
  4063. w2 = w2s
  4064. w3 = w3s
  4065. w4 = w4s
  4066. endif
  4067. !c
  4068. !c modify critical cloud workfunction by cloud base vertical velocity
  4069. !c
  4070. if(pdot(i).le.w4) then
  4071. acrtfct(i) = (pdot(i) - w4) / (w3 - w4)
  4072. elseif(pdot(i).ge.-w4) then
  4073. acrtfct(i) = - (pdot(i) + w4) / (w4 - w3)
  4074. else
  4075. acrtfct(i) = 0.
  4076. endif
  4077. val1 = -1.
  4078. acrtfct(i) = max(acrtfct(i),val1)
  4079. val2 = 1.
  4080. acrtfct(i) = min(acrtfct(i),val2)
  4081. acrtfct(i) = 1. - acrtfct(i)
  4082. !c
  4083. !c modify acrtfct(i) by colume mean rh if rhbar(i) is greater than 80 percent
  4084. !c
  4085. !c if(rhbar(i).ge..8) then
  4086. !c acrtfct(i) = acrtfct(i) * (.9 - min(rhbar(i),.9)) * 10.
  4087. !c endif
  4088. !c
  4089. !c modify adjustment time scale by cloud base vertical velocity
  4090. !c
  4091. val1=0.
  4092. dtconv(i) = dt2 + max((1800. - dt2),val1) * &
  4093. & (pdot(i) - w2) / (w1 - w2)
  4094. !c dtconv(i) = max(dtconv(i), dt2)
  4095. !c dtconv(i) = 1800. * (pdot(i) - w2) / (w1 - w2)
  4096. dtconv(i) = max(dtconv(i),dtmin)
  4097. dtconv(i) = min(dtconv(i),dtmax)
  4098. !c
  4099. endif
  4100. enddo
  4101. !c
  4102. !c--- large scale forcing
  4103. !c
  4104. xmbmx1=-1.e20
  4105. do i= 1, im
  4106. if(cnvflg(i)) then
  4107. fld(i)=(aa1(i)-acrt(i)* acrtfct(i))/dtconv(i)
  4108. if(fld(i).le.0.) cnvflg(i) = .false.
  4109. endif
  4110. if(cnvflg(i)) then
  4111. !c xaa0(i) = max(xaa0(i),0.)
  4112. xk(i) = (xaa0(i) - aa1(i)) / mbdt
  4113. if(xk(i).ge.0.) cnvflg(i) = .false.
  4114. endif
  4115. !c
  4116. !c--- kernel, cloud base mass flux
  4117. !c
  4118. if(cnvflg(i)) then
  4119. xmb(i) = -fld(i) / xk(i)
  4120. xmb(i) = min(xmb(i),xmbmax(i))
  4121. xmbmx1=max(xmbmx1,xmb(i))
  4122. endif
  4123. enddo
  4124. ! if(xmbmx1.gt.0.4)print*,'qingfu test xmbmx1=',xmbmx1
  4125. !!
  4126. totflg = .true.
  4127. do i=1,im
  4128. totflg = totflg .and. (.not. cnvflg(i))
  4129. enddo
  4130. if(totflg) return
  4131. !!
  4132. !c
  4133. !c restore to,qo,uo,vo to t1,q1,u1,v1 in case convection stops
  4134. !c
  4135. do k = 1, km
  4136. do i = 1, im
  4137. if (cnvflg(i) .and. k .le. kmax(i)) then
  4138. to(i,k) = t1(i,k)
  4139. qo(i,k) = q1(i,k)
  4140. uo(i,k) = u1(i,k)
  4141. vo(i,k) = v1(i,k)
  4142. qeso(i,k) = 0.01 * fpvs(t1(i,k)) ! fpvs is in pa
  4143. qeso(i,k) = eps * qeso(i,k) / (pfld(i,k) + epsm1*qeso(i,k))
  4144. val = 1.e-8
  4145. qeso(i,k) = max(qeso(i,k), val )
  4146. endif
  4147. enddo
  4148. enddo
  4149. !c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  4150. !c
  4151. !c--- feedback: simply the changes from the cloud with unit mass flux
  4152. !c--- multiplied by the mass flux necessary to keep the
  4153. !c--- equilibrium with the larger-scale.
  4154. !c
  4155. do i = 1, im
  4156. delhbar(i) = 0.
  4157. delqbar(i) = 0.
  4158. deltbar(i) = 0.
  4159. delubar(i) = 0.
  4160. delvbar(i) = 0.
  4161. qcond(i) = 0.
  4162. enddo
  4163. do k = 1, km
  4164. do i = 1, im
  4165. if (cnvflg(i) .and. k .le. kmax(i)) then
  4166. if(k.le.ktcon(i)) then
  4167. dellat = (dellah(i,k) - hvap * dellaq(i,k)) / cp
  4168. t1(i,k) = t1(i,k) + dellat * xmb(i) * dt2
  4169. q1(i,k) = q1(i,k) + dellaq(i,k) * xmb(i) * dt2
  4170. tem = 1./rcs(i)
  4171. u1(i,k) = u1(i,k) + dellau(i,k) * xmb(i) * dt2 * tem
  4172. v1(i,k) = v1(i,k) + dellav(i,k) * xmb(i) * dt2 * tem
  4173. dp = 1000. * del(i,k)
  4174. delhbar(i) = delhbar(i) + dellah(i,k)*xmb(i)*dp/g
  4175. delqbar(i) = delqbar(i) + dellaq(i,k)*xmb(i)*dp/g
  4176. deltbar(i) = deltbar(i) + dellat*xmb(i)*dp/g
  4177. delubar(i) = delubar(i) + dellau(i,k)*xmb(i)*dp/g
  4178. delvbar(i) = delvbar(i) + dellav(i,k)*xmb(i)*dp/g
  4179. endif
  4180. endif
  4181. enddo
  4182. enddo
  4183. do k = 1, km
  4184. do i = 1, im
  4185. if (cnvflg(i) .and. k .le. kmax(i)) then
  4186. if(k.le.ktcon(i)) then
  4187. qeso(i,k) = 0.01 * fpvs(t1(i,k)) ! fpvs is in pa
  4188. qeso(i,k) = eps * qeso(i,k)/(pfld(i,k) + epsm1*qeso(i,k))
  4189. val = 1.e-8
  4190. qeso(i,k) = max(qeso(i,k), val )
  4191. endif
  4192. endif
  4193. enddo
  4194. enddo
  4195. !c
  4196. do i = 1, im
  4197. rntot(i) = 0.
  4198. delqev(i) = 0.
  4199. delq2(i) = 0.
  4200. flg(i) = cnvflg(i)
  4201. enddo
  4202. do k = km, 1, -1
  4203. do i = 1, im
  4204. if (cnvflg(i) .and. k .le. kmax(i)) then
  4205. if(k.lt.ktcon(i)) then
  4206. aup = 1.
  4207. if(k.le.kb(i)) aup = 0.
  4208. adw = 1.
  4209. if(k.ge.jmin(i)) adw = 0.
  4210. rain = aup * pwo(i,k) + adw * edto(i) * pwdo(i,k)
  4211. rntot(i) = rntot(i) + rain * xmb(i) * .001 * dt2
  4212. endif
  4213. endif
  4214. enddo
  4215. enddo
  4216. do k = km, 1, -1
  4217. do i = 1, im
  4218. if (k .le. kmax(i)) then
  4219. deltv(i) = 0.
  4220. delq(i) = 0.
  4221. qevap(i) = 0.
  4222. if(cnvflg(i).and.k.lt.ktcon(i)) then
  4223. aup = 1.
  4224. if(k.le.kb(i)) aup = 0.
  4225. adw = 1.
  4226. if(k.ge.jmin(i)) adw = 0.
  4227. rain = aup * pwo(i,k) + adw * edto(i) * pwdo(i,k)
  4228. rn(i) = rn(i) + rain * xmb(i) * .001 * dt2
  4229. endif
  4230. if(flg(i).and.k.lt.ktcon(i)) then
  4231. evef = edt(i) * evfact
  4232. if(slimsk(i).eq.1.) evef=edt(i) * evfactl
  4233. ! if(slimsk(i).eq.1.) evef=.07
  4234. !c if(slimsk(i).ne.1.) evef = 0.
  4235. qcond(i) = evef * (q1(i,k) - qeso(i,k)) &
  4236. & / (1. + el2orc * qeso(i,k) / t1(i,k)**2)
  4237. dp = 1000. * del(i,k)
  4238. if(rn(i).gt.0..and.qcond(i).lt.0.) then
  4239. qevap(i) = -qcond(i) * (1.-exp(-.32*sqrt(dt2*rn(i))))
  4240. qevap(i) = min(qevap(i), rn(i)*1000.*g/dp)
  4241. delq2(i) = delqev(i) + .001 * qevap(i) * dp / g
  4242. endif
  4243. if(rn(i).gt.0..and.qcond(i).lt.0..and. &
  4244. & delq2(i).gt.rntot(i)) then
  4245. qevap(i) = 1000.* g * (rntot(i) - delqev(i)) / dp
  4246. flg(i) = .false.
  4247. endif
  4248. if(rn(i).gt.0..and.qevap(i).gt.0.) then
  4249. q1(i,k) = q1(i,k) + qevap(i)
  4250. t1(i,k) = t1(i,k) - elocp * qevap(i)
  4251. rn(i) = rn(i) - .001 * qevap(i) * dp / g
  4252. deltv(i) = - elocp*qevap(i)/dt2
  4253. delq(i) = + qevap(i)/dt2
  4254. delqev(i) = delqev(i) + .001*dp*qevap(i)/g
  4255. endif
  4256. dellaq(i,k) = dellaq(i,k) + delq(i) / xmb(i)
  4257. delqbar(i) = delqbar(i) + delq(i)*dp/g
  4258. deltbar(i) = deltbar(i) + deltv(i)*dp/g
  4259. endif
  4260. endif
  4261. enddo
  4262. enddo
  4263. !cj
  4264. ! do i = 1, im
  4265. ! if(me.eq.31.and.cnvflg(i)) then
  4266. ! if(cnvflg(i)) then
  4267. ! print *, ' deep delhbar, delqbar, deltbar = ',
  4268. ! & delhbar(i),hvap*delqbar(i),cp*deltbar(i)
  4269. ! print *, ' deep delubar, delvbar = ',delubar(i),delvbar(i)
  4270. ! print *, ' precip =', hvap*rn(i)*1000./dt2
  4271. ! print*,'pdif= ',pfld(i,kbcon(i))-pfld(i,ktcon(i))
  4272. ! endif
  4273. ! enddo
  4274. !c
  4275. !c precipitation rate converted to actual precip
  4276. !c in unit of m instead of kg
  4277. !c
  4278. do i = 1, im
  4279. if(cnvflg(i)) then
  4280. !c
  4281. !c in the event of upper level rain evaporation and lower level downdraft
  4282. !c moistening, rn can become negative, in this case, we back out of the
  4283. !c heating and the moistening
  4284. !c
  4285. if(rn(i).lt.0..and..not.flg(i)) rn(i) = 0.
  4286. if(rn(i).le.0.) then
  4287. rn(i) = 0.
  4288. else
  4289. ktop(i) = ktcon(i)
  4290. kbot(i) = kbcon(i)
  4291. kcnv(i) = 1
  4292. cldwrk(i) = aa1(i)
  4293. endif
  4294. endif
  4295. enddo
  4296. !c
  4297. !c cloud water
  4298. !c
  4299. if (ncloud.gt.0) then
  4300. !
  4301. val1=1.0
  4302. val2=0.0
  4303. do k = 1, km
  4304. do i = 1, im
  4305. if (cnvflg(i) .and. rn(i).gt.0.) then
  4306. if (k.gt.kb(i).and.k.le.ktcon(i)) then
  4307. tem = dellal(i,k) * xmb(i) * dt2
  4308. tem1 = max(val2, min(val1, (tcr-t1(i,k))*tcrf))
  4309. if (ql(i,k,2) .gt. -999.0) then
  4310. ql(i,k,1) = ql(i,k,1) + tem * tem1 ! ice
  4311. ql(i,k,2) = ql(i,k,2) + tem *(1.0-tem1) ! water
  4312. else
  4313. ql(i,k,1) = ql(i,k,1) + tem
  4314. endif
  4315. endif
  4316. endif
  4317. enddo
  4318. enddo
  4319. !
  4320. endif
  4321. !c
  4322. do k = 1, km
  4323. do i = 1, im
  4324. if(cnvflg(i).and.rn(i).le.0.) then
  4325. if (k .le. kmax(i)) then
  4326. t1(i,k) = to(i,k)
  4327. q1(i,k) = qo(i,k)
  4328. u1(i,k) = uo(i,k)
  4329. v1(i,k) = vo(i,k)
  4330. endif
  4331. endif
  4332. enddo
  4333. enddo
  4334. !
  4335. ! hchuang code change
  4336. !
  4337. ! do k = 1, km
  4338. ! do i = 1, im
  4339. ! if(cnvflg(i).and.rn(i).gt.0.) then
  4340. ! if(k.ge.kb(i) .and. k.lt.ktop(i)) then
  4341. ! ud_mf(i,k) = eta(i,k) * xmb(i) * dt2
  4342. ! endif
  4343. ! endif
  4344. ! enddo
  4345. ! enddo
  4346. ! do i = 1, im
  4347. ! if(cnvflg(i).and.rn(i).gt.0.) then
  4348. ! k = ktop(i)-1
  4349. ! dt_mf(i,k) = ud_mf(i,k)
  4350. ! endif
  4351. ! enddo
  4352. ! do k = 1, km
  4353. ! do i = 1, im
  4354. ! if(cnvflg(i).and.rn(i).gt.0.) then
  4355. ! if(k.ge.1 .and. k.le.jmin(i)) then
  4356. ! dd_mf(i,k) = edto(i) * etad(i,k) * xmb(i) * dt2
  4357. ! endif
  4358. ! endif
  4359. ! enddo
  4360. ! enddo
  4361. !!
  4362. return
  4363. end subroutine sascnvn
  4364. subroutine shalcnv(im,ix,km,jcap,delt,del,prsl,ps,phil,ql, &
  4365. & q1,t1,u1,v1,rcs,rn,kbot,ktop,kcnv,slimsk, &
  4366. & dot,ncloud,hpbl,heat,evap,pgcon)
  4367. !
  4368. use MODULE_GFS_machine , only : kind_phys
  4369. use MODULE_GFS_funcphys , only : fpvs
  4370. use MODULE_GFS_physcons, grav => con_g, cp => con_cp, hvap => con_hvap &
  4371. &, rv => con_rv, fv => con_fvirt, t0c => con_t0c &
  4372. &, rd => con_rd, cvap => con_cvap, cliq => con_cliq &
  4373. &, eps => con_eps, epsm1 => con_epsm1
  4374. implicit none
  4375. !
  4376. integer im, ix, km, jcap, ncloud, &
  4377. & kbot(im), ktop(im), kcnv(im)
  4378. real(kind=kind_phys) delt
  4379. real(kind=kind_phys) ps(im), del(ix,km), prsl(ix,km), &
  4380. & ql(ix,km,2),q1(ix,km), t1(ix,km), &
  4381. & u1(ix,km), v1(ix,km), rcs(im), &
  4382. & rn(im), slimsk(im), &
  4383. & dot(ix,km), phil(ix,km), hpbl(im), &
  4384. & heat(im), evap(im)
  4385. ! &, ud_mf(im,km),dt_mf(im,km)
  4386. real ud_mf(im,km),dt_mf(im,km)
  4387. !
  4388. integer i,j,indx, jmn, k, kk, latd, lond, km1
  4389. integer kpbl(im)
  4390. !
  4391. real(kind=kind_phys) c0, cpoel, dellat, delta, &
  4392. & desdt, deta, detad, dg, &
  4393. & dh, dhh, dlnsig, dp, &
  4394. & dq, dqsdp, dqsdt, dt, &
  4395. & dt2, dtmax, dtmin, dv1h, &
  4396. & dv1q, dv2h, dv2q, dv1u, &
  4397. & dv1v, dv2u, dv2v, dv3q, &
  4398. & dv3h, dv3u, dv3v, clam, &
  4399. & dz, dz1, e1, &
  4400. & el2orc, elocp, aafac, &
  4401. & es, etah, h1, dthk, &
  4402. & evef, evfact, evfactl, fact1, &
  4403. & fact2, factor, fjcap, &
  4404. & g, gamma, pprime, betaw, &
  4405. & qlk, qrch, qs, c1, &
  4406. & rain, rfact, shear, tem1, &
  4407. & tem2, terr, val, val1, &
  4408. & val2, w1, w1l, w1s, &
  4409. & w2, w2l, w2s, w3, &
  4410. & w3l, w3s, w4, w4l, &
  4411. & w4s, tem, ptem, ptem1, &
  4412. & pgcon
  4413. !
  4414. integer kb(im), kbcon(im), kbcon1(im), &
  4415. & ktcon(im), ktcon1(im), &
  4416. & kbm(im), kmax(im)
  4417. !
  4418. real(kind=kind_phys) aa1(im), &
  4419. & delhbar(im), delq(im), delq2(im), &
  4420. & delqbar(im), delqev(im), deltbar(im), &
  4421. & deltv(im), edt(im), &
  4422. & wstar(im), sflx(im), &
  4423. & pdot(im), po(im,km), &
  4424. & qcond(im), qevap(im), hmax(im), &
  4425. & rntot(im), vshear(im), &
  4426. & xlamud(im), xmb(im), xmbmax(im), &
  4427. & delubar(im), delvbar(im)
  4428. !c
  4429. real(kind=kind_phys) cincr, cincrmax, cincrmin
  4430. !cc
  4431. !c physical parameters
  4432. parameter(g=grav)
  4433. parameter(cpoel=cp/hvap,elocp=hvap/cp, &
  4434. & el2orc=hvap*hvap/(rv*cp))
  4435. parameter(terr=0.,c0=.002,c1=5.e-4,delta=fv)
  4436. parameter(fact1=(cvap-cliq)/rv,fact2=hvap/rv-fact1*t0c)
  4437. parameter(cincrmax=180.,cincrmin=120.,dthk=25.)
  4438. parameter(h1=0.33333333)
  4439. !c local variables and arrays
  4440. real(kind=kind_phys) pfld(im,km), to(im,km), qo(im,km), &
  4441. & uo(im,km), vo(im,km), qeso(im,km)
  4442. !c cloud water
  4443. real(kind=kind_phys) qlko_ktcon(im), dellal(im,km), &
  4444. & dbyo(im,km), zo(im,km), xlamue(im,km), &
  4445. & heo(im,km), heso(im,km), &
  4446. & dellah(im,km), dellaq(im,km), &
  4447. & dellau(im,km), dellav(im,km), hcko(im,km), &
  4448. & ucko(im,km), vcko(im,km), qcko(im,km), &
  4449. & eta(im,km), zi(im,km), pwo(im,km), &
  4450. & tx1(im)
  4451. !
  4452. logical totflg, cnvflg(im), flg(im)
  4453. !
  4454. real(kind=kind_phys) tf, tcr, tcrf
  4455. parameter (tf=233.16, tcr=263.16, tcrf=1.0/(tcr-tf))
  4456. !
  4457. !c-----------------------------------------------------------------------
  4458. !
  4459. km1 = km - 1
  4460. !c
  4461. !c compute surface buoyancy flux
  4462. !c
  4463. do i=1,im
  4464. sflx(i) = heat(i)+fv*t1(i,1)*evap(i)
  4465. enddo
  4466. !c
  4467. !c initialize arrays
  4468. !c
  4469. do i=1,im
  4470. cnvflg(i) = .true.
  4471. if(kcnv(i).eq.1) cnvflg(i) = .false.
  4472. if(sflx(i).le.0.) cnvflg(i) = .false.
  4473. if(cnvflg(i)) then
  4474. kbot(i)=km+1
  4475. ktop(i)=0
  4476. endif
  4477. rn(i)=0.
  4478. kbcon(i)=km
  4479. ktcon(i)=1
  4480. kb(i)=km
  4481. pdot(i) = 0.
  4482. qlko_ktcon(i) = 0.
  4483. edt(i) = 0.
  4484. aa1(i) = 0.
  4485. vshear(i) = 0.
  4486. enddo
  4487. ! hchuang code change
  4488. do k = 1, km
  4489. do i = 1, im
  4490. ud_mf(i,k) = 0.
  4491. dt_mf(i,k) = 0.
  4492. enddo
  4493. enddo
  4494. !!
  4495. totflg = .true.
  4496. do i=1,im
  4497. totflg = totflg .and. (.not. cnvflg(i))
  4498. enddo
  4499. if(totflg) return
  4500. !!
  4501. !c
  4502. dt2 = delt
  4503. val = 1200.
  4504. dtmin = max(dt2, val )
  4505. val = 3600.
  4506. dtmax = max(dt2, val )
  4507. !c model tunable parameters are all here
  4508. clam = .3
  4509. aafac = .1
  4510. betaw = .03
  4511. !c evef = 0.07
  4512. evfact = 0.3
  4513. evfactl = 0.3
  4514. !
  4515. fjcap = (float(jcap) / 126.) ** 2
  4516. val = 1.
  4517. fjcap = max(fjcap,val)
  4518. w1l = -8.e-3
  4519. w2l = -4.e-2
  4520. w3l = -5.e-3
  4521. w4l = -5.e-4
  4522. w1s = -2.e-4
  4523. w2s = -2.e-3
  4524. w3s = -1.e-3
  4525. w4s = -2.e-5
  4526. !c
  4527. !c define top layer for search of the downdraft originating layer
  4528. !c and the maximum thetae for updraft
  4529. !c
  4530. do i=1,im
  4531. kbm(i) = km
  4532. kmax(i) = km
  4533. tx1(i) = 1.0 / ps(i)
  4534. enddo
  4535. !
  4536. do k = 1, km
  4537. do i=1,im
  4538. if (prsl(i,k)*tx1(i) .gt. 0.70) kbm(i) = k + 1
  4539. if (prsl(i,k)*tx1(i) .gt. 0.60) kmax(i) = k + 1
  4540. enddo
  4541. enddo
  4542. do i=1,im
  4543. kbm(i) = min(kbm(i),kmax(i))
  4544. enddo
  4545. !c
  4546. !!c hydrostatic height assume zero terr and compute
  4547. !c updraft entrainment rate as an inverse function of height
  4548. !c
  4549. do k = 1, km
  4550. do i=1,im
  4551. zo(i,k) = phil(i,k) / g
  4552. enddo
  4553. enddo
  4554. do k = 1, km1
  4555. do i=1,im
  4556. zi(i,k) = 0.5*(zo(i,k)+zo(i,k+1))
  4557. xlamue(i,k) = clam / zi(i,k)
  4558. enddo
  4559. enddo
  4560. do i=1,im
  4561. xlamue(i,km) = xlamue(i,km1)
  4562. enddo
  4563. !c
  4564. !c pbl height
  4565. !c
  4566. do i=1,im
  4567. flg(i) = cnvflg(i)
  4568. kpbl(i)= 1
  4569. enddo
  4570. do k = 2, km1
  4571. do i=1,im
  4572. if (flg(i).and.zo(i,k).le.hpbl(i)) then
  4573. kpbl(i) = k
  4574. else
  4575. flg(i) = .false.
  4576. endif
  4577. enddo
  4578. enddo
  4579. do i=1,im
  4580. kpbl(i)= min(kpbl(i),kbm(i))
  4581. enddo
  4582. !c
  4583. !c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  4584. !c convert surface pressure to mb from cb
  4585. !c
  4586. do k = 1, km
  4587. do i = 1, im
  4588. if (cnvflg(i) .and. k .le. kmax(i)) then
  4589. pfld(i,k) = prsl(i,k) * 10.0
  4590. eta(i,k) = 1.
  4591. hcko(i,k) = 0.
  4592. qcko(i,k) = 0.
  4593. ucko(i,k) = 0.
  4594. vcko(i,k) = 0.
  4595. dbyo(i,k) = 0.
  4596. pwo(i,k) = 0.
  4597. dellal(i,k) = 0.
  4598. to(i,k) = t1(i,k)
  4599. qo(i,k) = q1(i,k)
  4600. uo(i,k) = u1(i,k) * rcs(i)
  4601. vo(i,k) = v1(i,k) * rcs(i)
  4602. endif
  4603. enddo
  4604. enddo
  4605. !c
  4606. !c column variables
  4607. !c p is pressure of the layer (mb)
  4608. !c t is temperature at t-dt (k)..tn
  4609. !c q is mixing ratio at t-dt (kg/kg)..qn
  4610. !c to is temperature at t+dt (k)... this is after advection and turbulan
  4611. !c qo is mixing ratio at t+dt (kg/kg)..q1
  4612. !c
  4613. do k = 1, km
  4614. do i=1,im
  4615. if (cnvflg(i) .and. k .le. kmax(i)) then
  4616. qeso(i,k) = 0.01 * fpvs(to(i,k)) ! fpvs is in pa
  4617. qeso(i,k) = eps * qeso(i,k) / (pfld(i,k) + epsm1*qeso(i,k))
  4618. val1 = 1.e-8
  4619. qeso(i,k) = max(qeso(i,k), val1)
  4620. val2 = 1.e-10
  4621. qo(i,k) = max(qo(i,k), val2 )
  4622. ! qo(i,k) = min(qo(i,k),qeso(i,k))
  4623. ! tvo(i,k) = to(i,k) + delta * to(i,k) * qo(i,k)
  4624. endif
  4625. enddo
  4626. enddo
  4627. !c
  4628. !c compute moist static energy
  4629. !c
  4630. do k = 1, km
  4631. do i=1,im
  4632. if (cnvflg(i) .and. k .le. kmax(i)) then
  4633. ! tem = g * zo(i,k) + cp * to(i,k)
  4634. tem = phil(i,k) + cp * to(i,k)
  4635. heo(i,k) = tem + hvap * qo(i,k)
  4636. heso(i,k) = tem + hvap * qeso(i,k)
  4637. !c heo(i,k) = min(heo(i,k),heso(i,k))
  4638. endif
  4639. enddo
  4640. enddo
  4641. !c
  4642. !c determine level with largest moist static energy within pbl
  4643. !c this is the level where updraft starts
  4644. !c
  4645. do i=1,im
  4646. if (cnvflg(i)) then
  4647. hmax(i) = heo(i,1)
  4648. kb(i) = 1
  4649. endif
  4650. enddo
  4651. do k = 2, km
  4652. do i=1,im
  4653. if (cnvflg(i).and.k.le.kpbl(i)) then
  4654. if(heo(i,k).gt.hmax(i)) then
  4655. kb(i) = k
  4656. hmax(i) = heo(i,k)
  4657. endif
  4658. endif
  4659. enddo
  4660. enddo
  4661. !c
  4662. do k = 1, km1
  4663. do i=1,im
  4664. if (cnvflg(i) .and. k .le. kmax(i)-1) then
  4665. dz = .5 * (zo(i,k+1) - zo(i,k))
  4666. dp = .5 * (pfld(i,k+1) - pfld(i,k))
  4667. es = 0.01 * fpvs(to(i,k+1)) ! fpvs is in pa
  4668. pprime = pfld(i,k+1) + epsm1 * es
  4669. qs = eps * es / pprime
  4670. dqsdp = - qs / pprime
  4671. desdt = es * (fact1 / to(i,k+1) + fact2 / (to(i,k+1)**2))
  4672. dqsdt = qs * pfld(i,k+1) * desdt / (es * pprime)
  4673. gamma = el2orc * qeso(i,k+1) / (to(i,k+1)**2)
  4674. dt = (g * dz + hvap * dqsdp * dp) / (cp * (1. + gamma))
  4675. dq = dqsdt * dt + dqsdp * dp
  4676. to(i,k) = to(i,k+1) + dt
  4677. qo(i,k) = qo(i,k+1) + dq
  4678. po(i,k) = .5 * (pfld(i,k) + pfld(i,k+1))
  4679. endif
  4680. enddo
  4681. enddo
  4682. !
  4683. do k = 1, km1
  4684. do i=1,im
  4685. if (cnvflg(i) .and. k .le. kmax(i)-1) then
  4686. qeso(i,k) = 0.01 * fpvs(to(i,k)) ! fpvs is in pa
  4687. qeso(i,k) = eps * qeso(i,k) / (po(i,k) + epsm1*qeso(i,k))
  4688. val1 = 1.e-8
  4689. qeso(i,k) = max(qeso(i,k), val1)
  4690. val2 = 1.e-10
  4691. qo(i,k) = max(qo(i,k), val2 )
  4692. ! qo(i,k) = min(qo(i,k),qeso(i,k))
  4693. heo(i,k) = .5 * g * (zo(i,k) + zo(i,k+1)) + &
  4694. & cp * to(i,k) + hvap * qo(i,k)
  4695. heso(i,k) = .5 * g * (zo(i,k) + zo(i,k+1)) + &
  4696. & cp * to(i,k) + hvap * qeso(i,k)
  4697. uo(i,k) = .5 * (uo(i,k) + uo(i,k+1))
  4698. vo(i,k) = .5 * (vo(i,k) + vo(i,k+1))
  4699. endif
  4700. enddo
  4701. enddo
  4702. !c
  4703. !c look for the level of free convection as cloud base
  4704. !c
  4705. do i=1,im
  4706. flg(i) = cnvflg(i)
  4707. if(flg(i)) kbcon(i) = kmax(i)
  4708. enddo
  4709. do k = 2, km1
  4710. do i=1,im
  4711. if (flg(i).and.k.lt.kbm(i)) then
  4712. if(k.gt.kb(i).and.heo(i,kb(i)).gt.heso(i,k)) then
  4713. kbcon(i) = k
  4714. flg(i) = .false.
  4715. endif
  4716. endif
  4717. enddo
  4718. enddo
  4719. !c
  4720. do i=1,im
  4721. if(cnvflg(i)) then
  4722. if(kbcon(i).eq.kmax(i)) cnvflg(i) = .false.
  4723. endif
  4724. enddo
  4725. !!
  4726. totflg = .true.
  4727. do i=1,im
  4728. totflg = totflg .and. (.not. cnvflg(i))
  4729. enddo
  4730. if(totflg) return
  4731. !!
  4732. !c
  4733. !c determine critical convective inhibition
  4734. !c as a function of vertical velocity at cloud base.
  4735. !c
  4736. do i=1,im
  4737. if(cnvflg(i)) then
  4738. pdot(i) = 10.* dot(i,kbcon(i))
  4739. endif
  4740. enddo
  4741. do i=1,im
  4742. if(cnvflg(i)) then
  4743. if(slimsk(i).eq.1.) then
  4744. w1 = w1l
  4745. w2 = w2l
  4746. w3 = w3l
  4747. w4 = w4l
  4748. else
  4749. w1 = w1s
  4750. w2 = w2s
  4751. w3 = w3s
  4752. w4 = w4s
  4753. endif
  4754. if(pdot(i).le.w4) then
  4755. ptem = (pdot(i) - w4) / (w3 - w4)
  4756. elseif(pdot(i).ge.-w4) then
  4757. ptem = - (pdot(i) + w4) / (w4 - w3)
  4758. else
  4759. ptem = 0.
  4760. endif
  4761. val1 = -1.
  4762. ptem = max(ptem,val1)
  4763. val2 = 1.
  4764. ptem = min(ptem,val2)
  4765. ptem = 1. - ptem
  4766. ptem1= .5*(cincrmax-cincrmin)
  4767. cincr = cincrmax - ptem * ptem1
  4768. tem1 = pfld(i,kb(i)) - pfld(i,kbcon(i))
  4769. if(tem1.gt.cincr) then
  4770. cnvflg(i) = .false.
  4771. endif
  4772. endif
  4773. enddo
  4774. !!
  4775. totflg = .true.
  4776. do i=1,im
  4777. totflg = totflg .and. (.not. cnvflg(i))
  4778. enddo
  4779. if(totflg) return
  4780. !!
  4781. !c
  4782. !c assume the detrainment rate for the updrafts to be same as
  4783. !c the entrainment rate at cloud base
  4784. !c
  4785. do i = 1, im
  4786. if(cnvflg(i)) then
  4787. xlamud(i) = xlamue(i,kbcon(i))
  4788. endif
  4789. enddo
  4790. !c
  4791. !c determine updraft mass flux for the subcloud layers
  4792. !c
  4793. do k = km1, 1, -1
  4794. do i = 1, im
  4795. if (cnvflg(i)) then
  4796. if(k.lt.kbcon(i).and.k.ge.kb(i)) then
  4797. dz = zi(i,k+1) - zi(i,k)
  4798. ptem = 0.5*(xlamue(i,k)+xlamue(i,k+1))-xlamud(i)
  4799. eta(i,k) = eta(i,k+1) / (1. + ptem * dz)
  4800. endif
  4801. endif
  4802. enddo
  4803. enddo
  4804. !c
  4805. !c compute mass flux above cloud base
  4806. !c
  4807. do k = 2, km1
  4808. do i = 1, im
  4809. if(cnvflg(i))then
  4810. if(k.gt.kbcon(i).and.k.lt.kmax(i)) then
  4811. dz = zi(i,k) - zi(i,k-1)
  4812. ptem = 0.5*(xlamue(i,k)+xlamue(i,k-1))-xlamud(i)
  4813. eta(i,k) = eta(i,k-1) * (1 + ptem * dz)
  4814. endif
  4815. endif
  4816. enddo
  4817. enddo
  4818. !c
  4819. !c compute updraft cloud property
  4820. !c
  4821. do i = 1, im
  4822. if(cnvflg(i)) then
  4823. indx = kb(i)
  4824. hcko(i,indx) = heo(i,indx)
  4825. ucko(i,indx) = uo(i,indx)
  4826. vcko(i,indx) = vo(i,indx)
  4827. endif
  4828. enddo
  4829. !c
  4830. do k = 2, km1
  4831. do i = 1, im
  4832. if (cnvflg(i)) then
  4833. if(k.gt.kb(i).and.k.lt.kmax(i)) then
  4834. dz = zi(i,k) - zi(i,k-1)
  4835. tem = 0.5 * (xlamue(i,k)+xlamue(i,k-1)) * dz
  4836. tem1 = 0.5 * xlamud(i) * dz
  4837. factor = 1. + tem - tem1
  4838. ptem = 0.5 * tem + pgcon
  4839. ptem1= 0.5 * tem - pgcon
  4840. hcko(i,k) = ((1.-tem1)*hcko(i,k-1)+tem*0.5* &
  4841. & (heo(i,k)+heo(i,k-1)))/factor
  4842. ucko(i,k) = ((1.-tem1)*ucko(i,k-1)+ptem*uo(i,k) &
  4843. & +ptem1*uo(i,k-1))/factor
  4844. vcko(i,k) = ((1.-tem1)*vcko(i,k-1)+ptem*vo(i,k) &
  4845. & +ptem1*vo(i,k-1))/factor
  4846. dbyo(i,k) = hcko(i,k) - heso(i,k)
  4847. endif
  4848. endif
  4849. enddo
  4850. enddo
  4851. !c
  4852. !c taking account into convection inhibition due to existence of
  4853. !c dry layers below cloud base
  4854. !c
  4855. do i=1,im
  4856. flg(i) = cnvflg(i)
  4857. kbcon1(i) = kmax(i)
  4858. enddo
  4859. do k = 2, km1
  4860. do i=1,im
  4861. if (flg(i).and.k.lt.kbm(i)) then
  4862. if(k.ge.kbcon(i).and.dbyo(i,k).gt.0.) then
  4863. kbcon1(i) = k
  4864. flg(i) = .false.
  4865. endif
  4866. endif
  4867. enddo
  4868. enddo
  4869. do i=1,im
  4870. if(cnvflg(i)) then
  4871. if(kbcon1(i).eq.kmax(i)) cnvflg(i) = .false.
  4872. endif
  4873. enddo
  4874. do i=1,im
  4875. if(cnvflg(i)) then
  4876. tem = pfld(i,kbcon(i)) - pfld(i,kbcon1(i))
  4877. if(tem.gt.dthk) then
  4878. cnvflg(i) = .false.
  4879. endif
  4880. endif
  4881. enddo
  4882. !!
  4883. totflg = .true.
  4884. do i = 1, im
  4885. totflg = totflg .and. (.not. cnvflg(i))
  4886. enddo
  4887. if(totflg) return
  4888. !!
  4889. !c
  4890. !c determine first guess cloud top as the level of zero buoyancy
  4891. !c limited to the level of sigma=0.7
  4892. !c
  4893. do i = 1, im
  4894. flg(i) = cnvflg(i)
  4895. if(flg(i)) ktcon(i) = kbm(i)
  4896. enddo
  4897. do k = 2, km1
  4898. do i=1,im
  4899. if (flg(i).and.k .lt. kbm(i)) then
  4900. if(k.gt.kbcon1(i).and.dbyo(i,k).lt.0.) then
  4901. ktcon(i) = k
  4902. flg(i) = .false.
  4903. endif
  4904. endif
  4905. enddo
  4906. enddo
  4907. !c
  4908. !c turn off shallow convection if cloud top is less than pbl top
  4909. !c
  4910. ! do i=1,im
  4911. ! if(cnvflg(i)) then
  4912. ! kk = kpbl(i)+1
  4913. ! if(ktcon(i).le.kk) cnvflg(i) = .false.
  4914. ! endif
  4915. ! enddo
  4916. !!
  4917. ! totflg = .true.
  4918. ! do i = 1, im
  4919. ! totflg = totflg .and. (.not. cnvflg(i))
  4920. ! enddo
  4921. ! if(totflg) return
  4922. !!
  4923. !c
  4924. !c specify upper limit of mass flux at cloud base
  4925. !c
  4926. do i = 1, im
  4927. if(cnvflg(i)) then
  4928. ! xmbmax(i) = .1
  4929. !
  4930. k = kbcon(i)
  4931. dp = 1000. * del(i,k)
  4932. xmbmax(i) = dp / (g * dt2)
  4933. !
  4934. ! tem = dp / (g * dt2)
  4935. ! xmbmax(i) = min(tem, xmbmax(i))
  4936. endif
  4937. enddo
  4938. !c
  4939. !c compute cloud moisture property and precipitation
  4940. !c
  4941. do i = 1, im
  4942. if (cnvflg(i)) then
  4943. aa1(i) = 0.
  4944. qcko(i,kb(i)) = qo(i,kb(i))
  4945. endif
  4946. enddo
  4947. do k = 2, km1
  4948. do i = 1, im
  4949. if (cnvflg(i)) then
  4950. if(k.gt.kb(i).and.k.lt.ktcon(i)) then
  4951. dz = zi(i,k) - zi(i,k-1)
  4952. gamma = el2orc * qeso(i,k) / (to(i,k)**2)
  4953. qrch = qeso(i,k) &
  4954. & + gamma * dbyo(i,k) / (hvap * (1. + gamma))
  4955. !cj
  4956. tem = 0.5 * (xlamue(i,k)+xlamue(i,k-1)) * dz
  4957. tem1 = 0.5 * xlamud(i) * dz
  4958. factor = 1. + tem - tem1
  4959. qcko(i,k) = ((1.-tem1)*qcko(i,k-1)+tem*0.5* &
  4960. & (qo(i,k)+qo(i,k-1)))/factor
  4961. !cj
  4962. dq = eta(i,k) * (qcko(i,k) - qrch)
  4963. !c
  4964. ! rhbar(i) = rhbar(i) + qo(i,k) / qeso(i,k)
  4965. !c
  4966. !c below lfc check if there is excess moisture to release latent heat
  4967. !c
  4968. if(k.ge.kbcon(i).and.dq.gt.0.) then
  4969. etah = .5 * (eta(i,k) + eta(i,k-1))
  4970. if(ncloud.gt.0.) then
  4971. dp = 1000. * del(i,k)
  4972. qlk = dq / (eta(i,k) + etah * (c0 + c1) * dz)
  4973. dellal(i,k) = etah * c1 * dz * qlk * g / dp
  4974. else
  4975. qlk = dq / (eta(i,k) + etah * c0 * dz)
  4976. endif
  4977. aa1(i) = aa1(i) - dz * g * qlk
  4978. qcko(i,k)= qlk + qrch
  4979. pwo(i,k) = etah * c0 * dz * qlk
  4980. endif
  4981. endif
  4982. endif
  4983. enddo
  4984. enddo
  4985. !c
  4986. !c calculate cloud work function
  4987. !c
  4988. do k = 2, km1
  4989. do i = 1, im
  4990. if (cnvflg(i)) then
  4991. if(k.ge.kbcon(i).and.k.lt.ktcon(i)) then
  4992. dz1 = zo(i,k+1) - zo(i,k)
  4993. gamma = el2orc * qeso(i,k) / (to(i,k)**2)
  4994. rfact = 1. + delta * cp * gamma &
  4995. & * to(i,k) / hvap
  4996. aa1(i) = aa1(i) + &
  4997. & dz1 * (g / (cp * to(i,k))) &
  4998. & * dbyo(i,k) / (1. + gamma) &
  4999. & * rfact
  5000. val = 0.
  5001. aa1(i)=aa1(i)+ &
  5002. & dz1 * g * delta * &
  5003. & max(val,(qeso(i,k) - qo(i,k)))
  5004. endif
  5005. endif
  5006. enddo
  5007. enddo
  5008. do i = 1, im
  5009. if(cnvflg(i).and.aa1(i).le.0.) cnvflg(i) = .false.
  5010. enddo
  5011. !!
  5012. totflg = .true.
  5013. do i=1,im
  5014. totflg = totflg .and. (.not. cnvflg(i))
  5015. enddo
  5016. if(totflg) return
  5017. !!
  5018. !c
  5019. !c estimate the onvective overshooting as the level
  5020. !c where the [aafac * cloud work function] becomes zero,
  5021. !c which is the final cloud top
  5022. !c limited to the level of sigma=0.7
  5023. !c
  5024. do i = 1, im
  5025. if (cnvflg(i)) then
  5026. aa1(i) = aafac * aa1(i)
  5027. endif
  5028. enddo
  5029. !c
  5030. do i = 1, im
  5031. flg(i) = cnvflg(i)
  5032. ktcon1(i) = kbm(i)
  5033. enddo
  5034. do k = 2, km1
  5035. do i = 1, im
  5036. if (flg(i)) then
  5037. if(k.ge.ktcon(i).and.k.lt.kbm(i)) then
  5038. dz1 = zo(i,k+1) - zo(i,k)
  5039. gamma = el2orc * qeso(i,k) / (to(i,k)**2)
  5040. rfact = 1. + delta * cp * gamma &
  5041. & * to(i,k) / hvap
  5042. aa1(i) = aa1(i) + &
  5043. & dz1 * (g / (cp * to(i,k))) &
  5044. & * dbyo(i,k) / (1. + gamma) &
  5045. & * rfact
  5046. if(aa1(i).lt.0.) then
  5047. ktcon1(i) = k
  5048. flg(i) = .false.
  5049. endif
  5050. endif
  5051. endif
  5052. enddo
  5053. enddo
  5054. !c
  5055. !c compute cloud moisture property, detraining cloud water
  5056. !c and precipitation in overshooting layers
  5057. !c
  5058. do k = 2, km1
  5059. do i = 1, im
  5060. if (cnvflg(i)) then
  5061. if(k.ge.ktcon(i).and.k.lt.ktcon1(i)) then
  5062. dz = zi(i,k) - zi(i,k-1)
  5063. gamma = el2orc * qeso(i,k) / (to(i,k)**2)
  5064. qrch = qeso(i,k) &
  5065. & + gamma * dbyo(i,k) / (hvap * (1. + gamma))
  5066. !cj
  5067. tem = 0.5 * (xlamue(i,k)+xlamue(i,k-1)) * dz
  5068. tem1 = 0.5 * xlamud(i) * dz
  5069. factor = 1. + tem - tem1
  5070. qcko(i,k) = ((1.-tem1)*qcko(i,k-1)+tem*0.5* &
  5071. & (qo(i,k)+qo(i,k-1)))/factor
  5072. !cj
  5073. dq = eta(i,k) * (qcko(i,k) - qrch)
  5074. !c
  5075. !c check if there is excess moisture to release latent heat
  5076. !c
  5077. if(dq.gt.0.) then
  5078. etah = .5 * (eta(i,k) + eta(i,k-1))
  5079. if(ncloud.gt.0.) then
  5080. dp = 1000. * del(i,k)
  5081. qlk = dq / (eta(i,k) + etah * (c0 + c1) * dz)
  5082. dellal(i,k) = etah * c1 * dz * qlk * g / dp
  5083. else
  5084. qlk = dq / (eta(i,k) + etah * c0 * dz)
  5085. endif
  5086. qcko(i,k) = qlk + qrch
  5087. pwo(i,k) = etah * c0 * dz * qlk
  5088. endif
  5089. endif
  5090. endif
  5091. enddo
  5092. enddo
  5093. !c
  5094. !c exchange ktcon with ktcon1
  5095. !c
  5096. do i = 1, im
  5097. if(cnvflg(i)) then
  5098. kk = ktcon(i)
  5099. ktcon(i) = ktcon1(i)
  5100. ktcon1(i) = kk
  5101. endif
  5102. enddo
  5103. !c
  5104. !c this section is ready for cloud water
  5105. !c
  5106. if(ncloud.gt.0) then
  5107. !c
  5108. !c compute liquid and vapor separation at cloud top
  5109. !c
  5110. do i = 1, im
  5111. if(cnvflg(i)) then
  5112. k = ktcon(i) - 1
  5113. gamma = el2orc * qeso(i,k) / (to(i,k)**2)
  5114. qrch = qeso(i,k) &
  5115. & + gamma * dbyo(i,k) / (hvap * (1. + gamma))
  5116. dq = qcko(i,k) - qrch
  5117. !c
  5118. !c check if there is excess moisture to release latent heat
  5119. !c
  5120. if(dq.gt.0.) then
  5121. qlko_ktcon(i) = dq
  5122. qcko(i,k) = qrch
  5123. endif
  5124. endif
  5125. enddo
  5126. endif
  5127. !!c
  5128. !c--- compute precipitation efficiency in terms of windshear
  5129. !c
  5130. do i = 1, im
  5131. if(cnvflg(i)) then
  5132. vshear(i) = 0.
  5133. endif
  5134. enddo
  5135. do k = 2, km
  5136. do i = 1, im
  5137. if (cnvflg(i)) then
  5138. if(k.gt.kb(i).and.k.le.ktcon(i)) then
  5139. shear= sqrt((uo(i,k)-uo(i,k-1)) ** 2 &
  5140. & + (vo(i,k)-vo(i,k-1)) ** 2)
  5141. vshear(i) = vshear(i) + shear
  5142. endif
  5143. endif
  5144. enddo
  5145. enddo
  5146. do i = 1, im
  5147. if(cnvflg(i)) then
  5148. vshear(i) = 1.e3 * vshear(i) / (zi(i,ktcon(i))-zi(i,kb(i)))
  5149. e1=1.591-.639*vshear(i) &
  5150. & +.0953*(vshear(i)**2)-.00496*(vshear(i)**3)
  5151. edt(i)=1.-e1
  5152. val = .9
  5153. edt(i) = min(edt(i),val)
  5154. val = .0
  5155. edt(i) = max(edt(i),val)
  5156. endif
  5157. enddo
  5158. !c
  5159. !c--- what would the change be, that a cloud with unit mass
  5160. !c--- will do to the environment?
  5161. !c
  5162. do k = 1, km
  5163. do i = 1, im
  5164. if(cnvflg(i) .and. k .le. kmax(i)) then
  5165. dellah(i,k) = 0.
  5166. dellaq(i,k) = 0.
  5167. dellau(i,k) = 0.
  5168. dellav(i,k) = 0.
  5169. endif
  5170. enddo
  5171. enddo
  5172. !c
  5173. !c--- changed due to subsidence and entrainment
  5174. !c
  5175. do k = 2, km1
  5176. do i = 1, im
  5177. if (cnvflg(i)) then
  5178. if(k.gt.kb(i).and.k.lt.ktcon(i)) then
  5179. dp = 1000. * del(i,k)
  5180. dz = zi(i,k) - zi(i,k-1)
  5181. !c
  5182. dv1h = heo(i,k)
  5183. dv2h = .5 * (heo(i,k) + heo(i,k-1))
  5184. dv3h = heo(i,k-1)
  5185. dv1q = qo(i,k)
  5186. dv2q = .5 * (qo(i,k) + qo(i,k-1))
  5187. dv3q = qo(i,k-1)
  5188. dv1u = uo(i,k)
  5189. dv2u = .5 * (uo(i,k) + uo(i,k-1))
  5190. dv3u = uo(i,k-1)
  5191. dv1v = vo(i,k)
  5192. dv2v = .5 * (vo(i,k) + vo(i,k-1))
  5193. dv3v = vo(i,k-1)
  5194. !c
  5195. tem = 0.5 * (xlamue(i,k)+xlamue(i,k-1))
  5196. tem1 = xlamud(i)
  5197. !cj
  5198. dellah(i,k) = dellah(i,k) + &
  5199. & ( eta(i,k)*dv1h - eta(i,k-1)*dv3h &
  5200. & - tem*eta(i,k-1)*dv2h*dz &
  5201. & + tem1*eta(i,k-1)*.5*(hcko(i,k)+hcko(i,k-1))*dz &
  5202. & ) *g/dp
  5203. !cj
  5204. dellaq(i,k) = dellaq(i,k) + &
  5205. & ( eta(i,k)*dv1q - eta(i,k-1)*dv3q &
  5206. & - tem*eta(i,k-1)*dv2q*dz &
  5207. & + tem1*eta(i,k-1)*.5*(qcko(i,k)+qcko(i,k-1))*dz &
  5208. & ) *g/dp
  5209. !cj
  5210. dellau(i,k) = dellau(i,k) + &
  5211. & ( eta(i,k)*dv1u - eta(i,k-1)*dv3u &
  5212. & - tem*eta(i,k-1)*dv2u*dz &
  5213. & + tem1*eta(i,k-1)*.5*(ucko(i,k)+ucko(i,k-1))*dz &
  5214. & - pgcon*eta(i,k-1)*(dv1u-dv3u) &
  5215. & ) *g/dp
  5216. !cj
  5217. dellav(i,k) = dellav(i,k) + &
  5218. & ( eta(i,k)*dv1v - eta(i,k-1)*dv3v &
  5219. & - tem*eta(i,k-1)*dv2v*dz &
  5220. & + tem1*eta(i,k-1)*.5*(vcko(i,k)+vcko(i,k-1))*dz &
  5221. & - pgcon*eta(i,k-1)*(dv1v-dv3v) &
  5222. & ) *g/dp
  5223. !cj
  5224. endif
  5225. endif
  5226. enddo
  5227. enddo
  5228. !c
  5229. !c------- cloud top
  5230. !c
  5231. do i = 1, im
  5232. if(cnvflg(i)) then
  5233. indx = ktcon(i)
  5234. dp = 1000. * del(i,indx)
  5235. dv1h = heo(i,indx-1)
  5236. dellah(i,indx) = eta(i,indx-1) * &
  5237. & (hcko(i,indx-1) - dv1h) * g / dp
  5238. dv1q = qo(i,indx-1)
  5239. dellaq(i,indx) = eta(i,indx-1) * &
  5240. & (qcko(i,indx-1) - dv1q) * g / dp
  5241. dv1u = uo(i,indx-1)
  5242. dellau(i,indx) = eta(i,indx-1) * &
  5243. & (ucko(i,indx-1) - dv1u) * g / dp
  5244. dv1v = vo(i,indx-1)
  5245. dellav(i,indx) = eta(i,indx-1) * &
  5246. & (vcko(i,indx-1) - dv1v) * g / dp
  5247. !c
  5248. !c cloud water
  5249. !c
  5250. dellal(i,indx) = eta(i,indx-1) * &
  5251. & qlko_ktcon(i) * g / dp
  5252. endif
  5253. enddo
  5254. !c
  5255. !c mass flux at cloud base for shallow convection
  5256. !c (Grant, 2001)
  5257. !c
  5258. do i= 1, im
  5259. if(cnvflg(i)) then
  5260. k = kbcon(i)
  5261. ! ptem = g*sflx(i)*zi(i,k)/t1(i,1)
  5262. ptem = g*sflx(i)*hpbl(i)/t1(i,1)
  5263. wstar(i) = ptem**h1
  5264. tem = po(i,k)*100. / (rd*t1(i,k))
  5265. xmb(i) = betaw*tem*wstar(i)
  5266. xmb(i) = min(xmb(i),xmbmax(i))
  5267. endif
  5268. enddo
  5269. !c
  5270. !c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  5271. !c
  5272. do k = 1, km
  5273. do i = 1, im
  5274. if (cnvflg(i) .and. k .le. kmax(i)) then
  5275. qeso(i,k) = 0.01 * fpvs(t1(i,k)) ! fpvs is in pa
  5276. qeso(i,k) = eps * qeso(i,k) / (pfld(i,k) + epsm1*qeso(i,k))
  5277. val = 1.e-8
  5278. qeso(i,k) = max(qeso(i,k), val )
  5279. endif
  5280. enddo
  5281. enddo
  5282. !c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  5283. !c
  5284. do i = 1, im
  5285. delhbar(i) = 0.
  5286. delqbar(i) = 0.
  5287. deltbar(i) = 0.
  5288. delubar(i) = 0.
  5289. delvbar(i) = 0.
  5290. qcond(i) = 0.
  5291. enddo
  5292. do k = 1, km
  5293. do i = 1, im
  5294. if (cnvflg(i)) then
  5295. if(k.gt.kb(i).and.k.le.ktcon(i)) then
  5296. dellat = (dellah(i,k) - hvap * dellaq(i,k)) / cp
  5297. t1(i,k) = t1(i,k) + dellat * xmb(i) * dt2
  5298. q1(i,k) = q1(i,k) + dellaq(i,k) * xmb(i) * dt2
  5299. tem = 1./rcs(i)
  5300. u1(i,k) = u1(i,k) + dellau(i,k) * xmb(i) * dt2 * tem
  5301. v1(i,k) = v1(i,k) + dellav(i,k) * xmb(i) * dt2 * tem
  5302. dp = 1000. * del(i,k)
  5303. delhbar(i) = delhbar(i) + dellah(i,k)*xmb(i)*dp/g
  5304. delqbar(i) = delqbar(i) + dellaq(i,k)*xmb(i)*dp/g
  5305. deltbar(i) = deltbar(i) + dellat*xmb(i)*dp/g
  5306. delubar(i) = delubar(i) + dellau(i,k)*xmb(i)*dp/g
  5307. delvbar(i) = delvbar(i) + dellav(i,k)*xmb(i)*dp/g
  5308. endif
  5309. endif
  5310. enddo
  5311. enddo
  5312. do k = 1, km
  5313. do i = 1, im
  5314. if (cnvflg(i)) then
  5315. if(k.gt.kb(i).and.k.le.ktcon(i)) then
  5316. qeso(i,k) = 0.01 * fpvs(t1(i,k)) ! fpvs is in pa
  5317. qeso(i,k) = eps * qeso(i,k)/(pfld(i,k) + epsm1*qeso(i,k))
  5318. val = 1.e-8
  5319. qeso(i,k) = max(qeso(i,k), val )
  5320. endif
  5321. endif
  5322. enddo
  5323. enddo
  5324. !c
  5325. do i = 1, im
  5326. rntot(i) = 0.
  5327. delqev(i) = 0.
  5328. delq2(i) = 0.
  5329. flg(i) = cnvflg(i)
  5330. enddo
  5331. do k = km, 1, -1
  5332. do i = 1, im
  5333. if (cnvflg(i)) then
  5334. if(k.lt.ktcon(i).and.k.gt.kb(i)) then
  5335. rntot(i) = rntot(i) + pwo(i,k) * xmb(i) * .001 * dt2
  5336. endif
  5337. endif
  5338. enddo
  5339. enddo
  5340. !c
  5341. !c evaporating rain
  5342. !c
  5343. do k = km, 1, -1
  5344. do i = 1, im
  5345. if (k .le. kmax(i)) then
  5346. deltv(i) = 0.
  5347. delq(i) = 0.
  5348. qevap(i) = 0.
  5349. if(cnvflg(i)) then
  5350. if(k.lt.ktcon(i).and.k.gt.kb(i)) then
  5351. rn(i) = rn(i) + pwo(i,k) * xmb(i) * .001 * dt2
  5352. endif
  5353. endif
  5354. if(flg(i).and.k.lt.ktcon(i)) then
  5355. evef = edt(i) * evfact
  5356. if(slimsk(i).eq.1.) evef=edt(i) * evfactl
  5357. ! if(slimsk(i).eq.1.) evef=.07
  5358. !c if(slimsk(i).ne.1.) evef = 0.
  5359. qcond(i) = evef * (q1(i,k) - qeso(i,k)) &
  5360. & / (1. + el2orc * qeso(i,k) / t1(i,k)**2)
  5361. dp = 1000. * del(i,k)
  5362. if(rn(i).gt.0..and.qcond(i).lt.0.) then
  5363. qevap(i) = -qcond(i) * (1.-exp(-.32*sqrt(dt2*rn(i))))
  5364. qevap(i) = min(qevap(i), rn(i)*1000.*g/dp)
  5365. delq2(i) = delqev(i) + .001 * qevap(i) * dp / g
  5366. endif
  5367. if(rn(i).gt.0..and.qcond(i).lt.0..and. &
  5368. & delq2(i).gt.rntot(i)) then
  5369. qevap(i) = 1000.* g * (rntot(i) - delqev(i)) / dp
  5370. flg(i) = .false.
  5371. endif
  5372. if(rn(i).gt.0..and.qevap(i).gt.0.) then
  5373. tem = .001 * dp / g
  5374. tem1 = qevap(i) * tem
  5375. if(tem1.gt.rn(i)) then
  5376. qevap(i) = rn(i) / tem
  5377. rn(i) = 0.
  5378. else
  5379. rn(i) = rn(i) - tem1
  5380. endif
  5381. q1(i,k) = q1(i,k) + qevap(i)
  5382. t1(i,k) = t1(i,k) - elocp * qevap(i)
  5383. deltv(i) = - elocp*qevap(i)/dt2
  5384. delq(i) = + qevap(i)/dt2
  5385. delqev(i) = delqev(i) + .001*dp*qevap(i)/g
  5386. endif
  5387. dellaq(i,k) = dellaq(i,k) + delq(i) / xmb(i)
  5388. delqbar(i) = delqbar(i) + delq(i)*dp/g
  5389. deltbar(i) = deltbar(i) + deltv(i)*dp/g
  5390. endif
  5391. endif
  5392. enddo
  5393. enddo
  5394. !cj
  5395. ! do i = 1, im
  5396. ! if(me.eq.31.and.cnvflg(i)) then
  5397. ! if(cnvflg(i)) then
  5398. ! print *, ' shallow delhbar, delqbar, deltbar = ',
  5399. ! & delhbar(i),hvap*delqbar(i),cp*deltbar(i)
  5400. ! print *, ' shallow delubar, delvbar = ',delubar(i),delvbar(i)
  5401. ! print *, ' precip =', hvap*rn(i)*1000./dt2
  5402. ! print*,'pdif= ',pfld(i,kbcon(i))-pfld(i,ktcon(i))
  5403. ! endif
  5404. ! enddo
  5405. !cj
  5406. do i = 1, im
  5407. if(cnvflg(i)) then
  5408. if(rn(i).lt.0..or..not.flg(i)) rn(i) = 0.
  5409. ktop(i) = ktcon(i)
  5410. kbot(i) = kbcon(i)
  5411. kcnv(i) = 0
  5412. endif
  5413. enddo
  5414. !c
  5415. !c cloud water
  5416. !c
  5417. if (ncloud.gt.0) then
  5418. !
  5419. val1 = 1.0
  5420. val2 = 0.
  5421. do k = 1, km1
  5422. do i = 1, im
  5423. if (cnvflg(i)) then
  5424. if (k.gt.kb(i).and.k.le.ktcon(i)) then
  5425. tem = dellal(i,k) * xmb(i) * dt2
  5426. ! tem1 = max(0.0, min(1.0, (tcr-t1(i,k))*tcrf))
  5427. tem1 = max(val2, min(val1, (tcr-t1(i,k))*tcrf))
  5428. if (ql(i,k,2) .gt. -999.0) then
  5429. ql(i,k,1) = ql(i,k,1) + tem * tem1 ! ice
  5430. ql(i,k,2) = ql(i,k,2) + tem *(1.0-tem1) ! water
  5431. else
  5432. ql(i,k,1) = ql(i,k,1) + tem
  5433. endif
  5434. endif
  5435. endif
  5436. enddo
  5437. enddo
  5438. !
  5439. endif
  5440. !
  5441. ! hchuang code change
  5442. !
  5443. do k = 1, km
  5444. do i = 1, im
  5445. if(cnvflg(i)) then
  5446. if(k.ge.kb(i) .and. k.lt.ktop(i)) then
  5447. ud_mf(i,k) = eta(i,k) * xmb(i) * dt2
  5448. endif
  5449. endif
  5450. enddo
  5451. enddo
  5452. do i = 1, im
  5453. if(cnvflg(i)) then
  5454. k = ktop(i)-1
  5455. dt_mf(i,k) = ud_mf(i,k)
  5456. endif
  5457. enddo
  5458. !!
  5459. return
  5460. end subroutine shalcnv
  5461. END MODULE module_cu_sas