phys/module_mp_morr_two_moment.F

   1 !WRF:MODEL_LAYER:PHYSICS
   2 !
   3
   4 ! THIS MODULE CONTAINS THE TWO-MOMENT MICROPHYSICS CODE DESCRIBED BY
   5 !     MORRISON ET AL. (2009, MWR)
   6
   7 ! CHANGES FOR V3.2, RELATIVE TO MOST RECENT (BUG-FIX) CODE FOR V3.1
   8
   9 ! 1) ADDED ACCELERATED MELTING OF GRAUPEL/SNOW DUE TO COLLISION WITH RAIN, FOLLOWING LIN ET AL. (1983)
  10 ! 2) INCREASED MINIMUM LAMBDA FOR RAIN, AND ADDED RAIN DROP BREAKUP FOLLOWING MODIFIED VERSION
  11 !     OF VERLINDE AND COTTON (1993)
  12 ! 3) CHANGE MINIMUM ALLOWED MIXING RATIOS IN DRY CONDITIONS (RH < 90%), THIS IMPROVES RADAR REFLECTIIVITY
  13 !     IN LOW REFLECTIVITY REGIONS
  14 ! 4) BUG FIX TO MAXIMUM ALLOWED PARTICLE FALLSPEEDS AS A FUNCTION OF AIR DENSITY
  15 ! 5) BUG FIX TO CALCULATION OF LIQUID WATER SATURATION VAPOR PRESSURE (CHANGE IS VERY MINOR)
  16 ! 6) INCLUDE WRF CONSTANTS PER SUGGESTION OF JIMY
  17
  18 ! bug fix, 5/12/10
  19 ! 7) bug fix for saturation vapor pressure in low pressure, to avoid division by zero
  20 ! 8) include 'EP2' WRF constant for saturation mixing ratio calculation, instead of hardwire constant
  21
  22 ! CHANGES FOR V3.3
  23 ! 1) MODIFICATION FOR COUPLING WITH WRF-CHEM (PREDICTED DROPLET NUMBER CONCENTRATION) AS AN OPTION
  24 ! 2) MODIFY FALLSPEED BELOW THE LOWEST LEVEL OF PRECIPITATION, WHICH PREVENTS
  25 !      POTENTIAL FOR SPURIOUS ACCUMULATION OF PRECIPITATION DURING SUB-STEPPING FOR SEDIMENTATION
  26 ! 3) BUG FIX TO LATENT HEAT RELEASE DUE TO COLLISIONS OF CLOUD ICE WITH RAIN
  27 ! 4) CLEAN UP OF COMMENTS IN THE CODE
  28
  29 ! additional minor bug fixes and small changes, 5/30/2011
  30 ! minor revisions by A. Ackerman April 2011:
  31 ! 1) replaced kinematic with dynamic viscosity
  32 ! 2) replaced scaling by air density for cloud droplet sedimentation
  33 !    with viscosity-dependent Stokes expression
  34 ! 3) use Ikawa and Saito (1991) air-density scaling for cloud ice
  35 ! 4) corrected typo in 2nd digit of ventilation constant F2R
  36
  37 ! additional fixes:
  38 ! 5) TEMPERATURE FOR ACCELERATED MELTING DUE TO COLLIIONS OF SNOW AND GRAUPEL
  39 !    WITH RAIN SHOULD USE CELSIUS, NOT KELVIN (BUG REPORTED BY K. VAN WEVERBERG)
  40 ! 6) NPRACS IS NOT SUBTRACTED FROM SNOW NUMBER CONCENTRATION, SINCE
  41 !    DECREASE IN SNOW NUMBER IS ALREADY ACCOUNTED FOR BY NSMLTS
  42 ! 7) fix for switch for running w/o graupel/hail (cloud ice and snow only)
  43
  44 ! hm bug fix 3/16/12
  45
  46 ! 1) very minor change to limits on autoconversion source of rain number when cloud water is depleted
  47
  48 ! WRFV3.5
  49 ! hm/A. Ackerman bug fix 11/08/12
  50
  51 ! 1) for accelerated melting from collisions, should use rain mass collected by snow, not snow mass
  52 !    collected by rain
  53 ! 2) minor changes to some comments
  54 ! 3) reduction of maximum-allowed ice concentration from 10 cm-3 to 0.3
  55 !    cm-3. This was done to address the problem of excessive and persistent
  56 !    anvil cirrus produced by the scheme.
  57
  58 ! CHANGES FOR WRFV3.5.1
  59 ! 1) added output for snow+cloud ice and graupel time step and accumulated
  60 !    surface precipitation
  61 ! 2) bug fix to option w/o graupel/hail (IGRAUP = 1), include PRACI, PGSACW,
  62 !    and PGRACS as sources for snow instead of graupel/hail, bug reported by
  63 !    Hailong Wang (PNNL)
  64 ! 3) very minor fix to immersion freezing rate formulation (negligible impact)
  65 ! 4) clarifications to code comments
  66 ! 5) minor change to shedding of rain, remove limit so that the number of
  67 !    collected drops can smaller than number of shed drops
  68 ! 6) change of specific heat of liquid water from 4218 to 4187 J/kg/K
  69
  70 ! CHANGES FOR WRFV3.6.1
  71 ! 1) minor bug fix to melting of snow and graupel, an extra factor of air density (RHO) was removed
  72 !    from the calculation of PSMLT and PGMLT
  73 ! 2) redundant initialization of PSMLT (non answer-changing)
  74
  75 ! CHANGES FOR WRFV3.8.1
  76 ! 1) changes and cleanup of code comments
  77 ! 2) correction to universal gas constant (very small change)
  78
  79 ! CHANGES FOR WRFV4.3
  80 ! 1) fix to saturation vapor pressure polysvp to work at T < -80 C
  81 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  82
  83 ! THIS SCHEME IS A BULK DOUBLE-MOMENT SCHEME THAT PREDICTS MIXING
  84 ! RATIOS AND NUMBER CONCENTRATIONS OF FIVE HYDROMETEOR SPECIES:
  85 ! CLOUD DROPLETS, CLOUD (SMALL) ICE, RAIN, SNOW, AND GRAUPEL/HAIL.
  86
  87 MODULE MODULE_MP_MORR_TWO_MOMENT
  88    USE     module_wrf_error
  89    USE module_mp_radar
  90
  91 !  USE WRF PHYSICS CONSTANTS
  92   use module_model_constants, ONLY: CP, G, R => r_d, RV => r_v, EP_2
  93
  94 !  USE module_state_description
  95
  96    IMPLICIT NONE
  97
  98    REAL, PARAMETER :: PI = 3.1415926535897932384626434
  99    REAL, PARAMETER :: xxx = 0.9189385332046727417803297
 100
 101    PUBLIC  ::  MP_MORR_TWO_MOMENT
 102    PUBLIC  ::  POLYSVP
 103
 104    PRIVATE :: GAMMA, DERF1
 105    PRIVATE :: PI, xxx
 106    PRIVATE :: MORR_TWO_MOMENT_MICRO
 107
 108 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 109 ! SWITCHES FOR MICROPHYSICS SCHEME
 110 ! IACT = 1, USE POWER-LAW CCN SPECTRA, NCCN = CS^K
 111 ! IACT = 2, USE LOGNORMAL AEROSOL SIZE DIST TO DERIVE CCN SPECTRA
 112 ! IACT = 3, ACTIVATION CALCULATED IN MODULE_MIXACTIVATE
 113
 114      INTEGER, PRIVATE ::  IACT
 115
 116 ! INUM = 0, PREDICT DROPLET CONCENTRATION
 117 ! INUM = 1, ASSUME CONSTANT DROPLET CONCENTRATION
 118 ! !!!NOTE: PREDICTED DROPLET CONCENTRATION NOT AVAILABLE IN THIS VERSION
 119 ! CONTACT HUGH MORRISON (morrison@ucar.edu) FOR FURTHER INFORMATION
 120
 121      INTEGER, PRIVATE ::  INUM
 122
 123 ! FOR INUM = 1, SET CONSTANT DROPLET CONCENTRATION (CM-3)
 124      REAL, PRIVATE ::      NDCNST
 125
 126 ! SWITCH FOR LIQUID-ONLY RUN
 127 ! ILIQ = 0, INCLUDE ICE
 128 ! ILIQ = 1, LIQUID ONLY, NO ICE
 129
 130      INTEGER, PRIVATE ::  ILIQ
 131
 132 ! SWITCH FOR ICE NUCLEATION
 133 ! INUC = 0, USE FORMULA FROM RASMUSSEN ET AL. 2002 (MID-LATITUDE)
 134 !      = 1, USE MPACE OBSERVATIONS
 135
 136      INTEGER, PRIVATE ::  INUC
 137
 138 ! IBASE = 1, NEGLECT DROPLET ACTIVATION AT LATERAL CLOUD EDGES DUE TO
 139 !             UNRESOLVED ENTRAINMENT AND MIXING, ACTIVATE
 140 !             AT CLOUD BASE OR IN REGION WITH LITTLE CLOUD WATER USING
 141 !             NON-EQULIBRIUM SUPERSATURATION,
 142 !             IN CLOUD INTERIOR ACTIVATE USING EQUILIBRIUM SUPERSATURATION
 143 ! IBASE = 2, ASSUME DROPLET ACTIVATION AT LATERAL CLOUD EDGES DUE TO
 144 !             UNRESOLVED ENTRAINMENT AND MIXING DOMINATES,
 145 !             ACTIVATE DROPLETS EVERYWHERE IN THE CLOUD USING NON-EQUILIBRIUM
 146 !             SUPERSATURATION, BASED ON THE
 147 !             LOCAL SUB-GRID AND/OR GRID-SCALE VERTICAL VELOCITY
 148 !             AT THE GRID POINT
 149
 150 ! NOTE: ONLY USED FOR PREDICTED DROPLET CONCENTRATION (INUM = 0) IN NON-WRF-CHEM VERSION OF CODE
 151
 152      INTEGER, PRIVATE ::  IBASE
 153
 154 ! INCLUDE SUB-GRID VERTICAL VELOCITY IN DROPLET ACTIVATION
 155 ! ISUB = 0, INCLUDE SUB-GRID W (RECOMMENDED FOR LOWER RESOLUTION)
 156 ! ISUB = 1, EXCLUDE SUB-GRID W, ONLY USE GRID-SCALE W
 157
 158 ! NOTE: ONLY USED FOR PREDICTED DROPLET CONCENTRATION (INUM = 0) IN NON-WRF-CHEM VERSION OF CODE
 159
 160      INTEGER, PRIVATE ::  ISUB
 161
 162 ! SWITCH FOR GRAUPEL/NO GRAUPEL
 163 ! IGRAUP = 0, INCLUDE GRAUPEL
 164 ! IGRAUP = 1, NO GRAUPEL
 165
 166      INTEGER, PRIVATE ::  IGRAUP
 167
 168 ! HM ADDED NEW OPTION FOR HAIL
 169 ! SWITCH FOR HAIL/GRAUPEL
 170 ! IHAIL = 0, DENSE PRECIPITATING ICE IS GRAUPEL
 171 ! IHAIL = 1, DENSE PRECIPITATING GICE IS HAIL
 172
 173      INTEGER, PRIVATE ::  IHAIL
 174
 175 ! CLOUD MICROPHYSICS CONSTANTS
 176
 177      REAL, PRIVATE ::      AI,AC,AS,AR,AG ! 'A' PARAMETER IN FALLSPEED-DIAM RELATIONSHIP
 178      REAL, PRIVATE ::      BI,BC,BS,BR,BG ! 'B' PARAMETER IN FALLSPEED-DIAM RELATIONSHIP
 179 !     REAL, PRIVATE ::      R           ! GAS CONSTANT FOR AIR
 180 !     REAL, PRIVATE ::      RV          ! GAS CONSTANT FOR WATER VAPOR
 181 !     REAL, PRIVATE ::      CP          ! SPECIFIC HEAT AT CONSTANT PRESSURE FOR DRY AIR
 182      REAL, PRIVATE ::      RHOSU       ! STANDARD AIR DENSITY AT 850 MB
 183      REAL, PRIVATE ::      RHOW        ! DENSITY OF LIQUID WATER
 184      REAL, PRIVATE ::      RHOI        ! BULK DENSITY OF CLOUD ICE
 185      REAL, PRIVATE ::      RHOSN       ! BULK DENSITY OF SNOW
 186      REAL, PRIVATE ::      RHOG        ! BULK DENSITY OF GRAUPEL
 187      REAL, PRIVATE ::      AIMM        ! PARAMETER IN BIGG IMMERSION FREEZING
 188      REAL, PRIVATE ::      BIMM        ! PARAMETER IN BIGG IMMERSION FREEZING
 189      REAL, PRIVATE ::      ECR         ! COLLECTION EFFICIENCY BETWEEN DROPLETS/RAIN AND SNOW/RAIN
 190      REAL, PRIVATE ::      DCS         ! THRESHOLD SIZE FOR CLOUD ICE AUTOCONVERSION
 191      REAL, PRIVATE ::      MI0         ! INITIAL SIZE OF NUCLEATED CRYSTAL
 192      REAL, PRIVATE ::      MG0         ! MASS OF EMBRYO GRAUPEL
 193      REAL, PRIVATE ::      F1S         ! VENTILATION PARAMETER FOR SNOW
 194      REAL, PRIVATE ::      F2S         ! VENTILATION PARAMETER FOR SNOW
 195      REAL, PRIVATE ::      F1R         ! VENTILATION PARAMETER FOR RAIN
 196      REAL, PRIVATE ::      F2R         ! VENTILATION PARAMETER FOR RAIN
 197 !     REAL, PRIVATE ::      G           ! GRAVITATIONAL ACCELERATION
 198      REAL, PRIVATE ::      QSMALL      ! SMALLEST ALLOWED HYDROMETEOR MIXING RATIO
 199      REAL, PRIVATE ::      CI,DI,CS,DS,CG,DG ! SIZE DISTRIBUTION PARAMETERS FOR CLOUD ICE, SNOW, GRAUPEL
 200      REAL, PRIVATE ::      EII         ! COLLECTION EFFICIENCY, ICE-ICE COLLISIONS
 201      REAL, PRIVATE ::      ECI         ! COLLECTION EFFICIENCY, ICE-DROPLET COLLISIONS
 202      REAL, PRIVATE ::      RIN     ! RADIUS OF CONTACT NUCLEI (M)
 203 ! hm, add for V3.2
 204      REAL, PRIVATE ::      CPW     ! SPECIFIC HEAT OF LIQUID WATER
 205
 206 ! CCN SPECTRA FOR IACT = 1
 207
 208      REAL, PRIVATE ::      C1     ! 'C' IN NCCN = CS^K (CM-3)
 209      REAL, PRIVATE ::      K1     ! 'K' IN NCCN = CS^K
 210
 211 ! AEROSOL PARAMETERS FOR IACT = 2
 212
 213      REAL, PRIVATE ::      MW      ! MOLECULAR WEIGHT WATER (KG/MOL)
 214      REAL, PRIVATE ::      OSM     ! OSMOTIC COEFFICIENT
 215      REAL, PRIVATE ::      VI      ! NUMBER OF ION DISSOCIATED IN SOLUTION
 216      REAL, PRIVATE ::      EPSM    ! AEROSOL SOLUBLE FRACTION
 217      REAL, PRIVATE ::      RHOA    ! AEROSOL BULK DENSITY (KG/M3)
 218      REAL, PRIVATE ::      MAP     ! MOLECULAR WEIGHT AEROSOL (KG/MOL)
 219      REAL, PRIVATE ::      MA      ! MOLECULAR WEIGHT OF 'AIR' (KG/MOL)
 220      REAL, PRIVATE ::      RR      ! UNIVERSAL GAS CONSTANT
 221      REAL, PRIVATE ::      BACT    ! ACTIVATION PARAMETER
 222      REAL, PRIVATE ::      RM1     ! GEOMETRIC MEAN RADIUS, MODE 1 (M)
 223      REAL, PRIVATE ::      RM2     ! GEOMETRIC MEAN RADIUS, MODE 2 (M)
 224      REAL, PRIVATE ::      NANEW1  ! TOTAL AEROSOL CONCENTRATION, MODE 1 (M^-3)
 225      REAL, PRIVATE ::      NANEW2  ! TOTAL AEROSOL CONCENTRATION, MODE 2 (M^-3)
 226      REAL, PRIVATE ::      SIG1    ! STANDARD DEVIATION OF AEROSOL S.D., MODE 1
 227      REAL, PRIVATE ::      SIG2    ! STANDARD DEVIATION OF AEROSOL S.D., MODE 2
 228      REAL, PRIVATE ::      F11     ! CORRECTION FACTOR FOR ACTIVATION, MODE 1
 229      REAL, PRIVATE ::      F12     ! CORRECTION FACTOR FOR ACTIVATION, MODE 1
 230      REAL, PRIVATE ::      F21     ! CORRECTION FACTOR FOR ACTIVATION, MODE 2
 231      REAL, PRIVATE ::      F22     ! CORRECTION FACTOR FOR ACTIVATION, MODE 2
 232      REAL, PRIVATE ::      MMULT   ! MASS OF SPLINTERED ICE PARTICLE
 233      REAL, PRIVATE ::      LAMMAXI,LAMMINI,LAMMAXR,LAMMINR,LAMMAXS,LAMMINS,LAMMAXG,LAMMING
 234
 235 ! CONSTANTS TO IMPROVE EFFICIENCY
 236
 237      REAL, PRIVATE :: CONS1,CONS2,CONS3,CONS4,CONS5,CONS6,CONS7,CONS8,CONS9,CONS10
 238      REAL, PRIVATE :: CONS11,CONS12,CONS13,CONS14,CONS15,CONS16,CONS17,CONS18,CONS19,CONS20
 239      REAL, PRIVATE :: CONS21,CONS22,CONS23,CONS24,CONS25,CONS26,CONS27,CONS28,CONS29,CONS30
 240      REAL, PRIVATE :: CONS31,CONS32,CONS33,CONS34,CONS35,CONS36,CONS37,CONS38,CONS39,CONS40
 241      REAL, PRIVATE :: CONS41
 242
 243
 244 CONTAINS
 245
 246 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 247 SUBROUTINE MORR_TWO_MOMENT_INIT(morr_rimed_ice) ! RAS
 248 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 249 ! THIS SUBROUTINE INITIALIZES ALL PHYSICAL CONSTANTS AMND PARAMETERS
 250 ! NEEDED BY THE MICROPHYSICS SCHEME.
 251 ! NEEDS TO BE CALLED AT FIRST TIME STEP, PRIOR TO CALL TO MAIN MICROPHYSICS INTERFACE
 252 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 253
 254       IMPLICIT NONE
 255
 256       INTEGER, INTENT(IN):: morr_rimed_ice ! RAS
 257
 258       integer n,i
 259
 260 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 261 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 262
 263 ! THE FOLLOWING PARAMETERS ARE USER-DEFINED SWITCHES AND NEED TO BE
 264 ! SET PRIOR TO CODE COMPILATION
 265
 266 ! INUM IS AUTOMATICALLY SET TO 0 FOR WRF-CHEM BELOW,
 267 ! ALLOWING PREDICTION OF DROPLET CONCENTRATION
 268 ! THUS, THIS PARAMETER SHOULD NOT BE CHANGED HERE
 269 ! AND SHOULD BE LEFT TO 1
 270
 271       INUM = 1
 272
 273 ! SET CONSTANT DROPLET CONCENTRATION (UNITS OF CM-3)
 274 ! IF NO COUPLING WITH WRF-CHEM
 275
 276       NDCNST = 250.
 277
 278 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 279 ! NOTE, THE FOLLOWING OPTIONS RELATED TO DROPLET ACTIVATION
 280 ! (IACT, IBASE, ISUB) ARE NOT AVAILABLE IN CURRENT VERSION
 281 ! FOR WRF-CHEM, DROPLET ACTIVATION IS PERFORMED
 282 ! IN 'MIX_ACTIVATE', NOT IN MICROPHYSICS SCHEME
 283
 284
 285 ! IACT = 1, USE POWER-LAW CCN SPECTRA, NCCN = CS^K
 286 ! IACT = 2, USE LOGNORMAL AEROSOL SIZE DIST TO DERIVE CCN SPECTRA
 287
 288       IACT = 2
 289
 290 ! IBASE = 1, NEGLECT DROPLET ACTIVATION AT LATERAL CLOUD EDGES DUE TO
 291 !             UNRESOLVED ENTRAINMENT AND MIXING, ACTIVATE
 292 !             AT CLOUD BASE OR IN REGION WITH LITTLE CLOUD WATER USING
 293 !             NON-EQULIBRIUM SUPERSATURATION ASSUMING NO INITIAL CLOUD WATER,
 294 !             IN CLOUD INTERIOR ACTIVATE USING EQUILIBRIUM SUPERSATURATION
 295 ! IBASE = 2, ASSUME DROPLET ACTIVATION AT LATERAL CLOUD EDGES DUE TO
 296 !             UNRESOLVED ENTRAINMENT AND MIXING DOMINATES,
 297 !             ACTIVATE DROPLETS EVERYWHERE IN THE CLOUD USING NON-EQUILIBRIUM
 298 !             SUPERSATURATION ASSUMING NO INITIAL CLOUD WATER, BASED ON THE
 299 !             LOCAL SUB-GRID AND/OR GRID-SCALE VERTICAL VELOCITY
 300 !             AT THE GRID POINT
 301
 302 ! NOTE: ONLY USED FOR PREDICTED DROPLET CONCENTRATION (INUM = 0)
 303
 304       IBASE = 2
 305
 306 ! INCLUDE SUB-GRID VERTICAL VELOCITY (standard deviation of w) IN DROPLET ACTIVATION
 307 ! ISUB = 0, INCLUDE SUB-GRID W (RECOMMENDED FOR LOWER RESOLUTION)
 308 ! currently, sub-grid w is constant of 0.5 m/s (not coupled with PBL/turbulence scheme)
 309 ! ISUB = 1, EXCLUDE SUB-GRID W, ONLY USE GRID-SCALE W
 310
 311 ! NOTE: ONLY USED FOR PREDICTED DROPLET CONCENTRATION (INUM = 0)
 312
 313       ISUB = 0
 314 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 315
 316
 317 ! SWITCH FOR LIQUID-ONLY RUN
 318 ! ILIQ = 0, INCLUDE ICE
 319 ! ILIQ = 1, LIQUID ONLY, NO ICE
 320
 321       ILIQ = 0
 322
 323 ! SWITCH FOR ICE NUCLEATION
 324 ! INUC = 0, USE FORMULA FROM RASMUSSEN ET AL. 2002 (MID-LATITUDE)
 325 !      = 1, USE MPACE OBSERVATIONS (ARCTIC ONLY)
 326
 327       INUC = 0
 328
 329 ! SWITCH FOR GRAUPEL/HAIL NO GRAUPEL/HAIL
 330 ! IGRAUP = 0, INCLUDE GRAUPEL/HAIL
 331 ! IGRAUP = 1, NO GRAUPEL/HAIL
 332
 333       IGRAUP = 0
 334
 335 ! HM ADDED 11/7/07
 336 ! SWITCH FOR HAIL/GRAUPEL
 337 ! IHAIL = 0, DENSE PRECIPITATING ICE IS GRAUPEL
 338 ! IHAIL = 1, DENSE PRECIPITATING ICE IS HAIL
 339 ! NOTE ---> RECOMMEND IHAIL = 1 FOR CONTINENTAL DEEP CONVECTION
 340
 341       !IHAIL = 0 !changed to namelist option (morr_rimed_ice) by RAS
 342       ! Check if namelist option is feasible, otherwise default to graupel - RAS
 343       IF (morr_rimed_ice .eq. 1) THEN
 344          IHAIL = 1
 345       ELSE
 346          IHAIL = 0
 347       ENDIF
 348
 349 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 350 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 351 ! SET PHYSICAL CONSTANTS
 352
 353 ! FALLSPEED PARAMETERS (V=AD^B)
 354          AI = 700.
 355          AC = 3.E7
 356          AS = 11.72
 357          AR = 841.99667
 358          BI = 1.
 359          BC = 2.
 360          BS = 0.41
 361          BR = 0.8
 362          IF (IHAIL.EQ.0) THEN
 363          AG = 19.3
 364          BG = 0.37
 365          ELSE ! (MATSUN AND HUGGINS 1980)
 366          AG = 114.5
 367          BG = 0.5
 368          END IF
 369
 370 ! CONSTANTS AND PARAMETERS
 371 !         R = 287.15
 372 !         RV = 461.5
 373 !         CP = 1005.
 374          RHOSU = 85000./(287.15*273.15)
 375          RHOW = 997.
 376          RHOI = 500.
 377          RHOSN = 100.
 378          IF (IHAIL.EQ.0) THEN
 379          RHOG = 400.
 380          ELSE
 381          RHOG = 900.
 382          END IF
 383          AIMM = 0.66
 384          BIMM = 100.
 385          ECR = 1.
 386          DCS = 125.E-6
 387          MI0 = 4./3.*PI*RHOI*(10.E-6)**3
 388          MG0 = 1.6E-10
 389          F1S = 0.86
 390          F2S = 0.28
 391          F1R = 0.78
 392 !         F2R = 0.32
 393 ! fix 053011
 394          F2R = 0.308
 395 !         G = 9.806
 396          QSMALL = 1.E-14
 397          EII = 0.1
 398          ECI = 0.7
 399 ! HM, ADD FOR V3.2
 400 ! hm, 7/23/13
 401 !         CPW = 4218.
 402          CPW = 4187.
 403
 404 ! SIZE DISTRIBUTION PARAMETERS
 405
 406          CI = RHOI*PI/6.
 407          DI = 3.
 408          CS = RHOSN*PI/6.
 409          DS = 3.
 410          CG = RHOG*PI/6.
 411          DG = 3.
 412
 413 ! RADIUS OF CONTACT NUCLEI
 414          RIN = 0.1E-6
 415
 416          MMULT = 4./3.*PI*RHOI*(5.E-6)**3
 417
 418 ! SIZE LIMITS FOR LAMBDA
 419
 420          LAMMAXI = 1./1.E-6
 421          LAMMINI = 1./(2.*DCS+100.E-6)
 422          LAMMAXR = 1./20.E-6
 423 !         LAMMINR = 1./500.E-6
 424          LAMMINR = 1./2800.E-6
 425
 426          LAMMAXS = 1./10.E-6
 427          LAMMINS = 1./2000.E-6
 428          LAMMAXG = 1./20.E-6
 429          LAMMING = 1./2000.E-6
 430
 431 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 432 ! note: these parameters only used by the non-wrf-chem version of the
 433 !       scheme with predicted droplet number
 434
 435 ! CCN SPECTRA FOR IACT = 1
 436
 437 ! MARITIME
 438 ! MODIFIED FROM RASMUSSEN ET AL. 2002
 439 ! NCCN = C*S^K, NCCN IS IN CM-3, S IS SUPERSATURATION RATIO IN %
 440
 441               K1 = 0.4
 442               C1 = 120.
 443
 444 ! CONTINENTAL
 445
 446 !              K1 = 0.5
 447 !              C1 = 1000.
 448
 449 ! AEROSOL ACTIVATION PARAMETERS FOR IACT = 2
 450 ! PARAMETERS CURRENTLY SET FOR AMMONIUM SULFATE
 451
 452          MW = 0.018
 453          OSM = 1.
 454          VI = 3.
 455          EPSM = 0.7
 456          RHOA = 1777.
 457          MAP = 0.132
 458          MA = 0.0284
 459 ! hm fix 6/23/16
 460 !         RR = 8.3187
 461          RR = 8.3145
 462          BACT = VI*OSM*EPSM*MW*RHOA/(MAP*RHOW)
 463
 464 ! AEROSOL SIZE DISTRIBUTION PARAMETERS CURRENTLY SET FOR MPACE
 465 ! (see morrison et al. 2007, JGR)
 466 ! MODE 1
 467
 468          RM1 = 0.052E-6
 469          SIG1 = 2.04
 470          NANEW1 = 72.2E6
 471          F11 = 0.5*EXP(2.5*(LOG(SIG1))**2)
 472          F21 = 1.+0.25*LOG(SIG1)
 473
 474 ! MODE 2
 475
 476          RM2 = 1.3E-6
 477          SIG2 = 2.5
 478          NANEW2 = 1.8E6
 479          F12 = 0.5*EXP(2.5*(LOG(SIG2))**2)
 480          F22 = 1.+0.25*LOG(SIG2)
 481 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 482
 483 ! CONSTANTS FOR EFFICIENCY
 484
 485          CONS1=GAMMA(1.+DS)*CS
 486          CONS2=GAMMA(1.+DG)*CG
 487          CONS3=GAMMA(4.+BS)/6.
 488          CONS4=GAMMA(4.+BR)/6.
 489          CONS5=GAMMA(1.+BS)
 490          CONS6=GAMMA(1.+BR)
 491          CONS7=GAMMA(4.+BG)/6.
 492          CONS8=GAMMA(1.+BG)
 493          CONS9=GAMMA(5./2.+BR/2.)
 494          CONS10=GAMMA(5./2.+BS/2.)
 495          CONS11=GAMMA(5./2.+BG/2.)
 496          CONS12=GAMMA(1.+DI)*CI
 497          CONS13=GAMMA(BS+3.)*PI/4.*ECI
 498          CONS14=GAMMA(BG+3.)*PI/4.*ECI
 499          CONS15=-1108.*EII*PI**((1.-BS)/3.)*RHOSN**((-2.-BS)/3.)/(4.*720.)
 500          CONS16=GAMMA(BI+3.)*PI/4.*ECI
 501          CONS17=4.*2.*3.*RHOSU*PI*ECI*ECI*GAMMA(2.*BS+2.)/(8.*(RHOG-RHOSN))
 502          CONS18=RHOSN*RHOSN
 503          CONS19=RHOW*RHOW
 504          CONS20=20.*PI*PI*RHOW*BIMM
 505          CONS21=4./(DCS*RHOI)
 506          CONS22=PI*RHOI*DCS**3/6.
 507          CONS23=PI/4.*EII*GAMMA(BS+3.)
 508          CONS24=PI/4.*ECR*GAMMA(BR+3.)
 509          CONS25=PI*PI/24.*RHOW*ECR*GAMMA(BR+6.)
 510          CONS26=PI/6.*RHOW
 511          CONS27=GAMMA(1.+BI)
 512          CONS28=GAMMA(4.+BI)/6.
 513          CONS29=4./3.*PI*RHOW*(25.E-6)**3
 514          CONS30=4./3.*PI*RHOW
 515          CONS31=PI*PI*ECR*RHOSN
 516          CONS32=PI/2.*ECR
 517          CONS33=PI*PI*ECR*RHOG
 518          CONS34=5./2.+BR/2.
 519          CONS35=5./2.+BS/2.
 520          CONS36=5./2.+BG/2.
 521          CONS37=4.*PI*1.38E-23/(6.*PI*RIN)
 522          CONS38=PI*PI/3.*RHOW
 523          CONS39=PI*PI/36.*RHOW*BIMM
 524          CONS40=PI/6.*BIMM
 525          CONS41=PI*PI*ECR*RHOW
 526
 527 !+---+-----------------------------------------------------------------+
 528 !..Set these variables needed for computing radar reflectivity.  These
 529 !.. get used within radar_init to create other variables used in the
 530 !.. radar module.
 531
 532          xam_r = PI*RHOW/6.
 533          xbm_r = 3.
 534          xmu_r = 0.
 535          xam_s = CS
 536          xbm_s = DS
 537          xmu_s = 0.
 538          xam_g = CG
 539          xbm_g = DG
 540          xmu_g = 0.
 541
 542          call radar_init
 543 !+---+-----------------------------------------------------------------+
 544
 545
 546 END SUBROUTINE MORR_TWO_MOMENT_INIT
 547
 548 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 549 ! THIS SUBROUTINE IS MAIN INTERFACE WITH THE TWO-MOMENT MICROPHYSICS SCHEME
 550 ! THIS INTERFACE TAKES IN 3D VARIABLES FROM DRIVER MODEL, CONVERTS TO 1D FOR
 551 ! CALL TO THE MAIN MICROPHYSICS SUBROUTINE (SUBROUTINE MORR_TWO_MOMENT_MICRO)
 552 ! WHICH OPERATES ON 1D VERTICAL COLUMNS.
 553 ! 1D VARIABLES FROM THE MAIN MICROPHYSICS SUBROUTINE ARE THEN REASSIGNED BACK TO 3D FOR OUTPUT
 554 ! BACK TO DRIVER MODEL USING THIS INTERFACE.
 555 ! MICROPHYSICS TENDENCIES ARE ADDED TO VARIABLES HERE BEFORE BEING PASSED BACK TO DRIVER MODEL.
 556
 557 ! THIS CODE WAS WRITTEN BY HUGH MORRISON (NCAR) AND SLAVA TATARSKII (GEORGIA TECH).
 558
 559 ! FOR QUESTIONS, CONTACT: HUGH MORRISON, E-MAIL: MORRISON@UCAR.EDU, PHONE:303-497-8916
 560
 561 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 562
 563 SUBROUTINE MP_MORR_TWO_MOMENT(ITIMESTEP,                       &
 564                 TH, QV, QC, QR, QI, QS, QG, NI, NS, NR, NG, &
 565                 RHO, PII, P, DT_IN, DZ, HT, W,          &
 566                 RAINNC, RAINNCV, SR,                    &
 567                 SNOWNC,SNOWNCV,GRAUPELNC,GRAUPELNCV,    & ! hm added 7/13/13
 568                 refl_10cm, diagflag, do_radar_ref,      & ! GT added for reflectivity calcs
 569                 qrcuten, qscuten, qicuten               & ! hm added
 570                ,F_QNDROP, qndrop                        & ! hm added, wrf-chem
 571                ,IDS,IDE, JDS,JDE, KDS,KDE               & ! domain dims
 572                ,IMS,IME, JMS,JME, KMS,KME               & ! memory dims
 573                ,ITS,ITE, JTS,JTE, KTS,KTE               & ! tile   dims            )
 574 !jdf               ,C2PREC3D,CSED3D,ISED3D,SSED3D,GSED3D,RSED3D & ! HM ADD, WRF-CHEM
 575                ,wetscav_on, rainprod, evapprod                      &
 576                    ,QLSINK,PRECR,PRECI,PRECS,PRECG &        ! HM ADD, WRF-CHEM
 577                                             )
 578
 579 ! QV - water vapor mixing ratio (kg/kg)
 580 ! QC - cloud water mixing ratio (kg/kg)
 581 ! QR - rain water mixing ratio (kg/kg)
 582 ! QI - cloud ice mixing ratio (kg/kg)
 583 ! QS - snow mixing ratio (kg/kg)
 584 ! QG - graupel mixing ratio (KG/KG)
 585 ! NI - cloud ice number concentration (1/kg)
 586 ! NS - Snow Number concentration (1/kg)
 587 ! NR - Rain Number concentration (1/kg)
 588 ! NG - Graupel number concentration (1/kg)
 589 ! NOTE: RHO AND HT NOT USED BY THIS SCHEME AND DO NOT NEED TO BE PASSED INTO SCHEME!!!!
 590 ! P - AIR PRESSURE (PA)
 591 ! W - VERTICAL AIR VELOCITY (M/S)
 592 ! TH - POTENTIAL TEMPERATURE (K)
 593 ! PII - exner function - used to convert potential temp to temp
 594 ! DZ - difference in height over interface (m)
 595 ! DT_IN - model time step (sec)
 596 ! ITIMESTEP - time step counter
 597 ! RAINNC - accumulated grid-scale precipitation (mm)
 598 ! RAINNCV - one time step grid scale precipitation (mm/time step)
 599 ! SNOWNC - accumulated grid-scale snow plus cloud ice (mm)
 600 ! SNOWNCV - one time step grid scale snow plus cloud ice (mm/time step)
 601 ! GRAUPELNC - accumulated grid-scale graupel (mm)
 602 ! GRAUPELNCV - one time step grid scale graupel (mm/time step)
 603 ! SR - one time step mass ratio of snow to total precip
 604 ! qrcuten, rain tendency from parameterized cumulus convection
 605 ! qscuten, snow tendency from parameterized cumulus convection
 606 ! qicuten, cloud ice tendency from parameterized cumulus convection
 607
 608 ! variables below currently not in use, not coupled to PBL or radiation codes
 609 ! TKE - turbulence kinetic energy (m^2 s-2), NEEDED FOR DROPLET ACTIVATION (SEE CODE BELOW)
 610 ! NCTEND - droplet concentration tendency from pbl (kg-1 s-1)
 611 ! NCTEND - CLOUD ICE concentration tendency from pbl (kg-1 s-1)
 612 ! KZH - heat eddy diffusion coefficient from YSU scheme (M^2 S-1), NEEDED FOR DROPLET ACTIVATION (SEE CODE BELOW)
 613 ! EFFCS - CLOUD DROPLET EFFECTIVE RADIUS OUTPUT TO RADIATION CODE (micron)
 614 ! EFFIS - CLOUD DROPLET EFFECTIVE RADIUS OUTPUT TO RADIATION CODE (micron)
 615 ! HM, ADDED FOR WRF-CHEM COUPLING
 616 ! QLSINK - TENDENCY OF CLOUD WATER TO RAIN, SNOW, GRAUPEL (KG/KG/S)
 617 ! CSED,ISED,SSED,GSED,RSED - SEDIMENTATION FLUXES (KG/M^2/S) FOR CLOUD WATER, ICE, SNOW, GRAUPEL, RAIN
 618 ! PRECI,PRECS,PRECG,PRECR - SEDIMENTATION FLUXES (KG/M^2/S) FOR ICE, SNOW, GRAUPEL, RAIN
 619
 620 ! rainprod - total tendency of conversion of cloud water/ice and graupel to rain (kg kg-1 s-1)
 621 ! evapprod - tendency of evaporation of rain (kg kg-1 s-1)
 622
 623 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 624
 625 ! reflectivity currently not included!!!!
 626 ! REFL_10CM - CALCULATED RADAR REFLECTIVITY AT 10 CM (DBZ)
 627 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 628
 629 ! EFFC - DROPLET EFFECTIVE RADIUS (MICRON)
 630 ! EFFR - RAIN EFFECTIVE RADIUS (MICRON)
 631 ! EFFS - SNOW EFFECTIVE RADIUS (MICRON)
 632 ! EFFI - CLOUD ICE EFFECTIVE RADIUS (MICRON)
 633
 634 ! ADDITIONAL OUTPUT FROM MICRO - SEDIMENTATION TENDENCIES, NEEDED FOR LIQUID-ICE STATIC ENERGY
 635
 636 ! QGSTEN - GRAUPEL SEDIMENTATION TEND (KG/KG/S)
 637 ! QRSTEN - RAIN SEDIMENTATION TEND (KG/KG/S)
 638 ! QISTEN - CLOUD ICE SEDIMENTATION TEND (KG/KG/S)
 639 ! QNISTEN - SNOW SEDIMENTATION TEND (KG/KG/S)
 640 ! QCSTEN - CLOUD WATER SEDIMENTATION TEND (KG/KG/S)
 641
 642 ! WVAR - STANDARD DEVIATION OF SUB-GRID VERTICAL VELOCITY (M/S)
 643
 644    IMPLICIT NONE
 645
 646    INTEGER,      INTENT(IN   )    ::   ids, ide, jds, jde, kds, kde , &
 647                                        ims, ime, jms, jme, kms, kme , &
 648                                        its, ite, jts, jte, kts, kte
 649 ! Temporary changed from INOUT to IN
 650
 651    REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT):: &
 652                           qv, qc, qr, qi, qs, qg, ni, ns, nr, TH, NG
 653 !jdf                      qndrop ! hm added, wrf-chem
 654    REAL, DIMENSION(ims:ime, kms:kme, jms:jme), optional,INTENT(INOUT):: qndrop
 655 !jdf  REAL, DIMENSION(ims:ime, kms:kme, jms:jme),INTENT(INOUT):: CSED3D, &
 656    REAL, DIMENSION(ims:ime, kms:kme, jms:jme), optional,INTENT(INOUT):: QLSINK, &
 657                           rainprod, evapprod, &
 658                           PRECI,PRECS,PRECG,PRECR ! HM, WRF-CHEM
 659 !, effcs, effis
 660
 661    REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(IN):: &
 662                           pii, p, dz, rho, w !, tke, nctend, nitend,kzh
 663    REAL, INTENT(IN):: dt_in
 664    INTEGER, INTENT(IN):: ITIMESTEP
 665
 666    REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT):: &
 667                           RAINNC, RAINNCV, SR, &
 668 ! hm added 7/13/13
 669                           SNOWNC,SNOWNCV,GRAUPELNC,GRAUPELNCV
 670
 671    REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT)::       &  ! GT
 672                           refl_10cm
 673
 674    REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN) ::       ht
 675
 676    LOGICAL, optional, INTENT(IN) :: wetscav_on
 677
 678    ! LOCAL VARIABLES
 679
 680    REAL, DIMENSION(its:ite, kts:kte, jts:jte)::                     &
 681                       effi, effs, effr, EFFG
 682
 683    REAL, DIMENSION(its:ite, kts:kte, jts:jte)::                     &
 684                       T, WVAR, EFFC
 685
 686    REAL, DIMENSION(kts:kte) ::                                                                &
 687                             QC_TEND1D, QI_TEND1D, QNI_TEND1D, QR_TEND1D,                      &
 688                             NI_TEND1D, NS_TEND1D, NR_TEND1D,                                  &
 689                             QC1D, QI1D, QR1D,NI1D, NS1D, NR1D, QS1D,                          &
 690                             T_TEND1D,QV_TEND1D, T1D, QV1D, P1D, W1D, WVAR1D,         &
 691                             EFFC1D, EFFI1D, EFFS1D, EFFR1D,DZ1D,   &
 692    ! HM ADD GRAUPEL
 693                             QG_TEND1D, NG_TEND1D, QG1D, NG1D, EFFG1D, &
 694
 695 ! ADD SEDIMENTATION TENDENCIES (UNITS OF KG/KG/S)
 696                             QGSTEN,QRSTEN, QISTEN, QNISTEN, QCSTEN, &
 697 ! ADD CUMULUS TENDENCIES
 698                             QRCU1D, QSCU1D, QICU1D
 699
 700 ! add cumulus tendencies
 701
 702    REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(IN):: &
 703       qrcuten, qscuten, qicuten
 704
 705   LOGICAL, INTENT(IN), OPTIONAL ::                F_QNDROP  ! wrf-chem
 706   LOGICAL :: flag_qndrop  ! wrf-chem
 707   integer :: iinum ! wrf-chem
 708
 709 ! wrf-chem
 710    REAL, DIMENSION(kts:kte) :: nc1d, nc_tend1d,C2PREC,CSED,ISED,SSED,GSED,RSED
 711    REAL, DIMENSION(kts:kte) :: rainprod1d, evapprod1d
 712 ! HM add reflectivity
 713    REAL, DIMENSION(kts:kte) :: dBZ
 714
 715    REAL PRECPRT1D, SNOWRT1D, SNOWPRT1D, GRPLPRT1D ! hm added 7/13/13
 716
 717    INTEGER I,K,J
 718
 719    REAL DT
 720
 721    LOGICAL, OPTIONAL, INTENT(IN) :: diagflag
 722    INTEGER, OPTIONAL, INTENT(IN) :: do_radar_ref
 723
 724    LOGICAL :: has_wetscav
 725
 726 ! below for wrf-chem
 727    flag_qndrop = .false.
 728    IF ( PRESENT ( f_qndrop ) ) flag_qndrop = f_qndrop
 729 !!!!!!!!!!!!!!!!!!!!!!
 730
 731    IF( PRESENT( wetscav_on ) ) THEN
 732      has_wetscav = wetscav_on
 733    ELSE
 734      has_wetscav = .false.
 735    ENDIF
 736
 737    ! Initialize tendencies (all set to 0) and transfer
 738    ! array to local variables
 739    DT = DT_IN
 740
 741    DO I=ITS,ITE
 742    DO J=JTS,JTE
 743    DO K=KTS,KTE
 744        T(I,K,J)        = TH(i,k,j)*PII(i,k,j)
 745
 746 ! NOTE: WVAR NOT CURRENTLY USED IN CODE !!!!!!!!!!
 747 ! currently assign wvar to 0.5 m/s (not coupled with PBL scheme)
 748
 749        WVAR(I,K,J)     = 0.5
 750
 751 ! currently mixing of number concentrations also is neglected (not coupled with PBL schemes)
 752
 753    END DO
 754    END DO
 755    END DO
 756
 757    do i=its,ite      ! i loop (east-west)
 758    do j=jts,jte      ! j loop (north-south)
 759    !
 760    ! Transfer 3D arrays into 1D for microphysical calculations
 761    !
 762
 763 ! hm , initialize 1d tendency arrays to zero
 764
 765       do k=kts,kte   ! k loop (vertical)
 766
 767           QC_TEND1D(k)  = 0.
 768           QI_TEND1D(k)  = 0.
 769           QNI_TEND1D(k) = 0.
 770           QR_TEND1D(k)  = 0.
 771           NI_TEND1D(k)  = 0.
 772           NS_TEND1D(k)  = 0.
 773           NR_TEND1D(k)  = 0.
 774           T_TEND1D(k)   = 0.
 775           QV_TEND1D(k)  = 0.
 776           nc_tend1d(k) = 0. ! wrf-chem
 777
 778           QC1D(k)       = QC(i,k,j)
 779           QI1D(k)       = QI(i,k,j)
 780           QS1D(k)       = QS(i,k,j)
 781           QR1D(k)       = QR(i,k,j)
 782
 783           NI1D(k)       = NI(i,k,j)
 784
 785           NS1D(k)       = NS(i,k,j)
 786           NR1D(k)       = NR(i,k,j)
 787 ! HM ADD GRAUPEL
 788           QG1D(K)       = QG(I,K,j)
 789           NG1D(K)       = NG(I,K,j)
 790           QG_TEND1D(K)  = 0.
 791           NG_TEND1D(K)  = 0.
 792
 793           T1D(k)        = T(i,k,j)
 794           QV1D(k)       = QV(i,k,j)
 795           P1D(k)        = P(i,k,j)
 796           DZ1D(k)       = DZ(i,k,j)
 797           W1D(k)        = W(i,k,j)
 798           WVAR1D(k)     = WVAR(i,k,j)
 799 ! add cumulus tendencies, already decoupled
 800           qrcu1d(k)     = qrcuten(i,k,j)
 801           qscu1d(k)     = qscuten(i,k,j)
 802           qicu1d(k)     = qicuten(i,k,j)
 803       end do  !jdf added this
 804 ! below for wrf-chem
 805    IF (flag_qndrop .AND. PRESENT( qndrop )) THEN
 806       iact = 3
 807       DO k = kts, kte
 808          nc1d(k)=qndrop(i,k,j)
 809          iinum=0
 810       ENDDO
 811    ELSE
 812       DO k = kts, kte
 813          nc1d(k)=0. ! temporary placeholder, set to constant in microphysics subroutine
 814          iinum=1
 815       ENDDO
 816    ENDIF
 817
 818 !jdf  end do
 819
 820       call MORR_TWO_MOMENT_MICRO(QC_TEND1D, QI_TEND1D, QNI_TEND1D, QR_TEND1D,            &
 821        NI_TEND1D, NS_TEND1D, NR_TEND1D,                                                  &
 822        QC1D, QI1D, QS1D, QR1D,NI1D, NS1D, NR1D,                                          &
 823        T_TEND1D,QV_TEND1D, T1D, QV1D, P1D, DZ1D, W1D, WVAR1D,                   &
 824        PRECPRT1D,SNOWRT1D,                                                               &
 825        SNOWPRT1D,GRPLPRT1D,                 & ! hm added 7/13/13
 826        EFFC1D,EFFI1D,EFFS1D,EFFR1D,DT,                                                   &
 827                                             IMS,IME, JMS,JME, KMS,KME,                   &
 828                                             ITS,ITE, JTS,JTE, KTS,KTE,                   & ! HM ADD GRAUPEL
 829                                     QG_TEND1D,NG_TEND1D,QG1D,NG1D,EFFG1D, &
 830                                     qrcu1d, qscu1d, qicu1d, &
 831 ! ADD SEDIMENTATION TENDENCIES
 832                                   QGSTEN,QRSTEN,QISTEN,QNISTEN,QCSTEN, &
 833                                   nc1d, nc_tend1d, iinum, C2PREC,CSED,ISED,SSED,GSED,RSED & !wrf-chem
 834 #if (WRF_CHEM == 1)
 835                                   ,has_wetscav,rainprod1d, evapprod1d & !wrf-chem
 836 #endif
 837                        )
 838
 839    !
 840    ! Transfer 1D arrays back into 3D arrays
 841    !
 842       do k=kts,kte
 843
 844 ! hm, add tendencies to update global variables
 845 ! HM, TENDENCIES FOR Q AND N NOW ADDED IN M2005MICRO, SO WE
 846 ! ONLY NEED TO TRANSFER 1D VARIABLES BACK TO 3D
 847
 848           QC(i,k,j)        = QC1D(k)
 849           QI(i,k,j)        = QI1D(k)
 850           QS(i,k,j)        = QS1D(k)
 851           QR(i,k,j)        = QR1D(k)
 852           NI(i,k,j)        = NI1D(k)
 853           NS(i,k,j)        = NS1D(k)
 854           NR(i,k,j)        = NR1D(k)
 855           QG(I,K,j)        = QG1D(K)
 856           NG(I,K,j)        = NG1D(K)
 857
 858           T(i,k,j)         = T1D(k)
 859           TH(I,K,J)        = T(i,k,j)/PII(i,k,j) ! CONVERT TEMP BACK TO POTENTIAL TEMP
 860           QV(i,k,j)        = QV1D(k)
 861
 862           EFFC(i,k,j)      = EFFC1D(k)
 863           EFFI(i,k,j)      = EFFI1D(k)
 864           EFFS(i,k,j)      = EFFS1D(k)
 865           EFFR(i,k,j)      = EFFR1D(k)
 866           EFFG(I,K,j)      = EFFG1D(K)
 867
 868 ! wrf-chem
 869           IF (flag_qndrop .AND. PRESENT( qndrop )) THEN
 870              qndrop(i,k,j) = nc1d(k)
 871 !jdf         CSED3D(I,K,J) = CSED(K)
 872           END IF
 873           IF ( PRESENT( QLSINK ) ) THEN
 874              if(qc(i,k,j)>1.e-10) then
 875                 QLSINK(I,K,J)  = C2PREC(K)/QC(I,K,J)
 876              else
 877                 QLSINK(I,K,J)  = 0.0
 878              endif
 879           END IF
 880           IF ( PRESENT( PRECR ) ) PRECR(I,K,J) = RSED(K)
 881           IF ( PRESENT( PRECI ) ) PRECI(I,K,J) = ISED(K)
 882           IF ( PRESENT( PRECS ) ) PRECS(I,K,J) = SSED(K)
 883           IF ( PRESENT( PRECG ) ) PRECG(I,K,J) = GSED(K)
 884 ! EFFECTIVE RADIUS FOR RADIATION CODE (currently not coupled)
 885 ! HM, ADD LIMIT TO PREVENT BLOWING UP OPTICAL PROPERTIES, 8/18/07
 886 !          EFFCS(I,K,J)     = MIN(EFFC(I,K,J),50.)
 887 !          EFFCS(I,K,J)     = MAX(EFFCS(I,K,J),1.)
 888 !          EFFIS(I,K,J)     = MIN(EFFI(I,K,J),130.)
 889 !          EFFIS(I,K,J)     = MAX(EFFIS(I,K,J),13.)
 890
 891 #if ( WRF_CHEM == 1)
 892            IF ( has_wetscav ) THEN
 893              IF ( PRESENT( rainprod ) ) rainprod(i,k,j) = rainprod1d(k)
 894              IF ( PRESENT( evapprod ) ) evapprod(i,k,j) = evapprod1d(k)
 895            ENDIF
 896 #endif
 897
 898       end do
 899
 900 ! hm modified so that m2005 precip variables correctly match wrf precip variables
 901       RAINNC(i,j) = RAINNC(I,J)+PRECPRT1D
 902       RAINNCV(i,j) = PRECPRT1D
 903 ! hm, added 7/13/13
 904       SNOWNC(i,j) = SNOWNC(I,J)+SNOWPRT1D
 905       SNOWNCV(i,j) = SNOWPRT1D
 906       GRAUPELNC(i,j) = GRAUPELNC(I,J)+GRPLPRT1D
 907       GRAUPELNCV(i,j) = GRPLPRT1D
 908       SR(i,j) = SNOWRT1D/(PRECPRT1D+1.E-12)
 909
 910 !+---+-----------------------------------------------------------------+
 911          IF ( PRESENT (diagflag) ) THEN
 912          if (diagflag .and. do_radar_ref == 1) then
 913           call refl10cm_hm (qv1d, qr1d, nr1d, qs1d, ns1d, qg1d, ng1d,   &
 914                       t1d, p1d, dBZ, kts, kte, i, j)
 915           do k = kts, kte
 916              refl_10cm(i,k,j) = MAX(-35., dBZ(k))
 917           enddo
 918          endif
 919          ENDIF
 920 !+---+-----------------------------------------------------------------+
 921
 922    end do
 923    end do
 924
 925 END SUBROUTINE MP_MORR_TWO_MOMENT
 926
 927 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
 928 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
 929       SUBROUTINE MORR_TWO_MOMENT_MICRO(QC3DTEN,QI3DTEN,QNI3DTEN,QR3DTEN,         &
 930        NI3DTEN,NS3DTEN,NR3DTEN,QC3D,QI3D,QNI3D,QR3D,NI3D,NS3D,NR3D,              &
 931        T3DTEN,QV3DTEN,T3D,QV3D,PRES,DZQ,W3D,WVAR,PRECRT,SNOWRT,            &
 932        SNOWPRT,GRPLPRT,                & ! hm added 7/13/13
 933        EFFC,EFFI,EFFS,EFFR,DT,                                                   &
 934                                             IMS,IME, JMS,JME, KMS,KME,           &
 935                                             ITS,ITE, JTS,JTE, KTS,KTE,           & ! ADD GRAUPEL
 936                         QG3DTEN,NG3DTEN,QG3D,NG3D,EFFG,qrcu1d,qscu1d, qicu1d,    &
 937                         QGSTEN,QRSTEN,QISTEN,QNISTEN,QCSTEN, &
 938                         nc3d,nc3dten,iinum, & ! wrf-chem
 939                                 c2prec,CSED,ISED,SSED,GSED,RSED  &  ! hm added, wrf-chem
 940 #if (WRF_CHEM == 1)
 941         ,has_wetscav,rainprod, evapprod &
 942 #endif
 943                         )
 944
 945 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
 946 ! THIS PROGRAM IS THE MAIN TWO-MOMENT MICROPHYSICS SUBROUTINE DESCRIBED BY
 947 ! MORRISON ET AL. 2005 JAS AND MORRISON ET AL. 2009 MWR
 948
 949 ! THIS SCHEME IS A BULK DOUBLE-MOMENT SCHEME THAT PREDICTS MIXING
 950 ! RATIOS AND NUMBER CONCENTRATIONS OF FIVE HYDROMETEOR SPECIES:
 951 ! CLOUD DROPLETS, CLOUD (SMALL) ICE, RAIN, SNOW, AND GRAUPEL/HAIL.
 952
 953 ! CODE STRUCTURE: MAIN SUBROUTINE IS 'MORR_TWO_MOMENT'. ALSO INCLUDED IN THIS FILE IS
 954 ! 'FUNCTION POLYSVP', 'FUNCTION DERF1', AND
 955 ! 'FUNCTION GAMMA'.
 956
 957 ! NOTE: THIS SUBROUTINE USES 1D ARRAY IN VERTICAL (COLUMN), EVEN THOUGH VARIABLES ARE CALLED '3D'......
 958
 959 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
 960
 961 ! DECLARATIONS
 962
 963       IMPLICIT NONE
 964
 965 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
 966 ! THESE VARIABLES BELOW MUST BE LINKED WITH THE MAIN MODEL.
 967 ! DEFINE ARRAY SIZES
 968
 969 ! INPUT NUMBER OF GRID CELLS
 970
 971 ! INPUT/OUTPUT PARAMETERS                                 ! DESCRIPTION (UNITS)
 972       INTEGER, INTENT( IN)  :: IMS,IME, JMS,JME, KMS,KME,          &
 973                                ITS,ITE, JTS,JTE, KTS,KTE
 974
 975       REAL, DIMENSION(KTS:KTE) ::  QC3DTEN            ! CLOUD WATER MIXING RATIO TENDENCY (KG/KG/S)
 976       REAL, DIMENSION(KTS:KTE) ::  QI3DTEN            ! CLOUD ICE MIXING RATIO TENDENCY (KG/KG/S)
 977       REAL, DIMENSION(KTS:KTE) ::  QNI3DTEN           ! SNOW MIXING RATIO TENDENCY (KG/KG/S)
 978       REAL, DIMENSION(KTS:KTE) ::  QR3DTEN            ! RAIN MIXING RATIO TENDENCY (KG/KG/S)
 979       REAL, DIMENSION(KTS:KTE) ::  NI3DTEN            ! CLOUD ICE NUMBER CONCENTRATION (1/KG/S)
 980       REAL, DIMENSION(KTS:KTE) ::  NS3DTEN            ! SNOW NUMBER CONCENTRATION (1/KG/S)
 981       REAL, DIMENSION(KTS:KTE) ::  NR3DTEN            ! RAIN NUMBER CONCENTRATION (1/KG/S)
 982       REAL, DIMENSION(KTS:KTE) ::  QC3D               ! CLOUD WATER MIXING RATIO (KG/KG)
 983       REAL, DIMENSION(KTS:KTE) ::  QI3D               ! CLOUD ICE MIXING RATIO (KG/KG)
 984       REAL, DIMENSION(KTS:KTE) ::  QNI3D              ! SNOW MIXING RATIO (KG/KG)
 985       REAL, DIMENSION(KTS:KTE) ::  QR3D               ! RAIN MIXING RATIO (KG/KG)
 986       REAL, DIMENSION(KTS:KTE) ::  NI3D               ! CLOUD ICE NUMBER CONCENTRATION (1/KG)
 987       REAL, DIMENSION(KTS:KTE) ::  NS3D               ! SNOW NUMBER CONCENTRATION (1/KG)
 988       REAL, DIMENSION(KTS:KTE) ::  NR3D               ! RAIN NUMBER CONCENTRATION (1/KG)
 989       REAL, DIMENSION(KTS:KTE) ::  T3DTEN             ! TEMPERATURE TENDENCY (K/S)
 990       REAL, DIMENSION(KTS:KTE) ::  QV3DTEN            ! WATER VAPOR MIXING RATIO TENDENCY (KG/KG/S)
 991       REAL, DIMENSION(KTS:KTE) ::  T3D                ! TEMPERATURE (K)
 992       REAL, DIMENSION(KTS:KTE) ::  QV3D               ! WATER VAPOR MIXING RATIO (KG/KG)
 993       REAL, DIMENSION(KTS:KTE) ::  PRES               ! ATMOSPHERIC PRESSURE (PA)
 994       REAL, DIMENSION(KTS:KTE) ::  DZQ                ! DIFFERENCE IN HEIGHT ACROSS LEVEL (m)
 995       REAL, DIMENSION(KTS:KTE) ::  W3D                ! GRID-SCALE VERTICAL VELOCITY (M/S)
 996       REAL, DIMENSION(KTS:KTE) ::  WVAR               ! SUB-GRID VERTICAL VELOCITY (M/S)
 997 ! below for wrf-chem
 998       REAL, DIMENSION(KTS:KTE) ::  nc3d
 999       REAL, DIMENSION(KTS:KTE) ::  nc3dten
1000       integer, intent(in) :: iinum
1001
1002 ! HM ADDED GRAUPEL VARIABLES
1003       REAL, DIMENSION(KTS:KTE) ::  QG3DTEN            ! GRAUPEL MIX RATIO TENDENCY (KG/KG/S)
1004       REAL, DIMENSION(KTS:KTE) ::  NG3DTEN            ! GRAUPEL NUMB CONC TENDENCY (1/KG/S)
1005       REAL, DIMENSION(KTS:KTE) ::  QG3D            ! GRAUPEL MIX RATIO (KG/KG)
1006       REAL, DIMENSION(KTS:KTE) ::  NG3D            ! GRAUPEL NUMBER CONC (1/KG)
1007
1008 ! HM, ADD 1/16/07, SEDIMENTATION TENDENCIES FOR MIXING RATIO
1009
1010       REAL, DIMENSION(KTS:KTE) ::  QGSTEN            ! GRAUPEL SED TEND (KG/KG/S)
1011       REAL, DIMENSION(KTS:KTE) ::  QRSTEN            ! RAIN SED TEND (KG/KG/S)
1012       REAL, DIMENSION(KTS:KTE) ::  QISTEN            ! CLOUD ICE SED TEND (KG/KG/S)
1013       REAL, DIMENSION(KTS:KTE) ::  QNISTEN           ! SNOW SED TEND (KG/KG/S)
1014       REAL, DIMENSION(KTS:KTE) ::  QCSTEN            ! CLOUD WAT SED TEND (KG/KG/S)
1015
1016 ! hm add cumulus tendencies for precip
1017         REAL, DIMENSION(KTS:KTE) ::   qrcu1d
1018         REAL, DIMENSION(KTS:KTE) ::   qscu1d
1019         REAL, DIMENSION(KTS:KTE) ::   qicu1d
1020
1021 ! OUTPUT VARIABLES
1022
1023         REAL PRECRT                ! TOTAL PRECIP PER TIME STEP (mm)
1024         REAL SNOWRT                ! SNOW PER TIME STEP (mm)
1025 ! hm added 7/13/13
1026         REAL SNOWPRT      ! TOTAL CLOUD ICE PLUS SNOW PER TIME STEP (mm)
1027         REAL GRPLPRT      ! TOTAL GRAUPEL PER TIME STEP (mm)
1028
1029         REAL, DIMENSION(KTS:KTE) ::   EFFC            ! DROPLET EFFECTIVE RADIUS (MICRON)
1030         REAL, DIMENSION(KTS:KTE) ::   EFFI            ! CLOUD ICE EFFECTIVE RADIUS (MICRON)
1031         REAL, DIMENSION(KTS:KTE) ::   EFFS            ! SNOW EFFECTIVE RADIUS (MICRON)
1032         REAL, DIMENSION(KTS:KTE) ::   EFFR            ! RAIN EFFECTIVE RADIUS (MICRON)
1033         REAL, DIMENSION(KTS:KTE) ::   EFFG            ! GRAUPEL EFFECTIVE RADIUS (MICRON)
1034
1035 ! MODEL INPUT PARAMETERS (FORMERLY IN COMMON BLOCKS)
1036
1037         REAL DT         ! MODEL TIME STEP (SEC)
1038
1039 #if (WRF_CHEM == 1)
1040       LOGICAL, INTENT(IN) :: has_wetscav
1041 #endif
1042
1043
1044 !.....................................................................................................
1045 ! LOCAL VARIABLES: ALL PARAMETERS BELOW ARE LOCAL TO SCHEME AND DON'T NEED TO COMMUNICATE WITH THE
1046 ! REST OF THE MODEL.
1047
1048 ! SIZE PARAMETER VARIABLES
1049
1050      REAL, DIMENSION(KTS:KTE) :: LAMC          ! SLOPE PARAMETER FOR DROPLETS (M-1)
1051      REAL, DIMENSION(KTS:KTE) :: LAMI          ! SLOPE PARAMETER FOR CLOUD ICE (M-1)
1052      REAL, DIMENSION(KTS:KTE) :: LAMS          ! SLOPE PARAMETER FOR SNOW (M-1)
1053      REAL, DIMENSION(KTS:KTE) :: LAMR          ! SLOPE PARAMETER FOR RAIN (M-1)
1054      REAL, DIMENSION(KTS:KTE) :: LAMG          ! SLOPE PARAMETER FOR GRAUPEL (M-1)
1055      REAL, DIMENSION(KTS:KTE) :: CDIST1        ! PSD PARAMETER FOR DROPLETS
1056      REAL, DIMENSION(KTS:KTE) :: N0I           ! INTERCEPT PARAMETER FOR CLOUD ICE (KG-1 M-1)
1057      REAL, DIMENSION(KTS:KTE) :: N0S           ! INTERCEPT PARAMETER FOR SNOW (KG-1 M-1)
1058      REAL, DIMENSION(KTS:KTE) :: N0RR          ! INTERCEPT PARAMETER FOR RAIN (KG-1 M-1)
1059      REAL, DIMENSION(KTS:KTE) :: N0G           ! INTERCEPT PARAMETER FOR GRAUPEL (KG-1 M-1)
1060      REAL, DIMENSION(KTS:KTE) :: PGAM          ! SPECTRAL SHAPE PARAMETER FOR DROPLETS
1061
1062 ! MICROPHYSICAL PROCESSES
1063
1064      REAL, DIMENSION(KTS:KTE) ::  NSUBC     ! LOSS OF NC DURING EVAP
1065      REAL, DIMENSION(KTS:KTE) ::  NSUBI     ! LOSS OF NI DURING SUB.
1066      REAL, DIMENSION(KTS:KTE) ::  NSUBS     ! LOSS OF NS DURING SUB.
1067      REAL, DIMENSION(KTS:KTE) ::  NSUBR     ! LOSS OF NR DURING EVAP
1068      REAL, DIMENSION(KTS:KTE) ::  PRD       ! DEP CLOUD ICE
1069      REAL, DIMENSION(KTS:KTE) ::  PRE       ! EVAP OF RAIN
1070      REAL, DIMENSION(KTS:KTE) ::  PRDS      ! DEP SNOW
1071      REAL, DIMENSION(KTS:KTE) ::  NNUCCC    ! CHANGE N DUE TO CONTACT FREEZ DROPLETS
1072      REAL, DIMENSION(KTS:KTE) ::  MNUCCC    ! CHANGE Q DUE TO CONTACT FREEZ DROPLETS
1073      REAL, DIMENSION(KTS:KTE) ::  PRA       ! ACCRETION DROPLETS BY RAIN
1074      REAL, DIMENSION(KTS:KTE) ::  PRC       ! AUTOCONVERSION DROPLETS
1075      REAL, DIMENSION(KTS:KTE) ::  PCC       ! COND/EVAP DROPLETS
1076      REAL, DIMENSION(KTS:KTE) ::  NNUCCD    ! CHANGE N FREEZING AEROSOL (PRIM ICE NUCLEATION)
1077      REAL, DIMENSION(KTS:KTE) ::  MNUCCD    ! CHANGE Q FREEZING AEROSOL (PRIM ICE NUCLEATION)
1078      REAL, DIMENSION(KTS:KTE) ::  MNUCCR    ! CHANGE Q DUE TO CONTACT FREEZ RAIN
1079      REAL, DIMENSION(KTS:KTE) ::  NNUCCR    ! CHANGE N DUE TO CONTACT FREEZ RAIN
1080      REAL, DIMENSION(KTS:KTE) ::  NPRA      ! CHANGE IN N DUE TO DROPLET ACC BY RAIN
1081      REAL, DIMENSION(KTS:KTE) ::  NRAGG     ! SELF-COLLECTION/BREAKUP OF RAIN
1082      REAL, DIMENSION(KTS:KTE) ::  NSAGG     ! SELF-COLLECTION OF SNOW
1083      REAL, DIMENSION(KTS:KTE) ::  NPRC      ! CHANGE NC AUTOCONVERSION DROPLETS
1084      REAL, DIMENSION(KTS:KTE) ::  NPRC1      ! CHANGE NR AUTOCONVERSION DROPLETS
1085      REAL, DIMENSION(KTS:KTE) ::  PRAI      ! CHANGE Q ACCRETION CLOUD ICE BY SNOW
1086      REAL, DIMENSION(KTS:KTE) ::  PRCI      ! CHANGE Q AUTOCONVERSIN CLOUD ICE TO SNOW
1087      REAL, DIMENSION(KTS:KTE) ::  PSACWS    ! CHANGE Q DROPLET ACCRETION BY SNOW
1088      REAL, DIMENSION(KTS:KTE) ::  NPSACWS   ! CHANGE N DROPLET ACCRETION BY SNOW
1089      REAL, DIMENSION(KTS:KTE) ::  PSACWI    ! CHANGE Q DROPLET ACCRETION BY CLOUD ICE
1090      REAL, DIMENSION(KTS:KTE) ::  NPSACWI   ! CHANGE N DROPLET ACCRETION BY CLOUD ICE
1091      REAL, DIMENSION(KTS:KTE) ::  NPRCI     ! CHANGE N AUTOCONVERSION CLOUD ICE BY SNOW
1092      REAL, DIMENSION(KTS:KTE) ::  NPRAI     ! CHANGE N ACCRETION CLOUD ICE
1093      REAL, DIMENSION(KTS:KTE) ::  NMULTS    ! ICE MULT DUE TO RIMING DROPLETS BY SNOW
1094      REAL, DIMENSION(KTS:KTE) ::  NMULTR    ! ICE MULT DUE TO RIMING RAIN BY SNOW
1095      REAL, DIMENSION(KTS:KTE) ::  QMULTS    ! CHANGE Q DUE TO ICE MULT DROPLETS/SNOW
1096      REAL, DIMENSION(KTS:KTE) ::  QMULTR    ! CHANGE Q DUE TO ICE RAIN/SNOW
1097      REAL, DIMENSION(KTS:KTE) ::  PRACS     ! CHANGE Q RAIN-SNOW COLLECTION
1098      REAL, DIMENSION(KTS:KTE) ::  NPRACS    ! CHANGE N RAIN-SNOW COLLECTION
1099      REAL, DIMENSION(KTS:KTE) ::  PCCN      ! CHANGE Q DROPLET ACTIVATION
1100      REAL, DIMENSION(KTS:KTE) ::  PSMLT     ! CHANGE Q MELTING SNOW TO RAIN
1101      REAL, DIMENSION(KTS:KTE) ::  EVPMS     ! CHNAGE Q MELTING SNOW EVAPORATING
1102      REAL, DIMENSION(KTS:KTE) ::  NSMLTS    ! CHANGE N MELTING SNOW
1103      REAL, DIMENSION(KTS:KTE) ::  NSMLTR    ! CHANGE N MELTING SNOW TO RAIN
1104 ! HM ADDED 12/13/06
1105      REAL, DIMENSION(KTS:KTE) ::  PIACR     ! CHANGE QR, ICE-RAIN COLLECTION
1106      REAL, DIMENSION(KTS:KTE) ::  NIACR     ! CHANGE N, ICE-RAIN COLLECTION
1107      REAL, DIMENSION(KTS:KTE) ::  PRACI     ! CHANGE QI, ICE-RAIN COLLECTION
1108      REAL, DIMENSION(KTS:KTE) ::  PIACRS     ! CHANGE QR, ICE RAIN COLLISION, ADDED TO SNOW
1109      REAL, DIMENSION(KTS:KTE) ::  NIACRS     ! CHANGE N, ICE RAIN COLLISION, ADDED TO SNOW
1110      REAL, DIMENSION(KTS:KTE) ::  PRACIS     ! CHANGE QI, ICE RAIN COLLISION, ADDED TO SNOW
1111      REAL, DIMENSION(KTS:KTE) ::  EPRD      ! SUBLIMATION CLOUD ICE
1112      REAL, DIMENSION(KTS:KTE) ::  EPRDS     ! SUBLIMATION SNOW
1113 ! HM ADDED GRAUPEL PROCESSES
1114      REAL, DIMENSION(KTS:KTE) ::  PRACG    ! CHANGE IN Q COLLECTION RAIN BY GRAUPEL
1115      REAL, DIMENSION(KTS:KTE) ::  PSACWG    ! CHANGE IN Q COLLECTION DROPLETS BY GRAUPEL
1116      REAL, DIMENSION(KTS:KTE) ::  PGSACW    ! CONVERSION Q TO GRAUPEL DUE TO COLLECTION DROPLETS BY SNOW
1117      REAL, DIMENSION(KTS:KTE) ::  PGRACS    ! CONVERSION Q TO GRAUPEL DUE TO COLLECTION RAIN BY SNOW
1118      REAL, DIMENSION(KTS:KTE) ::  PRDG    ! DEP OF GRAUPEL
1119      REAL, DIMENSION(KTS:KTE) ::  EPRDG    ! SUB OF GRAUPEL
1120      REAL, DIMENSION(KTS:KTE) ::  EVPMG    ! CHANGE Q MELTING OF GRAUPEL AND EVAPORATION
1121      REAL, DIMENSION(KTS:KTE) ::  PGMLT    ! CHANGE Q MELTING OF GRAUPEL
1122      REAL, DIMENSION(KTS:KTE) ::  NPRACG    ! CHANGE N COLLECTION RAIN BY GRAUPEL
1123      REAL, DIMENSION(KTS:KTE) ::  NPSACWG    ! CHANGE N COLLECTION DROPLETS BY GRAUPEL
1124      REAL, DIMENSION(KTS:KTE) ::  NSCNG    ! CHANGE N CONVERSION TO GRAUPEL DUE TO COLLECTION DROPLETS BY SNOW
1125      REAL, DIMENSION(KTS:KTE) ::  NGRACS    ! CHANGE N CONVERSION TO GRAUPEL DUE TO COLLECTION RAIN BY SNOW
1126      REAL, DIMENSION(KTS:KTE) ::  NGMLTG    ! CHANGE N MELTING GRAUPEL
1127      REAL, DIMENSION(KTS:KTE) ::  NGMLTR    ! CHANGE N MELTING GRAUPEL TO RAIN
1128      REAL, DIMENSION(KTS:KTE) ::  NSUBG    ! CHANGE N SUB/DEP OF GRAUPEL
1129      REAL, DIMENSION(KTS:KTE) ::  PSACR    ! CONVERSION DUE TO COLL OF SNOW BY RAIN
1130      REAL, DIMENSION(KTS:KTE) ::  NMULTG    ! ICE MULT DUE TO ACC DROPLETS BY GRAUPEL
1131      REAL, DIMENSION(KTS:KTE) ::  NMULTRG    ! ICE MULT DUE TO ACC RAIN BY GRAUPEL
1132      REAL, DIMENSION(KTS:KTE) ::  QMULTG    ! CHANGE Q DUE TO ICE MULT DROPLETS/GRAUPEL
1133      REAL, DIMENSION(KTS:KTE) ::  QMULTRG    ! CHANGE Q DUE TO ICE MULT RAIN/GRAUPEL
1134
1135 ! TIME-VARYING ATMOSPHERIC PARAMETERS
1136
1137      REAL, DIMENSION(KTS:KTE) ::   KAP   ! THERMAL CONDUCTIVITY OF AIR
1138      REAL, DIMENSION(KTS:KTE) ::   EVS   ! SATURATION VAPOR PRESSURE
1139      REAL, DIMENSION(KTS:KTE) ::   EIS   ! ICE SATURATION VAPOR PRESSURE
1140      REAL, DIMENSION(KTS:KTE) ::   QVS   ! SATURATION MIXING RATIO
1141      REAL, DIMENSION(KTS:KTE) ::   QVI   ! ICE SATURATION MIXING RATIO
1142      REAL, DIMENSION(KTS:KTE) ::   QVQVS ! SAUTRATION RATIO
1143      REAL, DIMENSION(KTS:KTE) ::   QVQVSI! ICE SATURAION RATIO
1144      REAL, DIMENSION(KTS:KTE) ::   DV    ! DIFFUSIVITY OF WATER VAPOR IN AIR
1145      REAL, DIMENSION(KTS:KTE) ::   XXLS  ! LATENT HEAT OF SUBLIMATION
1146      REAL, DIMENSION(KTS:KTE) ::   XXLV  ! LATENT HEAT OF VAPORIZATION
1147      REAL, DIMENSION(KTS:KTE) ::   CPM   ! SPECIFIC HEAT AT CONST PRESSURE FOR MOIST AIR
1148      REAL, DIMENSION(KTS:KTE) ::   MU    ! VISCOCITY OF AIR
1149      REAL, DIMENSION(KTS:KTE) ::   SC    ! SCHMIDT NUMBER
1150      REAL, DIMENSION(KTS:KTE) ::   XLF   ! LATENT HEAT OF FREEZING
1151      REAL, DIMENSION(KTS:KTE) ::   RHO   ! AIR DENSITY
1152      REAL, DIMENSION(KTS:KTE) ::   AB    ! CORRECTION TO CONDENSATION RATE DUE TO LATENT HEATING
1153      REAL, DIMENSION(KTS:KTE) ::   ABI    ! CORRECTION TO DEPOSITION RATE DUE TO LATENT HEATING
1154
1155 ! TIME-VARYING MICROPHYSICS PARAMETERS
1156
1157      REAL, DIMENSION(KTS:KTE) ::   DAP    ! DIFFUSIVITY OF AEROSOL
1158      REAL    NACNT                    ! NUMBER OF CONTACT IN
1159      REAL    FMULT                    ! TEMP.-DEP. PARAMETER FOR RIME-SPLINTERING
1160      REAL    COFFI                    ! ICE AUTOCONVERSION PARAMETER
1161
1162 ! FALL SPEED WORKING VARIABLES (DEFINED IN CODE)
1163
1164       REAL, DIMENSION(KTS:KTE) ::    DUMI,DUMR,DUMFNI,DUMG,DUMFNG
1165       REAL UNI, UMI,UMR
1166       REAL, DIMENSION(KTS:KTE) ::    FR, FI, FNI,FG,FNG
1167       REAL RGVM
1168       REAL, DIMENSION(KTS:KTE) ::   FALOUTR,FALOUTI,FALOUTNI
1169       REAL FALTNDR,FALTNDI,FALTNDNI,RHO2
1170       REAL, DIMENSION(KTS:KTE) ::   DUMQS,DUMFNS
1171       REAL UMS,UNS
1172       REAL, DIMENSION(KTS:KTE) ::   FS,FNS, FALOUTS,FALOUTNS,FALOUTG,FALOUTNG
1173       REAL FALTNDS,FALTNDNS,UNR,FALTNDG,FALTNDNG
1174       REAL, DIMENSION(KTS:KTE) ::    DUMC,DUMFNC
1175       REAL UNC,UMC,UNG,UMG
1176       REAL, DIMENSION(KTS:KTE) ::   FC,FALOUTC,FALOUTNC
1177       REAL FALTNDC,FALTNDNC
1178       REAL, DIMENSION(KTS:KTE) ::   FNC,DUMFNR,FALOUTNR
1179       REAL FALTNDNR
1180       REAL, DIMENSION(KTS:KTE) ::   FNR
1181
1182 ! FALL-SPEED PARAMETER 'A' WITH AIR DENSITY CORRECTION
1183
1184       REAL, DIMENSION(KTS:KTE) ::    AIN,ARN,ASN,ACN,AGN
1185
1186 ! EXTERNAL FUNCTION CALL RETURN VARIABLES
1187
1188 !      REAL GAMMA,      ! EULER GAMMA FUNCTION
1189 !      REAL POLYSVP,    ! SAT. PRESSURE FUNCTION
1190 !      REAL DERF1        ! ERROR FUNCTION
1191
1192 ! DUMMY VARIABLES
1193
1194      REAL DUM,DUM1,DUM2,DUMT,DUMQV,DUMQSS,DUMQSI,DUMS
1195
1196 ! PROGNOSTIC SUPERSATURATION
1197
1198      REAL DQSDT    ! CHANGE OF SAT. MIX. RAT. WITH TEMPERATURE
1199      REAL DQSIDT   ! CHANGE IN ICE SAT. MIXING RAT. WITH T
1200      REAL EPSI     ! 1/PHASE REL. TIME (SEE M2005), ICE
1201      REAL EPSS     ! 1/PHASE REL. TIME (SEE M2005), SNOW
1202      REAL EPSR     ! 1/PHASE REL. TIME (SEE M2005), RAIN
1203      REAL EPSG     ! 1/PHASE REL. TIME (SEE M2005), GRAUPEL
1204
1205 ! NEW DROPLET ACTIVATION VARIABLES
1206      REAL TAUC     ! PHASE REL. TIME (SEE M2005), DROPLETS
1207      REAL TAUR     ! PHASE REL. TIME (SEE M2005), RAIN
1208      REAL TAUI     ! PHASE REL. TIME (SEE M2005), CLOUD ICE
1209      REAL TAUS     ! PHASE REL. TIME (SEE M2005), SNOW
1210      REAL TAUG     ! PHASE REL. TIME (SEE M2005), GRAUPEL
1211      REAL DUMACT,DUM3
1212
1213 ! COUNTING/INDEX VARIABLES
1214
1215      INTEGER K,NSTEP,N ! ,I
1216
1217 ! LTRUE IS ONLY USED TO SPEED UP THE CODE !!
1218 ! LTRUE, SWITCH = 0, NO HYDROMETEORS IN COLUMN,
1219 !               = 1, HYDROMETEORS IN COLUMN
1220
1221       INTEGER LTRUE
1222
1223 ! DROPLET ACTIVATION/FREEZING AEROSOL
1224
1225
1226      REAL    CT      ! DROPLET ACTIVATION PARAMETER
1227      REAL    TEMP1   ! DUMMY TEMPERATURE
1228      REAL    SAT1    ! DUMMY SATURATION
1229      REAL    SIGVL   ! SURFACE TENSION LIQ/VAPOR
1230      REAL    KEL     ! KELVIN PARAMETER
1231      REAL    KC2     ! TOTAL ICE NUCLEATION RATE
1232
1233        REAL CRY,KRY   ! AEROSOL ACTIVATION PARAMETERS
1234
1235 ! MORE WORKING/DUMMY VARIABLES
1236
1237      REAL DUMQI,DUMNI,DC0,DS0,DG0
1238      REAL DUMQC,DUMQR,RATIO,SUM_DEP,FUDGEF
1239
1240 ! EFFECTIVE VERTICAL VELOCITY  (M/S)
1241      REAL WEF
1242
1243 ! WORKING PARAMETERS FOR ICE NUCLEATION
1244
1245       REAL ANUC,BNUC
1246
1247 ! WORKING PARAMETERS FOR AEROSOL ACTIVATION
1248
1249         REAL AACT,GAMM,GG,PSI,ETA1,ETA2,SM1,SM2,SMAX,UU1,UU2,ALPHA
1250
1251 ! DUMMY SIZE DISTRIBUTION PARAMETERS
1252
1253         REAL DLAMS,DLAMR,DLAMI,DLAMC,DLAMG,LAMMAX,LAMMIN
1254
1255         INTEGER IDROP
1256
1257 ! FOR WRF-CHEM
1258         REAL, DIMENSION(KTS:KTE)::C2PREC,CSED,ISED,SSED,GSED,RSED
1259 #if (WRF_CHEM == 1)
1260     REAL, DIMENSION(KTS:KTE), INTENT(INOUT) :: rainprod, evapprod
1261 #endif
1262     REAL, DIMENSION(KTS:KTE)                :: tqimelt ! melting of cloud ice (tendency)
1263
1264 ! comment lines for wrf-chem since these are intent(in) in that case
1265 !       REAL, DIMENSION(KTS:KTE) ::  NC3DTEN            ! CLOUD DROPLET NUMBER CONCENTRATION (1/KG/S)
1266 !       REAL, DIMENSION(KTS:KTE) ::  NC3D               ! CLOUD DROPLET NUMBER CONCENTRATION (1/KG)
1267
1268 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1269
1270 ! SET LTRUE INITIALLY TO 0
1271
1272          LTRUE = 0
1273
1274 ! ATMOSPHERIC PARAMETERS THAT VARY IN TIME AND HEIGHT
1275          DO K = KTS,KTE
1276
1277 ! NC3DTEN LOCAL ARRAY INITIALIZED
1278                NC3DTEN(K) = 0.
1279 ! INITIALIZE VARIABLES FOR WRF-CHEM OUTPUT TO ZERO
1280
1281                 C2PREC(K)=0.
1282                 CSED(K)=0.
1283                 ISED(K)=0.
1284                 SSED(K)=0.
1285                 GSED(K)=0.
1286                 RSED(K)=0.
1287
1288 #if (WRF_CHEM == 1)
1289          rainprod(K) = 0.
1290          evapprod(K) = 0.
1291          tqimelt(K)  = 0.
1292          PRC(K)      = 0.
1293          PRA(K)      = 0.
1294 #endif
1295
1296 ! LATENT HEAT OF VAPORATION
1297
1298             XXLV(K) = 3.1484E6-2370.*T3D(K)
1299
1300 ! LATENT HEAT OF SUBLIMATION
1301
1302             XXLS(K) = 3.15E6-2370.*T3D(K)+0.3337E6
1303
1304             CPM(K) = CP*(1.+0.887*QV3D(K))
1305
1306 ! SATURATION VAPOR PRESSURE AND MIXING RATIO
1307
1308 ! hm, add fix for low pressure, 5/12/10
1309
1310             EVS(K) = min(0.99*pres(k),POLYSVP(T3D(K),0))   ! PA
1311             EIS(K) = min(0.99*pres(k),POLYSVP(T3D(K),1))   ! PA
1312
1313 ! MAKE SURE ICE SATURATION DOESN'T EXCEED WATER SAT. NEAR FREEZING
1314
1315             IF (EIS(K).GT.EVS(K)) EIS(K) = EVS(K)
1316
1317             QVS(K) = EP_2*EVS(K)/(PRES(K)-EVS(K))
1318             QVI(K) = EP_2*EIS(K)/(PRES(K)-EIS(K))
1319
1320             QVQVS(K) = QV3D(K)/QVS(K)
1321             QVQVSI(K) = QV3D(K)/QVI(K)
1322
1323 ! AIR DENSITY
1324
1325             RHO(K) = PRES(K)/(R*T3D(K))
1326
1327 ! ADD NUMBER CONCENTRATION DUE TO CUMULUS TENDENCY
1328 ! ASSUME N0 ASSOCIATED WITH CUMULUS PARAM RAIN IS 10^7 M^-4
1329 ! ASSUME N0 ASSOCIATED WITH CUMULUS PARAM SNOW IS 2 X 10^7 M^-4
1330 ! FOR DETRAINED CLOUD ICE, ASSUME MEAN VOLUME DIAM OF 80 MICRON
1331
1332             IF (QRCU1D(K).GE.1.E-10) THEN
1333             DUM=1.8e5*(QRCU1D(K)*DT/(PI*RHOW*RHO(K)**3))**0.25
1334             NR3D(K)=NR3D(K)+DUM
1335             END IF
1336             IF (QSCU1D(K).GE.1.E-10) THEN
1337             DUM=3.e5*(QSCU1D(K)*DT/(CONS1*RHO(K)**3))**(1./(DS+1.))
1338             NS3D(K)=NS3D(K)+DUM
1339             END IF
1340             IF (QICU1D(K).GE.1.E-10) THEN
1341             DUM=QICU1D(K)*DT/(CI*(80.E-6)**DI)
1342             NI3D(K)=NI3D(K)+DUM
1343             END IF
1344
1345 ! AT SUBSATURATION, REMOVE SMALL AMOUNTS OF CLOUD/PRECIP WATER
1346 ! hm modify 7/0/09 change limit to 1.e-8
1347
1348              IF (QVQVS(K).LT.0.9) THEN
1349                IF (QR3D(K).LT.1.E-8) THEN
1350                   QV3D(K)=QV3D(K)+QR3D(K)
1351                   T3D(K)=T3D(K)-QR3D(K)*XXLV(K)/CPM(K)
1352                   QR3D(K)=0.
1353                END IF
1354                IF (QC3D(K).LT.1.E-8) THEN
1355                   QV3D(K)=QV3D(K)+QC3D(K)
1356                   T3D(K)=T3D(K)-QC3D(K)*XXLV(K)/CPM(K)
1357                   QC3D(K)=0.
1358                END IF
1359              END IF
1360
1361              IF (QVQVSI(K).LT.0.9) THEN
1362                IF (QI3D(K).LT.1.E-8) THEN
1363                   QV3D(K)=QV3D(K)+QI3D(K)
1364                   T3D(K)=T3D(K)-QI3D(K)*XXLS(K)/CPM(K)
1365                   QI3D(K)=0.
1366                END IF
1367                IF (QNI3D(K).LT.1.E-8) THEN
1368                   QV3D(K)=QV3D(K)+QNI3D(K)
1369                   T3D(K)=T3D(K)-QNI3D(K)*XXLS(K)/CPM(K)
1370                   QNI3D(K)=0.
1371                END IF
1372                IF (QG3D(K).LT.1.E-8) THEN
1373                   QV3D(K)=QV3D(K)+QG3D(K)
1374                   T3D(K)=T3D(K)-QG3D(K)*XXLS(K)/CPM(K)
1375                   QG3D(K)=0.
1376                END IF
1377              END IF
1378
1379 ! HEAT OF FUSION
1380
1381             XLF(K) = XXLS(K)-XXLV(K)
1382
1383 !..................................................................
1384 ! IF MIXING RATIO < QSMALL SET MIXING RATIO AND NUMBER CONC TO ZERO
1385
1386        IF (QC3D(K).LT.QSMALL) THEN
1387          QC3D(K) = 0.
1388          NC3D(K) = 0.
1389          EFFC(K) = 0.
1390        END IF
1391        IF (QR3D(K).LT.QSMALL) THEN
1392          QR3D(K) = 0.
1393          NR3D(K) = 0.
1394          EFFR(K) = 0.
1395        END IF
1396        IF (QI3D(K).LT.QSMALL) THEN
1397          QI3D(K) = 0.
1398          NI3D(K) = 0.
1399          EFFI(K) = 0.
1400        END IF
1401        IF (QNI3D(K).LT.QSMALL) THEN
1402          QNI3D(K) = 0.
1403          NS3D(K) = 0.
1404          EFFS(K) = 0.
1405        END IF
1406        IF (QG3D(K).LT.QSMALL) THEN
1407          QG3D(K) = 0.
1408          NG3D(K) = 0.
1409          EFFG(K) = 0.
1410        END IF
1411
1412 ! INITIALIZE SEDIMENTATION TENDENCIES FOR MIXING RATIO
1413
1414       QRSTEN(K) = 0.
1415       QISTEN(K) = 0.
1416       QNISTEN(K) = 0.
1417       QCSTEN(K) = 0.
1418       QGSTEN(K) = 0.
1419
1420 !..................................................................
1421 ! MICROPHYSICS PARAMETERS VARYING IN TIME/HEIGHT
1422
1423 ! fix 053011
1424             MU(K) = 1.496E-6*T3D(K)**1.5/(T3D(K)+120.)
1425
1426 ! FALL SPEED WITH DENSITY CORRECTION (HEYMSFIELD AND BENSSEMER 2006)
1427
1428             DUM = (RHOSU/RHO(K))**0.54
1429
1430 ! fix 053011
1431 !            AIN(K) = DUM*AI
1432 ! AA revision 4/1/11: Ikawa and Saito 1991 air-density correction
1433             AIN(K) = (RHOSU/RHO(K))**0.35*AI
1434             ARN(K) = DUM*AR
1435             ASN(K) = DUM*AS
1436 !            ACN(K) = DUM*AC
1437 ! AA revision 4/1/11: temperature-dependent Stokes fall speed
1438             ACN(K) = G*RHOW/(18.*MU(K))
1439 ! HM ADD GRAUPEL 8/28/06
1440             AGN(K) = DUM*AG
1441
1442 !hm 4/7/09 bug fix, initialize lami to prevent later division by zero
1443             LAMI(K)=0.
1444
1445 !..................................
1446 ! IF THERE IS NO CLOUD/PRECIP WATER, AND IF SUBSATURATED, THEN SKIP MICROPHYSICS
1447 ! FOR THIS LEVEL
1448
1449             IF (QC3D(K).LT.QSMALL.AND.QI3D(K).LT.QSMALL.AND.QNI3D(K).LT.QSMALL &
1450                  .AND.QR3D(K).LT.QSMALL.AND.QG3D(K).LT.QSMALL) THEN
1451                  IF (T3D(K).LT.273.15.AND.QVQVSI(K).LT.0.999) GOTO 200
1452                  IF (T3D(K).GE.273.15.AND.QVQVS(K).LT.0.999) GOTO 200
1453             END IF
1454
1455 ! THERMAL CONDUCTIVITY FOR AIR
1456
1457 ! fix 053011
1458             KAP(K) = 1.414E3*MU(K)
1459
1460 ! DIFFUSIVITY OF WATER VAPOR
1461
1462             DV(K) = 8.794E-5*T3D(K)**1.81/PRES(K)
1463
1464 ! SCHMIT NUMBER
1465
1466 ! fix 053011
1467             SC(K) = MU(K)/(RHO(K)*DV(K))
1468
1469 ! PSYCHOMETIC CORRECTIONS
1470
1471 ! RATE OF CHANGE SAT. MIX. RATIO WITH TEMPERATURE
1472
1473             DUM = (RV*T3D(K)**2)
1474
1475             DQSDT = XXLV(K)*QVS(K)/DUM
1476             DQSIDT =  XXLS(K)*QVI(K)/DUM
1477
1478             ABI(K) = 1.+DQSIDT*XXLS(K)/CPM(K)
1479             AB(K) = 1.+DQSDT*XXLV(K)/CPM(K)
1480
1481 !
1482 !.....................................................................
1483 !.....................................................................
1484 ! CASE FOR TEMPERATURE ABOVE FREEZING
1485
1486             IF (T3D(K).GE.273.15) THEN
1487
1488 !......................................................................
1489 !HM ADD, ALLOW FOR CONSTANT DROPLET NUMBER
1490 ! INUM = 0, PREDICT DROPLET NUMBER
1491 ! INUM = 1, SET CONSTANT DROPLET NUMBER
1492
1493          IF (iinum.EQ.1) THEN
1494 ! CONVERT NDCNST FROM CM-3 TO KG-1
1495             NC3D(K)=NDCNST*1.E6/RHO(K)
1496          END IF
1497
1498 ! GET SIZE DISTRIBUTION PARAMETERS
1499
1500 ! MELT VERY SMALL SNOW AND GRAUPEL MIXING RATIOS, ADD TO RAIN
1501        IF (QNI3D(K).LT.1.E-6) THEN
1502           QR3D(K)=QR3D(K)+QNI3D(K)
1503           NR3D(K)=NR3D(K)+NS3D(K)
1504           T3D(K)=T3D(K)-QNI3D(K)*XLF(K)/CPM(K)
1505           QNI3D(K) = 0.
1506           NS3D(K) = 0.
1507        END IF
1508        IF (QG3D(K).LT.1.E-6) THEN
1509           QR3D(K)=QR3D(K)+QG3D(K)
1510           NR3D(K)=NR3D(K)+NG3D(K)
1511           T3D(K)=T3D(K)-QG3D(K)*XLF(K)/CPM(K)
1512           QG3D(K) = 0.
1513           NG3D(K) = 0.
1514        END IF
1515
1516        IF (QC3D(K).LT.QSMALL.AND.QNI3D(K).LT.1.E-8.AND.QR3D(K).LT.QSMALL.AND.QG3D(K).LT.1.E-8) GOTO 300
1517
1518 ! MAKE SURE NUMBER CONCENTRATIONS AREN'T NEGATIVE
1519
1520       NS3D(K) = MAX(0.,NS3D(K))
1521       NC3D(K) = MAX(0.,NC3D(K))
1522       NR3D(K) = MAX(0.,NR3D(K))
1523       NG3D(K) = MAX(0.,NG3D(K))
1524
1525 !......................................................................
1526 ! RAIN
1527
1528       IF (QR3D(K).GE.QSMALL) THEN
1529       LAMR(K) = (PI*RHOW*NR3D(K)/QR3D(K))**(1./3.)
1530       N0RR(K) = NR3D(K)*LAMR(K)
1531
1532 ! CHECK FOR SLOPE
1533
1534 ! ADJUST VARS
1535
1536       IF (LAMR(K).LT.LAMMINR) THEN
1537
1538       LAMR(K) = LAMMINR
1539
1540       N0RR(K) = LAMR(K)**4*QR3D(K)/(PI*RHOW)
1541
1542       NR3D(K) = N0RR(K)/LAMR(K)
1543       ELSE IF (LAMR(K).GT.LAMMAXR) THEN
1544       LAMR(K) = LAMMAXR
1545       N0RR(K) = LAMR(K)**4*QR3D(K)/(PI*RHOW)
1546
1547       NR3D(K) = N0RR(K)/LAMR(K)
1548       END IF
1549       END IF
1550
1551 !......................................................................
1552 ! CLOUD DROPLETS
1553
1554 ! MARTIN ET AL. (1994) FORMULA FOR PGAM
1555
1556       IF (QC3D(K).GE.QSMALL) THEN
1557
1558          DUM = PRES(K)/(287.15*T3D(K))
1559          PGAM(K)=0.0005714*(NC3D(K)/1.E6*DUM)+0.2714
1560          PGAM(K)=1./(PGAM(K)**2)-1.
1561          PGAM(K)=MAX(PGAM(K),2.)
1562          PGAM(K)=MIN(PGAM(K),10.)
1563
1564 ! CALCULATE LAMC
1565
1566       LAMC(K) = (CONS26*NC3D(K)*GAMMA(PGAM(K)+4.)/   &
1567                  (QC3D(K)*GAMMA(PGAM(K)+1.)))**(1./3.)
1568
1569 ! LAMMIN, 60 MICRON DIAMETER
1570 ! LAMMAX, 1 MICRON
1571
1572       LAMMIN = (PGAM(K)+1.)/60.E-6
1573       LAMMAX = (PGAM(K)+1.)/1.E-6
1574
1575       IF (LAMC(K).LT.LAMMIN) THEN
1576       LAMC(K) = LAMMIN
1577
1578       NC3D(K) = EXP(3.*LOG(LAMC(K))+LOG(QC3D(K))+              &
1579                 LOG(GAMMA(PGAM(K)+1.))-LOG(GAMMA(PGAM(K)+4.)))/CONS26
1580       ELSE IF (LAMC(K).GT.LAMMAX) THEN
1581       LAMC(K) = LAMMAX
1582
1583       NC3D(K) = EXP(3.*LOG(LAMC(K))+LOG(QC3D(K))+              &
1584                 LOG(GAMMA(PGAM(K)+1.))-LOG(GAMMA(PGAM(K)+4.)))/CONS26
1585
1586       END IF
1587
1588       END IF
1589
1590 !......................................................................
1591 ! SNOW
1592
1593       IF (QNI3D(K).GE.QSMALL) THEN
1594       LAMS(K) = (CONS1*NS3D(K)/QNI3D(K))**(1./DS)
1595       N0S(K) = NS3D(K)*LAMS(K)
1596
1597 ! CHECK FOR SLOPE
1598
1599 ! ADJUST VARS
1600
1601       IF (LAMS(K).LT.LAMMINS) THEN
1602       LAMS(K) = LAMMINS
1603       N0S(K) = LAMS(K)**4*QNI3D(K)/CONS1
1604
1605       NS3D(K) = N0S(K)/LAMS(K)
1606
1607       ELSE IF (LAMS(K).GT.LAMMAXS) THEN
1608
1609       LAMS(K) = LAMMAXS
1610       N0S(K) = LAMS(K)**4*QNI3D(K)/CONS1
1611
1612       NS3D(K) = N0S(K)/LAMS(K)
1613       END IF
1614       END IF
1615
1616 !......................................................................
1617 ! GRAUPEL
1618
1619       IF (QG3D(K).GE.QSMALL) THEN
1620       LAMG(K) = (CONS2*NG3D(K)/QG3D(K))**(1./DG)
1621       N0G(K) = NG3D(K)*LAMG(K)
1622
1623 ! ADJUST VARS
1624
1625       IF (LAMG(K).LT.LAMMING) THEN
1626       LAMG(K) = LAMMING
1627       N0G(K) = LAMG(K)**4*QG3D(K)/CONS2
1628
1629       NG3D(K) = N0G(K)/LAMG(K)
1630
1631       ELSE IF (LAMG(K).GT.LAMMAXG) THEN
1632
1633       LAMG(K) = LAMMAXG
1634       N0G(K) = LAMG(K)**4*QG3D(K)/CONS2
1635
1636       NG3D(K) = N0G(K)/LAMG(K)
1637       END IF
1638       END IF
1639
1640 !.....................................................................
1641 ! ZERO OUT PROCESS RATES
1642
1643             PRC(K) = 0.
1644             NPRC(K) = 0.
1645             NPRC1(K) = 0.
1646             PRA(K) = 0.
1647             NPRA(K) = 0.
1648             NRAGG(K) = 0.
1649             NSMLTS(K) = 0.
1650             NSMLTR(K) = 0.
1651             EVPMS(K) = 0.
1652             PCC(K) = 0.
1653             PRE(K) = 0.
1654             NSUBC(K) = 0.
1655             NSUBR(K) = 0.
1656             PRACG(K) = 0.
1657             NPRACG(K) = 0.
1658             PSMLT(K) = 0.
1659             PGMLT(K) = 0.
1660             EVPMG(K) = 0.
1661             PRACS(K) = 0.
1662             NPRACS(K) = 0.
1663             NGMLTG(K) = 0.
1664             NGMLTR(K) = 0.
1665
1666 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1667 ! CALCULATION OF MICROPHYSICAL PROCESS RATES, T > 273.15 K
1668
1669 !.................................................................
1670 !.......................................................................
1671 ! AUTOCONVERSION OF CLOUD LIQUID WATER TO RAIN
1672 ! FORMULA FROM BEHENG (1994)
1673 ! USING NUMERICAL SIMULATION OF STOCHASTIC COLLECTION EQUATION
1674 ! AND INITIAL CLOUD DROPLET SIZE DISTRIBUTION SPECIFIED
1675 ! AS A GAMMA DISTRIBUTION
1676
1677 ! USE MINIMUM VALUE OF 1.E-6 TO PREVENT FLOATING POINT ERROR
1678
1679          IF (QC3D(K).GE.1.E-6) THEN
1680
1681 ! HM ADD 12/13/06, REPLACE WITH NEWER FORMULA
1682 ! FROM KHAIROUTDINOV AND KOGAN 2000, MWR
1683
1684                 PRC(K)=1350.*QC3D(K)**2.47*  &
1685            (NC3D(K)/1.e6*RHO(K))**(-1.79)
1686
1687 ! note: nprc1 is change in Nr,
1688 ! nprc is change in Nc
1689
1690         NPRC1(K) = PRC(K)/CONS29
1691         NPRC(K) = PRC(K)/(QC3D(k)/NC3D(K))
1692
1693 ! hm bug fix 3/20/12
1694                 NPRC(K) = MIN(NPRC(K),NC3D(K)/DT)
1695                 NPRC1(K) = MIN(NPRC1(K),NPRC(K))
1696
1697          END IF
1698
1699 !.......................................................................
1700 ! HM ADD 12/13/06, COLLECTION OF SNOW BY RAIN ABOVE FREEZING
1701 ! FORMULA FROM IKAWA AND SAITO (1991)
1702
1703          IF (QR3D(K).GE.1.E-8.AND.QNI3D(K).GE.1.E-8) THEN
1704
1705             UMS = ASN(K)*CONS3/(LAMS(K)**BS)
1706             UMR = ARN(K)*CONS4/(LAMR(K)**BR)
1707             UNS = ASN(K)*CONS5/LAMS(K)**BS
1708             UNR = ARN(K)*CONS6/LAMR(K)**BR
1709
1710 ! SET REASLISTIC LIMITS ON FALLSPEEDS
1711
1712 ! bug fix, 10/08/09
1713             dum=(rhosu/rho(k))**0.54
1714             UMS=MIN(UMS,1.2*dum)
1715             UNS=MIN(UNS,1.2*dum)
1716             UMR=MIN(UMR,9.1*dum)
1717             UNR=MIN(UNR,9.1*dum)
1718
1719 ! hm fix, 2/12/13
1720 ! for above freezing conditions to get accelerated melting of snow,
1721 ! we need collection of rain by snow (following Lin et al. 1983)
1722 !            PRACS(K) = CONS31*(((1.2*UMR-0.95*UMS)**2+              &
1723 !                  0.08*UMS*UMR)**0.5*RHO(K)*                     &
1724 !                 N0RR(K)*N0S(K)/LAMS(K)**3*                    &
1725 !                  (5./(LAMS(K)**3*LAMR(K))+                    &
1726 !                  2./(LAMS(K)**2*LAMR(K)**2)+                  &
1727 !                  0.5/(LAMS(K)*LAMR(K)**3)))
1728
1729             PRACS(K) = CONS41*(((1.2*UMR-0.95*UMS)**2+                   &
1730                   0.08*UMS*UMR)**0.5*RHO(K)*                      &
1731                   N0RR(K)*N0S(K)/LAMR(K)**3*                              &
1732                   (5./(LAMR(K)**3*LAMS(K))+                    &
1733                   2./(LAMR(K)**2*LAMS(K)**2)+                  &
1734                   0.5/(LAMR(k)*LAMS(k)**3)))
1735
1736 ! fix 053011, npracs no longer subtracted from snow
1737 !            NPRACS(K) = CONS32*RHO(K)*(1.7*(UNR-UNS)**2+            &
1738 !                0.3*UNR*UNS)**0.5*N0RR(K)*N0S(K)*              &
1739 !                (1./(LAMR(K)**3*LAMS(K))+                      &
1740 !                 1./(LAMR(K)**2*LAMS(K)**2)+                   &
1741 !                 1./(LAMR(K)*LAMS(K)**3))
1742
1743          END IF
1744
1745 ! ADD COLLECTION OF GRAUPEL BY RAIN ABOVE FREEZING
1746 ! ASSUME ALL RAIN COLLECTION BY GRAUPEL ABOVE FREEZING IS SHED
1747 ! ASSUME SHED DROPS ARE 1 MM IN SIZE
1748
1749          IF (QR3D(K).GE.1.E-8.AND.QG3D(K).GE.1.E-8) THEN
1750
1751             UMG = AGN(K)*CONS7/(LAMG(K)**BG)
1752             UMR = ARN(K)*CONS4/(LAMR(K)**BR)
1753             UNG = AGN(K)*CONS8/LAMG(K)**BG
1754             UNR = ARN(K)*CONS6/LAMR(K)**BR
1755
1756 ! SET REASLISTIC LIMITS ON FALLSPEEDS
1757 ! bug fix, 10/08/09
1758             dum=(rhosu/rho(k))**0.54
1759             UMG=MIN(UMG,20.*dum)
1760             UNG=MIN(UNG,20.*dum)
1761             UMR=MIN(UMR,9.1*dum)
1762             UNR=MIN(UNR,9.1*dum)
1763
1764 ! PRACG IS MIXING RATIO OF RAIN PER SEC COLLECTED BY GRAUPEL/HAIL
1765             PRACG(K) = CONS41*(((1.2*UMR-0.95*UMG)**2+                   &
1766                   0.08*UMG*UMR)**0.5*RHO(K)*                      &
1767                   N0RR(K)*N0G(K)/LAMR(K)**3*                              &
1768                   (5./(LAMR(K)**3*LAMG(K))+                    &
1769                   2./(LAMR(K)**2*LAMG(K)**2)+                              &
1770                                   0.5/(LAMR(k)*LAMG(k)**3)))
1771
1772 ! ASSUME 1 MM DROPS ARE SHED, GET NUMBER SHED PER SEC
1773
1774             DUM = PRACG(K)/5.2E-7
1775
1776             NPRACG(K) = CONS32*RHO(K)*(1.7*(UNR-UNG)**2+            &
1777                 0.3*UNR*UNG)**0.5*N0RR(K)*N0G(K)*              &
1778                 (1./(LAMR(K)**3*LAMG(K))+                      &
1779                  1./(LAMR(K)**2*LAMG(K)**2)+                   &
1780                  1./(LAMR(K)*LAMG(K)**3))
1781
1782 ! hm 7/15/13, remove limit so that the number of collected drops can smaller than
1783 ! number of shed drops
1784 !            NPRACG(K)=MAX(NPRACG(K)-DUM,0.)
1785             NPRACG(K)=NPRACG(K)-DUM
1786
1787             END IF
1788
1789 !.......................................................................
1790 ! ACCRETION OF CLOUD LIQUID WATER BY RAIN
1791 ! CONTINUOUS COLLECTION EQUATION WITH
1792 ! GRAVITATIONAL COLLECTION KERNEL, DROPLET FALL SPEED NEGLECTED
1793
1794          IF (QR3D(K).GE.1.E-8 .AND. QC3D(K).GE.1.E-8) THEN
1795
1796 ! 12/13/06 HM ADD, REPLACE WITH NEWER FORMULA FROM
1797 ! KHAIROUTDINOV AND KOGAN 2000, MWR
1798
1799            DUM=(QC3D(K)*QR3D(K))
1800            PRA(K) = 67.*(DUM)**1.15
1801            NPRA(K) = PRA(K)/(QC3D(K)/NC3D(K))
1802
1803          END IF
1804 !.......................................................................
1805 ! SELF-COLLECTION OF RAIN DROPS
1806 ! FROM BEHENG(1994)
1807 ! FROM NUMERICAL SIMULATION OF THE STOCHASTIC COLLECTION EQUATION
1808 ! AS DESCRINED ABOVE FOR AUTOCONVERSION
1809
1810          IF (QR3D(K).GE.1.E-8) THEN
1811 ! include breakup add 10/09/09
1812             dum1=300.e-6
1813             if (1./lamr(k).lt.dum1) then
1814             dum=1.
1815             else if (1./lamr(k).ge.dum1) then
1816             dum=2.-exp(2300.*(1./lamr(k)-dum1))
1817             end if
1818 !            NRAGG(K) = -8.*NR3D(K)*QR3D(K)*RHO(K)
1819             NRAGG(K) = -5.78*dum*NR3D(K)*QR3D(K)*RHO(K)
1820          END IF
1821
1822 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1823 ! CALCULATE EVAP OF RAIN (RUTLEDGE AND HOBBS 1983)
1824
1825       IF (QR3D(K).GE.QSMALL) THEN
1826         EPSR = 2.*PI*N0RR(K)*RHO(K)*DV(K)*                           &
1827                    (F1R/(LAMR(K)*LAMR(K))+                       &
1828                     F2R*(ARN(K)*RHO(K)/MU(K))**0.5*                      &
1829                     SC(K)**(1./3.)*CONS9/                   &
1830                 (LAMR(K)**CONS34))
1831       ELSE
1832       EPSR = 0.
1833       END IF
1834
1835 ! NO CONDENSATION ONTO RAIN, ONLY EVAP ALLOWED
1836
1837            IF (QV3D(K).LT.QVS(K)) THEN
1838               PRE(K) = EPSR*(QV3D(K)-QVS(K))/AB(K)
1839               PRE(K) = MIN(PRE(K),0.)
1840            ELSE
1841               PRE(K) = 0.
1842            END IF
1843
1844 !.......................................................................
1845 ! MELTING OF SNOW
1846
1847 ! SNOW MAY PERSITS ABOVE FREEZING, FORMULA FROM RUTLEDGE AND HOBBS, 1984
1848 ! IF WATER SUPERSATURATION, SNOW MELTS TO FORM RAIN
1849
1850           IF (QNI3D(K).GE.1.E-8) THEN
1851
1852 ! fix 053011
1853 ! HM, MODIFY FOR V3.2, ADD ACCELERATED MELTING DUE TO COLLISION WITH RAIN
1854 !             DUM = -CPW/XLF(K)*T3D(K)*PRACS(K)
1855              DUM = -CPW/XLF(K)*(T3D(K)-273.15)*PRACS(K)
1856
1857 ! hm fix 1/20/15
1858 !             PSMLT(K)=2.*PI*N0S(K)*KAP(K)*(273.15-T3D(K))/       &
1859 !                    XLF(K)*RHO(K)*(F1S/(LAMS(K)*LAMS(K))+        &
1860 !                    F2S*(ASN(K)*RHO(K)/MU(K))**0.5*                      &
1861 !                    SC(K)**(1./3.)*CONS10/                   &
1862 !                   (LAMS(K)**CONS35))+DUM
1863              PSMLT(K)=2.*PI*N0S(K)*KAP(K)*(273.15-T3D(K))/       &
1864                     XLF(K)*(F1S/(LAMS(K)*LAMS(K))+        &
1865                     F2S*(ASN(K)*RHO(K)/MU(K))**0.5*                      &
1866                     SC(K)**(1./3.)*CONS10/                   &
1867                    (LAMS(K)**CONS35))+DUM
1868
1869 ! IN WATER SUBSATURATION, SNOW MELTS AND EVAPORATES
1870
1871       IF (QVQVS(K).LT.1.) THEN
1872         EPSS = 2.*PI*N0S(K)*RHO(K)*DV(K)*                            &
1873                    (F1S/(LAMS(K)*LAMS(K))+                       &
1874                     F2S*(ASN(K)*RHO(K)/MU(K))**0.5*                      &
1875                     SC(K)**(1./3.)*CONS10/                   &
1876                (LAMS(K)**CONS35))
1877 ! hm fix 8/4/08
1878         EVPMS(K) = (QV3D(K)-QVS(K))*EPSS/AB(K)
1879         EVPMS(K) = MAX(EVPMS(K),PSMLT(K))
1880         PSMLT(K) = PSMLT(K)-EVPMS(K)
1881       END IF
1882       END IF
1883
1884 !.......................................................................
1885 ! MELTING OF GRAUPEL
1886
1887 ! GRAUPEL MAY PERSITS ABOVE FREEZING, FORMULA FROM RUTLEDGE AND HOBBS, 1984
1888 ! IF WATER SUPERSATURATION, GRAUPEL MELTS TO FORM RAIN
1889
1890           IF (QG3D(K).GE.1.E-8) THEN
1891
1892 ! fix 053011
1893 ! HM, MODIFY FOR V3.2, ADD ACCELERATED MELTING DUE TO COLLISION WITH RAIN
1894 !             DUM = -CPW/XLF(K)*T3D(K)*PRACG(K)
1895              DUM = -CPW/XLF(K)*(T3D(K)-273.15)*PRACG(K)
1896
1897 ! hm fix 1/20/15
1898 !             PGMLT(K)=2.*PI*N0G(K)*KAP(K)*(273.15-T3D(K))/              &
1899 !                    XLF(K)*RHO(K)*(F1S/(LAMG(K)*LAMG(K))+                &
1900 !                    F2S*(AGN(K)*RHO(K)/MU(K))**0.5*                      &
1901 !                    SC(K)**(1./3.)*CONS11/                   &
1902 !                   (LAMG(K)**CONS36))+DUM
1903              PGMLT(K)=2.*PI*N0G(K)*KAP(K)*(273.15-T3D(K))/               &
1904                     XLF(K)*(F1S/(LAMG(K)*LAMG(K))+                &
1905                     F2S*(AGN(K)*RHO(K)/MU(K))**0.5*                      &
1906                     SC(K)**(1./3.)*CONS11/                   &
1907                    (LAMG(K)**CONS36))+DUM
1908
1909 ! IN WATER SUBSATURATION, GRAUPEL MELTS AND EVAPORATES
1910
1911       IF (QVQVS(K).LT.1.) THEN
1912         EPSG = 2.*PI*N0G(K)*RHO(K)*DV(K)*                                &
1913                    (F1S/(LAMG(K)*LAMG(K))+                               &
1914                     F2S*(AGN(K)*RHO(K)/MU(K))**0.5*                      &
1915                     SC(K)**(1./3.)*CONS11/                   &
1916                (LAMG(K)**CONS36))
1917 ! hm fix 8/4/08
1918         EVPMG(K) = (QV3D(K)-QVS(K))*EPSG/AB(K)
1919         EVPMG(K) = MAX(EVPMG(K),PGMLT(K))
1920         PGMLT(K) = PGMLT(K)-EVPMG(K)
1921       END IF
1922       END IF
1923
1924 ! HM, V3.2
1925 ! RESET PRACG AND PRACS TO ZERO, THIS IS DONE BECAUSE THERE IS NO
1926 ! TRANSFER OF MASS FROM SNOW AND GRAUPEL TO RAIN DIRECTLY FROM COLLECTION
1927 ! ABOVE FREEZING, IT IS ONLY USED FOR ENHANCEMENT OF MELTING AND SHEDDING
1928
1929       PRACG(K) = 0.
1930       PRACS(K) = 0.
1931
1932 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1933
1934 ! FOR CLOUD ICE, ONLY PROCESSES OPERATING AT T > 273.15 IS
1935 ! MELTING, WHICH IS ALREADY CONSERVED DURING PROCESS
1936 ! CALCULATION
1937
1938 ! CONSERVATION OF QC
1939
1940       DUM = (PRC(K)+PRA(K))*DT
1941
1942       IF (DUM.GT.QC3D(K).AND.QC3D(K).GE.QSMALL) THEN
1943
1944         RATIO = QC3D(K)/DUM
1945
1946         PRC(K) = PRC(K)*RATIO
1947         PRA(K) = PRA(K)*RATIO
1948
1949         END IF
1950
1951 ! CONSERVATION OF SNOW
1952
1953         DUM = (-PSMLT(K)-EVPMS(K)+PRACS(K))*DT
1954
1955         IF (DUM.GT.QNI3D(K).AND.QNI3D(K).GE.QSMALL) THEN
1956
1957 ! NO SOURCE TERMS FOR SNOW AT T > FREEZING
1958         RATIO = QNI3D(K)/DUM
1959
1960         PSMLT(K) = PSMLT(K)*RATIO
1961         EVPMS(K) = EVPMS(K)*RATIO
1962         PRACS(K) = PRACS(K)*RATIO
1963
1964         END IF
1965
1966 ! CONSERVATION OF GRAUPEL
1967
1968         DUM = (-PGMLT(K)-EVPMG(K)+PRACG(K))*DT
1969
1970         IF (DUM.GT.QG3D(K).AND.QG3D(K).GE.QSMALL) THEN
1971
1972 ! NO SOURCE TERM FOR GRAUPEL ABOVE FREEZING
1973         RATIO = QG3D(K)/DUM
1974
1975         PGMLT(K) = PGMLT(K)*RATIO
1976         EVPMG(K) = EVPMG(K)*RATIO
1977         PRACG(K) = PRACG(K)*RATIO
1978
1979         END IF
1980
1981 ! CONSERVATION OF QR
1982 ! HM 12/13/06, ADDED CONSERVATION OF RAIN SINCE PRE IS NEGATIVE
1983
1984         DUM = (-PRACS(K)-PRACG(K)-PRE(K)-PRA(K)-PRC(K)+PSMLT(K)+PGMLT(K))*DT
1985
1986         IF (DUM.GT.QR3D(K).AND.QR3D(K).GE.QSMALL) THEN
1987
1988         RATIO = (QR3D(K)/DT+PRACS(K)+PRACG(K)+PRA(K)+PRC(K)-PSMLT(K)-PGMLT(K))/ &
1989                         (-PRE(K))
1990         PRE(K) = PRE(K)*RATIO
1991
1992         END IF
1993
1994 !....................................
1995
1996       QV3DTEN(K) = QV3DTEN(K)+(-PRE(K)-EVPMS(K)-EVPMG(K))
1997
1998       T3DTEN(K) = T3DTEN(K)+(PRE(K)*XXLV(K)+(EVPMS(K)+EVPMG(K))*XXLS(K)+&
1999                     (PSMLT(K)+PGMLT(K)-PRACS(K)-PRACG(K))*XLF(K))/CPM(K)
2000
2001       QC3DTEN(K) = QC3DTEN(K)+(-PRA(K)-PRC(K))
2002       QR3DTEN(K) = QR3DTEN(K)+(PRE(K)+PRA(K)+PRC(K)-PSMLT(K)-PGMLT(K)+PRACS(K)+PRACG(K))
2003       QNI3DTEN(K) = QNI3DTEN(K)+(PSMLT(K)+EVPMS(K)-PRACS(K))
2004       QG3DTEN(K) = QG3DTEN(K)+(PGMLT(K)+EVPMG(K)-PRACG(K))
2005 ! fix 053011
2006 !      NS3DTEN(K) = NS3DTEN(K)-NPRACS(K)
2007 ! HM, bug fix 5/12/08, npracg is subtracted from nr not ng
2008 !      NG3DTEN(K) = NG3DTEN(K)
2009       NC3DTEN(K) = NC3DTEN(K)+ (-NPRA(K)-NPRC(K))
2010       NR3DTEN(K) = NR3DTEN(K)+ (NPRC1(K)+NRAGG(K)-NPRACG(K))
2011
2012 ! HM ADD, WRF-CHEM, ADD TENDENCIES FOR C2PREC
2013
2014         C2PREC(K) = PRA(K)+PRC(K)
2015       IF (PRE(K).LT.0.) THEN
2016          DUM = PRE(K)*DT/QR3D(K)
2017            DUM = MAX(-1.,DUM)
2018          NSUBR(K) = DUM*NR3D(K)/DT
2019       END IF
2020
2021         IF (EVPMS(K)+PSMLT(K).LT.0.) THEN
2022          DUM = (EVPMS(K)+PSMLT(K))*DT/QNI3D(K)
2023            DUM = MAX(-1.,DUM)
2024          NSMLTS(K) = DUM*NS3D(K)/DT
2025         END IF
2026         IF (PSMLT(K).LT.0.) THEN
2027           DUM = PSMLT(K)*DT/QNI3D(K)
2028           DUM = MAX(-1.0,DUM)
2029           NSMLTR(K) = DUM*NS3D(K)/DT
2030         END IF
2031         IF (EVPMG(K)+PGMLT(K).LT.0.) THEN
2032          DUM = (EVPMG(K)+PGMLT(K))*DT/QG3D(K)
2033            DUM = MAX(-1.,DUM)
2034          NGMLTG(K) = DUM*NG3D(K)/DT
2035         END IF
2036         IF (PGMLT(K).LT.0.) THEN
2037           DUM = PGMLT(K)*DT/QG3D(K)
2038           DUM = MAX(-1.0,DUM)
2039           NGMLTR(K) = DUM*NG3D(K)/DT
2040         END IF
2041
2042          NS3DTEN(K) = NS3DTEN(K)+(NSMLTS(K))
2043          NG3DTEN(K) = NG3DTEN(K)+(NGMLTG(K))
2044          NR3DTEN(K) = NR3DTEN(K)+(NSUBR(K)-NSMLTR(K)-NGMLTR(K))
2045
2046  300  CONTINUE
2047
2048 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
2049 ! NOW CALCULATE SATURATION ADJUSTMENT TO CONDENSE EXTRA VAPOR ABOVE
2050 ! WATER SATURATION
2051
2052       DUMT = T3D(K)+DT*T3DTEN(K)
2053       DUMQV = QV3D(K)+DT*QV3DTEN(K)
2054 ! hm, add fix for low pressure, 5/12/10
2055       dum=min(0.99*pres(k),POLYSVP(DUMT,0))
2056       DUMQSS = EP_2*dum/(PRES(K)-dum)
2057       DUMQC = QC3D(K)+DT*QC3DTEN(K)
2058       DUMQC = MAX(DUMQC,0.)
2059
2060 ! SATURATION ADJUSTMENT FOR LIQUID
2061
2062       DUMS = DUMQV-DUMQSS
2063       PCC(K) = DUMS/(1.+XXLV(K)**2*DUMQSS/(CPM(K)*RV*DUMT**2))/DT
2064       IF (PCC(K)*DT+DUMQC.LT.0.) THEN
2065            PCC(K) = -DUMQC/DT
2066       END IF
2067
2068       QV3DTEN(K) = QV3DTEN(K)-PCC(K)
2069       T3DTEN(K) = T3DTEN(K)+PCC(K)*XXLV(K)/CPM(K)
2070       QC3DTEN(K) = QC3DTEN(K)+PCC(K)
2071
2072 #if (WRF_CHEM == 1)
2073       IF( has_wetscav ) THEN
2074          evapprod(k) = - PRE(K) - EVPMS(K) - EVPMG(K)
2075          rainprod(k) = PRA(K) + PRC(K) + tqimelt(K)
2076       END IF
2077 #endif
2078
2079 !.......................................................................
2080 ! ACTIVATION OF CLOUD DROPLETS
2081 ! ACTIVATION OF DROPLET CURRENTLY NOT CALCULATED
2082 ! DROPLET CONCENTRATION IS SPECIFIED !!!!!
2083
2084 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
2085 ! SUBLIMATE, MELT, OR EVAPORATE NUMBER CONCENTRATION
2086 ! THIS FORMULATION ASSUMES 1:1 RATIO BETWEEN MASS LOSS AND
2087 ! LOSS OF NUMBER CONCENTRATION
2088
2089 !     IF (PCC(K).LT.0.) THEN
2090 !        DUM = PCC(K)*DT/QC3D(K)
2091 !           DUM = MAX(-1.,DUM)
2092 !        NSUBC(K) = DUM*NC3D(K)/DT
2093 !     END IF
2094
2095 ! UPDATE TENDENCIES
2096
2097 !        NC3DTEN(K) = NC3DTEN(K)+NSUBC(K)
2098
2099 !.....................................................................
2100 !.....................................................................
2101          ELSE  ! TEMPERATURE < 273.15
2102
2103 !......................................................................
2104 !HM ADD, ALLOW FOR CONSTANT DROPLET NUMBER
2105 ! INUM = 0, PREDICT DROPLET NUMBER
2106 ! INUM = 1, SET CONSTANT DROPLET NUMBER
2107
2108          IF (iinum.EQ.1) THEN
2109 ! CONVERT NDCNST FROM CM-3 TO KG-1
2110             NC3D(K)=NDCNST*1.E6/RHO(K)
2111          END IF
2112
2113 ! CALCULATE SIZE DISTRIBUTION PARAMETERS
2114 ! MAKE SURE NUMBER CONCENTRATIONS AREN'T NEGATIVE
2115
2116       NI3D(K) = MAX(0.,NI3D(K))
2117       NS3D(K) = MAX(0.,NS3D(K))
2118       NC3D(K) = MAX(0.,NC3D(K))
2119       NR3D(K) = MAX(0.,NR3D(K))
2120       NG3D(K) = MAX(0.,NG3D(K))
2121
2122 !......................................................................
2123 ! CLOUD ICE
2124
2125       IF (QI3D(K).GE.QSMALL) THEN
2126          LAMI(K) = (CONS12*                 &
2127               NI3D(K)/QI3D(K))**(1./DI)
2128          N0I(K) = NI3D(K)*LAMI(K)
2129
2130 ! CHECK FOR SLOPE
2131
2132 ! ADJUST VARS
2133
2134       IF (LAMI(K).LT.LAMMINI) THEN
2135
2136       LAMI(K) = LAMMINI
2137
2138       N0I(K) = LAMI(K)**4*QI3D(K)/CONS12
2139
2140       NI3D(K) = N0I(K)/LAMI(K)
2141       ELSE IF (LAMI(K).GT.LAMMAXI) THEN
2142       LAMI(K) = LAMMAXI
2143       N0I(K) = LAMI(K)**4*QI3D(K)/CONS12
2144
2145       NI3D(K) = N0I(K)/LAMI(K)
2146       END IF
2147       END IF
2148
2149 !......................................................................
2150 ! RAIN
2151
2152       IF (QR3D(K).GE.QSMALL) THEN
2153       LAMR(K) = (PI*RHOW*NR3D(K)/QR3D(K))**(1./3.)
2154       N0RR(K) = NR3D(K)*LAMR(K)
2155
2156 ! CHECK FOR SLOPE
2157
2158 ! ADJUST VARS
2159
2160       IF (LAMR(K).LT.LAMMINR) THEN
2161
2162       LAMR(K) = LAMMINR
2163
2164       N0RR(K) = LAMR(K)**4*QR3D(K)/(PI*RHOW)
2165
2166       NR3D(K) = N0RR(K)/LAMR(K)
2167       ELSE IF (LAMR(K).GT.LAMMAXR) THEN
2168       LAMR(K) = LAMMAXR
2169       N0RR(K) = LAMR(K)**4*QR3D(K)/(PI*RHOW)
2170
2171       NR3D(K) = N0RR(K)/LAMR(K)
2172       END IF
2173       END IF
2174
2175 !......................................................................
2176 ! CLOUD DROPLETS
2177
2178 ! MARTIN ET AL. (1994) FORMULA FOR PGAM
2179
2180       IF (QC3D(K).GE.QSMALL) THEN
2181
2182          DUM = PRES(K)/(287.15*T3D(K))
2183          PGAM(K)=0.0005714*(NC3D(K)/1.E6*DUM)+0.2714
2184          PGAM(K)=1./(PGAM(K)**2)-1.
2185          PGAM(K)=MAX(PGAM(K),2.)
2186          PGAM(K)=MIN(PGAM(K),10.)
2187
2188 ! CALCULATE LAMC
2189
2190       LAMC(K) = (CONS26*NC3D(K)*GAMMA(PGAM(K)+4.)/   &
2191                  (QC3D(K)*GAMMA(PGAM(K)+1.)))**(1./3.)
2192
2193 ! LAMMIN, 60 MICRON DIAMETER
2194 ! LAMMAX, 1 MICRON
2195
2196       LAMMIN = (PGAM(K)+1.)/60.E-6
2197       LAMMAX = (PGAM(K)+1.)/1.E-6
2198
2199       IF (LAMC(K).LT.LAMMIN) THEN
2200       LAMC(K) = LAMMIN
2201
2202       NC3D(K) = EXP(3.*LOG(LAMC(K))+LOG(QC3D(K))+              &
2203                 LOG(GAMMA(PGAM(K)+1.))-LOG(GAMMA(PGAM(K)+4.)))/CONS26
2204       ELSE IF (LAMC(K).GT.LAMMAX) THEN
2205       LAMC(K) = LAMMAX
2206       NC3D(K) = EXP(3.*LOG(LAMC(K))+LOG(QC3D(K))+              &
2207                 LOG(GAMMA(PGAM(K)+1.))-LOG(GAMMA(PGAM(K)+4.)))/CONS26
2208
2209       END IF
2210
2211 ! TO CALCULATE DROPLET FREEZING
2212
2213         CDIST1(K) = NC3D(K)/GAMMA(PGAM(K)+1.)
2214
2215       END IF
2216
2217 !......................................................................
2218 ! SNOW
2219
2220       IF (QNI3D(K).GE.QSMALL) THEN
2221       LAMS(K) = (CONS1*NS3D(K)/QNI3D(K))**(1./DS)
2222       N0S(K) = NS3D(K)*LAMS(K)
2223
2224 ! CHECK FOR SLOPE
2225
2226 ! ADJUST VARS
2227
2228       IF (LAMS(K).LT.LAMMINS) THEN
2229       LAMS(K) = LAMMINS
2230       N0S(K) = LAMS(K)**4*QNI3D(K)/CONS1
2231
2232       NS3D(K) = N0S(K)/LAMS(K)
2233
2234       ELSE IF (LAMS(K).GT.LAMMAXS) THEN
2235
2236       LAMS(K) = LAMMAXS
2237       N0S(K) = LAMS(K)**4*QNI3D(K)/CONS1
2238
2239       NS3D(K) = N0S(K)/LAMS(K)
2240       END IF
2241       END IF
2242
2243 !......................................................................
2244 ! GRAUPEL
2245
2246       IF (QG3D(K).GE.QSMALL) THEN
2247       LAMG(K) = (CONS2*NG3D(K)/QG3D(K))**(1./DG)
2248       N0G(K) = NG3D(K)*LAMG(K)
2249
2250 ! CHECK FOR SLOPE
2251
2252 ! ADJUST VARS
2253
2254       IF (LAMG(K).LT.LAMMING) THEN
2255       LAMG(K) = LAMMING
2256       N0G(K) = LAMG(K)**4*QG3D(K)/CONS2
2257
2258       NG3D(K) = N0G(K)/LAMG(K)
2259
2260       ELSE IF (LAMG(K).GT.LAMMAXG) THEN
2261
2262       LAMG(K) = LAMMAXG
2263       N0G(K) = LAMG(K)**4*QG3D(K)/CONS2
2264
2265       NG3D(K) = N0G(K)/LAMG(K)
2266       END IF
2267       END IF
2268
2269 !.....................................................................
2270 ! ZERO OUT PROCESS RATES
2271
2272             MNUCCC(K) = 0.
2273             NNUCCC(K) = 0.
2274             PRC(K) = 0.
2275             NPRC(K) = 0.
2276             NPRC1(K) = 0.
2277             NSAGG(K) = 0.
2278             PSACWS(K) = 0.
2279             NPSACWS(K) = 0.
2280             PSACWI(K) = 0.
2281             NPSACWI(K) = 0.
2282             PRACS(K) = 0.
2283             NPRACS(K) = 0.
2284             NMULTS(K) = 0.
2285             QMULTS(K) = 0.
2286             NMULTR(K) = 0.
2287             QMULTR(K) = 0.
2288             NMULTG(K) = 0.
2289             QMULTG(K) = 0.
2290             NMULTRG(K) = 0.
2291             QMULTRG(K) = 0.
2292             MNUCCR(K) = 0.
2293             NNUCCR(K) = 0.
2294             PRA(K) = 0.
2295             NPRA(K) = 0.
2296             NRAGG(K) = 0.
2297             PRCI(K) = 0.
2298             NPRCI(K) = 0.
2299             PRAI(K) = 0.
2300             NPRAI(K) = 0.
2301             NNUCCD(K) = 0.
2302             MNUCCD(K) = 0.
2303             PCC(K) = 0.
2304             PRE(K) = 0.
2305             PRD(K) = 0.
2306             PRDS(K) = 0.
2307             EPRD(K) = 0.
2308             EPRDS(K) = 0.
2309             NSUBC(K) = 0.
2310             NSUBI(K) = 0.
2311             NSUBS(K) = 0.
2312             NSUBR(K) = 0.
2313             PIACR(K) = 0.
2314             NIACR(K) = 0.
2315             PRACI(K) = 0.
2316             PIACRS(K) = 0.
2317             NIACRS(K) = 0.
2318             PRACIS(K) = 0.
2319 ! HM: ADD GRAUPEL PROCESSES
2320             PRACG(K) = 0.
2321             PSACR(K) = 0.
2322             PSACWG(K) = 0.
2323             PGSACW(K) = 0.
2324             PGRACS(K) = 0.
2325             PRDG(K) = 0.
2326             EPRDG(K) = 0.
2327             NPRACG(K) = 0.
2328             NPSACWG(K) = 0.
2329             NSCNG(K) = 0.
2330             NGRACS(K) = 0.
2331             NSUBG(K) = 0.
2332
2333 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
2334 ! CALCULATION OF MICROPHYSICAL PROCESS RATES
2335 ! ACCRETION/AUTOCONVERSION/FREEZING/MELTING/COAG.
2336 !.......................................................................
2337 ! FREEZING OF CLOUD DROPLETS
2338 ! ONLY ALLOWED BELOW -4 C
2339         IF (QC3D(K).GE.QSMALL .AND. T3D(K).LT.269.15) THEN
2340
2341 ! NUMBER OF CONTACT NUCLEI (M^-3) FROM MEYERS ET AL., 1992
2342 ! FACTOR OF 1000 IS TO CONVERT FROM L^-1 TO M^-3
2343
2344 ! MEYERS CURVE
2345
2346            NACNT = EXP(-2.80+0.262*(273.15-T3D(K)))*1000.
2347
2348 ! COOPER CURVE
2349 !        NACNT =  5.*EXP(0.304*(273.15-T3D(K)))
2350
2351 ! FLECTHER
2352 !     NACNT = 0.01*EXP(0.6*(273.15-T3D(K)))
2353
2354 ! CONTACT FREEZING
2355
2356 ! MEAN FREE PATH
2357
2358             DUM = 7.37*T3D(K)/(288.*10.*PRES(K))/100.
2359
2360 ! EFFECTIVE DIFFUSIVITY OF CONTACT NUCLEI
2361 ! BASED ON BROWNIAN DIFFUSION
2362
2363             DAP(K) = CONS37*T3D(K)*(1.+DUM/RIN)/MU(K)
2364
2365            MNUCCC(K) = CONS38*DAP(K)*NACNT*EXP(LOG(CDIST1(K))+   &
2366                    LOG(GAMMA(PGAM(K)+5.))-4.*LOG(LAMC(K)))
2367            NNUCCC(K) = 2.*PI*DAP(K)*NACNT*CDIST1(K)*           &
2368                     GAMMA(PGAM(K)+2.)/                         &
2369                     LAMC(K)
2370
2371 ! IMMERSION FREEZING (BIGG 1953)
2372
2373 !           MNUCCC(K) = MNUCCC(K)+CONS39*                   &
2374 !                  EXP(LOG(CDIST1(K))+LOG(GAMMA(7.+PGAM(K)))-6.*LOG(LAMC(K)))*             &
2375 !                   EXP(AIMM*(273.15-T3D(K)))
2376
2377 !           NNUCCC(K) = NNUCCC(K)+                                  &
2378 !            CONS40*EXP(LOG(CDIST1(K))+LOG(GAMMA(PGAM(K)+4.))-3.*LOG(LAMC(K)))              &
2379 !                *EXP(AIMM*(273.15-T3D(K)))
2380
2381 ! hm 7/15/13 fix for consistency w/ original formula
2382            MNUCCC(K) = MNUCCC(K)+CONS39*                   &
2383                   EXP(LOG(CDIST1(K))+LOG(GAMMA(7.+PGAM(K)))-6.*LOG(LAMC(K)))*             &
2384                    (EXP(AIMM*(273.15-T3D(K)))-1.)
2385
2386            NNUCCC(K) = NNUCCC(K)+                                  &
2387             CONS40*EXP(LOG(CDIST1(K))+LOG(GAMMA(PGAM(K)+4.))-3.*LOG(LAMC(K)))              &
2388                 *(EXP(AIMM*(273.15-T3D(K)))-1.)
2389
2390 ! PUT IN A CATCH HERE TO PREVENT DIVERGENCE BETWEEN NUMBER CONC. AND
2391 ! MIXING RATIO, SINCE STRICT CONSERVATION NOT CHECKED FOR NUMBER CONC
2392
2393            NNUCCC(K) = MIN(NNUCCC(K),NC3D(K)/DT)
2394
2395         END IF
2396
2397 !.................................................................
2398 !.......................................................................
2399 ! AUTOCONVERSION OF CLOUD LIQUID WATER TO RAIN
2400 ! FORMULA FROM BEHENG (1994)
2401 ! USING NUMERICAL SIMULATION OF STOCHASTIC COLLECTION EQUATION
2402 ! AND INITIAL CLOUD DROPLET SIZE DISTRIBUTION SPECIFIED
2403 ! AS A GAMMA DISTRIBUTION
2404
2405 ! USE MINIMUM VALUE OF 1.E-6 TO PREVENT FLOATING POINT ERROR
2406
2407          IF (QC3D(K).GE.1.E-6) THEN
2408
2409 ! HM ADD 12/13/06, REPLACE WITH NEWER FORMULA
2410 ! FROM KHAIROUTDINOV AND KOGAN 2000, MWR
2411
2412                 PRC(K)=1350.*QC3D(K)**2.47*  &
2413            (NC3D(K)/1.e6*RHO(K))**(-1.79)
2414
2415 ! note: nprc1 is change in Nr,
2416 ! nprc is change in Nc
2417
2418         NPRC1(K) = PRC(K)/CONS29
2419         NPRC(K) = PRC(K)/(QC3D(K)/NC3D(K))
2420
2421 ! hm bug fix 3/20/12
2422                 NPRC(K) = MIN(NPRC(K),NC3D(K)/DT)
2423                 NPRC1(K) = MIN(NPRC1(K),NPRC(K))
2424
2425          END IF
2426
2427 !.......................................................................
2428 ! SELF-COLLECTION OF DROPLET NOT INCLUDED IN KK2000 SCHEME
2429
2430 ! SNOW AGGREGATION FROM PASSARELLI, 1978, USED BY REISNER, 1998
2431 ! THIS IS HARD-WIRED FOR BS = 0.4 FOR NOW
2432
2433          IF (QNI3D(K).GE.1.E-8) THEN
2434              NSAGG(K) = CONS15*ASN(K)*RHO(K)**            &
2435             ((2.+BS)/3.)*QNI3D(K)**((2.+BS)/3.)*                  &
2436             (NS3D(K)*RHO(K))**((4.-BS)/3.)/                       &
2437             (RHO(K))
2438          END IF
2439
2440 !.......................................................................
2441 ! ACCRETION OF CLOUD DROPLETS ONTO SNOW/GRAUPEL
2442 ! HERE USE CONTINUOUS COLLECTION EQUATION WITH
2443 ! SIMPLE GRAVITATIONAL COLLECTION KERNEL IGNORING
2444
2445 ! SNOW
2446
2447          IF (QNI3D(K).GE.1.E-8 .AND. QC3D(K).GE.QSMALL) THEN
2448
2449            PSACWS(K) = CONS13*ASN(K)*QC3D(K)*RHO(K)*               &
2450                   N0S(K)/                        &
2451                   LAMS(K)**(BS+3.)
2452            NPSACWS(K) = CONS13*ASN(K)*NC3D(K)*RHO(K)*              &
2453                   N0S(K)/                        &
2454                   LAMS(K)**(BS+3.)
2455
2456          END IF
2457
2458 !............................................................................
2459 ! COLLECTION OF CLOUD WATER BY GRAUPEL
2460
2461          IF (QG3D(K).GE.1.E-8 .AND. QC3D(K).GE.QSMALL) THEN
2462
2463            PSACWG(K) = CONS14*AGN(K)*QC3D(K)*RHO(K)*               &
2464                   N0G(K)/                        &
2465                   LAMG(K)**(BG+3.)
2466            NPSACWG(K) = CONS14*AGN(K)*NC3D(K)*RHO(K)*              &
2467                   N0G(K)/                        &
2468                   LAMG(K)**(BG+3.)
2469             END IF
2470
2471 !.......................................................................
2472 ! HM, ADD 12/13/06
2473 ! CLOUD ICE COLLECTING DROPLETS, ASSUME THAT CLOUD ICE MEAN DIAM > 100 MICRON
2474 ! BEFORE RIMING CAN OCCUR
2475 ! ASSUME THAT RIME COLLECTED ON CLOUD ICE DOES NOT LEAD
2476 ! TO HALLET-MOSSOP SPLINTERING
2477
2478          IF (QI3D(K).GE.1.E-8 .AND. QC3D(K).GE.QSMALL) THEN
2479
2480 ! PUT IN SIZE DEPENDENT COLLECTION EFFICIENCY BASED ON STOKES LAW
2481 ! FROM THOMPSON ET AL. 2004, MWR
2482
2483             IF (1./LAMI(K).GE.100.E-6) THEN
2484
2485            PSACWI(K) = CONS16*AIN(K)*QC3D(K)*RHO(K)*               &
2486                   N0I(K)/                        &
2487                   LAMI(K)**(BI+3.)
2488            NPSACWI(K) = CONS16*AIN(K)*NC3D(K)*RHO(K)*              &
2489                   N0I(K)/                        &
2490                   LAMI(K)**(BI+3.)
2491            END IF
2492          END IF
2493
2494 !.......................................................................
2495 ! ACCRETION OF RAIN WATER BY SNOW
2496 ! FORMULA FROM IKAWA AND SAITO, 1991, USED BY REISNER ET AL, 1998
2497
2498          IF (QR3D(K).GE.1.E-8.AND.QNI3D(K).GE.1.E-8) THEN
2499
2500             UMS = ASN(K)*CONS3/(LAMS(K)**BS)
2501             UMR = ARN(K)*CONS4/(LAMR(K)**BR)
2502             UNS = ASN(K)*CONS5/LAMS(K)**BS
2503             UNR = ARN(K)*CONS6/LAMR(K)**BR
2504
2505 ! SET REASLISTIC LIMITS ON FALLSPEEDS
2506
2507 ! bug fix, 10/08/09
2508             dum=(rhosu/rho(k))**0.54
2509             UMS=MIN(UMS,1.2*dum)
2510             UNS=MIN(UNS,1.2*dum)
2511             UMR=MIN(UMR,9.1*dum)
2512             UNR=MIN(UNR,9.1*dum)
2513
2514             PRACS(K) = CONS41*(((1.2*UMR-0.95*UMS)**2+                   &
2515                   0.08*UMS*UMR)**0.5*RHO(K)*                      &
2516                   N0RR(K)*N0S(K)/LAMR(K)**3*                              &
2517                   (5./(LAMR(K)**3*LAMS(K))+                    &
2518                   2./(LAMR(K)**2*LAMS(K)**2)+                  &
2519                   0.5/(LAMR(k)*LAMS(k)**3)))
2520
2521             NPRACS(K) = CONS32*RHO(K)*(1.7*(UNR-UNS)**2+            &
2522                 0.3*UNR*UNS)**0.5*N0RR(K)*N0S(K)*              &
2523                 (1./(LAMR(K)**3*LAMS(K))+                      &
2524                  1./(LAMR(K)**2*LAMS(K)**2)+                   &
2525                  1./(LAMR(K)*LAMS(K)**3))
2526
2527 ! MAKE SURE PRACS DOESN'T EXCEED TOTAL RAIN MIXING RATIO
2528 ! AS THIS MAY OTHERWISE RESULT IN TOO MUCH TRANSFER OF WATER DURING
2529 ! RIME-SPLINTERING
2530
2531             PRACS(K) = MIN(PRACS(K),QR3D(K)/DT)
2532
2533 ! COLLECTION OF SNOW BY RAIN - NEEDED FOR GRAUPEL CONVERSION CALCULATIONS
2534 ! ONLY CALCULATE IF SNOW AND RAIN MIXING RATIOS EXCEED 0.1 G/KG
2535
2536 ! HM MODIFY FOR WRFV3.1
2537 !            IF (IHAIL.EQ.0) THEN
2538             IF (QNI3D(K).GE.0.1E-3.AND.QR3D(K).GE.0.1E-3) THEN
2539             PSACR(K) = CONS31*(((1.2*UMR-0.95*UMS)**2+              &
2540                   0.08*UMS*UMR)**0.5*RHO(K)*                     &
2541                  N0RR(K)*N0S(K)/LAMS(K)**3*                               &
2542                   (5./(LAMS(K)**3*LAMR(K))+                    &
2543                   2./(LAMS(K)**2*LAMR(K)**2)+                  &
2544                   0.5/(LAMS(K)*LAMR(K)**3)))
2545             END IF
2546 !            END IF
2547
2548          END IF
2549
2550 !.......................................................................
2551
2552 ! COLLECTION OF RAINWATER BY GRAUPEL, FROM IKAWA AND SAITO 1990,
2553 ! USED BY REISNER ET AL 1998
2554          IF (QR3D(K).GE.1.E-8.AND.QG3D(K).GE.1.E-8) THEN
2555
2556             UMG = AGN(K)*CONS7/(LAMG(K)**BG)
2557             UMR = ARN(K)*CONS4/(LAMR(K)**BR)
2558             UNG = AGN(K)*CONS8/LAMG(K)**BG
2559             UNR = ARN(K)*CONS6/LAMR(K)**BR
2560
2561 ! SET REASLISTIC LIMITS ON FALLSPEEDS
2562 ! bug fix, 10/08/09
2563             dum=(rhosu/rho(k))**0.54
2564             UMG=MIN(UMG,20.*dum)
2565             UNG=MIN(UNG,20.*dum)
2566             UMR=MIN(UMR,9.1*dum)
2567             UNR=MIN(UNR,9.1*dum)
2568
2569             PRACG(K) = CONS41*(((1.2*UMR-0.95*UMG)**2+                   &
2570                   0.08*UMG*UMR)**0.5*RHO(K)*                      &
2571                   N0RR(K)*N0G(K)/LAMR(K)**3*                              &
2572                   (5./(LAMR(K)**3*LAMG(K))+                    &
2573                   2./(LAMR(K)**2*LAMG(K)**2)+                              &
2574                                   0.5/(LAMR(k)*LAMG(k)**3)))
2575
2576             NPRACG(K) = CONS32*RHO(K)*(1.7*(UNR-UNG)**2+            &
2577                 0.3*UNR*UNG)**0.5*N0RR(K)*N0G(K)*              &
2578                 (1./(LAMR(K)**3*LAMG(K))+                      &
2579                  1./(LAMR(K)**2*LAMG(K)**2)+                   &
2580                  1./(LAMR(K)*LAMG(K)**3))
2581
2582 ! MAKE SURE PRACG DOESN'T EXCEED TOTAL RAIN MIXING RATIO
2583 ! AS THIS MAY OTHERWISE RESULT IN TOO MUCH TRANSFER OF WATER DURING
2584 ! RIME-SPLINTERING
2585
2586             PRACG(K) = MIN(PRACG(K),QR3D(K)/DT)
2587
2588             END IF
2589
2590 !.......................................................................
2591 ! RIME-SPLINTERING - SNOW
2592 ! HALLET-MOSSOP (1974)
2593 ! NUMBER OF SPLINTERS FORMED IS BASED ON MASS OF RIMED WATER
2594
2595 ! DUM1 = MASS OF INDIVIDUAL SPLINTERS
2596
2597 ! HM ADD THRESHOLD SNOW AND DROPLET MIXING RATIO FOR RIME-SPLINTERING
2598 ! TO LIMIT RIME-SPLINTERING IN STRATIFORM CLOUDS
2599 ! THESE THRESHOLDS CORRESPOND WITH GRAUPEL THRESHOLDS IN RH 1984
2600
2601 !v1.4
2602          IF (QNI3D(K).GE.0.1E-3) THEN
2603          IF (QC3D(K).GE.0.5E-3.OR.QR3D(K).GE.0.1E-3) THEN
2604          IF (PSACWS(K).GT.0..OR.PRACS(K).GT.0.) THEN
2605             IF (T3D(K).LT.270.16 .AND. T3D(K).GT.265.16) THEN
2606
2607                IF (T3D(K).GT.270.16) THEN
2608                   FMULT = 0.
2609                ELSE IF (T3D(K).LE.270.16.AND.T3D(K).GT.268.16)  THEN
2610                   FMULT = (270.16-T3D(K))/2.
2611                ELSE IF (T3D(K).GE.265.16.AND.T3D(K).LE.268.16)   THEN
2612                   FMULT = (T3D(K)-265.16)/3.
2613                ELSE IF (T3D(K).LT.265.16) THEN
2614                   FMULT = 0.
2615                END IF
2616
2617 ! 1000 IS TO CONVERT FROM KG TO G
2618
2619 ! SPLINTERING FROM DROPLETS ACCRETED ONTO SNOW
2620
2621                IF (PSACWS(K).GT.0.) THEN
2622                   NMULTS(K) = 35.E4*PSACWS(K)*FMULT*1000.
2623                   QMULTS(K) = NMULTS(K)*MMULT
2624
2625 ! CONSTRAIN SO THAT TRANSFER OF MASS FROM SNOW TO ICE CANNOT BE MORE MASS
2626 ! THAN WAS RIMED ONTO SNOW
2627
2628                   QMULTS(K) = MIN(QMULTS(K),PSACWS(K))
2629                   PSACWS(K) = PSACWS(K)-QMULTS(K)
2630
2631                END IF
2632
2633 ! RIMING AND SPLINTERING FROM ACCRETED RAINDROPS
2634
2635                IF (PRACS(K).GT.0.) THEN
2636                    NMULTR(K) = 35.E4*PRACS(K)*FMULT*1000.
2637                    QMULTR(K) = NMULTR(K)*MMULT
2638
2639 ! CONSTRAIN SO THAT TRANSFER OF MASS FROM SNOW TO ICE CANNOT BE MORE MASS
2640 ! THAN WAS RIMED ONTO SNOW
2641
2642                    QMULTR(K) = MIN(QMULTR(K),PRACS(K))
2643
2644                    PRACS(K) = PRACS(K)-QMULTR(K)
2645
2646                END IF
2647
2648             END IF
2649          END IF
2650          END IF
2651          END IF
2652
2653 !.......................................................................
2654 ! RIME-SPLINTERING - GRAUPEL
2655 ! HALLET-MOSSOP (1974)
2656 ! NUMBER OF SPLINTERS FORMED IS BASED ON MASS OF RIMED WATER
2657
2658 ! DUM1 = MASS OF INDIVIDUAL SPLINTERS
2659
2660 ! HM ADD THRESHOLD SNOW MIXING RATIO FOR RIME-SPLINTERING
2661 ! TO LIMIT RIME-SPLINTERING IN STRATIFORM CLOUDS
2662
2663 !         IF (IHAIL.EQ.0) THEN
2664 ! v1.4
2665          IF (QG3D(K).GE.0.1E-3) THEN
2666          IF (QC3D(K).GE.0.5E-3.OR.QR3D(K).GE.0.1E-3) THEN
2667          IF (PSACWG(K).GT.0..OR.PRACG(K).GT.0.) THEN
2668             IF (T3D(K).LT.270.16 .AND. T3D(K).GT.265.16) THEN
2669
2670                IF (T3D(K).GT.270.16) THEN
2671                   FMULT = 0.
2672                ELSE IF (T3D(K).LE.270.16.AND.T3D(K).GT.268.16)  THEN
2673                   FMULT = (270.16-T3D(K))/2.
2674                ELSE IF (T3D(K).GE.265.16.AND.T3D(K).LE.268.16)   THEN
2675                   FMULT = (T3D(K)-265.16)/3.
2676                ELSE IF (T3D(K).LT.265.16) THEN
2677                   FMULT = 0.
2678                END IF
2679
2680 ! 1000 IS TO CONVERT FROM KG TO G
2681
2682 ! SPLINTERING FROM DROPLETS ACCRETED ONTO GRAUPEL
2683
2684                IF (PSACWG(K).GT.0.) THEN
2685                   NMULTG(K) = 35.E4*PSACWG(K)*FMULT*1000.
2686                   QMULTG(K) = NMULTG(K)*MMULT
2687
2688 ! CONSTRAIN SO THAT TRANSFER OF MASS FROM GRAUPEL TO ICE CANNOT BE MORE MASS
2689 ! THAN WAS RIMED ONTO GRAUPEL
2690
2691                   QMULTG(K) = MIN(QMULTG(K),PSACWG(K))
2692                   PSACWG(K) = PSACWG(K)-QMULTG(K)
2693
2694                END IF
2695
2696 ! RIMING AND SPLINTERING FROM ACCRETED RAINDROPS
2697
2698                IF (PRACG(K).GT.0.) THEN
2699                    NMULTRG(K) = 35.E4*PRACG(K)*FMULT*1000.
2700                    QMULTRG(K) = NMULTRG(K)*MMULT
2701
2702 ! CONSTRAIN SO THAT TRANSFER OF MASS FROM GRAUPEL TO ICE CANNOT BE MORE MASS
2703 ! THAN WAS RIMED ONTO GRAUPEL
2704
2705                    QMULTRG(K) = MIN(QMULTRG(K),PRACG(K))
2706                    PRACG(K) = PRACG(K)-QMULTRG(K)
2707
2708                END IF
2709                END IF
2710                END IF
2711             END IF
2712             END IF
2713 !         END IF
2714
2715 !........................................................................
2716 ! CONVERSION OF RIMED CLOUD WATER ONTO SNOW TO GRAUPEL/HAIL
2717
2718 !           IF (IHAIL.EQ.0) THEN
2719            IF (PSACWS(K).GT.0.) THEN
2720 ! ONLY ALLOW CONVERSION IF QNI > 0.1 AND QC > 0.5 G/KG FOLLOWING RUTLEDGE AND HOBBS (1984)
2721               IF (QNI3D(K).GE.0.1E-3.AND.QC3D(K).GE.0.5E-3) THEN
2722
2723 ! PORTION OF RIMING CONVERTED TO GRAUPEL (REISNER ET AL. 1998, ORIGINALLY IS1991)
2724              PGSACW(K) = MIN(PSACWS(K),CONS17*DT*N0S(K)*QC3D(K)*QC3D(K)* &
2725                           ASN(K)*ASN(K)/ &
2726                            (RHO(K)*LAMS(K)**(2.*BS+2.)))
2727
2728 ! MIX RAT CONVERTED INTO GRAUPEL AS EMBRYO (REISNER ET AL. 1998, ORIG M1990)
2729              DUM = MAX(RHOSN/(RHOG-RHOSN)*PGSACW(K),0.)
2730
2731 ! NUMBER CONCENTRAITON OF EMBRYO GRAUPEL FROM RIMING OF SNOW
2732              NSCNG(K) = DUM/MG0*RHO(K)
2733 ! LIMIT MAX NUMBER CONVERTED TO SNOW NUMBER
2734              NSCNG(K) = MIN(NSCNG(K),NS3D(K)/DT)
2735
2736 ! PORTION OF RIMING LEFT FOR SNOW
2737              PSACWS(K) = PSACWS(K) - PGSACW(K)
2738              END IF
2739            END IF
2740
2741 ! CONVERSION OF RIMED RAINWATER ONTO SNOW CONVERTED TO GRAUPEL
2742
2743            IF (PRACS(K).GT.0.) THEN
2744 ! ONLY ALLOW CONVERSION IF QNI > 0.1 AND QR > 0.1 G/KG FOLLOWING RUTLEDGE AND HOBBS (1984)
2745               IF (QNI3D(K).GE.0.1E-3.AND.QR3D(K).GE.0.1E-3) THEN
2746 ! PORTION OF COLLECTED RAINWATER CONVERTED TO GRAUPEL (REISNER ET AL. 1998)
2747               DUM = CONS18*(4./LAMS(K))**3*(4./LAMS(K))**3 &
2748                    /(CONS18*(4./LAMS(K))**3*(4./LAMS(K))**3+ &
2749                    CONS19*(4./LAMR(K))**3*(4./LAMR(K))**3)
2750               DUM=MIN(DUM,1.)
2751               DUM=MAX(DUM,0.)
2752               PGRACS(K) = (1.-DUM)*PRACS(K)
2753             NGRACS(K) = (1.-DUM)*NPRACS(K)
2754 ! LIMIT MAX NUMBER CONVERTED TO MIN OF EITHER RAIN OR SNOW NUMBER CONCENTRATION
2755             NGRACS(K) = MIN(NGRACS(K),NR3D(K)/DT)
2756             NGRACS(K) = MIN(NGRACS(K),NS3D(K)/DT)
2757
2758 ! AMOUNT LEFT FOR SNOW PRODUCTION
2759             PRACS(K) = PRACS(K) - PGRACS(K)
2760             NPRACS(K) = NPRACS(K) - NGRACS(K)
2761 ! CONVERSION TO GRAUPEL DUE TO COLLECTION OF SNOW BY RAIN
2762             PSACR(K)=PSACR(K)*(1.-DUM)
2763             END IF
2764            END IF
2765 !           END IF
2766
2767 !.......................................................................
2768 ! FREEZING OF RAIN DROPS
2769 ! FREEZING ALLOWED BELOW -4 C
2770
2771          IF (T3D(K).LT.269.15.AND.QR3D(K).GE.QSMALL) THEN
2772
2773 ! IMMERSION FREEZING (BIGG 1953)
2774 !            MNUCCR(K) = CONS20*NR3D(K)*EXP(AIMM*(273.15-T3D(K)))/LAMR(K)**3 &
2775 !                 /LAMR(K)**3
2776
2777 !            NNUCCR(K) = PI*NR3D(K)*BIMM*EXP(AIMM*(273.15-T3D(K)))/LAMR(K)**3
2778
2779 ! hm fix 7/15/13 for consistency w/ original formula
2780             MNUCCR(K) = CONS20*NR3D(K)*(EXP(AIMM*(273.15-T3D(K)))-1.)/LAMR(K)**3 &
2781                  /LAMR(K)**3
2782
2783             NNUCCR(K) = PI*NR3D(K)*BIMM*(EXP(AIMM*(273.15-T3D(K)))-1.)/LAMR(K)**3
2784
2785 ! PREVENT DIVERGENCE BETWEEN MIXING RATIO AND NUMBER CONC
2786             NNUCCR(K) = MIN(NNUCCR(K),NR3D(K)/DT)
2787
2788          END IF
2789
2790 !.......................................................................
2791 ! ACCRETION OF CLOUD LIQUID WATER BY RAIN
2792 ! CONTINUOUS COLLECTION EQUATION WITH
2793 ! GRAVITATIONAL COLLECTION KERNEL, DROPLET FALL SPEED NEGLECTED
2794
2795          IF (QR3D(K).GE.1.E-8 .AND. QC3D(K).GE.1.E-8) THEN
2796
2797 ! 12/13/06 HM ADD, REPLACE WITH NEWER FORMULA FROM
2798 ! KHAIROUTDINOV AND KOGAN 2000, MWR
2799
2800            DUM=(QC3D(K)*QR3D(K))
2801            PRA(K) = 67.*(DUM)**1.15
2802            NPRA(K) = PRA(K)/(QC3D(K)/NC3D(K))
2803
2804          END IF
2805 !.......................................................................
2806 ! SELF-COLLECTION OF RAIN DROPS
2807 ! FROM BEHENG(1994)
2808 ! FROM NUMERICAL SIMULATION OF THE STOCHASTIC COLLECTION EQUATION
2809 ! AS DESCRINED ABOVE FOR AUTOCONVERSION
2810
2811          IF (QR3D(K).GE.1.E-8) THEN
2812 ! include breakup add 10/09/09
2813             dum1=300.e-6
2814             if (1./lamr(k).lt.dum1) then
2815             dum=1.
2816             else if (1./lamr(k).ge.dum1) then
2817             dum=2.-exp(2300.*(1./lamr(k)-dum1))
2818             end if
2819 !            NRAGG(K) = -8.*NR3D(K)*QR3D(K)*RHO(K)
2820             NRAGG(K) = -5.78*dum*NR3D(K)*QR3D(K)*RHO(K)
2821          END IF
2822
2823 !.......................................................................
2824 ! AUTOCONVERSION OF CLOUD ICE TO SNOW
2825 ! FOLLOWING HARRINGTON ET AL. (1995) WITH MODIFICATION
2826 ! HERE IT IS ASSUMED THAT AUTOCONVERSION CAN ONLY OCCUR WHEN THE
2827 ! ICE IS GROWING, I.E. IN CONDITIONS OF ICE SUPERSATURATION
2828
2829          IF (QI3D(K).GE.1.E-8 .AND.QVQVSI(K).GE.1.) THEN
2830
2831 !           COFFI = 2./LAMI(K)
2832 !           IF (COFFI.GE.DCS) THEN
2833               NPRCI(K) = CONS21*(QV3D(K)-QVI(K))*RHO(K)                         &
2834                 *N0I(K)*EXP(-LAMI(K)*DCS)*DV(K)/ABI(K)
2835               PRCI(K) = CONS22*NPRCI(K)
2836               NPRCI(K) = MIN(NPRCI(K),NI3D(K)/DT)
2837
2838 !           END IF
2839          END IF
2840
2841 !.......................................................................
2842 ! ACCRETION OF CLOUD ICE BY SNOW
2843 ! FOR THIS CALCULATION, IT IS ASSUMED THAT THE VS >> VI
2844 ! AND DS >> DI FOR CONTINUOUS COLLECTION
2845
2846          IF (QNI3D(K).GE.1.E-8 .AND. QI3D(K).GE.QSMALL) THEN
2847             PRAI(K) = CONS23*ASN(K)*QI3D(K)*RHO(K)*N0S(K)/     &
2848                      LAMS(K)**(BS+3.)
2849             NPRAI(K) = CONS23*ASN(K)*NI3D(K)*                                       &
2850                   RHO(K)*N0S(K)/                                 &
2851                   LAMS(K)**(BS+3.)
2852             NPRAI(K)=MIN(NPRAI(K),NI3D(K)/DT)
2853          END IF
2854
2855 !.......................................................................
2856 ! HM, ADD 12/13/06, COLLISION OF RAIN AND ICE TO PRODUCE SNOW OR GRAUPEL
2857 ! FOLLOWS REISNER ET AL. 1998
2858 ! ASSUMED FALLSPEED AND SIZE OF ICE CRYSTAL << THAN FOR RAIN
2859
2860          IF (QR3D(K).GE.1.E-8.AND.QI3D(K).GE.1.E-8.AND.T3D(K).LE.273.15) THEN
2861
2862 ! ALLOW GRAUPEL FORMATION FROM RAIN-ICE COLLISIONS ONLY IF RAIN MIXING RATIO > 0.1 G/KG,
2863 ! OTHERWISE ADD TO SNOW
2864
2865             IF (QR3D(K).GE.0.1E-3) THEN
2866             NIACR(K)=CONS24*NI3D(K)*N0RR(K)*ARN(K) &
2867                 /LAMR(K)**(BR+3.)*RHO(K)
2868             PIACR(K)=CONS25*NI3D(K)*N0RR(K)*ARN(K) &
2869                 /LAMR(K)**(BR+3.)/LAMR(K)**3*RHO(K)
2870             PRACI(K)=CONS24*QI3D(K)*N0RR(K)*ARN(K)/ &
2871                 LAMR(K)**(BR+3.)*RHO(K)
2872             NIACR(K)=MIN(NIACR(K),NR3D(K)/DT)
2873             NIACR(K)=MIN(NIACR(K),NI3D(K)/DT)
2874             ELSE
2875             NIACRS(K)=CONS24*NI3D(K)*N0RR(K)*ARN(K) &
2876                 /LAMR(K)**(BR+3.)*RHO(K)
2877             PIACRS(K)=CONS25*NI3D(K)*N0RR(K)*ARN(K) &
2878                 /LAMR(K)**(BR+3.)/LAMR(K)**3*RHO(K)
2879             PRACIS(K)=CONS24*QI3D(K)*N0RR(K)*ARN(K)/ &
2880                 LAMR(K)**(BR+3.)*RHO(K)
2881             NIACRS(K)=MIN(NIACRS(K),NR3D(K)/DT)
2882             NIACRS(K)=MIN(NIACRS(K),NI3D(K)/DT)
2883             END IF
2884          END IF
2885
2886 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
2887 ! NUCLEATION OF CLOUD ICE FROM HOMOGENEOUS AND HETEROGENEOUS FREEZING ON AEROSOL
2888
2889          IF (INUC.EQ.0) THEN
2890
2891 ! add threshold according to Greg Thomspon
2892
2893          if ((QVQVS(K).GE.0.999.and.T3D(K).le.265.15).or. &
2894               QVQVSI(K).ge.1.08) then
2895
2896 ! hm, modify dec. 5, 2006, replace with cooper curve
2897       kc2 = 0.005*exp(0.304*(273.15-T3D(K)))*1000. ! convert from L-1 to m-3
2898 ! limit to 500 L-1
2899       kc2 = min(kc2,500.e3)
2900       kc2=MAX(kc2/rho(k),0.)  ! convert to kg-1
2901
2902           IF (KC2.GT.NI3D(K)+NS3D(K)+NG3D(K)) THEN
2903              NNUCCD(K) = (KC2-NI3D(K)-NS3D(K)-NG3D(K))/DT
2904              MNUCCD(K) = NNUCCD(K)*MI0
2905           END IF
2906
2907           END IF
2908
2909           ELSE IF (INUC.EQ.1) THEN
2910
2911           IF (T3D(K).LT.273.15.AND.QVQVSI(K).GT.1.) THEN
2912
2913              KC2 = 0.16*1000./RHO(K)  ! CONVERT FROM L-1 TO KG-1
2914           IF (KC2.GT.NI3D(K)+NS3D(K)+NG3D(K)) THEN
2915              NNUCCD(K) = (KC2-NI3D(K)-NS3D(K)-NG3D(K))/DT
2916              MNUCCD(K) = NNUCCD(K)*MI0
2917           END IF
2918           END IF
2919
2920          END IF
2921
2922 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
2923
2924  101      CONTINUE
2925
2926 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
2927 ! CALCULATE EVAP/SUB/DEP TERMS FOR QI,QNI,QR
2928
2929 ! NO VENTILATION FOR CLOUD ICE
2930
2931         IF (QI3D(K).GE.QSMALL) THEN
2932
2933          EPSI = 2.*PI*N0I(K)*RHO(K)*DV(K)/(LAMI(K)*LAMI(K))
2934
2935       ELSE
2936          EPSI = 0.
2937       END IF
2938
2939       IF (QNI3D(K).GE.QSMALL) THEN
2940         EPSS = 2.*PI*N0S(K)*RHO(K)*DV(K)*                            &
2941                    (F1S/(LAMS(K)*LAMS(K))+                       &
2942                     F2S*(ASN(K)*RHO(K)/MU(K))**0.5*                      &
2943                     SC(K)**(1./3.)*CONS10/                   &
2944                (LAMS(K)**CONS35))
2945       ELSE
2946       EPSS = 0.
2947       END IF
2948
2949       IF (QG3D(K).GE.QSMALL) THEN
2950         EPSG = 2.*PI*N0G(K)*RHO(K)*DV(K)*                                &
2951                    (F1S/(LAMG(K)*LAMG(K))+                               &
2952                     F2S*(AGN(K)*RHO(K)/MU(K))**0.5*                      &
2953                     SC(K)**(1./3.)*CONS11/                   &
2954                (LAMG(K)**CONS36))
2955
2956
2957       ELSE
2958       EPSG = 0.
2959       END IF
2960
2961       IF (QR3D(K).GE.QSMALL) THEN
2962         EPSR = 2.*PI*N0RR(K)*RHO(K)*DV(K)*                           &
2963                    (F1R/(LAMR(K)*LAMR(K))+                       &
2964                     F2R*(ARN(K)*RHO(K)/MU(K))**0.5*                      &
2965                     SC(K)**(1./3.)*CONS9/                   &
2966                 (LAMR(K)**CONS34))
2967       ELSE
2968       EPSR = 0.
2969       END IF
2970
2971 ! ONLY INCLUDE REGION OF ICE SIZE DIST < DCS
2972 ! DUM IS FRACTION OF D*N(D) < DCS
2973
2974 ! LOGIC BELOW FOLLOWS THAT OF HARRINGTON ET AL. 1995 (JAS)
2975               IF (QI3D(K).GE.QSMALL) THEN
2976               DUM=(1.-EXP(-LAMI(K)*DCS)*(1.+LAMI(K)*DCS))
2977               PRD(K) = EPSI*(QV3D(K)-QVI(K))/ABI(K)*DUM
2978               ELSE
2979               DUM=0.
2980               END IF
2981 ! ADD DEPOSITION IN TAIL OF ICE SIZE DIST TO SNOW IF SNOW IS PRESENT
2982               IF (QNI3D(K).GE.QSMALL) THEN
2983               PRDS(K) = EPSS*(QV3D(K)-QVI(K))/ABI(K)+ &
2984                 EPSI*(QV3D(K)-QVI(K))/ABI(K)*(1.-DUM)
2985 ! OTHERWISE ADD TO CLOUD ICE
2986               ELSE
2987               PRD(K) = PRD(K)+EPSI*(QV3D(K)-QVI(K))/ABI(K)*(1.-DUM)
2988               END IF
2989 ! VAPOR DPEOSITION ON GRAUPEL
2990               PRDG(K) = EPSG*(QV3D(K)-QVI(K))/ABI(K)
2991
2992 ! NO CONDENSATION ONTO RAIN, ONLY EVAP
2993
2994            IF (QV3D(K).LT.QVS(K)) THEN
2995               PRE(K) = EPSR*(QV3D(K)-QVS(K))/AB(K)
2996               PRE(K) = MIN(PRE(K),0.)
2997            ELSE
2998               PRE(K) = 0.
2999            END IF
3000
3001 ! MAKE SURE NOT PUSHED INTO ICE SUPERSAT/SUBSAT
3002 ! FORMULA FROM REISNER 2 SCHEME
3003
3004            DUM = (QV3D(K)-QVI(K))/DT
3005
3006            FUDGEF = 0.9999
3007            SUM_DEP = PRD(K)+PRDS(K)+MNUCCD(K)+PRDG(K)
3008
3009            IF( (DUM.GT.0. .AND. SUM_DEP.GT.DUM*FUDGEF) .OR.                      &
3010                (DUM.LT.0. .AND. SUM_DEP.LT.DUM*FUDGEF) ) THEN
3011                MNUCCD(K) = FUDGEF*MNUCCD(K)*DUM/SUM_DEP
3012                PRD(K) = FUDGEF*PRD(K)*DUM/SUM_DEP
3013                PRDS(K) = FUDGEF*PRDS(K)*DUM/SUM_DEP
3014                PRDG(K) = FUDGEF*PRDG(K)*DUM/SUM_DEP
3015            ENDIF
3016
3017 ! IF CLOUD ICE/SNOW/GRAUPEL VAP DEPOSITION IS NEG, THEN ASSIGN TO SUBLIMATION PROCESSES
3018
3019            IF (PRD(K).LT.0.) THEN
3020               EPRD(K)=PRD(K)
3021               PRD(K)=0.
3022            END IF
3023            IF (PRDS(K).LT.0.) THEN
3024               EPRDS(K)=PRDS(K)
3025               PRDS(K)=0.
3026            END IF
3027            IF (PRDG(K).LT.0.) THEN
3028               EPRDG(K)=PRDG(K)
3029               PRDG(K)=0.
3030            END IF
3031
3032 !.......................................................................
3033 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
3034
3035 ! CONSERVATION OF WATER
3036 ! THIS IS ADOPTED LOOSELY FROM MM5 RESINER CODE. HOWEVER, HERE WE
3037 ! ONLY ADJUST PROCESSES THAT ARE NEGATIVE, RATHER THAN ALL PROCESSES.
3038
3039 ! IF MIXING RATIOS LESS THAN QSMALL, THEN NO DEPLETION OF WATER
3040 ! THROUGH MICROPHYSICAL PROCESSES, SKIP CONSERVATION
3041
3042 ! NOTE: CONSERVATION CHECK NOT APPLIED TO NUMBER CONCENTRATION SPECIES. ADDITIONAL CATCH
3043 ! BELOW WILL PREVENT NEGATIVE NUMBER CONCENTRATION
3044 ! FOR EACH MICROPHYSICAL PROCESS WHICH PROVIDES A SOURCE FOR NUMBER, THERE IS A CHECK
3045 ! TO MAKE SURE THAT CAN'T EXCEED TOTAL NUMBER OF DEPLETED SPECIES WITH THE TIME
3046 ! STEP
3047
3048 !****SENSITIVITY - NO ICE
3049
3050       IF (ILIQ.EQ.1) THEN
3051       MNUCCC(K)=0.
3052       NNUCCC(K)=0.
3053       MNUCCR(K)=0.
3054       NNUCCR(K)=0.
3055       MNUCCD(K)=0.
3056       NNUCCD(K)=0.
3057       END IF
3058
3059 ! ****SENSITIVITY - NO GRAUPEL
3060       IF (IGRAUP.EQ.1) THEN
3061             PRACG(K) = 0.
3062             PSACR(K) = 0.
3063             PSACWG(K) = 0.
3064             PRDG(K) = 0.
3065             EPRDG(K) = 0.
3066             EVPMG(K) = 0.
3067             PGMLT(K) = 0.
3068             NPRACG(K) = 0.
3069             NPSACWG(K) = 0.
3070             NSCNG(K) = 0.
3071             NGRACS(K) = 0.
3072             NSUBG(K) = 0.
3073             NGMLTG(K) = 0.
3074             NGMLTR(K) = 0.
3075 ! fix 053011
3076             PIACRS(K)=PIACRS(K)+PIACR(K)
3077             PIACR(K) = 0.
3078 ! fix 070713
3079             PRACIS(K)=PRACIS(K)+PRACI(K)
3080             PRACI(K) = 0.
3081             PSACWS(K)=PSACWS(K)+PGSACW(K)
3082             PGSACW(K) = 0.
3083             PRACS(K)=PRACS(K)+PGRACS(K)
3084             PGRACS(K) = 0.
3085        END IF
3086
3087 ! CONSERVATION OF QC
3088
3089       DUM = (PRC(K)+PRA(K)+MNUCCC(K)+PSACWS(K)+PSACWI(K)+QMULTS(K)+PSACWG(K)+PGSACW(K)+QMULTG(K))*DT
3090
3091       IF (DUM.GT.QC3D(K).AND.QC3D(K).GE.QSMALL) THEN
3092         RATIO = QC3D(K)/DUM
3093
3094         PRC(K) = PRC(K)*RATIO
3095         PRA(K) = PRA(K)*RATIO
3096         MNUCCC(K) = MNUCCC(K)*RATIO
3097         PSACWS(K) = PSACWS(K)*RATIO
3098         PSACWI(K) = PSACWI(K)*RATIO
3099         QMULTS(K) = QMULTS(K)*RATIO
3100         QMULTG(K) = QMULTG(K)*RATIO
3101         PSACWG(K) = PSACWG(K)*RATIO
3102         PGSACW(K) = PGSACW(K)*RATIO
3103         END IF
3104
3105 ! CONSERVATION OF QI
3106
3107       DUM = (-PRD(K)-MNUCCC(K)+PRCI(K)+PRAI(K)-QMULTS(K)-QMULTG(K)-QMULTR(K)-QMULTRG(K) &
3108                 -MNUCCD(K)+PRACI(K)+PRACIS(K)-EPRD(K)-PSACWI(K))*DT
3109
3110       IF (DUM.GT.QI3D(K).AND.QI3D(K).GE.QSMALL) THEN
3111
3112         RATIO = (QI3D(K)/DT+PRD(K)+MNUCCC(K)+QMULTS(K)+QMULTG(K)+QMULTR(K)+QMULTRG(K)+ &
3113                      MNUCCD(K)+PSACWI(K))/ &
3114                       (PRCI(K)+PRAI(K)+PRACI(K)+PRACIS(K)-EPRD(K))
3115
3116         PRCI(K) = PRCI(K)*RATIO
3117         PRAI(K) = PRAI(K)*RATIO
3118         PRACI(K) = PRACI(K)*RATIO
3119         PRACIS(K) = PRACIS(K)*RATIO
3120         EPRD(K) = EPRD(K)*RATIO
3121
3122         END IF
3123
3124 ! CONSERVATION OF QR
3125
3126       DUM=((PRACS(K)-PRE(K))+(QMULTR(K)+QMULTRG(K)-PRC(K))+(MNUCCR(K)-PRA(K))+ &
3127              PIACR(K)+PIACRS(K)+PGRACS(K)+PRACG(K))*DT
3128
3129       IF (DUM.GT.QR3D(K).AND.QR3D(K).GE.QSMALL) THEN
3130
3131         RATIO = (QR3D(K)/DT+PRC(K)+PRA(K))/ &
3132              (-PRE(K)+QMULTR(K)+QMULTRG(K)+PRACS(K)+MNUCCR(K)+PIACR(K)+PIACRS(K)+PGRACS(K)+PRACG(K))
3133
3134         PRE(K) = PRE(K)*RATIO
3135         PRACS(K) = PRACS(K)*RATIO
3136         QMULTR(K) = QMULTR(K)*RATIO
3137         QMULTRG(K) = QMULTRG(K)*RATIO
3138         MNUCCR(K) = MNUCCR(K)*RATIO
3139         PIACR(K) = PIACR(K)*RATIO
3140         PIACRS(K) = PIACRS(K)*RATIO
3141         PGRACS(K) = PGRACS(K)*RATIO
3142         PRACG(K) = PRACG(K)*RATIO
3143
3144         END IF
3145
3146 ! CONSERVATION OF QNI
3147 ! CONSERVATION FOR GRAUPEL SCHEME
3148
3149         IF (IGRAUP.EQ.0) THEN
3150
3151       DUM = (-PRDS(K)-PSACWS(K)-PRAI(K)-PRCI(K)-PRACS(K)-EPRDS(K)+PSACR(K)-PIACRS(K)-PRACIS(K))*DT
3152
3153       IF (DUM.GT.QNI3D(K).AND.QNI3D(K).GE.QSMALL) THEN
3154
3155         RATIO = (QNI3D(K)/DT+PRDS(K)+PSACWS(K)+PRAI(K)+PRCI(K)+PRACS(K)+PIACRS(K)+PRACIS(K))/(-EPRDS(K)+PSACR(K))
3156
3157        EPRDS(K) = EPRDS(K)*RATIO
3158        PSACR(K) = PSACR(K)*RATIO
3159
3160        END IF
3161
3162 ! FOR NO GRAUPEL, NEED TO INCLUDE FREEZING OF RAIN FOR SNOW
3163        ELSE IF (IGRAUP.EQ.1) THEN
3164
3165       DUM = (-PRDS(K)-PSACWS(K)-PRAI(K)-PRCI(K)-PRACS(K)-EPRDS(K)+PSACR(K)-PIACRS(K)-PRACIS(K)-MNUCCR(K))*DT
3166
3167       IF (DUM.GT.QNI3D(K).AND.QNI3D(K).GE.QSMALL) THEN
3168
3169        RATIO = (QNI3D(K)/DT+PRDS(K)+PSACWS(K)+PRAI(K)+PRCI(K)+PRACS(K)+PIACRS(K)+PRACIS(K)+MNUCCR(K))/(-EPRDS(K)+PSACR(K))
3170
3171        EPRDS(K) = EPRDS(K)*RATIO
3172        PSACR(K) = PSACR(K)*RATIO
3173
3174        END IF
3175
3176        END IF
3177
3178 ! CONSERVATION OF QG
3179
3180       DUM = (-PSACWG(K)-PRACG(K)-PGSACW(K)-PGRACS(K)-PRDG(K)-MNUCCR(K)-EPRDG(K)-PIACR(K)-PRACI(K)-PSACR(K))*DT
3181
3182       IF (DUM.GT.QG3D(K).AND.QG3D(K).GE.QSMALL) THEN
3183
3184         RATIO = (QG3D(K)/DT+PSACWG(K)+PRACG(K)+PGSACW(K)+PGRACS(K)+PRDG(K)+MNUCCR(K)+PSACR(K)+&
3185                   PIACR(K)+PRACI(K))/(-EPRDG(K))
3186
3187        EPRDG(K) = EPRDG(K)*RATIO
3188
3189       END IF
3190
3191 ! TENDENCIES
3192
3193       QV3DTEN(K) = QV3DTEN(K)+(-PRE(K)-PRD(K)-PRDS(K)-MNUCCD(K)-EPRD(K)-EPRDS(K)-PRDG(K)-EPRDG(K))
3194
3195 ! BUG FIX HM, 3/1/11, INCLUDE PIACR AND PIACRS
3196       T3DTEN(K) = T3DTEN(K)+(PRE(K)                                 &
3197                *XXLV(K)+(PRD(K)+PRDS(K)+                            &
3198                 MNUCCD(K)+EPRD(K)+EPRDS(K)+PRDG(K)+EPRDG(K))*XXLS(K)+         &
3199                (PSACWS(K)+PSACWI(K)+MNUCCC(K)+MNUCCR(K)+                      &
3200                 QMULTS(K)+QMULTG(K)+QMULTR(K)+QMULTRG(K)+PRACS(K) &
3201                 +PSACWG(K)+PRACG(K)+PGSACW(K)+PGRACS(K)+PIACR(K)+PIACRS(K))*XLF(K))/CPM(K)
3202
3203       QC3DTEN(K) = QC3DTEN(K)+                                      &
3204                  (-PRA(K)-PRC(K)-MNUCCC(K)+PCC(K)-                  &
3205                   PSACWS(K)-PSACWI(K)-QMULTS(K)-QMULTG(K)-PSACWG(K)-PGSACW(K))
3206       QI3DTEN(K) = QI3DTEN(K)+                                      &
3207          (PRD(K)+EPRD(K)+PSACWI(K)+MNUCCC(K)-PRCI(K)-                                 &
3208                   PRAI(K)+QMULTS(K)+QMULTG(K)+QMULTR(K)+QMULTRG(K)+MNUCCD(K)-PRACI(K)-PRACIS(K))
3209       QR3DTEN(K) = QR3DTEN(K)+                                      &
3210                  (PRE(K)+PRA(K)+PRC(K)-PRACS(K)-MNUCCR(K)-QMULTR(K)-QMULTRG(K) &
3211              -PIACR(K)-PIACRS(K)-PRACG(K)-PGRACS(K))
3212
3213       IF (IGRAUP.EQ.0) THEN
3214
3215       QNI3DTEN(K) = QNI3DTEN(K)+                                    &
3216            (PRAI(K)+PSACWS(K)+PRDS(K)+PRACS(K)+PRCI(K)+EPRDS(K)-PSACR(K)+PIACRS(K)+PRACIS(K))
3217       NS3DTEN(K) = NS3DTEN(K)+(NSAGG(K)+NPRCI(K)-NSCNG(K)-NGRACS(K)+NIACRS(K))
3218       QG3DTEN(K) = QG3DTEN(K)+(PRACG(K)+PSACWG(K)+PGSACW(K)+PGRACS(K)+ &
3219                     PRDG(K)+EPRDG(K)+MNUCCR(K)+PIACR(K)+PRACI(K)+PSACR(K))
3220       NG3DTEN(K) = NG3DTEN(K)+(NSCNG(K)+NGRACS(K)+NNUCCR(K)+NIACR(K))
3221
3222 ! FOR NO GRAUPEL, NEED TO INCLUDE FREEZING OF RAIN FOR SNOW
3223       ELSE IF (IGRAUP.EQ.1) THEN
3224
3225       QNI3DTEN(K) = QNI3DTEN(K)+                                    &
3226            (PRAI(K)+PSACWS(K)+PRDS(K)+PRACS(K)+PRCI(K)+EPRDS(K)-PSACR(K)+PIACRS(K)+PRACIS(K)+MNUCCR(K))
3227       NS3DTEN(K) = NS3DTEN(K)+(NSAGG(K)+NPRCI(K)-NSCNG(K)-NGRACS(K)+NIACRS(K)+NNUCCR(K))
3228
3229       END IF
3230
3231       NC3DTEN(K) = NC3DTEN(K)+(-NNUCCC(K)-NPSACWS(K)                &
3232             -NPRA(K)-NPRC(K)-NPSACWI(K)-NPSACWG(K))
3233
3234       NI3DTEN(K) = NI3DTEN(K)+                                      &
3235        (NNUCCC(K)-NPRCI(K)-NPRAI(K)+NMULTS(K)+NMULTG(K)+NMULTR(K)+NMULTRG(K)+ &
3236                NNUCCD(K)-NIACR(K)-NIACRS(K))
3237
3238       NR3DTEN(K) = NR3DTEN(K)+(NPRC1(K)-NPRACS(K)-NNUCCR(K)      &
3239                    +NRAGG(K)-NIACR(K)-NIACRS(K)-NPRACG(K)-NGRACS(K))
3240
3241 ! HM ADD, WRF-CHEM, ADD TENDENCIES FOR C2PREC
3242
3243         C2PREC(K) = PRA(K)+PRC(K)+PSACWS(K)+QMULTS(K)+QMULTG(K)+PSACWG(K)+ &
3244        PGSACW(K)+MNUCCC(K)+PSACWI(K)
3245 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
3246 ! NOW CALCULATE SATURATION ADJUSTMENT TO CONDENSE EXTRA VAPOR ABOVE
3247 ! WATER SATURATION
3248
3249       DUMT = T3D(K)+DT*T3DTEN(K)
3250       DUMQV = QV3D(K)+DT*QV3DTEN(K)
3251 ! hm, add fix for low pressure, 5/12/10
3252       dum=min(0.99*pres(k),POLYSVP(DUMT,0))
3253       DUMQSS = EP_2*dum/(PRES(K)-dum)
3254       DUMQC = QC3D(K)+DT*QC3DTEN(K)
3255       DUMQC = MAX(DUMQC,0.)
3256
3257 ! SATURATION ADJUSTMENT FOR LIQUID
3258
3259       DUMS = DUMQV-DUMQSS
3260       PCC(K) = DUMS/(1.+XXLV(K)**2*DUMQSS/(CPM(K)*RV*DUMT**2))/DT
3261       IF (PCC(K)*DT+DUMQC.LT.0.) THEN
3262            PCC(K) = -DUMQC/DT
3263       END IF
3264
3265       QV3DTEN(K) = QV3DTEN(K)-PCC(K)
3266       T3DTEN(K) = T3DTEN(K)+PCC(K)*XXLV(K)/CPM(K)
3267       QC3DTEN(K) = QC3DTEN(K)+PCC(K)
3268
3269 !.......................................................................
3270 ! ACTIVATION OF CLOUD DROPLETS
3271 ! ACTIVATION OF DROPLET CURRENTLY NOT CALCULATED
3272 ! DROPLET CONCENTRATION IS SPECIFIED !!!!!
3273
3274 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
3275 ! SUBLIMATE, MELT, OR EVAPORATE NUMBER CONCENTRATION
3276 ! THIS FORMULATION ASSUMES 1:1 RATIO BETWEEN MASS LOSS AND
3277 ! LOSS OF NUMBER CONCENTRATION
3278
3279 !     IF (PCC(K).LT.0.) THEN
3280 !        DUM = PCC(K)*DT/QC3D(K)
3281 !           DUM = MAX(-1.,DUM)
3282 !        NSUBC(K) = DUM*NC3D(K)/DT
3283 !     END IF
3284
3285       IF (EPRD(K).LT.0.) THEN
3286          DUM = EPRD(K)*DT/QI3D(K)
3287             DUM = MAX(-1.,DUM)
3288          NSUBI(K) = DUM*NI3D(K)/DT
3289       END IF
3290       IF (EPRDS(K).LT.0.) THEN
3291          DUM = EPRDS(K)*DT/QNI3D(K)
3292            DUM = MAX(-1.,DUM)
3293          NSUBS(K) = DUM*NS3D(K)/DT
3294       END IF
3295       IF (PRE(K).LT.0.) THEN
3296          DUM = PRE(K)*DT/QR3D(K)
3297            DUM = MAX(-1.,DUM)
3298          NSUBR(K) = DUM*NR3D(K)/DT
3299       END IF
3300       IF (EPRDG(K).LT.0.) THEN
3301          DUM = EPRDG(K)*DT/QG3D(K)
3302            DUM = MAX(-1.,DUM)
3303          NSUBG(K) = DUM*NG3D(K)/DT
3304       END IF
3305
3306 !        nsubr(k)=0.
3307 !        nsubs(k)=0.
3308 !        nsubg(k)=0.
3309
3310 ! UPDATE TENDENCIES
3311
3312 !        NC3DTEN(K) = NC3DTEN(K)+NSUBC(K)
3313          NI3DTEN(K) = NI3DTEN(K)+NSUBI(K)
3314          NS3DTEN(K) = NS3DTEN(K)+NSUBS(K)
3315          NG3DTEN(K) = NG3DTEN(K)+NSUBG(K)
3316          NR3DTEN(K) = NR3DTEN(K)+NSUBR(K)
3317
3318 #if (WRF_CHEM == 1)
3319       IF( has_wetscav ) THEN
3320          evapprod(k) = - PRE(K) - EPRDS(K) - EPRDG(K)
3321          rainprod(k) = PRA(K) + PRC(K) + PSACWS(K) + PSACWG(K) + PGSACW(K) &
3322                        + PRAI(K) + PRCI(K) + PRACI(K) + PRACIS(K)  &
3323                        + PRDS(K) + PRDG(K)
3324       ENDIF
3325 #endif
3326
3327          END IF !!!!!! TEMPERATURE
3328
3329 ! SWITCH LTRUE TO 1, SINCE HYDROMETEORS ARE PRESENT
3330          LTRUE = 1
3331
3332  200     CONTINUE
3333
3334         END DO
3335
3336 ! INITIALIZE PRECIP AND SNOW RATES
3337       PRECRT = 0.
3338       SNOWRT = 0.
3339 ! hm added 7/13/13
3340       SNOWPRT = 0.
3341       GRPLPRT = 0.
3342
3343 ! IF THERE ARE NO HYDROMETEORS, THEN SKIP TO END OF SUBROUTINE
3344
3345         IF (LTRUE.EQ.0) GOTO 400
3346
3347 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
3348 !.......................................................................
3349 ! CALCULATE SEDIMENATION
3350 ! THE NUMERICS HERE FOLLOW FROM REISNER ET AL. (1998)
3351 ! FALLOUT TERMS ARE CALCULATED ON SPLIT TIME STEPS TO ENSURE NUMERICAL
3352 ! STABILITY, I.E. COURANT# < 1
3353
3354 !.......................................................................
3355
3356       NSTEP = 1
3357
3358       DO K = KTE,KTS,-1
3359
3360         DUMI(K) = QI3D(K)+QI3DTEN(K)*DT
3361         DUMQS(K) = QNI3D(K)+QNI3DTEN(K)*DT
3362         DUMR(K) = QR3D(K)+QR3DTEN(K)*DT
3363         DUMFNI(K) = NI3D(K)+NI3DTEN(K)*DT
3364         DUMFNS(K) = NS3D(K)+NS3DTEN(K)*DT
3365         DUMFNR(K) = NR3D(K)+NR3DTEN(K)*DT
3366         DUMC(K) = QC3D(K)+QC3DTEN(K)*DT
3367         DUMFNC(K) = NC3D(K)+NC3DTEN(K)*DT
3368         DUMG(K) = QG3D(K)+QG3DTEN(K)*DT
3369         DUMFNG(K) = NG3D(K)+NG3DTEN(K)*DT
3370
3371 ! SWITCH FOR CONSTANT DROPLET NUMBER
3372         IF (iinum.EQ.1) THEN
3373         DUMFNC(K) = NC3D(K)
3374         END IF
3375
3376 ! GET DUMMY LAMDA FOR SEDIMENTATION CALCULATIONS
3377
3378 ! MAKE SURE NUMBER CONCENTRATIONS ARE POSITIVE
3379       DUMFNI(K) = MAX(0.,DUMFNI(K))
3380       DUMFNS(K) = MAX(0.,DUMFNS(K))
3381       DUMFNC(K) = MAX(0.,DUMFNC(K))
3382       DUMFNR(K) = MAX(0.,DUMFNR(K))
3383       DUMFNG(K) = MAX(0.,DUMFNG(K))
3384
3385 !......................................................................
3386 ! CLOUD ICE
3387
3388       IF (DUMI(K).GE.QSMALL) THEN
3389         DLAMI = (CONS12*DUMFNI(K)/DUMI(K))**(1./DI)
3390         DLAMI=MAX(DLAMI,LAMMINI)
3391         DLAMI=MIN(DLAMI,LAMMAXI)
3392       END IF
3393 !......................................................................
3394 ! RAIN
3395
3396       IF (DUMR(K).GE.QSMALL) THEN
3397         DLAMR = (PI*RHOW*DUMFNR(K)/DUMR(K))**(1./3.)
3398         DLAMR=MAX(DLAMR,LAMMINR)
3399         DLAMR=MIN(DLAMR,LAMMAXR)
3400       END IF
3401 !......................................................................
3402 ! CLOUD DROPLETS
3403
3404       IF (DUMC(K).GE.QSMALL) THEN
3405          DUM = PRES(K)/(287.15*T3D(K))
3406          PGAM(K)=0.0005714*(NC3D(K)/1.E6*DUM)+0.2714
3407          PGAM(K)=1./(PGAM(K)**2)-1.
3408          PGAM(K)=MAX(PGAM(K),2.)
3409          PGAM(K)=MIN(PGAM(K),10.)
3410
3411         DLAMC = (CONS26*DUMFNC(K)*GAMMA(PGAM(K)+4.)/(DUMC(K)*GAMMA(PGAM(K)+1.)))**(1./3.)
3412         LAMMIN = (PGAM(K)+1.)/60.E-6
3413         LAMMAX = (PGAM(K)+1.)/1.E-6
3414         DLAMC=MAX(DLAMC,LAMMIN)
3415         DLAMC=MIN(DLAMC,LAMMAX)
3416       END IF
3417 !......................................................................
3418 ! SNOW
3419
3420       IF (DUMQS(K).GE.QSMALL) THEN
3421         DLAMS = (CONS1*DUMFNS(K)/ DUMQS(K))**(1./DS)
3422         DLAMS=MAX(DLAMS,LAMMINS)
3423         DLAMS=MIN(DLAMS,LAMMAXS)
3424       END IF
3425 !......................................................................
3426 ! GRAUPEL
3427
3428       IF (DUMG(K).GE.QSMALL) THEN
3429         DLAMG = (CONS2*DUMFNG(K)/ DUMG(K))**(1./DG)
3430         DLAMG=MAX(DLAMG,LAMMING)
3431         DLAMG=MIN(DLAMG,LAMMAXG)
3432       END IF
3433
3434 !......................................................................
3435 ! CALCULATE NUMBER-WEIGHTED AND MASS-WEIGHTED TERMINAL FALL SPEEDS
3436
3437 ! CLOUD WATER
3438
3439       IF (DUMC(K).GE.QSMALL) THEN
3440       UNC =  ACN(K)*GAMMA(1.+BC+PGAM(K))/ (DLAMC**BC*GAMMA(PGAM(K)+1.))
3441       UMC = ACN(K)*GAMMA(4.+BC+PGAM(K))/  (DLAMC**BC*GAMMA(PGAM(K)+4.))
3442       ELSE
3443       UMC = 0.
3444       UNC = 0.
3445       END IF
3446
3447       IF (DUMI(K).GE.QSMALL) THEN
3448       UNI =  AIN(K)*CONS27/DLAMI**BI
3449       UMI = AIN(K)*CONS28/(DLAMI**BI)
3450       ELSE
3451       UMI = 0.
3452       UNI = 0.
3453       END IF
3454
3455       IF (DUMR(K).GE.QSMALL) THEN
3456       UNR = ARN(K)*CONS6/DLAMR**BR
3457       UMR = ARN(K)*CONS4/(DLAMR**BR)
3458       ELSE
3459       UMR = 0.
3460       UNR = 0.
3461       END IF
3462
3463       IF (DUMQS(K).GE.QSMALL) THEN
3464       UMS = ASN(K)*CONS3/(DLAMS**BS)
3465       UNS = ASN(K)*CONS5/DLAMS**BS
3466       ELSE
3467       UMS = 0.
3468       UNS = 0.
3469       END IF
3470
3471       IF (DUMG(K).GE.QSMALL) THEN
3472       UMG = AGN(K)*CONS7/(DLAMG**BG)
3473       UNG = AGN(K)*CONS8/DLAMG**BG
3474       ELSE
3475       UMG = 0.
3476       UNG = 0.
3477       END IF
3478
3479 ! SET REALISTIC LIMITS ON FALLSPEED
3480
3481 ! bug fix, 10/08/09
3482         dum=(rhosu/rho(k))**0.54
3483         UMS=MIN(UMS,1.2*dum)
3484         UNS=MIN(UNS,1.2*dum)
3485 ! fix 053011
3486 ! fix for correction by AA 4/6/11
3487         UMI=MIN(UMI,1.2*(rhosu/rho(k))**0.35)
3488         UNI=MIN(UNI,1.2*(rhosu/rho(k))**0.35)
3489         UMR=MIN(UMR,9.1*dum)
3490         UNR=MIN(UNR,9.1*dum)
3491         UMG=MIN(UMG,20.*dum)
3492         UNG=MIN(UNG,20.*dum)
3493
3494       FR(K) = UMR
3495       FI(K) = UMI
3496       FNI(K) = UNI
3497       FS(K) = UMS
3498       FNS(K) = UNS
3499       FNR(K) = UNR
3500       FC(K) = UMC
3501       FNC(K) = UNC
3502       FG(K) = UMG
3503       FNG(K) = UNG
3504
3505 ! V3.3 MODIFY FALLSPEED BELOW LEVEL OF PRECIP
3506
3507         IF (K.LE.KTE-1) THEN
3508         IF (FR(K).LT.1.E-10) THEN
3509         FR(K)=FR(K+1)
3510         END IF
3511         IF (FI(K).LT.1.E-10) THEN
3512         FI(K)=FI(K+1)
3513         END IF
3514         IF (FNI(K).LT.1.E-10) THEN
3515         FNI(K)=FNI(K+1)
3516         END IF
3517         IF (FS(K).LT.1.E-10) THEN
3518         FS(K)=FS(K+1)
3519         END IF
3520         IF (FNS(K).LT.1.E-10) THEN
3521         FNS(K)=FNS(K+1)
3522         END IF
3523         IF (FNR(K).LT.1.E-10) THEN
3524         FNR(K)=FNR(K+1)
3525         END IF
3526         IF (FC(K).LT.1.E-10) THEN
3527         FC(K)=FC(K+1)
3528         END IF
3529         IF (FNC(K).LT.1.E-10) THEN
3530         FNC(K)=FNC(K+1)
3531         END IF
3532         IF (FG(K).LT.1.E-10) THEN
3533         FG(K)=FG(K+1)
3534         END IF
3535         IF (FNG(K).LT.1.E-10) THEN
3536         FNG(K)=FNG(K+1)
3537         END IF
3538         END IF ! K LE KTE-1
3539
3540 ! CALCULATE NUMBER OF SPLIT TIME STEPS
3541
3542       RGVM = MAX(FR(K),FI(K),FS(K),FC(K),FNI(K),FNR(K),FNS(K),FNC(K),FG(K),FNG(K))
3543 ! VVT CHANGED IFIX -> INT (GENERIC FUNCTION)
3544       NSTEP = MAX(INT(RGVM*DT/DZQ(K)+1.),NSTEP)
3545
3546 ! MULTIPLY VARIABLES BY RHO
3547       DUMR(k) = DUMR(k)*RHO(K)
3548       DUMI(k) = DUMI(k)*RHO(K)
3549       DUMFNI(k) = DUMFNI(K)*RHO(K)
3550       DUMQS(k) = DUMQS(K)*RHO(K)
3551       DUMFNS(k) = DUMFNS(K)*RHO(K)
3552       DUMFNR(k) = DUMFNR(K)*RHO(K)
3553       DUMC(k) = DUMC(K)*RHO(K)
3554       DUMFNC(k) = DUMFNC(K)*RHO(K)
3555       DUMG(k) = DUMG(K)*RHO(K)
3556       DUMFNG(k) = DUMFNG(K)*RHO(K)
3557
3558       END DO
3559
3560       DO N = 1,NSTEP
3561
3562       DO K = KTS,KTE
3563       FALOUTR(K) = FR(K)*DUMR(K)
3564       FALOUTI(K) = FI(K)*DUMI(K)
3565       FALOUTNI(K) = FNI(K)*DUMFNI(K)
3566       FALOUTS(K) = FS(K)*DUMQS(K)
3567       FALOUTNS(K) = FNS(K)*DUMFNS(K)
3568       FALOUTNR(K) = FNR(K)*DUMFNR(K)
3569       FALOUTC(K) = FC(K)*DUMC(K)
3570       FALOUTNC(K) = FNC(K)*DUMFNC(K)
3571       FALOUTG(K) = FG(K)*DUMG(K)
3572       FALOUTNG(K) = FNG(K)*DUMFNG(K)
3573       END DO
3574
3575 ! TOP OF MODEL
3576
3577       K = KTE
3578       FALTNDR = FALOUTR(K)/DZQ(k)
3579       FALTNDI = FALOUTI(K)/DZQ(k)
3580       FALTNDNI = FALOUTNI(K)/DZQ(k)
3581       FALTNDS = FALOUTS(K)/DZQ(k)
3582       FALTNDNS = FALOUTNS(K)/DZQ(k)
3583       FALTNDNR = FALOUTNR(K)/DZQ(k)
3584       FALTNDC = FALOUTC(K)/DZQ(k)
3585       FALTNDNC = FALOUTNC(K)/DZQ(k)
3586       FALTNDG = FALOUTG(K)/DZQ(k)
3587       FALTNDNG = FALOUTNG(K)/DZQ(k)
3588 ! ADD FALLOUT TERMS TO EULERIAN TENDENCIES
3589
3590       QRSTEN(K) = QRSTEN(K)-FALTNDR/NSTEP/RHO(k)
3591       QISTEN(K) = QISTEN(K)-FALTNDI/NSTEP/RHO(k)
3592       NI3DTEN(K) = NI3DTEN(K)-FALTNDNI/NSTEP/RHO(k)
3593       QNISTEN(K) = QNISTEN(K)-FALTNDS/NSTEP/RHO(k)
3594       NS3DTEN(K) = NS3DTEN(K)-FALTNDNS/NSTEP/RHO(k)
3595       NR3DTEN(K) = NR3DTEN(K)-FALTNDNR/NSTEP/RHO(k)
3596       QCSTEN(K) = QCSTEN(K)-FALTNDC/NSTEP/RHO(k)
3597       NC3DTEN(K) = NC3DTEN(K)-FALTNDNC/NSTEP/RHO(k)
3598       QGSTEN(K) = QGSTEN(K)-FALTNDG/NSTEP/RHO(k)
3599       NG3DTEN(K) = NG3DTEN(K)-FALTNDNG/NSTEP/RHO(k)
3600
3601       DUMR(K) = DUMR(K)-FALTNDR*DT/NSTEP
3602       DUMI(K) = DUMI(K)-FALTNDI*DT/NSTEP
3603       DUMFNI(K) = DUMFNI(K)-FALTNDNI*DT/NSTEP
3604       DUMQS(K) = DUMQS(K)-FALTNDS*DT/NSTEP
3605       DUMFNS(K) = DUMFNS(K)-FALTNDNS*DT/NSTEP
3606       DUMFNR(K) = DUMFNR(K)-FALTNDNR*DT/NSTEP
3607       DUMC(K) = DUMC(K)-FALTNDC*DT/NSTEP
3608       DUMFNC(K) = DUMFNC(K)-FALTNDNC*DT/NSTEP
3609       DUMG(K) = DUMG(K)-FALTNDG*DT/NSTEP
3610       DUMFNG(K) = DUMFNG(K)-FALTNDNG*DT/NSTEP
3611
3612       DO K = KTE-1,KTS,-1
3613       FALTNDR = (FALOUTR(K+1)-FALOUTR(K))/DZQ(K)
3614       FALTNDI = (FALOUTI(K+1)-FALOUTI(K))/DZQ(K)
3615       FALTNDNI = (FALOUTNI(K+1)-FALOUTNI(K))/DZQ(K)
3616       FALTNDS = (FALOUTS(K+1)-FALOUTS(K))/DZQ(K)
3617       FALTNDNS = (FALOUTNS(K+1)-FALOUTNS(K))/DZQ(K)
3618       FALTNDNR = (FALOUTNR(K+1)-FALOUTNR(K))/DZQ(K)
3619       FALTNDC = (FALOUTC(K+1)-FALOUTC(K))/DZQ(K)
3620       FALTNDNC = (FALOUTNC(K+1)-FALOUTNC(K))/DZQ(K)
3621       FALTNDG = (FALOUTG(K+1)-FALOUTG(K))/DZQ(K)
3622       FALTNDNG = (FALOUTNG(K+1)-FALOUTNG(K))/DZQ(K)
3623
3624 ! ADD FALLOUT TERMS TO EULERIAN TENDENCIES
3625
3626       QRSTEN(K) = QRSTEN(K)+FALTNDR/NSTEP/RHO(k)
3627       QISTEN(K) = QISTEN(K)+FALTNDI/NSTEP/RHO(k)
3628       NI3DTEN(K) = NI3DTEN(K)+FALTNDNI/NSTEP/RHO(k)
3629       QNISTEN(K) = QNISTEN(K)+FALTNDS/NSTEP/RHO(k)
3630       NS3DTEN(K) = NS3DTEN(K)+FALTNDNS/NSTEP/RHO(k)
3631       NR3DTEN(K) = NR3DTEN(K)+FALTNDNR/NSTEP/RHO(k)
3632       QCSTEN(K) = QCSTEN(K)+FALTNDC/NSTEP/RHO(k)
3633       NC3DTEN(K) = NC3DTEN(K)+FALTNDNC/NSTEP/RHO(k)
3634       QGSTEN(K) = QGSTEN(K)+FALTNDG/NSTEP/RHO(k)
3635       NG3DTEN(K) = NG3DTEN(K)+FALTNDNG/NSTEP/RHO(k)
3636
3637       DUMR(K) = DUMR(K)+FALTNDR*DT/NSTEP
3638       DUMI(K) = DUMI(K)+FALTNDI*DT/NSTEP
3639       DUMFNI(K) = DUMFNI(K)+FALTNDNI*DT/NSTEP
3640       DUMQS(K) = DUMQS(K)+FALTNDS*DT/NSTEP
3641       DUMFNS(K) = DUMFNS(K)+FALTNDNS*DT/NSTEP
3642       DUMFNR(K) = DUMFNR(K)+FALTNDNR*DT/NSTEP
3643       DUMC(K) = DUMC(K)+FALTNDC*DT/NSTEP
3644       DUMFNC(K) = DUMFNC(K)+FALTNDNC*DT/NSTEP
3645       DUMG(K) = DUMG(K)+FALTNDG*DT/NSTEP
3646       DUMFNG(K) = DUMFNG(K)+FALTNDNG*DT/NSTEP
3647
3648 ! FOR WRF-CHEM, NEED PRECIP RATES (UNITS OF KG/M^2/S)
3649           CSED(K)=CSED(K)+FALOUTC(K)/NSTEP
3650           ISED(K)=ISED(K)+FALOUTI(K)/NSTEP
3651           SSED(K)=SSED(K)+FALOUTS(K)/NSTEP
3652           GSED(K)=GSED(K)+FALOUTG(K)/NSTEP
3653           RSED(K)=RSED(K)+FALOUTR(K)/NSTEP
3654       END DO
3655
3656 ! GET PRECIPITATION AND SNOWFALL ACCUMULATION DURING THE TIME STEP
3657 ! FACTOR OF 1000 CONVERTS FROM M TO MM, BUT DIVISION BY DENSITY
3658 ! OF LIQUID WATER CANCELS THIS FACTOR OF 1000
3659
3660         PRECRT = PRECRT+(FALOUTR(KTS)+FALOUTC(KTS)+FALOUTS(KTS)+FALOUTI(KTS)+FALOUTG(KTS))  &
3661                      *DT/NSTEP
3662         SNOWRT = SNOWRT+(FALOUTS(KTS)+FALOUTI(KTS)+FALOUTG(KTS))*DT/NSTEP
3663 ! hm added 7/13/13
3664         SNOWPRT = SNOWPRT+(FALOUTI(KTS)+FALOUTS(KTS))*DT/NSTEP
3665         GRPLPRT = GRPLPRT+(FALOUTG(KTS))*DT/NSTEP
3666
3667       END DO
3668
3669         DO K=KTS,KTE
3670
3671 ! ADD ON SEDIMENTATION TENDENCIES FOR MIXING RATIO TO REST OF TENDENCIES
3672
3673         QR3DTEN(K)=QR3DTEN(K)+QRSTEN(K)
3674         QI3DTEN(K)=QI3DTEN(K)+QISTEN(K)
3675         QC3DTEN(K)=QC3DTEN(K)+QCSTEN(K)
3676         QG3DTEN(K)=QG3DTEN(K)+QGSTEN(K)
3677         QNI3DTEN(K)=QNI3DTEN(K)+QNISTEN(K)
3678
3679 ! PUT ALL CLOUD ICE IN SNOW CATEGORY IF MEAN DIAMETER EXCEEDS 2 * dcs
3680
3681 !hm 4/7/09 bug fix
3682 !        IF (QI3D(K).GE.QSMALL.AND.T3D(K).LT.273.15) THEN
3683         IF (QI3D(K).GE.QSMALL.AND.T3D(K).LT.273.15.AND.LAMI(K).GE.1.E-10) THEN
3684         IF (1./LAMI(K).GE.2.*DCS) THEN
3685            QNI3DTEN(K) = QNI3DTEN(K)+QI3D(K)/DT+ QI3DTEN(K)
3686            NS3DTEN(K) = NS3DTEN(K)+NI3D(K)/DT+   NI3DTEN(K)
3687            QI3DTEN(K) = -QI3D(K)/DT
3688            NI3DTEN(K) = -NI3D(K)/DT
3689         END IF
3690         END IF
3691
3692 ! hm add tendencies here, then call sizeparameter
3693 ! to ensure consisitency between mixing ratio and number concentration
3694
3695           QC3D(k)        = QC3D(k)+QC3DTEN(k)*DT
3696           QI3D(k)        = QI3D(k)+QI3DTEN(k)*DT
3697           QNI3D(k)        = QNI3D(k)+QNI3DTEN(k)*DT
3698           QR3D(k)        = QR3D(k)+QR3DTEN(k)*DT
3699           NC3D(k)        = NC3D(k)+NC3DTEN(k)*DT
3700           NI3D(k)        = NI3D(k)+NI3DTEN(k)*DT
3701           NS3D(k)        = NS3D(k)+NS3DTEN(k)*DT
3702           NR3D(k)        = NR3D(k)+NR3DTEN(k)*DT
3703
3704           IF (IGRAUP.EQ.0) THEN
3705           QG3D(k)        = QG3D(k)+QG3DTEN(k)*DT
3706           NG3D(k)        = NG3D(k)+NG3DTEN(k)*DT
3707           END IF
3708
3709 ! ADD TEMPERATURE AND WATER VAPOR TENDENCIES FROM MICROPHYSICS
3710           T3D(K)         = T3D(K)+T3DTEN(k)*DT
3711           QV3D(K)        = QV3D(K)+QV3DTEN(k)*DT
3712
3713 ! SATURATION VAPOR PRESSURE AND MIXING RATIO
3714
3715 ! hm, add fix for low pressure, 5/12/10
3716             EVS(K) = min(0.99*pres(k),POLYSVP(T3D(K),0))   ! PA
3717             EIS(K) = min(0.99*pres(k),POLYSVP(T3D(K),1))   ! PA
3718
3719 ! MAKE SURE ICE SATURATION DOESN'T EXCEED WATER SAT. NEAR FREEZING
3720
3721             IF (EIS(K).GT.EVS(K)) EIS(K) = EVS(K)
3722
3723             QVS(K) = EP_2*EVS(K)/(PRES(K)-EVS(K))
3724             QVI(K) = EP_2*EIS(K)/(PRES(K)-EIS(K))
3725
3726             QVQVS(K) = QV3D(K)/QVS(K)
3727             QVQVSI(K) = QV3D(K)/QVI(K)
3728
3729 ! AT SUBSATURATION, REMOVE SMALL AMOUNTS OF CLOUD/PRECIP WATER
3730 ! hm 7/9/09 change limit to 1.e-8
3731
3732              IF (QVQVS(K).LT.0.9) THEN
3733                IF (QR3D(K).LT.1.E-8) THEN
3734                   QV3D(K)=QV3D(K)+QR3D(K)
3735                   T3D(K)=T3D(K)-QR3D(K)*XXLV(K)/CPM(K)
3736                   QR3D(K)=0.
3737                END IF
3738                IF (QC3D(K).LT.1.E-8) THEN
3739                   QV3D(K)=QV3D(K)+QC3D(K)
3740                   T3D(K)=T3D(K)-QC3D(K)*XXLV(K)/CPM(K)
3741                   QC3D(K)=0.
3742                END IF
3743              END IF
3744
3745              IF (QVQVSI(K).LT.0.9) THEN
3746                IF (QI3D(K).LT.1.E-8) THEN
3747                   QV3D(K)=QV3D(K)+QI3D(K)
3748                   T3D(K)=T3D(K)-QI3D(K)*XXLS(K)/CPM(K)
3749                   QI3D(K)=0.
3750                END IF
3751                IF (QNI3D(K).LT.1.E-8) THEN
3752                   QV3D(K)=QV3D(K)+QNI3D(K)
3753                   T3D(K)=T3D(K)-QNI3D(K)*XXLS(K)/CPM(K)
3754                   QNI3D(K)=0.
3755                END IF
3756                IF (QG3D(K).LT.1.E-8) THEN
3757                   QV3D(K)=QV3D(K)+QG3D(K)
3758                   T3D(K)=T3D(K)-QG3D(K)*XXLS(K)/CPM(K)
3759                   QG3D(K)=0.
3760                END IF
3761              END IF
3762
3763 !..................................................................
3764 ! IF MIXING RATIO < QSMALL SET MIXING RATIO AND NUMBER CONC TO ZERO
3765
3766        IF (QC3D(K).LT.QSMALL) THEN
3767          QC3D(K) = 0.
3768          NC3D(K) = 0.
3769          EFFC(K) = 0.
3770        END IF
3771        IF (QR3D(K).LT.QSMALL) THEN
3772          QR3D(K) = 0.
3773          NR3D(K) = 0.
3774          EFFR(K) = 0.
3775        END IF
3776        IF (QI3D(K).LT.QSMALL) THEN
3777          QI3D(K) = 0.
3778          NI3D(K) = 0.
3779          EFFI(K) = 0.
3780        END IF
3781        IF (QNI3D(K).LT.QSMALL) THEN
3782          QNI3D(K) = 0.
3783          NS3D(K) = 0.
3784          EFFS(K) = 0.
3785        END IF
3786        IF (QG3D(K).LT.QSMALL) THEN
3787          QG3D(K) = 0.
3788          NG3D(K) = 0.
3789          EFFG(K) = 0.
3790        END IF
3791
3792 !..................................
3793 ! IF THERE IS NO CLOUD/PRECIP WATER, THEN SKIP CALCULATIONS
3794
3795             IF (QC3D(K).LT.QSMALL.AND.QI3D(K).LT.QSMALL.AND.QNI3D(K).LT.QSMALL &
3796                  .AND.QR3D(K).LT.QSMALL.AND.QG3D(K).LT.QSMALL) GOTO 500
3797
3798 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
3799 ! CALCULATE INSTANTANEOUS PROCESSES
3800
3801 ! ADD MELTING OF CLOUD ICE TO FORM RAIN
3802
3803         IF (QI3D(K).GE.QSMALL.AND.T3D(K).GE.273.15) THEN
3804            QR3D(K) = QR3D(K)+QI3D(K)
3805            T3D(K) = T3D(K)-QI3D(K)*XLF(K)/CPM(K)
3806 #if (WRF_CHEM == 1)
3807            tqimelt(K)=QI3D(K)/DT
3808 #endif
3809            QI3D(K) = 0.
3810            NR3D(K) = NR3D(K)+NI3D(K)
3811            NI3D(K) = 0.
3812         END IF
3813
3814 ! ****SENSITIVITY - NO ICE
3815         IF (ILIQ.EQ.1) GOTO 778
3816
3817 ! HOMOGENEOUS FREEZING OF CLOUD WATER
3818
3819         IF (T3D(K).LE.233.15.AND.QC3D(K).GE.QSMALL) THEN
3820            QI3D(K)=QI3D(K)+QC3D(K)
3821            T3D(K)=T3D(K)+QC3D(K)*XLF(K)/CPM(K)
3822            QC3D(K)=0.
3823            NI3D(K)=NI3D(K)+NC3D(K)
3824            NC3D(K)=0.
3825         END IF
3826
3827 ! HOMOGENEOUS FREEZING OF RAIN
3828
3829         IF (IGRAUP.EQ.0) THEN
3830
3831         IF (T3D(K).LE.233.15.AND.QR3D(K).GE.QSMALL) THEN
3832            QG3D(K) = QG3D(K)+QR3D(K)
3833            T3D(K) = T3D(K)+QR3D(K)*XLF(K)/CPM(K)
3834            QR3D(K) = 0.
3835            NG3D(K) = NG3D(K)+ NR3D(K)
3836            NR3D(K) = 0.
3837         END IF
3838
3839         ELSE IF (IGRAUP.EQ.1) THEN
3840
3841         IF (T3D(K).LE.233.15.AND.QR3D(K).GE.QSMALL) THEN
3842            QNI3D(K) = QNI3D(K)+QR3D(K)
3843            T3D(K) = T3D(K)+QR3D(K)*XLF(K)/CPM(K)
3844            QR3D(K) = 0.
3845            NS3D(K) = NS3D(K)+NR3D(K)
3846            NR3D(K) = 0.
3847         END IF
3848
3849         END IF
3850
3851  778    CONTINUE
3852
3853 ! MAKE SURE NUMBER CONCENTRATIONS AREN'T NEGATIVE
3854
3855       NI3D(K) = MAX(0.,NI3D(K))
3856       NS3D(K) = MAX(0.,NS3D(K))
3857       NC3D(K) = MAX(0.,NC3D(K))
3858       NR3D(K) = MAX(0.,NR3D(K))
3859       NG3D(K) = MAX(0.,NG3D(K))
3860
3861 !......................................................................
3862 ! CLOUD ICE
3863
3864       IF (QI3D(K).GE.QSMALL) THEN
3865          LAMI(K) = (CONS12*                 &
3866               NI3D(K)/QI3D(K))**(1./DI)
3867
3868 ! CHECK FOR SLOPE
3869
3870 ! ADJUST VARS
3871
3872       IF (LAMI(K).LT.LAMMINI) THEN
3873
3874       LAMI(K) = LAMMINI
3875
3876       N0I(K) = LAMI(K)**4*QI3D(K)/CONS12
3877
3878       NI3D(K) = N0I(K)/LAMI(K)
3879       ELSE IF (LAMI(K).GT.LAMMAXI) THEN
3880       LAMI(K) = LAMMAXI
3881       N0I(K) = LAMI(K)**4*QI3D(K)/CONS12
3882
3883       NI3D(K) = N0I(K)/LAMI(K)
3884       END IF
3885       END IF
3886
3887 !......................................................................
3888 ! RAIN
3889
3890       IF (QR3D(K).GE.QSMALL) THEN
3891       LAMR(K) = (PI*RHOW*NR3D(K)/QR3D(K))**(1./3.)
3892
3893 ! CHECK FOR SLOPE
3894
3895 ! ADJUST VARS
3896
3897       IF (LAMR(K).LT.LAMMINR) THEN
3898
3899       LAMR(K) = LAMMINR
3900
3901       N0RR(K) = LAMR(K)**4*QR3D(K)/(PI*RHOW)
3902
3903       NR3D(K) = N0RR(K)/LAMR(K)
3904       ELSE IF (LAMR(K).GT.LAMMAXR) THEN
3905       LAMR(K) = LAMMAXR
3906       N0RR(K) = LAMR(K)**4*QR3D(K)/(PI*RHOW)
3907
3908       NR3D(K) = N0RR(K)/LAMR(K)
3909       END IF
3910
3911       END IF
3912
3913 !......................................................................
3914 ! CLOUD DROPLETS
3915
3916 ! MARTIN ET AL. (1994) FORMULA FOR PGAM
3917
3918       IF (QC3D(K).GE.QSMALL) THEN
3919
3920          DUM = PRES(K)/(287.15*T3D(K))
3921          PGAM(K)=0.0005714*(NC3D(K)/1.E6*DUM)+0.2714
3922          PGAM(K)=1./(PGAM(K)**2)-1.
3923          PGAM(K)=MAX(PGAM(K),2.)
3924          PGAM(K)=MIN(PGAM(K),10.)
3925
3926 ! CALCULATE LAMC
3927
3928       LAMC(K) = (CONS26*NC3D(K)*GAMMA(PGAM(K)+4.)/   &
3929                  (QC3D(K)*GAMMA(PGAM(K)+1.)))**(1./3.)
3930
3931 ! LAMMIN, 60 MICRON DIAMETER
3932 ! LAMMAX, 1 MICRON
3933
3934       LAMMIN = (PGAM(K)+1.)/60.E-6
3935       LAMMAX = (PGAM(K)+1.)/1.E-6
3936
3937       IF (LAMC(K).LT.LAMMIN) THEN
3938       LAMC(K) = LAMMIN
3939       NC3D(K) = EXP(3.*LOG(LAMC(K))+LOG(QC3D(K))+              &
3940                 LOG(GAMMA(PGAM(K)+1.))-LOG(GAMMA(PGAM(K)+4.)))/CONS26
3941
3942       ELSE IF (LAMC(K).GT.LAMMAX) THEN
3943       LAMC(K) = LAMMAX
3944       NC3D(K) = EXP(3.*LOG(LAMC(K))+LOG(QC3D(K))+              &
3945                 LOG(GAMMA(PGAM(K)+1.))-LOG(GAMMA(PGAM(K)+4.)))/CONS26
3946
3947       END IF
3948
3949       END IF
3950
3951 !......................................................................
3952 ! SNOW
3953
3954       IF (QNI3D(K).GE.QSMALL) THEN
3955       LAMS(K) = (CONS1*NS3D(K)/QNI3D(K))**(1./DS)
3956
3957 ! CHECK FOR SLOPE
3958
3959 ! ADJUST VARS
3960
3961       IF (LAMS(K).LT.LAMMINS) THEN
3962       LAMS(K) = LAMMINS
3963       N0S(K) = LAMS(K)**4*QNI3D(K)/CONS1
3964
3965       NS3D(K) = N0S(K)/LAMS(K)
3966
3967       ELSE IF (LAMS(K).GT.LAMMAXS) THEN
3968
3969       LAMS(K) = LAMMAXS
3970       N0S(K) = LAMS(K)**4*QNI3D(K)/CONS1
3971       NS3D(K) = N0S(K)/LAMS(K)
3972       END IF
3973
3974       END IF
3975
3976 !......................................................................
3977 ! GRAUPEL
3978
3979       IF (QG3D(K).GE.QSMALL) THEN
3980       LAMG(K) = (CONS2*NG3D(K)/QG3D(K))**(1./DG)
3981
3982 ! CHECK FOR SLOPE
3983
3984 ! ADJUST VARS
3985
3986       IF (LAMG(K).LT.LAMMING) THEN
3987       LAMG(K) = LAMMING
3988       N0G(K) = LAMG(K)**4*QG3D(K)/CONS2
3989
3990       NG3D(K) = N0G(K)/LAMG(K)
3991
3992       ELSE IF (LAMG(K).GT.LAMMAXG) THEN
3993
3994       LAMG(K) = LAMMAXG
3995       N0G(K) = LAMG(K)**4*QG3D(K)/CONS2
3996
3997       NG3D(K) = N0G(K)/LAMG(K)
3998       END IF
3999
4000       END IF
4001
4002  500  CONTINUE
4003
4004 ! CALCULATE EFFECTIVE RADIUS
4005
4006       IF (QI3D(K).GE.QSMALL) THEN
4007          EFFI(K) = 3./LAMI(K)/2.*1.E6
4008       ELSE
4009          EFFI(K) = 25.
4010       END IF
4011
4012       IF (QNI3D(K).GE.QSMALL) THEN
4013          EFFS(K) = 3./LAMS(K)/2.*1.E6
4014       ELSE
4015          EFFS(K) = 25.
4016       END IF
4017
4018       IF (QR3D(K).GE.QSMALL) THEN
4019          EFFR(K) = 3./LAMR(K)/2.*1.E6
4020       ELSE
4021          EFFR(K) = 25.
4022       END IF
4023
4024       IF (QC3D(K).GE.QSMALL) THEN
4025       EFFC(K) = GAMMA(PGAM(K)+4.)/                        &
4026              GAMMA(PGAM(K)+3.)/LAMC(K)/2.*1.E6
4027       ELSE
4028       EFFC(K) = 25.
4029       END IF
4030
4031       IF (QG3D(K).GE.QSMALL) THEN
4032          EFFG(K) = 3./LAMG(K)/2.*1.E6
4033       ELSE
4034          EFFG(K) = 25.
4035       END IF
4036
4037 ! HM ADD 1/10/06, ADD UPPER BOUND ON ICE NUMBER, THIS IS NEEDED
4038 ! TO PREVENT VERY LARGE ICE NUMBER DUE TO HOMOGENEOUS FREEZING
4039 ! OF DROPLETS, ESPECIALLY WHEN INUM = 1, SET MAX AT 10 CM-3
4040 !          NI3D(K) = MIN(NI3D(K),10.E6/RHO(K))
4041 ! HM, 12/28/12, LOWER MAXIMUM ICE CONCENTRATION TO ADDRESS PROBLEM
4042 ! OF EXCESSIVE AND PERSISTENT ANVIL
4043 ! NOTE: THIS MAY CHANGE/REDUCE SENSITIVITY TO AEROSOL/CCN CONCENTRATION
4044           NI3D(K) = MIN(NI3D(K),0.3E6/RHO(K))
4045
4046 ! ADD BOUND ON DROPLET NUMBER - CANNOT EXCEED AEROSOL CONCENTRATION
4047           IF (iinum.EQ.0.AND.IACT.EQ.2) THEN
4048           NC3D(K) = MIN(NC3D(K),(NANEW1+NANEW2)/RHO(K))
4049           END IF
4050 ! SWITCH FOR CONSTANT DROPLET NUMBER
4051           IF (iinum.EQ.1) THEN
4052 ! CHANGE NDCNST FROM CM-3 TO KG-1
4053              NC3D(K) = NDCNST*1.E6/RHO(K)
4054           END IF
4055
4056       END DO !!! K LOOP
4057
4058  400         CONTINUE
4059
4060 ! ALL DONE !!!!!!!!!!!
4061       RETURN
4062       END SUBROUTINE MORR_TWO_MOMENT_MICRO
4063
4064 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
4065
4066       REAL FUNCTION POLYSVP (T,TYPE)
4067
4068 !-------------------------------------------
4069
4070 !  COMPUTE SATURATION VAPOR PRESSURE
4071
4072 !  POLYSVP RETURNED IN UNITS OF PA.
4073 !  T IS INPUT IN UNITS OF K.
4074 !  TYPE REFERS TO SATURATION WITH RESPECT TO LIQUID (0) OR ICE (1)
4075
4076 ! REPLACE GOFF-GRATCH WITH FASTER FORMULATION FROM FLATAU ET AL. 1992, TABLE 4 (RIGHT-HAND COLUMN)
4077
4078       IMPLICIT NONE
4079
4080       REAL DUM
4081       REAL T
4082       INTEGER TYPE
4083 ! ice
4084       real a0i,a1i,a2i,a3i,a4i,a5i,a6i,a7i,a8i
4085       data a0i,a1i,a2i,a3i,a4i,a5i,a6i,a7i,a8i /&
4086         6.11147274, 0.503160820, 0.188439774e-1, &
4087         0.420895665e-3, 0.615021634e-5,0.602588177e-7, &
4088         0.385852041e-9, 0.146898966e-11, 0.252751365e-14/
4089
4090 ! liquid
4091       real a0,a1,a2,a3,a4,a5,a6,a7,a8
4092
4093 ! V1.7
4094       data a0,a1,a2,a3,a4,a5,a6,a7,a8 /&
4095         6.11239921, 0.443987641, 0.142986287e-1, &
4096         0.264847430e-3, 0.302950461e-5, 0.206739458e-7, &
4097         0.640689451e-10,-0.952447341e-13,-0.976195544e-15/
4098       real dt
4099
4100 ! ICE
4101
4102       IF (TYPE.EQ.1) THEN
4103
4104 !         POLYSVP = 10.**(-9.09718*(273.16/T-1.)-3.56654*                &
4105 !          LOG10(273.16/T)+0.876793*(1.-T/273.16)+                                              &
4106 !          LOG10(6.1071))*100.
4107 ! hm 11/16/20, use Goff-Gratch for T < 195.8 K and Flatau et al. equal or above 195.8 K
4108       if (t.ge.195.8) then
4109          dt=t-273.15
4110          polysvp = a0i + dt*(a1i+dt*(a2i+dt*(a3i+dt*(a4i+dt*(a5i+dt*(a6i+dt*(a7i+a8i*dt)))))))
4111          polysvp = polysvp*100.
4112       else
4113
4114          polysvp = 10.**(-9.09718*(273.16/t-1.)-3.56654* &
4115           alog10(273.16/t)+0.876793*(1.-t/273.16)+ &
4116           alog10(6.1071))*100.
4117
4118       end if
4119
4120       END IF
4121
4122 ! LIQUID
4123
4124       IF (TYPE.EQ.0) THEN
4125
4126 !         POLYSVP = 10.**(-7.90298*(373.16/T-1.)+                        &
4127 !             5.02808*LOG10(373.16/T)-                                                                  &
4128 !             1.3816E-7*(10**(11.344*(1.-T/373.16))-1.)+                                &
4129 !             8.1328E-3*(10**(-3.49149*(373.16/T-1.))-1.)+                              &
4130 !             LOG10(1013.246))*100.
4131 ! hm 11/16/20, use Goff-Gratch for T < 202.0 K and Flatau et al. equal or above 202.0 K
4132       if (t.ge.202.0) then
4133          dt = t-273.15
4134          polysvp = a0 + dt*(a1+dt*(a2+dt*(a3+dt*(a4+dt*(a5+dt*(a6+dt*(a7+a8*dt)))))))
4135          polysvp = polysvp*100.
4136       else
4137
4138 ! note: uncertain below -70 C, but produces physical values (non-negative) unlike flatau
4139          polysvp = 10.**(-7.90298*(373.16/t-1.)+ &
4140              5.02808*alog10(373.16/t)- &
4141              1.3816e-7*(10**(11.344*(1.-t/373.16))-1.)+ &
4142              8.1328e-3*(10**(-3.49149*(373.16/t-1.))-1.)+ &
4143              alog10(1013.246))*100.
4144       end if
4145
4146       END IF
4147
4148
4149       END FUNCTION POLYSVP
4150
4151 !------------------------------------------------------------------------------
4152
4153       REAL FUNCTION GAMMA(X)
4154 !----------------------------------------------------------------------
4155 !
4156 ! THIS ROUTINE CALCULATES THE GAMMA FUNCTION FOR A REAL ARGUMENT X.
4157 !   COMPUTATION IS BASED ON AN ALGORITHM OUTLINED IN REFERENCE 1.
4158 !   THE PROGRAM USES RATIONAL FUNCTIONS THAT APPROXIMATE THE GAMMA
4159 !   FUNCTION TO AT LEAST 20 SIGNIFICANT DECIMAL DIGITS.  COEFFICIENTS
4160 !   FOR THE APPROXIMATION OVER THE INTERVAL (1,2) ARE UNPUBLISHED.
4161 !   THOSE FOR THE APPROXIMATION FOR X .GE. 12 ARE FROM REFERENCE 2.
4162 !   THE ACCURACY ACHIEVED DEPENDS ON THE ARITHMETIC SYSTEM, THE
4163 !   COMPILER, THE INTRINSIC FUNCTIONS, AND PROPER SELECTION OF THE
4164 !   MACHINE-DEPENDENT CONSTANTS.
4165 !
4166 !
4167 !*******************************************************************
4168 !*******************************************************************
4169 !
4170 ! EXPLANATION OF MACHINE-DEPENDENT CONSTANTS
4171 !
4172 ! BETA   - RADIX FOR THE FLOATING-POINT REPRESENTATION
4173 ! MAXEXP - THE SMALLEST POSITIVE POWER OF BETA THAT OVERFLOWS
4174 ! XBIG   - THE LARGEST ARGUMENT FOR WHICH GAMMA(X) IS REPRESENTABLE
4175 !          IN THE MACHINE, I.E., THE SOLUTION TO THE EQUATION
4176 !                  GAMMA(XBIG) = BETA**MAXEXP
4177 ! XINF   - THE LARGEST MACHINE REPRESENTABLE FLOATING-POINT NUMBER;
4178 !          APPROXIMATELY BETA**MAXEXP
4179 ! EPS    - THE SMALLEST POSITIVE FLOATING-POINT NUMBER SUCH THAT
4180 !          1.0+EPS .GT. 1.0
4181 ! XMININ - THE SMALLEST POSITIVE FLOATING-POINT NUMBER SUCH THAT
4182 !          1/XMININ IS MACHINE REPRESENTABLE
4183 !
4184 !     APPROXIMATE VALUES FOR SOME IMPORTANT MACHINES ARE:
4185 !
4186 !                            BETA       MAXEXP        XBIG
4187 !
4188 ! CRAY-1         (S.P.)        2         8191        966.961
4189 ! CYBER 180/855
4190 !   UNDER NOS    (S.P.)        2         1070        177.803
4191 ! IEEE (IBM/XT,
4192 !   SUN, ETC.)   (S.P.)        2          128        35.040
4193 ! IEEE (IBM/XT,
4194 !   SUN, ETC.)   (D.P.)        2         1024        171.624
4195 ! IBM 3033       (D.P.)       16           63        57.574
4196 ! VAX D-FORMAT   (D.P.)        2          127        34.844
4197 ! VAX G-FORMAT   (D.P.)        2         1023        171.489
4198 !
4199 !                            XINF         EPS        XMININ
4200 !
4201 ! CRAY-1         (S.P.)   5.45E+2465   7.11E-15    1.84E-2466
4202 ! CYBER 180/855
4203 !   UNDER NOS    (S.P.)   1.26E+322    3.55E-15    3.14E-294
4204 ! IEEE (IBM/XT,
4205 !   SUN, ETC.)   (S.P.)   3.40E+38     1.19E-7     1.18E-38
4206 ! IEEE (IBM/XT,
4207 !   SUN, ETC.)   (D.P.)   1.79D+308    2.22D-16    2.23D-308
4208 ! IBM 3033       (D.P.)   7.23D+75     2.22D-16    1.39D-76
4209 ! VAX D-FORMAT   (D.P.)   1.70D+38     1.39D-17    5.88D-39
4210 ! VAX G-FORMAT   (D.P.)   8.98D+307    1.11D-16    1.12D-308
4211 !
4212 !*******************************************************************
4213 !*******************************************************************
4214 !
4215 ! ERROR RETURNS
4216 !
4217 !  THE PROGRAM RETURNS THE VALUE XINF FOR SINGULARITIES OR
4218 !     WHEN OVERFLOW WOULD OCCUR.  THE COMPUTATION IS BELIEVED
4219 !     TO BE FREE OF UNDERFLOW AND OVERFLOW.
4220 !
4221 !
4222 !  INTRINSIC FUNCTIONS REQUIRED ARE:
4223 !
4224 !     INT, DBLE, EXP, LOG, REAL, SIN
4225 !
4226 !
4227 ! REFERENCES:  AN OVERVIEW OF SOFTWARE DEVELOPMENT FOR SPECIAL
4228 !              FUNCTIONS   W. J. CODY, LECTURE NOTES IN MATHEMATICS,
4229 !              506, NUMERICAL ANALYSIS DUNDEE, 1975, G. A. WATSON
4230 !              (ED.), SPRINGER VERLAG, BERLIN, 1976.
4231 !
4232 !              COMPUTER APPROXIMATIONS, HART, ET. AL., WILEY AND
4233 !              SONS, NEW YORK, 1968.
4234 !
4235 !  LATEST MODIFICATION: OCTOBER 12, 1989
4236 !
4237 !  AUTHORS: W. J. CODY AND L. STOLTZ
4238 !           APPLIED MATHEMATICS DIVISION
4239 !           ARGONNE NATIONAL LABORATORY
4240 !           ARGONNE, IL 60439
4241 !
4242 !----------------------------------------------------------------------
4243       implicit none
4244       INTEGER I,N
4245       LOGICAL PARITY
4246       REAL                                                          &
4247           CONV,EPS,FACT,HALF,ONE,RES,SUM,TWELVE,                    &
4248           TWO,X,XBIG,XDEN,XINF,XMININ,XNUM,Y,Y1,YSQ,Z,ZERO
4249       REAL, DIMENSION(7) :: C
4250       REAL, DIMENSION(8) :: P
4251       REAL, DIMENSION(8) :: Q
4252 !----------------------------------------------------------------------
4253 !  MATHEMATICAL CONSTANTS
4254 !----------------------------------------------------------------------
4255       DATA ONE,HALF,TWELVE,TWO,ZERO/1.0E0,0.5E0,12.0E0,2.0E0,0.0E0/
4256
4257
4258 !----------------------------------------------------------------------
4259 !  MACHINE DEPENDENT PARAMETERS
4260 !----------------------------------------------------------------------
4261       DATA XBIG,XMININ,EPS/35.040E0,1.18E-38,1.19E-7/,XINF/3.4E38/
4262 !----------------------------------------------------------------------
4263 !  NUMERATOR AND DENOMINATOR COEFFICIENTS FOR RATIONAL MINIMAX
4264 !     APPROXIMATION OVER (1,2).
4265 !----------------------------------------------------------------------
4266       DATA P/-1.71618513886549492533811E+0,2.47656508055759199108314E+1,  &
4267              -3.79804256470945635097577E+2,6.29331155312818442661052E+2,  &
4268              8.66966202790413211295064E+2,-3.14512729688483675254357E+4,  &
4269              -3.61444134186911729807069E+4,6.64561438202405440627855E+4/
4270       DATA Q/-3.08402300119738975254353E+1,3.15350626979604161529144E+2,  &
4271              -1.01515636749021914166146E+3,-3.10777167157231109440444E+3, &
4272               2.25381184209801510330112E+4,4.75584627752788110767815E+3,  &
4273             -1.34659959864969306392456E+5,-1.15132259675553483497211E+5/
4274 !----------------------------------------------------------------------
4275 !  COEFFICIENTS FOR MINIMAX APPROXIMATION OVER (12, INF).
4276 !----------------------------------------------------------------------
4277       DATA C/-1.910444077728E-03,8.4171387781295E-04,                      &
4278            -5.952379913043012E-04,7.93650793500350248E-04,                                 &
4279            -2.777777777777681622553E-03,8.333333333333333331554247E-02,    &
4280             5.7083835261E-03/
4281 !----------------------------------------------------------------------
4282 !  STATEMENT FUNCTIONS FOR CONVERSION BETWEEN INTEGER AND FLOAT
4283 !----------------------------------------------------------------------
4284       CONV(I) = REAL(I)
4285       PARITY=.FALSE.
4286       FACT=ONE
4287       N=0
4288       Y=X
4289       IF(Y.LE.ZERO)THEN
4290 !----------------------------------------------------------------------
4291 !  ARGUMENT IS NEGATIVE
4292 !----------------------------------------------------------------------
4293         Y=-X
4294         Y1=AINT(Y)
4295         RES=Y-Y1
4296         IF(RES.NE.ZERO)THEN
4297           IF(Y1.NE.AINT(Y1*HALF)*TWO)PARITY=.TRUE.
4298           FACT=-PI/SIN(PI*RES)
4299           Y=Y+ONE
4300         ELSE
4301           RES=XINF
4302           GOTO 900
4303         ENDIF
4304       ENDIF
4305 !----------------------------------------------------------------------
4306 !  ARGUMENT IS POSITIVE
4307 !----------------------------------------------------------------------
4308       IF(Y.LT.EPS)THEN
4309 !----------------------------------------------------------------------
4310 !  ARGUMENT .LT. EPS
4311 !----------------------------------------------------------------------
4312         IF(Y.GE.XMININ)THEN
4313           RES=ONE/Y
4314         ELSE
4315           RES=XINF
4316           GOTO 900
4317         ENDIF
4318       ELSEIF(Y.LT.TWELVE)THEN
4319         Y1=Y
4320         IF(Y.LT.ONE)THEN
4321 !----------------------------------------------------------------------
4322 !  0.0 .LT. ARGUMENT .LT. 1.0
4323 !----------------------------------------------------------------------
4324           Z=Y
4325           Y=Y+ONE
4326         ELSE
4327 !----------------------------------------------------------------------
4328 !  1.0 .LT. ARGUMENT .LT. 12.0, REDUCE ARGUMENT IF NECESSARY
4329 !----------------------------------------------------------------------
4330           N=INT(Y)-1
4331           Y=Y-CONV(N)
4332           Z=Y-ONE
4333         ENDIF
4334 !----------------------------------------------------------------------
4335 !  EVALUATE APPROXIMATION FOR 1.0 .LT. ARGUMENT .LT. 2.0
4336 !----------------------------------------------------------------------
4337         XNUM=ZERO
4338         XDEN=ONE
4339         DO I=1,8
4340           XNUM=(XNUM+P(I))*Z
4341           XDEN=XDEN*Z+Q(I)
4342         END DO
4343         RES=XNUM/XDEN+ONE
4344         IF(Y1.LT.Y)THEN
4345 !----------------------------------------------------------------------
4346 !  ADJUST RESULT FOR CASE  0.0 .LT. ARGUMENT .LT. 1.0
4347 !----------------------------------------------------------------------
4348           RES=RES/Y1
4349         ELSEIF(Y1.GT.Y)THEN
4350 !----------------------------------------------------------------------
4351 !  ADJUST RESULT FOR CASE  2.0 .LT. ARGUMENT .LT. 12.0
4352 !----------------------------------------------------------------------
4353           DO I=1,N
4354             RES=RES*Y
4355             Y=Y+ONE
4356           END DO
4357         ENDIF
4358       ELSE
4359 !----------------------------------------------------------------------
4360 !  EVALUATE FOR ARGUMENT .GE. 12.0,
4361 !----------------------------------------------------------------------
4362         IF(Y.LE.XBIG)THEN
4363           YSQ=Y*Y
4364           SUM=C(7)
4365           DO I=1,6
4366             SUM=SUM/YSQ+C(I)
4367           END DO
4368           SUM=SUM/Y-Y+xxx
4369           SUM=SUM+(Y-HALF)*LOG(Y)
4370           RES=EXP(SUM)
4371         ELSE
4372           RES=XINF
4373           GOTO 900
4374         ENDIF
4375       ENDIF
4376 !----------------------------------------------------------------------
4377 !  FINAL ADJUSTMENTS AND RETURN
4378 !----------------------------------------------------------------------
4379       IF(PARITY)RES=-RES
4380       IF(FACT.NE.ONE)RES=FACT/RES
4381   900 GAMMA=RES
4382       RETURN
4383 ! ---------- LAST LINE OF GAMMA ----------
4384       END FUNCTION GAMMA
4385
4386
4387       REAL FUNCTION DERF1(X)
4388       IMPLICIT NONE
4389       REAL X
4390       REAL, DIMENSION(0 : 64) :: A, B
4391       REAL W,T,Y
4392       INTEGER K,I
4393       DATA A/                                                 &
4394          0.00000000005958930743E0, -0.00000000113739022964E0, &
4395          0.00000001466005199839E0, -0.00000016350354461960E0, &
4396          0.00000164610044809620E0, -0.00001492559551950604E0, &
4397          0.00012055331122299265E0, -0.00085483269811296660E0, &
4398          0.00522397762482322257E0, -0.02686617064507733420E0, &
4399          0.11283791670954881569E0, -0.37612638903183748117E0, &
4400          1.12837916709551257377E0,                                &
4401          0.00000000002372510631E0, -0.00000000045493253732E0, &
4402          0.00000000590362766598E0, -0.00000006642090827576E0, &
4403          0.00000067595634268133E0, -0.00000621188515924000E0, &
4404          0.00005103883009709690E0, -0.00037015410692956173E0, &
4405          0.00233307631218880978E0, -0.01254988477182192210E0, &
4406          0.05657061146827041994E0, -0.21379664776456006580E0, &
4407          0.84270079294971486929E0,                                                        &
4408          0.00000000000949905026E0, -0.00000000018310229805E0, &
4409          0.00000000239463074000E0, -0.00000002721444369609E0, &
4410          0.00000028045522331686E0, -0.00000261830022482897E0, &
4411          0.00002195455056768781E0, -0.00016358986921372656E0, &
4412          0.00107052153564110318E0, -0.00608284718113590151E0, &
4413          0.02986978465246258244E0, -0.13055593046562267625E0, &
4414          0.67493323603965504676E0,                                                        &
4415          0.00000000000382722073E0, -0.00000000007421598602E0, &
4416          0.00000000097930574080E0, -0.00000001126008898854E0, &
4417          0.00000011775134830784E0, -0.00000111992758382650E0, &
4418          0.00000962023443095201E0, -0.00007404402135070773E0, &
4419          0.00050689993654144881E0, -0.00307553051439272889E0, &
4420          0.01668977892553165586E0, -0.08548534594781312114E0, &
4421          0.56909076642393639985E0,                                                        &
4422          0.00000000000155296588E0, -0.00000000003032205868E0, &
4423          0.00000000040424830707E0, -0.00000000471135111493E0, &
4424          0.00000005011915876293E0, -0.00000048722516178974E0, &
4425          0.00000430683284629395E0, -0.00003445026145385764E0, &
4426          0.00024879276133931664E0, -0.00162940941748079288E0, &
4427          0.00988786373932350462E0, -0.05962426839442303805E0, &
4428          0.49766113250947636708E0 /
4429       DATA (B(I), I = 0, 12) /                                  &
4430          -0.00000000029734388465E0,  0.00000000269776334046E0,  &
4431          -0.00000000640788827665E0, -0.00000001667820132100E0,  &
4432          -0.00000021854388148686E0,  0.00000266246030457984E0,  &
4433           0.00001612722157047886E0, -0.00025616361025506629E0,  &
4434           0.00015380842432375365E0,  0.00815533022524927908E0,  &
4435          -0.01402283663896319337E0, -0.19746892495383021487E0,  &
4436           0.71511720328842845913E0 /
4437       DATA (B(I), I = 13, 25) /                                 &
4438          -0.00000000001951073787E0, -0.00000000032302692214E0,  &
4439           0.00000000522461866919E0,  0.00000000342940918551E0,  &
4440          -0.00000035772874310272E0,  0.00000019999935792654E0,  &
4441           0.00002687044575042908E0, -0.00011843240273775776E0,  &
4442          -0.00080991728956032271E0,  0.00661062970502241174E0,  &
4443           0.00909530922354827295E0, -0.20160072778491013140E0,  &
4444           0.51169696718727644908E0 /
4445       DATA (B(I), I = 26, 38) /                                 &
4446          0.00000000003147682272E0, -0.00000000048465972408E0,   &
4447          0.00000000063675740242E0,  0.00000003377623323271E0,   &
4448         -0.00000015451139637086E0, -0.00000203340624738438E0,   &
4449          0.00001947204525295057E0,  0.00002854147231653228E0,   &
4450         -0.00101565063152200272E0,  0.00271187003520095655E0,   &
4451          0.02328095035422810727E0, -0.16725021123116877197E0,   &
4452          0.32490054966649436974E0 /
4453       DATA (B(I), I = 39, 51) /                                 &
4454          0.00000000002319363370E0, -0.00000000006303206648E0,   &
4455         -0.00000000264888267434E0,  0.00000002050708040581E0,   &
4456          0.00000011371857327578E0, -0.00000211211337219663E0,   &
4457          0.00000368797328322935E0,  0.00009823686253424796E0,   &
4458         -0.00065860243990455368E0, -0.00075285814895230877E0,   &
4459          0.02585434424202960464E0, -0.11637092784486193258E0,   &
4460          0.18267336775296612024E0 /
4461       DATA (B(I), I = 52, 64) /                                 &
4462         -0.00000000000367789363E0,  0.00000000020876046746E0,   &
4463         -0.00000000193319027226E0, -0.00000000435953392472E0,   &
4464          0.00000018006992266137E0, -0.00000078441223763969E0,   &
4465         -0.00000675407647949153E0,  0.00008428418334440096E0,   &
4466         -0.00017604388937031815E0, -0.00239729611435071610E0,   &
4467          0.02064129023876022970E0, -0.06905562880005864105E0,   &
4468          0.09084526782065478489E0 /
4469       W = ABS(X)
4470       IF (W .LT. 2.2D0) THEN
4471           T = W * W
4472           K = INT(T)
4473           T = T - K
4474           K = K * 13
4475           Y = ((((((((((((A(K) * T + A(K + 1)) * T +              &
4476               A(K + 2)) * T + A(K + 3)) * T + A(K + 4)) * T +     &
4477               A(K + 5)) * T + A(K + 6)) * T + A(K + 7)) * T +     &
4478               A(K + 8)) * T + A(K + 9)) * T + A(K + 10)) * T +    &
4479               A(K + 11)) * T + A(K + 12)) * W
4480       ELSE IF (W .LT. 6.9D0) THEN
4481           K = INT(W)
4482           T = W - K
4483           K = 13 * (K - 2)
4484           Y = (((((((((((B(K) * T + B(K + 1)) * T +               &
4485               B(K + 2)) * T + B(K + 3)) * T + B(K + 4)) * T +     &
4486               B(K + 5)) * T + B(K + 6)) * T + B(K + 7)) * T +     &
4487               B(K + 8)) * T + B(K + 9)) * T + B(K + 10)) * T +    &
4488               B(K + 11)) * T + B(K + 12)
4489           Y = Y * Y
4490           Y = Y * Y
4491           Y = Y * Y
4492           Y = 1 - Y * Y
4493       ELSE
4494           Y = 1
4495       END IF
4496       IF (X .LT. 0) Y = -Y
4497       DERF1 = Y
4498       END FUNCTION DERF1
4499
4500 !+---+-----------------------------------------------------------------+
4501
4502       subroutine refl10cm_hm (qv1d, qr1d, nr1d, qs1d, ns1d, qg1d, ng1d, &
4503                       t1d, p1d, dBZ, kts, kte, ii, jj)
4504
4505       IMPLICIT NONE
4506
4507 !..Sub arguments
4508       INTEGER, INTENT(IN):: kts, kte, ii, jj
4509       REAL, DIMENSION(kts:kte), INTENT(IN)::                            &
4510                       qv1d, qr1d, nr1d, qs1d, ns1d, qg1d, ng1d, t1d, p1d
4511       REAL, DIMENSION(kts:kte), INTENT(INOUT):: dBZ
4512
4513 !..Local variables
4514       REAL, DIMENSION(kts:kte):: temp, pres, qv, rho
4515       REAL, DIMENSION(kts:kte):: rr, nr, rs, ns, rg, ng
4516
4517       DOUBLE PRECISION, DIMENSION(kts:kte):: ilamr, ilamg, ilams
4518       DOUBLE PRECISION, DIMENSION(kts:kte):: N0_r, N0_g, N0_s
4519       DOUBLE PRECISION:: lamr, lamg, lams
4520       LOGICAL, DIMENSION(kts:kte):: L_qr, L_qs, L_qg
4521
4522       REAL, DIMENSION(kts:kte):: ze_rain, ze_snow, ze_graupel
4523       DOUBLE PRECISION:: fmelt_s, fmelt_g
4524       DOUBLE PRECISION:: cback, x, eta, f_d
4525
4526       INTEGER:: i, k, k_0, kbot, n
4527       LOGICAL:: melti
4528
4529 !+---+
4530
4531       do k = kts, kte
4532          dBZ(k) = -35.0
4533       enddo
4534
4535 !+---+-----------------------------------------------------------------+
4536 !..Put column of data into local arrays.
4537 !+---+-----------------------------------------------------------------+
4538       do k = kts, kte
4539          temp(k) = t1d(k)
4540          qv(k) = MAX(1.E-10, qv1d(k))
4541          pres(k) = p1d(k)
4542          rho(k) = 0.622*pres(k)/(R*temp(k)*(qv(k)+0.622))
4543
4544          if (qr1d(k) .gt. 1.E-9) then
4545             rr(k) = qr1d(k)*rho(k)
4546             nr(k) = nr1d(k)*rho(k)
4547             lamr = (xam_r*xcrg(3)*xorg2*nr(k)/rr(k))**xobmr
4548             ilamr(k) = 1./lamr
4549             N0_r(k) = nr(k)*xorg2*lamr**xcre(2)
4550             L_qr(k) = .true.
4551          else
4552             rr(k) = 1.E-12
4553             nr(k) = 1.E-12
4554             L_qr(k) = .false.
4555          endif
4556
4557          if (qs1d(k) .gt. 1.E-9) then
4558             rs(k) = qs1d(k)*rho(k)
4559             ns(k) = ns1d(k)*rho(k)
4560             lams = (xam_s*xcsg(3)*xosg2*ns(k)/rs(k))**xobms
4561             ilams(k) = 1./lams
4562             N0_s(k) = ns(k)*xosg2*lams**xcse(2)
4563             L_qs(k) = .true.
4564          else
4565             rs(k) = 1.E-12
4566             ns(k) = 1.E-12
4567             L_qs(k) = .false.
4568          endif
4569
4570          if (qg1d(k) .gt. 1.E-9) then
4571             rg(k) = qg1d(k)*rho(k)
4572             ng(k) = ng1d(k)*rho(k)
4573             lamg = (xam_g*xcgg(3)*xogg2*ng(k)/rg(k))**xobmg
4574             ilamg(k) = 1./lamg
4575             N0_g(k) = ng(k)*xogg2*lamg**xcge(2)
4576             L_qg(k) = .true.
4577          else
4578             rg(k) = 1.E-12
4579             ng(k) = 1.E-12
4580             L_qg(k) = .false.
4581          endif
4582       enddo
4583
4584 !+---+-----------------------------------------------------------------+
4585 !..Locate K-level of start of melting (k_0 is level above).
4586 !+---+-----------------------------------------------------------------+
4587       melti = .false.
4588       k_0 = kts
4589       do k = kte-1, kts, -1
4590          if ( (temp(k).gt.273.15) .and. L_qr(k)                         &
4591                                   .and. (L_qs(k+1).or.L_qg(k+1)) ) then
4592             k_0 = MAX(k+1, k_0)
4593             melti=.true.
4594             goto 195
4595          endif
4596       enddo
4597  195  continue
4598
4599 !+---+-----------------------------------------------------------------+
4600 !..Assume Rayleigh approximation at 10 cm wavelength. Rain (all temps)
4601 !.. and non-water-coated snow and graupel when below freezing are
4602 !.. simple. Integrations of m(D)*m(D)*N(D)*dD.
4603 !+---+-----------------------------------------------------------------+
4604
4605       do k = kts, kte
4606          ze_rain(k) = 1.e-22
4607          ze_snow(k) = 1.e-22
4608          ze_graupel(k) = 1.e-22
4609          if (L_qr(k)) ze_rain(k) = N0_r(k)*xcrg(4)*ilamr(k)**xcre(4)
4610          if (L_qs(k)) ze_snow(k) = (0.176/0.93) * (6.0/PI)*(6.0/PI)     &
4611                                  * (xam_s/900.0)*(xam_s/900.0)          &
4612                                  * N0_s(k)*xcsg(4)*ilams(k)**xcse(4)
4613          if (L_qg(k)) ze_graupel(k) = (0.176/0.93) * (6.0/PI)*(6.0/PI)  &
4614                                     * (xam_g/900.0)*(xam_g/900.0)       &
4615                                     * N0_g(k)*xcgg(4)*ilamg(k)**xcge(4)
4616       enddo
4617
4618 !+---+-----------------------------------------------------------------+
4619 !..Special case of melting ice (snow/graupel) particles.  Assume the
4620 !.. ice is surrounded by the liquid water.  Fraction of meltwater is
4621 !.. extremely simple based on amount found above the melting level.
4622 !.. Uses code from Uli Blahak (rayleigh_soak_wetgraupel and supporting
4623 !.. routines).
4624 !+---+-----------------------------------------------------------------+
4625
4626       if (melti .and. k_0.ge.kts+1) then
4627        do k = k_0-1, kts, -1
4628
4629 !..Reflectivity contributed by melting snow
4630           if (L_qs(k) .and. L_qs(k_0) ) then
4631            fmelt_s = MAX(0.005d0, MIN(1.0d0-rs(k)/rs(k_0), 0.99d0))
4632            eta = 0.d0
4633            lams = 1./ilams(k)
4634            do n = 1, nrbins
4635               x = xam_s * xxDs(n)**xbm_s
4636               call rayleigh_soak_wetgraupel (x,DBLE(xocms),DBLE(xobms), &
4637                     fmelt_s, melt_outside_s, m_w_0, m_i_0, lamda_radar, &
4638                     CBACK, mixingrulestring_s, matrixstring_s,          &
4639                     inclusionstring_s, hoststring_s,                    &
4640                     hostmatrixstring_s, hostinclusionstring_s)
4641               f_d = N0_s(k)*xxDs(n)**xmu_s * DEXP(-lams*xxDs(n))
4642               eta = eta + f_d * CBACK * simpson(n) * xdts(n)
4643            enddo
4644            ze_snow(k) = SNGL(lamda4 / (pi5 * K_w) * eta)
4645           endif
4646
4647
4648 !..Reflectivity contributed by melting graupel
4649
4650           if (L_qg(k) .and. L_qg(k_0) ) then
4651            fmelt_g = MAX(0.005d0, MIN(1.0d0-rg(k)/rg(k_0), 0.99d0))
4652            eta = 0.d0
4653            lamg = 1./ilamg(k)
4654            do n = 1, nrbins
4655               x = xam_g * xxDg(n)**xbm_g
4656               call rayleigh_soak_wetgraupel (x,DBLE(xocmg),DBLE(xobmg), &
4657                     fmelt_g, melt_outside_g, m_w_0, m_i_0, lamda_radar, &
4658                     CBACK, mixingrulestring_g, matrixstring_g,          &
4659                     inclusionstring_g, hoststring_g,                    &
4660                     hostmatrixstring_g, hostinclusionstring_g)
4661               f_d = N0_g(k)*xxDg(n)**xmu_g * DEXP(-lamg*xxDg(n))
4662               eta = eta + f_d * CBACK * simpson(n) * xdtg(n)
4663            enddo
4664            ze_graupel(k) = SNGL(lamda4 / (pi5 * K_w) * eta)
4665           endif
4666
4667        enddo
4668       endif
4669
4670       do k = kte, kts, -1
4671          dBZ(k) = 10.*log10((ze_rain(k)+ze_snow(k)+ze_graupel(k))*1.d18)
4672       enddo
4673
4674
4675       end subroutine refl10cm_hm
4676
4677 !+---+-----------------------------------------------------------------+
4678
4679 END MODULE module_mp_morr_two_moment
4680 !+---+-----------------------------------------------------------------+
4681 !+---+-----------------------------------------------------------------+