Fix omp reproducible

src/gen_modules_diag.F90

@@ -3254,7 +3254,11 @@ subroutine dvd_add_difflux_bhvisc(do_SDdvd, tr_num, dvd_tot, tr, trstar, gamma0_
     ! first round 
 !$OMP PARALLEL DEFAULT(SHARED) PRIVATE(edge, nz, ednodes, edelem, elnodes_l, elnodes_r, &
 !$OMP nu1, nl1, du, dv, dt, len, vi)
-!$OMP DO    
+#if defined(__openmp_reproducible)
+!$OMP DO ORDERED
+#else
+!$OMP DO
+#endif
     do edge=1, myDim_edge2D!+eDim_edge2D
         ! skip boundary edges only consider inner edges 
         if (myList_edge2D(edge) > edge2D_in) cycle
@@ -3301,7 +3305,11 @@ subroutine dvd_add_difflux_bhvisc(do_SDdvd, tr_num, dvd_tot, tr, trstar, gamma0_
 
     !___________________________________________________________________________
     ! second round:
+#if defined(__openmp_reproducible)
+!$OMP DO ORDERED
+#else
 !$OMP DO
+#endif
     do edge=1, myDim_edge2D!+eDim_edge2D
         ! skip boundary edges only consider inner edges 
         if (myList_edge2D(edge) > edge2D_in) cycle

src/gen_support.F90

@@ -337,13 +337,16 @@ subroutine integrate_nod_2D(data, int2D, partit, mesh)
 #if !defined(__openmp_reproducible)
 !$OMP PARALLEL DEFAULT(SHARED) PRIVATE(row)
 !$OMP DO REDUCTION (+: lval)
-#endif
   do row=1, myDim_nod2D
      lval=lval+data(row)*areasvol(ulevels_nod2D(row),row)
   end do
-#if !defined(__openmp_reproducible) 
 !$OMP END DO
 !$OMP END PARALLEL
+#else
+  ! Use serial computation for reproducible results
+  do row=1, myDim_nod2D
+     lval=lval+data(row)*areasvol(ulevels_nod2D(row),row)
+  end do
 #endif
   int2D=0.0_WP
   call MPI_AllREDUCE(lval, int2D, 1, MPI_DOUBLE_PRECISION, MPI_SUM, &
@@ -374,21 +377,32 @@ subroutine integrate_nod_3D(data, int3D, partit, mesh)
 #include "associate_mesh_ass.h" 
 
   lval=0.0_WP
+#if defined(__openmp_reproducible)
+!$OMP PARALLEL DEFAULT(SHARED) PRIVATE(row, k, lval_row)
+!$OMP DO ORDERED
+#else
 !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(row, k, lval_row) REDUCTION(+: lval)
+#endif
   do row=1, myDim_nod2D
      lval_row = 0.
      do k=ulevels_nod2D(row), nlevels_nod2D(row)-1
         lval_row=lval_row+data(k, row)*areasvol(k,row)*hnode_new(k,row)  ! --> TEST_cavity
      end do
 #if defined(__openmp_reproducible)
 !$OMP ORDERED
-#endif
+!$OMP ATOMIC UPDATE
      lval = lval + lval_row
-#if defined(__openmp_reproducible)
 !$OMP END ORDERED
+#else
+     lval = lval + lval_row
 #endif
   end do
+#if defined(__openmp_reproducible)
+!$OMP END DO
+!$OMP END PARALLEL
+#else
 !$OMP END PARALLEL DO
+#endif
 
   int3D=0.0_WP
   call MPI_AllREDUCE(lval, int3D, 1, MPI_DOUBLE_PRECISION, MPI_SUM, &
@@ -575,7 +589,7 @@ FUNCTION omp_min_max_sum2(arr, pos11, pos12, pos21, pos22, what, partit, nan)
   character(3),  intent(in)   :: what
   real(kind=WP), optional     :: nan !to be implemented upon the need (for masked arrays)
   real(kind=WP)               :: omp_min_max_sum2
-  real(kind=WP)               :: val, vmasked, val_part(pos11:pos12)
+  real(kind=WP)               :: val, vmasked, val_part(pos21:pos22)
   integer                     :: i, j
 
 
@@ -663,7 +677,12 @@ subroutine integrate_elem_3D(data, int3D, partit, mesh)
 #include "associate_mesh_ass.h" 
 
   lval=0.0_WP
+#if defined(__openmp_reproducible)
+!$OMP PARALLEL DEFAULT(SHARED) PRIVATE(row, k, lval_row)
+!$OMP DO ORDERED
+#else
 !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(row, k, lval_row) REDUCTION(+: lval)
+#endif
   do row=1, myDim_elem2D
      if(elem2D_nodes(1, row) > myDim_nod2D) cycle
      lval_row = 0.
@@ -672,13 +691,19 @@ subroutine integrate_elem_3D(data, int3D, partit, mesh)
      end do
 #if defined(__openmp_reproducible)
 !$OMP ORDERED
-#endif
+!$OMP ATOMIC UPDATE
      lval = lval + lval_row
-#if defined(__openmp_reproducible)
 !$OMP END ORDERED
+#else
+     lval = lval + lval_row
 #endif
   end do
+#if defined(__openmp_reproducible)
+!$OMP END DO
+!$OMP END PARALLEL
+#else
 !$OMP END PARALLEL DO
+#endif
 
   int3D=0.0_WP
   call MPI_AllREDUCE(lval, int3D, 1, MPI_DOUBLE_PRECISION, MPI_SUM, &
@@ -707,7 +732,12 @@ subroutine integrate_elem_2D(data, int2D, partit, mesh)
 #include "associate_mesh_ass.h" 
 
   lval=0.0_WP
+#if defined(__openmp_reproducible)
+!$OMP PARALLEL DEFAULT(SHARED) PRIVATE(row)
+!$OMP DO ORDERED
+#else
 !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(row) REDUCTION(+: lval)
+#endif
   do row=1, myDim_elem2D
      if(elem2D_nodes(1, row) > myDim_nod2D) cycle
 #if defined(__openmp_reproducible)
@@ -718,7 +748,12 @@ subroutine integrate_elem_2D(data, int2D, partit, mesh)
 !$OMP END ORDERED
 #endif
   end do
+#if defined(__openmp_reproducible)
+!$OMP END DO
+!$OMP END PARALLEL
+#else
 !$OMP END PARALLEL DO
+#endif
 
   int2D=0.0_WP
   call MPI_AllREDUCE(lval, int2D, 1, MPI_DOUBLE_PRECISION, MPI_SUM, &

src/icb_step.F90

@@ -567,7 +567,7 @@ subroutine iceberg_step1(ice, mesh, partit, dynamics, ib, height_ib_single,lengt
      case(2) 
       area_ib_tot = length_ib_single*width_ib_single*scaling(ib)
     end select
-!$OMP PARALLEL DEFAULT(SHARED) PRIVATE(idx, area_ib_tot)
+!$OMP PARALLEL DEFAULT(SHARED) PRIVATE(idx) REDUCTION(+:area_ib_tot)
 !$OMP DO
     do idx = 1, size(elem_block)
         if (elem_block(idx) == iceberg_elem) then
@@ -790,7 +790,7 @@ subroutine iceberg_step2(mesh, partit,arr, elem_from_block, ib, height_ib_single
       case(2) 
        area_ib_tot = length_ib_single*width_ib_single*scaling(ib)
      end select
-!$OMP PARALLEL DEFAULT(SHARED) PRIVATE(idx, area_ib_tot)
+!$OMP PARALLEL DEFAULT(SHARED) PRIVATE(idx) REDUCTION(+:area_ib_tot)
 !$OMP DO
      do idx = 1, size(elem_block_red)
          if (elem_block_red(idx) == iceberg_elem) then

src/ice_EVP.F90

@@ -227,7 +227,11 @@ subroutine stress2rhs(ice, partit, mesh)
 #endif
 
 #ifndef ENABLE_OPENACC
+#if defined(__openmp_reproducible)
+!$OMP DO ORDERED
+#else
 !$OMP DO
+#endif
 #else
     !$ACC PARALLEL LOOP GANG VECTOR DEFAULT(PRESENT)
 #if !defined(DISABLE_OPENACC_ATOMICS)
@@ -716,7 +720,11 @@ subroutine EVPdynamics(ice, partit, mesh)
         ! apply sea ice velocity boundary condition
 
 #ifndef ENABLE_OPENACC
+#if defined(__openmp_reproducible)
+!$OMP DO ORDERED
+#else
 !$OMP DO
+#endif
 #else
         ! With the binary data of np2 goes only inside the first if
         !$ACC PARALLEL LOOP GANG VECTOR DEFAULT(PRESENT)

src/ice_fct.F90

@@ -759,7 +759,11 @@ subroutine ice_fem_fct(tr_array_id, ice, partit, mesh)
 #endif
 
 #ifndef ENABLE_OPENACC
+#if defined(__openmp_reproducible)
+!$OMP DO ORDERED
+#else
 !$OMP DO
+#endif
 #else
 #if !defined(DISABLE_OPENACC_ATOMICS)
     !$ACC PARALLEL LOOP GANG VECTOR PRIVATE(elnodes) DEFAULT(PRESENT)
@@ -901,7 +905,11 @@ subroutine ice_fem_fct(tr_array_id, ice, partit, mesh)
 #endif
 
 #ifndef ENABLE_OPENACC
+#if defined(__openmp_reproducible)
+!$OMP DO ORDERED
+#else
 !$OMP DO
+#endif
 #else
 #if !defined(DISABLE_OPENACC_ATOMICS)
         !$ACC PARALLEL LOOP GANG VECTOR PRIVATE(elnodes) DEFAULT(PRESENT)
@@ -959,7 +967,11 @@ subroutine ice_fem_fct(tr_array_id, ice, partit, mesh)
         end do
 #ifndef ENABLE_OPENACC
 !$OMP END DO
+#if defined(__openmp_reproducible)
+!$OMP DO ORDERED
+#else
 !$OMP DO
+#endif
 #else
         !$ACC END PARALLEL LOOP
 #if !defined(DISABLE_OPENACC_ATOMICS)
@@ -1018,7 +1030,11 @@ subroutine ice_fem_fct(tr_array_id, ice, partit, mesh)
         end do
 #ifndef ENABLE_OPENACC
 !$OMP END DO
+#if defined(__openmp_reproducible)
+!$OMP DO ORDERED
+#else
 !$OMP DO
+#endif
 #else
         !$ACC END PARALLEL LOOP
 #if !defined(DISABLE_OPENACC_ATOMICS)
@@ -1078,7 +1094,11 @@ subroutine ice_fem_fct(tr_array_id, ice, partit, mesh)
         end do
 #ifndef ENABLE_OPENACC
 !$OMP END DO
+#if defined(__openmp_reproducible)
+!$OMP DO ORDERED
+#else
 !$OMP DO
+#endif
 #else
         !$ACC END PARALLEL LOOP
 #endif
@@ -1166,8 +1186,12 @@ SUBROUTINE ice_mass_matrix_fill(ice, partit, mesh)
     mass_matrix => ice%work%fct_massmatrix(:)
     !
     ! a)
-!$OMP PARALLEL DEFAULT(SHARED) PRIVATE(n, k, row, elem, elnodes, q, offset, ipos, aa)
+!$OMP PARALLEL DEFAULT(SHARED) PRIVATE(n, k, row, elem, elnodes, q, offset, ipos, aa, flag, iflag)
+#if defined(__openmp_reproducible)
+!$OMP DO ORDERED
+#else
 !$OMP DO
+#endif
     DO elem=1,myDim_elem2D
         elnodes=elem2D_nodes(:,elem)
 
@@ -1326,7 +1350,11 @@ subroutine ice_TG_rhs_div(ice, partit, mesh)
 
 #ifndef ENABLE_OPENACC
 !$OMP PARALLEL DEFAULT(SHARED) PRIVATE(diff, entries, um, vm, vol, dx, dy, n, q, row, elem, elnodes, c1, c2, c3, c4, cx1, cx2, cx3, cx4, entries2)
+#if defined(__openmp_reproducible)
+!$OMP DO ORDERED
+#else
 !$OMP DO
+#endif
 #else
 #if !defined(DISABLE_OPENACC_ATOMICS)
     !$ACC PARALLEL LOOP GANG VECTOR PRIVATE(elnodes, dx, dy, entries, entries2) DEFAULT(PRESENT)

src/ice_maEVP.F90

@@ -245,7 +245,11 @@ subroutine ssh2rhs(ice, partit, mesh)
     !_____________________________________________________________________________
     ! use floating sea ice for zlevel and zstar
     if (use_floatice .and.  .not. trim(which_ale)=='linfs') then
+#if defined(__openmp_reproducible)
+!$OMP DO ORDERED
+#else
 !$OMP DO
+#endif
         do elem=1,myDim_elem2d
             elnodes=elem2D_nodes(:,elem)
             !_______________________________________________________________________
@@ -285,7 +289,11 @@ subroutine ssh2rhs(ice, partit, mesh)
         end do
 !$OMP END DO
     else
+#if defined(__openmp_reproducible)
+!$OMP DO ORDERED
+#else
 !$OMP DO
+#endif
         do elem=1,myDim_elem2d
             elnodes=elem2D_nodes(:,elem)
             !_______________________________________________________________________
@@ -371,7 +379,11 @@ subroutine stress2rhs_m(ice, partit, mesh)
 !$OMP END PARALLEL DO
 
 !$OMP PARALLEL DEFAULT(SHARED) PRIVATE(elem, elnodes, k, row, dx, dy, vol, mf, aa, bb, mass, cluster_area, elevation_elem)
+#if defined(__openmp_reproducible)
+!$OMP DO ORDERED
+#else
 !$OMP DO
+#endif
     do elem=1,myDim_elem2d
         elnodes=elem2D_nodes(:,elem)
         !_______________________________________________________________________
@@ -545,7 +557,11 @@ subroutine EVPdynamics_m(ice, partit, mesh)
     ! use floating sea ice for zlevel and zstar
 !$OMP PARALLEL DEFAULT(SHARED) PRIVATE(el, elnodes, vol, dx, dy, p_ice, n, bb, aa)
     if (use_floatice .and.  .not. trim(which_ale)=='linfs') then
+#if defined(__openmp_reproducible)
+!$OMP DO ORDERED
+#else
 !$OMP DO
+#endif
         do el=1,myDim_elem2d
             elnodes=elem2D_nodes(:,el)
 
@@ -588,7 +604,11 @@ subroutine EVPdynamics_m(ice, partit, mesh)
     !_____________________________________________________________________________
     ! use levitating sea ice for linfs, zlevel and zstar
     else
+#if defined(__openmp_reproducible)
+!$OMP DO ORDERED
+#else
 !$OMP DO
+#endif
         do el=1,myDim_elem2d
             elnodes=elem2D_nodes(:,el)
             !_______________________________________________________________________
@@ -685,7 +705,11 @@ subroutine EVPdynamics_m(ice, partit, mesh)
         ! SD, 30.07.2014
         !_______________________________________________________________________
 !$OMP PARALLEL DEFAULT(SHARED) PRIVATE(el, i, ed, row, elnodes, dx, dy, meancos, eps1, eps2, delta, pressure, umod, drag, rhsu, rhsv, det, n)
+#if defined(__openmp_reproducible)
+!$OMP DO ORDERED
+#else
 !$OMP DO
+#endif
         do el=1,myDim_elem2D
             if (ulevels(el)>1) cycle
 

src/oce_adv_tra_driver.F90

@@ -131,7 +131,11 @@ subroutine do_oce_adv_tra(dt, vel, w, wi, we, tr_num, dynamics, tracers, partit,
 #endif
 
 #ifndef ENABLE_OPENACC
+#if defined(__openmp_reproducible)
+!$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(e, enodes, el, nl1, nu1, nl2, nu2, nu12, nl12, nz) ORDERED
+#else
 !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(e, enodes, el, nl1, nu1, nl2, nu2, nu12, nl12, nz)
+#endif
 #else
 #if !defined(DISABLE_OPENACC_ATOMICS)
         !$ACC PARALLEL LOOP GANG PRIVATE(enodes, el) DEFAULT(PRESENT) VECTOR_LENGTH(acc_vl)
@@ -493,7 +497,11 @@ subroutine oce_tra_adv_flux2dtracer(dt, dttf_h, dttf_v, flux_h, flux_v, partit,
 #endif
     ! Horizontal
 #ifndef ENABLE_OPENACC
+#if defined(__openmp_reproducible)
+!$OMP DO ORDERED
+#else
 !$OMP DO
+#endif
 #else
 #if !defined(DISABLE_OPENACC_ATOMICS)
     !$ACC PARALLEL LOOP GANG PRIVATE(enodes, el) DEFAULT(PRESENT) VECTOR_LENGTH(acc_vl)

src/oce_adv_tra_fct.F90

@@ -267,7 +267,7 @@ subroutine oce_tra_adv_fct(dt, ttf, lo, adf_h, adf_v, fct_ttf_min, fct_ttf_max,
 #endif
     !Vertical
 #ifndef ENABLE_OPENACC
-!$OMP DO
+!$OMP DO 
 #else
     !$ACC PARALLEL LOOP GANG DEFAULT(PRESENT) VECTOR_LENGTH(acc_vl)
 #endif
@@ -288,7 +288,11 @@ subroutine oce_tra_adv_fct(dt, ttf, lo, adf_h, adf_v, fct_ttf_min, fct_ttf_max,
 #endif
 
 #ifndef ENABLE_OPENACC
+#if defined(__openmp_reproducible)
+!$OMP DO ORDERED
+#else
 !$OMP DO
+#endif
 #else
     !Horizontal
 #if !defined(DISABLE_OPENACC_ATOMICS)