PageRenderTime 72ms CodeModel.GetById 19ms RepoModel.GetById 0ms app.codeStats 1ms

/wrfv2_fire/phys/module_mp_gsfcgce.F

http://github.com/jbeezley/wrf-fire
FORTRAN Legacy | 2838 lines | 1694 code | 322 blank | 822 comment | 147 complexity | 4686bc158b95f5e42c5d6615ed17ee47 MD5 | raw file
Possible License(s): AGPL-1.0 
    

  1. !WRF:MODEL_LAYER:PHYSICS
    2. 

  2. !
    3. 

  3. 

  4. MODULE module_mp_gsfcgce
    5. 

  5. 

  6.    USE     module_wrf_error
    7. 

  7. 

  8. !JJS 1/3/2008     vvvvv
    9. 

  9. 

  10. !  common /bt/
    11. 

  11.    REAL,    PRIVATE ::          rd1,  rd2,   al,   cp
    12. 

  12. 

  13. !  common /cont/
    14. 

  14.    REAL,    PRIVATE ::          c38, c358, c610, c149, &
    15. 

  15.                                c879, c172, c409,  c76, &
    16. 

  16.                                c218, c580, c141
    17. 

  17. !  common /b3cs/
    18. 

  18.    REAL,    PRIVATE ::           ag,   bg,   as,   bs, &
    19. 

  19.                                  aw,   bw,  bgh,  bgq, &
    20. 

  20.                                 bsh,  bsq,  bwh,  bwq
    21. 

  21. 

  22. !  common /size/
    23. 

  23.    REAL,    PRIVATE ::          tnw,  tns,  tng,       &
    24. 

  24.                                roqs, roqg, roqr
    25. 

  25. 

  26. !  common /bterv/
    27. 

  27.    REAL,    PRIVATE ::          zrc,  zgc,  zsc,       &
    28. 

  28.                                 vrc,  vgc,  vsc
    29. 

  29. 

  30. !  common /bsnw/
    31. 

  31.    REAL,    PRIVATE ::          alv,  alf,  als,   t0,   t00,     &
    32. 

  32.                                 avc,  afc,  asc,  rn1,  bnd1,     &
    33. 

  33.                                 rn2, bnd2,  rn3,  rn4,   rn5,     &
    34. 

  34.                                 rn6,  rn7,  rn8,  rn9,  rn10,     &
    35. 

  35.                               rn101,rn10a, rn11,rn11a,  rn12
    36. 

  36. 

  37.    REAL,    PRIVATE ::         rn14, rn15,rn15a, rn16,  rn17,     &
    38. 

  38.                               rn17a,rn17b,rn17c, rn18, rn18a,     &
    39. 

  39.                                rn19,rn19a,rn19b, rn20, rn20a,     &
    40. 

  40.                               rn20b, bnd3, rn21, rn22,  rn23,     &
    41. 

  41.                               rn23a,rn23b, rn25,rn30a, rn30b,     &
    42. 

  42.                               rn30c, rn31, beta, rn32
    43. 

  43. 

  44.    REAL,    PRIVATE, DIMENSION( 31 ) ::    rn12a, rn12b, rn13, rn25a
    45. 

  45. 

  46. !  common /rsnw1/
    47. 

  47.    REAL,    PRIVATE ::         rn10b, rn10c, rnn191, rnn192,  rn30,     &
    48. 

  48.                              rnn30a,  rn33,  rn331,  rn332
    49. 

  49. 

  50. !
    51. 

  51.    REAL,    PRIVATE, DIMENSION( 31 )  ::      aa1,  aa2
    52. 

  52.    DATA aa1/.7939e-7, .7841e-6, .3369e-5, .4336e-5, .5285e-5,     &
    53. 

  53.            .3728e-5, .1852e-5, .2991e-6, .4248e-6, .7434e-6,     &
    54. 

  54.            .1812e-5, .4394e-5, .9145e-5, .1725e-4, .3348e-4,     &
    55. 

  55.            .1725e-4, .9175e-5, .4412e-5, .2252e-5, .9115e-6,     &
    56. 

  56.            .4876e-6, .3473e-6, .4758e-6, .6306e-6, .8573e-6,     &
    57. 

  57.            .7868e-6, .7192e-6, .6513e-6, .5956e-6, .5333e-6,     &
    58. 

  58.            .4834e-6/
    59. 

  59.    DATA aa2/.4006, .4831, .5320, .5307, .5319,      &
    60. 

  60.            .5249, .4888, .3894, .4047, .4318,      &
    61. 

  61.            .4771, .5183, .5463, .5651, .5813,      &
    62. 

  62.            .5655, .5478, .5203, .4906, .4447,      &
    63. 

  63.            .4126, .3960, .4149, .4320, .4506,      &
    64. 

  64.            .4483, .4460, .4433, .4413, .4382,      &
    65. 

  65.            .4361/
    66. 

  66. 

  67. !JJS 1/3/2008     ^^^^^
    68. 

  68. 

  69. CONTAINS
    70. 

  70. 

  71. !-------------------------------------------------------------------
    72. 

  72. !  NASA/GSFC GCE
    73. 

  73. !  Tao et al, 2001, Meteo. & Atmos. Phy., 97-137
    74. 

  74. !-------------------------------------------------------------------
    75. 

  75. !  SUBROUTINE gsfcgce(  th, th_old,                                 &
    76. 

  76.   SUBROUTINE gsfcgce(  th,                                         &
    77. 

  77.                        qv, ql,                                     &
    78. 

  78.                        qr, qi,                                     &
    79. 

  79.                        qs,                                         &
    80. 

  80. !                       qvold, qlold,                               &
    81. 

  81. !                       qrold, qiold,                               &
    82. 

  82. !                       qsold,                                      &
    83. 

  83.                        rho, pii, p, dt_in, z,                      &
    84. 

  84.                        ht, dz8w, grav,                             &
    85. 

  85.                        rhowater, rhosnow,                          &
    86. 

  86.                        itimestep,                                  &
    87. 

  87.                        ids,ide, jds,jde, kds,kde,                  & ! domain dims
    88. 

  88.                        ims,ime, jms,jme, kms,kme,                  & ! memory dims
    89. 

  89.                        its,ite, jts,jte, kts,kte,                  & ! tile   dims
    90. 

  90.                        rainnc, rainncv,                            &
    91. 

  91.                        snownc, snowncv, sr,                        &
    92. 

  92.                        graupelnc, graupelncv,                      &
    93. 

  93. !                       f_qg, qg, pgold,                            &
    94. 

  94.                        f_qg, qg,                                   &
    95. 

  95.                        ihail, ice2                                 &
    96. 

  96.                                                                    )
    97. 

  97. 

  98. !-------------------------------------------------------------------
    99. 

  99.   IMPLICIT NONE
    100. 

  100. !-------------------------------------------------------------------
    101. 

  101. !
    102. 

  102. ! JJS 2/15/2005
    103. 

  103. !
    104. 

  104.   INTEGER,      INTENT(IN   )    ::   ids,ide, jds,jde, kds,kde , &
    105. 

  105.                                       ims,ime, jms,jme, kms,kme , &
    106. 

  106.                                       its,ite, jts,jte, kts,kte 
    107. 

  107.   INTEGER,      INTENT(IN   )    ::   itimestep, ihail, ice2 
    108. 

  108. 

  109.   REAL, DIMENSION( ims:ime , kms:kme , jms:jme ),                 &
    110. 

  110.         INTENT(INOUT) ::                                          &
    111. 

  111.                                                               th, &
    112. 

  112.                                                               qv, &
    113. 

  113.                                                               ql, &
    114. 

  114.                                                               qr, &
    115. 

  115.                                                               qi, &
    116. 

  116.                                                               qs, &
    117. 

  117.                                                               qg
    118. 

  118. !
    119. 

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

  120.         INTENT(IN   ) ::                                          &
    121. 

  121. !                                                         th_old, &
    122. 

  122. !                                                          qvold, &
    123. 

  123. !                                                          qlold, &
    124. 

  124. !                                                          qrold, &
    125. 

  125. !                                                          qiold, &
    126. 

  126. !                                                          qsold, &
    127. 

  127. !                                                          qgold, &
    128. 

  128.                                                              rho, &
    129. 

  129.                                                              pii, &
    130. 

  130.                                                                p, &
    131. 

  131.                                                             dz8w, &
    132. 

  132.                                                                z
    133. 

  133. 

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

  135.         INTENT(INOUT) ::                               rainnc,    &
    136. 

  136.                                                        rainncv,   &
    137. 

  137.                                                        snownc,    &   
    138. 

  138.                                                        snowncv,   &
    139. 

  139.                                                        sr,        &
    140. 

  140.                                                        graupelnc, &
    141. 

  141.                                                        graupelncv 
    142. 

  142. 

  143.   REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN) ::       ht
    144. 

  144. 

  145.   REAL, INTENT(IN   ) ::                                   dt_in, &
    146. 

  146.                                                             grav, &
    147. 

  147.                                                         rhowater, &
    148. 

  148.                                                          rhosnow
    149. 

  149. 

  150.   LOGICAL, INTENT(IN), OPTIONAL :: F_QG
    151. 

  151. 

  152. !  LOCAL VAR
    153. 

  153. 

  154. 

  155. !jjs  INTEGER             ::                            min_q, max_q
    156. 

  156. 

  157. !jjs  REAL, DIMENSION( its:ite , jts:jte )                            &
    158. 

  158. !jjs                               ::        rain, snow, graupel,ice
    159. 

  159. 

  160. !
    161. 

  161. !  INTEGER :: IHAIL, itaobraun, ice2, istatmin, new_ice_sat, id
    162. 

  162.   INTEGER ::  itaobraun, istatmin, new_ice_sat, id
    163. 

  163. 

  164.   INTEGER :: i, j, k
    165. 

  165.   INTEGER :: iskip, ih, icount, ibud, i24h 
    166. 

  166.   REAL    :: hour
    167. 

  167.   REAL    , PARAMETER :: cmin=1.e-20
    168. 

  168.   REAL    :: dth, dqv, dqrest, dqall, dqall1, rhotot, a1, a2 
    169. 

  169. !  REAL, DIMENSION( its:ite , kts:kte , jts:jte ) ::                   &
    170. 

  170. !                         th1, qv1, ql1, qr1, qi1, qs1, qg1
    171. 

  171.  
    172. 

  172.   LOGICAL :: flag_qg
    173. 

  173. 

  174. !
    175. 

  175. !c  ihail = 0    for graupel, for tropical region
    176. 

  176. !c  ihail = 1    for hail, for mid-lat region
    177. 

  177. 

  178. ! itaobraun: 0 for Tao's constantis, 1 for Braun's constants
    179. 

  179. !c        if ( itaobraun.eq.1 ) --> betah=0.5*beta=-.46*0.5=-0.23;   cn0=1.e-6
    180. 

  180. !c        if ( itaobraun.eq.0 ) --> betah=0.5*beta=-.6*0.5=-0.30;    cn0=1.e-8
    181. 

  181.    itaobraun = 1
    182. 

  182. 

  183. !c  ice2 = 0    for 3 ice --- ice, snow and graupel/hail
    184. 

  184. !c  ice2 = 1    for 2 ice --- ice and snow only
    185. 

  185. !c  ice2 = 2    for 2 ice --- ice and graupel only, use ihail = 0 only
    186. 

  186. !c  ice2 = 3    for 0 ice --- no ice, warm only
    187. 

  187. 

  188. !  if (ice2 .eq. 2) ihail = 0
    189. 

  189. 

  190.   i24h=nint(86400./dt_in)
    191. 

  191.   if (mod(itimestep,i24h).eq.1) then
    192. 

  192.      write(6,*) 'ihail=',ihail,'  ice2=',ice2
    193. 

  193.      if (ice2.eq.0) then
    194. 

  194.         write(6,*) 'Running 3-ice scheme in GSFCGCE with'
    195. 

  195.         if (ihail.eq.0) then 
    196. 

  196.            write(6,*) '     ice, snow and graupel'
    197. 

  197.         else if (ihail.eq.1) then
    198. 

  198.                 write(6,*) '     ice, snow and hail'
    199. 

  199.         else
    200. 

  200.              write(6,*) 'ihail has to be either 1 or 0'
    201. 

  201.              stop
    202. 

  202.         endif !ihail
    203. 

  203.      else if (ice2.eq.1) then
    204. 

  204.              write(6,*) 'Running 2-ice scheme in GSFCGCE with'
    205. 

  205.              write(6,*) '     ice and snow'
    206. 

  206.      else if (ice2.eq.2) then
    207. 

  207.              write(6,*) 'Running 2-ice scheme in GSFCGCE with'
    208. 

  208.              write(6,*) '     ice and graupel'
    209. 

  209.      else if (ice2.eq.3) then
    210. 

  210.              write(6,*) 'Running warm rain only scheme in GSFCGCE without any ice'
    211. 

  211.      else
    212. 

  212.              write(6,*) 'gsfcgce_2ice in namelist.input has to be 0, 1, 2, or 3'
    213. 

  213.              stop
    214. 

  214.      endif !ice2
    215. 

  215.   endif !itimestep
    216. 

  216. 

  217. !c  new_ice_sat = 0, 1 or 2 
    218. 

  218.     new_ice_sat = 2 
    219. 

  219. 

  220. !c istatmin
    221. 

  221.     istatmin = 180
    222. 

  222. 

  223. !c id = 0  without in-line staticstics
    224. 

  224. !c id = 1  with in-line staticstics
    225. 

  225.     id = 0
    226. 

  226. 

  227. !c ibud = 0 no calculation of dth, dqv, dqrest and dqall
    228. 

  228. !c ibud = 1 yes
    229. 

  229.     ibud = 0
    230. 

  230. 

  231. !jjs   dt=dt_in
    232. 

  232. !jjs   rhoe_s=1.29
    233. 

  233. !
    234. 

  234. !   IF (P_QI .lt. P_FIRST_SCALAR .or. P_QS .lt. P_FIRST_SCALAR) THEN
    235. 

  235. !      CALL wrf_error_fatal3 ( "module_mp_lin.b" , 130 ,  'module_mp_lin: Improper use of Lin et al scheme; no ice phase. Please chose another one.')
    236. 

  236. !   ENDIF
    237. 

  237. 

  238. ! calculte fallflux and precipiation in MKS system
    239. 

  239. 

  240.    call fall_flux(dt_in, qr, qi, qs, qg, p,                   &
    241. 

  241.                       rho, z, dz8w, ht, rainnc,               &
    242. 

  242.                       rainncv, grav,itimestep,                &
    243. 

  243.                       rhowater, rhosnow,                      &
    244. 

  244.                       snownc, snowncv, sr,                    &
    245. 

  245.                       graupelnc, graupelncv,                  &
    246. 

  246.                       ihail, ice2,                            &
    247. 

  247.                       ims,ime, jms,jme, kms,kme,              & ! memory dims
    248. 

  248.                       its,ite, jts,jte, kts,kte               ) ! tile   dims
    249. 

  249. !-----------------------------------------------------------------------
    250. 

  250. 

  251. !c  set up constants used internally in GCE
    252. 

  252. 

  253.    call consat_s (ihail, itaobraun)
    254. 

  254. 

  255. 

  256. !c Negative values correction
    257. 

  257. 

  258.    iskip = 1
    259. 

  259.  
    260. 

  260.    if (iskip.eq.0) then
    261. 

  261.       call negcor(qv,rho,dz8w,ims,ime,jms,jme,kms,kme, &
    262. 

  262.                            itimestep,1,             &
    263. 

  263.                            its,ite,jts,jte,kts,kte)
    264. 

  264.       call negcor(ql,rho,dz8w,ims,ime,jms,jme,kms,kme, &
    265. 

  265.                            itimestep,2,             &
    266. 

  266.                            its,ite,jts,jte,kts,kte)
    267. 

  267.       call negcor(qr,rho,dz8w,ims,ime,jms,jme,kms,kme, &
    268. 

  268.                            itimestep,3,             &
    269. 

  269.                            its,ite,jts,jte,kts,kte)
    270. 

  270.       call negcor(qi,rho,dz8w,ims,ime,jms,jme,kms,kme, &
    271. 

  271.                            itimestep,4,             &
    272. 

  272.                            its,ite,jts,jte,kts,kte)
    273. 

  273.       call negcor(qs,rho,dz8w,ims,ime,jms,jme,kms,kme, &
    274. 

  274.                            itimestep,5,             &
    275. 

  275.                            its,ite,jts,jte,kts,kte)
    276. 

  276.       call negcor(qg,rho,dz8w,ims,ime,jms,jme,kms,kme, &
    277. 

  277.                            itimestep,6,             &
    278. 

  278.                            its,ite,jts,jte,kts,kte)
    279. 

  279. !   else if (mod(itimestep,i24h).eq.1) then
    280. 

  280. !      print *,'no neg correction in mp at timestep=',itimestep
    281. 

  281.    endif ! iskip
    282. 

  282. 

  283. !c microphysics in GCE
    284. 

  284. 

  285.    call SATICEL_S( dt_in, IHAIL, itaobraun, ICE2, istatmin,     &
    286. 

  286.                    new_ice_sat, id,                             &
    287. 

  287. !                   th, th_old, qv, ql, qr,                      &
    288. 

  288.                    th, qv, ql, qr,                      &
    289. 

  289.                    qi, qs, qg,                                  &
    290. 

  290. !                   qvold, qlold, qrold,                         &
    291. 

  291. !                   qiold, qsold, qgold,                         &
    292. 

  292.                    rho, pii, p, itimestep,                      & 
    293. 

  293.                    ids,ide, jds,jde, kds,kde,                   & ! domain dims
    294. 

  294.                    ims,ime, jms,jme, kms,kme,                   & ! memory dims
    295. 

  295.                    its,ite, jts,jte, kts,kte                    & ! tile   dims
    296. 

  296.                                                                 ) 
    297. 

  297. 

  298.    END SUBROUTINE gsfcgce
    299. 

  299. 

  300. !-----------------------------------------------------------------------
    301. 

  301.    SUBROUTINE fall_flux ( dt, qr, qi, qs, qg, p,              &
    302. 

  302.                       rho, z, dz8w, topo, rainnc,             &
    303. 

  303.                       rainncv, grav, itimestep,               &
    304. 

  304.                       rhowater, rhosnow,                      &
    305. 

  305.                       snownc, snowncv, sr,                    &
    306. 

  306.                       graupelnc, graupelncv,                  &
    307. 

  307.                       ihail, ice2,                            &
    308. 

  308.                       ims,ime, jms,jme, kms,kme,              & ! memory dims
    309. 

  309.                       its,ite, jts,jte, kts,kte               ) ! tile   dims
    310. 

  310. !-----------------------------------------------------------------------
    311. 

  311. ! adopted from Jiun-Dar Chern's codes for Purdue Regional Model
    312. 

  312. ! adopted by Jainn J. Shi, 6/10/2005
    313. 

  313. !-----------------------------------------------------------------------
    314. 

  314. 

  315.   IMPLICIT NONE
    316. 

  316.   INTEGER, INTENT(IN   )               :: ihail, ice2,                &
    317. 

  317.                                           ims,ime, jms,jme, kms,kme,  &
    318. 

  318.                                           its,ite, jts,jte, kts,kte 
    319. 

  319.   INTEGER, INTENT(IN   )               :: itimestep
    320. 

  320.   REAL,    DIMENSION( ims:ime , kms:kme , jms:jme ),                  &
    321. 

  321.            INTENT(INOUT)               :: qr, qi, qs, qg       
    322. 

  322.   REAL,    DIMENSION( ims:ime , jms:jme ),                            &
    323. 

  323.            INTENT(INOUT)               :: rainnc, rainncv,            &
    324. 

  324.                                           snownc, snowncv, sr,        &
    325. 

  325.                                           graupelnc, graupelncv
    326. 

  326.   REAL,    DIMENSION( ims:ime , kms:kme , jms:jme ),                  &
    327. 

  327.            INTENT(IN   )               :: rho, z, dz8w, p     
    328. 

  328. 

  329.   REAL,    INTENT(IN   )               :: dt, grav, rhowater, rhosnow
    330. 

  330. 

  331. 

  332.   REAL,    DIMENSION( ims:ime , jms:jme ),                            &
    333. 

  333.            INTENT(IN   )               :: topo   
    334. 

  334. 

  335. ! temperary vars
    336. 

  336. 

  337.   REAL,    DIMENSION( kts:kte )           :: sqrhoz
    338. 

  338.   REAL                                    :: tmp1, term0
    339. 

  339.   REAL                                :: pptrain, pptsnow,        &
    340. 

  340.                                          pptgraul, pptice
    341. 

  341.   REAL,    DIMENSION( kts:kte )       :: qrz, qiz, qsz, qgz,      &
    342. 

  342.                                          zz, dzw, prez, rhoz,     &
    343. 

  343.                                          orhoz
    344. 

  344. 

  345. 

  346.    INTEGER                    :: k, i, j
    347. 

  347. !
    348. 

  348. 

  349.   REAL, DIMENSION( kts:kte )    :: vtr, vts, vtg, vti
    350. 

  350. 

  351.   REAL                          :: dtb, pi, consta, constc, gambp4,    &
    352. 

  352.                                    gamdp4, gam4pt5, gam4bbar
    353. 

  353. 

  354. !  Lin
    355. 

  355.    REAL    , PARAMETER ::     xnor = 8.0e6
    356. 

  356. !   REAL    , PARAMETER ::     xnos = 3.0e6
    357. 

  357.    REAL    , PARAMETER ::     xnos = 1.6e7   ! Tao's value
    358. 

  358.    REAL    , PARAMETER ::                              &
    359. 

  359. !             constb = 0.8, constd = 0.25, o6 = 1./6.,           &
    360. 

  360.              constb = 0.8, constd = 0.11, o6 = 1./6.,           &
    361. 

  361.              cdrag = 0.6
    362. 

  362. !  Lin
    363. 

  363. !  REAL    , PARAMETER ::     xnoh = 4.0e4
    364. 

  364.   REAL    , PARAMETER ::     xnoh = 2.0e5    ! Tao's value
    365. 

  365.   REAL    , PARAMETER ::     rhohail = 917.
    366. 

  366. 

  367. ! Hobbs
    368. 

  368.   REAL    , PARAMETER ::     xnog = 4.0e6
    369. 

  369.   REAL    , PARAMETER ::     rhograul = 400.
    370. 

  370.   REAL    , PARAMETER ::     abar = 19.3, bbar = 0.37,      &
    371. 

  371.                                       p0 = 1.0e5
    372. 

  372. 

  373.   REAL    , PARAMETER ::     rhoe_s = 1.29
    374. 

  374. 

  375. ! for terminal velocity flux
    376. 

  376.   INTEGER                       :: min_q, max_q
    377. 

  377.   REAL                          :: t_del_tv, del_tv, flux, fluxin, fluxout ,tmpqrz
    378. 

  378.   LOGICAL                       :: notlast
    379. 

  379. 

  380. !  if (itimestep.eq.1) then
    381. 

  381. !     write(6, *) 'in fall_flux'
    382. 

  382. !     write(6, *) 'ims=', ims, '  ime=', ime
    383. 

  383. !     write(6, *) 'jms=', jms, '  jme=', jme
    384. 

  384. !     write(6, *) 'kms=', kms, '  kme=', kme
    385. 

  385. !     write(6, *) 'its=', its, '  ite=', ite
    386. 

  386. !     write(6, *) 'jts=', jts, '  jte=', jte
    387. 

  387. !     write(6, *) 'kts=', kts, '  kte=', kte
    388. 

  388. !     write(6, *) 'dt=', dt
    389. 

  389. !     write(6, *) 'ihail=', ihail
    390. 

  390. !     write(6, *) 'ICE2=', ICE2
    391. 

  391. !     write(6, *) 'dt=', dt
    392. 

  392. !   endif 
    393. 

  393. 

  394. !-----------------------------------------------------------------------
    395. 

  395. !  This program calculates precipitation fluxes due to terminal velocities.
    396. 

  396. !-----------------------------------------------------------------------
    397. 

  397. 

  398.    dtb=dt
    399. 

  399.    pi=acos(-1.)
    400. 

  400.    consta=2115.0*0.01**(1-constb)
    401. 

  401. !   constc=152.93*0.01**(1-constd)
    402. 

  402.    constc=78.63*0.01**(1-constd)
    403. 

  403. 

  404. !  Gamma function
    405. 

  405.    gambp4=ggamma(constb+4.)
    406. 

  406.    gamdp4=ggamma(constd+4.)
    407. 

  407.    gam4pt5=ggamma(4.5)
    408. 

  408.    gam4bbar=ggamma(4.+bbar)
    409. 

  409. 

  410. !***********************************************************************
    411. 

  411. ! Calculate precipitation fluxes due to terminal velocities.
    412. 

  412. !***********************************************************************
    413. 

  413. !
    414. 

  414. !- Calculate termianl velocity (vt?)  of precipitation q?z
    415. 

  415. !- Find maximum vt? to determine the small delta t
    416. 

  416. 

  417.  j_loop:  do j = jts, jte
    418. 

  418.  i_loop:  do i = its, ite
    419. 

  419. 

  420.    pptrain = 0.
    421. 

  421.    pptsnow = 0.
    422. 

  422.    pptgraul = 0.
    423. 

  423.    pptice  = 0.
    424. 

  424. 

  425.    do k = kts, kte
    426. 

  426.       qrz(k)=qr(i,k,j)
    427. 

  427.       rhoz(k)=rho(i,k,j)
    428. 

  428.       orhoz(k)=1./rhoz(k)
    429. 

  429.       prez(k)=p(i,k,j)
    430. 

  430.       sqrhoz(k)=sqrt(rhoe_s/rhoz(k))
    431. 

  431.       zz(k)=z(i,k,j)
    432. 

  432.       dzw(k)=dz8w(i,k,j)
    433. 

  433.    enddo !k
    434. 

  434. 

  435.       DO k = kts, kte
    436. 

  436.          qiz(k)=qi(i,k,j)
    437. 

  437.       ENDDO
    438. 

  438. 

  439.       DO k = kts, kte
    440. 

  440.          qsz(k)=qs(i,k,j)
    441. 

  441.       ENDDO
    442. 

  442. 

  443.    IF (ice2 .eq. 0) THEN
    444. 

  444.       DO k = kts, kte
    445. 

  445.          qgz(k)=qg(i,k,j)
    446. 

  446.       ENDDO
    447. 

  447.    ELSE
    448. 

  448.       DO k = kts, kte
    449. 

  449.          qgz(k)=0.
    450. 

  450.       ENDDO
    451. 

  451.    ENDIF
    452. 

  452. 

  453. 

  454. !
    455. 

  455. !-- rain
    456. 

  456. !
    457. 

  457.     t_del_tv=0.
    458. 

  458.     del_tv=dtb
    459. 

  459.     notlast=.true.
    460. 

  460.     DO while (notlast)
    461. 

  461. !
    462. 

  462.       min_q=kte
    463. 

  463.       max_q=kts-1
    464. 

  464. !
    465. 

  465.       do k=kts,kte-1
    466. 

  466.          if (qrz(k) .gt. 1.0e-8) then
    467. 

  467.             min_q=min0(min_q,k)
    468. 

  468.             max_q=max0(max_q,k)
    469. 

  469.             tmp1=sqrt(pi*rhowater*xnor/rhoz(k)/qrz(k))
    470. 

  470.             tmp1=sqrt(tmp1)
    471. 

  471.             vtr(k)=consta*gambp4*sqrhoz(k)/tmp1**constb
    472. 

  472.             vtr(k)=vtr(k)/6.
    473. 

  473.             if (k .eq. 1) then
    474. 

  474.                del_tv=amin1(del_tv,0.9*(zz(k)-topo(i,j))/vtr(k))
    475. 

  475.             else
    476. 

  476.                del_tv=amin1(del_tv,0.9*(zz(k)-zz(k-1))/vtr(k))
    477. 

  477.             endif
    478. 

  478.          else
    479. 

  479.             vtr(k)=0.
    480. 

  480.          endif
    481. 

  481.       enddo
    482. 

  482. 

  483.       if (max_q .ge. min_q) then
    484. 

  484. !
    485. 

  485. !- Check if the summation of the small delta t >=  big delta t
    486. 

  486. !             (t_del_tv)          (del_tv)             (dtb)
    487. 

  487. 

  488.          t_del_tv=t_del_tv+del_tv
    489. 

  489. !
    490. 

  490.          if ( t_del_tv .ge. dtb ) then
    491. 

  491.               notlast=.false.
    492. 

  492.               del_tv=dtb+del_tv-t_del_tv
    493. 

  493.          endif
    494. 

  494. 

  495. ! use small delta t to calculate the qrz flux
    496. 

  496. ! termi is the qrz flux pass in the grid box through the upper boundary
    497. 

  497. ! termo is the qrz flux pass out the grid box through the lower boundary
    498. 

  498. !
    499. 

  499.          fluxin=0.
    500. 

  500.          do k=max_q,min_q,-1
    501. 

  501.             fluxout=rhoz(k)*vtr(k)*qrz(k)
    502. 

  502.             flux=(fluxin-fluxout)/rhoz(k)/dzw(k)
    503. 

  503. !            tmpqrz=qrz(k)
    504. 

  504.             qrz(k)=qrz(k)+del_tv*flux
    505. 

  505.             qrz(k)=amax1(0.,qrz(k))
    506. 

  506.             qr(i,k,j)=qrz(k)
    507. 

  507.             fluxin=fluxout
    508. 

  508.          enddo
    509. 

  509.          if (min_q .eq. 1) then
    510. 

  510.             pptrain=pptrain+fluxin*del_tv
    511. 

  511.          else
    512. 

  512.             qrz(min_q-1)=qrz(min_q-1)+del_tv*  &
    513. 

  513.                           fluxin/rhoz(min_q-1)/dzw(min_q-1)
    514. 

  514.             qr(i,min_q-1,j)=qrz(min_q-1)
    515. 

  515.          endif
    516. 

  516. !
    517. 

  517.       else
    518. 

  518.          notlast=.false.
    519. 

  519.       endif
    520. 

  520.     ENDDO
    521. 

  521. 

  522. !
    523. 

  523. !-- snow
    524. 

  524. !
    525. 

  525.     t_del_tv=0.
    526. 

  526.     del_tv=dtb
    527. 

  527.     notlast=.true.
    528. 

  528. 

  529.     DO while (notlast)
    530. 

  530. !
    531. 

  531.       min_q=kte
    532. 

  532.       max_q=kts-1
    533. 

  533. !
    534. 

  534.       do k=kts,kte-1
    535. 

  535.          if (qsz(k) .gt. 1.0e-8) then
    536. 

  536.             min_q=min0(min_q,k)
    537. 

  537.             max_q=max0(max_q,k)
    538. 

  538.             tmp1=sqrt(pi*rhosnow*xnos/rhoz(k)/qsz(k))
    539. 

  539.             tmp1=sqrt(tmp1)
    540. 

  540.             vts(k)=constc*gamdp4*sqrhoz(k)/tmp1**constd
    541. 

  541.             vts(k)=vts(k)/6.
    542. 

  542.             if (k .eq. 1) then
    543. 

  543.                del_tv=amin1(del_tv,0.9*(zz(k)-topo(i,j))/vts(k))
    544. 

  544.             else
    545. 

  545.                del_tv=amin1(del_tv,0.9*(zz(k)-zz(k-1))/vts(k))
    546. 

  546.             endif
    547. 

  547.          else
    548. 

  548.             vts(k)=0.
    549. 

  549.          endif
    550. 

  550.       enddo
    551. 

  551. 

  552.       if (max_q .ge. min_q) then
    553. 

  553. !
    554. 

  554. !
    555. 

  555. !- Check if the summation of the small delta t >=  big delta t
    556. 

  556. !             (t_del_tv)          (del_tv)             (dtb)
    557. 

  557. 

  558.          t_del_tv=t_del_tv+del_tv
    559. 

  559. 

  560.          if ( t_del_tv .ge. dtb ) then
    561. 

  561.               notlast=.false.
    562. 

  562.               del_tv=dtb+del_tv-t_del_tv
    563. 

  563.          endif
    564. 

  564. 

  565. ! use small delta t to calculate the qsz flux
    566. 

  566. ! termi is the qsz flux pass in the grid box through the upper boundary
    567. 

  567. ! termo is the qsz flux pass out the grid box through the lower boundary
    568. 

  568. !
    569. 

  569.          fluxin=0.
    570. 

  570.          do k=max_q,min_q,-1
    571. 

  571.             fluxout=rhoz(k)*vts(k)*qsz(k)
    572. 

  572.             flux=(fluxin-fluxout)/rhoz(k)/dzw(k)
    573. 

  573.             qsz(k)=qsz(k)+del_tv*flux
    574. 

  574.             qsz(k)=amax1(0.,qsz(k))
    575. 

  575.             qs(i,k,j)=qsz(k)
    576. 

  576.             fluxin=fluxout
    577. 

  577.          enddo
    578. 

  578.          if (min_q .eq. 1) then
    579. 

  579.             pptsnow=pptsnow+fluxin*del_tv
    580. 

  580.          else
    581. 

  581.             qsz(min_q-1)=qsz(min_q-1)+del_tv*  &
    582. 

  582.                          fluxin/rhoz(min_q-1)/dzw(min_q-1)
    583. 

  583.             qs(i,min_q-1,j)=qsz(min_q-1)
    584. 

  584.          endif
    585. 

  585. !
    586. 

  586.       else
    587. 

  587.          notlast=.false.
    588. 

  588.       endif
    589. 

  589. 

  590.     ENDDO
    591. 

  591. 

  592. !
    593. 

  593. !   ice2=0 --- with hail/graupel 
    594. 

  594. !   ice2=1 --- without hail/graupel 
    595. 

  595. !
    596. 

  596.   if (ice2.eq.0) then 
    597. 

  597. !
    598. 

  598. !-- If IHAIL=1, use hail.
    599. 

  599. !-- If IHAIL=0, use graupel.
    600. 

  600. !
    601. 

  601. !    if (ihail .eq. 1) then
    602. 

  602. !       xnog = xnoh
    603. 

  603. !       rhograul = rhohail
    604. 

  604. !    endif
    605. 

  605. 

  606.     t_del_tv=0.
    607. 

  607.     del_tv=dtb
    608. 

  608.     notlast=.true.
    609. 

  609. !
    610. 

  610.     DO while (notlast)
    611. 

  611. !
    612. 

  612.       min_q=kte
    613. 

  613.       max_q=kts-1
    614. 

  614. !
    615. 

  615.       do k=kts,kte-1
    616. 

  616.          if (qgz(k) .gt. 1.0e-8) then
    617. 

  617.             if (ihail .eq. 1) then
    618. 

  618. !  for hail, based on Lin et al (1983)
    619. 

  619.               min_q=min0(min_q,k)
    620. 

  620.               max_q=max0(max_q,k)
    621. 

  621.               tmp1=sqrt(pi*rhohail*xnoh/rhoz(k)/qgz(k))
    622. 

  622.               tmp1=sqrt(tmp1)
    623. 

  623.               term0=sqrt(4.*grav*rhohail/3./rhoz(k)/cdrag)
    624. 

  624.               vtg(k)=gam4pt5*term0*sqrt(1./tmp1)
    625. 

  625.               vtg(k)=vtg(k)/6.
    626. 

  626.               if (k .eq. 1) then
    627. 

  627.                  del_tv=amin1(del_tv,0.9*(zz(k)-topo(i,j))/vtg(k))
    628. 

  628.               else
    629. 

  629.                  del_tv=amin1(del_tv,0.9*(zz(k)-zz(k-1))/vtg(k))
    630. 

  630.               endif !k
    631. 

  631.             else
    632. 

  632. ! added by JJS
    633. 

  633. ! for graupel, based on RH (1984)
    634. 

  634.               min_q=min0(min_q,k)
    635. 

  635.               max_q=max0(max_q,k)
    636. 

  636.               tmp1=sqrt(pi*rhograul*xnog/rhoz(k)/qgz(k))
    637. 

  637.               tmp1=sqrt(tmp1)
    638. 

  638.               tmp1=tmp1**bbar
    639. 

  639.               tmp1=1./tmp1
    640. 

  640.               term0=abar*gam4bbar/6.
    641. 

  641.               vtg(k)=term0*tmp1*(p0/prez(k))**0.4
    642. 

  642.               if (k .eq. 1) then
    643. 

  643.                  del_tv=amin1(del_tv,0.9*(zz(k)-topo(i,j))/vtg(k))
    644. 

  644.               else
    645. 

  645.                  del_tv=amin1(del_tv,0.9*(zz(k)-zz(k-1))/vtg(k))
    646. 

  646.               endif !k
    647. 

  647.             endif !ihail
    648. 

  648.          else
    649. 

  649.             vtg(k)=0.
    650. 

  650.          endif !qgz
    651. 

  651.       enddo !k
    652. 

  652. 

  653.       if (max_q .ge. min_q) then
    654. 

  654. !
    655. 

  655. !
    656. 

  656. !- Check if the summation of the small delta t >=  big delta t
    657. 

  657. !             (t_del_tv)          (del_tv)             (dtb)
    658. 

  658. 

  659.          t_del_tv=t_del_tv+del_tv
    660. 

  660. 

  661.          if ( t_del_tv .ge. dtb ) then
    662. 

  662.               notlast=.false.
    663. 

  663.               del_tv=dtb+del_tv-t_del_tv
    664. 

  664.          endif
    665. 

  665. 

  666. ! use small delta t to calculate the qgz flux
    667. 

  667. ! termi is the qgz flux pass in the grid box through the upper boundary
    668. 

  668. ! termo is the qgz flux pass out the grid box through the lower boundary
    669. 

  669. !
    670. 

  670.          fluxin=0.
    671. 

  671.          do k=max_q,min_q,-1
    672. 

  672.             fluxout=rhoz(k)*vtg(k)*qgz(k)
    673. 

  673.             flux=(fluxin-fluxout)/rhoz(k)/dzw(k)
    674. 

  674.             qgz(k)=qgz(k)+del_tv*flux
    675. 

  675.             qgz(k)=amax1(0.,qgz(k))
    676. 

  676.             qg(i,k,j)=qgz(k)
    677. 

  677.             fluxin=fluxout
    678. 

  678.          enddo
    679. 

  679.          if (min_q .eq. 1) then
    680. 

  680.             pptgraul=pptgraul+fluxin*del_tv
    681. 

  681.          else
    682. 

  682.             qgz(min_q-1)=qgz(min_q-1)+del_tv*  &
    683. 

  683.                          fluxin/rhoz(min_q-1)/dzw(min_q-1)
    684. 

  684.             qg(i,min_q-1,j)=qgz(min_q-1)
    685. 

  685.          endif
    686. 

  686. !
    687. 

  687.       else
    688. 

  688.          notlast=.false.
    689. 

  689.       endif
    690. 

  690. !
    691. 

  691.     ENDDO
    692. 

  692.  ENDIF !ice2
    693. 

  693. !
    694. 

  694. !-- cloud ice  (03/21/02) follow Vaughan T.J. Phillips at GFDL
    695. 

  695. !
    696. 

  696. 

  697.     t_del_tv=0.
    698. 

  698.     del_tv=dtb
    699. 

  699.     notlast=.true.
    700. 

  700. !
    701. 

  701.     DO while (notlast)
    702. 

  702. !
    703. 

  703.       min_q=kte
    704. 

  704.       max_q=kts-1
    705. 

  705. !
    706. 

  706.       do k=kts,kte-1
    707. 

  707.          if (qiz(k) .gt. 1.0e-8) then
    708. 

  708.             min_q=min0(min_q,k)
    709. 

  709.             max_q=max0(max_q,k)
    710. 

  710.             vti(k)= 3.29 * (rhoz(k)* qiz(k))** 0.16  ! Heymsfield and Donner
    711. 

  711.             if (k .eq. 1) then
    712. 

  712.                del_tv=amin1(del_tv,0.9*(zz(k)-topo(i,j))/vti(k))
    713. 

  713.             else
    714. 

  714.                del_tv=amin1(del_tv,0.9*(zz(k)-zz(k-1))/vti(k))
    715. 

  715.             endif
    716. 

  716.          else
    717. 

  717.             vti(k)=0.
    718. 

  718.          endif
    719. 

  719.       enddo
    720. 

  720. 

  721.       if (max_q .ge. min_q) then
    722. 

  722. !
    723. 

  723. !
    724. 

  724. !- Check if the summation of the small delta t >=  big delta t
    725. 

  725. !             (t_del_tv)          (del_tv)             (dtb)
    726. 

  726. 

  727.          t_del_tv=t_del_tv+del_tv
    728. 

  728. 

  729.          if ( t_del_tv .ge. dtb ) then
    730. 

  730.               notlast=.false.
    731. 

  731.               del_tv=dtb+del_tv-t_del_tv
    732. 

  732.          endif
    733. 

  733. 

  734. ! use small delta t to calculate the qiz flux
    735. 

  735. ! termi is the qiz flux pass in the grid box through the upper boundary
    736. 

  736. ! termo is the qiz flux pass out the grid box through the lower boundary
    737. 

  737. !
    738. 

  738. 

  739.          fluxin=0.
    740. 

  740.          do k=max_q,min_q,-1
    741. 

  741.             fluxout=rhoz(k)*vti(k)*qiz(k)
    742. 

  742.             flux=(fluxin-fluxout)/rhoz(k)/dzw(k)
    743. 

  743.             qiz(k)=qiz(k)+del_tv*flux
    744. 

  744.             qiz(k)=amax1(0.,qiz(k))
    745. 

  745.             qi(i,k,j)=qiz(k)
    746. 

  746.             fluxin=fluxout
    747. 

  747.          enddo
    748. 

  748.          if (min_q .eq. 1) then
    749. 

  749.             pptice=pptice+fluxin*del_tv
    750. 

  750.          else
    751. 

  751.             qiz(min_q-1)=qiz(min_q-1)+del_tv*  &
    752. 

  752.                          fluxin/rhoz(min_q-1)/dzw(min_q-1)
    753. 

  753.             qi(i,min_q-1,j)=qiz(min_q-1)
    754. 

  754.          endif
    755. 

  755. !
    756. 

  756.       else
    757. 

  757.          notlast=.false.
    758. 

  758.       endif
    759. 

  759. !
    760. 

  760.    ENDDO !notlast
    761. 

  761. 

  762. !   prnc(i,j)=prnc(i,j)+pptrain
    763. 

  763. !   psnowc(i,j)=psnowc(i,j)+pptsnow
    764. 

  764. !   pgrauc(i,j)=pgrauc(i,j)+pptgraul
    765. 

  765. !   picec(i,j)=picec(i,j)+pptice
    766. 

  766. !                     
    767. 

  767. 

  768. !   write(6,*) 'i=',i,' j=',j,'   ', pptrain, pptsnow, pptgraul, pptice
    769. 

  769. !   call flush(6)
    770. 

  770. 

  771.    snowncv(i,j) = pptsnow
    772. 

  772.    snownc(i,j) = snownc(i,j) + pptsnow
    773. 

  773.    graupelncv(i,j) = pptgraul
    774. 

  774.    graupelnc(i,j) = graupelnc(i,j) + pptgraul 
    775. 

  775.    RAINNCV(i,j) = pptrain + pptsnow + pptgraul + pptice                 
    776. 

  776.    RAINNC(i,j)  = RAINNC(i,j) + pptrain + pptsnow + pptgraul + pptice
    777. 

  777.    sr(i,j) = 0.
    778. 

  778.    if (RAINNCV(i,j) .gt. 0.) sr(i,j) = (pptsnow + pptgraul + pptice) / RAINNCV(i,j) 
    779. 

  779. 

  780.   ENDDO i_loop
    781. 

  781.   ENDDO j_loop
    782. 

  782. 

  783. !  if (itimestep.eq.6480) then
    784. 

  784. !     write(51,*) 'in the end of fallflux, itimestep=',itimestep
    785. 

  785. !     do j=jts,jte
    786. 

  786. !        do i=its,ite 
    787. 

  787. !           if (rainnc(i,j).gt.400.) then
    788. 

  788. !              write(50,*) 'i=',i,' j=',j,' rainnc=',rainnc
    789. 

  789. !           endif
    790. 

  790. !        enddo
    791. 

  791. !     enddo
    792. 

  792. !  endif
    793. 

  793. 

  794.   END SUBROUTINE fall_flux
    795. 

  795. 

  796. !----------------------------------------------------------------
    797. 

  797.    REAL FUNCTION ggamma(X)
    798. 

  798. !----------------------------------------------------------------
    799. 

  799.    IMPLICIT NONE
    800. 

  800. !----------------------------------------------------------------
    801. 

  801.       REAL, INTENT(IN   ) :: x
    802. 

  802.       REAL, DIMENSION(8)  :: B
    803. 

  803.       INTEGER             ::j, K1
    804. 

  804.       REAL                ::PF, G1TO2 ,TEMP
    805. 

  805. 

  806.       DATA B/-.577191652,.988205891,-.897056937,.918206857,  &
    807. 

  807.              -.756704078,.482199394,-.193527818,.035868343/
    808. 

  808. 

  809.       PF=1.
    810. 

  810.       TEMP=X
    811. 

  811.       DO 10 J=1,200
    812. 

  812.       IF (TEMP .LE. 2) GO TO 20
    813. 

  813.       TEMP=TEMP-1.
    814. 

  814.    10 PF=PF*TEMP
    815. 

  815.   100 FORMAT(//,5X,'module_gsfcgce: INPUT TO GAMMA FUNCTION TOO LARGE, X=',E12.5)
    816. 

  816.       WRITE(wrf_err_message,100)X
    817. 

  817.       CALL wrf_error_fatal(wrf_err_message)
    818. 

  818.    20 G1TO2=1.
    819. 

  819.       TEMP=TEMP - 1.
    820. 

  820.       DO 30 K1=1,8
    821. 

  821.    30 G1TO2=G1TO2 + B(K1)*TEMP**K1
    822. 

  822.       ggamma=PF*G1TO2
    823. 

  823. 

  824.       END FUNCTION ggamma
    825. 

  825. 

  826. !-----------------------------------------------------------------------
    827. 

  827. !c Correction of negative values  
    828. 

  828.    SUBROUTINE negcor ( X, rho, dz8w,                         &
    829. 

  829.                       ims,ime, jms,jme, kms,kme,              & ! memory dims
    830. 

  830.                       itimestep, ics,                         &
    831. 

  831.                       its,ite, jts,jte, kts,kte               ) ! tile   dims
    832. 

  832. !-----------------------------------------------------------------------
    833. 

  833.   REAL, DIMENSION( ims:ime , kms:kme , jms:jme ),                 &
    834. 

  834.         INTENT(INOUT) ::                                     X   
    835. 

  835.   REAL, DIMENSION( ims:ime , kms:kme , jms:jme ),                 &
    836. 

  836.         INTENT(IN   ) ::                              rho, dz8w  
    837. 

  837.   integer, INTENT(IN   ) ::                           itimestep, ics 
    838. 

  838. 

  839. !c Local variables
    840. 

  840. !  REAL, DIMENSION( kts:kte ) ::  Y1, Y2
    841. 

  841.   REAL   ::   A0, A1, A2
    842. 

  842. 

  843.   A1=0.
    844. 

  844.   A2=0.
    845. 

  845.   do k=kts,kte
    846. 

  846.      do j=jts,jte
    847. 

  847.         do i=its,ite
    848. 

  848.         A1=A1+max(X(i,k,j), 0.)*rho(i,k,j)*dz8w(i,k,j)
    849. 

  849.         A2=A2+max(-X(i,k,j), 0.)*rho(i,k,j)*dz8w(i,k,j)
    850. 

  850.         enddo
    851. 

  851.      enddo
    852. 

  852.   enddo
    853. 

  853. 

  854. !  A1=0.0
    855. 

  855. !  A2=0.0
    856. 

  856. !  do k=kts,kte
    857. 

  857. !     A1=A1+Y1(k)
    858. 

  858. !     A2=A2+Y2(k)
    859. 

  859. !  enddo
    860. 

  860. 

  861.   A0=0.0
    862. 

  862. 

  863.   if (A1.NE.0.0.and.A1.GT.A2) then 
    864. 

  864.      A0=(A1-A2)/A1
    865. 

  865. 

  866.   if (mod(itimestep,540).eq.0) then
    867. 

  867.      if (ics.eq.1) then
    868. 

  868.         write(61,*) 'kms=',kms,'  kme=',kme,'  kts=',kts,'  kte=',kte
    869. 

  869.         write(61,*) 'jms=',jms,'  jme=',jme,'  jts=',jts,'  jte=',jte 
    870. 

  870.         write(61,*) 'ims=',ims,'  ime=',ime,'  its=',its,'  ite=',ite 
    871. 

  871.      endif 
    872. 

  872.      if (ics.eq.1) then
    873. 

  873.          write(61,*) 'qv timestep=',itimestep
    874. 

  874.          write(61,*) '  A1=',A1,'   A2=',A2,'   A0=',A0
    875. 

  875.      else if (ics.eq.2) then
    876. 

  876.              write(61,*) 'ql timestep=',itimestep
    877. 

  877.              write(61,*) '  A1=',A1,'   A2=',A2,'   A0=',A0
    878. 

  878.      else if (ics.eq.3) then
    879. 

  879.              write(61,*) 'qr timestep=',itimestep
    880. 

  880.              write(61,*) '  A1=',A1,'   A2=',A2,'   A0=',A0
    881. 

  881.      else if (ics.eq.4) then
    882. 

  882.              write(61,*) 'qi timestep=',itimestep
    883. 

  883.              write(61,*) '  A1=',A1,'   A2=',A2,'   A0=',A0
    884. 

  884.      else if (ics.eq.5) then
    885. 

  885.              write(61,*) 'qs timestep=',itimestep
    886. 

  886.              write(61,*) '  A1=',A1,'   A2=',A2,'   A0=',A0
    887. 

  887.      else if (ics.eq.6) then
    888. 

  888.              write(61,*) 'qg timestep=',itimestep
    889. 

  889.              write(61,*) '  A1=',A1,'   A2=',A2,'   A0=',A0
    890. 

  890.      else
    891. 

  891.              write(61,*) 'wrong cloud specieis number'
    892. 

  892.      endif 
    893. 

  893.   endif 
    894. 

  894. 

  895.      do k=kts,kte
    896. 

  896.         do j=jts,jte
    897. 

  897.            do i=its,ite
    898. 

  898.            X(i,k,j)=A0*AMAX1(X(i,k,j), 0.0)
    899. 

  899.            enddo
    900. 

  900.         enddo
    901. 

  901.      enddo
    902. 

  902.   endif
    903. 

  903. 

  904.   END SUBROUTINE negcor
    905. 

  905. 

  906.  SUBROUTINE consat_s (ihail,itaobraun)
    907. 

  907. 

  908. !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    909. 

  909. !                                                                      c
    910. 

  910. !   Tao, W.-K., and J. Simpson, 1989: Modeling study of a tropical     c
    911. 

  911. !   squall-type convective line. J. Atmos. Sci., 46, 177-202.          c
    912. 

  912. !                                                                      c
    913. 

  913. !   Tao, W.-K., J. Simpson and M. McCumber, 1989: An ice-water         c
    914. 

  914. !   saturation adjustment. Mon. Wea. Rev., 117, 231-235.               c
    915. 

  915. !                                                                      c
    916. 

  916. !   Tao, W.-K., and J. Simpson, 1993: The Goddard Cumulus Ensemble     c
    917. 

  917. !   Model. Part I: Model description. Terrestrial, Atmospheric and     c
    918. 

  918. !   Oceanic Sciences, 4, 35-72.                                        c
    919. 

  919. !                                                                      c
    920. 

  920. !   Tao, W.-K., J. Simpson, D. Baker, S. Braun, M.-D. Chou, B.         c
    921. 

  921. !   Ferrier,D. Johnson, A. Khain, S. Lang,  B. Lynn, C.-L. Shie,       c
    922. 

  922. !   D. Starr, C.-H. Sui, Y. Wang and P. Wetzel, 2003: Microphysics,    c
    923. 

  923. !   radiation and surface processes in the Goddard Cumulus Ensemble    c
    924. 

  924. !   (GCE) model, A Special Issue on Non-hydrostatic Mesoscale          c
    925. 

  925. !   Modeling, Meteorology and Atmospheric Physics, 82, 97-137.         c
    926. 

  926. !                                                                      c
    927. 

  927. !   Lang, S., W.-K. Tao, R. Cifelli, W. Olson, J. Halverson, S.        c
    928. 

  928. !   Rutledge, and J. Simpson, 2007: Improving simulations of           c
    929. 

  929. !   convective system from TRMM LBA: Easterly and Westerly regimes.    c
    930. 

  930. !   J. Atmos. Sci., 64, 1141-1164.                                     c
    931. 

  931. !                                                                      c
    932. 

  932. !   Coded by Tao (1989-2003), modified by S. Lang (2006/07)            c
    933. 

  933. !                                                                      c
    934. 

  934. !   Implemented into WRF  by Roger Shi 2006/2007                       c
    935. 

  935. !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    936. 

  936. 

  937. !        itaobraun=0   ! see Tao and Simpson (1993)
    938. 

  938. !        itaobraun=1   ! see Tao et al. (2003)
    939. 

  939. 

  940.  integer :: itaobraun
    941. 

  941.  real    :: cn0
    942. 

  942. 

  943. !JJS 1/3/2008  vvvvv
    944. 

  944. !JJS   the following common blocks have been moved to the top of
    945. 

  945. !JJS   module_mp_gsfcgce_driver_instat.F
    946. 

  946. !
    947. 

  947. ! real,   dimension (1:31) ::  a1, a2
    948. 

  948. ! data a1/.7939e-7,.7841e-6,.3369e-5,.4336e-5,.5285e-5,.3728e-5, &
    949. 

  949. !         .1852e-5,.2991e-6,.4248e-6,.7434e-6,.1812e-5,.4394e-5,.9145e-5, &
    950. 

  950. !         .1725e-4,.3348e-4,.1725e-4,.9175e-5,.4412e-5,.2252e-5,.9115e-6, &
    951. 

  951. !         .4876e-6,.3473e-6,.4758e-6,.6306e-6,.8573e-6,.7868e-6,.7192e-6, &
    952. 

  952. !         .6513e-6,.5956e-6,.5333e-6,.4834e-6/
    953. 

  953. ! data a2/.4006,.4831,.5320,.5307,.5319,.5249,.4888,.3894,.4047, &
    954. 

  954. !         .4318,.4771,.5183,.5463,.5651,.5813,.5655,.5478,.5203,.4906, &
    955. 

  955. !         .4447,.4126,.3960,.4149,.4320,.4506,.4483,.4460,.4433,.4413, &
    956. 

  956. !         .4382,.4361/
    957. 

  957. !JJS 1/3/2008  ^^^^^
    958. 

  958. 

  959. 

  960. !     ******************************************************************
    961. 

  961. !JJS
    962. 

  962.       al = 2.5e10
    963. 

  963.       cp = 1.004e7
    964. 

  964.       rd1 = 1.e-3
    965. 

  965.       rd2 = 2.2
    966. 

  966. !JJS
    967. 

  967.       cpi=4.*atan(1.)
    968. 

  968.       cpi2=cpi*cpi
    969. 

  969.       grvt=980.
    970. 

  970.       cd1=6.e-1
    971. 

  971.       cd2=4.*grvt/(3.*cd1)
    972. 

  972.       tca=2.43e3
    973. 

  973.       dwv=.226
    974. 

  974.       dva=1.718e-4
    975. 

  975.       amw=18.016
    976. 

  976.       ars=8.314e7
    977. 

  977.       scv=2.2904487
    978. 

  978.       t0=273.16
    979. 

  979.       t00=238.16
    980. 

  980.       alv=2.5e10
    981. 

  981.       alf=3.336e9
    982. 

  982.       als=2.8336e10
    983. 

  983.       avc=alv/cp
    984. 

  984.       afc=alf/cp
    985. 

  985.       asc=als/cp
    986. 

  986.       rw=4.615e6
    987. 

  987.       cw=4.187e7
    988. 

  988.       ci=2.093e7
    989. 

  989.       c76=7.66
    990. 

  990.       c358=35.86
    991. 

  991.       c172=17.26939
    992. 

  992.       c409=4098.026
    993. 

  993.       c218=21.87456
    994. 

  994.       c580=5807.695
    995. 

  995.       c610=6.1078e3
    996. 

  996.       c149=1.496286e-5
    997. 

  997.       c879=8.794142
    998. 

  998.       c141=1.4144354e7
    999. 

  999. !***   DEFINE THE COEFFICIENTS USED IN TERMINAL VELOCITY
    1000. 

  1000. !***   DEFINE THE DENSITY AND SIZE DISTRIBUTION OF PRECIPITATION
    1001. 

  1001. !**********   HAIL OR GRAUPEL PARAMETERS   **********
    1002. 

  1002.       if(ihail .eq. 1) then
    1003. 

  1003.          roqg=.9
    1004. 

  1004.          ag=sqrt(cd2*roqg)
    1005. 

  1005.          bg=.5
    1006. 

  1006.          tng=.002
    1007. 

  1007.       else
    1008. 

  1008.          roqg=.4
    1009. 

  1009.          ag=351.2
    1010. 

  1010. !        AG=372.3 ! if ice913=1 6/15/02 tao's
    1011. 

  1011.          bg=.37
    1012. 

  1012.          tng=.04
    1013. 

  1013.       endif
    1014. 

  1014. !**********         SNOW PARAMETERS        **********
    1015. 

  1015. !ccshie 6/15/02 tao's
    1016. 

  1016. !      TNS=1.
    1017. 

  1017. !      TNS=.08 ! if ice913=1, tao's
    1018. 

  1018.       tns=.16 ! if ice913=0, tao's
    1019. 

  1019.       roqs=.1
    1020. 

  1020. !      AS=152.93
    1021. 

  1021.       as=78.63
    1022. 

  1022. !      BS=.25
    1023. 

  1023.       bs=.11
    1024. 

  1024. !**********         RAIN PARAMETERS        **********
    1025. 

  1025.       aw=2115.
    1026. 

  1026.       bw=.8
    1027. 

  1027.       roqr=1.
    1028. 

  1028.       tnw=.08
    1029. 

  1029. !*****************************************************************
    1030. 

  1030.       bgh=.5*bg
    1031. 

  1031.       bsh=.5*bs
    1032. 

  1032.       bwh=.5*bw
    1033. 

  1033.       bgq=.25*bg
    1034. 

  1034.       bsq=.25*bs
    1035. 

  1035.       bwq=.25*bw
    1036. 

  1036. !**********GAMMA FUNCTION CALCULATIONS*************
    1037. 

  1037.       ga3b  = gammagce(3.+bw)
    1038. 

  1038.       ga4b  = gammagce(4.+bw)
    1039. 

  1039.       ga6b  = gammagce(6.+bw)
    1040. 

  1040.       ga5bh = gammagce((5.+bw)/2.)
    1041. 

  1041.       ga3g  = gammagce(3.+bg)
    1042. 

  1042.       ga4g  = gammagce(4.+bg)
    1043. 

  1043.       ga5gh = gammagce((5.+bg)/2.)
    1044. 

  1044.       ga3d  = gammagce(3.+bs)
    1045. 

  1045.       ga4d  = gammagce(4.+bs)
    1046. 

  1046.       ga5dh = gammagce((5.+bs)/2.)
    1047. 

  1047. !CCCCC        LIN ET AL., 1983 OR LORD ET AL., 1984   CCCCCCCCCCCCCCCCC
    1048. 

  1048.       ac1=aw
    1049. 

  1049. !JJS
    1050. 

  1050.       ac2=ag
    1051. 

  1051.       ac3=as
    1052. 

  1052. !JJS
    1053. 

  1053.       bc1=bw
    1054. 

  1054.       cc1=as
    1055. 

  1055.       dc1=bs
    1056. 

  1056.       zrc=(cpi*roqr*tnw)**0.25
    1057. 

  1057.       zsc=(cpi*roqs*tns)**0.25
    1058. 

  1058.       zgc=(cpi*roqg*tng)**0.25
    1059. 

  1059.       vrc=aw*ga4b/(6.*zrc**bw)
    1060. 

  1060.       vsc=as*ga4d/(6.*zsc**bs)
    1061. 

  1061.       vgc=ag*ga4g/(6.*zgc**bg)
    1062. 

  1062. !     ****************************
    1063. 

  1063. !     RN1=1.E-3
    1064. 

  1064.       rn1=9.4e-15 ! 6/15/02 tao's
    1065. 

  1065.       bnd1=6.e-4
    1066. 

  1066.       rn2=1.e-3
    1067. 

  1067. !     BND2=1.25E-3
    1068. 

  1068. !     BND2=1.5E-3 ! if ice913=1 6/15/02 tao's
    1069. 

  1069.       bnd2=2.0e-3 ! if ice913=0 6/15/02 tao's
    1070. 

  1070.       rn3=.25*cpi*tns*cc1*ga3d
    1071. 

  1071.       esw=1.
    1072. 

  1072.       rn4=.25*cpi*esw*tns*cc1*ga3d
    1073. 

  1073. !     ERI=1.
    1074. 

  1074.       eri=.1  ! 6/17/02 tao's ice913=0 (not 1)
    1075. 

  1075.       rn5=.25*cpi*eri*tnw*ac1*ga3b
    1076. 

  1076. !     AMI=1./(24.*4.19E-10)
    1077. 

  1077.       ami=1./(24.*6.e-9) ! 6/15/02 tao's
    1078. 

  1078.       rn6=cpi2*eri*tnw*ac1*roqr*ga6b*ami
    1079. 

  1079. !     ESR=1. ! also if ice913=1 for tao's
    1080. 

  1080.       esr=.5 ! 6/15/02 for ice913=0 tao's
    1081. 

  1081.       rn7=cpi2*esr*tnw*tns*roqs
    1082. 

  1082.       esr=1.
    1083. 

  1083.       rn8=cpi2*esr*tnw*tns*roqr
    1084. 

  1084.       rn9=cpi2*tns*tng*roqs
    1085. 

  1085.       rn10=2.*cpi*tns
    1086. 

  1086.       rn101=.31*ga5dh*sqrt(cc1)
    1087. 

  1087.       rn10a=als*als/rw
    1088. 

  1088. !JJS
    1089. 

  1089.        rn10b=alv/tca
    1090. 

  1090.        rn10c=ars/(dwv*amw)
    1091. 

  1091. !JJS
    1092. 

  1092.       rn11=2.*cpi*tns/alf
    1093. 

  1093.       rn11a=cw/alf
    1094. 

  1094. !     AMI50=1.51e-7
    1095. 

  1095.       ami50=3.84e-6 ! 6/15/02 tao's
    1096. 

  1096. !     AMI40=2.41e-8
    1097. 

  1097.       ami40=3.08e-8 ! 6/15/02 tao's
    1098. 

  1098.       eiw=1.
    1099. 

  1099. !     UI50=20.
    1100. 

  1100.       ui50=100. ! 6/15/02 tao's
    1101. 

  1101.       ri50=2.*5.e-3
    1102. 

  1102.       cmn=1.05e-15
    1103. 

  1103.       rn12=cpi*eiw*ui50*ri50**2
    1104. 

  1104. 

  1105.       do 10 k=1,31
    1106. 

  1106.          y1=1.-aa2(k)
    1107. 

  1107.          rn13(k)=aa1(k)*y1/(ami50**y1-ami40**y1)
    1108. 

  1108.          rn12a(k)=rn13(k)/ami50
    1109. 

  1109.          rn12b(k)=aa1(k)*ami50**aa2(k)
    1110. 

  1110.          rn25a(k)=aa1(k)*cmn**aa2(k)
    1111. 

  1111.    10 continue
    1112. 

  1112. 

  1113.       egw=1.
    1114. 

  1114.       rn14=.25*cpi*egw*tng*ga3g*ag
    1115. 

  1115.       egi=.1
    1116. 

  1116.       rn15=.25*cpi*egi*tng*ga3g*ag
    1117. 

  1117.       egi=1.
    1118. 

  1118.       rn15a=.25*cpi*egi*tng*ga3g*ag
    1119. 

  1119.       egr=1.
    1120. 

  1120.       rn16=cpi2*egr*tng*tnw*roqr
    1121. 

  1121.       rn17=2.*cpi*tng
    1122. 

  1122.       rn17a=.31*ga5gh*sqrt(ag)
    1123. 

  1123.       rn17b=cw-ci
    1124. 

  1124.       rn17c=cw
    1125. 

  1125.       apri=.66
    1126. 

  1126.       bpri=1.e-4
    1127. 

  1127.       bpri=0.5*bpri ! 6/17/02 tao's
    1128. 

  1128.       rn18=20.*cpi2*bpri*tnw*roqr
    1129. 

  1129.       rn18a=apri
    1130. 

  1130.       rn19=2.*cpi*tng/alf
    1131. 

  1131.       rn19a=.31*ga5gh*sqrt(ag)
    1132. 

  1132.       rn19b=cw/alf
    1133. 

  1133. !
    1134. 

  1134.        rnn191=.78
    1135. 

  1135.        rnn192=.31*ga5gh*sqrt(ac2/dva)
    1136. 

  1136. !
    1137. 

  1137.       rn20=2.*cpi*tng
    1138. 

  1138.       rn20a=als*als/rw
    1139. 

  1139.       rn20b=.31*ga5gh*sqrt(ag)
    1140. 

  1140.       bnd3=2.e-3
    1141. 

  1141.       rn21=1.e3*1.569e-12/0.15
    1142. 

  1142.       erw=1.
    1143. 

  1143.       rn22=.25*cpi*erw*ac1*tnw*ga3b
    1144. 

  1144.       rn23=2.*cpi*tnw
    1145. 

  1145.       rn23a=.31*ga5bh*sqrt(ac1)
    1146. 

  1146.       rn23b=alv*alv/rw
    1147. 

  1147. 

  1148. 

  1149. !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    1150. 

  1150. !cc
    1151. 

  1151. !cc        "c0" in routine      "consat" (2d), "consatrh" (3d)
    1152. 

  1152. !cc        if ( itaobraun.eq.1 ) --> betah=0.5*beta=-.46*0.5=-0.23;   cn0=1.e-6
    1153. 

  1153. !cc        if ( itaobraun.eq.0 ) --> betah=0.5*beta=-.6*0.5=-0.30;    cn0=1.e-8
    1154. 

  1154. 

  1155.        if (itaobraun .eq. 0) then
    1156. 

  1156.          cn0=1.e-8
    1157. 

  1157.          beta=-.6
    1158. 

  1158.        elseif (itaobraun .eq. 1) then
    1159. 

  1159.          cn0=1.e-6
    1160. 

  1160.          beta=-.46
    1161. 

  1161.        endif
    1162. 

  1162. !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    1163. 

  1163. !      CN0=1.E-6
    1164. 

  1164. !      CN0=1.E-8 ! 6/15/02 tao's
    1165. 

  1165. !      BETA=-.46
    1166. 

  1166. !      BETA=-.6  ! 6/15/02 tao's
    1167. 

    1168. 

  1168.       rn25=cn0
    1169. 

  1169.       rn30a=alv*als*amw/(tca*ars)
    1170. 

  1170.       rn30b=alv/tca
    1171. 

  1171.       rn30c=ars/(dwv*amw)
    1172. 

  1172.       rn31=1.e-17
    1173. 

  1173. 

  1174.       rn32=4.*51.545e-4
    1175. 

  1175. !
    1176. 

  1176.       rn30=2.*cpi*tng
    1177. 

  1177.       rnn30a=alv*alv*amw/(tca*ars)
    1178. 

  1178. !
    1179. 

  1179.       rn33=4.*tns
    1180. 

  1180.        rn331=.65
    1181. 

  1181.        rn332=.44*sqrt(ac3/dva)*ga5dh
    1182. 

  1182. !
    1183. 

  1183. 

  1184.     return
    1185. 

  1185.  END SUBROUTINE consat_s
    1186. 

  1186. 

  1187.  SUBROUTINE saticel_s (dt, ihail, itaobraun, ice2, istatmin, &
    1188. 

  1188.                        new_ice_sat, id, &
    1189. 

  1189.                        ptwrf, qvwrf, qlwrf, qrwrf, &
    1190. 

  1190.                        qiwrf, qswrf, qgwrf, &
    1191. 

  1191.                        rho_mks, pi_mks, p0_mks,itimestep, &
    1192. 

  1192.                        ids,ide, jds,jde, kds,kde, &
    1193. 

  1193.                        ims,ime, jms,jme, kms,kme, &
    1194. 

  1194.                        its,ite, jts,jte, kts,kte  &
    1195. 

  1195.                            )
    1196. 

  1196.     IMPLICIT NONE
    1197. 

  1197. !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    1198. 

  1198. !                                                                      c
    1199. 

  1199. !   Tao, W.-K., and J. Simpson, 1989: Modeling study of a tropical     c
    1200. 

  1200. !   squall-type convective line. J. Atmos. Sci., 46, 177-202.          c
    1201. 

  1201. !                                                                      c
    1202. 

  1202. !   Tao, W.-K., J. Simpson and M. McCumber, 1989: An ice-water         c
    1203. 

  1203. !   saturation adjustment. Mon. Wea. Rev., 117, 231-235.               c
    1204. 

  1204. !                                                                      c
    1205. 

  1205. !   Tao, W.-K., and J. Simpson, 1993: The Goddard Cumulus Ensemble     c
    1206. 

  1206. !   Model. Part I: Model description. Terrestrial, Atmospheric and     c
    1207. 

  1207. !   Oceanic Sciences, 4, 35-72.                                        c
    1208. 

  1208. !                                                                      c
    1209. 

  1209. !   Tao, W.-K., J. Simpson, D. Baker, S. Braun, M.-D. Chou, B.         c
    1210. 

  1210. !   Ferrier,D. Johnson, A. Khain, S. Lang,  B. Lynn, C.-L. Shie,       c
    1211. 

  1211. !   D. Starr, C.-H. Sui, Y. Wang and P. Wetzel, 2003: Microphysics,    c
    1212. 

  1212. !   radiation and surface processes in the Goddard Cumulus Ensemble    c
    1213. 

  1213. !   (GCE) model, A Special Issue on Non-hydrostatic Mesoscale          c
    1214. 

  1214. !   Modeling, Meteorology and Atmospheric Physics, 82, 97-137.         c
    1215. 

  1215. !                                                                      c
    1216. 

  1216. !   Lang, S., W.-K. Tao, R. Cifelli, W. Olson, J. Halverson, S.        c
    1217. 

  1217. !   Rutledge, and J. Simpson, 2007: Improving simulations of           c
    1218. 

  1218. !   convective system from TRMM LBA: Easterly and Westerly regimes.    c
    1219. 

  1219. !   J. Atmos. Sci., 64, 1141-1164.                                     c
    1220. 

  1220. !                                                                      c
    1221. 

  1221. !   Tao, W.-K., J. J. Shi,  S. Lang, C. Peters-Lidard, A. Hou, S.      c
    1222. 

  1222. !   Braun, and J. Simpson, 2007: New, improved bulk-microphysical      c
    1223. 

  1223. !   schemes for studying precipitation processes in WRF. Part I:       c
    1224. 

  1224. !   Comparisons with other schemes. to appear on Mon. Wea. Rev.        C
    1225. 

  1225. !                                                                      c
    1226. 

  1226. !   Coded by Tao (1989-2003), modified by S. Lang (2006/07)            c
    1227. 

  1227. !                                                                      c
    1228. 

  1228. !   Implemented into WRF  by Roger Shi 2006/2007                       c
    1229. 

  1229. !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    1230. 

  1230. !
    1231. 

  1231. !      COMPUTE ICE PHASE MICROPHYSICS AND SATURATION PROCESSES
    1232. 

  1232. !
    1233. 

  1233.   integer,    parameter ::  nt=2880, nt2=2*nt 
    1234. 

  1234. 

  1235. !cc   using scott braun's way for pint, pidep computations
    1236. 

  1236.   integer  ::   itaobraun,ice2,ihail,new_ice_sat,id,istatmin
    1237. 

  1237.   integer  ::   itimestep
    1238. 

  1238.   real     ::   tairccri, cn0, dt
    1239. 

  1239. !cc
    1240. 

  1240. 

  1241. !JJS      common/bxyz/ n,isec,nran,kt1,kt2
    1242. 

  1242. !JJS      common/option/ lipps,ijkadv,istatmin,iwater,itoga,imlifting,lin,
    1243. 

  1243. !JJS     1   irf,iadvh,irfg,ismg,id
    1244. 

  1244. 

  1245. !JJS      common/timestat/ ndt_stat,idq
    1246. 

  1246. !JJS      common/iice/ new_ice_sat
    1247. 

  1247. !JJS      common/bt/ dt,d2t,rijl2,dts,f5,rd1,rd2,bound,al,cp,ra,ck,ce,eps,
    1248. 

  1248. !JJS     1    psfc,fcor,sec,aminut,rdt
    1249. 

  1249. 

  1250. !JJS   the following common blocks have been moved to the top of 
    1251. 

  1251. !JJS   module_mp_gsfcgce_driver_instat.F
    1252. 

  1252. 

  1253. !      common/bt/ rd1,rd2,al,cp
    1254. 

  1254. !
    1255. 

  1255. !
    1256. 

  1256. !      common/bterv/ zrc,zgc,zsc,vrc,vgc,vsc
    1257. 

  1257. !      common/size/ tnw,tns,tng,roqs,roqg,roqr
    1258. 

  1258. !      common/cont/ c38,c358,c610,c149,c879,c172,c409,c76,c218,c580,c141
    1259. 

  1259. !        common/b3cs/ ag,bg,as,bs,aw,bw,bgh,bgq,bsh,bsq,bwh,bwq
    1260. 

  1260. !      common/bsnw/ alv,alf,als,t0,t00,avc,afc,asc,rn1,bnd1,rn2,bnd2, &
    1261. 

  1261. !         rn3,rn4,rn5,rn6,rn7,rn8,rn9,rn10,rn101,rn10a,rn11,rn11a, &
    1262. 

  1262. !         rn12,rn12a(31),rn12b(31),rn13(31),rn14,rn15,rn15a,rn16,rn17, &
    1263. 

  1263. !         rn17a,rn17b,rn17c,rn18,rn18a,rn19,rn19a,rn19b,rn20,rn20a,rn20b, &
    1264. 

  1264. !         bnd3,rn21,rn22,rn23,rn23a,rn23b,rn25,rn25a(31),rn30a,rn30b, &
    1265. 

  1265. !         rn30c,rn31,beta,rn32
    1266. 

  1266. !      common/rsnw1/ rn10b,rn10c,rnn191,rnn192,rn30,rnn30a,rn33,rn331, &
    1267. 

  1267. !                    rn332
    1268. 

  1268. !JJS
    1269. 

  1269. 

  1270.   integer ids,ide,jds,jde,kds,kde
    1271. 

  1271.   integer ims,ime,jms,jme,kms,kme
    1272. 

  1272.   integer its,ite,jts,jte,kts,kte
    1273. 

  1273.   integer i,j,k, kp
    1274. 

  1274. 

  1275.   real :: a0 ,a1 ,a2 ,afcp ,alvr ,ami100 ,ami40 ,ami50 ,ascp ,avcp ,betah &
    1276. 

  1276.    ,bg3 ,bgh5 ,bs3 ,bs6 ,bsh5 ,bw3 ,bw6 ,bwh5 ,cmin ,cmin1 ,cmin2 ,cp409 &
    1277. 

  1277.    ,cp580 ,cs580 ,cv409 ,d2t ,del ,dwvp ,ee1 ,ee2 ,f00 ,f2 ,f3 ,ft ,fv0 ,fvs &
    1278. 

  1278.    ,pi0 ,pir ,pr0 ,qb0 ,r00 ,r0s ,r101f ,r10ar ,r10t ,r11at ,r11rt ,r12r ,r14f &
    1279. 

  1279.    ,r14r ,r15af ,r15ar ,r15f ,r15r ,r16r ,r17aq ,r17as ,r17r ,r18r ,r19aq ,r19as &
    1280. 

  1280.    ,r19bt ,r19rt ,r20bq ,r20bs ,r20t ,r22f ,r23af ,r23br ,r23t ,r25a ,r25rt ,r2ice &
    1281. 

  1281.    ,r31r ,r32rt ,r3f ,r4f ,r5f ,r6f ,r7r ,r8r ,r9r ,r_nci ,rft ,rijl2 ,rp0 ,rr0 &
    1282. 

  1282.    ,rrq ,rrs ,rt0 ,scc ,sccc ,sddd ,see ,seee ,sfff ,smmm ,ssss ,tb0 ,temp ,ucog &
    1283. 

  1283.    ,ucor ,ucos ,uwet ,vgcf ,vgcr ,vrcf ,vscf ,zgr ,zrr ,zsr
    1284. 

  1284. 

  1285. 

  1286.   real, dimension (its:ite,jts:jte,kts:kte) ::  fv
    1287. 

  1287.   real, dimension (its:ite,jts:jte,kts:kte) ::  dpt, dqv
    1288. 

  1288.   real, dimension (its:ite,jts:jte,kts:kte) ::  qcl, qrn,      &
    1289. 

  1289.                                                 qci, qcs, qcg
    1290. 

  1290. !JJS 10/16/06    vvvv
    1291. 

  1291. !      real dpt1(ims:ime,jms:jme,kms:kme)
    1292. 

  1292. !      real dqv1(ims:ime,jms:jme,kms:kme),
    1293. 

  1293. !     1             qcl1(ims:ime,jms:jme,kms:kme)
    1294. 

  1294. !      real qrn1(ims:ime,jms:jme,kms:kme),
    1295. 

  1295. !     1             qci1(ims:ime,jms:jme,kms:kme)
    1296. 

  1296. !      real qcs1(ims:ime,jms:jme,kms:kme),
    1297. 

  1297. !     1             qcg1(ims:ime,jms:jme,kms:kme)
    1298. 

  1298. !JJS 10/16/06    ^^^^
    1299. 

  1299. 

  1300. !JJS
    1301. 

  1301. 

  1302.   real, dimension (ims:ime, kms:kme, jms:jme) ::  ptwrf, qvwrf 
    1303. 

  1303.   real, dimension (ims:ime, kms:kme, jms:jme) ::  qlwrf, qrwrf,        &
    1304. 

  1304.                                                   qiwrf, qswrf, qgwrf
    1305. 

  1305. !JJS 10/16/06    vvvv
    1306. 

  1306. !      real ptwrfold(ims:ime, kms:kme, jms:jme)
    1307. 

  1307. !      real qvwrfold(ims:ime, kms:kme, jms:jme),
    1308. 

  1308. !     1             qlwrfold(ims:ime, kms:kme, jms:jme)
    1309. 

  1309. !      real qrwrfold(ims:ime, kms:kme, jms:jme),
    1310. 

  1310. !     1             qiwrfold(ims:ime, kms:kme, jms:jme)
    1311. 

  1311. !      real qswrfold(ims:ime, kms:kme, jms:jme),
    1312. 

  1312. !     1             qgwrfold(ims:ime, kms:kme, jms:jme)
    1313. 

  1313. !JJS 10/16/06    ^^^^
    1314. 

  1314. 

  1315. !JJS in MKS
    1316. 

  1316.   real, dimension (ims:ime, kms:kme, jms:jme) ::  rho_mks
    1317. 

  1317.   real, dimension (ims:ime, kms:kme, jms:jme) ::  pi_mks
    1318. 

  1318.   real, dimension (ims:ime, kms:kme, jms:jme) ::  p0_mks
    1319. 

  1319. !JJS
    1320. 

  1320. !  real, dimension (its:ite,jts:jte,kts:kte) ::  ww1
    1321. 

  1321. !  real, dimension (its:ite,jts:jte,kts:kte) ::  rsw
    1322. 

  1322. !  real, dimension (its:ite,jts:jte,kts:kte) ::  rlw
    1323. 

  1323. 

  1324. !JJS      COMMON /BADV/
    1325. 

  1325.   real, dimension (its:ite,jts:jte) ::        &
    1326. 

  1326.            vg,      zg,       &
    1327. 

  1327.            ps,      pg,       &
    1328. 

  1328.           prn,     psn,       &
    1329. 

  1329.         pwacs,   wgacr,       &
    1330. 

  1330.         pidep,    pint,       &
    1331. 

  1331.           qsi,     ssi,       &
    1332. 

  1332.           esi,     esw,       &
    1333. 

  1333.           qsw,      pr,       &
    1334. 

  1334.           ssw,   pihom,       &
    1335. 

  1335.          pidw,   pimlt,       &
    1336. 

  1336.         psaut,   qracs,       &
    1337. 

  1337.         psaci,   psacw,       &
    1338. 

  1338.         qsacw,   praci,       &
    1339. 

  1339.         pmlts,   pmltg,       &
    1340. 

  1340.         asss,       y1,    y2
    1341. 

  1341. !JJS       Y2(its:ite,jts:jte),    DDE(NB)
    1342. 

  1342. 

  1343. !JJS      COMMON/BSAT/
    1344. 

  1344.   real, dimension (its:ite,jts:jte) ::        &
    1345. 

  1345.         praut,   pracw,       &
    1346. 

  1346.          psfw,    psfi,       &
    1347. 

  1347.         dgacs,   dgacw,       &
    1348. 

  1348.         dgaci,   dgacr,       &
    1349. 

  1349.         pgacs,   wgacs,       &
    1350. 

  1350.         qgacw,   wgaci,       &
    1351. 

  1351.         qgacr,   pgwet,       &
    1352. 

  1352.         pgaut,   pracs,       &
    1353. 

  1353.         psacr,   qsacr,       &
    1354. 

  1354.          pgfr,   psmlt,       &
    1355. 

  1355.         pgmlt,   psdep,       &
    1356. 

  1356.         pgdep,   piacr,       &
    1357. 

  1357.            y5,     scv,       &
    1358. 

  1358.           tca,     dwv,       &
    1359. 

  1359.           egs,      y3,       &
    1360. 

  1360.            y4,     ddb
    1361. 

  1361. 

  1362. !JJS      COMMON/BSAT1/
    1363. 

  1363.   real, dimension (its:ite,jts:jte) ::        &
    1364. 

  1364.            pt,      qv,       &
    1365. 

  1365.            qc,      qr,       &
    1366. 

  1366.            qi,      qs,       &
    1367. 

  1367.            qg,    tair,       &
    1368. 

  1368.         tairc,   rtair,       &
    1369. 

  1369.           dep,      dd,       &
    1370. 

  1370.           dd1,     qvs,       &
    1371. 

  1371.            dm,      rq,       &
    1372. 

  1372.         rsub1,     col,       &
    1373. 

  1373.           cnd,     ern,       &
    1374. 

  1374.          dlt1,    dlt2,       &
    1375. 

  1375.          dlt3,    dlt4,       &
    1376. 

  1376.            zr,      vr,       &
    1377. 

  1377.            zs,      vs,       &
    1378. 

  1378.                  pssub,       &
    1379. 

  1379.         pgsub,     dda
    1380. 

  1380. 

  1381. !JJS      COMMON/B5/
    1382. 

  1382.   real, dimension (its:ite,jts:jte,kts:kte) ::  rho
    1383. 

  1383.   real, dimension (kts:kte) ::                 & 
    1384. 

  1384.            tb,      qb,    rho1,              &
    1385. 

  1385.            ta,      qa,     ta1,     qa1,     &
    1386. 

  1386.          coef,      z1,      z2,      z3,     &
    1387. 

  1387.            am,     am1,      ub,      vb,     &
    1388. 

  1388.            wb,     ub1,     vb1,    rrho,     &
    1389. 

  1389.         rrho1,     wbx
    1390. 

  1390. 

  1391. !JJS      COMMON/B6/
    1392. 

  1392.   real, dimension (its:ite,jts:jte,kts:kte) ::  p0, pi, f0
    1393. 

  1393.   real, dimension (kts:kte) ::    & 
    1394. 

  1394.            fd,      fe,        &
    1395. 

  1395.            st,      sv,        &
    1396. 

  1396.            sq,      sc,        &
    1397. 

  1397.            se,     sqa
    1398. 

  1398. 

  1399. !JJS      COMMON/BRH1/
    1400. 

  1400.   real, dimension (kts:kte) ::    & 
    1401. 

  1401.          srro,    qrro,    sqc,    sqr,    &
    1402. 

  1402.           sqi,     sqs,    sqg,   stqc,    &
    1403. 

  1403.          stqr,    stqi,   stqs,   stqg
    1404. 

  1404.   real, dimension (nt) ::    & 
    1405. 

  1405.           tqc,     tqr,    tqi,    tqs,    tqg
    1406. 

  1406. 

  1407. !JJS     common/bls/ y0(nx,ny),ts0new(nx,ny),qss0new(nx,ny)
    1408. 

  1408. 

  1409. !JJS      COMMON/BLS/
    1410. 

  1410.   real, dimension (ims:ime,jms:jme) ::     &
    1411. 

  1411.            y0,     ts0,   qss0
    1412. 

  1412. 

  1413. !JJS      COMMON/BI/ IT(its:ite,jts:jte), ICS(its:ite,jts:jte,4)
    1414. 

  1414.   integer, dimension (its:ite,jts:jte) ::        it  
    1415. 

  1415.   integer, dimension (its:ite,jts:jte, 4) ::    ics 
    1416. 

  1416. 

  1417.   integer :: i24h
    1418. 

  1418.   integer :: iwarm
    1419. 

  1419.   real :: r2is, r2ig
    1420. 

  1420.   
    1421. 

  1421. 

  1422. !JJS      COMMON/MICRO/
    1423. 

  1423. !  real, dimension (ims:ime,kms:kme,jms:jme)  ::    dbz 
    1424. 

  1424. 

  1425. !23456789012345678901234567890123456789012345678901234567890123456789012
    1426. 

  1426. 

  1427. !
    1428. 

  1428. !JJS 1/3/2008  vvvvv
    1429. 

  1429. !JJS   the following common blocks have been moved to the top of
    1430. 

  1430. !JJS   module_mp_gsfcgce_driver.F
    1431. 

  1431. 

  1432. !  real, dimension (31)   ::      aa1,  aa2
    1433. 

  1433. !  data aa1/.7939e-7, .7841e-6, .3369e-5, .4336e-5, .5285e-5,     &
    1434. 

  1434. !           .3728e-5, .1852e-5, .2991e-6, .4248e-6, .7434e-6,     &
    1435. 

  1435. !           .1812e-5, .4394e-5, .9145e-5, .1725e-4, .3348e-4,     &
    1436. 

  1436. !           .1725e-4, .9175e-5, .4412e-5, .2252e-5, .9115e-6,     &
    1437. 

  1437. !           .4876e-6, .3473e-6, .4758e-6, .6306e-6, .8573e-6,     &
    1438. 

  1438. !           .7868e-6, .7192e-6, .6513e-6, .5956e-6, .5333e-6,     &
    1439. 

  1439. !           .4834e-6/
    1440. 

  1440. !  data aa2/.4006, .4831, .5320, .5307, .5319,      &
    1441. 

  1441. !           .5249, .4888, .3894, .4047, .4318,      &
    1442. 

  1442. !           .4771, .5183, .5463, .5651, .5813,      &
    1443. 

  1443. !           .5655, .5478, .5203, .4906, .4447,      &
    1444. 

  1444. !           .4126, .3960, .4149, .4320, .4506,      &
    1445. 

  1445. !           .4483, .4460, .4433, .4413, .4382,      &
    1446. 

  1446. !           .4361/
    1447. 

  1447. 

  1448. !JJS 1/3/2008  ^^^^^
    1449. 

  1449. 

  1450. !
    1451. 

  1451. !jm 20090220      save
    1452. 

  1452. 

  1453. !    i24h=nint(86400./dt) 
    1454. 

  1454. !    if (mod(itimestep,i24h).eq.1) then
    1455. 

  1455. !       write(6, *) 'ims=', ims, '  ime=', ime
    1456. 

  1456. !       write(6, *) 'jms=', jms, '  jme=', jme
    1457. 

  1457. !       write(6, *) 'kms=', kms, '  kme=', kme
    1458. 

  1458. !       write(6, *) 'its=', its, '  ite=', ite
    1459. 

  1459. !       write(6, *) 'jts=', jts, '  jte=', jte
    1460. 

  1460. !       write(6, *) 'kts=', kts, '  kte=', kte
    1461. 

  1461. !       write(6, *) '    ihail=', ihail
    1462. 

  1462. !       write(6, *) 'itaobraun=',itaobraun
    1463. 

  1463. !       write(6, *) '    ice2=', ICE2
    1464. 

  1464. !       write(6, *) 'istatmin=',istatmin
    1465. 

  1465. !       write(6, *) 'new_ice_sat=', new_ice_sat
    1466. 

  1466. !       write(6, *) 'id=', id
    1467. 

  1467. !       write(6, *) 'dt=', dt
    1468. 

  1468. !    endif
    1469. 

  1469. 

  1470. !JJS  convert from mks to cgs, and move from WRF grid to GCE grid
    1471. 

  1471.       do k=kts,kte
    1472. 

  1472.          do j=jts,jte
    1473. 

  1473.          do i=its,ite
    1474. 

  1474.          rho(i,j,k)=rho_mks(i,k,j)*0.001
    1475. 

  1475.          p0(i,j,k)=p0_mks(i,k,j)*10.0
    1476. 

  1476.          pi(i,j,k)=pi_mks(i,k,j)
    1477. 

  1477.          dpt(i,j,k)=ptwrf(i,k,j)
    1478. 

  1478.          dqv(i,j,k)=qvwrf(i,k,j)
    1479. 

  1479.          qcl(i,j,k)=qlwrf(i,k,j)
    1480. 

  1480.          qrn(i,j,k)=qrwrf(i,k,j)
    1481. 

  1481.          qci(i,j,k)=qiwrf(i,k,j)
    1482. 

  1482.          qcs(i,j,k)=qswrf(i,k,j)
    1483. 

  1483.          qcg(i,j,k)=qgwrf(i,k,j)
    1484. 

  1484. !JJS 10/16/06    vvvv
    1485. 

  1485. !         dpt1(i,j,k)=ptwrfold(i,k,j)
    1486. 

  1486. !         dqv1(i,j,k)=qvwrfold(i,k,j)
    1487. 

  1487. !         qcl1(i,j,k)=qlwrfold(i,k,j)
    1488. 

  1488. !         qrn1(i,j,k)=qrwrfold(i,k,j)
    1489. 

  1489. !         qci1(i,j,k)=qiwrfold(i,k,j)
    1490. 

  1490. !         qcs1(i,j,k)=qswrfold(i,k,j)
    1491. 

  1491. !         qcg1(i,j,k)=qgwrfold(i,k,j)
    1492. 

  1492. !JJS 10/16/06     ^^^^
    1493. 

  1493.          enddo !i
    1494. 

  1494.          enddo !j
    1495. 

  1495.       enddo !k
    1496. 

  1496. 

  1497.       do k=kts,kte
    1498. 

  1498.          do j=jts,jte
    1499. 

  1499.          do i=its,ite
    1500. 

  1500.          fv(i,j,k)=sqrt(rho(i,j,2)/rho(i,j,k))
    1501. 

  1501.          enddo !i
    1502. 

  1502.          enddo !j
    1503. 

  1503.       enddo !k
    1504. 

  1504. !JJS
    1505. 

  1505. 

  1506. !
    1507. 

  1507. !     ******   THREE CLASSES OF ICE-PHASE   (LIN ET AL, 1983)  *********
    1508. 

  1508. 

  1509. !JJS       D22T=D2T
    1510. 

  1510. !JJS       IF(IJKADV .EQ. 0) THEN
    1511. 

  1511. !JJS         D2T=D2T
    1512. 

  1512. !JJS       ELSE
    1513. 

  1513.          d2t=dt
    1514. 

  1514. !JJS       ENDIF
    1515. 

  1515. !
    1516. 

  1516. !       itaobraun=0 ! original pint and pidep & see Tao and Simpson 1993
    1517. 

  1517.         itaobraun=1 ! see Tao et al. (2003)
    1518. 

  1518. !
    1519. 

  1519.        if ( itaobraun.eq.0 ) then
    1520. 

  1520.           cn0=1.e-8
    1521. 

  1521. !c        beta=-.6
    1522. 

  1522.        elseif ( itaobraun.eq.1 ) then
    1523. 

  1523.           cn0=1.e-6
    1524. 

  1524. !         cn0=1.e-8  ! special
    1525. 

  1525. !c        beta=-.46
    1526. 

  1526.        endif
    1527. 

  1527. !C  TAO 2007 START
    1528. 

  1528. !   ICE2=0 ! default, 3ice with loud ice, snow and graupel
    1529. 

  1529. !              r2is=1., r2ig=1.
    1530. 

  1530. !   ICE2=1 ! 2ice with cloud ice and snow (no graupel) - r2iceg=1, r2ice=0.
    1531. 

  1531. !              r2is=1., r2ig=0.
    1532. 

  1532. !   ICE2=2 ! 2ice with cloud ice and graupel (no snow) - r2ice=1, r2iceg=0.
    1533. 

  1533. !              r2is=0., r2ig=1.
    1534. 

  1534. !c
    1535. 

  1535. !        r2ice=1.
    1536. 

  1536. !        r2iceg=1.
    1537. 

  1537.          r2ig=1.
    1538. 

  1538.          r2is=1.
    1539. 

  1539.           if (ice2 .eq. 1) then
    1540. 

  1540. !              r2ice=0.
    1541. 

  1541. !              r2iceg=1.
    1542. 

  1542.               r2ig=0.
    1543. 

  1543.               r2is=1.
    1544. 

  1544.           endif
    1545. 

  1545.           if (ice2 .eq. 2) then
    1546. 

  1546. !              r2ice=1.
    1547. 

  1547. !              r2iceg=0.
    1548. 

  1548.               r2ig=1.
    1549. 

  1549.               r2is=0.
    1550. 

  1550.           endif
    1551. 

  1551. !C  TAO 2007 END
    1552. 

  1552.      
    1553. 

  1553. !JJS  10/7/2008
    1554. 

  1554. !   ICE2=3 ! no ice, warm rain only
    1555. 

  1555.     iwarm = 0
    1556. 

  1556.     if (ice2 .eq. 3 ) iwarm = 1
    1557. 

  1557. 

  1558. 

  1559. 

  1560.       cmin=1.e-19
    1561. 

  1561.       cmin1=1.e-20
    1562. 

  1562.       cmin2=1.e-12
    1563. 

  1563.       ucor=3071.29/tnw**0.75
    1564. 

  1564.       ucos=687.97*roqs**0.25/tns**0.75
    1565. 

  1565.       ucog=687.97*roqg**0.25/tng**0.75
    1566. 

  1566.       uwet=4.464**0.95
    1567. 

  1567. 

  1568.       rijl2 = 1. / (ide-ids) / (jde-jds)
    1569. 

  1569. 

  1570. !JJScap $doacross local(j,i)
    1571. 

  1571. 

  1572. !JJS      DO 1 J=1,JMAX
    1573. 

  1573. !JJS      DO 1 I=1,IMAX
    1574. 

  1574.        do j=jts,jte
    1575. 

  1575.           do i=its,ite
    1576. 

  1576.           it(i,j)=1
    1577. 

  1577.           enddo
    1578. 

  1578.        enddo
    1579. 

  1579. 

  1580.       f2=rd1*d2t
    1581. 

  1581.       f3=rd2*d2t
    1582. 

  1582. 

  1583.       ft=dt/d2t
    1584. 

  1584.       rft=rijl2*ft
    1585. 

  1585.       a0=.5*istatmin*rijl2
    1586. 

  1586.       rt0=1./(t0-t00)
    1587. 

  1587.       bw3=bw+3.
    1588. 

  1588.       bs3=bs+3.
    1589. 

  1589.       bg3=bg+3.
    1590. 

  1590.       bsh5=2.5+bsh
    1591. 

  1591.       bgh5=2.5+bgh
    1592. 

  1592.       bwh5=2.5+bwh
    1593. 

  1593.       bw6=bw+6.
    1594. 

  1594.       bs6=bs+6.
    1595. 

  1595.       betah=.5*beta
    1596. 

  1596.       r10t=rn10*d2t
    1597. 

  1597.       r11at=rn11a*d2t
    1598. 

  1598.       r19bt=rn19b*d2t
    1599. 

  1599.       r20t=-rn20*d2t
    1600. 

  1600.       r23t=-rn23*d2t
    1601. 

  1601.       r25a=rn25
    1602. 

  1602. 

  1603. !     ami50 for use in PINT
    1604. 

  1604.        ami50=3.76e-8
    1605. 

  1605.        ami100=1.51e-7
    1606. 

  1606.        ami40=2.41e-8
    1607. 

  1607. 

  1608. !C    ******************************************************************
    1609. 

  1609. 

  1610. !JJS      DO 1000 K=2,kles
    1611. 

  1611.       do 1000 k=kts,kte
    1612. 

  1612.          kp=k+1
    1613. 

  1613. !JJS         tb0=ta1(k)
    1614. 

  1614. !JJS         qb0=qa1(k)
    1615. 

  1615.          tb0=0.
    1616. 

  1616.          qb0=0.
    1617. 

  1617. 

  1618.       do 2000 j=jts,jte
    1619. 

  1619.          do 2000 i=its,ite
    1620. 

  1620. 

  1621.          rp0=3.799052e3/p0(i,j,k)
    1622. 

  1622.          pi0=pi(i,j,k)
    1623. 

  1623.          pir=1./(pi(i,j,k))
    1624. 

  1624.          pr0=1./p0(i,j,k)
    1625. 

  1625.          r00=rho(i,j,k)
    1626. 

  1626.          r0s=sqrt(rho(i,j,k))
    1627. 

  1627. !JJS         RR0=RRHO(K)
    1628. 

  1628.          rr0=1./rho(i,j,k)
    1629. 

  1629. !JJS         RRS=SRRO(K)
    1630. 

  1630.          rrs=sqrt(rr0)
    1631. 

  1631. !JJS         RRQ=QRRO(K)
    1632. 

  1632.          rrq=sqrt(rrs)
    1633. 

  1633.          f0(i,j,k)=al/cp/pi(i,j,k)
    1634. 

  1634.          f00=f0(i,j,k)
    1635. 

  1635.          fv0=fv(i,j,k)
    1636. 

  1636.          fvs=sqrt(fv(i,j,k))
    1637. 

  1637.          zrr=1.e5*zrc*rrq
    1638. 

  1638.          zsr=1.e5*zsc*rrq
    1639. 

  1639.          zgr=1.e5*zgc*rrq
    1640. 

  1640.          cp409=c409*pi0
    1641. 

  1641.          cv409=c409*avc
    1642. 

  1642.          cp580=c580*pi0
    1643. 

  1643.          cs580=c580*asc
    1644. 

  1644.          alvr=r00*alv
    1645. 

  1645.          afcp=afc*pir
    1646. 

  1646.          avcp=avc*pir
    1647. 

  1647.          ascp=asc*pir
    1648. 

  1648.          vrcf=vrc*fv0
    1649. 

  1649.          vscf=vsc*fv0
    1650. 

  1650.          vgcf=vgc*fv0
    1651. 

  1651.          vgcr=vgc*rrs
    1652. 

  1652.          dwvp=c879*pr0
    1653. 

  1653.          r3f=rn3*fv0
    1654. 

  1654.          r4f=rn4*fv0
    1655. 

  1655.          r5f=rn5*fv0
    1656. 

  1656.          r6f=rn6*fv0
    1657. 

  1657.          r7r=rn7*rr0
    1658. 

  1658.          r8r=rn8*rr0
    1659. 

  1659.          r9r=rn9*rr0
    1660. 

  1660.          r101f=rn101*fvs
    1661. 

  1661.          r10ar=rn10a*r00
    1662. 

  1662.          r11rt=rn11*rr0*d2t
    1663. 

  1663.          r12r=rn12*r00
    1664. 

  1664.          r14r=rn14*rrs
    1665. 

  1665.          r14f=rn14*fv0
    1666. 

  1666.          r15r=rn15*rrs
    1667. 

  1667.          r15ar=rn15a*rrs
    1668. 

  1668.          r15f=rn15*fv0
    1669. 

  1669.          r15af=rn15a*fv0
    1670. 

  1670.          r16r=rn16*rr0
    1671. 

  1671.          r17r=rn17*rr0
    1672. 

  1672.          r17aq=rn17a*rrq
    1673. 

  1673.          r17as=rn17a*fvs
    1674. 

  1674.          r18r=rn18*rr0
    1675. 

  1675.          r19rt=rn19*rr0*d2t
    1676. 

  1676.          r19aq=rn19a*rrq
    1677. 

  1677.          r19as=rn19a*fvs
    1678. 

  1678.          r20bq=rn20b*rrq
    1679. 

  1679.          r20bs=rn20b*fvs
    1680. 

  1680.          r22f=rn22*fv0
    1681. 

  1681.          r23af=rn23a*fvs
    1682. 

  1682.          r23br=rn23b*r00
    1683. 

  1683.          r25rt=rn25*rr0*d2t
    1684. 

  1684.          r31r=rn31*rr0
    1685. 

  1685.          r32rt=rn32*d2t*rrs
    1686. 

  1686. 

  1687. !JJS       DO 100 J=3,JLES
    1688. 

  1688. !JJS       DO 100 I=3,ILES
    1689. 

  1689.         pt(i,j)=dpt(i,j,k)
    1690. 

  1690.         qv(i,j)=dqv(i,j,k)
    1691. 

  1691.         qc(i,j)=qcl(i,j,k)
    1692. 

  1692.         qr(i,j)=qrn(i,j,k)
    1693. 

  1693.         qi(i,j)=qci(i,j,k)
    1694. 

  1694.         qs(i,j)=qcs(i,j,k)
    1695. 

  1695.         qg(i,j)=qcg(i,j,k)
    1696. 

  1696. !        IF (QV(I,J)+QB0 .LE. 0.) QV(I,J)=-QB0
    1697. 

  1697.          if (qc(i,j) .le.  cmin1) qc(i,j)=0.0
    1698. 

  1698.          if (qr(i,j) .le.  cmin1) qr(i,j)=0.0
    1699. 

  1699.          if (qi(i,j) .le.  cmin1) qi(i,j)=0.0
    1700. 

  1700.          if (qs(i,j) .le.  cmin1) qs(i,j)=0.0
    1701. 

  1701.          if (qg(i,j) .le.  cmin1) qg(i,j)=0.0
    1702. 

  1702.         tair(i,j)=(pt(i,j)+tb0)*pi0
    1703. 

  1703.         tairc(i,j)=tair(i,j)-t0
    1704. 

  1704.          zr(i,j)=zrr
    1705. 

  1705.          zs(i,j)=zsr
    1706. 

  1706.          zg(i,j)=zgr
    1707. 

  1707.          vr(i,j)=0.0
    1708. 

  1708.          vs(i,j)=0.0
    1709. 

  1709.          vg(i,j)=0.0
    1710. 

  1710. 

  1711. !JJS 10/7/2008     vvvvv
    1712. 

  1712.     IF (IWARM .EQ. 1) THEN
    1713. 

  1713. !JJS   for calculating processes related to warm rain only
    1714. 

  1714.                 qi(i,j)=0.0
    1715. 

  1715.                 qs(i,j)=0.0
    1716. 

  1716.                 qg(i,j)=0.0
    1717. 

  1717.                 dep(i,j)=0.
    1718. 

  1718.                 pint(i,j)=0.
    1719. 

  1719.                 psdep(i,j)=0.
    1720. 

  1720.                 pgdep(i,j)=0.
    1721. 

  1721.                 dd1(i,j)=0.
    1722. 

  1722.                 pgsub(i,j)=0.
    1723. 

  1723.                 psmlt(i,j)=0.
    1724. 

  1724.                 pgmlt(i,j)=0.
    1725. 

  1725.                 pimlt(i,j)=0.
    1726. 

  1726.                 psacw(i,j)=0.
    1727. 

  1727.                 piacr(i,j)=0.
    1728. 

  1728.                 psfw(i,j)=0.
    1729. 

  1729.                 pgfr(i,j)=0.
    1730. 

  1730.                 dgacw(i,j)=0.
    1731. 

  1731.                 dgacr(i,j)=0.
    1732. 

  1732.                 psacr(i,j)=0.
    1733. 

  1733.                 wgacr(i,j)=0.
    1734. 

  1734.                 pihom(i,j)=0.
    1735. 

  1735.                 pidw(i,j)=0.
    1736. 

  1736. 

  1737.                 if (qr(i,j) .gt. cmin1) then
    1738. 

  1738.                    dd(i,j)=r00*qr(i,j)
    1739. 

  1739.                    y1(i,j)=dd(i,j)**.25
    1740. 

  1740.                    zr(i,j)=zrc/y1(i,j)
    1741. 

  1741.                    vr(i,j)=max(vrcf*dd(i,j)**bwq, 0.)
    1742. 

  1742.                 endif
    1743. 

  1743. 

  1744. !* 21 * PRAUT   AUTOCONVERSION OF QC TO QR                        **21**
    1745. 

  1745. !* 22 * PRACW : ACCRETION OF QC BY QR                             **22**
    1746. 

  1746.                 pracw(i,j)=0.
    1747. 

  1747.                 praut(i,j)=0.0
    1748. 

  1748.                 pracw(i,j)=r22f*qc(i,j)/zr(i,j)**bw3
    1749. 

  1749.                 y1(i,j)=qc(i,j)-bnd3
    1750. 

  1750.                 if (y1(i,j).gt.0.0) then
    1751. 

  1751.                     praut(i,j)=r00*y1(i,j)*y1(i,j)/(1.2e-4+rn21/y1(i,j))
    1752. 

  1752.                  endif
    1753. 

  1753. 

  1754. !C********   HANDLING THE NEGATIVE CLOUD WATER (QC)    ******************
    1755. 

  1755.                  Y1(I,J)=QC(I,J)/D2T
    1756. 

  1756.                  PRAUT(I,J)=MIN(Y1(I,J), PRAUT(I,J))
    1757. 

  1757.                  PRACW(I,J)=MIN(Y1(I,J), PRACW(I,J))
    1758. 

  1758.                  Y1(I,J)=(PRAUT(I,J)+PRACW(I,J))*D2T
    1759. 

  1759.                
    1760. 

  1760.                if (qc(i,j) .lt. y1(i,j) .and. y1(i,j) .ge. cmin2) then
    1761. 

  1761.                    y2(i,j)=qc(i,j)/(y1(i,j)+cmin2)
    1762. 

  1762.                    praut(i,j)=praut(i,j)*y2(i,j)
    1763. 

  1763.                    pracw(i,j)=pracw(i,j)*y2(i,j)
    1764. 

  1764.                    qc(i,j)=0.0
    1765. 

  1765.                else
    1766. 

  1766.                   qc(i,j)=qc(i,j)-y1(i,j)
    1767. 

  1767.                endif
    1768. 

  1768.                
    1769. 

  1769.                PR(I,J)=(PRAUT(I,J)+PRACW(I,J))*D2T
    1770. 

  1770.                QR(I,J)=QR(I,J)+PR(I,J)
    1771. 

  1771.                         
    1772. 

  1772. !*****   TAO ET AL (1989) SATURATION TECHNIQUE  ***********************
    1773. 

  1773.            
    1774. 

  1774.            cnd(i,j)=0.0
    1775. 

  1775.            tair(i,j)=(pt(i,j)+tb0)*pi0
    1776. 

  1776.               y1(i,j)=1./(tair(i,j)-c358)
    1777. 

  1777.               qsw(i,j)=rp0*exp(c172-c409*y1(i,j))
    1778. 

  1778.               dd(i,j)=cp409*y1(i,j)*y1(i,j)
    1779. 

  1779.               dm(i,j)=qv(i,j)+qb0-qsw(i,j)
    1780. 

  1780.               cnd(i,j)=dm(i,j)/(1.+avcp*dd(i,j)*qsw(i,j))
    1781. 

  1781. !c    ******   condensation or evaporation of qc  ******
    1782. 

  1782.               cnd(i,j)=max(-qc(i,j), cnd(i,j))
    1783. 

  1783.                          pt(i,j)=pt(i,j)+avcp*cnd(i,j)
    1784. 

  1784.              qv(i,j)=qv(i,j)-cnd(i,j)
    1785. 

  1785.                          qc(i,j)=qc(i,j)+cnd(i,j)
    1786. 

  1786. 

  1787. !C     ******   EVAPORATION   ******
    1788. 

  1788. !* 23 * ERN : EVAPORATION OF QR (SUBSATURATION)                   **23**
    1789. 

  1789.             ern(i,j)=0.0
    1790. 

  1790. 

  1791.             if(qr(i,j).gt.0.0) then
    1792. 

  1792.                tair(i,j)=(pt(i,j)+tb0)*pi0
    1793. 

  1793.                rtair(i,j)=1./(tair(i,j)-c358)
    1794. 

  1794.                qsw(i,j)=rp0*exp(c172-c409*rtair(i,j))
    1795. 

  1795.                ssw(i,j)=(qv(i,j)+qb0)/qsw(i,j)-1.0
    1796. 

  1796.                dm(i,j)=qv(i,j)+qb0-qsw(i,j)
    1797. 

  1797.                rsub1(i,j)=cv409*qsw(i,j)*rtair(i,j)*rtair(i,j)
    1798. 

  1798.                dd1(i,j)=max(-dm(i,j)/(1.+rsub1(i,j)), 0.0)
    1799. 

  1799.                y1(i,j)=.78/zr(i,j)**2+r23af*scv(i,j)/zr(i,j)**bwh5
    1800. 

  1800.                y2(i,j)=r23br/(tca(i,j)*tair(i,j)**2)+1./(dwv(i,j) &
    1801. 

  1801.                        *qsw(i,j))
    1802. 

  1802. !cccc
    1803. 

  1803.                ern(i,j)=r23t*ssw(i,j)*y1(i,j)/y2(i,j)
    1804. 

  1804.                ern(i,j)=min(dd1(i,j),qr(i,j),max(ern(i,j),0.))
    1805. 

  1805.                pt(i,j)=pt(i,j)-avcp*ern(i,j)
    1806. 

  1806.                qv(i,j)=qv(i,j)+ern(i,j)
    1807. 

  1807.                qr(i,j)=qr(i,j)-ern(i,j)
    1808. 

  1808.             endif
    1809. 

  1809. 

  1810.        ELSE       ! part of if (iwarm.eq.1) then
    1811. 

  1811. !JJS 10/7/2008     ^^^^^
    1812. 

  1812. 

  1813. !JJS   for calculating processes related to both ice and warm rain
    1814. 

  1814. 

  1815. !     ***   COMPUTE ZR,ZS,ZG,VR,VS,VG      *****************************
    1816. 

  1816. 

  1817.             if (qr(i,j) .gt. cmin1) then
    1818. 

  1818.                dd(i,j)=r00*qr(i,j)
    1819. 

  1819.                y1(i,j)=dd(i,j)**.25
    1820. 

  1820.                zr(i,j)=zrc/y1(i,j)
    1821. 

  1821.                vr(i,j)=max(vrcf*dd(i,j)**bwq, 0.)
    1822. 

  1822.             endif
    1823. 

  1823. 

  1824.             if (qs(i,j) .gt. cmin1) then
    1825. 

  1825.                dd(i,j)=r00*qs(i,j)
    1826. 

  1826.                y1(i,j)=dd(i,j)**.25
    1827. 

  1827.                zs(i,j)=zsc/y1(i,j)
    1828. 

  1828.                vs(i,j)=max(vscf*dd(i,j)**bsq, 0.)
    1829. 

  1829.             endif
    1830. 

  1830. 

  1831.             if (qg(i,j) .gt. cmin1) then
    1832. 

  1832.                dd(i,j)=r00*qg(i,j)
    1833. 

  1833.                y1(i,j)=dd(i,j)**.25
    1834. 

  1834.                zg(i,j)=zgc/y1(i,j)
    1835. 

  1835.                if(ihail .eq. 1) then
    1836. 

  1836.                   vg(i,j)=max(vgcr*dd(i,j)**bgq, 0.)
    1837. 

  1837.                else
    1838. 

  1838.                   vg(i,j)=max(vgcf*dd(i,j)**bgq, 0.)
    1839. 

  1839.                endif
    1840. 

  1840.             endif
    1841. 

  1841. 

  1842.             if (qr(i,j) .le. cmin2) vr(i,j)=0.0
    1843. 

  1843.             if (qs(i,j) .le. cmin2) vs(i,j)=0.0
    1844. 

  1844.             if (qg(i,j) .le. cmin2) vg(i,j)=0.0
    1845. 

  1845. 

  1846. !     ******************************************************************
    1847. 

  1847. !     ***   Y1 : DYNAMIC VISCOSITY OF AIR (U)
    1848. 

  1848. !     ***   DWV : DIFFUSIVITY OF WATER VAPOR IN AIR (PI)
    1849. 

  1849. !     ***   TCA : THERMAL CONDUCTIVITY OF AIR (KA)
    1850. 

  1850. !     ***   Y2 : KINETIC VISCOSITY (V)
    1851. 

  1851. 

  1852.             y1(i,j)=c149*tair(i,j)**1.5/(tair(i,j)+120.)
    1853. 

  1853.             dwv(i,j)=dwvp*tair(i,j)**1.81
    1854. 

  1854.             tca(i,j)=c141*y1(i,j)
    1855. 

  1855.             scv(i,j)=1./((rr0*y1(i,j))**.1666667*dwv(i,j)**.3333333)
    1856. 

  1856. !JJS  100    CONTINUE
    1857. 

  1857. 

  1858. !*  1 * PSAUT : AUTOCONVERSION OF QI TO QS                        ***1**
    1859. 

  1859. !*  3 * PSACI : ACCRETION OF QI TO QS                             ***3**
    1860. 

  1860. !*  4 * PSACW : ACCRETION OF QC BY QS (RIMING) (QSACW FOR PSMLT)  ***4**
    1861. 

  1861. !*  5 * PRACI : ACCRETION OF QI BY QR                             ***5**
    1862. 

  1862. !*  6 * PIACR : ACCRETION OF QR OR QG BY QI                       ***6**
    1863. 

  1863. 

  1864. !JJS         DO 125 J=3,JLES
    1865. 

  1865. !JJS         DO 125 I=3,ILES
    1866. 

  1866.             psaut(i,j)=0.0
    1867. 

  1867.             psaci(i,j)=0.0
    1868. 

  1868.             praci(i,j)=0.0
    1869. 

  1869.             piacr(i,j)=0.0
    1870. 

  1870.             psacw(i,j)=0.0
    1871. 

  1871.             qsacw(i,j)=0.0
    1872. 

  1872.             dd(i,j)=1./zs(i,j)**bs3
    1873. 

  1873. 

  1874.             if (tair(i,j).lt.t0) then
    1875. 

  1875.                esi(i,j)=exp(.025*tairc(i,j))
    1876. 

  1876.                psaut(i,j)=r2is*max(rn1*esi(i,j)*(qi(i,j)-bnd1) ,0.0)
    1877. 

  1877.                psaci(i,j)=r2is*r3f*esi(i,j)*qi(i,j)*dd(i,j)
    1878. 

  1878. !JJS 3/30/06
    1879. 

  1879. !    to cut water to snow accretion by half
    1880. 

  1880. !               PSACW(I,J)=R4F*QC(I,J)*DD(I,J)
    1881. 

  1881.                psacw(i,j)=r2is*0.5*r4f*qc(i,j)*dd(i,j)
    1882. 

  1882. !JJS 3/30/06
    1883. 

  1883.                praci(i,j)=r2is*r5f*qi(i,j)/zr(i,j)**bw3
    1884. 

  1884.                piacr(i,j)=r2is*r6f*qi(i,j)*(zr(i,j)**(-bw6))
    1885. 

  1885. !JJS               PIACR(I,J)=R6F*QI(I,J)/ZR(I,J)**BW6
    1886. 

  1886.             else
    1887. 

  1887.                qsacw(i,j)=r2is*r4f*qc(i,j)*dd(i,j)
    1888. 

  1888.             endif
    1889. 

  1889. 

  1890. !* 21 * PRAUT   AUTOCONVERSION OF QC TO QR                        **21**
    1891. 

  1891. !* 22 * PRACW : ACCRETION OF QC BY QR                             **22**
    1892. 

  1892. 

  1893.             pracw(i,j)=r22f*qc(i,j)/zr(i,j)**bw3
    1894. 

  1894.             praut(i,j)=0.0
    1895. 

  1895.             y1(i,j)=qc(i,j)-bnd3
    1896. 

  1896.             if (y1(i,j).gt.0.0) then
    1897. 

  1897.                praut(i,j)=r00*y1(i,j)*y1(i,j)/(1.2e-4+rn21/y1(i,j))
    1898. 

  1898.             endif
    1899. 

  1899. 

  1900. !* 12 * PSFW : BERGERON PROCESSES FOR QS (KOENING, 1971)          **12**
    1901. 

  1901. !* 13 * PSFI : BERGERON PROCESSES FOR QS                          **13**
    1902. 

  1902. 

  1903.             psfw(i,j)=0.0
    1904. 

  1904.             psfi(i,j)=0.0
    1905. 

  1905.             pidep(i,j)=0.0
    1906. 

  1906. 

  1907.             if(tair(i,j).lt.t0.and.qi(i,j).gt.cmin) then
    1908. 

  1908.                y1(i,j)=max( min(tairc(i,j), -1.), -31.)
    1909. 

  1909.                it(i,j)=int(abs(y1(i,j)))
    1910. 

  1910.                y1(i,j)=rn12a(it(i,j))
    1911. 

  1911.                y2(i,j)=rn12b(it(i,j))
    1912. 

  1912. !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    1913. 

  1913.           psfw(i,j)=r2is* &
    1914. 

  1914.                     max(d2t*y1(i,j)*(y2(i,j)+r12r*qc(i,j))*qi(i,j),0.0)
    1915. 

  1915.                rtair(i,j)=1./(tair(i,j)-c76)
    1916. 

  1916.                y2(i,j)=exp(c218-c580*rtair(i,j))
    1917. 

  1917.                qsi(i,j)=rp0*y2(i,j)
    1918. 

  1918.                esi(i,j)=c610*y2(i,j)
    1919. 

  1919.                ssi(i,j)=(qv(i,j)+qb0)/qsi(i,j)-1.
    1920. 

  1920.                r_nci=min(1.e-6*exp(-.46*tairc(i,j)),1.)
    1921. 

  1921. !              R_NCI=min(1.e-8*EXP(-.6*TAIRC(I,J)),1.) ! use Tao's
    1922. 

  1922.                dm(i,j)=max( (qv(i,j)+qb0-qsi(i,j)), 0.)
    1923. 

  1923.                rsub1(i,j)=cs580*qsi(i,j)*rtair(i,j)*rtair(i,j)
    1924. 

  1924.                y3(i,j)=1./tair(i,j)
    1925. 

  1925.           dd(i,j)=y3(i,j)*(rn30a*y3(i,j)-rn30b)+rn30c*tair(i,j)/esi(i,j)
    1926. 

  1926.                y1(i,j)=206.18*ssi(i,j)/dd(i,j)
    1927. 

  1927.                pidep(i,j)=y1(i,j)*sqrt(r_nci*qi(i,j)/r00)
    1928. 

  1928.                dep(i,j)=dm(i,j)/(1.+rsub1(i,j))/d2t
    1929. 

  1929.                if(dm(i,j).gt.cmin2) then
    1930. 

  1930.                   a2=1.
    1931. 

  1931.                 if(pidep(i,j).gt.dep(i,j).and.pidep(i,j).gt.cmin2) then
    1932. 

  1932.                      a2=dep(i,j)/pidep(i,j)
    1933. 

  1933.                      pidep(i,j)=dep(i,j)
    1934. 

  1934.                 endif
    1935. 

  1935.                   psfi(i,j)=r2is*a2*.5*qi(i,j)*y1(i,j)/(sqrt(ami100) &
    1936. 

  1936.                           -sqrt(ami40))
    1937. 

  1937.                   elseif(dm(i,j).lt.-cmin2) then
    1938. 

  1938. !
    1939. 

  1939. !        SUBLIMATION TERMS USED ONLY WHEN SATURATION ADJUSTMENT FOR ICE
    1940. 

  1940. !        IS TURNED OFF
    1941. 

  1941. !
    1942. 

  1942.                   pidep(i,j)=0.
    1943. 

  1943.                   psfi(i,j)=0.
    1944. 

  1944.                else
    1945. 

  1945.                   pidep(i,j)=0.
    1946. 

  1946.                   psfi(i,j)=0.
    1947. 

  1947.                endif
    1948. 

  1948.             endif
    1949. 

  1949. 

  1950. !TTT***** QG=QG+MIN(PGDRY,PGWET)
    1951. 

  1951. !*  9 * PGACS : ACCRETION OF QS BY QG (DGACS,WGACS: DRY AND WET)  ***9**
    1952. 

  1952. !* 14 * DGACW : ACCRETION OF QC BY QG (QGACW FOR PGMLT)           **14**
    1953. 

  1953. !* 16 * DGACR : ACCRETION OF QR TO QG (QGACR FOR PGMLT)           **16**
    1954. 

  1954. 

  1955.             if(qc(i,j)+qr(i,j).lt.1.e-4) then
    1956. 

  1956.                ee1=.01
    1957. 

  1957.               else
    1958. 

  1958.                  ee1=1.
    1959. 

  1959.               endif
    1960. 

  1960.             ee2=0.09
    1961. 

  1961.             egs(i,j)=ee1*exp(ee2*tairc(i,j))
    1962. 

  1962. !            EGS(I,J)=0.1 ! 6/15/02 tao's
    1963. 

  1963.             if (tair(i,j).ge.t0) egs(i,j)=1.0
    1964. 

  1964.             y1(i,j)=abs(vg(i,j)-vs(i,j))
    1965. 

  1965.             y2(i,j)=zs(i,j)*zg(i,j)
    1966. 

  1966.             y3(i,j)=5./y2(i,j)
    1967. 

  1967.             y4(i,j)=.08*y3(i,j)*y3(i,j)
    1968. 

  1968.             y5(i,j)=.05*y3(i,j)*y4(i,j)
    1969. 

  1969.             dd(i,j)=y1(i,j)*(y3(i,j)/zs(i,j)**5+y4(i,j)/zs(i,j)**3 &
    1970. 

  1970.                     +y5(i,j)/zs(i,j))
    1971. 

  1971.             pgacs(i,j)=r2ig*r2is*r9r*egs(i,j)*dd(i,j)
    1972. 

  1972. !JJS 1/3/06 from Steve and Chunglin
    1973. 

  1973.             if (ihail.eq.1) then
    1974. 

  1974.                dgacs(i,j)=pgacs(i,j)
    1975. 

  1975.             else
    1976. 

  1976.                dgacs(i,j)=0.
    1977. 

  1977.             endif
    1978. 

  1978. !JJS 1/3/06 from Steve and Chunglin
    1979. 

  1979.             wgacs(i,j)=r2ig*r2is*r9r*dd(i,j)
    1980. 

  1980. !            WGACS(I,J)=0.  ! 6/15/02 tao's
    1981. 

  1981.             y1(i,j)=1./zg(i,j)**bg3
    1982. 

  1982. 

  1983.             if(ihail .eq. 1) then
    1984. 

  1984.                dgacw(i,j)=r2ig*max(r14r*qc(i,j)*y1(i,j), 0.0)
    1985. 

  1985.             else
    1986. 

  1986.                dgacw(i,j)=r2ig*max(r14f*qc(i,j)*y1(i,j), 0.0)
    1987. 

  1987.             endif
    1988. 

  1988. 

  1989.             qgacw(i,j)=dgacw(i,j)
    1990. 

  1990.             y1(i,j)=abs(vg(i,j)-vr(i,j))
    1991. 

  1991.             y2(i,j)=zr(i,j)*zg(i,j)
    1992. 

  1992.             y3(i,j)=5./y2(i,j)
    1993. 

  1993.             y4(i,j)=.08*y3(i,j)*y3(i,j)
    1994. 

  1994.             y5(i,j)=.05*y3(i,j)*y4(i,j)
    1995. 

  1995.             dd(i,j)=r16r*y1(i,j)*(y3(i,j)/zr(i,j)**5+y4(i,j)/zr(i,j)**3 &
    1996. 

  1996.                     +y5(i,j)/zr(i,j))
    1997. 

  1997.             dgacr(i,j)=r2ig*max(dd(i,j), 0.0)
    1998. 

  1998.             qgacr(i,j)=dgacr(i,j)
    1999. 

  1999. 

  2000.             if (tair(i,j).ge.t0) then
    2001. 

  2001.                dgacs(i,j)=0.0
    2002. 

  2002.                wgacs(i,j)=0.0
    2003. 

  2003.                dgacw(i,j)=0.0
    2004. 

  2004.                dgacr(i,j)=0.0
    2005. 

  2005.             else
    2006. 

  2006.                pgacs(i,j)=0.0
    2007. 

  2007.                qgacw(i,j)=0.0
    2008. 

  2008.                qgacr(i,j)=0.0
    2009. 

  2009.             endif
    2010. 

  2010. 

  2011. !*******PGDRY : DGACW+DGACI+DGACR+DGACS                           ******
    2012. 

  2012. !* 15 * DGACI : ACCRETION OF QI BY QG (WGACI FOR WET GROWTH)      **15**
    2013. 

  2013. !* 17 * PGWET : WET GROWTH OF QG                                  **17**
    2014. 

  2014. 

  2015.             dgaci(i,j)=0.0
    2016. 

  2016.             wgaci(i,j)=0.0
    2017. 

  2017.             pgwet(i,j)=0.0
    2018. 

  2018. 

  2019.             if (tair(i,j).lt.t0) then
    2020. 

  2020.                y1(i,j)=qi(i,j)/zg(i,j)**bg3
    2021. 

  2021.                if (ihail.eq.1) then
    2022. 

  2022.                   dgaci(i,j)=r2ig*r15r*y1(i,j)
    2023. 

  2023.                   wgaci(i,j)=r2ig*r15ar*y1(i,j)
    2024. 

  2024. !                  WGACI(I,J)=0.  ! 6/15/02 tao's
    2025. 

  2025.                else
    2026. 

  2026. 

  2027. !JJS                  DGACI(I,J)=r2ig*R15F*Y1(I,J)
    2028. 

  2028.                    dgaci(i,j)=0.
    2029. 

  2029.                   wgaci(i,j)=r2ig*r15af*y1(i,j)
    2030. 

  2030. !                  WGACI(I,J)=0.  ! 6/15/02 tao's
    2031. 

  2031.                endif
    2032. 

  2032. !
    2033. 

  2033.                if (tairc(i,j).ge.-50.) then
    2034. 

  2034.                 if (alf+rn17c*tairc(i,j) .eq. 0.) then
    2035. 

  2035.                    write(91,*) itimestep, i,j,k, alf, rn17c, tairc(i,j)
    2036. 

  2036.                 endif
    2037. 

  2037.                 y1(i,j)=1./(alf+rn17c*tairc(i,j))
    2038. 

  2038.                 if (ihail.eq.1) then
    2039. 

  2039.                    y3(i,j)=.78/zg(i,j)**2+r17aq*scv(i,j)/zg(i,j)**bgh5
    2040. 

  2040.                 else
    2041. 

  2041.                    y3(i,j)=.78/zg(i,j)**2+r17as*scv(i,j)/zg(i,j)**bgh5
    2042. 

  2042.                 endif
    2043. 

  2043.                 y4(i,j)=alvr*dwv(i,j)*(rp0-(qv(i,j)+qb0)) &
    2044. 

  2044.                         -tca(i,j)*tairc(i,j)
    2045. 

  2045.                 dd(i,j)=y1(i,j)*(r17r*y4(i,j)*y3(i,j) &
    2046. 

  2046.                        +(wgaci(i,j)+wgacs(i,j))*(alf+rn17b*tairc(i,j)))
    2047. 

  2047.                 pgwet(i,j)=r2ig*max(dd(i,j), 0.0)
    2048. 

  2048.                endif
    2049. 

  2049.             endif
    2050. 

  2050. !JJS  125    CONTINUE
    2051. 

  2051. 

  2052. !********   HANDLING THE NEGATIVE CLOUD WATER (QC)    ******************
    2053. 

  2053. !********   HANDLING THE NEGATIVE CLOUD ICE (QI)      ******************
    2054. 

  2054. 

  2055. !JJS         DO 150 J=3,JLES
    2056. 

  2056. !JJS         DO 150 I=3,ILES
    2057. 

  2057. 

  2058.             y1(i,j)=qc(i,j)/d2t
    2059. 

  2059.             psacw(i,j)=min(y1(i,j), psacw(i,j))
    2060. 

  2060.             praut(i,j)=min(y1(i,j), praut(i,j))
    2061. 

  2061.             pracw(i,j)=min(y1(i,j), pracw(i,j))
    2062. 

  2062.             psfw(i,j)= min(y1(i,j), psfw(i,j))
    2063. 

  2063.             dgacw(i,j)=min(y1(i,j), dgacw(i,j))
    2064. 

  2064.             qsacw(i,j)=min(y1(i,j), qsacw(i,j))
    2065. 

  2065.             qgacw(i,j)=min(y1(i,j), qgacw(i,j))
    2066. 

  2066. 

  2067.             y1(i,j)=(psacw(i,j)+praut(i,j)+pracw(i,j)+psfw(i,j) &
    2068. 

  2068.                     +dgacw(i,j)+qsacw(i,j)+qgacw(i,j))*d2t
    2069. 

  2069.             qc(i,j)=qc(i,j)-y1(i,j)
    2070. 

  2070. 

  2071.             if (qc(i,j) .lt. 0.0) then
    2072. 

  2072.                a1=1.
    2073. 

  2073.                if (y1(i,j) .ne. 0.0) a1=qc(i,j)/y1(i,j)+1.
    2074. 

  2074.                psacw(i,j)=psacw(i,j)*a1
    2075. 

  2075.                praut(i,j)=praut(i,j)*a1
    2076. 

  2076.                pracw(i,j)=pracw(i,j)*a1
    2077. 

  2077.                psfw(i,j)=psfw(i,j)*a1
    2078. 

  2078.                dgacw(i,j)=dgacw(i,j)*a1
    2079. 

  2079.                qsacw(i,j)=qsacw(i,j)*a1
    2080. 

  2080.                qgacw(i,j)=qgacw(i,j)*a1
    2081. 

  2081.                qc(i,j)=0.0
    2082. 

  2082.             endif
    2083. 

  2083. !c
    2084. 

  2084. !
    2085. 

  2085. !******** SHED PROCESS (WGACR=PGWET-DGACW-WGACI-WGACS)
    2086. 

  2086. !c
    2087. 

  2087.             wgacr(i,j)=pgwet(i,j)-dgacw(i,j)-wgaci(i,j)-wgacs(i,j)
    2088. 

  2088.             y2(i,j)=dgacw(i,j)+dgaci(i,j)+dgacr(i,j)+dgacs(i,j)
    2089. 

  2089.             if (pgwet(i,j).ge.y2(i,j)) then
    2090. 

  2090.                wgacr(i,j)=0.0
    2091. 

  2091.                wgaci(i,j)=0.0
    2092. 

  2092.                wgacs(i,j)=0.0
    2093. 

  2093.             else
    2094. 

  2094.                dgacr(i,j)=0.0
    2095. 

  2095.                dgaci(i,j)=0.0
    2096. 

  2096.                dgacs(i,j)=0.0
    2097. 

  2097.             endif
    2098. 

  2098. !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    2099. 

  2099. !c
    2100. 

  2100.             y1(i,j)=qi(i,j)/d2t
    2101. 

  2101.             psaut(i,j)=min(y1(i,j), psaut(i,j))
    2102. 

  2102.             psaci(i,j)=min(y1(i,j), psaci(i,j))
    2103. 

  2103.             praci(i,j)=min(y1(i,j), praci(i,j))
    2104. 

  2104.             psfi(i,j)= min(y1(i,j), psfi(i,j))
    2105. 

  2105.             dgaci(i,j)=min(y1(i,j), dgaci(i,j))
    2106. 

  2106.             wgaci(i,j)=min(y1(i,j), wgaci(i,j))
    2107. 

  2107. !
    2108. 

  2108.             y2(i,j)=(psaut(i,j)+psaci(i,j)+praci(i,j)+psfi(i,j) &
    2109. 

  2109.                    +dgaci(i,j)+wgaci(i,j))*d2t
    2110. 

  2110.             qi(i,j)=qi(i,j)-y2(i,j)+pidep(i,j)*d2t
    2111. 

  2111. 

  2112.             if (qi(i,j).lt.0.0) then
    2113. 

  2113.                a2=1.
    2114. 

  2114.                if (y2(i,j) .ne. 0.0) a2=qi(i,j)/y2(i,j)+1.
    2115. 

  2115.                psaut(i,j)=psaut(i,j)*a2
    2116. 

  2116.                psaci(i,j)=psaci(i,j)*a2
    2117. 

  2117.                praci(i,j)=praci(i,j)*a2
    2118. 

  2118.                psfi(i,j)=psfi(i,j)*a2
    2119. 

  2119.                dgaci(i,j)=dgaci(i,j)*a2
    2120. 

  2120.                wgaci(i,j)=wgaci(i,j)*a2
    2121. 

  2121.                qi(i,j)=0.0
    2122. 

  2122.             endif
    2123. 

  2123. !
    2124. 

  2124.             dlt3(i,j)=0.0
    2125. 

  2125.             dlt2(i,j)=0.0
    2126. 

  2126. !
    2127. 

  2127. 

  2128. !            DLT4(I,J)=1.0
    2129. 

  2129. !            if(qc(i,j) .gt. 5.e-4) dlt4(i,j)=0.0
    2130. 

  2130. !            if(qs(i,j) .le. 1.e-4) dlt4(i,j)=1.0
    2131. 

  2131. !
    2132. 

  2132. !            IF (TAIR(I,J).ge.T0) THEN
    2133. 

  2133. !               DLT4(I,J)=0.0
    2134. 

  2134. !            ENDIF
    2135. 

  2135. 

  2136.             if (tair(i,j).lt.t0) then
    2137. 

  2137.                if (qr(i,j).lt.1.e-4) then
    2138. 

  2138.                   dlt3(i,j)=1.0
    2139. 

  2139.                   dlt2(i,j)=1.0
    2140. 

  2140.                endif
    2141. 

  2141.                if (qs(i,j).ge.1.e-4) then
    2142. 

  2142.                   dlt2(i,j)=0.0
    2143. 

  2143.                endif
    2144. 

  2144.             endif
    2145. 

  2145. 

  2146.             if (ice2 .eq. 1) then
    2147. 

  2147.                   dlt3(i,j)=1.0
    2148. 

  2148.                   dlt2(i,j)=1.0
    2149. 

  2149.             endif
    2150. 

  2150. !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    2151. 

  2151.             pr(i,j)=(qsacw(i,j)+praut(i,j)+pracw(i,j)+qgacw(i,j))*d2t
    2152. 

  2152.             ps(i,j)=(psaut(i,j)+psaci(i,j)+psacw(i,j)+psfw(i,j) &
    2153. 

  2153.                     +psfi(i,j)+dlt3(i,j)*praci(i,j))*d2t
    2154. 

  2154. !           PS(I,J)=(PSAUT(I,J)+PSACI(I,J)+dlt4(i,j)*PSACW(I,J)
    2155. 

  2155. !    1              +PSFW(I,J)+PSFI(I,J)+DLT3(I,J)*PRACI(I,J))*D2T
    2156. 

  2156.             pg(i,j)=((1.-dlt3(i,j))*praci(i,j)+dgaci(i,j)+wgaci(i,j) &
    2157. 

  2157.                     +dgacw(i,j))*d2t
    2158. 

  2158. !           PG(I,J)=((1.-DLT3(I,J))*PRACI(I,J)+DGACI(I,J)+WGACI(I,J)
    2159. 

  2159. !    1              +DGACW(I,J)+(1.-dlt4(i,j))*PSACW(I,J))*D2T
    2160. 

  2160. 

  2161. !JJS  150    CONTINUE
    2162. 

  2162. 

  2163. !*  7 * PRACS : ACCRETION OF QS BY QR                             ***7**
    2164. 

  2164. !*  8 * PSACR : ACCRETION OF QR BY QS (QSACR FOR PSMLT)           ***8**
    2165. 

  2165. 

  2166. !JJS         DO 175 J=3,JLES
    2167. 

  2167. !JJS         DO 175 I=3,ILES
    2168. 

  2168. 

  2169.             y1(i,j)=abs(vr(i,j)-vs(i,j))
    2170. 

  2170.             y2(i,j)=zr(i,j)*zs(i,j)
    2171. 

  2171.             y3(i,j)=5./y2(i,j)
    2172. 

  2172.             y4(i,j)=.08*y3(i,j)*y3(i,j)
    2173. 

  2173.             y5(i,j)=.05*y3(i,j)*y4(i,j)
    2174. 

  2174.             pracs(i,j)=r2ig*r2is*r7r*y1(i,j)*(y3(i,j)/zs(i,j)**5 &
    2175. 

  2175.                       +y4(i,j)/zs(i,j)**3+y5(i,j)/zs(i,j))
    2176. 

  2176.             psacr(i,j)=r2is*r8r*y1(i,j)*(y3(i,j)/zr(i,j)**5 &
    2177. 

  2177.                       +y4(i,j)/zr(i,j)**3+y5(i,j)/zr(i,j))
    2178. 

  2178.             qsacr(i,j)=psacr(i,j)
    2179. 

  2179. 

  2180.             if (tair(i,j).ge.t0) then
    2181. 

  2181.                pracs(i,j)=0.0
    2182. 

  2182.                psacr(i,j)=0.0
    2183. 

  2183.             else
    2184. 

  2184.                qsacr(i,j)=0.0
    2185. 

  2185.             endif
    2186. 

  2186. !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    2187. 

  2187. !*  2 * PGAUT : AUTOCONVERSION OF QS TO QG                        ***2**
    2188. 

  2188. !* 18 * PGFR : FREEZING OF QR TO QG                               **18**
    2189. 

  2189. 

  2190.             pgaut(i,j)=0.0
    2191. 

  2191.             pgfr(i,j)=0.0
    2192. 

  2192. 

  2193.             if (tair(i,j) .lt. t0) then
    2194. 

  2194. !               Y1(I,J)=EXP(.09*TAIRC(I,J))
    2195. 

  2195. !               PGAUT(I,J)=r2is*max(RN2*Y1(I,J)*(QS(I,J)-BND2), 0.0)
    2196. 

  2196. !         IF(IHAIL.EQ.1) PGAUT(I,J)=max(RN2*Y1(I,J)*(QS(I,J)-BND2),0.0)
    2197. 

  2197.                y2(i,j)=exp(rn18a*(t0-tair(i,j)))
    2198. 

  2198. !JJS              PGFR(I,J)=r2ig*max(R18R*(Y2(I,J)-1.)/ZR(I,J)**7., 0.0)
    2199. 

  2199. !              pgfr(i,j)=r2ice*max(r18r*(y2(i,j)-1.)* &
    2200. 

  2200. !                                    (zr(i,j)**(-7.)), 0.0)
    2201. 

  2201. !        modify to prevent underflow on some computers (JD)
    2202. 

  2202.                temp = 1./zr(i,j)
    2203. 

  2203.                temp = temp*temp*temp*temp*temp*temp*temp
    2204. 

  2204.                pgfr(i,j)=r2ig*max(r18r*(y2(i,j)-1.)* &
    2205. 

  2205.                                     temp, 0.0)
    2206. 

  2206.             endif
    2207. 

  2207. 

  2208. !JJS  175    CONTINUE
    2209. 

  2209. 

  2210. !********   HANDLING THE NEGATIVE RAIN WATER (QR)    *******************
    2211. 

  2211. !********   HANDLING THE NEGATIVE SNOW (QS)          *******************
    2212. 

  2212. 

  2213. !JJS         DO 200 J=3,JLES
    2214. 

  2214. !JJS         DO 200 I=3,ILES
    2215. 

  2215. 

  2216.             y1(i,j)=qr(i,j)/d2t
    2217. 

  2217.             y2(i,j)=-qg(i,j)/d2t
    2218. 

  2218.             piacr(i,j)=min(y1(i,j), piacr(i,j))
    2219. 

  2219.             dgacr(i,j)=min(y1(i,j), dgacr(i,j))
    2220. 

  2220.             wgacr(i,j)=min(y1(i,j), wgacr(i,j))
    2221. 

  2221.             wgacr(i,j)=max(y2(i,j), wgacr(i,j))
    2222. 

  2222.             psacr(i,j)=min(y1(i,j), psacr(i,j))
    2223. 

  2223.             pgfr(i,j)= min(y1(i,j), pgfr(i,j))
    2224. 

  2224.             del=0.
    2225. 

  2225.             if(wgacr(i,j) .lt. 0.) del=1.
    2226. 

  2226.             y1(i,j)=(piacr(i,j)+dgacr(i,j)+(1.-del)*wgacr(i,j) &
    2227. 

  2227.                     +psacr(i,j)+pgfr(i,j))*d2t
    2228. 

  2228.             qr(i,j)=qr(i,j)+pr(i,j)-y1(i,j)-del*wgacr(i,j)*d2t
    2229. 

  2229.             if (qr(i,j) .lt. 0.0) then
    2230. 

  2230.                a1=1.
    2231. 

  2231.                if(y1(i,j) .ne. 0.) a1=qr(i,j)/y1(i,j)+1.
    2232. 

  2232.                piacr(i,j)=piacr(i,j)*a1
    2233. 

  2233.                dgacr(i,j)=dgacr(i,j)*a1
    2234. 

  2234.                if (wgacr(i,j).gt.0.) wgacr(i,j)=wgacr(i,j)*a1
    2235. 

  2235.                pgfr(i,j)=pgfr(i,j)*a1
    2236. 

  2236.                psacr(i,j)=psacr(i,j)*a1
    2237. 

  2237.                qr(i,j)=0.0
    2238. 

  2238.             endif
    2239. 

  2239. !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    2240. 

  2240.             prn(i,j)=d2t*((1.-dlt3(i,j))*piacr(i,j)+dgacr(i,j) &
    2241. 

  2241.                      +wgacr(i,j)+(1.-dlt2(i,j))*psacr(i,j)+pgfr(i,j))
    2242. 

  2242.             ps(i,j)=ps(i,j)+d2t*(dlt3(i,j)*piacr(i,j) &
    2243. 

  2243.                     +dlt2(i,j)*psacr(i,j))
    2244. 

  2244.             pracs(i,j)=(1.-dlt2(i,j))*pracs(i,j)
    2245. 

  2245.             y1(i,j)=qs(i,j)/d2t
    2246. 

  2246.             pgacs(i,j)=min(y1(i,j), pgacs(i,j))
    2247. 

  2247.             dgacs(i,j)=min(y1(i,j), dgacs(i,j))
    2248. 

  2248.             wgacs(i,j)=min(y1(i,j), wgacs(i,j))
    2249. 

  2249.             pgaut(i,j)=min(y1(i,j), pgaut(i,j))
    2250. 

  2250.             pracs(i,j)=min(y1(i,j), pracs(i,j))
    2251. 

  2251.             psn(i,j)=d2t*(pgacs(i,j)+dgacs(i,j)+wgacs(i,j) &
    2252. 

  2252.                      +pgaut(i,j)+pracs(i,j))
    2253. 

  2253.             qs(i,j)=qs(i,j)+ps(i,j)-psn(i,j)
    2254. 

  2254. 

  2255.             if (qs(i,j).lt.0.0) then
    2256. 

  2256.                a2=1.
    2257. 

  2257.                if (psn(i,j) .ne. 0.0) a2=qs(i,j)/psn(i,j)+1.
    2258. 

  2258.                pgacs(i,j)=pgacs(i,j)*a2
    2259. 

  2259.                dgacs(i,j)=dgacs(i,j)*a2
    2260. 

  2260.                wgacs(i,j)=wgacs(i,j)*a2
    2261. 

  2261.                pgaut(i,j)=pgaut(i,j)*a2
    2262. 

  2262.                pracs(i,j)=pracs(i,j)*a2
    2263. 

  2263.                psn(i,j)=psn(i,j)*a2
    2264. 

  2264.                qs(i,j)=0.0
    2265. 

  2265.             endif
    2266. 

  2266. !
    2267. 

  2267. !C           PSN(I,J)=D2T*(PGACS(I,J)+DGACS(I,J)+WGACS(I,J)
    2268. 

  2268. !c                    +PGAUT(I,J)+PRACS(I,J))
    2269. 

  2269.             y2(i,j)=d2t*(psacw(i,j)+psfw(i,j)+dgacw(i,j)+piacr(i,j) &
    2270. 

  2270.                     +dgacr(i,j)+wgacr(i,j)+psacr(i,j)+pgfr(i,j))
    2271. 

  2271.             pt(i,j)=pt(i,j)+afcp*y2(i,j)
    2272. 

  2272.             qg(i,j)=qg(i,j)+pg(i,j)+prn(i,j)+psn(i,j)
    2273. 

  2273. 

  2274. !JJS  200    CONTINUE
    2275. 

  2275. 

  2276. !* 11 * PSMLT : MELTING OF QS                                     **11**
    2277. 

  2277. !* 19 * PGMLT : MELTING OF QG TO QR                               **19**
    2278. 

  2278. 

  2279. !JJS         DO 225 J=3,JLES
    2280. 

  2280. !JJS         DO 225 I=3,ILES
    2281. 

  2281. 

  2282.             psmlt(i,j)=0.0
    2283. 

  2283.             pgmlt(i,j)=0.0
    2284. 

  2284.             tair(i,j)=(pt(i,j)+tb0)*pi0
    2285. 

  2285. 

  2286.             if (tair(i,j).ge.t0) then
    2287. 

  2287.                tairc(i,j)=tair(i,j)-t0
    2288. 

  2288.                y1(i,j)=tca(i,j)*tairc(i,j)-alvr*dwv(i,j) &
    2289. 

  2289.                                *(rp0-(qv(i,j)+qb0))
    2290. 

  2290.                y2(i,j)=.78/zs(i,j)**2+r101f*scv(i,j)/zs(i,j)**bsh5
    2291. 

  2291.                dd(i,j)=r11rt*y1(i,j)*y2(i,j)+r11at*tairc(i,j) &
    2292. 

  2292.                        *(qsacw(i,j)+qsacr(i,j))
    2293. 

  2293.                psmlt(i,j)=r2is*max(0.0, min(dd(i,j), qs(i,j)))
    2294. 

  2294. 

  2295.                if(ihail.eq.1) then
    2296. 

  2296.                   y3(i,j)=.78/zg(i,j)**2+r19aq*scv(i,j)/zg(i,j)**bgh5
    2297. 

  2297.                else
    2298. 

  2298.                   y3(i,j)=.78/zg(i,j)**2+r19as*scv(i,j)/zg(i,j)**bgh5
    2299. 

  2299.                endif
    2300. 

  2300. 

  2301.                dd1(i,j)=r19rt*y1(i,j)*y3(i,j)+r19bt*tairc(i,j) &
    2302. 

  2302.                         *(qgacw(i,j)+qgacr(i,j))
    2303. 

  2303.                pgmlt(i,j)=r2ig*max(0.0, min(dd1(i,j), qg(i,j)))
    2304. 

  2304.                pt(i,j)=pt(i,j)-afcp*(psmlt(i,j)+pgmlt(i,j))
    2305. 

  2305.                qr(i,j)=qr(i,j)+psmlt(i,j)+pgmlt(i,j)
    2306. 

  2306.                qs(i,j)=qs(i,j)-psmlt(i,j)
    2307. 

  2307.                qg(i,j)=qg(i,j)-pgmlt(i,j)
    2308. 

  2308.             endif
    2309. 

  2309. !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    2310. 

  2310. !* 24 * PIHOM : HOMOGENEOUS FREEZING OF QC TO QI (T < T00)        **24**
    2311. 

  2311. !* 25 * PIDW : DEPOSITION GROWTH OF QC TO QI ( T0 < T <= T00)     **25**
    2312. 

  2312. !* 26 * PIMLT : MELTING OF QI TO QC (T >= T0)                     **26**
    2313. 

  2313. 

  2314.             if (qc(i,j).le.cmin1) qc(i,j)=0.0
    2315. 

  2315.             if (qi(i,j).le.cmin1) qi(i,j)=0.0
    2316. 

  2316.             tair(i,j)=(pt(i,j)+tb0)*pi0
    2317. 

  2317. 

  2318.             if(tair(i,j).le.t00) then
    2319. 

  2319.                pihom(i,j)=qc(i,j)
    2320. 

  2320.             else
    2321. 

  2321.                pihom(i,j)=0.0
    2322. 

  2322.             endif
    2323. 

  2323.             if(tair(i,j).ge.t0) then
    2324. 

  2324.                pimlt(i,j)=qi(i,j)
    2325. 

  2325.             else
    2326. 

  2326.                pimlt(i,j)=0.0
    2327. 

  2327.             endif
    2328. 

  2328.             pidw(i,j)=0.0
    2329. 

  2329. 

  2330.             if (tair(i,j).lt.t0 .and. tair(i,j).gt.t00) then
    2331. 

  2331.                tairc(i,j)=tair(i,j)-t0
    2332. 

  2332.                y1(i,j)=max( min(tairc(i,j), -1.), -31.)
    2333. 

  2333.                it(i,j)=int(abs(y1(i,j)))
    2334. 

  2334.                y2(i,j)=aa1(it(i,j))
    2335. 

  2335.                y3(i,j)=aa2(it(i,j))
    2336. 

  2336.                y4(i,j)=exp(abs(beta*tairc(i,j)))
    2337. 

  2337.                y5(i,j)=(r00*qi(i,j)/(r25a*y4(i,j)))**y3(i,j)
    2338. 

  2338.                pidw(i,j)=min(r25rt*y2(i,j)*y4(i,j)*y5(i,j), qc(i,j))
    2339. 

  2339.             endif
    2340. 

  2340. 

  2341.             y1(i,j)=pihom(i,j)-pimlt(i,j)+pidw(i,j)
    2342. 

  2342.             pt(i,j)=pt(i,j)+afcp*y1(i,j)+ascp*(pidep(i,j))*d2t
    2343. 

  2343.             qv(i,j)=qv(i,j)-(pidep(i,j))*d2t
    2344. 

  2344.             qc(i,j)=qc(i,j)-y1(i,j)
    2345. 

  2345.             qi(i,j)=qi(i,j)+y1(i,j)
    2346. 

  2346. 

  2347. !* 31 * PINT  : INITIATION OF QI                                  **31**
    2348. 

  2348. !* 32 * PIDEP : DEPOSITION OF QI                                  **32**
    2349. 

  2349. !
    2350. 

  2350. !     CALCULATION OF PINT USES DIFFERENT VALUES OF THE INTERCEPT AND SLOPE FOR
    2351. 

  2351. !     THE FLETCHER EQUATION. ALSO, ONLY INITIATE MORE ICE IF THE NEW NUMBER
    2352. 

  2352. !     CONCENTRATION EXCEEDS THAT ALREADY PRESENT.
    2353. 

  2353. !* 31 * pint  : initiation of qi                                  **31**
    2354. 

  2354. !* 32 * pidep : deposition of qi                                  **32**
    2355. 

  2355.            pint(i,j)=0.0
    2356. 

  2356. !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    2357. 

  2357.         if ( itaobraun.eq.1 ) then
    2358. 

  2358.             tair(i,j)=(pt(i,j)+tb0)*pi0
    2359. 

  2359.             if (tair(i,j) .lt. t0) then
    2360. 

  2360. !             if (qi(i,j) .le. cmin) qi(i,j)=0.
    2361. 

  2361.               if (qi(i,j) .le. cmin2) qi(i,j)=0.
    2362. 

  2362.                tairc(i,j)=tair(i,j)-t0
    2363. 

  2363.                rtair(i,j)=1./(tair(i,j)-c76)
    2364. 

  2364.                y2(i,j)=exp(c218-c580*rtair(i,j))
    2365. 

  2365.               qsi(i,j)=rp0*y2(i,j)
    2366. 

  2366.                esi(i,j)=c610*y2(i,j)
    2367. 

  2367.               ssi(i,j)=(qv(i,j)+qb0)/qsi(i,j)-1.
    2368. 

  2368.                         ami50=3.76e-8
    2369. 

  2369. 

  2370. !ccif ( itaobraun.eq.1 ) --> betah=0.5*beta=-.46*0.5=-0.23;   cn0=1.e-6
    2371. 

  2371. !ccif ( itaobraun.eq.0 ) --> betah=0.5*beta=-.6*0.5=-0.30;    cn0=1.e-8
    2372. 

  2372. 

  2373.              y1(i,j)=1./tair(i,j)
    2374. 

  2374. 

  2375. !cc insert a restriction on ice collection that ice collection
    2376. 

  2376. !cc should be stopped at -30 c (with cn0=1.e-6, beta=-.46)
    2377. 

  2377. 

  2378.              tairccri=tairc(i,j)          ! in degree c
    2379. 

  2379.              if(tairccri.le.-30.) tairccri=-30.
    2380. 

  2380. 

  2381. !            y2(i,j)=exp(betah*tairc(i,j))
    2382. 

  2382.              y2(i,j)=exp(betah*tairccri)
    2383. 

  2383.              y3(i,j)=sqrt(qi(i,j))
    2384. 

  2384.              dd(i,j)=y1(i,j)*(rn10a*y1(i,j)-rn10b)+rn10c*tair(i,j) &
    2385. 

  2385.                                                 /esi(i,j)
    2386. 

  2386.           pidep(i,j)=max(r32rt*ssi(i,j)*y2(i,j)*y3(i,j)/dd(i,j), 0.e0)
    2387. 

  2387. 

  2388.            r_nci=min(cn0*exp(beta*tairc(i,j)),1.)
    2389. 

  2389. !          r_nci=min(1.e-6*exp(-.46*tairc(i,j)),1.)
    2390. 

  2390. 

  2391.            dd(i,j)=max(1.e-9*r_nci/r00-qi(i,j)*1.e-9/ami50, 0.)
    2392. 

  2392.                 dm(i,j)=max( (qv(i,j)+qb0-qsi(i,j)), 0.0)
    2393. 

  2393.                 rsub1(i,j)=cs580*qsi(i,j)*rtair(i,j)*rtair(i,j)
    2394. 

  2394.               dep(i,j)=dm(i,j)/(1.+rsub1(i,j))
    2395. 

  2395.               pint(i,j)=max(min(dd(i,j), dm(i,j)), 0.)
    2396. 

  2396. 

  2397. !             pint(i,j)=min(pint(i,j), dep(i,j))
    2398. 

  2398.               pint(i,j)=min(pint(i,j)+pidep(i,j), dep(i,j))
    2399. 

  2399. 

  2400. !              if (pint(i,j) .le. cmin) pint(i,j)=0.
    2401. 

  2401.                if (pint(i,j) .le. cmin2) pint(i,j)=0.
    2402. 

  2402.               pt(i,j)=pt(i,j)+ascp*pint(i,j)
    2403. 

  2403.               qv(i,j)=qv(i,j)-pint(i,j)
    2404. 

  2404.               qi(i,j)=qi(i,j)+pint(i,j)
    2405. 

  2405.            endif
    2406. 

  2406.         endif  ! if ( itaobraun.eq.1 )
    2407. 

  2407. !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    2408. 

  2408. !
    2409. 

  2409. !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    2410. 

  2410.         if ( itaobraun.eq.0 ) then
    2411. 

  2411.              tair(i,j)=(pt(i,j)+tb0)*pi0
    2412. 

  2412.              if (tair(i,j) .lt. t0) then
    2413. 

  2413.                if (qi(i,j) .le. cmin1) qi(i,j)=0.
    2414. 

  2414.                tairc(i,j)=tair(i,j)-t0
    2415. 

  2415.                dd(i,j)=r31r*exp(beta*tairc(i,j))
    2416. 

  2416.                rtair(i,j)=1./(tair(i,j)-c76)
    2417. 

  2417.                y2(i,j)=exp(c218-c580*rtair(i,j))
    2418. 

  2418.                qsi(i,j)=rp0*y2(i,j)
    2419. 

  2419.                esi(i,j)=c610*y2(i,j)
    2420. 

  2420.                ssi(i,j)=(qv(i,j)+qb0)/qsi(i,j)-1.
    2421. 

  2421.                dm(i,j)=max( (qv(i,j)+qb0-qsi(i,j)), 0.)
    2422. 

  2422.                rsub1(i,j)=cs580*qsi(i,j)*rtair(i,j)*rtair(i,j)
    2423. 

  2423.                dep(i,j)=dm(i,j)/(1.+rsub1(i,j))
    2424. 

  2424.               pint(i,j)=max(min(dd(i,j), dm(i,j)), 0.)
    2425. 

  2425.                y1(i,j)=1./tair(i,j)
    2426. 

  2426.                y2(i,j)=exp(betah*tairc(i,j))
    2427. 

  2427.                y3(i,j)=sqrt(qi(i,j))
    2428. 

  2428.                dd(i,j)=y1(i,j)*(rn10a*y1(i,j)-rn10b) &
    2429. 

  2429.                      +rn10c*tair(i,j)/esi(i,j)
    2430. 

  2430.              pidep(i,j)=max(r32rt*ssi(i,j)*y2(i,j)*y3(i,j)/dd(i,j), 0.)
    2431. 

  2431.               pint(i,j)=pint(i,j)+pidep(i,j)
    2432. 

  2432.               pint(i,j)=min(pint(i,j),dep(i,j))
    2433. 

  2433. !c          if (pint(i,j) .le. cmin2) pint(i,j)=0.
    2434. 

  2434.              pt(i,j)=pt(i,j)+ascp*pint(i,j)
    2435. 

  2435.              qv(i,j)=qv(i,j)-pint(i,j)
    2436. 

  2436.              qi(i,j)=qi(i,j)+pint(i,j)
    2437. 

  2437.             endif
    2438. 

  2438.         endif  ! if ( itaobraun.eq.0 )
    2439. 

  2439. !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    2440. 

  2440. 

  2441. !JJS  225    CONTINUE
    2442. 

  2442. 

  2443. !*****   TAO ET AL (1989) SATURATION TECHNIQUE  ***********************
    2444. 

  2444. 

  2445.          if (new_ice_sat .eq. 0) then
    2446. 

  2446. 

  2447. !JJS            DO 250 J=3,JLES
    2448. 

  2448. !JJS            DO 250 I=3,ILES
    2449. 

  2449.                tair(i,j)=(pt(i,j)+tb0)*pi0
    2450. 

  2450.                cnd(i,j)=rt0*(tair(i,j)-t00)
    2451. 

  2451.                dep(i,j)=rt0*(t0-tair(i,j))
    2452. 

  2452.                y1(i,j)=1./(tair(i,j)-c358)
    2453. 

  2453.                y2(i,j)=1./(tair(i,j)-c76)
    2454. 

  2454.                qsw(i,j)=rp0*exp(c172-c409*y1(i,j))
    2455. 

  2455.                qsi(i,j)=rp0*exp(c218-c580*y2(i,j))
    2456. 

  2456.                dd(i,j)=cp409*y1(i,j)*y1(i,j)
    2457. 

  2457.                dd1(i,j)=cp580*y2(i,j)*y2(i,j)
    2458. 

  2458.                if (qc(i,j).le.cmin) qc(i,j)=cmin
    2459. 

  2459.                if (qi(i,j).le.cmin) qi(i,j)=cmin
    2460. 

  2460.                if (tair(i,j).ge.t0) then
    2461. 

  2461.                   dep(i,j)=0.0
    2462. 

  2462.                   cnd(i,j)=1.
    2463. 

  2463.                   qi(i,j)=0.0
    2464. 

  2464.                endif
    2465. 

  2465. 

  2466.                if (tair(i,j).lt.t00) then
    2467. 

  2467.                   cnd(i,j)=0.0
    2468. 

  2468.                   dep(i,j)=1.
    2469. 

  2469.                   qc(i,j)=0.0
    2470. 

  2470.                endif
    2471. 

  2471. 

  2472.                y5(i,j)=avcp*cnd(i,j)+ascp*dep(i,j)
    2473. 

  2473. !               if (qc(i,j) .ge. cmin .or. qi(i,j) .ge. cmin) then
    2474. 

  2474.                y1(i,j)=qc(i,j)*qsw(i,j)/(qc(i,j)+qi(i,j))
    2475. 

  2475.                y2(i,j)=qi(i,j)*qsi(i,j)/(qc(i,j)+qi(i,j))
    2476. 

  2476.                y4(i,j)=dd(i,j)*y1(i,j)+dd1(i,j)*y2(i,j)
    2477. 

  2477.                qvs(i,j)=y1(i,j)+y2(i,j)
    2478. 

  2478.                rsub1(i,j)=(qv(i,j)+qb0-qvs(i,j))/(1.+y4(i,j)*y5(i,j))
    2479. 

  2479.                cnd(i,j)=cnd(i,j)*rsub1(i,j)
    2480. 

  2480.                dep(i,j)=dep(i,j)*rsub1(i,j)
    2481. 

  2481.                if (qc(i,j).le.cmin) qc(i,j)=0.
    2482. 

  2482.                if (qi(i,j).le.cmin) qi(i,j)=0.
    2483. 

  2483. !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    2484. 

  2484. !c    ******   condensation or evaporation of qc  ******
    2485. 

  2485. 

  2486.                cnd(i,j)=max(-qc(i,j),cnd(i,j))
    2487. 

  2487. 

  2488. !c    ******   deposition or sublimation of qi    ******
    2489. 

  2489. 

  2490.                dep(i,j)=max(-qi(i,j),dep(i,j))
    2491. 

  2491. 

  2492.                pt(i,j)=pt(i,j)+avcp*cnd(i,j)+ascp*dep(i,j)
    2493. 

  2493.                qv(i,j)=qv(i,j)-cnd(i,j)-dep(i,j)
    2494. 

  2494.                qc(i,j)=qc(i,j)+cnd(i,j)
    2495. 

  2495.                qi(i,j)=qi(i,j)+dep(i,j)
    2496. 

  2496. !JJS  250       continue
    2497. 

  2497.          endif
    2498. 

  2498. 

  2499.          if (new_ice_sat .eq. 1) then
    2500. 

  2500. 

  2501. !JJS            DO J=3,JLES
    2502. 

  2502. !JJS            DO I=3,ILES
    2503. 

  2503. 

  2504.                tair(i,j)=(pt(i,j)+tb0)*pi0
    2505. 

  2505.                cnd(i,j)=rt0*(tair(i,j)-t00)
    2506. 

  2506.                dep(i,j)=rt0*(t0-tair(i,j))
    2507. 

  2507.                y1(i,j)=1./(tair(i,j)-c358)
    2508. 

  2508.                y2(i,j)=1./(tair(i,j)-c76)
    2509. 

  2509.                qsw(i,j)=rp0*exp(c172-c409*y1(i,j))
    2510. 

  2510.                qsi(i,j)=rp0*exp(c218-c580*y2(i,j))
    2511. 

  2511.                dd(i,j)=cp409*y1(i,j)*y1(i,j)
    2512. 

  2512.                dd1(i,j)=cp580*y2(i,j)*y2(i,j)
    2513. 

  2513.                y5(i,j)=avcp*cnd(i,j)+ascp*dep(i,j)
    2514. 

  2514.                y1(i,j)=rt0*(tair(i,j)-t00)*qsw(i,j)
    2515. 

  2515.                y2(i,j)=rt0*(t0-tair(i,j))*qsi(i,j)
    2516. 

  2516. !               IF (QC(I,J).LE.CMIN) QC(I,J)=CMIN
    2517. 

  2517. !               IF (QI(I,J).LE.CMIN) QI(I,J)=CMIN
    2518. 

  2518. 

  2519.                if (tair(i,j).ge.t0) then
    2520. 

  2520. !                 QI(I,J)=0.0
    2521. 

  2521.                   dep(i,j)=0.0
    2522. 

  2522.                   cnd(i,j)=1.
    2523. 

  2523.                   y2(i,j)=0.
    2524. 

  2524.                   y1(i,j)=qsw(i,j)
    2525. 

  2525.                endif
    2526. 

  2526.                if (tair(i,j).lt.t00) then
    2527. 

  2527.                   cnd(i,j)=0.0
    2528. 

  2528.                   dep(i,j)=1.
    2529. 

  2529.                   y2(i,j)=qsi(i,j)
    2530. 

  2530.                   y1(i,j)=0.
    2531. 

  2531. !                 QC(I,J)=0.0
    2532. 

  2532.                endif
    2533. 

  2533. 

  2534. !            Y1(I,J)=QC(I,J)*QSW(I,J)/(QC(I,J)+QI(I,J))
    2535. 

  2535. !            Y2(I,J)=QI(I,J)*QSI(I,J)/(QC(I,J)+QI(I,J))
    2536. 

  2536. 

  2537.                y4(i,j)=dd(i,j)*y1(i,j)+dd1(i,j)*y2(i,j)
    2538. 

  2538.                qvs(i,j)=y1(i,j)+y2(i,j)
    2539. 

  2539.                rsub1(i,j)=(qv(i,j)+qb0-qvs(i,j))/(1.+y4(i,j)*y5(i,j))
    2540. 

  2540.                cnd(i,j)=cnd(i,j)*rsub1(i,j)
    2541. 

  2541.                dep(i,j)=dep(i,j)*rsub1(i,j)
    2542. 

  2542. !               IF (QC(I,J).LE.CMIN) QC(I,J)=0.
    2543. 

  2543. !               IF (QI(I,J).LE.CMIN) QI(I,J)=0.
    2544. 

  2544. 

  2545. !C    ******   CONDENSATION OR EVAPORATION OF QC  ******
    2546. 

  2546. 

  2547.                cnd(i,j)=max(-qc(i,j),cnd(i,j))
    2548. 

  2548. 

  2549. !C    ******   DEPOSITION OR SUBLIMATION OF QI    ******
    2550. 

  2550. 

  2551.                dep(i,j)=max(-qi(i,j),dep(i,j))
    2552. 

  2552. 

  2553.                pt(i,j)=pt(i,j)+avcp*cnd(i,j)+ascp*dep(i,j)
    2554. 

  2554.                qv(i,j)=qv(i,j)-cnd(i,j)-dep(i,j)
    2555. 

  2555.                qc(i,j)=qc(i,j)+cnd(i,j)
    2556. 

  2556.                qi(i,j)=qi(i,j)+dep(i,j)
    2557. 

  2557. !JJS            ENDDO
    2558. 

  2558. !JJS            ENDDO
    2559. 

  2559.          endif
    2560. 

  2560. 

  2561. !c
    2562. 

  2562. !
    2563. 

  2563.           if (new_ice_sat .eq. 2) then
    2564. 

  2564. !JJS          do j=3,jles
    2565. 

  2565. !JJS             do i=3,iles
    2566. 

  2566.           dep(i,j)=0.0
    2567. 

  2567.           cnd(i,j)=0.0
    2568. 

  2568.           tair(i,j)=(pt(i,j)+tb0)*pi0
    2569. 

  2569.           if (tair(i,j) .ge. 253.16) then
    2570. 

  2570.               y1(i,j)=1./(tair(i,j)-c358)
    2571. 

  2571.               qsw(i,j)=rp0*exp(c172-c409*y1(i,j))
    2572. 

  2572.               dd(i,j)=cp409*y1(i,j)*y1(i,j)
    2573. 

  2573.               dm(i,j)=qv(i,j)+qb0-qsw(i,j)
    2574. 

  2574.               cnd(i,j)=dm(i,j)/(1.+avcp*dd(i,j)*qsw(i,j))
    2575. 

  2575. !c    ******   condensation or evaporation of qc  ******
    2576. 

  2576.               cnd(i,j)=max(-qc(i,j), cnd(i,j))
    2577. 

  2577.              pt(i,j)=pt(i,j)+avcp*cnd(i,j)
    2578. 

  2578.              qv(i,j)=qv(i,j)-cnd(i,j)
    2579. 

  2579.              qc(i,j)=qc(i,j)+cnd(i,j)
    2580. 

  2580.          endif
    2581. 

  2581.           if (tair(i,j) .le. 258.16) then
    2582. 

  2582. !c             cnd(i,j)=0.0
    2583. 

  2583.            y2(i,j)=1./(tair(i,j)-c76)
    2584. 

  2584.            qsi(i,j)=rp0*exp(c218-c580*y2(i,j))
    2585. 

  2585.           dd1(i,j)=cp580*y2(i,j)*y2(i,j)
    2586. 

  2586.          dep(i,j)=(qv(i,j)+qb0-qsi(i,j))/(1.+ascp*dd1(i,j)*qsi(i,j))
    2587. 

  2587. !c    ******   deposition or sublimation of qi    ******
    2588. 

  2588.              dep(i,j)=max(-qi(i,j),dep(i,j))
    2589. 

  2589.              pt(i,j)=pt(i,j)+ascp*dep(i,j)
    2590. 

  2590.              qv(i,j)=qv(i,j)-dep(i,j)
    2591. 

  2591.              qi(i,j)=qi(i,j)+dep(i,j)
    2592. 

  2592.          endif
    2593. 

  2593. !JJS       enddo
    2594. 

  2594. !JJS       enddo
    2595. 

  2595.       endif
    2596. 

  2596. 

  2597. !c
    2598. 

  2598. !
    2599. 

  2599. !* 10 * PSDEP : DEPOSITION OR SUBLIMATION OF QS                   **10**
    2600. 

  2600. !* 20 * PGSUB : SUBLIMATION OF QG                                 **20**
    2601. 

  2601. 

  2602. !JJS         DO 275 J=3,JLES
    2603. 

  2603. !JJS         DO 275 I=3,ILES
    2604. 

  2604. 

  2605.             psdep(i,j)=0.0
    2606. 

  2606.             pgdep(i,j)=0.0
    2607. 

  2607.             pssub(i,j)=0.0
    2608. 

  2608.             pgsub(i,j)=0.0
    2609. 

  2609.             tair(i,j)=(pt(i,j)+tb0)*pi0
    2610. 

  2610. 

  2611.             if(tair(i,j).lt.t0) then
    2612. 

  2612.                if(qs(i,j).lt.cmin1) qs(i,j)=0.0
    2613. 

  2613.                if(qg(i,j).lt.cmin1) qg(i,j)=0.0
    2614. 

  2614.                rtair(i,j)=1./(tair(i,j)-c76)
    2615. 

  2615.                qsi(i,j)=rp0*exp(c218-c580*rtair(i,j))
    2616. 

  2616.                ssi(i,j)=(qv(i,j)+qb0)/qsi(i,j)-1.
    2617. 

  2617. !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    2618. 

  2618.                y1(i,j)=r10ar/(tca(i,j)*tair(i,j)**2)+1./(dwv(i,j) &
    2619. 

  2619.                       *qsi(i,j))
    2620. 

  2620.                y2(i,j)=.78/zs(i,j)**2+r101f*scv(i,j)/zs(i,j)**bsh5
    2621. 

  2621.                psdep(i,j)=r10t*ssi(i,j)*y2(i,j)/y1(i,j)
    2622. 

  2622.                pssub(i,j)=psdep(i,j)
    2623. 

  2623.                psdep(i,j)=r2is*max(psdep(i,j), 0.)
    2624. 

  2624.                pssub(i,j)=r2is*max(-qs(i,j), min(pssub(i,j), 0.))
    2625. 

  2625. 

  2626.                if(ihail.eq.1) then
    2627. 

  2627.                   y2(i,j)=.78/zg(i,j)**2+r20bq*scv(i,j)/zg(i,j)**bgh5
    2628. 

  2628.                else
    2629. 

  2629.                   y2(i,j)=.78/zg(i,j)**2+r20bs*scv(i,j)/zg(i,j)**bgh5
    2630. 

  2630.                endif
    2631. 

  2631. 

  2632.                pgsub(i,j)=r2ig*r20t*ssi(i,j)*y2(i,j)/y1(i,j)
    2633. 

  2633.                dm(i,j)=qv(i,j)+qb0-qsi(i,j)
    2634. 

  2634.                rsub1(i,j)=cs580*qsi(i,j)*rtair(i,j)*rtair(i,j)
    2635. 

  2635. 

  2636. !     ********   DEPOSITION OR SUBLIMATION OF QS  **********************
    2637. 

  2637. 

  2638.                y1(i,j)=dm(i,j)/(1.+rsub1(i,j))
    2639. 

  2639.                psdep(i,j)=r2is*min(psdep(i,j),max(y1(i,j),0.))
    2640. 

  2640.                y2(i,j)=min(y1(i,j),0.)
    2641. 

  2641.                pssub(i,j)=r2is*max(pssub(i,j),y2(i,j))
    2642. 

  2642. 

  2643. !     ********   SUBLIMATION OF QG   ***********************************
    2644. 

  2644. 

  2645.                dd(i,j)=max((-y2(i,j)-qs(i,j)), 0.)
    2646. 

  2646.               pgsub(i,j)=r2ig*min(dd(i,j), qg(i,j), max(pgsub(i,j),0.))
    2647. 

  2647. 

  2648.                if(qc(i,j)+qi(i,j).gt.1.e-5) then
    2649. 

  2649.                   dlt1(i,j)=1.
    2650. 

  2650.                else
    2651. 

  2651.                   dlt1(i,j)=0.
    2652. 

  2652.                endif
    2653. 

  2653. 

  2654.                psdep(i,j)=dlt1(i,j)*psdep(i,j)
    2655. 

  2655.                pssub(i,j)=(1.-dlt1(i,j))*pssub(i,j)
    2656. 

  2656.                pgsub(i,j)=(1.-dlt1(i,j))*pgsub(i,j)
    2657. 

  2657. 

  2658.                pt(i,j)=pt(i,j)+ascp*(psdep(i,j)+pssub(i,j)-pgsub(i,j))
    2659. 

  2659.                qv(i,j)=qv(i,j)+pgsub(i,j)-pssub(i,j)-psdep(i,j)
    2660. 

  2660.                qs(i,j)=qs(i,j)+psdep(i,j)+pssub(i,j)
    2661. 

  2661.                qg(i,j)=qg(i,j)-pgsub(i,j)
    2662. 

  2662.             endif
    2663. 

  2663. 

  2664. !* 23 * ERN : EVAPORATION OF QR (SUBSATURATION)                   **23**
    2665. 

  2665. 

  2666.             ern(i,j)=0.0
    2667. 

  2667. 

  2668.             if(qr(i,j).gt.0.0) then
    2669. 

  2669.                tair(i,j)=(pt(i,j)+tb0)*pi0
    2670. 

  2670.                rtair(i,j)=1./(tair(i,j)-c358)
    2671. 

  2671.                qsw(i,j)=rp0*exp(c172-c409*rtair(i,j))
    2672. 

  2672.                ssw(i,j)=(qv(i,j)+qb0)/qsw(i,j)-1.0
    2673. 

  2673.                dm(i,j)=qv(i,j)+qb0-qsw(i,j)
    2674. 

  2674.                rsub1(i,j)=cv409*qsw(i,j)*rtair(i,j)*rtair(i,j)
    2675. 

  2675.                dd1(i,j)=max(-dm(i,j)/(1.+rsub1(i,j)), 0.0)
    2676. 

  2676.                y1(i,j)=.78/zr(i,j)**2+r23af*scv(i,j)/zr(i,j)**bwh5
    2677. 

  2677.                y2(i,j)=r23br/(tca(i,j)*tair(i,j)**2)+1./(dwv(i,j) &
    2678. 

  2678.                        *qsw(i,j))
    2679. 

  2679. !cccc
    2680. 

  2680.                ern(i,j)=r23t*ssw(i,j)*y1(i,j)/y2(i,j)
    2681. 

  2681.                ern(i,j)=min(dd1(i,j),qr(i,j),max(ern(i,j),0.))
    2682. 

  2682.                pt(i,j)=pt(i,j)-avcp*ern(i,j)
    2683. 

  2683.                qv(i,j)=qv(i,j)+ern(i,j)
    2684. 

  2684.                qr(i,j)=qr(i,j)-ern(i,j)
    2685. 

  2685.             endif
    2686. 

  2686. 

  2687. !JJS 10/7/2008     vvvvv
    2688. 

  2688.     ENDIF    ! part of if (iwarm.eq.1) then
    2689. 

  2689. !JJS 10/7/2008     ^^^^^
    2690. 

  2690. 

  2691. !            IF (QV(I,J)+QB0 .LE. 0.) QV(I,J)=-QB0
    2692. 

  2692.             if (qc(i,j) .le. cmin1) qc(i,j)=0.
    2693. 

  2693.             if (qr(i,j) .le. cmin1) qr(i,j)=0.
    2694. 

  2694.             if (qi(i,j) .le. cmin1) qi(i,j)=0.
    2695. 

  2695.             if (qs(i,j) .le. cmin1) qs(i,j)=0.
    2696. 

  2696.             if (qg(i,j) .le. cmin1) qg(i,j)=0.
    2697. 

  2697.             dpt(i,j,k)=pt(i,j)
    2698. 

  2698.             dqv(i,j,k)=qv(i,j)
    2699. 

  2699.             qcl(i,j,k)=qc(i,j)
    2700. 

  2700.             qrn(i,j,k)=qr(i,j)
    2701. 

  2701.             qci(i,j,k)=qi(i,j)
    2702. 

  2702.             qcs(i,j,k)=qs(i,j)
    2703. 

  2703.             qcg(i,j,k)=qg(i,j)
    2704. 

  2704. 

  2705. !JJS  275    CONTINUE
    2706. 

  2706. 

  2707.          scc=0.
    2708. 

  2708.          see=0.
    2709. 

  2709. 

  2710. !         DO 110 J=3,JLES
    2711. 

  2711. !         DO 110 I=3,ILES
    2712. 

  2712. !            DD(I,J)=MAX(-CND(I,J), 0.)
    2713. 

  2713. !            CND(I,J)=MAX(CND(I,J), 0.)
    2714. 

  2714. !            DD1(I,J)=MAX(-DEP(I,J), 0.)
    2715. 

  2715. 

  2716. !ccshie 2/21/02 shie follow tao
    2717. 

  2717. !CC for reference    QI(I,J)=QI(I,J)-Y2(I,J)+PIDEP(I,J)*D2T
    2718. 

  2718. !CC for reference    QV(I,J)=QV(I,J)-(PIDEP(I,J))*D2T
    2719. 

  2719. 

  2720. !c            DEP(I,J)=MAX(DEP(I,J), 0.)
    2721. 

  2721. !            DEP(I,J)=MAX(DEP(I,J), 0.)+PIDEP(I,J)*D2T
    2722. 

  2722. !            SCC=SCC+CND(I,J)
    2723. 

  2723. !            SEE=SEE+DD(I,J)+ERN(I,J)
    2724. 

  2724. 

  2725. !  110    CONTINUE
    2726. 

  2726. 

  2727. !         SC(K)=SCC+SC(K)
    2728. 

  2728. !         SE(K)=SEE+SE(K)
    2729. 

  2729. 

  2730. !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    2731. 

  2731. !c     henry:  please take a look  (start)
    2732. 

  2732. !JJS modified by JJS on 5/1/2007  vvvvv
    2733. 

  2733. 

  2734. !JJS       do 305 j=3,jles
    2735. 

  2735. !JJS       do 305 i=3,iles
    2736. 

  2736.             dd(i,j)=max(-cnd(i,j), 0.)
    2737. 

  2737.             cnd(i,j)=max(cnd(i,j), 0.)
    2738. 

  2738.             dd1(i,j)=max(-dep(i,j), 0.)+pidep(i,j)*d2t
    2739. 

  2739.             dep(i,j)=max(dep(i,j), 0.)
    2740. 

  2740. !JJS  305  continue
    2741. 

  2741. 

  2742. !JJS       do 306 j=3,jles
    2743. 

  2743. !JJS       do 306 i=3,iles
    2744. 

  2744. !JJS              scc=scc+cnd(i,j)
    2745. 

  2745. !JJS              see=see+(dd(i,j)+ern(i,j))
    2746. 

  2746. !
    2747. 

  2747. !JJS            sddd=sddd+(dep(i,j)+pint(i,j)+psdep(i,j)+pgdep(i,j))
    2748. 

  2748. !JJS          ssss=ssss+dd1(i,j)
    2749. 

  2749. !JJS
    2750. 

  2750. !            shhh=shhh+rsw(i,j,k)*d2t
    2751. 

  2751. !            sccc=sccc+rlw(i,j,k)*d2t
    2752. 

  2752. !jjs
    2753. 

  2753. !JJS              smmm=smmm+(psmlt(i,j)+pgmlt(i,j)+pimlt(i,j))
    2754. 

  2754. !JJS              sfff=sfff+d2t*(psacw(i,j)+piacr(i,j)+psfw(i,j)
    2755. 

  2755. !JJS     1         +pgfr(i,j)+dgacw(i,j)+dgacr(i,j)+psacr(i,j))
    2756. 

  2756. !JJS     2        -qracs(i,j)+pihom(i,j)+pidw(i,j)
    2757. 

  2757. 

  2758. 

  2759.               sccc=cnd(i,j)
    2760. 

  2760.               seee=dd(i,j)+ern(i,j)
    2761. 

  2761.               sddd=dep(i,j)+pint(i,j)+psdep(i,j)+pgdep(i,j)
    2762. 

  2762.               ssss=dd1(i,j) + pgsub(i,j)
    2763. 

  2763.               smmm=psmlt(i,j)+pgmlt(i,j)+pimlt(i,j)
    2764. 

  2764.               sfff=d2t*(psacw(i,j)+piacr(i,j)+psfw(i,j) &
    2765. 

  2765.                +pgfr(i,j)+dgacw(i,j)+dgacr(i,j)+psacr(i,j) &
    2766. 

  2766.                +wgacr(i,j))+pihom(i,j)+pidw(i,j)
    2767. 

  2767. 

  2768. !           physc(i,k,j) = avcp * sccc / d2t
    2769. 

  2769. !           physe(i,k,j) = avcp * seee / d2t
    2770. 

  2770. !           physd(i,k,j) = ascp * sddd / d2t
    2771. 

  2771. !           physs(i,k,j) = ascp * ssss / d2t
    2772. 

  2772. !           physf(i,k,j) = afcp * sfff / d2t
    2773. 

  2773. !           physm(i,k,j) = afcp * smmm / d2t
    2774. 

  2774. !           physc(i,k,j) = physc(i,k,j) + avcp * sccc 
    2775. 

  2775. !           physe(i,k,j) = physc(i,k,j) + avcp * seee 
    2776. 

  2776. !           physd(i,k,j) = physd(i,k,j) + ascp * sddd 
    2777. 

  2777. !           physs(i,k,j) = physs(i,k,j) + ascp * ssss 
    2778. 

  2778. !           physf(i,k,j) = physf(i,k,j) + afcp * sfff 
    2779. 

  2779. !           physm(i,k,j) = physm(i,k,j) + afcp * smmm 
    2780. 

  2780. 

  2781. !JJS modified by JJS on 5/1/2007  ^^^^^
    2782. 

  2782. 

  2783.  2000 continue
    2784. 

  2784. 

  2785.  1000 continue
    2786. 

  2786. 

  2787. !JJS  ****************************************************************
    2788. 

  2788. !JJS  convert from GCE grid back to WRF grid
    2789. 

  2789.       do k=kts,kte
    2790. 

  2790.          do j=jts,jte
    2791. 

  2791.          do i=its,ite
    2792. 

  2792.          ptwrf(i,k,j) = dpt(i,j,k)
    2793. 

  2793.          qvwrf(i,k,j) = dqv(i,j,k)
    2794. 

  2794.          qlwrf(i,k,j) = qcl(i,j,k)
    2795. 

  2795.          qrwrf(i,k,j) = qrn(i,j,k)
    2796. 

  2796.          qiwrf(i,k,j) = qci(i,j,k)
    2797. 

  2797.          qswrf(i,k,j) = qcs(i,j,k)
    2798. 

  2798.          qgwrf(i,k,j) = qcg(i,j,k)
    2799. 

  2799.          enddo !i
    2800. 

  2800.          enddo !j
    2801. 

  2801.       enddo !k
    2802. 

  2802. 

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

  2804. 

  2805.       return
    2806. 

  2806.  END SUBROUTINE saticel_s
    2807. 

  2807. 

  2808. !JJS
    2809. 

  2809. !JJS      REAL FUNCTION GAMMA(X)
    2810. 

  2810. !JJS        Y=GAMMLN(X)
    2811. 

  2811. !JJS        GAMMA=EXP(Y)
    2812. 

  2812. !JJS      RETURN
    2813. 

  2813. !JJS      END
    2814. 

  2814. !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    2815. 

  2815. !JJS      real function GAMMLN (xx)
    2816. 

  2816.   real function gammagce (xx)
    2817. 

  2817. !**********************************************************************
    2818. 

  2818.   real*8 cof(6),stp,half,one,fpf,x,tmp,ser
    2819. 

  2819.   data cof,stp /  76.18009173,-86.50532033,24.01409822, &
    2820. 

  2820.      -1.231739516,.120858003e-2,-.536382e-5, 2.50662827465 /
    2821. 

  2821.   data half,one,fpf / .5, 1., 5.5 /
    2822. 

  2822. !
    2823. 

  2823.       x=xx-one
    2824. 

  2824.       tmp=x+fpf
    2825. 

  2825.       tmp=(x+half)*log(tmp)-tmp
    2826. 

  2826.       ser=one
    2827. 

  2827.       do  j=1,6
    2828. 

  2828.          x=x+one
    2829. 

  2829.         ser=ser+cof(j)/x
    2830. 

  2830.       enddo !j
    2831. 

  2831.       gammln=tmp+log(stp*ser)
    2832. 

  2832. !JJS
    2833. 

  2833.       gammagce=exp(gammln)
    2834. 

  2834. !JJS
    2835. 

  2835.       return
    2836. 

  2836.  END FUNCTION gammagce
    2837. 

  2837. 

  2838. END MODULE  module_mp_gsfcgce
    2839. 

  2839.