phys/module_mp_gsfcgce.F

   1 !WRF:MODEL_LAYER:PHYSICS
   2 !
   3
   4 MODULE module_mp_gsfcgce
   5
   6    USE     module_wrf_error
   7    USE module_mp_radar
   8
   9 !JJS 1/3/2008     vvvvv
  10
  11 !  common /bt/
  12    REAL,    PRIVATE ::          rd1,  rd2,   al,   cp
  13
  14 !  common /cont/
  15    REAL,    PRIVATE ::          c38, c358, c610, c149, &
  16                                c879, c172, c409,  c76, &
  17                                c218, c580, c141
  18 !  common /b3cs/
  19    REAL,    PRIVATE ::           ag,   bg,   as,   bs, &
  20                                  aw,   bw,  bgh,  bgq, &
  21                                 bsh,  bsq,  bwh,  bwq
  22
  23 !  common /size/
  24    REAL,    PRIVATE ::          tnw,  tns,  tng,       &
  25                                roqs, roqg, roqr
  26
  27 !  common /bterv/
  28    REAL,    PRIVATE ::          zrc,  zgc,  zsc,       &
  29                                 vrc,  vgc,  vsc
  30
  31 !  common /bsnw/
  32    REAL,    PRIVATE ::          alv,  alf,  als,   t0,   t00,     &
  33                                 avc,  afc,  asc,  rn1,  bnd1,     &
  34                                 rn2, bnd2,  rn3,  rn4,   rn5,     &
  35                                 rn6,  rn7,  rn8,  rn9,  rn10,     &
  36                               rn101,rn10a, rn11,rn11a,  rn12
  37
  38    REAL,    PRIVATE ::         rn14, rn15,rn15a, rn16,  rn17,     &
  39                               rn17a,rn17b,rn17c, rn18, rn18a,     &
  40                                rn19,rn19a,rn19b, rn20, rn20a,     &
  41                               rn20b, bnd3, rn21, rn22,  rn23,     &
  42                               rn23a,rn23b, rn25,rn30a, rn30b,     &
  43                               rn30c, rn31, beta, rn32
  44
  45    REAL,    PRIVATE, DIMENSION( 31 ) ::    rn12a, rn12b, rn13, rn25a
  46
  47 !  common /rsnw1/
  48    REAL,    PRIVATE ::         rn10b, rn10c, rnn191, rnn192,  rn30,     &
  49                              rnn30a,  rn33,  rn331,  rn332
  50
  51 !
  52    REAL,    PRIVATE, DIMENSION( 31 )  ::      aa1,  aa2
  53    DATA aa1/.7939e-7, .7841e-6, .3369e-5, .4336e-5, .5285e-5,     &
  54            .3728e-5, .1852e-5, .2991e-6, .4248e-6, .7434e-6,     &
  55            .1812e-5, .4394e-5, .9145e-5, .1725e-4, .3348e-4,     &
  56            .1725e-4, .9175e-5, .4412e-5, .2252e-5, .9115e-6,     &
  57            .4876e-6, .3473e-6, .4758e-6, .6306e-6, .8573e-6,     &
  58            .7868e-6, .7192e-6, .6513e-6, .5956e-6, .5333e-6,     &
  59            .4834e-6/
  60    DATA aa2/.4006, .4831, .5320, .5307, .5319,      &
  61            .5249, .4888, .3894, .4047, .4318,      &
  62            .4771, .5183, .5463, .5651, .5813,      &
  63            .5655, .5478, .5203, .4906, .4447,      &
  64            .4126, .3960, .4149, .4320, .4506,      &
  65            .4483, .4460, .4433, .4413, .4382,      &
  66            .4361/
  67
  68 !+---+-----------------------------------------------------------------+
  69 !..The following 6 variables moved here to facilitate reflectivity
  70 !.. calculation similar to other MP schemes, because when they get
  71 !.. declared later in the code (now commented out), it makes things
  72 !.. more difficult to integreate with the radar code.
  73       REAL    , PARAMETER ::     xnor = 8.0e6
  74       REAL    , PARAMETER ::     xnos = 1.6e7
  75       REAL    , PARAMETER ::     xnoh = 2.0e5
  76       REAL    , PARAMETER ::     xnog = 4.0e6
  77       REAL    , PARAMETER ::     rhohail = 917.
  78       REAL    , PARAMETER ::     rhograul = 400.
  79 !+---+-----------------------------------------------------------------+
  80
  81 !JJS 1/3/2008     ^^^^^
  82
  83 CONTAINS
  84
  85 !-------------------------------------------------------------------
  86 !  NASA/GSFC GCE
  87 !  Tao et al, 2001, Meteo. & Atmos. Phy., 97-137
  88 !-------------------------------------------------------------------
  89 !  SUBROUTINE gsfcgce(  th, th_old,                                 &
  90   SUBROUTINE gsfcgce(  th,                                         &
  91                        qv, ql,                                     &
  92                        qr, qi,                                     &
  93                        qs,                                         &
  94 !                       qvold, qlold,                               &
  95 !                       qrold, qiold,                               &
  96 !                       qsold,                                      &
  97                        rho, pii, p, dt_in, z,                      &
  98                        ht, dz8w, grav,                             &
  99                        rhowater, rhosnow,                          &
 100                        itimestep,                                  &
 101                        ids,ide, jds,jde, kds,kde,                  & ! domain dims
 102                        ims,ime, jms,jme, kms,kme,                  & ! memory dims
 103                        its,ite, jts,jte, kts,kte,                  & ! tile   dims
 104                        rainnc, rainncv,                            &
 105                        snownc, snowncv, sr,                        &
 106                        graupelnc, graupelncv,                      &
 107                        refl_10cm, diagflag, do_radar_ref,          &
 108 !                       f_qg, qg, pgold,                            &
 109                        f_qg, qg,                                   &
 110                        ihail, ice2                                 &
 111                                                                    )
 112
 113 !-------------------------------------------------------------------
 114   IMPLICIT NONE
 115 !-------------------------------------------------------------------
 116 !
 117 ! JJS 2/15/2005
 118 !
 119   INTEGER,      INTENT(IN   )    ::   ids,ide, jds,jde, kds,kde , &
 120                                       ims,ime, jms,jme, kms,kme , &
 121                                       its,ite, jts,jte, kts,kte
 122   INTEGER,      INTENT(IN   )    ::   itimestep, ihail, ice2
 123
 124   REAL, DIMENSION( ims:ime , kms:kme , jms:jme ),                 &
 125         INTENT(INOUT) ::                                          &
 126                                                               th, &
 127                                                               qv, &
 128                                                               ql, &
 129                                                               qr, &
 130                                                               qi, &
 131                                                               qs, &
 132                                                               qg
 133 !
 134   REAL, DIMENSION( ims:ime , kms:kme , jms:jme ),                 &
 135         INTENT(IN   ) ::                                          &
 136 !                                                         th_old, &
 137 !                                                          qvold, &
 138 !                                                          qlold, &
 139 !                                                          qrold, &
 140 !                                                          qiold, &
 141 !                                                          qsold, &
 142 !                                                          qgold, &
 143                                                              rho, &
 144                                                              pii, &
 145                                                                p, &
 146                                                             dz8w, &
 147                                                                z
 148
 149   REAL, DIMENSION( ims:ime , jms:jme ),                           &
 150         INTENT(INOUT) ::                               rainnc,    &
 151                                                        rainncv,   &
 152                                                        snownc,    &
 153                                                        snowncv,   &
 154                                                        sr,        &
 155                                                        graupelnc, &
 156                                                        graupelncv
 157
 158 !+---+-----------------------------------------------------------------+
 159   REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT)::           &  ! GT
 160                                                        refl_10cm
 161   LOGICAL, OPTIONAL, INTENT(IN) :: diagflag
 162   INTEGER, OPTIONAL, INTENT(IN) :: do_radar_ref
 163 !+---+-----------------------------------------------------------------+
 164
 165   REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN) ::       ht
 166
 167   REAL, INTENT(IN   ) ::                                   dt_in, &
 168                                                             grav, &
 169                                                         rhowater, &
 170                                                          rhosnow
 171
 172   LOGICAL, INTENT(IN), OPTIONAL :: F_QG
 173
 174 !  LOCAL VAR
 175
 176
 177 !jjs  INTEGER             ::                            min_q, max_q
 178
 179 !jjs  REAL, DIMENSION( its:ite , jts:jte )                            &
 180 !jjs                               ::        rain, snow, graupel,ice
 181
 182 !
 183 !  INTEGER :: IHAIL, itaobraun, ice2, istatmin, new_ice_sat, id
 184   INTEGER ::  itaobraun, istatmin, new_ice_sat, id
 185
 186   INTEGER :: i, j, k
 187   INTEGER :: iskip, ih, icount, ibud, i24h
 188   REAL    :: hour
 189   REAL    , PARAMETER :: cmin=1.e-20
 190   REAL    :: dth, dqv, dqrest, dqall, dqall1, rhotot, a1, a2
 191 !  REAL, DIMENSION( its:ite , kts:kte , jts:jte ) ::                   &
 192 !                         th1, qv1, ql1, qr1, qi1, qs1, qg1
 193
 194   LOGICAL :: flag_qg
 195
 196 !+---+-----------------------------------------------------------------+
 197
 198       INTEGER:: NCALL=0
 199
 200       IF (NCALL .EQ. 0) THEN
 201 !..Set these variables needed for computing radar reflectivity.  These
 202 !.. get used within radar_init to create other variables used in the
 203 !.. radar module.
 204          xam_r = 3.14159*rhowater/6.
 205          xbm_r = 3.
 206          xmu_r = 0.
 207          xam_s = 3.14159*rhosnow/6.
 208          xbm_s = 3.
 209          xmu_s = 0.
 210          if (ihail .eq. 1) then
 211             xam_g = 3.14159*rhohail/6.
 212          else
 213             xam_g = 3.14159*rhograul/6.
 214          endif
 215          xbm_g = 3.
 216          xmu_g = 0.
 217
 218          call radar_init
 219          NCALL = 1
 220       ENDIF
 221 !+---+-----------------------------------------------------------------+
 222
 223 !
 224 !c  ihail = 0    for graupel, for tropical region
 225 !c  ihail = 1    for hail, for mid-lat region
 226
 227 ! itaobraun: 0 for Tao's constantis, 1 for Braun's constants
 228 !c        if ( itaobraun.eq.1 ) --> betah=0.5*beta=-.46*0.5=-0.23;   cn0=1.e-6
 229 !c        if ( itaobraun.eq.0 ) --> betah=0.5*beta=-.6*0.5=-0.30;    cn0=1.e-8
 230    itaobraun = 1
 231
 232 !c  ice2 = 0    for 3 ice --- ice, snow and graupel/hail
 233 !c  ice2 = 1    for 2 ice --- ice and snow only
 234 !c  ice2 = 2    for 2 ice --- ice and graupel only, use ihail = 0 only
 235 !c  ice2 = 3    for 0 ice --- no ice, warm only
 236
 237 !  if (ice2 .eq. 2) ihail = 0
 238
 239   i24h=nint(86400./dt_in)
 240   if (mod(itimestep,i24h).eq.1) then
 241      write(6,*) 'ihail=',ihail,'  ice2=',ice2
 242      if (ice2.eq.0) then
 243         write(6,*) 'Running 3-ice scheme in GSFCGCE with'
 244         if (ihail.eq.0) then
 245            write(6,*) '     ice, snow and graupel'
 246         else if (ihail.eq.1) then
 247                 write(6,*) '     ice, snow and hail'
 248         else
 249              write(6,*) 'ihail has to be either 1 or 0'
 250              stop
 251         endif !ihail
 252      else if (ice2.eq.1) then
 253              write(6,*) 'Running 2-ice scheme in GSFCGCE with'
 254              write(6,*) '     ice and snow'
 255      else if (ice2.eq.2) then
 256              write(6,*) 'Running 2-ice scheme in GSFCGCE with'
 257              write(6,*) '     ice and graupel'
 258      else if (ice2.eq.3) then
 259              write(6,*) 'Running warm rain only scheme in GSFCGCE without any ice'
 260      else
 261              write(6,*) 'gsfcgce_2ice in namelist.input has to be 0, 1, 2, or 3'
 262              stop
 263      endif !ice2
 264   endif !itimestep
 265
 266 !c  new_ice_sat = 0, 1 or 2
 267     new_ice_sat = 2
 268
 269 !c istatmin
 270     istatmin = 180
 271
 272 !c id = 0  without in-line staticstics
 273 !c id = 1  with in-line staticstics
 274     id = 0
 275
 276 !c ibud = 0 no calculation of dth, dqv, dqrest and dqall
 277 !c ibud = 1 yes
 278     ibud = 0
 279
 280 !jjs   dt=dt_in
 281 !jjs   rhoe_s=1.29
 282 !
 283 !   IF (P_QI .lt. P_FIRST_SCALAR .or. P_QS .lt. P_FIRST_SCALAR) THEN
 284 !      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.')
 285 !   ENDIF
 286
 287 ! calculte fallflux and precipiation in MKS system
 288
 289    call fall_flux(dt_in, qr, qi, qs, qg, p,                   &
 290                       rho, z, dz8w, ht, rainnc,               &
 291                       rainncv, grav,itimestep,                &
 292                       rhowater, rhosnow,                      &
 293                       snownc, snowncv, sr,                    &
 294                       graupelnc, graupelncv,                  &
 295                       ihail, ice2,                            &
 296                       ims,ime, jms,jme, kms,kme,              & ! memory dims
 297                       its,ite, jts,jte, kts,kte               ) ! tile   dims
 298 !-----------------------------------------------------------------------
 299
 300 !c  set up constants used internally in GCE
 301
 302    call consat_s (ihail, itaobraun)
 303
 304
 305 !c Negative values correction
 306
 307    iskip = 1
 308
 309    if (iskip.eq.0) then
 310       call negcor(qv,rho,dz8w,ims,ime,jms,jme,kms,kme, &
 311                            itimestep,1,             &
 312                            its,ite,jts,jte,kts,kte)
 313       call negcor(ql,rho,dz8w,ims,ime,jms,jme,kms,kme, &
 314                            itimestep,2,             &
 315                            its,ite,jts,jte,kts,kte)
 316       call negcor(qr,rho,dz8w,ims,ime,jms,jme,kms,kme, &
 317                            itimestep,3,             &
 318                            its,ite,jts,jte,kts,kte)
 319       call negcor(qi,rho,dz8w,ims,ime,jms,jme,kms,kme, &
 320                            itimestep,4,             &
 321                            its,ite,jts,jte,kts,kte)
 322       call negcor(qs,rho,dz8w,ims,ime,jms,jme,kms,kme, &
 323                            itimestep,5,             &
 324                            its,ite,jts,jte,kts,kte)
 325       call negcor(qg,rho,dz8w,ims,ime,jms,jme,kms,kme, &
 326                            itimestep,6,             &
 327                            its,ite,jts,jte,kts,kte)
 328 !   else if (mod(itimestep,i24h).eq.1) then
 329 !      print *,'no neg correction in mp at timestep=',itimestep
 330    endif ! iskip
 331
 332 !c microphysics in GCE
 333
 334    call SATICEL_S( dt_in, IHAIL, itaobraun, ICE2, istatmin,     &
 335                    new_ice_sat, id,                             &
 336 !                   th, th_old, qv, ql, qr,                      &
 337                    th, qv, ql, qr,                      &
 338                    qi, qs, qg,                                  &
 339 !                   qvold, qlold, qrold,                         &
 340 !                   qiold, qsold, qgold,                         &
 341                    rho, pii, p, itimestep,                      &
 342                    refl_10cm, diagflag, do_radar_ref,           & ! GT added for reflectivity calcs
 343                    ids,ide, jds,jde, kds,kde,                   & ! domain dims
 344                    ims,ime, jms,jme, kms,kme,                   & ! memory dims
 345                    its,ite, jts,jte, kts,kte                    & ! tile   dims
 346                                                                 )
 347
 348
 349    END SUBROUTINE gsfcgce
 350
 351 !-----------------------------------------------------------------------
 352    SUBROUTINE fall_flux ( dt, qr, qi, qs, qg, p,              &
 353                       rho, z, dz8w, topo, rainnc,             &
 354                       rainncv, grav, itimestep,               &
 355                       rhowater, rhosnow,                      &
 356                       snownc, snowncv, sr,                    &
 357                       graupelnc, graupelncv,                  &
 358                       ihail, ice2,                            &
 359                       ims,ime, jms,jme, kms,kme,              & ! memory dims
 360                       its,ite, jts,jte, kts,kte               ) ! tile   dims
 361 !-----------------------------------------------------------------------
 362 ! adopted from Jiun-Dar Chern's codes for Purdue Regional Model
 363 ! adopted by Jainn J. Shi, 6/10/2005
 364 !-----------------------------------------------------------------------
 365
 366   IMPLICIT NONE
 367   INTEGER, INTENT(IN   )               :: ihail, ice2,                &
 368                                           ims,ime, jms,jme, kms,kme,  &
 369                                           its,ite, jts,jte, kts,kte
 370   INTEGER, INTENT(IN   )               :: itimestep
 371   REAL,    DIMENSION( ims:ime , kms:kme , jms:jme ),                  &
 372            INTENT(INOUT)               :: qr, qi, qs, qg
 373   REAL,    DIMENSION( ims:ime , jms:jme ),                            &
 374            INTENT(INOUT)               :: rainnc, rainncv,            &
 375                                           snownc, snowncv, sr,        &
 376                                           graupelnc, graupelncv
 377   REAL,    DIMENSION( ims:ime , kms:kme , jms:jme ),                  &
 378            INTENT(IN   )               :: rho, z, dz8w, p
 379
 380   REAL,    INTENT(IN   )               :: dt, grav, rhowater, rhosnow
 381
 382
 383   REAL,    DIMENSION( ims:ime , jms:jme ),                            &
 384            INTENT(IN   )               :: topo
 385
 386 ! temperary vars
 387
 388   REAL,    DIMENSION( kts:kte )           :: sqrhoz
 389   REAL                                    :: tmp1, term0
 390   REAL                                :: pptrain, pptsnow,        &
 391                                          pptgraul, pptice
 392   REAL,    DIMENSION( kts:kte )       :: qrz, qiz, qsz, qgz,      &
 393                                          zz, dzw, prez, rhoz,     &
 394                                          orhoz
 395
 396
 397    INTEGER                    :: k, i, j
 398 !
 399
 400   REAL, DIMENSION( kts:kte )    :: vtr, vts, vtg, vti
 401
 402   REAL                          :: dtb, pi, consta, constc, gambp4,    &
 403                                    gamdp4, gam4pt5, gam4bbar
 404
 405 !  Lin
 406 !-GT  REAL    , PARAMETER ::     xnor = 8.0e6
 407 !   REAL    , PARAMETER ::     xnos = 3.0e6
 408 !-GT   REAL    , PARAMETER ::     xnos = 1.6e7   ! Tao's value
 409    REAL    , PARAMETER ::                              &
 410 !             constb = 0.8, constd = 0.25, o6 = 1./6.,           &
 411              constb = 0.8, constd = 0.11, o6 = 1./6.,           &
 412              cdrag = 0.6
 413 !  Lin
 414 !  REAL    , PARAMETER ::     xnoh = 4.0e4
 415 !-GT  REAL    , PARAMETER ::     xnoh = 2.0e5    ! Tao's value
 416 !-GT  REAL    , PARAMETER ::     rhohail = 917.
 417
 418 ! Hobbs
 419 !-GT  REAL    , PARAMETER ::     xnog = 4.0e6
 420 !-GT  REAL    , PARAMETER ::     rhograul = 400.
 421   REAL    , PARAMETER ::     abar = 19.3, bbar = 0.37,      &
 422                                       p0 = 1.0e5
 423
 424   REAL    , PARAMETER ::     rhoe_s = 1.29
 425
 426 ! for terminal velocity flux
 427   INTEGER                       :: min_q, max_q
 428   REAL                          :: t_del_tv, del_tv, flux, fluxin, fluxout ,tmpqrz
 429   LOGICAL                       :: notlast
 430
 431 !  if (itimestep.eq.1) then
 432 !     write(6, *) 'in fall_flux'
 433 !     write(6, *) 'ims=', ims, '  ime=', ime
 434 !     write(6, *) 'jms=', jms, '  jme=', jme
 435 !     write(6, *) 'kms=', kms, '  kme=', kme
 436 !     write(6, *) 'its=', its, '  ite=', ite
 437 !     write(6, *) 'jts=', jts, '  jte=', jte
 438 !     write(6, *) 'kts=', kts, '  kte=', kte
 439 !     write(6, *) 'dt=', dt
 440 !     write(6, *) 'ihail=', ihail
 441 !     write(6, *) 'ICE2=', ICE2
 442 !     write(6, *) 'dt=', dt
 443 !   endif
 444
 445 !-----------------------------------------------------------------------
 446 !  This program calculates precipitation fluxes due to terminal velocities.
 447 !-----------------------------------------------------------------------
 448
 449    dtb=dt
 450    pi=acos(-1.)
 451    consta=2115.0*0.01**(1-constb)
 452 !   constc=152.93*0.01**(1-constd)
 453    constc=78.63*0.01**(1-constd)
 454
 455 !  Gamma function
 456    gambp4=ggamma(constb+4.)
 457    gamdp4=ggamma(constd+4.)
 458    gam4pt5=ggamma(4.5)
 459    gam4bbar=ggamma(4.+bbar)
 460
 461 !***********************************************************************
 462 ! Calculate precipitation fluxes due to terminal velocities.
 463 !***********************************************************************
 464 !
 465 !- Calculate termianl velocity (vt?)  of precipitation q?z
 466 !- Find maximum vt? to determine the small delta t
 467
 468  j_loop:  do j = jts, jte
 469  i_loop:  do i = its, ite
 470
 471    pptrain = 0.
 472    pptsnow = 0.
 473    pptgraul = 0.
 474    pptice  = 0.
 475
 476    do k = kts, kte
 477       qrz(k)=qr(i,k,j)
 478       rhoz(k)=rho(i,k,j)
 479       orhoz(k)=1./rhoz(k)
 480       prez(k)=p(i,k,j)
 481       sqrhoz(k)=sqrt(rhoe_s/rhoz(k))
 482       zz(k)=z(i,k,j)
 483       dzw(k)=dz8w(i,k,j)
 484    enddo !k
 485
 486       DO k = kts, kte
 487          qiz(k)=qi(i,k,j)
 488       ENDDO
 489
 490       DO k = kts, kte
 491          qsz(k)=qs(i,k,j)
 492       ENDDO
 493
 494    IF (ice2 .eq. 0) THEN
 495       DO k = kts, kte
 496          qgz(k)=qg(i,k,j)
 497       ENDDO
 498    ELSE
 499       DO k = kts, kte
 500          qgz(k)=0.
 501       ENDDO
 502    ENDIF
 503
 504
 505 !
 506 !-- rain
 507 !
 508     t_del_tv=0.
 509     del_tv=dtb
 510     notlast=.true.
 511     DO while (notlast)
 512 !
 513       min_q=kte
 514       max_q=kts-1
 515 !
 516       do k=kts,kte-1
 517          if (qrz(k) .gt. 1.0e-8) then
 518             min_q=min0(min_q,k)
 519             max_q=max0(max_q,k)
 520             tmp1=sqrt(pi*rhowater*xnor/rhoz(k)/qrz(k))
 521             tmp1=sqrt(tmp1)
 522             vtr(k)=consta*gambp4*sqrhoz(k)/tmp1**constb
 523             vtr(k)=vtr(k)/6.
 524             if (k .eq. 1) then
 525                del_tv=amin1(del_tv,0.9*(zz(k)-topo(i,j))/vtr(k))
 526             else
 527                del_tv=amin1(del_tv,0.9*(zz(k)-zz(k-1))/vtr(k))
 528             endif
 529          else
 530             vtr(k)=0.
 531          endif
 532       enddo
 533
 534       if (max_q .ge. min_q) then
 535 !
 536 !- Check if the summation of the small delta t >=  big delta t
 537 !             (t_del_tv)          (del_tv)             (dtb)
 538
 539          t_del_tv=t_del_tv+del_tv
 540 !
 541          if ( t_del_tv .ge. dtb ) then
 542               notlast=.false.
 543               del_tv=dtb+del_tv-t_del_tv
 544          endif
 545
 546 ! use small delta t to calculate the qrz flux
 547 ! termi is the qrz flux pass in the grid box through the upper boundary
 548 ! termo is the qrz flux pass out the grid box through the lower boundary
 549 !
 550          fluxin=0.
 551          do k=max_q,min_q,-1
 552             fluxout=rhoz(k)*vtr(k)*qrz(k)
 553             flux=(fluxin-fluxout)/rhoz(k)/dzw(k)
 554 !            tmpqrz=qrz(k)
 555             qrz(k)=qrz(k)+del_tv*flux
 556             qrz(k)=amax1(0.,qrz(k))
 557             qr(i,k,j)=qrz(k)
 558             fluxin=fluxout
 559          enddo
 560          if (min_q .eq. 1) then
 561             pptrain=pptrain+fluxin*del_tv
 562          else
 563             qrz(min_q-1)=qrz(min_q-1)+del_tv*  &
 564                           fluxin/rhoz(min_q-1)/dzw(min_q-1)
 565             qr(i,min_q-1,j)=qrz(min_q-1)
 566          endif
 567 !
 568       else
 569          notlast=.false.
 570       endif
 571     ENDDO
 572
 573 !
 574 !-- snow
 575 !
 576     t_del_tv=0.
 577     del_tv=dtb
 578     notlast=.true.
 579
 580     DO while (notlast)
 581 !
 582       min_q=kte
 583       max_q=kts-1
 584 !
 585       do k=kts,kte-1
 586          if (qsz(k) .gt. 1.0e-8) then
 587             min_q=min0(min_q,k)
 588             max_q=max0(max_q,k)
 589             tmp1=sqrt(pi*rhosnow*xnos/rhoz(k)/qsz(k))
 590             tmp1=sqrt(tmp1)
 591             vts(k)=constc*gamdp4*sqrhoz(k)/tmp1**constd
 592             vts(k)=vts(k)/6.
 593             if (k .eq. 1) then
 594                del_tv=amin1(del_tv,0.9*(zz(k)-topo(i,j))/vts(k))
 595             else
 596                del_tv=amin1(del_tv,0.9*(zz(k)-zz(k-1))/vts(k))
 597             endif
 598          else
 599             vts(k)=0.
 600          endif
 601       enddo
 602
 603       if (max_q .ge. min_q) then
 604 !
 605 !
 606 !- Check if the summation of the small delta t >=  big delta t
 607 !             (t_del_tv)          (del_tv)             (dtb)
 608
 609          t_del_tv=t_del_tv+del_tv
 610
 611          if ( t_del_tv .ge. dtb ) then
 612               notlast=.false.
 613               del_tv=dtb+del_tv-t_del_tv
 614          endif
 615
 616 ! use small delta t to calculate the qsz flux
 617 ! termi is the qsz flux pass in the grid box through the upper boundary
 618 ! termo is the qsz flux pass out the grid box through the lower boundary
 619 !
 620          fluxin=0.
 621          do k=max_q,min_q,-1
 622             fluxout=rhoz(k)*vts(k)*qsz(k)
 623             flux=(fluxin-fluxout)/rhoz(k)/dzw(k)
 624             qsz(k)=qsz(k)+del_tv*flux
 625             qsz(k)=amax1(0.,qsz(k))
 626             qs(i,k,j)=qsz(k)
 627             fluxin=fluxout
 628          enddo
 629          if (min_q .eq. 1) then
 630             pptsnow=pptsnow+fluxin*del_tv
 631          else
 632             qsz(min_q-1)=qsz(min_q-1)+del_tv*  &
 633                          fluxin/rhoz(min_q-1)/dzw(min_q-1)
 634             qs(i,min_q-1,j)=qsz(min_q-1)
 635          endif
 636 !
 637       else
 638          notlast=.false.
 639       endif
 640
 641     ENDDO
 642
 643 !
 644 !   ice2=0 --- with hail/graupel
 645 !   ice2=1 --- without hail/graupel
 646 !
 647   if (ice2.eq.0) then
 648 !
 649 !-- If IHAIL=1, use hail.
 650 !-- If IHAIL=0, use graupel.
 651 !
 652 !    if (ihail .eq. 1) then
 653 !       xnog = xnoh
 654 !       rhograul = rhohail
 655 !    endif
 656
 657     t_del_tv=0.
 658     del_tv=dtb
 659     notlast=.true.
 660 !
 661     DO while (notlast)
 662 !
 663       min_q=kte
 664       max_q=kts-1
 665 !
 666       do k=kts,kte-1
 667          if (qgz(k) .gt. 1.0e-8) then
 668             if (ihail .eq. 1) then
 669 !  for hail, based on Lin et al (1983)
 670               min_q=min0(min_q,k)
 671               max_q=max0(max_q,k)
 672               tmp1=sqrt(pi*rhohail*xnoh/rhoz(k)/qgz(k))
 673               tmp1=sqrt(tmp1)
 674               term0=sqrt(4.*grav*rhohail/3./rhoz(k)/cdrag)
 675               vtg(k)=gam4pt5*term0*sqrt(1./tmp1)
 676               vtg(k)=vtg(k)/6.
 677               if (k .eq. 1) then
 678                  del_tv=amin1(del_tv,0.9*(zz(k)-topo(i,j))/vtg(k))
 679               else
 680                  del_tv=amin1(del_tv,0.9*(zz(k)-zz(k-1))/vtg(k))
 681               endif !k
 682             else
 683 ! added by JJS
 684 ! for graupel, based on RH (1984)
 685               min_q=min0(min_q,k)
 686               max_q=max0(max_q,k)
 687               tmp1=sqrt(pi*rhograul*xnog/rhoz(k)/qgz(k))
 688               tmp1=sqrt(tmp1)
 689               tmp1=tmp1**bbar
 690               tmp1=1./tmp1
 691               term0=abar*gam4bbar/6.
 692               vtg(k)=term0*tmp1*(p0/prez(k))**0.4
 693               if (k .eq. 1) then
 694                  del_tv=amin1(del_tv,0.9*(zz(k)-topo(i,j))/vtg(k))
 695               else
 696                  del_tv=amin1(del_tv,0.9*(zz(k)-zz(k-1))/vtg(k))
 697               endif !k
 698             endif !ihail
 699          else
 700             vtg(k)=0.
 701          endif !qgz
 702       enddo !k
 703
 704       if (max_q .ge. min_q) then
 705 !
 706 !
 707 !- Check if the summation of the small delta t >=  big delta t
 708 !             (t_del_tv)          (del_tv)             (dtb)
 709
 710          t_del_tv=t_del_tv+del_tv
 711
 712          if ( t_del_tv .ge. dtb ) then
 713               notlast=.false.
 714               del_tv=dtb+del_tv-t_del_tv
 715          endif
 716
 717 ! use small delta t to calculate the qgz flux
 718 ! termi is the qgz flux pass in the grid box through the upper boundary
 719 ! termo is the qgz flux pass out the grid box through the lower boundary
 720 !
 721          fluxin=0.
 722          do k=max_q,min_q,-1
 723             fluxout=rhoz(k)*vtg(k)*qgz(k)
 724             flux=(fluxin-fluxout)/rhoz(k)/dzw(k)
 725             qgz(k)=qgz(k)+del_tv*flux
 726             qgz(k)=amax1(0.,qgz(k))
 727             qg(i,k,j)=qgz(k)
 728             fluxin=fluxout
 729          enddo
 730          if (min_q .eq. 1) then
 731             pptgraul=pptgraul+fluxin*del_tv
 732          else
 733             qgz(min_q-1)=qgz(min_q-1)+del_tv*  &
 734                          fluxin/rhoz(min_q-1)/dzw(min_q-1)
 735             qg(i,min_q-1,j)=qgz(min_q-1)
 736          endif
 737 !
 738       else
 739          notlast=.false.
 740       endif
 741 !
 742     ENDDO
 743  ENDIF !ice2
 744 !
 745 !-- cloud ice  (03/21/02) follow Vaughan T.J. Phillips at GFDL
 746 !
 747
 748     t_del_tv=0.
 749     del_tv=dtb
 750     notlast=.true.
 751 !
 752     DO while (notlast)
 753 !
 754       min_q=kte
 755       max_q=kts-1
 756 !
 757       do k=kts,kte-1
 758          if (qiz(k) .gt. 1.0e-8) then
 759             min_q=min0(min_q,k)
 760             max_q=max0(max_q,k)
 761             vti(k)= 3.29 * (rhoz(k)* qiz(k))** 0.16  ! Heymsfield and Donner
 762             if (k .eq. 1) then
 763                del_tv=amin1(del_tv,0.9*(zz(k)-topo(i,j))/vti(k))
 764             else
 765                del_tv=amin1(del_tv,0.9*(zz(k)-zz(k-1))/vti(k))
 766             endif
 767          else
 768             vti(k)=0.
 769          endif
 770       enddo
 771
 772       if (max_q .ge. min_q) then
 773 !
 774 !
 775 !- Check if the summation of the small delta t >=  big delta t
 776 !             (t_del_tv)          (del_tv)             (dtb)
 777
 778          t_del_tv=t_del_tv+del_tv
 779
 780          if ( t_del_tv .ge. dtb ) then
 781               notlast=.false.
 782               del_tv=dtb+del_tv-t_del_tv
 783          endif
 784
 785 ! use small delta t to calculate the qiz flux
 786 ! termi is the qiz flux pass in the grid box through the upper boundary
 787 ! termo is the qiz flux pass out the grid box through the lower boundary
 788 !
 789
 790          fluxin=0.
 791          do k=max_q,min_q,-1
 792             fluxout=rhoz(k)*vti(k)*qiz(k)
 793             flux=(fluxin-fluxout)/rhoz(k)/dzw(k)
 794             qiz(k)=qiz(k)+del_tv*flux
 795             qiz(k)=amax1(0.,qiz(k))
 796             qi(i,k,j)=qiz(k)
 797             fluxin=fluxout
 798          enddo
 799          if (min_q .eq. 1) then
 800             pptice=pptice+fluxin*del_tv
 801          else
 802             qiz(min_q-1)=qiz(min_q-1)+del_tv*  &
 803                          fluxin/rhoz(min_q-1)/dzw(min_q-1)
 804             qi(i,min_q-1,j)=qiz(min_q-1)
 805          endif
 806 !
 807       else
 808          notlast=.false.
 809       endif
 810 !
 811    ENDDO !notlast
 812
 813 !   prnc(i,j)=prnc(i,j)+pptrain
 814 !   psnowc(i,j)=psnowc(i,j)+pptsnow
 815 !   pgrauc(i,j)=pgrauc(i,j)+pptgraul
 816 !   picec(i,j)=picec(i,j)+pptice
 817 !
 818
 819 !   write(6,*) 'i=',i,' j=',j,'   ', pptrain, pptsnow, pptgraul, pptice
 820 !   flush(6)
 821
 822    snowncv(i,j) = pptsnow
 823    snownc(i,j) = snownc(i,j) + pptsnow
 824    graupelncv(i,j) = pptgraul
 825    graupelnc(i,j) = graupelnc(i,j) + pptgraul
 826    RAINNCV(i,j) = pptrain + pptsnow + pptgraul + pptice
 827    RAINNC(i,j)  = RAINNC(i,j) + pptrain + pptsnow + pptgraul + pptice
 828    sr(i,j) = 0.
 829    if (RAINNCV(i,j) .gt. 0.) sr(i,j) = (pptsnow + pptgraul + pptice) / RAINNCV(i,j)
 830
 831   ENDDO i_loop
 832   ENDDO j_loop
 833
 834 !  if (itimestep.eq.6480) then
 835 !     write(51,*) 'in the end of fallflux, itimestep=',itimestep
 836 !     do j=jts,jte
 837 !        do i=its,ite
 838 !           if (rainnc(i,j).gt.400.) then
 839 !              write(50,*) 'i=',i,' j=',j,' rainnc=',rainnc
 840 !           endif
 841 !        enddo
 842 !     enddo
 843 !  endif
 844
 845   END SUBROUTINE fall_flux
 846
 847 !----------------------------------------------------------------
 848    REAL FUNCTION ggamma(X)
 849 !----------------------------------------------------------------
 850    IMPLICIT NONE
 851 !----------------------------------------------------------------
 852       REAL, INTENT(IN   ) :: x
 853       REAL, DIMENSION(8)  :: B
 854       INTEGER             ::j, K1
 855       REAL                ::PF, G1TO2 ,TEMP
 856
 857       DATA B/-.577191652,.988205891,-.897056937,.918206857,  &
 858              -.756704078,.482199394,-.193527818,.035868343/
 859
 860       PF=1.
 861       TEMP=X
 862       DO 10 J=1,200
 863       IF (TEMP .LE. 2) GO TO 20
 864       TEMP=TEMP-1.
 865    10 PF=PF*TEMP
 866   100 FORMAT(//,5X,'module_gsfcgce: INPUT TO GAMMA FUNCTION TOO LARGE, X=',E12.5)
 867       WRITE(wrf_err_message,100)X
 868       CALL wrf_error_fatal(wrf_err_message)
 869    20 G1TO2=1.
 870       TEMP=TEMP - 1.
 871       DO 30 K1=1,8
 872    30 G1TO2=G1TO2 + B(K1)*TEMP**K1
 873       ggamma=PF*G1TO2
 874
 875       END FUNCTION ggamma
 876
 877 !-----------------------------------------------------------------------
 878 !c Correction of negative values
 879    SUBROUTINE negcor ( X, rho, dz8w,                         &
 880                       ims,ime, jms,jme, kms,kme,              & ! memory dims
 881                       itimestep, ics,                         &
 882                       its,ite, jts,jte, kts,kte               ) ! tile   dims
 883 !-----------------------------------------------------------------------
 884   REAL, DIMENSION( ims:ime , kms:kme , jms:jme ),                 &
 885         INTENT(INOUT) ::                                     X
 886   REAL, DIMENSION( ims:ime , kms:kme , jms:jme ),                 &
 887         INTENT(IN   ) ::                              rho, dz8w
 888   integer, INTENT(IN   ) ::                           itimestep, ics
 889
 890 !c Local variables
 891 !  REAL, DIMENSION( kts:kte ) ::  Y1, Y2
 892   REAL   ::   A0, A1, A2
 893
 894   A1=0.
 895   A2=0.
 896   do k=kts,kte
 897      do j=jts,jte
 898         do i=its,ite
 899         A1=A1+max(X(i,k,j), 0.)*rho(i,k,j)*dz8w(i,k,j)
 900         A2=A2+max(-X(i,k,j), 0.)*rho(i,k,j)*dz8w(i,k,j)
 901         enddo
 902      enddo
 903   enddo
 904
 905 !  A1=0.0
 906 !  A2=0.0
 907 !  do k=kts,kte
 908 !     A1=A1+Y1(k)
 909 !     A2=A2+Y2(k)
 910 !  enddo
 911
 912   A0=0.0
 913
 914   if (A1.NE.0.0.and.A1.GT.A2) then
 915      A0=(A1-A2)/A1
 916
 917   if (mod(itimestep,540).eq.0) then
 918      if (ics.eq.1) then
 919         write(61,*) 'kms=',kms,'  kme=',kme,'  kts=',kts,'  kte=',kte
 920         write(61,*) 'jms=',jms,'  jme=',jme,'  jts=',jts,'  jte=',jte
 921         write(61,*) 'ims=',ims,'  ime=',ime,'  its=',its,'  ite=',ite
 922      endif
 923      if (ics.eq.1) then
 924          write(61,*) 'qv timestep=',itimestep
 925          write(61,*) '  A1=',A1,'   A2=',A2,'   A0=',A0
 926      else if (ics.eq.2) then
 927              write(61,*) 'ql timestep=',itimestep
 928              write(61,*) '  A1=',A1,'   A2=',A2,'   A0=',A0
 929      else if (ics.eq.3) then
 930              write(61,*) 'qr timestep=',itimestep
 931              write(61,*) '  A1=',A1,'   A2=',A2,'   A0=',A0
 932      else if (ics.eq.4) then
 933              write(61,*) 'qi timestep=',itimestep
 934              write(61,*) '  A1=',A1,'   A2=',A2,'   A0=',A0
 935      else if (ics.eq.5) then
 936              write(61,*) 'qs timestep=',itimestep
 937              write(61,*) '  A1=',A1,'   A2=',A2,'   A0=',A0
 938      else if (ics.eq.6) then
 939              write(61,*) 'qg timestep=',itimestep
 940              write(61,*) '  A1=',A1,'   A2=',A2,'   A0=',A0
 941      else
 942              write(61,*) 'wrong cloud specieis number'
 943      endif
 944   endif
 945
 946      do k=kts,kte
 947         do j=jts,jte
 948            do i=its,ite
 949            X(i,k,j)=A0*AMAX1(X(i,k,j), 0.0)
 950            enddo
 951         enddo
 952      enddo
 953   endif
 954
 955   END SUBROUTINE negcor
 956
 957  SUBROUTINE consat_s (ihail,itaobraun)
 958
 959 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 960 !                                                                      c
 961 !   Tao, W.-K., and J. Simpson, 1989: Modeling study of a tropical     c
 962 !   squall-type convective line. J. Atmos. Sci., 46, 177-202.          c
 963 !                                                                      c
 964 !   Tao, W.-K., J. Simpson and M. McCumber, 1989: An ice-water         c
 965 !   saturation adjustment. Mon. Wea. Rev., 117, 231-235.               c
 966 !                                                                      c
 967 !   Tao, W.-K., and J. Simpson, 1993: The Goddard Cumulus Ensemble     c
 968 !   Model. Part I: Model description. Terrestrial, Atmospheric and     c
 969 !   Oceanic Sciences, 4, 35-72.                                        c
 970 !                                                                      c
 971 !   Tao, W.-K., J. Simpson, D. Baker, S. Braun, M.-D. Chou, B.         c
 972 !   Ferrier,D. Johnson, A. Khain, S. Lang,  B. Lynn, C.-L. Shie,       c
 973 !   D. Starr, C.-H. Sui, Y. Wang and P. Wetzel, 2003: Microphysics,    c
 974 !   radiation and surface processes in the Goddard Cumulus Ensemble    c
 975 !   (GCE) model, A Special Issue on Non-hydrostatic Mesoscale          c
 976 !   Modeling, Meteorology and Atmospheric Physics, 82, 97-137.         c
 977 !                                                                      c
 978 !   Lang, S., W.-K. Tao, R. Cifelli, W. Olson, J. Halverson, S.        c
 979 !   Rutledge, and J. Simpson, 2007: Improving simulations of           c
 980 !   convective system from TRMM LBA: Easterly and Westerly regimes.    c
 981 !   J. Atmos. Sci., 64, 1141-1164.                                     c
 982 !                                                                      c
 983 !   Coded by Tao (1989-2003), modified by S. Lang (2006/07)            c
 984 !                                                                      c
 985 !   Implemented into WRF  by Roger Shi 2006/2007                       c
 986 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 987
 988 !        itaobraun=0   ! see Tao and Simpson (1993)
 989 !        itaobraun=1   ! see Tao et al. (2003)
 990
 991  integer :: itaobraun
 992  real    :: cn0
 993
 994 !JJS 1/3/2008  vvvvv
 995 !JJS   the following common blocks have been moved to the top of
 996 !JJS   module_mp_gsfcgce_driver_instat.F
 997 !
 998 ! real,   dimension (1:31) ::  a1, a2
 999 ! data a1/.7939e-7,.7841e-6,.3369e-5,.4336e-5,.5285e-5,.3728e-5, &
1000 !         .1852e-5,.2991e-6,.4248e-6,.7434e-6,.1812e-5,.4394e-5,.9145e-5, &
1001 !         .1725e-4,.3348e-4,.1725e-4,.9175e-5,.4412e-5,.2252e-5,.9115e-6, &
1002 !         .4876e-6,.3473e-6,.4758e-6,.6306e-6,.8573e-6,.7868e-6,.7192e-6, &
1003 !         .6513e-6,.5956e-6,.5333e-6,.4834e-6/
1004 ! data a2/.4006,.4831,.5320,.5307,.5319,.5249,.4888,.3894,.4047, &
1005 !         .4318,.4771,.5183,.5463,.5651,.5813,.5655,.5478,.5203,.4906, &
1006 !         .4447,.4126,.3960,.4149,.4320,.4506,.4483,.4460,.4433,.4413, &
1007 !         .4382,.4361/
1008 !JJS 1/3/2008  ^^^^^
1009
1010
1011 !     ******************************************************************
1012 !JJS
1013       al = 2.5e10
1014       cp = 1.004e7
1015       rd1 = 1.e-3
1016       rd2 = 2.2
1017 !JJS
1018       cpi=4.*atan(1.)
1019       cpi2=cpi*cpi
1020       grvt=980.
1021       cd1=6.e-1
1022       cd2=4.*grvt/(3.*cd1)
1023       tca=2.43e3
1024       dwv=.226
1025       dva=1.718e-4
1026       amw=18.016
1027       ars=8.314e7
1028       scv=2.2904487
1029       t0=273.16
1030       t00=238.16
1031       alv=2.5e10
1032       alf=3.336e9
1033       als=2.8336e10
1034       avc=alv/cp
1035       afc=alf/cp
1036       asc=als/cp
1037       rw=4.615e6
1038       cw=4.187e7
1039       ci=2.093e7
1040       c76=7.66
1041       c358=35.86
1042       c172=17.26939
1043       c409=4098.026
1044       c218=21.87456
1045       c580=5807.695
1046       c610=6.1078e3
1047       c149=1.496286e-5
1048       c879=8.794142
1049       c141=1.4144354e7
1050 !***   DEFINE THE COEFFICIENTS USED IN TERMINAL VELOCITY
1051 !***   DEFINE THE DENSITY AND SIZE DISTRIBUTION OF PRECIPITATION
1052 !**********   HAIL OR GRAUPEL PARAMETERS   **********
1053       if(ihail .eq. 1) then
1054          roqg=.9
1055          ag=sqrt(cd2*roqg)
1056          bg=.5
1057          tng=.002
1058       else
1059          roqg=.4
1060          ag=351.2
1061 !        AG=372.3 ! if ice913=1 6/15/02 tao's
1062          bg=.37
1063          tng=.04
1064       endif
1065 !**********         SNOW PARAMETERS        **********
1066 !ccshie 6/15/02 tao's
1067 !      TNS=1.
1068 !      TNS=.08 ! if ice913=1, tao's
1069       tns=.16 ! if ice913=0, tao's
1070       roqs=.1
1071 !      AS=152.93
1072       as=78.63
1073 !      BS=.25
1074       bs=.11
1075 !**********         RAIN PARAMETERS        **********
1076       aw=2115.
1077       bw=.8
1078       roqr=1.
1079       tnw=.08
1080 !*****************************************************************
1081       bgh=.5*bg
1082       bsh=.5*bs
1083       bwh=.5*bw
1084       bgq=.25*bg
1085       bsq=.25*bs
1086       bwq=.25*bw
1087 !**********GAMMA FUNCTION CALCULATIONS*************
1088       ga3b  = gammagce(3.+bw)
1089       ga4b  = gammagce(4.+bw)
1090       ga6b  = gammagce(6.+bw)
1091       ga5bh = gammagce((5.+bw)/2.)
1092       ga3g  = gammagce(3.+bg)
1093       ga4g  = gammagce(4.+bg)
1094       ga5gh = gammagce((5.+bg)/2.)
1095       ga3d  = gammagce(3.+bs)
1096       ga4d  = gammagce(4.+bs)
1097       ga5dh = gammagce((5.+bs)/2.)
1098 !CCCCC        LIN ET AL., 1983 OR LORD ET AL., 1984   CCCCCCCCCCCCCCCCC
1099       ac1=aw
1100 !JJS
1101       ac2=ag
1102       ac3=as
1103 !JJS
1104       bc1=bw
1105       cc1=as
1106       dc1=bs
1107       zrc=(cpi*roqr*tnw)**0.25
1108       zsc=(cpi*roqs*tns)**0.25
1109       zgc=(cpi*roqg*tng)**0.25
1110       vrc=aw*ga4b/(6.*zrc**bw)
1111       vsc=as*ga4d/(6.*zsc**bs)
1112       vgc=ag*ga4g/(6.*zgc**bg)
1113 !     ****************************
1114 !     RN1=1.E-3
1115       rn1=9.4e-15 ! 6/15/02 tao's
1116       bnd1=6.e-4
1117       rn2=1.e-3
1118 !     BND2=1.25E-3
1119 !     BND2=1.5E-3 ! if ice913=1 6/15/02 tao's
1120       bnd2=2.0e-3 ! if ice913=0 6/15/02 tao's
1121       rn3=.25*cpi*tns*cc1*ga3d
1122       esw=1.
1123       rn4=.25*cpi*esw*tns*cc1*ga3d
1124 !     ERI=1.
1125       eri=.1  ! 6/17/02 tao's ice913=0 (not 1)
1126       rn5=.25*cpi*eri*tnw*ac1*ga3b
1127 !     AMI=1./(24.*4.19E-10)
1128       ami=1./(24.*6.e-9) ! 6/15/02 tao's
1129       rn6=cpi2*eri*tnw*ac1*roqr*ga6b*ami
1130 !     ESR=1. ! also if ice913=1 for tao's
1131       esr=.5 ! 6/15/02 for ice913=0 tao's
1132       rn7=cpi2*esr*tnw*tns*roqs
1133       esr=1.
1134       rn8=cpi2*esr*tnw*tns*roqr
1135       rn9=cpi2*tns*tng*roqs
1136       rn10=2.*cpi*tns
1137       rn101=.31*ga5dh*sqrt(cc1)
1138       rn10a=als*als/rw
1139 !JJS
1140        rn10b=alv/tca
1141        rn10c=ars/(dwv*amw)
1142 !JJS
1143       rn11=2.*cpi*tns/alf
1144       rn11a=cw/alf
1145 !     AMI50=1.51e-7
1146       ami50=3.84e-6 ! 6/15/02 tao's
1147 !     AMI40=2.41e-8
1148       ami40=3.08e-8 ! 6/15/02 tao's
1149       eiw=1.
1150 !     UI50=20.
1151       ui50=100. ! 6/15/02 tao's
1152       ri50=2.*5.e-3
1153       cmn=1.05e-15
1154       rn12=cpi*eiw*ui50*ri50**2
1155
1156       do 10 k=1,31
1157          y1=1.-aa2(k)
1158          rn13(k)=aa1(k)*y1/(ami50**y1-ami40**y1)
1159          rn12a(k)=rn13(k)/ami50
1160          rn12b(k)=aa1(k)*ami50**aa2(k)
1161          rn25a(k)=aa1(k)*cmn**aa2(k)
1162    10 continue
1163
1164       egw=1.
1165       rn14=.25*cpi*egw*tng*ga3g*ag
1166       egi=.1
1167       rn15=.25*cpi*egi*tng*ga3g*ag
1168       egi=1.
1169       rn15a=.25*cpi*egi*tng*ga3g*ag
1170       egr=1.
1171       rn16=cpi2*egr*tng*tnw*roqr
1172       rn17=2.*cpi*tng
1173       rn17a=.31*ga5gh*sqrt(ag)
1174       rn17b=cw-ci
1175       rn17c=cw
1176       apri=.66
1177       bpri=1.e-4
1178       bpri=0.5*bpri ! 6/17/02 tao's
1179       rn18=20.*cpi2*bpri*tnw*roqr
1180       rn18a=apri
1181       rn19=2.*cpi*tng/alf
1182       rn19a=.31*ga5gh*sqrt(ag)
1183       rn19b=cw/alf
1184 !
1185        rnn191=.78
1186        rnn192=.31*ga5gh*sqrt(ac2/dva)
1187 !
1188       rn20=2.*cpi*tng
1189       rn20a=als*als/rw
1190       rn20b=.31*ga5gh*sqrt(ag)
1191       bnd3=2.e-3
1192       rn21=1.e3*1.569e-12/0.15
1193       erw=1.
1194       rn22=.25*cpi*erw*ac1*tnw*ga3b
1195       rn23=2.*cpi*tnw
1196       rn23a=.31*ga5bh*sqrt(ac1)
1197       rn23b=alv*alv/rw
1198
1199
1200 !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1201 !cc
1202 !cc        "c0" in routine      "consat" (2d), "consatrh" (3d)
1203 !cc        if ( itaobraun.eq.1 ) --> betah=0.5*beta=-.46*0.5=-0.23;   cn0=1.e-6
1204 !cc        if ( itaobraun.eq.0 ) --> betah=0.5*beta=-.6*0.5=-0.30;    cn0=1.e-8
1205
1206        if (itaobraun .eq. 0) then
1207          cn0=1.e-8
1208          beta=-.6
1209        elseif (itaobraun .eq. 1) then
1210          cn0=1.e-6
1211          beta=-.46
1212        endif
1213 !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1214 !      CN0=1.E-6
1215 !      CN0=1.E-8 ! 6/15/02 tao's
1216 !      BETA=-.46
1217 !      BETA=-.6  ! 6/15/02 tao's
1218
1219       rn25=cn0
1220       rn30a=alv*als*amw/(tca*ars)
1221       rn30b=alv/tca
1222       rn30c=ars/(dwv*amw)
1223       rn31=1.e-17
1224
1225       rn32=4.*51.545e-4
1226 !
1227       rn30=2.*cpi*tng
1228       rnn30a=alv*alv*amw/(tca*ars)
1229 !
1230       rn33=4.*tns
1231        rn331=.65
1232        rn332=.44*sqrt(ac3/dva)*ga5dh
1233 !
1234
1235     return
1236  END SUBROUTINE consat_s
1237
1238  SUBROUTINE saticel_s (dt, ihail, itaobraun, ice2, istatmin, &
1239                        new_ice_sat, id, &
1240                        ptwrf, qvwrf, qlwrf, qrwrf, &
1241                        qiwrf, qswrf, qgwrf, &
1242                        rho_mks, pi_mks, p0_mks,itimestep, &
1243                        refl_10cm, diagflag, do_radar_ref,           & ! GT added for reflectivity calcs
1244                        ids,ide, jds,jde, kds,kde, &
1245                        ims,ime, jms,jme, kms,kme, &
1246                        its,ite, jts,jte, kts,kte  &
1247                            )
1248     IMPLICIT NONE
1249 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1250 !                                                                      c
1251 !   Tao, W.-K., and J. Simpson, 1989: Modeling study of a tropical     c
1252 !   squall-type convective line. J. Atmos. Sci., 46, 177-202.          c
1253 !                                                                      c
1254 !   Tao, W.-K., J. Simpson and M. McCumber, 1989: An ice-water         c
1255 !   saturation adjustment. Mon. Wea. Rev., 117, 231-235.               c
1256 !                                                                      c
1257 !   Tao, W.-K., and J. Simpson, 1993: The Goddard Cumulus Ensemble     c
1258 !   Model. Part I: Model description. Terrestrial, Atmospheric and     c
1259 !   Oceanic Sciences, 4, 35-72.                                        c
1260 !                                                                      c
1261 !   Tao, W.-K., J. Simpson, D. Baker, S. Braun, M.-D. Chou, B.         c
1262 !   Ferrier,D. Johnson, A. Khain, S. Lang,  B. Lynn, C.-L. Shie,       c
1263 !   D. Starr, C.-H. Sui, Y. Wang and P. Wetzel, 2003: Microphysics,    c
1264 !   radiation and surface processes in the Goddard Cumulus Ensemble    c
1265 !   (GCE) model, A Special Issue on Non-hydrostatic Mesoscale          c
1266 !   Modeling, Meteorology and Atmospheric Physics, 82, 97-137.         c
1267 !                                                                      c
1268 !   Lang, S., W.-K. Tao, R. Cifelli, W. Olson, J. Halverson, S.        c
1269 !   Rutledge, and J. Simpson, 2007: Improving simulations of           c
1270 !   convective system from TRMM LBA: Easterly and Westerly regimes.    c
1271 !   J. Atmos. Sci., 64, 1141-1164.                                     c
1272 !                                                                      c
1273 !   Tao, W.-K., J. J. Shi,  S. Lang, C. Peters-Lidard, A. Hou, S.      c
1274 !   Braun, and J. Simpson, 2007: New, improved bulk-microphysical      c
1275 !   schemes for studying precipitation processes in WRF. Part I:       c
1276 !   Comparisons with other schemes. to appear on Mon. Wea. Rev.        C
1277 !                                                                      c
1278 !   Coded by Tao (1989-2003), modified by S. Lang (2006/07)            c
1279 !                                                                      c
1280 !   Implemented into WRF  by Roger Shi 2006/2007                       c
1281 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1282 !
1283 !      COMPUTE ICE PHASE MICROPHYSICS AND SATURATION PROCESSES
1284 !
1285   integer,    parameter ::  nt=2880, nt2=2*nt
1286
1287 !cc   using scott braun's way for pint, pidep computations
1288   integer  ::   itaobraun,ice2,ihail,new_ice_sat,id,istatmin
1289   integer  ::   itimestep
1290   real     ::   tairccri, cn0, dt
1291 !cc
1292
1293 !JJS      common/bxyz/ n,isec,nran,kt1,kt2
1294 !JJS      common/option/ lipps,ijkadv,istatmin,iwater,itoga,imlifting,lin,
1295 !JJS     1   irf,iadvh,irfg,ismg,id
1296
1297 !JJS      common/timestat/ ndt_stat,idq
1298 !JJS      common/iice/ new_ice_sat
1299 !JJS      common/bt/ dt,d2t,rijl2,dts,f5,rd1,rd2,bound,al,cp,ra,ck,ce,eps,
1300 !JJS     1    psfc,fcor,sec,aminut,rdt
1301
1302 !JJS   the following common blocks have been moved to the top of
1303 !JJS   module_mp_gsfcgce_driver_instat.F
1304
1305 !      common/bt/ rd1,rd2,al,cp
1306 !
1307 !
1308 !      common/bterv/ zrc,zgc,zsc,vrc,vgc,vsc
1309 !      common/size/ tnw,tns,tng,roqs,roqg,roqr
1310 !      common/cont/ c38,c358,c610,c149,c879,c172,c409,c76,c218,c580,c141
1311 !        common/b3cs/ ag,bg,as,bs,aw,bw,bgh,bgq,bsh,bsq,bwh,bwq
1312 !      common/bsnw/ alv,alf,als,t0,t00,avc,afc,asc,rn1,bnd1,rn2,bnd2, &
1313 !         rn3,rn4,rn5,rn6,rn7,rn8,rn9,rn10,rn101,rn10a,rn11,rn11a, &
1314 !         rn12,rn12a(31),rn12b(31),rn13(31),rn14,rn15,rn15a,rn16,rn17, &
1315 !         rn17a,rn17b,rn17c,rn18,rn18a,rn19,rn19a,rn19b,rn20,rn20a,rn20b, &
1316 !         bnd3,rn21,rn22,rn23,rn23a,rn23b,rn25,rn25a(31),rn30a,rn30b, &
1317 !         rn30c,rn31,beta,rn32
1318 !      common/rsnw1/ rn10b,rn10c,rnn191,rnn192,rn30,rnn30a,rn33,rn331, &
1319 !                    rn332
1320 !JJS
1321
1322   integer ids,ide,jds,jde,kds,kde
1323   integer ims,ime,jms,jme,kms,kme
1324   integer its,ite,jts,jte,kts,kte
1325   integer i,j,k, kp
1326
1327   real :: a0 ,a1 ,a2 ,afcp ,alvr ,ami100 ,ami40 ,ami50 ,ascp ,avcp ,betah &
1328    ,bg3 ,bgh5 ,bs3 ,bs6 ,bsh5 ,bw3 ,bw6 ,bwh5 ,cmin ,cmin1 ,cmin2 ,cp409 &
1329    ,cp580 ,cs580 ,cv409 ,d2t ,del ,dwvp ,ee1 ,ee2 ,f00 ,f2 ,f3 ,ft ,fv0 ,fvs &
1330    ,pi0 ,pir ,pr0 ,qb0 ,r00 ,r0s ,r101f ,r10ar ,r10t ,r11at ,r11rt ,r12r ,r14f &
1331    ,r14r ,r15af ,r15ar ,r15f ,r15r ,r16r ,r17aq ,r17as ,r17r ,r18r ,r19aq ,r19as &
1332    ,r19bt ,r19rt ,r20bq ,r20bs ,r20t ,r22f ,r23af ,r23br ,r23t ,r25a ,r25rt ,r2ice &
1333    ,r31r ,r32rt ,r3f ,r4f ,r5f ,r6f ,r7r ,r8r ,r9r ,r_nci ,rft ,rijl2 ,rp0 ,rr0 &
1334    ,rrq ,rrs ,rt0 ,scc ,sccc ,sddd ,see ,seee ,sfff ,smmm ,ssss ,tb0 ,temp ,ucog &
1335    ,ucor ,ucos ,uwet ,vgcf ,vgcr ,vrcf ,vscf ,zgr ,zrr ,zsr
1336
1337
1338   real, dimension (its:ite,jts:jte,kts:kte) ::  fv
1339   real, dimension (its:ite,jts:jte,kts:kte) ::  dpt, dqv
1340   real, dimension (its:ite,jts:jte,kts:kte) ::  qcl, qrn,      &
1341                                                 qci, qcs, qcg
1342 !JJS 10/16/06    vvvv
1343 !      real dpt1(ims:ime,jms:jme,kms:kme)
1344 !      real dqv1(ims:ime,jms:jme,kms:kme),
1345 !     1             qcl1(ims:ime,jms:jme,kms:kme)
1346 !      real qrn1(ims:ime,jms:jme,kms:kme),
1347 !     1             qci1(ims:ime,jms:jme,kms:kme)
1348 !      real qcs1(ims:ime,jms:jme,kms:kme),
1349 !     1             qcg1(ims:ime,jms:jme,kms:kme)
1350 !JJS 10/16/06    ^^^^
1351
1352 !JJS
1353
1354   real, dimension (ims:ime, kms:kme, jms:jme) ::  ptwrf, qvwrf
1355   real, dimension (ims:ime, kms:kme, jms:jme) ::  qlwrf, qrwrf,        &
1356                                                   qiwrf, qswrf, qgwrf
1357 !JJS 10/16/06    vvvv
1358 !      real ptwrfold(ims:ime, kms:kme, jms:jme)
1359 !      real qvwrfold(ims:ime, kms:kme, jms:jme),
1360 !     1             qlwrfold(ims:ime, kms:kme, jms:jme)
1361 !      real qrwrfold(ims:ime, kms:kme, jms:jme),
1362 !     1             qiwrfold(ims:ime, kms:kme, jms:jme)
1363 !      real qswrfold(ims:ime, kms:kme, jms:jme),
1364 !     1             qgwrfold(ims:ime, kms:kme, jms:jme)
1365 !JJS 10/16/06    ^^^^
1366
1367 !JJS in MKS
1368   real, dimension (ims:ime, kms:kme, jms:jme) ::  rho_mks
1369   real, dimension (ims:ime, kms:kme, jms:jme) ::  pi_mks
1370   real, dimension (ims:ime, kms:kme, jms:jme) ::  p0_mks
1371 !JJS
1372 !  real, dimension (its:ite,jts:jte,kts:kte) ::  ww1
1373 !  real, dimension (its:ite,jts:jte,kts:kte) ::  rsw
1374 !  real, dimension (its:ite,jts:jte,kts:kte) ::  rlw
1375
1376 !JJS      COMMON /BADV/
1377   real, dimension (its:ite,jts:jte) ::        &
1378            vg,      zg,       &
1379            ps,      pg,       &
1380           prn,     psn,       &
1381         pwacs,   wgacr,       &
1382         pidep,    pint,       &
1383           qsi,     ssi,       &
1384           esi,     esw,       &
1385           qsw,      pr,       &
1386           ssw,   pihom,       &
1387          pidw,   pimlt,       &
1388         psaut,   qracs,       &
1389         psaci,   psacw,       &
1390         qsacw,   praci,       &
1391         pmlts,   pmltg,       &
1392         asss,       y1,    y2
1393 !JJS       Y2(its:ite,jts:jte),    DDE(NB)
1394
1395 !JJS      COMMON/BSAT/
1396   real, dimension (its:ite,jts:jte) ::        &
1397         praut,   pracw,       &
1398          psfw,    psfi,       &
1399         dgacs,   dgacw,       &
1400         dgaci,   dgacr,       &
1401         pgacs,   wgacs,       &
1402         qgacw,   wgaci,       &
1403         qgacr,   pgwet,       &
1404         pgaut,   pracs,       &
1405         psacr,   qsacr,       &
1406          pgfr,   psmlt,       &
1407         pgmlt,   psdep,       &
1408         pgdep,   piacr,       &
1409            y5,     scv,       &
1410           tca,     dwv,       &
1411           egs,      y3,       &
1412            y4,     ddb
1413
1414 !JJS      COMMON/BSAT1/
1415   real, dimension (its:ite,jts:jte) ::        &
1416            pt,      qv,       &
1417            qc,      qr,       &
1418            qi,      qs,       &
1419            qg,    tair,       &
1420         tairc,   rtair,       &
1421           dep,      dd,       &
1422           dd1,     qvs,       &
1423            dm,      rq,       &
1424         rsub1,     col,       &
1425           cnd,     ern,       &
1426          dlt1,    dlt2,       &
1427          dlt3,    dlt4,       &
1428            zr,      vr,       &
1429            zs,      vs,       &
1430                  pssub,       &
1431         pgsub,     dda
1432
1433 !JJS      COMMON/B5/
1434   real, dimension (its:ite,jts:jte,kts:kte) ::  rho
1435   real, dimension (kts:kte) ::                 &
1436            tb,      qb,    rho1,              &
1437            ta,      qa,     ta1,     qa1,     &
1438          coef,      z1,      z2,      z3,     &
1439            am,     am1,      ub,      vb,     &
1440            wb,     ub1,     vb1,    rrho,     &
1441         rrho1,     wbx
1442
1443 !JJS      COMMON/B6/
1444   real, dimension (its:ite,jts:jte,kts:kte) ::  p0, pi, f0
1445   real, dimension (kts:kte) ::    &
1446            fd,      fe,        &
1447            st,      sv,        &
1448            sq,      sc,        &
1449            se,     sqa
1450
1451 !JJS      COMMON/BRH1/
1452   real, dimension (kts:kte) ::    &
1453          srro,    qrro,    sqc,    sqr,    &
1454           sqi,     sqs,    sqg,   stqc,    &
1455          stqr,    stqi,   stqs,   stqg
1456   real, dimension (nt) ::    &
1457           tqc,     tqr,    tqi,    tqs,    tqg
1458
1459 !JJS     common/bls/ y0(nx,ny),ts0new(nx,ny),qss0new(nx,ny)
1460
1461 !JJS      COMMON/BLS/
1462   real, dimension (ims:ime,jms:jme) ::     &
1463            y0,     ts0,   qss0
1464
1465 !JJS      COMMON/BI/ IT(its:ite,jts:jte), ICS(its:ite,jts:jte,4)
1466   integer, dimension (its:ite,jts:jte) ::        it
1467   integer, dimension (its:ite,jts:jte, 4) ::    ics
1468
1469   integer :: i24h
1470   integer :: iwarm
1471   real :: r2is, r2ig
1472
1473
1474 !JJS      COMMON/MICRO/
1475 !  real, dimension (ims:ime,kms:kme,jms:jme)  ::    dbz
1476
1477 !23456789012345678901234567890123456789012345678901234567890123456789012
1478
1479 !
1480 !JJS 1/3/2008  vvvvv
1481 !JJS   the following common blocks have been moved to the top of
1482 !JJS   module_mp_gsfcgce_driver.F
1483
1484 !  real, dimension (31)   ::      aa1,  aa2
1485 !  data aa1/.7939e-7, .7841e-6, .3369e-5, .4336e-5, .5285e-5,     &
1486 !           .3728e-5, .1852e-5, .2991e-6, .4248e-6, .7434e-6,     &
1487 !           .1812e-5, .4394e-5, .9145e-5, .1725e-4, .3348e-4,     &
1488 !           .1725e-4, .9175e-5, .4412e-5, .2252e-5, .9115e-6,     &
1489 !           .4876e-6, .3473e-6, .4758e-6, .6306e-6, .8573e-6,     &
1490 !           .7868e-6, .7192e-6, .6513e-6, .5956e-6, .5333e-6,     &
1491 !           .4834e-6/
1492 !  data aa2/.4006, .4831, .5320, .5307, .5319,      &
1493 !           .5249, .4888, .3894, .4047, .4318,      &
1494 !           .4771, .5183, .5463, .5651, .5813,      &
1495 !           .5655, .5478, .5203, .4906, .4447,      &
1496 !           .4126, .3960, .4149, .4320, .4506,      &
1497 !           .4483, .4460, .4433, .4413, .4382,      &
1498 !           .4361/
1499
1500 !JJS 1/3/2008  ^^^^^
1501
1502 !+---+-----------------------------------------------------------------+
1503   REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT):: refl_10cm  ! GT
1504
1505
1506       REAL, DIMENSION(kts:kte):: qv1d, t1d, p1d, qr1d, qs1d, qg1d, dBZ
1507       LOGICAL, OPTIONAL, INTENT(IN) :: diagflag
1508       INTEGER, OPTIONAL, INTENT(IN) :: do_radar_ref
1509 !+---+-----------------------------------------------------------------+
1510 !
1511 !jm 20090220      save
1512
1513 !    i24h=nint(86400./dt)
1514 !    if (mod(itimestep,i24h).eq.1) then
1515 !       write(6, *) 'ims=', ims, '  ime=', ime
1516 !       write(6, *) 'jms=', jms, '  jme=', jme
1517 !       write(6, *) 'kms=', kms, '  kme=', kme
1518 !       write(6, *) 'its=', its, '  ite=', ite
1519 !       write(6, *) 'jts=', jts, '  jte=', jte
1520 !       write(6, *) 'kts=', kts, '  kte=', kte
1521 !       write(6, *) '    ihail=', ihail
1522 !       write(6, *) 'itaobraun=',itaobraun
1523 !       write(6, *) '    ice2=', ICE2
1524 !       write(6, *) 'istatmin=',istatmin
1525 !       write(6, *) 'new_ice_sat=', new_ice_sat
1526 !       write(6, *) 'id=', id
1527 !       write(6, *) 'dt=', dt
1528 !    endif
1529
1530 !JJS  convert from mks to cgs, and move from WRF grid to GCE grid
1531       do k=kts,kte
1532          do j=jts,jte
1533          do i=its,ite
1534          rho(i,j,k)=rho_mks(i,k,j)*0.001
1535          p0(i,j,k)=p0_mks(i,k,j)*10.0
1536          pi(i,j,k)=pi_mks(i,k,j)
1537          dpt(i,j,k)=ptwrf(i,k,j)
1538          dqv(i,j,k)=qvwrf(i,k,j)
1539          qcl(i,j,k)=qlwrf(i,k,j)
1540          qrn(i,j,k)=qrwrf(i,k,j)
1541          qci(i,j,k)=qiwrf(i,k,j)
1542          qcs(i,j,k)=qswrf(i,k,j)
1543          qcg(i,j,k)=qgwrf(i,k,j)
1544 !JJS 10/16/06    vvvv
1545 !         dpt1(i,j,k)=ptwrfold(i,k,j)
1546 !         dqv1(i,j,k)=qvwrfold(i,k,j)
1547 !         qcl1(i,j,k)=qlwrfold(i,k,j)
1548 !         qrn1(i,j,k)=qrwrfold(i,k,j)
1549 !         qci1(i,j,k)=qiwrfold(i,k,j)
1550 !         qcs1(i,j,k)=qswrfold(i,k,j)
1551 !         qcg1(i,j,k)=qgwrfold(i,k,j)
1552 !JJS 10/16/06     ^^^^
1553          enddo !i
1554          enddo !j
1555       enddo !k
1556
1557       do k=kts,kte
1558          do j=jts,jte
1559          do i=its,ite
1560          fv(i,j,k)=sqrt(rho(i,j,2)/rho(i,j,k))
1561          enddo !i
1562          enddo !j
1563       enddo !k
1564 !JJS
1565
1566 !
1567 !     ******   THREE CLASSES OF ICE-PHASE   (LIN ET AL, 1983)  *********
1568
1569 !JJS       D22T=D2T
1570 !JJS       IF(IJKADV .EQ. 0) THEN
1571 !JJS         D2T=D2T
1572 !JJS       ELSE
1573          d2t=dt
1574 !JJS       ENDIF
1575 !
1576 !       itaobraun=0 ! original pint and pidep & see Tao and Simpson 1993
1577         itaobraun=1 ! see Tao et al. (2003)
1578 !
1579        if ( itaobraun.eq.0 ) then
1580           cn0=1.e-8
1581 !c        beta=-.6
1582        elseif ( itaobraun.eq.1 ) then
1583           cn0=1.e-6
1584 !         cn0=1.e-8  ! special
1585 !c        beta=-.46
1586        endif
1587 !C  TAO 2007 START
1588 !   ICE2=0 ! default, 3ice with loud ice, snow and graupel
1589 !              r2is=1., r2ig=1.
1590 !   ICE2=1 ! 2ice with cloud ice and snow (no graupel) - r2iceg=1, r2ice=0.
1591 !              r2is=1., r2ig=0.
1592 !   ICE2=2 ! 2ice with cloud ice and graupel (no snow) - r2ice=1, r2iceg=0.
1593 !              r2is=0., r2ig=1.
1594 !c
1595 !        r2ice=1.
1596 !        r2iceg=1.
1597          r2ig=1.
1598          r2is=1.
1599           if (ice2 .eq. 1) then
1600 !              r2ice=0.
1601 !              r2iceg=1.
1602               r2ig=0.
1603               r2is=1.
1604           endif
1605           if (ice2 .eq. 2) then
1606 !              r2ice=1.
1607 !              r2iceg=0.
1608               r2ig=1.
1609               r2is=0.
1610           endif
1611 !C  TAO 2007 END
1612
1613 !JJS  10/7/2008
1614 !   ICE2=3 ! no ice, warm rain only
1615     iwarm = 0
1616     if (ice2 .eq. 3 ) iwarm = 1
1617
1618
1619
1620       cmin=1.e-19
1621       cmin1=1.e-20
1622       cmin2=1.e-12
1623       ucor=3071.29/tnw**0.75
1624       ucos=687.97*roqs**0.25/tns**0.75
1625       ucog=687.97*roqg**0.25/tng**0.75
1626       uwet=4.464**0.95
1627
1628       rijl2 = 1. / (ide-ids) / (jde-jds)
1629
1630 !JJScap $doacross local(j,i)
1631
1632 !JJS      DO 1 J=1,JMAX
1633 !JJS      DO 1 I=1,IMAX
1634        do j=jts,jte
1635           do i=its,ite
1636           it(i,j)=1
1637           enddo
1638        enddo
1639
1640       f2=rd1*d2t
1641       f3=rd2*d2t
1642
1643       ft=dt/d2t
1644       rft=rijl2*ft
1645       a0=.5*istatmin*rijl2
1646       rt0=1./(t0-t00)
1647       bw3=bw+3.
1648       bs3=bs+3.
1649       bg3=bg+3.
1650       bsh5=2.5+bsh
1651       bgh5=2.5+bgh
1652       bwh5=2.5+bwh
1653       bw6=bw+6.
1654       bs6=bs+6.
1655       betah=.5*beta
1656       r10t=rn10*d2t
1657       r11at=rn11a*d2t
1658       r19bt=rn19b*d2t
1659       r20t=-rn20*d2t
1660       r23t=-rn23*d2t
1661       r25a=rn25
1662
1663 !     ami50 for use in PINT
1664        ami50=3.76e-8
1665        ami100=1.51e-7
1666        ami40=2.41e-8
1667
1668 !C    ******************************************************************
1669
1670 !JJS      DO 1000 K=2,kles
1671       do 1000 k=kts,kte
1672          kp=k+1
1673 !JJS         tb0=ta1(k)
1674 !JJS         qb0=qa1(k)
1675          tb0=0.
1676          qb0=0.
1677
1678       do 2000 j=jts,jte
1679          do 2000 i=its,ite
1680
1681          rp0=3.799052e3/p0(i,j,k)
1682          pi0=pi(i,j,k)
1683          pir=1./(pi(i,j,k))
1684          pr0=1./p0(i,j,k)
1685          r00=rho(i,j,k)
1686          r0s=sqrt(rho(i,j,k))
1687 !JJS         RR0=RRHO(K)
1688          rr0=1./rho(i,j,k)
1689 !JJS         RRS=SRRO(K)
1690          rrs=sqrt(rr0)
1691 !JJS         RRQ=QRRO(K)
1692          rrq=sqrt(rrs)
1693          f0(i,j,k)=al/cp/pi(i,j,k)
1694          f00=f0(i,j,k)
1695          fv0=fv(i,j,k)
1696          fvs=sqrt(fv(i,j,k))
1697          zrr=1.e5*zrc*rrq
1698          zsr=1.e5*zsc*rrq
1699          zgr=1.e5*zgc*rrq
1700          cp409=c409*pi0
1701          cv409=c409*avc
1702          cp580=c580*pi0
1703          cs580=c580*asc
1704          alvr=r00*alv
1705          afcp=afc*pir
1706          avcp=avc*pir
1707          ascp=asc*pir
1708          vrcf=vrc*fv0
1709          vscf=vsc*fv0
1710          vgcf=vgc*fv0
1711          vgcr=vgc*rrs
1712          dwvp=c879*pr0
1713          r3f=rn3*fv0
1714          r4f=rn4*fv0
1715          r5f=rn5*fv0
1716          r6f=rn6*fv0
1717          r7r=rn7*rr0
1718          r8r=rn8*rr0
1719          r9r=rn9*rr0
1720          r101f=rn101*fvs
1721          r10ar=rn10a*r00
1722          r11rt=rn11*rr0*d2t
1723          r12r=rn12*r00
1724          r14r=rn14*rrs
1725          r14f=rn14*fv0
1726          r15r=rn15*rrs
1727          r15ar=rn15a*rrs
1728          r15f=rn15*fv0
1729          r15af=rn15a*fv0
1730          r16r=rn16*rr0
1731          r17r=rn17*rr0
1732          r17aq=rn17a*rrq
1733          r17as=rn17a*fvs
1734          r18r=rn18*rr0
1735          r19rt=rn19*rr0*d2t
1736          r19aq=rn19a*rrq
1737          r19as=rn19a*fvs
1738          r20bq=rn20b*rrq
1739          r20bs=rn20b*fvs
1740          r22f=rn22*fv0
1741          r23af=rn23a*fvs
1742          r23br=rn23b*r00
1743          r25rt=rn25*rr0*d2t
1744          r31r=rn31*rr0
1745          r32rt=rn32*d2t*rrs
1746
1747 !JJS       DO 100 J=3,JLES
1748 !JJS       DO 100 I=3,ILES
1749         pt(i,j)=dpt(i,j,k)
1750         qv(i,j)=dqv(i,j,k)
1751         qc(i,j)=qcl(i,j,k)
1752         qr(i,j)=qrn(i,j,k)
1753         qi(i,j)=qci(i,j,k)
1754         qs(i,j)=qcs(i,j,k)
1755         qg(i,j)=qcg(i,j,k)
1756 !        IF (QV(I,J)+QB0 .LE. 0.) QV(I,J)=-QB0
1757          if (qc(i,j) .le.  cmin1) qc(i,j)=0.0
1758          if (qr(i,j) .le.  cmin1) qr(i,j)=0.0
1759          if (qi(i,j) .le.  cmin1) qi(i,j)=0.0
1760          if (qs(i,j) .le.  cmin1) qs(i,j)=0.0
1761          if (qg(i,j) .le.  cmin1) qg(i,j)=0.0
1762         tair(i,j)=(pt(i,j)+tb0)*pi0
1763         tairc(i,j)=tair(i,j)-t0
1764          zr(i,j)=zrr
1765          zs(i,j)=zsr
1766          zg(i,j)=zgr
1767          vr(i,j)=0.0
1768          vs(i,j)=0.0
1769          vg(i,j)=0.0
1770
1771 !JJS 10/7/2008     vvvvv
1772     IF (IWARM .EQ. 1) THEN
1773 !JJS   for calculating processes related to warm rain only
1774                 qi(i,j)=0.0
1775                 qs(i,j)=0.0
1776                 qg(i,j)=0.0
1777                 dep(i,j)=0.
1778                 pint(i,j)=0.
1779                 psdep(i,j)=0.
1780                 pgdep(i,j)=0.
1781                 dd1(i,j)=0.
1782                 pgsub(i,j)=0.
1783                 psmlt(i,j)=0.
1784                 pgmlt(i,j)=0.
1785                 pimlt(i,j)=0.
1786                 psacw(i,j)=0.
1787                 piacr(i,j)=0.
1788                 psfw(i,j)=0.
1789                 pgfr(i,j)=0.
1790                 dgacw(i,j)=0.
1791                 dgacr(i,j)=0.
1792                 psacr(i,j)=0.
1793                 wgacr(i,j)=0.
1794                 pihom(i,j)=0.
1795                 pidw(i,j)=0.
1796
1797                 if (qr(i,j) .gt. cmin1) then
1798                    dd(i,j)=r00*qr(i,j)
1799                    y1(i,j)=dd(i,j)**.25
1800                    zr(i,j)=zrc/y1(i,j)
1801                    vr(i,j)=max(vrcf*dd(i,j)**bwq, 0.)
1802                 endif
1803
1804 !* 21 * PRAUT   AUTOCONVERSION OF QC TO QR                        **21**
1805 !* 22 * PRACW : ACCRETION OF QC BY QR                             **22**
1806                 pracw(i,j)=0.
1807                 praut(i,j)=0.0
1808                 pracw(i,j)=r22f*qc(i,j)/zr(i,j)**bw3
1809                 y1(i,j)=qc(i,j)-bnd3
1810                 if (y1(i,j).gt.0.0) then
1811                     praut(i,j)=r00*y1(i,j)*y1(i,j)/(1.2e-4+rn21/y1(i,j))
1812                  endif
1813
1814 !C********   HANDLING THE NEGATIVE CLOUD WATER (QC)    ******************
1815                  Y1(I,J)=QC(I,J)/D2T
1816                  PRAUT(I,J)=MIN(Y1(I,J), PRAUT(I,J))
1817                  PRACW(I,J)=MIN(Y1(I,J), PRACW(I,J))
1818                  Y1(I,J)=(PRAUT(I,J)+PRACW(I,J))*D2T
1819
1820                if (qc(i,j) .lt. y1(i,j) .and. y1(i,j) .ge. cmin2) then
1821                    y2(i,j)=qc(i,j)/(y1(i,j)+cmin2)
1822                    praut(i,j)=praut(i,j)*y2(i,j)
1823                    pracw(i,j)=pracw(i,j)*y2(i,j)
1824                    qc(i,j)=0.0
1825                else
1826                   qc(i,j)=qc(i,j)-y1(i,j)
1827                endif
1828
1829                PR(I,J)=(PRAUT(I,J)+PRACW(I,J))*D2T
1830                QR(I,J)=QR(I,J)+PR(I,J)
1831
1832 !*****   TAO ET AL (1989) SATURATION TECHNIQUE  ***********************
1833
1834            cnd(i,j)=0.0
1835            tair(i,j)=(pt(i,j)+tb0)*pi0
1836               y1(i,j)=1./(tair(i,j)-c358)
1837               qsw(i,j)=rp0*exp(c172-c409*y1(i,j))
1838               dd(i,j)=cp409*y1(i,j)*y1(i,j)
1839               dm(i,j)=qv(i,j)+qb0-qsw(i,j)
1840               cnd(i,j)=dm(i,j)/(1.+avcp*dd(i,j)*qsw(i,j))
1841 !c    ******   condensation or evaporation of qc  ******
1842               cnd(i,j)=max(-qc(i,j), cnd(i,j))
1843                          pt(i,j)=pt(i,j)+avcp*cnd(i,j)
1844              qv(i,j)=qv(i,j)-cnd(i,j)
1845                          qc(i,j)=qc(i,j)+cnd(i,j)
1846
1847 !C     ******   EVAPORATION   ******
1848 !* 23 * ERN : EVAPORATION OF QR (SUBSATURATION)                   **23**
1849             ern(i,j)=0.0
1850
1851             if(qr(i,j).gt.0.0) then
1852                tair(i,j)=(pt(i,j)+tb0)*pi0
1853                rtair(i,j)=1./(tair(i,j)-c358)
1854                qsw(i,j)=rp0*exp(c172-c409*rtair(i,j))
1855                ssw(i,j)=(qv(i,j)+qb0)/qsw(i,j)-1.0
1856                dm(i,j)=qv(i,j)+qb0-qsw(i,j)
1857                rsub1(i,j)=cv409*qsw(i,j)*rtair(i,j)*rtair(i,j)
1858                dd1(i,j)=max(-dm(i,j)/(1.+rsub1(i,j)), 0.0)
1859                y1(i,j)=.78/zr(i,j)**2+r23af*scv(i,j)/zr(i,j)**bwh5
1860                y2(i,j)=r23br/(tca(i,j)*tair(i,j)**2)+1./(dwv(i,j) &
1861                        *qsw(i,j))
1862 !cccc
1863                ern(i,j)=r23t*ssw(i,j)*y1(i,j)/y2(i,j)
1864                ern(i,j)=min(dd1(i,j),qr(i,j),max(ern(i,j),0.))
1865                pt(i,j)=pt(i,j)-avcp*ern(i,j)
1866                qv(i,j)=qv(i,j)+ern(i,j)
1867                qr(i,j)=qr(i,j)-ern(i,j)
1868             endif
1869
1870        ELSE       ! part of if (iwarm.eq.1) then
1871 !JJS 10/7/2008     ^^^^^
1872
1873 !JJS   for calculating processes related to both ice and warm rain
1874
1875 !     ***   COMPUTE ZR,ZS,ZG,VR,VS,VG      *****************************
1876
1877             if (qr(i,j) .gt. cmin1) then
1878                dd(i,j)=r00*qr(i,j)
1879                y1(i,j)=dd(i,j)**.25
1880                zr(i,j)=zrc/y1(i,j)
1881                vr(i,j)=max(vrcf*dd(i,j)**bwq, 0.)
1882             endif
1883
1884             if (qs(i,j) .gt. cmin1) then
1885                dd(i,j)=r00*qs(i,j)
1886                y1(i,j)=dd(i,j)**.25
1887                zs(i,j)=zsc/y1(i,j)
1888                vs(i,j)=max(vscf*dd(i,j)**bsq, 0.)
1889             endif
1890
1891             if (qg(i,j) .gt. cmin1) then
1892                dd(i,j)=r00*qg(i,j)
1893                y1(i,j)=dd(i,j)**.25
1894                zg(i,j)=zgc/y1(i,j)
1895                if(ihail .eq. 1) then
1896                   vg(i,j)=max(vgcr*dd(i,j)**bgq, 0.)
1897                else
1898                   vg(i,j)=max(vgcf*dd(i,j)**bgq, 0.)
1899                endif
1900             endif
1901
1902             if (qr(i,j) .le. cmin2) vr(i,j)=0.0
1903             if (qs(i,j) .le. cmin2) vs(i,j)=0.0
1904             if (qg(i,j) .le. cmin2) vg(i,j)=0.0
1905
1906 !     ******************************************************************
1907 !     ***   Y1 : DYNAMIC VISCOSITY OF AIR (U)
1908 !     ***   DWV : DIFFUSIVITY OF WATER VAPOR IN AIR (PI)
1909 !     ***   TCA : THERMAL CONDUCTIVITY OF AIR (KA)
1910 !     ***   Y2 : KINETIC VISCOSITY (V)
1911
1912             y1(i,j)=c149*tair(i,j)**1.5/(tair(i,j)+120.)
1913             dwv(i,j)=dwvp*tair(i,j)**1.81
1914             tca(i,j)=c141*y1(i,j)
1915             scv(i,j)=1./((rr0*y1(i,j))**.1666667*dwv(i,j)**.3333333)
1916 !JJS  100    CONTINUE
1917
1918 !*  1 * PSAUT : AUTOCONVERSION OF QI TO QS                        ***1**
1919 !*  3 * PSACI : ACCRETION OF QI TO QS                             ***3**
1920 !*  4 * PSACW : ACCRETION OF QC BY QS (RIMING) (QSACW FOR PSMLT)  ***4**
1921 !*  5 * PRACI : ACCRETION OF QI BY QR                             ***5**
1922 !*  6 * PIACR : ACCRETION OF QR OR QG BY QI                       ***6**
1923
1924 !JJS         DO 125 J=3,JLES
1925 !JJS         DO 125 I=3,ILES
1926             psaut(i,j)=0.0
1927             psaci(i,j)=0.0
1928             praci(i,j)=0.0
1929             piacr(i,j)=0.0
1930             psacw(i,j)=0.0
1931             qsacw(i,j)=0.0
1932             dd(i,j)=1./zs(i,j)**bs3
1933
1934             if (tair(i,j).lt.t0) then
1935                esi(i,j)=exp(.025*tairc(i,j))
1936                psaut(i,j)=r2is*max(rn1*esi(i,j)*(qi(i,j)-bnd1) ,0.0)
1937                psaci(i,j)=r2is*r3f*esi(i,j)*qi(i,j)*dd(i,j)
1938 !JJS 3/30/06
1939 !    to cut water to snow accretion by half
1940 !               PSACW(I,J)=R4F*QC(I,J)*DD(I,J)
1941                psacw(i,j)=r2is*0.5*r4f*qc(i,j)*dd(i,j)
1942 !JJS 3/30/06
1943                praci(i,j)=r2is*r5f*qi(i,j)/zr(i,j)**bw3
1944                piacr(i,j)=r2is*r6f*qi(i,j)*(zr(i,j)**(-bw6))
1945 !JJS               PIACR(I,J)=R6F*QI(I,J)/ZR(I,J)**BW6
1946             else
1947                qsacw(i,j)=r2is*r4f*qc(i,j)*dd(i,j)
1948             endif
1949
1950 !* 21 * PRAUT   AUTOCONVERSION OF QC TO QR                        **21**
1951 !* 22 * PRACW : ACCRETION OF QC BY QR                             **22**
1952
1953             pracw(i,j)=r22f*qc(i,j)/zr(i,j)**bw3
1954             praut(i,j)=0.0
1955             y1(i,j)=qc(i,j)-bnd3
1956             if (y1(i,j).gt.0.0) then
1957                praut(i,j)=r00*y1(i,j)*y1(i,j)/(1.2e-4+rn21/y1(i,j))
1958             endif
1959
1960 !* 12 * PSFW : BERGERON PROCESSES FOR QS (KOENING, 1971)          **12**
1961 !* 13 * PSFI : BERGERON PROCESSES FOR QS                          **13**
1962
1963             psfw(i,j)=0.0
1964             psfi(i,j)=0.0
1965             pidep(i,j)=0.0
1966
1967             if(tair(i,j).lt.t0.and.qi(i,j).gt.cmin) then
1968                y1(i,j)=max( min(tairc(i,j), -1.), -31.)
1969                it(i,j)=int(abs(y1(i,j)))
1970                y1(i,j)=rn12a(it(i,j))
1971                y2(i,j)=rn12b(it(i,j))
1972 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1973           psfw(i,j)=r2is* &
1974                     max(d2t*y1(i,j)*(y2(i,j)+r12r*qc(i,j))*qi(i,j),0.0)
1975                rtair(i,j)=1./(tair(i,j)-c76)
1976                y2(i,j)=exp(c218-c580*rtair(i,j))
1977                qsi(i,j)=rp0*y2(i,j)
1978                esi(i,j)=c610*y2(i,j)
1979                ssi(i,j)=(qv(i,j)+qb0)/qsi(i,j)-1.
1980                r_nci=min(1.e-6*exp(-.46*tairc(i,j)),1.)
1981 !              R_NCI=min(1.e-8*EXP(-.6*TAIRC(I,J)),1.) ! use Tao's
1982                dm(i,j)=max( (qv(i,j)+qb0-qsi(i,j)), 0.)
1983                rsub1(i,j)=cs580*qsi(i,j)*rtair(i,j)*rtair(i,j)
1984                y3(i,j)=1./tair(i,j)
1985           dd(i,j)=y3(i,j)*(rn30a*y3(i,j)-rn30b)+rn30c*tair(i,j)/esi(i,j)
1986                y1(i,j)=206.18*ssi(i,j)/dd(i,j)
1987                pidep(i,j)=y1(i,j)*sqrt(r_nci*qi(i,j)/r00)
1988                dep(i,j)=dm(i,j)/(1.+rsub1(i,j))/d2t
1989                if(dm(i,j).gt.cmin2) then
1990                   a2=1.
1991                 if(pidep(i,j).gt.dep(i,j).and.pidep(i,j).gt.cmin2) then
1992                      a2=dep(i,j)/pidep(i,j)
1993                      pidep(i,j)=dep(i,j)
1994                 endif
1995                   psfi(i,j)=r2is*a2*.5*qi(i,j)*y1(i,j)/(sqrt(ami100) &
1996                           -sqrt(ami40))
1997                   elseif(dm(i,j).lt.-cmin2) then
1998 !
1999 !        SUBLIMATION TERMS USED ONLY WHEN SATURATION ADJUSTMENT FOR ICE
2000 !        IS TURNED OFF
2001 !
2002                   pidep(i,j)=0.
2003                   psfi(i,j)=0.
2004                else
2005                   pidep(i,j)=0.
2006                   psfi(i,j)=0.
2007                endif
2008             endif
2009
2010 !TTT***** QG=QG+MIN(PGDRY,PGWET)
2011 !*  9 * PGACS : ACCRETION OF QS BY QG (DGACS,WGACS: DRY AND WET)  ***9**
2012 !* 14 * DGACW : ACCRETION OF QC BY QG (QGACW FOR PGMLT)           **14**
2013 !* 16 * DGACR : ACCRETION OF QR TO QG (QGACR FOR PGMLT)           **16**
2014
2015             if(qc(i,j)+qr(i,j).lt.1.e-4) then
2016                ee1=.01
2017               else
2018                  ee1=1.
2019               endif
2020             ee2=0.09
2021             egs(i,j)=ee1*exp(ee2*tairc(i,j))
2022 !            EGS(I,J)=0.1 ! 6/15/02 tao's
2023             if (tair(i,j).ge.t0) egs(i,j)=1.0
2024             y1(i,j)=abs(vg(i,j)-vs(i,j))
2025             y2(i,j)=zs(i,j)*zg(i,j)
2026             y3(i,j)=5./y2(i,j)
2027             y4(i,j)=.08*y3(i,j)*y3(i,j)
2028             y5(i,j)=.05*y3(i,j)*y4(i,j)
2029             dd(i,j)=y1(i,j)*(y3(i,j)/zs(i,j)**5+y4(i,j)/zs(i,j)**3 &
2030                     +y5(i,j)/zs(i,j))
2031             pgacs(i,j)=r2ig*r2is*r9r*egs(i,j)*dd(i,j)
2032 !JJS 1/3/06 from Steve and Chunglin
2033             if (ihail.eq.1) then
2034                dgacs(i,j)=pgacs(i,j)
2035             else
2036                dgacs(i,j)=0.
2037             endif
2038 !JJS 1/3/06 from Steve and Chunglin
2039             wgacs(i,j)=r2ig*r2is*r9r*dd(i,j)
2040 !            WGACS(I,J)=0.  ! 6/15/02 tao's
2041             y1(i,j)=1./zg(i,j)**bg3
2042
2043             if(ihail .eq. 1) then
2044                dgacw(i,j)=r2ig*max(r14r*qc(i,j)*y1(i,j), 0.0)
2045             else
2046                dgacw(i,j)=r2ig*max(r14f*qc(i,j)*y1(i,j), 0.0)
2047             endif
2048
2049             qgacw(i,j)=dgacw(i,j)
2050             y1(i,j)=abs(vg(i,j)-vr(i,j))
2051             y2(i,j)=zr(i,j)*zg(i,j)
2052             y3(i,j)=5./y2(i,j)
2053             y4(i,j)=.08*y3(i,j)*y3(i,j)
2054             y5(i,j)=.05*y3(i,j)*y4(i,j)
2055             dd(i,j)=r16r*y1(i,j)*(y3(i,j)/zr(i,j)**5+y4(i,j)/zr(i,j)**3 &
2056                     +y5(i,j)/zr(i,j))
2057             dgacr(i,j)=r2ig*max(dd(i,j), 0.0)
2058             qgacr(i,j)=dgacr(i,j)
2059
2060             if (tair(i,j).ge.t0) then
2061                dgacs(i,j)=0.0
2062                wgacs(i,j)=0.0
2063                dgacw(i,j)=0.0
2064                dgacr(i,j)=0.0
2065             else
2066                pgacs(i,j)=0.0
2067                qgacw(i,j)=0.0
2068                qgacr(i,j)=0.0
2069             endif
2070
2071 !*******PGDRY : DGACW+DGACI+DGACR+DGACS                           ******
2072 !* 15 * DGACI : ACCRETION OF QI BY QG (WGACI FOR WET GROWTH)      **15**
2073 !* 17 * PGWET : WET GROWTH OF QG                                  **17**
2074
2075             dgaci(i,j)=0.0
2076             wgaci(i,j)=0.0
2077             pgwet(i,j)=0.0
2078
2079             if (tair(i,j).lt.t0) then
2080                y1(i,j)=qi(i,j)/zg(i,j)**bg3
2081                if (ihail.eq.1) then
2082                   dgaci(i,j)=r2ig*r15r*y1(i,j)
2083                   wgaci(i,j)=r2ig*r15ar*y1(i,j)
2084 !                  WGACI(I,J)=0.  ! 6/15/02 tao's
2085                else
2086
2087 !JJS                  DGACI(I,J)=r2ig*R15F*Y1(I,J)
2088                    dgaci(i,j)=0.
2089                   wgaci(i,j)=r2ig*r15af*y1(i,j)
2090 !                  WGACI(I,J)=0.  ! 6/15/02 tao's
2091                endif
2092 !
2093                if (tairc(i,j).ge.-50.) then
2094                 if (alf+rn17c*tairc(i,j) .eq. 0.) then
2095                    write(91,*) itimestep, i,j,k, alf, rn17c, tairc(i,j)
2096                 endif
2097                 y1(i,j)=1./(alf+rn17c*tairc(i,j))
2098                 if (ihail.eq.1) then
2099                    y3(i,j)=.78/zg(i,j)**2+r17aq*scv(i,j)/zg(i,j)**bgh5
2100                 else
2101                    y3(i,j)=.78/zg(i,j)**2+r17as*scv(i,j)/zg(i,j)**bgh5
2102                 endif
2103                 y4(i,j)=alvr*dwv(i,j)*(rp0-(qv(i,j)+qb0)) &
2104                         -tca(i,j)*tairc(i,j)
2105                 dd(i,j)=y1(i,j)*(r17r*y4(i,j)*y3(i,j) &
2106                        +(wgaci(i,j)+wgacs(i,j))*(alf+rn17b*tairc(i,j)))
2107                 pgwet(i,j)=r2ig*max(dd(i,j), 0.0)
2108                endif
2109             endif
2110 !JJS  125    CONTINUE
2111
2112 !********   HANDLING THE NEGATIVE CLOUD WATER (QC)    ******************
2113 !********   HANDLING THE NEGATIVE CLOUD ICE (QI)      ******************
2114
2115 !JJS         DO 150 J=3,JLES
2116 !JJS         DO 150 I=3,ILES
2117
2118             y1(i,j)=qc(i,j)/d2t
2119             psacw(i,j)=min(y1(i,j), psacw(i,j))
2120             praut(i,j)=min(y1(i,j), praut(i,j))
2121             pracw(i,j)=min(y1(i,j), pracw(i,j))
2122             psfw(i,j)= min(y1(i,j), psfw(i,j))
2123             dgacw(i,j)=min(y1(i,j), dgacw(i,j))
2124             qsacw(i,j)=min(y1(i,j), qsacw(i,j))
2125             qgacw(i,j)=min(y1(i,j), qgacw(i,j))
2126
2127             y1(i,j)=(psacw(i,j)+praut(i,j)+pracw(i,j)+psfw(i,j) &
2128                     +dgacw(i,j)+qsacw(i,j)+qgacw(i,j))*d2t
2129             qc(i,j)=qc(i,j)-y1(i,j)
2130
2131             if (qc(i,j) .lt. 0.0) then
2132                a1=1.
2133                if (y1(i,j) .ne. 0.0) a1=qc(i,j)/y1(i,j)+1.
2134                psacw(i,j)=psacw(i,j)*a1
2135                praut(i,j)=praut(i,j)*a1
2136                pracw(i,j)=pracw(i,j)*a1
2137                psfw(i,j)=psfw(i,j)*a1
2138                dgacw(i,j)=dgacw(i,j)*a1
2139                qsacw(i,j)=qsacw(i,j)*a1
2140                qgacw(i,j)=qgacw(i,j)*a1
2141                qc(i,j)=0.0
2142             endif
2143 !c
2144 !
2145 !******** SHED PROCESS (WGACR=PGWET-DGACW-WGACI-WGACS)
2146 !c
2147             wgacr(i,j)=pgwet(i,j)-dgacw(i,j)-wgaci(i,j)-wgacs(i,j)
2148             y2(i,j)=dgacw(i,j)+dgaci(i,j)+dgacr(i,j)+dgacs(i,j)
2149             if (pgwet(i,j).ge.y2(i,j)) then
2150                wgacr(i,j)=0.0
2151                wgaci(i,j)=0.0
2152                wgacs(i,j)=0.0
2153             else
2154                dgacr(i,j)=0.0
2155                dgaci(i,j)=0.0
2156                dgacs(i,j)=0.0
2157             endif
2158 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2159 !c
2160             y1(i,j)=qi(i,j)/d2t
2161             psaut(i,j)=min(y1(i,j), psaut(i,j))
2162             psaci(i,j)=min(y1(i,j), psaci(i,j))
2163             praci(i,j)=min(y1(i,j), praci(i,j))
2164             psfi(i,j)= min(y1(i,j), psfi(i,j))
2165             dgaci(i,j)=min(y1(i,j), dgaci(i,j))
2166             wgaci(i,j)=min(y1(i,j), wgaci(i,j))
2167 !
2168             y2(i,j)=(psaut(i,j)+psaci(i,j)+praci(i,j)+psfi(i,j) &
2169                    +dgaci(i,j)+wgaci(i,j))*d2t
2170             qi(i,j)=qi(i,j)-y2(i,j)+pidep(i,j)*d2t
2171
2172             if (qi(i,j).lt.0.0) then
2173                a2=1.
2174                if (y2(i,j) .ne. 0.0) a2=qi(i,j)/y2(i,j)+1.
2175                psaut(i,j)=psaut(i,j)*a2
2176                psaci(i,j)=psaci(i,j)*a2
2177                praci(i,j)=praci(i,j)*a2
2178                psfi(i,j)=psfi(i,j)*a2
2179                dgaci(i,j)=dgaci(i,j)*a2
2180                wgaci(i,j)=wgaci(i,j)*a2
2181                qi(i,j)=0.0
2182             endif
2183 !
2184             dlt3(i,j)=0.0
2185             dlt2(i,j)=0.0
2186 !
2187
2188 !            DLT4(I,J)=1.0
2189 !            if(qc(i,j) .gt. 5.e-4) dlt4(i,j)=0.0
2190 !            if(qs(i,j) .le. 1.e-4) dlt4(i,j)=1.0
2191 !
2192 !            IF (TAIR(I,J).ge.T0) THEN
2193 !               DLT4(I,J)=0.0
2194 !            ENDIF
2195
2196             if (tair(i,j).lt.t0) then
2197                if (qr(i,j).lt.1.e-4) then
2198                   dlt3(i,j)=1.0
2199                   dlt2(i,j)=1.0
2200                endif
2201                if (qs(i,j).ge.1.e-4) then
2202                   dlt2(i,j)=0.0
2203                endif
2204             endif
2205
2206             if (ice2 .eq. 1) then
2207                   dlt3(i,j)=1.0
2208                   dlt2(i,j)=1.0
2209             endif
2210 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2211             pr(i,j)=(qsacw(i,j)+praut(i,j)+pracw(i,j)+qgacw(i,j))*d2t
2212             ps(i,j)=(psaut(i,j)+psaci(i,j)+psacw(i,j)+psfw(i,j) &
2213                     +psfi(i,j)+dlt3(i,j)*praci(i,j))*d2t
2214 !           PS(I,J)=(PSAUT(I,J)+PSACI(I,J)+dlt4(i,j)*PSACW(I,J)
2215 !    1              +PSFW(I,J)+PSFI(I,J)+DLT3(I,J)*PRACI(I,J))*D2T
2216             pg(i,j)=((1.-dlt3(i,j))*praci(i,j)+dgaci(i,j)+wgaci(i,j) &
2217                     +dgacw(i,j))*d2t
2218 !           PG(I,J)=((1.-DLT3(I,J))*PRACI(I,J)+DGACI(I,J)+WGACI(I,J)
2219 !    1              +DGACW(I,J)+(1.-dlt4(i,j))*PSACW(I,J))*D2T
2220
2221 !JJS  150    CONTINUE
2222
2223 !*  7 * PRACS : ACCRETION OF QS BY QR                             ***7**
2224 !*  8 * PSACR : ACCRETION OF QR BY QS (QSACR FOR PSMLT)           ***8**
2225
2226 !JJS         DO 175 J=3,JLES
2227 !JJS         DO 175 I=3,ILES
2228
2229             y1(i,j)=abs(vr(i,j)-vs(i,j))
2230             y2(i,j)=zr(i,j)*zs(i,j)
2231             y3(i,j)=5./y2(i,j)
2232             y4(i,j)=.08*y3(i,j)*y3(i,j)
2233             y5(i,j)=.05*y3(i,j)*y4(i,j)
2234             pracs(i,j)=r2ig*r2is*r7r*y1(i,j)*(y3(i,j)/zs(i,j)**5 &
2235                       +y4(i,j)/zs(i,j)**3+y5(i,j)/zs(i,j))
2236             psacr(i,j)=r2is*r8r*y1(i,j)*(y3(i,j)/zr(i,j)**5 &
2237                       +y4(i,j)/zr(i,j)**3+y5(i,j)/zr(i,j))
2238             qsacr(i,j)=psacr(i,j)
2239
2240             if (tair(i,j).ge.t0) then
2241                pracs(i,j)=0.0
2242                psacr(i,j)=0.0
2243             else
2244                qsacr(i,j)=0.0
2245             endif
2246 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2247 !*  2 * PGAUT : AUTOCONVERSION OF QS TO QG                        ***2**
2248 !* 18 * PGFR : FREEZING OF QR TO QG                               **18**
2249
2250             pgaut(i,j)=0.0
2251             pgfr(i,j)=0.0
2252
2253             if (tair(i,j) .lt. t0) then
2254 !               Y1(I,J)=EXP(.09*TAIRC(I,J))
2255 !               PGAUT(I,J)=r2is*max(RN2*Y1(I,J)*(QS(I,J)-BND2), 0.0)
2256 !         IF(IHAIL.EQ.1) PGAUT(I,J)=max(RN2*Y1(I,J)*(QS(I,J)-BND2),0.0)
2257                y2(i,j)=exp(rn18a*(t0-tair(i,j)))
2258 !JJS              PGFR(I,J)=r2ig*max(R18R*(Y2(I,J)-1.)/ZR(I,J)**7., 0.0)
2259 !              pgfr(i,j)=r2ice*max(r18r*(y2(i,j)-1.)* &
2260 !                                    (zr(i,j)**(-7.)), 0.0)
2261 !        modify to prevent underflow on some computers (JD)
2262                temp = 1./zr(i,j)
2263                temp = temp*temp*temp*temp*temp*temp*temp
2264                pgfr(i,j)=r2ig*max(r18r*(y2(i,j)-1.)* &
2265                                     temp, 0.0)
2266             endif
2267
2268 !JJS  175    CONTINUE
2269
2270 !********   HANDLING THE NEGATIVE RAIN WATER (QR)    *******************
2271 !********   HANDLING THE NEGATIVE SNOW (QS)          *******************
2272
2273 !JJS         DO 200 J=3,JLES
2274 !JJS         DO 200 I=3,ILES
2275
2276             y1(i,j)=qr(i,j)/d2t
2277             y2(i,j)=-qg(i,j)/d2t
2278             piacr(i,j)=min(y1(i,j), piacr(i,j))
2279             dgacr(i,j)=min(y1(i,j), dgacr(i,j))
2280             wgacr(i,j)=min(y1(i,j), wgacr(i,j))
2281             wgacr(i,j)=max(y2(i,j), wgacr(i,j))
2282             psacr(i,j)=min(y1(i,j), psacr(i,j))
2283             pgfr(i,j)= min(y1(i,j), pgfr(i,j))
2284             del=0.
2285             if(wgacr(i,j) .lt. 0.) del=1.
2286             y1(i,j)=(piacr(i,j)+dgacr(i,j)+(1.-del)*wgacr(i,j) &
2287                     +psacr(i,j)+pgfr(i,j))*d2t
2288             qr(i,j)=qr(i,j)+pr(i,j)-y1(i,j)-del*wgacr(i,j)*d2t
2289             if (qr(i,j) .lt. 0.0) then
2290                a1=1.
2291                if(y1(i,j) .ne. 0.) a1=qr(i,j)/y1(i,j)+1.
2292                piacr(i,j)=piacr(i,j)*a1
2293                dgacr(i,j)=dgacr(i,j)*a1
2294                if (wgacr(i,j).gt.0.) wgacr(i,j)=wgacr(i,j)*a1
2295                pgfr(i,j)=pgfr(i,j)*a1
2296                psacr(i,j)=psacr(i,j)*a1
2297                qr(i,j)=0.0
2298             endif
2299 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2300             prn(i,j)=d2t*((1.-dlt3(i,j))*piacr(i,j)+dgacr(i,j) &
2301                      +wgacr(i,j)+(1.-dlt2(i,j))*psacr(i,j)+pgfr(i,j))
2302             ps(i,j)=ps(i,j)+d2t*(dlt3(i,j)*piacr(i,j) &
2303                     +dlt2(i,j)*psacr(i,j))
2304             pracs(i,j)=(1.-dlt2(i,j))*pracs(i,j)
2305             y1(i,j)=qs(i,j)/d2t
2306             pgacs(i,j)=min(y1(i,j), pgacs(i,j))
2307             dgacs(i,j)=min(y1(i,j), dgacs(i,j))
2308             wgacs(i,j)=min(y1(i,j), wgacs(i,j))
2309             pgaut(i,j)=min(y1(i,j), pgaut(i,j))
2310             pracs(i,j)=min(y1(i,j), pracs(i,j))
2311             psn(i,j)=d2t*(pgacs(i,j)+dgacs(i,j)+wgacs(i,j) &
2312                      +pgaut(i,j)+pracs(i,j))
2313             qs(i,j)=qs(i,j)+ps(i,j)-psn(i,j)
2314
2315             if (qs(i,j).lt.0.0) then
2316                a2=1.
2317                if (psn(i,j) .ne. 0.0) a2=qs(i,j)/psn(i,j)+1.
2318                pgacs(i,j)=pgacs(i,j)*a2
2319                dgacs(i,j)=dgacs(i,j)*a2
2320                wgacs(i,j)=wgacs(i,j)*a2
2321                pgaut(i,j)=pgaut(i,j)*a2
2322                pracs(i,j)=pracs(i,j)*a2
2323                psn(i,j)=psn(i,j)*a2
2324                qs(i,j)=0.0
2325             endif
2326 !
2327 !C           PSN(I,J)=D2T*(PGACS(I,J)+DGACS(I,J)+WGACS(I,J)
2328 !c                    +PGAUT(I,J)+PRACS(I,J))
2329             y2(i,j)=d2t*(psacw(i,j)+psfw(i,j)+dgacw(i,j)+piacr(i,j) &
2330                     +dgacr(i,j)+wgacr(i,j)+psacr(i,j)+pgfr(i,j))
2331             pt(i,j)=pt(i,j)+afcp*y2(i,j)
2332             qg(i,j)=qg(i,j)+pg(i,j)+prn(i,j)+psn(i,j)
2333
2334 !JJS  200    CONTINUE
2335
2336 !* 11 * PSMLT : MELTING OF QS                                     **11**
2337 !* 19 * PGMLT : MELTING OF QG TO QR                               **19**
2338
2339 !JJS         DO 225 J=3,JLES
2340 !JJS         DO 225 I=3,ILES
2341
2342             psmlt(i,j)=0.0
2343             pgmlt(i,j)=0.0
2344             tair(i,j)=(pt(i,j)+tb0)*pi0
2345
2346             if (tair(i,j).ge.t0) then
2347                tairc(i,j)=tair(i,j)-t0
2348                y1(i,j)=tca(i,j)*tairc(i,j)-alvr*dwv(i,j) &
2349                                *(rp0-(qv(i,j)+qb0))
2350                y2(i,j)=.78/zs(i,j)**2+r101f*scv(i,j)/zs(i,j)**bsh5
2351                dd(i,j)=r11rt*y1(i,j)*y2(i,j)+r11at*tairc(i,j) &
2352                        *(qsacw(i,j)+qsacr(i,j))
2353                psmlt(i,j)=r2is*max(0.0, min(dd(i,j), qs(i,j)))
2354
2355                if(ihail.eq.1) then
2356                   y3(i,j)=.78/zg(i,j)**2+r19aq*scv(i,j)/zg(i,j)**bgh5
2357                else
2358                   y3(i,j)=.78/zg(i,j)**2+r19as*scv(i,j)/zg(i,j)**bgh5
2359                endif
2360
2361                dd1(i,j)=r19rt*y1(i,j)*y3(i,j)+r19bt*tairc(i,j) &
2362                         *(qgacw(i,j)+qgacr(i,j))
2363                pgmlt(i,j)=r2ig*max(0.0, min(dd1(i,j), qg(i,j)))
2364                pt(i,j)=pt(i,j)-afcp*(psmlt(i,j)+pgmlt(i,j))
2365                qr(i,j)=qr(i,j)+psmlt(i,j)+pgmlt(i,j)
2366                qs(i,j)=qs(i,j)-psmlt(i,j)
2367                qg(i,j)=qg(i,j)-pgmlt(i,j)
2368             endif
2369 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2370 !* 24 * PIHOM : HOMOGENEOUS FREEZING OF QC TO QI (T < T00)        **24**
2371 !* 25 * PIDW : DEPOSITION GROWTH OF QC TO QI ( T0 < T <= T00)     **25**
2372 !* 26 * PIMLT : MELTING OF QI TO QC (T >= T0)                     **26**
2373
2374             if (qc(i,j).le.cmin1) qc(i,j)=0.0
2375             if (qi(i,j).le.cmin1) qi(i,j)=0.0
2376             tair(i,j)=(pt(i,j)+tb0)*pi0
2377
2378             if(tair(i,j).le.t00) then
2379                pihom(i,j)=qc(i,j)
2380             else
2381                pihom(i,j)=0.0
2382             endif
2383             if(tair(i,j).ge.t0) then
2384                pimlt(i,j)=qi(i,j)
2385             else
2386                pimlt(i,j)=0.0
2387             endif
2388             pidw(i,j)=0.0
2389
2390             if (tair(i,j).lt.t0 .and. tair(i,j).gt.t00) then
2391                tairc(i,j)=tair(i,j)-t0
2392                y1(i,j)=max( min(tairc(i,j), -1.), -31.)
2393                it(i,j)=int(abs(y1(i,j)))
2394                y2(i,j)=aa1(it(i,j))
2395                y3(i,j)=aa2(it(i,j))
2396                y4(i,j)=exp(abs(beta*tairc(i,j)))
2397                y5(i,j)=(r00*qi(i,j)/(r25a*y4(i,j)))**y3(i,j)
2398                pidw(i,j)=min(r25rt*y2(i,j)*y4(i,j)*y5(i,j), qc(i,j))
2399             endif
2400
2401             y1(i,j)=pihom(i,j)-pimlt(i,j)+pidw(i,j)
2402             pt(i,j)=pt(i,j)+afcp*y1(i,j)+ascp*(pidep(i,j))*d2t
2403             qv(i,j)=qv(i,j)-(pidep(i,j))*d2t
2404             qc(i,j)=qc(i,j)-y1(i,j)
2405             qi(i,j)=qi(i,j)+y1(i,j)
2406
2407 !* 31 * PINT  : INITIATION OF QI                                  **31**
2408 !* 32 * PIDEP : DEPOSITION OF QI                                  **32**
2409 !
2410 !     CALCULATION OF PINT USES DIFFERENT VALUES OF THE INTERCEPT AND SLOPE FOR
2411 !     THE FLETCHER EQUATION. ALSO, ONLY INITIATE MORE ICE IF THE NEW NUMBER
2412 !     CONCENTRATION EXCEEDS THAT ALREADY PRESENT.
2413 !* 31 * pint  : initiation of qi                                  **31**
2414 !* 32 * pidep : deposition of qi                                  **32**
2415            pint(i,j)=0.0
2416 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2417         if ( itaobraun.eq.1 ) then
2418             tair(i,j)=(pt(i,j)+tb0)*pi0
2419             if (tair(i,j) .lt. t0) then
2420 !             if (qi(i,j) .le. cmin) qi(i,j)=0.
2421               if (qi(i,j) .le. cmin2) qi(i,j)=0.
2422                tairc(i,j)=tair(i,j)-t0
2423                rtair(i,j)=1./(tair(i,j)-c76)
2424                y2(i,j)=exp(c218-c580*rtair(i,j))
2425               qsi(i,j)=rp0*y2(i,j)
2426                esi(i,j)=c610*y2(i,j)
2427               ssi(i,j)=(qv(i,j)+qb0)/qsi(i,j)-1.
2428                         ami50=3.76e-8
2429
2430 !ccif ( itaobraun.eq.1 ) --> betah=0.5*beta=-.46*0.5=-0.23;   cn0=1.e-6
2431 !ccif ( itaobraun.eq.0 ) --> betah=0.5*beta=-.6*0.5=-0.30;    cn0=1.e-8
2432
2433              y1(i,j)=1./tair(i,j)
2434
2435 !cc insert a restriction on ice collection that ice collection
2436 !cc should be stopped at -30 c (with cn0=1.e-6, beta=-.46)
2437
2438              tairccri=tairc(i,j)          ! in degree c
2439              if(tairccri.le.-30.) tairccri=-30.
2440
2441 !            y2(i,j)=exp(betah*tairc(i,j))
2442              y2(i,j)=exp(betah*tairccri)
2443              y3(i,j)=sqrt(qi(i,j))
2444              dd(i,j)=y1(i,j)*(rn10a*y1(i,j)-rn10b)+rn10c*tair(i,j) &
2445                                                 /esi(i,j)
2446           pidep(i,j)=max(r32rt*ssi(i,j)*y2(i,j)*y3(i,j)/dd(i,j), 0.e0)
2447
2448            r_nci=min(cn0*exp(beta*tairc(i,j)),1.)
2449 !          r_nci=min(1.e-6*exp(-.46*tairc(i,j)),1.)
2450
2451            dd(i,j)=max(1.e-9*r_nci/r00-qi(i,j)*1.e-9/ami50, 0.)
2452                 dm(i,j)=max( (qv(i,j)+qb0-qsi(i,j)), 0.0)
2453                 rsub1(i,j)=cs580*qsi(i,j)*rtair(i,j)*rtair(i,j)
2454               dep(i,j)=dm(i,j)/(1.+rsub1(i,j))
2455               pint(i,j)=max(min(dd(i,j), dm(i,j)), 0.)
2456
2457 !             pint(i,j)=min(pint(i,j), dep(i,j))
2458               pint(i,j)=min(pint(i,j)+pidep(i,j), dep(i,j))
2459
2460 !              if (pint(i,j) .le. cmin) pint(i,j)=0.
2461                if (pint(i,j) .le. cmin2) pint(i,j)=0.
2462               pt(i,j)=pt(i,j)+ascp*pint(i,j)
2463               qv(i,j)=qv(i,j)-pint(i,j)
2464               qi(i,j)=qi(i,j)+pint(i,j)
2465            endif
2466         endif  ! if ( itaobraun.eq.1 )
2467 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2468 !
2469 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2470         if ( itaobraun.eq.0 ) then
2471              tair(i,j)=(pt(i,j)+tb0)*pi0
2472              if (tair(i,j) .lt. t0) then
2473                if (qi(i,j) .le. cmin1) qi(i,j)=0.
2474                tairc(i,j)=tair(i,j)-t0
2475                dd(i,j)=r31r*exp(beta*tairc(i,j))
2476                rtair(i,j)=1./(tair(i,j)-c76)
2477                y2(i,j)=exp(c218-c580*rtair(i,j))
2478                qsi(i,j)=rp0*y2(i,j)
2479                esi(i,j)=c610*y2(i,j)
2480                ssi(i,j)=(qv(i,j)+qb0)/qsi(i,j)-1.
2481                dm(i,j)=max( (qv(i,j)+qb0-qsi(i,j)), 0.)
2482                rsub1(i,j)=cs580*qsi(i,j)*rtair(i,j)*rtair(i,j)
2483                dep(i,j)=dm(i,j)/(1.+rsub1(i,j))
2484               pint(i,j)=max(min(dd(i,j), dm(i,j)), 0.)
2485                y1(i,j)=1./tair(i,j)
2486                y2(i,j)=exp(betah*tairc(i,j))
2487                y3(i,j)=sqrt(qi(i,j))
2488                dd(i,j)=y1(i,j)*(rn10a*y1(i,j)-rn10b) &
2489                      +rn10c*tair(i,j)/esi(i,j)
2490              pidep(i,j)=max(r32rt*ssi(i,j)*y2(i,j)*y3(i,j)/dd(i,j), 0.)
2491               pint(i,j)=pint(i,j)+pidep(i,j)
2492               pint(i,j)=min(pint(i,j),dep(i,j))
2493 !c          if (pint(i,j) .le. cmin2) pint(i,j)=0.
2494              pt(i,j)=pt(i,j)+ascp*pint(i,j)
2495              qv(i,j)=qv(i,j)-pint(i,j)
2496              qi(i,j)=qi(i,j)+pint(i,j)
2497             endif
2498         endif  ! if ( itaobraun.eq.0 )
2499 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2500
2501 !JJS  225    CONTINUE
2502
2503 !*****   TAO ET AL (1989) SATURATION TECHNIQUE  ***********************
2504
2505          if (new_ice_sat .eq. 0) then
2506
2507 !JJS            DO 250 J=3,JLES
2508 !JJS            DO 250 I=3,ILES
2509                tair(i,j)=(pt(i,j)+tb0)*pi0
2510                cnd(i,j)=rt0*(tair(i,j)-t00)
2511                dep(i,j)=rt0*(t0-tair(i,j))
2512                y1(i,j)=1./(tair(i,j)-c358)
2513                y2(i,j)=1./(tair(i,j)-c76)
2514                qsw(i,j)=rp0*exp(c172-c409*y1(i,j))
2515                qsi(i,j)=rp0*exp(c218-c580*y2(i,j))
2516                dd(i,j)=cp409*y1(i,j)*y1(i,j)
2517                dd1(i,j)=cp580*y2(i,j)*y2(i,j)
2518                if (qc(i,j).le.cmin) qc(i,j)=cmin
2519                if (qi(i,j).le.cmin) qi(i,j)=cmin
2520                if (tair(i,j).ge.t0) then
2521                   dep(i,j)=0.0
2522                   cnd(i,j)=1.
2523                   qi(i,j)=0.0
2524                endif
2525
2526                if (tair(i,j).lt.t00) then
2527                   cnd(i,j)=0.0
2528                   dep(i,j)=1.
2529                   qc(i,j)=0.0
2530                endif
2531
2532                y5(i,j)=avcp*cnd(i,j)+ascp*dep(i,j)
2533 !               if (qc(i,j) .ge. cmin .or. qi(i,j) .ge. cmin) then
2534                y1(i,j)=qc(i,j)*qsw(i,j)/(qc(i,j)+qi(i,j))
2535                y2(i,j)=qi(i,j)*qsi(i,j)/(qc(i,j)+qi(i,j))
2536                y4(i,j)=dd(i,j)*y1(i,j)+dd1(i,j)*y2(i,j)
2537                qvs(i,j)=y1(i,j)+y2(i,j)
2538                rsub1(i,j)=(qv(i,j)+qb0-qvs(i,j))/(1.+y4(i,j)*y5(i,j))
2539                cnd(i,j)=cnd(i,j)*rsub1(i,j)
2540                dep(i,j)=dep(i,j)*rsub1(i,j)
2541                if (qc(i,j).le.cmin) qc(i,j)=0.
2542                if (qi(i,j).le.cmin) qi(i,j)=0.
2543 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2544 !c    ******   condensation or evaporation of qc  ******
2545
2546                cnd(i,j)=max(-qc(i,j),cnd(i,j))
2547
2548 !c    ******   deposition or sublimation of qi    ******
2549
2550                dep(i,j)=max(-qi(i,j),dep(i,j))
2551
2552                pt(i,j)=pt(i,j)+avcp*cnd(i,j)+ascp*dep(i,j)
2553                qv(i,j)=qv(i,j)-cnd(i,j)-dep(i,j)
2554                qc(i,j)=qc(i,j)+cnd(i,j)
2555                qi(i,j)=qi(i,j)+dep(i,j)
2556 !JJS  250       continue
2557          endif
2558
2559          if (new_ice_sat .eq. 1) then
2560
2561 !JJS            DO J=3,JLES
2562 !JJS            DO I=3,ILES
2563
2564                tair(i,j)=(pt(i,j)+tb0)*pi0
2565                cnd(i,j)=rt0*(tair(i,j)-t00)
2566                dep(i,j)=rt0*(t0-tair(i,j))
2567                y1(i,j)=1./(tair(i,j)-c358)
2568                y2(i,j)=1./(tair(i,j)-c76)
2569                qsw(i,j)=rp0*exp(c172-c409*y1(i,j))
2570                qsi(i,j)=rp0*exp(c218-c580*y2(i,j))
2571                dd(i,j)=cp409*y1(i,j)*y1(i,j)
2572                dd1(i,j)=cp580*y2(i,j)*y2(i,j)
2573                y5(i,j)=avcp*cnd(i,j)+ascp*dep(i,j)
2574                y1(i,j)=rt0*(tair(i,j)-t00)*qsw(i,j)
2575                y2(i,j)=rt0*(t0-tair(i,j))*qsi(i,j)
2576 !               IF (QC(I,J).LE.CMIN) QC(I,J)=CMIN
2577 !               IF (QI(I,J).LE.CMIN) QI(I,J)=CMIN
2578
2579                if (tair(i,j).ge.t0) then
2580 !                 QI(I,J)=0.0
2581                   dep(i,j)=0.0
2582                   cnd(i,j)=1.
2583                   y2(i,j)=0.
2584                   y1(i,j)=qsw(i,j)
2585                endif
2586                if (tair(i,j).lt.t00) then
2587                   cnd(i,j)=0.0
2588                   dep(i,j)=1.
2589                   y2(i,j)=qsi(i,j)
2590                   y1(i,j)=0.
2591 !                 QC(I,J)=0.0
2592                endif
2593
2594 !            Y1(I,J)=QC(I,J)*QSW(I,J)/(QC(I,J)+QI(I,J))
2595 !            Y2(I,J)=QI(I,J)*QSI(I,J)/(QC(I,J)+QI(I,J))
2596
2597                y4(i,j)=dd(i,j)*y1(i,j)+dd1(i,j)*y2(i,j)
2598                qvs(i,j)=y1(i,j)+y2(i,j)
2599                rsub1(i,j)=(qv(i,j)+qb0-qvs(i,j))/(1.+y4(i,j)*y5(i,j))
2600                cnd(i,j)=cnd(i,j)*rsub1(i,j)
2601                dep(i,j)=dep(i,j)*rsub1(i,j)
2602 !               IF (QC(I,J).LE.CMIN) QC(I,J)=0.
2603 !               IF (QI(I,J).LE.CMIN) QI(I,J)=0.
2604
2605 !C    ******   CONDENSATION OR EVAPORATION OF QC  ******
2606
2607                cnd(i,j)=max(-qc(i,j),cnd(i,j))
2608
2609 !C    ******   DEPOSITION OR SUBLIMATION OF QI    ******
2610
2611                dep(i,j)=max(-qi(i,j),dep(i,j))
2612
2613                pt(i,j)=pt(i,j)+avcp*cnd(i,j)+ascp*dep(i,j)
2614                qv(i,j)=qv(i,j)-cnd(i,j)-dep(i,j)
2615                qc(i,j)=qc(i,j)+cnd(i,j)
2616                qi(i,j)=qi(i,j)+dep(i,j)
2617 !JJS            ENDDO
2618 !JJS            ENDDO
2619          endif
2620
2621 !c
2622 !
2623           if (new_ice_sat .eq. 2) then
2624 !JJS          do j=3,jles
2625 !JJS             do i=3,iles
2626           dep(i,j)=0.0
2627           cnd(i,j)=0.0
2628           tair(i,j)=(pt(i,j)+tb0)*pi0
2629           if (tair(i,j) .ge. 253.16) then
2630               y1(i,j)=1./(tair(i,j)-c358)
2631               qsw(i,j)=rp0*exp(c172-c409*y1(i,j))
2632               dd(i,j)=cp409*y1(i,j)*y1(i,j)
2633               dm(i,j)=qv(i,j)+qb0-qsw(i,j)
2634               cnd(i,j)=dm(i,j)/(1.+avcp*dd(i,j)*qsw(i,j))
2635 !c    ******   condensation or evaporation of qc  ******
2636               cnd(i,j)=max(-qc(i,j), cnd(i,j))
2637              pt(i,j)=pt(i,j)+avcp*cnd(i,j)
2638              qv(i,j)=qv(i,j)-cnd(i,j)
2639              qc(i,j)=qc(i,j)+cnd(i,j)
2640          endif
2641           if (tair(i,j) .le. 258.16) then
2642 !c             cnd(i,j)=0.0
2643            y2(i,j)=1./(tair(i,j)-c76)
2644            qsi(i,j)=rp0*exp(c218-c580*y2(i,j))
2645           dd1(i,j)=cp580*y2(i,j)*y2(i,j)
2646          dep(i,j)=(qv(i,j)+qb0-qsi(i,j))/(1.+ascp*dd1(i,j)*qsi(i,j))
2647 !c    ******   deposition or sublimation of qi    ******
2648              dep(i,j)=max(-qi(i,j),dep(i,j))
2649              pt(i,j)=pt(i,j)+ascp*dep(i,j)
2650              qv(i,j)=qv(i,j)-dep(i,j)
2651              qi(i,j)=qi(i,j)+dep(i,j)
2652          endif
2653 !JJS       enddo
2654 !JJS       enddo
2655       endif
2656
2657 !c
2658 !
2659 !* 10 * PSDEP : DEPOSITION OR SUBLIMATION OF QS                   **10**
2660 !* 20 * PGSUB : SUBLIMATION OF QG                                 **20**
2661
2662 !JJS         DO 275 J=3,JLES
2663 !JJS         DO 275 I=3,ILES
2664
2665             psdep(i,j)=0.0
2666             pgdep(i,j)=0.0
2667             pssub(i,j)=0.0
2668             pgsub(i,j)=0.0
2669             tair(i,j)=(pt(i,j)+tb0)*pi0
2670
2671             if(tair(i,j).lt.t0) then
2672                if(qs(i,j).lt.cmin1) qs(i,j)=0.0
2673                if(qg(i,j).lt.cmin1) qg(i,j)=0.0
2674                rtair(i,j)=1./(tair(i,j)-c76)
2675                qsi(i,j)=rp0*exp(c218-c580*rtair(i,j))
2676                ssi(i,j)=(qv(i,j)+qb0)/qsi(i,j)-1.
2677 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2678                y1(i,j)=r10ar/(tca(i,j)*tair(i,j)**2)+1./(dwv(i,j) &
2679                       *qsi(i,j))
2680                y2(i,j)=.78/zs(i,j)**2+r101f*scv(i,j)/zs(i,j)**bsh5
2681                psdep(i,j)=r10t*ssi(i,j)*y2(i,j)/y1(i,j)
2682                pssub(i,j)=psdep(i,j)
2683                psdep(i,j)=r2is*max(psdep(i,j), 0.)
2684                pssub(i,j)=r2is*max(-qs(i,j), min(pssub(i,j), 0.))
2685
2686                if(ihail.eq.1) then
2687                   y2(i,j)=.78/zg(i,j)**2+r20bq*scv(i,j)/zg(i,j)**bgh5
2688                else
2689                   y2(i,j)=.78/zg(i,j)**2+r20bs*scv(i,j)/zg(i,j)**bgh5
2690                endif
2691
2692                pgsub(i,j)=r2ig*r20t*ssi(i,j)*y2(i,j)/y1(i,j)
2693                dm(i,j)=qv(i,j)+qb0-qsi(i,j)
2694                rsub1(i,j)=cs580*qsi(i,j)*rtair(i,j)*rtair(i,j)
2695
2696 !     ********   DEPOSITION OR SUBLIMATION OF QS  **********************
2697
2698                y1(i,j)=dm(i,j)/(1.+rsub1(i,j))
2699                psdep(i,j)=r2is*min(psdep(i,j),max(y1(i,j),0.))
2700                y2(i,j)=min(y1(i,j),0.)
2701                pssub(i,j)=r2is*max(pssub(i,j),y2(i,j))
2702
2703 !     ********   SUBLIMATION OF QG   ***********************************
2704
2705                dd(i,j)=max((-y2(i,j)-qs(i,j)), 0.)
2706               pgsub(i,j)=r2ig*min(dd(i,j), qg(i,j), max(pgsub(i,j),0.))
2707
2708                if(qc(i,j)+qi(i,j).gt.1.e-5) then
2709                   dlt1(i,j)=1.
2710                else
2711                   dlt1(i,j)=0.
2712                endif
2713
2714                psdep(i,j)=dlt1(i,j)*psdep(i,j)
2715                pssub(i,j)=(1.-dlt1(i,j))*pssub(i,j)
2716                pgsub(i,j)=(1.-dlt1(i,j))*pgsub(i,j)
2717
2718                pt(i,j)=pt(i,j)+ascp*(psdep(i,j)+pssub(i,j)-pgsub(i,j))
2719                qv(i,j)=qv(i,j)+pgsub(i,j)-pssub(i,j)-psdep(i,j)
2720                qs(i,j)=qs(i,j)+psdep(i,j)+pssub(i,j)
2721                qg(i,j)=qg(i,j)-pgsub(i,j)
2722             endif
2723
2724 !* 23 * ERN : EVAPORATION OF QR (SUBSATURATION)                   **23**
2725
2726             ern(i,j)=0.0
2727
2728             if(qr(i,j).gt.0.0) then
2729                tair(i,j)=(pt(i,j)+tb0)*pi0
2730                rtair(i,j)=1./(tair(i,j)-c358)
2731                qsw(i,j)=rp0*exp(c172-c409*rtair(i,j))
2732                ssw(i,j)=(qv(i,j)+qb0)/qsw(i,j)-1.0
2733                dm(i,j)=qv(i,j)+qb0-qsw(i,j)
2734                rsub1(i,j)=cv409*qsw(i,j)*rtair(i,j)*rtair(i,j)
2735                dd1(i,j)=max(-dm(i,j)/(1.+rsub1(i,j)), 0.0)
2736                y1(i,j)=.78/zr(i,j)**2+r23af*scv(i,j)/zr(i,j)**bwh5
2737                y2(i,j)=r23br/(tca(i,j)*tair(i,j)**2)+1./(dwv(i,j) &
2738                        *qsw(i,j))
2739 !cccc
2740                ern(i,j)=r23t*ssw(i,j)*y1(i,j)/y2(i,j)
2741                ern(i,j)=min(dd1(i,j),qr(i,j),max(ern(i,j),0.))
2742                pt(i,j)=pt(i,j)-avcp*ern(i,j)
2743                qv(i,j)=qv(i,j)+ern(i,j)
2744                qr(i,j)=qr(i,j)-ern(i,j)
2745             endif
2746
2747 !JJS 10/7/2008     vvvvv
2748     ENDIF    ! part of if (iwarm.eq.1) then
2749 !JJS 10/7/2008     ^^^^^
2750
2751 !            IF (QV(I,J)+QB0 .LE. 0.) QV(I,J)=-QB0
2752             if (qc(i,j) .le. cmin1) qc(i,j)=0.
2753             if (qr(i,j) .le. cmin1) qr(i,j)=0.
2754             if (qi(i,j) .le. cmin1) qi(i,j)=0.
2755             if (qs(i,j) .le. cmin1) qs(i,j)=0.
2756             if (qg(i,j) .le. cmin1) qg(i,j)=0.
2757             dpt(i,j,k)=pt(i,j)
2758             dqv(i,j,k)=qv(i,j)
2759             qcl(i,j,k)=qc(i,j)
2760             qrn(i,j,k)=qr(i,j)
2761             qci(i,j,k)=qi(i,j)
2762             qcs(i,j,k)=qs(i,j)
2763             qcg(i,j,k)=qg(i,j)
2764
2765 !JJS  275    CONTINUE
2766
2767          scc=0.
2768          see=0.
2769
2770 !         DO 110 J=3,JLES
2771 !         DO 110 I=3,ILES
2772 !            DD(I,J)=MAX(-CND(I,J), 0.)
2773 !            CND(I,J)=MAX(CND(I,J), 0.)
2774 !            DD1(I,J)=MAX(-DEP(I,J), 0.)
2775
2776 !ccshie 2/21/02 shie follow tao
2777 !CC for reference    QI(I,J)=QI(I,J)-Y2(I,J)+PIDEP(I,J)*D2T
2778 !CC for reference    QV(I,J)=QV(I,J)-(PIDEP(I,J))*D2T
2779
2780 !c            DEP(I,J)=MAX(DEP(I,J), 0.)
2781 !            DEP(I,J)=MAX(DEP(I,J), 0.)+PIDEP(I,J)*D2T
2782 !            SCC=SCC+CND(I,J)
2783 !            SEE=SEE+DD(I,J)+ERN(I,J)
2784
2785 !  110    CONTINUE
2786
2787 !         SC(K)=SCC+SC(K)
2788 !         SE(K)=SEE+SE(K)
2789
2790 !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2791 !c     henry:  please take a look  (start)
2792 !JJS modified by JJS on 5/1/2007  vvvvv
2793
2794 !JJS       do 305 j=3,jles
2795 !JJS       do 305 i=3,iles
2796             dd(i,j)=max(-cnd(i,j), 0.)
2797             cnd(i,j)=max(cnd(i,j), 0.)
2798             dd1(i,j)=max(-dep(i,j), 0.)+pidep(i,j)*d2t
2799             dep(i,j)=max(dep(i,j), 0.)
2800 !JJS  305  continue
2801
2802 !JJS       do 306 j=3,jles
2803 !JJS       do 306 i=3,iles
2804 !JJS              scc=scc+cnd(i,j)
2805 !JJS              see=see+(dd(i,j)+ern(i,j))
2806 !
2807 !JJS            sddd=sddd+(dep(i,j)+pint(i,j)+psdep(i,j)+pgdep(i,j))
2808 !JJS          ssss=ssss+dd1(i,j)
2809 !JJS
2810 !            shhh=shhh+rsw(i,j,k)*d2t
2811 !            sccc=sccc+rlw(i,j,k)*d2t
2812 !jjs
2813 !JJS              smmm=smmm+(psmlt(i,j)+pgmlt(i,j)+pimlt(i,j))
2814 !JJS              sfff=sfff+d2t*(psacw(i,j)+piacr(i,j)+psfw(i,j)
2815 !JJS     1         +pgfr(i,j)+dgacw(i,j)+dgacr(i,j)+psacr(i,j))
2816 !JJS     2        -qracs(i,j)+pihom(i,j)+pidw(i,j)
2817
2818
2819               sccc=cnd(i,j)
2820               seee=dd(i,j)+ern(i,j)
2821               sddd=dep(i,j)+pint(i,j)+psdep(i,j)+pgdep(i,j)
2822               ssss=dd1(i,j) + pgsub(i,j)
2823               smmm=psmlt(i,j)+pgmlt(i,j)+pimlt(i,j)
2824               sfff=d2t*(psacw(i,j)+piacr(i,j)+psfw(i,j) &
2825                +pgfr(i,j)+dgacw(i,j)+dgacr(i,j)+psacr(i,j) &
2826                +wgacr(i,j))+pihom(i,j)+pidw(i,j)
2827
2828 !           physc(i,k,j) = avcp * sccc / d2t
2829 !           physe(i,k,j) = avcp * seee / d2t
2830 !           physd(i,k,j) = ascp * sddd / d2t
2831 !           physs(i,k,j) = ascp * ssss / d2t
2832 !           physf(i,k,j) = afcp * sfff / d2t
2833 !           physm(i,k,j) = afcp * smmm / d2t
2834 !           physc(i,k,j) = physc(i,k,j) + avcp * sccc
2835 !           physe(i,k,j) = physc(i,k,j) + avcp * seee
2836 !           physd(i,k,j) = physd(i,k,j) + ascp * sddd
2837 !           physs(i,k,j) = physs(i,k,j) + ascp * ssss
2838 !           physf(i,k,j) = physf(i,k,j) + afcp * sfff
2839 !           physm(i,k,j) = physm(i,k,j) + afcp * smmm
2840
2841 !JJS modified by JJS on 5/1/2007  ^^^^^
2842
2843  2000 continue
2844
2845  1000 continue
2846
2847 !JJS  ****************************************************************
2848 !JJS  convert from GCE grid back to WRF grid
2849       do k=kts,kte
2850          do j=jts,jte
2851          do i=its,ite
2852          ptwrf(i,k,j) = dpt(i,j,k)
2853          qvwrf(i,k,j) = dqv(i,j,k)
2854          qlwrf(i,k,j) = qcl(i,j,k)
2855          qrwrf(i,k,j) = qrn(i,j,k)
2856          qiwrf(i,k,j) = qci(i,j,k)
2857          qswrf(i,k,j) = qcs(i,j,k)
2858          qgwrf(i,k,j) = qcg(i,j,k)
2859          enddo !i
2860          enddo !j
2861       enddo !k
2862
2863 !     ****************************************************************
2864
2865 !+---+-----------------------------------------------------------------+
2866          IF ( PRESENT (diagflag) ) THEN
2867          if (diagflag .and. do_radar_ref == 1) then
2868             do j=jts,jte
2869             do i=its,ite
2870                DO K=kts,kte
2871                   t1d(k)=ptwrf(i,k,j)*pi_mks(i,k,j)
2872                   p1d(k)=p0_mks(i,k,j)
2873                   qv1d(k)=qvwrf(i,k,j)
2874                   qr1d(k)=qrwrf(i,k,j)
2875                ENDDO
2876                if (ice2.eq.0) then
2877                   DO K=kts,kte
2878                      qs1d(k)=qswrf(i,k,j)
2879                      qg1d(k)=qgwrf(i,k,j)
2880                   ENDDO
2881                elseif (ice2.eq.1) then
2882                   DO K=kts,kte
2883                      qs1d(k)=qswrf(i,k,j)
2884                   ENDDO
2885                elseif (ice2.eq.2) then
2886                   DO K=kts,kte
2887                      qs1d(k)=0.
2888                      qg1d(k)=qgwrf(i,k,j)
2889                   ENDDO
2890                elseif (ice2.eq.3) then
2891                   DO K=kts,kte
2892                      qs1d(k)=0.
2893                      qg1d(k)=0.
2894                   ENDDO
2895                endif
2896                call refl10cm_gsfc (qv1d, qr1d, qs1d, qg1d,             &
2897                        t1d, p1d, dBZ, kts, kte, i, j, ihail)
2898                do k = kts, kte
2899                   refl_10cm(i,k,j) = MAX(-35., dBZ(k))
2900                enddo
2901             enddo
2902             enddo
2903          endif
2904          ENDIF
2905 !+---+-----------------------------------------------------------------+
2906
2907       return
2908  END SUBROUTINE saticel_s
2909
2910 !JJS
2911 !JJS      REAL FUNCTION GAMMA(X)
2912 !JJS        Y=GAMMLN(X)
2913 !JJS        GAMMA=EXP(Y)
2914 !JJS      RETURN
2915 !JJS      END
2916 !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2917 !JJS      real function GAMMLN (xx)
2918   real function gammagce (xx)
2919 !**********************************************************************
2920   real*8 cof(6),stp,half,one,fpf,x,tmp,ser
2921   data cof,stp /  76.18009173,-86.50532033,24.01409822, &
2922      -1.231739516,.120858003e-2,-.536382e-5, 2.50662827465 /
2923   data half,one,fpf / .5, 1., 5.5 /
2924 !
2925       x=xx-one
2926       tmp=x+fpf
2927       tmp=(x+half)*log(tmp)-tmp
2928       ser=one
2929       do  j=1,6
2930          x=x+one
2931         ser=ser+cof(j)/x
2932       enddo !j
2933       gammln=tmp+log(stp*ser)
2934 !JJS
2935       gammagce=exp(gammln)
2936 !JJS
2937       return
2938  END FUNCTION gammagce
2939
2940 !+---+-----------------------------------------------------------------+
2941
2942       subroutine refl10cm_gsfc (qv1d, qr1d, qs1d, qg1d,                 &
2943                        t1d, p1d, dBZ, kts, kte, ii, jj, ihail)
2944
2945       IMPLICIT NONE
2946
2947 !..Sub arguments
2948       INTEGER, INTENT(IN):: kts, kte, ii, jj, ihail
2949       REAL, DIMENSION(kts:kte), INTENT(IN)::                            &
2950                       qv1d, qr1d, qs1d, qg1d, t1d, p1d
2951       REAL, DIMENSION(kts:kte), INTENT(INOUT):: dBZ
2952
2953 !..Local variables
2954       REAL, DIMENSION(kts:kte):: temp, pres, qv, rho
2955       REAL, DIMENSION(kts:kte):: rr, rs, rg
2956
2957       DOUBLE PRECISION, DIMENSION(kts:kte):: ilamr, ilams, ilamg
2958       DOUBLE PRECISION, DIMENSION(kts:kte):: N0_r, N0_s, N0_g
2959       DOUBLE PRECISION:: lamr, lams, lamg
2960       LOGICAL, DIMENSION(kts:kte):: L_qr, L_qs, L_qg
2961
2962       REAL, DIMENSION(kts:kte):: ze_rain, ze_snow, ze_graupel
2963       DOUBLE PRECISION:: fmelt_s, fmelt_g
2964
2965       INTEGER:: i, k, k_0, kbot, n
2966       LOGICAL:: melti
2967
2968       DOUBLE PRECISION:: cback, x, eta, f_d
2969       REAL, PARAMETER:: R=287.
2970       REAL, PARAMETER:: PIx=3.1415926536
2971
2972 !+---+
2973
2974       do k = kts, kte
2975          dBZ(k) = -35.0
2976       enddo
2977
2978 !+---+-----------------------------------------------------------------+
2979 !..Put column of data into local arrays.
2980 !+---+-----------------------------------------------------------------+
2981       do k = kts, kte
2982          temp(k) = t1d(k)
2983          qv(k) = MAX(1.E-10, qv1d(k))
2984          pres(k) = p1d(k)
2985          rho(k) = 0.622*pres(k)/(R*temp(k)*(qv(k)+0.622))
2986
2987          if (qr1d(k) .gt. 1.E-9) then
2988             rr(k) = qr1d(k)*rho(k)
2989             N0_r(k) = xnor
2990             lamr = (xam_r*xcrg(3)*N0_r(k)/rr(k))**(1./xcre(1))
2991             ilamr(k) = 1./lamr
2992             L_qr(k) = .true.
2993          else
2994             rr(k) = 1.E-12
2995             L_qr(k) = .false.
2996          endif
2997
2998          if (qs1d(k) .gt. 1.E-9) then
2999             rs(k) = qs1d(k)*rho(k)
3000             N0_s(k) = xnos
3001             lams = (xam_s*xcsg(3)*N0_s(k)/rs(k))**(1./xcse(1))
3002             ilams(k) = 1./lams
3003             L_qs(k) = .true.
3004          else
3005             rs(k) = 1.E-12
3006             L_qs(k) = .false.
3007          endif
3008
3009          if (qg1d(k) .gt. 1.E-9) then
3010             rg(k) = qg1d(k)*rho(k)
3011             if (ihail.eq.1) then
3012                N0_g(k) = xnoh
3013             else
3014                N0_g(k) = xnog
3015             endif
3016             lamg = (xam_g*xcgg(3)*N0_g(k)/rg(k))**(1./xcge(1))
3017             ilamg(k) = 1./lamg
3018             L_qg(k) = .true.
3019          else
3020             rg(k) = 1.E-12
3021             L_qg(k) = .false.
3022          endif
3023       enddo
3024
3025 !+---+-----------------------------------------------------------------+
3026 !..Locate K-level of start of melting (k_0 is level above).
3027 !+---+-----------------------------------------------------------------+
3028       melti = .false.
3029       k_0 = kts
3030       do k = kte-1, kts, -1
3031          if ( (temp(k).gt.273.15) .and. L_qr(k)                         &
3032                                   .and. (L_qs(k+1).or.L_qg(k+1)) ) then
3033             k_0 = MAX(k+1, k_0)
3034             melti=.true.
3035             goto 195
3036          endif
3037       enddo
3038  195  continue
3039
3040 !+---+-----------------------------------------------------------------+
3041 !..Assume Rayleigh approximation at 10 cm wavelength. Rain (all temps)
3042 !.. and non-water-coated snow and graupel when below freezing are
3043 !.. simple. Integrations of m(D)*m(D)*N(D)*dD.
3044 !+---+-----------------------------------------------------------------+
3045
3046       do k = kts, kte
3047          ze_rain(k) = 1.e-22
3048          ze_snow(k) = 1.e-22
3049          ze_graupel(k) = 1.e-22
3050          if (L_qr(k)) ze_rain(k) = N0_r(k)*xcrg(4)*ilamr(k)**xcre(4)
3051          if (L_qs(k)) ze_snow(k) = (0.176/0.93) * (6.0/PIx)*(6.0/PIx)     &
3052                                  * (xam_s/900.0)*(xam_s/900.0)          &
3053                                  * N0_s(k)*xcsg(4)*ilams(k)**xcse(4)
3054          if (L_qg(k)) ze_graupel(k) = (0.176/0.93) * (6.0/PIx)*(6.0/PIx)  &
3055                                     * (xam_g/900.0)*(xam_g/900.0)       &
3056                                     * N0_g(k)*xcgg(4)*ilamg(k)**xcge(4)
3057       enddo
3058
3059
3060 !+---+-----------------------------------------------------------------+
3061 !..Special case of melting ice (snow/graupel) particles.  Assume the
3062 !.. ice is surrounded by the liquid water.  Fraction of meltwater is
3063 !.. extremely simple based on amount found above the melting level.
3064 !.. Uses code from Uli Blahak (rayleigh_soak_wetgraupel and supporting
3065 !.. routines).
3066 !+---+-----------------------------------------------------------------+
3067
3068       if (melti .and. k_0.ge.kts+1) then
3069        do k = k_0-1, kts, -1
3070
3071 !..Reflectivity contributed by melting snow
3072           if (L_qs(k) .and. L_qs(k_0) ) then
3073            fmelt_s = MAX(0.005d0, MIN(1.0d0-rs(k)/rs(k_0), 0.99d0))
3074            eta = 0.d0
3075            lams = 1./ilams(k)
3076            do n = 1, nrbins
3077               x = xam_s * xxDs(n)**xbm_s
3078               call rayleigh_soak_wetgraupel (x,DBLE(xocms),DBLE(xobms), &
3079                     fmelt_s, melt_outside_s, m_w_0, m_i_0, lamda_radar, &
3080                     CBACK, mixingrulestring_s, matrixstring_s,          &
3081                     inclusionstring_s, hoststring_s,                    &
3082                     hostmatrixstring_s, hostinclusionstring_s)
3083               f_d = N0_s(k)*xxDs(n)**xmu_s * DEXP(-lams*xxDs(n))
3084               eta = eta + f_d * CBACK * simpson(n) * xdts(n)
3085            enddo
3086            ze_snow(k) = SNGL(lamda4 / (pi5 * K_w) * eta)
3087           endif
3088
3089 !..Reflectivity contributed by melting graupel
3090
3091           if (L_qg(k) .and. L_qg(k_0) ) then
3092            fmelt_g = MAX(0.005d0, MIN(1.0d0-rg(k)/rg(k_0), 0.99d0))
3093            eta = 0.d0
3094            lamg = 1./ilamg(k)
3095            do n = 1, nrbins
3096               x = xam_g * xxDg(n)**xbm_g
3097               call rayleigh_soak_wetgraupel (x,DBLE(xocmg),DBLE(xobmg), &
3098                     fmelt_g, melt_outside_g, m_w_0, m_i_0, lamda_radar, &
3099                     CBACK, mixingrulestring_g, matrixstring_g,          &
3100                     inclusionstring_g, hoststring_g,                    &
3101                     hostmatrixstring_g, hostinclusionstring_g)
3102               f_d = N0_g(k)*xxDg(n)**xmu_g * DEXP(-lamg*xxDg(n))
3103               eta = eta + f_d * CBACK * simpson(n) * xdtg(n)
3104            enddo
3105            ze_graupel(k) = SNGL(lamda4 / (pi5 * K_w) * eta)
3106           endif
3107
3108        enddo
3109       endif
3110
3111       do k = kte, kts, -1
3112          dBZ(k) = 10.*log10((ze_rain(k)+ze_snow(k)+ze_graupel(k))*1.d18)
3113       enddo
3114
3115       end subroutine refl10cm_gsfc
3116
3117 !+---+-----------------------------------------------------------------+
3118
3119 END MODULE  module_mp_gsfcgce