/wrfv2_fire/phys/module_mp_milbrandt2mom.F

Large files files are truncated, but you can click here to view the full file

!_______________________________________________________________________________!
!                                                                               !
!  This module contains the microphysics sub-driver for the 2-moment version of !
!  the Milbrandt-Yau (2005, JAS) microphysics scheme, along with all associated !
!  subprograms.  The main subroutine, 'mp_milbrandt2mom_main', is essentially   !
!  directly from the RPN-CMC physics library of the Canadian GEM model.  It is  !
!  called by the wrapper 'mp_milbrandt2mom_driver' which makes the necessary    !
!  adjustments to the calling parameters for the interface to WRF.              !
!                                                                               !
!  For questions, bug reports, etc. pertaining to the scheme, or to request     !
!  updates to the code (before the next offical WRF release) please contact     !
!  Jason Milbrandt (Environment Canada) at jason.milbrandt@ec.gc.ca             !
!                                                                               !
!  Last modified:  2011-03-02                                                   !
!_______________________________________________________________________________!

module my_fncs_mod

!==============================================================================!
!  The following functions are used by the schemes in the multimoment package. !
!                                                                              !
!  Package version:  2.19.0      (internal bookkeeping)                        !
!  Last modified  :  2009-04-27                                                !
!==============================================================================!

   implicit none

   private
   public  :: NccnFNC,SxFNC,gamma,gammaDP,gser,gammln,gammp,cfg,gamminc

   contains

!==============================================================================!

 REAL FUNCTION NccnFNC(Win,Tin,Pin,CCNtype)

!---------------------------------------------------------------------------!
! This function returns number concentration (activated aerosols) as a
! function of w,T,p, based on polynomial approximations of detailed
! approach using a hypergeometric function, following Cohard and Pinty (2000a).
!---------------------------------------------------------------------------!

  IMPLICIT NONE

! PASSING PARAMETERS:
  real,    intent(in) :: Win, Tin, Pin
  integer, intent(in) :: CCNtype

! LOCAL PARAMETERS:
  real :: T,p,x,y,a,b,c,d,e,f,g,h,T2,T3,T4,x2,x3,x4,p2

  x= log10(Win*100.);   x2= x*x;  x3= x2*x;  x4= x2*x2
  T= Tin - 273.15;      T2= T*T;  T3= T2*T;  T4= T2*T2
  p= Pin*0.01;          p2= p*p

  if (CCNtype==1) then  !Maritime

     a= 1.47e-9*T4 -6.944e-8*T3 -9.933e-7*T2 +2.7278e-4*T -6.6853e-4
     b=-1.41e-8*T4 +6.662e-7*T3 +4.483e-6*T2 -2.0479e-3*T +4.0823e-2
     c= 5.12e-8*T4 -2.375e-6*T3 +4.268e-6*T2 +3.9681e-3*T -3.2356e-1
     d=-8.25e-8*T4 +3.629e-6*T3 -4.044e-5*T2 +2.1846e-3*T +9.1227e-1
     e= 5.02e-8*T4 -1.973e-6*T3 +3.944e-5*T2 -9.0734e-3*T +1.1256e0
     f= -1.424e-6*p2 +3.631e-3*p -1.986
     g= -0.0212*x4 +0.1765*x3 -0.3770*x2 -0.2200*x +1.0081
     h= 2.47e-6*T3 -3.654e-5*T2 +2.3327e-3*T +0.1938
     y= a*x4 + b*x3 + c*x2 + d*x + e + f*g*h
     NccnFNC= 10.**min(2.,max(0.,y)) *1.e6                ![m-3]

  else if (CCNtype==2) then  !Continental

     a= 0.
     b= 0.
     c=-2.112e-9*T4 +3.9836e-8*T3 +2.3703e-6*T2 -1.4542e-4*T -0.0698
     d=-4.210e-8*T4 +5.5745e-7*T3 +1.8460e-5*T2 +9.6078e-4*T +0.7120
     e= 1.434e-7*T4 -1.6455e-6*T3 -4.3334e-5*T2 -7.6720e-3*T +1.0056
     f= 1.340e-6*p2 -3.5114e-3*p  +1.9453
     g= 4.226e-3*x4 -5.6012e-3*x3 -8.7846e-2*x2 +2.7435e-2*x +0.9932
     h= 5.811e-9*T4 +1.5589e-7*T3 -3.8623e-5*T2 +1.4471e-3*T +0.1496
     y= a*x4 +b*x3 +c*x2 + d*x + e + (f*g*h)
     NccnFNC= 10.**max(0.,y) *1.e6

  else

    print*, '*** STOPPED in MODULE ### NccnFNC  *** '
    print*, '    Parameter CCNtype incorrectly specified'
    stop

  endif

 END FUNCTION NccnFNC
!======================================================================!

   real FUNCTION SxFNC(Win,Tin,Pin,Qsw,Qsi,CCNtype,WRT)

!---------------------------------------------------------------------------!
! This function returns the peak supersaturation achieved during ascent with
! activation of CCN aerosols as a function of w,T,p, based on polynomial
! approximations of detailed approach using a hypergeometric function,
! following Cohard and Pinty (2000a).
!---------------------------------------------------------------------------!

 IMPLICIT NONE

! PASSING PARAMETERS:
  integer, intent(IN) :: WRT
  integer, intent(IN) :: CCNtype
  real,    intent(IN) :: Win, Tin, Pin, Qsw, Qsi

! LOCAL PARAMETERS:
  real   ::  Si,Sw,Qv,T,p,x,a,b,c,d,f,g,h,Pcorr,T2corr,T2,T3,T4,x2,x3,x4,p2
  real, parameter :: TRPL= 273.15

  x= log10(max(Win,1.e-20)*100.);   x2= x*x;  x3= x2*x;  x4= x2*x2
  T= Tin;                           T2= T*T;  T3= T2*T;  T4= T2*T2
  p= Pin*0.01;                      p2= p*p

  if (CCNtype==1) then  !Maritime

     a= -5.109e-7*T4 -3.996e-5*T3 -1.066e-3*T2 -1.273e-2*T +0.0659
     b=  2.014e-6*T4 +1.583e-4*T3 +4.356e-3*T2 +4.943e-2*T -0.1538
     c= -2.037e-6*T4 -1.625e-4*T3 -4.541e-3*T2 -5.118e-2*T +0.1428
     d=  3.812e-7*T4 +3.065e-5*T3 +8.795e-4*T2 +9.440e-3*T +6.14e-3
     f= -2.012e-6*p2 + 4.1913e-3*p    - 1.785e0
     g=  2.832e-1*x3 -5.6990e-1*x2 +5.1105e-1*x -4.1747e-4
     h=  1.173e-6*T3 +3.2174e-5*T2 -6.8832e-4*T +6.7888e-2
     Pcorr= f*g*h
     T2corr= 0.9581-4.449e-3*T-2.016e-4*T2-3.307e-6*T3-1.725e-8*T4

  else if (CCNtype==2) then  !Continental [computed for -35<T<-5C]

     a=  3.80e-5*T2 +1.65e-4*T +9.88e-2
     b= -7.38e-5*T2 -2.53e-3*T -3.23e-1
     c=  8.39e-5*T2 +3.96e-3*T +3.50e-1
     d= -1.88e-6*T2 -1.33e-3*T -3.73e-2
     f= -1.9761e-6*p2 + 4.1473e-3*p - 1.771e0
     g=  0.1539*x4 -0.5575*x3 +0.9262*x2 -0.3498*x -0.1293
     h=-8.035e-9*T4+3.162e-7*T3+1.029e-5*T2-5.931e-4*T+5.62e-2
     Pcorr= f*g*h
     T2corr= 0.98888-5.0525e-4*T-1.7598e-5*T2-8.3308e-8*T3

  else

    print*, '*** STOPPED in MODULE ### SxFNC  *** '
    print*, '    Parameter CCNtype incorrectly specified'
    stop

  endif

  Sw= (a*x3 + b*x2 +c*x + d) + Pcorr
  Sw= 1. + 0.01*Sw
  Qv= Qsw*Sw
  Si= Qv/Qsi
  Si= Si*T2corr
  if (WRT.eq.1) then
     SxFNC= Sw
  else
     SxFNC= Si
  endif
  if (Win.le.0.) SxFNC= 1.

 END function SxFNC
!======================================================================!

 real FUNCTION gamma(xx)

!  Modified from "Numerical Recipes"

  IMPLICIT NONE

! PASSING PARAMETERS:
  real, intent(IN) :: xx

! LOCAL PARAMETERS:
  integer  :: j
  real*8   :: ser,stp,tmp,x,y,cof(6),gammadp


  SAVE cof,stp
  DATA cof,stp/76.18009172947146d0,-86.50532032941677d0,               &
       24.01409824083091d0,-1.231739572450155d0,.1208650973866179d-2,  &
       -.5395239384953d-5,2.5066282746310005d0/
  x=dble(xx)
  y=x
  tmp=x+5.5d0
  tmp=(x+0.5d0)*log(tmp)-tmp
  ser=1.000000000190015d0
! do j=1,6   !original
  do j=1,4
!!do j=1,3   !gives result to within ~ 3 %
     y=y+1.d0
     ser=ser+cof(j)/y
  enddo
  gammadp=tmp+log(stp*ser/x)
  gammadp= exp(gammadp)

#if (DWORDSIZE == 8 && RWORDSIZE == 8)
  gamma  =      gammadp
#elif (DWORDSIZE == 8 && RWORDSIZE == 4)
  gamma  = sngl(gammadp)
#else
     This is a temporary hack assuming double precision is 8 bytes.
#endif

 END FUNCTION gamma
!======================================================================!
! ! !
! ! ! -- USED BY DIAGNOSTIC-ALPHA DOUBLE-MOMENT (SINGLE-PRECISION) VERSION --
! ! !      FOR FUTURE VERSIONS OF M-Y PACKAGE WITH, THIS S/R CAN BE USED
! ! !
! ! !  real FUNCTION diagAlpha(Dm,x)
! ! !
! ! !   IMPLICIT NONE
! ! !
! ! !   integer :: x
! ! !   real    :: Dm
! ! !   real, dimension(5) :: c1,c2,c3,c4
! ! !   real, parameter    :: pi = 3.14159265
! ! !   real, parameter    :: alphaMAX= 80.e0
! ! !   data c1 /19.0, 12.0, 4.5, 5.5, 3.7/
! ! !   data c2 / 0.6,  0.7, 0.5, 0.7, 0.3/
! ! !   data c3 / 1.8,  1.7, 5.0, 4.5, 9.0/
! ! !   data c4 /17.0, 11.0, 5.5, 8.5, 6.5/
! ! !   diagAlpha= c1(x)*tanh(c2(x)*(1.e3*Dm-c3(x)))+c4(x)
! ! !   if (x==5.and.Dm>0.008) diagAlpha= 1.e3*Dm-2.6
! ! !   diagAlpha= min(diagAlpha, alphaMAX)
! ! !
! ! !  END function diagAlpha
! ! !
! ! ! !======================================================================!
! ! !
! ! ! -- USED BY DIAGNOSTIC-ALPHA DOUBLE-MOMENT (SINGLE-PRECISION) VERSION --
! ! !      FOR FUTURE VERSIONS OF M-Y PACKAGE WITH, THIS S/R CAN BE USED
! ! !
! ! !  real FUNCTION solveAlpha(Q,N,Z,Cx,rho)
! ! !
! ! !  IMPLICIT NONE
! ! !
! ! ! ! PASSING PARAMETERS:
! ! !   real, intent(IN) :: Q, N, Z, Cx, rho
! ! !
! ! ! ! LOCAL PARAMETERS:
! ! !   real             :: a,g,a1,g1,g2,tmp1
! ! !   integer          :: i
! ! !   real, parameter  :: alphaMax= 40.
! ! !   real, parameter  :: epsQ    = 1.e-14
! ! !   real, parameter  :: epsN    = 1.e-3
! ! !   real, parameter  :: epsZ    = 1.e-32
! ! !
! ! ! !  Q         mass mixing ratio
! ! ! !  N         total concentration
! ! ! !  Z         reflectivity
! ! ! !  Cx        (pi/6)*RHOx
! ! ! !  rho       air density
! ! ! !  a         alpha (returned as solveAlpha)
! ! ! !  g         function g(a)= [(6+a)(5+a)(4+a)]/[(3+a)(2+a)(1+a)],
! ! ! !              where g = (Cx/(rho*Q))**2.*(Z*N)
! ! !
! ! !
! ! !   if (Q==0. .or. N==0. .or. Z==0. .or. Cx==0. .or. rho==0.) then
! ! !   ! For testing/debugging only; this module should never be called
! ! !   ! if the above condition is true.
! ! !     print*,'*** STOPPED in MODULE ### solveAlpha *** '
! ! !     print*,'*** : ',Q,N,Z,Cx*1.9099,rho
! ! !     stop
! ! !   endif
! ! !
! ! !   IF (Q>epsQ .and. N>epsN .and. Z>epsZ ) THEN
! ! !
! ! !      tmp1= Cx/(rho*Q)
! ! !      g   = tmp1*Z*tmp1*N    ! g = (Z*N)*[Cx / (rho*Q)]^2
! ! !
! ! !  !Note: The above order avoids OVERFLOW, since tmp1*tmp1 is very large
! ! !
! ! ! !----------------------------------------------------------!
! ! ! ! !Solve alpha numerically: (brute-force; for testing only)
! ! ! !      a= 0.
! ! ! !      g2= 999.
! ! ! !      do i=0,4000
! ! ! !         a1= i*0.01
! ! ! !         g1= (6.+a1)*(5.+a1)*(4.+a1)/((3.+a1)*(2.+a1)*(1.+a1))
! ! ! !         if(abs(g-g1)<abs(g-g2)) then
! ! ! !            a = a1
! ! ! !            g2= g1
! ! ! !         endif
! ! ! !      enddo
! ! ! !----------------------------------------------------------!
! ! !
! ! ! !Piecewise-polynomial approximation of g(a) to solve for a:  [2004-11-29]
! ! !      if (g>=20.) then
! ! !        a= 0.
! ! !      else
! ! !        g2= g*g
! ! !        if (g<20.  .and.g>=13.31) a= 3.3638e-3*g2 - 1.7152e-1*g + 2.0857e+0
! ! !        if (g<13.31.and.g>=7.123) a= 1.5900e-2*g2 - 4.8202e-1*g + 4.0108e+0
! ! !        if (g<7.123.and.g>=4.200) a= 1.0730e-1*g2 - 1.7481e+0*g + 8.4246e+0
! ! !        if (g<4.200.and.g>=2.946) a= 5.9070e-1*g2 - 5.7918e+0*g + 1.6919e+1
! ! !        if (g<2.946.and.g>=1.793) a= 4.3966e+0*g2 - 2.6659e+1*g + 4.5477e+1
! ! !        if (g<1.793.and.g>=1.405) a= 4.7552e+1*g2 - 1.7958e+2*g + 1.8126e+2
! ! !        if (g<1.405.and.g>=1.230) a= 3.0889e+2*g2 - 9.0854e+2*g + 6.8995e+2
! ! !        if (g<1.230) a= alphaMax
! ! !      endif
! ! !
! ! !      solveAlpha= max(0.,min(a,alphaMax))
! ! !
! ! !   ELSE
! ! !
! ! !      solveAlpha= 0.
! ! !
! ! !   ENDIF
! ! !
! ! !  END FUNCTION solveAlpha
!======================================================================!

 FUNCTION gammaDP(xx)

!  Modified from "Numerical Recipes"

  IMPLICIT NONE

! PASSING PARAMETERS:
  DOUBLE PRECISION, INTENT(IN) :: xx

! LOCAL PARAMETERS:
  DOUBLE PRECISION  :: gammaDP
  INTEGER  :: j
  DOUBLE PRECISION  :: ser,stp,tmp,x,y,cof(6)


  SAVE cof,stp
  DATA cof,stp/76.18009172947146d0,-86.50532032941677d0,               &
       24.01409824083091d0,-1.231739572450155d0,.1208650973866179d-2,  &
       -.5395239384953d-5,2.5066282746310005d0/
  x=xx
  y=x
  tmp=x+5.5d0
  tmp=(x+0.5d0)*log(tmp)-tmp
  ser=1.000000000190015d0
! do j=1,6   !original
  do j=1,4
!!do j=1,3   !gives result to within ~ 3 %
     y=y+1.d0
     ser=ser+cof(j)/y
  enddo
  gammaDP=tmp+log(stp*ser/x)
  gammaDP= exp(gammaDP)

 END FUNCTION gammaDP
!======================================================================!

 SUBROUTINE gser(gamser,a,x,gln)

! USES gammln

!   Returns the incomplete gamma function P(a,x) evaluated by its series
!   representation as gamser.  Also returns GAMMA(a) as gln.

 implicit none

 integer :: itmax
 real    :: a,gamser,gln,x,eps
 parameter (itmax=100, eps=3.e-7)
 integer :: n
 real :: ap,de1,summ

 gln=gammln(a)
 if(x.le.0.)then
    if(x.lt.0.)pause 'x <0 in gser'
    gamser=0.
    return
 endif
 ap=a
 summ=1./a
 de1=summ
 do n=1,itmax
    ap=ap+1.
    de1=de1*x/ap
    summ=summ+de1
    if(abs(de1).lt.abs(summ)*eps) goto 1
 enddo
 pause 'a too large, itmax too small in gser'
1 gamser=summ*exp(-x+a*log(x)-gln)
 return

END SUBROUTINE gser
!======================================================================!

 real FUNCTION gammln(xx)

!  Returns value of ln(GAMMA(xx)) for xx>0
!   (modified from "Numerical Recipes")

  IMPLICIT NONE

! PASSING PARAMETERS:
  real, intent(IN) :: xx

! LOCAL PARAMETERS:
  integer  :: j
  real*8   :: ser,stp,tmp,x,y,cof(6)

  SAVE cof,stp
  DATA cof,stp/76.18009172947146d0,-86.50532032941677d0,               &
       24.01409824083091d0,-1.231739572450155d0,.1208650973866179d-2,  &
       -.5395239384953d-5,2.5066282746310005d0/
  x=dble(xx)
  y=x
  tmp=x+5.5d0
  tmp=(x+0.5d0)*log(tmp)-tmp
  ser=1.000000000190015d0
  do j=1,6   !original
!  do j=1,4
     y=y+1.d0
     ser=ser+cof(j)/y
  enddo
#if (DWORDSIZE == 8 && RWORDSIZE == 8)
  gammln=       tmp+log(stp*ser/x)
#elif (DWORDSIZE == 8 && RWORDSIZE == 4)
  gammln= sngl( tmp+log(stp*ser/x)  )
#else
     This is a temporary hack assuming double precision is 8 bytes.
#endif

 END FUNCTION gammln
!======================================================================!

 real FUNCTION gammp(a,x)

! USES gcf,gser

! Returns the incomplete gamma function P(a,x)

 implicit none

 real :: a,x,gammcf,gamser,gln

 if(x.lt.0..or.a.le.0.) pause 'bad arguments in gammq'
 if(x.lt.a+1.)then
    call gser(gamser,a,x,gln)
    gammp=gamser
 else
    call cfg(gammcf,a,x,gln)
    gammp=1.-gammcf
 endif
 return

 END FUNCTION gammp
!======================================================================!

 SUBROUTINE cfg(gammcf,a,x,gln)

! USES gammln

! Returns the incomplete gamma function (Q(a,x) evaluated by tis continued fraction
! representation as gammcf.  Also returns ln(GAMMA(a)) as gln.  ITMAX is the maximum
! allowed number of iterations; EPS is the relative accuracy; FPMIN is a number near
! the smallest representable floating-point number.

 implicit none

 integer :: i,itmax
 real    :: a,gammcf,gln,x,eps,fpmin
 real    :: an,b,c,d,de1,h
 parameter (itmax=100,eps=3.e-7)

 gln=gammln(a)
 b=x+1.-a
 c=1./fpmin
 d=1./b
 h=d
 do i= 1,itmax
   an=-i*(i-a)
   b=b+2.
   d=an*d+b
   if(abs(d).lt.fpmin)d=fpmin
   c=b+an/c
 if(abs(c).lt.fpmin) c=fpmin
   d=1./d
   de1=d*c
   h=h*de1
   if(abs(de1-1.).lt.eps) goto 1
 enddo
 pause 'a too large, itmax too small in gcf'
1 gammcf=exp(-x+a*log(x)-gln)*h
 return

END SUBROUTINE cfg
!======================================================================!

 real FUNCTION gamminc(p,xmax)

! USES gammp, gammln
! Returns incomplete gamma function, gamma(p,xmax)= P(p,xmax)*GAMMA(p)
 real :: p,xmax
 gamminc= gammp(p,xmax)*exp(gammln(p))

 end FUNCTION gamminc

!======================================================================!
!  real function x_tothe_y(x,y)
!
!     implicit none
!     real, intent(in) :: x,y
!     x_tothe_y= exp(y*log(x))
!
!  end function x_tothe_y
!======================================================================!

end module my_fncs_mod

!________________________________________________________________________________________!

module my_sedi_mod

!================================================================================!
!  The following subroutines are used by the schemes in the multimoment package. !
!                                                                                !
!  Package version:  2.19.0     (internal bookkeeping)                           !
!  Last modified  :  2011-01-07                                                  !
!================================================================================!

   implicit none

  private
  public :: SEDI_main_1b,SEDI_main_2,countColumns

   contains

!=====================================================================================!
 SUBROUTINE SEDI_main_2(QX,NX,cat,Q,T,DE,iDE,gamfact,epsQ,epsN,afx,bfx,cmx,dmx,      &
                        ckQx1,ckQx2,ckQx4,LXP,ni,nk,VxMax,DxMax,dt,DZ,massFlux,      &
                        ktop_sedi,GRAV,massFlux3D)

!-------------------------------------------------------------------------------------!
!  DOUBLE-MOMENT version of sedimentation subroutine for categories whose
!  fall velocity equation is V(D) = gamfact * afx * D^bfx
!-------------------------------------------------------------------------------------!

! Passing parameters:
!
!  VAR   Description
!  ---   ------------
!  QX    mass mixing ratio of category x
!  NX    number concentration of category x
!  cat:  hydrometeor category:
!   1     rain
!   2     ice
!   3     snow
!   4     graupel
!   5     hail
!-------------------------------------------------------------------------------------!

  use my_fncs_mod

  implicit none

! PASSING PARAMETERS:
  real, dimension(:,:), intent(inout) :: QX,NX,Q,T
  real, dimension(:),    intent(out)  :: massFlux
  real, optional, dimension(:,:), intent(out) :: massFlux3D
  real, dimension(:,:), intent(in)    :: DE,iDE,DZ
  real, intent(in) :: epsQ,epsN,VxMax,LXP,afx,bfx,cmx,dmx,ckQx1,ckQx2,ckQx4,DxMax,dt,GRAV
  integer, dimension(:), intent(in)   :: ktop_sedi
  integer, intent(in)                 :: ni,nk,cat

! LOCAL PARAMETERS:
  logical                :: slabHASmass,locallim,QxPresent
  integer                :: nnn,a,i,k,counter,l,km1,kp1,ks,kw,idzmin
  integer, dimension(nk) :: flim_Q,flim_N
  integer, dimension(ni) :: activeColumn,npassx,ke
  real                   :: VqMax,VnMax,iLAMx,iLAMxB0,tmp1,tmp2,tmp3,Dx,iDxMax,icmx,     &
                            VincFact,ratio_Vn2Vq,zmax_Q,zmax_N,tempo,idmx,Nos_Thompson,  &
                            No_s,iLAMs
  real, dimension(ni,nk) :: VVQ,VVN,RHOQX,gamfact
  real, dimension(ni)    :: dzMIN,dtx,VxMaxx
  real, dimension(nk)    :: vp_Q,vp_N,zt_Q,zt_N,zb_Q,zb_N,dzi,Q_star,N_star
  real, dimension(0:nk)  :: zz
  real, parameter        :: epsilon = 1.e-2
  real, parameter        :: thrd    = 1./3.
  real, parameter        :: sxth    = 1./6.
  real, parameter        :: CoMAX   = 2.0

!-------------------------------------------------------------------------------------!

   massFlux = 0.

  !Factor to estimate increased V from size-sorting:
  ! - this factor should be higher for categories with more time-splitting, since Vmax
  !   increases after each sedimentation split step (to be tuned)
   VincFact = 1.
   if (present(massFlux3D)) massFlux3D= 0.  !(for use in MAIN for diagnostics)

  !Determine for which slabs and columns sedimentation should be computes:
   call countColumns(QX,ni,nk,epsQ,counter,activeColumn,ktop_sedi)

   ratio_Vn2Vq= ckQx2/ckQx1
   iDxMax= 1./DxMax
   icmx  = 1./cmx
   idmx  = 1./dmx
   ks    = nk
   ke    = ktop_sedi  !(i-array) - formerly ke=1; now depends on max. level with hydrometeor
   kw    = -1         !direction of vertical leveling; -1 implies nk is bottom

   VVQ  = 0.
   VVN  = 0.
   VqMax= 0.
   VnMax= 0.

   DO a= 1,counter
      i= activeColumn(a)

      VVQ(i,:) = 0.
      do k= ktop_sedi(i),nk  !formerly do k= 1,nk
         QxPresent =  (QX(i,k)>epsQ .and. NX(i,k)>epsN)
         if (QxPresent) VVQ(i,k)= calcVV()*ckQx1
         if (present(massFlux3D)) massFlux3D(i,k)= VVQ(i,k)*DE(i,k)*QX(i,k)  !(for use in MAIN)
      enddo  !k-loop
      Vxmaxx(i)= min( VxMax, maxval(VVQ(i,:))*VincFact )

     !note: dzMIN is min. value in column (not necessarily lowest layer in general)
      dzMIN(i) = minval(DZ(i,:))
      npassx(i)= max(1, nint( dt*Vxmaxx(i)/(CoMAX*dzMIN(i)) ))
      dtx(i)   = dt/float(npassx(i))

!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -!
      DO nnn= 1,npassx(i)

         locallim = (nnn==1)

         do k= ktop_sedi(i),nk  !formerly do k= 1,nk
           RHOQX(i,k) = DE(i,k)*QX(i,k)
           QxPresent  = (QX(i,k)>epsQ .and. NX(i,k)>epsN)
           if (QxPresent) then
              if (locallim) then     !to avoid re-computing VVQ on first pass
                 VVQ(i,k)= -VVQ(i,k)
              else
                 VVQ(i,k)= -calcVV()*ckQx1
              endif
              VVN(i,k)= VVQ(i,k)*ratio_Vn2Vq
              VqMax   = max(VxMAX,-VVQ(i,k))
              VnMax   = max(VxMAX,-VVN(i,k))
           else
              VVQ(i,k)= 0.
              VVN(i,k)= 0.
              VqMax   = 0.
              VnMax   = 0.
           endif
         enddo  !k-loop

        !sum instantaneous surface mass flux at each split step: (for division later)
         massFlux(i)= massFlux(i) - VVQ(i,nk)*DE(i,nk)*QX(i,nk)

     !-- Perform single split sedimentation step:
     !   (formerly by calls to s/r 'blg4sedi', a modified [JM] version of 'blg2.ftn')
         zz(ks)= 0.
         do k= ks,ke(i),kw
            zz(k+kw)= zz(k)+dz(i,k)
            dzi(k)  = 1./dz(i,k)
            vp_Q(k) = 0.
            vp_N(k) = 0.
         enddo

         do k=ks,ke(i),kw
            zb_Q(k)= zz(k) + VVQ(i,k)*dtx(i)
            zb_N(k)= zz(k) + VVN(i,k)*dtx(i)
         enddo

         zt_Q(ke(i))= zb_Q(ke(i)) + dz(i,ke(i))
         zt_N(ke(i))= zb_N(ke(i)) + dz(i,ke(i))
         do k= ks,ke(i)-kw,kw
            zb_Q(k)= min(zb_Q(k+kw)-epsilon*dz(i,k), zz(k)+VVQ(i,k)*dtx(i))
            zb_N(k)= min(zb_N(k+kw)-epsilon*dz(i,k), zz(k)+VVN(i,k)*dtx(i))
            zt_Q(k)= zb_Q(k+kw)
            zt_N(k)= zb_N(k+kw)
         enddo

         do k=ks,ke(i),kw    !formerly k=1,nk
            Q_star(k)= RHOQX(i,k)*dz(i,k)/(zt_Q(k)-zb_Q(k))
            N_star(k)=    NX(i,k)*dz(i,k)/(zt_N(k)-zb_N(k))
         enddo

         if (locallim) then
            zmax_Q= abs(VqMax*dtx(i))
            zmax_N= abs(VnMax*dtx(i))
            do l=ks,ke(i),kw
               flim_Q(l)= l
               flim_N(l)= l
               do k= l,ke(i),kw
                  if (zmax_Q.ge.zz(k)-zz(l+kw)) flim_Q(l)= k
                  if (zmax_N.ge.zz(k)-zz(l+kw)) flim_N(l)= k
               enddo
            enddo
         endif

         do l=ks,ke(i),kw
            do k=l,flim_Q(l),kw
               vp_Q(l)= vp_Q(l) + Q_star(k)*max(0.,min(zz(l+kw),zt_Q(k))-max(zz(l),zb_Q(k)))
            enddo
            do k=l,flim_N(l),kw
               vp_N(l)= vp_N(l) + N_star(k)*max(0.,min(zz(l+kw),zt_N(k))-max(zz(l),zb_N(k)))
            enddo
         enddo

         do k=ks,ke(i),kw
            RHOQX(i,k)= vp_Q(k)*dzi(k)
               NX(i,k)= vp_N(k)*dzi(k)
         enddo
     !--

         do k= ktop_sedi(i),nk  !formerly do k= 1,nk
           QX(i,k)= RHOQX(i,k)*iDE(i,k)

         !Prevent levels with zero N and nonzero Q and size-limiter:
           QxPresent=  (QX(i,k)>epsQ .and. NX(i,k)>epsN)
           if (QxPresent) then    !size limiter
              Dx= (DE(i,k)*QX(i,k)/(NX(i,k)*cmx))**idmx
              if (cat==1 .and. Dx>3.e-3) then
                 tmp1   =  Dx-3.e-3;   tmp1= tmp1*tmp1
                 tmp2   = (Dx/DxMAX);  tmp2= tmp2*tmp2*tmp2
                 NX(i,k)= NX(i,k)*max((1.+2.e4*tmp1),tmp2)
              else
                 NX(i,k)= NX(i,k)*(max(Dx,DxMAX)*iDxMAX)**dmx   !impose Dx_max
              endif
           else   !here, "QxPresent" implies correlated QX and NX
              Q(i,k) = Q(i,k) + QX(i,k)
              T(i,k) = T(i,k) - LXP*QX(i,k)   !LCP for rain; LSP for i,s,g,h
              QX(i,k)= 0.
              NX(i,k)= 0.
           endif

         enddo

       ENDDO  !nnn-loop
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -!
      !compute average mass flux during the full time step: (used to compute the
      !instantaneous sedimentation rate [liq. equiv. volume flux] in the main s/r)
       massFlux(i)= massFlux(i)/float(npassx(i))

    ENDDO  !a(i)-loop
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -!

CONTAINS

   real function calcVV()
   !Calculates portion of moment-weighted fall velocities
      iLAMx   = ((QX(i,k)*DE(i,k)/NX(i,k))*ckQx4)**idmx
      iLAMxB0 = iLAMx**bfx
      calcVV  = gamfact(i,k)*iLAMxB0
   end function calcVV

 END SUBROUTINE SEDI_main_2

!=====================================================================================!
 SUBROUTINE SEDI_main_1b(QX,cat,T,DE,iDE,gamfact,epsQ,afx,bfx,icmx,dmx,ckQx1,ckQx4, &
                         ni,nk,VxMax,DxMax,dt,DZ,massFlux,No_x,ktop_sedi,GRAV,      &
                         massFlux3D)

!-------------------------------------------------------------------------------------!
!  SINGLE-MOMENT version of sedimentation subroutine for categories whose
!  fall velocity equation is V(D) = gamfact * afx * D^bfx
!-------------------------------------------------------------------------------------!

! Passing parameters:
!
!  VAR   Description
!  ---   ------------
!  QX    mass mixing ratio of category x
!  cat:  hydrometeor category:
!   1     rain
!   2     ice
!   3     snow
!   4     graupel
!   5     hail
!-------------------------------------------------------------------------------------!

  use my_fncs_mod

  implicit none

! PASSING PARAMETERS:
  real, dimension(:,:), intent(inout) :: QX,T
  real, dimension(:),    intent(out)   :: massFlux
  real, optional, dimension(:,:), intent(out) :: massFlux3D
  real, dimension(:,:), intent(in)    :: DE,iDE,DZ
  real,    intent(in)    :: epsQ,VxMax,afx,bfx,icmx,dmx,ckQx1,ckQx4,DxMax,dt,GRAV,No_x
  integer, dimension(:), intent(in) :: ktop_sedi
  integer, intent(in)    :: ni,nk,cat !,ktop_sedi

! LOCAL PARAMETERS:
  logical                :: slabHASmass,locallim,QxPresent
  integer                :: nnn,a,i,k,counter,l,km1,kp1,ks,kw,idzmin !,ke
  integer, dimension(nk) :: flim_Q
  integer, dimension(ni) :: activeColumn,npassx,ke
  real                   :: VqMax,iLAMx,iLAMxB0,tmp1,tmp2,Dx,iDxMax,VincFact,NX,iNo_x,   &
                            zmax_Q,zmax_N,tempo
  real, dimension(ni,nk) :: VVQ,RHOQX,gamfact
  real, dimension(ni)    :: dzMIN,dtx,VxMaxx
  real, dimension(nk)    :: vp_Q,zt_Q,zb_Q,dzi,Q_star
  real, dimension(0:nk)  :: zz
  real, parameter        :: epsilon = 1.e-2
  real, parameter        :: thrd  = 1./3.
  real, parameter        :: sxth  = 1./6.
  real, parameter        :: CoMAX = 2.0

!-------------------------------------------------------------------------------------!

   massFlux= 0.
  !Factor to estimate increased V from size-sorting:
  ! - this factor should be higher for categories with more time-splitting, since Vmax
  !   increases after each sedimentation split step (to be tuned)
   VincFact= 1.
   if (present(massFlux3D)) massFlux3D= 0.  !(for use in MAIN for diagnostics)

  !Determine for which slabs and columns sedimentation should be computes:
   call countColumns(QX,ni,nk,epsQ,counter,activeColumn,ktop_sedi)
   iNo_x = 1./No_x
   iDxMax= 1./DxMax
   ks    = nk
   ke    = ktop_sedi  !(i-array) - old: ke=1
   kw    = -1         !direction of vertical leveling

   VVQ  = 0.
   VqMax= 0.

   DO a= 1,counter
      i= activeColumn(a)

      VVQ(i,:) = 0.
      do k= ktop_sedi(i),nk  !do k= 1,nk
         QxPresent =  (QX(i,k)>epsQ)
!        if (QxPresent) VVQ(i,k)= calcVV()*ckQx1

         if (QxPresent) then
            !ice:
              if (cat==2) then
                 NX    = 5.*exp(0.304*(273.15-max(233.,T(i,k))))
                 iLAMx = (ckQx4*QX(i,k)*DE(i,k)/NX)**thrd
            !snow:
              else if (cat==3) then
                 iNo_x = 1./min(2.e+8, 2.e+6*exp(-0.12*min(-0.001,T(i,k)-273.15)))
                 iLAMx = sqrt(sqrt(QX(i,k)*DE(i,k)*icmx*sxth*iNo_x))
            !rain, graupel, hail:
              else
                 iLAMx = sqrt(sqrt(QX(i,k)*DE(i,k)*icmx*sxth*iNo_x))
              endif
              VVQ(i,k) = -gamfact(i,k)*ckQx1*iLAMx**bfx
         !    VqMax    = max(VxMAX,-VVQ(i,k))
         endif
         if (present(massFlux3D)) massFlux3D(i,k)= -VVQ(i,k)*DE(i,k)*QX(i,k)  !(for use in MAIN)

      enddo  !k-loop
      Vxmaxx(i)= min( VxMax, maxval(VVQ(i,:))*VincFact )

     !note: dzMIN is min. value in column (not necessarily lowest layer in general)
      dzMIN(i) = minval(DZ(i,:))
      npassx(i)= max(1, nint( dt*Vxmaxx(i)/(CoMAX*dzMIN(i)) ))
      dtx(i)   = dt/float(npassx(i))

!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -!
      DO nnn= 1,npassx(i)

         locallim = (nnn==1)

         do k= ktop_sedi(i),nk  !do k= 1,nk
           RHOQX(i,k) = DE(i,k)*QX(i,k)
           QxPresent  = (QX(i,k)>epsQ)
            if (QxPresent) then
             !ice:
               if (cat==2) then
                  NX    = 5.*exp(0.304*(273.15-max(233.,T(i,k))))
                  iLAMx = (ckQx4*QX(i,k)*DE(i,k)/NX)**thrd
             !snow:
               else if (cat==3) then
                  iNo_x = 1./min(2.e+8, 2.e+6*exp(-0.12*min(-0.001,T(i,k)-273.15)))
                  iLAMx = sqrt(sqrt(QX(i,k)*DE(i,k)*icmx*sxth*iNo_x))
             !rain, graupel, hail:
               else
                  iLAMx = sqrt(sqrt(QX(i,k)*DE(i,k)*icmx*sxth*iNo_x))
               endif
               VVQ(i,k) = -gamfact(i,k)*ckQx1*iLAMx**bfx
               VqMax    = max(VxMAX,-VVQ(i,k))
            endif

         enddo  !k-loop

     !-- Perform single split sedimentation step:  (formerly by calls to s/r 'blg4sedi')
         zz(ks)= 0.
         do k= ks,ke(i),kw
            zz(k+kw)= zz(k)+dz(i,k)
            dzi(k)  = 1./dz(i,k)
            vp_Q(k) = 0.
         enddo

         do k=ks,ke(i),kw
            zb_Q(k)= zz(k) + VVQ(i,k)*dtx(i)
         enddo

         zt_Q(ke(i))= zb_Q(ke(i)) + dz(i,ke(i))
         do k= ks,ke(i)-kw,kw
            zb_Q(k)= min(zb_Q(k+kw)-epsilon*dz(i,k), zz(k)+VVQ(i,k)*dtx(i))
            zt_Q(k)= zb_Q(k+kw)
         enddo

         do k=ks,ke(i),kw    !k=1,nk
            Q_star(k)= RHOQX(i,k)*dz(i,k)/(zt_Q(k)-zb_Q(k))
         enddo

         if (locallim) then
            zmax_Q= abs(VqMax*dtx(i))
            do l=ks,ke(i),kw
               flim_Q(l)= l
               do k= l,ke(i),kw
                  if (zmax_Q.ge.zz(k)-zz(l+kw)) flim_Q(l)= k
               enddo
            enddo
         endif

         do l=ks,ke(i),kw
            do k=l,flim_Q(l),kw
               vp_Q(l)= vp_Q(l) + Q_star(k)*max(0.,min(zz(l+kw),zt_Q(k))-max(zz(l),zb_Q(k)))
            enddo
         enddo

         do k=ks,ke(i),kw
            RHOQX(i,k)= vp_Q(k)*dzi(k)
         enddo
     !--

         do k= ktop_sedi(i),nk  ! do k= 1,nk
           QX(i,k)= RHOQX(i,k)*iDE(i,k)
         enddo

        !sum instantaneous flux at each split step: (for division later)
         massFlux(i)= massFlux(i) - VVQ(i,nk)*DE(i,nk)*QX(i,nk)

       ENDDO  !nnn-loop
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -!
      !compute average flux during the full time step: (this will be used to compute
      ! the instantaneous sedimentation rate [volume flux] in the main s/r)
       massFlux(i)= massFlux(i)/float(npassx(i))

    ENDDO  !a(i)-loop
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -!

 END SUBROUTINE SEDI_main_1b

!=====================================================================================!
 SUBROUTINE countColumns(QX,ni,nk,minQX,counter,activeColumn,ktop_sedi)

! Searches the hydrometeor array QX(ni,nk) for non-zero (>minQX) values.
! Returns the array if i-indices (activeColumn) for the columns (i)
! which contain at least one non-zero value, as well as the number of such
! columns (counter).

  implicit none

!PASSING PARAMETERS:
  integer, intent(in)                   :: ni,nk !,ktop_sedi
  integer, dimension(:), intent(in)     :: ktop_sedi
  integer,                 intent(out)  :: counter
  integer, dimension(:), intent(out)    :: activeColumn
  real,    dimension(:,:), intent(in)   :: QX
  real,    intent(in)                   :: minQX

!LOCAL PARAMETERS:
  integer                               :: i !,k
  integer, dimension(ni)                :: k

!    k= ktop_sedi-1  !  k=0
   counter     = 0
   activeColumn= 0

   do i=1,ni
      k(i)= ktop_sedi(i)-1  !  k=0
      do
         k(i)=k(i)+1
         if (QX(i,k(i))>minQX) then
            counter=counter+1
            activeColumn(counter)=i
            k(i)=0
            exit
         else
            if (k(i)==nk) then
               k(i)=0
               exit
            endif
         endif
      enddo
   enddo  !i-loop

 END SUBROUTINE countColumns

!=====================================================================================!

end module my_sedi_mod

!________________________________________________________________________________________!

module my_dmom_mod

  implicit none

  private
  public :: mp_milbrandt2mom_main

  contains

!_______________________________________________________________________________________!

 SUBROUTINE mp_milbrandt2mom_main(W_omega,T,Q,QC,QR,QI,QN,QG,QH,NC,NR,NY,NN,NG,NH,PS,TM,  &
     QM,QCM,QRM,QIM,QNM,QGM,QHM,NCM,NRM,NYM,NNM,NGM,NHM,PSM,S,RT_rn1,RT_rn2,RT_fr1,RT_fr2,&
     RT_sn1,RT_sn2,RT_sn3,RT_pe1,RT_pe2,RT_peL,RT_snd,GZ,T_TEND,Q_TEND,QCTEND,QRTEND,     &
     QITEND,QNTEND,QGTEND,QHTEND,NCTEND,NRTEND,NYTEND,NNTEND,NGTEND,NHTEND,dt,NI,N,NK,    &
     J,KOUNT,CCNtype,precipDiag_ON,sedi_ON,warmphase_ON,autoconv_ON,icephase_ON,snow_ON,  &
     initN,dblMom_c,dblMom_r,dblMom_i,dblMom_s,dblMom_g,dblMom_h,Dm_c,Dm_r,Dm_i,Dm_s,     &
     Dm_g,Dm_h,ZET,ZEC,SLW,VIS,VIS1,VIS2,VIS3,h_CB,h_ML1,h_ML2,h_SN,SS01,SS02,SS03,SS04,  &
     SS05,SS06,SS07,SS08,SS09,SS10,SS11,SS12,SS13,SS14,SS15,SS16,SS17,SS18,SS19,SS20)


  use my_fncs_mod
  use my_sedi_mod
!--WRF:
    use module_model_constants, ONLY: CPD => cp, CPV => cpv, RGASD => r_d, RGASV => r_v, &
        EPS1 => EP_2, DELTA => EP_1, CAPPA => rcp, GRAV => g, CHLC => XLV, CHLF => XLF
!==

  implicit none

!CALLING PARAMETERS:
  integer,               intent(in)    :: NI,NK,N,J,KOUNT,CCNtype
  real,                  intent(in)    :: dt
  real, dimension(:),    intent(in)    :: PS,PSM
  real, dimension(:),    intent(out)   :: h_CB,h_ML1,h_ML2,h_SN
  real, dimension(:),    intent(out)   :: RT_rn1,RT_rn2,RT_fr1,RT_fr2,RT_sn1,RT_sn2,   &
                                          RT_sn3,RT_pe1,RT_pe2,RT_peL,ZEC,RT_snd
  real, dimension(:,:),  intent(in)    :: W_omega,S,GZ
  real, dimension(:,:),  intent(inout) :: T,Q,QC,QR,QI,QN,QG,QH,NC,NR,NY,NN,NG,NH,     &
        TM,QM,QCM,QRM,QIM,QNM,QGM,QHM,NCM,NRM,NYM,NNM,NGM,NHM
  real, dimension(:,:),  intent(out)   :: T_TEND,QCTEND,QRTEND,QITEND,QNTEND,          &
        QGTEND,QHTEND,Q_TEND,NCTEND,NRTEND,NYTEND,NNTEND,NGTEND,NHTEND,ZET,Dm_c,       &
        Dm_r,Dm_i,Dm_s,Dm_g,Dm_h,SLW,VIS,VIS1,VIS2,VIS3,SS01,SS02,SS03,SS04,SS05,SS06, &
        SS07,SS08,SS09,SS10,SS11,SS12,SS13,SS14,SS15,SS16,SS17,SS18,SS19,SS20
  logical,               intent(in)    :: dblMom_c,dblMom_r,dblMom_i,dblMom_s,         &
        dblMom_g,dblMom_h,precipDiag_ON,sedi_ON,icephase_ON,snow_ON,warmphase_ON,      &
        autoconv_ON,initN

!_______________________________________________________________________________________
!                                                                                       !
!                    Milbrandt-Yau Multimoment Bulk Microphysics Scheme                 !
!                              - double-moment version   -                              !
!_______________________________________________________________________________________!
!  Package version:   2.19.0      (internal bookkeeping)                                !
!  Last modified  :   2011-03-02                                                        !
!_______________________________________________________________________________________!
!
!  Author:
!       J. Milbrandt, McGill University (August 2004)
!
!  Major revisions:
!
!  001  J. Milbrandt  (Dec 2006) - Converted the full Milbrandt-Yau (2005) multimoment
!        (RPN)                     scheme to an efficient fixed-dispersion double-moment
!                                  version
!  002  J. Milbrandt  (Mar 2007) - Added options for single-moment/double-moment for
!                                  each hydrometeor category
!  003  J. Milbrandt  (Feb 2008) - Modified single-moment version for use in GEM-LAM-2.5
!  004  J. Milbrandt  (Nov 2008) - Modified double-moment version for use in 2010 Vancouver
!                                  Olympics GEM-LAM configuration
!  005  J. Milbrandt  (Aug 2009) - Modified (dmom) for PHY_v5.0.4, for use in V2010 system:
!                                  + reduced ice/snow capacitance to C=0.25D (from C=0.5D)
!                                  + added diagnostic fields (VIS, levels, etc.)
!                                  + added constraints to snow size distribution (No_s and
!                                    LAMDA_s limits, plus changed m-D parameters
!                                  + modified solid-to-liquid ratio calculation, based on
!                                    volume flux (and other changes)
!                                  + added back sedimentation of ice category
!                                  + modified condition for conversion of graupel to hail
!                                  + corrected bug it diagnostic "ice pellets" vs. "hail"
!                                  + minor bug corrections (uninitialized values, etc.)
!  006  J. Milbrandt  (Jan 2011) - Bug fixes and minor code clean-up from PHY_v5.1.3 version
!                                  + corrected latent heat constants in thermodynamic functions
!                                    (ABi and ABw) for sublimation and evaporation
!                                  + properly initialized variables No_g and No_h
!                                  + changed max ice crystal size (fallspeed) to 5 mm (2 m s-1)
!                                  + imposed maximum ice number concentration of 1.e+7 m-3
!                                  + removed unused supersaturation reduction
!
!  Object:
!          Computes changes to the temperature, water vapor mixing ratio, and the
!          mixing ratios and total number concentrations of six hydrometeor species
!          resulting from cloud microphysical interactions at saturated grid points.
!          Liquid and solid surface precipitation rates from sedimenting hydrometeor
!          categories are also computed.
!
!          This subroutine and the associated modules form the single/double-moment
!          switchable verion of the multimoment bulk microphysics package, the full
!          version of which is described in the references below.
!
!  References:  Milbrandt and Yau, (2005a), J. Atmos. Sci., vol.62, 3051-3064
!               --------- and ---, (2005b), J. Atmos. Sci., vol.62, 3065-3081
!               (and references therein)
!
!  Please report bugs to:  jason.milbrandt@ec.gc.ca
!_______________________________________________________________________________________!
!
! Arguments:         Description:                                         Units:
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -!
!            - Input -
!
! NI                 number of x-dir points (in local subdomain)
! NK                 number of vertical levels
! N                  not used (to be removed)
! J                  y-dir index (local subdomain)
! KOUNT              current model time step number
! dt                 model time step                                      [s]
! CCNtype            switch for airmass type
!                      1 = maritime                   --> N_c =  80 cm-3  (1-moment cloud)
!                      2 = continental 1              --> N_c = 200 cm-3     "       "
!                      3 = continental 2  (polluted)  --> N_c = 500 cm-3     "       "
!                      4 = land-sea-mask-dependent (TBA)
! W_omega            vertical velocity                                    [Pa s-1]
! S                  sigma (=p/p_sfc)
! GZ                 geopotential
! dblMom_(x)         logical switch for double(T)-single(F)-moment for category (x)
! precipDiag_ON      logical switch, .F. to suppress calc. of sfc precip types
! sedi_ON            logical switch, .F. to suppress sedimentation
! warmphase_ON       logical switch, .F. to suppress warm-phase (Part II)
! autoconv_ON        logical switch, .F. to supppress autoconversion (cld->rn)
! icephase_ON        logical switch, .F. to suppress ice-phase (Part I)
! snow_ON            logical switch, .F. to suppress snow initiation
!
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -!
!            - Input/Output -
!
! T                  air temperature at time (t*)                         [K]
! TM                 air temperature at time (t-dt)                       [K]
! Q                  water vapor mixing ratio at (t*)                     [kg kg-1]
! QM                 water vapor mixing ratio at (t-dt)                   [kg kg-1]
! PS                 surface pressure at time (t*)                        [Pa]
! PSM                surface pressure at time (t-dt)                      [Pa]
!
!  For x = (C,R,I,N,G,H):  C = cloud
!                          R = rain
!                          I = ice (pristine) [except 'NY', not 'NI']
!                          N = snow
!                          G = graupel
!                          H = hail
!
! Q(x)               mixing ratio for hydrometeor x at (t*)               [kg kg-1]
! Q(x)M              mixing ratio for hydrometeor x at (t-dt)             [kg kg-1]
! N(x)               total number concentration for hydrometeor x  (t*)   [m-3]
! N(x)M              total number concentration for hydrometeor x  (t-dt) [m-3]
!
! Note:  The arrays "VM" (e.g. variables TM,QM,QCM etc.) are declared as INTENT(INOUT)
!        such that their values are modified in the code [VM = 0.5*(VM + V)].
!        This is to approxiate the values at time level (t), which are needed by
!        this routine but are unavailable to the PHYSICS.  The new values are discared
!        by the calling routine ('vkuocon6.ftn').  However, care should be taken with
!        interfacing with other modelling systems.  For GEM/MC2, it does not matter if
!        VM is modified since the calling module passes back only the tendencies
!        (VTEND) to the model.

!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -!
!            - Output -
!
! Q_TEND             tendency for water vapor mixing ratio                [kg kg-1 s-1]
! T_TEND             tendency for air temperature                         [K s-1]
! Q(x)TEND           tendency for mixing ratio for hydrometeor x          [kg kg-1 s-1]
! N(x)TEND           tendency for number concentration for hydrometeor x  [m-3 s-1]
! Dm_(x)             mean-mass diameter for hydrometeor x                 [m]
! H_CB               height of cloud base                                 [m]
! h_ML1              height of first melting level from ground            [m]
! h_ML2              height of first melting level from top               [m]
! h_SN               height of snow level                                 [m]
! RT_rn1             precipitation rate (at sfc) of liquid rain           [m+3 m-2 s-1]
! RT_rn2             precipitation rate (at sfc) of liquid drizzle        [m+3 m-2 s-1]
! RT_fr1             precipitation rate (at sfc) of freezing rain         [m+3 m-2 s-1]
! RT_fr2             precipitation rate (at sfc) of freezing drizzle      [m+3 m-2 s-1]
! RT_sn1             precipitation rate (at sfc) of ice crystals (liq-eq) [m+3 m-2 s-1]
! RT_sn2             precipitation rate (at sfc) of snow    (liq-equiv)   [m+3 m-2 s-1]
! RT_sn3             precipitation rate (at sfc) of graupel (liq-equiv)   [m+3 m-2 s-1]
! RT_snd             precipitation rate (at sfc) of snow    (frozen)      [m+3 m-2 s-1]
! RT_pe1             precipitation rate (at sfc) of ice pellets (liq-eq)  [m+3 m-2 s-1]
! RT_pe2             precipitation rate (at sfc) of hail (total; liq-eq)  [m+3 m-2 s-1]
! RT_peL             precipitation rate (at sfc) of hail (large only)     [m+3 m-2 s-1]
! SSxx               S/S terms (for testing purposes)
! SLW                supercooled liquid water content                     [kg m-3]
! VIS                visibility resulting from fog, rain, snow            [m]
! VIS1               visibility component through liquid cloud (fog)      [m]
! VIS2               visibility component through rain                    [m]
! VIS3               visibility component through snow                    [m]
! ZET                total equivalent radar reflectivity                  [dBZ]
! ZEC                composite (column-max) of ZET                        [dBZ]
!_______________________________________________________________________________________!


!LOCAL VARIABLES:

 !Variables to count active grid points:
  logical :: log1,log2,log3,log4,doneK,rainPresent,calcDiag,CB_found,ML_found,      &
             SN_found
  logical, dimension(size(QC,dim=1),size(QC,dim=2)) :: activePoint
  integer, dimension(size(QC,dim=1)) :: ktop_sedi
  integer :: i,k,niter,ll,start

  real    :: tmp1,tmp2,tmp3,tmp4,tmp5,tmp6,tmp7,tmp8,tmp9,tmp10,                    &
       VDmax,NNUmax,X,D,DEL,QREVP,NuDEPSOR,NuCONTA,NuCONTB,NuCONTC,iMUkin,Ecg,Erg,  &
       NuCONT,GG,Na,Tcc,F1,F2,Kdiff,PSIa,Kn,source,sink,sour,ratio,qvs0,Kstoke,     &
       DELqvs,ft,esi,Si,Simax,Vq,Vn,Vz,LAMr,No_r_DM,No_i,No_s,No_g,No_h,D_sll,      &
       iABi,ABw,VENTr,VENTs,VENTg,VENTi,VENTh,Cdiff,Ka,MUdyn,MUkin,DEo,Ng_tail,     &
       gam,ScTHRD,Tc,mi,ff,Ec,Ntr,Dho,DMrain,Ech,DMice,DMsnow,DMgrpl,DMhail,        &
       ssat,Swmax,dey,Esh,Eii,Eis,Ess,Eig,Eih,FRAC,JJ,Dirg,Dirh,Dsrs,Dsrg,Dsrh,     &
       Dgrg,Dgrh,SIGc,L,TAU,DrAUT,DrINIT,Di,Ds,Dg,Dh,qFact,nFact,Ki,Rz,NgCNgh,      &
       vr0,vi0,vs0,vg0,vh0,Dc,Dr,QCLcs,QCLrs,QCLis,QCLcg,QCLrg,QCLig,NhCNgh,        &
       QCLch,QCLrh,QCLsh,QMLir,QMLsr,QMLgr,QMLhr,QCLih,QVDvg,QVDvh,QSHhr,           &
       QFZci,QNUvi,QVDvi,QCNis,QCNis1,QCNis2,QCLir,QCLri,QCNsg,QCLsr,QCNgh,         &
       QCLgr,QHwet,QVDvs,QFZrh,QIMsi,QIMgi,NMLhr,NVDvh,NCLir,NCLri,NCLrh,           &
       NCLch,NCLsr,NCLirg,NCLirh,NrFZrh,NhFZrh,NCLsrs,NCLsrg,NCLsrh,NCLgrg,         &
       NCLgrh,NVDvg,NMLgr,NiCNis,NsCNis,NVDvs,NMLsr,NCLsh,NCLss,NNUvi,NFZci,NVDvi,  &
       NCLis,NCLig,NCLih,NMLir,NCLrs,NCNsg,NCLcs,NCLcg,NIMsi,NIMgi,NCLgr,NCLrg,     &
       NSHhr,RCAUTR,RCACCR,CCACCR,CCSCOC,CCAUTR,CRSCOR,ALFx,des_pmlt,Ecs,des,ides,  &
       LAMx,iLAMx,iLAMxB0,Dx,ffx,iLAMc,iNCM,iNRM,iNYM,iNNM,iNGM,iLAMs_D3,           &
       iLAMg,iLAMg2,iLAMgB0,iLAMgB1,iLAMgB2,iLAMh,iLAMhB0,iLAMhB1,iLAMhB2,iNHM,     &
       iLAMi,iLAMi2,iLAMi3,iLAMi4,iLAMi5,iLAMiB0,iLAMiB1,iLAMiB2,iLAMr6,iLAMh2,     &
       iLAMs,iLAMs2,iLAMsB0,iLAMsB1,iLAMsB2,iLAMr,iLAMr2,iLAMr3,iLAMr4,iLAMr5,      &
       iLAMc2,iLAMc3,iLAMc4,iLAMc5,iLAMc6,iQCM,iQRM,iQIM,iQNM,iQGM,iQHM,iEih,iEsh,  &
       N_c,N_r,N_i,N_s,N_g,N_h,fluxV_i,fluxV_g,fluxV_s,rhos_mlt,fracLiq

 !Variables that only need to be calulated on the first step (and saved):
  real, save :: idt,iMUc,cmr,cmi,cms,cmg,cmh,icmr,icmi,icmg,icms,icmh,idew,idei,    &
       ideh,ideg,GC1,imso,icexc9,cexr1,cexr2,cexr3,No_s_SM,No_r,idms,imgo,icexs2,   &
       cexr4,cexr5,cexr6,cexr9,icexr9,ckQr1,ckQr2,ckQr3,ckQi1,ckQi2,ckQi3,ckQi4,    &
       icexi9,ckQs1,ckQs2,cexs1,cexs2,ckQg1,ckQg2,ckQg4,ckQh1,ckQh2,ckQh4,GR37,dms, &
       LCP,LFP,LSP,ck5,ck6,PI2,PIov4,PIov6,CHLS,iCHLF,cxr,cxi,Gzr,Gzi,Gzs,Gzg,Gzh,  &
       N_c_SM,iGC1,GC2,GC3,GC4,GC5,iGC5,GC6,GC7,GC8,GC11,GC12,GC13,GC14,iGR34,mso,  &
       GC15,GR1,GR3,GR13,GR14,GR15,GR17,GR31,iGR31,GR32,GR33,GR34,GR35,GR36,GI4,    &
       GI6,GI20,GI21,GI22,GI31,GI32,GI33,GI34,GI35,iGI31,GI11,GI36,GI37,GI40,iGG34, &
       GS09,GS11,GS12,GS13,iGS20,GS31,iGS31,GS32,GS33,GS34,GS35,GS36,GS40,iGS40,    &
       GS50,GG09,GG11,GG12,GG13,GG31,iGG31,GG32,GG33,GG34,GG35,GG36,GG40,iGG99,GH09,&
       GH11,GH12,GH13,GH31,GH32,GH33,GH40,GR50,GG50,iGH34,GH50,iGH99,iGH31,iGS34,   &
       iGS20_D3,GS40_D3,cms_D3,eds,fds,rfact_FvFm

!Size distribution parameters:
  real, parameter :: MUc      =  3.    !shape parameter for cloud
  real, parameter :: alpha_c  =  1.    !shape parameter for cloud
  real, parameter :: alpha_r  =  0.    !shape parameter for rain
  real, parameter :: alpha_i  =  0.    !shape parameter for ice
  real, parameter :: alpha_s  =  0.    !shape parameter for snow
  real, parameter :: alpha_g  =  0.    !shape parameter for graupel
  real, parameter :: alpha_h  =  0.    !shape parameter for hail
  real, parameter :: No_s_max =  1.e+8 !max. allowable intercept for snow [m-4]
  real, parameter :: lamdas_min= 500.  !min. allowable LAMDA_s [m-1]

 !For single-moment:
  real, parameter :: No_r_SM  =  1.e+7  !intercept parameter for rain    [m-4]
  real, parameter :: No_g_SM  =  4.e+6  !intercept parameter for graupel [m-4]
  real, parameter :: No_h_SM  =  1.e+5  !intercept parameter for hail    [m-4]
  !note: No_s = f(T), rather than a fixed value
  !------------------------------------!
  ! Symbol convention: (dist. params.) ! MY05: Milbrandt & Yau, 2005a,b (JAS)
  !       MY05    F94       CP00       ! F94:  Ferrier, 1994            (JAS)
  !       ------  --------  ------     ! CP00: Cohard & Pinty, 2000a,b  (QJGR)
  !       ALFx    ALPHAx    MUx-1      !
  !       MUx     (1)       ALPHAx     !
  !       ALFx+1  ALPHAx+1  MUx        !
  !-------------------…