PageRenderTime 1385ms CodeModel.GetById 39ms RepoModel.GetById 1ms app.codeStats 1ms

/wrfv2_fire/phys/module_cu_kfeta.F

http://github.com/jbeezley/wrf-fire
FORTRAN Legacy | 3174 lines | 2002 code | 116 blank | 1056 comment | 19 complexity | 83c977654d6529edf366d4c5881efc42 MD5 | raw file
Possible License(s): AGPL-1.0

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

    

  1. MODULE module_cu_kfeta
    2. 

  2. 

  3.    USE module_wrf_error
    4. 

  4. 

  5. !
    6. 

  6. !  V3.3: A new trigger function is added based Ma and Tan (2009):
    7. 

  7. !   Ma, L.-M. and Z.-M. Tan, 2009: Improving the behavior of 
    8. 

  8. !   the cumulus parameterization for tropical cyclone prediction: 
    9. 

  9. !   Convection trigger. Atmospheric Research, 92, 190 - 211.
    10. 

  10. !
    11. 

  11. !--------------------------------------------------------------------
    12. 

  12. ! Lookup table variables:
    13. 

  13.       INTEGER, PARAMETER :: KFNT=250,KFNP=220
    14. 

  14.       REAL, DIMENSION(KFNT,KFNP),PRIVATE, SAVE :: TTAB,QSTAB
    15. 

  15.       REAL, DIMENSION(KFNP),PRIVATE, SAVE :: THE0K
    16. 

  16.       REAL, DIMENSION(200),PRIVATE, SAVE :: ALU
    17. 

  17.       REAL, PRIVATE, SAVE :: RDPR,RDTHK,PLUTOP
    18. 

  18. ! Note:  KF Lookup table is used by subroutines KF_eta_PARA, TPMIX2,
    19. 

  19. !        TPMIX2DD, ENVIRTHT
    20. 

  20. ! End of Lookup table variables:
    21. 

  21. 

  22. CONTAINS
    23. 

  23. 

  24.    SUBROUTINE KF_eta_CPS(                                    &
    25. 

  25.               ids,ide, jds,jde, kds,kde                      &
    26. 

  26.              ,ims,ime, jms,jme, kms,kme                      &
    27. 

  27.              ,its,ite, jts,jte, kts,kte                      &
    28. 

  28.              ,trigger                                        &
    29. 

  29.              ,DT,KTAU,DX,CUDT,CURR_SECS,ADAPT_STEP_FLAG      &
    30. 

  30.              ,CUDTACTTIME                                    & 
    31. 

  31.              ,rho,RAINCV,PRATEC,NCA                          &
    32. 

  32.              ,U,V,TH,T,W,dz8w,Pcps,pi                        &
    33. 

  33.              ,W0AVG,XLV0,XLV1,XLS0,XLS1,CP,R,G,EP1           &
    34. 

  34.              ,EP2,SVP1,SVP2,SVP3,SVPT0                       &
    35. 

  35.              ,STEPCU,CU_ACT_FLAG,warm_rain,CUTOP,CUBOT       &
    36. 

  36.              ,QV                                             &
    37. 

  37.             ! optionals
    38. 

  38.              ,F_QV    ,F_QC    ,F_QR    ,F_QI    ,F_QS       &
    39. 

  39.              ,RTHCUTEN,RQVCUTEN,RQCCUTEN,RQRCUTEN            &
    40. 

  40.              ,RQICUTEN,RQSCUTEN, RQVFTEN                     &
    41. 

  41.                                                              )
    42. 

  42. !
    43. 

  43. !-------------------------------------------------------------
    44. 

  44.    IMPLICIT NONE
    45. 

  45. !-------------------------------------------------------------
    46. 

  46.    INTEGER,      INTENT(IN   ) ::                            &
    47. 

  47.                                   ids,ide, jds,jde, kds,kde, &
    48. 

  48.                                   ims,ime, jms,jme, kms,kme, &
    49. 

  49.                                   its,ite, jts,jte, kts,kte
    50. 

  50. 

  51.    INTEGER,      INTENT(IN   ) :: trigger
    52. 

  52.    INTEGER,      INTENT(IN   ) :: STEPCU
    53. 

  53.    LOGICAL,      INTENT(IN   ) :: warm_rain
    54. 

  54. 

  55.    REAL,         INTENT(IN   ) :: XLV0,XLV1,XLS0,XLS1
    56. 

  56.    REAL,         INTENT(IN   ) :: CP,R,G,EP1,EP2
    57. 

  57.    REAL,         INTENT(IN   ) :: SVP1,SVP2,SVP3,SVPT0
    58. 

  58. 

  59.    INTEGER,      INTENT(IN   ) :: KTAU           
    60. 

  60. 

  61.    REAL,  DIMENSION( ims:ime , kms:kme , jms:jme )         , &
    62. 

  62.           INTENT(IN   ) ::                                   &
    63. 

  63.                                                           U, &
    64. 

  64.                                                           V, &
    65. 

  65.                                                           W, &
    66. 

  66.                                                          TH, &
    67. 

  67.                                                           T, &
    68. 

  68.                                                          QV, &
    69. 

  69.                                                        dz8w, &
    70. 

  70.                                                        Pcps, &
    71. 

  71.                                                         rho, &
    72. 

  72.                                                          pi
    73. 

  73. !
    74. 

  74.    REAL,  DIMENSION( ims:ime , kms:kme , jms:jme )         , &
    75. 

  75.           INTENT(INOUT) ::                                   &
    76. 

  76.                                                       W0AVG
    77. 

  77. 

  78.    REAL,  INTENT(IN   ) :: DT, DX
    79. 

  79.    REAL,  INTENT(IN   ) :: CUDT
    80. 

  80.    REAL,  INTENT(IN   ) :: CURR_SECS
    81. 

  81.    LOGICAL,OPTIONAL,INTENT(IN   ) :: ADAPT_STEP_FLAG
    82. 

  82.    REAL,  INTENT (INOUT) :: CUDTACTTIME       
    83. 

  83. !
    84. 

  84.    REAL, DIMENSION( ims:ime , jms:jme ),                     &
    85. 

  85.           INTENT(INOUT) ::                           RAINCV
    86. 

  86. 

  87.    REAL,    DIMENSION( ims:ime , jms:jme ),                  &
    88. 

  88.           INTENT(INOUT) ::                           PRATEC
    89. 

  89. 

  90.    REAL,    DIMENSION( ims:ime , jms:jme ),                  &
    91. 

  91.             INTENT(INOUT) ::                            NCA
    92. 

  92. 

  93.    REAL, DIMENSION( ims:ime , jms:jme ),                     &
    94. 

  94.           INTENT(OUT) ::                              CUBOT, &
    95. 

  95.                                                       CUTOP    
    96. 

  96. 

  97.    LOGICAL, DIMENSION( ims:ime , jms:jme ),                  &
    98. 

  98.           INTENT(INOUT) :: CU_ACT_FLAG
    99. 

  99. 

  100. !
    101. 

  101. ! Optional arguments
    102. 

  102. !
    103. 

  103. 

  104.    REAL, DIMENSION( ims:ime , kms:kme , jms:jme ),           &
    105. 

  105.          OPTIONAL,                                           &
    106. 

  106.          INTENT(INOUT) ::                                    &
    107. 

  107.                                                    RTHCUTEN, &
    108. 

  108.                                                    RQVCUTEN, &
    109. 

  109.                                                    RQCCUTEN, &
    110. 

  110.                                                    RQRCUTEN, &
    111. 

  111.                                                    RQICUTEN, &
    112. 

  112.                                                    RQSCUTEN, &
    113. 

  113.                                                    RQVFTEN
    114. 

  114. !
    115. 

  115. ! Flags relating to the optional tendency arrays declared above
    116. 

  116. ! Models that carry the optional tendencies will provdide the
    117. 

  117. ! optional arguments at compile time; these flags all the model
    118. 

  118. ! to determine at run-time whether a particular tracer is in
    119. 

  119. ! use or not.
    120. 

  120. !
    121. 

  121.    LOGICAL, OPTIONAL ::                                      &
    122. 

  122.                                                    F_QV      &
    123. 

  123.                                                   ,F_QC      &
    124. 

  124.                                                   ,F_QR      &
    125. 

  125.                                                   ,F_QI      &
    126. 

  126.                                                   ,F_QS
    127. 

  127. 

  128. 

  129. ! LOCAL VARS
    130. 

  130. 

  131.    LOGICAL :: flag_qr, flag_qi, flag_qs
    132. 

  132. 

  133.    REAL, DIMENSION( kts:kte ) ::                             &
    134. 

  134.                                                         U1D, &
    135. 

  135.                                                         V1D, &
    136. 

  136.                                                         T1D, &
    137. 

  137.                                                        DZ1D, &
    138. 

  138.                                                        QV1D, &
    139. 

  139.                                                         P1D, &
    140. 

  140.                                                       RHO1D, &
    141. 

  141.                                                   tpart_v1D, &
    142. 

  142.                                                   tpart_h1D, &
    143. 

  143.                                                     W0AVG1D
    144. 

  144. 

  145.    REAL, DIMENSION( kts:kte )::                              &
    146. 

  146.                                                        DQDT, &
    147. 

  147.                                                       DQIDT, &
    148. 

  148.                                                       DQCDT, &
    149. 

  149.                                                       DQRDT, &
    150. 

  150.                                                       DQSDT, &
    151. 

  151.                                                        DTDT
    152. 

  152. 

  153.   REAL, DIMENSION (its-1:ite+1,kts:kte,jts-1:jte+1) ::  aveh_t, aveh_q
    154. 

  154.   REAL, DIMENSION (its:ite,kts:kte,jts:jte) :: aveh_qmax, aveh_qmin
    155. 

  155.   REAL, DIMENSION (its:ite,kts:kte,jts:jte) ::  avev_t, avev_q 
    156. 

  156.   REAL, DIMENSION (its:ite,kts:kte,jts:jte) :: avev_qmax, avev_qmin
    157. 

  157.   REAL, DIMENSION (its:ite,kts:kte,jts:jte) ::  coef_v, coef_h, tpart_h, tpart_v
    158. 

  158.   INTEGER :: ii,jj,kk
    159. 

  159. 

  160.   REAL :: ttop
    161. 

  161.   REAL, DIMENSION (kts:kte)  :: z0
    162. 

  162. 

  163.    REAL    ::         TST,tv,PRS,RHOE,W0,SCR1,DXSQ,tmp
    164. 

  164.    integer :: ibegh,iendh,jbegh,jendh
    165. 

  165.    integer :: istart,iend,jstart,jend
    166. 

  166.    INTEGER :: i,j,k,NTST
    167. 

  167.    REAL    :: lastdt = -1.0
    168. 

  168.    REAL    :: W0AVGfctr, W0fctr, W0den
    169. 

  169.    LOGICAL :: run_param , doing_adapt_dt , decided
    170. 

  170.    
    171. 

  171. !
    172. 

  172.    DXSQ=DX*DX
    173. 

  173. 

  174. !----------------------
    175. 

  175.    NTST=STEPCU
    176. 

  176.    TST=float(NTST*2)
    177. 

  177.    flag_qr = .FALSE.
    178. 

  178.    flag_qi = .FALSE.
    179. 

  179.    flag_qs = .FALSE.
    180. 

  180.    IF ( PRESENT(F_QR) ) flag_qr = F_QR
    181. 

  181.    IF ( PRESENT(F_QI) ) flag_qi = F_QI
    182. 

  182.    IF ( PRESENT(F_QS) ) flag_qs = F_QS
    183. 

  183. !
    184. 

  184.    if (lastdt < 0) then
    185. 

  185.       lastdt = dt
    186. 

  186.    endif
    187. 

  187.    
    188. 

  188.    if (ADAPT_STEP_FLAG) then
    189. 

  189.       W0AVGfctr = 2 * MAX(CUDT*60,dt) - dt
    190. 

  190.       W0fctr = dt
    191. 

  191.       W0den = 2 * MAX(CUDT*60,dt)
    192. 

  192.    else
    193. 

  193.       W0AVGfctr = (TST-1.)
    194. 

  194.       W0fctr = 1.
    195. 

  195.       W0den = TST
    196. 

  196.    endif
    197. 

  197. 

  198.   DO J = jts,jte
    199. 

  199.       DO K=kts,kte
    200. 

  200.          DO I= its,ite
    201. 

  201. !            SCR1=-5.0E-4*G*rho(I,K,J)*(w(I,K,J)+w(I,K+1,J))
    202. 

  202. !            TV=T(I,K,J)*(1.+EP1*QV(I,K,J))
    203. 

  203. !            RHOE=Pcps(I,K,J)/(R*TV)
    204. 

  204. !            W0=-101.9368*SCR1/RHOE
    205. 

  205.             W0=0.5*(w(I,K,J)+w(I,K+1,J))
    206. 

  206. 

  207. !           Old:            
    208. 

  208. !
    209. 

  209. !            W0AVG(I,K,J)=(W0AVG(I,K,J)*(TST-1.)+W0)/TST            
    210. 

  210. !
    211. 

  211. !           New, to support adaptive time step:
    212. 

  212. !
    213. 

  213.             W0AVG(I,K,J) = ( W0AVG(I,K,J) * W0AVGfctr + W0 * W0fctr ) / W0den
    214. 

  214.          ENDDO
    215. 

  215.       ENDDO
    216. 

  216.    ENDDO
    217. 

  217.    lastdt = dt
    218. 

  218. 

  219. !  Initialization for adaptive time step.
    220. 

  220. 

  221.    doing_adapt_dt = .FALSE.
    222. 

  222.    IF ( PRESENT(adapt_step_flag) ) THEN
    223. 

  223.       IF ( adapt_step_flag ) THEN
    224. 

  224.          doing_adapt_dt = .TRUE.
    225. 

  225.          IF ( cudtacttime .EQ. 0. ) THEN
    226. 

  226.             cudtacttime = curr_secs + cudt*60.
    227. 

  227.          END IF
    228. 

  228.       END IF
    229. 

  229.    END IF
    230. 

  230. 

  231. !  Do we run through this scheme or not?
    232. 

  232. 

  233. !    Test 1:  If this is the initial model time, then yes.
    234. 

  234. !                KTAU=1
    235. 

  235. !    Test 2:  If the user asked for the cumulus to be run every time step, then yes.
    236. 

  236. !                CUDT=0 or STEPCU=1
    237. 

  237. !    Test 3:  If not adaptive dt, and this is on the requested cumulus frequency, then yes.
    238. 

  238. !                MOD(KTAU,NST)=0
    239. 

  239. !    Test 4:  If using adaptive dt and the current time is past the last requested activate cumulus time, then yes.
    240. 

  240. !                CURR_SECS >= CUDTACTTIME
    241. 

  241. 

  242. !  If we do run through the scheme, we set the flag run_param to TRUE and we set the decided flag
    243. 

  243. !  to TRUE.  The decided flag says that one of these tests was able to say "yes", run the scheme.
    244. 

  244. !  We only proceed to other tests if the previous tests all have left decided as FALSE.
    245. 

  245. 

  246. !  If we set run_param to TRUE and this is adaptive time stepping, we set the time to the next
    247. 

  247. !  cumulus run.
    248. 

  248. 

  249.    decided = .FALSE.
    250. 

  250.    run_param = .FALSE.
    251. 

  251.    IF ( ( .NOT. decided ) .AND. &
    252. 

  252.         ( ktau .EQ. 1 ) ) THEN
    253. 

  253.       run_param   = .TRUE.
    254. 

  254.       decided     = .TRUE.
    255. 

  255.    END IF
    256. 

  256. 

  257.    IF ( ( .NOT. decided ) .AND. &
    258. 

  258.         ( ( cudt .EQ. 0. ) .OR. ( stepcu .EQ. 1 ) ) ) THEN
    259. 

  259.       run_param   = .TRUE.
    260. 

  260.       decided     = .TRUE.
    261. 

  261.    END IF
    262. 

  262. 

  263.    IF ( ( .NOT. decided ) .AND. &
    264. 

  264.         ( .NOT. doing_adapt_dt ) .AND. &
    265. 

  265.         ( MOD(ktau,ntst) .EQ. 0 ) ) THEN
    266. 

  266.       run_param   = .TRUE.
    267. 

  267.       decided     = .TRUE.
    268. 

  268.    END IF
    269. 

  269. 

  270.    IF ( ( .NOT. decided ) .AND. &
    271. 

  271.         ( doing_adapt_dt ) .AND. &
    272. 

  272.         ( curr_secs .GE. cudtacttime ) ) THEN
    273. 

  273.       run_param   = .TRUE.
    274. 

  274.       decided     = .TRUE.
    275. 

  275.       cudtacttime = curr_secs + cudt*60
    276. 

  276.    END IF
    277. 

  277. 

  278.    if (run_param) then
    279. 

  279. 

  280. ! New trigger function
    281. 

  281.      IF (trigger.eq.2) THEN         
    282. 

  282. !
    283. 

  283. ! calculate 9-point average of moisture advection and temperature using halo (Horizontal)
    284. 

  284. !
    285. 

  285.      aveh_t=-999   ! horizontal 9-point ave
    286. 

  286.      aveh_q=-999
    287. 

  287.      avev_t=0   ! vertical 3-level ave
    288. 

  288.      avev_q=0
    289. 

  289.      avev_qmax=0
    290. 

  290.      avev_qmin=0
    291. 

  291.      aveh_qmax=0
    292. 

  292.      aveh_qmin=0
    293. 

  293.      tpart_h=0
    294. 

  294.      tpart_v=0
    295. 

  295.      coef_h=0
    296. 

  296.      coef_v=0
    297. 

  297.      ibegh=max(its-1, ids+1)   ! start from 2
    298. 

  298.      jbegh=max(jts-1, jds+1)
    299. 

  299.      iendh=min(ite+1, ide-2)   ! end at ide-2
    300. 

  300.      jendh=min(jte+1, jde-2)
    301. 

  301.         DO J = jbegh,jendh
    302. 

  302.         DO K = kts,kte
    303. 

  303.         DO I = ibegh,iendh
    304. 

  304.           aveh_t(i,k,j)=(T(i-1,k,j-1)+T(i-1,k,j)  +T(i-1,k,j+1)+ &
    305. 

  305.                          T(i,k,j-1)   +T(i,k,j)   +T(i,k,j+1)+         &
    306. 

  306.                          T(i+1,k,j-1) +T(i+1,k,j) +T(i+1,k,j+1))/9.
    307. 

  307.           aveh_q(i,k,j)=(rqvften(i-1,k,j-1)+rqvften(i-1,k,j)  +rqvften(i-1,k,j+1)+ &
    308. 

  308.                          rqvften(i,k,j-1)   +rqvften(i,k,j)   +rqvften(i,k,j+1)+         &
    309. 

  309.                          rqvften(i+1,k,j-1) +rqvften(i+1,k,j) +rqvften(i+1,k,j+1))/9.
    310. 

  310.         ENDDO
    311. 

  311.         ENDDO
    312. 

  312.         ENDDO
    313. 

  313. ! boundary value ( all processors will do the following? Or just those processsors handling sub-area including boundary)
    314. 

  314.         DO K = kts,kte
    315. 

  315.            DO J = jts-1,jte+1
    316. 

  316.             DO I = its-1,ite+1
    317. 

  317. 

  318.             if(i.eq.ids) then
    319. 

  319.             aveh_t(i,k,j)=aveh_t(i+1,k,j)
    320. 

  320.             aveh_q(i,k,j)=aveh_q(i+1,k,j)
    321. 

  321.             elseif(i.eq.ide-1) then
    322. 

  322.             aveh_t(i,k,j)=aveh_t(i-1,k,j)
    323. 

  323.             aveh_q(i,k,j)=aveh_q(i-1,k,j)
    324. 

  324.             endif
    325. 

  325. 

  326.             if(j.eq.jds) then
    327. 

  327.              aveh_t(i,k,j)=aveh_t(i,k,j+1)
    328. 

  328.              aveh_q(i,k,j)=aveh_q(i,k,j+1)
    329. 

  329.             elseif(j.eq.jde-1) then
    330. 

  330.             aveh_t(i,k,j)=aveh_t(i,k,j-1)
    331. 

  331.             aveh_q(i,k,j)=aveh_q(i,k,j-1)
    332. 

  332.             endif
    333. 

  333. 

  334.             if(j.eq.jds.and.i.eq.ids) then
    335. 

  335.             aveh_q(i,k,j)=aveh_q(i+1,k,j+1) 
    336. 

  336.             aveh_t(i,k,j)=aveh_t(i+1,k,j+1) 
    337. 

  337.             endif
    338. 

  338. 

  339.             if(j.eq.jde-1.and.i.eq.ids) then
    340. 

  340.             aveh_q(i,k,j)=aveh_q(i+1,k,j-1) 
    341. 

  341.             aveh_t(i,k,j)=aveh_t(i+1,k,j-1) 
    342. 

  342.             endif
    343. 

  343. 

  344.             if(j.eq.jde-1.and.i.eq.ide-1) then
    345. 

  345.             aveh_q(i,k,j)=aveh_q(i-1,k,j-1) 
    346. 

  346.             aveh_t(i,k,j)=aveh_t(i-1,k,j-1) 
    347. 

  347.             endif
    348. 

  348. 

  349.             if(j.eq.jds.and.i.eq.ide-1) then
    350. 

  350.             aveh_q(i,k,j)=aveh_q(i-1,k,j+1) 
    351. 

  351.             aveh_t(i,k,j)=aveh_t(i-1,k,j+1) 
    352. 

  352.             endif
    353. 

  353. 

  354.             ENDDO
    355. 

  355.            ENDDO
    356. 

  356.         ENDDO
    357. 

  357. ! search for max/min moisture advection in 9-point range, calculate horizontal T-perturbation (tpart_h)
    358. 

  358.      istart=max(its, ids+1)   ! start from 2
    359. 

  359.      jstart=max(jts, jds+1)
    360. 

  360.      iend=min(ite, ide-2)   ! end at ide-2
    361. 

  361.      jend=min(jte, jde-2)
    362. 

  362.         DO K = kts,kte
    363. 

  363.         DO J = jstart,jend
    364. 

  364.         DO I = istart,iend
    365. 

  365.            aveh_qmax(i,k,j)=aveh_q(i,k,j)
    366. 

  366.            aveh_qmin(i,k,j)=aveh_q(i,k,j)
    367. 

  367.           DO ii=-1, 1
    368. 

  368.            DO jj=-1,1
    369. 

  369.              if(aveh_q(i+II,k,j+JJ).gt.aveh_qmax(i,k,j)) aveh_qmax(i,k,j)=aveh_q(i+II,k,j+JJ)
    370. 

  370.              if(aveh_q(i+II,k,j+JJ).lt.aveh_qmin(i,k,j)) aveh_qmin(i,k,j)=aveh_q(i+II,k,j+JJ)
    371. 

  371.            ENDDO
    372. 

  372.           ENDDO 
    373. 

  373.           if(aveh_qmax(i,k,j).gt.aveh_qmin(i,k,j))then
    374. 

  374.           coef_h(i,k,j)=(aveh_q(i,k,j)-aveh_qmin(i,k,j))/(aveh_qmax(i,k,j)-aveh_qmin(i,k,j))
    375. 

  375.           else
    376. 

  376.           coef_h(i,k,j)=0.
    377. 

  377.           endif
    378. 

  378.           coef_h(i,k,j)=amin1(coef_h(i,k,j),1.0)
    379. 

  379.           coef_h(i,k,j)=amax1(coef_h(i,k,j),0.0)
    380. 

  380.           tpart_h(i,k,j)=coef_h(i,k,j)*(T(i,k,j)-aveh_t(i,k,j))
    381. 

  381.         ENDDO
    382. 

  382.         ENDDO
    383. 

  383.         ENDDO
    384. 

  384.     89 continue 
    385. 

  385. ! vertical 3-layer calculation
    386. 

  386.         DO J = jts, jte
    387. 

  387.         DO I = its, ite
    388. 

  388.           z0(1) = 0.5 * dz8w(i,1,j)
    389. 

  389.           DO K = 2, kte
    390. 

  390.             Z0(K) = Z0(K-1) + .5 * (DZ8W(i,K,j) + DZ8W(i,K-1,j))
    391. 

  391.           ENDDO
    392. 

  392.         DO K = kts+1,kte-1
    393. 

  393. 	  ttop = t(i,k,j) + ((t(i,k,j) - t(i,k+1,j)) / (z0(k) - z0(k+1))) * (z0(k)-z0(k-1))
    394. 

  394.           avev_t(i,k,j)=(T(i,k-1,j) + T(i,k,j) + ttop)/3.
    395. 

  395. !         avev_t(i,k,j)=(T(i,k-1,j)+T(i,k,j) + T(i,k+1,j))/3.
    396. 

  396.           avev_q(i,k,j)=(rqvften(i,k-1,j)+rqvften(i,k,j) + rqvften(i,k+1,j))/3.
    397. 

  397.         ENDDO
    398. 

  398.           avev_t(i,kts,j)=avev_t(i,kts+1,j)   ! lowest level value, is it the same as avev_t(i,kds,j)=avev_t(i,kds+1,j)?
    399. 

  399.           avev_q(i,kts,j)=avev_q(i,kts+1,j)   
    400. 

  400.           avev_t(i,kte,j)=avev_t(i,kte-1,j)   ! highest level value
    401. 

  401.           avev_q(i,kte,j)=avev_q(i,kte-1,j)   
    402. 

  402.         ENDDO
    403. 

  403.         ENDDO
    404. 

  404. ! max /min value
    405. 

  405.         DO J = jts, jte
    406. 

  406.         DO I = its, ite
    407. 

  407.         DO K = kts+1,kte-1
    408. 

  408.           avev_qmax(i,k,j)=avev_q(i,k,j)
    409. 

  409.           avev_qmin(i,k,j)=avev_q(i,k,j)
    410. 

  410.          DO kk=-1,1
    411. 

  411.          if(avev_q(i,k+kk,j).gt.avev_qmax(i,k,j)) avev_qmax(i,k,j)=avev_q(i,k+kk,j) 
    412. 

  412.          if(avev_q(i,k+kk,j).lt.avev_qmin(i,k,j)) avev_qmin(i,k,j)=avev_q(i,k+kk,j) 
    413. 

  413.          ENDDO
    414. 

  414.          if(avev_qmax(i,k,j).gt.avev_qmin(i,k,j)) then
    415. 

  415.          coef_v(i,k,j)=(avev_q(i,k,j)-avev_qmin(i,k,j))/(avev_qmax(i,k,j)-avev_qmin(i,k,j))
    416. 

  416.          else
    417. 

  417.          coef_v(i,k,j)=0
    418. 

  418.          endif
    419. 

  419.          tpart_v(i,k,j)=coef_v(i,k,j)*(T(i,k,j)-avev_t(i,k,j))
    420. 

  420.         ENDDO
    421. 

  421.          tpart_v(i,kts,j)= tpart_v(i,kts+1,j)    ! lowest level
    422. 

  422.          tpart_v(i,kte,j)= tpart_v(i,kte-1,j)    ! highest level
    423. 

  423.         ENDDO
    424. 

  424.         ENDDO
    425. 

  425.      ENDIF       ! endif (trigger.eq.2)          
    426. 

  426. !
    427. 

  427.      DO J = jts,jte
    428. 

  428.      DO I= its,ite
    429. 

  429.         CU_ACT_FLAG(i,j) = .true.
    430. 

  430.      ENDDO
    431. 

  431.      ENDDO
    432. 

  432. 

  433.      DO J = jts,jte
    434. 

  434.        DO I=its,ite
    435. 

  435.           
    436. 

  436. 

  437.          IF ( NCA(I,J) .ge. 0.5*DT ) then
    438. 

  438.             CU_ACT_FLAG(i,j) = .false.
    439. 

  439.          ELSE
    440. 

  440. 

  441.             DO k=kts,kte
    442. 

  442.                DQDT(k)=0.
    443. 

  443.                DQIDT(k)=0.
    444. 

  444.                DQCDT(k)=0.
    445. 

  445.                DQRDT(k)=0.
    446. 

  446.                DQSDT(k)=0.
    447. 

  447.                DTDT(k)=0.
    448. 

  448.             ENDDO
    449. 

  449.             RAINCV(I,J)=0.
    450. 

  450.             CUTOP(I,J)=KTS
    451. 

  451.             CUBOT(I,J)=KTE+1
    452. 

  452.             PRATEC(I,J)=0.
    453. 

  453. !
    454. 

  454. ! assign vars from 3D to 1D
    455. 

  455. 

  456.             DO K=kts,kte
    457. 

  457.                U1D(K) =U(I,K,J)
    458. 

  458.                V1D(K) =V(I,K,J)
    459. 

  459.                T1D(K) =T(I,K,J)
    460. 

  460.                RHO1D(K) =rho(I,K,J)
    461. 

  461.                QV1D(K)=QV(I,K,J)
    462. 

  462.                P1D(K) =Pcps(I,K,J)
    463. 

  463.                W0AVG1D(K) =W0AVG(I,K,J)
    464. 

  464.                DZ1D(k)=dz8w(I,K,J)
    465. 

  465.                IF (trigger.eq.2) THEN
    466. 

  466.                   tpart_h1D(K) =tpart_h(I,K,J)
    467. 

  467.                   tpart_v1D(K) =tpart_v(I,K,J)
    468. 

  468.                ELSE
    469. 

  469.                   tpart_h1D(K) = 0.
    470. 

  470.                   tpart_v1D(K) = 0.
    471. 

  471.                ENDIF
    472. 

  472.             ENDDO
    473. 

  473.             CALL KF_eta_PARA(I, J,                  &
    474. 

  474.                  U1D,V1D,T1D,QV1D,P1D,DZ1D,W0AVG1D, &
    475. 

  475.                  tpart_h1D,tpart_v1D,               &
    476. 

  476.                  trigger,                           &
    477. 

  477.                  DT,DX,DXSQ,RHO1D,                  &
    478. 

  478.                  XLV0,XLV1,XLS0,XLS1,CP,R,G,        &
    479. 

  479.                  EP2,SVP1,SVP2,SVP3,SVPT0,          &
    480. 

  480.                  DQDT,DQIDT,DQCDT,DQRDT,DQSDT,DTDT, &
    481. 

  481.                  RAINCV,PRATEC,NCA,                 &
    482. 

  482.                  flag_QI,flag_QS,warm_rain,         &
    483. 

  483.                  CUTOP,CUBOT,CUDT,                  &
    484. 

  484.                  ids,ide, jds,jde, kds,kde,         &
    485. 

  485.                  ims,ime, jms,jme, kms,kme,         &
    486. 

  486.                  its,ite, jts,jte, kts,kte)
    487. 

  487.             IF(PRESENT(rthcuten).AND.PRESENT(rqvcuten)) THEN
    488. 

  488.               DO K=kts,kte
    489. 

  489.                  RTHCUTEN(I,K,J)=DTDT(K)/pi(I,K,J)
    490. 

  490.                  RQVCUTEN(I,K,J)=DQDT(K)
    491. 

  491.               ENDDO
    492. 

  492.             ENDIF
    493. 

  493. 

  494.             IF(PRESENT(rqrcuten).AND.PRESENT(rqccuten)) THEN
    495. 

  495.               IF( F_QR )THEN
    496. 

  496.                 DO K=kts,kte
    497. 

  497.                    RQRCUTEN(I,K,J)=DQRDT(K)
    498. 

  498.                    RQCCUTEN(I,K,J)=DQCDT(K)
    499. 

  499.                 ENDDO
    500. 

  500.               ELSE
    501. 

  501. ! This is the case for Eta microphysics without 3d rain field
    502. 

  502.                 DO K=kts,kte
    503. 

  503.                    RQRCUTEN(I,K,J)=0.
    504. 

  504.                    RQCCUTEN(I,K,J)=DQRDT(K)+DQCDT(K)
    505. 

  505.                 ENDDO
    506. 

  506.               ENDIF
    507. 

  507.             ENDIF
    508. 

  508. 

  509. !......     QSTEN STORES GRAUPEL TENDENCY IF IT EXISTS, OTHERISE SNOW (V2)
    510. 

  510. 

  511. 

  512.             IF(PRESENT( rqicuten )) THEN
    513. 

  513.               IF ( F_QI ) THEN
    514. 

  514.                 DO K=kts,kte
    515. 

  515.                    RQICUTEN(I,K,J)=DQIDT(K)
    516. 

  516.                 ENDDO
    517. 

  517.               ENDIF
    518. 

  518.             ENDIF
    519. 

  519. 

  520.             IF(PRESENT( rqscuten )) THEN
    521. 

  521.               IF ( F_QS ) THEN
    522. 

  522.                 DO K=kts,kte
    523. 

  523.                    RQSCUTEN(I,K,J)=DQSDT(K)
    524. 

  524.                 ENDDO
    525. 

  525.               ENDIF
    526. 

  526.             ENDIF
    527. 

  527. !
    528. 

  528.          ENDIF 
    529. 

  529.        ENDDO     ! i-loop
    530. 

  530.      ENDDO       ! j-loop
    531. 

  531.    ENDIF         ! run_param
    532. 

  532. !
    533. 

  533.    END SUBROUTINE KF_eta_CPS
    534. 

  534. ! ****************************************************************************
    535. 

  535. !-----------------------------------------------------------
    536. 

  536.    SUBROUTINE KF_eta_PARA (I, J,                           &
    537. 

  537.                       U0,V0,T0,QV0,P0,DZQ,W0AVG1D,         &
    538. 

  538.                       TPART_H0,TPART_V0,                   &
    539. 

  539.                       trigger,                             &
    540. 

  540.                       DT,DX,DXSQ,rhoe,                     &
    541. 

  541.                       XLV0,XLV1,XLS0,XLS1,CP,R,G,          &
    542. 

  542.                       EP2,SVP1,SVP2,SVP3,SVPT0,            &
    543. 

  543.                       DQDT,DQIDT,DQCDT,DQRDT,DQSDT,DTDT,   &
    544. 

  544.                       RAINCV,PRATEC,NCA,                   &
    545. 

  545.                       F_QI,F_QS,warm_rain,                 &
    546. 

  546.                       CUTOP,CUBOT,CUDT,                    &
    547. 

  547.                       ids,ide, jds,jde, kds,kde,           &
    548. 

  548.                       ims,ime, jms,jme, kms,kme,           &
    549. 

  549.                       its,ite, jts,jte, kts,kte)
    550. 

  550. !-----------------------------------------------------------
    551. 

  551. !***** The KF scheme that is currently used in experimental runs of EMCs 
    552. 

  552. !***** Eta model....jsk 8/00
    553. 

  553. !
    554. 

  554.       IMPLICIT NONE
    555. 

  555. !-----------------------------------------------------------
    556. 

  556.       INTEGER, INTENT(IN   ) :: ids,ide, jds,jde, kds,kde, &
    557. 

  557.                                 ims,ime, jms,jme, kms,kme, &
    558. 

  558.                                 its,ite, jts,jte, kts,kte, &
    559. 

  559.                                 I,J
    560. 

  560.           ! ,P_QI,P_QS,P_FIRST_SCALAR
    561. 

  561.       INTEGER, INTENT(IN   ) ::  trigger
    562. 

  562. 

  563.       LOGICAL, INTENT(IN   ) :: F_QI, F_QS
    564. 

  564. 

  565.       LOGICAL, INTENT(IN   ) :: warm_rain
    566. 

  566. !
    567. 

  567.       REAL, DIMENSION( kts:kte ),                          &
    568. 

  568.             INTENT(IN   ) ::                           U0, &
    569. 

  569.                                                        V0, &
    570. 

  570.                                                  TPART_H0, &
    571. 

  571.                                                  TPART_V0, &
    572. 

  572.                                                        T0, &
    573. 

  573.                                                       QV0, &
    574. 

  574.                                                        P0, &
    575. 

  575.                                                      rhoe, &
    576. 

  576.                                                       DZQ, &
    577. 

  577.                                                   W0AVG1D
    578. 

  578. !
    579. 

  579.       REAL,  INTENT(IN   ) :: DT,DX,DXSQ
    580. 

  580. !
    581. 

  581. 

  582.       REAL,  INTENT(IN   ) :: XLV0,XLV1,XLS0,XLS1,CP,R,G
    583. 

  583.       REAL,  INTENT(IN   ) :: EP2,SVP1,SVP2,SVP3,SVPT0
    584. 

  584. 

  585. !
    586. 

  586.       REAL, DIMENSION( kts:kte ), INTENT(INOUT) ::         &
    587. 

  587.                                                      DQDT, &
    588. 

  588.                                                     DQIDT, &
    589. 

  589.                                                     DQCDT, &
    590. 

  590.                                                     DQRDT, &
    591. 

  591.                                                     DQSDT, &
    592. 

  592.                                                      DTDT
    593. 

  593. 

  594.       REAL,    DIMENSION( ims:ime , jms:jme ),             &
    595. 

  595.             INTENT(INOUT) ::                          NCA
    596. 

  596. 

  597.       REAL, DIMENSION( ims:ime , jms:jme ),                &
    598. 

  598.             INTENT(INOUT) ::                       RAINCV
    599. 

  599. 

  600.       REAL, DIMENSION( ims:ime , jms:jme ),                &
    601. 

  601.             INTENT(INOUT) ::                       PRATEC
    602. 

  602. 

  603.       REAL, DIMENSION( ims:ime , jms:jme ),                &
    604. 

  604.             INTENT(OUT) ::                          CUBOT, &
    605. 

  605.                                                     CUTOP
    606. 

  606.       REAL,  INTENT(IN   ) :: CUDT
    607. 

  607. !
    608. 

  608. !...DEFINE LOCAL VARIABLES...
    609. 

  609. !
    610. 

  610.       REAL, DIMENSION( kts:kte ) ::                        &
    611. 

  611.             Q0,Z0,TV0,TU,TVU,QU,TZ,TVD,                    &
    612. 

  612.             QD,QES,THTES,TG,TVG,QG,WU,WD,W0,EMS,EMSD,      &
    613. 

  613.             UMF,UER,UDR,DMF,DER,DDR,UMF2,UER2,             &
    614. 

  614.             UDR2,DMF2,DER2,DDR2,DZA,THTA0,THETEE,          &
    615. 

  615.             THTAU,THETEU,THTAD,THETED,QLIQ,QICE,           &
    616. 

  616.             QLQOUT,QICOUT,PPTLIQ,PPTICE,DETLQ,DETIC,       &
    617. 

  617.             DETLQ2,DETIC2,RATIO,RATIO2
    618. 

  618. 

  619. 

  620.       REAL, DIMENSION( kts:kte ) ::                        &
    621. 

  621.             DOMGDP,EXN,TVQU,DP,RH,EQFRC,WSPD,              &
    622. 

  622.             QDT,FXM,THTAG,THPA,THFXOUT,                    &
    623. 

  623.             THFXIN,QPA,QFXOUT,QFXIN,QLPA,QLFXIN,           &
    624. 

  624.             QLFXOUT,QIPA,QIFXIN,QIFXOUT,QRPA,              &
    625. 

  625.             QRFXIN,QRFXOUT,QSPA,QSFXIN,QSFXOUT,            &
    626. 

  626.             QL0,QLG,QI0,QIG,QR0,QRG,QS0,QSG
    627. 

  627. 

  628. 

  629.       REAL, DIMENSION( kts:kte+1 ) :: OMG
    630. 

  630.       REAL, DIMENSION( kts:kte ) :: RAINFB,SNOWFB
    631. 

  631.       REAL, DIMENSION( kts:kte ) ::                        &
    632. 

  632.             CLDHGT,QSD,DILFRC,DDILFRC,TKE,TGU,QGU,THTEEG
    633. 

  633. 

  634. ! LOCAL VARS
    635. 

  635. 

  636.       REAL    :: P00,T00,RLF,RHIC,RHBC,PIE,         &
    637. 

  637.                  TTFRZ,TBFRZ,C5,RATE
    638. 

  638.       REAL    :: GDRY,ROCP,ALIQ,BLIQ,                      &
    639. 

  639.                  CLIQ,DLIQ
    640. 

  640.       REAL    :: FBFRC,P300,DPTHMX,THMIX,QMIX,ZMIX,PMIX,   &
    641. 

  641.                  ROCPQ,TMIX,EMIX,TLOG,TDPT,TLCL,TVLCL,     &
    642. 

  642.                  CPORQ,PLCL,ES,DLP,TENV,QENV,TVEN,TVBAR,   &
    643. 

  643.                  ZLCL,WKL,WABS,TRPPT,WSIGNE,DTLCL,GDT,WLCL,&
    644. 

  644.                  TVAVG,QESE,WTW,RHOLCL,AU0,VMFLCL,UPOLD,   &
    645. 

  645.                  UPNEW,ABE,WKLCL,TTEMP,FRC1,   &
    646. 

  646.                  QNEWIC,RL,R1,QNWFRZ,EFFQ,BE,BOTERM,ENTERM,&
    647. 

  647.                  DZZ,UDLBE,REI,EE2,UD2,TTMP,F1,F2,         &
    648. 

  648.                  THTTMP,QTMP,TMPLIQ,TMPICE,TU95,TU10,EE1,  &
    649. 

  649.                  UD1,DPTT,QNEWLQ,DUMFDP,EE,TSAT,           &
    650. 

  650.                  THTA,VCONV,TIMEC,SHSIGN,VWS,PEF, &
    651. 

  651.                  CBH,RCBH,PEFCBH,PEFF,PEFF2,TDER,THTMIN,   &
    652. 

  652.                  DTMLTD,QS,TADVEC,DPDD,FRC,DPT,RDD,A1,     &
    653. 

  653.                  DSSDT,DTMP,T1RH,QSRH,PPTFLX,CPR,CNDTNF,   &
    654. 

  654.                  UPDINC,AINCM2,DEVDMF,PPR,RCED,DPPTDF,     &
    655. 

  655.                  DMFLFS,DMFLFS2,RCED2,DDINC,AINCMX,AINCM1, &
    656. 

  656.                  AINC,TDER2,PPTFL2,FABE,STAB,DTT,DTT1,     &
    657. 

  657.                  DTIME,TMA,TMB,TMM,BCOEFF,ACOEFF,QVDIFF,   &
    658. 

  658.                  TOPOMG,CPM,DQ,ABEG,DABE,DFDA,FRC2,DR,     &
    659. 

  659.                  UDFRC,TUC,QGS,RH0,RHG,QINIT,QFNL,ERR2,    &
    660. 

  660.                  RELERR,RLC,RLS,RNC,FABEOLD,AINCOLD,UEFRC, &
    661. 

  661.                  DDFRC,TDC,DEFRC,RHBAR,DMFFRC,DPMIN,DILBE
    662. 

  662.    REAL    ::    ASTRT,TP,VALUE,AINTRP,TKEMAX,QFRZ,&
    663. 

  663.                  QSS,PPTMLT,DTMELT,RHH,EVAC,BINC
    664. 

  664. !
    665. 

  665.       INTEGER :: INDLU,NU,NUCHM,NNN,KLFS
    666. 

  666.    REAL    :: CHMIN,PM15,CHMAX,DTRH,RAD,DPPP
    667. 

  667.    REAL    :: TVDIFF,DTTOT,ABSOMG,ABSOMGTC,FRDP
    668. 

  668. 

  669.       INTEGER :: KX,K,KL
    670. 

  670. !
    671. 

  671.       INTEGER :: NCHECK
    672. 

  672.       INTEGER, DIMENSION (kts:kte) :: KCHECK
    673. 

  673. 

  674.       INTEGER :: ISTOP,ML,L5,KMIX,LOW,                     &
    675. 

  675.                  LC,MXLAYR,LLFC,NLAYRS,NK,                 &
    676. 

  676.                  KPBL,KLCL,LCL,LET,IFLAG,                  &
    677. 

  677.                  NK1,LTOP,NJ,LTOP1,                        &
    678. 

  678.                  LTOPM1,LVF,KSTART,KMIN,LFS,               &
    679. 

  679.                  ND,NIC,LDB,LDT,ND1,NDK,                   &
    680. 

  680.                  NM,LMAX,NCOUNT,NOITR,                     &
    681. 

  681.                  NSTEP,NTC,NCHM,ISHALL,NSHALL
    682. 

  682.       LOGICAL :: IPRNT
    683. 

  683.       REAL :: u00,qslcl,rhlcl,dqssdt    !jfb
    684. 

  684.       CHARACTER*1024 message
    685. 

  685. !
    686. 

  686.       DATA P00,T00/1.E5,273.16/
    687. 

  687.       DATA RLF/3.339E5/
    688. 

  688.       DATA RHIC,RHBC/1.,0.90/
    689. 

  689.       DATA PIE,TTFRZ,TBFRZ,C5/3.141592654,268.16,248.16,1.0723E-3/
    690. 

  690.       DATA RATE/0.03/   ! wrf default
    691. 

  691. !      DATA RATE/0.01/  ! value used in NRCM
    692. 

  692. !      DATA RATE/0.001/  ! effectively turn off autoconversion
    693. 

  693. !-----------------------------------------------------------
    694. 

  694.       IPRNT=.FALSE.
    695. 

  695.       GDRY=-G/CP
    696. 

  696.       ROCP=R/CP
    697. 

  697.       NSHALL = 0
    698. 

  698.       KL=kte
    699. 

  699.       KX=kte
    700. 

  700. !
    701. 

  701. !     ALIQ = 613.3
    702. 

  702. !     BLIQ = 17.502
    703. 

  703. !     CLIQ = 4780.8
    704. 

  704. !     DLIQ = 32.19
    705. 

  705.       ALIQ = SVP1*1000.
    706. 

  706.       BLIQ = SVP2
    707. 

  707.       CLIQ = SVP2*SVPT0
    708. 

  708.       DLIQ = SVP3
    709. 

  709. !
    710. 

  710. !
    711. 

  711. !****************************************************************************
    712. 

  712. !                                                      ! PPT FB MODS
    713. 

  713. !...OPTION TO FEED CONVECTIVELY GENERATED RAINWATER    ! PPT FB MODS
    714. 

  714. !...INTO GRID-RESOLVED RAINWATER (OR SNOW/GRAUPEL)     ! PPT FB MODS
    715. 

  715. !...FIELD.  "FBFRC" IS THE FRACTION OF AVAILABLE       ! PPT FB MODS
    716. 

  716. !...PRECIPITATION TO BE FED BACK (0.0 - 1.0)...        ! PPT FB MODS
    717. 

  717.       FBFRC=0.0                                        ! PPT FB MODS
    718. 

  718. !...mods to allow shallow convection...
    719. 

  719.       NCHM = 0
    720. 

  720.       ISHALL = 0
    721. 

  721.       DPMIN = 5.E3
    722. 

  722. !...
    723. 

  723.       P300=P0(1)-30000.
    724. 

  724. !
    725. 

  725. !...PRESSURE PERTURBATION TERM IS ONLY DEFINED AT MID-POINT OF
    726. 

  726. !...VERTICAL LAYERS...SINCE TOTAL PRESSURE IS NEEDED AT THE TOP AND
    727. 

  727. !...BOTTOM OF LAYERS BELOW, DO AN INTERPOLATION...
    728. 

  728. !
    729. 

  729. !...INPUT A VERTICAL SOUNDING ... NOTE THAT MODEL LAYERS ARE NUMBERED
    730. 

  730. !...FROM BOTTOM-UP IN THE KF SCHEME...
    731. 

  731. !
    732. 

  732.       ML=0 
    733. 

  733. !SUE  tmprpsb=1./PSB(I,J)
    734. 

  734. !SUE  CELL=PTOP*tmprpsb
    735. 

  735. !
    736. 

  736.       DO K=1,KX
    737. 

  737. !
    738. 

  738. !  Saturation vapor pressure (ES) is calculated following Buck (1981)
    739. 

  739. !...IF Q0 IS ABOVE SATURATION VALUE, REDUCE IT TO SATURATION LEVEL...
    740. 

  740. !
    741. 

  741.          ES=ALIQ*EXP((BLIQ*T0(K)-CLIQ)/(T0(K)-DLIQ))
    742. 

  742.          QES(K)=0.622*ES/(P0(K)-ES)
    743. 

  743.          Q0(K)=AMIN1(QES(K),QV0(K))
    744. 

  744.          Q0(K)=AMAX1(0.000001,Q0(K))
    745. 

  745.          QL0(K)=0.
    746. 

  746.          QI0(K)=0.
    747. 

  747.          QR0(K)=0.
    748. 

  748.          QS0(K)=0.
    749. 

  749.          RH(K) = Q0(K)/QES(K)
    750. 

  750.          DILFRC(K) = 1.
    751. 

  751.          TV0(K)=T0(K)*(1.+0.608*Q0(K))
    752. 

  752. !        RHOE(K)=P0(K)/(R*TV0(K))
    753. 

  753. !   DP IS THE PRESSURE INTERVAL BETWEEN FULL SIGMA LEVELS...
    754. 

  754.          DP(K)=rhoe(k)*g*DZQ(k)
    755. 

  755. ! IF Turbulent Kinetic Energy (TKE) is available from turbulent mixing scheme
    756. 

  756. ! use it for shallow convection...For now, assume it is not available....
    757. 

  757. !         TKE(K) = Q2(I,J,NK)
    758. 

  758.          TKE(K) = 0.
    759. 

  759.          CLDHGT(K) = 0.
    760. 

  760. !        IF(P0(K).GE.500E2)L5=K
    761. 

  761.          IF(P0(K).GE.0.5*P0(1))L5=K
    762. 

  762.          IF(P0(K).GE.P300)LLFC=K
    763. 

  763.       ENDDO
    764. 

  764. !
    765. 

  765. !...DZQ IS DZ BETWEEN SIGMA SURFACES, DZA IS DZ BETWEEN MODEL HALF LEVEL
    766. 

  766.         Z0(1)=.5*DZQ(1)
    767. 

  767. !cdir novector
    768. 

  768.         DO K=2,KL
    769. 

  769.           Z0(K)=Z0(K-1)+.5*(DZQ(K)+DZQ(K-1))
    770. 

  770.           DZA(K-1)=Z0(K)-Z0(K-1)
    771. 

  771.         ENDDO   
    772. 

  772.         DZA(KL)=0.
    773. 

  773. !
    774. 

  774. !
    775. 

  775. !  To save time, specify a pressure interval to move up in sequential
    776. 

  776. !  check of different ~50 mb deep groups of adjacent model layers in
    777. 

  777. !  the process of identifying updraft source layer (USL).  Note that 
    778. 

  778. !  this search is terminated as soon as a buoyant parcel is found and 
    779. 

  779. !  this parcel can produce a cloud greater than specifed minimum depth
    780. 

  780. !  (CHMIN)...For now, set interval at 15 mb...
    781. 

  781. !
    782. 

  782.        NCHECK = 1
    783. 

  783.        KCHECK(NCHECK)=1
    784. 

  784.        PM15 = P0(1)-15.E2
    785. 

  785.        DO K=2,LLFC
    786. 

  786.          IF(P0(K).LT.PM15)THEN
    787. 

  787.            NCHECK = NCHECK+1
    788. 

  788.            KCHECK(NCHECK) = K
    789. 

  789.            PM15 = PM15-15.E2
    790. 

  790.          ENDIF
    791. 

  791.        ENDDO
    792. 

  792. !
    793. 

  793.        NU=0
    794. 

  794.        NUCHM=0
    795. 

  795. usl:   DO
    796. 

  796.            NU = NU+1
    797. 

  797.            IF(NU.GT.NCHECK)THEN 
    798. 

  798.              IF(ISHALL.EQ.1)THEN
    799. 

  799.                CHMAX = 0.
    800. 

  800.                NCHM = 0
    801. 

  801.                DO NK = 1,NCHECK
    802. 

  802.                  NNN=KCHECK(NK)
    803. 

  803.                  IF(CLDHGT(NNN).GT.CHMAX)THEN
    804. 

  804.                    NCHM = NNN
    805. 

  805.                    NUCHM = NK
    806. 

  806.                    CHMAX = CLDHGT(NNN)
    807. 

  807.                  ENDIF
    808. 

  808.                ENDDO
    809. 

  809.                NU = NUCHM-1
    810. 

  810.                FBFRC=1.
    811. 

  811.                CYCLE usl
    812. 

  812.              ELSE
    813. 

  813.                RETURN
    814. 

  814.              ENDIF
    815. 

  815.            ENDIF      
    816. 

  816.            KMIX = KCHECK(NU)
    817. 

  817.            LOW=KMIX
    818. 

  818. !...
    819. 

  819.            LC = LOW
    820. 

  820. !
    821. 

  821. !...ASSUME THAT IN ORDER TO SUPPORT A DEEP UPDRAFT YOU NEED A LAYER OF
    822. 

  822. !...UNSTABLE AIR AT LEAST 50 mb DEEP...TO APPROXIMATE THIS, ISOLATE A
    823. 

  823. !...GROUP OF ADJACENT INDIVIDUAL MODEL LAYERS, WITH THE BASE AT LEVEL
    824. 

  824. !...LC, SUCH THAT THE COMBINED DEPTH OF THESE LAYERS IS AT LEAST 50 mb..
    825. 

  825. !   
    826. 

  826.            NLAYRS=0
    827. 

  827.            DPTHMX=0.
    828. 

  828.            NK=LC-1
    829. 

  829.            IF ( NK+1 .LT. KTS ) THEN
    830. 

  830.              WRITE(message,*)'WOULD GO OFF BOTTOM: KF_ETA_PARA I,J,NK',I,J,NK
    831. 

  831.              CALL wrf_message (TRIM(message)) 
    832. 

  832.            ELSE
    833. 

  833.              DO 
    834. 

  834.                NK=NK+1   
    835. 

  835.                IF ( NK .GT. KTE ) THEN
    836. 

  836.                  WRITE(message,*)'WOULD GO OFF TOP: KF_ETA_PARA I,J,DPTHMX,DPMIN',I,J,DPTHMX,DPMIN
    837. 

  837.                  CALL wrf_message (TRIM(message))
    838. 

  838.                  EXIT
    839. 

  839.                ENDIF
    840. 

  840.                DPTHMX=DPTHMX+DP(NK)
    841. 

  841.                NLAYRS=NLAYRS+1
    842. 

  842.                IF(DPTHMX.GT.DPMIN)THEN
    843. 

  843.                  EXIT 
    844. 

  844.                ENDIF
    845. 

  845.              END DO    
    846. 

  846.            ENDIF
    847. 

  847.            IF(DPTHMX.LT.DPMIN)THEN 
    848. 

  848.              RETURN
    849. 

  849.            ENDIF
    850. 

  850.            KPBL=LC+NLAYRS-1   
    851. 

  851. !
    852. 

  852. !...********************************************************
    853. 

  853. !...for computational simplicity without much loss in accuracy,
    854. 

  854. !...mix temperature instead of theta for evaluating convective
    855. 

  855. !...initiation (triggering) potential...
    856. 

  856. !          THMIX=0.
    857. 

  857.            TMIX=0.
    858. 

  858.            QMIX=0.
    859. 

  859.            ZMIX=0.
    860. 

  860.            PMIX=0.
    861. 

  861. !
    862. 

  862. !...FIND THE THERMODYNAMIC CHARACTERISTICS OF THE LAYER BY
    863. 

  863. !...MASS-WEIGHTING THE CHARACTERISTICS OF THE INDIVIDUAL MODEL
    864. 

  864. !...LAYERS...
    865. 

  865. !
    866. 

  866. !cdir novector
    867. 

  867.            DO NK=LC,KPBL
    868. 

  868.              TMIX=TMIX+DP(NK)*T0(NK)
    869. 

  869.              QMIX=QMIX+DP(NK)*Q0(NK)
    870. 

  870.              ZMIX=ZMIX+DP(NK)*Z0(NK)
    871. 

  871.              PMIX=PMIX+DP(NK)*P0(NK)
    872. 

  872.            ENDDO   
    873. 

  873. !         THMIX=THMIX/DPTHMX
    874. 

  874.           TMIX=TMIX/DPTHMX
    875. 

  875.           QMIX=QMIX/DPTHMX
    876. 

  876.           ZMIX=ZMIX/DPTHMX
    877. 

  877.           PMIX=PMIX/DPTHMX
    878. 

  878.           EMIX=QMIX*PMIX/(0.622+QMIX)
    879. 

  879. !
    880. 

  880. !...FIND THE TEMPERATURE OF THE MIXTURE AT ITS LCL...
    881. 

  881. !
    882. 

  882. !        TLOG=ALOG(EMIX/ALIQ)
    883. 

  883. ! ...calculate dewpoint using lookup table...
    884. 

  884. !
    885. 

  885.           astrt=1.e-3
    886. 

  886.           ainc=0.075
    887. 

  887.           a1=emix/aliq
    888. 

  888.           tp=(a1-astrt)/ainc
    889. 

  889.           indlu=int(tp)+1
    890. 

  890.           value=(indlu-1)*ainc+astrt
    891. 

  891.           aintrp=(a1-value)/ainc
    892. 

  892.           tlog=aintrp*alu(indlu+1)+(1-aintrp)*alu(indlu)
    893. 

  893.           TDPT=(CLIQ-DLIQ*TLOG)/(BLIQ-TLOG)
    894. 

  894.           TLCL=TDPT-(.212+1.571E-3*(TDPT-T00)-4.36E-4*(TMIX-T00))*(TMIX-TDPT)
    895. 

  895.           TLCL=AMIN1(TLCL,TMIX)
    896. 

  896.           TVLCL=TLCL*(1.+0.608*QMIX)
    897. 

  897.           ZLCL = ZMIX+(TLCL-TMIX)/GDRY
    898. 

  898.           NK = LC-1
    899. 

  899.           DO 
    900. 

  900.             NK = NK+1
    901. 

  901.             KLCL=NK
    902. 

  902.             IF(ZLCL.LE.Z0(NK) .or. NK.GT.KL)THEN
    903. 

  903.               EXIT
    904. 

  904.             ENDIF 
    905. 

  905.           ENDDO   
    906. 

  906.           IF(NK.GT.KL)THEN
    907. 

  907.             RETURN  
    908. 

  908.           ENDIF
    909. 

  909.           K=KLCL-1
    910. 

  910. ! calculate DLP using Z instead of log(P)
    911. 

  911.           DLP=(ZLCL-Z0(K))/(Z0(KLCL)-Z0(K))
    912. 

  912. !     
    913. 

  913. !...ESTIMATE ENVIRONMENTAL TEMPERATURE AND MIXING RATIO AT THE LCL...
    914. 

  914. !     
    915. 

  915.           TENV=T0(K)+(T0(KLCL)-T0(K))*DLP
    916. 

  916.           QENV=Q0(K)+(Q0(KLCL)-Q0(K))*DLP
    917. 

  917.           TVEN=TENV*(1.+0.608*QENV)
    918. 

  918. !     
    919. 

  919. !...CHECK TO SEE IF CLOUD IS BUOYANT USING FRITSCH-CHAPPELL TRIGGER
    920. 

  920. !...FUNCTION DESCRIBED IN KAIN AND FRITSCH (1992)...W0 IS AN
    921. 

  921. !...APROXIMATE VALUE FOR THE RUNNING-MEAN GRID-SCALE VERTICAL
    922. 

  922. !...VELOCITY, WHICH GIVES SMOOTHER FIELDS OF CONVECTIVE INITIATION
    923. 

  923. !...THAN THE INSTANTANEOUS VALUE...FORMULA RELATING TEMPERATURE
    924. 

  924. !...PERTURBATION TO VERTICAL VELOCITY HAS BEEN USED WITH THE MOST
    925. 

  925. !...SUCCESS AT GRID LENGTHS NEAR 25 km.  FOR DIFFERENT GRID-LENGTHS,
    926. 

  926. !...ADJUST VERTICAL VELOCITY TO EQUIVALENT VALUE FOR 25 KM GRID
    927. 

  927. !...LENGTH, ASSUMING LINEAR DEPENDENCE OF W ON GRID LENGTH...
    928. 

  928. !     
    929. 

  929.           IF(ZLCL.LT.2.E3)THEN        ! Kain (2004) Eq. 2
    930. 

  930.             WKLCL=0.02*ZLCL/2.E3
    931. 

  931.           ELSE
    932. 

  932.             WKLCL=0.02                ! units of m/s
    933. 

  933.           ENDIF
    934. 

  934.           WKL=(W0AVG1D(K)+(W0AVG1D(KLCL)-W0AVG1D(K))*DLP)*DX/25.E3-WKLCL
    935. 

  935.           IF(WKL.LT.0.0001)THEN
    936. 

  936.             DTLCL=0.
    937. 

  937.           ELSE 
    938. 

  938.             DTLCL=4.64*WKL**0.33      ! Kain (2004) Eq. 1
    939. 

  939.           ENDIF
    940. 

  940. 

  941. ! New trigger function
    942. 

  942.            IF(trigger.eq.2) then
    943. 

  943.               DTLCL=amax1(TPART_H0(KLCL)+TPART_V0(KLCL),0.0)
    944. 

  944.            ENDIF
    945. 

  945. ! end for trigger function
    946. 

  946. !
    947. 

  947.         dtrh = 0.
    948. 

  948. 	if (trigger  .eq. 3) then
    949. 

  949. !...for ETA model, give parcel an extra temperature perturbation based
    950. 

  950. !...the threshold RH for condensation (U00)...
    951. 

  951. ! as described in Narita and Ohmori (2007, 12th Mesoscale Conf.
    952. 

  952. !
    953. 

  953. !...for now, just assume U00=0.75...
    954. 

  954. !...!!!!!! for MM5, SET DTRH = 0. !!!!!!!!
    955. 

  955.           U00 = 0.75
    956. 

  956.           IF(U00.lt.1.)THEN
    957. 

  957.             QSLCL=QES(K)+(QES(KLCL)-QES(K))*DLP
    958. 

  958.             RHLCL = QENV/QSLCL
    959. 

  959.             DQSSDT = QMIX*(CLIQ-BLIQ*DLIQ)/((TLCL-DLIQ)*(TLCL-DLIQ))
    960. 

  960.             IF(RHLCL.ge.0.75 .and. RHLCL.le.0.95)then
    961. 

  961.               DTRH = 0.25*(RHLCL-0.75)*QMIX/DQSSDT
    962. 

  962.             ELSEIF(RHLCL.GT.0.95)THEN
    963. 

  963.               DTRH = (1./RHLCL-1.)*QMIX/DQSSDT
    964. 

  964.             ELSE
    965. 

  965.                DTRH = 0.
    966. 

  966.             ENDIF
    967. 

  967.           ENDIF   
    968. 

  968.         endif   ! trigger 3
    969. 

  969. !         IF(ISHALL.EQ.1)IPRNT=.TRUE.
    970. 

  970. !         IPRNT=.TRUE.
    971. 

  971. !         IF(TLCL+DTLCL.GT.TENV)GOTO 45
    972. 

  972.  
    973. 

  973. trigger2:  IF(TLCL+DTLCL+DTRH.LT.TENV)THEN   
    974. 

  974. !
    975. 

  975. ! Parcel not buoyant, CYCLE back to start of trigger and evaluate next potential USL...
    976. 

  976. !
    977. 

  977.             CYCLE usl
    978. 

  978. !
    979. 

  979.           ELSE                            ! Parcel is buoyant, determine updraft
    980. 

  980. !     
    981. 

  981. !...CONVECTIVE TRIGGERING CRITERIA HAS BEEN SATISFIED...COMPUTE
    982. 

  982. !...EQUIVALENT POTENTIAL TEMPERATURE
    983. 

  983. !...(THETEU) AND VERTICAL VELOCITY OF THE RISING PARCEL AT THE LCL...
    984. 

  984. !     
    985. 

  985.             CALL ENVIRTHT(PMIX,TMIX,QMIX,THETEU(K),ALIQ,BLIQ,CLIQ,DLIQ)
    986. 

  986. !
    987. 

  987. !...modify calculation of initial parcel vertical velocity...jsk 11/26/97
    988. 

  988. !
    989. 

  989.             DTTOT = DTLCL+DTRH
    990. 

  990.             IF(DTTOT.GT.1.E-4)THEN
    991. 

  991.               GDT=2.*G*DTTOT*500./TVEN     ! Kain (2004) Eq. 3  (sort of)
    992. 

  992.               WLCL=1.+0.5*SQRT(GDT)
    993. 

  993.               WLCL = AMIN1(WLCL,3.)
    994. 

  994.             ELSE
    995. 

  995.               WLCL=1.
    996. 

  996.             ENDIF
    997. 

  997.             PLCL=P0(K)+(P0(KLCL)-P0(K))*DLP
    998. 

  998.             WTW=WLCL*WLCL
    999. 

  999. !
    1000. 

  1000.             TVLCL=TLCL*(1.+0.608*QMIX)
    1001. 

  1001.             RHOLCL=PLCL/(R*TVLCL)
    1002. 

  1002. !        
    1003. 

  1003.             LCL=KLCL
    1004. 

  1004.             LET=LCL
    1005. 

  1005. ! make RAD a function of background vertical velocity...    (Kain (2004) Eq. 6)
    1006. 

  1006.             IF(WKL.LT.0.)THEN
    1007. 

  1007.               RAD = 1000.
    1008. 

  1008.             ELSEIF(WKL.GT.0.1)THEN
    1009. 

  1009.               RAD = 2000.
    1010. 

  1010.             ELSE
    1011. 

  1011.               RAD = 1000.+1000*WKL/0.1
    1012. 

  1012.             ENDIF
    1013. 

  1013. !     
    1014. 

  1014. !*******************************************************************
    1015. 

  1015. !                                                                  *
    1016. 

  1016. !                 COMPUTE UPDRAFT PROPERTIES                       *
    1017. 

  1017. !                                                                  *
    1018. 

  1018. !*******************************************************************
    1019. 

  1019. !     
    1020. 

  1020. !     
    1021. 

  1021. !...
    1022. 

  1022. !...ESTIMATE INITIAL UPDRAFT MASS FLUX (UMF(K))...
    1023. 

  1023. !     
    1024. 

  1024.             WU(K)=WLCL
    1025. 

  1025.             AU0=0.01*DXSQ
    1026. 

  1026.             UMF(K)=RHOLCL*AU0
    1027. 

  1027.             VMFLCL=UMF(K)
    1028. 

  1028.             UPOLD=VMFLCL
    1029. 

  1029.             UPNEW=UPOLD
    1030. 

  1030. !     
    1031. 

  1031. !...RATIO2 IS THE DEGREE OF GLACIATION IN THE CLOUD (0 TO 1),
    1032. 

  1032. !...UER IS THE ENVIR ENTRAINMENT RATE, ABE IS AVAILABLE
    1033. 

  1033. !...BUOYANT ENERGY, TRPPT IS THE TOTAL RATE OF PRECIPITATION
    1034. 

  1034. !...PRODUCTION...
    1035. 

  1035. !     
    1036. 

  1036.             RATIO2(K)=0.
    1037. 

  1037.             UER(K)=0.
    1038. 

  1038.             ABE=0.
    1039. 

  1039.             TRPPT=0.
    1040. 

  1040.             TU(K)=TLCL
    1041. 

  1041.             TVU(K)=TVLCL
    1042. 

  1042.             QU(K)=QMIX
    1043. 

  1043.             EQFRC(K)=1.
    1044. 

  1044.             QLIQ(K)=0.
    1045. 

  1045.             QICE(K)=0.
    1046. 

  1046.             QLQOUT(K)=0.
    1047. 

  1047.             QICOUT(K)=0.
    1048. 

  1048.             DETLQ(K)=0.
    1049. 

  1049.             DETIC(K)=0.
    1050. 

  1050.             PPTLIQ(K)=0.
    1051. 

  1051.             PPTICE(K)=0.
    1052. 

  1052.             IFLAG=0
    1053. 

  1053. !     
    1054. 

  1054. !...TTEMP IS USED DURING CALCULATION OF THE LINEAR GLACIATION
    1055. 

  1055. !...PROCESS; IT IS INITIALLY SET TO THE TEMPERATURE AT WHICH
    1056. 

  1056. !...FREEZING IS SPECIFIED TO BEGIN.  WITHIN THE GLACIATION
    1057. 

  1057. !...INTERVAL, IT IS SET EQUAL TO THE UPDRAFT TEMP AT THE
    1058. 

  1058. !...PREVIOUS MODEL LEVEL...
    1059. 

  1059. !     
    1060. 

  1060.             TTEMP=TTFRZ
    1061. 

  1061. !     
    1062. 

  1062. !...ENTER THE LOOP FOR UPDRAFT CALCULATIONS...CALCULATE UPDRAFT TEMP,
    1063. 

  1063. !...MIXING RATIO, VERTICAL MASS FLUX, LATERAL DETRAINMENT OF MASS AND
    1064. 

  1064. !...MOISTURE, PRECIPITATION RATES AT EACH MODEL LEVEL...
    1065. 

  1065. !     
    1066. 

  1066. !    **1 variables indicate the bottom of a model layer
    1067. 

  1067. !    **2 variables indicate the top of a model layer
    1068. 

  1068. !     
    1069. 

  1069.             EE1=1.
    1070. 

  1070.             UD1=0.
    1071. 

  1071.             REI = 0.
    1072. 

  1072.             DILBE = 0.
    1073. 

  1073. updraft:    DO NK=K,KL-1
    1074. 

  1074.               NK1=NK+1
    1075. 

  1075.               RATIO2(NK1)=RATIO2(NK)
    1076. 

  1076.               FRC1=0.
    1077. 

  1077.               TU(NK1)=T0(NK1)
    1078. 

  1078.               THETEU(NK1)=THETEU(NK)
    1079. 

  1079.               QU(NK1)=QU(NK)
    1080. 

  1080.               QLIQ(NK1)=QLIQ(NK)
    1081. 

  1081.               QICE(NK1)=QICE(NK)
    1082. 

  1082.               call tpmix2(p0(nk1),theteu(nk1),tu(nk1),qu(nk1),qliq(nk1),        &
    1083. 

  1083.                      qice(nk1),qnewlq,qnewic,XLV1,XLV0)
    1084. 

  1084. !
    1085. 

  1085. !
    1086. 

  1086. !...CHECK TO SEE IF UPDRAFT TEMP IS ABOVE THE TEMPERATURE AT WHICH
    1087. 

  1087. !...GLACIATION IS ASSUMED TO INITIATE; IF IT IS, CALCULATE THE
    1088. 

  1088. !...FRACTION OF REMAINING LIQUID WATER TO FREEZE...TTFRZ IS THE
    1089. 

  1089. !...TEMP AT WHICH FREEZING BEGINS, TBFRZ THE TEMP BELOW WHICH ALL
    1090. 

  1090. !...LIQUID WATER IS FROZEN AT EACH LEVEL...
    1091. 

  1091. !
    1092. 

  1092.               IF(TU(NK1).LE.TTFRZ)THEN
    1093. 

  1093.                 IF(TU(NK1).GT.TBFRZ)THEN
    1094. 

  1094.                   IF(TTEMP.GT.TTFRZ)TTEMP=TTFRZ
    1095. 

  1095.                   FRC1=(TTEMP-TU(NK1))/(TTEMP-TBFRZ)
    1096. 

  1096.                 ELSE
    1097. 

  1097.                   FRC1=1.
    1098. 

  1098.                   IFLAG=1
    1099. 

  1099.                 ENDIF
    1100. 

  1100.                 TTEMP=TU(NK1)
    1101. 

  1101. !
    1102. 

  1102. !  DETERMINE THE EFFECTS OF LIQUID WATER FREEZING WHEN TEMPERATURE
    1103. 

  1103. !...IS BELOW TTFRZ...
    1104. 

  1104. !
    1105. 

  1105.                 QFRZ = (QLIQ(NK1)+QNEWLQ)*FRC1
    1106. 

  1106.                 QNEWIC=QNEWIC+QNEWLQ*FRC1
    1107. 

  1107.                 QNEWLQ=QNEWLQ-QNEWLQ*FRC1
    1108. 

  1108.                 QICE(NK1) = QICE(NK1)+QLIQ(NK1)*FRC1
    1109. 

  1109.                 QLIQ(NK1) = QLIQ(NK1)-QLIQ(NK1)*FRC1
    1110. 

  1110.                 CALL DTFRZNEW(TU(NK1),P0(NK1),THETEU(NK1),QU(NK1),QFRZ,         &
    1111. 

  1111.                           QICE(NK1),ALIQ,BLIQ,CLIQ,DLIQ)
    1112. 

  1112.               ENDIF
    1113. 

  1113.               TVU(NK1)=TU(NK1)*(1.+0.608*QU(NK1))
    1114. 

  1114. !
    1115. 

  1115. !  CALCULATE UPDRAFT VERTICAL VELOCITY AND PRECIPITATION FALLOUT...
    1116. 

  1116. !
    1117. 

  1117.               IF(NK.EQ.K)THEN
    1118. 

  1118.                 BE=(TVLCL+TVU(NK1))/(TVEN+TV0(NK1))-1.
    1119. 

  1119.                 BOTERM=2.*(Z0(NK1)-ZLCL)*G*BE/1.5
    1120. 

  1120.                 DZZ=Z0(NK1)-ZLCL
    1121. 

  1121.               ELSE
    1122. 

  1122.                 BE=(TVU(NK)+TVU(NK1))/(TV0(NK)+TV0(NK1))-1.
    1123. 

  1123.                 BOTERM=2.*DZA(NK)*G*BE/1.5
    1124. 

  1124.                 DZZ=DZA(NK)
    1125. 

  1125.               ENDIF
    1126. 

  1126.               ENTERM=2.*REI*WTW/UPOLD
    1127. 

  1127. 

  1128.               CALL CONDLOAD(QLIQ(NK1),QICE(NK1),WTW,DZZ,BOTERM,ENTERM,      &
    1129. 

  1129.                         RATE,QNEWLQ,QNEWIC,QLQOUT(NK1),QICOUT(NK1),G)
    1130. 

  1130. !
    1131. 

  1131. !...IF VERT VELOCITY IS LESS THAN ZERO, EXIT THE UPDRAFT LOOP AND,
    1132. 

  1132. !...IF CLOUD IS TALL ENOUGH, FINALIZE UPDRAFT CALCULATIONS...
    1133. 

  1133. !
    1134. 

  1134.               IF(WTW.LT.1.E-3)THEN
    1135. 

  1135.                 EXIT
    1136. 

  1136.               ELSE
    1137. 

  1137.                 WU(NK1)=SQRT(WTW)
    1138. 

  1138.               ENDIF
    1139. 

  1139. !...Calculate value of THETA-E in environment to entrain into updraft...
    1140. 

  1140. !
    1141. 

  1141.               CALL ENVIRTHT(P0(NK1),T0(NK1),Q0(NK1),THETEE(NK1),ALIQ,BLIQ,CLIQ,DLIQ)
    1142. 

  1142. !
    1143. 

  1143. !...REI IS THE RATE OF ENVIRONMENTAL INFLOW...
    1144. 

  1144. !
    1145. 

  1145.               REI=VMFLCL*DP(NK1)*0.03/RAD          !    KF (1990) Eq. 1; Kain (2004) Eq. 5
    1146. 

  1146.               TVQU(NK1)=TU(NK1)*(1.+0.608*QU(NK1)-QLIQ(NK1)-QICE(NK1))
    1147. 

  1147.               IF(NK.EQ.K)THEN
    1148. 

  1148.                 DILBE=((TVLCL+TVQU(NK1))/(TVEN+TV0(NK1))-1.)*DZZ
    1149. 

  1149.               ELSE
    1150. 

  1150.                 DILBE=((TVQU(NK)+TVQU(NK1))/(TV0(NK)+TV0(NK1))-1.)*DZZ
    1151. 

  1151.               ENDIF
    1152. 

  1152.               IF(DILBE.GT.0.)ABE=ABE+DILBE*G
    1153. 

  1153. !
    1154. 

  1154. !...IF CLOUD PARCELS ARE VIRTUALLY COLDER THAN THE ENVIRONMENT, MINIMAL 
    1155. 

  1155. !...ENTRAINMENT (0.5*REI) IS IMPOSED...
    1156. 

  1156. !
    1157. 

  1157.               IF(TVQU(NK1).LE.TV0(NK1))THEN    ! Entrain/Detrain IF BLOCK
    1158. 

  1158.                 EE2=0.5       ! Kain (2004)  Eq. 4
    1159. 

  1159.                 UD2=1.
    1160. 

  1160.                 EQFRC(NK1)=0.
    1161. 

  1161.               ELSE
    1162. 

  1162.                 LET=NK1
    1163. 

  1163.                 TTMP=TVQU(NK1)
    1164. 

  1164. !
    1165. 

  1165. !...DETERMINE THE CRITICAL MIXED FRACTION OF UPDRAFT AND ENVIRONMENTAL AIR...
    1166. 

  1166. !
    1167. 

  1167.                 F1=0.95
    1168. 

  1168.                 F2=1.-F1
    1169. 

  1169.                 THTTMP=F1*THETEE(NK1)+F2*THETEU(NK1)
    1170. 

  1170.                 QTMP=F1*Q0(NK1)+F2*QU(NK1)
    1171. 

  1171.                 TMPLIQ=F2*QLIQ(NK1)
    1172. 

  1172.                 TMPICE=F2*QICE(NK1)
    1173. 

  1173.                 call tpmix2(p0(nk1),thttmp,ttmp,qtmp,tmpliq,tmpice,        &
    1174. 

  1174.                            qnewlq,qnewic,XLV1,XLV0)
    1175. 

  1175.                 TU95=TTMP*(1.+0.608*QTMP-TMPLIQ-TMPICE)
    1176. 

  1176.                 IF(TU95.GT.TV0(NK1))THEN
    1177. 

  1177.                   EE2=1.
    1178. 

  1178.                   UD2=0.
    1179. 

  1179.                   EQFRC(NK1)=1.0
    1180. 

  1180.                 ELSE
    1181. 

  1181.                   F1=0.10
    1182. 

  1182.                   F2=1.-F1
    1183. 

  1183.                   THTTMP=F1*THETEE(NK1)+F2*THETEU(NK1)
    1184. 

  1184.                   QTMP=F1*Q0(NK1)+F2*QU(NK1)
    1185. 

  1185.                   TMPLIQ=F2*QLIQ(NK1)
    1186. 

  1186.                   TMPICE=F2*QICE(NK1)
    1187. 

  1187.                   call tpmix2(p0(nk1),thttmp,ttmp,qtmp,tmpliq,tmpice,        &
    1188. 

  1188.                                qnewlq,qnewic,XLV1,XLV0)
    1189. 

  1189.                   TU10=TTMP*(1.+0.608*QTMP-TMPLIQ-TMPICE)
    1190. 

  1190.                   TVDIFF = ABS(TU10-TVQU(NK1))
    1191. 

  1191.                   IF(TVDIFF.LT.1.e-3)THEN
    1192. 

  1192.                     EE2=1.
    1193. 

  1193.                     UD2=0.
    1194. 

  1194.                     EQFRC(NK1)=1.0
    1195. 

  1195.                   ELSE
    1196. 

  1196.                     EQFRC(NK1)=(TV0(NK1)-TVQU(NK1))*F1/(TU10-TVQU(NK1))
    1197. 

  1197.                     EQFRC(NK1)=AMAX1(0.0,EQFRC(NK1))
    1198. 

  1198.                     EQFRC(NK1)=AMIN1(1.0,EQFRC(NK1))
    1199. 

  1199.                     IF(EQFRC(NK1).EQ.1)THEN
    1200. 

  1200.                       EE2=1.
    1201. 

  1201.                       UD2=0.
    1202. 

  1202.                     ELSEIF(EQFRC(NK1).EQ.0.)THEN
    1203. 

  1203.                       EE2=0.
    1204. 

  1204.                       UD2=1.
    1205. 

  1205.                     ELSE
    1206. 

  1206. !
    1207. 

  1207. !...SUBROUTINE PROF5 INTEGRATES OVER THE GAUSSIAN DIST TO DETERMINE THE
    1208. 

  1208. !   FRACTIONAL ENTRAINMENT AND DETRAINMENT RATES...
    1209. 

  1209. !
    1210. 

  1210.                       CALL PROF5(EQFRC(NK1),EE2,UD2)
    1211. 

  1211.                     ENDIF
    1212. 

  1212.                   ENDIF
    1213. 

  1213.                 ENDIF
    1214. 

  1214.               ENDIF                            ! End of Entrain/Detrain IF BLOCK
    1215. 

  1215. !
    1216. 

  1216. !
    1217. 

  1217. !...NET ENTRAINMENT AND DETRAINMENT RATES ARE GIVEN BY THE AVERAGE FRACTIONAL
    1218. 

  1218. !   VALUES IN THE LAYER...
    1219. 

  1219. !
    1220. 

  1220.               EE2 = AMAX1(EE2,0.5)
    1221. 

  1221.               UD2 = 1.5*UD2
    1222. 

  1222.               UER(NK1)=0.5*REI*(EE1+EE2)
    1223. 

  1223.               UDR(NK1)=0.5*REI*(UD1+UD2)
    1224. 

  1224. !
    1225. 

  1225. !...IF THE CALCULATED UPDRAFT DETRAINMENT RATE IS GREATER THAN THE TOTAL
    1226. 

  1226. !   UPDRAFT MASS FLUX, ALL CLOUD MASS DETRAINS, EXIT UPDRAFT CALCULATIONS...
    1227. 

  1227. !
    1228. 

  1228.               IF(UMF(NK)-UDR(NK1).LT.10.)THEN
    1229. 

  1229. !
    1230. 

  1230. !...IF THE CALCULATED DETRAINED MASS FLUX IS GREATER THAN THE TOTAL UPD MASS
    1231. 

  1231. !   FLUX, IMPOSE TOTAL DETRAINMENT OF UPDRAFT MASS AT THE PREVIOUS MODEL LVL..
    1232. 

  1232. !   First, correct ABE calculation if needed...
    1233. 

  1233. !
    1234. 

  1234.                 IF(DILBE.GT.0.)THEN
    1235. 

  1235.                   ABE=ABE-DILBE*G
    1236. 

  1236.                 ENDIF
    1237. 

  1237.                 LET=NK
    1238. 

  1238. !               WRITE(98,1015)P0(NK1)/100.
    1239. 

  1239.                 EXIT 
    1240. 

  1240.               ELSE
    1241. 

  1241.                 EE1=EE2
    1242. 

  1242.                 UD1=UD2
    1243. 

  1243.                 UPOLD=UMF(NK)-UDR(NK1)
    1244. 

  1244.                 UPNEW=UPOLD+UER(NK1)
    1245. 

  1245.                 UMF(NK1)=UPNEW
    1246. 

  1246.                 DILFRC(NK1) = UPNEW/UPOLD
    1247. 

  1247. !
    1248. 

  1248. !...DETLQ AND DETIC ARE THE RATES OF DETRAINMENT OF LIQUID AND
    1249. 

  1249. !...ICE IN THE DETRAINING UPDRAFT MASS...
    1250. 

  1250. !
    1251. 

  1251.                 DETLQ(NK1)=QLIQ(NK1)*UDR(NK1)
    1252. 

  1252.                 DETIC(NK1)=QICE(NK1)*UDR(NK1)
    1253. 

  1253.                 QDT(NK1)=QU(NK1)
    1254. 

  1254.                 QU(NK1)=(UPOLD*QU(NK1)+UER(NK1)*Q0(NK1))/UPNEW
    1255. 

  1255.                 THETEU(NK1)=(THETEU(NK1)*UPOLD+THETEE(NK1)*UER(NK1))/UPNEW
    1256. 

  1256.                 QLIQ(NK1)=QLIQ(NK1)*UPOLD/UPNEW
    1257. 

  1257.                 QICE(NK1)=QICE(NK1)*UPOLD/UPNEW
    1258. 

  1258. !
    1259. 

  1259. !...PPTLIQ IS THE RATE OF GENERATION (FALLOUT) OF
    1260. 

  1260. !...LIQUID PRECIP AT A GIVEN MODEL LVL, PPTICE THE SAME FOR ICE,
    1261. 

  1261. !...TRPPT IS THE TOTAL RATE OF PRODUCTION OF PRECIP UP TO THE
    1262. 

  1262. !...CURRENT MODEL LEVEL...
    1263. 

  1263. !
    1264. 

  1264.                 PPTLIQ(NK1)=QLQOUT(NK1)*UMF(NK)
    1265. 

  1265.                 PPTICE(NK1)=QICOUT(NK1)*UMF(NK)
    1266. 

  1266. !
    1267. 

  1267.                 TRPPT=TRPPT+PPTLIQ(NK1)+PPTICE(NK1)
    1268. 

  1268.                 IF(NK1.LE.KPBL)UER(NK1)=UER(NK1)+VMFLCL*DP(NK1)/DPTHMX
    1269. 

  1269.               ENDIF
    1270. 

  1270. !
    1271. 

  1271.             END DO updraft
    1272. 

  1272. !
    1273. 

  1273. !...CHECK CLOUD DEPTH...IF CLOUD IS TALL ENOUGH, ESTIMATE THE EQUILIBRIUM
    1274. 

  1274. !   TEMPERATURE LEVEL (LET) AND ADJUST MASS FLUX PROFILE AT CLOUD TOP SO
    1275. 

  1275. !   THAT MASS FLUX DECREASES TO ZERO AS A LINEAR FUNCTION OF PRESSURE BETWEEN
    1276. 

  1276. !   THE LET AND CLOUD TOP...
    1277. 

  1277. !     
    1278. 

  1278. !...LTOP IS THE MODEL LEVEL JUST BELOW THE LEVEL AT WHICH VERTICAL VELOCITY
    1279. 

  1279. !   FIRST BECOMES NEGATIVE...
    1280. 

  1280. !     
    1281. 

  1281.             LTOP=NK
    1282. 

  1282.             CLDHGT(LC)=Z0(LTOP)-ZLCL 
    1283. 

  1283. !
    1284. 

  1284. !...Instead of using the same minimum cloud height (for deep convection)
    1285. 

  1285. !...everywhere, try specifying minimum cloud depth as a function of TLCL...
    1286. 

  1286. !
    1287. 

  1287. !   Kain (2004)  Eq. 7
    1288. 

  1288. !
    1289. 

  1289.             IF(TLCL.GT.293.)THEN
    1290. 

  1290.               CHMIN = 4.E3
    1291. 

  1291.             ELSEIF(TLCL.LE.293. .and. TLCL.GE.273)THEN
    1292. 

  1292.               CHMIN = 2.E3 + 100.*(TLCL-273.)
    1293. 

  1293.             ELSEIF(TLCL.LT.273.)THEN
    1294. 

  1294.               CHMIN = 2.E3
    1295. 

  1295.             ENDIF
    1296. 

  1296. 

  1297. !     
    1298. 

  1298. !...If cloud top height is less than the specified minimum for deep 
    1299. 

  1299. !...convection, save value to consider this level as source for 
    1300. 

  1300. !...shallow convection, go back up to check next level...
    1301. 

  1301. !     
    1302. 

  1302. !...Try specifying minimum cloud depth as a function of TLCL...
    1303. 

  1303. !
    1304. 

  1304. !
    1305. 

  1305. !...DO NOT ALLOW ANY CLOUD FROM THIS LAYER IF:
    1306. 

  1306. !
    1307. 

  1307. !...            1.) if there is no CAPE, or 
    1308. 

  1308. !...            2.) cloud top is at model level just above LCL, or
    1309. 

  1309. !...            3.) cloud top is within updraft source layer, or
    1310. 

  1310. !...            4.) cloud-top detrainment layer begins within 
    1311. 

  1311. !...                updraft source layer.
    1312. 

  1312. !
    1313. 

  1313.             IF(LTOP.LE.KLCL .or. LTOP.LE.KPBL .or. LET+1.LE.KPBL)THEN  ! No Convection Allowed
    1314. 

  1314.               CLDHGT(LC)=0.
    1315. 

  1315.               DO NK=K,LTOP
    1316. 

  1316.                 UMF(NK)=0.
    1317. 

  1317.                 UDR(NK)=0.
    1318. 

  1318.                 UER(NK)=0.
    1319. 

  1319.                 DETLQ(NK)=0.
    1320. 

  1320.                 DETIC(NK)=0.
    1321. 

  1321.                 PPTLIQ(NK)=0.
    1322. 

  1322.                 PPTICE(NK)=0.
    1323. 

  1323.               ENDDO
    1324. 

  1324. !        
    1325. 

  1325.             ELSEIF(CLDHGT(LC).GT.CHMIN .and. ABE.GT.1)THEN      ! Deep Convection allowed
    1326. 

  1326.               ISHALL=0
    1327. 

  1327.               EXIT usl
    1328. 

  1328.             ELSE
    1329. 

  1329. !
    1330. 

  1330. !...TO DISALLOW SHALLOW CONVECTION, COMMENT OUT NEXT LINE !!!!!!!!
    1331. 

  1331.               ISHALL = 1
    1332. 

  1332.               IF(NU.EQ.NUCHM)THEN
    1333. 

  1333.                 EXIT usl               ! Shallow Convection from this layer
    1334. 

  1334.               ELSE
    1335. 

  1335. ! Remember this layer (by virtue of non-zero CLDHGT) as potential shallow-cloud layer
    1336. 

  1336.                 DO NK=K,LTOP
    1337. 

  1337.                   UMF(NK)=0.
    1338. 

  1338.                   UDR(NK)=0.
    1339. 

  1339.                   UER(NK)=0.
    1340. 

  1340.                   DETLQ(NK)=0.
    1341. 

  1341.                   DETIC(NK)=0.
    1342. 

  1342.                   PPTLIQ(NK)=0.
    1343. 

  1343.                   PPTICE(NK)=0.
    1344. 

  1344.                 ENDDO
    1345. 

  1345.               ENDIF
    1346. 

  1346.             ENDIF
    1347. 

  1347.           ENDIF trigger2
    1348. 

  1348.         END DO usl
    1349. 

  1349.     IF(ISHALL.EQ.1)THEN
    1350. 

  1350.       KSTART=MAX0(KPBL,KLCL)
    1351. 

  1351.       LET=KSTART
    1352. 

  1352.     endif
    1353. 

  1353. !     
    1354. 

  1354. !...IF THE LET AND LTOP ARE THE SAME, DETRAIN ALL OF THE UPDRAFT MASS FL
    1355. 

  1355. !   THIS LEVEL...
    1356. 

  1356. !     
    1357. 

  1357.     IF(LET.EQ.LTOP)THEN
    1358. 

  1358.       UDR(LTOP)=UMF(LTOP)+UDR(LTOP)-UER(LTOP)
    1359. 

  1359.       DETLQ(LTOP)=QLIQ(LTOP)*UDR(LTOP)*UPNEW/UPOLD
    1360. 

  1360.       DETIC(LTOP)=QICE(LTOP)*UDR(LTOP)*UPNEW/UPOLD
    1361. 

  1361.       UER(LTOP)=0.
    1362. 

  1362.       UMF(LTOP)=0.
    1363. 

  1363.     ELSE 
    1364. 

  1364. !     
    1365. 

  1365. !   BEGIN TOTAL DETRAINMENT AT THE LEVEL ABOVE THE LET...
    1366. 

  1366. !     
    1367. 

  1367.       DPTT=0.
    1368. 

  1368.       DO NJ=LET+1,LTOP
    1369. 

  1369.         DPTT=DPTT+DP(NJ)
    1370. 

  1370.       ENDDO
    1371. 

  1371.       DUMFDP=UMF(LET)/DPTT
    1372. 

  1372. !     
    1373. 

  1373. !...ADJUST MASS FLUX PROFILES, DETRAINMENT RATES, AND PRECIPITATION FALL
    1374. 

  1374. !   RATES TO REFLECT THE LINEAR DECREASE IN MASS FLX BETWEEN THE LET AND
    1375. 

  1375. !     
    1376. 

  1376.       DO NK=LET+1,LTOP
    1377. 

  1377. !
    1378. 

  1378. !...entrainment is allowed at every level except for LTOP, so disallow
    1379. 

  1379. !...entrainment at LTOP and adjust entrainment rates between LET and LTOP
    1380. 

  1380. !...so the the dilution factor due to entyrianment is not changed but 
    1381. 

  1381. !...the actual entrainment rate will change due due forced total 
    1382. 

  1382. !...detrainment in this layer...
    1383. 

  1383. !
    1384. 

  1384.         IF(NK.EQ.LTOP)THEN
    1385. 

  1385.           UDR(NK) = UMF(NK-1)
    1386. 

  1386.           UER(NK) = 0.
    1387. 

  1387.           DETLQ(NK) = UDR(NK)*QLIQ(N…
    1388.

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