!WRF:MODEL_LAYER:PHYSICS
!
MODULE module_mp_kessler 1
CONTAINS
!----------------------------------------------------------------
SUBROUTINE kessler( t, qv, qc, qr, rho, pii & 1
,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