| blob | be7fd33ce6efd9a661853c33439a060e932117ec |
1 !WRF:MODEL_LAYER:PHYSICS
2 !
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
18 ! bug fix, 5/12/10
19 ! 7) bug fix for saturation vapor pressure in low pressure, to avoid division by zero
20 ! 8) include 'EP2' WRF constant for saturation mixing ratio calculation, instead of hardwire constant
22 ! CHANGES FOR V3.3
23 ! 1) MODIFICATION FOR COUPLING WITH WRF-CHEM (PREDICTED DROPLET NUMBER CONCENTRATION) AS AN OPTION
24 ! 2) MODIFY FALLSPEED BELOW THE LOWEST LEVEL OF PRECIPITATION, WHICH PREVENTS
25 ! POTENTIAL FOR SPURIOUS ACCUMULATION OF PRECIPITATION DURING SUB-STEPPING FOR SEDIMENTATION
26 ! 3) BUG FIX TO LATENT HEAT RELEASE DUE TO COLLISIONS OF CLOUD ICE WITH RAIN
27 ! 4) CLEAN UP OF COMMENTS IN THE CODE
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
37 ! additional fixes:
38 ! 5) TEMPERATURE FOR ACCELERATED MELTING DUE TO COLLIIONS OF SNOW AND GRAUPEL
39 ! WITH RAIN SHOULD USE CELSIUS, NOT KELVIN (BUG REPORTED BY K. VAN WEVERBERG)
40 ! 6) NPRACS IS NOT SUBTRACTED FROM SNOW NUMBER CONCENTRATION, SINCE
41 ! DECREASE IN SNOW NUMBER IS ALREADY ACCOUNTED FOR BY NSMLTS
42 ! 7) fix for switch for running w/o graupel/hail (cloud ice and snow only)
44 ! hm bug fix 3/16/12
46 ! 1) very minor change to limits on autoconversion source of rain number when cloud water is depleted
48 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
50 ! THIS SCHEME IS A BULK DOUBLE-MOMENT SCHEME THAT PREDICTS MIXING
51 ! RATIOS AND NUMBER CONCENTRATIONS OF FIVE HYDROMETEOR SPECIES:
52 ! CLOUD DROPLETS, CLOUD (SMALL) ICE, RAIN, SNOW, AND GRAUPEL.
54 MODULE MODULE_MP_MORR_TWO_MOMENT
55 USE module_wrf_error
56 ! USE module_utility, ONLY: WRFU_Clock, WRFU_Alarm ! GT
57 ! USE module_domain, ONLY : HISTORY_ALARM, Is_alarm_tstep ! GT
59 ! USE WRF PHYSICS CONSTANTS
60 use module_model_constants, ONLY: CP, G, R => r_d, RV => r_v, EP_2
62 ! USE module_state_description
64 IMPLICIT NONE
66 REAL, PARAMETER :: PI = 3.1415926535897932384626434
67 REAL, PARAMETER :: SQRTPI = 0.9189385332046727417803297
69 PUBLIC :: MP_MORR_TWO_MOMENT
70 PUBLIC :: POLYSVP
72 PRIVATE :: GAMMA, DERF1
73 PRIVATE :: PI, SQRTPI
74 PRIVATE :: MORR_TWO_MOMENT_MICRO
76 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
77 ! SWITCHES FOR MICROPHYSICS SCHEME
78 ! IACT = 1, USE POWER-LAW CCN SPECTRA, NCCN = CS^K
79 ! IACT = 2, USE LOGNORMAL AEROSOL SIZE DIST TO DERIVE CCN SPECTRA
80 ! IACT = 3, ACTIVATION CALCULATED IN MODULE_MIXACTIVATE
82 INTEGER, PRIVATE :: IACT
84 ! INUM = 0, PREDICT DROPLET CONCENTRATION
85 ! INUM = 1, ASSUME CONSTANT DROPLET CONCENTRATION
86 ! !!!NOTE: PREDICTED DROPLET CONCENTRATION NOT AVAILABLE IN THIS VERSION
87 ! CONTACT HUGH MORRISON (morrison@ucar.edu) FOR FURTHER INFORMATION
89 INTEGER, PRIVATE :: INUM
91 ! FOR INUM = 1, SET CONSTANT DROPLET CONCENTRATION (CM-3)
92 REAL, PRIVATE :: NDCNST
94 ! SWITCH FOR LIQUID-ONLY RUN
95 ! ILIQ = 0, INCLUDE ICE
96 ! ILIQ = 1, LIQUID ONLY, NO ICE
98 INTEGER, PRIVATE :: ILIQ
100 ! SWITCH FOR ICE NUCLEATION
101 ! INUC = 0, USE FORMULA FROM RASMUSSEN ET AL. 2002 (MID-LATITUDE)
102 ! = 1, USE MPACE OBSERVATIONS
104 INTEGER, PRIVATE :: INUC
106 ! IBASE = 1, NEGLECT DROPLET ACTIVATION AT LATERAL CLOUD EDGES DUE TO
107 ! UNRESOLVED ENTRAINMENT AND MIXING, ACTIVATE
108 ! AT CLOUD BASE OR IN REGION WITH LITTLE CLOUD WATER USING
109 ! NON-EQULIBRIUM SUPERSATURATION,
110 ! IN CLOUD INTERIOR ACTIVATE USING EQUILIBRIUM SUPERSATURATION
111 ! IBASE = 2, ASSUME DROPLET ACTIVATION AT LATERAL CLOUD EDGES DUE TO
112 ! UNRESOLVED ENTRAINMENT AND MIXING DOMINATES,
113 ! ACTIVATE DROPLETS EVERYWHERE IN THE CLOUD USING NON-EQUILIBRIUM
114 ! SUPERSATURATION, BASED ON THE
115 ! LOCAL SUB-GRID AND/OR GRID-SCALE VERTICAL VELOCITY
116 ! AT THE GRID POINT
118 ! NOTE: ONLY USED FOR PREDICTED DROPLET CONCENTRATION (INUM = 0)
120 INTEGER, PRIVATE :: IBASE
122 ! INCLUDE SUB-GRID VERTICAL VELOCITY IN DROPLET ACTIVATION
123 ! ISUB = 0, INCLUDE SUB-GRID W (RECOMMENDED FOR LOWER RESOLUTION)
124 ! ISUB = 1, EXCLUDE SUB-GRID W, ONLY USE GRID-SCALE W
126 INTEGER, PRIVATE :: ISUB
128 ! SWITCH FOR GRAUPEL/NO GRAUPEL
129 ! IGRAUP = 0, INCLUDE GRAUPEL
130 ! IGRAUP = 1, NO GRAUPEL
132 INTEGER, PRIVATE :: IGRAUP
134 ! HM ADDED NEW OPTION FOR HAIL
135 ! SWITCH FOR HAIL/GRAUPEL
136 ! IHAIL = 0, DENSE PRECIPITATING ICE IS GRAUPEL
137 ! IHAIL = 1, DENSE PRECIPITATING GICE IS HAIL
139 INTEGER, PRIVATE :: IHAIL
141 ! CLOUD MICROPHYSICS CONSTANTS
143 REAL, PRIVATE :: AI,AC,AS,AR,AG ! 'A' PARAMETER IN FALLSPEED-DIAM RELATIONSHIP
144 REAL, PRIVATE :: BI,BC,BS,BR,BG ! 'B' PARAMETER IN FALLSPEED-DIAM RELATIONSHIP
145 ! REAL, PRIVATE :: R ! GAS CONSTANT FOR AIR
146 ! REAL, PRIVATE :: RV ! GAS CONSTANT FOR WATER VAPOR
147 ! REAL, PRIVATE :: CP ! SPECIFIC HEAT AT CONSTANT PRESSURE FOR DRY AIR
148 REAL, PRIVATE :: RHOSU ! STANDARD AIR DENSITY AT 850 MB
149 REAL, PRIVATE :: RHOW ! DENSITY OF LIQUID WATER
150 REAL, PRIVATE :: RHOI ! BULK DENSITY OF CLOUD ICE
151 REAL, PRIVATE :: RHOSN ! BULK DENSITY OF SNOW
152 REAL, PRIVATE :: RHOG ! BULK DENSITY OF GRAUPEL
153 REAL, PRIVATE :: AIMM ! PARAMETER IN BIGG IMMERSION FREEZING
154 REAL, PRIVATE :: BIMM ! PARAMETER IN BIGG IMMERSION FREEZING
155 REAL, PRIVATE :: ECR ! COLLECTION EFFICIENCY BETWEEN DROPLETS/RAIN AND SNOW/RAIN
156 REAL, PRIVATE :: DCS ! THRESHOLD SIZE FOR CLOUD ICE AUTOCONVERSION
157 REAL, PRIVATE :: MI0 ! INITIAL SIZE OF NUCLEATED CRYSTAL
158 REAL, PRIVATE :: MG0 ! MASS OF EMBRYO GRAUPEL
159 REAL, PRIVATE :: F1S ! VENTILATION PARAMETER FOR SNOW
160 REAL, PRIVATE :: F2S ! VENTILATION PARAMETER FOR SNOW
161 REAL, PRIVATE :: F1R ! VENTILATION PARAMETER FOR RAIN
162 REAL, PRIVATE :: F2R ! VENTILATION PARAMETER FOR RAIN
163 ! REAL, PRIVATE :: G ! GRAVITATIONAL ACCELERATION
164 REAL, PRIVATE :: QSMALL ! SMALLEST ALLOWED HYDROMETEOR MIXING RATIO
165 REAL, PRIVATE :: CI,DI,CS,DS,CG,DG ! SIZE DISTRIBUTION PARAMETERS FOR CLOUD ICE, SNOW, GRAUPEL
166 REAL, PRIVATE :: EII ! COLLECTION EFFICIENCY, ICE-ICE COLLISIONS
167 REAL, PRIVATE :: ECI ! COLLECTION EFFICIENCY, ICE-DROPLET COLLISIONS
168 REAL, PRIVATE :: RIN ! RADIUS OF CONTACT NUCLEI (M)
169 ! hm, add for V3.2
170 REAL, PRIVATE :: CPW ! SPECIFIC HEAT OF LIQUID WATER
172 ! CCN SPECTRA FOR IACT = 1
174 REAL, PRIVATE :: C1 ! 'C' IN NCCN = CS^K (CM-3)
175 REAL, PRIVATE :: K1 ! 'K' IN NCCN = CS^K
177 ! AEROSOL PARAMETERS FOR IACT = 2
179 REAL, PRIVATE :: MW ! MOLECULAR WEIGHT WATER (KG/MOL)
180 REAL, PRIVATE :: OSM ! OSMOTIC COEFFICIENT
181 REAL, PRIVATE :: VI ! NUMBER OF ION DISSOCIATED IN SOLUTION
182 REAL, PRIVATE :: EPSM ! AEROSOL SOLUBLE FRACTION
183 REAL, PRIVATE :: RHOA ! AEROSOL BULK DENSITY (KG/M3)
184 REAL, PRIVATE :: MAP ! MOLECULAR WEIGHT AEROSOL (KG/MOL)
185 REAL, PRIVATE :: MA ! MOLECULAR WEIGHT OF 'AIR' (KG/MOL)
186 REAL, PRIVATE :: RR ! UNIVERSAL GAS CONSTANT
187 REAL, PRIVATE :: BACT ! ACTIVATION PARAMETER
188 REAL, PRIVATE :: RM1 ! GEOMETRIC MEAN RADIUS, MODE 1 (M)
189 REAL, PRIVATE :: RM2 ! GEOMETRIC MEAN RADIUS, MODE 2 (M)
190 REAL, PRIVATE :: NANEW1 ! TOTAL AEROSOL CONCENTRATION, MODE 1 (M^-3)
191 REAL, PRIVATE :: NANEW2 ! TOTAL AEROSOL CONCENTRATION, MODE 2 (M^-3)
192 REAL, PRIVATE :: SIG1 ! STANDARD DEVIATION OF AEROSOL S.D., MODE 1
193 REAL, PRIVATE :: SIG2 ! STANDARD DEVIATION OF AEROSOL S.D., MODE 2
194 REAL, PRIVATE :: F11 ! CORRECTION FACTOR FOR ACTIVATION, MODE 1
195 REAL, PRIVATE :: F12 ! CORRECTION FACTOR FOR ACTIVATION, MODE 1
196 REAL, PRIVATE :: F21 ! CORRECTION FACTOR FOR ACTIVATION, MODE 2
197 REAL, PRIVATE :: F22 ! CORRECTION FACTOR FOR ACTIVATION, MODE 2
198 REAL, PRIVATE :: MMULT ! MASS OF SPLINTERED ICE PARTICLE
199 REAL, PRIVATE :: LAMMAXI,LAMMINI,LAMMAXR,LAMMINR,LAMMAXS,LAMMINS,LAMMAXG,LAMMING
201 ! CONSTANTS TO IMPROVE EFFICIENCY
203 REAL, PRIVATE :: CONS1,CONS2,CONS3,CONS4,CONS5,CONS6,CONS7,CONS8,CONS9,CONS10
204 REAL, PRIVATE :: CONS11,CONS12,CONS13,CONS14,CONS15,CONS16,CONS17,CONS18,CONS19,CONS20
205 REAL, PRIVATE :: CONS21,CONS22,CONS23,CONS24,CONS25,CONS26,CONS27,CONS28,CONS29,CONS30
206 REAL, PRIVATE :: CONS31,CONS32,CONS33,CONS34,CONS35,CONS36,CONS37,CONS38,CONS39,CONS40
207 REAL, PRIVATE :: CONS41
209 CONTAINS
211 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
212 SUBROUTINE MORR_TWO_MOMENT_INIT
213 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
214 ! THIS SUBROUTINE INITIALIZES ALL PHYSICAL CONSTANTS AMND PARAMETERS
215 ! NEEDED BY THE MICROPHYSICS SCHEME.
216 ! NEEDS TO BE CALLED AT FIRST TIME STEP, PRIOR TO CALL TO MAIN MICROPHYSICS INTERFACE
217 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
219 IMPLICIT NONE
221 integer n,i
223 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
224 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
226 ! THE FOLLOWING PARAMETERS ARE USER-DEFINED SWITCHES AND NEED TO BE
227 ! SET PRIOR TO CODE COMPILATION
229 ! INUM IS AUTOMATICALLY SET TO 0 FOR WRF-CHEM BELOW,
230 ! ALLOWING PREDICTION OF DROPLET CONCENTRATION
231 ! THUS, THIS PARAMETER SHOULD NOT BE CHANGED HERE
232 ! AND SHOULD BE LEFT TO 1
234 INUM = 1
236 ! SET CONSTANT DROPLET CONCENTRATION (UNITS OF CM-3)
237 ! IF NO COUPLING WITH WRF-CHEM
239 NDCNST = 250.
241 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
242 ! NOTE, THE FOLLOWING OPTIONS RELATED TO DROPLET ACTIVATION
243 ! (IACT, IBASE, ISUB) ARE NOT AVAILABLE IN CURRENT VERSION
244 ! FOR WRF-CHEM, DROPLET ACTIVATION IS PERFORMED
245 ! IN 'MIX_ACTIVATE', NOT IN MICROPHYSICS SCHEME
248 ! IACT = 1, USE POWER-LAW CCN SPECTRA, NCCN = CS^K
249 ! IACT = 2, USE LOGNORMAL AEROSOL SIZE DIST TO DERIVE CCN SPECTRA
251 IACT = 2
253 ! IBASE = 1, NEGLECT DROPLET ACTIVATION AT LATERAL CLOUD EDGES DUE TO
254 ! UNRESOLVED ENTRAINMENT AND MIXING, ACTIVATE
255 ! AT CLOUD BASE OR IN REGION WITH LITTLE CLOUD WATER USING
256 ! NON-EQULIBRIUM SUPERSATURATION ASSUMING NO INITIAL CLOUD WATER,
257 ! IN CLOUD INTERIOR ACTIVATE USING EQUILIBRIUM SUPERSATURATION
258 ! IBASE = 2, ASSUME DROPLET ACTIVATION AT LATERAL CLOUD EDGES DUE TO
259 ! UNRESOLVED ENTRAINMENT AND MIXING DOMINATES,
260 ! ACTIVATE DROPLETS EVERYWHERE IN THE CLOUD USING NON-EQUILIBRIUM
261 ! SUPERSATURATION ASSUMING NO INITIAL CLOUD WATER, BASED ON THE
262 ! LOCAL SUB-GRID AND/OR GRID-SCALE VERTICAL VELOCITY
263 ! AT THE GRID POINT
265 ! NOTE: ONLY USED FOR PREDICTED DROPLET CONCENTRATION (INUM = 0)
267 IBASE = 2
269 ! INCLUDE SUB-GRID VERTICAL VELOCITY (standard deviation of w) IN DROPLET ACTIVATION
270 ! ISUB = 0, INCLUDE SUB-GRID W (RECOMMENDED FOR LOWER RESOLUTION)
271 ! currently, sub-grid w is constant of 0.5 m/s (not coupled with PBL/turbulence scheme)
272 ! ISUB = 1, EXCLUDE SUB-GRID W, ONLY USE GRID-SCALE W
274 ! NOTE: ONLY USED FOR PREDICTED DROPLET CONCENTRATION (INUM = 0)
276 ISUB = 0
277 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
280 ! SWITCH FOR LIQUID-ONLY RUN
281 ! ILIQ = 0, INCLUDE ICE
282 ! ILIQ = 1, LIQUID ONLY, NO ICE
284 ILIQ = 0
286 ! SWITCH FOR ICE NUCLEATION
287 ! INUC = 0, USE FORMULA FROM RASMUSSEN ET AL. 2002 (MID-LATITUDE)
288 ! = 1, USE MPACE OBSERVATIONS (ARCTIC ONLY)
290 INUC = 0
292 ! SWITCH FOR GRAUPEL/HAIL NO GRAUPEL/HAIL
293 ! IGRAUP = 0, INCLUDE GRAUPEL/HAIL
294 ! IGRAUP = 1, NO GRAUPEL/HAIL
296 IGRAUP = 0
298 ! HM ADDED 11/7/07
299 ! SWITCH FOR HAIL/GRAUPEL
300 ! IHAIL = 0, DENSE PRECIPITATING ICE IS GRAUPEL
301 ! IHAIL = 1, DENSE PRECIPITATING ICE IS HAIL
302 ! NOTE ---> RECOMMEND IHAIL = 1 FOR CONTINENTAL DEEP CONVECTION
304 IHAIL = 0
306 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
307 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
308 ! SET PHYSICAL CONSTANTS
310 ! FALLSPEED PARAMETERS (V=AD^B)
311 AI = 700.
312 AC = 3.E7
313 AS = 11.72
314 AR = 841.99667
315 BI = 1.
316 BC = 2.
317 BS = 0.41
318 BR = 0.8
319 IF (IHAIL.EQ.0) THEN
320 AG = 19.3
321 BG = 0.37
322 ELSE ! (MATSUN AND HUGGINS 1980)
323 AG = 114.5
324 BG = 0.5
325 END IF
327 ! CONSTANTS AND PARAMETERS
328 ! R = 287.15
329 ! RV = 461.5
330 ! CP = 1005.
331 RHOSU = 85000./(287.15*273.15)
332 RHOW = 997.
333 RHOI = 500.
334 RHOSN = 100.
335 IF (IHAIL.EQ.0) THEN
336 RHOG = 400.
337 ELSE
338 RHOG = 900.
339 END IF
340 AIMM = 0.66
341 BIMM = 100.
342 ECR = 1.
343 DCS = 125.E-6
344 MI0 = 4./3.*PI*RHOI*(10.E-6)**3
345 MG0 = 1.6E-10
346 F1S = 0.86
347 F2S = 0.28
348 F1R = 0.78
349 ! F2R = 0.32
350 ! fix 053011
351 F2R = 0.308
352 ! G = 9.806
353 QSMALL = 1.E-14
354 EII = 0.1
355 ECI = 0.7
356 ! HM, ADD FOR V3.2
357 CPW = 4218.
359 ! SIZE DISTRIBUTION PARAMETERS
361 CI = RHOI*PI/6.
362 DI = 3.
363 CS = RHOSN*PI/6.
364 DS = 3.
365 CG = RHOG*PI/6.
366 DG = 3.
368 ! RADIUS OF CONTACT NUCLEI
369 RIN = 0.1E-6
371 MMULT = 4./3.*PI*RHOI*(5.E-6)**3
373 ! SIZE LIMITS FOR LAMBDA
375 LAMMAXI = 1./1.E-6
376 LAMMINI = 1./(2.*DCS+100.E-6)
377 LAMMAXR = 1./20.E-6
378 ! LAMMINR = 1./500.E-6
379 LAMMINR = 1./2800.E-6
381 LAMMAXS = 1./10.E-6
382 LAMMINS = 1./2000.E-6
383 LAMMAXG = 1./20.E-6
384 LAMMING = 1./2000.E-6
386 ! CCN SPECTRA FOR IACT = 1
388 ! MARITIME
389 ! MODIFIED FROM RASMUSSEN ET AL. 2002
390 ! NCCN = C*S^K, NCCN IS IN CM-3, S IS SUPERSATURATION RATIO IN %
392 K1 = 0.4
393 C1 = 120.
395 ! CONTINENTAL
397 ! K1 = 0.5
398 ! C1 = 1000.
400 ! AEROSOL ACTIVATION PARAMETERS FOR IACT = 2
401 ! PARAMETERS CURRENTLY SET FOR AMMONIUM SULFATE
403 MW = 0.018
404 OSM = 1.
405 VI = 3.
406 EPSM = 0.7
407 RHOA = 1777.
408 MAP = 0.132
409 MA = 0.0284
410 RR = 8.3187
411 BACT = VI*OSM*EPSM*MW*RHOA/(MAP*RHOW)
413 ! AEROSOL SIZE DISTRIBUTION PARAMETERS CURRENTLY SET FOR MPACE
414 ! (see morrison et al. 2007, JGR)
415 ! MODE 1
417 RM1 = 0.052E-6
418 SIG1 = 2.04
419 NANEW1 = 72.2E6
420 F11 = 0.5*EXP(2.5*(LOG(SIG1))**2)
421 F21 = 1.+0.25*LOG(SIG1)
423 ! MODE 2
425 RM2 = 1.3E-6
426 SIG2 = 2.5
427 NANEW2 = 1.8E6
428 F12 = 0.5*EXP(2.5*(LOG(SIG2))**2)
429 F22 = 1.+0.25*LOG(SIG2)
431 ! CONSTANTS FOR EFFICIENCY
433 CONS1=GAMMA(1.+DS)*CS
434 CONS2=GAMMA(1.+DG)*CG
435 CONS3=GAMMA(4.+BS)/6.
436 CONS4=GAMMA(4.+BR)/6.
437 CONS5=GAMMA(1.+BS)
438 CONS6=GAMMA(1.+BR)
439 CONS7=GAMMA(4.+BG)/6.
440 CONS8=GAMMA(1.+BG)
441 CONS9=GAMMA(5./2.+BR/2.)
442 CONS10=GAMMA(5./2.+BS/2.)
443 CONS11=GAMMA(5./2.+BG/2.)
444 CONS12=GAMMA(1.+DI)*CI
445 CONS13=GAMMA(BS+3.)*PI/4.*ECI
446 CONS14=GAMMA(BG+3.)*PI/4.*ECI
447 CONS15=-1108.*EII*PI**((1.-BS)/3.)*RHOSN**((-2.-BS)/3.)/(4.*720.)
448 CONS16=GAMMA(BI+3.)*PI/4.*ECI
449 CONS17=4.*2.*3.*RHOSU*PI*ECI*ECI*GAMMA(2.*BS+2.)/(8.*(RHOG-RHOSN))
450 CONS18=RHOSN*RHOSN
451 CONS19=RHOW*RHOW
452 CONS20=20.*PI*PI*RHOW*BIMM
453 CONS21=4./(DCS*RHOI)
454 CONS22=PI*RHOI*DCS**3/6.
455 CONS23=PI/4.*EII*GAMMA(BS+3.)
456 CONS24=PI/4.*ECR*GAMMA(BR+3.)
457 CONS25=PI*PI/24.*RHOW*ECR*GAMMA(BR+6.)
458 CONS26=PI/6.*RHOW
459 CONS27=GAMMA(1.+BI)
460 CONS28=GAMMA(4.+BI)/6.
461 CONS29=4./3.*PI*RHOW*(25.E-6)**3
462 CONS30=4./3.*PI*RHOW
463 CONS31=PI*PI*ECR*RHOSN
464 CONS32=PI/2.*ECR
465 CONS33=PI*PI*ECR*RHOG
466 CONS34=5./2.+BR/2.
467 CONS35=5./2.+BS/2.
468 CONS36=5./2.+BG/2.
469 CONS37=4.*PI*1.38E-23/(6.*PI*RIN)
470 CONS38=PI*PI/3.*RHOW
471 CONS39=PI*PI/36.*RHOW*BIMM
472 CONS40=PI/6.*BIMM
473 CONS41=PI*PI*ECR*RHOW
475 END SUBROUTINE MORR_TWO_MOMENT_INIT
477 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
478 ! THIS SUBROUTINE IS MAIN INTERFACE WITH THE TWO-MOMENT MICROPHYSICS SCHEME
479 ! THIS INTERFACE TAKES IN 3D VARIABLES FROM DRIVER MODEL, CONVERTS TO 1D FOR
480 ! CALL TO THE MAIN MICROPHYSICS SUBROUTINE (SUBROUTINE MORR_TWO_MOMENT_MICRO)
481 ! WHICH OPERATES ON 1D VERTICAL COLUMNS.
482 ! 1D VARIABLES FROM THE MAIN MICROPHYSICS SUBROUTINE ARE THEN REASSIGNED BACK TO 3D FOR OUTPUT
483 ! BACK TO DRIVER MODEL USING THIS INTERFACE.
484 ! MICROPHYSICS TENDENCIES ARE ADDED TO VARIABLES HERE BEFORE BEING PASSED BACK TO DRIVER MODEL.
486 ! THIS CODE WAS WRITTEN BY HUGH MORRISON (NCAR) AND SLAVA TATARSKII (GEORGIA TECH).
488 ! FOR QUESTIONS, CONTACT: HUGH MORRISON, E-MAIL: MORRISON@UCAR.EDU, PHONE:303-497-8916
490 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
492 SUBROUTINE MP_MORR_TWO_MOMENT(ITIMESTEP, &
493 TH, QV, QC, QR, QI, QS, QG, NI, NS, NR, NG, &
494 RHO, PII, P, DT_IN, DZ, HT, W, &
495 RAINNC, RAINNCV, SR, &
496 qrcuten, qscuten, qicuten, mu & ! hm added
497 ,F_QNDROP, qndrop & ! hm added, wrf-chem
498 ,IDS,IDE, JDS,JDE, KDS,KDE & ! domain dims
499 ,IMS,IME, JMS,JME, KMS,KME & ! memory dims
500 ,ITS,ITE, JTS,JTE, KTS,KTE & ! tile dims )
501 !jdf ,C2PREC3D,CSED3D,ISED3D,SSED3D,GSED3D,RSED3D & ! HM ADD, WRF-CHEM
502 ,QLSINK,PRECR,PRECI,PRECS,PRECG & ! HM ADD, WRF-CHEM
503 )
505 ! QV - water vapor mixing ratio (kg/kg)
506 ! QC - cloud water mixing ratio (kg/kg)
507 ! QR - rain water mixing ratio (kg/kg)
508 ! QI - cloud ice mixing ratio (kg/kg)
509 ! QS - snow mixing ratio (kg/kg)
510 ! QG - graupel mixing ratio (KG/KG)
511 ! NI - cloud ice number concentration (1/kg)
512 ! NS - Snow Number concentration (1/kg)
513 ! NR - Rain Number concentration (1/kg)
514 ! NG - Graupel number concentration (1/kg)
515 ! NOTE: RHO AND HT NOT USED BY THIS SCHEME AND DO NOT NEED TO BE PASSED INTO SCHEME!!!!
516 ! P - AIR PRESSURE (PA)
517 ! W - VERTICAL AIR VELOCITY (M/S)
518 ! TH - POTENTIAL TEMPERATURE (K)
519 ! PII - exner function - used to convert potential temp to temp
520 ! DZ - difference in height over interface (m)
521 ! DT_IN - model time step (sec)
522 ! ITIMESTEP - time step counter
523 ! RAINNC - accumulated grid-scale precipitation (mm)
524 ! RAINNCV - one time step grid scale precipitation (mm/time step)
525 ! SR - one time step mass ratio of snow to total precip
526 ! qrcuten, rain tendency from parameterized cumulus convection
527 ! qscuten, snow tendency from parameterized cumulus convection
528 ! qicuten, cloud ice tendency from parameterized cumulus convection
530 ! variables below currently not in use, not coupled to PBL or radiation codes
531 ! TKE - turbulence kinetic energy (m^2 s-2), NEEDED FOR DROPLET ACTIVATION (SEE CODE BELOW)
532 ! NCTEND - droplet concentration tendency from pbl (kg-1 s-1)
533 ! NCTEND - CLOUD ICE concentration tendency from pbl (kg-1 s-1)
534 ! KZH - heat eddy diffusion coefficient from YSU scheme (M^2 S-1), NEEDED FOR DROPLET ACTIVATION (SEE CODE BELOW)
535 ! EFFCS - CLOUD DROPLET EFFECTIVE RADIUS OUTPUT TO RADIATION CODE (micron)
536 ! EFFIS - CLOUD DROPLET EFFECTIVE RADIUS OUTPUT TO RADIATION CODE (micron)
537 ! HM, ADDED FOR WRF-CHEM COUPLING
538 ! QLSINK - TENDENCY OF CLOUD WATER TO RAIN, SNOW, GRAUPEL (KG/KG/S)
539 ! CSED,ISED,SSED,GSED,RSED - SEDIMENTATION FLUXES (KG/M^2/S) FOR CLOUD WATER, ICE, SNOW, GRAUPEL, RAIN
540 ! PRECI,PRECS,PRECG,PRECR - SEDIMENTATION FLUXES (KG/M^2/S) FOR ICE, SNOW, GRAUPEL, RAIN
541 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
543 ! reflectivity currently not included!!!!
544 ! REFL_10CM - CALCULATED RADAR REFLECTIVITY AT 10 CM (DBZ)
545 !................................
546 ! GRID_CLOCK, GRID_ALARMS - parameters to limit radar reflectivity calculation only when needed
547 ! otherwise radar reflectivity calculation every time step is too slow
548 ! only needed for coupling with WRF, see code below for details
549 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
551 ! EFFC - DROPLET EFFECTIVE RADIUS (MICRON)
552 ! EFFR - RAIN EFFECTIVE RADIUS (MICRON)
553 ! EFFS - SNOW EFFECTIVE RADIUS (MICRON)
554 ! EFFI - CLOUD ICE EFFECTIVE RADIUS (MICRON)
556 ! ADDITIONAL OUTPUT FROM MICRO - SEDIMENTATION TENDENCIES, NEEDED FOR LIQUID-ICE STATIC ENERGY
558 ! QGSTEN - GRAUPEL SEDIMENTATION TEND (KG/KG/S)
559 ! QRSTEN - RAIN SEDIMENTATION TEND (KG/KG/S)
560 ! QISTEN - CLOUD ICE SEDIMENTATION TEND (KG/KG/S)
561 ! QNISTEN - SNOW SEDIMENTATION TEND (KG/KG/S)
562 ! QCSTEN - CLOUD WATER SEDIMENTATION TEND (KG/KG/S)
564 ! WVAR - STANDARD DEVIATION OF SUB-GRID VERTICAL VELOCITY (M/S)
566 IMPLICIT NONE
568 INTEGER, INTENT(IN ) :: ids, ide, jds, jde, kds, kde , &
569 ims, ime, jms, jme, kms, kme , &
570 its, ite, jts, jte, kts, kte
571 ! Temporary changed from INOUT to IN
573 REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT):: &
574 qv, qc, qr, qi, qs, qg, ni, ns, nr, TH, NG
575 !jdf qndrop ! hm added, wrf-chem
576 REAL, DIMENSION(ims:ime, kms:kme, jms:jme), optional,INTENT(INOUT):: qndrop
577 !jdf REAL, DIMENSION(ims:ime, kms:kme, jms:jme),INTENT(INOUT):: CSED3D, &
578 REAL, DIMENSION(ims:ime, kms:kme, jms:jme), optional,INTENT(INOUT):: QLSINK, &
579 PRECI,PRECS,PRECG,PRECR ! HM, WRF-CHEM
580 !, effcs, effis
582 REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(IN):: &
583 pii, p, dz, rho, w !, tke, nctend, nitend,kzh
584 REAL, INTENT(IN):: dt_in
585 INTEGER, INTENT(IN):: ITIMESTEP
587 REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT):: &
588 RAINNC, RAINNCV, SR
590 ! REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT):: & ! GT
591 ! refl_10cm
593 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN) :: ht
595 ! TYPE (WRFU_Clock):: grid_clock ! GT
596 ! TYPE (WRFU_Alarm), POINTER:: grid_alarms(:) ! GT
598 ! LOCAL VARIABLES
600 REAL, DIMENSION(its:ite, kts:kte, jts:jte):: &
601 effi, effs, effr, EFFG
603 REAL, DIMENSION(its:ite, kts:kte, jts:jte):: &
604 T, WVAR, EFFC
606 REAL, DIMENSION(kts:kte) :: &
607 QC_TEND1D, QI_TEND1D, QNI_TEND1D, QR_TEND1D, &
608 NI_TEND1D, NS_TEND1D, NR_TEND1D, &
609 QC1D, QI1D, QR1D,NI1D, NS1D, NR1D, QS1D, &
610 T_TEND1D,QV_TEND1D, T1D, QV1D, P1D, W1D, WVAR1D, &
611 EFFC1D, EFFI1D, EFFS1D, EFFR1D,DZ1D, &
612 ! HM ADD GRAUPEL
613 QG_TEND1D, NG_TEND1D, QG1D, NG1D, EFFG1D, &
615 ! ADD SEDIMENTATION TENDENCIES (UNITS OF KG/KG/S)
616 QGSTEN,QRSTEN, QISTEN, QNISTEN, QCSTEN, &
617 ! ADD CUMULUS TENDENCIES
618 QRCU1D, QSCU1D, QICU1D
620 ! add cumulus tendencies
622 REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(IN):: &
623 qrcuten, qscuten, qicuten
624 REAL, DIMENSION(ims:ime, jms:jme), INTENT(IN):: &
625 mu
627 LOGICAL, INTENT(IN), OPTIONAL :: F_QNDROP ! wrf-chem
628 LOGICAL :: flag_qndrop ! wrf-chem
629 integer :: iinum ! wrf-chem
631 ! wrf-chem
632 REAL, DIMENSION(kts:kte) :: nc1d, nc_tend1d,C2PREC,CSED,ISED,SSED,GSED,RSED
633 ! HM add reflectivity
634 ! dbz
636 REAL PRECPRT1D, SNOWRT1D
638 INTEGER I,K,J
640 REAL DT
642 ! LOGICAL:: dBZ_tstep ! GT
643 ! below for wrf-chem
644 flag_qndrop = .false.
645 IF ( PRESENT ( f_qndrop ) ) flag_qndrop = f_qndrop
646 !!!!!!!!!!!!!!!!!!!!!!
648 ! Initialize tendencies (all set to 0) and transfer
649 ! array to local variables
650 DT = DT_IN
652 DO I=ITS,ITE
653 DO J=JTS,JTE
654 DO K=KTS,KTE
655 T(I,K,J) = TH(i,k,j)*PII(i,k,j)
657 ! NOTE: WVAR NOT CURRENTLY USED IN CODE !!!!!!!!!!
658 ! currently assign wvar to 0.5 m/s (not coupled with PBL scheme)
660 WVAR(I,K,J) = 0.5
662 ! currently mixing of number concentrations also is neglected (not coupled with PBL schemes)
664 END DO
665 END DO
666 END DO
668 do i=its,ite ! i loop (east-west)
669 do j=jts,jte ! j loop (north-south)
670 !
671 ! Transfer 3D arrays into 1D for microphysical calculations
672 !
674 ! hm , initialize 1d tendency arrays to zero
676 do k=kts,kte ! k loop (vertical)
678 QC_TEND1D(k) = 0.
679 QI_TEND1D(k) = 0.
680 QNI_TEND1D(k) = 0.
681 QR_TEND1D(k) = 0.
682 NI_TEND1D(k) = 0.
683 NS_TEND1D(k) = 0.
684 NR_TEND1D(k) = 0.
685 T_TEND1D(k) = 0.
686 QV_TEND1D(k) = 0.
687 nc_tend1d(k) = 0. ! wrf-chem
689 QC1D(k) = QC(i,k,j)
690 QI1D(k) = QI(i,k,j)
691 QS1D(k) = QS(i,k,j)
692 QR1D(k) = QR(i,k,j)
694 NI1D(k) = NI(i,k,j)
696 NS1D(k) = NS(i,k,j)
697 NR1D(k) = NR(i,k,j)
698 ! HM ADD GRAUPEL
699 QG1D(K) = QG(I,K,j)
700 NG1D(K) = NG(I,K,j)
701 QG_TEND1D(K) = 0.
702 NG_TEND1D(K) = 0.
704 T1D(k) = T(i,k,j)
705 QV1D(k) = QV(i,k,j)
706 P1D(k) = P(i,k,j)
707 DZ1D(k) = DZ(i,k,j)
708 W1D(k) = W(i,k,j)
709 WVAR1D(k) = WVAR(i,k,j)
710 ! add cumulus tendencies, decouple from mu
711 qrcu1d(k) = qrcuten(i,k,j)/mu(i,j)
712 qscu1d(k) = qscuten(i,k,j)/mu(i,j)
713 qicu1d(k) = qicuten(i,k,j)/mu(i,j)
714 end do !jdf added this
715 ! below for wrf-chem
716 IF (flag_qndrop .AND. PRESENT( qndrop )) THEN
717 iact = 3
718 DO k = kts, kte
719 nc1d(k)=qndrop(i,k,j)
720 iinum=0
721 ENDDO
722 ELSE
723 DO k = kts, kte
724 nc1d(k)=0. ! temporary placeholder, set to constant in microphysics subroutine
725 iinum=1
726 ENDDO
727 ENDIF
729 !jdf end do
731 call MORR_TWO_MOMENT_MICRO(QC_TEND1D, QI_TEND1D, QNI_TEND1D, QR_TEND1D, &
732 NI_TEND1D, NS_TEND1D, NR_TEND1D, &
733 QC1D, QI1D, QS1D, QR1D,NI1D, NS1D, NR1D, &
734 T_TEND1D,QV_TEND1D, T1D, QV1D, P1D, DZ1D, W1D, WVAR1D, &
735 PRECPRT1D,SNOWRT1D, &
736 EFFC1D,EFFI1D,EFFS1D,EFFR1D,DT, &
737 IMS,IME, JMS,JME, KMS,KME, &
738 ITS,ITE, JTS,JTE, KTS,KTE, & ! HM ADD GRAUPEL
739 QG_TEND1D,NG_TEND1D,QG1D,NG1D,EFFG1D, &
740 qrcu1d, qscu1d, qicu1d, &
741 ! ADD SEDIMENTATION TENDENCIES
742 QGSTEN,QRSTEN,QISTEN,QNISTEN,QCSTEN, &
743 nc1d, nc_tend1d, iinum, C2PREC,CSED,ISED,SSED,GSED,RSED & !wrf-chem
744 )
746 !
747 ! Transfer 1D arrays back into 3D arrays
748 !
749 do k=kts,kte
751 ! hm, add tendencies to update global variables
752 ! HM, TENDENCIES FOR Q AND N NOW ADDED IN M2005MICRO, SO WE
753 ! ONLY NEED TO TRANSFER 1D VARIABLES BACK TO 3D
755 QC(i,k,j) = QC1D(k)
756 QI(i,k,j) = QI1D(k)
757 QS(i,k,j) = QS1D(k)
758 QR(i,k,j) = QR1D(k)
759 NI(i,k,j) = NI1D(k)
760 NS(i,k,j) = NS1D(k)
761 NR(i,k,j) = NR1D(k)
762 QG(I,K,j) = QG1D(K)
763 NG(I,K,j) = NG1D(K)
765 T(i,k,j) = T1D(k)
766 TH(I,K,J) = T(i,k,j)/PII(i,k,j) ! CONVERT TEMP BACK TO POTENTIAL TEMP
767 QV(i,k,j) = QV1D(k)
769 EFFC(i,k,j) = EFFC1D(k)
770 EFFI(i,k,j) = EFFI1D(k)
771 EFFS(i,k,j) = EFFS1D(k)
772 EFFR(i,k,j) = EFFR1D(k)
773 EFFG(I,K,j) = EFFG1D(K)
775 ! wrf-chem
776 IF (flag_qndrop .AND. PRESENT( qndrop )) THEN
777 qndrop(i,k,j) = nc1d(k)
778 !jdf CSED3D(I,K,J) = CSED(K)
779 END IF
780 IF ( PRESENT( QLSINK ) ) THEN
781 if(qc(i,k,j)>1.e-10) then
782 QLSINK(I,K,J) = C2PREC(K)/QC(I,K,J)
783 else
784 QLSINK(I,K,J) = 0.0
785 endif
786 END IF
787 IF ( PRESENT( PRECR ) ) PRECR(I,K,J) = RSED(K)
788 IF ( PRESENT( PRECI ) ) PRECI(I,K,J) = ISED(K)
789 IF ( PRESENT( PRECS ) ) PRECS(I,K,J) = SSED(K)
790 IF ( PRESENT( PRECG ) ) PRECG(I,K,J) = GSED(K)
791 ! EFFECTIVE RADIUS FOR RADIATION CODE (currently not coupled)
792 ! HM, ADD LIMIT TO PREVENT BLOWING UP OPTICAL PROPERTIES, 8/18/07
793 ! EFFCS(I,K,J) = MIN(EFFC(I,K,J),50.)
794 ! EFFCS(I,K,J) = MAX(EFFCS(I,K,J),1.)
795 ! EFFIS(I,K,J) = MIN(EFFI(I,K,J),130.)
796 ! EFFIS(I,K,J) = MAX(EFFIS(I,K,J),13.)
798 end do
800 ! hm modified so that m2005 precip variables correctly match wrf precip variables
801 RAINNC(i,j) = RAINNC(I,J)+PRECPRT1D
802 RAINNCV(i,j) = PRECPRT1D
803 SR(i,j) = SNOWRT1D/(PRECPRT1D+1.E-12)
805 end do
806 end do
808 END SUBROUTINE MP_MORR_TWO_MOMENT
810 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
811 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
812 SUBROUTINE MORR_TWO_MOMENT_MICRO(QC3DTEN,QI3DTEN,QNI3DTEN,QR3DTEN, &
813 NI3DTEN,NS3DTEN,NR3DTEN,QC3D,QI3D,QNI3D,QR3D,NI3D,NS3D,NR3D, &
814 T3DTEN,QV3DTEN,T3D,QV3D,PRES,DZQ,W3D,WVAR,PRECRT,SNOWRT, &
815 EFFC,EFFI,EFFS,EFFR,DT, &
816 IMS,IME, JMS,JME, KMS,KME, &
817 ITS,ITE, JTS,JTE, KTS,KTE, & ! ADD GRAUPEL
818 QG3DTEN,NG3DTEN,QG3D,NG3D,EFFG,qrcu1d,qscu1d, qicu1d, &
819 QGSTEN,QRSTEN,QISTEN,QNISTEN,QCSTEN, &
820 nc3d,nc3dten,iinum, & ! wrf-chem
821 c2prec,CSED,ISED,SSED,GSED,RSED & ! hm added, wrf-chem
822 )
824 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
825 ! THIS PROGRAM IS THE MAIN TWO-MOMENT MICROPHYSICS SUBROUTINE DESCRIBED BY
826 ! MORRISON ET AL. 2005 JAS; MORRISON AND PINTO 2005 JAS.
827 ! ADDITIONAL CHANGES ARE DESCRIBED IN DETAIL BY MORRISON, THOMPSON, TATARSKII (MWR, SUBMITTED)
829 ! THIS SCHEME IS A BULK DOUBLE-MOMENT SCHEME THAT PREDICTS MIXING
830 ! RATIOS AND NUMBER CONCENTRATIONS OF FIVE HYDROMETEOR SPECIES:
831 ! CLOUD DROPLETS, CLOUD (SMALL) ICE, RAIN, SNOW, AND GRAUPEL.
833 ! CODE STRUCTURE: MAIN SUBROUTINE IS 'MORR_TWO_MOMENT'. ALSO INCLUDED IN THIS FILE IS
834 ! 'FUNCTION POLYSVP', 'FUNCTION DERF1', AND
835 ! 'FUNCTION GAMMA'.
837 ! NOTE: THIS SUBROUTINE USES 1D ARRAY IN VERTICAL (COLUMN), EVEN THOUGH VARIABLES ARE CALLED '3D'......
839 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
841 ! DECLARATIONS
843 IMPLICIT NONE
845 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
846 ! THESE VARIABLES BELOW MUST BE LINKED WITH THE MAIN MODEL.
847 ! DEFINE ARRAY SIZES
849 ! INPUT NUMBER OF GRID CELLS
851 ! INPUT/OUTPUT PARAMETERS ! DESCRIPTION (UNITS)
852 INTEGER, INTENT( IN) :: IMS,IME, JMS,JME, KMS,KME, &
853 ITS,ITE, JTS,JTE, KTS,KTE
855 REAL, DIMENSION(KTS:KTE) :: QC3DTEN ! CLOUD WATER MIXING RATIO TENDENCY (KG/KG/S)
856 REAL, DIMENSION(KTS:KTE) :: QI3DTEN ! CLOUD ICE MIXING RATIO TENDENCY (KG/KG/S)
857 REAL, DIMENSION(KTS:KTE) :: QNI3DTEN ! SNOW MIXING RATIO TENDENCY (KG/KG/S)
858 REAL, DIMENSION(KTS:KTE) :: QR3DTEN ! RAIN MIXING RATIO TENDENCY (KG/KG/S)
859 REAL, DIMENSION(KTS:KTE) :: NI3DTEN ! CLOUD ICE NUMBER CONCENTRATION (1/KG/S)
860 REAL, DIMENSION(KTS:KTE) :: NS3DTEN ! SNOW NUMBER CONCENTRATION (1/KG/S)
861 REAL, DIMENSION(KTS:KTE) :: NR3DTEN ! RAIN NUMBER CONCENTRATION (1/KG/S)
862 REAL, DIMENSION(KTS:KTE) :: QC3D ! CLOUD WATER MIXING RATIO (KG/KG)
863 REAL, DIMENSION(KTS:KTE) :: QI3D ! CLOUD ICE MIXING RATIO (KG/KG)
864 REAL, DIMENSION(KTS:KTE) :: QNI3D ! SNOW MIXING RATIO (KG/KG)
865 REAL, DIMENSION(KTS:KTE) :: QR3D ! RAIN MIXING RATIO (KG/KG)
866 REAL, DIMENSION(KTS:KTE) :: NI3D ! CLOUD ICE NUMBER CONCENTRATION (1/KG)
867 REAL, DIMENSION(KTS:KTE) :: NS3D ! SNOW NUMBER CONCENTRATION (1/KG)
868 REAL, DIMENSION(KTS:KTE) :: NR3D ! RAIN NUMBER CONCENTRATION (1/KG)
869 REAL, DIMENSION(KTS:KTE) :: T3DTEN ! TEMPERATURE TENDENCY (K/S)
870 REAL, DIMENSION(KTS:KTE) :: QV3DTEN ! WATER VAPOR MIXING RATIO TENDENCY (KG/KG/S)
871 REAL, DIMENSION(KTS:KTE) :: T3D ! TEMPERATURE (K)
872 REAL, DIMENSION(KTS:KTE) :: QV3D ! WATER VAPOR MIXING RATIO (KG/KG)
873 REAL, DIMENSION(KTS:KTE) :: PRES ! ATMOSPHERIC PRESSURE (PA)
874 REAL, DIMENSION(KTS:KTE) :: DZQ ! DIFFERENCE IN HEIGHT ACROSS LEVEL (m)
875 REAL, DIMENSION(KTS:KTE) :: W3D ! GRID-SCALE VERTICAL VELOCITY (M/S)
876 REAL, DIMENSION(KTS:KTE) :: WVAR ! SUB-GRID VERTICAL VELOCITY (M/S)
877 ! below for wrf-chem
878 REAL, DIMENSION(KTS:KTE) :: nc3d
879 REAL, DIMENSION(KTS:KTE) :: nc3dten
880 integer, intent(in) :: iinum
882 ! HM ADDED GRAUPEL VARIABLES
883 REAL, DIMENSION(KTS:KTE) :: QG3DTEN ! GRAUPEL MIX RATIO TENDENCY (KG/KG/S)
884 REAL, DIMENSION(KTS:KTE) :: NG3DTEN ! GRAUPEL NUMB CONC TENDENCY (1/KG/S)
885 REAL, DIMENSION(KTS:KTE) :: QG3D ! GRAUPEL MIX RATIO (KG/KG)
886 REAL, DIMENSION(KTS:KTE) :: NG3D ! GRAUPEL NUMBER CONC (1/KG)
888 ! HM, ADD 1/16/07, SEDIMENTATION TENDENCIES FOR MIXING RATIO
890 REAL, DIMENSION(KTS:KTE) :: QGSTEN ! GRAUPEL SED TEND (KG/KG/S)
891 REAL, DIMENSION(KTS:KTE) :: QRSTEN ! RAIN SED TEND (KG/KG/S)
892 REAL, DIMENSION(KTS:KTE) :: QISTEN ! CLOUD ICE SED TEND (KG/KG/S)
893 REAL, DIMENSION(KTS:KTE) :: QNISTEN ! SNOW SED TEND (KG/KG/S)
894 REAL, DIMENSION(KTS:KTE) :: QCSTEN ! CLOUD WAT SED TEND (KG/KG/S)
896 ! hm add cumulus tendencies for precip
897 REAL, DIMENSION(KTS:KTE) :: qrcu1d
898 REAL, DIMENSION(KTS:KTE) :: qscu1d
899 REAL, DIMENSION(KTS:KTE) :: qicu1d
901 ! OUTPUT VARIABLES
903 REAL PRECRT ! TOTAL PRECIP PER TIME STEP (mm)
904 REAL SNOWRT ! SNOW PER TIME STEP (mm)
906 REAL, DIMENSION(KTS:KTE) :: EFFC ! DROPLET EFFECTIVE RADIUS (MICRON)
907 REAL, DIMENSION(KTS:KTE) :: EFFI ! CLOUD ICE EFFECTIVE RADIUS (MICRON)
908 REAL, DIMENSION(KTS:KTE) :: EFFS ! SNOW EFFECTIVE RADIUS (MICRON)
909 REAL, DIMENSION(KTS:KTE) :: EFFR ! RAIN EFFECTIVE RADIUS (MICRON)
910 REAL, DIMENSION(KTS:KTE) :: EFFG ! GRAUPEL EFFECTIVE RADIUS (MICRON)
912 ! MODEL INPUT PARAMETERS (FORMERLY IN COMMON BLOCKS)
914 REAL DT ! MODEL TIME STEP (SEC)
916 !.....................................................................................................
917 ! LOCAL VARIABLES: ALL PARAMETERS BELOW ARE LOCAL TO SCHEME AND DON'T NEED TO COMMUNICATE WITH THE
918 ! REST OF THE MODEL.
920 ! SIZE PARAMETER VARIABLES
922 REAL, DIMENSION(KTS:KTE) :: LAMC ! SLOPE PARAMETER FOR DROPLETS (M-1)
923 REAL, DIMENSION(KTS:KTE) :: LAMI ! SLOPE PARAMETER FOR CLOUD ICE (M-1)
924 REAL, DIMENSION(KTS:KTE) :: LAMS ! SLOPE PARAMETER FOR SNOW (M-1)
925 REAL, DIMENSION(KTS:KTE) :: LAMR ! SLOPE PARAMETER FOR RAIN (M-1)
926 REAL, DIMENSION(KTS:KTE) :: LAMG ! SLOPE PARAMETER FOR GRAUPEL (M-1)
927 REAL, DIMENSION(KTS:KTE) :: CDIST1 ! PSD PARAMETER FOR DROPLETS
928 REAL, DIMENSION(KTS:KTE) :: N0I ! INTERCEPT PARAMETER FOR CLOUD ICE (KG-1 M-1)
929 REAL, DIMENSION(KTS:KTE) :: N0S ! INTERCEPT PARAMETER FOR SNOW (KG-1 M-1)
930 REAL, DIMENSION(KTS:KTE) :: N0RR ! INTERCEPT PARAMETER FOR RAIN (KG-1 M-1)
931 REAL, DIMENSION(KTS:KTE) :: N0G ! INTERCEPT PARAMETER FOR GRAUPEL (KG-1 M-1)
932 REAL, DIMENSION(KTS:KTE) :: PGAM ! SPECTRAL SHAPE PARAMETER FOR DROPLETS
934 ! MICROPHYSICAL PROCESSES
936 REAL, DIMENSION(KTS:KTE) :: NSUBC ! LOSS OF NC DURING EVAP
937 REAL, DIMENSION(KTS:KTE) :: NSUBI ! LOSS OF NI DURING SUB.
938 REAL, DIMENSION(KTS:KTE) :: NSUBS ! LOSS OF NS DURING SUB.
939 REAL, DIMENSION(KTS:KTE) :: NSUBR ! LOSS OF NR DURING EVAP
940 REAL, DIMENSION(KTS:KTE) :: PRD ! DEP CLOUD ICE
941 REAL, DIMENSION(KTS:KTE) :: PRE ! EVAP OF RAIN
942 REAL, DIMENSION(KTS:KTE) :: PRDS ! DEP SNOW
943 REAL, DIMENSION(KTS:KTE) :: NNUCCC ! CHANGE N DUE TO CONTACT FREEZ DROPLETS
944 REAL, DIMENSION(KTS:KTE) :: MNUCCC ! CHANGE Q DUE TO CONTACT FREEZ DROPLETS
945 REAL, DIMENSION(KTS:KTE) :: PRA ! ACCRETION DROPLETS BY RAIN
946 REAL, DIMENSION(KTS:KTE) :: PRC ! AUTOCONVERSION DROPLETS
947 REAL, DIMENSION(KTS:KTE) :: PCC ! COND/EVAP DROPLETS
948 REAL, DIMENSION(KTS:KTE) :: NNUCCD ! CHANGE N FREEZING AEROSOL (PRIM ICE NUCLEATION)
949 REAL, DIMENSION(KTS:KTE) :: MNUCCD ! CHANGE Q FREEZING AEROSOL (PRIM ICE NUCLEATION)
950 REAL, DIMENSION(KTS:KTE) :: MNUCCR ! CHANGE Q DUE TO CONTACT FREEZ RAIN
951 REAL, DIMENSION(KTS:KTE) :: NNUCCR ! CHANGE N DUE TO CONTACT FREEZ RAIN
952 REAL, DIMENSION(KTS:KTE) :: NPRA ! CHANGE IN N DUE TO DROPLET ACC BY RAIN
953 REAL, DIMENSION(KTS:KTE) :: NRAGG ! SELF-COLLECTION OF RAIN
954 REAL, DIMENSION(KTS:KTE) :: NSAGG ! SELF-COLLECTION OF SNOW
955 REAL, DIMENSION(KTS:KTE) :: NPRC ! CHANGE NC AUTOCONVERSION DROPLETS
956 REAL, DIMENSION(KTS:KTE) :: NPRC1 ! CHANGE NR AUTOCONVERSION DROPLETS
957 REAL, DIMENSION(KTS:KTE) :: PRAI ! CHANGE Q AUTOCONVERSION CLOUD ICE
958 REAL, DIMENSION(KTS:KTE) :: PRCI ! CHANGE Q ACCRETION CLOUD ICE BY SNOW
959 REAL, DIMENSION(KTS:KTE) :: PSACWS ! CHANGE Q DROPLET ACCRETION BY SNOW
960 REAL, DIMENSION(KTS:KTE) :: NPSACWS ! CHANGE N DROPLET ACCRETION BY SNOW
961 REAL, DIMENSION(KTS:KTE) :: PSACWI ! CHANGE Q DROPLET ACCRETION BY CLOUD ICE
962 REAL, DIMENSION(KTS:KTE) :: NPSACWI ! CHANGE N DROPLET ACCRETION BY CLOUD ICE
963 REAL, DIMENSION(KTS:KTE) :: NPRCI ! CHANGE N AUTOCONVERSION CLOUD ICE BY SNOW
964 REAL, DIMENSION(KTS:KTE) :: NPRAI ! CHANGE N ACCRETION CLOUD ICE
965 REAL, DIMENSION(KTS:KTE) :: NMULTS ! ICE MULT DUE TO RIMING DROPLETS BY SNOW
966 REAL, DIMENSION(KTS:KTE) :: NMULTR ! ICE MULT DUE TO RIMING RAIN BY SNOW
967 REAL, DIMENSION(KTS:KTE) :: QMULTS ! CHANGE Q DUE TO ICE MULT DROPLETS/SNOW
968 REAL, DIMENSION(KTS:KTE) :: QMULTR ! CHANGE Q DUE TO ICE RAIN/SNOW
969 REAL, DIMENSION(KTS:KTE) :: PRACS ! CHANGE Q RAIN-SNOW COLLECTION
970 REAL, DIMENSION(KTS:KTE) :: NPRACS ! CHANGE N RAIN-SNOW COLLECTION
971 REAL, DIMENSION(KTS:KTE) :: PCCN ! CHANGE Q DROPLET ACTIVATION
972 REAL, DIMENSION(KTS:KTE) :: PSMLT ! CHANGE Q MELTING SNOW TO RAIN
973 REAL, DIMENSION(KTS:KTE) :: EVPMS ! CHNAGE Q MELTING SNOW EVAPORATING
974 REAL, DIMENSION(KTS:KTE) :: NSMLTS ! CHANGE N MELTING SNOW
975 REAL, DIMENSION(KTS:KTE) :: NSMLTR ! CHANGE N MELTING SNOW TO RAIN
976 ! HM ADDED 12/13/06
977 REAL, DIMENSION(KTS:KTE) :: PIACR ! CHANGE QR, ICE-RAIN COLLECTION
978 REAL, DIMENSION(KTS:KTE) :: NIACR ! CHANGE N, ICE-RAIN COLLECTION
979 REAL, DIMENSION(KTS:KTE) :: PRACI ! CHANGE QI, ICE-RAIN COLLECTION
980 REAL, DIMENSION(KTS:KTE) :: PIACRS ! CHANGE QR, ICE RAIN COLLISION, ADDED TO SNOW
981 REAL, DIMENSION(KTS:KTE) :: NIACRS ! CHANGE N, ICE RAIN COLLISION, ADDED TO SNOW
982 REAL, DIMENSION(KTS:KTE) :: PRACIS ! CHANGE QI, ICE RAIN COLLISION, ADDED TO SNOW
983 REAL, DIMENSION(KTS:KTE) :: EPRD ! SUBLIMATION CLOUD ICE
984 REAL, DIMENSION(KTS:KTE) :: EPRDS ! SUBLIMATION SNOW
985 ! HM ADDED GRAUPEL PROCESSES
986 REAL, DIMENSION(KTS:KTE) :: PRACG ! CHANGE IN Q COLLECTION RAIN BY GRAUPEL
987 REAL, DIMENSION(KTS:KTE) :: PSACWG ! CHANGE IN Q COLLECTION DROPLETS BY GRAUPEL
988 REAL, DIMENSION(KTS:KTE) :: PGSACW ! CONVERSION Q TO GRAUPEL DUE TO COLLECTION DROPLETS BY SNOW
989 REAL, DIMENSION(KTS:KTE) :: PGRACS ! CONVERSION Q TO GRAUPEL DUE TO COLLECTION RAIN BY SNOW
990 REAL, DIMENSION(KTS:KTE) :: PRDG ! DEP OF GRAUPEL
991 REAL, DIMENSION(KTS:KTE) :: EPRDG ! SUB OF GRAUPEL
992 REAL, DIMENSION(KTS:KTE) :: EVPMG ! CHANGE Q MELTING OF GRAUPEL AND EVAPORATION
993 REAL, DIMENSION(KTS:KTE) :: PGMLT ! CHANGE Q MELTING OF GRAUPEL
994 REAL, DIMENSION(KTS:KTE) :: NPRACG ! CHANGE N COLLECTION RAIN BY GRAUPEL
995 REAL, DIMENSION(KTS:KTE) :: NPSACWG ! CHANGE N COLLECTION DROPLETS BY GRAUPEL
996 REAL, DIMENSION(KTS:KTE) :: NSCNG ! CHANGE N CONVERSION TO GRAUPEL DUE TO COLLECTION DROPLETS BY SNOW
997 REAL, DIMENSION(KTS:KTE) :: NGRACS ! CHANGE N CONVERSION TO GRAUPEL DUE TO COLLECTION RAIN BY SNOW
998 REAL, DIMENSION(KTS:KTE) :: NGMLTG ! CHANGE N MELTING GRAUPEL
999 REAL, DIMENSION(KTS:KTE) :: NGMLTR ! CHANGE N MELTING GRAUPEL TO RAIN
1000 REAL, DIMENSION(KTS:KTE) :: NSUBG ! CHANGE N SUB/DEP OF GRAUPEL
1001 REAL, DIMENSION(KTS:KTE) :: PSACR ! CONVERSION DUE TO COLL OF SNOW BY RAIN
1002 REAL, DIMENSION(KTS:KTE) :: NMULTG ! ICE MULT DUE TO ACC DROPLETS BY GRAUPEL
1003 REAL, DIMENSION(KTS:KTE) :: NMULTRG ! ICE MULT DUE TO ACC RAIN BY GRAUPEL
1004 REAL, DIMENSION(KTS:KTE) :: QMULTG ! CHANGE Q DUE TO ICE MULT DROPLETS/GRAUPEL
1005 REAL, DIMENSION(KTS:KTE) :: QMULTRG ! CHANGE Q DUE TO ICE MULT RAIN/GRAUPEL
1007 ! TIME-VARYING ATMOSPHERIC PARAMETERS
1009 REAL, DIMENSION(KTS:KTE) :: KAP ! THERMAL CONDUCTIVITY OF AIR
1010 REAL, DIMENSION(KTS:KTE) :: EVS ! SATURATION VAPOR PRESSURE
1011 REAL, DIMENSION(KTS:KTE) :: EIS ! ICE SATURATION VAPOR PRESSURE
1012 REAL, DIMENSION(KTS:KTE) :: QVS ! SATURATION MIXING RATIO
1013 REAL, DIMENSION(KTS:KTE) :: QVI ! ICE SATURATION MIXING RATIO
1014 REAL, DIMENSION(KTS:KTE) :: QVQVS ! SAUTRATION RATIO
1015 REAL, DIMENSION(KTS:KTE) :: QVQVSI! ICE SATURAION RATIO
1016 REAL, DIMENSION(KTS:KTE) :: DV ! DIFFUSIVITY OF WATER VAPOR IN AIR
1017 REAL, DIMENSION(KTS:KTE) :: XXLS ! LATENT HEAT OF SUBLIMATION
1018 REAL, DIMENSION(KTS:KTE) :: XXLV ! LATENT HEAT OF VAPORIZATION
1019 REAL, DIMENSION(KTS:KTE) :: CPM ! SPECIFIC HEAT AT CONST PRESSURE FOR MOIST AIR
1020 REAL, DIMENSION(KTS:KTE) :: MU ! VISCOCITY OF AIR
1021 REAL, DIMENSION(KTS:KTE) :: SC ! SCHMIDT NUMBER
1022 REAL, DIMENSION(KTS:KTE) :: XLF ! LATENT HEAT OF FREEZING
1023 REAL, DIMENSION(KTS:KTE) :: RHO ! AIR DENSITY
1024 REAL, DIMENSION(KTS:KTE) :: AB ! CORRECTION TO CONDENSATION RATE DUE TO LATENT HEATING
1025 REAL, DIMENSION(KTS:KTE) :: ABI ! CORRECTION TO DEPOSITION RATE DUE TO LATENT HEATING
1027 ! TIME-VARYING MICROPHYSICS PARAMETERS
1029 REAL, DIMENSION(KTS:KTE) :: DAP ! DIFFUSIVITY OF AEROSOL
1030 REAL NACNT ! NUMBER OF CONTACT IN
1031 REAL FMULT ! TEMP.-DEP. PARAMETER FOR RIME-SPLINTERING
1032 REAL COFFI ! ICE AUTOCONVERSION PARAMETER
1034 ! FALL SPEED WORKING VARIABLES (DEFINED IN CODE)
1036 REAL, DIMENSION(KTS:KTE) :: DUMI,DUMR,DUMFNI,DUMG,DUMFNG
1037 REAL UNI, UMI,UMR
1038 REAL, DIMENSION(KTS:KTE) :: FR, FI, FNI,FG,FNG
1039 REAL RGVM
1040 REAL, DIMENSION(KTS:KTE) :: FALOUTR,FALOUTI,FALOUTNI
1041 REAL FALTNDR,FALTNDI,FALTNDNI,RHO2
1042 REAL, DIMENSION(KTS:KTE) :: DUMQS,DUMFNS
1043 REAL UMS,UNS
1044 REAL, DIMENSION(KTS:KTE) :: FS,FNS, FALOUTS,FALOUTNS,FALOUTG,FALOUTNG
1045 REAL FALTNDS,FALTNDNS,UNR,FALTNDG,FALTNDNG
1046 REAL, DIMENSION(KTS:KTE) :: DUMC,DUMFNC
1047 REAL UNC,UMC,UNG,UMG
1048 REAL, DIMENSION(KTS:KTE) :: FC,FALOUTC,FALOUTNC
1049 REAL FALTNDC,FALTNDNC
1050 REAL, DIMENSION(KTS:KTE) :: FNC,DUMFNR,FALOUTNR
1051 REAL FALTNDNR
1052 REAL, DIMENSION(KTS:KTE) :: FNR
1054 ! FALL-SPEED PARAMETER 'A' WITH AIR DENSITY CORRECTION
1056 REAL, DIMENSION(KTS:KTE) :: AIN,ARN,ASN,ACN,AGN
1058 ! EXTERNAL FUNCTION CALL RETURN VARIABLES
1060 ! REAL GAMMA, ! EULER GAMMA FUNCTION
1061 ! REAL POLYSVP, ! SAT. PRESSURE FUNCTION
1062 ! REAL DERF1 ! ERROR FUNCTION
1064 ! DUMMY VARIABLES
1066 REAL DUM,DUM1,DUM2,DUMT,DUMQV,DUMQSS,DUMQSI,DUMS
1068 ! PROGNOSTIC SUPERSATURATION
1070 REAL DQSDT ! CHANGE OF SAT. MIX. RAT. WITH TEMPERATURE
1071 REAL DQSIDT ! CHANGE IN ICE SAT. MIXING RAT. WITH T
1072 REAL EPSI ! 1/PHASE REL. TIME (SEE M2005), ICE
1073 REAL EPSS ! 1/PHASE REL. TIME (SEE M2005), SNOW
1074 REAL EPSR ! 1/PHASE REL. TIME (SEE M2005), RAIN
1075 REAL EPSG ! 1/PHASE REL. TIME (SEE M2005), GRAUPEL
1077 ! NEW DROPLET ACTIVATION VARIABLES
1078 REAL TAUC ! PHASE REL. TIME (SEE M2005), DROPLETS
1079 REAL TAUR ! PHASE REL. TIME (SEE M2005), RAIN
1080 REAL TAUI ! PHASE REL. TIME (SEE M2005), CLOUD ICE
1081 REAL TAUS ! PHASE REL. TIME (SEE M2005), SNOW
1082 REAL TAUG ! PHASE REL. TIME (SEE M2005), GRAUPEL
1083 REAL DUMACT,DUM3
1085 ! COUNTING/INDEX VARIABLES
1087 INTEGER K,NSTEP,N ! ,I
1089 ! LTRUE IS ONLY USED TO SPEED UP THE CODE !!
1090 ! LTRUE, SWITCH = 0, NO HYDROMETEORS IN COLUMN,
1091 ! = 1, HYDROMETEORS IN COLUMN
1093 INTEGER LTRUE
1095 ! DROPLET ACTIVATION/FREEZING AEROSOL
1098 REAL CT ! DROPLET ACTIVATION PARAMETER
1099 REAL TEMP1 ! DUMMY TEMPERATURE
1100 REAL SAT1 ! DUMMY SATURATION
1101 REAL SIGVL ! SURFACE TENSION LIQ/VAPOR
1102 REAL KEL ! KELVIN PARAMETER
1103 REAL KC2 ! TOTAL ICE NUCLEATION RATE
1105 REAL CRY,KRY ! AEROSOL ACTIVATION PARAMETERS
1107 ! MORE WORKING/DUMMY VARIABLES
1109 REAL DUMQI,DUMNI,DC0,DS0,DG0
1110 REAL DUMQC,DUMQR,RATIO,SUM_DEP,FUDGEF
1112 ! EFFECTIVE VERTICAL VELOCITY (M/S)
1113 REAL WEF
1115 ! WORKING PARAMETERS FOR ICE NUCLEATION
1117 REAL ANUC,BNUC
1119 ! WORKING PARAMETERS FOR AEROSOL ACTIVATION
1121 REAL AACT,GAMM,GG,PSI,ETA1,ETA2,SM1,SM2,SMAX,UU1,UU2,ALPHA
1123 ! DUMMY SIZE DISTRIBUTION PARAMETERS
1125 REAL DLAMS,DLAMR,DLAMI,DLAMC,DLAMG,LAMMAX,LAMMIN
1127 INTEGER IDROP
1129 ! FOR WRF-CHEM
1130 REAL, DIMENSION(KTS:KTE)::C2PREC,CSED,ISED,SSED,GSED,RSED
1132 ! comment lines for wrf-chem since these are intent(in) in that case
1133 ! REAL, DIMENSION(KTS:KTE) :: NC3DTEN ! CLOUD DROPLET NUMBER CONCENTRATION (1/KG/S)
1134 ! REAL, DIMENSION(KTS:KTE) :: NC3D ! CLOUD DROPLET NUMBER CONCENTRATION (1/KG)
1136 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1138 ! SET LTRUE INITIALLY TO 0
1140 LTRUE = 0
1142 ! ATMOSPHERIC PARAMETERS THAT VARY IN TIME AND HEIGHT
1143 DO K = KTS,KTE
1145 ! NC3DTEN LOCAL ARRAY INITIALIZED
1146 NC3DTEN(K) = 0.
1147 ! INITIALIZE VARIABLES FOR WRF-CHEM OUTPUT TO ZERO
1149 C2PREC(K)=0.
1150 CSED(K)=0.
1151 ISED(K)=0.
1152 SSED(K)=0.
1153 GSED(K)=0.
1154 RSED(K)=0.
1156 ! LATENT HEAT OF VAPORATION
1158 XXLV(K) = 3.1484E6-2370.*T3D(K)
1160 ! LATENT HEAT OF SUBLIMATION
1162 XXLS(K) = 3.15E6-2370.*T3D(K)+0.3337E6
1164 CPM(K) = CP*(1.+0.887*QV3D(K))
1166 ! SATURATION VAPOR PRESSURE AND MIXING RATIO
1168 ! hm, add fix for low pressure, 5/12/10
1169 EVS(K) = min(0.99*pres(k),POLYSVP(T3D(K),0)) ! PA
1170 EIS(K) = min(0.99*pres(k),POLYSVP(T3D(K),1)) ! PA
1172 ! MAKE SURE ICE SATURATION DOESN'T EXCEED WATER SAT. NEAR FREEZING
1174 IF (EIS(K).GT.EVS(K)) EIS(K) = EVS(K)
1176 QVS(K) = EP_2*EVS(K)/(PRES(K)-EVS(K))
1177 QVI(K) = EP_2*EIS(K)/(PRES(K)-EIS(K))
1179 QVQVS(K) = QV3D(K)/QVS(K)
1180 QVQVSI(K) = QV3D(K)/QVI(K)
1182 ! AIR DENSITY
1184 RHO(K) = PRES(K)/(R*T3D(K))
1186 ! ADD NUMBER CONCENTRATION DUE TO CUMULUS TENDENCY
1187 ! ASSUME N0 ASSOCIATED WITH CUMULUS PARAM RAIN IS 10^7 M^-4
1188 ! ASSUME N0 ASSOCIATED WITH CUMULUS PARAM SNOW IS 2 X 10^7 M^-4
1189 ! FOR DETRAINED CLOUD ICE, ASSUME MEAN VOLUME DIAM OF 80 MICRON
1191 IF (QRCU1D(K).GE.1.E-10) THEN
1192 DUM=1.8e5*(QRCU1D(K)*DT/(PI*RHOW*RHO(K)**3))**0.25
1193 NR3D(K)=NR3D(K)+DUM
1194 END IF
1195 IF (QSCU1D(K).GE.1.E-10) THEN
1196 DUM=3.e5*(QSCU1D(K)*DT/(CONS1*RHO(K)**3))**(1./(DS+1.))
1197 NS3D(K)=NS3D(K)+DUM
1198 END IF
1199 IF (QICU1D(K).GE.1.E-10) THEN
1200 DUM=QICU1D(K)*DT/(CI*(80.E-6)**DI)
1201 NI3D(K)=NI3D(K)+DUM
1202 END IF
1204 ! AT SUBSATURATION, REMOVE SMALL AMOUNTS OF CLOUD/PRECIP WATER
1205 ! hm modify 7/0/09 change limit to 1.e-8
1207 IF (QVQVS(K).LT.0.9) THEN
1208 IF (QR3D(K).LT.1.E-8) THEN
1209 QV3D(K)=QV3D(K)+QR3D(K)
1210 T3D(K)=T3D(K)-QR3D(K)*XXLV(K)/CPM(K)
1211 QR3D(K)=0.
1212 END IF
1213 IF (QC3D(K).LT.1.E-8) THEN
1214 QV3D(K)=QV3D(K)+QC3D(K)
1215 T3D(K)=T3D(K)-QC3D(K)*XXLV(K)/CPM(K)
1216 QC3D(K)=0.
1217 END IF
1218 END IF
1220 IF (QVQVSI(K).LT.0.9) THEN
1221 IF (QI3D(K).LT.1.E-8) THEN
1222 QV3D(K)=QV3D(K)+QI3D(K)
1223 T3D(K)=T3D(K)-QI3D(K)*XXLS(K)/CPM(K)
1224 QI3D(K)=0.
1225 END IF
1226 IF (QNI3D(K).LT.1.E-8) THEN
1227 QV3D(K)=QV3D(K)+QNI3D(K)
1228 T3D(K)=T3D(K)-QNI3D(K)*XXLS(K)/CPM(K)
1229 QNI3D(K)=0.
1230 END IF
1231 IF (QG3D(K).LT.1.E-8) THEN
1232 QV3D(K)=QV3D(K)+QG3D(K)
1233 T3D(K)=T3D(K)-QG3D(K)*XXLS(K)/CPM(K)
1234 QG3D(K)=0.
1235 END IF
1236 END IF
1238 ! HEAT OF FUSION
1240 XLF(K) = XXLS(K)-XXLV(K)
1242 !..................................................................
1243 ! IF MIXING RATIO < QSMALL SET MIXING RATIO AND NUMBER CONC TO ZERO
1245 IF (QC3D(K).LT.QSMALL) THEN
1246 QC3D(K) = 0.
1247 NC3D(K) = 0.
1248 EFFC(K) = 0.
1249 END IF
1250 IF (QR3D(K).LT.QSMALL) THEN
1251 QR3D(K) = 0.
1252 NR3D(K) = 0.
1253 EFFR(K) = 0.
1254 END IF
1255 IF (QI3D(K).LT.QSMALL) THEN
1256 QI3D(K) = 0.
1257 NI3D(K) = 0.
1258 EFFI(K) = 0.
1259 END IF
1260 IF (QNI3D(K).LT.QSMALL) THEN
1261 QNI3D(K) = 0.
1262 NS3D(K) = 0.
1263 EFFS(K) = 0.
1264 END IF
1265 IF (QG3D(K).LT.QSMALL) THEN
1266 QG3D(K) = 0.
1267 NG3D(K) = 0.
1268 EFFG(K) = 0.
1269 END IF
1271 ! INITIALIZE SEDIMENTATION TENDENCIES FOR MIXING RATIO
1273 QRSTEN(K) = 0.
1274 QISTEN(K) = 0.
1275 QNISTEN(K) = 0.
1276 QCSTEN(K) = 0.
1277 QGSTEN(K) = 0.
1279 !..................................................................
1280 ! MICROPHYSICS PARAMETERS VARYING IN TIME/HEIGHT
1282 ! fix 053011
1283 MU(K) = 1.496E-6*T3D(K)**1.5/(T3D(K)+120.)
1285 ! FALL SPEED WITH DENSITY CORRECTION (HEYMSFIELD AND BENSSEMER 2006)
1287 DUM = (RHOSU/RHO(K))**0.54
1289 ! fix 053011
1290 ! AIN(K) = DUM*AI
1291 ! AA revision 4/1/11: Ikawa and Saito 1991 air-density correction
1292 AIN(K) = (RHOSU/RHO(K))**0.35*AI
1293 ARN(K) = DUM*AR
1294 ASN(K) = DUM*AS
1295 ! ACN(K) = DUM*AC
1296 ! AA revision 4/1/11: temperature-dependent Stokes fall speed
1297 ACN(K) = G*RHOW/(18.*MU(K))
1298 ! HM ADD GRAUPEL 8/28/06
1299 AGN(K) = DUM*AG
1301 !hm 4/7/09 bug fix, initialize lami to prevent later division by zero
1302 LAMI(K)=0.
1304 !..................................
1305 ! IF THERE IS NO CLOUD/PRECIP WATER, AND IF SUBSATURATED, THEN SKIP MICROPHYSICS
1306 ! FOR THIS LEVEL
1308 IF (QC3D(K).LT.QSMALL.AND.QI3D(K).LT.QSMALL.AND.QNI3D(K).LT.QSMALL &
1309 .AND.QR3D(K).LT.QSMALL.AND.QG3D(K).LT.QSMALL) THEN
1310 IF (T3D(K).LT.273.15.AND.QVQVSI(K).LT.0.999) GOTO 200
1311 IF (T3D(K).GE.273.15.AND.QVQVS(K).LT.0.999) GOTO 200
1312 END IF
1314 ! THERMAL CONDUCTIVITY FOR AIR
1316 ! fix 053011
1317 KAP(K) = 1.414E3*MU(K)
1319 ! DIFFUSIVITY OF WATER VAPOR
1321 DV(K) = 8.794E-5*T3D(K)**1.81/PRES(K)
1323 ! SCHMIT NUMBER
1325 ! fix 053011
1326 SC(K) = MU(K)/(RHO(K)*DV(K))
1328 ! PSYCHOMETIC CORRECTIONS
1330 ! RATE OF CHANGE SAT. MIX. RATIO WITH TEMPERATURE
1332 DUM = (RV*T3D(K)**2)
1334 DQSDT = XXLV(K)*QVS(K)/DUM
1335 DQSIDT = XXLS(K)*QVI(K)/DUM
1337 ABI(K) = 1.+DQSIDT*XXLS(K)/CPM(K)
1338 AB(K) = 1.+DQSDT*XXLV(K)/CPM(K)
1340 !
1341 !.....................................................................
1342 !.....................................................................
1343 ! CASE FOR TEMPERATURE ABOVE FREEZING
1345 IF (T3D(K).GE.273.15) THEN
1347 !......................................................................
1348 !HM ADD, ALLOW FOR CONSTANT DROPLET NUMBER
1349 ! INUM = 0, PREDICT DROPLET NUMBER
1350 ! INUM = 1, SET CONSTANT DROPLET NUMBER
1352 IF (iinum.EQ.1) THEN
1353 ! CONVERT NDCNST FROM CM-3 TO KG-1
1354 NC3D(K)=NDCNST*1.E6/RHO(K)
1355 END IF
1357 ! GET SIZE DISTRIBUTION PARAMETERS
1359 ! MELT VERY SMALL SNOW AND GRAUPEL MIXING RATIOS, ADD TO RAIN
1360 IF (QNI3D(K).LT.1.E-6) THEN
1361 QR3D(K)=QR3D(K)+QNI3D(K)
1362 NR3D(K)=NR3D(K)+NS3D(K)
1363 T3D(K)=T3D(K)-QNI3D(K)*XLF(K)/CPM(K)
1364 QNI3D(K) = 0.
1365 NS3D(K) = 0.
1366 END IF
1367 IF (QG3D(K).LT.1.E-6) THEN
1368 QR3D(K)=QR3D(K)+QG3D(K)
1369 NR3D(K)=NR3D(K)+NG3D(K)
1370 T3D(K)=T3D(K)-QG3D(K)*XLF(K)/CPM(K)
1371 QG3D(K) = 0.
1372 NG3D(K) = 0.
1373 END IF
1375 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
1377 ! MAKE SURE NUMBER CONCENTRATIONS AREN'T NEGATIVE
1379 NS3D(K) = MAX(0.,NS3D(K))
1380 NC3D(K) = MAX(0.,NC3D(K))
1381 NR3D(K) = MAX(0.,NR3D(K))
1382 NG3D(K) = MAX(0.,NG3D(K))
1384 !......................................................................
1385 ! RAIN
1387 IF (QR3D(K).GE.QSMALL) THEN
1388 LAMR(K) = (PI*RHOW*NR3D(K)/QR3D(K))**(1./3.)
1389 N0RR(K) = NR3D(K)*LAMR(K)
1391 ! CHECK FOR SLOPE
1393 ! ADJUST VARS
1395 IF (LAMR(K).LT.LAMMINR) THEN
1397 LAMR(K) = LAMMINR
1399 N0RR(K) = LAMR(K)**4*QR3D(K)/(PI*RHOW)
1401 NR3D(K) = N0RR(K)/LAMR(K)
1402 ELSE IF (LAMR(K).GT.LAMMAXR) THEN
1403 LAMR(K) = LAMMAXR
1404 N0RR(K) = LAMR(K)**4*QR3D(K)/(PI*RHOW)
1406 NR3D(K) = N0RR(K)/LAMR(K)
1407 END IF
1408 END IF
1410 !......................................................................
1411 ! CLOUD DROPLETS
1413 ! MARTIN ET AL. (1994) FORMULA FOR PGAM
1415 IF (QC3D(K).GE.QSMALL) THEN
1417 DUM = PRES(K)/(287.15*T3D(K))
1418 PGAM(K)=0.0005714*(NC3D(K)/1.E6*DUM)+0.2714
1419 PGAM(K)=1./(PGAM(K)**2)-1.
1420 PGAM(K)=MAX(PGAM(K),2.)
1421 PGAM(K)=MIN(PGAM(K),10.)
1423 ! CALCULATE LAMC
1425 LAMC(K) = (CONS26*NC3D(K)*GAMMA(PGAM(K)+4.)/ &
1426 (QC3D(K)*GAMMA(PGAM(K)+1.)))**(1./3.)
1428 ! LAMMIN, 60 MICRON DIAMETER
1429 ! LAMMAX, 1 MICRON
1431 LAMMIN = (PGAM(K)+1.)/60.E-6
1432 LAMMAX = (PGAM(K)+1.)/1.E-6
1434 IF (LAMC(K).LT.LAMMIN) THEN
1435 LAMC(K) = LAMMIN
1437 NC3D(K) = EXP(3.*LOG(LAMC(K))+LOG(QC3D(K))+ &
1438 LOG(GAMMA(PGAM(K)+1.))-LOG(GAMMA(PGAM(K)+4.)))/CONS26
1439 ELSE IF (LAMC(K).GT.LAMMAX) THEN
1440 LAMC(K) = LAMMAX
1442 NC3D(K) = EXP(3.*LOG(LAMC(K))+LOG(QC3D(K))+ &
1443 LOG(GAMMA(PGAM(K)+1.))-LOG(GAMMA(PGAM(K)+4.)))/CONS26
1445 END IF
1447 END IF
1449 !......................................................................
1450 ! SNOW
1452 IF (QNI3D(K).GE.QSMALL) THEN
1453 LAMS(K) = (CONS1*NS3D(K)/QNI3D(K))**(1./DS)
1454 N0S(K) = NS3D(K)*LAMS(K)
1456 ! CHECK FOR SLOPE
1458 ! ADJUST VARS
1460 IF (LAMS(K).LT.LAMMINS) THEN
1461 LAMS(K) = LAMMINS
1462 N0S(K) = LAMS(K)**4*QNI3D(K)/CONS1
1464 NS3D(K) = N0S(K)/LAMS(K)
1466 ELSE IF (LAMS(K).GT.LAMMAXS) THEN
1468 LAMS(K) = LAMMAXS
1469 N0S(K) = LAMS(K)**4*QNI3D(K)/CONS1
1471 NS3D(K) = N0S(K)/LAMS(K)
1472 END IF
1473 END IF
1475 !......................................................................
1476 ! GRAUPEL
1478 IF (QG3D(K).GE.QSMALL) THEN
1479 LAMG(K) = (CONS2*NG3D(K)/QG3D(K))**(1./DG)
1480 N0G(K) = NG3D(K)*LAMG(K)
1482 ! ADJUST VARS
1484 IF (LAMG(K).LT.LAMMING) THEN
1485 LAMG(K) = LAMMING
1486 N0G(K) = LAMG(K)**4*QG3D(K)/CONS2
1488 NG3D(K) = N0G(K)/LAMG(K)
1490 ELSE IF (LAMG(K).GT.LAMMAXG) THEN
1492 LAMG(K) = LAMMAXG
1493 N0G(K) = LAMG(K)**4*QG3D(K)/CONS2
1495 NG3D(K) = N0G(K)/LAMG(K)
1496 END IF
1497 END IF
1499 !.....................................................................
1500 ! ZERO OUT PROCESS RATES
1502 PRC(K) = 0.
1503 NPRC(K) = 0.
1504 NPRC1(K) = 0.
1505 PRA(K) = 0.
1506 NPRA(K) = 0.
1507 NRAGG(K) = 0.
1508 PSMLT(K) = 0.
1509 NSMLTS(K) = 0.
1510 NSMLTR(K) = 0.
1511 EVPMS(K) = 0.
1512 PCC(K) = 0.
1513 PRE(K) = 0.
1514 NSUBC(K) = 0.
1515 NSUBR(K) = 0.
1516 PRACG(K) = 0.
1517 NPRACG(K) = 0.
1518 PSMLT(K) = 0.
1519 EVPMS(K) = 0.
1520 PGMLT(K) = 0.
1521 EVPMG(K) = 0.
1522 PRACS(K) = 0.
1523 NPRACS(K) = 0.
1524 NGMLTG(K) = 0.
1525 NGMLTR(K) = 0.
1527 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1528 ! CALCULATION OF MICROPHYSICAL PROCESS RATES, T > 273.15 K
1530 !.................................................................
1531 !.......................................................................
1532 ! AUTOCONVERSION OF CLOUD LIQUID WATER TO RAIN
1533 ! FORMULA FROM BEHENG (1994)
1534 ! USING NUMERICAL SIMULATION OF STOCHASTIC COLLECTION EQUATION
1535 ! AND INITIAL CLOUD DROPLET SIZE DISTRIBUTION SPECIFIED
1536 ! AS A GAMMA DISTRIBUTION
1538 ! USE MINIMUM VALUE OF 1.E-6 TO PREVENT FLOATING POINT ERROR
1540 IF (QC3D(K).GE.1.E-6) THEN
1542 ! HM ADD 12/13/06, REPLACE WITH NEWER FORMULA
1543 ! FROM KHAIROUTDINOV AND KOGAN 2000, MWR
1545 PRC(K)=1350.*QC3D(K)**2.47* &
1546 (NC3D(K)/1.e6*RHO(K))**(-1.79)
1548 ! note: nprc1 is change in Nr,
1549 ! nprc is change in Nc
1551 NPRC1(K) = PRC(K)/CONS29
1552 NPRC(K) = PRC(K)/(QC3D(k)/NC3D(K))
1554 ! hm bug fix 3/20/12
1555 NPRC(K) = MIN(NPRC(K),NC3D(K)/DT)
1556 NPRC1(K) = MIN(NPRC1(K),NPRC(K))
1558 END IF
1560 !.......................................................................
1561 ! HM ADD 12/13/06, COLLECTION OF SNOW BY RAIN ABOVE FREEZING
1562 ! FORMULA FROM IKAWA AND SAITO (1991)
1564 IF (QR3D(K).GE.1.E-8.AND.QNI3D(K).GE.1.E-8) THEN
1566 UMS = ASN(K)*CONS3/(LAMS(K)**BS)
1567 UMR = ARN(K)*CONS4/(LAMR(K)**BR)
1568 UNS = ASN(K)*CONS5/LAMS(K)**BS
1569 UNR = ARN(K)*CONS6/LAMR(K)**BR
1571 ! SET REASLISTIC LIMITS ON FALLSPEEDS
1573 ! bug fix, 10/08/09
1574 dum=(rhosu/rho(k))**0.54
1575 UMS=MIN(UMS,1.2*dum)
1576 UNS=MIN(UNS,1.2*dum)
1577 UMR=MIN(UMR,9.1*dum)
1578 UNR=MIN(UNR,9.1*dum)
1580 PRACS(K) = CONS31*(((1.2*UMR-0.95*UMS)**2+ &
1581 0.08*UMS*UMR)**0.5*RHO(K)* &
1582 N0RR(K)*N0S(K)/LAMS(K)**3* &
1583 (5./(LAMS(K)**3*LAMR(K))+ &
1584 2./(LAMS(K)**2*LAMR(K)**2)+ &
1585 0.5/(LAMS(K)*LAMR(K)**3)))
1587 ! fix 053011, npracs no longer subtracted from snow
1588 ! NPRACS(K) = CONS32*RHO(K)*(1.7*(UNR-UNS)**2+ &
1589 ! 0.3*UNR*UNS)**0.5*N0RR(K)*N0S(K)* &
1590 ! (1./(LAMR(K)**3*LAMS(K))+ &
1591 ! 1./(LAMR(K)**2*LAMS(K)**2)+ &
1592 ! 1./(LAMR(K)*LAMS(K)**3))
1594 END IF
1596 ! ADD COLLECTION OF GRAUPEL BY RAIN ABOVE FREEZING
1597 ! ASSUME ALL RAIN COLLECTION BY GRAUPEL ABOVE FREEZING IS SHED
1598 ! ASSUME SHED DROPS ARE 1 MM IN SIZE
1600 IF (QR3D(K).GE.1.E-8.AND.QG3D(K).GE.1.E-8) THEN
1602 UMG = AGN(K)*CONS7/(LAMG(K)**BG)
1603 UMR = ARN(K)*CONS4/(LAMR(K)**BR)
1604 UNG = AGN(K)*CONS8/LAMG(K)**BG
1605 UNR = ARN(K)*CONS6/LAMR(K)**BR
1607 ! SET REASLISTIC LIMITS ON FALLSPEEDS
1608 ! bug fix, 10/08/09
1609 dum=(rhosu/rho(k))**0.54
1610 UMG=MIN(UMG,20.*dum)
1611 UNG=MIN(UNG,20.*dum)
1612 UMR=MIN(UMR,9.1*dum)
1613 UNR=MIN(UNR,9.1*dum)
1615 ! PRACG IS MIXING RATIO OF RAIN PER SEC COLLECTED BY GRAUPEL/HAIL
1616 PRACG(K) = CONS41*(((1.2*UMR-0.95*UMG)**2+ &
1617 0.08*UMG*UMR)**0.5*RHO(K)* &
1618 N0RR(K)*N0G(K)/LAMR(K)**3* &
1619 (5./(LAMR(K)**3*LAMG(K))+ &
1620 2./(LAMR(K)**2*LAMG(K)**2)+ &
1621 0.5/(LAMR(k)*LAMG(k)**3)))
1623 ! ASSUME 1 MM DROPS ARE SHED, GET NUMBER SHED PER SEC
1625 DUM = PRACG(K)/5.2E-7
1627 NPRACG(K) = CONS32*RHO(K)*(1.7*(UNR-UNG)**2+ &
1628 0.3*UNR*UNG)**0.5*N0RR(K)*N0G(K)* &
1629 (1./(LAMR(K)**3*LAMG(K))+ &
1630 1./(LAMR(K)**2*LAMG(K)**2)+ &
1631 1./(LAMR(K)*LAMG(K)**3))
1633 NPRACG(K)=MAX(NPRACG(K)-DUM,0.)
1635 END IF
1637 !.......................................................................
1638 ! ACCRETION OF CLOUD LIQUID WATER BY RAIN
1639 ! CONTINUOUS COLLECTION EQUATION WITH
1640 ! GRAVITATIONAL COLLECTION KERNEL, DROPLET FALL SPEED NEGLECTED
1642 IF (QR3D(K).GE.1.E-8 .AND. QC3D(K).GE.1.E-8) THEN
1644 ! 12/13/06 HM ADD, REPLACE WITH NEWER FORMULA FROM
1645 ! KHAIROUTDINOV AND KOGAN 2000, MWR
1647 DUM=(QC3D(K)*QR3D(K))
1648 PRA(K) = 67.*(DUM)**1.15
1649 NPRA(K) = PRA(K)/(QC3D(K)/NC3D(K))
1651 END IF
1652 !.......................................................................
1653 ! SELF-COLLECTION OF RAIN DROPS
1654 ! FROM BEHENG(1994)
1655 ! FROM NUMERICAL SIMULATION OF THE STOCHASTIC COLLECTION EQUATION
1656 ! AS DESCRINED ABOVE FOR AUTOCONVERSION
1658 IF (QR3D(K).GE.1.E-8) THEN
1659 ! include breakup add 10/09/09
1660 dum1=300.e-6
1661 if (1./lamr(k).lt.dum1) then
1662 dum=1.
1663 else if (1./lamr(k).ge.dum1) then
1664 dum=2.-exp(2300.*(1./lamr(k)-dum1))
1665 end if
1666 ! NRAGG(K) = -8.*NR3D(K)*QR3D(K)*RHO(K)
1667 NRAGG(K) = -5.78*dum*NR3D(K)*QR3D(K)*RHO(K)
1668 END IF
1670 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1671 ! CALCULATE EVAP OF RAIN (RUTLEDGE AND HOBBS 1983)
1673 IF (QR3D(K).GE.QSMALL) THEN
1674 EPSR = 2.*PI*N0RR(K)*RHO(K)*DV(K)* &
1675 (F1R/(LAMR(K)*LAMR(K))+ &
1676 F2R*(ARN(K)*RHO(K)/MU(K))**0.5* &
1677 SC(K)**(1./3.)*CONS9/ &
1678 (LAMR(K)**CONS34))
1679 ELSE
1680 EPSR = 0.
1681 END IF
1683 ! NO CONDENSATION ONTO RAIN, ONLY EVAP ALLOWED
1685 IF (QV3D(K).LT.QVS(K)) THEN
1686 PRE(K) = EPSR*(QV3D(K)-QVS(K))/AB(K)
1687 PRE(K) = MIN(PRE(K),0.)
1688 ELSE
1689 PRE(K) = 0.
1690 END IF
1692 !.......................................................................
1693 ! MELTING OF SNOW
1695 ! SNOW MAY PERSITS ABOVE FREEZING, FORMULA FROM RUTLEDGE AND HOBBS, 1984
1696 ! IF WATER SUPERSATURATION, SNOW MELTS TO FORM RAIN
1698 IF (QNI3D(K).GE.1.E-8) THEN
1700 ! fix 053011
1701 ! HM, MODIFY FOR V3.2, ADD ACCELERATED MELTING DUE TO COLLISION WITH RAIN
1702 ! DUM = -CPW/XLF(K)*T3D(K)*PRACS(K)
1703 DUM = -CPW/XLF(K)*(T3D(K)-273.15)*PRACS(K)
1705 PSMLT(K)=2.*PI*N0S(K)*KAP(K)*(273.15-T3D(K))/ &
1706 XLF(K)*RHO(K)*(F1S/(LAMS(K)*LAMS(K))+ &
1707 F2S*(ASN(K)*RHO(K)/MU(K))**0.5* &
1708 SC(K)**(1./3.)*CONS10/ &
1709 (LAMS(K)**CONS35))+DUM
1711 ! IN WATER SUBSATURATION, SNOW MELTS AND EVAPORATES
1713 IF (QVQVS(K).LT.1.) THEN
1714 EPSS = 2.*PI*N0S(K)*RHO(K)*DV(K)* &
1715 (F1S/(LAMS(K)*LAMS(K))+ &
1716 F2S*(ASN(K)*RHO(K)/MU(K))**0.5* &
1717 SC(K)**(1./3.)*CONS10/ &
1718 (LAMS(K)**CONS35))
1719 ! hm fix 8/4/08
1720 EVPMS(K) = (QV3D(K)-QVS(K))*EPSS/AB(K)
1721 EVPMS(K) = MAX(EVPMS(K),PSMLT(K))
1722 PSMLT(K) = PSMLT(K)-EVPMS(K)
1723 END IF
1724 END IF
1726 !.......................................................................
1727 ! MELTING OF GRAUPEL
1729 ! GRAUPEL MAY PERSITS ABOVE FREEZING, FORMULA FROM RUTLEDGE AND HOBBS, 1984
1730 ! IF WATER SUPERSATURATION, GRAUPEL MELTS TO FORM RAIN
1732 IF (QG3D(K).GE.1.E-8) THEN
1734 ! fix 053011
1735 ! HM, MODIFY FOR V3.2, ADD ACCELERATED MELTING DUE TO COLLISION WITH RAIN
1736 ! DUM = -CPW/XLF(K)*T3D(K)*PRACG(K)
1737 DUM = -CPW/XLF(K)*(T3D(K)-273.15)*PRACG(K)
1739 PGMLT(K)=2.*PI*N0G(K)*KAP(K)*(273.15-T3D(K))/ &
1740 XLF(K)*RHO(K)*(F1S/(LAMG(K)*LAMG(K))+ &
1741 F2S*(AGN(K)*RHO(K)/MU(K))**0.5* &
1742 SC(K)**(1./3.)*CONS11/ &
1743 (LAMG(K)**CONS36))+DUM
1745 ! IN WATER SUBSATURATION, GRAUPEL MELTS AND EVAPORATES
1747 IF (QVQVS(K).LT.1.) THEN
1748 EPSG = 2.*PI*N0G(K)*RHO(K)*DV(K)* &
1749 (F1S/(LAMG(K)*LAMG(K))+ &
1750 F2S*(AGN(K)*RHO(K)/MU(K))**0.5* &
1751 SC(K)**(1./3.)*CONS11/ &
1752 (LAMG(K)**CONS36))
1753 ! hm fix 8/4/08
1754 EVPMG(K) = (QV3D(K)-QVS(K))*EPSG/AB(K)
1755 EVPMG(K) = MAX(EVPMG(K),PGMLT(K))
1756 PGMLT(K) = PGMLT(K)-EVPMG(K)
1757 END IF
1758 END IF
1760 ! HM, V3.2
1761 ! RESET PRACG AND PRACS TO ZERO, THIS IS DONE BECAUSE THERE IS NO
1762 ! TRANSFER OF MASS FROM SNOW AND GRAUPEL TO RAIN DIRECTLY FROM COLLECTION
1763 ! ABOVE FREEZING, IT IS ONLY USED FOR ENHANCEMENT OF MELTING AND SHEDDING
1765 PRACG(K) = 0.
1766 PRACS(K) = 0.
1768 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1770 ! FOR CLOUD ICE, ONLY PROCESSES OPERATING AT T > 273.15 IS
1771 ! MELTING, WHICH IS ALREADY CONSERVED DURING PROCESS
1772 ! CALCULATION
1774 ! CONSERVATION OF QC
1776 DUM = (PRC(K)+PRA(K))*DT
1778 IF (DUM.GT.QC3D(K).AND.QC3D(K).GE.QSMALL) THEN
1780 RATIO = QC3D(K)/DUM
1782 PRC(K) = PRC(K)*RATIO
1783 PRA(K) = PRA(K)*RATIO
1785 END IF
1787 ! CONSERVATION OF SNOW
1789 DUM = (-PSMLT(K)-EVPMS(K)+PRACS(K))*DT
1791 IF (DUM.GT.QNI3D(K).AND.QNI3D(K).GE.QSMALL) THEN
1793 ! NO SOURCE TERMS FOR SNOW AT T > FREEZING
1794 RATIO = QNI3D(K)/DUM
1796 PSMLT(K) = PSMLT(K)*RATIO
1797 EVPMS(K) = EVPMS(K)*RATIO
1798 PRACS(K) = PRACS(K)*RATIO
1800 END IF
1802 ! CONSERVATION OF GRAUPEL
1804 DUM = (-PGMLT(K)-EVPMG(K)+PRACG(K))*DT
1806 IF (DUM.GT.QG3D(K).AND.QG3D(K).GE.QSMALL) THEN
1808 ! NO SOURCE TERM FOR GRAUPEL ABOVE FREEZING
1809 RATIO = QG3D(K)/DUM
1811 PGMLT(K) = PGMLT(K)*RATIO
1812 EVPMG(K) = EVPMG(K)*RATIO
1813 PRACG(K) = PRACG(K)*RATIO
1815 END IF
1817 ! CONSERVATION OF QR
1818 ! HM 12/13/06, ADDED CONSERVATION OF RAIN SINCE PRE IS NEGATIVE
1820 DUM = (-PRACS(K)-PRACG(K)-PRE(K)-PRA(K)-PRC(K)+PSMLT(K)+PGMLT(K))*DT
1822 IF (DUM.GT.QR3D(K).AND.QR3D(K).GE.QSMALL) THEN
1824 RATIO = (QR3D(K)/DT+PRACS(K)+PRACG(K)+PRA(K)+PRC(K)-PSMLT(K)-PGMLT(K))/ &
1825 (-PRE(K))
1826 PRE(K) = PRE(K)*RATIO
1828 END IF
1830 !....................................
1832 QV3DTEN(K) = QV3DTEN(K)+(-PRE(K)-EVPMS(K)-EVPMG(K))
1834 T3DTEN(K) = T3DTEN(K)+(PRE(K)*XXLV(K)+(EVPMS(K)+EVPMG(K))*XXLS(K)+&
1835 (PSMLT(K)+PGMLT(K)-PRACS(K)-PRACG(K))*XLF(K))/CPM(K)
1837 QC3DTEN(K) = QC3DTEN(K)+(-PRA(K)-PRC(K))
1838 QR3DTEN(K) = QR3DTEN(K)+(PRE(K)+PRA(K)+PRC(K)-PSMLT(K)-PGMLT(K)+PRACS(K)+PRACG(K))
1839 QNI3DTEN(K) = QNI3DTEN(K)+(PSMLT(K)+EVPMS(K)-PRACS(K))
1840 QG3DTEN(K) = QG3DTEN(K)+(PGMLT(K)+EVPMG(K)-PRACG(K))
1841 ! fix 053011
1842 ! NS3DTEN(K) = NS3DTEN(K)-NPRACS(K)
1843 ! HM, bug fix 5/12/08, npracg is subtracted from nr not ng
1844 ! NG3DTEN(K) = NG3DTEN(K)
1845 NC3DTEN(K) = NC3DTEN(K)+ (-NPRA(K)-NPRC(K))
1846 NR3DTEN(K) = NR3DTEN(K)+ (NPRC1(K)+NRAGG(K)-NPRACG(K))
1848 ! HM ADD, WRF-CHEM, ADD TENDENCIES FOR C2PREC
1850 C2PREC(K) = PRA(K)+PRC(K)
1851 IF (PRE(K).LT.0.) THEN
1852 DUM = PRE(K)*DT/QR3D(K)
1853 DUM = MAX(-1.,DUM)
1854 NSUBR(K) = DUM*NR3D(K)/DT
1855 END IF
1857 IF (EVPMS(K)+PSMLT(K).LT.0.) THEN
1858 DUM = (EVPMS(K)+PSMLT(K))*DT/QNI3D(K)
1859 DUM = MAX(-1.,DUM)
1860 NSMLTS(K) = DUM*NS3D(K)/DT
1861 END IF
1862 IF (PSMLT(K).LT.0.) THEN
1863 DUM = PSMLT(K)*DT/QNI3D(K)
1864 DUM = MAX(-1.0,DUM)
1865 NSMLTR(K) = DUM*NS3D(K)/DT
1866 END IF
1867 IF (EVPMG(K)+PGMLT(K).LT.0.) THEN
1868 DUM = (EVPMG(K)+PGMLT(K))*DT/QG3D(K)
1869 DUM = MAX(-1.,DUM)
1870 NGMLTG(K) = DUM*NG3D(K)/DT
1871 END IF
1872 IF (PGMLT(K).LT.0.) THEN
1873 DUM = PGMLT(K)*DT/QG3D(K)
1874 DUM = MAX(-1.0,DUM)
1875 NGMLTR(K) = DUM*NG3D(K)/DT
1876 END IF
1878 NS3DTEN(K) = NS3DTEN(K)+(NSMLTS(K))
1879 NG3DTEN(K) = NG3DTEN(K)+(NGMLTG(K))
1880 NR3DTEN(K) = NR3DTEN(K)+(NSUBR(K)-NSMLTR(K)-NGMLTR(K))
1882 300 CONTINUE
1884 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1885 ! NOW CALCULATE SATURATION ADJUSTMENT TO CONDENSE EXTRA VAPOR ABOVE
1886 ! WATER SATURATION
1888 DUMT = T3D(K)+DT*T3DTEN(K)
1889 DUMQV = QV3D(K)+DT*QV3DTEN(K)
1890 ! hm, add fix for low pressure, 5/12/10
1891 dum=min(0.99*pres(k),POLYSVP(DUMT,0))
1892 DUMQSS = EP_2*dum/(PRES(K)-dum)
1893 DUMQC = QC3D(K)+DT*QC3DTEN(K)
1894 DUMQC = MAX(DUMQC,0.)
1896 ! SATURATION ADJUSTMENT FOR LIQUID
1898 DUMS = DUMQV-DUMQSS
1899 PCC(K) = DUMS/(1.+XXLV(K)**2*DUMQSS/(CPM(K)*RV*DUMT**2))/DT
1900 IF (PCC(K)*DT+DUMQC.LT.0.) THEN
1901 PCC(K) = -DUMQC/DT
1902 END IF
1904 QV3DTEN(K) = QV3DTEN(K)-PCC(K)
1905 T3DTEN(K) = T3DTEN(K)+PCC(K)*XXLV(K)/CPM(K)
1906 QC3DTEN(K) = QC3DTEN(K)+PCC(K)
1908 !.......................................................................
1909 ! ACTIVATION OF CLOUD DROPLETS
1910 ! ACTIVATION OF DROPLET CURRENTLY NOT CALCULATED
1911 ! DROPLET CONCENTRATION IS SPECIFIED !!!!!
1913 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1914 ! SUBLIMATE, MELT, OR EVAPORATE NUMBER CONCENTRATION
1915 ! THIS FORMULATION ASSUMES 1:1 RATIO BETWEEN MASS LOSS AND
1916 ! LOSS OF NUMBER CONCENTRATION
1918 ! IF (PCC(K).LT.0.) THEN
1919 ! DUM = PCC(K)*DT/QC3D(K)
1920 ! DUM = MAX(-1.,DUM)
1921 ! NSUBC(K) = DUM*NC3D(K)/DT
1922 ! END IF
1924 ! UPDATE TENDENCIES
1926 ! NC3DTEN(K) = NC3DTEN(K)+NSUBC(K)
1928 !.....................................................................
1929 !.....................................................................
1930 ELSE ! TEMPERATURE < 273.15
1932 !......................................................................
1933 !HM ADD, ALLOW FOR CONSTANT DROPLET NUMBER
1934 ! INUM = 0, PREDICT DROPLET NUMBER
1935 ! INUM = 1, SET CONSTANT DROPLET NUMBER
1937 IF (iinum.EQ.1) THEN
1938 ! CONVERT NDCNST FROM CM-3 TO KG-1
1939 NC3D(K)=NDCNST*1.E6/RHO(K)
1940 END IF
1942 ! CALCULATE SIZE DISTRIBUTION PARAMETERS
1943 ! MAKE SURE NUMBER CONCENTRATIONS AREN'T NEGATIVE
1945 NI3D(K) = MAX(0.,NI3D(K))
1946 NS3D(K) = MAX(0.,NS3D(K))
1947 NC3D(K) = MAX(0.,NC3D(K))
1948 NR3D(K) = MAX(0.,NR3D(K))
1949 NG3D(K) = MAX(0.,NG3D(K))
1951 !......................................................................
1952 ! CLOUD ICE
1954 IF (QI3D(K).GE.QSMALL) THEN
1955 LAMI(K) = (CONS12* &
1956 NI3D(K)/QI3D(K))**(1./DI)
1957 N0I(K) = NI3D(K)*LAMI(K)
1959 ! CHECK FOR SLOPE
1961 ! ADJUST VARS
1963 IF (LAMI(K).LT.LAMMINI) THEN
1965 LAMI(K) = LAMMINI
1967 N0I(K) = LAMI(K)**4*QI3D(K)/CONS12
1969 NI3D(K) = N0I(K)/LAMI(K)
1970 ELSE IF (LAMI(K).GT.LAMMAXI) THEN
1971 LAMI(K) = LAMMAXI
1972 N0I(K) = LAMI(K)**4*QI3D(K)/CONS12
1974 NI3D(K) = N0I(K)/LAMI(K)
1975 END IF
1976 END IF
1978 !......................................................................
1979 ! RAIN
1981 IF (QR3D(K).GE.QSMALL) THEN
1982 LAMR(K) = (PI*RHOW*NR3D(K)/QR3D(K))**(1./3.)
1983 N0RR(K) = NR3D(K)*LAMR(K)
1985 ! CHECK FOR SLOPE
1987 ! ADJUST VARS
1989 IF (LAMR(K).LT.LAMMINR) THEN
1991 LAMR(K) = LAMMINR
1993 N0RR(K) = LAMR(K)**4*QR3D(K)/(PI*RHOW)
1995 NR3D(K) = N0RR(K)/LAMR(K)
1996 ELSE IF (LAMR(K).GT.LAMMAXR) THEN
1997 LAMR(K) = LAMMAXR
1998 N0RR(K) = LAMR(K)**4*QR3D(K)/(PI*RHOW)
2000 NR3D(K) = N0RR(K)/LAMR(K)
2001 END IF
2002 END IF
2004 !......................................................................
2005 ! CLOUD DROPLETS
2007 ! MARTIN ET AL. (1994) FORMULA FOR PGAM
2009 IF (QC3D(K).GE.QSMALL) THEN
2011 DUM = PRES(K)/(287.15*T3D(K))
2012 PGAM(K)=0.0005714*(NC3D(K)/1.E6*DUM)+0.2714
2013 PGAM(K)=1./(PGAM(K)**2)-1.
2014 PGAM(K)=MAX(PGAM(K),2.)
2015 PGAM(K)=MIN(PGAM(K),10.)
2017 ! CALCULATE LAMC
2019 LAMC(K) = (CONS26*NC3D(K)*GAMMA(PGAM(K)+4.)/ &
2020 (QC3D(K)*GAMMA(PGAM(K)+1.)))**(1./3.)
2022 ! LAMMIN, 60 MICRON DIAMETER
2023 ! LAMMAX, 1 MICRON
2025 LAMMIN = (PGAM(K)+1.)/60.E-6
2026 LAMMAX = (PGAM(K)+1.)/1.E-6
2028 IF (LAMC(K).LT.LAMMIN) THEN
2029 LAMC(K) = LAMMIN
2031 NC3D(K) = EXP(3.*LOG(LAMC(K))+LOG(QC3D(K))+ &
2032 LOG(GAMMA(PGAM(K)+1.))-LOG(GAMMA(PGAM(K)+4.)))/CONS26
2033 ELSE IF (LAMC(K).GT.LAMMAX) THEN
2034 LAMC(K) = LAMMAX
2035 NC3D(K) = EXP(3.*LOG(LAMC(K))+LOG(QC3D(K))+ &
2036 LOG(GAMMA(PGAM(K)+1.))-LOG(GAMMA(PGAM(K)+4.)))/CONS26
2038 END IF
2040 ! TO CALCULATE DROPLET FREEZING
2042 CDIST1(K) = NC3D(K)/GAMMA(PGAM(K)+1.)
2044 END IF
2046 !......................................................................
2047 ! SNOW
2049 IF (QNI3D(K).GE.QSMALL) THEN
2050 LAMS(K) = (CONS1*NS3D(K)/QNI3D(K))**(1./DS)
2051 N0S(K) = NS3D(K)*LAMS(K)
2053 ! CHECK FOR SLOPE
2055 ! ADJUST VARS
2057 IF (LAMS(K).LT.LAMMINS) THEN
2058 LAMS(K) = LAMMINS
2059 N0S(K) = LAMS(K)**4*QNI3D(K)/CONS1
2061 NS3D(K) = N0S(K)/LAMS(K)
2063 ELSE IF (LAMS(K).GT.LAMMAXS) THEN
2065 LAMS(K) = LAMMAXS
2066 N0S(K) = LAMS(K)**4*QNI3D(K)/CONS1
2068 NS3D(K) = N0S(K)/LAMS(K)
2069 END IF
2070 END IF
2072 !......................................................................
2073 ! GRAUPEL
2075 IF (QG3D(K).GE.QSMALL) THEN
2076 LAMG(K) = (CONS2*NG3D(K)/QG3D(K))**(1./DG)
2077 N0G(K) = NG3D(K)*LAMG(K)
2079 ! CHECK FOR SLOPE
2081 ! ADJUST VARS
2083 IF (LAMG(K).LT.LAMMING) THEN
2084 LAMG(K) = LAMMING
2085 N0G(K) = LAMG(K)**4*QG3D(K)/CONS2
2087 NG3D(K) = N0G(K)/LAMG(K)
2089 ELSE IF (LAMG(K).GT.LAMMAXG) THEN
2091 LAMG(K) = LAMMAXG
2092 N0G(K) = LAMG(K)**4*QG3D(K)/CONS2
2094 NG3D(K) = N0G(K)/LAMG(K)
2095 END IF
2096 END IF
2098 !.....................................................................
2099 ! ZERO OUT PROCESS RATES
2101 MNUCCC(K) = 0.
2102 NNUCCC(K) = 0.
2103 PRC(K) = 0.
2104 NPRC(K) = 0.
2105 NPRC1(K) = 0.
2106 NSAGG(K) = 0.
2107 PSACWS(K) = 0.
2108 NPSACWS(K) = 0.
2109 PSACWI(K) = 0.
2110 NPSACWI(K) = 0.
2111 PRACS(K) = 0.
2112 NPRACS(K) = 0.
2113 NMULTS(K) = 0.
2114 QMULTS(K) = 0.
2115 NMULTR(K) = 0.
2116 QMULTR(K) = 0.
2117 NMULTG(K) = 0.
2118 QMULTG(K) = 0.
2119 NMULTRG(K) = 0.
2120 QMULTRG(K) = 0.
2121 MNUCCR(K) = 0.
2122 NNUCCR(K) = 0.
2123 PRA(K) = 0.
2124 NPRA(K) = 0.
2125 NRAGG(K) = 0.
2126 PRCI(K) = 0.
2127 NPRCI(K) = 0.
2128 PRAI(K) = 0.
2129 NPRAI(K) = 0.
2130 NNUCCD(K) = 0.
2131 MNUCCD(K) = 0.
2132 PCC(K) = 0.
2133 PRE(K) = 0.
2134 PRD(K) = 0.
2135 PRDS(K) = 0.
2136 EPRD(K) = 0.
2137 EPRDS(K) = 0.
2138 NSUBC(K) = 0.
2139 NSUBI(K) = 0.
2140 NSUBS(K) = 0.
2141 NSUBR(K) = 0.
2142 PIACR(K) = 0.
2143 NIACR(K) = 0.
2144 PRACI(K) = 0.
2145 PIACRS(K) = 0.
2146 NIACRS(K) = 0.
2147 PRACIS(K) = 0.
2148 ! HM: ADD GRAUPEL PROCESSES
2149 PRACG(K) = 0.
2150 PSACR(K) = 0.
2151 PSACWG(K) = 0.
2152 PGSACW(K) = 0.
2153 PGRACS(K) = 0.
2154 PRDG(K) = 0.
2155 EPRDG(K) = 0.
2156 NPRACG(K) = 0.
2157 NPSACWG(K) = 0.
2158 NSCNG(K) = 0.
2159 NGRACS(K) = 0.
2160 NSUBG(K) = 0.
2162 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
2163 ! CALCULATION OF MICROPHYSICAL PROCESS RATES
2164 ! ACCRETION/AUTOCONVERSION/FREEZING/MELTING/COAG.
2165 !.......................................................................
2166 ! FREEZING OF CLOUD DROPLETS
2167 ! ONLY ALLOWED BELOW -4 C
2168 IF (QC3D(K).GE.QSMALL .AND. T3D(K).LT.269.15) THEN
2170 ! NUMBER OF CONTACT NUCLEI (M^-3) FROM MEYERS ET AL., 1992
2171 ! FACTOR OF 1000 IS TO CONVERT FROM L^-1 TO M^-3
2173 ! MEYERS CURVE
2175 NACNT = EXP(-2.80+0.262*(273.15-T3D(K)))*1000.
2177 ! COOPER CURVE
2178 ! NACNT = 5.*EXP(0.304*(273.15-T3D(K)))
2180 ! FLECTHER
2181 ! NACNT = 0.01*EXP(0.6*(273.15-T3D(K)))
2183 ! CONTACT FREEZING
2185 ! MEAN FREE PATH
2187 DUM = 7.37*T3D(K)/(288.*10.*PRES(K))/100.
2189 ! EFFECTIVE DIFFUSIVITY OF CONTACT NUCLEI
2190 ! BASED ON BROWNIAN DIFFUSION
2192 DAP(K) = CONS37*T3D(K)*(1.+DUM/RIN)/MU(K)
2194 MNUCCC(K) = CONS38*DAP(K)*NACNT*EXP(LOG(CDIST1(K))+ &
2195 LOG(GAMMA(PGAM(K)+5.))-4.*LOG(LAMC(K)))
2196 NNUCCC(K) = 2.*PI*DAP(K)*NACNT*CDIST1(K)* &
2197 GAMMA(PGAM(K)+2.)/ &
2198 LAMC(K)
2200 ! IMMERSION FREEZING (BIGG 1953)
2202 MNUCCC(K) = MNUCCC(K)+CONS39* &
2203 EXP(LOG(CDIST1(K))+LOG(GAMMA(7.+PGAM(K)))-6.*LOG(LAMC(K)))* &
2204 EXP(AIMM*(273.15-T3D(K)))
2206 NNUCCC(K) = NNUCCC(K)+ &
2207 CONS40*EXP(LOG(CDIST1(K))+LOG(GAMMA(PGAM(K)+4.))-3.*LOG(LAMC(K))) &
2208 *EXP(AIMM*(273.15-T3D(K)))
2210 ! PUT IN A CATCH HERE TO PREVENT DIVERGENCE BETWEEN NUMBER CONC. AND
2211 ! MIXING RATIO, SINCE STRICT CONSERVATION NOT CHECKED FOR NUMBER CONC
2213 NNUCCC(K) = MIN(NNUCCC(K),NC3D(K)/DT)
2215 END IF
2217 !.................................................................
2218 !.......................................................................
2219 ! AUTOCONVERSION OF CLOUD LIQUID WATER TO RAIN
2220 ! FORMULA FROM BEHENG (1994)
2221 ! USING NUMERICAL SIMULATION OF STOCHASTIC COLLECTION EQUATION
2222 ! AND INITIAL CLOUD DROPLET SIZE DISTRIBUTION SPECIFIED
2223 ! AS A GAMMA DISTRIBUTION
2225 ! USE MINIMUM VALUE OF 1.E-6 TO PREVENT FLOATING POINT ERROR
2227 IF (QC3D(K).GE.1.E-6) THEN
2229 ! HM ADD 12/13/06, REPLACE WITH NEWER FORMULA
2230 ! FROM KHAIROUTDINOV AND KOGAN 2000, MWR
2232 PRC(K)=1350.*QC3D(K)**2.47* &
2233 (NC3D(K)/1.e6*RHO(K))**(-1.79)
2235 ! note: nprc1 is change in Nr,
2236 ! nprc is change in Nc
2238 NPRC1(K) = PRC(K)/CONS29
2239 NPRC(K) = PRC(K)/(QC3D(K)/NC3D(K))
2241 ! hm bug fix 3/20/12
2242 NPRC(K) = MIN(NPRC(K),NC3D(K)/DT)
2243 NPRC1(K) = MIN(NPRC1(K),NPRC(K))
2245 END IF
2247 !.......................................................................
2248 ! SELF-COLLECTION OF DROPLET NOT INCLUDED IN KK2000 SCHEME
2250 ! SNOW AGGREGATION FROM PASSARELLI, 1978, USED BY REISNER, 1998
2251 ! THIS IS HARD-WIRED FOR BS = 0.4 FOR NOW
2253 IF (QNI3D(K).GE.1.E-8) THEN
2254 NSAGG(K) = CONS15*ASN(K)*RHO(K)** &
2255 ((2.+BS)/3.)*QNI3D(K)**((2.+BS)/3.)* &
2256 (NS3D(K)*RHO(K))**((4.-BS)/3.)/ &
2257 (RHO(K))
2258 END IF
2260 !.......................................................................
2261 ! ACCRETION OF CLOUD DROPLETS ONTO SNOW/GRAUPEL
2262 ! HERE USE CONTINUOUS COLLECTION EQUATION WITH
2263 ! SIMPLE GRAVITATIONAL COLLECTION KERNEL IGNORING
2265 ! SNOW
2267 IF (QNI3D(K).GE.1.E-8 .AND. QC3D(K).GE.QSMALL) THEN
2269 PSACWS(K) = CONS13*ASN(K)*QC3D(K)*RHO(K)* &
2270 N0S(K)/ &
2271 LAMS(K)**(BS+3.)
2272 NPSACWS(K) = CONS13*ASN(K)*NC3D(K)*RHO(K)* &
2273 N0S(K)/ &
2274 LAMS(K)**(BS+3.)
2276 END IF
2278 !............................................................................
2279 ! COLLECTION OF CLOUD WATER BY GRAUPEL
2281 IF (QG3D(K).GE.1.E-8 .AND. QC3D(K).GE.QSMALL) THEN
2283 PSACWG(K) = CONS14*AGN(K)*QC3D(K)*RHO(K)* &
2284 N0G(K)/ &
2285 LAMG(K)**(BG+3.)
2286 NPSACWG(K) = CONS14*AGN(K)*NC3D(K)*RHO(K)* &
2287 N0G(K)/ &
2288 LAMG(K)**(BG+3.)
2289 END IF
2291 !.......................................................................
2292 ! HM, ADD 12/13/06
2293 ! CLOUD ICE COLLECTING DROPLETS, ASSUME THAT CLOUD ICE MEAN DIAM > 100 MICRON
2294 ! BEFORE RIMING CAN OCCUR
2295 ! ASSUME THAT RIME COLLECTED ON CLOUD ICE DOES NOT LEAD
2296 ! TO HALLET-MOSSOP SPLINTERING
2298 IF (QI3D(K).GE.1.E-8 .AND. QC3D(K).GE.QSMALL) THEN
2300 ! PUT IN SIZE DEPENDENT COLLECTION EFFICIENCY BASED ON STOKES LAW
2301 ! FROM THOMPSON ET AL. 2004, MWR
2303 IF (1./LAMI(K).GE.100.E-6) THEN
2305 PSACWI(K) = CONS16*AIN(K)*QC3D(K)*RHO(K)* &
2306 N0I(K)/ &
2307 LAMI(K)**(BI+3.)
2308 NPSACWI(K) = CONS16*AIN(K)*NC3D(K)*RHO(K)* &
2309 N0I(K)/ &
2310 LAMI(K)**(BI+3.)
2311 END IF
2312 END IF
2314 !.......................................................................
2315 ! ACCRETION OF RAIN WATER BY SNOW
2316 ! FORMULA FROM IKAWA AND SAITO, 1991, USED BY REISNER ET AL, 1998
2318 IF (QR3D(K).GE.1.E-8.AND.QNI3D(K).GE.1.E-8) THEN
2320 UMS = ASN(K)*CONS3/(LAMS(K)**BS)
2321 UMR = ARN(K)*CONS4/(LAMR(K)**BR)
2322 UNS = ASN(K)*CONS5/LAMS(K)**BS
2323 UNR = ARN(K)*CONS6/LAMR(K)**BR
2325 ! SET REASLISTIC LIMITS ON FALLSPEEDS
2327 ! bug fix, 10/08/09
2328 dum=(rhosu/rho(k))**0.54
2329 UMS=MIN(UMS,1.2*dum)
2330 UNS=MIN(UNS,1.2*dum)
2331 UMR=MIN(UMR,9.1*dum)
2332 UNR=MIN(UNR,9.1*dum)
2334 PRACS(K) = CONS41*(((1.2*UMR-0.95*UMS)**2+ &
2335 0.08*UMS*UMR)**0.5*RHO(K)* &
2336 N0RR(K)*N0S(K)/LAMR(K)**3* &
2337 (5./(LAMR(K)**3*LAMS(K))+ &
2338 2./(LAMR(K)**2*LAMS(K)**2)+ &
2339 0.5/(LAMR(k)*LAMS(k)**3)))
2341 NPRACS(K) = CONS32*RHO(K)*(1.7*(UNR-UNS)**2+ &
2342 0.3*UNR*UNS)**0.5*N0RR(K)*N0S(K)* &
2343 (1./(LAMR(K)**3*LAMS(K))+ &
2344 1./(LAMR(K)**2*LAMS(K)**2)+ &
2345 1./(LAMR(K)*LAMS(K)**3))
2347 ! MAKE SURE PRACS DOESN'T EXCEED TOTAL RAIN MIXING RATIO
2348 ! AS THIS MAY OTHERWISE RESULT IN TOO MUCH TRANSFER OF WATER DURING
2349 ! RIME-SPLINTERING
2351 PRACS(K) = MIN(PRACS(K),QR3D(K)/DT)
2353 ! COLLECTION OF SNOW BY RAIN - NEEDED FOR GRAUPEL CONVERSION CALCULATIONS
2354 ! ONLY CALCULATE IF SNOW AND RAIN MIXING RATIOS EXCEED 0.1 G/KG
2356 ! ASSUME COLLECTION OF SNOW BY RAIN PRODUCES GRAUPEL NOT HAIL
2358 ! HM MODIFY FOR WRFV3.1
2359 ! IF (IHAIL.EQ.0) THEN
2360 IF (QNI3D(K).GE.0.1E-3.AND.QR3D(K).GE.0.1E-3) THEN
2361 PSACR(K) = CONS31*(((1.2*UMR-0.95*UMS)**2+ &
2362 0.08*UMS*UMR)**0.5*RHO(K)* &
2363 N0RR(K)*N0S(K)/LAMS(K)**3* &
2364 (5./(LAMS(K)**3*LAMR(K))+ &
2365 2./(LAMS(K)**2*LAMR(K)**2)+ &
2366 0.5/(LAMS(K)*LAMR(K)**3)))
2367 END IF
2368 ! END IF
2370 END IF
2372 !.......................................................................
2374 ! COLLECTION OF RAINWATER BY GRAUPEL, FROM IKAWA AND SAITO 1990,
2375 ! USED BY REISNER ET AL 1998
2376 IF (QR3D(K).GE.1.E-8.AND.QG3D(K).GE.1.E-8) THEN
2378 UMG = AGN(K)*CONS7/(LAMG(K)**BG)
2379 UMR = ARN(K)*CONS4/(LAMR(K)**BR)
2380 UNG = AGN(K)*CONS8/LAMG(K)**BG
2381 UNR = ARN(K)*CONS6/LAMR(K)**BR
2383 ! SET REASLISTIC LIMITS ON FALLSPEEDS
2384 ! bug fix, 10/08/09
2385 dum=(rhosu/rho(k))**0.54
2386 UMG=MIN(UMG,20.*dum)
2387 UNG=MIN(UNG,20.*dum)
2388 UMR=MIN(UMR,9.1*dum)
2389 UNR=MIN(UNR,9.1*dum)
2391 PRACG(K) = CONS41*(((1.2*UMR-0.95*UMG)**2+ &
2392 0.08*UMG*UMR)**0.5*RHO(K)* &
2393 N0RR(K)*N0G(K)/LAMR(K)**3* &
2394 (5./(LAMR(K)**3*LAMG(K))+ &
2395 2./(LAMR(K)**2*LAMG(K)**2)+ &
2396 0.5/(LAMR(k)*LAMG(k)**3)))
2398 NPRACG(K) = CONS32*RHO(K)*(1.7*(UNR-UNG)**2+ &
2399 0.3*UNR*UNG)**0.5*N0RR(K)*N0G(K)* &
2400 (1./(LAMR(K)**3*LAMG(K))+ &
2401 1./(LAMR(K)**2*LAMG(K)**2)+ &
2402 1./(LAMR(K)*LAMG(K)**3))
2404 ! MAKE SURE PRACG DOESN'T EXCEED TOTAL RAIN MIXING RATIO
2405 ! AS THIS MAY OTHERWISE RESULT IN TOO MUCH TRANSFER OF WATER DURING
2406 ! RIME-SPLINTERING
2408 PRACG(K) = MIN(PRACG(K),QR3D(K)/DT)
2410 END IF
2412 !.......................................................................
2413 ! RIME-SPLINTERING - SNOW
2414 ! HALLET-MOSSOP (1974)
2415 ! NUMBER OF SPLINTERS FORMED IS BASED ON MASS OF RIMED WATER
2417 ! DUM1 = MASS OF INDIVIDUAL SPLINTERS
2419 ! HM ADD THRESHOLD SNOW AND DROPLET MIXING RATIO FOR RIME-SPLINTERING
2420 ! TO LIMIT RIME-SPLINTERING IN STRATIFORM CLOUDS
2421 ! THESE THRESHOLDS CORRESPOND WITH GRAUPEL THRESHOLDS IN RH 1984
2423 !v1.4
2424 IF (QNI3D(K).GE.0.1E-3) THEN
2425 IF (QC3D(K).GE.0.5E-3.OR.QR3D(K).GE.0.1E-3) THEN
2426 IF (PSACWS(K).GT.0..OR.PRACS(K).GT.0.) THEN
2427 IF (T3D(K).LT.270.16 .AND. T3D(K).GT.265.16) THEN
2429 IF (T3D(K).GT.270.16) THEN
2430 FMULT = 0.
2431 ELSE IF (T3D(K).LE.270.16.AND.T3D(K).GT.268.16) THEN
2432 FMULT = (270.16-T3D(K))/2.
2433 ELSE IF (T3D(K).GE.265.16.AND.T3D(K).LE.268.16) THEN
2434 FMULT = (T3D(K)-265.16)/3.
2435 ELSE IF (T3D(K).LT.265.16) THEN
2436 FMULT = 0.
2437 END IF
2439 ! 1000 IS TO CONVERT FROM KG TO G
2441 ! SPLINTERING FROM DROPLETS ACCRETED ONTO SNOW
2443 IF (PSACWS(K).GT.0.) THEN
2444 NMULTS(K) = 35.E4*PSACWS(K)*FMULT*1000.
2445 QMULTS(K) = NMULTS(K)*MMULT
2447 ! CONSTRAIN SO THAT TRANSFER OF MASS FROM SNOW TO ICE CANNOT BE MORE MASS
2448 ! THAN WAS RIMED ONTO SNOW
2450 QMULTS(K) = MIN(QMULTS(K),PSACWS(K))
2451 PSACWS(K) = PSACWS(K)-QMULTS(K)
2453 END IF
2455 ! RIMING AND SPLINTERING FROM ACCRETED RAINDROPS
2457 IF (PRACS(K).GT.0.) THEN
2458 NMULTR(K) = 35.E4*PRACS(K)*FMULT*1000.
2459 QMULTR(K) = NMULTR(K)*MMULT
2461 ! CONSTRAIN SO THAT TRANSFER OF MASS FROM SNOW TO ICE CANNOT BE MORE MASS
2462 ! THAN WAS RIMED ONTO SNOW
2464 QMULTR(K) = MIN(QMULTR(K),PRACS(K))
2466 PRACS(K) = PRACS(K)-QMULTR(K)
2468 END IF
2470 END IF
2471 END IF
2472 END IF
2473 END IF
2475 !.......................................................................
2476 ! RIME-SPLINTERING - GRAUPEL
2477 ! HALLET-MOSSOP (1974)
2478 ! NUMBER OF SPLINTERS FORMED IS BASED ON MASS OF RIMED WATER
2480 ! DUM1 = MASS OF INDIVIDUAL SPLINTERS
2482 ! HM ADD THRESHOLD SNOW MIXING RATIO FOR RIME-SPLINTERING
2483 ! TO LIMIT RIME-SPLINTERING IN STRATIFORM CLOUDS
2485 ! ONLY CALCULATE FOR GRAUPEL NOT HAIL
2486 ! IF (IHAIL.EQ.0) THEN
2487 ! v1.4
2488 IF (QG3D(K).GE.0.1E-3) THEN
2489 IF (QC3D(K).GE.0.5E-3.OR.QR3D(K).GE.0.1E-3) THEN
2490 IF (PSACWG(K).GT.0..OR.PRACG(K).GT.0.) THEN
2491 IF (T3D(K).LT.270.16 .AND. T3D(K).GT.265.16) THEN
2493 IF (T3D(K).GT.270.16) THEN
2494 FMULT = 0.
2495 ELSE IF (T3D(K).LE.270.16.AND.T3D(K).GT.268.16) THEN
2496 FMULT = (270.16-T3D(K))/2.
2497 ELSE IF (T3D(K).GE.265.16.AND.T3D(K).LE.268.16) THEN
2498 FMULT = (T3D(K)-265.16)/3.
2499 ELSE IF (T3D(K).LT.265.16) THEN
2500 FMULT = 0.
2501 END IF
2503 ! 1000 IS TO CONVERT FROM KG TO G
2505 ! SPLINTERING FROM DROPLETS ACCRETED ONTO GRAUPEL
2507 IF (PSACWG(K).GT.0.) THEN
2508 NMULTG(K) = 35.E4*PSACWG(K)*FMULT*1000.
2509 QMULTG(K) = NMULTG(K)*MMULT
2511 ! CONSTRAIN SO THAT TRANSFER OF MASS FROM GRAUPEL TO ICE CANNOT BE MORE MASS
2512 ! THAN WAS RIMED ONTO GRAUPEL
2514 QMULTG(K) = MIN(QMULTG(K),PSACWG(K))
2515 PSACWG(K) = PSACWG(K)-QMULTG(K)
2517 END IF
2519 ! RIMING AND SPLINTERING FROM ACCRETED RAINDROPS
2521 IF (PRACG(K).GT.0.) THEN
2522 NMULTRG(K) = 35.E4*PRACG(K)*FMULT*1000.
2523 QMULTRG(K) = NMULTRG(K)*MMULT
2525 ! CONSTRAIN SO THAT TRANSFER OF MASS FROM GRAUPEL TO ICE CANNOT BE MORE MASS
2526 ! THAN WAS RIMED ONTO GRAUPEL
2528 QMULTRG(K) = MIN(QMULTRG(K),PRACG(K))
2529 PRACG(K) = PRACG(K)-QMULTRG(K)
2531 END IF
2532 END IF
2533 END IF
2534 END IF
2535 END IF
2536 ! END IF
2538 !........................................................................
2539 ! CONVERSION OF RIMED CLOUD WATER ONTO SNOW TO GRAUPEL
2540 ! ASSUME CONVERTED SNOW FORMS GRAUPEL NOT HAIL
2541 ! HAIL ASSUMED TO ONLY FORM BY FREEZING OF RAIN
2542 ! OR COLLISIONS OF RAIN WITH CLOUD ICE
2544 ! IF (IHAIL.EQ.0) THEN
2545 IF (PSACWS(K).GT.0.) THEN
2546 ! ONLY ALLOW CONVERSION IF QNI > 0.1 AND QC > 0.5 G/KG FOLLOWING RUTLEDGE AND HOBBS (1984)
2547 IF (QNI3D(K).GE.0.1E-3.AND.QC3D(K).GE.0.5E-3) THEN
2549 ! PORTION OF RIMING CONVERTED TO GRAUPEL (REISNER ET AL. 1998, ORIGINALLY IS1991)
2550 PGSACW(K) = MIN(PSACWS(K),CONS17*DT*N0S(K)*QC3D(K)*QC3D(K)* &
2551 ASN(K)*ASN(K)/ &
2552 (RHO(K)*LAMS(K)**(2.*BS+2.)))
2554 ! MIX RAT CONVERTED INTO GRAUPEL AS EMBRYO (REISNER ET AL. 1998, ORIG M1990)
2555 DUM = MAX(RHOSN/(RHOG-RHOSN)*PGSACW(K),0.)
2557 ! NUMBER CONCENTRAITON OF EMBRYO GRAUPEL FROM RIMING OF SNOW
2558 NSCNG(K) = DUM/MG0*RHO(K)
2559 ! LIMIT MAX NUMBER CONVERTED TO SNOW NUMBER
2560 NSCNG(K) = MIN(NSCNG(K),NS3D(K)/DT)
2562 ! PORTION OF RIMING LEFT FOR SNOW
2563 PSACWS(K) = PSACWS(K) - PGSACW(K)
2564 END IF
2565 END IF
2567 ! CONVERSION OF RIMED RAINWATER ONTO SNOW CONVERTED TO GRAUPEL
2569 IF (PRACS(K).GT.0.) THEN
2570 ! ONLY ALLOW CONVERSION IF QNI > 0.1 AND QR > 0.1 G/KG FOLLOWING RUTLEDGE AND HOBBS (1984)
2571 IF (QNI3D(K).GE.0.1E-3.AND.QR3D(K).GE.0.1E-3) THEN
2572 ! PORTION OF COLLECTED RAINWATER CONVERTED TO GRAUPEL (REISNER ET AL. 1998)
2573 DUM = CONS18*(4./LAMS(K))**3*(4./LAMS(K))**3 &
2574 /(CONS18*(4./LAMS(K))**3*(4./LAMS(K))**3+ &
2575 CONS19*(4./LAMR(K))**3*(4./LAMR(K))**3)
2576 DUM=MIN(DUM,1.)
2577 DUM=MAX(DUM,0.)
2578 PGRACS(K) = (1.-DUM)*PRACS(K)
2579 NGRACS(K) = (1.-DUM)*NPRACS(K)
2580 ! LIMIT MAX NUMBER CONVERTED TO MIN OF EITHER RAIN OR SNOW NUMBER CONCENTRATION
2581 NGRACS(K) = MIN(NGRACS(K),NR3D(K)/DT)
2582 NGRACS(K) = MIN(NGRACS(K),NS3D(K)/DT)
2584 ! AMOUNT LEFT FOR SNOW PRODUCTION
2585 PRACS(K) = PRACS(K) - PGRACS(K)
2586 NPRACS(K) = NPRACS(K) - NGRACS(K)
2587 ! CONVERSION TO GRAUPEL DUE TO COLLECTION OF SNOW BY RAIN
2588 PSACR(K)=PSACR(K)*(1.-DUM)
2589 END IF
2590 END IF
2591 ! END IF
2593 !.......................................................................
2594 ! FREEZING OF RAIN DROPS
2595 ! FREEZING ALLOWED BELOW -4 C
2597 IF (T3D(K).LT.269.15.AND.QR3D(K).GE.QSMALL) THEN
2599 ! IMMERSION FREEZING (BIGG 1953)
2600 MNUCCR(K) = CONS20*NR3D(K)*EXP(AIMM*(273.15-T3D(K)))/LAMR(K)**3 &
2601 /LAMR(K)**3
2603 NNUCCR(K) = PI*NR3D(K)*BIMM*EXP(AIMM*(273.15-T3D(K)))/LAMR(K)**3
2605 ! PREVENT DIVERGENCE BETWEEN MIXING RATIO AND NUMBER CONC
2606 NNUCCR(K) = MIN(NNUCCR(K),NR3D(K)/DT)
2608 END IF
2610 !.......................................................................
2611 ! ACCRETION OF CLOUD LIQUID WATER BY RAIN
2612 ! CONTINUOUS COLLECTION EQUATION WITH
2613 ! GRAVITATIONAL COLLECTION KERNEL, DROPLET FALL SPEED NEGLECTED
2615 IF (QR3D(K).GE.1.E-8 .AND. QC3D(K).GE.1.E-8) THEN
2617 ! 12/13/06 HM ADD, REPLACE WITH NEWER FORMULA FROM
2618 ! KHAIROUTDINOV AND KOGAN 2000, MWR
2620 DUM=(QC3D(K)*QR3D(K))
2621 PRA(K) = 67.*(DUM)**1.15
2622 NPRA(K) = PRA(K)/(QC3D(K)/NC3D(K))
2624 END IF
2625 !.......................................................................
2626 ! SELF-COLLECTION OF RAIN DROPS
2627 ! FROM BEHENG(1994)
2628 ! FROM NUMERICAL SIMULATION OF THE STOCHASTIC COLLECTION EQUATION
2629 ! AS DESCRINED ABOVE FOR AUTOCONVERSION
2631 IF (QR3D(K).GE.1.E-8) THEN
2632 ! include breakup add 10/09/09
2633 dum1=300.e-6
2634 if (1./lamr(k).lt.dum1) then
2635 dum=1.
2636 else if (1./lamr(k).ge.dum1) then
2637 dum=2.-exp(2300.*(1./lamr(k)-dum1))
2638 end if
2639 ! NRAGG(K) = -8.*NR3D(K)*QR3D(K)*RHO(K)
2640 NRAGG(K) = -5.78*dum*NR3D(K)*QR3D(K)*RHO(K)
2641 END IF
2643 !.......................................................................
2644 ! AUTOCONVERSION OF CLOUD ICE TO SNOW
2645 ! FOLLOWING HARRINGTON ET AL. (1995) WITH MODIFICATION
2646 ! HERE IT IS ASSUMED THAT AUTOCONVERSION CAN ONLY OCCUR WHEN THE
2647 ! ICE IS GROWING, I.E. IN CONDITIONS OF ICE SUPERSATURATION
2649 IF (QI3D(K).GE.1.E-8 .AND.QVQVSI(K).GE.1.) THEN
2651 ! COFFI = 2./LAMI(K)
2652 ! IF (COFFI.GE.DCS) THEN
2653 NPRCI(K) = CONS21*(QV3D(K)-QVI(K))*RHO(K) &
2654 *N0I(K)*EXP(-LAMI(K)*DCS)*DV(K)/ABI(K)
2655 PRCI(K) = CONS22*NPRCI(K)
2656 NPRCI(K) = MIN(NPRCI(K),NI3D(K)/DT)
2658 ! END IF
2659 END IF
2661 !.......................................................................
2662 ! ACCRETION OF CLOUD ICE BY SNOW
2663 ! FOR THIS CALCULATION, IT IS ASSUMED THAT THE VS >> VI
2664 ! AND DS >> DI FOR CONTINUOUS COLLECTION
2666 IF (QNI3D(K).GE.1.E-8 .AND. QI3D(K).GE.QSMALL) THEN
2667 PRAI(K) = CONS23*ASN(K)*QI3D(K)*RHO(K)*N0S(K)/ &
2668 LAMS(K)**(BS+3.)
2669 NPRAI(K) = CONS23*ASN(K)*NI3D(K)* &
2670 RHO(K)*N0S(K)/ &
2671 LAMS(K)**(BS+3.)
2672 NPRAI(K)=MIN(NPRAI(K),NI3D(K)/DT)
2673 END IF
2675 !.......................................................................
2676 ! HM, ADD 12/13/06, COLLISION OF RAIN AND ICE TO PRODUCE SNOW OR GRAUPEL
2677 ! FOLLOWS REISNER ET AL. 1998
2678 ! ASSUMED FALLSPEED AND SIZE OF ICE CRYSTAL << THAN FOR RAIN
2680 IF (QR3D(K).GE.1.E-8.AND.QI3D(K).GE.1.E-8.AND.T3D(K).LE.273.15) THEN
2682 ! ALLOW GRAUPEL FORMATION FROM RAIN-ICE COLLISIONS ONLY IF RAIN MIXING RATIO > 0.1 G/KG,
2683 ! OTHERWISE ADD TO SNOW
2685 IF (QR3D(K).GE.0.1E-3) THEN
2686 NIACR(K)=CONS24*NI3D(K)*N0RR(K)*ARN(K) &
2687 /LAMR(K)**(BR+3.)*RHO(K)
2688 PIACR(K)=CONS25*NI3D(K)*N0RR(K)*ARN(K) &
2689 /LAMR(K)**(BR+3.)/LAMR(K)**3*RHO(K)
2690 PRACI(K)=CONS24*QI3D(K)*N0RR(K)*ARN(K)/ &
2691 LAMR(K)**(BR+3.)*RHO(K)
2692 NIACR(K)=MIN(NIACR(K),NR3D(K)/DT)
2693 NIACR(K)=MIN(NIACR(K),NI3D(K)/DT)
2694 ELSE
2695 NIACRS(K)=CONS24*NI3D(K)*N0RR(K)*ARN(K) &
2696 /LAMR(K)**(BR+3.)*RHO(K)
2697 PIACRS(K)=CONS25*NI3D(K)*N0RR(K)*ARN(K) &
2698 /LAMR(K)**(BR+3.)/LAMR(K)**3*RHO(K)
2699 PRACIS(K)=CONS24*QI3D(K)*N0RR(K)*ARN(K)/ &
2700 LAMR(K)**(BR+3.)*RHO(K)
2701 NIACRS(K)=MIN(NIACRS(K),NR3D(K)/DT)
2702 NIACRS(K)=MIN(NIACRS(K),NI3D(K)/DT)
2703 END IF
2704 END IF
2706 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
2707 ! NUCLEATION OF CLOUD ICE FROM HOMOGENEOUS AND HETEROGENEOUS FREEZING ON AEROSOL
2709 IF (INUC.EQ.0) THEN
2711 ! add threshold according to Greg Thomspon
2713 if ((QVQVS(K).GE.0.999.and.T3D(K).le.265.15).or. &
2714 QVQVSI(K).ge.1.08) then
2716 ! hm, modify dec. 5, 2006, replace with cooper curve
2717 kc2 = 0.005*exp(0.304*(273.15-T3D(K)))*1000. ! convert from L-1 to m-3
2718 ! limit to 500 L-1
2719 kc2 = min(kc2,500.e3)
2720 kc2=MAX(kc2/rho(k),0.) ! convert to kg-1
2722 IF (KC2.GT.NI3D(K)+NS3D(K)+NG3D(K)) THEN
2723 NNUCCD(K) = (KC2-NI3D(K)-NS3D(K)-NG3D(K))/DT
2724 MNUCCD(K) = NNUCCD(K)*MI0
2725 END IF
2727 END IF
2729 ELSE IF (INUC.EQ.1) THEN
2731 IF (T3D(K).LT.273.15.AND.QVQVSI(K).GT.1.) THEN
2733 KC2 = 0.16*1000./RHO(K) ! CONVERT FROM L-1 TO KG-1
2734 IF (KC2.GT.NI3D(K)+NS3D(K)+NG3D(K)) THEN
2735 NNUCCD(K) = (KC2-NI3D(K)-NS3D(K)-NG3D(K))/DT
2736 MNUCCD(K) = NNUCCD(K)*MI0
2737 END IF
2738 END IF
2740 END IF
2742 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
2744 101 CONTINUE
2746 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
2747 ! CALCULATE EVAP/SUB/DEP TERMS FOR QI,QNI,QR
2749 ! NO VENTILATION FOR CLOUD ICE
2751 IF (QI3D(K).GE.QSMALL) THEN
2753 EPSI = 2.*PI*N0I(K)*RHO(K)*DV(K)/(LAMI(K)*LAMI(K))
2755 ELSE
2756 EPSI = 0.
2757 END IF
2759 IF (QNI3D(K).GE.QSMALL) THEN
2760 EPSS = 2.*PI*N0S(K)*RHO(K)*DV(K)* &
2761 (F1S/(LAMS(K)*LAMS(K))+ &
2762 F2S*(ASN(K)*RHO(K)/MU(K))**0.5* &
2763 SC(K)**(1./3.)*CONS10/ &
2764 (LAMS(K)**CONS35))
2765 ELSE
2766 EPSS = 0.
2767 END IF
2769 IF (QG3D(K).GE.QSMALL) THEN
2770 EPSG = 2.*PI*N0G(K)*RHO(K)*DV(K)* &
2771 (F1S/(LAMG(K)*LAMG(K))+ &
2772 F2S*(AGN(K)*RHO(K)/MU(K))**0.5* &
2773 SC(K)**(1./3.)*CONS11/ &
2774 (LAMG(K)**CONS36))
2777 ELSE
2778 EPSG = 0.
2779 END IF
2781 IF (QR3D(K).GE.QSMALL) THEN
2782 EPSR = 2.*PI*N0RR(K)*RHO(K)*DV(K)* &
2783 (F1R/(LAMR(K)*LAMR(K))+ &
2784 F2R*(ARN(K)*RHO(K)/MU(K))**0.5* &
2785 SC(K)**(1./3.)*CONS9/ &
2786 (LAMR(K)**CONS34))
2787 ELSE
2788 EPSR = 0.
2789 END IF
2791 ! ONLY INCLUDE REGION OF ICE SIZE DIST < DCS
2792 ! DUM IS FRACTION OF D*N(D) < DCS
2794 ! LOGIC BELOW FOLLOWS THAT OF HARRINGTON ET AL. 1995 (JAS)
2795 IF (QI3D(K).GE.QSMALL) THEN
2796 DUM=(1.-EXP(-LAMI(K)*DCS)*(1.+LAMI(K)*DCS))
2797 PRD(K) = EPSI*(QV3D(K)-QVI(K))/ABI(K)*DUM
2798 ELSE
2799 DUM=0.
2800 END IF
2801 ! ADD DEPOSITION IN TAIL OF ICE SIZE DIST TO SNOW IF SNOW IS PRESENT
2802 IF (QNI3D(K).GE.QSMALL) THEN
2803 PRDS(K) = EPSS*(QV3D(K)-QVI(K))/ABI(K)+ &
2804 EPSI*(QV3D(K)-QVI(K))/ABI(K)*(1.-DUM)
2805 ! OTHERWISE ADD TO CLOUD ICE
2806 ELSE
2807 PRD(K) = PRD(K)+EPSI*(QV3D(K)-QVI(K))/ABI(K)*(1.-DUM)
2808 END IF
2809 ! VAPOR DPEOSITION ON GRAUPEL
2810 PRDG(K) = EPSG*(QV3D(K)-QVI(K))/ABI(K)
2812 ! NO CONDENSATION ONTO RAIN, ONLY EVAP
2814 IF (QV3D(K).LT.QVS(K)) THEN
2815 PRE(K) = EPSR*(QV3D(K)-QVS(K))/AB(K)
2816 PRE(K) = MIN(PRE(K),0.)
2817 ELSE
2818 PRE(K) = 0.
2819 END IF
2821 ! MAKE SURE NOT PUSHED INTO ICE SUPERSAT/SUBSAT
2822 ! FORMULA FROM REISNER 2 SCHEME
2824 DUM = (QV3D(K)-QVI(K))/DT
2826 FUDGEF = 0.9999
2827 SUM_DEP = PRD(K)+PRDS(K)+MNUCCD(K)+PRDG(K)
2829 IF( (DUM.GT.0. .AND. SUM_DEP.GT.DUM*FUDGEF) .OR. &
2830 (DUM.LT.0. .AND. SUM_DEP.LT.DUM*FUDGEF) ) THEN
2831 MNUCCD(K) = FUDGEF*MNUCCD(K)*DUM/SUM_DEP
2832 PRD(K) = FUDGEF*PRD(K)*DUM/SUM_DEP
2833 PRDS(K) = FUDGEF*PRDS(K)*DUM/SUM_DEP
2834 PRDG(K) = FUDGEF*PRDG(K)*DUM/SUM_DEP
2835 ENDIF
2837 ! IF CLOUD ICE/SNOW/GRAUPEL VAP DEPOSITION IS NEG, THEN ASSIGN TO SUBLIMATION PROCESSES
2839 IF (PRD(K).LT.0.) THEN
2840 EPRD(K)=PRD(K)
2841 PRD(K)=0.
2842 END IF
2843 IF (PRDS(K).LT.0.) THEN
2844 EPRDS(K)=PRDS(K)
2845 PRDS(K)=0.
2846 END IF
2847 IF (PRDG(K).LT.0.) THEN
2848 EPRDG(K)=PRDG(K)
2849 PRDG(K)=0.
2850 END IF
2852 !.......................................................................
2853 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
2855 ! CONSERVATION OF WATER
2856 ! THIS IS ADOPTED LOOSELY FROM MM5 RESINER CODE. HOWEVER, HERE WE
2857 ! ONLY ADJUST PROCESSES THAT ARE NEGATIVE, RATHER THAN ALL PROCESSES.
2859 ! IF MIXING RATIOS LESS THAN QSMALL, THEN NO DEPLETION OF WATER
2860 ! THROUGH MICROPHYSICAL PROCESSES, SKIP CONSERVATION
2862 ! NOTE: CONSERVATION CHECK NOT APPLIED TO NUMBER CONCENTRATION SPECIES. ADDITIONAL CATCH
2863 ! BELOW WILL PREVENT NEGATIVE NUMBER CONCENTRATION
2864 ! FOR EACH MICROPHYSICAL PROCESS WHICH PROVIDES A SOURCE FOR NUMBER, THERE IS A CHECK
2865 ! TO MAKE SURE THAT CAN'T EXCEED TOTAL NUMBER OF DEPLETED SPECIES WITH THE TIME
2866 ! STEP
2868 !****SENSITIVITY - NO ICE
2870 IF (ILIQ.EQ.1) THEN
2871 MNUCCC(K)=0.
2872 NNUCCC(K)=0.
2873 MNUCCR(K)=0.
2874 NNUCCR(K)=0.
2875 MNUCCD(K)=0.
2876 NNUCCD(K)=0.
2877 END IF
2879 ! ****SENSITIVITY - NO GRAUPEL
2880 IF (IGRAUP.EQ.1) THEN
2881 PRACG(K) = 0.
2882 PSACR(K) = 0.
2883 PSACWG(K) = 0.
2884 PGSACW(K) = 0.
2885 PGRACS(K) = 0.
2886 PRDG(K) = 0.
2887 EPRDG(K) = 0.
2888 EVPMG(K) = 0.
2889 PGMLT(K) = 0.
2890 NPRACG(K) = 0.
2891 NPSACWG(K) = 0.
2892 NSCNG(K) = 0.
2893 NGRACS(K) = 0.
2894 NSUBG(K) = 0.
2895 NGMLTG(K) = 0.
2896 NGMLTR(K) = 0.
2897 ! fix 053011
2898 PIACRS(K)=PIACRS(K)+PIACR(K)
2899 PIACR(K) = 0.
2900 END IF
2902 ! CONSERVATION OF QC
2904 DUM = (PRC(K)+PRA(K)+MNUCCC(K)+PSACWS(K)+PSACWI(K)+QMULTS(K)+PSACWG(K)+PGSACW(K)+QMULTG(K))*DT
2906 IF (DUM.GT.QC3D(K).AND.QC3D(K).GE.QSMALL) THEN
2907 RATIO = QC3D(K)/DUM
2909 PRC(K) = PRC(K)*RATIO
2910 PRA(K) = PRA(K)*RATIO
2911 MNUCCC(K) = MNUCCC(K)*RATIO
2912 PSACWS(K) = PSACWS(K)*RATIO
2913 PSACWI(K) = PSACWI(K)*RATIO
2914 QMULTS(K) = QMULTS(K)*RATIO
2915 QMULTG(K) = QMULTG(K)*RATIO
2916 PSACWG(K) = PSACWG(K)*RATIO
2917 PGSACW(K) = PGSACW(K)*RATIO
2918 END IF
2920 ! CONSERVATION OF QI
2922 DUM = (-PRD(K)-MNUCCC(K)+PRCI(K)+PRAI(K)-QMULTS(K)-QMULTG(K)-QMULTR(K)-QMULTRG(K) &
2923 -MNUCCD(K)+PRACI(K)+PRACIS(K)-EPRD(K)-PSACWI(K))*DT
2925 IF (DUM.GT.QI3D(K).AND.QI3D(K).GE.QSMALL) THEN
2927 RATIO = (QI3D(K)/DT+PRD(K)+MNUCCC(K)+QMULTS(K)+QMULTG(K)+QMULTR(K)+QMULTRG(K)+ &
2928 MNUCCD(K)+PSACWI(K))/ &
2929 (PRCI(K)+PRAI(K)+PRACI(K)+PRACIS(K)-EPRD(K))
2931 PRCI(K) = PRCI(K)*RATIO
2932 PRAI(K) = PRAI(K)*RATIO
2933 PRACI(K) = PRACI(K)*RATIO
2934 PRACIS(K) = PRACIS(K)*RATIO
2935 EPRD(K) = EPRD(K)*RATIO
2937 END IF
2939 ! CONSERVATION OF QR
2941 DUM=((PRACS(K)-PRE(K))+(QMULTR(K)+QMULTRG(K)-PRC(K))+(MNUCCR(K)-PRA(K))+ &
2942 PIACR(K)+PIACRS(K)+PGRACS(K)+PRACG(K))*DT
2944 IF (DUM.GT.QR3D(K).AND.QR3D(K).GE.QSMALL) THEN
2946 RATIO = (QR3D(K)/DT+PRC(K)+PRA(K))/ &
2947 (-PRE(K)+QMULTR(K)+QMULTRG(K)+PRACS(K)+MNUCCR(K)+PIACR(K)+PIACRS(K)+PGRACS(K)+PRACG(K))
2949 PRE(K) = PRE(K)*RATIO
2950 PRACS(K) = PRACS(K)*RATIO
2951 QMULTR(K) = QMULTR(K)*RATIO
2952 QMULTRG(K) = QMULTRG(K)*RATIO
2953 MNUCCR(K) = MNUCCR(K)*RATIO
2954 PIACR(K) = PIACR(K)*RATIO
2955 PIACRS(K) = PIACRS(K)*RATIO
2956 PGRACS(K) = PGRACS(K)*RATIO
2957 PRACG(K) = PRACG(K)*RATIO
2959 END IF
2961 ! CONSERVATION OF QNI
2962 ! CONSERVATION FOR GRAUPEL SCHEME
2964 IF (IGRAUP.EQ.0) THEN
2966 DUM = (-PRDS(K)-PSACWS(K)-PRAI(K)-PRCI(K)-PRACS(K)-EPRDS(K)+PSACR(K)-PIACRS(K)-PRACIS(K))*DT
2968 IF (DUM.GT.QNI3D(K).AND.QNI3D(K).GE.QSMALL) THEN
2970 RATIO = (QNI3D(K)/DT+PRDS(K)+PSACWS(K)+PRAI(K)+PRCI(K)+PRACS(K)+PIACRS(K)+PRACIS(K))/(-EPRDS(K)+PSACR(K))
2972 EPRDS(K) = EPRDS(K)*RATIO
2973 PSACR(K) = PSACR(K)*RATIO
2975 END IF
2977 ! FOR NO GRAUPEL, NEED TO INCLUDE FREEZING OF RAIN FOR SNOW
2978 ELSE IF (IGRAUP.EQ.1) THEN
2980 DUM = (-PRDS(K)-PSACWS(K)-PRAI(K)-PRCI(K)-PRACS(K)-EPRDS(K)+PSACR(K)-PIACRS(K)-PRACIS(K)-MNUCCR(K))*DT
2982 IF (DUM.GT.QNI3D(K).AND.QNI3D(K).GE.QSMALL) THEN
2984 RATIO = (QNI3D(K)/DT+PRDS(K)+PSACWS(K)+PRAI(K)+PRCI(K)+PRACS(K)+PIACRS(K)+PRACIS(K)+MNUCCR(K))/(-EPRDS(K)+PSACR(K))
2986 EPRDS(K) = EPRDS(K)*RATIO
2987 PSACR(K) = PSACR(K)*RATIO
2989 END IF
2991 END IF
2993 ! CONSERVATION OF QG
2995 DUM = (-PSACWG(K)-PRACG(K)-PGSACW(K)-PGRACS(K)-PRDG(K)-MNUCCR(K)-EPRDG(K)-PIACR(K)-PRACI(K)-PSACR(K))*DT
2997 IF (DUM.GT.QG3D(K).AND.QG3D(K).GE.QSMALL) THEN
2999 RATIO = (QG3D(K)/DT+PSACWG(K)+PRACG(K)+PGSACW(K)+PGRACS(K)+PRDG(K)+MNUCCR(K)+PSACR(K)+&
3000 PIACR(K)+PRACI(K))/(-EPRDG(K))
3002 EPRDG(K) = EPRDG(K)*RATIO
3004 END IF
3006 ! TENDENCIES
3008 QV3DTEN(K) = QV3DTEN(K)+(-PRE(K)-PRD(K)-PRDS(K)-MNUCCD(K)-EPRD(K)-EPRDS(K)-PRDG(K)-EPRDG(K))
3010 ! BUG FIX HM, 3/1/11, INCLUDE PIACR AND PIACRS
3011 T3DTEN(K) = T3DTEN(K)+(PRE(K) &
3012 *XXLV(K)+(PRD(K)+PRDS(K)+ &
3013 MNUCCD(K)+EPRD(K)+EPRDS(K)+PRDG(K)+EPRDG(K))*XXLS(K)+ &
3014 (PSACWS(K)+PSACWI(K)+MNUCCC(K)+MNUCCR(K)+ &
3015 QMULTS(K)+QMULTG(K)+QMULTR(K)+QMULTRG(K)+PRACS(K) &
3016 +PSACWG(K)+PRACG(K)+PGSACW(K)+PGRACS(K)+PIACR(K)+PIACRS(K))*XLF(K))/CPM(K)
3018 QC3DTEN(K) = QC3DTEN(K)+ &
3019 (-PRA(K)-PRC(K)-MNUCCC(K)+PCC(K)- &
3020 PSACWS(K)-PSACWI(K)-QMULTS(K)-QMULTG(K)-PSACWG(K)-PGSACW(K))
3021 QI3DTEN(K) = QI3DTEN(K)+ &
3022 (PRD(K)+EPRD(K)+PSACWI(K)+MNUCCC(K)-PRCI(K)- &
3023 PRAI(K)+QMULTS(K)+QMULTG(K)+QMULTR(K)+QMULTRG(K)+MNUCCD(K)-PRACI(K)-PRACIS(K))
3024 QR3DTEN(K) = QR3DTEN(K)+ &
3025 (PRE(K)+PRA(K)+PRC(K)-PRACS(K)-MNUCCR(K)-QMULTR(K)-QMULTRG(K) &
3026 -PIACR(K)-PIACRS(K)-PRACG(K)-PGRACS(K))
3028 IF (IGRAUP.EQ.0) THEN
3030 QNI3DTEN(K) = QNI3DTEN(K)+ &
3031 (PRAI(K)+PSACWS(K)+PRDS(K)+PRACS(K)+PRCI(K)+EPRDS(K)-PSACR(K)+PIACRS(K)+PRACIS(K))
3032 NS3DTEN(K) = NS3DTEN(K)+(NSAGG(K)+NPRCI(K)-NSCNG(K)-NGRACS(K)+NIACRS(K))
3033 QG3DTEN(K) = QG3DTEN(K)+(PRACG(K)+PSACWG(K)+PGSACW(K)+PGRACS(K)+ &
3034 PRDG(K)+EPRDG(K)+MNUCCR(K)+PIACR(K)+PRACI(K)+PSACR(K))
3035 NG3DTEN(K) = NG3DTEN(K)+(NSCNG(K)+NGRACS(K)+NNUCCR(K)+NIACR(K))
3037 ! FOR NO GRAUPEL, NEED TO INCLUDE FREEZING OF RAIN FOR SNOW
3038 ELSE IF (IGRAUP.EQ.1) THEN
3040 QNI3DTEN(K) = QNI3DTEN(K)+ &
3041 (PRAI(K)+PSACWS(K)+PRDS(K)+PRACS(K)+PRCI(K)+EPRDS(K)-PSACR(K)+PIACRS(K)+PRACIS(K)+MNUCCR(K))
3042 NS3DTEN(K) = NS3DTEN(K)+(NSAGG(K)+NPRCI(K)-NSCNG(K)-NGRACS(K)+NIACRS(K)+NNUCCR(K))
3044 END IF
3046 NC3DTEN(K) = NC3DTEN(K)+(-NNUCCC(K)-NPSACWS(K) &
3047 -NPRA(K)-NPRC(K)-NPSACWI(K)-NPSACWG(K))
3049 NI3DTEN(K) = NI3DTEN(K)+ &
3050 (NNUCCC(K)-NPRCI(K)-NPRAI(K)+NMULTS(K)+NMULTG(K)+NMULTR(K)+NMULTRG(K)+ &
3051 NNUCCD(K)-NIACR(K)-NIACRS(K))
3053 NR3DTEN(K) = NR3DTEN(K)+(NPRC1(K)-NPRACS(K)-NNUCCR(K) &
3054 +NRAGG(K)-NIACR(K)-NIACRS(K)-NPRACG(K)-NGRACS(K))
3056 ! HM ADD, WRF-CHEM, ADD TENDENCIES FOR C2PREC
3058 C2PREC(K) = PRA(K)+PRC(K)+PSACWS(K)+QMULTS(K)+QMULTG(K)+PSACWG(K)+ &
3059 PGSACW(K)+MNUCCC(K)+PSACWI(K)
3060 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
3061 ! NOW CALCULATE SATURATION ADJUSTMENT TO CONDENSE EXTRA VAPOR ABOVE
3062 ! WATER SATURATION
3064 DUMT = T3D(K)+DT*T3DTEN(K)
3065 DUMQV = QV3D(K)+DT*QV3DTEN(K)
3066 ! hm, add fix for low pressure, 5/12/10
3067 dum=min(0.99*pres(k),POLYSVP(DUMT,0))
3068 DUMQSS = EP_2*dum/(PRES(K)-dum)
3069 DUMQC = QC3D(K)+DT*QC3DTEN(K)
3070 DUMQC = MAX(DUMQC,0.)
3072 ! SATURATION ADJUSTMENT FOR LIQUID
3074 DUMS = DUMQV-DUMQSS
3075 PCC(K) = DUMS/(1.+XXLV(K)**2*DUMQSS/(CPM(K)*RV*DUMT**2))/DT
3076 IF (PCC(K)*DT+DUMQC.LT.0.) THEN
3077 PCC(K) = -DUMQC/DT
3078 END IF
3080 QV3DTEN(K) = QV3DTEN(K)-PCC(K)
3081 T3DTEN(K) = T3DTEN(K)+PCC(K)*XXLV(K)/CPM(K)
3082 QC3DTEN(K) = QC3DTEN(K)+PCC(K)
3084 !.......................................................................
3085 ! ACTIVATION OF CLOUD DROPLETS
3086 ! ACTIVATION OF DROPLET CURRENTLY NOT CALCULATED
3087 ! DROPLET CONCENTRATION IS SPECIFIED !!!!!
3089 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
3090 ! SUBLIMATE, MELT, OR EVAPORATE NUMBER CONCENTRATION
3091 ! THIS FORMULATION ASSUMES 1:1 RATIO BETWEEN MASS LOSS AND
3092 ! LOSS OF NUMBER CONCENTRATION
3094 ! IF (PCC(K).LT.0.) THEN
3095 ! DUM = PCC(K)*DT/QC3D(K)
3096 ! DUM = MAX(-1.,DUM)
3097 ! NSUBC(K) = DUM*NC3D(K)/DT
3098 ! END IF
3100 IF (EPRD(K).LT.0.) THEN
3101 DUM = EPRD(K)*DT/QI3D(K)
3102 DUM = MAX(-1.,DUM)
3103 NSUBI(K) = DUM*NI3D(K)/DT
3104 END IF
3105 IF (EPRDS(K).LT.0.) THEN
3106 DUM = EPRDS(K)*DT/QNI3D(K)
3107 DUM = MAX(-1.,DUM)
3108 NSUBS(K) = DUM*NS3D(K)/DT
3109 END IF
3110 IF (PRE(K).LT.0.) THEN
3111 DUM = PRE(K)*DT/QR3D(K)
3112 DUM = MAX(-1.,DUM)
3113 NSUBR(K) = DUM*NR3D(K)/DT
3114 END IF
3115 IF (EPRDG(K).LT.0.) THEN
3116 DUM = EPRDG(K)*DT/QG3D(K)
3117 DUM = MAX(-1.,DUM)
3118 NSUBG(K) = DUM*NG3D(K)/DT
3119 END IF
3121 ! nsubr(k)=0.
3122 ! nsubs(k)=0.
3123 ! nsubg(k)=0.
3125 ! UPDATE TENDENCIES
3127 ! NC3DTEN(K) = NC3DTEN(K)+NSUBC(K)
3128 NI3DTEN(K) = NI3DTEN(K)+NSUBI(K)
3129 NS3DTEN(K) = NS3DTEN(K)+NSUBS(K)
3130 NG3DTEN(K) = NG3DTEN(K)+NSUBG(K)
3131 NR3DTEN(K) = NR3DTEN(K)+NSUBR(K)
3133 END IF !!!!!! TEMPERATURE
3135 ! SWITCH LTRUE TO 1, SINCE HYDROMETEORS ARE PRESENT
3136 LTRUE = 1
3138 200 CONTINUE
3140 END DO
3142 ! INITIALIZE PRECIP AND SNOW RATES
3143 PRECRT = 0.
3144 SNOWRT = 0.
3146 ! IF THERE ARE NO HYDROMETEORS, THEN SKIP TO END OF SUBROUTINE
3148 IF (LTRUE.EQ.0) GOTO 400
3150 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
3151 !.......................................................................
3152 ! CALCULATE SEDIMENATION
3153 ! THE NUMERICS HERE FOLLOW FROM REISNER ET AL. (1998)
3154 ! FALLOUT TERMS ARE CALCULATED ON SPLIT TIME STEPS TO ENSURE NUMERICAL
3155 ! STABILITY, I.E. COURANT# < 1
3157 !.......................................................................
3159 NSTEP = 1
3161 DO K = KTE,KTS,-1
3163 DUMI(K) = QI3D(K)+QI3DTEN(K)*DT
3164 DUMQS(K) = QNI3D(K)+QNI3DTEN(K)*DT
3165 DUMR(K) = QR3D(K)+QR3DTEN(K)*DT
3166 DUMFNI(K) = NI3D(K)+NI3DTEN(K)*DT
3167 DUMFNS(K) = NS3D(K)+NS3DTEN(K)*DT
3168 DUMFNR(K) = NR3D(K)+NR3DTEN(K)*DT
3169 DUMC(K) = QC3D(K)+QC3DTEN(K)*DT
3170 DUMFNC(K) = NC3D(K)+NC3DTEN(K)*DT
3171 DUMG(K) = QG3D(K)+QG3DTEN(K)*DT
3172 DUMFNG(K) = NG3D(K)+NG3DTEN(K)*DT
3174 ! SWITCH FOR CONSTANT DROPLET NUMBER
3175 IF (iinum.EQ.1) THEN
3176 DUMFNC(K) = NC3D(K)
3177 END IF
3179 ! GET DUMMY LAMDA FOR SEDIMENTATION CALCULATIONS
3181 ! MAKE SURE NUMBER CONCENTRATIONS ARE POSITIVE
3182 DUMFNI(K) = MAX(0.,DUMFNI(K))
3183 DUMFNS(K) = MAX(0.,DUMFNS(K))
3184 DUMFNC(K) = MAX(0.,DUMFNC(K))
3185 DUMFNR(K) = MAX(0.,DUMFNR(K))
3186 DUMFNG(K) = MAX(0.,DUMFNG(K))
3188 !......................................................................
3189 ! CLOUD ICE
3191 IF (DUMI(K).GE.QSMALL) THEN
3192 DLAMI = (CONS12*DUMFNI(K)/DUMI(K))**(1./DI)
3193 DLAMI=MAX(DLAMI,LAMMINI)
3194 DLAMI=MIN(DLAMI,LAMMAXI)
3195 END IF
3196 !......................................................................
3197 ! RAIN
3199 IF (DUMR(K).GE.QSMALL) THEN
3200 DLAMR = (PI*RHOW*DUMFNR(K)/DUMR(K))**(1./3.)
3201 DLAMR=MAX(DLAMR,LAMMINR)
3202 DLAMR=MIN(DLAMR,LAMMAXR)
3203 END IF
3204 !......................................................................
3205 ! CLOUD DROPLETS
3207 IF (DUMC(K).GE.QSMALL) THEN
3208 DUM = PRES(K)/(287.15*T3D(K))
3209 PGAM(K)=0.0005714*(NC3D(K)/1.E6*DUM)+0.2714
3210 PGAM(K)=1./(PGAM(K)**2)-1.
3211 PGAM(K)=MAX(PGAM(K),2.)
3212 PGAM(K)=MIN(PGAM(K),10.)
3214 DLAMC = (CONS26*DUMFNC(K)*GAMMA(PGAM(K)+4.)/(DUMC(K)*GAMMA(PGAM(K)+1.)))**(1./3.)
3215 LAMMIN = (PGAM(K)+1.)/60.E-6
3216 LAMMAX = (PGAM(K)+1.)/1.E-6
3217 DLAMC=MAX(DLAMC,LAMMIN)
3218 DLAMC=MIN(DLAMC,LAMMAX)
3219 END IF
3220 !......................................................................
3221 ! SNOW
3223 IF (DUMQS(K).GE.QSMALL) THEN
3224 DLAMS = (CONS1*DUMFNS(K)/ DUMQS(K))**(1./DS)
3225 DLAMS=MAX(DLAMS,LAMMINS)
3226 DLAMS=MIN(DLAMS,LAMMAXS)
3227 END IF
3228 !......................................................................
3229 ! GRAUPEL
3231 IF (DUMG(K).GE.QSMALL) THEN
3232 DLAMG = (CONS2*DUMFNG(K)/ DUMG(K))**(1./DG)
3233 DLAMG=MAX(DLAMG,LAMMING)
3234 DLAMG=MIN(DLAMG,LAMMAXG)
3235 END IF
3237 !......................................................................
3238 ! CALCULATE NUMBER-WEIGHTED AND MASS-WEIGHTED TERMINAL FALL SPEEDS
3240 ! CLOUD WATER
3242 IF (DUMC(K).GE.QSMALL) THEN
3243 UNC = ACN(K)*GAMMA(1.+BC+PGAM(K))/ (DLAMC**BC*GAMMA(PGAM(K)+1.))
3244 UMC = ACN(K)*GAMMA(4.+BC+PGAM(K))/ (DLAMC**BC*GAMMA(PGAM(K)+4.))
3245 ELSE
3246 UMC = 0.
3247 UNC = 0.
3248 END IF
3250 IF (DUMI(K).GE.QSMALL) THEN
3251 UNI = AIN(K)*CONS27/DLAMI**BI
3252 UMI = AIN(K)*CONS28/(DLAMI**BI)
3253 ELSE
3254 UMI = 0.
3255 UNI = 0.
3256 END IF
3258 IF (DUMR(K).GE.QSMALL) THEN
3259 UNR = ARN(K)*CONS6/DLAMR**BR
3260 UMR = ARN(K)*CONS4/(DLAMR**BR)
3261 ELSE
3262 UMR = 0.
3263 UNR = 0.
3264 END IF
3266 IF (DUMQS(K).GE.QSMALL) THEN
3267 UMS = ASN(K)*CONS3/(DLAMS**BS)
3268 UNS = ASN(K)*CONS5/DLAMS**BS
3269 ELSE
3270 UMS = 0.
3271 UNS = 0.
3272 END IF
3274 IF (DUMG(K).GE.QSMALL) THEN
3275 UMG = AGN(K)*CONS7/(DLAMG**BG)
3276 UNG = AGN(K)*CONS8/DLAMG**BG
3277 ELSE
3278 UMG = 0.
3279 UNG = 0.
3280 END IF
3282 ! SET REALISTIC LIMITS ON FALLSPEED
3284 ! bug fix, 10/08/09
3285 dum=(rhosu/rho(k))**0.54
3286 UMS=MIN(UMS,1.2*dum)
3287 UNS=MIN(UNS,1.2*dum)
3288 ! fix 053011
3289 ! fix for correction by AA 4/6/11
3290 UMI=MIN(UMI,1.2*(rhosu/rho(k))**0.35)
3291 UNI=MIN(UNI,1.2*(rhosu/rho(k))**0.35)
3292 UMR=MIN(UMR,9.1*dum)
3293 UNR=MIN(UNR,9.1*dum)
3294 UMG=MIN(UMG,20.*dum)
3295 UNG=MIN(UNG,20.*dum)
3297 FR(K) = UMR
3298 FI(K) = UMI
3299 FNI(K) = UNI
3300 FS(K) = UMS
3301 FNS(K) = UNS
3302 FNR(K) = UNR
3303 FC(K) = UMC
3304 FNC(K) = UNC
3305 FG(K) = UMG
3306 FNG(K) = UNG
3308 ! V3.3 MODIFY FALLSPEED BELOW LEVEL OF PRECIP
3310 IF (K.LE.KTE-1) THEN
3311 IF (FR(K).LT.1.E-10) THEN
3312 FR(K)=FR(K+1)
3313 END IF
3314 IF (FI(K).LT.1.E-10) THEN
3315 FI(K)=FI(K+1)
3316 END IF
3317 IF (FNI(K).LT.1.E-10) THEN
3318 FNI(K)=FNI(K+1)
3319 END IF
3320 IF (FS(K).LT.1.E-10) THEN
3321 FS(K)=FS(K+1)
3322 END IF
3323 IF (FNS(K).LT.1.E-10) THEN
3324 FNS(K)=FNS(K+1)
3325 END IF
3326 IF (FNR(K).LT.1.E-10) THEN
3327 FNR(K)=FNR(K+1)
3328 END IF
3329 IF (FC(K).LT.1.E-10) THEN
3330 FC(K)=FC(K+1)
3331 END IF
3332 IF (FNC(K).LT.1.E-10) THEN
3333 FNC(K)=FNC(K+1)
3334 END IF
3335 IF (FG(K).LT.1.E-10) THEN
3336 FG(K)=FG(K+1)
3337 END IF
3338 IF (FNG(K).LT.1.E-10) THEN
3339 FNG(K)=FNG(K+1)
3340 END IF
3341 END IF ! K LE KTE-1
3343 ! CALCULATE NUMBER OF SPLIT TIME STEPS
3345 RGVM = MAX(FR(K),FI(K),FS(K),FC(K),FNI(K),FNR(K),FNS(K),FNC(K),FG(K),FNG(K))
3346 ! VVT CHANGED IFIX -> INT (GENERIC FUNCTION)
3347 NSTEP = MAX(INT(RGVM*DT/DZQ(K)+1.),NSTEP)
3349 ! MULTIPLY VARIABLES BY RHO
3350 DUMR(k) = DUMR(k)*RHO(K)
3351 DUMI(k) = DUMI(k)*RHO(K)
3352 DUMFNI(k) = DUMFNI(K)*RHO(K)
3353 DUMQS(k) = DUMQS(K)*RHO(K)
3354 DUMFNS(k) = DUMFNS(K)*RHO(K)
3355 DUMFNR(k) = DUMFNR(K)*RHO(K)
3356 DUMC(k) = DUMC(K)*RHO(K)
3357 DUMFNC(k) = DUMFNC(K)*RHO(K)
3358 DUMG(k) = DUMG(K)*RHO(K)
3359 DUMFNG(k) = DUMFNG(K)*RHO(K)
3361 END DO
3363 DO N = 1,NSTEP
3365 DO K = KTS,KTE
3366 FALOUTR(K) = FR(K)*DUMR(K)
3367 FALOUTI(K) = FI(K)*DUMI(K)
3368 FALOUTNI(K) = FNI(K)*DUMFNI(K)
3369 FALOUTS(K) = FS(K)*DUMQS(K)
3370 FALOUTNS(K) = FNS(K)*DUMFNS(K)
3371 FALOUTNR(K) = FNR(K)*DUMFNR(K)
3372 FALOUTC(K) = FC(K)*DUMC(K)
3373 FALOUTNC(K) = FNC(K)*DUMFNC(K)
3374 FALOUTG(K) = FG(K)*DUMG(K)
3375 FALOUTNG(K) = FNG(K)*DUMFNG(K)
3376 END DO
3378 ! TOP OF MODEL
3380 K = KTE
3381 FALTNDR = FALOUTR(K)/DZQ(k)
3382 FALTNDI = FALOUTI(K)/DZQ(k)
3383 FALTNDNI = FALOUTNI(K)/DZQ(k)
3384 FALTNDS = FALOUTS(K)/DZQ(k)
3385 FALTNDNS = FALOUTNS(K)/DZQ(k)
3386 FALTNDNR = FALOUTNR(K)/DZQ(k)
3387 FALTNDC = FALOUTC(K)/DZQ(k)
3388 FALTNDNC = FALOUTNC(K)/DZQ(k)
3389 FALTNDG = FALOUTG(K)/DZQ(k)
3390 FALTNDNG = FALOUTNG(K)/DZQ(k)
3391 ! ADD FALLOUT TERMS TO EULERIAN TENDENCIES
3393 QRSTEN(K) = QRSTEN(K)-FALTNDR/NSTEP/RHO(k)
3394 QISTEN(K) = QISTEN(K)-FALTNDI/NSTEP/RHO(k)
3395 NI3DTEN(K) = NI3DTEN(K)-FALTNDNI/NSTEP/RHO(k)
3396 QNISTEN(K) = QNISTEN(K)-FALTNDS/NSTEP/RHO(k)
3397 NS3DTEN(K) = NS3DTEN(K)-FALTNDNS/NSTEP/RHO(k)
3398 NR3DTEN(K) = NR3DTEN(K)-FALTNDNR/NSTEP/RHO(k)
3399 QCSTEN(K) = QCSTEN(K)-FALTNDC/NSTEP/RHO(k)
3400 NC3DTEN(K) = NC3DTEN(K)-FALTNDNC/NSTEP/RHO(k)
3401 QGSTEN(K) = QGSTEN(K)-FALTNDG/NSTEP/RHO(k)
3402 NG3DTEN(K) = NG3DTEN(K)-FALTNDNG/NSTEP/RHO(k)
3404 DUMR(K) = DUMR(K)-FALTNDR*DT/NSTEP
3405 DUMI(K) = DUMI(K)-FALTNDI*DT/NSTEP
3406 DUMFNI(K) = DUMFNI(K)-FALTNDNI*DT/NSTEP
3407 DUMQS(K) = DUMQS(K)-FALTNDS*DT/NSTEP
3408 DUMFNS(K) = DUMFNS(K)-FALTNDNS*DT/NSTEP
3409 DUMFNR(K) = DUMFNR(K)-FALTNDNR*DT/NSTEP
3410 DUMC(K) = DUMC(K)-FALTNDC*DT/NSTEP
3411 DUMFNC(K) = DUMFNC(K)-FALTNDNC*DT/NSTEP
3412 DUMG(K) = DUMG(K)-FALTNDG*DT/NSTEP
3413 DUMFNG(K) = DUMFNG(K)-FALTNDNG*DT/NSTEP
3415 DO K = KTE-1,KTS,-1
3416 FALTNDR = (FALOUTR(K+1)-FALOUTR(K))/DZQ(K)
3417 FALTNDI = (FALOUTI(K+1)-FALOUTI(K))/DZQ(K)
3418 FALTNDNI = (FALOUTNI(K+1)-FALOUTNI(K))/DZQ(K)
3419 FALTNDS = (FALOUTS(K+1)-FALOUTS(K))/DZQ(K)
3420 FALTNDNS = (FALOUTNS(K+1)-FALOUTNS(K))/DZQ(K)
3421 FALTNDNR = (FALOUTNR(K+1)-FALOUTNR(K))/DZQ(K)
3422 FALTNDC = (FALOUTC(K+1)-FALOUTC(K))/DZQ(K)
3423 FALTNDNC = (FALOUTNC(K+1)-FALOUTNC(K))/DZQ(K)
3424 FALTNDG = (FALOUTG(K+1)-FALOUTG(K))/DZQ(K)
3425 FALTNDNG = (FALOUTNG(K+1)-FALOUTNG(K))/DZQ(K)
3427 ! ADD FALLOUT TERMS TO EULERIAN TENDENCIES
3429 QRSTEN(K) = QRSTEN(K)+FALTNDR/NSTEP/RHO(k)
3430 QISTEN(K) = QISTEN(K)+FALTNDI/NSTEP/RHO(k)
3431 NI3DTEN(K) = NI3DTEN(K)+FALTNDNI/NSTEP/RHO(k)
3432 QNISTEN(K) = QNISTEN(K)+FALTNDS/NSTEP/RHO(k)
3433 NS3DTEN(K) = NS3DTEN(K)+FALTNDNS/NSTEP/RHO(k)
3434 NR3DTEN(K) = NR3DTEN(K)+FALTNDNR/NSTEP/RHO(k)
3435 QCSTEN(K) = QCSTEN(K)+FALTNDC/NSTEP/RHO(k)
3436 NC3DTEN(K) = NC3DTEN(K)+FALTNDNC/NSTEP/RHO(k)
3437 QGSTEN(K) = QGSTEN(K)+FALTNDG/NSTEP/RHO(k)
3438 NG3DTEN(K) = NG3DTEN(K)+FALTNDNG/NSTEP/RHO(k)
3440 DUMR(K) = DUMR(K)+FALTNDR*DT/NSTEP
3441 DUMI(K) = DUMI(K)+FALTNDI*DT/NSTEP
3442 DUMFNI(K) = DUMFNI(K)+FALTNDNI*DT/NSTEP
3443 DUMQS(K) = DUMQS(K)+FALTNDS*DT/NSTEP
3444 DUMFNS(K) = DUMFNS(K)+FALTNDNS*DT/NSTEP
3445 DUMFNR(K) = DUMFNR(K)+FALTNDNR*DT/NSTEP
3446 DUMC(K) = DUMC(K)+FALTNDC*DT/NSTEP
3447 DUMFNC(K) = DUMFNC(K)+FALTNDNC*DT/NSTEP
3448 DUMG(K) = DUMG(K)+FALTNDG*DT/NSTEP
3449 DUMFNG(K) = DUMFNG(K)+FALTNDNG*DT/NSTEP
3451 ! FOR WRF-CHEM, NEED PRECIP RATES (UNITS OF KG/M^2/S)
3452 CSED(K)=CSED(K)+FALOUTC(K)/NSTEP
3453 ISED(K)=ISED(K)+FALOUTI(K)/NSTEP
3454 SSED(K)=SSED(K)+FALOUTS(K)/NSTEP
3455 GSED(K)=GSED(K)+FALOUTG(K)/NSTEP
3456 RSED(K)=RSED(K)+FALOUTR(K)/NSTEP
3457 END DO
3459 ! GET PRECIPITATION AND SNOWFALL ACCUMULATION DURING THE TIME STEP
3460 ! FACTOR OF 1000 CONVERTS FROM M TO MM, BUT DIVISION BY DENSITY
3461 ! OF LIQUID WATER CANCELS THIS FACTOR OF 1000
3463 PRECRT = PRECRT+(FALOUTR(KTS)+FALOUTC(KTS)+FALOUTS(KTS)+FALOUTI(KTS)+FALOUTG(KTS)) &
3464 *DT/NSTEP
3465 SNOWRT = SNOWRT+(FALOUTS(KTS)+FALOUTI(KTS)+FALOUTG(KTS))*DT/NSTEP
3467 END DO
3469 DO K=KTS,KTE
3471 ! ADD ON SEDIMENTATION TENDENCIES FOR MIXING RATIO TO REST OF TENDENCIES
3473 QR3DTEN(K)=QR3DTEN(K)+QRSTEN(K)
3474 QI3DTEN(K)=QI3DTEN(K)+QISTEN(K)
3475 QC3DTEN(K)=QC3DTEN(K)+QCSTEN(K)
3476 QG3DTEN(K)=QG3DTEN(K)+QGSTEN(K)
3477 QNI3DTEN(K)=QNI3DTEN(K)+QNISTEN(K)
3479 ! PUT ALL CLOUD ICE IN SNOW CATEGORY IF MEAN DIAMETER EXCEEDS 2 * dcs
3481 !hm 4/7/09 bug fix
3482 ! IF (QI3D(K).GE.QSMALL.AND.T3D(K).LT.273.15) THEN
3483 IF (QI3D(K).GE.QSMALL.AND.T3D(K).LT.273.15.AND.LAMI(K).GE.1.E-10) THEN
3484 IF (1./LAMI(K).GE.2.*DCS) THEN
3485 QNI3DTEN(K) = QNI3DTEN(K)+QI3D(K)/DT+ QI3DTEN(K)
3486 NS3DTEN(K) = NS3DTEN(K)+NI3D(K)/DT+ NI3DTEN(K)
3487 QI3DTEN(K) = -QI3D(K)/DT
3488 NI3DTEN(K) = -NI3D(K)/DT
3489 END IF
3490 END IF
3492 ! hm add tendencies here, then call sizeparameter
3493 ! to ensure consisitency between mixing ratio and number concentration
3495 QC3D(k) = QC3D(k)+QC3DTEN(k)*DT
3496 QI3D(k) = QI3D(k)+QI3DTEN(k)*DT
3497 QNI3D(k) = QNI3D(k)+QNI3DTEN(k)*DT
3498 QR3D(k) = QR3D(k)+QR3DTEN(k)*DT
3499 NC3D(k) = NC3D(k)+NC3DTEN(k)*DT
3500 NI3D(k) = NI3D(k)+NI3DTEN(k)*DT
3501 NS3D(k) = NS3D(k)+NS3DTEN(k)*DT
3502 NR3D(k) = NR3D(k)+NR3DTEN(k)*DT
3504 IF (IGRAUP.EQ.0) THEN
3505 QG3D(k) = QG3D(k)+QG3DTEN(k)*DT
3506 NG3D(k) = NG3D(k)+NG3DTEN(k)*DT
3507 END IF
3509 ! ADD TEMPERATURE AND WATER VAPOR TENDENCIES FROM MICROPHYSICS
3510 T3D(K) = T3D(K)+T3DTEN(k)*DT
3511 QV3D(K) = QV3D(K)+QV3DTEN(k)*DT
3513 ! SATURATION VAPOR PRESSURE AND MIXING RATIO
3515 ! hm, add fix for low pressure, 5/12/10
3516 EVS(K) = min(0.99*pres(k),POLYSVP(T3D(K),0)) ! PA
3517 EIS(K) = min(0.99*pres(k),POLYSVP(T3D(K),1)) ! PA
3519 ! MAKE SURE ICE SATURATION DOESN'T EXCEED WATER SAT. NEAR FREEZING
3521 IF (EIS(K).GT.EVS(K)) EIS(K) = EVS(K)
3523 QVS(K) = EP_2*EVS(K)/(PRES(K)-EVS(K))
3524 QVI(K) = EP_2*EIS(K)/(PRES(K)-EIS(K))
3526 QVQVS(K) = QV3D(K)/QVS(K)
3527 QVQVSI(K) = QV3D(K)/QVI(K)
3529 ! AT SUBSATURATION, REMOVE SMALL AMOUNTS OF CLOUD/PRECIP WATER
3530 ! hm 7/9/09 change limit to 1.e-8
3532 IF (QVQVS(K).LT.0.9) THEN
3533 IF (QR3D(K).LT.1.E-8) THEN
3534 QV3D(K)=QV3D(K)+QR3D(K)
3535 T3D(K)=T3D(K)-QR3D(K)*XXLV(K)/CPM(K)
3536 QR3D(K)=0.
3537 END IF
3538 IF (QC3D(K).LT.1.E-8) THEN
3539 QV3D(K)=QV3D(K)+QC3D(K)
3540 T3D(K)=T3D(K)-QC3D(K)*XXLV(K)/CPM(K)
3541 QC3D(K)=0.
3542 END IF
3543 END IF
3545 IF (QVQVSI(K).LT.0.9) THEN
3546 IF (QI3D(K).LT.1.E-8) THEN
3547 QV3D(K)=QV3D(K)+QI3D(K)
3548 T3D(K)=T3D(K)-QI3D(K)*XXLS(K)/CPM(K)
3549 QI3D(K)=0.
3550 END IF
3551 IF (QNI3D(K).LT.1.E-8) THEN
3552 QV3D(K)=QV3D(K)+QNI3D(K)
3553 T3D(K)=T3D(K)-QNI3D(K)*XXLS(K)/CPM(K)
3554 QNI3D(K)=0.
3555 END IF
3556 IF (QG3D(K).LT.1.E-8) THEN
3557 QV3D(K)=QV3D(K)+QG3D(K)
3558 T3D(K)=T3D(K)-QG3D(K)*XXLS(K)/CPM(K)
3559 QG3D(K)=0.
3560 END IF
3561 END IF
3563 !..................................................................
3564 ! IF MIXING RATIO < QSMALL SET MIXING RATIO AND NUMBER CONC TO ZERO
3566 IF (QC3D(K).LT.QSMALL) THEN
3567 QC3D(K) = 0.
3568 NC3D(K) = 0.
3569 EFFC(K) = 0.
3570 END IF
3571 IF (QR3D(K).LT.QSMALL) THEN
3572 QR3D(K) = 0.
3573 NR3D(K) = 0.
3574 EFFR(K) = 0.
3575 END IF
3576 IF (QI3D(K).LT.QSMALL) THEN
3577 QI3D(K) = 0.
3578 NI3D(K) = 0.
3579 EFFI(K) = 0.
3580 END IF
3581 IF (QNI3D(K).LT.QSMALL) THEN
3582 QNI3D(K) = 0.
3583 NS3D(K) = 0.
3584 EFFS(K) = 0.
3585 END IF
3586 IF (QG3D(K).LT.QSMALL) THEN
3587 QG3D(K) = 0.
3588 NG3D(K) = 0.
3589 EFFG(K) = 0.
3590 END IF
3592 !..................................
3593 ! IF THERE IS NO CLOUD/PRECIP WATER, THEN SKIP CALCULATIONS
3595 IF (QC3D(K).LT.QSMALL.AND.QI3D(K).LT.QSMALL.AND.QNI3D(K).LT.QSMALL &
3596 .AND.QR3D(K).LT.QSMALL.AND.QG3D(K).LT.QSMALL) GOTO 500
3598 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
3599 ! CALCULATE INSTANTANEOUS PROCESSES
3601 ! ADD MELTING OF CLOUD ICE TO FORM RAIN
3603 IF (QI3D(K).GE.QSMALL.AND.T3D(K).GE.273.15) THEN
3604 QR3D(K) = QR3D(K)+QI3D(K)
3605 T3D(K) = T3D(K)-QI3D(K)*XLF(K)/CPM(K)
3606 QI3D(K) = 0.
3607 NR3D(K) = NR3D(K)+NI3D(K)
3608 NI3D(K) = 0.
3609 END IF
3611 ! ****SENSITIVITY - NO ICE
3612 IF (ILIQ.EQ.1) GOTO 778
3614 ! HOMOGENEOUS FREEZING OF CLOUD WATER
3616 IF (T3D(K).LE.233.15.AND.QC3D(K).GE.QSMALL) THEN
3617 QI3D(K)=QI3D(K)+QC3D(K)
3618 T3D(K)=T3D(K)+QC3D(K)*XLF(K)/CPM(K)
3619 QC3D(K)=0.
3620 NI3D(K)=NI3D(K)+NC3D(K)
3621 NC3D(K)=0.
3622 END IF
3624 ! HOMOGENEOUS FREEZING OF RAIN
3626 IF (IGRAUP.EQ.0) THEN
3628 IF (T3D(K).LE.233.15.AND.QR3D(K).GE.QSMALL) THEN
3629 QG3D(K) = QG3D(K)+QR3D(K)
3630 T3D(K) = T3D(K)+QR3D(K)*XLF(K)/CPM(K)
3631 QR3D(K) = 0.
3632 NG3D(K) = NG3D(K)+ NR3D(K)
3633 NR3D(K) = 0.
3634 END IF
3636 ELSE IF (IGRAUP.EQ.1) THEN
3638 IF (T3D(K).LE.233.15.AND.QR3D(K).GE.QSMALL) THEN
3639 QNI3D(K) = QNI3D(K)+QR3D(K)
3640 T3D(K) = T3D(K)+QR3D(K)*XLF(K)/CPM(K)
3641 QR3D(K) = 0.
3642 NS3D(K) = NS3D(K)+NR3D(K)
3643 NR3D(K) = 0.
3644 END IF
3646 END IF
3648 778 CONTINUE
3650 ! MAKE SURE NUMBER CONCENTRATIONS AREN'T NEGATIVE
3652 NI3D(K) = MAX(0.,NI3D(K))
3653 NS3D(K) = MAX(0.,NS3D(K))
3654 NC3D(K) = MAX(0.,NC3D(K))
3655 NR3D(K) = MAX(0.,NR3D(K))
3656 NG3D(K) = MAX(0.,NG3D(K))
3658 !......................................................................
3659 ! CLOUD ICE
3661 IF (QI3D(K).GE.QSMALL) THEN
3662 LAMI(K) = (CONS12* &
3663 NI3D(K)/QI3D(K))**(1./DI)
3665 ! CHECK FOR SLOPE
3667 ! ADJUST VARS
3669 IF (LAMI(K).LT.LAMMINI) THEN
3671 LAMI(K) = LAMMINI
3673 N0I(K) = LAMI(K)**4*QI3D(K)/CONS12
3675 NI3D(K) = N0I(K)/LAMI(K)
3676 ELSE IF (LAMI(K).GT.LAMMAXI) THEN
3677 LAMI(K) = LAMMAXI
3678 N0I(K) = LAMI(K)**4*QI3D(K)/CONS12
3680 NI3D(K) = N0I(K)/LAMI(K)
3681 END IF
3682 END IF
3684 !......................................................................
3685 ! RAIN
3687 IF (QR3D(K).GE.QSMALL) THEN
3688 LAMR(K) = (PI*RHOW*NR3D(K)/QR3D(K))**(1./3.)
3690 ! CHECK FOR SLOPE
3692 ! ADJUST VARS
3694 IF (LAMR(K).LT.LAMMINR) THEN
3696 LAMR(K) = LAMMINR
3698 N0RR(K) = LAMR(K)**4*QR3D(K)/(PI*RHOW)
3700 NR3D(K) = N0RR(K)/LAMR(K)
3701 ELSE IF (LAMR(K).GT.LAMMAXR) THEN
3702 LAMR(K) = LAMMAXR
3703 N0RR(K) = LAMR(K)**4*QR3D(K)/(PI*RHOW)
3705 NR3D(K) = N0RR(K)/LAMR(K)
3706 END IF
3708 END IF
3710 !......................................................................
3711 ! CLOUD DROPLETS
3713 ! MARTIN ET AL. (1994) FORMULA FOR PGAM
3715 IF (QC3D(K).GE.QSMALL) THEN
3717 DUM = PRES(K)/(287.15*T3D(K))
3718 PGAM(K)=0.0005714*(NC3D(K)/1.E6*DUM)+0.2714
3719 PGAM(K)=1./(PGAM(K)**2)-1.
3720 PGAM(K)=MAX(PGAM(K),2.)
3721 PGAM(K)=MIN(PGAM(K),10.)
3723 ! CALCULATE LAMC
3725 LAMC(K) = (CONS26*NC3D(K)*GAMMA(PGAM(K)+4.)/ &
3726 (QC3D(K)*GAMMA(PGAM(K)+1.)))**(1./3.)
3728 ! LAMMIN, 60 MICRON DIAMETER
3729 ! LAMMAX, 1 MICRON
3731 LAMMIN = (PGAM(K)+1.)/60.E-6
3732 LAMMAX = (PGAM(K)+1.)/1.E-6
3734 IF (LAMC(K).LT.LAMMIN) THEN
3735 LAMC(K) = LAMMIN
3736 NC3D(K) = EXP(3.*LOG(LAMC(K))+LOG(QC3D(K))+ &
3737 LOG(GAMMA(PGAM(K)+1.))-LOG(GAMMA(PGAM(K)+4.)))/CONS26
3739 ELSE IF (LAMC(K).GT.LAMMAX) THEN
3740 LAMC(K) = LAMMAX
3741 NC3D(K) = EXP(3.*LOG(LAMC(K))+LOG(QC3D(K))+ &
3742 LOG(GAMMA(PGAM(K)+1.))-LOG(GAMMA(PGAM(K)+4.)))/CONS26
3744 END IF
3746 END IF
3748 !......................................................................
3749 ! SNOW
3751 IF (QNI3D(K).GE.QSMALL) THEN
3752 LAMS(K) = (CONS1*NS3D(K)/QNI3D(K))**(1./DS)
3754 ! CHECK FOR SLOPE
3756 ! ADJUST VARS
3758 IF (LAMS(K).LT.LAMMINS) THEN
3759 LAMS(K) = LAMMINS
3760 N0S(K) = LAMS(K)**4*QNI3D(K)/CONS1
3762 NS3D(K) = N0S(K)/LAMS(K)
3764 ELSE IF (LAMS(K).GT.LAMMAXS) THEN
3766 LAMS(K) = LAMMAXS
3767 N0S(K) = LAMS(K)**4*QNI3D(K)/CONS1
3768 NS3D(K) = N0S(K)/LAMS(K)
3769 END IF
3771 END IF
3773 !......................................................................
3774 ! GRAUPEL
3776 IF (QG3D(K).GE.QSMALL) THEN
3777 LAMG(K) = (CONS2*NG3D(K)/QG3D(K))**(1./DG)
3779 ! CHECK FOR SLOPE
3781 ! ADJUST VARS
3783 IF (LAMG(K).LT.LAMMING) THEN
3784 LAMG(K) = LAMMING
3785 N0G(K) = LAMG(K)**4*QG3D(K)/CONS2
3787 NG3D(K) = N0G(K)/LAMG(K)
3789 ELSE IF (LAMG(K).GT.LAMMAXG) THEN
3791 LAMG(K) = LAMMAXG
3792 N0G(K) = LAMG(K)**4*QG3D(K)/CONS2
3794 NG3D(K) = N0G(K)/LAMG(K)
3795 END IF
3797 END IF
3799 500 CONTINUE
3801 ! CALCULATE EFFECTIVE RADIUS
3803 IF (QI3D(K).GE.QSMALL) THEN
3804 EFFI(K) = 3./LAMI(K)/2.*1.E6
3805 ELSE
3806 EFFI(K) = 25.
3807 END IF
3809 IF (QNI3D(K).GE.QSMALL) THEN
3810 EFFS(K) = 3./LAMS(K)/2.*1.E6
3811 ELSE
3812 EFFS(K) = 25.
3813 END IF
3815 IF (QR3D(K).GE.QSMALL) THEN
3816 EFFR(K) = 3./LAMR(K)/2.*1.E6
3817 ELSE
3818 EFFR(K) = 25.
3819 END IF
3821 IF (QC3D(K).GE.QSMALL) THEN
3822 EFFC(K) = GAMMA(PGAM(K)+4.)/ &
3823 GAMMA(PGAM(K)+3.)/LAMC(K)/2.*1.E6
3824 ELSE
3825 EFFC(K) = 25.
3826 END IF
3828 IF (QG3D(K).GE.QSMALL) THEN
3829 EFFG(K) = 3./LAMG(K)/2.*1.E6
3830 ELSE
3831 EFFG(K) = 25.
3832 END IF
3834 ! HM ADD 1/10/06, ADD UPPER BOUND ON ICE NUMBER, THIS IS NEEDED
3835 ! TO PREVENT VERY LARGE ICE NUMBER DUE TO HOMOGENEOUS FREEZING
3836 ! OF DROPLETS, ESPECIALLY WHEN INUM = 1, SET MAX AT 10 CM-3
3837 NI3D(K) = MIN(NI3D(K),10.E6/RHO(K))
3838 ! ADD BOUND ON DROPLET NUMBER - CANNOT EXCEED AEROSOL CONCENTRATION
3839 IF (iinum.EQ.0.AND.IACT.EQ.2) THEN
3840 NC3D(K) = MIN(NC3D(K),(NANEW1+NANEW2)/RHO(K))
3841 END IF
3842 ! SWITCH FOR CONSTANT DROPLET NUMBER
3843 IF (iinum.EQ.1) THEN
3844 ! CHANGE NDCNST FROM CM-3 TO KG-1
3845 NC3D(K) = NDCNST*1.E6/RHO(K)
3846 END IF
3848 END DO !!! K LOOP
3850 400 CONTINUE
3852 ! ALL DONE !!!!!!!!!!!
3853 RETURN
3854 END SUBROUTINE MORR_TWO_MOMENT_MICRO
3856 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
3858 REAL FUNCTION POLYSVP (T,TYPE)
3860 !-------------------------------------------
3862 ! COMPUTE SATURATION VAPOR PRESSURE
3864 ! POLYSVP RETURNED IN UNITS OF PA.
3865 ! T IS INPUT IN UNITS OF K.
3866 ! TYPE REFERS TO SATURATION WITH RESPECT TO LIQUID (0) OR ICE (1)
3868 ! REPLACE GOFF-GRATCH WITH FASTER FORMULATION FROM FLATAU ET AL. 1992, TABLE 4 (RIGHT-HAND COLUMN)
3870 IMPLICIT NONE
3872 REAL DUM
3873 REAL T
3874 INTEGER TYPE
3875 ! ice
3876 real a0i,a1i,a2i,a3i,a4i,a5i,a6i,a7i,a8i
3877 data a0i,a1i,a2i,a3i,a4i,a5i,a6i,a7i,a8i /&
3878 6.11147274, 0.503160820, 0.188439774e-1, &
3879 0.420895665e-3, 0.615021634e-5,0.602588177e-7, &
3880 0.385852041e-9, 0.146898966e-11, 0.252751365e-14/
3882 ! liquid
3883 real a0,a1,a2,a3,a4,a5,a6,a7,a8
3885 ! V1.7
3886 data a0,a1,a2,a3,a4,a5,a6,a7,a8 /&
3887 6.11239921, 0.443987641, 0.142986287e-1, &
3888 0.264847430e-3, 0.302950461e-5, 0.206739458e-7, &
3889 0.640689451e-10,-0.952447341e-13,-0.976195544e-15/
3890 real dt
3892 ! ICE
3894 IF (TYPE.EQ.1) THEN
3896 ! POLYSVP = 10.**(-9.09718*(273.16/T-1.)-3.56654* &
3897 ! LOG10(273.16/T)+0.876793*(1.-T/273.16)+ &
3898 ! LOG10(6.1071))*100.
3901 dt = max(-80.,t-273.16)
3902 polysvp = a0i + dt*(a1i+dt*(a2i+dt*(a3i+dt*(a4i+dt*(a5i+dt*(a6i+dt*(a7i+a8i*dt)))))))
3903 polysvp = polysvp*100.
3905 END IF
3907 ! LIQUID
3909 IF (TYPE.EQ.0) THEN
3911 dt = max(-80.,t-273.16)
3912 polysvp = a0 + dt*(a1+dt*(a2+dt*(a3+dt*(a4+dt*(a5+dt*(a6+dt*(a7+a8*dt)))))))
3913 polysvp = polysvp*100.
3915 ! POLYSVP = 10.**(-7.90298*(373.16/T-1.)+ &
3916 ! 5.02808*LOG10(373.16/T)- &
3917 ! 1.3816E-7*(10**(11.344*(1.-T/373.16))-1.)+ &
3918 ! 8.1328E-3*(10**(-3.49149*(373.16/T-1.))-1.)+ &
3919 ! LOG10(1013.246))*100.
3921 END IF
3924 END FUNCTION POLYSVP
3926 !------------------------------------------------------------------------------
3928 REAL FUNCTION GAMMA(X)
3929 !----------------------------------------------------------------------
3930 !
3931 ! THIS ROUTINE CALCULATES THE GAMMA FUNCTION FOR A REAL ARGUMENT X.
3932 ! COMPUTATION IS BASED ON AN ALGORITHM OUTLINED IN REFERENCE 1.
3933 ! THE PROGRAM USES RATIONAL FUNCTIONS THAT APPROXIMATE THE GAMMA
3934 ! FUNCTION TO AT LEAST 20 SIGNIFICANT DECIMAL DIGITS. COEFFICIENTS
3935 ! FOR THE APPROXIMATION OVER THE INTERVAL (1,2) ARE UNPUBLISHED.
3936 ! THOSE FOR THE APPROXIMATION FOR X .GE. 12 ARE FROM REFERENCE 2.
3937 ! THE ACCURACY ACHIEVED DEPENDS ON THE ARITHMETIC SYSTEM, THE
3938 ! COMPILER, THE INTRINSIC FUNCTIONS, AND PROPER SELECTION OF THE
3939 ! MACHINE-DEPENDENT CONSTANTS.
3940 !
3941 !
3942 !*******************************************************************
3943 !*******************************************************************
3944 !
3945 ! EXPLANATION OF MACHINE-DEPENDENT CONSTANTS
3946 !
3947 ! BETA - RADIX FOR THE FLOATING-POINT REPRESENTATION
3948 ! MAXEXP - THE SMALLEST POSITIVE POWER OF BETA THAT OVERFLOWS
3949 ! XBIG - THE LARGEST ARGUMENT FOR WHICH GAMMA(X) IS REPRESENTABLE
3950 ! IN THE MACHINE, I.E., THE SOLUTION TO THE EQUATION
3951 ! GAMMA(XBIG) = BETA**MAXEXP
3952 ! XINF - THE LARGEST MACHINE REPRESENTABLE FLOATING-POINT NUMBER;
3953 ! APPROXIMATELY BETA**MAXEXP
3954 ! EPS - THE SMALLEST POSITIVE FLOATING-POINT NUMBER SUCH THAT
3955 ! 1.0+EPS .GT. 1.0
3956 ! XMININ - THE SMALLEST POSITIVE FLOATING-POINT NUMBER SUCH THAT
3957 ! 1/XMININ IS MACHINE REPRESENTABLE
3958 !
3959 ! APPROXIMATE VALUES FOR SOME IMPORTANT MACHINES ARE:
3960 !
3961 ! BETA MAXEXP XBIG
3962 !
3963 ! CRAY-1 (S.P.) 2 8191 966.961
3964 ! CYBER 180/855
3965 ! UNDER NOS (S.P.) 2 1070 177.803
3966 ! IEEE (IBM/XT,
3967 ! SUN, ETC.) (S.P.) 2 128 35.040
3968 ! IEEE (IBM/XT,
3969 ! SUN, ETC.) (D.P.) 2 1024 171.624
3970 ! IBM 3033 (D.P.) 16 63 57.574
3971 ! VAX D-FORMAT (D.P.) 2 127 34.844
3972 ! VAX G-FORMAT (D.P.) 2 1023 171.489
3973 !
3974 ! XINF EPS XMININ
3975 !
3976 ! CRAY-1 (S.P.) 5.45E+2465 7.11E-15 1.84E-2466
3977 ! CYBER 180/855
3978 ! UNDER NOS (S.P.) 1.26E+322 3.55E-15 3.14E-294
3979 ! IEEE (IBM/XT,
3980 ! SUN, ETC.) (S.P.) 3.40E+38 1.19E-7 1.18E-38
3981 ! IEEE (IBM/XT,
3982 ! SUN, ETC.) (D.P.) 1.79D+308 2.22D-16 2.23D-308
3983 ! IBM 3033 (D.P.) 7.23D+75 2.22D-16 1.39D-76
3984 ! VAX D-FORMAT (D.P.) 1.70D+38 1.39D-17 5.88D-39
3985 ! VAX G-FORMAT (D.P.) 8.98D+307 1.11D-16 1.12D-308
3986 !
3987 !*******************************************************************
3988 !*******************************************************************
3989 !
3990 ! ERROR RETURNS
3991 !
3992 ! THE PROGRAM RETURNS THE VALUE XINF FOR SINGULARITIES OR
3993 ! WHEN OVERFLOW WOULD OCCUR. THE COMPUTATION IS BELIEVED
3994 ! TO BE FREE OF UNDERFLOW AND OVERFLOW.
3995 !
3996 !
3997 ! INTRINSIC FUNCTIONS REQUIRED ARE:
3998 !
3999 ! INT, DBLE, EXP, LOG, REAL, SIN
4000 !
4001 !
4002 ! REFERENCES: AN OVERVIEW OF SOFTWARE DEVELOPMENT FOR SPECIAL
4003 ! FUNCTIONS W. J. CODY, LECTURE NOTES IN MATHEMATICS,
4004 ! 506, NUMERICAL ANALYSIS DUNDEE, 1975, G. A. WATSON
4005 ! (ED.), SPRINGER VERLAG, BERLIN, 1976.
4006 !
4007 ! COMPUTER APPROXIMATIONS, HART, ET. AL., WILEY AND
4008 ! SONS, NEW YORK, 1968.
4009 !
4010 ! LATEST MODIFICATION: OCTOBER 12, 1989
4011 !
4012 ! AUTHORS: W. J. CODY AND L. STOLTZ
4013 ! APPLIED MATHEMATICS DIVISION
4014 ! ARGONNE NATIONAL LABORATORY
4015 ! ARGONNE, IL 60439
4016 !
4017 !----------------------------------------------------------------------
4018 implicit none
4019 INTEGER I,N
4020 LOGICAL PARITY
4021 REAL &
4022 CONV,EPS,FACT,HALF,ONE,RES,SUM,TWELVE, &
4023 TWO,X,XBIG,XDEN,XINF,XMININ,XNUM,Y,Y1,YSQ,Z,ZERO
4024 REAL, DIMENSION(7) :: C
4025 REAL, DIMENSION(8) :: P
4026 REAL, DIMENSION(8) :: Q
4027 !----------------------------------------------------------------------
4028 ! MATHEMATICAL CONSTANTS
4029 !----------------------------------------------------------------------
4030 DATA ONE,HALF,TWELVE,TWO,ZERO/1.0E0,0.5E0,12.0E0,2.0E0,0.0E0/
4033 !----------------------------------------------------------------------
4034 ! MACHINE DEPENDENT PARAMETERS
4035 !----------------------------------------------------------------------
4036 DATA XBIG,XMININ,EPS/35.040E0,1.18E-38,1.19E-7/,XINF/3.4E38/
4037 !----------------------------------------------------------------------
4038 ! NUMERATOR AND DENOMINATOR COEFFICIENTS FOR RATIONAL MINIMAX
4039 ! APPROXIMATION OVER (1,2).
4040 !----------------------------------------------------------------------
4041 DATA P/-1.71618513886549492533811E+0,2.47656508055759199108314E+1, &
4042 -3.79804256470945635097577E+2,6.29331155312818442661052E+2, &
4043 8.66966202790413211295064E+2,-3.14512729688483675254357E+4, &
4044 -3.61444134186911729807069E+4,6.64561438202405440627855E+4/
4045 DATA Q/-3.08402300119738975254353E+1,3.15350626979604161529144E+2, &
4046 -1.01515636749021914166146E+3,-3.10777167157231109440444E+3, &
4047 2.25381184209801510330112E+4,4.75584627752788110767815E+3, &
4048 -1.34659959864969306392456E+5,-1.15132259675553483497211E+5/
4049 !----------------------------------------------------------------------
4050 ! COEFFICIENTS FOR MINIMAX APPROXIMATION OVER (12, INF).
4051 !----------------------------------------------------------------------
4052 DATA C/-1.910444077728E-03,8.4171387781295E-04, &
4053 -5.952379913043012E-04,7.93650793500350248E-04, &
4054 -2.777777777777681622553E-03,8.333333333333333331554247E-02, &
4055 5.7083835261E-03/
4056 !----------------------------------------------------------------------
4057 ! STATEMENT FUNCTIONS FOR CONVERSION BETWEEN INTEGER AND FLOAT
4058 !----------------------------------------------------------------------
4059 CONV(I) = REAL(I)
4060 PARITY=.FALSE.
4061 FACT=ONE
4062 N=0
4063 Y=X
4064 IF(Y.LE.ZERO)THEN
4065 !----------------------------------------------------------------------
4066 ! ARGUMENT IS NEGATIVE
4067 !----------------------------------------------------------------------
4068 Y=-X
4069 Y1=AINT(Y)
4070 RES=Y-Y1
4071 IF(RES.NE.ZERO)THEN
4072 IF(Y1.NE.AINT(Y1*HALF)*TWO)PARITY=.TRUE.
4073 FACT=-PI/SIN(PI*RES)
4074 Y=Y+ONE
4075 ELSE
4076 RES=XINF
4077 GOTO 900
4078 ENDIF
4079 ENDIF
4080 !----------------------------------------------------------------------
4081 ! ARGUMENT IS POSITIVE
4082 !----------------------------------------------------------------------
4083 IF(Y.LT.EPS)THEN
4084 !----------------------------------------------------------------------
4085 ! ARGUMENT .LT. EPS
4086 !----------------------------------------------------------------------
4087 IF(Y.GE.XMININ)THEN
4088 RES=ONE/Y
4089 ELSE
4090 RES=XINF
4091 GOTO 900
4092 ENDIF
4093 ELSEIF(Y.LT.TWELVE)THEN
4094 Y1=Y
4095 IF(Y.LT.ONE)THEN
4096 !----------------------------------------------------------------------
4097 ! 0.0 .LT. ARGUMENT .LT. 1.0
4098 !----------------------------------------------------------------------
4099 Z=Y
4100 Y=Y+ONE
4101 ELSE
4102 !----------------------------------------------------------------------
4103 ! 1.0 .LT. ARGUMENT .LT. 12.0, REDUCE ARGUMENT IF NECESSARY
4104 !----------------------------------------------------------------------
4105 N=INT(Y)-1
4106 Y=Y-CONV(N)
4107 Z=Y-ONE
4108 ENDIF
4109 !----------------------------------------------------------------------
4110 ! EVALUATE APPROXIMATION FOR 1.0 .LT. ARGUMENT .LT. 2.0
4111 !----------------------------------------------------------------------
4112 XNUM=ZERO
4113 XDEN=ONE
4114 DO I=1,8
4115 XNUM=(XNUM+P(I))*Z
4116 XDEN=XDEN*Z+Q(I)
4117 END DO
4118 RES=XNUM/XDEN+ONE
4119 IF(Y1.LT.Y)THEN
4120 !----------------------------------------------------------------------
4121 ! ADJUST RESULT FOR CASE 0.0 .LT. ARGUMENT .LT. 1.0
4122 !----------------------------------------------------------------------
4123 RES=RES/Y1
4124 ELSEIF(Y1.GT.Y)THEN
4125 !----------------------------------------------------------------------
4126 ! ADJUST RESULT FOR CASE 2.0 .LT. ARGUMENT .LT. 12.0
4127 !----------------------------------------------------------------------
4128 DO I=1,N
4129 RES=RES*Y
4130 Y=Y+ONE
4131 END DO
4132 ENDIF
4133 ELSE
4134 !----------------------------------------------------------------------
4135 ! EVALUATE FOR ARGUMENT .GE. 12.0,
4136 !----------------------------------------------------------------------
4137 IF(Y.LE.XBIG)THEN
4138 YSQ=Y*Y
4139 SUM=C(7)
4140 DO I=1,6
4141 SUM=SUM/YSQ+C(I)
4142 END DO
4143 SUM=SUM/Y-Y+SQRTPI
4144 SUM=SUM+(Y-HALF)*LOG(Y)
4145 RES=EXP(SUM)
4146 ELSE
4147 RES=XINF
4148 GOTO 900
4149 ENDIF
4150 ENDIF
4151 !----------------------------------------------------------------------
4152 ! FINAL ADJUSTMENTS AND RETURN
4153 !----------------------------------------------------------------------
4154 IF(PARITY)RES=-RES
4155 IF(FACT.NE.ONE)RES=FACT/RES
4156 900 GAMMA=RES
4157 RETURN
4158 ! ---------- LAST LINE OF GAMMA ----------
4159 END FUNCTION GAMMA
4162 REAL FUNCTION DERF1(X)
4163 IMPLICIT NONE
4164 REAL X
4165 REAL, DIMENSION(0 : 64) :: A, B
4166 REAL W,T,Y
4167 INTEGER K,I
4168 DATA A/ &
4169 0.00000000005958930743E0, -0.00000000113739022964E0, &
4170 0.00000001466005199839E0, -0.00000016350354461960E0, &
4171 0.00000164610044809620E0, -0.00001492559551950604E0, &
4172 0.00012055331122299265E0, -0.00085483269811296660E0, &
4173 0.00522397762482322257E0, -0.02686617064507733420E0, &
4174 0.11283791670954881569E0, -0.37612638903183748117E0, &
4175 1.12837916709551257377E0, &
4176 0.00000000002372510631E0, -0.00000000045493253732E0, &
4177 0.00000000590362766598E0, -0.00000006642090827576E0, &
4178 0.00000067595634268133E0, -0.00000621188515924000E0, &
4179 0.00005103883009709690E0, -0.00037015410692956173E0, &
4180 0.00233307631218880978E0, -0.01254988477182192210E0, &
4181 0.05657061146827041994E0, -0.21379664776456006580E0, &
4182 0.84270079294971486929E0, &
4183 0.00000000000949905026E0, -0.00000000018310229805E0, &
4184 0.00000000239463074000E0, -0.00000002721444369609E0, &
4185 0.00000028045522331686E0, -0.00000261830022482897E0, &
4186 0.00002195455056768781E0, -0.00016358986921372656E0, &
4187 0.00107052153564110318E0, -0.00608284718113590151E0, &
4188 0.02986978465246258244E0, -0.13055593046562267625E0, &
4189 0.67493323603965504676E0, &
4190 0.00000000000382722073E0, -0.00000000007421598602E0, &
4191 0.00000000097930574080E0, -0.00000001126008898854E0, &
4192 0.00000011775134830784E0, -0.00000111992758382650E0, &
4193 0.00000962023443095201E0, -0.00007404402135070773E0, &
4194 0.00050689993654144881E0, -0.00307553051439272889E0, &
4195 0.01668977892553165586E0, -0.08548534594781312114E0, &
4196 0.56909076642393639985E0, &
4197 0.00000000000155296588E0, -0.00000000003032205868E0, &
4198 0.00000000040424830707E0, -0.00000000471135111493E0, &
4199 0.00000005011915876293E0, -0.00000048722516178974E0, &
4200 0.00000430683284629395E0, -0.00003445026145385764E0, &
4201 0.00024879276133931664E0, -0.00162940941748079288E0, &
4202 0.00988786373932350462E0, -0.05962426839442303805E0, &
4203 0.49766113250947636708E0 /
4204 DATA (B(I), I = 0, 12) / &
4205 -0.00000000029734388465E0, 0.00000000269776334046E0, &
4206 -0.00000000640788827665E0, -0.00000001667820132100E0, &
4207 -0.00000021854388148686E0, 0.00000266246030457984E0, &
4208 0.00001612722157047886E0, -0.00025616361025506629E0, &
4209 0.00015380842432375365E0, 0.00815533022524927908E0, &
4210 -0.01402283663896319337E0, -0.19746892495383021487E0, &
4211 0.71511720328842845913E0 /
4212 DATA (B(I), I = 13, 25) / &
4213 -0.00000000001951073787E0, -0.00000000032302692214E0, &
4214 0.00000000522461866919E0, 0.00000000342940918551E0, &
4215 -0.00000035772874310272E0, 0.00000019999935792654E0, &
4216 0.00002687044575042908E0, -0.00011843240273775776E0, &
4217 -0.00080991728956032271E0, 0.00661062970502241174E0, &
4218 0.00909530922354827295E0, -0.20160072778491013140E0, &
4219 0.51169696718727644908E0 /
4220 DATA (B(I), I = 26, 38) / &
4221 0.00000000003147682272E0, -0.00000000048465972408E0, &
4222 0.00000000063675740242E0, 0.00000003377623323271E0, &
4223 -0.00000015451139637086E0, -0.00000203340624738438E0, &
4224 0.00001947204525295057E0, 0.00002854147231653228E0, &
4225 -0.00101565063152200272E0, 0.00271187003520095655E0, &
4226 0.02328095035422810727E0, -0.16725021123116877197E0, &
4227 0.32490054966649436974E0 /
4228 DATA (B(I), I = 39, 51) / &
4229 0.00000000002319363370E0, -0.00000000006303206648E0, &
4230 -0.00000000264888267434E0, 0.00000002050708040581E0, &
4231 0.00000011371857327578E0, -0.00000211211337219663E0, &
4232 0.00000368797328322935E0, 0.00009823686253424796E0, &
4233 -0.00065860243990455368E0, -0.00075285814895230877E0, &
4234 0.02585434424202960464E0, -0.11637092784486193258E0, &
4235 0.18267336775296612024E0 /
4236 DATA (B(I), I = 52, 64) / &
4237 -0.00000000000367789363E0, 0.00000000020876046746E0, &
4238 -0.00000000193319027226E0, -0.00000000435953392472E0, &
4239 0.00000018006992266137E0, -0.00000078441223763969E0, &
4240 -0.00000675407647949153E0, 0.00008428418334440096E0, &
4241 -0.00017604388937031815E0, -0.00239729611435071610E0, &
4242 0.02064129023876022970E0, -0.06905562880005864105E0, &
4243 0.09084526782065478489E0 /
4244 W = ABS(X)
4245 IF (W .LT. 2.2D0) THEN
4246 T = W * W
4247 K = INT(T)
4248 T = T - K
4249 K = K * 13
4250 Y = ((((((((((((A(K) * T + A(K + 1)) * T + &
4251 A(K + 2)) * T + A(K + 3)) * T + A(K + 4)) * T + &
4252 A(K + 5)) * T + A(K + 6)) * T + A(K + 7)) * T + &
4253 A(K + 8)) * T + A(K + 9)) * T + A(K + 10)) * T + &
4254 A(K + 11)) * T + A(K + 12)) * W
4255 ELSE IF (W .LT. 6.9D0) THEN
4256 K = INT(W)
4257 T = W - K
4258 K = 13 * (K - 2)
4259 Y = (((((((((((B(K) * T + B(K + 1)) * T + &
4260 B(K + 2)) * T + B(K + 3)) * T + B(K + 4)) * T + &
4261 B(K + 5)) * T + B(K + 6)) * T + B(K + 7)) * T + &
4262 B(K + 8)) * T + B(K + 9)) * T + B(K + 10)) * T + &
4263 B(K + 11)) * T + B(K + 12)
4264 Y = Y * Y
4265 Y = Y * Y
4266 Y = Y * Y
4267 Y = 1 - Y * Y
4268 ELSE
4269 Y = 1
4270 END IF
4271 IF (X .LT. 0) Y = -Y
4272 DERF1 = Y
4273 END FUNCTION DERF1
4275 !+---+-----------------------------------------------------------------+
4277 END MODULE module_mp_morr_two_moment