Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

module_mp_morr_two_moment.F @ 3094

Last change on this file since 3094 was 2759, checked in by aslmd, 2 years ago
adding unmodified code from WRFV3.0.1.1, expurged from useless data +1M size
File size: 136.8 KB

Line
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 = 250.
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 = 0
242
243	! SWITCH FOR GRAUPEL/HAIL NO GRAUPEL/HAIL
244	! IGRAUP = 0, INCLUDE GRAUPEL/HAIL
245	! IGRAUP = 1, NO GRAUPEL/HAIL
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 = 461.5
281	CP = 1005.
282	RHOSU = 85000./(287.15*273.15)
283	RHOW = 997.
284	RHOI = 500.
285	RHOSN = 100.
286	IF (IHAIL.EQ.0) THEN
287	RHOG = 400.
288	ELSE
289	RHOG = 900.
290	END IF
291	AIMM = 0.66
292	BIMM = 100.
293	ECR = 1.
294	DCS = 125.E-6
295	MI0 = 4./3.PIRHOI(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.PIRHOI(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 = VIOSMEPSMMWRHOA/(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.5EXP(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.5EXP(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.EIIPI*((1.-BS)/3.)RHOSN*((-2.-BS)/3.)/(4.720.)
393	CONS16=GAMMA(BI+3.)PI/4.ECI
394	CONS17=4.2.3.RHOSUPIECIECIGAMMA(2.BS+2.)/(8.*(RHOG-RHOSN))
395	CONS18=RHOSN*RHOSN
396	CONS19=RHOW*RHOW
397	CONS20=20.PIPIRHOWBIMM
398	CONS21=4./(DCS*RHOI)
399	CONS22=PIRHOIDCS**3/6.
400	CONS23=PI/4.EIIGAMMA(BS+3.)
401	CONS24=PI/4.ECRGAMMA(BR+3.)
402	CONS25=PIPI/24.RHOWECRGAMMA(BR+6.)
403	CONS26=PI/6.*RHOW
404	CONS27=GAMMA(1.+BI)
405	CONS28=GAMMA(4.+BI)/6.
406	CONS29=4./3.PIRHOW(25.E-6)*3
407	CONS30=4./3.PIRHOW
408	CONS31=PIPIECR*RHOSN
409	CONS32=PI/2.*ECR
410	CONS33=PIPIECR*RHOG
411	CONS34=5./2.+BR/2.
412	CONS35=5./2.+BS/2.
413	CONS36=5./2.+BG/2.
414	CONS37=4.PI1.38E-23/(6.PIRIN)
415	CONS38=PIPI/3.RHOW
416	CONS39=PIPI/36.RHOW*BIMM
417	CONS40=PI/6.*BIMM
418	CONS41=PIPIECR*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(its:ite, kts:kte, jts:jte):: &
538	effi, effs, effr, EFFG
539
540	REAL, DIMENSION(its:ite, kts:kte, jts:jte):: &
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(KTS:KTE) :: QC3DTEN ! CLOUD WATER MIXING RATIO TENDENCY (KG/KG/S)
753	REAL, DIMENSION(KTS:KTE) :: QI3DTEN ! CLOUD ICE MIXING RATIO TENDENCY (KG/KG/S)
754	REAL, DIMENSION(KTS:KTE) :: QNI3DTEN ! SNOW MIXING RATIO TENDENCY (KG/KG/S)
755	REAL, DIMENSION(KTS:KTE) :: QR3DTEN ! RAIN MIXING RATIO TENDENCY (KG/KG/S)
756	REAL, DIMENSION(KTS:KTE) :: NI3DTEN ! CLOUD ICE NUMBER CONCENTRATION (1/KG/S)
757	REAL, DIMENSION(KTS:KTE) :: NS3DTEN ! SNOW NUMBER CONCENTRATION (1/KG/S)
758	REAL, DIMENSION(KTS:KTE) :: NR3DTEN ! RAIN NUMBER CONCENTRATION (1/KG/S)
759	REAL, DIMENSION(KTS:KTE) :: QC3D ! CLOUD WATER MIXING RATIO (KG/KG)
760	REAL, DIMENSION(KTS:KTE) :: QI3D ! CLOUD ICE MIXING RATIO (KG/KG)
761	REAL, DIMENSION(KTS:KTE) :: QNI3D ! SNOW MIXING RATIO (KG/KG)
762	REAL, DIMENSION(KTS:KTE) :: QR3D ! RAIN MIXING RATIO (KG/KG)
763	REAL, DIMENSION(KTS:KTE) :: NI3D ! CLOUD ICE NUMBER CONCENTRATION (1/KG)
764	REAL, DIMENSION(KTS:KTE) :: NS3D ! SNOW NUMBER CONCENTRATION (1/KG)
765	REAL, DIMENSION(KTS:KTE) :: NR3D ! RAIN NUMBER CONCENTRATION (1/KG)
766	REAL, DIMENSION(KTS:KTE) :: T3DTEN ! TEMPERATURE TENDENCY (K/S)
767	REAL, DIMENSION(KTS:KTE) :: QV3DTEN ! WATER VAPOR MIXING RATIO TENDENCY (KG/KG/S)
768	REAL, DIMENSION(KTS:KTE) :: T3D ! TEMPERATURE (K)
769	REAL, DIMENSION(KTS:KTE) :: QV3D ! WATER VAPOR MIXING RATIO (KG/KG)
770	REAL, DIMENSION(KTS:KTE) :: PRES ! ATMOSPHERIC PRESSURE (PA)
771	REAL, DIMENSION(KTS:KTE) :: DZQ ! DIFFERENCE IN HEIGHT ACROSS LEVEL (m)
772	REAL, DIMENSION(KTS:KTE) :: W3D ! GRID-SCALE VERTICAL VELOCITY (M/S)
773	REAL, DIMENSION(KTS:KTE) :: WVAR ! SUB-GRID VERTICAL VELOCITY (M/S)
774
775	! HM ADDED GRAUPEL VARIABLES
776	REAL, DIMENSION(KTS:KTE) :: QG3DTEN ! GRAUPEL MIX RATIO TENDENCY (KG/KG/S)
777	REAL, DIMENSION(KTS:KTE) :: NG3DTEN ! GRAUPEL NUMB CONC TENDENCY (1/KG/S)
778	REAL, DIMENSION(KTS:KTE) :: QG3D ! GRAUPEL MIX RATIO (KG/KG)
779	REAL, DIMENSION(KTS:KTE) :: NG3D ! GRAUPEL NUMBER CONC (1/KG)
780
781	! HM, ADD 1/16/07, SEDIMENTATION TENDENCIES FOR MIXING RATIO
782
783	REAL, DIMENSION(KTS:KTE) :: QGSTEN ! GRAUPEL SED TEND (KG/KG/S)
784	REAL, DIMENSION(KTS:KTE) :: QRSTEN ! RAIN SED TEND (KG/KG/S)
785	REAL, DIMENSION(KTS:KTE) :: QISTEN ! CLOUD ICE SED TEND (KG/KG/S)
786	REAL, DIMENSION(KTS:KTE) :: QNISTEN ! SNOW SED TEND (KG/KG/S)
787	REAL, DIMENSION(KTS:KTE) :: QCSTEN ! CLOUD WAT SED TEND (KG/KG/S)
788
789	! hm add cumulus tendencies for precip
790	REAL, DIMENSION(KTS:KTE) :: qrcu1d
791	REAL, DIMENSION(KTS:KTE) :: qscu1d
792	REAL, DIMENSION(KTS:KTE) :: qicu1d
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(KTS:KTE) :: EFFC ! DROPLET EFFECTIVE RADIUS (MICRON)
800	REAL, DIMENSION(KTS:KTE) :: EFFI ! CLOUD ICE EFFECTIVE RADIUS (MICRON)
801	REAL, DIMENSION(KTS:KTE) :: EFFS ! SNOW EFFECTIVE RADIUS (MICRON)
802	REAL, DIMENSION(KTS:KTE) :: EFFR ! RAIN EFFECTIVE RADIUS (MICRON)
803	REAL, DIMENSION(KTS:KTE) :: EFFG ! GRAUPEL EFFECTIVE RADIUS (MICRON)
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(KTS:KTE) :: LAMC ! SLOPE PARAMETER FOR DROPLETS (M-1)
816	REAL, DIMENSION(KTS:KTE) :: LAMI ! SLOPE PARAMETER FOR CLOUD ICE (M-1)
817	REAL, DIMENSION(KTS:KTE) :: LAMS ! SLOPE PARAMETER FOR SNOW (M-1)
818	REAL, DIMENSION(KTS:KTE) :: LAMR ! SLOPE PARAMETER FOR RAIN (M-1)
819	REAL, DIMENSION(KTS:KTE) :: LAMG ! SLOPE PARAMETER FOR GRAUPEL (M-1)
820	REAL, DIMENSION(KTS:KTE) :: CDIST1 ! PSD PARAMETER FOR DROPLETS
821	REAL, DIMENSION(KTS:KTE) :: N0I ! INTERCEPT PARAMETER FOR CLOUD ICE (KG-1 M-1)
822	REAL, DIMENSION(KTS:KTE) :: N0S ! INTERCEPT PARAMETER FOR SNOW (KG-1 M-1)
823	REAL, DIMENSION(KTS:KTE) :: N0RR ! INTERCEPT PARAMETER FOR RAIN (KG-1 M-1)
824	REAL, DIMENSION(KTS:KTE) :: N0G ! INTERCEPT PARAMETER FOR GRAUPEL (KG-1 M-1)
825	REAL, DIMENSION(KTS:KTE) :: PGAM ! SPECTRAL SHAPE PARAMETER FOR DROPLETS
826
827	! MICROPHYSICAL PROCESSES
828
829	REAL, DIMENSION(KTS:KTE) :: NSUBC ! LOSS OF NC DURING EVAP
830	REAL, DIMENSION(KTS:KTE) :: NSUBI ! LOSS OF NI DURING SUB.
831	REAL, DIMENSION(KTS:KTE) :: NSUBS ! LOSS OF NS DURING SUB.
832	REAL, DIMENSION(KTS:KTE) :: NSUBR ! LOSS OF NR DURING EVAP
833	REAL, DIMENSION(KTS:KTE) :: PRD ! DEP CLOUD ICE
834	REAL, DIMENSION(KTS:KTE) :: PRE ! EVAP OF RAIN
835	REAL, DIMENSION(KTS:KTE) :: PRDS ! DEP SNOW
836	REAL, DIMENSION(KTS:KTE) :: NNUCCC ! CHANGE N DUE TO CONTACT FREEZ DROPLETS
837	REAL, DIMENSION(KTS:KTE) :: MNUCCC ! CHANGE Q DUE TO CONTACT FREEZ DROPLETS
838	REAL, DIMENSION(KTS:KTE) :: PRA ! ACCRETION DROPLETS BY RAIN
839	REAL, DIMENSION(KTS:KTE) :: PRC ! AUTOCONVERSION DROPLETS
840	REAL, DIMENSION(KTS:KTE) :: PCC ! COND/EVAP DROPLETS
841	REAL, DIMENSION(KTS:KTE) :: NNUCCD ! CHANGE N FREEZING AEROSOL (PRIM ICE NUCLEATION)
842	REAL, DIMENSION(KTS:KTE) :: MNUCCD ! CHANGE Q FREEZING AEROSOL (PRIM ICE NUCLEATION)
843	REAL, DIMENSION(KTS:KTE) :: MNUCCR ! CHANGE Q DUE TO CONTACT FREEZ RAIN
844	REAL, DIMENSION(KTS:KTE) :: NNUCCR ! CHANGE N DUE TO CONTACT FREEZ RAIN
845	REAL, DIMENSION(KTS:KTE) :: NPRA ! CHANGE IN N DUE TO DROPLET ACC BY RAIN
846	REAL, DIMENSION(KTS:KTE) :: NRAGG ! SELF-COLLECTION OF RAIN
847	REAL, DIMENSION(KTS:KTE) :: NSAGG ! SELF-COLLECTION OF SNOW
848	REAL, DIMENSION(KTS:KTE) :: NPRC ! CHANGE NC AUTOCONVERSION DROPLETS
849	REAL, DIMENSION(KTS:KTE) :: NPRC1 ! CHANGE NR AUTOCONVERSION DROPLETS
850	REAL, DIMENSION(KTS:KTE) :: PRAI ! CHANGE Q AUTOCONVERSION CLOUD ICE
851	REAL, DIMENSION(KTS:KTE) :: PRCI ! CHANGE Q ACCRETION CLOUD ICE BY SNOW
852	REAL, DIMENSION(KTS:KTE) :: PSACWS ! CHANGE Q DROPLET ACCRETION BY SNOW
853	REAL, DIMENSION(KTS:KTE) :: NPSACWS ! CHANGE N DROPLET ACCRETION BY SNOW
854	REAL, DIMENSION(KTS:KTE) :: PSACWI ! CHANGE Q DROPLET ACCRETION BY CLOUD ICE
855	REAL, DIMENSION(KTS:KTE) :: NPSACWI ! CHANGE N DROPLET ACCRETION BY CLOUD ICE
856	REAL, DIMENSION(KTS:KTE) :: NPRCI ! CHANGE N AUTOCONVERSION CLOUD ICE BY SNOW
857	REAL, DIMENSION(KTS:KTE) :: NPRAI ! CHANGE N ACCRETION CLOUD ICE
858	REAL, DIMENSION(KTS:KTE) :: NMULTS ! ICE MULT DUE TO RIMING DROPLETS BY SNOW
859	REAL, DIMENSION(KTS:KTE) :: NMULTR ! ICE MULT DUE TO RIMING RAIN BY SNOW
860	REAL, DIMENSION(KTS:KTE) :: QMULTS ! CHANGE Q DUE TO ICE MULT DROPLETS/SNOW
861	REAL, DIMENSION(KTS:KTE) :: QMULTR ! CHANGE Q DUE TO ICE RAIN/SNOW
862	REAL, DIMENSION(KTS:KTE) :: PRACS ! CHANGE Q RAIN-SNOW COLLECTION
863	REAL, DIMENSION(KTS:KTE) :: NPRACS ! CHANGE N RAIN-SNOW COLLECTION
864	REAL, DIMENSION(KTS:KTE) :: PCCN ! CHANGE Q DROPLET ACTIVATION
865	REAL, DIMENSION(KTS:KTE) :: PSMLT ! CHANGE Q MELTING SNOW TO RAIN
866	REAL, DIMENSION(KTS:KTE) :: EVPMS ! CHNAGE Q MELTING SNOW EVAPORATING
867	REAL, DIMENSION(KTS:KTE) :: NSMLTS ! CHANGE N MELTING SNOW
868	REAL, DIMENSION(KTS:KTE) :: NSMLTR ! CHANGE N MELTING SNOW TO RAIN
869	! HM ADDED 12/13/06
870	REAL, DIMENSION(KTS:KTE) :: PIACR ! CHANGE QR, ICE-RAIN COLLECTION
871	REAL, DIMENSION(KTS:KTE) :: NIACR ! CHANGE N, ICE-RAIN COLLECTION
872	REAL, DIMENSION(KTS:KTE) :: PRACI ! CHANGE QI, ICE-RAIN COLLECTION
873	REAL, DIMENSION(KTS:KTE) :: PIACRS ! CHANGE QR, ICE RAIN COLLISION, ADDED TO SNOW
874	REAL, DIMENSION(KTS:KTE) :: NIACRS ! CHANGE N, ICE RAIN COLLISION, ADDED TO SNOW
875	REAL, DIMENSION(KTS:KTE) :: PRACIS ! CHANGE QI, ICE RAIN COLLISION, ADDED TO SNOW
876	REAL, DIMENSION(KTS:KTE) :: EPRD ! SUBLIMATION CLOUD ICE
877	REAL, DIMENSION(KTS:KTE) :: EPRDS ! SUBLIMATION SNOW
878	! HM ADDED GRAUPEL PROCESSES
879	REAL, DIMENSION(KTS:KTE) :: PRACG ! CHANGE IN Q COLLECTION RAIN BY GRAUPEL
880	REAL, DIMENSION(KTS:KTE) :: PSACWG ! CHANGE IN Q COLLECTION DROPLETS BY GRAUPEL
881	REAL, DIMENSION(KTS:KTE) :: PGSACW ! CONVERSION Q TO GRAUPEL DUE TO COLLECTION DROPLETS BY SNOW
882	REAL, DIMENSION(KTS:KTE) :: PGRACS ! CONVERSION Q TO GRAUPEL DUE TO COLLECTION RAIN BY SNOW
883	REAL, DIMENSION(KTS:KTE) :: PRDG ! DEP OF GRAUPEL
884	REAL, DIMENSION(KTS:KTE) :: EPRDG ! SUB OF GRAUPEL
885	REAL, DIMENSION(KTS:KTE) :: EVPMG ! CHANGE Q MELTING OF GRAUPEL AND EVAPORATION
886	REAL, DIMENSION(KTS:KTE) :: PGMLT ! CHANGE Q MELTING OF GRAUPEL
887	REAL, DIMENSION(KTS:KTE) :: NPRACG ! CHANGE N COLLECTION RAIN BY GRAUPEL
888	REAL, DIMENSION(KTS:KTE) :: NPSACWG ! CHANGE N COLLECTION DROPLETS BY GRAUPEL
889	REAL, DIMENSION(KTS:KTE) :: NSCNG ! CHANGE N CONVERSION TO GRAUPEL DUE TO COLLECTION DROPLETS BY SNOW
890	REAL, DIMENSION(KTS:KTE) :: NGRACS ! CHANGE N CONVERSION TO GRAUPEL DUE TO COLLECTION RAIN BY SNOW
891	REAL, DIMENSION(KTS:KTE) :: NGMLTG ! CHANGE N MELTING GRAUPEL
892	REAL, DIMENSION(KTS:KTE) :: NGMLTR ! CHANGE N MELTING GRAUPEL TO RAIN
893	REAL, DIMENSION(KTS:KTE) :: NSUBG ! CHANGE N SUB/DEP OF GRAUPEL
894	REAL, DIMENSION(KTS:KTE) :: PSACR ! CONVERSION DUE TO COLL OF SNOW BY RAIN
895	REAL, DIMENSION(KTS:KTE) :: NMULTG ! ICE MULT DUE TO ACC DROPLETS BY GRAUPEL
896	REAL, DIMENSION(KTS:KTE) :: NMULTRG ! ICE MULT DUE TO ACC RAIN BY GRAUPEL
897	REAL, DIMENSION(KTS:KTE) :: QMULTG ! CHANGE Q DUE TO ICE MULT DROPLETS/GRAUPEL
898	REAL, DIMENSION(KTS:KTE) :: QMULTRG ! CHANGE Q DUE TO ICE MULT RAIN/GRAUPEL
899
900	! TIME-VARYING ATMOSPHERIC PARAMETERS
901
902	REAL, DIMENSION(KTS:KTE) :: KAP ! THERMAL CONDUCTIVITY OF AIR
903	REAL, DIMENSION(KTS:KTE) :: EVS ! SATURATION VAPOR PRESSURE
904	REAL, DIMENSION(KTS:KTE) :: EIS ! ICE SATURATION VAPOR PRESSURE
905	REAL, DIMENSION(KTS:KTE) :: QVS ! SATURATION MIXING RATIO
906	REAL, DIMENSION(KTS:KTE) :: QVI ! ICE SATURATION MIXING RATIO
907	REAL, DIMENSION(KTS:KTE) :: QVQVS ! SAUTRATION RATIO
908	REAL, DIMENSION(KTS:KTE) :: QVQVSI! ICE SATURAION RATIO
909	REAL, DIMENSION(KTS:KTE) :: DV ! DIFFUSIVITY OF WATER VAPOR IN AIR
910	REAL, DIMENSION(KTS:KTE) :: XXLS ! LATENT HEAT OF SUBLIMATION
911	REAL, DIMENSION(KTS:KTE) :: XXLV ! LATENT HEAT OF VAPORIZATION
912	REAL, DIMENSION(KTS:KTE) :: CPM ! SPECIFIC HEAT AT CONST PRESSURE FOR MOIST AIR
913	REAL, DIMENSION(KTS:KTE) :: MU ! VISCOCITY OF AIR
914	REAL, DIMENSION(KTS:KTE) :: SC ! SCHMIDT NUMBER
915	REAL, DIMENSION(KTS:KTE) :: XLF ! LATENT HEAT OF FREEZING
916	REAL, DIMENSION(KTS:KTE) :: RHO ! AIR DENSITY
917	REAL, DIMENSION(KTS:KTE) :: AB ! CORRECTION TO CONDENSATION RATE DUE TO LATENT HEATING
918	REAL, DIMENSION(KTS:KTE) :: ABI ! CORRECTION TO DEPOSITION RATE DUE TO LATENT HEATING
919
920	! TIME-VARYING MICROPHYSICS PARAMETERS
921
922	REAL, DIMENSION(KTS:KTE) :: DAP ! DIFFUSIVITY OF AEROSOL
923	REAL NACNT ! NUMBER OF CONTACT IN
924	REAL FMULT ! TEMP.-DEP. PARAMETER FOR RIME-SPLINTERING
925	REAL COFFI ! ICE AUTOCONVERSION PARAMETER
926
927	! FALL SPEED WORKING VARIABLES (DEFINED IN CODE)
928
929	REAL, DIMENSION(KTS:KTE) :: DUMI,DUMR,DUMFNI,DUMG,DUMFNG
930	REAL UNI, UMI,UMR
931	REAL, DIMENSION(KTS:KTE) :: FR, FI, FNI,FG,FNG
932	REAL RGVM
933	REAL, DIMENSION(KTS:KTE) :: FALOUTR,FALOUTI,FALOUTNI
934	REAL FALTNDR,FALTNDI,FALTNDNI,RHO2
935	REAL, DIMENSION(KTS:KTE) :: DUMQS,DUMFNS
936	REAL UMS,UNS
937	REAL, DIMENSION(KTS:KTE) :: FS,FNS, FALOUTS,FALOUTNS,FALOUTG,FALOUTNG
938	REAL FALTNDS,FALTNDNS,UNR,FALTNDG,FALTNDNG
939	REAL, DIMENSION(KTS:KTE) :: DUMC,DUMFNC
940	REAL UNC,UMC,UNG,UMG
941	REAL, DIMENSION(KTS:KTE) :: FC,FALOUTC,FALOUTNC
942	REAL FALTNDC,FALTNDNC
943	REAL, DIMENSION(KTS:KTE) :: FNC,DUMFNR,FALOUTNR
944	REAL FALTNDNR
945	REAL, DIMENSION(KTS:KTE) :: FNR(KMS:KME)
946
947	! FALL-SPEED PARAMETER 'A' WITH AIR DENSITY CORRECTION
948
949	REAL, DIMENSION(KTS:KTE) :: 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(KTS:KTE) :: NC3DTEN ! CLOUD DROPLET NUMBER CONCENTRATION (1/KG/S)
1026	REAL, DIMENSION(KTS:KTE) :: 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.887QV3D(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/(PIRHOWRHO(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/(CONS1RHO(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-6T3D(K)*1.5/(T3D(K)+120.)
1185
1186	! KAP(K) = 1.414E31.496E-6T3D(K)**1.5/(T3D(K)+120.)
1187	KAP(K) = 1.414E3*DUM
1188
1189	! DIFFUSIVITY OF WATER VAPOR
1190
1191	DV(K) = 8.794E-5T3D(K)*1.81/PRES(K)
1192
1193	! VISCOSITY OF AIR
1194	! SCHMIT NUMBER
1195
1196	! MU(K) = 1.496E-6T3D(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 = (RVT3D(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) = (PIRHOWNR3D(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)*4QR3D(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)*4QR3D(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) = (CONS26NC3D(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) = (CONS1NS3D(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)*4QNI3D(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)*4QNI3D(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) = (CONS2NG3D(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)*4QG3D(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)*4QG3D(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.e6RHO(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.2UMR-0.95UMS)*2+ &
1448	0.08UMSUMR)*0.5RHO(K)* &
1449	N0RR(K)N0S(K)/LAMS(K)3 &
1450	(5./(LAMS(K)*3LAMR(K))+ &
1451	2./(LAMS(K)*2LAMR(K)**2)+ &
1452	0.5/(LAMS(K)LAMR(K)*3)))
1453
1454	NPRACS(K) = CONS32RHO(K)(1.7(UNR-UNS)*2+ &
1455	0.3UNRUNS)*0.5N0RR(K)N0S(K) &
1456	(1./(LAMR(K)*3LAMS(K))+ &
1457	1./(LAMR(K)*2LAMS(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.2UMR-0.95UMG)*2+ &
1481	0.08UMGUMR)*0.5RHO(K)* &
1482	N0RR(K)N0G(K)/LAMR(K)3 &
1483	(5./(LAMR(K)*3LAMG(K))+ &
1484	2./(LAMR(K)*2LAMG(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) = CONS32RHO(K)(1.7(UNR-UNG)*2+ &
1492	0.3UNRUNG)*0.5N0RR(K)N0G(K) &
1493	(1./(LAMR(K)*3LAMG(K))+ &
1494	1./(LAMR(K)*2LAMG(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.PIN0RR(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.PIN0S(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.PIN0S(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.PIN0S(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.PIN0G(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.PIN0G(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.PIN0G(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	! HM, bug fix 5/12/08, npracg is subtracted from nr not ng
1679	! NG3DTEN(K) = NG3DTEN(K)
1680	NC3DTEN(K) = NC3DTEN(K)+ (-NPRA(K)-NPRC(K))
1681	NR3DTEN(K) = NR3DTEN(K)+ (NPRC1(K)+NRAGG(K)-NPRACG(K))
1682
1683	IF (PRE(K).LT.0.) THEN
1684	DUM = PRE(K)*DT/QR3D(K)
1685	DUM = MAX(-1.,DUM)
1686	NSUBR(K) = DUM*NR3D(K)/DT
1687	END IF
1688
1689	IF (EVPMS(K)+PSMLT(K).LT.0.) THEN
1690	DUM = (EVPMS(K)+PSMLT(K))*DT/QNI3D(K)
1691	DUM = MAX(-1.,DUM)
1692	NSMLTS(K) = DUM*NS3D(K)/DT
1693	END IF
1694	IF (PSMLT(K).LT.0.) THEN
1695	DUM = PSMLT(K)*DT/QNI3D(K)
1696	DUM = MAX(-1.0,DUM)
1697	NSMLTR(K) = DUM*NS3D(K)/DT
1698	END IF
1699	IF (EVPMG(K)+PGMLT(K).LT.0.) THEN
1700	DUM = (EVPMG(K)+PGMLT(K))*DT/QG3D(K)
1701	DUM = MAX(-1.,DUM)
1702	NGMLTG(K) = DUM*NG3D(K)/DT
1703	END IF
1704	IF (PGMLT(K).LT.0.) THEN
1705	DUM = PGMLT(K)*DT/QG3D(K)
1706	DUM = MAX(-1.0,DUM)
1707	NGMLTR(K) = DUM*NG3D(K)/DT
1708	END IF
1709
1710	NS3DTEN(K) = NS3DTEN(K)+(NSMLTS(K))
1711	NG3DTEN(K) = NG3DTEN(K)+(NGMLTG(K))
1712	NR3DTEN(K) = NR3DTEN(K)+(NSUBR(K)-NSMLTR(K)-NGMLTR(K))
1713
1714	300 CONTINUE
1715
1716	!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1717	! NOW CALCULATE SATURATION ADJUSTMENT TO CONDENSE EXTRA VAPOR ABOVE
1718	! WATER SATURATION
1719
1720	DUMT = T3D(K)+DT*T3DTEN(K)
1721	DUMQV = QV3D(K)+DT*QV3DTEN(K)
1722	DUMQSS = 0.622*POLYSVP(DUMT,0)/ (PRES(K)-POLYSVP(DUMT,0))
1723	DUMQC = QC3D(K)+DT*QC3DTEN(K)
1724	DUMQC = MAX(DUMQC,0.)
1725
1726	! SATURATION ADJUSTMENT FOR LIQUID
1727
1728	DUMS = DUMQV-DUMQSS
1729	PCC(K) = DUMS/(1.+XXLV(K)*2DUMQSS/(CPM(K)RVDUMT**2))/DT
1730	IF (PCC(K)*DT+DUMQC.LT.0.) THEN
1731	PCC(K) = -DUMQC/DT
1732	END IF
1733
1734	QV3DTEN(K) = QV3DTEN(K)-PCC(K)
1735	T3DTEN(K) = T3DTEN(K)+PCC(K)*XXLV(K)/CPM(K)
1736	QC3DTEN(K) = QC3DTEN(K)+PCC(K)
1737
1738	!.......................................................................
1739	! ACTIVATION OF CLOUD DROPLETS
1740	! ACTIVATION OF DROPLET CURRENTLY NOT CALCULATED
1741	! DROPLET CONCENTRATION IS SPECIFIED !!!!!
1742
1743	!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1744	! SUBLIMATE, MELT, OR EVAPORATE NUMBER CONCENTRATION
1745	! THIS FORMULATION ASSUMES 1:1 RATIO BETWEEN MASS LOSS AND
1746	! LOSS OF NUMBER CONCENTRATION
1747
1748	! IF (PCC(K).LT.0.) THEN
1749	! DUM = PCC(K)*DT/QC3D(K)
1750	! DUM = MAX(-1.,DUM)
1751	! NSUBC(K) = DUM*NC3D(K)/DT
1752	! END IF
1753
1754	! UPDATE TENDENCIES
1755
1756	! NC3DTEN(K) = NC3DTEN(K)+NSUBC(K)
1757
1758	!.....................................................................
1759	!.....................................................................
1760	ELSE ! TEMPERATURE < 273.15
1761
1762	!......................................................................
1763	!HM ADD, ALLOW FOR CONSTANT DROPLET NUMBER
1764	! INUM = 0, PREDICT DROPLET NUMBER
1765	! INUM = 1, SET CONSTANT DROPLET NUMBER
1766
1767	! IF (INUM.EQ.1) THEN
1768	! CONVERT NDCNST FROM CM-3 TO KG-1
1769	NC3D(K)=NDCNST*1.E6/RHO(K)
1770	! END IF
1771
1772	! CALCULATE SIZE DISTRIBUTION PARAMETERS
1773	! MAKE SURE NUMBER CONCENTRATIONS AREN'T NEGATIVE
1774
1775	NI3D(K) = MAX(0.,NI3D(K))
1776	NS3D(K) = MAX(0.,NS3D(K))
1777	NC3D(K) = MAX(0.,NC3D(K))
1778	NR3D(K) = MAX(0.,NR3D(K))
1779	NG3D(K) = MAX(0.,NG3D(K))
1780
1781	!......................................................................
1782	! CLOUD ICE
1783
1784	IF (QI3D(K).GE.QSMALL) THEN
1785	LAMI(K) = (CONS12* &
1786	NI3D(K)/QI3D(K))**(1./DI)
1787	N0I(K) = NI3D(K)*LAMI(K)
1788
1789	! CHECK FOR SLOPE
1790
1791	! ADJUST VARS
1792
1793	IF (LAMI(K).LT.LAMMINI) THEN
1794
1795	LAMI(K) = LAMMINI
1796
1797	N0I(K) = LAMI(K)*4QI3D(K)/CONS12
1798
1799	NI3D(K) = N0I(K)/LAMI(K)
1800	ELSE IF (LAMI(K).GT.LAMMAXI) THEN
1801	LAMI(K) = LAMMAXI
1802	N0I(K) = LAMI(K)*4QI3D(K)/CONS12
1803
1804	NI3D(K) = N0I(K)/LAMI(K)
1805	END IF
1806	END IF
1807
1808	!......................................................................
1809	! RAIN
1810
1811	IF (QR3D(K).GE.QSMALL) THEN
1812	LAMR(K) = (PIRHOWNR3D(K)/QR3D(K))**(1./3.)
1813	N0RR(K) = NR3D(K)*LAMR(K)
1814
1815	! CHECK FOR SLOPE
1816
1817	! ADJUST VARS
1818
1819	IF (LAMR(K).LT.LAMMINR) THEN
1820
1821	LAMR(K) = LAMMINR
1822
1823	N0RR(K) = LAMR(K)*4QR3D(K)/(PI*RHOW)
1824
1825	NR3D(K) = N0RR(K)/LAMR(K)
1826	ELSE IF (LAMR(K).GT.LAMMAXR) THEN
1827	LAMR(K) = LAMMAXR
1828	N0RR(K) = LAMR(K)*4QR3D(K)/(PI*RHOW)
1829
1830	NR3D(K) = N0RR(K)/LAMR(K)
1831	END IF
1832	END IF
1833
1834	!......................................................................
1835	! CLOUD DROPLETS
1836
1837	! MARTIN ET AL. (1994) FORMULA FOR PGAM
1838
1839	IF (QC3D(K).GE.QSMALL) THEN
1840
1841	DUM = PRES(K)/(287.15*T3D(K))
1842	PGAM(K)=0.0005714*(NC3D(K)/1.E6/DUM)+0.2714
1843	PGAM(K)=1./(PGAM(K)**2)-1.
1844	PGAM(K)=MAX(PGAM(K),2.)
1845	PGAM(K)=MIN(PGAM(K),10.)
1846
1847	! CALCULATE LAMC
1848
1849	LAMC(K) = (CONS26NC3D(K)GAMMA(PGAM(K)+4.)/ &
1850	(QC3D(K)GAMMA(PGAM(K)+1.)))*(1./3.)
1851
1852	! LAMMIN, 60 MICRON DIAMETER
1853	! LAMMAX, 1 MICRON
1854
1855	LAMMIN = (PGAM(K)+1.)/60.E-6
1856	LAMMAX = (PGAM(K)+1.)/1.E-6
1857
1858	IF (LAMC(K).LT.LAMMIN) THEN
1859	LAMC(K) = LAMMIN
1860
1861	NC3D(K) = EXP(3.*LOG(LAMC(K))+LOG(QC3D(K))+ &
1862	LOG(GAMMA(PGAM(K)+1.))-LOG(GAMMA(PGAM(K)+4.)))/CONS26
1863	ELSE IF (LAMC(K).GT.LAMMAX) THEN
1864	LAMC(K) = LAMMAX
1865	NC3D(K) = EXP(3.*LOG(LAMC(K))+LOG(QC3D(K))+ &
1866	LOG(GAMMA(PGAM(K)+1.))-LOG(GAMMA(PGAM(K)+4.)))/CONS26
1867
1868	END IF
1869
1870	! TO CALCULATE DROPLET FREEZING
1871
1872	CDIST1(K) = NC3D(K)/GAMMA(PGAM(K)+1.)
1873
1874	END IF
1875
1876	!......................................................................
1877	! SNOW
1878
1879	IF (QNI3D(K).GE.QSMALL) THEN
1880	LAMS(K) = (CONS1NS3D(K)/QNI3D(K))*(1./DS)
1881	N0S(K) = NS3D(K)*LAMS(K)
1882
1883	! CHECK FOR SLOPE
1884
1885	! ADJUST VARS
1886
1887	IF (LAMS(K).LT.LAMMINS) THEN
1888	LAMS(K) = LAMMINS
1889	N0S(K) = LAMS(K)*4QNI3D(K)/CONS1
1890
1891	NS3D(K) = N0S(K)/LAMS(K)
1892
1893	ELSE IF (LAMS(K).GT.LAMMAXS) THEN
1894
1895	LAMS(K) = LAMMAXS
1896	N0S(K) = LAMS(K)*4QNI3D(K)/CONS1
1897
1898	NS3D(K) = N0S(K)/LAMS(K)
1899	END IF
1900	END IF
1901
1902	!......................................................................
1903	! GRAUPEL
1904
1905	IF (QG3D(K).GE.QSMALL) THEN
1906	LAMG(K) = (CONS2NG3D(K)/QG3D(K))*(1./DG)
1907	N0G(K) = NG3D(K)*LAMG(K)
1908
1909	! CHECK FOR SLOPE
1910
1911	! ADJUST VARS
1912
1913	IF (LAMG(K).LT.LAMMING) THEN
1914	LAMG(K) = LAMMING
1915	N0G(K) = LAMG(K)*4QG3D(K)/CONS2
1916
1917	NG3D(K) = N0G(K)/LAMG(K)
1918
1919	ELSE IF (LAMG(K).GT.LAMMAXG) THEN
1920
1921	LAMG(K) = LAMMAXG
1922	N0G(K) = LAMG(K)*4QG3D(K)/CONS2
1923
1924	NG3D(K) = N0G(K)/LAMG(K)
1925	END IF
1926	END IF
1927
1928	!.....................................................................
1929	! ZERO OUT PROCESS RATES
1930
1931	MNUCCC(K) = 0.
1932	NNUCCC(K) = 0.
1933	PRC(K) = 0.
1934	NPRC(K) = 0.
1935	NPRC1(K) = 0.
1936	NSAGG(K) = 0.
1937	PSACWS(K) = 0.
1938	NPSACWS(K) = 0.
1939	PSACWI(K) = 0.
1940	NPSACWI(K) = 0.
1941	PRACS(K) = 0.
1942	NPRACS(K) = 0.
1943	NMULTS(K) = 0.
1944	QMULTS(K) = 0.
1945	NMULTR(K) = 0.
1946	QMULTR(K) = 0.
1947	NMULTG(K) = 0.
1948	QMULTG(K) = 0.
1949	NMULTRG(K) = 0.
1950	QMULTRG(K) = 0.
1951	MNUCCR(K) = 0.
1952	NNUCCR(K) = 0.
1953	PRA(K) = 0.
1954	NPRA(K) = 0.
1955	NRAGG(K) = 0.
1956	PRCI(K) = 0.
1957	NPRCI(K) = 0.
1958	PRAI(K) = 0.
1959	NPRAI(K) = 0.
1960	NNUCCD(K) = 0.
1961	MNUCCD(K) = 0.
1962	PCC(K) = 0.
1963	PRE(K) = 0.
1964	PRD(K) = 0.
1965	PRDS(K) = 0.
1966	EPRD(K) = 0.
1967	EPRDS(K) = 0.
1968	NSUBC(K) = 0.
1969	NSUBI(K) = 0.
1970	NSUBS(K) = 0.
1971	NSUBR(K) = 0.
1972	PIACR(K) = 0.
1973	NIACR(K) = 0.
1974	PRACI(K) = 0.
1975	PIACRS(K) = 0.
1976	NIACRS(K) = 0.
1977	PRACIS(K) = 0.
1978	! HM: ADD GRAUPEL PROCESSES
1979	PRACG(K) = 0.
1980	PSACR(K) = 0.
1981	PSACWG(K) = 0.
1982	PGSACW(K) = 0.
1983	PGRACS(K) = 0.
1984	PRDG(K) = 0.
1985	EPRDG(K) = 0.
1986	NPRACG(K) = 0.
1987	NPSACWG(K) = 0.
1988	NSCNG(K) = 0.
1989	NGRACS(K) = 0.
1990	NSUBG(K) = 0.
1991
1992	!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1993	! CALCULATION OF MICROPHYSICAL PROCESS RATES
1994	! ACCRETION/AUTOCONVERSION/FREEZING/MELTING/COAG.
1995	!.......................................................................
1996	! FREEZING OF CLOUD DROPLETS
1997	! ONLY ALLOWED BELOW -4 C
1998	IF (QC3D(K).GE.QSMALL .AND. T3D(K).LT.269.15) THEN
1999
2000	! NUMBER OF CONTACT NUCLEI (M^-3) FROM MEYERS ET AL., 1992
2001	! FACTOR OF 1000 IS TO CONVERT FROM L^-1 TO M^-3
2002
2003	! MEYERS CURVE
2004
2005	NACNT = EXP(-2.80+0.262(273.15-T3D(K)))1000.
2006
2007	! COOPER CURVE
2008	! NACNT = 5.EXP(0.304(273.15-T3D(K)))
2009
2010	! FLECTHER
2011	! NACNT = 0.01EXP(0.6(273.15-T3D(K)))
2012
2013	! CONTACT FREEZING
2014
2015	! MEAN FREE PATH
2016
2017	DUM = 7.37T3D(K)/(288.10.*PRES(K))/100.
2018
2019	! EFFECTIVE DIFFUSIVITY OF CONTACT NUCLEI
2020	! BASED ON BROWNIAN DIFFUSION
2021
2022	DAP(K) = CONS37T3D(K)(1.+DUM/RIN)/MU(K)
2023
2024	MNUCCC(K) = CONS38DAP(K)NACNT*EXP(LOG(CDIST1(K))+ &
2025	LOG(GAMMA(PGAM(K)+5.))-4.*LOG(LAMC(K)))
2026	NNUCCC(K) = 2.PIDAP(K)NACNTCDIST1(K)* &
2027	GAMMA(PGAM(K)+2.)/ &
2028	LAMC(K)
2029
2030	! IMMERSION FREEZING (BIGG 1953)
2031
2032	MNUCCC(K) = MNUCCC(K)+CONS39* &
2033	EXP(LOG(CDIST1(K))+LOG(GAMMA(7.+PGAM(K)))-6.LOG(LAMC(K))) &
2034	EXP(AIMM*(273.15-T3D(K)))
2035
2036	NNUCCC(K) = NNUCCC(K)+ &
2037	CONS40EXP(LOG(CDIST1(K))+LOG(GAMMA(PGAM(K)+4.))-3.LOG(LAMC(K))) &
2038	EXP(AIMM(273.15-T3D(K)))
2039
2040	! PUT IN A CATCH HERE TO PREVENT DIVERGENCE BETWEEN NUMBER CONC. AND
2041	! MIXING RATIO, SINCE STRICT CONSERVATION NOT CHECKED FOR NUMBER CONC
2042
2043	NNUCCC(K) = MIN(NNUCCC(K),NC3D(K)/DT)
2044
2045	END IF
2046
2047	!.................................................................
2048	!.......................................................................
2049	! AUTOCONVERSION OF CLOUD LIQUID WATER TO RAIN
2050	! FORMULA FROM BEHENG (1994)
2051	! USING NUMERICAL SIMULATION OF STOCHASTIC COLLECTION EQUATION
2052	! AND INITIAL CLOUD DROPLET SIZE DISTRIBUTION SPECIFIED
2053	! AS A GAMMA DISTRIBUTION
2054
2055	! USE MINIMUM VALUE OF 1.E-6 TO PREVENT FLOATING POINT ERROR
2056
2057	IF (QC3D(K).GE.1.E-6) THEN
2058
2059	! HM ADD 12/13/06, REPLACE WITH NEWER FORMULA
2060	! FROM KHAIROUTDINOV AND KOGAN 2000, MWR
2061
2062	PRC(K)=1350.QC3D(K)2.47 &
2063	(NC3D(K)/1.e6RHO(K))*(-1.79)
2064
2065	! note: nprc1 is change in Nr,
2066	! nprc is change in Nc
2067
2068	NPRC1(K) = PRC(K)/CONS29
2069	NPRC(K) = PRC(K)/(QC3D(K)/NC3D(K))
2070
2071	NPRC(K) = MIN(NPRC(K),NC3D(K)/DT)
2072
2073	END IF
2074
2075	!.......................................................................
2076	! SELF-COLLECTION OF DROPLET NOT INCLUDED IN KK2000 SCHEME
2077
2078	! SNOW AGGREGATION FROM PASSARELLI, 1978, USED BY REISNER, 1998
2079	! THIS IS HARD-WIRED FOR BS = 0.4 FOR NOW
2080
2081	IF (QNI3D(K).GE.1.E-8) THEN
2082	NSAGG(K) = CONS15ASN(K)RHO(K)** &
2083	((2.+BS)/3.)QNI3D(K)((2.+BS)/3.) &
2084	(NS3D(K)RHO(K))*((4.-BS)/3.)/ &
2085	(RHO(K))
2086	END IF
2087
2088	!.......................................................................
2089	! ACCRETION OF CLOUD DROPLETS ONTO SNOW/GRAUPEL
2090	! HERE USE CONTINUOUS COLLECTION EQUATION WITH
2091	! SIMPLE GRAVITATIONAL COLLECTION KERNEL IGNORING
2092
2093	! SNOW
2094
2095	IF (QNI3D(K).GE.1.E-8 .AND. QC3D(K).GE.QSMALL) THEN
2096
2097	PSACWS(K) = CONS13ASN(K)QC3D(K)RHO(K) &
2098	N0S(K)/ &
2099	LAMS(K)**(BS+3.)
2100	NPSACWS(K) = CONS13ASN(K)NC3D(K)RHO(K) &
2101	N0S(K)/ &
2102	LAMS(K)**(BS+3.)
2103
2104	END IF
2105
2106	!............................................................................
2107	! COLLECTION OF CLOUD WATER BY GRAUPEL
2108
2109	IF (QG3D(K).GE.1.E-8 .AND. QC3D(K).GE.QSMALL) THEN
2110
2111	PSACWG(K) = CONS14AGN(K)QC3D(K)RHO(K) &
2112	N0G(K)/ &
2113	LAMG(K)**(BG+3.)
2114	NPSACWG(K) = CONS14AGN(K)NC3D(K)RHO(K) &
2115	N0G(K)/ &
2116	LAMG(K)**(BG+3.)
2117	END IF
2118
2119	!.......................................................................
2120	! HM, ADD 12/13/06
2121	! CLOUD ICE COLLECTING DROPLETS, ASSUME THAT CLOUD ICE MEAN DIAM > 100 MICRON
2122	! BEFORE RIMING CAN OCCUR
2123	! ASSUME THAT RIME COLLECTED ON CLOUD ICE DOES NOT LEAD
2124	! TO HALLET-MOSSOP SPLINTERING
2125
2126	IF (QI3D(K).GE.1.E-8 .AND. QC3D(K).GE.QSMALL) THEN
2127
2128	! PUT IN SIZE DEPENDENT COLLECTION EFFICIENCY BASED ON STOKES LAW
2129	! FROM THOMPSON ET AL. 2004, MWR
2130
2131	IF (1./LAMI(K).GE.100.E-6) THEN
2132
2133	PSACWI(K) = CONS16AIN(K)QC3D(K)RHO(K) &
2134	N0I(K)/ &
2135	LAMI(K)**(BI+3.)
2136	NPSACWI(K) = CONS16AIN(K)NC3D(K)RHO(K) &
2137	N0I(K)/ &
2138	LAMI(K)**(BI+3.)
2139	END IF
2140	END IF
2141
2142	!.......................................................................
2143	! ACCRETION OF RAIN WATER BY SNOW
2144	! FORMULA FROM IKAWA AND SAITO, 1991, USED BY REISNER ET AL, 1998
2145
2146	IF (QR3D(K).GE.1.E-8.AND.QNI3D(K).GE.1.E-8) THEN
2147
2148	UMS = ASN(K)CONS3/(LAMS(K)*BS)
2149	UMR = ARN(K)CONS4/(LAMR(K)*BR)
2150	UNS = ASN(K)CONS5/LAMS(K)*BS
2151	UNR = ARN(K)CONS6/LAMR(K)*BR
2152
2153	! SET REASLISTIC LIMITS ON FALLSPEEDS
2154	UMS=MIN(UMS,1.2)
2155	UNS=MIN(UNS,1.2)
2156	UMR=MIN(UMR,9.1)
2157	UNR=MIN(UNR,9.1)
2158
2159	PRACS(K) = CONS41(((1.2UMR-0.95UMS)*2+ &
2160	0.08UMSUMR)*0.5RHO(K)* &
2161	N0RR(K)N0S(K)/LAMR(K)3 &
2162	(5./(LAMR(K)*3LAMS(K))+ &
2163	2./(LAMR(K)*2LAMS(K)**2)+ &
2164	0.5/(LAMR(k)LAMS(k)*3)))
2165
2166	NPRACS(K) = CONS32RHO(K)(1.7(UNR-UNS)*2+ &
2167	0.3UNRUNS)*0.5N0RR(K)N0S(K) &
2168	(1./(LAMR(K)*3LAMS(K))+ &
2169	1./(LAMR(K)*2LAMS(K)**2)+ &
2170	1./(LAMR(K)LAMS(K)*3))
2171
2172	! MAKE SURE PRACS DOESN'T EXCEED TOTAL RAIN MIXING RATIO
2173	! AS THIS MAY OTHERWISE RESULT IN TOO MUCH TRANSFER OF WATER DURING
2174	! RIME-SPLINTERING
2175
2176	PRACS(K) = MIN(PRACS(K),QR3D(K)/DT)
2177
2178	! COLLECTION OF SNOW BY RAIN - NEEDED FOR GRAUPEL CONVERSION CALCULATIONS
2179	! ONLY CALCULATE IF SNOW AND RAIN MIXING RATIOS EXCEED 0.1 G/KG
2180
2181	! ASSUME COLLECTION OF SNOW BY RAIN PRODUCES GRAUPEL NOT HAIL
2182
2183	IF (IHAIL.EQ.0) THEN
2184	IF (QNI3D(K).GE.0.1E-3.AND.QR3D(K).GE.0.1E-3) THEN
2185	PSACR(K) = CONS31(((1.2UMR-0.95UMS)*2+ &
2186	0.08UMSUMR)*0.5RHO(K)* &
2187	N0RR(K)N0S(K)/LAMS(K)3 &
2188	(5./(LAMS(K)*3LAMR(K))+ &
2189	2./(LAMS(K)*2LAMR(K)**2)+ &
2190	0.5/(LAMS(K)LAMR(K)*3)))
2191	END IF
2192	END IF
2193
2194	END IF
2195
2196	!.......................................................................
2197
2198	! COLLECTION OF RAINWATER BY GRAUPEL, FROM IKAWA AND SAITO 1990,
2199	! USED BY REISNER ET AL 1998
2200	IF (QR3D(K).GE.1.E-8.AND.QG3D(K).GE.1.E-8) THEN
2201
2202	UMG = AGN(K)CONS7/(LAMG(K)*BG)
2203	UMR = ARN(K)CONS4/(LAMR(K)*BR)
2204	UNG = AGN(K)CONS8/LAMG(K)*BG
2205	UNR = ARN(K)CONS6/LAMR(K)*BR
2206
2207	! SET REASLISTIC LIMITS ON FALLSPEEDS
2208	UMG=MIN(UMG,20.)
2209	UNG=MIN(UNG,20.)
2210	UMR=MIN(UMR,9.1)
2211	UNR=MIN(UNR,9.1)
2212
2213	PRACG(K) = CONS41(((1.2UMR-0.95UMG)*2+ &
2214	0.08UMGUMR)*0.5RHO(K)* &
2215	N0RR(K)N0G(K)/LAMR(K)3 &
2216	(5./(LAMR(K)*3LAMG(K))+ &
2217	2./(LAMR(K)*2LAMG(K)**2)+ &
2218	0.5/(LAMR(k)LAMG(k)*3)))
2219
2220	NPRACG(K) = CONS32RHO(K)(1.7(UNR-UNG)*2+ &
2221	0.3UNRUNG)*0.5N0RR(K)N0G(K) &
2222	(1./(LAMR(K)*3LAMG(K))+ &
2223	1./(LAMR(K)*2LAMG(K)**2)+ &
2224	1./(LAMR(K)LAMG(K)*3))
2225
2226	! MAKE SURE PRACG DOESN'T EXCEED TOTAL RAIN MIXING RATIO
2227	! AS THIS MAY OTHERWISE RESULT IN TOO MUCH TRANSFER OF WATER DURING
2228	! RIME-SPLINTERING
2229
2230	PRACG(K) = MIN(PRACG(K),QR3D(K)/DT)
2231
2232	END IF
2233
2234	!.......................................................................
2235	! RIME-SPLINTERING - SNOW
2236	! HALLET-MOSSOP (1974)
2237	! NUMBER OF SPLINTERS FORMED IS BASED ON MASS OF RIMED WATER
2238
2239	! DUM1 = MASS OF INDIVIDUAL SPLINTERS
2240
2241	! HM ADD THRESHOLD SNOW AND DROPLET MIXING RATIO FOR RIME-SPLINTERING
2242	! TO LIMIT RIME-SPLINTERING IN STRATIFORM CLOUDS
2243	! THESE THRESHOLDS CORRESPOND WITH GRAUPEL THRESHOLDS IN RH 1984
2244
2245	IF (QNI3D(K).GE.0.1E-3.AND.QC3D(K).GE.0.5E-3) THEN
2246	IF (PSACWS(K).GT.0..OR.PRACS(K).GT.0.) THEN
2247	IF (T3D(K).LT.270.16 .AND. T3D(K).GT.265.16) THEN
2248
2249	IF (T3D(K).GT.270.16) THEN
2250	FMULT = 0.
2251	ELSE IF (T3D(K).LE.270.16.AND.T3D(K).GT.268.16) THEN
2252	FMULT = (270.16-T3D(K))/2.
2253	ELSE IF (T3D(K).GE.265.16.AND.T3D(K).LE.268.16) THEN
2254	FMULT = (T3D(K)-265.16)/3.
2255	ELSE IF (T3D(K).LT.265.16) THEN
2256	FMULT = 0.
2257	END IF
2258
2259	! 1000 IS TO CONVERT FROM KG TO G
2260
2261	! SPLINTERING FROM DROPLETS ACCRETED ONTO SNOW
2262
2263	IF (PSACWS(K).GT.0.) THEN
2264	NMULTS(K) = 35.E4PSACWS(K)FMULT*1000.
2265	QMULTS(K) = NMULTS(K)*MMULT
2266
2267	! CONSTRAIN SO THAT TRANSFER OF MASS FROM SNOW TO ICE CANNOT BE MORE MASS
2268	! THAN WAS RIMED ONTO SNOW
2269
2270	QMULTS(K) = MIN(QMULTS(K),PSACWS(K))
2271	PSACWS(K) = PSACWS(K)-QMULTS(K)
2272
2273	END IF
2274
2275	! RIMING AND SPLINTERING FROM ACCRETED RAINDROPS
2276
2277	IF (PRACS(K).GT.0.) THEN
2278	NMULTR(K) = 35.E4PRACS(K)FMULT*1000.
2279	QMULTR(K) = NMULTR(K)*MMULT
2280
2281	! CONSTRAIN SO THAT TRANSFER OF MASS FROM SNOW TO ICE CANNOT BE MORE MASS
2282	! THAN WAS RIMED ONTO SNOW
2283
2284	QMULTR(K) = MIN(QMULTR(K),PRACS(K))
2285
2286	PRACS(K) = PRACS(K)-QMULTR(K)
2287
2288	END IF
2289
2290	END IF
2291	END IF
2292	END IF
2293
2294	!.......................................................................
2295	! RIME-SPLINTERING - GRAUPEL
2296	! HALLET-MOSSOP (1974)
2297	! NUMBER OF SPLINTERS FORMED IS BASED ON MASS OF RIMED WATER
2298
2299	! DUM1 = MASS OF INDIVIDUAL SPLINTERS
2300
2301	! HM ADD THRESHOLD SNOW MIXING RATIO FOR RIME-SPLINTERING
2302	! TO LIMIT RIME-SPLINTERING IN STRATIFORM CLOUDS
2303
2304	! ONLY CALCULATE FOR GRAUPEL NOT HAIL
2305	IF (IHAIL.EQ.0) THEN
2306	IF (PSACWG(K).GT.0..OR.PRACG(K).GT.0.) THEN
2307	IF (T3D(K).LT.270.16 .AND. T3D(K).GT.265.16) THEN
2308
2309	IF (T3D(K).GT.270.16) THEN
2310	FMULT = 0.
2311	ELSE IF (T3D(K).LE.270.16.AND.T3D(K).GT.268.16) THEN
2312	FMULT = (270.16-T3D(K))/2.
2313	ELSE IF (T3D(K).GE.265.16.AND.T3D(K).LE.268.16) THEN
2314	FMULT = (T3D(K)-265.16)/3.
2315	ELSE IF (T3D(K).LT.265.16) THEN
2316	FMULT = 0.
2317	END IF
2318
2319	! 1000 IS TO CONVERT FROM KG TO G
2320
2321	! SPLINTERING FROM DROPLETS ACCRETED ONTO GRAUPEL
2322
2323	IF (PSACWG(K).GT.0.) THEN
2324	NMULTG(K) = 35.E4PSACWG(K)FMULT*1000.
2325	QMULTG(K) = NMULTG(K)*MMULT
2326
2327	! CONSTRAIN SO THAT TRANSFER OF MASS FROM GRAUPEL TO ICE CANNOT BE MORE MASS
2328	! THAN WAS RIMED ONTO GRAUPEL
2329
2330	QMULTG(K) = MIN(QMULTG(K),PSACWG(K))
2331	PSACWG(K) = PSACWG(K)-QMULTG(K)
2332
2333	END IF
2334
2335	! RIMING AND SPLINTERING FROM ACCRETED RAINDROPS
2336
2337	IF (PRACG(K).GT.0.) THEN
2338	NMULTRG(K) = 35.E4PRACG(K)FMULT*1000.
2339	QMULTRG(K) = NMULTRG(K)*MMULT
2340
2341	! CONSTRAIN SO THAT TRANSFER OF MASS FROM GRAUPEL TO ICE CANNOT BE MORE MASS
2342	! THAN WAS RIMED ONTO GRAUPEL
2343
2344	QMULTRG(K) = MIN(QMULTRG(K),PRACG(K))
2345	PRACG(K) = PRACG(K)-QMULTRG(K)
2346
2347	END IF
2348
2349	END IF
2350	END IF
2351	END IF
2352
2353	!........................................................................
2354	! CONVERSION OF RIMED CLOUD WATER ONTO SNOW TO GRAUPEL
2355	! ASSUME CONVERTED SNOW FORMS GRAUPEL NOT HAIL
2356	! HAIL ASSUMED TO ONLY FORM BY FREEZING OF RAIN
2357	! OR COLLISIONS OF RAIN WITH CLOUD ICE
2358
2359	IF (IHAIL.EQ.0) THEN
2360	IF (PSACWS(K).GT.0.) THEN
2361	! ONLY ALLOW CONVERSION IF QNI > 0.1 AND QC > 0.5 G/KG FOLLOWING RUTLEDGE AND HOBBS (1984)
2362	IF (QNI3D(K).GE.0.1E-3.AND.QC3D(K).GE.0.5E-3) THEN
2363
2364	! PORTION OF RIMING CONVERTED TO GRAUPEL (REISNER ET AL. 1998, ORIGINALLY IS1991)
2365	PGSACW(K) = MIN(PSACWS(K),CONS17DTN0S(K)QC3D(K)QC3D(K)* &
2366	ASN(K)*ASN(K)/ &
2367	(RHO(K)LAMS(K)(2.BS+2.)))
2368
2369	! MIX RAT CONVERTED INTO GRAUPEL AS EMBRYO (REISNER ET AL. 1998, ORIG M1990)
2370	DUM = MAX(RHOSN/(RHOG-RHOSN)*PGSACW(K),0.)
2371
2372	! NUMBER CONCENTRAITON OF EMBRYO GRAUPEL FROM RIMING OF SNOW
2373	NSCNG(K) = DUM/MG0*RHO(K)
2374	! LIMIT MAX NUMBER CONVERTED TO SNOW NUMBER
2375	NSCNG(K) = MIN(NSCNG(K),NS3D(K)/DT)
2376
2377	! PORTION OF RIMING LEFT FOR SNOW
2378	PSACWS(K) = PSACWS(K) - PGSACW(K)
2379	END IF
2380	END IF
2381
2382	! CONVERSION OF RIMED RAINWATER ONTO SNOW CONVERTED TO GRAUPEL
2383
2384	IF (PRACS(K).GT.0.) THEN
2385	! ONLY ALLOW CONVERSION IF QNI > 0.1 AND QR > 0.1 G/KG FOLLOWING RUTLEDGE AND HOBBS (1984)
2386	IF (QNI3D(K).GE.0.1E-3.AND.QR3D(K).GE.0.1E-3) THEN
2387	! PORTION OF COLLECTED RAINWATER CONVERTED TO GRAUPEL (REISNER ET AL. 1998)
2388	DUM = CONS18(4./LAMS(K))3(4./LAMS(K))**3 &
2389	/(CONS18(4./LAMS(K))3(4./LAMS(K))**3+ &
2390	CONS19(4./LAMR(K))3(4./LAMR(K))**3)
2391	DUM=MIN(DUM,1.)
2392	DUM=MAX(DUM,0.)
2393	PGRACS(K) = (1.-DUM)*PRACS(K)
2394	NGRACS(K) = (1.-DUM)*NPRACS(K)
2395	! LIMIT MAX NUMBER CONVERTED TO MIN OF EITHER RAIN OR SNOW NUMBER CONCENTRATION
2396	NGRACS(K) = MIN(NGRACS(K),NR3D(K)/DT)
2397	NGRACS(K) = MIN(NGRACS(K),NS3D(K)/DT)
2398
2399	! AMOUNT LEFT FOR SNOW PRODUCTION
2400	PRACS(K) = PRACS(K) - PGRACS(K)
2401	NPRACS(K) = NPRACS(K) - NGRACS(K)
2402	! CONVERSION TO GRAUPEL DUE TO COLLECTION OF SNOW BY RAIN
2403	PSACR(K)=PSACR(K)*(1.-DUM)
2404	END IF
2405	END IF
2406	END IF
2407
2408	!.......................................................................
2409	! FREEZING OF RAIN DROPS
2410	! FREEZING ALLOWED BELOW -4 C
2411
2412	IF (T3D(K).LT.269.15.AND.QR3D(K).GE.QSMALL) THEN
2413
2414	! IMMERSION FREEZING (BIGG 1953)
2415	MNUCCR(K) = CONS20NR3D(K)EXP(AIMM(273.15-T3D(K)))/LAMR(K)*3 &
2416	/LAMR(K)**3
2417
2418	NNUCCR(K) = PINR3D(K)BIMMEXP(AIMM(273.15-T3D(K)))/LAMR(K)**3
2419
2420	! PREVENT DIVERGENCE BETWEEN MIXING RATIO AND NUMBER CONC
2421	NNUCCR(K) = MIN(NNUCCR(K),NR3D(K)/DT)
2422
2423	END IF
2424
2425	!.......................................................................
2426	! ACCRETION OF CLOUD LIQUID WATER BY RAIN
2427	! CONTINUOUS COLLECTION EQUATION WITH
2428	! GRAVITATIONAL COLLECTION KERNEL, DROPLET FALL SPEED NEGLECTED
2429
2430	IF (QR3D(K).GE.1.E-8 .AND. QC3D(K).GE.1.E-8) THEN
2431
2432	! 12/13/06 HM ADD, REPLACE WITH NEWER FORMULA FROM
2433	! KHAIROUTDINOV AND KOGAN 2000, MWR
2434
2435	DUM=(QC3D(K)*QR3D(K))
2436	PRA(K) = 67.(DUM)*1.15
2437	NPRA(K) = PRA(K)/(QC3D(K)/NC3D(K))
2438
2439	END IF
2440	!.......................................................................
2441	! SELF-COLLECTION OF RAIN DROPS
2442	! FROM BEHENG(1994)
2443	! FROM NUMERICAL SIMULATION OF THE STOCHASTIC COLLECTION EQUATION
2444	! AS DESCRINED ABOVE FOR AUTOCONVERSION
2445
2446	IF (QR3D(K).GE.1.E-8) THEN
2447	NRAGG(K) = -8.NR3D(K)QR3D(K)*RHO(K)
2448	END IF
2449
2450	!.......................................................................
2451	! AUTOCONVERSION OF CLOUD ICE TO SNOW
2452	! FOLLOWING HARRINGTON ET AL. (1995) WITH MODIFICATION
2453	! HERE IT IS ASSUMED THAT AUTOCONVERSION CAN ONLY OCCUR WHEN THE
2454	! ICE IS GROWING, I.E. IN CONDITIONS OF ICE SUPERSATURATION
2455
2456	IF (QI3D(K).GE.1.E-8 .AND.QVQVSI(K).GE.1.) THEN
2457
2458	! COFFI = 2./LAMI(K)
2459	! IF (COFFI.GE.DCS) THEN
2460	NPRCI(K) = CONS21(QV3D(K)-QVI(K))RHO(K) &
2461	N0I(K)EXP(-LAMI(K)DCS)DV(K)/ABI(K)
2462	PRCI(K) = CONS22*NPRCI(K)
2463	NPRCI(K) = MIN(NPRCI(K),NI3D(K)/DT)
2464
2465	! END IF
2466	END IF
2467
2468	!.......................................................................
2469	! ACCRETION OF CLOUD ICE BY SNOW
2470	! FOR THIS CALCULATION, IT IS ASSUMED THAT THE VS >> VI
2471	! AND DS >> DI FOR CONTINUOUS COLLECTION
2472
2473	IF (QNI3D(K).GE.1.E-8 .AND. QI3D(K).GE.QSMALL) THEN
2474	PRAI(K) = CONS23ASN(K)QI3D(K)RHO(K)N0S(K)/ &
2475	LAMS(K)**(BS+3.)
2476	NPRAI(K) = CONS23ASN(K)NI3D(K)* &
2477	RHO(K)*N0S(K)/ &
2478	LAMS(K)**(BS+3.)
2479	NPRAI(K)=MIN(NPRAI(K),NI3D(K)/DT)
2480	END IF
2481
2482	!.......................................................................
2483	! HM, ADD 12/13/06, COLLISION OF RAIN AND ICE TO PRODUCE SNOW OR GRAUPEL
2484	! FOLLOWS REISNER ET AL. 1998
2485	! ASSUMED FALLSPEED AND SIZE OF ICE CRYSTAL << THAN FOR RAIN
2486
2487	IF (QR3D(K).GE.1.E-8.AND.QI3D(K).GE.1.E-8.AND.T3D(K).LE.273.15) THEN
2488
2489	! ALLOW GRAUPEL FORMATION FROM RAIN-ICE COLLISIONS ONLY IF RAIN MIXING RATIO > 0.1 G/KG,
2490	! OTHERWISE ADD TO SNOW
2491
2492	IF (QR3D(K).GE.0.1E-3) THEN
2493	NIACR(K)=CONS24NI3D(K)N0RR(K)*ARN(K) &
2494	/LAMR(K)*(BR+3.)RHO(K)
2495	PIACR(K)=CONS25NI3D(K)N0RR(K)*ARN(K) &
2496	/LAMR(K)(BR+3.)/LAMR(K)3*RHO(K)
2497	PRACI(K)=CONS24QI3D(K)N0RR(K)*ARN(K)/ &
2498	LAMR(K)*(BR+3.)RHO(K)
2499	NIACR(K)=MIN(NIACR(K),NR3D(K)/DT)
2500	NIACR(K)=MIN(NIACR(K),NI3D(K)/DT)
2501	ELSE
2502	NIACRS(K)=CONS24NI3D(K)N0RR(K)*ARN(K) &
2503	/LAMR(K)*(BR+3.)RHO(K)
2504	PIACRS(K)=CONS25NI3D(K)N0RR(K)*ARN(K) &
2505	/LAMR(K)(BR+3.)/LAMR(K)3*RHO(K)
2506	PRACIS(K)=CONS24QI3D(K)N0RR(K)*ARN(K)/ &
2507	LAMR(K)*(BR+3.)RHO(K)
2508	NIACRS(K)=MIN(NIACRS(K),NR3D(K)/DT)
2509	NIACRS(K)=MIN(NIACRS(K),NI3D(K)/DT)
2510	END IF
2511	END IF
2512
2513	!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
2514	! NUCLEATION OF CLOUD ICE FROM HOMOGENEOUS AND HETEROGENEOUS FREEZING ON AEROSOL
2515
2516	IF (INUC.EQ.0) THEN
2517
2518	! FREEZING OF AEROSOL ONLY ALLOWED BELOW -5 C
2519	! AND ABOVE DELIQUESCENCE THRESHOLD OF 80%
2520	! AND ABOVE ICE SATURATION
2521
2522	! add threshold according to Greg Thomspon
2523
2524	if ((QVQVS(K).GE.0.999.and.T3D(K).le.265.15).or. &
2525	QVQVSI(K).ge.1.08) then
2526
2527	! hm, modify dec. 5, 2006, replace with cooper curve
2528	kc2 = 0.005exp(0.304(273.15-T3D(K)))*1000. ! convert from L-1 to m-3
2529	! limit to 500 L-1
2530	kc2 = min(kc2,500.e3)
2531	kc2=MAX(kc2/rho(k),0.) ! convert to kg-1
2532
2533	IF (KC2.GT.NI3D(K)+NS3D(K)+NG3D(K)) THEN
2534	NNUCCD(K) = (KC2-NI3D(K)-NS3D(K)-NG3D(K))/DT
2535	MNUCCD(K) = NNUCCD(K)*MI0
2536	END IF
2537
2538	END IF
2539
2540	ELSE IF (INUC.EQ.1) THEN
2541
2542	IF (T3D(K).LT.273.15.AND.QVQVSI(K).GT.1.) THEN
2543
2544	KC2 = 0.16*1000./RHO(K) ! CONVERT FROM L-1 TO KG-1
2545	IF (KC2.GT.NI3D(K)+NS3D(K)+NG3D(K)) THEN
2546	NNUCCD(K) = (KC2-NI3D(K)-NS3D(K)-NG3D(K))/DT
2547	MNUCCD(K) = NNUCCD(K)*MI0
2548	END IF
2549	END IF
2550
2551	END IF
2552
2553	!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
2554
2555	101 CONTINUE
2556
2557	!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
2558	! CALCULATE EVAP/SUB/DEP TERMS FOR QI,QNI,QR
2559
2560	! NO VENTILATION FOR CLOUD ICE
2561
2562	IF (QI3D(K).GE.QSMALL) THEN
2563
2564	EPSI = 2.PIN0I(K)RHO(K)DV(K)/(LAMI(K)*LAMI(K))
2565
2566	ELSE
2567	EPSI = 0.
2568	END IF
2569
2570	IF (QNI3D(K).GE.QSMALL) THEN
2571	EPSS = 2.PIN0S(K)RHO(K)DV(K)* &
2572	(F1S/(LAMS(K)*LAMS(K))+ &
2573	F2S(ASN(K)RHO(K)/MU(K))*0.5 &
2574	SC(K)*(1./3.)CONS10/ &
2575	(LAMS(K)**CONS35))
2576	ELSE
2577	EPSS = 0.
2578	END IF
2579
2580	IF (QG3D(K).GE.QSMALL) THEN
2581	EPSG = 2.PIN0G(K)RHO(K)DV(K)* &
2582	(F1S/(LAMG(K)*LAMG(K))+ &
2583	F2S(AGN(K)RHO(K)/MU(K))*0.5 &
2584	SC(K)*(1./3.)CONS11/ &
2585	(LAMG(K)**CONS36))
2586
2587
2588	ELSE
2589	EPSG = 0.
2590	END IF
2591
2592	IF (QR3D(K).GE.QSMALL) THEN
2593	EPSR = 2.PIN0RR(K)RHO(K)DV(K)* &
2594	(F1R/(LAMR(K)*LAMR(K))+ &
2595	F2R(ARN(K)RHO(K)/MU(K))*0.5 &
2596	SC(K)*(1./3.)CONS9/ &
2597	(LAMR(K)**CONS34))
2598	ELSE
2599	EPSR = 0.
2600	END IF
2601
2602	! ONLY INCLUDE REGION OF ICE SIZE DIST < DCS
2603	! DUM IS FRACTION OF D*N(D) < DCS
2604
2605	! LOGIC BELOW FOLLOWS THAT OF HARRINGTON ET AL. 1995 (JAS)
2606	IF (QI3D(K).GE.QSMALL) THEN
2607	DUM=(1.-EXP(-LAMI(K)DCS)(1.+LAMI(K)*DCS))
2608	PRD(K) = EPSI(QV3D(K)-QVI(K))/ABI(K)DUM
2609	ELSE
2610	DUM=0.
2611	END IF
2612	! ADD DEPOSITION IN TAIL OF ICE SIZE DIST TO SNOW IF SNOW IS PRESENT
2613	IF (QNI3D(K).GE.QSMALL) THEN
2614	PRDS(K) = EPSS*(QV3D(K)-QVI(K))/ABI(K)+ &
2615	EPSI(QV3D(K)-QVI(K))/ABI(K)(1.-DUM)
2616	! OTHERWISE ADD TO CLOUD ICE
2617	ELSE
2618	PRD(K) = PRD(K)+EPSI(QV3D(K)-QVI(K))/ABI(K)(1.-DUM)
2619	END IF
2620	! VAPOR DPEOSITION ON GRAUPEL
2621	PRDG(K) = EPSG*(QV3D(K)-QVI(K))/ABI(K)
2622
2623	! NO CONDENSATION ONTO RAIN, ONLY EVAP
2624
2625	IF (QV3D(K).LT.QVS(K)) THEN
2626	PRE(K) = EPSR*(QV3D(K)-QVS(K))/AB(K)
2627	PRE(K) = MIN(PRE(K),0.)
2628	ELSE
2629	PRE(K) = 0.
2630	END IF
2631
2632	! MAKE SURE NOT PUSHED INTO ICE SUPERSAT/SUBSAT
2633	! FORMULA FROM REISNER 2 SCHEME
2634
2635	DUM = (QV3D(K)-QVI(K))/DT
2636
2637	FUDGEF = 0.9999
2638	SUM_DEP = PRD(K)+PRDS(K)+MNUCCD(K)+PRDG(K)
2639
2640	IF( (DUM.GT.0. .AND. SUM_DEP.GT.DUM*FUDGEF) .OR. &
2641	(DUM.LT.0. .AND. SUM_DEP.LT.DUM*FUDGEF) ) THEN
2642	MNUCCD(K) = FUDGEFMNUCCD(K)DUM/SUM_DEP
2643	PRD(K) = FUDGEFPRD(K)DUM/SUM_DEP
2644	PRDS(K) = FUDGEFPRDS(K)DUM/SUM_DEP
2645	PRDG(K) = FUDGEFPRDG(K)DUM/SUM_DEP
2646	ENDIF
2647
2648	! IF CLOUD ICE/SNOW/GRAUPEL VAP DEPOSITION IS NEG, THEN ASSIGN TO SUBLIMATION PROCESSES
2649
2650	IF (PRD(K).LT.0.) THEN
2651	EPRD(K)=PRD(K)
2652	PRD(K)=0.
2653	END IF
2654	IF (PRDS(K).LT.0.) THEN
2655	EPRDS(K)=PRDS(K)
2656	PRDS(K)=0.
2657	END IF
2658	IF (PRDG(K).LT.0.) THEN
2659	EPRDG(K)=PRDG(K)
2660	PRDG(K)=0.
2661	END IF
2662
2663	!.......................................................................
2664	!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
2665
2666	! CONSERVATION OF WATER
2667	! THIS IS ADOPTED LOOSELY FROM MM5 RESINER CODE. HOWEVER, HERE WE
2668	! ONLY ADJUST PROCESSES THAT ARE NEGATIVE, RATHER THAN ALL PROCESSES.
2669
2670	! IF MIXING RATIOS LESS THAN QSMALL, THEN NO DEPLETION OF WATER
2671	! THROUGH MICROPHYSICAL PROCESSES, SKIP CONSERVATION
2672
2673	! NOTE: CONSERVATION CHECK NOT APPLIED TO NUMBER CONCENTRATION SPECIES. ADDITIONAL CATCH
2674	! BELOW WILL PREVENT NEGATIVE NUMBER CONCENTRATION
2675	! FOR EACH MICROPHYSICAL PROCESS WHICH PROVIDES A SOURCE FOR NUMBER, THERE IS A CHECK
2676	! TO MAKE SURE THAT CAN'T EXCEED TOTAL NUMBER OF DEPLETED SPECIES WITH THE TIME
2677	! STEP
2678
2679	!****SENSITIVITY - NO ICE
2680
2681	IF (ILIQ.EQ.1) THEN
2682	MNUCCC(K)=0.
2683	NNUCCC(K)=0.
2684	MNUCCR(K)=0.
2685	NNUCCR(K)=0.
2686	MNUCCD(K)=0.
2687	NNUCCD(K)=0.
2688	END IF
2689
2690	! ****SENSITIVITY - NO GRAUPEL
2691	IF (IGRAUP.EQ.1) THEN
2692	PRACG(K) = 0.
2693	PSACR(K) = 0.
2694	PSACWG(K) = 0.
2695	PGSACW(K) = 0.
2696	PGRACS(K) = 0.
2697	PRDG(K) = 0.
2698	EPRDG(K) = 0.
2699	EVPMG(K) = 0.
2700	PGMLT(K) = 0.
2701	NPRACG(K) = 0.
2702	NPSACWG(K) = 0.
2703	NSCNG(K) = 0.
2704	NGRACS(K) = 0.
2705	NSUBG(K) = 0.
2706	NGMLTG(K) = 0.
2707	NGMLTR(K) = 0.
2708	END IF
2709
2710	! CONSERVATION OF QC
2711
2712	DUM = (PRC(K)+PRA(K)+MNUCCC(K)+PSACWS(K)+PSACWI(K)+QMULTS(K)+PSACWG(K)+PGSACW(K)+QMULTG(K))*DT
2713
2714	IF (DUM.GT.QC3D(K).AND.QC3D(K).GE.QSMALL) THEN
2715	RATIO = QC3D(K)/DUM
2716
2717	PRC(K) = PRC(K)*RATIO
2718	PRA(K) = PRA(K)*RATIO
2719	MNUCCC(K) = MNUCCC(K)*RATIO
2720	PSACWS(K) = PSACWS(K)*RATIO
2721	PSACWI(K) = PSACWI(K)*RATIO
2722	QMULTS(K) = QMULTS(K)*RATIO
2723	QMULTG(K) = QMULTG(K)*RATIO
2724	PSACWG(K) = PSACWG(K)*RATIO
2725	PGSACW(K) = PGSACW(K)*RATIO
2726	END IF
2727
2728	! CONSERVATION OF QI
2729
2730	DUM = (-PRD(K)-MNUCCC(K)+PRCI(K)+PRAI(K)-QMULTS(K)-QMULTG(K)-QMULTR(K)-QMULTRG(K) &
2731	-MNUCCD(K)+PRACI(K)+PRACIS(K)-EPRD(K)-PSACWI(K))*DT
2732
2733	IF (DUM.GT.QI3D(K).AND.QI3D(K).GE.QSMALL) THEN
2734
2735	RATIO = (QI3D(K)/DT+PRD(K)+MNUCCC(K)+QMULTS(K)+QMULTG(K)+QMULTR(K)+QMULTRG(K)+ &
2736	MNUCCD(K)+PSACWI(K))/ &
2737	(PRCI(K)+PRAI(K)+PRACI(K)+PRACIS(K)-EPRD(K))
2738
2739	PRCI(K) = PRCI(K)*RATIO
2740	PRAI(K) = PRAI(K)*RATIO
2741	PRACI(K) = PRACI(K)*RATIO
2742	PRACIS(K) = PRACIS(K)*RATIO
2743	EPRD(K) = EPRD(K)*RATIO
2744
2745	END IF
2746
2747	! CONSERVATION OF QR
2748
2749	DUM=((PRACS(K)-PRE(K))+(QMULTR(K)+QMULTRG(K)-PRC(K))+(MNUCCR(K)-PRA(K))+ &
2750	PIACR(K)+PIACRS(K)+PGRACS(K)+PRACG(K))*DT
2751
2752	IF (DUM.GT.QR3D(K).AND.QR3D(K).GE.QSMALL) THEN
2753
2754	RATIO = (QR3D(K)/DT+PRC(K)+PRA(K))/ &
2755	(-PRE(K)+QMULTR(K)+QMULTRG(K)+PRACS(K)+MNUCCR(K)+PIACR(K)+PIACRS(K)+PGRACS(K)+PRACG(K))
2756
2757	PRE(K) = PRE(K)*RATIO
2758	PRACS(K) = PRACS(K)*RATIO
2759	QMULTR(K) = QMULTR(K)*RATIO
2760	QMULTRG(K) = QMULTRG(K)*RATIO
2761	MNUCCR(K) = MNUCCR(K)*RATIO
2762	PIACR(K) = PIACR(K)*RATIO
2763	PIACRS(K) = PIACRS(K)*RATIO
2764	PGRACS(K) = PGRACS(K)*RATIO
2765	PRACG(K) = PRACG(K)*RATIO
2766
2767	END IF
2768
2769	! CONSERVATION OF QNI
2770	! CONSERVATION FOR GRAUPEL SCHEME
2771
2772	IF (IGRAUP.EQ.0) THEN
2773
2774	DUM = (-PRDS(K)-PSACWS(K)-PRAI(K)-PRCI(K)-PRACS(K)-EPRDS(K)+PSACR(K)-PIACRS(K)-PRACIS(K))*DT
2775
2776	IF (DUM.GT.QNI3D(K).AND.QNI3D(K).GE.QSMALL) THEN
2777
2778	RATIO = (QNI3D(K)/DT+PRDS(K)+PSACWS(K)+PRAI(K)+PRCI(K)+PRACS(K)+PIACRS(K)+PRACIS(K))/(-EPRDS(K)+PSACR(K))
2779
2780	EPRDS(K) = EPRDS(K)*RATIO
2781	PSACR(K) = PSACR(K)*RATIO
2782
2783	END IF
2784
2785	! FOR NO GRAUPEL, NEED TO INCLUDE FREEZING OF RAIN FOR SNOW
2786	ELSE IF (IGRAUP.EQ.1) THEN
2787
2788	DUM = (-PRDS(K)-PSACWS(K)-PRAI(K)-PRCI(K)-PRACS(K)-EPRDS(K)+PSACR(K)-PIACRS(K)-PRACIS(K)-MNUCCR(K))*DT
2789
2790	IF (DUM.GT.QNI3D(K).AND.QNI3D(K).GE.QSMALL) THEN
2791
2792	RATIO = (QNI3D(K)/DT+PRDS(K)+PSACWS(K)+PRAI(K)+PRCI(K)+PRACS(K)+PIACRS(K)+PRACIS(K)+MNUCCR(K))/(-EPRDS(K)+PSACR(K))
2793
2794	EPRDS(K) = EPRDS(K)*RATIO
2795	PSACR(K) = PSACR(K)*RATIO
2796
2797	END IF
2798
2799	END IF
2800
2801	! CONSERVATION OF QG
2802
2803	DUM = (-PSACWG(K)-PRACG(K)-PGSACW(K)-PGRACS(K)-PRDG(K)-MNUCCR(K)-EPRDG(K)-PIACR(K)-PRACI(K)-PSACR(K))*DT
2804
2805	IF (DUM.GT.QG3D(K).AND.QG3D(K).GE.QSMALL) THEN
2806
2807	RATIO = (QG3D(K)/DT+PSACWG(K)+PRACG(K)+PGSACW(K)+PGRACS(K)+PRDG(K)+MNUCCR(K)+PSACR(K)+&
2808	PIACR(K)+PRACI(K))/(-EPRDG(K))
2809
2810	EPRDG(K) = EPRDG(K)*RATIO
2811
2812	END IF
2813
2814	! TENDENCIES
2815
2816	QV3DTEN(K) = QV3DTEN(K)+(-PRE(K)-PRD(K)-PRDS(K)-MNUCCD(K)-EPRD(K)-EPRDS(K)-PRDG(K)-EPRDG(K))
2817	T3DTEN(K) = T3DTEN(K)+(PRE(K) &
2818	*XXLV(K)+(PRD(K)+PRDS(K)+ &
2819	MNUCCD(K)+EPRD(K)+EPRDS(K)+PRDG(K)+EPRDG(K))*XXLS(K)+ &
2820	(PSACWS(K)+PSACWI(K)+MNUCCC(K)+MNUCCR(K)+ &
2821	QMULTS(K)+QMULTG(K)+QMULTR(K)+QMULTRG(K)+PRACS(K) &
2822	+PSACWG(K)+PRACG(K)+PGSACW(K)+PGRACS(K))*XLF(K))/CPM(K)
2823
2824	QC3DTEN(K) = QC3DTEN(K)+ &
2825	(-PRA(K)-PRC(K)-MNUCCC(K)+PCC(K)- &
2826	PSACWS(K)-PSACWI(K)-QMULTS(K)-QMULTG(K)-PSACWG(K)-PGSACW(K))
2827	QI3DTEN(K) = QI3DTEN(K)+ &
2828	(PRD(K)+EPRD(K)+PSACWI(K)+MNUCCC(K)-PRCI(K)- &
2829	PRAI(K)+QMULTS(K)+QMULTG(K)+QMULTR(K)+QMULTRG(K)+MNUCCD(K)-PRACI(K)-PRACIS(K))
2830	QR3DTEN(K) = QR3DTEN(K)+ &
2831	(PRE(K)+PRA(K)+PRC(K)-PRACS(K)-MNUCCR(K)-QMULTR(K)-QMULTRG(K) &
2832	-PIACR(K)-PIACRS(K)-PRACG(K)-PGRACS(K))
2833
2834	IF (IGRAUP.EQ.0) THEN
2835
2836	QNI3DTEN(K) = QNI3DTEN(K)+ &
2837	(PRAI(K)+PSACWS(K)+PRDS(K)+PRACS(K)+PRCI(K)+EPRDS(K)-PSACR(K)+PIACRS(K)+PRACIS(K))
2838	NS3DTEN(K) = NS3DTEN(K)+(NSAGG(K)+NPRCI(K)-NSCNG(K)-NGRACS(K)+NIACRS(K))
2839	QG3DTEN(K) = QG3DTEN(K)+(PRACG(K)+PSACWG(K)+PGSACW(K)+PGRACS(K)+ &
2840	PRDG(K)+EPRDG(K)+MNUCCR(K)+PIACR(K)+PRACI(K)+PSACR(K))
2841	NG3DTEN(K) = NG3DTEN(K)+(NSCNG(K)+NGRACS(K)+NNUCCR(K)+NIACR(K))
2842
2843	! FOR NO GRAUPEL, NEED TO INCLUDE FREEZING OF RAIN FOR SNOW
2844	ELSE IF (IGRAUP.EQ.1) THEN
2845
2846	QNI3DTEN(K) = QNI3DTEN(K)+ &
2847	(PRAI(K)+PSACWS(K)+PRDS(K)+PRACS(K)+PRCI(K)+EPRDS(K)-PSACR(K)+PIACRS(K)+PRACIS(K)+MNUCCR(K))
2848	NS3DTEN(K) = NS3DTEN(K)+(NSAGG(K)+NPRCI(K)-NSCNG(K)-NGRACS(K)+NIACRS(K)+NNUCCR(K))
2849
2850	END IF
2851
2852	NC3DTEN(K) = NC3DTEN(K)+(-NNUCCC(K)-NPSACWS(K) &
2853	-NPRA(K)-NPRC(K)-NPSACWI(K)-NPSACWG(K))
2854
2855	NI3DTEN(K) = NI3DTEN(K)+ &
2856	(NNUCCC(K)-NPRCI(K)-NPRAI(K)+NMULTS(K)+NMULTG(K)+NMULTR(K)+NMULTRG(K)+ &
2857	NNUCCD(K)-NIACR(K)-NIACRS(K))
2858
2859	NR3DTEN(K) = NR3DTEN(K)+(NPRC1(K)-NPRACS(K)-NNUCCR(K) &
2860	+NRAGG(K)-NIACR(K)-NIACRS(K)-NPRACG(K)-NGRACS(K))
2861
2862	!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
2863	! NOW CALCULATE SATURATION ADJUSTMENT TO CONDENSE EXTRA VAPOR ABOVE
2864	! WATER SATURATION
2865
2866	DUMT = T3D(K)+DT*T3DTEN(K)
2867	DUMQV = QV3D(K)+DT*QV3DTEN(K)
2868	DUMQSS = 0.622*POLYSVP(DUMT,0)/ (PRES(K)-POLYSVP(DUMT,0))
2869	DUMQC = QC3D(K)+DT*QC3DTEN(K)
2870	DUMQC = MAX(DUMQC,0.)
2871
2872	! SATURATION ADJUSTMENT FOR LIQUID
2873
2874	DUMS = DUMQV-DUMQSS
2875	PCC(K) = DUMS/(1.+XXLV(K)*2DUMQSS/(CPM(K)RVDUMT**2))/DT
2876	IF (PCC(K)*DT+DUMQC.LT.0.) THEN
2877	PCC(K) = -DUMQC/DT
2878	END IF
2879
2880	QV3DTEN(K) = QV3DTEN(K)-PCC(K)
2881	T3DTEN(K) = T3DTEN(K)+PCC(K)*XXLV(K)/CPM(K)
2882	QC3DTEN(K) = QC3DTEN(K)+PCC(K)
2883
2884	!.......................................................................
2885	! ACTIVATION OF CLOUD DROPLETS
2886	! ACTIVATION OF DROPLET CURRENTLY NOT CALCULATED
2887	! DROPLET CONCENTRATION IS SPECIFIED !!!!!
2888
2889	!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
2890	! SUBLIMATE, MELT, OR EVAPORATE NUMBER CONCENTRATION
2891	! THIS FORMULATION ASSUMES 1:1 RATIO BETWEEN MASS LOSS AND
2892	! LOSS OF NUMBER CONCENTRATION
2893
2894	! IF (PCC(K).LT.0.) THEN
2895	! DUM = PCC(K)*DT/QC3D(K)
2896	! DUM = MAX(-1.,DUM)
2897	! NSUBC(K) = DUM*NC3D(K)/DT
2898	! END IF
2899
2900	IF (EPRD(K).LT.0.) THEN
2901	DUM = EPRD(K)*DT/QI3D(K)
2902	DUM = MAX(-1.,DUM)
2903	NSUBI(K) = DUM*NI3D(K)/DT
2904	END IF
2905	IF (EPRDS(K).LT.0.) THEN
2906	DUM = EPRDS(K)*DT/QNI3D(K)
2907	DUM = MAX(-1.,DUM)
2908	NSUBS(K) = DUM*NS3D(K)/DT
2909	END IF
2910	IF (PRE(K).LT.0.) THEN
2911	DUM = PRE(K)*DT/QR3D(K)
2912	DUM = MAX(-1.,DUM)
2913	NSUBR(K) = DUM*NR3D(K)/DT
2914	END IF
2915	IF (EPRDG(K).LT.0.) THEN
2916	DUM = EPRDG(K)*DT/QG3D(K)
2917	DUM = MAX(-1.,DUM)
2918	NSUBG(K) = DUM*NG3D(K)/DT
2919	END IF
2920
2921	! nsubr(k)=0.
2922	! nsubs(k)=0.
2923	! nsubg(k)=0.
2924
2925	! UPDATE TENDENCIES
2926
2927	! NC3DTEN(K) = NC3DTEN(K)+NSUBC(K)
2928	NI3DTEN(K) = NI3DTEN(K)+NSUBI(K)
2929	NS3DTEN(K) = NS3DTEN(K)+NSUBS(K)
2930	NG3DTEN(K) = NG3DTEN(K)+NSUBG(K)
2931	NR3DTEN(K) = NR3DTEN(K)+NSUBR(K)
2932
2933	END IF !!!!!! TEMPERATURE
2934
2935	! SWITCH LTRUE TO 1, SINCE HYDROMETEORS ARE PRESENT
2936	LTRUE = 1
2937
2938	200 CONTINUE
2939
2940	END DO
2941
2942	! INITIALIZE PRECIP AND SNOW RATES
2943	PRECRT = 0.
2944	SNOWRT = 0.
2945
2946	! IF THERE ARE NO HYDROMETEORS, THEN SKIP TO END OF SUBROUTINE
2947
2948	IF (LTRUE.EQ.0) GOTO 400
2949
2950	!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
2951	!.......................................................................
2952	! CALCULATE SEDIMENATION
2953	! THE NUMERICS HERE FOLLOW FROM REISNER ET AL. (1998)
2954	! FALLOUT TERMS ARE CALCULATED ON SPLIT TIME STEPS TO ENSURE NUMERICAL
2955	! STABILITY, I.E. COURANT# < 1
2956
2957	!.......................................................................
2958
2959	NSTEP = 1
2960
2961	DO K = KTS,KTE
2962
2963	DUMI(K) = QI3D(K)+QI3DTEN(K)*DT
2964	DUMQS(K) = QNI3D(K)+QNI3DTEN(K)*DT
2965	DUMR(K) = QR3D(K)+QR3DTEN(K)*DT
2966	DUMFNI(K) = NI3D(K)+NI3DTEN(K)*DT
2967	DUMFNS(K) = NS3D(K)+NS3DTEN(K)*DT
2968	DUMFNR(K) = NR3D(K)+NR3DTEN(K)*DT
2969	DUMC(K) = QC3D(K)+QC3DTEN(K)*DT
2970	DUMFNC(K) = NC3D(K)+NC3DTEN(K)*DT
2971	DUMG(K) = QG3D(K)+QG3DTEN(K)*DT
2972	DUMFNG(K) = NG3D(K)+NG3DTEN(K)*DT
2973
2974	! SWITCH FOR CONSTANT DROPLET NUMBER
2975	! IF (INUM.EQ.1) THEN
2976	DUMFNC(K) = NC3D(K)
2977	! END IF
2978
2979	! GET DUMMY LAMDA FOR SEDIMENTATION CALCULATIONS
2980
2981	! MAKE SURE NUMBER CONCENTRATIONS ARE POSITIVE
2982	DUMFNI(K) = MAX(0.,DUMFNI(K))
2983	DUMFNS(K) = MAX(0.,DUMFNS(K))
2984	DUMFNC(K) = MAX(0.,DUMFNC(K))
2985	DUMFNR(K) = MAX(0.,DUMFNR(K))
2986	DUMFNG(K) = MAX(0.,DUMFNG(K))
2987
2988	!......................................................................
2989	! CLOUD ICE
2990
2991	IF (DUMI(K).GE.QSMALL) THEN
2992	DLAMI = (CONS12DUMFNI(K)/DUMI(K))*(1./DI)
2993	DLAMI=MAX(DLAMI,LAMMINI)
2994	DLAMI=MIN(DLAMI,LAMMAXI)
2995	END IF
2996	!......................................................................
2997	! RAIN
2998
2999	IF (DUMR(K).GE.QSMALL) THEN
3000	DLAMR = (PIRHOWDUMFNR(K)/DUMR(K))**(1./3.)
3001	DLAMR=MAX(DLAMR,LAMMINR)
3002	DLAMR=MIN(DLAMR,LAMMAXR)
3003	END IF
3004	!......................................................................
3005	! CLOUD DROPLETS
3006
3007	IF (DUMC(K).GE.QSMALL) THEN
3008	DUM = PRES(K)/(287.15*T3D(K))
3009	PGAM(K)=0.0005714*(NC3D(K)/1.E6/DUM)+0.2714
3010	PGAM(K)=1./(PGAM(K)**2)-1.
3011	PGAM(K)=MAX(PGAM(K),2.)
3012	PGAM(K)=MIN(PGAM(K),10.)
3013
3014	DLAMC = (CONS26DUMFNC(K)GAMMA(PGAM(K)+4.)/(DUMC(K)GAMMA(PGAM(K)+1.)))*(1./3.)
3015	LAMMIN = (PGAM(K)+1.)/60.E-6
3016	LAMMAX = (PGAM(K)+1.)/1.E-6
3017	DLAMC=MAX(DLAMC,LAMMIN)
3018	DLAMC=MIN(DLAMC,LAMMAX)
3019	END IF
3020	!......................................................................
3021	! SNOW
3022
3023	IF (DUMQS(K).GE.QSMALL) THEN
3024	DLAMS = (CONS1DUMFNS(K)/ DUMQS(K))*(1./DS)
3025	DLAMS=MAX(DLAMS,LAMMINS)
3026	DLAMS=MIN(DLAMS,LAMMAXS)
3027	END IF
3028	!......................................................................
3029	! GRAUPEL
3030
3031	IF (DUMG(K).GE.QSMALL) THEN
3032	DLAMG = (CONS2DUMFNG(K)/ DUMG(K))*(1./DG)
3033	DLAMG=MAX(DLAMG,LAMMING)
3034	DLAMG=MIN(DLAMG,LAMMAXG)
3035	END IF
3036
3037	!......................................................................
3038	! CALCULATE NUMBER-WEIGHTED AND MASS-WEIGHTED TERMINAL FALL SPEEDS
3039
3040	! CLOUD WATER
3041
3042	IF (DUMC(K).GE.QSMALL) THEN
3043	UNC = ACN(K)GAMMA(1.+BC+PGAM(K))/ (DLAMCBCGAMMA(PGAM(K)+1.))
3044	UMC = ACN(K)GAMMA(4.+BC+PGAM(K))/ (DLAMCBCGAMMA(PGAM(K)+4.))
3045	ELSE
3046	UMC = 0.
3047	UNC = 0.
3048	END IF
3049
3050	IF (DUMI(K).GE.QSMALL) THEN
3051	UNI = AIN(K)CONS27/DLAMI*BI
3052	UMI = AIN(K)CONS28/(DLAMI*BI)
3053	ELSE
3054	UMI = 0.
3055	UNI = 0.
3056	END IF
3057
3058	IF (DUMR(K).GE.QSMALL) THEN
3059	UNR = ARN(K)CONS6/DLAMR*BR
3060	UMR = ARN(K)CONS4/(DLAMR*BR)
3061	ELSE
3062	UMR = 0.
3063	UNR = 0.
3064	END IF
3065
3066	IF (DUMQS(K).GE.QSMALL) THEN
3067	UMS = ASN(K)CONS3/(DLAMS*BS)
3068	UNS = ASN(K)CONS5/DLAMS*BS
3069	ELSE
3070	UMS = 0.
3071	UNS = 0.
3072	END IF
3073
3074	IF (DUMG(K).GE.QSMALL) THEN
3075	UMG = AGN(K)CONS7/(DLAMG*BG)
3076	UNG = AGN(K)CONS8/DLAMG*BG
3077	ELSE
3078	UMG = 0.
3079	UNG = 0.
3080	END IF
3081
3082	! SET REALISTIC LIMITS ON FALLSPEED
3083
3084	UMS=MIN(UMS,1.2)
3085	UNS=MIN(UNS,1.2)
3086	UMI=MIN(UMI,1.2)
3087	UNI=MIN(UNI,1.2)
3088	UMR=MIN(UMR,9.1)
3089	UNR=MIN(UNR,9.1)
3090	UMG=MIN(UMG,20.)
3091	UNG=MIN(UNG,20.)
3092
3093	FR(K) = UMR
3094	FI(K) = UMI
3095	FNI(K) = UNI
3096	FS(K) = UMS
3097	FNS(K) = UNS
3098	FNR(K) = UNR
3099	FC(K) = UMC
3100	FNC(K) = UNC
3101	FG(K) = UMG
3102	FNG(K) = UNG
3103
3104	! CALCULATE NUMBER OF SPLIT TIME STEPS
3105
3106	RGVM = MAX(FR(K),FI(K),FS(K),FC(K),FNI(K),FNR(K),FNS(K),FNC(K),FG(K),FNG(K))
3107	! VVT CHANGED IFIX -> INT (GENERIC FUNCTION)
3108	NSTEP = MAX(INT(RGVM*DT/DZQ(K)+1.),NSTEP)
3109
3110	! MULTIPLY VARIABLES BY RHO
3111	DUMR(k) = DUMR(k)*RHO(K)
3112	DUMI(k) = DUMI(k)*RHO(K)
3113	DUMFNI(k) = DUMFNI(K)*RHO(K)
3114	DUMQS(k) = DUMQS(K)*RHO(K)
3115	DUMFNS(k) = DUMFNS(K)*RHO(K)
3116	DUMFNR(k) = DUMFNR(K)*RHO(K)
3117	DUMC(k) = DUMC(K)*RHO(K)
3118	DUMFNC(k) = DUMFNC(K)*RHO(K)
3119	DUMG(k) = DUMG(K)*RHO(K)
3120	DUMFNG(k) = DUMFNG(K)*RHO(K)
3121
3122	END DO
3123
3124	DO N = 1,NSTEP
3125
3126	DO K = KTS,KTE
3127	FALOUTR(K) = FR(K)*DUMR(K)
3128	FALOUTI(K) = FI(K)*DUMI(K)
3129	FALOUTNI(K) = FNI(K)*DUMFNI(K)
3130	FALOUTS(K) = FS(K)*DUMQS(K)
3131	FALOUTNS(K) = FNS(K)*DUMFNS(K)
3132	FALOUTNR(K) = FNR(K)*DUMFNR(K)
3133	FALOUTC(K) = FC(K)*DUMC(K)
3134	FALOUTNC(K) = FNC(K)*DUMFNC(K)
3135	FALOUTG(K) = FG(K)*DUMG(K)
3136	FALOUTNG(K) = FNG(K)*DUMFNG(K)
3137	END DO
3138
3139	! TOP OF MODEL
3140
3141	K = KTE
3142	FALTNDR = FALOUTR(K)/DZQ(k)
3143	FALTNDI = FALOUTI(K)/DZQ(k)
3144	FALTNDNI = FALOUTNI(K)/DZQ(k)
3145	FALTNDS = FALOUTS(K)/DZQ(k)
3146	FALTNDNS = FALOUTNS(K)/DZQ(k)
3147	FALTNDNR = FALOUTNR(K)/DZQ(k)
3148	FALTNDC = FALOUTC(K)/DZQ(k)
3149	FALTNDNC = FALOUTNC(K)/DZQ(k)
3150	FALTNDG = FALOUTG(K)/DZQ(k)
3151	FALTNDNG = FALOUTNG(K)/DZQ(k)
3152	! ADD FALLOUT TERMS TO EULERIAN TENDENCIES
3153
3154	QRSTEN(K) = QRSTEN(K)-FALTNDR/NSTEP/RHO(k)
3155	QISTEN(K) = QISTEN(K)-FALTNDI/NSTEP/RHO(k)
3156	NI3DTEN(K) = NI3DTEN(K)-FALTNDNI/NSTEP/RHO(k)
3157	QNISTEN(K) = QNISTEN(K)-FALTNDS/NSTEP/RHO(k)
3158	NS3DTEN(K) = NS3DTEN(K)-FALTNDNS/NSTEP/RHO(k)
3159	NR3DTEN(K) = NR3DTEN(K)-FALTNDNR/NSTEP/RHO(k)
3160	QCSTEN(K) = QCSTEN(K)-FALTNDC/NSTEP/RHO(k)
3161	NC3DTEN(K) = NC3DTEN(K)-FALTNDNC/NSTEP/RHO(k)
3162	QGSTEN(K) = QGSTEN(K)-FALTNDG/NSTEP/RHO(k)
3163	NG3DTEN(K) = NG3DTEN(K)-FALTNDNG/NSTEP/RHO(k)
3164
3165	DUMR(K) = DUMR(K)-FALTNDR*DT/NSTEP
3166	DUMI(K) = DUMI(K)-FALTNDI*DT/NSTEP
3167	DUMFNI(K) = DUMFNI(K)-FALTNDNI*DT/NSTEP
3168	DUMQS(K) = DUMQS(K)-FALTNDS*DT/NSTEP
3169	DUMFNS(K) = DUMFNS(K)-FALTNDNS*DT/NSTEP
3170	DUMFNR(K) = DUMFNR(K)-FALTNDNR*DT/NSTEP
3171	DUMC(K) = DUMC(K)-FALTNDC*DT/NSTEP
3172	DUMFNC(K) = DUMFNC(K)-FALTNDNC*DT/NSTEP
3173	DUMG(K) = DUMG(K)-FALTNDG*DT/NSTEP
3174	DUMFNG(K) = DUMFNG(K)-FALTNDNG*DT/NSTEP
3175
3176	DO K = KTE-1,KTS,-1
3177	FALTNDR = (FALOUTR(K+1)-FALOUTR(K))/DZQ(K)
3178	FALTNDI = (FALOUTI(K+1)-FALOUTI(K))/DZQ(K)
3179	FALTNDNI = (FALOUTNI(K+1)-FALOUTNI(K))/DZQ(K)
3180	FALTNDS = (FALOUTS(K+1)-FALOUTS(K))/DZQ(K)
3181	FALTNDNS = (FALOUTNS(K+1)-FALOUTNS(K))/DZQ(K)
3182	FALTNDNR = (FALOUTNR(K+1)-FALOUTNR(K))/DZQ(K)
3183	FALTNDC = (FALOUTC(K+1)-FALOUTC(K))/DZQ(K)
3184	FALTNDNC = (FALOUTNC(K+1)-FALOUTNC(K))/DZQ(K)
3185	FALTNDG = (FALOUTG(K+1)-FALOUTG(K))/DZQ(K)
3186	FALTNDNG = (FALOUTNG(K+1)-FALOUTNG(K))/DZQ(K)
3187
3188	! ADD FALLOUT TERMS TO EULERIAN TENDENCIES
3189
3190	QRSTEN(K) = QRSTEN(K)+FALTNDR/NSTEP/RHO(k)
3191	QISTEN(K) = QISTEN(K)+FALTNDI/NSTEP/RHO(k)
3192	NI3DTEN(K) = NI3DTEN(K)+FALTNDNI/NSTEP/RHO(k)
3193	QNISTEN(K) = QNISTEN(K)+FALTNDS/NSTEP/RHO(k)
3194	NS3DTEN(K) = NS3DTEN(K)+FALTNDNS/NSTEP/RHO(k)
3195	NR3DTEN(K) = NR3DTEN(K)+FALTNDNR/NSTEP/RHO(k)
3196	QCSTEN(K) = QCSTEN(K)+FALTNDC/NSTEP/RHO(k)
3197	NC3DTEN(K) = NC3DTEN(K)+FALTNDNC/NSTEP/RHO(k)
3198	QGSTEN(K) = QGSTEN(K)+FALTNDG/NSTEP/RHO(k)
3199	NG3DTEN(K) = NG3DTEN(K)+FALTNDNG/NSTEP/RHO(k)
3200
3201	DUMR(K) = DUMR(K)+FALTNDR*DT/NSTEP
3202	DUMI(K) = DUMI(K)+FALTNDI*DT/NSTEP
3203	DUMFNI(K) = DUMFNI(K)+FALTNDNI*DT/NSTEP
3204	DUMQS(K) = DUMQS(K)+FALTNDS*DT/NSTEP
3205	DUMFNS(K) = DUMFNS(K)+FALTNDNS*DT/NSTEP
3206	DUMFNR(K) = DUMFNR(K)+FALTNDNR*DT/NSTEP
3207	DUMC(K) = DUMC(K)+FALTNDC*DT/NSTEP
3208	DUMFNC(K) = DUMFNC(K)+FALTNDNC*DT/NSTEP
3209	DUMG(K) = DUMG(K)+FALTNDG*DT/NSTEP
3210	DUMFNG(K) = DUMFNG(K)+FALTNDNG*DT/NSTEP
3211
3212	END DO
3213
3214	! GET PRECIPITATION AND SNOWFALL ACCUMULATION DURING THE TIME STEP
3215	! FACTOR OF 1000 CONVERTS FROM M TO MM, BUT DIVISION BY DENSITY
3216	! OF LIQUID WATER CANCELS THIS FACTOR OF 1000
3217
3218	PRECRT = PRECRT+(FALOUTR(KTS)+FALOUTC(KTS)+FALOUTS(KTS)+FALOUTI(KTS)+FALOUTG(KTS)) &
3219	*DT/NSTEP
3220	SNOWRT = SNOWRT+(FALOUTS(KTS)+FALOUTI(KTS)+FALOUTG(KTS))*DT/NSTEP
3221
3222	END DO
3223
3224	DO K=KTS,KTE
3225
3226	! ADD ON SEDIMENTATION TENDENCIES FOR MIXING RATIO TO REST OF TENDENCIES
3227
3228	QR3DTEN(K)=QR3DTEN(K)+QRSTEN(K)
3229	QI3DTEN(K)=QI3DTEN(K)+QISTEN(K)
3230	QC3DTEN(K)=QC3DTEN(K)+QCSTEN(K)
3231	QG3DTEN(K)=QG3DTEN(K)+QGSTEN(K)
3232	QNI3DTEN(K)=QNI3DTEN(K)+QNISTEN(K)
3233
3234	! PUT ALL CLOUD ICE IN SNOW CATEGORY IF MEAN DIAMETER EXCEEDS 2 * dcs
3235
3236	IF (QI3D(K).GE.QSMALL.AND.T3D(K).LT.273.15) THEN
3237	IF (1./LAMI(K).GE.2.*DCS) THEN
3238	QNI3DTEN(K) = QNI3DTEN(K)+QI3D(K)/DT+ QI3DTEN(K)
3239	NS3DTEN(K) = NS3DTEN(K)+NI3D(K)/DT+ NI3DTEN(K)
3240	QI3DTEN(K) = -QI3D(K)/DT
3241	NI3DTEN(K) = -NI3D(K)/DT
3242	END IF
3243	END IF
3244
3245	! hm add tendencies here, then call sizeparameter
3246	! to ensure consisitency between mixing ratio and number concentration
3247
3248	QC3D(k) = QC3D(k)+QC3DTEN(k)*DT
3249	QI3D(k) = QI3D(k)+QI3DTEN(k)*DT
3250	QNI3D(k) = QNI3D(k)+QNI3DTEN(k)*DT
3251	QR3D(k) = QR3D(k)+QR3DTEN(k)*DT
3252	! NC3D(k) = NC3D(k)+NC3DTEN(k)*DT
3253	NI3D(k) = NI3D(k)+NI3DTEN(k)*DT
3254	NS3D(k) = NS3D(k)+NS3DTEN(k)*DT
3255	NR3D(k) = NR3D(k)+NR3DTEN(k)*DT
3256
3257	IF (IGRAUP.EQ.0) THEN
3258	QG3D(k) = QG3D(k)+QG3DTEN(k)*DT
3259	NG3D(k) = NG3D(k)+NG3DTEN(k)*DT
3260	END IF
3261
3262	! ADD TEMPERATURE AND WATER VAPOR TENDENCIES FROM MICROPHYSICS
3263	T3D(K) = T3D(K)+T3DTEN(k)*DT
3264	QV3D(K) = QV3D(K)+QV3DTEN(k)*DT
3265
3266	! SATURATION VAPOR PRESSURE AND MIXING RATIO
3267
3268	EVS(K) = POLYSVP(T3D(K),0) ! PA
3269	EIS(K) = POLYSVP(T3D(K),1) ! PA
3270
3271	! MAKE SURE ICE SATURATION DOESN'T EXCEED WATER SAT. NEAR FREEZING
3272
3273	IF (EIS(K).GT.EVS(K)) EIS(K) = EVS(K)
3274
3275	QVS(K) = .622*EVS(K)/(PRES(K)-EVS(K))
3276	QVI(K) = .622*EIS(K)/(PRES(K)-EIS(K))
3277
3278	QVQVS(K) = QV3D(K)/QVS(K)
3279	QVQVSI(K) = QV3D(K)/QVI(K)
3280
3281	! AT SUBSATURATION, REMOVE SMALL AMOUNTS OF CLOUD/PRECIP WATER
3282
3283	IF (QVQVS(K).LT.0.9) THEN
3284	IF (QR3D(K).LT.1.E-6) THEN
3285	QV3D(K)=QV3D(K)+QR3D(K)
3286	T3D(K)=T3D(K)-QR3D(K)*XXLV(K)/CPM(K)
3287	QR3D(K)=0.
3288	END IF
3289	IF (QC3D(K).LT.1.E-6) THEN
3290	QV3D(K)=QV3D(K)+QC3D(K)
3291	T3D(K)=T3D(K)-QC3D(K)*XXLV(K)/CPM(K)
3292	QC3D(K)=0.
3293	END IF
3294	END IF
3295
3296	IF (QVQVSI(K).LT.0.9) THEN
3297	IF (QI3D(K).LT.1.E-6) THEN
3298	QV3D(K)=QV3D(K)+QI3D(K)
3299	T3D(K)=T3D(K)-QI3D(K)*XXLS(K)/CPM(K)
3300	QI3D(K)=0.
3301	END IF
3302	IF (QNI3D(K).LT.1.E-6) THEN
3303	QV3D(K)=QV3D(K)+QNI3D(K)
3304	T3D(K)=T3D(K)-QNI3D(K)*XXLS(K)/CPM(K)
3305	QNI3D(K)=0.
3306	END IF
3307	IF (QG3D(K).LT.1.E-6) THEN
3308	QV3D(K)=QV3D(K)+QG3D(K)
3309	T3D(K)=T3D(K)-QG3D(K)*XXLS(K)/CPM(K)
3310	QG3D(K)=0.
3311	END IF
3312	END IF
3313
3314	!..................................................................
3315	! IF MIXING RATIO < QSMALL SET MIXING RATIO AND NUMBER CONC TO ZERO
3316
3317	IF (QC3D(K).LT.QSMALL) THEN
3318	QC3D(K) = 0.
3319	NC3D(K) = 0.
3320	EFFC(K) = 0.
3321	END IF
3322	IF (QR3D(K).LT.QSMALL) THEN
3323	QR3D(K) = 0.
3324	NR3D(K) = 0.
3325	EFFR(K) = 0.
3326	END IF
3327	IF (QI3D(K).LT.QSMALL) THEN
3328	QI3D(K) = 0.
3329	NI3D(K) = 0.
3330	EFFI(K) = 0.
3331	END IF
3332	IF (QNI3D(K).LT.QSMALL) THEN
3333	QNI3D(K) = 0.
3334	NS3D(K) = 0.
3335	EFFS(K) = 0.
3336	END IF
3337	IF (QG3D(K).LT.QSMALL) THEN
3338	QG3D(K) = 0.
3339	NG3D(K) = 0.
3340	EFFG(K) = 0.
3341	END IF
3342
3343	!..................................
3344	! IF THERE IS NO CLOUD/PRECIP WATER, THEN SKIP CALCULATIONS
3345
3346	IF (QC3D(K).LT.QSMALL.AND.QI3D(K).LT.QSMALL.AND.QNI3D(K).LT.QSMALL &
3347	.AND.QR3D(K).LT.QSMALL.AND.QG3D(K).LT.QSMALL) GOTO 500
3348
3349	!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
3350	! CALCULATE INSTANTANEOUS PROCESSES
3351
3352	! ADD MELTING OF CLOUD ICE TO FORM RAIN
3353
3354	IF (QI3D(K).GE.QSMALL.AND.T3D(K).GE.273.15) THEN
3355	QR3D(K) = QR3D(K)+QI3D(K)
3356	T3D(K) = T3D(K)-QI3D(K)*XLF(K)/CPM(K)
3357	QI3D(K) = 0.
3358	NR3D(K) = NR3D(K)+NI3D(K)
3359	NI3D(K) = 0.
3360	END IF
3361
3362	! ****SENSITIVITY - NO ICE
3363	IF (ILIQ.EQ.1) GOTO 778
3364
3365	! HOMOGENEOUS FREEZING OF CLOUD WATER
3366
3367	IF (T3D(K).LE.233.15.AND.QC3D(K).GE.QSMALL) THEN
3368	QI3D(K)=QI3D(K)+QC3D(K)
3369	T3D(K)=T3D(K)+QC3D(K)*XLF(K)/CPM(K)
3370	QC3D(K)=0.
3371	NI3D(K)=NI3D(K)+NC3D(K)
3372	NC3D(K)=0.
3373	END IF
3374
3375	! HOMOGENEOUS FREEZING OF RAIN
3376
3377	IF (IGRAUP.EQ.0) THEN
3378
3379	IF (T3D(K).LE.233.15.AND.QR3D(K).GE.QSMALL) THEN
3380	QG3D(K) = QG3D(K)+QR3D(K)
3381	T3D(K) = T3D(K)+QR3D(K)*XLF(K)/CPM(K)
3382	QR3D(K) = 0.
3383	NG3D(K) = NG3D(K)+ NR3D(K)
3384	NR3D(K) = 0.
3385	END IF
3386
3387	ELSE IF (IGRAUP.EQ.1) THEN
3388
3389	IF (T3D(K).LE.233.15.AND.QR3D(K).GE.QSMALL) THEN
3390	QNI3D(K) = QNI3D(K)+QR3D(K)
3391	T3D(K) = T3D(K)+QR3D(K)*XLF(K)/CPM(K)
3392	QR3D(K) = 0.
3393	NS3D(K) = NS3D(K)+NR3D(K)
3394	NR3D(K) = 0.
3395	END IF
3396
3397	END IF
3398
3399	778 CONTINUE
3400
3401	! MAKE SURE NUMBER CONCENTRATIONS AREN'T NEGATIVE
3402
3403	NI3D(K) = MAX(0.,NI3D(K))
3404	NS3D(K) = MAX(0.,NS3D(K))
3405	NC3D(K) = MAX(0.,NC3D(K))
3406	NR3D(K) = MAX(0.,NR3D(K))
3407	NG3D(K) = MAX(0.,NG3D(K))
3408
3409	!......................................................................
3410	! CLOUD ICE
3411
3412	IF (QI3D(K).GE.QSMALL) THEN
3413	LAMI(K) = (CONS12* &
3414	NI3D(K)/QI3D(K))**(1./DI)
3415
3416	! CHECK FOR SLOPE
3417
3418	! ADJUST VARS
3419
3420	IF (LAMI(K).LT.LAMMINI) THEN
3421
3422	LAMI(K) = LAMMINI
3423
3424	N0I(K) = LAMI(K)*4QI3D(K)/CONS12
3425
3426	NI3D(K) = N0I(K)/LAMI(K)
3427	ELSE IF (LAMI(K).GT.LAMMAXI) THEN
3428	LAMI(K) = LAMMAXI
3429	N0I(K) = LAMI(K)*4QI3D(K)/CONS12
3430
3431	NI3D(K) = N0I(K)/LAMI(K)
3432	END IF
3433	END IF
3434
3435	!......................................................................
3436	! RAIN
3437
3438	IF (QR3D(K).GE.QSMALL) THEN
3439	LAMR(K) = (PIRHOWNR3D(K)/QR3D(K))**(1./3.)
3440
3441	! CHECK FOR SLOPE
3442
3443	! ADJUST VARS
3444
3445	IF (LAMR(K).LT.LAMMINR) THEN
3446
3447	LAMR(K) = LAMMINR
3448
3449	N0RR(K) = LAMR(K)*4QR3D(K)/(PI*RHOW)
3450
3451	NR3D(K) = N0RR(K)/LAMR(K)
3452	ELSE IF (LAMR(K).GT.LAMMAXR) THEN
3453	LAMR(K) = LAMMAXR
3454	N0RR(K) = LAMR(K)*4QR3D(K)/(PI*RHOW)
3455
3456	NR3D(K) = N0RR(K)/LAMR(K)
3457	END IF
3458
3459	END IF
3460
3461	!......................................................................
3462	! CLOUD DROPLETS
3463
3464	! MARTIN ET AL. (1994) FORMULA FOR PGAM
3465
3466	IF (QC3D(K).GE.QSMALL) THEN
3467
3468	DUM = PRES(K)/(287.15*T3D(K))
3469	PGAM(K)=0.0005714*(NC3D(K)/1.E6/DUM)+0.2714
3470	PGAM(K)=1./(PGAM(K)**2)-1.
3471	PGAM(K)=MAX(PGAM(K),2.)
3472	PGAM(K)=MIN(PGAM(K),10.)
3473
3474	! CALCULATE LAMC
3475
3476	LAMC(K) = (CONS26NC3D(K)GAMMA(PGAM(K)+4.)/ &
3477	(QC3D(K)GAMMA(PGAM(K)+1.)))*(1./3.)
3478
3479	! LAMMIN, 60 MICRON DIAMETER
3480	! LAMMAX, 1 MICRON
3481
3482	LAMMIN = (PGAM(K)+1.)/60.E-6
3483	LAMMAX = (PGAM(K)+1.)/1.E-6
3484
3485	IF (LAMC(K).LT.LAMMIN) THEN
3486	LAMC(K) = LAMMIN
3487	NC3D(K) = EXP(3.*LOG(LAMC(K))+LOG(QC3D(K))+ &
3488	LOG(GAMMA(PGAM(K)+1.))-LOG(GAMMA(PGAM(K)+4.)))/CONS26
3489
3490	ELSE IF (LAMC(K).GT.LAMMAX) THEN
3491	LAMC(K) = LAMMAX
3492	NC3D(K) = EXP(3.*LOG(LAMC(K))+LOG(QC3D(K))+ &
3493	LOG(GAMMA(PGAM(K)+1.))-LOG(GAMMA(PGAM(K)+4.)))/CONS26
3494
3495	END IF
3496
3497	END IF
3498
3499	!......................................................................
3500	! SNOW
3501
3502	IF (QNI3D(K).GE.QSMALL) THEN
3503	LAMS(K) = (CONS1NS3D(K)/QNI3D(K))*(1./DS)
3504
3505	! CHECK FOR SLOPE
3506
3507	! ADJUST VARS
3508
3509	IF (LAMS(K).LT.LAMMINS) THEN
3510	LAMS(K) = LAMMINS
3511	N0S(K) = LAMS(K)*4QNI3D(K)/CONS1
3512
3513	NS3D(K) = N0S(K)/LAMS(K)
3514
3515	ELSE IF (LAMS(K).GT.LAMMAXS) THEN
3516
3517	LAMS(K) = LAMMAXS
3518	N0S(K) = LAMS(K)*4QNI3D(K)/CONS1
3519	NS3D(K) = N0S(K)/LAMS(K)
3520	END IF
3521
3522	END IF
3523
3524	!......................................................................
3525	! GRAUPEL
3526
3527	IF (QG3D(K).GE.QSMALL) THEN
3528	LAMG(K) = (CONS2NG3D(K)/QG3D(K))*(1./DG)
3529
3530	! CHECK FOR SLOPE
3531
3532	! ADJUST VARS
3533
3534	IF (LAMG(K).LT.LAMMING) THEN
3535	LAMG(K) = LAMMING
3536	N0G(K) = LAMG(K)*4QG3D(K)/CONS2
3537
3538	NG3D(K) = N0G(K)/LAMG(K)
3539
3540	ELSE IF (LAMG(K).GT.LAMMAXG) THEN
3541
3542	LAMG(K) = LAMMAXG
3543	N0G(K) = LAMG(K)*4QG3D(K)/CONS2
3544
3545	NG3D(K) = N0G(K)/LAMG(K)
3546	END IF
3547
3548	END IF
3549
3550	500 CONTINUE
3551
3552	! CALCULATE EFFECTIVE RADIUS
3553
3554	IF (QI3D(K).GE.QSMALL) THEN
3555	EFFI(K) = 3./LAMI(K)/2.*1.E6
3556	ELSE
3557	EFFI(K) = 25.
3558	END IF
3559
3560	IF (QNI3D(K).GE.QSMALL) THEN
3561	EFFS(K) = 3./LAMS(K)/2.*1.E6
3562	ELSE
3563	EFFS(K) = 25.
3564	END IF
3565
3566	IF (QR3D(K).GE.QSMALL) THEN
3567	EFFR(K) = 3./LAMR(K)/2.*1.E6
3568	ELSE
3569	EFFR(K) = 25.
3570	END IF
3571
3572	IF (QC3D(K).GE.QSMALL) THEN
3573	EFFC(K) = GAMMA(PGAM(K)+4.)/ &
3574	GAMMA(PGAM(K)+3.)/LAMC(K)/2.*1.E6
3575	ELSE
3576	EFFC(K) = 25.
3577	END IF
3578
3579	IF (QG3D(K).GE.QSMALL) THEN
3580	EFFG(K) = 3./LAMG(K)/2.*1.E6
3581	ELSE
3582	EFFG(K) = 25.
3583	END IF
3584
3585	! HM ADD 1/10/06, ADD UPPER BOUND ON ICE NUMBER, THIS IS NEEDED
3586	! TO PREVENT VERY LARGE ICE NUMBER DUE TO HOMOGENEOUS FREEZING
3587	! OF DROPLETS, ESPECIALLY WHEN INUM = 1, SET MAX AT 10 CM-3
3588	NI3D(K) = MIN(NI3D(K),10.E6/RHO(K))
3589	! ADD BOUND ON DROPLET NUMBER - CANNOT EXCEED AEROSOL CONCENTRATION
3590	IF (INUM.EQ.0.AND.IACT.EQ.2) THEN
3591	NC3D(K) = MIN(NC3D(K),(NANEW1+NANEW2)/RHO(K))
3592	END IF
3593	! SWITCH FOR CONSTANT DROPLET NUMBER
3594	! IF (INUM.EQ.1) THEN
3595	! CHANGE NDCNST FROM CM-3 TO KG-1
3596	NC3D(K) = NDCNST*1.E6/RHO(K)
3597	! END IF
3598
3599	END DO !!! K LOOP
3600
3601	400 CONTINUE
3602
3603	! ALL DONE !!!!!!!!!!!
3604	RETURN
3605	END SUBROUTINE MORR_TWO_MOMENT_MICRO
3606
3607	!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
3608
3609	REAL FUNCTION POLYSVP (T,TYPE)
3610
3611	!-------------------------------------------
3612
3613	! COMPUTE SATURATION VAPOR PRESSURE
3614
3615	! POLYSVP RETURNED IN UNITS OF PA.
3616	! T IS INPUT IN UNITS OF K.
3617	! TYPE REFERS TO SATURATION WITH RESPECT TO LIQUID (0) OR ICE (1)
3618
3619	! REPLACE GOFF-GRATCH WITH FASTER FORMULATION FROM MARAT KHROUTDINOV
3620
3621	IMPLICIT NONE
3622
3623	REAL DUM
3624	REAL T
3625	INTEGER TYPE
3626	! ice
3627	real a0i,a1i,a2i,a3i,a4i,a5i,a6i,a7i,a8i
3628	data a0i,a1i,a2i,a3i,a4i,a5i,a6i,a7i,a8i /&
3629	6.11147274, 0.503160820, 0.188439774e-1, &
3630	0.420895665e-3, 0.615021634e-5,0.602588177e-7, &
3631	0.385852041e-9, 0.146898966e-11, 0.252751365e-14/
3632
3633	! liquid
3634	real a0,a1,a2,a3,a4,a5,a6,a7,a8
3635	data a0,a1,a2,a3,a4,a5,a6,a7,a8 /&
3636	6.105851, 0.4440316, 0.1430341e-1, &
3637	0.2641412e-3, 0.2995057e-5, 0.2031998e-7, &
3638	0.6936113e-10, 0.2564861e-13,-0.3704404e-15/
3639	real dt
3640
3641	! ICE
3642
3643	IF (TYPE.EQ.1) THEN
3644
3645	! POLYSVP = 10.*(-9.09718(273.16/T-1.)-3.56654* &
3646	! LOG10(273.16/T)+0.876793*(1.-T/273.16)+ &
3647	! LOG10(6.1071))*100.
3648
3649
3650	dt = max(-80.,t-273.16)
3651	polysvp = a0i + dt(a1i+dt(a2i+dt(a3i+dt(a4i+dt(a5i+dt(a6i+dt(a7i+a8idt)))))))
3652	polysvp = polysvp*100.
3653
3654	END IF
3655
3656	! LIQUID
3657
3658	IF (TYPE.EQ.0) THEN
3659
3660	dt = max(-80.,t-273.16)
3661	polysvp = a0 + dt(a1+dt(a2+dt(a3+dt(a4+dt(a5+dt(a6+dt(a7+a8dt)))))))
3662	polysvp = polysvp*100.
3663
3664	! POLYSVP = 10.*(-7.90298(373.16/T-1.)+ &
3665	! 5.02808*LOG10(373.16/T)- &
3666	! 1.3816E-7(10(11.344(1.-T/373.16))-1.)+ &
3667	! 8.1328E-3(10(-3.49149(373.16/T-1.))-1.)+ &
3668	! LOG10(1013.246))*100.
3669
3670	END IF
3671
3672
3673	END FUNCTION POLYSVP
3674
3675	!------------------------------------------------------------------------------
3676
3677	REAL FUNCTION GAMMA(X)
3678	!----------------------------------------------------------------------
3679	!
3680	! THIS ROUTINE CALCULATES THE GAMMA FUNCTION FOR A REAL ARGUMENT X.
3681	! COMPUTATION IS BASED ON AN ALGORITHM OUTLINED IN REFERENCE 1.
3682	! THE PROGRAM USES RATIONAL FUNCTIONS THAT APPROXIMATE THE GAMMA
3683	! FUNCTION TO AT LEAST 20 SIGNIFICANT DECIMAL DIGITS. COEFFICIENTS
3684	! FOR THE APPROXIMATION OVER THE INTERVAL (1,2) ARE UNPUBLISHED.
3685	! THOSE FOR THE APPROXIMATION FOR X .GE. 12 ARE FROM REFERENCE 2.
3686	! THE ACCURACY ACHIEVED DEPENDS ON THE ARITHMETIC SYSTEM, THE
3687	! COMPILER, THE INTRINSIC FUNCTIONS, AND PROPER SELECTION OF THE
3688	! MACHINE-DEPENDENT CONSTANTS.
3689	!
3690	!
3691	!*******************************************************************
3692	!*******************************************************************
3693	!
3694	! EXPLANATION OF MACHINE-DEPENDENT CONSTANTS
3695	!
3696	! BETA - RADIX FOR THE FLOATING-POINT REPRESENTATION
3697	! MAXEXP - THE SMALLEST POSITIVE POWER OF BETA THAT OVERFLOWS
3698	! XBIG - THE LARGEST ARGUMENT FOR WHICH GAMMA(X) IS REPRESENTABLE
3699	! IN THE MACHINE, I.E., THE SOLUTION TO THE EQUATION
3700	! GAMMA(XBIG) = BETA**MAXEXP
3701	! XINF - THE LARGEST MACHINE REPRESENTABLE FLOATING-POINT NUMBER;
3702	! APPROXIMATELY BETA**MAXEXP
3703	! EPS - THE SMALLEST POSITIVE FLOATING-POINT NUMBER SUCH THAT
3704	! 1.0+EPS .GT. 1.0
3705	! XMININ - THE SMALLEST POSITIVE FLOATING-POINT NUMBER SUCH THAT
3706	! 1/XMININ IS MACHINE REPRESENTABLE
3707	!
3708	! APPROXIMATE VALUES FOR SOME IMPORTANT MACHINES ARE:
3709	!
3710	! BETA MAXEXP XBIG
3711	!
3712	! CRAY-1 (S.P.) 2 8191 966.961
3713	! CYBER 180/855
3714	! UNDER NOS (S.P.) 2 1070 177.803
3715	! IEEE (IBM/XT,
3716	! SUN, ETC.) (S.P.) 2 128 35.040
3717	! IEEE (IBM/XT,
3718	! SUN, ETC.) (D.P.) 2 1024 171.624
3719	! IBM 3033 (D.P.) 16 63 57.574
3720	! VAX D-FORMAT (D.P.) 2 127 34.844
3721	! VAX G-FORMAT (D.P.) 2 1023 171.489
3722	!
3723	! XINF EPS XMININ
3724	!
3725	! CRAY-1 (S.P.) 5.45E+2465 7.11E-15 1.84E-2466
3726	! CYBER 180/855
3727	! UNDER NOS (S.P.) 1.26E+322 3.55E-15 3.14E-294
3728	! IEEE (IBM/XT,
3729	! SUN, ETC.) (S.P.) 3.40E+38 1.19E-7 1.18E-38
3730	! IEEE (IBM/XT,
3731	! SUN, ETC.) (D.P.) 1.79D+308 2.22D-16 2.23D-308
3732	! IBM 3033 (D.P.) 7.23D+75 2.22D-16 1.39D-76
3733	! VAX D-FORMAT (D.P.) 1.70D+38 1.39D-17 5.88D-39
3734	! VAX G-FORMAT (D.P.) 8.98D+307 1.11D-16 1.12D-308
3735	!
3736	!*******************************************************************
3737	!*******************************************************************
3738	!
3739	! ERROR RETURNS
3740	!
3741	! THE PROGRAM RETURNS THE VALUE XINF FOR SINGULARITIES OR
3742	! WHEN OVERFLOW WOULD OCCUR. THE COMPUTATION IS BELIEVED
3743	! TO BE FREE OF UNDERFLOW AND OVERFLOW.
3744	!
3745	!
3746	! INTRINSIC FUNCTIONS REQUIRED ARE:
3747	!
3748	! INT, DBLE, EXP, LOG, REAL, SIN
3749	!
3750	!
3751	! REFERENCES: AN OVERVIEW OF SOFTWARE DEVELOPMENT FOR SPECIAL
3752	! FUNCTIONS W. J. CODY, LECTURE NOTES IN MATHEMATICS,
3753	! 506, NUMERICAL ANALYSIS DUNDEE, 1975, G. A. WATSON
3754	! (ED.), SPRINGER VERLAG, BERLIN, 1976.
3755	!
3756	! COMPUTER APPROXIMATIONS, HART, ET. AL., WILEY AND
3757	! SONS, NEW YORK, 1968.
3758	!
3759	! LATEST MODIFICATION: OCTOBER 12, 1989
3760	!
3761	! AUTHORS: W. J. CODY AND L. STOLTZ
3762	! APPLIED MATHEMATICS DIVISION
3763	! ARGONNE NATIONAL LABORATORY
3764	! ARGONNE, IL 60439
3765	!
3766	!----------------------------------------------------------------------
3767	implicit none
3768	INTEGER I,N
3769	LOGICAL PARITY
3770	REAL &
3771	CONV,EPS,FACT,HALF,ONE,RES,SUM,TWELVE, &
3772	TWO,X,XBIG,XDEN,XINF,XMININ,XNUM,Y,Y1,YSQ,Z,ZERO
3773	REAL, DIMENSION(7) :: C
3774	REAL, DIMENSION(8) :: P
3775	REAL, DIMENSION(8) :: Q
3776	!----------------------------------------------------------------------
3777	! MATHEMATICAL CONSTANTS
3778	!----------------------------------------------------------------------
3779	DATA ONE,HALF,TWELVE,TWO,ZERO/1.0E0,0.5E0,12.0E0,2.0E0,0.0E0/
3780
3781
3782	!----------------------------------------------------------------------
3783	! MACHINE DEPENDENT PARAMETERS
3784	!----------------------------------------------------------------------
3785	DATA XBIG,XMININ,EPS/35.040E0,1.18E-38,1.19E-7/,XINF/3.4E38/
3786	!----------------------------------------------------------------------
3787	! NUMERATOR AND DENOMINATOR COEFFICIENTS FOR RATIONAL MINIMAX
3788	! APPROXIMATION OVER (1,2).
3789	!----------------------------------------------------------------------
3790	DATA P/-1.71618513886549492533811E+0,2.47656508055759199108314E+1, &
3791	-3.79804256470945635097577E+2,6.29331155312818442661052E+2, &
3792	8.66966202790413211295064E+2,-3.14512729688483675254357E+4, &
3793	-3.61444134186911729807069E+4,6.64561438202405440627855E+4/
3794	DATA Q/-3.08402300119738975254353E+1,3.15350626979604161529144E+2, &
3795	-1.01515636749021914166146E+3,-3.10777167157231109440444E+3, &
3796	2.25381184209801510330112E+4,4.75584627752788110767815E+3, &
3797	-1.34659959864969306392456E+5,-1.15132259675553483497211E+5/
3798	!----------------------------------------------------------------------
3799	! COEFFICIENTS FOR MINIMAX APPROXIMATION OVER (12, INF).
3800	!----------------------------------------------------------------------
3801	DATA C/-1.910444077728E-03,8.4171387781295E-04, &
3802	-5.952379913043012E-04,7.93650793500350248E-04, &
3803	-2.777777777777681622553E-03,8.333333333333333331554247E-02, &
3804	5.7083835261E-03/
3805	!----------------------------------------------------------------------
3806	! STATEMENT FUNCTIONS FOR CONVERSION BETWEEN INTEGER AND FLOAT
3807	!----------------------------------------------------------------------
3808	CONV(I) = REAL(I)
3809	PARITY=.FALSE.
3810	FACT=ONE
3811	N=0
3812	Y=X
3813	IF(Y.LE.ZERO)THEN
3814	!----------------------------------------------------------------------
3815	! ARGUMENT IS NEGATIVE
3816	!----------------------------------------------------------------------
3817	Y=-X
3818	Y1=AINT(Y)
3819	RES=Y-Y1
3820	IF(RES.NE.ZERO)THEN
3821	IF(Y1.NE.AINT(Y1HALF)TWO)PARITY=.TRUE.
3822	FACT=-PI/SIN(PI*RES)
3823	Y=Y+ONE
3824	ELSE
3825	RES=XINF
3826	GOTO 900
3827	ENDIF
3828	ENDIF
3829	!----------------------------------------------------------------------
3830	! ARGUMENT IS POSITIVE
3831	!----------------------------------------------------------------------
3832	IF(Y.LT.EPS)THEN
3833	!----------------------------------------------------------------------
3834	! ARGUMENT .LT. EPS
3835	!----------------------------------------------------------------------
3836	IF(Y.GE.XMININ)THEN
3837	RES=ONE/Y
3838	ELSE
3839	RES=XINF
3840	GOTO 900
3841	ENDIF
3842	ELSEIF(Y.LT.TWELVE)THEN
3843	Y1=Y
3844	IF(Y.LT.ONE)THEN
3845	!----------------------------------------------------------------------
3846	! 0.0 .LT. ARGUMENT .LT. 1.0
3847	!----------------------------------------------------------------------
3848	Z=Y
3849	Y=Y+ONE
3850	ELSE
3851	!----------------------------------------------------------------------
3852	! 1.0 .LT. ARGUMENT .LT. 12.0, REDUCE ARGUMENT IF NECESSARY
3853	!----------------------------------------------------------------------
3854	N=INT(Y)-1
3855	Y=Y-CONV(N)
3856	Z=Y-ONE
3857	ENDIF
3858	!----------------------------------------------------------------------
3859	! EVALUATE APPROXIMATION FOR 1.0 .LT. ARGUMENT .LT. 2.0
3860	!----------------------------------------------------------------------
3861	XNUM=ZERO
3862	XDEN=ONE
3863	DO I=1,8
3864	XNUM=(XNUM+P(I))*Z
3865	XDEN=XDEN*Z+Q(I)
3866	END DO
3867	RES=XNUM/XDEN+ONE
3868	IF(Y1.LT.Y)THEN
3869	!----------------------------------------------------------------------
3870	! ADJUST RESULT FOR CASE 0.0 .LT. ARGUMENT .LT. 1.0
3871	!----------------------------------------------------------------------
3872	RES=RES/Y1
3873	ELSEIF(Y1.GT.Y)THEN
3874	!----------------------------------------------------------------------
3875	! ADJUST RESULT FOR CASE 2.0 .LT. ARGUMENT .LT. 12.0
3876	!----------------------------------------------------------------------
3877	DO I=1,N
3878	RES=RES*Y
3879	Y=Y+ONE
3880	END DO
3881	ENDIF
3882	ELSE
3883	!----------------------------------------------------------------------
3884	! EVALUATE FOR ARGUMENT .GE. 12.0,
3885	!----------------------------------------------------------------------
3886	IF(Y.LE.XBIG)THEN
3887	YSQ=Y*Y
3888	SUM=C(7)
3889	DO I=1,6
3890	SUM=SUM/YSQ+C(I)
3891	END DO
3892	SUM=SUM/Y-Y+SQRTPI
3893	SUM=SUM+(Y-HALF)*LOG(Y)
3894	RES=EXP(SUM)
3895	ELSE
3896	RES=XINF
3897	GOTO 900
3898	ENDIF
3899	ENDIF
3900	!----------------------------------------------------------------------
3901	! FINAL ADJUSTMENTS AND RETURN
3902	!----------------------------------------------------------------------
3903	IF(PARITY)RES=-RES
3904	IF(FACT.NE.ONE)RES=FACT/RES
3905	900 GAMMA=RES
3906	RETURN
3907	! ---------- LAST LINE OF GAMMA ----------
3908	END FUNCTION GAMMA
3909
3910
3911	REAL FUNCTION DERF1(X)
3912	IMPLICIT NONE
3913	REAL X
3914	REAL, DIMENSION(0 : 64) :: A, B
3915	REAL W,T,Y
3916	INTEGER K,I
3917	DATA A/ &
3918	0.00000000005958930743E0, -0.00000000113739022964E0, &
3919	0.00000001466005199839E0, -0.00000016350354461960E0, &
3920	0.00000164610044809620E0, -0.00001492559551950604E0, &
3921	0.00012055331122299265E0, -0.00085483269811296660E0, &
3922	0.00522397762482322257E0, -0.02686617064507733420E0, &
3923	0.11283791670954881569E0, -0.37612638903183748117E0, &
3924	1.12837916709551257377E0, &
3925	0.00000000002372510631E0, -0.00000000045493253732E0, &
3926	0.00000000590362766598E0, -0.00000006642090827576E0, &
3927	0.00000067595634268133E0, -0.00000621188515924000E0, &
3928	0.00005103883009709690E0, -0.00037015410692956173E0, &
3929	0.00233307631218880978E0, -0.01254988477182192210E0, &
3930	0.05657061146827041994E0, -0.21379664776456006580E0, &
3931	0.84270079294971486929E0, &
3932	0.00000000000949905026E0, -0.00000000018310229805E0, &
3933	0.00000000239463074000E0, -0.00000002721444369609E0, &
3934	0.00000028045522331686E0, -0.00000261830022482897E0, &
3935	0.00002195455056768781E0, -0.00016358986921372656E0, &
3936	0.00107052153564110318E0, -0.00608284718113590151E0, &
3937	0.02986978465246258244E0, -0.13055593046562267625E0, &
3938	0.67493323603965504676E0, &
3939	0.00000000000382722073E0, -0.00000000007421598602E0, &
3940	0.00000000097930574080E0, -0.00000001126008898854E0, &
3941	0.00000011775134830784E0, -0.00000111992758382650E0, &
3942	0.00000962023443095201E0, -0.00007404402135070773E0, &
3943	0.00050689993654144881E0, -0.00307553051439272889E0, &
3944	0.01668977892553165586E0, -0.08548534594781312114E0, &
3945	0.56909076642393639985E0, &
3946	0.00000000000155296588E0, -0.00000000003032205868E0, &
3947	0.00000000040424830707E0, -0.00000000471135111493E0, &
3948	0.00000005011915876293E0, -0.00000048722516178974E0, &
3949	0.00000430683284629395E0, -0.00003445026145385764E0, &
3950	0.00024879276133931664E0, -0.00162940941748079288E0, &
3951	0.00988786373932350462E0, -0.05962426839442303805E0, &
3952	0.49766113250947636708E0 /
3953	DATA (B(I), I = 0, 12) / &
3954	-0.00000000029734388465E0, 0.00000000269776334046E0, &
3955	-0.00000000640788827665E0, -0.00000001667820132100E0, &
3956	-0.00000021854388148686E0, 0.00000266246030457984E0, &
3957	0.00001612722157047886E0, -0.00025616361025506629E0, &
3958	0.00015380842432375365E0, 0.00815533022524927908E0, &
3959	-0.01402283663896319337E0, -0.19746892495383021487E0, &
3960	0.71511720328842845913E0 /
3961	DATA (B(I), I = 13, 25) / &
3962	-0.00000000001951073787E0, -0.00000000032302692214E0, &
3963	0.00000000522461866919E0, 0.00000000342940918551E0, &
3964	-0.00000035772874310272E0, 0.00000019999935792654E0, &
3965	0.00002687044575042908E0, -0.00011843240273775776E0, &
3966	-0.00080991728956032271E0, 0.00661062970502241174E0, &
3967	0.00909530922354827295E0, -0.20160072778491013140E0, &
3968	0.51169696718727644908E0 /
3969	DATA (B(I), I = 26, 38) / &
3970	0.00000000003147682272E0, -0.00000000048465972408E0, &
3971	0.00000000063675740242E0, 0.00000003377623323271E0, &
3972	-0.00000015451139637086E0, -0.00000203340624738438E0, &
3973	0.00001947204525295057E0, 0.00002854147231653228E0, &
3974	-0.00101565063152200272E0, 0.00271187003520095655E0, &
3975	0.02328095035422810727E0, -0.16725021123116877197E0, &
3976	0.32490054966649436974E0 /
3977	DATA (B(I), I = 39, 51) / &
3978	0.00000000002319363370E0, -0.00000000006303206648E0, &
3979	-0.00000000264888267434E0, 0.00000002050708040581E0, &
3980	0.00000011371857327578E0, -0.00000211211337219663E0, &
3981	0.00000368797328322935E0, 0.00009823686253424796E0, &
3982	-0.00065860243990455368E0, -0.00075285814895230877E0, &
3983	0.02585434424202960464E0, -0.11637092784486193258E0, &
3984	0.18267336775296612024E0 /
3985	DATA (B(I), I = 52, 64) / &
3986	-0.00000000000367789363E0, 0.00000000020876046746E0, &
3987	-0.00000000193319027226E0, -0.00000000435953392472E0, &
3988	0.00000018006992266137E0, -0.00000078441223763969E0, &
3989	-0.00000675407647949153E0, 0.00008428418334440096E0, &
3990	-0.00017604388937031815E0, -0.00239729611435071610E0, &
3991	0.02064129023876022970E0, -0.06905562880005864105E0, &
3992	0.09084526782065478489E0 /
3993	W = ABS(X)
3994	IF (W .LT. 2.2D0) THEN
3995	T = W * W
3996	K = INT(T)
3997	T = T - K
3998	K = K * 13
3999	Y = ((((((((((((A(K) * T + A(K + 1)) * T + &
4000	A(K + 2)) * T + A(K + 3)) * T + A(K + 4)) * T + &
4001	A(K + 5)) * T + A(K + 6)) * T + A(K + 7)) * T + &
4002	A(K + 8)) * T + A(K + 9)) * T + A(K + 10)) * T + &
4003	A(K + 11)) * T + A(K + 12)) * W
4004	ELSE IF (W .LT. 6.9D0) THEN
4005	K = INT(W)
4006	T = W - K
4007	K = 13 * (K - 2)
4008	Y = (((((((((((B(K) * T + B(K + 1)) * T + &
4009	B(K + 2)) * T + B(K + 3)) * T + B(K + 4)) * T + &
4010	B(K + 5)) * T + B(K + 6)) * T + B(K + 7)) * T + &
4011	B(K + 8)) * T + B(K + 9)) * T + B(K + 10)) * T + &
4012	B(K + 11)) * T + B(K + 12)
4013	Y = Y * Y
4014	Y = Y * Y
4015	Y = Y * Y
4016	Y = 1 - Y * Y
4017	ELSE
4018	Y = 1
4019	END IF
4020	IF (X .LT. 0) Y = -Y
4021	DERF1 = Y
4022	END FUNCTION DERF1
4023
4024	!+---+-----------------------------------------------------------------+
4025
4026	END MODULE module_mp_morr_two_moment

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format