#ifdef _ACCEL
# include "module_mp_wsm5_accel.F"
#else
#if ( RWORDSIZE == 4 )
# define VREC vsrec
# define VSQRT vssqrt
#else
# define VREC vrec
# define VSQRT vsqrt
#endif
!Including inline expansion statistical function
MODULE module_mp_wsm5 2
!
USE module_utility, ONLY: WRFU_Clock, WRFU_Alarm
USE module_domain, ONLY : HISTORY_ALARM, Is_alarm_tstep
USE module_mp_radar
!
REAL, PARAMETER, PRIVATE :: dtcldcr = 120. ! maximum time step for minor loops
REAL, PARAMETER, PRIVATE :: n0r = 8.e6 ! intercept parameter rain
REAL, PARAMETER, PRIVATE :: avtr = 841.9 ! a constant for terminal velocity of rain
REAL, PARAMETER, PRIVATE :: bvtr = 0.8 ! a constant for terminal velocity of rain
REAL, PARAMETER, PRIVATE :: r0 = .8e-5 ! 8 microm in contrast to 10 micro m
REAL, PARAMETER, PRIVATE :: peaut = .55 ! collection efficiency
REAL, PARAMETER, PRIVATE :: xncr = 3.e8 ! maritime cloud in contrast to 3.e8 in tc80
REAL, PARAMETER, PRIVATE :: xmyu = 1.718e-5 ! the dynamic viscosity kgm-1s-1
REAL, PARAMETER, PRIVATE :: avts = 11.72 ! a constant for terminal velocity of snow
REAL, PARAMETER, PRIVATE :: bvts = .41 ! a constant for terminal velocity of snow
REAL, PARAMETER, PRIVATE :: n0smax = 1.e11 ! maximum n0s (t=-90C unlimited)
REAL, PARAMETER, PRIVATE :: lamdarmax = 8.e4 ! limited maximum value for slope parameter of rain
REAL, PARAMETER, PRIVATE :: lamdasmax = 1.e5 ! limited maximum value for slope parameter of snow
REAL, PARAMETER, PRIVATE :: lamdagmax = 6.e4 ! limited maximum value for slope parameter of graupel
REAL, PARAMETER, PRIVATE :: dicon = 11.9 ! constant for the cloud-ice diamter
REAL, PARAMETER, PRIVATE :: dimax = 500.e-6 ! limited maximum value for the cloud-ice diamter
REAL, PARAMETER, PRIVATE :: n0s = 2.e6 ! temperature dependent intercept parameter snow
REAL, PARAMETER, PRIVATE :: alpha = .12 ! .122 exponen factor for n0s
REAL, PARAMETER, PRIVATE :: pfrz1 = 100. ! constant in Biggs freezing
REAL, PARAMETER, PRIVATE :: pfrz2 = 0.66 ! constant in Biggs freezing
REAL, PARAMETER, PRIVATE :: qcrmin = 1.e-9 ! minimun values for qr, qs, and qg
REAL, PARAMETER, PRIVATE :: eacrc = 1.0 ! Snow/cloud-water collection efficiency
REAL, SAVE :: &
qc0, qck1,bvtr1,bvtr2,bvtr3,bvtr4,g1pbr, &
g3pbr,g4pbr,g5pbro2,pvtr,eacrr,pacrr, &
precr1,precr2,xmmax,roqimax,bvts1, &
bvts2,bvts3,bvts4,g1pbs,g3pbs,g4pbs, &
g5pbso2,pvts,pacrs,precs1,precs2,pidn0r, &
pidn0s,xlv1,pacrc,pi, &
rslopermax,rslopesmax,rslopegmax, &
rsloperbmax,rslopesbmax,rslopegbmax, &
rsloper2max,rslopes2max,rslopeg2max, &
rsloper3max,rslopes3max,rslopeg3max
!
! Specifies code-inlining of fpvs function in WSM52D below. JM 20040507
!
CONTAINS
!===================================================================
!
SUBROUTINE wsm5(th, q, qc, qr, qi, qs & 2,6
,den, pii, p, delz &
,delt,g, cpd, cpv, rd, rv, t0c &
,ep1, ep2, qmin &
,XLS, XLV0, XLF0, den0, denr &
,cliq,cice,psat &
,rain, rainncv &
,snow, snowncv &
,sr &
,refl_10cm, diagflag, do_radar_ref &
,ids,ide, jds,jde, kds,kde &
,ims,ime, jms,jme, kms,kme &
,its,ite, jts,jte, kts,kte &
)
!-------------------------------------------------------------------
IMPLICIT NONE
!-------------------------------------------------------------------
!
INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde , &
ims,ime, jms,jme, kms,kme , &
its,ite, jts,jte, kts,kte
REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
INTENT(INOUT) :: &
th, &
q, &
qc, &
qi, &
qr, &
qs
REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
INTENT(IN ) :: &
den, &
pii, &
p, &
delz
REAL, INTENT(IN ) :: delt, &
g, &
rd, &
rv, &
t0c, &
den0, &
cpd, &
cpv, &
ep1, &
ep2, &
qmin, &
XLS, &
XLV0, &
XLF0, &
cliq, &
cice, &
psat, &
denr
REAL, DIMENSION( ims:ime , jms:jme ), &
INTENT(INOUT) :: rain, &
rainncv, &
sr
!+---+-----------------------------------------------------------------+
REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT):: & ! GT
refl_10cm
!+---+-----------------------------------------------------------------+
REAL, DIMENSION( ims:ime , jms:jme ), OPTIONAL, &
INTENT(INOUT) :: snow, &
snowncv
! LOCAL VAR
#ifdef XEON_OPTIMIZED_WSM5
# include "mic-wsm5-3-5-locvar.h"
#else
REAL, DIMENSION( its:ite , kts:kte ) :: t
REAL, DIMENSION( its:ite , kts:kte, 2 ) :: qci, qrs
CHARACTER*256 :: emess
INTEGER :: i,j,k
#endif
!+---+-----------------------------------------------------------------+
REAL, DIMENSION(kts:kte):: qv1d, t1d, p1d, qr1d, qs1d, dBZ
LOGICAL, OPTIONAL, INTENT(IN) :: diagflag
INTEGER, OPTIONAL, INTENT(IN) :: do_radar_ref
!+---+-----------------------------------------------------------------+
!+---+-----------------------------------------------------------------+
#ifndef XEON_OPTIMIZED_WSM5
DO j=jts,jte
DO k=kts,kte
DO i=its,ite
t(i,k)=th(i,k,j)*pii(i,k,j)
qci(i,k,1) = qc(i,k,j)
qci(i,k,2) = qi(i,k,j)
qrs(i,k,1) = qr(i,k,j)
qrs(i,k,2) = qs(i,k,j)
ENDDO
ENDDO
! Sending array starting locations of optional variables may cause
! troubles, so we explicitly change the call.
CALL wsm52D(t, q(ims,kms,j), qci, qrs &
,den(ims,kms,j) &
,p(ims,kms,j), delz(ims,kms,j) &
,delt,g, cpd, cpv, rd, rv, t0c &
,ep1, ep2, qmin &
,XLS, XLV0, XLF0, den0, denr &
,cliq,cice,psat &
,j &
,rain(ims,j),rainncv(ims,j) &
,sr(ims,j) &
,ids,ide, jds,jde, kds,kde &
,ims,ime, jms,jme, kms,kme &
,its,ite, jts,jte, kts,kte &
,snow,snowncv &
)
DO K=kts,kte
DO I=its,ite
th(i,k,j)=t(i,k)/pii(i,k,j)
qc(i,k,j) = qci(i,k,1)
qi(i,k,j) = qci(i,k,2)
qr(i,k,j) = qrs(i,k,1)
qs(i,k,j) = qrs(i,k,2)
ENDDO
ENDDO
!+---+-----------------------------------------------------------------+
IF ( PRESENT (diagflag) ) THEN
if (diagflag .and. do_radar_ref == 1) then
! WRITE(emess,*)'calling refl10cm_wsm5 ',its, jts
! CALL wrf_debug ( 0, emess )
DO I=its,ite
DO K=kts,kte
t1d(k)=th(i,k,j)*pii(i,k,j)
p1d(k)=p(i,k,j)
qv1d(k)=q(i,k,j)
qr1d(k)=qr(i,k,j)
qs1d(k)=qs(i,k,j)
ENDDO
call refl10cm_wsm5 (qv1d, qr1d, qs1d, &
t1d, p1d, dBZ, kts, kte, i, j)
do k = kts, kte
refl_10cm(i,k,j) = MAX(-35., dBZ(k))
enddo
ENDDO
endif
ENDIF
!+---+-----------------------------------------------------------------+
ENDDO
#else
! Code to call the XEON-optimized WSM5 is included here
! this also contains OpenMP directives, so they are
! removed from the driver in the case where WSM5SCHEME
! is selected.
# include "mic-wsm5-3-5-callsite.h"
IF ( PRESENT (diagflag) ) THEN
IF (diagflag .and. do_radar_ref == 1) then
!$OMP PARALLEL DO &
!$OMP PRIVATE ( i,j,k,t1d,p1d,qv1d,qr1d,qs1d,dBZ )
DO j=jts,jte
DO I=its,ite
DO K=kts,kte
t1d(k)=th(i,k,j)*pii(i,k,j)
p1d(k)=p(i,k,j)
qv1d(k)=q(i,k,j)
qr1d(k)=qr(i,k,j)
qs1d(k)=qs(i,k,j)
ENDDO
call refl10cm_wsm5 (qv1d, qr1d, qs1d, &
t1d, p1d, dBZ, kts, kte, i, j)
DO k = kts, kte
refl_10cm(i,k,j) = MAX(-35., dBZ(k))
ENDDO
ENDDO
ENDDO
!$OMP END PARALLEL DO
ENDIF
ENDIF
#endif
END SUBROUTINE wsm5
#ifndef XEON_OPTIMIZED_WSM5
!===================================================================
!
SUBROUTINE wsm52D(t, q & 3,10
,qci, qrs, den, p, delz &
,delt,g, cpd, cpv, rd, rv, t0c &
,ep1, ep2, qmin &
,XLS, XLV0, XLF0, den0, denr &
,cliq,cice,psat &
,lat &
,rain,rainncv &
,sr &
,ids,ide, jds,jde, kds,kde &
,ims,ime, jms,jme, kms,kme &
,its,ite, jts,jte, kts,kte &
,snow,snowncv &
)
!-------------------------------------------------------------------
IMPLICIT NONE
!-------------------------------------------------------------------
!
! This code is a 5-class mixed ice microphyiscs scheme (WSM5) of the
! Single-Moment MicroPhyiscs (WSMMP). The WSMMP assumes that ice nuclei
! number concentration is a function of temperature, and seperate assumption
! is developed, in which ice crystal number concentration is a function
! of ice amount. A theoretical background of the ice-microphysics and related
! processes in the WSMMPs are described in Hong et al. (2004).
! Production terms in the WSM6 scheme are described in Hong and Lim (2006).
! All units are in m.k.s. and source/sink terms in kgkg-1s-1.
!
! WSM5 cloud scheme
!
! Coded by Song-You Hong (Yonsei Univ.)
! Jimy Dudhia (NCAR) and Shu-Hua Chen (UC Davis)
! Summer 2002
!
! Implemented by Song-You Hong (Yonsei Univ.) and Jimy Dudhia (NCAR)
! Summer 2003
!
! Reference) Hong, Dudhia, Chen (HDC, 2004) Mon. Wea. Rev.
! Rutledge, Hobbs (RH83, 1983) J. Atmos. Sci.
! Hong and Lim (HL, 2006) J. Korean Meteor. Soc.
!
INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde , &
ims,ime, jms,jme, kms,kme , &
its,ite, jts,jte, kts,kte, &
lat
REAL, DIMENSION( its:ite , kts:kte ), &
INTENT(INOUT) :: &
t
REAL, DIMENSION( its:ite , kts:kte, 2 ), &
INTENT(INOUT) :: &
qci, &
qrs
REAL, DIMENSION( ims:ime , kms:kme ), &
INTENT(INOUT) :: &
q
REAL, DIMENSION( ims:ime , kms:kme ), &
INTENT(IN ) :: &
den, &
p, &
delz
REAL, INTENT(IN ) :: delt, &
g, &
cpd, &
cpv, &
t0c, &
den0, &
rd, &
rv, &
ep1, &
ep2, &
qmin, &
XLS, &
XLV0, &
XLF0, &
cliq, &
cice, &
psat, &
denr
REAL, DIMENSION( ims:ime ), &
INTENT(INOUT) :: rain, &
rainncv, &
sr
REAL, DIMENSION( ims:ime, jms:jme ), OPTIONAL, &
INTENT(INOUT) :: snow, &
snowncv
! LOCAL VAR
REAL, DIMENSION( its:ite , kts:kte , 2) :: &
rh, &
qs, &
rslope, &
rslope2, &
rslope3, &
rslopeb, &
qrs_tmp, &
falk, &
fall, &
work1
REAL, DIMENSION( its:ite , kts:kte ) :: &
falkc, &
fallc, &
xl, &
cpm, &
denfac, &
xni, &
denqrs1, &
denqrs2, &
denqci, &
n0sfac, &
work2, &
workr, &
works, &
work1c, &
work2c
REAL, DIMENSION( its:ite , kts:kte ) :: &
den_tmp, &
delz_tmp
REAL, DIMENSION( its:ite ) :: &
delqrs1, &
delqrs2, &
delqi
REAL, DIMENSION( its:ite , kts:kte ) :: &
pigen, &
pidep, &
psdep, &
praut, &
psaut, &
prevp, &
psevp, &
pracw, &
psacw, &
psaci, &
pcond, &
psmlt
INTEGER, DIMENSION( its:ite ) :: &
mstep, &
numdt
REAL, DIMENSION(its:ite) :: tstepsnow
REAL, DIMENSION(its:ite) :: rmstep
REAL dtcldden, rdelz, rdtcld
LOGICAL, DIMENSION( its:ite ) :: flgcld
#define WSM_NO_CONDITIONAL_IN_VECTOR
#ifdef WSM_NO_CONDITIONAL_IN_VECTOR
REAL, DIMENSION(its:ite) :: xal, xbl
#endif
REAL :: &
cpmcal, xlcal, diffus, &
viscos, xka, venfac, conden, diffac, &
x, y, z, a, b, c, d, e, &
qdt, holdrr, holdrs, supcol, supcolt, pvt, &
coeres, supsat, dtcld, xmi, eacrs, satdt, &
vt2i,vt2s,acrfac, &
qimax, diameter, xni0, roqi0, &
fallsum, fallsum_qsi, xlwork2, factor, source, &
value, xlf, pfrzdtc, pfrzdtr, supice, holdc, holdci
! variables for optimization
REAL, DIMENSION( its:ite ) :: tvec1
REAL :: temp
INTEGER :: i, j, k, mstepmax, &
iprt, latd, lond, loop, loops, ifsat, n, idim, kdim
! Temporaries used for inlining fpvs function
REAL :: dldti, xb, xai, tr, xbi, xa, hvap, cvap, hsub, dldt, ttp
REAL :: logtr
!
!=================================================================
! compute internal functions
!
cpmcal(x) = cpd*(1.-max(x,qmin))+max(x,qmin)*cpv
xlcal(x) = xlv0-xlv1*(x-t0c)
!----------------------------------------------------------------
! diffus: diffusion coefficient of the water vapor
! viscos: kinematic viscosity(m2s-1)
! diffus(x,y) = 8.794e-5 * exp(log(x)*(1.81)) / y ! 8.794e-5*x**1.81/y
! viscos(x,y) = 1.496e-6 * (x*sqrt(x)) /(x+120.)/y ! 1.496e-6*x**1.5/(x+120.)/y
! xka(x,y) = 1.414e3*viscos(x,y)*y
! diffac(a,b,c,d,e) = d*a*a/(xka(c,d)*rv*c*c)+1./(e*diffus(c,b))
! venfac(a,b,c) = exp(log((viscos(b,c)/diffus(b,a)))*((.3333333))) &
! /sqrt(viscos(b,c))*sqrt(sqrt(den0/c))
! conden(a,b,c,d,e) = (max(b,qmin)-c)/(1.+d*d/(rv*e)*c/(a*a))
!
!----------------------------------------------------------------
idim = ite-its+1
kdim = kte-kts+1
!
!----------------------------------------------------------------
! paddint 0 for negative values generated by dynamics
!
do k = kts, kte
do i = its, ite
qci(i,k,1) = max(qci(i,k,1),0.0)
qrs(i,k,1) = max(qrs(i,k,1),0.0)
qci(i,k,2) = max(qci(i,k,2),0.0)
qrs(i,k,2) = max(qrs(i,k,2),0.0)
enddo
enddo
!
!----------------------------------------------------------------
! latent heat for phase changes and heat capacity. neglect the
! changes during microphysical process calculation
! emanuel(1994)
!
do k = kts, kte
do i = its, ite
cpm(i,k) = cpmcal(q(i,k))
xl(i,k) = xlcal(t(i,k))
enddo
enddo
do k = kts, kte
do i = its, ite
delz_tmp(i,k) = delz(i,k)
den_tmp(i,k) = den(i,k)
enddo
enddo
!
!----------------------------------------------------------------
! initialize the surface rain, snow
!
do i = its, ite
rainncv(i) = 0.
if(PRESENT (snowncv) .AND. PRESENT (snow)) snowncv(i,lat) = 0.
sr(i) = 0.
! new local array to catch step snow
tstepsnow(i) = 0.
enddo
!
!----------------------------------------------------------------
! compute the minor time steps.
!
loops = max(nint(delt/dtcldcr),1)
dtcld = delt/loops
if(delt.le.dtcldcr) dtcld = delt
!
do loop = 1,loops
!
!----------------------------------------------------------------
! initialize the large scale variables
!
do i = its, ite
mstep(i) = 1
flgcld(i) = .true.
enddo
!
! do k = kts, kte
! do i = its, ite
! denfac(i,k) = sqrt(den0/den(i,k))
! enddo
! enddo
do k = kts, kte
CALL VREC( tvec1(its), den(its,k), ite-its+1)
do i = its, ite
tvec1(i) = tvec1(i)*den0
enddo
CALL VSQRT( denfac(its,k), tvec1(its), ite-its+1)
enddo
!
! Inline expansion for fpvs
! qs(i,k,1) = fpvs(t(i,k),0,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
! qs(i,k,2) = fpvs(t(i,k),1,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
hsub = xls
hvap = xlv0
cvap = cpv
ttp=t0c+0.01
dldt=cvap-cliq
xa=-dldt/rv
xb=xa+hvap/(rv*ttp)
dldti=cvap-cice
xai=-dldti/rv
xbi=xai+hsub/(rv*ttp)
! this is for compilers where the conditional inhibits vectorization
#ifdef WSM_NO_CONDITIONAL_IN_VECTOR
do k = kts, kte
do i = its, ite
if(t(i,k).lt.ttp) then
xal(i) = xai
xbl(i) = xbi
else
xal(i) = xa
xbl(i) = xb
endif
enddo
do i = its, ite
tr=ttp/t(i,k)
logtr=log(tr)
qs(i,k,1)=psat*exp(logtr*(xa)+xb*(1.-tr))
qs(i,k,1) = min(qs(i,k,1),0.99*p(i,k))
qs(i,k,1) = ep2 * qs(i,k,1) / (p(i,k) - qs(i,k,1))
qs(i,k,1) = max(qs(i,k,1),qmin)
rh(i,k,1) = max(q(i,k) / qs(i,k,1),qmin)
qs(i,k,2)=psat*exp(logtr*(xal(i))+xbl(i)*(1.-tr))
qs(i,k,2) = min(qs(i,k,2),0.99*p(i,k))
qs(i,k,2) = ep2 * qs(i,k,2) / (p(i,k) - qs(i,k,2))
qs(i,k,2) = max(qs(i,k,2),qmin)
rh(i,k,2) = max(q(i,k) / qs(i,k,2),qmin)
#else
do k = kts, kte
do i = its, ite
tr=ttp/t(i,k)
logtr=log(tr)
qs(i,k,1)=psat*exp(logtr*(xa)+xb*(1.-tr))
qs(i,k,1) = min(qs(i,k,1),0.99*p(i,k))
qs(i,k,1) = ep2 * qs(i,k,1) / (p(i,k) - qs(i,k,1))
qs(i,k,1) = max(qs(i,k,1),qmin)
rh(i,k,1) = max(q(i,k) / qs(i,k,1),qmin)
if(t(i,k).lt.ttp) then
qs(i,k,2)=psat*exp(logtr*(xai)+xbi*(1.-tr))
else
qs(i,k,2)=psat*exp(logtr*(xa)+xb*(1.-tr))
endif
qs(i,k,2) = min(qs(i,k,2),0.99*p(i,k))
qs(i,k,2) = ep2 * qs(i,k,2) / (p(i,k) - qs(i,k,2))
qs(i,k,2) = max(qs(i,k,2),qmin)
rh(i,k,2) = max(q(i,k) / qs(i,k,2),qmin)
enddo
enddo
#endif
enddo
enddo
!
!----------------------------------------------------------------
! initialize the variables for microphysical physics
!
!
do k = kts, kte
do i = its, ite
prevp(i,k) = 0.
psdep(i,k) = 0.
praut(i,k) = 0.
psaut(i,k) = 0.
pracw(i,k) = 0.
psaci(i,k) = 0.
psacw(i,k) = 0.
pigen(i,k) = 0.
pidep(i,k) = 0.
pcond(i,k) = 0.
psmlt(i,k) = 0.
psevp(i,k) = 0.
falk(i,k,1) = 0.
falk(i,k,2) = 0.
fall(i,k,1) = 0.
fall(i,k,2) = 0.
fallc(i,k) = 0.
falkc(i,k) = 0.
xni(i,k) = 1.e3
enddo
enddo
!-------------------------------------------------------------
! Ni: ice crystal number concentraiton [HDC 5c]
!-------------------------------------------------------------
do k = kts, kte
do i = its, ite
temp = (den(i,k)*max(qci(i,k,2),qmin))
temp = sqrt(sqrt(temp*temp*temp))
xni(i,k) = min(max(5.38e7*temp,1.e3),1.e6)
enddo
enddo
!
!----------------------------------------------------------------
! compute the fallout term:
! first, vertical terminal velosity for minor loops
!----------------------------------------------------------------
do k = kts, kte
do i = its, ite
qrs_tmp(i,k,1) = qrs(i,k,1)
qrs_tmp(i,k,2) = qrs(i,k,2)
enddo
enddo
call slope_wsm5(qrs_tmp,den_tmp,denfac,t,rslope,rslopeb,rslope2,rslope3, &
work1,its,ite,kts,kte)
!
do k = kte, kts, -1
do i = its, ite
workr(i,k) = work1(i,k,1)
works(i,k) = work1(i,k,2)
denqrs1(i,k) = den(i,k)*qrs(i,k,1)
denqrs2(i,k) = den(i,k)*qrs(i,k,2)
if(qrs(i,k,1).le.0.0) workr(i,k) = 0.0
if(qrs(i,k,2).le.0.0) works(i,k) = 0.0
enddo
enddo
call nislfv_rain_plm(idim,kdim,den_tmp,denfac,t,delz_tmp,workr,denqrs1, &
delqrs1,dtcld,1,1)
call nislfv_rain_plm(idim,kdim,den_tmp,denfac,t,delz_tmp,works,denqrs2, &
delqrs2,dtcld,2,1)
do k = kts, kte
do i = its, ite
qrs(i,k,1) = max(denqrs1(i,k)/den(i,k),0.)
qrs(i,k,2) = max(denqrs2(i,k)/den(i,k),0.)
fall(i,k,1) = denqrs1(i,k)*workr(i,k)/delz(i,k)
fall(i,k,2) = denqrs2(i,k)*works(i,k)/delz(i,k)
enddo
enddo
do i = its, ite
fall(i,1,1) = delqrs1(i)/delz(i,1)/dtcld
fall(i,1,2) = delqrs2(i)/delz(i,1)/dtcld
enddo
do k = kts, kte
do i = its, ite
qrs_tmp(i,k,1) = qrs(i,k,1)
qrs_tmp(i,k,2) = qrs(i,k,2)
enddo
enddo
call slope_wsm5(qrs_tmp,den_tmp,denfac,t,rslope,rslopeb,rslope2,rslope3, &
work1,its,ite,kts,kte)
do k = kte, kts, -1
do i = its, ite
supcol = t0c-t(i,k)
n0sfac(i,k) = max(min(exp(alpha*supcol),n0smax/n0s),1.)
if(t(i,k).gt.t0c.and.qrs(i,k,2).gt.0.) then
!----------------------------------------------------------------
! psmlt: melting of snow [HL A33] [RH83 A25]
! (T>T0: S->R)
!----------------------------------------------------------------
xlf = xlf0
! work2(i,k)= venfac(p(i,k),t(i,k),den(i,k))
work2(i,k)= (exp(log(((1.496e-6*((t(i,k))*sqrt(t(i,k))) &
/((t(i,k))+120.)/(den(i,k)))/(8.794e-5 &
*exp(log(t(i,k))*(1.81))/p(i,k)))) &
*((.3333333)))/sqrt((1.496e-6*((t(i,k)) &
*sqrt(t(i,k)))/((t(i,k))+120.)/(den(i,k)))) &
*sqrt(sqrt(den0/(den(i,k)))))
coeres = rslope2(i,k,2)*sqrt(rslope(i,k,2)*rslopeb(i,k,2))
psmlt(i,k) = (1.414e3*(1.496e-6*((t(i,k))*sqrt(t(i,k))) &
/((t(i,k))+120.)/(den(i,k)) )*(den(i,k))) &
/xlf*(t0c-t(i,k))*pi/2. &
*n0sfac(i,k)*(precs1*rslope2(i,k,2)+precs2 &
*work2(i,k)*coeres)
psmlt(i,k) = min(max(psmlt(i,k)*dtcld/mstep(i), &
-qrs(i,k,2)/mstep(i)),0.)
qrs(i,k,2) = qrs(i,k,2) + psmlt(i,k)
qrs(i,k,1) = qrs(i,k,1) - psmlt(i,k)
t(i,k) = t(i,k) + xlf/cpm(i,k)*psmlt(i,k)
endif
enddo
enddo
!---------------------------------------------------------------
! Vice [ms-1] : fallout of ice crystal [HDC 5a]
!---------------------------------------------------------------
do k = kte, kts, -1
do i = its, ite
if(qci(i,k,2).le.0.) then
work1c(i,k) = 0.
else
xmi = den(i,k)*qci(i,k,2)/xni(i,k)
diameter = max(min(dicon * sqrt(xmi),dimax), 1.e-25)
work1c(i,k) = 1.49e4*exp(log(diameter)*(1.31))
endif
enddo
enddo
!
! forward semi-laglangian scheme (JH), PCM (piecewise constant), (linear)
!
do k = kte, kts, -1
do i = its, ite
denqci(i,k) = den(i,k)*qci(i,k,2)
enddo
enddo
call nislfv_rain_plm(idim,kdim,den_tmp,denfac,t,delz_tmp,work1c,denqci, &
delqi,dtcld,1,0)
do k = kts, kte
do i = its, ite
qci(i,k,2) = max(denqci(i,k)/den(i,k),0.)
enddo
enddo
do i = its, ite
fallc(i,1) = delqi(i)/delz(i,1)/dtcld
enddo
!
!----------------------------------------------------------------
! rain (unit is mm/sec;kgm-2s-1: /1000*delt ===> m)==> mm for wrf
!
do i = its, ite
fallsum = fall(i,1,1)+fall(i,1,2)+fallc(i,1)
fallsum_qsi = fall(i,1,2)+fallc(i,1)
if(fallsum.gt.0.) then
rainncv(i) = fallsum*delz(i,1)/denr*dtcld*1000. + rainncv(i)
rain(i) = fallsum*delz(i,1)/denr*dtcld*1000. + rain(i)
endif
if(fallsum_qsi.gt.0.) then
tstepsnow(i) = fallsum_qsi*delz(i,kts)/denr*dtcld*1000. &
+tstepsnow(i)
IF ( PRESENT (snowncv) .AND. PRESENT (snow)) THEN
snowncv(i,lat) = fallsum_qsi*delz(i,kts)/denr*dtcld*1000. &
+snowncv(i,lat)
snow(i,lat) = fallsum_qsi*delz(i,kts)/denr*dtcld*1000. + snow(i,lat)
ENDIF
endif
! if(fallsum.gt.0.)sr(i)=snowncv(i,lat)/(rainncv(i)+1.e-12)
if(fallsum.gt.0.)sr(i)=tstepsnow(i)/(rainncv(i)+1.e-12)
enddo
!
!---------------------------------------------------------------
! pimlt: instantaneous melting of cloud ice [HL A47] [RH83 A28]
! (T>T0: I->C)
!---------------------------------------------------------------
do k = kts, kte
do i = its, ite
supcol = t0c-t(i,k)
xlf = xls-xl(i,k)
if(supcol.lt.0.) xlf = xlf0
if(supcol.lt.0.and.qci(i,k,2).gt.0.) then
qci(i,k,1) = qci(i,k,1) + qci(i,k,2)
t(i,k) = t(i,k) - xlf/cpm(i,k)*qci(i,k,2)
qci(i,k,2) = 0.
endif
!---------------------------------------------------------------
! pihmf: homogeneous freezing of cloud water below -40c [HL A45]
! (T<-40C: C->I)
!---------------------------------------------------------------
if(supcol.gt.40..and.qci(i,k,1).gt.0.) then
qci(i,k,2) = qci(i,k,2) + qci(i,k,1)
t(i,k) = t(i,k) + xlf/cpm(i,k)*qci(i,k,1)
qci(i,k,1) = 0.
endif
!---------------------------------------------------------------
! pihtf: heterogeneous freezing of cloud water [HL A44]
! (T0>T>-40C: C->I)
!---------------------------------------------------------------
if(supcol.gt.0..and.qci(i,k,1).gt.0.) then
supcolt=min(supcol,50.)
! pfrzdtc = min(pfrz1*(exp(pfrz2*supcol)-1.) &
! *den(i,k)/denr/xncr*qci(i,k,1)**2*dtcld,qci(i,k,1))
pfrzdtc = min(pfrz1*(exp(pfrz2*supcolt)-1.) &
*den(i,k)/denr/xncr*qci(i,k,1)*qci(i,k,1)*dtcld,qci(i,k,1))
qci(i,k,2) = qci(i,k,2) + pfrzdtc
t(i,k) = t(i,k) + xlf/cpm(i,k)*pfrzdtc
qci(i,k,1) = qci(i,k,1)-pfrzdtc
endif
!---------------------------------------------------------------
! psfrz: freezing of rain water [HL A20] [LFO 45]
! (T<T0, R->S)
!---------------------------------------------------------------
if(supcol.gt.0..and.qrs(i,k,1).gt.0.) then
supcolt=min(supcol,50.)
! pfrzdtr = min(20.*pi**2*pfrz1*n0r*denr/den(i,k) &
! *(exp(pfrz2*supcol)-1.)*rslope(i,k,1)**7*dtcld, &
! qrs(i,k,1))
temp = rslope(i,k,1)
temp = temp*temp*temp*temp*temp*temp*temp
pfrzdtr = min(20.*(pi*pi)*pfrz1*n0r*denr/den(i,k) &
*(exp(pfrz2*supcolt)-1.)*temp*dtcld, &
qrs(i,k,1))
qrs(i,k,2) = qrs(i,k,2) + pfrzdtr
t(i,k) = t(i,k) + xlf/cpm(i,k)*pfrzdtr
qrs(i,k,1) = qrs(i,k,1)-pfrzdtr
endif
enddo
enddo
!
!----------------------------------------------------------------
! update the slope parameters for microphysics computation
!
do k = kts, kte
do i = its, ite
qrs_tmp(i,k,1) = qrs(i,k,1)
qrs_tmp(i,k,2) = qrs(i,k,2)
enddo
enddo
call slope_wsm5(qrs_tmp,den_tmp,denfac,t,rslope,rslopeb,rslope2,rslope3, &
work1,its,ite,kts,kte)
!------------------------------------------------------------------
! work1: the thermodynamic term in the denominator associated with
! heat conduction and vapor diffusion
! (ry88, y93, h85)
! work2: parameter associated with the ventilation effects(y93)
!
do k = kts, kte
do i = its, ite
! work1(i,k,1) = diffac(xl(i,k),p(i,k),t(i,k),den(i,k),qs(i,k,1))
work1(i,k,1) = ((((den(i,k))*(xl(i,k))*(xl(i,k)))*((t(i,k))+120.) &
*(den(i,k)))/(1.414e3*(1.496e-6*((t(i,k))*sqrt(t(i,k))))&
*(den(i,k))*(rv*(t(i,k))*(t(i,k))))) &
+ p(i,k)/((qs(i,k,1))*(8.794e-5*exp(log(t(i,k))*(1.81))))
! work1(i,k,2) = diffac(xls,p(i,k),t(i,k),den(i,k),qs(i,k,2))
work1(i,k,2) = ((((den(i,k))*(xls)*(xls))*((t(i,k))+120.)*(den(i,k)))&
/(1.414e3*(1.496e-6*((t(i,k))*sqrt(t(i,k))))*(den(i,k)) &
*(rv*(t(i,k))*(t(i,k)))) &
+ p(i,k)/(qs(i,k,2)*(8.794e-5*exp(log(t(i,k))*(1.81)))))
! work2(i,k) = venfac(p(i,k),t(i,k),den(i,k))
work2(i,k) = (exp(.3333333*log(((1.496e-6 * ((t(i,k))*sqrt(t(i,k)))) &
*p(i,k))/(((t(i,k))+120.)*den(i,k)*(8.794e-5 &
*exp(log(t(i,k))*(1.81))))))*sqrt(sqrt(den0/(den(i,k))))) &
/sqrt((1.496e-6*((t(i,k))*sqrt(t(i,k)))) &
/(((t(i,k))+120.)*den(i,k)))
enddo
enddo
!
!===============================================================
!
! warm rain processes
!
! - follows the processes in RH83 and LFO except for autoconcersion
!
!===============================================================
!
do k = kts, kte
do i = its, ite
supsat = max(q(i,k),qmin)-qs(i,k,1)
satdt = supsat/dtcld
!---------------------------------------------------------------
! praut: auto conversion rate from cloud to rain [HDC 16]
! (C->R)
!---------------------------------------------------------------
if(qci(i,k,1).gt.qc0) then
praut(i,k) = qck1*exp(log(qci(i,k,1))*((7./3.)))
praut(i,k) = min(praut(i,k),qci(i,k,1)/dtcld)
endif
!---------------------------------------------------------------
! pracw: accretion of cloud water by rain [HL A40] [LFO 51]
! (C->R)
!---------------------------------------------------------------
if(qrs(i,k,1).gt.qcrmin.and.qci(i,k,1).gt.qmin) then
pracw(i,k) = min(pacrr*rslope3(i,k,1)*rslopeb(i,k,1) &
*qci(i,k,1)*denfac(i,k),qci(i,k,1)/dtcld)
endif
!---------------------------------------------------------------
! prevp: evaporation/condensation rate of rain [HDC 14]
! (V->R or R->V)
!---------------------------------------------------------------
if(qrs(i,k,1).gt.0.) then
coeres = rslope2(i,k,1)*sqrt(rslope(i,k,1)*rslopeb(i,k,1))
prevp(i,k) = (rh(i,k,1)-1.)*(precr1*rslope2(i,k,1) &
+precr2*work2(i,k)*coeres)/work1(i,k,1)
if(prevp(i,k).lt.0.) then
prevp(i,k) = max(prevp(i,k),-qrs(i,k,1)/dtcld)
prevp(i,k) = max(prevp(i,k),satdt/2)
else
prevp(i,k) = min(prevp(i,k),satdt/2)
endif
endif
enddo
enddo
!
!===============================================================
!
! cold rain processes
!
! - follows the revised ice microphysics processes in HDC
! - the processes same as in RH83 and RH84 and LFO behave
! following ice crystal hapits defined in HDC, inclduing
! intercept parameter for snow (n0s), ice crystal number
! concentration (ni), ice nuclei number concentration
! (n0i), ice diameter (d)
!
!===============================================================
!
rdtcld = 1./dtcld
do k = kts, kte
do i = its, ite
supcol = t0c-t(i,k)
n0sfac(i,k) = max(min(exp(alpha*supcol),n0smax/n0s),1.)
supsat = max(q(i,k),qmin)-qs(i,k,2)
satdt = supsat/dtcld
ifsat = 0
!-------------------------------------------------------------
! Ni: ice crystal number concentraiton [HDC 5c]
!-------------------------------------------------------------
! xni(i,k) = min(max(5.38e7*(den(i,k) &
! *max(qci(i,k,2),qmin))**0.75,1.e3),1.e6)
temp = (den(i,k)*max(qci(i,k,2),qmin))
temp = sqrt(sqrt(temp*temp*temp))
xni(i,k) = min(max(5.38e7*temp,1.e3),1.e6)
eacrs = exp(0.07*(-supcol))
!
if(supcol.gt.0) then
if(qrs(i,k,2).gt.qcrmin.and.qci(i,k,2).gt.qmin) then
xmi = den(i,k)*qci(i,k,2)/xni(i,k)
diameter = min(dicon * sqrt(xmi),dimax)
vt2i = 1.49e4*diameter**1.31
vt2s = pvts*rslopeb(i,k,2)*denfac(i,k)
!-------------------------------------------------------------
! psaci: Accretion of cloud ice by rain [HDC 10]
! (T<T0: I->S)
!-------------------------------------------------------------
acrfac = 2.*rslope3(i,k,2)+2.*diameter*rslope2(i,k,2) &
+diameter**2*rslope(i,k,2)
psaci(i,k) = pi*qci(i,k,2)*eacrs*n0s*n0sfac(i,k) &
*abs(vt2s-vt2i)*acrfac/4.
endif
endif
!-------------------------------------------------------------
! psacw: Accretion of cloud water by snow [HL A7] [LFO 24]
! (T<T0: C->S, and T>=T0: C->R)
!-------------------------------------------------------------
if(qrs(i,k,2).gt.qcrmin.and.qci(i,k,1).gt.qmin) then
psacw(i,k) = min(pacrc*n0sfac(i,k)*rslope3(i,k,2) &
*rslopeb(i,k,2)*qci(i,k,1)*denfac(i,k) &
! ,qci(i,k,1)/dtcld)
,qci(i,k,1)*rdtcld)
endif
if(supcol .gt. 0) then
!-------------------------------------------------------------
! pidep: Deposition/Sublimation rate of ice [HDC 9]
! (T<T0: V->I or I->V)
!-------------------------------------------------------------
if(qci(i,k,2).gt.0.and.ifsat.ne.1) then
xmi = den(i,k)*qci(i,k,2)/xni(i,k)
diameter = dicon * sqrt(xmi)
pidep(i,k) = 4.*diameter*xni(i,k)*(rh(i,k,2)-1.)/work1(i,k,2)
supice = satdt-prevp(i,k)
if(pidep(i,k).lt.0.) then
! pidep(i,k) = max(max(pidep(i,k),satdt/2),supice)
! pidep(i,k) = max(pidep(i,k),-qci(i,k,2)/dtcld)
pidep(i,k) = max(max(pidep(i,k),satdt*.5),supice)
pidep(i,k) = max(pidep(i,k),-qci(i,k,2)*rdtcld)
else
! pidep(i,k) = min(min(pidep(i,k),satdt/2),supice)
pidep(i,k) = min(min(pidep(i,k),satdt*.5),supice)
endif
if(abs(prevp(i,k)+pidep(i,k)).ge.abs(satdt)) ifsat = 1
endif
!-------------------------------------------------------------
! psdep: deposition/sublimation rate of snow [HDC 14]
! (V->S or S->V)
!-------------------------------------------------------------
if(qrs(i,k,2).gt.0..and.ifsat.ne.1) then
coeres = rslope2(i,k,2)*sqrt(rslope(i,k,2)*rslopeb(i,k,2))
psdep(i,k) = (rh(i,k,2)-1.)*n0sfac(i,k) &
*(precs1*rslope2(i,k,2)+precs2 &
*work2(i,k)*coeres)/work1(i,k,2)
supice = satdt-prevp(i,k)-pidep(i,k)
if(psdep(i,k).lt.0.) then
! psdep(i,k) = max(psdep(i,k),-qrs(i,k,2)/dtcld)
! psdep(i,k) = max(max(psdep(i,k),satdt/2),supice)
psdep(i,k) = max(psdep(i,k),-qrs(i,k,2)*rdtcld)
psdep(i,k) = max(max(psdep(i,k),satdt*.5),supice)
else
! psdep(i,k) = min(min(psdep(i,k),satdt/2),supice)
psdep(i,k) = min(min(psdep(i,k),satdt*.5),supice)
endif
if(abs(prevp(i,k)+pidep(i,k)+psdep(i,k)).ge.abs(satdt)) &
ifsat = 1
endif
!-------------------------------------------------------------
! pigen: generation(nucleation) of ice from vapor [HL A50] [HDC 7-8]
! (T<T0: V->I)
!-------------------------------------------------------------
if(supsat.gt.0.and.ifsat.ne.1) then
supice = satdt-prevp(i,k)-pidep(i,k)-psdep(i,k)
xni0 = 1.e3*exp(0.1*supcol)
roqi0 = 4.92e-11*exp(log(xni0)*(1.33))
pigen(i,k) = max(0.,(roqi0/den(i,k)-max(qci(i,k,2),0.)) &
! /dtcld)
*rdtcld)
pigen(i,k) = min(min(pigen(i,k),satdt),supice)
endif
!
!-------------------------------------------------------------
! psaut: conversion(aggregation) of ice to snow [HDC 12]
! (T<T0: I->S)
!-------------------------------------------------------------
if(qci(i,k,2).gt.0.) then
qimax = roqimax/den(i,k)
! psaut(i,k) = max(0.,(qci(i,k,2)-qimax)/dtcld)
psaut(i,k) = max(0.,(qci(i,k,2)-qimax)*rdtcld)
endif
endif
!-------------------------------------------------------------
! psevp: Evaporation of melting snow [HL A35] [RH83 A27]
! (T>T0: S->V)
!-------------------------------------------------------------
if(supcol.lt.0.) then
if(qrs(i,k,2).gt.0..and.rh(i,k,1).lt.1.) &
psevp(i,k) = psdep(i,k)*work1(i,k,2)/work1(i,k,1)
! psevp(i,k) = min(max(psevp(i,k),-qrs(i,k,2)/dtcld),0.)
psevp(i,k) = min(max(psevp(i,k),-qrs(i,k,2)*rdtcld),0.)
endif
enddo
enddo
!
!
!----------------------------------------------------------------
! check mass conservation of generation terms and feedback to the
! large scale
!
do k = kts, kte
do i = its, ite
if(t(i,k).le.t0c) then
!
! cloud water
!
value = max(qmin,qci(i,k,1))
source = (praut(i,k)+pracw(i,k)+psacw(i,k))*dtcld
if (source.gt.value) then
factor = value/source
praut(i,k) = praut(i,k)*factor
pracw(i,k) = pracw(i,k)*factor
psacw(i,k) = psacw(i,k)*factor
endif
!
! cloud ice
!
value = max(qmin,qci(i,k,2))
source = (psaut(i,k)+psaci(i,k)-pigen(i,k)-pidep(i,k))*dtcld
if (source.gt.value) then
factor = value/source
psaut(i,k) = psaut(i,k)*factor
psaci(i,k) = psaci(i,k)*factor
pigen(i,k) = pigen(i,k)*factor
pidep(i,k) = pidep(i,k)*factor
endif
!
! rain
!
!
value = max(qmin,qrs(i,k,1))
source = (-praut(i,k)-pracw(i,k)-prevp(i,k))*dtcld
if (source.gt.value) then
factor = value/source
praut(i,k) = praut(i,k)*factor
pracw(i,k) = pracw(i,k)*factor
prevp(i,k) = prevp(i,k)*factor
endif
!
! snow
!
value = max(qmin,qrs(i,k,2))
source = (-psdep(i,k)-psaut(i,k)-psaci(i,k)-psacw(i,k))*dtcld
if (source.gt.value) then
factor = value/source
psdep(i,k) = psdep(i,k)*factor
psaut(i,k) = psaut(i,k)*factor
psaci(i,k) = psaci(i,k)*factor
psacw(i,k) = psacw(i,k)*factor
endif
!
work2(i,k)=-(prevp(i,k)+psdep(i,k)+pigen(i,k)+pidep(i,k))
! update
q(i,k) = q(i,k)+work2(i,k)*dtcld
qci(i,k,1) = max(qci(i,k,1)-(praut(i,k)+pracw(i,k) &
+psacw(i,k))*dtcld,0.)
qrs(i,k,1) = max(qrs(i,k,1)+(praut(i,k)+pracw(i,k) &
+prevp(i,k))*dtcld,0.)
qci(i,k,2) = max(qci(i,k,2)-(psaut(i,k)+psaci(i,k) &
-pigen(i,k)-pidep(i,k))*dtcld,0.)
qrs(i,k,2) = max(qrs(i,k,2)+(psdep(i,k)+psaut(i,k) &
+psaci(i,k)+psacw(i,k))*dtcld,0.)
xlf = xls-xl(i,k)
xlwork2 = -xls*(psdep(i,k)+pidep(i,k)+pigen(i,k)) &
-xl(i,k)*prevp(i,k)-xlf*psacw(i,k)
t(i,k) = t(i,k)-xlwork2/cpm(i,k)*dtcld
else
!
! cloud water
!
value = max(qmin,qci(i,k,1))
source=(praut(i,k)+pracw(i,k)+psacw(i,k))*dtcld
if (source.gt.value) then
factor = value/source
praut(i,k) = praut(i,k)*factor
pracw(i,k) = pracw(i,k)*factor
psacw(i,k) = psacw(i,k)*factor
endif
!
! rain
!
value = max(qmin,qrs(i,k,1))
source = (-praut(i,k)-pracw(i,k)-prevp(i,k)-psacw(i,k))*dtcld
if (source.gt.value) then
factor = value/source
praut(i,k) = praut(i,k)*factor
pracw(i,k) = pracw(i,k)*factor
prevp(i,k) = prevp(i,k)*factor
psacw(i,k) = psacw(i,k)*factor
endif
!
! snow
!
value = max(qcrmin,qrs(i,k,2))
source=(-psevp(i,k))*dtcld
if (source.gt.value) then
factor = value/source
psevp(i,k) = psevp(i,k)*factor
endif
work2(i,k)=-(prevp(i,k)+psevp(i,k))
! update
q(i,k) = q(i,k)+work2(i,k)*dtcld
qci(i,k,1) = max(qci(i,k,1)-(praut(i,k)+pracw(i,k) &
+psacw(i,k))*dtcld,0.)
qrs(i,k,1) = max(qrs(i,k,1)+(praut(i,k)+pracw(i,k) &
+prevp(i,k) +psacw(i,k))*dtcld,0.)
qrs(i,k,2) = max(qrs(i,k,2)+psevp(i,k)*dtcld,0.)
xlf = xls-xl(i,k)
xlwork2 = -xl(i,k)*(prevp(i,k)+psevp(i,k))
t(i,k) = t(i,k)-xlwork2/cpm(i,k)*dtcld
endif
enddo
enddo
!
! Inline expansion for fpvs
! qs(i,k,1) = fpvs(t(i,k),0,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
! qs(i,k,2) = fpvs(t(i,k),1,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
hsub = xls
hvap = xlv0
cvap = cpv
ttp=t0c+0.01
dldt=cvap-cliq
xa=-dldt/rv
xb=xa+hvap/(rv*ttp)
dldti=cvap-cice
xai=-dldti/rv
xbi=xai+hsub/(rv*ttp)
do k = kts, kte
do i = its, ite
tr=ttp/t(i,k)
logtr = log(tr)
qs(i,k,1)=psat*exp(logtr*(xa)+xb*(1.-tr))
qs(i,k,1) = min(qs(i,k,1),0.99*p(i,k))
qs(i,k,1) = ep2 * qs(i,k,1) / (p(i,k) - qs(i,k,1))
qs(i,k,1) = max(qs(i,k,1),qmin)
enddo
enddo
!
!----------------------------------------------------------------
! pcond: condensational/evaporational rate of cloud water [HL A46] [RH83 A6]
! if there exists additional water vapor condensated/if
! evaporation of cloud water is not enough to remove subsaturation
!
do k = kts, kte
do i = its, ite
! work1(i,k,1) = conden(t(i,k),q(i,k),qs(i,k,1),xl(i,k),cpm(i,k))
work1(i,k,1) = ((max(q(i,k),qmin)-(qs(i,k,1)))/(1.+(xl(i,k)) &
*(xl(i,k))/(rv*(cpm(i,k)))*(qs(i,k,1)) &
/((t(i,k))*(t(i,k)))))
work2(i,k) = qci(i,k,1)+work1(i,k,1)
pcond(i,k) = min(max(work1(i,k,1)/dtcld,0.),max(q(i,k),0.)/dtcld)
if(qci(i,k,1).gt.0..and.work1(i,k,1).lt.0.) &
pcond(i,k) = max(work1(i,k,1),-qci(i,k,1))/dtcld
q(i,k) = q(i,k)-pcond(i,k)*dtcld
qci(i,k,1) = max(qci(i,k,1)+pcond(i,k)*dtcld,0.)
t(i,k) = t(i,k)+pcond(i,k)*xl(i,k)/cpm(i,k)*dtcld
enddo
enddo
!
!
!----------------------------------------------------------------
! padding for small values
!
do k = kts, kte
do i = its, ite
if(qci(i,k,1).le.qmin) qci(i,k,1) = 0.0
if(qci(i,k,2).le.qmin) qci(i,k,2) = 0.0
enddo
enddo
enddo ! big loops
END SUBROUTINE wsm52d
! ...................................................................
!--------------------------------------------------------------------------
REAL FUNCTION fpvs(t,ice,rd,rv,cvap,cliq,cice,hvap,hsub,psat,t0c) 9
!--------------------------------------------------------------------------
IMPLICIT NONE
!--------------------------------------------------------------------------
REAL t,rd,rv,cvap,cliq,cice,hvap,hsub,psat,t0c,dldt,xa,xb,dldti, &
xai,xbi,ttp,tr
INTEGER ice
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
ttp=t0c+0.01
dldt=cvap-cliq
xa=-dldt/rv
xb=xa+hvap/(rv*ttp)
dldti=cvap-cice
xai=-dldti/rv
xbi=xai+hsub/(rv*ttp)
tr=ttp/t
if(t.lt.ttp.and.ice.eq.1) then
fpvs=psat*exp(log(tr)*(xai))*exp(xbi*(1.-tr))
else
fpvs=psat*exp(log(tr)*(xa))*exp(xb*(1.-tr))
endif
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
END FUNCTION fpvs
!------------------------------------------------------------------------------
subroutine slope_wsm5(qrs,den,denfac,t,rslope,rslopeb,rslope2,rslope3, & 3
vt,its,ite,kts,kte)
IMPLICIT NONE
INTEGER :: its,ite, jts,jte, kts,kte
REAL, DIMENSION( its:ite , kts:kte,2) :: &
qrs, &
rslope, &
rslopeb, &
rslope2, &
rslope3, &
vt
REAL, DIMENSION( its:ite , kts:kte) :: &
den, &
denfac, &
t
REAL, PARAMETER :: t0c = 273.15
REAL, DIMENSION( its:ite , kts:kte ) :: &
n0sfac
REAL :: lamdar, lamdas, x, y, z, supcol
integer :: i, j, k
!----------------------------------------------------------------
! size distributions: (x=mixing ratio, y=air density):
! valid for mixing ratio > 1.e-9 kg/kg.
lamdar(x,y)= sqrt(sqrt(pidn0r/(x*y))) ! (pidn0r/(x*y))**.25
lamdas(x,y,z)= sqrt(sqrt(pidn0s*z/(x*y))) ! (pidn0s*z/(x*y))**.25
!
do k = kts, kte
do i = its, ite
supcol = t0c-t(i,k)
!---------------------------------------------------------------
! n0s: Intercept parameter for snow [m-4] [HDC 6]
!---------------------------------------------------------------
n0sfac(i,k) = max(min(exp(alpha*supcol),n0smax/n0s),1.)
if(qrs(i,k,1).le.qcrmin)then
rslope(i,k,1) = rslopermax
rslopeb(i,k,1) = rsloperbmax
rslope2(i,k,1) = rsloper2max
rslope3(i,k,1) = rsloper3max
else
rslope(i,k,1) = 1./lamdar(qrs(i,k,1),den(i,k))
rslopeb(i,k,1) = exp(log(rslope(i,k,1))*(bvtr))
rslope2(i,k,1) = rslope(i,k,1)*rslope(i,k,1)
rslope3(i,k,1) = rslope2(i,k,1)*rslope(i,k,1)
endif
if(qrs(i,k,2).le.qcrmin)then
rslope(i,k,2) = rslopesmax
rslopeb(i,k,2) = rslopesbmax
rslope2(i,k,2) = rslopes2max
rslope3(i,k,2) = rslopes3max
else
rslope(i,k,2) = 1./lamdas(qrs(i,k,2),den(i,k),n0sfac(i,k))
rslopeb(i,k,2) = exp(log(rslope(i,k,2))*(bvts))
rslope2(i,k,2) = rslope(i,k,2)*rslope(i,k,2)
rslope3(i,k,2) = rslope2(i,k,2)*rslope(i,k,2)
endif
vt(i,k,1) = pvtr*rslopeb(i,k,1)*denfac(i,k)
vt(i,k,2) = pvts*rslopeb(i,k,2)*denfac(i,k)
if(qrs(i,k,1).le.0.0) vt(i,k,1) = 0.0
if(qrs(i,k,2).le.0.0) vt(i,k,2) = 0.0
enddo
enddo
END subroutine slope_wsm5
!-----------------------------------------------------------------------------
subroutine slope_rain(qrs,den,denfac,t,rslope,rslopeb,rslope2,rslope3, & 6
vt,its,ite,kts,kte)
IMPLICIT NONE
INTEGER :: its,ite, jts,jte, kts,kte
REAL, DIMENSION( its:ite , kts:kte) :: &
qrs, &
rslope, &
rslopeb, &
rslope2, &
rslope3, &
vt, &
den, &
denfac, &
t
REAL, PARAMETER :: t0c = 273.15
REAL, DIMENSION( its:ite , kts:kte ) :: &
n0sfac
REAL :: lamdar, x, y, z, supcol
integer :: i, j, k
!----------------------------------------------------------------
! size distributions: (x=mixing ratio, y=air density):
! valid for mixing ratio > 1.e-9 kg/kg.
lamdar(x,y)= sqrt(sqrt(pidn0r/(x*y))) ! (pidn0r/(x*y))**.25
!
do k = kts, kte
do i = its, ite
if(qrs(i,k).le.qcrmin)then
rslope(i,k) = rslopermax
rslopeb(i,k) = rsloperbmax
rslope2(i,k) = rsloper2max
rslope3(i,k) = rsloper3max
else
rslope(i,k) = 1./lamdar(qrs(i,k),den(i,k))
rslopeb(i,k) = rslope(i,k)**bvtr
rslope2(i,k) = rslope(i,k)*rslope(i,k)
rslope3(i,k) = rslope2(i,k)*rslope(i,k)
endif
vt(i,k) = pvtr*rslopeb(i,k)*denfac(i,k)
if(qrs(i,k).le.0.0) vt(i,k) = 0.0
enddo
enddo
END subroutine slope_rain
!------------------------------------------------------------------------------
subroutine slope_snow(qrs,den,denfac,t,rslope,rslopeb,rslope2,rslope3, & 4
vt,its,ite,kts,kte)
IMPLICIT NONE
INTEGER :: its,ite, jts,jte, kts,kte
REAL, DIMENSION( its:ite , kts:kte) :: &
qrs, &
rslope, &
rslopeb, &
rslope2, &
rslope3, &
vt, &
den, &
denfac, &
t
REAL, PARAMETER :: t0c = 273.15
REAL, DIMENSION( its:ite , kts:kte ) :: &
n0sfac
REAL :: lamdas, x, y, z, supcol
integer :: i, j, k
!----------------------------------------------------------------
! size distributions: (x=mixing ratio, y=air density):
! valid for mixing ratio > 1.e-9 kg/kg.
lamdas(x,y,z)= sqrt(sqrt(pidn0s*z/(x*y))) ! (pidn0s*z/(x*y))**.25
!
do k = kts, kte
do i = its, ite
supcol = t0c-t(i,k)
!---------------------------------------------------------------
! n0s: Intercept parameter for snow [m-4] [HDC 6]
!---------------------------------------------------------------
n0sfac(i,k) = max(min(exp(alpha*supcol),n0smax/n0s),1.)
if(qrs(i,k).le.qcrmin)then
rslope(i,k) = rslopesmax
rslopeb(i,k) = rslopesbmax
rslope2(i,k) = rslopes2max
rslope3(i,k) = rslopes3max
else
rslope(i,k) = 1./lamdas(qrs(i,k),den(i,k),n0sfac(i,k))
rslopeb(i,k) = rslope(i,k)**bvts
rslope2(i,k) = rslope(i,k)*rslope(i,k)
rslope3(i,k) = rslope2(i,k)*rslope(i,k)
endif
vt(i,k) = pvts*rslopeb(i,k)*denfac(i,k)
if(qrs(i,k).le.0.0) vt(i,k) = 0.0
enddo
enddo
END subroutine slope_snow
!-------------------------------------------------------------------
SUBROUTINE nislfv_rain_plm(im,km,denl,denfacl,tkl,dzl,wwl,rql,precip,dt,id,iter) 9,5
!-------------------------------------------------------------------
!
! for non-iteration semi-Lagrangain forward advection for cloud
! with mass conservation and positive definite advection
! 2nd order interpolation with monotonic piecewise linear method
! this routine is under assumption of decfl < 1 for semi_Lagrangian
!
! dzl depth of model layer in meter
! wwl terminal velocity at model layer m/s
! rql cloud density*mixing ration
! precip precipitation
! dt time step
! id kind of precip: 0 test case; 1 raindrop 2: snow
! iter how many time to guess mean terminal velocity: 0 pure forward.
! 0 : use departure wind for advection
! 1 : use mean wind for advection
! > 1 : use mean wind after iter-1 iterations
!
! author: hann-ming henry juang <henry.juang@noaa.gov>
! implemented by song-you hong
!
implicit none
integer im,km,id
real dt
real dzl(im,km),wwl(im,km),rql(im,km),precip(im)
real denl(im,km),denfacl(im,km),tkl(im,km)
!
integer i,k,n,m,kk,kb,kt,iter
real tl,tl2,qql,dql,qqd
real th,th2,qqh,dqh
real zsum,qsum,dim,dip,c1,con1,fa1,fa2
real allold, allnew, zz, dzamin, cflmax, decfl
real dz(km), ww(km), qq(km), wd(km), wa(km), was(km)
real den(km), denfac(km), tk(km)
real wi(km+1), zi(km+1), za(km+1)
real qn(km), qr(km),tmp(km),tmp1(km),tmp2(km),tmp3(km)
real dza(km+1), qa(km+1), qmi(km+1), qpi(km+1)
!
precip(:) = 0.0
!
i_loop : do i=1,im
! -----------------------------------
dz(:) = dzl(i,:)
qq(:) = rql(i,:)
ww(:) = wwl(i,:)
den(:) = denl(i,:)
denfac(:) = denfacl(i,:)
tk(:) = tkl(i,:)
! skip for no precipitation for all layers
allold = 0.0
do k=1,km
allold = allold + qq(k)
enddo
if(allold.le.0.0) then
cycle i_loop
endif
!
! compute interface values
zi(1)=0.0
do k=1,km
zi(k+1) = zi(k)+dz(k)
enddo
!
! save departure wind
wd(:) = ww(:)
n=1
100 continue
! plm is 2nd order, we can use 2nd order wi or 3rd order wi
! 2nd order interpolation to get wi
wi(1) = ww(1)
wi(km+1) = ww(km)
do k=2,km
wi(k) = (ww(k)*dz(k-1)+ww(k-1)*dz(k))/(dz(k-1)+dz(k))
enddo
! 3rd order interpolation to get wi
fa1 = 9./16.
fa2 = 1./16.
wi(1) = ww(1)
wi(2) = 0.5*(ww(2)+ww(1))
do k=3,km-1
wi(k) = fa1*(ww(k)+ww(k-1))-fa2*(ww(k+1)+ww(k-2))
enddo
wi(km) = 0.5*(ww(km)+ww(km-1))
wi(km+1) = ww(km)
!
! terminate of top of raingroup
do k=2,km
if( ww(k).eq.0.0 ) wi(k)=ww(k-1)
enddo
!
! diffusivity of wi
con1 = 0.05
do k=km,1,-1
decfl = (wi(k+1)-wi(k))*dt/dz(k)
if( decfl .gt. con1 ) then
wi(k) = wi(k+1) - con1*dz(k)/dt
endif
enddo
! compute arrival point
do k=1,km+1
za(k) = zi(k) - wi(k)*dt
enddo
!
do k=1,km
dza(k) = za(k+1)-za(k)
enddo
dza(km+1) = zi(km+1) - za(km+1)
!
! computer deformation at arrival point
do k=1,km
qa(k) = qq(k)*dz(k)/dza(k)
qr(k) = qa(k)/den(k)
enddo
qa(km+1) = 0.0
! call maxmin(km,1,qa,' arrival points ')
!
! compute arrival terminal velocity, and estimate mean terminal velocity
! then back to use mean terminal velocity
if( n.le.iter ) then
if (id.eq.1) then
call slope_rain(qr,den,denfac,tk,tmp,tmp1,tmp2,tmp3,wa,1,1,1,km)
else
call slope_snow(qr,den,denfac,tk,tmp,tmp1,tmp2,tmp3,wa,1,1,1,km)
endif
if( n.ge.2 ) wa(1:km)=0.5*(wa(1:km)+was(1:km))
do k=1,km
!#ifdef DEBUG
! print*,' slope_wsm3 ',qr(k)*1000.,den(k),denfac(k),tk(k),tmp(k),tmp1(k),tmp2(k),ww(k),wa(k)
!#endif
! mean wind is average of departure and new arrival winds
ww(k) = 0.5* ( wd(k)+wa(k) )
enddo
was(:) = wa(:)
n=n+1
go to 100
endif
!
! estimate values at arrival cell interface with monotone
do k=2,km
dip=(qa(k+1)-qa(k))/(dza(k+1)+dza(k))
dim=(qa(k)-qa(k-1))/(dza(k-1)+dza(k))
if( dip*dim.le.0.0 ) then
qmi(k)=qa(k)
qpi(k)=qa(k)
else
qpi(k)=qa(k)+0.5*(dip+dim)*dza(k)
qmi(k)=2.0*qa(k)-qpi(k)
if( qpi(k).lt.0.0 .or. qmi(k).lt.0.0 ) then
qpi(k) = qa(k)
qmi(k) = qa(k)
endif
endif
enddo
qpi(1)=qa(1)
qmi(1)=qa(1)
qmi(km+1)=qa(km+1)
qpi(km+1)=qa(km+1)
!
! interpolation to regular point
qn = 0.0
kb=1
kt=1
intp : do k=1,km
kb=max(kb-1,1)
kt=max(kt-1,1)
! find kb and kt
if( zi(k).ge.za(km+1) ) then
exit intp
else
find_kb : do kk=kb,km
if( zi(k).le.za(kk+1) ) then
kb = kk
exit find_kb
else
cycle find_kb
endif
enddo find_kb
find_kt : do kk=kt,km
if( zi(k+1).le.za(kk) ) then
kt = kk
exit find_kt
else
cycle find_kt
endif
enddo find_kt
kt = kt - 1
! compute q with piecewise constant method
if( kt.eq.kb ) then
tl=(zi(k)-za(kb))/dza(kb)
th=(zi(k+1)-za(kb))/dza(kb)
tl2=tl*tl
th2=th*th
qqd=0.5*(qpi(kb)-qmi(kb))
qqh=qqd*th2+qmi(kb)*th
qql=qqd*tl2+qmi(kb)*tl
qn(k) = (qqh-qql)/(th-tl)
else if( kt.gt.kb ) then
tl=(zi(k)-za(kb))/dza(kb)
tl2=tl*tl
qqd=0.5*(qpi(kb)-qmi(kb))
qql=qqd*tl2+qmi(kb)*tl
dql = qa(kb)-qql
zsum = (1.-tl)*dza(kb)
qsum = dql*dza(kb)
if( kt-kb.gt.1 ) then
do m=kb+1,kt-1
zsum = zsum + dza(m)
qsum = qsum + qa(m) * dza(m)
enddo
endif
th=(zi(k+1)-za(kt))/dza(kt)
th2=th*th
qqd=0.5*(qpi(kt)-qmi(kt))
dqh=qqd*th2+qmi(kt)*th
zsum = zsum + th*dza(kt)
qsum = qsum + dqh*dza(kt)
qn(k) = qsum/zsum
endif
cycle intp
endif
!
enddo intp
!
! rain out
sum_precip: do k=1,km
if( za(k).lt.0.0 .and. za(k+1).lt.0.0 ) then
precip(i) = precip(i) + qa(k)*dza(k)
cycle sum_precip
else if ( za(k).lt.0.0 .and. za(k+1).ge.0.0 ) then
precip(i) = precip(i) + qa(k)*(0.0-za(k))
exit sum_precip
endif
exit sum_precip
enddo sum_precip
!
! replace the new values
rql(i,:) = qn(:)
!
! ----------------------------------
enddo i_loop
!
END SUBROUTINE nislfv_rain_plm
# else
# include "mic-wsm5-3-5-code.h"
# endif
! end of #ifndef XEON_OPTIMIZED_WSM5
! remainder of routines are common to original and MIC version
!+---+-----------------------------------------------------------------+
subroutine refl10cm_wsm5 (qv1d, qr1d, qs1d, & 2,1
t1d, p1d, dBZ, kts, kte, ii, jj)
IMPLICIT NONE
!..Sub arguments
INTEGER, INTENT(IN):: kts, kte, ii, jj
REAL, DIMENSION(kts:kte), INTENT(IN):: &
qv1d, qr1d, qs1d, t1d, p1d
REAL, DIMENSION(kts:kte), INTENT(INOUT):: dBZ
!..Local variables
REAL, DIMENSION(kts:kte):: temp, pres, qv, rho
REAL, DIMENSION(kts:kte):: rr, rs
REAL:: temp_C
DOUBLE PRECISION, DIMENSION(kts:kte):: ilamr, ilams
DOUBLE PRECISION, DIMENSION(kts:kte):: N0_r, N0_s
DOUBLE PRECISION:: lamr, lams
LOGICAL, DIMENSION(kts:kte):: L_qr, L_qs
REAL, DIMENSION(kts:kte):: ze_rain, ze_snow
DOUBLE PRECISION:: fmelt_s
INTEGER:: i, k, k_0, kbot, n
LOGICAL:: melti
DOUBLE PRECISION:: cback, x, eta, f_d
REAL, PARAMETER:: R=287.
!+---+
do k = kts, kte
dBZ(k) = -35.0
enddo
!+---+-----------------------------------------------------------------+
!..Put column of data into local arrays.
!+---+-----------------------------------------------------------------+
do k = kts, kte
temp(k) = t1d(k)
temp_C = min(-0.001, temp(K)-273.15)
qv(k) = MAX(1.E-10, qv1d(k))
pres(k) = p1d(k)
rho(k) = 0.622*pres(k)/(R*temp(k)*(qv(k)+0.622))
if (qr1d(k) .gt. 1.E-9) then
rr(k) = qr1d(k)*rho(k)
N0_r(k) = n0r
lamr = (xam_r*xcrg(3)*N0_r(k)/rr(k))**(1./xcre(1))
ilamr(k) = 1./lamr
L_qr(k) = .true.
else
rr(k) = 1.E-12
L_qr(k) = .false.
endif
if (qs1d(k) .gt. 1.E-9) then
rs(k) = qs1d(k)*rho(k)
N0_s(k) = min(n0smax, n0s*exp(-alpha*temp_C))
lams = (xam_s*xcsg(3)*N0_s(k)/rs(k))**(1./xcse(1))
ilams(k) = 1./lams
L_qs(k) = .true.
else
rs(k) = 1.E-12
L_qs(k) = .false.
endif
enddo
!+---+-----------------------------------------------------------------+
!..Locate K-level of start of melting (k_0 is level above).
!+---+-----------------------------------------------------------------+
melti = .false.
k_0 = kts
do k = kte-1, kts, -1
if ( (temp(k).gt.273.15) .and. L_qr(k) .and. L_qs(k+1) ) then
k_0 = MAX(k+1, k_0)
melti=.true.
goto 195
endif
enddo
195 continue
!+---+-----------------------------------------------------------------+
!..Assume Rayleigh approximation at 10 cm wavelength. Rain (all temps)
!.. and non-water-coated snow and graupel when below freezing are
!.. simple. Integrations of m(D)*m(D)*N(D)*dD.
!+---+-----------------------------------------------------------------+
do k = kts, kte
ze_rain(k) = 1.e-22
ze_snow(k) = 1.e-22
if (L_qr(k)) ze_rain(k) = N0_r(k)*xcrg(4)*ilamr(k)**xcre(4)
if (L_qs(k)) ze_snow(k) = (0.176/0.93) * (6.0/PI)*(6.0/PI) &
* (xam_s/900.0)*(xam_s/900.0) &
* N0_s(k)*xcsg(4)*ilams(k)**xcse(4)
enddo
!+---+-----------------------------------------------------------------+
!..Special case of melting ice (snow/graupel) particles. Assume the
!.. ice is surrounded by the liquid water. Fraction of meltwater is
!.. extremely simple based on amount found above the melting level.
!.. Uses code from Uli Blahak (rayleigh_soak_wetgraupel and supporting
!.. routines).
!+---+-----------------------------------------------------------------+
if (melti .and. k_0.ge.kts+1) then
do k = k_0-1, kts, -1
!..Reflectivity contributed by melting snow
if (L_qs(k) .and. L_qs(k_0) ) then
fmelt_s = MAX(0.005d0, MIN(1.0d0-rs(k)/rs(k_0), 0.99d0))
eta = 0.d0
lams = 1./ilams(k)
do n = 1, nrbins
x = xam_s * xxDs(n)**xbm_s
call rayleigh_soak_wetgraupel (x,DBLE(xocms),DBLE(xobms), &
fmelt_s, melt_outside_s, m_w_0, m_i_0, lamda_radar, &
CBACK, mixingrulestring_s, matrixstring_s, &
inclusionstring_s, hoststring_s, &
hostmatrixstring_s, hostinclusionstring_s)
f_d = N0_s(k)*xxDs(n)**xmu_s * DEXP(-lams*xxDs(n))
eta = eta + f_d * CBACK * simpson(n) * xdts(n)
enddo
ze_snow(k) = SNGL(lamda4 / (pi5 * K_w) * eta)
endif
enddo
endif
do k = kte, kts, -1
dBZ(k) = 10.*log10((ze_rain(k)+ze_snow(k))*1.d18)
enddo
end subroutine refl10cm_wsm5
!+---+-----------------------------------------------------------------+
!-------------------------------------------------------------------
SUBROUTINE wsm5init(den0,denr,dens,cl,cpv,allowed_to_read) 1,17
!-------------------------------------------------------------------
IMPLICIT NONE
!-------------------------------------------------------------------
!.... constants which may not be tunable
REAL, INTENT(IN) :: den0,denr,dens,cl,cpv
LOGICAL, INTENT(IN) :: allowed_to_read
!
pi = 4.*atan(1.)
xlv1 = cl-cpv
!
qc0 = 4./3.*pi*denr*r0**3*xncr/den0 ! 0.419e-3 -- .61e-3
qck1 = .104*9.8*peaut/(xncr*denr)**(1./3.)/xmyu*den0**(4./3.) ! 7.03
!
bvtr1 = 1.+bvtr
bvtr2 = 2.5+.5*bvtr
bvtr3 = 3.+bvtr
bvtr4 = 4.+bvtr
g1pbr = rgmma(bvtr1)
g3pbr = rgmma(bvtr3)
g4pbr = rgmma(bvtr4) ! 17.837825
g5pbro2 = rgmma(bvtr2) ! 1.8273
pvtr = avtr*g4pbr/6.
eacrr = 1.0
pacrr = pi*n0r*avtr*g3pbr*.25*eacrr
precr1 = 2.*pi*n0r*.78
precr2 = 2.*pi*n0r*.31*avtr**.5*g5pbro2
xmmax = (dimax/dicon)**2
roqimax = 2.08e22*dimax**8
!
bvts1 = 1.+bvts
bvts2 = 2.5+.5*bvts
bvts3 = 3.+bvts
bvts4 = 4.+bvts
g1pbs = rgmma(bvts1) !.8875
g3pbs = rgmma(bvts3)
g4pbs = rgmma(bvts4) ! 12.0786
g5pbso2 = rgmma(bvts2)
pvts = avts*g4pbs/6.
pacrs = pi*n0s*avts*g3pbs*.25
precs1 = 4.*n0s*.65
precs2 = 4.*n0s*.44*avts**.5*g5pbso2
pidn0r = pi*denr*n0r
pidn0s = pi*dens*n0s
pacrc = pi*n0s*avts*g3pbs*.25*eacrc
!
rslopermax = 1./lamdarmax
rslopesmax = 1./lamdasmax
rsloperbmax = rslopermax ** bvtr
rslopesbmax = rslopesmax ** bvts
rsloper2max = rslopermax * rslopermax
rslopes2max = rslopesmax * rslopesmax
rsloper3max = rsloper2max * rslopermax
rslopes3max = rslopes2max * rslopesmax
!
!+---+-----------------------------------------------------------------+
!..Set these variables needed for computing radar reflectivity. These
!.. get used within radar_init to create other variables used in the
!.. radar module.
xam_r = PI*denr/6.
xbm_r = 3.
xmu_r = 0.
xam_s = PI*dens/6.
xbm_s = 3.
xmu_s = 0.
xam_g = PI*dens/6. ! These 3 variables for graupel are set but unused.
xbm_g = 3.
xmu_g = 0.
call radar_init
!+---+-----------------------------------------------------------------+
END SUBROUTINE wsm5init
! ...................................................................
REAL FUNCTION rgmma(x) 74
!-------------------------------------------------------------------
IMPLICIT NONE
!-------------------------------------------------------------------
! rgmma function: use infinite product form
REAL :: euler
PARAMETER (euler=0.577215664901532)
REAL :: x, y
INTEGER :: i
if(x.eq.1.)then
rgmma=0.
else
rgmma=x*exp(euler*x)
do i=1,10000
y=float(i)
rgmma=rgmma*(1.000+x/y)*exp(-x/y)
enddo
rgmma=1./rgmma
endif
END FUNCTION rgmma
END MODULE module_mp_wsm5
#endif