Philippe 27/10/2016: OpenACC: use BITREP library to get same results on CPU and GPU

src/MNH/turb.f90

@@ -372,6 +372,10 @@ USE MODI_ETHETA
 !
 USE MODI_SECOND_MNH
 !
+#ifdef MNH_BITREP
+USE MODI_BITREP
+#endif
+!
 !!  use, intrinsic :: ISO_C_BINDING
 !
 IMPLICIT NONE
@@ -662,20 +666,20 @@ END DO
 DO JRR = 2+KRRL,1+KRRL+KRRI                ! loop on the solid components   
   ZCP(:,:,:)  = ZCP(:,:,:)  + XCI * PRT(:,:,:,JRR)
 END DO
-!$acc end kernels
 !
 !*      2.2 Exner function at t
 !
 !PW: "BUG" PGI : results different on CPU and GPU due to the power function
 !See http://www.pgroup.com/userforum/viewtopic.php?t=5364&sid=ec7b762b17fb9bb3332a47f0db57af55
-!!$acc kernels
+!Use of own functions allow bit-reproducible results
+#ifndef MNH_BITREP
 ZEXN(:,:,:) = (PPABST(:,:,:)/XP00) ** (XRD/XCPD)
-!!$acc end kernels
-!$acc update device(ZEXN)
+#else
+ZEXN(:,:,:) = BR_POW(PPABST(:,:,:)/XP00,XRD/XCPD)
+#endif
 !
 !*      2.3 dissipative heating coeff a t
 !
-!$acc kernels
 ZCOEF_DISS(:,:,:) = 1/(ZCP(:,:,:) * ZEXN(:,:,:)) 
 !
 !
@@ -1515,21 +1519,20 @@ REAL, DIMENSION(SIZE(PEXN,1),SIZE(PEXN,2),SIZE(PEXN,3)) :: ZDRVSATDT
 !
 !$acc kernels
   PLOCPEXN(:,:,:) = ( PLTT + (XCPV-PC) *  (PT(:,:,:)-XTT) ) / PCP(:,:,:)
-!$acc end kernels
 !
 !*      1.2 Saturation vapor pressure at t
 !
 !PW: "BUG" PGI : results different on CPU and GPU due to the EXP and LOG functions
 !See http://www.pgroup.com/userforum/viewtopic.php?t=5364&sid=ec7b762b17fb9bb3332a47f0db57af55
-!$acc update self(PT)
-!!!$acc kernels
+!Use of own functions allow bit-reproducible results
+#ifndef MNH_BITREP
   ZRVSAT(:,:,:) =  EXP( PALP - PBETA/PT(:,:,:) - PGAM*ALOG( PT(:,:,:) ) )
-!!!$acc end kernels
+#else
+  ZRVSAT(:,:,:) =  BR_EXP( PALP - PBETA/PT(:,:,:) - PGAM*BR_LOG( PT(:,:,:) ) )
+#endif
 !
 !*      1.3 saturation  mixing ratio at t
 !
-!$acc update device(ZRVSAT)
-!$acc kernels
   ZRVSAT(:,:,:) =  ZRVSAT(:,:,:) * ZEPS / ( PPABST(:,:,:) - ZRVSAT(:,:,:) )
 !
 !*      1.4 compute the saturation mixing ratio derivative (rvs')
@@ -1696,25 +1699,33 @@ IF ( HTURBDIM /= '1DIM' ) THEN  ! 3D turbulence scheme
 !$acc end kernels
 #endif
   ELSE
+!PW: "BUG" PGI : results different on CPU and GPU due to the power function
+!See http://www.pgroup.com/userforum/viewtopic.php?t=5364&sid=ec7b762b17fb9bb3332a47f0db57af55
+!Use of own functions allow bit-reproducible results
+#ifndef MNH_BITREP
+!
 #ifndef _OPENACC
     PLM(:,:,:) = (PLM(:,:,:)*MXF(PDXX(:,:,:))*MYF(PDYY(:,:,:)) ) ** (1./3.)
 #else
-#if 0
-!PW: "BUG" PGI : results different on CPU and GPU due to the power function
-!See http://www.pgroup.com/userforum/viewtopic.php?t=5364&sid=ec7b762b17fb9bb3332a47f0db57af55
     CALL MXF_DEVICE(PDXX,ZTMP1_DEVICE)
     CALL MYF_DEVICE(PDYY,ZTMP2_DEVICE)
-!$acc update self(ZTMP1_DEVICE,ZTMP2_DEVICE)
-!!$acc kernels
+!$acc kernels
     PLM(:,:,:) = (PLM(:,:,:)*ZTMP1_DEVICE*ZTMP2_DEVICE ) ** (1./3.)
-!!$acc end kernels
-!$acc update device(PLM)
+!$acc end kernels
+#endif !_OPENACC
+!
+#else !MNH_BITREP
+!
+#ifndef _OPENACC
+    PLM(:,:,:) = BR_POW( PLM(:,:,:)*MXF(PDXX(:,:,:))*MYF(PDYY(:,:,:)) , 1./3. )
 #else
-!$acc update self(PLM)
-    PLM(:,:,:) = (PLM(:,:,:)*MXF(PDXX(:,:,:))*MYF(PDYY(:,:,:)) ) ** (1./3.)
-!$acc update device(PLM)
-#endif
-#endif
+    CALL MXF_DEVICE(PDXX,ZTMP1_DEVICE)
+    CALL MYF_DEVICE(PDYY,ZTMP2_DEVICE)
+!$acc kernels
+    PLM(:,:,:) = BR_POW( PLM(:,:,:)*ZTMP1_DEVICE    *ZTMP2_DEVICE     , 1./3. )
+!$acc end kernels
+#endif !_OPENACC
+#endif !MNH_BITREP
   END IF
 END IF
 !   compute a mixing length limited by the stability