CALVIS_GSD.f

!**********************************************************************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%)
!   aextc55--aerosol extinction coefficient (m**-1)
!
!
!   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)]
!
!   The 'aextc55' field is 3-D and is derived from the 'aod_3d' field
!   by dividing by dz, the vertical thickness of that model level in 'm'.
!   This can be handled as a 2-D field if needed to save resources.
!
! 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
!   2015-17        S. Benjamin, T. Smirnova - modifications for RH-based clear-air vis
!   2017-12        R. Ahmadov, Steve Albers - addition for attenuation from aerosols
!                              (not related to water vapor or RH at this point)
!   2021-05        Wen Meng        Unify CONST1 and VISRH. 
!   2021-05        Wen Meng  - Add checking for undefined points invloved in computation
!   2021-08        Wen Meng  - Restrict divided by 0.
!                           
!------------------------------------------------------------------
!

      use vrbls3d, only: qqw, qqi, qqs, qqr, qqg, t, pmid, q, u, v, extcof55, aextc55
      use params_mod, only: h1, d608, rd
      use ctlblk_mod, only: jm, im, jsta_2l, jend_2u, lm, modelname, spval

      implicit none

      integer :: j, i, k, ll
      integer :: method
      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
!------------------------------------------------------------------

! Method used for clear-air visibility with extension for aerosols
      method = 3 
!                RH-only method (1), 
!                Aerosol method (2), 
!                Smoke added to RH method for clear air  (3)
!                   3 - option to add reducted visibility from smoke-based aerosols.
     
      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

! - snow-based vis attenuation - coefficient 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
!    Rasmussen et al. 2003,  J. App. Meteor.
!    Snow Nowcasting Using a Real-Time Correlation of Radar Reflectivity with Snow Gauge Accumulation
      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
      const1= 3.000

! visibility with respect to RH is
!   calculated from optical depth linearly
!   related to RH as follows:

!    vis = 60 exp (-2.5 * (RH-15)/80)
!       changed on 3/14/01

!  coefficient of 3 gives visibility of 5 km
!    at 95% RH

! Total visibility is minimum of vis-rh  (developed by Benjamin, Brown, Smirnova)
!   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
        vis(i,j)=spval
! -checking undedined points
        if(t(i,j,lm)<spval .and. u(i,j,lm)<spval .and. v(i,j,lm)<spval &
           .and. pmid(i,j,lm)<spval) then 
!  - take max hydrometeor mixing ratios in lowest 3 levels (lowest 13 hPa, 100m with RAP/HRRR
      qrain = 0.
      qsnow = 0.
      qgraupel = 0.
      qclw = 0.
      qclice = 0.

          do k = 1,3
               ll=lm-k+1
            if(qqw(i,j,ll)<spval)qclw = max(qclw, qqw(i,j,ll) )
            if(qqi(i,j,ll)<spval)qclice = max(qclice, qqi(i,j,ll) )
            if(qqs(i,j,ll)<spval)qsnow = max(qsnow, qqs(i,j,ll) )
            if(qqr(i,j,ll)<spval)qrain = max(qrain, qqr(i,j,ll) )
            if(qqg(i,j,ll)<spval)qgraupel  = max(qgraupel, qqg(i,j,ll) )
! - compute relative humidity
        tx=t(i,j,ll)-273.15
        rhb(i,j,ll)=0.
        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)))))))))
        if(abs(pol) > 0.) THEN
        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)
       ENDIF

          enddo

!  - take max RH of levels 1 and 2 near the sfc
          rhmax = max(rhb(i,j,lm),rhb(i,j,lm-1))
          qrh = max(0.0,min(0.8,(rhmax/100.-0.15)))

!tgs 23 feb 2017 - increase of base value to 90 km to reduce attenuation
!                  from RH for clear-air visibility.  (i.e., increase clear-air vis overall)
       visrh = 90. * 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 25 hPa for HRRR/RAP)
         shear = sqrt( (u(i,j,lm-3)-u(i,j,lm))**2         &
                     +(v(i,j,lm-3)-v(i,j,lm))**2  )

        shear_fac = min(1.,max(0.,(shear-4.)/2.) )
        if (visrh<10.) visrh = visrh + (10.-visrh)*    &
           shear_fac

        if (shear>4.) shear4_cnt = shear4_cnt +1
        if (shear>5.) shear5_cnt = shear5_cnt +1
        if (shear>6.) shear8_cnt = shear8_cnt +1

        if (shear>4..and.visrh<10)                  &
          shear4_cnt_lowvis = shear4_cnt_lowvis +1
        if (shear>5..and.visrh<10)                  &
          shear5_cnt_lowvis = shear5_cnt_lowvis +1
        if (shear>6..and.visrh<10)                  &
          shear8_cnt_lowvis = shear8_cnt_lowvis +1

        if (visrh<10.) visrh10_cnt = visrh10_cnt+1
        if (czen(i,j)<0.) night_cnt = night_cnt + 1
        if (czen(i,j)<0.1) lowsun_cnt = lowsun_cnt + 1

        tv=t(i,j,lm)*(h1+d608*q(i,j,lm))

        rhoair=pmid(i,j,lm)/(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.) )

         coef_snow = coeffp_dry*(1.-temp_fac)                   &
                   + coeffp_wet* temp_fac    

          if (t(i,j,lm)< 270. .and. temp_fac==1.)          &
             write (6,*) 'Problem w/ temp_fac - calvis'

! Key calculation of attenuation from each hydrometeor type (cloud, snow, graupel, rain, ice)
        betav=coeffc*concfc**exponfc                            &
             + coef_snow*concfp**exponfp                        &
             + coeflc*conclc**exponlc + coeflp*conclp**exponlp    &
             + coeffg*concfg**exponfg  +1.e-10

! Addition of attenuation from aerosols if option selected
        if(method == 2 .or. method == 3)then ! aerosol method
            betav = betav + aextc55(i,j,lm)*1000.
        endif

       if (i==290 .and. j==112) then
         write (6,*) 'BETAV, extcof55 =',betav,extcof55(i,j,lm)
       end if

!  Calculation of visibility based on hydrometeor and aerosols.  (RH effect not yet included.)
        vis(i,j)=min(90.,const1/(betav+extcof55(i,j,lm)))      ! max of 90km

        if (vis(i,j)<vis_min) vis_min = vis(i,j)
        if (visrh<visrh_min) visrh_min = visrh

        if (visrh<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==290 .and. j==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==290 .and. j==112) then
         write (6,*) 'visrh, vis =',visrh, vis(i,j)
       end if

       if(method == 1 .or. method == 3)then ! RH method (if lower vis)
         vis(i,j) = min(vis(i,j),visrh)
       endif

        if (vis(i,j)<1.) vis1km_cnt = vis1km_cnt + 1
        if (vis(i,j)<3.) vis3km_cnt = vis3km_cnt + 1
        if (vis(i,j)<5.) vis5km_cnt = vis5km_cnt + 1
! convert vis from km to [m]
        vis(i,j) = vis(i,j) * 1000.

        endif !end checking undefined points
      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