Merge branch 'brian-oneill/ocean/openacc-tend-thick' into next (PR #4792

)

GPU port of thickness tendency with OpenACC

GPU port of ocn_tend_thick and all subprocesses. Changes include:
* replace mesh pool variables and arrays with counterparts in ocn_mesh
  module
* remove meshPool from ocn_tend_thick subroutine calls in time integrator
  routines
* resolve OpenACC error in mpas_ocn_diagnostics.F by copying needed arrays
  onto device in ocn_compute_land_ice_flux_input_fields routine
* add OpenACC directives to loops in ocn_tend_thick and associated
  subroutines

[BFB]

components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_rk4.F

@@ -890,7 +890,7 @@ subroutine ocn_time_integrator_rk4_compute_thick_tends(block, dt, rkSubstepWeigh
             vertAleTransportTop, err)
       endif
 
-      call ocn_tend_thick(tendPool, forcingPool, meshPool)
+      call ocn_tend_thick(tendPool, forcingPool)
 
 
    end subroutine ocn_time_integrator_rk4_compute_thick_tends!}}}
@@ -1035,7 +1035,7 @@ subroutine ocn_time_integrator_rk4_compute_tends(block, dt, rkWeight, err)!{{{
             vertAleTransportTop, err)
       endif
 
-      call ocn_tend_thick(tendPool, forcingPool, meshPool)
+      call ocn_tend_thick(tendPool, forcingPool)
 
       if (config_filter_btr_mode) then
           call ocn_filter_btr_mode_tend_vel(tendPool, provisStatePool, meshPool, 1)

components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_si.F

@@ -2612,7 +2612,7 @@ subroutine ocn_time_integrator_si(domain, dt)!{{{
 #endif
          call mpas_timer_stop('thick vert trans vel top')
 
-         call ocn_tend_thick(tendPool, forcingPool, meshPool)
+         call ocn_tend_thick(tendPool, forcingPool)
 
          call mpas_timer_stop('si thick tend')
 

components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_split.F

@@ -1516,7 +1516,7 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{
 #endif
          call mpas_timer_stop('thick vert trans vel top')
 
-         call ocn_tend_thick(tendPool, forcingPool, meshPool)
+         call ocn_tend_thick(tendPool, forcingPool)
 
          call mpas_timer_stop('se thick tend')
 

components/mpas-ocean/src/shared/mpas_ocn_diagnostics.F

@@ -2658,6 +2658,8 @@ subroutine ocn_compute_land_ice_flux_input_fields(layerThickness, normalVelocity
 
       ! Compute top drag
 #ifdef MPAS_OPENACC
+      !$acc enter data copyin(landIceFraction)
+
       !$acc parallel loop present(cellsOnEdge, kineticEnergyCell, minLevelCell, minLevelEdgeBot) &
       !$acc    present(landIceFraction, topDrag, normalVelocity)
 #else
@@ -2808,6 +2810,8 @@ subroutine ocn_compute_land_ice_flux_input_fields(layerThickness, normalVelocity
 
       ! modify the spatially-varying attenuation coefficient where there is land ice
 #ifdef MPAS_OPENACC
+      !$acc enter data copyin(landIceMask)
+
       !$acc parallel loop present(landIceMask, sfcFlxAttCoeff)
 #else
       !$omp parallel

components/mpas-ocean/src/shared/mpas_ocn_frazil_forcing.F

@@ -27,6 +27,7 @@ module ocn_frazil_forcing
    use mpas_timer
    use ocn_constants
    use ocn_config
+   use ocn_mesh
    use ocn_diagnostics_variables
    use ocn_equation_of_state
 
@@ -134,14 +135,13 @@ end subroutine ocn_frazil_forcing_tracers!}}}
 !
 !-----------------------------------------------------------------------
 
-   subroutine ocn_frazil_forcing_layer_thickness(meshPool, forcingPool, layerThicknessTend, err)!{{{
+   subroutine ocn_frazil_forcing_layer_thickness(forcingPool, layerThicknessTend, err)!{{{
 
       !-----------------------------------------------------------------
       !
       ! input variables
       !
       !-----------------------------------------------------------------
-      type (mpas_pool_type), intent(in) :: meshPool !< Input: mesh information
       type (mpas_pool_type), intent(in) :: forcingPool !< Input: Forcing information
 
       !-----------------------------------------------------------------
@@ -166,8 +166,6 @@ subroutine ocn_frazil_forcing_layer_thickness(meshPool, forcingPool, layerThickn
       !-----------------------------------------------------------------
 
       integer :: iCell, k, nCells
-      integer, dimension(:), pointer :: nCellsArray
-      integer, dimension(:), pointer :: minLevelCell, maxLevelCell
       real (kind=RKIND), dimension(:,:), pointer :: frazilLayerThicknessTendency
 
       err = 0
@@ -176,24 +174,35 @@ subroutine ocn_frazil_forcing_layer_thickness(meshPool, forcingPool, layerThickn
 
       call mpas_timer_start("frazil thickness tendency")
 
-      call mpas_pool_get_dimension(meshPool, 'nCellsArray', nCellsArray)
-      call mpas_pool_get_array(meshPool, 'minLevelCell', minLevelCell)
-      call mpas_pool_get_array(meshPool, 'maxLevelCell', maxLevelCell)
       call mpas_pool_get_array(forcingPool, 'frazilLayerThicknessTendency', frazilLayerThicknessTendency)
 
       ! frazil fields are needed only over 0 and 1 halos
-      nCells = nCellsArray( 2 )
+      nCells = nCellsHalo( 1 )
 
       ! Build surface fluxes at cell centers
+#ifdef MPAS_OPENACC
+      !$acc enter data copyin(frazilLayerThicknessTendency)
+
+      !$acc parallel loop &
+      !$acc     present(layerThicknessTend, minLevelCell, maxLevelCell, frazilLayerThicknessTendency) &
+      !$acc     private(k)
+#else
       !$omp parallel
       !$omp do schedule(runtime) private(k)
+#endif
       do iCell = 1, nCells
         do k = minLevelCell(iCell), maxLevelCell(iCell)
           layerThicknessTend(k,iCell) = layerThicknessTend(k,iCell) + frazilLayerThicknessTendency(k,iCell)
         end do
       end do
+#ifndef MPAS_OPENACC
       !$omp end do
       !$omp end parallel
+#endif
+
+#ifdef MPAS_OPENACC
+      !$acc exit data delete(frazilLayerThicknessTendency)
+#endif
 
       call mpas_timer_stop("frazil thickness tendency")
 

components/mpas-ocean/src/shared/mpas_ocn_surface_bulk_forcing.F

@@ -279,14 +279,7 @@ end subroutine ocn_surface_bulk_forcing_vel!}}}
 !
 !-----------------------------------------------------------------------
 
-   subroutine ocn_surface_bulk_forcing_thick(meshPool, forcingPool, surfaceThicknessFlux, surfaceThicknessFluxRunoff, err)!{{{
-
-      !-----------------------------------------------------------------
-      !
-      ! input variables
-      !
-      !-----------------------------------------------------------------
-      type (mpas_pool_type), intent(in) :: meshPool !< Input: mesh information
+   subroutine ocn_surface_bulk_forcing_thick(forcingPool, surfaceThicknessFlux, surfaceThicknessFluxRunoff, err)!{{{
 
       !-----------------------------------------------------------------
       !
@@ -313,7 +306,6 @@ subroutine ocn_surface_bulk_forcing_thick(meshPool, forcingPool, surfaceThicknes
       !-----------------------------------------------------------------
 
       integer :: iCell, nCells
-      integer, dimension(:), pointer :: nCellsArray
 
       real (kind=RKIND), dimension(:), pointer :: evaporationFlux, snowFlux
       real (kind=RKIND), dimension(:), pointer :: seaIceFreshWaterFlux, icebergFreshWaterFlux, riverRunoffFlux, iceRunoffFlux
@@ -325,8 +317,6 @@ subroutine ocn_surface_bulk_forcing_thick(meshPool, forcingPool, surfaceThicknes
 
       call mpas_timer_start("bulk_thick", .false.)
 
-      call mpas_pool_get_dimension(meshPool, 'nCellsArray', nCellsArray)
-
       call mpas_pool_get_array(forcingPool, 'evaporationFlux', evaporationFlux)
       call mpas_pool_get_array(forcingPool, 'snowFlux', snowFlux)
       call mpas_pool_get_array(forcingPool, 'seaIceFreshWaterFlux', seaIceFreshWaterFlux)
@@ -335,18 +325,34 @@ subroutine ocn_surface_bulk_forcing_thick(meshPool, forcingPool, surfaceThicknes
       call mpas_pool_get_array(forcingPool, 'iceRunoffFlux', iceRunoffFlux)
       call mpas_pool_get_array(forcingPool, 'rainFlux', rainFlux)
 
-      nCells = nCellsArray( 3 )
+      nCells = nCellsHalo( 2 )
 
       ! Build surface fluxes at cell centers
+#ifdef MPAS_OPENACC
+      !$acc enter data copyin(evaporationFlux, snowFlux, seaIceFreshWaterFlux, icebergFreshWaterFlux, &
+      !$acc     riverRunoffFlux, iceRunoffFlux, rainFlux)
+
+      !$acc parallel loop &
+      !$acc     present(surfaceThicknessFlux, surfaceThicknessFluxRunoff, evaporationFlux, snowFlux, &
+      !$acc     seaIceFreshWaterFlux, icebergFreshWaterFlux, riverRunoffFlux, iceRunoffFlux, rainFlux)
+#else
       !$omp parallel
       !$omp do schedule(runtime)
+#endif
       do iCell = 1, nCells
         surfaceThicknessFlux(iCell) = surfaceThicknessFlux(iCell) + ( snowFlux(iCell) + rainFlux(iCell) + evaporationFlux(iCell) &
                                     + seaIceFreshWaterFlux(iCell) + icebergFreshWaterFlux(iCell) + iceRunoffFlux(iCell) ) / rho_sw
         surfaceThicknessFluxRunoff(iCell) = riverRunoffFlux(iCell) / rho_sw
       end do
+#ifndef MPAS_OPENACC
       !$omp end do
       !$omp end parallel
+#endif
+
+#ifdef MPAS_OPENACC
+      !$acc exit data delete(evaporationFlux, snowFlux, seaIceFreshWaterFlux, icebergFreshWaterFlux, &
+      !$acc     riverRunoffFlux, iceRunoffFlux, rainFlux)
+#endif
 
       call mpas_timer_stop("bulk_thick")
 

components/mpas-ocean/src/shared/mpas_ocn_surface_land_ice_fluxes.F

@@ -245,14 +245,7 @@ end subroutine ocn_surface_land_ice_fluxes_vel!}}}
 !
 !-----------------------------------------------------------------------
 
-   subroutine ocn_surface_land_ice_fluxes_thick(meshPool, forcingPool, surfaceThicknessFlux, err)!{{{
-
-      !-----------------------------------------------------------------
-      !
-      ! input variables
-      !
-      !-----------------------------------------------------------------
-      type (mpas_pool_type), intent(in) :: meshPool !< Input: mesh information
+   subroutine ocn_surface_land_ice_fluxes_thick(forcingPool, surfaceThicknessFlux, err)!{{{
 
       !-----------------------------------------------------------------
       !
@@ -277,7 +270,6 @@ subroutine ocn_surface_land_ice_fluxes_thick(meshPool, forcingPool, surfaceThick
       !-----------------------------------------------------------------
 
       integer :: iCell
-      integer, pointer :: nCells
 
       real (kind=RKIND), dimension(:), pointer :: landIceFreshwaterFlux
 
@@ -287,19 +279,28 @@ subroutine ocn_surface_land_ice_fluxes_thick(meshPool, forcingPool, surfaceThick
 
       call mpas_timer_start("land_ice_thick")
 
-      call mpas_pool_get_dimension(meshPool, 'nCells', nCells)
-
       call mpas_pool_get_array(forcingPool, 'landIceFreshwaterFlux', landIceFreshwaterFlux)
 
       ! Build surface fluxes at cell centers
+#ifdef MPAS_OPENACC
+      !$acc enter data copyin(landIceFreshwaterFlux)
+
+      !$acc parallel loop present(surfaceThicknessFlux, landIceFreshwaterFlux)
+#else
       !$omp parallel
       !$omp do schedule(runtime)
-      do iCell = 1, nCells
+#endif
+      do iCell = 1, nCellsAll
         surfaceThicknessFlux(iCell) = surfaceThicknessFlux(iCell) + landIceFreshwaterFlux(iCell) / rho_sw
       end do
+#ifndef MPAS_OPENACC
       !$omp end do
       !$omp end parallel
+#endif
 
+#ifdef MPAS_OPENACC
+      !$acc exit data delete(landIceFreshwaterFlux)
+#endif
       call mpas_timer_stop("land_ice_thick")
 
    end subroutine ocn_surface_land_ice_fluxes_thick!}}}

components/mpas-ocean/src/shared/mpas_ocn_tendency.F

@@ -111,14 +111,7 @@ module ocn_tendency
 !
 !-----------------------------------------------------------------------
 
-   subroutine ocn_tend_thick(tendPool, forcingPool, meshPool)!{{{
-
-      !-----------------------------------------------------------------
-      ! input variables
-      !-----------------------------------------------------------------
-
-      type (mpas_pool_type), intent(in) :: &
-         meshPool            !< [in] Mesh information (for pass-thru)
+   subroutine ocn_tend_thick(tendPool, forcingPool)!{{{
 
       !-----------------------------------------------------------------
       ! output variables
@@ -166,60 +159,73 @@ subroutine ocn_tend_thick(tendPool, forcingPool, meshPool)!{{{
       ! layer thickness tendency: 
       !    initialize to zero and start accumulating tendency terms
       !
+#ifdef MPAS_OPENACC
+      !$acc enter data create(tendThick, surfaceThicknessFlux, surfaceThicknessFluxRunoff)
 
+      !$acc parallel loop &
+      !$acc     present(tendThick, surfaceThicknessFlux, surfaceThicknessFluxRunoff) &
+      !$acc     private(k)
+#else
       !$omp parallel
       !$omp do schedule(runtime) private(k)
+#endif
       do iCell = 1, nCellsAll
          surfaceThicknessFlux(iCell) = 0.0_RKIND
          surfaceThicknessFluxRunoff(iCell) = 0.0_RKIND
          do k=1,nVertLevels
             tendThick(k, iCell) = 0.0_RKIND
          end do
       end do
+#ifndef MPAS_OPENACC
       !$omp end do
       !$omp end parallel
+#endif
 
       ! If turned off, return with zero fluxes, tendencies
       ! Otherwise, start time and call routines to accumulate
       if (config_disable_thick_all_tend) return
       call mpas_timer_start("ocn_tend_thick")
 
       ! Compute surface mass flux array from bulk forcing
-      call ocn_surface_bulk_forcing_thick(meshPool, forcingPool, &
+      call ocn_surface_bulk_forcing_thick(forcingPool, &
                                           surfaceThicknessFlux, &
                                           surfaceThicknessFluxRunoff, err)
 
       ! Compute surface thickness flux from land ice
-      call ocn_surface_land_ice_fluxes_thick(meshPool, forcingPool, &
+      call ocn_surface_land_ice_fluxes_thick(forcingPool, &
                                              surfaceThicknessFlux, err)
 
       !
       ! Compute horizontal advection term -\nabla\cdot ( hu)
       ! See Ringler et al. (2010) jcp paper, eqn 19, 21, and fig. 3.
       ! for explanation of divergence operator.
       !
-      call ocn_thick_hadv_tend(meshPool, normalTransportVelocity, &
+      call ocn_thick_hadv_tend(normalTransportVelocity, &
                                layerThickEdge, tendThick, err)
 
       ! Compute vertical advection term -d/dz(hw)
-      call ocn_thick_vadv_tend(meshPool, vertAleTransportTop, &
+      call ocn_thick_vadv_tend(vertAleTransportTop, &
                                tendThick, err)
 
       ! Compute surface flux tendency
-      call ocn_thick_surface_flux_tend(meshPool, fractionAbsorbed, &
+      call ocn_thick_surface_flux_tend(fractionAbsorbed, &
                                        fractionAbsorbedRunoff, &
                                        surfaceThicknessFlux, &
                                        surfaceThicknessFluxRunoff, &
                                        tendThick, err)
 
       ! Compute contribution from frazil ice formation
-      call ocn_frazil_forcing_layer_thickness(meshPool, forcingPool, &
+      call ocn_frazil_forcing_layer_thickness(forcingPool, &
                                               tendThick, err)
 
       ! Compute thickness change due to tidal forcing
-      call ocn_tidal_forcing_layer_thickness(meshPool, forcingPool, &
+      call ocn_tidal_forcing_layer_thickness(forcingPool, &
                                              tendThick, err)
 
+#ifdef MPAS_OPENACC
+      !$acc exit data copyout(tendThick, surfaceThicknessFlux, surfaceThicknessFluxRunoff)
+#endif
+
       call mpas_timer_stop("ocn_tend_thick")
 
    !--------------------------------------------------------------------