[Dart-dev] DART/branches Revision: 12832

[dart at ucar.edu]

nancy at ucar.edu
2018-09-18 10:54:31 -0600 (Tue, 18 Sep 2018)
239
updated dewpoint equation which avoids the very slight
chance of dividing by 0.  happened once in 11 years of
use, but it generates NaNs in that case as the result of
the forward operator.  updated docs as well.  code and
docs from jeff.




Modified: DART/branches/rma_trunk/observations/forward_operators/obs_def_dew_point_mod.f90
===================================================================
--- DART/branches/rma_trunk/observations/forward_operators/obs_def_dew_point_mod.f90	2018-09-18 16:51:41 UTC (rev 12831)
+++ DART/branches/rma_trunk/observations/forward_operators/obs_def_dew_point_mod.f90	2018-09-18 16:54:31 UTC (rev 12832)
@@ -142,9 +142,10 @@
 real(r8) :: qv(ens_size)            ! water vapor mixing ratio (kg/kg)
 real(r8) :: e_mb(ens_size)          ! water vapor pressure (mb)
 real(r8),   PARAMETER :: e_min = 0.001_r8 ! threshold for minimum vapor pressure (mb),
-                                          !   to avoid problems near zero in Bolton's equation
+                                          !   to avoid problems near zero vapor pressure in Bolton's equation
 real(r8) :: p_Pa(ens_size)          ! pressure (Pa)
 real(r8) :: p_mb(ens_size)          ! pressure (mb)
+real(r8) :: log_term(ens_size)      ! Intermediate term in computation of dewpoint temperature
 !> @todo make strings longer
 character(len=129) :: errstring
 logical  :: return_now
@@ -193,9 +194,21 @@
 
    !------------------------------------------------------------------------------
    !  Use Bolton's approximation to compute dewpoint.
+   !  Bolton, David, 1980: The Computation of Equivalent Potential Temperature. 
+   !  Monthly Weather Review. 108 (7): 1046-1053.
+   !  Bolton does not explicitly have this formula, but the pieces are there.
+   !  He uses the three constant values used here.
+   !  Other authors, and the Wikipedia page on dewpoint, provide this formula
+   !  explicitly and suggest 6.1121 instead of 6.112.
    !------------------------------------------------------------------------------
 
-   td = t_kelvin + (243.5_r8 / ((17.67_r8 / log(e_mb/6.112_r8)) - 1.0_r8) )
+   ! The following expression can fail numerically for dewpoints very close to 0 C
+   !td = t_kelvin + (243.5_r8 / ((17.67_r8 / log(e_mb/6.112_r8)) - 1.0_r8) )
+
+   ! A numerically robust formula that avoids the failure near dewpoints of 0 C
+   log_term = log(e_mb / 6.112_r8)
+   td = t_kelvin + 243.5_r8 * log_term / (17.67_r8 - log_term)
+
 elsewhere
    td = missing_r8
 end where

Modified: DART/branches/rma_trunk/observations/forward_operators/obs_def_dew_point_mod.html
===================================================================
--- DART/branches/rma_trunk/observations/forward_operators/obs_def_dew_point_mod.html	2018-09-18 16:51:41 UTC (rev 12831)
+++ DART/branches/rma_trunk/observations/forward_operators/obs_def_dew_point_mod.html	2018-09-18 16:54:31 UTC (rev 12832)
@@ -41,6 +41,8 @@
 <br>
 Revision 2801 implements a more robust method 
 (based on Bolton's Approximation) for calculating dew point.
+This has been further revised to avoid a numerical instability that could lead
+to failed forward operators for dewpoints almost exactly 0 C.
 </P>
 
 <!--==================================================================-->
@@ -137,7 +139,7 @@
 <div class="top">[<a href="#">top</a>]</div><hr />
 <H2>REFERENCES</H2>
 <ol>
-<li> none </li>
+<li> Bolton, David, 1980: The Computation of Equivalent Potential Temperature. Monthly Weather Review, 108, 1046-1053. </li>
 </ol>
 
 <!--==================================================================-->