Update ocean BGC threading

To make multi-threaded runs BFB with single-threaded runs in coupled cases.

Author: amametjanov

src/core_ocean/mode_forward/mpas_ocn_time_integration_rk4.F

@@ -439,7 +439,6 @@ subroutine ocn_time_integrator_rk4(domain, dt)!{{{
                     * highFreqThicknessTend(:, iCell)
               end do
               !$omp end do
-              call mpas_threading_barrier()
               block => block % next
            end do
 
@@ -642,7 +641,6 @@ subroutine ocn_time_integrator_rk4(domain, dt)!{{{
             normalTransportVelocity(:, iEdge) = normalVelocityNew(:, iEdge)
          end do
          !$omp end do
-         call mpas_threading_barrier()
 
          ! Compute normalGMBolusVelocity and the tracer transport velocity
          if (config_use_standardGM) then
@@ -1161,7 +1159,6 @@ subroutine ocn_time_integrator_rk4_diagnostic_update(block, dt, rkWeight, err)!{
          normalTransportVelocity(:, iEdge) = normalVelocityProvis(:, iEdge)
       end do
       !$omp end do
-      call mpas_threading_barrier()
 
       ! Compute normalGMBolusVelocity, relativeSlope and RediDiffVertCoef if respective flags are turned on
       if (config_use_standardGM) then
@@ -1176,7 +1173,6 @@ subroutine ocn_time_integrator_rk4_diagnostic_update(block, dt, rkWeight, err)!{
          end do
          !$omp end do
       end if
-      call mpas_threading_barrier()
       ! ------------------------------------------------------------------
       ! End: Accumulating various parametrizations of the transport velocity
       ! ------------------------------------------------------------------

src/core_ocean/mode_forward/mpas_ocn_time_integration_split.F

@@ -560,7 +560,7 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{
 
                allocate(uTemp(nVertLevels))
 
-               !$omp do schedule(runtime) private(cell1, cell2, k, normalThicknessFluxSum, thicknessSum)
+               !$omp do schedule(runtime)
                do iEdge = 1, nEdges
                   cell1 = cellsOnEdge(1,iEdge)
                   cell2 = cellsOnEdge(2,iEdge)
@@ -813,7 +813,7 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{
                      nEdges = nEdgesPtr
                      nEdges = nEdgesArray( edgeHaloComputeCounter )
 
-                     !$omp do schedule(runtime) private(cell1, cell2, CoriolisTerm, i, eoe)
+                     !$omp do schedule(runtime)
                      do iEdge = 1, nEdges
 
                         temp_mask = edgeMask(1, iEdge)
@@ -1025,7 +1025,7 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{
 
                    nEdges = nEdgesArray( edgeHaloComputeCounter )
 
-                   !$omp do schedule(runtime) private(cell1, cell2, eoe, CoriolisTerm, i, sshCell1, sshCell2)
+                   !$omp do schedule(runtime)
                    do iEdge = 1, nEdges
 
                      ! asarje: added to avoid redundant computations based on mask
@@ -1338,7 +1338,7 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{
                   useVelocityCorrection = 0
                endif
 
-               !$omp do schedule(runtime) private(k, normalThicknessFluxSum, thicknessSum, normalVelocityCorrection)
+               !$omp do schedule(runtime)
                do iEdge = 1, nEdges
 
                   ! velocity for normalVelocityCorrectionection is normalBarotropicVelocity + normalBaroclinicVelocity + uBolus

src/core_ocean/shared/mpas_ocn_diagnostics.F

@@ -285,12 +285,13 @@ subroutine ocn_diagnostic_solve(dt, statePool, forcingPool, meshPool, diagnostic
       !    used -1e34 so error clearly occurs if these values are used.
       !
 
-      !$omp single
+      !$omp master
       normalVelocity(:,nEdges+1) = -1e34
       layerThickness(:,nCells+1) = -1e34
       activeTracers(indexTemperature,:,nCells+1) = -1e34
       activeTracers(indexSalinity,:,nCells+1) = -1e34
-      !$omp end single
+      !$omp end master
+      call mpas_threading_barrier()
 
       call ocn_relativeVorticity_circulation(relativeVorticity, circulation, meshPool, normalVelocity, err)
 
@@ -320,7 +321,7 @@ subroutine ocn_diagnostic_solve(dt, statePool, forcingPool, meshPool, diagnostic
       !
       allocate(div_hu(nVertLevels),div_huTransport(nVertLevels),div_huGMBolus(nVertLevels))
 
-      !$omp do schedule(runtime) private(invAreaCell1, iEdge, r_tmp, i, k)
+      !$omp do schedule(runtime)
       do iCell = 1, nCells
          divergence(:, iCell) = 0.0_RKIND
          kineticEnergyCell(:, iCell) = 0.0_RKIND
@@ -364,7 +365,7 @@ subroutine ocn_diagnostic_solve(dt, statePool, forcingPool, meshPool, diagnostic
 
       nEdges = nEdgesArray( 2 )
 
-      !$omp do schedule(runtime) private(eoe, i, k)
+      !$omp do schedule(runtime)
       do iEdge = 1, nEdges
          tangentialVelocity(:, iEdge) = 0.0_RKIND
          ! Compute v (tangential) velocities
@@ -621,7 +622,7 @@ subroutine ocn_diagnostic_solve(dt, statePool, forcingPool, meshPool, diagnostic
 
         allocate(pTop(nVertLevels))
 
-        !$omp do schedule(runtime) private(k)
+        !$omp do schedule(runtime)
         do iCell = 1, nCells
 
            ! assume atmospheric pressure at the surface is zero for now.
@@ -752,7 +753,7 @@ subroutine ocn_diagnostic_solve(dt, statePool, forcingPool, meshPool, diagnostic
 
         nCells = nCellsArray( 1 )
 
-        !$omp do schedule(runtime) private(surfaceLayerDepth, sumSurfaceLayer, k, rSurfaceLayer)
+        !$omp do schedule(runtime)
         do iCell=1,nCells
           surfaceLayerDepth = config_cvmix_kpp_surface_layer_averaging
           sumSurfaceLayer=0.0_RKIND
@@ -1168,7 +1169,7 @@ subroutine ocn_filter_btr_mode_vel(statePool, diagnosticsPool, meshPool, timeLev
 
       call mpas_pool_get_dimension(meshPool, 'nEdges', nEdges)
 
-      !$omp do schedule(runtime) private(normalThicknessFluxSum, thicknessSum, k, vertSum)
+      !$omp do schedule(runtime)
       do iEdge = 1, nEdges
 
         ! thicknessSum is initialized outside the loop because on land boundaries
@@ -1189,8 +1190,6 @@ subroutine ocn_filter_btr_mode_vel(statePool, diagnosticsPool, meshPool, timeLev
       enddo ! iEdge
       !$omp end do
 
-      call mpas_threading_barrier()
-
       call mpas_timer_stop("ocn_filter_btr_mode_vel")
 
    end subroutine ocn_filter_btr_mode_vel!}}}
@@ -1239,7 +1238,7 @@ subroutine ocn_filter_btr_mode_tend_vel(tendPool, statePool, diagnosticsPool, me
 
       call mpas_pool_get_dimension(meshPool, 'nEdges', nEdges)
 
-      !$omp do schedule(runtime) private(normalThicknessFluxSum, thicknessSum, vertSum, k)
+      !$omp do schedule(runtime)
       do iEdge = 1, nEdges
 
         ! thicknessSum is initialized outside the loop because on land boundaries
@@ -1450,8 +1449,7 @@ subroutine ocn_compute_KPP_input_fields(statePool, forcingPool, meshPool, diagno
                                          densitySurfaceDisplaced, err, thermalExpansionCoeff, salineContractionCoeff, &
                                          timeLevel)
 
-      !$omp do schedule(runtime) private(invAreaCell, deltaVelocitySquared, i, iEdge, factor, delU2, fracAbsorbed, &
-      !$omp                              fracAbsorbedRunoff)
+      !$omp do schedule(runtime)
       do iCell = 1, nCells
        invAreaCell = 1.0_RKIND / areaCell(iCell)
 
@@ -1514,8 +1512,6 @@ subroutine ocn_compute_KPP_input_fields(statePool, forcingPool, meshPool, diagno
       enddo
       !$omp end do
 
-      call mpas_threading_barrier()
-
       ! deallocate scratch space
       call mpas_deallocate_scratch_field(thermalExpansionCoeffField, .true.)
       call mpas_deallocate_scratch_field(salineContractionCoeffField, .true.)
@@ -1740,7 +1736,7 @@ subroutine ocn_compute_land_ice_flux_input_fields(meshPool, statePool, &
 
 
       ! average temperature and salinity over horizontal neighbors and the sub-ice-shelf boundary layer
-      !$omp do schedule(runtime) private(blThickness, iLevel, dz)
+      !$omp do schedule(runtime)
       do iCell = 1, nCells
          blThickness = 0.0_RKIND
          blTempScratch(iCell) = 0.0_RKIND

src/core_ocean/shared/mpas_ocn_equation_of_state_jm.F

@@ -309,7 +309,7 @@ subroutine ocn_equation_of_state_jm_density(meshPool, scratchPool, nCells, k_dis
          enddo
       endif
 
-      !$omp do schedule(runtime) private(k, DRDT0, DKDT, DRHODT, DRDS0, DKDS, DRHODS)
+      !$omp do schedule(runtime)
       do iCell=1,nCells
          if (displacement_type == 'surfaceDisplaced') then
            if(present(tracersSurfaceLayerValue)) then

src/core_ocean/shared/mpas_ocn_frazil_forcing.F

@@ -451,9 +451,7 @@ subroutine ocn_frazil_forcing_build_arrays(domain, meshPool, forcingPool, diagno
       nCells = nCellsArray( 2 )
 
       ! loop over all columns
-      !$omp do schedule(runtime) private(kBottomFrazil, k, columnTemperatureMin, sumNewFrazilIceThickness, &
-      !$omp                              oceanFreezingTemperature, potential, freezingEnergy, meltingEnergy, &
-      !$omp                              newFrazilIceThickness, meltedFrazilIceThickness)
+      !$omp do schedule(runtime)
       do iCell=1,nCells
 
         underLandIce = .false.

src/core_ocean/shared/mpas_ocn_gm.F

@@ -418,11 +418,10 @@ subroutine ocn_gm_compute_Bolus_velocity(diagnosticsPool, meshPool, scratchPool)
       ! use relativeSlopeTaperingCell as a temporary space for smoothing of relativeSlopeTopOfCell
       relativeSlopeTaperingCell = relativeSlopeTopOfCell
       do iter = 1, 5
-         !$omp barrier
 
          nCells = nCellsArray( 2 )
 
-         !$omp do schedule(runtime) private(k, rtmp)
+         !$omp do schedule(runtime)
          do iCell=1,nCells
            relativeSlopeTaperingCell(1, iCell) = 0.0_RKIND
            relativeSlopeTaperingCell(maxLevelCell(iCell):nVertLevels, iCell) = 0.0_RKIND
@@ -523,7 +522,7 @@ subroutine ocn_gm_compute_Bolus_velocity(diagnosticsPool, meshPool, scratchPool)
 
       nEdges = nEdgesArray( 3 )
 
-      !$omp do schedule(runtime) private(cell1, cell2, k, BruntVaisalaFreqTopEdge, N)
+      !$omp do schedule(runtime)
       do iEdge = 1, nEdges
          cell1 = cellsOnEdge(1,iEdge)
          cell2 = cellsOnEdge(2,iEdge)

src/core_ocean/shared/mpas_ocn_sea_ice.F

@@ -135,8 +135,7 @@ subroutine ocn_sea_ice_formation(meshPool, indexTemperature, indexSalinity, laye
 
       density_ice = rho_ice
 
-      !$omp do schedule(runtime) private(maxLevel, netEnergyChange, k, freezingTemp, availableEnergyChange, energyChange, &
-      !$omp                              temperatureChange, thicknessChange, iceThicknessChange, iTracer)
+      !$omp do schedule(runtime)
       do iCell = 1, nCellsSolve ! Check performance of these two loop definitions
 !     do iCell = nCellsSolve, 1, -1
          maxLevel = min(maxLevelCell(iCell), verticalLevelCap)

src/core_ocean/shared/mpas_ocn_surface_land_ice_fluxes.F

@@ -491,7 +491,7 @@ subroutine ocn_surface_land_ice_fluxes_build_arrays(meshPool, diagnosticsPool, &
       nCells = nCellsArray( size(nCellsArray) )
 
       if(isomipOn) then
-         !$omp do schedule(runtime) private(heatFlux)
+         !$omp do schedule(runtime)
          do iCell = 1, nCells
             if (landIceMask(iCell) == 0) cycle
 
@@ -795,7 +795,7 @@ subroutine compute_melt_fluxes( &
             .and. present(kappa_land_ice)
     Tlatent = latent_heat_fusion_mks/cp_sw
 
-    !$omp do schedule(runtime) private(iceHeatFluxCoeff, nu, iceDeltaT, T0, transferVelocityRatio, a, b, c)
+    !$omp do schedule(runtime)
     do iCell = 1, nCells
       if (mask(iCell) == 0) cycle
 
@@ -952,7 +952,7 @@ subroutine compute_HJ99_melt_fluxes( &
 
     err = 0
     cpRatio = cp_land_ice/cp_sw
-    !$omp do schedule(runtime) private(T0, transferVelocityRatio, Tlatent, eta, TlatentStar, a, b, c)
+    !$omp do schedule(runtime) reduction(+:err)
     do iCell = 1, nCells
       if (mask(iCell) == 0) cycle
 
@@ -975,9 +975,7 @@ subroutine compute_HJ99_melt_fluxes( &
       outInterfaceSalinity(iCell) = (-b + sqrt(b**2 - 4.0_RKIND*a*c*oceanSalinity(iCell)))/(2.0_RKIND*a)
       if (outInterfaceSalinity(iCell) .le. 0.0_RKIND) then
          err = 1
-         call mpas_log_write( &
-            'interfaceSalinity is negative.', &
-            MPAS_LOG_CRIT)
+         call mpas_log_write('interfaceSalinity is negative.', MPAS_LOG_CRIT)
       end if
       outInterfaceTemperature(iCell) = dTf_dS*outInterfaceSalinity(iCell)+T0
 

src/core_ocean/shared/mpas_ocn_tendency.F

@@ -628,6 +628,9 @@ subroutine ocn_tend_tracer(tendPool, statePool, forcingPool, diagnosticsPool, me
                                                         tracerGroupSurfaceFluxRemoved, dt, layerThickness, err)
                end if
 
+
+               !$omp master
+
                !
                ! compute ecosystem source-sink tendencies and net surface fluxes
                ! NOTE: must be called before ocn_tracer_surface_flux_tend
@@ -673,6 +676,9 @@ subroutine ocn_tend_tracer(tendPool, statePool, forcingPool, diagnosticsPool, me
                      nTracerGroup, nCellsSolve, zMid, indexTemperature, indexSalinity, tracerGroupSurfaceFlux, err)
                endif
 
+               !$omp end master
+               call mpas_threading_barrier()
+
                !
                ! ocean surface restoring
                !
@@ -974,7 +980,7 @@ subroutine ocn_tend_freq_filtered_thickness(tendPool, statePool, diagnosticsPool
 
       allocate(div_hu(nVertLevels))
 
-      !$omp do schedule(runtime) private(div_hu_btr, invAreaCell, i, iEdge, k, totalThickness, flux)
+      !$omp do schedule(runtime)
       do iCell = 1, nCells
         tend_lowFreqDivergence(:, iCell) = 0.0_RKIND
         tend_highFreqThickness(:, iCell) = 0.0_RKIND

src/core_ocean/shared/mpas_ocn_thick_ale.F

@@ -147,7 +147,7 @@ subroutine ocn_ALE_thickness(meshPool, verticalMeshPool, SSH, ALE_thickness, err
       !
       ! ALE thickness alteration due to SSH (z-star)
       !
-      !$omp do schedule(runtime) private(kMax, thicknessSum, k)
+      !$omp do schedule(runtime)
       do iCell = 1, nCells
          kMax = maxLevelCell(iCell)
 
@@ -182,7 +182,7 @@ subroutine ocn_ALE_thickness(meshPool, verticalMeshPool, SSH, ALE_thickness, err
       !
       if (config_use_min_max_thickness) then
 
-         !$omp do schedule(runtime) private(kMax, remainder, k, newThickness)
+         !$omp do schedule(runtime)
          do iCell = 1, nCells
             kMax = maxLevelCell(iCell)