CALVIS.f

!**********************************************************************c
      SUBROUTINE calvis(QV,QC,QR,QI,QS,TT,PP,VIS)
!
!   This routine computes horizontal visibility at the
!   surface or lowest model layer, from qc, qr, qi, and qs.  
!   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)
!   tt--temperature              (k)
!   pp--pressure                 (Pa)
!
!   If iice=0:
!      qprc=qr     qrain=qr and qclw=qc if T>0C
!      qcld=qc          =0          =0  if T<0C
!                  qsnow=qs and qclice=qc  if T<0C
!                       =0            =0   if T>0C
!   If iice=1:
!      qprc=qr+qs   qrain=qr and qclw=qc
!      qcld=qc+qi   qsnow=qs and qclice=qc
!
!   Independent of the above definitions, the scheme can use different
!   assumptions of the state of hydrometeors:
!        meth='d': qprc is all frozen if T<0, liquid if T>0
!        meth='b': Bocchieri scheme used to determine whether qprc
!           is rain or snow. A temperature assumption is used to
!           determine whether qcld is liquid or frozen.
!        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)
!
!   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
!
!   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)]
!
! 2021-05  Wen Meng  -Add checking for undfined points invloved in 
!                     computation.
!------------------------------------------------------------------
    use params_mod, only: h1, d608, rd
    use ctlblk_mod, only: jsta, jend, im, jsta_2l, jend_2u, spval
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    implicit none
!

      real,dimension(IM,jsta_2l:jend_2u),intent(in)    :: qv,qc,qr,qi,qs,tt,pp
      real,dimension(IM,jsta_2l:jend_2u),intent(inout) :: vis
 
      CHARACTER meth*1
      real celkel,tice,coeflc,coeflp,coeffc,coeffp,exponlc,      &
           exponlp,const1,rhoice,rhowat,qprc,qcld,qrain,qsnow,   &
           qclw,qclice,tv,rhoair,vovermd,conclc,concld,concfc,   &
           concfd,betav,exponfc,exponfp,conclp,concfp
      integer i,j
!------------------------------------------------------------------
!------------------------------------------------------------------
      celkel=273.15
      tice=celkel-10.
      coeflc=144.7
      coeflp=2.24
      coeffc=327.8
      coeffp=10.36
      exponlc=0.8800
      exponlp=0.7500
      exponfc=1.0000
      exponfp=0.7776
      const1=-log(.02)
      rhoice=917.
      rhowat=1000.
!
      DO j=jsta,jend
      DO i=1,im
        vis(i,j)=spval
!       IF(IICE==0)THEN
!         QPRC=QR
!         QCLD=QC
!         IF(TT<CELKEL)THEN
!           QRAIN=0.
!           QSNOW=QPRC
!           QCLW=0.
!           QCLICE=QCLD
!         ELSE
!           QRAIN=QPRC
!           QSNOW=0.
!           QCLW=QCLD
!           QCLICE=0.
!         ENDIF
!       ELSE
       IF (qr(i,j) < spval .and. qs(i,j) < spval .and. &
           qc(i,j) < spval .and. qi(i,j) < spval .and. &
           tt(i,j) < spval .and. qv(i,j) < spval .and. &
           pp(i,j) < spval) THEN
          qprc=qr(i,j)+qs(i,j)
          qcld=qc(i,j)+qi(i,j)
          qrain=qr(i,j)
          qsnow=qs(i,j)
          qclw=qc(i,j)
          qclice=qi(i,j)
!       ENDIF
!       TV=VIRTUAL(TT,QV)
        tv=tt(i,j)*(h1+d608*qv(i,j))
        rhoair=pp(i,j)/(rd*tv)
!       IF(METH=='D')THEN
!         IF(TT<CELKEL)THEN
!           VOVERMD=(1.+QV)/RHOAIR+(QPRC+QCLD)/RHOICE
!           CONCLC = 0.
!           CONCLP = 0.
!           CONCFC = QCLD/VOVERMD*1000.
!           CONCFP = QPRC/VOVERMD*1000.
!         ELSE
!           VOVERMD=(1.+QV)/RHOAIR+(QPRC+QCLD)/RHOWAT
!           CONCLC = QCLD/VOVERMD*1000.
!           CONCLP = QPRC/VOVERMD*1000.
!           CONCFC = 0.
!           CONCFP = 0.
!         ENDIF
!       ELSEIF(METH=='B')THEN
!         IF(TT<TICE)THEN
!           VOVERMD=(1.+QV)/RHOAIR+(QPRC+QCLD)/RHOICE
!           CONCLC = 0.
!           CONCLP = 0.
!           CONCFC = QCLD/VOVERMD*1000.
!           CONCFP = QPRC/VOVERMD*1000.
!         ELSEIF(PRSNOW>=50.)THEN
!           VOVERMD=(1.+QV)/RHOAIR+QPRC/RHOICE+QCLD/RHOWAT
!           CONCLC = QCLD/VOVERMD*1000.
!           CONCLP = 0.
!           CONCFC = 0.
!           CONCFP = QPRC/VOVERMD*1000.
!         ELSE
!           VOVERMD=(1.+QV)/RHOAIR+(QPRC+QCLD)/RHOWAT
!           CONCLC = QCLD/VOVERMD*1000.
!           CONCLP = QPRC/VOVERMD*1000.
!           CONCFC = 0.
!           CONCFP = 0.
!         ENDIF
!       ELSEIF(METH=='R')THEN
          vovermd=(1.+qv(i,j))/rhoair+(qclw+qrain)/rhowat+        &
                  (qclice+qsnow)/rhoice
          conclc = max(0., qclw/vovermd*1000.)
          conclp = max(0., qrain/vovermd*1000.)
          concfc = max(0., qclice/vovermd*1000.)
          concfp = max(0., qsnow/vovermd*1000.)
!       ENDIF
        betav=coeffc*concfc**exponfc+coeffp*concfp**exponfp        &
             +coeflc*conclc**exponlc+coeflp*conclp**exponlp        &
             +1.e-10
! CHANGED GSM 3-10-00 -->  no point in distinguishing values
!       above 20 km, so make that value the max (prev max was 80)
!        VIS(I,J)=1.E3*MIN(20.,CONST1/BETAV)   ! max of 20km
! Chuang: Per Geoff, the max visibility was changed to be cosistent with visibility ceiling in obs
        vis(i,j) = 1.e3*min(24.135,const1/betav)   ! change max to be consistent with obs
        ENDIF
      ENDDO
      ENDDO
!
      RETURN
      END