| blob | e7b1a3383807ea82acdacc509fd09ee686784e3b |
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 !_______________________________________________________________________________!
15 ! !
16 ! Package version: 2.25.2_beta_04 (internal bookkeeping) !
17 ! Last modified: 2015-02-11 !
18 !_______________________________________________________________________________!
20 MODULE my2_mod
22 private
23 public :: mp_milbrandt2mom_main
24 private :: NccnFNC,SxFNC,gamma,gser,gammln,gammp,cfg,gamminc,polysvp,qsat,check_values
25 private :: sedi_wrapper_2,sedi_1D,count_columns,des_OF_Ds,Dm_x,iLAMDA_x,N_Cooper,Nos_Thompson
27 CONTAINS
29 !==============================================================================!
31 REAL FUNCTION NccnFNC(Win,Tin,Pin,CCNtype)
33 !---------------------------------------------------------------------------!
34 ! This function returns number concentration (activated aerosols) as a
35 ! function of w,T,p, based on polynomial approximations of detailed
36 ! approach using a hypergeometric function, following Cohard and Pinty (2000a).
37 !---------------------------------------------------------------------------!
39 IMPLICIT NONE
41 ! PASSING PARAMETERS:
42 real, intent(in) :: Win, Tin, Pin
43 integer, intent(in) :: CCNtype
45 ! LOCAL PARAMETERS:
46 real :: T,p,x,y,a,b,c,d,e,f,g,h,T2,T3,T4,x2,x3,x4,p2
48 x= log10(Win*100.); x2= x*x; x3= x2*x; x4= x2*x2
49 T= Tin - 273.15; T2= T*T; T3= T2*T; T4= T2*T2
50 p= Pin*0.01; p2= p*p
52 if (CCNtype==1) then !Maritime
54 a= 1.47e-9*T4 -6.944e-8*T3 -9.933e-7*T2 +2.7278e-4*T -6.6853e-4
55 b=-1.41e-8*T4 +6.662e-7*T3 +4.483e-6*T2 -2.0479e-3*T +4.0823e-2
56 c= 5.12e-8*T4 -2.375e-6*T3 +4.268e-6*T2 +3.9681e-3*T -3.2356e-1
57 d=-8.25e-8*T4 +3.629e-6*T3 -4.044e-5*T2 +2.1846e-3*T +9.1227e-1
58 e= 5.02e-8*T4 -1.973e-6*T3 +3.944e-5*T2 -9.0734e-3*T +1.1256e0
59 f= -1.424e-6*p2 +3.631e-3*p -1.986
60 g= -0.0212*x4 +0.1765*x3 -0.3770*x2 -0.2200*x +1.0081
61 h= 2.47e-6*T3 -3.654e-5*T2 +2.3327e-3*T +0.1938
62 y= a*x4 + b*x3 + c*x2 + d*x + e + f*g*h
63 NccnFNC= 10.**min(2.,max(0.,y)) *1.e6 ![m-3]
65 else if (CCNtype==2) then !Continental
67 a= 0.
68 b= 0.
69 c=-2.112e-9*T4 +3.9836e-8*T3 +2.3703e-6*T2 -1.4542e-4*T -0.0698
70 d=-4.210e-8*T4 +5.5745e-7*T3 +1.8460e-5*T2 +9.6078e-4*T +0.7120
71 e= 1.434e-7*T4 -1.6455e-6*T3 -4.3334e-5*T2 -7.6720e-3*T +1.0056
72 f= 1.340e-6*p2 -3.5114e-3*p +1.9453
73 g= 4.226e-3*x4 -5.6012e-3*x3 -8.7846e-2*x2 +2.7435e-2*x +0.9932
74 h= 5.811e-9*T4 +1.5589e-7*T3 -3.8623e-5*T2 +1.4471e-3*T +0.1496
75 y= a*x4 +b*x3 +c*x2 + d*x + e + (f*g*h)
76 NccnFNC= 10.**max(0.,y) *1.e6
78 else
80 print*, '*** STOPPED in MODULE ### NccnFNC *** '
81 print*, ' Parameter CCNtype incorrectly specified'
82 stop
84 endif
86 END FUNCTION NccnFNC
87 !======================================================================!
89 real FUNCTION SxFNC(Win,Tin,Pin,Qsw,Qsi,CCNtype,WRT)
91 !---------------------------------------------------------------------------!
92 ! This function returns the peak supersaturation achieved during ascent with
93 ! activation of CCN aerosols as a function of w,T,p, based on polynomial
94 ! approximations of detailed approach using a hypergeometric function,
95 ! following Cohard and Pinty (2000a).
96 !---------------------------------------------------------------------------!
98 IMPLICIT NONE
100 ! PASSING PARAMETERS:
101 integer, intent(IN) :: WRT
102 integer, intent(IN) :: CCNtype
103 real, intent(IN) :: Win, Tin, Pin, Qsw, Qsi
105 ! LOCAL PARAMETERS:
106 real :: Si,Sw,Qv,T,p,x,a,b,c,d,f,g,h,Pcorr,T2corr,T2,T3,T4,x2,x3,x4,p2
107 real, parameter :: TRPL= 273.15
109 x= log10(max(Win,1.e-20)*100.); x2= x*x; x3= x2*x; x4= x2*x2
110 T= Tin; T2= T*T; T3= T2*T; T4= T2*T2
111 p= Pin*0.01; p2= p*p
113 if (CCNtype==1) then !Maritime
115 a= -5.109e-7*T4 -3.996e-5*T3 -1.066e-3*T2 -1.273e-2*T +0.0659
116 b= 2.014e-6*T4 +1.583e-4*T3 +4.356e-3*T2 +4.943e-2*T -0.1538
117 c= -2.037e-6*T4 -1.625e-4*T3 -4.541e-3*T2 -5.118e-2*T +0.1428
118 d= 3.812e-7*T4 +3.065e-5*T3 +8.795e-4*T2 +9.440e-3*T +6.14e-3
119 f= -2.012e-6*p2 + 4.1913e-3*p - 1.785e0
120 g= 2.832e-1*x3 -5.6990e-1*x2 +5.1105e-1*x -4.1747e-4
121 h= 1.173e-6*T3 +3.2174e-5*T2 -6.8832e-4*T +6.7888e-2
122 Pcorr= f*g*h
123 T2corr= 0.9581-4.449e-3*T-2.016e-4*T2-3.307e-6*T3-1.725e-8*T4
125 else if (CCNtype==2) then !Continental [computed for -35<T<-5C]
127 a= 3.80e-5*T2 +1.65e-4*T +9.88e-2
128 b= -7.38e-5*T2 -2.53e-3*T -3.23e-1
129 c= 8.39e-5*T2 +3.96e-3*T +3.50e-1
130 d= -1.88e-6*T2 -1.33e-3*T -3.73e-2
131 f= -1.9761e-6*p2 + 4.1473e-3*p - 1.771e0
132 g= 0.1539*x4 -0.5575*x3 +0.9262*x2 -0.3498*x -0.1293
133 h=-8.035e-9*T4+3.162e-7*T3+1.029e-5*T2-5.931e-4*T+5.62e-2
134 Pcorr= f*g*h
135 T2corr= 0.98888-5.0525e-4*T-1.7598e-5*T2-8.3308e-8*T3
137 else
139 print*, '*** STOPPED in MODULE ### SxFNC *** '
140 print*, ' Parameter CCNtype incorrectly specified'
141 stop
143 endif
145 Sw= (a*x3 + b*x2 +c*x + d) + Pcorr
146 Sw= 1. + 0.01*Sw
147 Qv= Qsw*Sw
148 Si= Qv/Qsi
149 Si= Si*T2corr
150 if (WRT.eq.1) then
151 SxFNC= Sw
152 else
153 SxFNC= Si
154 endif
155 if (Win.le.0.) SxFNC= 1.
157 END function SxFNC
158 !======================================================================!
160 real FUNCTION gamma(xx)
162 ! Modified from "Numerical Recipes"
164 IMPLICIT NONE
166 ! PASSING PARAMETERS:
167 real, intent(IN) :: xx
169 ! LOCAL PARAMETERS:
170 integer :: j
171 real*8 :: ser,stp,tmp,x,y,cof(6),gammadp
174 SAVE cof,stp
175 DATA cof,stp/76.18009172947146d0,-86.50532032941677d0, &
176 24.01409824083091d0,-1.231739572450155d0,.1208650973866179d-2, &
177 -.5395239384953d-5,2.5066282746310005d0/
178 x=dble(xx)
179 y=x
180 tmp=x+5.5d0
181 tmp=(x+0.5d0)*log(tmp)-tmp
182 ser=1.000000000190015d0
183 ! do j=1,6 !original
184 do j=1,4
185 !!do j=1,3 !gives result to within ~ 3 %
186 y=y+1.d0
187 ser=ser+cof(j)/y
188 enddo
189 gammadp=tmp+log(stp*ser/x)
190 gammadp= exp(gammadp)
192 !!GEM:
193 gamma = gammadp
195 END FUNCTION gamma
196 !======================================================================!
198 SUBROUTINE gser(gamser,a,x,gln)
200 ! USES gammln
202 ! Returns the incomplete gamma function P(a,x) evaluated by its series
203 ! representation as gamser. Also returns GAMMA(a) as gln.
205 implicit none
207 integer :: itmax
208 real :: a,gamser,gln,x,eps
209 parameter (itmax=100, eps=3.e-7)
210 integer :: n
211 real :: ap,de1,summ
213 gln=gammln(a)
214 if(x.le.0.)then
215 gamser=0.
216 return
217 endif
218 ap=a
219 summ=1./a
220 de1=summ
221 do n=1,itmax
222 ap=ap+1.
223 de1=de1*x/ap
224 summ=summ+de1
225 if(abs(de1).lt.abs(summ)*eps) goto 1
226 enddo
227 1 gamser=summ*exp(-x+a*log(x)-gln)
228 return
230 END SUBROUTINE gser
231 !======================================================================!
233 real FUNCTION gammln(xx)
235 ! Returns value of ln(GAMMA(xx)) for xx>0
236 ! (modified from "Numerical Recipes")
238 IMPLICIT NONE
240 ! PASSING PARAMETERS:
241 real, intent(IN) :: xx
243 ! LOCAL PARAMETERS:
244 integer :: j
245 real*8 :: ser,stp,tmp,x,y,cof(6)
247 SAVE cof,stp
248 DATA cof,stp/76.18009172947146d0,-86.50532032941677d0, &
249 24.01409824083091d0,-1.231739572450155d0,.1208650973866179d-2, &
250 -.5395239384953d-5,2.5066282746310005d0/
251 x=dble(xx)
252 y=x
253 tmp=x+5.5d0
254 tmp=(x+0.5d0)*log(tmp)-tmp
255 ser=1.000000000190015d0
256 do j=1,6 !original
257 ! do j=1,4
258 y=y+1.d0
259 ser=ser+cof(j)/y
260 enddo
262 !! GEM:
263 gammln=tmp+log(stp*ser/x)
265 END FUNCTION gammln
266 !======================================================================!
268 real FUNCTION gammp(a,x)
270 ! USES gcf,gser
272 ! Returns the incomplete gamma function P(a,x)
274 implicit none
276 real :: a,x,gammcf,gamser,gln
278 if(x.lt.a+1.)then
279 call gser(gamser,a,x,gln)
280 gammp=gamser
281 else
282 call cfg(gammcf,a,x,gln)
283 gammp=1.-gammcf
284 endif
285 return
287 END FUNCTION gammp
288 !======================================================================!
290 SUBROUTINE cfg(gammcf,a,x,gln)
292 ! USES gammln
294 ! Returns the incomplete gamma function (Q(a,x) evaluated by tis continued fraction
295 ! representation as gammcf. Also returns ln(GAMMA(a)) as gln. ITMAX is the maximum
296 ! allowed number of iterations; EPS is the relative accuracy; FPMIN is a number near
297 ! the smallest representable floating-point number.
299 implicit none
301 integer :: i,itmax
302 real :: a,gammcf,gln,x,eps,fpmin
303 real :: an,b,c,d,de1,h
304 parameter (itmax=100,eps=3.e-7)
306 gln=gammln(a)
307 b=x+1.-a
308 c=1./fpmin
309 d=1./b
310 h=d
311 do i= 1,itmax
312 an=-i*(i-a)
313 b=b+2.
314 d=an*d+b
315 if(abs(d).lt.fpmin)d=fpmin
316 c=b+an/c
317 if(abs(c).lt.fpmin) c=fpmin
318 d=1./d
319 de1=d*c
320 h=h*de1
321 if(abs(de1-1.).lt.eps) goto 1
322 enddo
323 1 gammcf=exp(-x+a*log(x)-gln)*h
324 return
326 END SUBROUTINE cfg
327 !======================================================================!
329 real FUNCTION gamminc(p,xmax)
331 ! USES gammp, gammln
332 ! Returns incomplete gamma function, gamma(p,xmax)= P(p,xmax)*GAMMA(p)
333 real :: p,xmax
334 gamminc= gammp(p,xmax)*exp(gammln(p))
336 end FUNCTION gamminc
337 !======================================================================!
339 real function polysvp(T,TYPE)
341 !--------------------------------------------------------------
342 ! Taken from 'module_mp_morr_two_moment.F' (WRFV3.4)
344 ! COMPUTE SATURATION VAPOR PRESSURE
346 ! POLYSVP RETURNED IN UNITS OF PA.
347 ! T IS INPUT IN UNITS OF K.
348 ! TYPE REFERS TO SATURATION WITH RESPECT TO LIQUID (0) OR ICE (1)
350 ! REPLACE GOFF-GRATCH WITH FASTER FORMULATION FROM FLATAU ET AL. 1992,
351 ! TABLE 4 (RIGHT-HAND COLUMN)
352 !--------------------------------------------------------------
354 IMPLICIT NONE
356 REAL DUM
357 REAL T
358 INTEGER TYPE
359 ! ice
360 real a0i,a1i,a2i,a3i,a4i,a5i,a6i,a7i,a8i
361 data a0i,a1i,a2i,a3i,a4i,a5i,a6i,a7i,a8i /&
362 6.11147274, 0.503160820, 0.188439774e-1, &
363 0.420895665e-3, 0.615021634e-5,0.602588177e-7, &
364 0.385852041e-9, 0.146898966e-11, 0.252751365e-14/
366 ! liquid
367 real a0,a1,a2,a3,a4,a5,a6,a7,a8
369 ! V1.7
370 data a0,a1,a2,a3,a4,a5,a6,a7,a8 /&
371 6.11239921, 0.443987641, 0.142986287e-1, &
372 0.264847430e-3, 0.302950461e-5, 0.206739458e-7, &
373 0.640689451e-10,-0.952447341e-13,-0.976195544e-15/
374 real dt
376 ! ICE
378 IF (TYPE.EQ.1) THEN
380 ! POLYSVP = 10.**(-9.09718*(273.16/T-1.)-3.56654* &
381 ! LOG10(273.16/T)+0.876793*(1.-T/273.16)+ &
382 ! LOG10(6.1071))*100.
385 dt = max(-80.,t-273.16)
386 polysvp = a0i + dt*(a1i+dt*(a2i+dt*(a3i+dt*(a4i+dt*(a5i+dt*(a6i+dt*(a7i+a8i*dt)))))))
387 polysvp = polysvp*100.
389 END IF
391 ! LIQUID
393 IF (TYPE.EQ.0) THEN
395 dt = max(-80.,t-273.16)
396 polysvp = a0 + dt*(a1+dt*(a2+dt*(a3+dt*(a4+dt*(a5+dt*(a6+dt*(a7+a8*dt)))))))
397 polysvp = polysvp*100.
399 ! POLYSVP = 10.**(-7.90298*(373.16/T-1.)+ &
400 ! 5.02808*LOG10(373.16/T)- &
401 ! 1.3816E-7*(10**(11.344*(1.-T/373.16))-1.)+ &
402 ! 8.1328E-3*(10**(-3.49149*(373.16/T-1.))-1.)+ &
403 ! LOG10(1013.246))*100.
405 END IF
407 end function polysvp
409 !==============================================================================!
410 real function qsat(temp,pres,wtype)
412 !-----------------------------------------------------------------------------
413 ! Returns the saturation mixing ratio [kg kg-1], as a function of temperature
414 ! pressure, with respect to liquid water [wtype=0] or ice [wtype=1], by calling
415 ! function POLYSVP to obtain the saturation vapor pressure.
417 ! 2013-08-06
418 !-----------------------------------------------------------------------------
420 implicit none
422 !Calling parameters:
423 real, intent(in) :: temp !temperature [K]
424 real, intent(in) :: pres !pressure [Pa]
425 integer, intent(in) :: wtype !0=liquid water; 1=ice
427 !Local variables:
428 real :: tmp1
430 tmp1 = polysvp(temp,wtype) !esat [Pa], wrt liquid (Flatau formulation)
431 qsat = 0.622*tmp1/(pres-tmp1)
433 end function qsat
435 !==============================================================================!
436 subroutine check_values(Qv,T,Qc,Qr,Qi,Qn,Qg,Qh,Nc,Nr,Ny,Nn,Ng,Nh,epsQ,epsN, &
437 check_consistency,force_abort,source_ind)
439 !-----------------------------------------------------------------------------
440 ! Checks current values of prognotic variables for reasonable values and
441 ! stops and prints values if they are out of specified allowable ranges.
442 !
443 ! 'trapAll = 1' means include trap for inconsistency in hydrometeor moments;
444 ! otherwise, only trap for Q, T, and negative Qx, Nx. This option is here
445 ! to allow for Q<epsQ.and.N>epsN or Q>epsQ.and.N<epsN which can be produced
446 ! at the leading edges due to sedimentation and whose values are accpetable
447 ! since Dmax is imposed in SEDI (so one does not necessarily want to trap
448 ! for inconsistency after sedimentation has been called).
449 !
450 ! The value 'source_ind' indicates the approximate line number in 'my2_main'
451 ! from where 'check_values' was called before it resulted in a trap.
452 !
453 ! 2013-06-21
454 !-----------------------------------------------------------------------------
456 implicit none
458 !Calling parameters:
459 real, dimension(:,:), intent(in) :: Qv,T,Qc,Qr,Qi,Qn,Qg,Qh,Nc,Nr,Ny,Nn,Ng,Nh
460 real, intent(in) :: epsQ,epsN
461 integer, intent(in) :: source_ind
462 logical, intent(in) :: force_abort !.TRUE. = forces abort if value violation is detected
463 logical, intent(in) :: check_consistency !.TRUE. = check for sign consistency between Qx and Nx
465 !Local variables:
466 real, parameter :: T_low = 173.
467 real, parameter :: T_high = 323.
468 real, parameter :: Q_high = 1.e-1
469 real, parameter :: N_high = 1.e+20
470 real :: x_low
471 integer :: i,k,ni,nk
472 logical :: trap
474 if (source_ind == 100) then
475 x_low = -1.e+30 !for call at beginning of main
476 else
477 x_low = 0.
478 endif
480 ni = size(Qv,dim=1)
481 nk = size(Qv,dim=2)
483 trap = .false.
484 do k = 1,nk
485 do i = 1,ni
487 if (.not.(T(i,k)>T_low .and. T(i,k)<T_high)) then
488 print*,'** WARNING IN MICRO **'
489 print*,'** src,i,k,T: ',source_ind,i,k,T(i,k)
490 trap = .true.
491 endif
492 if (.not.(Qv(i,k)>=0. .and. Qv(i,k)<Q_high)) then
493 print*,'** WARNING IN MICRO **'
494 print*,'** src,i,k,Qv (HU): ',source_ind,i,k,Qv(i,k)
495 endif
496 !-- NAN check:
497 if (.not.(Qc(i,k)>=x_low .and. Qc(i,k)<Q_high .and. &
498 Qr(i,k)>=x_low .and. Qr(i,k)<Q_high .and. &
499 Qi(i,k)>=x_low .and. Qi(i,k)<Q_high .and. &
500 Qn(i,k)>=x_low .and. Qn(i,k)<Q_high .and. &
501 Qg(i,k)>=x_low .and. Qg(i,k)<Q_high .and. &
502 Qh(i,k)>=x_low .and. Qh(i,k)<Q_high )) then
503 print*,'** WARNING IN MICRO **'
504 print*, '** src,i,k,Qc,Qr,Qi,Qn,Qg,Qh: ',source_ind,i,k,Qc(i,k),Qr(i,k), &
505 Qi(i,k),Qn(i,k),Qg(i,k),Qh(i,k)
506 trap = .true.
507 endif
508 if (.not.(Nc(i,k)>=x_low .and. Nc(i,k)<N_high .and. &
509 Nr(i,k)>=x_low .and. Nr(i,k)<N_high .and. &
510 Ny(i,k)>=x_low .and. Ny(i,k)<N_high .and. &
511 Nn(i,k)>=x_low .and. Nn(i,k)<N_high .and. &
512 Ng(i,k)>=x_low .and. Ng(i,k)<N_high .and. &
513 Nh(i,k)>=x_low .and. Nh(i,k)<N_high )) then
514 print*,'** WARNING IN MICRO **'
515 print*, '** src,i,k,Nc,Nr,Ny,Nn,Ng,Nh: ',source_ind,i,k,Nc(i,k),Nr(i,k), &
516 Ny(i,k),Nn(i,k),Ng(i,k),Nh(i,k)
517 trap = .true.
518 endif
519 !==
520 if (check_consistency) then
521 if ((Qc(i,k)>epsQ.and.Nc(i,k)<epsN) .or. (Qc(i,k)<epsQ.and.Nc(i,k)>epsN)) then
522 print*,'** WARNING IN MICRO **'
523 print*, '** src,i,k,Qc,Nc: ',source_ind,i,k,Qc(i,k),Nc(i,k)
524 trap = .true.
525 endif
526 if ((Qr(i,k)>epsQ.and.Nr(i,k)<epsN) .or. (Qr(i,k)<epsQ.and.Nr(i,k)>epsN)) then
527 print*,'** WARNING IN MICRO **'
528 print*, '** src,i,k,Qr,Nr: ',source_ind,i,k,Qr(i,k),Nr(i,k)
529 trap = .true.
530 endif
531 if ((Qi(i,k)>epsQ.and.Ny(i,k)<epsN) .or. (Qi(i,k)<epsQ.and.Ny(i,k)>epsN)) then
532 print*,'** WARNING IN MICRO **'
533 print*, '** src,i,k,Qi,Ny: ',source_ind,i,k,Qi(i,k),Ny(i,k)
534 trap = .true.
535 endif
536 if ((Qn(i,k)>epsQ.and.Nn(i,k)<epsN) .or. (Qn(i,k)<epsQ.and.Nn(i,k)>epsN)) then
537 print*,'** WARNING IN MICRO **'
538 print*, '** src,i,k,Qn,Nn: ',source_ind,i,k,Qn(i,k),Nn(i,k)
539 trap = .true.
540 endif
541 if ((Qg(i,k)>epsQ.and.Ng(i,k)<epsN) .or. (Qg(i,k)<epsQ.and.Ng(i,k)>epsN)) then
542 print*,'** WARNING IN MICRO **'
543 print*, '** src,i,k,Qg,Ng: ',source_ind,i,k,Qg(i,k),Ng(i,k)
544 trap = .true.
545 endif
546 if ((Qh(i,k)>epsQ.and.Nh(i,k)<epsN) .or. (Qh(i,k)<epsQ.and.Nh(i,k)>epsN)) then
547 print*,'** WARNING IN MICRO **'
548 print*, '** src,i,k,Qh,Nh: ',source_ind,i,k,Qh(i,k),Nh(i,k)
549 trap = .true.
550 endif
551 endif !if (check_consistency)
552 enddo
553 enddo
555 if (trap .and. force_abort) then
556 print*,'** DEBUG TRAP IN MICRO, s/r CHECK_VALUES -- source: ',source_ind
557 if (source_ind/=100) stop
558 endif
560 end subroutine check_values
563 !=====================================================================================!
564 SUBROUTINE sedi_wrapper_2(QX,NX,cat,epsQ,epsQ_sedi,epsN,dmx,ni,VxMax,DxMax,dt, &
565 massFlux_bot,kdir,kbot,ktop_sedi,GRAV,zheight,nk,DE,iDE,iDP, &
566 DZ,iDZ,gamfact,kount,afx_in,bfx_in,cmx_in,ckQx1_in,ckQx2_in,ckQx4_in)
568 !-------------------------------------------------------------------------------------!
569 ! Wrapper for s/r SEDI, for computation on all vertical levels. Called from MY2_MAIN.
570 !-------------------------------------------------------------------------------------!
572 implicit none
574 ! PASSING PARAMETERS:
575 real, dimension(:,:), intent(inout) :: QX,NX
576 real, dimension(:), intent(out) :: massFlux_bot
577 real, dimension(:,:), intent(in) :: zheight, DE,iDE,iDP,DZ,iDZ,gamfact
578 real, intent(in) :: epsQ,epsQ_sedi,epsN,VxMax,dmx,DxMax,dt,GRAV
579 real, intent(in), optional :: afx_in,bfx_in,cmx_in,ckQx1_in,ckQx2_in,ckQx4_in
580 integer, dimension(:), intent(in) :: ktop_sedi
581 integer, intent(in) :: ni,cat,kbot,kdir,nk,kount
583 ! LOCAL VARIABLES:
584 integer, dimension(size(QX,dim=1)) :: activeColumn,ktop
585 integer :: counter
586 integer :: a,i,k
589 massFlux_bot = 0.
591 ktop = ktop_sedi !(i-array) - for complete column, ktop(:)=1 (GEM) or =nk (WRF)
592 call count_columns(QX,ni,epsQ_sedi,counter,activeColumn,kdir,kbot,ktop)
594 DO a = 1,counter
595 i= activeColumn(a)
596 !From here, all sedi calcs are done for each column i
598 call sedi_1D(QX(i,:),NX(i,:),cat,DE(i,:),iDE(i,:),iDP(i,:),gamfact(i,:),epsQ,epsN, &
599 dmx,VxMax,DxMax,dt,DZ(i,:),iDZ(i,:),massFlux_bot(i),kdir,kbot,ktop(i), &
600 GRAV,afx_in=afx_in,bfx_in=bfx_in,cmx_in=cmx_in,ckQx1_in=ckQx1_in, &
601 ckQx2_in=ckQx2_in,ckQx4_in=ckQx4_in)
603 ENDDO !a-loop
605 END SUBROUTINE sedi_wrapper_2
607 !=====================================================================================!
608 SUBROUTINE sedi_1D(QX1d,NX1d,cat,DE1d,iDE1d,iDP1d,gamfact1d,epsQ,epsN,dmx,VxMax,DxMax, &
609 dt,DZ1d,iDZ1d,massFlux_bot,kdir,kbot,ktop,GRAV,afx_in,bfx_in,cmx_in, &
610 ckQx1_in,ckQx2_in,ckQx4_in,BX1d,epsB)
612 implicit none
614 ! PASSING PARAMETERS:
615 real, dimension(:), intent(inout), optional :: BX1d
616 real, dimension(:), intent(inout) :: QX1d,NX1d
617 real, dimension(:), intent(in) :: gamfact1d
618 real, intent(out) :: massFlux_bot
619 real, dimension(:), intent(in) :: DE1d,iDE1d,iDP1d,DZ1d,iDZ1d
620 real, intent(in) :: epsQ,epsN,VxMax,dmx,DxMax,dt,GRAV
621 real, optional, intent(in) :: afx_in,bfx_in,cmx_in,ckQx1_in,ckQx2_in, &
622 ckQx4_in,epsB
623 integer, intent(in) :: cat,kbot,kdir
624 integer, intent(in) :: ktop
626 ! LOCAL PARAMETERS:
627 integer :: npassx
628 real, dimension(size(QX1d,dim=1)) :: VVQ,VVN
629 real :: dzMIN,dtx,VxMaxx
630 logical :: firstPass,QxPresent,BX_present
631 integer :: nnn,i,k,l,km1,kp1,idzmin,kk
632 real :: VqMax,VnMax,iLAMx,iLAMxB0,tmp1,tmp2,tmp3,Dx, &
633 iDxMax,icmx,VincFact,ratio_Vn2Vq,zmax_Q,zmax_N, &
634 idmx
635 real :: alpha_x,afx,bfx,cmx,ckQx1,ckQx2,ckQx4
637 real, parameter :: thrd = 1./3.
638 real, parameter :: sxth = 1./6.
639 ! real, parameter :: CoMAX = 0.5 !0.8
640 real, parameter :: CoMAX = 0.8
641 real, parameter :: PIov6 = 3.14159265*sxth
642 !-------------------------------------------------------------------------------------!
644 BX_present = present(BX1d)
646 !for rain, ice, snow, hail:
647 ! ! if (.not. (cat==4 .and. BX_present)) then
648 afx = afx_in
649 bfx = bfx_in
650 cmx = cmx_in
651 icmx = 1./cmx
652 ckQx1 = ckQx1_in
653 ckQx2 = ckQx2_in
654 ckQx4 = ckQx4_in
655 ratio_Vn2Vq = ckQx2/ckQx1
656 ! ! endif
658 massFlux_bot = 0.
659 iDxMax = 1./DxMax
660 idmx = 1./dmx
661 VVQ = 0.
662 VVN = 0.
663 VqMax = 0.
664 VnMax = 0.
665 VVQ(:) = 0.
667 do k= kbot,ktop,kdir
668 QxPresent = (QX1d(k)>epsQ .and. NX1d(k)>epsN)
669 if (QxPresent) VVQ(k)= VV_Q()
670 enddo
672 Vxmaxx= min( VxMax, maxval(VVQ(:)))
673 if (kdir==1) then
674 dzMIN = minval(DZ1d) !WRF (to be tested)
675 else
676 dzMIN = minval(DZ1d(ktop:kbot+kdir)) !GEM
677 endif
678 npassx= max(1, nint( dt*Vxmaxx/(CoMAX*dzMIN) ))
681 dtx = dt/float(npassx)
683 !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -!
684 DO nnn= 1,npassx
686 firstPass = (nnn==1)
687 do k= kbot,ktop,kdir
688 QxPresent = (QX1d(k)>epsQ .and. NX1d(k)>epsN)
689 if (QxPresent) then
690 if (firstPass) then !to avoid re-computing VVQ on first pass
691 VVQ(k)= -VVQ(k)
692 else
693 VVQ(k)= -VV_Q()
694 endif
695 !--
696 !control excessive size-sorting for hail:
697 if (cat==5) then
698 tmp1 = (icmx*QX1d(k)/NX1d(k))**thrd !Dmh
699 tmp2 = min(50., 0.1*(1000.*tmp1)) !mu = const*Dmh [mm]
700 ratio_Vn2Vq = ((3.+tmp2)*(2.+tmp2)*(1.+tmp2))/((3.+bfx+tmp2)* &
701 (2.+bfx+tmp2)*(1.+bfx+tmp2))
702 endif
703 !==
704 VVN(k) = VVQ(k)*ratio_Vn2Vq
705 VqMax = max(VxMAX,-VVQ(k))
706 VnMax = max(VxMAX,-VVN(k))
707 else
708 VVQ(k) = 0.
709 VVN(k) = 0.
710 VqMax = 0.
711 VnMax = 0.
712 endif
713 enddo !k-loop
715 massFlux_bot= massFlux_bot - VVQ(kbot)*DE1d(kbot)*QX1d(kbot)
716 do k= kbot,ktop,kdir
717 QX1d(k)= QX1d(k) + dtx*iDE1d(k)*(-DE1d(k+kdir)*QX1d(k+kdir)*VVQ(k+kdir) + &
718 DE1d(k)*QX1d(k)*VVQ(k))*iDZ1d(k+kdir)
719 NX1d(k)= NX1d(k) + dtx*(-NX1d(k+kdir)*VVN(k+kdir) + NX1d(k)*VVN(k))*iDZ1d(k+kdir)
720 QX1d(k) = max( QX1d(k), 0.)
721 NX1d(k) = max( NX1d(k), 0.)
722 enddo
723 !
724 if (BX_present) then
725 do k= kbot,ktop,kdir
726 BX1d(k)= BX1d(k) + dtx*iDE1d(k)*(-DE1d(k+kdir)*BX1d(k+kdir)*VVQ(k+kdir) + &
727 DE1d(k)*BX1d(k)*VVQ(k))*iDZ1d(k+kdir)
728 BX1d(k) = max( BX1d(k), 0.)
729 enddo
730 endif
731 !--
733 do k= kbot,ktop,kdir
735 !Prescribe NX if QX>0.and.NX=0:
736 if (QX1d(k)>epsQ .and. NX1d(k)<epsN) then
737 !determine first level above with NX>epsN
738 do kk = k+kdir,ktop,kdir
739 !note: the following condition should normally be satisfied immediately;
740 ! that is, the next level up should contain NX>epsN
741 if (NX1d(kk)>=epsN) exit
742 enddo
743 !prescribe new NX:
744 ! note: if no kk with NX>epsN found [i.e. if kk=ktop at this point] then
745 ! epsN is prescribed; this will then be modified via size-limiter below
746 NX1d(k) = max(epsN,NX1d(kk))
747 endif
749 !Impose size-limiter / drop-breakup:
750 if (QX1d(k)>epsQ .and. NX1d(k)>epsN) then
751 Dx= (DE1d(k)*QX1d(k)/(NX1d(k)*cmx))**idmx
752 if (cat==1 .and. Dx>3.e-3) then
753 NX1d(k)= NX1d(k)*max((1.+2.e4*(Dx-3.e-3)**2),(Dx*iDxMAX)**3)
754 else
755 NX1d(k)= NX1d(k)*(max(Dx,DxMAX)*iDxMAX)**dmx !impose Dx_max
756 endif
757 endif
759 enddo
761 ENDDO !nnn-loop
762 !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -!
763 !Compute average mass flux during the full time step: (used to compute the
764 !instantaneous sedimentation rate [liq. equiv. volume flux] in the main s/r)
765 massFlux_bot= massFlux_bot/float(npassx)
767 !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -!
769 CONTAINS
771 real function VV_Q()
772 !Calculates Q-weighted fall velocity
773 iLAMx = ((QX1d(k)*DE1d(k)/NX1d(k))*ckQx4)**idmx
774 iLAMxB0 = iLAMx**bfx
775 VV_Q = gamfact1d(k)*iLAMxB0*ckQx1
776 end function VV_Q
778 real function VV_Qg()
779 !Calculates Q-weighted fall velocity
780 ! iLAMx is already computed in 'calc_grpl_params' (for graupel only)
781 iLAMxB0 = iLAMx**bfx
782 VV_Qg = gamfact1d(k)*iLAMxB0*ckQx1
783 end function VV_Qg
785 END SUBROUTINE sedi_1D
787 !=====================================================================================!
788 SUBROUTINE count_columns(QX,ni,minQX,counter,activeColumn,kdir,kbot,ktop)
790 !--------------------------------------------------------------------------
791 ! Searches the hydrometeor array QX(ni,nk) for non-zero (>minQX) values.
792 ! Returns the array if i-indices (activeColumn) for the column indices (i) which
793 ! contain at least one non-zero value, the number of such columns (counter),
794 ! and the k-indices of the maximum level to compute sedimentation.
795 !--------------------------------------------------------------------------
797 implicit none
799 !PASSING PARAMETERS:
800 real, dimension(:,:),intent(in) :: QX ! mixing ratio
801 real, intent(in) :: minQX ! mixing ratio threshold
802 integer, intent(in) :: ni ! total number of columns (input)
803 integer, intent(in) :: kbot ! k index of lowest level
804 integer, intent(in) :: kdir ! -1 of k=1 is top (GEM); 1 if k=1 is bottom (WRF)
805 integer, dimension(:), intent(inout):: ktop ! IN: array of highest level to look at; OUT: array of highest level with QX>epsQ
806 integer, intent(out) :: counter ! number of columns containing at least one QX>epsQ
807 integer, dimension(:), intent(out) :: activeColumn ! array of i-indices with columns containing at least one QX>epsQ
809 !LOCAL PARAMETERS:
810 integer :: i
811 integer, dimension(size(QX,dim=1)) :: k
813 counter = 0
814 activeColumn = 0
816 !Note: k_top(i) must be at least one level higher than the level with non-zero Qx
818 do i=1,ni
819 k(i)= ktop(i)
820 do
821 k(i)=k(i)-kdir !step 1 level downward (towards lowest-level k)
822 if (QX(i,k(i))>minQX) then
823 counter=counter+1
824 activeColumn(counter)=i
825 ktop(i)= k(i) !set ktop(k) to highest level with QX>minQX
826 exit
827 else
828 if (k(i)==kbot) then
829 ktop(i) = kbot
830 exit
831 endif
832 endif
833 enddo
834 enddo
836 END SUBROUTINE count_columns
837 !=======================================================================================!
839 !_______________________________________________________________________________________!
841 SUBROUTINE mp_milbrandt2mom_main(WZ,T,Q,QC,QR,QI,QN,QG,QH,NC,NR,NY,NN,NG,NH,PS, &
842 sigma,RT_rn1,RT_rn2,RT_fr1,RT_fr2,RT_sn1,RT_sn2,RT_sn3,RT_pe1,RT_pe2,RT_peL,RT_snd, &
843 dt,NI,NK,J,KOUNT,CCNtype,precipDiag_ON,sedi_ON,warmphase_ON,autoconv_ON,icephase_ON, &
844 snow_ON,Dm_c,Dm_r,Dm_i,Dm_s,Dm_g,Dm_h,ZET,ZEC,SS,nk_bottom)
846 implicit none
848 !CALLING PARAMETERS:
849 integer, intent(in) :: NI,NK,J,KOUNT,CCNtype
850 real, intent(in) :: dt
851 real, dimension(:), intent(in) :: PS
852 real, dimension(:), intent(out) :: RT_rn1,RT_rn2,RT_fr1,RT_fr2,RT_sn1,RT_sn2, &
853 RT_sn3,RT_pe1,RT_pe2,RT_peL,ZEC,RT_snd
854 real, dimension(:,:), intent(in) :: WZ,sigma
855 real, dimension(:,:), intent(inout) :: T,Q,QC,QR,QI,QN,QG,QH,NC,NR,NY,NN,NG,NH
856 real, dimension(:,:), intent(out) :: ZET,Dm_c,Dm_r,Dm_i,Dm_s,Dm_g,Dm_h
857 real, dimension(:,:,:),intent(out) :: SS
858 logical, intent(in) :: precipDiag_ON,sedi_ON,icephase_ON,snow_ON, &
859 warmphase_ON,autoconv_ON,nk_BOTTOM
861 !_______________________________________________________________________________________
862 ! !
863 ! Milbrandt-Yau Multimoment Bulk Microphysics Scheme !
864 ! - double-moment version - !
865 !_______________________________________________________________________________________!
866 !
867 ! Author:
868 ! J. Milbrandt, McGill University (August 2004)
869 !
870 ! Major revisions:
871 !
872 ! 001 J. Milbrandt (Dec 2006) - Converted the full Milbrandt-Yau (2005) multimoment
873 ! (RPN) scheme to an efficient fixed-dispersion double-moment
874 ! version
875 ! 002 J. Milbrandt (Mar 2007) - Added options for single-moment/double-moment for
876 ! each hydrometeor category
877 ! 003 J. Milbrandt (Feb 2008) - Modified single-moment version for use in GEM-LAM-2.5
878 ! 004 J. Milbrandt (Nov 2008) - Modified double-moment version for use in 2010 Vancouver
879 ! Olympics GEM-LAM configuration
880 ! 005 J. Milbrandt (Aug 2009) - Modified (dmom) for PHY_v5.0.4, for use in V2010 system:
881 ! + reduced ice/snow capacitance to C=0.25D (from C=0.5D)
882 ! + added diagnostic fields (VIS, levels, etc.)
883 ! + added constraints to snow size distribution (No_s and
884 ! LAMDA_s limits, plus changed m-D parameters
885 ! + modified solid-to-liquid ratio calculation, based on
886 ! volume flux (and other changes)
887 ! + added back sedimentation of ice category
888 ! + modified condition for conversion of graupel to hail
889 ! + corrected bug it diagnostic "ice pellets" vs. "hail"
890 ! + minor bug corrections (uninitialized values, etc.)
891 ! 006 J. Milbrandt (Jan 2011) - Bug fixes and minor code clean-up from PHY_v5.1.3 version
892 ! + corrected latent heat constants in thermodynamic functions
893 ! (ABi and ABw) for sublimation and evaporation
894 ! + properly initialized variables No_g and No_h
895 ! + changed max ice crystal size (fallspeed) to 5 mm (2 m s-1)
896 ! + imposed maximum ice number concentration of 1.e+7 m-3
897 ! + removed unused supersaturation reduction
898 !
899 ! Object:
900 ! Computes changes to the temperature, water vapor mixing ratio, and the
901 ! mixing ratios and total number concentrations of six hydrometeor species
902 ! resulting from cloud microphysical interactions at saturated grid points.
903 ! Liquid and solid surface precipitation rates from sedimenting hydrometeor
904 ! categories are also computed.
905 !
906 ! This subroutine and the associated modules form the single/double-moment
907 ! switchable verion of the multimoment bulk microphysics package, the full
908 ! version of which is described in the references below.
909 !
910 ! References: Milbrandt and Yau, (2005a), J. Atmos. Sci., vol.62, 3051-3064
911 ! --------- and ---, (2005b), J. Atmos. Sci., vol.62, 3065-3081
912 ! (and references therein)
913 !
914 ! Please report bugs to: jason.milbrandt@ec.gc.ca
915 !_______________________________________________________________________________________!
916 !
917 ! Arguments: Description: Units:
918 !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -!
919 ! - Input -
920 !
921 ! NI number of x-dir points (in local subdomain)
922 ! NK number of vertical levels
923 ! N not used (to be removed)
924 ! J y-dir index (local subdomain)
925 ! KOUNT current model time step number
926 ! dt model time step [s]
927 ! CCNtype switch for airmass type
928 ! 1 = maritime --> N_c = 80 cm-3 (1-moment cloud)
929 ! 2 = continental 1 --> N_c = 200 cm-3 " "
930 ! 3 = continental 2 (polluted) --> N_c = 500 cm-3 " "
931 ! 4 = land-sea-mask-dependent (TBA)
932 ! WZ vertical velocity [m s-1]
933 ! sigma sigma = p/p_sfc
934 ! precipDiag_ON logical switch, .F. to suppress calc. of sfc precip types
935 ! sedi_ON logical switch, .F. to suppress sedimentation
936 ! warmphase_ON logical switch, .F. to suppress warm-phase (Part II)
937 ! autoconv_ON logical switch, .F. to supppress autoconversion (cld->rn)
938 ! icephase_ON logical switch, .F. to suppress ice-phase (Part I)
939 ! snow_ON logical switch, .F. to suppress snow initiation
940 ! nk_BOTTOM logical switch, .T. for nk at bottom (GEM, 1dkin); .F. for nk at top (WRF, 2dkin)
941 !
942 !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -!
943 ! - Input/Output -
944 !
945 ! T air temperature at time (t*) [K]
946 ! Q water vapor mixing ratio at (t*) [kg kg-1]
947 ! PS surface pressure at time (t*) [Pa]
948 !
949 ! For x = (C,R,I,N,G,H): C = cloud
950 ! R = rain
951 ! I = ice (pristine) [except 'NY', not 'NI']
952 ! N = snow
953 ! G = graupel
954 ! H = hail
955 !
956 ! Q(x) mixing ratio for hydrometeor x at (t*) [kg kg-1]
957 ! N(x) total number concentration for hydrometeor x (t*) [m-3]
958 !
959 !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -!
960 ! - Output -
961 !
962 ! Dm_(x) mean-mass diameter for hydrometeor x [m]
963 ! RT_rn1 precipitation rate (at sfc) of liquid rain [m+3 m-2 s-1]
964 ! RT_rn2 precipitation rate (at sfc) of liquid drizzle [m+3 m-2 s-1]
965 ! RT_fr1 precipitation rate (at sfc) of freezing rain [m+3 m-2 s-1]
966 ! RT_fr2 precipitation rate (at sfc) of freezing drizzle [m+3 m-2 s-1]
967 ! RT_sn1 precipitation rate (at sfc) of ice crystals (liq-eq) [m+3 m-2 s-1]
968 ! RT_sn2 precipitation rate (at sfc) of snow (liq-equiv) [m+3 m-2 s-1]
969 ! RT_sn3 precipitation rate (at sfc) of graupel (liq-equiv) [m+3 m-2 s-1]
970 ! RT_snd precipitation rate (at sfc) of snow (frozen) [m+3 m-2 s-1]
971 ! RT_pe1 precipitation rate (at sfc) of ice pellets (liq-eq) [m+3 m-2 s-1]
972 ! RT_pe2 precipitation rate (at sfc) of hail (total; liq-eq) [m+3 m-2 s-1]
973 ! RT_peL precipitation rate (at sfc) of hail (large only) [m+3 m-2 s-1]
974 ! SS(i,k,n) array (n) for 3D diagnostic output (e.g. S/S term)
975 ! ZET total equivalent radar reflectivity [dBZ]
976 ! ZEC composite (column-max) of ZET [dBZ]
977 !_______________________________________________________________________________________!
980 !LOCAL VARIABLES:
982 !Variables to count active grid points:
983 logical :: log1,log2,log3,log4,doneK,rainPresent,calcDiag
984 logical, dimension(size(QC,dim=1),size(QC,dim=2)) :: activePoint
985 integer, dimension(size(QC,dim=1)) :: ktop_sedi
986 integer :: i,k,niter,ll,start,kskip_1,ktop,kbot,kdir
988 real :: tmp1,tmp2,tmp3,tmp4,tmp5,tmp6,tmp7,tmp8,tmp9,tmp10, &
989 VDmax,NNUmax,X,D,DEL,QREVP,NuDEPSOR,NuCONTA,NuCONTB,NuCONTC,iMUkin,Ecg,Erg, &
990 NuCONT,GG,Na,Tcc,F1,F2,Kdiff,PSIa,Kn,source,sink,sour,ratio,qvs0,Kstoke, &
991 DELqvs,ft,esi,Si,Simax,Vq,Vn,Vz,LAMr,No_r_DM,No_i,No_s,No_g,No_h,D_sll, &
992 iABi,ABw,VENTr,VENTs,VENTg,VENTi,VENTh,Cdiff,Ka,MUdyn,MUkin,Ng_tail, &
993 gam,ScTHRD,Tc,mi,ff,Ec,Ntr,Dho,DMrain,Ech,DMice,DMsnow,DMgrpl,DMhail, &
994 ssat,Swmax,dey,Esh,Eii,Eis,Ess,Eig,Eih,FRAC,JJ,Dirg,Dirh,Dsrs,Dsrg,Dsrh, &
995 Dgrg,Dgrh,SIGc,L,TAU,DrAUT,DrINIT,Di,Ds,Dg,Dh,qFact,nFact,Ki,Rz,NgCNgh, &
996 vr0,vi0,vs0,vg0,vh0,Dc,Dr,QCLcs,QCLrs,QCLis,QCLcg,QCLrg,QCLig,NhCNgh, &
997 QCLch,QCLrh,QCLsh,QMLir,QMLsr,QMLgr,QMLhr,QCLih,QVDvg,QVDvh,QSHhr, &
998 QFZci,QNUvi,QVDvi,QCNis,QCNis1,QCNis2,QCLir,QCLri,QCNsg,QCLsr,QCNgh, &
999 QCLgr,QHwet,QVDvs,QFZrh,QIMsi,QIMgi,NMLhr,NVDvh,NCLir,NCLri,NCLrh, &
1000 NCLch,NCLsr,NCLirg,NCLirh,NrFZrh,NhFZrh,NCLsrs,NCLsrg,NCLsrh,NCLgrg, &
1001 NCLgrh,NVDvg,NMLgr,NiCNis,NsCNis,NVDvs,NMLsr,NCLsh,NCLss,NNUvi,NFZci,NVDvi, &
1002 NCLis,NCLig,NCLih,NMLir,NCLrs,NCNsg,NCLcs,NCLcg,NIMsi,NIMgi,NCLgr,NCLrg, &
1003 NSHhr,RCAUTR,RCACCR,CCACCR,CCSCOC,CCAUTR,CRSCOR,ALFx,des_pmlt,Ecs,des,ides, &
1004 LAMx,iLAMx,iLAMxB0,Dx,ffx,iLAMc,iNC,iNR,iNY,iNN,iNG,iLAMs_D3, &
1005 iLAMg,iLAMg2,iLAMgB0,iLAMgB1,iLAMgB2,iLAMh,iLAMhB0,iLAMhB1,iLAMhB2,iNH, &
1006 iLAMi,iLAMi2,iLAMi3,iLAMi4,iLAMi5,iLAMiB0,iLAMiB1,iLAMiB2,iLAMr6,iLAMh2, &
1007 iLAMs,iLAMs2,iLAMsB0,iLAMsB1,iLAMsB2,iLAMr,iLAMr2,iLAMr3,iLAMr4,iLAMr5, &
1008 iLAMc2,iLAMc3,iLAMc4,iLAMc5,iLAMc6,iQC,iQR,iQI,iQN,iQG,iQH,iEih,iEsh, &
1009 N_c,N_r,N_i,N_s,N_g,N_h,fluxV_i,fluxV_g,fluxV_s,rhos_mlt,fracLiq
1011 !Variables that only need to be calulated on the first step (and saved):
1012 real, save :: idt,iMUc,cmr,cmi,cms,cmg,cmh,icmr,icmi,icmg,icms,icmh,idew,idei, &
1013 ideh,ideg,GC1,imso,icexc9,cexr1,cexr2,cexr3,No_s_SM,No_r,idms,imgo,icexs2, &
1014 cexr4,cexr5,cexr6,cexr9,icexr9,ckQr1,ckQr2,ckQr3,ckQi1,ckQi2,ckQi3,ckQi4, &
1015 icexi9,ckQs1,ckQs2,cexs1,cexs2,ckQg1,ckQg2,ckQg4,ckQh1,ckQh2,ckQh4,GR37,dms, &
1016 LCP,LFP,LSP,ck5,ck6,PI2,PIov4,PIov6,CHLS,iCHLF,cxr,cxi,Gzr,Gzi,Gzs,Gzg,Gzh, &
1017 N_c_SM,iGC1,GC2,GC3,GC4,GC5,iGC5,GC6,GC7,GC8,GC11,GC12,GC13,GC14,iGR34,mso, &
1018 GC15,GR1,GR3,GR13,GR14,GR15,GR17,GR31,iGR31,GR32,GR33,GR34,GR35,GR36,GI4, &
1019 GI6,GI20,GI21,GI22,GI31,GI32,GI33,GI34,GI35,iGI31,GI11,GI36,GI37,GI40,iGG34, &
1020 GS09,GS11,GS12,GS13,iGS20,GS31,iGS31,GS32,GS33,GS34,GS35,GS36,GS40,iGS40, &
1021 GS50,GG09,GG11,GG12,GG13,GG31,iGG31,GG32,GG33,GG34,GG35,GG36,GG40,iGG99,GH09,&
1022 GH11,GH12,GH13,GH31,GH32,GH33,GH40,GR50,GG50,iGH34,GH50,iGH99,iGH31,iGS34, &
1023 iGS20_D3,GS40_D3,cms_D3,eds,fds,rfact_FvFm
1025 !Size distribution parameters:
1026 real, parameter :: MUc = 3. !shape parameter for cloud
1027 real, parameter :: alpha_c = 1. !shape parameter for cloud
1028 real, parameter :: alpha_r = 0. !shape parameter for rain
1029 real, parameter :: alpha_i = 0. !shape parameter for ice
1030 real, parameter :: alpha_s = 0. !shape parameter for snow
1031 real, parameter :: alpha_g = 0. !shape parameter for graupel
1032 real, parameter :: alpha_h = 0. !shape parameter for hail
1033 real, parameter :: No_s_max = 1.e+8 !max. allowable intercept for snow [m-4]
1034 real, parameter :: lamdas_min= 500. !min. allowable LAMDA_s [m-1]
1036 !For single-moment:
1037 real, parameter :: No_r_SM = 1.e+7 !intercept parameter for rain [m-4]
1038 real, parameter :: No_g_SM = 4.e+6 !intercept parameter for graupel [m-4]
1039 real, parameter :: No_h_SM = 1.e+5 !intercept parameter for hail [m-4]
1040 !note: No_s = f(T), rather than a fixed value
1041 !------------------------------------!
1042 ! Symbol convention: (dist. params.) ! MY05: Milbrandt & Yau, 2005a,b (JAS)
1043 ! MY05 F94 CP00 ! F94: Ferrier, 1994 (JAS)
1044 ! ------ -------- ------ ! CP00: Cohard & Pinty, 2000a,b (QJGR)
1045 ! ALFx ALPHAx MUx-1 !
1046 ! MUx (1) ALPHAx !
1047 ! ALFx+1 ALPHAx+1 MUx !
1048 !------------------------------------!
1049 ! Note: The symbols for MU and ALPHA are REVERSED from that of CP2000a,b
1050 ! Explicit appearance of MUr = 1. has been removed.
1052 ! Fallspeed parameters:
1053 real, parameter :: afr= 149.100, bfr= 0.5000 !Tripoloi and Cotton (1980)
1054 real, parameter :: afi= 71.340, bfi= 0.6635 !Ferrier (1994)
1055 real, parameter :: afs= 11.720, bfs= 0.4100 !Locatelli and Hobbs (1974)
1056 real, parameter :: afg= 19.300, bfg= 0.3700 !Ferrier (1994)
1057 real, parameter :: afh= 206.890, bfh= 0.6384 !Ferrier (1994)
1058 !options:
1059 !real, parameter :: afs= 8.996, bfs= 0.4200 !Ferrier (1994)
1060 !real, parameter :: afg= 6.4800, bfg= 0.2400 !LH74 (grpl-like snow of lump type)
1062 real, parameter :: epsQ = 1.e-14 !kg kg-1, min. allowable mixing ratio
1063 real, parameter :: epsN = 1.e-3 !m-3, min. allowable number concentration
1064 real, parameter :: epsQ2 = 1.e-6 !kg kg-1, mixing ratio threshold for diagnostics
1065 real, parameter :: epsVIS= 1. !m, min. allowable visibility
1067 real, parameter :: iLAMmin1= 1.e-6 !min. iLAMx (prevents underflow in Nox and VENTx calcs)
1068 real, parameter :: iLAMmin2= 1.e-10 !min. iLAMx (prevents underflow in Nox and VENTx calcs)
1069 real, parameter :: eps = 1.e-32
1071 real, parameter :: deg = 400., mgo= 1.6e-10
1072 real, parameter :: deh = 900.
1073 real, parameter :: dei = 500., mio=1.e-12, Nti0=1.e3
1074 real, parameter :: dew = 1000.
1075 real, parameter :: desFix= 100. !used for snowSpherical = .true.
1076 real, parameter :: desMax= 500.
1077 real, parameter :: Dso = 125.e-6 ![m]; embryo snow diameter (mean-volume particle)
1078 real, parameter :: dmr = 3., dmi= 3., dmg= 3., dmh= 3.
1080 ! NOTE: VxMAX below are the max.allowable mass-weighted fallspeeds for sedimentation.
1081 ! Thus, Vx corresponds to DxMAX (at sea-level) times the max. density factor, GAM.
1082 ! [GAMmax=sqrt(DEo/DEmin)=sqrt(1.25/0.4)~2.] e.g. VrMAX = 2.*8.m/s = 16.m/s
1083 real, parameter :: DrMax= 5.e-3, VrMax= 16., epsQr_sedi= 1.e-8
1084 real, parameter :: DiMax= 5.e-3, ViMax= 2., epsQi_sedi= 1.e-10
1085 real, parameter :: DsMax= 5.e-3, VsMax= 4., epsQs_sedi= 1.e-8
1086 real, parameter :: DgMax= 2.e-3, VgMax= 6., epsQg_sedi= 1.e-8
1087 real, parameter :: DhMax= 80.e-3, VhMax= 25., epsQh_sedi= 1.e-10
1089 real, parameter :: CPW = 4218. ![J kg-1 K-1] specific heat capacity of water
1090 real, parameter :: DEo = 1.225 ![kg m-3] reference air density
1091 real, parameter :: thrd = 1./3.
1092 real, parameter :: sixth = 0.5*thrd
1093 real, parameter :: Ers = 1., Eci= 1. !collection efficiencies, Exy, between categories x and y
1094 real, parameter :: Eri = 1., Erh= 1.
1095 real, parameter :: Xdisp = 0.25 !dispersion of the fall velocity of ice
1096 real, parameter :: aa11 = 9.44e15, aa22= 5.78e3, Rh= 41.e-6
1097 real, parameter :: Avx = 0.78, Bvx= 0.30 !ventilation coefficients [F94 (B.36)]
1098 real, parameter :: Abigg = 0.66, Bbigg= 100. !parameters in probabilistic freezing
1099 real, parameter :: fdielec = 4.464 !ratio of dielectric factor, |K|w**2/|K|i**2
1100 real, parameter :: zfact = 1.e+18 !conversion factor for m-3 to mm2 m-6 for Ze
1101 real, parameter :: minZET = -99. ![dBZ] min threshold for ZET
1102 real, parameter :: maxVIS = 99.e+3 ![m] max. allowable VIS (visibility)
1103 real, parameter :: Drshed = 0.001 ![m] mean diam. of drop shed during wet growth
1104 real, parameter :: SIGcTHRS = 15.e-6 !threshold cld std.dev. before autoconversion
1105 real, parameter :: KK1 = 3.03e3 !parameter in Long (1974) kernel
1106 real, parameter :: KK2 = 2.59e15 !parameter in Long (1974) kernel
1107 real, parameter :: Dhh = 82.e-6 ![m] diameter that rain hump first appears
1108 real, parameter :: zMax_sedi = 20000. ![m] maximum height to compute sedimentation
1109 real, parameter :: Dr_large = 200.e-6 ![m] size threshold to distinguish rain/drizzle for precip rates
1110 real, parameter :: Ds_large = 200.e-6 ![m] size threshold to distinguish snow/snow-grains for precip rates
1111 real, parameter :: Dh_large = 1.0e-2 ![m] size threshold for "large" hail precipitation rate
1112 real, parameter :: Dh_min = 1.0e-3 ![m] size threhsold for below which hail converts to graupel
1113 real, parameter :: Dr_3cmpThrs = 2.5e-3 ![m] size threshold for hail production from 3-comp freezing
1114 real, parameter :: w_CNgh = 3. ![m s-1] vertical motion threshold for CNgh
1115 real, parameter :: Ngh_crit = 1.e+0 ![m-3] critical graupel concentration for CNgh
1116 real, parameter :: Tc_FZrh = -10. !temp-threshold (C) for FZrh
1117 real, parameter :: CNsgThres = 1.0 !threshold for CLcs/VDvs ratio for CNsg
1118 real, parameter :: capFact_i = 0.5 !capacitace factor for ice (C= 0.5*D*capFact_i)
1119 real, parameter :: capFact_s = 0.5 !capacitace factor for snow (C= 0.5*D*capFact_s)
1120 real, parameter :: Fv_Dsmin = 125.e-6 ![m] min snow size to compute volume flux
1121 real, parameter :: Fv_Dsmax = 0.008 ![m] max snow size to compute volume flux
1122 real, parameter :: Ni_max = 1.e+7 ![m-3] max ice crystal concentration
1123 real, parameter :: satw_peak = 1.01 !assumed max. peak saturation w.r.t. water (for calc of Simax)
1125 !------------------------------------------------------------------------------!
1126 !-- For GEM:
1127 !#include "consphy.cdk"
1129 !-- For WRF + kin_1d + kin_2d:
1130 real, parameter :: CPI =.21153e+4 !J kg-1 K-1; specific heat capacity of ice
1131 real, parameter :: TRPL =.27316e+3 !K; triple point of water
1132 real, parameter :: TCDK =.27315e+3 !conversion from kelvin to celsius
1133 real, parameter :: RAUW =.1e+4 !density of liquid H2O
1134 real, parameter :: EPS2 =.3780199778986 !1 - EPS1
1135 real, parameter :: TGL =.27316e+3 !K; ice temperature in the atmosphere
1136 real, parameter :: CONSOL =.1367e+4 !W m-2; solar constant
1137 real, parameter :: RAYT =.637122e+7 !M; mean radius of the earth
1138 real, parameter :: STEFAN =.566948e-7 !J m-2 s-1 K-4; Stefan-Boltzmann constant
1139 real, parameter :: PI =.314159265359e+1 !PI constant = ACOS(-1)
1140 real, parameter :: OMEGA =.7292e-4 !s-1; angular speed of rotation of the earth
1141 real, parameter :: KNAMS =.514791 !conversion from knots to m/s
1142 real, parameter :: STLO =.6628486583943e-3 !K s2 m-2; Schuman-Newell Lapse Rate
1143 real, parameter :: KARMAN =.35 !Von Karman constant
1144 real, parameter :: RIC =.2 !Critical Richardson number
1145 !-- For kin_1d + kin_2d (exclude from WRF):
1146 real, parameter :: CHLC =.2501e+7 !J kg-1; latent heat of condensation
1147 real, parameter :: CHLF =.334e+6 !J kg-1; latent heat of fusion
1148 real, parameter :: CPD =.100546e+4 !J K-1 kg-1; specific heat of dry air
1149 real, parameter :: CPV =.186946e+4 !J K-1 kg-1; specific heat of water vapour
1150 real, parameter :: RGASD =.28705e+3 !J K-1 kg-1; gas constant for dry air
1151 real, parameter :: RGASV =.46151e+3 !J K-1 kg-1; gas constant for water vapour
1152 real, parameter :: EPS1 =.62194800221014 !RGASD/RGASV
1153 real, parameter :: DELTA =.6077686814144 !1/EPS1 - 1
1154 real, parameter :: CAPPA =.28549121795 !RGASD/CPD
1155 real, parameter :: GRAV =.980616e+1 !M s-2; gravitational acceleration
1156 !==
1157 !------------------------------------------------------------------------------!
1159 ! Constants used for contact ice nucleation:
1160 real, parameter :: LAMa0 = 6.6e-8 ![m] mean free path at T0 and p0 [W95_eqn58]
1161 real, parameter :: T0 = 293.15 ![K] ref. temp.
1162 real, parameter :: p0 = 101325. ![Pa] ref. pres.
1163 real, parameter :: Ra = 1.e-6 ![m] aerosol (IN) radius [M92 p.713; W95_eqn60]
1164 real, parameter :: kBoltz = 1.381e-23 !Boltzmann's constant
1165 real, parameter :: KAPa = 5.39e5 !aerosol thermal conductivity
1167 !Test switches:
1168 logical, parameter :: DEBUG_ON = .false. !.true. to switch on debugging checks/traps throughout code
1169 logical, parameter :: DEBUG_abort = .true. !.true. will result in forced abort in s/r 'check_values'
1170 logical, parameter :: iceDep_ON = .true. !.false. to suppress depositional growth of ice
1171 logical, parameter :: grpl_ON = .true. !.false. to suppress graupel initiation
1172 logical, parameter :: hail_ON = .true. !.false. to suppress hail initiation
1173 logical, parameter :: rainAccr_ON = .true. ! rain accretion and self-collection ON/OFF
1174 logical, parameter :: snowSpherical = .false. !.true.: m(D)=(pi/6)*const_des*D^3 | .false.: m(D)= 0.069*D^2
1175 integer, parameter :: primIceNucl = 1 !1= Meyers+contact ; 2= Cooper
1176 real, parameter :: outfreq = 60. !frequency to compute output diagnostics [s]
1178 real, dimension(size(QC,dim=1),size(QC,dim=2)) :: DE,iDE,iDP,QSW,QSI,DZ,iDZ,zz,VQQ, &
1179 gamfact,pres,zheight,QC_in,QR_in,NC_in,NR_in
1180 real, dimension(size(QC,dim=1)) :: fluxM_r,fluxM_i,fluxM_s,fluxM_g, &
1181 fluxM_h,dum
1182 integer, dimension(size(QC,dim=1)) :: activeColumn
1184 !-- for use with sedimentation on subset of levels:
1185 !integer :: k_sub,nk_sub,nk_skip
1186 !integer :: status !for allocate/deallocate statements (0 for success)
1187 !integer, allocatable, dimension(:) :: kfull,kskip
1188 !integer, allocatable, dimension(:,:) :: iint
1189 !real, dimension(:,:), allocatable :: DE_sub,iDE_sub,iDP_sub,pres_sub, &
1190 !== DZ_sub,zheight_sub,iDZ_sub,gamfact_sub
1193 !==================================================================================!
1195 !----------------------------------------------------------------------------------!
1196 ! PART 1: Prelimiary Calculations !
1197 !----------------------------------------------------------------------------------!
1199 !Switch on here later, once it is certain that calling routine is not supposed to
1200 !pass negative values of tracers.
1201 if (DEBUG_ON) call check_values(Q,T,QC,QR,QI,QN,QG,QH,NC,NR,NY,NN,NG,NH,epsQ,epsN,.false.,DEBUG_abort,100)
1203 if (nk_BOTTOM) then
1204 ! !GEM / kin_1d:
1205 ktop = 1 !k of top level
1206 kbot = nk !k of bottom level
1207 kdir = -1 !direction of vertical leveling (k: 1=top, nk=bottom)
1208 else
1209 !WRF / kin_2d: (assuming no array flipping in wrapper)
1210 ktop = nk !k of top level
1211 kbot = 1 !k of bottom level
1212 kdir = 1 !direction of vertical leveling (k: 1=bottom, nk=top)
1213 endif
1215 do k= kbot,ktop,kdir
1216 pres(:,k)= PS(:)*sigma(:,k) !air pressure [Pa]
1217 do i=1,ni
1218 QSW(i,k) = qsat(T(i,k),pres(i,k),0) !wrt. liquid water
1219 QSI(i,k) = qsat(T(i,k),pres(i,k),1) !wrt. ice
1220 enddo
1221 enddo
1223 !Air density:
1224 DE = pres/(RGASD*T)
1225 iDE = 1./DE
1227 !Convert N from #/kg to #/m3:
1228 NC = NC*DE
1229 NR = NR*DE
1230 NY = NY*DE
1231 NN = NN*DE
1232 NG = NG*DE
1233 NH = NH*DE
1235 ! The SS(i,k,n) array is passed to 'vkuocon6' where it is converted into individual
1236 ! arrays [a_ss01(i,k)] and then passed to the volatile bus for output as 3-D diagnostic
1237 ! output variables, for testing purposes. For example, to output the
1238 ! instantanous value of the deposition rate, add 'SS(i,k,1) = QVDvi' in the
1239 ! appropriate place. It can then be output as a 3-D physics variable by adding
1240 ! SS01 to the sortie_p list in 'outcfgs.out'
1241 SS= 0.
1243 !Compute diagnostic values only every 'outfreq' minutes:
1244 !calcDiag= (mod(DT*float(KOUNT),outfreq)==0.)
1245 calcDiag = .true. !compute diagnostics every step (for time-series output)
1247 !#### These need only to be computed once per model integration:
1248 ! (note: These variables must be declared with the SAVE attribute)
1250 ! if (KOUNT==0) then
1251 !*** For restarts, these values are not saved. Therefore, the condition statement
1252 ! must be modified to something like: IF (MOD(Step_rsti,KOUNT).eq.0) THEN
1253 ! in order that these be computed only on the first step of a given restart.
1254 ! (...to be done. For now, changing condition to IF(TRUE) to compute at each step.)
1257 if (.TRUE.) then
1259 PI2 = PI*2.
1260 PIov4 = 0.25*PI
1261 PIov6 = PI*sixth
1262 CHLS = CHLC+CHLF !J k-1; latent heat of sublimation
1263 LCP = CHLC/CPD
1264 LFP = CHLF/CPD
1265 iCHLF = 1./CHLF
1266 LSP = LCP+LFP
1267 ck5 = 4098.170*LCP
1268 ck6 = 5806.485*LSP
1269 idt = 1./dt
1270 imgo = 1./mgo
1271 idew = 1./dew
1272 idei = 1./dei
1273 ideg = 1./deg
1274 ideh = 1./deh
1276 !Constants based on size distribution parameters:
1278 ! Mass parameters [ m(D) = cD^d ]
1279 cmr = PIov6*dew; icmr= 1./cmr
1280 cmi = 440.; icmi= 1./cmi
1281 cmg = PIov6*deg; icmg= 1./cmg
1282 cmh = PIov6*deh; icmh= 1./cmh
1284 cms_D3 = PIov6*desFix !used for snowSpherical = .T. or .F.
1285 if (snowSpherical) then
1286 cms = cms_D3
1287 dms = 3.
1288 else
1289 ! cms = 0.0690; dms = 2.000 !Cox, 1988 (QJRMS)
1290 cms = 0.1597; dms = 2.078 !Brandes et al., 2007 (JAMC)
1291 endif
1292 icms = 1./cms
1293 idms = 1./dms
1294 mso = cms*Dso**dms
1295 imso = 1./mso
1296 !bulk density parameters: [rho(D) = eds*D^fds]
1297 ! These are implied by the mass-diameter parameters, by computing the bulk
1298 ! density of a sphere with the equaivalent mass.
1299 ! e.g. m(D) = cD^d = (pi/6)rhoD^3 and solve for rho(D)
1300 eds = cms/PIov6
1301 fds = dms-3.
1302 if (fds/=-1. .and..not.snowSpherical) GS50= gamma(1.+fds+alpha_s)
1304 ! Cloud:
1305 iMUc = 1./MUc
1306 GC1 = gamma(alpha_c+1.0)
1307 iGC1 = 1./GC1
1308 GC2 = gamma(alpha_c+1.+3.0*iMUc) !i.e. gamma(alf + 4)
1309 GC3 = gamma(alpha_c+1.+6.0*iMUc) !i.e. gamma(alf + 7)
1310 GC4 = gamma(alpha_c+1.+9.0*iMUc) !i.e. gamma(alf + 10)
1311 GC11 = gamma(1.0*iMUc+1.0+alpha_c)
1312 GC12 = gamma(2.0*iMUc+1.0+alpha_c)
1313 GC5 = gamma(1.0+alpha_c)
1314 iGC5 = 1./GC5
1315 GC6 = gamma(1.0+alpha_c+1.0*iMUc)
1316 GC7 = gamma(1.0+alpha_c+2.0*iMUc)
1317 GC8 = gamma(1.0+alpha_c+3.0*iMUc)
1318 GC13 = gamma(3.0*iMUc+1.0+alpha_c)
1319 GC14 = gamma(4.0*iMUc+1.0+alpha_c)
1320 GC15 = gamma(5.0*iMUc+1.0+alpha_c)
1321 icexc9 = 1./(GC2*iGC1*PIov6*dew)
1322 !specify cloud droplet number concentration [m-3] based on 'CCNtype' (1-moment):
1323 if (CCNtype==1) then
1324 N_c_SM = 0.8e+8 !maritime
1325 elseif (CCNtype==2) then
1326 N_c_SM = 2.0e+8 !continental 1
1327 elseif (CCNtype==3) then
1328 N_c_SM = 5.0e+8 !continental 2 (polluted)
1329 else
1330 N_c_SM = 2.0e+8 !default (cont1), if 'CCNtype' specified incorrectly
1331 endif
1333 ! Rain:
1334 cexr1 = 1.+alpha_r+dmr+bfr
1335 cexr2 = 1.+alpha_r+dmr
1336 GR17 = gamma(2.5+alpha_r+0.5*bfr)
1337 GR31 = gamma(1.+alpha_r)
1338 iGR31 = 1./GR31
1339 GR32 = gamma(2.+alpha_r)
1340 GR33 = gamma(3.+alpha_r)
1341 GR34 = gamma(4.+alpha_r)
1342 iGR34 = 1./GR34
1343 GR35 = gamma(5.+alpha_r)
1344 GR36 = gamma(6.+alpha_r)
1345 GR37 = gamma(7.+alpha_r)
1346 GR50 = (No_r_SM*GR31)**0.75 !for 1-moment or Nr-initialization
1347 cexr5 = 2.+alpha_r
1348 cexr6 = 2.5+alpha_r+0.5*bfr
1349 cexr9 = cmr*GR34*iGR31; icexr9= 1./cexr9
1350 cexr3 = 1.+bfr+alpha_r
1351 cexr4 = 1.+alpha_r
1352 ckQr1 = afr*gamma(1.+alpha_r+dmr+bfr)/gamma(1.+alpha_r+dmr)
1353 ckQr2 = afr*gamma(1.+alpha_r+bfr)*GR31
1354 ckQr3 = afr*gamma(7.+alpha_r+bfr)/GR37
1356 ! Ice:
1357 GI4 = gamma(alpha_i+dmi+bfi)
1358 GI6 = gamma(2.5+bfi*0.5+alpha_i)
1359 GI11 = gamma(1.+bfi+alpha_i)
1360 GI20 = gamma(0.+bfi+1.+alpha_i)
1361 GI21 = gamma(1.+bfi+1.+alpha_i)
1362 GI22 = gamma(2.+bfi+1.+alpha_i)
1363 GI31 = gamma(1.+alpha_i)
1364 iGI31 = 1./GI31
1365 GI32 = gamma(2.+alpha_i)
1366 GI33 = gamma(3.+alpha_i)
1367 GI34 = gamma(4.+alpha_i)
1368 GI35 = gamma(5.+alpha_i)
1369 GI36 = gamma(6.+alpha_i)
1370 GI40 = gamma(1.+alpha_i+dmi)
1371 icexi9 = 1./(cmi*gamma(1.+alpha_i+dmi)*iGI31)
1372 ckQi1 = afi*gamma(1.+alpha_i+dmi+bfi)/GI40
1373 ckQi2 = afi*GI11*iGI31
1374 ckQi4 = 1./(cmi*GI40*iGI31)
1376 ! Snow:
1377 cexs1 = 2.5+0.5*bfs+alpha_s
1378 cexs2 = 1.+alpha_s+dms
1379 icexs2 = 1./cexs2
1380 GS09 = gamma(2.5+bfs*0.5+alpha_s)
1381 GS11 = gamma(1.+bfs+alpha_s)
1382 GS12 = gamma(2.+bfs+alpha_s)
1383 GS13 = gamma(3.+bfs+alpha_s)
1384 GS31 = gamma(1.+alpha_s)
1385 iGS31 = 1./GS31
1386 GS32 = gamma(2.+alpha_s)
1387 GS33 = gamma(3.+alpha_s)
1388 GS34 = gamma(4.+alpha_s)
1389 iGS34 = 1./GS34
1390 GS35 = gamma(5.+alpha_s)
1391 GS36 = gamma(6.+alpha_s)
1392 GS40 = gamma(1.+alpha_s+dms)
1393 iGS40 = 1./GS40
1394 iGS20 = 1./(GS40*iGS31*cms)
1395 ckQs1 = afs*gamma(1.+alpha_s+dms+bfs)*iGS40
1396 ckQs2 = afs*GS11*iGS31
1397 GS40_D3 = gamma(1.+alpha_s+3.)
1398 iGS20_D3= 1./(GS40_D3*iGS31*cms_D3)
1399 rfact_FvFm= PIov6*icms*gamma(4.+bfs+alpha_s)/gamma(1.+dms+bfs+alpha_s)
1401 ! Graupel:
1402 GG09 = gamma(2.5+0.5*bfg+alpha_g)
1403 GG11 = gamma(1.+bfg+alpha_g)
1404 GG12 = gamma(2.+bfg+alpha_g)
1405 GG13 = gamma(3.+bfg+alpha_g)
1406 GG31 = gamma(1.+alpha_g)
1407 iGG31 = 1./GG31
1408 GG32 = gamma(2.+alpha_g)
1409 GG33 = gamma(3.+alpha_g)
1410 GG34 = gamma(4.+alpha_g)
1411 iGG34 = 1./GG34
1412 GG35 = gamma(5.+alpha_g)
1413 GG36 = gamma(6.+alpha_g)
1414 GG40 = gamma(1.+alpha_g+dmg)
1415 iGG99 = 1./(GG40*iGG31*cmg)
1416 GG50 = (No_g_SM*GG31)**0.75 !for 1-moment only
1417 ckQg1 = afg*gamma(1.+alpha_g+dmg+bfg)/GG40
1418 ckQg2 = afg*GG11*iGG31
1419 ckQg4 = 1./(cmg*GG40*iGG31)
1421 ! Hail:
1422 GH09 = gamma(2.5+bfh*0.5+alpha_h)
1423 GH11 = gamma(1.+bfh+alpha_h)
1424 GH12 = gamma(2.+bfh+alpha_h)
1425 GH13 = gamma(3.+bfh+alpha_h)
1426 GH31 = gamma(1.+alpha_h)
1427 iGH31 = 1./GH31
1428 GH32 = gamma(2.+alpha_h)
1429 GH33 = gamma(3.+alpha_h)
1430 iGH34 = 1./gamma(4.+alpha_h)
1431 GH40 = gamma(1.+alpha_h+dmh)
1432 iGH99 = 1./(GH40*iGH31*cmh)
1433 GH50 = (No_h_SM*GH31)**0.75 !for 1-moment only
1434 ckQh1 = afh*gamma(1.+alpha_h+dmh+bfh)/GH40
1435 ckQh2 = afh*GH11*iGH31
1436 ckQh4 = 1./(cmh*GH40*iGH31)
1438 endif !if (KOUNT=0)
1439 !####
1441 !=======================================================================================!
1442 ! if (DEBUG_ON) call check_values(Q,T,QC,QR,QI,QN,QG,QH,NC,NR,NY,NN,NG,NH,epsQ,epsN,.false.,DEBUG_abort,200)
1444 !--- Ensure consistency between moments:
1445 do k= kbot,ktop,kdir
1446 do i= 1,ni
1448 tmp1 = QSW(i,k)/max(Q(i,k),1.e-20) !saturation w.r.t. water
1449 tmp2 = QSI(i,k)/max(Q(i,k),1.e-20) !saturation w.r.t. ice
1450 !NOTE: Hydrometeor (Qx,Nx) is clipped if Qx is tiny or if Qx is small
1451 ! and is expected to completely evaporate/sublimate in one time step
1452 ! anyway. (To avoid creating mass from a parallel universe, only
1453 ! positive Qx values are added to water vapor upon clipping.)
1454 ! ** RH thresholds for clipping need to be tuned (especially for graupel
1455 ! and hail, for which it may be preferable to reduce threshold)
1458 !cloud:
1459 if (QC(i,k)>epsQ .and. NC(i,k)<epsN) then
1460 NC(i,k) = N_c_SM
1461 elseif (QC(i,k)<=epsQ .or. (QC(i,k)<epsQ2 .and. tmp1<0.90)) then
1462 tmp3 = max(0., QC(i,k))
1463 Q(i,k) = Q(i,k) + tmp3
1464 T(i,k) = T(i,k) - LCP*tmp3
1465 QC(i,k) = 0.
1466 NC(i,k) = 0.
1467 endif
1469 !rain
1470 if (QR(i,k)>epsQ .and. NR(i,k)<epsN) then
1471 NR(i,k) = (No_r_SM*GR31)**(3./(4.+alpha_r))*(GR31*iGR34*DE(i,k)*QR(i,k)* &
1472 icmr)**((1.+alpha_r)/(4.+alpha_r))
1473 elseif (QR(i,k)<=epsQ .or. (QR(i,k)<epsQ2 .and. tmp1<0.90)) then
1474 tmp3 = max(0., QR(i,k))
1475 Q(i,k) = Q(i,k) + tmp3
1476 T(i,k) = T(i,k) - LCP*tmp3
1477 QR(i,k) = 0.
1478 NR(i,k) = 0.
1479 endif
1481 !ice:
1482 if (QI(i,k)>epsQ .and. NY(i,k)<epsN) then
1483 ! NY(i,k) = N_Cooper(TRPL,T(i,k))
1484 NY(i,k) = max(2.*epsN, N_Cooper(TRPL,T(i,k)) )
1485 elseif (QI(i,k)<=epsQ .or. (QI(i,k)<epsQ2 .and. tmp2<0.80)) then
1486 tmp3 = max(0., QI(i,k))
1487 Q(i,k) = Q(i,k) + tmp3
1488 T(i,k) = T(i,k) - LSP*tmp3
1489 QI(i,k) = 0.
1490 NY(i,k) = 0.
1491 endif
1493 !snow:
1494 if (QN(i,k)>epsQ .and. NN(i,k)<epsN) then
1495 No_s = Nos_Thompson(TRPL,T(i,k))
1496 NN(i,k) = (No_s*GS31)**(dms*icexs2)*(GS31*iGS40*icms*DE(i,k)*QN(i,k))** &
1497 ((1.+alpha_s)*icexs2)
1498 elseif (QN(i,k)<=epsQ .or. (QN(i,k)<epsQ2 .and. tmp2<0.80)) then
1499 tmp3 = max(0., QN(i,k))
1500 Q(i,k) = Q(i,k) + tmp3
1501 T(i,k) = T(i,k) - LSP*tmp3
1502 QN(i,k) = 0.
1503 NN(i,k) = 0.
1504 endif
1506 !grpl:
1507 if (QG(i,k)>epsQ .and. NG(i,k)<epsN) then
1508 NG(i,k) = (No_g_SM*GG31)**(3./(4.+alpha_g))*(GG31*iGG34*DE(i,k)*QG(i,k)* &
1509 icmg)**((1.+alpha_g)/(4.+alpha_g))
1510 elseif (QG(i,k)<=epsQ .or. (QG(i,k)<epsQ2 .and. tmp2<0.80)) then
1511 tmp3 = max(0., QG(i,k))
1512 Q(i,k) = Q(i,k) + tmp3
1513 T(i,k) = T(i,k) - LSP*tmp3
1514 QG(i,k) = 0.
1515 NG(i,k) = 0.
1516 endif
1518 !hail:
1519 if (QH(i,k)>epsQ .and. NH(i,k)<epsN) then
1520 NH(i,k) = (No_h_SM*GH31)**(3./(4.+alpha_h))*(GH31*iGH34*DE(i,k)*QH(i,k)* &
1521 icmh)**((1.+alpha_h)/(4.+alpha_h))
1522 elseif (QH(i,k)<=epsQ .or. (QH(i,k)<epsQ2 .and. tmp2<0.80)) then
1523 tmp3 = max(0., QH(i,k))
1524 Q(i,k) = Q(i,k) + tmp3
1525 T(i,k) = T(i,k) - LSP*tmp3
1526 QH(i,k) = 0.
1527 NH(i,k) = 0.
1528 endif
1530 enddo !i-loop
1531 enddo !k-loop
1532 !===
1534 if (DEBUG_ON) call check_values(Q,T,QC,QR,QI,QN,QG,QH,NC,NR,NY,NN,NG,NH,epsQ,epsN,.true.,DEBUG_abort,300)
1536 !Temporarily store arrays at time (t*) in order to compute warm-rain coalescence
1537 ! equations (Part 3a):
1538 QC_in = QC
1539 QR_in = QR
1540 NC_in = NC
1541 NR_in = NR
1544 !Air-density factor (for fall velocity computations):
1545 do i= 1,ni
1546 gamfact(i,:) = sqrt(DEo/(DE(i,:)))
1547 enddo
1549 !Pressure difference between levels (used for sedimentation)
1550 iDP(:,kbot) = 1./(PS(:)-pres(:,kbot))
1551 do k = kbot+kdir,ktop,kdir
1552 iDP(:,k) = 1./(pres(:,k-kdir)-pres(:,k))
1553 enddo
1555 !Compute thickness of layers for sedimentation calculation: (optional use if height coordiates)
1556 ! note - from 'cldoptx4.ftn': dz(i,k)= dp(i,k)/aird(i,k)*rec_grav
1557 iDZ = DE*GRAV*iDP
1558 DZ = 1./iDZ
1560 !Compute height above surface (zheight) and height above lowest prog level (zz):
1561 zheight(:,kbot)= DZ(:,kbot)
1562 zz(:,kbot)= 0. !<-- define height of lowest prognosic level to be 0. (for diagnostics only)
1563 do k = kbot+kdir,ktop,kdir
1564 zheight(:,k) = zheight(:,k-kdir) + DZ(:,k)
1565 zz(:,k) = zz(:,k-kdir) + DZ(:,k)
1566 enddo
1568 !Determine the upper-most level in each column to which to compute sedimentation:
1569 !ktop_sedi= ktop+kdir !<Optional (in place of code below) - to compute sedimentation at all levels
1570 ktop_sedi= 0
1571 do i=1,ni
1572 do k= ktop,kbot,-kdir
1573 ktop_sedi(i)= k
1574 if (zheight(i,k)<zMax_sedi) exit
1575 enddo
1576 enddo
1578 !----------------------------------------------------------------------------------!
1580 !----------------------------------------------------------------------------------!
1581 ! End of Preliminary Calculation section (Part 1) !
1582 !----------------------------------------------------------------------------------!
1584 if (DEBUG_ON) call check_values(Q,T,QC,QR,QI,QN,QG,QH,NC,NR,NY,NN,NG,NH,epsQ,epsN,.true.,DEBUG_abort,400)
1586 !----------------------------------------------------------------------------------!
1587 ! PART 2: Cold Microphysics Processes !
1588 !----------------------------------------------------------------------------------!
1590 ! Determine the active grid points (i.e. those which scheme should treat):
1591 activePoint = .false.
1592 DO k= ktop-kdir,kbot,-kdir
1593 DO i=1,ni
1594 log1= ((QI(i,k)+QG(i,k)+QN(i,k)+QH(i,k))<epsQ) !no solid (i,g,s,h)
1595 log2= ((QC(i,k)+QR(i,k)) <epsQ) !no liquid (c,r)
1596 log3= ((T(i,k)>TRPL) .and. log1) !T>0C & no i,g,s,h
1597 log4= log1.and.log2.and.(Q(i,k)<QSI(i,k)) !no sol. or liq.; subsat(i)
1598 if (.not.( log3 .or. log4 ) .and. icephase_ON) then
1599 activePoint(i,k)= .true.
1600 endif
1601 ENDDO
1602 ENDDO
1604 ! Size distribution parameters:
1605 ! Note: + 'thrd' should actually be '1/dmx'(but dmx=3 for all categories x)
1606 ! + If Qx=0, LAMx etc. are never be used in any calculations
1607 ! (If Qc=0, CLcy etc. will never be calculated. iLAMx is set to 0
1608 ! to avoid possible problems due to bugs.)
1610 DO k= ktop-kdir,kbot,-kdir
1611 DO i= 1,ni
1612 IF (activePoint(i,k)) THEN
1614 Tc= T(i,k)-TRPL
1615 if (Tc<-120. .or. Tc>50.) then
1616 print*, '***WARNING*** -- In MICROPHYSICS -- Ambient Temp.(C),step,i,k:',Tc,kount,i,k
1617 !stop
1618 endif
1619 Cdiff = (2.2157e-5+0.0155e-5*Tc)*1.e5/pres(i,k)
1620 MUdyn = 1.72e-5*(393./(T(i,k)+120.))*(T(i,k)/TRPL)**1.5 !RYp.102
1621 MUkin = MUdyn*iDE(i,k)
1622 iMUkin= 1./MUkin
1623 ScTHRD= (MUkin/Cdiff)**thrd ! i.e. Sc^(1/3)
1624 Ka = 2.3971e-2 + 0.0078e-2*Tc !therm.cond.(air)
1625 Kdiff = (9.1018e-11*T(i,k)*T(i,k)+8.8197e-8*T(i,k)-(1.0654e-5)) !therm.diff.(air)
1626 gam = gamfact(i,k)
1628 !Collection efficiencies:
1629 Eis = min(0.05*exp(0.1*Tc),1.) !Ferrier, 1995 (Table 1)
1630 Eig = min(0.01*exp(0.1*Tc),1.) !dry (Eig=1.0 for wet growth)
1631 Eii = 0.1*Eis
1632 Ess = Eis; Eih = Eig; Esh = Eig
1633 iEih = 1./Eih
1634 iEsh = 1./Esh
1635 !note: Eri=Ers=Erh=1. (constant parameters)
1636 ! - Ecs is computed in CLcs section
1637 ! - Ech is computed in CLch section
1638 ! - Ecg is computed in CLcg section
1639 ! - Erg is computed in CLrg section
1641 qvs0 = qsat(TRPL,pres(i,k),0) !sat.mix.ratio at 0C
1642 DELqvs = qvs0-(Q(i,k))
1644 ! Cloud:
1645 if (QC(i,k)>epsQ) then
1646 iQC = 1./QC(i,k)
1647 iNC = 1./NC(i,k)
1648 Dc = Dm_x(DE(i,k),QC(i,k),iNC,icmr,thrd)
1650 iLAMc = iLAMDA_x(DE(i,k),QC(i,k),iNC,icexc9,thrd)
1651 iLAMc2 = iLAMc *iLAMc
1652 iLAMc3 = iLAMc2*iLAMc
1653 iLAMc4 = iLAMc2*iLAMc2
1654 iLAMc5 = iLAMc3*iLAMc2
1655 else
1656 Dc = 0.; iLAMc3= 0.
1657 iLAMc = 0.; iLAMc4= 0.
1658 iLAMc2 = 0.; iLAMc5= 0.
1659 endif
1661 ! Rain:
1662 if (QR(i,k)>epsQ) then
1663 iQR = 1./QR(i,k)
1664 iNR = 1./NR(i,k)
1665 Dr = Dm_x(DE(i,k),QR(i,k),iNR,icmr,thrd)
1666 iLAMr = max( iLAMmin1, iLAMDA_x(DE(i,k),QR(i,k),iNR,icexr9,thrd) )
1667 tmp1 = 1./iLAMr
1668 iLAMr2 = iLAMr**2
1669 iLAMr3 = iLAMr**3
1670 iLAMr4 = iLAMr**4
1671 iLAMr5 = iLAMr**5
1672 if (Dr>40.e-6) then
1673 vr0 = gamfact(i,k)*ckQr1*iLAMr**bfr
1674 else
1675 vr0 = 0.
1676 endif
1677 else
1678 iLAMr = 0.; Dr = 0.; vr0 = 0.
1679 iLAMr2 = 0.; iLAMr3= 0.; iLAMr4= 0.; iLAMr5 = 0.
1680 endif
1682 ! Ice:
1683 if (QI(i,k)>epsQ) then
1684 iQI = 1./QI(i,k)
1685 iNY = 1./NY(i,k)
1686 iLAMi = max( iLAMmin2, iLAMDA_x(DE(i,k),QI(i,k),iNY,icexi9,thrd) )
1687 iLAMi2 = iLAMi**2
1688 iLAMi3 = iLAMi**3
1689 iLAMi4 = iLAMi**4
1690 iLAMi5 = iLAMi**5
1691 iLAMiB0= iLAMi**(bfi)
1692 iLAMiB1= iLAMi**(bfi+1.)
1693 iLAMiB2= iLAMi**(bfi+2.)
1694 vi0 = gamfact(i,k)*ckQi1*iLAMiB0
1695 Di = Dm_x(DE(i,k),QI(i,k),iNY,icmi,thrd)
1696 else
1697 iLAMi = 0.; vi0 = 0.; Di = 0.
1698 iLAMi2 = 0.; iLAMi3 = 0.; iLAMi4 = 0.; iLAMi5= 0.
1699 iLAMiB0= 0.; iLAMiB1= 0.; iLAMiB2= 0.
1700 endif
1702 ! Snow:
1703 if (QN(i,k)>epsQ) then
1704 iQN = 1./QN(i,k)
1705 iNN = 1./NN(i,k)
1706 iLAMs = max( iLAMmin2, iLAMDA_x(DE(i,k),QN(i,k),iNN,iGS20,idms) )
1707 iLAMs_D3= max(iLAMmin2, iLAMDA_x(DE(i,k),QN(i,k),iNN,iGS20_D3,thrd) )
1708 iLAMs2 = iLAMs**2
1709 iLAMsB0= iLAMs**(bfs)
1710 iLAMsB1= iLAMs**(bfs+1.)
1711 iLAMsB2= iLAMs**(bfs+2.)
1712 vs0 = gamfact(i,k)*ckQs1*iLAMsB0
1713 Ds = min(DsMax, Dm_x(DE(i,k),QN(i,k),iNN,icms,idms))
1714 if (snowSpherical) then
1715 des = desFix
1716 else
1717 des = des_OF_Ds(Ds,desMax,eds,fds)
1718 endif
1719 !!-- generalized equations (any alpha_s):
1720 ! No_s = (NN(i,k))*iGS31/iLAMs**(1.+alpha_s)
1721 ! VENTs = Avx*GS32*iLAMs**(2.+alpha_s)+Bvx*ScTHRD*sqrt(gam*afs*iMUkin)* &
1722 !!-- GS09*iLAMs**(2.5+0.5*bfs+alpha_s)
1723 !The following equations for No_s and VENTs is based on m(D)=(pi/6)*100.*D**3 for snow.
1724 ! Strict application of m(D)=c*D**2 would require re-derivation using implied
1725 ! definition of D as the MAXIMUM DIMENSION of an ellipsoid, rather than a sphere.
1726 ! For simplicity, the m-D^3 relation is applied -- used for VDvs and MLsr only.
1727 !No_s= NN(i,k)*iGS31/iLAMs !optimized for alpha_s=0
1728 No_s= NN(i,k)*iGS31/iLAMs_D3 !based on m-D^3 (consistent with VENTs, below)
1729 VENTs= Avx*GS32*iLAMs_D3**2. + Bvx*ScTHRD*sqrt(gamfact(i,k)*afs*iMUkin)*GS09* &
1730 iLAMs_D3**cexs1
1731 else
1732 iLAMs = 0.; vs0 = 0.; Ds = 0.; iLAMs2= 0.
1733 iLAMsB0= 0.; iLAMsB1= 0.; iLAMsB1= 0.
1734 des = desFix !used for 3-component freezing if QN=0 (even for snowSpherical=.F.)
1735 endif
1736 ides = 1./des
1739 ! Graupel:
1740 if (QG(i,k)>epsQ) then
1741 iQG = 1./QG(i,k)
1742 iNG = 1./NG(i,k)
1743 iLAMg = max( iLAMmin1, iLAMDA_x(DE(i,k),QG(i,k),iNG,iGG99,thrd) )
1744 iLAMg2 = iLAMg**2
1745 iLAMgB0= iLAMg**(bfg)
1746 iLAMgB1= iLAMg**(bfg+1.)
1747 iLAMgB2= iLAMg**(bfg+2.)
1748 !No_g = (NG(i,k))*iGG31/iLAMg**(1.+alpha_g)
1749 No_g= NG(i,k)*iGG31/iLAMg !optimized for alpha_g=0
1750 vg0 = gamfact(i,k)*ckQg1*iLAMgB0
1751 Dg = Dm_x(DE(i,k),QG(i,k),iNG,icmg,thrd)
1752 else
1753 iLAMg = 0.; vg0 = 0.; Dg = 0.; No_g = 0.
1754 iLAMg2 = 0.; iLAMgB0= 0.; iLAMgB1= 0.; iLAMgB1= 0.
1755 endif
1757 ! Hail:
1758 if (QH(i,k)>epsQ) then
1759 iQH = 1./QH(i,k)
1760 iNH = 1./NH(i,k)
1761 iLAMh = max( iLAMmin1, iLAMDA_x(DE(i,k),QH(i,k),iNH,iGH99,thrd) )
1762 iLAMh2 = iLAMh**2
1763 iLAMhB0= iLAMh**(bfh)
1764 iLAMhB1= iLAMh**(bfh+1.)
1765 iLAMhB2= iLAMh**(bfh+2.)
1766 No_h= NH(i,k)*iGH31/iLAMh**(1.+alpha_h)
1767 vh0 = gamfact(i,k)*ckQh1*iLAMhB0
1768 Dh = Dm_x(DE(i,k),QH(i,k),iNH,icmh,thrd)
1769 else
1770 iLAMh = 0.; vh0 = 0.; Dh = 0.; No_h= 0.
1771 iLAMhB0= 0.; iLAMhB1= 0.; iLAMhB1= 0.
1772 endif
1773 !------
1775 !Calculating ice-phase source/sink terms:
1777 ! Initialize all source terms to zero:
1778 QNUvi=0.; QVDvi=0.; QVDvs=0.; QVDvg=0.; QVDvh=0.
1779 QCLcs=0.; QCLcg=0.; QCLch=0.; QFZci=0.; QCLri=0.; QMLsr=0.
1780 QCLrs=0.; QCLrg=0.; QMLgr=0.; QCLrh=0.; QMLhr=0.; QFZrh=0.
1781 QMLir=0.; QCLsr=0.; QCLsh=0.; QCLgr=0.; QCNgh=0.
1782 QCNis=0.; QCLir=0.; QCLis=0.; QCLih=0.
1783 QIMsi=0.; QIMgi=0.; QCNsg=0.; QHwet=0.
1785 NCLcs= 0.; NCLcg=0.; NCLch=0.; NFZci=0.; NMLhr=0.; NhCNgh=0.
1786 NCLri= 0.; NCLrs=0.; NCLrg=0.; NCLrh=0.; NMLsr=0.; NMLgr=0.
1787 NMLir= 0.; NSHhr=0.; NNUvi=0.; NVDvi=0.; NVDvh=0.; QCLig=0.
1788 NCLir= 0.; NCLis=0.; NCLig=0.; NCLih=0.; NIMsi=0.; NIMgi=0.
1789 NiCNis=0.; NsCNis=0.; NVDvs=0.; NCNsg=0.; NCLgr=0.; NCLsrh=0.
1790 NCLss= 0.; NCLsr=0.; NCLsh=0.; NCLsrs=0.; NCLgrg=0.; NgCNgh=0.
1791 NVDvg= 0.; NCLirg=0.; NCLsrg=0.; NCLgrh=0.; NrFZrh=0.; NhFZrh=0.
1792 NCLirh=0.
1794 Dirg=0.; Dirh=0.; Dsrs= 0.; Dsrg= 0.; Dsrh= 0.; Dgrg=0.; Dgrh=0.
1796 !-------------------------------------------------------------------------------------------!
1798 Si = Q(i,k)/QSI(i,k)
1799 iABi = 1./( CHLS*CHLS/(Ka*RGASV*T(i,k)**2) + 1./(DE(i,k)*(QSI(i,k))*Cdiff) )
1801 ! COLLECTION by snow, graupel, hail:
1802 ! (i.e. wet or dry ice-categories [=> excludes ice crystals])
1804 ! Collection by SNOW:
1805 if (QN(i,k)>epsQ) then
1806 ! cloud:
1807 if (QC(i,k)>epsQ) then
1809 !Approximation of Ecs based on Pruppacher & Klett (1997) Fig. 14-11
1810 Ecs= min(Dc,30.e-6)*3.333e+4*sqrt(min(Ds,1.e-3)*1.e+3)
1811 QCLcs= dt*gam*afs*cmr*Ecs*PIov4*iDE(i,k)*(NC(i,k)*NN(i,k))*iGC5*iGS31* &
1812 (GC13*GS13*iLAMc3*iLAMsB2+2.*GC14*GS12*iLAMc4*iLAMsB1+GC15*GS11* &
1813 iLAMc5*iLAMsB0)
1815 NCLcs= dt*gam*afs*PIov4*Ecs*(NC(i,k)*NN(i,k))*iGC5*iGS31*(GC5*GS13* &
1816 iLAMsB2+2.*GC11*GS12*iLAMc*iLAMsB1+GC12*GS11*iLAMc2*iLAMsB0)
1818 !continuous collection: (alternative; gives values ~0.95 of SCE [above])
1819 !QCLcs= dt*gam*Ecs*PIov4*afs*QC(i,k)*NN(i,k)*iLAMs**(2.+bfs)*GS13*iGS31
1820 !NCLcs= QCLcs*NC(i,k)/QC(i,k)
1822 !Correction factor for non-spherical snow [D = maximum dimension] which
1823 !changes projected area: [assumption: A=0.50*D**2 (vs. A=(PI/4)*D**2)]
1824 ! note: Strictly speaking, this correction should only be applied to
1825 ! continuous growth approximation for cloud. [factor = 0.50/(pi/4)]
1826 if (.not. snowSpherical) then
1827 tmp1 = 0.6366 !factor = 0.50/(pi/4)
1828 QCLcs= tmp1*QCLcs
1829 NCLcs= tmp1*NCLcs
1830 endif
1832 QCLcs= min(QCLcs, QC(i,k))
1833 NCLcs= min(NCLcs, NC(i,k))
1834 else
1835 QCLcs= 0.; NCLcs= 0.
1836 endif
1838 ! ice:
1839 if (QI(i,k)>epsQ) then
1840 tmp1= vs0-vi0
1841 tmp3= sqrt(tmp1*tmp1+0.04*vs0*vi0)
1843 QCLis= dt*cmi*iDE(i,k)*PI*6.*Eis*(NY(i,k)*NN(i,k))*tmp3*iGI31*iGS31*(0.5* &
1844 iLAMs2*iLAMi3+2.*iLAMs*iLAMi4+5.*iLAMi5)
1846 NCLis= dt*PIov4*Eis*(NY(i,k)*NN(i,k))*GI31*GS31*tmp3*(GI33*GS31*iLAMi2+ &
1847 2.*GI32*GS32*iLAMi*iLAMs+GI31*GS33*iLAMs2)
1849 QCLis= min(QCLis, (QI(i,k)))
1850 NCLis= min(QCLis*(NY(i,k)*iQI), NCLis)
1851 else
1852 QCLis= 0.; NCLis= 0.
1853 endif
1855 !snow: (i.e. self-collection [aggregation])
1856 NCLss= dt*0.93952*Ess*(DE(i,k)*(QN(i,k)))**((2.+bfs)*thrd)*(NN(i,k))** &
1857 ((4.-bfs)*thrd)
1858 !Note: 0.91226 = I(bfs)*afs*PI^((1-bfs)/3)*des^((-2-bfs)/3); I(bfs=0.41)=1138
1859 ! 0.93952 = I(bfs)*afs*PI^((1-bfs)/3)*des^((-2-bfs)/3); I(bfs=0.42)=1172
1860 ! [interpolated from 3rd-order polynomial approx. of values given in RRB98;
1861 ! see eqn(A.35)]
1862 NCLss= min(NCLss, 0.5*(NN(i,k)))
1864 else
1865 QCLcs= 0.; NCLcs= 0.; QCLis= 0.; NCLis= 0.; NCLss= 0.
1866 endif
1868 ! Collection by GRAUPEL:
1869 if (QG(i,k)>epsQ) then
1871 ! cloud:
1872 if (QC(i,k)>epsQ) then
1874 !(parameterization of Ecg based on Cober and List, 1993 [JAS])
1875 Kstoke = dew*vg0*Dc*Dc/(9.*MUdyn*Dg)
1876 Kstoke = max(1.5,min(10.,Kstoke))
1877 Ecg = 0.55*log10(2.51*Kstoke)
1879 QCLcg= dt*gam*afg*cmr*Ecg*PIov4*iDE(i,k)*(NC(i,k)*NG(i,k))*iGC5*iGG31* &
1880 (GC13*GG13*iLAMc3*iLAMgB2+ 2.*GC14*GG12*iLAMc4*iLAMgB1+GC15*GG11* &
1881 iLAMc5*iLAMgB0)
1883 NCLcg= dt*gam*afg*PIov4*Ecg*(NC(i,k)*NG(i,k))*iGC5*iGG31*(GC5*GG13* &
1884 iLAMgB2+2.*GC11*GG12*iLAMc*iLAMgB1+GC12*GG11*iLAMc2*iLAMgB0)
1886 QCLcg= min(QCLcg, (QC(i,k)))
1887 NCLcg= min(NCLcg, (NC(i,k)))
1888 else
1889 QCLcg= 0.; NCLcg= 0.
1890 endif
1892 ! ice:
1893 if (QI(i,k)>epsQ) then
1894 tmp1= vg0-vi0
1895 tmp3= sqrt(tmp1*tmp1+0.04*vg0*vi0)
1897 QCLig= dt*cmi*iDE(i,k)*PI*6.*Eig*(NY(i,k)*NG(i,k))*tmp3*iGI31*iGG31*(0.5* &
1898 iLAMg2*iLAMi3+2.*iLAMg*iLAMi4+5.*iLAMi5)
1899 NCLig= dt*PIov4*Eig*(NY(i,k)*NG(i,k))*GI31*GG31*tmp3*(GI33*GG31*iLAMi2+ &
1900 2.*GI32*GG32*iLAMi*iLAMg+GI31*GG33*iLAMg2)
1902 QCLig= min(QCLig, (QI(i,k)))
1903 NCLig= min(QCLig*(NY(i,k)*iQI), NCLig)
1904 else
1905 QCLig= 0.; NCLig= 0.
1906 endif
1908 !Deposition/sublimation:
1909 VENTg= Avx*GG32*iLAMg*iLAMg+Bvx*ScTHRD*sqrt(gam*afg*iMUkin)*GG09*iLAMg** &
1910 (2.5+0.5*bfg+alpha_g)
1911 ! QVDvg = dt*iDE(i,k)*iABi*(PI2*(Si-1.)*No_g*VENTg - CHLS*CHLF/(Ka*RGASV* &
1912 ! T(i,k)**2)*QCLcg*idt)
1913 QVDvg = dt*iDE(i,k)*iABi*(PI2*(Si-1.)*No_g*VENTg) !neglect accretion term
1914 ! Prevent overdepletion of vapor:
1915 VDmax = (Q(i,k)-QSI(i,k))/(1.+ck6*QSI(i,k)/(T(i,k)-7.66)**2) !KY97_A.33
1916 if(Si>=1.) then
1917 QVDvg= min(max(QVDvg,0.),VDmax)
1918 else
1919 if (VDmax<0.) QVDvg= max(QVDvg,VDmax)
1920 !IF prevents subl.(QVDvs<0 at t) changing to dep.(VDmax>0 at t*)
1921 endif
1922 !NVDvg = -min(0.,NG(i,k)*iQG*QVDvg) !assume slope does not change during sublimation (pos. quantity)
1923 NVDvg = 0. !assume number does not change during sublimation
1925 else
1926 QCLcg= 0.; QCLrg= 0.; QCLig= 0.
1927 NCLcg= 0.; NCLrg= 0.; NCLig= 0.
1928 endif
1930 ! Collection by HAIL:
1931 if (QH(i,k)>epsQ) then
1933 ! cloud:
1934 if (QC(i,k)>epsQ) then
1935 Ech = exp(-8.68e-7*Dc**(-1.6)*Dh) !Ziegler (1985) A24
1937 QCLch= dt*gam*afh*cmr*Ech*PIov4*iDE(i,k)*(NC(i,k)*NH(i,k))*iGC5*iGH31* &
1938 (GC13*GH13*iLAMc3*iLAMhB2+2.*GC14*GH12*iLAMc4*iLAMhB1+GC15*GH11* &
1939 iLAMc5*iLAMhB0)
1941 NCLch= dt*gam*afh*PIov4*Ech*(NC(i,k)*NH(i,k))*iGC5*iGH31*(GC5*GH13* &
1942 iLAMhB2+2.*GC11*GH12*iLAMc*iLAMhB1+GC12*GH11*iLAMc2*iLAMhB0)
1944 QCLch= min(QCLch, QC(i,k))
1945 NCLch= min(NCLch, NC(i,k))
1946 else
1947 QCLch= 0.; NCLch= 0.
1948 endif
1950 ! rain:
1951 if (QR(i,k)>epsQ) then
1952 tmp1= vh0-vr0
1953 tmp3= sqrt(tmp1*tmp1+0.04*vh0*vr0)
1954 QCLrh= dt*cmr*Erh*PIov4*iDE(i,k)*(NH(i,k)*NR(i,k))*iGR31*iGH31*tmp3* &
1955 (GR36*GH31*iLAMr5+2.*GR35*GH32*iLAMr4*iLAMh+GR34*GH33*iLAMr3*iLAMh2)
1957 NCLrh= dt*PIov4*Erh*(NH(i,k)*NR(i,k))*iGR31*iGH31*tmp3*(GR33*GH31* &
1958 iLAMr2+2.*GR32*GH32*iLAMr*iLAMh+GR31*GH33*iLAMh2)
1960 QCLrh= min(QCLrh, QR(i,k))
1961 NCLrh= min(NCLrh, QCLrh*(NR(i,k)*iQR))
1962 else
1963 QCLrh= 0.; NCLrh= 0.
1964 endif
1966 ! ice:
1967 if (QI(i,k)>epsQ) then
1968 tmp1 = vh0-vi0
1969 tmp3 = sqrt(tmp1*tmp1+0.04*vh0*vi0)
1971 QCLih= dt*cmi*iDE(i,k)*PI*6.*Eih*(NY(i,k)*NH(i,k))*tmp3*iGI31*iGH31*(0.5* &
1972 iLAMh2*iLAMi3+2.*iLAMh*iLAMi4+5.*iLAMi5)
1974 NCLih= dt*PIov4*Eih*(NY(i,k)*NH(i,k))*GI31*GH31*tmp3*(GI33*GH31*iLAMi2+ &
1975 2.*GI32*GH32*iLAMi*iLAMh+GI31*GH33*iLAMh2)
1977 QCLih= min(QCLih, QI(i,k))
1978 NCLih= min(QCLih*(NY(i,k)*iQI), NCLih)
1979 else
1980 QCLih= 0.; NCLih= 0.
1981 endif
1983 ! snow:
1984 if (QN(i,k)>epsQ) then
1985 tmp1 = vh0-vs0
1986 tmp3 = sqrt(tmp1*tmp1+0.04*vh0*vs0)
1987 tmp4 = iLAMs2*iLAMs2
1989 if (snowSpherical) then
1990 !hardcoded for dms=3:
1991 QCLsh= dt*cms*iDE(i,k)*PI*6.*Esh*(NN(i,k)*NH(i,k))*tmp3*iGS31*iGH31* &
1992 (0.5*iLAMh2*iLAMs2*iLAMs+2.*iLAMh*tmp4+5.*tmp4*iLAMs)
1993 else
1994 !hardcoded for dms=2:
1995 QCLsh= dt*cms*iDE(i,k)*PI*0.25*Esh*tmp3*NN(i,k)*NH(i,k)*iGS31*iGH31* &
1996 (GH33*GS33*iLAMh**2.*iLAMs**2. + 2.*GH32*GS34*iLAMh*iLAMs**3. + &
1997 GH31*GS35*iLAMs**4.)
1998 endif
2000 NCLsh= dt*PIov4*Esh*(NN(i,k)*NH(i,k))*GS31*GH31*tmp3*(GS33*GH31*iLAMs2+ &
2001 2.*GS32*GH32*iLAMs*iLAMh+GS31*GH33*iLAMh2)
2003 QCLsh= min(QCLsh, (QN(i,k)))
2004 NCLsh= min((NN(i,k)*iQN)*QCLsh, NCLsh, (NN(i,k)))
2005 else
2006 QCLsh= 0.; NCLsh= 0.
2007 endif
2009 !wet growth:
2010 VENTh= Avx*GH32*iLAMh**(2.+alpha_h) + Bvx*ScTHRD*sqrt(gam*afh*iMUkin)*GH09* &
2011 iLAMh**(2.5+0.5*bfh+alpha_h)
2012 QHwet= max(0., dt*PI2*(DE(i,k)*CHLC*Cdiff*DELqvs-Ka*Tc)*No_h*iDE(i,k)/(CHLF+ &
2013 CPW*Tc)*VENTh+(QCLih*iEih+QCLsh*iEsh)*(1.-CPI*Tc/(CHLF+CPW*Tc)) )
2015 !Deposition/sublimation:
2016 ! QVDvh = dt*iDE(i,k)*iABi*(PI2*(Si-1.)*No_h*VENTh - CHLS*CHLF/(Ka*RGASV* &
2017 ! T(i,k)**2)*QCLch*idt)
2018 QVDvh = dt*iDE(i,k)*iABi*(PI2*(Si-1.)*No_h*VENTh) !neglect acretion term
2019 !prevent overdepletion of vapor:
2020 VDmax = (Q(i,k)-QSI(i,k))/(1.+ck6*(QSI(i,k))/(T(i,k)-7.66)**2) !KY97_A.33 ** USED BY OTHERS; COULD BE PUT ABOVE
2021 if(Si>=1.) then
2022 QVDvh= min(max(QVDvh,0.),VDmax)
2023 else
2024 if (VDmax<0.) QVDvh= max(QVDvh,VDmax) !prevents subl.(QVDvs<0 at t) changing to dep.(VDmax>0 at t*)
2025 endif
2026 ! NVDvh= -min(0.,NH(i,k)*iQH*QVDvh) !assume SLOPE does not change during sublimation (pos. quantity)
2027 NVDvh= 0. !assume NUMBER does not change during sublimation
2029 else
2030 QCLch= 0.; QCLrh= 0.; QCLih= 0.; QCLsh= 0.; QHwet= 0.
2031 NCLch= 0.; NCLrh= 0.; NCLsh= 0.; NCLih= 0.
2032 endif
2034 IF (T(i,k)>TRPL .and. warmphase_ON) THEN
2035 !**********!
2036 ! T > To !
2037 !**********!
2039 ! MELTING of frozen particles:
2040 ! ICE:
2041 QMLir = QI(i,k) !all pristine ice melts in one time step
2042 QI(i,k)= 0.
2043 NMLir = NY(i,k)
2045 ! SNOW:
2046 if (QN(i,k)>epsQ) then
2047 QMLsr= dt*(PI2*iDE(i,k)*iCHLF*No_s*VENTs*(Ka*Tc-CHLC*Cdiff*DELqvs) + CPW* &
2048 iCHLF*Tc*(QCLcs+QCLrs)*idt)
2049 QMLsr= min(max(QMLsr,0.), QN(i,k))
2050 NMLsr= NN(i,k)*iQN*QMLsr
2051 else
2052 QMLsr= 0.; NMLsr= 0.
2053 endif
2055 ! GRAUPEL:
2056 if (QG(i,k)>epsQ) then
2057 QMLgr= dt*(PI2*iDE(i,k)*iCHLF*No_g*VENTg*(Ka*Tc-CHLC*Cdiff*DELqvs) + CPW* &
2058 iCHLF*Tc*(QCLcg+QCLrg)*idt)
2059 QMLgr= min(max(QMLgr,0.), QG(i,k))
2060 NMLgr= NG(i,k)*iQG*QMLgr
2061 else
2062 QMLgr= 0.; NMLgr= 0.
2063 endif
2065 ! HAIL:
2066 if (QH(i,k)>epsQ.and.Tc>5.) then
2067 VENTh= Avx*GH32*iLAMh**(2.+alpha_h) + Bvx*ScTHRD*sqrt(gam*afh*iMUkin)*GH09* &
2068 iLAMh**(2.5+0.5*bfh+alpha_h)
2069 QMLhr= dt*(PI2*iDE(i,k)*iCHLF*No_h*VENTh*(Ka*Tc-CHLC*Cdiff*DELqvs) + CPW/ &
2070 CHLF*Tc*(QCLch+QCLrh)*idt)
2071 QMLhr= min(max(QMLhr,0.), QH(i,k))
2072 NMLhr= NH(i,k)*iQH*QMLhr
2073 if(QCLrh>0.) NMLhr= NMLhr*0.1 !Prevents problems when hail is ML & CL
2074 else
2075 QMLhr= 0.; NMLhr= 0.
2076 endif
2078 ! Cold (sub-zero) source/sink terms:
2079 QNUvi= 0.; QFZci= 0.; QVDvi= 0.; QVDvs= 0.
2080 QCLis= 0.; QCNis1=0.; QCNis2=0.; QCLri= 0.
2081 QCNgh= 0.; QIMsi= 0.; QIMgi= 0.; QCLir= 0.
2082 QCLrs= 0.; QCLgr= 0.; QCLrg= 0.; QCNis= 0.
2083 QCNsg= 0.; QCLsr= 0.
2085 NNUvi= 0.; NFZci= 0.; NCLgr= 0.; NCLrg= 0.; NgCNgh= 0.
2086 NCLis= 0.; NVDvi= 0.; NVDvs= 0.; NCLri= 0.; NCLsr= 0.
2087 NCNsg= 0.; NhCNgh= 0.; NiCNis=0.; NsCNis=0.
2088 NIMsi= 0.; NIMgi= 0.; NCLir= 0.; NCLrs= 0.
2090 ELSE !----------!
2091 ! T < To !
2092 !----------!
2094 ! Warm-air-only source/sink terms:
2095 QMLir= 0.; QMLsr= 0.; QMLgr= 0.; QMLhr= 0.
2096 NMLir= 0.; NMLsr= 0.; NMLgr= 0.; NMLhr= 0.
2098 !Probabilistic freezing (Bigg) of rain:
2099 if (Tc<Tc_FZrh .and. QR(i,k)>epsQ .and. hail_ON) then
2100 !note: - (Tc<-10.C) condition is based on Pruppacher-Klett (1997) Fig. 9-41
2101 ! - Small raindrops will freeze to hail. However, if after all S/S terms
2102 ! are added Dh<Dh_min, then hail will be converted to graupel. Thus,
2103 ! probabilistic freezing of small rain is effectively a source of graupel.
2104 NrFZrh= -dt*Bbigg*(exp(Abigg*Tc)-1.)*DE(i,k)*QR(i,k)*idew
2105 Rz= 1. !N and Z (and Q) are conserved for FZrh with triple-moment
2106 ! The Rz factor serves to conserve reflectivity when a rain distribution
2107 ! converts to an distribution with a different shape parameter, alpha.
2108 ! (e.g. when rain freezes to hail) The factor Rz non-conserves N while
2109 ! acting to conserve Z for double-moment. See Ferrier, 1994 App. D)
2110 ! Rz= (gamma(7.d0+alpha_h)*GH31*GR34*GR34)/(GR36(i,k)*GR31* &
2111 ! gamma(4.d0+alpha_h)*gamma(4.d0+alpha_h))
2112 NhFZrh= Rz*NrFZrh
2113 QFZrh = NrFZrh*(QR(i,k)*iNR)
2114 else
2115 QFZrh= 0.; NrFZrh= 0.; NhFZrh= 0.
2116 endif
2118 !--------!
2119 ! ICE: !
2120 !--------!
2121 ! Homogeneous freezing of cloud to ice:
2122 if (QC(i,k)>epsQ) then
2123 tmp2 = Tc*Tc; tmp3= tmp2*Tc; tmp4= tmp2*tmp2
2124 JJ = (10.**max(-20.,(-606.3952-52.6611*Tc-1.7439*tmp2-0.0265*tmp3- &
2125 1.536e-4*tmp4)))
2126 tmp1 = 1.e6*(DE(i,k)*(QC(i,k)*iNC)*icmr) !i.e. Dc[cm]**3
2127 FRAC = 1.-exp(-JJ*PIov6*tmp1*dt)
2128 if (Tc>-30.) FRAC= 0.
2129 if (Tc<-50.) FRAC= 1.
2130 QFZci= FRAC*QC(i,k)
2131 NFZci= FRAC*NC(i,k)
2132 else
2133 QFZci= 0.; NFZci= 0.
2134 endif
2136 !Primary ice nucleation:
2137 NNUvi= 0.; QNUvi= 0.
2138 if (primIceNucl==1) then
2140 NuDEPSOR= 0.; NuCONT= 0.
2141 if (QSI(i,k)>1.e-20) then
2142 Simax = min(Si, satw_peak*QSW(i,k)/QSI(i,k))
2143 else
2144 Simax = 0.
2145 endif
2146 tmp1 = T(i,k)-7.66
2147 NNUmax = max(0., DE(i,k)/mio*(Q(i,k)-QSI(i,k))/(1.+ck6*(QSI(i,k)/(tmp1* &
2148 tmp1))))
2149 !Deposition/sorption nucleation:
2150 if (Tc<-5. .and. Si>1.) then
2151 NuDEPSOR= max(0., 1.e3*exp(12.96*(Simax-1.)-0.639)-(NY(i,k))) !Meyers(1992)
2152 endif
2153 !Contact nucleation:
2154 if (QC(i,k)>epsQ .and. Tc<-2. .and. WZ(i,k)>0.001) then
2155 GG = 1.*idew/(RGASV*(T(i,k))/((QSW(i,k)*pres(i,k))/EPS1)/ &
2156 Cdiff+CHLC/Ka/(T(i,k))*(CHLC/RGASV/(T(i,k))-1.)) !CP00a
2157 Swmax = SxFNC(WZ(i,k),Tc,pres(i,k),QSW(i,k),QSI(i,k),CCNtype,1)
2158 if (QSW(i,k)>1.e-20) then
2159 ssat= min((Q(i,k)/QSW(i,k)), Swmax) -1.
2160 else
2161 ssat= 0.
2162 endif
2163 Tcc = Tc + GG*ssat*CHLC/Kdiff !C86_eqn64
2164 Na = exp(4.11-0.262*Tcc) !W95_eqn60/M92_2.6
2165 Kn = LAMa0*(T(i,k))*p0/(T0*pres(i,k)*Ra) !W95_eqn59
2166 PSIa = -kBoltz*Tcc/(6.*pi*Ra*MUdyn)*(1.+Kn) !W95_eqn58
2167 ft = 0.4*(1.+1.45*Kn+0.4*Kn*exp(-1./Kn))*(Ka+2.5*Kn*KAPa)/ &
2168 (1.+3.*Kn)/(2.*Ka+5.*KAPa*Kn+KAPa) !W95_eqn57
2169 Dc = (DE(i,k)*(QC(i,k)*iNC)*icmr)**thrd
2170 F1 = PI2*Dc*Na*(NC(i,k)) !W95_eqn55
2171 F2 = Ka/pres(i,k)*(Tc-Tcc) !W95_eqn56
2172 NuCONTA= -F1*F2*RGASV*(T(i,k))/CHLC*iDE(i,k) !diffusiophoresis
2173 NuCONTB= F1*F2*ft*iDE(i,k) !thermeophoresis
2174 NuCONTC= F1*PSIa !Brownian diffusion
2175 NuCONT = max(0.,(NuCONTA+NuCONTB+NuCONTC)*dt)
2176 endif
2177 !Total primary ice nucleation:
2178 if (icephase_ON) then
2179 NNUvi= min(NNUmax, NuDEPSOR + NuCONT )
2180 QNUvi= mio*iDE(i,k)*NNUvi
2181 QNUvi= min(QNUvi,(Q(i,k)))
2182 endif
2184 elseif (primIceNucl==2) then
2185 if (Tc<-5. .and. Si>1.08) then !following Thompson etal (2006)
2186 NNUvi= max(N_Cooper(TRPL,T(i,k))-NY(i,k),0.)
2187 QNUvi= min(mio*iDE(i,k)*NNUvi, Q(i,k))
2188 endif
2189 !elseif (primIceNucl==3) then
2190 !! (for alternative [future] ice nucleation parameterizations)
2191 ! NNUvi=...
2192 ! QNUvi=...
2193 endif !if (primIceNucl==1)
2197 IF (QI(i,k)>epsQ) THEN
2199 !Deposition/sublimation:
2200 ! No_i = NY(i,k)*iGI31/iLAMi**(1.+alpha_i)
2201 ! VENTi= Avx*GI32*iLAMi**(2.+alpha_i)+Bvx*ScTHRD*sqrt(gam*afi*iMUkin)*GI6* &
2202 ! iLAMi**(2.5+0.5*bfi+alpha_i)
2203 No_i = NY(i,k)*iGI31/iLAMi !optimized for alpha_i=0
2204 VENTi= Avx*GI32*iLAMi*iLAMi+Bvx*ScTHRD*sqrt(gam*afi*iMUkin)*GI6*iLAMi** &
2205 (2.5+0.5*bfi+alpha_i)
2206 !Note: ice crystal capacitance is implicitly C = 0.5*D*capFact_i
2207 ! QVDvi= dt*capFact_i*iABi*(PI2*(Si-1.)*No_i*VENTi)
2208 QVDvi = dt*iDE(i,k)*capFact_i*iABi*(PI2*(Si-1.)*No_i*VENTi)
2210 ! Prevent overdepletion of vapor:
2211 VDmax = (Q(i,k)-QSI(i,k))/(1.+ck6*(QSI(i,k))/(T(i,k)-7.66)**2)
2212 if(Si>=1.) then
2213 QVDvi= min(max(QVDvi,0.),VDmax)
2214 else
2215 if (VDmax<0.) QVDvi= max(QVDvi,VDmax)
2216 !IF prevents subl.(QVDvi<0 at t) changing to dep.(VDmax>0 at t*) 2005-06-28
2217 endif
2218 if (.not. iceDep_ON) QVDvi= 0. !suppresses depositional growth
2219 NVDvi= min(0., (NY(i,k)*iQI)*QVDvi) !dNi/dt=0 for deposition
2221 ! Conversion to snow:
2222 if (QI(i,k)+QVDvi>epsQ .and. NY(i,k)+NVDvi>epsN) then
2223 tmp5 = iLAMi !hold value
2224 tmp6 = No_i !hold value
2225 !estimate ice PSD after VDvi (if there were no CNis):
2226 tmp1 = QI(i,k) + QVDvi
2227 tmp2 = NY(i,k) + NVDvi
2228 tmp3 = 1./tmp2
2229 iLAMi = max( iLAMmin2, iLAMDA_x(DE(i,k),tmp1,tmp3,icexi9,thrd) )
2230 No_i = tmp2*iGI31/iLAMi !optimized for alpha_i=0
2231 !compute number and mass of ice converted to snow as the integral from
2232 ! Dso to INF of Ni(D)dD and m(D)Ni(D)dD, respectively:
2233 tmp4 = exp(-Dso/iLAMi)
2234 NiCNis = No_i*iLAMi*tmp4
2235 NsCNis = NiCNis
2236 QCNis = cmi*No_i*tmp4*(Dso**3*iLAMi + 3.*Dso**2*iLAMi**2 + 6.*Dso* &
2237 iLAMi**3 + 6.*iLAMi**4)
2238 iLAMi = tmp5 !(restore value)
2239 No_i = tmp6 !(restore value)
2240 endif
2242 if (.not.(snow_ON)) then !Suppress SNOW initiation (for testing only)
2243 QCNis = 0.
2244 NiCNis = 0.
2245 NsCNis = 0.
2246 endif
2248 ! 3-component freezing (collisions with rain):
2249 if (QR(i,k)>epsQ .and. QI(i,k)>epsQ) then
2250 tmp1 = vr0-vi0
2251 tmp3 = sqrt(tmp1*tmp1+0.04*vr0*vi0)
2253 QCLir= dt*cmi*Eri*PIov4*iDE(i,k)*(NR(i,k)*NY(i,k))*iGI31*iGR31*tmp3* &
2254 (GI36*GR31*iLAMi5+2.*GI35*GR32*iLAMi4*iLAMr+GI34*GR33*iLAMi3* &
2255 iLAMr2)
2257 NCLri= dt*PIov4*Eri*(NR(i,k)*NY(i,k))*iGI31*iGR31*tmp3*(GI33*GR31* &
2258 iLAMi2+2.*GI32*GR32*iLAMi*iLAMr+GI31*GR33*iLAMr2)
2260 QCLri= dt*cmr*Eri*PIov4*iDE(i,k)*(NY(i,k)*NR(i,k))*iGR31*iGI31*tmp3* &
2261 (GR36*GI31 *iLAMr5+2.*GR35*GI32*iLAMr4*iLAMi+GR34*GI33*iLAMr3* &
2262 iLAMi2)
2264 !note: For explicit eqns, both NCLri and NCLir are mathematically identical)
2265 NCLir= min(QCLir*(NY(i,k)*iQI), NCLri)
2266 QCLri= min(QCLri, (QR(i,k))); QCLir= min(QCLir, (QI(i,k)))
2267 NCLri= min(NCLri, (NR(i,k))); NCLir= min(NCLir, (NY(i,k)))
2269 !Determine destination of 3-comp.freezing:
2270 tmp1= max(Di,Dr)
2271 dey= (dei*Di*Di*Di+dew*Dr*Dr*Dr)/(tmp1*tmp1*tmp1)
2272 if (dey>0.5*(deg+deh) .and. Dr>Dr_3cmpThrs .and. hail_ON) then
2273 Dirg= 0.; Dirh= 1.
2274 else
2275 Dirg= 1.; Dirh= 0.
2276 endif
2277 if (.not. grpl_ON) Dirg= 0.
2279 else
2280 QCLir= 0.; NCLir= 0.; QCLri= 0.
2281 NCLri= 0.; Dirh = 0.; Dirg= 0.
2282 endif
2284 ! Rime-splintering (ice multiplication):
2285 ff= 0.
2286 if(Tc>=-8..and.Tc<=-5.) ff= 3.5e8*(Tc +8.)*thrd
2287 if(Tc> -5..and.Tc< -3.) ff= 3.5e8*(-3.-Tc)*0.5
2288 NIMsi= DE(i,k)*ff*QCLcs
2289 NIMgi= DE(i,k)*ff*QCLcg
2290 QIMsi= mio*iDE(i,k)*NIMsi
2291 QIMgi= mio*iDE(i,k)*NIMgi
2293 ELSE
2295 QVDvi= 0.; QCNis= 0.
2296 QIMsi= 0.; QIMgi= 0.; QCLri= 0.; QCLir= 0.
2297 NVDvi= 0.; NCLir= 0.; NIMsi= 0.
2298 NiCNis=0.; NsCNis=0.; NIMgi= 0.; NCLri= 0.
2300 ENDIF
2301 !---------!
2302 ! SNOW: !
2303 !---------!
2304 IF (QN(i,k)>epsQ) THEN
2306 !Deposition/sublimation:
2307 !note: - snow crystal capacitance is implicitly C = 0.5*D*capFact_s
2308 ! - No_s and VENTs are computed above
2309 ! QVDvs = dt*capFact_s*iABi*(PI2*(Si-1.)*No_s*VENTs - CHLS*CHLF/(Ka*RGASV* &
2310 ! T(i,k)*T(i,k))*QCLcs*idt)
2311 QVDvs = dt*iDE(i,k)*capFact_s*iABi*(PI2*(Si-1.)*No_s*VENTs - CHLS*CHLF/(Ka* &
2312 RGASV*T(i,k)**2)*QCLcs*idt)
2314 ! Prevent overdepletion of vapor:
2315 VDmax = (Q(i,k)-QSI(i,k))/(1.+ck6*(QSI(i,k))/(T(i,k)-7.66)**2)
2316 if(Si>=1.) then
2317 QVDvs= min(max(QVDvs,0.),VDmax)
2318 else
2319 if (VDmax<0.) QVDvs= max(QVDvs,VDmax)
2320 !IF prevents subl.(QVDvs<0 at t) changing to dep.(VDmax>0 at t*)
2321 endif
2322 NVDvs= -min(0.,(NN(i,k)*iQN)*QVDvs) !pos. quantity
2324 ! Conversion to graupel:
2325 if (QCLcs>0. .and. QCLcs>CNsgThres*QVDvs .and. grpl_ON) then
2326 tmp1 = 100. !tuning factor: controls amount of mass either added or removed
2327 !from snow (QN) during partial conversion to graupel. [If QCLcs/QN > 1/tmp1,
2328 !then some snow mass will be converted to graupel (in addition to rime mass).]
2329 QCNsg = min( QN(i,k)+QCLcs, QCLcs*(tmp1*QCLcs/QN(i,k)) )
2330 !calculate NCNsg: [explicit logic]
2331 !mgo = DE(i,k)*(QN(i,k)+QCLsg)/NN(i,k) !mean-mass of new graupel
2332 !NCNsg = DE(i,k)*QCNsg/mgo
2333 !calculate NCNsg: [optimized; substituting mgo and factoring]
2334 NCNsg = DE(i,k)*QCNsg/(QN(i,k)+QCLcs)
2335 else
2336 QCNsg = 0.
2337 NCNsg = 0.
2338 endif
2340 ! 3-component freezing (collisions with rain):
2341 if (QR(i,k)>epsQ .and. QN(i,k)>epsQ .and. Tc<-5.) then
2342 tmp1 = vs0-vr0
2343 tmp2 = sqrt(tmp1*tmp1+0.04*vs0*vr0)
2344 tmp6 = iLAMs2*iLAMs2*iLAMs
2346 QCLrs= dt*cmr*Ers*PIov4*iDE(i,k)*NN(i,k)*NR(i,k)*iGR31*iGS31*tmp2* &
2347 (GR36*GS31*iLAMr5+2.*GR35*GS32*iLAMr4*iLAMs+GR34*GS33*iLAMr3* &
2348 iLAMs2)
2350 NCLrs= dt*0.25e0*PI*Ers*(NN(i,k)*NR(i,k))*iGR31*iGS31*tmp2*(GR33* &
2351 GS31*iLAMr2+2.*GR32*GS32*iLAMr*iLAMs+GR31*GS33*iLAMs2)
2353 if (snowSpherical) then
2354 !hardcoded for dms=3:
2355 QCLsr= dt*cms*Ers*PIov4*iDE(i,k)*(NR(i,k)*NN(i,k))*iGS31*iGR31* &
2356 tmp2*(GS36*GR31*tmp6+2.*GS35*GR32*iLAMs2*iLAMs2*iLAMr+GS34* &
2357 GR33*iLAMs2*iLAMs*iLAMr2)
2358 else
2359 !hardcoded for dms=2:
2360 QCLsr= dt*cms*iDE(i,k)*PI*0.25*ERS*tmp2*NN(i,k)*NR(i,k)*iGS31* &
2361 iGR31*(GR33*GS33*iLAMr**2.*iLAMs**2. + 2.*GR32*GS34*iLAMr* &
2362 iLAMs**3. +GR31*GS35*iLAMs**4.)
2363 endif
2365 !note: For explicit eqns, NCLsr = NCLrs
2366 NCLsr= min(QCLsr*(NN(i,k)*iQN), NCLrs)
2367 QCLrs= min(QCLrs, QR(i,k)); QCLsr= min(QCLsr, QN(i,k))
2368 NCLrs= min(NCLrs, NR(i,k)); NCLsr= min(NCLsr, NN(i,k))
2370 ! Determine destination of 3-comp.freezing:
2371 Dsrs= 0.; Dsrg= 0.; Dsrh= 0.
2372 tmp1= max(Ds,Dr)
2373 tmp2= tmp1*tmp1*tmp1
2374 dey = (des*Ds*Ds*Ds + dew*Dr*Dr*Dr)/tmp2
2375 if (dey<=0.5*(des+deg) ) Dsrs= 1. !snow
2376 if (dey >0.5*(des+deg) .and. dey<0.5*(deg+deh)) Dsrg= 1. !graupel
2377 if (dey>=0.5*(deg+deh)) then
2378 Dsrh= 1. !hail
2379 if (.not.hail_ON .or. Dr<Dr_3cmpThrs) then
2380 Dsrg= 1.; Dsrh= 0. !graupel
2381 endif
2382 endif
2383 if (.not. grpl_ON) Dsrg=0.
2385 else
2386 QCLrs= 0.; QCLsr= 0.; NCLrs= 0.; NCLsr= 0.
2387 endif
2389 ELSE
2391 QVDvs= 0.; QCLcs= 0.; QCNsg= 0.; QCLsr= 0.; QCLrs= 0.
2392 NVDvs= 0.; NCLcs= 0.; NCLsr= 0.; NCLrs= 0.; NCNsg= 0.
2394 ENDIF
2395 !------------!
2396 ! GRAUPEL: !
2397 !------------!
2398 IF (QG(i,k)>epsQ .and. NG(i,k)>epsN) THEN
2400 !Conversion to hail: (D_sll given by S-L limit)
2401 if ( (QCLcg+QCLrg)>0. .and. hail_ON ) then
2402 ! D_sll = 0.01*(exp(min(20.,-Tc/(1.1e4*DE(i,k)*(QC(i,k)+QR(i,k))-1.3e3*DE(i,k)*QI(i,k)+1.)))-1.)
2403 ! D_sll = 2.0*D_sll !correction factor [error Ziegler (1985), as per Young (1993)]
2404 ! D_sll = 2.*0.01*(exp(min(20.,-Tc/(1.1e4*DE(i,k)*(QC(i,k)+QR(i,k)) + 1.)))-1.)
2405 tmp1 = 1.1e4*DE(i,k)*(QC(i,k)+QR(i,k)) + 1.
2406 tmp1 = max(1.,tmp1) !to prevent div-by-zero
2407 D_sll = 2.*0.01*(exp(min(20.,-Tc/tmp1))-1.)
2408 D_sll = min(1., max(0.0001,D_sll)) !smallest D_sll=0.1mm; largest=1m
2409 tmp1 = iLAMg !hold value
2410 tmp2 = No_g !hold value
2411 !estimate PSD after accretion (and before conversion to hail): (assume inverse-exponential)
2412 tmp3 = QG(i,k) + QCLcg + QCLrg
2413 iLAMg = exp(thrd*log(DE(i,k)*tmp3/(NG(i,k)*6*cmg)))
2414 No_g = NG(i,k)/iLAMg
2415 tmp4 = exp(-D_sll/iLAMg)
2416 Ng_tail = No_g*iLAMg*tmp4
2417 if (Ng_tail > Ngh_crit) then
2418 NgCNgh= min(NG(i,k), Ng_tail)
2419 QCNgh = min(QG(i,k), cmg*No_g*tmp4*(D_sll**3*iLAMg + 3.*D_sll**2* &
2420 iLAMg**2 + 6.*D_sll*iLAMg**3 + 6.*iLAMg**4) )
2421 Rz= 1.
2422 ! The Rz factor (/=1) serves to conserve reflectivity when graupel
2423 ! converts to hail with a a different shape parameter, alpha.
2424 ! (See Ferrier, 1994 App. D).
2425 NhCNgh = Rz*NgCNgh
2426 else
2427 QCNgh = 0.
2428 NgCNgh = 0.
2429 NhCNgh = 0.
2430 endif
2431 iLAMg = tmp1 !restore value
2432 No_g = tmp2 !restore value
2433 endif
2435 !3-component freezing (collisions with rain):
2436 ! if (QR(i,k)>epsQ) then
2437 if (QR(i,k)>epsQ .and. Tc<-5.) then
2438 tmp1 = vg0-vr0
2439 tmp2 = sqrt(tmp1*tmp1 + 0.04*vg0*vr0)
2440 tmp8 = iLAMg2*iLAMg ! iLAMg**3
2441 tmp9 = tmp8*iLAMg ! iLAMg**4
2442 tmp10= tmp9*iLAMg ! iLAMg**5
2444 !(parameterization of Erg based on Cober and List, 1993 [JAS])
2445 Kstoke = dew*abs(vg0-vr0)*Dr*Dr/(9.*MUdyn*Dg)
2446 Kstoke = max(1.5,min(10.,Kstoke))
2447 Erg = 0.55*log10(2.51*Kstoke)
2449 QCLrg= dt*cmr*Erg*PIov4*iDE(i,k)*(NG(i,k)*NR(i,k))*iGR31*iGG31*tmp2* &
2450 (GR36*GG31*iLAMr5+2.*GR35*GG32*iLAMr4*iLAMg+GR34*GG33*iLAMr3* &
2451 iLAMg2)
2453 NCLrg= dt*PIov4*Erg*(NG(i,k)*NR(i,k))*iGR31*iGG31*tmp2*(GR33*GG31* &
2454 iLAMr2+2.*GR32*GG32*iLAMr*iLAMg+GR31*GG33*iLAMg2)
2456 QCLgr= dt*cmg*Erg*PIov4*iDE(i,k)*(NR(i,k)*NG(i,k))*iGG31*iGR31*tmp2* &
2457 (GG36*GR31*tmp10+2.*GG35*GR32*tmp9*iLAMr+GG34*GR33*tmp8*iLAMr2)
2459 !(note: For explicit eqns, NCLgr= NCLrg)
2460 NCLgr= min(NCLrg, QCLgr*(NG(i,k)*iQG))
2461 QCLrg= min(QCLrg, QR(i,k)); QCLgr= min(QCLgr, QG(i,k))
2462 NCLrg= min(NCLrg, NR(i,k)); NCLgr= min(NCLgr, NG(i,k))
2464 ! Determine destination of 3-comp.freezing:
2465 tmp1= max(Dg,Dr)
2466 tmp2= tmp1*tmp1*tmp1
2467 dey = (deg*Dg*Dg*Dg + dew*Dr*Dr*Dr)/tmp2
2468 if (dey>0.5*(deg+deh) .and. Dr>Dr_3cmpThrs .and. hail_ON) then
2469 Dgrg= 0.; Dgrh= 1.
2470 else
2471 Dgrg= 1.; Dgrh= 0.
2472 endif
2473 else
2474 QCLgr= 0.; QCLrg= 0.; NCLgr= 0.; NCLrg= 0.
2475 endif
2477 ELSE
2479 QCNgh= 0.; QCLgr= 0.; QCLrg= 0.; NgCNgh= 0.
2480 NhCNgh= 0.; NCLgr= 0.; NCLrg= 0.
2482 ENDIF
2483 !---------!
2484 ! HAIL: !
2485 !---------!
2486 IF (QH(i,k)>epsQ) THEN
2488 !Wet growth:
2489 if (QHwet<(QCLch+QCLrh+QCLih+QCLsh) .and. Tc>-40.) then
2490 QCLih= min(QCLih*iEih, QI(i,k)) !change Eih to 1. in CLih
2491 NCLih= min(NCLih*iEih, NY(i,k)) ! " "
2492 QCLsh= min(QCLsh*iEsh, QN(i,k)) !change Esh to 1. in CLsh
2493 NCLsh= min(NCLsh*iEsh, NN(i,k)) ! " "
2494 tmp3 = QCLrh
2495 QCLrh= QHwet-(QCLch+QCLih+QCLsh) !actual QCLrh minus QSHhr
2496 QSHhr= tmp3-QCLrh !QSHhr used here only
2497 NSHhr= DE(i,k)*QSHhr/(cmr*Drshed*Drshed*Drshed)
2498 else
2499 NSHhr= 0.
2500 endif
2501 ELSE
2502 NSHhr= 0.
2503 ENDIF
2505 ENDIF ! ( if Tc<0C Block )
2507 !----- Prevent mass transfer from accretion during melting: ---!
2508 ! (only if exepcted NX (X=s,g,h) < 0 after source/sinks added)
2509 ! The purpose is to prevent mass tranfer from cloud to X and then
2510 ! to vapor (during "prevent overdepletion" code) if all of X
2511 ! would otherwise be completely depleted (e.g. due to melting).
2513 !estimate NN after S/S:
2514 tmp1= NN(i,k) +NsCNis -NVDvs -NCNsg -NMLsr -NCLss -NCLsr -NCLsh +NCLsrs
2515 if (tmp1<epsN) then
2516 QCLcs = 0.
2517 NCLcs = 0.
2518 QCLrs = 0.
2519 NCLrs = 0.
2520 if (Dsrs==1.) then
2521 QCLrs = 0.
2522 QCLsr = 0.
2523 endif
2524 endif
2525 !estimate NG after S/S:
2526 tmp2= NG(i,k) +NCNsg -NCLgr -NVDvg -NMLgr +NCLirg +NCLsrg +NCLgrg -NgCNgh
2527 if (tmp2<epsN) then
2528 QCLcg = 0.
2529 NCLcg = 0.
2530 QCLrg = 0.
2531 NCLrg = 0.
2532 if (Dirg==1.) then
2533 QCLri = 0.
2534 QCLir = 0.
2535 endif
2536 if (Dsrg==1.) then
2537 QCLrs = 0.
2538 QCLsr = 0.
2539 endif
2540 endif
2541 !estimate NH after S/S:
2542 tmp3= NH(i,k) +NhFZrh +NhCNgh -NMLhr -NVDvh +NCLirh +NCLsrh +NCLgrh
2543 if (tmp3<epsN) then
2544 QCLch = 0.
2545 NCLch = 0.
2546 QCLrh = 0.
2547 NCLrh = 0.
2548 if (Dirh==1.) then
2549 QCLri = 0.
2550 QCLir = 0.
2551 endif
2552 if (Dsrh==1.) then
2553 QCLrs = 0.
2554 QCLsr = 0.
2555 endif
2556 endif
2557 !=====
2559 !------------ End of source/sink term calculation -------------!
2561 !-- Adjustment of source/sink terms to prevent overdepletion: --!
2562 do niter= 1,2
2564 ! (1) Vapor:
2565 source= Q(i,k) +dim(-QVDvi,0.)+dim(-QVDvs,0.)+dim(-QVDvg,0.)+dim(-QVDvh,0.)
2566 sink = QNUvi+dim(QVDvi,0.)+dim(QVDvs,0.)
2567 sour = max(source,0.)
2568 if(sink>sour) then
2569 ratio= sour/sink
2570 QNUvi= ratio*QNUvi; NNUvi= ratio*NNUvi
2571 if(QVDvi>0.) then
2572 QVDvi= ratio*QVDvi; NVDvi= ratio*NVDvi
2573 endif
2574 if(QVDvs>0.) then
2575 QVDvs=ratio*QVDvs; NVDvs=ratio*NVDvs
2576 endif
2577 QVDvg= ratio*QVDvg; NVDvg= ratio*NVDvg
2578 QVDvh= ratio*QVDvh; NVDvh= ratio*NVDvh
2579 endif
2581 ! (2) Cloud:
2582 source= QC(i,k)
2583 sink = QCLcs+QCLcg+QCLch+QFZci
2584 sour = max(source,0.)
2585 if(sink>sour) then
2586 ratio= sour/sink
2587 QFZci= ratio*QFZci; NFZci= ratio*NFZci
2588 QCLcs= ratio*QCLcs; NCLcs= ratio*NCLcs
2589 QCLcg= ratio*QCLcg; NCLcg= ratio*NCLcg
2590 QCLch= ratio*QCLch; NCLch= ratio*NCLch
2591 endif
2593 ! (3) Rain:
2594 source= QR(i,k)+QMLsr+QMLgr+QMLhr+QMLir
2595 sink = QCLri+QCLrs+QCLrg+QCLrh+QFZrh
2596 sour = max(source,0.)
2597 if(sink>sour) then
2598 ratio= sour/sink
2599 QCLrg= ratio*QCLrg; QCLri= ratio*QCLri; NCLri= ratio*NCLri
2600 QCLrs= ratio*QCLrs; NCLrs= ratio*NCLrs
2601 NCLrg= ratio*NCLrg; QCLrh= ratio*QCLrh; NCLrh= ratio*NCLrh
2602 QFZrh= ratio*QFZrh; NrFZrh=ratio*NrFZrh; NhFZrh=ratio*NhFZrh
2603 if (ratio==0.) then
2604 Dirg= 0.; Dirh= 0.; Dgrg= 0.; Dgrh= 0.
2605 Dsrs= 0.; Dsrg= 0.; Dsrh= 0.
2606 endif
2607 endif
2609 ! (4) Ice:
2610 source= QI(i,k)+QNUvi+dim(QVDvi,0.)+QFZci
2611 sink = QCNis+QCLir+dim(-QVDvi,0.)+QCLis+QCLig+QCLih+QMLir
2612 sour = max(source,0.)
2613 if(sink>sour) then
2614 ratio= sour/sink
2615 QMLir= ratio*QMLir; NMLir= ratio*NMLir
2616 if (QVDvi<0.) then
2617 QVDvi= ratio*QVDvi; NVDvi= ratio*NVDvi
2618 endif
2619 QCNis= ratio*QCNis; NiCNis= ratio*NiCNis; NsCNis= ratio*NsCNis
2620 QCLir= ratio*QCLir; NCLir= ratio*NCLir; QCLig= ratio*QCLig
2621 QCLis= ratio*QCLis; NCLis= ratio*NCLis
2622 QCLih= ratio*QCLih; NCLih= ratio*NCLih
2623 if (ratio==0.) then
2624 Dirg= 0.; Dirh= 0.
2625 endif
2626 endif
2628 ! (5) Snow:
2629 source= QN(i,k)+QCNis+dim(QVDvs,0.)+QCLis+Dsrs*(QCLrs+QCLsr)+QCLcs
2630 sink = dim(-QVDvs,0.)+QCNsg+QMLsr+QCLsr+QCLsh
2631 sour = max(source,0.)
2632 if(sink>sour) then
2633 ratio= sour/sink
2634 if(QVDvs<=0.) then
2635 QVDvs= ratio*QVDvs; NVDvs= ratio*NVDvs
2636 endif
2637 QCNsg= ratio*QCNsg; NCNsg= ratio*NCNsg; QMLsr= ratio*QMLsr
2638 NMLsr= ratio*NMLsr; QCLsr= ratio*QCLsr; NCLsr= ratio*NCLsr
2639 QCLsh= ratio*QCLsh; NCLsh= ratio*NCLsh
2640 if (ratio==0.) then
2641 Dsrs= 0.; Dsrg= 0.; Dsrh= 0.
2642 endif
2643 endif
2645 ! (6) Graupel:
2646 source= QG(i,k)+QCNsg+dim(QVDvg,0.)+Dirg*(QCLri+QCLir)+Dgrg*(QCLrg+QCLgr)+ &
2647 QCLcg+Dsrg*(QCLrs+QCLsr)+QCLig
2648 sink = dim(-QVDvg,0.)+QMLgr+QCNgh+QCLgr
2649 sour = max(source,0.)
2650 if(sink>sour) then
2651 ratio= sour/sink
2652 QVDvg= ratio*QVDvg; NVDvg= ratio*NVDvg; QMLgr = ratio*QMLgr
2653 NMLgr= ratio*NMLgr; QCNgh= ratio*QCNgh; NgCNgh= ratio*NgCNgh
2654 QCLgr= ratio*QCLgr; NCLgr= ratio*NCLgr; NhCNgh= ratio*NhCNgh
2655 if (ratio==0.) then
2656 Dgrg= 0.; Dgrh= 0.
2657 endif
2658 endif
2660 ! (7) Hail:
2661 source= QH(i,k)+dim(QVDvh,0.)+QCLch+QCLrh+Dirh*(QCLri+QCLir)+QCLih+QCLsh+ &
2662 Dsrh*(QCLrs+QCLsr)+QCNgh+Dgrh*(QCLrg+QCLgr)+QFZrh
2663 sink = dim(-QVDvh,0.)+QMLhr
2664 sour = max(source,0.)
2665 if(sink>sour) then
2666 ratio= sour/sink
2667 QVDvh= ratio*QVDvh; NVDvh= ratio*NVDvh
2668 QMLhr= ratio*QMLhr; NMLhr= ratio*NMLhr
2669 endif
2671 enddo
2672 !--------------- End of source/sink term adjustment ------------------!
2674 !Compute N-tendencies for destination categories of 3-comp.freezing:
2675 NCLirg= 0.; NCLirh= 0.; NCLsrs= 0.; NCLsrg= 0.
2676 NCLsrh= 0.; NCLgrg= 0.; NCLgrh= 0.
2678 if (QCLir+QCLri>0.) then
2679 tmp1 = max(Dr,Di)
2680 tmp2 = tmp1*tmp1*tmp1*PIov6
2681 NCLirg= Dirg*DE(i,k)*(QCLir+QCLri)/(deg*tmp2)
2682 NCLirh= Dirh*DE(i,k)*(QCLir+QCLri)/(deh*tmp2)
2683 endif
2685 if (QCLsr+QCLrs>0.) then
2686 tmp1 = max(Dr,Ds)
2687 tmp2 = tmp1*tmp1*tmp1*PIov6
2688 NCLsrs= Dsrs*DE(i,k)*(QCLsr+QCLrs)/(des*tmp2)
2689 NCLsrg= Dsrg*DE(i,k)*(QCLsr+QCLrs)/(deg*tmp2)
2690 NCLsrh= Dsrh*DE(i,k)*(QCLsr+QCLrs)/(deh*tmp2)
2691 endif
2693 if (QCLgr+QCLrg>0.) then
2694 tmp1 = max(Dr,Dg)
2695 tmp2 = tmp1*tmp1*tmp1*PIov6
2696 NCLgrg= Dgrg*DE(i,k)*(QCLgr+QCLrg)/(deg*tmp2)
2697 NCLgrh= Dgrh*DE(i,k)*(QCLgr+QCLrg)/(deh*tmp2)
2698 endif
2700 !========================================================================!
2701 ! Add all source/sink terms to all predicted moments: !
2702 !========================================================================!
2704 !Diagnostic S/S terms: (to facilitate output of 3D variables for diagnostics)
2705 !e.g. SS(i,k,1)= QVDvs*idt (for depositional growth rate of snow, kg kg-1 s-1)
2707 ! Q-Source/Sink Terms:
2708 Q(i,k) = Q(i,k) -QNUvi -QVDvi -QVDvs -QVDvg -QVDvh
2709 QC(i,k)= QC(i,k) -QCLcs -QCLcg -QCLch -QFZci
2710 QR(i,k)= QR(i,k) -QCLri +QMLsr -QCLrs -QCLrg +QMLgr -QCLrh +QMLhr -QFZrh +QMLir
2711 QI(i,k)= QI(i,k) +QNUvi +QVDvi +QFZci -QCNis -QCLir -QCLis -QCLig &
2712 -QMLir -QCLih +QIMsi +QIMgi
2713 QG(i,k)= QG(i,k) +QCNsg +QVDvg +QCLcg -QCLgr-QMLgr -QCNgh -QIMgi +QCLig &
2714 +Dirg*(QCLri+QCLir) +Dgrg*(QCLrg+QCLgr) +Dsrg*(QCLrs+QCLsr)
2715 QN(i,k)= QN(i,k) +QCNis +QVDvs +QCLcs -QCNsg -QMLsr -QIMsi -QCLsr +QCLis -QCLsh &
2716 +Dsrs*(QCLrs+QCLsr)
2717 QH(i,k)= QH(i,k) +Dirh*(QCLri+QCLir) -QMLhr +QVDvh +QCLch +Dsrh*(QCLrs+QCLsr) &
2718 +QCLih +QCLsh +QFZrh +QCLrh +QCNgh +Dgrh*(QCLrg+QCLgr)
2720 ! N-Source/Sink Terms:
2721 NC(i,k)= NC(i,k) -NCLcs -NCLcg -NCLch -NFZci
2722 NR(i,k)= NR(i,k) -NCLri -NCLrs -NCLrg -NCLrh +NMLsr +NMLgr +NMLhr -NrFZrh +NMLir &
2723 +NSHhr
2724 NY(i,k)= NY(i,k) +NNUvi +NVDvi +NFZci -NCLir -NCLis -NCLig -NCLih -NMLir +NIMsi &
2725 +NIMgi -NiCNis
2726 NN(i,k)= NN(i,k) +NsCNis -NVDvs -NCNsg -NMLsr -NCLss -NCLsr -NCLsh +NCLsrs
2727 NG(i,k)= NG(i,k) +NCNsg -NCLgr -NVDvg -NMLgr +NCLirg +NCLsrg +NCLgrg -NgCNgh
2728 NH(i,k)= NH(i,k) +NhFZrh +NhCNgh -NMLhr -NVDvh +NCLirh +NCLsrh +NCLgrh
2730 T(i,k)= T(i,k) +LFP*(QCLri+QCLcs+QCLrs+QFZci-QMLsr+QCLcg+QCLrg-QMLir-QMLgr &
2731 -QMLhr+QCLch+QCLrh+QFZrh) +LSP*(QNUvi+QVDvi+QVDvs+QVDvg+QVDvh)
2733 !--- Ensure consistency between moments: (prescribe NX (and BG) in case of inconsistency)
2734 tmp1= pres(i,k)/(RGASD*T(i,k)) !updated air density
2735 !note: In Part 2a (warm-rain coalescence), air density at initial time (before
2736 ! modification of T due to ice-phase) is used; hence, DE is not recomputed here.
2738 !cloud:
2739 if (QC(i,k)>epsQ .and. NC(i,k)<epsN) then
2740 NC(i,k)= N_c_SM
2741 elseif (QC(i,k)<=epsQ) then
2742 Q(i,k) = Q(i,k) + QC(i,k)
2743 T(i,k) = T(i,k) - LCP*QC(i,k)
2744 QC(i,k) = 0.
2745 NC(i,k) = 0.
2746 endif
2748 !rain
2749 if (QR(i,k)>epsQ .and. NR(i,k)<epsN) then
2750 NR(i,k)= (No_r_SM*GR31)**(3./(4.+alpha_r))*(GR31*iGR34*tmp1*QR(i,k)* &
2751 icmr)**((1.+alpha_r)/(4.+alpha_r))
2752 elseif (QR(i,k)<=epsQ) then
2753 Q(i,k) = Q(i,k) + QR(i,k)
2754 T(i,k) = T(i,k) - LCP*QR(i,k)
2755 QR(i,k) = 0.
2756 NR(i,k) = 0.
2757 endif
2759 !ice:
2760 if (QI(i,k)>epsQ .and. NY(i,k)<epsN) then
2761 ! NY(i,k)= N_Cooper(TRPL,T(i,k))
2762 NY(i,k) = max(2.*epsN, N_Cooper(TRPL,T(i,k)) )
2763 elseif (QI(i,k)<=epsQ) then
2764 Q(i,k) = Q(i,k) + QI(i,k)
2765 T(i,k) = T(i,k) - LSP*QI(i,k)
2766 QI(i,k) = 0.
2767 NY(i,k) = 0.
2768 endif
2770 !snow:
2771 if (QN(i,k)>epsQ .and. NN(i,k)<epsN) then
2772 No_s= Nos_Thompson(TRPL,T(i,k))
2773 NN(i,k)= (No_s*GS31)**(dms*icexs2)*(GS31*iGS40*icms*tmp1*QN(i,k))**((1.+ &
2774 alpha_s)*icexs2)
2775 elseif (QN(i,k)<=epsQ) then
2776 Q(i,k) = Q(i,k) + QN(i,k)
2777 T(i,k) = T(i,k) - LSP*QN(i,k)
2778 QN(i,k) = 0.
2779 NN(i,k) = 0.
2780 endif
2782 !grpl:
2783 if (QG(i,k)>epsQ .and. NG(i,k)<epsN) then
2784 NG(i,k) = (No_g_SM*GG31)**(3./(4.+alpha_g))*(GG31*iGG34*tmp1*QG(i,k)* &
2785 icmg)**((1.+alpha_g)/(4.+alpha_g))
2786 elseif (QG(i,k)<=epsQ) then
2787 Q(i,k) = Q(i,k) + QG(i,k)
2788 T(i,k) = T(i,k) - LSP*QG(i,k)
2789 QG(i,k) = 0.
2790 NG(i,k) = 0.
2791 endif
2793 !hail:
2794 if (QH(i,k)>epsQ .and. NH(i,k)<epsN) then
2795 NH(i,k)= (No_h_SM*GH31)**(3./(4.+alpha_h))*(GH31*iGH34*tmp1*QH(i,k)* &
2796 icmh)**((1.+alpha_h)/(4.+alpha_h))
2797 elseif (QH(i,k)<=epsQ) then
2798 Q(i,k) = Q(i,k) + QH(i,k)
2799 T(i,k) = T(i,k) - LSP*QH(i,k)
2800 QH(i,k) = 0.
2801 NH(i,k) = 0.
2802 endif
2803 if (QH(i,k)>epsQ .and. NH(i,k)>epsN) then
2804 !transfer small hail to graupel:
2805 tmp1 = 1./NH(i,k)
2806 Dh = Dm_x(DE(i,k),QH(i,k),tmp1,icmh,thrd)
2807 if (Dh < Dh_min) then
2808 QG(i,k) = QG(i,k) + QH(i,k)
2809 NG(i,k) = NG(i,k) + NH(i,k)
2810 QH(i,k) = 0.
2811 NH(i,k) = 0.
2812 endif
2813 endif
2814 !===
2815 Q(i,k)= max(Q(i,k),0.)
2816 NY(i,k)= min(NY(i,k), Ni_max)
2818 !-----
2819 if ( T(i,k)<173. .or. T(i,k)>323.) then !** DEBUG **
2820 print*, '** STOPPING IN MICROPHYSICS: (Part 2, end) **' !** DEBUG **
2821 print*, '** i,k,T [K]: ',i,k,T(i,k) !** DEBUG **
2822 stop !** DEBUG **
2823 endif !** DEBUG **
2824 !=====
2826 ENDIF !if (activePoint)
2827 ENDDO
2828 ENDDO
2830 !----------------------------------------------------------------------------------!
2831 ! End of ice phase microphysics (Part 2) !
2832 !----------------------------------------------------------------------------------!
2834 if (DEBUG_ON) call check_values(Q,T,QC,QR,QI,QN,QG,QH,NC,NR,NY,NN,NG,NH,epsQ,epsN,.true.,DEBUG_abort,450)
2836 !----------------------------------------------------------------------------------!
2837 ! PART 3: Warm Microphysics Processes !
2838 ! !
2839 ! Equations for warm-rain coalescence based on Cohard and Pinty (2000a,b; QJRMS) !
2840 ! Condensation/evaportaion equations based on Kong and Yau (1997; Atmos-Ocean) !
2841 ! Equations for rain reflectivity (ZR) based on Milbrandt and Yau (2005b; JAS) !
2842 !----------------------------------------------------------------------------------!
2844 ! Part 3a - Warm-rain Coallescence:
2846 IF (warmphase_ON) THEN
2848 DO k= ktop-kdir,kbot,-kdir
2849 DO i= 1,ni
2851 RCAUTR= 0.; CCACCR= 0.; Dc= 0.; iLAMc= 0.; L = 0.
2852 RCACCR= 0.; CCSCOC= 0.; Dr= 0.; iLAMr= 0.; TAU= 0.
2853 CCAUTR= 0.; CRSCOR= 0.; SIGc= 0.; DrINIT= 0.
2854 iLAMc3= 0.; iLAMc6= 0.; iLAMr3= 0.; iLAMr6= 0.
2856 rainPresent= (QR_in(i,k)>epsQ .and. NR_in(i,k)>epsN)
2858 if (QC_in(i,k)>epsQ .and. NC_in(i,k)>epsN) then
2859 iNC = 1./NC_in(i,k)
2860 iLAMc = iLAMDA_x(DE(i,k),QC_in(i,k),iNC,icexc9,thrd)
2861 iLAMc3= iLAMc*iLAMc*iLAMc
2862 iLAMc6= iLAMc3*iLAMc3
2863 Dc = iLAMc*(GC2*iGC1)**thrd
2864 SIGc = iLAMc*( GC3*iGC1- (GC2*iGC1)*(GC2*iGC1) )**sixth
2865 L = 0.027*DE(i,k)*QC_in(i,k)*(6.25e18*SIGc*SIGc*SIGc*Dc-0.4)
2866 if (SIGc>SIGcTHRS) TAU= 3.7/(DE(i,k)*(QC_in(i,k))*(0.5e6*SIGc-7.5))
2867 endif
2869 if (rainPresent) then
2870 iNR = 1./NR_in(i,k)
2871 Dr = Dm_x(DE(i,k),QR_in(i,k),iNR,icmr,thrd)
2873 !Drop-size limiter [prevents initially large drops]
2874 if (Dr>3.e-3) then
2875 tmp1 = (Dr-3.e-3)
2876 tmp2 = (Dr/DrMAX)
2877 tmp3 = tmp2*tmp2*tmp2
2878 NR_in(i,k)= NR_in(i,k)*max((1.+2.e4*tmp1*tmp1),tmp3)
2879 iNR = 1./NR_in(i,k)
2880 Dr = Dm_x(DE(i,k),QR_in(i,k),iNR,icmr,thrd)
2881 endif
2882 iLAMr = iLAMDA_x(DE(i,k),QR_in(i,k),iNR,icexr9,thrd)
2883 iLAMr3 = iLAMr*iLAMr*iLAMr
2884 iLAMr6 = iLAMr3*iLAMr3
2885 endif
2887 ! Autoconversion:
2888 if (QC_in(i,k)>epsQ .and. SIGc>SIGcTHRS .and. autoconv_ON) then
2889 RCAUTR= min( max(L/TAU,0.), QC(i,k)*idt )
2890 DrINIT= max(83.e-6, 12.6e-4/(0.5e6*SIGc-3.5)) !initiation regime Dr
2891 DrAUT = max(DrINIT, Dr) !init. or feeding DrAUT
2892 CCAUTR= RCAUTR*DE(i,k)/(cmr*DrAUT*DrAUT*DrAUT)
2894 ! ---------------------------------------------------------------------------- !
2895 ! NOTE: The formulation for CCAUTR here (dNr/dt|initiation) does NOT follow
2896 ! eqn (18) in CP2000a, but rather it comes from the F90 code provided
2897 ! by J-P Pinty (subroutine: 'rain_c2r2.f90').
2898 ! (See notes: 2001-10-17; 2001-10-22)
2899 !
2900 ! Similarly, the condition for the activation of accretion and self-
2901 ! collection depends on whether or not autoconversion is in the feeding
2902 ! regime (see notes 2002-01-07). This is apparent in the F90 code, but
2903 ! NOT in CP2000a.
2904 ! ---------------------------------------------------------------------------- !
2906 ! cloud self-collection: (dNc/dt_autoconversion) {CP eqn(25)}
2907 CCSCOC= min(KK2*NC_in(i,k)*NC_in(i,k)*GC3*iGC1*iLAMc6, NC_in(i,k)*idt) !{CP00a eqn(25)}
2908 endif
2910 ! Accretion, rain self-collection, and collisional breakup:
2911 if (((QR_in(i,k))>1.2*max(L,0.)*iDE(i,k).or.Dr>max(5.e-6,DrINIT)).and.rainAccr_ON &
2912 .and. rainPresent) then
2914 ! Accretion: !{CP00a eqn(22)}
2915 if (QC_in(i,k)>epsQ.and.L>0.) then
2916 if (Dr.ge.100.e-6) then
2917 CCACCR = KK1*(NC_in(i,k)*NR_in(i,k))*(GC2*iGC1*iLAMc3+GR34*iGR31*iLAMr3)
2918 RCACCR = cmr*iDE(i,k)*KK1*(NC_in(i,k)*NR_in(i,k))*iLAMc3*(GC3*iGC1*iLAMc3+ &
2919 GC2*iGC1*GR34*iGR31*iLAMr3)
2920 else
2921 CCACCR = KK2*(NC_in(i,k)*NR_in(i,k))*(GC3*iGC1*iLAMc6+GR37*iGR31*iLAMr6)
2923 ! RCACCR= cmr*iDE(i,k)*KK2*(NC(i,k)*NR(i,k))*iLAMc3* &
2924 ! (GC4*iGR31*iLAMc6+GC2*iGC1*GR37*iGR31*iLAMr6)
2925 !++ The following calculation of RCACCR avoids overflow:
2926 tmp1 = cmr*iDE(i,k)
2927 tmp2 = KK2*(NC_in(i,k)*NR_in(i,k))*iLAMc3
2928 RCACCR = tmp1 * tmp2
2929 tmp1 = GC4*iGR31
2930 tmp1 = (tmp1)*iLAMc6
2931 tmp2 = GC2*iGC1
2932 tmp2 = tmp2*GR37*iGR31
2933 tmp2 = (tmp2)*iLAMr6
2934 RCACCR = RCACCR * (tmp1 + tmp2)
2935 !++
2936 endif
2937 CCACCR = min(CCACCR,(NC(i,k))*idt)
2938 RCACCR = min(RCACCR,(QC(i,k))*idt)
2939 endif
2941 !Rain self-collection:
2942 tmp1= NR_in(i,k)*NR_in(i,k)
2943 if (Dr.ge.100.e-6) then
2944 CRSCOR= KK1*tmp1*GR34*iGR31*iLAMr3 !{CP00a eqn(24)}
2945 else
2946 CRSCOR= KK2*tmp1*GR37*iGR31*iLAMr6 !{CP00a eqn(25)}
2947 endif
2948 !Raindrop breakup: !{CP00a eqn(26)}
2949 Ec= 1.
2950 ! if (Dr > 300.e-6) Ec = 2.-exp(2300.*(Dr-300.e-6))
2951 if (iLAMr > 300.e-6) Ec = 2.-exp(2300.*(iLAMr-300.e-6)) !(assumes alpha_r=0)
2953 CRSCOR= min(Ec*CRSCOR,(0.5*NR(i,k))*idt) !0.5 prevents depletion of NR
2955 endif !accretion/self-collection/breakup
2957 ! Prevent overdepletion of cloud:
2958 source= QC(i,k)
2959 sink = (RCAUTR+RCACCR)*dt
2960 if (sink>source) then
2961 ratio = source/sink
2962 RCAUTR= ratio*RCAUTR
2963 RCACCR= ratio*RCACCR
2964 CCACCR= ratio*CCACCR
2965 endif
2967 ! Apply tendencies:
2968 QC(i,k)= max(0., QC(i,k)+(-RCAUTR-RCACCR)*dt )
2969 QR(i,k)= max(0., QR(i,k)+( RCAUTR+RCACCR)*dt )
2970 NC(i,k)= max(0., NC(i,k)+(-CCACCR-CCSCOC)*dt )
2971 NR(i,k)= max(0., NR(i,k)+( CCAUTR-CRSCOR)*dt )
2973 if (QR(i,k)>epsQ .and. NR(i,k)>epsN) then
2974 iNR = 1./NR(i,k)
2975 Dr = Dm_x(DE(i,k),QR(i,k),iNR,icmr,thrd)
2976 if (Dr>3.e-3) then
2977 tmp1= (Dr-3.e-3); tmp2= tmp1*tmp1
2978 tmp3= (Dr/DrMAX); tmp4= tmp3*tmp3*tmp3
2979 NR(i,k)= NR(i,k)*(max((1.+2.e4*tmp2),tmp4))
2980 elseif (Dr<Dhh) then
2981 !Convert small raindrops to cloud:
2982 QC(i,k)= QC(i,k) + QR(i,k)
2983 NC(i,k)= NC(i,k) + NR(i,k)
2984 QR(i,k)= 0.; NR(i,k)= 0.
2985 endif
2986 else
2987 QR(i,k)= 0.; NR(i,k)= 0.
2988 endif !(Qr,Nr>eps)
2990 ! Part 3b - Condensation/Evaporation:
2992 QSW(i,k) = qsat(T(i,k),pres(i,k),0) !Flatau formulation
2993 if (QSW(i,k)>1.e-20) then
2994 ssat = Q(i,k)/QSW(i,k)-1. !supersaturation ratio
2995 else
2996 ssat = 0.
2997 endif
2998 X = Q(i,k)-QSW(i,k) !saturation excess (deficit)
2999 !adjustment for latent heating during cond/evap
3000 ! X = X/(1.+ck5*QSW(i,k)/(T(i,k)-35.86)**2) !orig (KY97)
3001 X = X / ( 1.+ ((3.1484e6-2370.*T(i,k))**2 * QSW(i,k))/( (1005.*(1.+ & !morr2mom
3002 0.887*Q(i,k))) *461.5*T(i,k)**2 ) )
3003 X = max(X, -QC(i,k)) !ensure no overdepletion of QC
3004 QC(i,k) = QC(i,k) + X
3005 Q(i,k) = Q(i,k) - X
3006 T(i,k) = T(i,k) + LCP*X
3008 if (X>0.) then
3009 !nucleation of cloud droplets:
3010 if (WZ(i,k)>0.001) then
3011 !condensation and non-negligible upward motion:
3012 !note: WZ threshold of 1 mm/s is to overflow problem in NccnFNC, which
3013 ! uses a polynomial approximation that is invalid for tiny WZ.
3014 NC(i,k) = max(NC(i,k), NccnFNC(WZ(i,k),T(i,k),pres(i,k),CCNtype))
3015 else
3016 !condensation and negible or downward vertical motion:
3017 NC(i,k) = max(NC(i,k), N_c_SM)
3018 endif
3019 else
3020 if (QC(i,k)>epsQ) then
3021 !partial evaporation of cloud droplets:
3022 NC(i,k) = max(0., NC(i,k) + X*NC(i,k)/max(QC(i,k),epsQ) ) !(dNc/dt)|evap
3023 else
3024 NC(i,k) = 0.
3025 endif
3026 endif
3027 if (QC(i,k)>epsQ .and. NC(i,k)<epsN) NC(i,k)= N_c_SM !prevents non-zero_Q & zero_N
3029 !rain evaporation:
3030 !saturation adjustment:
3031 QSW(i,k) = qsat(T(i,k),pres(i,k),0) !Flatau formulation
3032 if (Q(i,k)<QSW(i,k) .and. QR(i,k)>epsQ .and. NR(i,k)>epsN) then
3033 ssat = Q(i,k)/QSW(i,k)-1.
3034 Tc = T(i,k)-TRPL
3035 Cdiff = max(1.62e-5, (2.2157e-5 + 0.0155e-5*Tc)) *1.e5/pres(i,k)
3036 MUdyn = max(1.51e-5, (1.7153e-5 + 0.0050e-5*Tc))
3037 Ka = max(2.07e-2, (2.3971e-2 + 0.0078e-2*Tc))
3038 MUkin = MUdyn*iDE(i,k)
3039 iMUkin = 1./MUkin
3040 ScTHRD = (MUkin/Cdiff)**thrd
3041 X = QSW(i,k) - Q(i,k) !saturation exesss(deficit)
3042 !adjustment for latent cooling during evaporation:
3043 ! X = X/(1.+ck5*QSW(i,k)/(T(i,k)-35.86)**2) !orig (KY97)
3044 X = X / ( 1.+ ((3.1484e6-2370.*T(i,k))**2 * QSW(i,k))/( (1005.*(1.+ & !morr2mom
3045 0.887*Q(i,k))) *461.5*T(i,k)**2 ) )
3046 DE(i,k) = pres(i,k)/(RGASD*T(i,k)) !recompute air density (with updated T)
3047 iDE(i,k) = 1./DE(i,k)
3048 gam = sqrt(DEo*iDE(i,k))
3049 tmp1 = 1./NR(i,k)
3050 iLAMr = iLAMDA_x(DE(i,k),QR(i,k),tmp1,icexr9,thrd)
3051 LAMr = 1./iLAMr
3052 !note: The following coding of 'No_r=...' prevents overflow:
3053 !No_r = NR(i,k)*LAMr**(1.+alpha_r))*iGR31
3054 No_r = dble(NR(i,k))*dble(LAMr)**dble(1.+alpha_r)*iGR31
3055 !note: There is an error in MY05a_eq(8) for VENTx (corrected in code)
3056 VENTr = Avx*GR32*iLAMr**cexr5 + Bvx*ScTHRD*sqrt(gam*afr*iMUkin)*GR17*iLAMr**cexr6
3057 ABw = CHLC**2/(Ka*RGASV*T(i,k)**2)+1./(DE(i,k)*(QSW(i,k))*Cdiff)
3058 QREVP = min( QR(i,k), -dt*(iDE(i,k)*PI2*ssat*No_r*VENTr/ABw) )
3059 tmp1 = QR(i,k) !value of QR before update, used in NR update eqn
3060 T(i,k) = T(i,k) - LCP*QREVP
3061 Q(i,k) = Q(i,k) + QREVP
3062 QR(i,k) = QR(i,k) - QREVP
3063 NR(i,k) = max(0., NR(i,k) - QREVP*NR(i,k)/tmp1)
3064 !Protect against negative values due to overdepletion:
3065 if (QR(i,k)<epsQ .or. NR(i,k)<epsN) then
3066 Q(i,k) = Q(i,k) + QR(i,k)
3067 T(i,k) = T(i,k) - QR(i,k)*LCP
3068 QR(i,k) = 0.
3069 NR(i,k) = 0.
3070 endif
3071 endif
3073 !homogeneous freezing of cloud:
3074 Tc = T(i,k) - TRPL
3075 if (QC(i,k)>epsQ .and. Tc<-30. .and. icephase_ON) then
3077 !-detailed:
3078 ! Tc = T(i,k) - TRPL
3079 ! JJ = (10.**max(-20.,(-606.3952-52.6611*Tc-1.7439*Tc**2-0.0265*Tc**3- &
3080 ! 1.536e-4*Tc**4)))
3081 ! tmp1 = 1.e6*(DE(i,k)*(QC(i,k)/NC(i,k))*icmr) !i.e. Dc[cm]**3
3082 ! FRAC = 1.-exp(-JJ*PIov6*tmp1*dt)
3083 ! if (Tc>-30.) FRAC= 0.
3084 ! if (Tc<-50.) FRAC= 1.
3085 !-simplified:
3086 if (Tc<-35.) then
3087 FRAC = 1.
3088 else
3089 FRAC = 0.
3090 endif
3091 !=
3092 QFZci = FRAC*QC(i,k)
3093 NFZci = FRAC*NC(i,k)
3094 QC(i,k) = QC(i,k) - QFZci
3095 NC(i,k) = NC(i,k) - NFZci
3096 QI(i,k) = QI(i,k) + QFZci
3097 NY(i,k) = NY(i,k) + NFZci
3098 T(i,k) = T(i,k) + LFP*QFZci
3100 if (QC(i,k)>epsQ .and. NC(i,k)<epsN) then
3101 NC(i,k)= N_c_SM
3102 elseif (QC(i,k)<=epsQ) then
3103 Q(i,k) = Q(i,k) + QC(i,k)
3104 T(i,k) = T(i,k) - LCP*QC(i,k)
3105 QC(i,k) = 0.
3106 NC(i,k) = 0.
3107 endif
3109 if (QI(i,k)>epsQ .and. NY(i,k)<epsN) then
3110 ! NY(i,k)= N_Cooper(TRPL,T(i,k))
3111 NY(i,k) = max(2.*epsN, N_Cooper(TRPL,T(i,k)) )
3112 elseif (QI(i,k)<=epsQ) then
3113 Q(i,k) = Q(i,k) + QI(i,k)
3114 T(i,k) = T(i,k) - LSP*QI(i,k)
3115 QI(i,k) = 0.
3116 NY(i,k) = 0.
3117 endif
3119 endif
3121 !==
3123 ENDDO
3124 ENDDO !cond/evap [k-loop]
3126 ENDIF !if warmphase_ON
3128 !----------------------------------------------------------------------------------!
3129 ! End of warm-phase microphysics (Part 3) !
3130 !----------------------------------------------------------------------------------!
3132 if (DEBUG_ON) call check_values(Q,T,QC,QR,QI,QN,QG,QH,NC,NR,NY,NN,NG,NH,epsQ,epsN,.true.,DEBUG_abort,500)
3134 !----------------------------------------------------------------------------------!
3135 ! PART 4: Sedimentation !
3136 !----------------------------------------------------------------------------------!
3138 IF (sedi_ON) THEN
3140 fluxM_r= 0.; fluxM_i= 0.; fluxM_s= 0.; fluxM_g= 0.; fluxM_h= 0.
3141 RT_rn1 = 0.; RT_rn2 = 0.; RT_fr1 = 0.; RT_fr2 = 0.; RT_sn1 = 0.
3142 RT_sn2 = 0.; RT_sn3 = 0.; RT_pe1 = 0.; RT_pe2 = 0.; RT_peL = 0.
3144 !---- For sedimentation on all levels:
3145 call sedi_wrapper_2(QR,NR,1,epsQ,epsQr_sedi,epsN,dmr,ni,VrMax,DrMax,dt,fluxM_r,kdir, &
3146 kbot,ktop_sedi,GRAV,zheight,nk,DE,iDE,iDP,DZ,iDZ,gamfact,kount, &
3147 afr,bfr,cmr,ckQr1,ckQr2,icexr9)
3148 if (DEBUG_ON) call check_values(Q,T,QC,QR,QI,QN,QG,QH,NC,NR,NY,NN,NG,NH,epsQ,epsN,.false.,DEBUG_abort,610)
3150 call sedi_wrapper_2(QI,NY,2,epsQ,epsQi_sedi,epsN,dmi,ni,ViMax,DiMax,dt,fluxM_i,kdir, &
3151 kbot,ktop_sedi,GRAV,zheight,nk,DE,iDE,iDP,DZ,iDZ,gamfact,kount, &
3152 afi,bfi,cmi,ckQi1,ckQi2,ckQi4)
3154 if (DEBUG_ON) call check_values(Q,T,QC,QR,QI,QN,QG,QH,NC,NR,NY,NN,NG,NH,epsQ,epsN,.false.,DEBUG_abort,620)
3156 call sedi_wrapper_2(QN,NN,3,epsQ,epsQs_sedi,epsN,dms,ni,VsMax,DsMax,dt,fluxM_s,kdir, &
3157 kbot,ktop_sedi,GRAV,zheight,nk,DE,iDE,iDP,DZ,iDZ,gamfact,kount, &
3158 afs,bfs,cms,ckQs1,ckQs2,iGS20)
3160 if (DEBUG_ON) call check_values(Q,T,QC,QR,QI,QN,QG,QH,NC,NR,NY,NN,NG,NH,epsQ,epsN,.false.,DEBUG_abort,630)
3162 call sedi_wrapper_2(QG,NG,4,epsQ,epsQg_sedi,epsN,dmg,ni,VgMax,DgMax,dt,fluxM_g,kdir, &
3163 kbot,ktop_sedi,GRAV,zheight,nk,DE,iDE,iDP,DZ,iDZ,gamfact,kount, &
3164 afg,bfg,cmg,ckQg1,ckQg2,ckQg4)
3166 if (DEBUG_ON) call check_values(Q,T,QC,QR,QI,QN,QG,QH,NC,NR,NY,NN,NG,NH,epsQ,epsN,.false.,DEBUG_abort,640)
3168 call sedi_wrapper_2(QH,NH,5,epsQ,epsQh_sedi,epsN,dmh,ni,VhMax,DhMax,dt,fluxM_h,kdir, &
3169 kbot,ktop_sedi,GRAV,zheight,nk,DE,iDE,iDP,DZ,iDZ,gamfact,kount, &
3170 afh,bfh,cmh,ckQh1,ckQh2,ckQh4)
3171 !====
3174 if (DEBUG_ON) call check_values(Q,T,QC,QR,QI,QN,QG,QH,NC,NR,NY,NN,NG,NH,epsQ,epsN,.false.,DEBUG_abort,650)
3176 !--- Impose constraints on size distribution parameters ---!
3177 do k= ktop,kbot,-kdir
3178 do i= 1,ni
3179 !snow:
3180 if (QN(i,k)>epsQ .and. NN(i,k)>epsN) then
3181 !Impose No_s max for snow: (assumes alpha_s=0.)
3182 tmp1 = 1./NN(i,k)
3183 iLAMs = max( iLAMmin2, iLAMDA_x(DE(i,k),QN(i,k),tmp1,iGS20,idms) )
3184 tmp1 = min(NN(i,k)/iLAMs,No_s_max) !min. No_s
3185 NN(i,k)= tmp1**(dms/(1.+dms))*(iGS20*DE(i,k)*QN(i,k))**(1./(1.+dms)) !impose Nos_max
3186 !Impose LAMDAs_min (by increasing LAMDAs if it is below LAMDAs_min2 [2xLAMDAs_min])
3187 tmp1 = 1./NN(i,k)
3188 iLAMs = max( iLAMmin2, iLAMDA_x(DE(i,k),QN(i,k),tmp1,iGS20,idms) )
3189 tmp2 = 1./iLAMs !LAMs before adjustment
3190 !adjust value of lamdas_min to be applied:
3191 ! This adjusts for multiple corrections (each time step). The factor 0.6 was obtained by
3192 ! trial-and-error to ultimately give reasonable LAMs profiles, smooth and with min LAMs~lamdas_min
3193 tmp4 = 0.6*lamdas_min
3194 tmp5 = 2.*tmp4
3195 tmp3 = tmp2 + tmp4*(max(0.,tmp5-tmp2)/tmp5)**2. !LAMs after adjustment
3196 tmp3 = max(tmp3, lamdas_min) !final correction
3197 NN(i,k)= NN(i,k)*(tmp3*iLAMs)**dms !re-compute NN after LAMs adjustment
3198 endif
3199 enddo !i-loop
3200 enddo !k-loop
3201 !===
3203 !Compute melted (liquid-equivalent) volume fluxes [m3 (liquid) m-2 (sfc area) s-1]:
3204 ! (note: For other precipitation schemes in RPN-CMC physics, this is computed in 'vkuocon6.ftn')
3205 RT_rn1 = fluxM_r *idew
3206 RT_sn1 = fluxM_i *idew
3207 RT_sn2 = fluxM_s *idew
3208 RT_sn3 = fluxM_g *idew
3209 RT_pe1 = fluxM_h *idew
3211 if (DEBUG_ON) call check_values(Q,T,QC,QR,QI,QN,QG,QH,NC,NR,NY,NN,NG,NH,epsQ,epsN,.false.,DEBUG_abort,700)
3213 !----
3214 !Compute sum of solid (unmelted) volume fluxes [m3 (bulk hydrometeor) m-2 (sfc area) s-1]:
3215 !(this is the precipitation rate for UNmelted total snow [i+s+g])
3217 ! Note: In 'calcdiagn.ftn', the total solid precipitation (excluding hail) SN is computed
3218 ! from the sum of the liq-eq.vol fluxes, tss_sn1 + tss_sn2 + tss_sn3. With the
3219 ! accumulation of SND (in 'calcdiag.ftn'), the solid-to-liquid ratio for the total
3220 ! accumulated "snow" (i+s+g) can be compute as SND/SN. Likewise, the instantaneous
3221 ! solid-to-liquid ratio of solid precipitation is computed (in 'calcdiag.ftn') as
3222 ! RS2L = RSND/RSN.
3224 do i= 1,ni
3226 fluxV_i= fluxM_i(i)*idei
3227 fluxV_g= fluxM_g(i)*ideg
3228 !Compute unmelted volume flux for snow:
3229 ! note: This is based on the strict calculation of the volume flux, where vol=(pi/6)D^3,
3230 ! and remains in the integral calculation Fv = int[V(D)*vol(D)*N(D)*dD].
3231 ! For a constant density (ice and graupel), vol(D) = m(D)/dex, dex comes out of
3232 ! integral and Fv_x=Fm_x/dex
3233 ! Optimized for alpha_s = 0.
3234 if (QN(i,nk)>epsQ .and. NN(i,nk)>epsN .and. fluxM_s(i)>0.) then
3235 tmp2 = 1./NN(i,nk)
3236 tmp1 = 1./iLAMDA_x(DE(i,nk),QN(i,nk),tmp2,iGS20,idms) !LAMDA_s
3237 fluxV_s = fluxM_s(i)*rfact_FvFm*tmp1**(dms-3.)
3238 else
3239 fluxV_s = 0.
3240 endif
3242 !total solid unmelted volume flux, before accounting for partial melting:
3243 tmp1= fluxV_i + fluxV_g + fluxV_s
3245 !liquid-fraction of partially-melted solid precipitation:
3246 ! The physical premise is that if QR>0, QN+QI+QG>0, and T>0, then QR
3247 ! originates from melting and can be used to estimate the liquid portion
3248 ! of the partially-melted solid hydrometeor.
3249 tmp2= QR(i,nk) + QI(i,nk) + QN(i,nk) + QG(i,nk)
3250 if (T(i,nk)>TRPL .and. tmp2>epsQ) then
3251 fracLiq= QR(i,nk)/tmp2
3252 else
3253 fracLiq= 0.
3254 endif
3256 !Tend total volume flux towards the liquid-equivalent as the liquid-fraction increases to 1:
3257 tmp3= RT_sn1(i) + RT_sn2(i) + RT_sn3(i) !total liquid-equivalent volume flux (RSN, Fv_sol)
3258 RT_snd(i)= (1.-fracLiq)*tmp1 + fracLiq*tmp3 !total volume flux with partial melting (RSND, Fvsle_sol)
3259 !Note: Calculation of instantaneous solid-to-liquid ratio [RS2L = RSND/RSN]
3260 ! is based on the above quantities and is done on 'calcdiag.ftn'.
3262 enddo !i-loop
3263 !====
3265 !++++
3266 ! Diagnose surface precipitation types:
3267 !
3268 ! The following involves diagnostic conditions to determine surface precipitation rates
3269 ! for various precipitation elements defined in Canadian Meteorological Operational Internship
3270 ! Program module 3.1 (plus one for large hail) based on the sedimentation rates of the five
3271 ! sedimenting hydrometeor categories.
3272 !
3273 ! With the diagnostics shut off (precipDiag_ON=.false.), 5 rates are just the 5 category
3274 ! rates, with the other 6 rates just 0. The model output variables will have:
3275 !
3276 ! total liquid = RT_rn1 [RAIN]
3277 ! total solid = RT_sn1 [ICE] + RT_sn2 [SNOW] + RT_sn3 [GRAUPEL] + RT_pe1 [HAIL]
3278 !
3279 ! With the diagnostics on, the 5 sedimentation rates are partitioned into 9 rates,
3280 ! with the following model output variable:
3281 !
3282 ! total liquid = RT_rn1 [liquid rain] + RT_rn2 [liquid drizzle]
3283 ! total solid = RT_fr1 [freezing rain] + RT_fr2 [freezing drizzle] + RT_sn1 [ice crystals] +
3284 ! RT_sn2 [snow] + RT_sn3 [graupel] + RT_pe1 [ice pellets] + RT_pe2 [hail]
3285 !
3286 ! NOTE: - The above total liquid/solid rates are computed in 'calcdiag.ftn' (as R2/R4).
3287 ! - RT_peL [large hail] is a sub-set of RT_pe2 [hail] and is thus not added to the total
3288 ! solid precipitation.
3290 ! call tmg_start0(98,'mmCalcDIAG')
3292 IF (precipDiag_ON) THEN
3293 DO i= 1,ni
3295 DE(i,kbot)= pres(i,kbot)/(RGASD*T(i,kbot))
3297 !rain vs. drizzle:
3298 N_r= NR(i,nk)
3299 if (QR(i,kbot)>epsQ .and. N_r>epsN) then
3300 Dm_r(i,kbot)= (DE(i,nk)*icmr*QR(i,kbot)/N_r)**thrd
3301 if (Dm_r(i,kbot)>Dr_large) then !Dr_large is rain/drizzle size threshold
3302 RT_rn2(i)= RT_rn1(i); RT_rn1(i)= 0.
3303 endif
3304 endif
3306 !liquid vs. freezing rain or drizzle:
3307 if (T(i,nk)<TRPL) then
3308 RT_fr1(i)= RT_rn1(i); RT_rn1(i)= 0.
3309 RT_fr2(i)= RT_rn2(i); RT_rn2(i)= 0.
3310 endif
3312 !ice pellets vs. hail:
3313 if (T(i,nk)>(TRPL+5.0)) then
3314 !note: The condition (T_sfc<5C) for ice pellets is a simply proxy for the presence
3315 ! of a warm layer aloft, though which falling snow or graupel will melt to rain,
3316 ! over a sub-freezinglayer, where the rain will freeze into the 'hail' category
3317 RT_pe2(i)= RT_pe1(i); RT_pe1(i)= 0.
3318 endif
3320 !large hail:
3321 if (QH(i,kbot)>epsQ) then
3322 N_h = NH(i,kbot)
3323 tmp1 = 1./N_h
3324 Dm_h(i,kbot)= Dm_x(DE(i,kbot),QH(i,kbot),tmp1,icmh,thrd)
3325 if (DM_h(i,kbot)>Dh_large) RT_peL(i)= RT_pe2(i)
3326 !note: large hail (RT_peL) is a subset of the total hail (RT_pe2)
3327 endif
3329 ENDDO
3330 ENDIF !if (precipDiag_ON)
3331 !
3332 !++++
3334 ENDIF ! if (sedi_ON)
3336 where (Q<0.) Q= 0.
3338 !-----------------------------------------------------------------------------------!
3339 ! End of sedimentation calculations (Part 4) !
3340 !-----------------------------------------------------------------------------------!
3342 if (DEBUG_ON) call check_values(Q,T,QC,QR,QI,QN,QG,QH,NC,NR,NY,NN,NG,NH,epsQ,epsN,.false.,DEBUG_abort,800)
3344 !===================================================================================!
3345 ! End of microphysics scheme !
3346 !===================================================================================!
3348 !-----------------------------------------------------------------------------------!
3349 ! Calculation of diagnostic variables: !
3351 !Compute effective radii for cloud and ice (to be passed to radiation scheme):
3352 ! - based of definition, r_eff = M_r(3)/M_r(2) = 0.5*M_D(3)/M_D(2),
3353 ! where the pth moment w.r.t. diameter, M_D(p), is given by MY2005a, eqn (2)
3354 do k = kbot,ktop,kdir
3355 do i = 1,ni
3357 !cloud:
3358 if (QC(i,k)>epsQ .and. NC(i,k)>epsN) then
3359 !hardcoded for alpha_c = 1. and mu_c = 3.
3360 iNC = 1./NC(i,k)
3361 iLAMc = iLAMDA_x(DE(i,k),QC(i,k),iNC,icexc9,thrd)
3362 ! reff_c(i,k) = 0.664639*iLAMc
3363 else
3364 ! reff_c(i,k) = 0.
3365 endif
3367 !ice:
3368 if (QI(i,k)>epsQ .and. NY(i,k)>epsN) then
3369 !hardcoded for alpha_i = 0. and mu_i = 1.
3370 iNY = 1./NY(i,k)
3371 iLAMi = max( iLAMmin2, iLAMDA_x(DE(i,k),QI(i,k),iNY,icexi9,thrd) )
3372 ! reff_i(i,k) = 1.5*iLAMi
3373 else
3374 ! reff_i(i,k) = 0.
3375 endif
3377 enddo
3378 enddo
3381 IF (calcDiag) THEN
3383 !For reflectivity calculations:
3384 ZEC= minZET
3385 cxr= icmr*icmr !for rain
3386 cxi= 1./fdielec*icmr*icmr !for all frozen categories
3387 Gzr= (6.+alpha_r)*(5.+alpha_r)*(4.+alpha_r)/((3.+alpha_r)*(2.+alpha_r)*(1.+alpha_r))
3388 Gzi= (6.+alpha_i)*(5.+alpha_i)*(4.+alpha_i)/((3.+alpha_i)*(2.+alpha_i)*(1.+alpha_i))
3389 if (snowSpherical) then !dms=3
3390 Gzs= (6.+alpha_s)*(5.+alpha_s)*(4.+alpha_s)/((3.+alpha_s)*(2.+alpha_s)* &
3391 (1.+alpha_s))
3392 else !dms=2
3393 Gzs= (4.+alpha_s)*(3.+alpha_s)/((2.+alpha_s)*(1.+alpha_s))
3394 endif
3395 Gzg= (6.+alpha_g)*(5.+alpha_g)*(4.+alpha_g)/((3.+alpha_g)*(2.+alpha_g)*(1.+alpha_g))
3396 Gzh= (6.+alpha_h)*(5.+alpha_h)*(4.+alpha_h)/((3.+alpha_h)*(2.+alpha_h)*(1.+alpha_h))
3398 do k= ktop,kbot,-kdir
3399 do i= 1,ni
3400 DE(i,k)= pres(i,k)/(RGASD*T(i,k))
3401 tmp9= DE(i,k)*DE(i,k)
3403 N_c= NC(i,k)
3404 N_r= NR(i,k)
3405 N_i= NY(i,k)
3406 N_s= NN(i,k)
3407 N_g= NG(i,k)
3408 N_h= NH(i,k)
3410 !Total equivalent reflectivity: (units of [dBZ])
3411 tmp1= 0.; tmp2= 0.; tmp3= 0.; tmp4= 0.; tmp5= 0.
3412 if (QR(i,k)>epsQ .and. N_r>epsN) tmp1 = cxr*Gzr*tmp9*QR(i,k)*QR(i,k)/N_r
3413 if (QI(i,k)>epsQ .and. N_i>epsN) tmp2 = cxi*Gzi*tmp9*QI(i,k)*QI(i,k)/N_i
3414 if (QN(i,k)>epsQ .and. N_s>epsN) tmp3 = cxi*Gzs*tmp9*QN(i,k)*QN(i,k)/N_s
3415 if (QG(i,k)>epsQ .and. N_g>epsN) tmp4 = cxi*Gzg*tmp9*QG(i,k)*QG(i,k)/N_g
3416 if (QH(i,k)>epsQ .and. N_h>epsN) tmp5 = cxi*Gzh*tmp9*QH(i,k)*QH(i,k)/N_h
3417 !Modifiy dielectric constant for melting ice-phase categories:
3418 if ( T(i,k)>TRPL) then
3419 tmp2= tmp2*fdielec
3420 tmp3= tmp3*fdielec
3421 tmp4= tmp4*fdielec
3422 tmp5= tmp5*fdielec
3423 endif
3424 ZET(i,k) = tmp1 + tmp2 + tmp3 + tmp4 + tmp5 != Zr+Zi+Zs+Zg+Zh
3425 if (ZET(i,k)>0.) then
3426 ZET(i,k)= 10.*log10((ZET(i,k)*Zfact)) !convert to dBZ
3427 else
3428 ZET(i,k)= minZET
3429 endif
3430 ZET(i,k)= max(ZET(i,k),minZET)
3431 ZEC(i)= max(ZEC(i),ZET(i,k)) !composite (max in column) of ZET
3433 !Mean-mass diameters: (units of [m])
3434 Dm_c(i,k)= 0.; Dm_r(i,k)= 0.; Dm_i(i,k)= 0.
3435 Dm_s(i,k)= 0.; Dm_g(i,k)= 0.; Dm_h(i,k)= 0.
3436 if (QC(i,k)>epsQ.and.N_c>epsN) then
3437 tmp1 = 1./N_c
3438 Dm_c(i,k) = Dm_x(DE(i,k),QC(i,k),tmp1,icmr,thrd)
3439 endif
3440 if (QR(i,k)>epsQ.and.N_r>epsN) then
3441 tmp1 = 1./N_r
3442 Dm_r(i,k) = Dm_x(DE(i,k),QR(i,k),tmp1,icmr,thrd)
3443 endif
3444 if (QI(i,k)>epsQ.and.N_i>epsN) then
3445 tmp1 = 1./N_i
3446 Dm_i(i,k) = Dm_x(DE(i,k),QI(i,k),tmp1,icmi,thrd)
3447 endif
3448 if (QN(i,k)>epsQ.and.N_s>epsN) then
3449 tmp1 = 1./N_s
3450 Dm_s(i,k) = Dm_x(DE(i,k),QN(i,k),tmp1,icms,idms)
3451 endif
3452 if (QG(i,k)>epsQ.and.N_g>epsN) then
3453 tmp1 = 1./N_g
3454 Dm_g(i,k) = Dm_x(DE(i,k),QG(i,k),tmp1,icmg,thrd)
3455 endif
3456 if (QH(i,k)>epsQ.and.N_h>epsN) then
3457 tmp1 = 1./N_h
3458 Dm_h(i,k) = Dm_x(DE(i,k),QH(i,k),tmp1,icmh,thrd)
3459 endif
3461 enddo !i-loop
3462 enddo !k-loop
3464 ENDIF
3466 !-- Convert N from #/m3 to #/kg:
3467 iDE = (RGASD*T)/pres
3468 NC = NC*iDE
3469 NR = NR*iDE
3470 NY = NY*iDE
3471 NN = NN*iDE
3472 NG = NG*iDE
3473 NH = NH*iDE
3475 !-- Ensure than no negative final values:
3476 ! QC = max(QC, 0.); NC = max(NC, 0.)
3477 ! QR = max(QR, 0.); NR = max(NR, 0.)
3478 ! QI = max(QI, 0.); NY = max(NY, 0.)
3479 ! QN = max(QN, 0.); NN = max(NN, 0.)
3480 ! QG = max(QG, 0.); NG = max(NG, 0.)
3481 ! QH = max(QH, 0.); NH = max(NH, 0.)
3483 if (DEBUG_ON) call check_values(Q,T,QC,QR,QI,QN,QG,QH,NC,NR,NY,NN,NG,NH,epsQ,epsN,.false.,DEBUG_abort,900)
3485 END SUBROUTINE mp_milbrandt2mom_main
3487 !___________________________________________________________________________________!
3489 real function des_OF_Ds(Ds_local,desMax_local,eds_local,fds_local)
3490 !Computes density of equivalent-volume snow particle based on (pi/6*des)*Ds^3 = cms*Ds^dms
3491 real :: Ds_local,desMax_local,eds_local,fds_local
3492 ! des_OF_Ds= min(desMax_local, eds_local*Ds_local**fds_local)
3493 des_OF_Ds= min(desMax_local, eds_local*exp(fds_local*log(Ds_local))) !IBM optimization
3494 end function des_OF_Ds
3497 real function Dm_x(DE_local,QX_local,iNX_local,icmx_local,idmx_local)
3498 !Computes mean-mass diameter
3499 real :: DE_local,QX_local,iNX_local,icmx_local,idmx_local
3500 !Dm_x = (DE_local*QX_local*iNX_local*icmx_local)**idmx_local
3501 Dm_x = exp(idmx_local*log(DE_local*QX_local*iNX_local*icmx_local)) !IBM optimization
3502 end function Dm_x
3505 real function iLAMDA_x(DE_local,QX_local,iNX_local,icex_local,idmx_local)
3506 !Computes 1/LAMDA ("slope" parameter):
3507 real :: DE_local,QX_local,iNX_local,icex_local,idmx_local
3508 !iLAMDA_x = (DE_local*QX_local*iNX_local*icex_local)**idmx_local
3509 iLAMDA_x = exp(idmx_local*log(DE_local*QX_local*iNX_local*icex_local)) !IBM optimization
3510 end function
3513 real function N_Cooper(TRPL_local,T_local)
3514 !Computes total number concentration of ice as a function of temperature
3515 !according to parameterization of Cooper (1986):
3516 real :: TRPL_local,T_local
3517 N_Cooper= 5.*exp(0.304*(TRPL_local-max(233.,T_local)))
3518 end function N_Cooper
3520 real function Nos_Thompson(TRPL_local,T_local)
3521 !Computes the snow intercept, No_s, as a function of temperature
3522 !according to the parameterization of Thompson et al. (2004):
3523 real :: TRPL_local,T_local
3524 Nos_Thompson= min(2.e+8, 2.e+6*exp(-0.12*min(-0.001,T_local-TRPL_local)))
3525 end function Nos_Thompson
3527 !===================================================================================================!
3529 END MODULE my2_mod
3531 !________________________________________________________________________________________!
3533 MODULE module_mp_milbrandt2mom
3535 use module_wrf_error
3536 use my2_mod, ONLY: mp_milbrandt2mom_main
3538 implicit none
3540 ! To be done later. Currently, parameters are initialized in the main routine
3541 ! (at every time step).
3543 CONTAINS
3545 !----------------------------------------------------------------------------------------!
3546 SUBROUTINE milbrandt2mom_init
3548 ! To be done later. Currently, parameters are initialized in the main routine (at every time step).
3550 END SUBROUTINE milbrandt2mom_init
3552 !----------------------------------------------------------------------------------------!
3554 !+---------------------------------------------------------------------+
3555 ! This is a wrapper routine designed to transfer values from 3D to 2D. !
3556 !+---------------------------------------------------------------------+
3559 SUBROUTINE mp_milbrandt2mom_driver(qv, qc, qr, qi, qs, qg, qh, nc, nr, ni, ns, ng, &
3560 nh, th, pii, p, w, dz, dt_in, itimestep, p8w, &
3561 RAINNC, RAINNCV, SNOWNC, SNOWNCV, GRPLNC, GRPLNCV, &
3562 HAILNC, HAILNCV, SR, Zet, &
3563 ids,ide, jds,jde, kds,kde, & ! domain dims
3564 ims,ime, jms,jme, kms,kme, & ! memory dims
3565 its,ite, jts,jte, kts,kte) ! tile dims
3567 implicit none
3569 !Subroutine arguments:
3570 integer, intent(in):: ids,ide, jds,jde, kds,kde, &
3571 ims,ime, jms,jme, kms,kme, &
3572 its,ite, jts,jte, kts,kte
3573 real, dimension(ims:ime, kms:kme, jms:jme), intent(inout):: &
3574 qv,qc,qr,qi,qs,qg,qh,nc,nr,ni,ns,ng,nh,th,Zet
3575 real, dimension(ims:ime, kms:kme, jms:jme), intent(in):: &
3576 pii,p,w,dz
3577 real, dimension(ims:ime, kms:kme, jms:jme), intent(in):: p8w
3578 real, dimension(ims:ime, jms:jme), intent(inout):: &
3579 RAINNC,RAINNCV,SNOWNC,SNOWNCV,GRPLNC,GRPLNCV,HAILNC,HAILNCV, &
3580 SR
3581 real, intent(in):: dt_in
3582 integer, intent(in):: itimestep !, ccntype
3584 !Local variables:
3585 real, dimension(1:ite-its+1,1:kte-kts+1) :: t2d,p2d,sigma2d
3587 !tentatively local; to be passed out as output variables later
3588 real, dimension(1:ite-its+1,1:kte-kts+1) :: Dm_c,Dm_r,Dm_i,Dm_s,Dm_g,Dm_h
3589 real, dimension(1:ite-its+1,1:kte-kts+1,1:20) :: SS
3590 !note: 20 is temporarily hardcoded; eventually should allocate size with parameter
3591 !integer, parameter :: numSS = 20 ! number of diagnostic arrays [SS(i,k,:)]
3593 real, dimension(1:ite-its+1) :: rt_rn1,rt_rn2,rt_fr1,rt_fr2,rt_sn1,rt_sn2,rt_sn3, &
3594 rt_pe1,rt_pe2,rt_peL,rt_snd,ZEC,p_sfc
3595 real, dimension(ims:ime, kms:kme, jms:jme) :: t3d
3597 real :: dt,ms2mmstp
3598 integer :: i,j,k,i2d,k2d,i2d_max,k2d_max
3599 integer :: i_start, j_start, i_end, j_end, CCNtype
3600 logical :: precipDiag_ON,sedi_ON,warmphase_ON,autoconv_ON,icephase_ON,snow_ON
3602 real, parameter :: ms2mmh = 3.6e+6 !conversion factor for precipitation rates
3603 logical, parameter :: nk_BOTTOM = .false. !.F. for k=1 at bottom level (WRF)
3605 real, parameter :: statfreq = 300. !frequency (seconds) to output block stats
3607 !+---+
3608 i2d_max = ite-its+1
3609 k2d_max = kte-kts+1
3610 dt = dt_in
3611 ms2mmstp = 1.e+3*dt !conversion factor: m/2 to mm/step
3613 !--- temporary initialization (until variables are put as namelist options:
3614 ! CCNtype = 1. !maritime --> N_c = 80 cm-3 for dblMom_c = .F.
3615 CCNtype = 2. !continental --> N_c = 200 cm-3 for dblMom_c = .F.
3617 precipDiag_ON = .true.
3618 sedi_ON = .true.
3619 warmphase_ON = .true.
3620 autoconv_ON = .true.
3621 icephase_ON = .true.
3622 snow_ON = .true.
3623 !---
3625 RAINNCV(its:ite,jts:jte) = 0.
3626 SNOWNCV(its:ite,jts:jte) = 0.
3627 GRPLNCV(its:ite,jts:jte) = 0.
3628 HAILNCV(its:ite,jts:jte) = 0.
3629 SR(its:ite,jts:jte) = 0.
3631 !--run-time stats:
3632 ! if (mod(itimestep*dt,statfreq)==0.) then
3633 ! t3d = th*pii
3634 ! print*, 'Bloc Stats -- BEFORE micro call (my-2.25.1_b11)'
3635 ! write(6,'(a8,7e15.5)') 'Max qx: ',maxval(qc(its:ite,kts:kte,jts:jte)),maxval(qr(its:ite,kts:kte,jts:jte)),maxval(qi(its:ite,kts:kte,jts:jte)), &
3636 ! maxval(qs(its:ite,kts:kte,jts:jte)),maxval(qg(its:ite,kts:kte,jts:jte)),maxval(qh(its:ite,kts:kte,jts:jte))
3637 ! write(6,'(a8,6e15.5)') 'Max Nx: ',maxval(nc(its:ite,kts:kte,jts:jte)),maxval(nr(its:ite,kts:kte,jts:jte)),maxval(ni(its:ite,kts:kte,jts:jte)), &
3638 ! maxval(ns(its:ite,kts:kte,jts:jte)),maxval(ng(its:ite,kts:kte,jts:jte)),maxval(nh(its:ite,kts:kte,jts:jte))
3639 ! write(6,'(a8,7e15.5)') 'Min qx: ',minval(qc(its:ite,kts:kte,jts:jte)),minval(qr(its:ite,kts:kte,jts:jte)),minval(qi(its:ite,kts:kte,jts:jte)), &
3640 ! minval(qs(its:ite,kts:kte,jts:jte)),minval(qg(its:ite,kts:kte,jts:jte)),minval(qh(its:ite,kts:kte,jts:jte))
3641 ! write(6,'(a8,6e15.5)') 'Min Nx: ',minval(nc(its:ite,kts:kte,jts:jte)),minval(nr(its:ite,kts:kte,jts:jte)),minval(ni(its:ite,kts:kte,jts:jte)), &
3642 ! minval(ns(its:ite,kts:kte,jts:jte)),minval(ng(its:ite,kts:kte,jts:jte)),minval(nh(its:ite,kts:kte,jts:jte))
3643 ! write(6,'(a8,6e15.5)') 'W,T,qv: ',minval( w(its:ite,kts:kte,jts:jte)),maxval( w(its:ite,kts:kte,jts:jte)),minval(t3d(its:ite,kts:kte,jts:jte)),&
3644 ! maxval(t3d(its:ite,kts:kte,jts:jte)),minval(qv(its:ite,kts:kte,jts:jte)),maxval(qv(its:ite,kts:kte,jts:jte))
3645 ! endif
3647 do j = jts, jte
3649 t2d(:,:) = th(its:ite,kts:kte,j)*pii(its:ite,kts:kte,j)
3650 p2d(:,:) = p(its:ite,kts:kte,j)
3651 ! p_sfc(:) = p2d(:,k2d_max)
3652 p_sfc(:) = p8w(its:ite,kms, j)
3654 do i = its, ite
3655 i2d = i-its+1
3656 sigma2d(i2d,:) = p2d(i2d,:)/p_sfc(i2d)
3657 enddo
3658 call mp_milbrandt2mom_main(w(its:ite,kts:kte,j),t2d,qv(its:ite,kts:kte,j), &
3659 qc(its:ite,kts:kte,j),qr(its:ite,kts:kte,j),qi(its:ite,kts:kte,j), &
3660 qs(its:ite,kts:kte,j),qg(its:ite,kts:kte,j),qh(its:ite,kts:kte,j), &
3661 nc(its:ite,kts:kte,j),nr(its:ite,kts:kte,j),ni(its:ite,kts:kte,j), &
3662 ns(its:ite,kts:kte,j),ng(its:ite,kts:kte,j),nh(its:ite,kts:kte,j), &
3663 p_sfc,sigma2d,rt_rn1,rt_rn2,rt_fr1,rt_fr2,rt_sn1,rt_sn2,rt_sn3, &
3664 rt_pe1,rt_pe2,rt_peL,rt_snd,dt,i2d_max,k2d_max,j,itimestep,CCNtype, &
3665 precipDiag_ON,sedi_ON,warmphase_ON,autoconv_ON,icephase_ON,snow_ON, &
3666 Dm_c,Dm_r,Dm_i,Dm_s,Dm_g,Dm_h,Zet(its:ite,kts:kte,j),ZEC,SS,nk_BOTTOM)
3668 th(its:ite,kts:kte,j) = t2d(:,:)/pii(its:ite,kts:kte,j)
3670 !Convert individual precipitation rates (in m/s) to WRF precipitation fields (mm/step):
3671 RAINNCV(its:ite,j) = ( rt_rn1(:)+rt_rn2(:)+rt_fr1(:)+rt_fr2(:)+rt_sn1(:)+rt_sn2(:)+ &
3672 rt_sn3(:)+rt_pe1(:)+rt_pe2(:) )*ms2mmstp
3673 SNOWNCV(its:ite,j) = (rt_sn1(:) + rt_sn2(:))*ms2mmstp
3674 HAILNCV(its:ite,j) = (rt_pe1(:) + rt_pe2(:))*ms2mmstp
3675 GRPLNCV(its:ite,j) = rt_sn3(:)*ms2mmstp
3677 RAINNC(its:ite,j) = RAINNC(its:ite,j) + RAINNCV(its:ite,j)
3678 SNOWNC(its:ite,j) = SNOWNC(its:ite,j) + SNOWNCV(its:ite,j)
3679 HAILNC(its:ite,j) = HAILNC(its:ite,j) + HAILNCV(its:ite,j)
3680 GRPLNC(its:ite,j) = GRPLNC(its:ite,j) + GRPLNCV(its:ite,j)
3681 SR(its:ite,j) = (SNOWNCV(its:ite,j)+HAILNCV(its:ite,j)+GRPLNCV(its:ite,j))/(RAINNCV(its:ite,j)+1.e-12)
3684 enddo !j_loop
3686 !--run-time stats:
3687 ! if (mod(itimestep*dt,statfreq)==0.) then
3688 ! t3d = th*pii
3689 ! print*, 'Bloc Stats -- AFTER micro call; my_2.25.1-b11)'
3690 ! write(6,'(a8,7e15.5)') 'Max qx: ',maxval(qc(its:ite,kts:kte,jts:jte)),maxval(qr(its:ite,kts:kte,jts:jte)),maxval(qi(its:ite,kts:kte,jts:jte)), &
3691 ! maxval(qs(its:ite,kts:kte,jts:jte)),maxval(qg(its:ite,kts:kte,jts:jte)),maxval(qh(its:ite,kts:kte,jts:jte))
3692 ! write(6,'(a8,6e15.5)') 'Max Nx: ',maxval(nc(its:ite,kts:kte,jts:jte)),maxval(nr(its:ite,kts:kte,jts:jte)),maxval(ni(its:ite,kts:kte,jts:jte)), &
3693 ! maxval(ns(its:ite,kts:kte,jts:jte)),maxval(ng(its:ite,kts:kte,jts:jte)),maxval(nh(its:ite,kts:kte,jts:jte))
3694 ! write(6,'(a8,7e15.5)') 'Min qx: ',minval(qc(its:ite,kts:kte,jts:jte)),minval(qr(its:ite,kts:kte,jts:jte)),minval(qi(its:ite,kts:kte,jts:jte)), &
3695 ! minval(qs(its:ite,kts:kte,jts:jte)),minval(qg(its:ite,kts:kte,jts:jte)),minval(qh(its:ite,kts:kte,jts:jte))
3696 ! write(6,'(a8,6e15.5)') 'Min Nx: ',minval(nc(its:ite,kts:kte,jts:jte)),minval(nr(its:ite,kts:kte,jts:jte)),minval(ni(its:ite,kts:kte,jts:jte)), &
3697 ! minval(ns(its:ite,kts:kte,jts:jte)),minval(ng(its:ite,kts:kte,jts:jte)),minval(nh(its:ite,kts:kte,jts:jte))
3698 ! write(6,'(a8,6e15.5)') 'W,T,qv: ',minval( w(its:ite,kts:kte,jts:jte)),maxval( w(its:ite,kts:kte,jts:jte)),minval(t3d(its:ite,kts:kte,jts:jte)),&
3699 ! maxval(t3d(its:ite,kts:kte,jts:jte)),minval(qv(its:ite,kts:kte,jts:jte)),maxval(qv(its:ite,kts:kte,jts:jte))
3700 ! write(6,'(a8,2e15.5)') 'Zet : ',maxval(Zet(its:ite,kts:kte,jts:jte))
3701 ! endif
3703 END SUBROUTINE mp_milbrandt2mom_driver
3705 !+---+-----------------------------------------------------------------+
3706 !________________________________________________________________________________________!
3708 END MODULE module_mp_milbrandt2mom