UPP/upp-srw-v2.2.0/CALVIS__GSD_8f_source.html

!**********************************************************************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.

!   2021-10        Jesse Meng - 2D DECOMPOSITION

!

!------------------------------------------------------------------

!


      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,&

                                    ista_2l, iend_2u


      implicit none


      integer :: j, i, k, ll

      integer :: method

      real :: tx, pol, esx, es, e

      REAL VIS(ista_2l:iend_2u,jsta_2l:jend_2u)

      REAL RHB(ista_2l:iend_2u,jsta_2l:jend_2u,LM)

      REAL CZEN(ista_2l:iend_2u,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=ista_2l,iend_2u

        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


!  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

       vis(i,j) = zen_fac * vis(i,j) + (1.-zen_fac)*vis_night


       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