[Dart-dev] DART/branches Revision: 12831

[dart at ucar.edu]

nancy at ucar.edu
2018-09-18 10:51:41 -0600 (Tue, 18 Sep 2018)
176
updated code from jeff - numerically stable computation
and updated documentation.  note this is on the recam
branch and will be moving to rma_trunk and then on to
manhattan.




Modified: DART/branches/recam/observations/forward_operators/obs_def_dew_point_mod.f90
===================================================================
--- DART/branches/recam/observations/forward_operators/obs_def_dew_point_mod.f90	2018-09-18 14:07:08 UTC (rev 12830)
+++ DART/branches/recam/observations/forward_operators/obs_def_dew_point_mod.f90	2018-09-18 16:51:41 UTC (rev 12831)
@@ -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/recam/observations/forward_operators/obs_def_dew_point_mod.html
===================================================================
--- DART/branches/recam/observations/forward_operators/obs_def_dew_point_mod.html	2018-09-18 14:07:08 UTC (rev 12830)
+++ DART/branches/recam/observations/forward_operators/obs_def_dew_point_mod.html	2018-09-18 16:51:41 UTC (rev 12831)
@@ -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>
 
 <!--==================================================================-->