1 !WRF:MODEL_LAYER:PHYSICS
4 ! THIS MODULE CONTAINS THE TWO-MOMENT MICROPHYSICS CODE DESCRIBED BY
5 ! MORRISON ET AL. (2009, MWR)
7 ! CHANGES FOR V3.2, RELATIVE TO MOST RECENT (BUG-FIX) CODE FOR V3.1
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
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
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
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
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)
46 ! 1) very minor change to limits on autoconversion source of rain number when cloud water is depleted
49 ! hm/A. Ackerman bug fix 11/08/12
51 ! 1) for accelerated melting from collisions, should use rain mass collected by snow, not snow mass
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.
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
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
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)
75 ! CHANGES FOR WRFV3.8.1
76 ! 1) changes and cleanup of code comments
77 ! 2) correction to universal gas constant (very small change)
80 ! 1) fix to saturation vapor pressure polysvp to work at T < -80 C
81 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
84 ! TWG = Timothy Glotfelty, EPA
85 ! Adapted from WRFV3.8.1 Morrison Double Moment Scheme
86 ! Author's of Adapted Version: Timothy Glotfelty and Kiran Alapaty, EPA
87 ! CHANGES FOR AEROSOL Version
89 ! 1) Add Original Aerosol Activation Changes of Hugh Morrison, NCAR and Amy Solomon (NOAA) from WRFV3.1
91 ! 2) Add the Aerosol Activation Scheme of Abdul-Razzak and Ghan (2000) from the
92 ! Morrison-Gettelman Scheme and Song and Zhang (2011) microphysics treatment as an
93 ! aerosol activation option
95 ! 3) Add the Ice Nucleation scheme of Liu and Penner (2005) / Liu et al. (2007)
96 ! from the Morrison-Gettelman Scheme and Song and Zhang (2011) microphysics as an ice
99 ! 4) Introduce prescribed aerosol concentrations from CESM into the above
100 ! mentioned schemes to predict cloud droplet number concentrations and ice
101 ! number concentrations
103 ! 5) Output the cloud liquid, ice, and snow effective radii to be coupled to the
104 ! RRTMG radiation scheme.
106 ! 6) Tune the DCS parameter by increasing it the 350 um in order to get better
107 ! radaition performance after coupling.
108 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
110 ! THIS SCHEME IS A BULK DOUBLE-MOMENT SCHEME THAT PREDICTS MIXING
111 ! RATIOS AND NUMBER CONCENTRATIONS OF FIVE HYDROMETEOR SPECIES:
112 ! CLOUD DROPLETS, CLOUD (SMALL) ICE, RAIN, SNOW, AND GRAUPEL/HAIL.
114 MODULE MODULE_MP_MORR_TWO_MOMENT_AERO
118 ! USE WRF PHYSICS CONSTANTS
119 use module_model_constants, ONLY: CP, G, R => r_d, RV => r_v, EP_2, t0 ! TWG add
121 ! USE module_state_description
125 REAL, PARAMETER :: PI = 3.1415926535897932384626434
126 REAL, PARAMETER :: xxx = 0.9189385332046727417803297
127 REAL, PARAMETER :: real_tiny = TINY(1.0)
129 PUBLIC :: MP_MORR_TWO_MOMENT_AERO
132 PRIVATE :: GAMMA, DERF1
133 PRIVATE :: PI, xxx, real_tiny
134 PRIVATE :: MORR_TWO_MOMENT_MICRO
136 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
137 ! SWITCHES FOR MICROPHYSICS SCHEME
138 ! IACT = 1, USE POWER-LAW CCN SPECTRA, NCCN = CS^K
139 ! IACT = 2, USE LOGNORMAL AEROSOL SIZE DIST TO DERIVE CCN SPECTRA
140 ! IACT = 3, ACTIVATION CALCULATED IN MODULE_MIXACTIVATE
141 ! IACT = 4, TWG 2016 Activation of Prescribed CESM Aerosool
143 INTEGER, PRIVATE :: IACT
145 ! INUM = 0, PREDICT DROPLET CONCENTRATION
146 ! INUM = 1, ASSUME CONSTANT DROPLET CONCENTRATION
147 ! !!!NOTE: PREDICTED DROPLET CONCENTRATION NOT AVAILABLE IN THIS VERSION
148 ! CONTACT HUGH MORRISON (morrison@ucar.edu) FOR FURTHER INFORMATION
150 INTEGER, PRIVATE :: INUM
152 ! FOR INUM = 1, SET CONSTANT DROPLET CONCENTRATION (CM-3)
153 REAL, PRIVATE :: NDCNST
155 ! SWITCH FOR LIQUID-ONLY RUN
156 ! ILIQ = 0, INCLUDE ICE
157 ! ILIQ = 1, LIQUID ONLY, NO ICE
159 INTEGER, PRIVATE :: ILIQ
161 ! SWITCH FOR ICE NUCLEATION
162 ! INUC = 0, USE FORMULA FROM RASMUSSEN ET AL. 2002 (MID-LATITUDE)
163 ! = 1, USE MPACE OBSERVATIONS
164 ! = 2, USE Prescribed Aerosol Ice Nucleation ! TWG 2016 add
166 INTEGER, PRIVATE :: INUC
168 ! IBASE = 1, NEGLECT DROPLET ACTIVATION AT LATERAL CLOUD EDGES DUE TO
169 ! UNRESOLVED ENTRAINMENT AND MIXING, ACTIVATE
170 ! AT CLOUD BASE OR IN REGION WITH LITTLE CLOUD WATER USING
171 ! NON-EQULIBRIUM SUPERSATURATION,
172 ! IN CLOUD INTERIOR ACTIVATE USING EQUILIBRIUM SUPERSATURATION
173 ! IBASE = 2, ASSUME DROPLET ACTIVATION AT LATERAL CLOUD EDGES DUE TO
174 ! UNRESOLVED ENTRAINMENT AND MIXING DOMINATES,
175 ! ACTIVATE DROPLETS EVERYWHERE IN THE CLOUD USING NON-EQUILIBRIUM
176 ! SUPERSATURATION, BASED ON THE
177 ! LOCAL SUB-GRID AND/OR GRID-SCALE VERTICAL VELOCITY
180 ! NOTE: ONLY USED FOR PREDICTED DROPLET CONCENTRATION (INUM = 0) IN NON-WRF-CHEM VERSION OF CODE
182 INTEGER, PRIVATE :: IBASE
184 ! INCLUDE SUB-GRID VERTICAL VELOCITY IN DROPLET ACTIVATION
185 ! ISUB = 0, INCLUDE SUB-GRID W (RECOMMENDED FOR LOWER RESOLUTION)
186 ! ISUB = 1, EXCLUDE SUB-GRID W, ONLY USE GRID-SCALE W
188 ! NOTE: ONLY USED FOR PREDICTED DROPLET CONCENTRATION (INUM = 0) IN NON-WRF-CHEM VERSION OF CODE
190 INTEGER, PRIVATE :: ISUB
192 ! SWITCH FOR GRAUPEL/NO GRAUPEL
193 ! IGRAUP = 0, INCLUDE GRAUPEL
194 ! IGRAUP = 1, NO GRAUPEL
196 INTEGER, PRIVATE :: IGRAUP
198 ! HM ADDED NEW OPTION FOR HAIL
199 ! SWITCH FOR HAIL/GRAUPEL
200 ! IHAIL = 0, DENSE PRECIPITATING ICE IS GRAUPEL
201 ! IHAIL = 1, DENSE PRECIPITATING GICE IS HAIL
203 INTEGER, PRIVATE :: IHAIL
205 ! CLOUD MICROPHYSICS CONSTANTS
207 REAL, PRIVATE :: AI,AC,AS,AR,AG ! 'A' PARAMETER IN FALLSPEED-DIAM RELATIONSHIP
208 REAL, PRIVATE :: BI,BC,BS,BR,BG ! 'B' PARAMETER IN FALLSPEED-DIAM RELATIONSHIP
209 ! REAL, PRIVATE :: R ! GAS CONSTANT FOR AIR
210 ! REAL, PRIVATE :: RV ! GAS CONSTANT FOR WATER VAPOR
211 ! REAL, PRIVATE :: CP ! SPECIFIC HEAT AT CONSTANT PRESSURE FOR DRY AIR
212 REAL, PRIVATE :: RHOSU ! STANDARD AIR DENSITY AT 850 MB
213 REAL, PRIVATE :: RHOW ! DENSITY OF LIQUID WATER
214 REAL, PRIVATE :: RHOI ! BULK DENSITY OF CLOUD ICE
215 REAL, PRIVATE :: RHOSN ! BULK DENSITY OF SNOW
216 REAL, PRIVATE :: RHOG ! BULK DENSITY OF GRAUPEL
217 REAL, PRIVATE :: AIMM ! PARAMETER IN BIGG IMMERSION FREEZING
218 REAL, PRIVATE :: BIMM ! PARAMETER IN BIGG IMMERSION FREEZING
219 REAL, PRIVATE :: ECR ! COLLECTION EFFICIENCY BETWEEN DROPLETS/RAIN AND SNOW/RAIN
220 REAL, PRIVATE :: DCS ! THRESHOLD SIZE FOR CLOUD ICE AUTOCONVERSION
221 REAL, PRIVATE :: MI0 ! INITIAL SIZE OF NUCLEATED CRYSTAL
222 REAL, PRIVATE :: MG0 ! MASS OF EMBRYO GRAUPEL
223 REAL, PRIVATE :: F1S ! VENTILATION PARAMETER FOR SNOW
224 REAL, PRIVATE :: F2S ! VENTILATION PARAMETER FOR SNOW
225 REAL, PRIVATE :: F1R ! VENTILATION PARAMETER FOR RAIN
226 REAL, PRIVATE :: F2R ! VENTILATION PARAMETER FOR RAIN
227 ! REAL, PRIVATE :: G ! GRAVITATIONAL ACCELERATION
228 REAL, PRIVATE :: QSMALL ! SMALLEST ALLOWED HYDROMETEOR MIXING RATIO
229 REAL, PRIVATE :: CI,DI,CS,DS,CG,DG ! SIZE DISTRIBUTION PARAMETERS FOR CLOUD ICE, SNOW, GRAUPEL
230 REAL, PRIVATE :: EII ! COLLECTION EFFICIENCY, ICE-ICE COLLISIONS
231 REAL, PRIVATE :: ECI ! COLLECTION EFFICIENCY, ICE-DROPLET COLLISIONS
232 REAL, PRIVATE :: RIN ! RADIUS OF CONTACT NUCLEI (M)
234 REAL, PRIVATE :: CPW ! SPECIFIC HEAT OF LIQUID WATER
236 ! CCN SPECTRA FOR IACT = 1
238 REAL, PRIVATE :: C1 ! 'C' IN NCCN = CS^K (CM-3)
239 REAL, PRIVATE :: K1 ! 'K' IN NCCN = CS^K
241 ! AEROSOL PARAMETERS FOR IACT = 2
243 REAL, PRIVATE :: MW ! MOLECULAR WEIGHT WATER (KG/MOL)
244 REAL, PRIVATE :: OSM ! OSMOTIC COEFFICIENT
245 REAL, PRIVATE :: VI ! NUMBER OF ION DISSOCIATED IN SOLUTION
246 REAL, PRIVATE :: EPSM ! AEROSOL SOLUBLE FRACTION
247 REAL, PRIVATE :: RHOA ! AEROSOL BULK DENSITY (KG/M3)
248 REAL, PRIVATE :: MAP ! MOLECULAR WEIGHT AEROSOL (KG/MOL)
249 REAL, PRIVATE :: MA ! MOLECULAR WEIGHT OF 'AIR' (KG/MOL)
250 REAL, PRIVATE :: RR ! UNIVERSAL GAS CONSTANT
251 REAL, PRIVATE :: BACT ! ACTIVATION PARAMETER
252 REAL, PRIVATE :: RM1 ! GEOMETRIC MEAN RADIUS, MODE 1 (M)
253 REAL, PRIVATE :: RM2 ! GEOMETRIC MEAN RADIUS, MODE 2 (M)
254 REAL, PRIVATE :: NANEW1 ! TOTAL AEROSOL CONCENTRATION, MODE 1 (M^-3)
255 REAL, PRIVATE :: NANEW2 ! TOTAL AEROSOL CONCENTRATION, MODE 2 (M^-3)
256 REAL, PRIVATE :: SIG1 ! STANDARD DEVIATION OF AEROSOL S.D., MODE 1
257 REAL, PRIVATE :: SIG2 ! STANDARD DEVIATION OF AEROSOL S.D., MODE 2
258 REAL, PRIVATE :: F11 ! CORRECTION FACTOR FOR ACTIVATION, MODE 1
259 REAL, PRIVATE :: F12 ! CORRECTION FACTOR FOR ACTIVATION, MODE 1
260 REAL, PRIVATE :: F21 ! CORRECTION FACTOR FOR ACTIVATION, MODE 2
261 REAL, PRIVATE :: F22 ! CORRECTION FACTOR FOR ACTIVATION, MODE 2
262 REAL, PRIVATE :: MMULT ! MASS OF SPLINTERED ICE PARTICLE
263 REAL, PRIVATE :: LAMMAXI,LAMMINI,LAMMAXR,LAMMINR,LAMMAXS,LAMMINS,LAMMAXG,LAMMING
265 ! CONSTANTS TO IMPROVE EFFICIENCY
267 REAL, PRIVATE :: CONS1,CONS2,CONS3,CONS4,CONS5,CONS6,CONS7,CONS8,CONS9,CONS10
268 REAL, PRIVATE :: CONS11,CONS12,CONS13,CONS14,CONS15,CONS16,CONS17,CONS18,CONS19,CONS20
269 REAL, PRIVATE :: CONS21,CONS22,CONS23,CONS24,CONS25,CONS26,CONS27,CONS28,CONS29,CONS30
270 REAL, PRIVATE :: CONS31,CONS32,CONS33,CONS34,CONS35,CONS36,CONS37,CONS38,CONS39,CONS40
271 REAL, PRIVATE :: CONS41
273 !TWG 2016 Begin Add the following for prescribed aerosols
274 integer, parameter :: naer_cu = 10 !xsong 2013-08-22 !
275 real, private:: aten_pamdm
276 real, private:: alogsig_pamdm(naer_cu) ! natl log of geometric standard dev of aerosol
277 real, private:: exp45logsig_pamdm(naer_cu)
278 real, private:: argfactor_pamdm(naer_cu)
279 real, private:: amcube_pamdm(naer_cu) ! cube of dry mode radius (m)
280 real, private:: smcrit_pamdm(naer_cu) ! critical supersatuation for activation
281 real, private:: lnsm_pamdm(naer_cu) ! ln(smcrit_pamdm)
282 real, private:: amcubefactor_pamdm(naer_cu) ! factors for calculating mode radius
283 real, private:: smcritfactor_pamdm(naer_cu) ! factors for calculatingcritical supersaturation
284 real, private:: alog2_pamdm,alog3_pamdm,alogaten_pamdm,surften_pamdm
286 real, private:: f1_pamdm(naer_cu),f2_pamdm(naer_cu) ! abdul-razzak functions of width
287 real, private:: third_pamdm
288 integer, private:: sat_pamdm
289 real, private:: sq2_pamdm, arg_pamdm
290 integer, parameter:: psat_pamdm=7 ! number of supersaturations to calc ccn concentration
291 real, private:: super_pamdm(psat_pamdm)
292 real, private, parameter :: supersat_pamdm(psat_pamdm)= &! supersaturation (%) to determine ccn concentration
293 (/0.02,0.05,0.1,0.2,0.3,0.5,1.0/)
294 real, private:: ccnfact_pamdm(psat_pamdm,naer_cu) !factor to calculate diagnostic CCN
298 !xsong 2013-08-22---------------
300 character(len=20) aername(naer_cu)
301 REAL dryrad_aer(naer_cu)
302 REAL density_aer(naer_cu)
303 REAL hygro_aer(naer_cu)
304 REAL dispersion_aer(naer_cu)
305 REAL num_to_mass_aer(naer_cu)
307 !xsong 2013-08-22--------------------
308 data aername /"SULFATE","SEASALT2","DUST1","DUST2","DUST3","DUST4","OCPHO","BCPHO", &
310 data dryrad_aer /0.695E-7,0.200E-5,0.151E-5,0.151E-5,0.151E-5,0.151E-5, &
311 0.212E-7,0.118E-7,0.212E-7, 0.118E-7/
312 data density_aer /1770.,2200.,2600.,2600.,2600.,2600.,1800., &
314 data hygro_aer /0.507,1.160,0.140,0.140,0.140,0.140,0.100,0.100, &
316 data dispersion_aer /2.030,1.3732,1.900,1.900,1.900,1.900,2.240, &
318 data num_to_mass_aer /42097098109277080.,8626504211623.,3484000000000000.,213800000000000.,&
319 22050000000000.,3165000000000.,0.745645E+18,0.167226E+20,&
320 0.516216E+18,0.167226E+20/
321 !xsong 2013-08-22-----------------------
328 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
329 SUBROUTINE MORR_TWO_MOMENT_INIT_AERO(morr_rimed_ice) ! RAS
330 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
331 ! THIS SUBROUTINE INITIALIZES ALL PHYSICAL CONSTANTS AMND PARAMETERS
332 ! NEEDED BY THE MICROPHYSICS SCHEME.
333 ! NEEDS TO BE CALLED AT FIRST TIME STEP, PRIOR TO CALL TO MAIN MICROPHYSICS INTERFACE
334 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
337 SAVE !TWG 2017 add to avoid memory issues
339 INTEGER, INTENT(IN):: morr_rimed_ice ! RAS
341 integer n,i,m !TWG add
343 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
344 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
346 ! THE FOLLOWING PARAMETERS ARE USER-DEFINED SWITCHES AND NEED TO BE
347 ! SET PRIOR TO CODE COMPILATION
349 ! INUM IS AUTOMATICALLY SET TO 0 FOR WRF-CHEM BELOW,
350 ! ALLOWING PREDICTION OF DROPLET CONCENTRATION
351 ! THUS, THIS PARAMETER SHOULD NOT BE CHANGED HERE
352 ! AND SHOULD BE LEFT TO 1
356 ! SET CONSTANT DROPLET CONCENTRATION (UNITS OF CM-3)
357 ! IF NO COUPLING WITH WRF-CHEM
361 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
362 ! NOTE, THE FOLLOWING OPTIONS RELATED TO DROPLET ACTIVATION
363 ! (IACT, IBASE, ISUB) ARE NOT AVAILABLE IN CURRENT VERSION
364 ! FOR WRF-CHEM, DROPLET ACTIVATION IS PERFORMED
365 ! IN 'MIX_ACTIVATE', NOT IN MICROPHYSICS SCHEME
368 ! IACT = 1, USE POWER-LAW CCN SPECTRA, NCCN = CS^K
369 ! IACT = 2, USE LOGNORMAL AEROSOL SIZE DIST TO DERIVE CCN SPECTRA
373 ! IBASE = 1, NEGLECT DROPLET ACTIVATION AT LATERAL CLOUD EDGES DUE TO
374 ! UNRESOLVED ENTRAINMENT AND MIXING, ACTIVATE
375 ! AT CLOUD BASE OR IN REGION WITH LITTLE CLOUD WATER USING
376 ! NON-EQULIBRIUM SUPERSATURATION ASSUMING NO INITIAL CLOUD WATER,
377 ! IN CLOUD INTERIOR ACTIVATE USING EQUILIBRIUM SUPERSATURATION
378 ! IBASE = 2, ASSUME DROPLET ACTIVATION AT LATERAL CLOUD EDGES DUE TO
379 ! UNRESOLVED ENTRAINMENT AND MIXING DOMINATES,
380 ! ACTIVATE DROPLETS EVERYWHERE IN THE CLOUD USING NON-EQUILIBRIUM
381 ! SUPERSATURATION ASSUMING NO INITIAL CLOUD WATER, BASED ON THE
382 ! LOCAL SUB-GRID AND/OR GRID-SCALE VERTICAL VELOCITY
385 ! NOTE: ONLY USED FOR PREDICTED DROPLET CONCENTRATION (INUM = 0)
389 ! INCLUDE SUB-GRID VERTICAL VELOCITY (standard deviation of w) IN DROPLET ACTIVATION
390 ! ISUB = 0, INCLUDE SUB-GRID W (RECOMMENDED FOR LOWER RESOLUTION)
391 ! currently, sub-grid w is constant of 0.5 m/s (not coupled with PBL/turbulence scheme)
392 ! ISUB = 1, EXCLUDE SUB-GRID W, ONLY USE GRID-SCALE W
394 ! NOTE: ONLY USED FOR PREDICTED DROPLET CONCENTRATION (INUM = 0)
397 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
400 ! SWITCH FOR LIQUID-ONLY RUN
401 ! ILIQ = 0, INCLUDE ICE
402 ! ILIQ = 1, LIQUID ONLY, NO ICE
406 ! SWITCH FOR ICE NUCLEATION
407 ! INUC = 0, USE FORMULA FROM RASMUSSEN ET AL. 2002 (MID-LATITUDE)
408 ! = 1, USE MPACE OBSERVATIONS (ARCTIC ONLY)
412 ! SWITCH FOR GRAUPEL/HAIL NO GRAUPEL/HAIL
413 ! IGRAUP = 0, INCLUDE GRAUPEL/HAIL
414 ! IGRAUP = 1, NO GRAUPEL/HAIL
419 ! SWITCH FOR HAIL/GRAUPEL
420 ! IHAIL = 0, DENSE PRECIPITATING ICE IS GRAUPEL
421 ! IHAIL = 1, DENSE PRECIPITATING ICE IS HAIL
422 ! NOTE ---> RECOMMEND IHAIL = 1 FOR CONTINENTAL DEEP CONVECTION
424 !IHAIL = 0 !changed to namelist option (morr_rimed_ice) by RAS
425 ! Check if namelist option is feasible, otherwise default to graupel - RAS
426 IF (morr_rimed_ice .eq. 1) THEN
432 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
433 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
434 ! SET PHYSICAL CONSTANTS
436 ! FALLSPEED PARAMETERS (V=AD^B)
448 ELSE ! (MATSUN AND HUGGINS 1980)
453 ! CONSTANTS AND PARAMETERS
457 RHOSU = 85000./(287.15*273.15)
469 DCS = 350.E-6 !TWG Increase from 150.E-6
470 MI0 = 4./3.*PI*RHOI*(10.E-6)**3
487 ! SIZE DISTRIBUTION PARAMETERS
496 ! RADIUS OF CONTACT NUCLEI
499 MMULT = 4./3.*PI*RHOI*(5.E-6)**3
501 ! SIZE LIMITS FOR LAMBDA
502 ! TWG Feb 2017 Effective radii Test
504 LAMMINI = 1./(2.*DCS+100.E-6)
506 ! LAMMINI = 1./120.E-6
508 ! LAMMINR = 1./500.E-6
509 LAMMINR = 1./2800.E-6
512 LAMMINS = 1./2000.E-6
513 ! LAMMINS = 1./1000.E-6
515 LAMMING = 1./2000.E-6
517 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
518 ! note: these parameters only used by the non-wrf-chem version of the
519 ! scheme with predicted droplet number
521 ! CCN SPECTRA FOR IACT = 1
524 ! MODIFIED FROM RASMUSSEN ET AL. 2002
525 ! NCCN = C*S^K, NCCN IS IN CM-3, S IS SUPERSATURATION RATIO IN %
535 ! AEROSOL ACTIVATION PARAMETERS FOR IACT = 2
536 ! PARAMETERS CURRENTLY SET FOR AMMONIUM SULFATE
548 BACT = VI*OSM*EPSM*MW*RHOA/(MAP*RHOW)
550 ! AEROSOL SIZE DISTRIBUTION PARAMETERS CURRENTLY SET FOR MPACE
551 ! (see morrison et al. 2007, JGR)
557 F11 = 0.5*EXP(2.5*(LOG(SIG1))**2)
558 F21 = 1.+0.25*LOG(SIG1)
565 F12 = 0.5*EXP(2.5*(LOG(SIG2))**2)
566 F22 = 1.+0.25*LOG(SIG2)
567 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
569 ! CONSTANTS FOR EFFICIENCY
571 CONS1=GAMMA(1.+DS)*CS
572 CONS2=GAMMA(1.+DG)*CG
573 CONS3=GAMMA(4.+BS)/6.
574 CONS4=GAMMA(4.+BR)/6.
577 CONS7=GAMMA(4.+BG)/6.
579 CONS9=GAMMA(5./2.+BR/2.)
580 CONS10=GAMMA(5./2.+BS/2.)
581 CONS11=GAMMA(5./2.+BG/2.)
582 CONS12=GAMMA(1.+DI)*CI
583 CONS13=GAMMA(BS+3.)*PI/4.*ECI
584 CONS14=GAMMA(BG+3.)*PI/4.*ECI
585 CONS15=-1108.*EII*PI**((1.-BS)/3.)*RHOSN**((-2.-BS)/3.)/(4.*720.)
586 CONS16=GAMMA(BI+3.)*PI/4.*ECI
587 CONS17=4.*2.*3.*RHOSU*PI*ECI*ECI*GAMMA(2.*BS+2.)/(8.*(RHOG-RHOSN))
590 CONS20=20.*PI*PI*RHOW*BIMM
592 CONS22=PI*RHOI*DCS**3/6.
593 CONS23=PI/4.*EII*GAMMA(BS+3.)
594 CONS24=PI/4.*ECR*GAMMA(BR+3.)
595 CONS25=PI*PI/24.*RHOW*ECR*GAMMA(BR+6.)
598 CONS28=GAMMA(4.+BI)/6.
599 CONS29=4./3.*PI*RHOW*(25.E-6)**3
601 CONS31=PI*PI*ECR*RHOSN
603 CONS33=PI*PI*ECR*RHOG
607 CONS37=4.*PI*1.38E-23/(6.*PI*RIN)
609 CONS39=PI*PI/36.*RHOW*BIMM
611 CONS41=PI*PI*ECR*RHOW
613 !+---+-----------------------------------------------------------------+
614 !..Set these variables needed for computing radar reflectivity. These
615 !.. get used within radar_init to create other variables used in the
629 !+---+-----------------------------------------------------------------+
631 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
633 ! set parameters for prescribed CESM droplet activation, following abdul-razzak
637 ! mathematical constants
640 sq2_pamdm=1.41421356237
643 aten_pamdm=2.*MW*surften_pamdm/(RR*t0*RHOW)
644 alogaten_pamdm=log(aten_pamdm)
649 ! use only if width of size distribution is prescribed
650 alogsig_pamdm(m)=log(dispersion_aer(m))
651 exp45logsig_pamdm(m)=exp(4.5*alogsig_pamdm(m)*alogsig_pamdm(m))
652 argfactor_pamdm(m)=2./(3.*SQRT(2.)*alogsig_pamdm(m))
653 f1_pamdm(m)=0.5*exp(2.5*alogsig_pamdm(m)*alogsig_pamdm(m))
654 f2_pamdm(m)=1.+0.25*alogsig_pamdm(m)
655 amcubefactor_pamdm(m)=3./(4.*pi*exp45logsig_pamdm(m)*density_aer(m))
656 smcritfactor_pamdm(m)=2.*aten_pamdm*SQRT(aten_pamdm/(27.*max(1.e-10,hygro_aer(m))))
657 ! use only if mode radius of size distribution is prescribed
658 amcube_pamdm(m)=amcubefactor_pamdm(m)/(num_to_mass_aer(m))
659 ! use only if only one component per mode
660 if(hygro_aer(m).gt.1.e-10) then
661 smcrit_pamdm(m)=smcritfactor_pamdm(m)/SQRT(amcube_pamdm(m))
665 lnsm_pamdm(m)=log(smcrit_pamdm(m))
669 do sat_pamdm=1,psat_pamdm
670 super_pamdm(sat_pamdm) = 0.01*supersat_pamdm(sat_pamdm)
671 arg_pamdm=argfactor_pamdm(m)*log(smcrit_pamdm(m)/super_pamdm(sat_pamdm))
672 if(arg_pamdm.lt.2) then
673 if(arg_pamdm.lt.-2) then
674 ccnfact_pamdm(sat_pamdm,m)=1.e-6
676 ccnfact_pamdm(sat_pamdm,m)=1.e-6*0.5*(1-DERF1(arg_pamdm))
679 ccnfact_pamdm(sat_pamdm,m)=0.0
681 ! print*,'m',m,' sat ',sat_pamdm,' ccnfact',ccnfact_pamdm(sat_pamdm,m)
688 END SUBROUTINE MORR_TWO_MOMENT_INIT_AERO
690 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
691 ! THIS SUBROUTINE IS MAIN INTERFACE WITH THE TWO-MOMENT MICROPHYSICS SCHEME
692 ! THIS INTERFACE TAKES IN 3D VARIABLES FROM DRIVER MODEL, CONVERTS TO 1D FOR
693 ! CALL TO THE MAIN MICROPHYSICS SUBROUTINE (SUBROUTINE MORR_TWO_MOMENT_MICRO)
694 ! WHICH OPERATES ON 1D VERTICAL COLUMNS.
695 ! 1D VARIABLES FROM THE MAIN MICROPHYSICS SUBROUTINE ARE THEN REASSIGNED BACK TO 3D FOR OUTPUT
696 ! BACK TO DRIVER MODEL USING THIS INTERFACE.
697 ! MICROPHYSICS TENDENCIES ARE ADDED TO VARIABLES HERE BEFORE BEING PASSED BACK TO DRIVER MODEL.
699 ! THIS CODE WAS WRITTEN BY HUGH MORRISON (NCAR) AND SLAVA TATARSKII (GEORGIA TECH).
701 ! FOR QUESTIONS, CONTACT: HUGH MORRISON, E-MAIL: MORRISON@UCAR.EDU, PHONE:303-497-8916
703 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
705 SUBROUTINE MP_MORR_TWO_MOMENT_AERO(ITIMESTEP, &
706 TH, QV, QC, QR, QI, QS, QG, NI, NS, NR, NG, NC, & !TWG/amy add
707 RHO, PII, P, DT_IN, DZ, HT, W, KZH, &
708 RAINNC, RAINNCV, SR, &
709 SNOWNC,SNOWNCV,GRAUPELNC,GRAUPELNCV, & ! hm added 7/13/13
710 refl_10cm,mskf_refl_10cm,diagflag, do_radar_ref, & ! GT added for reflectivity calcs
711 qrcuten, qscuten, qicuten & ! hm added
712 ,F_QNDROP, qndrop & ! hm added, wrf-chem
713 ,IDS,IDE, JDS,JDE, KDS,KDE & ! domain dims
714 ,IMS,IME, JMS,JME, KMS,KME & ! memory dims
715 ,ITS,ITE, JTS,JTE, KTS,KTE & ! tile dimsi )
718 ,aercu_opt & ! TWG add
719 ,aercu_fct & ! TWG add
720 ,no_src_types_cu & ! TWG add
721 ,EFCG, EFIG, EFSG, WACT & !TWG add
722 ,CCN1_GS, CCN2_GS, CCN3_GS, CCN4_GS & !TWG add
723 ,CCN5_GS, CCN6_GS, CCN7_GS, NR_CU & ! TWG add
724 ,QR_CU, NS_CU, QS_CU, CU_UAF & ! TWG add
725 !jdf ,C2PREC3D,CSED3D,ISED3D,SSED3D,GSED3D,RSED3D & ! HM ADD, WRF-CHEM
726 ,wetscav_on, rainprod, evapprod &
727 ,QLSINK,PRECR,PRECI,PRECS,PRECG & ! HM ADD, WRF-CHEM
730 ! QV - water vapor mixing ratio (kg/kg)
731 ! QC - cloud water mixing ratio (kg/kg)
732 ! QR - rain water mixing ratio (kg/kg)
733 ! QI - cloud ice mixing ratio (kg/kg)
734 ! QS - snow mixing ratio (kg/kg)
735 ! QG - graupel mixing ratio (KG/KG)
736 ! NI - cloud ice number concentration (1/kg)
737 ! NS - Snow Number concentration (1/kg)
738 ! NR - Rain Number concentration (1/kg)
739 ! NG - Graupel number concentration (1/kg)
740 ! NOTE: RHO AND HT NOT USED BY THIS SCHEME AND DO NOT NEED TO BE PASSED INTO SCHEME!!!!
741 ! P - AIR PRESSURE (PA)
742 ! W - VERTICAL AIR VELOCITY (M/S)
743 ! TH - POTENTIAL TEMPERATURE (K)
744 ! PII - exner function - used to convert potential temp to temp
745 ! DZ - difference in height over interface (m)
746 ! DT_IN - model time step (sec)
747 ! ITIMESTEP - time step counter
748 ! RAINNC - accumulated grid-scale precipitation (mm)
749 ! RAINNCV - one time step grid scale precipitation (mm/time step)
750 ! SNOWNC - accumulated grid-scale snow plus cloud ice (mm)
751 ! SNOWNCV - one time step grid scale snow plus cloud ice (mm/time step)
752 ! GRAUPELNC - accumulated grid-scale graupel (mm)
753 ! GRAUPELNCV - one time step grid scale graupel (mm/time step)
754 ! SR - one time step mass ratio of snow to total precip
755 ! qrcuten, rain tendency from parameterized cumulus convection
756 ! qscuten, snow tendency from parameterized cumulus convection
757 ! qicuten, cloud ice tendency from parameterized cumulus convection
759 ! variables below currently not in use, not coupled to PBL or radiation codes
760 ! KZH - heat eddy diffusion coefficient from YSU scheme (exch_h, M^2 S-1), or other PBL that provides this
761 ! NEEDED FOR DROPLET ACTIVATION (SEE CODE BELOW)
762 ! EFFCS - CLOUD DROPLET EFFECTIVE RADIUS OUTPUT TO RADIATION CODE (micron)
763 ! EFFIS - CLOUD DROPLET EFFECTIVE RADIUS OUTPUT TO RADIATION CODE (micron)
764 ! HM, ADDED FOR WRF-CHEM COUPLING
765 ! QLSINK - TENDENCY OF CLOUD WATER TO RAIN, SNOW, GRAUPEL (KG/KG/S)
766 ! CSED,ISED,SSED,GSED,RSED - SEDIMENTATION FLUXES (KG/M^2/S) FOR CLOUD WATER, ICE, SNOW, GRAUPEL, RAIN
767 ! PRECI,PRECS,PRECG,PRECR - SEDIMENTATION FLUXES (KG/M^2/S) FOR ICE, SNOW, GRAUPEL, RAIN
769 ! rainprod - total tendency of conversion of cloud water/ice and graupel to rain (kg kg-1 s-1)
770 ! evapprod - tendency of evaporation of rain (kg kg-1 s-1)
772 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
774 ! reflectivity currently not included!!!!
775 ! REFL_10CM - CALCULATED RADAR REFLECTIVITY AT 10 CM (DBZ)
776 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
778 ! EFFC - DROPLET EFFECTIVE RADIUS (MICRON)
779 ! EFFR - RAIN EFFECTIVE RADIUS (MICRON)
780 ! EFFS - SNOW EFFECTIVE RADIUS (MICRON)
781 ! EFFI - CLOUD ICE EFFECTIVE RADIUS (MICRON)
783 ! ADDITIONAL OUTPUT FROM MICRO - SEDIMENTATION TENDENCIES, NEEDED FOR LIQUID-ICE STATIC ENERGY
785 ! QGSTEN - GRAUPEL SEDIMENTATION TEND (KG/KG/S)
786 ! QRSTEN - RAIN SEDIMENTATION TEND (KG/KG/S)
787 ! QISTEN - CLOUD ICE SEDIMENTATION TEND (KG/KG/S)
788 ! QNISTEN - SNOW SEDIMENTATION TEND (KG/KG/S)
789 ! QCSTEN - CLOUD WATER SEDIMENTATION TEND (KG/KG/S)
791 ! WVAR - STANDARD DEVIATION OF SUB-GRID VERTICAL VELOCITY (M/S)
794 SAVE !TWG 2017 add to avoid memeory issues
796 INTEGER, INTENT(IN ) :: ids, ide, jds, jde, kds, kde , &
797 ims, ime, jms, jme, kms, kme , &
798 its, ite, jts, jte, kts, kte
799 INTEGER, INTENT(IN ) :: PBL
800 ! Temporary changed from INOUT to IN
802 REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT):: &
803 qv, qc, qr, qi, qs, qg, ni, ns, nr, TH, NG, NC !TWG add
804 !jdf qndrop ! hm added, wrf-chem
805 REAL, DIMENSION(ims:ime, kms:kme, jms:jme), optional,INTENT(INOUT):: qndrop
806 !jdf REAL, DIMENSION(ims:ime, kms:kme, jms:jme),INTENT(INOUT):: CSED3D, &
807 REAL, DIMENSION(ims:ime, kms:kme, jms:jme), optional,INTENT(INOUT):: QLSINK, &
808 rainprod, evapprod, &
809 PRECI,PRECS,PRECG,PRECR ! HM, WRF-CHEM
812 REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(IN):: &
813 pii, p, dz, rho, w, kzh ! TWG/amy add
814 REAL, INTENT(IN):: dt_in
815 INTEGER, INTENT(IN):: ITIMESTEP
817 REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT):: &
818 RAINNC, RAINNCV, SR, &
820 SNOWNC,SNOWNCV,GRAUPELNC,GRAUPELNCV
822 REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT):: & ! GT
825 REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT):: & ! GT
828 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN) :: ht, CU_UAF !TWG add
830 REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(IN):: &
831 NR_CU,QR_CU,NS_CU,QS_CU
833 REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT):: &
834 EFCG, EFIG, EFSG, WACT, CCN1_GS, &
835 CCN2_GS, CCN3_GS, CCN4_GS, CCN5_GS, CCN6_GS, CCN7_GS
837 LOGICAL, optional, INTENT(IN) :: wetscav_on
841 REAL, DIMENSION(its:ite, kts:kte, jts:jte):: &
842 effi, effs, effr, EFFG
844 REAL, DIMENSION(its:ite, kts:kte, jts:jte):: &
847 REAL, DIMENSION(kts:kte) :: &
848 QC_TEND1D, QI_TEND1D, QNI_TEND1D, QR_TEND1D, &
849 NI_TEND1D, NS_TEND1D, NR_TEND1D, &
850 QC1D, QI1D, QR1D,NI1D, NS1D, NR1D, QS1D, &
851 T_TEND1D,QV_TEND1D, T1D, QV1D, P1D, W1D, WVAR1D, &
852 EFFC1D, EFFI1D, EFFS1D, EFFR1D,DZ1D, &
854 QG_TEND1D, NG_TEND1D, QG1D, NG1D, EFFG1D, &
856 ! ADD SEDIMENTATION TENDENCIES (UNITS OF KG/KG/S)
857 QGSTEN,QRSTEN, QISTEN, QNISTEN, QCSTEN, &
858 ! ADD CUMULUS TENDENCIES
859 QRCU1D, QSCU1D, QICU1D, &
861 NRMSCU1D, NSMSCU1D, QRMSCU1D, QSMSCU1D, & ! TWG add
863 NRRAD, QRRAD, NSRAD, QSRAD
865 ! add cumulus tendencies
867 REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(IN):: &
868 qrcuten, qscuten, qicuten
870 LOGICAL, INTENT(IN), OPTIONAL :: F_QNDROP ! wrf-chem
871 LOGICAL :: flag_qndrop ! wrf-chem
872 integer :: iinum ! wrf-chem
875 REAL, DIMENSION(kts:kte) :: nc1d, nc_tend1d,C2PREC,CSED,ISED,SSED,GSED,RSED
876 REAL, DIMENSION(kts:kte) :: rainprod1d, evapprod1d
877 ! HM add reflectivity
878 REAL, DIMENSION(kts:kte) :: dBZ,dBZ2 !TWG 2017 add
880 REAL PRECPRT1D, SNOWRT1D, SNOWPRT1D, GRPLPRT1D ! hm added 7/13/13
882 INTEGER I,K,J,L !TWG Add
886 LOGICAL, OPTIONAL, INTENT(IN) :: diagflag
887 INTEGER, OPTIONAL, INTENT(IN) :: do_radar_ref
889 ! TWG begin declare prescribed aerosol parameters
890 INTEGER, INTENT(IN ) :: no_src_types_cu !PSH/TWG
891 INTEGER, INTENT(IN ) :: aercu_opt
892 REAL, INTENT(IN ) :: aercu_fct
893 REAL, DIMENSION( ims:ime, kms:kme, jms:jme, no_src_types_cu), OPTIONAL, &
894 INTENT(INOUT) :: aerocu !PSH/TWG
896 REAL, DIMENSION(kts:kte, no_src_types_cu) :: maerosol, naer
899 LOGICAL :: has_wetscav
902 flag_qndrop = .false.
903 IF ( PRESENT ( f_qndrop ) ) flag_qndrop = f_qndrop
904 !!!!!!!!!!!!!!!!!!!!!!
906 IF ( PRESENT ( wetscav_on ) ) THEN
907 has_wetscav = wetscav_on
909 has_wetscav = .false.
912 ! Initialize tendencies (all set to 0) and transfer
913 ! array to local variables
919 T(I,K,J) = TH(i,k,j)*PII(i,k,j)
921 ! NOTE: WVAR NOT CURRENTLY USED IN CODE !!!!!!!!!!
922 ! currently assign wvar to 0.5 m/s (not coupled with PBL scheme)
924 ! wvar is the ST. DEV. OF sub-grid vertical velocity, used for calculating
927 ! WVAR COMES FROM EDDY DIFFUSION COEFFICIENT KZH (e.g. FROM YSU PBL SCHEME)
928 ! NOTE: IF MODEL HAS HIGH ENOUGH RESOLUTION TO RESOLVE UPDRAFTS, WVAR MAY
931 WVAR(I,K,J) = KZH(I,K+1,J)/20.
932 WVAR(I,K,J) = MAX(0.10,WVAR(I,K,J))
933 WVAR(I,K,J) = MIN(50.,WVAR(I,K,J))
935 ! currently mixing of number concentrations also is neglected (not coupled with PBL schemes)
941 do i=its,ite ! i loop (east-west)
942 do j=jts,jte ! j loop (north-south)
944 ! Transfer 3D arrays into 1D for microphysical calculations
947 ! hm , initialize 1d tendency arrays to zero
949 do k=kts,kte ! k loop (vertical)
960 nc_tend1d(k) = 0. ! wrf-chem
962 !TWG add Initalize CCN
970 NRMSCU1D(k) = NR_CU(i,k,j)
971 NSMSCU1D(k) = NS_CU(i,k,j)
972 QRMSCU1D(k) = QR_CU(i,k,j)
973 QSMSCU1D(k) = QS_CU(i,k,j)
1000 WVAR1D(k) = WVAR(i,k,j)
1001 ! add cumulus tendencies
1002 qrcu1d(k) = qrcuten(i,k,j)
1003 qscu1d(k) = qscuten(i,k,j)
1004 qicu1d(k) = qicuten(i,k,j)
1005 end do !jdf added this
1006 ! below for wrf-chem
1007 IF (flag_qndrop .AND. PRESENT( qndrop )) THEN
1010 nc1d(k)=qndrop(i,k,j)
1015 ! TWG 2016 comment this out temporarily
1016 ! nc1d(k)=0. ! temporary placeholder, set to constant in microphysics subroutine
1018 IF (aercu_opt.eq.2) THEN
1024 DO L=1,no_src_types_cu
1029 ! print*,'Doing TWG added activation'
1031 maerosol(k,1) = aercu_fct*aerocu(i,k,j,6)*1E-9
1032 naer(k,1) = 5.64259e13 * maerosol(k,1)**0.58
1034 maerosol(k,2) = aercu_fct*aerocu(i,k,j,5)*1E-9
1035 naer(k,2) = maerosol(k,2)*num_to_mass_aer(2)
1037 maerosol(k,3) = aercu_fct*1.44*aerocu(i,k,j,1)*1E-9
1038 naer(k,3) = maerosol(k,3)*num_to_mass_aer(3)
1040 maerosol(k,4) = aercu_fct*1.44*aerocu(i,k,j,2)*1E-9
1041 naer(k,4) = maerosol(k,4)*num_to_mass_aer(4)
1043 maerosol(k,5) = aercu_fct*1.44*aerocu(i,k,j,3)*1E-9
1044 naer(k,5) = maerosol(k,5)*num_to_mass_aer(5)
1046 maerosol(k,6) = aercu_fct*1.44*aerocu(i,k,j,4)*1E-9
1047 naer(k,6) = maerosol(k,6)*num_to_mass_aer(6)
1049 maerosol(k,7) = aercu_fct*1.54*aerocu(i,k,j,9)*1E-9
1050 naer(k,7) = maerosol(k,7)*num_to_mass_aer(7)
1052 maerosol(k,8) = aercu_fct*1.37*aerocu(i,k,j,7)*1E-9
1053 naer(k,8) = maerosol(k,8)*num_to_mass_aer(8)
1055 maerosol(k,9) = aercu_fct*1.25*aerocu(i,k,j,10)*1E-9
1056 naer(k,9) = maerosol(k,9)*num_to_mass_aer(9)
1058 maerosol(k,10) = aercu_fct*1.37*aerocu(i,k,j,8)*1E-9
1059 naer(k,10) = maerosol(k,10)*num_to_mass_aer(10)
1061 DO L=1,no_src_types_cu
1062 CCN1_GS(i,k,j) = CCN1_GS(i,k,j)+(naer(k,L)*ccnfact_pamdm(1,L))
1063 CCN2_GS(i,k,j) = CCN2_GS(i,k,j)+(naer(k,L)*ccnfact_pamdm(2,L))
1064 CCN3_GS(i,k,j) = CCN3_GS(i,k,j)+(naer(k,L)*ccnfact_pamdm(3,L))
1065 CCN4_GS(i,k,j) = CCN4_GS(i,k,j)+(naer(k,L)*ccnfact_pamdm(4,L))
1066 CCN5_GS(i,k,j) = CCN5_GS(i,k,j)+(naer(k,L)*ccnfact_pamdm(5,L))
1067 CCN6_GS(i,k,j) = CCN6_GS(i,k,j)+(naer(k,L)*ccnfact_pamdm(6,L))
1068 CCN7_GS(i,k,j) = CCN7_GS(i,k,j)+(naer(k,L)*ccnfact_pamdm(7,L))
1072 nc1d(k)=0. ! temporary placeholder, set to constant in microphysics
1082 call MORR_TWO_MOMENT_MICRO(QC_TEND1D, QI_TEND1D, QNI_TEND1D, QR_TEND1D, &
1083 NI_TEND1D, NS_TEND1D, NR_TEND1D, &
1084 QC1D, QI1D, QS1D, QR1D,NI1D, NS1D, NR1D, &
1085 T_TEND1D,QV_TEND1D, T1D, QV1D, P1D, DZ1D, W1D, WVAR1D, &
1086 PRECPRT1D,SNOWRT1D, &
1087 SNOWPRT1D,GRPLPRT1D, & ! hm added 7/13/13
1088 EFFC1D,EFFI1D,EFFS1D,EFFR1D,DT, &
1089 IMS,IME, JMS,JME, KMS,KME, &
1090 ITS,ITE, JTS,JTE, KTS,KTE, & ! HM ADD GRAUPEL
1091 QG_TEND1D,NG_TEND1D,QG1D,NG1D,EFFG1D, &
1092 qrcu1d, qscu1d, qicu1d, &
1093 ! ADD SEDIMENTATION TENDENCIES
1094 QGSTEN,QRSTEN,QISTEN,QNISTEN,QCSTEN, &
1095 nc1d, nc_tend1d, iinum, C2PREC,CSED,ISED,SSED,GSED,RSED & !wrf-chem
1096 ,no_src_types_cu,maerosol,naer & !TWG add prescribed aerosol
1098 ,has_wetscav, rainprod1d, evapprod1d & !wrf-chem
1102 ! Transfer 1D arrays back into 3D arrays
1106 ! hm, add tendencies to update global variables
1107 ! HM, TENDENCIES FOR Q AND N NOW ADDED IN M2005MICRO, SO WE
1108 ! ONLY NEED TO TRANSFER 1D VARIABLES BACK TO 3D
1121 TH(I,K,J) = T(i,k,j)/PII(i,k,j) ! CONVERT TEMP BACK TO POTENTIAL TEMP
1124 EFFC(i,k,j) = EFFC1D(k)
1125 EFFI(i,k,j) = EFFI1D(k)
1126 EFFS(i,k,j) = EFFS1D(k)
1127 EFFR(i,k,j) = EFFR1D(k)
1128 EFFG(I,K,j) = EFFG1D(K)
1130 !TWG add to send effective radius to RRTMG
1131 IF (aercu_opt.gt.0) THEN
1132 EFCG(i,k,j) = MAX(2.49, MIN(EFFC1D(k), 50.))
1133 EFIG(i,k,j) = MAX(4.99, MIN(EFFI1D(k), 120.))
1134 EFSG(i,k,j) = MAX(9.99, MIN(EFFS1D(k), 999.))
1135 WACT(i,k,j) = WVAR(i,k,j)+W(i,k,j)
1139 IF (flag_qndrop .AND. PRESENT( qndrop )) THEN
1140 qndrop(i,k,j) = nc1d(k)
1141 !jdf CSED3D(I,K,J) = CSED(K)
1142 !TWG 2016 modify for prescribed aerosol
1145 !TWG end modification
1147 IF ( PRESENT( QLSINK ) ) THEN
1148 if(qc(i,k,j)>1.e-10) then
1149 QLSINK(I,K,J) = C2PREC(K)/QC(I,K,J)
1154 IF ( PRESENT( PRECR ) ) PRECR(I,K,J) = RSED(K)
1155 IF ( PRESENT( PRECI ) ) PRECI(I,K,J) = ISED(K)
1156 IF ( PRESENT( PRECS ) ) PRECS(I,K,J) = SSED(K)
1157 IF ( PRESENT( PRECG ) ) PRECG(I,K,J) = GSED(K)
1158 ! EFFECTIVE RADIUS FOR RADIATION CODE (currently not coupled)
1159 ! HM, ADD LIMIT TO PREVENT BLOWING UP OPTICAL PROPERTIES, 8/18/07
1160 ! EFFCS(I,K,J) = MIN(EFFC(I,K,J),50.)
1161 ! EFFCS(I,K,J) = MAX(EFFCS(I,K,J),1.)
1162 ! EFFIS(I,K,J) = MIN(EFFI(I,K,J),130.)
1163 ! EFFIS(I,K,J) = MAX(EFFIS(I,K,J),13.)
1165 #if ( WRF_CHEM == 1)
1166 IF ( has_wetscav ) THEN
1167 IF ( PRESENT( rainprod ) ) rainprod(i,k,j) = rainprod1d(k)
1168 IF ( PRESENT( evapprod ) ) evapprod(i,k,j) = evapprod1d(k)
1174 ! hm modified so that m2005 precip variables correctly match wrf precip variables
1175 RAINNC(i,j) = RAINNC(I,J)+PRECPRT1D
1176 RAINNCV(i,j) = PRECPRT1D
1178 SNOWNC(i,j) = SNOWNC(I,J)+SNOWPRT1D
1179 SNOWNCV(i,j) = SNOWPRT1D
1180 GRAUPELNC(i,j) = GRAUPELNC(I,J)+GRPLPRT1D
1181 GRAUPELNCV(i,j) = GRPLPRT1D
1182 SR(i,j) = SNOWRT1D/(PRECPRT1D+1.E-12)
1184 !+---+-----------------------------------------------------------------+
1185 IF ( PRESENT (diagflag) ) THEN
1186 if (diagflag .and. do_radar_ref == 1) then
1187 call refl10cm_hm (qv1d, qr1d, nr1d, qs1d, ns1d, qg1d, ng1d, &
1188 t1d, p1d, dBZ, kts, kte, i, j)
1190 refl_10cm(i,k,j) = MAX(-35., dBZ(k))
1193 NRRAD = nr1d + (NRMSCU1D*CU_UAF(i,j)) ! TWG add
1194 QRRAD = qr1d + (QRMSCU1D*CU_UAF(i,j)) ! TWG add
1195 NSRAD = ns1d + (NSMSCU1D*CU_UAF(i,j)) ! TWG add
1196 QSRAD = qs1d + (QSMSCU1D*CU_UAF(i,j)) ! TWG add
1197 call refl10cm_hm (qv1d, QRRAD, NRRAD, QSRAD, NSRAD, qg1d, ng1d, & ! Modified TWG for cumulus hydrometeors
1198 t1d, p1d, dBZ2, kts, kte, i, j)
1200 mskf_refl_10cm(i,k,j) = MAX(-35., dBZ2(k))
1204 !+---+-----------------------------------------------------------------+
1209 END SUBROUTINE MP_MORR_TWO_MOMENT_AERO
1211 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1212 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1213 SUBROUTINE MORR_TWO_MOMENT_MICRO(QC3DTEN,QI3DTEN,QNI3DTEN,QR3DTEN, &
1214 NI3DTEN,NS3DTEN,NR3DTEN,QC3D,QI3D,QNI3D,QR3D,NI3D,NS3D,NR3D, &
1215 T3DTEN,QV3DTEN,T3D,QV3D,PRES,DZQ,W3D,WVAR,PRECRT,SNOWRT, &
1216 SNOWPRT,GRPLPRT, & ! hm added 7/13/13
1217 EFFC,EFFI,EFFS,EFFR,DT, &
1218 IMS,IME, JMS,JME, KMS,KME, &
1219 ITS,ITE, JTS,JTE, KTS,KTE, & ! ADD GRAUPEL
1220 QG3DTEN,NG3DTEN,QG3D,NG3D,EFFG,qrcu1d,qscu1d, qicu1d, &
1221 QGSTEN,QRSTEN,QISTEN,QNISTEN,QCSTEN, &
1222 nc3d,nc3dten,iinum, & ! wrf-chem
1223 c2prec,CSED,ISED,SSED,GSED,RSED & ! hm added, wrf-chem
1224 ,no_src_types_cu,maerosol,naer & ! TWG add
1227 ,has_wetscav,rainprod, evapprod &
1231 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1232 ! THIS PROGRAM IS THE MAIN TWO-MOMENT MICROPHYSICS SUBROUTINE DESCRIBED BY
1233 ! MORRISON ET AL. 2005 JAS AND MORRISON ET AL. 2009 MWR
1235 ! THIS SCHEME IS A BULK DOUBLE-MOMENT SCHEME THAT PREDICTS MIXING
1236 ! RATIOS AND NUMBER CONCENTRATIONS OF FIVE HYDROMETEOR SPECIES:
1237 ! CLOUD DROPLETS, CLOUD (SMALL) ICE, RAIN, SNOW, AND GRAUPEL/HAIL.
1239 ! CODE STRUCTURE: MAIN SUBROUTINE IS 'MORR_TWO_MOMENT'. ALSO INCLUDED IN THIS FILE IS
1240 ! 'FUNCTION POLYSVP', 'FUNCTION DERF1', AND
1243 ! NOTE: THIS SUBROUTINE USES 1D ARRAY IN VERTICAL (COLUMN), EVEN THOUGH VARIABLES ARE CALLED '3D'......
1245 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1250 SAVE !TWG 2017 add to avoid memeory issues
1252 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1253 ! THESE VARIABLES BELOW MUST BE LINKED WITH THE MAIN MODEL.
1254 ! DEFINE ARRAY SIZES
1256 ! INPUT NUMBER OF GRID CELLS
1258 ! INPUT/OUTPUT PARAMETERS ! DESCRIPTION (UNITS)
1259 INTEGER, INTENT( IN) :: IMS,IME, JMS,JME, KMS,KME, &
1260 ITS,ITE, JTS,JTE, KTS,KTE
1262 REAL, DIMENSION(KTS:KTE) :: QC3DTEN ! CLOUD WATER MIXING RATIO TENDENCY (KG/KG/S)
1263 REAL, DIMENSION(KTS:KTE) :: QI3DTEN ! CLOUD ICE MIXING RATIO TENDENCY (KG/KG/S)
1264 REAL, DIMENSION(KTS:KTE) :: QNI3DTEN ! SNOW MIXING RATIO TENDENCY (KG/KG/S)
1265 REAL, DIMENSION(KTS:KTE) :: QR3DTEN ! RAIN MIXING RATIO TENDENCY (KG/KG/S)
1266 REAL, DIMENSION(KTS:KTE) :: NI3DTEN ! CLOUD ICE NUMBER CONCENTRATION (1/KG/S)
1267 REAL, DIMENSION(KTS:KTE) :: NS3DTEN ! SNOW NUMBER CONCENTRATION (1/KG/S)
1268 REAL, DIMENSION(KTS:KTE) :: NR3DTEN ! RAIN NUMBER CONCENTRATION (1/KG/S)
1269 REAL, DIMENSION(KTS:KTE) :: QC3D ! CLOUD WATER MIXING RATIO (KG/KG)
1270 REAL, DIMENSION(KTS:KTE) :: QI3D ! CLOUD ICE MIXING RATIO (KG/KG)
1271 REAL, DIMENSION(KTS:KTE) :: QNI3D ! SNOW MIXING RATIO (KG/KG)
1272 REAL, DIMENSION(KTS:KTE) :: QR3D ! RAIN MIXING RATIO (KG/KG)
1273 REAL, DIMENSION(KTS:KTE) :: NI3D ! CLOUD ICE NUMBER CONCENTRATION (1/KG)
1274 REAL, DIMENSION(KTS:KTE) :: NS3D ! SNOW NUMBER CONCENTRATION (1/KG)
1275 REAL, DIMENSION(KTS:KTE) :: NR3D ! RAIN NUMBER CONCENTRATION (1/KG)
1276 REAL, DIMENSION(KTS:KTE) :: T3DTEN ! TEMPERATURE TENDENCY (K/S)
1277 REAL, DIMENSION(KTS:KTE) :: QV3DTEN ! WATER VAPOR MIXING RATIO TENDENCY (KG/KG/S)
1278 REAL, DIMENSION(KTS:KTE) :: T3D ! TEMPERATURE (K)
1279 REAL, DIMENSION(KTS:KTE) :: QV3D ! WATER VAPOR MIXING RATIO (KG/KG)
1280 REAL, DIMENSION(KTS:KTE) :: PRES ! ATMOSPHERIC PRESSURE (PA)
1281 REAL, DIMENSION(KTS:KTE) :: DZQ ! DIFFERENCE IN HEIGHT ACROSS LEVEL (m)
1282 REAL, DIMENSION(KTS:KTE) :: W3D ! GRID-SCALE VERTICAL VELOCITY (M/S)
1283 REAL, DIMENSION(KTS:KTE) :: WVAR ! SUB-GRID VERTICAL VELOCITY (M/S)
1284 ! below for wrf-chem
1285 REAL, DIMENSION(KTS:KTE) :: nc3d
1286 REAL, DIMENSION(KTS:KTE) :: nc3dten
1287 integer, intent(in) :: iinum
1289 ! HM ADDED GRAUPEL VARIABLES
1290 REAL, DIMENSION(KTS:KTE) :: QG3DTEN ! GRAUPEL MIX RATIO TENDENCY (KG/KG/S)
1291 REAL, DIMENSION(KTS:KTE) :: NG3DTEN ! GRAUPEL NUMB CONC TENDENCY (1/KG/S)
1292 REAL, DIMENSION(KTS:KTE) :: QG3D ! GRAUPEL MIX RATIO (KG/KG)
1293 REAL, DIMENSION(KTS:KTE) :: NG3D ! GRAUPEL NUMBER CONC (1/KG)
1295 ! HM, ADD 1/16/07, SEDIMENTATION TENDENCIES FOR MIXING RATIO
1297 REAL, DIMENSION(KTS:KTE) :: QGSTEN ! GRAUPEL SED TEND (KG/KG/S)
1298 REAL, DIMENSION(KTS:KTE) :: QRSTEN ! RAIN SED TEND (KG/KG/S)
1299 REAL, DIMENSION(KTS:KTE) :: QISTEN ! CLOUD ICE SED TEND (KG/KG/S)
1300 REAL, DIMENSION(KTS:KTE) :: QNISTEN ! SNOW SED TEND (KG/KG/S)
1301 REAL, DIMENSION(KTS:KTE) :: QCSTEN ! CLOUD WAT SED TEND (KG/KG/S)
1303 ! hm add cumulus tendencies for precip
1304 REAL, DIMENSION(KTS:KTE) :: qrcu1d
1305 REAL, DIMENSION(KTS:KTE) :: qscu1d
1306 REAL, DIMENSION(KTS:KTE) :: qicu1d
1310 REAL PRECRT ! TOTAL PRECIP PER TIME STEP (mm)
1311 REAL SNOWRT ! SNOW PER TIME STEP (mm)
1313 REAL SNOWPRT ! TOTAL CLOUD ICE PLUS SNOW PER TIME STEP (mm)
1314 REAL GRPLPRT ! TOTAL GRAUPEL PER TIME STEP (mm)
1316 REAL, DIMENSION(KTS:KTE) :: EFFC ! DROPLET EFFECTIVE RADIUS (MICRON)
1317 REAL, DIMENSION(KTS:KTE) :: EFFI ! CLOUD ICE EFFECTIVE RADIUS (MICRON)
1318 REAL, DIMENSION(KTS:KTE) :: EFFS ! SNOW EFFECTIVE RADIUS (MICRON)
1319 REAL, DIMENSION(KTS:KTE) :: EFFR ! RAIN EFFECTIVE RADIUS (MICRON)
1320 REAL, DIMENSION(KTS:KTE) :: EFFG ! GRAUPEL EFFECTIVE RADIUS (MICRON)
1322 ! MODEL INPUT PARAMETERS (FORMERLY IN COMMON BLOCKS)
1324 REAL DT ! MODEL TIME STEP (SEC)
1326 !.....................................................................................................
1327 ! LOCAL VARIABLES: ALL PARAMETERS BELOW ARE LOCAL TO SCHEME AND DON'T NEED TO COMMUNICATE WITH THE
1328 ! REST OF THE MODEL.
1330 ! SIZE PARAMETER VARIABLES
1332 REAL, DIMENSION(KTS:KTE) :: LAMC ! SLOPE PARAMETER FOR DROPLETS (M-1)
1333 REAL, DIMENSION(KTS:KTE) :: LAMI ! SLOPE PARAMETER FOR CLOUD ICE (M-1)
1334 REAL, DIMENSION(KTS:KTE) :: LAMS ! SLOPE PARAMETER FOR SNOW (M-1)
1335 REAL, DIMENSION(KTS:KTE) :: LAMR ! SLOPE PARAMETER FOR RAIN (M-1)
1336 REAL, DIMENSION(KTS:KTE) :: LAMG ! SLOPE PARAMETER FOR GRAUPEL (M-1)
1337 REAL, DIMENSION(KTS:KTE) :: CDIST1 ! PSD PARAMETER FOR DROPLETS
1338 REAL, DIMENSION(KTS:KTE) :: N0I ! INTERCEPT PARAMETER FOR CLOUD ICE (KG-1 M-1)
1339 REAL, DIMENSION(KTS:KTE) :: N0S ! INTERCEPT PARAMETER FOR SNOW (KG-1 M-1)
1340 REAL, DIMENSION(KTS:KTE) :: N0RR ! INTERCEPT PARAMETER FOR RAIN (KG-1 M-1)
1341 REAL, DIMENSION(KTS:KTE) :: N0G ! INTERCEPT PARAMETER FOR GRAUPEL (KG-1 M-1)
1342 REAL, DIMENSION(KTS:KTE) :: PGAM ! SPECTRAL SHAPE PARAMETER FOR DROPLETS
1344 ! MICROPHYSICAL PROCESSES
1346 REAL, DIMENSION(KTS:KTE) :: NSUBC ! LOSS OF NC DURING EVAP
1347 REAL, DIMENSION(KTS:KTE) :: NSUBI ! LOSS OF NI DURING SUB.
1348 REAL, DIMENSION(KTS:KTE) :: NSUBS ! LOSS OF NS DURING SUB.
1349 REAL, DIMENSION(KTS:KTE) :: NSUBR ! LOSS OF NR DURING EVAP
1350 REAL, DIMENSION(KTS:KTE) :: PRD ! DEP CLOUD ICE
1351 REAL, DIMENSION(KTS:KTE) :: PRE ! EVAP OF RAIN
1352 REAL, DIMENSION(KTS:KTE) :: PRDS ! DEP SNOW
1353 REAL, DIMENSION(KTS:KTE) :: NNUCCC ! CHANGE N DUE TO CONTACT FREEZ DROPLETS
1354 REAL, DIMENSION(KTS:KTE) :: MNUCCC ! CHANGE Q DUE TO CONTACT FREEZ DROPLETS
1355 REAL, DIMENSION(KTS:KTE) :: PRA ! ACCRETION DROPLETS BY RAIN
1356 REAL, DIMENSION(KTS:KTE) :: PRC ! AUTOCONVERSION DROPLETS
1357 REAL, DIMENSION(KTS:KTE) :: PCC ! COND/EVAP DROPLETS
1358 REAL, DIMENSION(KTS:KTE) :: NNUCCD ! CHANGE N FREEZING AEROSOL (PRIM ICE NUCLEATION)
1359 REAL, DIMENSION(KTS:KTE) :: MNUCCD ! CHANGE Q FREEZING AEROSOL (PRIM ICE NUCLEATION)
1360 REAL, DIMENSION(KTS:KTE) :: MNUCCR ! CHANGE Q DUE TO CONTACT FREEZ RAIN
1361 REAL, DIMENSION(KTS:KTE) :: NNUCCR ! CHANGE N DUE TO CONTACT FREEZ RAIN
1362 REAL, DIMENSION(KTS:KTE) :: NPRA ! CHANGE IN N DUE TO DROPLET ACC BY RAIN
1363 REAL, DIMENSION(KTS:KTE) :: NRAGG ! SELF-COLLECTION/BREAKUP OF RAIN
1364 REAL, DIMENSION(KTS:KTE) :: NSAGG ! SELF-COLLECTION OF SNOW
1365 REAL, DIMENSION(KTS:KTE) :: NPRC ! CHANGE NC AUTOCONVERSION DROPLETS
1366 REAL, DIMENSION(KTS:KTE) :: NPRC1 ! CHANGE NR AUTOCONVERSION DROPLETS
1367 REAL, DIMENSION(KTS:KTE) :: PRAI ! CHANGE Q ACCRETION CLOUD ICE BY SNOW
1368 REAL, DIMENSION(KTS:KTE) :: PRCI ! CHANGE Q AUTOCONVERSIN CLOUD ICE TO SNOW
1369 REAL, DIMENSION(KTS:KTE) :: PSACWS ! CHANGE Q DROPLET ACCRETION BY SNOW
1370 REAL, DIMENSION(KTS:KTE) :: NPSACWS ! CHANGE N DROPLET ACCRETION BY SNOW
1371 REAL, DIMENSION(KTS:KTE) :: PSACWI ! CHANGE Q DROPLET ACCRETION BY CLOUD ICE
1372 REAL, DIMENSION(KTS:KTE) :: NPSACWI ! CHANGE N DROPLET ACCRETION BY CLOUD ICE
1373 REAL, DIMENSION(KTS:KTE) :: NPRCI ! CHANGE N AUTOCONVERSION CLOUD ICE BY SNOW
1374 REAL, DIMENSION(KTS:KTE) :: NPRAI ! CHANGE N ACCRETION CLOUD ICE
1375 REAL, DIMENSION(KTS:KTE) :: NMULTS ! ICE MULT DUE TO RIMING DROPLETS BY SNOW
1376 REAL, DIMENSION(KTS:KTE) :: NMULTR ! ICE MULT DUE TO RIMING RAIN BY SNOW
1377 REAL, DIMENSION(KTS:KTE) :: QMULTS ! CHANGE Q DUE TO ICE MULT DROPLETS/SNOW
1378 REAL, DIMENSION(KTS:KTE) :: QMULTR ! CHANGE Q DUE TO ICE RAIN/SNOW
1379 REAL, DIMENSION(KTS:KTE) :: PRACS ! CHANGE Q RAIN-SNOW COLLECTION
1380 REAL, DIMENSION(KTS:KTE) :: NPRACS ! CHANGE N RAIN-SNOW COLLECTION
1381 REAL, DIMENSION(KTS:KTE) :: PCCN ! CHANGE Q DROPLET ACTIVATION
1382 REAL, DIMENSION(KTS:KTE) :: PSMLT ! CHANGE Q MELTING SNOW TO RAIN
1383 REAL, DIMENSION(KTS:KTE) :: EVPMS ! CHNAGE Q MELTING SNOW EVAPORATING
1384 REAL, DIMENSION(KTS:KTE) :: NSMLTS ! CHANGE N MELTING SNOW
1385 REAL, DIMENSION(KTS:KTE) :: NSMLTR ! CHANGE N MELTING SNOW TO RAIN
1387 REAL, DIMENSION(KTS:KTE) :: PIACR ! CHANGE QR, ICE-RAIN COLLECTION
1388 REAL, DIMENSION(KTS:KTE) :: NIACR ! CHANGE N, ICE-RAIN COLLECTION
1389 REAL, DIMENSION(KTS:KTE) :: PRACI ! CHANGE QI, ICE-RAIN COLLECTION
1390 REAL, DIMENSION(KTS:KTE) :: PIACRS ! CHANGE QR, ICE RAIN COLLISION, ADDED TO SNOW
1391 REAL, DIMENSION(KTS:KTE) :: NIACRS ! CHANGE N, ICE RAIN COLLISION, ADDED TO SNOW
1392 REAL, DIMENSION(KTS:KTE) :: PRACIS ! CHANGE QI, ICE RAIN COLLISION, ADDED TO SNOW
1393 REAL, DIMENSION(KTS:KTE) :: EPRD ! SUBLIMATION CLOUD ICE
1394 REAL, DIMENSION(KTS:KTE) :: EPRDS ! SUBLIMATION SNOW
1395 ! HM ADDED GRAUPEL PROCESSES
1396 REAL, DIMENSION(KTS:KTE) :: PRACG ! CHANGE IN Q COLLECTION RAIN BY GRAUPEL
1397 REAL, DIMENSION(KTS:KTE) :: PSACWG ! CHANGE IN Q COLLECTION DROPLETS BY GRAUPEL
1398 REAL, DIMENSION(KTS:KTE) :: PGSACW ! CONVERSION Q TO GRAUPEL DUE TO COLLECTION DROPLETS BY SNOW
1399 REAL, DIMENSION(KTS:KTE) :: PGRACS ! CONVERSION Q TO GRAUPEL DUE TO COLLECTION RAIN BY SNOW
1400 REAL, DIMENSION(KTS:KTE) :: PRDG ! DEP OF GRAUPEL
1401 REAL, DIMENSION(KTS:KTE) :: EPRDG ! SUB OF GRAUPEL
1402 REAL, DIMENSION(KTS:KTE) :: EVPMG ! CHANGE Q MELTING OF GRAUPEL AND EVAPORATION
1403 REAL, DIMENSION(KTS:KTE) :: PGMLT ! CHANGE Q MELTING OF GRAUPEL
1404 REAL, DIMENSION(KTS:KTE) :: NPRACG ! CHANGE N COLLECTION RAIN BY GRAUPEL
1405 REAL, DIMENSION(KTS:KTE) :: NPSACWG ! CHANGE N COLLECTION DROPLETS BY GRAUPEL
1406 REAL, DIMENSION(KTS:KTE) :: NSCNG ! CHANGE N CONVERSION TO GRAUPEL DUE TO COLLECTION DROPLETS BY SNOW
1407 REAL, DIMENSION(KTS:KTE) :: NGRACS ! CHANGE N CONVERSION TO GRAUPEL DUE TO COLLECTION RAIN BY SNOW
1408 REAL, DIMENSION(KTS:KTE) :: NGMLTG ! CHANGE N MELTING GRAUPEL
1409 REAL, DIMENSION(KTS:KTE) :: NGMLTR ! CHANGE N MELTING GRAUPEL TO RAIN
1410 REAL, DIMENSION(KTS:KTE) :: NSUBG ! CHANGE N SUB/DEP OF GRAUPEL
1411 REAL, DIMENSION(KTS:KTE) :: PSACR ! CONVERSION DUE TO COLL OF SNOW BY RAIN
1412 REAL, DIMENSION(KTS:KTE) :: NMULTG ! ICE MULT DUE TO ACC DROPLETS BY GRAUPEL
1413 REAL, DIMENSION(KTS:KTE) :: NMULTRG ! ICE MULT DUE TO ACC RAIN BY GRAUPEL
1414 REAL, DIMENSION(KTS:KTE) :: QMULTG ! CHANGE Q DUE TO ICE MULT DROPLETS/GRAUPEL
1415 REAL, DIMENSION(KTS:KTE) :: QMULTRG ! CHANGE Q DUE TO ICE MULT RAIN/GRAUPEL
1417 ! TIME-VARYING ATMOSPHERIC PARAMETERS
1419 REAL, DIMENSION(KTS:KTE) :: KAP ! THERMAL CONDUCTIVITY OF AIR
1420 REAL, DIMENSION(KTS:KTE) :: EVS ! SATURATION VAPOR PRESSURE
1421 REAL, DIMENSION(KTS:KTE) :: EIS ! ICE SATURATION VAPOR PRESSURE
1422 REAL, DIMENSION(KTS:KTE) :: QVS ! SATURATION MIXING RATIO
1423 REAL, DIMENSION(KTS:KTE) :: QVI ! ICE SATURATION MIXING RATIO
1424 REAL, DIMENSION(KTS:KTE) :: QVQVS ! SAUTRATION RATIO
1425 REAL, DIMENSION(KTS:KTE) :: QVQVSI! ICE SATURAION RATIO
1426 REAL, DIMENSION(KTS:KTE) :: DV ! DIFFUSIVITY OF WATER VAPOR IN AIR
1427 REAL, DIMENSION(KTS:KTE) :: XXLS ! LATENT HEAT OF SUBLIMATION
1428 REAL, DIMENSION(KTS:KTE) :: XXLV ! LATENT HEAT OF VAPORIZATION
1429 REAL, DIMENSION(KTS:KTE) :: CPM ! SPECIFIC HEAT AT CONST PRESSURE FOR MOIST AIR
1430 REAL, DIMENSION(KTS:KTE) :: MU ! VISCOCITY OF AIR
1431 REAL, DIMENSION(KTS:KTE) :: SC ! SCHMIDT NUMBER
1432 REAL, DIMENSION(KTS:KTE) :: XLF ! LATENT HEAT OF FREEZING
1433 REAL, DIMENSION(KTS:KTE) :: RHO ! AIR DENSITY
1434 REAL, DIMENSION(KTS:KTE) :: AB ! CORRECTION TO CONDENSATION RATE DUE TO LATENT HEATING
1435 REAL, DIMENSION(KTS:KTE) :: ABI ! CORRECTION TO DEPOSITION RATE DUE TO LATENT HEATING
1437 ! TIME-VARYING MICROPHYSICS PARAMETERS
1439 REAL, DIMENSION(KTS:KTE) :: DAP ! DIFFUSIVITY OF AEROSOL
1440 REAL NACNT ! NUMBER OF CONTACT IN
1441 REAL FMULT ! TEMP.-DEP. PARAMETER FOR RIME-SPLINTERING
1442 REAL COFFI ! ICE AUTOCONVERSION PARAMETER
1444 ! FALL SPEED WORKING VARIABLES (DEFINED IN CODE)
1446 REAL, DIMENSION(KTS:KTE) :: DUMI,DUMR,DUMFNI,DUMG,DUMFNG
1448 REAL, DIMENSION(KTS:KTE) :: FR, FI, FNI,FG,FNG
1450 REAL, DIMENSION(KTS:KTE) :: FALOUTR,FALOUTI,FALOUTNI
1451 REAL FALTNDR,FALTNDI,FALTNDNI,RHO2
1452 REAL, DIMENSION(KTS:KTE) :: DUMQS,DUMFNS
1454 REAL, DIMENSION(KTS:KTE) :: FS,FNS, FALOUTS,FALOUTNS,FALOUTG,FALOUTNG
1455 REAL FALTNDS,FALTNDNS,UNR,FALTNDG,FALTNDNG
1456 REAL, DIMENSION(KTS:KTE) :: DUMC,DUMFNC
1457 REAL UNC,UMC,UNG,UMG
1458 REAL, DIMENSION(KTS:KTE) :: FC,FALOUTC,FALOUTNC
1459 REAL FALTNDC,FALTNDNC
1460 REAL, DIMENSION(KTS:KTE) :: FNC,DUMFNR,FALOUTNR
1462 REAL, DIMENSION(KTS:KTE) :: FNR
1464 ! FALL-SPEED PARAMETER 'A' WITH AIR DENSITY CORRECTION
1466 REAL, DIMENSION(KTS:KTE) :: AIN,ARN,ASN,ACN,AGN
1468 ! EXTERNAL FUNCTION CALL RETURN VARIABLES
1470 ! REAL GAMMA, ! EULER GAMMA FUNCTION
1471 ! REAL POLYSVP, ! SAT. PRESSURE FUNCTION
1472 ! REAL DERF1 ! ERROR FUNCTION
1476 REAL DUM,DUM1,DUM2,DUMT,DUMQV,DUMQSS,DUMQSI,DUMS
1478 ! PROGNOSTIC SUPERSATURATION
1480 REAL DQSDT ! CHANGE OF SAT. MIX. RAT. WITH TEMPERATURE
1481 REAL DQSIDT ! CHANGE IN ICE SAT. MIXING RAT. WITH T
1482 REAL EPSI ! 1/PHASE REL. TIME (SEE M2005), ICE
1483 REAL EPSS ! 1/PHASE REL. TIME (SEE M2005), SNOW
1484 REAL EPSR ! 1/PHASE REL. TIME (SEE M2005), RAIN
1485 REAL EPSG ! 1/PHASE REL. TIME (SEE M2005), GRAUPEL
1487 ! NEW DROPLET ACTIVATION VARIABLES
1488 REAL TAUC ! PHASE REL. TIME (SEE M2005), DROPLETS
1489 REAL TAUR ! PHASE REL. TIME (SEE M2005), RAIN
1490 REAL TAUI ! PHASE REL. TIME (SEE M2005), CLOUD ICE
1491 REAL TAUS ! PHASE REL. TIME (SEE M2005), SNOW
1492 REAL TAUG ! PHASE REL. TIME (SEE M2005), GRAUPEL
1495 ! COUNTING/INDEX VARIABLES
1497 INTEGER K,NSTEP,N,L !TWG add ! ,I
1499 ! LTRUE IS ONLY USED TO SPEED UP THE CODE !!
1500 ! LTRUE, SWITCH = 0, NO HYDROMETEORS IN COLUMN,
1501 ! = 1, HYDROMETEORS IN COLUMN
1505 ! DROPLET ACTIVATION/FREEZING AEROSOL
1508 REAL CT ! DROPLET ACTIVATION PARAMETER
1509 REAL TEMP1 ! DUMMY TEMPERATURE
1510 REAL SAT1 ! DUMMY SATURATION
1511 REAL SIGVL ! SURFACE TENSION LIQ/VAPOR
1512 REAL KEL ! KELVIN PARAMETER
1513 REAL KC2 ! TOTAL ICE NUCLEATION RATE
1515 REAL KC2H ! ICE NUCLEATED FROM HOMOGENOUS FREEZING OF SULFATE
1516 REAL KC2IM ! ICE NUCLEATED FROM IMMERSION FREEZING
1517 REAL KC2D ! ICE NUCLEATED FROM DEPOSITION NUCLEATION
1518 REAL KC2DM ! ICE NUCLEATED FROM MIXED PHASE DEPOSITION NUCLEATION
1520 REAL CRY,KRY ! AEROSOL ACTIVATION PARAMETERS
1522 ! MORE WORKING/DUMMY VARIABLES
1524 REAL DUMQI,DUMNI,DC0,DS0,DG0
1525 REAL DUMQC,DUMQR,RATIO,SUM_DEP,FUDGEF
1527 ! EFFECTIVE VERTICAL VELOCITY (M/S)
1530 ! WORKING PARAMETERS FOR ICE NUCLEATION
1534 ! WORKING PARAMETERS FOR AEROSOL ACTIVATION
1536 REAL AACT,GAMM,GG,PSI,ETA1,ETA2,SM1,SM2,SMAX,UU1,UU2,ALPHA
1538 ! DUMMY SIZE DISTRIBUTION PARAMETERS
1540 REAL DLAMS,DLAMR,DLAMI,DLAMC,DLAMG,LAMMAX,LAMMIN
1544 ! TWG 2016 Begin FOR Prescribed Aerosol
1545 INTEGER no_src_types_cu
1546 REAL, DIMENSION(KTS:KTE,no_src_types_cu) :: maerosol,naer
1547 REAL, DIMENSION(no_src_types_cu) :: maero,naero
1552 REAL, DIMENSION(KTS:KTE)::C2PREC,CSED,ISED,SSED,GSED,RSED
1554 REAL, DIMENSION(KTS:KTE), INTENT(INOUT) :: rainprod, evapprod
1555 LOGICAL, INTENT(IN) :: has_wetscav
1557 REAL, DIMENSION(KTS:KTE) :: tqimelt ! melting of cloud ice (tendency)
1559 ! comment lines for wrf-chem since these are intent(in) in that case
1560 ! REAL, DIMENSION(KTS:KTE) :: NC3DTEN ! CLOUD DROPLET NUMBER CONCENTRATION (1/KG/S)
1561 ! REAL, DIMENSION(KTS:KTE) :: NC3D ! CLOUD DROPLET NUMBER CONCENTRATION (1/KG)
1563 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1565 ! SET LTRUE INITIALLY TO 0
1569 ! ATMOSPHERIC PARAMETERS THAT VARY IN TIME AND HEIGHT
1572 ! NC3DTEN LOCAL ARRAY INITIALIZED
1574 ! TWG/amy comment out this initialization
1575 ! INITIALIZE VARIABLES FOR WRF-CHEM OUTPUT TO ZERO
1593 !Make single point aerosol number and mass
1594 DO L=1,no_src_types_cu
1595 maero(L)=maerosol(K,L)
1601 ! LATENT HEAT OF VAPORATION
1603 XXLV(K) = 3.1484E6-2370.*T3D(K)
1605 ! LATENT HEAT OF SUBLIMATION
1607 XXLS(K) = 3.15E6-2370.*T3D(K)+0.3337E6
1609 CPM(K) = CP*(1.+0.887*QV3D(K))
1611 ! SATURATION VAPOR PRESSURE AND MIXING RATIO
1613 ! hm, add fix for low pressure, 5/12/10
1614 EVS(K) = min(0.99*pres(k),POLYSVP(T3D(K),0)) ! PA
1615 EIS(K) = min(0.99*pres(k),POLYSVP(T3D(K),1)) ! PA
1617 ! MAKE SURE ICE SATURATION DOESN'T EXCEED WATER SAT. NEAR FREEZING
1619 IF (EIS(K).GT.EVS(K)) EIS(K) = EVS(K)
1621 QVS(K) = EP_2*EVS(K)/(PRES(K)-EVS(K))
1622 QVI(K) = EP_2*EIS(K)/(PRES(K)-EIS(K))
1624 QVQVS(K) = QV3D(K)/QVS(K)
1625 QVQVSI(K) = QV3D(K)/QVI(K)
1629 RHO(K) = PRES(K)/(R*T3D(K))
1631 ! ADD NUMBER CONCENTRATION DUE TO CUMULUS TENDENCY
1632 ! ASSUME N0 ASSOCIATED WITH CUMULUS PARAM RAIN IS 10^7 M^-4
1633 ! ASSUME N0 ASSOCIATED WITH CUMULUS PARAM SNOW IS 2 X 10^7 M^-4
1634 ! FOR DETRAINED CLOUD ICE, ASSUME MEAN VOLUME DIAM OF 80 MICRON
1636 IF (QRCU1D(K).GE.1.E-10) THEN
1637 DUM=1.8e5*(QRCU1D(K)*DT/(PI*RHOW*RHO(K)**3))**0.25
1640 IF (QSCU1D(K).GE.1.E-10) THEN
1641 DUM=3.e5*(QSCU1D(K)*DT/(CONS1*RHO(K)**3))**(1./(DS+1.))
1644 IF (QICU1D(K).GE.1.E-10) THEN
1645 DUM=QICU1D(K)*DT/(CI*(80.E-6)**DI)
1649 ! AT SUBSATURATION, REMOVE SMALL AMOUNTS OF CLOUD/PRECIP WATER
1650 ! hm modify 7/0/09 change limit to 1.e-8
1652 IF (QVQVS(K).LT.0.9) THEN
1653 IF (QR3D(K).LT.1.E-8) THEN
1654 QV3D(K)=QV3D(K)+QR3D(K)
1655 T3D(K)=T3D(K)-QR3D(K)*XXLV(K)/CPM(K)
1658 IF (QC3D(K).LT.1.E-8) THEN
1659 QV3D(K)=QV3D(K)+QC3D(K)
1660 T3D(K)=T3D(K)-QC3D(K)*XXLV(K)/CPM(K)
1665 IF (QVQVSI(K).LT.0.9) THEN
1666 IF (QI3D(K).LT.1.E-8) THEN
1667 QV3D(K)=QV3D(K)+QI3D(K)
1668 T3D(K)=T3D(K)-QI3D(K)*XXLS(K)/CPM(K)
1671 IF (QNI3D(K).LT.1.E-8) THEN
1672 QV3D(K)=QV3D(K)+QNI3D(K)
1673 T3D(K)=T3D(K)-QNI3D(K)*XXLS(K)/CPM(K)
1676 IF (QG3D(K).LT.1.E-8) THEN
1677 QV3D(K)=QV3D(K)+QG3D(K)
1678 T3D(K)=T3D(K)-QG3D(K)*XXLS(K)/CPM(K)
1685 XLF(K) = XXLS(K)-XXLV(K)
1687 !..................................................................
1688 ! IF MIXING RATIO < QSMALL SET MIXING RATIO AND NUMBER CONC TO ZERO
1690 IF (QC3D(K).LT.QSMALL) THEN
1695 IF (QR3D(K).LT.QSMALL) THEN
1700 IF (QI3D(K).LT.QSMALL) THEN
1705 IF (QNI3D(K).LT.QSMALL) THEN
1710 IF (QG3D(K).LT.QSMALL) THEN
1716 ! INITIALIZE SEDIMENTATION TENDENCIES FOR MIXING RATIO
1724 !..................................................................
1725 ! MICROPHYSICS PARAMETERS VARYING IN TIME/HEIGHT
1728 MU(K) = 1.496E-6*T3D(K)**1.5/(T3D(K)+120.)
1730 ! FALL SPEED WITH DENSITY CORRECTION (HEYMSFIELD AND BENSSEMER 2006)
1732 DUM = (RHOSU/RHO(K))**0.54
1736 ! AA revision 4/1/11: Ikawa and Saito 1991 air-density correction
1737 AIN(K) = (RHOSU/RHO(K))**0.35*AI
1741 ! AA revision 4/1/11: temperature-dependent Stokes fall speed
1742 ACN(K) = G*RHOW/(18.*MU(K))
1743 ! HM ADD GRAUPEL 8/28/06
1746 !hm 4/7/09 bug fix, initialize lami to prevent later division by zero
1749 !..................................
1750 ! IF THERE IS NO CLOUD/PRECIP WATER, AND IF SUBSATURATED, THEN SKIP MICROPHYSICS
1753 IF (QC3D(K).LT.QSMALL.AND.QI3D(K).LT.QSMALL.AND.QNI3D(K).LT.QSMALL &
1754 .AND.QR3D(K).LT.QSMALL.AND.QG3D(K).LT.QSMALL) THEN
1755 IF (T3D(K).LT.273.15.AND.QVQVSI(K).LT.0.999) GOTO 200
1756 IF (T3D(K).GE.273.15.AND.QVQVS(K).LT.0.999) GOTO 200
1759 ! THERMAL CONDUCTIVITY FOR AIR
1762 KAP(K) = 1.414E3*MU(K)
1764 ! DIFFUSIVITY OF WATER VAPOR
1766 DV(K) = 8.794E-5*T3D(K)**1.81/PRES(K)
1771 SC(K) = MU(K)/(RHO(K)*DV(K))
1773 ! PSYCHOMETIC CORRECTIONS
1775 ! RATE OF CHANGE SAT. MIX. RATIO WITH TEMPERATURE
1777 DUM = (RV*T3D(K)**2)
1779 DQSDT = XXLV(K)*QVS(K)/DUM
1780 DQSIDT = XXLS(K)*QVI(K)/DUM
1782 ABI(K) = 1.+DQSIDT*XXLS(K)/CPM(K)
1783 AB(K) = 1.+DQSDT*XXLV(K)/CPM(K)
1786 !.....................................................................
1787 !.....................................................................
1788 ! CASE FOR TEMPERATURE ABOVE FREEZING
1790 IF (T3D(K).GE.273.15) THEN
1792 !......................................................................
1793 !HM ADD, ALLOW FOR CONSTANT DROPLET NUMBER
1794 ! INUM = 0, PREDICT DROPLET NUMBER
1795 ! INUM = 1, SET CONSTANT DROPLET NUMBER
1797 IF (iinum.EQ.1) THEN
1798 ! CONVERT NDCNST FROM CM-3 TO KG-1
1799 NC3D(K)=NDCNST*1.E6/RHO(K)
1802 ! GET SIZE DISTRIBUTION PARAMETERS
1804 ! MELT VERY SMALL SNOW AND GRAUPEL MIXING RATIOS, ADD TO RAIN
1805 IF (QNI3D(K).LT.1.E-6) THEN
1806 QR3D(K)=QR3D(K)+QNI3D(K)
1807 NR3D(K)=NR3D(K)+NS3D(K)
1808 T3D(K)=T3D(K)-QNI3D(K)*XLF(K)/CPM(K)
1812 IF (QG3D(K).LT.1.E-6) THEN
1813 QR3D(K)=QR3D(K)+QG3D(K)
1814 NR3D(K)=NR3D(K)+NG3D(K)
1815 T3D(K)=T3D(K)-QG3D(K)*XLF(K)/CPM(K)
1820 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
1822 ! MAKE SURE NUMBER CONCENTRATIONS AREN'T NEGATIVE
1824 NS3D(K) = MAX(0.,NS3D(K))
1825 NC3D(K) = MAX(0.,NC3D(K))
1826 NR3D(K) = MAX(0.,NR3D(K))
1827 NG3D(K) = MAX(0.,NG3D(K))
1829 !......................................................................
1832 IF (QR3D(K).GE.QSMALL) THEN
1833 LAMR(K) = (PI*RHOW*NR3D(K)/QR3D(K))**(1./3.)
1834 N0RR(K) = NR3D(K)*LAMR(K)
1840 IF (LAMR(K).LT.LAMMINR) THEN
1844 N0RR(K) = LAMR(K)**4*QR3D(K)/(PI*RHOW)
1846 NR3D(K) = N0RR(K)/LAMR(K)
1847 ELSE IF (LAMR(K).GT.LAMMAXR) THEN
1849 N0RR(K) = LAMR(K)**4*QR3D(K)/(PI*RHOW)
1851 NR3D(K) = N0RR(K)/LAMR(K)
1855 !......................................................................
1858 ! MARTIN ET AL. (1994) FORMULA FOR PGAM
1860 IF (QC3D(K).GE.QSMALL) THEN
1862 DUM = PRES(K)/(287.15*T3D(K))
1863 PGAM(K)=0.0005714*(NC3D(K)/1.E6*DUM)+0.2714
1864 PGAM(K)=1./(PGAM(K)**2)-1.
1865 PGAM(K)=MAX(PGAM(K),2.)
1866 PGAM(K)=MIN(PGAM(K),10.)
1870 LAMC(K) = (CONS26*NC3D(K)*GAMMA(PGAM(K)+4.)/ &
1871 (QC3D(K)*GAMMA(PGAM(K)+1.)))**(1./3.)
1873 ! LAMMIN, 60 MICRON DIAMETER
1876 LAMMIN = (PGAM(K)+1.)/60.E-6
1877 LAMMAX = (PGAM(K)+1.)/1.E-6
1879 IF (LAMC(K).LT.LAMMIN) THEN
1882 NC3D(K) = EXP(3.*LOG(LAMC(K))+LOG(QC3D(K))+ &
1883 LOG(GAMMA(PGAM(K)+1.))-LOG(GAMMA(PGAM(K)+4.)))/CONS26
1884 ELSE IF (LAMC(K).GT.LAMMAX) THEN
1887 NC3D(K) = EXP(3.*LOG(LAMC(K))+LOG(QC3D(K))+ &
1888 LOG(GAMMA(PGAM(K)+1.))-LOG(GAMMA(PGAM(K)+4.)))/CONS26
1894 !......................................................................
1897 IF (QNI3D(K).GE.QSMALL) THEN
1898 LAMS(K) = (CONS1*NS3D(K)/QNI3D(K))**(1./DS)
1899 N0S(K) = NS3D(K)*LAMS(K)
1905 IF (LAMS(K).LT.LAMMINS) THEN
1907 N0S(K) = LAMS(K)**4*QNI3D(K)/CONS1
1909 NS3D(K) = N0S(K)/LAMS(K)
1911 ELSE IF (LAMS(K).GT.LAMMAXS) THEN
1914 N0S(K) = LAMS(K)**4*QNI3D(K)/CONS1
1916 NS3D(K) = N0S(K)/LAMS(K)
1920 !......................................................................
1923 IF (QG3D(K).GE.QSMALL) THEN
1924 LAMG(K) = (CONS2*NG3D(K)/QG3D(K))**(1./DG)
1925 N0G(K) = NG3D(K)*LAMG(K)
1929 IF (LAMG(K).LT.LAMMING) THEN
1931 N0G(K) = LAMG(K)**4*QG3D(K)/CONS2
1933 NG3D(K) = N0G(K)/LAMG(K)
1935 ELSE IF (LAMG(K).GT.LAMMAXG) THEN
1938 N0G(K) = LAMG(K)**4*QG3D(K)/CONS2
1940 NG3D(K) = N0G(K)/LAMG(K)
1944 !.....................................................................
1945 ! ZERO OUT PROCESS RATES
1970 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1971 ! CALCULATION OF MICROPHYSICAL PROCESS RATES, T > 273.15 K
1973 !.................................................................
1974 !.......................................................................
1975 ! AUTOCONVERSION OF CLOUD LIQUID WATER TO RAIN
1976 ! FORMULA FROM BEHENG (1994)
1977 ! USING NUMERICAL SIMULATION OF STOCHASTIC COLLECTION EQUATION
1978 ! AND INITIAL CLOUD DROPLET SIZE DISTRIBUTION SPECIFIED
1979 ! AS A GAMMA DISTRIBUTION
1981 ! USE MINIMUM VALUE OF 1.E-6 TO PREVENT FLOATING POINT ERROR
1983 IF (QC3D(K).GE.1.E-6) THEN
1985 ! HM ADD 12/13/06, REPLACE WITH NEWER FORMULA
1986 ! FROM KHAIROUTDINOV AND KOGAN 2000, MWR
1988 PRC(K)=1350.*QC3D(K)**2.47* &
1989 (NC3D(K)/1.e6*RHO(K))**(-1.79)
1991 ! note: nprc1 is change in Nr,
1992 ! nprc is change in Nc
1994 NPRC1(K) = PRC(K)/CONS29
1995 NPRC(K) = PRC(K)/(QC3D(k)/NC3D(K))
1997 ! hm bug fix 3/20/12
1998 NPRC(K) = MIN(NPRC(K),NC3D(K)/DT)
1999 NPRC1(K) = MIN(NPRC1(K),NPRC(K))
2003 !.......................................................................
2004 ! HM ADD 12/13/06, COLLECTION OF SNOW BY RAIN ABOVE FREEZING
2005 ! FORMULA FROM IKAWA AND SAITO (1991)
2007 IF (QR3D(K).GE.1.E-8.AND.QNI3D(K).GE.1.E-8) THEN
2009 UMS = ASN(K)*CONS3/(LAMS(K)**BS)
2010 UMR = ARN(K)*CONS4/(LAMR(K)**BR)
2011 UNS = ASN(K)*CONS5/LAMS(K)**BS
2012 UNR = ARN(K)*CONS6/LAMR(K)**BR
2014 ! SET REASLISTIC LIMITS ON FALLSPEEDS
2017 dum=(rhosu/rho(k))**0.54
2018 UMS=MIN(UMS,1.2*dum)
2019 UNS=MIN(UNS,1.2*dum)
2020 UMR=MIN(UMR,9.1*dum)
2021 UNR=MIN(UNR,9.1*dum)
2024 ! for above freezing conditions to get accelerated melting of snow,
2025 ! we need collection of rain by snow (following Lin et al. 1983)
2026 ! PRACS(K) = CONS31*(((1.2*UMR-0.95*UMS)**2+ &
2027 ! 0.08*UMS*UMR)**0.5*RHO(K)* &
2028 ! N0RR(K)*N0S(K)/LAMS(K)**3* &
2029 ! (5./(LAMS(K)**3*LAMR(K))+ &
2030 ! 2./(LAMS(K)**2*LAMR(K)**2)+ &
2031 ! 0.5/(LAMS(K)*LAMR(K)**3)))
2033 PRACS(K) = CONS41*(((1.2*UMR-0.95*UMS)**2+ &
2034 0.08*UMS*UMR)**0.5*RHO(K)* &
2035 N0RR(K)*N0S(K)/LAMR(K)**3* &
2036 (5./(LAMR(K)**3*LAMS(K))+ &
2037 2./(LAMR(K)**2*LAMS(K)**2)+ &
2038 0.5/(LAMR(k)*LAMS(k)**3)))
2040 ! fix 053011, npracs no longer subtracted from snow
2041 ! NPRACS(K) = CONS32*RHO(K)*(1.7*(UNR-UNS)**2+ &
2042 ! 0.3*UNR*UNS)**0.5*N0RR(K)*N0S(K)* &
2043 ! (1./(LAMR(K)**3*LAMS(K))+ &
2044 ! 1./(LAMR(K)**2*LAMS(K)**2)+ &
2045 ! 1./(LAMR(K)*LAMS(K)**3))
2049 ! ADD COLLECTION OF GRAUPEL BY RAIN ABOVE FREEZING
2050 ! ASSUME ALL RAIN COLLECTION BY GRAUPEL ABOVE FREEZING IS SHED
2051 ! ASSUME SHED DROPS ARE 1 MM IN SIZE
2053 IF (QR3D(K).GE.1.E-8.AND.QG3D(K).GE.1.E-8) THEN
2055 UMG = AGN(K)*CONS7/(LAMG(K)**BG)
2056 UMR = ARN(K)*CONS4/(LAMR(K)**BR)
2057 UNG = AGN(K)*CONS8/LAMG(K)**BG
2058 UNR = ARN(K)*CONS6/LAMR(K)**BR
2060 ! SET REASLISTIC LIMITS ON FALLSPEEDS
2062 dum=(rhosu/rho(k))**0.54
2063 UMG=MIN(UMG,20.*dum)
2064 UNG=MIN(UNG,20.*dum)
2065 UMR=MIN(UMR,9.1*dum)
2066 UNR=MIN(UNR,9.1*dum)
2068 ! PRACG IS MIXING RATIO OF RAIN PER SEC COLLECTED BY GRAUPEL/HAIL
2069 PRACG(K) = CONS41*(((1.2*UMR-0.95*UMG)**2+ &
2070 0.08*UMG*UMR)**0.5*RHO(K)* &
2071 N0RR(K)*N0G(K)/LAMR(K)**3* &
2072 (5./(LAMR(K)**3*LAMG(K))+ &
2073 2./(LAMR(K)**2*LAMG(K)**2)+ &
2074 0.5/(LAMR(k)*LAMG(k)**3)))
2076 ! ASSUME 1 MM DROPS ARE SHED, GET NUMBER SHED PER SEC
2078 DUM = PRACG(K)/5.2E-7
2080 NPRACG(K) = CONS32*RHO(K)*(1.7*(UNR-UNG)**2+ &
2081 0.3*UNR*UNG)**0.5*N0RR(K)*N0G(K)* &
2082 (1./(LAMR(K)**3*LAMG(K))+ &
2083 1./(LAMR(K)**2*LAMG(K)**2)+ &
2084 1./(LAMR(K)*LAMG(K)**3))
2086 ! hm 7/15/13, remove limit so that the number of collected drops can smaller than
2087 ! number of shed drops
2088 ! NPRACG(K)=MAX(NPRACG(K)-DUM,0.)
2089 NPRACG(K)=NPRACG(K)-DUM
2093 !.......................................................................
2094 ! ACCRETION OF CLOUD LIQUID WATER BY RAIN
2095 ! CONTINUOUS COLLECTION EQUATION WITH
2096 ! GRAVITATIONAL COLLECTION KERNEL, DROPLET FALL SPEED NEGLECTED
2098 IF (QR3D(K).GE.1.E-8 .AND. QC3D(K).GE.1.E-8) THEN
2100 ! 12/13/06 HM ADD, REPLACE WITH NEWER FORMULA FROM
2101 ! KHAIROUTDINOV AND KOGAN 2000, MWR
2103 DUM=(QC3D(K)*QR3D(K))
2104 PRA(K) = 67.*(DUM)**1.15
2105 NPRA(K) = PRA(K)/(QC3D(K)/NC3D(K))
2108 !.......................................................................
2109 ! SELF-COLLECTION OF RAIN DROPS
2111 ! FROM NUMERICAL SIMULATION OF THE STOCHASTIC COLLECTION EQUATION
2112 ! AS DESCRINED ABOVE FOR AUTOCONVERSION
2114 IF (QR3D(K).GE.1.E-8) THEN
2115 ! include breakup add 10/09/09
2117 if (1./lamr(k).lt.dum1) then
2119 else if (1./lamr(k).ge.dum1) then
2120 dum=2.-exp(2300.*(1./lamr(k)-dum1))
2122 ! NRAGG(K) = -8.*NR3D(K)*QR3D(K)*RHO(K)
2123 NRAGG(K) = -5.78*dum*NR3D(K)*QR3D(K)*RHO(K)
2126 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
2127 ! CALCULATE EVAP OF RAIN (RUTLEDGE AND HOBBS 1983)
2129 IF (QR3D(K).GE.QSMALL) THEN
2130 EPSR = 2.*PI*N0RR(K)*RHO(K)*DV(K)* &
2131 (F1R/(LAMR(K)*LAMR(K))+ &
2132 F2R*(ARN(K)*RHO(K)/MU(K))**0.5* &
2133 SC(K)**(1./3.)*CONS9/ &
2139 ! NO CONDENSATION ONTO RAIN, ONLY EVAP ALLOWED
2141 IF (QV3D(K).LT.QVS(K)) THEN
2142 PRE(K) = EPSR*(QV3D(K)-QVS(K))/AB(K)
2143 PRE(K) = MIN(PRE(K),0.)
2148 !.......................................................................
2151 ! SNOW MAY PERSITS ABOVE FREEZING, FORMULA FROM RUTLEDGE AND HOBBS, 1984
2152 ! IF WATER SUPERSATURATION, SNOW MELTS TO FORM RAIN
2154 IF (QNI3D(K).GE.1.E-8) THEN
2157 ! HM, MODIFY FOR V3.2, ADD ACCELERATED MELTING DUE TO COLLISION WITH RAIN
2158 ! DUM = -CPW/XLF(K)*T3D(K)*PRACS(K)
2159 DUM = -CPW/XLF(K)*(T3D(K)-273.15)*PRACS(K)
2162 ! PSMLT(K)=2.*PI*N0S(K)*KAP(K)*(273.15-T3D(K))/ &
2163 ! XLF(K)*RHO(K)*(F1S/(LAMS(K)*LAMS(K))+ &
2164 ! F2S*(ASN(K)*RHO(K)/MU(K))**0.5* &
2165 ! SC(K)**(1./3.)*CONS10/ &
2166 ! (LAMS(K)**CONS35))+DUM
2167 PSMLT(K)=2.*PI*N0S(K)*KAP(K)*(273.15-T3D(K))/ &
2168 XLF(K)*(F1S/(LAMS(K)*LAMS(K))+ &
2169 F2S*(ASN(K)*RHO(K)/MU(K))**0.5* &
2170 SC(K)**(1./3.)*CONS10/ &
2171 (LAMS(K)**CONS35))+DUM
2173 ! IN WATER SUBSATURATION, SNOW MELTS AND EVAPORATES
2175 IF (QVQVS(K).LT.1.) THEN
2176 EPSS = 2.*PI*N0S(K)*RHO(K)*DV(K)* &
2177 (F1S/(LAMS(K)*LAMS(K))+ &
2178 F2S*(ASN(K)*RHO(K)/MU(K))**0.5* &
2179 SC(K)**(1./3.)*CONS10/ &
2182 EVPMS(K) = (QV3D(K)-QVS(K))*EPSS/AB(K)
2183 EVPMS(K) = MAX(EVPMS(K),PSMLT(K))
2184 PSMLT(K) = PSMLT(K)-EVPMS(K)
2188 !.......................................................................
2189 ! MELTING OF GRAUPEL
2191 ! GRAUPEL MAY PERSITS ABOVE FREEZING, FORMULA FROM RUTLEDGE AND HOBBS, 1984
2192 ! IF WATER SUPERSATURATION, GRAUPEL MELTS TO FORM RAIN
2194 IF (QG3D(K).GE.1.E-8) THEN
2197 ! HM, MODIFY FOR V3.2, ADD ACCELERATED MELTING DUE TO COLLISION WITH RAIN
2198 ! DUM = -CPW/XLF(K)*T3D(K)*PRACG(K)
2199 DUM = -CPW/XLF(K)*(T3D(K)-273.15)*PRACG(K)
2202 ! PGMLT(K)=2.*PI*N0G(K)*KAP(K)*(273.15-T3D(K))/ &
2203 ! XLF(K)*RHO(K)*(F1S/(LAMG(K)*LAMG(K))+ &
2204 ! F2S*(AGN(K)*RHO(K)/MU(K))**0.5* &
2205 ! SC(K)**(1./3.)*CONS11/ &
2206 ! (LAMG(K)**CONS36))+DUM
2207 PGMLT(K)=2.*PI*N0G(K)*KAP(K)*(273.15-T3D(K))/ &
2208 XLF(K)*(F1S/(LAMG(K)*LAMG(K))+ &
2209 F2S*(AGN(K)*RHO(K)/MU(K))**0.5* &
2210 SC(K)**(1./3.)*CONS11/ &
2211 (LAMG(K)**CONS36))+DUM
2213 ! IN WATER SUBSATURATION, GRAUPEL MELTS AND EVAPORATES
2215 IF (QVQVS(K).LT.1.) THEN
2216 EPSG = 2.*PI*N0G(K)*RHO(K)*DV(K)* &
2217 (F1S/(LAMG(K)*LAMG(K))+ &
2218 F2S*(AGN(K)*RHO(K)/MU(K))**0.5* &
2219 SC(K)**(1./3.)*CONS11/ &
2222 EVPMG(K) = (QV3D(K)-QVS(K))*EPSG/AB(K)
2223 EVPMG(K) = MAX(EVPMG(K),PGMLT(K))
2224 PGMLT(K) = PGMLT(K)-EVPMG(K)
2229 ! RESET PRACG AND PRACS TO ZERO, THIS IS DONE BECAUSE THERE IS NO
2230 ! TRANSFER OF MASS FROM SNOW AND GRAUPEL TO RAIN DIRECTLY FROM COLLECTION
2231 ! ABOVE FREEZING, IT IS ONLY USED FOR ENHANCEMENT OF MELTING AND SHEDDING
2236 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
2238 ! FOR CLOUD ICE, ONLY PROCESSES OPERATING AT T > 273.15 IS
2239 ! MELTING, WHICH IS ALREADY CONSERVED DURING PROCESS
2242 ! CONSERVATION OF QC
2244 DUM = (PRC(K)+PRA(K))*DT
2246 IF (DUM.GT.QC3D(K).AND.QC3D(K).GE.QSMALL) THEN
2250 PRC(K) = PRC(K)*RATIO
2251 PRA(K) = PRA(K)*RATIO
2255 ! CONSERVATION OF SNOW
2257 DUM = (-PSMLT(K)-EVPMS(K)+PRACS(K))*DT
2259 IF (DUM.GT.QNI3D(K).AND.QNI3D(K).GE.QSMALL) THEN
2261 ! NO SOURCE TERMS FOR SNOW AT T > FREEZING
2262 RATIO = QNI3D(K)/DUM
2264 PSMLT(K) = PSMLT(K)*RATIO
2265 EVPMS(K) = EVPMS(K)*RATIO
2266 PRACS(K) = PRACS(K)*RATIO
2270 ! CONSERVATION OF GRAUPEL
2272 DUM = (-PGMLT(K)-EVPMG(K)+PRACG(K))*DT
2274 IF (DUM.GT.QG3D(K).AND.QG3D(K).GE.QSMALL) THEN
2276 ! NO SOURCE TERM FOR GRAUPEL ABOVE FREEZING
2279 PGMLT(K) = PGMLT(K)*RATIO
2280 EVPMG(K) = EVPMG(K)*RATIO
2281 PRACG(K) = PRACG(K)*RATIO
2285 ! CONSERVATION OF QR
2286 ! HM 12/13/06, ADDED CONSERVATION OF RAIN SINCE PRE IS NEGATIVE
2288 DUM = (-PRACS(K)-PRACG(K)-PRE(K)-PRA(K)-PRC(K)+PSMLT(K)+PGMLT(K))*DT
2290 IF (DUM.GT.QR3D(K).AND.QR3D(K).GE.QSMALL) THEN
2292 RATIO = (QR3D(K)/DT+PRACS(K)+PRACG(K)+PRA(K)+PRC(K)-PSMLT(K)-PGMLT(K))/ &
2294 PRE(K) = PRE(K)*RATIO
2298 !....................................
2300 QV3DTEN(K) = QV3DTEN(K)+(-PRE(K)-EVPMS(K)-EVPMG(K))
2302 T3DTEN(K) = T3DTEN(K)+(PRE(K)*XXLV(K)+(EVPMS(K)+EVPMG(K))*XXLS(K)+&
2303 (PSMLT(K)+PGMLT(K)-PRACS(K)-PRACG(K))*XLF(K))/CPM(K)
2305 QC3DTEN(K) = QC3DTEN(K)+(-PRA(K)-PRC(K))
2306 QR3DTEN(K) = QR3DTEN(K)+(PRE(K)+PRA(K)+PRC(K)-PSMLT(K)-PGMLT(K)+PRACS(K)+PRACG(K))
2307 QNI3DTEN(K) = QNI3DTEN(K)+(PSMLT(K)+EVPMS(K)-PRACS(K))
2308 QG3DTEN(K) = QG3DTEN(K)+(PGMLT(K)+EVPMG(K)-PRACG(K))
2310 ! NS3DTEN(K) = NS3DTEN(K)-NPRACS(K)
2311 ! HM, bug fix 5/12/08, npracg is subtracted from nr not ng
2312 ! NG3DTEN(K) = NG3DTEN(K)
2313 NC3DTEN(K) = NC3DTEN(K)+ (-NPRA(K)-NPRC(K))
2314 NR3DTEN(K) = NR3DTEN(K)+ (NPRC1(K)+NRAGG(K)-NPRACG(K))
2316 ! HM ADD, WRF-CHEM, ADD TENDENCIES FOR C2PREC
2318 C2PREC(K) = PRA(K)+PRC(K)
2319 IF (PRE(K).LT.0.) THEN
2320 DUM = PRE(K)*DT/QR3D(K)
2322 NSUBR(K) = DUM*NR3D(K)/DT
2325 IF (EVPMS(K)+PSMLT(K).LT.0.) THEN
2326 DUM = (EVPMS(K)+PSMLT(K))*DT/QNI3D(K)
2328 NSMLTS(K) = DUM*NS3D(K)/DT
2330 IF (PSMLT(K).LT.0.) THEN
2331 DUM = PSMLT(K)*DT/QNI3D(K)
2333 NSMLTR(K) = DUM*NS3D(K)/DT
2335 IF (EVPMG(K)+PGMLT(K).LT.0.) THEN
2336 DUM = (EVPMG(K)+PGMLT(K))*DT/QG3D(K)
2338 NGMLTG(K) = DUM*NG3D(K)/DT
2340 IF (PGMLT(K).LT.0.) THEN
2341 DUM = PGMLT(K)*DT/QG3D(K)
2343 NGMLTR(K) = DUM*NG3D(K)/DT
2346 NS3DTEN(K) = NS3DTEN(K)+(NSMLTS(K))
2347 NG3DTEN(K) = NG3DTEN(K)+(NGMLTG(K))
2348 NR3DTEN(K) = NR3DTEN(K)+(NSUBR(K)-NSMLTR(K)-NGMLTR(K))
2352 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
2353 ! NOW CALCULATE SATURATION ADJUSTMENT TO CONDENSE EXTRA VAPOR ABOVE
2356 DUMT = T3D(K)+DT*T3DTEN(K)
2357 DUMQV = QV3D(K)+DT*QV3DTEN(K)
2358 ! hm, add fix for low pressure, 5/12/10
2359 dum=min(0.99*pres(k),POLYSVP(DUMT,0))
2360 DUMQSS = EP_2*dum/(PRES(K)-dum)
2361 DUMQC = QC3D(K)+DT*QC3DTEN(K)
2362 DUMQC = MAX(DUMQC,0.)
2364 ! SATURATION ADJUSTMENT FOR LIQUID
2367 PCC(K) = DUMS/(1.+XXLV(K)**2*DUMQSS/(CPM(K)*RV*DUMT**2))/DT
2368 IF (PCC(K)*DT+DUMQC.LT.0.) THEN
2372 QV3DTEN(K) = QV3DTEN(K)-PCC(K)
2373 T3DTEN(K) = T3DTEN(K)+PCC(K)*XXLV(K)/CPM(K)
2374 QC3DTEN(K) = QC3DTEN(K)+PCC(K)
2377 IF ( has_wetscav ) THEN
2378 evapprod(k) = - PRE(K) - EVPMS(K) - EVPMG(K)
2379 rainprod(k) = PRA(K) + PRC(K) + tqimelt(K)
2383 !.......................................................................
2384 ! ACTIVATION OF CLOUD DROPLETS
2385 ! DROPLET CONCENTRATION IS SPECIFIED !!!!!
2387 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
2388 !!TWG/amy begin Morrison Activation Code
2391 IF (QC3D(K)+QC3DTEN(K)*DT.GE.QSMALL) THEN
2393 ! EFFECTIVE VERTICAL VELOCITY (M/S)
2396 ! ADD SUB-GRID VERTICAL VELOCITY
2397 DUM = W3D(K)+WVAR(K)
2399 ! ASSUME MINIMUM EFF. SUB-GRID VELOCITY 0.10 M/S
2402 ELSE IF (ISUB.EQ.1) THEN
2405 ! ONLY ACTIVATE IN REGIONS OF UPWARD MOTION
2406 IF (DUM.GE.0.001) THEN
2408 IF (IBASE.EQ.1) THEN
2410 ! ACTIVATE ONLY IF THERE IS LITTLE CLOUD WATER
2411 ! OR IF AT CLOUD BASE, OR AT LOWEST MODEL LEVEL (K=1)
2415 IF (QC3D(K)+QC3DTEN(K)*DT.LE.0.05E-3/RHO(K)) THEN
2420 ELSE IF (K.GE.2) THEN
2421 IF (QC3D(K)+QC3DTEN(K)*DT.GT.0.05E-3/RHO(K).AND. &
2422 QC3D(K-1)+QC3DTEN(K-1)*DT.LE.0.05E-3/RHO(K-1)) THEN
2427 IF (IDROP.EQ.1) THEN
2428 ! ACTIVATE AT CLOUD BASE OR REGIONS WITH VERY LITTLE LIQ WATER
2431 ! USE ROGERS AND YAU (1989) TO RELATE NUMBER ACTIVATED TO W
2432 ! BASED ON TWOMEY 1959
2434 DUM=DUM*100. ! CONVERT FROM M/S TO CM/S
2435 DUM2 = 0.88*C1**(2./(K1+2.))*(7.E-2*DUM**1.5)**(K1/(K1+2.))
2436 DUM2=DUM2*1.E6 ! CONVERT FROM CM-3 TO M-3
2437 DUM2=DUM2/RHO(K) ! CONVERT FROM M-3 TO KG-1
2438 DUM2 = (DUM2-NC3D(K))/DT
2440 NC3DTEN(K) = NC3DTEN(K)+DUM2
2442 ELSE IF (IACT.EQ.2) THEN
2443 ! DROPLET ACTIVATION FROM ABDUL-RAZZAK AND GHAN (2000)
2445 SIGVL = 0.0761-1.55E-4*(T3D(K)-273.15)
2446 AACT = 2.*MW/(RHOW*RR)*SIGVL/T3D(K)
2447 ALPHA = G*MW*XXLV(K)/(CPM(K)*RR*T3D(K)**2)-G*MA/(RR*T3D(K))
2448 GAMM = RR*T3D(K)/(EVS(K)*MW)+MW*XXLV(K)**2/(CPM(K)*PRES(K)*MA*T3D(K))
2450 GG =1./(RHOW*RR*T3D(K)/(EVS(K)*DV(K)*MW)+XXLV(K)*RHOW/(KAP(K)*T3D(K))*(XXLV(K)*MW/ &
2453 PSI = 2./3.*(ALPHA*DUM/GG)**0.5*AACT
2455 ETA1 = (ALPHA*DUM/GG)**1.5/(2.*PI*RHOW*GAMM*NANEW1)
2456 ETA2 = (ALPHA*DUM/GG)**1.5/(2.*PI*RHOW*GAMM*NANEW2)
2458 SM1 = 2./BACT**0.5*(AACT/(3.*RM1))**1.5
2459 SM2 = 2./BACT**0.5*(AACT/(3.*RM2))**1.5
2461 DUM1=1./SM1**2*(F11*(PSI/ETA1)**1.5+F21*(SM1**2/(ETA1+3.*PSI))**0.75)
2462 DUM2=1./SM2**2*(F12*(PSI/ETA2)**1.5+F22*(SM2**2/(ETA2+3.*PSI))**0.75)
2464 SMAX = 1./(DUM1+DUM2)**0.5
2466 UU1 = 2.*LOG(SM1/SMAX)/(4.242*LOG(SIG1))
2467 UU2 = 2.*LOG(SM2/SMAX)/(4.242*LOG(SIG2))
2468 DUM1 = NANEW1/2.*(1.-DERF1(UU1))
2469 DUM2 = NANEW2/2.*(1.-DERF1(UU2))
2471 DUM2 = (DUM1+DUM2)/RHO(K) !CONVERT TO KG-1
2473 ! MAKE SURE THIS VALUE ISN'T GREATER THAN TOTAL NUMBER OF AEROSOL
2474 DUM2 = MIN((NANEW1+NANEW2)/RHO(K),DUM2)
2476 DUM2 = (DUM2-NC3D(K))/DT
2478 NC3DTEN(K) = NC3DTEN(K)+DUM2
2481 ELSE IF (IACT.EQ.4) THEN
2482 DUM = W3D(K)+WVAR(K)
2483 call mdm_prescribed_activate(DUM,T3D(K),RHO(K), &
2484 naero, naer_cu,naer_cu, maero, &
2485 dispersion_aer,hygro_aer, density_aer, DUM2, XXLV(K))
2487 DUM2 = (DUM2-NC3D(K))/DT
2489 NC3DTEN(K) = NC3DTEN(K)+DUM2
2493 !.............................................................................
2494 ELSE IF (IDROP.EQ.0) THEN
2495 ! ACTIVATE IN CLOUD INTERIOR
2496 ! FIND EQUILIBRIUM SUPERSATURATION
2498 TAUC=1./(2.*PI*RHO(k)*DV(K)*NC3D(K)*(PGAM(K)+1.)/LAMC(K))
2499 IF (EPSR.GT.1.E-8) THEN
2506 ! DUM3=(QVS(K)*RHO(K)/(PRES(K)-EVS(K))+DQSDT/CP)*G*DUM
2507 DUM3=(-QVS(K)*RHO(K)/(PRES(K)-EVS(K))+DQSDT/CP)*G*DUM
2508 DUM3=DUM3*TAUC*TAUR/(TAUC+TAUR)
2510 IF (DUM3/QVS(K).GE.1.E-6) THEN
2513 ! FIND MAXIMUM ALLOWED ACTIVATION WITH NON-EQUILIBRIUM SS
2515 DUM=DUM*100. ! CONVERT FROM M/S TO CM/S
2516 DUMACT = 0.88*C1**(2./(K1+2.))*(7.E-2*DUM**1.5)**(K1/(K1+2.))
2518 ! USE POWER LAW CCN SPECTRA
2520 ! CONVERT FROM ABSOLUTE SUPERSATURATION TO SUPERSATURATION RATIO IN %
2521 DUM3=DUM3/QVS(K)*100.
2524 ! MAKE SURE VALUE DOESN'T EXCEED THAT FOR NON-EQUILIBRIUM SS
2525 DUM2=MIN(DUM2,DUMACT)
2526 DUM2=DUM2*1.E6 ! CONVERT FROM CM-3 TO M-3
2527 DUM2=DUM2/RHO(K) ! CONVERT FROM M-3 TO KG-1
2528 DUM2 = (DUM2-NC3D(K))/DT
2530 NC3DTEN(K) = NC3DTEN(K)+DUM2
2532 ELSE IF (IACT.EQ.2) THEN
2534 ! FIND MAXIMUM ALLOWED ACTIVATION WITH NON-EQUILIBRIUM SS
2536 SIGVL = 0.0761-1.55E-4*(T3D(K)-273.15)
2537 AACT = 2.*MW/(RHOW*RR)*SIGVL/T3D(K)
2538 ! ALPHA = G*MW*XXLV(K)/(CPM(K)*RR*T3D(K)**2)-G*MA/(RR*T3D(K))
2540 ! RR*T3D(K)/(EVS(K)*MW)+MW*XXLV(K)**2/(CPM(K)*PRES(K)*MA*T3D(K))
2542 GG =1./(RHOW*RR*T3D(K)/(EVS(K)*DV(K)*MW)+XXLV(K)*RHOW/(KAP(K)*T3D(K))*(XXLV(K)*MW/ &
2545 PSI = 2./3.*(ALPHA*DUM/GG)**0.5*AACT
2547 ETA1 = (ALPHA*DUM/GG)**1.5/(2.*PI*RHOW*GAMM*NANEW1)
2548 ETA2 = (ALPHA*DUM/GG)**1.5/(2.*PI*RHOW*GAMM*NANEW2)
2550 SM1 = 2./BACT**0.5*(AACT/(3.*RM1))**1.5
2551 SM2 = 2./BACT**0.5*(AACT/(3.*RM2))**1.5
2553 DUM1=1./SM1**2*(F11*(PSI/ETA1)**1.5+F21*(SM1**2/(ETA1+3.*PSI))**0.75)
2554 DUM2=1./SM2**2*(F12*(PSI/ETA2)**1.5+F22*(SM2**2/(ETA2+3.*PSI))**0.75)
2556 SMAX = 1./(DUM1+DUM2)**0.5
2558 UU1 = 2.*LOG(SM1/SMAX)/(4.242*LOG(SIG1))
2559 UU2 = 2.*LOG(SM2/SMAX)/(4.242*LOG(SIG2))
2560 DUM1 = NANEW1/2.*(1.-DERF1(UU1))
2561 DUM2 = NANEW2/2.*(1.-DERF1(UU2))
2563 DUM2 = (DUM1+DUM2)/RHO(K) !CONVERT TO KG-1
2565 ! MAKE SURE THIS VALUE ISN'T GREATER THAN TOTAL NUMBER OF AEROSOL
2567 DUM2 = (DUM1+DUM2)/RHO(K) !CONVERT TO KG-1
2569 ! MAKE SURE THIS VALUE ISN'T GREATER THAN TOTAL NUMBER OF AEROSOL
2571 DUMACT = MIN((NANEW1+NANEW2)/RHO(K),DUM2)
2573 ! USE LOGNORMAL AEROSOL
2574 SIGVL = 0.0761-1.55E-4*(T3D(K)-273.15)
2575 AACT = 2.*MW/(RHOW*RR)*SIGVL/T3D(K)
2577 SM1 = 2./BACT**0.5*(AACT/(3.*RM1))**1.5
2578 SM2 = 2./BACT**0.5*(AACT/(3.*RM2))**1.5
2580 ! GET SUPERSATURATION RATIO FROM ABSOLUTE SUPERSATURATION
2583 UU1 = 2.*LOG(SM1/SMAX)/(4.242*LOG(SIG1))
2584 UU2 = 2.*LOG(SM2/SMAX)/(4.242*LOG(SIG2))
2585 DUM1 = NANEW1/2.*(1.-DERF1(UU1))
2586 DUM2 = NANEW2/2.*(1.-DERF1(UU2))
2588 DUM2 = (DUM1+DUM2)/RHO(K) !CONVERT TO KG-1
2590 ! MAKE SURE THIS VALUE ISN'T GREATER THAN TOTAL NUMBER OF AEROSOL
2592 DUM2 = MIN((NANEW1+NANEW2)/RHO(K),DUM2)
2594 ! MAKE SURE ISN'T GREATER THAN NON-EQUIL. SS
2595 DUM2=MIN(DUM2,DUMACT)
2597 DUM2 = (DUM2-NC3D(K))/DT
2599 NC3DTEN(K) = NC3DTEN(K)+DUM2
2601 ELSE IF (IACT.EQ.4) THEN
2602 DUM = W3D(K)+WVAR(K)
2603 call mdm_prescribed_activate(DUM,T3D(K),RHO(K), &
2604 naero, naer_cu,naer_cu, maero, &
2605 dispersion_aer,hygro_aer, density_aer, DUM2, XXLV(K))
2607 DUM2 = (DUM2-NC3D(K))/DT
2609 NC3DTEN(K) = NC3DTEN(K)+DUM2
2612 END IF ! DUM3/QVS > 1.E-6
2615 !.......................................................................
2616 ELSE IF (IBASE.EQ.2) THEN
2619 ! USE ROGERS AND YAU (1989) TO RELATE NUMBER ACTIVATED TO W
2620 ! BASED ON TWOMEY 1959
2622 DUM=DUM*100. ! CONVERT FROM M/S TO CM/S
2623 DUM2 = 0.88*C1**(2./(K1+2.))*(7.E-2*DUM**1.5)**(K1/(K1+2.))
2624 DUM2=DUM2*1.E6 ! CONVERT FROM CM-3 TO M-3
2625 DUM2=DUM2/RHO(K) ! CONVERT FROM M-3 TO KG-1
2626 DUM2 = (DUM2-NC3D(K))/DT
2628 NC3DTEN(K) = NC3DTEN(K)+DUM2
2630 ELSE IF (IACT.EQ.2) THEN
2632 SIGVL = 0.0761-1.55E-4*(T3D(K)-273.15)
2633 AACT = 2.*MW/(RHOW*RR)*SIGVL/T3D(K)
2634 ALPHA = G*MW*XXLV(K)/(CPM(K)*RR*T3D(K)**2)-G*MA/(RR*T3D(K))
2635 GAMM = RR*T3D(K)/(EVS(K)*MW)+MW*XXLV(K)**2/(CPM(K)*PRES(K)*MA*T3D(K))
2637 GG =1./(RHOW*RR*T3D(K)/(EVS(K)*DV(K)*MW)+XXLV(K)*RHOW/(KAP(K)*T3D(K))*(XXLV(K)*MW/ &
2640 PSI = 2./3.*(ALPHA*DUM/GG)**0.5*AACT
2642 ETA1 = (ALPHA*DUM/GG)**1.5/(2.*PI*RHOW*GAMM*NANEW1)
2643 ETA2 = (ALPHA*DUM/GG)**1.5/(2.*PI*RHOW*GAMM*NANEW2)
2645 SM1 = 2./BACT**0.5*(AACT/(3.*RM1))**1.5
2646 SM2 = 2./BACT**0.5*(AACT/(3.*RM2))**1.5
2648 DUM1=1./SM1**2*(F11*(PSI/ETA1)**1.5+F21*(SM1**2/(ETA1+3.*PSI))**0.75)
2649 DUM2=1./SM2**2*(F12*(PSI/ETA2)**1.5+F22*(SM2**2/(ETA2+3.*PSI))**0.75)
2651 SMAX = 1./(DUM1+DUM2)**0.5
2653 UU1 = 2.*LOG(SM1/SMAX)/(4.242*LOG(SIG1))
2654 UU2 = 2.*LOG(SM2/SMAX)/(4.242*LOG(SIG2))
2655 DUM1 = NANEW1/2.*(1.-DERF1(UU1))
2656 DUM2 = NANEW2/2.*(1.-DERF1(UU2))
2658 DUM2 = (DUM1+DUM2)/RHO(K) !CONVERT TO KG-1
2660 ! MAKE SURE THIS VALUE ISN'T GREATER THAN TOTAL NUMBER OF AEROSOL
2662 DUM2 = MIN((NANEW1+NANEW2)/RHO(K),DUM2)
2664 DUM2 = (DUM2-NC3D(K))/DT
2666 NC3DTEN(K) = NC3DTEN(K)+DUM2
2669 ELSE IF (IACT.EQ.4) THEN
2670 DUM = W3D(K)+WVAR(K)
2671 call mdm_prescribed_activate(DUM,T3D(K),RHO(K), &
2672 naero, naer_cu,naer_cu, maero, &
2673 dispersion_aer,hygro_aer, density_aer, DUM2, XXLV(K))
2675 DUM2 = (DUM2-NC3D(K))/DT
2677 NC3DTEN(K) = NC3DTEN(K)+DUM2
2682 END IF ! QC3D > QSMALL
2688 ! SUBLIMATE, MELT, OR EVAPORATE NUMBER CONCENTRATION
2689 ! THIS FORMULATION ASSUMES 1:1 RATIO BETWEEN MASS LOSS AND
2690 ! LOSS OF NUMBER CONCENTRATION
2692 ! IF (PCC(K).LT.0.) THEN
2693 ! DUM = PCC(K)*DT/QC3D(K)
2694 ! DUM = MAX(-1.,DUM)
2695 ! NSUBC(K) = DUM*NC3D(K)/DT
2700 ! NC3DTEN(K) = NC3DTEN(K)+NSUBC(K)
2702 !.....................................................................
2703 !.....................................................................
2704 ELSE ! TEMPERATURE < 273.15
2706 !......................................................................
2707 !HM ADD, ALLOW FOR CONSTANT DROPLET NUMBER
2708 ! INUM = 0, PREDICT DROPLET NUMBER
2709 ! INUM = 1, SET CONSTANT DROPLET NUMBER
2711 IF (iinum.EQ.1) THEN
2712 ! CONVERT NDCNST FROM CM-3 TO KG-1
2713 NC3D(K)=NDCNST*1.E6/RHO(K)
2716 ! CALCULATE SIZE DISTRIBUTION PARAMETERS
2717 ! MAKE SURE NUMBER CONCENTRATIONS AREN'T NEGATIVE
2719 NI3D(K) = MAX(0.,NI3D(K))
2720 NS3D(K) = MAX(0.,NS3D(K))
2721 NC3D(K) = MAX(0.,NC3D(K))
2722 NR3D(K) = MAX(0.,NR3D(K))
2723 NG3D(K) = MAX(0.,NG3D(K))
2725 !......................................................................
2728 IF (QI3D(K).GE.QSMALL) THEN
2729 LAMI(K) = (CONS12* &
2730 NI3D(K)/QI3D(K))**(1./DI)
2731 N0I(K) = NI3D(K)*LAMI(K)
2737 IF (LAMI(K).LT.LAMMINI) THEN
2741 N0I(K) = LAMI(K)**4*QI3D(K)/CONS12
2743 NI3D(K) = N0I(K)/LAMI(K)
2744 ELSE IF (LAMI(K).GT.LAMMAXI) THEN
2746 N0I(K) = LAMI(K)**4*QI3D(K)/CONS12
2748 NI3D(K) = N0I(K)/LAMI(K)
2752 !......................................................................
2755 IF (QR3D(K).GE.QSMALL) THEN
2756 LAMR(K) = (PI*RHOW*NR3D(K)/QR3D(K))**(1./3.)
2757 N0RR(K) = NR3D(K)*LAMR(K)
2763 IF (LAMR(K).LT.LAMMINR) THEN
2767 N0RR(K) = LAMR(K)**4*QR3D(K)/(PI*RHOW)
2769 NR3D(K) = N0RR(K)/LAMR(K)
2770 ELSE IF (LAMR(K).GT.LAMMAXR) THEN
2772 N0RR(K) = LAMR(K)**4*QR3D(K)/(PI*RHOW)
2774 NR3D(K) = N0RR(K)/LAMR(K)
2778 !......................................................................
2781 ! MARTIN ET AL. (1994) FORMULA FOR PGAM
2783 IF (QC3D(K).GE.QSMALL) THEN
2785 DUM = PRES(K)/(287.15*T3D(K))
2786 PGAM(K)=0.0005714*(NC3D(K)/1.E6*DUM)+0.2714
2787 PGAM(K)=1./(PGAM(K)**2)-1.
2788 PGAM(K)=MAX(PGAM(K),2.)
2789 PGAM(K)=MIN(PGAM(K),10.)
2793 LAMC(K) = (CONS26*NC3D(K)*GAMMA(PGAM(K)+4.)/ &
2794 (QC3D(K)*GAMMA(PGAM(K)+1.)))**(1./3.)
2796 ! LAMMIN, 60 MICRON DIAMETER
2799 LAMMIN = (PGAM(K)+1.)/60.E-6
2800 LAMMAX = (PGAM(K)+1.)/1.E-6
2802 IF (LAMC(K).LT.LAMMIN) THEN
2805 NC3D(K) = EXP(3.*LOG(LAMC(K))+LOG(QC3D(K))+ &
2806 LOG(GAMMA(PGAM(K)+1.))-LOG(GAMMA(PGAM(K)+4.)))/CONS26
2807 ELSE IF (LAMC(K).GT.LAMMAX) THEN
2809 NC3D(K) = EXP(3.*LOG(LAMC(K))+LOG(QC3D(K))+ &
2810 LOG(GAMMA(PGAM(K)+1.))-LOG(GAMMA(PGAM(K)+4.)))/CONS26
2814 ! TO CALCULATE DROPLET FREEZING
2816 CDIST1(K) = NC3D(K)/GAMMA(PGAM(K)+1.)
2820 !......................................................................
2823 IF (QNI3D(K).GE.QSMALL) THEN
2824 LAMS(K) = (CONS1*NS3D(K)/QNI3D(K))**(1./DS)
2825 N0S(K) = NS3D(K)*LAMS(K)
2831 IF (LAMS(K).LT.LAMMINS) THEN
2833 N0S(K) = LAMS(K)**4*QNI3D(K)/CONS1
2835 NS3D(K) = N0S(K)/LAMS(K)
2837 ELSE IF (LAMS(K).GT.LAMMAXS) THEN
2840 N0S(K) = LAMS(K)**4*QNI3D(K)/CONS1
2842 NS3D(K) = N0S(K)/LAMS(K)
2846 !......................................................................
2849 IF (QG3D(K).GE.QSMALL) THEN
2850 LAMG(K) = (CONS2*NG3D(K)/QG3D(K))**(1./DG)
2851 N0G(K) = NG3D(K)*LAMG(K)
2857 IF (LAMG(K).LT.LAMMING) THEN
2859 N0G(K) = LAMG(K)**4*QG3D(K)/CONS2
2861 NG3D(K) = N0G(K)/LAMG(K)
2863 ELSE IF (LAMG(K).GT.LAMMAXG) THEN
2866 N0G(K) = LAMG(K)**4*QG3D(K)/CONS2
2868 NG3D(K) = N0G(K)/LAMG(K)
2872 !.....................................................................
2873 ! ZERO OUT PROCESS RATES
2922 ! HM: ADD GRAUPEL PROCESSES
2936 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
2937 ! CALCULATION OF MICROPHYSICAL PROCESS RATES
2938 ! ACCRETION/AUTOCONVERSION/FREEZING/MELTING/COAG.
2939 !.......................................................................
2940 ! FREEZING OF CLOUD DROPLETS
2941 ! ONLY ALLOWED BELOW -4 C
2942 IF (QC3D(K).GE.QSMALL .AND. T3D(K).LT.269.15) THEN
2944 ! NUMBER OF CONTACT NUCLEI (M^-3) FROM MEYERS ET AL., 1992
2945 ! FACTOR OF 1000 IS TO CONVERT FROM L^-1 TO M^-3
2949 NACNT = EXP(-2.80+0.262*(273.15-T3D(K)))*1000.
2952 ! NACNT = 5.*EXP(0.304*(273.15-T3D(K)))
2955 ! NACNT = 0.01*EXP(0.6*(273.15-T3D(K)))
2961 DUM = 7.37*T3D(K)/(288.*10.*PRES(K))/100.
2963 ! EFFECTIVE DIFFUSIVITY OF CONTACT NUCLEI
2964 ! BASED ON BROWNIAN DIFFUSION
2966 DAP(K) = CONS37*T3D(K)*(1.+DUM/RIN)/MU(K)
2968 MNUCCC(K) = CONS38*DAP(K)*NACNT*EXP(LOG(CDIST1(K))+ &
2969 LOG(GAMMA(PGAM(K)+5.))-4.*LOG(LAMC(K)))
2970 NNUCCC(K) = 2.*PI*DAP(K)*NACNT*CDIST1(K)* &
2971 GAMMA(PGAM(K)+2.)/ &
2974 ! IMMERSION FREEZING (BIGG 1953)
2976 ! MNUCCC(K) = MNUCCC(K)+CONS39* &
2977 ! EXP(LOG(CDIST1(K))+LOG(GAMMA(7.+PGAM(K)))-6.*LOG(LAMC(K)))* &
2978 ! EXP(AIMM*(273.15-T3D(K)))
2980 ! NNUCCC(K) = NNUCCC(K)+ &
2981 ! CONS40*EXP(LOG(CDIST1(K))+LOG(GAMMA(PGAM(K)+4.))-3.*LOG(LAMC(K))) &
2982 ! *EXP(AIMM*(273.15-T3D(K)))
2984 ! hm 7/15/13 fix for consistency w/ original formula
2985 MNUCCC(K) = MNUCCC(K)+CONS39* &
2986 EXP(LOG(CDIST1(K))+LOG(GAMMA(7.+PGAM(K)))-6.*LOG(LAMC(K)))* &
2987 (EXP(AIMM*(273.15-T3D(K)))-1.)
2989 NNUCCC(K) = NNUCCC(K)+ &
2990 CONS40*EXP(LOG(CDIST1(K))+LOG(GAMMA(PGAM(K)+4.))-3.*LOG(LAMC(K))) &
2991 *(EXP(AIMM*(273.15-T3D(K)))-1.)
2993 ! PUT IN A CATCH HERE TO PREVENT DIVERGENCE BETWEEN NUMBER CONC. AND
2994 ! MIXING RATIO, SINCE STRICT CONSERVATION NOT CHECKED FOR NUMBER CONC
2996 NNUCCC(K) = MIN(NNUCCC(K),NC3D(K)/DT)
3000 !.................................................................
3001 !.......................................................................
3002 ! AUTOCONVERSION OF CLOUD LIQUID WATER TO RAIN
3003 ! FORMULA FROM BEHENG (1994)
3004 ! USING NUMERICAL SIMULATION OF STOCHASTIC COLLECTION EQUATION
3005 ! AND INITIAL CLOUD DROPLET SIZE DISTRIBUTION SPECIFIED
3006 ! AS A GAMMA DISTRIBUTION
3008 ! USE MINIMUM VALUE OF 1.E-6 TO PREVENT FLOATING POINT ERROR
3010 IF (QC3D(K).GE.1.E-6) THEN
3012 ! HM ADD 12/13/06, REPLACE WITH NEWER FORMULA
3013 ! FROM KHAIROUTDINOV AND KOGAN 2000, MWR
3015 PRC(K)=1350.*QC3D(K)**2.47* &
3016 (NC3D(K)/1.e6*RHO(K))**(-1.79)
3018 ! note: nprc1 is change in Nr,
3019 ! nprc is change in Nc
3021 NPRC1(K) = PRC(K)/CONS29
3022 NPRC(K) = PRC(K)/(QC3D(K)/NC3D(K))
3024 ! hm bug fix 3/20/12
3025 NPRC(K) = MIN(NPRC(K),NC3D(K)/DT)
3026 NPRC1(K) = MIN(NPRC1(K),NPRC(K))
3030 !.......................................................................
3031 ! SELF-COLLECTION OF DROPLET NOT INCLUDED IN KK2000 SCHEME
3033 ! SNOW AGGREGATION FROM PASSARELLI, 1978, USED BY REISNER, 1998
3034 ! THIS IS HARD-WIRED FOR BS = 0.4 FOR NOW
3036 IF (QNI3D(K).GE.1.E-8) THEN
3037 NSAGG(K) = CONS15*ASN(K)*RHO(K)** &
3038 ((2.+BS)/3.)*QNI3D(K)**((2.+BS)/3.)* &
3039 (NS3D(K)*RHO(K))**((4.-BS)/3.)/ &
3043 !.......................................................................
3044 ! ACCRETION OF CLOUD DROPLETS ONTO SNOW/GRAUPEL
3045 ! HERE USE CONTINUOUS COLLECTION EQUATION WITH
3046 ! SIMPLE GRAVITATIONAL COLLECTION KERNEL IGNORING
3050 IF (QNI3D(K).GE.1.E-8 .AND. QC3D(K).GE.QSMALL) THEN
3052 PSACWS(K) = CONS13*ASN(K)*QC3D(K)*RHO(K)* &
3055 NPSACWS(K) = CONS13*ASN(K)*NC3D(K)*RHO(K)* &
3061 !............................................................................
3062 ! COLLECTION OF CLOUD WATER BY GRAUPEL
3064 IF (QG3D(K).GE.1.E-8 .AND. QC3D(K).GE.QSMALL) THEN
3066 PSACWG(K) = CONS14*AGN(K)*QC3D(K)*RHO(K)* &
3069 NPSACWG(K) = CONS14*AGN(K)*NC3D(K)*RHO(K)* &
3074 !.......................................................................
3076 ! CLOUD ICE COLLECTING DROPLETS, ASSUME THAT CLOUD ICE MEAN DIAM > 100 MICRON
3077 ! BEFORE RIMING CAN OCCUR
3078 ! ASSUME THAT RIME COLLECTED ON CLOUD ICE DOES NOT LEAD
3079 ! TO HALLET-MOSSOP SPLINTERING
3081 IF (QI3D(K).GE.1.E-8 .AND. QC3D(K).GE.QSMALL) THEN
3083 ! PUT IN SIZE DEPENDENT COLLECTION EFFICIENCY BASED ON STOKES LAW
3084 ! FROM THOMPSON ET AL. 2004, MWR
3086 IF (1./LAMI(K).GE.100.E-6) THEN
3088 PSACWI(K) = CONS16*AIN(K)*QC3D(K)*RHO(K)* &
3091 NPSACWI(K) = CONS16*AIN(K)*NC3D(K)*RHO(K)* &
3097 !.......................................................................
3098 ! ACCRETION OF RAIN WATER BY SNOW
3099 ! FORMULA FROM IKAWA AND SAITO, 1991, USED BY REISNER ET AL, 1998
3101 IF (QR3D(K).GE.1.E-8.AND.QNI3D(K).GE.1.E-8) THEN
3103 UMS = ASN(K)*CONS3/(LAMS(K)**BS)
3104 UMR = ARN(K)*CONS4/(LAMR(K)**BR)
3105 UNS = ASN(K)*CONS5/LAMS(K)**BS
3106 UNR = ARN(K)*CONS6/LAMR(K)**BR
3108 ! SET REASLISTIC LIMITS ON FALLSPEEDS
3111 dum=(rhosu/rho(k))**0.54
3112 UMS=MIN(UMS,1.2*dum)
3113 UNS=MIN(UNS,1.2*dum)
3114 UMR=MIN(UMR,9.1*dum)
3115 UNR=MIN(UNR,9.1*dum)
3117 PRACS(K) = CONS41*(((1.2*UMR-0.95*UMS)**2+ &
3118 0.08*UMS*UMR)**0.5*RHO(K)* &
3119 N0RR(K)*N0S(K)/LAMR(K)**3* &
3120 (5./(LAMR(K)**3*LAMS(K))+ &
3121 2./(LAMR(K)**2*LAMS(K)**2)+ &
3122 0.5/(LAMR(k)*LAMS(k)**3)))
3124 NPRACS(K) = CONS32*RHO(K)*(1.7*(UNR-UNS)**2+ &
3125 0.3*UNR*UNS)**0.5*N0RR(K)*N0S(K)* &
3126 (1./(LAMR(K)**3*LAMS(K))+ &
3127 1./(LAMR(K)**2*LAMS(K)**2)+ &
3128 1./(LAMR(K)*LAMS(K)**3))
3130 ! MAKE SURE PRACS DOESN'T EXCEED TOTAL RAIN MIXING RATIO
3131 ! AS THIS MAY OTHERWISE RESULT IN TOO MUCH TRANSFER OF WATER DURING
3134 PRACS(K) = MIN(PRACS(K),QR3D(K)/DT)
3136 ! COLLECTION OF SNOW BY RAIN - NEEDED FOR GRAUPEL CONVERSION CALCULATIONS
3137 ! ONLY CALCULATE IF SNOW AND RAIN MIXING RATIOS EXCEED 0.1 G/KG
3139 ! HM MODIFY FOR WRFV3.1
3140 ! IF (IHAIL.EQ.0) THEN
3141 IF (QNI3D(K).GE.0.1E-3.AND.QR3D(K).GE.0.1E-3) THEN
3142 PSACR(K) = CONS31*(((1.2*UMR-0.95*UMS)**2+ &
3143 0.08*UMS*UMR)**0.5*RHO(K)* &
3144 N0RR(K)*N0S(K)/LAMS(K)**3* &
3145 (5./(LAMS(K)**3*LAMR(K))+ &
3146 2./(LAMS(K)**2*LAMR(K)**2)+ &
3147 0.5/(LAMS(K)*LAMR(K)**3)))
3153 !.......................................................................
3155 ! COLLECTION OF RAINWATER BY GRAUPEL, FROM IKAWA AND SAITO 1990,
3156 ! USED BY REISNER ET AL 1998
3157 IF (QR3D(K).GE.1.E-8.AND.QG3D(K).GE.1.E-8) THEN
3159 UMG = AGN(K)*CONS7/(LAMG(K)**BG)
3160 UMR = ARN(K)*CONS4/(LAMR(K)**BR)
3161 UNG = AGN(K)*CONS8/LAMG(K)**BG
3162 UNR = ARN(K)*CONS6/LAMR(K)**BR
3164 ! SET REASLISTIC LIMITS ON FALLSPEEDS
3166 dum=(rhosu/rho(k))**0.54
3167 UMG=MIN(UMG,20.*dum)
3168 UNG=MIN(UNG,20.*dum)
3169 UMR=MIN(UMR,9.1*dum)
3170 UNR=MIN(UNR,9.1*dum)
3172 PRACG(K) = CONS41*(((1.2*UMR-0.95*UMG)**2+ &
3173 0.08*UMG*UMR)**0.5*RHO(K)* &
3174 N0RR(K)*N0G(K)/LAMR(K)**3* &
3175 (5./(LAMR(K)**3*LAMG(K))+ &
3176 2./(LAMR(K)**2*LAMG(K)**2)+ &
3177 0.5/(LAMR(k)*LAMG(k)**3)))
3179 NPRACG(K) = CONS32*RHO(K)*(1.7*(UNR-UNG)**2+ &
3180 0.3*UNR*UNG)**0.5*N0RR(K)*N0G(K)* &
3181 (1./(LAMR(K)**3*LAMG(K))+ &
3182 1./(LAMR(K)**2*LAMG(K)**2)+ &
3183 1./(LAMR(K)*LAMG(K)**3))
3185 ! MAKE SURE PRACG DOESN'T EXCEED TOTAL RAIN MIXING RATIO
3186 ! AS THIS MAY OTHERWISE RESULT IN TOO MUCH TRANSFER OF WATER DURING
3189 PRACG(K) = MIN(PRACG(K),QR3D(K)/DT)
3193 !.......................................................................
3194 ! RIME-SPLINTERING - SNOW
3195 ! HALLET-MOSSOP (1974)
3196 ! NUMBER OF SPLINTERS FORMED IS BASED ON MASS OF RIMED WATER
3198 ! DUM1 = MASS OF INDIVIDUAL SPLINTERS
3200 ! HM ADD THRESHOLD SNOW AND DROPLET MIXING RATIO FOR RIME-SPLINTERING
3201 ! TO LIMIT RIME-SPLINTERING IN STRATIFORM CLOUDS
3202 ! THESE THRESHOLDS CORRESPOND WITH GRAUPEL THRESHOLDS IN RH 1984
3205 IF (QNI3D(K).GE.0.1E-3) THEN
3206 IF (QC3D(K).GE.0.5E-3.OR.QR3D(K).GE.0.1E-3) THEN
3207 IF (PSACWS(K).GT.0..OR.PRACS(K).GT.0.) THEN
3208 IF (T3D(K).LT.270.16 .AND. T3D(K).GT.265.16) THEN
3210 IF (T3D(K).GT.270.16) THEN
3212 ELSE IF (T3D(K).LE.270.16.AND.T3D(K).GT.268.16) THEN
3213 FMULT = (270.16-T3D(K))/2.
3214 ELSE IF (T3D(K).GE.265.16.AND.T3D(K).LE.268.16) THEN
3215 FMULT = (T3D(K)-265.16)/3.
3216 ELSE IF (T3D(K).LT.265.16) THEN
3220 ! 1000 IS TO CONVERT FROM KG TO G
3222 ! SPLINTERING FROM DROPLETS ACCRETED ONTO SNOW
3224 IF (PSACWS(K).GT.0.) THEN
3225 NMULTS(K) = 35.E4*PSACWS(K)*FMULT*1000.
3226 QMULTS(K) = NMULTS(K)*MMULT
3228 ! CONSTRAIN SO THAT TRANSFER OF MASS FROM SNOW TO ICE CANNOT BE MORE MASS
3229 ! THAN WAS RIMED ONTO SNOW
3231 QMULTS(K) = MIN(QMULTS(K),PSACWS(K))
3232 PSACWS(K) = PSACWS(K)-QMULTS(K)
3236 ! RIMING AND SPLINTERING FROM ACCRETED RAINDROPS
3238 IF (PRACS(K).GT.0.) THEN
3239 NMULTR(K) = 35.E4*PRACS(K)*FMULT*1000.
3240 QMULTR(K) = NMULTR(K)*MMULT
3242 ! CONSTRAIN SO THAT TRANSFER OF MASS FROM SNOW TO ICE CANNOT BE MORE MASS
3243 ! THAN WAS RIMED ONTO SNOW
3245 QMULTR(K) = MIN(QMULTR(K),PRACS(K))
3247 PRACS(K) = PRACS(K)-QMULTR(K)
3256 !.......................................................................
3257 ! RIME-SPLINTERING - GRAUPEL
3258 ! HALLET-MOSSOP (1974)
3259 ! NUMBER OF SPLINTERS FORMED IS BASED ON MASS OF RIMED WATER
3261 ! DUM1 = MASS OF INDIVIDUAL SPLINTERS
3263 ! HM ADD THRESHOLD SNOW MIXING RATIO FOR RIME-SPLINTERING
3264 ! TO LIMIT RIME-SPLINTERING IN STRATIFORM CLOUDS
3266 ! IF (IHAIL.EQ.0) THEN
3268 IF (QG3D(K).GE.0.1E-3) THEN
3269 IF (QC3D(K).GE.0.5E-3.OR.QR3D(K).GE.0.1E-3) THEN
3270 IF (PSACWG(K).GT.0..OR.PRACG(K).GT.0.) THEN
3271 IF (T3D(K).LT.270.16 .AND. T3D(K).GT.265.16) THEN
3273 IF (T3D(K).GT.270.16) THEN
3275 ELSE IF (T3D(K).LE.270.16.AND.T3D(K).GT.268.16) THEN
3276 FMULT = (270.16-T3D(K))/2.
3277 ELSE IF (T3D(K).GE.265.16.AND.T3D(K).LE.268.16) THEN
3278 FMULT = (T3D(K)-265.16)/3.
3279 ELSE IF (T3D(K).LT.265.16) THEN
3283 ! 1000 IS TO CONVERT FROM KG TO G
3285 ! SPLINTERING FROM DROPLETS ACCRETED ONTO GRAUPEL
3287 IF (PSACWG(K).GT.0.) THEN
3288 NMULTG(K) = 35.E4*PSACWG(K)*FMULT*1000.
3289 QMULTG(K) = NMULTG(K)*MMULT
3291 ! CONSTRAIN SO THAT TRANSFER OF MASS FROM GRAUPEL TO ICE CANNOT BE MORE MASS
3292 ! THAN WAS RIMED ONTO GRAUPEL
3294 QMULTG(K) = MIN(QMULTG(K),PSACWG(K))
3295 PSACWG(K) = PSACWG(K)-QMULTG(K)
3299 ! RIMING AND SPLINTERING FROM ACCRETED RAINDROPS
3301 IF (PRACG(K).GT.0.) THEN
3302 NMULTRG(K) = 35.E4*PRACG(K)*FMULT*1000.
3303 QMULTRG(K) = NMULTRG(K)*MMULT
3305 ! CONSTRAIN SO THAT TRANSFER OF MASS FROM GRAUPEL TO ICE CANNOT BE MORE MASS
3306 ! THAN WAS RIMED ONTO GRAUPEL
3308 QMULTRG(K) = MIN(QMULTRG(K),PRACG(K))
3309 PRACG(K) = PRACG(K)-QMULTRG(K)
3318 !........................................................................
3319 ! CONVERSION OF RIMED CLOUD WATER ONTO SNOW TO GRAUPEL/HAIL
3321 ! IF (IHAIL.EQ.0) THEN
3322 IF (PSACWS(K).GT.0.) THEN
3323 ! ONLY ALLOW CONVERSION IF QNI > 0.1 AND QC > 0.5 G/KG FOLLOWING RUTLEDGE AND HOBBS (1984)
3324 IF (QNI3D(K).GE.0.1E-3.AND.QC3D(K).GE.0.5E-3) THEN
3326 ! PORTION OF RIMING CONVERTED TO GRAUPEL (REISNER ET AL. 1998, ORIGINALLY IS1991)
3327 PGSACW(K) = MIN(PSACWS(K),CONS17*DT*N0S(K)*QC3D(K)*QC3D(K)* &
3329 (RHO(K)*LAMS(K)**(2.*BS+2.)))
3331 ! MIX RAT CONVERTED INTO GRAUPEL AS EMBRYO (REISNER ET AL. 1998, ORIG M1990)
3332 DUM = MAX(RHOSN/(RHOG-RHOSN)*PGSACW(K),0.)
3334 ! NUMBER CONCENTRAITON OF EMBRYO GRAUPEL FROM RIMING OF SNOW
3335 NSCNG(K) = DUM/MG0*RHO(K)
3336 ! LIMIT MAX NUMBER CONVERTED TO SNOW NUMBER
3337 NSCNG(K) = MIN(NSCNG(K),NS3D(K)/DT)
3339 ! PORTION OF RIMING LEFT FOR SNOW
3340 PSACWS(K) = PSACWS(K) - PGSACW(K)
3344 ! CONVERSION OF RIMED RAINWATER ONTO SNOW CONVERTED TO GRAUPEL
3346 IF (PRACS(K).GT.0.) THEN
3347 ! ONLY ALLOW CONVERSION IF QNI > 0.1 AND QR > 0.1 G/KG FOLLOWING RUTLEDGE AND HOBBS (1984)
3348 IF (QNI3D(K).GE.0.1E-3.AND.QR3D(K).GE.0.1E-3) THEN
3349 ! PORTION OF COLLECTED RAINWATER CONVERTED TO GRAUPEL (REISNER ET AL. 1998)
3350 DUM = CONS18*(4./LAMS(K))**3*(4./LAMS(K))**3 &
3351 /(CONS18*(4./LAMS(K))**3*(4./LAMS(K))**3+ &
3352 CONS19*(4./LAMR(K))**3*(4./LAMR(K))**3)
3355 PGRACS(K) = (1.-DUM)*PRACS(K)
3356 NGRACS(K) = (1.-DUM)*NPRACS(K)
3357 ! LIMIT MAX NUMBER CONVERTED TO MIN OF EITHER RAIN OR SNOW NUMBER CONCENTRATION
3358 NGRACS(K) = MIN(NGRACS(K),NR3D(K)/DT)
3359 NGRACS(K) = MIN(NGRACS(K),NS3D(K)/DT)
3361 ! AMOUNT LEFT FOR SNOW PRODUCTION
3362 PRACS(K) = PRACS(K) - PGRACS(K)
3363 NPRACS(K) = NPRACS(K) - NGRACS(K)
3364 ! CONVERSION TO GRAUPEL DUE TO COLLECTION OF SNOW BY RAIN
3365 PSACR(K)=PSACR(K)*(1.-DUM)
3370 !.......................................................................
3371 ! FREEZING OF RAIN DROPS
3372 ! FREEZING ALLOWED BELOW -4 C
3374 IF (T3D(K).LT.269.15.AND.QR3D(K).GE.QSMALL) THEN
3376 ! IMMERSION FREEZING (BIGG 1953)
3377 ! MNUCCR(K) = CONS20*NR3D(K)*EXP(AIMM*(273.15-T3D(K)))/LAMR(K)**3 &
3380 ! NNUCCR(K) = PI*NR3D(K)*BIMM*EXP(AIMM*(273.15-T3D(K)))/LAMR(K)**3
3382 ! hm fix 7/15/13 for consistency w/ original formula
3383 MNUCCR(K) = CONS20*NR3D(K)*(EXP(AIMM*(273.15-T3D(K)))-1.)/LAMR(K)**3 &
3386 NNUCCR(K) = PI*NR3D(K)*BIMM*(EXP(AIMM*(273.15-T3D(K)))-1.)/LAMR(K)**3
3388 ! PREVENT DIVERGENCE BETWEEN MIXING RATIO AND NUMBER CONC
3389 NNUCCR(K) = MIN(NNUCCR(K),NR3D(K)/DT)
3393 !.......................................................................
3394 ! ACCRETION OF CLOUD LIQUID WATER BY RAIN
3395 ! CONTINUOUS COLLECTION EQUATION WITH
3396 ! GRAVITATIONAL COLLECTION KERNEL, DROPLET FALL SPEED NEGLECTED
3398 IF (QR3D(K).GE.1.E-8 .AND. QC3D(K).GE.1.E-8) THEN
3400 ! 12/13/06 HM ADD, REPLACE WITH NEWER FORMULA FROM
3401 ! KHAIROUTDINOV AND KOGAN 2000, MWR
3403 DUM=(QC3D(K)*QR3D(K))
3404 PRA(K) = 67.*(DUM)**1.15
3405 NPRA(K) = PRA(K)/(QC3D(K)/NC3D(K))
3408 !.......................................................................
3409 ! SELF-COLLECTION OF RAIN DROPS
3411 ! FROM NUMERICAL SIMULATION OF THE STOCHASTIC COLLECTION EQUATION
3412 ! AS DESCRINED ABOVE FOR AUTOCONVERSION
3414 IF (QR3D(K).GE.1.E-8) THEN
3415 ! include breakup add 10/09/09
3417 if (1./lamr(k).lt.dum1) then
3419 else if (1./lamr(k).ge.dum1) then
3420 dum=2.-exp(2300.*(1./lamr(k)-dum1))
3422 ! NRAGG(K) = -8.*NR3D(K)*QR3D(K)*RHO(K)
3423 NRAGG(K) = -5.78*dum*NR3D(K)*QR3D(K)*RHO(K)
3426 !.......................................................................
3427 ! AUTOCONVERSION OF CLOUD ICE TO SNOW
3428 ! FOLLOWING HARRINGTON ET AL. (1995) WITH MODIFICATION
3429 ! HERE IT IS ASSUMED THAT AUTOCONVERSION CAN ONLY OCCUR WHEN THE
3430 ! ICE IS GROWING, I.E. IN CONDITIONS OF ICE SUPERSATURATION
3432 IF (QI3D(K).GE.1.E-8 .AND.QVQVSI(K).GE.1.) THEN
3434 ! COFFI = 2./LAMI(K)
3435 ! IF (COFFI.GE.DCS) THEN
3436 NPRCI(K) = CONS21*(QV3D(K)-QVI(K))*RHO(K) &
3437 *N0I(K)*EXP(-LAMI(K)*DCS)*DV(K)/ABI(K)
3438 PRCI(K) = CONS22*NPRCI(K)
3439 NPRCI(K) = MIN(NPRCI(K),NI3D(K)/DT)
3444 !.......................................................................
3445 ! ACCRETION OF CLOUD ICE BY SNOW
3446 ! FOR THIS CALCULATION, IT IS ASSUMED THAT THE VS >> VI
3447 ! AND DS >> DI FOR CONTINUOUS COLLECTION
3449 IF (QNI3D(K).GE.1.E-8 .AND. QI3D(K).GE.QSMALL) THEN
3450 PRAI(K) = CONS23*ASN(K)*QI3D(K)*RHO(K)*N0S(K)/ &
3452 NPRAI(K) = CONS23*ASN(K)*NI3D(K)* &
3455 NPRAI(K)=MIN(NPRAI(K),NI3D(K)/DT)
3458 !.......................................................................
3459 ! HM, ADD 12/13/06, COLLISION OF RAIN AND ICE TO PRODUCE SNOW OR GRAUPEL
3460 ! FOLLOWS REISNER ET AL. 1998
3461 ! ASSUMED FALLSPEED AND SIZE OF ICE CRYSTAL << THAN FOR RAIN
3463 IF (QR3D(K).GE.1.E-8.AND.QI3D(K).GE.1.E-8.AND.T3D(K).LE.273.15) THEN
3465 ! ALLOW GRAUPEL FORMATION FROM RAIN-ICE COLLISIONS ONLY IF RAIN MIXING RATIO > 0.1 G/KG,
3466 ! OTHERWISE ADD TO SNOW
3468 IF (QR3D(K).GE.0.1E-3) THEN
3469 NIACR(K)=CONS24*NI3D(K)*N0RR(K)*ARN(K) &
3470 /LAMR(K)**(BR+3.)*RHO(K)
3471 PIACR(K)=CONS25*NI3D(K)*N0RR(K)*ARN(K) &
3472 /LAMR(K)**(BR+3.)/LAMR(K)**3*RHO(K)
3473 PRACI(K)=CONS24*QI3D(K)*N0RR(K)*ARN(K)/ &
3474 LAMR(K)**(BR+3.)*RHO(K)
3475 NIACR(K)=MIN(NIACR(K),NR3D(K)/DT)
3476 NIACR(K)=MIN(NIACR(K),NI3D(K)/DT)
3478 NIACRS(K)=CONS24*NI3D(K)*N0RR(K)*ARN(K) &
3479 /LAMR(K)**(BR+3.)*RHO(K)
3480 PIACRS(K)=CONS25*NI3D(K)*N0RR(K)*ARN(K) &
3481 /LAMR(K)**(BR+3.)/LAMR(K)**3*RHO(K)
3482 PRACIS(K)=CONS24*QI3D(K)*N0RR(K)*ARN(K)/ &
3483 LAMR(K)**(BR+3.)*RHO(K)
3484 NIACRS(K)=MIN(NIACRS(K),NR3D(K)/DT)
3485 NIACRS(K)=MIN(NIACRS(K),NI3D(K)/DT)
3489 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
3490 ! NUCLEATION OF CLOUD ICE FROM HOMOGENEOUS AND HETEROGENEOUS FREEZING ON AEROSOL
3500 ! add threshold according to Greg Thomspon
3502 if ((QVQVS(K).GE.0.999.and.T3D(K).le.265.15).or. &
3503 QVQVSI(K).ge.1.08) then
3505 ! hm, modify dec. 5, 2006, replace with cooper curve
3506 kc2 = 0.005*exp(0.304*(273.15-T3D(K)))*1000. ! convert from L-1 to m-3
3508 kc2 = min(kc2,500.e3)
3509 kc2=MAX(kc2/rho(k),0.) ! convert to kg-1
3511 IF (KC2.GT.NI3D(K)+NS3D(K)+NG3D(K)) THEN
3512 NNUCCD(K) = (KC2-NI3D(K)-NS3D(K)-NG3D(K))/DT
3513 MNUCCD(K) = NNUCCD(K)*MI0
3518 ELSE IF (INUC.EQ.1) THEN
3520 IF (T3D(K).LT.273.15.AND.QVQVSI(K).GT.1.) THEN
3522 KC2 = 0.16*1000./RHO(K) ! CONVERT FROM L-1 TO KG-1
3523 IF (KC2.GT.NI3D(K)+NS3D(K)+NG3D(K)) THEN
3524 NNUCCD(K) = (KC2-NI3D(K)-NS3D(K)-NG3D(K))/DT
3525 MNUCCD(K) = NNUCCD(K)*MI0
3529 ELSE IF (INUC.EQ.2) THEN
3531 ! IF (T3D(K).LT.273.15.AND.QVQVSI(K).GT.1.) THEN
3532 if ((QVQVS(K).GE.0.999.and.T3D(K).le.265.15).or. &
3533 QVQVSI(K).ge.1.08) then
3535 DUM = W3D(K)+WVAR(K)
3536 call mdm_prescribed_nucleati(DUM,T3D(K),QVQVS(K),QVQVSI(K),QC3D(K),RHO(K), & !TWG Fix Ice Nucleation add QVQVSI
3541 IF (KC2.GT.NI3D(K)+NS3D(K)+NG3D(K)) THEN
3542 NNUCCD(K) = (KC2-NI3D(K)-NS3D(K)-NG3D(K))/DT
3543 NNUCCD(K) = max(NNUCCD(K),0.0)
3544 MNUCCD(K) = NNUCCD(K)*MI0
3550 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
3554 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
3555 ! CALCULATE EVAP/SUB/DEP TERMS FOR QI,QNI,QR
3557 ! NO VENTILATION FOR CLOUD ICE
3559 IF (QI3D(K).GE.QSMALL) THEN
3561 EPSI = 2.*PI*N0I(K)*RHO(K)*DV(K)/(LAMI(K)*LAMI(K))
3567 IF (QNI3D(K).GE.QSMALL) THEN
3568 EPSS = 2.*PI*N0S(K)*RHO(K)*DV(K)* &
3569 (F1S/(LAMS(K)*LAMS(K))+ &
3570 F2S*(ASN(K)*RHO(K)/MU(K))**0.5* &
3571 SC(K)**(1./3.)*CONS10/ &
3577 IF (QG3D(K).GE.QSMALL) THEN
3578 EPSG = 2.*PI*N0G(K)*RHO(K)*DV(K)* &
3579 (F1S/(LAMG(K)*LAMG(K))+ &
3580 F2S*(AGN(K)*RHO(K)/MU(K))**0.5* &
3581 SC(K)**(1./3.)*CONS11/ &
3589 IF (QR3D(K).GE.QSMALL) THEN
3590 EPSR = 2.*PI*N0RR(K)*RHO(K)*DV(K)* &
3591 (F1R/(LAMR(K)*LAMR(K))+ &
3592 F2R*(ARN(K)*RHO(K)/MU(K))**0.5* &
3593 SC(K)**(1./3.)*CONS9/ &
3599 ! ONLY INCLUDE REGION OF ICE SIZE DIST < DCS
3600 ! DUM IS FRACTION OF D*N(D) < DCS
3602 ! LOGIC BELOW FOLLOWS THAT OF HARRINGTON ET AL. 1995 (JAS)
3603 IF (QI3D(K).GE.QSMALL) THEN
3604 DUM=(1.-EXP(-LAMI(K)*DCS)*(1.+LAMI(K)*DCS))
3605 PRD(K) = EPSI*(QV3D(K)-QVI(K))/ABI(K)*DUM
3609 ! ADD DEPOSITION IN TAIL OF ICE SIZE DIST TO SNOW IF SNOW IS PRESENT
3610 IF (QNI3D(K).GE.QSMALL) THEN
3611 PRDS(K) = EPSS*(QV3D(K)-QVI(K))/ABI(K)+ &
3612 EPSI*(QV3D(K)-QVI(K))/ABI(K)*(1.-DUM)
3613 ! OTHERWISE ADD TO CLOUD ICE
3615 PRD(K) = PRD(K)+EPSI*(QV3D(K)-QVI(K))/ABI(K)*(1.-DUM)
3617 ! VAPOR DPEOSITION ON GRAUPEL
3618 PRDG(K) = EPSG*(QV3D(K)-QVI(K))/ABI(K)
3620 ! NO CONDENSATION ONTO RAIN, ONLY EVAP
3622 IF (QV3D(K).LT.QVS(K)) THEN
3623 PRE(K) = EPSR*(QV3D(K)-QVS(K))/AB(K)
3624 PRE(K) = MIN(PRE(K),0.)
3629 ! MAKE SURE NOT PUSHED INTO ICE SUPERSAT/SUBSAT
3630 ! FORMULA FROM REISNER 2 SCHEME
3632 DUM = (QV3D(K)-QVI(K))/DT
3635 SUM_DEP = PRD(K)+PRDS(K)+MNUCCD(K)+PRDG(K)
3637 IF( (DUM.GT.0. .AND. SUM_DEP.GT.DUM*FUDGEF) .OR. &
3638 (DUM.LT.0. .AND. SUM_DEP.LT.DUM*FUDGEF) ) THEN
3639 MNUCCD(K) = FUDGEF*MNUCCD(K)*DUM/SUM_DEP
3640 PRD(K) = FUDGEF*PRD(K)*DUM/SUM_DEP
3641 PRDS(K) = FUDGEF*PRDS(K)*DUM/SUM_DEP
3642 PRDG(K) = FUDGEF*PRDG(K)*DUM/SUM_DEP
3645 ! IF CLOUD ICE/SNOW/GRAUPEL VAP DEPOSITION IS NEG, THEN ASSIGN TO SUBLIMATION PROCESSES
3647 IF (PRD(K).LT.0.) THEN
3651 IF (PRDS(K).LT.0.) THEN
3655 IF (PRDG(K).LT.0.) THEN
3660 !.......................................................................
3661 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
3663 ! CONSERVATION OF WATER
3664 ! THIS IS ADOPTED LOOSELY FROM MM5 RESINER CODE. HOWEVER, HERE WE
3665 ! ONLY ADJUST PROCESSES THAT ARE NEGATIVE, RATHER THAN ALL PROCESSES.
3667 ! IF MIXING RATIOS LESS THAN QSMALL, THEN NO DEPLETION OF WATER
3668 ! THROUGH MICROPHYSICAL PROCESSES, SKIP CONSERVATION
3670 ! NOTE: CONSERVATION CHECK NOT APPLIED TO NUMBER CONCENTRATION SPECIES. ADDITIONAL CATCH
3671 ! BELOW WILL PREVENT NEGATIVE NUMBER CONCENTRATION
3672 ! FOR EACH MICROPHYSICAL PROCESS WHICH PROVIDES A SOURCE FOR NUMBER, THERE IS A CHECK
3673 ! TO MAKE SURE THAT CAN'T EXCEED TOTAL NUMBER OF DEPLETED SPECIES WITH THE TIME
3676 !****SENSITIVITY - NO ICE
3687 ! ****SENSITIVITY - NO GRAUPEL
3688 IF (IGRAUP.EQ.1) THEN
3704 PIACRS(K)=PIACRS(K)+PIACR(K)
3707 PRACIS(K)=PRACIS(K)+PRACI(K)
3709 PSACWS(K)=PSACWS(K)+PGSACW(K)
3711 PRACS(K)=PRACS(K)+PGRACS(K)
3715 ! CONSERVATION OF QC
3717 DUM = (PRC(K)+PRA(K)+MNUCCC(K)+PSACWS(K)+PSACWI(K)+QMULTS(K)+PSACWG(K)+PGSACW(K)+QMULTG(K))*DT
3719 IF (DUM.GT.QC3D(K).AND.QC3D(K).GE.QSMALL) THEN
3722 PRC(K) = PRC(K)*RATIO
3723 PRA(K) = PRA(K)*RATIO
3724 MNUCCC(K) = MNUCCC(K)*RATIO
3725 PSACWS(K) = PSACWS(K)*RATIO
3726 PSACWI(K) = PSACWI(K)*RATIO
3727 QMULTS(K) = QMULTS(K)*RATIO
3728 QMULTG(K) = QMULTG(K)*RATIO
3729 PSACWG(K) = PSACWG(K)*RATIO
3730 PGSACW(K) = PGSACW(K)*RATIO
3733 ! CONSERVATION OF QI
3735 DUM = (-PRD(K)-MNUCCC(K)+PRCI(K)+PRAI(K)-QMULTS(K)-QMULTG(K)-QMULTR(K)-QMULTRG(K) &
3736 -MNUCCD(K)+PRACI(K)+PRACIS(K)-EPRD(K)-PSACWI(K))*DT
3738 IF (DUM.GT.QI3D(K).AND.QI3D(K).GE.QSMALL) THEN
3740 RATIO = (QI3D(K)/DT+PRD(K)+MNUCCC(K)+QMULTS(K)+QMULTG(K)+QMULTR(K)+QMULTRG(K)+ &
3741 MNUCCD(K)+PSACWI(K))/ &
3742 (PRCI(K)+PRAI(K)+PRACI(K)+PRACIS(K)-EPRD(K))
3744 PRCI(K) = PRCI(K)*RATIO
3745 PRAI(K) = PRAI(K)*RATIO
3746 PRACI(K) = PRACI(K)*RATIO
3747 PRACIS(K) = PRACIS(K)*RATIO
3748 EPRD(K) = EPRD(K)*RATIO
3752 ! CONSERVATION OF QR
3754 DUM=((PRACS(K)-PRE(K))+(QMULTR(K)+QMULTRG(K)-PRC(K))+(MNUCCR(K)-PRA(K))+ &
3755 PIACR(K)+PIACRS(K)+PGRACS(K)+PRACG(K))*DT
3757 IF (DUM.GT.QR3D(K).AND.QR3D(K).GE.QSMALL) THEN
3759 RATIO = (QR3D(K)/DT+PRC(K)+PRA(K))/ &
3760 (-PRE(K)+QMULTR(K)+QMULTRG(K)+PRACS(K)+MNUCCR(K)+PIACR(K)+PIACRS(K)+PGRACS(K)+PRACG(K))
3762 PRE(K) = PRE(K)*RATIO
3763 PRACS(K) = PRACS(K)*RATIO
3764 QMULTR(K) = QMULTR(K)*RATIO
3765 QMULTRG(K) = QMULTRG(K)*RATIO
3766 MNUCCR(K) = MNUCCR(K)*RATIO
3767 PIACR(K) = PIACR(K)*RATIO
3768 PIACRS(K) = PIACRS(K)*RATIO
3769 PGRACS(K) = PGRACS(K)*RATIO
3770 PRACG(K) = PRACG(K)*RATIO
3774 ! CONSERVATION OF QNI
3775 ! CONSERVATION FOR GRAUPEL SCHEME
3777 IF (IGRAUP.EQ.0) THEN
3779 DUM = (-PRDS(K)-PSACWS(K)-PRAI(K)-PRCI(K)-PRACS(K)-EPRDS(K)+PSACR(K)-PIACRS(K)-PRACIS(K))*DT
3781 IF (DUM.GT.QNI3D(K).AND.QNI3D(K).GE.QSMALL) THEN
3783 RATIO = (QNI3D(K)/DT+PRDS(K)+PSACWS(K)+PRAI(K)+PRCI(K)+PRACS(K)+PIACRS(K)+PRACIS(K))/(-EPRDS(K)+PSACR(K))
3785 EPRDS(K) = EPRDS(K)*RATIO
3786 PSACR(K) = PSACR(K)*RATIO
3790 ! FOR NO GRAUPEL, NEED TO INCLUDE FREEZING OF RAIN FOR SNOW
3791 ELSE IF (IGRAUP.EQ.1) THEN
3793 DUM = (-PRDS(K)-PSACWS(K)-PRAI(K)-PRCI(K)-PRACS(K)-EPRDS(K)+PSACR(K)-PIACRS(K)-PRACIS(K)-MNUCCR(K))*DT
3795 IF (DUM.GT.QNI3D(K).AND.QNI3D(K).GE.QSMALL) THEN
3797 RATIO = (QNI3D(K)/DT+PRDS(K)+PSACWS(K)+PRAI(K)+PRCI(K)+PRACS(K)+PIACRS(K)+PRACIS(K)+MNUCCR(K))/(-EPRDS(K)+PSACR(K))
3799 EPRDS(K) = EPRDS(K)*RATIO
3800 PSACR(K) = PSACR(K)*RATIO
3806 ! CONSERVATION OF QG
3808 DUM = (-PSACWG(K)-PRACG(K)-PGSACW(K)-PGRACS(K)-PRDG(K)-MNUCCR(K)-EPRDG(K)-PIACR(K)-PRACI(K)-PSACR(K))*DT
3810 IF (DUM.GT.QG3D(K).AND.QG3D(K).GE.QSMALL) THEN
3812 RATIO = (QG3D(K)/DT+PSACWG(K)+PRACG(K)+PGSACW(K)+PGRACS(K)+PRDG(K)+MNUCCR(K)+PSACR(K)+&
3813 PIACR(K)+PRACI(K))/(-EPRDG(K))
3815 EPRDG(K) = EPRDG(K)*RATIO
3821 QV3DTEN(K) = QV3DTEN(K)+(-PRE(K)-PRD(K)-PRDS(K)-MNUCCD(K)-EPRD(K)-EPRDS(K)-PRDG(K)-EPRDG(K))
3823 ! BUG FIX HM, 3/1/11, INCLUDE PIACR AND PIACRS
3824 T3DTEN(K) = T3DTEN(K)+(PRE(K) &
3825 *XXLV(K)+(PRD(K)+PRDS(K)+ &
3826 MNUCCD(K)+EPRD(K)+EPRDS(K)+PRDG(K)+EPRDG(K))*XXLS(K)+ &
3827 (PSACWS(K)+PSACWI(K)+MNUCCC(K)+MNUCCR(K)+ &
3828 QMULTS(K)+QMULTG(K)+QMULTR(K)+QMULTRG(K)+PRACS(K) &
3829 +PSACWG(K)+PRACG(K)+PGSACW(K)+PGRACS(K)+PIACR(K)+PIACRS(K))*XLF(K))/CPM(K)
3831 QC3DTEN(K) = QC3DTEN(K)+ &
3832 (-PRA(K)-PRC(K)-MNUCCC(K)+PCC(K)- &
3833 PSACWS(K)-PSACWI(K)-QMULTS(K)-QMULTG(K)-PSACWG(K)-PGSACW(K))
3834 QI3DTEN(K) = QI3DTEN(K)+ &
3835 (PRD(K)+EPRD(K)+PSACWI(K)+MNUCCC(K)-PRCI(K)- &
3836 PRAI(K)+QMULTS(K)+QMULTG(K)+QMULTR(K)+QMULTRG(K)+MNUCCD(K)-PRACI(K)-PRACIS(K))
3837 QR3DTEN(K) = QR3DTEN(K)+ &
3838 (PRE(K)+PRA(K)+PRC(K)-PRACS(K)-MNUCCR(K)-QMULTR(K)-QMULTRG(K) &
3839 -PIACR(K)-PIACRS(K)-PRACG(K)-PGRACS(K))
3841 IF (IGRAUP.EQ.0) THEN
3843 QNI3DTEN(K) = QNI3DTEN(K)+ &
3844 (PRAI(K)+PSACWS(K)+PRDS(K)+PRACS(K)+PRCI(K)+EPRDS(K)-PSACR(K)+PIACRS(K)+PRACIS(K))
3845 NS3DTEN(K) = NS3DTEN(K)+(NSAGG(K)+NPRCI(K)-NSCNG(K)-NGRACS(K)+NIACRS(K))
3846 QG3DTEN(K) = QG3DTEN(K)+(PRACG(K)+PSACWG(K)+PGSACW(K)+PGRACS(K)+ &
3847 PRDG(K)+EPRDG(K)+MNUCCR(K)+PIACR(K)+PRACI(K)+PSACR(K))
3848 NG3DTEN(K) = NG3DTEN(K)+(NSCNG(K)+NGRACS(K)+NNUCCR(K)+NIACR(K))
3850 ! FOR NO GRAUPEL, NEED TO INCLUDE FREEZING OF RAIN FOR SNOW
3851 ELSE IF (IGRAUP.EQ.1) THEN
3853 QNI3DTEN(K) = QNI3DTEN(K)+ &
3854 (PRAI(K)+PSACWS(K)+PRDS(K)+PRACS(K)+PRCI(K)+EPRDS(K)-PSACR(K)+PIACRS(K)+PRACIS(K)+MNUCCR(K))
3855 NS3DTEN(K) = NS3DTEN(K)+(NSAGG(K)+NPRCI(K)-NSCNG(K)-NGRACS(K)+NIACRS(K)+NNUCCR(K))
3859 NC3DTEN(K) = NC3DTEN(K)+(-NNUCCC(K)-NPSACWS(K) &
3860 -NPRA(K)-NPRC(K)-NPSACWI(K)-NPSACWG(K))
3862 NI3DTEN(K) = NI3DTEN(K)+ &
3863 (NNUCCC(K)-NPRCI(K)-NPRAI(K)+NMULTS(K)+NMULTG(K)+NMULTR(K)+NMULTRG(K)+ &
3864 NNUCCD(K)-NIACR(K)-NIACRS(K))
3866 NR3DTEN(K) = NR3DTEN(K)+(NPRC1(K)-NPRACS(K)-NNUCCR(K) &
3867 +NRAGG(K)-NIACR(K)-NIACRS(K)-NPRACG(K)-NGRACS(K))
3869 ! HM ADD, WRF-CHEM, ADD TENDENCIES FOR C2PREC
3871 C2PREC(K) = PRA(K)+PRC(K)+PSACWS(K)+QMULTS(K)+QMULTG(K)+PSACWG(K)+ &
3872 PGSACW(K)+MNUCCC(K)+PSACWI(K)
3873 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
3874 ! NOW CALCULATE SATURATION ADJUSTMENT TO CONDENSE EXTRA VAPOR ABOVE
3877 DUMT = T3D(K)+DT*T3DTEN(K)
3878 DUMQV = QV3D(K)+DT*QV3DTEN(K)
3879 ! hm, add fix for low pressure, 5/12/10
3880 dum=min(0.99*pres(k),POLYSVP(DUMT,0))
3881 DUMQSS = EP_2*dum/(PRES(K)-dum)
3882 DUMQC = QC3D(K)+DT*QC3DTEN(K)
3883 DUMQC = MAX(DUMQC,0.)
3885 ! SATURATION ADJUSTMENT FOR LIQUID
3888 PCC(K) = DUMS/(1.+XXLV(K)**2*DUMQSS/(CPM(K)*RV*DUMT**2))/DT
3889 IF (PCC(K)*DT+DUMQC.LT.0.) THEN
3893 QV3DTEN(K) = QV3DTEN(K)-PCC(K)
3894 T3DTEN(K) = T3DTEN(K)+PCC(K)*XXLV(K)/CPM(K)
3895 QC3DTEN(K) = QC3DTEN(K)+PCC(K)
3897 !.......................................................................
3898 ! ACTIVATION OF CLOUD DROPLETS
3899 ! ACTIVATION OF DROPLET CURRENTLY NOT CALCULATED
3900 ! DROPLET CONCENTRATION IS SPECIFIED !!!!!
3902 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
3903 !!TWG/amy added code to predict droplet concentration back in
3906 IF (QC3D(K)+QC3DTEN(K)*DT.GE.QSMALL) THEN
3908 ! EFFECTIVE VERTICAL VELOCITY (M/S)
3911 ! ADD SUB-GRID VERTICAL VELOCITY
3912 DUM = W3D(K)+WVAR(K)
3914 ! ASSUME MINIMUM EFF. SUB-GRID VELOCITY 0.10 M/S
3917 ELSE IF (ISUB.EQ.1) THEN
3921 ! ONLY ACTIVATE IN REGIONS OF UPWARD MOTION
3922 IF (DUM.GE.0.001) THEN
3924 IF (IBASE.EQ.1) THEN
3926 ! ACTIVATE ONLY IF THERE IS LITTLE CLOUD WATER
3927 ! OR IF AT CLOUD BASE, OR AT LOWEST MODEL LEVEL (K=1)
3931 IF (QC3D(K)+QC3DTEN(K)*DT.LE.0.05E-3/RHO(K)) THEN
3936 ELSE IF (K.GE.2) THEN
3937 IF (QC3D(K)+QC3DTEN(K)*DT.GT.0.05E-3/RHO(K).AND. &
3938 QC3D(K-1)+QC3DTEN(K-1)*DT.LE.0.05E-3/RHO(K-1)) THEN
3943 IF (IDROP.EQ.1) THEN
3944 ! ACTIVATE AT CLOUD BASE OR REGIONS WITH VERY LITTLE LIQ WATER
3947 ! USE ROGERS AND YAU (1989) TO RELATE NUMBER ACTIVATED TO W
3948 ! BASED ON TWOMEY 1959
3950 DUM=DUM*100. ! CONVERT FROM M/S TO CM/S
3951 DUM2 = 0.88*C1**(2./(K1+2.))*(7.E-2*DUM**1.5)**(K1/(K1+2.))
3952 DUM2=DUM2*1.E6 ! CONVERT FROM CM-3 TO M-3
3953 DUM2=DUM2/RHO(K) ! CONVERT FROM M-3 TO KG-1
3954 DUM2 = (DUM2-NC3D(K))/DT
3956 NC3DTEN(K) = NC3DTEN(K)+DUM2
3958 ELSE IF (IACT.EQ.2) THEN
3959 ! DROPLET ACTIVATION FROM ABDUL-RAZZAK AND GHAN (2000)
3961 SIGVL = 0.0761-1.55E-4*(T3D(K)-273.15)
3962 AACT = 2.*MW/(RHOW*RR)*SIGVL/T3D(K)
3963 ALPHA = G*MW*XXLV(K)/(CPM(K)*RR*T3D(K)**2)-G*MA/(RR*T3D(K))
3964 GAMM = RR*T3D(K)/(EVS(K)*MW)+MW*XXLV(K)**2/(CPM(K)*PRES(K)*MA*T3D(K))
3966 GG =1./(RHOW*RR*T3D(K)/(EVS(K)*DV(K)*MW)+XXLV(K)*RHOW/(KAP(K)*T3D(K))*(XXLV(K)*MW/ &
3969 PSI = 2./3.*(ALPHA*DUM/GG)**0.5*AACT
3971 ETA1 = (ALPHA*DUM/GG)**1.5/(2.*PI*RHOW*GAMM*NANEW1)
3972 ETA2 = (ALPHA*DUM/GG)**1.5/(2.*PI*RHOW*GAMM*NANEW2)
3974 SM1 = 2./BACT**0.5*(AACT/(3.*RM1))**1.5
3975 SM2 = 2./BACT**0.5*(AACT/(3.*RM2))**1.5
3977 DUM1 =1./SM1**2*(F11*(PSI/ETA1)**1.5+F21*(SM1**2/(ETA1+3.*PSI))**0.75)
3978 DUM2 =1./SM2**2*(F12*(PSI/ETA2)**1.5+F22*(SM2**2/(ETA2+3.*PSI))**0.75)
3980 SMAX = 1./(DUM1+DUM2)**0.5
3982 UU1 = 2.*LOG(SM1/SMAX)/(4.242*LOG(SIG1))
3983 UU2 = 2.*LOG(SM2/SMAX)/(4.242*LOG(SIG2))
3984 DUM1 = NANEW1/2.*(1.-DERF1(UU1))
3985 DUM2 = NANEW2/2.*(1.-DERF1(UU2))
3987 DUM2 = (DUM1+DUM2)/RHO(K) !CONVERT TO KG-1
3989 ! MAKE SURE THIS VALUE ISN'T GREATER THAN TOTAL NUMBER OF AEROSOL
3991 DUM2 = MIN((NANEW1+NANEW2)/RHO(K),DUM2)
3993 DUM2 = (DUM2-NC3D(K))/DT
3995 NC3DTEN(K) = NC3DTEN(K)+DUM2
3997 ELSE IF (IACT.EQ.4) THEN
3998 DUM = W3D(K)+WVAR(K)
3999 call mdm_prescribed_activate(DUM,T3D(K),RHO(K), &
4000 naero, naer_cu,naer_cu, maero, &
4001 dispersion_aer,hygro_aer, density_aer, DUM2, XXLV(K))
4003 DUM2 = (DUM2-NC3D(K))/DT
4005 NC3DTEN(K) = NC3DTEN(K)+DUM2
4009 !.............................................................................
4010 ELSE IF (IDROP.EQ.0) THEN
4011 ! ACTIVATE IN CLOUD INTERIOR
4012 ! FIND EQUILIBRIUM SUPERSATURATION
4014 TAUC=1./(2.*PI*RHO(k)*DV(K)*NC3D(K)*(PGAM(K)+1.)/LAMC(K))
4015 IF (EPSR.GT.1.E-8) THEN
4020 !!amy taui,taus,taug lines added in v3
4021 IF (EPSI.GT.1.E-8) THEN
4026 IF (EPSS.GT.1.E-8) THEN
4031 IF (EPSG.GT.1.E-8) THEN
4037 ! EQUILIBRIUM SS INCLUDING BERGERON EFFECT
4038 !!amy added taui,taus,taug to these lines in v3
4040 ! DUM3=(QVS(K)*RHO(K)/(PRES(K)-EVS(K))+DQSDT/CP)*G*DUM
4041 DUM3=(-QVS(K)*RHO(K)/(PRES(K)-EVS(K))+DQSDT/CP)*G*DUM
4043 DUM3=(DUM3*TAUC*TAUR*TAUI*TAUS*TAUG- &
4044 (QVS(K)-QVI(K))*(TAUC*TAUR*TAUI*TAUG+TAUC*TAUR*TAUS*TAUG+TAUC*TAUR*TAUI*TAUS))/&
4045 (TAUC*TAUR*TAUI*TAUG+TAUC*TAUR*TAUS*TAUG+TAUC*TAUR*TAUI*TAUS+ &
4046 TAUR*TAUI*TAUS*TAUG+TAUC*TAUI*TAUS*TAUG)
4048 IF (DUM3/QVS(K).GE.1.E-6) THEN
4051 ! FIND MAXIMUM ALLOWED ACTIVATION WITH NON-EQUILIBRIUM SS
4053 DUM=DUM*100. ! CONVERT FROM M/S TO CM/S
4054 DUMACT = 0.88*C1**(2./(K1+2.))*(7.E-2*DUM**1.5)**(K1/(K1+2.))
4056 ! USE POWER LAW CCN SPECTRA
4058 ! CONVERT FROM ABSOLUTE SUPERSATURATION TO SUPERSATURATION RATIO IN %
4059 DUM3=DUM3/QVS(K)*100.
4062 ! MAKE SURE VALUE DOESN'T EXCEED THAT FOR NON-EQUILIBRIUM SS
4063 DUM2=MIN(DUM2,DUMACT)
4064 DUM2=DUM2*1.E6 ! CONVERT FROM CM-3 TO M-3
4065 DUM2=DUM2/RHO(K) ! CONVERT FROM M-3 TO KG-1
4066 DUM2 = (DUM2-NC3D(K))/DT
4068 NC3DTEN(K) = NC3DTEN(K)+DUM2
4070 ELSE IF (IACT.EQ.2) THEN
4072 ! FIND MAXIMUM ALLOWED ACTIVATION WITH NON-EQUILIBRIUM SS
4074 SIGVL = 0.0761-1.55E-4*(T3D(K)-273.15)
4075 AACT = 2.*MW/(RHOW*RR)*SIGVL/T3D(K)
4076 ALPHA = G*MW*XXLV(K)/(CPM(K)*RR*T3D(K)**2)-G*MA/(RR*T3D(K))
4077 GAMM = RR*T3D(K)/(EVS(K)*MW)+MW*XXLV(K)**2/(CPM(K)*PRES(K)*MA*T3D(K))
4079 GG = 1./(RHOW*RR*T3D(K)/(EVS(K)*DV(K)*MW)+XXLV(K)*RHOW/(KAP(K)*T3D(K))*(XXLV(K)*MW/ &
4082 PSI = 2./3.*(ALPHA*DUM/GG)**0.5*AACT
4084 ETA1 = (ALPHA*DUM/GG)**1.5/(2.*PI*RHOW*GAMM*NANEW1)
4085 ETA2 = (ALPHA*DUM/GG)**1.5/(2.*PI*RHOW*GAMM*NANEW2)
4087 SM1 = 2./BACT**0.5*(AACT/(3.*RM1))**1.5
4088 SM2 = 2./BACT**0.5*(AACT/(3.*RM2))**1.5
4090 DUM1 =1./SM1**2*(F11*(PSI/ETA1)**1.5+F21*(SM1**2/(ETA1+3.*PSI))**0.75)
4091 DUM2 =1./SM2**2*(F12*(PSI/ETA2)**1.5+F22*(SM2**2/(ETA2+3.*PSI))**0.75)
4093 SMAX = 1./(DUM1+DUM2)**0.5
4095 UU1 = 2.*LOG(SM1/SMAX)/(4.242*LOG(SIG1))
4096 UU2 = 2.*LOG(SM2/SMAX)/(4.242*LOG(SIG2))
4097 DUM1 = NANEW1/2.*(1.-DERF1(UU1))
4098 DUM2 = NANEW2/2.*(1.-DERF1(UU2))
4100 DUM2 = (DUM1+DUM2)/RHO(K) !CONVERT TO KG-1
4102 ! MAKE SURE THIS VALUE ISN'T GREATER THAN TOTAL NUMBER OF AEROSOL
4104 DUMACT = MIN((NANEW1+NANEW2)/RHO(K),DUM2)
4106 ! USE LOGNORMAL AEROSOL
4107 SIGVL = 0.0761-1.55E-4*(T3D(K)-273.15)
4108 AACT = 2.*MW/(RHOW*RR)*SIGVL/T3D(K)
4110 SM1 = 2./BACT**0.5*(AACT/(3.*RM1))**1.5
4111 SM2 = 2./BACT**0.5*(AACT/(3.*RM2))**1.5
4113 ! GET SUPERSATURATION RATIO FROM ABSOLUTE SUPERSATURATION
4116 UU1 = 2.*LOG(SM1/SMAX)/(4.242*LOG(SIG1))
4117 UU2 = 2.*LOG(SM2/SMAX)/(4.242*LOG(SIG2))
4118 DUM1 = NANEW1/2.*(1.-DERF1(UU1))
4119 DUM2 = NANEW2/2.*(1.-DERF1(UU2))
4121 DUM2 = (DUM1+DUM2)/RHO(K) !CONVERT TO KG-1
4123 ! MAKE SURE THIS VALUE ISN'T GREATER THAN TOTAL NUMBER OF AEROSOL
4125 DUM2 = MIN((NANEW1+NANEW2)/RHO(K),DUM2)
4127 ! MAKE SURE ISN'T GREATER THAN NON-EQUIL. SS
4128 DUM2=MIN(DUM2,DUMACT)
4130 DUM2 = (DUM2-NC3D(K))/DT
4132 NC3DTEN(K) = NC3DTEN(K)+DUM2
4134 ELSE IF (IACT.EQ.4) THEN
4135 DUM = W3D(K)+WVAR(K)
4136 call mdm_prescribed_activate(DUM,T3D(K),RHO(K), &
4137 naero, naer_cu,naer_cu, maero, &
4138 dispersion_aer,hygro_aer, density_aer, DUM2, XXLV(K))
4140 DUM2 = (DUM2-NC3D(K))/DT
4142 NC3DTEN(K) = NC3DTEN(K)+DUM2
4145 END IF ! DUM3/QVS > 1.E-6
4148 !.......................................................................
4149 ELSE IF (IBASE.EQ.2) THEN
4152 ! USE ROGERS AND YAU (1989) TO RELATE NUMBER ACTIVATED TO W
4153 ! BASED ON TWOMEY 1959
4155 DUM=DUM*100. ! CONVERT FROM M/S TO CM/S
4156 DUM2 = 0.88*C1**(2./(K1+2.))*(7.E-2*DUM**1.5)**(K1/(K1+2.))
4157 DUM2=DUM2*1.E6 ! CONVERT FROM CM-3 TO M-3
4158 DUM2=DUM2/RHO(K) ! CONVERT FROM M-3 TO KG-1
4159 DUM2 = (DUM2-NC3D(K))/DT
4161 NC3DTEN(K) = NC3DTEN(K)+DUM2
4163 ELSE IF (IACT.EQ.2) THEN
4165 SIGVL = 0.0761-1.55E-4*(T3D(K)-273.15)
4166 AACT = 2.*MW/(RHOW*RR)*SIGVL/T3D(K)
4167 ALPHA = G*MW*XXLV(K)/(CPM(K)*RR*T3D(K)**2)-G*MA/(RR*T3D(K))
4168 GAMM = RR*T3D(K)/(EVS(K)*MW)+MW*XXLV(K)**2/(CPM(K)*PRES(K)*MA*T3D(K))
4170 GG =1./(RHOW*RR*T3D(K)/(EVS(K)*DV(K)*MW)+XXLV(K)*RHOW/(KAP(K)*T3D(K))*(XXLV(K)*MW/ &
4173 PSI = 2./3.*(ALPHA*DUM/GG)**0.5*AACT
4175 ETA1 = (ALPHA*DUM/GG)**1.5/(2.*PI*RHOW*GAMM*NANEW1)
4176 ETA2 = (ALPHA*DUM/GG)**1.5/(2.*PI*RHOW*GAMM*NANEW2)
4178 SM1 = 2./BACT**0.5*(AACT/(3.*RM1))**1.5
4179 SM2 = 2./BACT**0.5*(AACT/(3.*RM2))**1.5
4181 DUM1 =1./SM1**2*(F11*(PSI/ETA1)**1.5+F21*(SM1**2/(ETA1+3.*PSI))**0.75)
4182 DUM2 =1./SM2**2*(F12*(PSI/ETA2)**1.5+F22*(SM2**2/(ETA2+3.*PSI))**0.75)
4184 SMAX = 1./(DUM1+DUM2)**0.5
4186 UU1 = 2.*LOG(SM1/SMAX)/(4.242*LOG(SIG1))
4187 UU2 = 2.*LOG(SM2/SMAX)/(4.242*LOG(SIG2))
4188 DUM1 = NANEW1/2.*(1.-DERF1(UU1))
4189 DUM2 = NANEW2/2.*(1.-DERF1(UU2))
4191 DUM2 = (DUM1+DUM2)/RHO(K) !CONVERT TO KG-1
4193 ! MAKE SURE THIS VALUE ISN'T GREATER THAN TOTAL NUMBER OF AEROSOL
4195 DUM2 = MIN((NANEW1+NANEW2)/RHO(K),DUM2)
4197 DUM2 = (DUM2-NC3D(K))/DT
4199 NC3DTEN(K) = NC3DTEN(K)+DUM2
4201 ELSE IF (IACT.EQ.4) THEN
4202 DUM = W3D(K)+WVAR(K)
4203 call mdm_prescribed_activate(DUM,T3D(K),RHO(K), &
4204 naero, naer_cu,naer_cu, maero, &
4205 dispersion_aer,hygro_aer, density_aer, DUM2, XXLV(K))
4207 DUM2 = (DUM2-NC3D(K))/DT
4209 NC3DTEN(K) = NC3DTEN(K)+DUM2
4214 END IF ! QC3D > QSMALL
4219 ! SUBLIMATE, MELT, OR EVAPORATE NUMBER CONCENTRATION
4220 ! THIS FORMULATION ASSUMES 1:1 RATIO BETWEEN MASS LOSS AND
4221 ! LOSS OF NUMBER CONCENTRATION
4223 ! IF (PCC(K).LT.0.) THEN
4224 ! DUM = PCC(K)*DT/QC3D(K)
4225 ! DUM = MAX(-1.,DUM)
4226 ! NSUBC(K) = DUM*NC3D(K)/DT
4229 IF (EPRD(K).LT.0.) THEN
4230 DUM = EPRD(K)*DT/QI3D(K)
4232 NSUBI(K) = DUM*NI3D(K)/DT
4234 IF (EPRDS(K).LT.0.) THEN
4235 DUM = EPRDS(K)*DT/QNI3D(K)
4237 NSUBS(K) = DUM*NS3D(K)/DT
4239 IF (PRE(K).LT.0.) THEN
4240 DUM = PRE(K)*DT/QR3D(K)
4242 NSUBR(K) = DUM*NR3D(K)/DT
4244 IF (EPRDG(K).LT.0.) THEN
4245 DUM = EPRDG(K)*DT/QG3D(K)
4247 NSUBG(K) = DUM*NG3D(K)/DT
4256 ! NC3DTEN(K) = NC3DTEN(K)+NSUBC(K)
4257 NI3DTEN(K) = NI3DTEN(K)+NSUBI(K)
4258 NS3DTEN(K) = NS3DTEN(K)+NSUBS(K)
4259 NG3DTEN(K) = NG3DTEN(K)+NSUBG(K)
4260 NR3DTEN(K) = NR3DTEN(K)+NSUBR(K)
4263 IF ( has_wetscav ) THEN
4264 evapprod(k) = - PRE(K) - EPRDS(K) - EPRDG(K)
4265 rainprod(k) = PRA(K) + PRC(K) + PSACWS(K) + PSACWG(K) + PGSACW(K) &
4266 + PRAI(K) + PRCI(K) + PRACI(K) + PRACIS(K) + &
4271 END IF !!!!!! TEMPERATURE
4273 ! SWITCH LTRUE TO 1, SINCE HYDROMETEORS ARE PRESENT
4280 ! INITIALIZE PRECIP AND SNOW RATES
4287 ! IF THERE ARE NO HYDROMETEORS, THEN SKIP TO END OF SUBROUTINE
4289 IF (LTRUE.EQ.0) GOTO 400
4291 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
4292 !.......................................................................
4293 ! CALCULATE SEDIMENATION
4294 ! THE NUMERICS HERE FOLLOW FROM REISNER ET AL. (1998)
4295 ! FALLOUT TERMS ARE CALCULATED ON SPLIT TIME STEPS TO ENSURE NUMERICAL
4296 ! STABILITY, I.E. COURANT# < 1
4298 !.......................................................................
4304 DUMI(K) = QI3D(K)+QI3DTEN(K)*DT
4305 DUMQS(K) = QNI3D(K)+QNI3DTEN(K)*DT
4306 DUMR(K) = QR3D(K)+QR3DTEN(K)*DT
4307 DUMFNI(K) = NI3D(K)+NI3DTEN(K)*DT
4308 DUMFNS(K) = NS3D(K)+NS3DTEN(K)*DT
4309 DUMFNR(K) = NR3D(K)+NR3DTEN(K)*DT
4310 DUMC(K) = QC3D(K)+QC3DTEN(K)*DT
4311 DUMFNC(K) = NC3D(K)+NC3DTEN(K)*DT
4312 DUMG(K) = QG3D(K)+QG3DTEN(K)*DT
4313 DUMFNG(K) = NG3D(K)+NG3DTEN(K)*DT
4315 ! SWITCH FOR CONSTANT DROPLET NUMBER
4316 IF (iinum.EQ.1) THEN
4320 ! GET DUMMY LAMDA FOR SEDIMENTATION CALCULATIONS
4322 ! MAKE SURE NUMBER CONCENTRATIONS ARE POSITIVE
4323 DUMFNI(K) = MAX(0.,DUMFNI(K))
4324 DUMFNS(K) = MAX(0.,DUMFNS(K))
4325 DUMFNC(K) = MAX(0.,DUMFNC(K))
4326 DUMFNR(K) = MAX(0.,DUMFNR(K))
4327 DUMFNG(K) = MAX(0.,DUMFNG(K))
4329 !......................................................................
4332 IF (DUMI(K).GE.QSMALL) THEN
4333 DLAMI = (CONS12*DUMFNI(K)/DUMI(K))**(1./DI)
4334 DLAMI=MAX(DLAMI,LAMMINI)
4335 DLAMI=MIN(DLAMI,LAMMAXI)
4337 !......................................................................
4340 IF (DUMR(K).GE.QSMALL) THEN
4341 DLAMR = (PI*RHOW*DUMFNR(K)/DUMR(K))**(1./3.)
4342 DLAMR=MAX(DLAMR,LAMMINR)
4343 DLAMR=MIN(DLAMR,LAMMAXR)
4345 !......................................................................
4348 IF (DUMC(K).GE.QSMALL) THEN
4349 DUM = PRES(K)/(287.15*T3D(K))
4350 PGAM(K)=0.0005714*(NC3D(K)/1.E6*DUM)+0.2714
4351 PGAM(K)=1./(PGAM(K)**2)-1.
4352 PGAM(K)=MAX(PGAM(K),2.)
4353 PGAM(K)=MIN(PGAM(K),10.)
4355 DLAMC = (CONS26*DUMFNC(K)*GAMMA(PGAM(K)+4.)/(DUMC(K)*GAMMA(PGAM(K)+1.)))**(1./3.)
4356 LAMMIN = (PGAM(K)+1.)/60.E-6
4357 LAMMAX = (PGAM(K)+1.)/1.E-6
4358 DLAMC=MAX(DLAMC,LAMMIN)
4359 DLAMC=MIN(DLAMC,LAMMAX)
4361 !......................................................................
4364 IF (DUMQS(K).GE.QSMALL) THEN
4365 DLAMS = (CONS1*DUMFNS(K)/ DUMQS(K))**(1./DS)
4366 DLAMS=MAX(DLAMS,LAMMINS)
4367 DLAMS=MIN(DLAMS,LAMMAXS)
4369 !......................................................................
4372 IF (DUMG(K).GE.QSMALL) THEN
4373 DLAMG = (CONS2*DUMFNG(K)/ DUMG(K))**(1./DG)
4374 DLAMG=MAX(DLAMG,LAMMING)
4375 DLAMG=MIN(DLAMG,LAMMAXG)
4378 !......................................................................
4379 ! CALCULATE NUMBER-WEIGHTED AND MASS-WEIGHTED TERMINAL FALL SPEEDS
4383 IF (DUMC(K).GE.QSMALL) THEN
4384 UNC = ACN(K)*GAMMA(1.+BC+PGAM(K))/ (DLAMC**BC*GAMMA(PGAM(K)+1.))
4385 UMC = ACN(K)*GAMMA(4.+BC+PGAM(K))/ (DLAMC**BC*GAMMA(PGAM(K)+4.))
4391 IF (DUMI(K).GE.QSMALL) THEN
4392 UNI = AIN(K)*CONS27/DLAMI**BI
4393 UMI = AIN(K)*CONS28/(DLAMI**BI)
4399 IF (DUMR(K).GE.QSMALL) THEN
4400 UNR = ARN(K)*CONS6/DLAMR**BR
4401 UMR = ARN(K)*CONS4/(DLAMR**BR)
4407 IF (DUMQS(K).GE.QSMALL) THEN
4408 UMS = ASN(K)*CONS3/(DLAMS**BS)
4409 UNS = ASN(K)*CONS5/DLAMS**BS
4415 IF (DUMG(K).GE.QSMALL) THEN
4416 UMG = AGN(K)*CONS7/(DLAMG**BG)
4417 UNG = AGN(K)*CONS8/DLAMG**BG
4423 ! SET REALISTIC LIMITS ON FALLSPEED
4426 dum=(rhosu/rho(k))**0.54
4427 UMS=MIN(UMS,1.2*dum)
4428 UNS=MIN(UNS,1.2*dum)
4430 ! fix for correction by AA 4/6/11
4431 UMI=MIN(UMI,1.2*(rhosu/rho(k))**0.35)
4432 UNI=MIN(UNI,1.2*(rhosu/rho(k))**0.35)
4433 UMR=MIN(UMR,9.1*dum)
4434 UNR=MIN(UNR,9.1*dum)
4435 UMG=MIN(UMG,20.*dum)
4436 UNG=MIN(UNG,20.*dum)
4449 ! V3.3 MODIFY FALLSPEED BELOW LEVEL OF PRECIP
4451 IF (K.LE.KTE-1) THEN
4452 IF (FR(K).LT.1.E-10) THEN
4455 IF (FI(K).LT.1.E-10) THEN
4458 IF (FNI(K).LT.1.E-10) THEN
4461 IF (FS(K).LT.1.E-10) THEN
4464 IF (FNS(K).LT.1.E-10) THEN
4467 IF (FNR(K).LT.1.E-10) THEN
4470 IF (FC(K).LT.1.E-10) THEN
4473 IF (FNC(K).LT.1.E-10) THEN
4476 IF (FG(K).LT.1.E-10) THEN
4479 IF (FNG(K).LT.1.E-10) THEN
4484 ! CALCULATE NUMBER OF SPLIT TIME STEPS
4486 RGVM = MAX(FR(K),FI(K),FS(K),FC(K),FNI(K),FNR(K),FNS(K),FNC(K),FG(K),FNG(K))
4487 ! VVT CHANGED IFIX -> INT (GENERIC FUNCTION)
4488 NSTEP = MAX(INT(RGVM*DT/DZQ(K)+1.),NSTEP)
4490 ! MULTIPLY VARIABLES BY RHO
4491 DUMR(k) = DUMR(k)*RHO(K)
4492 DUMI(k) = DUMI(k)*RHO(K)
4493 DUMFNI(k) = DUMFNI(K)*RHO(K)
4494 DUMQS(k) = DUMQS(K)*RHO(K)
4495 DUMFNS(k) = DUMFNS(K)*RHO(K)
4496 DUMFNR(k) = DUMFNR(K)*RHO(K)
4497 DUMC(k) = DUMC(K)*RHO(K)
4498 DUMFNC(k) = DUMFNC(K)*RHO(K)
4499 DUMG(k) = DUMG(K)*RHO(K)
4500 DUMFNG(k) = DUMFNG(K)*RHO(K)
4507 FALOUTR(K) = FR(K)*DUMR(K)
4508 FALOUTI(K) = FI(K)*DUMI(K)
4509 FALOUTNI(K) = FNI(K)*DUMFNI(K)
4510 FALOUTS(K) = FS(K)*DUMQS(K)
4511 FALOUTNS(K) = FNS(K)*DUMFNS(K)
4512 FALOUTNR(K) = FNR(K)*DUMFNR(K)
4513 FALOUTC(K) = FC(K)*DUMC(K)
4514 FALOUTNC(K) = FNC(K)*DUMFNC(K)
4515 FALOUTG(K) = FG(K)*DUMG(K)
4516 FALOUTNG(K) = FNG(K)*DUMFNG(K)
4522 FALTNDR = FALOUTR(K)/DZQ(k)
4523 FALTNDI = FALOUTI(K)/DZQ(k)
4524 FALTNDNI = FALOUTNI(K)/DZQ(k)
4525 FALTNDS = FALOUTS(K)/DZQ(k)
4526 FALTNDNS = FALOUTNS(K)/DZQ(k)
4527 FALTNDNR = FALOUTNR(K)/DZQ(k)
4528 FALTNDC = FALOUTC(K)/DZQ(k)
4529 FALTNDNC = FALOUTNC(K)/DZQ(k)
4530 FALTNDG = FALOUTG(K)/DZQ(k)
4531 FALTNDNG = FALOUTNG(K)/DZQ(k)
4532 ! ADD FALLOUT TERMS TO EULERIAN TENDENCIES
4534 QRSTEN(K) = QRSTEN(K)-FALTNDR/NSTEP/RHO(k)
4535 QISTEN(K) = QISTEN(K)-FALTNDI/NSTEP/RHO(k)
4536 NI3DTEN(K) = NI3DTEN(K)-FALTNDNI/NSTEP/RHO(k)
4537 QNISTEN(K) = QNISTEN(K)-FALTNDS/NSTEP/RHO(k)
4538 NS3DTEN(K) = NS3DTEN(K)-FALTNDNS/NSTEP/RHO(k)
4539 NR3DTEN(K) = NR3DTEN(K)-FALTNDNR/NSTEP/RHO(k)
4540 QCSTEN(K) = QCSTEN(K)-FALTNDC/NSTEP/RHO(k)
4541 NC3DTEN(K) = NC3DTEN(K)-FALTNDNC/NSTEP/RHO(k)
4542 QGSTEN(K) = QGSTEN(K)-FALTNDG/NSTEP/RHO(k)
4543 NG3DTEN(K) = NG3DTEN(K)-FALTNDNG/NSTEP/RHO(k)
4545 DUMR(K) = DUMR(K)-FALTNDR*DT/NSTEP
4546 DUMI(K) = DUMI(K)-FALTNDI*DT/NSTEP
4547 DUMFNI(K) = DUMFNI(K)-FALTNDNI*DT/NSTEP
4548 DUMQS(K) = DUMQS(K)-FALTNDS*DT/NSTEP
4549 DUMFNS(K) = DUMFNS(K)-FALTNDNS*DT/NSTEP
4550 DUMFNR(K) = DUMFNR(K)-FALTNDNR*DT/NSTEP
4551 DUMC(K) = DUMC(K)-FALTNDC*DT/NSTEP
4552 DUMFNC(K) = DUMFNC(K)-FALTNDNC*DT/NSTEP
4553 DUMG(K) = DUMG(K)-FALTNDG*DT/NSTEP
4554 DUMFNG(K) = DUMFNG(K)-FALTNDNG*DT/NSTEP
4557 FALTNDR = (FALOUTR(K+1)-FALOUTR(K))/DZQ(K)
4558 FALTNDI = (FALOUTI(K+1)-FALOUTI(K))/DZQ(K)
4559 FALTNDNI = (FALOUTNI(K+1)-FALOUTNI(K))/DZQ(K)
4560 FALTNDS = (FALOUTS(K+1)-FALOUTS(K))/DZQ(K)
4561 FALTNDNS = (FALOUTNS(K+1)-FALOUTNS(K))/DZQ(K)
4562 FALTNDNR = (FALOUTNR(K+1)-FALOUTNR(K))/DZQ(K)
4563 FALTNDC = (FALOUTC(K+1)-FALOUTC(K))/DZQ(K)
4564 FALTNDNC = (FALOUTNC(K+1)-FALOUTNC(K))/DZQ(K)
4565 FALTNDG = (FALOUTG(K+1)-FALOUTG(K))/DZQ(K)
4566 FALTNDNG = (FALOUTNG(K+1)-FALOUTNG(K))/DZQ(K)
4568 ! ADD FALLOUT TERMS TO EULERIAN TENDENCIES
4570 QRSTEN(K) = QRSTEN(K)+FALTNDR/NSTEP/RHO(k)
4571 QISTEN(K) = QISTEN(K)+FALTNDI/NSTEP/RHO(k)
4572 NI3DTEN(K) = NI3DTEN(K)+FALTNDNI/NSTEP/RHO(k)
4573 QNISTEN(K) = QNISTEN(K)+FALTNDS/NSTEP/RHO(k)
4574 NS3DTEN(K) = NS3DTEN(K)+FALTNDNS/NSTEP/RHO(k)
4575 NR3DTEN(K) = NR3DTEN(K)+FALTNDNR/NSTEP/RHO(k)
4576 QCSTEN(K) = QCSTEN(K)+FALTNDC/NSTEP/RHO(k)
4577 NC3DTEN(K) = NC3DTEN(K)+FALTNDNC/NSTEP/RHO(k)
4578 QGSTEN(K) = QGSTEN(K)+FALTNDG/NSTEP/RHO(k)
4579 NG3DTEN(K) = NG3DTEN(K)+FALTNDNG/NSTEP/RHO(k)
4581 DUMR(K) = DUMR(K)+FALTNDR*DT/NSTEP
4582 DUMI(K) = DUMI(K)+FALTNDI*DT/NSTEP
4583 DUMFNI(K) = DUMFNI(K)+FALTNDNI*DT/NSTEP
4584 DUMQS(K) = DUMQS(K)+FALTNDS*DT/NSTEP
4585 DUMFNS(K) = DUMFNS(K)+FALTNDNS*DT/NSTEP
4586 DUMFNR(K) = DUMFNR(K)+FALTNDNR*DT/NSTEP
4587 DUMC(K) = DUMC(K)+FALTNDC*DT/NSTEP
4588 DUMFNC(K) = DUMFNC(K)+FALTNDNC*DT/NSTEP
4589 DUMG(K) = DUMG(K)+FALTNDG*DT/NSTEP
4590 DUMFNG(K) = DUMFNG(K)+FALTNDNG*DT/NSTEP
4592 ! FOR WRF-CHEM, NEED PRECIP RATES (UNITS OF KG/M^2/S)
4593 CSED(K)=CSED(K)+FALOUTC(K)/NSTEP
4594 ISED(K)=ISED(K)+FALOUTI(K)/NSTEP
4595 SSED(K)=SSED(K)+FALOUTS(K)/NSTEP
4596 GSED(K)=GSED(K)+FALOUTG(K)/NSTEP
4597 RSED(K)=RSED(K)+FALOUTR(K)/NSTEP
4600 ! GET PRECIPITATION AND SNOWFALL ACCUMULATION DURING THE TIME STEP
4601 ! FACTOR OF 1000 CONVERTS FROM M TO MM, BUT DIVISION BY DENSITY
4602 ! OF LIQUID WATER CANCELS THIS FACTOR OF 1000
4604 PRECRT = PRECRT+(FALOUTR(KTS)+FALOUTC(KTS)+FALOUTS(KTS)+FALOUTI(KTS)+FALOUTG(KTS)) &
4606 SNOWRT = SNOWRT+(FALOUTS(KTS)+FALOUTI(KTS)+FALOUTG(KTS))*DT/NSTEP
4608 SNOWPRT = SNOWPRT+(FALOUTI(KTS)+FALOUTS(KTS))*DT/NSTEP
4609 GRPLPRT = GRPLPRT+(FALOUTG(KTS))*DT/NSTEP
4615 ! ADD ON SEDIMENTATION TENDENCIES FOR MIXING RATIO TO REST OF TENDENCIES
4617 QR3DTEN(K)=QR3DTEN(K)+QRSTEN(K)
4618 QI3DTEN(K)=QI3DTEN(K)+QISTEN(K)
4619 QC3DTEN(K)=QC3DTEN(K)+QCSTEN(K)
4620 QG3DTEN(K)=QG3DTEN(K)+QGSTEN(K)
4621 QNI3DTEN(K)=QNI3DTEN(K)+QNISTEN(K)
4623 ! PUT ALL CLOUD ICE IN SNOW CATEGORY IF MEAN DIAMETER EXCEEDS 2 * dcs
4626 ! IF (QI3D(K).GE.QSMALL.AND.T3D(K).LT.273.15) THEN
4627 IF (QI3D(K).GE.QSMALL.AND.T3D(K).LT.273.15.AND.LAMI(K).GE.1.E-10) THEN
4628 IF (1./LAMI(K).GE.2.*DCS) THEN
4629 QNI3DTEN(K) = QNI3DTEN(K)+QI3D(K)/DT+ QI3DTEN(K)
4630 NS3DTEN(K) = NS3DTEN(K)+NI3D(K)/DT+ NI3DTEN(K)
4631 QI3DTEN(K) = -QI3D(K)/DT
4632 NI3DTEN(K) = -NI3D(K)/DT
4636 ! hm add tendencies here, then call sizeparameter
4637 ! to ensure consisitency between mixing ratio and number concentration
4639 QC3D(k) = QC3D(k)+QC3DTEN(k)*DT
4640 QI3D(k) = QI3D(k)+QI3DTEN(k)*DT
4641 QNI3D(k) = QNI3D(k)+QNI3DTEN(k)*DT
4642 QR3D(k) = QR3D(k)+QR3DTEN(k)*DT
4643 NC3D(k) = NC3D(k)+NC3DTEN(k)*DT
4644 NI3D(k) = NI3D(k)+NI3DTEN(k)*DT
4645 NS3D(k) = NS3D(k)+NS3DTEN(k)*DT
4646 NR3D(k) = NR3D(k)+NR3DTEN(k)*DT
4648 IF (IGRAUP.EQ.0) THEN
4649 QG3D(k) = QG3D(k)+QG3DTEN(k)*DT
4650 NG3D(k) = NG3D(k)+NG3DTEN(k)*DT
4653 ! ADD TEMPERATURE AND WATER VAPOR TENDENCIES FROM MICROPHYSICS
4654 T3D(K) = T3D(K)+T3DTEN(k)*DT
4655 QV3D(K) = QV3D(K)+QV3DTEN(k)*DT
4657 ! SATURATION VAPOR PRESSURE AND MIXING RATIO
4659 ! hm, add fix for low pressure, 5/12/10
4660 EVS(K) = min(0.99*pres(k),POLYSVP(T3D(K),0)) ! PA
4661 EIS(K) = min(0.99*pres(k),POLYSVP(T3D(K),1)) ! PA
4663 ! MAKE SURE ICE SATURATION DOESN'T EXCEED WATER SAT. NEAR FREEZING
4665 IF (EIS(K).GT.EVS(K)) EIS(K) = EVS(K)
4667 QVS(K) = EP_2*EVS(K)/(PRES(K)-EVS(K))
4668 QVI(K) = EP_2*EIS(K)/(PRES(K)-EIS(K))
4670 QVQVS(K) = QV3D(K)/QVS(K)
4671 QVQVSI(K) = QV3D(K)/QVI(K)
4673 ! AT SUBSATURATION, REMOVE SMALL AMOUNTS OF CLOUD/PRECIP WATER
4674 ! hm 7/9/09 change limit to 1.e-8
4676 IF (QVQVS(K).LT.0.9) THEN
4677 IF (QR3D(K).LT.1.E-8) THEN
4678 QV3D(K)=QV3D(K)+QR3D(K)
4679 T3D(K)=T3D(K)-QR3D(K)*XXLV(K)/CPM(K)
4682 IF (QC3D(K).LT.1.E-8) THEN
4683 QV3D(K)=QV3D(K)+QC3D(K)
4684 T3D(K)=T3D(K)-QC3D(K)*XXLV(K)/CPM(K)
4689 IF (QVQVSI(K).LT.0.9) THEN
4690 IF (QI3D(K).LT.1.E-8) THEN
4691 QV3D(K)=QV3D(K)+QI3D(K)
4692 T3D(K)=T3D(K)-QI3D(K)*XXLS(K)/CPM(K)
4695 IF (QNI3D(K).LT.1.E-8) THEN
4696 QV3D(K)=QV3D(K)+QNI3D(K)
4697 T3D(K)=T3D(K)-QNI3D(K)*XXLS(K)/CPM(K)
4700 IF (QG3D(K).LT.1.E-8) THEN
4701 QV3D(K)=QV3D(K)+QG3D(K)
4702 T3D(K)=T3D(K)-QG3D(K)*XXLS(K)/CPM(K)
4707 !..................................................................
4708 ! IF MIXING RATIO < QSMALL SET MIXING RATIO AND NUMBER CONC TO ZERO
4710 IF (QC3D(K).LT.QSMALL) THEN
4715 IF (QR3D(K).LT.QSMALL) THEN
4720 IF (QI3D(K).LT.QSMALL) THEN
4725 IF (QNI3D(K).LT.QSMALL) THEN
4730 IF (QG3D(K).LT.QSMALL) THEN
4736 !..................................
4737 ! IF THERE IS NO CLOUD/PRECIP WATER, THEN SKIP CALCULATIONS
4739 IF (QC3D(K).LT.QSMALL.AND.QI3D(K).LT.QSMALL.AND.QNI3D(K).LT.QSMALL &
4740 .AND.QR3D(K).LT.QSMALL.AND.QG3D(K).LT.QSMALL) GOTO 500
4742 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
4743 ! CALCULATE INSTANTANEOUS PROCESSES
4745 ! ADD MELTING OF CLOUD ICE TO FORM RAIN
4747 IF (QI3D(K).GE.QSMALL.AND.T3D(K).GE.273.15) THEN
4748 QR3D(K) = QR3D(K)+QI3D(K)
4749 T3D(K) = T3D(K)-QI3D(K)*XLF(K)/CPM(K)
4751 tqimelt(K)=QI3D(K)/DT
4754 NR3D(K) = NR3D(K)+NI3D(K)
4758 ! ****SENSITIVITY - NO ICE
4759 IF (ILIQ.EQ.1) GOTO 778
4761 ! HOMOGENEOUS FREEZING OF CLOUD WATER
4763 IF (T3D(K).LE.233.15.AND.QC3D(K).GE.QSMALL) THEN
4764 QI3D(K)=QI3D(K)+QC3D(K)
4765 T3D(K)=T3D(K)+QC3D(K)*XLF(K)/CPM(K)
4767 NI3D(K)=NI3D(K)+NC3D(K)
4771 ! HOMOGENEOUS FREEZING OF RAIN
4773 IF (IGRAUP.EQ.0) THEN
4775 IF (T3D(K).LE.233.15.AND.QR3D(K).GE.QSMALL) THEN
4776 QG3D(K) = QG3D(K)+QR3D(K)
4777 T3D(K) = T3D(K)+QR3D(K)*XLF(K)/CPM(K)
4779 NG3D(K) = NG3D(K)+ NR3D(K)
4783 ELSE IF (IGRAUP.EQ.1) THEN
4785 IF (T3D(K).LE.233.15.AND.QR3D(K).GE.QSMALL) THEN
4786 QNI3D(K) = QNI3D(K)+QR3D(K)
4787 T3D(K) = T3D(K)+QR3D(K)*XLF(K)/CPM(K)
4789 NS3D(K) = NS3D(K)+NR3D(K)
4797 ! MAKE SURE NUMBER CONCENTRATIONS AREN'T NEGATIVE
4799 NI3D(K) = MAX(0.,NI3D(K))
4800 NS3D(K) = MAX(0.,NS3D(K))
4801 NC3D(K) = MAX(0.,NC3D(K))
4802 NR3D(K) = MAX(0.,NR3D(K))
4803 NG3D(K) = MAX(0.,NG3D(K))
4805 !......................................................................
4808 IF (QI3D(K).GE.QSMALL) THEN
4809 LAMI(K) = (CONS12* &
4810 NI3D(K)/QI3D(K))**(1./DI)
4816 IF (LAMI(K).LT.LAMMINI) THEN
4820 N0I(K) = LAMI(K)**4*QI3D(K)/CONS12
4822 NI3D(K) = N0I(K)/LAMI(K)
4823 ELSE IF (LAMI(K).GT.LAMMAXI) THEN
4825 N0I(K) = LAMI(K)**4*QI3D(K)/CONS12
4827 NI3D(K) = N0I(K)/LAMI(K)
4831 !......................................................................
4834 IF (QR3D(K).GE.QSMALL) THEN
4835 LAMR(K) = (PI*RHOW*NR3D(K)/QR3D(K))**(1./3.)
4841 IF (LAMR(K).LT.LAMMINR) THEN
4845 N0RR(K) = LAMR(K)**4*QR3D(K)/(PI*RHOW)
4847 NR3D(K) = N0RR(K)/LAMR(K)
4848 ELSE IF (LAMR(K).GT.LAMMAXR) THEN
4850 N0RR(K) = LAMR(K)**4*QR3D(K)/(PI*RHOW)
4852 NR3D(K) = N0RR(K)/LAMR(K)
4857 !......................................................................
4860 ! MARTIN ET AL. (1994) FORMULA FOR PGAM
4862 IF (QC3D(K).GE.QSMALL) THEN
4864 DUM = PRES(K)/(287.15*T3D(K))
4865 PGAM(K)=0.0005714*(NC3D(K)/1.E6*DUM)+0.2714
4866 PGAM(K)=1./(PGAM(K)**2)-1.
4867 PGAM(K)=MAX(PGAM(K),2.)
4868 PGAM(K)=MIN(PGAM(K),10.)
4872 LAMC(K) = (CONS26*NC3D(K)*GAMMA(PGAM(K)+4.)/ &
4873 (QC3D(K)*GAMMA(PGAM(K)+1.)))**(1./3.)
4875 ! LAMMIN, 60 MICRON DIAMETER
4878 LAMMIN = (PGAM(K)+1.)/60.E-6
4879 LAMMAX = (PGAM(K)+1.)/1.E-6
4881 IF (LAMC(K).LT.LAMMIN) THEN
4883 NC3D(K) = EXP(3.*LOG(LAMC(K))+LOG(QC3D(K))+ &
4884 LOG(GAMMA(PGAM(K)+1.))-LOG(GAMMA(PGAM(K)+4.)))/CONS26
4886 ELSE IF (LAMC(K).GT.LAMMAX) THEN
4888 NC3D(K) = EXP(3.*LOG(LAMC(K))+LOG(QC3D(K))+ &
4889 LOG(GAMMA(PGAM(K)+1.))-LOG(GAMMA(PGAM(K)+4.)))/CONS26
4895 !......................................................................
4898 IF (QNI3D(K).GE.QSMALL) THEN
4899 LAMS(K) = (CONS1*NS3D(K)/QNI3D(K))**(1./DS)
4905 IF (LAMS(K).LT.LAMMINS) THEN
4907 N0S(K) = LAMS(K)**4*QNI3D(K)/CONS1
4909 NS3D(K) = N0S(K)/LAMS(K)
4911 ELSE IF (LAMS(K).GT.LAMMAXS) THEN
4914 N0S(K) = LAMS(K)**4*QNI3D(K)/CONS1
4915 NS3D(K) = N0S(K)/LAMS(K)
4920 !......................................................................
4923 IF (QG3D(K).GE.QSMALL) THEN
4924 LAMG(K) = (CONS2*NG3D(K)/QG3D(K))**(1./DG)
4930 IF (LAMG(K).LT.LAMMING) THEN
4932 N0G(K) = LAMG(K)**4*QG3D(K)/CONS2
4934 NG3D(K) = N0G(K)/LAMG(K)
4936 ELSE IF (LAMG(K).GT.LAMMAXG) THEN
4939 N0G(K) = LAMG(K)**4*QG3D(K)/CONS2
4941 NG3D(K) = N0G(K)/LAMG(K)
4948 ! CALCULATE EFFECTIVE RADIUS
4950 IF (QI3D(K).GE.QSMALL) THEN
4951 EFFI(K) = 3./LAMI(K)/2.*1.E6
4953 !EFFI(K) = 25. !TWG change for consistency with RRTMG
4957 IF (QNI3D(K).GE.QSMALL) THEN
4958 EFFS(K) = 3./LAMS(K)/2.*1.E6
4963 IF (QR3D(K).GE.QSMALL) THEN
4964 EFFR(K) = 3./LAMR(K)/2.*1.E6
4969 IF (QC3D(K).GE.QSMALL) THEN
4970 EFFC(K) = GAMMA(PGAM(K)+4.)/ &
4971 GAMMA(PGAM(K)+3.)/LAMC(K)/2.*1.E6
4977 IF (QG3D(K).GE.QSMALL) THEN
4978 EFFG(K) = 3./LAMG(K)/2.*1.E6
4983 ! HM ADD 1/10/06, ADD UPPER BOUND ON ICE NUMBER, THIS IS NEEDED
4984 ! TO PREVENT VERY LARGE ICE NUMBER DUE TO HOMOGENEOUS FREEZING
4985 ! OF DROPLETS, ESPECIALLY WHEN INUM = 1, SET MAX AT 10 CM-3
4986 ! NI3D(K) = MIN(NI3D(K),10.E6/RHO(K))
4987 ! HM, 12/28/12, LOWER MAXIMUM ICE CONCENTRATION TO ADDRESS PROBLEM
4988 ! OF EXCESSIVE AND PERSISTENT ANVIL
4989 ! NOTE: THIS MAY CHANGE/REDUCE SENSITIVITY TO AEROSOL/CCN CONCENTRATION
4990 NI3D(K) = MIN(NI3D(K),0.3E6/RHO(K))
4992 ! ADD BOUND ON DROPLET NUMBER - CANNOT EXCEED AEROSOL CONCENTRATION
4993 IF (iinum.EQ.0.AND.IACT.EQ.2) THEN
4994 NC3D(K) = MIN(NC3D(K),(NANEW1+NANEW2)/RHO(K))
4996 ! SWITCH FOR CONSTANT DROPLET NUMBER
4997 IF (iinum.EQ.1) THEN
4998 ! CHANGE NDCNST FROM CM-3 TO KG-1
4999 NC3D(K) = NDCNST*1.E6/RHO(K)
5006 ! ALL DONE !!!!!!!!!!!
5008 END SUBROUTINE MORR_TWO_MOMENT_MICRO
5010 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
5012 REAL FUNCTION POLYSVP (T,TYPE)
5014 !-------------------------------------------
5016 ! COMPUTE SATURATION VAPOR PRESSURE
5018 ! POLYSVP RETURNED IN UNITS OF PA.
5019 ! T IS INPUT IN UNITS OF K.
5020 ! TYPE REFERS TO SATURATION WITH RESPECT TO LIQUID (0) OR ICE (1)
5022 ! REPLACE GOFF-GRATCH WITH FASTER FORMULATION FROM FLATAU ET AL. 1992, TABLE 4 (RIGHT-HAND COLUMN)
5025 SAVE !TWG 2017 add to avoid memory issues
5031 real a0i,a1i,a2i,a3i,a4i,a5i,a6i,a7i,a8i
5032 data a0i,a1i,a2i,a3i,a4i,a5i,a6i,a7i,a8i /&
5033 6.11147274, 0.503160820, 0.188439774e-1, &
5034 0.420895665e-3, 0.615021634e-5,0.602588177e-7, &
5035 0.385852041e-9, 0.146898966e-11, 0.252751365e-14/
5038 real a0,a1,a2,a3,a4,a5,a6,a7,a8
5041 data a0,a1,a2,a3,a4,a5,a6,a7,a8 /&
5042 6.11239921, 0.443987641, 0.142986287e-1, &
5043 0.264847430e-3, 0.302950461e-5, 0.206739458e-7, &
5044 0.640689451e-10,-0.952447341e-13,-0.976195544e-15/
5051 ! POLYSVP = 10.**(-9.09718*(273.16/T-1.)-3.56654* &
5052 ! LOG10(273.16/T)+0.876793*(1.-T/273.16)+ &
5053 ! LOG10(6.1071))*100.
5054 ! hm 11/16/20, use Goff-Gratch for T < 195.8 K and Flatau et al. equal or above 195.8 K
5055 if (t.ge.195.8) then
5057 polysvp = a0i + dt*(a1i+dt*(a2i+dt*(a3i+dt*(a4i+dt*(a5i+dt*(a6i+dt*(a7i+a8i*dt)))))))
5058 polysvp = polysvp*100.
5061 polysvp = 10.**(-9.09718*(273.16/t-1.)-3.56654* &
5062 alog10(273.16/t)+0.876793*(1.-t/273.16)+ &
5063 alog10(6.1071))*100.
5073 ! POLYSVP = 10.**(-7.90298*(373.16/T-1.)+ &
5074 ! 5.02808*LOG10(373.16/T)- &
5075 ! 1.3816E-7*(10**(11.344*(1.-T/373.16))-1.)+ &
5076 ! 8.1328E-3*(10**(-3.49149*(373.16/T-1.))-1.)+ &
5077 ! LOG10(1013.246))*100.
5078 ! hm 11/16/20, use Goff-Gratch for T < 202.0 K and Flatau et al. equal or above 202.0 K
5079 if (t.ge.202.0) then
5081 polysvp = a0 + dt*(a1+dt*(a2+dt*(a3+dt*(a4+dt*(a5+dt*(a6+dt*(a7+a8*dt)))))))
5082 polysvp = polysvp*100.
5085 ! note: uncertain below -70 C, but produces physical values (non-negative) unlike flatau
5086 polysvp = 10.**(-7.90298*(373.16/t-1.)+ &
5087 5.02808*alog10(373.16/t)- &
5088 1.3816e-7*(10**(11.344*(1.-t/373.16))-1.)+ &
5089 8.1328e-3*(10**(-3.49149*(373.16/t-1.))-1.)+ &
5090 alog10(1013.246))*100.
5096 END FUNCTION POLYSVP
5098 !------------------------------------------------------------------------------
5100 REAL FUNCTION GAMMA(X)
5101 !----------------------------------------------------------------------
5103 ! THIS ROUTINE CALCULATES THE GAMMA FUNCTION FOR A REAL ARGUMENT X.
5104 ! COMPUTATION IS BASED ON AN ALGORITHM OUTLINED IN REFERENCE 1.
5105 ! THE PROGRAM USES RATIONAL FUNCTIONS THAT APPROXIMATE THE GAMMA
5106 ! FUNCTION TO AT LEAST 20 SIGNIFICANT DECIMAL DIGITS. COEFFICIENTS
5107 ! FOR THE APPROXIMATION OVER THE INTERVAL (1,2) ARE UNPUBLISHED.
5108 ! THOSE FOR THE APPROXIMATION FOR X .GE. 12 ARE FROM REFERENCE 2.
5109 ! THE ACCURACY ACHIEVED DEPENDS ON THE ARITHMETIC SYSTEM, THE
5110 ! COMPILER, THE INTRINSIC FUNCTIONS, AND PROPER SELECTION OF THE
5111 ! MACHINE-DEPENDENT CONSTANTS.
5114 !*******************************************************************
5115 !*******************************************************************
5117 ! EXPLANATION OF MACHINE-DEPENDENT CONSTANTS
5119 ! BETA - RADIX FOR THE FLOATING-POINT REPRESENTATION
5120 ! MAXEXP - THE SMALLEST POSITIVE POWER OF BETA THAT OVERFLOWS
5121 ! XBIG - THE LARGEST ARGUMENT FOR WHICH GAMMA(X) IS REPRESENTABLE
5122 ! IN THE MACHINE, I.E., THE SOLUTION TO THE EQUATION
5123 ! GAMMA(XBIG) = BETA**MAXEXP
5124 ! XINF - THE LARGEST MACHINE REPRESENTABLE FLOATING-POINT NUMBER;
5125 ! APPROXIMATELY BETA**MAXEXP
5126 ! EPS - THE SMALLEST POSITIVE FLOATING-POINT NUMBER SUCH THAT
5128 ! XMININ - THE SMALLEST POSITIVE FLOATING-POINT NUMBER SUCH THAT
5129 ! 1/XMININ IS MACHINE REPRESENTABLE
5131 ! APPROXIMATE VALUES FOR SOME IMPORTANT MACHINES ARE:
5135 ! CRAY-1 (S.P.) 2 8191 966.961
5137 ! UNDER NOS (S.P.) 2 1070 177.803
5139 ! SUN, ETC.) (S.P.) 2 128 35.040
5141 ! SUN, ETC.) (D.P.) 2 1024 171.624
5142 ! IBM 3033 (D.P.) 16 63 57.574
5143 ! VAX D-FORMAT (D.P.) 2 127 34.844
5144 ! VAX G-FORMAT (D.P.) 2 1023 171.489
5148 ! CRAY-1 (S.P.) 5.45E+2465 7.11E-15 1.84E-2466
5150 ! UNDER NOS (S.P.) 1.26E+322 3.55E-15 3.14E-294
5152 ! SUN, ETC.) (S.P.) 3.40E+38 1.19E-7 1.18E-38
5154 ! SUN, ETC.) (D.P.) 1.79D+308 2.22D-16 2.23D-308
5155 ! IBM 3033 (D.P.) 7.23D+75 2.22D-16 1.39D-76
5156 ! VAX D-FORMAT (D.P.) 1.70D+38 1.39D-17 5.88D-39
5157 ! VAX G-FORMAT (D.P.) 8.98D+307 1.11D-16 1.12D-308
5159 !*******************************************************************
5160 !*******************************************************************
5164 ! THE PROGRAM RETURNS THE VALUE XINF FOR SINGULARITIES OR
5165 ! WHEN OVERFLOW WOULD OCCUR. THE COMPUTATION IS BELIEVED
5166 ! TO BE FREE OF UNDERFLOW AND OVERFLOW.
5169 ! INTRINSIC FUNCTIONS REQUIRED ARE:
5171 ! INT, DBLE, EXP, LOG, REAL, SIN
5174 ! REFERENCES: AN OVERVIEW OF SOFTWARE DEVELOPMENT FOR SPECIAL
5175 ! FUNCTIONS W. J. CODY, LECTURE NOTES IN MATHEMATICS,
5176 ! 506, NUMERICAL ANALYSIS DUNDEE, 1975, G. A. WATSON
5177 ! (ED.), SPRINGER VERLAG, BERLIN, 1976.
5179 ! COMPUTER APPROXIMATIONS, HART, ET. AL., WILEY AND
5180 ! SONS, NEW YORK, 1968.
5182 ! LATEST MODIFICATION: OCTOBER 12, 1989
5184 ! AUTHORS: W. J. CODY AND L. STOLTZ
5185 ! APPLIED MATHEMATICS DIVISION
5186 ! ARGONNE NATIONAL LABORATORY
5189 !----------------------------------------------------------------------
5194 CONV,EPS,FACT,HALF,ONE,RES,SUM,TWELVE, &
5195 TWO,X,XBIG,XDEN,XINF,XMININ,XNUM,Y,Y1,YSQ,Z,ZERO
5196 REAL, DIMENSION(7) :: C
5197 REAL, DIMENSION(8) :: P
5198 REAL, DIMENSION(8) :: Q
5199 !----------------------------------------------------------------------
5200 ! MATHEMATICAL CONSTANTS
5201 !----------------------------------------------------------------------
5202 DATA ONE,HALF,TWELVE,TWO,ZERO/1.0E0,0.5E0,12.0E0,2.0E0,0.0E0/
5205 !----------------------------------------------------------------------
5206 ! MACHINE DEPENDENT PARAMETERS
5207 !----------------------------------------------------------------------
5208 DATA XBIG,XMININ,EPS/35.040E0,1.18E-38,1.19E-7/,XINF/3.4E38/
5209 !----------------------------------------------------------------------
5210 ! NUMERATOR AND DENOMINATOR COEFFICIENTS FOR RATIONAL MINIMAX
5211 ! APPROXIMATION OVER (1,2).
5212 !----------------------------------------------------------------------
5213 DATA P/-1.71618513886549492533811E+0,2.47656508055759199108314E+1, &
5214 -3.79804256470945635097577E+2,6.29331155312818442661052E+2, &
5215 8.66966202790413211295064E+2,-3.14512729688483675254357E+4, &
5216 -3.61444134186911729807069E+4,6.64561438202405440627855E+4/
5217 DATA Q/-3.08402300119738975254353E+1,3.15350626979604161529144E+2, &
5218 -1.01515636749021914166146E+3,-3.10777167157231109440444E+3, &
5219 2.25381184209801510330112E+4,4.75584627752788110767815E+3, &
5220 -1.34659959864969306392456E+5,-1.15132259675553483497211E+5/
5221 !----------------------------------------------------------------------
5222 ! COEFFICIENTS FOR MINIMAX APPROXIMATION OVER (12, INF).
5223 !----------------------------------------------------------------------
5224 DATA C/-1.910444077728E-03,8.4171387781295E-04, &
5225 -5.952379913043012E-04,7.93650793500350248E-04, &
5226 -2.777777777777681622553E-03,8.333333333333333331554247E-02, &
5228 !----------------------------------------------------------------------
5229 ! STATEMENT FUNCTIONS FOR CONVERSION BETWEEN INTEGER AND FLOAT
5230 !----------------------------------------------------------------------
5237 !----------------------------------------------------------------------
5238 ! ARGUMENT IS NEGATIVE
5239 !----------------------------------------------------------------------
5244 IF(Y1.NE.AINT(Y1*HALF)*TWO)PARITY=.TRUE.
5245 FACT=-PI/SIN(PI*RES)
5252 !----------------------------------------------------------------------
5253 ! ARGUMENT IS POSITIVE
5254 !----------------------------------------------------------------------
5256 !----------------------------------------------------------------------
5258 !----------------------------------------------------------------------
5265 ELSEIF(Y.LT.TWELVE)THEN
5268 !----------------------------------------------------------------------
5269 ! 0.0 .LT. ARGUMENT .LT. 1.0
5270 !----------------------------------------------------------------------
5274 !----------------------------------------------------------------------
5275 ! 1.0 .LT. ARGUMENT .LT. 12.0, REDUCE ARGUMENT IF NECESSARY
5276 !----------------------------------------------------------------------
5281 !----------------------------------------------------------------------
5282 ! EVALUATE APPROXIMATION FOR 1.0 .LT. ARGUMENT .LT. 2.0
5283 !----------------------------------------------------------------------
5292 !----------------------------------------------------------------------
5293 ! ADJUST RESULT FOR CASE 0.0 .LT. ARGUMENT .LT. 1.0
5294 !----------------------------------------------------------------------
5297 !----------------------------------------------------------------------
5298 ! ADJUST RESULT FOR CASE 2.0 .LT. ARGUMENT .LT. 12.0
5299 !----------------------------------------------------------------------
5306 !----------------------------------------------------------------------
5307 ! EVALUATE FOR ARGUMENT .GE. 12.0,
5308 !----------------------------------------------------------------------
5316 SUM=SUM+(Y-HALF)*LOG(Y)
5323 !----------------------------------------------------------------------
5324 ! FINAL ADJUSTMENTS AND RETURN
5325 !----------------------------------------------------------------------
5327 IF(FACT.NE.ONE)RES=FACT/RES
5330 ! ---------- LAST LINE OF GAMMA ----------
5334 REAL FUNCTION DERF1(X)
5337 REAL, DIMENSION(0 : 64) :: A, B
5341 0.00000000005958930743E0, -0.00000000113739022964E0, &
5342 0.00000001466005199839E0, -0.00000016350354461960E0, &
5343 0.00000164610044809620E0, -0.00001492559551950604E0, &
5344 0.00012055331122299265E0, -0.00085483269811296660E0, &
5345 0.00522397762482322257E0, -0.02686617064507733420E0, &
5346 0.11283791670954881569E0, -0.37612638903183748117E0, &
5347 1.12837916709551257377E0, &
5348 0.00000000002372510631E0, -0.00000000045493253732E0, &
5349 0.00000000590362766598E0, -0.00000006642090827576E0, &
5350 0.00000067595634268133E0, -0.00000621188515924000E0, &
5351 0.00005103883009709690E0, -0.00037015410692956173E0, &
5352 0.00233307631218880978E0, -0.01254988477182192210E0, &
5353 0.05657061146827041994E0, -0.21379664776456006580E0, &
5354 0.84270079294971486929E0, &
5355 0.00000000000949905026E0, -0.00000000018310229805E0, &
5356 0.00000000239463074000E0, -0.00000002721444369609E0, &
5357 0.00000028045522331686E0, -0.00000261830022482897E0, &
5358 0.00002195455056768781E0, -0.00016358986921372656E0, &
5359 0.00107052153564110318E0, -0.00608284718113590151E0, &
5360 0.02986978465246258244E0, -0.13055593046562267625E0, &
5361 0.67493323603965504676E0, &
5362 0.00000000000382722073E0, -0.00000000007421598602E0, &
5363 0.00000000097930574080E0, -0.00000001126008898854E0, &
5364 0.00000011775134830784E0, -0.00000111992758382650E0, &
5365 0.00000962023443095201E0, -0.00007404402135070773E0, &
5366 0.00050689993654144881E0, -0.00307553051439272889E0, &
5367 0.01668977892553165586E0, -0.08548534594781312114E0, &
5368 0.56909076642393639985E0, &
5369 0.00000000000155296588E0, -0.00000000003032205868E0, &
5370 0.00000000040424830707E0, -0.00000000471135111493E0, &
5371 0.00000005011915876293E0, -0.00000048722516178974E0, &
5372 0.00000430683284629395E0, -0.00003445026145385764E0, &
5373 0.00024879276133931664E0, -0.00162940941748079288E0, &
5374 0.00988786373932350462E0, -0.05962426839442303805E0, &
5375 0.49766113250947636708E0 /
5376 DATA (B(I), I = 0, 12) / &
5377 -0.00000000029734388465E0, 0.00000000269776334046E0, &
5378 -0.00000000640788827665E0, -0.00000001667820132100E0, &
5379 -0.00000021854388148686E0, 0.00000266246030457984E0, &
5380 0.00001612722157047886E0, -0.00025616361025506629E0, &
5381 0.00015380842432375365E0, 0.00815533022524927908E0, &
5382 -0.01402283663896319337E0, -0.19746892495383021487E0, &
5383 0.71511720328842845913E0 /
5384 DATA (B(I), I = 13, 25) / &
5385 -0.00000000001951073787E0, -0.00000000032302692214E0, &
5386 0.00000000522461866919E0, 0.00000000342940918551E0, &
5387 -0.00000035772874310272E0, 0.00000019999935792654E0, &
5388 0.00002687044575042908E0, -0.00011843240273775776E0, &
5389 -0.00080991728956032271E0, 0.00661062970502241174E0, &
5390 0.00909530922354827295E0, -0.20160072778491013140E0, &
5391 0.51169696718727644908E0 /
5392 DATA (B(I), I = 26, 38) / &
5393 0.00000000003147682272E0, -0.00000000048465972408E0, &
5394 0.00000000063675740242E0, 0.00000003377623323271E0, &
5395 -0.00000015451139637086E0, -0.00000203340624738438E0, &
5396 0.00001947204525295057E0, 0.00002854147231653228E0, &
5397 -0.00101565063152200272E0, 0.00271187003520095655E0, &
5398 0.02328095035422810727E0, -0.16725021123116877197E0, &
5399 0.32490054966649436974E0 /
5400 DATA (B(I), I = 39, 51) / &
5401 0.00000000002319363370E0, -0.00000000006303206648E0, &
5402 -0.00000000264888267434E0, 0.00000002050708040581E0, &
5403 0.00000011371857327578E0, -0.00000211211337219663E0, &
5404 0.00000368797328322935E0, 0.00009823686253424796E0, &
5405 -0.00065860243990455368E0, -0.00075285814895230877E0, &
5406 0.02585434424202960464E0, -0.11637092784486193258E0, &
5407 0.18267336775296612024E0 /
5408 DATA (B(I), I = 52, 64) / &
5409 -0.00000000000367789363E0, 0.00000000020876046746E0, &
5410 -0.00000000193319027226E0, -0.00000000435953392472E0, &
5411 0.00000018006992266137E0, -0.00000078441223763969E0, &
5412 -0.00000675407647949153E0, 0.00008428418334440096E0, &
5413 -0.00017604388937031815E0, -0.00239729611435071610E0, &
5414 0.02064129023876022970E0, -0.06905562880005864105E0, &
5415 0.09084526782065478489E0 /
5417 IF (W .LT. 2.2D0) THEN
5422 Y = ((((((((((((A(K) * T + A(K + 1)) * T + &
5423 A(K + 2)) * T + A(K + 3)) * T + A(K + 4)) * T + &
5424 A(K + 5)) * T + A(K + 6)) * T + A(K + 7)) * T + &
5425 A(K + 8)) * T + A(K + 9)) * T + A(K + 10)) * T + &
5426 A(K + 11)) * T + A(K + 12)) * W
5427 ELSE IF (W .LT. 6.9D0) THEN
5431 Y = (((((((((((B(K) * T + B(K + 1)) * T + &
5432 B(K + 2)) * T + B(K + 3)) * T + B(K + 4)) * T + &
5433 B(K + 5)) * T + B(K + 6)) * T + B(K + 7)) * T + &
5434 B(K + 8)) * T + B(K + 9)) * T + B(K + 10)) * T + &
5435 B(K + 11)) * T + B(K + 12)
5443 IF (X .LT. 0) Y = -Y
5447 !+---+-----------------------------------------------------------------+
5449 subroutine refl10cm_hm (qv1d, qr1d, nr1d, qs1d, ns1d, qg1d, ng1d, &
5450 t1d, p1d, dBZ, kts, kte, ii, jj)
5453 SAVE !TWG 2017 add to avoid memeory issues
5456 INTEGER, INTENT(IN):: kts, kte, ii, jj
5457 REAL, DIMENSION(kts:kte), INTENT(IN):: &
5458 qv1d, qr1d, nr1d, qs1d, ns1d, qg1d, ng1d, t1d, p1d
5459 REAL, DIMENSION(kts:kte), INTENT(INOUT):: dBZ
5462 REAL, DIMENSION(kts:kte):: temp, pres, qv, rho
5463 REAL, DIMENSION(kts:kte):: rr, nr, rs, ns, rg, ng
5465 DOUBLE PRECISION, DIMENSION(kts:kte):: ilamr, ilamg, ilams
5466 DOUBLE PRECISION, DIMENSION(kts:kte):: N0_r, N0_g, N0_s
5467 DOUBLE PRECISION:: lamr, lamg, lams
5468 LOGICAL, DIMENSION(kts:kte):: L_qr, L_qs, L_qg
5470 REAL, DIMENSION(kts:kte):: ze_rain, ze_snow, ze_graupel
5471 DOUBLE PRECISION:: fmelt_s, fmelt_g
5472 DOUBLE PRECISION:: cback, x, eta, f_d
5474 INTEGER:: i, k, k_0, kbot, n
5483 !+---+-----------------------------------------------------------------+
5484 !..Put column of data into local arrays.
5485 !+---+-----------------------------------------------------------------+
5488 qv(k) = MAX(1.E-10, qv1d(k))
5490 rho(k) = 0.622*pres(k)/(R*temp(k)*(qv(k)+0.622))
5492 if (qr1d(k) .gt. 1.E-9) then
5493 rr(k) = qr1d(k)*rho(k)
5494 nr(k) = nr1d(k)*rho(k)
5495 lamr = (xam_r*xcrg(3)*xorg2*nr(k)/rr(k))**xobmr
5497 N0_r(k) = nr(k)*xorg2*lamr**xcre(2)
5505 if (qs1d(k) .gt. 1.E-9) then
5506 rs(k) = qs1d(k)*rho(k)
5507 ns(k) = ns1d(k)*rho(k)
5508 lams = (xam_s*xcsg(3)*xosg2*ns(k)/rs(k))**xobms
5510 N0_s(k) = ns(k)*xosg2*lams**xcse(2)
5518 if (qg1d(k) .gt. 1.E-9) then
5519 rg(k) = qg1d(k)*rho(k)
5520 ng(k) = ng1d(k)*rho(k)
5521 lamg = (xam_g*xcgg(3)*xogg2*ng(k)/rg(k))**xobmg
5523 N0_g(k) = ng(k)*xogg2*lamg**xcge(2)
5532 !+---+-----------------------------------------------------------------+
5533 !..Locate K-level of start of melting (k_0 is level above).
5534 !+---+-----------------------------------------------------------------+
5537 do k = kte-1, kts, -1
5538 if ( (temp(k).gt.273.15) .and. L_qr(k) &
5539 .and. (L_qs(k+1).or.L_qg(k+1)) ) then
5547 !+---+-----------------------------------------------------------------+
5548 !..Assume Rayleigh approximation at 10 cm wavelength. Rain (all temps)
5549 !.. and non-water-coated snow and graupel when below freezing are
5550 !.. simple. Integrations of m(D)*m(D)*N(D)*dD.
5551 !+---+-----------------------------------------------------------------+
5556 ze_graupel(k) = 1.e-22
5557 if (L_qr(k)) ze_rain(k) = N0_r(k)*xcrg(4)*ilamr(k)**xcre(4)
5558 if (L_qs(k)) ze_snow(k) = (0.176/0.93) * (6.0/PI)*(6.0/PI) &
5559 * (xam_s/900.0)*(xam_s/900.0) &
5560 * N0_s(k)*xcsg(4)*ilams(k)**xcse(4)
5561 if (L_qg(k)) ze_graupel(k) = (0.176/0.93) * (6.0/PI)*(6.0/PI) &
5562 * (xam_g/900.0)*(xam_g/900.0) &
5563 * N0_g(k)*xcgg(4)*ilamg(k)**xcge(4)
5566 !+---+-----------------------------------------------------------------+
5567 !..Special case of melting ice (snow/graupel) particles. Assume the
5568 !.. ice is surrounded by the liquid water. Fraction of meltwater is
5569 !.. extremely simple based on amount found above the melting level.
5570 !.. Uses code from Uli Blahak (rayleigh_soak_wetgraupel and supporting
5572 !+---+-----------------------------------------------------------------+
5574 if (melti .and. k_0.ge.kts+1) then
5575 do k = k_0-1, kts, -1
5577 !..Reflectivity contributed by melting snow
5578 if (L_qs(k) .and. L_qs(k_0) ) then
5579 fmelt_s = MAX(0.005d0, MIN(1.0d0-rs(k)/rs(k_0), 0.99d0))
5583 x = xam_s * xxDs(n)**xbm_s
5584 call rayleigh_soak_wetgraupel (x,DBLE(xocms),DBLE(xobms), &
5585 fmelt_s, melt_outside_s, m_w_0, m_i_0, lamda_radar, &
5586 CBACK, mixingrulestring_s, matrixstring_s, &
5587 inclusionstring_s, hoststring_s, &
5588 hostmatrixstring_s, hostinclusionstring_s)
5589 f_d = N0_s(k)*xxDs(n)**xmu_s * DEXP(-lams*xxDs(n))
5590 eta = eta + f_d * CBACK * simpson(n) * xdts(n)
5592 ze_snow(k) = SNGL(lamda4 / (pi5 * K_w) * eta)
5596 !..Reflectivity contributed by melting graupel
5598 if (L_qg(k) .and. L_qg(k_0) ) then
5599 fmelt_g = MAX(0.005d0, MIN(1.0d0-rg(k)/rg(k_0), 0.99d0))
5603 x = xam_g * xxDg(n)**xbm_g
5604 call rayleigh_soak_wetgraupel (x,DBLE(xocmg),DBLE(xobmg), &
5605 fmelt_g, melt_outside_g, m_w_0, m_i_0, lamda_radar, &
5606 CBACK, mixingrulestring_g, matrixstring_g, &
5607 inclusionstring_g, hoststring_g, &
5608 hostmatrixstring_g, hostinclusionstring_g)
5609 f_d = N0_g(k)*xxDg(n)**xmu_g * DEXP(-lamg*xxDg(n))
5610 eta = eta + f_d * CBACK * simpson(n) * xdtg(n)
5612 ze_graupel(k) = SNGL(lamda4 / (pi5 * K_w) * eta)
5619 dBZ(k) = 10.*log10((ze_rain(k)+ze_snow(k)+ze_graupel(k))*1.d18)
5623 end subroutine refl10cm_hm
5625 !+---+-----------------------------------------------------------------+
5626 ! TWG 2016 Begin Add activation subroutine for prescribed aerosols
5628 subroutine mdm_prescribed_activate(wbar, tair, rhoair, &
5629 naero, pmode, nmode, maero, sigman, hygro, rhodry, nact,latvap)
5630 ! calculates number, surface, and mass fraction of aerosols activated as
5632 ! calculates flux of cloud droplets, surface area, and aerosol mass into
5634 ! assumes an internal mixture within each of up to pmode multiple aerosol
5636 ! a gaussiam spectrum of updrafts can be treated.
5640 ! Abdul-Razzak and Ghan, A parameterization of aerosol activation.
5641 ! 2. Multiple aerosol types. J. Geophys. Res., 105, 6837-6844.
5643 USE module_model_constants, only: svp1,svp2,svp3,ep_2
5652 integer pmode,ptype ! dimension of modes, types in modes
5653 real wbar ! grid cell mean vertical velocity (m/s)
5654 real tair ! air temperature (K)
5655 real rhoair ! air density (kg/m3)
5656 real naero(pmode) ! aerosol number concentration (/m3)
5657 integer nmode ! number of aerosol modes
5658 real maero(pmode) ! aerosol mass concentration (kg/m3)
5659 real rhodry(pmode) ! density of aerosol material
5660 real sigman(pmode) ! geometric standard deviation of aerosol size distribution
5661 real hygro(pmode) ! hygroscopicity of aerosol mode
5666 real nact ! number fraction of aerosols activated
5669 ! real derf,derfc, erf_alt
5670 integer, parameter:: nx=200
5673 real p0 ! reference pressure (Pa)
5677 real :: volc(naer_cu) ! total aerosol volume concentration (m3/m3)
5678 real tmass ! total aerosol mass concentration (g/cm3)
5679 real rm ! number mode radius of aerosol at max supersat (cm)
5680 real pres ! pressure (Pa)
5681 real path ! mean free path (m)
5682 real diff ! diffusivity (m2/s)
5683 real conduct ! thermal conductivity (Joule/m/sec/deg)
5685 real qs ! water vapor saturation mixing ratio
5686 real dqsdt ! change in qs with temperature
5687 real dqsdp ! change in qs with pressure
5688 real gloc ! thermodynamic function (m2/s)
5690 real :: eta(naer_cu)
5691 real :: smc(naer_cu)
5692 real lnsmax ! ln(smax)
5698 real rlo,rhi,xint1,xint2,xint3,xint4
5703 real xmincoeff,xcut,volcut,surfcut
5704 real z,z1,z2,wf1,wf2,zf1,zf2,gf1,gf2,gf
5705 real :: etafactor1,etafactor2max
5706 real :: etafactor2(naer_cu)
5711 real :: amcubeloc(naer_cu)
5712 real :: lnsmloc(naer_cu)
5718 !print*,'Doing MSKF-AERO Activation MDM'
5720 if(nmode.eq.1.and.naero(1).lt.1.e-20)return
5722 if(wbar.le.0.)return
5725 diff0=0.211e-4*(p0/pres)*(tair/t0)**1.94
5726 conduct0=(5.69+0.017*(tair-t0))*4.186e2*1.e-5 ! convert to J/m/s/deg
5727 es=1000.*svp1*exp( svp2*(tair-t0)/(tair-svp3) )
5728 qs=ep_2*es/(pres-es)
5729 dqsdt=latvap/(RV*tair*tair)*qs
5730 alpha_cesm=G*(latvap/(CP*RV*tair*tair)-1./(R*tair))
5731 gammaloc=(1+latvap/CP*dqsdt)/(rhoair*qs)
5732 ! growth coefficent Abdul-Razzak & Ghan 1998 eqn 16
5733 ! should depend on mean radius of mode to account for gas kinetic effects
5734 gloc=1./(RHOW/(diff0*rhoair*qs) &
5735 +latvap*RHOW/(conduct0*tair)*(latvap/(RV*tair)-1.))
5737 beta=4.*pi*RHOW*gloc*gammaloc
5738 etafactor2max=1.e10/(alpha_cesm*wbar)**1.5 ! this should make eta big if na is very small.
5741 ! internal mixture of aerosols
5742 volc(m)=maero(m)/(rhodry(m)) ! only if variable size dist
5743 if(volc(m).gt.real_tiny.and.naero(m).gt.real_tiny)then
5744 etafactor2(m)=1./(naero(m)*beta*sqrtg) !fixed or variable size dist
5745 ! number mode radius (m)
5746 amcubeloc(m)=(3.*volc(m)/(4.*pi*exp45logsig_pamdm(m)*naero(m))) !only if variable size dist
5747 smc(m)=smcrit_pamdm(m) ! only for prescribed size dist
5752 if(hygro(m).gt.1.e-10)then ! loop only if variable size dist
5753 smc(m)=2.*aten_pamdm*sqrt(aten_pamdm/(27.*hygro(m)*amcubeloc(m)))
5760 etafactor2(m)=etafactor2max ! this should make eta big if na is very small.
5762 lnsmloc(m)=log(smc(m)) ! only if variable size dist
5771 zeta=2.*sqrtalw*aten_pamdm/(3.*sqrtg)
5772 etafactor1=2.*alw*sqrtalw
5775 eta(m)=etafactor1*etafactor2(m)
5778 call mdm_prescribed_maxsat(zeta,eta,nmode,smc,smax)
5781 xmincoeff=alogaten_pamdm-2.*third_pamdm*(lnsmax-alog2_pamdm)-alog3_pamdm
5785 x=2*(lnsmloc(m)-lnsmax)/(3*sq2_pamdm*alogsig_pamdm(m))
5786 nact=nact+0.5*(1.-DERF1(x))*naero(m) ! danger erf
5788 nact=nact/rhoair ! convert from #/m3 to #/kg
5791 end subroutine mdm_prescribed_activate
5793 subroutine mdm_prescribed_maxsat(zeta,eta,nmode,smc,smax)
5795 ! calculates maximum supersaturation for multiple
5796 ! competing aerosol modes.
5798 ! Abdul-Razzak and Ghan, A parameterization of aerosol activation.
5799 ! 2. Multiple aerosol types. J. Geophys. Res., 105, 6837-6844.
5803 integer nmode ! number of modes
5804 real :: smc(:) ! critical supersaturation for number mode radius
5807 real smax ! maximum supersaturation
5808 integer m ! mode index
5812 if(zeta.gt.1.e5*eta(m).or.smc(m)*smc(m).gt.1.e5*eta(m))then
5813 ! weak forcing. essentially none activated
5816 ! significant activation of this mode. calc activation all modes.
5827 if(eta(m).gt.1.e-20)then
5828 g1=sqrt(zeta/eta(m))
5830 g2=smc(m)/sqrt(eta(m)+3*zeta)
5833 sum=sum+(f1_pamdm(m)*g1+f2_pamdm(m)*g2)/(smc(m)*smc(m))
5844 end subroutine mdm_prescribed_maxsat
5846 subroutine mdm_prescribed_nucleati(wbar, tair, relhum, relhumi, qc, rhoair, &
5847 naero, naer_all, nuci &
5848 , onihf, oniimm, onidep, onimey)
5850 !---------------------------------------------------------------
5852 ! The parameterization of ice nucleation.
5854 ! Method: The current method is based on Liu & Penner (2005)
5855 ! It related the ice nucleation with the aerosol number, temperature and the
5856 ! updraft velocity. It includes homogeneous freezing of sulfate, immersion
5857 ! freezing of soot, and Meyers et al. (1992) deposition nucleation
5859 ! Authors: Xiaohong Liu, 01/2005, modifications by A. Gettelman 2009-2010
5860 !----------------------------------------------------------------
5865 real :: wbar ! grid cell mean vertical velocity (m/s)
5866 real :: tair ! temperature (K)
5867 real :: relhum ! relative humidity with respoect to water
5868 real :: relhumi ! relative humidity with respective to ice
5870 real :: qc ! liquid water mixing ratio (kg/kg)
5871 real :: rhoair ! air density (kg/m3)
5872 real :: naero(naer_all) ! aerosol number concentration (/m3)
5877 real :: nuci ! ice number nucleated (#/kg)
5878 real :: onihf ! nucleated number from homogeneous freezing of so4
5879 real :: oniimm ! nucleated number from immersion freezing
5880 real :: onidep ! nucleated number from deposition nucleation
5881 real :: onimey ! nucleated number from deposition nucleation (meyers: mixed phase)
5886 real so4_num ! so4 aerosol number(#/cm^3)
5887 real soot_num ! soot (hydrophilic) aerosol number (#/cm^3)
5888 real dst1_num,dst2_num,dst3_num,dst4_num ! dust aerosol number (#/cm^3)
5889 real dst_num ! total dust aerosol number (#/cm^3)
5890 real nihf ! nucleated number from homogeneous freezing of so4
5891 real niimm ! nucleated number from immersion freezing
5892 real nidep ! nucleated number from deposition nucleation
5893 real nimey ! nucleated number from deposition nucleation (meyers)
5894 real n1,ni ! nucleated number
5895 real tc,A,B,C,regm ! work variable
5896 real esl,esi,deles ! work variable
5899 real dmc,ssmc ! variables for modal scheme.
5902 !print*,'Doing MSKF Aerosol Ice Nucleation MDM'
5912 !For modal aerosols, assume for the upper troposphere:
5913 ! soot = accumulation mode
5914 ! sulfate = aiken mode
5915 ! dust = coarse mode
5916 ! since modal has internal mixtures.
5918 so4_num=naero(1)*1.0e-6 ! #/cm^3
5919 soot_num=naero(10)*1.0e-6 !#/cm^3
5920 dst1_num=naero(3)*1.0e-6 !#/cm^3
5921 dst2_num=naero(4)*1.0e-6 !#/cm^3
5922 dst3_num=naero(5)*1.0e-6 !#/cm^3
5923 dst4_num=naero(6)*1.0e-6 !#/cm^3
5925 dst_num =dst1_num+dst2_num+dst3_num+dst4_num
5926 ! no soot nucleation for now.
5937 if(so4_num.ge.1.0e-10 .and. (soot_num+dst_num).ge.1.0e-10 ) then
5941 if((tc.le.-35.0).and.((relhumi*subgrid).ge.1.2))then !Make saturation consistent with Morrison Scheme
5942 !< regm => T in Eq.10 of Liu et al., J. Climate, 2007>
5943 A = -1.4938 * log(soot_num+dst_num) + 12.884
5944 B = -10.41 * log(soot_num+dst_num) - 67.69
5945 regm = A * log(wbar) + B
5947 if(tc.gt.regm) then ! heterogeneous nucleation only
5948 if(tc.lt.-40. .and. wbar.gt.1.) then ! exclude T<-40 & W>1m/s from hetero. nucleation
5949 call mdm_prescribed_hf(tc,wbar,relhum,subgrid,so4_num,nihf)
5954 call mdm_prescribed_hetero(tc,wbar,soot_num+dst_num,niimm,nidep)
5958 elseif (tc.lt.regm-5.) then ! homogeneous nucleation only
5959 call mdm_prescribed_hf(tc,wbar,relhum,subgrid,so4_num,nihf)
5963 else ! transition between homogeneous and heterogeneous: interpolate in-between
5964 if(tc.lt.-40. .and. wbar.gt.1.) then ! exclude T<-40 & W>1m/s from hetero. nucleation
5965 call mdm_prescribed_hf(tc,wbar,relhum,subgrid,so4_num,nihf)
5970 call mdm_prescribed_hf(regm-5.,wbar,relhum,subgrid,so4_num,nihf)
5971 call mdm_prescribed_hetero(regm,wbar,soot_num+dst_num,niimm,nidep)
5972 if(nihf.le.(niimm+nidep)) then
5975 n1=(niimm+nidep)*((niimm+nidep)/nihf)**((tc-regm)/5.)
5986 ! deposition/condensation nucleation in mixed clouds (-37<T<0C) (Meyers, 1992)
5987 !<Eq.12 of Liu et al., J. Climate, 2007
5988 ! Nid(L-1)*1.e-3 => Nid(m-3)
5989 ! Question: RHi=RHw*esl/esi
5991 if(tc.lt.0. .and. tc.gt.-37. .and. qc.gt.1.e-12) then
5992 ! if(tc.lt.0. .and. tc.gt.-37) then
5993 ! esl = mdm_prescribed_polysvp(tair,0) ! over water in mixed clouds
5994 ! esi = mdm_prescribed_polysvp(tair,1) ! over ice
5995 ! deles = (esl - esi)
5996 ! deles = (relhum*esl - esi) TWG use original formulation
5997 if(relhumi.gt.1.5) then ! TWG add 50% ice superaturation constraint to avoid numerical issues.
5998 deles = 1.5 ! This constraint is double the environmental applicability of Meyers 1992 (25%)
6003 nimey=1.e-3*exp(12.96*(deles-1.0) - 0.639) ! TWG Fix meyers formulation
6009 if(nuci.gt.9999..or.nuci.lt.0.) then
6010 write(*, *) 'incorrect ice nucleation number'
6011 write(*, *) ni, tair, relhum, relhumi, wbar, nihf,niimm,nidep,nimey,dst2_num,dst3_num,dst4_num
6013 stop 'nuclei prbolem?'
6016 nuci=nuci*1.e+6/rhoair ! change unit from #/cm3 to #/kg
6017 onimey=nimey*1.e+6/rhoair
6018 onidep=nidep*1.e+6/rhoair
6019 oniimm=niimm*1.e+6/rhoair
6020 onihf=nihf*1.e+6/rhoair
6023 end subroutine mdm_prescribed_nucleati
6025 subroutine mdm_prescribed_hetero(T,ww,Ns,Nis,Nid)
6030 real A11,A12,A21,A22,B11,B12,B21,B22
6034 !---------------------------------------------------------------------
6045 !<Eq.11 of Liu et al., J. Climate, 2007>
6046 ! ice from immersion nucleation (cm-3)
6048 B = (A11+B11*log(Ns)) * log(ww) + (A12+B12*log(Ns))
6051 Nis = exp(A22) * Ns**B22 * exp(B*T) * ww**C
6053 Nid = 0.0 ! don't include deposition nucleation for cirrus clouds when T<-37C
6055 end subroutine mdm_prescribed_hetero
6057 subroutine mdm_prescribed_hf(T,ww,RH,subgrid,Na,Ni)
6060 real :: T, ww, RH, subgrid, Na
6061 real, intent(out) :: Ni
6063 real A1_fast,A21_fast,A22_fast,B1_fast,B21_fast,B22_fast
6064 real A2_fast,B2_fast
6065 real C1_fast,C2_fast,k1_fast,k2_fast
6066 real A1_slow,A2_slow,B1_slow,B2_slow,B3_slow
6067 real C1_slow,C2_slow,k1_slow,k2_slow
6072 !---------------------------------------------------------------------
6073 !<Table 1 of Liu et al., J. Climate, 2007>
6077 A21_fast =-1.6387 !(T>-64 deg)
6078 A22_fast =-6.045 !(T<=-64 deg)
6080 B21_fast =-0.042 !(T>-64 deg)
6081 B22_fast =-0.112 !(T<=-64 deg)
6095 !----------------------------
6096 !<Eq.6 of Liu et al., J. Climate, 2007
6097 ! w~m/s, T~degree C, RHw~% => RHw*0.01~fraction >
6098 !RHw xiaohong's parameter
6099 A = 6.0e-4*log(ww)+6.6e-3
6100 B = 6.0e-2*log(ww)+1.052
6101 C = 1.68 *log(ww)+129.35
6102 RHw=(A*T*T+B*T+C)*0.01
6104 if((T.le.-37.0) .and. ((RH*subgrid).ge.RHw)) then
6106 !<Eq.9 of Liu et al., J. Climate, 2007>
6107 regm = 6.07*log(ww)-55.0
6109 if(T.ge.regm) then ! fast-growth regime
6119 !<Eq.7 of Liu et al., J. Climate, 2007>
6120 k1_fast = exp(A2_fast + B2_fast*T + C2_fast*log(ww))
6121 k2_fast = A1_fast+B1_fast*T+C1_fast*log(ww)
6123 Ni = k1_fast*Na**(k2_fast)
6125 else ! slow-growth regime
6126 !<Eq.7 of Liu et al., J. Climate, 2007>
6127 k1_slow = exp(A2_slow + (B2_slow+B3_slow*log(ww))*T + C2_slow*log(ww))
6128 k2_slow = A1_slow+B1_slow*T+C1_slow*log(ww)
6130 Ni = k1_slow*Na**(k2_slow)
6136 end subroutine mdm_prescribed_hf
6138 function mdm_prescribed_polysvp (T,type)
6139 ! Compute saturation vapor pressure by using
6140 ! function from Goff and Gatch (1946)
6142 ! Polysvp returned in units of pa.
6143 ! T is input in units of K.
6144 ! type refers to saturation with respect to liquid (0) or ice (1)
6149 real T,mdm_prescribed_polysvp
6157 ! Goff Gatch equation (good down to -100 C)
6159 mdm_prescribed_polysvp = 10.**(-9.09718*(273.16/t-1.)-3.56654* &
6160 log10(273.16/t)+0.876793*(1.-t/273.16)+ &
6166 ! Goff Gatch equation, uncertain below -70 C
6169 mdm_prescribed_polysvp = 10.**(-7.90298*(373.16/t-1.)+ &
6170 5.02808*log10(373.16/t)- &
6171 1.3816e-7*(10.**(11.344*(1.-t/373.16))-1.)+ &
6172 8.1328e-3*(10.**(-3.49149*(373.16/t-1.))-1.)+ &
6173 log10(1013.246))*100.
6177 end function mdm_prescribed_polysvp
6180 END MODULE module_mp_morr_two_moment_aero
6181 !+---+-----------------------------------------------------------------+
6182 !+---+-----------------------------------------------------------------+