Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

module_mp_milbrandt2mom.F @ 1

Last change on this file since 1 was 1, checked in by lfita, 10 years ago

-- --- Opening of the WRF+LMDZ coupling repository --- -- -

WRF: version v3.3
LMDZ: version v1818

More details in:

File size: 184.1 KB

Line
1	!_______________________________________________________________________________!
2	! !
3	! This module contains the microphysics sub-driver for the 2-moment version of !
4	! the Milbrandt-Yau (2005, JAS) microphysics scheme, along with all associated !
5	! subprograms. The main subroutine, 'mp_milbrandt2mom_main', is essentially !
6	! directly from the RPN-CMC physics library of the Canadian GEM model. It is !
7	! called by the wrapper 'mp_milbrandt2mom_driver' which makes the necessary !
8	! adjustments to the calling parameters for the interface to WRF. !
9	! !
10	! For questions, bug reports, etc. pertaining to the scheme, or to request !
11	! updates to the code (before the next offical WRF release) please contact !
12	! Jason Milbrandt (Environment Canada) at jason.milbrandt@ec.gc.ca !
13	! !
14	! Last modified: 2011-03-02 !
15	!_______________________________________________________________________________!
16
17	module my_fncs_mod
18
19	!==============================================================================!
20	! The following functions are used by the schemes in the multimoment package. !
21	! !
22	! Package version: 2.19.0 (internal bookkeeping) !
23	! Last modified : 2009-04-27 !
24	!==============================================================================!
25
26	implicit none
27
28	private
29	public :: NccnFNC,SxFNC,gamma,gammaDP,gser,gammln,gammp,cfg,gamminc
30
31	contains
32
33	!==============================================================================!
34
35	REAL FUNCTION NccnFNC(Win,Tin,Pin,CCNtype)
36
37	!---------------------------------------------------------------------------!
38	! This function returns number concentration (activated aerosols) as a
39	! function of w,T,p, based on polynomial approximations of detailed
40	! approach using a hypergeometric function, following Cohard and Pinty (2000a).
41	!---------------------------------------------------------------------------!
42
43	IMPLICIT NONE
44
45	! PASSING PARAMETERS:
46	real, intent(in) :: Win, Tin, Pin
47	integer, intent(in) :: CCNtype
48
49	! LOCAL PARAMETERS:
50	real :: T,p,x,y,a,b,c,d,e,f,g,h,T2,T3,T4,x2,x3,x4,p2
51
52	x= log10(Win100.); x2= xx; x3= x2x; x4= x2x2
53	T= Tin - 273.15; T2= TT; T3= T2T; T4= T2*T2
54	p= Pin0.01; p2= pp
55
56	if (CCNtype==1) then !Maritime
57
58	a= 1.47e-9T4 -6.944e-8T3 -9.933e-7T2 +2.7278e-4T -6.6853e-4
59	b=-1.41e-8T4 +6.662e-7T3 +4.483e-6T2 -2.0479e-3T +4.0823e-2
60	c= 5.12e-8T4 -2.375e-6T3 +4.268e-6T2 +3.9681e-3T -3.2356e-1
61	d=-8.25e-8T4 +3.629e-6T3 -4.044e-5T2 +2.1846e-3T +9.1227e-1
62	e= 5.02e-8T4 -1.973e-6T3 +3.944e-5T2 -9.0734e-3T +1.1256e0
63	f= -1.424e-6p2 +3.631e-3p -1.986
64	g= -0.0212x4 +0.1765x3 -0.3770x2 -0.2200x +1.0081
65	h= 2.47e-6T3 -3.654e-5T2 +2.3327e-3*T +0.1938
66	y= ax4 + bx3 + cx2 + dx + e + fgh
67	NccnFNC= 10.*min(2.,max(0.,y)) 1.e6 ![m-3]
68
69	else if (CCNtype==2) then !Continental
70
71	a= 0.
72	b= 0.
73	c=-2.112e-9T4 +3.9836e-8T3 +2.3703e-6T2 -1.4542e-4T -0.0698
74	d=-4.210e-8T4 +5.5745e-7T3 +1.8460e-5T2 +9.6078e-4T +0.7120
75	e= 1.434e-7T4 -1.6455e-6T3 -4.3334e-5T2 -7.6720e-3T +1.0056
76	f= 1.340e-6p2 -3.5114e-3p +1.9453
77	g= 4.226e-3x4 -5.6012e-3x3 -8.7846e-2x2 +2.7435e-2x +0.9932
78	h= 5.811e-9T4 +1.5589e-7T3 -3.8623e-5T2 +1.4471e-3T +0.1496
79	y= ax4 +bx3 +cx2 + dx + e + (fgh)
80	NccnFNC= 10.*max(0.,y) 1.e6
81
82	else
83
84	print, ' STOPPED in MODULE ### NccnFNC * '
85	print*, ' Parameter CCNtype incorrectly specified'
86	stop
87
88	endif
89
90	END FUNCTION NccnFNC
91	!======================================================================!
92
93	real FUNCTION SxFNC(Win,Tin,Pin,Qsw,Qsi,CCNtype,WRT)
94
95	!---------------------------------------------------------------------------!
96	! This function returns the peak supersaturation achieved during ascent with
97	! activation of CCN aerosols as a function of w,T,p, based on polynomial
98	! approximations of detailed approach using a hypergeometric function,
99	! following Cohard and Pinty (2000a).
100	!---------------------------------------------------------------------------!
101
102	IMPLICIT NONE
103
104	! PASSING PARAMETERS:
105	integer, intent(IN) :: WRT
106	integer, intent(IN) :: CCNtype
107	real, intent(IN) :: Win, Tin, Pin, Qsw, Qsi
108
109	! LOCAL PARAMETERS:
110	real :: Si,Sw,Qv,T,p,x,a,b,c,d,f,g,h,Pcorr,T2corr,T2,T3,T4,x2,x3,x4,p2
111	real, parameter :: TRPL= 273.15
112
113	x= log10(max(Win,1.e-20)100.); x2= xx; x3= x2x; x4= x2x2
114	T= Tin; T2= TT; T3= T2T; T4= T2*T2
115	p= Pin0.01; p2= pp
116
117	if (CCNtype==1) then !Maritime
118
119	a= -5.109e-7T4 -3.996e-5T3 -1.066e-3T2 -1.273e-2T +0.0659
120	b= 2.014e-6T4 +1.583e-4T3 +4.356e-3T2 +4.943e-2T -0.1538
121	c= -2.037e-6T4 -1.625e-4T3 -4.541e-3T2 -5.118e-2T +0.1428
122	d= 3.812e-7T4 +3.065e-5T3 +8.795e-4T2 +9.440e-3T +6.14e-3
123	f= -2.012e-6p2 + 4.1913e-3p - 1.785e0
124	g= 2.832e-1x3 -5.6990e-1x2 +5.1105e-1*x -4.1747e-4
125	h= 1.173e-6T3 +3.2174e-5T2 -6.8832e-4*T +6.7888e-2
126	Pcorr= fgh
127	T2corr= 0.9581-4.449e-3T-2.016e-4T2-3.307e-6T3-1.725e-8T4
128
129	else if (CCNtype==2) then !Continental [computed for -35<T<-5C]
130
131	a= 3.80e-5T2 +1.65e-4T +9.88e-2
132	b= -7.38e-5T2 -2.53e-3T -3.23e-1
133	c= 8.39e-5T2 +3.96e-3T +3.50e-1
134	d= -1.88e-6T2 -1.33e-3T -3.73e-2
135	f= -1.9761e-6p2 + 4.1473e-3p - 1.771e0
136	g= 0.1539x4 -0.5575x3 +0.9262x2 -0.3498x -0.1293
137	h=-8.035e-9T4+3.162e-7T3+1.029e-5T2-5.931e-4T+5.62e-2
138	Pcorr= fgh
139	T2corr= 0.98888-5.0525e-4T-1.7598e-5T2-8.3308e-8*T3
140
141	else
142
143	print, ' STOPPED in MODULE ### SxFNC * '
144	print*, ' Parameter CCNtype incorrectly specified'
145	stop
146
147	endif
148
149	Sw= (ax3 + bx2 +c*x + d) + Pcorr
150	Sw= 1. + 0.01*Sw
151	Qv= Qsw*Sw
152	Si= Qv/Qsi
153	Si= Si*T2corr
154	if (WRT.eq.1) then
155	SxFNC= Sw
156	else
157	SxFNC= Si
158	endif
159	if (Win.le.0.) SxFNC= 1.
160
161	END function SxFNC
162	!======================================================================!
163
164	real FUNCTION gamma(xx)
165
166	! Modified from "Numerical Recipes"
167
168	IMPLICIT NONE
169
170	! PASSING PARAMETERS:
171	real, intent(IN) :: xx
172
173	! LOCAL PARAMETERS:
174	integer :: j
175	real*8 :: ser,stp,tmp,x,y,cof(6),gammadp
176
177
178	SAVE cof,stp
179	DATA cof,stp/76.18009172947146d0,-86.50532032941677d0, &
180	24.01409824083091d0,-1.231739572450155d0,.1208650973866179d-2, &
181	-.5395239384953d-5,2.5066282746310005d0/
182	x=dble(xx)
183	y=x
184	tmp=x+5.5d0
185	tmp=(x+0.5d0)*log(tmp)-tmp
186	ser=1.000000000190015d0
187	! do j=1,6 !original
188	do j=1,4
189	!!do j=1,3 !gives result to within ~ 3 %
190	y=y+1.d0
191	ser=ser+cof(j)/y
192	enddo
193	gammadp=tmp+log(stp*ser/x)
194	gammadp= exp(gammadp)
195
196	#if (DWORDSIZE == 8 && RWORDSIZE == 8)
197	gamma = gammadp
198	#elif (DWORDSIZE == 8 && RWORDSIZE == 4)
199	gamma = sngl(gammadp)
200	#else
201	This is a temporary hack assuming double precision is 8 bytes.
202	#endif
203
204	END FUNCTION gamma
205	!======================================================================!
206	! ! !
207	! ! ! -- USED BY DIAGNOSTIC-ALPHA DOUBLE-MOMENT (SINGLE-PRECISION) VERSION --
208	! ! ! FOR FUTURE VERSIONS OF M-Y PACKAGE WITH, THIS S/R CAN BE USED
209	! ! !
210	! ! ! real FUNCTION diagAlpha(Dm,x)
211	! ! !
212	! ! ! IMPLICIT NONE
213	! ! !
214	! ! ! integer :: x
215	! ! ! real :: Dm
216	! ! ! real, dimension(5) :: c1,c2,c3,c4
217	! ! ! real, parameter :: pi = 3.14159265
218	! ! ! real, parameter :: alphaMAX= 80.e0
219	! ! ! data c1 /19.0, 12.0, 4.5, 5.5, 3.7/
220	! ! ! data c2 / 0.6, 0.7, 0.5, 0.7, 0.3/
221	! ! ! data c3 / 1.8, 1.7, 5.0, 4.5, 9.0/
222	! ! ! data c4 /17.0, 11.0, 5.5, 8.5, 6.5/
223	! ! ! diagAlpha= c1(x)tanh(c2(x)(1.e3*Dm-c3(x)))+c4(x)
224	! ! ! if (x==5.and.Dm>0.008) diagAlpha= 1.e3*Dm-2.6
225	! ! ! diagAlpha= min(diagAlpha, alphaMAX)
226	! ! !
227	! ! ! END function diagAlpha
228	! ! !
229	! ! ! !======================================================================!
230	! ! !
231	! ! ! -- USED BY DIAGNOSTIC-ALPHA DOUBLE-MOMENT (SINGLE-PRECISION) VERSION --
232	! ! ! FOR FUTURE VERSIONS OF M-Y PACKAGE WITH, THIS S/R CAN BE USED
233	! ! !
234	! ! ! real FUNCTION solveAlpha(Q,N,Z,Cx,rho)
235	! ! !
236	! ! ! IMPLICIT NONE
237	! ! !
238	! ! ! ! PASSING PARAMETERS:
239	! ! ! real, intent(IN) :: Q, N, Z, Cx, rho
240	! ! !
241	! ! ! ! LOCAL PARAMETERS:
242	! ! ! real :: a,g,a1,g1,g2,tmp1
243	! ! ! integer :: i
244	! ! ! real, parameter :: alphaMax= 40.
245	! ! ! real, parameter :: epsQ = 1.e-14
246	! ! ! real, parameter :: epsN = 1.e-3
247	! ! ! real, parameter :: epsZ = 1.e-32
248	! ! !
249	! ! ! ! Q mass mixing ratio
250	! ! ! ! N total concentration
251	! ! ! ! Z reflectivity
252	! ! ! ! Cx (pi/6)*RHOx
253	! ! ! ! rho air density
254	! ! ! ! a alpha (returned as solveAlpha)
255	! ! ! ! g function g(a)= [(6+a)(5+a)(4+a)]/[(3+a)(2+a)(1+a)],
256	! ! ! ! where g = (Cx/(rhoQ))2.(Z*N)
257	! ! !
258	! ! !
259	! ! ! if (Q==0. .or. N==0. .or. Z==0. .or. Cx==0. .or. rho==0.) then
260	! ! ! ! For testing/debugging only; this module should never be called
261	! ! ! ! if the above condition is true.
262	! ! ! print,' STOPPED in MODULE ### solveAlpha * '
263	! ! ! print,'** : ',Q,N,Z,Cx*1.9099,rho
264	! ! ! stop
265	! ! ! endif
266	! ! !
267	! ! ! IF (Q>epsQ .and. N>epsN .and. Z>epsZ ) THEN
268	! ! !
269	! ! ! tmp1= Cx/(rho*Q)
270	! ! ! g = tmp1Ztmp1N ! g = (ZN)[Cx / (rhoQ)]^2
271	! ! !
272	! ! ! !Note: The above order avoids OVERFLOW, since tmp1*tmp1 is very large
273	! ! !
274	! ! ! !----------------------------------------------------------!
275	! ! ! ! !Solve alpha numerically: (brute-force; for testing only)
276	! ! ! ! a= 0.
277	! ! ! ! g2= 999.
278	! ! ! ! do i=0,4000
279	! ! ! ! a1= i*0.01
280	! ! ! ! g1= (6.+a1)(5.+a1)(4.+a1)/((3.+a1)(2.+a1)(1.+a1))
281	! ! ! ! if(abs(g-g1)<abs(g-g2)) then
282	! ! ! ! a = a1
283	! ! ! ! g2= g1
284	! ! ! ! endif
285	! ! ! ! enddo
286	! ! ! !----------------------------------------------------------!
287	! ! !
288	! ! ! !Piecewise-polynomial approximation of g(a) to solve for a: [2004-11-29]
289	! ! ! if (g>=20.) then
290	! ! ! a= 0.
291	! ! ! else
292	! ! ! g2= g*g
293	! ! ! if (g<20. .and.g>=13.31) a= 3.3638e-3g2 - 1.7152e-1g + 2.0857e+0
294	! ! ! if (g<13.31.and.g>=7.123) a= 1.5900e-2g2 - 4.8202e-1g + 4.0108e+0
295	! ! ! if (g<7.123.and.g>=4.200) a= 1.0730e-1g2 - 1.7481e+0g + 8.4246e+0
296	! ! ! if (g<4.200.and.g>=2.946) a= 5.9070e-1g2 - 5.7918e+0g + 1.6919e+1
297	! ! ! if (g<2.946.and.g>=1.793) a= 4.3966e+0g2 - 2.6659e+1g + 4.5477e+1
298	! ! ! if (g<1.793.and.g>=1.405) a= 4.7552e+1g2 - 1.7958e+2g + 1.8126e+2
299	! ! ! if (g<1.405.and.g>=1.230) a= 3.0889e+2g2 - 9.0854e+2g + 6.8995e+2
300	! ! ! if (g<1.230) a= alphaMax
301	! ! ! endif
302	! ! !
303	! ! ! solveAlpha= max(0.,min(a,alphaMax))
304	! ! !
305	! ! ! ELSE
306	! ! !
307	! ! ! solveAlpha= 0.
308	! ! !
309	! ! ! ENDIF
310	! ! !
311	! ! ! END FUNCTION solveAlpha
312	!======================================================================!
313
314	FUNCTION gammaDP(xx)
315
316	! Modified from "Numerical Recipes"
317
318	IMPLICIT NONE
319
320	! PASSING PARAMETERS:
321	DOUBLE PRECISION, INTENT(IN) :: xx
322
323	! LOCAL PARAMETERS:
324	DOUBLE PRECISION :: gammaDP
325	INTEGER :: j
326	DOUBLE PRECISION :: ser,stp,tmp,x,y,cof(6)
327
328
329	SAVE cof,stp
330	DATA cof,stp/76.18009172947146d0,-86.50532032941677d0, &
331	24.01409824083091d0,-1.231739572450155d0,.1208650973866179d-2, &
332	-.5395239384953d-5,2.5066282746310005d0/
333	x=xx
334	y=x
335	tmp=x+5.5d0
336	tmp=(x+0.5d0)*log(tmp)-tmp
337	ser=1.000000000190015d0
338	! do j=1,6 !original
339	do j=1,4
340	!!do j=1,3 !gives result to within ~ 3 %
341	y=y+1.d0
342	ser=ser+cof(j)/y
343	enddo
344	gammaDP=tmp+log(stp*ser/x)
345	gammaDP= exp(gammaDP)
346
347	END FUNCTION gammaDP
348	!======================================================================!
349
350	SUBROUTINE gser(gamser,a,x,gln)
351
352	! USES gammln
353
354	! Returns the incomplete gamma function P(a,x) evaluated by its series
355	! representation as gamser. Also returns GAMMA(a) as gln.
356
357	implicit none
358
359	integer :: itmax
360	real :: a,gamser,gln,x,eps
361	parameter (itmax=100, eps=3.e-7)
362	integer :: n
363	real :: ap,de1,summ
364
365	gln=gammln(a)
366	if(x.le.0.)then
367	if(x.lt.0.)pause 'x <0 in gser'
368	gamser=0.
369	return
370	endif
371	ap=a
372	summ=1./a
373	de1=summ
374	do n=1,itmax
375	ap=ap+1.
376	de1=de1*x/ap
377	summ=summ+de1
378	if(abs(de1).lt.abs(summ)*eps) goto 1
379	enddo
380	pause 'a too large, itmax too small in gser'
381	1 gamser=summexp(-x+alog(x)-gln)
382	return
383
384	END SUBROUTINE gser
385	!======================================================================!
386
387	real FUNCTION gammln(xx)
388
389	! Returns value of ln(GAMMA(xx)) for xx>0
390	! (modified from "Numerical Recipes")
391
392	IMPLICIT NONE
393
394	! PASSING PARAMETERS:
395	real, intent(IN) :: xx
396
397	! LOCAL PARAMETERS:
398	integer :: j
399	real*8 :: ser,stp,tmp,x,y,cof(6)
400
401	SAVE cof,stp
402	DATA cof,stp/76.18009172947146d0,-86.50532032941677d0, &
403	24.01409824083091d0,-1.231739572450155d0,.1208650973866179d-2, &
404	-.5395239384953d-5,2.5066282746310005d0/
405	x=dble(xx)
406	y=x
407	tmp=x+5.5d0
408	tmp=(x+0.5d0)*log(tmp)-tmp
409	ser=1.000000000190015d0
410	do j=1,6 !original
411	! do j=1,4
412	y=y+1.d0
413	ser=ser+cof(j)/y
414	enddo
415	#if (DWORDSIZE == 8 && RWORDSIZE == 8)
416	gammln= tmp+log(stp*ser/x)
417	#elif (DWORDSIZE == 8 && RWORDSIZE == 4)
418	gammln= sngl( tmp+log(stp*ser/x) )
419	#else
420	This is a temporary hack assuming double precision is 8 bytes.
421	#endif
422
423	END FUNCTION gammln
424	!======================================================================!
425
426	real FUNCTION gammp(a,x)
427
428	! USES gcf,gser
429
430	! Returns the incomplete gamma function P(a,x)
431
432	implicit none
433
434	real :: a,x,gammcf,gamser,gln
435
436	if(x.lt.0..or.a.le.0.) pause 'bad arguments in gammq'
437	if(x.lt.a+1.)then
438	call gser(gamser,a,x,gln)
439	gammp=gamser
440	else
441	call cfg(gammcf,a,x,gln)
442	gammp=1.-gammcf
443	endif
444	return
445
446	END FUNCTION gammp
447	!======================================================================!
448
449	SUBROUTINE cfg(gammcf,a,x,gln)
450
451	! USES gammln
452
453	! Returns the incomplete gamma function (Q(a,x) evaluated by tis continued fraction
454	! representation as gammcf. Also returns ln(GAMMA(a)) as gln. ITMAX is the maximum
455	! allowed number of iterations; EPS is the relative accuracy; FPMIN is a number near
456	! the smallest representable floating-point number.
457
458	implicit none
459
460	integer :: i,itmax
461	real :: a,gammcf,gln,x,eps,fpmin
462	real :: an,b,c,d,de1,h
463	parameter (itmax=100,eps=3.e-7)
464
465	gln=gammln(a)
466	b=x+1.-a
467	c=1./fpmin
468	d=1./b
469	h=d
470	do i= 1,itmax
471	an=-i*(i-a)
472	b=b+2.
473	d=an*d+b
474	if(abs(d).lt.fpmin)d=fpmin
475	c=b+an/c
476	if(abs(c).lt.fpmin) c=fpmin
477	d=1./d
478	de1=d*c
479	h=h*de1
480	if(abs(de1-1.).lt.eps) goto 1
481	enddo
482	pause 'a too large, itmax too small in gcf'
483	1 gammcf=exp(-x+alog(x)-gln)h
484	return
485
486	END SUBROUTINE cfg
487	!======================================================================!
488
489	real FUNCTION gamminc(p,xmax)
490
491	! USES gammp, gammln
492	! Returns incomplete gamma function, gamma(p,xmax)= P(p,xmax)*GAMMA(p)
493	real :: p,xmax
494	gamminc= gammp(p,xmax)*exp(gammln(p))
495
496	end FUNCTION gamminc
497
498	!======================================================================!
499	! real function x_tothe_y(x,y)
500	!
501	! implicit none
502	! real, intent(in) :: x,y
503	! x_tothe_y= exp(y*log(x))
504	!
505	! end function x_tothe_y
506	!======================================================================!
507
508	end module my_fncs_mod
509
510	!________________________________________________________________________________________!
511
512	module my_sedi_mod
513
514	!================================================================================!
515	! The following subroutines are used by the schemes in the multimoment package. !
516	! !
517	! Package version: 2.19.0 (internal bookkeeping) !
518	! Last modified : 2011-01-07 !
519	!================================================================================!
520
521	implicit none
522
523	private
524	public :: SEDI_main_1b,SEDI_main_2,countColumns
525
526	contains
527
528	!=====================================================================================!
529	SUBROUTINE SEDI_main_2(QX,NX,cat,Q,T,DE,iDE,gamfact,epsQ,epsN,afx,bfx,cmx,dmx, &
530	ckQx1,ckQx2,ckQx4,LXP,ni,nk,VxMax,DxMax,dt,DZ,massFlux, &
531	ktop_sedi,GRAV,massFlux3D)
532
533	!-------------------------------------------------------------------------------------!
534	! DOUBLE-MOMENT version of sedimentation subroutine for categories whose
535	! fall velocity equation is V(D) = gamfact * afx * D^bfx
536	!-------------------------------------------------------------------------------------!
537
538	! Passing parameters:
539	!
540	! VAR Description
541	! --- ------------
542	! QX mass mixing ratio of category x
543	! NX number concentration of category x
544	! cat: hydrometeor category:
545	! 1 rain
546	! 2 ice
547	! 3 snow
548	! 4 graupel
549	! 5 hail
550	!-------------------------------------------------------------------------------------!
551
552	use my_fncs_mod
553
554	implicit none
555
556	! PASSING PARAMETERS:
557	real, dimension(:,:), intent(inout) :: QX,NX,Q,T
558	real, dimension(:), intent(out) :: massFlux
559	real, optional, dimension(:,:), intent(out) :: massFlux3D
560	real, dimension(:,:), intent(in) :: DE,iDE,DZ
561	real, intent(in) :: epsQ,epsN,VxMax,LXP,afx,bfx,cmx,dmx,ckQx1,ckQx2,ckQx4,DxMax,dt,GRAV
562	integer, dimension(:), intent(in) :: ktop_sedi
563	integer, intent(in) :: ni,nk,cat
564
565	! LOCAL PARAMETERS:
566	logical :: slabHASmass,locallim,QxPresent
567	integer :: nnn,a,i,k,counter,l,km1,kp1,ks,kw,idzmin
568	integer, dimension(nk) :: flim_Q,flim_N
569	integer, dimension(ni) :: activeColumn,npassx,ke
570	real :: VqMax,VnMax,iLAMx,iLAMxB0,tmp1,tmp2,tmp3,Dx,iDxMax,icmx, &
571	VincFact,ratio_Vn2Vq,zmax_Q,zmax_N,tempo,idmx,Nos_Thompson, &
572	No_s,iLAMs
573	real, dimension(ni,nk) :: VVQ,VVN,RHOQX,gamfact
574	real, dimension(ni) :: dzMIN,dtx,VxMaxx
575	real, dimension(nk) :: vp_Q,vp_N,zt_Q,zt_N,zb_Q,zb_N,dzi,Q_star,N_star
576	real, dimension(0:nk) :: zz
577	real, parameter :: epsilon = 1.e-2
578	real, parameter :: thrd = 1./3.
579	real, parameter :: sxth = 1./6.
580	real, parameter :: CoMAX = 2.0
581
582	!-------------------------------------------------------------------------------------!
583
584	massFlux = 0.
585
586	!Factor to estimate increased V from size-sorting:
587	! - this factor should be higher for categories with more time-splitting, since Vmax
588	! increases after each sedimentation split step (to be tuned)
589	VincFact = 1.
590	if (present(massFlux3D)) massFlux3D= 0. !(for use in MAIN for diagnostics)
591
592	!Determine for which slabs and columns sedimentation should be computes:
593	call countColumns(QX,ni,nk,epsQ,counter,activeColumn,ktop_sedi)
594
595	ratio_Vn2Vq= ckQx2/ckQx1
596	iDxMax= 1./DxMax
597	icmx = 1./cmx
598	idmx = 1./dmx
599	ks = nk
600	ke = ktop_sedi !(i-array) - formerly ke=1; now depends on max. level with hydrometeor
601	kw = -1 !direction of vertical leveling; -1 implies nk is bottom
602
603	VVQ = 0.
604	VVN = 0.
605	VqMax= 0.
606	VnMax= 0.
607
608	DO a= 1,counter
609	i= activeColumn(a)
610
611	VVQ(i,:) = 0.
612	do k= ktop_sedi(i),nk !formerly do k= 1,nk
613	QxPresent = (QX(i,k)>epsQ .and. NX(i,k)>epsN)
614	if (QxPresent) VVQ(i,k)= calcVV()*ckQx1
615	if (present(massFlux3D)) massFlux3D(i,k)= VVQ(i,k)DE(i,k)QX(i,k) !(for use in MAIN)
616	enddo !k-loop
617	Vxmaxx(i)= min( VxMax, maxval(VVQ(i,:))*VincFact )
618
619	!note: dzMIN is min. value in column (not necessarily lowest layer in general)
620	dzMIN(i) = minval(DZ(i,:))
621	npassx(i)= max(1, nint( dtVxmaxx(i)/(CoMAXdzMIN(i)) ))
622	dtx(i) = dt/float(npassx(i))
623
624	!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -!
625	DO nnn= 1,npassx(i)
626
627	locallim = (nnn==1)
628
629	do k= ktop_sedi(i),nk !formerly do k= 1,nk
630	RHOQX(i,k) = DE(i,k)*QX(i,k)
631	QxPresent = (QX(i,k)>epsQ .and. NX(i,k)>epsN)
632	if (QxPresent) then
633	if (locallim) then !to avoid re-computing VVQ on first pass
634	VVQ(i,k)= -VVQ(i,k)
635	else
636	VVQ(i,k)= -calcVV()*ckQx1
637	endif
638	VVN(i,k)= VVQ(i,k)*ratio_Vn2Vq
639	VqMax = max(VxMAX,-VVQ(i,k))
640	VnMax = max(VxMAX,-VVN(i,k))
641	else
642	VVQ(i,k)= 0.
643	VVN(i,k)= 0.
644	VqMax = 0.
645	VnMax = 0.
646	endif
647	enddo !k-loop
648
649	!sum instantaneous surface mass flux at each split step: (for division later)
650	massFlux(i)= massFlux(i) - VVQ(i,nk)DE(i,nk)QX(i,nk)
651
652	!-- Perform single split sedimentation step:
653	! (formerly by calls to s/r 'blg4sedi', a modified [JM] version of 'blg2.ftn')
654	zz(ks)= 0.
655	do k= ks,ke(i),kw
656	zz(k+kw)= zz(k)+dz(i,k)
657	dzi(k) = 1./dz(i,k)
658	vp_Q(k) = 0.
659	vp_N(k) = 0.
660	enddo
661
662	do k=ks,ke(i),kw
663	zb_Q(k)= zz(k) + VVQ(i,k)*dtx(i)
664	zb_N(k)= zz(k) + VVN(i,k)*dtx(i)
665	enddo
666
667	zt_Q(ke(i))= zb_Q(ke(i)) + dz(i,ke(i))
668	zt_N(ke(i))= zb_N(ke(i)) + dz(i,ke(i))
669	do k= ks,ke(i)-kw,kw
670	zb_Q(k)= min(zb_Q(k+kw)-epsilondz(i,k), zz(k)+VVQ(i,k)dtx(i))
671	zb_N(k)= min(zb_N(k+kw)-epsilondz(i,k), zz(k)+VVN(i,k)dtx(i))
672	zt_Q(k)= zb_Q(k+kw)
673	zt_N(k)= zb_N(k+kw)
674	enddo
675
676	do k=ks,ke(i),kw !formerly k=1,nk
677	Q_star(k)= RHOQX(i,k)*dz(i,k)/(zt_Q(k)-zb_Q(k))
678	N_star(k)= NX(i,k)*dz(i,k)/(zt_N(k)-zb_N(k))
679	enddo
680
681	if (locallim) then
682	zmax_Q= abs(VqMax*dtx(i))
683	zmax_N= abs(VnMax*dtx(i))
684	do l=ks,ke(i),kw
685	flim_Q(l)= l
686	flim_N(l)= l
687	do k= l,ke(i),kw
688	if (zmax_Q.ge.zz(k)-zz(l+kw)) flim_Q(l)= k
689	if (zmax_N.ge.zz(k)-zz(l+kw)) flim_N(l)= k
690	enddo
691	enddo
692	endif
693
694	do l=ks,ke(i),kw
695	do k=l,flim_Q(l),kw
696	vp_Q(l)= vp_Q(l) + Q_star(k)*max(0.,min(zz(l+kw),zt_Q(k))-max(zz(l),zb_Q(k)))
697	enddo
698	do k=l,flim_N(l),kw
699	vp_N(l)= vp_N(l) + N_star(k)*max(0.,min(zz(l+kw),zt_N(k))-max(zz(l),zb_N(k)))
700	enddo
701	enddo
702
703	do k=ks,ke(i),kw
704	RHOQX(i,k)= vp_Q(k)*dzi(k)
705	NX(i,k)= vp_N(k)*dzi(k)
706	enddo
707	!--
708
709	do k= ktop_sedi(i),nk !formerly do k= 1,nk
710	QX(i,k)= RHOQX(i,k)*iDE(i,k)
711
712	!Prevent levels with zero N and nonzero Q and size-limiter:
713	QxPresent= (QX(i,k)>epsQ .and. NX(i,k)>epsN)
714	if (QxPresent) then !size limiter
715	Dx= (DE(i,k)QX(i,k)/(NX(i,k)cmx))**idmx
716	if (cat==1 .and. Dx>3.e-3) then
717	tmp1 = Dx-3.e-3; tmp1= tmp1*tmp1
718	tmp2 = (Dx/DxMAX); tmp2= tmp2tmp2tmp2
719	NX(i,k)= NX(i,k)max((1.+2.e4tmp1),tmp2)
720	else
721	NX(i,k)= NX(i,k)(max(Dx,DxMAX)iDxMAX)**dmx !impose Dx_max
722	endif
723	else !here, "QxPresent" implies correlated QX and NX
724	Q(i,k) = Q(i,k) + QX(i,k)
725	T(i,k) = T(i,k) - LXP*QX(i,k) !LCP for rain; LSP for i,s,g,h
726	QX(i,k)= 0.
727	NX(i,k)= 0.
728	endif
729
730	enddo
731
732	ENDDO !nnn-loop
733	!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -!
734	!compute average mass flux during the full time step: (used to compute the
735	!instantaneous sedimentation rate [liq. equiv. volume flux] in the main s/r)
736	massFlux(i)= massFlux(i)/float(npassx(i))
737
738	ENDDO !a(i)-loop
739	!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -!
740
741	CONTAINS
742
743	real function calcVV()
744	!Calculates portion of moment-weighted fall velocities
745	iLAMx = ((QX(i,k)DE(i,k)/NX(i,k))ckQx4)**idmx
746	iLAMxB0 = iLAMx**bfx
747	calcVV = gamfact(i,k)*iLAMxB0
748	end function calcVV
749
750	END SUBROUTINE SEDI_main_2
751
752	!=====================================================================================!
753	SUBROUTINE SEDI_main_1b(QX,cat,T,DE,iDE,gamfact,epsQ,afx,bfx,icmx,dmx,ckQx1,ckQx4, &
754	ni,nk,VxMax,DxMax,dt,DZ,massFlux,No_x,ktop_sedi,GRAV, &
755	massFlux3D)
756
757	!-------------------------------------------------------------------------------------!
758	! SINGLE-MOMENT version of sedimentation subroutine for categories whose
759	! fall velocity equation is V(D) = gamfact * afx * D^bfx
760	!-------------------------------------------------------------------------------------!
761
762	! Passing parameters:
763	!
764	! VAR Description
765	! --- ------------
766	! QX mass mixing ratio of category x
767	! cat: hydrometeor category:
768	! 1 rain
769	! 2 ice
770	! 3 snow
771	! 4 graupel
772	! 5 hail
773	!-------------------------------------------------------------------------------------!
774
775	use my_fncs_mod
776
777	implicit none
778
779	! PASSING PARAMETERS:
780	real, dimension(:,:), intent(inout) :: QX,T
781	real, dimension(:), intent(out) :: massFlux
782	real, optional, dimension(:,:), intent(out) :: massFlux3D
783	real, dimension(:,:), intent(in) :: DE,iDE,DZ
784	real, intent(in) :: epsQ,VxMax,afx,bfx,icmx,dmx,ckQx1,ckQx4,DxMax,dt,GRAV,No_x
785	integer, dimension(:), intent(in) :: ktop_sedi
786	integer, intent(in) :: ni,nk,cat !,ktop_sedi
787
788	! LOCAL PARAMETERS:
789	logical :: slabHASmass,locallim,QxPresent
790	integer :: nnn,a,i,k,counter,l,km1,kp1,ks,kw,idzmin !,ke
791	integer, dimension(nk) :: flim_Q
792	integer, dimension(ni) :: activeColumn,npassx,ke
793	real :: VqMax,iLAMx,iLAMxB0,tmp1,tmp2,Dx,iDxMax,VincFact,NX,iNo_x, &
794	zmax_Q,zmax_N,tempo
795	real, dimension(ni,nk) :: VVQ,RHOQX,gamfact
796	real, dimension(ni) :: dzMIN,dtx,VxMaxx
797	real, dimension(nk) :: vp_Q,zt_Q,zb_Q,dzi,Q_star
798	real, dimension(0:nk) :: zz
799	real, parameter :: epsilon = 1.e-2
800	real, parameter :: thrd = 1./3.
801	real, parameter :: sxth = 1./6.
802	real, parameter :: CoMAX = 2.0
803
804	!-------------------------------------------------------------------------------------!
805
806	massFlux= 0.
807	!Factor to estimate increased V from size-sorting:
808	! - this factor should be higher for categories with more time-splitting, since Vmax
809	! increases after each sedimentation split step (to be tuned)
810	VincFact= 1.
811	if (present(massFlux3D)) massFlux3D= 0. !(for use in MAIN for diagnostics)
812
813	!Determine for which slabs and columns sedimentation should be computes:
814	call countColumns(QX,ni,nk,epsQ,counter,activeColumn,ktop_sedi)
815	iNo_x = 1./No_x
816	iDxMax= 1./DxMax
817	ks = nk
818	ke = ktop_sedi !(i-array) - old: ke=1
819	kw = -1 !direction of vertical leveling
820
821	VVQ = 0.
822	VqMax= 0.
823
824	DO a= 1,counter
825	i= activeColumn(a)
826
827	VVQ(i,:) = 0.
828	do k= ktop_sedi(i),nk !do k= 1,nk
829	QxPresent = (QX(i,k)>epsQ)
830	! if (QxPresent) VVQ(i,k)= calcVV()*ckQx1
831
832	if (QxPresent) then
833	!ice:
834	if (cat==2) then
835	NX = 5.exp(0.304(273.15-max(233.,T(i,k))))
836	iLAMx = (ckQx4QX(i,k)DE(i,k)/NX)**thrd
837	!snow:
838	else if (cat==3) then
839	iNo_x = 1./min(2.e+8, 2.e+6exp(-0.12min(-0.001,T(i,k)-273.15)))
840	iLAMx = sqrt(sqrt(QX(i,k)DE(i,k)icmxsxthiNo_x))
841	!rain, graupel, hail:
842	else
843	iLAMx = sqrt(sqrt(QX(i,k)DE(i,k)icmxsxthiNo_x))
844	endif
845	VVQ(i,k) = -gamfact(i,k)ckQx1iLAMx**bfx
846	! VqMax = max(VxMAX,-VVQ(i,k))
847	endif
848	if (present(massFlux3D)) massFlux3D(i,k)= -VVQ(i,k)DE(i,k)QX(i,k) !(for use in MAIN)
849
850	enddo !k-loop
851	Vxmaxx(i)= min( VxMax, maxval(VVQ(i,:))*VincFact )
852
853	!note: dzMIN is min. value in column (not necessarily lowest layer in general)
854	dzMIN(i) = minval(DZ(i,:))
855	npassx(i)= max(1, nint( dtVxmaxx(i)/(CoMAXdzMIN(i)) ))
856	dtx(i) = dt/float(npassx(i))
857
858	!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -!
859	DO nnn= 1,npassx(i)
860
861	locallim = (nnn==1)
862
863	do k= ktop_sedi(i),nk !do k= 1,nk
864	RHOQX(i,k) = DE(i,k)*QX(i,k)
865	QxPresent = (QX(i,k)>epsQ)
866	if (QxPresent) then
867	!ice:
868	if (cat==2) then
869	NX = 5.exp(0.304(273.15-max(233.,T(i,k))))
870	iLAMx = (ckQx4QX(i,k)DE(i,k)/NX)**thrd
871	!snow:
872	else if (cat==3) then
873	iNo_x = 1./min(2.e+8, 2.e+6exp(-0.12min(-0.001,T(i,k)-273.15)))
874	iLAMx = sqrt(sqrt(QX(i,k)DE(i,k)icmxsxthiNo_x))
875	!rain, graupel, hail:
876	else
877	iLAMx = sqrt(sqrt(QX(i,k)DE(i,k)icmxsxthiNo_x))
878	endif
879	VVQ(i,k) = -gamfact(i,k)ckQx1iLAMx**bfx
880	VqMax = max(VxMAX,-VVQ(i,k))
881	endif
882
883	enddo !k-loop
884
885	!-- Perform single split sedimentation step: (formerly by calls to s/r 'blg4sedi')
886	zz(ks)= 0.
887	do k= ks,ke(i),kw
888	zz(k+kw)= zz(k)+dz(i,k)
889	dzi(k) = 1./dz(i,k)
890	vp_Q(k) = 0.
891	enddo
892
893	do k=ks,ke(i),kw
894	zb_Q(k)= zz(k) + VVQ(i,k)*dtx(i)
895	enddo
896
897	zt_Q(ke(i))= zb_Q(ke(i)) + dz(i,ke(i))
898	do k= ks,ke(i)-kw,kw
899	zb_Q(k)= min(zb_Q(k+kw)-epsilondz(i,k), zz(k)+VVQ(i,k)dtx(i))
900	zt_Q(k)= zb_Q(k+kw)
901	enddo
902
903	do k=ks,ke(i),kw !k=1,nk
904	Q_star(k)= RHOQX(i,k)*dz(i,k)/(zt_Q(k)-zb_Q(k))
905	enddo
906
907	if (locallim) then
908	zmax_Q= abs(VqMax*dtx(i))
909	do l=ks,ke(i),kw
910	flim_Q(l)= l
911	do k= l,ke(i),kw
912	if (zmax_Q.ge.zz(k)-zz(l+kw)) flim_Q(l)= k
913	enddo
914	enddo
915	endif
916
917	do l=ks,ke(i),kw
918	do k=l,flim_Q(l),kw
919	vp_Q(l)= vp_Q(l) + Q_star(k)*max(0.,min(zz(l+kw),zt_Q(k))-max(zz(l),zb_Q(k)))
920	enddo
921	enddo
922
923	do k=ks,ke(i),kw
924	RHOQX(i,k)= vp_Q(k)*dzi(k)
925	enddo
926	!--
927
928	do k= ktop_sedi(i),nk ! do k= 1,nk
929	QX(i,k)= RHOQX(i,k)*iDE(i,k)
930	enddo
931
932	!sum instantaneous flux at each split step: (for division later)
933	massFlux(i)= massFlux(i) - VVQ(i,nk)DE(i,nk)QX(i,nk)
934
935	ENDDO !nnn-loop
936	!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -!
937	!compute average flux during the full time step: (this will be used to compute
938	! the instantaneous sedimentation rate [volume flux] in the main s/r)
939	massFlux(i)= massFlux(i)/float(npassx(i))
940
941	ENDDO !a(i)-loop
942	!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -!
943
944	END SUBROUTINE SEDI_main_1b
945
946	!=====================================================================================!
947	SUBROUTINE countColumns(QX,ni,nk,minQX,counter,activeColumn,ktop_sedi)
948
949	! Searches the hydrometeor array QX(ni,nk) for non-zero (>minQX) values.
950	! Returns the array if i-indices (activeColumn) for the columns (i)
951	! which contain at least one non-zero value, as well as the number of such
952	! columns (counter).
953
954	implicit none
955
956	!PASSING PARAMETERS:
957	integer, intent(in) :: ni,nk !,ktop_sedi
958	integer, dimension(:), intent(in) :: ktop_sedi
959	integer, intent(out) :: counter
960	integer, dimension(:), intent(out) :: activeColumn
961	real, dimension(:,:), intent(in) :: QX
962	real, intent(in) :: minQX
963
964	!LOCAL PARAMETERS:
965	integer :: i !,k
966	integer, dimension(ni) :: k
967
968	! k= ktop_sedi-1 ! k=0
969	counter = 0
970	activeColumn= 0
971
972	do i=1,ni
973	k(i)= ktop_sedi(i)-1 ! k=0
974	do
975	k(i)=k(i)+1
976	if (QX(i,k(i))>minQX) then
977	counter=counter+1
978	activeColumn(counter)=i
979	k(i)=0
980	exit
981	else
982	if (k(i)==nk) then
983	k(i)=0
984	exit
985	endif
986	endif
987	enddo
988	enddo !i-loop
989
990	END SUBROUTINE countColumns
991
992	!=====================================================================================!
993
994	end module my_sedi_mod
995
996	!________________________________________________________________________________________!
997
998	module my_dmom_mod
999
1000	implicit none
1001
1002	private
1003	public :: mp_milbrandt2mom_main
1004
1005	contains
1006
1007	!_______________________________________________________________________________________!
1008
1009	SUBROUTINE mp_milbrandt2mom_main(W_omega,T,Q,QC,QR,QI,QN,QG,QH,NC,NR,NY,NN,NG,NH,PS,TM, &
1010	QM,QCM,QRM,QIM,QNM,QGM,QHM,NCM,NRM,NYM,NNM,NGM,NHM,PSM,S,RT_rn1,RT_rn2,RT_fr1,RT_fr2,&
1011	RT_sn1,RT_sn2,RT_sn3,RT_pe1,RT_pe2,RT_peL,RT_snd,GZ,T_TEND,Q_TEND,QCTEND,QRTEND, &
1012	QITEND,QNTEND,QGTEND,QHTEND,NCTEND,NRTEND,NYTEND,NNTEND,NGTEND,NHTEND,dt,NI,N,NK, &
1013	J,KOUNT,CCNtype,precipDiag_ON,sedi_ON,warmphase_ON,autoconv_ON,icephase_ON,snow_ON, &
1014	initN,dblMom_c,dblMom_r,dblMom_i,dblMom_s,dblMom_g,dblMom_h,Dm_c,Dm_r,Dm_i,Dm_s, &
1015	Dm_g,Dm_h,ZET,ZEC,SLW,VIS,VIS1,VIS2,VIS3,h_CB,h_ML1,h_ML2,h_SN,SS01,SS02,SS03,SS04, &
1016	SS05,SS06,SS07,SS08,SS09,SS10,SS11,SS12,SS13,SS14,SS15,SS16,SS17,SS18,SS19,SS20)
1017
1018
1019	use my_fncs_mod
1020	use my_sedi_mod
1021	!--WRF:
1022	use module_model_constants, ONLY: CPD => cp, CPV => cpv, RGASD => r_d, RGASV => r_v, &
1023	EPS1 => EP_2, DELTA => EP_1, CAPPA => rcp, GRAV => g, CHLC => XLV, CHLF => XLF
1024	!==
1025
1026	implicit none
1027
1028	!CALLING PARAMETERS:
1029	integer, intent(in) :: NI,NK,N,J,KOUNT,CCNtype
1030	real, intent(in) :: dt
1031	real, dimension(:), intent(in) :: PS,PSM
1032	real, dimension(:), intent(out) :: h_CB,h_ML1,h_ML2,h_SN
1033	real, dimension(:), intent(out) :: RT_rn1,RT_rn2,RT_fr1,RT_fr2,RT_sn1,RT_sn2, &
1034	RT_sn3,RT_pe1,RT_pe2,RT_peL,ZEC,RT_snd
1035	real, dimension(:,:), intent(in) :: W_omega,S,GZ
1036	real, dimension(:,:), intent(inout) :: T,Q,QC,QR,QI,QN,QG,QH,NC,NR,NY,NN,NG,NH, &
1037	TM,QM,QCM,QRM,QIM,QNM,QGM,QHM,NCM,NRM,NYM,NNM,NGM,NHM
1038	real, dimension(:,:), intent(out) :: T_TEND,QCTEND,QRTEND,QITEND,QNTEND, &
1039	QGTEND,QHTEND,Q_TEND,NCTEND,NRTEND,NYTEND,NNTEND,NGTEND,NHTEND,ZET,Dm_c, &
1040	Dm_r,Dm_i,Dm_s,Dm_g,Dm_h,SLW,VIS,VIS1,VIS2,VIS3,SS01,SS02,SS03,SS04,SS05,SS06, &
1041	SS07,SS08,SS09,SS10,SS11,SS12,SS13,SS14,SS15,SS16,SS17,SS18,SS19,SS20
1042	logical, intent(in) :: dblMom_c,dblMom_r,dblMom_i,dblMom_s, &
1043	dblMom_g,dblMom_h,precipDiag_ON,sedi_ON,icephase_ON,snow_ON,warmphase_ON, &
1044	autoconv_ON,initN
1045
1046	!_______________________________________________________________________________________
1047	! !
1048	! Milbrandt-Yau Multimoment Bulk Microphysics Scheme !
1049	! - double-moment version - !
1050	!_______________________________________________________________________________________!
1051	! Package version: 2.19.0 (internal bookkeeping) !
1052	! Last modified : 2011-03-02 !
1053	!_______________________________________________________________________________________!
1054	!
1055	! Author:
1056	! J. Milbrandt, McGill University (August 2004)
1057	!
1058	! Major revisions:
1059	!
1060	! 001 J. Milbrandt (Dec 2006) - Converted the full Milbrandt-Yau (2005) multimoment
1061	! (RPN) scheme to an efficient fixed-dispersion double-moment
1062	! version
1063	! 002 J. Milbrandt (Mar 2007) - Added options for single-moment/double-moment for
1064	! each hydrometeor category
1065	! 003 J. Milbrandt (Feb 2008) - Modified single-moment version for use in GEM-LAM-2.5
1066	! 004 J. Milbrandt (Nov 2008) - Modified double-moment version for use in 2010 Vancouver
1067	! Olympics GEM-LAM configuration
1068	! 005 J. Milbrandt (Aug 2009) - Modified (dmom) for PHY_v5.0.4, for use in V2010 system:
1069	! + reduced ice/snow capacitance to C=0.25D (from C=0.5D)
1070	! + added diagnostic fields (VIS, levels, etc.)
1071	! + added constraints to snow size distribution (No_s and
1072	! LAMDA_s limits, plus changed m-D parameters
1073	! + modified solid-to-liquid ratio calculation, based on
1074	! volume flux (and other changes)
1075	! + added back sedimentation of ice category
1076	! + modified condition for conversion of graupel to hail
1077	! + corrected bug it diagnostic "ice pellets" vs. "hail"
1078	! + minor bug corrections (uninitialized values, etc.)
1079	! 006 J. Milbrandt (Jan 2011) - Bug fixes and minor code clean-up from PHY_v5.1.3 version
1080	! + corrected latent heat constants in thermodynamic functions
1081	! (ABi and ABw) for sublimation and evaporation
1082	! + properly initialized variables No_g and No_h
1083	! + changed max ice crystal size (fallspeed) to 5 mm (2 m s-1)
1084	! + imposed maximum ice number concentration of 1.e+7 m-3
1085	! + removed unused supersaturation reduction
1086	!
1087	! Object:
1088	! Computes changes to the temperature, water vapor mixing ratio, and the
1089	! mixing ratios and total number concentrations of six hydrometeor species
1090	! resulting from cloud microphysical interactions at saturated grid points.
1091	! Liquid and solid surface precipitation rates from sedimenting hydrometeor
1092	! categories are also computed.
1093	!
1094	! This subroutine and the associated modules form the single/double-moment
1095	! switchable verion of the multimoment bulk microphysics package, the full
1096	! version of which is described in the references below.
1097	!
1098	! References: Milbrandt and Yau, (2005a), J. Atmos. Sci., vol.62, 3051-3064
1099	! --------- and ---, (2005b), J. Atmos. Sci., vol.62, 3065-3081
1100	! (and references therein)
1101	!
1102	! Please report bugs to: jason.milbrandt@ec.gc.ca
1103	!_______________________________________________________________________________________!
1104	!
1105	! Arguments: Description: Units:
1106	!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -!
1107	! - Input -
1108	!
1109	! NI number of x-dir points (in local subdomain)
1110	! NK number of vertical levels
1111	! N not used (to be removed)
1112	! J y-dir index (local subdomain)
1113	! KOUNT current model time step number
1114	! dt model time step [s]
1115	! CCNtype switch for airmass type
1116	! 1 = maritime --> N_c = 80 cm-3 (1-moment cloud)
1117	! 2 = continental 1 --> N_c = 200 cm-3 " "
1118	! 3 = continental 2 (polluted) --> N_c = 500 cm-3 " "
1119	! 4 = land-sea-mask-dependent (TBA)
1120	! W_omega vertical velocity [Pa s-1]
1121	! S sigma (=p/p_sfc)
1122	! GZ geopotential
1123	! dblMom_(x) logical switch for double(T)-single(F)-moment for category (x)
1124	! precipDiag_ON logical switch, .F. to suppress calc. of sfc precip types
1125	! sedi_ON logical switch, .F. to suppress sedimentation
1126	! warmphase_ON logical switch, .F. to suppress warm-phase (Part II)
1127	! autoconv_ON logical switch, .F. to supppress autoconversion (cld->rn)
1128	! icephase_ON logical switch, .F. to suppress ice-phase (Part I)
1129	! snow_ON logical switch, .F. to suppress snow initiation
1130	!
1131	!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -!
1132	! - Input/Output -
1133	!
1134	! T air temperature at time (t*) [K]
1135	! TM air temperature at time (t-dt) [K]
1136	! Q water vapor mixing ratio at (t*) [kg kg-1]
1137	! QM water vapor mixing ratio at (t-dt) [kg kg-1]
1138	! PS surface pressure at time (t*) [Pa]
1139	! PSM surface pressure at time (t-dt) [Pa]
1140	!
1141	! For x = (C,R,I,N,G,H): C = cloud
1142	! R = rain
1143	! I = ice (pristine) [except 'NY', not 'NI']
1144	! N = snow
1145	! G = graupel
1146	! H = hail
1147	!
1148	! Q(x) mixing ratio for hydrometeor x at (t*) [kg kg-1]
1149	! Q(x)M mixing ratio for hydrometeor x at (t-dt) [kg kg-1]
1150	! N(x) total number concentration for hydrometeor x (t*) [m-3]
1151	! N(x)M total number concentration for hydrometeor x (t-dt) [m-3]
1152	!
1153	! Note: The arrays "VM" (e.g. variables TM,QM,QCM etc.) are declared as INTENT(INOUT)
1154	! such that their values are modified in the code [VM = 0.5*(VM + V)].
1155	! This is to approxiate the values at time level (t), which are needed by
1156	! this routine but are unavailable to the PHYSICS. The new values are discared
1157	! by the calling routine ('vkuocon6.ftn'). However, care should be taken with
1158	! interfacing with other modelling systems. For GEM/MC2, it does not matter if
1159	! VM is modified since the calling module passes back only the tendencies
1160	! (VTEND) to the model.
1161
1162	!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -!
1163	! - Output -
1164	!
1165	! Q_TEND tendency for water vapor mixing ratio [kg kg-1 s-1]
1166	! T_TEND tendency for air temperature [K s-1]
1167	! Q(x)TEND tendency for mixing ratio for hydrometeor x [kg kg-1 s-1]
1168	! N(x)TEND tendency for number concentration for hydrometeor x [m-3 s-1]
1169	! Dm_(x) mean-mass diameter for hydrometeor x [m]
1170	! H_CB height of cloud base [m]
1171	! h_ML1 height of first melting level from ground [m]
1172	! h_ML2 height of first melting level from top [m]
1173	! h_SN height of snow level [m]
1174	! RT_rn1 precipitation rate (at sfc) of liquid rain [m+3 m-2 s-1]
1175	! RT_rn2 precipitation rate (at sfc) of liquid drizzle [m+3 m-2 s-1]
1176	! RT_fr1 precipitation rate (at sfc) of freezing rain [m+3 m-2 s-1]
1177	! RT_fr2 precipitation rate (at sfc) of freezing drizzle [m+3 m-2 s-1]
1178	! RT_sn1 precipitation rate (at sfc) of ice crystals (liq-eq) [m+3 m-2 s-1]
1179	! RT_sn2 precipitation rate (at sfc) of snow (liq-equiv) [m+3 m-2 s-1]
1180	! RT_sn3 precipitation rate (at sfc) of graupel (liq-equiv) [m+3 m-2 s-1]
1181	! RT_snd precipitation rate (at sfc) of snow (frozen) [m+3 m-2 s-1]
1182	! RT_pe1 precipitation rate (at sfc) of ice pellets (liq-eq) [m+3 m-2 s-1]
1183	! RT_pe2 precipitation rate (at sfc) of hail (total; liq-eq) [m+3 m-2 s-1]
1184	! RT_peL precipitation rate (at sfc) of hail (large only) [m+3 m-2 s-1]
1185	! SSxx S/S terms (for testing purposes)
1186	! SLW supercooled liquid water content [kg m-3]
1187	! VIS visibility resulting from fog, rain, snow [m]
1188	! VIS1 visibility component through liquid cloud (fog) [m]
1189	! VIS2 visibility component through rain [m]
1190	! VIS3 visibility component through snow [m]
1191	! ZET total equivalent radar reflectivity [dBZ]
1192	! ZEC composite (column-max) of ZET [dBZ]
1193	!_______________________________________________________________________________________!
1194
1195
1196	!LOCAL VARIABLES:
1197
1198	!Variables to count active grid points:
1199	logical :: log1,log2,log3,log4,doneK,rainPresent,calcDiag,CB_found,ML_found, &
1200	SN_found
1201	logical, dimension(size(QC,dim=1),size(QC,dim=2)) :: activePoint
1202	integer, dimension(size(QC,dim=1)) :: ktop_sedi
1203	integer :: i,k,niter,ll,start
1204
1205	real :: tmp1,tmp2,tmp3,tmp4,tmp5,tmp6,tmp7,tmp8,tmp9,tmp10, &
1206	VDmax,NNUmax,X,D,DEL,QREVP,NuDEPSOR,NuCONTA,NuCONTB,NuCONTC,iMUkin,Ecg,Erg, &
1207	NuCONT,GG,Na,Tcc,F1,F2,Kdiff,PSIa,Kn,source,sink,sour,ratio,qvs0,Kstoke, &
1208	DELqvs,ft,esi,Si,Simax,Vq,Vn,Vz,LAMr,No_r_DM,No_i,No_s,No_g,No_h,D_sll, &
1209	iABi,ABw,VENTr,VENTs,VENTg,VENTi,VENTh,Cdiff,Ka,MUdyn,MUkin,DEo,Ng_tail, &
1210	gam,ScTHRD,Tc,mi,ff,Ec,Ntr,Dho,DMrain,Ech,DMice,DMsnow,DMgrpl,DMhail, &
1211	ssat,Swmax,dey,Esh,Eii,Eis,Ess,Eig,Eih,FRAC,JJ,Dirg,Dirh,Dsrs,Dsrg,Dsrh, &
1212	Dgrg,Dgrh,SIGc,L,TAU,DrAUT,DrINIT,Di,Ds,Dg,Dh,qFact,nFact,Ki,Rz,NgCNgh, &
1213	vr0,vi0,vs0,vg0,vh0,Dc,Dr,QCLcs,QCLrs,QCLis,QCLcg,QCLrg,QCLig,NhCNgh, &
1214	QCLch,QCLrh,QCLsh,QMLir,QMLsr,QMLgr,QMLhr,QCLih,QVDvg,QVDvh,QSHhr, &
1215	QFZci,QNUvi,QVDvi,QCNis,QCNis1,QCNis2,QCLir,QCLri,QCNsg,QCLsr,QCNgh, &
1216	QCLgr,QHwet,QVDvs,QFZrh,QIMsi,QIMgi,NMLhr,NVDvh,NCLir,NCLri,NCLrh, &
1217	NCLch,NCLsr,NCLirg,NCLirh,NrFZrh,NhFZrh,NCLsrs,NCLsrg,NCLsrh,NCLgrg, &
1218	NCLgrh,NVDvg,NMLgr,NiCNis,NsCNis,NVDvs,NMLsr,NCLsh,NCLss,NNUvi,NFZci,NVDvi, &
1219	NCLis,NCLig,NCLih,NMLir,NCLrs,NCNsg,NCLcs,NCLcg,NIMsi,NIMgi,NCLgr,NCLrg, &
1220	NSHhr,RCAUTR,RCACCR,CCACCR,CCSCOC,CCAUTR,CRSCOR,ALFx,des_pmlt,Ecs,des,ides, &
1221	LAMx,iLAMx,iLAMxB0,Dx,ffx,iLAMc,iNCM,iNRM,iNYM,iNNM,iNGM,iLAMs_D3, &
1222	iLAMg,iLAMg2,iLAMgB0,iLAMgB1,iLAMgB2,iLAMh,iLAMhB0,iLAMhB1,iLAMhB2,iNHM, &
1223	iLAMi,iLAMi2,iLAMi3,iLAMi4,iLAMi5,iLAMiB0,iLAMiB1,iLAMiB2,iLAMr6,iLAMh2, &
1224	iLAMs,iLAMs2,iLAMsB0,iLAMsB1,iLAMsB2,iLAMr,iLAMr2,iLAMr3,iLAMr4,iLAMr5, &
1225	iLAMc2,iLAMc3,iLAMc4,iLAMc5,iLAMc6,iQCM,iQRM,iQIM,iQNM,iQGM,iQHM,iEih,iEsh, &
1226	N_c,N_r,N_i,N_s,N_g,N_h,fluxV_i,fluxV_g,fluxV_s,rhos_mlt,fracLiq
1227
1228	!Variables that only need to be calulated on the first step (and saved):
1229	real, save :: idt,iMUc,cmr,cmi,cms,cmg,cmh,icmr,icmi,icmg,icms,icmh,idew,idei, &
1230	ideh,ideg,GC1,imso,icexc9,cexr1,cexr2,cexr3,No_s_SM,No_r,idms,imgo,icexs2, &
1231	cexr4,cexr5,cexr6,cexr9,icexr9,ckQr1,ckQr2,ckQr3,ckQi1,ckQi2,ckQi3,ckQi4, &
1232	icexi9,ckQs1,ckQs2,cexs1,cexs2,ckQg1,ckQg2,ckQg4,ckQh1,ckQh2,ckQh4,GR37,dms, &
1233	LCP,LFP,LSP,ck5,ck6,PI2,PIov4,PIov6,CHLS,iCHLF,cxr,cxi,Gzr,Gzi,Gzs,Gzg,Gzh, &
1234	N_c_SM,iGC1,GC2,GC3,GC4,GC5,iGC5,GC6,GC7,GC8,GC11,GC12,GC13,GC14,iGR34,mso, &
1235	GC15,GR1,GR3,GR13,GR14,GR15,GR17,GR31,iGR31,GR32,GR33,GR34,GR35,GR36,GI4, &
1236	GI6,GI20,GI21,GI22,GI31,GI32,GI33,GI34,GI35,iGI31,GI11,GI36,GI37,GI40,iGG34, &
1237	GS09,GS11,GS12,GS13,iGS20,GS31,iGS31,GS32,GS33,GS34,GS35,GS36,GS40,iGS40, &
1238	GS50,GG09,GG11,GG12,GG13,GG31,iGG31,GG32,GG33,GG34,GG35,GG36,GG40,iGG99,GH09,&
1239	GH11,GH12,GH13,GH31,GH32,GH33,GH40,GR50,GG50,iGH34,GH50,iGH99,iGH31,iGS34, &
1240	iGS20_D3,GS40_D3,cms_D3,eds,fds,rfact_FvFm
1241
1242	!Size distribution parameters:
1243	real, parameter :: MUc = 3. !shape parameter for cloud
1244	real, parameter :: alpha_c = 1. !shape parameter for cloud
1245	real, parameter :: alpha_r = 0. !shape parameter for rain
1246	real, parameter :: alpha_i = 0. !shape parameter for ice
1247	real, parameter :: alpha_s = 0. !shape parameter for snow
1248	real, parameter :: alpha_g = 0. !shape parameter for graupel
1249	real, parameter :: alpha_h = 0. !shape parameter for hail
1250	real, parameter :: No_s_max = 1.e+8 !max. allowable intercept for snow [m-4]
1251	real, parameter :: lamdas_min= 500. !min. allowable LAMDA_s [m-1]
1252
1253	!For single-moment:
1254	real, parameter :: No_r_SM = 1.e+7 !intercept parameter for rain [m-4]
1255	real, parameter :: No_g_SM = 4.e+6 !intercept parameter for graupel [m-4]
1256	real, parameter :: No_h_SM = 1.e+5 !intercept parameter for hail [m-4]
1257	!note: No_s = f(T), rather than a fixed value
1258	!------------------------------------!
1259	! Symbol convention: (dist. params.) ! MY05: Milbrandt & Yau, 2005a,b (JAS)
1260	! MY05 F94 CP00 ! F94: Ferrier, 1994 (JAS)
1261	! ------ -------- ------ ! CP00: Cohard & Pinty, 2000a,b (QJGR)
1262	! ALFx ALPHAx MUx-1 !
1263	! MUx (1) ALPHAx !
1264	! ALFx+1 ALPHAx+1 MUx !
1265	!------------------------------------!
1266	! Note: The symbols for MU and ALPHA are REVERSED from that of CP2000a,b
1267	! Explicit appearance of MUr = 1. has been removed.
1268
1269	! Fallspeed parameters:
1270	real, parameter :: afr= 149.100, bfr= 0.5000 !Tripoloi and Cotton (1980)
1271	real, parameter :: afi= 71.340, bfi= 0.6635 !Ferrier (1994)
1272	real, parameter :: afs= 11.720, bfs= 0.4100 !Locatelli and Hobbs (1974)
1273	real, parameter :: afg= 19.300, bfg= 0.3700 !Ferrier (1994)
1274	real, parameter :: afh= 206.890, bfh= 0.6384 !Ferrier (1994)
1275	!options:
1276	!real, parameter :: afs= 8.996, bfs= 0.4200 !Ferrier (1994)
1277	!real, parameter :: afg= 6.4800, bfg= 0.2400 !LH74 (grpl-like snow of lump type)
1278
1279	real, parameter :: epsQ = 1.e-14 !kg kg-1, min. allowable mixing ratio
1280	real, parameter :: epsN = 1.e-3 !m-3, min. allowable number concentration
1281	real, parameter :: epsQ2 = 1.e-6 !kg kg-1, mixing ratio threshold for diagnostics
1282	real, parameter :: epsVIS= 1. !m, min. allowable visibility
1283
1284	real, parameter :: iLAMmin1= 1.e-6 !min. iLAMx (prevents underflow in Nox and VENTx calcs)
1285	real, parameter :: iLAMmin2= 1.e-10 !min. iLAMx (prevents underflow in Nox and VENTx calcs)
1286	real, parameter :: eps = 1.e-32
1287	real, parameter :: k1 = 0.001
1288	real, parameter :: k2 = 0.0005
1289	real, parameter :: k3 = 2.54
1290	real, parameter :: CPW = 4218., CPI=2093.
1291
1292	real, parameter :: deg = 400., mgo= 1.6e-10
1293	real, parameter :: deh = 900.
1294	real, parameter :: dei = 500., mio=1.e-12, Nti0=1.e3
1295	real, parameter :: dew = 1000.
1296	real, parameter :: desFix= 100. !used for snowSpherical = .true.
1297	real, parameter :: desMax= 500.
1298	real, parameter :: Dso = 125.e-6 ![m]; embryo snow diameter (mean-volume particle)
1299	real, parameter :: dmr = 3., dmi= 3., dmg= 3., dmh= 3.
1300
1301	! NOTE: VxMAX below are the max.allowable mass-weighted fallspeeds for sedimentation.
1302	! Thus, Vx corresponds to DxMAX (at sea-level) times the max. density factor, GAM.
1303	! [GAMmax=sqrt(DEo/DEmin)=sqrt(1.25/0.4)~2.] e.g. VrMAX = 2.*8.m/s = 16.m/s
1304	real, parameter :: DrMax= 5.e-3, VrMax= 16., epsQr_sedi= 1.e-8
1305	real, parameter :: DiMax= 5.e-3, ViMax= 2., epsQi_sedi= 1.e-10
1306	real, parameter :: DsMax= 5.e-3, VsMax= 2., epsQs_sedi= 1.e-8
1307	real, parameter :: DgMax= 50.e-3, VgMax= 8., epsQg_sedi= 1.e-8
1308	real, parameter :: DhMax= 80.e-3, VhMax= 25., epsQh_sedi= 1.e-10
1309
1310	real, parameter :: thrd = 1./3.
1311	real, parameter :: sixth = 0.5*thrd
1312	real, parameter :: Ers = 1., Eci= 1. !collection efficiencies, Exy, between categories x and y
1313	real, parameter :: Eri = 1., Erh= 1.
1314	real, parameter :: Xdisp = 0.25 !dispersion of the fall velocity of ice
1315	real, parameter :: aa11 = 9.44e15, aa22= 5.78e3, Rh= 41.e-6
1316	real, parameter :: Avx = 0.78, Bvx= 0.30 !ventilation coefficients [F94 (B.36)]
1317	real, parameter :: Abigg = 0.66, Bbigg= 100. !parameters in probabilistic freezing
1318	real, parameter :: fdielec = 4.464 !ratio of dielectric factor, \|K\|w2/\|K\|i2
1319	real, parameter :: zfact = 1.e+18 !conversion factor for m-3 to mm2 m-6 for Ze
1320	real, parameter :: minZET = -99. ![dBZ] min threshold for ZET
1321	real, parameter :: maxVIS = 99.e+3 ![m] max. allowable VIS (visibility)
1322	real, parameter :: Drshed = 0.001 ![m] mean diam. of drop shed during wet growth
1323	real, parameter :: SIGcTHRS = 15.e-6 !threshold cld std.dev. before autoconversion
1324	real, parameter :: KK1 = 3.03e3 !parameter in Long (1974) kernel
1325	real, parameter :: KK2 = 2.59e15 !parameter in Long (1974) kernel
1326	real, parameter :: Dhh = 82.e-6 ![m] diameter that rain hump first appears
1327	real, parameter :: gzMax_sedi = 200000. !GZ value below which sedimentation is computed
1328	real, parameter :: Dr_large = 200.e-6 ![m] size threshold to distinguish rain/drizzle for precip rates
1329	real, parameter :: Ds_large = 200.e-6 ![m] size threshold to distinguish snow/snow-grains for precip rates
1330	real, parameter :: Dh_large = 1.0e-2 ![m] size threshold for "large" hail precipitation rate
1331	real, parameter :: Dh_min = 5.0e-3 ![m] size threhsold for below which hail converts to graupel
1332	real, parameter :: Dr_3cmpThrs = 2.5e-3 ![m] size threshold for hail production from 3-comp freezing
1333	real, parameter :: w_CNgh = 3. ![m s-1] vertical motion threshold for CNgh
1334	! real, parameter :: r_CNgh = 0.05 !Dg/Dho ratio threshold for CNgh
1335	real, parameter :: Ngh_crit = 0.01 ![m-3] critical graupel concentration for CNgh
1336	real, parameter :: Tc_FZrh = -10. !temp-threshold (C) for FZrh
1337	real, parameter :: CNsgThres = 1.0 !threshold for CLcs/VDvs ratio for CNsg
1338	real, parameter :: capFact_i = 0.5 !capacitace factor for ice (C= 0.5DcapFact_i)
1339	real, parameter :: capFact_s = 0.5 !capacitace factor for snow (C= 0.5DcapFact_s)
1340	real, parameter :: noVal_h_XX = -1. !non-value indicator for h_CB, h_ML1, h_ML2, h_SN
1341	real, parameter :: minSnowSize = 1.e-4 ![m] snow size threshold to compute h_SN
1342	real, parameter :: Fv_Dsmin = 125.e-6 ![m] min snow size to compute volume flux
1343	real, parameter :: Fv_Dsmax = 0.008 ![m] max snow size to compute volume flux
1344	real, parameter :: Ni_max = 1.e+7 ![m-3] max ice crystal concentration
1345
1346	!-- For GEM:
1347	!#include "consphy.cdk"
1348	!#include "dintern.cdk"
1349	!#include "fintern.cdk"
1350
1351	!-- For WRF:
1352	!------------------------------------------------------------------------------!
1353	!#include "consphy.cdk"
1354	! real, parameter :: CPD =.100546e+4 !J K-1 kg-1; specific heat of dry air
1355	! real, parameter :: CPV =.186946e+4 !J K-1 kg-1; specific heat of water vapour
1356	! real, parameter :: RGASD =.28705e+3 !J K-1 kg-1; gas constant for dry air
1357	! real, parameter :: RGASV =.46151e+3 !J K-1 kg-1; gas constant for water vapour
1358	real, parameter :: TRPL =.27316e+3 !K; triple point of water
1359	real, parameter :: TCDK =.27315e+3 !conversion from kelvin to celsius
1360	real, parameter :: RAUW =.1e+4 !density of liquid H2O
1361	! real, parameter :: EPS1 =.62194800221014 !RGASD/RGASV
1362	real, parameter :: EPS2 =.3780199778986 !1 - EPS1
1363	! real, parameter :: DELTA =.6077686814144 !1/EPS1 - 1
1364	! real, parameter :: CAPPA =.28549121795 !RGASD/CPD
1365	real, parameter :: TGL =.27316e+3 !K; ice temperature in the atmosphere
1366	real, parameter :: CONSOL =.1367e+4 !W m-2; solar constant
1367	! real, parameter :: GRAV =.980616e+1 !M s-2; gravitational acceleration
1368	real, parameter :: RAYT =.637122e+7 !M; mean radius of the earth
1369	real, parameter :: STEFAN =.566948e-7 !J m-2 s-1 K-4; Stefan-Boltzmann constant
1370	real, parameter :: PI =.314159265359e+1 !PI constant = ACOS(-1)
1371	real, parameter :: OMEGA =.7292e-4 !s-1; angular speed of rotation of the earth
1372	real, parameter :: KNAMS =.514791 !conversion from knots to m/s
1373	real, parameter :: STLO =.6628486583943e-3 !K s2 m-2; Schuman-Newell Lapse Rate
1374	real, parameter :: KARMAN =.35 !Von Karman constant
1375	real, parameter :: RIC =.2 !Critical Richardson number
1376	! real, parameter :: CHLC =.2501e+7 !J kg-1; latent heat of condensation
1377	! real, parameter :: CHLF =.334e+6 !J kg-1; latent heat of fusion
1378
1379	!------------------------------------------------------------------------------!
1380	!#include "dintern.cdk"
1381	REAL TTT, PRS, QQQ, EEE, TVI, QST, QQH
1382	REAL T00, PR0, TF, PF,FFF , DDFF
1383	REAL QSM , DLEMX
1384	REAL*8 FOEW,FODLE,FOQST,FODQS,FOEFQ,FOQFE,FOTVT,FOTTV,FOHR
1385	REAL*8 FOLV,FOLS,FOPOIT,FOPOIP,FOTTVH,FOTVHT
1386	REAL*8 FOEWA,FODLA,FOQSA,FODQA,FOHRA
1387	REAL*8 FESI,FDLESI,FESMX,FDLESMX,FQSMX,FDQSMX
1388
1389	!------------------------------------------------------------------------------!
1390	!#include "fintern.cdk"
1391	! DEFINITION DES FONCTIONS THERMODYNAMIQUES DE BASE
1392	! POUR LES CONSTANTES, UTILISER LE COMMON /CONSPHY/
1393	! NOTE: TOUTES LES FONCTIONS TRAVAILLENT AVEC LES UNITES S.I.
1394	! I.E. TTT EN DEG K, PRS EN PA, QQQ EN KG/KG
1395	! * N. BRUNET - MAI 90 *
1396	! * REVISION 01 - MAI 94 - N. BRUNET
1397	! NOUVELLE VERSION POUR FAIBLES PRESSIONS
1398	! * REVISION 02 - AOUT 2000 - J-P TOVIESSI
1399	! CALCUL EN REAL*8
1400	! * REVISION 03 - SEPT 2000 - N. BRUNET
1401	! AJOUT DE NOUVELLES FONCTIONS
1402	! * REVISION 04 - JANV 2000 - J. MAILHOT
1403	! FONCTIONS EN PHASE MIXTE
1404	! * REVISION 05 - DEC 2001 - G. LEMAY
1405	! DOUBLE PRECISION POUR PHASE MIXTE
1406	! * REVISION 06 - AVR 2002 - A. PLANTE
1407	! AJOUT DES NOUVELLES FONCTIONS FOTTVH ET FOTVHT
1408	!
1409	! FONCTION DE TENSION DE VAPEUR SATURANTE (TETENS) - EW OU EI SELON TT
1410	FOEW(TTT) = 610.78D0*DEXP( DMIN1(DSIGN(17.269D0, &
1411	DBLE(TTT)-DBLE(TRPL)),DSIGN &
1412	(21.875D0,DBLE(TTT)-DBLE(TRPL)))*DABS(DBLE(TTT)-DBLE(TRPL))/ &
1413	(DBLE(TTT)-35.86D0+DMAX1(0.D0,DSIGN &
1414	(28.2D0,DBLE(TRPL)-DBLE(TTT)))))
1415	!
1416	! FONCTION CALCULANT LA DERIVEE SELON T DE LN EW (OU LN EI)
1417	FODLE(TTT)=(4097.93D0+DMAX1(0.D0,DSIGN(1709.88D0, &
1418	DBLE(TRPL)-DBLE(TTT)))) &
1419	/((DBLE(TTT)-35.86D0+DMAX1(0.D0,DSIGN(28.2D0, &
1420	DBLE(TRPL)-DBLE(TTT))))*(DBLE(TTT)-35.86D0+DMAX1(0.D0 &
1421	,DSIGN(28.2D0,DBLE(TRPL)-DBLE(TTT)))))
1422	!
1423	! FONCTION CALCULANT L'HUMIDITE SPECIFIQUE SATURANTE (QSAT)
1424	FOQST(TTT,PRS) = DBLE(EPS1)/(DMAX1(1.D0,DBLE(PRS)/FOEW(TTT))- &
1425	DBLE(EPS2))
1426	!
1427	! FONCTION CALCULANT LA DERIVEE DE QSAT SELON T
1428	FODQS(QST,TTT)=DBLE(QST)(1.D0+DBLE(DELTA)DBLE(QST))*FODLE(TTT)
1429	! QST EST LA SORTIE DE FOQST
1430	!
1431	! FONCTION CALCULANT TENSION VAP (EEE) FN DE HUM SP (QQQ) ET PRS
1432	FOEFQ(QQQ,PRS) = DMIN1(DBLE(PRS),(DBLE(QQQ)*DBLE(PRS)) / &
1433	(DBLE(EPS1) + DBLE(EPS2)*DBLE(QQQ)))
1434	!
1435	! FONCTION CALCULANT HUM SP (QQQ) DE TENS. VAP (EEE) ET PRES (PRS)
1436	FOQFE(EEE,PRS) = DMIN1(1.D0,DBLE(EPS1)*DBLE(EEE)/(DBLE(PRS)- &
1437	DBLE(EPS2)*DBLE(EEE)))
1438	!
1439	! FONCTION CALCULANT TEMP VIRT. (TVI) DE TEMP (TTT) ET HUM SP (QQQ)
1440	FOTVT(TTT,QQQ) = DBLE(TTT) * (1.0D0 + DBLE(DELTA)*DBLE(QQQ))
1441
1442	! FONCTION CALCULANT TEMP VIRT. (TVI) DE TEMP (TTT), HUM SP (QQQ) ET
1443	! MASSE SP DES HYDROMETEORES.
1444	FOTVHT(TTT,QQQ,QQH) = DBLE(TTT) * &
1445	(1.0D0 + DBLE(DELTA)*DBLE(QQQ) - DBLE(QQH))
1446	!
1447	! FONCTION CALCULANT TTT DE TEMP VIRT. (TVI) ET HUM SP (QQQ)
1448	FOTTV(TVI,QQQ) = DBLE(TVI) / (1.0D0 + DBLE(DELTA)*DBLE(QQQ))
1449
1450	! FONCTION CALCULANT TTT DE TEMP VIRT. (TVI), HUM SP (QQQ) ET
1451	! MASSE SP DES HYDROMETEORES (QQH)
1452	FOTTVH(TVI,QQQ,QQH) = DBLE(TVI) / &
1453	(1.0D0 + DBLE(DELTA)*DBLE(QQQ) - DBLE(QQH))
1454	!
1455	! FONCTION CALCULANT HUM REL DE HUM SP (QQQ), TEMP (TTT) ET PRES (PRS)
1456	! HR = E/ESAT
1457	#if (DWORDSIZE == 8 && RWORDSIZE == 8)
1458	FOHR(QQQ,TTT,PRS) = MIN( PRS ,FOEFQ(QQQ,PRS)) / FOEW(TTT)
1459	#elif (DWORDSIZE == 8 && RWORDSIZE == 4)
1460	FOHR(QQQ,TTT,PRS) = MIN(DBLE(PRS),FOEFQ(QQQ,PRS)) / FOEW(TTT)
1461	#else
1462	This is a temporary hack assuming double precision is 8 bytes.
1463	#endif
1464	!
1465	! FONCTION CALCULANT LA CHALEUR LATENTE DE CONDENSATION
1466	FOLV(TTT) =DBLE(CHLC) - 2317.D0*(DBLE(TTT)-DBLE(TRPL))
1467	!
1468	! FONCTION CALCULANT LA CHALEUR LATENTE DE SUBLIMATION
1469	FOLS(TTT) = DBLE(CHLC)+DBLE(CHLF)+(DBLE(CPV)- &
1470	(7.24D0DBLE(TTT)+128.4D0))(DBLE(TTT)-DBLE(TRPL))
1471	!
1472	! FONCTION RESOLVANT L'EQN. DE POISSON POUR LA TEMPERATURE
1473	! NOTE: SI PF=1000*100, "FOPOIT" DONNE LE THETA STANDARD
1474	FOPOIT(T00,PR0,PF)=DBLE(T00)(DBLE(PR0)/DBLE(PF))* &
1475	(-DBLE(CAPPA))
1476	!
1477	! FONCTION RESOLVANT L'EQN. DE POISSON POUR LA PRESSION
1478	FOPOIP(T00,TF,PR0)=DBLE(PR0)*DEXP(-(DLOG(DBLE(T00)/DBLE(TF))/ &
1479	DBLE(CAPPA)))
1480	!
1481	! LES 5 FONCTIONS SUIVANTES SONT VALIDES DANS LE CONTEXTE OU ON
1482	! NE DESIRE PAS TENIR COMPTE DE LA PHASE GLACE DANS LES CALCULS
1483	! DE SATURATION.
1484	! FONCTION DE VAPEUR SATURANTE (TETENS)
1485	FOEWA(TTT)=610.78D0DEXP(17.269D0(DBLE(TTT)-DBLE(TRPL))/ &
1486	(DBLE(TTT)-35.86D0))
1487	! FONCTION CALCULANT LA DERIVEE SELON T DE LN EW
1488	FODLA(TTT)=17.269D0(DBLE(TRPL)-35.86D0)/(DBLE(TTT)-35.86D0)*2
1489	! FONCTION CALCULANT L'HUMIDITE SPECIFIQUE SATURANTE
1490	FOQSA(TTT,PRS)=DBLE(EPS1)/(DMAX1(1.D0,DBLE(PRS)/FOEWA(TTT))- &
1491	DBLE(EPS2))
1492	! FONCTION CALCULANT LA DERIVEE DE QSAT SELON T
1493	FODQA(QST,TTT)=DBLE(QST)(1.D0+DBLE(DELTA)DBLE(QST))*FODLA(TTT)
1494	! FONCTION CALCULANT L'HUMIDITE RELATIVE
1495	#if (DWORDSIZE == 8 && RWORDSIZE == 8)
1496	FOHRA(QQQ,TTT,PRS)=MIN( PRS ,FOEFQ(QQQ,PRS))/FOEWA(TTT)
1497	#elif (DWORDSIZE == 8 && RWORDSIZE == 4)
1498	FOHRA(QQQ,TTT,PRS)=MIN(DBLE(PRS),FOEFQ(QQQ,PRS))/FOEWA(TTT)
1499	#else
1500	This is a temporary hack assuming double precision is 8 bytes.
1501	#endif
1502	!
1503	!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1504	!
1505	! Definition of basic thermodynamic functions in mixed-phase mode
1506	! FFF is the fraction of ice and DDFF its derivative w/r to T
1507	! NOTE: S.I. units are used
1508	! i.e. TTT in deg K, PRS in Pa
1509	! * J. Mailhot - Jan. 2000 *
1510	!
1511	! Saturation calculations in presence of liquid phase only
1512	! Function for saturation vapor pressure (TETENS)
1513	FESI(TTT)=610.78D0DEXP(21.875D0(DBLE(TTT)-DBLE(TRPL))/ &
1514	(DBLE(TTT)-7.66D0) )
1515	FDLESI(TTT)=21.875D0(DBLE(TRPL)-7.66D0)/(DBLE(TTT)-7.66D0)*2
1516	FESMX(TTT,FFF) = (1.D0-DBLE(FFF))FOEWA(TTT)+DBLE(FFF)FESI(TTT)
1517	FDLESMX(TTT,FFF,DDFF) = ( (1.D0-DBLE(FFF))FOEWA(TTT)FODLA(TTT) &
1518	+ DBLE(FFF)FESI(TTT)FDLESI(TTT) &
1519	+ DBLE(DDFF)*(FESI(TTT)-FOEWA(TTT)) )/FESMX(TTT,FFF)
1520	FQSMX(TTT,PRS,FFF) = DBLE(EPS1)/ &
1521	(DMAX1(1.D0,DBLE(PRS)/FESMX(TTT,FFF) ) - DBLE(EPS2) )
1522	FDQSMX(QSM,DLEMX) = DBLE(QSM ) (1.D0 + DBLE(DELTA) DBLE(QSM ) ) &
1523	* DBLE(DLEMX )
1524	!
1525	! ! !------------------------------------------------------------------------------!
1526	!*** END of Replace 3 #includes (for WRF) *
1527
1528	! Constants used for contact ice nucleation:
1529	real, parameter :: LAMa0 = 6.6e-8 ![m] mean free path at T0 and p0 [W95_eqn58]
1530	real, parameter :: T0 = 293.15 ![K] ref. temp.
1531	real, parameter :: p0 = 101325. ![Pa] ref. pres.
1532	real, parameter :: Ra = 1.e-6 ![m] aerosol (IN) radius [M92 p.713; W95_eqn60]
1533	real, parameter :: kBoltz = 1.381e-23 !Boltzmann's constant
1534	real, parameter :: KAPa = 5.39e5 !aerosol thermal conductivity
1535
1536	!Test switches:
1537	logical, parameter :: iceDep_ON = .true. !.false. to suppress depositional growth of ice
1538	logical, parameter :: grpl_ON = .true. !.false. to suppress graupel initiation
1539	logical, parameter :: hail_ON = .true. !.false. to suppress hail initiation
1540	logical, parameter :: rainAccr_ON = .true. ! rain accretion and self-collection ON/OFF
1541	logical, parameter :: snowSpherical = .false. !.true.: m(D)=(pi/6)const_desD^3 \| .false.: m(D)= 0.069*D^2
1542	integer, parameter :: primIceNucl = 1 !1= Meyers+contact ; 2= Cooper
1543	real, parameter :: outfreq = 60. !frequency to compute output diagnostics [s]
1544
1545	!Passed as physics namelist parameters:
1546	! logical, parameter :: precipDiag_ON = .true. !.false. to suppress calc. of sfc precip types
1547	! logical, parameter :: sedi_ON = .true. !.false. to suppress sedimentation
1548	! logical, parameter :: warmphase_ON = .true. !.false. to suppress warm-phase (Part II)
1549	! logical, parameter :: autoconv_ON = .true. ! autoconversion ON/OFF
1550	! logical, parameter :: icephase_ON = .true. !.false. to suppress ice-phase (Part I)
1551	! logical, parameter :: snow_ON = .true. !.false. to suppress snow initiation
1552	! logical, parameter :: initN = .true. !.true. to initialize Nx of Qx>0 and Nx=0
1553
1554	real, dimension(size(QC,dim=1),size(QC,dim=2)) :: DE,iDE,DP,QSS,QSW,QSI,WZ,DZ,RHOQX,FLIM, &
1555	VQQ,gamfact,gamfact_r,massFlux3D_r,massFlux3D_s
1556	real, dimension(size(QC,dim=1)) :: fluxM_r,fluxM_i,fluxM_s,fluxM_g,fluxM_h, &
1557	HPS,dum
1558	integer, dimension(size(QC,dim=1)) :: activeColumn
1559
1560
1561	!==================================================================================!
1562
1563	!----------------------------------------------------------------------------------!
1564	! PART 1: Prelimiary Calculations !
1565	!----------------------------------------------------------------------------------!
1566
1567	!-------------
1568	!Convert N from #/kg to #/m3:
1569	do k= 1,nk
1570	do i= 1,ni
1571	tmp1= S(i,k)PSM(i)/(RGASDTM(i,k)) !air density at time (t-1)
1572	tmp2= S(i,k)PS(i)/(RGASDT(i,k)) !air density at time (*)
1573
1574	NCM(i,k)= NCM(i,k)tmp1; NC(i,k)= NC(i,k)tmp2
1575	NRM(i,k)= NRM(i,k)tmp1; NR(i,k)= NR(i,k)tmp2
1576	NYM(i,k)= NYM(i,k)tmp1; NY(i,k)= NY(i,k)tmp2
1577	NNM(i,k)= NNM(i,k)tmp1; NN(i,k)= NN(i,k)tmp2
1578	NGM(i,k)= NGM(i,k)tmp1; NG(i,k)= NG(i,k)tmp2
1579	NHM(i,k)= NHM(i,k)tmp1; NH(i,k)= NH(i,k)tmp2
1580	enddo
1581	enddo
1582	!=============
1583
1584	! The SSxx arrays are for passed to the volatile bus for output as 3-D diagnostic
1585	! output variables, for testing purposes. For example, to output the
1586	! instantanous value of the deposition rate, add 'SS01(i,k) = QVDvi' in the
1587	! appropriate place. It can then be output as a 3-D physics variable by adding
1588	! it to the sortie_p list in 'outcfgs.out'
1589
1590	SS01= 0.; SS02= 0.; SS03= 0.; SS04= 0.; SS05= 0.; SS06= 0.; SS07= 0.; SS08= 0.
1591	SS09= 0.; SS10= 0.; SS11= 0.; SS12= 0.; SS13= 0.; SS14= 0.; SS15= 0.; SS16= 0.
1592	SS17= 0.; SS18= 0.; SS19= 0.; SS20= 0.
1593
1594	!Determine the upper-most level in each column to which to compute sedimentation:
1595	ktop_sedi= 0
1596	do i=1,ni
1597	do k=1,nk
1598	ktop_sedi(i)= k
1599	if (GZ(i,k)<gzMax_sedi) exit
1600	enddo
1601	enddo
1602
1603	!Compute diagnostic values only every 'outfreq' minutes:
1604	!calcDiag= (mod(DT*float(KOUNT),outfreq)==0.)
1605	calcDiag = .true. !compute diagnostics every step (for time-series output)
1606
1607	!#### These need only to be computed once per model integration:
1608	! (note: These variables must be declared with the SAVE attribute)
1609
1610	! if (KOUNT==0) then
1611	!*** For restarts, these values are not saved. Therefore, the condition statement
1612	! must be modified to something like: IF (MOD(Step_rsti,KOUNT).eq.0) THEN
1613	! in order that these be computed only on the first step of a given restart.
1614	! (...to be done. For now, changing condition to IF(TRUE) to compute at each step.)
1615
1616
1617	if (.TRUE.) then
1618
1619	PI2 = PI*2.
1620	PIov4 = 0.25*PI
1621	PIov6 = PI*sixth
1622	CHLS = CHLC+CHLF !J k-1; latent heat of sublimation
1623	LCP = CHLC/CPD
1624	LFP = CHLF/CPD
1625	iCHLF = 1./CHLF
1626	LSP = LCP+LFP
1627	ck5 = 4098.170*LCP
1628	ck6 = 5806.485*LSP
1629	idt = 1./dt
1630	imgo = 1./mgo
1631	idew = 1./dew
1632	idei = 1./dei
1633	ideg = 1./deg
1634	ideh = 1./deh
1635
1636	!Constants based on size distribution parameters:
1637
1638	! Mass parameters [ m(D) = cD^d ]
1639	cmr = PIov6*dew; icmr= 1./cmr
1640	cmi = 440.; icmi= 1./cmi
1641	cmg = PIov6*deg; icmg= 1./cmg
1642	cmh = PIov6*deh; icmh= 1./cmh
1643
1644	cms_D3 = PIov6*desFix !used for snowSpherical = .T. or .F.
1645	if (snowSpherical) then
1646	cms = cms_D3
1647	dms = 3.
1648	else
1649	! cms = 0.0690; dms = 2.000 !Cox, 1988 (QJRMS)
1650	cms = 0.1597; dms = 2.078 !Brandes et al., 2007 (JAMC)
1651	endif
1652	icms = 1./cms
1653	idms = 1./dms
1654	mso = cmsDso*dms
1655	imso = 1./mso
1656	!bulk density parameters: [rho(D) = eds*D^fds]
1657	! These are implied by the mass-diameter parameters, by computing the bulk
1658	! density of a sphere with the equaivalent mass.
1659	! e.g. m(D) = cD^d = (pi/6)rhoD^3 and solve for rho(D)
1660	eds = cms/PIov6
1661	fds = dms-3.
1662	if (fds/=-1. .and..not.snowSpherical) GS50= gamma(1.+fds+alpha_s)
1663
1664	! Cloud:
1665	iMUc = 1./MUc
1666	GC1 = gamma(alpha_c+1.0)
1667	iGC1 = 1./GC1
1668	GC2 = gamma(alpha_c+1.+3.0*iMUc) !i.e. gamma(alf + 4)
1669	GC3 = gamma(alpha_c+1.+6.0*iMUc) !i.e. gamma(alf + 7)
1670	GC4 = gamma(alpha_c+1.+9.0*iMUc) !i.e. gamma(alf + 10)
1671	GC11 = gamma(1.0*iMUc+1.0+alpha_c)
1672	GC12 = gamma(2.0*iMUc+1.0+alpha_c)
1673	GC5 = gamma(1.0+alpha_c)
1674	iGC5 = 1./GC5
1675	GC6 = gamma(1.0+alpha_c+1.0*iMUc)
1676	GC7 = gamma(1.0+alpha_c+2.0*iMUc)
1677	GC8 = gamma(1.0+alpha_c+3.0*iMUc)
1678	GC13 = gamma(3.0*iMUc+1.0+alpha_c)
1679	GC14 = gamma(4.0*iMUc+1.0+alpha_c)
1680	GC15 = gamma(5.0*iMUc+1.0+alpha_c)
1681	icexc9 = 1./(GC2iGC1PIov6*dew)
1682	!specify cloud droplet number concentration [m-3] based on 'CCNtype' (1-moment):
1683	if (CCNtype==1) then
1684	N_c_SM = 0.8e+8 !maritime
1685	elseif (CCNtype==2) then
1686	N_c_SM = 2.0e+8 !continental 1
1687	elseif (CCNtype==3) then
1688	N_c_SM = 5.0e+8 !continental 2 (polluted)
1689	else
1690	N_c_SM = 2.0e+8 !default (cont1), if 'CCNtype' specified incorrectly
1691	endif
1692
1693	! Rain:
1694	cexr1 = 1.+alpha_r+dmr+bfr
1695	cexr2 = 1.+alpha_r+dmr
1696	GR17 = gamma(2.5+alpha_r+0.5*bfr)
1697	GR31 = gamma(1.+alpha_r)
1698	iGR31 = 1./GR31
1699	GR32 = gamma(2.+alpha_r)
1700	GR33 = gamma(3.+alpha_r)
1701	GR34 = gamma(4.+alpha_r)
1702	iGR34 = 1./GR34
1703	GR35 = gamma(5.+alpha_r)
1704	GR36 = gamma(6.+alpha_r)
1705	GR37 = gamma(7.+alpha_r)
1706	GR50 = (No_r_SMGR31)*0.75 !for 1-moment or Nr-initialization
1707	cexr5 = 2.+alpha_r
1708	cexr6 = 2.5+alpha_r+0.5*bfr
1709	cexr9 = cmrGR34iGR31; icexr9= 1./cexr9
1710	cexr3 = 1.+bfr+alpha_r
1711	cexr4 = 1.+alpha_r
1712	ckQr1 = afr*gamma(1.+alpha_r+dmr+bfr)/gamma(1.+alpha_r+dmr)
1713	ckQr2 = afrgamma(1.+alpha_r+bfr)GR31
1714	ckQr3 = afr*gamma(7.+alpha_r+bfr)/GR37
1715	if (.not.dblMom_r) then
1716	No_r = No_r_SM
1717	endif
1718
1719	! Ice:
1720	GI4 = gamma(alpha_i+dmi+bfi)
1721	GI6 = gamma(2.5+bfi*0.5+alpha_i)
1722	GI11 = gamma(1.+bfi+alpha_i)
1723	GI20 = gamma(0.+bfi+1.+alpha_i)
1724	GI21 = gamma(1.+bfi+1.+alpha_i)
1725	GI22 = gamma(2.+bfi+1.+alpha_i)
1726	GI31 = gamma(1.+alpha_i)
1727	iGI31 = 1./GI31
1728	GI32 = gamma(2.+alpha_i)
1729	GI33 = gamma(3.+alpha_i)
1730	GI34 = gamma(4.+alpha_i)
1731	GI35 = gamma(5.+alpha_i)
1732	GI36 = gamma(6.+alpha_i)
1733	GI40 = gamma(1.+alpha_i+dmi)
1734	icexi9 = 1./(cmigamma(1.+alpha_i+dmi)iGI31)
1735	ckQi1 = afi*gamma(1.+alpha_i+dmi+bfi)/GI40
1736	ckQi2 = afiGI11iGI31
1737	ckQi4 = 1./(cmiGI40iGI31)
1738
1739	! Snow:
1740	cexs1 = 2.5+0.5*bfs+alpha_s
1741	cexs2 = 1.+alpha_s+dms
1742	icexs2 = 1./cexs2
1743	GS09 = gamma(2.5+bfs*0.5+alpha_s)
1744	GS11 = gamma(1.+bfs+alpha_s)
1745	GS12 = gamma(2.+bfs+alpha_s)
1746	GS13 = gamma(3.+bfs+alpha_s)
1747	GS31 = gamma(1.+alpha_s)
1748	iGS31 = 1./GS31
1749	GS32 = gamma(2.+alpha_s)
1750	GS33 = gamma(3.+alpha_s)
1751	GS34 = gamma(4.+alpha_s)
1752	iGS34 = 1./GS34
1753	GS35 = gamma(5.+alpha_s)
1754	GS36 = gamma(6.+alpha_s)
1755	GS40 = gamma(1.+alpha_s+dms)
1756	iGS40 = 1./GS40
1757	iGS20 = 1./(GS40iGS31cms)
1758	ckQs1 = afsgamma(1.+alpha_s+dms+bfs)iGS40
1759	ckQs2 = afsGS11iGS31
1760	GS40_D3 = gamma(1.+alpha_s+3.)
1761	iGS20_D3= 1./(GS40_D3iGS31cms_D3)
1762	rfact_FvFm= PIov6icmsgamma(4.+bfs+alpha_s)/gamma(1.+dms+bfs+alpha_s)
1763
1764	! Graupel:
1765	GG09 = gamma(2.5+0.5*bfg+alpha_g)
1766	GG11 = gamma(1.+bfg+alpha_g)
1767	GG12 = gamma(2.+bfg+alpha_g)
1768	GG13 = gamma(3.+bfg+alpha_g)
1769	GG31 = gamma(1.+alpha_g)
1770	iGG31 = 1./GG31
1771	GG32 = gamma(2.+alpha_g)
1772	GG33 = gamma(3.+alpha_g)
1773	GG34 = gamma(4.+alpha_g)
1774	iGG34 = 1./GG34
1775	GG35 = gamma(5.+alpha_g)
1776	GG36 = gamma(6.+alpha_g)
1777	GG40 = gamma(1.+alpha_g+dmg)
1778	iGG99 = 1./(GG40iGG31cmg)
1779	GG50 = (No_g_SMGG31)*0.75 !for 1-moment only
1780	ckQg1 = afg*gamma(1.+alpha_g+dmg+bfg)/GG40
1781	ckQg2 = afgGG11iGG31
1782	ckQg4 = 1./(cmgGG40iGG31)
1783
1784	! Hail:
1785	GH09 = gamma(2.5+bfh*0.5+alpha_h)
1786	GH11 = gamma(1.+bfh+alpha_h)
1787	GH12 = gamma(2.+bfh+alpha_h)
1788	GH13 = gamma(3.+bfh+alpha_h)
1789	GH31 = gamma(1.+alpha_h)
1790	iGH31 = 1./GH31
1791	GH32 = gamma(2.+alpha_h)
1792	GH33 = gamma(3.+alpha_h)
1793	iGH34 = 1./gamma(4.+alpha_h)
1794	GH40 = gamma(1.+alpha_h+dmh)
1795	iGH99 = 1./(GH40iGH31cmh)
1796	GH50 = (No_h_SMGH31)*0.75 !for 1-moment only
1797	ckQh1 = afh*gamma(1.+alpha_h+dmh+bfh)/GH40
1798	ckQh2 = afhGH11iGH31
1799	ckQh4 = 1./(cmhGH40iGH31)
1800
1801	endif !if (KOUNT=0)
1802	!####
1803
1804	!=======================================================================================!
1805
1806	!Compute thickness of layers for sedimentation calcuation:
1807	! (note; 'GZ' passed in is geopotential, not geopotential height)
1808	tmp1= 1./GRAV
1809	do k=2,nk
1810	DZ(:,k)= (GZ(:,k-1)-GZ(:,k))*tmp1
1811	enddo
1812	DZ(:,1)= DZ(:,2)
1813
1814	! Temporarily store arrays at time (t*) in order to compute (at the end of subroutine)
1815	! the final VXTEND as VXTEND = ( VX{t+1} - VX{t*} )/dt :
1816	T_TEND = T ; Q_TEND = Q
1817	QCTEND = QC; QRTEND = QR; QITEND = QI; QNTEND = QN; QGTEND = QG; QHTEND = QH
1818	NCTEND = NC; NRTEND = NR; NYTEND = NY; NNTEND = NN; NGTEND = NG; NHTEND = NH
1819
1820	!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -!
1821	! Initialize Nx if Qx>0 and Nx=0: (for nesting from 1-moment to 2-moment):
1822	IF (initN) THEN
1823	do k= 1,nk
1824	do i= 1,ni
1825	tmp1= S(i,k)PSM(i)/(RGASDTM(i,k)) !air density at time (t-1)
1826	tmp2= S(i,k)PS(i)/(RGASDT(i,k)) !air density at time (*)
1827
1828	!cloud:
1829	if (QCM(i,k)>epsQ .and. NCM(i,k)<epsN) &
1830	NCM(i,k)= N_c_SM
1831	if (QC(i,k)>epsQ .and. NC(i,k)<epsN) &
1832	NC(i,k) = N_c_SM
1833	!rain
1834	if (QRM(i,k)>epsQ .and. NRM(i,k)<epsN) &
1835	NRM(i,k)= (No_r_SMGR31)(3./(4.+alpha_r))(GR31iGR34tmp1QRM(i,k) &
1836	icmr)**((1.+alpha_r)/(4.+alpha_r))
1837	if (QR(i,k)>epsQ .and. NR(i,k)<epsN) &
1838	NR(i,k)= (No_r_SMGR31)(3./(4.+alpha_r))(GR31iGR34tmp2QR(i,k) &
1839	icmr)**((1.+alpha_r)/(4.+alpha_r))
1840	!ice:
1841	if (QIM(i,k)>epsQ .and. NYM(i,k)<epsN) &
1842	NYM(i,k)= N_Cooper(TRPL,TM(i,k))
1843	if (QI(i,k)>epsQ .and. NY(i,k)<epsN) &
1844	NY(i,k)= N_Cooper(TRPL,T(i,k))
1845	!snow:
1846	if (QNM(i,k)>epsQ .and. NNM(i,k)<epsN) then
1847	No_s= Nos_Thompson(TRPL,TM(i,k))
1848	NNM(i,k)= (No_sGS31)(dmsicexs2)(GS31iGS40icmstmp1QNM(i,k))* &
1849	((1.+alpha_s)*icexs2)
1850	endif
1851	if (QN(i,k)>epsQ .and. NN(i,k)<epsN) then
1852	No_s= Nos_Thompson(TRPL,T(i,k))
1853	NN(i,k)= (No_sGS31)(dmsicexs2)(GS31iGS40icmstmp2QN(i,k))* &
1854	((1.+alpha_s)*icexs2)
1855	endif
1856	!grpl:
1857	if (QGM(i,k)>epsQ .and. NGM(i,k)<epsN) &
1858	NGM(i,k)= (No_g_SMGG31)(3./(4.+alpha_g))(GG31iGG34tmp1QGM(i,k) &
1859	icmg)**((1.+alpha_g)/(4.+alpha_g))
1860	if (QG(i,k)>epsQ .and. NG(i,k)<epsN) &
1861	NG(i,k)= (No_g_SMGG31)(3./(4.+alpha_g))(GG31iGG34tmp2QG(i,k) &
1862	icmg)**((1.+alpha_g)/(4.+alpha_g))
1863	!hail:
1864	if (QHM(i,k)>epsQ .and. NHM(i,k)<epsN) &
1865	NHM(i,k)= (No_h_SMGH31)(3./(4.+alpha_h))(GH31iGH34tmp1QHM(i,k) &
1866	icmh)**((1.+alpha_h)/(4.+alpha_h))
1867	if (QH(i,k)>epsQ .and. NH(i,k)<epsN) &
1868	NH(i,k)= (No_h_SMGH31)(3./(4.+alpha_h))(GH31iGH34tmp2QH(i,k) &
1869	icmh)**((1.+alpha_h)/(4.+alpha_h))
1870
1871	enddo !i-loop
1872	enddo !k-loop
1873	ENDIF !N-initialization
1874
1875	!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -!
1876	! Clip all moments to zero if one or more corresponding category moments are less than
1877	! the minimum allowable value:
1878	! (Note: Clipped mass is added back to water vapor field to conserve total mass)
1879	do k= 1,nk
1880	do i= 1,ni
1881
1882	IF (dblMom_c) THEN
1883	if(QC(i,k)<epsQ .or. NC(i,k)<epsN) then
1884	Q(i,k) = Q(i,k) + QC(i,k)
1885	QC(i,k)= 0.; NC(i,k)= 0.
1886	endif
1887	if(QCM(i,k)<epsQ .or. NCM(i,k)<epsN) then
1888	QM(i,k) = QM(i,k) + QCM(i,k)
1889	QCM(i,k)= 0.; NCM(i,k)= 0.
1890	endif
1891	ELSE
1892	if(QC(i,k)<epsQ) then
1893	Q(i,k) = Q(i,k) + QC(i,k)
1894	QC(i,k)= 0.
1895	endif
1896	if(QCM(i,k)<epsQ) then
1897	QM(i,k) = QM(i,k) + QCM(i,k)
1898	QCM(i,k)= 0.
1899	endif
1900	ENDIF
1901
1902	IF (dblMom_r) THEN
1903	if (QR(i,k)<epsQ .or. NR(i,k)<epsN) then
1904	Q(i,k) = Q(i,k) + QR(i,k)
1905	QR(i,k)= 0.; NR(i,k)= 0.
1906	endif
1907	if (QRM(i,k)<epsQ .or. NRM(i,k)<epsN) then
1908	QM(i,k) = QM(i,k) + QRM(i,k)
1909	QRM(i,k)= 0.; NRM(i,k)= 0.
1910	endif
1911	ELSE
1912	if (QR(i,k)<epsQ) then
1913	Q(i,k) = Q(i,k) + QR(i,k)
1914	QR(i,k)= 0.
1915	endif
1916	if (QRM(i,k)<epsQ) then
1917	QM(i,k) = QM(i,k) + QRM(i,k)
1918	QRM(i,k)= 0.
1919	endif
1920	ENDIF
1921
1922	IF (dblMom_i) THEN
1923	if (QI(i,k)<epsQ .or. NY(i,k)<epsN) then
1924	Q(i,k) = Q(i,k) + QI(i,k)
1925	QI(i,k)= 0.; NY(i,k)= 0.
1926	endif
1927	if (QIM(i,k)<epsQ .or. NYM(i,k)<epsN) then
1928	QM(i,k) = QM(i,k) + QIM(i,k)
1929	QIM(i,k)= 0.; NYM(i,k)= 0.
1930	endif
1931	ELSE
1932	if (QI(i,k)<epsQ) then
1933	Q(i,k) = Q(i,k) + QI(i,k)
1934	QI(i,k)= 0.
1935	endif
1936	if (QIM(i,k)<epsQ) then
1937	QM(i,k) = QM(i,k) + QIM(i,k)
1938	QIM(i,k)= 0.
1939	endif
1940	ENDIF
1941
1942	IF (dblMom_s) THEN
1943	if (QN(i,k)<epsQ .or. NN(i,k)<epsN) then
1944	Q(i,k) = Q(i,k) + QN(i,k)
1945	QN(i,k)= 0.; NN(i,k)= 0.
1946	endif
1947	if (QNM(i,k)<epsQ .or. NNM(i,k)<epsN) then
1948	QM(i,k) = QM(i,k) + QNM(i,k)
1949	QNM(i,k)= 0.; NNM(i,k)= 0.
1950	endif
1951	ELSE
1952	if (QN(i,k)<epsQ) then
1953	Q(i,k) = Q(i,k) + QN(i,k)
1954	QN(i,k)= 0.
1955	endif
1956	if (QNM(i,k)<epsQ) then
1957	QM(i,k) = QM(i,k) + QNM(i,k)
1958	QNM(i,k)= 0.
1959	endif
1960	ENDIF
1961
1962	IF (dblMom_g) THEN
1963	if (QG(i,k)<epsQ .or. NG(i,k)<epsN) then
1964	Q(i,k) = Q(i,k) + QG(i,k)
1965	QG(i,k)= 0.; NG(i,k)= 0.
1966	endif
1967	if (QGM(i,k)<epsQ .or. NGM(i,k)<epsN) then
1968	QM(i,k) = QM(i,k) + QGM(i,k)
1969	QGM(i,k)= 0.; NGM(i,k)= 0.
1970	endif
1971	ELSE
1972	if (QG(i,k)<epsQ) then
1973	Q(i,k) = Q(i,k) + QG(i,k)
1974	QG(i,k)= 0.
1975	endif
1976	if (QGM(i,k)<epsQ) then
1977	QM(i,k) = QM(i,k) + QGM(i,k)
1978	QGM(i,k)= 0.
1979	endif
1980	ENDIF
1981
1982	IF (dblMom_h) THEN
1983	if (QH(i,k)<epsQ .or. NH(i,k)<epsN) then
1984	Q(i,k) = Q(i,k) + QH(i,k)
1985	QH(i,k)= 0.; NH(i,k)= 0.
1986	endif
1987	if (QHM(i,k)<epsQ .or. NHM(i,k)<epsN) then
1988	QM(i,k) = QM(i,k) + QHM(i,k)
1989	QHM(i,k)= 0.; NHM(i,k)= 0.
1990	endif
1991	ELSE
1992	if (QH(i,k)<epsQ) then
1993	Q(i,k) = Q(i,k) + QH(i,k)
1994	QH(i,k)= 0.
1995	endif
1996	if (QHM(i,k)<epsQ) then
1997	QM(i,k) = QM(i,k) + QHM(i,k)
1998	QHM(i,k)= 0.
1999	endif
2000	ENDIF
2001
2002	enddo !i-loop
2003	enddo !k-loop; (clipping)
2004	QM = max(QM,0.)
2005	Q = max(Q ,0.)
2006
2007	!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -!
2008
2009	! Approximate values at time {t}:
2010	! [ ave. of values at {*} (advected, but no physics tendency added) and {t-dt} ]:
2011	HPS= 0.5(PSM+PS); TM = 0.5(TM + T); QM = 0.5*(QM + Q)
2012	QCM= 0.5(QCM+QC); QRM= 0.5(QRM+QR); QIM= 0.5*(QIM+QI)
2013	QNM= 0.5(QNM+QN); QGM= 0.5(QGM+QG); QHM= 0.5*(QHM+QH)
2014
2015	if (dblMom_c) NCM= 0.5*(NCM+NC)
2016	if (dblMom_r) NRM= 0.5*(NRM+NR)
2017	if (dblMom_i) NYM= 0.5*(NYM+NY)
2018	if (dblMom_s) NNM= 0.5*(NNM+NN)
2019	if (dblMom_g) NGM= 0.5*(NGM+NG)
2020	if (dblMom_h) NHM= 0.5*(NHM+NH)
2021
2022	do k=1,nk
2023	do i=1,ni
2024	!WRF:
2025	#if (DWORDSIZE == 8 && RWORDSIZE == 8)
2026	QSW(i,k)= FOQSA(TM(i,k),HPS(i)*S(i,k)) !wrt. liquid water at (t)
2027	QSS(i,k)= FOQST( T(i,k), PS(i)S(i,k)) !wrt. ice surface at ()
2028	QSI(i,k)= FOQST(TM(i,k),HPS(i)*S(i,k)) !wrt. ice surface at (t)
2029	#elif (DWORDSIZE == 8 && RWORDSIZE == 4)
2030	QSW(i,k)= sngl(FOQSA(TM(i,k),HPS(i)*S(i,k))) !wrt. liquid water at (t)
2031	QSS(i,k)= sngl(FOQST( T(i,k), PS(i)S(i,k))) !wrt. ice surface at ()
2032	QSI(i,k)= sngl(FOQST(TM(i,k),HPS(i)*S(i,k))) !wrt. ice surface at (t)
2033	#else
2034	!! This is a temporary hack assuming double precision is 8 bytes.
2035	#endif
2036	!Air density at time (t)
2037	DE(i,k) = S(i,k)HPS(i)/(RGASDTM(i,k)) !air density at time (t)
2038	iDE(i,k)= 1./DE(i,k)
2039	enddo
2040	enddo
2041
2042	do i= 1,ni
2043	!Air-density factor: (for fall velocity computations)
2044	DEo = DE(i,nk)
2045	gamfact(i,:) = sqrt(DEo/(DE(i,:)))
2046	gamfact_r(i,:)= sqrt( 1./(DE(i,:)))
2047	!Convert 'W_omega' (on thermodynamic levels) to 'w' (on momentum):
2048	do k= 2,nk-1
2049	WZ(i,k)= -0.5/(DE(i,k)GRAV)(W_omega(i,k-1)+W_omega(i,k+1))
2050	enddo
2051	WZ(i,1) = -0.5/(DE(i,1) GRAV)W_omega(i,1)
2052	WZ(i,nk)= -0.5/(DE(i,nk)GRAV)W_omega(i,nk)
2053	enddo
2054
2055	!----------------------------------------------------------------------------------!
2056	! End of Preliminary Calculation section (Part 1) !
2057	!----------------------------------------------------------------------------------!
2058
2059	!----------------------------------------------------------------------------------!
2060	! PART 2: Cold Microphysics Processes !
2061	!----------------------------------------------------------------------------------!
2062
2063	! Determine the active grid points (i.e. those which scheme should treat):
2064	activePoint = .false.
2065	DO k=2,nk
2066	DO i=1,ni
2067	log1= ((QIM(i,k)+QGM(i,k)+QNM(i,k)+QHM(i,k))<epsQ) !no solid (i,g,s,h)
2068	log2= ((QCM(i,k)+QRM(i,k)) <epsQ) !no liquid (c,r)
2069	log3= ((TM(i,k)>TRPL) .and. log1) !T>0C & no i,g,s,h
2070	log4= log1.and.log2.and.(QM(i,k)<QSI(i,k)) !no sol. or liq.; subsat(i)
2071	if (.not.( log3 .or. log4 ) .and. icephase_ON) then
2072	activePoint(i,k)= .true.
2073	endif
2074	ENDDO
2075	ENDDO
2076
2077	! Size distribution parameters:
2078	! Note: + 'thrd' should actually be '1/dmx'(but dmx=3 for all categories x)
2079	! + If Qx=0, LAMx etc. are never be used in any calculations
2080	! (If Qc=0, CLcy etc. will never be calculated. iLAMx is set to 0
2081	! to avoid possible problems due to bugs.)
2082
2083	DO k= 2,nk !Main loop for Part 2
2084	DO i= 1,ni
2085	IF (activePoint(i,k)) THEN
2086
2087
2088	Tc= TM(i,k)-TRPL
2089	if (Tc<-120. .or. Tc>50.) &
2090	print, 'WARNING* -- In MICROPHYSICS -- Ambient Temp.(C):',Tc
2091	Cdiff = (2.2157e-5+0.0155e-5Tc)1.e5/(S(i,k)*HPS(i))
2092	MUdyn = 1.72e-5(393./(TM(i,k)+120.))(TM(i,k)/TRPL)**1.5 !RYp.102
2093	MUkin = MUdyn*iDE(i,k)
2094	iMUkin= 1./MUkin
2095	ScTHRD= (MUkin/Cdiff)**thrd ! i.e. Sc^(1/3)
2096	Ka = 2.3971e-2 + 0.0078e-2*Tc !therm.cond.(air)
2097	Kdiff = (9.1018e-11TM(i,k)TM(i,k)+8.8197e-8*TM(i,k)-(1.0654e-5)) !therm.diff.(air)
2098	gam = gamfact(i,k)
2099
2100	!Collection efficiencies:
2101	Eis = min(0.05exp(0.1Tc),1.) !Ferrier, 1995 (Table 1)
2102	Eig = min(0.01exp(0.1Tc),1.) !dry (Eig=1.0 for wet growth)
2103	Eii = 0.1*Eis
2104	Ess = Eis; Eih = Eig; Esh = Eig
2105	iEih = 1./Eih
2106	iEsh = 1./Esh
2107	!note: Eri=Ers=Erh=1. (constant parameters)
2108	! - Ecs is computed in CLcs section
2109	! - Ech is computed in CLch section
2110	! - Ecg is computed in CLcg section
2111	! - Erg is computed in CLrg section
2112
2113	!WRF:
2114	#if (DWORDSIZE == 8 && RWORDSIZE == 8)
2115	qvs0 = FOQSA(TRPL,HPS(i)*S(i,k)) !sat.mix.ratio at 0C
2116	#elif (DWORDSIZE == 8 && RWORDSIZE == 4)
2117	qvs0 = sngl(FOQSA(TRPL,HPS(i)*S(i,k))) !sat.mix.ratio at 0C
2118	#else
2119	!! This is a temporary hack assuming double precision is 8 bytes.
2120	#endif
2121	DELqvs= qvs0-(QM(i,k))
2122
2123	! Cloud:
2124	if (QCM(i,k)>epsQ) then
2125	if (.not. dblMom_c) NCM(i,k)= N_c_SM
2126	iQCM = 1./QCM(i,k)
2127	iNCM = 1./NCM(i,k)
2128	Dc = Dm_x(DE(i,k),QCM(i,k),iNCM,icmr,thrd)
2129
2130	iLAMc = iLAMDA_x(DE(i,k),QCM(i,k),iNCM,icexc9,thrd)
2131	iLAMc2 = iLAMc *iLAMc
2132	iLAMc3 = iLAMc2*iLAMc
2133	iLAMc4 = iLAMc2*iLAMc2
2134	iLAMc5 = iLAMc3*iLAMc2
2135	else
2136	Dc = 0.; iLAMc3= 0.
2137	iLAMc = 0.; iLAMc4= 0.
2138	iLAMc2 = 0.; iLAMc5= 0.
2139	endif
2140
2141	! Rain:
2142	if (QRM(i,k)>epsQ) then
2143	if (.not. dblMom_r) NRM(i,k)= GR50sqrt(sqrt(GR31iGR34DE(i,k)QRM(i,k)*icmr))
2144	iQRM = 1./QRM(i,k)
2145	iNRM = 1./NRM(i,k)
2146	Dr = Dm_x(DE(i,k),QRM(i,k),iNRM,icmr,thrd)
2147	iLAMr = max( iLAMmin1, iLAMDA_x(DE(i,k),QRM(i,k),iNRM,icexr9,thrd) )
2148	tmp1 = 1./iLAMr
2149	iLAMr2 = iLAMr *iLAMr
2150	iLAMr3 = iLAMr2*iLAMr
2151	iLAMr4 = iLAMr2*iLAMr2
2152	iLAMr5 = iLAMr3*iLAMr2
2153	if (Dr>40.e-6) then
2154	vr0 = gamfact_r(i,k)ckQr1iLAMr**bfr
2155	else
2156	vr0 = 0.
2157	endif
2158	else
2159	iLAMr = 0.; Dr = 0.; vr0 = 0.
2160	iLAMr2 = 0.; iLAMr3= 0.; iLAMr4= 0.; iLAMr5 = 0.
2161	endif
2162
2163	! Ice:
2164	if (QIM(i,k)>epsQ) then
2165	if (.not. dblMom_i) NYM(i,k)= N_Cooper(TRPL,TM(i,k))
2166
2167	iQIM = 1./QIM(i,k)
2168	iNYM = 1./NYM(i,k)
2169	iLAMi = max( iLAMmin2, iLAMDA_x(DE(i,k),QIM(i,k),iNYM,icexi9,thrd) )
2170	iLAMi2 = iLAMi *iLAMi
2171	iLAMi3 = iLAMi2*iLAMi
2172	iLAMi4 = iLAMi2*iLAMi2
2173	iLAMi5 = iLAMi3*iLAMi2
2174	iLAMiB0= iLAMi**(bfi)
2175	iLAMiB1= iLAMi**(bfi+1.)
2176	iLAMiB2= iLAMi**(bfi+2.)
2177	vi0 = gamfact(i,k)ckQi1iLAMiB0
2178	Di = Dm_x(DE(i,k),QIM(i,k),iNYM,icmi,thrd)
2179	else
2180	iLAMi = 0.; vi0 = 0.; Di = 0.
2181	iLAMi2 = 0.; iLAMi3 = 0.; iLAMi4 = 0.; iLAMi5= 0.
2182	iLAMiB0= 0.; iLAMiB1= 0.; iLAMiB2= 0.
2183	endif
2184
2185	! Snow:
2186	if (QNM(i,k)>epsQ) then
2187	if (.not.dblMom_s) then
2188	No_s_SM = Nos_Thompson(TRPL,TM(i,k))
2189	NNM(i,k)= (No_sGS31)(dmsicexs2)(GS31iGS40icmsDE(i,k)QNM(i,k))* &
2190	((1.+alpha_s)*icexs2)
2191	endif
2192	iQNM = 1./QNM(i,k)
2193	iNNM = 1./NNM(i,k)
2194	iLAMs = max( iLAMmin2, iLAMDA_x(DE(i,k),QNM(i,k),iNNM,iGS20,idms) )
2195	iLAMs_D3= max(iLAMmin2, iLAMDA_x(DE(i,k),QNM(i,k),iNNM,iGS20_D3,thrd) )
2196	iLAMs2 = iLAMs*iLAMs
2197	iLAMsB0= iLAMs**(bfs)
2198	iLAMsB1= iLAMs**(bfs+1.)
2199	iLAMsB2= iLAMs**(bfs+2.)
2200	vs0 = gamfact(i,k)ckQs1iLAMsB0
2201	Ds = min(DsMax, Dm_x(DE(i,k),QNM(i,k),iNNM,icms,idms))
2202	if (snowSpherical) then
2203	des = desFix
2204	else
2205	des = des_OF_Ds(Ds,desMax,eds,fds)
2206	endif
2207	!!-- generalized equations (any alpha_s):
2208	! No_s = (NNM(i,k))iGS31/iLAMs*(1.+alpha_s)
2209	! VENTs = AvxGS32iLAMs*(2.+alpha_s)+BvxScTHRDsqrt(gamafsiMUkin) &
2210	!!-- GS09iLAMs(2.5+0.5bfs+alpha_s)
2211	!The following equations for No_s and VENTs is based on m(D)=(pi/6)100.D**3 for snow.
2212	! Strict application of m(D)=cD*2 would require re-derivation using implied
2213	! definition of D as the MAXIMUM DIMENSION of an ellipsoid, rather than a sphere.
2214	! For simplicity, the m-D^3 relation is applied -- used for VDvs and MLsr only.
2215	if (dblMom_s) then
2216	!No_s= NNM(i,k)*iGS31/iLAMs !optimized for alpha_s=0
2217	No_s= NNM(i,k)*iGS31/iLAMs_D3 !based on m-D^3 (consistent with VENTs, below)
2218	else
2219	No_s= No_s_SM
2220	endif
2221	VENTs= AvxGS32iLAMs_D3*2. + BvxScTHRDsqrt(gamfact(i,k)afsiMUkin)GS09* &
2222	iLAMs_D3**cexs1
2223	else
2224	iLAMs = 0.; vs0 = 0.; Ds = 0.; iLAMs2= 0.
2225	iLAMsB0= 0.; iLAMsB1= 0.; iLAMsB1= 0.
2226	des = desFix !used for 3-component freezing if QNM=0 (even for snowSpherical=.F.)
2227	endif
2228	ides = 1./des
2229
2230
2231	! Graupel:
2232	if (QGM(i,k)>epsQ) then
2233	if (.not.dblMom_g) NGM(i,k)= GG50sqrt(sqrt(GG31GG34DE(i,k)QGM(i,k)*icmg))
2234	iQGM = 1./QGM(i,k)
2235	iNGM = 1./NGM(i,k)
2236	iLAMg = max( iLAMmin1, iLAMDA_x(DE(i,k),QGM(i,k),iNGM,iGG99,thrd) )
2237	iLAMg2 = iLAMg *iLAMg
2238	iLAMgB0= iLAMg**(bfg)
2239	iLAMgB1= iLAMg**(bfg+1.)
2240	iLAMgB2= iLAMg**(bfg+2.)
2241	if (dblMom_g) then
2242	!No_g = (NGM(i,k))iGG31/iLAMg*(1.+alpha_g)
2243	No_g= NGM(i,k)*iGG31/iLAMg !optimized for alpha_g=0
2244	else
2245	No_g= No_g_SM
2246	endif
2247	vg0 = gamfact(i,k)ckQg1iLAMgB0
2248	Dg = Dm_x(DE(i,k),QGM(i,k),iNGM,icmg,thrd)
2249	else
2250	iLAMg = 0.; vg0 = 0.; Dg = 0.; No_g = 0.
2251	iLAMg2 = 0.; iLAMgB0= 0.; iLAMgB1= 0.; iLAMgB1= 0.
2252	endif
2253
2254	! Hail:
2255	if (QHM(i,k)>epsQ) then
2256	if (.not.dblMom_h) NHM(i,k)= GH50sqrt(sqrt(GH31iGH34DE(i,k)QHM(i,k)*icmh))
2257	iQHM = 1./QHM(i,k)
2258	iNHM = 1./NHM(i,k)
2259	iLAMh = max( iLAMmin1, iLAMDA_x(DE(i,k),QHM(i,k),iNHM,iGH99,thrd) )
2260	iLAMh2 = iLAMh*iLAMh
2261	iLAMhB0= iLAMh**(bfh)
2262	iLAMhB1= iLAMh**(bfh+1.)
2263	iLAMhB2= iLAMh**(bfh+2.)
2264	if (dblMom_h) then
2265	No_h= NHM(i,k)iGH31/iLAMh*(1.+alpha_h)
2266	else
2267	No_h= No_h_SM
2268	endif
2269	vh0 = gamfact(i,k)ckQh1iLAMhB0
2270	Dh = Dm_x(DE(i,k),QHM(i,k),iNHM,icmh,thrd)
2271	else
2272	iLAMh = 0.; vh0 = 0.; Dh = 0.; No_h= 0.
2273	iLAMhB0= 0.; iLAMhB1= 0.; iLAMhB1= 0.
2274	endif
2275	!------
2276
2277	!Calculating ice-phase source/sink terms:
2278
2279	! Initialize all source terms to zero:
2280	QNUvi=0.; QVDvi=0.; QVDvs=0.; QVDvg=0.; QVDvh=0.
2281	QCLcs=0.; QCLcg=0.; QCLch=0.; QFZci=0.; QCLri=0.; QMLsr=0.
2282	QCLrs=0.; QCLrg=0.; QMLgr=0.; QCLrh=0.; QMLhr=0.; QFZrh=0.
2283	QMLir=0.; QCLsr=0.; QCLsh=0.; QCLgr=0.; QCNgh=0.
2284	QCNis=0.; QCLir=0.; QCLis=0.; QCLih=0.
2285	QIMsi=0.; QIMgi=0.; QCNsg=0.; QHwet=0.
2286
2287	NCLcs= 0.; NCLcg=0.; NCLch=0.; NFZci=0.; NMLhr=0.; NhCNgh=0.
2288	NCLri= 0.; NCLrs=0.; NCLrg=0.; NCLrh=0.; NMLsr=0.; NMLgr=0.
2289	NMLir= 0.; NSHhr=0.; NNUvi=0.; NVDvi=0.; NVDvh=0.; QCLig=0.
2290	NCLir= 0.; NCLis=0.; NCLig=0.; NCLih=0.; NIMsi=0.; NIMgi=0.
2291	NiCNis=0.; NsCNis=0.; NVDvs=0.; NCNsg=0.; NCLgr=0.; NCLsrh=0.
2292	NCLss= 0.; NCLsr=0.; NCLsh=0.; NCLsrs=0.; NCLgrg=0.; NgCNgh=0.
2293	NVDvg= 0.; NCLirg=0.; NCLsrg=0.; NCLgrh=0.; NrFZrh=0.; NhFZrh=0.
2294	NCLirh=0.
2295
2296	Dirg=0.; Dirh=0.; Dsrs= 0.; Dsrg= 0.; Dsrh= 0.; Dgrg=0.; Dgrh=0.
2297
2298	!-------------------------------------------------------------------------------------------!
2299
2300	! COLLECTION by snow, graupel, hail:
2301	! (i.e. wet or dry ice-categories [=> excludes ice crystals])
2302
2303	! Collection by SNOW:
2304	if (QNM(i,k)>epsQ) then
2305	! cloud:
2306	if (QCM(i,k)>epsQ) then
2307
2308	!Approximation of Ecs based on Pruppacher & Klett (1997) Fig. 14-11
2309	Ecs= min(Dc,30.e-6)3.333e+4sqrt(min(Ds,1.e-3)*1.e+3)
2310	QCLcs= dtgamafscmrEcsPIov4iDE(i,k)(NCM(i,k)NNM(i,k))iGC5iGS31* &
2311	(GC13GS13iLAMc3iLAMsB2+2.GC14GS12iLAMc4iLAMsB1+GC15GS11* &
2312	iLAMc5*iLAMsB0)
2313
2314	NCLcs= dtgamafsPIov4Ecs(NCM(i,k)NNM(i,k))iGC5iGS31(GC5GS13* &
2315	iLAMsB2+2.GC11GS12iLAMciLAMsB1+GC12GS11iLAMc2*iLAMsB0)
2316
2317	!continuous collection: (alternative; gives values ~0.95 of SCE [above])
2318	!QCLcs= dtgamEcsPIov4afsQCM(i,k)NNM(i,k)iLAMs(2.+bfs)GS13*iGS31
2319	!NCLcs= QCLcs*NCM(i,k)/QCM(i,k)
2320
2321	!Correction factor for non-spherical snow [D = maximum dimension] which
2322	!changes projected area: [assumption: A=0.50D2 (vs. A=(PI/4)D**2)]
2323	! note: Strictly speaking, this correction should only be applied to
2324	! continuous growth approximation for cloud. [factor = 0.50/(pi/4)]
2325	if (.not. snowSpherical) then
2326	tmp1 = 0.6366 !factor = 0.50/(pi/4)
2327	QCLcs= tmp1*QCLcs
2328	NCLcs= tmp1*NCLcs
2329	endif
2330
2331	QCLcs= min(QCLcs, QCM(i,k))
2332	NCLcs= min(NCLcs, NCM(i,k))
2333	else
2334	QCLcs= 0.; NCLcs= 0.
2335	endif
2336
2337	! ice:
2338	if (QIM(i,k)>epsQ) then
2339	tmp1= vs0-vi0
2340	tmp3= sqrt(tmp1tmp1+0.04vs0*vi0)
2341
2342	QCLis= dtcmiiDE(i,k)PI6.Eis(NYM(i,k)NNM(i,k))tmp3iGI31iGS31(0.5 &
2343	iLAMs2iLAMi3+2.iLAMsiLAMi4+5.iLAMi5)
2344
2345	NCLis= dtPIov4Eis(NYM(i,k)NNM(i,k))GI31GS31tmp3(GI33GS31iLAMi2+ &
2346	2.GI32GS32iLAMiiLAMs+GI31GS33iLAMs2)
2347
2348	QCLis= min(QCLis, (QIM(i,k)))
2349	NCLis= min(QCLis(NYM(i,k)iQIM), NCLis)
2350	else
2351	QCLis= 0.; NCLis= 0.
2352	endif
2353
2354	if (dblMom_s) then
2355	!snow: (i.e. self-collection [aggregation])
2356	NCLss= dt0.93952Ess(DE(i,k)(QNM(i,k)))*((2.+bfs)thrd)(NNM(i,k))* &
2357	((4.-bfs)*thrd)
2358	!Note: 0.91226 = I(bfs)afsPI^((1-bfs)/3)*des^((-2-bfs)/3); I(bfs=0.41)=1138
2359	! 0.93952 = I(bfs)afsPI^((1-bfs)/3)*des^((-2-bfs)/3); I(bfs=0.42)=1172
2360	! [interpolated from 3rd-order polynomial approx. of values given in RRB98;
2361	! see eqn(A.35)]
2362	NCLss= min(NCLss, 0.5*(NNM(i,k)))
2363	endif
2364
2365	else
2366	QCLcs= 0.; NCLcs= 0.; QCLis= 0.; NCLis= 0.; NCLss= 0.
2367	endif
2368
2369	! Collection by GRAUPEL:
2370	if (QGM(i,k)>epsQ) then
2371
2372	! cloud:
2373	if (QCM(i,k)>epsQ) then
2374
2375	!(parameterization of Ecg based on Cober and List, 1993 [JAS])
2376	Kstoke = dewvg0DcDc/(9.MUdyn*Dg)
2377	Kstoke = max(1.5,min(10.,Kstoke))
2378	Ecg = 0.55log10(2.51Kstoke)
2379
2380	QCLcg= dtgamafgcmrEcgPIov4iDE(i,k)(NCM(i,k)NGM(i,k))iGC5iGG31* &
2381	(GC13GG13iLAMc3iLAMgB2+ 2.GC14GG12iLAMc4iLAMgB1+GC15GG11* &
2382	iLAMc5*iLAMgB0)
2383
2384	NCLcg= dtgamafgPIov4Ecg(NCM(i,k)NGM(i,k))iGC5iGG31(GC5GG13* &
2385	iLAMgB2+2.GC11GG12iLAMciLAMgB1+GC12GG11iLAMc2*iLAMgB0)
2386
2387	QCLcg= min(QCLcg, (QCM(i,k)))
2388	NCLcg= min(NCLcg, (NCM(i,k)))
2389	else
2390	QCLcg= 0.; NCLcg= 0.
2391	endif
2392
2393	! ice:
2394	if (QIM(i,k)>epsQ) then
2395	tmp1= vg0-vi0
2396	tmp3= sqrt(tmp1tmp1+0.04vg0*vi0)
2397
2398	QCLig= dtcmiiDE(i,k)PI6.Eig(NYM(i,k)NGM(i,k))tmp3iGI31iGG31(0.5 &
2399	iLAMg2iLAMi3+2.iLAMgiLAMi4+5.iLAMi5)
2400	NCLig= dtPIov4Eig(NYM(i,k)NGM(i,k))GI31GG31tmp3(GI33GG31iLAMi2+ &
2401	2.GI32GG32iLAMiiLAMg+GI31GG33iLAMg2)
2402
2403	QCLig= min(QCLig, (QIM(i,k)))
2404	NCLig= min(QCLig(NYM(i,k)iQIM), NCLig)
2405	else
2406	QCLig= 0.; NCLig= 0.
2407	endif
2408
2409	else
2410	QCLcg= 0.; QCLrg= 0.; QCLig= 0.
2411	NCLcg= 0.; NCLrg= 0.; NCLig= 0.
2412	endif
2413
2414	! Collection by HAIL:
2415	if (QHM(i,k)>epsQ) then
2416
2417	! cloud:
2418	if (QCM(i,k)>epsQ) then
2419	Ech = exp(-8.68e-7Dc(-1.6)Dh) !Ziegler (1985) A24
2420
2421	QCLch= dtgamafhcmrEchPIov4iDE(i,k)(NCM(i,k)NHM(i,k))iGC5iGH31* &
2422	(GC13GH13iLAMc3iLAMhB2+2.GC14GH12iLAMc4iLAMhB1+GC15GH11* &
2423	iLAMc5*iLAMhB0)
2424
2425	NCLch= dtgamafhPIov4Ech(NCM(i,k)NHM(i,k))iGC5iGH31(GC5GH13* &
2426	iLAMhB2+2.GC11GH12iLAMciLAMhB1+GC12GH11iLAMc2*iLAMhB0)
2427
2428	QCLch= min(QCLch, QCM(i,k))
2429	NCLch= min(NCLch, NCM(i,k))
2430	else
2431	QCLch= 0.; NCLch= 0.
2432	endif
2433
2434	! rain:
2435	if (QRM(i,k)>epsQ) then
2436	tmp1= vh0-vr0
2437	tmp3= sqrt(tmp1tmp1+0.04vh0*vr0)
2438	QCLrh= dtcmrErhPIov4iDE(i,k)(NHM(i,k)NRM(i,k))iGR31iGH31tmp3 &
2439	(GR36GH31iLAMr5+2.GR35GH32iLAMr4iLAMh+GR34GH33iLAMr3*iLAMh2)
2440
2441	NCLrh= dtPIov4Erh(NHM(i,k)NRM(i,k))iGR31iGH31tmp3(GR33GH31 &
2442	iLAMr2+2.GR32GH32iLAMriLAMh+GR31GH33iLAMh2)
2443
2444	QCLrh= min(QCLrh, QRM(i,k))
2445	NCLrh= min(NCLrh, QCLrh(NRM(i,k)iQRM))
2446	else
2447	QCLrh= 0.; NCLrh= 0.
2448	endif
2449
2450	! ice:
2451	if (QIM(i,k)>epsQ) then
2452	tmp1 = vh0-vi0
2453	tmp3 = sqrt(tmp1tmp1+0.04vh0*vi0)
2454
2455	QCLih= dtcmiiDE(i,k)PI6.Eih(NYM(i,k)NHM(i,k))tmp3iGI31iGH31(0.5 &
2456	iLAMh2iLAMi3+2.iLAMhiLAMi4+5.iLAMi5)
2457
2458	NCLih= dtPIov4Eih(NYM(i,k)NHM(i,k))GI31GH31tmp3(GI33GH31iLAMi2+ &
2459	2.GI32GH32iLAMiiLAMh+GI31GH33iLAMh2)
2460
2461	QCLih= min(QCLih, QIM(i,k))
2462	NCLih= min(QCLih(NYM(i,k)iQIM), NCLih)
2463	else
2464	QCLih= 0.; NCLih= 0.
2465	endif
2466
2467	! snow:
2468	if (QNM(i,k)>epsQ) then
2469	tmp1 = vh0-vs0
2470	tmp3 = sqrt(tmp1tmp1+0.04vh0*vs0)
2471	tmp4 = iLAMs2*iLAMs2
2472
2473	if (snowSpherical) then
2474	!hardcoded for dms=3:
2475	QCLsh= dtcmsiDE(i,k)PI6.Esh(NNM(i,k)NHM(i,k))tmp3iGS31iGH31* &
2476	(0.5iLAMh2iLAMs2iLAMs+2.iLAMhtmp4+5.tmp4*iLAMs)
2477	else
2478	!hardcoded for dms=2:
2479	QCLsh= dtcmsiDE(i,k)PI0.25Eshtmp3NNM(i,k)NHM(i,k)iGS31iGH31* &
2480	(GH33GS33iLAMh*2.iLAMs*2. + 2.GH32GS34iLAMhiLAMs*3. + &
2481	GH31GS35iLAMs**4.)
2482	endif
2483
2484	NCLsh= dtPIov4Esh(NNM(i,k)NHM(i,k))GS31GH31tmp3(GS33GH31iLAMs2+ &
2485	2.GS32GH32iLAMsiLAMh+GS31GH33iLAMh2)
2486
2487	QCLsh= min(QCLsh, (QNM(i,k)))
2488	NCLsh= min((NNM(i,k)iQNM)QCLsh, NCLsh, (NNM(i,k)))
2489	else
2490	QCLsh= 0.; NCLsh= 0.
2491	endif
2492
2493	!wet growth:
2494	VENTh= AvxGH32iLAMh*(2.+alpha_h) + BvxScTHRDsqrt(gamafhiMUkin)GH09* &
2495	iLAMh*(2.5+0.5bfh+alpha_h)
2496	QHwet= max(0., dtPI2(DE(i,k)CHLCCdiffDELqvs-KaTc)No_hiDE(i,k)/(CHLF+ &
2497	CPWTc)VENTh+(QCLihiEih+QCLshiEsh)(1.-CPITc/(CHLF+CPW*Tc)) )
2498
2499	else
2500	QCLch= 0.; QCLrh= 0.; QCLih= 0.; QCLsh= 0.; QHwet= 0.
2501	NCLch= 0.; NCLrh= 0.; NCLsh= 0.; NCLih= 0.
2502	endif
2503
2504	IF (TM(i,k)>TRPL .and. warmphase_ON) THEN
2505	!**********!
2506	! T > To !
2507	!**********!
2508
2509	! MELTING of frozen particles:
2510	! ICE:
2511	QMLir = QIM(i,k) !all pristine ice melts in one time step
2512	QIM(i,k)= 0.
2513	NMLir = NYM(i,k)
2514
2515	! SNOW:
2516	if (QNM(i,k)>epsQ) then
2517	QMLsr= dt(PI2iDE(i,k)iCHLFNo_sVENTs(KaTc-CHLCCdiffDELqvs) + CPW &
2518	iCHLFTc(QCLcs+QCLrs)*idt)
2519	QMLsr= min(max(QMLsr,0.), QNM(i,k))
2520	NMLsr= NNM(i,k)iQNMQMLsr
2521	else
2522	QMLsr= 0.; NMLsr= 0.
2523	endif
2524
2525	! GRAUPEL:
2526	if (QGM(i,k)>epsQ) then
2527	VENTg= AvxGG32iLAMgiLAMg+BvxScTHRDsqrt(gamafgiMUkin)GG09iLAMg* &
2528	(2.5+0.5*bfg+alpha_g)
2529	QMLgr= dt(PI2iDE(i,k)iCHLFNo_gVENTg(KaTc-CHLCCdiffDELqvs) + CPW &
2530	iCHLFTc(QCLcg+QCLrg)*idt)
2531	QMLgr= min(max(QMLgr,0.), QGM(i,k))
2532	NMLgr= NGM(i,k)iQGMQMLgr
2533	else
2534	QMLgr= 0.; NMLgr= 0.
2535	endif
2536
2537	! HAIL:
2538	if (QHM(i,k)>epsQ.and.Tc>5.) then
2539	VENTh= AvxGH32iLAMh*(2.+alpha_h) + BvxScTHRDsqrt(gamafhiMUkin)GH09* &
2540	iLAMh*(2.5+0.5bfh+alpha_h)
2541	QMLhr= dt(PI2iDE(i,k)iCHLFNo_hVENTh(KaTc-CHLCCdiff*DELqvs) + CPW/ &
2542	CHLFTc(QCLch+QCLrh)*idt)
2543	QMLhr= min(max(QMLhr,0.), QHM(i,k))
2544	NMLhr= NHM(i,k)iQHMQMLhr
2545	if(QCLrh>0.) NMLhr= NMLhr*0.1 !Prevents problems when hail is ML & CL
2546	else
2547	QMLhr= 0.; NMLhr= 0.
2548	endif
2549
2550	! Cold (sub-zero) source/sink terms:
2551	QNUvi= 0.; QFZci= 0.; QVDvi= 0.; QVDvs= 0.; QVDvg= 0.
2552	QCLis= 0.; QCNis1=0.; QCNis2=0.
2553	QCNgh= 0.; QIMsi= 0.; QIMgi= 0.; QCLir= 0.; QCLri= 0.
2554	QCLrs= 0.; QCLgr= 0.; QCLrg= 0.; QCNis= 0.; QVDvh= 0.
2555	QCNsg= 0.; QCLsr= 0.
2556
2557	NNUvi= 0.; NFZci= 0.; NCLgr= 0.; NCLrg= 0.; NgCNgh= 0.
2558	NCLis= 0.; NVDvi= 0.; NVDvs= 0.; NVDvg= 0.; NVDvh= 0.
2559	NCNsg= 0.; NhCNgh= 0.; NiCNis=0.; NsCNis=0.; NCLrs= 0.
2560	NIMsi= 0.; NIMgi= 0.; NCLir= 0.; NCLri= 0.; NCLsr= 0.
2561
2562	ELSE
2563	!----------!
2564	! T < To !
2565	!----------!
2566	tmp1 = 1./QSI(i,k)
2567	Si = QM(i,k) *tmp1
2568	tmp2 = TM(i,k)*TM(i,k)
2569	iABi = 1./( CHLSCHLS/(KaRGASVtmp2) + 1./(DE(i,k)(QSI(i,k))*Cdiff) )
2570
2571	! Warm-air-only source/sink terms:
2572	QMLir= 0.; QMLsr= 0.; QMLgr= 0.; QMLhr= 0.
2573	NMLir= 0.; NMLsr= 0.; NMLgr= 0.; NMLhr= 0.
2574
2575	!Probabilistic freezing (Bigg) of rain:
2576	if (Tc<Tc_FZrh .and. QRM(i,k)>epsQ .and. hail_ON) then
2577	!note: - (Tc<-10.C) condition is based on Pruppacher-Klett (1997) Fig. 9-41
2578	! - Small raindrops will freeze to hail. However, if after all S/S terms
2579	! are added Dh<Dh_min, then hail will be converted to graupel. Thus,
2580	! probabilistic freezing of small rain is effectively a source of graupel.
2581	NrFZrh= -dtBbigg(exp(AbiggTc)-1.)DE(i,k)QRM(i,k)idew
2582	Rz= 1. !N and Z (and Q) are conserved for FZrh with triple-moment
2583	! The Rz factor serves to conserve reflectivity when a rain distribution
2584	! converts to an distribution with a different shape parameter, alpha.
2585	! (e.g. when rain freezes to hail) The factor Rz non-conserves N while
2586	! acting to conserve Z for double-moment. See Ferrier, 1994 App. D)
2587	! Rz= (gamma(7.d0+alpha_h)GH31GR34GR34)/(GR36(i,k)GR31* &
2588	! gamma(4.d0+alpha_h)*gamma(4.d0+alpha_h))
2589	NhFZrh= Rz*NrFZrh
2590	QFZrh = NrFZrh(QRM(i,k)iNRM)
2591	else
2592	QFZrh= 0.; NrFZrh= 0.; NhFZrh= 0.
2593	endif
2594
2595	!--------!
2596	! ICE: !
2597	!--------!
2598	! Homogeneous freezing of cloud to ice:
2599	if (dblMom_c) then
2600	if (QCM(i,k)>epsQ) then
2601	tmp2 = TcTc; tmp3= tmp2Tc; tmp4= tmp2*tmp2
2602	JJ = (10.*max(-20.,(-606.3952-52.6611Tc-1.7439tmp2-0.0265tmp3- &
2603	1.536e-4*tmp4)))
2604	tmp1 = 1.e6(DE(i,k)(QCM(i,k)iNCM)icmr) !i.e. Dc[cm]**3
2605	FRAC = 1.-exp(-JJPIov6tmp1*dt)
2606	if (Tc>-30.) FRAC= 0.
2607	if (Tc<-50.) FRAC= 1.
2608	QFZci= FRAC*QCM(i,k)
2609	NFZci= FRAC*NCM(i,k)
2610	else
2611	QFZci= 0.; NFZci= 0.
2612	endif
2613	else
2614	!Homogeneous freezing of cloud to ice: (simplified)
2615	if (QCM(i,k)>epsQ .and. Tc<-35.) then
2616	FRAC= 1. !if T<-35
2617	QFZci= FRAC*QCM(i,k)
2618	NFZci= FRAC*N_c_SM
2619	else
2620	QFZci= 0.; NFZci= 0.
2621	endif
2622	endif
2623
2624	if (dblMom_i) then
2625	!Primary ice nucleation:
2626	NNUvi= 0.; QNUvi= 0.
2627	if (primIceNucl==1) then
2628
2629	NuDEPSOR= 0.; NuCONT= 0.
2630	Simax = min(Si, SxFNC(WZ(i,k),Tc,HPS(i)*S(i,k),QSW(i,k),QSI(i,k),CCNtype, &
2631	2))
2632	tmp1 = T(i,k)-7.66
2633	NNUmax = max(0., DE(i,k)/mio(Q(i,k)-QSS(i,k))/(1.+ck6(QSS(i,k)/(tmp1* &
2634	tmp1))))
2635	!Deposition/sorption nucleation:
2636	if (Tc<-5. .and. Si>1.) then
2637	NuDEPSOR= max(0., 1.e3exp(12.96(Simax-1.)-0.639)-(NYM(i,k))) !Meyers(1992)
2638	endif
2639	!Contact nucleation:
2640	if (QCM(i,k)>epsQ .and. Tc<-2.) then
2641	GG = 1.idew/(RGASV(TM(i,k))/((QSW(i,k)HPS(i)S(i,k))/EPS1)/ &
2642	Cdiff+CHLC/Ka/(TM(i,k))*(CHLC/RGASV/(TM(i,k))-1.)) !CP00a
2643	Swmax = SxFNC(WZ(i,k),Tc,HPS(i)*S(i,k),QSW(i,k),QSI(i,k),CCNtype,1)
2644	ssat = min((QM(i,k)/QSW(i,k)), Swmax) -1.
2645	Tcc = Tc + GGssatCHLC/Kdiff !C86_eqn64
2646	Na = exp(4.11-0.262*Tcc) !W95_eqn60/M92_2.6
2647	Kn = LAMa0(TM(i,k))p0/(T0(HPS(i)S(i,k))*Ra) !W95_eqn59
2648	PSIa = -kBoltzTcc/(6.piRaMUdyn)*(1.+Kn) !W95_eqn58
2649	ft = 0.4(1.+1.45Kn+0.4Knexp(-1./Kn))(Ka+2.5Kn*KAPa)/ &
2650	(1.+3.Kn)/(2.Ka+5.KAPaKn+KAPa) !W95_eqn57
2651	Dc = (DE(i,k)(QCM(i,k)iNCM)icmr)*thrd
2652	F1 = PI2DcNa*(NCM(i,k)) !W95_eqn55
2653	F2 = Ka/(HPS(i)S(i,k))(Tc-Tcc) !W95_eqn56
2654	NuCONTA= -F1F2RGASV(TM(i,k))/CHLCiDE(i,k) !diffusiophoresis
2655	NuCONTB= F1F2ft*iDE(i,k) !thermeophoresis
2656	NuCONTC= F1*PSIa !Brownian diffusion
2657	NuCONT = max(0.,(NuCONTA+NuCONTB+NuCONTC)*dt)
2658	endif
2659	!Total primary ice nucleation:
2660	if (icephase_ON) then
2661	NNUvi= min(NNUmax, NuDEPSOR + NuCONT )
2662	QNUvi= mioiDE(i,k)NNUvi
2663	QNUvi= min(QNUvi,(Q(i,k)))
2664	endif
2665
2666	elseif (primIceNucl==2) then
2667	if (Tc<-5. .and. Si>1.08) then !following Thompson etal (2006)
2668	NNUvi= max(N_Cooper(TRPL,T(i,k))-NYM(i,k),0.)
2669	QNUvi= min(mioiDE(i,k)NNUvi, Q(i,k))
2670	endif
2671	!elseif (primIceNucl==3) then
2672	!! (for alternative [future] ice nucleation parameterizations)
2673	! NNUvi=...
2674	! QNUvi=...
2675	endif !if (primIceNucl==1)
2676
2677	else !dblMom_i
2678	!Ice initiation (single-moment):
2679	if (QIM(i,k)<=epsQ .and. Tc<-5. .and. Si>1.08) then !following Thompson etal (2006)
2680	NNUvi = N_Cooper(TRPL,T(i,k))
2681	QNUvi= mioiDE(i,k)NNUvi
2682	QNUvi= min(QNUvi,Q(i,k))
2683	endif
2684	endif !dblMom_i
2685
2686
2687	IF (QIM(i,k)>epsQ) THEN
2688
2689	!Deposition/sublimation:
2690	! No_i = NYM(i,k)iGI31/iLAMi*(1.+alpha_i)
2691	! VENTi= AvxGI32iLAMi*(2.+alpha_i)+BvxScTHRDsqrt(gamafiiMUkin)GI6* &
2692	! iLAMi*(2.5+0.5bfi+alpha_i)
2693	No_i = NYM(i,k)*iGI31/iLAMi !optimized for alpha_i=0
2694	VENTi= AvxGI32iLAMiiLAMi+BvxScTHRDsqrt(gamafiiMUkin)GI6iLAMi* &
2695	(2.5+0.5*bfi+alpha_i)
2696	!Note: ice crystal capacitance is implicitly C = 0.5DcapFact_i
2697	QVDvi= dtcapFact_iiABi(PI2(Si-1.)No_iVENTi)
2698
2699	! Prevent overdepletion of vapor:
2700	tmp1 = T(i,k)-7.66
2701	VDmax = (Q(i,k)-QSS(i,k))/(1.+ck6(QSS(i,k))/(tmp1tmp1))
2702	if(Si>=1.) then
2703	QVDvi= min(max(QVDvi,0.),VDmax)
2704	else
2705	if (VDmax<0.) QVDvi= max(QVDvi,VDmax)
2706	!IF prevents subl.(QVDvi<0 at t) changing to dep.(VDmax>0 at t*) 2005-06-28
2707	endif
2708	if (.not. iceDep_ON) QVDvi= 0. !suppresses depositional growth
2709	NVDvi= min(0., (NYM(i,k)iQIM)QVDvi) !dNi/dt=0 for deposition
2710
2711	! Conversion to snow:
2712	! +depostion of ice:
2713	mi= DE(i,k)(QIM(i,k)iNYM)
2714	if (mi<=0.5mso.and.abs(0.5mso-mi)>1.e-20) then
2715	QCNis1= (mi/(mso-mi))*QVDvi
2716	else
2717	QCNis1= QVDvi + (1.-0.5mso/mi)QIM(i,k)
2718	endif
2719	QCNis1= max(0., QCNis1)
2720	! +aggregation of ice:
2721	if(Di<0.5*Dso) then
2722	Ki = PIov6DiDivi0Eii*Xdisp
2723	tmp1 = log(Di/Dso)
2724	tmp2 = tmp1tmp1tmp1
2725	QCNis2= -dt0.5(QIM(i,k)NYM(i,k))Ki/tmp2
2726	else
2727	Ki= 0.; QCNis2= 0.
2728	endif
2729	! +total conversion rate:
2730	QCNis = QCNis1 + QCNis2
2731	NsCNis= DE(i,k)imsoQCNis !source for snow (Ns)
2732	NiCNis= (DE(i,k)imsoQCNis1 + 0.5KiNYM(i,k)*NYM(i,k)) !sink for ice (Ni)
2733	NiCNis= min(NiCNis, NYM(i,k)*0.1) !Prevents overdepl. of NY when final QI>0
2734
2735	if (.not.(snow_ON)) then
2736	QCNis= 0.; NiCNis= 0.; NsCNis= 0. !Suppress SNOW initiation
2737	endif
2738
2739	! 3-component freezing (collisions with rain):
2740	if (QRM(i,k)>epsQ .and. QIM(i,k)>epsQ) then
2741	tmp1 = vr0-vi0
2742	tmp3 = sqrt(tmp1tmp1+0.04vr0*vi0)
2743
2744	QCLir= dtcmiEriPIov4iDE(i,k)(NRM(i,k)NYM(i,k))iGI31iGR31tmp3 &
2745	(GI36GR31iLAMi5+2.GI35GR32iLAMi4iLAMr+GI34GR33iLAMi3* &
2746	iLAMr2)
2747
2748	NCLri= dtPIov4Eri(NRM(i,k)NYM(i,k))iGI31iGR31tmp3(GI33GR31 &
2749	iLAMi2+2.GI32GR32iLAMiiLAMr+GI31GR33iLAMr2)
2750
2751	QCLri= dtcmrEriPIov4iDE(i,k)(NYM(i,k)NRM(i,k))iGR31iGI31tmp3 &
2752	(GR36GI31 iLAMr5+2.GR35GI32iLAMr4iLAMi+GR34GI33iLAMr3* &
2753	iLAMi2)
2754
2755	!note: For explicit eqns, both NCLri and NCLir are mathematically identical)
2756	NCLir= min(QCLir(NYM(i,k)iQIM), NCLri)
2757	QCLri= min(QCLri, (QRM(i,k))); QCLir= min(QCLir, (QIM(i,k)))
2758	NCLri= min(NCLri, (NRM(i,k))); NCLir= min(NCLir, (NYM(i,k)))
2759
2760	!Determine destination of 3-comp.freezing:
2761	tmp1= max(Di,Dr)
2762	dey= (deiDiDiDi+dewDrDrDr)/(tmp1tmp1tmp1)
2763	if (dey>0.5*(deg+deh) .and. Dr>Dr_3cmpThrs .and. hail_ON) then
2764	Dirg= 0.; Dirh= 1.
2765	else
2766	Dirg= 1.; Dirh= 0.
2767	endif
2768	if (.not. grpl_ON) Dirg= 0.
2769
2770	else
2771	QCLir= 0.; NCLir= 0.; QCLri= 0.
2772	NCLri= 0.; Dirh = 0.; Dirg= 0.
2773	endif
2774
2775	! Rime-splintering (ice multiplication):
2776	ff= 0.
2777	if(Tc>=-8..and.Tc<=-5.) ff= 3.5e8(Tc +8.)thrd
2778	if(Tc> -5..and.Tc< -3.) ff= 3.5e8(-3.-Tc)0.5
2779	NIMsi= DE(i,k)ffQCLcs
2780	NIMgi= DE(i,k)ffQCLcg
2781	QIMsi= mioiDE(i,k)NIMsi
2782	QIMgi= mioiDE(i,k)NIMgi
2783
2784	ELSE
2785
2786	QVDvi= 0.; QCNis= 0.
2787	QIMsi= 0.; QIMgi= 0.; QCLri= 0.; QCLir= 0.
2788	NVDvi= 0.; NCLir= 0.; NIMsi= 0.
2789	NiCNis=0.; NsCNis=0.; NIMgi= 0.; NCLri= 0.
2790
2791	ENDIF
2792	!---------!
2793	! SNOW: !
2794	!---------!
2795	IF (QNM(i,k)>epsQ) THEN
2796
2797	!Deposition/sublimation:
2798	!note: - snow crystal capacitance is implicitly C = 0.5DcapFact_s
2799	! - No_s and VENTs are computed above
2800	QVDvs = dtcapFact_siABi(PI2(Si-1.)No_sVENTs - CHLSCHLF/(KaRGASV* &
2801	TM(i,k)TM(i,k))QCLcs*idt)
2802
2803	! Prevent overdepletion of vapor:
2804	tmp1 = T(i,k)-7.66
2805	VDmax = (Q(i,k)-QSS(i,k))/(1.+ck6(QSS(i,k))/(tmp1tmp1)) !KY97_A.33
2806	if(Si>=1.) then
2807	QVDvs= min(max(QVDvs,0.),VDmax)
2808	else
2809	if (VDmax<0.) QVDvs= max(QVDvs,VDmax)
2810	!IF prevents subl.(QVDvs<0 at t) changing to dep.(VDmax>0 at t*)
2811	endif
2812	NVDvs= -min(0.,(NNM(i,k)iQNM)QVDvs) !pos. quantity
2813
2814	! Conversion to graupel:
2815	if (QCLcs>CNsgThresQVDvs .and. 0.99deg>des) then
2816	!note: The (deg>des) condition equates to (Ds>330microns) for m(D)=0.069D^2
2817	! relation for snow, which implies a variable bulk density. The physical
2818	! assumption in the QCNsg equation is that snow converts to graupel due
2819	! to densification from riming.
2820	! The 0.99 is to prevent overflow if des~deg
2821	QCNsg= (deg/(deg-des))*QCLcs
2822	else
2823	QCNsg= 0.
2824	endif
2825	if (.not. grpl_ON) QCNsg= 0.
2826	NCNsg= DE(i,k)imgoQCNsg
2827	NCNsg= min(NCNsg, (0.5NNM(i,k)iQNM)*QCNsg) !Prevents incorrect Ns-depletion
2828
2829	! 3-component freezing (collisions with rain):
2830	if (QRM(i,k)>epsQ .and. QNM(i,k)>epsQ .and. Tc<-5.) then
2831	tmp1 = vs0-vr0
2832	tmp2 = sqrt(tmp1tmp1+0.04vs0*vr0)
2833	tmp6 = iLAMs2iLAMs2iLAMs
2834
2835	QCLrs= dtcmrErsPIov4iDE(i,k)NNM(i,k)NRM(i,k)iGR31iGS31tmp2 &
2836	(GR36GS31iLAMr5+2.GR35GS32iLAMr4iLAMs+GR34GS33iLAMr3* &
2837	iLAMs2)
2838
2839	NCLrs= dt0.25e0PIErs(NNM(i,k)NRM(i,k))iGR31iGS31tmp2(GR33 &
2840	GS31iLAMr2+2.GR32GS32iLAMriLAMs+GR31GS33*iLAMs2)
2841
2842	if (snowSpherical) then
2843	!hardcoded for dms=3:
2844	QCLsr= dtcmsErsPIov4iDE(i,k)(NRM(i,k)NNM(i,k))iGS31iGR31* &
2845	tmp2(GS36GR31tmp6+2.GS35GR32iLAMs2iLAMs2iLAMr+GS34* &
2846	GR33iLAMs2iLAMs*iLAMr2)
2847	else
2848	!hardcoded for dms=2:
2849	QCLsr= dtcmsiDE(i,k)PI0.25ERStmp2NNM(i,k)NRM(i,k)iGS31 &
2850	iGR31(GR33GS33iLAMr2.iLAMs*2. + 2.GR32GS34iLAMr* &
2851	iLAMs*3. +GR31GS35iLAMs*4.)
2852	endif
2853
2854	!note: For explicit eqns, NCLsr = NCLrs
2855	NCLsr= min(QCLsr(NNM(i,k)iQNM), NCLrs)
2856	QCLrs= min(QCLrs, QRM(i,k)); QCLsr= min(QCLsr, QNM(i,k))
2857	NCLrs= min(NCLrs, NRM(i,k)); NCLsr= min(NCLsr, NNM(i,k))
2858
2859	! Determine destination of 3-comp.freezing:
2860	Dsrs= 0.; Dsrg= 0.; Dsrh= 0.
2861	tmp1= max(Ds,Dr)
2862	tmp2= tmp1tmp1tmp1
2863	dey = (desDsDsDs + dewDrDrDr)/tmp2
2864	if (dey<=0.5*(des+deg) ) Dsrs= 1. !snow
2865	if (dey >0.5(des+deg) .and. dey<0.5(deg+deh)) Dsrg= 1. !graupel
2866	if (dey>=0.5*(deg+deh)) then
2867	Dsrh= 1. !hail
2868	if (.not.hail_ON .or. Dr<Dr_3cmpThrs) then
2869	Dsrg= 1.; Dsrh= 0. !graupel
2870	endif
2871	endif
2872	if (.not. grpl_ON) Dsrg=0.
2873	else
2874	QCLrs= 0.; QCLsr= 0.; NCLrs= 0.; NCLsr= 0.
2875	endif
2876
2877	ELSE
2878
2879	QVDvs= 0.; QCLcs= 0.; QCNsg= 0.; QCLsr= 0.; QCLrs= 0.
2880	NVDvs= 0.; NCLcs= 0.; NCLsr= 0.; NCLrs= 0.; NCNsg= 0.
2881
2882	ENDIF
2883	!------------!
2884	! GRAUPEL: !
2885	!------------!
2886	IF (QGM(i,k)>epsQ) THEN
2887
2888	!Conversion to hail: (D_sll given by S-L limit)
2889	if (WZ(i,k)>w_CNgh .and. hail_ON) then
2890	D_sll = 0.01(exp(min(20.,-Tc/(1.1e4DE(i,k)(QCM(i,k)+QRM(i,k))-1.3e3 &
2891	DE(i,k)*(QIM(i,k))+1.)))-1.)
2892	!Add correction factor: [to account error in equation of Ziegler (1985), as per Young (1993)]
2893	D_sll = 2.0*D_sll
2894	D_sll = min(1., max(0.0001,D_sll)) !smallest D_sll=0.1mm; largest=1m
2895
2896	!Old approach: (pre-my-2.15.0)
2897	! ratio= Dg/D_sll
2898	! if (ratio>r_CNgh) then
2899	! QCNgh= (0.5ratio)(QCLcg+QCLrg+QCLig)
2900	! QCNgh= min(QCNgh,(QGM(i,k))+QCLcg+QCLrg+QCLig)
2901	! NCNgh= DE(i,k)QCNghicmh/(D_sllD_sllD_sll)
2902	! else
2903	! QCNgh= 0.
2904	! NCNgh= 0.
2905	! endif
2906	!New approach:
2907	tmp1 = exp(-D_sll/iLAMg)
2908	Ng_tail = No_giLAMgtmp1 !integral(Dsll,inf) of N(D)dD
2909	if (Ng_tail > Ngh_crit) then
2910	QCNgh = idtcmgNo_gtmp1(D_sll*3.iLAMg + 3.D_sll2.iLAMg**2. &
2911	+ 6.D_slliLAMg*3. + 6.iLAMg**4.)
2912	NgCNgh= idtNo_giLAMg*tmp1
2913	Rz= 1.
2914	!---
2915	! The Rz factor (<>1) serves to conserve reflectivity when graupel
2916	! converts to hail with a a different shape parameter, alpha.
2917	! The factor Rz non-conserves N while acting to conserve Z for
2918	! double-moment. See Ferrier, 1994 App. D). However, Rz=1 is
2919	! used since it is deemed more important to conserve concentration
2920	! than reflectivity (see Milbrandt and McTaggart-Cowan, 2010 JAS).
2921	!---
2922	! Code to conserve total reflectivity:
2923	! if (QHM(i,k)>epsQ) then
2924	! Rz= (gamma(7.+alpha_h)GH31GG34*2.)/(GG36GG31GH34*2.)
2925	! else
2926	! Rz= 1.
2927	! endif
2928	!---
2929	NhCNgh= Rz*NgCNgh
2930	else
2931	QCNgh = 0.; NgCNgh = 0.; NhCNgh = 0.
2932	endif
2933	endif
2934
2935	!3-component freezing (collisions with rain):
2936	if (QRM(i,k)>epsQ) then
2937	tmp1 = vg0-vr0
2938	tmp2 = sqrt(tmp1tmp1 + 0.04vg0*vr0)
2939	tmp8 = iLAMg2iLAMg ! iLAMg*3
2940	tmp9 = tmp8iLAMg ! iLAMg*4
2941	tmp10= tmp9iLAMg ! iLAMg*5
2942
2943	!(parameterization of Erg based on Cober and List, 1993 [JAS])
2944	Kstoke = dewabs(vg0-vr0)DrDr/(9.MUdyn*Dg)
2945	Kstoke = max(1.5,min(10.,Kstoke))
2946	Erg = 0.55log10(2.51Kstoke)
2947
2948	QCLrg= dtcmrErgPIov4iDE(i,k)(NGM(i,k)NRM(i,k))iGR31iGG31tmp2 &
2949	(GR36GG31iLAMr5+2.GR35GG32iLAMr4iLAMg+GR34GG33iLAMr3* &
2950	iLAMg2)
2951
2952	NCLrg= dtPIov4Erg(NGM(i,k)NRM(i,k))iGR31iGG31tmp2(GR33GG31 &
2953	iLAMr2+2.GR32GG32iLAMriLAMg+GR31GG33iLAMg2)
2954
2955	QCLgr= dtcmgErgPIov4iDE(i,k)(NRM(i,k)NGM(i,k))iGG31iGR31tmp2 &
2956	(GG36GR31tmp10+2.GG35GR32tmp9iLAMr+GG34GR33tmp8*iLAMr2)
2957
2958	!(note: For explicit eqns, NCLgr= NCLrg)
2959	NCLgr= min(NCLrg, QCLgr(NGM(i,k)iQGM))
2960	QCLrg= min(QCLrg, QRM(i,k)); QCLgr= min(QCLgr, QGM(i,k))
2961	NCLrg= min(NCLrg, NRM(i,k)); NCLgr= min(NCLgr, NGM(i,k))
2962
2963	! Determine destination of 3-comp.freezing:
2964	tmp1= max(Dg,Dr)
2965	tmp2= tmp1tmp1tmp1
2966	dey = (degDgDgDg + dewDrDrDr)/tmp2
2967	if (dey>0.5*(deg+deh) .and. Dr>Dr_3cmpThrs .and. hail_ON) then
2968	Dgrg= 0.; Dgrh= 1.
2969	else
2970	Dgrg= 1.; Dgrh= 0.
2971	endif
2972	else
2973	QCLgr= 0.; QCLrg= 0.; NCLgr= 0.; NCLrg= 0.
2974	endif
2975
2976	ELSE
2977
2978	QVDvg= 0.; QCNgh= 0.; QCLgr= 0.; QCLrg= 0.; NgCNgh= 0.
2979	NVDvg= 0.; NhCNgh= 0.; NCLgr= 0.; NCLrg= 0.
2980
2981	ENDIF
2982	!---------!
2983	! HAIL: !
2984	!---------!
2985	IF (QHM(i,k)>epsQ) THEN
2986
2987	!Wet growth:
2988	if (QHwet<(QCLch+QCLrh+QCLih+QCLsh) .and. Tc>-40.) then
2989	QCLih= min(QCLih*iEih, QIM(i,k)) !change Eih to 1. in CLih
2990	NCLih= min(NCLih*iEih, NYM(i,k)) ! " "
2991	QCLsh= min(QCLsh*iEsh, QNM(i,k)) !change Esh to 1. in CLsh
2992	NCLsh= min(NCLsh*iEsh, NNM(i,k)) ! " "
2993	tmp3 = QCLrh
2994	QCLrh= QHwet-(QCLch+QCLih+QCLsh) !actual QCLrh minus QSHhr
2995	QSHhr= tmp3-QCLrh !QSHhr used here only
2996	NSHhr= DE(i,k)QSHhr/(cmrDrshedDrshedDrshed)
2997	else
2998	NSHhr= 0.
2999	endif
3000	ELSE
3001	QVDvh= 0.; NVDvh= 0.; NSHhr= 0.
3002	ENDIF
3003
3004	ENDIF ! ( if Tc<0C Block )
3005
3006	!------------ End of source/sink term calculation -------------!
3007
3008	!-- Adjustment of source/sink terms to prevent overdepletion: --!
3009	do niter= 1,2
3010
3011	! (1) Vapor:
3012	source= Q(i,k) +dim(-QVDvi,0.)+dim(-QVDvs,0.)+dim(-QVDvg,0.)+dim(-QVDvh,0.)
3013	sink = QNUvi+dim(QVDvi,0.)+dim(QVDvs,0.)
3014	sour = max(source,0.)
3015	if(sink>sour) then
3016	ratio= sour/sink
3017	QNUvi= ratioQNUvi; NNUvi= ratioNNUvi
3018	if(QVDvi>0.) then
3019	QVDvi= ratioQVDvi; NVDvi= ratioNVDvi
3020	endif
3021	if(QVDvs>0.) then
3022	QVDvs=ratioQVDvs; NVDvs=ratioNVDvs
3023	endif
3024	QVDvg= ratioQVDvg; NVDvg= ratioNVDvg
3025	QVDvh= ratioQVDvh; NVDvh= ratioNVDvh
3026	endif
3027
3028	! (2) Cloud:
3029	source= QC(i,k)
3030	sink = QCLcs+QCLcg+QCLch+QFZci
3031	sour = max(source,0.)
3032	if(sink>sour) then
3033	ratio= sour/sink
3034	QFZci= ratioQFZci; NFZci= ratioNFZci
3035	QCLcs= ratioQCLcs; NCLcs= ratioNCLcs
3036	QCLcg= ratioQCLcg; NCLcg= ratioNCLcg
3037	QCLch= ratioQCLch; NCLch= ratioNCLch
3038	endif
3039
3040	! (3) Rain:
3041	source= QR(i,k)+QMLsr+QMLgr+QMLhr+QMLir
3042	sink = QCLri+QCLrs+QCLrg+QCLrh+QFZrh
3043	sour = max(source,0.)
3044	if(sink>sour) then
3045	ratio= sour/sink
3046	QCLrg= ratioQCLrg; QCLri= ratioQCLri; NCLri= ratio*NCLri
3047	QCLrs= ratioQCLrs; NCLrs= ratioNCLrs; QCLrg= ratio*QCLrg
3048	NCLrg= ratioNCLrg; QCLrh= ratioQCLrh; NCLrh= ratio*NCLrh
3049	QFZrh= ratioQFZrh; NrFZrh=ratioNrFZrh; NhFZrh=ratio*NhFZrh
3050	if (ratio==0.) then
3051	Dirg= 0.; Dirh= 0.; Dgrg= 0.; Dgrh= 0.
3052	Dsrs= 0.; Dsrg= 0.; Dsrh= 0.
3053	endif
3054	endif
3055
3056	! (4) Ice:
3057	source= QI(i,k)+QNUvi+dim(QVDvi,0.)+QFZci
3058	sink = QCNis+QCLir+dim(-QVDvi,0.)+QCLis+QCLig+QCLih+QMLir
3059	sour = max(source,0.)
3060	if(sink>sour) then
3061	ratio= sour/sink
3062	QMLir= ratioQMLir; NMLir= ratioNMLir
3063	if (QVDvi<0.) then
3064	QVDvi= ratioQVDvi; NVDvi= ratioNVDvi
3065	endif
3066	QCNis= ratioQCNis; NiCNis= ratioNiCNis; NsCNis= ratio*NsCNis
3067	QCLir= ratioQCLir; NCLir= ratioNCLir; QCLig= ratio*QCLig
3068	QCLis= ratioQCLis; NCLis= ratioNCLis
3069	QCLih= ratioQCLih; NCLih= ratioNCLih
3070	if (ratio==0.) then
3071	Dirg= 0.; Dirh= 0.
3072	endif
3073	endif
3074
3075	! (5) Snow:
3076	source= QN(i,k)+QCNis+dim(QVDvs,0.)+QCLis+Dsrs*(QCLrs+QCLsr)+QCLcs
3077	sink = dim(-QVDvs,0.)+QCNsg+QMLsr+QCLsr+QCLsh
3078	sour = max(source,0.)
3079	if(sink>sour) then
3080	ratio= sour/sink
3081	if(QVDvs<=0.) then
3082	QVDvs= ratioQVDvs; NVDvs= ratioNVDvs
3083	endif
3084	QCNsg= ratioQCNsg; NCNsg= ratioNCNsg; QMLsr= ratio*QMLsr
3085	NMLsr= ratioNMLsr; QCLsr= ratioQCLsr; NCLsr= ratio*NCLsr
3086	QCLsh= ratioQCLsh; NCLsh= ratioNCLsh
3087	if (ratio==0.) then
3088	Dsrs= 0.; Dsrg= 0.; Dsrh= 0.
3089	endif
3090	endif
3091
3092	! (6) Graupel:
3093	source= QG(i,k)+QCNsg+dim(QVDvg,0.)+Dirg(QCLri+QCLir)+Dgrg(QCLrg+QCLgr)+ &
3094	QCLcg+Dsrg*(QCLrs+QCLsr)+QCLig
3095	sink = dim(-QVDvg,0.)+QMLgr+QCNgh+QCLgr
3096	sour = max(source,0.)
3097	if(sink>sour) then
3098	ratio= sour/sink
3099	QVDvg= ratioQVDvg; NVDvg= ratioNVDvg; QMLgr = ratio*QMLgr
3100	NMLgr= ratioNMLgr; QCNgh= ratioQCNgh; NgCNgh= ratio*NgCNgh
3101	QCLgr= ratioQCLgr; NCLgr= ratioNCLgr; NhCNgh= ratio*NhCNgh
3102	if (ratio==0.) then
3103	Dgrg= 0.; Dgrh= 0.
3104	endif
3105	endif
3106
3107	! (7) Hail:
3108	source= QH(i,k)+dim(QVDvh,0.)+QCLch+QCLrh+Dirh*(QCLri+QCLir)+QCLih+QCLsh+ &
3109	Dsrh(QCLrs+QCLsr)+QCNgh+Dgrh(QCLrg+QCLgr)+QFZrh
3110	sink = dim(-QVDvh,0.)+QMLhr
3111	sour = max(source,0.)
3112	if(sink>sour) then
3113	ratio= sour/sink
3114	QVDvh= ratioQVDvh; NVDvh= ratioNVDvh
3115	QMLhr= ratioQMLhr; NMLhr= ratioNMLhr
3116	endif
3117
3118	enddo
3119	!--------------- End of source/sink term adjustment ------------------!
3120
3121	!Compute N-tendencies for destination categories of 3-comp.freezing:
3122	NCLirg= 0.; NCLirh= 0.; NCLsrs= 0.; NCLsrg= 0.
3123	NCLsrh= 0.; NCLgrg= 0.; NCLgrh= 0.
3124
3125	if (QCLir+QCLri>0.) then
3126	tmp1 = max(Dr,Di)
3127	tmp2 = tmp1tmp1tmp1*PIov6
3128	NCLirg= DirgDE(i,k)(QCLir+QCLri)/(deg*tmp2)
3129	NCLirh= DirhDE(i,k)(QCLir+QCLri)/(deh*tmp2)
3130	endif
3131
3132	if (QCLsr+QCLrs>0.) then
3133	tmp1 = max(Dr,Ds)
3134	tmp2 = tmp1tmp1tmp1*PIov6
3135	NCLsrs= DsrsDE(i,k)(QCLsr+QCLrs)/(des*tmp2)
3136	NCLsrg= DsrgDE(i,k)(QCLsr+QCLrs)/(deg*tmp2)
3137	NCLsrh= DsrhDE(i,k)(QCLsr+QCLrs)/(deh*tmp2)
3138	endif
3139
3140	if (QCLgr+QCLrg>0.) then
3141	tmp1 = max(Dr,Dg)
3142	tmp2 = tmp1tmp1tmp1*PIov6
3143	NCLgrg= DgrgDE(i,k)(QCLgr+QCLrg)/(deg*tmp2)
3144	NCLgrh= DgrhDE(i,k)(QCLgr+QCLrg)/(deh*tmp2)
3145	endif
3146
3147	!========================================================================!
3148	! Add all source/sink terms to all predicted moments: !
3149	!========================================================================!
3150
3151	!Diagnostic S/S terms: (to facilitate output of 3D variables for diagnostics)
3152	!SS01(i,k)= QVDvs*idt (e.g., for depositional growth rate of snow, kg kg-1 s-1)
3153
3154	! Q-Source/Sink Terms:
3155	Q(i,k) = Q(i,k) -QNUvi -QVDvi -QVDvs -QVDvg -QVDvh
3156	QC(i,k)= QC(i,k) -QCLcs -QCLcg -QCLch -QFZci
3157	QR(i,k)= QR(i,k) -QCLri +QMLsr -QCLrs -QCLrg +QMLgr -QCLrh +QMLhr -QFZrh +QMLir
3158	QI(i,k)= QI(i,k) +QNUvi +QVDvi +QFZci -QCNis -QCLir -QCLis -QCLig &
3159	-QMLir -QCLih +QIMsi +QIMgi
3160	QG(i,k)= QG(i,k) +QCNsg +QVDvg +QCLcg -QCLgr-QMLgr -QCNgh -QIMgi +QCLig &
3161	+Dirg(QCLri+QCLir) +Dgrg(QCLrg+QCLgr) +Dsrg*(QCLrs+QCLsr)
3162	QN(i,k)= QN(i,k) +QCNis +QVDvs +QCLcs -QCNsg -QMLsr -QIMsi -QCLsr +QCLis -QCLsh &
3163	+Dsrs*(QCLrs+QCLsr)
3164	QH(i,k)= QH(i,k) +Dirh(QCLri+QCLir) -QMLhr +QVDvh +QCLch +Dsrh(QCLrs+QCLsr) &
3165	+QCLih +QCLsh +QFZrh +QCLrh +QCNgh +Dgrh*(QCLrg+QCLgr)
3166
3167	! N-Source/Sink Terms:
3168	if (dblMom_c) NC(i,k)= NC(i,k) -NCLcs -NCLcg -NCLch -NFZci
3169	if (dblMom_r) NR(i,k)= NR(i,k) -NCLri -NCLrs -NCLrg -NCLrh +NMLsr +NMLgr +NMLhr &
3170	-NrFZrh +NMLir +NSHhr
3171	if (dblMom_i) NY(i,k)= NY(i,k) +NNUvi +NVDvi +NFZci -NCLir -NCLis -NCLig -NCLih &
3172	-NMLir +NIMsi +NIMgi -NiCNis
3173	if (dblMom_s) NN(i,k)= NN(i,k) +NsCNis -NVDvs -NCNsg -NMLsr -NCLss -NCLsr -NCLsh &
3174	+NCLsrs
3175	if (dblMom_g) NG(i,k)= NG(i,k) +NCNsg -NCLgr -NVDvg -NMLgr +NCLirg +NCLsrg &
3176	+NCLgrg -NgCNgh
3177	if (dblMom_h) NH(i,k)= NH(i,k) +NhFZrh +NhCNgh -NMLhr -NVDvh +NCLirh +NCLsrh &
3178	+NCLgrh
3179
3180	T(i,k)= T(i,k) +LFP*(QCLri+QCLcs+QCLrs+QFZci-QMLsr+QCLcg+QCLrg-QMLir-QMLgr &
3181	-QMLhr+QCLch+QCLrh+QFZrh) +LSP*(QNUvi+QVDvi+QVDvs+QVDvg+QVDvh)
3182
3183	!Prevent overdepletion:
3184	IF (dblMom_c) THEN
3185	if(QC(i,k)<epsQ .or. NC(i,k)<epsN) then
3186	Q(i,k) = Q(i,k) + QC(i,k)
3187	T(i,k) = T(i,k) - LCP*QC(i,k)
3188	QC(i,k)= 0.; NC(i,k)= 0.
3189	endif
3190	ELSE
3191	if(QC(i,k)<epsQ) then
3192	Q(i,k) = Q(i,k) + QC(i,k)
3193	T(i,k) = T(i,k) - LCP*QC(i,k)
3194	QC(i,k)= 0.
3195	endif
3196	ENDIF
3197
3198	IF (dblMom_r) THEN
3199	if (QR(i,k)<epsQ .or. NR(i,k)<epsN) then
3200	Q(i,k) = Q(i,k) + QR(i,k)
3201	T(i,k) = T(i,k) - LCP*QR(i,k)
3202	QR(i,k)= 0.; NR(i,k)= 0.
3203	endif
3204	ELSE
3205	if (QR(i,k)<epsQ) then
3206	Q(i,k) = Q(i,k) + QR(i,k)
3207	T(i,k) = T(i,k) - LCP*QR(i,k)
3208	QR(i,k)= 0.
3209	endif
3210	ENDIF
3211
3212	IF (dblMom_i) THEN
3213	if (QI(i,k)<epsQ .or. NY(i,k)<epsN) then
3214	Q(i,k) = Q(i,k) + QI(i,k)
3215	T(i,k) = T(i,k) - LSP*QI(i,k)
3216	QI(i,k)= 0.; NY(i,k)= 0.
3217	endif
3218	ELSE
3219	if (QI(i,k)<epsQ) then
3220	Q(i,k) = Q(i,k) + QI(i,k)
3221	T(i,k) = T(i,k) - LSP*QI(i,k)
3222	QI(i,k)= 0.
3223	endif
3224	ENDIF
3225
3226	IF (dblMom_s) THEN
3227	if (QN(i,k)<epsQ .or. NN(i,k)<epsN) then
3228	Q(i,k) = Q(i,k) + QN(i,k)
3229	T(i,k) = T(i,k) - LSP*QN(i,k)
3230	QN(i,k)= 0.; NN(i,k)= 0.
3231	endif
3232	ELSE
3233	if (QN(i,k)<epsQ) then
3234	Q(i,k) = Q(i,k) + QN(i,k)
3235	T(i,k) = T(i,k) - LSP*QN(i,k)
3236	QN(i,k)= 0.
3237	endif
3238	ENDIF
3239
3240	IF (dblMom_g) THEN
3241	if (QG(i,k)<epsQ .or. NG(i,k)<epsN) then
3242	Q(i,k) = Q(i,k) + QG(i,k)
3243	T(i,k) = T(i,k) - LSP*QG(i,k)
3244	QG(i,k)= 0.; NG(i,k)= 0.
3245	endif
3246	ELSE
3247	if (QG(i,k)<epsQ) then
3248	Q(i,k) = Q(i,k) + QG(i,k)
3249	T(i,k) = T(i,k) - LSP*QG(i,k)
3250	QG(i,k)= 0.
3251	endif
3252	ENDIF
3253
3254	IF (dblMom_h) THEN
3255	if (QH(i,k)<epsQ .or. NH(i,k)<epsN) then
3256	Q(i,k) = Q(i,k) + QH(i,k)
3257	T(i,k) = T(i,k) - LSP*QH(i,k)
3258	QH(i,k)= 0.; NH(i,k)= 0.
3259	else if (QH(i,k)>epsQ .and. NH(i,k)>epsN) then
3260	!Conversion to graupel of hail is small:
3261	Dh= (DE(i,k)QH(i,k)/NH(i,k)icmh)**thrd
3262	if (Dh<Dh_min) then
3263	QG(i,k)= QG(i,k) + QH(i,k)
3264	NG(i,k)= NG(i,k) + NH(i,k)
3265	QH(i,k)= 0.; NH(i,k)= 0.
3266	endif
3267	endif
3268	ELSE
3269	if (QH(i,k)<epsQ) then
3270	Q(i,k) = Q(i,k) + QH(i,k)
3271	T(i,k) = T(i,k) - LSP*QH(i,k)
3272	QH(i,k)= 0.
3273	endif
3274	ENDIF
3275	Q(i,k)= max(Q(i,k),0.)
3276	NY(i,k)= min(NY(i,k), Ni_max)
3277
3278	ENDIF !if (activePoint)
3279	ENDDO
3280	ENDDO
3281
3282	!----------------------------------------------------------------------------------!
3283	! End of ice phase microphysics (Part 2) !
3284	!----------------------------------------------------------------------------------!
3285
3286	!----------------------------------------------------------------------------------!
3287	! PART 3: Warm Microphysics Processes !
3288	! !
3289	! Equations for warm-rain coalescence based on Cohard and Pinty (2000a,b; QJRMS) !
3290	! Condensation/evaportaion equations based on Kong and Yau (1997; Atmos-Ocean) !
3291	! Equations for rain reflectivity (ZR) based on Milbrandt and Yau (2005b; JAS) !
3292	!----------------------------------------------------------------------------------!
3293
3294	! Part 3a - Warm-rain Coallescence:
3295
3296	IF (warmphase_ON) THEN
3297
3298	DO k= 2,nk
3299	DO i= 1,ni
3300
3301	RCAUTR= 0.; CCACCR= 0.; Dc= 0.; iLAMc= 0.; L = 0.
3302	RCACCR= 0.; CCSCOC= 0.; Dr= 0.; iLAMr= 0.; TAU= 0.
3303	CCAUTR= 0.; CRSCOR= 0.; SIGc= 0.; DrINIT= 0.
3304	iLAMc3= 0.; iLAMc6= 0.; iLAMr3= 0.; iLAMr6= 0.
3305
3306	if (dblMom_r) then
3307	rainPresent= (QRM(i,k)>epsQ .and. NRM(i,k)>epsN)
3308	else
3309	rainPresent= (QRM(i,k)>epsQ)
3310	endif
3311
3312	if (.not. dblMom_c) NCM(i,k)= N_c_SM
3313	if (QCM(i,k)>epsQ .and. NCM(i,k)>epsN) then
3314	iLAMc = iLAMDA_x(DE(i,k),QCM(i,k),1./NCM(i,k),icexc9,thrd)
3315	iLAMc3= iLAMciLAMciLAMc
3316	iLAMc6= iLAMc3*iLAMc3
3317	Dc = iLAMc(GC2iGC1)**thrd
3318	SIGc = iLAMc( GC3iGC1- (GC2iGC1)(GC2iGC1) )*sixth
3319	L = 0.027DE(i,k)QCM(i,k)(6.25e18SIGcSIGcSIGc*Dc-0.4)
3320	if (SIGc>SIGcTHRS) TAU= 3.7/(DE(i,k)(QCM(i,k))(0.5e6*SIGc-7.5))
3321	endif
3322
3323	if (rainPresent) then
3324	if (dblMom_r) then
3325	Dr = Dm_x(DE(i,k),QRM(i,k),1./NRM(i,k),icmr,thrd)
3326	!Drop-size limiter [prevents initially large drops from melted hail]
3327	if (Dr>3.e-3) then
3328	tmp1 = (Dr-3.e-3); tmp2= (Dr/DrMAX); tmp3= tmp2tmp2tmp2
3329	NRM(i,k)= NRM(i,k)max((1.+2.e4tmp1*tmp1),tmp3)
3330	tmp1 = DE(i,k)QRM(i,k)icmr
3331	Dr = (tmp1/NRM(i,k))**thrd
3332	endif
3333	else
3334	NRM(i,k)= GR50sqrt(sqrt(GR31iGR34DE(i,k)QRM(i,k)*icmr))
3335	Dr = Dm_x(DE(i,k),QRM(i,k),1./NRM(i,k),icmr,thrd)
3336	endif
3337	iLAMr = iLAMDA_x(DE(i,k),QRM(i,k),1./NRM(i,k),icexr9,thrd)
3338	iLAMr3= iLAMriLAMriLAMr
3339	iLAMr6= iLAMr3*iLAMr3
3340	endif
3341
3342	! Autoconversion:
3343	if (QCM(i,k)>epsQ .and. SIGc>SIGcTHRS .and. autoconv_ON) then
3344	RCAUTR= min( max(L/TAU,0.), QCM(i,k)*idt )
3345	DrINIT= max(83.e-6, 12.6e-4/(0.5e6*SIGc-3.5)) !initiation regime Dr
3346	DrAUT = max(DrINIT, Dr) !init. or feeding DrAUT
3347	CCAUTR= RCAUTRDE(i,k)/(cmrDrAUTDrAUTDrAUT)
3348
3349	! ---------------------------------------------------------------------------- !
3350	! NOTE: The formulation for CCAUTR here (dNr/dt\|initiation) does NOT follow
3351	! eqn (18) in CP2000a, but rather it comes from the F90 code provided
3352	! by J-P Pinty (subroutine: 'rain_c2r2.f90').
3353	! (See notes: 2001-10-17; 2001-10-22)
3354	!
3355	! Similarly, the condition for the activation of accretion and self-
3356	! collection depends on whether or not autoconversion is in the feeding
3357	! regime (see notes 2002-01-07). This is apparent in the F90 code, but
3358	! NOT in CP2000a.
3359	! ---------------------------------------------------------------------------- !
3360
3361	! cloud self-collection: (dNc/dt_autoconversion) {CP eqn(25)}
3362	if (dblMom_c) CCSCOC= min(KK2NCM(i,k)NCM(i,k)GC3iGC1iLAMc6, NCM(i,k) &
3363	idt) !{CP00a eqn(25)}
3364	endif
3365
3366	! Accretion, rain self-collection, and collisional breakup:
3367	if (((QRM(i,k))>1.2max(L,0.)iDE(i,k).or.Dr>max(5.e-6,DrINIT)).and.rainAccr_ON &
3368	.and. rainPresent) then
3369
3370	! Accretion: !{CP00a eqn(22)}
3371	if (QCM(i,k)>epsQ.and.L>0.) then
3372	if (Dr.ge.100.e-6) then
3373	CCACCR = KK1(NCM(i,k)NRM(i,k))(GC2iGC1iLAMc3+GR34iGR31*iLAMr3)
3374	RCACCR = cmriDE(i,k)KK1(NCM(i,k)NRM(i,k))iLAMc3(GC3iGC1iLAMc3+ &
3375	GC2iGC1GR34iGR31iLAMr3)
3376	else
3377	CCACCR = KK2(NCM(i,k)NRM(i,k))(GC3iGC1iLAMc6+GR37iGR31*iLAMr6)
3378
3379	! RCACCR= cmriDE(i,k)KK2(NCM(i,k)NRM(i,k))iLAMc3 &
3380	! (GC4iGR31iLAMc6+GC2iGC1GR37iGR31iLAMr6)
3381	!++ The following calculation of RCACCR avoids overflow:
3382	tmp1 = cmr*iDE(i,k)
3383	tmp2 = KK2(NCM(i,k)NRM(i,k))*iLAMc3
3384	RCACCR = tmp1 * tmp2
3385	tmp1 = GC4*iGR31
3386	tmp1 = (tmp1)*iLAMc6
3387	tmp2 = GC2*iGC1
3388	tmp2 = tmp2GR37iGR31
3389	tmp2 = (tmp2)*iLAMr6
3390	RCACCR = RCACCR * (tmp1 + tmp2)
3391	!++
3392	endif
3393	CCACCR = min(CCACCR,(NC(i,k))*idt)
3394	RCACCR = min(RCACCR,(QC(i,k))*idt)
3395	endif
3396
3397	if (dblMom_r) then
3398	!Rain self-collection:
3399	tmp1= NRM(i,k)*NRM(i,k)
3400	if (Dr.ge.100.e-6) then
3401	CRSCOR= KK1tmp1GR34iGR31iLAMr3 !{CP00a eqn(24)}
3402	else
3403	CRSCOR= KK2tmp1GR37iGR31iLAMr6 !{CP00a eqn(25)}
3404	endif
3405	!Raindrop breakup: !{CP00a eqn(26)}
3406	Ec= 1.
3407	if (Dr >= 600.e-6) Ec= exp(-2.5e3*(Dr-6.e-4))
3408	if (Dr >= 2000.e-6) Ec= 0.
3409	CRSCOR= min(EcCRSCOR,(0.5NR(i,k))*idt) !0.5 prevents depletion of NR
3410	endif
3411
3412	endif !accretion/self-collection/breakup
3413
3414	! Prevent overdepletion of cloud:
3415	source= QC(i,k)
3416	sink = (RCAUTR+RCACCR)*dt
3417	if (sink>source) then
3418	ratio = source/sink
3419	RCAUTR= ratio*RCAUTR
3420	RCACCR= ratio*RCACCR
3421	CCACCR= ratio*CCACCR
3422	endif
3423
3424	! Apply tendencies:
3425	QC(i,k)= max(0., QC(i,k)+(-RCAUTR-RCACCR)*dt )
3426	QR(i,k)= max(0., QR(i,k)+( RCAUTR+RCACCR)*dt )
3427	if (dblMom_c) NC(i,k)= max(0., NC(i,k)+(-CCACCR-CCSCOC)*dt )
3428	if (dblMom_r) NR(i,k)= max(0., NR(i,k)+( CCAUTR-CRSCOR)*dt )
3429
3430	if (dblMom_r) then
3431	if (QR(i,k)>epsQ .and. NR(i,k)>epsN) then
3432	Dr = Dm_x(DE(i,k),QR(i,k),1./NR(i,k),icmr,thrd)
3433	if (Dr>3.e-3) then
3434	tmp1= (Dr-3.e-3); tmp2= tmp1*tmp1
3435	tmp3= (Dr/DrMAX); tmp4= tmp3tmp3tmp3
3436	NR(i,k)= NR(i,k)(max((1.+2.e4tmp2),tmp4))
3437	elseif (Dr<Dhh) then
3438	!Convert small raindrops to cloud:
3439	QC(i,k)= QC(i,k) + QR(i,k)
3440	NC(i,k)= NC(i,k) + NR(i,k)
3441	QR(i,k)= 0.; NR(i,k)= 0.
3442	endif
3443	else
3444	QR(i,k)= 0.; NR(i,k)= 0.
3445	endif !(Qr,Nr>eps)
3446	endif
3447
3448	ENDDO
3449	ENDDO
3450
3451	! Part 3b - Condensation/Evaporation:
3452
3453	DO k=1,nk
3454	DO i=1,ni
3455
3456	DEo = DE(i,nk)
3457	gam = sqrt(DEo*iDE(i,k))
3458	#if (DWORDSIZE == 8 && RWORDSIZE == 8)
3459	QSS(i,k)= FOQSA(T(i,k), PS(i)*S(i,k)) ! Re-calculates QS with new T (w.r.t. liquid)
3460	#elif (DWORDSIZE == 8 && RWORDSIZE == 4)
3461	QSS(i,k)= sngl(FOQSA(T(i,k), PS(i)*S(i,k))) ! Re-calculates QS with new T (w.r.t. liquid)
3462	#else
3463	This is a temporary hack assuming double precision is 8 bytes.
3464	#endif
3465	ssat = Q(i,k)/QSS(i,k)-1.
3466	Tc = T(i,k)-TRPL
3467	Cdiff = max(1.62e-5, (2.2157e-5 + 0.0155e-5Tc)) 1.e5/(S(i,k)*PS(i))
3468	MUdyn = max(1.51e-5, (1.7153e-5 + 0.0050e-5*Tc))
3469	MUkin = MUdyn*iDE(i,k)
3470	iMUkin = 1./MUkin
3471	Ka = max(2.07e-2, (2.3971e-2 + 0.0078e-2*Tc))
3472	ScTHRD = (MUkin/Cdiff)**thrd ! i.e. Sc^(1/3)
3473
3474	!Condensation/evaporation:
3475	! Capacity of evap/cond in one time step is determined by saturation
3476	! adjustment technique [Kong and Yau, 1997 App.A]. Equation for rain evaporation rate
3477	! comes from Cohard and Pinty, 2000a. Explicit condensation rate is not considered
3478	! (as it is in Ziegler, 1985), but rather complete removal of supersaturation is assumed.
3479
3480	X= Q(i,k)-QSS(i,k)
3481	if (dblMom_r) then
3482	rainPresent= (QR(i,k)>epsQ .and. NR(i,k)>epsN)
3483	else
3484	rainPresent= (QR(i,k)>epsQ)
3485	endif
3486	IF(X>0. .or. QC(i,k)>epsQ .or. rainPresent) THEN
3487	tmp1 = T(i,k)-35.86
3488	X = X/(1.+ck5QSS(i,k)/(tmp1tmp1))
3489	if (X<(-QC(i,k))) then
3490	D= 0.
3491	if(rainPresent) then
3492	if(QM(i,k)<QSW(i,k)) then
3493	MUkin = (1.715e-5+5.e-8Tc)iDE(i,k)
3494	iMUkin= 1./MUkin
3495	if (dblMom_r) then
3496	Dr = Dm_x(DE(i,k),QR(i,k),1./NR(i,k),icmr,thrd)
3497	iLAMr= iLAMDA_x(DE(i,k),QR(i,k),1./NR(i,k),icexr9,thrd)
3498	LAMr = 1./iLAMr
3499	!note: The following coding of 'No_r=...' prevents overflow:
3500	!No_r_DM= NR(i,k)LAMr(1.+alpha_r))iGR31
3501	No_r_DM= sngl(dble(NR(i,k))dble(LAMr)dble(1.+alpha_r))iGR31
3502	No_r = No_r_DM
3503	else
3504	iLAMr = sqrt(sqrt( (QR(i,k)DE(i,k))/(GR34cmr*No_r) ))
3505	!note: No_r= No_r_SM is already done (in Part 1)
3506	endif
3507	!note: There is an error in MY05a_eq(8) for VENTx (corrected in code)
3508	VENTr= AvxGR32iLAMr*cexr5 + BvxScTHRDsqrt(gamafriMUkin)GR17* &
3509	iLAMr**cexr6
3510	ABw = CHLCCHLC/(KaRGASVT(i,k)T(i,k))+1./(DE(i,k)(QSS(i,k)) &
3511	Cdiff)
3512	QREVP= -dt(PI2ssatNo_rVENTr/ABw)
3513	!! QREVP= 0. !to suppress evaporation of rain
3514	if ((QR(i,k))>QREVP) then !Note: QREVP is [(dQ/dt)*dt]
3515	DEL= -QREVP
3516	else
3517	DEL= -QR(i,k)
3518	endif
3519	D= max(X+QC(i,k), DEL)
3520	endif !QM< QSM
3521	endif !QR<eps & NR<eps
3522	X= D - QC(i,k)
3523
3524	QR(i,k)= QR(i,k) + D
3525	if (QR(i,k)>0. .and. dblMom_r) &
3526	NR(i,k)= max(0.,NR(i,k)+D*NR(i,k)/QR(i,k)) !(dNr/dt)\|evap
3527	! The above expression of (dNr/dt)\|evap is from Ferrier, 1994.
3528	! In CP2000a, Nr is not affected by evap. (except if Qr goes to zero).
3529	QC(i,k)= 0.; NC(i,k)= 0.
3530	T(i,k) = T(i,k) + LCP*X
3531	Q(i,k) = Q(i,k) - X
3532
3533	else ![if(X >= -QC)]
3534
3535	! Nucleation of cloud droplets:
3536	if (ssat>0. .and. WZ(i,k)>0. .and. dblMom_c) &
3537	NC(i,k)= max(NC(i,k),NccnFNC(WZ(i,k),TM(i,k),HPS(i)*S(i,k),CCNtype))
3538
3539	! All supersaturation is removed (condensed onto cloud field).
3540	T(i,k) = T(i,k) + LCP*X
3541	Q(i,k) = Q(i,k) - X
3542	QC(i,k) = QC(i,k) + X
3543	if (dblMom_c) then
3544	if (X<0.) then
3545	if (QC(i,k)>0.) then
3546	NC(i,k)= max(0., NC(i,k) + X*NC(i,k)/QC(i,k) ) !(dNc/dt)\|evap
3547	else
3548	NC(i,k)= 0.
3549	endif
3550	endif
3551	if (QC(i,k)>0..and.NC(i,k)==0.) NC(i,k)= 1.e7 !prevents non-zero_Q & zero_N
3552	endif
3553
3554	endif
3555
3556	ENDIF
3557
3558	!Protect against negative values due to overdepletion:
3559	if (dblMom_r) then
3560	if (QR(i,k)<epsQ.or.NR(i,k)<epsN) then
3561	Q(i,k) = Q(i,k) + QR(i,k)
3562	T(i,k) = T(i,k) - QR(i,k)*LCP
3563	QR(i,k)= 0.; NR(i,k)= 0.
3564	endif
3565	else
3566	if (QR(i,k)<epsQ) then
3567	Q(i,k) = Q(i,k) + QR(i,k)
3568	T(i,k) = T(i,k) - QR(i,k)*LCP
3569	QR(i,k)= 0.
3570	endif
3571	endif
3572
3573	ENDDO
3574	ENDDO !cond/evap [k-loop]
3575
3576	ENDIF !if warmphase_ON
3577
3578	!----------------------------------------------------------------------------------!
3579	! End of warm-phase microphysics (Part 3) !
3580	!----------------------------------------------------------------------------------!
3581
3582	!----------------------------------------------------------------------------------!
3583	! PART 4: Sedimentation !
3584	!----------------------------------------------------------------------------------!
3585
3586	!----------------------------------------------------------------------------------!
3587	! Sedimentation is computed using a modified version of the box-Lagrangian !
3588	! scheme. Sedimentation is only computed for columns containing non-zero !
3589	! hydrometeor quantities (at at least one level). !
3590	!----------------------------------------------------------------------------------!
3591
3592	IF (sedi_ON) THEN
3593
3594	fluxM_r= 0.; fluxM_i= 0.; fluxM_s= 0.; fluxM_g= 0.; fluxM_h= 0.
3595	RT_rn1 = 0.; RT_rn2 = 0.; RT_fr1 = 0.; RT_fr2 = 0.; RT_sn1 = 0.
3596	RT_sn2 = 0.; RT_sn3 = 0.; RT_pe1 = 0.; RT_pe2 = 0.; RT_peL = 0.
3597
3598	!-- RAIN sedimentation:
3599	if (DblMom_r) then
3600	call SEDI_main_2(QR,NR,1,Q,T,DE,iDE,gamfact_r,epsQr_sedi,epsN,afr,bfr,cmr,dmr, &
3601	ckQr1,ckQr2,icexr9,LCP,ni,nk,VrMax,DrMax,dt,DZ,fluxM_r,ktop_sedi, &
3602	GRAV,massFlux3D=massFlux3D_r)
3603	else !if DblMom_r
3604	call SEDI_main_1b(QR,1,T,DE,iDE,gamfact_r,epsQr_sedi,afr,bfr,icmr,dmr,ckQr1, &
3605	icexr9,ni,nk,VrMax,DrMax,dt,DZ,fluxM_r,No_r_SM,ktop_sedi,GRAV, &
3606	massFlux3D=massFlux3D_r)
3607	endif !if DblMom_r
3608
3609
3610	!-- ICE sedimentation:
3611	if (DblMom_i) then
3612	call SEDI_main_2(QI,NY,2,Q,T,DE,iDE,gamfact,epsQi_sedi,epsN,afi,bfi,cmi,dmi,ckQi1, &
3613	ckQi2,ckQi4,LSP,ni,nk,ViMax,DiMax,dt,DZ,fluxM_i,ktop_sedi,GRAV)
3614	else
3615	call SEDI_main_1b(QI,2,T,DE,iDE,gamfact,epsQi_sedi,afi,bfi,icmi,dmi,ckQi1,ckQi4, &
3616	ni,nk,ViMax,DiMax,dt,DZ,fluxM_i,-99.,ktop_sedi,GRAV)
3617	endif
3618
3619
3620	!-- SNOW sedimentation:
3621	if (DblMom_s) then
3622	call SEDI_main_2(QN,NN,3,Q,T,DE,iDE,gamfact,epsQs_sedi,epsN,afs,bfs,cms,dms,ckQs1, &
3623	ckQs2,iGS20,LSP,ni,nk,VsMax,DsMax,dt,DZ,fluxM_s,ktop_sedi,GRAV, &
3624	massFlux3D=massFlux3D_s)
3625	else
3626	call SEDI_main_1b(QN,3,T,DE,iDE,gamfact,epsQs_sedi,afs,bfs,icms,dms,ckQs1,iGS20, &
3627	ni,nk,VsMax,DsMax,dt,DZ,fluxM_s,-99.,ktop_sedi,GRAV,massFlux3D= &
3628	massFlux3D_s)
3629	endif
3630
3631
3632	!-- GRAUPEL sedimentation:
3633	if (DblMom_g) then
3634	call SEDI_main_2(QG,NG,4,Q,T,DE,iDE,gamfact,epsQg_sedi,epsN,afg,bfg,cmg,dmg,ckQg1, &
3635	ckQg2,ckQg4,LSP,ni,nk,VgMax,DgMax,dt,DZ,fluxM_g,ktop_sedi,GRAV)
3636	else
3637	call SEDI_main_1b(QG,4,T,DE,iDE,gamfact,epsQg_sedi,afg,bfg,icmg,dmg,ckQg1,ckQg4, &
3638	ni,nk,VgMax,DgMax,dt,DZ,fluxM_g,No_g_SM,ktop_sedi,GRAV)
3639	endif
3640
3641
3642	!-- HAIL sedimentation:
3643	if (DblMom_h) then
3644	call SEDI_main_2(QH,NH,5,Q,T,DE,iDE,gamfact,epsQh_sedi,epsN,afh,bfh,cmh,dmh,ckQh1, &
3645	ckQh2,ckQh4,LSP,ni,nk,VhMax,DhMax,dt,DZ,fluxM_h,ktop_sedi,GRAV)
3646	else
3647	call SEDI_main_1b(QH,5,T,DE,iDE,gamfact,epsQh_sedi,afh,bfh,icmh,dmh,ckQh1,ckQh4, &
3648	ni,nk,VhMax,DhMax,dt,DZ,fluxM_h,No_h_SM,ktop_sedi,GRAV)
3649	endif
3650
3651	!======= End of sedimentation for each category ========!
3652
3653	!--- Impose constraints on size distribution parameters ---!
3654	do k= 1,nk
3655	do i= 1,ni
3656
3657	!snow:
3658	if (QN(i,k)>epsQ .and. NN(i,k)>epsN) then
3659
3660	!Impose No_s max for snow: (assumes alpha_s=0.)
3661	iLAMs = max( iLAMmin2, iLAMDA_x(DE(i,k),QN(i,k), 1./NN(i,k),iGS20,idms) )
3662	tmp1 = min(NN(i,k)/iLAMs,No_s_max) !min. No_s
3663	NN(i,k)= tmp1*(dms/(1.+dms))(iGS20DE(i,k)QN(i,k))**(1./(1.+dms)) !impose Nos_max
3664
3665	!Impose LAMDAs_min (by increasing LAMDAs if it is below LAMDAs_min2 [2xLAMDAs_min])
3666	iLAMs = max( iLAMmin2, iLAMDA_x(DE(i,k),QN(i,k),1./NN(i,k),iGS20,idms) )
3667	tmp2 = 1./iLAMs !LAMs before adjustment
3668	!adjust value of lamdas_min to be applied:
3669	! This adjusts for multiple corrections (each time step). The factor 0.6 was obtained by
3670	! trial-and-error to ultimately give reasonable LAMs profiles, smooth and with min LAMs~lamdas_min
3671	tmp4 = 0.6*lamdas_min
3672	tmp5 = 2.*tmp4
3673	tmp3 = tmp2 + tmp4(max(0.,tmp5-tmp2)/tmp5)*2. !LAMs after adjustment
3674	tmp3 = max(tmp3, lamdas_min) !final correction
3675	NN(i,k)= NN(i,k)(tmp3iLAMs)**dms !re-compute NN after LAMs adjustment
3676	endif
3677
3678	enddo !i-loop
3679	enddo !k-loop
3680	!===
3681
3682	!Compute melted (liquid-equivalent) volume fluxes [m3 (liquid) m-2 (sfc area) s-1]:
3683	! (note: For other precipitation schemes in RPN-CMC physics, this is computed in 'vkuocon6.ftn')
3684	RT_rn1 = fluxM_r *idew
3685	RT_sn1 = fluxM_i *idew
3686	RT_sn2 = fluxM_s *idew
3687	RT_sn3 = fluxM_g *idew
3688	RT_pe1 = fluxM_h *idew
3689
3690	!----
3691	!Compute sum of solid (unmelted) volume fluxes [m3 (bulk hydrometeor) m-2 (sfc area) s-1]:
3692	!(this is the precipitation rate for UNmelted total snow [i+s+g])
3693
3694	! Note: In 'calcdiagn.ftn', the total solid precipitation (excluding hail) SN is computed
3695	! from the sum of the liq-eq.vol fluxes, tss_sn1 + tss_sn2 + tss_sn3. With the
3696	! accumulation of SND (in 'calcdiag.ftn'), the solid-to-liquid ratio for the total
3697	! accumulated "snow" (i+s+g) can be compute as SND/SN. Likewise, the instantaneous
3698	! solid-to-liquid ratio of solid precipitation is computed (in 'calcdiag.ftn') as
3699	! RS2L = RSND/RSN.
3700
3701	do i= 1,ni
3702
3703	fluxV_i= fluxM_i(i)*idei
3704	fluxV_g= fluxM_g(i)*ideg
3705	!Compute unmelted volume flux for snow:
3706	! note: This is based on the strict calculation of the volume flux, where vol=(pi/6)D^3,
3707	! and remains in the integral calculation Fv = int[V(D)vol(D)N(D)*dD].
3708	! For a constant density (ice and graupel), vol(D) = m(D)/dex, dex comes out of
3709	! integral and Fv_x=Fm_x/dex
3710	! Optimized for alpha_s = 0.
3711	if (QN(i,nk)>epsQ .and. NN(i,nk)>epsN .and. fluxM_s(i)>0.) then
3712	tmp1= 1./iLAMDA_x(DE(i,nk),QN(i,nk),1./NN(i,nk),iGS20,idms) !LAMDA_s
3713	fluxV_s= fluxM_s(i)rfact_FvFmtmp1**(dms-3.)
3714	else
3715	fluxV_s=0.
3716	endif
3717
3718	!total solid unmelted volume flux, before accounting for partial melting:
3719	tmp1= fluxV_i + fluxV_g + fluxV_s
3720
3721	!liquid-fraction of partially-melted solid precipitation:
3722	! The physical premise is that if QR>0, QN+QI+QG>0, and T>0, then QR
3723	! originates from melting and can be used to estimate the liquid portion
3724	! of the partially-melted solid hydrometeor.
3725	tmp2= QR(i,nk) + QI(i,nk) + QN(i,nk) + QG(i,nk)
3726	if (T(i,nk)>TRPL .and. tmp2>epsQ) then
3727	fracLiq= QR(i,nk)/tmp2
3728	else
3729	fracLiq= 0.
3730	endif
3731
3732	!Tend total volume flux towards the liquid-equivalent as the liquid-fraction increases to 1:
3733	tmp3= RT_sn1(i) + RT_sn2(i) + RT_sn3(i) !total liquid-equivalent volume flux (RSN, Fv_sol)
3734	RT_snd(i)= (1.-fracLiq)tmp1 + fracLiqtmp3 !total volume flux with partial melting (RSND, Fvsle_sol)
3735	!Note: Calculation of instantaneous solid-to-liquid ratio [RS2L = RSND/RSN]
3736	! is based on the above quantities and is done on 'calcdiag.ftn'.
3737
3738	enddo !i-loop
3739	!====
3740
3741	!++++
3742	! Diagnose surface precipitation types:
3743	!
3744	! The following involves diagnostic conditions to determine surface precipitation rates
3745	! for various precipitation elements defined in Canadian Meteorological Operational Internship
3746	! Program module 3.1 (plus one for large hail) based on the sedimentation rates of the five
3747	! sedimenting hydrometeor categories.
3748	!
3749	! With the diagnostics shut off (precipDiag_ON=.false.), 5 rates are just the 5 category
3750	! rates, with the other 6 rates just 0. The model output variables will have:
3751	!
3752	! total liquid = RT_rn1 [RAIN]
3753	! total solid = RT_sn1 [ICE] + RT_sn2 [SNOW] + RT_sn3 [GRAUPEL] + RT_pe1 [HAIL]
3754	!
3755	! With the diagnostics on, the 5 sedimentation rates are partitioned into 9 rates,
3756	! with the following model output variable:
3757	!
3758	! total liquid = RT_rn1 [liquid rain] + RT_rn2 [liquid drizzle]
3759	! total solid = RT_fr1 [freezing rain] + RT_fr2 [freezing drizzle] + RT_sn1 [ice crystals] +
3760	! RT_sn2 [snow] + RT_sn3 [graupel] + RT_pe1 [ice pellets] + RT_pe2 [hail]
3761	!
3762	! NOTE: - The above total liquid/solid rates are computed in 'calcdiag.ftn' (as R2/R4).
3763	! - RT_peL [large hail] is a sub-set of RT_pe2 [hail] and is thus not added to the total
3764	! solid precipitation.
3765
3766	! call tmg_start0(98,'mmCalcDIAG')
3767
3768	IF (precipDiag_ON) THEN
3769	DO i= 1,ni
3770
3771	DE(i,nk)= S(i,nk)PS(i)/(RGASDT(i,nk))
3772
3773	!rain vs. drizzle:
3774	if (DblMom_r) then
3775	N_r= NR(i,nk)
3776	else
3777	N_r= (No_rGR31)(3./(4.+alpha_r))(GR31iGR34DE(i,nk)QR(i,nk)icmr)** &
3778	((1.+alpha_r)/(4.+alpha_r)) !i.e. NR = f(No_r,QR)
3779	endif
3780	if (QR(i,nk)>epsQ .and. N_r>epsN) then
3781	Dm_r(i,nk)= (DE(i,nk)icmrQR(i,nk)/N_r)**thrd
3782	if (Dm_r(i,nk)>Dr_large) then !Dr_large is rain/drizzle size threshold
3783	RT_rn2(i)= RT_rn1(i); RT_rn1(i)= 0.
3784	endif
3785	endif
3786
3787	!liquid vs. freezing rain or drizzle:
3788	if (T(i,nk)<TRPL) then
3789	RT_fr1(i)= RT_rn1(i); RT_rn1(i)= 0.
3790	RT_fr2(i)= RT_rn2(i); RT_rn2(i)= 0.
3791	endif
3792
3793	!ice pellets vs. hail:
3794	if (T(i,nk)>(TRPL+5.0)) then
3795	!note: The condition (T_sfc<5C) for ice pellets is a simply proxy for the presence
3796	! of a warm layer aloft, though which falling snow or graupel will melt to rain,
3797	! over a sub-freezinglayer, where the rain will freeze into the 'hail' category
3798	RT_pe2(i)= RT_pe1(i); RT_pe1(i)= 0.
3799	endif
3800
3801	!large hail:
3802	if (QH(i,nk)>epsQ) then
3803	if (DblMom_h) then
3804	N_h= NH(i,nk)
3805	else
3806	N_h= (No_h_SMGH31)(3./(4.+alpha_h))(GH31iGH34DE(i,nk)QH(i,nk) &
3807	icmh)**((1.+alpha_h)/(4.+alpha_h)) !i.e. Nh = f(No_h,Qh)
3808	endif
3809	Dm_h(i,nk)= Dm_x(DE(i,nk),QH(i,nk),1./N_h,icmh,thrd)
3810	if (DM_h(i,nk)>Dh_large) RT_peL(i)= RT_pe2(i)
3811	!note: large hail (RT_peL) is a subset of the total hail (RT_pe2)
3812	endif
3813
3814	ENDDO
3815	ENDIF !if (precipDiag_ON)
3816	!
3817	!++++
3818
3819	ELSE
3820
3821	massFlux3D_r= 0.
3822	massFlux3D_s= 0.
3823
3824	ENDIF ! if (sedi_ON)
3825
3826	where (Q<0.) Q= 0.
3827
3828	!-----------------------------------------------------------------------------------!
3829	! End of sedimentation calculations (Part 4) !
3830	!-----------------------------------------------------------------------------------!
3831
3832
3833	!===================================================================================!
3834	! End of microphysics scheme !
3835	!===================================================================================!
3836
3837	!-----------------------------------------------------------------------------------!
3838	! Calculation of diagnostic output variables: !
3839
3840	IF (calcDiag) THEN
3841
3842	!For reflectivity calculations:
3843	ZEC= minZET
3844	cxr= icmr*icmr !for rain
3845	cxi= 1./fdielecicmricmr !for all frozen categories
3846	Gzr= (6.+alpha_r)(5.+alpha_r)(4.+alpha_r)/((3.+alpha_r)(2.+alpha_r)(1.+alpha_r))
3847	Gzi= (6.+alpha_i)(5.+alpha_i)(4.+alpha_i)/((3.+alpha_i)(2.+alpha_i)(1.+alpha_i))
3848	if (snowSpherical) then !dms=3
3849	Gzs= (6.+alpha_s)(5.+alpha_s)(4.+alpha_s)/((3.+alpha_s)(2.+alpha_s) &
3850	(1.+alpha_s))
3851	else !dms=2
3852	Gzs= (4.+alpha_s)(3.+alpha_s)/((2.+alpha_s)(1.+alpha_s))
3853	endif
3854	Gzg= (6.+alpha_g)(5.+alpha_g)(4.+alpha_g)/((3.+alpha_g)(2.+alpha_g)(1.+alpha_g))
3855	Gzh= (6.+alpha_h)(5.+alpha_h)(4.+alpha_h)/((3.+alpha_h)(2.+alpha_h)(1.+alpha_h))
3856
3857	do k= 1,nk
3858	do i= 1,ni
3859	DE(i,k)= S(i,k)PS(i)/(RGASDT(i,k))
3860	tmp9= DE(i,k)*DE(i,k)
3861
3862	!Compute N_x for single-moment categories:
3863	if (DblMom_c) then
3864	N_c= NC(i,k)
3865	else
3866	N_c= N_c_SM
3867	endif
3868	if (DblMom_r) then
3869	N_r= NR(i,k)
3870	else
3871	N_r= (No_r_SMGR31)(3./(4.+alpha_r))(GR31iGR34DE(i,k)QR(i,k)icmr)** &
3872	((1.+alpha_r)/(4.+alpha_r)) !i.e. NR = f(No_r,QR)
3873	endif
3874	if (DblMom_i) then
3875	N_i= NY(i,k)
3876	else
3877	N_i= N_Cooper(TRPL,T(i,k))
3878	endif
3879	if (DblMom_s) then
3880	N_s= NN(i,k)
3881	else
3882	No_s= Nos_Thompson(TRPL,T(i,k))
3883	N_s = (No_sGS31)(dms/(1.+dms+alpha_s))(GS31iGS34DE(i,k)QN(i,k) &
3884	icms)**((1.+alpha_s)/(1.+dms+alpha_s))
3885	endif
3886	if (DblMom_g) then
3887	N_g= NG(i,k)
3888	else
3889	N_g= (No_g_SMGG31)(3./(4.+alpha_g))(GG31GG34DE(i,k)QG(i,k)icmg)** &
3890	((1.+alpha_g)/(4.+alpha_g)) !i.e. NX = f(No_x,QX)
3891	endif
3892	if (DblMom_h) then
3893	N_h= NH(i,k)
3894	else
3895	N_h= (No_h_SMGH31)(3./(4.+alpha_h))(GH31iGH34DE(i,k)QH(i,k)icmh)** &
3896	((1.+alpha_h)/(4.+alpha_h)) !i.e. NX = f(No_x,QX)
3897	endif
3898
3899	!Total equivalent reflectivity: (units of [dBZ])
3900	tmp1= 0.; tmp2= 0.; tmp3= 0.; tmp4= 0.; tmp5= 0.
3901	if (QR(i,k)>epsQ .and. N_r>epsN) tmp1 = cxrGzrtmp9QR(i,k)QR(i,k)/N_r
3902	if (QI(i,k)>epsQ .and. N_i>epsN) tmp2 = cxiGzitmp9QI(i,k)QI(i,k)/N_i
3903	if (QN(i,k)>epsQ .and. N_s>epsN) tmp3 = cxiGzstmp9QN(i,k)QN(i,k)/N_s
3904	if (QG(i,k)>epsQ .and. N_g>epsN) tmp4 = cxiGzgtmp9QG(i,k)QG(i,k)/N_g
3905	if (QH(i,k)>epsQ .and. N_h>epsN) tmp5 = cxiGzhtmp9QH(i,k)QH(i,k)/N_h
3906	!Modifiy dielectric constant for melting ice-phase categories:
3907	if ( T(i,k)>TRPL) then
3908	tmp2= tmp2*fdielec
3909	tmp3= tmp3*fdielec
3910	tmp4= tmp4*fdielec
3911	tmp5= tmp5*fdielec
3912	endif
3913	ZET(i,k) = tmp1 + tmp2 + tmp3 + tmp4 + tmp5 != Zr+Zi+Zs+Zg+Zh
3914	if (ZET(i,k)>0.) then
3915	ZET(i,k)= 10.log10((ZET(i,k)Zfact)) !convert to dBZ
3916	else
3917	ZET(i,k)= minZET
3918	endif
3919	ZET(i,k)= max(ZET(i,k),minZET)
3920	ZEC(i)= max(ZEC(i),ZET(i,k)) !composite (max in column) of ZET
3921
3922	!Mean-mass diameters: (units of [m])
3923	Dm_c(i,k)= 0.; Dm_r(i,k)= 0.; Dm_i(i,k)= 0.
3924	Dm_s(i,k)= 0.; Dm_g(i,k)= 0.; Dm_h(i,k)= 0.
3925	if(QC(i,k)>epsQ.and.N_c>epsN) Dm_c(i,k)=Dm_x(DE(i,k),QC(i,k),1./N_c,icmr,thrd)
3926	if(QR(i,k)>epsQ.and.N_r>epsN) Dm_r(i,k)=Dm_x(DE(i,k),QR(i,k),1./N_r,icmr,thrd)
3927	if(QI(i,k)>epsQ.and.N_i>epsN) Dm_i(i,k)=Dm_x(DE(i,k),QI(i,k),1./N_i,icmi,thrd)
3928	if(QN(i,k)>epsQ.and.N_s>epsN) Dm_s(i,k)=Dm_x(DE(i,k),QN(i,k),1./N_s,icms,idms)
3929	if(QG(i,k)>epsQ.and.N_g>epsN) Dm_g(i,k)=Dm_x(DE(i,k),QG(i,k),1./N_g,icmg,thrd)
3930	if(QH(i,k)>epsQ.and.N_h>epsN) Dm_h(i,k)=Dm_x(DE(i,k),QH(i,k),1./N_h,icmh,thrd)
3931
3932	!Supercooled liquid water:
3933	SLW(i,k)= 0.
3934	if (T(i,k)<TRPL) SLW(i,k)= DE(i,k)*(QC(i,k)+QR(i,k)) !(units of [kg/m3])
3935
3936	!Visibility:
3937	!VIS1: component through liquid cloud (fog) [m]
3938	! (following parameterization of Gultepe and Milbrandt, 2007)
3939	tmp1= QC(i,k)DE(i,k)1.e+3 !LWC [g m-3]
3940	tmp2= N_c*1.e-6 !Nc [cm-3]
3941	if (tmp1>0.005 .and. tmp2>1.) then
3942	VIS1(i,k)= max(epsVIS,1000.(1.13(tmp1tmp2)*(-0.51))) !based on FRAM [GM2007, eqn (4)
3943	!VIS1(i,k)= max(epsVIS,min(maxVIS, (tmp1tmp2)*(-0.65))) !based on RACE [GM2007, eqn (3)
3944	else
3945	VIS1(i,k)= 3.*maxVIS !gets set to maxVIS after calc. of VIS
3946	endif
3947
3948	!VIS2: component through rain !based on Gultepe and Milbrandt, 2008, Table 2 eqn (1)
3949	tmp1= massFlux3D_r(i,k)idew3.6e+6 !rain rate [mm h-1]
3950	if (tmp1>0.01) then
3951	VIS2(i,k)= max(epsVIS,1000.(-4.12tmp1**0.176+9.01)) ![m]
3952	else
3953	VIS2(i,k)= 3.*maxVIS
3954	endif
3955
3956	!VIS3: component through snow !based on Gultepe and Milbrandt, 2008, Table 2 eqn (6)
3957	tmp1= massFlux3D_s(i,k)idew3.6e+6 !snow rate, liq-eq [mm h-1]
3958	if (tmp1>0.01) then
3959	VIS3(i,k)= max(epsVIS,1000.(1.10tmp1**(-0.701))) ![m]
3960	else
3961	VIS3(i,k)= 3.*maxVIS
3962	endif
3963
3964	!VIS: visibility due to reduction from all components 1, 2, and 3
3965	! (based on sum of extinction coefficients and Koschmieders's Law)
3966	VIS(i,k) = min(maxVIS, 1./(1./VIS1(i,k) + 1./VIS2(i,k) + 1./VIS3(i,k)))
3967	VIS1(i,k)= min(maxVIS, VIS1(i,k))
3968	VIS2(i,k)= min(maxVIS, VIS2(i,k))
3969	VIS3(i,k)= min(maxVIS, VIS3(i,k))
3970
3971	enddo !i-loop
3972	enddo !k-loop
3973
3974	!Diagnostic levels:
3975	h_CB = noVal_h_XX !height (AGL) of cloud base
3976	h_SN = noVal_h_XX !height (AGL) of snow level [conventional snow (not just QN>0.)]
3977	h_ML1= noVal_h_XX !height (AGL) of melting level [first 0C isotherm from ground]
3978	h_ML2= noVal_h_XX !height (AGL) of melting level [first 0C isotherm from top]
3979	! note: h_ML2 = h_ML1 implies only 1 melting level
3980	tmp1= 1./GRAV
3981	do i= 1,ni
3982	CB_found= .false.; SN_found= .false.; ML_found= .false.
3983	do k= nk,2,-1
3984	!cloud base:
3985	if ((QC(i,k)>epsQ2.or.QI(i,k)>epsQ2) .and. .not.CB_found) then
3986	h_CB(i) = GZ(i,k)*tmp1
3987	CB_found= .true.
3988	endif
3989	!snow level:
3990	if ( ((QN(i,k)>epsQ2 .and. Dm_s(i,k)>minSnowSize) .or. &
3991	(QG(i,k)>epsQ2 .and. Dm_g(i,k)>minSnowSize)) .and. .not.SN_found) then
3992	h_SN(i) = GZ(i,k)*tmp1
3993	SN_found= .true.
3994	endif
3995	!melting level: (height of lowest 0C isotherm)
3996	if (T(i,k)>TRPL .and. T(i,k-1)<TRPL .and. .not.ML_found) then
3997	h_ML1(i) = GZ(i,k)*tmp1
3998	ML_found= .true.
3999	endif
4000	enddo
4001	enddo
4002	do i= 1,ni
4003	ML_found= .false. !from top
4004	do k= 2,nk
4005	!melting level: (height of highest 0C isotherm)
4006	if (T(i,k)>TRPL .and. T(i,k-1)<TRPL .and. .not.ML_found) then
4007	h_ML2(i) = GZ(i,k)*tmp1
4008	ML_found= .true.
4009	endif
4010	enddo
4011	enddo
4012
4013	ENDIF
4014	! !
4015
4016	!-------------
4017	!Convert N from #/m3 to #/kg:
4018	! note: - at this point, NX is updated NX (at t+1); NXTEND is NX before S/S (at t*)
4019	! - NXM is no longer used (it does not need a unit conversion)
4020	do k= 1,nk
4021	DE(:,k) = S(:,k)PS(:)/(RGASDT(:,k)) !air density at time (t)
4022	iDE(:,k)= 1./DE(:,k)
4023	enddo
4024	NC= NCiDE; NCTEND= NCTENDiDE
4025	NR= NRiDE; NRTEND= NRTENDiDE
4026	NY= NYiDE; NYTEND= NYTENDiDE
4027	NN= NNiDE; NNTEND= NNTENDiDE
4028	NG= NGiDE; NGTEND= NGTENDiDE
4029	NH= NHiDE; NHTEND= NHTENDiDE
4030	!=============
4031
4032	!-----------------------------------------------------------------------------------!
4033	! Compute the tendencies of T, Q, QC, etc. (to be passed back to model dynamics) !
4034	! and reset the fields to their initial (saved) values at time {*}: !
4035
4036	do k= 1,nk
4037	do i= 1,ni
4038
4039	tmp1=T_TEND(i,k); T_TEND(i,k)=(T(i,k) -T_TEND(i,k))*iDT; T(i,k) = tmp1
4040	tmp1=Q_TEND(i,k); Q_TEND(i,k)=(Q(i,k) -Q_TEND(i,k))*iDT; Q(i,k) = tmp1
4041	tmp1=QCTEND(i,k); QCTEND(i,k)=(QC(i,k)-QCTEND(i,k))*iDT; QC(i,k)= tmp1
4042	tmp1=QRTEND(i,k); QRTEND(i,k)=(QR(i,k)-QRTEND(i,k))*iDT; QR(i,k)= tmp1
4043	tmp1=QITEND(i,k); QITEND(i,k)=(QI(i,k)-QITEND(i,k))*iDT; QI(i,k)= tmp1
4044	tmp1=QNTEND(i,k); QNTEND(i,k)=(QN(i,k)-QNTEND(i,k))*iDT; QN(i,k)= tmp1
4045	tmp1=QGTEND(i,k); QGTEND(i,k)=(QG(i,k)-QGTEND(i,k))*iDT; QG(i,k)= tmp1
4046	tmp1=QHTEND(i,k); QHTEND(i,k)=(QH(i,k)-QHTEND(i,k))*iDT; QH(i,k)= tmp1
4047
4048	if (DblMom_c) then
4049	tmp1=NCTEND(i,k); NCTEND(i,k)=(NC(i,k)-NCTEND(i,k))*iDT; NC(i,k)= tmp1
4050	endif
4051	if (DblMom_r) then
4052	tmp1=NRTEND(i,k); NRTEND(i,k)=(NR(i,k)-NRTEND(i,k))*iDT; NR(i,k)= tmp1
4053	endif
4054	if (DblMom_i) then
4055	tmp1=NYTEND(i,k); NYTEND(i,k)=(NY(i,k)-NYTEND(i,k))*iDT; NY(i,k)= tmp1
4056	endif
4057	if (DblMom_s) then
4058	tmp1=NNTEND(i,k); NNTEND(i,k)=(NN(i,k)-NNTEND(i,k))*iDT; NN(i,k)= tmp1
4059	endif
4060	if (DblMom_g) then
4061	tmp1=NGTEND(i,k); NGTEND(i,k)=(NG(i,k)-NGTEND(i,k))*iDT; NG(i,k)= tmp1
4062	endif
4063	if (DblMom_h) then
4064	tmp1=NHTEND(i,k); NHTEND(i,k)=(NH(i,k)-NHTEND(i,k))*iDT; NH(i,k)= tmp1
4065	endif
4066
4067	enddo
4068	enddo
4069	! !
4070	!-----------------------------------------------------------------------------------!
4071	END SUBROUTINE mp_milbrandt2mom_main
4072	!___________________________________________________________________________________!
4073
4074	real function des_OF_Ds(Ds_local,desMax_local,eds_local,fds_local)
4075	!Computes density of equivalent-volume snow particle based on (pi/6des)Ds^3 = cms*Ds^dms
4076	real :: Ds_local,desMax_local,eds_local,fds_local
4077	! des_OF_Ds= min(desMax_local, eds_localDs_local*fds_local)
4078	des_OF_Ds= min(desMax_local, eds_localexp(fds_locallog(Ds_local))) !IBM optimization
4079	end function des_OF_Ds
4080
4081
4082	real function Dm_x(DE_local,QX_local,iNX_local,icmx_local,idmx_local)
4083	!Computes mean-mass diameter
4084	real :: DE_local,QX_local,iNX_local,icmx_local,idmx_local
4085	!Dm_x = (DE_localQX_localiNX_localicmx_local)*idmx_local
4086	Dm_x = exp(idmx_locallog(DE_localQX_localiNX_localicmx_local)) !IBM optimization
4087	end function Dm_x
4088
4089
4090	real function iLAMDA_x(DE_local,QX_local,iNX_local,icex_local,idmx_local)
4091	!Computes 1/LAMDA ("slope" parameter):
4092	real :: DE_local,QX_local,iNX_local,icex_local,idmx_local
4093	!iLAMDA_x = (DE_localQX_localiNX_localicex_local)*idmx_local
4094	iLAMDA_x = exp(idmx_locallog(DE_localQX_localiNX_localicex_local)) !IBM optimization
4095	end function
4096
4097
4098	real function N_Cooper(TRPL_local,T_local)
4099	!Computes total number concentration of ice as a function of temperature
4100	!according to parameterization of Cooper (1986):
4101	real :: TRPL_local,T_local
4102	N_Cooper= 5.exp(0.304(TRPL_local-max(233.,T_local)))
4103	end function N_Cooper
4104
4105	real function Nos_Thompson(TRPL_local,T_local)
4106	!Computes the snow intercept, No_s, as a function of temperature
4107	!according to the parameterization of Thompson et al. (2004):
4108	real :: TRPL_local,T_local
4109	Nos_Thompson= min(2.e+8, 2.e+6exp(-0.12min(-0.001,T_local-TRPL_local)))
4110	end function Nos_Thompson
4111
4112	!===================================================================================================!
4113
4114	END MODULE my_dmom_mod
4115
4116	!________________________________________________________________________________________!
4117
4118	MODULE module_mp_milbrandt2mom
4119
4120	use module_wrf_error
4121	use my_dmom_mod
4122
4123	implicit none
4124
4125	! To be done later. Currently, parameters are initialized in the main routine
4126	! (at every time step).
4127
4128	CONTAINS
4129
4130	!----------------------------------------------------------------------------------------!
4131	SUBROUTINE milbrandt2mom_init
4132
4133	! To be done later. Currently, parameters are initialized in the main routine (at every time step).
4134
4135	END SUBROUTINE milbrandt2mom_init
4136
4137	!----------------------------------------------------------------------------------------!
4138
4139	!+---------------------------------------------------------------------+
4140	! This is a wrapper routine designed to transfer values from 3D to 2D. !
4141	!+---------------------------------------------------------------------+
4142
4143
4144	SUBROUTINE mp_milbrandt2mom_driver(qv, qc, qr, qi, qs, qg, qh, nc, nr, ni, ns, ng, &
4145	nh, th, pii, p, w, dz, dt_in, itimestep, &
4146	RAINNC, RAINNCV, SNOWNC, SNOWNCV, GRPLNC, GRPLNCV, &
4147	! HAILNC, HAILNCV, SR, Zet, ccntype, &
4148	HAILNC, HAILNCV, SR, Zet, &
4149	ids,ide, jds,jde, kds,kde, & ! domain dims
4150	ims,ime, jms,jme, kms,kme, & ! memory dims
4151	its,ite, jts,jte, kts,kte) ! tile dims
4152
4153	implicit none
4154
4155	!Subroutine arguments:
4156	integer, intent(in):: ids,ide, jds,jde, kds,kde, &
4157	ims,ime, jms,jme, kms,kme, &
4158	its,ite, jts,jte, kts,kte
4159	real, dimension(ims:ime, kms:kme, jms:jme), intent(inout):: &
4160	qv,qc,qr,qi,qs,qg,qh,nc,nr,ni,ns,ng,nh,th,Zet
4161	real, dimension(ims:ime, kms:kme, jms:jme), intent(in):: &
4162	pii,p,w,dz
4163	real, dimension(ims:ime, jms:jme), intent(inout):: &
4164	RAINNC,RAINNCV,SNOWNC,SNOWNCV,GRPLNC,GRPLNCV,HAILNC,HAILNCV, &
4165	SR
4166	real, intent(in):: dt_in
4167	integer, intent(in):: itimestep !, ccntype
4168
4169	!Local variables:
4170	real, dimension(1:ite-its+1,1:kte-kts+1) :: t2d,qv2d,qc2d,qr2d,qi2d,qs2d,qg2d,qh2d,&
4171	nc2d,nr2d,ni2d,ns2d,ng2d,nh2d,p2d,dz2d,rho,irho,omega2d,t2d_m,qv2d_m,qc2d_m, &
4172	qr2d_m,qi2d_m,qs2d_m,qg2d_m,qh2d_m,nc2d_m,nr2d_m,ni2d_m,ns2d_m,ng2d_m,nh2d_m,&
4173	sigma2d,tmp01,tmp02,tmp03,tmp04,tmp05,tmp06,tmp07,tmp08,tmp09,tmp10,tmp11, &
4174	tmp12,tmp13,tmp14,tmp15,tmp16,tmp17,tmp18,gz2d,zet2d
4175
4176	!tentatively local; to be passed out as output variables later
4177	real, dimension(1:ite-its+1,1:kte-kts+1) :: Dm_c,Dm_r,Dm_i,Dm_s,Dm_g,Dm_h, &
4178	SLW,VIS,VIS1,VIS2,VIS3,SS01,SS02,SS03,SS04,SS05,SS06,SS07,SS08,SS09,SS10, &
4179	SS11,SS12,SS13,SS14,SS15,SS16,SS17,SS18,SS19,SS20,T_tend,Q_tend,QCtend, &
4180	QRtend,QItend,QStend,QGtend,QHtend,NCtend,NRtend,NItend,NStend,NGtend,NHtend
4181
4182	real, dimension(1:ite-its+1) :: rt_rn1,rt_rn2,rt_fr1,rt_fr2,rt_sn1,rt_sn2,rt_sn3, &
4183	rt_pe1,rt_pe2,rt_peL,rt_snd,ZEC,h_CB,h_ML1,h_ML2,h_SN,p_src
4184
4185	real :: dt,ms2mmstp
4186	real :: qc_max,qr_max,qs_max,qi_max,qg_max,qh_max,nc_max,nr_max,ns_max,ni_max, &
4187	ng_max,nh_max
4188	integer :: i,j,k,i2d,j2d,k2d,i2d_max,k2d_max
4189	integer :: imax_qc, imax_qr, imax_qi, imax_qs, imax_qg, imax_qh
4190	integer :: imax_nc, imax_nr, imax_ni, imax_ns, imax_ng, imax_nh
4191	integer :: jmax_qc, jmax_qr, jmax_qi, jmax_qs, jmax_qg, jmax_qh
4192	integer :: jmax_nc, jmax_nr, jmax_ni, jmax_ns, jmax_ng, jmax_nh
4193	integer :: kmax_qc, kmax_qr, kmax_qi, kmax_qs, kmax_qg, kmax_qh
4194	integer :: kmax_nc, kmax_nr, kmax_ni, kmax_ns, kmax_ng, kmax_nh
4195	integer :: i_start, j_start, i_end, j_end, CCNtype
4196	logical :: precipDiag_ON,sedi_ON,warmphase_ON,autoconv_ON,icephase_ON,snow_ON, &
4197	initN,dblMom_c,dblMom_r,dblMom_i,dblMom_s,dblMom_g,dblMom_h
4198
4199	real, parameter :: ms2mmh = 3.6e+6 !conversion factor for precipitation rates
4200	real, parameter :: R_d = 287.04 !gas constant for dry air
4201	character*512 :: mp_debug
4202	!+---+
4203	i2d_max = ite-its+1
4204	k2d_max = kte-kts+1
4205
4206	dt = dt_in
4207	ms2mmstp = 1.e+3*dt !conversion factor: m/2 to mm/step
4208
4209	!--- temporary initialization (until variables are put as namelist options:
4210	! CCNtype = 1. !maritime --> N_c = 80 cm-3 for dblMom_c = .F.
4211	CCNtype = 2. !continental --> N_c = 200 cm-3 for dblMom_c = .F.
4212
4213	precipDiag_ON = .true.; dblMom_c = .true.
4214	sedi_ON = .true.; dblMom_r = .true.
4215	warmphase_ON = .true.; dblMom_i = .true.
4216	autoconv_ON = .true.; dblMom_s = .true.
4217	icephase_ON = .true.; dblMom_g = .true.
4218	snow_ON = .true.; dblMom_h = .true.
4219	initN = .true.
4220	!---
4221
4222
4223	qc_max = 0.; nc_max = 0.
4224	qr_max = 0.; nr_max = 0.
4225	qi_max = 0.; ni_max = 0.
4226	qs_max = 0.; ns_max = 0.
4227	qg_max = 0.; ng_max = 0.
4228	qh_max = 0.; nh_max = 0.
4229
4230	imax_qc = 0; imax_nc = 0; jmax_qc = 0; jmax_nc = 0; kmax_qc = 0; kmax_nc = 0
4231	imax_qr = 0; imax_nr = 0; jmax_qr = 0; jmax_nr = 0; kmax_qr = 0; kmax_nr = 0
4232	imax_qi = 0; imax_ni = 0; jmax_qi = 0; jmax_ni = 0; kmax_qi = 0; kmax_ni = 0
4233	imax_qs = 0; imax_ns = 0; jmax_qs = 0; jmax_ns = 0; kmax_qs = 0; kmax_ns = 0
4234	imax_qg = 0; imax_ng = 0; jmax_qg = 0; jmax_ng = 0; kmax_qg = 0; kmax_ng = 0
4235	imax_qh = 0; imax_nh = 0; jmax_qh = 0; jmax_nh = 0; kmax_qh = 0; kmax_nh = 0
4236
4237	RAINNCV = 0.
4238	SNOWNCV = 0.
4239	GRPLNCV = 0.
4240	HAILNCV = 0.
4241	SR = 0.
4242
4243	do i = 1, 512
4244	mp_debug(i:i) = char(0)
4245	enddo
4246
4247	j_loop1: do j = jts, jte
4248
4249	j2d = j-jts+1 !index value for 2D arrays, to be passed to main micro scheme
4250
4251	i_loop1: do i = its, ite
4252
4253	i2d = i-its+1 !index value for 2D arrays, to be passed to main micro scheme
4254
4255	!Approximate geopotential:
4256	! (assumes lowest model level is at sea-level; acceptable for purposes of scheme)
4257	gz2d(i2d,kts)= 0.
4258	do k = kts+1, kte
4259	gz2d(i2d,k)= gz2d(i2d,k-1) + dz(i,k,j)*9.81
4260	enddo
4261
4262	k_loop1: do k = kts, kte
4263	k2d = k-kts+1 !index value for 2D arrays, to be passed to main micro scheme
4264
4265	!Note: The 3D number concentration variables (seen by WRF dynamics) are in units of 1/kg.
4266	! However, the 2D variables must be converted to units of 1/m3 (by multiplying by air
4267	! density) before being passed to the main subroutine mp_milbrandtsmom. They are then
4268	! converted back after the call, upon putting them back from 2D to 3D variables.
4269
4270	!Convert 3D to 2D arrays (etc.):
4271	t2d(i2d,k2d) = th(i,k,j)*pii(i,k,j)
4272	p2d(i2d,k2d) = p(i,k,j)
4273	dz2d(i2d,k2d) = dz(i,k,j)
4274	qv2d(i2d,k2d) = qv(i,k,j)
4275	!chen rho(i2d,k2d) = p2d(i2d,k2d)/(R_d*t2d(i2d,k2d))
4276	!chen omega2d(i2d,k2d)= -w(i,k,j)p2d(i2d,k2d)9.81
4277	rho(i2d,k2d) = p2d(i2d,k)/(R_d*t2d(i2d,k))
4278	omega2d(i2d,k2d)= -w(i,k,j)rho(i2d,k2d)9.81
4279	qc2d(i2d,k2d) = qc(i,k,j); nc2d(i2d,k2d) = nc(i,k,j)
4280	qi2d(i2d,k2d) = qi(i,k,j); ni2d(i2d,k2d) = ni(i,k,j)
4281	qr2d(i2d,k2d) = qr(i,k,j); nr2d(i2d,k2d) = nr(i,k,j)
4282	qs2d(i2d,k2d) = qs(i,k,j); ns2d(i2d,k2d) = ns(i,k,j)
4283	qg2d(i2d,k2d) = qg(i,k,j); ng2d(i2d,k2d) = ng(i,k,j)
4284	qh2d(i2d,k2d) = qh(i,k,j); nh2d(i2d,k2d) = nh(i,k,j)
4285	sigma2d(i2d,k2d)= p2d(i2d,k2d)/p2d(i2d,kte-kts+1)
4286
4287	enddo k_loop1
4288
4289	K_loop9: do k= kts, kte
4290	k2d = k-kts+1 !index value for 2D arrays, to be passed to main micro scheme
4291	sigma2d(i2d,k2d)= p2d(i2d,k2d)/p2d(i2d,kte-kts+1)
4292	enddo K_loop9
4293
4294	enddo i_loop1
4295
4296	p_src(:)= p2d(:,k2d_max)
4297
4298	!Flip arrays: (to conform to vertical leveling in GEM)
4299	! Note: This step (and the flipping back) could be avoided by changing the indexing
4300	! in the sedimentation subroutine. It is done this way to allow for directly
4301	! pasting the GEM code directly into this subdriver without having to change
4302	! the code.
4303	tmp01= omega2d; tmp02= t2d; tmp03= qv2d; tmp04= qc2d; tmp05=qr2d; tmp06=qi2d
4304	tmp07= qs2d; tmp08= qg2d; tmp09= qh2d; tmp10= nc2d; tmp11=nr2d; tmp12=ni2d
4305	tmp13= ns2d; tmp14= ng2d; tmp15= nh2d; tmp16= sigma2d; tmp17=dz2d; tmp18=gz2d
4306	do k = kts-1,kte-1
4307	k2d = k-kts+1
4308	omega2d(:,k2d+1)= tmp01(:,k2d_max-k2d)
4309	t2d(:,k2d+1) = tmp02(:,k2d_max-k2d)
4310	qv2d(:,k2d+1) = tmp03(:,k2d_max-k2d)
4311	qc2d(:,k2d+1) = tmp04(:,k2d_max-k2d)
4312	qr2d(:,k2d+1) = tmp05(:,k2d_max-k2d)
4313	qi2d(:,k2d+1) = tmp06(:,k2d_max-k2d)
4314	qs2d(:,k2d+1) = tmp07(:,k2d_max-k2d)
4315	qg2d(:,k2d+1) = tmp08(:,k2d_max-k2d)
4316	qh2d(:,k2d+1) = tmp09(:,k2d_max-k2d)
4317	nc2d(:,k2d+1) = tmp10(:,k2d_max-k2d)
4318	nr2d(:,k2d+1) = tmp11(:,k2d_max-k2d)
4319	ni2d(:,k2d+1) = tmp12(:,k2d_max-k2d)
4320	ns2d(:,k2d+1) = tmp13(:,k2d_max-k2d)
4321	ng2d(:,k2d+1) = tmp14(:,k2d_max-k2d)
4322	nh2d(:,k2d+1) = tmp15(:,k2d_max-k2d)
4323	sigma2d(:,k2d+1)= tmp16(:,k2d_max-k2d)
4324	dz2d(:,k2d+1) = tmp17(:,k2d_max-k2d)
4325	gz2d(:,k2d+1) = tmp18(:,k2d_max-k2d)
4326	enddo
4327
4328	!Copy 2d arrays xx2d to xx2d_m: (to facilitate inclusion of main milbrandt2mom
4329	! subroutine which uses arrays at two different time levels, for GEM model)
4330	t2d_m = t2d; qv2d_m = qv2d
4331	qc2d_m = qc2d; nc2d_m = nc2d
4332	qr2d_m = qr2d; nr2d_m = nr2d
4333	qi2d_m = qi2d; ni2d_m = ni2d
4334	qs2d_m = qs2d; ns2d_m = ns2d
4335	qg2d_m = qg2d; ng2d_m = ng2d
4336	qh2d_m = qh2d; nh2d_m = nh2d
4337	call mp_milbrandt2mom_main(omega2d,t2d,qv2d,qc2d,qr2d,qi2d,qs2d,qg2d,qh2d,nc2d, &
4338	nr2d,ni2d,ns2d,ng2d,nh2d,p_src,t2d_m,qv2d_m,qc2d_m,qr2d_m,qi2d_m,qs2d_m, &
4339	qg2d_m,qh2d_m,nc2d_m,nr2d_m,ni2d_m,ns2d_m,ng2d_m,nh2d_m,p_src,sigma2d, &
4340	rt_rn1,rt_rn2,rt_fr1,rt_fr2,rt_sn1,rt_sn2,rt_sn3,rt_pe1,rt_pe2,rt_peL,rt_snd,&
4341	gz2d,T_tend,Q_tend,QCtend,QRtend,QItend,QStend,QGtend,QHtend,NCtend,NRtend, &
4342	NItend,NStend,NGtend,NHtend,dt,i2d_max,1,k2d_max,j,itimestep,CCNtype,precipDiag_ON,&
4343	sedi_ON,warmphase_ON,autoconv_ON,icephase_ON,snow_ON,initN,dblMom_c,dblMom_r,&
4344	dblMom_i,dblMom_s,dblMom_g,dblMom_h,Dm_c,Dm_r,Dm_i,Dm_s,Dm_g,Dm_h,Zet2d,ZEC, &
4345	SLW,VIS,VIS1,VIS2,VIS3,h_CB,h_ML1,h_ML2,h_SN,SS01,SS02,SS03,SS04,SS05,SS06, &
4346	SS07,SS08,SS09,SS10,SS11,SS12,SS13,SS14,SS15,SS16,SS17,SS18,SS19,SS20)
4347
4348	!Add tendencies:
4349	t2d(:,:) = t2d(:,:) + T_tend(:,:)*dt
4350	qv2d(:,:)= qv2d(:,:) + Q_tend(:,:)*dt
4351	qc2d(:,:)= qc2d(:,:) + QCtend(:,:)dt; nc2d(:,:)= nc2d(:,:) + NCtend(:,:)dt
4352	qr2d(:,:)= qr2d(:,:) + QRtend(:,:)dt; nr2d(:,:)= nr2d(:,:) + NRtend(:,:)dt
4353	qi2d(:,:)= qi2d(:,:) + QItend(:,:)dt; ni2d(:,:)= ni2d(:,:) + NItend(:,:)dt
4354	qs2d(:,:)= qs2d(:,:) + QStend(:,:)dt; ns2d(:,:)= ns2d(:,:) + NStend(:,:)dt
4355	qg2d(:,:)= qg2d(:,:) + QGtend(:,:)dt; ng2d(:,:)= ng2d(:,:) + NGtend(:,:)dt
4356	qh2d(:,:)= qh2d(:,:) + QHtend(:,:)dt; nh2d(:,:)= nh2d(:,:) + NHtend(:,:)dt
4357
4358
4359	!Flip arrays back : (to conform to vertical leveling in WRF)
4360	tmp02= t2d; tmp03= qv2d; tmp04= qc2d; tmp05=qr2d; tmp06=qi2d
4361	tmp07= qs2d; tmp08= qg2d; tmp09= qh2d; tmp10= nc2d; tmp11=nr2d; tmp12=ni2d
4362	tmp13= ns2d; tmp14= ng2d; tmp15= nh2d; tmp16= Zet2d; tmp17=ss01; tmp18=ss02
4363	do k = kts-1,kte-1
4364	k2d = k-kts+1
4365	t2d(:,k2d+1) = tmp02(:,k2d_max-k2d)
4366	qv2d(:,k2d+1) = tmp03(:,k2d_max-k2d)
4367	qc2d(:,k2d+1) = tmp04(:,k2d_max-k2d)
4368	qr2d(:,k2d+1) = tmp05(:,k2d_max-k2d)
4369	qi2d(:,k2d+1) = tmp06(:,k2d_max-k2d)
4370	qs2d(:,k2d+1) = tmp07(:,k2d_max-k2d)
4371	qg2d(:,k2d+1) = tmp08(:,k2d_max-k2d)
4372	qh2d(:,k2d+1) = tmp09(:,k2d_max-k2d)
4373	nc2d(:,k2d+1) = tmp10(:,k2d_max-k2d)
4374	nr2d(:,k2d+1) = tmp11(:,k2d_max-k2d)
4375	ni2d(:,k2d+1) = tmp12(:,k2d_max-k2d)
4376	ns2d(:,k2d+1) = tmp13(:,k2d_max-k2d)
4377	ng2d(:,k2d+1) = tmp14(:,k2d_max-k2d)
4378	nh2d(:,k2d+1) = tmp15(:,k2d_max-k2d)
4379	Zet2d(:,k2d+1) = tmp16(:,k2d_max-k2d)
4380	enddo
4381	i_loop2: do i = its, ite
4382	i2d = i-its+1
4383
4384	!Convert individual precipitation rates (in m/s) to WRF precipitation fields:
4385	! note: RAINNC is not actually "rain"; it is the total precipitation.
4386	! The liquid precipitation is the total multiplied by the liquid fraction,
4387	! --> rain = RAINNC*(1-SR) (done elsewhere in WRF)
4388
4389	RAINNCV(i,j) = (rt_rn1(i2d)+rt_rn2(i2d)+rt_fr1(i2d)+rt_fr2(i2d)+rt_sn1(i2d)+ &
4390	rt_sn2(i2d)+rt_sn3(i2d)+rt_pe1(i2d)+rt_pe2(i2d))*ms2mmstp
4391	SNOWNCV(i,j) = (rt_sn1(i2d) + rt_sn2(i2d))*ms2mmstp
4392	HAILNCV(i,j) = (rt_pe1(i2d) + rt_pe2(i2d))*ms2mmstp
4393	GRPLNCV(i,j) = rt_sn3(i2d) *ms2mmstp
4394	RAINNC(i,j) = RAINNC(i,j) + RAINNCV(i,j)
4395	SNOWNC(i,j) = SNOWNC(i,j) + SNOWNCV(i,j)
4396	HAILNC(i,j) = HAILNC(i,j) + HAILNCV(i,j)
4397	GRPLNC(i,j) = GRPLNC(i,j) + GRPLNCV(i,j)
4398	SR(i,j) = (SNOWNCV(i,j)+HAILNCV(i,j)+GRPLNCV(i,j))/(RAINNCV(i,j)+1.e-12)
4399
4400	k_loop2: do k = kts, kte
4401	k2d = k-kts+1
4402	if(.not.(t2d(i2d,k2d)>=173.) .or. (t2d(i2d,k2d)>1000.)) then
4403	write(6,*)
4404	write(6,) ' Stopping in mp_milbrandt2mom_driver due to unrealistic temperature *'
4405	write(6,*) ' step: ',itimestep
4406	write(6,'(a5,5i5,8e15.5)') 'i,k: ',i,j,k,i2d,k2d,t2d(i2d,k2d),qv2d(i2d,k2d),qc2d(i2d,k2d),qr2d(i2d,k2d), &
4407	qi2d(i2d,k2d),qs2d(i2d,k2d),qg2d(i2d,k2d),qh2d(i2d,k2d)
4408	write(6,*)
4409	stop
4410	endif
4411
4412	!Convert back to 3D arrays (and change units of number concentrations back to kg-1):
4413	th(i,k,j) = t2d(i2d,k2d)/pii(i,k,j)
4414	qv(i,k,j) = qv2d(i2d,k2d)
4415	! irho(i,k) = R_d*t2d(i2d,k2d)/p2d(i2d,k2d)
4416	qc(i,k,j) = qc2d(i2d,k2d); nc(i,k,j) = nc2d(i2d,k2d)
4417	qi(i,k,j) = qi2d(i2d,k2d); ni(i,k,j) = ni2d(i2d,k2d)
4418	qr(i,k,j) = qr2d(i2d,k2d); nr(i,k,j) = nr2d(i2d,k2d)
4419	qs(i,k,j) = qs2d(i2d,k2d); ns(i,k,j) = ns2d(i2d,k2d)
4420	qg(i,k,j) = qg2d(i2d,k2d); ng(i,k,j) = ng2d(i2d,k2d)
4421	qh(i,k,j) = qh2d(i2d,k2d); nh(i,k,j) = nh2d(i2d,k2d)
4422	Zet(i,k,j)= Zet2d(i2d,k2d)
4423
4424	enddo k_loop2
4425	enddo i_loop2
4426
4427	enddo j_loop1
4428
4429	do i = 1, 256
4430	mp_debug(i:i) = char(0)
4431	enddo
4432
4433	END SUBROUTINE mp_milbrandt2mom_driver
4434
4435	!+---+-----------------------------------------------------------------+
4436	!________________________________________________________________________________________!
4437
4438	END MODULE module_mp_milbrandt2mom

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

Download in other formats:

Original Format