[Dart-dev] DART/branches Revision: 12686

[dart at ucar.edu]

jla at ucar.edu
2018-06-27 14:55:57 -0600 (Wed, 27 Jun 2018)
91
Added check for failed interpolation and pass back missing_r8 for those ensemble members. 



Modified: DART/branches/gigg/observations/forward_operators/obs_def_1d_state_mod.f90
===================================================================
--- DART/branches/gigg/observations/forward_operators/obs_def_1d_state_mod.f90	2018-06-27 20:16:28 UTC (rev 12685)
+++ DART/branches/gigg/observations/forward_operators/obs_def_1d_state_mod.f90	2018-06-27 20:55:57 UTC (rev 12686)
@@ -203,7 +203,9 @@
 ! for tracking/debug use if needed.
 
 if (is_ascii) then
-   write(ifile, *) power(powkey)
+   write(ifile, 11) power(powkey)
+   ! Avoid having an integer entered value get truncated to non-integer which would be disastrous for power
+   11 format(G25.16)
    write(ifile, *) powkey
 else
    write(ifile)    power(powkey)
@@ -394,7 +396,7 @@
 real(r8)              :: range, loc, bottom, dx, x
 integer               :: j !< loop variable
 integer               :: point_istatus(ens_size)
-logical               :: return_now ! used to return early if an interpoaltion fails
+logical               :: return_now ! used to return early if an interpolation fails
 real(r8), dimension(ens_size) :: sum, dist, weight, weight_sum
 
 if ( .not. module_initialized ) call initialize_module
@@ -492,17 +494,24 @@
 
 pow = power(powkey)
 
-! NOTE: WORRY ABOUT SMALL ROUND_OFF AWAY FROM INTEGER IN FILES?
 if(pow == int(pow)) then
    ! Integer power, just use standard definition
-   val = pval ** pow
+   where(istatus == 0) 
+      val = pval ** pow
+   elsewhere
+      val = missing_r8
+   endwhere
 else
    ! For non-integer powers, fix up values for negative bases
    do i = 1, ens_size
-      if(pval(i) >= 0.0_r8) then
-         val(i) = pval(i) ** pow
+      if(istatus(i) == 0) then
+         if(pval(i) >= 0.0_r8) then
+            val(i) = pval(i) ** pow
+         else
+            val(i) = - (- pval(i)) ** pow
+         endif
       else
-         val(i) = -1.0_r8 * (-1.0_r8 * pval(i)) ** pow
+         val(i) = missing_r8
       endif
    end do
 endif