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

/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 
    

  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(NK)*DILFRC(NK)
    1388. 

  1388.           DETIC(NK) = UDR(NK)*QICE(NK)*DILFRC(NK)
    1389. 

  1389.         ELSE
    1390. 

  1390.           UMF(NK)=UMF(NK-1)-DP(NK)*DUMFDP
    1391. 

  1391.           UER(NK)=UMF(NK)*(1.-1./DILFRC(NK))
    1392. 

  1392.           UDR(NK)=UMF(NK-1)-UMF(NK)+UER(NK)
    1393. 

  1393.           DETLQ(NK)=UDR(NK)*QLIQ(NK)*DILFRC(NK)
    1394. 

  1394.           DETIC(NK)=UDR(NK)*QICE(NK)*DILFRC(NK)
    1395. 

  1395.         ENDIF
    1396. 

  1396.         IF(NK.GE.LET+2)THEN
    1397. 

  1397.           TRPPT=TRPPT-PPTLIQ(NK)-PPTICE(NK)
    1398. 

  1398.           PPTLIQ(NK)=UMF(NK-1)*QLQOUT(NK)
    1399. 

  1399.           PPTICE(NK)=UMF(NK-1)*QICOUT(NK)
    1400. 

  1400.           TRPPT=TRPPT+PPTLIQ(NK)+PPTICE(NK)
    1401. 

  1401.         ENDIF
    1402. 

  1402.       ENDDO
    1403. 

  1403.     ENDIF
    1404. 

  1404. !     
    1405. 

  1405. ! Initialize some arrays below cloud base and above cloud top...
    1406. 

  1406. !
    1407. 

  1407.     DO NK=1,LTOP
    1408. 

  1408.       IF(T0(NK).GT.T00)ML=NK
    1409. 

  1409.     ENDDO
    1410. 

  1410.     DO NK=1,K
    1411. 

  1411.       IF(NK.GE.LC)THEN
    1412. 

  1412.         IF(NK.EQ.LC)THEN
    1413. 

  1413.           UMF(NK)=VMFLCL*DP(NK)/DPTHMX
    1414. 

  1414.           UER(NK)=VMFLCL*DP(NK)/DPTHMX
    1415. 

  1415.         ELSEIF(NK.LE.KPBL)THEN
    1416. 

  1416.           UER(NK)=VMFLCL*DP(NK)/DPTHMX
    1417. 

  1417.           UMF(NK)=UMF(NK-1)+UER(NK)
    1418. 

  1418.         ELSE
    1419. 

  1419.           UMF(NK)=VMFLCL
    1420. 

  1420.           UER(NK)=0.
    1421. 

  1421.         ENDIF
    1422. 

  1422.         TU(NK)=TMIX+(Z0(NK)-ZMIX)*GDRY
    1423. 

  1423.         QU(NK)=QMIX
    1424. 

  1424.         WU(NK)=WLCL
    1425. 

  1425.       ELSE
    1426. 

  1426.         TU(NK)=0.
    1427. 

  1427.         QU(NK)=0.
    1428. 

  1428.         UMF(NK)=0.
    1429. 

  1429.         WU(NK)=0.
    1430. 

  1430.         UER(NK)=0.
    1431. 

  1431.       ENDIF
    1432. 

  1432.       UDR(NK)=0.
    1433. 

  1433.       QDT(NK)=0.
    1434. 

  1434.       QLIQ(NK)=0.
    1435. 

  1435.       QICE(NK)=0.
    1436. 

  1436.       QLQOUT(NK)=0.
    1437. 

  1437.       QICOUT(NK)=0.
    1438. 

  1438.       PPTLIQ(NK)=0.
    1439. 

  1439.       PPTICE(NK)=0.
    1440. 

  1440.       DETLQ(NK)=0.
    1441. 

  1441.       DETIC(NK)=0.
    1442. 

  1442.       RATIO2(NK)=0.
    1443. 

  1443.       CALL ENVIRTHT(P0(NK),T0(NK),Q0(NK),THETEE(NK),ALIQ,BLIQ,CLIQ,DLIQ)
    1444. 

  1444.       EQFRC(NK)=1.0
    1445. 

  1445.     ENDDO
    1446. 

  1446. !     
    1447. 

  1447.       LTOP1=LTOP+1
    1448. 

  1448.       LTOPM1=LTOP-1
    1449. 

  1449. !     
    1450. 

  1450. !...DEFINE VARIABLES ABOVE CLOUD TOP...
    1451. 

  1451. !     
    1452. 

  1452.       DO NK=LTOP1,KX
    1453. 

  1453.         UMF(NK)=0.
    1454. 

  1454.         UDR(NK)=0.
    1455. 

  1455.         UER(NK)=0.
    1456. 

  1456.         QDT(NK)=0.
    1457. 

  1457.         QLIQ(NK)=0.
    1458. 

  1458.         QICE(NK)=0.
    1459. 

  1459.         QLQOUT(NK)=0.
    1460. 

  1460.         QICOUT(NK)=0.
    1461. 

  1461.         DETLQ(NK)=0.
    1462. 

  1462.         DETIC(NK)=0.
    1463. 

  1463.         PPTLIQ(NK)=0.
    1464. 

  1464.         PPTICE(NK)=0.
    1465. 

  1465.         IF(NK.GT.LTOP1)THEN
    1466. 

  1466.           TU(NK)=0.
    1467. 

  1467.           QU(NK)=0.
    1468. 

  1468.           WU(NK)=0.
    1469. 

  1469.         ENDIF
    1470. 

  1470.         THTA0(NK)=0.
    1471. 

  1471.         THTAU(NK)=0.
    1472. 

  1472.         EMS(NK)=0.
    1473. 

  1473.         EMSD(NK)=0.
    1474. 

  1474.         TG(NK)=T0(NK)
    1475. 

  1475.         QG(NK)=Q0(NK)
    1476. 

  1476.         QLG(NK)=0.
    1477. 

  1477.         QIG(NK)=0.
    1478. 

  1478.         QRG(NK)=0.
    1479. 

  1479.         QSG(NK)=0.
    1480. 

  1480.         OMG(NK)=0.
    1481. 

  1481.       ENDDO
    1482. 

  1482.         OMG(KX+1)=0.
    1483. 

  1483.         DO NK=1,LTOP
    1484. 

  1484.           EMS(NK)=DP(NK)*DXSQ/G
    1485. 

  1485.           EMSD(NK)=1./EMS(NK)
    1486. 

  1486. !     
    1487. 

  1487. !...INITIALIZE SOME VARIABLES TO BE USED LATER IN THE VERT ADVECTION SCHEME
    1488. 

  1488. !     
    1489. 

  1489.           EXN(NK)=(P00/P0(NK))**(0.2854*(1.-0.28*QDT(NK)))
    1490. 

  1490.           THTAU(NK)=TU(NK)*EXN(NK)
    1491. 

  1491.           EXN(NK)=(P00/P0(NK))**(0.2854*(1.-0.28*Q0(NK)))
    1492. 

  1492.           THTA0(NK)=T0(NK)*EXN(NK)
    1493. 

  1493.           DDILFRC(NK) = 1./DILFRC(NK)
    1494. 

  1494.           OMG(NK)=0.
    1495. 

  1495.         ENDDO
    1496. 

  1496. !     IF (XTIME.LT.10.)THEN
    1497. 

  1497. !      WRITE(98,1025)KLCL,ZLCL,DTLCL,LTOP,P0(LTOP),IFLAG,
    1498. 

  1498. !    * TMIX-T00,PMIX,QMIX,ABE
    1499. 

  1499. !      WRITE(98,1030)P0(LET)/100.,P0(LTOP)/100.,VMFLCL,PLCL/100.,
    1500. 

  1500. !    * WLCL,CLDHGT
    1501. 

  1501. !     ENDIF
    1502. 

  1502. !     
    1503. 

  1503. !...COMPUTE CONVECTIVE TIME SCALE(TIMEC). THE MEAN WIND AT THE LCL
    1504. 

  1504. !...AND MIDTROPOSPHERE IS USED.
    1505. 

  1505. !     
    1506. 

  1506.         WSPD(KLCL)=SQRT(U0(KLCL)*U0(KLCL)+V0(KLCL)*V0(KLCL))
    1507. 

  1507.         WSPD(L5)=SQRT(U0(L5)*U0(L5)+V0(L5)*V0(L5))
    1508. 

  1508.         WSPD(LTOP)=SQRT(U0(LTOP)*U0(LTOP)+V0(LTOP)*V0(LTOP))
    1509. 

  1509.         VCONV=.5*(WSPD(KLCL)+WSPD(L5))
    1510. 

  1510. !...for ETA model, DX is a function of location...
    1511. 

  1511. !       TIMEC=DX(I,J)/VCONV
    1512. 

  1512.         TIMEC=DX/VCONV
    1513. 

  1513.         TADVEC=TIMEC
    1514. 

  1514.         TIMEC=AMAX1(1800.,TIMEC)         ! 30 minutes  >= TIMEC <= 60 minutes
    1515. 

  1515.         TIMEC=AMIN1(3600.,TIMEC)
    1516. 

  1516.         IF(ISHALL.EQ.1)TIMEC=2400.       ! shallow convection TIMEC = 40 minutes
    1517. 

  1517.         NIC=NINT(TIMEC/DT)
    1518. 

  1518.         TIMEC=FLOAT(NIC)*DT
    1519. 

  1519. !     
    1520. 

  1520. !...COMPUTE WIND SHEAR AND PRECIPITATION EFFICIENCY.
    1521. 

  1521. !     
    1522. 

  1522.         IF(WSPD(LTOP).GT.WSPD(KLCL))THEN
    1523. 

  1523.           SHSIGN=1.
    1524. 

  1524.         ELSE
    1525. 

  1525.           SHSIGN=-1.
    1526. 

  1526.         ENDIF
    1527. 

  1527.         VWS=(U0(LTOP)-U0(KLCL))*(U0(LTOP)-U0(KLCL))+(V0(LTOP)-V0(KLCL))*   &
    1528. 

  1528.             (V0(LTOP)-V0(KLCL))
    1529. 

  1529.         VWS=1.E3*SHSIGN*SQRT(VWS)/(Z0(LTOP)-Z0(LCL))
    1530. 

  1530.         PEF=1.591+VWS*(-.639+VWS*(9.53E-2-VWS*4.96E-3))
    1531. 

  1531.         PEF=AMAX1(PEF,.2)
    1532. 

  1532.         PEF=AMIN1(PEF,.9)
    1533. 

  1533. !     
    1534. 

  1534. !...PRECIPITATION EFFICIENCY IS A FUNCTION OF THE HEIGHT OF CLOUD BASE.
    1535. 

  1535. !     
    1536. 

  1536.         CBH=(ZLCL-Z0(1))*3.281E-3
    1537. 

  1537.         IF(CBH.LT.3.)THEN
    1538. 

  1538.           RCBH=.02
    1539. 

  1539.         ELSE
    1540. 

  1540.           RCBH=.96729352+CBH*(-.70034167+CBH*(.162179896+CBH*(-            &
    1541. 

  1541.                1.2569798E-2+CBH*(4.2772E-4-CBH*5.44E-6))))
    1542. 

  1542.         ENDIF
    1543. 

  1543.         IF(CBH.GT.25)RCBH=2.4
    1544. 

  1544.         PEFCBH=1./(1.+RCBH)
    1545. 

  1545.         PEFCBH=AMIN1(PEFCBH,.9)
    1546. 

  1546. !     
    1547. 

  1547. !... MEAN PEF. IS USED TO COMPUTE RAINFALL.
    1548. 

  1548. !     
    1549. 

  1549.         PEFF=.5*(PEF+PEFCBH)
    1550. 

  1550.         PEFF2 = PEFF                                ! JSK MODS
    1551. 

  1551.        IF(IPRNT)THEN  
    1552. 

  1552. !         WRITE(98,1035)PEF,PEFCBH,LC,LET,WKL,VWS
    1553. 

  1553.          WRITE(message,1035)PEF,PEFCBH,LC,LET,WKL,VWS
    1554. 

  1554.          CALL wrf_message( message )
    1555. 

  1555. !       call flush(98)   
    1556. 

  1556.        endif     
    1557. 

  1557. !        WRITE(98,1035)PEF,PEFCBH,LC,LET,WKL,VWS
    1558. 

  1558. !*****************************************************************
    1559. 

  1559. !                                                                *
    1560. 

  1560. !                  COMPUTE DOWNDRAFT PROPERTIES                  *
    1561. 

  1561. !                                                                *
    1562. 

  1562. !*****************************************************************
    1563. 

  1563. !     
    1564. 

  1564. !     
    1565. 

  1565.        TDER=0.
    1566. 

  1566.  devap:IF(ISHALL.EQ.1)THEN
    1567. 

  1567.          LFS = 1
    1568. 

  1568.        ELSE
    1569. 

  1569. !
    1570. 

  1570. !...start downdraft about 150 mb above cloud base...
    1571. 

  1571. !
    1572. 

  1572. !        KSTART=MAX0(KPBL,KLCL)
    1573. 

  1573. !        KSTART=KPBL                                  ! Changed 7/23/99
    1574. 

  1574.          KSTART=KPBL+1                                ! Changed 7/23/99
    1575. 

  1575.          KLFS = LET-1
    1576. 

  1576.          DO NK = KSTART+1,KL
    1577. 

  1577.            DPPP = P0(KSTART)-P0(NK)
    1578. 

  1578. !          IF(DPPP.GT.200.E2)THEN
    1579. 

  1579.            IF(DPPP.GT.150.E2)THEN
    1580. 

  1580.              KLFS = NK
    1581. 

  1581.              EXIT 
    1582. 

  1582.            ENDIF
    1583. 

  1583.          ENDDO
    1584. 

  1584.          KLFS = MIN0(KLFS,LET-1)
    1585. 

  1585.          LFS = KLFS
    1586. 

  1586. !
    1587. 

  1587. !...if LFS is not at least 50 mb above cloud base (implying that the 
    1588. 

  1588. !...level of equil temp, LET, is just above cloud base) do not allow a
    1589. 

  1589. !...downdraft...
    1590. 

  1590. !
    1591. 

  1591.         IF((P0(KSTART)-P0(LFS)).GT.50.E2)THEN
    1592. 

  1592.           THETED(LFS) = THETEE(LFS)
    1593. 

  1593.           QD(LFS) = Q0(LFS)
    1594. 

  1594. !
    1595. 

  1595. !...call tpmix2dd to find wet-bulb temp, qv...
    1596. 

  1596. !
    1597. 

  1597.           call tpmix2dd(p0(lfs),theted(lfs),tz(lfs),qss,i,j)
    1598. 

  1598.           THTAD(LFS)=TZ(LFS)*(P00/P0(LFS))**(0.2854*(1.-0.28*QSS))
    1599. 

  1599. !     
    1600. 

  1600. !...TAKE A FIRST GUESS AT THE INITIAL DOWNDRAFT MASS FLUX...
    1601. 

  1601. !     
    1602. 

  1602.           TVD(LFS)=TZ(LFS)*(1.+0.608*QSS)
    1603. 

  1603.           RDD=P0(LFS)/(R*TVD(LFS))
    1604. 

  1604.           A1=(1.-PEFF)*AU0
    1605. 

  1605.           DMF(LFS)=-A1*RDD
    1606. 

  1606.           DER(LFS)=DMF(LFS)
    1607. 

  1607.           DDR(LFS)=0.
    1608. 

  1608.           RHBAR = RH(LFS)*DP(LFS)
    1609. 

  1609.           DPTT = DP(LFS)
    1610. 

  1610.           DO ND = LFS-1,KSTART,-1
    1611. 

  1611.             ND1 = ND+1
    1612. 

  1612.             DER(ND)=DER(LFS)*EMS(ND)/EMS(LFS)
    1613. 

  1613.             DDR(ND)=0.
    1614. 

  1614.             DMF(ND)=DMF(ND1)+DER(ND)
    1615. 

  1615.             THETED(ND)=(THETED(ND1)*DMF(ND1)+THETEE(ND)*DER(ND))/DMF(ND)
    1616. 

  1616.             QD(ND)=(QD(ND1)*DMF(ND1)+Q0(ND)*DER(ND))/DMF(ND)    
    1617. 

  1617.             DPTT = DPTT+DP(ND)
    1618. 

  1618.             RHBAR = RHBAR+RH(ND)*DP(ND)
    1619. 

  1619.           ENDDO
    1620. 

  1620.           RHBAR = RHBAR/DPTT
    1621. 

  1621.           DMFFRC = 2.*(1.-RHBAR)    ! Kain (2004) eq. 11
    1622. 

  1622.           DPDD = 0.
    1623. 

  1623. !...Calculate melting effect
    1624. 

  1624. !... first, compute total frozen precipitation generated...
    1625. 

  1625. !
    1626. 

  1626.           pptmlt = 0.
    1627. 

  1627.           DO NK = KLCL,LTOP
    1628. 

  1628.             PPTMLT = PPTMLT+PPTICE(NK)
    1629. 

  1629.           ENDDO
    1630. 

  1630.           if(lc.lt.ml)then
    1631. 

  1631. !...For now, calculate melting effect as if DMF = -UMF at KLCL, i.e., as
    1632. 

  1632. !...if DMFFRC=1.  Otherwise, for small DMFFRC, DTMELT gets too large!
    1633. 

  1633. !...12/14/98 jsk...
    1634. 

  1634.             DTMELT = RLF*PPTMLT/(CP*UMF(KLCL))
    1635. 

  1635.           else
    1636. 

  1636.             DTMELT = 0.
    1637. 

  1637.           endif
    1638. 

  1638.           LDT = MIN0(LFS-1,KSTART-1)
    1639. 

  1639. !
    1640. 

  1640.           call tpmix2dd(p0(kstart),theted(kstart),tz(kstart),qss,i,j)
    1641. 

  1641. !
    1642. 

  1642.           tz(kstart) = tz(kstart)-dtmelt
    1643. 

  1643.           ES=ALIQ*EXP((BLIQ*TZ(KSTART)-CLIQ)/(TZ(KSTART)-DLIQ))
    1644. 

  1644.           QSS=0.622*ES/(P0(KSTART)-ES)
    1645. 

  1645.           THETED(KSTART)=TZ(KSTART)*(1.E5/P0(KSTART))**(0.2854*(1.-0.28*QSS))*    &
    1646. 

  1646.                 EXP((3374.6525/TZ(KSTART)-2.5403)*QSS*(1.+0.81*QSS))
    1647. 

  1647. !....  
    1648. 

  1648.           LDT = MIN0(LFS-1,KSTART-1)
    1649. 

  1649.           DO ND = LDT,1,-1
    1650. 

  1650.             DPDD = DPDD+DP(ND)
    1651. 

  1651.             THETED(ND) = THETED(KSTART)
    1652. 

  1652.             QD(ND)     = QD(KSTART)       
    1653. 

  1653. !
    1654. 

  1654. !...call tpmix2dd to find wet bulb temp, saturation mixing ratio...
    1655. 

  1655. !
    1656. 

  1656.             call tpmix2dd(p0(nd),theted(nd),tz(nd),qss,i,j)
    1657. 

  1657.             qsd(nd) = qss
    1658. 

  1658. !
    1659. 

  1659. !...specify RH decrease of 20%/km in downdraft...
    1660. 

  1660. !
    1661. 

  1661.             RHH = 1.-0.2/1000.*(Z0(KSTART)-Z0(ND))
    1662. 

  1662. !
    1663. 

  1663. !...adjust downdraft TEMP, Q to specified RH:
    1664. 

  1664. !
    1665. 

  1665.             IF(RHH.LT.1.)THEN
    1666. 

  1666.               DSSDT=(CLIQ-BLIQ*DLIQ)/((TZ(ND)-DLIQ)*(TZ(ND)-DLIQ))
    1667. 

  1667.               RL=XLV0-XLV1*TZ(ND)
    1668. 

  1668.               DTMP=RL*QSS*(1.-RHH)/(CP+RL*RHH*QSS*DSSDT)
    1669. 

  1669.               T1RH=TZ(ND)+DTMP
    1670. 

  1670.               ES=RHH*ALIQ*EXP((BLIQ*T1RH-CLIQ)/(T1RH-DLIQ))
    1671. 

  1671.               QSRH=0.622*ES/(P0(ND)-ES)
    1672. 

  1672. !
    1673. 

  1673. !...CHECK TO SEE IF MIXING RATIO AT SPECIFIED RH IS LESS THAN ACTUAL
    1674. 

  1674. !...MIXING RATIO...IF SO, ADJUST TO GIVE ZERO EVAPORATION...
    1675. 

  1675. !
    1676. 

  1676.               IF(QSRH.LT.QD(ND))THEN
    1677. 

  1677.                 QSRH=QD(ND)
    1678. 

  1678.                 T1RH=TZ(ND)+(QSS-QSRH)*RL/CP
    1679. 

  1679.               ENDIF
    1680. 

  1680.               TZ(ND)=T1RH
    1681. 

  1681.               QSS=QSRH
    1682. 

  1682.               QSD(ND) = QSS
    1683. 

  1683.             ENDIF         
    1684. 

  1684.             TVD(nd) = tz(nd)*(1.+0.608*qsd(nd))
    1685. 

  1685.             IF(TVD(ND).GT.TV0(ND).OR.ND.EQ.1)THEN
    1686. 

  1686.               LDB=ND
    1687. 

  1687.               EXIT
    1688. 

  1688.             ENDIF
    1689. 

  1689.           ENDDO
    1690. 

  1690.           IF((P0(LDB)-P0(LFS)) .gt. 50.E2)THEN   ! minimum Downdraft depth! 
    1691. 

  1691.             DO ND=LDT,LDB,-1
    1692. 

  1692.               ND1 = ND+1
    1693. 

  1693.               DDR(ND) = -DMF(KSTART)*DP(ND)/DPDD
    1694. 

  1694.               DER(ND) = 0.
    1695. 

  1695.               DMF(ND) = DMF(ND1)+DDR(ND)
    1696. 

  1696.               TDER=TDER+(QSD(nd)-QD(ND))*DDR(ND)
    1697. 

  1697.               QD(ND)=QSD(nd)
    1698. 

  1698.               THTAD(ND)=TZ(ND)*(P00/P0(ND))**(0.2854*(1.-0.28*QD(ND)))
    1699. 

  1699.             ENDDO
    1700. 

  1700.           ENDIF
    1701. 

  1701.         ENDIF
    1702. 

  1702.       ENDIF devap
    1703. 

  1703. !
    1704. 

  1704. !...IF DOWNDRAFT DOES NOT EVAPORATE ANY WATER FOR SPECIFIED RELATIVE
    1705. 

  1705. !...HUMIDITY, NO DOWNDRAFT IS ALLOWED...
    1706. 

  1706. !
    1707. 

  1707. d_mf:   IF(TDER.LT.1.)THEN
    1708. 

  1708. !           WRITE(98,3004)I,J 
    1709. 

  1709. !3004       FORMAT(' ','No Downdraft!;  I=',I3,2X,'J=',I3,'ISHALL =',I2)
    1710. 

  1710.           PPTFLX=TRPPT
    1711. 

  1711.           CPR=TRPPT
    1712. 

  1712.           TDER=0.
    1713. 

  1713.           CNDTNF=0.
    1714. 

  1714.           UPDINC=1.
    1715. 

  1715.           LDB=LFS
    1716. 

  1716.           DO NDK=1,LTOP
    1717. 

  1717.             DMF(NDK)=0.
    1718. 

  1718.             DER(NDK)=0.
    1719. 

  1719.             DDR(NDK)=0.
    1720. 

  1720.             THTAD(NDK)=0.
    1721. 

  1721.             WD(NDK)=0.
    1722. 

  1722.             TZ(NDK)=0.
    1723. 

  1723.             QD(NDK)=0.
    1724. 

  1724.           ENDDO
    1725. 

  1725.           AINCM2=100.
    1726. 

  1726.         ELSE 
    1727. 

  1727.           DDINC = -DMFFRC*UMF(KLCL)/DMF(KSTART)
    1728. 

  1728.           UPDINC=1.
    1729. 

  1729.           IF(TDER*DDINC.GT.TRPPT)THEN
    1730. 

  1730.             DDINC = TRPPT/TDER
    1731. 

  1731.           ENDIF
    1732. 

  1732.           TDER = TDER*DDINC
    1733. 

  1733.           DO NK=LDB,LFS
    1734. 

  1734.             DMF(NK)=DMF(NK)*DDINC
    1735. 

  1735.             DER(NK)=DER(NK)*DDINC
    1736. 

  1736.             DDR(NK)=DDR(NK)*DDINC
    1737. 

  1737.           ENDDO
    1738. 

  1738.          CPR=TRPPT
    1739. 

  1739.          PPTFLX = TRPPT-TDER
    1740. 

  1740.          PEFF=PPTFLX/TRPPT
    1741. 

  1741.          IF(IPRNT)THEN
    1742. 

  1742. !           write(98,*)'PRECIP EFFICIENCY =',PEFF
    1743. 

  1743.            write(message,*)'PRECIP EFFICIENCY =',PEFF
    1744. 

  1744.            CALL wrf_message(message)
    1745. 

  1745. !          call flush(98)   
    1746. 

  1746.          ENDIF
    1747. 

  1747. !
    1748. 

  1748. !
    1749. 

  1749. !...ADJUST UPDRAFT MASS FLUX, MASS DETRAINMENT RATE, AND LIQUID WATER AN
    1750. 

  1750. !   DETRAINMENT RATES TO BE CONSISTENT WITH THE TRANSFER OF THE ESTIMATE
    1751. 

  1751. !   FROM THE UPDRAFT TO THE DOWNDRAFT AT THE LFS...
    1752. 

  1752. !     
    1753. 

  1753. !         DO NK=LC,LFS
    1754. 

  1754. !           UMF(NK)=UMF(NK)*UPDINC
    1755. 

  1755. !           UDR(NK)=UDR(NK)*UPDINC
    1756. 

  1756. !           UER(NK)=UER(NK)*UPDINC
    1757. 

  1757. !           PPTLIQ(NK)=PPTLIQ(NK)*UPDINC
    1758. 

  1758. !           PPTICE(NK)=PPTICE(NK)*UPDINC
    1759. 

  1759. !           DETLQ(NK)=DETLQ(NK)*UPDINC
    1760. 

  1760. !           DETIC(NK)=DETIC(NK)*UPDINC
    1761. 

  1761. !         ENDDO
    1762. 

  1762. !     
    1763. 

  1763. !...ZERO OUT THE ARRAYS FOR DOWNDRAFT DATA AT LEVELS ABOVE AND BELOW THE
    1764. 

  1764. !...DOWNDRAFT...
    1765. 

  1765. !     
    1766. 

  1766.          IF(LDB.GT.1)THEN
    1767. 

  1767.            DO NK=1,LDB-1
    1768. 

  1768.              DMF(NK)=0.
    1769. 

  1769.              DER(NK)=0.
    1770. 

  1770.              DDR(NK)=0.
    1771. 

  1771.              WD(NK)=0.
    1772. 

  1772.              TZ(NK)=0.
    1773. 

  1773.              QD(NK)=0.
    1774. 

  1774.              THTAD(NK)=0.
    1775. 

  1775.            ENDDO
    1776. 

  1776.          ENDIF
    1777. 

  1777.          DO NK=LFS+1,KX
    1778. 

  1778.            DMF(NK)=0.
    1779. 

  1779.            DER(NK)=0.
    1780. 

  1780.            DDR(NK)=0.
    1781. 

  1781.            WD(NK)=0.
    1782. 

  1782.            TZ(NK)=0.
    1783. 

  1783.            QD(NK)=0.
    1784. 

  1784.            THTAD(NK)=0.
    1785. 

  1785.          ENDDO
    1786. 

  1786.          DO NK=LDT+1,LFS-1
    1787. 

  1787.            TZ(NK)=0.
    1788. 

  1788.            QD(NK)=0.
    1789. 

  1789.            THTAD(NK)=0.
    1790. 

  1790.          ENDDO
    1791. 

  1791.        ENDIF d_mf
    1792. 

  1792. !
    1793. 

  1793. !...SET LIMITS ON THE UPDRAFT AND DOWNDRAFT MASS FLUXES SO THAT THE INFLOW
    1794. 

  1794. !   INTO CONVECTIVE DRAFTS FROM A GIVEN LAYER IS NO MORE THAN IS AVAILABLE
    1795. 

  1795. !   IN THAT LAYER INITIALLY...
    1796. 

  1796. !     
    1797. 

  1797.        AINCMX=1000.
    1798. 

  1798.        LMAX=MAX0(KLCL,LFS)
    1799. 

  1799.        DO NK=LC,LMAX
    1800. 

  1800.          IF((UER(NK)-DER(NK)).GT.1.e-3)THEN
    1801. 

  1801.            AINCM1=EMS(NK)/((UER(NK)-DER(NK))*TIMEC)
    1802. 

  1802.            AINCMX=AMIN1(AINCMX,AINCM1)
    1803. 

  1803.          ENDIF
    1804. 

  1804.        ENDDO
    1805. 

  1805.        AINC=1.
    1806. 

  1806.        IF(AINCMX.LT.AINC)AINC=AINCMX
    1807. 

  1807. !     
    1808. 

  1808. !...SAVE THE RELEVENT VARIABLES FOR A UNIT UPDRAFT AND DOWNDRAFT...THEY WILL 
    1809. 

  1809. !...BE ITERATIVELY ADJUSTED BY THE FACTOR AINC TO SATISFY THE STABILIZATION
    1810. 

  1810. !...CLOSURE...
    1811. 

  1811. !     
    1812. 

  1812.        TDER2=TDER
    1813. 

  1813.        PPTFL2=PPTFLX
    1814. 

  1814.        DO NK=1,LTOP
    1815. 

  1815.          DETLQ2(NK)=DETLQ(NK)
    1816. 

  1816.          DETIC2(NK)=DETIC(NK)
    1817. 

  1817.          UDR2(NK)=UDR(NK)
    1818. 

  1818.          UER2(NK)=UER(NK)
    1819. 

  1819.          DDR2(NK)=DDR(NK)
    1820. 

  1820.          DER2(NK)=DER(NK)
    1821. 

  1821.          UMF2(NK)=UMF(NK)
    1822. 

  1822.          DMF2(NK)=DMF(NK)
    1823. 

  1823.        ENDDO
    1824. 

  1824.        FABE=1.
    1825. 

  1825.        STAB=0.95
    1826. 

  1826.        NOITR=0
    1827. 

  1827.        ISTOP=0
    1828. 

  1828. !
    1829. 

  1829.         IF(ISHALL.EQ.1)THEN                              ! First for shallow convection
    1830. 

  1830. !
    1831. 

  1831. ! No iteration for shallow convection; if turbulent kinetic energy (TKE) is available
    1832. 

  1832. ! from a turbulence parameterization, scale cloud-base updraft mass flux as a function
    1833. 

  1833. ! of TKE, but for now, just specify shallow-cloud mass flux using TKEMAX = 5...
    1834. 

  1834. !
    1835. 

  1835. !...find the maximum TKE value between LC and KLCL...
    1836. 

  1836. !         TKEMAX = 0.
    1837. 

  1837.           TKEMAX = 5.
    1838. 

  1838. !          DO 173 K = LC,KLCL
    1839. 

  1839. !            NK = KX-K+1
    1840. 

  1840. !            TKEMAX = AMAX1(TKEMAX,Q2(I,J,NK))
    1841. 

  1841. ! 173      CONTINUE
    1842. 

  1842. !          TKEMAX = AMIN1(TKEMAX,10.)
    1843. 

  1843. !          TKEMAX = AMAX1(TKEMAX,5.)
    1844. 

  1844. !c         TKEMAX = 10.
    1845. 

  1845. !c...3_24_99...DPMIN was changed for shallow convection so that it is the
    1846. 

  1846. !c...          the same as for deep convection (5.E3).  Since this doubles
    1847. 

  1847. !c...          (roughly) the value of DPTHMX, add a factor of 0.5 to calcu-
    1848. 

  1848. !c...          lation of EVAC...
    1849. 

  1849. !c         EVAC  = TKEMAX*0.1
    1850. 

  1850.           EVAC  = 0.5*TKEMAX*0.1
    1851. 

  1851. !         AINC = 0.1*DPTHMX*DXIJ*DXIJ/(VMFLCL*G*TIMEC)
    1852. 

  1852. !          AINC = EVAC*DPTHMX*DX(I,J)*DX(I,J)/(VMFLCL*G*TIMEC)
    1853. 

  1853.           AINC = EVAC*DPTHMX*DXSQ/(VMFLCL*G*TIMEC)
    1854. 

  1854.           TDER=TDER2*AINC
    1855. 

  1855.           PPTFLX=PPTFL2*AINC
    1856. 

  1856.           DO NK=1,LTOP
    1857. 

  1857.             UMF(NK)=UMF2(NK)*AINC
    1858. 

  1858.             DMF(NK)=DMF2(NK)*AINC
    1859. 

  1859.             DETLQ(NK)=DETLQ2(NK)*AINC
    1860. 

  1860.             DETIC(NK)=DETIC2(NK)*AINC
    1861. 

  1861.             UDR(NK)=UDR2(NK)*AINC
    1862. 

  1862.             UER(NK)=UER2(NK)*AINC
    1863. 

  1863.             DER(NK)=DER2(NK)*AINC
    1864. 

  1864.             DDR(NK)=DDR2(NK)*AINC
    1865. 

  1865.           ENDDO
    1866. 

  1866.         ENDIF                                           ! Otherwise for deep convection
    1867. 

  1867. ! use iterative procedure to find mass fluxes...
    1868. 

  1868. iter:     DO NCOUNT=1,10
    1869. 

  1869. !     
    1870. 

  1870. !*****************************************************************
    1871. 

  1871. !                                                                *
    1872. 

  1872. !           COMPUTE PROPERTIES FOR COMPENSATIONAL SUBSIDENCE     *
    1873. 

  1873. !                                                                *
    1874. 

  1874. !*****************************************************************
    1875. 

  1875. !     
    1876. 

  1876. !...DETERMINE OMEGA VALUE NECESSARY AT TOP AND BOTTOM OF EACH LAYER TO
    1877. 

  1877. !...SATISFY MASS CONTINUITY...
    1878. 

  1878. !     
    1879. 

  1879.             DTT=TIMEC
    1880. 

  1880.             DO NK=1,LTOP
    1881. 

  1881.               DOMGDP(NK)=-(UER(NK)-DER(NK)-UDR(NK)-DDR(NK))*EMSD(NK)
    1882. 

  1882.               IF(NK.GT.1)THEN
    1883. 

  1883.                 OMG(NK)=OMG(NK-1)-DP(NK-1)*DOMGDP(NK-1)
    1884. 

  1884.                 ABSOMG = ABS(OMG(NK))
    1885. 

  1885.                 ABSOMGTC = ABSOMG*TIMEC
    1886. 

  1886.                 FRDP = 0.75*DP(NK-1)
    1887. 

  1887.                 IF(ABSOMGTC.GT.FRDP)THEN
    1888. 

  1888.                   DTT1 = FRDP/ABSOMG
    1889. 

  1889.                   DTT=AMIN1(DTT,DTT1)
    1890. 

  1890.                 ENDIF
    1891. 

  1891.               ENDIF
    1892. 

  1892.             ENDDO
    1893. 

  1893.             DO NK=1,LTOP
    1894. 

  1894.               THPA(NK)=THTA0(NK)
    1895. 

  1895.               QPA(NK)=Q0(NK)
    1896. 

  1896.               NSTEP=NINT(TIMEC/DTT+1)
    1897. 

  1897.               DTIME=TIMEC/FLOAT(NSTEP)
    1898. 

  1898.               FXM(NK)=OMG(NK)*DXSQ/G
    1899. 

  1899.             ENDDO
    1900. 

  1900. !     
    1901. 

  1901. !...DO AN UPSTREAM/FORWARD-IN-TIME ADVECTION OF THETA, QV...
    1902. 

  1902. !     
    1903. 

  1903.         DO NTC=1,NSTEP
    1904. 

  1904. !     
    1905. 

  1905. !...ASSIGN THETA AND Q VALUES AT THE TOP AND BOTTOM OF EACH LAYER BASED ON
    1906. 

  1906. !...SIGN OF OMEGA...
    1907. 

  1907. !     
    1908. 

  1908.             DO  NK=1,LTOP
    1909. 

  1909.               THFXIN(NK)=0.
    1910. 

  1910.               THFXOUT(NK)=0.
    1911. 

  1911.               QFXIN(NK)=0.
    1912. 

  1912.               QFXOUT(NK)=0.
    1913. 

  1913.             ENDDO
    1914. 

  1914.             DO NK=2,LTOP
    1915. 

  1915.               IF(OMG(NK).LE.0.)THEN
    1916. 

  1916.                 THFXIN(NK)=-FXM(NK)*THPA(NK-1)
    1917. 

  1917.                 QFXIN(NK)=-FXM(NK)*QPA(NK-1)
    1918. 

  1918.                 THFXOUT(NK-1)=THFXOUT(NK-1)+THFXIN(NK)
    1919. 

  1919.                 QFXOUT(NK-1)=QFXOUT(NK-1)+QFXIN(NK)
    1920. 

  1920.               ELSE
    1921. 

  1921.                 THFXOUT(NK)=FXM(NK)*THPA(NK)
    1922. 

  1922.                 QFXOUT(NK)=FXM(NK)*QPA(NK)
    1923. 

  1923.                 THFXIN(NK-1)=THFXIN(NK-1)+THFXOUT(NK)
    1924. 

  1924.                 QFXIN(NK-1)=QFXIN(NK-1)+QFXOUT(NK)
    1925. 

  1925.               ENDIF
    1926. 

  1926.             ENDDO
    1927. 

  1927. !     
    1928. 

  1928. !...UPDATE THE THETA AND QV VALUES AT EACH LEVEL...
    1929. 

  1929. !     
    1930. 

  1930.             DO NK=1,LTOP
    1931. 

  1931.               THPA(NK)=THPA(NK)+(THFXIN(NK)+UDR(NK)*THTAU(NK)+DDR(NK)*      &
    1932. 

  1932.                        THTAD(NK)-THFXOUT(NK)-(UER(NK)-DER(NK))*THTA0(NK))*  &
    1933. 

  1933.                        DTIME*EMSD(NK)
    1934. 

  1934.               QPA(NK)=QPA(NK)+(QFXIN(NK)+UDR(NK)*QDT(NK)+DDR(NK)*QD(NK)-    &
    1935. 

  1935.                       QFXOUT(NK)-(UER(NK)-DER(NK))*Q0(NK))*DTIME*EMSD(NK)
    1936. 

  1936.             ENDDO   
    1937. 

  1937.           ENDDO   
    1938. 

  1938.           DO NK=1,LTOP
    1939. 

  1939.             THTAG(NK)=THPA(NK)
    1940. 

  1940.             QG(NK)=QPA(NK)
    1941. 

  1941.           ENDDO
    1942. 

  1942. !     
    1943. 

  1943. !...CHECK TO SEE IF MIXING RATIO DIPS BELOW ZERO ANYWHERE;  IF SO, BORROW
    1944. 

  1944. !...MOISTURE FROM ADJACENT LAYERS TO BRING IT BACK UP ABOVE ZERO...
    1945. 

  1945. !     
    1946. 

  1946.         DO NK=1,LTOP
    1947. 

  1947.           IF(QG(NK).LT.0.)THEN
    1948. 

  1948.             IF(NK.EQ.1)THEN                             ! JSK MODS
    1949. 

  1949. !              PRINT *,' PROBLEM WITH KF SCHEME:  ' ! JSK MODS
    1950. 

  1950. !              PRINT *,'QG = 0 AT THE SURFACE!!!!!!!'    ! JSK MODS
    1951. 

  1951.               CALL wrf_error_fatal ( 'QG, QG(NK).LT.0') ! JSK MODS
    1952. 

  1952.             ENDIF                                       ! JSK MODS
    1953. 

  1953.             NK1=NK+1
    1954. 

  1954.             IF(NK.EQ.LTOP)THEN
    1955. 

  1955.               NK1=KLCL
    1956. 

  1956.             ENDIF
    1957. 

  1957.             TMA=QG(NK1)*EMS(NK1)
    1958. 

  1958.             TMB=QG(NK-1)*EMS(NK-1)
    1959. 

  1959.             TMM=(QG(NK)-1.E-9)*EMS(NK  )
    1960. 

  1960.             BCOEFF=-TMM/((TMA*TMA)/TMB+TMB)
    1961. 

  1961.             ACOEFF=BCOEFF*TMA/TMB
    1962. 

  1962.             TMB=TMB*(1.-BCOEFF)
    1963. 

  1963.             TMA=TMA*(1.-ACOEFF)
    1964. 

  1964.             IF(NK.EQ.LTOP)THEN
    1965. 

  1965.               QVDIFF=(QG(NK1)-TMA*EMSD(NK1))*100./QG(NK1)
    1966. 

  1966. !              IF(ABS(QVDIFF).GT.1.)THEN
    1967. 

  1967. !             PRINT *,'!!!WARNING!!! CLOUD BASE WATER VAPOR CHANGES BY ',     &
    1968. 

  1968. !                      QVDIFF,                                                &
    1969. 

  1969. !                     '% WHEN MOISTURE IS BORROWED TO PREVENT NEGATIVE ',     &
    1970. 

  1970. !                     'VALUES IN KAIN-FRITSCH'
    1971. 

  1971. !              ENDIF
    1972. 

  1972.             ENDIF
    1973. 

  1973.             QG(NK)=1.E-9
    1974. 

  1974.             QG(NK1)=TMA*EMSD(NK1)
    1975. 

  1975.             QG(NK-1)=TMB*EMSD(NK-1)
    1976. 

  1976.           ENDIF
    1977. 

  1977.         ENDDO
    1978. 

  1978.         TOPOMG=(UDR(LTOP)-UER(LTOP))*DP(LTOP)*EMSD(LTOP)
    1979. 

  1979.         IF(ABS(TOPOMG-OMG(LTOP)).GT.1.E-3)THEN
    1980. 

  1980. !       WRITE(99,*)'ERROR:  MASS DOES NOT BALANCE IN KF SCHEME;            &
    1981. 

  1981. !      TOPOMG, OMG =',TOPOMG,OMG(LTOP)
    1982. 

  1982. !      TOPOMG, OMG =',TOPOMG,OMG(LTOP)
    1983. 

  1983.           ISTOP=1
    1984. 

  1984.           IPRNT=.TRUE.
    1985. 

  1985.           EXIT iter
    1986. 

  1986.         ENDIF
    1987. 

  1987. !     
    1988. 

  1988. !...CONVERT THETA TO T...
    1989. 

  1989. !     
    1990. 

  1990.         DO NK=1,LTOP
    1991. 

  1991.           EXN(NK)=(P00/P0(NK))**(0.2854*(1.-0.28*QG(NK)))
    1992. 

  1992.           TG(NK)=THTAG(NK)/EXN(NK)
    1993. 

  1993.           TVG(NK)=TG(NK)*(1.+0.608*QG(NK))
    1994. 

  1994.         ENDDO
    1995. 

  1995.         IF(ISHALL.EQ.1)THEN
    1996. 

  1996.           EXIT iter
    1997. 

  1997.         ENDIF
    1998. 

  1998. !     
    1999. 

  1999. !*******************************************************************
    2000. 

  2000. !                                                                  *
    2001. 

  2001. !     COMPUTE NEW CLOUD AND CHANGE IN AVAILABLE BUOYANT ENERGY.    *
    2002. 

  2002. !                                                                  *
    2003. 

  2003. !*******************************************************************
    2004. 

  2004. !     
    2005. 

  2005. !...THE FOLLOWING COMPUTATIONS ARE SIMILAR TO THAT FOR UPDRAFT
    2006. 

  2006. !     
    2007. 

  2007. !        THMIX=0.
    2008. 

  2008.           TMIX=0.
    2009. 

  2009.           QMIX=0.
    2010. 

  2010. !
    2011. 

  2011. !...FIND THE THERMODYNAMIC CHARACTERISTICS OF THE LAYER BY
    2012. 

  2012. !...MASS-WEIGHTING THE CHARACTERISTICS OF THE INDIVIDUAL MODEL
    2013. 

  2013. !...LAYERS...
    2014. 

  2014. !
    2015. 

  2015.           DO NK=LC,KPBL
    2016. 

  2016.             TMIX=TMIX+DP(NK)*TG(NK)
    2017. 

  2017.             QMIX=QMIX+DP(NK)*QG(NK)  
    2018. 

  2018.           ENDDO
    2019. 

  2019.           TMIX=TMIX/DPTHMX
    2020. 

  2020.           QMIX=QMIX/DPTHMX
    2021. 

  2021.           ES=ALIQ*EXP((TMIX*BLIQ-CLIQ)/(TMIX-DLIQ))
    2022. 

  2022.           QSS=0.622*ES/(PMIX-ES)
    2023. 

  2023. !     
    2024. 

  2024. !...REMOVE SUPERSATURATION FOR DIAGNOSTIC PURPOSES, IF NECESSARY...
    2025. 

  2025. !     
    2026. 

  2026.           IF(QMIX.GT.QSS)THEN
    2027. 

  2027.             RL=XLV0-XLV1*TMIX
    2028. 

  2028.             CPM=CP*(1.+0.887*QMIX)
    2029. 

  2029.             DSSDT=QSS*(CLIQ-BLIQ*DLIQ)/((TMIX-DLIQ)*(TMIX-DLIQ))
    2030. 

  2030.             DQ=(QMIX-QSS)/(1.+RL*DSSDT/CPM)
    2031. 

  2031.             TMIX=TMIX+RL/CP*DQ
    2032. 

  2032.             QMIX=QMIX-DQ
    2033. 

  2033.             TLCL=TMIX
    2034. 

  2034.           ELSE
    2035. 

  2035.             QMIX=AMAX1(QMIX,0.)
    2036. 

  2036.             EMIX=QMIX*PMIX/(0.622+QMIX)
    2037. 

  2037.             astrt=1.e-3
    2038. 

  2038.             binc=0.075
    2039. 

  2039.             a1=emix/aliq
    2040. 

  2040.             tp=(a1-astrt)/binc
    2041. 

  2041.             indlu=int(tp)+1
    2042. 

  2042.             value=(indlu-1)*binc+astrt
    2043. 

  2043.             aintrp=(a1-value)/binc
    2044. 

  2044.             tlog=aintrp*alu(indlu+1)+(1-aintrp)*alu(indlu)
    2045. 

  2045.             TDPT=(CLIQ-DLIQ*TLOG)/(BLIQ-TLOG)
    2046. 

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

  2047.             TLCL=AMIN1(TLCL,TMIX)
    2048. 

  2048.           ENDIF
    2049. 

  2049.           TVLCL=TLCL*(1.+0.608*QMIX)
    2050. 

  2050.           ZLCL = ZMIX+(TLCL-TMIX)/GDRY
    2051. 

  2051.           DO NK = LC,KL
    2052. 

  2052.             KLCL=NK
    2053. 

  2053.             IF(ZLCL.LE.Z0(NK))THEN
    2054. 

  2054.               EXIT 
    2055. 

  2055.             ENDIF
    2056. 

  2056.           ENDDO
    2057. 

  2057.           K=KLCL-1
    2058. 

  2058.           DLP=(ZLCL-Z0(K))/(Z0(KLCL)-Z0(K))
    2059. 

  2059. !     
    2060. 

  2060. !...ESTIMATE ENVIRONMENTAL TEMPERATURE AND MIXING RATIO AT THE LCL...
    2061. 

  2061. !     
    2062. 

  2062.           TENV=TG(K)+(TG(KLCL)-TG(K))*DLP
    2063. 

  2063.           QENV=QG(K)+(QG(KLCL)-QG(K))*DLP
    2064. 

  2064.           TVEN=TENV*(1.+0.608*QENV)
    2065. 

  2065.           PLCL=P0(K)+(P0(KLCL)-P0(K))*DLP
    2066. 

  2066.           THETEU(K)=TMIX*(1.E5/PMIX)**(0.2854*(1.-0.28*QMIX))*             &
    2067. 

  2067.                   EXP((3374.6525/TLCL-2.5403)*QMIX*(1.+0.81*QMIX))
    2068. 

  2068. !     
    2069. 

  2069. !...COMPUTE ADJUSTED ABE(ABEG).
    2070. 

  2070. !     
    2071. 

  2071.           ABEG=0.
    2072. 

  2072.           DO NK=K,LTOPM1
    2073. 

  2073.             NK1=NK+1
    2074. 

  2074.             THETEU(NK1) = THETEU(NK)
    2075. 

  2075. !
    2076. 

  2076.             call tpmix2dd(p0(nk1),theteu(nk1),tgu(nk1),qgu(nk1),i,j)
    2077. 

  2077. !
    2078. 

  2078.             TVQU(NK1)=TGU(NK1)*(1.+0.608*QGU(NK1)-QLIQ(NK1)-QICE(NK1))
    2079. 

  2079.             IF(NK.EQ.K)THEN
    2080. 

  2080.               DZZ=Z0(KLCL)-ZLCL
    2081. 

  2081.               DILBE=((TVLCL+TVQU(NK1))/(TVEN+TVG(NK1))-1.)*DZZ
    2082. 

  2082.             ELSE
    2083. 

  2083.               DZZ=DZA(NK)
    2084. 

  2084.               DILBE=((TVQU(NK)+TVQU(NK1))/(TVG(NK)+TVG(NK1))-1.)*DZZ
    2085. 

  2085.             ENDIF
    2086. 

  2086.             IF(DILBE.GT.0.)ABEG=ABEG+DILBE*G
    2087. 

  2087. !
    2088. 

  2088. !...DILUTE BY ENTRAINMENT BY THE RATE AS ORIGINAL UPDRAFT...
    2089. 

  2089. !
    2090. 

  2090.             CALL ENVIRTHT(P0(NK1),TG(NK1),QG(NK1),THTEEG(NK1),ALIQ,BLIQ,CLIQ,DLIQ)
    2091. 

  2091.             THETEU(NK1)=THETEU(NK1)*DDILFRC(NK1)+THTEEG(NK1)*(1.-DDILFRC(NK1))
    2092. 

  2092.           ENDDO
    2093. 

  2093. !     
    2094. 

  2094. !...ASSUME AT LEAST 90% OF CAPE (ABE) IS REMOVED BY CONVECTION DURING
    2095. 

  2095. !...THE PERIOD TIMEC...
    2096. 

  2096. !     
    2097. 

  2097.           IF(NOITR.EQ.1)THEN
    2098. 

  2098. !         write(98,*)' '
    2099. 

  2099. !         write(98,*)'TAU, I, J, =',NTSD,I,J
    2100. 

  2100. !         WRITE(98,1060)FABE
    2101. 

  2101. !          GOTO 265
    2102. 

  2102.           EXIT iter
    2103. 

  2103.           ENDIF
    2104. 

  2104.           DABE=AMAX1(ABE-ABEG,0.1*ABE)
    2105. 

  2105.           FABE=ABEG/ABE
    2106. 

  2106.           IF(FABE.GT.1. .and. ISHALL.EQ.0)THEN
    2107. 

  2107. !          WRITE(98,*)'UPDRAFT/DOWNDRAFT COUPLET INCREASES CAPE AT THIS
    2108. 

  2108. !     *GRID POINT; NO CONVECTION ALLOWED!'
    2109. 

  2109.             RETURN  
    2110. 

  2110.           ENDIF
    2111. 

  2111.           IF(NCOUNT.NE.1)THEN
    2112. 

  2112.             IF(ABS(AINC-AINCOLD).LT.0.0001)THEN
    2113. 

  2113.               NOITR=1
    2114. 

  2114.               AINC=AINCOLD
    2115. 

  2115.               CYCLE iter
    2116. 

  2116.             ENDIF
    2117. 

  2117.             DFDA=(FABE-FABEOLD)/(AINC-AINCOLD)
    2118. 

  2118.             IF(DFDA.GT.0.)THEN
    2119. 

  2119.               NOITR=1
    2120. 

  2120.               AINC=AINCOLD
    2121. 

  2121.               CYCLE iter
    2122. 

  2122.             ENDIF
    2123. 

  2123.           ENDIF
    2124. 

  2124.           AINCOLD=AINC
    2125. 

  2125.           FABEOLD=FABE
    2126. 

  2126.           IF(AINC/AINCMX.GT.0.999.AND.FABE.GT.1.05-STAB)THEN
    2127. 

  2127. !           write(98,*)' '
    2128. 

  2128. !           write(98,*)'TAU, I, J, =',NTSD,I,J
    2129. 

  2129. !           WRITE(98,1055)FABE
    2130. 

  2130. !            GOTO 265
    2131. 

  2131.             EXIT
    2132. 

  2132.           ENDIF
    2133. 

  2133.           IF((FABE.LE.1.05-STAB.AND.FABE.GE.0.95-STAB) .or. NCOUNT.EQ.10)THEN
    2134. 

  2134.             EXIT iter
    2135. 

  2135.           ELSE
    2136. 

  2136.             IF(NCOUNT.GT.10)THEN
    2137. 

  2137. !             write(98,*)' '
    2138. 

  2138. !             write(98,*)'TAU, I, J, =',NTSD,I,J
    2139. 

  2139. !             WRITE(98,1060)FABE
    2140. 

  2140. !             GOTO 265
    2141. 

  2141.               EXIT
    2142. 

  2142.             ENDIF
    2143. 

  2143. !     
    2144. 

  2144. !...IF MORE THAN 10% OF THE ORIGINAL CAPE REMAINS, INCREASE THE CONVECTIVE
    2145. 

  2145. !...MASS FLUX BY THE FACTOR AINC:
    2146. 

  2146. !     
    2147. 

  2147.             IF(FABE.EQ.0.)THEN
    2148. 

  2148.               AINC=AINC*0.5
    2149. 

  2149.             ELSE
    2150. 

  2150.               IF(DABE.LT.1.e-4)THEN
    2151. 

  2151.                 NOITR=1
    2152. 

  2152.                 AINC=AINCOLD
    2153. 

  2153.                 CYCLE iter
    2154. 

  2154.               ELSE
    2155. 

  2155.                 AINC=AINC*STAB*ABE/DABE
    2156. 

  2156.               ENDIF
    2157. 

  2157.             ENDIF
    2158. 

  2158. !           AINC=AMIN1(AINCMX,AINC)
    2159. 

  2159.             AINC=AMIN1(AINCMX,AINC)
    2160. 

  2160. !...IF AINC BECOMES VERY SMALL, EFFECTS OF CONVECTION ! JSK MODS
    2161. 

  2161. !...WILL BE MINIMAL SO JUST IGNORE IT...              ! JSK MODS
    2162. 

  2162.             IF(AINC.LT.0.05)then
    2163. 

  2163.               RETURN                          ! JSK MODS
    2164. 

  2164.             ENDIF
    2165. 

  2165. !            AINC=AMAX1(AINC,0.05)                        ! JSK MODS
    2166. 

  2166.             TDER=TDER2*AINC
    2167. 

  2167.             PPTFLX=PPTFL2*AINC
    2168. 

  2168. !           IF (XTIME.LT.10.)THEN
    2169. 

  2169. !           WRITE(98,1080)LFS,LDB,LDT,TIMEC,TADVEC,NSTEP,NCOUNT,
    2170. 

  2170. !          *              FABEOLD,AINCOLD 
    2171. 

  2171. !           ENDIF
    2172. 

  2172.             DO NK=1,LTOP
    2173. 

  2173.               UMF(NK)=UMF2(NK)*AINC
    2174. 

  2174.               DMF(NK)=DMF2(NK)*AINC
    2175. 

  2175.               DETLQ(NK)=DETLQ2(NK)*AINC
    2176. 

  2176.               DETIC(NK)=DETIC2(NK)*AINC
    2177. 

  2177.               UDR(NK)=UDR2(NK)*AINC
    2178. 

  2178.               UER(NK)=UER2(NK)*AINC
    2179. 

  2179.               DER(NK)=DER2(NK)*AINC
    2180. 

  2180.               DDR(NK)=DDR2(NK)*AINC
    2181. 

  2181.             ENDDO
    2182. 

  2182. !     
    2183. 

  2183. !...GO BACK UP FOR ANOTHER ITERATION...
    2184. 

  2184. !     
    2185. 

  2185.           ENDIF
    2186. 

  2186.         ENDDO iter
    2187. 

  2187. !     
    2188. 

  2188. !...COMPUTE HYDROMETEOR TENDENCIES AS IS DONE FOR T, QV...
    2189. 

  2189. !     
    2190. 

  2190. !...FRC2 IS THE FRACTION OF TOTAL CONDENSATE      !  PPT FB MODS
    2191. 

  2191. !...GENERATED THAT GOES INTO PRECIPITIATION       !  PPT FB MODS
    2192. 

  2192. !
    2193. 

  2193. !  Redistribute hydormeteors according to the final mass-flux values:
    2194. 

  2194. !
    2195. 

  2195.         IF(CPR.GT.0.)THEN 
    2196. 

  2196.           FRC2=PPTFLX/(CPR*AINC)                    !  PPT FB MODS
    2197. 

  2197.         ELSE
    2198. 

  2198.            FRC2=0.
    2199. 

  2199.         ENDIF
    2200. 

  2200.         DO NK=1,LTOP
    2201. 

  2201.           QLPA(NK)=QL0(NK)
    2202. 

  2202.           QIPA(NK)=QI0(NK)
    2203. 

  2203.           QRPA(NK)=QR0(NK)
    2204. 

  2204.           QSPA(NK)=QS0(NK)
    2205. 

  2205.           RAINFB(NK)=PPTLIQ(NK)*AINC*FBFRC*FRC2   !  PPT FB MODS
    2206. 

  2206.           SNOWFB(NK)=PPTICE(NK)*AINC*FBFRC*FRC2   !  PPT FB MODS
    2207. 

  2207.         ENDDO
    2208. 

  2208.         DO NTC=1,NSTEP
    2209. 

  2209. !     
    2210. 

  2210. !...ASSIGN HYDROMETEORS CONCENTRATIONS AT THE TOP AND BOTTOM OF EACH LAYER
    2211. 

  2211. !...BASED ON THE SIGN OF OMEGA...
    2212. 

  2212. !     
    2213. 

  2213.           DO NK=1,LTOP
    2214. 

  2214.             QLFXIN(NK)=0.
    2215. 

  2215.             QLFXOUT(NK)=0.
    2216. 

  2216.             QIFXIN(NK)=0.
    2217. 

  2217.             QIFXOUT(NK)=0.
    2218. 

  2218.             QRFXIN(NK)=0.
    2219. 

  2219.             QRFXOUT(NK)=0.
    2220. 

  2220.             QSFXIN(NK)=0.
    2221. 

  2221.             QSFXOUT(NK)=0.
    2222. 

  2222.           ENDDO   
    2223. 

  2223.           DO NK=2,LTOP
    2224. 

  2224.             IF(OMG(NK).LE.0.)THEN
    2225. 

  2225.               QLFXIN(NK)=-FXM(NK)*QLPA(NK-1)
    2226. 

  2226.               QIFXIN(NK)=-FXM(NK)*QIPA(NK-1)
    2227. 

  2227.               QRFXIN(NK)=-FXM(NK)*QRPA(NK-1)
    2228. 

  2228.               QSFXIN(NK)=-FXM(NK)*QSPA(NK-1)
    2229. 

  2229.               QLFXOUT(NK-1)=QLFXOUT(NK-1)+QLFXIN(NK)
    2230. 

  2230.               QIFXOUT(NK-1)=QIFXOUT(NK-1)+QIFXIN(NK)
    2231. 

  2231.               QRFXOUT(NK-1)=QRFXOUT(NK-1)+QRFXIN(NK)
    2232. 

  2232.               QSFXOUT(NK-1)=QSFXOUT(NK-1)+QSFXIN(NK)
    2233. 

  2233.             ELSE
    2234. 

  2234.               QLFXOUT(NK)=FXM(NK)*QLPA(NK)
    2235. 

  2235.               QIFXOUT(NK)=FXM(NK)*QIPA(NK)
    2236. 

  2236.               QRFXOUT(NK)=FXM(NK)*QRPA(NK)
    2237. 

  2237.               QSFXOUT(NK)=FXM(NK)*QSPA(NK)
    2238. 

  2238.               QLFXIN(NK-1)=QLFXIN(NK-1)+QLFXOUT(NK)
    2239. 

  2239.               QIFXIN(NK-1)=QIFXIN(NK-1)+QIFXOUT(NK)
    2240. 

  2240.               QRFXIN(NK-1)=QRFXIN(NK-1)+QRFXOUT(NK)
    2241. 

  2241.               QSFXIN(NK-1)=QSFXIN(NK-1)+QSFXOUT(NK)
    2242. 

  2242.             ENDIF
    2243. 

  2243.           ENDDO   
    2244. 

  2244. !     
    2245. 

  2245. !...UPDATE THE HYDROMETEOR CONCENTRATION VALUES AT EACH LEVEL...
    2246. 

  2246. !     
    2247. 

  2247.           DO NK=1,LTOP
    2248. 

  2248.             QLPA(NK)=QLPA(NK)+(QLFXIN(NK)+DETLQ(NK)-QLFXOUT(NK))*DTIME*EMSD(NK)
    2249. 

  2249.             QIPA(NK)=QIPA(NK)+(QIFXIN(NK)+DETIC(NK)-QIFXOUT(NK))*DTIME*EMSD(NK)
    2250. 

  2250.             QRPA(NK)=QRPA(NK)+(QRFXIN(NK)-QRFXOUT(NK)+RAINFB(NK))*DTIME*EMSD(NK)         !  PPT FB MODS
    2251. 

  2251.             QSPA(NK)=QSPA(NK)+(QSFXIN(NK)-QSFXOUT(NK)+SNOWFB(NK))*DTIME*EMSD(NK)         !  PPT FB MODS
    2252. 

  2252.           ENDDO     
    2253. 

  2253.         ENDDO
    2254. 

  2254.         DO NK=1,LTOP
    2255. 

  2255.           QLG(NK)=QLPA(NK)
    2256. 

  2256.           QIG(NK)=QIPA(NK)
    2257. 

  2257.           QRG(NK)=QRPA(NK)
    2258. 

  2258.           QSG(NK)=QSPA(NK)
    2259. 

  2259.         ENDDO   
    2260. 

  2260. !
    2261. 

  2261. !...CLEAN THINGS UP, CALCULATE CONVECTIVE FEEDBACK TENDENCIES FOR THIS
    2262. 

  2262. !...GRID POINT...
    2263. 

  2263. !     
    2264. 

  2264. !     IF (XTIME.LT.10.)THEN
    2265. 

  2265. !     WRITE(98,1080)LFS,LDB,LDT,TIMEC,TADVEC,NSTEP,NCOUNT,FABE,AINC 
    2266. 

  2266. !     ENDIF
    2267. 

  2267.        IF(IPRNT)THEN  
    2268. 

  2268. !         WRITE(98,1080)LFS,LDB,LDT,TIMEC,TADVEC,NSTEP,NCOUNT,FABE,AINC
    2269. 

  2269.          WRITE(message,1080)LFS,LDB,LDT,TIMEC,TADVEC,NSTEP,NCOUNT,FABE,AINC
    2270. 

  2270.          CALL wrf_message(message)
    2271. 

  2271. !        call flush(98)   
    2272. 

  2272.        endif  
    2273. 

  2273. !     
    2274. 

  2274. !...SEND FINAL PARAMETERIZED VALUES TO OUTPUT FILES...
    2275. 

  2275. !     
    2276. 

  2276. !297   IF(IPRNT)then 
    2277. 

  2277.        IF(IPRNT)then 
    2278. 

  2278. !    if(I.eq.16 .and. J.eq.41)then
    2279. 

  2279. !      IF(ISTOP.EQ.1)THEN
    2280. 

  2280.          write(98,*)
    2281. 

  2281. !        write(98,*)'At t(h), I, J =',float(NTSD)*72./3600.,I,J
    2282. 

  2282.          write(message,*)'P(LC), DTP, WKL, WKLCL =',p0(LC)/100.,       &
    2283. 

  2283.                      TLCL+DTLCL+dtrh-TENV,WKL,WKLCL
    2284. 

  2284.          call wrf_message(message)
    2285. 

  2285.          write(message,*)'TLCL, DTLCL, DTRH, TENV =',TLCL,DTLCL,       &
    2286. 

  2286.                       DTRH,TENV   
    2287. 

  2287.          call wrf_message(message)
    2288. 

  2288.          WRITE(message,1025)KLCL,ZLCL,DTLCL,LTOP,P0(LTOP),IFLAG,       &
    2289. 

  2289.          TMIX-T00,PMIX,QMIX,ABE
    2290. 

  2290.          call wrf_message(message)
    2291. 

  2291.          WRITE(message,1030)P0(LET)/100.,P0(LTOP)/100.,VMFLCL,PLCL/100.,  &
    2292. 

  2292.          WLCL,CLDHGT(LC)
    2293. 

  2293.          call wrf_message(message)
    2294. 

  2294.          WRITE(message,1035)PEF,PEFCBH,LC,LET,WKL,VWS 
    2295. 

  2295.          call wrf_message(message)
    2296. 

  2296.          write(message,*)'PRECIP EFFICIENCY =',PEFF 
    2297. 

  2297.          call wrf_message(message)
    2298. 

  2298.       WRITE(message,1080)LFS,LDB,LDT,TIMEC,TADVEC,NSTEP,NCOUNT,FABE,AINC
    2299. 

  2299.          call wrf_message(message)
    2300. 

  2300. !      ENDIF
    2301. 

  2301. !!!!! HERE !!!!!!!
    2302. 

  2302.            WRITE(message,1070)'  P  ','   DP ',' DT K/D ',' DR K/D ','   OMG  ',        &
    2303. 

  2303.           ' DOMGDP ','   UMF  ','   UER  ','   UDR  ','   DMF  ','   DER  '        &
    2304. 

  2304.           ,'   DDR  ','   EMS  ','    W0  ','  DETLQ ',' DETIC '
    2305. 

  2305.          call wrf_message(message)
    2306. 

  2306.            write(message,*)'just before DO 300...'
    2307. 

  2307.          call wrf_message(message)
    2308. 

  2308. !          call flush(98)
    2309. 

  2309.            DO NK=1,LTOP
    2310. 

  2310.              K=LTOP-NK+1
    2311. 

  2311.              DTT=(TG(K)-T0(K))*86400./TIMEC
    2312. 

  2312.              RL=XLV0-XLV1*TG(K)
    2313. 

  2313.              DR=-(QG(K)-Q0(K))*RL*86400./(TIMEC*CP)
    2314. 

  2314.              UDFRC=UDR(K)*TIMEC*EMSD(K)
    2315. 

  2315.              UEFRC=UER(K)*TIMEC*EMSD(K)
    2316. 

  2316.              DDFRC=DDR(K)*TIMEC*EMSD(K)
    2317. 

  2317.              DEFRC=-DER(K)*TIMEC*EMSD(K)
    2318. 

  2318.              WRITE(message,1075)P0(K)/100.,DP(K)/100.,DTT,DR,OMG(K),DOMGDP(K)*1.E4,       &
    2319. 

  2319.              UMF(K)/1.E6,UEFRC,UDFRC,DMF(K)/1.E6,DEFRC,DDFRC,EMS(K)/1.E11,           &
    2320. 

  2320.              W0AVG1D(K)*1.E2,DETLQ(K)*TIMEC*EMSD(K)*1.E3,DETIC(K)*                   &
    2321. 

  2321.              TIMEC*EMSD(K)*1.E3
    2322. 

  2322.          call wrf_message(message)
    2323. 

  2323.            ENDDO
    2324. 

  2324.            WRITE(message,1085)'K','P','Z','T0','TG','DT','TU','TD','Q0','QG',             &
    2325. 

  2325.                   'DQ','QU','QD','QLG','QIG','QRG','QSG','RH0','RHG'
    2326. 

  2326.          call wrf_message(message)
    2327. 

  2327.            DO NK=1,KL
    2328. 

  2328.              K=KX-NK+1
    2329. 

  2329.              DTT=TG(K)-T0(K)
    2330. 

  2330.              TUC=TU(K)-T00
    2331. 

  2331.              IF(K.LT.LC.OR.K.GT.LTOP)TUC=0.
    2332. 

  2332.              TDC=TZ(K)-T00
    2333. 

  2333.              IF((K.LT.LDB.OR.K.GT.LDT).AND.K.NE.LFS)TDC=0.
    2334. 

  2334.              IF(T0(K).LT.T00)THEN
    2335. 

  2335.                ES=ALIQ*EXP((BLIQ*TG(K)-CLIQ)/(TG(K)-DLIQ))
    2336. 

  2336.              ELSE
    2337. 

  2337.                ES=ALIQ*EXP((BLIQ*TG(K)-CLIQ)/(TG(K)-DLIQ))
    2338. 

  2338.              ENDIF  
    2339. 

  2339.              QGS=ES*0.622/(P0(K)-ES)
    2340. 

  2340.              RH0=Q0(K)/QES(K)
    2341. 

  2341.              RHG=QG(K)/QGS
    2342. 

  2342.              WRITE(message,1090)K,P0(K)/100.,Z0(K),T0(K)-T00,TG(K)-T00,DTT,TUC,            &
    2343. 

  2343.              TDC,Q0(K)*1000.,QG(K)*1000.,(QG(K)-Q0(K))*1000.,QU(K)*                   &
    2344. 

  2344.              1000.,QD(K)*1000.,QLG(K)*1000.,QIG(K)*1000.,QRG(K)*1000.,                &
    2345. 

  2345.              QSG(K)*1000.,RH0,RHG
    2346. 

  2346.          call wrf_message(message)
    2347. 

  2347.            ENDDO
    2348. 

  2348. !     
    2349. 

  2349. !...IF CALCULATIONS ABOVE SHOW AN ERROR IN THE MASS BUDGET, PRINT OUT A
    2350. 

  2350. !...TO BE USED LATER FOR DIAGNOSTIC PURPOSES, THEN ABORT RUN...
    2351. 

  2351. !     
    2352. 

  2352. !         IF(ISTOP.EQ.1 .or. ISHALL.EQ.1)THEN
    2353. 

  2353. 

  2354. !         IF(ISHALL.NE.1)THEN
    2355. 

  2355. !            write(98,4421)i,j,iyr,imo,idy,ihr,imn
    2356. 

  2356. !           write(98)i,j,iyr,imo,idy,ihr,imn,kl
    2357. 

  2357. ! 4421       format(7i4)
    2358. 

  2358. !            write(98,4422)kl
    2359. 

  2359. ! 4422       format(i6) 
    2360. 

  2360.             DO 310 NK = 1,KL
    2361. 

  2361.               k = kl - nk + 1
    2362. 

  2362.               write(98,4455) p0(k)/100.,t0(k)-273.16,q0(k)*1000.,       &
    2363. 

  2363.                        u0(k),v0(k),W0AVG1D(K),dp(k),tke(k)
    2364. 

  2364. !             write(98) p0,t0,q0,u0,v0,w0,dp,tke
    2365. 

  2365. !           WRITE(98,1115)Z0(K),P0(K)/100.,T0(K)-273.16,Q0(K)*1000.,
    2366. 

  2366. !    *               U0(K),V0(K),DP(K)/100.,W0AVG(I,J,K)
    2367. 

  2367.  310        CONTINUE
    2368. 

  2368.             IF(ISTOP.EQ.1)THEN
    2369. 

  2369.               CALL wrf_error_fatal ( 'KAIN-FRITSCH, istop=1, diags' )
    2370. 

  2370.             ENDIF
    2371. 

  2371. !         ENDIF
    2372. 

  2372.   4455  format(8f11.3) 
    2373. 

  2373.        ENDIF
    2374. 

  2374.         CNDTNF=(1.-EQFRC(LFS))*(QLIQ(LFS)+QICE(LFS))*DMF(LFS)
    2375. 

  2375.         PRATEC(I,J)=PPTFLX*(1.-FBFRC)/DXSQ
    2376. 

  2376.         RAINCV(I,J)=DT*PRATEC(I,J)     !  PPT FB MODS
    2377. 

  2377. !        RAINCV(I,J)=.1*.5*DT*PPTFLX/DXSQ               !  PPT FB MODS
    2378. 

  2378. !         RNC=0.1*TIMEC*PPTFLX/DXSQ
    2379. 

  2379.         RNC=RAINCV(I,J)*NIC
    2380. 

  2380.        IF(ISHALL.EQ.0.AND.IPRNT)write (98,909)I,J,RNC
    2381. 

  2381. 

  2382. !     WRITE(98,1095)CPR*AINC,TDER+PPTFLX+CNDTNF
    2383. 

  2383. !     
    2384. 

  2384. !  EVALUATE MOISTURE BUDGET...
    2385. 

  2385. !     
    2386. 

  2386. 

  2387.         QINIT=0.
    2388. 

  2388.         QFNL=0.
    2389. 

  2389.         DPT=0.
    2390. 

  2390.         DO 315 NK=1,LTOP
    2391. 

  2391.           DPT=DPT+DP(NK)
    2392. 

  2392.           QINIT=QINIT+Q0(NK)*EMS(NK)
    2393. 

  2393.           QFNL=QFNL+QG(NK)*EMS(NK)
    2394. 

  2394.           QFNL=QFNL+(QLG(NK)+QIG(NK)+QRG(NK)+QSG(NK))*EMS(NK)
    2395. 

  2395.   315   CONTINUE
    2396. 

  2396.         QFNL=QFNL+PPTFLX*TIMEC*(1.-FBFRC)       !  PPT FB MODS
    2397. 

  2397. !        QFNL=QFNL+PPTFLX*TIMEC                 !  PPT FB MODS
    2398. 

  2398.         ERR2=(QFNL-QINIT)*100./QINIT
    2399. 

  2399.        IF(IPRNT)WRITE(98,1110)QINIT,QFNL,ERR2
    2400. 

  2400.       IF(ABS(ERR2).GT.0.05 .AND. ISTOP.EQ.0)THEN 
    2401. 

  2401. !       write(99,*)'!!!!!!!! MOISTURE BUDGET ERROR IN KFPARA !!!'
    2402. 

  2402. !       WRITE(99,1110)QINIT,QFNL,ERR2
    2403. 

  2403.         IPRNT=.TRUE.
    2404. 

  2404.         ISTOP=1
    2405. 

  2405.             write(98,4422)kl
    2406. 

  2406.  4422       format(i6)
    2407. 

  2407.             DO 311 NK = 1,KL
    2408. 

  2408.               k = kl - nk + 1
    2409. 

  2409. !             write(99,4455) p0(k)/100.,t0(k)-273.16,q0(k)*1000.,       &
    2410. 

  2410. !                      u0(k),v0(k),W0AVG1D(K),dp(k)
    2411. 

  2411. !             write(98) p0,t0,q0,u0,v0,w0,dp,tke
    2412. 

  2412. !           WRITE(98,1115)P0(K)/100.,T0(K)-273.16,Q0(K)*1000.,          &
    2413. 

  2413. !                    U0(K),V0(K),W0AVG1D(K),dp(k)/100.,tke(k)
    2414. 

  2414.             WRITE(98,4456)P0(K)/100.,T0(K)-273.16,Q0(K)*1000.,          &
    2415. 

  2415.                      U0(K),V0(K),W0AVG1D(K),dp(k)/100.,tke(k)
    2416. 

  2416.  311        CONTINUE
    2417. 

  2417. !           call flush(98)
    2418. 

  2418. 

  2419. !        GOTO 297
    2420. 

  2420. !         STOP 'QVERR'
    2421. 

  2421.       ENDIF
    2422. 

  2422.  1115 FORMAT (2X,F7.2,2X,F5.1,2X,F6.3,2(2X,F5.1),2X,F7.2,2X,F7.4)
    2423. 

  2423.  4456  format(8f12.3)
    2424. 

  2424.         IF(PPTFLX.GT.0.)THEN
    2425. 

  2425.           RELERR=ERR2*QINIT/(PPTFLX*TIMEC)
    2426. 

  2426.         ELSE
    2427. 

  2427.           RELERR=0.
    2428. 

  2428.         ENDIF
    2429. 

  2429.      IF(IPRNT)THEN
    2430. 

  2430.         WRITE(98,1120)RELERR
    2431. 

  2431.         WRITE(98,*)'TDER, CPR, TRPPT =',              &
    2432. 

  2432.           TDER,CPR*AINC,TRPPT*AINC
    2433. 

  2433.      ENDIF
    2434. 

  2434. !     
    2435. 

  2435. !...FEEDBACK TO RESOLVABLE SCALE TENDENCIES.
    2436. 

  2436. !     
    2437. 

  2437. !...IF THE ADVECTIVE TIME PERIOD (TADVEC) IS LESS THAN SPECIFIED MINIMUM
    2438. 

  2438. !...TIMEC, ALLOW FEEDBACK TO OCCUR ONLY DURING TADVEC...
    2439. 

  2439. !     
    2440. 

  2440.         IF(TADVEC.LT.TIMEC)NIC=NINT(TADVEC/DT)
    2441. 

  2441.         NCA(I,J) = REAL(NIC)*DT
    2442. 

  2442.         IF(ISHALL.EQ.1)THEN
    2443. 

  2443.            TIMEC = 2400.
    2444. 

  2444.            NCA(I,J) = CUDT*60.
    2445. 

  2445.            NSHALL = NSHALL+1
    2446. 

  2446.         ENDIF
    2447. 

  2447. 

  2448.         DO K=1,KX
    2449. 

  2449. !         IF(IMOIST(INEST).NE.2)THEN
    2450. 

  2450. !
    2451. 

  2451. !...IF HYDROMETEORS ARE NOT ALLOWED, THEY MUST BE EVAPORATED OR SUBLIMATED
    2452. 

  2452. !...AND FED BACK AS VAPOR, ALONG WITH ASSOCIATED CHANGES IN TEMPERATURE.
    2453. 

  2453. !...NOTE:  THIS WILL INTRODUCE CHANGES IN THE CONVECTIVE TEMPERATURE AND
    2454. 

  2454. !...WATER VAPOR FEEDBACK TENDENCIES AND MAY LEAD TO SUPERSATURATED VALUE
    2455. 

  2455. !...OF QG...
    2456. 

  2456. !
    2457. 

  2457. !           RLC=XLV0-XLV1*TG(K)
    2458. 

  2458. !           RLS=XLS0-XLS1*TG(K)
    2459. 

  2459. !           CPM=CP*(1.+0.887*QG(K))
    2460. 

  2460. !           TG(K)=TG(K)-(RLC*(QLG(K)+QRG(K))+RLS*(QIG(K)+QSG(K)))/CPM
    2461. 

  2461. !           QG(K)=QG(K)+(QLG(K)+QRG(K)+QIG(K)+QSG(K))
    2462. 

  2462. !           DQLDT(I,J,NK)=0.
    2463. 

  2463. !           DQIDT(I,J,NK)=0.
    2464. 

  2464. !           DQRDT(I,J,NK)=0.
    2465. 

  2465. !           DQSDT(I,J,NK)=0.
    2466. 

  2466. !         ELSE
    2467. 

  2467. !
    2468. 

  2468. !...IF ICE PHASE IS NOT ALLOWED, MELT ALL FROZEN HYDROMETEORS...
    2469. 

  2469. !
    2470. 

  2470.           IF(.NOT. F_QI .and. warm_rain)THEN
    2471. 

  2471. 

  2472.             CPM=CP*(1.+0.887*QG(K))
    2473. 

  2473.             TG(K)=TG(K)-(QIG(K)+QSG(K))*RLF/CPM
    2474. 

  2474.             DQCDT(K)=(QLG(K)+QIG(K)-QL0(K)-QI0(K))/TIMEC
    2475. 

  2475.             DQIDT(K)=0.
    2476. 

  2476.             DQRDT(K)=(QRG(K)+QSG(K)-QR0(K)-QS0(K))/TIMEC
    2477. 

  2477.             DQSDT(K)=0.
    2478. 

  2478.           ELSEIF(.NOT. F_QI .and. .not. warm_rain)THEN
    2479. 

  2479. !
    2480. 

  2480. !...IF ICE PHASE IS ALLOWED, BUT MIXED PHASE IS NOT, MELT FROZEN HYDROMETEORS
    2481. 

  2481. !...BELOW THE MELTING LEVEL, FREEZE LIQUID WATER ABOVE THE MELTING LEVEL
    2482. 

  2482. !
    2483. 

  2483.             CPM=CP*(1.+0.887*QG(K))
    2484. 

  2484.             IF(K.LE.ML)THEN
    2485. 

  2485.               TG(K)=TG(K)-(QIG(K)+QSG(K))*RLF/CPM
    2486. 

  2486.             ELSEIF(K.GT.ML)THEN
    2487. 

  2487.               TG(K)=TG(K)+(QLG(K)+QRG(K))*RLF/CPM
    2488. 

  2488.             ENDIF
    2489. 

  2489.             DQCDT(K)=(QLG(K)+QIG(K)-QL0(K)-QI0(K))/TIMEC
    2490. 

  2490.             DQIDT(K)=0.
    2491. 

  2491.             DQRDT(K)=(QRG(K)+QSG(K)-QR0(K)-QS0(K))/TIMEC
    2492. 

  2492.             DQSDT(K)=0.
    2493. 

  2493.           ELSEIF(F_QI) THEN
    2494. 

  2494. !
    2495. 

  2495. !...IF MIXED PHASE HYDROMETEORS ARE ALLOWED, FEED BACK CONVECTIVE TENDENCIES
    2496. 

  2496. !...OF HYDROMETEORS DIRECTLY...
    2497. 

  2497. !
    2498. 

  2498.             DQCDT(K)=(QLG(K)-QL0(K))/TIMEC
    2499. 

  2499.             DQIDT(K)=(QIG(K)-QI0(K))/TIMEC
    2500. 

  2500.             DQRDT(K)=(QRG(K)-QR0(K))/TIMEC
    2501. 

  2501.             IF (F_QS) THEN
    2502. 

  2502.                DQSDT(K)=(QSG(K)-QS0(K))/TIMEC
    2503. 

  2503.             ELSE
    2504. 

  2504.                DQIDT(K)=DQIDT(K)+(QSG(K)-QS0(K))/TIMEC
    2505. 

  2505.             ENDIF
    2506. 

  2506.           ELSE
    2507. 

  2507. !              PRINT *,'THIS COMBINATION OF IMOIST, IEXICE, IICE NOT ALLOWED!'
    2508. 

  2508.               CALL wrf_error_fatal ( 'KAIN-FRITSCH, THIS COMBINATION OF IMOIST, IEXICE, IICE NOT ALLOWED' )
    2509. 

  2509.           ENDIF
    2510. 

  2510.           DTDT(K)=(TG(K)-T0(K))/TIMEC
    2511. 

  2511.           DQDT(K)=(QG(K)-Q0(K))/TIMEC
    2512. 

  2512.         ENDDO
    2513. 

  2513.         PRATEC(I,J)=PPTFLX*(1.-FBFRC)/DXSQ
    2514. 

  2514.         RAINCV(I,J)=DT*PRATEC(I,J)
    2515. 

  2515. !        RAINCV(I,J)=.1*.5*DT*PPTFLX/DXSQ               !  PPT FB MODS
    2516. 

  2516. !         RNC=0.1*TIMEC*PPTFLX/DXSQ
    2517. 

  2517.         RNC=RAINCV(I,J)*NIC
    2518. 

  2518.  909     FORMAT('AT I, J =',i3,1x,i3,' CONVECTIVE RAINFALL =',F8.4,' mm')
    2519. 

  2519. !      write (98,909)I,J,RNC
    2520. 

  2520. !      write (6,909)I,J,RNC
    2521. 

  2521. !      WRITE(98,*)'at NTSD =',NTSD,',No. of KF points activated =',
    2522. 

  2522. !     *            NCCNT
    2523. 

  2523. !      call flush(98)
    2524. 

  2524. 1000  FORMAT(' ',10A8)
    2525. 

  2525. 1005  FORMAT(' ',F6.0,2X,F6.4,2X,F7.3,1X,F6.4,2X,4(F6.3,2X),2(F7.3,1X))
    2526. 

  2526. 1010  FORMAT(' ',' VERTICAL VELOCITY IS NEGATIVE AT ',F4.0,' MB')
    2527. 

  2527. 1015   FORMAT(' ','ALL REMAINING MASS DETRAINS BELOW ',F4.0,' MB')
    2528. 

  2528. 1025   FORMAT(5X,' KLCL=',I2,' ZLCL=',F7.1,'M',                         &
    2529. 

  2529.         ' DTLCL=',F5.2,' LTOP=',I2,' P0(LTOP)=',-2PF5.1,'MB FRZ LV=',   &
    2530. 

  2530.         I2,' TMIX=',0PF4.1,1X,'PMIX=',-2PF6.1,' QMIX=',3PF5.1,          &
    2531. 

  2531.         ' CAPE=',0PF7.1)
    2532. 

  2532. 1030   FORMAT(' ',' P0(LET) = ',F6.1,' P0(LTOP) = ',F6.1,' VMFLCL =',   &
    2533. 

  2533.       E12.3,' PLCL =',F6.1,' WLCL =',F6.3,' CLDHGT =',                  &
    2534. 

  2534.       F8.1)
    2535. 

  2535. 1035  FORMAT(1X,'PEF(WS)=',F4.2,'(CB)=',F4.2,'LC,LET=',2I3,'WKL='       &
    2536. 

  2536.       ,F6.3,'VWS=',F5.2)
    2537. 

  2537. !1055  FORMAT('*** DEGREE OF STABILIZATION =',F5.3,                  &
    2538. 

  2538. !      ', NO MORE MASS FLUX IS ALLOWED!')
    2539. 

  2539. !1060     FORMAT(' ITERATION DOES NOT CONVERGE TO GIVE THE SPECIFIED    &
    2540. 

  2540. !      &DEGREE OF STABILIZATION!  FABE= ',F6.4) 
    2541. 

  2541.  1070 FORMAT (16A8) 
    2542. 

  2542.  1075 FORMAT (F8.2,3(F8.2),2(F8.3),F8.2,2F8.3,F8.2,6F8.3) 
    2543. 

  2543.  1080 FORMAT(2X,'LFS,LDB,LDT =',3I3,' TIMEC, TADVEC, NSTEP=',           &
    2544. 

  2544.               2(1X,F5.0),I3,'NCOUNT, FABE, AINC=',I2,1X,F5.3,F6.2) 
    2545. 

  2545.  1085 FORMAT (A3,16A7,2A8) 
    2546. 

  2546.  1090 FORMAT (I3,F7.2,F7.0,10F7.2,4F7.3,2F8.3) 
    2547. 

  2547.  1095 FORMAT(' ','  PPT PRODUCTION RATE= ',F10.0,' TOTAL EVAP+PPT= ',F10.0)
    2548. 

  2548. 1105   FORMAT(' ','NET LATENT HEAT RELEASE =',E12.5,' ACTUAL HEATING =',&
    2549. 

  2549.        E12.5,' J/KG-S, DIFFERENCE = ',F9.3,'%')
    2550. 

  2550. 1110   FORMAT(' ','INITIAL WATER =',E12.5,' FINAL WATER =',E12.5,       &
    2551. 

  2551.        ' TOTAL WATER CHANGE =',F8.2,'%')
    2552. 

  2552. ! 1115 FORMAT (2X,F6.0,2X,F7.2,2X,F5.1,2X,F6.3,2(2X,F5.1),2X,F7.2,2X,F7.4)
    2553. 

  2553. 1120   FORMAT(' ','MOISTURE ERROR AS FUNCTION OF TOTAL PPT =',F9.3,'%')
    2554. 

  2554. !
    2555. 

  2555. !-----------------------------------------------------------------------
    2556. 

  2556. !--------------SAVE CLOUD TOP AND BOTTOM FOR RADIATION------------------
    2557. 

  2557. !-----------------------------------------------------------------------
    2558. 

  2558. !
    2559. 

  2559.       CUTOP(I,J)=REAL(LTOP)
    2560. 

  2560.       CUBOT(I,J)=REAL(LCL)
    2561. 

  2561. !
    2562. 

  2562. !-----------------------------------------------------------------------
    2563. 

  2563.    END SUBROUTINE  KF_eta_PARA
    2564. 

  2564. !********************************************************************
    2565. 

  2565. ! ***********************************************************************
    2566. 

  2566.    SUBROUTINE TPMIX2(p,thes,tu,qu,qliq,qice,qnewlq,qnewic,XLV1,XLV0)
    2567. 

  2567. !
    2568. 

  2568. ! Lookup table variables:
    2569. 

  2569. !     INTEGER, PARAMETER :: (KFNT=250,KFNP=220)
    2570. 

  2570. !     REAL, SAVE, DIMENSION(1:KFNT,1:KFNP) :: TTAB,QSTAB
    2571. 

  2571. !     REAL, SAVE, DIMENSION(1:KFNP) :: THE0K
    2572. 

  2572. !     REAL, SAVE, DIMENSION(1:200) :: ALU
    2573. 

  2573. !     REAL, SAVE :: RDPR,RDTHK,PLUTOP
    2574. 

  2574. ! End of Lookup table variables:
    2575. 

  2575. !-----------------------------------------------------------------------
    2576. 

  2576.    IMPLICIT NONE
    2577. 

  2577. !-----------------------------------------------------------------------
    2578. 

  2578.    REAL,         INTENT(IN   )   :: P,THES,XLV1,XLV0
    2579. 

  2579.    REAL,         INTENT(OUT  )   :: QNEWLQ,QNEWIC
    2580. 

  2580.    REAL,         INTENT(INOUT)   :: TU,QU,QLIQ,QICE
    2581. 

  2581.    REAL    ::    TP,QQ,BTH,TTH,PP,T00,T10,T01,T11,Q00,Q10,Q01,Q11,          &
    2582. 

  2582.                  TEMP,QS,QNEW,DQ,QTOT,RLL,CPP
    2583. 

  2583.    INTEGER ::    IPTB,ITHTB
    2584. 

  2584. !-----------------------------------------------------------------------
    2585. 

  2585. 

  2586. !c******** LOOKUP TABLE VARIABLES... ****************************
    2587. 

  2587. !      parameter(kfnt=250,kfnp=220)
    2588. 

  2588. !c
    2589. 

  2589. !      COMMON/KFLUT/ ttab(kfnt,kfnp),qstab(kfnt,kfnp),the0k(kfnp),
    2590. 

  2590. !     *              alu(200),rdpr,rdthk,plutop 
    2591. 

  2591. !C*************************************************************** 
    2592. 

  2592. !c
    2593. 

  2593. !c***********************************************************************
    2594. 

  2594. !c     scaling pressure and tt table index                         
    2595. 

  2595. !c***********************************************************************
    2596. 

  2596. !c
    2597. 

  2597.       tp=(p-plutop)*rdpr
    2598. 

  2598.       qq=tp-aint(tp)
    2599. 

  2599.       iptb=int(tp)+1
    2600. 

  2600. 

  2601. !
    2602. 

  2602. !***********************************************************************
    2603. 

  2603. !              base and scaling factor for the                           
    2604. 

  2604. !***********************************************************************
    2605. 

  2605. !
    2606. 

  2606. !  scaling the and tt table index                                        
    2607. 

  2607.       bth=(the0k(iptb+1)-the0k(iptb))*qq+the0k(iptb)
    2608. 

  2608.       tth=(thes-bth)*rdthk
    2609. 

  2609.       pp   =tth-aint(tth)
    2610. 

  2610.       ithtb=int(tth)+1
    2611. 

  2611.        IF(IPTB.GE.220 .OR. IPTB.LE.1 .OR. ITHTB.GE.250 .OR. ITHTB.LE.1)THEN
    2612. 

  2612.          write(98,*)'**** OUT OF BOUNDS *********'
    2613. 

  2613. !        call flush(98)
    2614. 

  2614.        ENDIF
    2615. 

  2615. !
    2616. 

  2616.       t00=ttab(ithtb  ,iptb  )
    2617. 

  2617.       t10=ttab(ithtb+1,iptb  )
    2618. 

  2618.       t01=ttab(ithtb  ,iptb+1)
    2619. 

  2619.       t11=ttab(ithtb+1,iptb+1)
    2620. 

  2620. !
    2621. 

  2621.       q00=qstab(ithtb  ,iptb  )
    2622. 

  2622.       q10=qstab(ithtb+1,iptb  )
    2623. 

  2623.       q01=qstab(ithtb  ,iptb+1)
    2624. 

  2624.       q11=qstab(ithtb+1,iptb+1)
    2625. 

  2625. !
    2626. 

  2626. !***********************************************************************
    2627. 

  2627. !              parcel temperature                                        
    2628. 

  2628. !***********************************************************************
    2629. 

  2629. !
    2630. 

  2630.       temp=(t00+(t10-t00)*pp+(t01-t00)*qq+(t00-t10-t01+t11)*pp*qq)
    2631. 

  2631. !
    2632. 

  2632.       qs=(q00+(q10-q00)*pp+(q01-q00)*qq+(q00-q10-q01+q11)*pp*qq)
    2633. 

  2633. !
    2634. 

  2634.       DQ=QS-QU
    2635. 

  2635.       IF(DQ.LE.0.)THEN
    2636. 

  2636.         QNEW=QU-QS
    2637. 

  2637.         QU=QS
    2638. 

  2638.       ELSE 
    2639. 

  2639. !
    2640. 

  2640. !   IF THE PARCEL IS SUBSATURATED, TEMPERATURE AND MIXING RATIO MUST BE
    2641. 

  2641. !   ADJUSTED...IF LIQUID WATER IS PRESENT, IT IS ALLOWED TO EVAPORATE
    2642. 

  2642. ! 
    2643. 

  2643.         QNEW=0.
    2644. 

  2644.         QTOT=QLIQ+QICE
    2645. 

  2645. !
    2646. 

  2646. !   IF THERE IS ENOUGH LIQUID OR ICE TO SATURATE THE PARCEL, TEMP STAYS AT ITS
    2647. 

  2647. !   WET BULB VALUE, VAPOR MIXING RATIO IS AT SATURATED LEVEL, AND THE MIXING
    2648. 

  2648. !   RATIOS OF LIQUID AND ICE ARE ADJUSTED TO MAKE UP THE ORIGINAL SATURATION
    2649. 

  2649. !   DEFICIT... OTHERWISE, ANY AVAILABLE LIQ OR ICE VAPORIZES AND APPROPRIATE
    2650. 

  2650. !   ADJUSTMENTS TO PARCEL TEMP; VAPOR, LIQUID, AND ICE MIXING RATIOS ARE MADE.
    2651. 

  2651. !
    2652. 

  2652. !...subsaturated values only occur in calculations involving various mixtures of
    2653. 

  2653. !...updraft and environmental air for estimation of entrainment and detrainment.
    2654. 

  2654. !...For these purposes, assume that reasonable estimates can be given using 
    2655. 

  2655. !...liquid water saturation calculations only - i.e., ignore the effect of the
    2656. 

  2656. !...ice phase in this process only...will not affect conservative properties...
    2657. 

  2657. !
    2658. 

  2658.         IF(QTOT.GE.DQ)THEN
    2659. 

  2659.           qliq=qliq-dq*qliq/(qtot+1.e-10)
    2660. 

  2660.           qice=qice-dq*qice/(qtot+1.e-10)
    2661. 

  2661.           QU=QS
    2662. 

  2662.         ELSE
    2663. 

  2663.           RLL=XLV0-XLV1*TEMP
    2664. 

  2664.           CPP=1004.5*(1.+0.89*QU)
    2665. 

  2665.           IF(QTOT.LT.1.E-10)THEN
    2666. 

  2666. !
    2667. 

  2667. !...IF NO LIQUID WATER OR ICE IS AVAILABLE, TEMPERATURE IS GIVEN BY:
    2668. 

  2668.             TEMP=TEMP+RLL*(DQ/(1.+DQ))/CPP
    2669. 

  2669.           ELSE
    2670. 

  2670. !
    2671. 

  2671. !...IF SOME LIQ WATER/ICE IS AVAILABLE, BUT NOT ENOUGH TO ACHIEVE SATURATION,
    2672. 

  2672. !   THE TEMPERATURE IS GIVEN BY:
    2673. 

  2673. !
    2674. 

  2674.             TEMP=TEMP+RLL*((DQ-QTOT)/(1+DQ-QTOT))/CPP
    2675. 

  2675.             QU=QU+QTOT
    2676. 

  2676.             QTOT=0.
    2677. 

  2677.             QLIQ=0.
    2678. 

  2678.             QICE=0.
    2679. 

  2679.           ENDIF
    2680. 

  2680.         ENDIF
    2681. 

  2681.       ENDIF
    2682. 

  2682.       TU=TEMP
    2683. 

  2683.       qnewlq=qnew
    2684. 

  2684.       qnewic=0.
    2685. 

  2685. !
    2686. 

  2686.    END SUBROUTINE TPMIX2
    2687. 

  2687. !******************************************************************************
    2688. 

  2688.       SUBROUTINE DTFRZNEW(TU,P,THTEU,QU,QFRZ,QICE,ALIQ,BLIQ,CLIQ,DLIQ)
    2689. 

  2689. !-----------------------------------------------------------------------
    2690. 

  2690.    IMPLICIT NONE
    2691. 

  2691. !-----------------------------------------------------------------------
    2692. 

  2692.    REAL,         INTENT(IN   )   :: P,QFRZ,ALIQ,BLIQ,CLIQ,DLIQ
    2693. 

  2693.    REAL,         INTENT(INOUT)   :: TU,THTEU,QU,QICE
    2694. 

  2694.    REAL    ::    RLC,RLS,RLF,CPP,A,DTFRZ,ES,QS,DQEVAP,PII
    2695. 

  2695. !-----------------------------------------------------------------------
    2696. 

  2696. !
    2697. 

  2697. !...ALLOW THE FREEZING OF LIQUID WATER IN THE UPDRAFT TO PROCEED AS AN 
    2698. 

  2698. !...APPROXIMATELY LINEAR FUNCTION OF TEMPERATURE IN THE TEMPERATURE RANGE 
    2699. 

  2699. !...TTFRZ TO TBFRZ...
    2700. 

  2700. !...FOR COLDER TEMPERATURES, FREEZE ALL LIQUID WATER...
    2701. 

  2701. !...THERMODYNAMIC PROPERTIES ARE STILL CALCULATED WITH RESPECT TO LIQUID WATER
    2702. 

  2702. !...TO ALLOW THE USE OF LOOKUP TABLE TO EXTRACT TMP FROM THETAE...
    2703. 

  2703. !
    2704. 

  2704.       RLC=2.5E6-2369.276*(TU-273.16)
    2705. 

  2705.       RLS=2833922.-259.532*(TU-273.16)
    2706. 

  2706.       RLF=RLS-RLC
    2707. 

  2707.       CPP=1004.5*(1.+0.89*QU)
    2708. 

  2708. !
    2709. 

  2709. !  A = D(es)/DT IS THAT CALCULATED FROM BUCK (1981) EMPERICAL FORMULAS
    2710. 

  2710. !  FOR SATURATION VAPOR PRESSURE...
    2711. 

  2711. !
    2712. 

  2712.       A=(CLIQ-BLIQ*DLIQ)/((TU-DLIQ)*(TU-DLIQ))
    2713. 

  2713.       DTFRZ = RLF*QFRZ/(CPP+RLS*QU*A)
    2714. 

  2714.       TU = TU+DTFRZ
    2715. 

  2715.       
    2716. 

  2716.       ES = ALIQ*EXP((BLIQ*TU-CLIQ)/(TU-DLIQ))
    2717. 

  2717.       QS = ES*0.622/(P-ES)
    2718. 

  2718. !
    2719. 

  2719. !...FREEZING WARMS THE AIR AND IT BECOMES UNSATURATED...ASSUME THAT SOME OF THE 
    2720. 

  2720. !...LIQUID WATER THAT IS AVAILABLE FOR FREEZING EVAPORATES TO MAINTAIN SATURA-
    2721. 

  2721. !...TION...SINCE THIS WATER HAS ALREADY BEEN TRANSFERRED TO THE ICE CATEGORY,
    2722. 

  2722. !...SUBTRACT IT FROM ICE CONCENTRATION, THEN SET UPDRAFT MIXING RATIO AT THE NEW
    2723. 

  2723. !...TEMPERATURE TO THE SATURATION VALUE...
    2724. 

  2724. !
    2725. 

  2725.       DQEVAP = QS-QU
    2726. 

  2726.       QICE = QICE-DQEVAP
    2727. 

  2727.       QU = QU+DQEVAP
    2728. 

  2728.       PII=(1.E5/P)**(0.2854*(1.-0.28*QU))
    2729. 

  2729.       THTEU=TU*PII*EXP((3374.6525/TU-2.5403)*QU*(1.+0.81*QU))
    2730. 

  2730. !
    2731. 

  2731.    END SUBROUTINE DTFRZNEW
    2732. 

  2732. ! --------------------------------------------------------------------------------
    2733. 

  2733. 

  2734.       SUBROUTINE CONDLOAD(QLIQ,QICE,WTW,DZ,BOTERM,ENTERM,RATE,QNEWLQ,           &
    2735. 

  2735.                           QNEWIC,QLQOUT,QICOUT,G)
    2736. 

  2736. 

  2737. !-----------------------------------------------------------------------
    2738. 

  2738.    IMPLICIT NONE
    2739. 

  2739. !-----------------------------------------------------------------------
    2740. 

  2740. !  9/18/88...THIS PRECIPITATION FALLOUT SCHEME IS BASED ON THE SCHEME US
    2741. 

  2741. !  BY OGURA AND CHO (1973).  LIQUID WATER FALLOUT FROM A PARCEL IS CAL-
    2742. 

  2742. !  CULATED USING THE EQUATION DQ=-RATE*Q*DT, BUT TO SIMULATE A QUASI-
    2743. 

  2743. !  CONTINUOUS PROCESS, AND TO ELIMINATE A DEPENDENCY ON VERTICAL
    2744. 

  2744. !  RESOLUTION THIS IS EXPRESSED AS Q=Q*EXP(-RATE*DZ).
    2745. 

  2745. 

  2746.       REAL, INTENT(IN   )   :: G
    2747. 

  2747.       REAL, INTENT(IN   )   :: DZ,BOTERM,ENTERM,RATE
    2748. 

  2748.       REAL, INTENT(INOUT)   :: QLQOUT,QICOUT,WTW,QLIQ,QICE,QNEWLQ,QNEWIC
    2749. 

  2749.       REAL :: QTOT,QNEW,QEST,G1,WAVG,CONV,RATIO3,OLDQ,RATIO4,DQ,PPTDRG
    2750. 

  2750. !
    2751. 

  2751.       QTOT=QLIQ+QICE                                                    
    2752. 

  2752.       QNEW=QNEWLQ+QNEWIC                                                
    2753. 

  2753. !                                                                       
    2754. 

  2754. !  ESTIMATE THE VERTICAL VELOCITY SO THAT AN AVERAGE VERTICAL VELOCITY 
    2755. 

  2755. !  BE CALCULATED TO ESTIMATE THE TIME REQUIRED FOR ASCENT BETWEEN MODEL 
    2756. 

  2756. !  LEVELS...                                                            
    2757. 

  2757. !                                                                       
    2758. 

  2758.       QEST=0.5*(QTOT+QNEW)                                              
    2759. 

  2759.       G1=WTW+BOTERM-ENTERM-2.*G*DZ*QEST/1.5                             
    2760. 

  2760.       IF(G1.LT.0.0)G1=0.                                                
    2761. 

  2761.       WAVG=0.5*(SQRT(WTW)+SQRT(G1))                                      
    2762. 

  2762.       CONV=RATE*DZ/WAVG               ! KF90  Eq. 9
    2763. 

  2763. !                                                                       
    2764. 

  2764. !  RATIO3 IS THE FRACTION OF LIQUID WATER IN FRESH CONDENSATE, RATIO4 IS
    2765. 

  2765. !  THE FRACTION OF LIQUID WATER IN THE TOTAL AMOUNT OF CONDENSATE INVOLV
    2766. 

  2766. !  IN THE PRECIPITATION PROCESS - NOTE THAT ONLY 60% OF THE FRESH CONDEN
    2767. 

  2767. !  SATE IS IS ALLOWED TO PARTICIPATE IN THE CONVERSION PROCESS...       
    2768. 

  2768. !                                                                       
    2769. 

  2769.       RATIO3=QNEWLQ/(QNEW+1.E-8)                                       
    2770. 

  2770. !     OLDQ=QTOT                                                         
    2771. 

  2771.       QTOT=QTOT+0.6*QNEW                                                
    2772. 

  2772.       OLDQ=QTOT                                                         
    2773. 

  2773.       RATIO4=(0.6*QNEWLQ+QLIQ)/(QTOT+1.E-8)                            
    2774. 

  2774.       QTOT=QTOT*EXP(-CONV)            ! KF90  Eq. 9                                              
    2775. 

  2775. !                                                                       
    2776. 

  2776. !  DETERMINE THE AMOUNT OF PRECIPITATION THAT FALLS OUT OF THE UPDRAFT  
    2777. 

  2777. !  PARCEL AT THIS LEVEL...                                              
    2778. 

  2778. !                                                                       
    2779. 

  2779.       DQ=OLDQ-QTOT                                                      
    2780. 

  2780.       QLQOUT=RATIO4*DQ                                                  
    2781. 

  2781.       QICOUT=(1.-RATIO4)*DQ                                             
    2782. 

  2782. !                                                                       
    2783. 

  2783. !  ESTIMATE THE MEAN LOAD OF CONDENSATE ON THE UPDRAFT IN THE LAYER, CAL
    2784. 

  2784. !  LATE VERTICAL VELOCITY                                               
    2785. 

  2785. !                                                                       
    2786. 

  2786.       PPTDRG=0.5*(OLDQ+QTOT-0.2*QNEW)                                   
    2787. 

  2787.       WTW=WTW+BOTERM-ENTERM-2.*G*DZ*PPTDRG/1.5                          
    2788. 

  2788.       IF(ABS(WTW).LT.1.E-4)WTW=1.E-4
    2789. 

  2789. !                                                                       
    2790. 

  2790. !  DETERMINE THE NEW LIQUID WATER AND ICE CONCENTRATIONS INCLUDING LOSSE
    2791. 

  2791. !  DUE TO PRECIPITATION AND GAINS FROM CONDENSATION...                  
    2792. 

  2792. !                                                                       
    2793. 

  2793.       QLIQ=RATIO4*QTOT+RATIO3*0.4*QNEW                                  
    2794. 

  2794.       QICE=(1.-RATIO4)*QTOT+(1.-RATIO3)*0.4*QNEW                        
    2795. 

  2795.       QNEWLQ=0.                                                         
    2796. 

  2796.       QNEWIC=0.                                                         
    2797. 

  2797. 

  2798.    END SUBROUTINE CONDLOAD
    2799. 

  2799. 

  2800. ! ----------------------------------------------------------------------
    2801. 

  2801.    SUBROUTINE PROF5(EQ,EE,UD)                                        
    2802. 

  2802. !
    2803. 

  2803. !***********************************************************************
    2804. 

  2804. !*****    GAUSSIAN TYPE MIXING PROFILE....******************************
    2805. 

  2805. !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
    2806. 

  2806. !  THIS SUBROUTINE INTEGRATES THE AREA UNDER THE CURVE IN THE GAUSSIAN  
    2807. 

  2807. !  DISTRIBUTION...THE NUMERICAL APPROXIMATION TO THE INTEGRAL IS TAKEN FROM
    2808. 

  2808. !  "HANDBOOK OF MATHEMATICAL FUNCTIONS WITH FORMULAS, GRAPHS AND MATHEMATICS TABLES"
    2809. 

  2809. !  ED. BY ABRAMOWITZ AND STEGUN, NATL BUREAU OF STANDARDS APPLIED
    2810. 

  2810. !  MATHEMATICS SERIES.  JUNE, 1964., MAY, 1968.                         
    2811. 

  2811. !                                     JACK KAIN                         
    2812. 

  2812. !                                     7/6/89                            
    2813. 

  2813. !  Solves for KF90 Eq. 2
    2814. 

  2814. !
    2815. 

  2815. !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
    2816. 

  2816. !-----------------------------------------------------------------------
    2817. 

  2817.    IMPLICIT NONE
    2818. 

  2818. !-----------------------------------------------------------------------
    2819. 

  2819.    REAL,         INTENT(IN   )   :: EQ
    2820. 

  2820.    REAL,         INTENT(INOUT)   :: EE,UD
    2821. 

  2821.    REAL ::       SQRT2P,A1,A2,A3,P,SIGMA,FE,X,Y,EY,E45,T1,T2,C1,C2
    2822. 

  2822. 

  2823.       DATA SQRT2P,A1,A2,A3,P,SIGMA,FE/2.506628,0.4361836,-0.1201676,       &
    2824. 

  2824.            0.9372980,0.33267,0.166666667,0.202765151/                        
    2825. 

  2825.       X=(EQ-0.5)/SIGMA                                                  
    2826. 

  2826.       Y=6.*EQ-3.                                                        
    2827. 

  2827.       EY=EXP(Y*Y/(-2))                                                  
    2828. 

  2828.       E45=EXP(-4.5)                                                     
    2829. 

  2829.       T2=1./(1.+P*ABS(Y))                                               
    2830. 

  2830.       T1=0.500498                                                       
    2831. 

  2831.       C1=A1*T1+A2*T1*T1+A3*T1*T1*T1                                     
    2832. 

  2832.       C2=A1*T2+A2*T2*T2+A3*T2*T2*T2                                     
    2833. 

  2833.       IF(Y.GE.0.)THEN                                                   
    2834. 

  2834.         EE=SIGMA*(0.5*(SQRT2P-E45*C1-EY*C2)+SIGMA*(E45-EY))-E45*EQ*EQ/2.
    2835. 

  2835.         UD=SIGMA*(0.5*(EY*C2-E45*C1)+SIGMA*(E45-EY))-E45*(0.5+EQ*EQ/2.-    &
    2836. 

  2836.            EQ)                                                          
    2837. 

  2837.       ELSE                                                              
    2838. 

  2838.         EE=SIGMA*(0.5*(EY*C2-E45*C1)+SIGMA*(E45-EY))-E45*EQ*EQ/2.       
    2839. 

  2839.         UD=SIGMA*(0.5*(SQRT2P-E45*C1-EY*C2)+SIGMA*(E45-EY))-E45*(0.5+EQ*   &
    2840. 

  2840.            EQ/2.-EQ)                                                    
    2841. 

  2841.       ENDIF                                                             
    2842. 

  2842.       EE=EE/FE                                                          
    2843. 

  2843.       UD=UD/FE                                                          
    2844. 

  2844. 

  2845.    END SUBROUTINE PROF5
    2846. 

  2846. 

  2847. ! ------------------------------------------------------------------------
    2848. 

  2848.    SUBROUTINE TPMIX2DD(p,thes,ts,qs,i,j)
    2849. 

  2849. !
    2850. 

  2850. ! Lookup table variables:
    2851. 

  2851. !     INTEGER, PARAMETER :: (KFNT=250,KFNP=220)
    2852. 

  2852. !     REAL, SAVE, DIMENSION(1:KFNT,1:KFNP) :: TTAB,QSTAB
    2853. 

  2853. !     REAL, SAVE, DIMENSION(1:KFNP) :: THE0K
    2854. 

  2854. !     REAL, SAVE, DIMENSION(1:200) :: ALU
    2855. 

  2855. !     REAL, SAVE :: RDPR,RDTHK,PLUTOP
    2856. 

  2856. ! End of Lookup table variables:
    2857. 

  2857. !-----------------------------------------------------------------------
    2858. 

  2858.    IMPLICIT NONE
    2859. 

  2859. !-----------------------------------------------------------------------
    2860. 

  2860.    REAL,         INTENT(IN   )   :: P,THES
    2861. 

  2861.    REAL,         INTENT(INOUT)   :: TS,QS
    2862. 

  2862.    INTEGER,      INTENT(IN   )   :: i,j     ! avail for debugging
    2863. 

  2863.    REAL    ::    TP,QQ,BTH,TTH,PP,T00,T10,T01,T11,Q00,Q10,Q01,Q11
    2864. 

  2864.    INTEGER ::    IPTB,ITHTB
    2865. 

  2865.    CHARACTER*256 :: MESS
    2866. 

  2866. !-----------------------------------------------------------------------
    2867. 

  2867. 

  2868. !
    2869. 

  2869. !******** LOOKUP TABLE VARIABLES (F77 format)... ****************************
    2870. 

  2870. !     parameter(kfnt=250,kfnp=220)
    2871. 

  2871. !
    2872. 

  2872. !     COMMON/KFLUT/ ttab(kfnt,kfnp),qstab(kfnt,kfnp),the0k(kfnp),        &
    2873. 

  2873. !                   alu(200),rdpr,rdthk,plutop 
    2874. 

  2874. !*************************************************************** 
    2875. 

  2875. !
    2876. 

  2876. !***********************************************************************
    2877. 

  2877. !     scaling pressure and tt table index                         
    2878. 

  2878. !***********************************************************************
    2879. 

  2879. !
    2880. 

  2880.       tp=(p-plutop)*rdpr
    2881. 

  2881.       qq=tp-aint(tp)
    2882. 

  2882.       iptb=int(tp)+1
    2883. 

  2883. !
    2884. 

  2884. !***********************************************************************
    2885. 

  2885. !              base and scaling factor for the                           
    2886. 

  2886. !***********************************************************************
    2887. 

  2887. !
    2888. 

  2888. !  scaling the and tt table index                                        
    2889. 

  2889.       bth=(the0k(iptb+1)-the0k(iptb))*qq+the0k(iptb)
    2890. 

  2890.       tth=(thes-bth)*rdthk
    2891. 

  2891.       pp   =tth-aint(tth)
    2892. 

  2892.       ithtb=int(tth)+1
    2893. 

  2893. !
    2894. 

  2894.       t00=ttab(ithtb  ,iptb  )
    2895. 

  2895.       t10=ttab(ithtb+1,iptb  )
    2896. 

  2896.       t01=ttab(ithtb  ,iptb+1)
    2897. 

  2897.       t11=ttab(ithtb+1,iptb+1)
    2898. 

  2898. !
    2899. 

  2899.       q00=qstab(ithtb  ,iptb  )
    2900. 

  2900.       q10=qstab(ithtb+1,iptb  )
    2901. 

  2901.       q01=qstab(ithtb  ,iptb+1)
    2902. 

  2902.       q11=qstab(ithtb+1,iptb+1)
    2903. 

  2903. !
    2904. 

  2904. !***********************************************************************
    2905. 

  2905. !              parcel temperature and saturation mixing ratio                                        
    2906. 

  2906. !***********************************************************************
    2907. 

  2907. !
    2908. 

  2908.       ts=(t00+(t10-t00)*pp+(t01-t00)*qq+(t00-t10-t01+t11)*pp*qq)
    2909. 

  2909. !
    2910. 

  2910.       qs=(q00+(q10-q00)*pp+(q01-q00)*qq+(q00-q10-q01+q11)*pp*qq)
    2911. 

  2911. !
    2912. 

  2912.    END SUBROUTINE TPMIX2DD
    2913. 

  2913. 

  2914. ! -----------------------------------------------------------------------
    2915. 

  2915.   SUBROUTINE ENVIRTHT(P1,T1,Q1,THT1,ALIQ,BLIQ,CLIQ,DLIQ)                       
    2916. 

  2916. !
    2917. 

  2917. !-----------------------------------------------------------------------
    2918. 

  2918.    IMPLICIT NONE
    2919. 

  2919. !-----------------------------------------------------------------------
    2920. 

  2920.    REAL,         INTENT(IN   )   :: P1,T1,Q1,ALIQ,BLIQ,CLIQ,DLIQ
    2921. 

  2921.    REAL,         INTENT(INOUT)   :: THT1
    2922. 

  2922.    REAL    ::    EE,TLOG,ASTRT,AINC,A1,TP,VALUE,AINTRP,TDPT,TSAT,THT,      &
    2923. 

  2923.                  T00,P00,C1,C2,C3,C4,C5
    2924. 

  2924.    INTEGER ::    INDLU
    2925. 

  2925. !-----------------------------------------------------------------------
    2926. 

  2926.       DATA T00,P00,C1,C2,C3,C4,C5/273.16,1.E5,3374.6525,2.5403,3114.834,   &
    2927. 

  2927.            0.278296,1.0723E-3/                                          
    2928. 

  2928. !                                                                       
    2929. 

  2929. !  CALCULATE ENVIRONMENTAL EQUIVALENT POTENTIAL TEMPERATURE...          
    2930. 

  2930. !                                                                       
    2931. 

  2931. ! NOTE: Calculations for mixed/ice phase no longer used...jsk 8/00
    2932. 

  2932. !        For example, KF90 Eq. 10 no longer used
    2933. 

  2933. !
    2934. 

  2934.       EE=Q1*P1/(0.622+Q1)                                             
    2935. 

  2935. !     TLOG=ALOG(EE/ALIQ)                                              
    2936. 

  2936. ! ...calculate LOG term using lookup table...
    2937. 

  2937. !
    2938. 

  2938.       astrt=1.e-3
    2939. 

  2939.       ainc=0.075
    2940. 

  2940.       a1=ee/aliq
    2941. 

  2941.       tp=(a1-astrt)/ainc
    2942. 

  2942.       indlu=int(tp)+1
    2943. 

  2943.       value=(indlu-1)*ainc+astrt
    2944. 

  2944.       aintrp=(a1-value)/ainc
    2945. 

  2945.       tlog=aintrp*alu(indlu+1)+(1-aintrp)*alu(indlu)
    2946. 

  2946. !
    2947. 

  2947.       TDPT=(CLIQ-DLIQ*TLOG)/(BLIQ-TLOG)                               
    2948. 

  2948.       TSAT=TDPT-(.212+1.571E-3*(TDPT-T00)-4.36E-4*(T1-T00))*(T1-TDPT) 
    2949. 

  2949.       THT=T1*(P00/P1)**(0.2854*(1.-0.28*Q1))                          
    2950. 

  2950.       THT1=THT*EXP((C1/TSAT-C2)*Q1*(1.+0.81*Q1))                      
    2951. 

  2951. !
    2952. 

  2952.   END SUBROUTINE ENVIRTHT                                                              
    2953. 

  2953. ! ***********************************************************************
    2954. 

  2954. !====================================================================
    2955. 

  2955.    SUBROUTINE kf_eta_init(RTHCUTEN,RQVCUTEN,RQCCUTEN,RQRCUTEN,      &
    2956. 

  2956.                      RQICUTEN,RQSCUTEN,NCA,W0AVG,P_QI,P_QS,         &
    2957. 

  2957.                      SVP1,SVP2,SVP3,SVPT0,                          &
    2958. 

  2958.                      P_FIRST_SCALAR,restart,allowed_to_read,        &
    2959. 

  2959.                      ids, ide, jds, jde, kds, kde,                  &
    2960. 

  2960.                      ims, ime, jms, jme, kms, kme,                  &
    2961. 

  2961.                      its, ite, jts, jte, kts, kte                   )
    2962. 

  2962. !--------------------------------------------------------------------
    2963. 

  2963.    IMPLICIT NONE
    2964. 

  2964. !--------------------------------------------------------------------
    2965. 

  2965.    LOGICAL , INTENT(IN)           ::  restart,allowed_to_read
    2966. 

  2966.    INTEGER , INTENT(IN)           ::  ids, ide, jds, jde, kds, kde, &
    2967. 

  2967.                                       ims, ime, jms, jme, kms, kme, &
    2968. 

  2968.                                       its, ite, jts, jte, kts, kte
    2969. 

  2969.    INTEGER , INTENT(IN)           ::  P_QI,P_QS,P_FIRST_SCALAR
    2970. 

  2970. 

  2971.    REAL,     DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(OUT) ::       &
    2972. 

  2972.                                                           RTHCUTEN, &
    2973. 

  2973.                                                           RQVCUTEN, &
    2974. 

  2974.                                                           RQCCUTEN, &
    2975. 

  2975.                                                           RQRCUTEN, &
    2976. 

  2976.                                                           RQICUTEN, &
    2977. 

  2977.                                                           RQSCUTEN
    2978. 

  2978. 

  2979.    REAL ,   DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(OUT) :: W0AVG
    2980. 

  2980. 

  2981.    REAL, DIMENSION( ims:ime , jms:jme ), INTENT(INOUT):: NCA
    2982. 

  2982. 

  2983.    INTEGER :: i, j, k, itf, jtf, ktf
    2984. 

  2984.    REAL, INTENT(IN)    :: SVP1,SVP2,SVP3,SVPT0
    2985. 

  2985. 

  2986.    jtf=min0(jte,jde-1)
    2987. 

  2987.    ktf=min0(kte,kde-1)
    2988. 

  2988.    itf=min0(ite,ide-1)
    2989. 

  2989. 

  2990.    IF(.not.restart)THEN
    2991. 

  2991. 

  2992.       DO j=jts,jtf
    2993. 

  2993.       DO k=kts,ktf
    2994. 

  2994.       DO i=its,itf
    2995. 

  2995.          RTHCUTEN(i,k,j)=0.
    2996. 

  2996.          RQVCUTEN(i,k,j)=0.
    2997. 

  2997.          RQCCUTEN(i,k,j)=0.
    2998. 

  2998.          RQRCUTEN(i,k,j)=0.
    2999. 

  2999.       ENDDO
    3000. 

  3000.       ENDDO
    3001. 

  3001.       ENDDO
    3002. 

  3002. 

  3003.       IF (P_QI .ge. P_FIRST_SCALAR) THEN
    3004. 

  3004.          DO j=jts,jtf
    3005. 

  3005.          DO k=kts,ktf
    3006. 

  3006.          DO i=its,itf
    3007. 

  3007.             RQICUTEN(i,k,j)=0.
    3008. 

  3008.          ENDDO
    3009. 

  3009.          ENDDO
    3010. 

  3010.          ENDDO
    3011. 

  3011.       ENDIF
    3012. 

  3012. 

  3013.       IF (P_QS .ge. P_FIRST_SCALAR) THEN
    3014. 

  3014.          DO j=jts,jtf
    3015. 

  3015.          DO k=kts,ktf
    3016. 

  3016.          DO i=its,itf
    3017. 

  3017.             RQSCUTEN(i,k,j)=0.
    3018. 

  3018.          ENDDO
    3019. 

  3019.          ENDDO
    3020. 

  3020.          ENDDO
    3021. 

  3021.       ENDIF
    3022. 

  3022. 

  3023.       DO j=jts,jtf
    3024. 

  3024.       DO i=its,itf
    3025. 

  3025.          NCA(i,j)=-100.
    3026. 

  3026.       ENDDO
    3027. 

  3027.       ENDDO
    3028. 

  3028. 

  3029.       DO j=jts,jtf
    3030. 

  3030.       DO k=kts,ktf
    3031. 

  3031.       DO i=its,itf
    3032. 

  3032.          W0AVG(i,k,j)=0.
    3033. 

  3033.       ENDDO
    3034. 

  3034.       ENDDO
    3035. 

  3035.       ENDDO
    3036. 

  3036. 

  3037.    endif
    3038. 

  3038.  
    3039. 

  3039.    CALL KF_LUTAB(SVP1,SVP2,SVP3,SVPT0)
    3040. 

  3040. 

  3041.    END SUBROUTINE kf_eta_init
    3042. 

  3042. 

  3043. !-------------------------------------------------------
    3044. 

  3044. 

  3045.       subroutine kf_lutab(SVP1,SVP2,SVP3,SVPT0)
    3046. 

  3046. !
    3047. 

  3047. !  This subroutine is a lookup table.
    3048. 

  3048. !  Given a series of series of saturation equivalent potential 
    3049. 

  3049. !  temperatures, the temperature is calculated.
    3050. 

  3050. !
    3051. 

  3051. !--------------------------------------------------------------------
    3052. 

  3052.    IMPLICIT NONE
    3053. 

  3053. !--------------------------------------------------------------------
    3054. 

  3054. ! Lookup table variables
    3055. 

  3055. !     INTEGER, SAVE, PARAMETER :: KFNT=250,KFNP=220
    3056. 

  3056. !     REAL, SAVE, DIMENSION(1:KFNT,1:KFNP) :: TTAB,QSTAB
    3057. 

  3057. !     REAL, SAVE, DIMENSION(1:KFNP) :: THE0K
    3058. 

  3058. !     REAL, SAVE, DIMENSION(1:200) :: ALU
    3059. 

  3059. !     REAL, SAVE :: RDPR,RDTHK,PLUTOP
    3060. 

  3060. ! End of Lookup table variables
    3061. 

  3061. 

  3062.      INTEGER :: KP,IT,ITCNT,I
    3063. 

  3063.      REAL :: DTH,TMIN,TOLER,PBOT,DPR,                               &
    3064. 

  3064.              TEMP,P,ES,QS,PI,THES,TGUES,THGUES,F0,T1,T0,THGS,F1,DT, &
    3065. 

  3065.              ASTRT,AINC,A1,THTGS
    3066. 

  3066. !    REAL    :: ALIQ,BLIQ,CLIQ,DLIQ,SVP1,SVP2,SVP3,SVPT0
    3067. 

  3067.      REAL    :: ALIQ,BLIQ,CLIQ,DLIQ
    3068. 

  3068.      REAL, INTENT(IN)    :: SVP1,SVP2,SVP3,SVPT0
    3069. 

  3069. !
    3070. 

  3070. ! equivalent potential temperature increment
    3071. 

  3071.       data dth/1./
    3072. 

  3072. ! minimum starting temp 
    3073. 

  3073.       data tmin/150./
    3074. 

  3074. ! tolerance for accuracy of temperature 
    3075. 

  3075.       data toler/0.001/
    3076. 

  3076. ! top pressure (pascals)
    3077. 

  3077.       plutop=5000.0
    3078. 

  3078. ! bottom pressure (pascals)
    3079. 

  3079.       pbot=110000.0
    3080. 

  3080. 

  3081.       ALIQ = SVP1*1000.
    3082. 

  3082.       BLIQ = SVP2
    3083. 

  3083.       CLIQ = SVP2*SVPT0
    3084. 

  3084.       DLIQ = SVP3
    3085. 

  3085. 

  3086. !
    3087. 

  3087. ! compute parameters
    3088. 

  3088. !
    3089. 

  3089. ! 1._over_(sat. equiv. theta increment)
    3090. 

  3090.       rdthk=1./dth
    3091. 

  3091. ! pressure increment
    3092. 

  3092. !
    3093. 

  3093.       DPR=(PBOT-PLUTOP)/REAL(KFNP-1)
    3094. 

  3094. !      dpr=(pbot-plutop)/REAL(kfnp-1)
    3095. 

  3095. ! 1._over_(pressure increment)
    3096. 

  3096.       rdpr=1./dpr
    3097. 

  3097. ! compute the spread of thes
    3098. 

  3098. !     thespd=dth*(kfnt-1)
    3099. 

  3099. !
    3100. 

  3100. ! calculate the starting sat. equiv. theta
    3101. 

  3101. !
    3102. 

  3102.       temp=tmin 
    3103. 

  3103.       p=plutop-dpr
    3104. 

  3104.       do kp=1,kfnp
    3105. 

  3105.         p=p+dpr
    3106. 

  3106.         es=aliq*exp((bliq*temp-cliq)/(temp-dliq))
    3107. 

  3107.         qs=0.622*es/(p-es)
    3108. 

  3108.         pi=(1.e5/p)**(0.2854*(1.-0.28*qs))
    3109. 

  3109.         the0k(kp)=temp*pi*exp((3374.6525/temp-2.5403)*qs*        &
    3110. 

  3110.                (1.+0.81*qs))
    3111. 

  3111.       enddo   
    3112. 

  3112. !
    3113. 

  3113. ! compute temperatures for each sat. equiv. potential temp.
    3114. 

  3114. !
    3115. 

  3115.       p=plutop-dpr
    3116. 

  3116.       do kp=1,kfnp
    3117. 

  3117.         thes=the0k(kp)-dth
    3118. 

  3118.         p=p+dpr
    3119. 

  3119.         do it=1,kfnt
    3120. 

  3120. ! define sat. equiv. pot. temp.
    3121. 

  3121.           thes=thes+dth
    3122. 

  3122. ! iterate to find temperature
    3123. 

  3123. ! find initial guess
    3124. 

  3124.           if(it.eq.1) then
    3125. 

  3125.             tgues=tmin
    3126. 

  3126.           else
    3127. 

  3127.             tgues=ttab(it-1,kp)
    3128. 

  3128.           endif
    3129. 

  3129.           es=aliq*exp((bliq*tgues-cliq)/(tgues-dliq))
    3130. 

  3130.           qs=0.622*es/(p-es)
    3131. 

  3131.           pi=(1.e5/p)**(0.2854*(1.-0.28*qs))
    3132. 

  3132.           thgues=tgues*pi*exp((3374.6525/tgues-2.5403)*qs*      &
    3133. 

  3133.                (1.+0.81*qs))
    3134. 

  3134.           f0=thgues-thes
    3135. 

  3135.           t1=tgues-0.5*f0
    3136. 

  3136.           t0=tgues
    3137. 

  3137.           itcnt=0
    3138. 

  3138. ! iteration loop
    3139. 

  3139.           do itcnt=1,11
    3140. 

  3140.             es=aliq*exp((bliq*t1-cliq)/(t1-dliq))
    3141. 

  3141.             qs=0.622*es/(p-es)
    3142. 

  3142.             pi=(1.e5/p)**(0.2854*(1.-0.28*qs))
    3143. 

  3143.             thtgs=t1*pi*exp((3374.6525/t1-2.5403)*qs*(1.+0.81*qs))
    3144. 

  3144.             f1=thtgs-thes
    3145. 

  3145.             if(abs(f1).lt.toler)then
    3146. 

  3146.               exit
    3147. 

  3147.             endif
    3148. 

  3148. !           itcnt=itcnt+1
    3149. 

  3149.             dt=f1*(t1-t0)/(f1-f0)
    3150. 

  3150.             t0=t1
    3151. 

  3151.             f0=f1
    3152. 

  3152.             t1=t1-dt
    3153. 

  3153.           enddo 
    3154. 

  3154.           ttab(it,kp)=t1 
    3155. 

  3155.           qstab(it,kp)=qs
    3156. 

  3156.         enddo
    3157. 

  3157.       enddo   
    3158. 

  3158. !
    3159. 

  3159. ! lookup table for tlog(emix/aliq)
    3160. 

  3160. !
    3161. 

  3161. ! set up intial values for lookup tables
    3162. 

  3162. !
    3163. 

  3163.        astrt=1.e-3
    3164. 

  3164.        ainc=0.075
    3165. 

  3165. !
    3166. 

  3166.        a1=astrt-ainc
    3167. 

  3167.        do i=1,200
    3168. 

  3168.          a1=a1+ainc
    3169. 

  3169.          alu(i)=alog(a1)
    3170. 

  3170.        enddo   
    3171. 

  3171. !
    3172. 

  3172.    END SUBROUTINE KF_LUTAB
    3173. 

  3173. 

  3174. END MODULE module_cu_kfeta
    3175. 

  3175.