Reintroducing visibility including blowing snow effect. (#59)

NOAA-GSL · Oct 4, 2023 · 4db9f8f · 4db9f8f
1 parent ed689d4
commit 4db9f8f
Show file tree

Hide file tree

Showing 5 changed files with 169 additions and 6 deletions.
diff --git a/parm/fv3lam_rrfs.xml b/parm/fv3lam_rrfs.xml
@@ -925,6 +925,12 @@
        <scale>6.0</scale>
        </param>
 
+       <param>
+       <shortname>GSD_BLSN_ON_SURFACE</shortname>
+       <pname>VISBSN</pname>
+       <scale>6.0</scale>
+       </param>
+
        <param>
        <shortname>HGT_ON_LVL_OF_ADIAB_COND_FROM_SFC</shortname>
        <pname>HGT</pname>
@@ -3406,6 +3412,12 @@
        <scale>6.0</scale>
        </param>
 
+       <param>
+       <shortname>GSD_BLSN_ON_SURFACE</shortname>
+       <pname>VISBSN</pname>
+       <scale>6.0</scale>
+       </param>
+
        <param>
        <shortname>REFC_ON_ENTIRE_ATMOS</shortname>
        <pname>REFC</pname>

diff --git a/parm/post_avblflds.xml b/parm/post_avblflds.xml
@@ -3857,6 +3857,15 @@
 
 <!--438-442-->
 
+    <param>
+       <post_avblfldidx>438</post_avblfldidx>
+       <shortname>GSD_BLSN_ON_SURFACE</shortname>
+       <longname>GSD_blsn on surface</longname>
+       <pname>VISBSN</pname>
+       <fixed_sfc1_type>surface</fixed_sfc1_type>
+       <scale>6.0</scale>
+    </param>
+
 <!--443-->
     <param>
        <post_avblfldidx>443</post_avblfldidx>

diff --git a/parm/postxconfig-NT-fv3lam_rrfs.txt b/parm/postxconfig-NT-fv3lam_rrfs.txt
@@ -1,7 +1,7 @@
 3
 3
-249
-306
+250
+307
 PRSLEV
 32769
 ncep_nco
@@ -4828,6 +4828,43 @@ surface
 ?
 ?
 ?
+438
+GSD_BLSN_ON_SURFACE
+GSD_blsn on surface
+1
+tmpl4_0
+VISBSN
+?
+?
+surface
+0
+?
+0
+?
+?
+0
+?
+0
+?
+?
+?
+0
+0.0
+0
+0.0
+?
+0
+0.0
+0
+0.0
+1
+6.0
+0
+0
+0
+?
+?
+?
 109
 HGT_ON_LVL_OF_ADIAB_COND_FROM_SFC
 ?
@@ -18460,6 +18497,43 @@ surface
 ?
 ?
 ?
+438
+GSD_BLSN_ON_SURFACE
+GSD_blsn on surface
+1
+tmpl4_0
+VISBSN
+?
+?
+surface
+0
+?
+0
+?
+?
+0
+?
+0
+?
+?
+?
+0
+0.0
+0
+0.0
+?
+0
+0.0
+0
+0.0
+1
+6.0
+0
+0
+0
+?
+?
+?
 252
 REFC_ON_ENTIRE_ATMOS
 ?

diff --git a/sorc/ncep_post.fd/CALVIS_GSD.f b/sorc/ncep_post.fd/CALVIS_GSD.f
@@ -1,5 +1,5 @@
 !**********************************************************************c
-      SUBROUTINE CALVIS_GSD(CZEN,VIS)
+      SUBROUTINE CALVIS_GSD(CZEN,VIS,METHOD_BLSN)
 
 ! SUBPROGRAM:    CALVIS      CALCULATE HORIZONTAL VISIBILITY   
 !
@@ -91,24 +91,30 @@ SUBROUTINE CALVIS_GSD(CZEN,VIS)
 !   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
+!   2023-02        Tim Corrie, Eric James - addition of attenuation for blowing snow
 !                           
 !------------------------------------------------------------------
 !
 
+      use vrbls2d, only: sno, si, ustar
       use vrbls3d, only: qqw, qqi, qqs, qqr, qqg, t, pmid, q, u, v, extcof55, aextc55
-      use params_mod, only: h1, d608, rd
+      use params_mod, only: h1, d608, rd, g
       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
+      integer :: method, method_blsn
       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 :: z, ustar_t, u_p, lamda, r_bar, alpha
+      real :: rho_sno
+      real :: z_r, Q_s, C_r, c_z, c_alpha, vis_blsn, BETABLSN
+
       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,  &
@@ -136,6 +142,7 @@ SUBROUTINE CALVIS_GSD(CZEN,VIS)
 !                Aerosol method (2), 
 !                Smoke added to RH method for clear air  (3)
 !                   3 - option to add reducted visibility from smoke-based aerosols.
+! method_blsn = 0 ! 1 for 'on', 0 for 'off'
      
       CELKEL     = 273.15
       TICE       = CELKEL-10.
@@ -300,6 +307,45 @@ SUBROUTINE CALVIS_GSD(CZEN,VIS)
           if (t(i,j,lm)< 270. .and. temp_fac==1.)          &
              write (6,*) 'Problem w/ temp_fac - calvis'
 
+! Key calculation of attenuation from blowing snow -- updated 5 August 2022 by Tim Corrie
+! Framework is from Letcher et al (2021)
+
+        ustar_t = 0.2
+        u_p = 2.8*ustar_t
+        lamda = 0.45
+        z = 2.0
+        alpha = 15.0
+        r_bar = 0.0002
+
+!        print *, i,j
+
+
+        if (si(i,j)<spval .and. si(i,j) .ge. 1.0) then
+            z_r = 1.6*(ustar(i,j)**2./(2.*g))
+            Q_s = max((0.68/ustar(i,j))*(RHOAIR/g)*(ustar(i,j)**2.-ustar_t**2.),0.0)
+            C_r = (Q_s/u_p)*(lamda*g/ustar(i,j)**2.)*exp(-lamda*z_r*g/ustar(i,j)**2.)
+            c_z = max(C_r * exp(-1.55*((0.05628*ustar(i,j))**-0.544 - z**-0.544)),1e-15)
+            c_alpha = alpha/(alpha+2) !simplified version of (6) in Letcher et al (2021)    
+            rho_sno = sno(i,j)/(si(i,j)/1.0e3)
+            rho_sno = rho_sno*2. + 10.*max(0.,rho_sno-0.15)
+            vis_blsn = (5.217*rho_sno*r_bar**1.011)/(1.82*c_z*c_alpha)
+            BETABLSN = 3.912/(vis_blsn/1000.0)
+            ! print to ensure quality
+            !print *, "z_r", z_r
+            !print *, "Q_s", Q_s
+            !print *, "C_r", C_r
+            !print *, "c_z", c_z
+            !print *, "c_alpha", c_alpha
+            !print *, "sno/SWE", sno(i,j)
+            !print *, "si/SNOD", si(i,j)/1.0e3
+            !print *, "rho_sno", rho_sno
+            !print *, "vis_blsn", vis_blsn
+            !print *, "BETABLSN", BETABLSN
+        else
+            BETABLSN = 0
+            !print *, "BETABLSN", BETABLSN
+        end if
+
 ! Key calculation of attenuation from each hydrometeor type (cloud, snow, graupel, rain, ice)
         BETAV=COEFFC*CONCFC**EXPONFC                            &
              + coef_SNOW*CONCFP**EXPONFP                        &
@@ -309,6 +355,9 @@ SUBROUTINE CALVIS_GSD(CZEN,VIS)
 ! Addition of attenuation from aerosols if option selected
         if(method == 2 .or. method == 3)then ! aerosol method
             BETAV = BETAV + aextc55(i,j,lm)*1000.
+            if(method_blsn .eq. 1) then ! BLSN method, updated 8 August 2022 by Tim Corrie
+                BETAV = BETAV + BETABLSN
+            endif
         endif
 
 !  Calculation of visibility based on hydrometeor and aerosols.  (RH effect not yet included.)

diff --git a/sorc/ncep_post.fd/MDLFLD.f b/sorc/ncep_post.fd/MDLFLD.f
@@ -3550,7 +3550,7 @@ SUBROUTINE MDLFLD
 ! --- GSD VISIBILITY
 !
       IF (IGET(410)>0) THEN
-        CALL CALVIS_GSD(CZEN,VIS)
+        CALL CALVIS_GSD(CZEN,VIS,0)
         DO J=JSTA,JEND
         DO I=ista,iend
           GRID1(I,J)=VIS(I,J)
@@ -3563,6 +3563,25 @@ SUBROUTINE MDLFLD
          datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=GRID1(ista:iend,jsta:jend)
         endif
       ENDIF
+
+!
+! --- GSD VISIBILITY WITH BLSN
+!
+      IF (IGET(438)>0) THEN
+        CALL CALVIS_GSD(CZEN,VIS,1)
+        DO J=JSTA,JEND
+        DO I=ista,iend
+          GRID1(I,J)=VIS(I,J)
+        END DO
+        END DO
+        if(grib=="grib2") then
+         cfld=cfld+1
+         fld_info(cfld)%ifld=IAVBLFLD(IGET(438))
+         fld_info(cfld)%lvl=LVLSXML(1,IGET(438))
+         datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=GRID1(ista:iend,jsta:jend)
+        endif
+      ENDIF
+
 !
 ! --- RADAR REFLECT - 1km
 !