!WRF:MODEL_LAYER:PHYSICS
!

MODULE module_mp_kessler

CONTAINS
!----------------------------------------------------------------
   SUBROUTINE kessler( t, qv, qc, qr, rho, pii                  &
                      ,dt_in, z, xlv, cp                        &
                      ,EP2,SVP1,SVP2,SVP3,SVPT0,rhowater        &
                      ,dz8w                                     &
                      ,RAINNC, RAINNCV                          &
                      ,ids,ide, jds,jde, kds,kde                & ! domain dims
                      ,ims,ime, jms,jme, kms,kme                & ! memory dims
                      ,its,ite, jts,jte, kts,kte                & ! tile   dims
                                                                )
!----------------------------------------------------------------
   IMPLICIT NONE
!----------------------------------------------------------------
   !  taken from the COMMAS code - WCS 10 May 1999.
   !  converted from FORTRAN 77 to 90, tiled, WCS 10 May 1999.
!----------------------------------------------------------------
   REAL    , PARAMETER ::  c1 = .001 
   REAL    , PARAMETER ::  c2 = .001 
   REAL    , PARAMETER ::  c3 = 2.2 
   REAL    , PARAMETER ::  c4 = .875 
   REAL    , PARAMETER ::  fudge = 1.0 
   REAL    , PARAMETER ::  mxfall = 10.0 
!----------------------------------------------------------------
   INTEGER,      INTENT(IN   )    :: ids,ide, jds,jde, kds,kde, &
                                     ims,ime, jms,jme, kms,kme, &
                                     its,ite, jts,jte, kts,kte
   REAL   ,      INTENT(IN   )    :: xlv, cp
   REAL   ,      INTENT(IN   )    :: EP2,SVP1,SVP2,SVP3,SVPT0
   REAL   ,      INTENT(IN   )    :: rhowater

   REAL, DIMENSION( ims:ime , kms:kme , jms:jme ),              &
         INTENT(INOUT) ::                                       &
                                                            t , &
                                                            qv, &
                                                            qc, &
                                                            qr

   REAL, DIMENSION( ims:ime , kms:kme , jms:jme ),              &
         INTENT(IN   ) ::                                       &
                                                           rho, &
                                                           pii, &
                                                          dz8w 

   REAL, DIMENSION( ims:ime , kms:kme , jms:jme ),              &
         INTENT(IN   ) ::                                    z

   REAL, INTENT(IN   ) :: dt_in

   REAL, DIMENSION( ims:ime , jms:jme ),                        &
         INTENT(INOUT) ::                               RAINNC, &
                                                       RAINNCV



   ! local variables

   REAL :: qrprod, ern, gam, rcgs, rcgsi
   REAL, DIMENSION( its:ite , kts:kte, jts:jte ) ::     prod
   REAL, DIMENSION(kts:kte) :: vt, prodk, vtden,rdzk,rhok,factor,rdzw
   INTEGER :: i,j,k
   INTEGER :: nfall, n, nfall_new
   REAL    :: qrr, pressure, temp, es, qvs, dz, dt
   REAL    :: f5, dtfall, rdz, product
   REAL    :: max_heating, max_condense, max_rain, maxqrp
   REAL    :: vtmax, ernmax, crmax, factorn, time_sediment
   REAL    :: qcr, factorr, ppt
   REAL, PARAMETER :: max_cr_sedimentation = 0.75
!----------------------------------------------------------------

   INTEGER :: imax, kmax

    dt = dt_in

!   f5 = 237.3 * 17.27 * 2.5e6 / cp 
    f5 = svp2*(svpt0-svp3)*xlv/cp
    ernmax = 0.
    maxqrp = -100.

!------------------------------------------------------------------------------
! parameters for the time split terminal advection
!------------------------------------------------------------------------------

      max_heating = 0.
      max_condense = 0.
      max_rain = 0.

!-----------------------------------------------------------------------------
! outer J loop for entire microphysics, outer i loop for sedimentation
!-----------------------------------------------------------------------------

  microphysics_outer_j_loop: DO j = jts, jte

  sedimentation_outer_i_loop: DO i = its,ite

!   vtmax = 0.
   crmax = 0.


!------------------------------------------------------------------------------
! Terminal velocity calculation and advection, set up coefficients and
! compute stable timestep
!------------------------------------------------------------------------------

   DO k = 1, kte
     prodk(k)   = qr(i,k,j)
     rhok(k) = rho(i,k,j)
     qrr = amax1(0.,qr(i,k,j)*0.001*rhok(k))
     vtden(k) = sqrt(rhok(1)/rhok(k))
     vt(k) = 36.34*(qrr**0.1364) * vtden(k)
!     vtmax = amax1(vt(k), vtmax)
     rdzw(k) = 1./dz8w(i,k,j)
     crmax = amax1(vt(k)*dt*rdzw(k),crmax)
   ENDDO
   DO k = 1, kte-1
     rdzk(k) = 1./(z(i,k+1,j) - z(i,k,j))
   ENDDO
   rdzk(kte) = 1./(z(i,kte,j) - z(i,kte-1,j))

   nfall = max(1,nint(0.5+crmax/max_cr_sedimentation))  ! courant number for big timestep.
   dtfall = dt / float(nfall)                           ! splitting so courant number for sedimentation
   time_sediment = dt                                   ! is stable

!------------------------------------------------------------------------------
! Terminal velocity calculation and advection
! Do a time split loop on this for stability.
!------------------------------------------------------------------------------

   column_sedimentation: DO WHILE ( nfall > 0 )

   time_sediment = time_sediment - dtfall
   DO k = 1, kte-1
     factor(k) = dtfall*rdzk(k)/rhok(k)
   ENDDO
   factor(kte) = dtfall*rdzk(kte)

   ppt=0.

      k = 1
      ppt=rhok(k)*prodk(k)*vt(k)*dtfall/rhowater
      RAINNCV(i,j)=ppt*1000.
      RAINNC(i,j)=RAINNC(i,j)+ppt*1000.  ! unit = mm
 
!------------------------------------------------------------------------------
! Time split loop, Fallout done with flux upstream
!------------------------------------------------------------------------------

      DO k = kts, kte-1
        prodk(k) = prodk(k) - factor(k)           &
                  * (rhok(k)*prodk(k)*vt(k)       &
                    -rhok(k+1)*prodk(k+1)*vt(k+1))
      ENDDO

      k = kte
      prodk(k) = prodk(k) - factor(k)*prodk(k)*vt(k)

!------------------------------------------------------------------------------
! compute new sedimentation velocity, and check/recompute new 
! sedimentation timestep if this isn't the last split step.
!------------------------------------------------------------------------------

      IF( nfall > 1 ) THEN ! this wasn't the last split sedimentation timestep

        nfall = nfall - 1
        crmax = 0.
        DO k = kts, kte 
          qrr = amax1(0.,prodk(k)*0.001*rhok(k))
          vt(k) = 36.34*(qrr**0.1364) * vtden(k)
!          vtmax = amax1(vt(k), vtmax)
          crmax = amax1(vt(k)*time_sediment*rdzw(k),crmax)
        ENDDO

        nfall_new = max(1,nint(0.5+crmax/max_cr_sedimentation))
        if (nfall_new /= nfall ) then
          nfall = nfall_new
          dtfall = time_sediment/nfall
        end if

      ELSE  ! this was the last timestep

        DO k=kts,kte
          prod(i,k,j) = prodk(k)
        ENDDO
        nfall = 0  ! exit condition for sedimentation loop

      END IF

   ENDDO column_sedimentation

   ENDDO sedimentation_outer_i_loop

!------------------------------------------------------------------------------
! Production of rain and deletion of qc
! Production of qc from supersaturation
! Evaporation of QR
!------------------------------------------------------------------------------

     DO k = kts, kte
     DO i = its, ite
       factorn = 1.0 / (1.+c3*dt*amax1(0.,qr(i,k,j))**c4)
       qrprod = qc(i,k,j) * (1.0 - factorn)           &
             + factorn*c1*dt*amax1(qc(i,k,j)-c2,0.)      
       rcgs = 0.001*rho(i,k,j)

       qc(i,k,j) = amax1(qc(i,k,j) - qrprod,0.)
       qr(i,k,j) = (qr(i,k,j) + prod(i,k,j)-qr(i,k,j))
       qr(i,k,j) = amax1(qr(i,k,j) + qrprod,0.)

       temp      = pii(i,k,j)*t(i,k,j)
       pressure = 1.000e+05 * (pii(i,k,j)**(1004./287.))
       gam = 2.5e+06/(1004.*pii(i,k,j))
!      qvs       = 380.*exp(17.27*(temp-273.)/(temp- 36.))/pressure
       es        = 1000.*svp1*exp(svp2*(temp-svpt0)/(temp-svp3))
       qvs       = ep2*es/(pressure-es)
!      prod(i,k,j) = (qv(i,k,j)-qvs) / (1.+qvs*f5/(temp-36.)**2)
       prod(i,k,j) = (qv(i,k,j)-qvs) / (1.+pressure/(pressure-es)*qvs*f5/(temp-svp3)**2)
       ern  = amin1(dt*(((1.6+124.9*(rcgs*qr(i,k,j))**.2046)   &
          *(rcgs*qr(i,k,j))**.525)/(2.55e8/(pressure*qvs)       &
          +5.4e5))*(dim(qvs,qv(i,k,j))/(rcgs*qvs)),             &
          amax1(-prod(i,k,j)-qc(i,k,j),0.),qr(i,k,j))

! Update all variables

       product = amax1(prod(i,k,j),-qc(i,k,j))
       t (i,k,j) = t(i,k,j) + gam*(product - ern)
       qv(i,k,j) = amax1(qv(i,k,j) - product + ern,0.)
       qc(i,k,j) =       qc(i,k,j) + product
       qr(i,k,j) = qr(i,k,j) - ern

     ENDDO
     ENDDO

  ENDDO  microphysics_outer_j_loop

  RETURN

  END SUBROUTINE kessler

END MODULE module_mp_kessler