!WRF:MODEL_LAYER:PHYSICS
!
MODULE module_mp_lin 1
USE module_wrf_error
!
REAL , PARAMETER, PRIVATE :: RH = 1.0
! REAL , PARAMETER, PRIVATE :: episp0 = 0.622*611.21
REAL , PARAMETER, PRIVATE :: xnor = 8.0e6
REAL , PARAMETER, PRIVATE :: xnos = 3.0e6
! Lin
! REAL , PARAMETER, PRIVATE :: xnog = 4.0e4
! REAL , PARAMETER, PRIVATE :: rhograul = 917.
! Hobbs
REAL , PARAMETER, PRIVATE :: xnog = 4.0e6
REAL , PARAMETER, PRIVATE :: rhograul = 400.
!
REAL , PARAMETER, PRIVATE :: &
qi0 = 1.0e-3, ql0 = 7.0e-4, qs0 = 6.0E-4, &
xmi50 = 4.8e-10, xmi40 = 2.46e-10, &
constb = 0.8, constd = 0.25, &
o6 = 1./6., cdrag = 0.6, &
avisc = 1.49628e-6, adiffwv = 8.7602e-5, &
axka = 1.4132e3, di50 = 1.0e-4, xmi = 4.19e-13, &
cw = 4.187e3, vf1s = 0.78, vf2s = 0.31, &
xni0 = 1.0e-2, xmnin = 1.05e-18, bni = 0.5, &
ci = 2.093e3
CONTAINS
!-------------------------------------------------------------------
! Lin et al., 1983, JAM, 1065-1092, and
! Rutledge and Hobbs, 1984, JAS, 2949-2972
!-------------------------------------------------------------------
SUBROUTINE lin_et_al(th, qv, ql, qr, qi, qs, qg, & 1,2
qrold, qsold, qgold, &
rho, pii, p, RAINNC, RAINNCV, &
dt_in, z, ht, dz8w, &
grav, cp, Rair, rvapor, &
XLS, XLV, XLF, rhowater, rhosnow, &
EP2,SVP1,SVP2,SVP3,SVPT0, &
P_QI,P_QS,P_QG, &
P_FIRST_SCALAR, &
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
!-------------------------------------------------------------------
!
! Shuhua 12/17/99
!
INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde , &
ims,ime, jms,jme, kms,kme , &
its,ite, jts,jte, kts,kte , &
P_QI,P_QS,P_QG, &
P_FIRST_SCALAR
REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
INTENT(INOUT) :: &
th, &
qv, &
ql, &
qr, &
qi, &
qs, &
qg
!
REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
INTENT(IN ) :: &
qrold, &
qsold, &
qgold, &
rho, &
pii, &
p, &
dz8w
REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
INTENT(IN ) :: z
REAL, DIMENSION( ims:ime , jms:jme ), &
INTENT(INOUT) :: RAINNC, &
RAINNCV
REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN) :: ht
REAL, INTENT(IN ) :: dt_in, &
grav, &
Rair, &
rvapor, &
cp, &
XLS, &
XLV, &
XLF, &
rhowater, &
rhosnow
REAL, INTENT(IN ) :: EP2,SVP1,SVP2,SVP3,SVPT0
! LOCAL VAR
INTEGER :: min_q, max_q
REAL, DIMENSION( its:ite , jts:jte ) &
:: rain, snow, graupel,ice
REAL, DIMENSION( kts:kte ) :: qvz, qlz, qrz, &
qiz, qsz, qgz, &
qrzold, qszold, &
qgzold, thz, &
tothz, rhoz, &
orhoz, sqrhoz, &
prez, zz, &
dzw
!
REAL :: dt, pptrain, pptsnow, pptgraul, rhoe_s, &
gindex, pptice
INTEGER :: i,j,k
!
dt=dt_in
rhoe_s=1.29
!
IF (P_QI .lt. P_FIRST_SCALAR .or. P_QS .lt. P_FIRST_SCALAR) THEN
CALL wrf_error_fatal ( 'module_mp_lin: Improper use of Lin et al scheme; no ice phase. Please chose another one.')
ENDIF
gindex=1.0
IF (P_QG .lt. P_FIRST_SCALAR) gindex=0.
j_loop: DO j = jts, jte
i_loop: DO i = its, ite
!
!- write data from 3-D to 1-D
!
DO k = kts, kte
qvz(k)=qv(i,k,j)
qlz(k)=ql(i,k,j)
qrz(k)=qr(i,k,j)
qrzold(k)=qrold(i,k,j)
thz(k)=th(i,k,j)
!
rhoz(k)=rho(i,k,j)
orhoz(k)=1./rhoz(k)
prez(k)=p(i,k,j)
sqrhoz(k)=sqrt(rhoe_s*orhoz(k))
tothz(k)=pii(i,k,j)
zz(k)=z(i,k,j)
dzw(k)=dz8w(i,k,j)
END DO
IF (P_QI .ge. P_FIRST_SCALAR) THEN
DO k = kts, kte
qiz(k)=qi(i,k,j)
ENDDO
ELSE
DO k = kts, kte
qiz(k)=0.
ENDDO
ENDIF
IF (P_QS .ge. P_FIRST_SCALAR) THEN
DO k = kts, kte
qsz(k)=qs(i,k,j)
qszold(k)=qsold(i,k,j)
ENDDO
ELSE
DO k = kts, kte
qsz(k)=0.
qszold(k)=0.
ENDDO
ENDIF
IF (P_QG .ge. P_FIRST_SCALAR) THEN
DO k = kts, kte
qgz(k)=qg(i,k,j)
qgzold(k)=qgold(i,k,j)
ENDDO
ELSE
DO k = kts, kte
qgz(k)=0.
qgzold(k)=0.
ENDDO
ENDIF
!
pptrain=0.
pptsnow=0.
pptgraul=0.
pptice=0.
CALL clphy1d( dt, qvz, qlz, qrz, qiz, qsz, qgz, &
qrzold, qszold, qgzold, &
thz, tothz, rhoz, orhoz, sqrhoz, &
prez, zz, dzw, ht(I,J), &
pptrain, pptsnow, pptgraul, pptice, &
grav, cp, Rair, rvapor, gindex, &
XLS, XLV, XLF, rhowater, rhosnow, &
EP2,SVP1,SVP2,SVP3,SVPT0, &
kts, kte, i, j )
!
! Precipitation from cloud microphysics -- only for one time step
!
! unit is transferred from m to mm
!
rain(i,j)=pptrain
snow(i,j)=pptsnow
graupel(i,j)=pptgraul
ice(i,j)=pptice
!
RAINNCV(i,j)= pptrain + pptsnow + pptgraul + pptice
RAINNC(i,j)=RAINNC(i,j) + pptrain + pptsnow + pptgraul + pptice
!
!- update data from 1-D back to 3-D
!
!
DO k = kts, kte
qv(i,k,j)=qvz(k)
ql(i,k,j)=qlz(k)
qr(i,k,j)=qrz(k)
th(i,k,j)=thz(k)
END DO
!
IF (P_QI .ge. P_FIRST_SCALAR) THEN
DO k = kts, kte
qi(i,k,j)=qiz(k)
ENDDO
ENDIF
IF (P_QS .ge. P_FIRST_SCALAR) THEN
DO k = kts, kte
qs(i,k,j)=qsz(k)
ENDDO
ENDIF
IF (P_QG .ge. P_FIRST_SCALAR) THEN
DO k = kts, kte
qg(i,k,j)=qgz(k)
ENDDO
ENDIF
!
ENDDO i_loop
ENDDO j_loop
END SUBROUTINE lin_et_al
!-----------------------------------------------------------------------
SUBROUTINE clphy1d(dt, qvz, qlz, qrz, qiz, qsz, qgz, & 1,16
qrzold, qszold, qgzold, &
thz, tothz, rho, orho, sqrho, &
prez, zz, dzw, zsfc, pptrain, pptsnow, pptgraul, &
pptice, grav, cp, Rair, rvapor, gindex, &
XLS, XLV, XLF, rhowater, rhosnow, &
EP2,SVP1,SVP2,SVP3,SVPT0, &
kts, kte, i, j )
!-----------------------------------------------------------------------
IMPLICIT NONE
!-----------------------------------------------------------------------
! This program handles the vertical 1-D cloud micphysics
!-----------------------------------------------------------------------
! avisc: constant in empirical formular for dynamic viscosity of air
! =1.49628e-6 [kg/m/s] = 1.49628e-5 [g/cm/s]
! adiffwv: constant in empirical formular for diffusivity of water
! vapor in air
! = 8.7602e-5 [kgm/s3] = 8.7602 [gcm/s3]
! axka: constant in empirical formular for thermal conductivity of air
! = 1.4132e3 [m2/s2/K] = 1.4132e7 [cm2/s2/K]
! qi0: mixing ratio threshold for cloud ice aggregation [kg/kg]
! xmi50: mass of a 50 micron ice crystal
! = 4.8e-10 [kg] =4.8e-7 [g]
! xmi40: mass of a 40 micron ice crystal
! = 2.46e-10 [kg] = 2.46e-7 [g]
! di50: diameter of a 50 micro (radius) ice crystal
! =1.0e-4 [m]
! xmi: mass of one cloud ice crystal
! =4.19e-13 [kg] = 4.19e-10 [g]
! oxmi=1.0/xmi
!
! xni0=1.0e-2 [m-3] The value given in Lin et al. is wrong.(see
! Hsie et al.(1980) and Rutledge and Hobbs(1983) )
! bni=0.5 [K-1]
! xmnin: mass of a natural ice nucleus
! = 1.05e-18 [kg] = 1.05e-15 [g] This values is suggested by
! Hsie et al. (1980)
! = 1.0e-12 [kg] suggested by Rutlegde and Hobbs (1983)
! rhowater: density of water=1.0 g/cm3=1000.0 kg/m3
! consta: constant in empirical formular for terminal
! velocity of raindrop
! =2115.0 [cm**(1-b)/s] = 2115.0*0.01**(1-b) [m**(1-b)/s]
! constb: constant in empirical formular for terminal
! velocity of raindrop
! =0.8
! xnor: intercept parameter of the raindrop size distribution
! = 0.08 cm-4 = 8.0e6 m-4
! araut: time sacle for autoconversion of cloud water to raindrops
! =1.0e-3 [s-1]
! ql0: mixing ratio threshold for cloud watercoalescence [kg/kg]
! vf1r: ventilation factors for rain =0.78
! vf2r: ventilation factors for rain =0.31
! rhosnow: density of snow=0.1 g/cm3=100.0 kg/m3
! constc: constant in empirical formular for terminal
! velocity of snow
! =152.93 [cm**(1-d)/s] = 152.93*0.01**(1-d) [m**(1-d)/s]
! constd: constant in empirical formular for terminal
! velocity of snow
! =0.25
! xnos: intercept parameter of the snowflake size distribution
! vf1s: ventilation factors for snow =0.78
! vf2s: ventilation factors for snow =0.31
!
!----------------------------------------------------------------------
INTEGER, INTENT(IN ) :: kts, kte, i, j
REAL, DIMENSION( kts:kte ), &
INTENT(INOUT) :: qvz, qlz, qrz, qiz, qsz, &
qgz, thz
REAL, DIMENSION( kts:kte ), &
INTENT(IN ) :: tothz, rho, orho, sqrho, &
prez, zz, qrzold, qszold, &
qgzold, dzw
REAL, INTENT(IN ) :: dt, grav, cp, Rair, rvapor, &
XLS, XLV, XLF, rhowater, &
rhosnow,EP2,SVP1,SVP2,SVP3,SVPT0
REAL, INTENT(INOUT) :: pptrain, pptsnow, pptgraul, pptice
REAL, INTENT(IN ) :: zsfc
! local vars
REAL :: obp4, bp3, bp5, bp6, odp4, &
dp3, dp5, dp5o2
! temperary vars
REAL :: tmp, tmp0, tmp1, tmp2,tmp3, &
tmp4,delta2,delta3, delta4, &
tmpa,tmpb,tmpc,tmpd,alpha1, &
qic, abi,abr, abg, odtberg, &
vti50,eiw,eri,esi,esr, esw, &
erw,delrs,term0,term1,araut, &
constg2, vf1r, vf2r,alpha2, &
Ap, Bp, egw, egi, egs, egr, &
constg, gdelta4, g1sdelt4, &
factor, tmp_r, tmp_s,tmp_g, &
qlpqi, rsat, a1, a2, xnin
INTEGER :: k
!
REAL, DIMENSION( kts:kte ) :: oprez, tem, temcc, theiz, qswz, &
qsiz, qvoqswz, qvoqsiz, qvzodt, &
qlzodt, qizodt, qszodt, qrzodt, &
qgzodt
REAL, DIMENSION( kts:kte ) :: psnow, psaut, psfw, psfi, praci, &
piacr, psaci, psacw, psdep, pssub, &
pracs, psacr, psmlt, psmltevp, &
prain, praut, pracw, prevp, pvapor, &
pclw, pladj, pcli, pimlt, pihom, &
pidw, piadj, pgraupel, pgaut, &
pgfr, pgacw, pgaci, pgacr, pgacs, &
pgacip,pgacrp,pgacsp,pgwet, pdry, &
pgsub, pgdep, pgmlt, pgmltevp, &
qschg, qgchg
!
REAL, DIMENSION( kts:kte ) :: qvsbar, rs0, viscmu, visc, diffwv, &
schmidt, xka
REAL, DIMENSION( kts:kte ) :: vtr, vts, vtg, &
vtrold, vtsold, vtgold, vtiold, &
xlambdar, xlambdas, xlambdag, &
olambdar, olambdas, olambdag
REAL :: episp0k, dtb, odtb, pi, pio4, &
pio6, oxLf, xLvocp, xLfocp, consta, &
constc, ocdrag, gambp4, gamdp4, &
gam4pt5, Cpor, oxmi, gambp3, gamdp3,&
gambp6, gam3pt5, gam2pt75, gambp5o2,&
gamdp5o2, cwoxlf, ocp, xni50, es
!
REAL :: qvmin=1.e-20
REAL :: gindex
REAL :: temc1,save1,save2,xni50mx
! for terminal velocity flux
INTEGER :: min_q, max_q
REAL :: t_del_tv, del_tv, flux, fluxin, fluxout ,tmpqrz
LOGICAL :: notlast
!
dtb=dt
odtb=1./dtb
pi=acos(-1.)
pio4=acos(-1.)/4.
pio6=acos(-1.)/6.
ocp=1./cp
oxLf=1./xLf
xLvocp=xLv/cp
xLfocp=xLf/cp
consta=2115.0*0.01**(1-constb)
constc=152.93*0.01**(1-constd)
ocdrag=1./Cdrag
! episp0k=RH*episp0
episp0k=RH*ep2*1000.*svp1
!
gambp4=ggamma(constb+4.)
gamdp4=ggamma(constd+4.)
gam4pt5=ggamma(4.5)
Cpor=cp/Rair
oxmi=1.0/xmi
gambp3=ggamma(constb+3.)
gamdp3=ggamma(constd+3.)
gambp6=ggamma(constb+6)
gam3pt5=ggamma(3.5)
gam2pt75=ggamma(2.75)
gambp5o2=ggamma((constb+5.)/2.)
gamdp5o2=ggamma((constd+5.)/2.)
cwoxlf=cw/xlf
delta2=0.
delta3=0.
delta4=0.
!
!-----------------------------------------------------------------------
! oprez 1./prez ( prez : pressure)
! qsw saturated mixing ratio on water surface
! qsi saturated mixing ratio on ice surface
! episp0k RH*e*saturated pressure at 273.15 K
! qvoqsw qv/qsw
! qvoqsi qv/qsi
! qvzodt qv/dt
! qlzodt ql/dt
! qizodt qi/dt
! qszodt qs/dt
! qrzodt qr/dt
! qgzodt qg/dt
!
! temcc temperature in dregee C
!
obp4=1.0/(constb+4.0)
bp3=constb+3.0
bp5=constb+5.0
bp6=constb+6.0
odp4=1.0/(constd+4.0)
dp3=constd+3.0
dp5=constd+5.0
dp5o2=0.5*(constd+5.0)
!
do k=kts,kte
oprez(k)=1./prez(k)
enddo
do k=kts,kte
qlz(k)=amax1( 0.0,qlz(k) )
qiz(k)=amax1( 0.0,qiz(k) )
qvz(k)=amax1( qvmin,qvz(k) )
qsz(k)=amax1( 0.0,qsz(k) )
qrz(k)=amax1( 0.0,qrz(k) )
qgz(k)=amax1( 0.0,qgz(k) )
!
tem(k)=thz(k)*tothz(k)
temcc(k)=tem(k)-273.15
!
! qswz(k)=episp0k*oprez(k)* &
! exp( svp2*temcc(k)/(tem(k)-svp3) )
es=1000.*svp1*exp( svp2*temcc(k)/(tem(k)-svp3) )
qswz(k)=ep2*es/(prez(k)-es)
if (tem(k) .lt. 273.15 ) then
! qsiz(k)=episp0k*oprez(k)* &
! exp( 21.8745584*(tem(k)-273.16)/(tem(k)-7.66) )
es=1000.*svp1*exp( 21.8745584*(tem(k)-273.16)/(tem(k)-7.66) )
qsiz(k)=ep2*es/(prez(k)-es)
if (temcc(k) .lt. -40.0) qswz(k)=qsiz(k)
else
qsiz(k)=qswz(k)
endif
!
qvoqswz(k)=qvz(k)/qswz(k)
qvoqsiz(k)=qvz(k)/qsiz(k)
qvzodt(k)=amax1( 0.0,odtb*qvz(k) )
qlzodt(k)=amax1( 0.0,odtb*qlz(k) )
qizodt(k)=amax1( 0.0,odtb*qiz(k) )
qszodt(k)=amax1( 0.0,odtb*qsz(k) )
qrzodt(k)=amax1( 0.0,odtb*qrz(k) )
qgzodt(k)=amax1( 0.0,odtb*qgz(k) )
theiz(k)=thz(k)+(xlvocp*qvz(k)-xlfocp*qiz(k))/tothz(k)
enddo
!
!
!-----------------------------------------------------------------------
! In this simple stable cloud parameterization scheme, only five
! forms of water substance (water vapor, cloud water, cloud ice,
! rain and snow are considered. The prognostic variables are total
! water (qp),cloud water (ql), and cloud ice (qi). Rain and snow are
! diagnosed following Nagata and Ogura, 1991, MWR, 1309-1337. Eq (A7).
! the micro physics are based on (1) Hsie et al.,1980, JAM, 950-977 ;
! (2) Lin et al., 1983, JAM, 1065-1092 ; (3) Rutledge and Hobbs, 1983,
! JAS, 1185-1206 ; (4) Rutledge and Hobbs, 1984, JAS, 2949-2972.
!-----------------------------------------------------------------------
!
! rhowater: density of water=1.0 g/cm3=1000.0 kg/m3
! rhosnow: density of snow=0.1 g/cm3=100.0 kg/m3
! xnor: intercept parameter of the raindrop size distribution
! = 0.08 cm-4 = 8.0e6 m-4
! xnos: intercept parameter of the snowflake size distribution
! = 0.03 cm-4 = 3.0e6 m-4
! xnog: intercept parameter of the graupel size distribution
! = 4.0e-4 cm-4 = 4.0e4 m-4
! consta: constant in empirical formular for terminal
! velocity of raindrop
! =2115.0 [cm**(1-b)/s] = 2115.0*0.01**(1-b) [m**(1-b)/s]
! constb: constant in empirical formular for terminal
! velocity of raindrop
! =0.8
! constc: constant in empirical formular for terminal
! velocity of snow
! =152.93 [cm**(1-d)/s] = 152.93*0.01**(1-d) [m**(1-d)/s]
! constd: constant in empirical formular for terminal
! velocity of snow
! =0.25
! avisc: constant in empirical formular for dynamic viscosity of air
! =1.49628e-6 [kg/m/s] = 1.49628e-5 [g/cm/s]
! adiffwv: constant in empirical formular for diffusivity of water
! vapor in air
! = 8.7602e-5 [kgm/s3] = 8.7602 [gcm/s3]
! axka: constant in empirical formular for thermal conductivity of air
! = 1.4132e3 [m2/s2/K] = 1.4132e7 [cm2/s2/K]
! qi0: mixing ratio threshold for cloud ice aggregation [kg/kg]
! = 1.0e-3 g/g = 1.0e-3 kg/gk
! ql0: mixing ratio threshold for cloud watercoalescence [kg/kg]
! = 2.0e-3 g/g = 2.0e-3 kg/gk
! qs0: mixing ratio threshold for snow aggregation
! = 6.0e-4 g/g = 6.0e-4 kg/gk
! xmi50: mass of a 50 micron ice crystal
! = 4.8e-10 [kg] =4.8e-7 [g]
! xmi40: mass of a 40 micron ice crystal
! = 2.46e-10 [kg] = 2.46e-7 [g]
! di50: diameter of a 50 micro (radius) ice crystal
! =1.0e-4 [m]
! xmi: mass of one cloud ice crystal
! =4.19e-13 [kg] = 4.19e-10 [g]
! oxmi=1.0/xmi
!
! if gindex=1.0 include graupel
! if gindex=0. no graupel
!
!
do k=kts,kte
psnow(k)=0.0
psaut(k)=0.0
psfw(k)=0.0
psfi(k)=0.0
praci(k)=0.0
piacr(k)=0.0
psaci(k)=0.0
psacw(k)=0.0
psdep(k)=0.0
pssub(k)=0.0
pracs(k)=0.0
psacr(k)=0.0
psmlt(k)=0.0
psmltevp(k)=0.0
!
prain(k)=0.0
praut(k)=0.0
pracw(k)=0.0
prevp(k)=0.0
!
pvapor(k)=0.0
!
pclw(k)=0.0
pladj(k)=0.0
!
pcli(k)=0.0
pimlt(k)=0.0
pihom(k)=0.0
pidw(k)=0.0
piadj(k)=0.0
enddo
!
!!! graupel
!
do k=kts,kte
pgraupel(k)=0.0
pgaut(k)=0.0
pgfr(k)=0.0
pgacw(k)=0.0
pgaci(k)=0.0
pgacr(k)=0.0
pgacs(k)=0.0
pgacip(k)=0.0
pgacrP(k)=0.0
pgacsp(k)=0.0
pgwet(k)=0.0
pdry(k)=0.0
pgsub(k)=0.0
pgdep(k)=0.0
pgmlt(k)=0.0
pgmltevp(k)=0.0
qschg(k)=0.
qgchg(k)=0.
enddo
!
!
! Set rs0=episp0*oprez(k)
! episp0=e*saturated pressure at 273.15 K
! e = 0.622
!
DO k=kts,kte
rs0(k)=ep2*1000.*svp1/(prez(k)-1000.*svp1)
END DO
!
!***********************************************************************
! Calculate precipitation fluxes due to terminal velocities.
!***********************************************************************
!
!- Calculate termianl velocity (vt?) of precipitation q?z
!- Find maximum vt? to determine the small delta t
!
!-- rain
!
t_del_tv=0.
del_tv=dtb
notlast=.true.
DO while (notlast)
!
min_q=kte
max_q=kts-1
!
do k=kts,kte-1
! if (qrzold(k) .gt. 1.0e-12) then
!! if (qrzold(k) .gt. 1.0e-8) then
if (qrz(k) .gt. 1.0e-8) then
min_q=min0(min_q,k)
max_q=max0(max_q,k)
!! tmp1=sqrt(pi*rhowater*xnor/rho(k)/qrzold(k))
tmp1=sqrt(pi*rhowater*xnor/rho(k)/qrz(k))
tmp1=sqrt(tmp1)
vtrold(k)=o6*consta*gambp4*sqrho(k)/tmp1**constb
if (k .eq. 1) then
del_tv=amin1(del_tv,0.9*(zz(k)-zsfc)/vtrold(k))
else
del_tv=amin1(del_tv,0.9*(zz(k)-zz(k-1))/vtrold(k))
endif
else
vtrold(k)=0.
endif
enddo
if (max_q .ge. min_q) then
!
!- Check if the summation of the small delta t >= big delta t
! (t_del_tv) (del_tv) (dtb)
t_del_tv=t_del_tv+del_tv
!
if ( t_del_tv .ge. dtb ) then
notlast=.false.
del_tv=dtb+del_tv-t_del_tv
endif
! use small delta t to calculate the q?z flux
! termi is the qrzold flux pass in the grid box through the upper boundary
! termo is the qrzold flux pass out the grid box through the lower boundary
!
fluxin=0.
do k=max_q,min_q,-1
! fluxout=rho(k)*vtrold(k)*qrzold(k)
fluxout=rho(k)*vtrold(k)*qrz(k)
flux=(fluxin-fluxout)/rho(k)/dzw(k)
tmpqrz=qrz(k)
qrz(k)=qrz(k)+del_tv*flux
fluxin=fluxout
enddo
if (min_q .eq. 1) then
pptrain=pptrain+fluxin*del_tv
else
qrz(min_q-1)=qrz(min_q-1)+del_tv* &
fluxin/rho(min_q-1)/dzw(min_q-1)
endif
!
else
notlast=.false.
endif
ENDDO
!
!-- snow
!
t_del_tv=0.
del_tv=dtb
notlast=.true.
DO while (notlast)
!
min_q=kte
max_q=kts-1
!
do k=kts,kte-1
! if (qszold(k) .gt. 1.0e-12) then
!! if (qszold(k) .gt. 1.0e-8) then
if (qsz(k) .gt. 1.0e-8) then
min_q=min0(min_q,k)
max_q=max0(max_q,k)
!! tmp1=sqrt(pi*rhosnow*xnos/rho(k)/qszold(k))
tmp1=sqrt(pi*rhosnow*xnos/rho(k)/qsz(k))
tmp1=sqrt(tmp1)
vtsold(k)=o6*constc*gamdp4*sqrho(k)/tmp1**constd
if (k .eq. 1) then
del_tv=amin1(del_tv,0.9*(zz(k)-zsfc)/vtsold(k))
else
del_tv=amin1(del_tv,0.9*(zz(k)-zz(k-1))/vtsold(k))
endif
else
vtsold(k)=0.
endif
enddo
if (max_q .ge. min_q) then
!
!
!- Check if the summation of the small delta t >= big delta t
! (t_del_tv) (del_tv) (dtb)
t_del_tv=t_del_tv+del_tv
if ( t_del_tv .ge. dtb ) then
notlast=.false.
del_tv=dtb+del_tv-t_del_tv
endif
! use small delta t to calculate the qszold flux
! termi is the qszold flux pass in the grid box through the upper boundary
! termo is the qszold flux pass out the grid box through the lower boundary
!
fluxin=0.
do k=max_q,min_q,-1
!! fluxout=rho(k)*vtsold(k)*qszold(k)
fluxout=rho(k)*vtsold(k)*qsz(k)
flux=(fluxin-fluxout)/rho(k)/dzw(k)
qsz(k)=qsz(k)+del_tv*flux
qsz(k)=amax1(0.,qsz(k))
fluxin=fluxout
enddo
if (min_q .eq. 1) then
pptsnow=pptsnow+fluxin*del_tv
else
qsz(min_q-1)=qsz(min_q-1)+del_tv* &
fluxin/rho(min_q-1)/dzw(min_q-1)
endif
!
else
notlast=.false.
endif
ENDDO
!
!-- grauupel
!
t_del_tv=0.
del_tv=dtb
notlast=.true.
!
DO while (notlast)
!
min_q=kte
max_q=kts-1
!
do k=kts,kte-1
! if (qgzold(k) .gt. 1.0e-12) then
!! if (qgzold(k) .gt. 1.0e-8) then
if (qgz(k) .gt. 1.0e-8) then
min_q=min0(min_q,k)
max_q=max0(max_q,k)
!! tmp1=sqrt(pi*rhograul*xnog/rho(k)/qgzold(k))
tmp1=sqrt(pi*rhograul*xnog/rho(k)/qgz(k))
tmp1=sqrt(tmp1)
term0=sqrt(4.*grav*rhograul*0.33334/rho(k)/cdrag)
vtgold(k)=o6*gam4pt5*term0*sqrt(1./tmp1)
if (k .eq. 1) then
del_tv=amin1(del_tv,0.9*(zz(k)-zsfc)/vtgold(k))
else
del_tv=amin1(del_tv,0.9*(zz(k)-zz(k-1))/vtgold(k))
endif
else
vtgold(k)=0.
endif
enddo
if (max_q .ge. min_q) then
!
!
!- Check if the summation of the small delta t >= big delta t
! (t_del_tv) (del_tv) (dtb)
t_del_tv=t_del_tv+del_tv
if ( t_del_tv .ge. dtb ) then
notlast=.false.
del_tv=dtb+del_tv-t_del_tv
endif
! use small delta t to calculate the qgzold flux
! termi is the qgzold flux pass in the grid box through the upper boundary
! termo is the qgzold flux pass out the grid box through the lower boundary
!
fluxin=0.
do k=max_q,min_q,-1
!! fluxout=rho(k)*vtgold(k)*qgzold(k)
fluxout=rho(k)*vtgold(k)*qgz(k)
flux=(fluxin-fluxout)/rho(k)/dzw(k)
qgz(k)=qgz(k)+del_tv*flux
qgz(k)=amax1(0.,qgz(k))
fluxin=fluxout
enddo
if (min_q .eq. 1) then
pptgraul=pptgraul+fluxin*del_tv
else
qgz(min_q-1)=qgz(min_q-1)+del_tv* &
fluxin/rho(min_q-1)/dzw(min_q-1)
endif
!
else
notlast=.false.
endif
!
ENDDO
!
!-- cloud ice (03/21/02) follow Vaughan T.J. Phillips at GFDL
!
t_del_tv=0.
del_tv=dtb
notlast=.true.
!
DO while (notlast)
!
min_q=kte
max_q=kts-1
!
do k=kts,kte-1
if (qiz(k) .gt. 1.0e-8) then
min_q=min0(min_q,k)
max_q=max0(max_q,k)
vtiold(k)= 3.29 * (rho(k)* qiz(k))** 0.16 ! Heymsfield and Donner
if (k .eq. 1) then
del_tv=amin1(del_tv,0.9*(zz(k)-zsfc)/vtiold(k))
else
del_tv=amin1(del_tv,0.9*(zz(k)-zz(k-1))/vtiold(k))
endif
else
vtiold(k)=0.
endif
enddo
if (max_q .ge. min_q) then
!
!
!- Check if the summation of the small delta t >= big delta t
! (t_del_tv) (del_tv) (dtb)
t_del_tv=t_del_tv+del_tv
if ( t_del_tv .ge. dtb ) then
notlast=.false.
del_tv=dtb+del_tv-t_del_tv
endif
! use small delta t to calculate the qgzold flux
! termi is the qgzold flux pass in the grid box through the upper boundary
! termo is the qgzold flux pass out the grid box through the lower boundary
!
fluxin=0.
do k=max_q,min_q,-1
fluxout=rho(k)*vtiold(k)*qiz(k)
flux=(fluxin-fluxout)/rho(k)/dzw(k)
qiz(k)=qiz(k)+del_tv*flux
qiz(k)=amax1(0.,qiz(k))
fluxin=fluxout
enddo
if (min_q .eq. 1) then
pptice=pptice+fluxin*del_tv
else
qiz(min_q-1)=qiz(min_q-1)+del_tv* &
fluxin/rho(min_q-1)/dzw(min_q-1)
endif
!
else
notlast=.false.
endif
!
ENDDO
! Microphpysics processes
!
DO 2000 k=kts,kte
!
!***********************************************************************
!***** diagnose mixing ratios (qrz,qsz), terminal *****
!***** velocities (vtr,vts), and slope parameters in size *****
!***** distribution(xlambdar,xlambdas) of rain and snow *****
!***** follows Nagata and Ogura, 1991, MWR, 1309-1337. Eq (A7) *****
!***********************************************************************
!
!**** assuming no cloud water can exist in the top two levels due to
!**** radiation consideration
!
!! if
!! unsaturated,
!! no cloud water, rain, ice, snow and graupel
!! then
!! skip these processes and jump to line 2000
!
!
tmp=qiz(k)+qlz(k)+qsz(k)+qrz(k)+qgz(k)*gindex
if( qvz(k)+qlz(k)+qiz(k) .lt. qsiz(k) &
.and. tmp .eq. 0.0 ) go to 2000
!! calculate terminal velocity of rain
!
! if (qrzold(k) .gt. 1.0e-12) then
!! if (qrzold(k) .gt. 1.0e-8) then
if (qrz(k) .gt. 1.0e-8) then
!! tmp1=sqrt(pi*rhowater*xnor*orho(k)/qrzold(k))
tmp1=sqrt(pi*rhowater*xnor*orho(k)/qrz(k))
xlambdar(k)=sqrt(tmp1)
olambdar(k)=1.0/xlambdar(k)
vtrold(k)=o6*consta*gambp4*sqrho(k)*olambdar(k)**constb
else
vtrold(k)=0.
olambdar(k)=0.
endif
!
! if (qrz(k) .gt. 1.0e-12) then
if (qrz(k) .gt. 1.0e-8) then
tmp1=sqrt(pi*rhowater*xnor*orho(k)/qrz(k))
xlambdar(k)=sqrt(tmp1)
olambdar(k)=1.0/xlambdar(k)
vtr(k)=o6*consta*gambp4*sqrho(k)*olambdar(k)**constb
else
vtr(k)=0.
olambdar(k)=0.
endif
!
!! calculate terminal velocity of snow
!
! if (qszold(k) .gt. 1.0e-12) then
!! if (qszold(k) .gt. 1.0e-8) then
if (qsz(k) .gt. 1.0e-8) then
!! tmp1=sqrt(pi*rhosnow*xnos*orho(k)/qszold(k))
tmp1=sqrt(pi*rhosnow*xnos*orho(k)/qsz(k))
xlambdas(k)=sqrt(tmp1)
olambdas(k)=1.0/xlambdas(k)
vtsold(k)=o6*constc*gamdp4*sqrho(k)*olambdas(k)**constd
else
vtsold(k)=0.
olambdas(k)=0.
endif
!
! if (qsz(k) .gt. 1.0e-12) then
if (qsz(k) .gt. 1.0e-8) then
tmp1=sqrt(pi*rhosnow*xnos*orho(k)/qsz(k))
xlambdas(k)=sqrt(tmp1)
olambdas(k)=1.0/xlambdas(k)
vts(k)=o6*constc*gamdp4*sqrho(k)*olambdas(k)**constd
else
vts(k)=0.
olambdas(k)=0.
endif
!
!! calculate terminal velocity of graupel
!
! if (qgzold(k) .gt. 1.0e-12) then
!! if (qgzold(k) .gt. 1.0e-8) then
if (qgz(k) .gt. 1.0e-8) then
!! tmp1=sqrt( pi*rhograul*xnog*orho(k)/qgzold(k))
tmp1=sqrt( pi*rhograul*xnog*orho(k)/qgz(k))
xlambdag(k)=sqrt(tmp1)
olambdag(k)=1.0/xlambdag(k)
term0=sqrt(4.*grav*rhograul*0.33334*orho(k)*ocdrag)
vtgold(k)=o6*gam4pt5*term0*sqrt(olambdag(k))
else
vtgold(k)=0.
olambdag(k)=0.
endif
!
! if (qgz(k) .gt. 1.0e-12) then
if (qgz(k) .gt. 1.0e-8) then
tmp1=sqrt( pi*rhograul*xnog*orho(k)/qgz(k))
xlambdag(k)=sqrt(tmp1)
olambdag(k)=1.0/xlambdag(k)
term0=sqrt(4.*grav*rhograul*0.33334*orho(k)*ocdrag)
vtg(k)=o6*gam4pt5*term0*sqrt(olambdag(k))
else
vtg(k)=0.
olambdag(k)=0.
endif
!
!***********************************************************************
!***** compute viscosity,difusivity,thermal conductivity, and ******
!***** Schmidt number ******
!***********************************************************************
!c------------------------------------------------------------------
!c viscmu: dynamic viscosity of air kg/m/s
!c visc: kinematic viscosity of air = viscmu/rho (m2/s)
!c avisc=1.49628e-6 kg/m/s=1.49628e-5 g/cm/s
!c viscmu=1.718e-5 kg/m/s in RH
!c diffwv: Diffusivity of water vapor in air
!c adiffwv = 8.7602e-5 (8.794e-5 in MM5) kgm/s3
!c = 8.7602 (8.794 in MM5) gcm/s3
!c diffwv(k)=2.26e-5 m2/s
!c schmidt: Schmidt number=visc/diffwv
!c xka: thermal conductivity of air J/m/s/K (Kgm/s3/K)
!c xka(k)=2.43e-2 J/m/s/K in RH
!c axka=1.4132e3 (1.414e3 in MM5) m2/s2/k = 1.4132e7 cm2/s2/k
!c------------------------------------------------------------------
viscmu(k)=avisc*tem(k)**1.5/(tem(k)+120.0)
visc(k)=viscmu(k)*orho(k)
diffwv(k)=adiffwv*tem(k)**1.81*oprez(k)
schmidt(k)=visc(k)/diffwv(k)
xka(k)=axka*viscmu(k)
if (tem(k) .lt. 273.15) then
!
!***********************************************************************
!********* snow production processes for T < 0 C **********
!***********************************************************************
!c
!c (1) ICE CRYSTAL AGGREGATION TO SNOW (Psaut): Lin (21)
!c! psaut=alpha1*(qi-qi0)
!c! alpha1=1.0e-3*exp(0.025*(T-T0))
!c
! alpha1=1.0e-3*exp( 0.025*temcc(k) )
alpha1=1.0e-3*exp( 0.025*temcc(k) )
!
if(temcc(k) .lt. -20.0) then
tmp1=-7.6+4.0*exp( -0.2443e-3*(abs(temcc(k))-20)**2.455 )
qic=1.0e-3*exp(tmp1)*orho(k)
else
qic=qi0
end if
!testing
! tmp1=amax1( 0.0,alpha1*(qiz(k)-qic) )
! psaut(k)=amin1( tmp1,qizodt(k) )
tmp1=odtb*(qiz(k)-qic)*(1.0-exp(-alpha1*dtb))
psaut(k)=amax1( 0.0,tmp1 )
!c
!c (2) BERGERON PROCESS TRANSFER OF CLOUD WATER TO SNOW (Psfw)
!c this process only considered when -31 C < T < 0 C
!c Lin (33) and Hsie (17)
!c
!c!
!c! parama1 and parama2 functions must be user supplied
!c!
! testing
if( qlz(k) .gt. 1.0e-10 ) then
temc1=amax1(-30.99,temcc(k))
! print*,'temc1',temc1,qlz(k)
a1=parama1( temc1 )
a2=parama2( temc1 )
tmp1=1.0-a2
!! change unit from cgs to mks
a1=a1*0.001**tmp1
!c! dtberg is the time needed for a crystal to grow from 40 to 50 um
!c ! odtberg=1.0/dtberg
odtberg=(a1*tmp1)/(xmi50**tmp1-xmi40**tmp1)
!
!c! compute terminal velocity of a 50 micron ice cystal
!
vti50=constc*di50**constd*sqrho(k)
!
eiw=1.0
save1=a1*xmi50**a2
save2=0.25*pi*eiw*rho(k)*di50*di50*vti50
!
tmp2=( save1 + save2*qlz(k) )
!
!! maximum number of 50 micron crystals limited by the amount
!! of supercool water
!
xni50mx=qlzodt(k)/tmp2
!
!! number of 50 micron crystals produced
!
!
xni50=qiz(k)*( 1.0-exp(-dtb*odtberg) )/xmi50
xni50=amin1(xni50,xni50mx)
!
tmp3=odtb*tmp2/save2*( 1.0-exp(-save2*xni50*dtb) )
psfw(k)=amin1( tmp3,qlzodt(k) )
!testing
! psfw(k)=0.
!0915 if( temcc(k).gt.-30.99 ) then
!0915 a1=parama1( temcc(k) )
!0915 a2=parama2( temcc(k) )
!0915 tmp1=1.0-a2
!! change unit from cgs to mks
!0915 a1=a1*0.001**tmp1
!c! dtberg is the time needed for a crystal to grow from 40 to 50 um
!c! odtberg=1.0/dtberg
!0915 odtberg=(a1*tmp1)/(xmi50**tmp1-xmi40**tmp1)
!c! number of 50 micron crystals produced
!0915 xni50=qiz(k)*dtb*odtberg/xmi50
!c! need to calculate the terminal velocity of a 50 micron
!c! ice cystal
!0915 vti50=constc*di50**constd*sqrho(k)
!0915 eiw=1.0
!0915 tmp2=xni50*( a1*xmi50**a2 + &
!0915 0.25*qlz(k)*pi*eiw*rho(k)*di50*di50*vti50 )
!0915 psfw(k)=amin1( tmp2,qlzodt(k) )
!0915 psfw(k)=0.
!c
!c (3) REDUCTION OF CLOUD ICE BY BERGERON PROCESS (Psfi): Lin (34)
!c this process only considered when -31 C < T < 0 C
!c
tmp1=xni50*xmi50-psfw(k)
psfi(k)=amin1(tmp1,qizodt(k))
! testing
! psfi(k)=0.
end if
!
!0915 tmp1=qiz(k)*odtberg
!0915 psfi(k)=amin1(tmp1,qizodt(k))
! testing
!0915 psfi(k)=0.
!0915 end if
!
if(qrz(k) .le. 0.0) go to 1000
!
! Processes (4) and (5) only need when qrz > 0.0
!
!c
!c (4) CLOUD ICE ACCRETION BY RAIN (Praci): Lin (25)
!c may produce snow or graupel
!c
eri=1.0
!0915 tmp1=qiz(k)*pio4*eri*xnor*consta*sqrho(k)
!0915 tmp2=tmp1*gambp3*olambdar(k)**bp3
!0915 praci(k)=amin1( tmp2,qizodt(k) )
save1=pio4*eri*xnor*consta*sqrho(k)
tmp1=save1*gambp3*olambdar(k)**bp3
praci(k)=qizodt(k)*( 1.0-exp(-tmp1*dtb) )
!c
!c (5) RAIN ACCRETION BY CLOUD ICE (Piacr): Lin (26)
!c
!0915 tmp2=tmp1*rho(k)*pio6*rhowater*gambp6*oxmi* &
!0915 olambdar(k)**bp6
!0915 piacr(k)=amin1( tmp2,qrzodt(k) )
tmp2=qiz(k)*save1*rho(k)*pio6*rhowater*gambp6*oxmi* &
olambdar(k)**bp6
piacr(k)=amin1( tmp2,qrzodt(k) )
!
1000 continue
!
if(qsz(k) .le. 0.0) go to 1200
!
! Compute the following processes only when qsz > 0.0
!
!c
!c (6) ICE CRYSTAL ACCRETION BY SNOW (Psaci): Lin (22)
!c
esi=exp( 0.025*temcc(k) )
save1=pio4*xnos*constc*gamdp3*sqrho(k)* &
olambdas(k)**dp3
tmp1=esi*save1
psaci(k)=qizodt(k)*( 1.0-exp(-tmp1*dtb) )
!0915 tmp1=pio4*xnos*constc*gamdp3*sqrho(k)* &
!0915 olambdas(k)**dp3
!0915 tmp2=qiz(k)*esi*tmp1
!0915 psaci(k)=amin1( tmp2,qizodt(k) )
!c
!c (7) CLOUD WATER ACCRETION BY SNOW (Psacw): Lin (24)
!c
esw=1.0
tmp1=esw*save1
psacw(k)=qlzodt(K)*( 1.0-exp(-tmp1*dtb) )
!0915 tmp2=qlz(k)*esw*tmp1
!0915 psacw(k)=amin1( tmp2,qlzodt(k) )
!c
!c (8) DEPOSITION/SUBLIMATION OF SNOW (Psdep/Pssub): Lin (31)
!c includes consideration of ventilation effect
!c
!c abi=2*pi*(Si-1)/rho/(A"+B")
!c
tmpa=rvapor*xka(k)*tem(k)*tem(k)
tmpb=xls*xls*rho(k)*qsiz(k)*diffwv(k)
tmpc=tmpa*qsiz(k)*diffwv(k)
abi=2.0*pi*(qvoqsiz(k)-1.0)*tmpc/(tmpa+tmpb)
!
!c vf1s,vf2s=ventilation factors for snow
!c vf1s=0.78,vf2s=0.31 in LIN
!
tmp1=constc*sqrho(k)*olambdas(k)**dp5/visc(k)
tmp2=abi*xnos*( vf1s*olambdas(k)*olambdas(k)+ &
vf2s*schmidt(k)**0.33334*gamdp5o2*sqrt(tmp1) )
tmp3=odtb*( qvz(k)-qsiz(k) )
!
if( tmp2 .le. 0.0) then
tmp2=amax1( tmp2,tmp3)
pssub(k)=amax1( tmp2,-qszodt(k) )
psdep(k)=0.0
else
psdep(k)=amin1( tmp2,tmp3 )
pssub(k)=0.0
end if
!0915 psdep(k)=amax1(0.0,tmp2)
!0915 pssub(k)=amin1(0.0,tmp2)
!0915 pssub(k)=amax1( pssub(k),-qszodt(k) )
!
if(qrz(k) .le. 0.0) go to 1200
!
! Compute processes (9) and (10) only when qsz > 0.0 and qrz > 0.0
!
!c
!c (9) ACCRETION OF SNOW BY RAIN (Pracs): Lin (27)
!c
esr=1.0
tmpa=olambdar(k)*olambdar(k)
tmpb=olambdas(k)*olambdas(k)
tmpc=olambdar(k)*olambdas(k)
tmp1=pi*pi*esr*xnor*xnos*abs( vtr(k)-vts(k) )*orho(k)
tmp2=tmpb*tmpb*olambdar(k)*(5.0*tmpb+2.0*tmpc+0.5*tmpa)
tmp3=tmp1*rhosnow*tmp2
pracs(k)=amin1( tmp3,qszodt(k) )
!c
!c (10) ACCRETION OF RAIN BY SNOW (Psacr): Lin (28)
!c
tmp3=tmpa*tmpa*olambdas(k)*(5.0*tmpa+2.0*tmpc+0.5*tmpb)
tmp4=tmp1*rhowater*tmp3
psacr(k)=amin1( tmp4,qrzodt(k) )
!
1200 continue
!
else
!
!***********************************************************************
!********* snow production processes for T > 0 C **********
!***********************************************************************
!
if (qsz(k) .le. 0.0) go to 1400
!c
!c (1) CLOUD WATER ACCRETION BY SNOW (Psacw): Lin (24)
!c
esw=1.0
tmp1=esw*pio4*xnos*constc*gamdp3*sqrho(k)* &
olambdas(k)**dp3
psacw(k)=qlzodt(k)*( 1.0-exp(-tmp1*dtb) )
!0915 tmp1=pio4*xnos*constc*gamdp3*sqrho(k)* &
!0915 olambdas(k)**dp3
!0915 tmp2=qlz(k)*esw*tmp1
!0915 psacw(k)=amin1( tmp2,qlzodt(k) )
!c
!c (2) ACCRETION OF RAIN BY SNOW (Psacr): Lin (28)
!c
esr=1.0
tmpa=olambdar(k)*olambdar(k)
tmpb=olambdas(k)*olambdas(k)
tmpc=olambdar(k)*olambdas(k)
tmp1=pi*pi*esr*xnor*xnos*abs( vtr(k)-vts(k) )*orho(k)
tmp2=tmpa*tmpa*olambdas(k)*(5.0*tmpa+2.0*tmpc+0.5*tmpb)
tmp3=tmp1*rhowater*tmp2
psacr(k)=amin1( tmp3,qrzodt(k) )
!c
!c (3) MELTING OF SNOW (Psmlt): Lin (32)
!c Psmlt is negative value
!
delrs=rs0(k)-qvz(k)
term1=2.0*pi*orho(k)*( xlv*diffwv(k)*rho(k)*delrs- &
xka(k)*temcc(k) )
tmp1=constc*sqrho(k)*olambdas(k)**dp5/visc(k)
tmp2=xnos*( vf1s*olambdas(k)*olambdas(k)+ &
vf2s*schmidt(k)**0.33334*gamdp5o2*sqrt(tmp1) )
tmp3=term1*oxlf*tmp2-cwoxlf*temcc(k)*( psacw(k)+psacr(k) )
tmp4=amin1(0.0,tmp3)
psmlt(k)=amax1( tmp4,-qszodt(k) )
!c
!c (4) EVAPORATION OF MELTING SNOW (Psmltevp): HR (A27)
!c but use Lin et al. coefficience
!c Psmltevp is a negative value
!c
tmpa=rvapor*xka(k)*tem(k)*tem(k)
tmpb=xlv*xlv*rho(k)*qswz(k)*diffwv(k)
tmpc=tmpa*qswz(k)*diffwv(k)
tmpd=amin1( 0.0,(qvoqswz(k)-0.90)*qswz(k)*odtb )
! abr=2.0*pi*(qvoqswz(k)-1.0)*tmpc/(tmpa+tmpb)
abr=2.0*pi*(qvoqswz(k)-0.90)*tmpc/(tmpa+tmpb)
!
!**** allow evaporation to occur when RH less than 90%
!**** here not using 100% because the evaporation cooling
!**** of temperature is not taking into account yet; hence,
!**** the qsw value is a little bit larger. This will avoid
!**** evaporation can generate cloud.
!
!c vf1s,vf2s=ventilation factors for snow
!c vf1s=0.78,vf2s=0.31 in LIN
!
tmp1=constc*sqrho(k)*olambdas(k)**dp5/visc(k)
tmp2=abr*xnos*( vf1s*olambdas(k)*olambdas(k)+ &
vf2s*schmidt(k)**0.33334*gamdp5o2*sqrt(tmp1) )
tmp3=amin1(0.0,tmp2)
tmp3=amax1( tmp3,tmpd )
psmltevp(k)=amax1( tmp3,-qszodt(k) )
1400 continue
!
end if
!***********************************************************************
!********* rain production processes **********
!***********************************************************************
!
!c
!c (1) AUTOCONVERSION OF RAIN (Praut): RH
!c araut=afa*rho
!c afa=0.001 Rate coefficient for autoconvergence
!c
!c araut=1.0e-3
!c
araut=0.001
!testing
! tmp1=amax1( 0.0,araut*(qlz(k)-ql0) )
! praut(k)=amin1( tmp1,qlzodt(k) )
tmp1=odtb*(qlz(k)-ql0)*( 1.0-exp(-araut*dtb) )
praut(k)=amax1( 0.0,tmp1 )
!c
!c (2) ACCRETION OF CLOUD WATER BY RAIN (Pracw): Lin (51)
!c
erw=1.0
! tmp1=qlz(k)*pio4*erw*xnor*consta*sqrho(k)
! tmp2=tmp1*gambp3*olambdar(k)**bp3
! pracw(k)=amin1( tmp2,qlzodt(k) )
tmp1=pio4*erw*xnor*consta*sqrho(k)* &
gambp3*olambdar(k)**bp3
pracw(k)=qlzodt(k)*( 1.0-exp(-tmp1*dtb) )
!c
!c (3) EVAPORATION OF RAIN (Prevp): Lin (52)
!c Prevp is negative value
!c
!c Sw=qvoqsw : saturation ratio
!c
tmpa=rvapor*xka(k)*tem(k)*tem(k)
tmpb=xlv*xlv*rho(k)*qswz(k)*diffwv(k)
tmpc=tmpa*qswz(k)*diffwv(k)
tmpd=amin1(0.0,(qvoqswz(k)-0.90)*qswz(k)*odtb)
!
! abr=2.0*pi*(qvoqswz(k)-1.0)*tmpc/(tmpa+tmpb)
abr=2.0*pi*(qvoqswz(k)-0.90)*tmpc/(tmpa+tmpb)
!
!c vf1r,vf2r=ventilation factors for rain
!c vf1r=0.78,vf2r=0.31 in RH, LIN and MM5
!
vf1r=0.78
vf2r=0.31
tmp1=consta*sqrho(k)*olambdar(k)**bp5/visc(k)
tmp2=abr*xnor*( vf1r*olambdar(k)*olambdar(k)+ &
vf2r*schmidt(k)**0.33334*gambp5o2*sqrt(tmp1) )
tmp3=amin1( 0.0,tmp2 )
tmp3=amax1( tmp3,tmpd )
prevp(k)=amax1( tmp3,-qrzodt(k) )
!
! if(iout .gt. 0) write(20,*)'tmp1,tmp2,tmp3=',tmp1,tmp2,tmp3
! if(iout .gt. 0) write(20,*)'qlz,qiz,qrz=',qlz(k),qiz(k),qrz(k)
! if(iout .gt. 0) write(20,*)'tem,qsz,qvz=',tem(k),qsz(k),qvz(k)
! if (gindex .eq. 0.) goto 900
!
if (tem(k) .lt. 273.15) then
!
!
!-- graupel
!***********************************************************************
!********* graupel production processes for T < 0 C **********
!***********************************************************************
!c
!c (1) AUTOCONVERSION OF SNOW TO FORM GRAUPEL (Pgaut): Lin (37)
!c pgaut=alpha2*(qsz-qs0)
!c qs0=6.0E-4
!c alpha2=1.0e-3*exp(0.09*temcc(k)) Lin (38)
!
alpha2=1.0e-3*exp(0.09*temcc(k))
!
! testing
! tmp1=alpha2*(qsz(k)-qs0)
! tmp1=amax1(0.0,tmp1)
! pgaut(k)=amin1( tmp1,qszodt(k) )
tmp1=odtb*(qsz(k)-qs0)*(1.0-exp(-alpha2*dtb))
pgaut(k)=amax1( 0.0,tmp1 )
!c
!c (2) FREEZING OF RAIN TO FORM GRAUPEL (Pgfr): Lin (45)
!c positive value
!c Constant in Bigg freezing Aplume=Ap=0.66 /k
!c Constant in raindrop freezing equ. Bplume=Bp=100./m/m/m/s
!
if (qrz(k) .gt. 1.e-8 ) then
Bp=100.
Ap=0.66
tmp1=olambdar(k)*olambdar(k)*olambdar(k)
tmp2=20.*pi*pi*Bp*xnor*rhowater*orho(k)* &
(exp(-Ap*temcc(k))-1.0)*tmp1*tmp1*olambdar(k)
Pgfr(k)=amin1( tmp2,qrzodt(k) )
else
Pgfr(k)=0
endif
!c
!c if (qgz(k) = 0.0) skip the other step below about graupel
!c
if (qgz(k) .eq. 0.0) goto 4000
!c
!c Comparing Pgwet(wet process) and Pdry(dry process),
!c we will pick up the small one.
!c
!c ---------------
!c | dry processes |
!c ---------------
!c
!c (3) ACCRETION OF CLOUD WATER BY GRAUPEL (Pgacw): Lin (40)
!c egw=1.0
!c Cdrag=0.6 drag coefficients for hairstone
!c constg=sqrt(4.*grav*rhograul*0.33334*orho(k)/Cdrag)
!c
egw=1.0
constg=sqrt(4.*grav*rhograul*0.33334*orho(k)*oCdrag)
tmp1=pio4*xnog*gam3pt5*constg*olambdag(k)**3.5
tmp2=qlz(k)*egw*tmp1
Pgacw(k)=amin1( tmp2,qlzodt(k) )
!c
!c (4) ACCRETION OF ICE CRYSTAL BY GRAUPEL (Pgaci): Lin (41)
!c egi=1. for wet growth
!c egi=0.1 for dry growth
!c
egi=0.1
tmp2=qiz(k)*egi*tmp1
pgaci(k)=amin1( tmp2,qizodt(k) )
!c
!c (5) ACCRETION OF SNOW BY GRAUPEL (Pgacs) : Lin (29)
!c Compute processes (6) only when qsz > 0.0 and qgz > 0.0
!c
egs=exp(0.09*temcc(k))
tmpa=olambdas(k)*olambdas(k)
tmpb=olambdag(k)*olambdag(k)
tmpc=olambdas(k)*olambdag(k)
tmp1=pi*pi*xnos*xnog*abs( vts(k)-vtg(k) )*orho(k)
tmp2=tmpa*tmpa*olambdag(k)*(5.0*tmpa+2.0*tmpc+0.5*tmpb)
tmp3=tmp1*egs*rhosnow*tmp2
Pgacs(k)=amin1( tmp3,qszodt(k) )
!c
!c (6) ACCRETION OF RAIN BY GRAUPEL (Pgacr): Lin (42)
!c Compute processes (6) only when qrz > 0.0 and qgz > 0.0
!c egr=1.
!c
egr=1.
tmpa=olambdar(k)*olambdar(k)
tmpb=olambdag(k)*olambdag(k)
tmpc=olambdar(k)*olambdag(k)
tmp1=pi*pi*xnor*xnog*abs( vtr(k)-vtg(k) )*orho(k)
tmp2=tmpa*tmpa*olambdag(k)*(5.0*tmpa+2.0*tmpc+0.5*tmpb)
tmp3=tmp1*egr*rhowater*tmp2
pgacr(k)=amin1( tmp3,qrzodt(k) )
!c
!c (7) Calculate total dry process effect Pdry(k)
!c
Pdry(k)=Pgacw(k)+pgaci(k)+Pgacs(k)+pgacr(k)
!c ---------------
!c | wet processes |
!c ---------------
!c
!c (3) ACCRETION OF ICE CRYSTAL BY GRAUPEL (Pgacip): Lin (41)
!c egi=1. for wet growth
!c egi=0.1 for dry growth
!c
tmp2=10.*pgaci(k)
pgacip(k)=amin1( tmp2,qizodt(k) )
!c
!c (4) ACCRETION OF SNOW BY GRAUPEL ((Pgacsp) : Lin (29)
!c Compute processes (6) only when qsz > 0.0 and qgz > 0.0
!c egs=exp(0.09*(tem(k)-273.15)) when T < 273.15 k
!c
tmp3=Pgacs(k)*1.0/egs
Pgacsp(k)=amin1( tmp3,qszodt(k) )
!c
!c (5) WET GROWTH OF GRAUPEL (Pgwet) : Lin (43)
!c may involve Pgacs or Pgaci and
!c must include PPgacw or Pgacr, or both.
!c ( The amount of Pgacw which is not able
!c to freeze is shed to rain. )
term0=constg*olambdag(k)**5.5/visc(k)
!c
!c vf1s,vf2s=ventilation factors for graupel
!c vf1s=0.78,vf2s=0.31 in LIN
!c Cdrag=0.6 drag coefficient for hairstone
!c constg2=vf1s*olambdag(k)*olambdag(k)+
!c vf2s*schmidt(k)**0.33334*gam2pt75*sqrt(term0)
delrs=rs0(k)-qvz(k)
tmp0=1./(xlf+cw*temcc(k))
tmp1=2.*pi*xnog*(rho(k)*xlv*diffwv(k)*delrs-xka(k)* &
temcc(k))*orho(k)*tmp0
constg2=vf1s*olambdag(k)*olambdag(k)+ &
vf2s*schmidt(k)**0.33334*gam2pt75*sqrt(term0)
tmp3=tmp1*constg2+(Pgacip(k)+Pgacsp(k))* &
(1-Ci*temcc(k)*tmp0)
tmp3=amax1(0.0,tmp3)
Pgwet(k)=amax1(tmp3,qlzodt(k)+qszodt(k)+qizodt(k) )
!c
!c Comparing Pgwet(wet process) and Pdry(dry process),
!c we will apply the small one.
!c if dry processes then delta4=1.0
!c if wet processes then delta4=0.0
!
if ( Pdry(k) .lt. Pgwet(k) ) then
delta4=1.0
else
delta4=0.0
endif
!c
!c
!c (6) Pgacrp(k)=Pgwet(k)-Pgacw(k)-Pgacip(k)-Pgacsp(k)
!c if Pgacrp(k) > 0. then some of the rain is frozen to hail
!c if Pgacrp(k) < 0. then some of the cloud water collected
!c by the hail is unable to freeze and is
!c shed as rain.
!c
Pgacrp(k)=Pgwet(k)-Pgacw(k)-Pgacip(k)-Pgacsp(k)
!c
!c (8) DEPOSITION/SUBLIMATION OF GRAUPEL (Pgdep/Pgsub): Lin (46)
!c includes ventilation effect
!c constg=sqrt(4.*grav*rhograul*0.33334*orho(k)/Cdrag)
!c constg2=vf1s*olambdag(k)*olambdag(k)+
!c vf2s*schmidt(k)**0.33334*gam2pt75*constg
!c
!c abg=2*pi*(Si-1)/rho/(A"+B")
!c
tmpa=rvapor*xka(k)*tem(k)*tem(k)
tmpb=xls*xls*rho(k)*qsiz(k)*diffwv(k)
tmpc=tmpa*qsiz(k)*diffwv(k)
abg=2.0*pi*(qvoqsiz(k)-1.0)*tmpc/(tmpa+tmpb)
!c
!c vf1s,vf2s=ventilation factors for graupel
!c vf1s=0.78,vf2s=0.31 in LIN
!c Cdrag=0.6 drag coefficient for hairstone
!c
tmp2=abg*xnog*constg2
pgdep(k)=amax1(0.0,tmp2)
pgsub(k)=amin1(0.0,tmp2)
pgsub(k)=amax1( pgsub(k),-qgzodt(k) )
4000 continue
else
!
!***********************************************************************
!********* graupel production processes for T > 0 C **********
!***********************************************************************
!
!c
!c (1) ACCRETION OF CLOUD WATER BY GRAUPEL (Pgacw): Lin (40)
!c egw=1.0
!c Cdrag=0.6 drag coefficients for hairstone
!c constg=sqrt(4.*grav*rhograul*0.33334*orho(k)/Cdrag)
egw=1.0
constg=sqrt(4.*grav*rhograul*0.33334*orho(k)*oCdrag)
tmp1=pio4*xnog*gam3pt5*constg*olambdag(k)**3.5
tmp2=qlz(k)*egw*tmp1
Pgacw(k)=amin1( tmp2,qlzodt(k) )
!c
!c (2) ACCRETION OF RAIN BY GRAUPEL (Pgacr): Lin (42)
!c Compute processes (5) only when qrz > 0.0 and qgz > 0.0
!c egr=1.
!c
egr=1.
tmpa=olambdar(k)*olambdar(k)
tmpb=olambdag(k)*olambdag(k)
tmpc=olambdar(k)*olambdag(k)
tmp1=pi*pi*xnor*xnog*abs( vtr(k)-vtg(k) )*orho(k)
tmp2=tmpa*tmpa*olambdag(k)*(5.0*tmpa+2.0*tmpc+0.5*tmpb)
tmp3=tmp1*egr*rhowater*tmp2
pgacr(k)=amin1( tmp3,qrzodt(k) )
!c
!c (3) GRAUPEL MELTING TO FORM RAIN (Pgmlt): Lin (47)
!c Pgmlt is negative value
!c constg=sqrt(4.*grav*rhograul*0.33334*orho(k)/Cdrag)
!c constg2=vf1s*olambdag(k)*olambdag(k)+
!c vf2s*schmidt(k)**0.33334*gam2pt75*constg
!c Cdrag=0.6 drag coefficients for hairstone
!
delrs=rs0(k)-qvz(k)
term1=2.0*pi*orho(k)*( xlv*diffwv(k)*rho(k)*delrs- &
xka(k)*temcc(k) )
term0=sqrt(4.*grav*rhograul*0.33334*orho(k)*ocdrag) &
*olambdag(k)**5.5/visc(k)
constg2=vf1s*olambdag(k)*olambdag(k)+ &
vf2s*schmidt(k)**0.33334*gam2pt75*sqrt(term0)
tmp2=xnog*constg2
tmp3=term1*oxlf*tmp2-cwoxlf*temcc(k)*( pgacw(k)+pgacr(k) )
tmp4=amin1(0.0,tmp3)
pgmlt(k)=amax1( tmp4,-qgzodt(k) )
!c
!c (4) EVAPORATION OF MELTING GRAUPEL (Pgmltevp) : HR (A19)
!c but use Lin et al. coefficience
!c Pgmltevp is a negative value
!c abg=2.0*pi*(qvoqsiz(k)-1.0)*tmpc/(tmpa+tmpb)
!c
tmpa=rvapor*xka(k)*tem(k)*tem(k)
tmpb=xlv*xlv*rho(k)*qswz(k)*diffwv(k)
tmpc=tmpa*qswz(k)*diffwv(k)
tmpd=amin1( 0.0,(qvoqswz(k)-0.90)*qswz(k)*odtb )
!c
!c abg=2*pi*(Si-1)/rho/(A"+B")
!c
abg=2.0*pi*(qvoqswz(k)-0.90)*tmpc/(tmpa+tmpb)
!
!**** allow evaporation to occur when RH less than 90%
!**** here not using 100% because the evaporation cooling
!**** of temperature is not taking into account yet; hence,
!**** the qgw value is a little bit larger. This will avoid
!**** evaporation can generate cloud.
!
!c vf1s,vf2s=ventilation factors for snow
!c vf1s=0.78,vf2s=0.31 in LIN
!c constg=sqrt(4.*grav*rhograul*0.33334*orho(k)/Cdrag)
!c constg2=vf1s*olambdag(k)*olambdag(k)+
!c vf2s*schmidt(k)**0.33334*gam2pt75*constg
!
tmp2=abg*xnog*constg2
tmp3=amin1(0.0,tmp2)
tmp3=amax1( tmp3,tmpd )
pgmltevp(k)=amax1( tmp3,-qgzodt(k) )
!c
!c (5) ACCRETION OF SNOW BY GRAUPEL (Pgacs) : Lin (29)
!c Compute processes (3) only when qsz > 0.0 and qgz > 0.0
!c egs=1.0
!c
egs=1.
tmpa=olambdas(k)*olambdas(k)
tmpb=olambdag(k)*olambdag(k)
tmpc=olambdas(k)*olambdag(k)
tmp1=pi*pi*xnos*xnog*abs( vts(k)-vtg(k) )*orho(k)
tmp2=tmpa*tmpa*olambdag(k)*(5.0*tmpa+2.0*tmpc+0.5*tmpb)
tmp3=tmp1*egs*rhosnow*tmp2
Pgacs(k)=amin1( tmp3,qszodt(k) )
endif
!
900 continue
!cc
!c
!c**********************************************************************
!c***** combine all processes together and avoid negative *****
!c***** water substances
!***********************************************************************
!c
if ( temcc(k) .lt. 0.0) then
!,delta4,1.-delta4
!c
!c gdelta4=gindex*delta4
!c g1sdelt4=gindex*(1.-delta4)
!c
gdelta4=gindex*delta4
g1sdelt4=gindex*(1.-delta4)
!c
!c combined water vapor depletions
!c
!cc graupel
tmp=psdep(k)+pgdep(k)*gindex
if ( tmp .gt. qvzodt(k) ) then
factor=qvzodt(k)/tmp
psdep(k)=psdep(k)*factor
pgdep(k)=pgdep(k)*factor*gindex
end if
!c
!c combined cloud water depletions
!c
tmp=praut(k)+psacw(k)+psfw(k)+pracw(k)+gindex*Pgacw(k)
if ( tmp .gt. qlzodt(k) ) then
factor=qlzodt(k)/tmp
praut(k)=praut(k)*factor
psacw(k)=psacw(k)*factor
psfw(k)=psfw(k)*factor
pracw(k)=pracw(k)*factor
!cc graupel
Pgacw(k)=Pgacw(k)*factor*gindex
end if
!c
!c combined cloud ice depletions
!c
tmp=psaut(k)+psaci(k)+praci(k)+psfi(k)+Pgaci(k)*gdelta4 &
+Pgacip(k)*g1sdelt4
if (tmp .gt. qizodt(k) ) then
factor=qizodt(k)/tmp
psaut(k)=psaut(k)*factor
psaci(k)=psaci(k)*factor
praci(k)=praci(k)*factor
psfi(k)=psfi(k)*factor
!cc graupel
Pgaci(k)=Pgaci(k)*factor*gdelta4
Pgacip(k)=Pgacip(k)*factor*g1sdelt4
endif
!c
!c combined all rain processes
!c
tmp_r=piacr(k)+psacr(k)-prevp(k)-praut(k)-pracw(k) &
+Pgfr(k)*gindex+Pgacr(k)*gdelta4 &
+Pgacrp(k)*g1sdelt4
if (tmp_r .gt. qrzodt(k) ) then
factor=qrzodt(k)/tmp_r
piacr(k)=piacr(k)*factor
psacr(k)=psacr(k)*factor
prevp(k)=prevp(k)*factor
!cc graupel
Pgfr(k)=Pgfr(k)*factor*gindex
Pgacr(k)=Pgacr(k)*factor*gdelta4
Pgacrp(k)=Pgacrp(k)*factor*g1sdelt4
endif
!c
!c if qrz < 1.0E-4 and qsz < 1.0E-4 then delta2=1.
!c (all Pracs and Psacr become to snow)
!c if qrz >= 1.0E-4 or qsz >= 1.0E-4 then delta2=0.
!c (all Pracs and Psacr become to graupel)
!c
if (qrz(k) .lt. 1.0E-4 .and. qsz(k) .lt. 1.0E-4) then
delta2=1.0
else
delta2=0.0
endif
!
!cc graupel
!c
!c if qrz(k) < 1.0e-4 then delta3=1. means praci(k) --> qs
!c piacr(k) --> qs
!c if qrz(k) > 1.0e-4 then delta3=0. means praci(k) --> qg
!c piacr(k) --> qg : Lin (20)
if (qrz(k) .lt. 1.0e-4) then
delta3=1.0
else
delta3=0.0
endif
!
!c
!c if gindex = 0.(no graupel) then delta2=1.0
!c delta3=1.0
!c
if (gindex .eq. 0.) then
delta2=1.0
delta3=1.0
endif
!
!c
!c combined all snow processes
!c
tmp_s=-pssub(k)-(psaut(k)+psaci(k)+psacw(k)+psfw(k)+ &
psfi(k)+praci(k)*delta3+piacr(k)*delta3+ &
psdep(k))+Pgaut(k)*gindex+Pgacs(k)*gdelta4+ &
Pgacsp(k)*g1sdelt4+Pracs(k)*(1.-delta2)- &
Psacr(k)*delta2
if ( tmp_s .gt. qszodt(k) ) then
factor=qszodt(k)/tmp_s
pssub(k)=pssub(k)*factor
Pracs(k)=Pracs(k)*factor
!cc graupel
Pgaut(k)=Pgaut(k)*factor*gindex
Pgacs(k)=Pgacs(k)*factor*gdelta4
Pgacsp(k)=Pgacsp(k)*factor*g1sdelt4
endif
!cc graupel
!
! if (gindex .eq. 0.) goto 998
!c
!c combined all graupel processes
!c
tmp_g=-pgaut(k)-pgfr(k)-Pgacw(k)*delta4-Pgaci(k)*delta4 &
-Pgacr(k)*delta4-Pgacs(k)*delta4 &
-pgwet(k)*(1.-delta4)-pgsub(k)-pgdep(k) &
-psacr(k)*(1-delta2)-Pracs(k)*(1-delta2) &
-praci(k)*(1-delta3)-piacr(k)*(1-delta3)
if (tmp_g .gt. qgzodt(k)) then
factor=qgzodt(k)/tmp_g
pgsub(k)=pgsub(k)*factor
endif
998 continue
!c
!c calculate new water substances, thetae, tem, and qvsbar
!c
!cc graupel
pvapor(k)=-pssub(k)-psdep(k)-prevp(k)-pgsub(k)*gindex &
-pgdep(k)*gindex
qvz(k)=amax1( qvmin,qvz(k)+dtb*pvapor(k) )
pclw(k)=-praut(k)-pracw(k)-psacw(k)-psfw(k)-pgacw(k)*gindex
qlz(k)=amax1( 0.0,qlz(k)+dtb*pclw(k) )
pcli(k)=-psaut(k)-psfi(k)-psaci(k)-praci(k)-pgaci(k)*gdelta4 &
-Pgacip(k)*g1sdelt4
qiz(k)=amax1( 0.0,qiz(k)+dtb*pcli(k) )
tmp_r=piacr(k)+psacr(k)-prevp(k)-praut(k)-pracw(k) &
+Pgfr(k)*gindex+Pgacr(k)*gdelta4 &
+Pgacrp(k)*g1sdelt4
232 format(i2,1x,6(f9.3,1x))
prain(k)=-tmp_r
qrz(k)=amax1( 0.0,qrz(k)+dtb*prain(k) )
tmp_s=-pssub(k)-(psaut(k)+psaci(k)+psacw(k)+psfw(k)+ &
psfi(k)+praci(k)*delta3+piacr(k)*delta3+ &
psdep(k))+Pgaut(k)*gindex+Pgacs(k)*gdelta4+ &
Pgacsp(k)*g1sdelt4+Pracs(k)*(1.-delta2)- &
Psacr(k)*delta2
psnow(k)=-tmp_s
qsz(k)=amax1( 0.0,qsz(k)+dtb*psnow(k) )
qschg(k)=qschg(k)+psnow(k)
qschg(k)=psnow(k)
!cc graupel
tmp_g=-pgaut(k)-pgfr(k)-Pgacw(k)*delta4-Pgaci(k)*delta4 &
-Pgacr(k)*delta4-Pgacs(k)*delta4 &
-pgwet(k)*(1.-delta4)-pgsub(k)-pgdep(k) &
-psacr(k)*(1-delta2)-Pracs(k)*(1-delta2) &
-praci(k)*(1-delta3)-piacr(k)*(1-delta3)
252 format(i2,1x,6(f12.9,1x))
262 format(i2,1x,7(f12.9,1x))
pgraupel(k)=-tmp_g
pgraupel(k)=pgraupel(k)*gindex
qgz(k)=amax1( 0.0,qgz(k)+dtb*pgraupel(k))
! qgchg(k)=qgchg(k)+pgraupel(k)
qgchg(k)=pgraupel(k)
qgz(k)=qgz(k)*gindex
tmp=ocp/tothz(k)*xLf*(qschg(k)+qgchg(k))
theiz(k)=theiz(k)+dtb*tmp
thz(k)=theiz(k)-(xLvocp*qvz(k)-xLfocp*qiz(k))/tothz(k)
tem(k)=thz(k)*tothz(k)
temcc(k)=tem(k)-273.15
if( temcc(k) .lt. -40.0 ) qswz(k)=qsiz(k)
qlpqi=qlz(k)+qiz(k)
if ( qlpqi .eq. 0.0 ) then
qvsbar(k)=qsiz(k)
else
qvsbar(k)=( qiz(k)*qsiz(k)+qlz(k)*qswz(k) )/qlpqi
endif
!
else
!c
!c combined cloud water depletions
!c
tmp=praut(k)+psacw(k)+pracw(k)+pgacw(k)*gindex
if ( tmp .gt. qlzodt(k) ) then
factor=qlzodt(k)/tmp
praut(k)=praut(k)*factor
psacw(k)=psacw(k)*factor
pracw(k)=pracw(k)*factor
!cc graupel
pgacw(k)=pgacw(k)*factor*gindex
end if
!c
!c combined all snow processes
!c
tmp_s=-(psmlt(k)+psmltevp(k))+Pgacs(k)*gindex
if (tmp_s .gt. qszodt(k) ) then
factor=qszodt(k)/tmp_s
psmlt(k)=psmlt(k)*factor
psmltevp(k)=psmltevp(k)*factor
!cc graupel
Pgacs(k)=Pgacs(k)*factor*gindex
endif
!c
!c
!cc graupel
!c
! if (gindex .eq. 0.) goto 997
!c
!c combined all graupel processes
!c
tmp_g=-pgmlt(k)-pgacs(k)-pgmltevp(k)
if (tmp_g .gt. qgzodt(k)) then
factor=qgzodt(k)/tmp_g
pgmltevp(k)=pgmltevp(k)*factor
pgmlt(k)=pgmlt(k)*factor
endif
!c
997 continue
!c
!c combined all rain processes
!c
tmp_r=-prevp(k)-(praut(k)+pracw(k)+psacw(k)-psmlt(k)) &
+pgmlt(k)*gindex-pgacw(k)*gindex
if (tmp_r .gt. qrzodt(k) ) then
factor=qrzodt(k)/tmp_r
prevp(k)=prevp(k)*factor
endif
!c
!c
!c calculate new water substances and thetae
!c
pvapor(k)=-psmltevp(k)-prevp(k)-pgmltevp(k)
qvz(k)=amax1( qvmin,qvz(k)+dtb*pvapor(k))
pclw(k)=-praut(k)-pracw(k)-psacw(k)-pgacw(k)*gindex
qlz(k)=amax1( 0.0,qlz(k)+dtb*pclw(k) )
pcli(k)=0.0
qiz(k)=amax1( 0.0,qiz(k)+dtb*pcli(k) )
tmp_r=-prevp(k)-(praut(k)+pracw(k)+psacw(k)-psmlt(k)) &
+pgmlt(k)*gindex-pgacw(k)*gindex
242 format(i2,1x,7(f9.6,1x))
prain(k)=-tmp_r
tmpqrz=qrz(k)
qrz(k)=amax1( 0.0,qrz(k)+dtb*prain(k) )
tmp_s=-(psmlt(k)+psmltevp(k))+Pgacs(k)*gindex
psnow(k)=-tmp_s
qsz(k)=amax1( 0.0,qsz(k)+dtb*psnow(k) )
! qschg(k)=qschg(k)+psnow(k)
qschg(k)=psnow(k)
!cc graupel
tmp_g=-pgmlt(k)-pgacs(k)-pgmltevp(k)
! write(*,272)k,pgmlt(k),pgacs(k),pgmltevp(k),
272 format(i2,1x,3(f12.9,1x))
pgraupel(k)=-tmp_g*gindex
qgz(k)=amax1( 0.0,qgz(k)+dtb*pgraupel(k))
! qgchg(k)=qgchg(k)+pgraupel(k)
qgchg(k)=pgraupel(k)
qgz(k)=qgz(k)*gindex
!
tmp=ocp/tothz(k)*xLf*(qschg(k)+qgchg(k))
theiz(k)=theiz(k)+dtb*tmp
thz(k)=theiz(k)-(xLvocp*qvz(k)-xLfocp*qiz(k))/tothz(k)
tem(k)=thz(k)*tothz(k)
temcc(k)=tem(k)-273.15
! qswz(k)=episp0k*oprez(k)* &
! exp( svp2*temcc(k)/(tem(k)-svp3) )
es=1000.*svp1*exp( svp2*temcc(k)/(tem(k)-svp3) )
qswz(k)=ep2*es/(prez(k)-es)
qsiz(k)=qswz(k)
qvsbar(k)=qswz(k)
!
end if
!
!***********************************************************************
!********** saturation adjustment **********
!***********************************************************************
!
! allow supersaturation exits linearly from 0% at 500 mb to 50%
! above 300 mb
! 5.0e-5=1.0/(500mb-300mb)
!
rsat=1.0+0.5*(50000.0-prez(k))*5.0e-5
rsat=amax1(1.0,rsat)
rsat=amin1(1.5,rsat)
rsat=1.0
if( qvz(k)+qlz(k)+qiz(k) .lt. rsat*qvsbar(k) ) then
!c
!c unsaturated
!c
qvz(k)=qvz(k)+qlz(k)+qiz(k)
qlz(k)=0.0
qiz(k)=0.0
thz(k)=theiz(k)-(xLvocp*qvz(k)-xLfocp*qiz(k))/tothz(k)
tem(k)=thz(k)*tothz(k)
temcc(k)=tem(k)-273.15
go to 1800
!
else
!c
!c saturated
!c
!
pladj(k)=qlz(k)
piadj(k)=qiz(k)
!
CALL satadj(qvz, qlz, qiz, prez, theiz, thz, tothz, kts, kte, &
k, xLvocp, xLfocp, episp0k, EP2,SVP1,SVP2,SVP3,SVPT0 )
!
pladj(k)=odtb*(qlz(k)-pladj(k))
piadj(k)=odtb*(qiz(k)-piadj(k))
!
pclw(k)=pclw(k)+pladj(k)
pcli(k)=pcli(k)+piadj(k)
pvapor(k)=pvapor(k)-( pladj(k)+piadj(k) )
!
thz(k)=theiz(k)-(xLvocp*qvz(k)-xLfocp*qiz(k))/tothz(k)
tem(k)=thz(k)*tothz(k)
temcc(k)=tem(k)-273.15
! qswz(k)=episp0k*oprez(k)* &
! exp( svp2*temcc(k)/(tem(k)-svp3) )
es=1000.*svp1*exp( svp2*temcc(k)/(tem(k)-svp3) )
qswz(k)=ep2*es/(prez(k)-es)
if (tem(k) .lt. 273.15 ) then
! qsiz(k)=episp0k*oprez(k)* &
! exp( 21.8745584*(tem(k)-273.16)/(tem(k)-7.66) )
es=1000.*svp1*exp( 21.8745584*(tem(k)-273.16)/(tem(k)-7.66) )
qsiz(k)=ep2*es/(prez(k)-es)
if (temcc(k) .lt. -40.0) qswz(k)=qsiz(k)
else
qsiz(k)=qswz(k)
endif
qlpqi=qlz(k)+qiz(k)
if ( qlpqi .eq. 0.0 ) then
qvsbar(k)=qsiz(k)
else
qvsbar(k)=( qiz(k)*qsiz(k)+qlz(k)*qswz(k) )/qlpqi
endif
end if
!
!***********************************************************************
!***** melting and freezing of cloud ice and cloud water *****
!***********************************************************************
qlpqi=qlz(k)+qiz(k)
if(qlpqi .le. 0.0) go to 1800
!
!c
!c (1) HOMOGENEOUS NUCLEATION WHEN T< -40 C (Pihom)
!c
if(temcc(k) .lt. -40.0) pihom(k)=qlz(k)*odtb
!c
!c (2) MELTING OF ICE CRYSTAL WHEN T> 0 C (Pimlt)
!c
if(temcc(k) .gt. 0.0) pimlt(k)=qiz(k)*odtb
!c
!c (3) PRODUCTION OF CLOUD ICE BY BERGERON PROCESS (Pidw): Hsie (p957)
!c this process only considered when -31 C < T < 0 C
!c
if(temcc(k) .lt. 0.0 .and. temcc(k) .gt. -31.0) then
!c!
!c! parama1 and parama2 functions must be user supplied
!c!
a1=parama1( temcc(k) )
a2=parama2( temcc(k) )
!! change unit from cgs to mks
a1=a1*0.001**(1.0-a2)
xnin=xni0*exp(-bni*temcc(k))
pidw(k)=xnin*orho(k)*(a1*xmnin**a2)
end if
!
pcli(k)=pcli(k)+pihom(k)-pimlt(k)+pidw(k)
pclw(k)=pclw(k)-pihom(k)+pimlt(k)-pidw(k)
qlz(k)=amax1( 0.0,qlz(k)+dtb*(-pihom(k)+pimlt(k)-pidw(k)) )
qiz(k)=amax1( 0.0,qiz(k)+dtb*(pihom(k)-pimlt(k)+pidw(k)) )
!
CALL satadj(qvz, qlz, qiz, prez, theiz, thz, tothz, kts, kte, &
k, xLvocp, xLfocp, episp0k ,EP2,SVP1,SVP2,SVP3,SVPT0)
thz(k)=theiz(k)-(xLvocp*qvz(k)-xLfocp*qiz(k))/tothz(k)
tem(k)=thz(k)*tothz(k)
temcc(k)=tem(k)-273.15
! qswz(k)=episp0k*oprez(k)* &
! exp( svp2*temcc(k)/(tem(k)-svp3) )
es=1000.*svp1*exp( svp2*temcc(k)/(tem(k)-svp3) )
qswz(k)=ep2*es/(prez(k)-es)
if (tem(k) .lt. 273.15 ) then
! qsiz(k)=episp0k*oprez(k)* &
! exp( 21.8745584*(tem(k)-273.16)/(tem(k)-7.66) )
es=1000.*svp1*exp( 21.8745584*(tem(k)-273.16)/(tem(k)-7.66) )
qsiz(k)=ep2*es/(prez(k)-es)
if (temcc(k) .lt. -40.0) qswz(k)=qsiz(k)
else
qsiz(k)=qswz(k)
endif
qlpqi=qlz(k)+qiz(k)
if ( qlpqi .eq. 0.0 ) then
qvsbar(k)=qsiz(k)
else
qvsbar(k)=( qiz(k)*qsiz(k)+qlz(k)*qswz(k) )/qlpqi
endif
1800 continue
!
!***********************************************************************
!********** integrate the productions of rain and snow **********
!***********************************************************************
!c
2000 continue
!---------------------------------------------------------------------
!
!***********************************************************************
!****** Write terms in cloud physics to time series dataset *****
!***********************************************************************
!
! open(unit=24,form='formatted',status='new',
! & file='cloud.dat')
!9030 format(10e12.6)
! write(24,*)'tmp'
! write(24,9030) (tem(k),k=kts+1,kte)
! write(24,*)'qiz'
! write(24,9030) (qiz(k),k=kts+1,kte)
! write(24,*)'qsz'
! write(24,9030) (qsz(k),k=kts+1,kte)
! write(24,*)'qrz'
! write(24,9030) (qrz(k),k=kts+1,kte)
! write(24,*)'qgz'
! write(24,9030) (qgz(k),k=kts+1,kte)
! write(24,*)'qvoqsw'
! write(24,9030) (qvoqswz(k),k=kts+1,kte)
! write(24,*)'qvoqsi'
! write(24,9030) (qvoqsiz(k),k=kts+1,kte)
! write(24,*)'vtr'
! write(24,9030) (vtr(k),k=kts+1,kte)
! write(24,*)'vts'
! write(24,9030) (vts(k),k=kts+1,kte)
! write(24,*)'vtg'
! write(24,9030) (vtg(k),k=kts+1,kte)
! write(24,*)'pclw'
! write(24,9030) (pclw(k),k=kts+1,kte)
! write(24,*)'pvapor'
! write(24,9030) (pvapor(k),k=kts+1,kte)
! write(24,*)'pcli'
! write(24,9030) (pcli(k),k=kts+1,kte)
! write(24,*)'pimlt'
! write(24,9030) (pimlt(k),k=kts+1,kte)
! write(24,*)'pihom'
! write(24,9030) (pihom(k),k=kts+1,kte)
! write(24,*)'pidw'
! write(24,9030) (pidw(k),k=kts+1,kte)
! write(24,*)'prain'
! write(24,9030) (prain(k),k=kts+1,kte)
! write(24,*)'praut'
! write(24,9030) (praut(k),k=kts+1,kte)
! write(24,*)'pracw'
! write(24,9030) (pracw(k),k=kts+1,kte)
! write(24,*)'prevp'
! write(24,9030) (prevp(k),k=kts+1,kte)
! write(24,*)'psnow'
! write(24,9030) (psnow(k),k=kts+1,kte)
! write(24,*)'psaut'
! write(24,9030) (psaut(k),k=kts+1,kte)
! write(24,*)'psfw'
! write(24,9030) (psfw(k),k=kts+1,kte)
! write(24,*)'psfi'
! write(24,9030) (psfi(k),k=kts+1,kte)
! write(24,*)'praci'
! write(24,9030) (praci(k),k=kts+1,kte)
! write(24,*)'piacr'
! write(24,9030) (piacr(k),k=kts+1,kte)
! write(24,*)'psaci'
! write(24,9030) (psaci(k),k=kts+1,kte)
! write(24,*)'psacw'
! write(24,9030) (psacw(k),k=kts+1,kte)
! write(24,*)'psdep'
! write(24,9030) (psdep(k),k=kts+1,kte)
! write(24,*)'pssub'
! write(24,9030) (pssub(k),k=kts+1,kte)
! write(24,*)'pracs'
! write(24,9030) (pracs(k),k=kts+1,kte)
! write(24,*)'psacr'
! write(24,9030) (psacr(k),k=kts+1,kte)
! write(24,*)'psmlt'
! write(24,9030) (psmlt(k),k=kts+1,kte)
! write(24,*)'psmltevp'
! write(24,9030) (psmltevp(k),k=kts+1,kte)
! write(24,*)'pladj'
! write(24,9030) (pladj(k),k=kts+1,kte)
! write(24,*)'piadj'
! write(24,9030) (piadj(k),k=kts+1,kte)
! write(24,*)'pgraupel'
! write(24,9030) (pgraupel(k),k=kts+1,kte)
! write(24,*)'pgaut'
! write(24,9030) (pgaut(k),k=kts+1,kte)
! write(24,*)'pgfr'
! write(24,9030) (pgfr(k),k=kts+1,kte)
! write(24,*)'pgacw'
! write(24,9030) (pgacw(k),k=kts+1,kte)
! write(24,*)'pgaci'
! write(24,9030) (pgaci(k),k=kts+1,kte)
! write(24,*)'pgacr'
! write(24,9030) (pgacr(k),k=kts+1,kte)
! write(24,*)'pgacs'
! write(24,9030) (pgacs(k),k=kts+1,kte)
! write(24,*)'pgacip'
! write(24,9030) (pgacip(k),k=kts+1,kte)
! write(24,*)'pgacrP'
! write(24,9030) (pgacrP(k),k=kts+1,kte)
! write(24,*)'pgacsp'
! write(24,9030) (pgacsp(k),k=kts+1,kte)
! write(24,*)'pgwet'
! write(24,9030) (pgwet(k),k=kts+1,kte)
! write(24,*)'pdry'
! write(24,9030) (pdry(k),k=kts+1,kte)
! write(24,*)'pgsub'
! write(24,9030) (pgsub(k),k=kts+1,kte)
! write(24,*)'pgdep'
! write(24,9030) (pgdep(k),k=kts+1,kte)
! write(24,*)'pgmlt'
! write(24,9030) (pgmlt(k),k=kts+1,kte)
! write(24,*)'pgmltevp'
! write(24,9030) (pgmltevp(k),k=kts+1,kte)
!**** below if qv < qvmin then qv=qvmin, ql=0.0, and qi=0.0
!
do k=kts+1,kte
if ( qvz(k) .lt. qvmin ) then
qlz(k)=0.0
qiz(k)=0.0
qvz(k)=amax1( qvmin,qvz(k)+qlz(k)+qiz(k) )
end if
enddo
!
END SUBROUTINE clphy1d
!---------------------------------------------------------------------
! SATURATED ADJUSTMENT
!---------------------------------------------------------------------
SUBROUTINE satadj(qvz, qlz, qiz, prez, theiz, thz, tothz, & 2
kts, kte, k, xLvocp, xLfocp, episp0k, EP2,SVP1,SVP2,SVP3,SVPT0)
!---------------------------------------------------------------------
IMPLICIT NONE
!---------------------------------------------------------------------
! This program use Newton's method for finding saturated temperature
! and saturation mixing ratio.
!
! In this saturation adjustment scheme we assume
! (1) the saturation mixing ratio is the mass weighted average of
! saturation values over liquid water (qsw), and ice (qsi)
! following Lord et al., 1984 and Tao, 1989
!
! (2) the percentage of cloud liquid and cloud ice will
! be fixed during the saturation calculation
!---------------------------------------------------------------------
!
INTEGER, INTENT(IN ) :: kts, kte, k
REAL, DIMENSION( kts:kte ), &
INTENT(INOUT) :: qvz, qlz, qiz
!
REAL, DIMENSION( kts:kte ), &
INTENT(IN ) :: prez, theiz, tothz
REAL, INTENT(IN ) :: xLvocp, xLfocp, episp0k
REAL, INTENT(IN ) :: EP2,SVP1,SVP2,SVP3,SVPT0
! LOCAL VARS
INTEGER :: n
REAL, DIMENSION( kts:kte ) :: thz, tem, temcc, qsiz, &
qswz, qvsbar
REAL :: qsat, qlpqi, ratql, t0, t1, tmp1, ratqi, tsat, absft, &
denom1, denom2, dqvsbar, ftsat, dftsat, qpz, qvold, &
gindex, es
!
!---------------------------------------------------------------------
thz(k)=theiz(k)-(xLvocp*qvz(k)-xLfocp*qiz(k))/tothz(k)
tem(k)=tothz(k)*thz(k)
if (tem(k) .gt. 273.15) then
! qsat=episp0k/prez(k)* &
! exp( svp2*(tem(k)-273.15)/(tem(k)-svp3) )
es=1000.*svp1*exp( svp2*(tem(k)-svpt0)/(tem(k)-svp3) )
qsat=ep2*es/(prez(k)-es)
else
qsat=episp0k/prez(k)* &
exp( 21.8745584*(tem(k)-273.15)/(tem(k)-7.66) )
end if
qpz=qvz(k)+qlz(k)+qiz(k)
qvold=qvz(k)
if (qpz .lt. qsat) then
qvz(k)=qpz
qiz(k)=0.0
qlz(k)=0.0
go to 400
end if
qlpqi=qlz(k)+qiz(k)
if( qlpqi .ge. 1.0e-5) then
ratql=qlz(k)/qlpqi
ratqi=qiz(k)/qlpqi
else
t0=273.15
! t1=233.15
t1=248.15
tmp1=( t0-tem(k) )/(t0-t1)
tmp1=amin1(1.0,tmp1)
tmp1=amax1(0.0,tmp1)
ratqi=tmp1
ratql=1.0-tmp1
end if
!
!
!-- saturation mixing ratios over water and ice
!-- at the outset we will follow Bolton 1980 MWR for
!-- the water and Murray JAS 1967 for the ice
!
!-- dqvsbar=d(qvsbar)/dT
!-- ftsat=F(Tsat)
!-- dftsat=d(F(T))/dT
!
! First guess of tsat
tsat=tem(k)
absft=1.0
!
do 200 n=1,20
denom1=1.0/(tsat-svp3)
denom2=1.0/(tsat-7.66)
! qswz(k)=episp0k/prez(k)* &
! exp( svp2*denom1*(tsat-273.15) )
es=1000.*svp1*exp( svp2*denom1*(tsat-svpt0) )
qswz(k)=ep2*es/(prez(k)-es)
if (tem(k) .lt. 273.15) then
! qsiz(k)=episp0k/prez(k)* &
! exp( 21.8745584*denom2*(tsat-273.15) )
es=1000.*svp1*exp( 21.8745584*denom2*(tsat-273.15) )
qsiz(k)=ep2*es/(prez(k)-es)
if (tem(k) .lt. 233.15) qswz(k)=qsiz(k)
else
qsiz(k)=qswz(k)
endif
qvsbar(k)=ratql*qswz(k)+ratqi*qsiz(k)
!
! if( absft .lt. 0.01 .and. n .gt. 3 ) go to 300
if( absft .lt. 0.01 ) go to 300
!
dqvsbar=ratql*qswz(k)*svp2*243.5*denom1*denom1+ &
ratqi*qsiz(k)*21.8745584*265.5*denom2*denom2
ftsat=tsat+(xlvocp+ratqi*xlfocp)*qvsbar(k)- &
tothz(k)*theiz(k)-xlfocp*ratqi*(qvz(k)+qlz(k)+qiz(k))
dftsat=1.0+(xlvocp+ratqi*xlfocp)*dqvsbar
tsat=tsat-ftsat/dftsat
absft=abs(ftsat)
200 continue
9020 format(1x,'point can not converge, absft,n=',e12.5,i5)
!
300 continue
if( qpz .gt. qvsbar(k) ) then
qvz(k)=qvsbar(k)
qiz(k)=ratqi*( qpz-qvz(k) )
qlz(k)=ratql*( qpz-qvz(k) )
else
qvz(k)=qpz
qiz(k)=0.0
qlz(k)=0.0
end if
400 continue
END SUBROUTINE satadj
!----------------------------------------------------------------
REAL FUNCTION parama1(temp) 2
!----------------------------------------------------------------
IMPLICIT NONE
!----------------------------------------------------------------
! This program calculate the parameter for crystal growth rate
! in Bergeron process
!----------------------------------------------------------------
REAL, INTENT (IN ) :: temp
REAL, DIMENSION(32) :: a1
INTEGER :: i1, i1p1
REAL :: ratio
data a1/0.100e-10,0.7939e-7,0.7841e-6,0.3369e-5,0.4336e-5, &
0.5285e-5,0.3728e-5,0.1852e-5,0.2991e-6,0.4248e-6, &
0.7434e-6,0.1812e-5,0.4394e-5,0.9145e-5,0.1725e-4, &
0.3348e-4,0.1725e-4,0.9175e-5,0.4412e-5,0.2252e-5, &
0.9115e-6,0.4876e-6,0.3473e-6,0.4758e-6,0.6306e-6, &
0.8573e-6,0.7868e-6,0.7192e-6,0.6513e-6,0.5956e-6, &
0.5333e-6,0.4834e-6/
i1=int(-temp)+1
i1p1=i1+1
ratio=-(temp)-float(i1-1)
parama1=a1(i1)+ratio*( a1(i1p1)-a1(i1) )
END FUNCTION parama1
!----------------------------------------------------------------
REAL FUNCTION parama2(temp) 2
!----------------------------------------------------------------
IMPLICIT NONE
!----------------------------------------------------------------
! This program calculate the parameter for crystal growth rate
! in Bergeron process
!----------------------------------------------------------------
REAL, INTENT (IN ) :: temp
REAL, DIMENSION(32) :: a2
INTEGER :: i1, i1p1
REAL :: ratio
data a2/0.0100,0.4006,0.4831,0.5320,0.5307,0.5319,0.5249, &
0.4888,0.3849,0.4047,0.4318,0.4771,0.5183,0.5463, &
0.5651,0.5813,0.5655,0.5478,0.5203,0.4906,0.4447, &
0.4126,0.3960,0.4149,0.4320,0.4506,0.4483,0.4460, &
0.4433,0.4413,0.4382,0.4361/
i1=int(-temp)+1
i1p1=i1+1
ratio=-(temp)-float(i1-1)
parama2=a2(i1)+ratio*( a2(i1p1)-a2(i1) )
END FUNCTION parama2
!----------------------------------------------------------------
REAL FUNCTION ggamma(X) 10,1
!----------------------------------------------------------------
IMPLICIT NONE
!----------------------------------------------------------------
REAL, INTENT(IN ) :: x
REAL, DIMENSION(8) :: B
INTEGER ::j, K1
REAL ::PF, G1TO2 ,TEMP
DATA B/-.577191652,.988205891,-.897056937,.918206857, &
-.756704078,.482199394,-.193527818,.035868343/
PF=1.
TEMP=X
DO 10 J=1,200
IF (TEMP .LE. 2) GO TO 20
TEMP=TEMP-1.
10 PF=PF*TEMP
100 FORMAT(//,5X,'module_mp_lin: INPUT TO GAMMA FUNCTION TOO LARGE, X=',E12.5)
WRITE(wrf_err_message,100)X
CALL wrf_error_fatal(wrf_err_message)
20 G1TO2=1.
TEMP=TEMP - 1.
DO 30 K1=1,8
30 G1TO2=G1TO2 + B(K1)*TEMP**K1
ggamma=PF*G1TO2
END FUNCTION ggamma
!----------------------------------------------------------------
END MODULE module_mp_lin