From cb6caf3c15cf65d541d54388978795107402ac89 Mon Sep 17 00:00:00 2001
From: Juan ESCOBAR <juan.escobar@aero.obs-mip.fr>
Date: Wed, 30 Mar 2022 14:34:42 +0200
Subject: [PATCH] Juan 30/03/2022:Rules.LXcray+ZSOLVER/: For Bit
 Reproductibility on Cray+OMP ,add key MNH_BITREP_OMP in
 Rules.LXcray.mk,BITREP/br_transcendentals.cpp,ZSOLVER/pressurez.f90

---
 src/LIB/BITREP/br_transcendentals.cpp | 41 +++++++++++++++++++++++++--
 src/Rules.LXcray.mk                   | 21 ++++++++++++--
 src/ZSOLVER/pressurez.f90             | 14 +++++----
 3 files changed, 66 insertions(+), 10 deletions(-)

diff --git a/src/LIB/BITREP/br_transcendentals.cpp b/src/LIB/BITREP/br_transcendentals.cpp
index 2bbbb4dae..f0ce340a7 100644
--- a/src/LIB/BITREP/br_transcendentals.cpp
+++ b/src/LIB/BITREP/br_transcendentals.cpp
@@ -59,7 +59,11 @@ double log(double);
  * HELPER FUNCTIONS *
  ********************/
 
+#ifdef MNH_BITREP_OMP
+#pragma omp declare target
+#else
 #pragma acc routine seq
+#endif
 double __internal_copysign_pos(double a, double b)
 {
     union {
@@ -263,7 +267,11 @@ double __internal_asin_kernel(double x)
   return r;
 }
 
+#ifdef MNH_BITREP_OMP
+#pragma omp declare target
+#else
 #pragma acc routine seq
+#endif
 double __internal_atan_kernel(double x)
 {
   double t, x2;
@@ -391,7 +399,11 @@ double acos(double x)
     return t0;
 }
 
+#ifdef MNH_BITREP_OMP
+#pragma omp declare target
+#else
 #pragma acc routine seq
+#endif
 double atan(double x)
 {
     double t0, t1;
@@ -665,8 +677,11 @@ double atanh(double x)
  **************/
 
 
-
+#ifdef MNH_BITREP_OMP
+#pragma omp declare target
+#else
 #pragma acc routine seq
+#endif
 double log(double x)
 {
     double m, f, g, u, v, tmp, q, ulo, log_lo, log_hi;
@@ -817,7 +832,11 @@ double __internal_exp_scale(double x, int i)
     return x;
 }
 
+#ifdef MNH_BITREP_OMP
+#pragma omp declare target
+#else
 #pragma acc routine seq
+#endif
 double __internal_exp_kernel(double x, int scale)
 { 
   double t, z;
@@ -838,7 +857,11 @@ double __internal_exp_kernel(double x, int scale)
   return z;
 }   
 
+#ifdef MNH_BITREP_OMP
+#pragma omp declare target
+#else
 #pragma acc routine seq
+#endif
 double exp(double x)
 {
   double t;
@@ -862,7 +885,9 @@ double exp(double x)
   }
   return t;
 }
-
+#ifdef MNH_BITREP_OMP
+#pragma omp end declare target
+#endif
 
 } // End of namespace bitrep
 
@@ -874,7 +899,11 @@ double br_cos  (double x) { return bitrep::cos  (x); }
 double br_tan  (double x) { return bitrep::tan  (x); }
 double br_asin (double x) { return bitrep::asin (x); }
 double br_acos (double x) { return bitrep::acos (x); }
+#ifdef MNH_BITREP_OMP
+#pragma omp declare target
+#else
 #pragma acc routine seq
+#endif
 double br_atan (double x) { return bitrep::atan (x); }
 double br_sinh (double x) { return bitrep::sinh (x); }
 double br_cosh (double x) { return bitrep::cosh (x); }
@@ -882,9 +911,17 @@ double br_tanh (double x) { return bitrep::tanh (x); }
 double br_asinh(double x) { return bitrep::asinh(x); }
 double br_acosh(double x) { return bitrep::acosh(x); }
 double br_atanh(double x) { return bitrep::atanh(x); }
+#ifdef MNH_BITREP_OMP
+#pragma omp declare target
+#else
 #pragma acc routine seq
+#endif
 double br_log  (double x) { return bitrep::log  (x); }
 double br_log1p(double x) { return bitrep::log1p(x); }
+#ifdef MNH_BITREP_OMP
+#pragma omp declare target
+#else
 #pragma acc routine seq
+#endif
 double br_exp  (double x) { return bitrep::exp  (x); }
 }
diff --git a/src/Rules.LXcray.mk b/src/Rules.LXcray.mk
index 10379a889..1a1633ee6 100644
--- a/src/Rules.LXcray.mk
+++ b/src/Rules.LXcray.mk
@@ -55,11 +55,12 @@ CFLAGS   += -g
 endif
 #
 ifeq "$(OPTLEVEL)" "OPENACC"
+MNH_BITREP_OMP=YES
 CPPFLAGS    += -DMNH_OPENACC -DMNH_GPUDIRECT
 OPT       = $(OPT_BASE) $(OPT_PERF2) $(OPT_OPENACC)
 OPT0      = $(OPT_BASE) $(OPT_PERF0) $(OPT_OPENACC)
 OPT_NOCB  = $(OPT_BASE) $(OPT_PERF2) $(OPT_OPENACC)
-#CXXFLAGS = -acc -Kieee -Mnofma $(OPT_OPENACC)
+#CXXFLAGS = -fopenmp
 #OBJS_REPROD= spll_mode_sum_ll.o
 #$(OBJS_REPROD) : OPT = $(OPT_BASE) $(OPT_PERF2) $(OPT_OPENACC) -Mvect=nosimd -Minfo=all -g
 #OBJS_O1_OPENACC= spll_ice4_tendencies.o spll_turb_ver_thermo_flux.o
@@ -67,6 +68,7 @@ OPT_NOCB  = $(OPT_BASE) $(OPT_PERF2) $(OPT_OPENACC)
 endif
 #
 ifeq "$(OPTLEVEL)" "OPENACCDEFONLY"
+MNH_BITREP_OMP=YES
 CPPFLAGS    += -DMNH_OPENACC -D_FAKEOPENACC
 OPT       = $(OPT_BASE) $(OPT_PERF2) $(OPT_NOOPENACC)  
 OPT0      = $(OPT_BASE) $(OPT_PERF0) $(OPT_NOOPENACC) 
@@ -105,6 +107,7 @@ endif
 FC = ftn
 FCFLAGS = -emf
 CC=cc
+CXX=CC
 export FC CC FCFLAGS
 ifeq "$(VER_MPI)" "MPIAUTO"
 F90 = mpif90
@@ -149,6 +152,19 @@ ifeq "$(MNH_BITREP)" "YES"
 CPPFLAGS_MNH += -DMNH_BITREP
 endif
 #
+# Test of bitrep with OMP compilation 
+#
+ifeq "$(MNH_BITREP_OMP)" "YES"
+CXXFLAGS = -fopenmp
+CPPFLAGS_MNH += -DMNH_BITREP_OMP
+DIR_MASTER += LIB/BITREP
+VPATH += LIB/BITREP
+OBJS_LISTE_MASTER += br_transcendentals.o
+LIBS += -lstdc++
+%.o : %.cpp
+	$(CXX) $(INC) $(CXXFLAGS) $(CPPFLAGS) -c $< -o $(OBJDIR)/$(*F).o
+endif
+#
 # LIBTOOLS flags
 #
 #if MNH_TOOLS exists => compile the tools
@@ -199,7 +215,8 @@ include Makefile.MESONH.mk
 ##########################################################
 # Juan & Maud 20/03/2008 --> Ifort 10.1.008 Bug O2 optimization
 #OPT_PERF1  =  -O1
-OBJS_O1= spll_schu.o spll_ps2str.o spll_p_abs.o spll_ini_one_way_n.o spll_urban_solar_abs.o spll_mode_ekf.o
+OBJS_O1= spll_schu.o spll_ps2str.o spll_ini_one_way_n.o spll_urban_solar_abs.o spll_mode_ekf.o
+#spll_p_abs.o
 $(OBJS_O1) : OPT = $(OPT_BASE) $(OPT_PERF1)
 OBJS_O0= spll_mode_gridproj.o spll_ini_dynamics.o spll_sunpos_n.o spll_average_diag.o
 # spll_ground_param_n.o
diff --git a/src/ZSOLVER/pressurez.f90 b/src/ZSOLVER/pressurez.f90
index 83a84ad32..4cca04485 100644
--- a/src/ZSOLVER/pressurez.f90
+++ b/src/ZSOLVER/pressurez.f90
@@ -263,7 +263,7 @@ USE MODE_MPPDB
 USE MODE_MSG
 USE MODE_SUM2_ll,     ONLY: GMAXLOC_ll
 !
-#ifdef MNH_BITREP
+#if defined(MNH_BITREP) || defined(MNH_BITREP_OMP)
 USE MODI_BITREP
 #endif
 USE MODI_CONJGRAD
@@ -648,7 +648,7 @@ IF(CEQNSYS=='MAE' .OR. CEQNSYS=='DUR') THEN
   ENDIF
   !
 #ifndef MNH_OPENACC   
-#ifndef MNH_BITREP
+#if !defined(MNH_BITREP) && !defined(MNH_BITREP_OMP)
   ZPHIT(:,:,:)=(PPABST(:,:,:)/XP00)**(XRD/XCPD)-PEXNREF(:,:,:)
 #else
   ZPHIT(:,:,:)=BR_POW(PPABST(:,:,:)/XP00,XRD/XCPD)-PEXNREF(:,:,:)
@@ -656,7 +656,7 @@ IF(CEQNSYS=='MAE' .OR. CEQNSYS=='DUR') THEN
 #else
   !$acc kernels
   DO CONCURRENT ( JI=1:IIU,JJ=1:IJU,JK=1:IKU )
-#ifndef MNH_BITREP
+#if !defined(MNH_BITREP) && !defined(MNH_BITREP_OMP)
     ZPHIT(JI,JJ,JK)=(PPABST(JI,JJ,JK)/XP00)**(XRD/XCPD)-PEXNREF(JI,JJ,JK)
 #else
     ZPHIT(JI,JJ,JK)=BR_POW((PPABST(JI,JJ,JK)/XP00),(XRD/XCPD))-PEXNREF(JI,JJ,JK)
@@ -1039,10 +1039,12 @@ IF ((ZMAX_ll > 1.E-12) .AND. KTCOUNT >0 ) THEN
 !
   IF(CEQNSYS=='MAE' .OR. CEQNSYS=='DUR') THEN
      !$acc kernels
-#ifndef MNH_BITREP
-     PPABST(:,:,:)=XP00*(ZPHIT+PEXNREF)**(XCPD/XRD)
+#if !defined(MNH_BITREP) && !defined(MNH_BITREP_OMP)
+     PPABST(:,:,:)=XP00*(ZPHIT(:,:,:)+PEXNREF(:,:,:))**(XCPD/XRD)
 #else
-     PPABST(:,:,:)=XP00*BR_POW((ZPHIT+PEXNREF),(XCPD/XRD))
+     DO CONCURRENT(JI=1:IIU,JJ=1:IJU,JK=1:IKU)
+        PPABST(JI,JJ,JK)=XP00*BR_POW((ZPHIT(JI,JJ,JK)+PEXNREF(JI,JJ,JK)),(XCPD/XRD))
+     END DO
 #endif
      !$acc end kernels
   ELSEIF(CEQNSYS=='LHE') THEN
-- 
GitLab