PageRenderTime 65ms CodeModel.GetById 24ms RepoModel.GetById 0ms app.codeStats 1ms

/wrfv2_fire/phys/module_cu_osas.F

http://github.com/jbeezley/wrf-fire
FORTRAN Legacy | 2586 lines | 1848 code | 99 blank | 639 comment | 53 complexity | e5e0d51bc82750324794a1aad415c4a7 MD5 | raw file
Possible License(s): AGPL-1.0

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

    

  1. !!
    2. 

  2. MODULE module_cu_osas 
    3. 

  3. 

  4. CONTAINS
    5. 

  5. 

  6. !-----------------------------------------------------------------
    7. 

  7.       SUBROUTINE CU_OSAS(DT,ITIMESTEP,STEPCU,                        &
    8. 

  8.                  RTHCUTEN,RQVCUTEN,RQCCUTEN,RQICUTEN,               &
    9. 

  9.                  RUCUTEN,RVCUTEN,                                   & ! gopal's doing for SAS
    10. 

  10.                  RAINCV,PRATEC,HTOP,HBOT,                           &
    11. 

  11.                  U3D,V3D,W,T3D,QV3D,QC3D,QI3D,PI3D,RHO3D,           &
    12. 

  12.                  DZ8W,PCPS,P8W,XLAND,CU_ACT_FLAG,                   &
    13. 

  13.                  P_QC,                                              & 
    14. 

  14.                  STORE_RAND,MOMMIX, & ! gopal's doing
    15. 

  15.                  P_QI,P_FIRST_SCALAR,                               & 
    16. 

  16.                  CUDT, CURR_SECS, ADAPT_STEP_FLAG,                  &
    17. 

  17.                  CUDTACTTIME,                                       & 
    18. 

  18.                  ids,ide, jds,jde, kds,kde,                         &
    19. 

  19.                  ims,ime, jms,jme, kms,kme,                         &
    20. 

  20.                  its,ite, jts,jte, kts,kte                          )
    21. 

  21. 

  22. !-------------------------------------------------------------------
    23. 

  23.       USE MODULE_GFS_MACHINE , ONLY : kind_phys
    24. 

  24.       USE MODULE_GFS_FUNCPHYS , ONLY : gfuncphys
    25. 

  25.       USE MODULE_GFS_PHYSCONS, grav => con_g, CP => con_CP, HVAP => con_HVAP  &
    26. 

  26.      &,             RV => con_RV, FV => con_fvirt, T0C => con_T0C       &
    27. 

  27.      &,             CVAP => con_CVAP, CLIQ => con_CLIQ                  & 
    28. 

  28.      &,             EPS => con_eps, EPSM1 => con_epsm1                  &
    29. 

  29.      &,             ROVCP => con_rocp, RD => con_rd
    30. 

  30. !-------------------------------------------------------------------
    31. 

  31.       IMPLICIT NONE
    32. 

  32. !-------------------------------------------------------------------
    33. 

  33. !-- U3D         3D u-velocity interpolated to theta points (m/s)
    34. 

  34. !-- V3D         3D v-velocity interpolated to theta points (m/s)
    35. 

  35. !-- TH3D	3D potential temperature (K)
    36. 

  36. !-- T3D         temperature (K)
    37. 

  37. !-- QV3D        3D water vapor mixing ratio (Kg/Kg)
    38. 

  38. !-- QC3D        3D cloud mixing ratio (Kg/Kg)
    39. 

  39. !-- QI3D        3D ice mixing ratio (Kg/Kg)
    40. 

  40. !-- P8w         3D pressure at full levels (Pa)
    41. 

  41. !-- Pcps        3D pressure (Pa)
    42. 

  42. !-- PI3D	3D exner function (dimensionless)
    43. 

  43. !-- rr3D	3D dry air density (kg/m^3)
    44. 

  44. !-- RUBLTEN     U tendency due to
    45. 

  45. !               PBL parameterization (m/s^2)
    46. 

  46. !-- RVBLTEN     V tendency due to
    47. 

  47. !               PBL parameterization (m/s^2)
    48. 

  48. !-- RTHBLTEN    Theta tendency due to
    49. 

  49. !               PBL parameterization (K/s)
    50. 

  50. !-- RQVBLTEN    Qv tendency due to
    51. 

  51. !               PBL parameterization (kg/kg/s)
    52. 

  52. !-- RQCBLTEN    Qc tendency due to
    53. 

  53. !               PBL parameterization (kg/kg/s)
    54. 

  54. !-- RQIBLTEN    Qi tendency due to
    55. 

  55. !               PBL parameterization (kg/kg/s)
    56. 

  56. !
    57. 

  57. !-- MOMMIX      MOMENTUM MIXING COEFFICIENT (can be set in the namelist)
    58. 

  58. !-- RUCUTEN     U tendency due to Cumulus Momentum Mixing (gopal's doing for SAS)
    59. 

  59. !-- RVCUTEN     V tendency due to Cumulus Momentum Mixing (gopal's doing for SAS)
    60. 

  60. !
    61. 

  61. !-- CP          heat capacity at constant pressure for dry air (J/kg/K)
    62. 

  62. !-- GRAV        acceleration due to gravity (m/s^2)
    63. 

  63. !-- ROVCP       R/CP
    64. 

  64. !-- RD          gas constant for dry air (J/kg/K)
    65. 

  65. !-- ROVG 	R/G
    66. 

  66. !-- P_QI	species index for cloud ice
    67. 

  67. !-- dz8w	dz between full levels (m)
    68. 

  68. !-- z		height above sea level (m)
    69. 

  69. !-- PSFC        pressure at the surface (Pa)
    70. 

  70. !-- UST		u* in similarity theory (m/s)
    71. 

  71. !-- PBL		PBL height (m)
    72. 

  72. !-- PSIM        similarity stability function for momentum
    73. 

  73. !-- PSIH        similarity stability function for heat
    74. 

  74. !-- HFX		upward heat flux at the surface (W/m^2)
    75. 

  75. !-- QFX		upward moisture flux at the surface (kg/m^2/s)
    76. 

  76. !-- TSK		surface temperature (K)
    77. 

  77. !-- GZ1OZ0      log(z/z0) where z0 is roughness length
    78. 

  78. !-- WSPD        wind speed at lowest model level (m/s)
    79. 

  79. !-- BR          bulk Richardson number in surface layer
    80. 

  80. !-- DT		time step (s)
    81. 

  81. !-- rvovrd      R_v divided by R_d (dimensionless)
    82. 

  82. !-- EP1         constant for virtual temperature (R_v/R_d - 1) (dimensionless)
    83. 

  83. !-- KARMAN      Von Karman constant
    84. 

  84. !-- ids         start index for i in domain
    85. 

  85. !-- ide         end index for i in domain
    86. 

  86. !-- jds         start index for j in domain
    87. 

  87. !-- jde         end index for j in domain
    88. 

  88. !-- kds         start index for k in domain
    89. 

  89. !-- kde         end index for k in domain
    90. 

  90. !-- ims         start index for i in memory
    91. 

  91. !-- ime         end index for i in memory
    92. 

  92. !-- jms         start index for j in memory
    93. 

  93. !-- jme         end index for j in memory
    94. 

  94. !-- kms         start index for k in memory
    95. 

  95. !-- kme         end index for k in memory
    96. 

  96. !-- its         start index for i in tile
    97. 

  97. !-- ite         end index for i in tile
    98. 

  98. !-- jts         start index for j in tile
    99. 

  99. !-- jte         end index for j in tile
    100. 

  100. !-- kts         start index for k in tile
    101. 

  101. !-- kte         end index for k in tile
    102. 

  102. !-------------------------------------------------------------------
    103. 

  103. 

  104.       INTEGER ::                        ICLDCK
    105. 

  105. 

  106.       INTEGER, INTENT(IN) ::            ids,ide, jds,jde, kds,kde,      &
    107. 

  107.                                         ims,ime, jms,jme, kms,kme,      &
    108. 

  108.                                         its,ite, jts,jte, kts,kte,      &
    109. 

  109.                                         ITIMESTEP,                      &     !NSTD
    110. 

  110.                                         P_FIRST_SCALAR,                 &
    111. 

  111.                                         P_QC,                           &
    112. 

  112.                                         P_QI,                           &
    113. 

  113.                                         STEPCU
    114. 

  114. 

  115.       REAL,    INTENT(IN) ::                                            &
    116. 

  116.                                         DT
    117. 

  117. 

  118. 

  119.       REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT) ::      &
    120. 

  120.                                         RQCCUTEN,                       &
    121. 

  121.                                         RQICUTEN,                       &
    122. 

  122.                                         RQVCUTEN,                       &
    123. 

  123.                                         RTHCUTEN
    124. 

  124.       REAL, DIMENSION(ims:ime, jms:jme, kms:kme), INTENT(INOUT) ::      &
    125. 

  125.                                         RUCUTEN,                        &  ! gopal's doing for SAS
    126. 

  126.                                         RVCUTEN                            ! gopal's doing for SAS 
    127. 

  127.       REAL, OPTIONAL,   INTENT(IN) ::    MOMMIX
    128. 

  128.       REAL, DIMENSION( ims:ime , jms:jme ), OPTIONAL,                   &
    129. 

  129.                          INTENT(IN) :: STORE_RAND
    130. 

  130. 

  131.       REAL,    DIMENSION(ims:ime, jms:jme), INTENT(IN) ::               &
    132. 

  132.                                         XLAND
    133. 

  133. 

  134.       REAL,    DIMENSION(ims:ime, jms:jme), INTENT(INOUT) ::            &
    135. 

  135.                                         RAINCV, PRATEC
    136. 

  136. 

  137.       REAL,    DIMENSION(ims:ime, jms:jme), INTENT(OUT) ::              &
    138. 

  138.                                         HBOT,                           &
    139. 

  139.                                         HTOP
    140. 

  140. 

  141.       LOGICAL, DIMENSION(IMS:IME,JMS:JME), INTENT(INOUT) ::             &
    142. 

  142.                                         CU_ACT_FLAG
    143. 

  143. 

  144. 

  145.       REAL,    DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(IN) ::      &
    146. 

  146.                                         DZ8W,                           &
    147. 

  147.                                         P8w,                            &
    148. 

  148.                                         Pcps,                           &
    149. 

  149.                                         PI3D,                           &
    150. 

  150.                                         QC3D,                           &
    151. 

  151.                                         QI3D,                           &
    152. 

  152.                                         QV3D,                           &
    153. 

  153.                                         RHO3D,                          &
    154. 

  154.                                         T3D,                            &
    155. 

  155.                                         U3D,                            &
    156. 

  156.                                         V3D,                            &
    157. 

  157.                                         W
    158. 

  158. 

  159. ! Adaptive time-step variables
    160. 

  160.       REAL,  INTENT(IN   ) :: CUDT
    161. 

  161.       REAL,  INTENT(IN   ) :: CURR_SECS
    162. 

  162.       LOGICAL,INTENT(IN   ) , OPTIONAL :: ADAPT_STEP_FLAG
    163. 

  163.       REAL,  INTENT (INOUT) :: CUDTACTTIME       
    164. 

  164. 

  165. !--------------------------- LOCAL VARS ------------------------------
    166. 

  166. 

  167.       REAL,    DIMENSION(ims:ime, jms:jme) ::                           &
    168. 

  168.                                         PSFC
    169. 

  169. 

  170. 

  171.       REAL     (kind=kind_phys) ::                                      &
    172. 

  172.                                         DELT,                           &
    173. 

  173.                                         DPSHC,                          &
    174. 

  174.                                         RDELT,                          &
    175. 

  175.                                         RSEED
    176. 

  176. 

  177.       REAL     (kind=kind_phys), DIMENSION(its:ite) ::                  &
    178. 

  178.                                         CLDWRK,                         &
    179. 

  179.                                         PS,                             &
    180. 

  180.                                         RCS,                            &
    181. 

  181.                                         RN,                             &
    182. 

  182.                                         SLIMSK,                         &
    183. 

  183.                                         XKT2
    184. 

  184. 

  185.       REAL     (kind=kind_phys), DIMENSION(its:ite, kts:kte+1) ::       &
    186. 

  186.                                         PRSI                            
    187. 

  187. 

  188.       REAL     (kind=kind_phys), DIMENSION(its:ite, kts:kte) ::         &
    189. 

  189.                                         DEL,                            &
    190. 

  190.                                         DOT,                            &
    191. 

  191.                                         PHIL,                           &
    192. 

  192.                                         PRSL,                           &
    193. 

  193.                                         PRSLK,                          &
    194. 

  194.                                         Q1,                             & 
    195. 

  195.                                         T1,                             & 
    196. 

  196.                                         U1,                             & 
    197. 

  197.                                         V1,                             & 
    198. 

  198.                                         ZI,                             & 
    199. 

  199.                                         ZL 
    200. 

  200. 

  201.       REAL     (kind=kind_phys), DIMENSION(its:ite, kts:kte, 2) ::      &
    202. 

  202.                                         QL 
    203. 

  203. 

  204.       INTEGER, DIMENSION(its:ite) ::                                    &
    205. 

  205.                                         KBOT,                           &
    206. 

  206.                                         KTOP,                           &
    207. 

  207.                                         KUO
    208. 

  208. 

  209.       INTEGER ::                                                        &
    210. 

  210.                                         I,                              &
    211. 

  211.                                         IGPVS,                          &
    212. 

  212.                                         IM,                             &
    213. 

  213.                                         J,                              &
    214. 

  214.                                         JCAP,                           &
    215. 

  215.                                         K,                              &
    216. 

  216.                                         KM,                             &
    217. 

  217.                                         KP,                             &
    218. 

  218.                                         KX,                             &
    219. 

  219.                                         NCLOUD 
    220. 

  220. 

  221.       LOGICAL :: run_param , doing_adapt_dt , decided
    222. 

  222. 

  223.       DATA IGPVS/0/
    224. 

  224. 

  225. !-----------------------------------------------------------------------
    226. 

  226. !
    227. 

  227. !***  CHECK TO SEE IF THIS IS A CONVECTION TIMESTEP
    228. 

  228. !
    229. 

  229. 

  230. !  Initialization for adaptive time step.
    231. 

  231. 

  232.    doing_adapt_dt = .FALSE.
    233. 

  233.    IF ( PRESENT(adapt_step_flag) ) THEN
    234. 

  234.       IF ( adapt_step_flag ) THEN
    235. 

  235.          doing_adapt_dt = .TRUE.
    236. 

  236.          IF ( cudtacttime .EQ. 0. ) THEN
    237. 

  237.             cudtacttime = curr_secs + cudt*60.
    238. 

  238.          END IF
    239. 

  239.       END IF
    240. 

  240.    END IF
    241. 

  241. 

  242. !  Do we run through this scheme or not?
    243. 

  243. 

  244. !    Test 1:  If this is the initial model time, then yes.
    245. 

  245. !                ITIMESTEP=1
    246. 

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

  247. !                CUDT=0 or STEPCU=1
    248. 

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

  249. !                MOD(ITIMESTEP,STEPCU)=0
    250. 

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

  251. !                CURR_SECS >= CUDTACTTIME
    252. 

  252. 

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

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

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

  256. 

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

  258. !  cumulus run.
    259. 

  259. 

  260.    decided = .FALSE.
    261. 

  261.    run_param = .FALSE.
    262. 

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

  263.         ( itimestep .EQ. 1 ) ) THEN
    264. 

  264.       run_param   = .TRUE.
    265. 

  265.       decided     = .TRUE.
    266. 

  266.    END IF
    267. 

  267. 

  268.    IF ( ( .NOT. decided ) .AND. &
    269. 

  269.         ( ( cudt .EQ. 0. ) .OR. ( stepcu .EQ. 1 ) ) ) THEN
    270. 

  270.       run_param   = .TRUE.
    271. 

  271.       decided     = .TRUE.
    272. 

  272.    END IF
    273. 

  273. 

  274.    IF ( ( .NOT. decided ) .AND. &
    275. 

  275.         ( .NOT. doing_adapt_dt ) .AND. &
    276. 

  276.         ( MOD(itimestep,stepcu) .EQ. 0 ) ) THEN
    277. 

  277.       run_param   = .TRUE.
    278. 

  278.       decided     = .TRUE.
    279. 

  279.    END IF
    280. 

  280. 

  281.    IF ( ( .NOT. decided ) .AND. &
    282. 

  282.         ( doing_adapt_dt ) .AND. &
    283. 

  283.         ( curr_secs .GE. cudtacttime ) ) THEN
    284. 

  284.       run_param   = .TRUE.
    285. 

  285.       decided     = .TRUE.
    286. 

  286.       cudtacttime = curr_secs + cudt*60
    287. 

  287.    END IF
    288. 

  288. 

  289. !-----------------------------------------------------------------------
    290. 

  290. 

  291. 

  292.    IF(run_param) THEN
    293. 

  293. 

  294.       DO J=JTS,JTE
    295. 

  295.          DO I=ITS,ITE
    296. 

  296.             CU_ACT_FLAG(I,J)=.TRUE.
    297. 

  297.          ENDDO
    298. 

  298.       ENDDO
    299. 

  299.  
    300. 

  300.       IM=ITE-ITS+1
    301. 

  301.       KX=KTE-KTS+1
    302. 

  302.       JCAP=126
    303. 

  303.       DPSHC=30_kind_phys
    304. 

  304.       DELT=DT*STEPCU
    305. 

  305.       RDELT=1./DELT
    306. 

  306.       NCLOUD=1
    307. 

  307. 

  308. 

  309.    DO J=jms,jme
    310. 

  310.      DO I=ims,ime
    311. 

  311.        PSFC(i,j)=P8w(i,kms,j)
    312. 

  312.      ENDDO
    313. 

  313.    ENDDO
    314. 

  314. 

  315.    if(igpvs.eq.0) CALL GFUNCPHYS
    316. 

  316.    igpvs=1
    317. 

  317. 

  318. !-------------  J LOOP (OUTER) --------------------------------------------------
    319. 

  319. 

  320.    DO J=jts,jte
    321. 

  321. 

  322. ! --------------- compute zi and zl -----------------------------------------
    323. 

  323.       DO i=its,ite
    324. 

  324.         ZI(I,KTS)=0.0
    325. 

  325.       ENDDO
    326. 

  326. 

  327.       DO k=kts+1,kte
    328. 

  328.         KM=K-1
    329. 

  329.         DO i=its,ite
    330. 

  330.           ZI(I,K)=ZI(I,KM)+dz8w(i,km,j)
    331. 

  331.         ENDDO
    332. 

  332.       ENDDO
    333. 

  333. 

  334.       DO k=kts+1,kte
    335. 

  335.         KM=K-1
    336. 

  336.         DO i=its,ite
    337. 

  337.           ZL(I,KM)=(ZI(I,K)+ZI(I,KM))*0.5
    338. 

  338.         ENDDO
    339. 

  339.       ENDDO
    340. 

  340. 

  341.       DO i=its,ite
    342. 

  342.         ZL(I,KTE)=2.*ZI(I,KTE)-ZL(I,KTE-1)
    343. 

  343.       ENDDO
    344. 

  344. 

  345. ! --------------- end compute zi and zl -------------------------------------
    346. 

  346. 

  347. !    Based on some important findings from Morris Bender, XKT2 was defined in
    348. 

  348. !    terms of random number instead of random number based cloud tops
    349. 

  349. !    Also, these random numbers are stored and are changed only once in
    350. 

  350. !    approximately 5 minutes interval now. This is gopal's doing for HWRF.
    351. 

    352. 

  352. !     call random_number(XKT2)
    353. 

  353. 

  354. #if (EM_CORE == 1)
    355. 

  355. !    XKT2 was defined in terms of random number instead of random number based cloud tops
    356. 

  356. !    ZCX   
    357. 

  357.      call init_random_seed()
    358. 

  358.      call random_number(XKT2)
    359. 

  359. #ifdef REGTEST
    360. 

  360. ! for regtest only
    361. 

  361.      xkt2 = 0.1
    362. 

  362. #endif
    363. 

  363. #endif
    364. 

  364. !   
    365. 

  365. #if (NMM_CORE == 1)
    366. 

  366.       DO i=its,ite
    367. 

  367.          XKT2(i) = STORE_RAND(i,j)
    368. 

  368.       ENDDO
    369. 

  369. #endif
    370. 

  370. 

  371.       DO i=its,ite
    372. 

  372.         PS(i)=PSFC(i,j)*.001
    373. 

  373.         RCS(i)=1.
    374. 

  374.         SLIMSK(i)=ABS(XLAND(i,j)-2.)
    375. 

  375.       ENDDO
    376. 

  376. 

  377.       DO i=its,ite
    378. 

  378.         PRSI(i,kts)=PS(i)
    379. 

  379.       ENDDO
    380. 

  380. 

  381.       DO k=kts,kte
    382. 

  382.         kp=k+1
    383. 

  383.         DO i=its,ite
    384. 

  384.           PRSL(I,K)=Pcps(i,k,j)*.001
    385. 

  385.           PHIL(I,K)=ZL(I,K)*GRAV
    386. 

  386.           DOT(i,k)=-5.0E-4*GRAV*rho3d(i,k,j)*(w(i,k,j)+w(i,kp,j))
    387. 

  387.         ENDDO
    388. 

  388.       ENDDO
    389. 

  389. 

  390.       DO k=kts,kte
    391. 

  391.         DO i=its,ite
    392. 

  392.           DEL(i,k)=PRSL(i,k)*GRAV/RD*dz8w(i,k,j)/T3D(i,k,j)
    393. 

  393.           U1(i,k)=U3D(i,k,j)
    394. 

  394.           V1(i,k)=V3D(i,k,j)
    395. 

  395.           Q1(i,k)=QV3D(i,k,j)/(1.+QV3D(i,k,j))
    396. 

  396.           T1(i,k)=T3D(i,k,j)
    397. 

  397.           QL(i,k,1)=QI3D(i,k,j)/(1.+QI3D(i,k,j))
    398. 

  398.           QL(i,k,2)=QC3D(i,k,j)/(1.+QC3D(i,k,j))
    399. 

  399.           PRSLK(I,K)=(PRSL(i,k)*.01)**ROVCP
    400. 

  400.         ENDDO
    401. 

  401.       ENDDO
    402. 

  402. 

  403.       DO k=kts+1,kte+1
    404. 

  404.         km=k-1
    405. 

  405.         DO i=its,ite
    406. 

  406.           PRSI(i,k)=PRSI(i,km)-del(i,km) 
    407. 

  407.         ENDDO
    408. 

  408.       ENDDO
    409. 

  409. 

  410. 

  411.       CALL OSASCNV(IM,IM,KX,JCAP,DELT,DEL,PRSL,PS,PHIL,                  &
    412. 

  412.                   QL,Q1,T1,U1,V1,RCS,CLDWRK,RN,KBOT,                    &
    413. 

  413.                   KTOP,KUO,SLIMSK,DOT,XKT2,NCLOUD) 
    414. 

  414. 

  415. !!!   make more like GFDL ... eliminate shallow convection.....
    416. 

  416. !!!   CALL SHALCV(IM,IM,KX,DELT,DEL,PRSI,PRSL,PRSLK,KUO,Q1,T1,DPSHC)
    417. 

  417. #if (EM_CORE == 1)
    418. 

  418.       CALL SHALCV(IM,IM,KX,DELT,DEL,PRSI,PRSL,PRSLK,KUO,Q1,T1,DPSHC)
    419. 

  419. #endif
    420. 

  420. 

  421.       DO I=ITS,ITE
    422. 

  422.         RAINCV(I,J)=RN(I)*1000./STEPCU
    423. 

  423.         PRATEC(I,J)=RN(I)*1000./(STEPCU * DT)
    424. 

  424.         HBOT(I,J)=KBOT(I)
    425. 

  425.         HTOP(I,J)=KTOP(I)
    426. 

  426.       ENDDO
    427. 

  427. 

  428.       DO K=KTS,KTE
    429. 

  429.         DO I=ITS,ITE
    430. 

  430.           RTHCUTEN(I,K,J)=(T1(I,K)-T3D(I,K,J))/PI3D(I,K,J)*RDELT
    431. 

  431.           RQVCUTEN(I,K,J)=(Q1(I,K)/(1.-q1(i,k))-QV3D(I,K,J))*RDELT
    432. 

  432.         ENDDO
    433. 

  433.       ENDDO
    434. 

  434. 

  435. !===============================================================================
    436. 

  436. !     ADD MOMENTUM MIXING TERM AS TENDENCIES. This is gopal's doing for SAS
    437. 

  437. !     MOMMIX is the reduction factor set to 0.7 by default. Because NMM has 
    438. 

  438. !     divergence damping term, a reducion factor for cumulum mixing may be
    439. 

  439. !     required otherwise storms were too weak.
    440. 

  440. !===============================================================================
    441. 

  441. !
    442. 

  442. #if (NMM_CORE == 1)
    443. 

  443.       DO K=KTS,KTE
    444. 

  444.         DO I=ITS,ITE
    445. 

  445.          RUCUTEN(I,J,K)=MOMMIX*(U1(I,K)-U3D(I,K,J))*RDELT
    446. 

  446.          RVCUTEN(I,J,K)=MOMMIX*(V1(I,K)-V3D(I,K,J))*RDELT
    447. 

  447.         ENDDO
    448. 

  448.       ENDDO
    449. 

  449. #endif
    450. 

  450. 

  451. 

  452.       IF(P_QC .ge. P_FIRST_SCALAR)THEN
    453. 

  453.         DO K=KTS,KTE
    454. 

  454.           DO I=ITS,ITE
    455. 

  455.             RQCCUTEN(I,K,J)=(QL(I,K,2)/(1.-ql(i,k,2))-QC3D(I,K,J))*RDELT
    456. 

  456.           ENDDO
    457. 

  457.         ENDDO
    458. 

  458.       ENDIF
    459. 

  459. 

  460.       IF(P_QI .ge. P_FIRST_SCALAR)THEN
    461. 

  461.         DO K=KTS,KTE
    462. 

  462.           DO I=ITS,ITE
    463. 

  463.             RQICUTEN(I,K,J)=(QL(I,K,1)/(1.-ql(i,k,1))-QI3D(I,K,J))*RDELT
    464. 

  464.           ENDDO
    465. 

  465.         ENDDO
    466. 

  466.       ENDIF
    467. 

  467. 

  468.    ENDDO    ! Outer most J loop
    469. 

  469. 

  470.    ENDIF
    471. 

  471. 

  472.    END SUBROUTINE CU_OSAS
    473. 

  473. 

  474. !====================================================================
    475. 

  475.    SUBROUTINE osasinit(RTHCUTEN,RQVCUTEN,RQCCUTEN,RQICUTEN,          &
    476. 

  476.                       RUCUTEN,RVCUTEN,                              &   ! gopal's doing for SAS
    477. 

  477.                       RESTART,P_QC,P_QI,P_FIRST_SCALAR,             &
    478. 

  478.                       allowed_to_read,                              &
    479. 

  479.                       ids, ide, jds, jde, kds, kde,                 &
    480. 

  480.                       ims, ime, jms, jme, kms, kme,                 &
    481. 

  481.                       its, ite, jts, jte, kts, kte                  )
    482. 

  482. !--------------------------------------------------------------------
    483. 

  483.    IMPLICIT NONE
    484. 

  484. !--------------------------------------------------------------------
    485. 

  485.    LOGICAL , INTENT(IN)           ::  allowed_to_read,restart
    486. 

  486.    INTEGER , INTENT(IN)           ::  ids, ide, jds, jde, kds, kde, &
    487. 

  487.                                       ims, ime, jms, jme, kms, kme, &
    488. 

  488.                                       its, ite, jts, jte, kts, kte
    489. 

  489.    INTEGER , INTENT(IN)           ::  P_FIRST_SCALAR, P_QI, P_QC
    490. 

  490. 

  491.    REAL,     DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(OUT) ::  &
    492. 

  492.                                                               RTHCUTEN, &
    493. 

  493.                                                               RQVCUTEN, &
    494. 

  494.                                                               RQCCUTEN, &
    495. 

  495.                                                               RQICUTEN
    496. 

  496.    REAL,     DIMENSION( ims:ime , jms:jme , kms:kme ) , INTENT(OUT) ::  &
    497. 

  497.                                                               RUCUTEN,  & ! gopal's doing for SAS
    498. 

  498.                                                               RVCUTEN   
    499. 

  499. 

  500.    INTEGER :: i, j, k, itf, jtf, ktf
    501. 

  501. 

  502.    jtf=min0(jte,jde-1)
    503. 

  503.    ktf=min0(kte,kde-1)
    504. 

  504.    itf=min0(ite,ide-1)
    505. 

  505. 

  506. #ifdef HWRF
    507. 

  507. !zhang's doing
    508. 

  508.    IF(.not.restart .or. .not.allowed_to_read)THEN
    509. 

  509. !end of zhang's doing
    510. 

  510. #else
    511. 

  511.    IF(.not.restart)THEN
    512. 

  512. #endif
    513. 

  513.      DO j=jts,jtf
    514. 

  514.      DO k=kts,ktf
    515. 

  515.      DO i=its,itf
    516. 

  516.        RTHCUTEN(i,k,j)=0.
    517. 

  517.        RQVCUTEN(i,k,j)=0.
    518. 

  518.        RUCUTEN(i,j,k)=0.   ! gopal's doing for SAS
    519. 

  519.        RVCUTEN(i,j,k)=0.    ! gopal's doing for SAS
    520. 

  520.      ENDDO
    521. 

  521.      ENDDO
    522. 

  522.      ENDDO
    523. 

  523. 

  524.      IF (P_QC .ge. P_FIRST_SCALAR) THEN
    525. 

  525.         DO j=jts,jtf
    526. 

  526.         DO k=kts,ktf
    527. 

  527.         DO i=its,itf
    528. 

  528.            RQCCUTEN(i,k,j)=0.
    529. 

  529.         ENDDO
    530. 

  530.         ENDDO
    531. 

  531.         ENDDO
    532. 

  532.      ENDIF
    533. 

  533. 

  534.      IF (P_QI .ge. P_FIRST_SCALAR) THEN
    535. 

  535.         DO j=jts,jtf
    536. 

  536.         DO k=kts,ktf
    537. 

  537.         DO i=its,itf
    538. 

  538.            RQICUTEN(i,k,j)=0.
    539. 

  539.         ENDDO
    540. 

  540.         ENDDO
    541. 

  541.         ENDDO
    542. 

  542.      ENDIF
    543. 

  543.    ENDIF
    544. 

  544. 

  545.       END SUBROUTINE osasinit
    546. 

  546. 

  547. ! ------------------------------------------------------------------------
    548. 

  548. 

  549.       SUBROUTINE OSASCNV(IM,IX,KM,JCAP,DELT,DEL,PRSL,PS,PHIL,QL,         &
    550. 

  550.      &       Q1,T1,U1,V1,RCS,CLDWRK,RN,KBOT,KTOP,KUO,SLIMSK,            &
    551. 

  551.      &       DOT,XKT2,ncloud)
    552. 

  552. !  for cloud water version
    553. 

  553. !     parameter(ncloud=0)
    554. 

  554. !     SUBROUTINE OSASCNV(KM,JCAP,DELT,DEL,SL,SLK,PS,QL,
    555. 

  555. !    &       Q1,T1,U1,V1,RCS,CLDWRK,RN,KBOT,KTOP,KUO,SLIMSK,
    556. 

  556. !    &       DOT,xkt2,ncloud)
    557. 

  557. !
    558. 

  558.       USE MODULE_GFS_MACHINE , ONLY : kind_phys
    559. 

  559.       USE MODULE_GFS_FUNCPHYS ,ONLY : fpvs
    560. 

  560.       USE MODULE_GFS_PHYSCONS, grav => con_g, CP => con_CP, HVAP => con_HVAP &
    561. 

  561.      &,             RV => con_RV, FV => con_fvirt, T0C => con_T0C       &
    562. 

  562.      &,             CVAP => con_CVAP, CLIQ => con_CLIQ                  &
    563. 

  563.      &,             EPS => con_eps, EPSM1 => con_epsm1
    564. 

  564. 

  565.       implicit none
    566. 

  566. !
    567. 

  567. !     include 'constant.h'
    568. 

  568. !
    569. 

  569.       integer            IM, IX,  KM, JCAP, ncloud,                     &
    570. 

  570.      &                   KBOT(IM), KTOP(IM), KUO(IM), J
    571. 

  571.       real(kind=kind_phys) DELT
    572. 

  572.       real(kind=kind_phys) PS(IM),      DEL(IX,KM),  PRSL(IX,KM),       &
    573. 

  573. !     real(kind=kind_phys)              DEL(IX,KM),  PRSL(IX,KM),
    574. 

  574.      &                     QL(IX,KM,2), Q1(IX,KM),   T1(IX,KM),         &
    575. 

  575.      &                     U1(IX,KM),   V1(IX,KM),   RCS(IM),           &
    576. 

  576.      &                     CLDWRK(IM),  RN(IM),      SLIMSK(IM),        &
    577. 

  577.      &                     DOT(IX,KM),  XKT2(IM),    PHIL(IX,KM)
    578. 

  578. !
    579. 

  579.       integer              I, INDX, jmn, k, knumb, latd, lond, km1
    580. 

  580. !
    581. 

  581.       real(kind=kind_phys) adw,     alpha,   alphal,  alphas,           &
    582. 

  582.      &                     aup,     beta,    betal,   betas,            &
    583. 

  583.      &                     c0,      cpoel,   dellat,  delta,            &
    584. 

  584.      &                     desdt,   deta,    detad,   dg,               &
    585. 

  585.      &                     dh,      dhh,     dlnsig,  dp,               &
    586. 

  586.      &                     dq,      dqsdp,   dqsdt,   dt,               &
    587. 

  587.      &                     dt2,     dtmax,   dtmin,   dv1,              &
    588. 

  588.      &                     dv1q,    dv2,     dv2q,    dv1u,             &
    589. 

  589.      &                     dv1v,    dv2u,    dv2v,    dv3u,             &
    590. 

  590.      &                     dv3v,    dv3,     dv3q,    dvq1,             &
    591. 

  591.      &                     dz,      dz1,     e1,      edtmax,           &
    592. 

  592.      &                     edtmaxl, edtmaxs, el2orc,  elocp,            &
    593. 

  593.      &                     es,      etah,                               &
    594. 

  594.      &                     evef,    evfact,  evfactl, fact1,            &
    595. 

  595.      &                     fact2,   factor,  fjcap,   fkm,              &
    596. 

  596.      &                     fuv,     g,       gamma,   onemf,            &
    597. 

  597.      &                     onemfu,  pdetrn,  pdpdwn,  pprime,           &
    598. 

  598.      &                     qc,      qlk,     qrch,    qs,               &
    599. 

  599.      &                     rain,    rfact,   shear,   tem1,             &
    600. 

  600.      &                     tem2,    terr,    val,     val1,             &
    601. 

  601.      &                     val2,    w1,      w1l,     w1s,              &
    602. 

  602.      &                     w2,      w2l,     w2s,     w3,               &
    603. 

  603.      &                     w3l,     w3s,     w4,      w4l,              & 
    604. 

  604.      &                     w4s,     xdby,    xpw,     xpwd,             & 
    605. 

  605.      &                     xqc,     xqrch,   xlambu,  mbdt,             &
    606. 

  606.      &                     tem
    607. 

  607. !
    608. 

  608. !
    609. 

  609.       integer              JMIN(IM), KB(IM), KBCON(IM), KBDTR(IM),      & 
    610. 

  610.      &                     KT2(IM),  KTCON(IM), LMIN(IM),               &
    611. 

  611.      &                     kbm(IM),  kbmax(IM), kmax(IM)
    612. 

  612. !
    613. 

  613.       real(kind=kind_phys) AA1(IM),     ACRT(IM),   ACRTFCT(IM),        & 
    614. 

  614.      &                     DELHBAR(IM), DELQ(IM),   DELQ2(IM),          &
    615. 

  615.      &                     DELQBAR(IM), DELQEV(IM), DELTBAR(IM),        &
    616. 

  616.      &                     DELTV(IM),   DTCONV(IM), EDT(IM),            &
    617. 

  617.      &                     EDTO(IM),    EDTX(IM),   FLD(IM),            &
    618. 

  618.      &                     HCDO(IM),    HKBO(IM),   HMAX(IM),           &
    619. 

  619.      &                     HMIN(IM),    HSBAR(IM),  UCDO(IM),           &
    620. 

  620.      &                     UKBO(IM),    VCDO(IM),   VKBO(IM),           &
    621. 

  621.      &                     PBCDIF(IM),  PDOT(IM),   PO(IM,KM),          &
    622. 

  622.      &                                  PWAVO(IM),  PWEVO(IM),          &
    623. 

  623. !    &                     PSFC(IM),    PWAVO(IM),  PWEVO(IM),          &
    624. 

  624.      &                     QCDO(IM),    QCOND(IM),  QEVAP(IM),          &
    625. 

  625.      &                     QKBO(IM),    RNTOT(IM),  VSHEAR(IM),         &
    626. 

  626.      &                     XAA0(IM),    XHCD(IM),   XHKB(IM),           & 
    627. 

  627.      &                     XK(IM),      XLAMB(IM),  XLAMD(IM),          &
    628. 

  628.      &                     XMB(IM),     XMBMAX(IM), XPWAV(IM),          &
    629. 

  629.      &                     XPWEV(IM),   XQCD(IM),   XQKB(IM)
    630. 

  630. !
    631. 

  631. !  PHYSICAL PARAMETERS
    632. 

  632.       PARAMETER(G=grav)
    633. 

  633.       PARAMETER(CPOEL=CP/HVAP,ELOCP=HVAP/CP,                            &
    634. 

  634.      &          EL2ORC=HVAP*HVAP/(RV*CP))
    635. 

  635.       PARAMETER(TERR=0.,C0=.002,DELTA=fv)
    636. 

  636.       PARAMETER(FACT1=(CVAP-CLIQ)/RV,FACT2=HVAP/RV-FACT1*T0C)
    637. 

  637. !  LOCAL VARIABLES AND ARRAYS
    638. 

  638.       real(kind=kind_phys) PFLD(IM,KM),    TO(IM,KM),     QO(IM,KM),    &
    639. 

  639.      &                     UO(IM,KM),      VO(IM,KM),     QESO(IM,KM)
    640. 

  640. !  cloud water
    641. 

  641.       real(kind=kind_phys) QLKO_KTCON(IM), DELLAL(IM),    TVO(IM,KM),   &
    642. 

  642.      &                     DBYO(IM,KM),    ZO(IM,KM),     SUMZ(IM,KM),  &
    643. 

  643.      &                     SUMH(IM,KM),    HEO(IM,KM),    HESO(IM,KM),  &
    644. 

  644.      &                     QRCD(IM,KM),    DELLAH(IM,KM), DELLAQ(IM,KM),&
    645. 

  645.      &                     DELLAU(IM,KM),  DELLAV(IM,KM), HCKO(IM,KM),  &
    646. 

  646.      &                     UCKO(IM,KM),    VCKO(IM,KM),   QCKO(IM,KM),  &
    647. 

  647.      &                     ETA(IM,KM),     ETAU(IM,KM),   ETAD(IM,KM),  &
    648. 

  648.      &                     QRCDO(IM,KM),   PWO(IM,KM),    PWDO(IM,KM),  &
    649. 

  649.      &                     RHBAR(IM),      TX1(IM)
    650. 

  650. !
    651. 

  651.       LOGICAL TOTFLG, CNVFLG(IM), DWNFLG(IM), DWNFLG2(IM), FLG(IM)
    652. 

  652. !
    653. 

  653.       real(kind=kind_phys) PCRIT(15), ACRITT(15), ACRIT(15)
    654. 

  654. !     SAVE PCRIT, ACRITT
    655. 

  655.       DATA PCRIT/850.,800.,750.,700.,650.,600.,550.,500.,450.,400.,     &
    656. 

  656.      &           350.,300.,250.,200.,150./
    657. 

  657.       DATA ACRITT/.0633,.0445,.0553,.0664,.075,.1082,.1521,.2216,       &
    658. 

  658.      &           .3151,.3677,.41,.5255,.7663,1.1686,1.6851/
    659. 

  659. !  GDAS DERIVED ACRIT
    660. 

  660. !     DATA ACRITT/.203,.515,.521,.566,.625,.665,.659,.688,              & 
    661. 

  661. !    &            .743,.813,.886,.947,1.138,1.377,1.896/
    662. 

  662. !
    663. 

  663.       real(kind=kind_phys) TF, TCR, TCRF, RZERO, RONE
    664. 

  664.       parameter (TF=233.16, TCR=263.16, TCRF=1.0/(TCR-TF))
    665. 

  665.       parameter (RZERO=0.0,RONE=1.0)
    666. 

  666. !-----------------------------------------------------------------------
    667. 

  667. !
    668. 

  668.       km1 = km - 1
    669. 

  669. !  INITIALIZE ARRAYS
    670. 

  670. !
    671. 

  671.       DO I=1,IM
    672. 

  672.         RN(I)=0.
    673. 

  673.         KBOT(I)=KM+1
    674. 

  674.         KTOP(I)=0
    675. 

  675.         KUO(I)=0
    676. 

  676.         CNVFLG(I) = .TRUE.
    677. 

  677.         DTCONV(I) = 3600.
    678. 

  678.         CLDWRK(I) = 0.
    679. 

  679.         PDOT(I) = 0.
    680. 

  680.         KT2(I) = 0
    681. 

  681.         QLKO_KTCON(I) = 0.
    682. 

  682.         DELLAL(I) = 0.
    683. 

  683.       ENDDO
    684. 

  684. !!
    685. 

  685.       DO K = 1, 15
    686. 

  686.         ACRIT(K) = ACRITT(K) * (975. - PCRIT(K))
    687. 

  687.       ENDDO
    688. 

  688.       DT2 = DELT
    689. 

  689. !cmr  dtmin = max(dt2,1200.)
    690. 

  690.       val   =         1200.
    691. 

  691.       dtmin = max(dt2, val )
    692. 

  692. !cmr  dtmax = max(dt2,3600.)
    693. 

  693.       val   =         3600.
    694. 

  694.       dtmax = max(dt2, val )
    695. 

  695. !  MODEL TUNABLE PARAMETERS ARE ALL HERE
    696. 

  696.       MBDT    = 10.
    697. 

  697.       EDTMAXl = .3
    698. 

  698.       EDTMAXs = .3
    699. 

  699.       ALPHAl  = .5
    700. 

  700.       ALPHAs  = .5
    701. 

  701.       BETAl   = .15
    702. 

  702.       betas   = .15
    703. 

  703.       BETAl   = .05
    704. 

  704.       betas   = .05
    705. 

  705. !     change for hurricane model
    706. 

  706.         BETAl = .5
    707. 

  707.         betas = .5
    708. 

  708. !     EVEF    = 0.07
    709. 

  709.       evfact  = 0.3
    710. 

  710.       evfactl = 0.3
    711. 

  711. !     change for hurricane model
    712. 

  712.          evfact = 0.6
    713. 

  713.          evfactl = .6
    714. 

  714. #if ( EM_CORE == 1 )
    715. 

  715. !  HAWAII TEST - ZCX
    716. 

  716.       ALPHAl  = .5
    717. 

  717.       ALPHAs  = .75
    718. 

  718.       BETAl   = .05
    719. 

  719.       betas   = .05
    720. 

  720.       evfact  = 0.5
    721. 

  721.       evfactl = 0.5
    722. 

  722. #endif
    723. 

  723.       PDPDWN  = 0.
    724. 

  724.       PDETRN  = 200.
    725. 

  725.       xlambu  = 1.e-4
    726. 

  726.       fjcap   = (float(jcap) / 126.) ** 2
    727. 

  727. !cmr  fjcap   = max(fjcap,1.)
    728. 

  728.       val     =           1.
    729. 

  729.       fjcap   = max(fjcap,val)
    730. 

  730.       fkm     = (float(km) / 28.) ** 2
    731. 

  731. !cmr  fkm     = max(fkm,1.)
    732. 

  732.       fkm     = max(fkm,val)
    733. 

  733.       W1l     = -8.E-3 
    734. 

  734.       W2l     = -4.E-2
    735. 

  735.       W3l     = -5.E-3 
    736. 

  736.       W4l     = -5.E-4
    737. 

  737.       W1s     = -2.E-4
    738. 

  738.       W2s     = -2.E-3
    739. 

  739.       W3s     = -1.E-3
    740. 

  740.       W4s     = -2.E-5
    741. 

  741. !CCCC IF(IM.EQ.384) THEN
    742. 

  742.         LATD  = 92
    743. 

  743.         lond  = 189
    744. 

  744. !CCCC ELSEIF(IM.EQ.768) THEN
    745. 

  745. !CCCC   LATD = 80
    746. 

  746. !CCCC ELSE
    747. 

  747. !CCCC   LATD = 0
    748. 

  748. !CCCC ENDIF
    749. 

  749. !
    750. 

  750. !  DEFINE TOP LAYER FOR SEARCH OF THE DOWNDRAFT ORIGINATING LAYER
    751. 

  751. !  AND THE MAXIMUM THETAE FOR UPDRAFT
    752. 

  752. !
    753. 

  753.       DO I=1,IM
    754. 

  754.         KBMAX(I) = KM
    755. 

  755.         KBM(I)   = KM
    756. 

  756.         KMAX(I)  = KM
    757. 

  757.         TX1(I)   = 1.0 / PS(I)
    758. 

  758.       ENDDO
    759. 

  759. !     
    760. 

  760.       DO K = 1, KM
    761. 

  761.         DO I=1,IM
    762. 

  762.           IF (prSL(I,K)*tx1(I) .GT. 0.45) KBMAX(I) = K + 1
    763. 

  763.           IF (prSL(I,K)*tx1(I) .GT. 0.70) KBM(I)   = K + 1
    764. 

  764.           IF (prSL(I,K)*tx1(I) .GT. 0.04) KMAX(I)  = MIN(KM,K + 1)
    765. 

  765.         ENDDO
    766. 

  766.       ENDDO
    767. 

  767.       DO I=1,IM
    768. 

  768.         KBMAX(I) = MIN(KBMAX(I),KMAX(I))
    769. 

  769.         KBM(I)   = MIN(KBM(I),KMAX(I))
    770. 

  770.       ENDDO
    771. 

  771. !
    772. 

  772. !   CONVERT SURFACE PRESSURE TO MB FROM CB
    773. 

  773. !
    774. 

  774. !!
    775. 

  775.       DO K = 1, KM
    776. 

  776.         DO I=1,IM
    777. 

  777.           if (K .le. kmax(i)) then
    778. 

  778.             PFLD(I,k) = PRSL(I,K) * 10.0
    779. 

  779.             PWO(I,k)  = 0.
    780. 

  780.             PWDO(I,k) = 0.
    781. 

  781.             TO(I,k)   = T1(I,k)
    782. 

  782.             QO(I,k)   = Q1(I,k)
    783. 

  783.             UO(I,k)   = U1(I,k)
    784. 

  784.             VO(I,k)   = V1(I,k)
    785. 

  785.             DBYO(I,k) = 0.
    786. 

  786.             SUMZ(I,k) = 0.
    787. 

  787.             SUMH(I,k) = 0.
    788. 

  788.           endif
    789. 

  789.         ENDDO
    790. 

  790.       ENDDO
    791. 

  791. 

  792. !
    793. 

  793. !  COLUMN VARIABLES
    794. 

  794. !  P IS PRESSURE OF THE LAYER (MB)
    795. 

  795. !  T IS TEMPERATURE AT T-DT (K)..TN
    796. 

  796. !  Q IS MIXING RATIO AT T-DT (KG/KG)..QN
    797. 

  797. !  TO IS TEMPERATURE AT T+DT (K)... THIS IS AFTER ADVECTION AND TURBULAN
    798. 

  798. !  QO IS MIXING RATIO AT T+DT (KG/KG)..Q1
    799. 

  799. !
    800. 

  800.       DO K = 1, KM
    801. 

  801.         DO I=1,IM
    802. 

  802.           if (k .le. kmax(i)) then
    803. 

  803.          !jfe        QESO(I,k) = 10. * FPVS(T1(I,k))
    804. 

  804.          !
    805. 

  805.             QESO(I,k) = 0.01 * fpvs(T1(I,K))      ! fpvs is in Pa
    806. 

  806.          !
    807. 

  807.             QESO(I,k) = EPS * QESO(I,k) / (PFLD(I,k) + EPSM1*QESO(I,k))
    808. 

  808.          !cmr        QESO(I,k) = MAX(QESO(I,k),1.E-8)
    809. 

  809.             val1      =             1.E-8
    810. 

  810.             QESO(I,k) = MAX(QESO(I,k), val1)
    811. 

  811.          !cmr        QO(I,k)   = max(QO(I,k),1.e-10)
    812. 

  812.             val2      =           1.e-10
    813. 

  813.             QO(I,k)   = max(QO(I,k), val2 )
    814. 

  814.          !           QO(I,k)   = MIN(QO(I,k),QESO(I,k))
    815. 

  815.             TVO(I,k)  = TO(I,k) + DELTA * TO(I,k) * QO(I,k)
    816. 

  816.           endif
    817. 

  817.         ENDDO
    818. 

  818.       ENDDO
    819. 

  819. 

  820. !
    821. 

  821. !  HYDROSTATIC HEIGHT ASSUME ZERO TERR
    822. 

  822. !
    823. 

  823.       DO K = 1, KM
    824. 

  824.         DO I=1,IM
    825. 

  825.           ZO(I,k) = PHIL(I,k) / G
    826. 

  826.         ENDDO
    827. 

  827.       ENDDO
    828. 

  828. !  COMPUTE MOIST STATIC ENERGY
    829. 

  829.       DO K = 1, KM
    830. 

  830.         DO I=1,IM
    831. 

  831.           if (K .le. kmax(i)) then
    832. 

  832. !           tem       = G * ZO(I,k) + CP * TO(I,k)
    833. 

  833.             tem       = PHIL(I,k) + CP * TO(I,k)
    834. 

  834.             HEO(I,k)  = tem  + HVAP * QO(I,k)
    835. 

  835.             HESO(I,k) = tem  + HVAP * QESO(I,k)
    836. 

  836. !           HEO(I,k)  = MIN(HEO(I,k),HESO(I,k))
    837. 

  837.           endif
    838. 

  838.         ENDDO
    839. 

  839.       ENDDO
    840. 

  840. !
    841. 

  841. !  DETERMINE LEVEL WITH LARGEST MOIST STATIC ENERGY
    842. 

  842. !  THIS IS THE LEVEL WHERE UPDRAFT STARTS
    843. 

  843. !
    844. 

  844.       DO I=1,IM
    845. 

  845.         HMAX(I) = HEO(I,1)
    846. 

  846.         KB(I) = 1
    847. 

  847.       ENDDO
    848. 

  848. !!
    849. 

  849.       DO K = 2, KM
    850. 

  850.         DO I=1,IM
    851. 

  851.           if (k .le. kbm(i)) then
    852. 

  852.             IF(HEO(I,k).GT.HMAX(I).AND.CNVFLG(I)) THEN
    853. 

  853.               KB(I)   = K
    854. 

  854.               HMAX(I) = HEO(I,k)
    855. 

  855.             ENDIF
    856. 

  856.           endif
    857. 

  857.         ENDDO
    858. 

  858.       ENDDO
    859. 

  859. !     DO K = 1, KMAX - 1
    860. 

  860. !         TOL(k) = .5 * (TO(I,k) + TO(I,k+1))
    861. 

  861. !         QOL(k) = .5 * (QO(I,k) + QO(I,k+1))
    862. 

  862. !         QESOL(I,k) = .5 * (QESO(I,k) + QESO(I,k+1))
    863. 

  863. !         HEOL(I,k) = .5 * (HEO(I,k) + HEO(I,k+1))
    864. 

  864. !         HESOL(I,k) = .5 * (HESO(I,k) + HESO(I,k+1))
    865. 

  865. !     ENDDO
    866. 

  866.       DO K = 1, KM1
    867. 

  867.         DO I=1,IM
    868. 

  868.           if (k .le. kmax(i)-1) then
    869. 

  869.             DZ      = .5 * (ZO(I,k+1) - ZO(I,k))
    870. 

  870.             DP      = .5 * (PFLD(I,k+1) - PFLD(I,k))
    871. 

  871. !jfe        ES      = 10. * FPVS(TO(I,k+1))
    872. 

  872. !
    873. 

  873.             ES      = 0.01 * fpvs(TO(I,K+1))      ! fpvs is in Pa
    874. 

  874. !
    875. 

  875.             PPRIME  = PFLD(I,k+1) + EPSM1 * ES
    876. 

  876.             QS      = EPS * ES / PPRIME
    877. 

  877.             DQSDP   = - QS / PPRIME
    878. 

  878.             DESDT   = ES * (FACT1 / TO(I,k+1) + FACT2 / (TO(I,k+1)**2))
    879. 

  879.             DQSDT   = QS * PFLD(I,k+1) * DESDT / (ES * PPRIME)
    880. 

  880.             GAMMA   = EL2ORC * QESO(I,k+1) / (TO(I,k+1)**2)
    881. 

  881.             DT      = (G * DZ + HVAP * DQSDP * DP) / (CP * (1. + GAMMA))
    882. 

  882.             DQ      = DQSDT * DT + DQSDP * DP
    883. 

  883.             TO(I,k) = TO(I,k+1) + DT
    884. 

  884.             QO(I,k) = QO(I,k+1) + DQ
    885. 

  885.             PO(I,k) = .5 * (PFLD(I,k) + PFLD(I,k+1))
    886. 

  886.           endif
    887. 

  887.         ENDDO
    888. 

  888.       ENDDO
    889. 

  889. !
    890. 

  890.       DO K = 1, KM1
    891. 

  891.         DO I=1,IM
    892. 

  892.           if (k .le. kmax(I)-1) then
    893. 

  893. !jfe        QESO(I,k) = 10. * FPVS(TO(I,k))
    894. 

  894. !
    895. 

  895.             QESO(I,k) = 0.01 * fpvs(TO(I,K))      ! fpvs is in Pa
    896. 

  896. !
    897. 

  897.             QESO(I,k) = EPS * QESO(I,k) / (PO(I,k) + EPSM1*QESO(I,k))
    898. 

  898. !cmr        QESO(I,k) = MAX(QESO(I,k),1.E-8)
    899. 

  899.             val1      =             1.E-8
    900. 

  900.             QESO(I,k) = MAX(QESO(I,k), val1)
    901. 

  901. !cmr        QO(I,k)   = max(QO(I,k),1.e-10)
    902. 

  902.             val2      =           1.e-10
    903. 

  903.             QO(I,k)   = max(QO(I,k), val2 )
    904. 

  904. !           QO(I,k)   = MIN(QO(I,k),QESO(I,k))
    905. 

  905.             HEO(I,k)  = .5 * G * (ZO(I,k) + ZO(I,k+1)) +                &
    906. 

  906.      &                  CP * TO(I,k) + HVAP * QO(I,k)
    907. 

  907.             HESO(I,k) = .5 * G * (ZO(I,k) + ZO(I,k+1)) +                & 
    908. 

  908.      &                  CP * TO(I,k) + HVAP * QESO(I,k)
    909. 

  909.             UO(I,k)   = .5 * (UO(I,k) + UO(I,k+1))
    910. 

  910.             VO(I,k)   = .5 * (VO(I,k) + VO(I,k+1))
    911. 

  911.           endif
    912. 

  912.         ENDDO
    913. 

  913.       ENDDO
    914. 

  914. !     k = kmax
    915. 

  915. !       HEO(I,k) = HEO(I,k)
    916. 

  916. !       hesol(k) = HESO(I,k)
    917. 

  917. !      IF(LAT.EQ.LATD.AND.lon.eq.lond.and.CNVFLG(I)) THEN
    918. 

  918. !        PRINT *, '   HEO ='
    919. 

  919. !        PRINT 6001, (HEO(I,K),K=1,KMAX)
    920. 

  920. !        PRINT *, '   HESO ='
    921. 

  921. !        PRINT 6001, (HESO(I,K),K=1,KMAX)
    922. 

  922. !        PRINT *, '   TO ='
    923. 

  923. !        PRINT 6002, (TO(I,K)-273.16,K=1,KMAX)
    924. 

  924. !        PRINT *, '   QO ='
    925. 

  925. !        PRINT 6003, (QO(I,K),K=1,KMAX)
    926. 

  926. !        PRINT *, '   QSO ='
    927. 

  927. !        PRINT 6003, (QESO(I,K),K=1,KMAX)
    928. 

  928. !      ENDIF
    929. 

  929. !
    930. 

  930. !  LOOK FOR CONVECTIVE CLOUD BASE AS THE LEVEL OF FREE CONVECTION
    931. 

  931. !
    932. 

  932.       DO I=1,IM
    933. 

  933.         IF(CNVFLG(I)) THEN
    934. 

  934.           INDX    = KB(I)
    935. 

  935.           HKBO(I) = HEO(I,INDX)
    936. 

  936.           QKBO(I) = QO(I,INDX)
    937. 

  937.           UKBO(I) = UO(I,INDX)
    938. 

  938.           VKBO(I) = VO(I,INDX)
    939. 

  939.         ENDIF
    940. 

  940.         FLG(I)    = CNVFLG(I)
    941. 

  941.         KBCON(I)  = KMAX(I)
    942. 

  942.       ENDDO
    943. 

  943. !!
    944. 

  944.       DO K = 1, KM
    945. 

  945.         DO I=1,IM
    946. 

  946.           if (k .le. kbmax(i)) then
    947. 

  947.             IF(FLG(I).AND.K.GT.KB(I)) THEN
    948. 

  948.               HSBAR(I)   = HESO(I,k)
    949. 

  949.               IF(HKBO(I).GT.HSBAR(I)) THEN
    950. 

  950.                 FLG(I)   = .FALSE.
    951. 

  951.                 KBCON(I) = K
    952. 

  952.               ENDIF
    953. 

  953.             ENDIF
    954. 

  954.           endif
    955. 

  955.         ENDDO
    956. 

  956.       ENDDO
    957. 

  957.       DO I=1,IM
    958. 

  958.         IF(CNVFLG(I)) THEN
    959. 

  959.           PBCDIF(I) = -PFLD(I,KBCON(I)) + PFLD(I,KB(I))
    960. 

  960.           PDOT(I)   = 10.* DOT(I,KBCON(I))
    961. 

  961.           IF(PBCDIF(I).GT.150.)    CNVFLG(I) = .FALSE.
    962. 

  962.           IF(KBCON(I).EQ.KMAX(I))  CNVFLG(I) = .FALSE.
    963. 

  963.         ENDIF
    964. 

  964.       ENDDO
    965. 

  965. !!
    966. 

  966.       TOTFLG = .TRUE.
    967. 

  967.       DO I=1,IM
    968. 

  968.         TOTFLG = TOTFLG .AND. (.NOT. CNVFLG(I))
    969. 

  969.       ENDDO
    970. 

  970.       IF(TOTFLG) RETURN
    971. 

  971. !  FOUND LFC, CAN DEFINE REST OF VARIABLES
    972. 

  972.  6001 FORMAT(2X,-2P10F12.2)
    973. 

  973.  6002 FORMAT(2X,10F12.2)
    974. 

  974.  6003 FORMAT(2X,3P10F12.2)
    975. 

  975. 

  976. !
    977. 

  977. !  DETERMINE ENTRAINMENT RATE BETWEEN KB AND KBCON
    978. 

  978. !
    979. 

  979.       DO I = 1, IM
    980. 

  980.         alpha = alphas
    981. 

  981.         if(SLIMSK(I).eq.1.) alpha = alphal
    982. 

  982.         IF(CNVFLG(I)) THEN
    983. 

  983.           IF(KB(I).EQ.1) THEN
    984. 

  984.             DZ = .5 * (ZO(I,KBCON(I)) + ZO(I,KBCON(I)-1)) - ZO(I,1)
    985. 

  985.           ELSE
    986. 

  986.             DZ = .5 * (ZO(I,KBCON(I)) + ZO(I,KBCON(I)-1))               &
    987. 

  987.      &         - .5 * (ZO(I,KB(I)) + ZO(I,KB(I)-1))
    988. 

  988.           ENDIF
    989. 

  989.           IF(KBCON(I).NE.KB(I)) THEN
    990. 

  990. !cmr        XLAMB(I) = -ALOG(ALPHA) / DZ
    991. 

  991.             XLAMB(I) = - LOG(ALPHA) / DZ
    992. 

  992.           ELSE
    993. 

  993.             XLAMB(I) = 0.
    994. 

  994.           ENDIF
    995. 

  995.         ENDIF
    996. 

  996.       ENDDO
    997. 

  997. !  DETERMINE UPDRAFT MASS FLUX
    998. 

  998.       DO K = 1, KM
    999. 

  999.         DO I = 1, IM
    1000. 

  1000.           if (k .le. kmax(i) .and. CNVFLG(I)) then
    1001. 

  1001.             ETA(I,k)  = 1.
    1002. 

  1002.             ETAU(I,k) = 1.
    1003. 

  1003.           ENDIF
    1004. 

  1004.         ENDDO
    1005. 

  1005.       ENDDO
    1006. 

  1006.       DO K = KM1, 2, -1
    1007. 

  1007.         DO I = 1, IM
    1008. 

  1008.           if (k .le. kbmax(i)) then
    1009. 

  1009.             IF(CNVFLG(I).AND.K.LT.KBCON(I).AND.K.GE.KB(I)) THEN
    1010. 

  1010.               DZ        = .5 * (ZO(I,k+1) - ZO(I,k-1))
    1011. 

  1011.               ETA(I,k)  = ETA(I,k+1) * EXP(-XLAMB(I) * DZ)
    1012. 

  1012.               ETAU(I,k) = ETA(I,k)
    1013. 

  1013.             ENDIF
    1014. 

  1014.           endif
    1015. 

  1015.         ENDDO
    1016. 

  1016.       ENDDO
    1017. 

  1017.       DO I = 1, IM
    1018. 

  1018.         IF(CNVFLG(I).AND.KB(I).EQ.1.AND.KBCON(I).GT.1) THEN
    1019. 

  1019.           DZ = .5 * (ZO(I,2) - ZO(I,1))
    1020. 

  1020.           ETA(I,1) = ETA(I,2) * EXP(-XLAMB(I) * DZ)
    1021. 

  1021.           ETAU(I,1) = ETA(I,1)
    1022. 

  1022.         ENDIF
    1023. 

  1023.       ENDDO
    1024. 

  1024. !
    1025. 

  1025. !  WORK UP UPDRAFT CLOUD PROPERTIES
    1026. 

  1026. !
    1027. 

  1027.       DO I = 1, IM
    1028. 

  1028.         IF(CNVFLG(I)) THEN
    1029. 

  1029.           INDX         = KB(I)
    1030. 

  1030.           HCKO(I,INDX) = HKBO(I)
    1031. 

  1031.           QCKO(I,INDX) = QKBO(I)
    1032. 

  1032.           UCKO(I,INDX) = UKBO(I)
    1033. 

  1033.           VCKO(I,INDX) = VKBO(I)
    1034. 

  1034.           PWAVO(I)     = 0.
    1035. 

  1035.         ENDIF
    1036. 

  1036.       ENDDO
    1037. 

  1037. !
    1038. 

  1038. !  CLOUD PROPERTY BELOW CLOUD BASE IS MODIFIED BY THE ENTRAINMENT PROCES
    1039. 

  1039. !
    1040. 

  1040.       DO K = 2, KM1
    1041. 

  1041.         DO I = 1, IM
    1042. 

  1042.           if (k .le. kmax(i)-1) then
    1043. 

  1043.             IF(CNVFLG(I).AND.K.GT.KB(I).AND.K.LE.KBCON(I)) THEN
    1044. 

  1044.               FACTOR = ETA(I,k-1) / ETA(I,k)
    1045. 

  1045.               ONEMF = 1. - FACTOR
    1046. 

  1046.               HCKO(I,k) = FACTOR * HCKO(I,k-1) + ONEMF *                &
    1047. 

  1047.      &                    .5 * (HEO(I,k) + HEO(I,k+1))
    1048. 

  1048.               UCKO(I,k) = FACTOR * UCKO(I,k-1) + ONEMF *                & 
    1049. 

  1049.      &                    .5 * (UO(I,k) + UO(I,k+1))
    1050. 

  1050.               VCKO(I,k) = FACTOR * VCKO(I,k-1) + ONEMF *                &
    1051. 

  1051.      &                    .5 * (VO(I,k) + VO(I,k+1))
    1052. 

  1052.               DBYO(I,k) = HCKO(I,k) - HESO(I,k)
    1053. 

  1053.             ENDIF
    1054. 

  1054.             IF(CNVFLG(I).AND.K.GT.KBCON(I)) THEN
    1055. 

  1055.               HCKO(I,k) = HCKO(I,k-1)
    1056. 

  1056.               UCKO(I,k) = UCKO(I,k-1)
    1057. 

  1057.               VCKO(I,k) = VCKO(I,k-1)
    1058. 

  1058.               DBYO(I,k) = HCKO(I,k) - HESO(I,k)
    1059. 

  1059.             ENDIF
    1060. 

  1060.           endif
    1061. 

  1061.         ENDDO
    1062. 

  1062.       ENDDO
    1063. 

  1063. !  DETERMINE CLOUD TOP
    1064. 

  1064.       DO I = 1, IM
    1065. 

  1065.         FLG(I) = CNVFLG(I)
    1066. 

  1066.         KTCON(I) = 1
    1067. 

  1067.       ENDDO
    1068. 

  1068. !     DO K = 2, KMAX
    1069. 

  1069. !       KK = KMAX - K + 1
    1070. 

  1070. !         IF(DBYO(I,kK).GE.0..AND.FLG(I).AND.KK.GT.KBCON(I)) THEN
    1071. 

  1071. !           KTCON(I) = KK + 1
    1072. 

  1072. !           FLG(I) = .FALSE.
    1073. 

  1073. !         ENDIF
    1074. 

  1074. !     ENDDO
    1075. 

  1075.       DO K = 2, KM
    1076. 

  1076.         DO I = 1, IM
    1077. 

  1077.           if (k .le. kmax(i)) then
    1078. 

  1078.             IF(DBYO(I,k).LT.0..AND.FLG(I).AND.K.GT.KBCON(I)) THEN
    1079. 

  1079.               KTCON(I) = K
    1080. 

  1080.               FLG(I) = .FALSE.
    1081. 

  1081.             ENDIF
    1082. 

  1082.           endif
    1083. 

  1083.         ENDDO
    1084. 

  1084.       ENDDO
    1085. 

  1085.       DO I = 1, IM
    1086. 

  1086.         IF(CNVFLG(I).AND.(PFLD(I,KBCON(I)) - PFLD(I,KTCON(I))).LT.150.) &
    1087. 

  1087.      &  CNVFLG(I) = .FALSE.
    1088. 

  1088.       ENDDO
    1089. 

  1089.       TOTFLG = .TRUE.
    1090. 

  1090.       DO I = 1, IM
    1091. 

  1091.         TOTFLG = TOTFLG .AND. (.NOT. CNVFLG(I))
    1092. 

  1092.       ENDDO
    1093. 

  1093.       IF(TOTFLG) RETURN
    1094. 

  1094. !
    1095. 

  1095. !  SEARCH FOR DOWNDRAFT ORIGINATING LEVEL ABOVE THETA-E MINIMUM
    1096. 

  1096. !
    1097. 

  1097.       DO I = 1, IM
    1098. 

  1098.         HMIN(I) = HEO(I,KBCON(I))
    1099. 

  1099.         LMIN(I) = KBMAX(I)
    1100. 

  1100.         JMIN(I) = KBMAX(I)
    1101. 

  1101.       ENDDO
    1102. 

  1102.       DO I = 1, IM
    1103. 

  1103.         DO K = KBCON(I), KBMAX(I)
    1104. 

  1104.           IF(HEO(I,k).LT.HMIN(I).AND.CNVFLG(I)) THEN
    1105. 

  1105.             LMIN(I) = K + 1
    1106. 

  1106.             HMIN(I) = HEO(I,k)
    1107. 

  1107.           ENDIF
    1108. 

  1108.         ENDDO
    1109. 

  1109.       ENDDO
    1110. 

  1110. !
    1111. 

  1111. !  Make sure that JMIN(I) is within the cloud
    1112. 

  1112. !
    1113. 

  1113.       DO I = 1, IM
    1114. 

  1114.         IF(CNVFLG(I)) THEN
    1115. 

  1115.           JMIN(I) = MIN(LMIN(I),KTCON(I)-1)
    1116. 

  1116.           XMBMAX(I) = .1
    1117. 

  1117.           JMIN(I) = MAX(JMIN(I),KBCON(I)+1)
    1118. 

  1118.         ENDIF
    1119. 

  1119.       ENDDO
    1120. 

  1120. !
    1121. 

  1121. !  ENTRAINING CLOUD
    1122. 

  1122. !
    1123. 

  1123.       do k = 2, km1
    1124. 

  1124.         DO I = 1, IM
    1125. 

  1125.           if (k .le. kmax(i)-1) then
    1126. 

  1126.             if(CNVFLG(I).and.k.gt.JMIN(I).and.k.le.KTCON(I)) THEN
    1127. 

  1127.               SUMZ(I,k) = SUMZ(I,k-1) + .5 * (ZO(I,k+1) - ZO(I,k-1))
    1128. 

  1128.               SUMH(I,k) = SUMH(I,k-1) + .5 * (ZO(I,k+1) - ZO(I,k-1))    &
    1129. 

  1129.      &                  * HEO(I,k)
    1130. 

  1130.             ENDIF
    1131. 

  1131.           endif
    1132. 

  1132.         enddo
    1133. 

  1133.       enddo
    1134. 

  1134. !!
    1135. 

  1135.       DO I = 1, IM
    1136. 

  1136.         IF(CNVFLG(I)) THEN
    1137. 

  1137. !         call random_number(XKT2)
    1138. 

  1138. !         call srand(fhour)
    1139. 

  1139. !         XKT2(I) = rand()
    1140. 

  1140.           KT2(I) = nint(XKT2(I)*float(KTCON(I)-JMIN(I))-.5)+JMIN(I)+1
    1141. 

  1141. !         KT2(I) = nint(sqrt(XKT2(I))*float(KTCON(I)-JMIN(I))-.5) + JMIN(I) + 1
    1142. 

  1142. !         KT2(I) = nint(ranf() *float(KTCON(I)-JMIN(I))-.5) + JMIN(I) + 1
    1143. 

  1143.           tem1 = (HCKO(I,JMIN(I)) - HESO(I,KT2(I)))
    1144. 

  1144.           tem2 = (SUMZ(I,KT2(I)) * HESO(I,KT2(I)) - SUMH(I,KT2(I)))
    1145. 

  1145.           if (abs(tem2) .gt. 0.000001) THEN
    1146. 

  1146.             XLAMB(I) = tem1 / tem2
    1147. 

  1147.           else
    1148. 

  1148.             CNVFLG(I) = .false.
    1149. 

  1149.           ENDIF
    1150. 

  1150. !         XLAMB(I) = (HCKO(I,JMIN(I)) - HESO(I,KT2(I)))
    1151. 

  1151. !    &          / (SUMZ(I,KT2(I)) * HESO(I,KT2(I)) - SUMH(I,KT2(I)))
    1152. 

  1152.           XLAMB(I) = max(XLAMB(I),RZERO)
    1153. 

  1153.           XLAMB(I) = min(XLAMB(I),2.3/SUMZ(I,KT2(I)))
    1154. 

  1154.         ENDIF
    1155. 

  1155.       ENDDO
    1156. 

  1156. !!
    1157. 

  1157.       DO I = 1, IM
    1158. 

  1158.        DWNFLG(I)  = CNVFLG(I)
    1159. 

  1159.        DWNFLG2(I) = CNVFLG(I)
    1160. 

  1160.        IF(CNVFLG(I)) THEN
    1161. 

  1161.         if(KT2(I).ge.KTCON(I)) DWNFLG(I) = .false.
    1162. 

  1162.       if(XLAMB(I).le.1.e-30.or.HCKO(I,JMIN(I))-HESO(I,KT2(I)).le.1.e-30)&
    1163. 

  1163.      &  DWNFLG(I) = .false.
    1164. 

  1164.         do k = JMIN(I), KT2(I)
    1165. 

  1165.           if(DWNFLG(I).and.HEO(I,k).gt.HESO(I,KT2(I))) DWNFLG(I)=.false.
    1166. 

  1166.         enddo
    1167. 

  1167. !       IF(CNVFLG(I).AND.(PFLD(KBCON(I))-PFLD(KTCON(I))).GT.PDETRN)
    1168. 

  1168. !    &     DWNFLG(I)=.FALSE.
    1169. 

  1169.         IF(CNVFLG(I).AND.(PFLD(I,KBCON(I))-PFLD(I,KTCON(I))).LT.PDPDWN) &
    1170. 

  1170.      &     DWNFLG2(I)=.FALSE.
    1171. 

  1171.        ENDIF
    1172. 

  1172.       ENDDO
    1173. 

  1173. !!
    1174. 

  1174.       DO K = 2, KM1
    1175. 

  1175.         DO I = 1, IM
    1176. 

  1176.           if (k .le. kmax(i)-1) then
    1177. 

  1177.             IF(DWNFLG(I).AND.K.GT.JMIN(I).AND.K.LE.KT2(I)) THEN
    1178. 

  1178.               DZ        = .5 * (ZO(I,k+1) - ZO(I,k-1))
    1179. 

  1179. !             ETA(I,k)  = ETA(I,k-1) * EXP( XLAMB(I) * DZ)
    1180. 

  1180. !  to simplify matter, we will take the linear approach here
    1181. 

  1181. !
    1182. 

  1182.               ETA(I,k)  = ETA(I,k-1) * (1. + XLAMB(I) * dz)
    1183. 

  1183.               ETAU(I,k) = ETAU(I,k-1) * (1. + (XLAMB(I)+xlambu) * dz)
    1184. 

  1184.             ENDIF
    1185. 

  1185.           endif
    1186. 

  1186.         ENDDO
    1187. 

  1187.       ENDDO
    1188. 

  1188. !!
    1189. 

  1189.       DO K = 2, KM1
    1190. 

  1190.         DO I = 1, IM
    1191. 

  1191.           if (k .le. kmax(i)-1) then
    1192. 

  1192. !           IF(.NOT.DWNFLG(I).AND.K.GT.JMIN(I).AND.K.LE.KT2(I)) THEN
    1193. 

  1193.             IF(.NOT.DWNFLG(I).AND.K.GT.JMIN(I).AND.K.LE.KTCON(I)) THEN
    1194. 

  1194.               DZ        = .5 * (ZO(I,k+1) - ZO(I,k-1))
    1195. 

  1195.               ETAU(I,k) = ETAU(I,k-1) * (1. + xlambu * dz)
    1196. 

  1196.             ENDIF
    1197. 

  1197.           endif
    1198. 

  1198.         ENDDO
    1199. 

  1199.       ENDDO
    1200. 

  1200. !      IF(LAT.EQ.LATD.AND.lon.eq.lond.and.CNVFLG(I)) THEN
    1201. 

  1201. !        PRINT *, ' LMIN(I), KT2(I)=', LMIN(I), KT2(I)
    1202. 

  1202. !        PRINT *, ' KBOT, KTOP, JMIN(I) =', KBCON(I), KTCON(I), JMIN(I)
    1203. 

  1203. !      ENDIF
    1204. 

  1204. !     IF(LAT.EQ.LATD.AND.lon.eq.lond) THEN
    1205. 

  1205. !       print *, ' xlamb =', xlamb
    1206. 

  1206. !       print *, ' eta =', (eta(k),k=1,KT2(I))
    1207. 

  1207. !       print *, ' ETAU =', (ETAU(I,k),k=1,KT2(I))
    1208. 

  1208. !       print *, ' HCKO =', (HCKO(I,k),k=1,KT2(I))
    1209. 

  1209. !       print *, ' SUMZ =', (SUMZ(I,k),k=1,KT2(I))
    1210. 

  1210. !       print *, ' SUMH =', (SUMH(I,k),k=1,KT2(I))
    1211. 

  1211. !     ENDIF
    1212. 

  1212.       DO I = 1, IM
    1213. 

  1213.         if(DWNFLG(I)) THEN
    1214. 

  1214.           KTCON(I) = KT2(I)
    1215. 

  1215.         ENDIF
    1216. 

  1216.       ENDDO
    1217. 

  1217. !
    1218. 

  1218. !  CLOUD PROPERTY ABOVE CLOUD Base IS MODIFIED BY THE DETRAINMENT PROCESS
    1219. 

  1219. !
    1220. 

  1220.       DO K = 2, KM1
    1221. 

  1221.         DO I = 1, IM
    1222. 

  1222.           if (k .le. kmax(i)-1) then
    1223. 

  1223. !jfe
    1224. 

  1224.             IF(CNVFLG(I).AND.K.GT.KBCON(I).AND.K.LE.KTCON(I)) THEN
    1225. 

  1225. !jfe      IF(K.GT.KBCON(I).AND.K.LE.KTCON(I)) THEN
    1226. 

  1226.               FACTOR    = ETA(I,k-1) / ETA(I,k)
    1227. 

  1227.               ONEMF     = 1. - FACTOR
    1228. 

  1228.               fuv       = ETAU(I,k-1) / ETAU(I,k)
    1229. 

  1229.               onemfu    = 1. - fuv
    1230. 

  1230.               HCKO(I,k) = FACTOR * HCKO(I,k-1) + ONEMF *                &
    1231. 

  1231.      &                    .5 * (HEO(I,k) + HEO(I,k+1))
    1232. 

  1232.               UCKO(I,k) = fuv * UCKO(I,k-1) + ONEMFu *                  &
    1233. 

  1233.      &                    .5 * (UO(I,k) + UO(I,k+1))
    1234. 

  1234.               VCKO(I,k) = fuv * VCKO(I,k-1) + ONEMFu *                  &
    1235. 

  1235.      &                    .5 * (VO(I,k) + VO(I,k+1))
    1236. 

  1236.               DBYO(I,k) = HCKO(I,k) - HESO(I,k)
    1237. 

  1237.             ENDIF
    1238. 

  1238.           endif
    1239. 

  1239.         ENDDO
    1240. 

  1240.       ENDDO
    1241. 

  1241. !      IF(LAT.EQ.LATD.AND.lon.eq.lond.and.CNVFLG(I)) THEN
    1242. 

  1242. !        PRINT *, ' UCKO=', (UCKO(I,k),k=KBCON(I)+1,KTCON(I))
    1243. 

  1243. !        PRINT *, ' uenv=', (.5*(UO(I,k)+UO(I,k-1)),k=KBCON(I)+1,KTCON(I))
    1244. 

  1244. !      ENDIF
    1245. 

  1245.       DO I = 1, IM
    1246. 

  1246.         if(CNVFLG(I).and.DWNFLG2(I).and.JMIN(I).le.KBCON(I))            &
    1247. 

  1247.      &     THEN
    1248. 

  1248.           CNVFLG(I) = .false.
    1249. 

  1249.           DWNFLG(I) = .false.
    1250. 

  1250.           DWNFLG2(I) = .false.
    1251. 

  1251.         ENDIF
    1252. 

  1252.       ENDDO
    1253. 

  1253. !!
    1254. 

  1254.       TOTFLG = .TRUE.
    1255. 

  1255.       DO I = 1, IM
    1256. 

  1256.         TOTFLG = TOTFLG .AND. (.NOT. CNVFLG(I))
    1257. 

  1257.       ENDDO
    1258. 

  1258.       IF(TOTFLG) RETURN
    1259. 

  1259. !!
    1260. 

  1260. !
    1261. 

  1261. !  COMPUTE CLOUD MOISTURE PROPERTY AND PRECIPITATION
    1262. 

  1262. !
    1263. 

  1263.       DO I = 1, IM
    1264. 

  1264.           AA1(I) = 0.
    1265. 

  1265.           RHBAR(I) = 0.
    1266. 

  1266.       ENDDO
    1267. 

  1267.       DO K = 1, KM
    1268. 

  1268.         DO I = 1, IM
    1269. 

  1269.           if (k .le. kmax(i)) then
    1270. 

  1270.             IF(CNVFLG(I).AND.K.GT.KB(I).AND.K.LT.KTCON(I)) THEN
    1271. 

  1271.               DZ = .5 * (ZO(I,k+1) - ZO(I,k-1))
    1272. 

  1272.               DZ1 = (ZO(I,k) - ZO(I,k-1))
    1273. 

  1273.               GAMMA = EL2ORC * QESO(I,k) / (TO(I,k)**2)
    1274. 

  1274.               QRCH = QESO(I,k)                                          &
    1275. 

  1275.      &             + GAMMA * DBYO(I,k) / (HVAP * (1. + GAMMA))
    1276. 

  1276.               FACTOR = ETA(I,k-1) / ETA(I,k)
    1277. 

  1277.               ONEMF = 1. - FACTOR
    1278. 

  1278.               QCKO(I,k) = FACTOR * QCKO(I,k-1) + ONEMF *                &
    1279. 

  1279.      &                    .5 * (QO(I,k) + QO(I,k+1))
    1280. 

  1280.               DQ = ETA(I,k) * QCKO(I,k) - ETA(I,k) * QRCH
    1281. 

  1281.               RHBAR(I) = RHBAR(I) + QO(I,k) / QESO(I,k)
    1282. 

  1282. !
    1283. 

  1283. !  BELOW LFC CHECK IF THERE IS EXCESS MOISTURE TO RELEASE LATENT HEAT
    1284. 

  1284. !
    1285. 

  1285.               IF(DQ.GT.0.) THEN
    1286. 

  1286.                 ETAH = .5 * (ETA(I,k) + ETA(I,k-1))
    1287. 

  1287.                 QLK = DQ / (ETA(I,k) + ETAH * C0 * DZ)
    1288. 

  1288.                 AA1(I) = AA1(I) - DZ1 * G * QLK
    1289. 

  1289.                 QC = QLK + QRCH
    1290. 

  1290.                 PWO(I,k) = ETAH * C0 * DZ * QLK
    1291. 

  1291.                 QCKO(I,k) = QC
    1292. 

  1292.                 PWAVO(I) = PWAVO(I) + PWO(I,k)
    1293. 

  1293.               ENDIF
    1294. 

  1294.             ENDIF
    1295. 

  1295.           endif
    1296. 

  1296.         ENDDO
    1297. 

  1297.       ENDDO
    1298. 

  1298.       DO I = 1, IM
    1299. 

  1299.         RHBAR(I) = RHBAR(I) / float(KTCON(I) - KB(I) - 1)
    1300. 

  1300.       ENDDO
    1301. 

  1301. !
    1302. 

  1302. !  this section is ready for cloud water
    1303. 

  1303. !
    1304. 

  1304.       if(ncloud.gt.0) THEN
    1305. 

  1305. !
    1306. 

  1306. !  compute liquid and vapor separation at cloud top
    1307. 

  1307. !
    1308. 

  1308.       DO I = 1, IM
    1309. 

  1309.         k = KTCON(I)
    1310. 

  1310.         IF(CNVFLG(I)) THEN
    1311. 

  1311.           GAMMA = EL2ORC * QESO(I,K) / (TO(I,K)**2)
    1312. 

  1312.           QRCH = QESO(I,K)                                              &
    1313. 

  1313.      &         + GAMMA * DBYO(I,K) / (HVAP * (1. + GAMMA))
    1314. 

  1314.           DQ = QCKO(I,K-1) - QRCH
    1315. 

  1315. !
    1316. 

  1316. !  CHECK IF THERE IS EXCESS MOISTURE TO RELEASE LATENT HEAT
    1317. 

  1317. !
    1318. 

  1318.           IF(DQ.GT.0.) THEN
    1319. 

  1319.             QLKO_KTCON(I) = dq
    1320. 

  1320.             QCKO(I,K-1) = QRCH
    1321. 

  1321.           ENDIF
    1322. 

  1322.         ENDIF
    1323. 

  1323.       ENDDO
    1324. 

  1324.       ENDIF
    1325. 

  1325. !
    1326. 

  1326. !  CALCULATE CLOUD WORK FUNCTION AT T+DT
    1327. 

  1327. !
    1328. 

  1328.       DO K = 1, KM
    1329. 

  1329.         DO I = 1, IM
    1330. 

  1330.           if (k .le. kmax(i)) then
    1331. 

  1331.             IF(CNVFLG(I).AND.K.GT.KBCON(I).AND.K.LE.KTCON(I)) THEN
    1332. 

  1332.               DZ1 = ZO(I,k) - ZO(I,k-1)
    1333. 

  1333.               GAMMA = EL2ORC * QESO(I,k-1) / (TO(I,k-1)**2)
    1334. 

  1334.               RFACT =  1. + DELTA * CP * GAMMA                          &
    1335. 

  1335.      &                 * TO(I,k-1) / HVAP
    1336. 

  1336.               AA1(I) = AA1(I) +                                         &
    1337. 

  1337.      &                 DZ1 * (G / (CP * TO(I,k-1)))                     &
    1338. 

  1338.      &                 * DBYO(I,k-1) / (1. + GAMMA)                     &
    1339. 

  1339.      &                 * RFACT
    1340. 

  1340.               val = 0.
    1341. 

  1341.               AA1(I)=AA1(I)+                                            &
    1342. 

  1342.      &                 DZ1 * G * DELTA *                                &
    1343. 

  1343. !cmr &                 MAX( 0.,(QESO(I,k-1) - QO(I,k-1)))               & 
    1344. 

  1344.      &                 MAX(val,(QESO(I,k-1) - QO(I,k-1)))
    1345. 

  1345.             ENDIF
    1346. 

  1346.           endif
    1347. 

  1347.         ENDDO
    1348. 

  1348.       ENDDO
    1349. 

  1349.       DO I = 1, IM
    1350. 

  1350.         IF(CNVFLG(I).AND.AA1(I).LE.0.) DWNFLG(I)  = .FALSE.
    1351. 

  1351.         IF(CNVFLG(I).AND.AA1(I).LE.0.) DWNFLG2(I) = .FALSE.
    1352. 

  1352.         IF(CNVFLG(I).AND.AA1(I).LE.0.) CNVFLG(I)  = .FALSE.
    1353. 

  1353.       ENDDO
    1354. 

  1354. !!
    1355. 

  1355.       TOTFLG = .TRUE.
    1356. 

  1356.       DO I = 1, IM
    1357. 

  1357.         TOTFLG = TOTFLG .AND. (.NOT. CNVFLG(I))
    1358. 

  1358.       ENDDO
    1359. 

  1359.       IF(TOTFLG) RETURN
    1360. 

  1360. !!
    1361. 

  1361. !cccc IF(LAT.EQ.LATD.AND.lon.eq.lond.and.CNVFLG(I)) THEN
    1362. 

  1362. !cccc   PRINT *, ' AA1(I) BEFORE DWNDRFT =', AA1(I)
    1363. 

  1363. !cccc ENDIF
    1364. 

  1364. !
    1365. 

  1365. !------- DOWNDRAFT CALCULATIONS
    1366. 

  1366. !
    1367. 

  1367. !
    1368. 

  1368. !--- DETERMINE DOWNDRAFT STRENGTH IN TERMS OF WINDSHEAR
    1369. 

  1369. !
    1370. 

  1370.       DO I = 1, IM
    1371. 

  1371.         IF(CNVFLG(I)) THEN
    1372. 

  1372.           VSHEAR(I) = 0.
    1373. 

  1373.         ENDIF
    1374. 

  1374.       ENDDO
    1375. 

  1375.       DO K = 1, KM
    1376. 

  1376.         DO I = 1, IM
    1377. 

  1377.           if (k .le. kmax(i)) then
    1378. 

  1378.             IF(K.GE.KB(I).AND.K.LE.KTCON(I).AND.CNVFLG(I)) THEN
    1379. 

  1379.               shear=rcs(I) * sqrt((UO(I,k+1)-UO(I,k)) ** 2              &
    1380. 

  1380.      &                          + (VO(I,k+1)-VO(I,k)) ** 2)
    1381. 

  1381.               VSHEAR(I) = VSHEAR(I) + SHEAR
    1382. 

  1382.             ENDIF
    1383. 

  1383.           endif
    1384. 

  1384.         ENDDO
    1385. 

  1385.       ENDDO
    1386. 

  1386.       DO I = 1, IM
    1387. 

  1387.         EDT(I) = 0.
    1388. 

  1388.         IF(CNVFLG(I)) THEN
    1389. 

  1389.           KNUMB = KTCON(I) - KB(I) + 1
    1390. 

  1390.           KNUMB = MAX(KNUMB,1)
    1391. 

  1391.           VSHEAR(I) = 1.E3 * VSHEAR(I) / (ZO(I,KTCON(I))-ZO(I,KB(I)))
    1392. 

  1392.           E1=1.591-.639*VSHEAR(I)                                       &
    1393. 

  1393.      &       +.0953*(VSHEAR(I)**2)-.00496*(VSHEAR(I)**3)
    1394. 

  1394.           EDT(I)=1.-E1
    1395. 

  1395. !cmr      EDT(I) = MIN(EDT(I),.9)
    1396. 

  1396.           val =         .9
    1397. 

  1397.           EDT(I) = MIN(EDT(I),val)
    1398. 

  1398. !cmr      EDT(I) = MAX(EDT(I),.0)
    1399. 

  1399.           val =         .0
    1400. 

  1400.           EDT(I) = MAX(EDT(I),val)
    1401. 

  1401.       …
    1402.

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