!WRF:MODEL_LAYER:PHYSICS
!
MODULE module_mp_gsfcgce 1
USE module_wrf_error
USE module_utility, ONLY: WRFU_Clock, WRFU_Alarm
USE module_domain, ONLY : HISTORY_ALARM, Is_alarm_tstep
USE module_mp_radar
!JJS 1/3/2008 vvvvv
! common /bt/
REAL, PRIVATE :: rd1, rd2, al, cp
! common /cont/
REAL, PRIVATE :: c38, c358, c610, c149, &
c879, c172, c409, c76, &
c218, c580, c141
! common /b3cs/
REAL, PRIVATE :: ag, bg, as, bs, &
aw, bw, bgh, bgq, &
bsh, bsq, bwh, bwq
! common /size/
REAL, PRIVATE :: tnw, tns, tng, &
roqs, roqg, roqr
! common /bterv/
REAL, PRIVATE :: zrc, zgc, zsc, &
vrc, vgc, vsc
! common /bsnw/
REAL, PRIVATE :: alv, alf, als, t0, t00, &
avc, afc, asc, rn1, bnd1, &
rn2, bnd2, rn3, rn4, rn5, &
rn6, rn7, rn8, rn9, rn10, &
rn101,rn10a, rn11,rn11a, rn12
REAL, PRIVATE :: rn14, rn15,rn15a, rn16, rn17, &
rn17a,rn17b,rn17c, rn18, rn18a, &
rn19,rn19a,rn19b, rn20, rn20a, &
rn20b, bnd3, rn21, rn22, rn23, &
rn23a,rn23b, rn25,rn30a, rn30b, &
rn30c, rn31, beta, rn32
REAL, PRIVATE, DIMENSION( 31 ) :: rn12a, rn12b, rn13, rn25a
! common /rsnw1/
REAL, PRIVATE :: rn10b, rn10c, rnn191, rnn192, rn30, &
rnn30a, rn33, rn331, rn332
!
REAL, PRIVATE, DIMENSION( 31 ) :: aa1, aa2
DATA aa1/.7939e-7, .7841e-6, .3369e-5, .4336e-5, .5285e-5, &
.3728e-5, .1852e-5, .2991e-6, .4248e-6, .7434e-6, &
.1812e-5, .4394e-5, .9145e-5, .1725e-4, .3348e-4, &
.1725e-4, .9175e-5, .4412e-5, .2252e-5, .9115e-6, &
.4876e-6, .3473e-6, .4758e-6, .6306e-6, .8573e-6, &
.7868e-6, .7192e-6, .6513e-6, .5956e-6, .5333e-6, &
.4834e-6/
DATA aa2/.4006, .4831, .5320, .5307, .5319, &
.5249, .4888, .3894, .4047, .4318, &
.4771, .5183, .5463, .5651, .5813, &
.5655, .5478, .5203, .4906, .4447, &
.4126, .3960, .4149, .4320, .4506, &
.4483, .4460, .4433, .4413, .4382, &
.4361/
!+---+-----------------------------------------------------------------+
!..The following 6 variables moved here to facilitate reflectivity
!.. calculation similar to other MP schemes, because when they get
!.. declared later in the code (now commented out), it makes things
!.. more difficult to integreate with the radar code.
REAL , PARAMETER :: xnor = 8.0e6
REAL , PARAMETER :: xnos = 1.6e7
REAL , PARAMETER :: xnoh = 2.0e5
REAL , PARAMETER :: xnog = 4.0e6
REAL , PARAMETER :: rhohail = 917.
REAL , PARAMETER :: rhograul = 400.
!+---+-----------------------------------------------------------------+
!JJS 1/3/2008 ^^^^^
CONTAINS
!-------------------------------------------------------------------
! NASA/GSFC GCE
! Tao et al, 2001, Meteo. & Atmos. Phy., 97-137
!-------------------------------------------------------------------
! SUBROUTINE gsfcgce( th, th_old, &
SUBROUTINE gsfcgce( th, & 1,10
qv, ql, &
qr, qi, &
qs, &
! qvold, qlold, &
! qrold, qiold, &
! qsold, &
rho, pii, p, dt_in, z, &
ht, dz8w, grav, &
rhowater, rhosnow, &
itimestep, &
ids,ide, jds,jde, kds,kde, & ! domain dims
ims,ime, jms,jme, kms,kme, & ! memory dims
its,ite, jts,jte, kts,kte, & ! tile dims
rainnc, rainncv, &
snownc, snowncv, sr, &
graupelnc, graupelncv, &
refl_10cm, diagflag, do_radar_ref, &
! f_qg, qg, pgold, &
f_qg, qg, &
ihail, ice2 &
)
!-------------------------------------------------------------------
IMPLICIT NONE
!-------------------------------------------------------------------
!
! JJS 2/15/2005
!
INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde , &
ims,ime, jms,jme, kms,kme , &
its,ite, jts,jte, kts,kte
INTEGER, INTENT(IN ) :: itimestep, ihail, ice2
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 ) :: &
! th_old, &
! qvold, &
! qlold, &
! qrold, &
! qiold, &
! qsold, &
! qgold, &
rho, &
pii, &
p, &
dz8w, &
z
REAL, DIMENSION( ims:ime , jms:jme ), &
INTENT(INOUT) :: rainnc, &
rainncv, &
snownc, &
snowncv, &
sr, &
graupelnc, &
graupelncv
!+---+-----------------------------------------------------------------+
REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT):: & ! GT
refl_10cm
LOGICAL, OPTIONAL, INTENT(IN) :: diagflag
INTEGER, OPTIONAL, INTENT(IN) :: do_radar_ref
!+---+-----------------------------------------------------------------+
REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN) :: ht
REAL, INTENT(IN ) :: dt_in, &
grav, &
rhowater, &
rhosnow
LOGICAL, INTENT(IN), OPTIONAL :: F_QG
! LOCAL VAR
!jjs INTEGER :: min_q, max_q
!jjs REAL, DIMENSION( its:ite , jts:jte ) &
!jjs :: rain, snow, graupel,ice
!
! INTEGER :: IHAIL, itaobraun, ice2, istatmin, new_ice_sat, id
INTEGER :: itaobraun, istatmin, new_ice_sat, id
INTEGER :: i, j, k
INTEGER :: iskip, ih, icount, ibud, i24h
REAL :: hour
REAL , PARAMETER :: cmin=1.e-20
REAL :: dth, dqv, dqrest, dqall, dqall1, rhotot, a1, a2
! REAL, DIMENSION( its:ite , kts:kte , jts:jte ) :: &
! th1, qv1, ql1, qr1, qi1, qs1, qg1
LOGICAL :: flag_qg
!+---+-----------------------------------------------------------------+
INTEGER:: NCALL=0
IF (NCALL .EQ. 0) THEN
!..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 = 3.14159*rhowater/6.
xbm_r = 3.
xmu_r = 0.
xam_s = 3.14159*rhosnow/6.
xbm_s = 3.
xmu_s = 0.
if (ihail .eq. 1) then
xam_g = 3.14159*rhohail/6.
else
xam_g = 3.14159*rhograul/6.
endif
xbm_g = 3.
xmu_g = 0.
call radar_init
NCALL = 1
ENDIF
!+---+-----------------------------------------------------------------+
!
!c ihail = 0 for graupel, for tropical region
!c ihail = 1 for hail, for mid-lat region
! itaobraun: 0 for Tao's constantis, 1 for Braun's constants
!c if ( itaobraun.eq.1 ) --> betah=0.5*beta=-.46*0.5=-0.23; cn0=1.e-6
!c if ( itaobraun.eq.0 ) --> betah=0.5*beta=-.6*0.5=-0.30; cn0=1.e-8
itaobraun = 1
!c ice2 = 0 for 3 ice --- ice, snow and graupel/hail
!c ice2 = 1 for 2 ice --- ice and snow only
!c ice2 = 2 for 2 ice --- ice and graupel only, use ihail = 0 only
!c ice2 = 3 for 0 ice --- no ice, warm only
! if (ice2 .eq. 2) ihail = 0
i24h=nint(86400./dt_in)
if (mod(itimestep,i24h).eq.1) then
write(6,*) 'ihail=',ihail,' ice2=',ice2
if (ice2.eq.0) then
write(6,*) 'Running 3-ice scheme in GSFCGCE with'
if (ihail.eq.0) then
write(6,*) ' ice, snow and graupel'
else if (ihail.eq.1) then
write(6,*) ' ice, snow and hail'
else
write(6,*) 'ihail has to be either 1 or 0'
stop
endif !ihail
else if (ice2.eq.1) then
write(6,*) 'Running 2-ice scheme in GSFCGCE with'
write(6,*) ' ice and snow'
else if (ice2.eq.2) then
write(6,*) 'Running 2-ice scheme in GSFCGCE with'
write(6,*) ' ice and graupel'
else if (ice2.eq.3) then
write(6,*) 'Running warm rain only scheme in GSFCGCE without any ice'
else
write(6,*) 'gsfcgce_2ice in namelist.input has to be 0, 1, 2, or 3'
stop
endif !ice2
endif !itimestep
!c new_ice_sat = 0, 1 or 2
new_ice_sat = 2
!c istatmin
istatmin = 180
!c id = 0 without in-line staticstics
!c id = 1 with in-line staticstics
id = 0
!c ibud = 0 no calculation of dth, dqv, dqrest and dqall
!c ibud = 1 yes
ibud = 0
!jjs dt=dt_in
!jjs rhoe_s=1.29
!
! IF (P_QI .lt. P_FIRST_SCALAR .or. P_QS .lt. P_FIRST_SCALAR) THEN
! CALL wrf_error_fatal3 ( "module_mp_lin.b" , 130 , 'module_mp_lin: Improper use of Lin et al scheme; no ice phase. Please chose another one.')
! ENDIF
! calculte fallflux and precipiation in MKS system
call fall_flux(dt_in, qr, qi, qs, qg, p, &
rho, z, dz8w, ht, rainnc, &
rainncv, grav,itimestep, &
rhowater, rhosnow, &
snownc, snowncv, sr, &
graupelnc, graupelncv, &
ihail, ice2, &
ims,ime, jms,jme, kms,kme, & ! memory dims
its,ite, jts,jte, kts,kte ) ! tile dims
!-----------------------------------------------------------------------
!c set up constants used internally in GCE
call consat_s (ihail, itaobraun)
!c Negative values correction
iskip = 1
if (iskip.eq.0) then
call negcor(qv,rho,dz8w,ims,ime,jms,jme,kms,kme, &
itimestep,1, &
its,ite,jts,jte,kts,kte)
call negcor(ql,rho,dz8w,ims,ime,jms,jme,kms,kme, &
itimestep,2, &
its,ite,jts,jte,kts,kte)
call negcor(qr,rho,dz8w,ims,ime,jms,jme,kms,kme, &
itimestep,3, &
its,ite,jts,jte,kts,kte)
call negcor(qi,rho,dz8w,ims,ime,jms,jme,kms,kme, &
itimestep,4, &
its,ite,jts,jte,kts,kte)
call negcor(qs,rho,dz8w,ims,ime,jms,jme,kms,kme, &
itimestep,5, &
its,ite,jts,jte,kts,kte)
call negcor(qg,rho,dz8w,ims,ime,jms,jme,kms,kme, &
itimestep,6, &
its,ite,jts,jte,kts,kte)
! else if (mod(itimestep,i24h).eq.1) then
! print *,'no neg correction in mp at timestep=',itimestep
endif ! iskip
!c microphysics in GCE
call SATICEL_S( dt_in, IHAIL, itaobraun, ICE2, istatmin, &
new_ice_sat, id, &
! th, th_old, qv, ql, qr, &
th, qv, ql, qr, &
qi, qs, qg, &
! qvold, qlold, qrold, &
! qiold, qsold, qgold, &
rho, pii, p, itimestep, &
refl_10cm, diagflag, do_radar_ref, & ! GT added for reflectivity calcs
! refl_10cm, grid_clock, grid_alarms, & ! GT added for reflectivity calcs
ids,ide, jds,jde, kds,kde, & ! domain dims
ims,ime, jms,jme, kms,kme, & ! memory dims
its,ite, jts,jte, kts,kte & ! tile dims
)
END SUBROUTINE gsfcgce
!-----------------------------------------------------------------------
SUBROUTINE fall_flux ( dt, qr, qi, qs, qg, p, & 1,4
rho, z, dz8w, topo, rainnc, &
rainncv, grav, itimestep, &
rhowater, rhosnow, &
snownc, snowncv, sr, &
graupelnc, graupelncv, &
ihail, ice2, &
ims,ime, jms,jme, kms,kme, & ! memory dims
its,ite, jts,jte, kts,kte ) ! tile dims
!-----------------------------------------------------------------------
! adopted from Jiun-Dar Chern's codes for Purdue Regional Model
! adopted by Jainn J. Shi, 6/10/2005
!-----------------------------------------------------------------------
IMPLICIT NONE
INTEGER, INTENT(IN ) :: ihail, ice2, &
ims,ime, jms,jme, kms,kme, &
its,ite, jts,jte, kts,kte
INTEGER, INTENT(IN ) :: itimestep
REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
INTENT(INOUT) :: qr, qi, qs, qg
REAL, DIMENSION( ims:ime , jms:jme ), &
INTENT(INOUT) :: rainnc, rainncv, &
snownc, snowncv, sr, &
graupelnc, graupelncv
REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
INTENT(IN ) :: rho, z, dz8w, p
REAL, INTENT(IN ) :: dt, grav, rhowater, rhosnow
REAL, DIMENSION( ims:ime , jms:jme ), &
INTENT(IN ) :: topo
! temperary vars
REAL, DIMENSION( kts:kte ) :: sqrhoz
REAL :: tmp1, term0
REAL :: pptrain, pptsnow, &
pptgraul, pptice
REAL, DIMENSION( kts:kte ) :: qrz, qiz, qsz, qgz, &
zz, dzw, prez, rhoz, &
orhoz
INTEGER :: k, i, j
!
REAL, DIMENSION( kts:kte ) :: vtr, vts, vtg, vti
REAL :: dtb, pi, consta, constc, gambp4, &
gamdp4, gam4pt5, gam4bbar
! Lin
!-GT REAL , PARAMETER :: xnor = 8.0e6
! REAL , PARAMETER :: xnos = 3.0e6
!-GT REAL , PARAMETER :: xnos = 1.6e7 ! Tao's value
REAL , PARAMETER :: &
! constb = 0.8, constd = 0.25, o6 = 1./6., &
constb = 0.8, constd = 0.11, o6 = 1./6., &
cdrag = 0.6
! Lin
! REAL , PARAMETER :: xnoh = 4.0e4
!-GT REAL , PARAMETER :: xnoh = 2.0e5 ! Tao's value
!-GT REAL , PARAMETER :: rhohail = 917.
! Hobbs
!-GT REAL , PARAMETER :: xnog = 4.0e6
!-GT REAL , PARAMETER :: rhograul = 400.
REAL , PARAMETER :: abar = 19.3, bbar = 0.37, &
p0 = 1.0e5
REAL , PARAMETER :: rhoe_s = 1.29
! for terminal velocity flux
INTEGER :: min_q, max_q
REAL :: t_del_tv, del_tv, flux, fluxin, fluxout ,tmpqrz
LOGICAL :: notlast
! if (itimestep.eq.1) then
! write(6, *) 'in fall_flux'
! write(6, *) 'ims=', ims, ' ime=', ime
! write(6, *) 'jms=', jms, ' jme=', jme
! write(6, *) 'kms=', kms, ' kme=', kme
! write(6, *) 'its=', its, ' ite=', ite
! write(6, *) 'jts=', jts, ' jte=', jte
! write(6, *) 'kts=', kts, ' kte=', kte
! write(6, *) 'dt=', dt
! write(6, *) 'ihail=', ihail
! write(6, *) 'ICE2=', ICE2
! write(6, *) 'dt=', dt
! endif
!-----------------------------------------------------------------------
! This program calculates precipitation fluxes due to terminal velocities.
!-----------------------------------------------------------------------
dtb=dt
pi=acos(-1.)
consta=2115.0*0.01**(1-constb)
! constc=152.93*0.01**(1-constd)
constc=78.63*0.01**(1-constd)
! Gamma function
gambp4=ggamma(constb+4.)
gamdp4=ggamma(constd+4.)
gam4pt5=ggamma(4.5)
gam4bbar=ggamma(4.+bbar)
!***********************************************************************
! Calculate precipitation fluxes due to terminal velocities.
!***********************************************************************
!
!- Calculate termianl velocity (vt?) of precipitation q?z
!- Find maximum vt? to determine the small delta t
j_loop: do j = jts, jte
i_loop: do i = its, ite
pptrain = 0.
pptsnow = 0.
pptgraul = 0.
pptice = 0.
do k = kts, kte
qrz(k)=qr(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/rhoz(k))
zz(k)=z(i,k,j)
dzw(k)=dz8w(i,k,j)
enddo !k
DO k = kts, kte
qiz(k)=qi(i,k,j)
ENDDO
DO k = kts, kte
qsz(k)=qs(i,k,j)
ENDDO
IF (ice2 .eq. 0) THEN
DO k = kts, kte
qgz(k)=qg(i,k,j)
ENDDO
ELSE
DO k = kts, kte
qgz(k)=0.
ENDDO
ENDIF
!
!-- 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 (qrz(k) .gt. 1.0e-8) then
min_q=min0(min_q,k)
max_q=max0(max_q,k)
tmp1=sqrt(pi*rhowater*xnor/rhoz(k)/qrz(k))
tmp1=sqrt(tmp1)
vtr(k)=consta*gambp4*sqrhoz(k)/tmp1**constb
vtr(k)=vtr(k)/6.
if (k .eq. 1) then
del_tv=amin1(del_tv,0.9*(zz(k)-topo(i,j))/vtr(k))
else
del_tv=amin1(del_tv,0.9*(zz(k)-zz(k-1))/vtr(k))
endif
else
vtr(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 qrz flux
! termi is the qrz flux pass in the grid box through the upper boundary
! termo is the qrz flux pass out the grid box through the lower boundary
!
fluxin=0.
do k=max_q,min_q,-1
fluxout=rhoz(k)*vtr(k)*qrz(k)
flux=(fluxin-fluxout)/rhoz(k)/dzw(k)
! tmpqrz=qrz(k)
qrz(k)=qrz(k)+del_tv*flux
qrz(k)=amax1(0.,qrz(k))
qr(i,k,j)=qrz(k)
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/rhoz(min_q-1)/dzw(min_q-1)
qr(i,min_q-1,j)=qrz(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 (qsz(k) .gt. 1.0e-8) then
min_q=min0(min_q,k)
max_q=max0(max_q,k)
tmp1=sqrt(pi*rhosnow*xnos/rhoz(k)/qsz(k))
tmp1=sqrt(tmp1)
vts(k)=constc*gamdp4*sqrhoz(k)/tmp1**constd
vts(k)=vts(k)/6.
if (k .eq. 1) then
del_tv=amin1(del_tv,0.9*(zz(k)-topo(i,j))/vts(k))
else
del_tv=amin1(del_tv,0.9*(zz(k)-zz(k-1))/vts(k))
endif
else
vts(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 qsz flux
! termi is the qsz flux pass in the grid box through the upper boundary
! termo is the qsz flux pass out the grid box through the lower boundary
!
fluxin=0.
do k=max_q,min_q,-1
fluxout=rhoz(k)*vts(k)*qsz(k)
flux=(fluxin-fluxout)/rhoz(k)/dzw(k)
qsz(k)=qsz(k)+del_tv*flux
qsz(k)=amax1(0.,qsz(k))
qs(i,k,j)=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/rhoz(min_q-1)/dzw(min_q-1)
qs(i,min_q-1,j)=qsz(min_q-1)
endif
!
else
notlast=.false.
endif
ENDDO
!
! ice2=0 --- with hail/graupel
! ice2=1 --- without hail/graupel
!
if (ice2.eq.0) then
!
!-- If IHAIL=1, use hail.
!-- If IHAIL=0, use graupel.
!
! if (ihail .eq. 1) then
! xnog = xnoh
! rhograul = rhohail
! endif
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 (qgz(k) .gt. 1.0e-8) then
if (ihail .eq. 1) then
! for hail, based on Lin et al (1983)
min_q=min0(min_q,k)
max_q=max0(max_q,k)
tmp1=sqrt(pi*rhohail*xnoh/rhoz(k)/qgz(k))
tmp1=sqrt(tmp1)
term0=sqrt(4.*grav*rhohail/3./rhoz(k)/cdrag)
vtg(k)=gam4pt5*term0*sqrt(1./tmp1)
vtg(k)=vtg(k)/6.
if (k .eq. 1) then
del_tv=amin1(del_tv,0.9*(zz(k)-topo(i,j))/vtg(k))
else
del_tv=amin1(del_tv,0.9*(zz(k)-zz(k-1))/vtg(k))
endif !k
else
! added by JJS
! for graupel, based on RH (1984)
min_q=min0(min_q,k)
max_q=max0(max_q,k)
tmp1=sqrt(pi*rhograul*xnog/rhoz(k)/qgz(k))
tmp1=sqrt(tmp1)
tmp1=tmp1**bbar
tmp1=1./tmp1
term0=abar*gam4bbar/6.
vtg(k)=term0*tmp1*(p0/prez(k))**0.4
if (k .eq. 1) then
del_tv=amin1(del_tv,0.9*(zz(k)-topo(i,j))/vtg(k))
else
del_tv=amin1(del_tv,0.9*(zz(k)-zz(k-1))/vtg(k))
endif !k
endif !ihail
else
vtg(k)=0.
endif !qgz
enddo !k
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 qgz flux
! termi is the qgz flux pass in the grid box through the upper boundary
! termo is the qgz flux pass out the grid box through the lower boundary
!
fluxin=0.
do k=max_q,min_q,-1
fluxout=rhoz(k)*vtg(k)*qgz(k)
flux=(fluxin-fluxout)/rhoz(k)/dzw(k)
qgz(k)=qgz(k)+del_tv*flux
qgz(k)=amax1(0.,qgz(k))
qg(i,k,j)=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/rhoz(min_q-1)/dzw(min_q-1)
qg(i,min_q-1,j)=qgz(min_q-1)
endif
!
else
notlast=.false.
endif
!
ENDDO
ENDIF !ice2
!
!-- 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)
vti(k)= 3.29 * (rhoz(k)* qiz(k))** 0.16 ! Heymsfield and Donner
if (k .eq. 1) then
del_tv=amin1(del_tv,0.9*(zz(k)-topo(i,j))/vti(k))
else
del_tv=amin1(del_tv,0.9*(zz(k)-zz(k-1))/vti(k))
endif
else
vti(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 qiz flux
! termi is the qiz flux pass in the grid box through the upper boundary
! termo is the qiz flux pass out the grid box through the lower boundary
!
fluxin=0.
do k=max_q,min_q,-1
fluxout=rhoz(k)*vti(k)*qiz(k)
flux=(fluxin-fluxout)/rhoz(k)/dzw(k)
qiz(k)=qiz(k)+del_tv*flux
qiz(k)=amax1(0.,qiz(k))
qi(i,k,j)=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/rhoz(min_q-1)/dzw(min_q-1)
qi(i,min_q-1,j)=qiz(min_q-1)
endif
!
else
notlast=.false.
endif
!
ENDDO !notlast
! prnc(i,j)=prnc(i,j)+pptrain
! psnowc(i,j)=psnowc(i,j)+pptsnow
! pgrauc(i,j)=pgrauc(i,j)+pptgraul
! picec(i,j)=picec(i,j)+pptice
!
! write(6,*) 'i=',i,' j=',j,' ', pptrain, pptsnow, pptgraul, pptice
! call flush(6)
snowncv(i,j) = pptsnow
snownc(i,j) = snownc(i,j) + pptsnow
graupelncv(i,j) = pptgraul
graupelnc(i,j) = graupelnc(i,j) + pptgraul
RAINNCV(i,j) = pptrain + pptsnow + pptgraul + pptice
RAINNC(i,j) = RAINNC(i,j) + pptrain + pptsnow + pptgraul + pptice
sr(i,j) = 0.
if (RAINNCV(i,j) .gt. 0.) sr(i,j) = (pptsnow + pptgraul + pptice) / RAINNCV(i,j)
ENDDO i_loop
ENDDO j_loop
! if (itimestep.eq.6480) then
! write(51,*) 'in the end of fallflux, itimestep=',itimestep
! do j=jts,jte
! do i=its,ite
! if (rainnc(i,j).gt.400.) then
! write(50,*) 'i=',i,' j=',j,' rainnc=',rainnc
! endif
! enddo
! enddo
! endif
END SUBROUTINE fall_flux
!----------------------------------------------------------------
REAL FUNCTION ggamma(X) 20,2
!----------------------------------------------------------------
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_gsfcgce: 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
!-----------------------------------------------------------------------
!c Correction of negative values
SUBROUTINE negcor ( X, rho, dz8w, & 6
ims,ime, jms,jme, kms,kme, & ! memory dims
itimestep, ics, &
its,ite, jts,jte, kts,kte ) ! tile dims
!-----------------------------------------------------------------------
REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
INTENT(INOUT) :: X
REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
INTENT(IN ) :: rho, dz8w
integer, INTENT(IN ) :: itimestep, ics
!c Local variables
! REAL, DIMENSION( kts:kte ) :: Y1, Y2
REAL :: A0, A1, A2
A1=0.
A2=0.
do k=kts,kte
do j=jts,jte
do i=its,ite
A1=A1+max(X(i,k,j), 0.)*rho(i,k,j)*dz8w(i,k,j)
A2=A2+max(-X(i,k,j), 0.)*rho(i,k,j)*dz8w(i,k,j)
enddo
enddo
enddo
! A1=0.0
! A2=0.0
! do k=kts,kte
! A1=A1+Y1(k)
! A2=A2+Y2(k)
! enddo
A0=0.0
if (A1.NE.0.0.and.A1.GT.A2) then
A0=(A1-A2)/A1
if (mod(itimestep,540).eq.0) then
if (ics.eq.1) then
write(61,*) 'kms=',kms,' kme=',kme,' kts=',kts,' kte=',kte
write(61,*) 'jms=',jms,' jme=',jme,' jts=',jts,' jte=',jte
write(61,*) 'ims=',ims,' ime=',ime,' its=',its,' ite=',ite
endif
if (ics.eq.1) then
write(61,*) 'qv timestep=',itimestep
write(61,*) ' A1=',A1,' A2=',A2,' A0=',A0
else if (ics.eq.2) then
write(61,*) 'ql timestep=',itimestep
write(61,*) ' A1=',A1,' A2=',A2,' A0=',A0
else if (ics.eq.3) then
write(61,*) 'qr timestep=',itimestep
write(61,*) ' A1=',A1,' A2=',A2,' A0=',A0
else if (ics.eq.4) then
write(61,*) 'qi timestep=',itimestep
write(61,*) ' A1=',A1,' A2=',A2,' A0=',A0
else if (ics.eq.5) then
write(61,*) 'qs timestep=',itimestep
write(61,*) ' A1=',A1,' A2=',A2,' A0=',A0
else if (ics.eq.6) then
write(61,*) 'qg timestep=',itimestep
write(61,*) ' A1=',A1,' A2=',A2,' A0=',A0
else
write(61,*) 'wrong cloud specieis number'
endif
endif
do k=kts,kte
do j=jts,jte
do i=its,ite
X(i,k,j)=A0*AMAX1(X(i,k,j), 0.0)
enddo
enddo
enddo
endif
END SUBROUTINE negcor
SUBROUTINE consat_s (ihail,itaobraun) 1,10
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
! c
! Tao, W.-K., and J. Simpson, 1989: Modeling study of a tropical c
! squall-type convective line. J. Atmos. Sci., 46, 177-202. c
! c
! Tao, W.-K., J. Simpson and M. McCumber, 1989: An ice-water c
! saturation adjustment. Mon. Wea. Rev., 117, 231-235. c
! c
! Tao, W.-K., and J. Simpson, 1993: The Goddard Cumulus Ensemble c
! Model. Part I: Model description. Terrestrial, Atmospheric and c
! Oceanic Sciences, 4, 35-72. c
! c
! Tao, W.-K., J. Simpson, D. Baker, S. Braun, M.-D. Chou, B. c
! Ferrier,D. Johnson, A. Khain, S. Lang, B. Lynn, C.-L. Shie, c
! D. Starr, C.-H. Sui, Y. Wang and P. Wetzel, 2003: Microphysics, c
! radiation and surface processes in the Goddard Cumulus Ensemble c
! (GCE) model, A Special Issue on Non-hydrostatic Mesoscale c
! Modeling, Meteorology and Atmospheric Physics, 82, 97-137. c
! c
! Lang, S., W.-K. Tao, R. Cifelli, W. Olson, J. Halverson, S. c
! Rutledge, and J. Simpson, 2007: Improving simulations of c
! convective system from TRMM LBA: Easterly and Westerly regimes. c
! J. Atmos. Sci., 64, 1141-1164. c
! c
! Coded by Tao (1989-2003), modified by S. Lang (2006/07) c
! c
! Implemented into WRF by Roger Shi 2006/2007 c
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
! itaobraun=0 ! see Tao and Simpson (1993)
! itaobraun=1 ! see Tao et al. (2003)
integer :: itaobraun
real :: cn0
!JJS 1/3/2008 vvvvv
!JJS the following common blocks have been moved to the top of
!JJS module_mp_gsfcgce_driver_instat.F
!
! real, dimension (1:31) :: a1, a2
! data a1/.7939e-7,.7841e-6,.3369e-5,.4336e-5,.5285e-5,.3728e-5, &
! .1852e-5,.2991e-6,.4248e-6,.7434e-6,.1812e-5,.4394e-5,.9145e-5, &
! .1725e-4,.3348e-4,.1725e-4,.9175e-5,.4412e-5,.2252e-5,.9115e-6, &
! .4876e-6,.3473e-6,.4758e-6,.6306e-6,.8573e-6,.7868e-6,.7192e-6, &
! .6513e-6,.5956e-6,.5333e-6,.4834e-6/
! data a2/.4006,.4831,.5320,.5307,.5319,.5249,.4888,.3894,.4047, &
! .4318,.4771,.5183,.5463,.5651,.5813,.5655,.5478,.5203,.4906, &
! .4447,.4126,.3960,.4149,.4320,.4506,.4483,.4460,.4433,.4413, &
! .4382,.4361/
!JJS 1/3/2008 ^^^^^
! ******************************************************************
!JJS
al = 2.5e10
cp = 1.004e7
rd1 = 1.e-3
rd2 = 2.2
!JJS
cpi=4.*atan(1.)
cpi2=cpi*cpi
grvt=980.
cd1=6.e-1
cd2=4.*grvt/(3.*cd1)
tca=2.43e3
dwv=.226
dva=1.718e-4
amw=18.016
ars=8.314e7
scv=2.2904487
t0=273.16
t00=238.16
alv=2.5e10
alf=3.336e9
als=2.8336e10
avc=alv/cp
afc=alf/cp
asc=als/cp
rw=4.615e6
cw=4.187e7
ci=2.093e7
c76=7.66
c358=35.86
c172=17.26939
c409=4098.026
c218=21.87456
c580=5807.695
c610=6.1078e3
c149=1.496286e-5
c879=8.794142
c141=1.4144354e7
!*** DEFINE THE COEFFICIENTS USED IN TERMINAL VELOCITY
!*** DEFINE THE DENSITY AND SIZE DISTRIBUTION OF PRECIPITATION
!********** HAIL OR GRAUPEL PARAMETERS **********
if(ihail .eq. 1) then
roqg=.9
ag=sqrt(cd2*roqg)
bg=.5
tng=.002
else
roqg=.4
ag=351.2
! AG=372.3 ! if ice913=1 6/15/02 tao's
bg=.37
tng=.04
endif
!********** SNOW PARAMETERS **********
!ccshie 6/15/02 tao's
! TNS=1.
! TNS=.08 ! if ice913=1, tao's
tns=.16 ! if ice913=0, tao's
roqs=.1
! AS=152.93
as=78.63
! BS=.25
bs=.11
!********** RAIN PARAMETERS **********
aw=2115.
bw=.8
roqr=1.
tnw=.08
!*****************************************************************
bgh=.5*bg
bsh=.5*bs
bwh=.5*bw
bgq=.25*bg
bsq=.25*bs
bwq=.25*bw
!**********GAMMA FUNCTION CALCULATIONS*************
ga3b = gammagce(3.+bw)
ga4b = gammagce(4.+bw)
ga6b = gammagce(6.+bw)
ga5bh = gammagce((5.+bw)/2.)
ga3g = gammagce(3.+bg)
ga4g = gammagce(4.+bg)
ga5gh = gammagce((5.+bg)/2.)
ga3d = gammagce(3.+bs)
ga4d = gammagce(4.+bs)
ga5dh = gammagce((5.+bs)/2.)
!CCCCC LIN ET AL., 1983 OR LORD ET AL., 1984 CCCCCCCCCCCCCCCCC
ac1=aw
!JJS
ac2=ag
ac3=as
!JJS
bc1=bw
cc1=as
dc1=bs
zrc=(cpi*roqr*tnw)**0.25
zsc=(cpi*roqs*tns)**0.25
zgc=(cpi*roqg*tng)**0.25
vrc=aw*ga4b/(6.*zrc**bw)
vsc=as*ga4d/(6.*zsc**bs)
vgc=ag*ga4g/(6.*zgc**bg)
! ****************************
! RN1=1.E-3
rn1=9.4e-15 ! 6/15/02 tao's
bnd1=6.e-4
rn2=1.e-3
! BND2=1.25E-3
! BND2=1.5E-3 ! if ice913=1 6/15/02 tao's
bnd2=2.0e-3 ! if ice913=0 6/15/02 tao's
rn3=.25*cpi*tns*cc1*ga3d
esw=1.
rn4=.25*cpi*esw*tns*cc1*ga3d
! ERI=1.
eri=.1 ! 6/17/02 tao's ice913=0 (not 1)
rn5=.25*cpi*eri*tnw*ac1*ga3b
! AMI=1./(24.*4.19E-10)
ami=1./(24.*6.e-9) ! 6/15/02 tao's
rn6=cpi2*eri*tnw*ac1*roqr*ga6b*ami
! ESR=1. ! also if ice913=1 for tao's
esr=.5 ! 6/15/02 for ice913=0 tao's
rn7=cpi2*esr*tnw*tns*roqs
esr=1.
rn8=cpi2*esr*tnw*tns*roqr
rn9=cpi2*tns*tng*roqs
rn10=2.*cpi*tns
rn101=.31*ga5dh*sqrt(cc1)
rn10a=als*als/rw
!JJS
rn10b=alv/tca
rn10c=ars/(dwv*amw)
!JJS
rn11=2.*cpi*tns/alf
rn11a=cw/alf
! AMI50=1.51e-7
ami50=3.84e-6 ! 6/15/02 tao's
! AMI40=2.41e-8
ami40=3.08e-8 ! 6/15/02 tao's
eiw=1.
! UI50=20.
ui50=100. ! 6/15/02 tao's
ri50=2.*5.e-3
cmn=1.05e-15
rn12=cpi*eiw*ui50*ri50**2
do 10 k=1,31
y1=1.-aa2(k)
rn13(k)=aa1(k)*y1/(ami50**y1-ami40**y1)
rn12a(k)=rn13(k)/ami50
rn12b(k)=aa1(k)*ami50**aa2(k)
rn25a(k)=aa1(k)*cmn**aa2(k)
10 continue
egw=1.
rn14=.25*cpi*egw*tng*ga3g*ag
egi=.1
rn15=.25*cpi*egi*tng*ga3g*ag
egi=1.
rn15a=.25*cpi*egi*tng*ga3g*ag
egr=1.
rn16=cpi2*egr*tng*tnw*roqr
rn17=2.*cpi*tng
rn17a=.31*ga5gh*sqrt(ag)
rn17b=cw-ci
rn17c=cw
apri=.66
bpri=1.e-4
bpri=0.5*bpri ! 6/17/02 tao's
rn18=20.*cpi2*bpri*tnw*roqr
rn18a=apri
rn19=2.*cpi*tng/alf
rn19a=.31*ga5gh*sqrt(ag)
rn19b=cw/alf
!
rnn191=.78
rnn192=.31*ga5gh*sqrt(ac2/dva)
!
rn20=2.*cpi*tng
rn20a=als*als/rw
rn20b=.31*ga5gh*sqrt(ag)
bnd3=2.e-3
rn21=1.e3*1.569e-12/0.15
erw=1.
rn22=.25*cpi*erw*ac1*tnw*ga3b
rn23=2.*cpi*tnw
rn23a=.31*ga5bh*sqrt(ac1)
rn23b=alv*alv/rw
!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!cc
!cc "c0" in routine "consat" (2d), "consatrh" (3d)
!cc if ( itaobraun.eq.1 ) --> betah=0.5*beta=-.46*0.5=-0.23; cn0=1.e-6
!cc if ( itaobraun.eq.0 ) --> betah=0.5*beta=-.6*0.5=-0.30; cn0=1.e-8
if (itaobraun .eq. 0) then
cn0=1.e-8
beta=-.6
elseif (itaobraun .eq. 1) then
cn0=1.e-6
beta=-.46
endif
!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
! CN0=1.E-6
! CN0=1.E-8 ! 6/15/02 tao's
! BETA=-.46
! BETA=-.6 ! 6/15/02 tao's
rn25=cn0
rn30a=alv*als*amw/(tca*ars)
rn30b=alv/tca
rn30c=ars/(dwv*amw)
rn31=1.e-17
rn32=4.*51.545e-4
!
rn30=2.*cpi*tng
rnn30a=alv*alv*amw/(tca*ars)
!
rn33=4.*tns
rn331=.65
rn332=.44*sqrt(ac3/dva)*ga5dh
!
return
END SUBROUTINE consat_s
SUBROUTINE saticel_s (dt, ihail, itaobraun, ice2, istatmin, & 1,1
new_ice_sat, id, &
ptwrf, qvwrf, qlwrf, qrwrf, &
qiwrf, qswrf, qgwrf, &
rho_mks, pi_mks, p0_mks,itimestep, &
refl_10cm, diagflag, do_radar_ref, & ! GT added for reflectivity calcs
! refl_10cm, grid_clock, grid_alarms, & ! GT added for reflectivity calcs
ids,ide, jds,jde, kds,kde, &
ims,ime, jms,jme, kms,kme, &
its,ite, jts,jte, kts,kte &
)
IMPLICIT NONE
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
! c
! Tao, W.-K., and J. Simpson, 1989: Modeling study of a tropical c
! squall-type convective line. J. Atmos. Sci., 46, 177-202. c
! c
! Tao, W.-K., J. Simpson and M. McCumber, 1989: An ice-water c
! saturation adjustment. Mon. Wea. Rev., 117, 231-235. c
! c
! Tao, W.-K., and J. Simpson, 1993: The Goddard Cumulus Ensemble c
! Model. Part I: Model description. Terrestrial, Atmospheric and c
! Oceanic Sciences, 4, 35-72. c
! c
! Tao, W.-K., J. Simpson, D. Baker, S. Braun, M.-D. Chou, B. c
! Ferrier,D. Johnson, A. Khain, S. Lang, B. Lynn, C.-L. Shie, c
! D. Starr, C.-H. Sui, Y. Wang and P. Wetzel, 2003: Microphysics, c
! radiation and surface processes in the Goddard Cumulus Ensemble c
! (GCE) model, A Special Issue on Non-hydrostatic Mesoscale c
! Modeling, Meteorology and Atmospheric Physics, 82, 97-137. c
! c
! Lang, S., W.-K. Tao, R. Cifelli, W. Olson, J. Halverson, S. c
! Rutledge, and J. Simpson, 2007: Improving simulations of c
! convective system from TRMM LBA: Easterly and Westerly regimes. c
! J. Atmos. Sci., 64, 1141-1164. c
! c
! Tao, W.-K., J. J. Shi, S. Lang, C. Peters-Lidard, A. Hou, S. c
! Braun, and J. Simpson, 2007: New, improved bulk-microphysical c
! schemes for studying precipitation processes in WRF. Part I: c
! Comparisons with other schemes. to appear on Mon. Wea. Rev. C
! c
! Coded by Tao (1989-2003), modified by S. Lang (2006/07) c
! c
! Implemented into WRF by Roger Shi 2006/2007 c
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!
! COMPUTE ICE PHASE MICROPHYSICS AND SATURATION PROCESSES
!
integer, parameter :: nt=2880, nt2=2*nt
!cc using scott braun's way for pint, pidep computations
integer :: itaobraun,ice2,ihail,new_ice_sat,id,istatmin
integer :: itimestep
real :: tairccri, cn0, dt
!cc
!JJS common/bxyz/ n,isec,nran,kt1,kt2
!JJS common/option/ lipps,ijkadv,istatmin,iwater,itoga,imlifting,lin,
!JJS 1 irf,iadvh,irfg,ismg,id
!JJS common/timestat/ ndt_stat,idq
!JJS common/iice/ new_ice_sat
!JJS common/bt/ dt,d2t,rijl2,dts,f5,rd1,rd2,bound,al,cp,ra,ck,ce,eps,
!JJS 1 psfc,fcor,sec,aminut,rdt
!JJS the following common blocks have been moved to the top of
!JJS module_mp_gsfcgce_driver_instat.F
! common/bt/ rd1,rd2,al,cp
!
!
! common/bterv/ zrc,zgc,zsc,vrc,vgc,vsc
! common/size/ tnw,tns,tng,roqs,roqg,roqr
! common/cont/ c38,c358,c610,c149,c879,c172,c409,c76,c218,c580,c141
! common/b3cs/ ag,bg,as,bs,aw,bw,bgh,bgq,bsh,bsq,bwh,bwq
! common/bsnw/ alv,alf,als,t0,t00,avc,afc,asc,rn1,bnd1,rn2,bnd2, &
! rn3,rn4,rn5,rn6,rn7,rn8,rn9,rn10,rn101,rn10a,rn11,rn11a, &
! rn12,rn12a(31),rn12b(31),rn13(31),rn14,rn15,rn15a,rn16,rn17, &
! rn17a,rn17b,rn17c,rn18,rn18a,rn19,rn19a,rn19b,rn20,rn20a,rn20b, &
! bnd3,rn21,rn22,rn23,rn23a,rn23b,rn25,rn25a(31),rn30a,rn30b, &
! rn30c,rn31,beta,rn32
! common/rsnw1/ rn10b,rn10c,rnn191,rnn192,rn30,rnn30a,rn33,rn331, &
! rn332
!JJS
integer ids,ide,jds,jde,kds,kde
integer ims,ime,jms,jme,kms,kme
integer its,ite,jts,jte,kts,kte
integer i,j,k, kp
real :: a0 ,a1 ,a2 ,afcp ,alvr ,ami100 ,ami40 ,ami50 ,ascp ,avcp ,betah &
,bg3 ,bgh5 ,bs3 ,bs6 ,bsh5 ,bw3 ,bw6 ,bwh5 ,cmin ,cmin1 ,cmin2 ,cp409 &
,cp580 ,cs580 ,cv409 ,d2t ,del ,dwvp ,ee1 ,ee2 ,f00 ,f2 ,f3 ,ft ,fv0 ,fvs &
,pi0 ,pir ,pr0 ,qb0 ,r00 ,r0s ,r101f ,r10ar ,r10t ,r11at ,r11rt ,r12r ,r14f &
,r14r ,r15af ,r15ar ,r15f ,r15r ,r16r ,r17aq ,r17as ,r17r ,r18r ,r19aq ,r19as &
,r19bt ,r19rt ,r20bq ,r20bs ,r20t ,r22f ,r23af ,r23br ,r23t ,r25a ,r25rt ,r2ice &
,r31r ,r32rt ,r3f ,r4f ,r5f ,r6f ,r7r ,r8r ,r9r ,r_nci ,rft ,rijl2 ,rp0 ,rr0 &
,rrq ,rrs ,rt0 ,scc ,sccc ,sddd ,see ,seee ,sfff ,smmm ,ssss ,tb0 ,temp ,ucog &
,ucor ,ucos ,uwet ,vgcf ,vgcr ,vrcf ,vscf ,zgr ,zrr ,zsr
real, dimension (its:ite,jts:jte,kts:kte) :: fv
real, dimension (its:ite,jts:jte,kts:kte) :: dpt, dqv
real, dimension (its:ite,jts:jte,kts:kte) :: qcl, qrn, &
qci, qcs, qcg
!JJS 10/16/06 vvvv
! real dpt1(ims:ime,jms:jme,kms:kme)
! real dqv1(ims:ime,jms:jme,kms:kme),
! 1 qcl1(ims:ime,jms:jme,kms:kme)
! real qrn1(ims:ime,jms:jme,kms:kme),
! 1 qci1(ims:ime,jms:jme,kms:kme)
! real qcs1(ims:ime,jms:jme,kms:kme),
! 1 qcg1(ims:ime,jms:jme,kms:kme)
!JJS 10/16/06 ^^^^
!JJS
real, dimension (ims:ime, kms:kme, jms:jme) :: ptwrf, qvwrf
real, dimension (ims:ime, kms:kme, jms:jme) :: qlwrf, qrwrf, &
qiwrf, qswrf, qgwrf
!JJS 10/16/06 vvvv
! real ptwrfold(ims:ime, kms:kme, jms:jme)
! real qvwrfold(ims:ime, kms:kme, jms:jme),
! 1 qlwrfold(ims:ime, kms:kme, jms:jme)
! real qrwrfold(ims:ime, kms:kme, jms:jme),
! 1 qiwrfold(ims:ime, kms:kme, jms:jme)
! real qswrfold(ims:ime, kms:kme, jms:jme),
! 1 qgwrfold(ims:ime, kms:kme, jms:jme)
!JJS 10/16/06 ^^^^
!JJS in MKS
real, dimension (ims:ime, kms:kme, jms:jme) :: rho_mks
real, dimension (ims:ime, kms:kme, jms:jme) :: pi_mks
real, dimension (ims:ime, kms:kme, jms:jme) :: p0_mks
!JJS
! real, dimension (its:ite,jts:jte,kts:kte) :: ww1
! real, dimension (its:ite,jts:jte,kts:kte) :: rsw
! real, dimension (its:ite,jts:jte,kts:kte) :: rlw
!JJS COMMON /BADV/
real, dimension (its:ite,jts:jte) :: &
vg, zg, &
ps, pg, &
prn, psn, &
pwacs, wgacr, &
pidep, pint, &
qsi, ssi, &
esi, esw, &
qsw, pr, &
ssw, pihom, &
pidw, pimlt, &
psaut, qracs, &
psaci, psacw, &
qsacw, praci, &
pmlts, pmltg, &
asss, y1, y2
!JJS Y2(its:ite,jts:jte), DDE(NB)
!JJS COMMON/BSAT/
real, dimension (its:ite,jts:jte) :: &
praut, pracw, &
psfw, psfi, &
dgacs, dgacw, &
dgaci, dgacr, &
pgacs, wgacs, &
qgacw, wgaci, &
qgacr, pgwet, &
pgaut, pracs, &
psacr, qsacr, &
pgfr, psmlt, &
pgmlt, psdep, &
pgdep, piacr, &
y5, scv, &
tca, dwv, &
egs, y3, &
y4, ddb
!JJS COMMON/BSAT1/
real, dimension (its:ite,jts:jte) :: &
pt, qv, &
qc, qr, &
qi, qs, &
qg, tair, &
tairc, rtair, &
dep, dd, &
dd1, qvs, &
dm, rq, &
rsub1, col, &
cnd, ern, &
dlt1, dlt2, &
dlt3, dlt4, &
zr, vr, &
zs, vs, &
pssub, &
pgsub, dda
!JJS COMMON/B5/
real, dimension (its:ite,jts:jte,kts:kte) :: rho
real, dimension (kts:kte) :: &
tb, qb, rho1, &
ta, qa, ta1, qa1, &
coef, z1, z2, z3, &
am, am1, ub, vb, &
wb, ub1, vb1, rrho, &
rrho1, wbx
!JJS COMMON/B6/
real, dimension (its:ite,jts:jte,kts:kte) :: p0, pi, f0
real, dimension (kts:kte) :: &
fd, fe, &
st, sv, &
sq, sc, &
se, sqa
!JJS COMMON/BRH1/
real, dimension (kts:kte) :: &
srro, qrro, sqc, sqr, &
sqi, sqs, sqg, stqc, &
stqr, stqi, stqs, stqg
real, dimension (nt) :: &
tqc, tqr, tqi, tqs, tqg
!JJS common/bls/ y0(nx,ny),ts0new(nx,ny),qss0new(nx,ny)
!JJS COMMON/BLS/
real, dimension (ims:ime,jms:jme) :: &
y0, ts0, qss0
!JJS COMMON/BI/ IT(its:ite,jts:jte), ICS(its:ite,jts:jte,4)
integer, dimension (its:ite,jts:jte) :: it
integer, dimension (its:ite,jts:jte, 4) :: ics
integer :: i24h
integer :: iwarm
real :: r2is, r2ig
!JJS COMMON/MICRO/
! real, dimension (ims:ime,kms:kme,jms:jme) :: dbz
!23456789012345678901234567890123456789012345678901234567890123456789012
!
!JJS 1/3/2008 vvvvv
!JJS the following common blocks have been moved to the top of
!JJS module_mp_gsfcgce_driver.F
! real, dimension (31) :: aa1, aa2
! data aa1/.7939e-7, .7841e-6, .3369e-5, .4336e-5, .5285e-5, &
! .3728e-5, .1852e-5, .2991e-6, .4248e-6, .7434e-6, &
! .1812e-5, .4394e-5, .9145e-5, .1725e-4, .3348e-4, &
! .1725e-4, .9175e-5, .4412e-5, .2252e-5, .9115e-6, &
! .4876e-6, .3473e-6, .4758e-6, .6306e-6, .8573e-6, &
! .7868e-6, .7192e-6, .6513e-6, .5956e-6, .5333e-6, &
! .4834e-6/
! data aa2/.4006, .4831, .5320, .5307, .5319, &
! .5249, .4888, .3894, .4047, .4318, &
! .4771, .5183, .5463, .5651, .5813, &
! .5655, .5478, .5203, .4906, .4447, &
! .4126, .3960, .4149, .4320, .4506, &
! .4483, .4460, .4433, .4413, .4382, &
! .4361/
!JJS 1/3/2008 ^^^^^
!+---+-----------------------------------------------------------------+
REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT):: refl_10cm ! GT
! TYPE (WRFU_Clock):: grid_clock
! TYPE (WRFU_Alarm), POINTER:: grid_alarms(:)
REAL, DIMENSION(kts:kte):: qv1d, t1d, p1d, qr1d, qs1d, qg1d, dBZ
LOGICAL, OPTIONAL, INTENT(IN) :: diagflag
INTEGER, OPTIONAL, INTENT(IN) :: do_radar_ref
!+---+-----------------------------------------------------------------+
!
!jm 20090220 save
! i24h=nint(86400./dt)
! if (mod(itimestep,i24h).eq.1) then
! write(6, *) 'ims=', ims, ' ime=', ime
! write(6, *) 'jms=', jms, ' jme=', jme
! write(6, *) 'kms=', kms, ' kme=', kme
! write(6, *) 'its=', its, ' ite=', ite
! write(6, *) 'jts=', jts, ' jte=', jte
! write(6, *) 'kts=', kts, ' kte=', kte
! write(6, *) ' ihail=', ihail
! write(6, *) 'itaobraun=',itaobraun
! write(6, *) ' ice2=', ICE2
! write(6, *) 'istatmin=',istatmin
! write(6, *) 'new_ice_sat=', new_ice_sat
! write(6, *) 'id=', id
! write(6, *) 'dt=', dt
! endif
!JJS convert from mks to cgs, and move from WRF grid to GCE grid
do k=kts,kte
do j=jts,jte
do i=its,ite
rho(i,j,k)=rho_mks(i,k,j)*0.001
p0(i,j,k)=p0_mks(i,k,j)*10.0
pi(i,j,k)=pi_mks(i,k,j)
dpt(i,j,k)=ptwrf(i,k,j)
dqv(i,j,k)=qvwrf(i,k,j)
qcl(i,j,k)=qlwrf(i,k,j)
qrn(i,j,k)=qrwrf(i,k,j)
qci(i,j,k)=qiwrf(i,k,j)
qcs(i,j,k)=qswrf(i,k,j)
qcg(i,j,k)=qgwrf(i,k,j)
!JJS 10/16/06 vvvv
! dpt1(i,j,k)=ptwrfold(i,k,j)
! dqv1(i,j,k)=qvwrfold(i,k,j)
! qcl1(i,j,k)=qlwrfold(i,k,j)
! qrn1(i,j,k)=qrwrfold(i,k,j)
! qci1(i,j,k)=qiwrfold(i,k,j)
! qcs1(i,j,k)=qswrfold(i,k,j)
! qcg1(i,j,k)=qgwrfold(i,k,j)
!JJS 10/16/06 ^^^^
enddo !i
enddo !j
enddo !k
do k=kts,kte
do j=jts,jte
do i=its,ite
fv(i,j,k)=sqrt(rho(i,j,2)/rho(i,j,k))
enddo !i
enddo !j
enddo !k
!JJS
!
! ****** THREE CLASSES OF ICE-PHASE (LIN ET AL, 1983) *********
!JJS D22T=D2T
!JJS IF(IJKADV .EQ. 0) THEN
!JJS D2T=D2T
!JJS ELSE
d2t=dt
!JJS ENDIF
!
! itaobraun=0 ! original pint and pidep & see Tao and Simpson 1993
itaobraun=1 ! see Tao et al. (2003)
!
if ( itaobraun.eq.0 ) then
cn0=1.e-8
!c beta=-.6
elseif ( itaobraun.eq.1 ) then
cn0=1.e-6
! cn0=1.e-8 ! special
!c beta=-.46
endif
!C TAO 2007 START
! ICE2=0 ! default, 3ice with loud ice, snow and graupel
! r2is=1., r2ig=1.
! ICE2=1 ! 2ice with cloud ice and snow (no graupel) - r2iceg=1, r2ice=0.
! r2is=1., r2ig=0.
! ICE2=2 ! 2ice with cloud ice and graupel (no snow) - r2ice=1, r2iceg=0.
! r2is=0., r2ig=1.
!c
! r2ice=1.
! r2iceg=1.
r2ig=1.
r2is=1.
if (ice2 .eq. 1) then
! r2ice=0.
! r2iceg=1.
r2ig=0.
r2is=1.
endif
if (ice2 .eq. 2) then
! r2ice=1.
! r2iceg=0.
r2ig=1.
r2is=0.
endif
!C TAO 2007 END
!JJS 10/7/2008
! ICE2=3 ! no ice, warm rain only
iwarm = 0
if (ice2 .eq. 3 ) iwarm = 1
cmin=1.e-19
cmin1=1.e-20
cmin2=1.e-12
ucor=3071.29/tnw**0.75
ucos=687.97*roqs**0.25/tns**0.75
ucog=687.97*roqg**0.25/tng**0.75
uwet=4.464**0.95
rijl2 = 1. / (ide-ids) / (jde-jds)
!JJScap $doacross local(j,i)
!JJS DO 1 J=1,JMAX
!JJS DO 1 I=1,IMAX
do j=jts,jte
do i=its,ite
it(i,j)=1
enddo
enddo
f2=rd1*d2t
f3=rd2*d2t
ft=dt/d2t
rft=rijl2*ft
a0=.5*istatmin*rijl2
rt0=1./(t0-t00)
bw3=bw+3.
bs3=bs+3.
bg3=bg+3.
bsh5=2.5+bsh
bgh5=2.5+bgh
bwh5=2.5+bwh
bw6=bw+6.
bs6=bs+6.
betah=.5*beta
r10t=rn10*d2t
r11at=rn11a*d2t
r19bt=rn19b*d2t
r20t=-rn20*d2t
r23t=-rn23*d2t
r25a=rn25
! ami50 for use in PINT
ami50=3.76e-8
ami100=1.51e-7
ami40=2.41e-8
!C ******************************************************************
!JJS DO 1000 K=2,kles
do 1000 k=kts,kte
kp=k+1
!JJS tb0=ta1(k)
!JJS qb0=qa1(k)
tb0=0.
qb0=0.
do 2000 j=jts,jte
do 2000 i=its,ite
rp0=3.799052e3/p0(i,j,k)
pi0=pi(i,j,k)
pir=1./(pi(i,j,k))
pr0=1./p0(i,j,k)
r00=rho(i,j,k)
r0s=sqrt(rho(i,j,k))
!JJS RR0=RRHO(K)
rr0=1./rho(i,j,k)
!JJS RRS=SRRO(K)
rrs=sqrt(rr0)
!JJS RRQ=QRRO(K)
rrq=sqrt(rrs)
f0(i,j,k)=al/cp/pi(i,j,k)
f00=f0(i,j,k)
fv0=fv(i,j,k)
fvs=sqrt(fv(i,j,k))
zrr=1.e5*zrc*rrq
zsr=1.e5*zsc*rrq
zgr=1.e5*zgc*rrq
cp409=c409*pi0
cv409=c409*avc
cp580=c580*pi0
cs580=c580*asc
alvr=r00*alv
afcp=afc*pir
avcp=avc*pir
ascp=asc*pir
vrcf=vrc*fv0
vscf=vsc*fv0
vgcf=vgc*fv0
vgcr=vgc*rrs
dwvp=c879*pr0
r3f=rn3*fv0
r4f=rn4*fv0
r5f=rn5*fv0
r6f=rn6*fv0
r7r=rn7*rr0
r8r=rn8*rr0
r9r=rn9*rr0
r101f=rn101*fvs
r10ar=rn10a*r00
r11rt=rn11*rr0*d2t
r12r=rn12*r00
r14r=rn14*rrs
r14f=rn14*fv0
r15r=rn15*rrs
r15ar=rn15a*rrs
r15f=rn15*fv0
r15af=rn15a*fv0
r16r=rn16*rr0
r17r=rn17*rr0
r17aq=rn17a*rrq
r17as=rn17a*fvs
r18r=rn18*rr0
r19rt=rn19*rr0*d2t
r19aq=rn19a*rrq
r19as=rn19a*fvs
r20bq=rn20b*rrq
r20bs=rn20b*fvs
r22f=rn22*fv0
r23af=rn23a*fvs
r23br=rn23b*r00
r25rt=rn25*rr0*d2t
r31r=rn31*rr0
r32rt=rn32*d2t*rrs
!JJS DO 100 J=3,JLES
!JJS DO 100 I=3,ILES
pt(i,j)=dpt(i,j,k)
qv(i,j)=dqv(i,j,k)
qc(i,j)=qcl(i,j,k)
qr(i,j)=qrn(i,j,k)
qi(i,j)=qci(i,j,k)
qs(i,j)=qcs(i,j,k)
qg(i,j)=qcg(i,j,k)
! IF (QV(I,J)+QB0 .LE. 0.) QV(I,J)=-QB0
if (qc(i,j) .le. cmin1) qc(i,j)=0.0
if (qr(i,j) .le. cmin1) qr(i,j)=0.0
if (qi(i,j) .le. cmin1) qi(i,j)=0.0
if (qs(i,j) .le. cmin1) qs(i,j)=0.0
if (qg(i,j) .le. cmin1) qg(i,j)=0.0
tair(i,j)=(pt(i,j)+tb0)*pi0
tairc(i,j)=tair(i,j)-t0
zr(i,j)=zrr
zs(i,j)=zsr
zg(i,j)=zgr
vr(i,j)=0.0
vs(i,j)=0.0
vg(i,j)=0.0
!JJS 10/7/2008 vvvvv
IF (IWARM .EQ. 1) THEN
!JJS for calculating processes related to warm rain only
qi(i,j)=0.0
qs(i,j)=0.0
qg(i,j)=0.0
dep(i,j)=0.
pint(i,j)=0.
psdep(i,j)=0.
pgdep(i,j)=0.
dd1(i,j)=0.
pgsub(i,j)=0.
psmlt(i,j)=0.
pgmlt(i,j)=0.
pimlt(i,j)=0.
psacw(i,j)=0.
piacr(i,j)=0.
psfw(i,j)=0.
pgfr(i,j)=0.
dgacw(i,j)=0.
dgacr(i,j)=0.
psacr(i,j)=0.
wgacr(i,j)=0.
pihom(i,j)=0.
pidw(i,j)=0.
if (qr(i,j) .gt. cmin1) then
dd(i,j)=r00*qr(i,j)
y1(i,j)=dd(i,j)**.25
zr(i,j)=zrc/y1(i,j)
vr(i,j)=max(vrcf*dd(i,j)**bwq, 0.)
endif
!* 21 * PRAUT AUTOCONVERSION OF QC TO QR **21**
!* 22 * PRACW : ACCRETION OF QC BY QR **22**
pracw(i,j)=0.
praut(i,j)=0.0
pracw(i,j)=r22f*qc(i,j)/zr(i,j)**bw3
y1(i,j)=qc(i,j)-bnd3
if (y1(i,j).gt.0.0) then
praut(i,j)=r00*y1(i,j)*y1(i,j)/(1.2e-4+rn21/y1(i,j))
endif
!C******** HANDLING THE NEGATIVE CLOUD WATER (QC) ******************
Y1(I,J)=QC(I,J)/D2T
PRAUT(I,J)=MIN(Y1(I,J), PRAUT(I,J))
PRACW(I,J)=MIN(Y1(I,J), PRACW(I,J))
Y1(I,J)=(PRAUT(I,J)+PRACW(I,J))*D2T
if (qc(i,j) .lt. y1(i,j) .and. y1(i,j) .ge. cmin2) then
y2(i,j)=qc(i,j)/(y1(i,j)+cmin2)
praut(i,j)=praut(i,j)*y2(i,j)
pracw(i,j)=pracw(i,j)*y2(i,j)
qc(i,j)=0.0
else
qc(i,j)=qc(i,j)-y1(i,j)
endif
PR(I,J)=(PRAUT(I,J)+PRACW(I,J))*D2T
QR(I,J)=QR(I,J)+PR(I,J)
!***** TAO ET AL (1989) SATURATION TECHNIQUE ***********************
cnd(i,j)=0.0
tair(i,j)=(pt(i,j)+tb0)*pi0
y1(i,j)=1./(tair(i,j)-c358)
qsw(i,j)=rp0*exp(c172-c409*y1(i,j))
dd(i,j)=cp409*y1(i,j)*y1(i,j)
dm(i,j)=qv(i,j)+qb0-qsw(i,j)
cnd(i,j)=dm(i,j)/(1.+avcp*dd(i,j)*qsw(i,j))
!c ****** condensation or evaporation of qc ******
cnd(i,j)=max(-qc(i,j), cnd(i,j))
pt(i,j)=pt(i,j)+avcp*cnd(i,j)
qv(i,j)=qv(i,j)-cnd(i,j)
qc(i,j)=qc(i,j)+cnd(i,j)
!C ****** EVAPORATION ******
!* 23 * ERN : EVAPORATION OF QR (SUBSATURATION) **23**
ern(i,j)=0.0
if(qr(i,j).gt.0.0) then
tair(i,j)=(pt(i,j)+tb0)*pi0
rtair(i,j)=1./(tair(i,j)-c358)
qsw(i,j)=rp0*exp(c172-c409*rtair(i,j))
ssw(i,j)=(qv(i,j)+qb0)/qsw(i,j)-1.0
dm(i,j)=qv(i,j)+qb0-qsw(i,j)
rsub1(i,j)=cv409*qsw(i,j)*rtair(i,j)*rtair(i,j)
dd1(i,j)=max(-dm(i,j)/(1.+rsub1(i,j)), 0.0)
y1(i,j)=.78/zr(i,j)**2+r23af*scv(i,j)/zr(i,j)**bwh5
y2(i,j)=r23br/(tca(i,j)*tair(i,j)**2)+1./(dwv(i,j) &
*qsw(i,j))
!cccc
ern(i,j)=r23t*ssw(i,j)*y1(i,j)/y2(i,j)
ern(i,j)=min(dd1(i,j),qr(i,j),max(ern(i,j),0.))
pt(i,j)=pt(i,j)-avcp*ern(i,j)
qv(i,j)=qv(i,j)+ern(i,j)
qr(i,j)=qr(i,j)-ern(i,j)
endif
ELSE ! part of if (iwarm.eq.1) then
!JJS 10/7/2008 ^^^^^
!JJS for calculating processes related to both ice and warm rain
! *** COMPUTE ZR,ZS,ZG,VR,VS,VG *****************************
if (qr(i,j) .gt. cmin1) then
dd(i,j)=r00*qr(i,j)
y1(i,j)=dd(i,j)**.25
zr(i,j)=zrc/y1(i,j)
vr(i,j)=max(vrcf*dd(i,j)**bwq, 0.)
endif
if (qs(i,j) .gt. cmin1) then
dd(i,j)=r00*qs(i,j)
y1(i,j)=dd(i,j)**.25
zs(i,j)=zsc/y1(i,j)
vs(i,j)=max(vscf*dd(i,j)**bsq, 0.)
endif
if (qg(i,j) .gt. cmin1) then
dd(i,j)=r00*qg(i,j)
y1(i,j)=dd(i,j)**.25
zg(i,j)=zgc/y1(i,j)
if(ihail .eq. 1) then
vg(i,j)=max(vgcr*dd(i,j)**bgq, 0.)
else
vg(i,j)=max(vgcf*dd(i,j)**bgq, 0.)
endif
endif
if (qr(i,j) .le. cmin2) vr(i,j)=0.0
if (qs(i,j) .le. cmin2) vs(i,j)=0.0
if (qg(i,j) .le. cmin2) vg(i,j)=0.0
! ******************************************************************
! *** Y1 : DYNAMIC VISCOSITY OF AIR (U)
! *** DWV : DIFFUSIVITY OF WATER VAPOR IN AIR (PI)
! *** TCA : THERMAL CONDUCTIVITY OF AIR (KA)
! *** Y2 : KINETIC VISCOSITY (V)
y1(i,j)=c149*tair(i,j)**1.5/(tair(i,j)+120.)
dwv(i,j)=dwvp*tair(i,j)**1.81
tca(i,j)=c141*y1(i,j)
scv(i,j)=1./((rr0*y1(i,j))**.1666667*dwv(i,j)**.3333333)
!JJS 100 CONTINUE
!* 1 * PSAUT : AUTOCONVERSION OF QI TO QS ***1**
!* 3 * PSACI : ACCRETION OF QI TO QS ***3**
!* 4 * PSACW : ACCRETION OF QC BY QS (RIMING) (QSACW FOR PSMLT) ***4**
!* 5 * PRACI : ACCRETION OF QI BY QR ***5**
!* 6 * PIACR : ACCRETION OF QR OR QG BY QI ***6**
!JJS DO 125 J=3,JLES
!JJS DO 125 I=3,ILES
psaut(i,j)=0.0
psaci(i,j)=0.0
praci(i,j)=0.0
piacr(i,j)=0.0
psacw(i,j)=0.0
qsacw(i,j)=0.0
dd(i,j)=1./zs(i,j)**bs3
if (tair(i,j).lt.t0) then
esi(i,j)=exp(.025*tairc(i,j))
psaut(i,j)=r2is*max(rn1*esi(i,j)*(qi(i,j)-bnd1) ,0.0)
psaci(i,j)=r2is*r3f*esi(i,j)*qi(i,j)*dd(i,j)
!JJS 3/30/06
! to cut water to snow accretion by half
! PSACW(I,J)=R4F*QC(I,J)*DD(I,J)
psacw(i,j)=r2is*0.5*r4f*qc(i,j)*dd(i,j)
!JJS 3/30/06
praci(i,j)=r2is*r5f*qi(i,j)/zr(i,j)**bw3
piacr(i,j)=r2is*r6f*qi(i,j)*(zr(i,j)**(-bw6))
!JJS PIACR(I,J)=R6F*QI(I,J)/ZR(I,J)**BW6
else
qsacw(i,j)=r2is*r4f*qc(i,j)*dd(i,j)
endif
!* 21 * PRAUT AUTOCONVERSION OF QC TO QR **21**
!* 22 * PRACW : ACCRETION OF QC BY QR **22**
pracw(i,j)=r22f*qc(i,j)/zr(i,j)**bw3
praut(i,j)=0.0
y1(i,j)=qc(i,j)-bnd3
if (y1(i,j).gt.0.0) then
praut(i,j)=r00*y1(i,j)*y1(i,j)/(1.2e-4+rn21/y1(i,j))
endif
!* 12 * PSFW : BERGERON PROCESSES FOR QS (KOENING, 1971) **12**
!* 13 * PSFI : BERGERON PROCESSES FOR QS **13**
psfw(i,j)=0.0
psfi(i,j)=0.0
pidep(i,j)=0.0
if(tair(i,j).lt.t0.and.qi(i,j).gt.cmin) then
y1(i,j)=max( min(tairc(i,j), -1.), -31.)
it(i,j)=int(abs(y1(i,j)))
y1(i,j)=rn12a(it(i,j))
y2(i,j)=rn12b(it(i,j))
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
psfw(i,j)=r2is* &
max(d2t*y1(i,j)*(y2(i,j)+r12r*qc(i,j))*qi(i,j),0.0)
rtair(i,j)=1./(tair(i,j)-c76)
y2(i,j)=exp(c218-c580*rtair(i,j))
qsi(i,j)=rp0*y2(i,j)
esi(i,j)=c610*y2(i,j)
ssi(i,j)=(qv(i,j)+qb0)/qsi(i,j)-1.
r_nci=min(1.e-6*exp(-.46*tairc(i,j)),1.)
! R_NCI=min(1.e-8*EXP(-.6*TAIRC(I,J)),1.) ! use Tao's
dm(i,j)=max( (qv(i,j)+qb0-qsi(i,j)), 0.)
rsub1(i,j)=cs580*qsi(i,j)*rtair(i,j)*rtair(i,j)
y3(i,j)=1./tair(i,j)
dd(i,j)=y3(i,j)*(rn30a*y3(i,j)-rn30b)+rn30c*tair(i,j)/esi(i,j)
y1(i,j)=206.18*ssi(i,j)/dd(i,j)
pidep(i,j)=y1(i,j)*sqrt(r_nci*qi(i,j)/r00)
dep(i,j)=dm(i,j)/(1.+rsub1(i,j))/d2t
if(dm(i,j).gt.cmin2) then
a2=1.
if(pidep(i,j).gt.dep(i,j).and.pidep(i,j).gt.cmin2) then
a2=dep(i,j)/pidep(i,j)
pidep(i,j)=dep(i,j)
endif
psfi(i,j)=r2is*a2*.5*qi(i,j)*y1(i,j)/(sqrt(ami100) &
-sqrt(ami40))
elseif(dm(i,j).lt.-cmin2) then
!
! SUBLIMATION TERMS USED ONLY WHEN SATURATION ADJUSTMENT FOR ICE
! IS TURNED OFF
!
pidep(i,j)=0.
psfi(i,j)=0.
else
pidep(i,j)=0.
psfi(i,j)=0.
endif
endif
!TTT***** QG=QG+MIN(PGDRY,PGWET)
!* 9 * PGACS : ACCRETION OF QS BY QG (DGACS,WGACS: DRY AND WET) ***9**
!* 14 * DGACW : ACCRETION OF QC BY QG (QGACW FOR PGMLT) **14**
!* 16 * DGACR : ACCRETION OF QR TO QG (QGACR FOR PGMLT) **16**
if(qc(i,j)+qr(i,j).lt.1.e-4) then
ee1=.01
else
ee1=1.
endif
ee2=0.09
egs(i,j)=ee1*exp(ee2*tairc(i,j))
! EGS(I,J)=0.1 ! 6/15/02 tao's
if (tair(i,j).ge.t0) egs(i,j)=1.0
y1(i,j)=abs(vg(i,j)-vs(i,j))
y2(i,j)=zs(i,j)*zg(i,j)
y3(i,j)=5./y2(i,j)
y4(i,j)=.08*y3(i,j)*y3(i,j)
y5(i,j)=.05*y3(i,j)*y4(i,j)
dd(i,j)=y1(i,j)*(y3(i,j)/zs(i,j)**5+y4(i,j)/zs(i,j)**3 &
+y5(i,j)/zs(i,j))
pgacs(i,j)=r2ig*r2is*r9r*egs(i,j)*dd(i,j)
!JJS 1/3/06 from Steve and Chunglin
if (ihail.eq.1) then
dgacs(i,j)=pgacs(i,j)
else
dgacs(i,j)=0.
endif
!JJS 1/3/06 from Steve and Chunglin
wgacs(i,j)=r2ig*r2is*r9r*dd(i,j)
! WGACS(I,J)=0. ! 6/15/02 tao's
y1(i,j)=1./zg(i,j)**bg3
if(ihail .eq. 1) then
dgacw(i,j)=r2ig*max(r14r*qc(i,j)*y1(i,j), 0.0)
else
dgacw(i,j)=r2ig*max(r14f*qc(i,j)*y1(i,j), 0.0)
endif
qgacw(i,j)=dgacw(i,j)
y1(i,j)=abs(vg(i,j)-vr(i,j))
y2(i,j)=zr(i,j)*zg(i,j)
y3(i,j)=5./y2(i,j)
y4(i,j)=.08*y3(i,j)*y3(i,j)
y5(i,j)=.05*y3(i,j)*y4(i,j)
dd(i,j)=r16r*y1(i,j)*(y3(i,j)/zr(i,j)**5+y4(i,j)/zr(i,j)**3 &
+y5(i,j)/zr(i,j))
dgacr(i,j)=r2ig*max(dd(i,j), 0.0)
qgacr(i,j)=dgacr(i,j)
if (tair(i,j).ge.t0) then
dgacs(i,j)=0.0
wgacs(i,j)=0.0
dgacw(i,j)=0.0
dgacr(i,j)=0.0
else
pgacs(i,j)=0.0
qgacw(i,j)=0.0
qgacr(i,j)=0.0
endif
!*******PGDRY : DGACW+DGACI+DGACR+DGACS ******
!* 15 * DGACI : ACCRETION OF QI BY QG (WGACI FOR WET GROWTH) **15**
!* 17 * PGWET : WET GROWTH OF QG **17**
dgaci(i,j)=0.0
wgaci(i,j)=0.0
pgwet(i,j)=0.0
if (tair(i,j).lt.t0) then
y1(i,j)=qi(i,j)/zg(i,j)**bg3
if (ihail.eq.1) then
dgaci(i,j)=r2ig*r15r*y1(i,j)
wgaci(i,j)=r2ig*r15ar*y1(i,j)
! WGACI(I,J)=0. ! 6/15/02 tao's
else
!JJS DGACI(I,J)=r2ig*R15F*Y1(I,J)
dgaci(i,j)=0.
wgaci(i,j)=r2ig*r15af*y1(i,j)
! WGACI(I,J)=0. ! 6/15/02 tao's
endif
!
if (tairc(i,j).ge.-50.) then
if (alf+rn17c*tairc(i,j) .eq. 0.) then
write(91,*) itimestep, i,j,k, alf, rn17c, tairc(i,j)
endif
y1(i,j)=1./(alf+rn17c*tairc(i,j))
if (ihail.eq.1) then
y3(i,j)=.78/zg(i,j)**2+r17aq*scv(i,j)/zg(i,j)**bgh5
else
y3(i,j)=.78/zg(i,j)**2+r17as*scv(i,j)/zg(i,j)**bgh5
endif
y4(i,j)=alvr*dwv(i,j)*(rp0-(qv(i,j)+qb0)) &
-tca(i,j)*tairc(i,j)
dd(i,j)=y1(i,j)*(r17r*y4(i,j)*y3(i,j) &
+(wgaci(i,j)+wgacs(i,j))*(alf+rn17b*tairc(i,j)))
pgwet(i,j)=r2ig*max(dd(i,j), 0.0)
endif
endif
!JJS 125 CONTINUE
!******** HANDLING THE NEGATIVE CLOUD WATER (QC) ******************
!******** HANDLING THE NEGATIVE CLOUD ICE (QI) ******************
!JJS DO 150 J=3,JLES
!JJS DO 150 I=3,ILES
y1(i,j)=qc(i,j)/d2t
psacw(i,j)=min(y1(i,j), psacw(i,j))
praut(i,j)=min(y1(i,j), praut(i,j))
pracw(i,j)=min(y1(i,j), pracw(i,j))
psfw(i,j)= min(y1(i,j), psfw(i,j))
dgacw(i,j)=min(y1(i,j), dgacw(i,j))
qsacw(i,j)=min(y1(i,j), qsacw(i,j))
qgacw(i,j)=min(y1(i,j), qgacw(i,j))
y1(i,j)=(psacw(i,j)+praut(i,j)+pracw(i,j)+psfw(i,j) &
+dgacw(i,j)+qsacw(i,j)+qgacw(i,j))*d2t
qc(i,j)=qc(i,j)-y1(i,j)
if (qc(i,j) .lt. 0.0) then
a1=1.
if (y1(i,j) .ne. 0.0) a1=qc(i,j)/y1(i,j)+1.
psacw(i,j)=psacw(i,j)*a1
praut(i,j)=praut(i,j)*a1
pracw(i,j)=pracw(i,j)*a1
psfw(i,j)=psfw(i,j)*a1
dgacw(i,j)=dgacw(i,j)*a1
qsacw(i,j)=qsacw(i,j)*a1
qgacw(i,j)=qgacw(i,j)*a1
qc(i,j)=0.0
endif
!c
!
!******** SHED PROCESS (WGACR=PGWET-DGACW-WGACI-WGACS)
!c
wgacr(i,j)=pgwet(i,j)-dgacw(i,j)-wgaci(i,j)-wgacs(i,j)
y2(i,j)=dgacw(i,j)+dgaci(i,j)+dgacr(i,j)+dgacs(i,j)
if (pgwet(i,j).ge.y2(i,j)) then
wgacr(i,j)=0.0
wgaci(i,j)=0.0
wgacs(i,j)=0.0
else
dgacr(i,j)=0.0
dgaci(i,j)=0.0
dgacs(i,j)=0.0
endif
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!c
y1(i,j)=qi(i,j)/d2t
psaut(i,j)=min(y1(i,j), psaut(i,j))
psaci(i,j)=min(y1(i,j), psaci(i,j))
praci(i,j)=min(y1(i,j), praci(i,j))
psfi(i,j)= min(y1(i,j), psfi(i,j))
dgaci(i,j)=min(y1(i,j), dgaci(i,j))
wgaci(i,j)=min(y1(i,j), wgaci(i,j))
!
y2(i,j)=(psaut(i,j)+psaci(i,j)+praci(i,j)+psfi(i,j) &
+dgaci(i,j)+wgaci(i,j))*d2t
qi(i,j)=qi(i,j)-y2(i,j)+pidep(i,j)*d2t
if (qi(i,j).lt.0.0) then
a2=1.
if (y2(i,j) .ne. 0.0) a2=qi(i,j)/y2(i,j)+1.
psaut(i,j)=psaut(i,j)*a2
psaci(i,j)=psaci(i,j)*a2
praci(i,j)=praci(i,j)*a2
psfi(i,j)=psfi(i,j)*a2
dgaci(i,j)=dgaci(i,j)*a2
wgaci(i,j)=wgaci(i,j)*a2
qi(i,j)=0.0
endif
!
dlt3(i,j)=0.0
dlt2(i,j)=0.0
!
! DLT4(I,J)=1.0
! if(qc(i,j) .gt. 5.e-4) dlt4(i,j)=0.0
! if(qs(i,j) .le. 1.e-4) dlt4(i,j)=1.0
!
! IF (TAIR(I,J).ge.T0) THEN
! DLT4(I,J)=0.0
! ENDIF
if (tair(i,j).lt.t0) then
if (qr(i,j).lt.1.e-4) then
dlt3(i,j)=1.0
dlt2(i,j)=1.0
endif
if (qs(i,j).ge.1.e-4) then
dlt2(i,j)=0.0
endif
endif
if (ice2 .eq. 1) then
dlt3(i,j)=1.0
dlt2(i,j)=1.0
endif
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
pr(i,j)=(qsacw(i,j)+praut(i,j)+pracw(i,j)+qgacw(i,j))*d2t
ps(i,j)=(psaut(i,j)+psaci(i,j)+psacw(i,j)+psfw(i,j) &
+psfi(i,j)+dlt3(i,j)*praci(i,j))*d2t
! PS(I,J)=(PSAUT(I,J)+PSACI(I,J)+dlt4(i,j)*PSACW(I,J)
! 1 +PSFW(I,J)+PSFI(I,J)+DLT3(I,J)*PRACI(I,J))*D2T
pg(i,j)=((1.-dlt3(i,j))*praci(i,j)+dgaci(i,j)+wgaci(i,j) &
+dgacw(i,j))*d2t
! PG(I,J)=((1.-DLT3(I,J))*PRACI(I,J)+DGACI(I,J)+WGACI(I,J)
! 1 +DGACW(I,J)+(1.-dlt4(i,j))*PSACW(I,J))*D2T
!JJS 150 CONTINUE
!* 7 * PRACS : ACCRETION OF QS BY QR ***7**
!* 8 * PSACR : ACCRETION OF QR BY QS (QSACR FOR PSMLT) ***8**
!JJS DO 175 J=3,JLES
!JJS DO 175 I=3,ILES
y1(i,j)=abs(vr(i,j)-vs(i,j))
y2(i,j)=zr(i,j)*zs(i,j)
y3(i,j)=5./y2(i,j)
y4(i,j)=.08*y3(i,j)*y3(i,j)
y5(i,j)=.05*y3(i,j)*y4(i,j)
pracs(i,j)=r2ig*r2is*r7r*y1(i,j)*(y3(i,j)/zs(i,j)**5 &
+y4(i,j)/zs(i,j)**3+y5(i,j)/zs(i,j))
psacr(i,j)=r2is*r8r*y1(i,j)*(y3(i,j)/zr(i,j)**5 &
+y4(i,j)/zr(i,j)**3+y5(i,j)/zr(i,j))
qsacr(i,j)=psacr(i,j)
if (tair(i,j).ge.t0) then
pracs(i,j)=0.0
psacr(i,j)=0.0
else
qsacr(i,j)=0.0
endif
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!* 2 * PGAUT : AUTOCONVERSION OF QS TO QG ***2**
!* 18 * PGFR : FREEZING OF QR TO QG **18**
pgaut(i,j)=0.0
pgfr(i,j)=0.0
if (tair(i,j) .lt. t0) then
! Y1(I,J)=EXP(.09*TAIRC(I,J))
! PGAUT(I,J)=r2is*max(RN2*Y1(I,J)*(QS(I,J)-BND2), 0.0)
! IF(IHAIL.EQ.1) PGAUT(I,J)=max(RN2*Y1(I,J)*(QS(I,J)-BND2),0.0)
y2(i,j)=exp(rn18a*(t0-tair(i,j)))
!JJS PGFR(I,J)=r2ig*max(R18R*(Y2(I,J)-1.)/ZR(I,J)**7., 0.0)
! pgfr(i,j)=r2ice*max(r18r*(y2(i,j)-1.)* &
! (zr(i,j)**(-7.)), 0.0)
! modify to prevent underflow on some computers (JD)
temp = 1./zr(i,j)
temp = temp*temp*temp*temp*temp*temp*temp
pgfr(i,j)=r2ig*max(r18r*(y2(i,j)-1.)* &
temp, 0.0)
endif
!JJS 175 CONTINUE
!******** HANDLING THE NEGATIVE RAIN WATER (QR) *******************
!******** HANDLING THE NEGATIVE SNOW (QS) *******************
!JJS DO 200 J=3,JLES
!JJS DO 200 I=3,ILES
y1(i,j)=qr(i,j)/d2t
y2(i,j)=-qg(i,j)/d2t
piacr(i,j)=min(y1(i,j), piacr(i,j))
dgacr(i,j)=min(y1(i,j), dgacr(i,j))
wgacr(i,j)=min(y1(i,j), wgacr(i,j))
wgacr(i,j)=max(y2(i,j), wgacr(i,j))
psacr(i,j)=min(y1(i,j), psacr(i,j))
pgfr(i,j)= min(y1(i,j), pgfr(i,j))
del=0.
if(wgacr(i,j) .lt. 0.) del=1.
y1(i,j)=(piacr(i,j)+dgacr(i,j)+(1.-del)*wgacr(i,j) &
+psacr(i,j)+pgfr(i,j))*d2t
qr(i,j)=qr(i,j)+pr(i,j)-y1(i,j)-del*wgacr(i,j)*d2t
if (qr(i,j) .lt. 0.0) then
a1=1.
if(y1(i,j) .ne. 0.) a1=qr(i,j)/y1(i,j)+1.
piacr(i,j)=piacr(i,j)*a1
dgacr(i,j)=dgacr(i,j)*a1
if (wgacr(i,j).gt.0.) wgacr(i,j)=wgacr(i,j)*a1
pgfr(i,j)=pgfr(i,j)*a1
psacr(i,j)=psacr(i,j)*a1
qr(i,j)=0.0
endif
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
prn(i,j)=d2t*((1.-dlt3(i,j))*piacr(i,j)+dgacr(i,j) &
+wgacr(i,j)+(1.-dlt2(i,j))*psacr(i,j)+pgfr(i,j))
ps(i,j)=ps(i,j)+d2t*(dlt3(i,j)*piacr(i,j) &
+dlt2(i,j)*psacr(i,j))
pracs(i,j)=(1.-dlt2(i,j))*pracs(i,j)
y1(i,j)=qs(i,j)/d2t
pgacs(i,j)=min(y1(i,j), pgacs(i,j))
dgacs(i,j)=min(y1(i,j), dgacs(i,j))
wgacs(i,j)=min(y1(i,j), wgacs(i,j))
pgaut(i,j)=min(y1(i,j), pgaut(i,j))
pracs(i,j)=min(y1(i,j), pracs(i,j))
psn(i,j)=d2t*(pgacs(i,j)+dgacs(i,j)+wgacs(i,j) &
+pgaut(i,j)+pracs(i,j))
qs(i,j)=qs(i,j)+ps(i,j)-psn(i,j)
if (qs(i,j).lt.0.0) then
a2=1.
if (psn(i,j) .ne. 0.0) a2=qs(i,j)/psn(i,j)+1.
pgacs(i,j)=pgacs(i,j)*a2
dgacs(i,j)=dgacs(i,j)*a2
wgacs(i,j)=wgacs(i,j)*a2
pgaut(i,j)=pgaut(i,j)*a2
pracs(i,j)=pracs(i,j)*a2
psn(i,j)=psn(i,j)*a2
qs(i,j)=0.0
endif
!
!C PSN(I,J)=D2T*(PGACS(I,J)+DGACS(I,J)+WGACS(I,J)
!c +PGAUT(I,J)+PRACS(I,J))
y2(i,j)=d2t*(psacw(i,j)+psfw(i,j)+dgacw(i,j)+piacr(i,j) &
+dgacr(i,j)+wgacr(i,j)+psacr(i,j)+pgfr(i,j))
pt(i,j)=pt(i,j)+afcp*y2(i,j)
qg(i,j)=qg(i,j)+pg(i,j)+prn(i,j)+psn(i,j)
!JJS 200 CONTINUE
!* 11 * PSMLT : MELTING OF QS **11**
!* 19 * PGMLT : MELTING OF QG TO QR **19**
!JJS DO 225 J=3,JLES
!JJS DO 225 I=3,ILES
psmlt(i,j)=0.0
pgmlt(i,j)=0.0
tair(i,j)=(pt(i,j)+tb0)*pi0
if (tair(i,j).ge.t0) then
tairc(i,j)=tair(i,j)-t0
y1(i,j)=tca(i,j)*tairc(i,j)-alvr*dwv(i,j) &
*(rp0-(qv(i,j)+qb0))
y2(i,j)=.78/zs(i,j)**2+r101f*scv(i,j)/zs(i,j)**bsh5
dd(i,j)=r11rt*y1(i,j)*y2(i,j)+r11at*tairc(i,j) &
*(qsacw(i,j)+qsacr(i,j))
psmlt(i,j)=r2is*max(0.0, min(dd(i,j), qs(i,j)))
if(ihail.eq.1) then
y3(i,j)=.78/zg(i,j)**2+r19aq*scv(i,j)/zg(i,j)**bgh5
else
y3(i,j)=.78/zg(i,j)**2+r19as*scv(i,j)/zg(i,j)**bgh5
endif
dd1(i,j)=r19rt*y1(i,j)*y3(i,j)+r19bt*tairc(i,j) &
*(qgacw(i,j)+qgacr(i,j))
pgmlt(i,j)=r2ig*max(0.0, min(dd1(i,j), qg(i,j)))
pt(i,j)=pt(i,j)-afcp*(psmlt(i,j)+pgmlt(i,j))
qr(i,j)=qr(i,j)+psmlt(i,j)+pgmlt(i,j)
qs(i,j)=qs(i,j)-psmlt(i,j)
qg(i,j)=qg(i,j)-pgmlt(i,j)
endif
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!* 24 * PIHOM : HOMOGENEOUS FREEZING OF QC TO QI (T < T00) **24**
!* 25 * PIDW : DEPOSITION GROWTH OF QC TO QI ( T0 < T <= T00) **25**
!* 26 * PIMLT : MELTING OF QI TO QC (T >= T0) **26**
if (qc(i,j).le.cmin1) qc(i,j)=0.0
if (qi(i,j).le.cmin1) qi(i,j)=0.0
tair(i,j)=(pt(i,j)+tb0)*pi0
if(tair(i,j).le.t00) then
pihom(i,j)=qc(i,j)
else
pihom(i,j)=0.0
endif
if(tair(i,j).ge.t0) then
pimlt(i,j)=qi(i,j)
else
pimlt(i,j)=0.0
endif
pidw(i,j)=0.0
if (tair(i,j).lt.t0 .and. tair(i,j).gt.t00) then
tairc(i,j)=tair(i,j)-t0
y1(i,j)=max( min(tairc(i,j), -1.), -31.)
it(i,j)=int(abs(y1(i,j)))
y2(i,j)=aa1(it(i,j))
y3(i,j)=aa2(it(i,j))
y4(i,j)=exp(abs(beta*tairc(i,j)))
y5(i,j)=(r00*qi(i,j)/(r25a*y4(i,j)))**y3(i,j)
pidw(i,j)=min(r25rt*y2(i,j)*y4(i,j)*y5(i,j), qc(i,j))
endif
y1(i,j)=pihom(i,j)-pimlt(i,j)+pidw(i,j)
pt(i,j)=pt(i,j)+afcp*y1(i,j)+ascp*(pidep(i,j))*d2t
qv(i,j)=qv(i,j)-(pidep(i,j))*d2t
qc(i,j)=qc(i,j)-y1(i,j)
qi(i,j)=qi(i,j)+y1(i,j)
!* 31 * PINT : INITIATION OF QI **31**
!* 32 * PIDEP : DEPOSITION OF QI **32**
!
! CALCULATION OF PINT USES DIFFERENT VALUES OF THE INTERCEPT AND SLOPE FOR
! THE FLETCHER EQUATION. ALSO, ONLY INITIATE MORE ICE IF THE NEW NUMBER
! CONCENTRATION EXCEEDS THAT ALREADY PRESENT.
!* 31 * pint : initiation of qi **31**
!* 32 * pidep : deposition of qi **32**
pint(i,j)=0.0
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
if ( itaobraun.eq.1 ) then
tair(i,j)=(pt(i,j)+tb0)*pi0
if (tair(i,j) .lt. t0) then
! if (qi(i,j) .le. cmin) qi(i,j)=0.
if (qi(i,j) .le. cmin2) qi(i,j)=0.
tairc(i,j)=tair(i,j)-t0
rtair(i,j)=1./(tair(i,j)-c76)
y2(i,j)=exp(c218-c580*rtair(i,j))
qsi(i,j)=rp0*y2(i,j)
esi(i,j)=c610*y2(i,j)
ssi(i,j)=(qv(i,j)+qb0)/qsi(i,j)-1.
ami50=3.76e-8
!ccif ( itaobraun.eq.1 ) --> betah=0.5*beta=-.46*0.5=-0.23; cn0=1.e-6
!ccif ( itaobraun.eq.0 ) --> betah=0.5*beta=-.6*0.5=-0.30; cn0=1.e-8
y1(i,j)=1./tair(i,j)
!cc insert a restriction on ice collection that ice collection
!cc should be stopped at -30 c (with cn0=1.e-6, beta=-.46)
tairccri=tairc(i,j) ! in degree c
if(tairccri.le.-30.) tairccri=-30.
! y2(i,j)=exp(betah*tairc(i,j))
y2(i,j)=exp(betah*tairccri)
y3(i,j)=sqrt(qi(i,j))
dd(i,j)=y1(i,j)*(rn10a*y1(i,j)-rn10b)+rn10c*tair(i,j) &
/esi(i,j)
pidep(i,j)=max(r32rt*ssi(i,j)*y2(i,j)*y3(i,j)/dd(i,j), 0.e0)
r_nci=min(cn0*exp(beta*tairc(i,j)),1.)
! r_nci=min(1.e-6*exp(-.46*tairc(i,j)),1.)
dd(i,j)=max(1.e-9*r_nci/r00-qi(i,j)*1.e-9/ami50, 0.)
dm(i,j)=max( (qv(i,j)+qb0-qsi(i,j)), 0.0)
rsub1(i,j)=cs580*qsi(i,j)*rtair(i,j)*rtair(i,j)
dep(i,j)=dm(i,j)/(1.+rsub1(i,j))
pint(i,j)=max(min(dd(i,j), dm(i,j)), 0.)
! pint(i,j)=min(pint(i,j), dep(i,j))
pint(i,j)=min(pint(i,j)+pidep(i,j), dep(i,j))
! if (pint(i,j) .le. cmin) pint(i,j)=0.
if (pint(i,j) .le. cmin2) pint(i,j)=0.
pt(i,j)=pt(i,j)+ascp*pint(i,j)
qv(i,j)=qv(i,j)-pint(i,j)
qi(i,j)=qi(i,j)+pint(i,j)
endif
endif ! if ( itaobraun.eq.1 )
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
if ( itaobraun.eq.0 ) then
tair(i,j)=(pt(i,j)+tb0)*pi0
if (tair(i,j) .lt. t0) then
if (qi(i,j) .le. cmin1) qi(i,j)=0.
tairc(i,j)=tair(i,j)-t0
dd(i,j)=r31r*exp(beta*tairc(i,j))
rtair(i,j)=1./(tair(i,j)-c76)
y2(i,j)=exp(c218-c580*rtair(i,j))
qsi(i,j)=rp0*y2(i,j)
esi(i,j)=c610*y2(i,j)
ssi(i,j)=(qv(i,j)+qb0)/qsi(i,j)-1.
dm(i,j)=max( (qv(i,j)+qb0-qsi(i,j)), 0.)
rsub1(i,j)=cs580*qsi(i,j)*rtair(i,j)*rtair(i,j)
dep(i,j)=dm(i,j)/(1.+rsub1(i,j))
pint(i,j)=max(min(dd(i,j), dm(i,j)), 0.)
y1(i,j)=1./tair(i,j)
y2(i,j)=exp(betah*tairc(i,j))
y3(i,j)=sqrt(qi(i,j))
dd(i,j)=y1(i,j)*(rn10a*y1(i,j)-rn10b) &
+rn10c*tair(i,j)/esi(i,j)
pidep(i,j)=max(r32rt*ssi(i,j)*y2(i,j)*y3(i,j)/dd(i,j), 0.)
pint(i,j)=pint(i,j)+pidep(i,j)
pint(i,j)=min(pint(i,j),dep(i,j))
!c if (pint(i,j) .le. cmin2) pint(i,j)=0.
pt(i,j)=pt(i,j)+ascp*pint(i,j)
qv(i,j)=qv(i,j)-pint(i,j)
qi(i,j)=qi(i,j)+pint(i,j)
endif
endif ! if ( itaobraun.eq.0 )
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!JJS 225 CONTINUE
!***** TAO ET AL (1989) SATURATION TECHNIQUE ***********************
if (new_ice_sat .eq. 0) then
!JJS DO 250 J=3,JLES
!JJS DO 250 I=3,ILES
tair(i,j)=(pt(i,j)+tb0)*pi0
cnd(i,j)=rt0*(tair(i,j)-t00)
dep(i,j)=rt0*(t0-tair(i,j))
y1(i,j)=1./(tair(i,j)-c358)
y2(i,j)=1./(tair(i,j)-c76)
qsw(i,j)=rp0*exp(c172-c409*y1(i,j))
qsi(i,j)=rp0*exp(c218-c580*y2(i,j))
dd(i,j)=cp409*y1(i,j)*y1(i,j)
dd1(i,j)=cp580*y2(i,j)*y2(i,j)
if (qc(i,j).le.cmin) qc(i,j)=cmin
if (qi(i,j).le.cmin) qi(i,j)=cmin
if (tair(i,j).ge.t0) then
dep(i,j)=0.0
cnd(i,j)=1.
qi(i,j)=0.0
endif
if (tair(i,j).lt.t00) then
cnd(i,j)=0.0
dep(i,j)=1.
qc(i,j)=0.0
endif
y5(i,j)=avcp*cnd(i,j)+ascp*dep(i,j)
! if (qc(i,j) .ge. cmin .or. qi(i,j) .ge. cmin) then
y1(i,j)=qc(i,j)*qsw(i,j)/(qc(i,j)+qi(i,j))
y2(i,j)=qi(i,j)*qsi(i,j)/(qc(i,j)+qi(i,j))
y4(i,j)=dd(i,j)*y1(i,j)+dd1(i,j)*y2(i,j)
qvs(i,j)=y1(i,j)+y2(i,j)
rsub1(i,j)=(qv(i,j)+qb0-qvs(i,j))/(1.+y4(i,j)*y5(i,j))
cnd(i,j)=cnd(i,j)*rsub1(i,j)
dep(i,j)=dep(i,j)*rsub1(i,j)
if (qc(i,j).le.cmin) qc(i,j)=0.
if (qi(i,j).le.cmin) qi(i,j)=0.
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!c ****** condensation or evaporation of qc ******
cnd(i,j)=max(-qc(i,j),cnd(i,j))
!c ****** deposition or sublimation of qi ******
dep(i,j)=max(-qi(i,j),dep(i,j))
pt(i,j)=pt(i,j)+avcp*cnd(i,j)+ascp*dep(i,j)
qv(i,j)=qv(i,j)-cnd(i,j)-dep(i,j)
qc(i,j)=qc(i,j)+cnd(i,j)
qi(i,j)=qi(i,j)+dep(i,j)
!JJS 250 continue
endif
if (new_ice_sat .eq. 1) then
!JJS DO J=3,JLES
!JJS DO I=3,ILES
tair(i,j)=(pt(i,j)+tb0)*pi0
cnd(i,j)=rt0*(tair(i,j)-t00)
dep(i,j)=rt0*(t0-tair(i,j))
y1(i,j)=1./(tair(i,j)-c358)
y2(i,j)=1./(tair(i,j)-c76)
qsw(i,j)=rp0*exp(c172-c409*y1(i,j))
qsi(i,j)=rp0*exp(c218-c580*y2(i,j))
dd(i,j)=cp409*y1(i,j)*y1(i,j)
dd1(i,j)=cp580*y2(i,j)*y2(i,j)
y5(i,j)=avcp*cnd(i,j)+ascp*dep(i,j)
y1(i,j)=rt0*(tair(i,j)-t00)*qsw(i,j)
y2(i,j)=rt0*(t0-tair(i,j))*qsi(i,j)
! IF (QC(I,J).LE.CMIN) QC(I,J)=CMIN
! IF (QI(I,J).LE.CMIN) QI(I,J)=CMIN
if (tair(i,j).ge.t0) then
! QI(I,J)=0.0
dep(i,j)=0.0
cnd(i,j)=1.
y2(i,j)=0.
y1(i,j)=qsw(i,j)
endif
if (tair(i,j).lt.t00) then
cnd(i,j)=0.0
dep(i,j)=1.
y2(i,j)=qsi(i,j)
y1(i,j)=0.
! QC(I,J)=0.0
endif
! Y1(I,J)=QC(I,J)*QSW(I,J)/(QC(I,J)+QI(I,J))
! Y2(I,J)=QI(I,J)*QSI(I,J)/(QC(I,J)+QI(I,J))
y4(i,j)=dd(i,j)*y1(i,j)+dd1(i,j)*y2(i,j)
qvs(i,j)=y1(i,j)+y2(i,j)
rsub1(i,j)=(qv(i,j)+qb0-qvs(i,j))/(1.+y4(i,j)*y5(i,j))
cnd(i,j)=cnd(i,j)*rsub1(i,j)
dep(i,j)=dep(i,j)*rsub1(i,j)
! IF (QC(I,J).LE.CMIN) QC(I,J)=0.
! IF (QI(I,J).LE.CMIN) QI(I,J)=0.
!C ****** CONDENSATION OR EVAPORATION OF QC ******
cnd(i,j)=max(-qc(i,j),cnd(i,j))
!C ****** DEPOSITION OR SUBLIMATION OF QI ******
dep(i,j)=max(-qi(i,j),dep(i,j))
pt(i,j)=pt(i,j)+avcp*cnd(i,j)+ascp*dep(i,j)
qv(i,j)=qv(i,j)-cnd(i,j)-dep(i,j)
qc(i,j)=qc(i,j)+cnd(i,j)
qi(i,j)=qi(i,j)+dep(i,j)
!JJS ENDDO
!JJS ENDDO
endif
!c
!
if (new_ice_sat .eq. 2) then
!JJS do j=3,jles
!JJS do i=3,iles
dep(i,j)=0.0
cnd(i,j)=0.0
tair(i,j)=(pt(i,j)+tb0)*pi0
if (tair(i,j) .ge. 253.16) then
y1(i,j)=1./(tair(i,j)-c358)
qsw(i,j)=rp0*exp(c172-c409*y1(i,j))
dd(i,j)=cp409*y1(i,j)*y1(i,j)
dm(i,j)=qv(i,j)+qb0-qsw(i,j)
cnd(i,j)=dm(i,j)/(1.+avcp*dd(i,j)*qsw(i,j))
!c ****** condensation or evaporation of qc ******
cnd(i,j)=max(-qc(i,j), cnd(i,j))
pt(i,j)=pt(i,j)+avcp*cnd(i,j)
qv(i,j)=qv(i,j)-cnd(i,j)
qc(i,j)=qc(i,j)+cnd(i,j)
endif
if (tair(i,j) .le. 258.16) then
!c cnd(i,j)=0.0
y2(i,j)=1./(tair(i,j)-c76)
qsi(i,j)=rp0*exp(c218-c580*y2(i,j))
dd1(i,j)=cp580*y2(i,j)*y2(i,j)
dep(i,j)=(qv(i,j)+qb0-qsi(i,j))/(1.+ascp*dd1(i,j)*qsi(i,j))
!c ****** deposition or sublimation of qi ******
dep(i,j)=max(-qi(i,j),dep(i,j))
pt(i,j)=pt(i,j)+ascp*dep(i,j)
qv(i,j)=qv(i,j)-dep(i,j)
qi(i,j)=qi(i,j)+dep(i,j)
endif
!JJS enddo
!JJS enddo
endif
!c
!
!* 10 * PSDEP : DEPOSITION OR SUBLIMATION OF QS **10**
!* 20 * PGSUB : SUBLIMATION OF QG **20**
!JJS DO 275 J=3,JLES
!JJS DO 275 I=3,ILES
psdep(i,j)=0.0
pgdep(i,j)=0.0
pssub(i,j)=0.0
pgsub(i,j)=0.0
tair(i,j)=(pt(i,j)+tb0)*pi0
if(tair(i,j).lt.t0) then
if(qs(i,j).lt.cmin1) qs(i,j)=0.0
if(qg(i,j).lt.cmin1) qg(i,j)=0.0
rtair(i,j)=1./(tair(i,j)-c76)
qsi(i,j)=rp0*exp(c218-c580*rtair(i,j))
ssi(i,j)=(qv(i,j)+qb0)/qsi(i,j)-1.
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
y1(i,j)=r10ar/(tca(i,j)*tair(i,j)**2)+1./(dwv(i,j) &
*qsi(i,j))
y2(i,j)=.78/zs(i,j)**2+r101f*scv(i,j)/zs(i,j)**bsh5
psdep(i,j)=r10t*ssi(i,j)*y2(i,j)/y1(i,j)
pssub(i,j)=psdep(i,j)
psdep(i,j)=r2is*max(psdep(i,j), 0.)
pssub(i,j)=r2is*max(-qs(i,j), min(pssub(i,j), 0.))
if(ihail.eq.1) then
y2(i,j)=.78/zg(i,j)**2+r20bq*scv(i,j)/zg(i,j)**bgh5
else
y2(i,j)=.78/zg(i,j)**2+r20bs*scv(i,j)/zg(i,j)**bgh5
endif
pgsub(i,j)=r2ig*r20t*ssi(i,j)*y2(i,j)/y1(i,j)
dm(i,j)=qv(i,j)+qb0-qsi(i,j)
rsub1(i,j)=cs580*qsi(i,j)*rtair(i,j)*rtair(i,j)
! ******** DEPOSITION OR SUBLIMATION OF QS **********************
y1(i,j)=dm(i,j)/(1.+rsub1(i,j))
psdep(i,j)=r2is*min(psdep(i,j),max(y1(i,j),0.))
y2(i,j)=min(y1(i,j),0.)
pssub(i,j)=r2is*max(pssub(i,j),y2(i,j))
! ******** SUBLIMATION OF QG ***********************************
dd(i,j)=max((-y2(i,j)-qs(i,j)), 0.)
pgsub(i,j)=r2ig*min(dd(i,j), qg(i,j), max(pgsub(i,j),0.))
if(qc(i,j)+qi(i,j).gt.1.e-5) then
dlt1(i,j)=1.
else
dlt1(i,j)=0.
endif
psdep(i,j)=dlt1(i,j)*psdep(i,j)
pssub(i,j)=(1.-dlt1(i,j))*pssub(i,j)
pgsub(i,j)=(1.-dlt1(i,j))*pgsub(i,j)
pt(i,j)=pt(i,j)+ascp*(psdep(i,j)+pssub(i,j)-pgsub(i,j))
qv(i,j)=qv(i,j)+pgsub(i,j)-pssub(i,j)-psdep(i,j)
qs(i,j)=qs(i,j)+psdep(i,j)+pssub(i,j)
qg(i,j)=qg(i,j)-pgsub(i,j)
endif
!* 23 * ERN : EVAPORATION OF QR (SUBSATURATION) **23**
ern(i,j)=0.0
if(qr(i,j).gt.0.0) then
tair(i,j)=(pt(i,j)+tb0)*pi0
rtair(i,j)=1./(tair(i,j)-c358)
qsw(i,j)=rp0*exp(c172-c409*rtair(i,j))
ssw(i,j)=(qv(i,j)+qb0)/qsw(i,j)-1.0
dm(i,j)=qv(i,j)+qb0-qsw(i,j)
rsub1(i,j)=cv409*qsw(i,j)*rtair(i,j)*rtair(i,j)
dd1(i,j)=max(-dm(i,j)/(1.+rsub1(i,j)), 0.0)
y1(i,j)=.78/zr(i,j)**2+r23af*scv(i,j)/zr(i,j)**bwh5
y2(i,j)=r23br/(tca(i,j)*tair(i,j)**2)+1./(dwv(i,j) &
*qsw(i,j))
!cccc
ern(i,j)=r23t*ssw(i,j)*y1(i,j)/y2(i,j)
ern(i,j)=min(dd1(i,j),qr(i,j),max(ern(i,j),0.))
pt(i,j)=pt(i,j)-avcp*ern(i,j)
qv(i,j)=qv(i,j)+ern(i,j)
qr(i,j)=qr(i,j)-ern(i,j)
endif
!JJS 10/7/2008 vvvvv
ENDIF ! part of if (iwarm.eq.1) then
!JJS 10/7/2008 ^^^^^
! IF (QV(I,J)+QB0 .LE. 0.) QV(I,J)=-QB0
if (qc(i,j) .le. cmin1) qc(i,j)=0.
if (qr(i,j) .le. cmin1) qr(i,j)=0.
if (qi(i,j) .le. cmin1) qi(i,j)=0.
if (qs(i,j) .le. cmin1) qs(i,j)=0.
if (qg(i,j) .le. cmin1) qg(i,j)=0.
dpt(i,j,k)=pt(i,j)
dqv(i,j,k)=qv(i,j)
qcl(i,j,k)=qc(i,j)
qrn(i,j,k)=qr(i,j)
qci(i,j,k)=qi(i,j)
qcs(i,j,k)=qs(i,j)
qcg(i,j,k)=qg(i,j)
!JJS 275 CONTINUE
scc=0.
see=0.
! DO 110 J=3,JLES
! DO 110 I=3,ILES
! DD(I,J)=MAX(-CND(I,J), 0.)
! CND(I,J)=MAX(CND(I,J), 0.)
! DD1(I,J)=MAX(-DEP(I,J), 0.)
!ccshie 2/21/02 shie follow tao
!CC for reference QI(I,J)=QI(I,J)-Y2(I,J)+PIDEP(I,J)*D2T
!CC for reference QV(I,J)=QV(I,J)-(PIDEP(I,J))*D2T
!c DEP(I,J)=MAX(DEP(I,J), 0.)
! DEP(I,J)=MAX(DEP(I,J), 0.)+PIDEP(I,J)*D2T
! SCC=SCC+CND(I,J)
! SEE=SEE+DD(I,J)+ERN(I,J)
! 110 CONTINUE
! SC(K)=SCC+SC(K)
! SE(K)=SEE+SE(K)
!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!c henry: please take a look (start)
!JJS modified by JJS on 5/1/2007 vvvvv
!JJS do 305 j=3,jles
!JJS do 305 i=3,iles
dd(i,j)=max(-cnd(i,j), 0.)
cnd(i,j)=max(cnd(i,j), 0.)
dd1(i,j)=max(-dep(i,j), 0.)+pidep(i,j)*d2t
dep(i,j)=max(dep(i,j), 0.)
!JJS 305 continue
!JJS do 306 j=3,jles
!JJS do 306 i=3,iles
!JJS scc=scc+cnd(i,j)
!JJS see=see+(dd(i,j)+ern(i,j))
!
!JJS sddd=sddd+(dep(i,j)+pint(i,j)+psdep(i,j)+pgdep(i,j))
!JJS ssss=ssss+dd1(i,j)
!JJS
! shhh=shhh+rsw(i,j,k)*d2t
! sccc=sccc+rlw(i,j,k)*d2t
!jjs
!JJS smmm=smmm+(psmlt(i,j)+pgmlt(i,j)+pimlt(i,j))
!JJS sfff=sfff+d2t*(psacw(i,j)+piacr(i,j)+psfw(i,j)
!JJS 1 +pgfr(i,j)+dgacw(i,j)+dgacr(i,j)+psacr(i,j))
!JJS 2 -qracs(i,j)+pihom(i,j)+pidw(i,j)
sccc=cnd(i,j)
seee=dd(i,j)+ern(i,j)
sddd=dep(i,j)+pint(i,j)+psdep(i,j)+pgdep(i,j)
ssss=dd1(i,j) + pgsub(i,j)
smmm=psmlt(i,j)+pgmlt(i,j)+pimlt(i,j)
sfff=d2t*(psacw(i,j)+piacr(i,j)+psfw(i,j) &
+pgfr(i,j)+dgacw(i,j)+dgacr(i,j)+psacr(i,j) &
+wgacr(i,j))+pihom(i,j)+pidw(i,j)
! physc(i,k,j) = avcp * sccc / d2t
! physe(i,k,j) = avcp * seee / d2t
! physd(i,k,j) = ascp * sddd / d2t
! physs(i,k,j) = ascp * ssss / d2t
! physf(i,k,j) = afcp * sfff / d2t
! physm(i,k,j) = afcp * smmm / d2t
! physc(i,k,j) = physc(i,k,j) + avcp * sccc
! physe(i,k,j) = physc(i,k,j) + avcp * seee
! physd(i,k,j) = physd(i,k,j) + ascp * sddd
! physs(i,k,j) = physs(i,k,j) + ascp * ssss
! physf(i,k,j) = physf(i,k,j) + afcp * sfff
! physm(i,k,j) = physm(i,k,j) + afcp * smmm
!JJS modified by JJS on 5/1/2007 ^^^^^
2000 continue
1000 continue
!JJS ****************************************************************
!JJS convert from GCE grid back to WRF grid
do k=kts,kte
do j=jts,jte
do i=its,ite
ptwrf(i,k,j) = dpt(i,j,k)
qvwrf(i,k,j) = dqv(i,j,k)
qlwrf(i,k,j) = qcl(i,j,k)
qrwrf(i,k,j) = qrn(i,j,k)
qiwrf(i,k,j) = qci(i,j,k)
qswrf(i,k,j) = qcs(i,j,k)
qgwrf(i,k,j) = qcg(i,j,k)
enddo !i
enddo !j
enddo !k
! ****************************************************************
!+---+-----------------------------------------------------------------+
IF ( PRESENT (diagflag) ) THEN
if (diagflag .and. do_radar_ref == 1) then
do j=jts,jte
do i=its,ite
DO K=kts,kte
t1d(k)=ptwrf(i,k,j)*pi_mks(i,k,j)
p1d(k)=p0_mks(i,k,j)
qv1d(k)=qvwrf(i,k,j)
qr1d(k)=qrwrf(i,k,j)
ENDDO
if (ice2.eq.0) then
DO K=kts,kte
qs1d(k)=qswrf(i,k,j)
qg1d(k)=qgwrf(i,k,j)
ENDDO
elseif (ice2.eq.1) then
DO K=kts,kte
qs1d(k)=qswrf(i,k,j)
ENDDO
elseif (ice2.eq.2) then
DO K=kts,kte
qs1d(k)=0.
qg1d(k)=qgwrf(i,k,j)
ENDDO
elseif (ice2.eq.3) then
DO K=kts,kte
qs1d(k)=0.
qg1d(k)=0.
ENDDO
endif
call refl10cm_gsfc (qv1d, qr1d, qs1d, qg1d, &
t1d, p1d, dBZ, kts, kte, i, j, ihail)
do k = kts, kte
refl_10cm(i,k,j) = MAX(-35., dBZ(k))
enddo
enddo
enddo
endif
ENDIF
!+---+-----------------------------------------------------------------+
return
END SUBROUTINE saticel_s
!JJS
!JJS REAL FUNCTION GAMMA(X)
!JJS Y=GAMMLN(X)
!JJS GAMMA=EXP(Y)
!JJS RETURN
!JJS END
!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!JJS real function GAMMLN (xx)
real function gammagce (xx) 10
!**********************************************************************
real*8 cof(6),stp,half,one,fpf,x,tmp,ser
data cof,stp / 76.18009173,-86.50532033,24.01409822, &
-1.231739516,.120858003e-2,-.536382e-5, 2.50662827465 /
data half,one,fpf / .5, 1., 5.5 /
!
x=xx-one
tmp=x+fpf
tmp=(x+half)*log(tmp)-tmp
ser=one
do j=1,6
x=x+one
ser=ser+cof(j)/x
enddo !j
gammln=tmp+log(stp*ser)
!JJS
gammagce=exp(gammln)
!JJS
return
END FUNCTION gammagce
!+---+-----------------------------------------------------------------+
subroutine refl10cm_gsfc (qv1d, qr1d, qs1d, qg1d, & 1,2
t1d, p1d, dBZ, kts, kte, ii, jj, ihail)
IMPLICIT NONE
!..Sub arguments
INTEGER, INTENT(IN):: kts, kte, ii, jj, ihail
REAL, DIMENSION(kts:kte), INTENT(IN):: &
qv1d, qr1d, qs1d, qg1d, 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, rg
DOUBLE PRECISION, DIMENSION(kts:kte):: ilamr, ilams, ilamg
DOUBLE PRECISION, DIMENSION(kts:kte):: N0_r, N0_s, N0_g
DOUBLE PRECISION:: lamr, lams, lamg
LOGICAL, DIMENSION(kts:kte):: L_qr, L_qs, L_qg
REAL, DIMENSION(kts:kte):: ze_rain, ze_snow, ze_graupel
DOUBLE PRECISION:: fmelt_s, fmelt_g
INTEGER:: i, k, k_0, kbot, n
LOGICAL:: melti
DOUBLE PRECISION:: cback, x, eta, f_d
REAL, PARAMETER:: R=287.
REAL, PARAMETER:: PIx=3.1415926536
!+---+
do k = kts, kte
dBZ(k) = -35.0
enddo
!+---+-----------------------------------------------------------------+
!..Put column of data into local arrays.
!+---+-----------------------------------------------------------------+
do k = kts, kte
temp(k) = t1d(k)
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) = xnor
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) = xnos
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
if (qg1d(k) .gt. 1.E-9) then
rg(k) = qg1d(k)*rho(k)
if (ihail.eq.1) then
N0_g(k) = xnoh
else
N0_g(k) = xnog
endif
lamg = (xam_g*xcgg(3)*N0_g(k)/rg(k))**(1./xcge(1))
ilamg(k) = 1./lamg
L_qg(k) = .true.
else
rg(k) = 1.E-12
L_qg(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).or.L_qg(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
ze_graupel(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/PIx)*(6.0/PIx) &
* (xam_s/900.0)*(xam_s/900.0) &
* N0_s(k)*xcsg(4)*ilams(k)**xcse(4)
if (L_qg(k)) ze_graupel(k) = (0.176/0.93) * (6.0/PIx)*(6.0/PIx) &
* (xam_g/900.0)*(xam_g/900.0) &
* N0_g(k)*xcgg(4)*ilamg(k)**xcge(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
!..Reflectivity contributed by melting graupel
if (L_qg(k) .and. L_qg(k_0) ) then
fmelt_g = MAX(0.005d0, MIN(1.0d0-rg(k)/rg(k_0), 0.99d0))
eta = 0.d0
lamg = 1./ilamg(k)
do n = 1, nrbins
x = xam_g * xxDg(n)**xbm_g
call rayleigh_soak_wetgraupel (x,DBLE(xocmg),DBLE(xobmg), &
fmelt_g, melt_outside_g, m_w_0, m_i_0, lamda_radar, &
CBACK, mixingrulestring_g, matrixstring_g, &
inclusionstring_g, hoststring_g, &
hostmatrixstring_g, hostinclusionstring_g)
f_d = N0_g(k)*xxDg(n)**xmu_g * DEXP(-lamg*xxDg(n))
eta = eta + f_d * CBACK * simpson(n) * xdtg(n)
enddo
ze_graupel(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)+ze_graupel(k))*1.d18)
enddo
end subroutine refl10cm_gsfc
!+---+-----------------------------------------------------------------+
END MODULE module_mp_gsfcgce