!**********************************************************************c
SUBROUTINE CALVIS_GSD(CZEN,VIS)
! SUBPROGRAM: CALVIS CALCULATE HORIZONTAL VISIBILITY
!
! PRGMMR: BENJAMIN, STAN ORG: NOAA/FSL DATE: 99-09-07
!
! ABSTRACT:
!
! Started with Stoelinga-Warner algorithm for hydrometeors only.
! Added coefficients for graupel.
! Added algorithm for clear-air RH-based visibility.
!
! This routine computes horizontal visibility (in km) at the
! surface or lowest model layer, from qc, qr, qi, qs, and qg.
! qv--water vapor mixing ratio (kg/kg)
! qc--cloud water mixing ratio (kg/kg)
! qr--rain water mixing ratio (kg/kg)
! qi--cloud ice mixing ratio (kg/kg)
! qs--snow mixing ratio (kg/kg)
! qg--graupel mixing ratio (kg/kg)
! u/v - u/v wind components (m/s)
! tb -- temperature (k)
! pp--pressure (Pa)
! rhb-- relative humidity (0-100%)
!
!
! Independent of the above definitions, the scheme can use different
! assumptions of the state of hydrometeors:
! meth='r': Uses the four mixing ratios qrain, qsnow, qclw,
! and qclice
!
! The routine uses the following
! expressions for extinction coefficient, beta (in km**-1),
! with C being the mass concentration (in g/m**3):
!
! cloud water: beta = 144.7 * C ** (0.8800)
! rain water: beta = 2.24 * C ** (0.7500)
! cloud ice: beta = 327.8 * C ** (1.0000)
! snow: beta = 10.36 * C ** (0.7776)
! graupel: beta = 8.0 * C ** (0.7500)
!
! These expressions were obtained from the following sources:
!
! for cloud water: from Kunkel (1984)
! for rainwater: from M-P dist'n, with No=8e6 m**-4 and
! rho_w=1000 kg/m**3
! for cloud ice: assume randomly oriented plates which follow
! mass-diameter relationship from Rutledge and Hobbs (1983)
! for snow: from Stallabrass (1985), assuming beta = -ln(.02)/vis
! for graupel: guestimate by John Brown and Stan Benjamin,
! similar to snow, but a smaller extinction coef seemed
! reasonable. 27 Aug 99
!
! The extinction coefficient for each water species present is
! calculated, and then all applicable betas are summed to yield
! a single beta. Then the following relationship is used to
! determine visibility (in km), where epsilon is the threshhold
! of contrast, usually taken to be .02:
!
! vis = -ln(epsilon)/beta [found in Kunkel (1984)]
!
! HISTORY
! PROGRAM HISTORY LOG:
! 99-05- Version from Eta model and from
! Mark Stoelinga and Tom Warner
! 99-09-07 S. Benjamin Modified for MM5 microphysics variables
! include graupel mixing ratio
! 99-09 S. Benjamin Added algorithm for RH-based clear-air
! visibility
! 00-08 S. Benjamin Added mods for base of 60km instead of 90km,
! max RH from lowest 2 levels instead of
! lev 2 only, max hydrometeor mix ratio
! from lowest 5 levs instead of lev 1 only
! Based on Schwartz stats and Smirnova et al
! paper, and on METAR verif started this week
! Dec 03 S. Benjamin - updates
! - day/night distinction for vis constants
! from Roy Rasmussen
! - low-level wind shear term
! - recommended by Evan Kuchera
!
!------------------------------------------------------------------
!
use vrbls3d, only: qqw, qqi, qqs, qqr, qqg, t, pmid, q, u, v, extcof55
use params_mod, only: h1, d608, rd
use ctlblk_mod, only: jm, im, jsta_2l, jend_2u, lm
implicit none
integer :: j, i, k, ll
real :: tx, pol, esx, es, e
REAL VIS(IM,jsta_2l:jend_2u) ,RHB(IM,jsta_2l:jend_2u,LM), CZEN(IM,jsta_2l:jend_2u)
real celkel,tice,coeflc,coeflp,coeffc,coeffp,coeffg
real exponlc,exponlp,exponfc,exponfp,exponfg,const1
real rhoice,rhowat,qrain,qsnow,qgraupel,qclw,qclice,tv,rhoair, &
vovermd,conclc,conclp,concfc,concfp,concfg,betav
real coeffp_dry, coeffp_wet, shear_fac, temp_fac
real coef_snow, shear
real coefrh,qrh,visrh
real rhmax,shear5_cnt, shear8_cnt
real shear5_cnt_lowvis, shear8_cnt_lowvis
real shear4_cnt, shear4_cnt_lowvis
integer night_cnt, lowsun_cnt
real visrh10_cnt, vis1km_cnt, visrh_lower
real vis3km_cnt
real vis5km_cnt
real vis_min, visrh_min
real vis_night, zen_fac
!------------------------------------------------------------------
CELKEL = 273.15
TICE = CELKEL-10.
COEFLC = 144.7
COEFLP = 2.24
COEFFC = 327.8
COEFFP = 10.36
! - Initialize counters
shear4_cnt = 0
shear5_cnt = 0
shear8_cnt = 0
shear4_cnt_lowvis = 0
shear5_cnt_lowvis = 0
shear8_cnt_lowvis = 0
night_cnt = 0
lowsun_cnt = 0
visrh10_cnt = 0
visrh_lower = 0
vis1km_cnt = 0
vis3km_cnt = 0
vis5km_cnt = 0
! - values from Roy Rasmussen - Dec 2003
! COEFFP_dry = 17.7
! COEFFP_wet = 4.18
! - modified number - Stan B. - Dec 2007
! after quick talks with John Brown and Ismail Gultepe
COEFFP_dry = 10.0
COEFFP_wet = 6.0
! COEFFg = 8.0
! - values from Roy Rasmussen - Dec 2003
COEFFg = 4.0
EXPONLC = 0.8800
EXPONLP = 0.7500
EXPONFC = 1.0000
! EXPONFP = 0.7776
! - new value from Roy Rasmussen - Dec 2003
EXPONFP = 1.0
EXPONFg = 0.75
! CONST1=-LOG(.02)
CONST1= 3.912
! visibility with respect to RH is
! calculated from optical depth linearly
! related to RH as follows:
! vis = 60 exp (-3 * (RH-15)/80)
! changed on 8/23/00
! vis = 60 exp (-2.5 * (RH-15)/80)
! changed on 3/14/01
! Previous algorithm gave vis of 3km, not 5 km
! at 95% RH, now fixed. Stan B.
! coefficient of 3 gives visibility of 5 km
! at 95% RH
! Total visibility is minimum of vis-rh
! and vis-hydrometeors from Stoelinga/Warner
RHOICE=917.
RHOWAT=1000.
vis_min = 1.e6
visrh_min = 1.e6
DO J=jsta_2l,jend_2u
DO I=1,IM
! - take max hydrometeor mixing ratios in lowest 25 mb (lowest 5 levels)
! - change - 3/8/01 - Stan B. - based on apparent underforecasting
! of visibility (vis too low from 20km RUC)
! - take max hydrometeor mixing ratios in lowest 15 mb (lowest 4 levels)
qrain = 0.
qsnow = 0.
qgraupel = 0.
qclw = 0.
qclice = 0.
!???? - Stan B. - questions as of 1/1/04
!???? - use mean of these values over
! lowest 2 levels
! do k = 1,4
! QRAIN = qrain+QR(I,J,k)
! QSNOW = qsnow+max(0.,QS(I,J,k) )
! Qgraupel= qgraupel+Qg(I,J,k)
! end do
do k = 1,3
LL=LM-k+1
QCLW = max(qclw, QQW(I,J,ll) )
QCLICE = max(qclice, QQI(I,J,ll) )
Qsnow = max(qsnow, QQS(I,J,ll) )
Qrain = max(qrain, QQR(I,J,ll) )
Qgraupel = max(qgraupel, QQG(I,J,ll) )
! - compute relative humidity
Tx=T(I,J,LL)-273.15
POL = 0.99999683 + TX*(-0.90826951E-02 + &
TX*(0.78736169E-04 + TX*(-0.61117958E-06 + &
TX*(0.43884187E-08 + TX*(-0.29883885E-10 + &
TX*(0.21874425E-12 + TX*(-0.17892321E-14 + &
TX*(0.11112018E-16 + TX*(-0.30994571E-19)))))))))
esx = 6.1078/POL**8
ES = esx
E = PMID(I,J,LL)/100.*Q(I,J,LL)/(0.62197+Q(I,J,LL)*0.37803)
RHB(I,J,LL) = 100.*AMIN1(1.,E/ES)
enddo
! qrain = max(0., qrain / 4.)
! qsnow = max(0., qsnow / 4.)
! qgraupel = max(0., qgraupel / 4.)
! qclw = max(0., qclw / 2.)
! qclice = max(0., qclice / 2.)
! - take max RH of levels 1 and 2 near the sfc
rhmax = max (rhb(i,j,lm),rhb(i,j,lm-1))
! rhmax = max (rhb(i,j,1),rhb(i,j,2))
! - vary RH coefficient between 75% up to max at 95%
! to give 5.4 km visibility at 95%.
! qrh = max(0.0,min(1.0,(rhb(i,j)/100.-0.15)/0.80))
qrh = max(0.0,min(0.8,(rhmax/100.-0.15)))
!15aug11 visrh = 80. * exp(-2.5*qrh)
visrh = 60. * exp(-2.5*qrh)
! -- add term to increase RH vis term for
! low-level wind shear increasing from 4 to 6 ms-1
! (using Evan Kuchera's paper as a guideline)
! -- calculate term for shear between levels 1 and 4
! (about 15 mb)
shear = sqrt( (u(i,j,lm-3)-u(i,j,lm))**2 &
+(v(i,j,lm-3)-v(i,j,lm))**2 )
! shear = sqrt( (u(i,j,4)-u(i,j,1))**2
! 1 +(v(i,j,4)-v(i,j,1))**2 )
! shear_fac = min(1.,max(0.,(shear-5.)/3.) )
shear_fac = min(1.,max(0.,(shear-4.)/2.) )
if (visrh.lt.10.) visrh = visrh + (10.-visrh)* &
shear_fac
if (shear.gt.4.) shear4_cnt = shear4_cnt +1
if (shear.gt.5.) shear5_cnt = shear5_cnt +1
if (shear.gt.6.) shear8_cnt = shear8_cnt +1
if (shear.gt.4..and.visrh.lt.10) &
shear4_cnt_lowvis = shear4_cnt_lowvis +1
if (shear.gt.5..and.visrh.lt.10) &
shear5_cnt_lowvis = shear5_cnt_lowvis +1
if (shear.gt.6..and.visrh.lt.10) &
shear8_cnt_lowvis = shear8_cnt_lowvis +1
if (visrh.lt.10.) visrh10_cnt = visrh10_cnt+1
!new
if (czen(i,j).lt.0.) night_cnt = night_cnt + 1
if (czen(i,j).lt.0.1) lowsun_cnt = lowsun_cnt + 1
TV=T(I,J,lm)*(H1+D608*Q(I,J,lm))
! tv = t(i,j,1)*(1. + 0.6078*q(i,j,1))
RHOAIR=PMID(I,J,lm)/(RD*TV)
! RHOAIR=PMID(I,J,1)/(RD*TV)
VOVERMD=(1.+Q(I,J,lm))/RHOAIR+(QCLW+QRAIN)/RHOWAT+ &
! VOVERMD=(1.+Q(I,J,1))/RHOAIR+(QCLW+QRAIN)/RHOWAT+
(qgraupel+QCLICE+QSNOW)/RHOICE
CONCLC=QCLW/VOVERMD*1000.
CONCLP=QRAIN/VOVERMD*1000.
CONCFC=QCLICE/VOVERMD*1000.
CONCFP=QSNOW/VOVERMD*1000.
CONCFg=Qgraupel/VOVERMD*1000.
temp_fac = min(1.,max((t(i,j,lm)-271.15),0.) )
! temp_fac = min(1.,max((t(i,j,1)-271.15),0.) )
coef_snow = coeffp_dry*(1.-temp_fac) &
+ coeffp_wet* temp_fac
if (t(i,j,lm).lt. 270. .and. temp_fac.eq.1.) &
! if (t(i,j,1).lt. 270. .and. temp_fac.eq.1.)
write (6,*) 'Problem w/ temp_fac - calvis'
! BETAV=COEFFC*CONCFC**EXPONFC+COEFFP*CONCFP**EXPONFP
BETAV=COEFFC*CONCFC**EXPONFC &
+ coef_SNOW*CONCFP**EXPONFP &
+ COEFLC*CONCLC**EXPONLC+COEFLP*CONCLP**EXPONLP &
+ coeffg*concfg**exponfg +1.E-10
if (i.eq.290 .and. j.eq.112) then
write (6,*) 'BETAV, extcof55 =',BETAV,extcof55(i,j,lm)
end if
VIS(I,J)=MIN(90.,CONST1/BETAV+extcof55(i,j,lm)) ! max of 90km
if (vis(i,j).lt.vis_min) vis_min = vis(i,j)
if (visrh.lt.visrh_min) visrh_min = visrh
if (visrh.lt.vis(i,j)) visrh_lower = visrh_lower + 1
! -- Dec 2003 - Roy Rasmussen (NCAR) expression for night vs. day vis
! 1.609 factor is number of km in mile.
vis_night = 1.69 * ((vis(i,j)/1.609)**0.86) * 1.609
zen_fac = min(0.1,max(czen(i,j),0.))/ 0.1
if (i.eq.290 .and. j.eq.112) then
write (6,*) 'zen_fac,vis_night, vis =',zen_fac,vis_night, vis(i,j)
end if
vis(i,j) = zen_fac * vis(i,j) + (1.-zen_fac)*vis_night
if (i.eq.290 .and. j.eq.112) then
write (6,*) 'visrh, vis =',visrh, vis(i,j)
end if
vis(i,j) = min(vis(i,j),visrh)
if (vis(i,j).lt.1.) vis1km_cnt = vis1km_cnt + 1
if (vis(i,j).lt.3.) vis3km_cnt = vis3km_cnt + 1
if (vis(i,j).lt.5.) vis5km_cnt = vis5km_cnt + 1
! convert to [m]
vis(i,j) = vis(i,j) * 1000.
ENDDO
ENDDO
write (6,*)
write (6,*) ' Visibility diagnostics follow: ------------'
write (6,*) ' -------------------------------------------'
write (6,*) &
' any vis / vis < 10 km '
write (6,*)'No. of grid pts with shear (lev4-1) > 4m/s', &
shear4_cnt, shear4_cnt_lowvis
write (6,*)'No. of grid pts with shear (lev4-1) > 5m/s', &
shear5_cnt, shear5_cnt_lowvis
write (6,*)'No. of grid pts with shear (lev4-1) > 6m/s', &
shear8_cnt, shear8_cnt_lowvis
write (6,*)
write (6,*)'No. of grid pts with vis-RH < 10 km', &
visrh10_cnt
write (6,*)'No. of grid pts with vis < 1 km', &
vis1km_cnt
write (6,*)'No. of grid pts with vis < 3 km', &
vis3km_cnt
write (6,*)'No. of grid pts with vis < 5 km', &
vis5km_cnt
write (6,*)
write (6,*)'Min vis-hydrometeor, vis-RH', vis_min, visrh_min
write (6,*)'No. of grid pts with visRH < vis(hydrometeor)', &
visrh_lower
write (6,*)'% grid pts with night/cos(zen) < 0.1', &
float(night_cnt)/float(IM*JM),float(lowsun_cnt)/ &
float(IM*JM)
write (6,*)
!
RETURN
END