Fixes momentum flux in gustiness parameterization

Merge singhbalwinder/atm/gust-fix (PR #1624)

This bug and its fix is found by Po-Lun Ma. The original code has a
bug,which makes the gustiness effect on momentum flux smaller than
what it should be. The effects of the bug on clouds and precipitation
are small, but has an important effect on the westerly anomalies in
tropical Pacific which is relevant to ENSO. The bug fix enhances the
momentum flux in the tropics, weakening the prevailing easterlies so
that the cold tongue won't extend too far into the west. The largest
signal comes from TWP, so it might also improve WWB to help ENSO.

[BFB] - Bit-For-Bit

src/share/util/shr_flux_mod.F90

@@ -189,6 +189,7 @@ SUBROUTINE shr_flux_atmOcn(nMax  ,zbot  ,ubot  ,vbot  ,thbot ,  prec_gust, gust_
    !--- local variables --------------------------------
    integer(IN) :: n      ! vector loop index
    real(R8)    :: vmag   ! surface wind magnitude   (m/s)
+   real(R8)    :: vmag_old   ! surface wind magnitude without gustiness (m/s)
    real(R8)    :: thvbot ! virtual temperature      (K)
    real(R8)    :: ssq    ! sea surface humidity     (kg/kg)
    real(R8)    :: delt   ! potential T difference   (K)
@@ -281,10 +282,18 @@ SUBROUTINE shr_flux_atmOcn(nMax  ,zbot  ,ubot  ,vbot  ,thbot ,  prec_gust, gust_
         ! old version
         !vmag   = max(umin, sqrt( (ubot(n)-us(n))**2 + (vbot(n)-vs(n))**2) )
 
-        !--- vmag+ugust (convective gustiness) Limit to a max precip 6 cm/day = 0.00069444 mm/s.
+        !--- vmag+ugust (convective gustiness) Limit to a max precip 6 cm/day = 0.00069444 m/s.
         !--- reverts to original formula if gust_fac=0 
-        vmag   = max(umin, sqrt( (ubot(n)-us(n))**2 + (vbot(n)-vs(n))**2) +  ugust(min(prec_gust(n),6.94444e-4_R8)))
 
+        !PMA saves vmag_old for taux tauy computation 
+
+        vmag_old    = max(umin, sqrt( (ubot(n)-us(n))**2 + (vbot(n)-vs(n))**2) )
+
+        if (gust_fac .gt. 1.e-12_R8) then
+         vmag       = max(umin, sqrt( (ubot(n)-us(n))**2 + (vbot(n)-vs(n))**2) + ugust(min(prec_gust(n),6.94444e-4_R8)))
+        else
+         vmag       = vmag_old
+        endif
 
         thvbot = thbot(n) * (1.0_R8 + loc_zvir * qbot(n)) ! virtual temp (K)
         ssq    = 0.98_R8 * qsat(ts(n)) / rbot(n)   ! sea surf hum (kg/kg)
@@ -377,8 +386,8 @@ SUBROUTINE shr_flux_atmOcn(nMax  ,zbot  ,ubot  ,vbot  ,thbot ,  prec_gust, gust_
         tau = rbot(n) * ustar * ustar 
 
         !--- momentum flux ---
-        taux(n) = tau * (ubot(n)-us(n)) / vmag 
-        tauy(n) = tau * (vbot(n)-vs(n)) / vmag 
+        taux(n) = tau * (ubot(n)-us(n)) / vmag_old !PMA uses vmag_old for taux 
+        tauy(n) = tau * (vbot(n)-vs(n)) / vmag_old !    tauy c20170620
 
         !--- heat flux ---
         sen (n) =          cp * tau * tstar / ustar