Merge pull request #245 from bensonr/2019.01_quadfix

Fix quad precision issue with hard on/off switch
NOAA-GFDL · Dec 23, 2019 · a722a24 · a722a24
2 parents 708b8d5 + e5abd22
commit a722a24
Show file tree

Hide file tree

Showing 4 changed files with 15 additions and 38 deletions.
diff --git a/include/fms_platform.h b/include/fms_platform.h
@@ -107,13 +107,14 @@ use,intrinsic :: iso_c_binding, only: c_double,c_float,c_int64_t, &
 
 
 !If you want to use quad-precision.
-#if defined(QUAD_PRECISION) && defined(__PGI)
-#error "cannot define QUAD_PRECISION and use PGI compiler"
-#elif !defined(QUAD_PRECISION) && defined(__PGI)
+! The NO_QUAD_PRECISION macro will be deprecated and removed at some future time.
+! Model code will rely solely upon the ENABLE_QUAD_PRECISION macro thereafer.
+#if defined(ENABLE_QUAD_PRECISION) 
+#undef NO_QUAD_PRECISION
+#else
 #define NO_QUAD_PRECISION
+#undef QUAD_KIND
 #define QUAD_KIND DOUBLE_KIND
-#else
-#undef NO_QUAD_PRECISION
 #endif
 
 

diff --git a/interpolator/interpolator.F90 b/interpolator/interpolator.F90
@@ -362,12 +362,12 @@ module interpolator_mod
 real ::  missing_value = -1.e10                                                       !< No description
 ! sjs integer :: itaum, itaup
 
-#ifdef NO_QUAD_PRECISION
-! 64-bit precision (kind=8)
- integer, parameter:: f_p = selected_real_kind(15)     !< 64-bit precision (kind=8)
-#else
+#ifdef ENABLE_QUAD_PRECISION
 ! Higher precision (kind=16) for grid geometrical factors:
  integer, parameter:: f_p = selected_real_kind(20)     !< Higher precision (kind=16) for grid geometrical factors
+#else
+! 64-bit precision (kind=8)
+ integer, parameter:: f_p = selected_real_kind(15)     !< 64-bit precision (kind=8)
 #endif
 
 logical :: read_all_on_init = .false.          !< No description

diff --git a/mosaic/create_xgrid.c b/mosaic/create_xgrid.c
@@ -2424,15 +2424,15 @@ double poly_ctrlon(const double x[], const double y[], int n, double clon)
     if( dphi2 > M_PI) dphi2 -= 2.0*M_PI;
     if( dphi2 <-M_PI) dphi2 += 2.0*M_PI;
 
-    if(abs(dphi2 -dphi1) < M_PI) {
+    if(fabs(dphi2 -dphi1) < M_PI) {
       ctrlon -= dphi * (dphi1*f1+dphi2*f2)/2.0;
     }
     else {
       if(dphi1 > 0.0)
         fac = M_PI;
       else
         fac = -M_PI;
-      fint = f1 + (f2-f1)*(fac-dphi1)/abs(dphi);
+      fint = f1 + (f2-f1)*(fac-dphi1)/fabs(dphi);
       ctrlon -= 0.5*dphi1*(dphi1-fac)*f1 - 0.5*dphi2*(dphi2+fac)*f2
           + 0.5*fac*(dphi1+dphi2)*fint;
     }

diff --git a/mosaic/mosaic_util.c b/mosaic/mosaic_util.c
@@ -541,31 +541,7 @@ double great_circle_area(int n, const double *x, const double *y, const double *
 double spherical_angle(const double *v1, const double *v2, const double *v3)
 {
   double angle;
-#ifdef NO_QUAD_PRECISION
-  double px, py, pz, qx, qy, qz, ddd;
-#ifndef SQRT_
-#define SQRT_ sqrt
-#else
-#error "SQRT_ Previously Defined"
-#endif  /* SQRT_ */
-#ifndef ABS_
-#define ABS_ fabsl
-#else
-#error "ABS_ Previously Defined"
-#endif  /* ABS_ */
-#else
   long double px, py, pz, qx, qy, qz, ddd;
-#ifndef SQRT_
-#define SQRT_ sqrtl
-#else
-#error "SQRT_ Previously Defined"
-#endif  /* SQRT_ */
-#ifndef ABS_
-#define ABS_ fabs
-#else
-#error "ABS_ Previously Defined"
-#endif  /* ABS_ */
-#endif
 
   /* vector product between v1 and v2 */
   px = v1[1]*v2[2] - v1[2]*v2[1];
@@ -580,9 +556,9 @@ double spherical_angle(const double *v1, const double *v2, const double *v3)
   if ( ddd <= 0.0 )
     angle = 0. ;
   else {
-    ddd = (px*qx+py*qy+pz*qz) / SQRT_(ddd);
-    if( ABS_(ddd-1) < EPSLN30 ) ddd = 1;
-    if( ABS_(ddd+1) < EPSLN30 ) ddd = -1;
+    ddd = (px*qx+py*qy+pz*qz) / sqrtl(ddd);
+    if( fabsl(ddd-1) < EPSLN30 ) ddd = 1;
+    if( fabsl(ddd+1) < EPSLN30 ) ddd = -1;
     if ( ddd>1. || ddd<-1. ) {
       /*FIX (lmh) to correctly handle co-linear points (angle near pi or 0) */
       if (ddd < 0.)