| blob | a225395ee3e02b5b450fa7d6a9392a82724f3385 |
1 !WRF:MODEL_LAYER:PHYSICS
2 !
4 ! THIS MODULE CONTAINS THE TWO-MOMENT MICROPHYSICS CODE DESCRIBED BY
5 ! MORRISON ET AL. 2005 JAS; MORRISON AND PINTO 2005 JAS.
6 ! ADDITIONAL CHANGES ARE DESCRIBED IN DETAIL BY MORRISON, THOMPSON, TATARSKII (MWR, SUBMITTED)
8 ! THIS SCHEME IS A BULK DOUBLE-MOMENT SCHEME THAT PREDICTS MIXING
9 ! RATIOS AND NUMBER CONCENTRATIONS OF FIVE HYDROMETEOR SPECIES:
10 ! CLOUD DROPLETS, CLOUD (SMALL) ICE, RAIN, SNOW, AND GRAUPEL.
12 MODULE MODULE_MP_MORR_TWO_MOMENT
13 USE module_wrf_error
14 ! USE module_utility, ONLY: WRFU_Clock, WRFU_Alarm ! GT
15 ! USE module_domain, ONLY : HISTORY_ALARM, Is_alarm_tstep ! GT
17 ! USE module_state_description
19 IMPLICIT NONE
21 REAL, PARAMETER :: PI = 3.1415926535897932384626434
22 REAL, PARAMETER :: SQRTPI = 0.9189385332046727417803297
24 PUBLIC :: MP_MORR_TWO_MOMENT
25 PUBLIC :: POLYSVP
27 PRIVATE :: GAMMA, DERF1
28 PRIVATE :: PI, SQRTPI
29 PRIVATE :: MORR_TWO_MOMENT_MICRO
31 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
32 ! SWITCHES FOR MICROPHYSICS SCHEME
33 ! IACT = 1, USE POWER-LAW CCN SPECTRA, NCCN = CS^K
34 ! IACT = 2, USE LOGNORMAL AEROSOL SIZE DIST TO DERIVE CCN SPECTRA
36 INTEGER, PRIVATE :: IACT
38 ! INUM = 0, PREDICT DROPLET CONCENTRATION
39 ! INUM = 1, ASSUME CONSTANT DROPLET CONCENTRATION
40 ! !!!NOTE: PREDICTED DROPLET CONCENTRATION NOT AVAILABLE IN THIS VERSION
41 ! CONTACT HUGH MORRISON (morrison@ucar.edu) FOR FURTHER INFORMATION
43 INTEGER, PRIVATE :: INUM
45 ! FOR INUM = 1, SET CONSTANT DROPLET CONCENTRATION (CM-3)
46 REAL, PRIVATE :: NDCNST
48 ! SWITCH FOR LIQUID-ONLY RUN
49 ! ILIQ = 0, INCLUDE ICE
50 ! ILIQ = 1, LIQUID ONLY, NO ICE
52 INTEGER, PRIVATE :: ILIQ
54 ! SWITCH FOR ICE NUCLEATION
55 ! INUC = 0, USE FORMULA FROM RASMUSSEN ET AL. 2002 (MID-LATITUDE)
56 ! = 1, USE MPACE OBSERVATIONS
58 INTEGER, PRIVATE :: INUC
60 ! IBASE = 1, NEGLECT DROPLET ACTIVATION AT LATERAL CLOUD EDGES DUE TO
61 ! UNRESOLVED ENTRAINMENT AND MIXING, ACTIVATE
62 ! AT CLOUD BASE OR IN REGION WITH LITTLE CLOUD WATER USING
63 ! NON-EQULIBRIUM SUPERSATURATION,
64 ! IN CLOUD INTERIOR ACTIVATE USING EQUILIBRIUM SUPERSATURATION
65 ! IBASE = 2, ASSUME DROPLET ACTIVATION AT LATERAL CLOUD EDGES DUE TO
66 ! UNRESOLVED ENTRAINMENT AND MIXING DOMINATES,
67 ! ACTIVATE DROPLETS EVERYWHERE IN THE CLOUD USING NON-EQUILIBRIUM
68 ! SUPERSATURATION, BASED ON THE
69 ! LOCAL SUB-GRID AND/OR GRID-SCALE VERTICAL VELOCITY
70 ! AT THE GRID POINT
72 ! NOTE: ONLY USED FOR PREDICTED DROPLET CONCENTRATION (INUM = 0)
74 INTEGER, PRIVATE :: IBASE
76 ! INCLUDE SUB-GRID VERTICAL VELOCITY IN DROPLET ACTIVATION
77 ! ISUB = 0, INCLUDE SUB-GRID W (RECOMMENDED FOR LOWER RESOLUTION)
78 ! ISUB = 1, EXCLUDE SUB-GRID W, ONLY USE GRID-SCALE W
80 INTEGER, PRIVATE :: ISUB
82 ! SWITCH FOR GRAUPEL/NO GRAUPEL
83 ! IGRAUP = 0, INCLUDE GRAUPEL
84 ! IGRAUP = 1, NO GRAUPEL
86 INTEGER, PRIVATE :: IGRAUP
88 ! HM ADDED NEW OPTION FOR HAIL
89 ! SWITCH FOR HAIL/GRAUPEL
90 ! IHAIL = 0, DENSE PRECIPITATING ICE IS GRAUPEL
91 ! IHAIL = 1, DENSE PRECIPITATING GICE IS HAIL
93 INTEGER, PRIVATE :: IHAIL
95 ! CLOUD MICROPHYSICS CONSTANTS
97 REAL, PRIVATE :: AI,AC,AS,AR,AG ! 'A' PARAMETER IN FALLSPEED-DIAM RELATIONSHIP
98 REAL, PRIVATE :: BI,BC,BS,BR,BG ! 'B' PARAMETER IN FALLSPEED-DIAM RELATIONSHIP
99 REAL, PRIVATE :: R ! GAS CONSTANT FOR AIR
100 REAL, PRIVATE :: RV ! GAS CONSTANT FOR WATER VAPOR
101 REAL, PRIVATE :: CP ! SPECIFIC HEAT AT CONSTANT PRESSURE FOR DRY AIR
102 REAL, PRIVATE :: RHOSU ! STANDARD AIR DENSITY AT 850 MB
103 REAL, PRIVATE :: RHOW ! DENSITY OF LIQUID WATER
104 REAL, PRIVATE :: RHOI ! BULK DENSITY OF CLOUD ICE
105 REAL, PRIVATE :: RHOSN ! BULK DENSITY OF SNOW
106 REAL, PRIVATE :: RHOG ! BULK DENSITY OF GRAUPEL
107 REAL, PRIVATE :: AIMM ! PARAMETER IN BIGG IMMERSION FREEZING
108 REAL, PRIVATE :: BIMM ! PARAMETER IN BIGG IMMERSION FREEZING
109 REAL, PRIVATE :: ECR ! COLLECTION EFFICIENCY BETWEEN DROPLETS/RAIN AND SNOW/RAIN
110 REAL, PRIVATE :: DCS ! THRESHOLD SIZE FOR CLOUD ICE AUTOCONVERSION
111 REAL, PRIVATE :: MI0 ! INITIAL SIZE OF NUCLEATED CRYSTAL
112 REAL, PRIVATE :: MG0 ! MASS OF EMBRYO GRAUPEL
113 REAL, PRIVATE :: F1S ! VENTILATION PARAMETER FOR SNOW
114 REAL, PRIVATE :: F2S ! VENTILATION PARAMETER FOR SNOW
115 REAL, PRIVATE :: F1R ! VENTILATION PARAMETER FOR RAIN
116 REAL, PRIVATE :: F2R ! VENTILATION PARAMETER FOR RAIN
117 REAL, PRIVATE :: G ! GRAVITATIONAL ACCELERATION
118 REAL, PRIVATE :: QSMALL ! SMALLEST ALLOWED HYDROMETEOR MIXING RATIO
119 REAL, PRIVATE :: CI,DI,CS,DS,CG,DG ! SIZE DISTRIBUTION PARAMETERS FOR CLOUD ICE, SNOW, GRAUPEL
120 REAL, PRIVATE :: EII ! COLLECTION EFFICIENCY, ICE-ICE COLLISIONS
121 REAL, PRIVATE :: ECI ! COLLECTION EFFICIENCY, ICE-DROPLET COLLISIONS
122 REAL, PRIVATE :: RIN ! RADIUS OF CONTACT NUCLEI (M)
124 ! CCN SPECTRA FOR IACT = 1
126 REAL, PRIVATE :: C1 ! 'C' IN NCCN = CS^K (CM-3)
127 REAL, PRIVATE :: K1 ! 'K' IN NCCN = CS^K
129 ! AEROSOL PARAMETERS FOR IACT = 2
131 REAL, PRIVATE :: MW ! MOLECULAR WEIGHT WATER (KG/MOL)
132 REAL, PRIVATE :: OSM ! OSMOTIC COEFFICIENT
133 REAL, PRIVATE :: VI ! NUMBER OF ION DISSOCIATED IN SOLUTION
134 REAL, PRIVATE :: EPSM ! AEROSOL SOLUBLE FRACTION
135 REAL, PRIVATE :: RHOA ! AEROSOL BULK DENSITY (KG/M3)
136 REAL, PRIVATE :: MAP ! MOLECULAR WEIGHT AEROSOL (KG/MOL)
137 REAL, PRIVATE :: MA ! MOLECULAR WEIGHT OF 'AIR' (KG/MOL)
138 REAL, PRIVATE :: RR ! UNIVERSAL GAS CONSTANT
139 REAL, PRIVATE :: BACT ! ACTIVATION PARAMETER
140 REAL, PRIVATE :: RM1 ! GEOMETRIC MEAN RADIUS, MODE 1 (M)
141 REAL, PRIVATE :: RM2 ! GEOMETRIC MEAN RADIUS, MODE 2 (M)
142 REAL, PRIVATE :: NANEW1 ! TOTAL AEROSOL CONCENTRATION, MODE 1 (M^-3)
143 REAL, PRIVATE :: NANEW2 ! TOTAL AEROSOL CONCENTRATION, MODE 2 (M^-3)
144 REAL, PRIVATE :: SIG1 ! STANDARD DEVIATION OF AEROSOL S.D., MODE 1
145 REAL, PRIVATE :: SIG2 ! STANDARD DEVIATION OF AEROSOL S.D., MODE 2
146 REAL, PRIVATE :: F11 ! CORRECTION FACTOR FOR ACTIVATION, MODE 1
147 REAL, PRIVATE :: F12 ! CORRECTION FACTOR FOR ACTIVATION, MODE 1
148 REAL, PRIVATE :: F21 ! CORRECTION FACTOR FOR ACTIVATION, MODE 2
149 REAL, PRIVATE :: F22 ! CORRECTION FACTOR FOR ACTIVATION, MODE 2
150 REAL, PRIVATE :: MMULT ! MASS OF SPLINTERED ICE PARTICLE
151 REAL, PRIVATE :: LAMMAXI,LAMMINI,LAMMAXR,LAMMINR,LAMMAXS,LAMMINS,LAMMAXG,LAMMING
153 ! CONSTANTS TO IMPROVE EFFICIENCY
155 REAL, PRIVATE :: CONS1,CONS2,CONS3,CONS4,CONS5,CONS6,CONS7,CONS8,CONS9,CONS10
156 REAL, PRIVATE :: CONS11,CONS12,CONS13,CONS14,CONS15,CONS16,CONS17,CONS18,CONS19,CONS20
157 REAL, PRIVATE :: CONS21,CONS22,CONS23,CONS24,CONS25,CONS26,CONS27,CONS28,CONS29,CONS30
158 REAL, PRIVATE :: CONS31,CONS32,CONS33,CONS34,CONS35,CONS36,CONS37,CONS38,CONS39,CONS40
159 REAL, PRIVATE :: CONS41
161 CONTAINS
163 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
164 SUBROUTINE MORR_TWO_MOMENT_INIT
165 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
166 ! THIS SUBROUTINE INITIALIZES ALL PHYSICAL CONSTANTS AMND PARAMETERS
167 ! NEEDED BY THE MICROPHYSICS SCHEME.
168 ! NEEDS TO BE CALLED AT FIRST TIME STEP, PRIOR TO CALL TO MAIN MICROPHYSICS INTERFACE
169 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
171 IMPLICIT NONE
173 integer n,i
175 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
176 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
178 ! THE FOLLOWING PARAMETERS ARE USER-DEFINED SWITCHES AND NEED TO BE
179 ! SET PRIOR TO CODE COMPILATION
181 ! INUM = 0, PREDICT DROPLET CONCENTRATION
182 ! INUM = 1, ASSUME CONSTANT DROPLET CONCENTRATION
183 ! !!!NOTE: PREDICTED DROPLET CONCENTRATION NOT AVAILABLE IN THIS VERSION
184 ! CONTACT HUGH MORRISON (morrison@ucar.edu) FOR FURTHER INFORMATION
185 ! INUM=1 ONLY IN CURRENT VERSION
187 INUM = 1
189 ! FOR INUM = 1, SET CONSTANT DROPLET CONCENTRATION (UNITS OF CM-3)
191 NDCNST = 250.
193 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
194 ! NOTE, FOLLOWING OPTIONS NOT AVAILABLE IN CURRENT VERSION
195 ! ONLY USED WHEN INUM=0
198 ! IACT = 1, USE POWER-LAW CCN SPECTRA, NCCN = CS^K
199 ! IACT = 2, USE LOGNORMAL AEROSOL SIZE DIST TO DERIVE CCN SPECTRA
200 ! NOTE: ONLY USED FOR PREDICTED DROPLET CONCENTRATION (INUM = 0)
202 IACT = 2
204 ! IBASE = 1, NEGLECT DROPLET ACTIVATION AT LATERAL CLOUD EDGES DUE TO
205 ! UNRESOLVED ENTRAINMENT AND MIXING, ACTIVATE
206 ! AT CLOUD BASE OR IN REGION WITH LITTLE CLOUD WATER USING
207 ! NON-EQULIBRIUM SUPERSATURATION ASSUMING NO INITIAL CLOUD WATER,
208 ! IN CLOUD INTERIOR ACTIVATE USING EQUILIBRIUM SUPERSATURATION
209 ! IBASE = 2, ASSUME DROPLET ACTIVATION AT LATERAL CLOUD EDGES DUE TO
210 ! UNRESOLVED ENTRAINMENT AND MIXING DOMINATES,
211 ! ACTIVATE DROPLETS EVERYWHERE IN THE CLOUD USING NON-EQUILIBRIUM
212 ! SUPERSATURATION ASSUMING NO INITIAL CLOUD WATER, BASED ON THE
213 ! LOCAL SUB-GRID AND/OR GRID-SCALE VERTICAL VELOCITY
214 ! AT THE GRID POINT
216 ! NOTE: ONLY USED FOR PREDICTED DROPLET CONCENTRATION (INUM = 0)
218 IBASE = 2
220 ! INCLUDE SUB-GRID VERTICAL VELOCITY (standard deviation of w) IN DROPLET ACTIVATION
221 ! ISUB = 0, INCLUDE SUB-GRID W (RECOMMENDED FOR LOWER RESOLUTION)
222 ! currently, sub-grid w is constant of 0.5 m/s (not coupled with PBL/turbulence scheme)
223 ! ISUB = 1, EXCLUDE SUB-GRID W, ONLY USE GRID-SCALE W
225 ! NOTE: ONLY USED FOR PREDICTED DROPLET CONCENTRATION (INUM = 0)
227 ISUB = 0
228 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
231 ! SWITCH FOR LIQUID-ONLY RUN
232 ! ILIQ = 0, INCLUDE ICE
233 ! ILIQ = 1, LIQUID ONLY, NO ICE
235 ILIQ = 0
237 ! SWITCH FOR ICE NUCLEATION
238 ! INUC = 0, USE FORMULA FROM RASMUSSEN ET AL. 2002 (MID-LATITUDE)
239 ! = 1, USE MPACE OBSERVATIONS (ARCTIC ONLY)
241 INUC = 0
243 ! SWITCH FOR GRAUPEL/HAIL NO GRAUPEL/HAIL
244 ! IGRAUP = 0, INCLUDE GRAUPEL/HAIL
245 ! IGRAUP = 1, NO GRAUPEL/HAIL
247 IGRAUP = 0
249 ! HM ADDED 11/7/07
250 ! SWITCH FOR HAIL/GRAUPEL
251 ! IHAIL = 0, DENSE PRECIPITATING ICE IS GRAUPEL
252 ! IHAIL = 1, DENSE PRECIPITATING ICE IS HAIL
253 ! NOTE ---> RECOMMEND IHAIL = 1 FOR CONTINENTAL DEEP CONVECTION
255 IHAIL = 0
257 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
258 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
259 ! SET PHYSICAL CONSTANTS
261 ! FALLSPEED PARAMETERS (V=AD^B)
262 AI = 700.
263 AC = 3.E7
264 AS = 11.72
265 AR = 841.99667
266 BI = 1.
267 BC = 2.
268 BS = 0.41
269 BR = 0.8
270 IF (IHAIL.EQ.0) THEN
271 AG = 19.3
272 BG = 0.37
273 ELSE ! (MATSUN AND HUGGINS 1980)
274 AG = 114.5
275 BG = 0.5
276 END IF
278 ! CONSTANTS AND PARAMETERS
279 R = 287.15
280 RV = 461.5
281 CP = 1005.
282 RHOSU = 85000./(287.15*273.15)
283 RHOW = 997.
284 RHOI = 500.
285 RHOSN = 100.
286 IF (IHAIL.EQ.0) THEN
287 RHOG = 400.
288 ELSE
289 RHOG = 900.
290 END IF
291 AIMM = 0.66
292 BIMM = 100.
293 ECR = 1.
294 DCS = 125.E-6
295 MI0 = 4./3.*PI*RHOI*(10.E-6)**3
296 MG0 = 1.6E-10
297 F1S = 0.86
298 F2S = 0.28
299 F1R = 0.78
300 F2R = 0.32
301 G = 9.806
302 QSMALL = 1.E-14
303 EII = 0.1
304 ECI = 0.7
306 ! SIZE DISTRIBUTION PARAMETERS
308 CI = RHOI*PI/6.
309 DI = 3.
310 CS = RHOSN*PI/6.
311 DS = 3.
312 CG = RHOG*PI/6.
313 DG = 3.
315 ! RADIUS OF CONTACT NUCLEI
316 RIN = 0.1E-6
318 MMULT = 4./3.*PI*RHOI*(5.E-6)**3
320 ! SIZE LIMITS FOR LAMBDA
322 LAMMAXI = 1./1.E-6
323 LAMMINI = 1./(2.*DCS+100.E-6)
324 LAMMAXR = 1./20.E-6
325 LAMMINR = 1./500.E-6
326 LAMMAXS = 1./10.E-6
327 LAMMINS = 1./2000.E-6
328 LAMMAXG = 1./20.E-6
329 LAMMING = 1./2000.E-6
331 ! CCN SPECTRA FOR IACT = 1
333 ! MARITIME
334 ! MODIFIED FROM RASMUSSEN ET AL. 2002
335 ! NCCN = C*S^K, NCCN IS IN CM-3, S IS SUPERSATURATION RATIO IN %
337 K1 = 0.4
338 C1 = 120.
340 ! CONTINENTAL
342 ! K1 = 0.5
343 ! C1 = 1000.
345 ! AEROSOL ACTIVATION PARAMETERS FOR IACT = 2
346 ! PARAMETERS CURRENTLY SET FOR AMMONIUM SULFATE
348 MW = 0.018
349 OSM = 1.
350 VI = 3.
351 EPSM = 0.7
352 RHOA = 1777.
353 MAP = 0.132
354 MA = 0.0284
355 RR = 8.3187
356 BACT = VI*OSM*EPSM*MW*RHOA/(MAP*RHOW)
358 ! AEROSOL SIZE DISTRIBUTION PARAMETERS CURRENTLY SET FOR MPACE
359 ! (see morrison et al. 2007, JGR)
360 ! MODE 1
362 RM1 = 0.052E-6
363 SIG1 = 2.04
364 NANEW1 = 72.2E6
365 F11 = 0.5*EXP(2.5*(LOG(SIG1))**2)
366 F21 = 1.+0.25*LOG(SIG1)
368 ! MODE 2
370 RM2 = 1.3E-6
371 SIG2 = 2.5
372 NANEW2 = 1.8E6
373 F12 = 0.5*EXP(2.5*(LOG(SIG2))**2)
374 F22 = 1.+0.25*LOG(SIG2)
376 ! CONSTANTS FOR EFFICIENCY
378 CONS1=GAMMA(1.+DS)*CS
379 CONS2=GAMMA(1.+DG)*CG
380 CONS3=GAMMA(4.+BS)/6.
381 CONS4=GAMMA(4.+BR)/6.
382 CONS5=GAMMA(1.+BS)
383 CONS6=GAMMA(1.+BR)
384 CONS7=GAMMA(4.+BG)/6.
385 CONS8=GAMMA(1.+BG)
386 CONS9=GAMMA(5./2.+BR/2.)
387 CONS10=GAMMA(5./2.+BS/2.)
388 CONS11=GAMMA(5./2.+BG/2.)
389 CONS12=GAMMA(1.+DI)*CI
390 CONS13=GAMMA(BS+3.)*PI/4.*ECI
391 CONS14=GAMMA(BG+3.)*PI/4.*ECI
392 CONS15=-1108.*EII*PI**((1.-BS)/3.)*RHOSN**((-2.-BS)/3.)/(4.*720.)
393 CONS16=GAMMA(BI+3.)*PI/4.*ECI
394 CONS17=4.*2.*3.*RHOSU*PI*ECI*ECI*GAMMA(2.*BS+2.)/(8.*(RHOG-RHOSN))
395 CONS18=RHOSN*RHOSN
396 CONS19=RHOW*RHOW
397 CONS20=20.*PI*PI*RHOW*BIMM
398 CONS21=4./(DCS*RHOI)
399 CONS22=PI*RHOI*DCS**3/6.
400 CONS23=PI/4.*EII*GAMMA(BS+3.)
401 CONS24=PI/4.*ECR*GAMMA(BR+3.)
402 CONS25=PI*PI/24.*RHOW*ECR*GAMMA(BR+6.)
403 CONS26=PI/6.*RHOW
404 CONS27=GAMMA(1.+BI)
405 CONS28=GAMMA(4.+BI)/6.
406 CONS29=4./3.*PI*RHOW*(25.E-6)**3
407 CONS30=4./3.*PI*RHOW
408 CONS31=PI*PI*ECR*RHOSN
409 CONS32=PI/2.*ECR
410 CONS33=PI*PI*ECR*RHOG
411 CONS34=5./2.+BR/2.
412 CONS35=5./2.+BS/2.
413 CONS36=5./2.+BG/2.
414 CONS37=4.*PI*1.38E-23/(6.*PI*RIN)
415 CONS38=PI*PI/3.*RHOW
416 CONS39=PI*PI/36.*RHOW*BIMM
417 CONS40=PI/6.*BIMM
418 CONS41=PI*PI*ECR*RHOW
420 END SUBROUTINE MORR_TWO_MOMENT_INIT
422 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
423 ! THIS SUBROUTINE IS MAIN INTERFACE WITH THE TWO-MOMENT MICROPHYSICS SCHEME
424 ! THIS INTERFACE TAKES IN 3D VARIABLES FROM DRIVER MODEL, CONVERTS TO 1D FOR
425 ! CALL TO THE MAIN MICROPHYSICS SUBROUTINE (SUBROUTINE MORR_TWO_MOMENT_MICRO)
426 ! WHICH OPERATES ON 1D VERTICAL COLUMNS.
427 ! 1D VARIABLES FROM THE MAIN MICROPHYSICS SUBROUTINE ARE THEN REASSIGNED BACK TO 3D FOR OUTPUT
428 ! BACK TO DRIVER MODEL USING THIS INTERFACE.
429 ! MICROPHYSICS TENDENCIES ARE ADDED TO VARIABLES HERE BEFORE BEING PASSED BACK TO DRIVER MODEL.
431 ! THIS CODE WAS WRITTEN BY HUGH MORRISON (NCAR) AND SLAVA TATARSKII (GEORGIA TECH).
433 ! FOR QUESTIONS, CONTACT: HUGH MORRISON, E-MAIL: MORRISON@UCAR.EDU, PHONE:303-497-8916
435 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
437 SUBROUTINE MP_MORR_TWO_MOMENT(ITIMESTEP, &
438 TH, QV, QC, QR, QI, QS, QG, NI, NS, NR, NG, &
439 RHO, PII, P, DT_IN, DZ, HT, W, &
440 RAINNC, RAINNCV, SR, &
441 qrcuten, qscuten, qicuten, mu & ! hm added
442 ,IDS,IDE, JDS,JDE, KDS,KDE & ! domain dims
443 ,IMS,IME, JMS,JME, KMS,KME & ! memory dims
444 ,ITS,ITE, JTS,JTE, KTS,KTE & ! tile dims )
445 )
447 ! QV - water vapor mixing ratio (kg/kg)
448 ! QC - cloud water mixing ratio (kg/kg)
449 ! QR - rain water mixing ratio (kg/kg)
450 ! QI - cloud ice mixing ratio (kg/kg)
451 ! QS - snow mixing ratio (kg/kg)
452 ! QG - graupel mixing ratio (KG/KG)
453 ! NI - cloud ice number concentration (1/kg)
454 ! NS - Snow Number concentration (1/kg)
455 ! NR - Rain Number concentration (1/kg)
456 ! NG - Graupel number concentration (1/kg)
457 ! NOTE: RHO AND HT NOT USED BY THIS SCHEME AND DO NOT NEED TO BE PASSED INTO SCHEME!!!!
458 ! P - AIR PRESSURE (PA)
459 ! W - VERTICAL AIR VELOCITY (M/S)
460 ! TH - POTENTIAL TEMPERATURE (K)
461 ! PII - exner function - used to convert potential temp to temp
462 ! DZ - difference in height over interface (m)
463 ! DT_IN - model time step (sec)
464 ! ITIMESTEP - time step counter
465 ! RAINNC - accumulated grid-scale precipitation (mm)
466 ! RAINNCV - one time step grid scale precipitation (mm/time step)
467 ! SR - one time step mass ratio of snow to total precip
468 ! qrcuten, rain tendency from parameterized cumulus convection
469 ! qscuten, snow tendency from parameterized cumulus convection
470 ! qicuten, cloud ice tendency from parameterized cumulus convection
472 ! variables below currently not in use, not coupled to PBL or radiation codes
473 ! TKE - turbulence kinetic energy (m^2 s-2), NEEDED FOR DROPLET ACTIVATION (SEE CODE BELOW)
474 ! NCTEND - droplet concentration tendency from pbl (kg-1 s-1)
475 ! NCTEND - CLOUD ICE concentration tendency from pbl (kg-1 s-1)
476 ! KZH - heat eddy diffusion coefficient from YSU scheme (M^2 S-1), NEEDED FOR DROPLET ACTIVATION (SEE CODE BELOW)
477 ! EFFCS - CLOUD DROPLET EFFECTIVE RADIUS OUTPUT TO RADIATION CODE (micron)
478 ! EFFIS - CLOUD DROPLET EFFECTIVE RADIUS OUTPUT TO RADIATION CODE (micron)
479 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
481 ! reflectivity currently not included!!!!
482 ! REFL_10CM - CALCULATED RADAR REFLECTIVITY AT 10 CM (DBZ)
483 !................................
484 ! GRID_CLOCK, GRID_ALARMS - parameters to limit radar reflectivity calculation only when needed
485 ! otherwise radar reflectivity calculation every time step is too slow
486 ! only needed for coupling with WRF, see code below for details
487 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
489 ! EFFC - DROPLET EFFECTIVE RADIUS (MICRON)
490 ! EFFR - RAIN EFFECTIVE RADIUS (MICRON)
491 ! EFFS - SNOW EFFECTIVE RADIUS (MICRON)
492 ! EFFI - CLOUD ICE EFFECTIVE RADIUS (MICRON)
494 ! ADDITIONAL OUTPUT FROM MICRO - SEDIMENTATION TENDENCIES, NEEDED FOR LIQUID-ICE STATIC ENERGY
496 ! QGSTEN - GRAUPEL SEDIMENTATION TEND (KG/KG/S)
497 ! QRSTEN - RAIN SEDIMENTATION TEND (KG/KG/S)
498 ! QISTEN - CLOUD ICE SEDIMENTATION TEND (KG/KG/S)
499 ! QNISTEN - SNOW SEDIMENTATION TEND (KG/KG/S)
500 ! QCSTEN - CLOUD WATER SEDIMENTATION TEND (KG/KG/S)
502 ! ADDITIONAL INPUT NEEDED BY MICRO
503 ! ********NOTE: WVAR IS SHOULD BE USED IN DROPLET ACTIVATION
504 ! FOR CASES WHEN UPDRAFT IS NOT RESOLVED, EITHER BECAUSE OF
505 ! LOW MODEL RESOLUTION OR CLOUD TYPE
506 ! WVAR - STANDARD DEVIATION OF SUB-GRID VERTICAL VELOCITY (M/S)
508 IMPLICIT NONE
510 INTEGER, INTENT(IN ) :: ids, ide, jds, jde, kds, kde , &
511 ims, ime, jms, jme, kms, kme , &
512 its, ite, jts, jte, kts, kte
513 ! Temporary changed from INOUT to IN
515 REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT):: &
516 qv, qc, qr, qi, qs, qg, ni, ns, nr, TH, NG
517 !, effcs, effis
519 REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(IN):: &
520 pii, p, dz, rho, w !, tke, nctend, nitend,kzh
521 REAL, INTENT(IN):: dt_in
522 INTEGER, INTENT(IN):: ITIMESTEP
524 REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT):: &
525 RAINNC, RAINNCV, SR
527 ! REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT):: & ! GT
528 ! refl_10cm
530 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN) :: ht
532 ! TYPE (WRFU_Clock):: grid_clock ! GT
533 ! TYPE (WRFU_Alarm), POINTER:: grid_alarms(:) ! GT
535 ! LOCAL VARIABLES
537 REAL, DIMENSION(its:ite, kts:kte, jts:jte):: &
538 effi, effs, effr, EFFG
540 REAL, DIMENSION(its:ite, kts:kte, jts:jte):: &
541 T, WVAR, EFFC
543 REAL, DIMENSION(kts:kte) :: &
544 QC_TEND1D, QI_TEND1D, QNI_TEND1D, QR_TEND1D, &
545 NI_TEND1D, NS_TEND1D, NR_TEND1D, &
546 QC1D, QI1D, QR1D,NI1D, NS1D, NR1D, QS1D, &
547 T_TEND1D,QV_TEND1D, T1D, QV1D, P1D, W1D, WVAR1D, &
548 EFFC1D, EFFI1D, EFFS1D, EFFR1D,DZ1D, &
549 ! HM ADD GRAUPEL
550 QG_TEND1D, NG_TEND1D, QG1D, NG1D, EFFG1D, &
552 ! ADD SEDIMENTATION TENDENCIES (UNITS OF KG/KG/S)
553 QGSTEN,QRSTEN, QISTEN, QNISTEN, QCSTEN, &
554 ! ADD CUMULUS TENDENCIES
555 QRCU1D, QSCU1D, QICU1D
557 ! add cumulus tendencies
559 REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(IN):: &
560 qrcuten, qscuten, qicuten
561 REAL, DIMENSION(ims:ime, jms:jme), INTENT(IN):: &
562 mu
564 ! HM add reflectivity
565 ! dbz
567 REAL PRECPRT1D, SNOWRT1D
569 INTEGER I,K,J
571 REAL DT
573 ! LOGICAL:: dBZ_tstep ! GT
575 ! Initialize tendencies (all set to 0) and transfer
576 ! array to local variables
577 DT = DT_IN
579 DO I=ITS,ITE
580 DO J=JTS,JTE
581 DO K=KTS,KTE
582 T(I,K,J) = TH(i,k,j)*PII(i,k,j)
584 ! wvar is the ST. DEV. OF sub-grid vertical velocity, used for calculating droplet
585 ! activation rates.
586 ! WVAR CAN BE DERIVED EITHER FROM PREDICTED TKE (AS IN MYJ PBL SCHEME),
587 ! OR FROM EDDY DIFFUSION COEFFICIENT KZH (AS IN YSU PBL SCHEME),
588 ! DEPENDING ON THE PARTICULAR pbl SCHEME DRIVER MODEL IS COUPLED WITH
589 ! NOTE: IF MODEL HAS HIGH ENOUGH RESOLUTION TO RESOLVE UPDRAFTS, WVAR MAY
590 ! NOT BE NEEDED
592 ! currently assign wvar to 0.5 m/s (not coupled with PBL scheme)
594 WVAR(I,K,J) = 0.5
596 ! currently mixing of number concentrations also is neglected (not coupled with PBL schemes)
598 END DO
599 END DO
600 END DO
602 do i=its,ite ! i loop (east-west)
603 do j=jts,jte ! j loop (north-south)
604 !
605 ! Transfer 3D arrays into 1D for microphysical calculations
606 !
608 ! hm , initialize 1d tendency arrays to zero
610 do k=kts,kte ! k loop (vertical)
612 QC_TEND1D(k) = 0.
613 QI_TEND1D(k) = 0.
614 QNI_TEND1D(k) = 0.
615 QR_TEND1D(k) = 0.
616 NI_TEND1D(k) = 0.
617 NS_TEND1D(k) = 0.
618 NR_TEND1D(k) = 0.
619 T_TEND1D(k) = 0.
620 QV_TEND1D(k) = 0.
622 QC1D(k) = QC(i,k,j)
623 QI1D(k) = QI(i,k,j)
624 QS1D(k) = QS(i,k,j)
625 QR1D(k) = QR(i,k,j)
627 NI1D(k) = NI(i,k,j)
629 NS1D(k) = NS(i,k,j)
630 NR1D(k) = NR(i,k,j)
631 ! HM ADD GRAUPEL
632 QG1D(K) = QG(I,K,j)
633 NG1D(K) = NG(I,K,j)
634 QG_TEND1D(K) = 0.
635 NG_TEND1D(K) = 0.
637 T1D(k) = T(i,k,j)
638 QV1D(k) = QV(i,k,j)
639 P1D(k) = P(i,k,j)
640 DZ1D(k) = DZ(i,k,j)
641 W1D(k) = W(i,k,j)
642 WVAR1D(k) = WVAR(i,k,j)
643 ! add cumulus tendencies, decouple from mu
644 qrcu1d(k) = qrcuten(i,k,j)/mu(i,j)
645 qscu1d(k) = qscuten(i,k,j)/mu(i,j)
646 qicu1d(k) = qicuten(i,k,j)/mu(i,j)
647 end do
649 call MORR_TWO_MOMENT_MICRO(QC_TEND1D, QI_TEND1D, QNI_TEND1D, QR_TEND1D, &
650 NI_TEND1D, NS_TEND1D, NR_TEND1D, &
651 QC1D, QI1D, QS1D, QR1D,NI1D, NS1D, NR1D, &
652 T_TEND1D,QV_TEND1D, T1D, QV1D, P1D, DZ1D, W1D, WVAR1D, &
653 PRECPRT1D,SNOWRT1D, &
654 EFFC1D,EFFI1D,EFFS1D,EFFR1D,DT, &
655 IMS,IME, JMS,JME, KMS,KME, &
656 ITS,ITE, JTS,JTE, KTS,KTE, & ! HM ADD GRAUPEL
657 QG_TEND1D,NG_TEND1D,QG1D,NG1D,EFFG1D, &
658 qrcu1d, qscu1d, qicu1d, &
659 ! ADD SEDIMENTATION TENDENCIES
660 QGSTEN,QRSTEN,QISTEN,QNISTEN,QCSTEN)
662 !
663 ! Transfer 1D arrays back into 3D arrays
664 !
665 do k=kts,kte
667 ! hm, add tendencies to update global variables
668 ! HM, TENDENCIES FOR Q AND N NOW ADDED IN M2005MICRO, SO WE
669 ! ONLY NEED TO TRANSFER 1D VARIABLES BACK TO 3D
671 QC(i,k,j) = QC1D(k)
672 QI(i,k,j) = QI1D(k)
673 QS(i,k,j) = QS1D(k)
674 QR(i,k,j) = QR1D(k)
675 NI(i,k,j) = NI1D(k)
676 NS(i,k,j) = NS1D(k)
677 NR(i,k,j) = NR1D(k)
678 QG(I,K,j) = QG1D(K)
679 NG(I,K,j) = NG1D(K)
681 T(i,k,j) = T1D(k)
682 TH(I,K,J) = T(i,k,j)/PII(i,k,j) ! CONVERT TEMP BACK TO POTENTIAL TEMP
683 QV(i,k,j) = QV1D(k)
685 EFFC(i,k,j) = EFFC1D(k)
686 EFFI(i,k,j) = EFFI1D(k)
687 EFFS(i,k,j) = EFFS1D(k)
688 EFFR(i,k,j) = EFFR1D(k)
689 EFFG(I,K,j) = EFFG1D(K)
691 ! EFFECTIVE RADIUS FOR RADIATION CODE (currently not coupled)
692 ! HM, ADD LIMIT TO PREVENT BLOWING UP OPTICAL PROPERTIES, 8/18/07
693 ! EFFCS(I,K,J) = MIN(EFFC(I,K,J),50.)
694 ! EFFCS(I,K,J) = MAX(EFFCS(I,K,J),1.)
695 ! EFFIS(I,K,J) = MIN(EFFI(I,K,J),130.)
696 ! EFFIS(I,K,J) = MAX(EFFIS(I,K,J),13.)
698 end do
700 ! hm modified so that m2005 precip variables correctly match wrf precip variables
701 RAINNC(i,j) = RAINNC(I,J)+PRECPRT1D
702 RAINNCV(i,j) = PRECPRT1D
703 SR(i,j) = SNOWRT1D/(PRECPRT1D+1.E-12)
705 end do
706 end do
708 END SUBROUTINE MP_MORR_TWO_MOMENT
710 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
711 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
712 SUBROUTINE MORR_TWO_MOMENT_MICRO(QC3DTEN,QI3DTEN,QNI3DTEN,QR3DTEN, &
713 NI3DTEN,NS3DTEN,NR3DTEN,QC3D,QI3D,QNI3D,QR3D,NI3D,NS3D,NR3D, &
714 T3DTEN,QV3DTEN,T3D,QV3D,PRES,DZQ,W3D,WVAR,PRECRT,SNOWRT, &
715 EFFC,EFFI,EFFS,EFFR,DT, &
716 IMS,IME, JMS,JME, KMS,KME, &
717 ITS,ITE, JTS,JTE, KTS,KTE, & ! ADD GRAUPEL
718 QG3DTEN,NG3DTEN,QG3D,NG3D,EFFG,qrcu1d,qscu1d, qicu1d, &
719 QGSTEN,QRSTEN,QISTEN,QNISTEN,QCSTEN)
721 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
722 ! THIS PROGRAM IS THE MAIN TWO-MOMENT MICROPHYSICS SUBROUTINE DESCRIBED BY
723 ! MORRISON ET AL. 2005 JAS; MORRISON AND PINTO 2005 JAS.
724 ! ADDITIONAL CHANGES ARE DESCRIBED IN DETAIL BY MORRISON, THOMPSON, TATARSKII (MWR, SUBMITTED)
726 ! THIS SCHEME IS A BULK DOUBLE-MOMENT SCHEME THAT PREDICTS MIXING
727 ! RATIOS AND NUMBER CONCENTRATIONS OF FIVE HYDROMETEOR SPECIES:
728 ! CLOUD DROPLETS, CLOUD (SMALL) ICE, RAIN, SNOW, AND GRAUPEL.
730 ! CODE STRUCTURE: MAIN SUBROUTINE IS 'MORR_TWO_MOMENT'. ALSO INCLUDED IN THIS FILE IS
731 ! 'FUNCTION POLYSVP', 'FUNCTION DERF1', AND
732 ! 'FUNCTION GAMMA'.
734 ! NOTE: THIS SUBROUTINE USES 1D ARRAY IN VERTICAL (COLUMN), EVEN THOUGH VARIABLES ARE CALLED '3D'......
736 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
738 ! DECLARATIONS
740 IMPLICIT NONE
742 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
743 ! THESE VARIABLES BELOW MUST BE LINKED WITH THE MAIN MODEL.
744 ! DEFINE ARRAY SIZES
746 ! INPUT NUMBER OF GRID CELLS
748 ! INPUT/OUTPUT PARAMETERS ! DESCRIPTION (UNITS)
749 INTEGER, INTENT( IN) :: IMS,IME, JMS,JME, KMS,KME, &
750 ITS,ITE, JTS,JTE, KTS,KTE
752 REAL, DIMENSION(KTS:KTE) :: QC3DTEN ! CLOUD WATER MIXING RATIO TENDENCY (KG/KG/S)
753 REAL, DIMENSION(KTS:KTE) :: QI3DTEN ! CLOUD ICE MIXING RATIO TENDENCY (KG/KG/S)
754 REAL, DIMENSION(KTS:KTE) :: QNI3DTEN ! SNOW MIXING RATIO TENDENCY (KG/KG/S)
755 REAL, DIMENSION(KTS:KTE) :: QR3DTEN ! RAIN MIXING RATIO TENDENCY (KG/KG/S)
756 REAL, DIMENSION(KTS:KTE) :: NI3DTEN ! CLOUD ICE NUMBER CONCENTRATION (1/KG/S)
757 REAL, DIMENSION(KTS:KTE) :: NS3DTEN ! SNOW NUMBER CONCENTRATION (1/KG/S)
758 REAL, DIMENSION(KTS:KTE) :: NR3DTEN ! RAIN NUMBER CONCENTRATION (1/KG/S)
759 REAL, DIMENSION(KTS:KTE) :: QC3D ! CLOUD WATER MIXING RATIO (KG/KG)
760 REAL, DIMENSION(KTS:KTE) :: QI3D ! CLOUD ICE MIXING RATIO (KG/KG)
761 REAL, DIMENSION(KTS:KTE) :: QNI3D ! SNOW MIXING RATIO (KG/KG)
762 REAL, DIMENSION(KTS:KTE) :: QR3D ! RAIN MIXING RATIO (KG/KG)
763 REAL, DIMENSION(KTS:KTE) :: NI3D ! CLOUD ICE NUMBER CONCENTRATION (1/KG)
764 REAL, DIMENSION(KTS:KTE) :: NS3D ! SNOW NUMBER CONCENTRATION (1/KG)
765 REAL, DIMENSION(KTS:KTE) :: NR3D ! RAIN NUMBER CONCENTRATION (1/KG)
766 REAL, DIMENSION(KTS:KTE) :: T3DTEN ! TEMPERATURE TENDENCY (K/S)
767 REAL, DIMENSION(KTS:KTE) :: QV3DTEN ! WATER VAPOR MIXING RATIO TENDENCY (KG/KG/S)
768 REAL, DIMENSION(KTS:KTE) :: T3D ! TEMPERATURE (K)
769 REAL, DIMENSION(KTS:KTE) :: QV3D ! WATER VAPOR MIXING RATIO (KG/KG)
770 REAL, DIMENSION(KTS:KTE) :: PRES ! ATMOSPHERIC PRESSURE (PA)
771 REAL, DIMENSION(KTS:KTE) :: DZQ ! DIFFERENCE IN HEIGHT ACROSS LEVEL (m)
772 REAL, DIMENSION(KTS:KTE) :: W3D ! GRID-SCALE VERTICAL VELOCITY (M/S)
773 REAL, DIMENSION(KTS:KTE) :: WVAR ! SUB-GRID VERTICAL VELOCITY (M/S)
775 ! HM ADDED GRAUPEL VARIABLES
776 REAL, DIMENSION(KTS:KTE) :: QG3DTEN ! GRAUPEL MIX RATIO TENDENCY (KG/KG/S)
777 REAL, DIMENSION(KTS:KTE) :: NG3DTEN ! GRAUPEL NUMB CONC TENDENCY (1/KG/S)
778 REAL, DIMENSION(KTS:KTE) :: QG3D ! GRAUPEL MIX RATIO (KG/KG)
779 REAL, DIMENSION(KTS:KTE) :: NG3D ! GRAUPEL NUMBER CONC (1/KG)
781 ! HM, ADD 1/16/07, SEDIMENTATION TENDENCIES FOR MIXING RATIO
783 REAL, DIMENSION(KTS:KTE) :: QGSTEN ! GRAUPEL SED TEND (KG/KG/S)
784 REAL, DIMENSION(KTS:KTE) :: QRSTEN ! RAIN SED TEND (KG/KG/S)
785 REAL, DIMENSION(KTS:KTE) :: QISTEN ! CLOUD ICE SED TEND (KG/KG/S)
786 REAL, DIMENSION(KTS:KTE) :: QNISTEN ! SNOW SED TEND (KG/KG/S)
787 REAL, DIMENSION(KTS:KTE) :: QCSTEN ! CLOUD WAT SED TEND (KG/KG/S)
789 ! hm add cumulus tendencies for precip
790 REAL, DIMENSION(KTS:KTE) :: qrcu1d
791 REAL, DIMENSION(KTS:KTE) :: qscu1d
792 REAL, DIMENSION(KTS:KTE) :: qicu1d
794 ! OUTPUT VARIABLES
796 REAL PRECRT ! TOTAL PRECIP PER TIME STEP (mm)
797 REAL SNOWRT ! SNOW PER TIME STEP (mm)
799 REAL, DIMENSION(KTS:KTE) :: EFFC ! DROPLET EFFECTIVE RADIUS (MICRON)
800 REAL, DIMENSION(KTS:KTE) :: EFFI ! CLOUD ICE EFFECTIVE RADIUS (MICRON)
801 REAL, DIMENSION(KTS:KTE) :: EFFS ! SNOW EFFECTIVE RADIUS (MICRON)
802 REAL, DIMENSION(KTS:KTE) :: EFFR ! RAIN EFFECTIVE RADIUS (MICRON)
803 REAL, DIMENSION(KTS:KTE) :: EFFG ! GRAUPEL EFFECTIVE RADIUS (MICRON)
805 ! MODEL INPUT PARAMETERS (FORMERLY IN COMMON BLOCKS)
807 REAL DT ! MODEL TIME STEP (SEC)
809 !.....................................................................................................
810 ! LOCAL VARIABLES: ALL PARAMETERS BELOW ARE LOCAL TO SCHEME AND DON'T NEED TO COMMUNICATE WITH THE
811 ! REST OF THE MODEL.
813 ! SIZE PARAMETER VARIABLES
815 REAL, DIMENSION(KTS:KTE) :: LAMC ! SLOPE PARAMETER FOR DROPLETS (M-1)
816 REAL, DIMENSION(KTS:KTE) :: LAMI ! SLOPE PARAMETER FOR CLOUD ICE (M-1)
817 REAL, DIMENSION(KTS:KTE) :: LAMS ! SLOPE PARAMETER FOR SNOW (M-1)
818 REAL, DIMENSION(KTS:KTE) :: LAMR ! SLOPE PARAMETER FOR RAIN (M-1)
819 REAL, DIMENSION(KTS:KTE) :: LAMG ! SLOPE PARAMETER FOR GRAUPEL (M-1)
820 REAL, DIMENSION(KTS:KTE) :: CDIST1 ! PSD PARAMETER FOR DROPLETS
821 REAL, DIMENSION(KTS:KTE) :: N0I ! INTERCEPT PARAMETER FOR CLOUD ICE (KG-1 M-1)
822 REAL, DIMENSION(KTS:KTE) :: N0S ! INTERCEPT PARAMETER FOR SNOW (KG-1 M-1)
823 REAL, DIMENSION(KTS:KTE) :: N0RR ! INTERCEPT PARAMETER FOR RAIN (KG-1 M-1)
824 REAL, DIMENSION(KTS:KTE) :: N0G ! INTERCEPT PARAMETER FOR GRAUPEL (KG-1 M-1)
825 REAL, DIMENSION(KTS:KTE) :: PGAM ! SPECTRAL SHAPE PARAMETER FOR DROPLETS
827 ! MICROPHYSICAL PROCESSES
829 REAL, DIMENSION(KTS:KTE) :: NSUBC ! LOSS OF NC DURING EVAP
830 REAL, DIMENSION(KTS:KTE) :: NSUBI ! LOSS OF NI DURING SUB.
831 REAL, DIMENSION(KTS:KTE) :: NSUBS ! LOSS OF NS DURING SUB.
832 REAL, DIMENSION(KTS:KTE) :: NSUBR ! LOSS OF NR DURING EVAP
833 REAL, DIMENSION(KTS:KTE) :: PRD ! DEP CLOUD ICE
834 REAL, DIMENSION(KTS:KTE) :: PRE ! EVAP OF RAIN
835 REAL, DIMENSION(KTS:KTE) :: PRDS ! DEP SNOW
836 REAL, DIMENSION(KTS:KTE) :: NNUCCC ! CHANGE N DUE TO CONTACT FREEZ DROPLETS
837 REAL, DIMENSION(KTS:KTE) :: MNUCCC ! CHANGE Q DUE TO CONTACT FREEZ DROPLETS
838 REAL, DIMENSION(KTS:KTE) :: PRA ! ACCRETION DROPLETS BY RAIN
839 REAL, DIMENSION(KTS:KTE) :: PRC ! AUTOCONVERSION DROPLETS
840 REAL, DIMENSION(KTS:KTE) :: PCC ! COND/EVAP DROPLETS
841 REAL, DIMENSION(KTS:KTE) :: NNUCCD ! CHANGE N FREEZING AEROSOL (PRIM ICE NUCLEATION)
842 REAL, DIMENSION(KTS:KTE) :: MNUCCD ! CHANGE Q FREEZING AEROSOL (PRIM ICE NUCLEATION)
843 REAL, DIMENSION(KTS:KTE) :: MNUCCR ! CHANGE Q DUE TO CONTACT FREEZ RAIN
844 REAL, DIMENSION(KTS:KTE) :: NNUCCR ! CHANGE N DUE TO CONTACT FREEZ RAIN
845 REAL, DIMENSION(KTS:KTE) :: NPRA ! CHANGE IN N DUE TO DROPLET ACC BY RAIN
846 REAL, DIMENSION(KTS:KTE) :: NRAGG ! SELF-COLLECTION OF RAIN
847 REAL, DIMENSION(KTS:KTE) :: NSAGG ! SELF-COLLECTION OF SNOW
848 REAL, DIMENSION(KTS:KTE) :: NPRC ! CHANGE NC AUTOCONVERSION DROPLETS
849 REAL, DIMENSION(KTS:KTE) :: NPRC1 ! CHANGE NR AUTOCONVERSION DROPLETS
850 REAL, DIMENSION(KTS:KTE) :: PRAI ! CHANGE Q AUTOCONVERSION CLOUD ICE
851 REAL, DIMENSION(KTS:KTE) :: PRCI ! CHANGE Q ACCRETION CLOUD ICE BY SNOW
852 REAL, DIMENSION(KTS:KTE) :: PSACWS ! CHANGE Q DROPLET ACCRETION BY SNOW
853 REAL, DIMENSION(KTS:KTE) :: NPSACWS ! CHANGE N DROPLET ACCRETION BY SNOW
854 REAL, DIMENSION(KTS:KTE) :: PSACWI ! CHANGE Q DROPLET ACCRETION BY CLOUD ICE
855 REAL, DIMENSION(KTS:KTE) :: NPSACWI ! CHANGE N DROPLET ACCRETION BY CLOUD ICE
856 REAL, DIMENSION(KTS:KTE) :: NPRCI ! CHANGE N AUTOCONVERSION CLOUD ICE BY SNOW
857 REAL, DIMENSION(KTS:KTE) :: NPRAI ! CHANGE N ACCRETION CLOUD ICE
858 REAL, DIMENSION(KTS:KTE) :: NMULTS ! ICE MULT DUE TO RIMING DROPLETS BY SNOW
859 REAL, DIMENSION(KTS:KTE) :: NMULTR ! ICE MULT DUE TO RIMING RAIN BY SNOW
860 REAL, DIMENSION(KTS:KTE) :: QMULTS ! CHANGE Q DUE TO ICE MULT DROPLETS/SNOW
861 REAL, DIMENSION(KTS:KTE) :: QMULTR ! CHANGE Q DUE TO ICE RAIN/SNOW
862 REAL, DIMENSION(KTS:KTE) :: PRACS ! CHANGE Q RAIN-SNOW COLLECTION
863 REAL, DIMENSION(KTS:KTE) :: NPRACS ! CHANGE N RAIN-SNOW COLLECTION
864 REAL, DIMENSION(KTS:KTE) :: PCCN ! CHANGE Q DROPLET ACTIVATION
865 REAL, DIMENSION(KTS:KTE) :: PSMLT ! CHANGE Q MELTING SNOW TO RAIN
866 REAL, DIMENSION(KTS:KTE) :: EVPMS ! CHNAGE Q MELTING SNOW EVAPORATING
867 REAL, DIMENSION(KTS:KTE) :: NSMLTS ! CHANGE N MELTING SNOW
868 REAL, DIMENSION(KTS:KTE) :: NSMLTR ! CHANGE N MELTING SNOW TO RAIN
869 ! HM ADDED 12/13/06
870 REAL, DIMENSION(KTS:KTE) :: PIACR ! CHANGE QR, ICE-RAIN COLLECTION
871 REAL, DIMENSION(KTS:KTE) :: NIACR ! CHANGE N, ICE-RAIN COLLECTION
872 REAL, DIMENSION(KTS:KTE) :: PRACI ! CHANGE QI, ICE-RAIN COLLECTION
873 REAL, DIMENSION(KTS:KTE) :: PIACRS ! CHANGE QR, ICE RAIN COLLISION, ADDED TO SNOW
874 REAL, DIMENSION(KTS:KTE) :: NIACRS ! CHANGE N, ICE RAIN COLLISION, ADDED TO SNOW
875 REAL, DIMENSION(KTS:KTE) :: PRACIS ! CHANGE QI, ICE RAIN COLLISION, ADDED TO SNOW
876 REAL, DIMENSION(KTS:KTE) :: EPRD ! SUBLIMATION CLOUD ICE
877 REAL, DIMENSION(KTS:KTE) :: EPRDS ! SUBLIMATION SNOW
878 ! HM ADDED GRAUPEL PROCESSES
879 REAL, DIMENSION(KTS:KTE) :: PRACG ! CHANGE IN Q COLLECTION RAIN BY GRAUPEL
880 REAL, DIMENSION(KTS:KTE) :: PSACWG ! CHANGE IN Q COLLECTION DROPLETS BY GRAUPEL
881 REAL, DIMENSION(KTS:KTE) :: PGSACW ! CONVERSION Q TO GRAUPEL DUE TO COLLECTION DROPLETS BY SNOW
882 REAL, DIMENSION(KTS:KTE) :: PGRACS ! CONVERSION Q TO GRAUPEL DUE TO COLLECTION RAIN BY SNOW
883 REAL, DIMENSION(KTS:KTE) :: PRDG ! DEP OF GRAUPEL
884 REAL, DIMENSION(KTS:KTE) :: EPRDG ! SUB OF GRAUPEL
885 REAL, DIMENSION(KTS:KTE) :: EVPMG ! CHANGE Q MELTING OF GRAUPEL AND EVAPORATION
886 REAL, DIMENSION(KTS:KTE) :: PGMLT ! CHANGE Q MELTING OF GRAUPEL
887 REAL, DIMENSION(KTS:KTE) :: NPRACG ! CHANGE N COLLECTION RAIN BY GRAUPEL
888 REAL, DIMENSION(KTS:KTE) :: NPSACWG ! CHANGE N COLLECTION DROPLETS BY GRAUPEL
889 REAL, DIMENSION(KTS:KTE) :: NSCNG ! CHANGE N CONVERSION TO GRAUPEL DUE TO COLLECTION DROPLETS BY SNOW
890 REAL, DIMENSION(KTS:KTE) :: NGRACS ! CHANGE N CONVERSION TO GRAUPEL DUE TO COLLECTION RAIN BY SNOW
891 REAL, DIMENSION(KTS:KTE) :: NGMLTG ! CHANGE N MELTING GRAUPEL
892 REAL, DIMENSION(KTS:KTE) :: NGMLTR ! CHANGE N MELTING GRAUPEL TO RAIN
893 REAL, DIMENSION(KTS:KTE) :: NSUBG ! CHANGE N SUB/DEP OF GRAUPEL
894 REAL, DIMENSION(KTS:KTE) :: PSACR ! CONVERSION DUE TO COLL OF SNOW BY RAIN
895 REAL, DIMENSION(KTS:KTE) :: NMULTG ! ICE MULT DUE TO ACC DROPLETS BY GRAUPEL
896 REAL, DIMENSION(KTS:KTE) :: NMULTRG ! ICE MULT DUE TO ACC RAIN BY GRAUPEL
897 REAL, DIMENSION(KTS:KTE) :: QMULTG ! CHANGE Q DUE TO ICE MULT DROPLETS/GRAUPEL
898 REAL, DIMENSION(KTS:KTE) :: QMULTRG ! CHANGE Q DUE TO ICE MULT RAIN/GRAUPEL
900 ! TIME-VARYING ATMOSPHERIC PARAMETERS
902 REAL, DIMENSION(KTS:KTE) :: KAP ! THERMAL CONDUCTIVITY OF AIR
903 REAL, DIMENSION(KTS:KTE) :: EVS ! SATURATION VAPOR PRESSURE
904 REAL, DIMENSION(KTS:KTE) :: EIS ! ICE SATURATION VAPOR PRESSURE
905 REAL, DIMENSION(KTS:KTE) :: QVS ! SATURATION MIXING RATIO
906 REAL, DIMENSION(KTS:KTE) :: QVI ! ICE SATURATION MIXING RATIO
907 REAL, DIMENSION(KTS:KTE) :: QVQVS ! SAUTRATION RATIO
908 REAL, DIMENSION(KTS:KTE) :: QVQVSI! ICE SATURAION RATIO
909 REAL, DIMENSION(KTS:KTE) :: DV ! DIFFUSIVITY OF WATER VAPOR IN AIR
910 REAL, DIMENSION(KTS:KTE) :: XXLS ! LATENT HEAT OF SUBLIMATION
911 REAL, DIMENSION(KTS:KTE) :: XXLV ! LATENT HEAT OF VAPORIZATION
912 REAL, DIMENSION(KTS:KTE) :: CPM ! SPECIFIC HEAT AT CONST PRESSURE FOR MOIST AIR
913 REAL, DIMENSION(KTS:KTE) :: MU ! VISCOCITY OF AIR
914 REAL, DIMENSION(KTS:KTE) :: SC ! SCHMIDT NUMBER
915 REAL, DIMENSION(KTS:KTE) :: XLF ! LATENT HEAT OF FREEZING
916 REAL, DIMENSION(KTS:KTE) :: RHO ! AIR DENSITY
917 REAL, DIMENSION(KTS:KTE) :: AB ! CORRECTION TO CONDENSATION RATE DUE TO LATENT HEATING
918 REAL, DIMENSION(KTS:KTE) :: ABI ! CORRECTION TO DEPOSITION RATE DUE TO LATENT HEATING
920 ! TIME-VARYING MICROPHYSICS PARAMETERS
922 REAL, DIMENSION(KTS:KTE) :: DAP ! DIFFUSIVITY OF AEROSOL
923 REAL NACNT ! NUMBER OF CONTACT IN
924 REAL FMULT ! TEMP.-DEP. PARAMETER FOR RIME-SPLINTERING
925 REAL COFFI ! ICE AUTOCONVERSION PARAMETER
927 ! FALL SPEED WORKING VARIABLES (DEFINED IN CODE)
929 REAL, DIMENSION(KTS:KTE) :: DUMI,DUMR,DUMFNI,DUMG,DUMFNG
930 REAL UNI, UMI,UMR
931 REAL, DIMENSION(KTS:KTE) :: FR, FI, FNI,FG,FNG
932 REAL RGVM
933 REAL, DIMENSION(KTS:KTE) :: FALOUTR,FALOUTI,FALOUTNI
934 REAL FALTNDR,FALTNDI,FALTNDNI,RHO2
935 REAL, DIMENSION(KTS:KTE) :: DUMQS,DUMFNS
936 REAL UMS,UNS
937 REAL, DIMENSION(KTS:KTE) :: FS,FNS, FALOUTS,FALOUTNS,FALOUTG,FALOUTNG
938 REAL FALTNDS,FALTNDNS,UNR,FALTNDG,FALTNDNG
939 REAL, DIMENSION(KTS:KTE) :: DUMC,DUMFNC
940 REAL UNC,UMC,UNG,UMG
941 REAL, DIMENSION(KTS:KTE) :: FC,FALOUTC,FALOUTNC
942 REAL FALTNDC,FALTNDNC
943 REAL, DIMENSION(KTS:KTE) :: FNC,DUMFNR,FALOUTNR
944 REAL FALTNDNR
945 REAL, DIMENSION(KTS:KTE) :: FNR(KMS:KME)
947 ! FALL-SPEED PARAMETER 'A' WITH AIR DENSITY CORRECTION
949 REAL, DIMENSION(KTS:KTE) :: AIN,ARN,ASN,ACN,AGN
951 ! EXTERNAL FUNCTION CALL RETURN VARIABLES
953 ! REAL GAMMA, ! EULER GAMMA FUNCTION
954 ! REAL POLYSVP, ! SAT. PRESSURE FUNCTION
955 ! REAL DERF1 ! ERROR FUNCTION
957 ! DUMMY VARIABLES
959 REAL DUM,DUM1,DUM2,DUMT,DUMQV,DUMQSS,DUMQSI,DUMS
961 ! PROGNOSTIC SUPERSATURATION
963 REAL DQSDT ! CHANGE OF SAT. MIX. RAT. WITH TEMPERATURE
964 REAL DQSIDT ! CHANGE IN ICE SAT. MIXING RAT. WITH T
965 REAL EPSI ! 1/PHASE REL. TIME (SEE M2005), ICE
966 REAL EPSS ! 1/PHASE REL. TIME (SEE M2005), SNOW
967 REAL EPSR ! 1/PHASE REL. TIME (SEE M2005), RAIN
968 REAL EPSG ! 1/PHASE REL. TIME (SEE M2005), GRAUPEL
970 ! NEW DROPLET ACTIVATION VARIABLES
971 REAL TAUC ! PHASE REL. TIME (SEE M2005), DROPLETS
972 REAL TAUR ! PHASE REL. TIME (SEE M2005), RAIN
973 REAL TAUI ! PHASE REL. TIME (SEE M2005), CLOUD ICE
974 REAL TAUS ! PHASE REL. TIME (SEE M2005), SNOW
975 REAL TAUG ! PHASE REL. TIME (SEE M2005), GRAUPEL
976 REAL DUMACT,DUM3
978 ! COUNTING/INDEX VARIABLES
980 INTEGER K,NSTEP,N ! ,I
982 ! LTRUE IS ONLY USED TO SPEED UP THE CODE !!
983 ! LTRUE, SWITCH = 0, NO HYDROMETEORS IN COLUMN,
984 ! = 1, HYDROMETEORS IN COLUMN
986 INTEGER LTRUE
988 ! DROPLET ACTIVATION/FREEZING AEROSOL
991 REAL CT ! DROPLET ACTIVATION PARAMETER
992 REAL TEMP1 ! DUMMY TEMPERATURE
993 REAL SAT1 ! DUMMY SATURATION
994 REAL SIGVL ! SURFACE TENSION LIQ/VAPOR
995 REAL KEL ! KELVIN PARAMETER
996 REAL KC2 ! TOTAL ICE NUCLEATION RATE
998 REAL CRY,KRY ! AEROSOL ACTIVATION PARAMETERS
1000 ! MORE WORKING/DUMMY VARIABLES
1002 REAL DUMQI,DUMNI,DC0,DS0,DG0
1003 REAL DUMQC,DUMQR,RATIO,SUM_DEP,FUDGEF
1005 ! EFFECTIVE VERTICAL VELOCITY (M/S)
1006 REAL WEF
1008 ! WORKING PARAMETERS FOR ICE NUCLEATION
1010 REAL ANUC,BNUC
1012 ! WORKING PARAMETERS FOR AEROSOL ACTIVATION
1014 REAL AACT,GAMM,GG,PSI,ETA1,ETA2,SM1,SM2,SMAX,UU1,UU2,ALPHA
1016 ! DUMMY SIZE DISTRIBUTION PARAMETERS
1018 REAL DLAMS,DLAMR,DLAMI,DLAMC,DLAMG,LAMMAX,LAMMIN
1020 INTEGER IDROP
1022 ! DROPLET CONCENTRATION AND ITS TENDENCY
1023 ! NOTE: CURRENTLY DROPLET CONCENTRATION IS SPECIFIED !!!!!
1024 ! TENDENCY OF NC IS CALCULATED BUT IT IS NOT USED !!!
1025 REAL, DIMENSION(KTS:KTE) :: NC3DTEN ! CLOUD DROPLET NUMBER CONCENTRATION (1/KG/S)
1026 REAL, DIMENSION(KTS:KTE) :: NC3D ! CLOUD DROPLET NUMBER CONCENTRATION (1/KG)
1028 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1030 ! SET LTRUE INITIALLY TO 0
1032 LTRUE = 0
1034 ! ATMOSPHERIC PARAMETERS THAT VARY IN TIME AND HEIGHT
1035 DO K = KTS,KTE
1037 ! LATENT HEAT OF VAPORATION
1039 XXLV(K) = 3.1484E6-2370.*T3D(K)
1041 ! LATENT HEAT OF SUBLIMATION
1043 XXLS(K) = 3.15E6-2370.*T3D(K)+0.3337E6
1045 CPM(K) = CP*(1.+0.887*QV3D(K))
1047 ! SATURATION VAPOR PRESSURE AND MIXING RATIO
1049 EVS(K) = POLYSVP(T3D(K),0) ! PA
1050 EIS(K) = POLYSVP(T3D(K),1) ! PA
1052 ! MAKE SURE ICE SATURATION DOESN'T EXCEED WATER SAT. NEAR FREEZING
1054 IF (EIS(K).GT.EVS(K)) EIS(K) = EVS(K)
1056 QVS(K) = .622*EVS(K)/(PRES(K)-EVS(K))
1057 QVI(K) = .622*EIS(K)/(PRES(K)-EIS(K))
1059 QVQVS(K) = QV3D(K)/QVS(K)
1060 QVQVSI(K) = QV3D(K)/QVI(K)
1062 ! AIR DENSITY
1064 RHO(K) = PRES(K)/(R*T3D(K))
1066 ! ADD NUMBER CONCENTRATION DUE TO CUMULUS TENDENCY
1067 ! ASSUME N0 ASSOCIATED WITH CUMULUS PARAM RAIN IS 10^7 M^-4
1068 ! ASSUME N0 ASSOCIATED WITH CUMULUS PARAM SNOW IS 2 X 10^7 M^-4
1069 ! FOR DETRAINED CLOUD ICE, ASSUME MEAN VOLUME DIAM OF 80 MICRON
1071 IF (QRCU1D(K).GE.1.E-10) THEN
1072 DUM=1.8e5*(QRCU1D(K)*DT/(PI*RHOW*RHO(K)**3))**0.25
1073 NR3D(K)=NR3D(K)+DUM
1074 END IF
1075 IF (QSCU1D(K).GE.1.E-10) THEN
1076 DUM=3.e5*(QSCU1D(K)*DT/(CONS1*RHO(K)**3))**(1./(DS+1.))
1077 NS3D(K)=NS3D(K)+DUM
1078 END IF
1079 IF (QICU1D(K).GE.1.E-10) THEN
1080 DUM=QICU1D(K)*DT/(CI*(80.E-6)**DI)
1081 NI3D(K)=NI3D(K)+DUM
1082 END IF
1084 ! AT SUBSATURATION, REMOVE SMALL AMOUNTS OF CLOUD/PRECIP WATER
1086 IF (QVQVS(K).LT.0.9) THEN
1087 IF (QR3D(K).LT.1.E-6) THEN
1088 QV3D(K)=QV3D(K)+QR3D(K)
1089 T3D(K)=T3D(K)-QR3D(K)*XXLV(K)/CPM(K)
1090 QR3D(K)=0.
1091 END IF
1092 IF (QC3D(K).LT.1.E-6) THEN
1093 QV3D(K)=QV3D(K)+QC3D(K)
1094 T3D(K)=T3D(K)-QC3D(K)*XXLV(K)/CPM(K)
1095 QC3D(K)=0.
1096 END IF
1097 END IF
1099 IF (QVQVSI(K).LT.0.9) THEN
1100 IF (QI3D(K).LT.1.E-6) THEN
1101 QV3D(K)=QV3D(K)+QI3D(K)
1102 T3D(K)=T3D(K)-QI3D(K)*XXLS(K)/CPM(K)
1103 QI3D(K)=0.
1104 END IF
1105 IF (QNI3D(K).LT.1.E-6) THEN
1106 QV3D(K)=QV3D(K)+QNI3D(K)
1107 T3D(K)=T3D(K)-QNI3D(K)*XXLS(K)/CPM(K)
1108 QNI3D(K)=0.
1109 END IF
1110 IF (QG3D(K).LT.1.E-6) THEN
1111 QV3D(K)=QV3D(K)+QG3D(K)
1112 T3D(K)=T3D(K)-QG3D(K)*XXLS(K)/CPM(K)
1113 QG3D(K)=0.
1114 END IF
1115 END IF
1117 ! HEAT OF FUSION
1119 XLF(K) = XXLS(K)-XXLV(K)
1121 !..................................................................
1122 ! IF MIXING RATIO < QSMALL SET MIXING RATIO AND NUMBER CONC TO ZERO
1124 IF (QC3D(K).LT.QSMALL) THEN
1125 QC3D(K) = 0.
1126 NC3D(K) = 0.
1127 EFFC(K) = 0.
1128 END IF
1129 IF (QR3D(K).LT.QSMALL) THEN
1130 QR3D(K) = 0.
1131 NR3D(K) = 0.
1132 EFFR(K) = 0.
1133 END IF
1134 IF (QI3D(K).LT.QSMALL) THEN
1135 QI3D(K) = 0.
1136 NI3D(K) = 0.
1137 EFFI(K) = 0.
1138 END IF
1139 IF (QNI3D(K).LT.QSMALL) THEN
1140 QNI3D(K) = 0.
1141 NS3D(K) = 0.
1142 EFFS(K) = 0.
1143 END IF
1144 IF (QG3D(K).LT.QSMALL) THEN
1145 QG3D(K) = 0.
1146 NG3D(K) = 0.
1147 EFFG(K) = 0.
1148 END IF
1150 ! INITIALIZE SEDIMENTATION TENDENCIES FOR MIXING RATIO
1152 QRSTEN(K) = 0.
1153 QISTEN(K) = 0.
1154 QNISTEN(K) = 0.
1155 QCSTEN(K) = 0.
1156 QGSTEN(K) = 0.
1158 !..................................................................
1159 ! MICROPHYSICS PARAMETERS VARYING IN TIME/HEIGHT
1161 ! FALL SPEED WITH DENSITY CORRECTION (HEYMSFIELD AND BENSSEMER 2006)
1163 DUM = (RHOSU/RHO(K))**0.54
1165 AIN(K) = DUM*AI
1166 ARN(K) = DUM*AR
1167 ASN(K) = DUM*AS
1168 ACN(K) = DUM*AC
1169 ! HM ADD GRAUPEL 8/28/06
1170 AGN(K) = DUM*AG
1172 !..................................
1173 ! IF THERE IS NO CLOUD/PRECIP WATER, AND IF SUBSATURATED, THEN SKIP MICROPHYSICS
1174 ! FOR THIS LEVEL
1176 IF (QC3D(K).LT.QSMALL.AND.QI3D(K).LT.QSMALL.AND.QNI3D(K).LT.QSMALL &
1177 .AND.QR3D(K).LT.QSMALL.AND.QG3D(K).LT.QSMALL) THEN
1178 IF (T3D(K).LT.273.15.AND.QVQVSI(K).LT.0.999) GOTO 200
1179 IF (T3D(K).GE.273.15.AND.QVQVS(K).LT.0.999) GOTO 200
1180 END IF
1182 ! THERMAL CONDUCTIVITY FOR AIR
1184 DUM = 1.496E-6*T3D(K)**1.5/(T3D(K)+120.)
1186 ! KAP(K) = 1.414E3*1.496E-6*T3D(K)**1.5/(T3D(K)+120.)
1187 KAP(K) = 1.414E3*DUM
1189 ! DIFFUSIVITY OF WATER VAPOR
1191 DV(K) = 8.794E-5*T3D(K)**1.81/PRES(K)
1193 ! VISCOSITY OF AIR
1194 ! SCHMIT NUMBER
1196 ! MU(K) = 1.496E-6*T3D(K)**1.5/(T3D(K)+120.)/RHO(K)
1197 MU(K) = DUM/RHO(K)
1198 SC(K) = MU(K)/DV(K)
1200 ! PSYCHOMETIC CORRECTIONS
1202 ! RATE OF CHANGE SAT. MIX. RATIO WITH TEMPERATURE
1204 DUM = (RV*T3D(K)**2)
1206 DQSDT = XXLV(K)*QVS(K)/DUM
1207 DQSIDT = XXLS(K)*QVI(K)/DUM
1209 ABI(K) = 1.+DQSIDT*XXLS(K)/CPM(K)
1210 AB(K) = 1.+DQSDT*XXLV(K)/CPM(K)
1212 !
1213 !.....................................................................
1214 !.....................................................................
1215 ! CASE FOR TEMPERATURE ABOVE FREEZING
1217 IF (T3D(K).GE.273.15) THEN
1219 !......................................................................
1220 !HM ADD, ALLOW FOR CONSTANT DROPLET NUMBER
1221 ! INUM = 0, PREDICT DROPLET NUMBER
1222 ! INUM = 1, SET CONSTANT DROPLET NUMBER
1224 ! IF (INUM.EQ.1) THEN
1225 ! CONVERT NDCNST FROM CM-3 TO KG-1
1226 NC3D(K)=NDCNST*1.E6/RHO(K)
1227 ! END IF
1229 ! GET SIZE DISTRIBUTION PARAMETERS
1231 ! MELT VERY SMALL SNOW AND GRAUPEL MIXING RATIOS, ADD TO RAIN
1232 IF (QNI3D(K).LT.1.E-6) THEN
1233 QR3D(K)=QR3D(K)+QNI3D(K)
1234 NR3D(K)=NR3D(K)+NS3D(K)
1235 T3D(K)=T3D(K)-QNI3D(K)*XLF(K)/CPM(K)
1236 QNI3D(K) = 0.
1237 NS3D(K) = 0.
1238 END IF
1239 IF (QG3D(K).LT.1.E-6) THEN
1240 QR3D(K)=QR3D(K)+QG3D(K)
1241 NR3D(K)=NR3D(K)+NG3D(K)
1242 T3D(K)=T3D(K)-QG3D(K)*XLF(K)/CPM(K)
1243 QG3D(K) = 0.
1244 NG3D(K) = 0.
1245 END IF
1247 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
1249 ! MAKE SURE NUMBER CONCENTRATIONS AREN'T NEGATIVE
1251 NS3D(K) = MAX(0.,NS3D(K))
1252 NC3D(K) = MAX(0.,NC3D(K))
1253 NR3D(K) = MAX(0.,NR3D(K))
1254 NG3D(K) = MAX(0.,NG3D(K))
1256 !......................................................................
1257 ! RAIN
1259 IF (QR3D(K).GE.QSMALL) THEN
1260 LAMR(K) = (PI*RHOW*NR3D(K)/QR3D(K))**(1./3.)
1261 N0RR(K) = NR3D(K)*LAMR(K)
1263 ! CHECK FOR SLOPE
1265 ! ADJUST VARS
1267 IF (LAMR(K).LT.LAMMINR) THEN
1269 LAMR(K) = LAMMINR
1271 N0RR(K) = LAMR(K)**4*QR3D(K)/(PI*RHOW)
1273 NR3D(K) = N0RR(K)/LAMR(K)
1274 ELSE IF (LAMR(K).GT.LAMMAXR) THEN
1275 LAMR(K) = LAMMAXR
1276 N0RR(K) = LAMR(K)**4*QR3D(K)/(PI*RHOW)
1278 NR3D(K) = N0RR(K)/LAMR(K)
1279 END IF
1280 END IF
1282 !......................................................................
1283 ! CLOUD DROPLETS
1285 ! MARTIN ET AL. (1994) FORMULA FOR PGAM
1287 IF (QC3D(K).GE.QSMALL) THEN
1289 DUM = PRES(K)/(287.15*T3D(K))
1290 PGAM(K)=0.0005714*(NC3D(K)/1.E6/DUM)+0.2714
1291 PGAM(K)=1./(PGAM(K)**2)-1.
1292 PGAM(K)=MAX(PGAM(K),2.)
1293 PGAM(K)=MIN(PGAM(K),10.)
1295 ! CALCULATE LAMC
1297 LAMC(K) = (CONS26*NC3D(K)*GAMMA(PGAM(K)+4.)/ &
1298 (QC3D(K)*GAMMA(PGAM(K)+1.)))**(1./3.)
1300 ! LAMMIN, 60 MICRON DIAMETER
1301 ! LAMMAX, 1 MICRON
1303 LAMMIN = (PGAM(K)+1.)/60.E-6
1304 LAMMAX = (PGAM(K)+1.)/1.E-6
1306 IF (LAMC(K).LT.LAMMIN) THEN
1307 LAMC(K) = LAMMIN
1309 NC3D(K) = EXP(3.*LOG(LAMC(K))+LOG(QC3D(K))+ &
1310 LOG(GAMMA(PGAM(K)+1.))-LOG(GAMMA(PGAM(K)+4.)))/CONS26
1311 ELSE IF (LAMC(K).GT.LAMMAX) THEN
1312 LAMC(K) = LAMMAX
1314 NC3D(K) = EXP(3.*LOG(LAMC(K))+LOG(QC3D(K))+ &
1315 LOG(GAMMA(PGAM(K)+1.))-LOG(GAMMA(PGAM(K)+4.)))/CONS26
1317 END IF
1319 END IF
1321 !......................................................................
1322 ! SNOW
1324 IF (QNI3D(K).GE.QSMALL) THEN
1325 LAMS(K) = (CONS1*NS3D(K)/QNI3D(K))**(1./DS)
1326 N0S(K) = NS3D(K)*LAMS(K)
1328 ! CHECK FOR SLOPE
1330 ! ADJUST VARS
1332 IF (LAMS(K).LT.LAMMINS) THEN
1333 LAMS(K) = LAMMINS
1334 N0S(K) = LAMS(K)**4*QNI3D(K)/CONS1
1336 NS3D(K) = N0S(K)/LAMS(K)
1338 ELSE IF (LAMS(K).GT.LAMMAXS) THEN
1340 LAMS(K) = LAMMAXS
1341 N0S(K) = LAMS(K)**4*QNI3D(K)/CONS1
1343 NS3D(K) = N0S(K)/LAMS(K)
1344 END IF
1345 END IF
1347 !......................................................................
1348 ! GRAUPEL
1350 IF (QG3D(K).GE.QSMALL) THEN
1351 LAMG(K) = (CONS2*NG3D(K)/QG3D(K))**(1./DG)
1352 N0G(K) = NG3D(K)*LAMG(K)
1354 ! ADJUST VARS
1356 IF (LAMG(K).LT.LAMMING) THEN
1357 LAMG(K) = LAMMING
1358 N0G(K) = LAMG(K)**4*QG3D(K)/CONS2
1360 NG3D(K) = N0G(K)/LAMG(K)
1362 ELSE IF (LAMG(K).GT.LAMMAXG) THEN
1364 LAMG(K) = LAMMAXG
1365 N0G(K) = LAMG(K)**4*QG3D(K)/CONS2
1367 NG3D(K) = N0G(K)/LAMG(K)
1368 END IF
1369 END IF
1371 !.....................................................................
1372 ! ZERO OUT PROCESS RATES
1374 PRC(K) = 0.
1375 NPRC(K) = 0.
1376 NPRC1(K) = 0.
1377 PRA(K) = 0.
1378 NPRA(K) = 0.
1379 NRAGG(K) = 0.
1380 PSMLT(K) = 0.
1381 NSMLTS(K) = 0.
1382 NSMLTR(K) = 0.
1383 EVPMS(K) = 0.
1384 PCC(K) = 0.
1385 PRE(K) = 0.
1386 NSUBC(K) = 0.
1387 NSUBR(K) = 0.
1388 PRACG(K) = 0.
1389 NPRACG(K) = 0.
1390 PSMLT(K) = 0.
1391 EVPMS(K) = 0.
1392 PGMLT(K) = 0.
1393 EVPMG(K) = 0.
1394 PRACS(K) = 0.
1395 NPRACS(K) = 0.
1396 NGMLTG(K) = 0.
1397 NGMLTR(K) = 0.
1399 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1400 ! CALCULATION OF MICROPHYSICAL PROCESS RATES, T > 273.15 K
1402 !.................................................................
1403 !.......................................................................
1404 ! AUTOCONVERSION OF CLOUD LIQUID WATER TO RAIN
1405 ! FORMULA FROM BEHENG (1994)
1406 ! USING NUMERICAL SIMULATION OF STOCHASTIC COLLECTION EQUATION
1407 ! AND INITIAL CLOUD DROPLET SIZE DISTRIBUTION SPECIFIED
1408 ! AS A GAMMA DISTRIBUTION
1410 ! USE MINIMUM VALUE OF 1.E-6 TO PREVENT FLOATING POINT ERROR
1412 IF (QC3D(K).GE.1.E-6) THEN
1414 ! HM ADD 12/13/06, REPLACE WITH NEWER FORMULA
1415 ! FROM KHAIROUTDINOV AND KOGAN 2000, MWR
1417 PRC(K)=1350.*QC3D(K)**2.47* &
1418 (NC3D(K)/1.e6*RHO(K))**(-1.79)
1420 ! note: nprc1 is change in Nr,
1421 ! nprc is change in Nc
1423 NPRC1(K) = PRC(K)/CONS29
1424 NPRC(K) = PRC(K)/(QC3D(k)/NC3D(K))
1426 NPRC(K) = MIN(NPRC(K),NC3D(K)/DT)
1428 END IF
1430 !.......................................................................
1431 ! HM ADD 12/13/06, COLLECTION OF SNOW BY RAIN ABOVE FREEZING
1432 ! FORMULA FROM IKAWA AND SAITO (1991)
1434 IF (QR3D(K).GE.1.E-8.AND.QNI3D(K).GE.1.E-8) THEN
1436 UMS = ASN(K)*CONS3/(LAMS(K)**BS)
1437 UMR = ARN(K)*CONS4/(LAMR(K)**BR)
1438 UNS = ASN(K)*CONS5/LAMS(K)**BS
1439 UNR = ARN(K)*CONS6/LAMR(K)**BR
1441 ! SET REASLISTIC LIMITS ON FALLSPEEDS
1442 UMS=MIN(UMS,1.2)
1443 UNS=MIN(UNS,1.2)
1444 UMR=MIN(UMR,9.1)
1445 UNR=MIN(UNR,9.1)
1447 PRACS(K) = CONS31*(((1.2*UMR-0.95*UMS)**2+ &
1448 0.08*UMS*UMR)**0.5*RHO(K)* &
1449 N0RR(K)*N0S(K)/LAMS(K)**3* &
1450 (5./(LAMS(K)**3*LAMR(K))+ &
1451 2./(LAMS(K)**2*LAMR(K)**2)+ &
1452 0.5/(LAMS(K)*LAMR(K)**3)))
1454 NPRACS(K) = CONS32*RHO(K)*(1.7*(UNR-UNS)**2+ &
1455 0.3*UNR*UNS)**0.5*N0RR(K)*N0S(K)* &
1456 (1./(LAMR(K)**3*LAMS(K))+ &
1457 1./(LAMR(K)**2*LAMS(K)**2)+ &
1458 1./(LAMR(K)*LAMS(K)**3))
1460 END IF
1462 ! ADD COLLECTION OF GRAUPEL BY RAIN ABOVE FREEZING
1463 ! ASSUME ALL RAIN COLLECTION BY GRAUPEL ABOVE FREEZING IS SHED
1464 ! ASSUME SHED DROPS ARE 1 MM IN SIZE
1466 IF (QR3D(K).GE.1.E-8.AND.QG3D(K).GE.1.E-8) THEN
1468 UMG = AGN(K)*CONS7/(LAMG(K)**BG)
1469 UMR = ARN(K)*CONS4/(LAMR(K)**BR)
1470 UNG = AGN(K)*CONS8/LAMG(K)**BG
1471 UNR = ARN(K)*CONS6/LAMR(K)**BR
1473 ! SET REASLISTIC LIMITS ON FALLSPEEDS
1474 UMG=MIN(UMG,20.)
1475 UNG=MIN(UNG,20.)
1476 UMR=MIN(UMR,9.1)
1477 UNR=MIN(UNR,9.1)
1479 ! DUM IS MIXING RATIO OF RAIN PER SEC COLLECTED BY GRAUPEL/HAIL
1480 DUM = CONS41*(((1.2*UMR-0.95*UMG)**2+ &
1481 0.08*UMG*UMR)**0.5*RHO(K)* &
1482 N0RR(K)*N0G(K)/LAMR(K)**3* &
1483 (5./(LAMR(K)**3*LAMG(K))+ &
1484 2./(LAMR(K)**2*LAMG(K)**2)+ &
1485 0.5/(LAMR(k)*LAMG(k)**3)))
1487 ! ASSUME 1 MM DROPS ARE SHED, GET NUMBER SHED PER SEC
1489 DUM = DUM/5.2E-7
1491 NPRACG(K) = CONS32*RHO(K)*(1.7*(UNR-UNG)**2+ &
1492 0.3*UNR*UNG)**0.5*N0RR(K)*N0G(K)* &
1493 (1./(LAMR(K)**3*LAMG(K))+ &
1494 1./(LAMR(K)**2*LAMG(K)**2)+ &
1495 1./(LAMR(K)*LAMG(K)**3))
1497 NPRACG(K)=MAX(NPRACG(K)-DUM,0.)
1499 END IF
1501 !.......................................................................
1502 ! ACCRETION OF CLOUD LIQUID WATER BY RAIN
1503 ! CONTINUOUS COLLECTION EQUATION WITH
1504 ! GRAVITATIONAL COLLECTION KERNEL, DROPLET FALL SPEED NEGLECTED
1506 IF (QR3D(K).GE.1.E-8 .AND. QC3D(K).GE.1.E-8) THEN
1508 ! 12/13/06 HM ADD, REPLACE WITH NEWER FORMULA FROM
1509 ! KHAIROUTDINOV AND KOGAN 2000, MWR
1511 DUM=(QC3D(K)*QR3D(K))
1512 PRA(K) = 67.*(DUM)**1.15
1513 NPRA(K) = PRA(K)/(QC3D(K)/NC3D(K))
1515 END IF
1516 !.......................................................................
1517 ! SELF-COLLECTION OF RAIN DROPS
1518 ! FROM BEHENG(1994)
1519 ! FROM NUMERICAL SIMULATION OF THE STOCHASTIC COLLECTION EQUATION
1520 ! AS DESCRINED ABOVE FOR AUTOCONVERSION
1522 IF (QR3D(K).GE.1.E-8) THEN
1523 NRAGG(K) = -8.*NR3D(K)*QR3D(K)*RHO(K)
1524 END IF
1526 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1527 ! CALCULATE EVAP OF RAIN (RUTLEDGE AND HOBBS 1983)
1529 IF (QR3D(K).GE.QSMALL) THEN
1530 EPSR = 2.*PI*N0RR(K)*RHO(K)*DV(K)* &
1531 (F1R/(LAMR(K)*LAMR(K))+ &
1532 F2R*(ARN(K)*RHO(K)/MU(K))**0.5* &
1533 SC(K)**(1./3.)*CONS9/ &
1534 (LAMR(K)**CONS34))
1535 ELSE
1536 EPSR = 0.
1537 END IF
1539 ! NO CONDENSATION ONTO RAIN, ONLY EVAP ALLOWED
1541 IF (QV3D(K).LT.QVS(K)) THEN
1542 PRE(K) = EPSR*(QV3D(K)-QVS(K))/AB(K)
1543 PRE(K) = MIN(PRE(K),0.)
1544 ELSE
1545 PRE(K) = 0.
1546 END IF
1548 !.......................................................................
1549 ! MELTING OF SNOW
1551 ! SNOW MAY PERSITS ABOVE FREEZING, FORMULA FROM RUTLEDGE AND HOBBS, 1984
1552 ! IF WATER SUPERSATURATION, SNOW MELTS TO FORM RAIN
1554 IF (QNI3D(K).GE.1.E-8) THEN
1556 PSMLT(K)=2.*PI*N0S(K)*KAP(K)*(273.15-T3D(K))/ &
1557 XLF(K)*RHO(K)*(F1S/(LAMS(K)*LAMS(K))+ &
1558 F2S*(ASN(K)*RHO(K)/MU(K))**0.5* &
1559 SC(K)**(1./3.)*CONS10/ &
1560 (LAMS(K)**CONS35))
1562 ! IN WATER SUBSATURATION, SNOW MELTS AND EVAPORATES
1564 IF (QVQVS(K).LT.1.) THEN
1565 EPSS = 2.*PI*N0S(K)*RHO(K)*DV(K)* &
1566 (F1S/(LAMS(K)*LAMS(K))+ &
1567 F2S*(ASN(K)*RHO(K)/MU(K))**0.5* &
1568 SC(K)**(1./3.)*CONS10/ &
1569 (LAMS(K)**CONS35))
1570 EVPMS(K) = 2.*PI*N0S(K)*(QV3D(K)-QVS(K))*EPSS/AB(K)
1571 EVPMS(K) = MAX(EVPMS(K),PSMLT(K))
1572 PSMLT(K) = PSMLT(K)-EVPMS(K)
1573 END IF
1574 END IF
1576 !.......................................................................
1577 ! MELTING OF GRAUPEL
1579 ! GRAUPEL MAY PERSITS ABOVE FREEZING, FORMULA FROM RUTLEDGE AND HOBBS, 1984
1580 ! IF WATER SUPERSATURATION, GRAUPEL MELTS TO FORM RAIN
1582 IF (QG3D(K).GE.1.E-8) THEN
1584 PGMLT(K)=2.*PI*N0G(K)*KAP(K)*(273.15-T3D(K))/ &
1585 XLF(K)*RHO(K)*(F1S/(LAMG(K)*LAMG(K))+ &
1586 F2S*(AGN(K)*RHO(K)/MU(K))**0.5* &
1587 SC(K)**(1./3.)*CONS11/ &
1588 (LAMG(K)**CONS36))
1590 ! IN WATER SUBSATURATION, GRAUPEL MELTS AND EVAPORATES
1592 IF (QVQVS(K).LT.1.) THEN
1593 EPSG = 2.*PI*N0G(K)*RHO(K)*DV(K)* &
1594 (F1S/(LAMG(K)*LAMG(K))+ &
1595 F2S*(AGN(K)*RHO(K)/MU(K))**0.5* &
1596 SC(K)**(1./3.)*CONS11/ &
1597 (LAMG(K)**CONS36))
1598 EVPMG(K) = 2.*PI*N0G(K)*(QV3D(K)-QVS(K))*EPSG/AB(K)
1599 EVPMG(K) = MAX(EVPMG(K),PGMLT(K))
1600 PGMLT(K) = PGMLT(K)-EVPMG(K)
1601 END IF
1602 END IF
1604 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1606 ! FOR CLOUD ICE, ONLY PROCESSES OPERATING AT T > 273.15 IS
1607 ! MELTING, WHICH IS ALREADY CONSERVED DURING PROCESS
1608 ! CALCULATION
1610 ! CONSERVATION OF QC
1612 DUM = (PRC(K)+PRA(K))*DT
1614 IF (DUM.GT.QC3D(K).AND.QC3D(K).GE.QSMALL) THEN
1616 RATIO = QC3D(K)/DUM
1618 PRC(K) = PRC(K)*RATIO
1619 PRA(K) = PRA(K)*RATIO
1621 END IF
1623 ! CONSERVATION OF SNOW
1625 DUM = (-PSMLT(K)-EVPMS(K)+PRACS(K))*DT
1627 IF (DUM.GT.QNI3D(K).AND.QNI3D(K).GE.QSMALL) THEN
1629 ! NO SOURCE TERMS FOR SNOW AT T > FREEZING
1630 RATIO = QNI3D(K)/DUM
1632 PSMLT(K) = PSMLT(K)*RATIO
1633 EVPMS(K) = EVPMS(K)*RATIO
1634 PRACS(K) = PRACS(K)*RATIO
1636 END IF
1638 ! CONSERVATION OF GRAUPEL
1640 DUM = (-PGMLT(K)-EVPMG(K)+PRACG(K))*DT
1642 IF (DUM.GT.QG3D(K).AND.QG3D(K).GE.QSMALL) THEN
1644 ! NO SOURCE TERM FOR GRAUPEL ABOVE FREEZING
1645 RATIO = QG3D(K)/DUM
1647 PGMLT(K) = PGMLT(K)*RATIO
1648 EVPMG(K) = EVPMG(K)*RATIO
1649 PRACG(K) = PRACG(K)*RATIO
1651 END IF
1653 ! CONSERVATION OF QR
1654 ! HM 12/13/06, ADDED CONSERVATION OF RAIN SINCE PRE IS NEGATIVE
1656 DUM = (-PRACS(K)-PRACG(K)-PRE(K)-PRA(K)-PRC(K)+PSMLT(K)+PGMLT(K))*DT
1658 IF (DUM.GT.QR3D(K).AND.QR3D(K).GE.QSMALL) THEN
1660 RATIO = (QR3D(K)/DT+PRACS(K)+PRACG(K)+PRA(K)+PRC(K)-PSMLT(K)-PGMLT(K))/ &
1661 (-PRE(K))
1662 PRE(K) = PRE(K)*RATIO
1664 END IF
1666 !....................................
1668 QV3DTEN(K) = QV3DTEN(K)+(-PRE(K)-EVPMS(K)-EVPMG(K))
1670 T3DTEN(K) = T3DTEN(K)+(PRE(K)*XXLV(K)+(EVPMS(K)+EVPMG(K))*XXLS(K)+&
1671 (PSMLT(K)+PGMLT(K)-PRACS(K)-PRACG(K))*XLF(K))/CPM(K)
1673 QC3DTEN(K) = QC3DTEN(K)+(-PRA(K)-PRC(K))
1674 QR3DTEN(K) = QR3DTEN(K)+(PRE(K)+PRA(K)+PRC(K)-PSMLT(K)-PGMLT(K)+PRACS(K)+PRACG(K))
1675 QNI3DTEN(K) = QNI3DTEN(K)+(PSMLT(K)+EVPMS(K)-PRACS(K))
1676 QG3DTEN(K) = QG3DTEN(K)+(PGMLT(K)+EVPMG(K)-PRACG(K))
1677 NS3DTEN(K) = NS3DTEN(K)-NPRACS(K)
1678 ! HM, bug fix 5/12/08, npracg is subtracted from nr not ng
1679 ! NG3DTEN(K) = NG3DTEN(K)
1680 NC3DTEN(K) = NC3DTEN(K)+ (-NPRA(K)-NPRC(K))
1681 NR3DTEN(K) = NR3DTEN(K)+ (NPRC1(K)+NRAGG(K)-NPRACG(K))
1683 IF (PRE(K).LT.0.) THEN
1684 DUM = PRE(K)*DT/QR3D(K)
1685 DUM = MAX(-1.,DUM)
1686 NSUBR(K) = DUM*NR3D(K)/DT
1687 END IF
1689 IF (EVPMS(K)+PSMLT(K).LT.0.) THEN
1690 DUM = (EVPMS(K)+PSMLT(K))*DT/QNI3D(K)
1691 DUM = MAX(-1.,DUM)
1692 NSMLTS(K) = DUM*NS3D(K)/DT
1693 END IF
1694 IF (PSMLT(K).LT.0.) THEN
1695 DUM = PSMLT(K)*DT/QNI3D(K)
1696 DUM = MAX(-1.0,DUM)
1697 NSMLTR(K) = DUM*NS3D(K)/DT
1698 END IF
1699 IF (EVPMG(K)+PGMLT(K).LT.0.) THEN
1700 DUM = (EVPMG(K)+PGMLT(K))*DT/QG3D(K)
1701 DUM = MAX(-1.,DUM)
1702 NGMLTG(K) = DUM*NG3D(K)/DT
1703 END IF
1704 IF (PGMLT(K).LT.0.) THEN
1705 DUM = PGMLT(K)*DT/QG3D(K)
1706 DUM = MAX(-1.0,DUM)
1707 NGMLTR(K) = DUM*NG3D(K)/DT
1708 END IF
1710 NS3DTEN(K) = NS3DTEN(K)+(NSMLTS(K))
1711 NG3DTEN(K) = NG3DTEN(K)+(NGMLTG(K))
1712 NR3DTEN(K) = NR3DTEN(K)+(NSUBR(K)-NSMLTR(K)-NGMLTR(K))
1714 300 CONTINUE
1716 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1717 ! NOW CALCULATE SATURATION ADJUSTMENT TO CONDENSE EXTRA VAPOR ABOVE
1718 ! WATER SATURATION
1720 DUMT = T3D(K)+DT*T3DTEN(K)
1721 DUMQV = QV3D(K)+DT*QV3DTEN(K)
1722 DUMQSS = 0.622*POLYSVP(DUMT,0)/ (PRES(K)-POLYSVP(DUMT,0))
1723 DUMQC = QC3D(K)+DT*QC3DTEN(K)
1724 DUMQC = MAX(DUMQC,0.)
1726 ! SATURATION ADJUSTMENT FOR LIQUID
1728 DUMS = DUMQV-DUMQSS
1729 PCC(K) = DUMS/(1.+XXLV(K)**2*DUMQSS/(CPM(K)*RV*DUMT**2))/DT
1730 IF (PCC(K)*DT+DUMQC.LT.0.) THEN
1731 PCC(K) = -DUMQC/DT
1732 END IF
1734 QV3DTEN(K) = QV3DTEN(K)-PCC(K)
1735 T3DTEN(K) = T3DTEN(K)+PCC(K)*XXLV(K)/CPM(K)
1736 QC3DTEN(K) = QC3DTEN(K)+PCC(K)
1738 !.......................................................................
1739 ! ACTIVATION OF CLOUD DROPLETS
1740 ! ACTIVATION OF DROPLET CURRENTLY NOT CALCULATED
1741 ! DROPLET CONCENTRATION IS SPECIFIED !!!!!
1743 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1744 ! SUBLIMATE, MELT, OR EVAPORATE NUMBER CONCENTRATION
1745 ! THIS FORMULATION ASSUMES 1:1 RATIO BETWEEN MASS LOSS AND
1746 ! LOSS OF NUMBER CONCENTRATION
1748 ! IF (PCC(K).LT.0.) THEN
1749 ! DUM = PCC(K)*DT/QC3D(K)
1750 ! DUM = MAX(-1.,DUM)
1751 ! NSUBC(K) = DUM*NC3D(K)/DT
1752 ! END IF
1754 ! UPDATE TENDENCIES
1756 ! NC3DTEN(K) = NC3DTEN(K)+NSUBC(K)
1758 !.....................................................................
1759 !.....................................................................
1760 ELSE ! TEMPERATURE < 273.15
1762 !......................................................................
1763 !HM ADD, ALLOW FOR CONSTANT DROPLET NUMBER
1764 ! INUM = 0, PREDICT DROPLET NUMBER
1765 ! INUM = 1, SET CONSTANT DROPLET NUMBER
1767 ! IF (INUM.EQ.1) THEN
1768 ! CONVERT NDCNST FROM CM-3 TO KG-1
1769 NC3D(K)=NDCNST*1.E6/RHO(K)
1770 ! END IF
1772 ! CALCULATE SIZE DISTRIBUTION PARAMETERS
1773 ! MAKE SURE NUMBER CONCENTRATIONS AREN'T NEGATIVE
1775 NI3D(K) = MAX(0.,NI3D(K))
1776 NS3D(K) = MAX(0.,NS3D(K))
1777 NC3D(K) = MAX(0.,NC3D(K))
1778 NR3D(K) = MAX(0.,NR3D(K))
1779 NG3D(K) = MAX(0.,NG3D(K))
1781 !......................................................................
1782 ! CLOUD ICE
1784 IF (QI3D(K).GE.QSMALL) THEN
1785 LAMI(K) = (CONS12* &
1786 NI3D(K)/QI3D(K))**(1./DI)
1787 N0I(K) = NI3D(K)*LAMI(K)
1789 ! CHECK FOR SLOPE
1791 ! ADJUST VARS
1793 IF (LAMI(K).LT.LAMMINI) THEN
1795 LAMI(K) = LAMMINI
1797 N0I(K) = LAMI(K)**4*QI3D(K)/CONS12
1799 NI3D(K) = N0I(K)/LAMI(K)
1800 ELSE IF (LAMI(K).GT.LAMMAXI) THEN
1801 LAMI(K) = LAMMAXI
1802 N0I(K) = LAMI(K)**4*QI3D(K)/CONS12
1804 NI3D(K) = N0I(K)/LAMI(K)
1805 END IF
1806 END IF
1808 !......................................................................
1809 ! RAIN
1811 IF (QR3D(K).GE.QSMALL) THEN
1812 LAMR(K) = (PI*RHOW*NR3D(K)/QR3D(K))**(1./3.)
1813 N0RR(K) = NR3D(K)*LAMR(K)
1815 ! CHECK FOR SLOPE
1817 ! ADJUST VARS
1819 IF (LAMR(K).LT.LAMMINR) THEN
1821 LAMR(K) = LAMMINR
1823 N0RR(K) = LAMR(K)**4*QR3D(K)/(PI*RHOW)
1825 NR3D(K) = N0RR(K)/LAMR(K)
1826 ELSE IF (LAMR(K).GT.LAMMAXR) THEN
1827 LAMR(K) = LAMMAXR
1828 N0RR(K) = LAMR(K)**4*QR3D(K)/(PI*RHOW)
1830 NR3D(K) = N0RR(K)/LAMR(K)
1831 END IF
1832 END IF
1834 !......................................................................
1835 ! CLOUD DROPLETS
1837 ! MARTIN ET AL. (1994) FORMULA FOR PGAM
1839 IF (QC3D(K).GE.QSMALL) THEN
1841 DUM = PRES(K)/(287.15*T3D(K))
1842 PGAM(K)=0.0005714*(NC3D(K)/1.E6/DUM)+0.2714
1843 PGAM(K)=1./(PGAM(K)**2)-1.
1844 PGAM(K)=MAX(PGAM(K),2.)
1845 PGAM(K)=MIN(PGAM(K),10.)
1847 ! CALCULATE LAMC
1849 LAMC(K) = (CONS26*NC3D(K)*GAMMA(PGAM(K)+4.)/ &
1850 (QC3D(K)*GAMMA(PGAM(K)+1.)))**(1./3.)
1852 ! LAMMIN, 60 MICRON DIAMETER
1853 ! LAMMAX, 1 MICRON
1855 LAMMIN = (PGAM(K)+1.)/60.E-6
1856 LAMMAX = (PGAM(K)+1.)/1.E-6
1858 IF (LAMC(K).LT.LAMMIN) THEN
1859 LAMC(K) = LAMMIN
1861 NC3D(K) = EXP(3.*LOG(LAMC(K))+LOG(QC3D(K))+ &
1862 LOG(GAMMA(PGAM(K)+1.))-LOG(GAMMA(PGAM(K)+4.)))/CONS26
1863 ELSE IF (LAMC(K).GT.LAMMAX) THEN
1864 LAMC(K) = LAMMAX
1865 NC3D(K) = EXP(3.*LOG(LAMC(K))+LOG(QC3D(K))+ &
1866 LOG(GAMMA(PGAM(K)+1.))-LOG(GAMMA(PGAM(K)+4.)))/CONS26
1868 END IF
1870 ! TO CALCULATE DROPLET FREEZING
1872 CDIST1(K) = NC3D(K)/GAMMA(PGAM(K)+1.)
1874 END IF
1876 !......................................................................
1877 ! SNOW
1879 IF (QNI3D(K).GE.QSMALL) THEN
1880 LAMS(K) = (CONS1*NS3D(K)/QNI3D(K))**(1./DS)
1881 N0S(K) = NS3D(K)*LAMS(K)
1883 ! CHECK FOR SLOPE
1885 ! ADJUST VARS
1887 IF (LAMS(K).LT.LAMMINS) THEN
1888 LAMS(K) = LAMMINS
1889 N0S(K) = LAMS(K)**4*QNI3D(K)/CONS1
1891 NS3D(K) = N0S(K)/LAMS(K)
1893 ELSE IF (LAMS(K).GT.LAMMAXS) THEN
1895 LAMS(K) = LAMMAXS
1896 N0S(K) = LAMS(K)**4*QNI3D(K)/CONS1
1898 NS3D(K) = N0S(K)/LAMS(K)
1899 END IF
1900 END IF
1902 !......................................................................
1903 ! GRAUPEL
1905 IF (QG3D(K).GE.QSMALL) THEN
1906 LAMG(K) = (CONS2*NG3D(K)/QG3D(K))**(1./DG)
1907 N0G(K) = NG3D(K)*LAMG(K)
1909 ! CHECK FOR SLOPE
1911 ! ADJUST VARS
1913 IF (LAMG(K).LT.LAMMING) THEN
1914 LAMG(K) = LAMMING
1915 N0G(K) = LAMG(K)**4*QG3D(K)/CONS2
1917 NG3D(K) = N0G(K)/LAMG(K)
1919 ELSE IF (LAMG(K).GT.LAMMAXG) THEN
1921 LAMG(K) = LAMMAXG
1922 N0G(K) = LAMG(K)**4*QG3D(K)/CONS2
1924 NG3D(K) = N0G(K)/LAMG(K)
1925 END IF
1926 END IF
1928 !.....................................................................
1929 ! ZERO OUT PROCESS RATES
1931 MNUCCC(K) = 0.
1932 NNUCCC(K) = 0.
1933 PRC(K) = 0.
1934 NPRC(K) = 0.
1935 NPRC1(K) = 0.
1936 NSAGG(K) = 0.
1937 PSACWS(K) = 0.
1938 NPSACWS(K) = 0.
1939 PSACWI(K) = 0.
1940 NPSACWI(K) = 0.
1941 PRACS(K) = 0.
1942 NPRACS(K) = 0.
1943 NMULTS(K) = 0.
1944 QMULTS(K) = 0.
1945 NMULTR(K) = 0.
1946 QMULTR(K) = 0.
1947 NMULTG(K) = 0.
1948 QMULTG(K) = 0.
1949 NMULTRG(K) = 0.
1950 QMULTRG(K) = 0.
1951 MNUCCR(K) = 0.
1952 NNUCCR(K) = 0.
1953 PRA(K) = 0.
1954 NPRA(K) = 0.
1955 NRAGG(K) = 0.
1956 PRCI(K) = 0.
1957 NPRCI(K) = 0.
1958 PRAI(K) = 0.
1959 NPRAI(K) = 0.
1960 NNUCCD(K) = 0.
1961 MNUCCD(K) = 0.
1962 PCC(K) = 0.
1963 PRE(K) = 0.
1964 PRD(K) = 0.
1965 PRDS(K) = 0.
1966 EPRD(K) = 0.
1967 EPRDS(K) = 0.
1968 NSUBC(K) = 0.
1969 NSUBI(K) = 0.
1970 NSUBS(K) = 0.
1971 NSUBR(K) = 0.
1972 PIACR(K) = 0.
1973 NIACR(K) = 0.
1974 PRACI(K) = 0.
1975 PIACRS(K) = 0.
1976 NIACRS(K) = 0.
1977 PRACIS(K) = 0.
1978 ! HM: ADD GRAUPEL PROCESSES
1979 PRACG(K) = 0.
1980 PSACR(K) = 0.
1981 PSACWG(K) = 0.
1982 PGSACW(K) = 0.
1983 PGRACS(K) = 0.
1984 PRDG(K) = 0.
1985 EPRDG(K) = 0.
1986 NPRACG(K) = 0.
1987 NPSACWG(K) = 0.
1988 NSCNG(K) = 0.
1989 NGRACS(K) = 0.
1990 NSUBG(K) = 0.
1992 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1993 ! CALCULATION OF MICROPHYSICAL PROCESS RATES
1994 ! ACCRETION/AUTOCONVERSION/FREEZING/MELTING/COAG.
1995 !.......................................................................
1996 ! FREEZING OF CLOUD DROPLETS
1997 ! ONLY ALLOWED BELOW -4 C
1998 IF (QC3D(K).GE.QSMALL .AND. T3D(K).LT.269.15) THEN
2000 ! NUMBER OF CONTACT NUCLEI (M^-3) FROM MEYERS ET AL., 1992
2001 ! FACTOR OF 1000 IS TO CONVERT FROM L^-1 TO M^-3
2003 ! MEYERS CURVE
2005 NACNT = EXP(-2.80+0.262*(273.15-T3D(K)))*1000.
2007 ! COOPER CURVE
2008 ! NACNT = 5.*EXP(0.304*(273.15-T3D(K)))
2010 ! FLECTHER
2011 ! NACNT = 0.01*EXP(0.6*(273.15-T3D(K)))
2013 ! CONTACT FREEZING
2015 ! MEAN FREE PATH
2017 DUM = 7.37*T3D(K)/(288.*10.*PRES(K))/100.
2019 ! EFFECTIVE DIFFUSIVITY OF CONTACT NUCLEI
2020 ! BASED ON BROWNIAN DIFFUSION
2022 DAP(K) = CONS37*T3D(K)*(1.+DUM/RIN)/MU(K)
2024 MNUCCC(K) = CONS38*DAP(K)*NACNT*EXP(LOG(CDIST1(K))+ &
2025 LOG(GAMMA(PGAM(K)+5.))-4.*LOG(LAMC(K)))
2026 NNUCCC(K) = 2.*PI*DAP(K)*NACNT*CDIST1(K)* &
2027 GAMMA(PGAM(K)+2.)/ &
2028 LAMC(K)
2030 ! IMMERSION FREEZING (BIGG 1953)
2032 MNUCCC(K) = MNUCCC(K)+CONS39* &
2033 EXP(LOG(CDIST1(K))+LOG(GAMMA(7.+PGAM(K)))-6.*LOG(LAMC(K)))* &
2034 EXP(AIMM*(273.15-T3D(K)))
2036 NNUCCC(K) = NNUCCC(K)+ &
2037 CONS40*EXP(LOG(CDIST1(K))+LOG(GAMMA(PGAM(K)+4.))-3.*LOG(LAMC(K))) &
2038 *EXP(AIMM*(273.15-T3D(K)))
2040 ! PUT IN A CATCH HERE TO PREVENT DIVERGENCE BETWEEN NUMBER CONC. AND
2041 ! MIXING RATIO, SINCE STRICT CONSERVATION NOT CHECKED FOR NUMBER CONC
2043 NNUCCC(K) = MIN(NNUCCC(K),NC3D(K)/DT)
2045 END IF
2047 !.................................................................
2048 !.......................................................................
2049 ! AUTOCONVERSION OF CLOUD LIQUID WATER TO RAIN
2050 ! FORMULA FROM BEHENG (1994)
2051 ! USING NUMERICAL SIMULATION OF STOCHASTIC COLLECTION EQUATION
2052 ! AND INITIAL CLOUD DROPLET SIZE DISTRIBUTION SPECIFIED
2053 ! AS A GAMMA DISTRIBUTION
2055 ! USE MINIMUM VALUE OF 1.E-6 TO PREVENT FLOATING POINT ERROR
2057 IF (QC3D(K).GE.1.E-6) THEN
2059 ! HM ADD 12/13/06, REPLACE WITH NEWER FORMULA
2060 ! FROM KHAIROUTDINOV AND KOGAN 2000, MWR
2062 PRC(K)=1350.*QC3D(K)**2.47* &
2063 (NC3D(K)/1.e6*RHO(K))**(-1.79)
2065 ! note: nprc1 is change in Nr,
2066 ! nprc is change in Nc
2068 NPRC1(K) = PRC(K)/CONS29
2069 NPRC(K) = PRC(K)/(QC3D(K)/NC3D(K))
2071 NPRC(K) = MIN(NPRC(K),NC3D(K)/DT)
2073 END IF
2075 !.......................................................................
2076 ! SELF-COLLECTION OF DROPLET NOT INCLUDED IN KK2000 SCHEME
2078 ! SNOW AGGREGATION FROM PASSARELLI, 1978, USED BY REISNER, 1998
2079 ! THIS IS HARD-WIRED FOR BS = 0.4 FOR NOW
2081 IF (QNI3D(K).GE.1.E-8) THEN
2082 NSAGG(K) = CONS15*ASN(K)*RHO(K)** &
2083 ((2.+BS)/3.)*QNI3D(K)**((2.+BS)/3.)* &
2084 (NS3D(K)*RHO(K))**((4.-BS)/3.)/ &
2085 (RHO(K))
2086 END IF
2088 !.......................................................................
2089 ! ACCRETION OF CLOUD DROPLETS ONTO SNOW/GRAUPEL
2090 ! HERE USE CONTINUOUS COLLECTION EQUATION WITH
2091 ! SIMPLE GRAVITATIONAL COLLECTION KERNEL IGNORING
2093 ! SNOW
2095 IF (QNI3D(K).GE.1.E-8 .AND. QC3D(K).GE.QSMALL) THEN
2097 PSACWS(K) = CONS13*ASN(K)*QC3D(K)*RHO(K)* &
2098 N0S(K)/ &
2099 LAMS(K)**(BS+3.)
2100 NPSACWS(K) = CONS13*ASN(K)*NC3D(K)*RHO(K)* &
2101 N0S(K)/ &
2102 LAMS(K)**(BS+3.)
2104 END IF
2106 !............................................................................
2107 ! COLLECTION OF CLOUD WATER BY GRAUPEL
2109 IF (QG3D(K).GE.1.E-8 .AND. QC3D(K).GE.QSMALL) THEN
2111 PSACWG(K) = CONS14*AGN(K)*QC3D(K)*RHO(K)* &
2112 N0G(K)/ &
2113 LAMG(K)**(BG+3.)
2114 NPSACWG(K) = CONS14*AGN(K)*NC3D(K)*RHO(K)* &
2115 N0G(K)/ &
2116 LAMG(K)**(BG+3.)
2117 END IF
2119 !.......................................................................
2120 ! HM, ADD 12/13/06
2121 ! CLOUD ICE COLLECTING DROPLETS, ASSUME THAT CLOUD ICE MEAN DIAM > 100 MICRON
2122 ! BEFORE RIMING CAN OCCUR
2123 ! ASSUME THAT RIME COLLECTED ON CLOUD ICE DOES NOT LEAD
2124 ! TO HALLET-MOSSOP SPLINTERING
2126 IF (QI3D(K).GE.1.E-8 .AND. QC3D(K).GE.QSMALL) THEN
2128 ! PUT IN SIZE DEPENDENT COLLECTION EFFICIENCY BASED ON STOKES LAW
2129 ! FROM THOMPSON ET AL. 2004, MWR
2131 IF (1./LAMI(K).GE.100.E-6) THEN
2133 PSACWI(K) = CONS16*AIN(K)*QC3D(K)*RHO(K)* &
2134 N0I(K)/ &
2135 LAMI(K)**(BI+3.)
2136 NPSACWI(K) = CONS16*AIN(K)*NC3D(K)*RHO(K)* &
2137 N0I(K)/ &
2138 LAMI(K)**(BI+3.)
2139 END IF
2140 END IF
2142 !.......................................................................
2143 ! ACCRETION OF RAIN WATER BY SNOW
2144 ! FORMULA FROM IKAWA AND SAITO, 1991, USED BY REISNER ET AL, 1998
2146 IF (QR3D(K).GE.1.E-8.AND.QNI3D(K).GE.1.E-8) THEN
2148 UMS = ASN(K)*CONS3/(LAMS(K)**BS)
2149 UMR = ARN(K)*CONS4/(LAMR(K)**BR)
2150 UNS = ASN(K)*CONS5/LAMS(K)**BS
2151 UNR = ARN(K)*CONS6/LAMR(K)**BR
2153 ! SET REASLISTIC LIMITS ON FALLSPEEDS
2154 UMS=MIN(UMS,1.2)
2155 UNS=MIN(UNS,1.2)
2156 UMR=MIN(UMR,9.1)
2157 UNR=MIN(UNR,9.1)
2159 PRACS(K) = CONS41*(((1.2*UMR-0.95*UMS)**2+ &
2160 0.08*UMS*UMR)**0.5*RHO(K)* &
2161 N0RR(K)*N0S(K)/LAMR(K)**3* &
2162 (5./(LAMR(K)**3*LAMS(K))+ &
2163 2./(LAMR(K)**2*LAMS(K)**2)+ &
2164 0.5/(LAMR(k)*LAMS(k)**3)))
2166 NPRACS(K) = CONS32*RHO(K)*(1.7*(UNR-UNS)**2+ &
2167 0.3*UNR*UNS)**0.5*N0RR(K)*N0S(K)* &
2168 (1./(LAMR(K)**3*LAMS(K))+ &
2169 1./(LAMR(K)**2*LAMS(K)**2)+ &
2170 1./(LAMR(K)*LAMS(K)**3))
2172 ! MAKE SURE PRACS DOESN'T EXCEED TOTAL RAIN MIXING RATIO
2173 ! AS THIS MAY OTHERWISE RESULT IN TOO MUCH TRANSFER OF WATER DURING
2174 ! RIME-SPLINTERING
2176 PRACS(K) = MIN(PRACS(K),QR3D(K)/DT)
2178 ! COLLECTION OF SNOW BY RAIN - NEEDED FOR GRAUPEL CONVERSION CALCULATIONS
2179 ! ONLY CALCULATE IF SNOW AND RAIN MIXING RATIOS EXCEED 0.1 G/KG
2181 ! ASSUME COLLECTION OF SNOW BY RAIN PRODUCES GRAUPEL NOT HAIL
2183 IF (IHAIL.EQ.0) THEN
2184 IF (QNI3D(K).GE.0.1E-3.AND.QR3D(K).GE.0.1E-3) THEN
2185 PSACR(K) = CONS31*(((1.2*UMR-0.95*UMS)**2+ &
2186 0.08*UMS*UMR)**0.5*RHO(K)* &
2187 N0RR(K)*N0S(K)/LAMS(K)**3* &
2188 (5./(LAMS(K)**3*LAMR(K))+ &
2189 2./(LAMS(K)**2*LAMR(K)**2)+ &
2190 0.5/(LAMS(K)*LAMR(K)**3)))
2191 END IF
2192 END IF
2194 END IF
2196 !.......................................................................
2198 ! COLLECTION OF RAINWATER BY GRAUPEL, FROM IKAWA AND SAITO 1990,
2199 ! USED BY REISNER ET AL 1998
2200 IF (QR3D(K).GE.1.E-8.AND.QG3D(K).GE.1.E-8) THEN
2202 UMG = AGN(K)*CONS7/(LAMG(K)**BG)
2203 UMR = ARN(K)*CONS4/(LAMR(K)**BR)
2204 UNG = AGN(K)*CONS8/LAMG(K)**BG
2205 UNR = ARN(K)*CONS6/LAMR(K)**BR
2207 ! SET REASLISTIC LIMITS ON FALLSPEEDS
2208 UMG=MIN(UMG,20.)
2209 UNG=MIN(UNG,20.)
2210 UMR=MIN(UMR,9.1)
2211 UNR=MIN(UNR,9.1)
2213 PRACG(K) = CONS41*(((1.2*UMR-0.95*UMG)**2+ &
2214 0.08*UMG*UMR)**0.5*RHO(K)* &
2215 N0RR(K)*N0G(K)/LAMR(K)**3* &
2216 (5./(LAMR(K)**3*LAMG(K))+ &
2217 2./(LAMR(K)**2*LAMG(K)**2)+ &
2218 0.5/(LAMR(k)*LAMG(k)**3)))
2220 NPRACG(K) = CONS32*RHO(K)*(1.7*(UNR-UNG)**2+ &
2221 0.3*UNR*UNG)**0.5*N0RR(K)*N0G(K)* &
2222 (1./(LAMR(K)**3*LAMG(K))+ &
2223 1./(LAMR(K)**2*LAMG(K)**2)+ &
2224 1./(LAMR(K)*LAMG(K)**3))
2226 ! MAKE SURE PRACG DOESN'T EXCEED TOTAL RAIN MIXING RATIO
2227 ! AS THIS MAY OTHERWISE RESULT IN TOO MUCH TRANSFER OF WATER DURING
2228 ! RIME-SPLINTERING
2230 PRACG(K) = MIN(PRACG(K),QR3D(K)/DT)
2232 END IF
2234 !.......................................................................
2235 ! RIME-SPLINTERING - SNOW
2236 ! HALLET-MOSSOP (1974)
2237 ! NUMBER OF SPLINTERS FORMED IS BASED ON MASS OF RIMED WATER
2239 ! DUM1 = MASS OF INDIVIDUAL SPLINTERS
2241 ! HM ADD THRESHOLD SNOW AND DROPLET MIXING RATIO FOR RIME-SPLINTERING
2242 ! TO LIMIT RIME-SPLINTERING IN STRATIFORM CLOUDS
2243 ! THESE THRESHOLDS CORRESPOND WITH GRAUPEL THRESHOLDS IN RH 1984
2245 IF (QNI3D(K).GE.0.1E-3.AND.QC3D(K).GE.0.5E-3) THEN
2246 IF (PSACWS(K).GT.0..OR.PRACS(K).GT.0.) THEN
2247 IF (T3D(K).LT.270.16 .AND. T3D(K).GT.265.16) THEN
2249 IF (T3D(K).GT.270.16) THEN
2250 FMULT = 0.
2251 ELSE IF (T3D(K).LE.270.16.AND.T3D(K).GT.268.16) THEN
2252 FMULT = (270.16-T3D(K))/2.
2253 ELSE IF (T3D(K).GE.265.16.AND.T3D(K).LE.268.16) THEN
2254 FMULT = (T3D(K)-265.16)/3.
2255 ELSE IF (T3D(K).LT.265.16) THEN
2256 FMULT = 0.
2257 END IF
2259 ! 1000 IS TO CONVERT FROM KG TO G
2261 ! SPLINTERING FROM DROPLETS ACCRETED ONTO SNOW
2263 IF (PSACWS(K).GT.0.) THEN
2264 NMULTS(K) = 35.E4*PSACWS(K)*FMULT*1000.
2265 QMULTS(K) = NMULTS(K)*MMULT
2267 ! CONSTRAIN SO THAT TRANSFER OF MASS FROM SNOW TO ICE CANNOT BE MORE MASS
2268 ! THAN WAS RIMED ONTO SNOW
2270 QMULTS(K) = MIN(QMULTS(K),PSACWS(K))
2271 PSACWS(K) = PSACWS(K)-QMULTS(K)
2273 END IF
2275 ! RIMING AND SPLINTERING FROM ACCRETED RAINDROPS
2277 IF (PRACS(K).GT.0.) THEN
2278 NMULTR(K) = 35.E4*PRACS(K)*FMULT*1000.
2279 QMULTR(K) = NMULTR(K)*MMULT
2281 ! CONSTRAIN SO THAT TRANSFER OF MASS FROM SNOW TO ICE CANNOT BE MORE MASS
2282 ! THAN WAS RIMED ONTO SNOW
2284 QMULTR(K) = MIN(QMULTR(K),PRACS(K))
2286 PRACS(K) = PRACS(K)-QMULTR(K)
2288 END IF
2290 END IF
2291 END IF
2292 END IF
2294 !.......................................................................
2295 ! RIME-SPLINTERING - GRAUPEL
2296 ! HALLET-MOSSOP (1974)
2297 ! NUMBER OF SPLINTERS FORMED IS BASED ON MASS OF RIMED WATER
2299 ! DUM1 = MASS OF INDIVIDUAL SPLINTERS
2301 ! HM ADD THRESHOLD SNOW MIXING RATIO FOR RIME-SPLINTERING
2302 ! TO LIMIT RIME-SPLINTERING IN STRATIFORM CLOUDS
2304 ! ONLY CALCULATE FOR GRAUPEL NOT HAIL
2305 IF (IHAIL.EQ.0) THEN
2306 IF (PSACWG(K).GT.0..OR.PRACG(K).GT.0.) THEN
2307 IF (T3D(K).LT.270.16 .AND. T3D(K).GT.265.16) THEN
2309 IF (T3D(K).GT.270.16) THEN
2310 FMULT = 0.
2311 ELSE IF (T3D(K).LE.270.16.AND.T3D(K).GT.268.16) THEN
2312 FMULT = (270.16-T3D(K))/2.
2313 ELSE IF (T3D(K).GE.265.16.AND.T3D(K).LE.268.16) THEN
2314 FMULT = (T3D(K)-265.16)/3.
2315 ELSE IF (T3D(K).LT.265.16) THEN
2316 FMULT = 0.
2317 END IF
2319 ! 1000 IS TO CONVERT FROM KG TO G
2321 ! SPLINTERING FROM DROPLETS ACCRETED ONTO GRAUPEL
2323 IF (PSACWG(K).GT.0.) THEN
2324 NMULTG(K) = 35.E4*PSACWG(K)*FMULT*1000.
2325 QMULTG(K) = NMULTG(K)*MMULT
2327 ! CONSTRAIN SO THAT TRANSFER OF MASS FROM GRAUPEL TO ICE CANNOT BE MORE MASS
2328 ! THAN WAS RIMED ONTO GRAUPEL
2330 QMULTG(K) = MIN(QMULTG(K),PSACWG(K))
2331 PSACWG(K) = PSACWG(K)-QMULTG(K)
2333 END IF
2335 ! RIMING AND SPLINTERING FROM ACCRETED RAINDROPS
2337 IF (PRACG(K).GT.0.) THEN
2338 NMULTRG(K) = 35.E4*PRACG(K)*FMULT*1000.
2339 QMULTRG(K) = NMULTRG(K)*MMULT
2341 ! CONSTRAIN SO THAT TRANSFER OF MASS FROM GRAUPEL TO ICE CANNOT BE MORE MASS
2342 ! THAN WAS RIMED ONTO GRAUPEL
2344 QMULTRG(K) = MIN(QMULTRG(K),PRACG(K))
2345 PRACG(K) = PRACG(K)-QMULTRG(K)
2347 END IF
2349 END IF
2350 END IF
2351 END IF
2353 !........................................................................
2354 ! CONVERSION OF RIMED CLOUD WATER ONTO SNOW TO GRAUPEL
2355 ! ASSUME CONVERTED SNOW FORMS GRAUPEL NOT HAIL
2356 ! HAIL ASSUMED TO ONLY FORM BY FREEZING OF RAIN
2357 ! OR COLLISIONS OF RAIN WITH CLOUD ICE
2359 IF (IHAIL.EQ.0) THEN
2360 IF (PSACWS(K).GT.0.) THEN
2361 ! ONLY ALLOW CONVERSION IF QNI > 0.1 AND QC > 0.5 G/KG FOLLOWING RUTLEDGE AND HOBBS (1984)
2362 IF (QNI3D(K).GE.0.1E-3.AND.QC3D(K).GE.0.5E-3) THEN
2364 ! PORTION OF RIMING CONVERTED TO GRAUPEL (REISNER ET AL. 1998, ORIGINALLY IS1991)
2365 PGSACW(K) = MIN(PSACWS(K),CONS17*DT*N0S(K)*QC3D(K)*QC3D(K)* &
2366 ASN(K)*ASN(K)/ &
2367 (RHO(K)*LAMS(K)**(2.*BS+2.)))
2369 ! MIX RAT CONVERTED INTO GRAUPEL AS EMBRYO (REISNER ET AL. 1998, ORIG M1990)
2370 DUM = MAX(RHOSN/(RHOG-RHOSN)*PGSACW(K),0.)
2372 ! NUMBER CONCENTRAITON OF EMBRYO GRAUPEL FROM RIMING OF SNOW
2373 NSCNG(K) = DUM/MG0*RHO(K)
2374 ! LIMIT MAX NUMBER CONVERTED TO SNOW NUMBER
2375 NSCNG(K) = MIN(NSCNG(K),NS3D(K)/DT)
2377 ! PORTION OF RIMING LEFT FOR SNOW
2378 PSACWS(K) = PSACWS(K) - PGSACW(K)
2379 END IF
2380 END IF
2382 ! CONVERSION OF RIMED RAINWATER ONTO SNOW CONVERTED TO GRAUPEL
2384 IF (PRACS(K).GT.0.) THEN
2385 ! ONLY ALLOW CONVERSION IF QNI > 0.1 AND QR > 0.1 G/KG FOLLOWING RUTLEDGE AND HOBBS (1984)
2386 IF (QNI3D(K).GE.0.1E-3.AND.QR3D(K).GE.0.1E-3) THEN
2387 ! PORTION OF COLLECTED RAINWATER CONVERTED TO GRAUPEL (REISNER ET AL. 1998)
2388 DUM = CONS18*(4./LAMS(K))**3*(4./LAMS(K))**3 &
2389 /(CONS18*(4./LAMS(K))**3*(4./LAMS(K))**3+ &
2390 CONS19*(4./LAMR(K))**3*(4./LAMR(K))**3)
2391 DUM=MIN(DUM,1.)
2392 DUM=MAX(DUM,0.)
2393 PGRACS(K) = (1.-DUM)*PRACS(K)
2394 NGRACS(K) = (1.-DUM)*NPRACS(K)
2395 ! LIMIT MAX NUMBER CONVERTED TO MIN OF EITHER RAIN OR SNOW NUMBER CONCENTRATION
2396 NGRACS(K) = MIN(NGRACS(K),NR3D(K)/DT)
2397 NGRACS(K) = MIN(NGRACS(K),NS3D(K)/DT)
2399 ! AMOUNT LEFT FOR SNOW PRODUCTION
2400 PRACS(K) = PRACS(K) - PGRACS(K)
2401 NPRACS(K) = NPRACS(K) - NGRACS(K)
2402 ! CONVERSION TO GRAUPEL DUE TO COLLECTION OF SNOW BY RAIN
2403 PSACR(K)=PSACR(K)*(1.-DUM)
2404 END IF
2405 END IF
2406 END IF
2408 !.......................................................................
2409 ! FREEZING OF RAIN DROPS
2410 ! FREEZING ALLOWED BELOW -4 C
2412 IF (T3D(K).LT.269.15.AND.QR3D(K).GE.QSMALL) THEN
2414 ! IMMERSION FREEZING (BIGG 1953)
2415 MNUCCR(K) = CONS20*NR3D(K)*EXP(AIMM*(273.15-T3D(K)))/LAMR(K)**3 &
2416 /LAMR(K)**3
2418 NNUCCR(K) = PI*NR3D(K)*BIMM*EXP(AIMM*(273.15-T3D(K)))/LAMR(K)**3
2420 ! PREVENT DIVERGENCE BETWEEN MIXING RATIO AND NUMBER CONC
2421 NNUCCR(K) = MIN(NNUCCR(K),NR3D(K)/DT)
2423 END IF
2425 !.......................................................................
2426 ! ACCRETION OF CLOUD LIQUID WATER BY RAIN
2427 ! CONTINUOUS COLLECTION EQUATION WITH
2428 ! GRAVITATIONAL COLLECTION KERNEL, DROPLET FALL SPEED NEGLECTED
2430 IF (QR3D(K).GE.1.E-8 .AND. QC3D(K).GE.1.E-8) THEN
2432 ! 12/13/06 HM ADD, REPLACE WITH NEWER FORMULA FROM
2433 ! KHAIROUTDINOV AND KOGAN 2000, MWR
2435 DUM=(QC3D(K)*QR3D(K))
2436 PRA(K) = 67.*(DUM)**1.15
2437 NPRA(K) = PRA(K)/(QC3D(K)/NC3D(K))
2439 END IF
2440 !.......................................................................
2441 ! SELF-COLLECTION OF RAIN DROPS
2442 ! FROM BEHENG(1994)
2443 ! FROM NUMERICAL SIMULATION OF THE STOCHASTIC COLLECTION EQUATION
2444 ! AS DESCRINED ABOVE FOR AUTOCONVERSION
2446 IF (QR3D(K).GE.1.E-8) THEN
2447 NRAGG(K) = -8.*NR3D(K)*QR3D(K)*RHO(K)
2448 END IF
2450 !.......................................................................
2451 ! AUTOCONVERSION OF CLOUD ICE TO SNOW
2452 ! FOLLOWING HARRINGTON ET AL. (1995) WITH MODIFICATION
2453 ! HERE IT IS ASSUMED THAT AUTOCONVERSION CAN ONLY OCCUR WHEN THE
2454 ! ICE IS GROWING, I.E. IN CONDITIONS OF ICE SUPERSATURATION
2456 IF (QI3D(K).GE.1.E-8 .AND.QVQVSI(K).GE.1.) THEN
2458 ! COFFI = 2./LAMI(K)
2459 ! IF (COFFI.GE.DCS) THEN
2460 NPRCI(K) = CONS21*(QV3D(K)-QVI(K))*RHO(K) &
2461 *N0I(K)*EXP(-LAMI(K)*DCS)*DV(K)/ABI(K)
2462 PRCI(K) = CONS22*NPRCI(K)
2463 NPRCI(K) = MIN(NPRCI(K),NI3D(K)/DT)
2465 ! END IF
2466 END IF
2468 !.......................................................................
2469 ! ACCRETION OF CLOUD ICE BY SNOW
2470 ! FOR THIS CALCULATION, IT IS ASSUMED THAT THE VS >> VI
2471 ! AND DS >> DI FOR CONTINUOUS COLLECTION
2473 IF (QNI3D(K).GE.1.E-8 .AND. QI3D(K).GE.QSMALL) THEN
2474 PRAI(K) = CONS23*ASN(K)*QI3D(K)*RHO(K)*N0S(K)/ &
2475 LAMS(K)**(BS+3.)
2476 NPRAI(K) = CONS23*ASN(K)*NI3D(K)* &
2477 RHO(K)*N0S(K)/ &
2478 LAMS(K)**(BS+3.)
2479 NPRAI(K)=MIN(NPRAI(K),NI3D(K)/DT)
2480 END IF
2482 !.......................................................................
2483 ! HM, ADD 12/13/06, COLLISION OF RAIN AND ICE TO PRODUCE SNOW OR GRAUPEL
2484 ! FOLLOWS REISNER ET AL. 1998
2485 ! ASSUMED FALLSPEED AND SIZE OF ICE CRYSTAL << THAN FOR RAIN
2487 IF (QR3D(K).GE.1.E-8.AND.QI3D(K).GE.1.E-8.AND.T3D(K).LE.273.15) THEN
2489 ! ALLOW GRAUPEL FORMATION FROM RAIN-ICE COLLISIONS ONLY IF RAIN MIXING RATIO > 0.1 G/KG,
2490 ! OTHERWISE ADD TO SNOW
2492 IF (QR3D(K).GE.0.1E-3) THEN
2493 NIACR(K)=CONS24*NI3D(K)*N0RR(K)*ARN(K) &
2494 /LAMR(K)**(BR+3.)*RHO(K)
2495 PIACR(K)=CONS25*NI3D(K)*N0RR(K)*ARN(K) &
2496 /LAMR(K)**(BR+3.)/LAMR(K)**3*RHO(K)
2497 PRACI(K)=CONS24*QI3D(K)*N0RR(K)*ARN(K)/ &
2498 LAMR(K)**(BR+3.)*RHO(K)
2499 NIACR(K)=MIN(NIACR(K),NR3D(K)/DT)
2500 NIACR(K)=MIN(NIACR(K),NI3D(K)/DT)
2501 ELSE
2502 NIACRS(K)=CONS24*NI3D(K)*N0RR(K)*ARN(K) &
2503 /LAMR(K)**(BR+3.)*RHO(K)
2504 PIACRS(K)=CONS25*NI3D(K)*N0RR(K)*ARN(K) &
2505 /LAMR(K)**(BR+3.)/LAMR(K)**3*RHO(K)
2506 PRACIS(K)=CONS24*QI3D(K)*N0RR(K)*ARN(K)/ &
2507 LAMR(K)**(BR+3.)*RHO(K)
2508 NIACRS(K)=MIN(NIACRS(K),NR3D(K)/DT)
2509 NIACRS(K)=MIN(NIACRS(K),NI3D(K)/DT)
2510 END IF
2511 END IF
2513 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
2514 ! NUCLEATION OF CLOUD ICE FROM HOMOGENEOUS AND HETEROGENEOUS FREEZING ON AEROSOL
2516 IF (INUC.EQ.0) THEN
2518 ! FREEZING OF AEROSOL ONLY ALLOWED BELOW -5 C
2519 ! AND ABOVE DELIQUESCENCE THRESHOLD OF 80%
2520 ! AND ABOVE ICE SATURATION
2522 ! add threshold according to Greg Thomspon
2524 if ((QVQVS(K).GE.0.999.and.T3D(K).le.265.15).or. &
2525 QVQVSI(K).ge.1.08) then
2527 ! hm, modify dec. 5, 2006, replace with cooper curve
2528 kc2 = 0.005*exp(0.304*(273.15-T3D(K)))*1000. ! convert from L-1 to m-3
2529 ! limit to 500 L-1
2530 kc2 = min(kc2,500.e3)
2531 kc2=MAX(kc2/rho(k),0.) ! convert to kg-1
2533 IF (KC2.GT.NI3D(K)+NS3D(K)+NG3D(K)) THEN
2534 NNUCCD(K) = (KC2-NI3D(K)-NS3D(K)-NG3D(K))/DT
2535 MNUCCD(K) = NNUCCD(K)*MI0
2536 END IF
2538 END IF
2540 ELSE IF (INUC.EQ.1) THEN
2542 IF (T3D(K).LT.273.15.AND.QVQVSI(K).GT.1.) THEN
2544 KC2 = 0.16*1000./RHO(K) ! CONVERT FROM L-1 TO KG-1
2545 IF (KC2.GT.NI3D(K)+NS3D(K)+NG3D(K)) THEN
2546 NNUCCD(K) = (KC2-NI3D(K)-NS3D(K)-NG3D(K))/DT
2547 MNUCCD(K) = NNUCCD(K)*MI0
2548 END IF
2549 END IF
2551 END IF
2553 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
2555 101 CONTINUE
2557 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
2558 ! CALCULATE EVAP/SUB/DEP TERMS FOR QI,QNI,QR
2560 ! NO VENTILATION FOR CLOUD ICE
2562 IF (QI3D(K).GE.QSMALL) THEN
2564 EPSI = 2.*PI*N0I(K)*RHO(K)*DV(K)/(LAMI(K)*LAMI(K))
2566 ELSE
2567 EPSI = 0.
2568 END IF
2570 IF (QNI3D(K).GE.QSMALL) THEN
2571 EPSS = 2.*PI*N0S(K)*RHO(K)*DV(K)* &
2572 (F1S/(LAMS(K)*LAMS(K))+ &
2573 F2S*(ASN(K)*RHO(K)/MU(K))**0.5* &
2574 SC(K)**(1./3.)*CONS10/ &
2575 (LAMS(K)**CONS35))
2576 ELSE
2577 EPSS = 0.
2578 END IF
2580 IF (QG3D(K).GE.QSMALL) THEN
2581 EPSG = 2.*PI*N0G(K)*RHO(K)*DV(K)* &
2582 (F1S/(LAMG(K)*LAMG(K))+ &
2583 F2S*(AGN(K)*RHO(K)/MU(K))**0.5* &
2584 SC(K)**(1./3.)*CONS11/ &
2585 (LAMG(K)**CONS36))
2588 ELSE
2589 EPSG = 0.
2590 END IF
2592 IF (QR3D(K).GE.QSMALL) THEN
2593 EPSR = 2.*PI*N0RR(K)*RHO(K)*DV(K)* &
2594 (F1R/(LAMR(K)*LAMR(K))+ &
2595 F2R*(ARN(K)*RHO(K)/MU(K))**0.5* &
2596 SC(K)**(1./3.)*CONS9/ &
2597 (LAMR(K)**CONS34))
2598 ELSE
2599 EPSR = 0.
2600 END IF
2602 ! ONLY INCLUDE REGION OF ICE SIZE DIST < DCS
2603 ! DUM IS FRACTION OF D*N(D) < DCS
2605 ! LOGIC BELOW FOLLOWS THAT OF HARRINGTON ET AL. 1995 (JAS)
2606 IF (QI3D(K).GE.QSMALL) THEN
2607 DUM=(1.-EXP(-LAMI(K)*DCS)*(1.+LAMI(K)*DCS))
2608 PRD(K) = EPSI*(QV3D(K)-QVI(K))/ABI(K)*DUM
2609 ELSE
2610 DUM=0.
2611 END IF
2612 ! ADD DEPOSITION IN TAIL OF ICE SIZE DIST TO SNOW IF SNOW IS PRESENT
2613 IF (QNI3D(K).GE.QSMALL) THEN
2614 PRDS(K) = EPSS*(QV3D(K)-QVI(K))/ABI(K)+ &
2615 EPSI*(QV3D(K)-QVI(K))/ABI(K)*(1.-DUM)
2616 ! OTHERWISE ADD TO CLOUD ICE
2617 ELSE
2618 PRD(K) = PRD(K)+EPSI*(QV3D(K)-QVI(K))/ABI(K)*(1.-DUM)
2619 END IF
2620 ! VAPOR DPEOSITION ON GRAUPEL
2621 PRDG(K) = EPSG*(QV3D(K)-QVI(K))/ABI(K)
2623 ! NO CONDENSATION ONTO RAIN, ONLY EVAP
2625 IF (QV3D(K).LT.QVS(K)) THEN
2626 PRE(K) = EPSR*(QV3D(K)-QVS(K))/AB(K)
2627 PRE(K) = MIN(PRE(K),0.)
2628 ELSE
2629 PRE(K) = 0.
2630 END IF
2632 ! MAKE SURE NOT PUSHED INTO ICE SUPERSAT/SUBSAT
2633 ! FORMULA FROM REISNER 2 SCHEME
2635 DUM = (QV3D(K)-QVI(K))/DT
2637 FUDGEF = 0.9999
2638 SUM_DEP = PRD(K)+PRDS(K)+MNUCCD(K)+PRDG(K)
2640 IF( (DUM.GT.0. .AND. SUM_DEP.GT.DUM*FUDGEF) .OR. &
2641 (DUM.LT.0. .AND. SUM_DEP.LT.DUM*FUDGEF) ) THEN
2642 MNUCCD(K) = FUDGEF*MNUCCD(K)*DUM/SUM_DEP
2643 PRD(K) = FUDGEF*PRD(K)*DUM/SUM_DEP
2644 PRDS(K) = FUDGEF*PRDS(K)*DUM/SUM_DEP
2645 PRDG(K) = FUDGEF*PRDG(K)*DUM/SUM_DEP
2646 ENDIF
2648 ! IF CLOUD ICE/SNOW/GRAUPEL VAP DEPOSITION IS NEG, THEN ASSIGN TO SUBLIMATION PROCESSES
2650 IF (PRD(K).LT.0.) THEN
2651 EPRD(K)=PRD(K)
2652 PRD(K)=0.
2653 END IF
2654 IF (PRDS(K).LT.0.) THEN
2655 EPRDS(K)=PRDS(K)
2656 PRDS(K)=0.
2657 END IF
2658 IF (PRDG(K).LT.0.) THEN
2659 EPRDG(K)=PRDG(K)
2660 PRDG(K)=0.
2661 END IF
2663 !.......................................................................
2664 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
2666 ! CONSERVATION OF WATER
2667 ! THIS IS ADOPTED LOOSELY FROM MM5 RESINER CODE. HOWEVER, HERE WE
2668 ! ONLY ADJUST PROCESSES THAT ARE NEGATIVE, RATHER THAN ALL PROCESSES.
2670 ! IF MIXING RATIOS LESS THAN QSMALL, THEN NO DEPLETION OF WATER
2671 ! THROUGH MICROPHYSICAL PROCESSES, SKIP CONSERVATION
2673 ! NOTE: CONSERVATION CHECK NOT APPLIED TO NUMBER CONCENTRATION SPECIES. ADDITIONAL CATCH
2674 ! BELOW WILL PREVENT NEGATIVE NUMBER CONCENTRATION
2675 ! FOR EACH MICROPHYSICAL PROCESS WHICH PROVIDES A SOURCE FOR NUMBER, THERE IS A CHECK
2676 ! TO MAKE SURE THAT CAN'T EXCEED TOTAL NUMBER OF DEPLETED SPECIES WITH THE TIME
2677 ! STEP
2679 !****SENSITIVITY - NO ICE
2681 IF (ILIQ.EQ.1) THEN
2682 MNUCCC(K)=0.
2683 NNUCCC(K)=0.
2684 MNUCCR(K)=0.
2685 NNUCCR(K)=0.
2686 MNUCCD(K)=0.
2687 NNUCCD(K)=0.
2688 END IF
2690 ! ****SENSITIVITY - NO GRAUPEL
2691 IF (IGRAUP.EQ.1) THEN
2692 PRACG(K) = 0.
2693 PSACR(K) = 0.
2694 PSACWG(K) = 0.
2695 PGSACW(K) = 0.
2696 PGRACS(K) = 0.
2697 PRDG(K) = 0.
2698 EPRDG(K) = 0.
2699 EVPMG(K) = 0.
2700 PGMLT(K) = 0.
2701 NPRACG(K) = 0.
2702 NPSACWG(K) = 0.
2703 NSCNG(K) = 0.
2704 NGRACS(K) = 0.
2705 NSUBG(K) = 0.
2706 NGMLTG(K) = 0.
2707 NGMLTR(K) = 0.
2708 END IF
2710 ! CONSERVATION OF QC
2712 DUM = (PRC(K)+PRA(K)+MNUCCC(K)+PSACWS(K)+PSACWI(K)+QMULTS(K)+PSACWG(K)+PGSACW(K)+QMULTG(K))*DT
2714 IF (DUM.GT.QC3D(K).AND.QC3D(K).GE.QSMALL) THEN
2715 RATIO = QC3D(K)/DUM
2717 PRC(K) = PRC(K)*RATIO
2718 PRA(K) = PRA(K)*RATIO
2719 MNUCCC(K) = MNUCCC(K)*RATIO
2720 PSACWS(K) = PSACWS(K)*RATIO
2721 PSACWI(K) = PSACWI(K)*RATIO
2722 QMULTS(K) = QMULTS(K)*RATIO
2723 QMULTG(K) = QMULTG(K)*RATIO
2724 PSACWG(K) = PSACWG(K)*RATIO
2725 PGSACW(K) = PGSACW(K)*RATIO
2726 END IF
2728 ! CONSERVATION OF QI
2730 DUM = (-PRD(K)-MNUCCC(K)+PRCI(K)+PRAI(K)-QMULTS(K)-QMULTG(K)-QMULTR(K)-QMULTRG(K) &
2731 -MNUCCD(K)+PRACI(K)+PRACIS(K)-EPRD(K)-PSACWI(K))*DT
2733 IF (DUM.GT.QI3D(K).AND.QI3D(K).GE.QSMALL) THEN
2735 RATIO = (QI3D(K)/DT+PRD(K)+MNUCCC(K)+QMULTS(K)+QMULTG(K)+QMULTR(K)+QMULTRG(K)+ &
2736 MNUCCD(K)+PSACWI(K))/ &
2737 (PRCI(K)+PRAI(K)+PRACI(K)+PRACIS(K)-EPRD(K))
2739 PRCI(K) = PRCI(K)*RATIO
2740 PRAI(K) = PRAI(K)*RATIO
2741 PRACI(K) = PRACI(K)*RATIO
2742 PRACIS(K) = PRACIS(K)*RATIO
2743 EPRD(K) = EPRD(K)*RATIO
2745 END IF
2747 ! CONSERVATION OF QR
2749 DUM=((PRACS(K)-PRE(K))+(QMULTR(K)+QMULTRG(K)-PRC(K))+(MNUCCR(K)-PRA(K))+ &
2750 PIACR(K)+PIACRS(K)+PGRACS(K)+PRACG(K))*DT
2752 IF (DUM.GT.QR3D(K).AND.QR3D(K).GE.QSMALL) THEN
2754 RATIO = (QR3D(K)/DT+PRC(K)+PRA(K))/ &
2755 (-PRE(K)+QMULTR(K)+QMULTRG(K)+PRACS(K)+MNUCCR(K)+PIACR(K)+PIACRS(K)+PGRACS(K)+PRACG(K))
2757 PRE(K) = PRE(K)*RATIO
2758 PRACS(K) = PRACS(K)*RATIO
2759 QMULTR(K) = QMULTR(K)*RATIO
2760 QMULTRG(K) = QMULTRG(K)*RATIO
2761 MNUCCR(K) = MNUCCR(K)*RATIO
2762 PIACR(K) = PIACR(K)*RATIO
2763 PIACRS(K) = PIACRS(K)*RATIO
2764 PGRACS(K) = PGRACS(K)*RATIO
2765 PRACG(K) = PRACG(K)*RATIO
2767 END IF
2769 ! CONSERVATION OF QNI
2770 ! CONSERVATION FOR GRAUPEL SCHEME
2772 IF (IGRAUP.EQ.0) THEN
2774 DUM = (-PRDS(K)-PSACWS(K)-PRAI(K)-PRCI(K)-PRACS(K)-EPRDS(K)+PSACR(K)-PIACRS(K)-PRACIS(K))*DT
2776 IF (DUM.GT.QNI3D(K).AND.QNI3D(K).GE.QSMALL) THEN
2778 RATIO = (QNI3D(K)/DT+PRDS(K)+PSACWS(K)+PRAI(K)+PRCI(K)+PRACS(K)+PIACRS(K)+PRACIS(K))/(-EPRDS(K)+PSACR(K))
2780 EPRDS(K) = EPRDS(K)*RATIO
2781 PSACR(K) = PSACR(K)*RATIO
2783 END IF
2785 ! FOR NO GRAUPEL, NEED TO INCLUDE FREEZING OF RAIN FOR SNOW
2786 ELSE IF (IGRAUP.EQ.1) THEN
2788 DUM = (-PRDS(K)-PSACWS(K)-PRAI(K)-PRCI(K)-PRACS(K)-EPRDS(K)+PSACR(K)-PIACRS(K)-PRACIS(K)-MNUCCR(K))*DT
2790 IF (DUM.GT.QNI3D(K).AND.QNI3D(K).GE.QSMALL) THEN
2792 RATIO = (QNI3D(K)/DT+PRDS(K)+PSACWS(K)+PRAI(K)+PRCI(K)+PRACS(K)+PIACRS(K)+PRACIS(K)+MNUCCR(K))/(-EPRDS(K)+PSACR(K))
2794 EPRDS(K) = EPRDS(K)*RATIO
2795 PSACR(K) = PSACR(K)*RATIO
2797 END IF
2799 END IF
2801 ! CONSERVATION OF QG
2803 DUM = (-PSACWG(K)-PRACG(K)-PGSACW(K)-PGRACS(K)-PRDG(K)-MNUCCR(K)-EPRDG(K)-PIACR(K)-PRACI(K)-PSACR(K))*DT
2805 IF (DUM.GT.QG3D(K).AND.QG3D(K).GE.QSMALL) THEN
2807 RATIO = (QG3D(K)/DT+PSACWG(K)+PRACG(K)+PGSACW(K)+PGRACS(K)+PRDG(K)+MNUCCR(K)+PSACR(K)+&
2808 PIACR(K)+PRACI(K))/(-EPRDG(K))
2810 EPRDG(K) = EPRDG(K)*RATIO
2812 END IF
2814 ! TENDENCIES
2816 QV3DTEN(K) = QV3DTEN(K)+(-PRE(K)-PRD(K)-PRDS(K)-MNUCCD(K)-EPRD(K)-EPRDS(K)-PRDG(K)-EPRDG(K))
2817 T3DTEN(K) = T3DTEN(K)+(PRE(K) &
2818 *XXLV(K)+(PRD(K)+PRDS(K)+ &
2819 MNUCCD(K)+EPRD(K)+EPRDS(K)+PRDG(K)+EPRDG(K))*XXLS(K)+ &
2820 (PSACWS(K)+PSACWI(K)+MNUCCC(K)+MNUCCR(K)+ &
2821 QMULTS(K)+QMULTG(K)+QMULTR(K)+QMULTRG(K)+PRACS(K) &
2822 +PSACWG(K)+PRACG(K)+PGSACW(K)+PGRACS(K))*XLF(K))/CPM(K)
2824 QC3DTEN(K) = QC3DTEN(K)+ &
2825 (-PRA(K)-PRC(K)-MNUCCC(K)+PCC(K)- &
2826 PSACWS(K)-PSACWI(K)-QMULTS(K)-QMULTG(K)-PSACWG(K)-PGSACW(K))
2827 QI3DTEN(K) = QI3DTEN(K)+ &
2828 (PRD(K)+EPRD(K)+PSACWI(K)+MNUCCC(K)-PRCI(K)- &
2829 PRAI(K)+QMULTS(K)+QMULTG(K)+QMULTR(K)+QMULTRG(K)+MNUCCD(K)-PRACI(K)-PRACIS(K))
2830 QR3DTEN(K) = QR3DTEN(K)+ &
2831 (PRE(K)+PRA(K)+PRC(K)-PRACS(K)-MNUCCR(K)-QMULTR(K)-QMULTRG(K) &
2832 -PIACR(K)-PIACRS(K)-PRACG(K)-PGRACS(K))
2834 IF (IGRAUP.EQ.0) THEN
2836 QNI3DTEN(K) = QNI3DTEN(K)+ &
2837 (PRAI(K)+PSACWS(K)+PRDS(K)+PRACS(K)+PRCI(K)+EPRDS(K)-PSACR(K)+PIACRS(K)+PRACIS(K))
2838 NS3DTEN(K) = NS3DTEN(K)+(NSAGG(K)+NPRCI(K)-NSCNG(K)-NGRACS(K)+NIACRS(K))
2839 QG3DTEN(K) = QG3DTEN(K)+(PRACG(K)+PSACWG(K)+PGSACW(K)+PGRACS(K)+ &
2840 PRDG(K)+EPRDG(K)+MNUCCR(K)+PIACR(K)+PRACI(K)+PSACR(K))
2841 NG3DTEN(K) = NG3DTEN(K)+(NSCNG(K)+NGRACS(K)+NNUCCR(K)+NIACR(K))
2843 ! FOR NO GRAUPEL, NEED TO INCLUDE FREEZING OF RAIN FOR SNOW
2844 ELSE IF (IGRAUP.EQ.1) THEN
2846 QNI3DTEN(K) = QNI3DTEN(K)+ &
2847 (PRAI(K)+PSACWS(K)+PRDS(K)+PRACS(K)+PRCI(K)+EPRDS(K)-PSACR(K)+PIACRS(K)+PRACIS(K)+MNUCCR(K))
2848 NS3DTEN(K) = NS3DTEN(K)+(NSAGG(K)+NPRCI(K)-NSCNG(K)-NGRACS(K)+NIACRS(K)+NNUCCR(K))
2850 END IF
2852 NC3DTEN(K) = NC3DTEN(K)+(-NNUCCC(K)-NPSACWS(K) &
2853 -NPRA(K)-NPRC(K)-NPSACWI(K)-NPSACWG(K))
2855 NI3DTEN(K) = NI3DTEN(K)+ &
2856 (NNUCCC(K)-NPRCI(K)-NPRAI(K)+NMULTS(K)+NMULTG(K)+NMULTR(K)+NMULTRG(K)+ &
2857 NNUCCD(K)-NIACR(K)-NIACRS(K))
2859 NR3DTEN(K) = NR3DTEN(K)+(NPRC1(K)-NPRACS(K)-NNUCCR(K) &
2860 +NRAGG(K)-NIACR(K)-NIACRS(K)-NPRACG(K)-NGRACS(K))
2862 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
2863 ! NOW CALCULATE SATURATION ADJUSTMENT TO CONDENSE EXTRA VAPOR ABOVE
2864 ! WATER SATURATION
2866 DUMT = T3D(K)+DT*T3DTEN(K)
2867 DUMQV = QV3D(K)+DT*QV3DTEN(K)
2868 DUMQSS = 0.622*POLYSVP(DUMT,0)/ (PRES(K)-POLYSVP(DUMT,0))
2869 DUMQC = QC3D(K)+DT*QC3DTEN(K)
2870 DUMQC = MAX(DUMQC,0.)
2872 ! SATURATION ADJUSTMENT FOR LIQUID
2874 DUMS = DUMQV-DUMQSS
2875 PCC(K) = DUMS/(1.+XXLV(K)**2*DUMQSS/(CPM(K)*RV*DUMT**2))/DT
2876 IF (PCC(K)*DT+DUMQC.LT.0.) THEN
2877 PCC(K) = -DUMQC/DT
2878 END IF
2880 QV3DTEN(K) = QV3DTEN(K)-PCC(K)
2881 T3DTEN(K) = T3DTEN(K)+PCC(K)*XXLV(K)/CPM(K)
2882 QC3DTEN(K) = QC3DTEN(K)+PCC(K)
2884 !.......................................................................
2885 ! ACTIVATION OF CLOUD DROPLETS
2886 ! ACTIVATION OF DROPLET CURRENTLY NOT CALCULATED
2887 ! DROPLET CONCENTRATION IS SPECIFIED !!!!!
2889 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
2890 ! SUBLIMATE, MELT, OR EVAPORATE NUMBER CONCENTRATION
2891 ! THIS FORMULATION ASSUMES 1:1 RATIO BETWEEN MASS LOSS AND
2892 ! LOSS OF NUMBER CONCENTRATION
2894 ! IF (PCC(K).LT.0.) THEN
2895 ! DUM = PCC(K)*DT/QC3D(K)
2896 ! DUM = MAX(-1.,DUM)
2897 ! NSUBC(K) = DUM*NC3D(K)/DT
2898 ! END IF
2900 IF (EPRD(K).LT.0.) THEN
2901 DUM = EPRD(K)*DT/QI3D(K)
2902 DUM = MAX(-1.,DUM)
2903 NSUBI(K) = DUM*NI3D(K)/DT
2904 END IF
2905 IF (EPRDS(K).LT.0.) THEN
2906 DUM = EPRDS(K)*DT/QNI3D(K)
2907 DUM = MAX(-1.,DUM)
2908 NSUBS(K) = DUM*NS3D(K)/DT
2909 END IF
2910 IF (PRE(K).LT.0.) THEN
2911 DUM = PRE(K)*DT/QR3D(K)
2912 DUM = MAX(-1.,DUM)
2913 NSUBR(K) = DUM*NR3D(K)/DT
2914 END IF
2915 IF (EPRDG(K).LT.0.) THEN
2916 DUM = EPRDG(K)*DT/QG3D(K)
2917 DUM = MAX(-1.,DUM)
2918 NSUBG(K) = DUM*NG3D(K)/DT
2919 END IF
2921 ! nsubr(k)=0.
2922 ! nsubs(k)=0.
2923 ! nsubg(k)=0.
2925 ! UPDATE TENDENCIES
2927 ! NC3DTEN(K) = NC3DTEN(K)+NSUBC(K)
2928 NI3DTEN(K) = NI3DTEN(K)+NSUBI(K)
2929 NS3DTEN(K) = NS3DTEN(K)+NSUBS(K)
2930 NG3DTEN(K) = NG3DTEN(K)+NSUBG(K)
2931 NR3DTEN(K) = NR3DTEN(K)+NSUBR(K)
2933 END IF !!!!!! TEMPERATURE
2935 ! SWITCH LTRUE TO 1, SINCE HYDROMETEORS ARE PRESENT
2936 LTRUE = 1
2938 200 CONTINUE
2940 END DO
2942 ! INITIALIZE PRECIP AND SNOW RATES
2943 PRECRT = 0.
2944 SNOWRT = 0.
2946 ! IF THERE ARE NO HYDROMETEORS, THEN SKIP TO END OF SUBROUTINE
2948 IF (LTRUE.EQ.0) GOTO 400
2950 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
2951 !.......................................................................
2952 ! CALCULATE SEDIMENATION
2953 ! THE NUMERICS HERE FOLLOW FROM REISNER ET AL. (1998)
2954 ! FALLOUT TERMS ARE CALCULATED ON SPLIT TIME STEPS TO ENSURE NUMERICAL
2955 ! STABILITY, I.E. COURANT# < 1
2957 !.......................................................................
2959 NSTEP = 1
2961 DO K = KTS,KTE
2963 DUMI(K) = QI3D(K)+QI3DTEN(K)*DT
2964 DUMQS(K) = QNI3D(K)+QNI3DTEN(K)*DT
2965 DUMR(K) = QR3D(K)+QR3DTEN(K)*DT
2966 DUMFNI(K) = NI3D(K)+NI3DTEN(K)*DT
2967 DUMFNS(K) = NS3D(K)+NS3DTEN(K)*DT
2968 DUMFNR(K) = NR3D(K)+NR3DTEN(K)*DT
2969 DUMC(K) = QC3D(K)+QC3DTEN(K)*DT
2970 DUMFNC(K) = NC3D(K)+NC3DTEN(K)*DT
2971 DUMG(K) = QG3D(K)+QG3DTEN(K)*DT
2972 DUMFNG(K) = NG3D(K)+NG3DTEN(K)*DT
2974 ! SWITCH FOR CONSTANT DROPLET NUMBER
2975 ! IF (INUM.EQ.1) THEN
2976 DUMFNC(K) = NC3D(K)
2977 ! END IF
2979 ! GET DUMMY LAMDA FOR SEDIMENTATION CALCULATIONS
2981 ! MAKE SURE NUMBER CONCENTRATIONS ARE POSITIVE
2982 DUMFNI(K) = MAX(0.,DUMFNI(K))
2983 DUMFNS(K) = MAX(0.,DUMFNS(K))
2984 DUMFNC(K) = MAX(0.,DUMFNC(K))
2985 DUMFNR(K) = MAX(0.,DUMFNR(K))
2986 DUMFNG(K) = MAX(0.,DUMFNG(K))
2988 !......................................................................
2989 ! CLOUD ICE
2991 IF (DUMI(K).GE.QSMALL) THEN
2992 DLAMI = (CONS12*DUMFNI(K)/DUMI(K))**(1./DI)
2993 DLAMI=MAX(DLAMI,LAMMINI)
2994 DLAMI=MIN(DLAMI,LAMMAXI)
2995 END IF
2996 !......................................................................
2997 ! RAIN
2999 IF (DUMR(K).GE.QSMALL) THEN
3000 DLAMR = (PI*RHOW*DUMFNR(K)/DUMR(K))**(1./3.)
3001 DLAMR=MAX(DLAMR,LAMMINR)
3002 DLAMR=MIN(DLAMR,LAMMAXR)
3003 END IF
3004 !......................................................................
3005 ! CLOUD DROPLETS
3007 IF (DUMC(K).GE.QSMALL) THEN
3008 DUM = PRES(K)/(287.15*T3D(K))
3009 PGAM(K)=0.0005714*(NC3D(K)/1.E6/DUM)+0.2714
3010 PGAM(K)=1./(PGAM(K)**2)-1.
3011 PGAM(K)=MAX(PGAM(K),2.)
3012 PGAM(K)=MIN(PGAM(K),10.)
3014 DLAMC = (CONS26*DUMFNC(K)*GAMMA(PGAM(K)+4.)/(DUMC(K)*GAMMA(PGAM(K)+1.)))**(1./3.)
3015 LAMMIN = (PGAM(K)+1.)/60.E-6
3016 LAMMAX = (PGAM(K)+1.)/1.E-6
3017 DLAMC=MAX(DLAMC,LAMMIN)
3018 DLAMC=MIN(DLAMC,LAMMAX)
3019 END IF
3020 !......................................................................
3021 ! SNOW
3023 IF (DUMQS(K).GE.QSMALL) THEN
3024 DLAMS = (CONS1*DUMFNS(K)/ DUMQS(K))**(1./DS)
3025 DLAMS=MAX(DLAMS,LAMMINS)
3026 DLAMS=MIN(DLAMS,LAMMAXS)
3027 END IF
3028 !......................................................................
3029 ! GRAUPEL
3031 IF (DUMG(K).GE.QSMALL) THEN
3032 DLAMG = (CONS2*DUMFNG(K)/ DUMG(K))**(1./DG)
3033 DLAMG=MAX(DLAMG,LAMMING)
3034 DLAMG=MIN(DLAMG,LAMMAXG)
3035 END IF
3037 !......................................................................
3038 ! CALCULATE NUMBER-WEIGHTED AND MASS-WEIGHTED TERMINAL FALL SPEEDS
3040 ! CLOUD WATER
3042 IF (DUMC(K).GE.QSMALL) THEN
3043 UNC = ACN(K)*GAMMA(1.+BC+PGAM(K))/ (DLAMC**BC*GAMMA(PGAM(K)+1.))
3044 UMC = ACN(K)*GAMMA(4.+BC+PGAM(K))/ (DLAMC**BC*GAMMA(PGAM(K)+4.))
3045 ELSE
3046 UMC = 0.
3047 UNC = 0.
3048 END IF
3050 IF (DUMI(K).GE.QSMALL) THEN
3051 UNI = AIN(K)*CONS27/DLAMI**BI
3052 UMI = AIN(K)*CONS28/(DLAMI**BI)
3053 ELSE
3054 UMI = 0.
3055 UNI = 0.
3056 END IF
3058 IF (DUMR(K).GE.QSMALL) THEN
3059 UNR = ARN(K)*CONS6/DLAMR**BR
3060 UMR = ARN(K)*CONS4/(DLAMR**BR)
3061 ELSE
3062 UMR = 0.
3063 UNR = 0.
3064 END IF
3066 IF (DUMQS(K).GE.QSMALL) THEN
3067 UMS = ASN(K)*CONS3/(DLAMS**BS)
3068 UNS = ASN(K)*CONS5/DLAMS**BS
3069 ELSE
3070 UMS = 0.
3071 UNS = 0.
3072 END IF
3074 IF (DUMG(K).GE.QSMALL) THEN
3075 UMG = AGN(K)*CONS7/(DLAMG**BG)
3076 UNG = AGN(K)*CONS8/DLAMG**BG
3077 ELSE
3078 UMG = 0.
3079 UNG = 0.
3080 END IF
3082 ! SET REALISTIC LIMITS ON FALLSPEED
3084 UMS=MIN(UMS,1.2)
3085 UNS=MIN(UNS,1.2)
3086 UMI=MIN(UMI,1.2)
3087 UNI=MIN(UNI,1.2)
3088 UMR=MIN(UMR,9.1)
3089 UNR=MIN(UNR,9.1)
3090 UMG=MIN(UMG,20.)
3091 UNG=MIN(UNG,20.)
3093 FR(K) = UMR
3094 FI(K) = UMI
3095 FNI(K) = UNI
3096 FS(K) = UMS
3097 FNS(K) = UNS
3098 FNR(K) = UNR
3099 FC(K) = UMC
3100 FNC(K) = UNC
3101 FG(K) = UMG
3102 FNG(K) = UNG
3104 ! CALCULATE NUMBER OF SPLIT TIME STEPS
3106 RGVM = MAX(FR(K),FI(K),FS(K),FC(K),FNI(K),FNR(K),FNS(K),FNC(K),FG(K),FNG(K))
3107 ! VVT CHANGED IFIX -> INT (GENERIC FUNCTION)
3108 NSTEP = MAX(INT(RGVM*DT/DZQ(K)+1.),NSTEP)
3110 ! MULTIPLY VARIABLES BY RHO
3111 DUMR(k) = DUMR(k)*RHO(K)
3112 DUMI(k) = DUMI(k)*RHO(K)
3113 DUMFNI(k) = DUMFNI(K)*RHO(K)
3114 DUMQS(k) = DUMQS(K)*RHO(K)
3115 DUMFNS(k) = DUMFNS(K)*RHO(K)
3116 DUMFNR(k) = DUMFNR(K)*RHO(K)
3117 DUMC(k) = DUMC(K)*RHO(K)
3118 DUMFNC(k) = DUMFNC(K)*RHO(K)
3119 DUMG(k) = DUMG(K)*RHO(K)
3120 DUMFNG(k) = DUMFNG(K)*RHO(K)
3122 END DO
3124 DO N = 1,NSTEP
3126 DO K = KTS,KTE
3127 FALOUTR(K) = FR(K)*DUMR(K)
3128 FALOUTI(K) = FI(K)*DUMI(K)
3129 FALOUTNI(K) = FNI(K)*DUMFNI(K)
3130 FALOUTS(K) = FS(K)*DUMQS(K)
3131 FALOUTNS(K) = FNS(K)*DUMFNS(K)
3132 FALOUTNR(K) = FNR(K)*DUMFNR(K)
3133 FALOUTC(K) = FC(K)*DUMC(K)
3134 FALOUTNC(K) = FNC(K)*DUMFNC(K)
3135 FALOUTG(K) = FG(K)*DUMG(K)
3136 FALOUTNG(K) = FNG(K)*DUMFNG(K)
3137 END DO
3139 ! TOP OF MODEL
3141 K = KTE
3142 FALTNDR = FALOUTR(K)/DZQ(k)
3143 FALTNDI = FALOUTI(K)/DZQ(k)
3144 FALTNDNI = FALOUTNI(K)/DZQ(k)
3145 FALTNDS = FALOUTS(K)/DZQ(k)
3146 FALTNDNS = FALOUTNS(K)/DZQ(k)
3147 FALTNDNR = FALOUTNR(K)/DZQ(k)
3148 FALTNDC = FALOUTC(K)/DZQ(k)
3149 FALTNDNC = FALOUTNC(K)/DZQ(k)
3150 FALTNDG = FALOUTG(K)/DZQ(k)
3151 FALTNDNG = FALOUTNG(K)/DZQ(k)
3152 ! ADD FALLOUT TERMS TO EULERIAN TENDENCIES
3154 QRSTEN(K) = QRSTEN(K)-FALTNDR/NSTEP/RHO(k)
3155 QISTEN(K) = QISTEN(K)-FALTNDI/NSTEP/RHO(k)
3156 NI3DTEN(K) = NI3DTEN(K)-FALTNDNI/NSTEP/RHO(k)
3157 QNISTEN(K) = QNISTEN(K)-FALTNDS/NSTEP/RHO(k)
3158 NS3DTEN(K) = NS3DTEN(K)-FALTNDNS/NSTEP/RHO(k)
3159 NR3DTEN(K) = NR3DTEN(K)-FALTNDNR/NSTEP/RHO(k)
3160 QCSTEN(K) = QCSTEN(K)-FALTNDC/NSTEP/RHO(k)
3161 NC3DTEN(K) = NC3DTEN(K)-FALTNDNC/NSTEP/RHO(k)
3162 QGSTEN(K) = QGSTEN(K)-FALTNDG/NSTEP/RHO(k)
3163 NG3DTEN(K) = NG3DTEN(K)-FALTNDNG/NSTEP/RHO(k)
3165 DUMR(K) = DUMR(K)-FALTNDR*DT/NSTEP
3166 DUMI(K) = DUMI(K)-FALTNDI*DT/NSTEP
3167 DUMFNI(K) = DUMFNI(K)-FALTNDNI*DT/NSTEP
3168 DUMQS(K) = DUMQS(K)-FALTNDS*DT/NSTEP
3169 DUMFNS(K) = DUMFNS(K)-FALTNDNS*DT/NSTEP
3170 DUMFNR(K) = DUMFNR(K)-FALTNDNR*DT/NSTEP
3171 DUMC(K) = DUMC(K)-FALTNDC*DT/NSTEP
3172 DUMFNC(K) = DUMFNC(K)-FALTNDNC*DT/NSTEP
3173 DUMG(K) = DUMG(K)-FALTNDG*DT/NSTEP
3174 DUMFNG(K) = DUMFNG(K)-FALTNDNG*DT/NSTEP
3176 DO K = KTE-1,KTS,-1
3177 FALTNDR = (FALOUTR(K+1)-FALOUTR(K))/DZQ(K)
3178 FALTNDI = (FALOUTI(K+1)-FALOUTI(K))/DZQ(K)
3179 FALTNDNI = (FALOUTNI(K+1)-FALOUTNI(K))/DZQ(K)
3180 FALTNDS = (FALOUTS(K+1)-FALOUTS(K))/DZQ(K)
3181 FALTNDNS = (FALOUTNS(K+1)-FALOUTNS(K))/DZQ(K)
3182 FALTNDNR = (FALOUTNR(K+1)-FALOUTNR(K))/DZQ(K)
3183 FALTNDC = (FALOUTC(K+1)-FALOUTC(K))/DZQ(K)
3184 FALTNDNC = (FALOUTNC(K+1)-FALOUTNC(K))/DZQ(K)
3185 FALTNDG = (FALOUTG(K+1)-FALOUTG(K))/DZQ(K)
3186 FALTNDNG = (FALOUTNG(K+1)-FALOUTNG(K))/DZQ(K)
3188 ! ADD FALLOUT TERMS TO EULERIAN TENDENCIES
3190 QRSTEN(K) = QRSTEN(K)+FALTNDR/NSTEP/RHO(k)
3191 QISTEN(K) = QISTEN(K)+FALTNDI/NSTEP/RHO(k)
3192 NI3DTEN(K) = NI3DTEN(K)+FALTNDNI/NSTEP/RHO(k)
3193 QNISTEN(K) = QNISTEN(K)+FALTNDS/NSTEP/RHO(k)
3194 NS3DTEN(K) = NS3DTEN(K)+FALTNDNS/NSTEP/RHO(k)
3195 NR3DTEN(K) = NR3DTEN(K)+FALTNDNR/NSTEP/RHO(k)
3196 QCSTEN(K) = QCSTEN(K)+FALTNDC/NSTEP/RHO(k)
3197 NC3DTEN(K) = NC3DTEN(K)+FALTNDNC/NSTEP/RHO(k)
3198 QGSTEN(K) = QGSTEN(K)+FALTNDG/NSTEP/RHO(k)
3199 NG3DTEN(K) = NG3DTEN(K)+FALTNDNG/NSTEP/RHO(k)
3201 DUMR(K) = DUMR(K)+FALTNDR*DT/NSTEP
3202 DUMI(K) = DUMI(K)+FALTNDI*DT/NSTEP
3203 DUMFNI(K) = DUMFNI(K)+FALTNDNI*DT/NSTEP
3204 DUMQS(K) = DUMQS(K)+FALTNDS*DT/NSTEP
3205 DUMFNS(K) = DUMFNS(K)+FALTNDNS*DT/NSTEP
3206 DUMFNR(K) = DUMFNR(K)+FALTNDNR*DT/NSTEP
3207 DUMC(K) = DUMC(K)+FALTNDC*DT/NSTEP
3208 DUMFNC(K) = DUMFNC(K)+FALTNDNC*DT/NSTEP
3209 DUMG(K) = DUMG(K)+FALTNDG*DT/NSTEP
3210 DUMFNG(K) = DUMFNG(K)+FALTNDNG*DT/NSTEP
3212 END DO
3214 ! GET PRECIPITATION AND SNOWFALL ACCUMULATION DURING THE TIME STEP
3215 ! FACTOR OF 1000 CONVERTS FROM M TO MM, BUT DIVISION BY DENSITY
3216 ! OF LIQUID WATER CANCELS THIS FACTOR OF 1000
3218 PRECRT = PRECRT+(FALOUTR(KTS)+FALOUTC(KTS)+FALOUTS(KTS)+FALOUTI(KTS)+FALOUTG(KTS)) &
3219 *DT/NSTEP
3220 SNOWRT = SNOWRT+(FALOUTS(KTS)+FALOUTI(KTS)+FALOUTG(KTS))*DT/NSTEP
3222 END DO
3224 DO K=KTS,KTE
3226 ! ADD ON SEDIMENTATION TENDENCIES FOR MIXING RATIO TO REST OF TENDENCIES
3228 QR3DTEN(K)=QR3DTEN(K)+QRSTEN(K)
3229 QI3DTEN(K)=QI3DTEN(K)+QISTEN(K)
3230 QC3DTEN(K)=QC3DTEN(K)+QCSTEN(K)
3231 QG3DTEN(K)=QG3DTEN(K)+QGSTEN(K)
3232 QNI3DTEN(K)=QNI3DTEN(K)+QNISTEN(K)
3234 ! PUT ALL CLOUD ICE IN SNOW CATEGORY IF MEAN DIAMETER EXCEEDS 2 * dcs
3236 IF (QI3D(K).GE.QSMALL.AND.T3D(K).LT.273.15) THEN
3237 IF (1./LAMI(K).GE.2.*DCS) THEN
3238 QNI3DTEN(K) = QNI3DTEN(K)+QI3D(K)/DT+ QI3DTEN(K)
3239 NS3DTEN(K) = NS3DTEN(K)+NI3D(K)/DT+ NI3DTEN(K)
3240 QI3DTEN(K) = -QI3D(K)/DT
3241 NI3DTEN(K) = -NI3D(K)/DT
3242 END IF
3243 END IF
3245 ! hm add tendencies here, then call sizeparameter
3246 ! to ensure consisitency between mixing ratio and number concentration
3248 QC3D(k) = QC3D(k)+QC3DTEN(k)*DT
3249 QI3D(k) = QI3D(k)+QI3DTEN(k)*DT
3250 QNI3D(k) = QNI3D(k)+QNI3DTEN(k)*DT
3251 QR3D(k) = QR3D(k)+QR3DTEN(k)*DT
3252 ! NC3D(k) = NC3D(k)+NC3DTEN(k)*DT
3253 NI3D(k) = NI3D(k)+NI3DTEN(k)*DT
3254 NS3D(k) = NS3D(k)+NS3DTEN(k)*DT
3255 NR3D(k) = NR3D(k)+NR3DTEN(k)*DT
3257 IF (IGRAUP.EQ.0) THEN
3258 QG3D(k) = QG3D(k)+QG3DTEN(k)*DT
3259 NG3D(k) = NG3D(k)+NG3DTEN(k)*DT
3260 END IF
3262 ! ADD TEMPERATURE AND WATER VAPOR TENDENCIES FROM MICROPHYSICS
3263 T3D(K) = T3D(K)+T3DTEN(k)*DT
3264 QV3D(K) = QV3D(K)+QV3DTEN(k)*DT
3266 ! SATURATION VAPOR PRESSURE AND MIXING RATIO
3268 EVS(K) = POLYSVP(T3D(K),0) ! PA
3269 EIS(K) = POLYSVP(T3D(K),1) ! PA
3271 ! MAKE SURE ICE SATURATION DOESN'T EXCEED WATER SAT. NEAR FREEZING
3273 IF (EIS(K).GT.EVS(K)) EIS(K) = EVS(K)
3275 QVS(K) = .622*EVS(K)/(PRES(K)-EVS(K))
3276 QVI(K) = .622*EIS(K)/(PRES(K)-EIS(K))
3278 QVQVS(K) = QV3D(K)/QVS(K)
3279 QVQVSI(K) = QV3D(K)/QVI(K)
3281 ! AT SUBSATURATION, REMOVE SMALL AMOUNTS OF CLOUD/PRECIP WATER
3283 IF (QVQVS(K).LT.0.9) THEN
3284 IF (QR3D(K).LT.1.E-6) THEN
3285 QV3D(K)=QV3D(K)+QR3D(K)
3286 T3D(K)=T3D(K)-QR3D(K)*XXLV(K)/CPM(K)
3287 QR3D(K)=0.
3288 END IF
3289 IF (QC3D(K).LT.1.E-6) THEN
3290 QV3D(K)=QV3D(K)+QC3D(K)
3291 T3D(K)=T3D(K)-QC3D(K)*XXLV(K)/CPM(K)
3292 QC3D(K)=0.
3293 END IF
3294 END IF
3296 IF (QVQVSI(K).LT.0.9) THEN
3297 IF (QI3D(K).LT.1.E-6) THEN
3298 QV3D(K)=QV3D(K)+QI3D(K)
3299 T3D(K)=T3D(K)-QI3D(K)*XXLS(K)/CPM(K)
3300 QI3D(K)=0.
3301 END IF
3302 IF (QNI3D(K).LT.1.E-6) THEN
3303 QV3D(K)=QV3D(K)+QNI3D(K)
3304 T3D(K)=T3D(K)-QNI3D(K)*XXLS(K)/CPM(K)
3305 QNI3D(K)=0.
3306 END IF
3307 IF (QG3D(K).LT.1.E-6) THEN
3308 QV3D(K)=QV3D(K)+QG3D(K)
3309 T3D(K)=T3D(K)-QG3D(K)*XXLS(K)/CPM(K)
3310 QG3D(K)=0.
3311 END IF
3312 END IF
3314 !..................................................................
3315 ! IF MIXING RATIO < QSMALL SET MIXING RATIO AND NUMBER CONC TO ZERO
3317 IF (QC3D(K).LT.QSMALL) THEN
3318 QC3D(K) = 0.
3319 NC3D(K) = 0.
3320 EFFC(K) = 0.
3321 END IF
3322 IF (QR3D(K).LT.QSMALL) THEN
3323 QR3D(K) = 0.
3324 NR3D(K) = 0.
3325 EFFR(K) = 0.
3326 END IF
3327 IF (QI3D(K).LT.QSMALL) THEN
3328 QI3D(K) = 0.
3329 NI3D(K) = 0.
3330 EFFI(K) = 0.
3331 END IF
3332 IF (QNI3D(K).LT.QSMALL) THEN
3333 QNI3D(K) = 0.
3334 NS3D(K) = 0.
3335 EFFS(K) = 0.
3336 END IF
3337 IF (QG3D(K).LT.QSMALL) THEN
3338 QG3D(K) = 0.
3339 NG3D(K) = 0.
3340 EFFG(K) = 0.
3341 END IF
3343 !..................................
3344 ! IF THERE IS NO CLOUD/PRECIP WATER, THEN SKIP CALCULATIONS
3346 IF (QC3D(K).LT.QSMALL.AND.QI3D(K).LT.QSMALL.AND.QNI3D(K).LT.QSMALL &
3347 .AND.QR3D(K).LT.QSMALL.AND.QG3D(K).LT.QSMALL) GOTO 500
3349 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
3350 ! CALCULATE INSTANTANEOUS PROCESSES
3352 ! ADD MELTING OF CLOUD ICE TO FORM RAIN
3354 IF (QI3D(K).GE.QSMALL.AND.T3D(K).GE.273.15) THEN
3355 QR3D(K) = QR3D(K)+QI3D(K)
3356 T3D(K) = T3D(K)-QI3D(K)*XLF(K)/CPM(K)
3357 QI3D(K) = 0.
3358 NR3D(K) = NR3D(K)+NI3D(K)
3359 NI3D(K) = 0.
3360 END IF
3362 ! ****SENSITIVITY - NO ICE
3363 IF (ILIQ.EQ.1) GOTO 778
3365 ! HOMOGENEOUS FREEZING OF CLOUD WATER
3367 IF (T3D(K).LE.233.15.AND.QC3D(K).GE.QSMALL) THEN
3368 QI3D(K)=QI3D(K)+QC3D(K)
3369 T3D(K)=T3D(K)+QC3D(K)*XLF(K)/CPM(K)
3370 QC3D(K)=0.
3371 NI3D(K)=NI3D(K)+NC3D(K)
3372 NC3D(K)=0.
3373 END IF
3375 ! HOMOGENEOUS FREEZING OF RAIN
3377 IF (IGRAUP.EQ.0) THEN
3379 IF (T3D(K).LE.233.15.AND.QR3D(K).GE.QSMALL) THEN
3380 QG3D(K) = QG3D(K)+QR3D(K)
3381 T3D(K) = T3D(K)+QR3D(K)*XLF(K)/CPM(K)
3382 QR3D(K) = 0.
3383 NG3D(K) = NG3D(K)+ NR3D(K)
3384 NR3D(K) = 0.
3385 END IF
3387 ELSE IF (IGRAUP.EQ.1) THEN
3389 IF (T3D(K).LE.233.15.AND.QR3D(K).GE.QSMALL) THEN
3390 QNI3D(K) = QNI3D(K)+QR3D(K)
3391 T3D(K) = T3D(K)+QR3D(K)*XLF(K)/CPM(K)
3392 QR3D(K) = 0.
3393 NS3D(K) = NS3D(K)+NR3D(K)
3394 NR3D(K) = 0.
3395 END IF
3397 END IF
3399 778 CONTINUE
3401 ! MAKE SURE NUMBER CONCENTRATIONS AREN'T NEGATIVE
3403 NI3D(K) = MAX(0.,NI3D(K))
3404 NS3D(K) = MAX(0.,NS3D(K))
3405 NC3D(K) = MAX(0.,NC3D(K))
3406 NR3D(K) = MAX(0.,NR3D(K))
3407 NG3D(K) = MAX(0.,NG3D(K))
3409 !......................................................................
3410 ! CLOUD ICE
3412 IF (QI3D(K).GE.QSMALL) THEN
3413 LAMI(K) = (CONS12* &
3414 NI3D(K)/QI3D(K))**(1./DI)
3416 ! CHECK FOR SLOPE
3418 ! ADJUST VARS
3420 IF (LAMI(K).LT.LAMMINI) THEN
3422 LAMI(K) = LAMMINI
3424 N0I(K) = LAMI(K)**4*QI3D(K)/CONS12
3426 NI3D(K) = N0I(K)/LAMI(K)
3427 ELSE IF (LAMI(K).GT.LAMMAXI) THEN
3428 LAMI(K) = LAMMAXI
3429 N0I(K) = LAMI(K)**4*QI3D(K)/CONS12
3431 NI3D(K) = N0I(K)/LAMI(K)
3432 END IF
3433 END IF
3435 !......................................................................
3436 ! RAIN
3438 IF (QR3D(K).GE.QSMALL) THEN
3439 LAMR(K) = (PI*RHOW*NR3D(K)/QR3D(K))**(1./3.)
3441 ! CHECK FOR SLOPE
3443 ! ADJUST VARS
3445 IF (LAMR(K).LT.LAMMINR) THEN
3447 LAMR(K) = LAMMINR
3449 N0RR(K) = LAMR(K)**4*QR3D(K)/(PI*RHOW)
3451 NR3D(K) = N0RR(K)/LAMR(K)
3452 ELSE IF (LAMR(K).GT.LAMMAXR) THEN
3453 LAMR(K) = LAMMAXR
3454 N0RR(K) = LAMR(K)**4*QR3D(K)/(PI*RHOW)
3456 NR3D(K) = N0RR(K)/LAMR(K)
3457 END IF
3459 END IF
3461 !......................................................................
3462 ! CLOUD DROPLETS
3464 ! MARTIN ET AL. (1994) FORMULA FOR PGAM
3466 IF (QC3D(K).GE.QSMALL) THEN
3468 DUM = PRES(K)/(287.15*T3D(K))
3469 PGAM(K)=0.0005714*(NC3D(K)/1.E6/DUM)+0.2714
3470 PGAM(K)=1./(PGAM(K)**2)-1.
3471 PGAM(K)=MAX(PGAM(K),2.)
3472 PGAM(K)=MIN(PGAM(K),10.)
3474 ! CALCULATE LAMC
3476 LAMC(K) = (CONS26*NC3D(K)*GAMMA(PGAM(K)+4.)/ &
3477 (QC3D(K)*GAMMA(PGAM(K)+1.)))**(1./3.)
3479 ! LAMMIN, 60 MICRON DIAMETER
3480 ! LAMMAX, 1 MICRON
3482 LAMMIN = (PGAM(K)+1.)/60.E-6
3483 LAMMAX = (PGAM(K)+1.)/1.E-6
3485 IF (LAMC(K).LT.LAMMIN) THEN
3486 LAMC(K) = LAMMIN
3487 NC3D(K) = EXP(3.*LOG(LAMC(K))+LOG(QC3D(K))+ &
3488 LOG(GAMMA(PGAM(K)+1.))-LOG(GAMMA(PGAM(K)+4.)))/CONS26
3490 ELSE IF (LAMC(K).GT.LAMMAX) THEN
3491 LAMC(K) = LAMMAX
3492 NC3D(K) = EXP(3.*LOG(LAMC(K))+LOG(QC3D(K))+ &
3493 LOG(GAMMA(PGAM(K)+1.))-LOG(GAMMA(PGAM(K)+4.)))/CONS26
3495 END IF
3497 END IF
3499 !......................................................................
3500 ! SNOW
3502 IF (QNI3D(K).GE.QSMALL) THEN
3503 LAMS(K) = (CONS1*NS3D(K)/QNI3D(K))**(1./DS)
3505 ! CHECK FOR SLOPE
3507 ! ADJUST VARS
3509 IF (LAMS(K).LT.LAMMINS) THEN
3510 LAMS(K) = LAMMINS
3511 N0S(K) = LAMS(K)**4*QNI3D(K)/CONS1
3513 NS3D(K) = N0S(K)/LAMS(K)
3515 ELSE IF (LAMS(K).GT.LAMMAXS) THEN
3517 LAMS(K) = LAMMAXS
3518 N0S(K) = LAMS(K)**4*QNI3D(K)/CONS1
3519 NS3D(K) = N0S(K)/LAMS(K)
3520 END IF
3522 END IF
3524 !......................................................................
3525 ! GRAUPEL
3527 IF (QG3D(K).GE.QSMALL) THEN
3528 LAMG(K) = (CONS2*NG3D(K)/QG3D(K))**(1./DG)
3530 ! CHECK FOR SLOPE
3532 ! ADJUST VARS
3534 IF (LAMG(K).LT.LAMMING) THEN
3535 LAMG(K) = LAMMING
3536 N0G(K) = LAMG(K)**4*QG3D(K)/CONS2
3538 NG3D(K) = N0G(K)/LAMG(K)
3540 ELSE IF (LAMG(K).GT.LAMMAXG) THEN
3542 LAMG(K) = LAMMAXG
3543 N0G(K) = LAMG(K)**4*QG3D(K)/CONS2
3545 NG3D(K) = N0G(K)/LAMG(K)
3546 END IF
3548 END IF
3550 500 CONTINUE
3552 ! CALCULATE EFFECTIVE RADIUS
3554 IF (QI3D(K).GE.QSMALL) THEN
3555 EFFI(K) = 3./LAMI(K)/2.*1.E6
3556 ELSE
3557 EFFI(K) = 25.
3558 END IF
3560 IF (QNI3D(K).GE.QSMALL) THEN
3561 EFFS(K) = 3./LAMS(K)/2.*1.E6
3562 ELSE
3563 EFFS(K) = 25.
3564 END IF
3566 IF (QR3D(K).GE.QSMALL) THEN
3567 EFFR(K) = 3./LAMR(K)/2.*1.E6
3568 ELSE
3569 EFFR(K) = 25.
3570 END IF
3572 IF (QC3D(K).GE.QSMALL) THEN
3573 EFFC(K) = GAMMA(PGAM(K)+4.)/ &
3574 GAMMA(PGAM(K)+3.)/LAMC(K)/2.*1.E6
3575 ELSE
3576 EFFC(K) = 25.
3577 END IF
3579 IF (QG3D(K).GE.QSMALL) THEN
3580 EFFG(K) = 3./LAMG(K)/2.*1.E6
3581 ELSE
3582 EFFG(K) = 25.
3583 END IF
3585 ! HM ADD 1/10/06, ADD UPPER BOUND ON ICE NUMBER, THIS IS NEEDED
3586 ! TO PREVENT VERY LARGE ICE NUMBER DUE TO HOMOGENEOUS FREEZING
3587 ! OF DROPLETS, ESPECIALLY WHEN INUM = 1, SET MAX AT 10 CM-3
3588 NI3D(K) = MIN(NI3D(K),10.E6/RHO(K))
3589 ! ADD BOUND ON DROPLET NUMBER - CANNOT EXCEED AEROSOL CONCENTRATION
3590 IF (INUM.EQ.0.AND.IACT.EQ.2) THEN
3591 NC3D(K) = MIN(NC3D(K),(NANEW1+NANEW2)/RHO(K))
3592 END IF
3593 ! SWITCH FOR CONSTANT DROPLET NUMBER
3594 ! IF (INUM.EQ.1) THEN
3595 ! CHANGE NDCNST FROM CM-3 TO KG-1
3596 NC3D(K) = NDCNST*1.E6/RHO(K)
3597 ! END IF
3599 END DO !!! K LOOP
3601 400 CONTINUE
3603 ! ALL DONE !!!!!!!!!!!
3604 RETURN
3605 END SUBROUTINE MORR_TWO_MOMENT_MICRO
3607 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
3609 REAL FUNCTION POLYSVP (T,TYPE)
3611 !-------------------------------------------
3613 ! COMPUTE SATURATION VAPOR PRESSURE
3615 ! POLYSVP RETURNED IN UNITS OF PA.
3616 ! T IS INPUT IN UNITS OF K.
3617 ! TYPE REFERS TO SATURATION WITH RESPECT TO LIQUID (0) OR ICE (1)
3619 ! REPLACE GOFF-GRATCH WITH FASTER FORMULATION FROM MARAT KHROUTDINOV
3621 IMPLICIT NONE
3623 REAL DUM
3624 REAL T
3625 INTEGER TYPE
3626 ! ice
3627 real a0i,a1i,a2i,a3i,a4i,a5i,a6i,a7i,a8i
3628 data a0i,a1i,a2i,a3i,a4i,a5i,a6i,a7i,a8i /&
3629 6.11147274, 0.503160820, 0.188439774e-1, &
3630 0.420895665e-3, 0.615021634e-5,0.602588177e-7, &
3631 0.385852041e-9, 0.146898966e-11, 0.252751365e-14/
3633 ! liquid
3634 real a0,a1,a2,a3,a4,a5,a6,a7,a8
3635 data a0,a1,a2,a3,a4,a5,a6,a7,a8 /&
3636 6.105851, 0.4440316, 0.1430341e-1, &
3637 0.2641412e-3, 0.2995057e-5, 0.2031998e-7, &
3638 0.6936113e-10, 0.2564861e-13,-0.3704404e-15/
3639 real dt
3641 ! ICE
3643 IF (TYPE.EQ.1) THEN
3645 ! POLYSVP = 10.**(-9.09718*(273.16/T-1.)-3.56654* &
3646 ! LOG10(273.16/T)+0.876793*(1.-T/273.16)+ &
3647 ! LOG10(6.1071))*100.
3650 dt = max(-80.,t-273.16)
3651 polysvp = a0i + dt*(a1i+dt*(a2i+dt*(a3i+dt*(a4i+dt*(a5i+dt*(a6i+dt*(a7i+a8i*dt)))))))
3652 polysvp = polysvp*100.
3654 END IF
3656 ! LIQUID
3658 IF (TYPE.EQ.0) THEN
3660 dt = max(-80.,t-273.16)
3661 polysvp = a0 + dt*(a1+dt*(a2+dt*(a3+dt*(a4+dt*(a5+dt*(a6+dt*(a7+a8*dt)))))))
3662 polysvp = polysvp*100.
3664 ! POLYSVP = 10.**(-7.90298*(373.16/T-1.)+ &
3665 ! 5.02808*LOG10(373.16/T)- &
3666 ! 1.3816E-7*(10**(11.344*(1.-T/373.16))-1.)+ &
3667 ! 8.1328E-3*(10**(-3.49149*(373.16/T-1.))-1.)+ &
3668 ! LOG10(1013.246))*100.
3670 END IF
3673 END FUNCTION POLYSVP
3675 !------------------------------------------------------------------------------
3677 REAL FUNCTION GAMMA(X)
3678 !----------------------------------------------------------------------
3679 !
3680 ! THIS ROUTINE CALCULATES THE GAMMA FUNCTION FOR A REAL ARGUMENT X.
3681 ! COMPUTATION IS BASED ON AN ALGORITHM OUTLINED IN REFERENCE 1.
3682 ! THE PROGRAM USES RATIONAL FUNCTIONS THAT APPROXIMATE THE GAMMA
3683 ! FUNCTION TO AT LEAST 20 SIGNIFICANT DECIMAL DIGITS. COEFFICIENTS
3684 ! FOR THE APPROXIMATION OVER THE INTERVAL (1,2) ARE UNPUBLISHED.
3685 ! THOSE FOR THE APPROXIMATION FOR X .GE. 12 ARE FROM REFERENCE 2.
3686 ! THE ACCURACY ACHIEVED DEPENDS ON THE ARITHMETIC SYSTEM, THE
3687 ! COMPILER, THE INTRINSIC FUNCTIONS, AND PROPER SELECTION OF THE
3688 ! MACHINE-DEPENDENT CONSTANTS.
3689 !
3690 !
3691 !*******************************************************************
3692 !*******************************************************************
3693 !
3694 ! EXPLANATION OF MACHINE-DEPENDENT CONSTANTS
3695 !
3696 ! BETA - RADIX FOR THE FLOATING-POINT REPRESENTATION
3697 ! MAXEXP - THE SMALLEST POSITIVE POWER OF BETA THAT OVERFLOWS
3698 ! XBIG - THE LARGEST ARGUMENT FOR WHICH GAMMA(X) IS REPRESENTABLE
3699 ! IN THE MACHINE, I.E., THE SOLUTION TO THE EQUATION
3700 ! GAMMA(XBIG) = BETA**MAXEXP
3701 ! XINF - THE LARGEST MACHINE REPRESENTABLE FLOATING-POINT NUMBER;
3702 ! APPROXIMATELY BETA**MAXEXP
3703 ! EPS - THE SMALLEST POSITIVE FLOATING-POINT NUMBER SUCH THAT
3704 ! 1.0+EPS .GT. 1.0
3705 ! XMININ - THE SMALLEST POSITIVE FLOATING-POINT NUMBER SUCH THAT
3706 ! 1/XMININ IS MACHINE REPRESENTABLE
3707 !
3708 ! APPROXIMATE VALUES FOR SOME IMPORTANT MACHINES ARE:
3709 !
3710 ! BETA MAXEXP XBIG
3711 !
3712 ! CRAY-1 (S.P.) 2 8191 966.961
3713 ! CYBER 180/855
3714 ! UNDER NOS (S.P.) 2 1070 177.803
3715 ! IEEE (IBM/XT,
3716 ! SUN, ETC.) (S.P.) 2 128 35.040
3717 ! IEEE (IBM/XT,
3718 ! SUN, ETC.) (D.P.) 2 1024 171.624
3719 ! IBM 3033 (D.P.) 16 63 57.574
3720 ! VAX D-FORMAT (D.P.) 2 127 34.844
3721 ! VAX G-FORMAT (D.P.) 2 1023 171.489
3722 !
3723 ! XINF EPS XMININ
3724 !
3725 ! CRAY-1 (S.P.) 5.45E+2465 7.11E-15 1.84E-2466
3726 ! CYBER 180/855
3727 ! UNDER NOS (S.P.) 1.26E+322 3.55E-15 3.14E-294
3728 ! IEEE (IBM/XT,
3729 ! SUN, ETC.) (S.P.) 3.40E+38 1.19E-7 1.18E-38
3730 ! IEEE (IBM/XT,
3731 ! SUN, ETC.) (D.P.) 1.79D+308 2.22D-16 2.23D-308
3732 ! IBM 3033 (D.P.) 7.23D+75 2.22D-16 1.39D-76
3733 ! VAX D-FORMAT (D.P.) 1.70D+38 1.39D-17 5.88D-39
3734 ! VAX G-FORMAT (D.P.) 8.98D+307 1.11D-16 1.12D-308
3735 !
3736 !*******************************************************************
3737 !*******************************************************************
3738 !
3739 ! ERROR RETURNS
3740 !
3741 ! THE PROGRAM RETURNS THE VALUE XINF FOR SINGULARITIES OR
3742 ! WHEN OVERFLOW WOULD OCCUR. THE COMPUTATION IS BELIEVED
3743 ! TO BE FREE OF UNDERFLOW AND OVERFLOW.
3744 !
3745 !
3746 ! INTRINSIC FUNCTIONS REQUIRED ARE:
3747 !
3748 ! INT, DBLE, EXP, LOG, REAL, SIN
3749 !
3750 !
3751 ! REFERENCES: AN OVERVIEW OF SOFTWARE DEVELOPMENT FOR SPECIAL
3752 ! FUNCTIONS W. J. CODY, LECTURE NOTES IN MATHEMATICS,
3753 ! 506, NUMERICAL ANALYSIS DUNDEE, 1975, G. A. WATSON
3754 ! (ED.), SPRINGER VERLAG, BERLIN, 1976.
3755 !
3756 ! COMPUTER APPROXIMATIONS, HART, ET. AL., WILEY AND
3757 ! SONS, NEW YORK, 1968.
3758 !
3759 ! LATEST MODIFICATION: OCTOBER 12, 1989
3760 !
3761 ! AUTHORS: W. J. CODY AND L. STOLTZ
3762 ! APPLIED MATHEMATICS DIVISION
3763 ! ARGONNE NATIONAL LABORATORY
3764 ! ARGONNE, IL 60439
3765 !
3766 !----------------------------------------------------------------------
3767 implicit none
3768 INTEGER I,N
3769 LOGICAL PARITY
3770 REAL &
3771 CONV,EPS,FACT,HALF,ONE,RES,SUM,TWELVE, &
3772 TWO,X,XBIG,XDEN,XINF,XMININ,XNUM,Y,Y1,YSQ,Z,ZERO
3773 REAL, DIMENSION(7) :: C
3774 REAL, DIMENSION(8) :: P
3775 REAL, DIMENSION(8) :: Q
3776 !----------------------------------------------------------------------
3777 ! MATHEMATICAL CONSTANTS
3778 !----------------------------------------------------------------------
3779 DATA ONE,HALF,TWELVE,TWO,ZERO/1.0E0,0.5E0,12.0E0,2.0E0,0.0E0/
3782 !----------------------------------------------------------------------
3783 ! MACHINE DEPENDENT PARAMETERS
3784 !----------------------------------------------------------------------
3785 DATA XBIG,XMININ,EPS/35.040E0,1.18E-38,1.19E-7/,XINF/3.4E38/
3786 !----------------------------------------------------------------------
3787 ! NUMERATOR AND DENOMINATOR COEFFICIENTS FOR RATIONAL MINIMAX
3788 ! APPROXIMATION OVER (1,2).
3789 !----------------------------------------------------------------------
3790 DATA P/-1.71618513886549492533811E+0,2.47656508055759199108314E+1, &
3791 -3.79804256470945635097577E+2,6.29331155312818442661052E+2, &
3792 8.66966202790413211295064E+2,-3.14512729688483675254357E+4, &
3793 -3.61444134186911729807069E+4,6.64561438202405440627855E+4/
3794 DATA Q/-3.08402300119738975254353E+1,3.15350626979604161529144E+2, &
3795 -1.01515636749021914166146E+3,-3.10777167157231109440444E+3, &
3796 2.25381184209801510330112E+4,4.75584627752788110767815E+3, &
3797 -1.34659959864969306392456E+5,-1.15132259675553483497211E+5/
3798 !----------------------------------------------------------------------
3799 ! COEFFICIENTS FOR MINIMAX APPROXIMATION OVER (12, INF).
3800 !----------------------------------------------------------------------
3801 DATA C/-1.910444077728E-03,8.4171387781295E-04, &
3802 -5.952379913043012E-04,7.93650793500350248E-04, &
3803 -2.777777777777681622553E-03,8.333333333333333331554247E-02, &
3804 5.7083835261E-03/
3805 !----------------------------------------------------------------------
3806 ! STATEMENT FUNCTIONS FOR CONVERSION BETWEEN INTEGER AND FLOAT
3807 !----------------------------------------------------------------------
3808 CONV(I) = REAL(I)
3809 PARITY=.FALSE.
3810 FACT=ONE
3811 N=0
3812 Y=X
3813 IF(Y.LE.ZERO)THEN
3814 !----------------------------------------------------------------------
3815 ! ARGUMENT IS NEGATIVE
3816 !----------------------------------------------------------------------
3817 Y=-X
3818 Y1=AINT(Y)
3819 RES=Y-Y1
3820 IF(RES.NE.ZERO)THEN
3821 IF(Y1.NE.AINT(Y1*HALF)*TWO)PARITY=.TRUE.
3822 FACT=-PI/SIN(PI*RES)
3823 Y=Y+ONE
3824 ELSE
3825 RES=XINF
3826 GOTO 900
3827 ENDIF
3828 ENDIF
3829 !----------------------------------------------------------------------
3830 ! ARGUMENT IS POSITIVE
3831 !----------------------------------------------------------------------
3832 IF(Y.LT.EPS)THEN
3833 !----------------------------------------------------------------------
3834 ! ARGUMENT .LT. EPS
3835 !----------------------------------------------------------------------
3836 IF(Y.GE.XMININ)THEN
3837 RES=ONE/Y
3838 ELSE
3839 RES=XINF
3840 GOTO 900
3841 ENDIF
3842 ELSEIF(Y.LT.TWELVE)THEN
3843 Y1=Y
3844 IF(Y.LT.ONE)THEN
3845 !----------------------------------------------------------------------
3846 ! 0.0 .LT. ARGUMENT .LT. 1.0
3847 !----------------------------------------------------------------------
3848 Z=Y
3849 Y=Y+ONE
3850 ELSE
3851 !----------------------------------------------------------------------
3852 ! 1.0 .LT. ARGUMENT .LT. 12.0, REDUCE ARGUMENT IF NECESSARY
3853 !----------------------------------------------------------------------
3854 N=INT(Y)-1
3855 Y=Y-CONV(N)
3856 Z=Y-ONE
3857 ENDIF
3858 !----------------------------------------------------------------------
3859 ! EVALUATE APPROXIMATION FOR 1.0 .LT. ARGUMENT .LT. 2.0
3860 !----------------------------------------------------------------------
3861 XNUM=ZERO
3862 XDEN=ONE
3863 DO I=1,8
3864 XNUM=(XNUM+P(I))*Z
3865 XDEN=XDEN*Z+Q(I)
3866 END DO
3867 RES=XNUM/XDEN+ONE
3868 IF(Y1.LT.Y)THEN
3869 !----------------------------------------------------------------------
3870 ! ADJUST RESULT FOR CASE 0.0 .LT. ARGUMENT .LT. 1.0
3871 !----------------------------------------------------------------------
3872 RES=RES/Y1
3873 ELSEIF(Y1.GT.Y)THEN
3874 !----------------------------------------------------------------------
3875 ! ADJUST RESULT FOR CASE 2.0 .LT. ARGUMENT .LT. 12.0
3876 !----------------------------------------------------------------------
3877 DO I=1,N
3878 RES=RES*Y
3879 Y=Y+ONE
3880 END DO
3881 ENDIF
3882 ELSE
3883 !----------------------------------------------------------------------
3884 ! EVALUATE FOR ARGUMENT .GE. 12.0,
3885 !----------------------------------------------------------------------
3886 IF(Y.LE.XBIG)THEN
3887 YSQ=Y*Y
3888 SUM=C(7)
3889 DO I=1,6
3890 SUM=SUM/YSQ+C(I)
3891 END DO
3892 SUM=SUM/Y-Y+SQRTPI
3893 SUM=SUM+(Y-HALF)*LOG(Y)
3894 RES=EXP(SUM)
3895 ELSE
3896 RES=XINF
3897 GOTO 900
3898 ENDIF
3899 ENDIF
3900 !----------------------------------------------------------------------
3901 ! FINAL ADJUSTMENTS AND RETURN
3902 !----------------------------------------------------------------------
3903 IF(PARITY)RES=-RES
3904 IF(FACT.NE.ONE)RES=FACT/RES
3905 900 GAMMA=RES
3906 RETURN
3907 ! ---------- LAST LINE OF GAMMA ----------
3908 END FUNCTION GAMMA
3911 REAL FUNCTION DERF1(X)
3912 IMPLICIT NONE
3913 REAL X
3914 REAL, DIMENSION(0 : 64) :: A, B
3915 REAL W,T,Y
3916 INTEGER K,I
3917 DATA A/ &
3918 0.00000000005958930743E0, -0.00000000113739022964E0, &
3919 0.00000001466005199839E0, -0.00000016350354461960E0, &
3920 0.00000164610044809620E0, -0.00001492559551950604E0, &
3921 0.00012055331122299265E0, -0.00085483269811296660E0, &
3922 0.00522397762482322257E0, -0.02686617064507733420E0, &
3923 0.11283791670954881569E0, -0.37612638903183748117E0, &
3924 1.12837916709551257377E0, &
3925 0.00000000002372510631E0, -0.00000000045493253732E0, &
3926 0.00000000590362766598E0, -0.00000006642090827576E0, &
3927 0.00000067595634268133E0, -0.00000621188515924000E0, &
3928 0.00005103883009709690E0, -0.00037015410692956173E0, &
3929 0.00233307631218880978E0, -0.01254988477182192210E0, &
3930 0.05657061146827041994E0, -0.21379664776456006580E0, &
3931 0.84270079294971486929E0, &
3932 0.00000000000949905026E0, -0.00000000018310229805E0, &
3933 0.00000000239463074000E0, -0.00000002721444369609E0, &
3934 0.00000028045522331686E0, -0.00000261830022482897E0, &
3935 0.00002195455056768781E0, -0.00016358986921372656E0, &
3936 0.00107052153564110318E0, -0.00608284718113590151E0, &
3937 0.02986978465246258244E0, -0.13055593046562267625E0, &
3938 0.67493323603965504676E0, &
3939 0.00000000000382722073E0, -0.00000000007421598602E0, &
3940 0.00000000097930574080E0, -0.00000001126008898854E0, &
3941 0.00000011775134830784E0, -0.00000111992758382650E0, &
3942 0.00000962023443095201E0, -0.00007404402135070773E0, &
3943 0.00050689993654144881E0, -0.00307553051439272889E0, &
3944 0.01668977892553165586E0, -0.08548534594781312114E0, &
3945 0.56909076642393639985E0, &
3946 0.00000000000155296588E0, -0.00000000003032205868E0, &
3947 0.00000000040424830707E0, -0.00000000471135111493E0, &
3948 0.00000005011915876293E0, -0.00000048722516178974E0, &
3949 0.00000430683284629395E0, -0.00003445026145385764E0, &
3950 0.00024879276133931664E0, -0.00162940941748079288E0, &
3951 0.00988786373932350462E0, -0.05962426839442303805E0, &
3952 0.49766113250947636708E0 /
3953 DATA (B(I), I = 0, 12) / &
3954 -0.00000000029734388465E0, 0.00000000269776334046E0, &
3955 -0.00000000640788827665E0, -0.00000001667820132100E0, &
3956 -0.00000021854388148686E0, 0.00000266246030457984E0, &
3957 0.00001612722157047886E0, -0.00025616361025506629E0, &
3958 0.00015380842432375365E0, 0.00815533022524927908E0, &
3959 -0.01402283663896319337E0, -0.19746892495383021487E0, &
3960 0.71511720328842845913E0 /
3961 DATA (B(I), I = 13, 25) / &
3962 -0.00000000001951073787E0, -0.00000000032302692214E0, &
3963 0.00000000522461866919E0, 0.00000000342940918551E0, &
3964 -0.00000035772874310272E0, 0.00000019999935792654E0, &
3965 0.00002687044575042908E0, -0.00011843240273775776E0, &
3966 -0.00080991728956032271E0, 0.00661062970502241174E0, &
3967 0.00909530922354827295E0, -0.20160072778491013140E0, &
3968 0.51169696718727644908E0 /
3969 DATA (B(I), I = 26, 38) / &
3970 0.00000000003147682272E0, -0.00000000048465972408E0, &
3971 0.00000000063675740242E0, 0.00000003377623323271E0, &
3972 -0.00000015451139637086E0, -0.00000203340624738438E0, &
3973 0.00001947204525295057E0, 0.00002854147231653228E0, &
3974 -0.00101565063152200272E0, 0.00271187003520095655E0, &
3975 0.02328095035422810727E0, -0.16725021123116877197E0, &
3976 0.32490054966649436974E0 /
3977 DATA (B(I), I = 39, 51) / &
3978 0.00000000002319363370E0, -0.00000000006303206648E0, &
3979 -0.00000000264888267434E0, 0.00000002050708040581E0, &
3980 0.00000011371857327578E0, -0.00000211211337219663E0, &
3981 0.00000368797328322935E0, 0.00009823686253424796E0, &
3982 -0.00065860243990455368E0, -0.00075285814895230877E0, &
3983 0.02585434424202960464E0, -0.11637092784486193258E0, &
3984 0.18267336775296612024E0 /
3985 DATA (B(I), I = 52, 64) / &
3986 -0.00000000000367789363E0, 0.00000000020876046746E0, &
3987 -0.00000000193319027226E0, -0.00000000435953392472E0, &
3988 0.00000018006992266137E0, -0.00000078441223763969E0, &
3989 -0.00000675407647949153E0, 0.00008428418334440096E0, &
3990 -0.00017604388937031815E0, -0.00239729611435071610E0, &
3991 0.02064129023876022970E0, -0.06905562880005864105E0, &
3992 0.09084526782065478489E0 /
3993 W = ABS(X)
3994 IF (W .LT. 2.2D0) THEN
3995 T = W * W
3996 K = INT(T)
3997 T = T - K
3998 K = K * 13
3999 Y = ((((((((((((A(K) * T + A(K + 1)) * T + &
4000 A(K + 2)) * T + A(K + 3)) * T + A(K + 4)) * T + &
4001 A(K + 5)) * T + A(K + 6)) * T + A(K + 7)) * T + &
4002 A(K + 8)) * T + A(K + 9)) * T + A(K + 10)) * T + &
4003 A(K + 11)) * T + A(K + 12)) * W
4004 ELSE IF (W .LT. 6.9D0) THEN
4005 K = INT(W)
4006 T = W - K
4007 K = 13 * (K - 2)
4008 Y = (((((((((((B(K) * T + B(K + 1)) * T + &
4009 B(K + 2)) * T + B(K + 3)) * T + B(K + 4)) * T + &
4010 B(K + 5)) * T + B(K + 6)) * T + B(K + 7)) * T + &
4011 B(K + 8)) * T + B(K + 9)) * T + B(K + 10)) * T + &
4012 B(K + 11)) * T + B(K + 12)
4013 Y = Y * Y
4014 Y = Y * Y
4015 Y = Y * Y
4016 Y = 1 - Y * Y
4017 ELSE
4018 Y = 1
4019 END IF
4020 IF (X .LT. 0) Y = -Y
4021 DERF1 = Y
4022 END FUNCTION DERF1
4024 !+---+-----------------------------------------------------------------+
4026 END MODULE module_mp_morr_two_moment