module_mp_morr_two_moment.F
References to this file elsewhere.
1 !WRF:MODEL_LAYER:PHYSICS
2 !
3
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)
7
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.
11
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
16
17 ! USE module_state_description
18
19 IMPLICIT NONE
20
21 REAL, PARAMETER :: PI = 3.1415926535897932384626434
22 REAL, PARAMETER :: SQRTPI = 0.9189385332046727417803297
23
24 PUBLIC :: MP_MORR_TWO_MOMENT
25 PUBLIC :: POLYSVP
26
27 PRIVATE :: GAMMA, DERF1
28 PRIVATE :: PI, SQRTPI
29 PRIVATE :: MORR_TWO_MOMENT_MICRO
30
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
35
36 INTEGER, PRIVATE :: IACT
37
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
42
43 INTEGER, PRIVATE :: INUM
44
45 ! FOR INUM = 1, SET CONSTANT DROPLET CONCENTRATION (CM-3)
46 REAL, PRIVATE :: NDCNST
47
48 ! SWITCH FOR LIQUID-ONLY RUN
49 ! ILIQ = 0, INCLUDE ICE
50 ! ILIQ = 1, LIQUID ONLY, NO ICE
51
52 INTEGER, PRIVATE :: ILIQ
53
54 ! SWITCH FOR ICE NUCLEATION
55 ! INUC = 0, USE FORMULA FROM RASMUSSEN ET AL. 2002 (MID-LATITUDE)
56 ! = 1, USE MPACE OBSERVATIONS
57
58 INTEGER, PRIVATE :: INUC
59
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
71
72 ! NOTE: ONLY USED FOR PREDICTED DROPLET CONCENTRATION (INUM = 0)
73
74 INTEGER, PRIVATE :: IBASE
75
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
79
80 INTEGER, PRIVATE :: ISUB
81
82 ! SWITCH FOR GRAUPEL/NO GRAUPEL
83 ! IGRAUP = 0, INCLUDE GRAUPEL
84 ! IGRAUP = 1, NO GRAUPEL
85
86 INTEGER, PRIVATE :: IGRAUP
87
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
92
93 INTEGER, PRIVATE :: IHAIL
94
95 ! CLOUD MICROPHYSICS CONSTANTS
96
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)
123
124 ! CCN SPECTRA FOR IACT = 1
125
126 REAL, PRIVATE :: C1 ! 'C' IN NCCN = CS^K (CM-3)
127 REAL, PRIVATE :: K1 ! 'K' IN NCCN = CS^K
128
129 ! AEROSOL PARAMETERS FOR IACT = 2
130
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
152
153 ! CONSTANTS TO IMPROVE EFFICIENCY
154
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
160
161 CONTAINS
162
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 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
170
171 IMPLICIT NONE
172
173 integer n,i
174
175 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
176 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
177
178 ! THE FOLLOWING PARAMETERS ARE USER-DEFINED SWITCHES AND NEED TO BE
179 ! SET PRIOR TO CODE COMPILATION
180
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
186
187 INUM = 1
188
189 ! FOR INUM = 1, SET CONSTANT DROPLET CONCENTRATION (UNITS OF CM-3)
190
191 NDCNST = 100.
192
193 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
194 ! NOTE, FOLLOWING OPTIONS NOT AVAILABLE IN CURRENT VERSION
195 ! ONLY USED WHEN INUM=0
196
197
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)
201
202 IACT = 2
203
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
215
216 ! NOTE: ONLY USED FOR PREDICTED DROPLET CONCENTRATION (INUM = 0)
217
218 IBASE = 2
219
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
224
225 ! NOTE: ONLY USED FOR PREDICTED DROPLET CONCENTRATION (INUM = 0)
226
227 ISUB = 0
228 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
229
230
231 ! SWITCH FOR LIQUID-ONLY RUN
232 ! ILIQ = 0, INCLUDE ICE
233 ! ILIQ = 1, LIQUID ONLY, NO ICE
234
235 ILIQ = 0
236
237 ! SWITCH FOR ICE NUCLEATION
238 ! INUC = 0, USE FORMULA FROM RASMUSSEN ET AL. 2002 (MID-LATITUDE)
239 ! = 1, USE MPACE OBSERVATIONS (ARCTIC ONLY)
240
241 INUC = 1
242
243 ! SWITCH FOR GRAUPEL/NO GRAUPEL
244 ! IGRAUP = 0, INCLUDE GRAUPEL
245 ! IGRAUP = 1, NO GRAUPEL
246
247 IGRAUP = 0
248
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
254
255 IHAIL = 0
256
257 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
258 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
259 ! SET PHYSICAL CONSTANTS
260
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
277
278 ! CONSTANTS AND PARAMETERS
279 R = 287.15
280 RV = 465.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
305
306 ! SIZE DISTRIBUTION PARAMETERS
307
308 CI = RHOI*PI/6.
309 DI = 3.
310 CS = RHOSN*PI/6.
311 DS = 3.
312 CG = RHOG*PI/6.
313 DG = 3.
314
315 ! RADIUS OF CONTACT NUCLEI
316 RIN = 0.1E-6
317
318 MMULT = 4./3.*PI*RHOI*(5.E-6)**3
319
320 ! SIZE LIMITS FOR LAMBDA
321
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
330
331 ! CCN SPECTRA FOR IACT = 1
332
333 ! MARITIME
334 ! MODIFIED FROM RASMUSSEN ET AL. 2002
335 ! NCCN = C*S^K, NCCN IS IN CM-3, S IS SUPERSATURATION RATIO IN %
336
337 K1 = 0.4
338 C1 = 120.
339
340 ! CONTINENTAL
341
342 ! K1 = 0.5
343 ! C1 = 1000.
344
345 ! AEROSOL ACTIVATION PARAMETERS FOR IACT = 2
346 ! PARAMETERS CURRENTLY SET FOR AMMONIUM SULFATE
347
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)
357
358 ! AEROSOL SIZE DISTRIBUTION PARAMETERS CURRENTLY SET FOR MPACE
359 ! (see morrison et al. 2007, JGR)
360 ! MODE 1
361
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)
367
368 ! MODE 2
369
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)
375
376 ! CONSTANTS FOR EFFICIENCY
377
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
419
420 END SUBROUTINE MORR_TWO_MOMENT_INIT
421
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.
430
431 ! THIS CODE WAS WRITTEN BY HUGH MORRISON (NCAR) AND SLAVA TATARSKII (GEORGIA TECH).
432
433 ! FOR QUESTIONS, CONTACT: HUGH MORRISON, E-MAIL: MORRISON@UCAR.EDU, PHONE:303-497-8916
434
435 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
436
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 )
446
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
471
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 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
480
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 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
488
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)
493
494 ! ADDITIONAL OUTPUT FROM MICRO - SEDIMENTATION TENDENCIES, NEEDED FOR LIQUID-ICE STATIC ENERGY
495
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)
501
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)
507
508 IMPLICIT NONE
509
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
514
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
518
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
523
524 REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT):: &
525 RAINNC, RAINNCV, SR
526
527 ! REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT):: & ! GT
528 ! refl_10cm
529
530 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN) :: ht
531
532 ! TYPE (WRFU_Clock):: grid_clock ! GT
533 ! TYPE (WRFU_Alarm), POINTER:: grid_alarms(:) ! GT
534
535 ! LOCAL VARIABLES
536
537 REAL, DIMENSION(ims:ime, kms:kme, jms:jme):: &
538 effi, effs, effr, EFFG
539
540 REAL, DIMENSION(ims:ime, kms:kme, jms:jme):: &
541 T, WVAR, EFFC
542
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, &
551
552 ! ADD SEDIMENTATION TENDENCIES (UNITS OF KG/KG/S)
553 QGSTEN,QRSTEN, QISTEN, QNISTEN, QCSTEN, &
554 ! ADD CUMULUS TENDENCIES
555 QRCU1D, QSCU1D, QICU1D
556
557 ! add cumulus tendencies
558
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
563
564 ! HM add reflectivity
565 ! dbz
566
567 REAL PRECPRT1D, SNOWRT1D
568
569 INTEGER I,K,J
570
571 REAL DT
572
573 ! LOGICAL:: dBZ_tstep ! GT
574
575 ! Initialize tendencies (all set to 0) and transfer
576 ! array to local variables
577 DT = DT_IN
578
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)
583
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
591
592 ! currently assign wvar to 0.5 m/s (not coupled with PBL scheme)
593
594 WVAR(I,K,J) = 0.5
595
596 ! currently mixing of number concentrations also is neglected (not coupled with PBL schemes)
597
598 END DO
599 END DO
600 END DO
601
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 !
607
608 ! hm , initialize 1d tendency arrays to zero
609
610 do k=kts,kte ! k loop (vertical)
611
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.
621
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)
626
627 NI1D(k) = NI(i,k,j)
628
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.
636
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
648
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)
661
662 !
663 ! Transfer 1D arrays back into 3D arrays
664 !
665 do k=kts,kte
666
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
670
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)
680
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)
684
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)
690
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.)
697
698 end do
699
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)
704
705 end do
706 end do
707
708 END SUBROUTINE MP_MORR_TWO_MOMENT
709
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)
720
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)
725
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.
729
730 ! CODE STRUCTURE: MAIN SUBROUTINE IS 'MORR_TWO_MOMENT'. ALSO INCLUDED IN THIS FILE IS
731 ! 'FUNCTION POLYSVP', 'FUNCTION DERF1', AND
732 ! 'FUNCTION GAMMA'.
733
734 ! NOTE: THIS SUBROUTINE USES 1D ARRAY IN VERTICAL (COLUMN), EVEN THOUGH VARIABLES ARE CALLED '3D'......
735
736 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
737
738 ! DECLARATIONS
739
740 IMPLICIT NONE
741
742 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
743 ! THESE VARIABLES BELOW MUST BE LINKED WITH THE MAIN MODEL.
744 ! DEFINE ARRAY SIZES
745
746 ! INPUT NUMBER OF GRID CELLS
747
748 ! INPUT/OUTPUT PARAMETERS ! DESCRIPTION (UNITS)
749 INTEGER, INTENT( IN) :: IMS,IME, JMS,JME, KMS,KME, &
750 ITS,ITE, JTS,JTE, KTS,KTE
751
752 REAL, DIMENSION(KMS:KME) :: QC3DTEN ! CLOUD WATER MIXING RATIO TENDENCY (KG/KG/S)
753 REAL, DIMENSION(KMS:KME) :: QI3DTEN ! CLOUD ICE MIXING RATIO TENDENCY (KG/KG/S)
754 REAL, DIMENSION(KMS:KME) :: QNI3DTEN ! SNOW MIXING RATIO TENDENCY (KG/KG/S)
755 REAL, DIMENSION(KMS:KME) :: QR3DTEN ! RAIN MIXING RATIO TENDENCY (KG/KG/S)
756 REAL, DIMENSION(KMS:KME) :: NI3DTEN ! CLOUD ICE NUMBER CONCENTRATION (1/KG/S)
757 REAL, DIMENSION(KMS:KME) :: NS3DTEN ! SNOW NUMBER CONCENTRATION (1/KG/S)
758 REAL, DIMENSION(KMS:KME) :: NR3DTEN ! RAIN NUMBER CONCENTRATION (1/KG/S)
759 REAL, DIMENSION(KMS:KME) :: QC3D ! CLOUD WATER MIXING RATIO (KG/KG)
760 REAL, DIMENSION(KMS:KME) :: QI3D ! CLOUD ICE MIXING RATIO (KG/KG)
761 REAL, DIMENSION(KMS:KME) :: QNI3D ! SNOW MIXING RATIO (KG/KG)
762 REAL, DIMENSION(KMS:KME) :: QR3D ! RAIN MIXING RATIO (KG/KG)
763 REAL, DIMENSION(KMS:KME) :: NI3D ! CLOUD ICE NUMBER CONCENTRATION (1/KG)
764 REAL, DIMENSION(KMS:KME) :: NS3D ! SNOW NUMBER CONCENTRATION (1/KG)
765 REAL, DIMENSION(KMS:KME) :: NR3D ! RAIN NUMBER CONCENTRATION (1/KG)
766 REAL, DIMENSION(KMS:KME) :: T3DTEN ! TEMPERATURE TENDENCY (K/S)
767 REAL, DIMENSION(KMS:KME) :: QV3DTEN ! WATER VAPOR MIXING RATIO TENDENCY (KG/KG/S)
768 REAL, DIMENSION(KMS:KME) :: T3D ! TEMPERATURE (K)
769 REAL, DIMENSION(KMS:KME) :: QV3D ! WATER VAPOR MIXING RATIO (KG/KG)
770 REAL, DIMENSION(KMS:KME) :: PRES ! ATMOSPHERIC PRESSURE (PA)
771 REAL, DIMENSION(KMS:KME) :: DZQ ! DIFFERENCE IN HEIGHT ACROSS LEVEL (m)
772 REAL, DIMENSION(KMS:KME) :: W3D ! GRID-SCALE VERTICAL VELOCITY (M/S)
773 REAL, DIMENSION(KMS:KME) :: WVAR ! SUB-GRID VERTICAL VELOCITY (M/S)
774
775 ! HM ADDED GRAUPEL VARIABLES
776 REAL, DIMENSION(KMS:KME) :: QG3DTEN ! GRAUPEL MIX RATIO TENDENCY (KG/KG/S)
777 REAL, DIMENSION(KMS:KME) :: NG3DTEN ! GRAUPEL NUMB CONC TENDENCY (1/KG/S)
778 REAL, DIMENSION(KMS:KME) :: QG3D ! GRAUPEL MIX RATIO (KG/KG)
779 REAL, DIMENSION(KMS:KME) :: NG3D ! GRAUPEL NUMBER CONC (1/KG)
780
781 ! HM, ADD 1/16/07, SEDIMENTATION TENDENCIES FOR MIXING RATIO
782
783 REAL, DIMENSION(KMS:KME) :: QGSTEN ! GRAUPEL SED TEND (KG/KG/S)
784 REAL, DIMENSION(KMS:KME) :: QRSTEN ! RAIN SED TEND (KG/KG/S)
785 REAL, DIMENSION(KMS:KME) :: QISTEN ! CLOUD ICE SED TEND (KG/KG/S)
786 REAL, DIMENSION(KMS:KME) :: QNISTEN ! SNOW SED TEND (KG/KG/S)
787 REAL, DIMENSION(KMS:KME) :: QCSTEN ! CLOUD WAT SED TEND (KG/KG/S)
788
789 ! hm add cumulus tendencies for precip
790 REAL, DIMENSION(KMS:KME) :: qrcu1d
791 REAL, DIMENSION(KMS:KME) :: qscu1d
792 REAL, DIMENSION(KMS:KME) :: qicu1d
793
794 ! OUTPUT VARIABLES
795
796 REAL PRECRT ! TOTAL PRECIP PER TIME STEP (mm)
797 REAL SNOWRT ! SNOW PER TIME STEP (mm)
798
799 REAL, DIMENSION(KMS:KME) :: EFFC ! DROPLET EFFECTIVE RADIUS (MICRON)
800 REAL, DIMENSION(KMS:KME) :: EFFI ! CLOUD ICE EFFECTIVE RADIUS (MICRON)
801 REAL, DIMENSION(KMS:KME) :: EFFS ! SNOW EFFECTIVE RADIUS (MICRON)
802 REAL, DIMENSION(KMS:KME) :: EFFR ! RAIN EFFECTIVE RADIUS (MICRON)
803 REAL, DIMENSION(KMS:KME) :: EFFG ! GRAUPEL EFFECTIVE RADIUS (MICRON)
804
805 ! MODEL INPUT PARAMETERS (FORMERLY IN COMMON BLOCKS)
806
807 REAL DT ! MODEL TIME STEP (SEC)
808
809 !.....................................................................................................
810 ! LOCAL VARIABLES: ALL PARAMETERS BELOW ARE LOCAL TO SCHEME AND DON'T NEED TO COMMUNICATE WITH THE
811 ! REST OF THE MODEL.
812
813 ! SIZE PARAMETER VARIABLES
814
815 REAL, DIMENSION(KMS:KME) :: LAMC ! SLOPE PARAMETER FOR DROPLETS (M-1)
816 REAL, DIMENSION(KMS:KME) :: LAMI ! SLOPE PARAMETER FOR CLOUD ICE (M-1)
817 REAL, DIMENSION(KMS:KME) :: LAMS ! SLOPE PARAMETER FOR SNOW (M-1)
818 REAL, DIMENSION(KMS:KME) :: LAMR ! SLOPE PARAMETER FOR RAIN (M-1)
819 REAL, DIMENSION(KMS:KME) :: LAMG ! SLOPE PARAMETER FOR GRAUPEL (M-1)
820 REAL, DIMENSION(KMS:KME) :: CDIST1 ! PSD PARAMETER FOR DROPLETS
821 REAL, DIMENSION(KMS:KME) :: N0I ! INTERCEPT PARAMETER FOR CLOUD ICE (KG-1 M-1)
822 REAL, DIMENSION(KMS:KME) :: N0S ! INTERCEPT PARAMETER FOR SNOW (KG-1 M-1)
823 REAL, DIMENSION(KMS:KME) :: N0RR ! INTERCEPT PARAMETER FOR RAIN (KG-1 M-1)
824 REAL, DIMENSION(KMS:KME) :: N0G ! INTERCEPT PARAMETER FOR GRAUPEL (KG-1 M-1)
825 REAL, DIMENSION(KMS:KME) :: PGAM ! SPECTRAL SHAPE PARAMETER FOR DROPLETS
826
827 ! MICROPHYSICAL PROCESSES
828
829 REAL, DIMENSION(KMS:KME) :: NSUBC ! LOSS OF NC DURING EVAP
830 REAL, DIMENSION(KMS:KME) :: NSUBI ! LOSS OF NI DURING SUB.
831 REAL, DIMENSION(KMS:KME) :: NSUBS ! LOSS OF NS DURING SUB.
832 REAL, DIMENSION(KMS:KME) :: NSUBR ! LOSS OF NR DURING EVAP
833 REAL, DIMENSION(KMS:KME) :: PRD ! DEP CLOUD ICE
834 REAL, DIMENSION(KMS:KME) :: PRE ! EVAP OF RAIN
835 REAL, DIMENSION(KMS:KME) :: PRDS ! DEP SNOW
836 REAL, DIMENSION(KMS:KME) :: NNUCCC ! CHANGE N DUE TO CONTACT FREEZ DROPLETS
837 REAL, DIMENSION(KMS:KME) :: MNUCCC ! CHANGE Q DUE TO CONTACT FREEZ DROPLETS
838 REAL, DIMENSION(KMS:KME) :: PRA ! ACCRETION DROPLETS BY RAIN
839 REAL, DIMENSION(KMS:KME) :: PRC ! AUTOCONVERSION DROPLETS
840 REAL, DIMENSION(KMS:KME) :: PCC ! COND/EVAP DROPLETS
841 REAL, DIMENSION(KMS:KME) :: NNUCCD ! CHANGE N FREEZING AEROSOL (PRIM ICE NUCLEATION)
842 REAL, DIMENSION(KMS:KME) :: MNUCCD ! CHANGE Q FREEZING AEROSOL (PRIM ICE NUCLEATION)
843 REAL, DIMENSION(KMS:KME) :: MNUCCR ! CHANGE Q DUE TO CONTACT FREEZ RAIN
844 REAL, DIMENSION(KMS:KME) :: NNUCCR ! CHANGE N DUE TO CONTACT FREEZ RAIN
845 REAL, DIMENSION(KMS:KME) :: NPRA ! CHANGE IN N DUE TO DROPLET ACC BY RAIN
846 REAL, DIMENSION(KMS:KME) :: NRAGG ! SELF-COLLECTION OF RAIN
847 REAL, DIMENSION(KMS:KME) :: NSAGG ! SELF-COLLECTION OF SNOW
848 REAL, DIMENSION(KMS:KME) :: NPRC ! CHANGE NC AUTOCONVERSION DROPLETS
849 REAL, DIMENSION(KMS:KME) :: NPRC1 ! CHANGE NR AUTOCONVERSION DROPLETS
850 REAL, DIMENSION(KMS:KME) :: PRAI ! CHANGE Q AUTOCONVERSION CLOUD ICE
851 REAL, DIMENSION(KMS:KME) :: PRCI ! CHANGE Q ACCRETION CLOUD ICE BY SNOW
852 REAL, DIMENSION(KMS:KME) :: PSACWS ! CHANGE Q DROPLET ACCRETION BY SNOW
853 REAL, DIMENSION(KMS:KME) :: NPSACWS ! CHANGE N DROPLET ACCRETION BY SNOW
854 REAL, DIMENSION(KMS:KME) :: PSACWI ! CHANGE Q DROPLET ACCRETION BY CLOUD ICE
855 REAL, DIMENSION(KMS:KME) :: NPSACWI ! CHANGE N DROPLET ACCRETION BY CLOUD ICE
856 REAL, DIMENSION(KMS:KME) :: NPRCI ! CHANGE N AUTOCONVERSION CLOUD ICE BY SNOW
857 REAL, DIMENSION(KMS:KME) :: NPRAI ! CHANGE N ACCRETION CLOUD ICE
858 REAL, DIMENSION(KMS:KME) :: NMULTS ! ICE MULT DUE TO RIMING DROPLETS BY SNOW
859 REAL, DIMENSION(KMS:KME) :: NMULTR ! ICE MULT DUE TO RIMING RAIN BY SNOW
860 REAL, DIMENSION(KMS:KME) :: QMULTS ! CHANGE Q DUE TO ICE MULT DROPLETS/SNOW
861 REAL, DIMENSION(KMS:KME) :: QMULTR ! CHANGE Q DUE TO ICE RAIN/SNOW
862 REAL, DIMENSION(KMS:KME) :: PRACS ! CHANGE Q RAIN-SNOW COLLECTION
863 REAL, DIMENSION(KMS:KME) :: NPRACS ! CHANGE N RAIN-SNOW COLLECTION
864 REAL, DIMENSION(KMS:KME) :: PCCN ! CHANGE Q DROPLET ACTIVATION
865 REAL, DIMENSION(KMS:KME) :: PSMLT ! CHANGE Q MELTING SNOW TO RAIN
866 REAL, DIMENSION(KMS:KME) :: EVPMS ! CHNAGE Q MELTING SNOW EVAPORATING
867 REAL, DIMENSION(KMS:KME) :: NSMLTS ! CHANGE N MELTING SNOW
868 REAL, DIMENSION(KMS:KME) :: NSMLTR ! CHANGE N MELTING SNOW TO RAIN
869 ! HM ADDED 12/13/06
870 REAL, DIMENSION(KMS:KME) :: PIACR ! CHANGE QR, ICE-RAIN COLLECTION
871 REAL, DIMENSION(KMS:KME) :: NIACR ! CHANGE N, ICE-RAIN COLLECTION
872 REAL, DIMENSION(KMS:KME) :: PRACI ! CHANGE QI, ICE-RAIN COLLECTION
873 REAL, DIMENSION(KMS:KME) :: PIACRS ! CHANGE QR, ICE RAIN COLLISION, ADDED TO SNOW
874 REAL, DIMENSION(KMS:KME) :: NIACRS ! CHANGE N, ICE RAIN COLLISION, ADDED TO SNOW
875 REAL, DIMENSION(KMS:KME) :: PRACIS ! CHANGE QI, ICE RAIN COLLISION, ADDED TO SNOW
876 REAL, DIMENSION(KMS:KME) :: EPRD ! SUBLIMATION CLOUD ICE
877 REAL, DIMENSION(KMS:KME) :: EPRDS ! SUBLIMATION SNOW
878 ! HM ADDED GRAUPEL PROCESSES
879 REAL, DIMENSION(KMS:KME) :: PRACG ! CHANGE IN Q COLLECTION RAIN BY GRAUPEL
880 REAL, DIMENSION(KMS:KME) :: PSACWG ! CHANGE IN Q COLLECTION DROPLETS BY GRAUPEL
881 REAL, DIMENSION(KMS:KME) :: PGSACW ! CONVERSION Q TO GRAUPEL DUE TO COLLECTION DROPLETS BY SNOW
882 REAL, DIMENSION(KMS:KME) :: PGRACS ! CONVERSION Q TO GRAUPEL DUE TO COLLECTION RAIN BY SNOW
883 REAL, DIMENSION(KMS:KME) :: PRDG ! DEP OF GRAUPEL
884 REAL, DIMENSION(KMS:KME) :: EPRDG ! SUB OF GRAUPEL
885 REAL, DIMENSION(KMS:KME) :: EVPMG ! CHANGE Q MELTING OF GRAUPEL AND EVAPORATION
886 REAL, DIMENSION(KMS:KME) :: PGMLT ! CHANGE Q MELTING OF GRAUPEL
887 REAL, DIMENSION(KMS:KME) :: NPRACG ! CHANGE N COLLECTION RAIN BY GRAUPEL
888 REAL, DIMENSION(KMS:KME) :: NPSACWG ! CHANGE N COLLECTION DROPLETS BY GRAUPEL
889 REAL, DIMENSION(KMS:KME) :: NSCNG ! CHANGE N CONVERSION TO GRAUPEL DUE TO COLLECTION DROPLETS BY SNOW
890 REAL, DIMENSION(KMS:KME) :: NGRACS ! CHANGE N CONVERSION TO GRAUPEL DUE TO COLLECTION RAIN BY SNOW
891 REAL, DIMENSION(KMS:KME) :: NGMLTG ! CHANGE N MELTING GRAUPEL
892 REAL, DIMENSION(KMS:KME) :: NGMLTR ! CHANGE N MELTING GRAUPEL TO RAIN
893 REAL, DIMENSION(KMS:KME) :: NSUBG ! CHANGE N SUB/DEP OF GRAUPEL
894 REAL, DIMENSION(KMS:KME) :: PSACR ! CONVERSION DUE TO COLL OF SNOW BY RAIN
895 REAL, DIMENSION(KMS:KME) :: NMULTG ! ICE MULT DUE TO ACC DROPLETS BY GRAUPEL
896 REAL, DIMENSION(KMS:KME) :: NMULTRG ! ICE MULT DUE TO ACC RAIN BY GRAUPEL
897 REAL, DIMENSION(KMS:KME) :: QMULTG ! CHANGE Q DUE TO ICE MULT DROPLETS/GRAUPEL
898 REAL, DIMENSION(KMS:KME) :: QMULTRG ! CHANGE Q DUE TO ICE MULT RAIN/GRAUPEL
899
900 ! TIME-VARYING ATMOSPHERIC PARAMETERS
901
902 REAL, DIMENSION(KMS:KME) :: KAP ! THERMAL CONDUCTIVITY OF AIR
903 REAL, DIMENSION(KMS:KME) :: EVS ! SATURATION VAPOR PRESSURE
904 REAL, DIMENSION(KMS:KME) :: EIS ! ICE SATURATION VAPOR PRESSURE
905 REAL, DIMENSION(KMS:KME) :: QVS ! SATURATION MIXING RATIO
906 REAL, DIMENSION(KMS:KME) :: QVI ! ICE SATURATION MIXING RATIO
907 REAL, DIMENSION(KMS:KME) :: QVQVS ! SAUTRATION RATIO
908 REAL, DIMENSION(KMS:KME) :: QVQVSI! ICE SATURAION RATIO
909 REAL, DIMENSION(KMS:KME) :: DV ! DIFFUSIVITY OF WATER VAPOR IN AIR
910 REAL, DIMENSION(KMS:KME) :: XXLS ! LATENT HEAT OF SUBLIMATION
911 REAL, DIMENSION(KMS:KME) :: XXLV ! LATENT HEAT OF VAPORIZATION
912 REAL, DIMENSION(KMS:KME) :: CPM ! SPECIFIC HEAT AT CONST PRESSURE FOR MOIST AIR
913 REAL, DIMENSION(KMS:KME) :: MU ! VISCOCITY OF AIR
914 REAL, DIMENSION(KMS:KME) :: SC ! SCHMIDT NUMBER
915 REAL, DIMENSION(KMS:KME) :: XLF ! LATENT HEAT OF FREEZING
916 REAL, DIMENSION(KMS:KME) :: RHO ! AIR DENSITY
917 REAL, DIMENSION(KMS:KME) :: AB ! CORRECTION TO CONDENSATION RATE DUE TO LATENT HEATING
918 REAL, DIMENSION(KMS:KME) :: ABI ! CORRECTION TO DEPOSITION RATE DUE TO LATENT HEATING
919
920 ! TIME-VARYING MICROPHYSICS PARAMETERS
921
922 REAL, DIMENSION(KMS:KME) :: 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
926
927 ! FALL SPEED WORKING VARIABLES (DEFINED IN CODE)
928
929 REAL, DIMENSION(KMS:KME) :: DUMI,DUMR,DUMFNI,DUMG,DUMFNG
930 REAL UNI, UMI,UMR
931 REAL, DIMENSION(KMS:KME) :: FR, FI, FNI,FG,FNG
932 REAL RGVM
933 REAL, DIMENSION(KMS:KME) :: FALOUTR,FALOUTI,FALOUTNI
934 REAL FALTNDR,FALTNDI,FALTNDNI,RHO2
935 REAL, DIMENSION(KMS:KME) :: DUMQS,DUMFNS
936 REAL UMS,UNS
937 REAL, DIMENSION(KMS:KME) :: FS,FNS, FALOUTS,FALOUTNS,FALOUTG,FALOUTNG
938 REAL FALTNDS,FALTNDNS,UNR,FALTNDG,FALTNDNG
939 REAL, DIMENSION(KMS:KME) :: DUMC,DUMFNC
940 REAL UNC,UMC,UNG,UMG
941 REAL, DIMENSION(KMS:KME) :: FC,FALOUTC,FALOUTNC
942 REAL FALTNDC,FALTNDNC
943 REAL, DIMENSION(KMS:KME) :: FNC,DUMFNR,FALOUTNR
944 REAL FALTNDNR
945 REAL, DIMENSION(KMS:KME) :: FNR(KMS:KME)
946
947 ! FALL-SPEED PARAMETER 'A' WITH AIR DENSITY CORRECTION
948
949 REAL, DIMENSION(KMS:KME) :: AIN,ARN,ASN,ACN,AGN
950
951 ! EXTERNAL FUNCTION CALL RETURN VARIABLES
952
953 ! REAL GAMMA, ! EULER GAMMA FUNCTION
954 ! REAL POLYSVP, ! SAT. PRESSURE FUNCTION
955 ! REAL DERF1 ! ERROR FUNCTION
956
957 ! DUMMY VARIABLES
958
959 REAL DUM,DUM1,DUM2,DUMT,DUMQV,DUMQSS,DUMQSI,DUMS
960
961 ! PROGNOSTIC SUPERSATURATION
962
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
969
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
977
978 ! COUNTING/INDEX VARIABLES
979
980 INTEGER K,NSTEP,N ! ,I
981
982 ! LTRUE IS ONLY USED TO SPEED UP THE CODE !!
983 ! LTRUE, SWITCH = 0, NO HYDROMETEORS IN COLUMN,
984 ! = 1, HYDROMETEORS IN COLUMN
985
986 INTEGER LTRUE
987
988 ! DROPLET ACTIVATION/FREEZING AEROSOL
989
990
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
997
998 REAL CRY,KRY ! AEROSOL ACTIVATION PARAMETERS
999
1000 ! MORE WORKING/DUMMY VARIABLES
1001
1002 REAL DUMQI,DUMNI,DC0,DS0,DG0
1003 REAL DUMQC,DUMQR,RATIO,SUM_DEP,FUDGEF
1004
1005 ! EFFECTIVE VERTICAL VELOCITY (M/S)
1006 REAL WEF
1007
1008 ! WORKING PARAMETERS FOR ICE NUCLEATION
1009
1010 REAL ANUC,BNUC
1011
1012 ! WORKING PARAMETERS FOR AEROSOL ACTIVATION
1013
1014 REAL AACT,GAMM,GG,PSI,ETA1,ETA2,SM1,SM2,SMAX,UU1,UU2,ALPHA
1015
1016 ! DUMMY SIZE DISTRIBUTION PARAMETERS
1017
1018 REAL DLAMS,DLAMR,DLAMI,DLAMC,DLAMG,LAMMAX,LAMMIN
1019
1020 INTEGER IDROP
1021
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(KMS:KME) :: NC3DTEN ! CLOUD DROPLET NUMBER CONCENTRATION (1/KG/S)
1026 REAL, DIMENSION(KMS:KME) :: NC3D ! CLOUD DROPLET NUMBER CONCENTRATION (1/KG)
1027
1028 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1029
1030 ! SET LTRUE INITIALLY TO 0
1031
1032 LTRUE = 0
1033
1034 ! ATMOSPHERIC PARAMETERS THAT VARY IN TIME AND HEIGHT
1035 DO K = KTS,KTE
1036
1037 ! LATENT HEAT OF VAPORATION
1038
1039 XXLV(K) = 3.1484E6-2370.*T3D(K)
1040
1041 ! LATENT HEAT OF SUBLIMATION
1042
1043 XXLS(K) = 3.15E6-2370.*T3D(K)+0.3337E6
1044
1045 CPM(K) = CP*(1.+0.887*QV3D(K))
1046
1047 ! SATURATION VAPOR PRESSURE AND MIXING RATIO
1048
1049 EVS(K) = POLYSVP(T3D(K),0) ! PA
1050 EIS(K) = POLYSVP(T3D(K),1) ! PA
1051
1052 ! MAKE SURE ICE SATURATION DOESN'T EXCEED WATER SAT. NEAR FREEZING
1053
1054 IF (EIS(K).GT.EVS(K)) EIS(K) = EVS(K)
1055
1056 QVS(K) = .622*EVS(K)/(PRES(K)-EVS(K))
1057 QVI(K) = .622*EIS(K)/(PRES(K)-EIS(K))
1058
1059 QVQVS(K) = QV3D(K)/QVS(K)
1060 QVQVSI(K) = QV3D(K)/QVI(K)
1061
1062 ! AIR DENSITY
1063
1064 RHO(K) = PRES(K)/(R*T3D(K))
1065
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
1070
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
1083
1084 ! AT SUBSATURATION, REMOVE SMALL AMOUNTS OF CLOUD/PRECIP WATER
1085
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
1098
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
1116
1117 ! HEAT OF FUSION
1118
1119 XLF(K) = XXLS(K)-XXLV(K)
1120
1121 !..................................................................
1122 ! IF MIXING RATIO < QSMALL SET MIXING RATIO AND NUMBER CONC TO ZERO
1123
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
1149
1150 ! INITIALIZE SEDIMENTATION TENDENCIES FOR MIXING RATIO
1151
1152 QRSTEN(K) = 0.
1153 QISTEN(K) = 0.
1154 QNISTEN(K) = 0.
1155 QCSTEN(K) = 0.
1156 QGSTEN(K) = 0.
1157
1158 !..................................................................
1159 ! MICROPHYSICS PARAMETERS VARYING IN TIME/HEIGHT
1160
1161 ! FALL SPEED WITH DENSITY CORRECTION (HEYMSFIELD AND BENSSEMER 2006)
1162
1163 DUM = (RHOSU/RHO(K))**0.54
1164
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
1171
1172 !..................................
1173 ! IF THERE IS NO CLOUD/PRECIP WATER, AND IF SUBSATURATED, THEN SKIP MICROPHYSICS
1174 ! FOR THIS LEVEL
1175
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
1181
1182 ! THERMAL CONDUCTIVITY FOR AIR
1183
1184 DUM = 1.496E-6*T3D(K)**1.5/(T3D(K)+120.)
1185
1186 ! KAP(K) = 1.414E3*1.496E-6*T3D(K)**1.5/(T3D(K)+120.)
1187 KAP(K) = 1.414E3*DUM
1188
1189 ! DIFFUSIVITY OF WATER VAPOR
1190
1191 DV(K) = 8.794E-5*T3D(K)**1.81/PRES(K)
1192
1193 ! VISCOSITY OF AIR
1194 ! SCHMIT NUMBER
1195
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)
1199
1200 ! PSYCHOMETIC CORRECTIONS
1201
1202 ! RATE OF CHANGE SAT. MIX. RATIO WITH TEMPERATURE
1203
1204 DUM = (RV*T3D(K)**2)
1205
1206 DQSDT = XXLV(K)*QVS(K)/DUM
1207 DQSIDT = XXLS(K)*QVI(K)/DUM
1208
1209 ABI(K) = 1.+DQSIDT*XXLS(K)/CPM(K)
1210 AB(K) = 1.+DQSDT*XXLV(K)/CPM(K)
1211
1212 !
1213 !.....................................................................
1214 !.....................................................................
1215 ! CASE FOR TEMPERATURE ABOVE FREEZING
1216
1217 IF (T3D(K).GE.273.15) THEN
1218
1219 !......................................................................
1220 !HM ADD, ALLOW FOR CONSTANT DROPLET NUMBER
1221 ! INUM = 0, PREDICT DROPLET NUMBER
1222 ! INUM = 1, SET CONSTANT DROPLET NUMBER
1223
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
1228
1229 ! GET SIZE DISTRIBUTION PARAMETERS
1230
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
1246
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
1248
1249 ! MAKE SURE NUMBER CONCENTRATIONS AREN'T NEGATIVE
1250
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))
1255
1256 !......................................................................
1257 ! RAIN
1258
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)
1262
1263 ! CHECK FOR SLOPE
1264
1265 ! ADJUST VARS
1266
1267 IF (LAMR(K).LT.LAMMINR) THEN
1268
1269 LAMR(K) = LAMMINR
1270
1271 N0RR(K) = LAMR(K)**4*QR3D(K)/(PI*RHOW)
1272
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)
1277
1278 NR3D(K) = N0RR(K)/LAMR(K)
1279 END IF
1280 END IF
1281
1282 !......................................................................
1283 ! CLOUD DROPLETS
1284
1285 ! MARTIN ET AL. (1994) FORMULA FOR PGAM
1286
1287 IF (QC3D(K).GE.QSMALL) THEN
1288
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.)
1294
1295 ! CALCULATE LAMC
1296
1297 LAMC(K) = (CONS26*NC3D(K)*GAMMA(PGAM(K)+4.)/ &
1298 (QC3D(K)*GAMMA(PGAM(K)+1.)))**(1./3.)
1299
1300 ! LAMMIN, 60 MICRON DIAMETER
1301 ! LAMMAX, 1 MICRON
1302
1303 LAMMIN = (PGAM(K)+1.)/60.E-6
1304 LAMMAX = (PGAM(K)+1.)/1.E-6
1305
1306 IF (LAMC(K).LT.LAMMIN) THEN
1307 LAMC(K) = LAMMIN
1308
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
1313
1314 NC3D(K) = EXP(3.*LOG(LAMC(K))+LOG(QC3D(K))+ &
1315 LOG(GAMMA(PGAM(K)+1.))-LOG(GAMMA(PGAM(K)+4.)))/CONS26
1316
1317 END IF
1318
1319 END IF
1320
1321 !......................................................................
1322 ! SNOW
1323
1324 IF (QNI3D(K).GE.QSMALL) THEN
1325 LAMS(K) = (CONS1*NS3D(K)/QNI3D(K))**(1./DS)
1326 N0S(K) = NS3D(K)*LAMS(K)
1327
1328 ! CHECK FOR SLOPE
1329
1330 ! ADJUST VARS
1331
1332 IF (LAMS(K).LT.LAMMINS) THEN
1333 LAMS(K) = LAMMINS
1334 N0S(K) = LAMS(K)**4*QNI3D(K)/CONS1
1335
1336 NS3D(K) = N0S(K)/LAMS(K)
1337
1338 ELSE IF (LAMS(K).GT.LAMMAXS) THEN
1339
1340 LAMS(K) = LAMMAXS
1341 N0S(K) = LAMS(K)**4*QNI3D(K)/CONS1
1342
1343 NS3D(K) = N0S(K)/LAMS(K)
1344 END IF
1345 END IF
1346
1347 !......................................................................
1348 ! GRAUPEL
1349
1350 IF (QG3D(K).GE.QSMALL) THEN
1351 LAMG(K) = (CONS2*NG3D(K)/QG3D(K))**(1./DG)
1352 N0G(K) = NG3D(K)*LAMG(K)
1353
1354 ! ADJUST VARS
1355
1356 IF (LAMG(K).LT.LAMMING) THEN
1357 LAMG(K) = LAMMING
1358 N0G(K) = LAMG(K)**4*QG3D(K)/CONS2
1359
1360 NG3D(K) = N0G(K)/LAMG(K)
1361
1362 ELSE IF (LAMG(K).GT.LAMMAXG) THEN
1363
1364 LAMG(K) = LAMMAXG
1365 N0G(K) = LAMG(K)**4*QG3D(K)/CONS2
1366
1367 NG3D(K) = N0G(K)/LAMG(K)
1368 END IF
1369 END IF
1370
1371 !.....................................................................
1372 ! ZERO OUT PROCESS RATES
1373
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.
1398
1399 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1400 ! CALCULATION OF MICROPHYSICAL PROCESS RATES, T > 273.15 K
1401
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
1409
1410 ! USE MINIMUM VALUE OF 1.E-6 TO PREVENT FLOATING POINT ERROR
1411
1412 IF (QC3D(K).GE.1.E-6) THEN
1413
1414 ! HM ADD 12/13/06, REPLACE WITH NEWER FORMULA
1415 ! FROM KHAIROUTDINOV AND KOGAN 2000, MWR
1416
1417 PRC(K)=1350.*QC3D(K)**2.47* &
1418 (NC3D(K)/1.e6*RHO(K))**(-1.79)
1419
1420 ! note: nprc1 is change in Nr,
1421 ! nprc is change in Nc
1422
1423 NPRC1(K) = PRC(K)/CONS29
1424 NPRC(K) = PRC(K)/(QC3D(k)/NC3D(K))
1425
1426 NPRC(K) = MIN(NPRC(K),NC3D(K)/DT)
1427
1428 END IF
1429
1430 !.......................................................................
1431 ! HM ADD 12/13/06, COLLECTION OF SNOW BY RAIN ABOVE FREEZING
1432 ! FORMULA FROM IKAWA AND SAITO (1991)
1433
1434 IF (QR3D(K).GE.1.E-8.AND.QNI3D(K).GE.1.E-8) THEN
1435
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
1440
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)
1446
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)))
1453
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))
1459
1460 END IF
1461
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
1465
1466 IF (QR3D(K).GE.1.E-8.AND.QG3D(K).GE.1.E-8) THEN
1467
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
1472
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)
1478
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)))
1486
1487 ! ASSUME 1 MM DROPS ARE SHED, GET NUMBER SHED PER SEC
1488
1489 DUM = DUM/5.2E-7
1490
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))
1496
1497 NPRACG(K)=MAX(NPRACG(K)-DUM,0.)
1498
1499 END IF
1500
1501 !.......................................................................
1502 ! ACCRETION OF CLOUD LIQUID WATER BY RAIN
1503 ! CONTINUOUS COLLECTION EQUATION WITH
1504 ! GRAVITATIONAL COLLECTION KERNEL, DROPLET FALL SPEED NEGLECTED
1505
1506 IF (QR3D(K).GE.1.E-8 .AND. QC3D(K).GE.1.E-8) THEN
1507
1508 ! 12/13/06 HM ADD, REPLACE WITH NEWER FORMULA FROM
1509 ! KHAIROUTDINOV AND KOGAN 2000, MWR
1510
1511 DUM=(QC3D(K)*QR3D(K))
1512 PRA(K) = 67.*(DUM)**1.15
1513 NPRA(K) = PRA(K)/(QC3D(K)/NC3D(K))
1514
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
1521
1522 IF (QR3D(K).GE.1.E-8) THEN
1523 NRAGG(K) = -8.*NR3D(K)*QR3D(K)*RHO(K)
1524 END IF
1525
1526 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1527 ! CALCULATE EVAP OF RAIN (RUTLEDGE AND HOBBS 1983)
1528
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
1538
1539 ! NO CONDENSATION ONTO RAIN, ONLY EVAP ALLOWED
1540
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
1547
1548 !.......................................................................
1549 ! MELTING OF SNOW
1550
1551 ! SNOW MAY PERSITS ABOVE FREEZING, FORMULA FROM RUTLEDGE AND HOBBS, 1984
1552 ! IF WATER SUPERSATURATION, SNOW MELTS TO FORM RAIN
1553
1554 IF (QNI3D(K).GE.1.E-8) THEN
1555
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))
1561
1562 ! IN WATER SUBSATURATION, SNOW MELTS AND EVAPORATES
1563
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
1575
1576 !.......................................................................
1577 ! MELTING OF GRAUPEL
1578
1579 ! GRAUPEL MAY PERSITS ABOVE FREEZING, FORMULA FROM RUTLEDGE AND HOBBS, 1984
1580 ! IF WATER SUPERSATURATION, GRAUPEL MELTS TO FORM RAIN
1581
1582 IF (QG3D(K).GE.1.E-8) THEN
1583
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))
1589
1590 ! IN WATER SUBSATURATION, GRAUPEL MELTS AND EVAPORATES
1591
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
1603
1604 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1605
1606 ! FOR CLOUD ICE, ONLY PROCESSES OPERATING AT T > 273.15 IS
1607 ! MELTING, WHICH IS ALREADY CONSERVED DURING PROCESS
1608 ! CALCULATION
1609
1610 ! CONSERVATION OF QC
1611
1612 DUM = (PRC(K)+PRA(K))*DT
1613
1614 IF (DUM.GT.QC3D(K).AND.QC3D(K).GE.QSMALL) THEN
1615
1616 RATIO = QC3D(K)/DUM
1617
1618 PRC(K) = PRC(K)*RATIO
1619 PRA(K) = PRA(K)*RATIO
1620
1621 END IF
1622
1623 ! CONSERVATION OF SNOW
1624
1625 DUM = (-PSMLT(K)-EVPMS(K)+PRACS(K))*DT
1626
1627 IF (DUM.GT.QNI3D(K).AND.QNI3D(K).GE.QSMALL) THEN
1628
1629 ! NO SOURCE TERMS FOR SNOW AT T > FREEZING
1630 RATIO = QNI3D(K)/DUM
1631
1632 PSMLT(K) = PSMLT(K)*RATIO
1633 EVPMS(K) = EVPMS(K)*RATIO
1634 PRACS(K) = PRACS(K)*RATIO
1635
1636 END IF
1637
1638 ! CONSERVATION OF GRAUPEL
1639
1640 DUM = (-PGMLT(K)-EVPMG(K)+PRACG(K))*DT
1641
1642 IF (DUM.GT.QG3D(K).AND.QG3D(K).GE.QSMALL) THEN
1643
1644 ! NO SOURCE TERM FOR GRAUPEL ABOVE FREEZING
1645 RATIO = QG3D(K)/DUM
1646
1647 PGMLT(K) = PGMLT(K)*RATIO
1648 EVPMG(K) = EVPMG(K)*RATIO
1649 PRACG(K) = PRACG(K)*RATIO
1650
1651 END IF
1652
1653 ! CONSERVATION OF QR
1654 ! HM 12/13/06, ADDED CONSERVATION OF RAIN SINCE PRE IS NEGATIVE
1655
1656 DUM = (-PRACS(K)-PRACG(K)-PRE(K)-PRA(K)-PRC(K)+PSMLT(K)+PGMLT(K))*DT
1657
1658 IF (DUM.GT.QR3D(K).AND.QR3D(K).GE.QSMALL) THEN
1659
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
1663
1664 END IF
1665
1666 !....................................
1667
1668 QV3DTEN(K) = QV3DTEN(K)+(-PRE(K)-EVPMS(K)-EVPMG(K))
1669
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)
1672
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 NG3DTEN(K) = NG3DTEN(K)-NPRACG(K)
1679 NC3DTEN(K) = NC3DTEN(K)+ (-NPRA(K)-NPRC(K))
1680 NR3DTEN(K) = NR3DTEN(K)+ (NPRC1(K)+NRAGG(K))
1681
1682 IF (PRE(K).LT.0.) THEN
1683 DUM = PRE(K)*DT/QR3D(K)
1684 DUM = MAX(-1.,DUM)
1685 NSUBR(K) = DUM*NR3D(K)/DT
1686 END IF
1687
1688 IF (EVPMS(K)+PSMLT(K).LT.0.) THEN
1689 DUM = (EVPMS(K)+PSMLT(K))*DT/QNI3D(K)
1690 DUM = MAX(-1.,DUM)
1691 NSMLTS(K) = DUM*NS3D(K)/DT
1692 END IF
1693 IF (PSMLT(K).LT.0.) THEN
1694 DUM = PSMLT(K)*DT/QNI3D(K)
1695 DUM = MAX(-1.0,DUM)
1696 NSMLTR(K) = DUM*NS3D(K)/DT
1697 END IF
1698 IF (EVPMG(K)+PGMLT(K).LT.0.) THEN
1699 DUM = (EVPMG(K)+PGMLT(K))*DT/QG3D(K)
1700 DUM = MAX(-1.,DUM)
1701 NGMLTG(K) = DUM*NG3D(K)/DT
1702 END IF
1703 IF (PGMLT(K).LT.0.) THEN
1704 DUM = PGMLT(K)*DT/QG3D(K)
1705 DUM = MAX(-1.0,DUM)
1706 NGMLTR(K) = DUM*NG3D(K)/DT
1707 END IF
1708
1709 NS3DTEN(K) = NS3DTEN(K)+(NSMLTS(K))
1710 NG3DTEN(K) = NG3DTEN(K)+(NGMLTG(K))
1711 NR3DTEN(K) = NR3DTEN(K)+(NSUBR(K)-NSMLTR(K)-NGMLTR(K))
1712
1713 300 CONTINUE
1714
1715 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1716 ! NOW CALCULATE SATURATION ADJUSTMENT TO CONDENSE EXTRA VAPOR ABOVE
1717 ! WATER SATURATION
1718
1719 DUMT = T3D(K)+DT*T3DTEN(K)
1720 DUMQV = QV3D(K)+DT*QV3DTEN(K)
1721 DUMQSS = 0.622*POLYSVP(DUMT,0)/ (PRES(K)-POLYSVP(DUMT,0))
1722 DUMQC = QC3D(K)+DT*QC3DTEN(K)
1723 DUMQC = MAX(DUMQC,0.)
1724
1725 ! SATURATION ADJUSTMENT FOR LIQUID
1726
1727 DUMS = DUMQV-DUMQSS
1728 PCC(K) = DUMS/(1.+XXLV(K)**2*DUMQSS/(CPM(K)*RV*DUMT**2))/DT
1729 IF (PCC(K)*DT+DUMQC.LT.0.) THEN
1730 PCC(K) = -DUMQC/DT
1731 END IF
1732
1733 QV3DTEN(K) = QV3DTEN(K)-PCC(K)
1734 T3DTEN(K) = T3DTEN(K)+PCC(K)*XXLV(K)/CPM(K)
1735 QC3DTEN(K) = QC3DTEN(K)+PCC(K)
1736
1737 !.......................................................................
1738 ! ACTIVATION OF CLOUD DROPLETS
1739 ! ACTIVATION OF DROPLET CURRENTLY NOT CALCULATED
1740 ! DROPLET CONCENTRATION IS SPECIFIED !!!!!
1741
1742 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1743 ! SUBLIMATE, MELT, OR EVAPORATE NUMBER CONCENTRATION
1744 ! THIS FORMULATION ASSUMES 1:1 RATIO BETWEEN MASS LOSS AND
1745 ! LOSS OF NUMBER CONCENTRATION
1746
1747 ! IF (PCC(K).LT.0.) THEN
1748 ! DUM = PCC(K)*DT/QC3D(K)
1749 ! DUM = MAX(-1.,DUM)
1750 ! NSUBC(K) = DUM*NC3D(K)/DT
1751 ! END IF
1752
1753 ! UPDATE TENDENCIES
1754
1755 ! NC3DTEN(K) = NC3DTEN(K)+NSUBC(K)
1756
1757 !.....................................................................
1758 !.....................................................................
1759 ELSE ! TEMPERATURE < 273.15
1760
1761 !......................................................................
1762 !HM ADD, ALLOW FOR CONSTANT DROPLET NUMBER
1763 ! INUM = 0, PREDICT DROPLET NUMBER
1764 ! INUM = 1, SET CONSTANT DROPLET NUMBER
1765
1766 ! IF (INUM.EQ.1) THEN
1767 ! CONVERT NDCNST FROM CM-3 TO KG-1
1768 NC3D(K)=NDCNST*1.E6/RHO(K)
1769 ! END IF
1770
1771 ! CALCULATE SIZE DISTRIBUTION PARAMETERS
1772 ! MAKE SURE NUMBER CONCENTRATIONS AREN'T NEGATIVE
1773
1774 NI3D(K) = MAX(0.,NI3D(K))
1775 NS3D(K) = MAX(0.,NS3D(K))
1776 NC3D(K) = MAX(0.,NC3D(K))
1777 NR3D(K) = MAX(0.,NR3D(K))
1778 NG3D(K) = MAX(0.,NG3D(K))
1779
1780 !......................................................................
1781 ! CLOUD ICE
1782
1783 IF (QI3D(K).GE.QSMALL) THEN
1784 LAMI(K) = (CONS12* &
1785 NI3D(K)/QI3D(K))**(1./DI)
1786 N0I(K) = NI3D(K)*LAMI(K)
1787
1788 ! CHECK FOR SLOPE
1789
1790 ! ADJUST VARS
1791
1792 IF (LAMI(K).LT.LAMMINI) THEN
1793
1794 LAMI(K) = LAMMINI
1795
1796 N0I(K) = LAMI(K)**4*QI3D(K)/CONS12
1797
1798 NI3D(K) = N0I(K)/LAMI(K)
1799 ELSE IF (LAMI(K).GT.LAMMAXI) THEN
1800 LAMI(K) = LAMMAXI
1801 N0I(K) = LAMI(K)**4*QI3D(K)/CONS12
1802
1803 NI3D(K) = N0I(K)/LAMI(K)
1804 END IF
1805 END IF
1806
1807 !......................................................................
1808 ! RAIN
1809
1810 IF (QR3D(K).GE.QSMALL) THEN
1811 LAMR(K) = (PI*RHOW*NR3D(K)/QR3D(K))**(1./3.)
1812 N0RR(K) = NR3D(K)*LAMR(K)
1813
1814 ! CHECK FOR SLOPE
1815
1816 ! ADJUST VARS
1817
1818 IF (LAMR(K).LT.LAMMINR) THEN
1819
1820 LAMR(K) = LAMMINR
1821
1822 N0RR(K) = LAMR(K)**4*QR3D(K)/(PI*RHOW)
1823
1824 NR3D(K) = N0RR(K)/LAMR(K)
1825 ELSE IF (LAMR(K).GT.LAMMAXR) THEN
1826 LAMR(K) = LAMMAXR
1827 N0RR(K) = LAMR(K)**4*QR3D(K)/(PI*RHOW)
1828
1829 NR3D(K) = N0RR(K)/LAMR(K)
1830 END IF
1831 END IF
1832
1833 !......................................................................
1834 ! CLOUD DROPLETS
1835
1836 ! MARTIN ET AL. (1994) FORMULA FOR PGAM
1837
1838 IF (QC3D(K).GE.QSMALL) THEN
1839
1840 DUM = PRES(K)/(287.15*T3D(K))
1841 PGAM(K)=0.0005714*(NC3D(K)/1.E6/DUM)+0.2714
1842 PGAM(K)=1./(PGAM(K)**2)-1.
1843 PGAM(K)=MAX(PGAM(K),2.)
1844 PGAM(K)=MIN(PGAM(K),10.)
1845
1846 ! CALCULATE LAMC
1847
1848 LAMC(K) = (CONS26*NC3D(K)*GAMMA(PGAM(K)+4.)/ &
1849 (QC3D(K)*GAMMA(PGAM(K)+1.)))**(1./3.)
1850
1851 ! LAMMIN, 60 MICRON DIAMETER
1852 ! LAMMAX, 1 MICRON
1853
1854 LAMMIN = (PGAM(K)+1.)/60.E-6
1855 LAMMAX = (PGAM(K)+1.)/1.E-6
1856
1857 IF (LAMC(K).LT.LAMMIN) THEN
1858 LAMC(K) = LAMMIN
1859
1860 NC3D(K) = EXP(3.*LOG(LAMC(K))+LOG(QC3D(K))+ &
1861 LOG(GAMMA(PGAM(K)+1.))-LOG(GAMMA(PGAM(K)+4.)))/CONS26
1862 ELSE IF (LAMC(K).GT.LAMMAX) THEN
1863 LAMC(K) = LAMMAX
1864 NC3D(K) = EXP(3.*LOG(LAMC(K))+LOG(QC3D(K))+ &
1865 LOG(GAMMA(PGAM(K)+1.))-LOG(GAMMA(PGAM(K)+4.)))/CONS26
1866
1867 END IF
1868
1869 ! TO CALCULATE DROPLET FREEZING
1870
1871 CDIST1(K) = NC3D(K)/GAMMA(PGAM(K)+1.)
1872
1873 END IF
1874
1875 !......................................................................
1876 ! SNOW
1877
1878 IF (QNI3D(K).GE.QSMALL) THEN
1879 LAMS(K) = (CONS1*NS3D(K)/QNI3D(K))**(1./DS)
1880 N0S(K) = NS3D(K)*LAMS(K)
1881
1882 ! CHECK FOR SLOPE
1883
1884 ! ADJUST VARS
1885
1886 IF (LAMS(K).LT.LAMMINS) THEN
1887 LAMS(K) = LAMMINS
1888 N0S(K) = LAMS(K)**4*QNI3D(K)/CONS1
1889
1890 NS3D(K) = N0S(K)/LAMS(K)
1891
1892 ELSE IF (LAMS(K).GT.LAMMAXS) THEN
1893
1894 LAMS(K) = LAMMAXS
1895 N0S(K) = LAMS(K)**4*QNI3D(K)/CONS1
1896
1897 NS3D(K) = N0S(K)/LAMS(K)
1898 END IF
1899 END IF
1900
1901 !......................................................................
1902 ! GRAUPEL
1903
1904 IF (QG3D(K).GE.QSMALL) THEN
1905 LAMG(K) = (CONS2*NG3D(K)/QG3D(K))**(1./DG)
1906 N0G(K) = NG3D(K)*LAMG(K)
1907
1908 ! CHECK FOR SLOPE
1909
1910 ! ADJUST VARS
1911
1912 IF (LAMG(K).LT.LAMMING) THEN
1913 LAMG(K) = LAMMING
1914 N0G(K) = LAMG(K)**4*QG3D(K)/CONS2
1915
1916 NG3D(K) = N0G(K)/LAMG(K)
1917
1918 ELSE IF (LAMG(K).GT.LAMMAXG) THEN
1919
1920 LAMG(K) = LAMMAXG
1921 N0G(K) = LAMG(K)**4*QG3D(K)/CONS2
1922
1923 NG3D(K) = N0G(K)/LAMG(K)
1924 END IF
1925 END IF
1926
1927 !.....................................................................
1928 ! ZERO OUT PROCESS RATES
1929
1930 MNUCCC(K) = 0.
1931 NNUCCC(K) = 0.
1932 PRC(K) = 0.
1933 NPRC(K) = 0.
1934 NPRC1(K) = 0.
1935 NSAGG(K) = 0.
1936 PSACWS(K) = 0.
1937 NPSACWS(K) = 0.
1938 PSACWI(K) = 0.
1939 NPSACWI(K) = 0.
1940 PRACS(K) = 0.
1941 NPRACS(K) = 0.
1942 NMULTS(K) = 0.
1943 QMULTS(K) = 0.
1944 NMULTR(K) = 0.
1945 QMULTR(K) = 0.
1946 NMULTG(K) = 0.
1947 QMULTG(K) = 0.
1948 NMULTRG(K) = 0.
1949 QMULTRG(K) = 0.
1950 MNUCCR(K) = 0.
1951 NNUCCR(K) = 0.
1952 PRA(K) = 0.
1953 NPRA(K) = 0.
1954 NRAGG(K) = 0.
1955 PRCI(K) = 0.
1956 NPRCI(K) = 0.
1957 PRAI(K) = 0.
1958 NPRAI(K) = 0.
1959 NNUCCD(K) = 0.
1960 MNUCCD(K) = 0.
1961 PCC(K) = 0.
1962 PRE(K) = 0.
1963 PRD(K) = 0.
1964 PRDS(K) = 0.
1965 EPRD(K) = 0.
1966 EPRDS(K) = 0.
1967 NSUBC(K) = 0.
1968 NSUBI(K) = 0.
1969 NSUBS(K) = 0.
1970 NSUBR(K) = 0.
1971 PIACR(K) = 0.
1972 NIACR(K) = 0.
1973 PRACI(K) = 0.
1974 PIACRS(K) = 0.
1975 NIACRS(K) = 0.
1976 PRACIS(K) = 0.
1977 ! HM: ADD GRAUPEL PROCESSES
1978 PRACG(K) = 0.
1979 PSACR(K) = 0.
1980 PSACWG(K) = 0.
1981 PGSACW(K) = 0.
1982 PGRACS(K) = 0.
1983 PRDG(K) = 0.
1984 EPRDG(K) = 0.
1985 NPRACG(K) = 0.
1986 NPSACWG(K) = 0.
1987 NSCNG(K) = 0.
1988 NGRACS(K) = 0.
1989 NSUBG(K) = 0.
1990
1991 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1992 ! CALCULATION OF MICROPHYSICAL PROCESS RATES
1993 ! ACCRETION/AUTOCONVERSION/FREEZING/MELTING/COAG.
1994 !.......................................................................
1995 ! FREEZING OF CLOUD DROPLETS
1996 ! ONLY ALLOWED BELOW -4 C
1997 IF (QC3D(K).GE.QSMALL .AND. T3D(K).LT.269.15) THEN
1998
1999 ! NUMBER OF CONTACT NUCLEI (M^-3) FROM MEYERS ET AL., 1992
2000 ! FACTOR OF 1000 IS TO CONVERT FROM L^-1 TO M^-3
2001
2002 ! MEYERS CURVE
2003
2004 NACNT = EXP(-2.80+0.262*(273.15-T3D(K)))*1000.
2005
2006 ! COOPER CURVE
2007 ! NACNT = 5.*EXP(0.304*(273.15-T3D(K)))
2008
2009 ! FLECTHER
2010 ! NACNT = 0.01*EXP(0.6*(273.15-T3D(K)))
2011
2012 ! CONTACT FREEZING
2013
2014 ! MEAN FREE PATH
2015
2016 DUM = 7.37*T3D(K)/(288.*10.*PRES(K))/100.
2017
2018 ! EFFECTIVE DIFFUSIVITY OF CONTACT NUCLEI
2019 ! BASED ON BROWNIAN DIFFUSION
2020
2021 DAP(K) = CONS37*T3D(K)*(1.+DUM/RIN)/MU(K)
2022
2023 MNUCCC(K) = CONS38*DAP(K)*NACNT*EXP(LOG(CDIST1(K))+ &
2024 LOG(GAMMA(PGAM(K)+5.))-4.*LOG(LAMC(K)))
2025 NNUCCC(K) = 2.*PI*DAP(K)*NACNT*CDIST1(K)* &
2026 GAMMA(PGAM(K)+2.)/ &
2027 LAMC(K)
2028
2029 ! IMMERSION FREEZING (BIGG 1953)
2030
2031 MNUCCC(K) = MNUCCC(K)+CONS39* &
2032 EXP(LOG(CDIST1(K))+LOG(GAMMA(7.+PGAM(K)))-6.*LOG(LAMC(K)))* &
2033 EXP(AIMM*(273.15-T3D(K)))
2034
2035 NNUCCC(K) = NNUCCC(K)+ &
2036 CONS40*EXP(LOG(CDIST1(K))+LOG(GAMMA(PGAM(K)+4.))-3.*LOG(LAMC(K))) &
2037 *EXP(AIMM*(273.15-T3D(K)))
2038
2039 ! PUT IN A CATCH HERE TO PREVENT DIVERGENCE BETWEEN NUMBER CONC. AND
2040 ! MIXING RATIO, SINCE STRICT CONSERVATION NOT CHECKED FOR NUMBER CONC
2041
2042 NNUCCC(K) = MIN(NNUCCC(K),NC3D(K)/DT)
2043
2044 END IF
2045
2046 !.................................................................
2047 !.......................................................................
2048 ! AUTOCONVERSION OF CLOUD LIQUID WATER TO RAIN
2049 ! FORMULA FROM BEHENG (1994)
2050 ! USING NUMERICAL SIMULATION OF STOCHASTIC COLLECTION EQUATION
2051 ! AND INITIAL CLOUD DROPLET SIZE DISTRIBUTION SPECIFIED
2052 ! AS A GAMMA DISTRIBUTION
2053
2054 ! USE MINIMUM VALUE OF 1.E-6 TO PREVENT FLOATING POINT ERROR
2055
2056 IF (QC3D(K).GE.1.E-6) THEN
2057
2058 ! HM ADD 12/13/06, REPLACE WITH NEWER FORMULA
2059 ! FROM KHAIROUTDINOV AND KOGAN 2000, MWR
2060
2061 PRC(K)=1350.*QC3D(K)**2.47* &
2062 (NC3D(K)/1.e6*RHO(K))**(-1.79)
2063
2064 ! note: nprc1 is change in Nr,
2065 ! nprc is change in Nc
2066
2067 NPRC1(K) = PRC(K)/CONS29
2068 NPRC(K) = PRC(K)/(QC3D(K)/NC3D(K))
2069
2070 NPRC(K) = MIN(NPRC(K),NC3D(K)/DT)
2071
2072 END IF
2073
2074 !.......................................................................
2075 ! SELF-COLLECTION OF DROPLET NOT INCLUDED IN KK2000 SCHEME
2076
2077 ! SNOW AGGREGATION FROM PASSARELLI, 1978, USED BY REISNER, 1998
2078 ! THIS IS HARD-WIRED FOR BS = 0.4 FOR NOW
2079
2080 IF (QNI3D(K).GE.1.E-8) THEN
2081 NSAGG(K) = CONS15*ASN(K)*RHO(K)** &
2082 ((2.+BS)/3.)*QNI3D(K)**((2.+BS)/3.)* &
2083 (NS3D(K)*RHO(K))**((4.-BS)/3.)/ &
2084 (RHO(K))
2085 END IF
2086
2087 !.......................................................................
2088 ! ACCRETION OF CLOUD DROPLETS ONTO SNOW/GRAUPEL
2089 ! HERE USE CONTINUOUS COLLECTION EQUATION WITH
2090 ! SIMPLE GRAVITATIONAL COLLECTION KERNEL IGNORING
2091
2092 ! SNOW
2093
2094 IF (QNI3D(K).GE.1.E-8 .AND. QC3D(K).GE.QSMALL) THEN
2095
2096 PSACWS(K) = CONS13*ASN(K)*QC3D(K)*RHO(K)* &
2097 N0S(K)/ &
2098 LAMS(K)**(BS+3.)
2099 NPSACWS(K) = CONS13*ASN(K)*NC3D(K)*RHO(K)* &
2100 N0S(K)/ &
2101 LAMS(K)**(BS+3.)
2102
2103 END IF
2104
2105 !............................................................................
2106 ! COLLECTION OF CLOUD WATER BY GRAUPEL
2107
2108 IF (QG3D(K).GE.1.E-8 .AND. QC3D(K).GE.QSMALL) THEN
2109
2110 PSACWG(K) = CONS14*AGN(K)*QC3D(K)*RHO(K)* &
2111 N0G(K)/ &
2112 LAMG(K)**(BG+3.)
2113 NPSACWG(K) = CONS14*AGN(K)*NC3D(K)*RHO(K)* &
2114 N0G(K)/ &
2115 LAMG(K)**(BG+3.)
2116 END IF
2117
2118 !.......................................................................
2119 ! HM, ADD 12/13/06
2120 ! CLOUD ICE COLLECTING DROPLETS, ASSUME THAT CLOUD ICE MEAN DIAM > 100 MICRON
2121 ! BEFORE RIMING CAN OCCUR
2122 ! ASSUME THAT RIME COLLECTED ON CLOUD ICE DOES NOT LEAD
2123 ! TO HALLET-MOSSOP SPLINTERING
2124
2125 IF (QI3D(K).GE.1.E-8 .AND. QC3D(K).GE.QSMALL) THEN
2126
2127 ! PUT IN SIZE DEPENDENT COLLECTION EFFICIENCY BASED ON STOKES LAW
2128 ! FROM THOMPSON ET AL. 2004, MWR
2129
2130 IF (1./LAMI(K).GE.100.E-6) THEN
2131
2132 PSACWI(K) = CONS16*AIN(K)*QC3D(K)*RHO(K)* &
2133 N0I(K)/ &
2134 LAMI(K)**(BI+3.)
2135 NPSACWI(K) = CONS16*AIN(K)*NC3D(K)*RHO(K)* &
2136 N0I(K)/ &
2137 LAMI(K)**(BI+3.)
2138 END IF
2139 END IF
2140
2141 !.......................................................................
2142 ! ACCRETION OF RAIN WATER BY SNOW
2143 ! FORMULA FROM IKAWA AND SAITO, 1991, USED BY REISNER ET AL, 1998
2144
2145 IF (QR3D(K).GE.1.E-8.AND.QNI3D(K).GE.1.E-8) THEN
2146
2147 UMS = ASN(K)*CONS3/(LAMS(K)**BS)
2148 UMR = ARN(K)*CONS4/(LAMR(K)**BR)
2149 UNS = ASN(K)*CONS5/LAMS(K)**BS
2150 UNR = ARN(K)*CONS6/LAMR(K)**BR
2151
2152 ! SET REASLISTIC LIMITS ON FALLSPEEDS
2153 UMS=MIN(UMS,1.2)
2154 UNS=MIN(UNS,1.2)
2155 UMR=MIN(UMR,9.1)
2156 UNR=MIN(UNR,9.1)
2157
2158 PRACS(K) = CONS41*(((1.2*UMR-0.95*UMS)**2+ &
2159 0.08*UMS*UMR)**0.5*RHO(K)* &
2160 N0RR(K)*N0S(K)/LAMR(K)**3* &
2161 (5./(LAMR(K)**3*LAMS(K))+ &
2162 2./(LAMR(K)**2*LAMS(K)**2)+ &
2163 0.5/(LAMR(k)*LAMS(k)**3)))
2164
2165 NPRACS(K) = CONS32*RHO(K)*(1.7*(UNR-UNS)**2+ &
2166 0.3*UNR*UNS)**0.5*N0RR(K)*N0S(K)* &
2167 (1./(LAMR(K)**3*LAMS(K))+ &
2168 1./(LAMR(K)**2*LAMS(K)**2)+ &
2169 1./(LAMR(K)*LAMS(K)**3))
2170
2171 ! MAKE SURE PRACS DOESN'T EXCEED TOTAL RAIN MIXING RATIO
2172 ! AS THIS MAY OTHERWISE RESULT IN TOO MUCH TRANSFER OF WATER DURING
2173 ! RIME-SPLINTERING
2174
2175 PRACS(K) = MIN(PRACS(K),QR3D(K)/DT)
2176
2177 ! COLLECTION OF SNOW BY RAIN - NEEDED FOR GRAUPEL CONVERSION CALCULATIONS
2178 ! ONLY CALCULATE IF SNOW AND RAIN MIXING RATIOS EXCEED 0.1 G/KG
2179
2180 ! ASSUME COLLECTION OF SNOW BY RAIN PRODUCES GRAUPEL NOT HAIL
2181
2182 IF (IHAIL.EQ.0) THEN
2183 IF (QNI3D(K).GE.0.1E-3.AND.QR3D(K).GE.0.1E-3) THEN
2184 PSACR(K) = CONS31*(((1.2*UMR-0.95*UMS)**2+ &
2185 0.08*UMS*UMR)**0.5*RHO(K)* &
2186 N0RR(K)*N0S(K)/LAMS(K)**3* &
2187 (5./(LAMS(K)**3*LAMR(K))+ &
2188 2./(LAMS(K)**2*LAMR(K)**2)+ &
2189 0.5/(LAMS(K)*LAMR(K)**3)))
2190 END IF
2191 END IF
2192
2193 END IF
2194
2195 !.......................................................................
2196
2197 ! COLLECTION OF RAINWATER BY GRAUPEL, FROM IKAWA AND SAITO 1990,
2198 ! USED BY REISNER ET AL 1998
2199 IF (QR3D(K).GE.1.E-8.AND.QG3D(K).GE.1.E-8) THEN
2200
2201 UMG = AGN(K)*CONS7/(LAMG(K)**BG)
2202 UMR = ARN(K)*CONS4/(LAMR(K)**BR)
2203 UNG = AGN(K)*CONS8/LAMG(K)**BG
2204 UNR = ARN(K)*CONS6/LAMR(K)**BR
2205
2206 ! SET REASLISTIC LIMITS ON FALLSPEEDS
2207 UMG=MIN(UMG,20.)
2208 UNG=MIN(UNG,20.)
2209 UMR=MIN(UMR,9.1)
2210 UNR=MIN(UNR,9.1)
2211
2212 PRACG(K) = CONS41*(((1.2*UMR-0.95*UMG)**2+ &
2213 0.08*UMG*UMR)**0.5*RHO(K)* &
2214 N0RR(K)*N0G(K)/LAMR(K)**3* &
2215 (5./(LAMR(K)**3*LAMG(K))+ &
2216 2./(LAMR(K)**2*LAMG(K)**2)+ &
2217 0.5/(LAMR(k)*LAMG(k)**3)))
2218
2219 NPRACG(K) = CONS32*RHO(K)*(1.7*(UNR-UNG)**2+ &
2220 0.3*UNR*UNG)**0.5*N0RR(K)*N0G(K)* &
2221 (1./(LAMR(K)**3*LAMG(K))+ &
2222 1./(LAMR(K)**2*LAMG(K)**2)+ &
2223 1./(LAMR(K)*LAMG(K)**3))
2224
2225 ! MAKE SURE PRACG DOESN'T EXCEED TOTAL RAIN MIXING RATIO
2226 ! AS THIS MAY OTHERWISE RESULT IN TOO MUCH TRANSFER OF WATER DURING
2227 ! RIME-SPLINTERING
2228
2229 PRACG(K) = MIN(PRACG(K),QR3D(K)/DT)
2230
2231 END IF
2232
2233 !.......................................................................
2234 ! RIME-SPLINTERING - SNOW
2235 ! HALLET-MOSSOP (1974)
2236 ! NUMBER OF SPLINTERS FORMED IS BASED ON MASS OF RIMED WATER
2237
2238 ! DUM1 = MASS OF INDIVIDUAL SPLINTERS
2239
2240 ! HM ADD THRESHOLD SNOW AND DROPLET MIXING RATIO FOR RIME-SPLINTERING
2241 ! TO LIMIT RIME-SPLINTERING IN STRATIFORM CLOUDS
2242 ! THESE THRESHOLDS CORRESPOND WITH GRAUPEL THRESHOLDS IN RH 1984
2243
2244 IF (QNI3D(K).GE.0.1E-3.AND.QC3D(K).GE.0.5E-3) THEN
2245 IF (PSACWS(K).GT.0..OR.PRACS(K).GT.0.) THEN
2246 IF (T3D(K).LT.270.16 .AND. T3D(K).GT.265.16) THEN
2247
2248 IF (T3D(K).GT.270.16) THEN
2249 FMULT = 0.
2250 ELSE IF (T3D(K).LE.270.16.AND.T3D(K).GT.268.16) THEN
2251 FMULT = (270.16-T3D(K))/2.
2252 ELSE IF (T3D(K).GE.265.16.AND.T3D(K).LE.268.16) THEN
2253 FMULT = (T3D(K)-265.16)/3.
2254 ELSE IF (T3D(K).LT.265.16) THEN
2255 FMULT = 0.
2256 END IF
2257
2258 ! 1000 IS TO CONVERT FROM KG TO G
2259
2260 ! SPLINTERING FROM DROPLETS ACCRETED ONTO SNOW
2261
2262 IF (PSACWS(K).GT.0.) THEN
2263 NMULTS(K) = 35.E4*PSACWS(K)*FMULT*1000.
2264 QMULTS(K) = NMULTS(K)*MMULT
2265
2266 ! CONSTRAIN SO THAT TRANSFER OF MASS FROM SNOW TO ICE CANNOT BE MORE MASS
2267 ! THAN WAS RIMED ONTO SNOW
2268
2269 QMULTS(K) = MIN(QMULTS(K),PSACWS(K))
2270 PSACWS(K) = PSACWS(K)-QMULTS(K)
2271
2272 END IF
2273
2274 ! RIMING AND SPLINTERING FROM ACCRETED RAINDROPS
2275
2276 IF (PRACS(K).GT.0.) THEN
2277 NMULTR(K) = 35.E4*PRACS(K)*FMULT*1000.
2278 QMULTR(K) = NMULTR(K)*MMULT
2279
2280 ! CONSTRAIN SO THAT TRANSFER OF MASS FROM SNOW TO ICE CANNOT BE MORE MASS
2281 ! THAN WAS RIMED ONTO SNOW
2282
2283 QMULTR(K) = MIN(QMULTR(K),PRACS(K))
2284
2285 PRACS(K) = PRACS(K)-QMULTR(K)
2286
2287 END IF
2288
2289 END IF
2290 END IF
2291 END IF
2292
2293 !.......................................................................
2294 ! RIME-SPLINTERING - GRAUPEL
2295 ! HALLET-MOSSOP (1974)
2296 ! NUMBER OF SPLINTERS FORMED IS BASED ON MASS OF RIMED WATER
2297
2298 ! DUM1 = MASS OF INDIVIDUAL SPLINTERS
2299
2300 ! HM ADD THRESHOLD SNOW MIXING RATIO FOR RIME-SPLINTERING
2301 ! TO LIMIT RIME-SPLINTERING IN STRATIFORM CLOUDS
2302
2303 ! ONLY CALCULATE FOR GRAUPEL NOT HAIL
2304 IF (IHAIL.EQ.0) THEN
2305 IF (PSACWG(K).GT.0..OR.PRACG(K).GT.0.) THEN
2306 IF (T3D(K).LT.270.16 .AND. T3D(K).GT.265.16) THEN
2307
2308 IF (T3D(K).GT.270.16) THEN
2309 FMULT = 0.
2310 ELSE IF (T3D(K).LE.270.16.AND.T3D(K).GT.268.16) THEN
2311 FMULT = (270.16-T3D(K))/2.
2312 ELSE IF (T3D(K).GE.265.16.AND.T3D(K).LE.268.16) THEN
2313 FMULT = (T3D(K)-265.16)/3.
2314 ELSE IF (T3D(K).LT.265.16) THEN
2315 FMULT = 0.
2316 END IF
2317
2318 ! 1000 IS TO CONVERT FROM KG TO G
2319
2320 ! SPLINTERING FROM DROPLETS ACCRETED ONTO GRAUPEL
2321
2322 IF (PSACWG(K).GT.0.) THEN
2323 NMULTG(K) = 35.E4*PSACWG(K)*FMULT*1000.
2324 QMULTG(K) = NMULTG(K)*MMULT
2325
2326 ! CONSTRAIN SO THAT TRANSFER OF MASS FROM GRAUPEL TO ICE CANNOT BE MORE MASS
2327 ! THAN WAS RIMED ONTO GRAUPEL
2328
2329 QMULTG(K) = MIN(QMULTG(K),PSACWG(K))
2330 PSACWG(K) = PSACWG(K)-QMULTG(K)
2331
2332 END IF
2333
2334 ! RIMING AND SPLINTERING FROM ACCRETED RAINDROPS
2335
2336 IF (PRACG(K).GT.0.) THEN
2337 NMULTRG(K) = 35.E4*PRACG(K)*FMULT*1000.
2338 QMULTRG(K) = NMULTRG(K)*MMULT
2339
2340 ! CONSTRAIN SO THAT TRANSFER OF MASS FROM GRAUPEL TO ICE CANNOT BE MORE MASS
2341 ! THAN WAS RIMED ONTO GRAUPEL
2342
2343 QMULTRG(K) = MIN(QMULTRG(K),PRACG(K))
2344 PRACG(K) = PRACG(K)-QMULTRG(K)
2345
2346 END IF
2347
2348 END IF
2349 END IF
2350 END IF
2351
2352 !........................................................................
2353 ! CONVERSION OF RIMED CLOUD WATER ONTO SNOW TO GRAUPEL
2354 ! ASSUME CONVERTED SNOW FORMS GRAUPEL NOT HAIL
2355 ! HAIL ASSUMED TO ONLY FORM BY FREEZING OF RAIN
2356 ! OR COLLISIONS OF RAIN WITH CLOUD ICE
2357
2358 IF (IHAIL.EQ.0) THEN
2359 IF (PSACWS(K).GT.0.) THEN
2360 ! ONLY ALLOW CONVERSION IF QNI > 0.1 AND QC > 0.5 G/KG FOLLOWING RUTLEDGE AND HOBBS (1984)
2361 IF (QNI3D(K).GE.0.1E-3.AND.QC3D(K).GE.0.5E-3) THEN
2362
2363 ! PORTION OF RIMING CONVERTED TO GRAUPEL (REISNER ET AL. 1998, ORIGINALLY IS1991)
2364 PGSACW(K) = MIN(PSACWS(K),CONS17*DT*N0S(K)*QC3D(K)*QC3D(K)* &
2365 ASN(K)*ASN(K)/ &
2366 (RHO(K)*LAMS(K)**(2.*BS+2.)))
2367
2368 ! MIX RAT CONVERTED INTO GRAUPEL AS EMBRYO (REISNER ET AL. 1998, ORIG M1990)
2369 DUM = MAX(RHOSN/(RHOG-RHOSN)*PGSACW(K),0.)
2370
2371 ! NUMBER CONCENTRAITON OF EMBRYO GRAUPEL FROM RIMING OF SNOW
2372 NSCNG(K) = DUM/MG0*RHO(K)
2373 ! LIMIT MAX NUMBER CONVERTED TO SNOW NUMBER
2374 NSCNG(K) = MIN(NSCNG(K),NS3D(K)/DT)
2375
2376 ! PORTION OF RIMING LEFT FOR SNOW
2377 PSACWS(K) = PSACWS(K) - PGSACW(K)
2378 END IF
2379 END IF
2380
2381 ! CONVERSION OF RIMED RAINWATER ONTO SNOW CONVERTED TO GRAUPEL
2382
2383 IF (PRACS(K).GT.0.) THEN
2384 ! ONLY ALLOW CONVERSION IF QNI > 0.1 AND QR > 0.1 G/KG FOLLOWING RUTLEDGE AND HOBBS (1984)
2385 IF (QNI3D(K).GE.0.1E-3.AND.QR3D(K).GE.0.1E-3) THEN
2386 ! PORTION OF COLLECTED RAINWATER CONVERTED TO GRAUPEL (REISNER ET AL. 1998)
2387 DUM = CONS18*(4./LAMS(K))**3*(4./LAMS(K))**3 &
2388 /(CONS18*(4./LAMS(K))**3*(4./LAMS(K))**3+ &
2389 CONS19*(4./LAMR(K))**3*(4./LAMR(K))**3)
2390 DUM=MIN(DUM,1.)
2391 DUM=MAX(DUM,0.)
2392 PGRACS(K) = (1.-DUM)*PRACS(K)
2393 NGRACS(K) = (1.-DUM)*NPRACS(K)
2394 ! LIMIT MAX NUMBER CONVERTED TO MIN OF EITHER RAIN OR SNOW NUMBER CONCENTRATION
2395 NGRACS(K) = MIN(NGRACS(K),NR3D(K)/DT)
2396 NGRACS(K) = MIN(NGRACS(K),NS3D(K)/DT)
2397
2398 ! AMOUNT LEFT FOR SNOW PRODUCTION
2399 PRACS(K) = PRACS(K) - PGRACS(K)
2400 NPRACS(K) = NPRACS(K) - NGRACS(K)
2401 ! CONVERSION TO GRAUPEL DUE TO COLLECTION OF SNOW BY RAIN
2402 PSACR(K)=PSACR(K)*(1.-DUM)
2403 END IF
2404 END IF
2405 END IF
2406
2407 !.......................................................................
2408 ! FREEZING OF RAIN DROPS
2409 ! FREEZING ALLOWED BELOW -4 C
2410
2411 IF (T3D(K).LT.269.15.AND.QR3D(K).GE.QSMALL) THEN
2412
2413 ! IMMERSION FREEZING (BIGG 1953)
2414 MNUCCR(K) = CONS20*NR3D(K)*EXP(AIMM*(273.15-T3D(K)))/LAMR(K)**3 &
2415 /LAMR(K)**3
2416
2417 NNUCCR(K) = PI*NR3D(K)*BIMM*EXP(AIMM*(273.15-T3D(K)))/LAMR(K)**3
2418
2419 ! PREVENT DIVERGENCE BETWEEN MIXING RATIO AND NUMBER CONC
2420 NNUCCR(K) = MIN(NNUCCR(K),NR3D(K)/DT)
2421
2422 END IF
2423
2424 !.......................................................................
2425 ! ACCRETION OF CLOUD LIQUID WATER BY RAIN
2426 ! CONTINUOUS COLLECTION EQUATION WITH
2427 ! GRAVITATIONAL COLLECTION KERNEL, DROPLET FALL SPEED NEGLECTED
2428
2429 IF (QR3D(K).GE.1.E-8 .AND. QC3D(K).GE.1.E-8) THEN
2430
2431 ! 12/13/06 HM ADD, REPLACE WITH NEWER FORMULA FROM
2432 ! KHAIROUTDINOV AND KOGAN 2000, MWR
2433
2434 DUM=(QC3D(K)*QR3D(K))
2435 PRA(K) = 67.*(DUM)**1.15
2436 NPRA(K) = PRA(K)/(QC3D(K)/NC3D(K))
2437
2438 END IF
2439 !.......................................................................
2440 ! SELF-COLLECTION OF RAIN DROPS
2441 ! FROM BEHENG(1994)
2442 ! FROM NUMERICAL SIMULATION OF THE STOCHASTIC COLLECTION EQUATION
2443 ! AS DESCRINED ABOVE FOR AUTOCONVERSION
2444
2445 IF (QR3D(K).GE.1.E-8) THEN
2446 NRAGG(K) = -8.*NR3D(K)*QR3D(K)*RHO(K)
2447 END IF
2448
2449 !.......................................................................
2450 ! AUTOCONVERSION OF CLOUD ICE TO SNOW
2451 ! FOLLOWING HARRINGTON ET AL. (1995) WITH MODIFICATION
2452 ! HERE IT IS ASSUMED THAT AUTOCONVERSION CAN ONLY OCCUR WHEN THE
2453 ! ICE IS GROWING, I.E. IN CONDITIONS OF ICE SUPERSATURATION
2454
2455 IF (QI3D(K).GE.1.E-8 .AND.QVQVSI(K).GE.1.) THEN
2456
2457 ! COFFI = 2./LAMI(K)
2458 ! IF (COFFI.GE.DCS) THEN
2459 NPRCI(K) = CONS21*(QV3D(K)-QVI(K))*RHO(K) &
2460 *N0I(K)*EXP(-LAMI(K)*DCS)*DV(K)/ABI(K)
2461 PRCI(K) = CONS22*NPRCI(K)
2462 NPRCI(K) = MIN(NPRCI(K),NI3D(K)/DT)
2463
2464 ! END IF
2465 END IF
2466
2467 !.......................................................................
2468 ! ACCRETION OF CLOUD ICE BY SNOW
2469 ! FOR THIS CALCULATION, IT IS ASSUMED THAT THE VS >> VI
2470 ! AND DS >> DI FOR CONTINUOUS COLLECTION
2471
2472 IF (QNI3D(K).GE.1.E-8 .AND. QI3D(K).GE.QSMALL) THEN
2473 PRAI(K) = CONS23*ASN(K)*QI3D(K)*RHO(K)*N0S(K)/ &
2474 LAMS(K)**(BS+3.)
2475 NPRAI(K) = CONS23*ASN(K)*NI3D(K)* &
2476 RHO(K)*N0S(K)/ &
2477 LAMS(K)**(BS+3.)
2478 NPRAI(K)=MIN(NPRAI(K),NI3D(K)/DT)
2479 END IF
2480
2481 !.......................................................................
2482 ! HM, ADD 12/13/06, COLLISION OF RAIN AND ICE TO PRODUCE SNOW OR GRAUPEL
2483 ! FOLLOWS REISNER ET AL. 1998
2484 ! ASSUMED FALLSPEED AND SIZE OF ICE CRYSTAL << THAN FOR RAIN
2485
2486 IF (QR3D(K).GE.1.E-8.AND.QI3D(K).GE.1.E-8.AND.T3D(K).LE.273.15) THEN
2487
2488 ! ALLOW GRAUPEL FORMATION FROM RAIN-ICE COLLISIONS ONLY IF RAIN MIXING RATIO > 0.1 G/KG,
2489 ! OTHERWISE ADD TO SNOW
2490
2491 IF (QR3D(K).GE.0.1E-3) THEN
2492 NIACR(K)=CONS24*NI3D(K)*N0RR(K)*ARN(K) &
2493 /LAMR(K)**(BR+3.)*RHO(K)
2494 PIACR(K)=CONS25*NI3D(K)*N0RR(K)*ARN(K) &
2495 /LAMR(K)**(BR+3.)/LAMR(K)**3*RHO(K)
2496 PRACI(K)=CONS24*QI3D(K)*N0RR(K)*ARN(K)/ &
2497 LAMR(K)**(BR+3.)*RHO(K)
2498 NIACR(K)=MIN(NIACR(K),NR3D(K)/DT)
2499 NIACR(K)=MIN(NIACR(K),NI3D(K)/DT)
2500 ELSE
2501 NIACRS(K)=CONS24*NI3D(K)*N0RR(K)*ARN(K) &
2502 /LAMR(K)**(BR+3.)*RHO(K)
2503 PIACRS(K)=CONS25*NI3D(K)*N0RR(K)*ARN(K) &
2504 /LAMR(K)**(BR+3.)/LAMR(K)**3*RHO(K)
2505 PRACIS(K)=CONS24*QI3D(K)*N0RR(K)*ARN(K)/ &
2506 LAMR(K)**(BR+3.)*RHO(K)
2507 NIACRS(K)=MIN(NIACRS(K),NR3D(K)/DT)
2508 NIACRS(K)=MIN(NIACRS(K),NI3D(K)/DT)
2509 END IF
2510 END IF
2511
2512 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
2513 ! NUCLEATION OF CLOUD ICE FROM HOMOGENEOUS AND HETEROGENEOUS FREEZING ON AEROSOL
2514
2515 IF (INUC.EQ.0) THEN
2516
2517 ! FREEZING OF AEROSOL ONLY ALLOWED BELOW -5 C
2518 ! AND ABOVE DELIQUESCENCE THRESHOLD OF 80%
2519 ! AND ABOVE ICE SATURATION
2520
2521 ! add threshold according to Greg Thomspon
2522
2523 if ((QVQVS(K).GE.0.999.and.T3D(K).le.265.15).or. &
2524 QVQVSI(K).ge.1.08) then
2525
2526 ! hm, modify dec. 5, 2006, replace with cooper curve
2527 kc2 = 0.005*exp(0.304*(273.15-T3D(K)))*1000. ! convert from L-1 to m-3
2528 ! limit to 500 L-1
2529 kc2 = min(kc2,500.e3)
2530 kc2=MAX(kc2/rho(k),0.) ! convert to kg-1
2531
2532 IF (KC2.GT.NI3D(K)+NS3D(K)+NG3D(K)) THEN
2533 NNUCCD(K) = (KC2-NI3D(K)-NS3D(K)-NG3D(K))/DT
2534 MNUCCD(K) = NNUCCD(K)*MI0
2535 END IF
2536
2537 END IF
2538
2539 ELSE IF (INUC.EQ.1) THEN
2540
2541 IF (T3D(K).LT.273.15.AND.QVQVSI(K).GT.1.) THEN
2542
2543 KC2 = 0.16*1000./RHO(K) ! CONVERT FROM L-1 TO KG-1
2544 IF (KC2.GT.NI3D(K)+NS3D(K)+NG3D(K)) THEN
2545 NNUCCD(K) = (KC2-NI3D(K)-NS3D(K)-NG3D(K))/DT
2546 MNUCCD(K) = NNUCCD(K)*MI0
2547 END IF
2548 END IF
2549
2550 END IF
2551
2552 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
2553
2554 101 CONTINUE
2555
2556 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
2557 ! CALCULATE EVAP/SUB/DEP TERMS FOR QI,QNI,QR
2558
2559 ! NO VENTILATION FOR CLOUD ICE
2560
2561 IF (QI3D(K).GE.QSMALL) THEN
2562
2563 EPSI = 2.*PI*N0I(K)*RHO(K)*DV(K)/(LAMI(K)*LAMI(K))
2564
2565 ELSE
2566 EPSI = 0.
2567 END IF
2568
2569 IF (QNI3D(K).GE.QSMALL) THEN
2570 EPSS = 2.*PI*N0S(K)*RHO(K)*DV(K)* &
2571 (F1S/(LAMS(K)*LAMS(K))+ &
2572 F2S*(ASN(K)*RHO(K)/MU(K))**0.5* &
2573 SC(K)**(1./3.)*CONS10/ &
2574 (LAMS(K)**CONS35))
2575 ELSE
2576 EPSS = 0.
2577 END IF
2578
2579 IF (QG3D(K).GE.QSMALL) THEN
2580 EPSG = 2.*PI*N0G(K)*RHO(K)*DV(K)* &
2581 (F1S/(LAMG(K)*LAMG(K))+ &
2582 F2S*(AGN(K)*RHO(K)/MU(K))**0.5* &
2583 SC(K)**(1./3.)*CONS11/ &
2584 (LAMG(K)**CONS36))
2585
2586
2587 ELSE
2588 EPSG = 0.
2589 END IF
2590
2591 IF (QR3D(K).GE.QSMALL) THEN
2592 EPSR = 2.*PI*N0RR(K)*RHO(K)*DV(K)* &
2593 (F1R/(LAMR(K)*LAMR(K))+ &
2594 F2R*(ARN(K)*RHO(K)/MU(K))**0.5* &
2595 SC(K)**(1./3.)*CONS9/ &
2596 (LAMR(K)**CONS34))
2597 ELSE
2598 EPSR = 0.
2599 END IF
2600
2601 ! ONLY INCLUDE REGION OF ICE SIZE DIST < DCS
2602 ! DUM IS FRACTION OF D*N(D) < DCS
2603
2604 ! LOGIC BELOW FOLLOWS THAT OF HARRINGTON ET AL. 1995 (JAS)
2605 IF (QI3D(K).GE.QSMALL) THEN
2606 DUM=(1.-EXP(-LAMI(K)*DCS)*(1.+LAMI(K)*DCS))
2607 PRD(K) = EPSI*(QV3D(K)-QVI(K))/ABI(K)*DUM
2608 ELSE
2609 DUM=0.
2610 END IF
2611 ! ADD DEPOSITION IN TAIL OF ICE SIZE DIST TO SNOW IF SNOW IS PRESENT
2612 IF (QNI3D(K).GE.QSMALL) THEN
2613 PRDS(K) = EPSS*(QV3D(K)-QVI(K))/ABI(K)+ &
2614 EPSI*(QV3D(K)-QVI(K))/ABI(K)*(1.-DUM)
2615 ! OTHERWISE ADD TO CLOUD ICE
2616 ELSE
2617 PRD(K) = PRD(K)+EPSI*(QV3D(K)-QVI(K))/ABI(K)*(1.-DUM)
2618 END IF
2619 ! VAPOR DPEOSITION ON GRAUPEL
2620 PRDG(K) = EPSG*(QV3D(K)-QVI(K))/ABI(K)
2621
2622 ! NO CONDENSATION ONTO RAIN, ONLY EVAP
2623
2624 IF (QV3D(K).LT.QVS(K)) THEN
2625 PRE(K) = EPSR*(QV3D(K)-QVS(K))/AB(K)
2626 PRE(K) = MIN(PRE(K),0.)
2627 ELSE
2628 PRE(K) = 0.
2629 END IF
2630
2631 ! MAKE SURE NOT PUSHED INTO ICE SUPERSAT/SUBSAT
2632 ! FORMULA FROM REISNER 2 SCHEME
2633
2634 DUM = (QV3D(K)-QVI(K))/DT
2635
2636 FUDGEF = 0.9999
2637 SUM_DEP = PRD(K)+PRDS(K)+MNUCCD(K)+PRDG(K)
2638
2639 IF( (DUM.GT.0. .AND. SUM_DEP.GT.DUM*FUDGEF) .OR. &
2640 (DUM.LT.0. .AND. SUM_DEP.LT.DUM*FUDGEF) ) THEN
2641 MNUCCD(K) = FUDGEF*MNUCCD(K)*DUM/SUM_DEP
2642 PRD(K) = FUDGEF*PRD(K)*DUM/SUM_DEP
2643 PRDS(K) = FUDGEF*PRDS(K)*DUM/SUM_DEP
2644 PRDG(K) = FUDGEF*PRDG(K)*DUM/SUM_DEP
2645 ENDIF
2646
2647 ! IF CLOUD ICE/SNOW/GRAUPEL VAP DEPOSITION IS NEG, THEN ASSIGN TO SUBLIMATION PROCESSES
2648
2649 IF (PRD(K).LT.0.) THEN
2650 EPRD(K)=PRD(K)
2651 PRD(K)=0.
2652 END IF
2653 IF (PRDS(K).LT.0.) THEN
2654 EPRDS(K)=PRDS(K)
2655 PRDS(K)=0.
2656 END IF
2657 IF (PRDG(K).LT.0.) THEN
2658 EPRDG(K)=PRDG(K)
2659 PRDG(K)=0.
2660 END IF
2661
2662 !.......................................................................
2663 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
2664
2665 ! CONSERVATION OF WATER
2666 ! THIS IS ADOPTED LOOSELY FROM MM5 RESINER CODE. HOWEVER, HERE WE
2667 ! ONLY ADJUST PROCESSES THAT ARE NEGATIVE, RATHER THAN ALL PROCESSES.
2668
2669 ! IF MIXING RATIOS LESS THAN QSMALL, THEN NO DEPLETION OF WATER
2670 ! THROUGH MICROPHYSICAL PROCESSES, SKIP CONSERVATION
2671
2672 ! NOTE: CONSERVATION CHECK NOT APPLIED TO NUMBER CONCENTRATION SPECIES. ADDITIONAL CATCH
2673 ! BELOW WILL PREVENT NEGATIVE NUMBER CONCENTRATION
2674 ! FOR EACH MICROPHYSICAL PROCESS WHICH PROVIDES A SOURCE FOR NUMBER, THERE IS A CHECK
2675 ! TO MAKE SURE THAT CAN'T EXCEED TOTAL NUMBER OF DEPLETED SPECIES WITH THE TIME
2676 ! STEP
2677
2678 !****SENSITIVITY - NO ICE
2679
2680 IF (ILIQ.EQ.1) THEN
2681 MNUCCC(K)=0.
2682 NNUCCC(K)=0.
2683 MNUCCR(K)=0.
2684 NNUCCR(K)=0.
2685 MNUCCD(K)=0.
2686 NNUCCD(K)=0.
2687 END IF
2688
2689 ! ****SENSITIVITY - NO GRAUPEL
2690 IF (IGRAUP.EQ.1) THEN
2691 PRACG(K) = 0.
2692 PSACR(K) = 0.
2693 PSACWG(K) = 0.
2694 PGSACW(K) = 0.
2695 PGRACS(K) = 0.
2696 PRDG(K) = 0.
2697 EPRDG(K) = 0.
2698 EVPMG(K) = 0.
2699 PGMLT(K) = 0.
2700 NPRACG(K) = 0.
2701 NPSACWG(K) = 0.
2702 NSCNG(K) = 0.
2703 NGRACS(K) = 0.
2704 NSUBG(K) = 0.
2705 NGMLTG(K) = 0.
2706 NGMLTR(K) = 0.
2707 END IF
2708
2709 ! CONSERVATION OF QC
2710
2711 DUM = (PRC(K)+PRA(K)+MNUCCC(K)+PSACWS(K)+PSACWI(K)+QMULTS(K)+PSACWG(K)+PGSACW(K)+QMULTG(K))*DT
2712
2713 IF (DUM.GT.QC3D(K).AND.QC3D(K).GE.QSMALL) THEN
2714 RATIO = QC3D(K)/DUM
2715
2716 PRC(K) = PRC(K)*RATIO
2717 PRA(K) = PRA(K)*RATIO
2718 MNUCCC(K) = MNUCCC(K)*RATIO
2719 PSACWS(K) = PSACWS(K)*RATIO
2720 PSACWI(K) = PSACWI(K)*RATIO
2721 QMULTS(K) = QMULTS(K)*RATIO
2722 QMULTG(K) = QMULTG(K)*RATIO
2723 PSACWG(K) = PSACWG(K)*RATIO
2724 PGSACW(K) = PGSACW(K)*RATIO
2725 END IF
2726
2727 ! CONSERVATION OF QI
2728
2729 DUM = (-PRD(K)-MNUCCC(K)+PRCI(K)+PRAI(K)-QMULTS(K)-QMULTG(K)-QMULTR(K)-QMULTRG(K) &
2730 -MNUCCD(K)+PRACI(K)+PRACIS(K)-EPRD(K)-PSACWI(K))*DT
2731
2732 IF (DUM.GT.QI3D(K).AND.QI3D(K).GE.QSMALL) THEN
2733
2734 RATIO = (QI3D(K)/DT+PRD(K)+MNUCCC(K)+QMULTS(K)+QMULTG(K)+QMULTR(K)+QMULTRG(K)+ &
2735 MNUCCD(K)+PSACWI(K))/ &
2736 (PRCI(K)+PRAI(K)+PRACI(K)+PRACIS(K)-EPRD(K))
2737
2738 PRCI(K) = PRCI(K)*RATIO
2739 PRAI(K) = PRAI(K)*RATIO
2740 PRACI(K) = PRACI(K)*RATIO
2741 PRACIS(K) = PRACIS(K)*RATIO
2742 EPRD(K) = EPRD(K)*RATIO
2743
2744 END IF
2745
2746 ! CONSERVATION OF QR
2747
2748 DUM=((PRACS(K)-PRE(K))+(QMULTR(K)+QMULTRG(K)-PRC(K))+(MNUCCR(K)-PRA(K))+ &
2749 PIACR(K)+PIACRS(K)+PGRACS(K)+PRACG(K))*DT
2750
2751 IF (DUM.GT.QR3D(K).AND.QR3D(K).GE.QSMALL) THEN
2752
2753 RATIO = (QR3D(K)/DT+PRC(K)+PRA(K))/ &
2754 (-PRE(K)+QMULTR(K)+QMULTRG(K)+PRACS(K)+MNUCCR(K)+PIACR(K)+PIACRS(K)+PGRACS(K)+PRACG(K))
2755
2756 PRE(K) = PRE(K)*RATIO
2757 PRACS(K) = PRACS(K)*RATIO
2758 QMULTR(K) = QMULTR(K)*RATIO
2759 QMULTRG(K) = QMULTRG(K)*RATIO
2760 MNUCCR(K) = MNUCCR(K)*RATIO
2761 PIACR(K) = PIACR(K)*RATIO
2762 PIACRS(K) = PIACRS(K)*RATIO
2763 PGRACS(K) = PGRACS(K)*RATIO
2764 PRACG(K) = PRACG(K)*RATIO
2765
2766 END IF
2767
2768 ! CONSERVATION OF QNI
2769 ! CONSERVATION FOR GRAUPEL SCHEME
2770
2771 IF (IGRAUP.EQ.0) THEN
2772
2773 DUM = (-PRDS(K)-PSACWS(K)-PRAI(K)-PRCI(K)-PRACS(K)-EPRDS(K)+PSACR(K)-PIACRS(K)-PRACIS(K))*DT
2774
2775 IF (DUM.GT.QNI3D(K).AND.QNI3D(K).GE.QSMALL) THEN
2776
2777 RATIO = (QNI3D(K)/DT+PRDS(K)+PSACWS(K)+PRAI(K)+PRCI(K)+PRACS(K)+PIACRS(K)+PRACIS(K))/(-EPRDS(K)+PSACR(K))
2778
2779 EPRDS(K) = EPRDS(K)*RATIO
2780 PSACR(K) = PSACR(K)*RATIO
2781
2782 END IF
2783
2784 ! FOR NO GRAUPEL, NEED TO INCLUDE FREEZING OF RAIN FOR SNOW
2785 ELSE IF (IGRAUP.EQ.1) THEN
2786
2787 DUM = (-PRDS(K)-PSACWS(K)-PRAI(K)-PRCI(K)-PRACS(K)-EPRDS(K)+PSACR(K)-PIACRS(K)-PRACIS(K)-MNUCCR(K))*DT
2788
2789 IF (DUM.GT.QNI3D(K).AND.QNI3D(K).GE.QSMALL) THEN
2790
2791 RATIO = (QNI3D(K)/DT+PRDS(K)+PSACWS(K)+PRAI(K)+PRCI(K)+PRACS(K)+PIACRS(K)+PRACIS(K)+MNUCCR(K))/(-EPRDS(K)+PSACR(K))
2792
2793 EPRDS(K) = EPRDS(K)*RATIO
2794 PSACR(K) = PSACR(K)*RATIO
2795
2796 END IF
2797
2798 END IF
2799
2800 ! CONSERVATION OF QG
2801
2802 DUM = (-PSACWG(K)-PRACG(K)-PGSACW(K)-PGRACS(K)-PRDG(K)-MNUCCR(K)-EPRDG(K)-PIACR(K)-PRACI(K)-PSACR(K))*DT
2803
2804 IF (DUM.GT.QG3D(K).AND.QG3D(K).GE.QSMALL) THEN
2805
2806 RATIO = (QG3D(K)/DT+PSACWG(K)+PRACG(K)+PGSACW(K)+PGRACS(K)+PRDG(K)+MNUCCR(K)+PSACR(K)+&
2807 PIACR(K)+PRACI(K))/(-EPRDG(K))
2808
2809 EPRDG(K) = EPRDG(K)*RATIO
2810
2811 END IF
2812
2813 ! TENDENCIES
2814
2815 QV3DTEN(K) = QV3DTEN(K)+(-PRE(K)-PRD(K)-PRDS(K)-MNUCCD(K)-EPRD(K)-EPRDS(K)-PRDG(K)-EPRDG(K))
2816 T3DTEN(K) = T3DTEN(K)+(PRE(K) &
2817 *XXLV(K)+(PRD(K)+PRDS(K)+ &
2818 MNUCCD(K)+EPRD(K)+EPRDS(K)+PRDG(K)+EPRDG(K))*XXLS(K)+ &
2819 (PSACWS(K)+PSACWI(K)+MNUCCC(K)+MNUCCR(K)+ &
2820 QMULTS(K)+QMULTG(K)+QMULTR(K)+QMULTRG(K)+PRACS(K) &
2821 +PSACWG(K)+PRACG(K)+PGSACW(K)+PGRACS(K))*XLF(K))/CPM(K)
2822
2823 QC3DTEN(K) = QC3DTEN(K)+ &
2824 (-PRA(K)-PRC(K)-MNUCCC(K)+PCC(K)- &
2825 PSACWS(K)-PSACWI(K)-QMULTS(K)-QMULTG(K)-PSACWG(K)-PGSACW(K))
2826 QI3DTEN(K) = QI3DTEN(K)+ &
2827 (PRD(K)+EPRD(K)+PSACWI(K)+MNUCCC(K)-PRCI(K)- &
2828 PRAI(K)+QMULTS(K)+QMULTG(K)+QMULTR(K)+QMULTRG(K)+MNUCCD(K)-PRACI(K)-PRACIS(K))
2829 QR3DTEN(K) = QR3DTEN(K)+ &
2830 (PRE(K)+PRA(K)+PRC(K)-PRACS(K)-MNUCCR(K)-QMULTR(K)-QMULTRG(K) &
2831 -PIACR(K)-PIACRS(K)-PRACG(K)-PGRACS(K))
2832
2833 IF (IGRAUP.EQ.0) THEN
2834
2835 QNI3DTEN(K) = QNI3DTEN(K)+ &
2836 (PRAI(K)+PSACWS(K)+PRDS(K)+PRACS(K)+PRCI(K)+EPRDS(K)-PSACR(K)+PIACRS(K)+PRACIS(K))
2837 NS3DTEN(K) = NS3DTEN(K)+(NSAGG(K)+NPRCI(K)-NSCNG(K)-NGRACS(K)+NIACRS(K))
2838 QG3DTEN(K) = QG3DTEN(K)+(PRACG(K)+PSACWG(K)+PGSACW(K)+PGRACS(K)+ &
2839 PRDG(K)+EPRDG(K)+MNUCCR(K)+PIACR(K)+PRACI(K)+PSACR(K))
2840 NG3DTEN(K) = NG3DTEN(K)+(NSCNG(K)+NGRACS(K)+NNUCCR(K)+NIACR(K))
2841
2842 ! FOR NO GRAUPEL, NEED TO INCLUDE FREEZING OF RAIN FOR SNOW
2843 ELSE IF (IGRAUP.EQ.1) THEN
2844
2845 QNI3DTEN(K) = QNI3DTEN(K)+ &
2846 (PRAI(K)+PSACWS(K)+PRDS(K)+PRACS(K)+PRCI(K)+EPRDS(K)-PSACR(K)+PIACRS(K)+PRACIS(K)+MNUCCR(K))
2847 NS3DTEN(K) = NS3DTEN(K)+(NSAGG(K)+NPRCI(K)-NSCNG(K)-NGRACS(K)+NIACRS(K)+NNUCCR(K))
2848
2849 END IF
2850
2851 NC3DTEN(K) = NC3DTEN(K)+(-NNUCCC(K)-NPSACWS(K) &
2852 -NPRA(K)-NPRC(K)-NPSACWI(K)-NPSACWG(K))
2853
2854 NI3DTEN(K) = NI3DTEN(K)+ &
2855 (NNUCCC(K)-NPRCI(K)-NPRAI(K)+NMULTS(K)+NMULTG(K)+NMULTR(K)+NMULTRG(K)+ &
2856 NNUCCD(K)-NIACR(K)-NIACRS(K))
2857
2858 NR3DTEN(K) = NR3DTEN(K)+(NPRC1(K)-NPRACS(K)-NNUCCR(K) &
2859 +NRAGG(K)-NIACR(K)-NIACRS(K)-NPRACG(K)-NGRACS(K))
2860
2861 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
2862 ! NOW CALCULATE SATURATION ADJUSTMENT TO CONDENSE EXTRA VAPOR ABOVE
2863 ! WATER SATURATION
2864
2865 DUMT = T3D(K)+DT*T3DTEN(K)
2866 DUMQV = QV3D(K)+DT*QV3DTEN(K)
2867 DUMQSS = 0.622*POLYSVP(DUMT,0)/ (PRES(K)-POLYSVP(DUMT,0))
2868 DUMQC = QC3D(K)+DT*QC3DTEN(K)
2869 DUMQC = MAX(DUMQC,0.)
2870
2871 ! SATURATION ADJUSTMENT FOR LIQUID
2872
2873 DUMS = DUMQV-DUMQSS
2874 PCC(K) = DUMS/(1.+XXLV(K)**2*DUMQSS/(CPM(K)*RV*DUMT**2))/DT
2875 IF (PCC(K)*DT+DUMQC.LT.0.) THEN
2876 PCC(K) = -DUMQC/DT
2877 END IF
2878
2879 QV3DTEN(K) = QV3DTEN(K)-PCC(K)
2880 T3DTEN(K) = T3DTEN(K)+PCC(K)*XXLV(K)/CPM(K)
2881 QC3DTEN(K) = QC3DTEN(K)+PCC(K)
2882
2883 !.......................................................................
2884 ! ACTIVATION OF CLOUD DROPLETS
2885 ! ACTIVATION OF DROPLET CURRENTLY NOT CALCULATED
2886 ! DROPLET CONCENTRATION IS SPECIFIED !!!!!
2887
2888 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
2889 ! SUBLIMATE, MELT, OR EVAPORATE NUMBER CONCENTRATION
2890 ! THIS FORMULATION ASSUMES 1:1 RATIO BETWEEN MASS LOSS AND
2891 ! LOSS OF NUMBER CONCENTRATION
2892
2893 ! IF (PCC(K).LT.0.) THEN
2894 ! DUM = PCC(K)*DT/QC3D(K)
2895 ! DUM = MAX(-1.,DUM)
2896 ! NSUBC(K) = DUM*NC3D(K)/DT
2897 ! END IF
2898
2899 IF (EPRD(K).LT.0.) THEN
2900 DUM = EPRD(K)*DT/QI3D(K)
2901 DUM = MAX(-1.,DUM)
2902 NSUBI(K) = DUM*NI3D(K)/DT
2903 END IF
2904 IF (EPRDS(K).LT.0.) THEN
2905 DUM = EPRDS(K)*DT/QNI3D(K)
2906 DUM = MAX(-1.,DUM)
2907 NSUBS(K) = DUM*NS3D(K)/DT
2908 END IF
2909 IF (PRE(K).LT.0.) THEN
2910 DUM = PRE(K)*DT/QR3D(K)
2911 DUM = MAX(-1.,DUM)
2912 NSUBR(K) = DUM*NR3D(K)/DT
2913 END IF
2914 IF (EPRDG(K).LT.0.) THEN
2915 DUM = EPRDG(K)*DT/QG3D(K)
2916 DUM = MAX(-1.,DUM)
2917 NSUBG(K) = DUM*NG3D(K)/DT
2918 END IF
2919
2920 ! nsubr(k)=0.
2921 ! nsubs(k)=0.
2922 ! nsubg(k)=0.
2923
2924 ! UPDATE TENDENCIES
2925
2926 ! NC3DTEN(K) = NC3DTEN(K)+NSUBC(K)
2927 NI3DTEN(K) = NI3DTEN(K)+NSUBI(K)
2928 NS3DTEN(K) = NS3DTEN(K)+NSUBS(K)
2929 NG3DTEN(K) = NG3DTEN(K)+NSUBG(K)
2930 NR3DTEN(K) = NR3DTEN(K)+NSUBR(K)
2931
2932 END IF !!!!!! TEMPERATURE
2933
2934 ! SWITCH LTRUE TO 1, SINCE HYDROMETEORS ARE PRESENT
2935 LTRUE = 1
2936
2937 200 CONTINUE
2938
2939 END DO
2940
2941 ! INITIALIZE PRECIP AND SNOW RATES
2942 PRECRT = 0.
2943 SNOWRT = 0.
2944
2945 ! IF THERE ARE NO HYDROMETEORS, THEN SKIP TO END OF SUBROUTINE
2946
2947 IF (LTRUE.EQ.0) GOTO 400
2948
2949 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
2950 !.......................................................................
2951 ! CALCULATE SEDIMENATION
2952 ! THE NUMERICS HERE FOLLOW FROM REISNER ET AL. (1998)
2953 ! FALLOUT TERMS ARE CALCULATED ON SPLIT TIME STEPS TO ENSURE NUMERICAL
2954 ! STABILITY, I.E. COURANT# < 1
2955
2956 !.......................................................................
2957
2958 NSTEP = 1
2959
2960 DO K = KTS,KTE
2961
2962 DUMI(K) = QI3D(K)+QI3DTEN(K)*DT
2963 DUMQS(K) = QNI3D(K)+QNI3DTEN(K)*DT
2964 DUMR(K) = QR3D(K)+QR3DTEN(K)*DT
2965 DUMFNI(K) = NI3D(K)+NI3DTEN(K)*DT
2966 DUMFNS(K) = NS3D(K)+NS3DTEN(K)*DT
2967 DUMFNR(K) = NR3D(K)+NR3DTEN(K)*DT
2968 DUMC(K) = QC3D(K)+QC3DTEN(K)*DT
2969 DUMFNC(K) = NC3D(K)+NC3DTEN(K)*DT
2970 DUMG(K) = QG3D(K)+QG3DTEN(K)*DT
2971 DUMFNG(K) = NG3D(K)+NG3DTEN(K)*DT
2972
2973 ! SWITCH FOR CONSTANT DROPLET NUMBER
2974 ! IF (INUM.EQ.1) THEN
2975 DUMFNC(K) = NC3D(K)
2976 ! END IF
2977
2978 ! GET DUMMY LAMDA FOR SEDIMENTATION CALCULATIONS
2979
2980 ! MAKE SURE NUMBER CONCENTRATIONS ARE POSITIVE
2981 DUMFNI(K) = MAX(0.,DUMFNI(K))
2982 DUMFNS(K) = MAX(0.,DUMFNS(K))
2983 DUMFNC(K) = MAX(0.,DUMFNC(K))
2984 DUMFNR(K) = MAX(0.,DUMFNR(K))
2985 DUMFNG(K) = MAX(0.,DUMFNG(K))
2986
2987 !......................................................................
2988 ! CLOUD ICE
2989
2990 IF (DUMI(K).GE.QSMALL) THEN
2991 DLAMI = (CONS12*DUMFNI(K)/DUMI(K))**(1./DI)
2992 DLAMI=MAX(DLAMI,LAMMINI)
2993 DLAMI=MIN(DLAMI,LAMMAXI)
2994 END IF
2995 !......................................................................
2996 ! RAIN
2997
2998 IF (DUMR(K).GE.QSMALL) THEN
2999 DLAMR = (PI*RHOW*DUMFNR(K)/DUMR(K))**(1./3.)
3000 DLAMR=MAX(DLAMR,LAMMINR)
3001 DLAMR=MIN(DLAMR,LAMMAXR)
3002 END IF
3003 !......................................................................
3004 ! CLOUD DROPLETS
3005
3006 IF (DUMC(K).GE.QSMALL) THEN
3007 DUM = PRES(K)/(287.15*T3D(K))
3008 PGAM(K)=0.0005714*(NC3D(K)/1.E6/DUM)+0.2714
3009 PGAM(K)=1./(PGAM(K)**2)-1.
3010 PGAM(K)=MAX(PGAM(K),2.)
3011 PGAM(K)=MIN(PGAM(K),10.)
3012
3013 DLAMC = (CONS26*DUMFNC(K)*GAMMA(PGAM(K)+4.)/(DUMC(K)*GAMMA(PGAM(K)+1.)))**(1./3.)
3014 LAMMIN = (PGAM(K)+1.)/60.E-6
3015 LAMMAX = (PGAM(K)+1.)/1.E-6
3016 DLAMC=MAX(DLAMC,LAMMIN)
3017 DLAMC=MIN(DLAMC,LAMMAX)
3018 END IF
3019 !......................................................................
3020 ! SNOW
3021
3022 IF (DUMQS(K).GE.QSMALL) THEN
3023 DLAMS = (CONS1*DUMFNS(K)/ DUMQS(K))**(1./DS)
3024 DLAMS=MAX(DLAMS,LAMMINS)
3025 DLAMS=MIN(DLAMS,LAMMAXS)
3026 END IF
3027 !......................................................................
3028 ! GRAUPEL
3029
3030 IF (DUMG(K).GE.QSMALL) THEN
3031 DLAMG = (CONS2*DUMFNG(K)/ DUMG(K))**(1./DG)
3032 DLAMG=MAX(DLAMG,LAMMING)
3033 DLAMG=MIN(DLAMG,LAMMAXG)
3034 END IF
3035
3036 !......................................................................
3037 ! CALCULATE NUMBER-WEIGHTED AND MASS-WEIGHTED TERMINAL FALL SPEEDS
3038
3039 ! CLOUD WATER
3040
3041 IF (DUMC(K).GE.QSMALL) THEN
3042 UNC = ACN(K)*GAMMA(1.+BC+PGAM(K))/ (DLAMC**BC*GAMMA(PGAM(K)+1.))
3043 UMC = ACN(K)*GAMMA(4.+BC+PGAM(K))/ (DLAMC**BC*GAMMA(PGAM(K)+4.))
3044 ELSE
3045 UMC = 0.
3046 UNC = 0.
3047 END IF
3048
3049 IF (DUMI(K).GE.QSMALL) THEN
3050 UNI = AIN(K)*CONS27/DLAMI**BI
3051 UMI = AIN(K)*CONS28/(DLAMI**BI)
3052 ELSE
3053 UMI = 0.
3054 UNI = 0.
3055 END IF
3056
3057 IF (DUMR(K).GE.QSMALL) THEN
3058 UNR = ARN(K)*CONS6/DLAMR**BR
3059 UMR = ARN(K)*CONS4/(DLAMR**BR)
3060 ELSE
3061 UMR = 0.
3062 UNR = 0.
3063 END IF
3064
3065 IF (DUMQS(K).GE.QSMALL) THEN
3066 UMS = ASN(K)*CONS3/(DLAMS**BS)
3067 UNS = ASN(K)*CONS5/DLAMS**BS
3068 ELSE
3069 UMS = 0.
3070 UNS = 0.
3071 END IF
3072
3073 IF (DUMG(K).GE.QSMALL) THEN
3074 UMG = AGN(K)*CONS7/(DLAMG**BG)
3075 UNG = AGN(K)*CONS8/DLAMG**BG
3076 ELSE
3077 UMG = 0.
3078 UNG = 0.
3079 END IF
3080
3081 ! SET REALISTIC LIMITS ON FALLSPEED
3082
3083 UMS=MIN(UMS,1.2)
3084 UNS=MIN(UNS,1.2)
3085 UMI=MIN(UMI,1.2)
3086 UNI=MIN(UNI,1.2)
3087 UMR=MIN(UMR,9.1)
3088 UNR=MIN(UNR,9.1)
3089 UMG=MIN(UMG,20.)
3090 UNG=MIN(UNG,20.)
3091
3092 FR(K) = UMR
3093 FI(K) = UMI
3094 FNI(K) = UNI
3095 FS(K) = UMS
3096 FNS(K) = UNS
3097 FNR(K) = UNR
3098 FC(K) = UMC
3099 FNC(K) = UNC
3100 FG(K) = UMG
3101 FNG(K) = UNG
3102
3103 ! CALCULATE NUMBER OF SPLIT TIME STEPS
3104
3105 RGVM = MAX(FR(K),FI(K),FS(K),FC(K),FNI(K),FNR(K),FNS(K),FNC(K),FG(K),FNG(K))
3106 ! VVT CHANGED IFIX -> INT (GENERIC FUNCTION)
3107 NSTEP = MAX(INT(RGVM*DT/DZQ(K)+1.),NSTEP)
3108
3109 ! MULTIPLY VARIABLES BY RHO
3110 DUMR(k) = DUMR(k)*RHO(K)
3111 DUMI(k) = DUMI(k)*RHO(K)
3112 DUMFNI(k) = DUMFNI(K)*RHO(K)
3113 DUMQS(k) = DUMQS(K)*RHO(K)
3114 DUMFNS(k) = DUMFNS(K)*RHO(K)
3115 DUMFNR(k) = DUMFNR(K)*RHO(K)
3116 DUMC(k) = DUMC(K)*RHO(K)
3117 DUMFNC(k) = DUMFNC(K)*RHO(K)
3118 DUMG(k) = DUMG(K)*RHO(K)
3119 DUMFNG(k) = DUMFNG(K)*RHO(K)
3120
3121 END DO
3122
3123 DO N = 1,NSTEP
3124
3125 DO K = KTS,KTE
3126 FALOUTR(K) = FR(K)*DUMR(K)
3127 FALOUTI(K) = FI(K)*DUMI(K)
3128 FALOUTNI(K) = FNI(K)*DUMFNI(K)
3129 FALOUTS(K) = FS(K)*DUMQS(K)
3130 FALOUTNS(K) = FNS(K)*DUMFNS(K)
3131 FALOUTNR(K) = FNR(K)*DUMFNR(K)
3132 FALOUTC(K) = FC(K)*DUMC(K)
3133 FALOUTNC(K) = FNC(K)*DUMFNC(K)
3134 FALOUTG(K) = FG(K)*DUMG(K)
3135 FALOUTNG(K) = FNG(K)*DUMFNG(K)
3136 END DO
3137
3138 ! TOP OF MODEL
3139
3140 K = KTE
3141 FALTNDR = FALOUTR(K)/DZQ(k)
3142 FALTNDI = FALOUTI(K)/DZQ(k)
3143 FALTNDNI = FALOUTNI(K)/DZQ(k)
3144 FALTNDS = FALOUTS(K)/DZQ(k)
3145 FALTNDNS = FALOUTNS(K)/DZQ(k)
3146 FALTNDNR = FALOUTNR(K)/DZQ(k)
3147 FALTNDC = FALOUTC(K)/DZQ(k)
3148 FALTNDNC = FALOUTNC(K)/DZQ(k)
3149 FALTNDG = FALOUTG(K)/DZQ(k)
3150 FALTNDNG = FALOUTNG(K)/DZQ(k)
3151 ! ADD FALLOUT TERMS TO EULERIAN TENDENCIES
3152
3153 QRSTEN(K) = QRSTEN(K)-FALTNDR/NSTEP/RHO(k)
3154 QISTEN(K) = QISTEN(K)-FALTNDI/NSTEP/RHO(k)
3155 NI3DTEN(K) = NI3DTEN(K)-FALTNDNI/NSTEP/RHO(k)
3156 QNISTEN(K) = QNISTEN(K)-FALTNDS/NSTEP/RHO(k)
3157 NS3DTEN(K) = NS3DTEN(K)-FALTNDNS/NSTEP/RHO(k)
3158 NR3DTEN(K) = NR3DTEN(K)-FALTNDNR/NSTEP/RHO(k)
3159 QCSTEN(K) = QCSTEN(K)-FALTNDC/NSTEP/RHO(k)
3160 NC3DTEN(K) = NC3DTEN(K)-FALTNDNC/NSTEP/RHO(k)
3161 QGSTEN(K) = QGSTEN(K)-FALTNDG/NSTEP/RHO(k)
3162 NG3DTEN(K) = NG3DTEN(K)-FALTNDNG/NSTEP/RHO(k)
3163
3164 DUMR(K) = DUMR(K)-FALTNDR*DT/NSTEP
3165 DUMI(K) = DUMI(K)-FALTNDI*DT/NSTEP
3166 DUMFNI(K) = DUMFNI(K)-FALTNDNI*DT/NSTEP
3167 DUMQS(K) = DUMQS(K)-FALTNDS*DT/NSTEP
3168 DUMFNS(K) = DUMFNS(K)-FALTNDNS*DT/NSTEP
3169 DUMFNR(K) = DUMFNR(K)-FALTNDNR*DT/NSTEP
3170 DUMC(K) = DUMC(K)-FALTNDC*DT/NSTEP
3171 DUMFNC(K) = DUMFNC(K)-FALTNDNC*DT/NSTEP
3172 DUMG(K) = DUMG(K)-FALTNDG*DT/NSTEP
3173 DUMFNG(K) = DUMFNG(K)-FALTNDNG*DT/NSTEP
3174
3175 DO K = KTE-1,KTS,-1
3176 FALTNDR = (FALOUTR(K+1)-FALOUTR(K))/DZQ(K)
3177 FALTNDI = (FALOUTI(K+1)-FALOUTI(K))/DZQ(K)
3178 FALTNDNI = (FALOUTNI(K+1)-FALOUTNI(K))/DZQ(K)
3179 FALTNDS = (FALOUTS(K+1)-FALOUTS(K))/DZQ(K)
3180 FALTNDNS = (FALOUTNS(K+1)-FALOUTNS(K))/DZQ(K)
3181 FALTNDNR = (FALOUTNR(K+1)-FALOUTNR(K))/DZQ(K)
3182 FALTNDC = (FALOUTC(K+1)-FALOUTC(K))/DZQ(K)
3183 FALTNDNC = (FALOUTNC(K+1)-FALOUTNC(K))/DZQ(K)
3184 FALTNDG = (FALOUTG(K+1)-FALOUTG(K))/DZQ(K)
3185 FALTNDNG = (FALOUTNG(K+1)-FALOUTNG(K))/DZQ(K)
3186
3187 ! ADD FALLOUT TERMS TO EULERIAN TENDENCIES
3188
3189 QRSTEN(K) = QRSTEN(K)+FALTNDR/NSTEP/RHO(k)
3190 QISTEN(K) = QISTEN(K)+FALTNDI/NSTEP/RHO(k)
3191 NI3DTEN(K) = NI3DTEN(K)+FALTNDNI/NSTEP/RHO(k)
3192 QNISTEN(K) = QNISTEN(K)+FALTNDS/NSTEP/RHO(k)
3193 NS3DTEN(K) = NS3DTEN(K)+FALTNDNS/NSTEP/RHO(k)
3194 NR3DTEN(K) = NR3DTEN(K)+FALTNDNR/NSTEP/RHO(k)
3195 QCSTEN(K) = QCSTEN(K)+FALTNDC/NSTEP/RHO(k)
3196 NC3DTEN(K) = NC3DTEN(K)+FALTNDNC/NSTEP/RHO(k)
3197 QGSTEN(K) = QGSTEN(K)+FALTNDG/NSTEP/RHO(k)
3198 NG3DTEN(K) = NG3DTEN(K)+FALTNDNG/NSTEP/RHO(k)
3199
3200 DUMR(K) = DUMR(K)+FALTNDR*DT/NSTEP
3201 DUMI(K) = DUMI(K)+FALTNDI*DT/NSTEP
3202 DUMFNI(K) = DUMFNI(K)+FALTNDNI*DT/NSTEP
3203 DUMQS(K) = DUMQS(K)+FALTNDS*DT/NSTEP
3204 DUMFNS(K) = DUMFNS(K)+FALTNDNS*DT/NSTEP
3205 DUMFNR(K) = DUMFNR(K)+FALTNDNR*DT/NSTEP
3206 DUMC(K) = DUMC(K)+FALTNDC*DT/NSTEP
3207 DUMFNC(K) = DUMFNC(K)+FALTNDNC*DT/NSTEP
3208 DUMG(K) = DUMG(K)+FALTNDG*DT/NSTEP
3209 DUMFNG(K) = DUMFNG(K)+FALTNDNG*DT/NSTEP
3210
3211 END DO
3212
3213 ! GET PRECIPITATION AND SNOWFALL ACCUMULATION DURING THE TIME STEP
3214 ! FACTOR OF 1000 CONVERTS FROM M TO MM, BUT DIVISION BY DENSITY
3215 ! OF LIQUID WATER CANCELS THIS FACTOR OF 1000
3216
3217 PRECRT = PRECRT+(FALOUTR(KTS)+FALOUTC(KTS)+FALOUTS(KTS)+FALOUTI(KTS)+FALOUTG(KTS)) &
3218 *DT/NSTEP
3219 SNOWRT = SNOWRT+(FALOUTS(KTS)+FALOUTI(KTS)+FALOUTG(KTS))*DT/NSTEP
3220
3221 END DO
3222
3223 DO K=KTS,KTE
3224
3225 ! ADD ON SEDIMENTATION TENDENCIES FOR MIXING RATIO TO REST OF TENDENCIES
3226
3227 QR3DTEN(K)=QR3DTEN(K)+QRSTEN(K)
3228 QI3DTEN(K)=QI3DTEN(K)+QISTEN(K)
3229 QC3DTEN(K)=QC3DTEN(K)+QCSTEN(K)
3230 QG3DTEN(K)=QG3DTEN(K)+QGSTEN(K)
3231 QNI3DTEN(K)=QNI3DTEN(K)+QNISTEN(K)
3232
3233 ! PUT ALL CLOUD ICE IN SNOW CATEGORY IF MEAN DIAMETER EXCEEDS 2 * dcs
3234
3235 IF (QI3D(K).GE.QSMALL.AND.T3D(K).LT.273.15) THEN
3236 IF (1./LAMI(K).GE.2.*DCS) THEN
3237 QNI3DTEN(K) = QNI3DTEN(K)+QI3D(K)/DT+ QI3DTEN(K)
3238 NS3DTEN(K) = NS3DTEN(K)+NI3D(K)/DT+ NI3DTEN(K)
3239 QI3DTEN(K) = -QI3D(K)/DT
3240 NI3DTEN(K) = -NI3D(K)/DT
3241 END IF
3242 END IF
3243
3244 ! hm add tendencies here, then call sizeparameter
3245 ! to ensure consisitency between mixing ratio and number concentration
3246
3247 QC3D(k) = QC3D(k)+QC3DTEN(k)*DT
3248 QI3D(k) = QI3D(k)+QI3DTEN(k)*DT
3249 QNI3D(k) = QNI3D(k)+QNI3DTEN(k)*DT
3250 QR3D(k) = QR3D(k)+QR3DTEN(k)*DT
3251 ! NC3D(k) = NC3D(k)+NC3DTEN(k)*DT
3252 NI3D(k) = NI3D(k)+NI3DTEN(k)*DT
3253 NS3D(k) = NS3D(k)+NS3DTEN(k)*DT
3254 NR3D(k) = NR3D(k)+NR3DTEN(k)*DT
3255
3256 IF (IGRAUP.EQ.0) THEN
3257 QG3D(k) = QG3D(k)+QG3DTEN(k)*DT
3258 NG3D(k) = NG3D(k)+NG3DTEN(k)*DT
3259 END IF
3260
3261 ! ADD TEMPERATURE AND WATER VAPOR TENDENCIES FROM MICROPHYSICS
3262 T3D(K) = T3D(K)+T3DTEN(k)*DT
3263 QV3D(K) = QV3D(K)+QV3DTEN(k)*DT
3264
3265 ! SATURATION VAPOR PRESSURE AND MIXING RATIO
3266
3267 EVS(K) = POLYSVP(T3D(K),0) ! PA
3268 EIS(K) = POLYSVP(T3D(K),1) ! PA
3269
3270 ! MAKE SURE ICE SATURATION DOESN'T EXCEED WATER SAT. NEAR FREEZING
3271
3272 IF (EIS(K).GT.EVS(K)) EIS(K) = EVS(K)
3273
3274 QVS(K) = .622*EVS(K)/(PRES(K)-EVS(K))
3275 QVI(K) = .622*EIS(K)/(PRES(K)-EIS(K))
3276
3277 QVQVS(K) = QV3D(K)/QVS(K)
3278 QVQVSI(K) = QV3D(K)/QVI(K)
3279
3280 ! AT SUBSATURATION, REMOVE SMALL AMOUNTS OF CLOUD/PRECIP WATER
3281
3282 IF (QVQVS(K).LT.0.9) THEN
3283 IF (QR3D(K).LT.1.E-6) THEN
3284 QV3D(K)=QV3D(K)+QR3D(K)
3285 T3D(K)=T3D(K)-QR3D(K)*XXLV(K)/CPM(K)
3286 QR3D(K)=0.
3287 END IF
3288 IF (QC3D(K).LT.1.E-6) THEN
3289 QV3D(K)=QV3D(K)+QC3D(K)
3290 T3D(K)=T3D(K)-QC3D(K)*XXLV(K)/CPM(K)
3291 QC3D(K)=0.
3292 END IF
3293 END IF
3294
3295 IF (QVQVSI(K).LT.0.9) THEN
3296 IF (QI3D(K).LT.1.E-6) THEN
3297 QV3D(K)=QV3D(K)+QI3D(K)
3298 T3D(K)=T3D(K)-QI3D(K)*XXLS(K)/CPM(K)
3299 QI3D(K)=0.
3300 END IF
3301 IF (QNI3D(K).LT.1.E-6) THEN
3302 QV3D(K)=QV3D(K)+QNI3D(K)
3303 T3D(K)=T3D(K)-QNI3D(K)*XXLS(K)/CPM(K)
3304 QNI3D(K)=0.
3305 END IF
3306 IF (QG3D(K).LT.1.E-6) THEN
3307 QV3D(K)=QV3D(K)+QG3D(K)
3308 T3D(K)=T3D(K)-QG3D(K)*XXLS(K)/CPM(K)
3309 QG3D(K)=0.
3310 END IF
3311 END IF
3312
3313 !..................................................................
3314 ! IF MIXING RATIO < QSMALL SET MIXING RATIO AND NUMBER CONC TO ZERO
3315
3316 IF (QC3D(K).LT.QSMALL) THEN
3317 QC3D(K) = 0.
3318 NC3D(K) = 0.
3319 EFFC(K) = 0.
3320 END IF
3321 IF (QR3D(K).LT.QSMALL) THEN
3322 QR3D(K) = 0.
3323 NR3D(K) = 0.
3324 EFFR(K) = 0.
3325 END IF
3326 IF (QI3D(K).LT.QSMALL) THEN
3327 QI3D(K) = 0.
3328 NI3D(K) = 0.
3329 EFFI(K) = 0.
3330 END IF
3331 IF (QNI3D(K).LT.QSMALL) THEN
3332 QNI3D(K) = 0.
3333 NS3D(K) = 0.
3334 EFFS(K) = 0.
3335 END IF
3336 IF (QG3D(K).LT.QSMALL) THEN
3337 QG3D(K) = 0.
3338 NG3D(K) = 0.
3339 EFFG(K) = 0.
3340 END IF
3341
3342 !..................................
3343 ! IF THERE IS NO CLOUD/PRECIP WATER, THEN SKIP CALCULATIONS
3344
3345 IF (QC3D(K).LT.QSMALL.AND.QI3D(K).LT.QSMALL.AND.QNI3D(K).LT.QSMALL &
3346 .AND.QR3D(K).LT.QSMALL.AND.QG3D(K).LT.QSMALL) GOTO 500
3347
3348 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
3349 ! CALCULATE INSTANTANEOUS PROCESSES
3350
3351 ! ADD MELTING OF CLOUD ICE TO FORM RAIN
3352
3353 IF (QI3D(K).GE.QSMALL.AND.T3D(K).GE.273.15) THEN
3354 QR3D(K) = QR3D(K)+QI3D(K)
3355 T3D(K) = T3D(K)-QI3D(K)*XLF(K)/CPM(K)
3356 QI3D(K) = 0.
3357 NR3D(K) = NR3D(K)+NI3D(K)
3358 NI3D(K) = 0.
3359 END IF
3360
3361 ! ****SENSITIVITY - NO ICE
3362 IF (ILIQ.EQ.1) GOTO 778
3363
3364 ! HOMOGENEOUS FREEZING OF CLOUD WATER
3365
3366 IF (T3D(K).LE.233.15.AND.QC3D(K).GE.QSMALL) THEN
3367 QI3D(K)=QI3D(K)+QC3D(K)
3368 T3D(K)=T3D(K)+QC3D(K)*XLF(K)/CPM(K)
3369 QC3D(K)=0.
3370 NI3D(K)=NI3D(K)+NC3D(K)
3371 NC3D(K)=0.
3372 END IF
3373
3374 ! HOMOGENEOUS FREEZING OF RAIN
3375
3376 IF (IGRAUP.EQ.0) THEN
3377
3378 IF (T3D(K).LE.233.15.AND.QR3D(K).GE.QSMALL) THEN
3379 QG3D(K) = QG3D(K)+QR3D(K)
3380 T3D(K) = T3D(K)+QR3D(K)*XLF(K)/CPM(K)
3381 QR3D(K) = 0.
3382 NG3D(K) = NG3D(K)+ NR3D(K)
3383 NR3D(K) = 0.
3384 END IF
3385
3386 ELSE IF (IGRAUP.EQ.1) THEN
3387
3388 IF (T3D(K).LE.233.15.AND.QR3D(K).GE.QSMALL) THEN
3389 QNI3D(K) = QNI3D(K)+QR3D(K)
3390 T3D(K) = T3D(K)+QR3D(K)*XLF(K)/CPM(K)
3391 QR3D(K) = 0.
3392 NS3D(K) = NS3D(K)+NR3D(K)
3393 NR3D(K) = 0.
3394 END IF
3395
3396 END IF
3397
3398 778 CONTINUE
3399
3400 ! MAKE SURE NUMBER CONCENTRATIONS AREN'T NEGATIVE
3401
3402 NI3D(K) = MAX(0.,NI3D(K))
3403 NS3D(K) = MAX(0.,NS3D(K))
3404 NC3D(K) = MAX(0.,NC3D(K))
3405 NR3D(K) = MAX(0.,NR3D(K))
3406 NG3D(K) = MAX(0.,NG3D(K))
3407
3408 !......................................................................
3409 ! CLOUD ICE
3410
3411 IF (QI3D(K).GE.QSMALL) THEN
3412 LAMI(K) = (CONS12* &
3413 NI3D(K)/QI3D(K))**(1./DI)
3414
3415 ! CHECK FOR SLOPE
3416
3417 ! ADJUST VARS
3418
3419 IF (LAMI(K).LT.LAMMINI) THEN
3420
3421 LAMI(K) = LAMMINI
3422
3423 N0I(K) = LAMI(K)**4*QI3D(K)/CONS12
3424
3425 NI3D(K) = N0I(K)/LAMI(K)
3426 ELSE IF (LAMI(K).GT.LAMMAXI) THEN
3427 LAMI(K) = LAMMAXI
3428 N0I(K) = LAMI(K)**4*QI3D(K)/CONS12
3429
3430 NI3D(K) = N0I(K)/LAMI(K)
3431 END IF
3432 END IF
3433
3434 !......................................................................
3435 ! RAIN
3436
3437 IF (QR3D(K).GE.QSMALL) THEN
3438 LAMR(K) = (PI*RHOW*NR3D(K)/QR3D(K))**(1./3.)
3439
3440 ! CHECK FOR SLOPE
3441
3442 ! ADJUST VARS
3443
3444 IF (LAMR(K).LT.LAMMINR) THEN
3445
3446 LAMR(K) = LAMMINR
3447
3448 N0RR(K) = LAMR(K)**4*QR3D(K)/(PI*RHOW)
3449
3450 NR3D(K) = N0RR(K)/LAMR(K)
3451 ELSE IF (LAMR(K).GT.LAMMAXR) THEN
3452 LAMR(K) = LAMMAXR
3453 N0RR(K) = LAMR(K)**4*QR3D(K)/(PI*RHOW)
3454
3455 NR3D(K) = N0RR(K)/LAMR(K)
3456 END IF
3457
3458 END IF
3459
3460 !......................................................................
3461 ! CLOUD DROPLETS
3462
3463 ! MARTIN ET AL. (1994) FORMULA FOR PGAM
3464
3465 IF (QC3D(K).GE.QSMALL) THEN
3466
3467 DUM = PRES(K)/(287.15*T3D(K))
3468 PGAM(K)=0.0005714*(NC3D(K)/1.E6/DUM)+0.2714
3469 PGAM(K)=1./(PGAM(K)**2)-1.
3470 PGAM(K)=MAX(PGAM(K),2.)
3471 PGAM(K)=MIN(PGAM(K),10.)
3472
3473 ! CALCULATE LAMC
3474
3475 LAMC(K) = (CONS26*NC3D(K)*GAMMA(PGAM(K)+4.)/ &
3476 (QC3D(K)*GAMMA(PGAM(K)+1.)))**(1./3.)
3477
3478 ! LAMMIN, 60 MICRON DIAMETER
3479 ! LAMMAX, 1 MICRON
3480
3481 LAMMIN = (PGAM(K)+1.)/60.E-6
3482 LAMMAX = (PGAM(K)+1.)/1.E-6
3483
3484 IF (LAMC(K).LT.LAMMIN) THEN
3485 LAMC(K) = LAMMIN
3486 NC3D(K) = EXP(3.*LOG(LAMC(K))+LOG(QC3D(K))+ &
3487 LOG(GAMMA(PGAM(K)+1.))-LOG(GAMMA(PGAM(K)+4.)))/CONS26
3488
3489 ELSE IF (LAMC(K).GT.LAMMAX) THEN
3490 LAMC(K) = LAMMAX
3491 NC3D(K) = EXP(3.*LOG(LAMC(K))+LOG(QC3D(K))+ &
3492 LOG(GAMMA(PGAM(K)+1.))-LOG(GAMMA(PGAM(K)+4.)))/CONS26
3493
3494 END IF
3495
3496 END IF
3497
3498 !......................................................................
3499 ! SNOW
3500
3501 IF (QNI3D(K).GE.QSMALL) THEN
3502 LAMS(K) = (CONS1*NS3D(K)/QNI3D(K))**(1./DS)
3503
3504 ! CHECK FOR SLOPE
3505
3506 ! ADJUST VARS
3507
3508 IF (LAMS(K).LT.LAMMINS) THEN
3509 LAMS(K) = LAMMINS
3510 N0S(K) = LAMS(K)**4*QNI3D(K)/CONS1
3511
3512 NS3D(K) = N0S(K)/LAMS(K)
3513
3514 ELSE IF (LAMS(K).GT.LAMMAXS) THEN
3515
3516 LAMS(K) = LAMMAXS
3517 N0S(K) = LAMS(K)**4*QNI3D(K)/CONS1
3518 NS3D(K) = N0S(K)/LAMS(K)
3519 END IF
3520
3521 END IF
3522
3523 !......................................................................
3524 ! GRAUPEL
3525
3526 IF (QG3D(K).GE.QSMALL) THEN
3527 LAMG(K) = (CONS2*NG3D(K)/QG3D(K))**(1./DG)
3528
3529 ! CHECK FOR SLOPE
3530
3531 ! ADJUST VARS
3532
3533 IF (LAMG(K).LT.LAMMING) THEN
3534 LAMG(K) = LAMMING
3535 N0G(K) = LAMG(K)**4*QG3D(K)/CONS2
3536
3537 NG3D(K) = N0G(K)/LAMG(K)
3538
3539 ELSE IF (LAMG(K).GT.LAMMAXG) THEN
3540
3541 LAMG(K) = LAMMAXG
3542 N0G(K) = LAMG(K)**4*QG3D(K)/CONS2
3543
3544 NG3D(K) = N0G(K)/LAMG(K)
3545 END IF
3546
3547 END IF
3548
3549 500 CONTINUE
3550
3551 ! CALCULATE EFFECTIVE RADIUS
3552
3553 IF (QI3D(K).GE.QSMALL) THEN
3554 EFFI(K) = 3./LAMI(K)/2.*1.E6
3555 ELSE
3556 EFFI(K) = 25.
3557 END IF
3558
3559 IF (QNI3D(K).GE.QSMALL) THEN
3560 EFFS(K) = 3./LAMS(K)/2.*1.E6
3561 ELSE
3562 EFFS(K) = 25.
3563 END IF
3564
3565 IF (QR3D(K).GE.QSMALL) THEN
3566 EFFR(K) = 3./LAMR(K)/2.*1.E6
3567 ELSE
3568 EFFR(K) = 25.
3569 END IF
3570
3571 IF (QC3D(K).GE.QSMALL) THEN
3572 EFFC(K) = GAMMA(PGAM(K)+4.)/ &
3573 GAMMA(PGAM(K)+3.)/LAMC(K)/2.*1.E6
3574 ELSE
3575 EFFC(K) = 25.
3576 END IF
3577
3578 IF (QG3D(K).GE.QSMALL) THEN
3579 EFFG(K) = 3./LAMG(K)/2.*1.E6
3580 ELSE
3581 EFFG(K) = 25.
3582 END IF
3583
3584 ! HM ADD 1/10/06, ADD UPPER BOUND ON ICE NUMBER, THIS IS NEEDED
3585 ! TO PREVENT VERY LARGE ICE NUMBER DUE TO HOMOGENEOUS FREEZING
3586 ! OF DROPLETS, ESPECIALLY WHEN INUM = 1, SET MAX AT 10 CM-3
3587 NI3D(K) = MIN(NI3D(K),10.E6/RHO(K))
3588 ! ADD BOUND ON DROPLET NUMBER - CANNOT EXCEED AEROSOL CONCENTRATION
3589 IF (INUM.EQ.0.AND.IACT.EQ.2) THEN
3590 NC3D(K) = MIN(NC3D(K),(NANEW1+NANEW2)/RHO(K))
3591 END IF
3592 ! SWITCH FOR CONSTANT DROPLET NUMBER
3593 ! IF (INUM.EQ.1) THEN
3594 ! CHANGE NDCNST FROM CM-3 TO KG-1
3595 NC3D(K) = NDCNST*1.E6/RHO(K)
3596 ! END IF
3597
3598 END DO !!! K LOOP
3599
3600 400 CONTINUE
3601
3602 ! ALL DONE !!!!!!!!!!!
3603 RETURN
3604 END SUBROUTINE MORR_TWO_MOMENT_MICRO
3605
3606 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
3607
3608 REAL FUNCTION POLYSVP (T,TYPE)
3609
3610 !-------------------------------------------
3611
3612 ! COMPUTE SATURATION VAPOR PRESSURE
3613
3614 ! POLYSVP RETURNED IN UNITS OF PA.
3615 ! T IS INPUT IN UNITS OF K.
3616 ! TYPE REFERS TO SATURATION WITH RESPECT TO LIQUID (0) OR ICE (1)
3617
3618 ! REPLACE GOFF-GRATCH WITH FASTER FORMULATION FROM MARAT KHROUTDINOV
3619
3620 IMPLICIT NONE
3621
3622 REAL DUM
3623 REAL T
3624 INTEGER TYPE
3625 ! ice
3626 real a0i,a1i,a2i,a3i,a4i,a5i,a6i,a7i,a8i
3627 data a0i,a1i,a2i,a3i,a4i,a5i,a6i,a7i,a8i /&
3628 6.11147274, 0.503160820, 0.188439774e-1, &
3629 0.420895665e-3, 0.615021634e-5,0.602588177e-7, &
3630 0.385852041e-9, 0.146898966e-11, 0.252751365e-14/
3631
3632 ! liquid
3633 real a0,a1,a2,a3,a4,a5,a6,a7,a8
3634 data a0,a1,a2,a3,a4,a5,a6,a7,a8 /&
3635 6.105851, 0.4440316, 0.1430341e-1, &
3636 0.2641412e-3, 0.2995057e-5, 0.2031998e-7, &
3637 0.6936113e-10, 0.2564861e-13,-0.3704404e-15/
3638 real dt
3639
3640 ! ICE
3641
3642 IF (TYPE.EQ.1) THEN
3643
3644 ! POLYSVP = 10.**(-9.09718*(273.16/T-1.)-3.56654* &
3645 ! LOG10(273.16/T)+0.876793*(1.-T/273.16)+ &
3646 ! LOG10(6.1071))*100.
3647
3648
3649 dt = max(-80.,t-273.16)
3650 polysvp = a0i + dt*(a1i+dt*(a2i+dt*(a3i+dt*(a4i+dt*(a5i+dt*(a6i+dt*(a7i+a8i*dt)))))))
3651 polysvp = polysvp*100.
3652
3653 END IF
3654
3655 ! LIQUID
3656
3657 IF (TYPE.EQ.0) THEN
3658
3659 dt = max(-80.,t-273.16)
3660 polysvp = a0 + dt*(a1+dt*(a2+dt*(a3+dt*(a4+dt*(a5+dt*(a6+dt*(a7+a8*dt)))))))
3661 polysvp = polysvp*100.
3662
3663 ! POLYSVP = 10.**(-7.90298*(373.16/T-1.)+ &
3664 ! 5.02808*LOG10(373.16/T)- &
3665 ! 1.3816E-7*(10**(11.344*(1.-T/373.16))-1.)+ &
3666 ! 8.1328E-3*(10**(-3.49149*(373.16/T-1.))-1.)+ &
3667 ! LOG10(1013.246))*100.
3668
3669 END IF
3670
3671
3672 END FUNCTION POLYSVP
3673
3674 !------------------------------------------------------------------------------
3675
3676 REAL FUNCTION GAMMA(X)
3677 !----------------------------------------------------------------------
3678 !
3679 ! THIS ROUTINE CALCULATES THE GAMMA FUNCTION FOR A REAL ARGUMENT X.
3680 ! COMPUTATION IS BASED ON AN ALGORITHM OUTLINED IN REFERENCE 1.
3681 ! THE PROGRAM USES RATIONAL FUNCTIONS THAT APPROXIMATE THE GAMMA
3682 ! FUNCTION TO AT LEAST 20 SIGNIFICANT DECIMAL DIGITS. COEFFICIENTS
3683 ! FOR THE APPROXIMATION OVER THE INTERVAL (1,2) ARE UNPUBLISHED.
3684 ! THOSE FOR THE APPROXIMATION FOR X .GE. 12 ARE FROM REFERENCE 2.
3685 ! THE ACCURACY ACHIEVED DEPENDS ON THE ARITHMETIC SYSTEM, THE
3686 ! COMPILER, THE INTRINSIC FUNCTIONS, AND PROPER SELECTION OF THE
3687 ! MACHINE-DEPENDENT CONSTANTS.
3688 !
3689 !
3690 !*******************************************************************
3691 !*******************************************************************
3692 !
3693 ! EXPLANATION OF MACHINE-DEPENDENT CONSTANTS
3694 !
3695 ! BETA - RADIX FOR THE FLOATING-POINT REPRESENTATION
3696 ! MAXEXP - THE SMALLEST POSITIVE POWER OF BETA THAT OVERFLOWS
3697 ! XBIG - THE LARGEST ARGUMENT FOR WHICH GAMMA(X) IS REPRESENTABLE
3698 ! IN THE MACHINE, I.E., THE SOLUTION TO THE EQUATION
3699 ! GAMMA(XBIG) = BETA**MAXEXP
3700 ! XINF - THE LARGEST MACHINE REPRESENTABLE FLOATING-POINT NUMBER;
3701 ! APPROXIMATELY BETA**MAXEXP
3702 ! EPS - THE SMALLEST POSITIVE FLOATING-POINT NUMBER SUCH THAT
3703 ! 1.0+EPS .GT. 1.0
3704 ! XMININ - THE SMALLEST POSITIVE FLOATING-POINT NUMBER SUCH THAT
3705 ! 1/XMININ IS MACHINE REPRESENTABLE
3706 !
3707 ! APPROXIMATE VALUES FOR SOME IMPORTANT MACHINES ARE:
3708 !
3709 ! BETA MAXEXP XBIG
3710 !
3711 ! CRAY-1 (S.P.) 2 8191 966.961
3712 ! CYBER 180/855
3713 ! UNDER NOS (S.P.) 2 1070 177.803
3714 ! IEEE (IBM/XT,
3715 ! SUN, ETC.) (S.P.) 2 128 35.040
3716 ! IEEE (IBM/XT,
3717 ! SUN, ETC.) (D.P.) 2 1024 171.624
3718 ! IBM 3033 (D.P.) 16 63 57.574
3719 ! VAX D-FORMAT (D.P.) 2 127 34.844
3720 ! VAX G-FORMAT (D.P.) 2 1023 171.489
3721 !
3722 ! XINF EPS XMININ
3723 !
3724 ! CRAY-1 (S.P.) 5.45E+2465 7.11E-15 1.84E-2466
3725 ! CYBER 180/855
3726 ! UNDER NOS (S.P.) 1.26E+322 3.55E-15 3.14E-294
3727 ! IEEE (IBM/XT,
3728 ! SUN, ETC.) (S.P.) 3.40E+38 1.19E-7 1.18E-38
3729 ! IEEE (IBM/XT,
3730 ! SUN, ETC.) (D.P.) 1.79D+308 2.22D-16 2.23D-308
3731 ! IBM 3033 (D.P.) 7.23D+75 2.22D-16 1.39D-76
3732 ! VAX D-FORMAT (D.P.) 1.70D+38 1.39D-17 5.88D-39
3733 ! VAX G-FORMAT (D.P.) 8.98D+307 1.11D-16 1.12D-308
3734 !
3735 !*******************************************************************
3736 !*******************************************************************
3737 !
3738 ! ERROR RETURNS
3739 !
3740 ! THE PROGRAM RETURNS THE VALUE XINF FOR SINGULARITIES OR
3741 ! WHEN OVERFLOW WOULD OCCUR. THE COMPUTATION IS BELIEVED
3742 ! TO BE FREE OF UNDERFLOW AND OVERFLOW.
3743 !
3744 !
3745 ! INTRINSIC FUNCTIONS REQUIRED ARE:
3746 !
3747 ! INT, DBLE, EXP, LOG, REAL, SIN
3748 !
3749 !
3750 ! REFERENCES: AN OVERVIEW OF SOFTWARE DEVELOPMENT FOR SPECIAL
3751 ! FUNCTIONS W. J. CODY, LECTURE NOTES IN MATHEMATICS,
3752 ! 506, NUMERICAL ANALYSIS DUNDEE, 1975, G. A. WATSON
3753 ! (ED.), SPRINGER VERLAG, BERLIN, 1976.
3754 !
3755 ! COMPUTER APPROXIMATIONS, HART, ET. AL., WILEY AND
3756 ! SONS, NEW YORK, 1968.
3757 !
3758 ! LATEST MODIFICATION: OCTOBER 12, 1989
3759 !
3760 ! AUTHORS: W. J. CODY AND L. STOLTZ
3761 ! APPLIED MATHEMATICS DIVISION
3762 ! ARGONNE NATIONAL LABORATORY
3763 ! ARGONNE, IL 60439
3764 !
3765 !----------------------------------------------------------------------
3766 implicit none
3767 INTEGER I,N
3768 LOGICAL PARITY
3769 REAL &
3770 CONV,EPS,FACT,HALF,ONE,RES,SUM,TWELVE, &
3771 TWO,X,XBIG,XDEN,XINF,XMININ,XNUM,Y,Y1,YSQ,Z,ZERO
3772 REAL, DIMENSION(7) :: C
3773 REAL, DIMENSION(8) :: P
3774 REAL, DIMENSION(8) :: Q
3775 !----------------------------------------------------------------------
3776 ! MATHEMATICAL CONSTANTS
3777 !----------------------------------------------------------------------
3778 DATA ONE,HALF,TWELVE,TWO,ZERO/1.0E0,0.5E0,12.0E0,2.0E0,0.0E0/
3779
3780
3781 !----------------------------------------------------------------------
3782 ! MACHINE DEPENDENT PARAMETERS
3783 !----------------------------------------------------------------------
3784 DATA XBIG,XMININ,EPS/35.040E0,1.18E-38,1.19E-7/,XINF/3.4E38/
3785 !----------------------------------------------------------------------
3786 ! NUMERATOR AND DENOMINATOR COEFFICIENTS FOR RATIONAL MINIMAX
3787 ! APPROXIMATION OVER (1,2).
3788 !----------------------------------------------------------------------
3789 DATA P/-1.71618513886549492533811E+0,2.47656508055759199108314E+1, &
3790 -3.79804256470945635097577E+2,6.29331155312818442661052E+2, &
3791 8.66966202790413211295064E+2,-3.14512729688483675254357E+4, &
3792 -3.61444134186911729807069E+4,6.64561438202405440627855E+4/
3793 DATA Q/-3.08402300119738975254353E+1,3.15350626979604161529144E+2, &
3794 -1.01515636749021914166146E+3,-3.10777167157231109440444E+3, &
3795 2.25381184209801510330112E+4,4.75584627752788110767815E+3, &
3796 -1.34659959864969306392456E+5,-1.15132259675553483497211E+5/
3797 !----------------------------------------------------------------------
3798 ! COEFFICIENTS FOR MINIMAX APPROXIMATION OVER (12, INF).
3799 !----------------------------------------------------------------------
3800 DATA C/-1.910444077728E-03,8.4171387781295E-04, &
3801 -5.952379913043012E-04,7.93650793500350248E-04, &
3802 -2.777777777777681622553E-03,8.333333333333333331554247E-02, &
3803 5.7083835261E-03/
3804 !----------------------------------------------------------------------
3805 ! STATEMENT FUNCTIONS FOR CONVERSION BETWEEN INTEGER AND FLOAT
3806 !----------------------------------------------------------------------
3807 CONV(I) = REAL(I)
3808 PARITY=.FALSE.
3809 FACT=ONE
3810 N=0
3811 Y=X
3812 IF(Y.LE.ZERO)THEN
3813 !----------------------------------------------------------------------
3814 ! ARGUMENT IS NEGATIVE
3815 !----------------------------------------------------------------------
3816 Y=-X
3817 Y1=AINT(Y)
3818 RES=Y-Y1
3819 IF(RES.NE.ZERO)THEN
3820 IF(Y1.NE.AINT(Y1*HALF)*TWO)PARITY=.TRUE.
3821 FACT=-PI/SIN(PI*RES)
3822 Y=Y+ONE
3823 ELSE
3824 RES=XINF
3825 GOTO 900
3826 ENDIF
3827 ENDIF
3828 !----------------------------------------------------------------------
3829 ! ARGUMENT IS POSITIVE
3830 !----------------------------------------------------------------------
3831 IF(Y.LT.EPS)THEN
3832 !----------------------------------------------------------------------
3833 ! ARGUMENT .LT. EPS
3834 !----------------------------------------------------------------------
3835 IF(Y.GE.XMININ)THEN
3836 RES=ONE/Y
3837 ELSE
3838 RES=XINF
3839 GOTO 900
3840 ENDIF
3841 ELSEIF(Y.LT.TWELVE)THEN
3842 Y1=Y
3843 IF(Y.LT.ONE)THEN
3844 !----------------------------------------------------------------------
3845 ! 0.0 .LT. ARGUMENT .LT. 1.0
3846 !----------------------------------------------------------------------
3847 Z=Y
3848 Y=Y+ONE
3849 ELSE
3850 !----------------------------------------------------------------------
3851 ! 1.0 .LT. ARGUMENT .LT. 12.0, REDUCE ARGUMENT IF NECESSARY
3852 !----------------------------------------------------------------------
3853 N=INT(Y)-1
3854 Y=Y-CONV(N)
3855 Z=Y-ONE
3856 ENDIF
3857 !----------------------------------------------------------------------
3858 ! EVALUATE APPROXIMATION FOR 1.0 .LT. ARGUMENT .LT. 2.0
3859 !----------------------------------------------------------------------
3860 XNUM=ZERO
3861 XDEN=ONE
3862 DO I=1,8
3863 XNUM=(XNUM+P(I))*Z
3864 XDEN=XDEN*Z+Q(I)
3865 END DO
3866 RES=XNUM/XDEN+ONE
3867 IF(Y1.LT.Y)THEN
3868 !----------------------------------------------------------------------
3869 ! ADJUST RESULT FOR CASE 0.0 .LT. ARGUMENT .LT. 1.0
3870 !----------------------------------------------------------------------
3871 RES=RES/Y1
3872 ELSEIF(Y1.GT.Y)THEN
3873 !----------------------------------------------------------------------
3874 ! ADJUST RESULT FOR CASE 2.0 .LT. ARGUMENT .LT. 12.0
3875 !----------------------------------------------------------------------
3876 DO I=1,N
3877 RES=RES*Y
3878 Y=Y+ONE
3879 END DO
3880 ENDIF
3881 ELSE
3882 !----------------------------------------------------------------------
3883 ! EVALUATE FOR ARGUMENT .GE. 12.0,
3884 !----------------------------------------------------------------------
3885 IF(Y.LE.XBIG)THEN
3886 YSQ=Y*Y
3887 SUM=C(7)
3888 DO I=1,6
3889 SUM=SUM/YSQ+C(I)
3890 END DO
3891 SUM=SUM/Y-Y+SQRTPI
3892 SUM=SUM+(Y-HALF)*LOG(Y)
3893 RES=EXP(SUM)
3894 ELSE
3895 RES=XINF
3896 GOTO 900
3897 ENDIF
3898 ENDIF
3899 !----------------------------------------------------------------------
3900 ! FINAL ADJUSTMENTS AND RETURN
3901 !----------------------------------------------------------------------
3902 IF(PARITY)RES=-RES
3903 IF(FACT.NE.ONE)RES=FACT/RES
3904 900 GAMMA=RES
3905 RETURN
3906 ! ---------- LAST LINE OF GAMMA ----------
3907 END FUNCTION GAMMA
3908
3909
3910 REAL FUNCTION DERF1(X)
3911 IMPLICIT NONE
3912 REAL X
3913 REAL, DIMENSION(0 : 64) :: A, B
3914 REAL W,T,Y
3915 INTEGER K,I
3916 DATA A/ &
3917 0.00000000005958930743E0, -0.00000000113739022964E0, &
3918 0.00000001466005199839E0, -0.00000016350354461960E0, &
3919 0.00000164610044809620E0, -0.00001492559551950604E0, &
3920 0.00012055331122299265E0, -0.00085483269811296660E0, &
3921 0.00522397762482322257E0, -0.02686617064507733420E0, &
3922 0.11283791670954881569E0, -0.37612638903183748117E0, &
3923 1.12837916709551257377E0, &
3924 0.00000000002372510631E0, -0.00000000045493253732E0, &
3925 0.00000000590362766598E0, -0.00000006642090827576E0, &
3926 0.00000067595634268133E0, -0.00000621188515924000E0, &
3927 0.00005103883009709690E0, -0.00037015410692956173E0, &
3928 0.00233307631218880978E0, -0.01254988477182192210E0, &
3929 0.05657061146827041994E0, -0.21379664776456006580E0, &
3930 0.84270079294971486929E0, &
3931 0.00000000000949905026E0, -0.00000000018310229805E0, &
3932 0.00000000239463074000E0, -0.00000002721444369609E0, &
3933 0.00000028045522331686E0, -0.00000261830022482897E0, &
3934 0.00002195455056768781E0, -0.00016358986921372656E0, &
3935 0.00107052153564110318E0, -0.00608284718113590151E0, &
3936 0.02986978465246258244E0, -0.13055593046562267625E0, &
3937 0.67493323603965504676E0, &
3938 0.00000000000382722073E0, -0.00000000007421598602E0, &
3939 0.00000000097930574080E0, -0.00000001126008898854E0, &
3940 0.00000011775134830784E0, -0.00000111992758382650E0, &
3941 0.00000962023443095201E0, -0.00007404402135070773E0, &
3942 0.00050689993654144881E0, -0.00307553051439272889E0, &
3943 0.01668977892553165586E0, -0.08548534594781312114E0, &
3944 0.56909076642393639985E0, &
3945 0.00000000000155296588E0, -0.00000000003032205868E0, &
3946 0.00000000040424830707E0, -0.00000000471135111493E0, &
3947 0.00000005011915876293E0, -0.00000048722516178974E0, &
3948 0.00000430683284629395E0, -0.00003445026145385764E0, &
3949 0.00024879276133931664E0, -0.00162940941748079288E0, &
3950 0.00988786373932350462E0, -0.05962426839442303805E0, &
3951 0.49766113250947636708E0 /
3952 DATA (B(I), I = 0, 12) / &
3953 -0.00000000029734388465E0, 0.00000000269776334046E0, &
3954 -0.00000000640788827665E0, -0.00000001667820132100E0, &
3955 -0.00000021854388148686E0, 0.00000266246030457984E0, &
3956 0.00001612722157047886E0, -0.00025616361025506629E0, &
3957 0.00015380842432375365E0, 0.00815533022524927908E0, &
3958 -0.01402283663896319337E0, -0.19746892495383021487E0, &
3959 0.71511720328842845913E0 /
3960 DATA (B(I), I = 13, 25) / &
3961 -0.00000000001951073787E0, -0.00000000032302692214E0, &
3962 0.00000000522461866919E0, 0.00000000342940918551E0, &
3963 -0.00000035772874310272E0, 0.00000019999935792654E0, &
3964 0.00002687044575042908E0, -0.00011843240273775776E0, &
3965 -0.00080991728956032271E0, 0.00661062970502241174E0, &
3966 0.00909530922354827295E0, -0.20160072778491013140E0, &
3967 0.51169696718727644908E0 /
3968 DATA (B(I), I = 26, 38) / &
3969 0.00000000003147682272E0, -0.00000000048465972408E0, &
3970 0.00000000063675740242E0, 0.00000003377623323271E0, &
3971 -0.00000015451139637086E0, -0.00000203340624738438E0, &
3972 0.00001947204525295057E0, 0.00002854147231653228E0, &
3973 -0.00101565063152200272E0, 0.00271187003520095655E0, &
3974 0.02328095035422810727E0, -0.16725021123116877197E0, &
3975 0.32490054966649436974E0 /
3976 DATA (B(I), I = 39, 51) / &
3977 0.00000000002319363370E0, -0.00000000006303206648E0, &
3978 -0.00000000264888267434E0, 0.00000002050708040581E0, &
3979 0.00000011371857327578E0, -0.00000211211337219663E0, &
3980 0.00000368797328322935E0, 0.00009823686253424796E0, &
3981 -0.00065860243990455368E0, -0.00075285814895230877E0, &
3982 0.02585434424202960464E0, -0.11637092784486193258E0, &
3983 0.18267336775296612024E0 /
3984 DATA (B(I), I = 52, 64) / &
3985 -0.00000000000367789363E0, 0.00000000020876046746E0, &
3986 -0.00000000193319027226E0, -0.00000000435953392472E0, &
3987 0.00000018006992266137E0, -0.00000078441223763969E0, &
3988 -0.00000675407647949153E0, 0.00008428418334440096E0, &
3989 -0.00017604388937031815E0, -0.00239729611435071610E0, &
3990 0.02064129023876022970E0, -0.06905562880005864105E0, &
3991 0.09084526782065478489E0 /
3992 W = ABS(X)
3993 IF (W .LT. 2.2D0) THEN
3994 T = W * W
3995 K = INT(T)
3996 T = T - K
3997 K = K * 13
3998 Y = ((((((((((((A(K) * T + A(K + 1)) * T + &
3999 A(K + 2)) * T + A(K + 3)) * T + A(K + 4)) * T + &
4000 A(K + 5)) * T + A(K + 6)) * T + A(K + 7)) * T + &
4001 A(K + 8)) * T + A(K + 9)) * T + A(K + 10)) * T + &
4002 A(K + 11)) * T + A(K + 12)) * W
4003 ELSE IF (W .LT. 6.9D0) THEN
4004 K = INT(W)
4005 T = W - K
4006 K = 13 * (K - 2)
4007 Y = (((((((((((B(K) * T + B(K + 1)) * T + &
4008 B(K + 2)) * T + B(K + 3)) * T + B(K + 4)) * T + &
4009 B(K + 5)) * T + B(K + 6)) * T + B(K + 7)) * T + &
4010 B(K + 8)) * T + B(K + 9)) * T + B(K + 10)) * T + &
4011 B(K + 11)) * T + B(K + 12)
4012 Y = Y * Y
4013 Y = Y * Y
4014 Y = Y * Y
4015 Y = 1 - Y * Y
4016 ELSE
4017 Y = 1
4018 END IF
4019 IF (X .LT. 0) Y = -Y
4020 DERF1 = Y
4021 END FUNCTION DERF1
4022
4023 !+---+-----------------------------------------------------------------+
4024
4025 END MODULE module_mp_morr_two_moment